一种基于S-K特征的钻孔应变数据异常提取方法与流程

本发明属于地震前兆观测数据异常检测领域，具体地来讲为一种基于s-k特征的钻孔应变数据异常提取方法。

背景技术

国的钻孔应变观测起步较早，特别是首创的水平四分量钻孔应变仪，设计理念具有国际领先水平。国家“十五”期间，这种四分量钻孔应变仪观测点发展较快，已经从原来屈指可数的几个观测点增加到50多个。其中一些观测点已经开始产出高质量的资料。中国地震局地震监测网络的首要目的是地震预测预报。地震预报研究的一个重要发展方向，是从“经验预报”向物理预报转变。连续地应变观测资料正是物理预报赖以进行的依据。地震预报的另一个发展方向，是解决短l临前兆问题。在数种主流大地测量观测方法(测震、gps和钻孔应变观测)中，钻孔应变仪也正是在数月至数小时时段上观测精度最高。除此之外，钻孔应变观测在避免地面干扰、少占用地等方面也具有优势。

cn106918836a公开了一种基于主成分分析的钻孔应变数据异常提取方法，首先是将同一台站的钻孔应变数据序列进行应变换算，对换算后的数据进行预处理；将预处理后的钻孔应变数据构造成一个矩阵；并对每一天的矩阵进行主成分分析，以获取每个矩阵的特征值和特征向量；将得到的特征值与计算出来的特征向量角度与地震事件相对应，以得到异常的检测结果。通过本发明能够有效的利用主成分分析的方法对钻孔应变数据进行分析，根据钻孔应变各测项的相关性，对可能的地震前兆异常进行提取。本发明钻孔应变数据异常提取方法，利用主成分分析中的特征值和特征向量角度分别将地壳的微弱变化表征出来；实现了在有较强背景干扰的情况下对钻孔应变数据异常的精确提取。

根据长期的观测资料和经验，一般认为短周期的高频信息在无其他干扰因素的条件下，多表现为一种服从正态分布的随机信号。从统计的角度分析，在去除趋势变化后，且无其他影响因素存在的条件下，观测资料变化值服从正态分布。偏度和峰度就是统计学中描述数据分布偏离正态分布的两个参数，基于这一特性，本发明提出了一种提取钻孔应变数据中异常的方法。到目前为止，尚未见有关于偏度和峰度的方法对钻孔应变数据进行异常提取的报道。

技术实现要素：

本发明所要解决的技术问题在于提供一种基于s-k特征的钻孔应变数据异常提取方法，能够有效的对钻孔应变数据进行分析，对可能的地震前兆异常进行提取。

本发明是这样实现的，一种基于s-k特征的钻孔应变数据异常提取方法，包括如下步骤：

a、对钻孔应变数据进行数据有效性的验证，若为有效，则进行下一步；

b、对钻孔应变数据进行面应变换算与调和分析去低频信息；

c、对处理后数据差分计算并制作样本数据；

d、计算每天样本数据的偏度和峰度；

e、将偏度和峰度画在s-k平面内；

f、输入s-k平面散点图；

g、定义一个地壳平静背景圆；

h、计算每天样本数据偏离背景的程度，定义为s-k偏差；

i、输出随时间变化的s-k偏差曲线。

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

进一步地，步骤b所述对钻孔应变数据进行应变换算式根据下式(2)将四分量钻孔应变观测数据换算成面应变sa：

sa(t)＝(s1(t)+s2(t)+s3(t)+s4(t))/2(2)。

进一步地，步骤b中所述调和分析用来去除数据中固体潮低频信息的周期项,其调和分析函数s(t)表达式为：

其中a0为时间序列的直流分量，m为谐波的次数，系数am,bm是权重因子，表示各次谐波对总序列的贡献。

进一步地，差分表达式为：

dsa(t)＝(sa(t+1)-s(t+1))-(sa(t)-s(t))，(4)

，按天分段制成n天的样本数据ds(n)，其中ds(n)表示差分后的面应变。

进一步地，计算每天样本数据的偏度和峰度，利用统计学上高阶累积量的算法计算，定义现有随机变量x，其均值是μ，方差是σ，则偏度和峰度表达式如下：

其中，随机变量x取自样本数据，sk为偏度，ku为峰度，正态分布的偏度为0，若数据分布是对称的，偏度为0，若偏度>0,则分布为右偏；若偏度<0，则分布为左偏，同时偏度的绝对值越大，说明分布的偏移程度越严重。

进一步地，定义一个地壳平静背景圆包括：s-k平面内所有的点中去掉一部分的极端点，剩余点求坐标参数的平均值作为圆点，坐标参数的标准差的平均值作为圆的半径。

本发明与现有技术相比，有益效果在于：本发明公开的基于偏度和峰度的钻孔应变数据异常提取方法，第一次从数据分布上分析钻孔应变数据，并第一次把偏度峰度和数据的物理意义对应，深入的了解了钻孔应变数据；第一次发现钻孔应变数据的偏度和峰度成抛物线关系，并定义了合适背景，计算了每天面应变与背景的随时间的偏差，提出了有效的提取分布异常的方法。

采用该方法能够对可能的地震前兆异常进行提取，根据以往发生地震前的数据分析，能否正确的提取地震前数据的异常。

附图说明

图1是基于s-k特征的钻孔应变数据异常提取方法流程图。

图2是姑咱地震前兆监测台站与汶川地震震中位置示意图。

图3是2007年1月1日到2009年7月1日的钻孔应变面应变的差分数据。

图4是2007年1月1日到2009年7月1日的钻孔应变面应变的差分数据的每天的偏度(a)峰度(b)。

图5是所有时间内的s-k平面图；

图中灰色的点代表所有点，橙色的点是2007年1月1日到2009年7月1日的(s，k)，放大图中的灰色阴影圆是地壳平静背景；

图6是2007年1月1日到2009年7月1日的s-k偏差曲线。

具体实施方式

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

实施例

针对汶川地震，以四川地区地震前兆监测台站中姑咱台的钻孔应变数据为例。姑咱台站与汶川地震震中位置如图2所示。该数据由yry四分量钻孔应变仪测得，采样一分钟一次，研究的时间段为2007年1月1日至2009年7月1日。参见图1所示：

步骤一、录入钻孔应变数据,并进行数据有效性的验证，是，进行下一步，录入钻孔应变数据是选取一个台站的四分量钻孔应变数据，制作成按照分钟值采样的时间序列。数据有效性验证是按照不同分量，具体为：

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

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

计算出的k值，选取自洽系数k在0.9以上的数据为有效数据，由k≥0.9可知姑咱台站2007年1月1日至2013年12月31日钻孔应变分钟值数据有效。

步骤二、对钻孔应变数据进行面应变换算与调和分析处理：对钻孔应变数据进行应变换算是根据公式(2)将四分量钻孔应变观测数据换算成面应变sa。

sa(t)＝(s1(t)+s2(t)+s3(t)+s4(t))/2，(2)

然后，应用调和分析是用来去除数据中固体潮等低频信息的周期项。其调和分析函数s(t)表达式为：

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

步骤三、对处理后数据取差分。其差分表达式为：

dsa(t)＝(sa(t+1)-s(t+1))-(sa(t)-s(t))，(4)

然后，按天分段制成n天的样本数据ds(n),ds(n)表示表示差分后的面应变，其中2007年的1月1日到2009年7月1日的样本数据如图3所示。

步骤四、计算每天样本数据的偏度、峰度。这两个参数是利用统计学上高阶累积量的算法计算的，现有随机变量x，其均值是μ，方差是σ，表达式如下：

这里的μ，σ分别是样本数据ds(n)的平均值和方差，sk为偏度，ku为峰度。针对面应变差分数据，如果当天的样本数据满足对称分布的的，当天偏度sk为0，若偏度>0，则认为分布为右偏，即地壳有受拉张的趋势；若偏度<0，则可认为分布为左偏，即地壳有受挤压的趋势，同时偏度的绝对值越大，说明分布的偏移程度越严重。当峰度ku>0，从形态上看，它相比于正态分布要更陡峭或尾部更厚，这个分布往往比正态分布的尾部具有更大的“质量”，即含有更多的极端值；而峰度ku<0,从形态山看，则它相比于正态分布更平缓或尾部更薄。其中姑咱台站2007年1月1日至2009年7月1日的样本数据的偏度峰度分别如图4(a)和图4(b)所示。

步骤五、将偏度、峰度画在s-k平面上，是将每天样本数据的偏度作为横坐标，峰度作为纵坐标，横纵坐标确定了一个点，将该点代表该天样本数据的分布状态。步骤五就是将所有的(s,k)画在一个平面。发现总体上，偏度和峰度呈现抛物线关系。将峰度看作是偏度的函数，再将峰度在偏度为零处进行泰勒展开，可以解释钻孔形变的面应变差分数据在s-k平面有一种成抛物线的非线性关系。事实上，不仅是面应变差分曲线，面应变原始曲线，观测曲线及观测曲线的差分形式均有此关系。

步骤六、输出s-k平面图，所有时间段内的s-k平面图如图5所示，灰色点是全部时间内的(s,k)，深黑点是汶川地震研究时间段的(s,k)，可以看到汶川前的数据大多分布在抛物线的左侧，意味着汶川地震前大多数天都受挤压力的影响。

步骤七、因为要判断异常，必然要有一个判断异常的标准。理论上，高斯分布的偏度峰度都为零。没有异常的平静天里，钻孔面应变的样本数据是满足准高斯分布的。所以，在s-k平面上，定义了一个地壳平静背景。这里将背景定义为一个圆，在所有的点中去掉一部分的极端点，剩余相对平稳的点求剩余点坐标参数的平均值作为圆点，剩余点坐标参数的标准差的平均值作为圆的半径。这样背景圆为：

图5上部分的小图是背景处的放大，灰色阴影是定义的地壳平静背景圆。

步骤八、计算每天偏离背景的程度。步骤七定义了一个背景圆，认为这个圆就是阈值，在圆内的点就是正常值，在圆外的点就是异常值。为了判断异常的状态，这里的状态是指既要判断异常点与圆的距离，又要判断异常点是在抛物线的左侧还是右侧。所以，定义了一个偏移d。d的表达式为：

根号下的部分代表异常的大小，前面乘以偏度的符号代表异常的方向。

步骤九、输出姑咱台站2007年1月1日至2009年7月1日的样本数据的s-k偏差随时间变化的偏移曲线。可以看出汶川地震前偏差一直为很大的负值，说明其观测地区地壳一直长期主要受到挤压力的作用。

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