CN102707312B - 一种检测地震信号的控制方法及装置 - Google Patents

Abstract

一种检测地震信号的控制方法及系统，包括：获取地震波的N个样点的离散序列；计算当前时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到所有时窗N个样点的峭度和偏斜度的数值得到峭度曲线和偏斜度曲线，判断地震波到时的位置。本发明解决了当前传统的基于地震波形峭度和偏斜度的检测方法存在运算量过大导致效率低的问题。

Description

一种检测地震信号的控制方法及装置

技术领域

本发明涉及地震信号采集领域，尤其涉及一种检测地震信号的控制方法及装置。

背景技术

地震勘探主要分采集、处理和解释三大部分，地震数据采集时需要拾取地震信号，地震波初至的拾取广泛应用在表层调查方法和地质构造的探测中,初至的拾取对矿物资源和石油资源的开发，都起着重要的作用。

在地震观测中，地震事件的检测识别及初至到时拾取是一项重要而又繁琐的工作，当前出现了一种基于地震波形峭度和偏斜度的检测方法。对于一个有限长度的离散实数序列{x_i}，其偏斜度S(skewness)和峭度K(kurtosis,又译峰度)定义为：

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

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

其中为{x_i}序列的平均值，σ_x为其标准差。

通过地震波形峭度和偏斜度检测地震波到时的方法具有计算稳定强健、实现简单等优点。

但是对于一个长度为N的地震波形序列{x_i}，求偏斜度和峭度的时窗长度为M，按照公式(1.1)和（1.2），沿地震波形滑动求取偏斜度和峭度，要调用（N-M）次公式，计算量为O(M·(N-M))，意味着是跟M*(N-M)成正比，这样对于海量的多道连续观测记录，系统的计算效率较低。

因此当前需要一种检测地震信号控制的技术方案在初至自动拾取地震信号时可以提高系统的效率。

发明内容

本发明所要解决的技术问题是提供一种检测地震信号的控制方法及装置，解决了当前传统的基于地震波形峭度和偏斜度的检测方法存在运算量过大导致效率低的问题。

为了解决上述问题，本发明提供了一种检测地震信号的控制方法，包括：

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

计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到所有时窗N个样点的峭度和偏斜度的数值；并通过得到的所有时窗N个样点的峭度和偏斜度的数值得到峭度曲线和偏斜度曲线；

通过得到的地震波的峭度曲线和偏斜度曲线，判断地震波到时的位置。

进一步地，上述方法还可包括：所述得到所有时窗N个样点的峭度K的数值的步骤，包括：

计算当前时窗长度为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)$ 得到当前时窗长度为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),$

得到新时窗的后，再计算新时窗的峭度K的数值，依次类推，直到得到所有时窗的峭度K的数值。

进一步地，上述方法还可包括：所述得到所有时窗N个样点的偏斜度S的数值的步骤，包括：

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

${\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;$

然后通过公式 $S=\frac{1}{(M-1)\·{\mathrm{\σ}}_x^3}({\mathrm{sm}}_3-3{\mathrm{sm}}_2\stackrel{\‾}{x}+2M{\stackrel{\‾}{x}}^3)$ 得到当前时窗长度为M的时窗的偏斜度S，

其中为：

$\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，公式为：

sm₁=sm₁+（r-f)，

sm₂=sm₂+（r²-f²），

sm₃=sm₃+（r³-f³），

得到新时窗的后，计算新时窗的偏斜度S的数值，依次类推，直到得到所有时窗的偏斜度S的数值。

进一步地，上述方法还可包括：所述通过得到的地震波的峭度曲线和偏斜度曲线，判断地震波到时的位置的步骤，包括：

通过地震波形峭度和偏斜度检测地震波到时来判断峭度曲线和偏斜度曲线的极大值，根据极大值点前曲线斜率最大的位置判断地震波到时的位置。

本发明还提供了一种检测地震信号的控制装置，包括：

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

计算模块，用于计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到所有时窗N个样点的峭度和偏斜度的数值；

存储模块，用于存储数据；及

判断模块，用于通过得到的所有时窗N个样点的峭度和偏斜度的数值得到峭度曲线和偏斜度曲线，并通过得到的地震波的峭度曲线和偏斜度曲线，判断地震波到时的位置。

进一步地，上述装置还可包括：所述计算模块得到所有时窗N个样点的峭度K的数值，是指：

所述计算模块计算当前时窗长度为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)$ 得到当前时窗长度为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),$

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

直到得到所有时窗的峭度K的数值。

进一步地，上述装置还可包括：所述计算模块得到所有时窗N个样点的偏斜度S的数值，是指：

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

${\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;$

然后通过公式 $S=\frac{1}{(M-1)\·{\mathrm{\σ}}_x^3}({\mathrm{sm}}_3-3{\mathrm{sm}}_2\stackrel{\‾}{x}+2M{\stackrel{\‾}{x}}^3)$ 得到当前时窗长度为M的时窗的偏斜度S，

其中为：

$\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，公式为：

sm₁=sm₁+（r-f)，

sm₂=sm₂+（r²-f²），

sm₃=sm₃+（r³-f³），

得到新时窗的后，计算新时窗的偏斜度S的数值，依次类推，直到得到所有时窗的偏斜度S的数值。

进一步地，上述装置还可包括：所述判断模块通过得到的地震波的峭度曲线和偏斜度曲线，判断地震波到时的位置，是指所述判断模块通过地震波形峭度和偏斜度检测地震波到时来判断峭度曲线和偏斜度曲线的极大值，根据极大值点前曲线斜率最大的位置判断地震波到时的位置。

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

附图说明

图1是是沿滑动时窗求取j次方和的过程示意图；

图2是本发明的检测地震信号的控制方法的流程图；

图3是本发明的检测地震信号的控制装置的结构示意图。

具体实施方式

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

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

通过地震波形峭度和偏斜度检测地震波到时来判断峭度曲线和偏斜度曲线的极大值，以极大值点前曲线斜率最大的位置作为地震波到时的位置，在本发明的地震信号采集的控制方法中，首先求峭度K，将1.1式展开并变换：

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

$=\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-4{\stackrel{\‾}{x}}^3{\mathrm{\Σ}}_{i=1}^M{x}_i+M\·{\stackrel{\‾}{x}}^4)$

$=\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)---\left(2.1\right)$

其中为{x_i}序列的平均值，σ_x为其标准差，一个长度为N的地震波形序列{x_i}，偏斜度和峭度的时窗长度为M；

将时窗内离散序列的算术和、平方和、立方和、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)---\left(2.2\right)$

其中为：

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

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

${\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)---\left(2.4\right)$

采用（2.2）-（2.3）式计算的优点是，当时窗移动时，计算新时窗的sm_j不必重复计算时窗内所有点的值，只须加上新进入时窗的值，并减去离开时窗的值即可。

其次，对于偏斜度S的计算，将1.2式展开，得到：

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

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

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

将上式中的求和项记为sm_j，得到：

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

其中平均值和标准差的计算与2.3式和2.4式相同。

下面具体说明计算过程：

对于N个样点的离散序列{x_i},计算峭度的时窗长度为M，计算峭度曲线的过程如下：

步骤1、计算第1个时窗的峭度，首先计算该时窗内样点的算术和、平方和、立方和、4次方和：

${\mathrm{sm}}_1={\mathrm{\Σ}}_{i=1}^M{x}_i---\left(3.1\right)$

${\mathrm{sm}}_2={\mathrm{\Σ}}_{i=1}^M{x}_i^2---\left(3.2\right)$

${\mathrm{sm}}_3={\mathrm{\Σ}}_{i=1}^M{x}_i^3---\left(3.3\right)$

${\mathrm{sm}}_4={\mathrm{\Σ}}_{i=1}^M{x}_i^4---\left(3.4\right)$

然后利用公式（2.2-2.4）计算得到第1个时窗的峭度；

步骤2、时窗向前移动一个样点，新时窗内各值的计算，须在上一个时窗的值基础上，减去上一个时窗第1个样点的值f，并加上当前时窗最后一个样点的值r，如图1所示，具体计算公式为：

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

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

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

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

得到新时窗的后，再利用公式(2.2)-(2.4)计算新时窗的峭度值。

图1是沿滑动时窗求取j次方和的过程示意图。实线标注的是当前时窗，点线标出的是上一个时窗，r为新进入当前时窗的点，f为刚脱离当前时窗的点。

步骤3、重复2的步骤直到所有时窗计算完成。

对于N个样点的离散序列{x_i},计算偏斜度的时窗长度为M，计算偏斜度曲线的过程与计算峭度类似，具体过程如下：

步骤1、计算第1个时窗的峭度，首先计算该时窗内样点的算术和、平方和、立方和：

${\mathrm{sm}}_1={\mathrm{\Σ}}_{i=1}^M{x}_i---\left(3.9\right)$

${\mathrm{sm}}_2={\mathrm{\Σ}}_{i=1}^M{x}_i^2---\left(3.10\right)$

${\mathrm{sm}}_3={\mathrm{\Σ}}_{i=1}^M{x}_i^3---\left(3.11\right)$

然后利用公式（2.3）-（2.5）式计算得到第1个时窗的偏斜度。步骤2、时窗向前移动一个样点，新时窗内各值的计算，须在上一个时窗的值基础上，减去上一个时窗第1个样点的值f，并加上当前时窗最后一个样点的值r，具体计算公式为：

sm₁=sm₁+（r-f)    （3.12）

sm₂=sm₂+（r²-f²） （3.13）

sm₃=sm₃+（r³-f³ ）（3.14）

得到新时窗的后，再利用公式(2.3)-(2.5)计算新时窗的偏斜度值。

步骤3、重复2的步骤直到所有时窗计算完成。

如图2所示，本发明的检测地震信号的控制方法，包括以下步骤：

步骤210、获取地震波的N个样点的离散序列{x_i}，当前时窗长度为M，其中N、M和i都为自然数；

步骤220、计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值；依次类推，得到所有时窗N个样点的峭度和偏斜度的数值，并通过得到的所有时窗N个样点的峭度和偏斜度的数值得到峭度曲线和偏斜度曲线；

步骤230、通过地震波形峭度和偏斜度检测地震波到时来判断峭度曲线和偏斜度曲线的极大值，根据极大值点前曲线斜率最大的位置判断地震波到时的位置。

如图3所示，本发明的检测地震信号的控制装置，包括：获取地震波样点模块、计算模块、存储模块及判断模块，其中，

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

计算模块，用于计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到所有时窗N个样点的峭度和偏斜度的数值；

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

判断模块，用于通过得到的所有时窗N个样点的峭度和偏斜度的数值得到峭度曲线和偏斜度曲线，并通过得到的地震波的峭度曲线和偏斜度曲线，判断地震波到时的位置。

所述计算模块得到所有时窗N个样点的峭度K的数值，是指：

所述计算模块计算当前时窗长度为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)$ 得到当前时窗长度为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),$

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

直到得到所有时窗的峭度K的数值。

所述计算模块得到所有时窗N个样点的偏斜度S的数值，是指：

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

${\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;$

然后通过公式 $S=\frac{1}{(M-1)\·{\mathrm{\σ}}_x^3}({\mathrm{sm}}_3-3{\mathrm{sm}}_2\stackrel{\‾}{x}+2M{\stackrel{\‾}{x}}^3)$ 得到当

前时窗长度为M的时窗的偏斜度S，

其中为：

$\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，公式为：

sm₁=sm₁+（r-f)，

sm₂=sm₂+（r²-f²），

sm₃=sm₃+（r³-f³），

得到新时窗的后，计算新时窗的偏斜度S的数值，依次类推，直到得到所有时窗的偏斜度S的数值。

所述判断模块通过得到的地震波的峭度曲线和偏斜度曲线，判断地震波到时的位置，是指所述判断模块通过地震波形峭度和偏斜度检测地震波到时来判断峭度曲线和偏斜度曲线的极大值，根据极大值点前曲线斜率最大的位置判断地震波到时的位置。

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

本方法与传统方法的实际数据计算效率进行比较，其中所采用的数据为柴达木盆地西部地区进行的微地震观测数据，共有13359条垂直分量记录。每条记录地震波形长度15秒，滑动时窗长度分别采用了1秒/1.5秒/2秒三种情况。两种方法的计算时间如下表所示。计算时间为所有地震记录计算峭度或偏斜度的cpu耗时，未计入硬盘IO时间，测试在一台cpu I7-920、主板x58机器上进行。

通过以上可以看出，本方法与传统方法相比，计算速度大约提高了20倍左右，明显提升了系统的效率。

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

Claims (2)

1.一种检测地震信号的控制方法，其特征在于，包括：

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

计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到所有时窗N个样点的峭度和偏斜度的数值；并通过得到的所有时窗N个样点的峭度和偏斜度的数值得到峭度曲线和偏斜度曲线；

其中得到所有时窗N个样点的峭度K的数值的步骤，包括：

计算当前时窗长度为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)$ 得到当前时窗长

度为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),$

得到新时窗的后，再计算新时窗的峭度K的数值，依次类推，直到得到所有时窗的峭度K的数值；

其中得到所有时窗N个样点的偏斜度S的数值的步骤，包括：

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

${\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;$

然后通过公式 $S=\frac{1}{(M-1)\·{\mathrm{\σ}}_x^3}({\mathrm{sm}}_3-3{\mathrm{sm}}_2\stackrel{\‾}{x}+2M{\stackrel{\‾}{x}}^3)$ 得到当前时窗长度为M的时窗的偏斜度S，

其中为：

$\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，公式为：

sm₁＝sm₁+(r-f)   ，

sm₂＝sm₂+(r²-f²)   ，

sm₃＝sm₃+(r³-f³)   ，

得到新时窗的后，计算新时窗的偏斜度S的数值，依次类推，直到得到所有时窗的偏斜度S的数值；

通过得到的地震波的峭度曲线和偏斜度曲线，根据地震波形峭度和偏斜度检测地震波到时来判断峭度曲线和偏斜度曲线的极大值，根据极大值点前曲线斜率最大的位置判断地震波到时的位置。

2.一种检测地震信号的控制装置，其特征在于，包括：

获取地震波样点模块，用于获取地震波的N个样点的离散序列{x_i}发送给计算模块和存储模块，当前时窗长度为M；

计算模块，用于计算当前时窗长度为M的时窗内样点的算术和、平方和、立方和、4次方和；将时窗向前移动一个样点，计算新时窗内样点的各值，是在上一个时窗的值基础上，减去上一个时窗第1个样点的值，并加上当前时窗最后一个样点的值，依次类推，得到所有时窗N个样点的峭度和偏斜度的数值发送给判断模块和存储模块；

其中计算模块根据以下方式得到所有时窗N个样点的峭度K的数值：

计算模块计算当前时窗长度为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)$ 得到当前时窗长

度为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),$

得到新时窗的后，再计算新时窗的峭度K的数值，依修次改类替换推页，直到得到所有时窗的峭度K的数值；

其中计算模块根据以下方式得到所有时窗N个样点的偏斜度S的数值：

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

${\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;$

然后通过公式 $S=\frac{1}{(M-1)\·{\mathrm{\σ}}_x^3}({\mathrm{sm}}_3-3{\mathrm{sm}}_2\stackrel{\‾}{x}+2M{\stackrel{\‾}{x}}^3)$ 得到当前时窗长度为M的时窗的偏斜度S，

其中为：

$\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，公式为：

sm₁＝sm₁+(r-f)   ，

sm₂＝sm₂+(r²-f²)   ，

sm₃＝sm₃+(r³-f³)   ，

得到新时窗的后，计算新时窗的偏斜度S的数值，依次类推，直到得到所有时窗的偏斜度S的数值；

存储模块，用于存储数据；及

判断模块，用于通过得到的所有时窗N个样点的峭度和偏斜度的数值得到峭度曲线和偏斜度曲线，并通过得到的地震波的峭度曲线和偏斜度曲线，根据地震波形峭度和偏斜度检测地震波到时来判断峭度曲线和偏斜度曲线的极大值，根据极大值点前曲线斜率最大的位置判断地震波到时的位置。