基于最小噪声分离的钻孔应变数据潮汐应变分量去除方法与流程

本发明涉及钻孔应变数据处理领域，特别涉及一种基于最小噪声分离的钻孔应变数据潮汐应变分量去除方法。

背景技术

地球表面变形和地壳内部的构造运动及其产生的各种地质灾害都与地壳应力作用密切相关，地壳应力应变状态的变化是导致断裂、褶皱乃至发生地震的最直接动因。钻孔应变观测通过对地层内部应变状态依时间连续变化的精细观测，发现并掌握地震应变前兆的(长)中期-短期-临震以及震后调整的时空分布与发展变化规律，是重要的地震前兆观测手段。钻孔应变观测由于其精度高能够清晰的记录到固体潮汐响应。固体潮汐响应会使钻孔应变数据产生潮汐应变分量，掩盖引起构造地震的地壳应变短周期高频信息。最小噪声分离方法理论上可以判定数据内在的维数，分离数据中的信号，可以解决去除钻孔应变数据中潮汐应变分量的问题。

cn106918836a公开了一种基于主成分分析的钻孔应变数据异常提取方法，利用调和分析的方法对钻孔应变数据的周期响应进行去除，并利用主成分分析中的特征值和特征向量分别的将地壳的微弱变化表征出来。

cn104390733a公开了一种地应力大小和方向的确定方法，在测量应力的环氧树脂空心包体上沿圆周方向设置三个应变花，相邻应变花的夹角为120°，每个应变花上有四个应变片，每个应变花上的应变片依次旋转45°设置；利用圆周方向上12个不同方向应变值，能够准确有效的计算出地应力和应力方向。

邱泽华等(用小波-超限率分析提取宁陕台汶川地震体应变异常，2012)采用高通滤波的方法提取钻孔应变数据的高频成分。并利用小波-超限率方法对宁陕台汶川地震提应变异常进行分析。李进武等(钻孔应变仪观测的面应变潮汐因子初步分析，2014)对钻孔应变差分数据进行计算，并分析了潮汐因子变化特征以及影响潮汐因子动态变化的因素。但到目前为止，尚未见运用最小噪声分离的方法对钻孔应变数据潮汐应变分量进行去除的报道。现有方法没有考虑各谐波的物理意义，缺乏针对性；有些需要计算理论固体潮，计算量大。

技术实现要素：

本发明所要解决的技术问题在于提供一种基于最小噪声分离的钻孔应变数据潮汐应变分量去除方法，有针对性的去除钻孔应变数据中潮汐应变分量。

本发明是这样实现的，一种基于最小噪声分离的钻孔应变数据潮汐应变分量去除方法，该方法包括：

a、将钻孔应变数据制作成样本矩阵y(n×1440)；

b、采用自适应窗宽滤波算法对钻孔应变数据样本矩阵y(n×1440)进行高频信号h(n×1440)估计；

c、计算估计的高频信号h(n×1440)的协方差矩阵ch(n×n)，并进行特征值分解；

d、利用高频信号的协方差矩阵ch(n×n)来构造调整矩阵p；

e、利用构造的调整矩阵p对钻孔应变数据样本矩阵的协方差矩阵cy(n×n)进行调整，得到调整后的钻孔应变数据样本矩阵的协方差矩阵cp(n×n)；

f、对调整后的钻孔应变数据样本矩阵的协方差矩阵cp(n×n)进行最小噪声分离变换并得到最小噪声分离变换成分，最小噪声分离变换成分所构成的矩阵表示为ψ；

g、计算最小噪声分离成分的谐波周期ψ(ω)；

h、选取谐波周期与潮汐应变分量谐波周期相对应的最小噪声分离成分进行重构；

i、计算去除潮汐应变分量的数据，并输出结果。

进一步地，所述步骤a包括：选取一组钻孔应变分钟值数据，将钻孔应变数据按照时间序列表示为：

x1＝[x1(1)，x1(2)，...，x1(1440)]，...，xn＝[xn(1)，xn(2)，...，xn(1440)]，得到样本矩阵y＝[x1，x2，x3，...xn]^t，样本矩阵y(n×1440)的表达式为：

进一步地，步骤b包括将制作好的样本矩阵y(n×1440)经过自适应窗宽滤波器进行滤波，估计出样本矩阵中的高频信号，自适应窗宽滤波器的表达式如下：

其中，o是滤波后数据，y是输入的样本数据；w为自适应窗宽，其表达式如下：

其中，wl为可选择的最小窗宽，wu为可选择的最大窗宽，δ(j)为样本矩阵y(n×1440)的二阶差分；

则估计的高频信号为：

h＝y-o(4)

其中，h表示的是h(n×1440)为估计的高频信号，y表示的是y(n×1440)为钻孔应变数据样本矩阵，o表示的是o(n×1440)为自适应窗宽滤波后的数据。

进一步地，步骤c包括计算估计的高频信号h(n×1440)的协方差矩阵ch(n×n)，协方差矩阵ch(n×n)元素γpq由公式(5)计算得到，

其中，和分别为第i行第p和第q分钟数据；和分别是n行数据的第p和第q分钟数据的均值；

对计算出的协方差矩阵ch(n×n)进行特征值分解，表达式如下：

式中，dh为特征值按照降序排列的对角矩阵；rh为特征值对角矩阵对应的特征向量矩阵，t为矩阵的转置。

进一步地，步骤d包括利用估计的高频信号h(n×1440)的协方差矩阵ch(n×n)来构造调整矩阵p，使得调整矩阵满足p^tchp＝e，其中e为单位矩阵，则调整矩阵的表达式如下：

p＝rhdh^-0.5(7)

其中，dh为特征值按照降序排列的对角矩阵；rh为特征值对角矩阵对应的特征向量矩阵。

进一步地，步骤e包括利用构造的调整矩阵p对钻孔应变数据样本矩阵的协方差矩阵cy(n×n)进行调整，其表达式为：cp＝p^tcyp，其中cp(n×n)为调整后的钻孔应变数据样本矩阵的协方差矩阵。

进一步地，步骤f包括对调整后的钻孔应变数据样本矩阵的协方差矩阵cp(n×n)做特征值分解，其表达式为：

其中，dp为是调整后的钻孔应变数据样本矩阵的协方差矩阵的特征值按照降序排列的对角阵，vp是调整后的钻孔应变数据样本矩阵的协方差矩阵的特征值对应的特征向量矩阵，最小噪声分离的变换矩阵r表示为：r＝pvp，其中p为调整矩阵，最小噪声分离变换成分所构成的矩阵ψ可表示为：

ψ＝yr＝y(pvp)＝[ψ1，ψ2，...，ψn]，n＝1，2，...，m，(9)

其中，y表示的是y(n×1440)为钻孔应变数据样本矩阵，ψ1，ψ2，...，ψn是最小噪声分离成分。

进一步地，步骤g包括利用傅立叶变换求取各最小噪声分离成分的谐波周期，其计算公式如下：

其中，ψn(t)为最小噪声分离成分，ψ(ω)为最小噪声分离成分的谐波周期。

进一步地，步骤h包括选取谐波周期与潮汐应变分量的谐波周期相对应的k个最小噪声分离成分进行重构，重构后的数据表达式如下：

ymnf＝ψr^-1＝[ψ1，ψ2，...，ψk，0，0，...，0]r^-1，1≤k≤m，(11)。

进一步地，步骤i包括利用钻孔应变数据样本矩阵y(n×1440)与最小噪声分离变换重构后的数据矩阵ymnf(n×1440)求去除潮汐应变分量的数据ys(n×1440)，其表达式如下：

ys＝y-ymnf(12)

其中，y(n×1440)是钻孔应变数据样本矩阵，ymnf(n×1440)是最小噪声分离变换重构后的数据矩阵。

本发明与现有技术相比，有益效果在于：

本发明可以判定数据内在的维数，利用钻孔应变数据中潮汐应变分量在频率上的特性，有针对性的去除潮汐应变分量，无需计算理论固体潮。

附图说明

图1为基于最小噪声分离的钻孔应变数据潮汐应变分量去除方法流程图；

图2为样本矩阵示意图；

图3为最小噪声分离成分图；

图4中(a)、(b)、(c)、(d)、(e)、(f)、(g)、(h)、(i)为9个与潮汐应变分量的谐波周期相对应的最小噪声分离成分的频谱图；

图5为利用选取的最小噪声分离成分重构后的数据图；

图6为去除了潮汐应变分量的时间序列数据曲线图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

以四川地区姑咱地震前兆监测台站的钻孔应变分钟值数据为例。该数据于2012年1月1日到2013年12月31日使用yry四分量钻孔应变仪测得。

a，录入钻孔应变数据s1，将钻孔应变每天的数据表示为：x1＝[x1(1)，x1(2)，...，x1(1440)]，...，x731＝[x731(1)，x731(2)，...，x731(1440)]，可以得到样本矩阵y＝[x1，x2，...，x731]^t，如图2所示。样本矩阵y(731×1440)的表达式为：

b，采用自适应窗宽滤波算法对钻孔应变数据样本矩阵进行高频信号估计是将制作好的样本矩阵y(731×1440)经过自适应窗宽滤波器进行滤波，估计出样本矩阵中的高频信号。自适应窗宽滤波器的表达式如下：

其中，o是滤波后数据，y是输入的样本数据。w为自适应窗宽，其表达式如下：

其中，wl为可选择的最小窗宽，wu为可选择的最大窗宽。δ(j)为样本矩阵y(731×1440)的二阶差分。

则估计的高频信号为：

h＝y-o(4)

其中，h(731×1440)是估计的高频信号，y(731×1440)是钻孔应变数据样本矩阵，o(731×1440)是自适应窗宽滤波后的数据。

c，计算估计的高频信号h(731×1440)的协方差矩阵ch(731×731)，协方差矩阵ch(731×731)元素γpq由公式(5)计算得到，

其中，和分别为第i行第p和第q分钟数据；和分别是n行数据的第p和第q分钟数据的均值。

对计算出的协方差矩阵ch(731×731)进行特征值分解，表达式如下：

式中，dh为特征值按照降序排列的对角矩阵；rh为特征值对角矩阵对应的特征向量矩阵，t为矩阵的转置。

d，构造调整矩阵，利用估计的高频信号h(731×1440)的协方差矩阵ch(731×731)来构造调整矩阵p，使得调整矩阵满足p^tchp＝e，其中e为单位矩阵。则调整矩阵的表达式如下：

p＝rhdh^-0.5(7)

其中，dh为特征值按照降序排列的对角矩阵；rh为特征值对角矩阵对应的特征向量矩阵。

e，利用构造的调整矩阵p对钻孔应变数据样本矩阵的协方差矩阵cy(731×731)进行调整，其表达式为：cp＝p^tcyp，其中cp(731×731)为调整后的钻孔应变数据样本协方差矩阵，cy(731×731)可利用公式(5)算出。

f，对调整后的钻孔应变数据样本矩阵的协方差矩阵cp(731×731)做特征值分解，其表达式为：

其中，dp为是调整后的钻孔应变数据样本矩阵的协方差矩阵的特征值按照降序排列的对角阵，vp是调整后的钻孔应变数据样本矩阵的协方差矩阵的特征值对应的特征向量矩阵。最小噪声分离的变换矩阵r可表示为：r＝pvp，其中p为调整矩阵。最小噪声分离变换成分所构成的矩阵ψ可表示为：

ψ＝yr＝y(pvp)＝[ψ1，ψ2，...，ψn]，n＝1，2，...，731，(9)

其中，y(731×1440)是钻孔应变数据样本矩阵，ψ1，ψ2，...，ψn是最小噪声分离成分。图3为最小噪声分离成分图。

g，利用傅立叶变换求取各最小噪声分离成分的谐波周期，其计算公式如下：

其中，ψn(t)为最小噪声分离成分，ψ(ω)为最小噪声分离成分的谐波周期。

h，选取谐波周期与潮汐应变分量的谐波周期相对应的9个最小噪声分离成分进行重构，图4为9个与潮汐应变分量的谐波周期相对应的最小噪声分离成分的频谱图。重构后的数据表达式如下：

ymnf＝ψr^-1＝[ψ1，ψ2，...，ψ9，0，0，...，0]r^-1，(11)

重构后的数据如图5所示。

i，求去除了潮汐应变分量的数据ys(731×1440)，其表达式如下：

ys＝y-ymnf

(12)

其中，y(731×1440)是钻孔应变数据样本矩阵，ymnf(731×1440)是最小噪声分离变换重构后的数据矩阵。ys(731×1440)是去除了潮汐应变分量的数据。去除了潮汐应变分量的时间序列数据如图6所示。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。