基于主成分分析的钻孔应变数据异常提取方法与流程

本发明涉及一种地震前兆观测数据异常检测方法，具体涉及一种基于主成分分析的钻孔应变数据异常提取方法。

背景技术：

地球科学是一门以观测和测量资料为基础的学科，对观测系统的不断完善是发展地球科学、减轻灾害的必由之路。地球表面变形和地壳内部的构造运动及其产生的各种地质灾害都与地壳应力作用密切相关，地壳应力状态的变化是导致断裂、褶皱乃至发生地震的最直接动因。钻孔应变观测通过对地层内部应变状态依时间连续变化的精细观测，发现和掌握地震应变前兆的(长)中期-短期-临震以及震后调整的时空分布与发展变化规律，是重要的地震前兆观测手段。钻孔应变观测凭借精度高、抗干扰能力强、多测项观测、仪器易于维护等特点成为了仅次于gps观测的前兆手段。

主成分分析方法理论上可以将原始信号中相对微弱的信号从较强干扰背景中识别出来，解决了钻孔应变观测中地壳变化微弱信号淹没在较强干扰背景中的问题。

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

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

邱泽华等(用小波-超限率分析提取宁陕台汶川地震体应变异常，2012)利用小波分解和改进的超限率分析方法对姑咱台站四分量钻孔应变数据进行研究，将细微的地震异常变化精细的提取出来。池顺良等(2013年庐山ms7.0地震的震前及临震应变异常，2013)利用庐山地震前后姑咱台站四分量钻孔应变观测资料，对地震同震信号进行分析。但到目前为止，尚未见运用主成分分析的方法对钻孔应变数据进行异常提取的报道。

技术实现要素：

本发明的目的就在于针对上述现有技术的不足，提供一种基于主成分分析的钻孔应变数据异常提取方法。

本发明的目的是通过以下技术方案实现的：

本发明利用主成分分析对钻孔应变数据进行异常提取，首先是将同一台站的钻孔应变数据序列进行应变换算，对换算后的数据进行预处理；将预处理后的钻孔应变数据构造成一个矩阵；并对每一天的矩阵进行主成分分析，以获取每个矩阵的特征值和特征向量；将得到的特征值与计算出来的特征向量角度与地震事件相对应，以得到异常的检测结果。通过本发明能够有效的利用主成分分析的方法对钻孔应变数据进行分析，根据钻孔应变各测项的相关性，对可能的地震前兆异常进行提取。

一种基于主成分分析的钻孔应变数据异常提取方法，包括以下步骤：

a、录入钻孔应变数据，并进行数据有效性的验证，“是”，进行下一步；

b、对钻孔应变数据进行应变换算；

c、数据预处理；

d、制作样本数据；

e、求样本数据的协方差矩阵；

f、求出协方差矩阵的特征值和特征向量；

g、求出特征值向量的角度；

h、输出特征值变化曲线；

i、输出特征向量角度变化曲线及统计直方图。

步骤a，所述的录入钻孔应变数据是选取一个地区不同台站的四分量钻孔应变数据，制作成按照分钟值采样的时间序列。数据有效性的验证是按照不同分量，数据可表示为s1,s2,s3,s4；根据由带孔平板弹性理论模型导出的圆孔径向变形与区域应力的关系式，对钻孔应变数据进行自洽分析，其关系式为：

s1+s3＝k(s2+s4)，(1)

选取自洽系数k在0.9以上的数据为有效数据。

步骤b，所述的对钻孔应变数据进行应变换算是根据公式(2)将四分量钻孔应变观测数据换算成两个剪应变s13、s24和一个面应变sa，

步骤c，所述的数据预处理是对数据进行调和分析，用来去除数据的周期项；其调和分析函数表达式为：

其中，a0为时间序列的直流分量，m为谐波的次数，系数am,bm是权重因子，表示各次谐波对总序列的贡献；原始数据减去拟合后的周期数据就是我们想要的预处理后数据。

步骤d，所述的制作样本数据是将预处理后台站每天的数据按时间序列表示为z1＝[z1(1),z1(2),...,z1(1440)],…,zn＝[zn(1),zn(2),...,zn(1440)]，可以得到样本矩阵y＝[z1,z2,z3,...,zn]^t。样本矩阵y的表达式为：

步骤e，样本数据y(n×1440)其协方差矩阵cy(n×n)的元素γpq由公式(5)计算得到，

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

步骤f，求出协方差矩阵的特征值是对矩阵y的协方差矩阵cy进行特征值分解，

cy＝r∧r^t，(6)

式中，∧(n×n)为大小降序排列的特征值对角矩阵；r(n×1)为与特征值对角矩阵相对应的特征向量矩阵。特征值表示为λ1,λ2,...,λn(λ1＞λ2＞...＞λn)。

步骤g，所述求出特征向量角度是设特征向量为r＝[v1v2...vn]^t，由图2所示，特征向量角度为：

与现有技术相比，本发明的有益效果在于：本发明基于主成分分析的钻孔应变数据异常提取方法，利用主成分分析中的特征值和特征向量角度分别将地壳的微弱变化表征出来；实现了在有较强背景干扰的情况下对钻孔应变数据异常的精确提取。

附图说明

图1为基于主成分分析的钻孔应变数据异常提取方法流程图；

图2为特征向量角度示意图；

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

图4为第一主成分特征值变化曲线；

图5为第一主成分特征向量角度变化示意图；

图6a为全部时间段的特征向量角度三维直方图；

图6b为无震期间特征向量角度三维直方图；

图6c为震前5个月地震前兆异常数据的特征向量角度三维直方图；

图6d为地震前后异常数据的特征向量角度三维直方图。

具体实施方式

下面结合实施实例对本发明进行详细说明。

针对芦山地震，以四川地区姑咱地震前兆监测台站的钻孔应变分钟值数据为例。姑咱台站与芦山地震位置如图3所示。该数据于2012年1月1日至2013年12月31日使用yry四分量钻孔应变仪测得。

a、录入姑咱台站2012年1月1日至2013年12月31日钻孔应变分钟值时间序列，按照北南分量、东西分量、北东分量、北西分量的顺序将数据分别记为s1、s2、s3、s4；根据由带孔平板弹性理论模型导出的圆孔径向变形与区域应力的关系式，对钻孔应变数据进行自洽分析，其关系式为：

s1+s3＝k(s2+s4)，(1)

计算出的k值，由k≥0.9可知姑咱台站2012年1月1日至2013年12月31日钻孔应变分钟值数据有效。

b、对钻孔应变数据进行应变换算是根据公式(2)将四分量钻孔应变观测数据换算成两个剪应变s13、s24和一个面应变sa；

c、数据预处理，即对换算后的数据进行调和分析，去掉应变数据的周期响应；调和分析函数表达式为：

其中a0为时间序列的直流分量，m为谐波的次数，系数am,bm是权重因子，表示各次谐波对总序列的贡献。原始数据减去拟合后的周期数据就是我们想要的预处理后数据。

d、将预处理后的钻孔应变每天的数据表示为

zn＝[zn(1),zn(2),...,zn(1440)],n＝1,2,3，

可以得到样本矩阵y＝[z1,z2,z3]^t。样本矩阵y的表达式为：

e、求样本数据y的协方差矩阵，样本数据y(3×1440)其协方差矩阵cy(3×3)的元素γpq由公式(5)计算得到，

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

f、对矩阵y的协方差矩阵cy进行特征值分解，

cy＝r∧r^t，(6)

式中，∧(3×3)为大小降序排列的特征值对角矩阵；r为与特征值对角矩阵相对应的特征向量矩阵，特征值表示为λ1,λ2,λ3(λ1＞λ2＞λ3)。

g、求出特征向量角度；以第一主成分特征向量为例设r1＝[v1,v2,v3]，根据图2所示，由公式(7)求出

第二、第三主成分特征向量同理可求。

h、得到特征值变化曲线，如图4所示，图中，实线为第一主成分特征值，---虚线为芦山地震时刻，矩形框中的为提取到的地震相关异常，可以看出特征值很好的表征了地震相关异常。

i、得到特征向量角度变化图及统计直方图，如图5，图6a-图6d所示。图5中空心圆点为特征向量角度，---虚线为芦山地震时刻，矩形框中的为提取到的地震相关异常。图6a-图6d为第一主成分特征向量角度变化三维直方图，从图6d可以看出，在异常部分的特征向量角度有明显的聚集现象，很好的提取到了地震前兆异常。