CN110618458B - 一种基于ica-ra的钻孔应变数据的分频段级联校正方法 - Google Patents

Abstract

本发明涉及钻孔应变领域，具体地而言为一种基于ICA‑RA的钻孔应变数据的分频段级联校正方法，目的是去除外界环境干扰引起的应变变化，首先是将同一观测站点的钻孔应变数据、气温数据、气压数据、降水数据和计算的理论固体潮数据进行4层小波分解；在分解后的5个频段中，逐频段进行多因素与应变的相关性分析；利用ICA‑RA级联校正方法，在各个频段内找到与钻孔应变数据相关性最大的影响因素。引入独立成分分析(ICA)方法，计算由于水位变化引起的应变变化并消除；对于非是气温、气压、和固体潮这类数据，利用回归分析(RA)方法，计算对应频段内主要因素产生的影响并校正。通过本发明能够有效的对钻孔应变数据各频段进行分析，对真实的应变情况有着重要的校正作用。

Description

一种基于ICA-RA的钻孔应变数据的分频段级联校正方法

技术领域

本发明涉及钻孔应变领域，具体地而言为一种基于ICA-RA的钻孔应变数据的分频段级联校正方法，目的是去除外界环境干扰引起的应变变化。

背景技术

国家“十五”期间，中国首创的四分量钻孔应变仪观测点发展较快，全国观测站点增加到50多个。在数种主流大地测量观测方法(测震、GPS和钻孔应变观测)中，钻孔应变仪填补了测震和GPS的测量盲区，在数月至数小时时段上观测精度最高。然而经过几十年的观测实践中发现，影响钻孔应变观测变化的因素很多，这就给提取有效的地下应变信息增加了难度。一般认为气象变化不是钻孔应变的观测的目标。消除掉这种气象影响是提高数据准确性的重要途径。中国地震局“十五”以来建设的钻孔应变观测点，很多都配备气象三要素(气温、气压和降雨量)观测，为后来的数据分析提供了必要的参考资料。

要消除这种影响，首先必须把这种影响的规律或者特征搞清楚。对于气温，一方面，温度和地面向深部的衰减很迅速；另一方面，温度向深部传导也很缓慢。所以，温度对于应变的影响是长周期的，一般可忽略不计。与气温不同，气压变化的主要影响是相对短周期的，且影响明显。周龙寿等(地壳应变场对气压短周期变化的响应，2008)和张凌空等(不同周期气压波对钻孔体应变仪观测结果的影响，2009)针对体应变仪，研究了不同周期的气压对应不同周期应变的影响。但需指出，目前对四分量钻孔应变仪观测数据的分析还很匮乏。在气象变化影响中，降水对钻孔应变的影响是最复杂的，对这个问题深入研究也很少。除去气象影响之外，钻孔应变数据还清晰的记录着固体潮汐的作用。

发明内容

本发明所要解决的技术问题在于提供一种基于ICA-RA的钻孔应变数据的分频段级联校正方法，解决四分量钻孔应变仪观测数据干扰引起的应变变化问题。

本发明是这样实现的，

一种基于ICA-RA的钻孔应变数据的分频段级联校正方法，该方法包括：

a、录入钻孔应变数据，录入同站点的气温、气压和水位数据共4组；

b、应用滑动平均法计算步骤a录入的4组数据的小时值数据；

c、计算同一站点位置的等时间跨度的理论固体潮数据；

d、对步骤b的4组数据的小时值数据和步骤c的理论固体潮数据，共5组数据均进行4层小波分解并重构，将数据分解到5个不同频段内，其中5个不同频段包括低频段作为第一频段，以及细节频段作为第二频段；

e、选取N＝1频段作为第一频段时的各因素作为初始状态；

f、计算各因素与对应频段应变的相关系数，并按降序排列；

g、判断是否存在一个相关系数大于0.8，是，进行步骤h；否，进行步骤p；

h、判断应变是否受水位数据影响最大，是，下一步；否，跳至步骤j，判断气温；

i、利用独立成分分析方法计算水位变化引起的应变变化并校正；

j、判断应变是否受气温数据影响最大，是，下一步；否，跳至步骤l，判断气压；

k、利用回归分析计算气温变化引起的应变变化并校正；

l、判断应变是否受气压数据影响最大，是，下一步；否，跳至步骤n，判断固体潮；

m、利用回归分析计算气压变化引起的应变变化并校正；

n、判断应变是否受固体潮数据影响最大，是，下一步；否，跳至步骤p；

o、利用回归分析计算固体潮变化引起的应变变化并校正；

p、判断是否最后一层已经校正结束，是，下一步；否，频段数N加1进入细节频段并重复步骤f～o；

q、叠加各频段内剩余的应变。

进一步地，步骤a所述的录入钻孔应变数据，录入同站点的气温、气压和水位数据具体为选取一个台站的四分量钻孔应变数据、气温数据、气压数据和水位数据，制作成按照分钟值采样的时间序列。

进一步地，步骤b所述的计算4组数据的小时值数据具体包括将分钟值采样的时间序列转换小时值，采用滑动窗平均法，其中第t个窗内的表达式如公式(1)，

n为窗长，i为分钟值采样点数，t代表小时时间序列；

将第t个窗内集合为：A_h＝[A_h(1),A_h(2),...,A_h(t)] (2)；

应变、气温、气压和水位数据代入到(1)和(2)式，转换后的应变、气温、气压和水位数据分别表示为S_h,T_h,P_h,W_h。

进一步地，步骤c所述的计算同一站点位置的等时间跨度的理论固体潮数据为计算序列长度为t的小时值理论固体潮，记为ET_h。

进一步地，步骤d中利用一维离散小波分析将每组数据分解4层并重构到5个不同频段中。

进一步地，步骤f所述的计算各因素与对应频段应变的相关系数，具体包括：选取时间范围(-t,+t)，逐步尝试所有可能的滞后时间，计算钻孔应变数据和另一组带滞后项的因素数据的相关系数，以使相关系数达到最大的滞后时间t_d为准。

进一步地，步骤i是利用独立成分分析计算水位变化引起的应变具体包括：

设m个信号源S＝[s₁,s₂,...,s_m]，n个观测数据，记为X＝[x₁,x₂,...,x_m],则

X＝AS (3)

式中A为m×n的混叠矩阵，A和S为未知，根据S的先验知识和多个观测信号X估计A和S；

寻找一个分离矩阵W，使Y＝WX,得到原信号的估计：

式中

表示对S的估计；

选取观测矩阵X，构造X，将X进行中心化，使它的均值为0；

再进行白化/球化处理，使观测矩阵X投影到新的子空间上变成白化向量Z：

Z＝W_nX (6)

其中，W_n是白化矩阵；

选择需要估计的分量个数m，迭代次数p＝1；

选择一个随机的初始权矢量W；

通过公式(7)为W_p赋值：

其中，

通过公式(6)更新W_p值，判断是否收敛，若W_p不收敛，转公式(4)：

p＝p+1,若p≤m，转公式(3)；

由公式(3)解得对原信号的估计值Y，Y为由于水位变化引起的应变变化；

用对应频段的应变减去水位引起的应变为校正后的应变。

进一步地，在相应频段按照影响应变的最主要因素，构造该因素的一元线性回归模型：

Yn_s＝β₀+β×Yn_x+ε

其中，n为1～4，x为T、P、或ET，其中T表示气温，P表示气压，ET表示固体潮，Yn_S是由于Yn_x引起的低频或细节频段变化，ε是随机误差项，β₀,β是模型参数，用b₀,b分别表示β₀,β通过样本数据的最小二乘估计，于是得到估计的一元线性回归方程为

其中，b₀是估计得回归直线在y轴上的截距，b是直线的斜率，

是Yn_S的估计值，表示由于因素的变动而引起的应变变化；

将影响去除，表示为

Y′n_S就是该频段内校正后的应变变化。

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

本发明公开的基于ICA-RA的钻孔应变数据的分频段级联校正方法，第一次提出在分频段分析影响钻孔应变数据的多可能因素，并通过相关性最大的判据确定了各频段影响应变变化的主要因素；第一次利用独立成分分析(ICA)的方法分析降水对应变的影响，并成功分离了由于降水引起的应变变化。第一次利用ICA-RA级联方法分析单一频段内的多影响因素。本发明根据不同影响因素的不同特征，提出了一种去除引起应变变化的因素影响的处理办法，有效的还原了钻孔应变的真实变化情况。

附图说明

图1是基于ICA-RA的钻孔应变数据的分频段级联校正方法流程图；

图2是姑咱地震前兆监测台站与芦山地震震中位置示意图；

图3是2012年1月1日到2013年6月30日的钻孔应变数据(a)，气温数据(b)，气压数据(c)，水位数据(d)和理论固体潮数(e)；

图4是在低频段去除水位影响前后的钻孔应变变化，其中，(a)含水位影响，(b)为去除含水位影响；

图5是在细节频段第一频段去除气压影响前后的钻孔应变变化(a)为含气压影响，(b)为去除气压影响；

图6是在细节频段第二频段和第三频段去除固体潮影响前后的钻孔应变变化(a)含固体潮，(b)为去除固体潮；

图7是去除以上因素影响的校正后的钻孔应变变化情况(a)为原始钻孔应变数据，(b)校正后的钻孔应变数据。

具体实施方式

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

参见图1，一种基于ICA-RA的钻孔应变数据的分频段级联校正方法，该方法包括：

a、录入钻孔应变数据，录入同站点的气温、气压和水位数据共4组；

b、应用滑动平均法计算步骤a录入的钻孔应变数据、同站点的气温、气压和水位数据的小时值数据；

c、计算同一站点位置的等时间跨度的理论固体潮数据；

d、对步骤b的4组数据的小时值数据和步骤c的理论固体潮数据，共5组数据均进行4层小波分解并重构，将数据分解到5个不同频段内，其中5个不同频段包括低频段作为第一频段，以及细节频段作为第二频段；

e、选取N＝1频段作为第一频段时的各因素作为初始状态；

f、计算各因素与对应频段应变的相关系数，并按降序排列；

g、判断是否存在一个相关系数大于0.8，是，进行步骤h；否，进行步骤p；

h、判断应变是否受水位数据影响最大，是，下一步；否，跳至步骤j，判断气温；

i、利用独立成分分析方法计算水位变化引起的应变变化并校正；

j、判断应变是否受气温数据影响最大，是，下一步；否，跳至步骤l，判断气压；

k、利用回归分析计算气温变化引起的应变变化并校正；

l、判断应变是否受气压数据影响最大，是，下一步；否，跳至步骤n，判断固体潮；

m、利用回归分析计算气压变化引起的应变变化并校正；

n、判断应变是否受固体潮数据影响最大，是，下一步；否，跳至步骤p；

o、利用回归分析计算固体潮变化引起的应变变化并校正；

p、判断是否最后一层已经校正结束，是，下一步；否，频段数N加1进入细节频段并重复步骤f～o；

q、叠加各频段内剩余的应变。

步骤a录入钻孔应变数据，录入同站点的气温、气压和水位数据具体为选取一个台站的四分量钻孔应变数据、气温数据、气压数据和水位数据，制作成按照分钟值采样的时间序列。这里用同站点的水位数据数据代替降雨量数据，因为同孔中的水位数据既包含了降水信息，也包含了地下水变化的信息。钻孔应变、气温、气压和水位数据可分别表示为S_m,T_m,P_m,W_m。

步骤b所述的计算4组数据的小时值数据具体包括将分钟值采样的时间序列转换小时值，采用滑动窗平均法，其中第t个窗内的表达式如公式(1)，

n为窗长，i为分钟值采样点数，t代表小时时间序列；

将第t个窗内集合为：A_h＝[A_h(1),A_h(2),...,A_h(t)] (2)；

应变、气温、气压和水位数据代入到(1)和(2)式，转换后的应变、气温、气压和水位数据分别表示为S_h,T_h,P_h,W_h。

步骤c所述的计算同一站点位置的等时间跨度的理论固体潮数据为计算序列长度为t的小时值理论固体潮，记为ET_h。

步骤d中利用一维离散小波分析将每组数据分解4层并重构到5个不同频段中。用小波方法分解观测时间序列，可以帮助我们分析不同频段中的影响因素和观测应变的变化情况，以及两者关系。以应变数据为例，应用db4小波基对S_h进行4层分解，得到A1_S,D1_S,D2_S,D3_S,D4_S分别是重构的低频段的应变，第一到第四层细节频段的应变。同样的，可以得到重构的低频气温和一到四层的细节频段的气温，记为A1_T,D1_T,D2_T,D3_T,D4_T，重构的低频气压和一到四层的细节频段的气压，记为A1_P,D1_P,D2_P,D3_P,D4_P，重构的低频水位和一到四层的细节频段的水位，记为A1_W,D1_W,D2_W,D3_W,D4_W，重构的低频固体潮和一到四层的细节频段的固体潮，记为A1_ET,D1_ET,D2_ET,D3_ET,D4_ET。

步骤f所述的计算各因素与对应频段应变的相关系数，即利用相关系数计算分析5组数据间的相关性。计算相关性的关键是确定恰当的滞后时间t_d，单位是小时。具体方法是，选取适当的时间范围(-t,+t)，逐步尝试所有可能的滞后时间，计算钻孔应变数据和另一组带滞后项的因素数据的相关系数，以使相关系数达到最大的滞后时间t_d为准，此时的两组数据相关性显著，可判定为该频段应变数据的主要影响因素，滞后时间为t_d。以应变数据和气压数据的第一层为例，应变数据为参考数据，两者的相关系数为

步骤g是判断是否存在一个相关系数大于0.8，是，进行步骤h，否，进行步骤p。步骤f的计算会产生4个相关系数。需要判断是否存在相关系数r在0.8以上的两组数据，认为此时是有效相关，即有影响。如果不存在任一r＞0.8，认为该频段应变没有收到外界干扰，校正时不考虑该层。然后进行步骤p。另外，此时最好将4个相关系数按降序排列，以便步骤h判断。

步骤h是判断该频段应变是否受水位数据影响最大，即判断在步骤f中的4个相关系数中，最大相关是否是应变与水位数据产生，是，进行下一步；否，跳至步骤j，判断气温因素；

步骤i是利用独立成分分析计算水位变化引起的应变，由于是降水的特殊性，有及时响应和延时响应，所以水位对应变的影响不是简单的线性相关可以解释的。这里我们引入盲源分离的思想，利用独立成分分析(Independent Component Analysis，ICA)的方法，将由于水位变化引起的应变变化分离出来。基本原理是利用信号的高阶统计特性，对多个观测信号进行盲源分离，估计出原始信号的近似值。

具体包括：

设m个信号源S＝[s₁,s₂,...,s_m]，n个观测数据，记为X＝[x₁,x₂,...,x_m],则

X＝AS (3)

式中A为m×n的混叠矩阵，A和S为未知，根据S的先验知识和多个观测信号X估计A和S；

寻找一个分离矩阵W，使Y＝WX,得到原信号的估计：

式中

表示对S的估计；

本发明采用基于负熵最大的fastICA算法，负熵定义如下：

N_g(Y)＝{E[g(Y)]-E[g(Y_gauss)]}² (5)

其中，E[·]为均值运算，g(·)为非线性函数，Y_gauss是与Y具有相同方差的高斯随机变量，以步骤d后的低频段应变与水位为例，fastICA算法流程如下：

(1)选取观测矩阵X，构造X，将X进行中心化，使它的均值为0；

再进行白化/球化处理，使观测矩阵X投影到新的子空间上变成白化向量Z：

Z＝W_nX (6)

其中，W_n是白化矩阵；

(2)选择需要估计的分量个数m，迭代次数p＝1；

(3)选择一个随机的初始权矢量W；

(4)通过公式(7)为W_p赋值：

其中，

(5)通过公式(6)更新W_p值，判断是否收敛，若W_p不收敛，转公式(4)：

(6)p＝p+1,若p≤m，转公式(3)；

(7)由公式(3)解得对原信号的估计值Y，Y为由于水位变化引起的应变变化；

用对应频段的应变减去水位引起的应变为校正后的应变；以第一频段应变与水位为例，即A1′_S＝A1_S-Y。

在相应频段按照影响应变的最主要因素，构造该因素的一元线性回归模型：

Yn_s＝β₀+β×Yn_x+ε (8)

其中，n为1～4，x为T、P、或ET，其中T表示气温，P表示气压，ET表示固体潮，Yn_S是由于Yn_x引起的低频或细节频段变化，ε是随机误差项，β₀,β是模型参数，用b₀,b分别表示β₀,β通过样本数据的最小二乘估计，于是得到估计的一元线性回归方程为

其中，b₀是估计得回归直线在y轴上的截距，b是直线的斜率，

是Yn_S的估计值，表示由于因素的变动而引起的应变变化；

将影响去除，表示为

Y′n_S就是该频段内校正后的应变变化。

实施例

针对芦山地震，以四川地区地震前兆监测台站中姑咱台的钻孔应变数据为例。姑咱台站与芦山地震震中位置如图2所示。该数据由YRY四分量钻孔应变仪测得，采样一分钟一次，时间段为2012年1月1日到2013年6月30日。

步骤一、录入姑咱台站钻孔应变数据，录入同站点的气温、气压和水位数据是选取一个台站的四分量钻孔应变数据、气温数据、气压数据和水位数据，制作成按照分钟值采样的时间序列。这里用同站点的水位数据代替降雨量数据，因为其与钻孔应变数据同孔，该点水位数据既包含了降水信息，也包含了地下水变化的信息。钻孔应变、气温、气压和水位数据可分别表示为S_m,T_m,P_m,W_m，共计四组数据。

步骤二、计算4组数据的小时值数据，将分钟值采样的时间序列转换小时值，采用滑动窗平均法，其中第t个窗内的表达式如公式1，以应变数据为例：

n为窗长，这里取60,i为分钟值采样点数。转换后的应变、气温、气压和水位数据可分别表示为S_h,T_h,P_h,W_h,以应变数据为例，S_h＝[S_h(1),S_h(2),...,S_h(t)]，t代表小时。

步骤三、计算姑咱台的2012年1月1日到2013年6月30日的理论固体潮数据，记为ET_h。五组数据分别为钻孔应变数据、同站点的气温、气压和水位数据的小时值数据以及理论固体潮数据如图3所示。

步骤四、对5组数据均进行4层小波分解是利用一维离散小波分析，将每组数据分解并重构到不同频带上。以应变数据为例，应用db4小波基对S_h进行4层分解，得到A1_S,D1_S,D2_S,D3_S,D4_S分别是重构的低频段的应变，第一到第四层细节频段的应变。同样的，可以得到重构的低频气温和一到四层的细节频段的气温，记为A1_T,D1_T,D2_T,D3_T,D4_T，重构的低频气压和一到四层的细节频段的气压，记为A1_P,D1_P,D2_P,D3_P,D4_P，重构的低频水位和一到四层的细节频段的水位，记为A1_W,D1_W,D2_W,D3_W,D4_W，重构的低频固体潮和一到四层的细节频段的固体潮，记为A1_ET,D1_ET,D2_ET,D3_ET,D4_ET。

步骤五、选取每种影响因素的第一频段数据即低频段数据，即选取N＝1频段时的各个因素,即选取A1_S,A1_T,A1_P,A1_W,A1_ET作为初始状态，以备下一步计算使用。

步骤六、计算各因素与对应频段应变的相关系数，即利用相关系数计算分析5组数据间的相关性。

计算相关性的关键是确定恰当的滞后时间t_d，单位是小时。具体方法是，选取适当的时间范围(-120,+120)，逐步尝试所有可能的滞后时间，计算钻孔应变数据和另一组带滞后项的数据的相关系数，以使相关系数达到最大的滞后时间t_d为准，此时的两组数据相关性显著，判定该因素是应变数据的主要影响，滞后时间为t_d。以应变数据和气压数据的第一频段为例，两者的相关系数为

步骤七、判断该频段是否存在一个相关系数大于0.8，是，进行步骤八，否，进行步骤十六。步骤六的计算会产生4个相关系数。需要判断是否存在相关系数r在0.8以上的两组数据，认为此时是有效相关。如果不存在任一r＞0.8，认为该频段应变没有收到外界干扰，校正时不考虑该层。然后频段数N加1，进行步骤十六。研究时间段内的各个频段内，4组数据与钻孔应变数据相关分析的对应滞后时间和相关系数如表1。

表1为钻孔应变，气温，气压，水位和理论固体潮做小波分解后，各个频段内，5组数据与钻孔应变数据相关分析的对应滞后时间和相关系数表。表示为(A,B)，A为最优的滞后时间td，B是此时的相关系数r。

表1

	气温	气压	水位	固体潮
					低频段	(-95,-0.9102)	(-50,-0.5584)	(120,-0.9561)	(80,-0.0925)
细节频段第1层	(80,0.8019)	(50,-0.8407)	(-15,0.5767)	(-80,0.7674)
					细节频段第2层	(120,-0.7033)	(25,-0.6719)	(0,-0.6323)	(-80,0.9396)
细节频段第3层	(-120,-0.4885)	(20,0.6251)	(0,-0.3048)	(-80,0.8258)
					细节频段第4层	(70,-0.0829)	(-50,0.2460)	(0,-0.2084)	(-80,-0.3389)

另外，此时优选地将4个相关系数按降序排列，以便步骤九判断。

步骤八、判断该频段应变是否受水位数据影响最大，即判断在步骤f中的4个相关系数，最大相关是否是应变与水位数据产生，这里第一频段应变与降水的相关系数为-0.9561，是最大，所以进行下一步。

步骤九、对于水位因素，利用独立成分分析分离其影响并去除。由于是降水的特殊性，有及时响应和延时响应，所以水位对应变的影响不是简单的线性相关可以解释的。这里引入盲源分离的思想，利用独立成分分析(Independent Component Analysis，ICA)的方法，将由于水位变化引起的应变变化分离出来。基本原理是利用信号的高阶统计特性，对多个观测信号进行盲源分离，估计出原始信号的近似值。

设m个信号源S＝[s₁,s₂,...,s_m]，n个观测数据，记为X＝[x₁,x₂,...,x_m],则

X＝AS (3)

式中A为m×n的混叠矩阵，A和S是未知的，ICA需要根据S的先验知识和多个观测信号X来估计A和S。因此，需要找到一个分离矩阵W，使Y＝WX,从而得到原信号的估计：

式中

表示对S的估计。

本发明采用基于负熵最大的fastICA算法，负熵定义如下：

N_g(Y)＝{E[g(Y)]-E[g(Y_gauss)]}² (5)

其中，E[·]为均值运算，g(·)为非线性函数，Y_gauss是与Y具有相同方差的高斯随机变量，以步骤d后的低频段应变与水位为例，fastICA算法流程如下：

(1)选取观测矩阵X，构造X＝[A1_S,A1_W]，将X进行中心化，使它的均值为0；再进行白化/球化处理，使观测矩阵X投影到新的子空间上变成白化向量Z。

Z＝W_nX (6)

其中，W_n是白化矩阵；

(2)选择需要估计的分量个数m,这里m＝2。迭代次数p＝1；

(3)选择一个随机的初始权矢量W；

(4)通过公式(7)为W_p赋值：

其中，

(5)通过公式(6)更新W_p值，判断是否收敛，若W_p不收敛，转(4)；

(6)p＝p+1,若p≤m，转(3)；

(7)由公式(3)解得对原信号的估计值Y。

这里Y就是由于水位变化引起的应变变化。用低频层的应变减去水位引起的应变就是校正后的应变，即A1′_S＝A1_S-Y。水位因素去除前后如图4。

由于第一频段中，水位是影响应变的最大因素，可直接跳到步骤十六，进行下一频段的判断，即回到步骤六。

第二频段中，步骤六与步骤七与第一频段判断方法相同，不重复赘述。步骤八中，判断该频段应变是否受水位数据影响最大。由表1可知，这里为否，该频段影响最大的是气压，相关系数为-0.8407，可直接跳到步骤十三。

步骤十三、利用回归分析(RA)方法计算气压变化引起的应变并校正。

第二频段中，即细节频段第一层中，气压是影响应变的最主要因素，那么可以构造如下的一元线性回归模型

D1_S＝β₀+β×D1_P+ε (12)

模型的线性部分反映了由于D1_P的变化而引起的D1_S的变化，ε_i是随机误差项，β₀,β是模型参数，b₀,b分别是β₀,β通过样本数据的最小二乘估计，于是得到估计的一元线性回归方程为

其中，b₀是估计得回归直线在y轴上的截距，b是直线的斜率，

是D1_S的估计值，即表示在第二层细节频段内，由于气温的变动而引起的气压变化。那么接下来，就可以将这种影响去除，表示为

D1′_S就是该频段内校正气压后的应变变化。气压因素去除前后如图5。

由于第二频段中，气压是影响应变的最大因素，可直接跳到步骤十六，进行下一频段的判断，即回到步骤六。

第三频段中，步骤六与步骤七与前两个频段相同，不重复赘述。步骤八中，判断该频段应变是否受水位数据影响最大。由表1可知，这里为否，该频段影响最大的是固体潮，相关系数为0.9396，可直接跳到步骤十五。

步骤十五是利用回归分析(RA)计算固体潮变化引起的应变并校正。校正方法与步骤十三一致，不重复赘述。该步校正表示为

D2′_ET就是该频段内校正固体潮后的应变变化。

由于第三,频段中，固体潮是影响应变的最大因素，可直接跳到步骤十六，进行下一频段的判断，即回到步骤六。

第四频段中，步骤六与步骤七与前三个频段相同，不重复赘述。步骤八中，判断该频段应变是否受水位数据影响最大。由表一可知，这里为否，该频段影响最大的是仍然固体潮，相关系数为0.8258，可直接跳到步骤十五。

步骤十五是利用回归分析(RA)计算固体潮变化引起的应变并校正。校正方法与步骤十三一致，不重复赘述。该步校正表示为

D3′_ET就是该频段内校正固体潮后的应变变化。第三层和第四层固体潮因素去除前后如图6。

由于第四频段中，固体潮是影响应变的最大因素，可直接跳到步骤十六，进行下一频段的判断，即回到步骤六。

第五频段中，由表1可知，不存在任一r＞0.8，认为该频段应变没有收到外界干扰，校正时不考虑该频段。直接跳到步骤十六。

步骤十六、判断是否最后一频段已经校正结束。第五频段已经是最后一频段，进行步骤十七。

步骤十七、将剩余的应变叠加得到去除固体潮和气象因素影响的应变变化就是将上述步骤中各频段剩余的应变相加，认为校正完毕。图7是去除以上因素影响的真实的钻孔应变变化情况。

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

Claims (8)

1.一种基于ICA-RA的钻孔应变数据的分频段级联校正方法，其特征在于，该方法包括：

a、录入钻孔应变数据，录入同站点的气温、气压和水位数据共4组；

b、应用滑动平均法计算步骤a录入的4组数据的小时值数据；

c、计算同一站点位置的等时间跨度的理论固体潮数据；

d、对步骤b的4组数据的小时值数据和步骤c的理论固体潮数据，共5组数据均进行4层小波分解并重构，将数据分解到5个不同频段内，其中5个不同频段包括低频段作为第一频段，以及细节频段作为第二频段；

e、选取N＝1频段作为第一频段时的各因素作为初始状态；

f、计算各因素与对应频段应变的相关系数，并按降序排列；

g、判断是否存在一个相关系数大于0.8，是，进行步骤h；否，进行步骤p；

h、判断应变是否受水位数据影响最大，是，下一步；否，跳至步骤j，判断气温；

i、利用独立成分分析方法计算水位变化引起的应变变化并校正；

j、判断应变是否受气温数据影响最大，是，下一步；否，跳至步骤l，判断气压；

k、利用回归分析计算气温变化引起的应变变化并校正；

l、判断应变是否受气压数据影响最大，是，下一步；否，跳至步骤n，判断固体潮；

m、利用回归分析计算气压变化引起的应变变化并校正；

n、判断应变是否受固体潮数据影响最大，是，下一步；否，跳至步骤p；

o、利用回归分析计算固体潮变化引起的应变变化并校正；

p、判断是否最后一层已经校正结束，是，下一步；否，频段数N加1进入细节频段并重复步骤f～o；

q、叠加各频段内剩余的应变。

2.按照权利要求1所述的方法，其特征在于，步骤a所述的录入钻孔应变数据，录入同站点的气温、气压和水位数据具体为选取一个台站的四分量钻孔应变数据、气温数据、气压数据和水位数据，制作成按照分钟值采样的时间序列。

3.按照权利要求2所述的方法，其特征在于，步骤b所述的计算4组数据的小时值数据具体包括将分钟值采样的时间序列转换小时值，采用滑动窗平均法，其中第T个窗内的表达式如公式(1)，

w为窗长，i为分钟值采样点数，t代表小时时间序列；

将第T个窗内集合为：A_h＝[A_h(1),A_h(2),...,A_h(T)] (2)；

应变、气温、气压和水位数据代入到(1)和(2)式，转换后的应变、气温、气压和水位数据分别表示为S_h,T_h,P_h,W_h。

4.按照权利要求1所述的方法，其特征在于，步骤c所述的计算同一站点位置的等时间跨度的理论固体潮数据为计算序列长度为N的小时值理论固体潮，记为ET_h。

5.按照权利要求1所述的方法，其特征在于，

步骤d中利用一维离散小波分析将每组数据分解4层并重构到5个不同频段中。

6.按照权利要求1所述的方法，其特征在于，

步骤f所述的计算各因素与对应频段应变的相关系数，具体包括：选取时间范围(-a,+a)，逐步尝试所有可能的滞后时间，计算钻孔应变数据和另一组带滞后项的因素数据的相关系数，以使相关系数达到最大的滞后时间t_d为准。

7.按照权利要求1所述的方法，其特征在于，

步骤i是利用独立成分分析计算水位变化引起的应变具体包括：

设m个信号源S＝[s₁,s₂,...,s_m]，n个观测数据，记为X＝[x₁,x₂,...,x_n],则

X＝AS (3)

式中A为m×n的混叠矩阵，A和S为未知，根据S的先验知识和多个观测信号X估计A和S；

寻找一个分离矩阵W，使Y＝WX,得到原信号的估计：

式中

表示对S的估计；

选取观测矩阵X，将观测矩阵X进行中心化，使它的均值为0；

再进行白化/球化处理，使观测矩阵X投影到新的子空间上变成白化向量Z：

Z＝W_nX (6)

其中，W_n是白化矩阵；

选择需要估计的信号源个数m，迭代次数p＝1；

选择一个随机的分离矩阵W作为初始权矢量；

通过公式(7)为W_p赋值：

其中，

通过公式(6)更新W_p值，判断是否收敛，若W_p不收敛，转公式(4)：

p＝p+1,若p≤m，转公式(3)；

由公式(3)和公式(4)解得对原信号的估计值Y，Y为由于水位变化引起的应变变化；

用对应频段的应变减去水位引起的应变为校正后的应变。

8.按照权利要求1所述的方法，其特征在于，在相应频段按照影响应变的最主要因素，构造该因素的一元线性回归模型：

Yf_s＝β₀+β×Yf_x+ε

其中，f为1～4，x为T、P或ET，其中T表示气温，P表示气压，ET表示固体潮，Yf_S是由于Yf_x引起的低频或细节频段变化，ε是随机误差项，β₀,β是模型参数，用b₀,b分别表示β₀,β通过样本数据的最小二乘估计，于是得到估计的一元线性回归方程为

其中，b₀是估计得回归直线在y轴上的截距，b是直线的斜率，

是Yf_S的估计值，表示由于因素的变动而引起的应变变化；

将影响去除，表示为

Y′f_S就是该频段内校正后的应变变化。

Images