CN102879813B - 一种微地震信号到时自动拾取的方法及装置 - Google Patents

Abstract

本发明公开了一种微地震信号到时自动拾取的方法及装置，包括：获取地震波的N个样点，通过得到的第i个样点前M个地震波形数据计算第i个样点处的峭度；设置得到的第i个样点处的峭度作为特征函数，通过特征函数的长短时均值比和短时均值判断微地震信号的初至位置。本发明解决了在致密砂岩油气藏和页岩气藏的开采中，如何对微地震进行到时自动拾取地震信号的问题。

Description

一种微地震信号到时自动拾取的方法及装置

技术领域

本发明涉及微地震信号采集领域，尤其涉及一种微地震信号到时自动拾取的方法及装置。

背景技术

近年来美国大规模的页岩气开采引发了各国对页岩气资源的关注。我国致密砂岩油气藏和页岩气藏资源丰富，在松辽、鄂尔多斯、四川及南方诸多盆地均有分布，这些致密砂岩油气藏和页岩气藏的天然气资源将成为我国重要的能源支柱。对致密砂岩油气藏和页岩气藏主要采用水平井分段水力压裂技术方式开采,但水力压裂使地层形成新的裂缝，产生微震。监测微地震的分布和特征可为水力压裂设计和压裂效果的监测提供依据，提高开采效益。

微地震记录的特点是频率高、信噪比低，因此微地震事件的自动识别和初至到时拾取对实现海量微地震数据的自动处理有重要意义。对于天然地震事件，已提出了多种自动识别方法。这些方法主要根据在地震记录中地震事件到达前后质点振动性质差别构建特征函数加以判断。例如，根据在时间域能量和能量变化构建特征函数提出的长短时均值比方法(STA/LTA法)；根据在时间-频率域能量变化提出的频谱比方法；根据地震信号到达前后地震波形数据统计性质的差别提出的AIC方法，还有基于地震波形偏斜度和峭度的PAISK方法；以及分形分维方法等。与天然地震的记录相比，微震的震级更小，通常在1级以下，而且信噪比更低。因此在致密砂岩油气藏和页岩气藏的开采中，如何对微地震进行到时自动拾取地震信号尤为重要。

发明内容

本发明所要解决的技术问题是提供一种微地震信号到时自动拾取的方法及装置，解决了在致密砂岩油气藏和页岩气藏的开采中，如何对微地震进行到时自动拾取地震信号的问题。

为了解决上述问题，本发明提供了一种微地震信号到时自动拾取的方法，包括：

获取地震波的N个样点，通过得到的第i个样点前M个地震波形数据计算第i个样点处的峭度K_i，其中N、M、i为自然数；

设置得到的第i个样点处的峭度K_i作为特征函数E_i，通过特征函数E_i的长短时均值比判断微地震信号的初至位置，其中，短时均值函数为：α_i=(1-C₃)·α_i-1+C₃·E_i；长时均值函数为：β_i=(1-C₄)·β_i-1+C₄·E_i；若地震信号的位置满足下列条件，则判定该位置为微地震信号的初至位置：

1)α_i≥r₁·β_i且α_i≥r₂；

2)满足条件1的位置记为n，由位置n向时间减小方向进行搜索K_i的数值，直到K_i≤K_i-1；

其中，C₃为短时均值的权系数、C₄为长时均值的权系数；r₁为长短时均值比阀值，r₂为峭度短时均值的阀值。

本发明还提供了一种微地震信号到时自动拾取的装置，包括：

获取地震波形峭度模块，用于获取地震波的N个样点，通过得到的第i个样点前M个地震波形数据计算第i个样点处的峭度K_i，其中N、M、i为自然数；及

判断初至位置模块，用于设置得到的第i个样点处的峭度K_i作为特征函数E_i，通过特征函数E_i的长短时均值比判断微地震信号的初至位置，其中，短时均值函数为：α_i=(1-C₃)·α_i-1+C₃·E_i；长时均值函数为：β_i=(1-C₄)·β_i-1+C₄·E_i；若地震信号的位置满足下列条件，则判定该位置为微地震信号的初至位置：

1)α_i≥r₁·β_i且α_i≥r₂；

2)满足条件1的位置记为n，由位置n向时间减小方向进行搜索K_i的数值，直到K_i≤K_i-1；

其中，C₃为短时均值的权系数、C₄为长时均值的权系数；r₁为长短时均值比阀值，r₂为峭度短时均值的阀值；及

存储模块，用于存储数据。

与现有技术相比，应用本发明，减少了系统的重复计算，将计算量减少为O(N)，意味着是与N成正比，极大地提高了计算效率，解决了当前传统的基于地震波形峭度和偏斜度的检测方法存在运算量过大导致效率低的问题。与原PAISK方法相比，本发明的到时拾取准确率提高约10个百分点，与传统的STA／LTA方法相比，准确率提高约5个百分点。

附图说明

图1是一地震记录的波形和峭度、以及峭度的长短时均值曲线的示意图；

图2是图1的初至附近的局部放大图；

图3是本发明的微地震信号到时自动拾取的方法的流程图；

图4是本发明的微地震信号到时自动拾取的装置的结构示意图。

具体实施方式

本发明的方法在具体实现中要通过系统中各设备之间信息交互来进行信息和/或数据的收集，并通过其内的控制器（可以是CPU等进行控制处理信息和/或数据，本发明对此不作任何限定），其间还可以通过各种存储器（可以是内存、硬盘或其他存储设备）进行信息和/或数据的储存和传送，本发明对此不作任何限定。

下面结合附图和具体实施方式对本发明作进一步说明。

传统的PAISK方法，即基于地震波形峭度和偏斜度的拾取初至方法，其中，对于一个有限长度的离散实数序列{x_i}，其偏斜度S(skewness)和峭度K(kurtosis)定义为：

$K=\frac{{\mathrm{\Σ}}_{i=1}^M{(x_i-\stackrel{\‾}{x})}^4}{(M-1)\·{\mathrm{\σ}}_x^4}-3---\left(1\right)$

$S=\frac{{\mathrm{\Σ}}_{i=1}^M{(x_i-\stackrel{\‾}{x})}^3}{(M-1)\·{\mathrm{\σ}}_x^3}---\left(2\right)$

其中为{x_i}序列的平均值，σ_x为其标准差，峭度的时窗长度为M。

从上可以看出，对于一个长度为N的地震波形序列{x_i}，峭度的时窗长度为M，按照公式(1)沿地震波形滑动求取峭度，要调用（N-M）次公式，计算量为O(M·(N-M))，对于海量的多道连续观测记录，其计算效率较低。本发明根据PAISK方法，提出移动时窗峭度的快速算法，将上述公式展开，将公式变换为求取离散实数序列的算术和、平方和、立方和、4次方和的线性组合形式，减少重复计算，计算量减少为O(N)，极大地提高了计算效率。

将（1）式展开并变换得到：

$K=\frac{1}{(M-1)\·{\mathrm{\σ}}_x^4}({\mathrm{\Σ}}_{i=1}^M{x}_i^4-4\stackrel{\‾}{x}{\mathrm{\Σ}}_{i=1}^M{x}_i^3+6{\stackrel{\‾}{x}}^2{\mathrm{\Σ}}_{i=1}^M{x}_i^2-3M\·{\stackrel{\‾}{x}}^4)-3---\left(3\right)$

将时窗内离散序列的算术和、平方和、立方和、4次方和分别记为sm_j,j=1,2,3,4，得到：

$K=\frac{1}{(M-1)\·{\mathrm{\σ}}_x}({\mathrm{sm}}_4-4\stackrel{\‾}{x}\·{\mathrm{sm}}_3+6{\stackrel{\‾}{x}}^2\·{\mathrm{sm}}_2-3M{\stackrel{\‾}{x}}^4)-3---\left(4\right)$

其中为：

$\stackrel{\‾}{x}=\frac{1}{M}{\mathrm{sm}}_1---\left(5\right)$

类似地，对于标准差σ_x，有下列关系：

${\mathrm{\σ}}_x={\left(\frac{1}{M}{\mathrm{\Σ}}_{i=1}^M{(x_i-\stackrel{\‾}{x})}^2\right)}^{\frac{1}{2}}={\left(\frac{1}{M}({\mathrm{sm}}_2-M\·{\stackrel{\‾}{x}}^2)\right)}^{\frac{1}{2}}---\left(6\right)$

采用（4）-（6）式计算的优点是，当时窗移动时计算新时窗的sm_j不必重复计算时窗内所有点的值，只须在上一个时窗各值的基础上，加上新进入时窗的值，并减去离开时窗的值即可。经实际数据计算测试，在时窗长度取0.5s、1.0s、1.5s的情况下，采用快速算法的计算速度分别提高了22倍、32倍、43倍。对于偏斜度的计算，同样可采用类似的方法提高速度。

传统的STA／LTA方法涉及参数过多，比较复杂，其参数意义及取值如表1所示。

表1.STA/LTA方法参数意义及取值

对微地震信号进行到时判别，PAISK方法采用峭度最大值点左侧斜率最大的点作为初至位置，但实际数据的结果显示，该点往往滞后于真实的初至位置，而峭度曲线开始起跳的位置更接近初至位置。在某些情况下，初至附近的峭度极值点也不是最大值点。针对这些情况，本发明将STA/LTA方法与PAISK法结合，以地震波形的峭度K作为特征函数，以特征函数的长短时均值比做为判断初至位置的依据。具体方法如下：

定义特征函数E_i为：

R_i=K_i，（7）

其中K_i为第i个样点处的峭度，K_i的值由第i个样点前M个地震波形数据计算。与STA／LTA的方法类似，长短时均值函数定义为：

α_i=(1-C₃)·α_i-1+C₃·E_i,（8）

β_i=(1-C₄)·β_i-1+C₄·E_i,（9）

满足下列条件的位置判定为初至：

1)α_i≥r₁·β_i并且α_i≥r₂

2)满足条件1的位置记为n，由n开始向前(时间减小方向)搜索K_i的值，直到满足K_i≤K_i-1。

图1为一个地震记录的波形和峭度、以及峭度的长短时均值曲线的示意图，其中，从上到下分别为：垂直分量地震记录波形；地震波形的峭度曲线；峭度曲线的长时均值；峭度曲线的短时均值。

图2为图1的初至附近的局部放大图，图2中初至附近的波形(实线)峭度曲线长时均值(虚线)、短时均值(点线)。波形曲线按照峭度曲线做了归一化处理。

上述公式中的C₃、C₄和r₁、r₂为常数。将上述方法应用于实际数据，采用DE全局搜索法搜索最佳参数，保持C₃=0.6,C₄=0.03不变，以人工拾取的初至作为参照，以时差在0.3s以内的记录所占百分比作为目标函数，自动搜索最佳的r₁和r₂的数值。所采用的实际地震数据为我国西部某地观测到的震中距范围在0.3至20公里、震级范围-0.5至1级的约13359个微震记录。当r₁=2.71,r₂=1.43时，自动识别的效果最佳。得到的结果如表2中所示。表2中还列出了传统的STA／LTA方法、AIC方法、以及PAISK方法的结果，对于STA／LTA方法同样采用了DE全局搜索方法搜索最佳拾取参数。从表2可以看出，本发明方法识别准确率较PAISK方法提高了约10％，比传统的STA／LTA方法提高了约5％。

表2.各种方法初至拾取情况一览表（其中计算耗时未计入硬盘I/O传送时间）

本方法主要有两个控制参数r₁和r₂，分别是长短时均值比阀值和峭度短时均值的阀值，在P波初至附近的极值点处，长短时的均值比和极值都随地震信号的信噪比增加而增加，峭度的极值尤其明显。而根据(1)式和(2)式，偏斜度和峭度均为无量纲量，因而对于不同类型的地震记录(速度、加速度、不同的仪器增益等)，该方法的拾取参数具有一定的普遍意义。

如图3所示，本发明的微地震信号到时自动拾取的方法，包括以下步骤：

步骤310、获取地震波的N个样点，通过得到的第i个样点前M个地震波形数据计算第i个样点处的峭度K_i，其中N、M、i为自然数；

所述获取地震波的N个样点，通过得到的第i个样点前M个地震波形数据计算第i个样点处的峭度K_i的步骤，包括：

获取地震波的N个样点的离散序列{x_i}，当前时窗长度为M；

计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到第i个样点处的峭度K_i。

所述计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到第i个样点处的峭度K_i的步骤，包括：

计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和，其中，将当前时窗长度为M的时窗内离散序列的算术和、平方和、立方和、4次方和分别记为sm_j,j=1,2,3,4，公式如下所示：

${\mathrm{sm}}_1={\mathrm{\Σ}}_{i=1}^M{x}_i,$

${\mathrm{sm}}_2={\mathrm{\Σ}}_{i=1}^M{x}_i^2,$

${\mathrm{sm}}_{\text{3}}={\mathrm{\Σ}}_{i=1}^M{x}_i^3,$

${\mathrm{sm}}_4={\mathrm{\Σ}}_{i=1}^M{x}_i^4;$

然后通过公式

$K=\frac{1}{(M-1)\·{\mathrm{\σ}}_x}({\mathrm{sm}}_4-4\stackrel{\‾}{x}\·{\mathrm{sm}}_3+6{\stackrel{\‾}{x}}^2\·{\mathrm{sm}}_2-3M{\stackrel{\‾}{x}}^4)-3$ 得到当前时窗长度为M的峭度K，

其中为：

$\stackrel{\‾}{x}=\frac{1}{M}{\mathrm{sm}}_1,$

标准差 ${\mathrm{\σ}}_x=\mathrm{sqrt}\left(\frac{1}{M}{\mathrm{\Σ}}_{i=1}^M{(x_i-\stackrel{\‾}{x})}^2\right)=\mathrm{sqrt}\left(\frac{1}{M}({\mathrm{sm}}_2-M\·{\stackrel{\‾}{x}}^2)\right);$

将时窗向前移动一个样点，计算新时窗内的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值f，并加上当前时窗最后一个样点的值r，公式如下所示：

${\mathrm{sm}}_1^n={\mathrm{sm}}_1+(r-f),$

${\mathrm{sm}}_2^n={\mathrm{sm}}_2+(r^2-f^2),$

${\mathrm{sm}}_3^n={\mathrm{sm}}_3+(r^3-f^3),$

${\mathrm{sm}}_4^n={\mathrm{sm}}_4+(r^4-f^4),$

得到新时窗的后，再计算新时窗的峭度的数值，依次类推，

得到第i个样点处的峭度K_i。

步骤320、设置得到的第i个样点处的峭度K_i作为特征函数E_i，通过特征函数E_i的长短时均值比和短时均值判断微地震信号的初至位置，其中，短时均值函数为：α_i=(1-C₃)·α_i-1+C₃·E_i；长时均值函数为：β_i=(1-C₄)·β_i-1+C₄·E_i；若地震信号的位置满足下列条件，则判定该位置为微地震信号的初至位置：

1)α_i≥r₁·β_i且α_i≥r₂；

2)满足条件1的位置记为n，由位置n向时间减小方向进行搜索K_i的数值，直到K_i≤K_i-1；

其中，C₃、C₄、r₁和r₂为常数，C₃为短时均值的权系数、C₄为长时均值的权系数，一般取C₃=0.6，C₄=0.03；r₁为长短时均值比阀值，r₂为峭度短时均值的阀值，应根据实际数据情况选择，对于本文中的实例数据，r₁=2.71,r₂=1.43

如图4所示，本发明的一种微地震信号到时自动拾取的装置，包括：获取地震波形峭度模块和判断初至位置模块，其中，

获取地震波形峭度模块，用于获取地震波的N个样点，通过得到的第i个样点前M个地震波形数据计算第i个样点处的峭度K_i，其中N、M、i为自然数；

判断初至位置模块，用于设置得到的第i个样点处的峭度K_i作为特征函数E_i，通过特征函数E_i的长短时均值比和短时均值判断微地震信号的初至位置，其中，短时均值函数为：α_i=(1-C₃)·α_i-1+C₃·E_i；长时均值函数为：β_i=(1-C₄)·β_i-1+C₄·E_i；若地震信号的位置满足下列条件，则判定该位置为微地震信号的初至位置：

1)α_i≥r₁·β_i且α_i≥r₂；

2)满足条件1的位置记为n，由位置n向时间减小方向进行搜索K_i的数值，直到K_i≤K_i-1；

其中，C₃、C₄、r₁和r₂为常数，C₃=0.6,C₄=0.03，对于本文的实例数据，r₁=2.71,r₂=1.43；及

存储模块，用于存储数据。

所述获取地震波形峭度模块获取地震波的N个样点，通过得到的第i个样点前M个地震波形数据计算第i个样点处的峭度K_i，是指：

所述获取地震波形峭度模块获取地震波的N个样点的离散序列{x_i}，当前时窗长度为M；计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到第i个样点处的峭度K_i。

所述获取地震波形峭度模块计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到第i个样点处的峭度K_i，是指：

所述获取地震波形峭度模块计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和，其中，将当前时窗长度为M的时窗内离散序列的算术和、平方和、立方和、4次方和分别记为sm_j,j=1,2,3,4，公式如下所示：

${\mathrm{sm}}_1={\mathrm{\Σ}}_{i=1}^M{x}_i,$

${\mathrm{sm}}_2={\mathrm{\Σ}}_{i=1}^M{x}_i^2,$

${\mathrm{sm}}_3={\mathrm{\Σ}}_{i=1}^M{x}_i^3,$

${\mathrm{sm}}_4={\mathrm{\Σ}}_{i=1}^M{x}_i^4,$

然后通过公式

$K=\frac{1}{(M-1)\·{\mathrm{\σ}}_x}({\mathrm{sm}}_4-4\stackrel{\‾}{x}\·{\mathrm{sm}}_3+6{\stackrel{\‾}{x}}^2\·{\mathrm{sm}}_2-3M{\stackrel{\‾}{x}}^4)-3$ 得到当前时窗长度为M的峭度K，

其中为：

$\stackrel{\‾}{x}=\frac{1}{M}{\mathrm{sm}}_1,$

标准差 ${\mathrm{\σ}}_x=\mathrm{sqrt}\left(\frac{1}{M}{\mathrm{\Σ}}_{i=1}^M{(x_i-\stackrel{\‾}{x})}^2\right)=\mathrm{sqrt}\left(\frac{1}{M}({\mathrm{sm}}_2-M\·{\stackrel{\‾}{x}}^2)\right);$

将时窗向前移动一个样点，计算新时窗内的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值f，并加上当前时窗最后一个样点的值r，公式如下所示：

${\mathrm{sm}}_1^n={\mathrm{sm}}_1+(r-f),$

${\mathrm{sm}}_2^n={\mathrm{sm}}_2+(r^2-f^2),$

${\mathrm{sm}}_3^n={\mathrm{sm}}_3+(r^3-f^3),$

${\mathrm{sm}}_4^n={\mathrm{sm}}_4+(r^4-f^4),$

得到新时窗的后，再计算新时窗的峭度的数值，依次类推，得到第i个样点处的峭度K_i。

下面通过实例对本发明作进一步说明。

通过传统的STA／LTA方法、AIC方法、以及PAISK方法与本发明方法拾取结果的对比，得到如下结论：

STA／LTA法历史悠久，算法稳定强健，目前仍被广泛使用。该方法可沿地震记录波形滑动进行，因而可用于地震信号的检测，但该方法涉及参数较多，不易掌控。

AIC方法在时窗正好包含了地震信号前后一段波形的情况下，能够得到非常好的结果，适用于已有震源先验信息的情形，例如流动地震观测中记录到的远震信号和区域地震信号，以及在已经用其他方法检测出地震信号所在时段后，再应用AIC方法判断到时。

PAISK方法原理简单实现容易，算法同样稳定强健，由于该方法采用了高阶统计量判断信号初至，具有一定的抗干扰能力。

本发明提出的方法，可方便地沿连续地震记录波形滑动进行，因而可用于地震信号的检测，对于本文中的实例数据，拾取误差小于0.3秒的记录占比达到83.8%，是上述方法中最高的；该方法涉及的参数较少，同时由于峭度的无量纲特性，因而所涉及参数的选取具有一定的普遍性。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (4)

1.一种微地震信号到时自动拾取的方法，其特征在于，包括：

获取地震波的N个样点，通过得到的第i个样点前M个地震波形数据计算第i个样点处的峭度K_i，其中N、M、i为自然数；其中具体包括步骤：获取地震波的N个样点的离散序列{x_i}，当前时窗长度为M；计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到第i个样点处的峭度K_i；

设置得到的第i个样点处的峭度K_i作为特征函数E_i，通过特征函数E_i的长短时均值比判断微地震信号的初至位置，其中，短时均值函数为：α_i＝(1-C₃)·α_i-1+C₃·E_i；长时均值函数为：β_i＝(1-C₄)·β_i-1+C₄·E_i；若地震信号的位置满足下列条件，则判定该位置为微地震信号的初至位置：

1)α_i≥r₁·β_i且α_i≥r₂；

2)满足条件1的位置记为n，由位置n向时间减小方向进行搜索K_i的数值，直到K_i≤K_i-1；

其中，C₃为短时均值的权系数、C₄为长时均值的权系数；r₁为长短时均值比阀值，r₂为峭度短时均值的阀值。

2.如权利要求1所述的方法，其特征在于，

所述计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到第i个样点处的峭度K_i的步骤，包括：

计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和，其中，将当前时窗长度为M的时窗内离散序列的算术和、平方和、立方和、4次方和分别记为sm_j,j＝1,2,3,4，公式如下所示：

${\mathrm{sm}}_1={\mathrm{\Σ}}_{i=1}^M{x}_i,$

${\mathrm{sm}}_2={\mathrm{\Σ}}_{i=1}^m{x}_i^2,$

${\mathrm{sm}}_3={\mathrm{\Σ}}_{i=1}^M{x}_i^3,$

${\mathrm{sm}}_4={\mathrm{\Σ}}_{i=1}^M{x}_i^4;$

然后通过公式

$K=\frac{1}{(M-1)\·{\mathrm{\σ}}_x}({\mathrm{sm}}_4-4\stackrel{\‾}{x}\·{\mathrm{sm}}_3+6{\stackrel{\‾}{x}}^2\·{\mathrm{sm}}_2-3M{\stackrel{\‾}{x}}^4)-3$ 得到当前时窗长度为M的峭度K，

其中为：

$\stackrel{\‾}{x}=\frac{1}{M}{\mathrm{sm}}_1,$

标准差 ${\mathrm{\σ}}_x=\mathrm{sqrt}\left(\frac{1}{M}{\mathrm{\Σ}}_{i=1}^M{(x_i-\stackrel{\‾}{x})}^2\right)=\mathrm{sqrt}\left(\frac{1}{M}({\mathrm{sm}}_2-M\·{\stackrel{\‾}{x}}^2)\right);$

将时窗向前移动一个样点，计算新时窗内的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值f，并加上当前时窗最后一个样点的值r，公式如下所示：

${\mathrm{sm}}_1^n={\mathrm{sm}}_1+(r-f),$

${\mathrm{sm}}_2^n={\mathrm{sm}}_2+(r^2-f^2),$

${\mathrm{sm}}_3^n={\mathrm{sm}}_3+(r^3-f^3),$

${\mathrm{sm}}_4^n={\mathrm{sm}}_4+(r^4-f^4),$

得到新时窗的后，再计算新时窗的峭度的数值，依次类推，得到第i个样点处的峭度K_i。

3.一种微地震信号到时自动拾取的装置，其特征在于，包括：

获取地震波形峭度模块，用于获取地震波的N个样点，通过得到的第i个样点前M个地震波形数据计算第i个样点处的峭度K_i，其中N、M、i为自然数，其中所述获取地震波形峭度模块获取地震波的N个样点，通过得到的第i个样点前M个地震波形数据计算第i个样点处的峭度K_i，是指：所述获取地震波形峭度模块获取地震波的N个样点的离散序列{x_i}，当前时窗长度为M；计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到第i个样点处的峭度K_i；及

判断初至位置模块，用于设置得到的第i个样点处的峭度K_i作为特征函数E_i，通过特征函数E_i的长短时均值比判断微地震信号的初至位置，其中，短时均值函数为：α_i＝(1-C₃)·α_i-1+C₃·E_i；长时均值函数为：β_i＝(1-C₄)·β_i-1+C₄·E_i；若地震信号的位置满足下列条件，则判定该位置为微地震信号的初至位置：

1)α_i≥r₁·β_i且α_i≥r₂；

2)满足条件1的位置记为n，由位置n向时间减小方向进行搜索K_i的数值，直到K_i≤K_i-1；

其中，C₃为短时均值的权系数、C₄为长时均值的权系数；r₁为长短时均值比阀值，r₂为峭度短时均值的阀值；及

存储模块，用于存储数据。

4.如权利要求3所述的装置，其特征在于，包括：

所述获取地震波形峭度模块计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到第i个样点处的峭度K_i，是指：

所述获取地震波形峭度模块计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和，其中，将当前时窗长度为M的时窗内离散序列的算术和、平方和、立方和、4次方和分别记为sm_j,j＝1,2,3,4，公式如下所示：

${\mathrm{sm}}_1={\mathrm{\Σ}}_{i=1}^M{x}_i,$

${\mathrm{sm}}_2={\mathrm{\Σ}}_{i=1}^m{x}_i^2,$

${\mathrm{sm}}_3={\mathrm{\Σ}}_{i=1}^M{x}_i^3,$

${\mathrm{sm}}_4={\mathrm{\Σ}}_{i=1}^M{x}_i^4;$

然后通过公式

$K=\frac{1}{(M-1)\·{\mathrm{\σ}}_x}({\mathrm{sm}}_4-4\stackrel{\‾}{x}\·{\mathrm{sm}}_3+6{\stackrel{\‾}{x}}^2\·{\mathrm{sm}}_2-3M{\stackrel{\‾}{x}}^4)-3$ 得到当前时窗长度为M的峭度K，

其中为：

$\stackrel{\‾}{x}=\frac{1}{M}{\mathrm{sm}}_1,$

标准差 ${\mathrm{\σ}}_x=\mathrm{sqrt}\left(\frac{1}{M}{\mathrm{\Σ}}_{i=1}^M{(x_i-\stackrel{\‾}{x})}^2\right)=\mathrm{sqrt}\left(\frac{1}{M}({\mathrm{sm}}_2-M\·{\stackrel{\‾}{x}}^2)\right);$

将时窗向前移动一个样点，计算新时窗内的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值f，并加上当前时窗最后一个样点的值r，公式如下所示：

${\mathrm{sm}}_1^n={\mathrm{sm}}_1+(r-f),$

${\mathrm{sm}}_2^n={\mathrm{sm}}_2+(r^2-f^2),$

${\mathrm{sm}}_3^n={\mathrm{sm}}_3+(r^3-f^3),$

${\mathrm{sm}}_4^n={\mathrm{sm}}_4+(r^4-f^4),$

得到新时窗的后，再计算新时窗的峭度的数值，依次类推，得到第i个样点处的峭度K_i。