CN109214092A - 基于希尔伯特黄变换的四分量钻孔应变数据异常提取方法 - Google Patents

Abstract

本发明涉及四分量钻孔应变数据异常提取领域，特别涉及一种基于希尔伯特黄变换的四分量钻孔应变数据异常提取方法，包括：录入钻孔应变四分量数据和气温数据、气压数据、钻孔水位数据、钻孔气压数据，并对钻孔应变四分量数据进行应变换算为面应变S_a；对应变换算后的钻孔应变数据进行集合经验模态分解；计算分解后的各本征模函数的谐波周期和气温、气压、钻孔水位、钻孔气压干扰因素的谐波周期；去除与干扰因素谐波周期相对应的本征模函数；对余下的本征模函数进行希尔伯特变换，并计算瞬时能量；输出瞬时能量图。能有针对性的提取到钻孔应变数据影响因素的特征，为影响因素的排除及提取与地壳应变变化相关的异常信息的问题提供支持。

Description

基于希尔伯特黄变换的四分量钻孔应变数据异常提取方法

技术领域

本发明涉及四分量钻孔应变数据异常提取领域，特别涉及一种基于希尔伯特黄变换的四分量钻孔应变数据异常提取方法。

背景技术

钻孔应变观测是通过对地层内部应变状态依时间连续变化的精细观测，发现并掌握震前、震中和震后地壳应变变化的时空分布与发展变化规律，对地震的预测预报工作有着重要的意义。钻孔应变观测由于其精度高，很容易受到环境因素的干扰，使与地壳应变变化相关的信息淹没在这些干扰中。钻孔应变数据主要的干扰因素有气温、气压、钻孔水位(降水量)、钻孔气压。希尔伯特黄变换在理论上可以将复杂的信号分解成为有限个本征模函数(IMF)，所分解出的本征模函数包含了原信号的不同时间尺度的局部特征信息，可以解决排除干扰因素并提取到与地壳应变变化相关的异常信息的问题。

CN106918836A公开了一种基于主成分分析的钻孔应变数据异常提取方法，根据数据之间相关性的特点，利用主成分分析中的特征值和特征向量分别的将地壳的微弱变化表征出来。

CN202452947U公开了一种四分量钻孔应变仪观测系统，包括井上设置的数据存储网络传输设备、井下设置的主四分量应变仪及连接井上井下的信号线缆，还包括井下设置的辅助四分量应变仪。通过在观测系统内设置一主、一辅两个四分量应变仪，增加了系统的可靠性，避免了因一个装置损坏就使得整个系统瘫痪的情形。

邱泽华等(用小波-超限率分析提取宁陕台汶川地震体应变异常，2012)采用高通滤波的方法提取四分量钻孔应变数据的高频成分。并利用小波-超限率方法对宁陕台汶川地震提应变异常进行提取分析。池顺良等(汶川、鲁甸、康定地震前应变数据由自洽到失洽的转变与地震成核，2014)利用地震前后自恰系数k值的变化，对汶川、鲁甸、康定地震前四分量钻孔应变数据由自洽到失洽的转变的异常现象进行了分析。但到目前为止，尚未见运用希尔伯特黄变换的方法对四分量钻孔应变数据进行异常提取的报道。然后现有技术中无法实现在数据处理的过程中有针对性的排除钻孔应变数据的影响因素干扰并提取与地壳应变变化相关的异常信息。

发明内容

本发明所要解决的技术问题在于提供一种基于基于希尔伯特黄变换的四分量钻孔应变数据异常提取方法,其能有针对性的提取到钻孔应变数据影响因素的特征，为影响因素的排除及提取与地壳应变变化相关的异常信息的问题提供支持。

本发明是这样实现的，

一种基于希尔伯特黄变换的四分量钻孔应变数据异常提取方法，该方法包括：

a、录入钻孔应变四分量数据和气温数据、气压数据、钻孔水位数据、钻孔气压数据，并对钻孔应变四分量数据进行应变换算为面应变Sa；

b、对应变换算后的钻孔应变数据进行集合经验模态分解；

c、计算分解后的各本征模函数的谐波周期和气温、气压、钻孔水位、钻孔气压干扰因素的谐波周期；

d、去除与干扰因素谐波周期相对应的本征模函数；

e、对余下的本征模函数进行希尔伯特变换，并计算瞬时能量；

f、输出瞬时能量图。

进一步地，步骤a中将钻孔应变四分量数据记为S₁、S₂、S₃、S₄，气温数据记为T，气压数据记为P，钻孔水位数据记为L，钻孔气压数据记为P_z，根据公式(1)将钻孔应变四分量数据换算为面应变S_a，

S_a＝(S₁+S₂+S₃+S₄)/2 (1)。

进一步地，所述的对应变换算后的钻孔应变数据进行集合经验模态分解是对应变换算后的钻孔应变数据进行集合经验模态分解，其步骤如下：

步骤1，向面应变S_a中加入正态分布白噪声，其关系式如下所示：

Y＝S_a+w (2)

其中，Y是加入白噪声后的数据，w是加入的白噪声数据；

步骤2，找出Y中所有的极大值点和极小值点，并用三次样条插值函数分别拟合形成上、下包络，三次样条插值函数的表达式为：

S_i(x)＝a_ix³+b_ix²+c_ix+d_i，(i＝1，2，...，n-1) (3)

并求取上、下包络的均值并记做M，则有：H＝Y-M,H为数据残差。

步骤3，持续对数据残差H进行分解直到H满足在所有时间内局部零极点和过零点的数目相等或相差一个，并且在任意时刻局部上、下包络的均值为零时停止分解，此时记H＝IMF，IMF为分解得到的本征模函数；

重复步骤1-3，得到所有的本征模函数IMF_n。

进一步地，步骤c利用傅里叶变换对各本征模函数和气温、气压、钻孔水位、钻孔气压干扰因素进行谐波周期的求取，其计算公式如下：

其中，IMF_n为各本征模函数，ψ(ω)为本征模函数的频率。

进一步地，步骤d中去除与干扰因素谐波周期相对应的本征模函数包括将计算出的各本征模函数的谐波周期与气温、气压、钻孔水位、钻孔气压干扰因素的谐波周期进行比较，去除与干扰因素谐波周期相对应的本征模函数。

进一步地，步骤e中希尔伯特黄变换后的数据表示为：

其中，A是瞬时幅值，ω是瞬时频率，瞬时能量E通过式(6)计算：

E(t)＝∑_ωIMF_n ²(t，ω)，n＝1，2，...k.(6)。

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

本发明能将复杂的钻孔应变数据分解为有限个本征模函数，所分解出来的各本征模函数分量包含了原钻孔应变数据的不同时间尺度的局部特征信号，并且其频率具有一定的物理意义。与现有技术相比，这种方法是直观的、直接的、后验的和自适应的，可以有效的排除干扰因素并提取到与地壳应变变化相关的异常信息。

附图说明

图1为基于希尔伯特黄变换的四分量钻孔应变数据异常提取方法的流程图；

图2为钻孔应变面应变数据及影响因素曲线图；

图3 IMF分解图(a)及其相应的谐波周期图(b)；

图4为影响因素谐波周期图，其中a为气温、b为钻孔气压、c为气压、d为钻孔水位；

图5为第6个IMF分量的瞬时能量图。

具体实施方式

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

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

步骤a，录入钻孔应变四分量数据并记为S₁、S₂、S₃、S₄，气温数据记为T，气压数据记为P，钻孔水位数据记为L，钻孔气压数据记为P_z。根据公式(1)将钻孔应变四分量数据换算为面应变S_a，

S_a＝(S₁+S₂+S₃+S₄)/2 (1)

图2为钻孔应变面应变数据及影响因素曲线图。

步骤b，首先向面应变S_a中加入正态分布白噪声，其关系式如下所示：

Y＝S_a+w (2)

其中，Y是加入白噪声后的数据，w是加入的白噪声数据。

找出Y中所有的极大值点和极小值点，并用三次样条插值函数分别拟合形成上、下包络，三次样条插值函数的表达式为：

S_i(x)＝a_ix³+b_ix²+c_ix+d_i，(i＝1，2，...，n-1) (3)

并求取上、下包络的均值并记做M，则有：H＝Y-M，这里H为数据残差。

持续对数据残差H进行分解直到H满足在所有时间内局部零极点和过零点的数目相等或相差一个，并且在任意时刻局部上、下包络的均值为零时停止分解。此时记H＝IMF，IMF即为分解得到的本征模函数。

重复步骤1-3，得到所有的本征模函数IMF_n。我们选择将数据分解成为15个本征模函数，最后一个本征模函数记为数据的趋势项res。

步骤c，利用傅里叶变换对各本征模函数和气温、气压、钻孔水位、钻孔气压等干扰因素进行谐波周期的求取，其计算公式如下：

其中，IMF_n为各本征模函数，ψ(ω)为本征模函数的频率。利用公式(4)分别求取气温、气压、钻孔水位、钻孔气压等干扰因素的谐波周期。图3为分解的各IMF_n图以及其对应的的谐波周期，其中IMF分解图(a)及其相应的谐波周期图(b)。图4为求取的影响因素谐波周期图。

步骤d，将计算出的各本征模函数的谐波周期与气温、气压、钻孔水位、钻孔气压等干扰因素的谐波周期进行比较：IMF₁到IMF₅的谐波周期小于1小时，我们将其视为仪器噪声影响；IMF₇到IMF₁₀的谐波周期与气温、气压、钻孔气压的谐波周期相似，我们将其视为气温、气压、钻孔气压的影响；IMF₁₁到IMF₁₃的谐波周期与钻孔水位的谐波周期相似，我们将其视为钻孔水位的影响；IMF₁₄和res表征的是数据的长周期信号和整体趋势信号。经过排除我们将IMF₆视为与地壳应变变化相关的信号。

步骤e，利用希尔伯特黄变换，IMF₆可以表示为：

其中，A是瞬时幅值，ω是瞬时频率。瞬时能量E的计算公式如下：

E(t)＝∑_ωIMF₆ ²(t，ω) (6)

步骤f，输出瞬时能量图，图5为第6个IMF分量的瞬时能量图。如图5所示，除了在地震之前几天出现能量波动外，在12年10月末到12年12月末也出现了能量的异常现象。瞬时能量很好的表征了大地震前四分量钻孔应变数据的异常现象。

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

Claims (6)

1.一种基于希尔伯特黄变换的四分量钻孔应变数据异常提取方法，其特征在于，该方法包括：

a、录入钻孔应变四分量数据和气温数据、气压数据、钻孔水位数据、钻孔气压数据，并对钻孔应变四分量数据进行应变换算为面应变S_a；

b、对应变换算后的钻孔应变数据进行集合经验模态分解；

c、计算分解后的各本征模函数的谐波周期和气温、气压、钻孔水位、钻孔气压干扰因素的谐波周期；

d、去除与干扰因素谐波周期相对应的本征模函数；

e、对余下的本征模函数进行希尔伯特变换，并计算瞬时能量；

f、输出瞬时能量图。

2.按照权利要求1所述的方法，其特征在于，步骤a中将钻孔应变四分量数据记为S₁、S₂、S₃、S₄，气温数据记为T，气压数据记为P，钻孔水位数据记为L，钻孔气压数据记为P_z，根据公式(1)将钻孔应变四分量数据换算为面应变S_a，

S_a＝(S₁+S₂+S₃+S₄)/2 (1)。

3.按照权利要求1所述的方法，其特征在于，所述的对应变换算后的钻孔应变数据进行集合经验模态分解是对应变换算后的钻孔应变数据进行集合经验模态分解，其步骤如下：

步骤1，向面应变S_a中加入正态分布白噪声，其关系式如下所示：

Y＝S_a+w (2)

其中，Y是加入白噪声后的数据，w是加入的白噪声数据；

步骤2，找出Y中所有的极大值点和极小值点，并用三次样条插值函数分别拟合形成上、下包络，三次样条插值函数的表达式为：

S_i(x)＝a_ix³+b_ix²+c_ix+d_i，(i＝1，2，...，n-1) (3)

并求取上、下包络的均值并记做M，则有：H＝Y-M,H为数据残差。

步骤3，持续对数据残差H进行分解直到H满足在所有时间内局部零极点和过零点的数目相等或相差一个，并且在任意时刻局部上、下包络的均值为零时停止分解，此时记H＝IMF，IMF为分解得到的本征模函数；

重复步骤1-3，得到所有的本征模函数IMF_n。

4.按照权利要求1所述的方法，其特征在于，步骤c利用傅里叶变换对各本征模函数和气温、气压、钻孔水位、钻孔气压干扰因素进行谐波周期的求取，其计算公式如下：

其中，IMF_n为各本征模函数，ψ(ω)为本征模函数的频率。

5.按照权利要求1所述的方法，其特征在于，步骤d中去除与干扰因素谐波周期相对应的本征模函数包括将计算出的各本征模函数的谐波周期与气温、气压、钻孔水位、钻孔气压干扰因素的谐波周期进行比较，去除与干扰因素谐波周期相对应的本征模函数。

6.按照权利要求1所述的方法，其特征在于，步骤e中希尔伯特黄变换后的数据表示为：

其中，A是瞬时幅值，ω是瞬时频率，瞬时能量E通过式(6)计算：

E(t)＝∑_ωIMF_n ²(t，ω)，n＝1，2，...k. (6)。