CN101334482B - 一种预测地震波中的多次波和一次波信号的方法 - Google Patents

Abstract

本发明公开了一种预测地震波中的多次波和一次波信号的方法，属于地震观测信号处理技术领域。本发明所提供的预测地震波中的多次波信号的方法，包括以下步骤：1)用式(9)的目标函数优化滤波器系数a(i)；2)将步骤1)中得到的优化滤波器系数a(i)代入式(1)，得到多次波信号。用地震观测信号减去所述多次波的信号还可以得到地震波中的一次波信号。本发明克服了传统的预测方法需要假设预测多次波和一次波是正交的限制，更有效地消除了地震信号中的多次波。

(1)；maxJ(a(i))＝[E{G(ε₀(n))}-E{G(v(n))}]²，i＝1，2，…，m，(9)。

Description

一种预测地震波中的多次波和一次波信号的方法

技术领域

本发明属于地震信号处理领域，特别涉及一种预测地震波中的多次波和一次波信号的方法。

背景技术

海上多次波一直以来是海洋石油地震勘探中的一个关键问题。由于海面与空气的分界面是一个良好的反射界面，具有较大的反射系数，所以当震源产生的地震波经海底或地层分界面反射后传播到海面时，就会再次反射形成多次反射波。多次波的存在会干扰甚至破坏偏移剖面上真实界面的成像，因此，一般认为多次波是一种相干噪声，它常常干扰对一次波的有效识别。因此，可靠地识别多次波并削弱或消除多次波的影响是地震资料处理中的一项重要任务。

现有的消除多次波的方法大致可分为两类(Weglen，A B.，1999，Multipleattenuation：An overview of recent advance and the road ahead：The LeadingEdge，18，40-44)，一类是基于一次波和多次波差异性的滤波方法，简称为滤波方法(Lu，W.，Zhang，X.，Li，Y.，2003，Multiple removal based on detectionand estimation of localized coherent signal：Geophysics，68，745-750Niu，B.，Sun，C.，Zhang，Z.，2001，Polynomial Radon transform：Chinese J.Geophys.(in Chinese)，44，263-271)；另一类是基于波动方程的预测相减法(Lu，W.，2006，Adaptive Multiple Subtraction Using independent component analysis：Geophysics，71，S179-184Spitz，S.，1999，Pattern recognition，spatialpredictivity，and subtraction of multiple events：The Leading Edge，18，55-59Verschuur，D.，Berkhout，A.，Wapenaar，C.，1992，Adapt ivesurface-related multiple elimination：Geophysics，57，1166-1177)。

预测反褶积的基本理论是线性预测，作为理论在40年代就已建立起来。常规预测反褶积算法是一种基于多次波信号具有周期性的滤波方法(Robinson，E.，1957，Predictive decomposition of seismic trace：Geophysics，22，767-779)，并广泛用于多次波消除。常规预测反褶积算法采用了维纳滤波(Peacock，K.L.，andTreitel，S.，1969，Predictive deconvolution-theory and practice：Geophysics，34，155-169)，并假设地震信号中的一次波和多次波是正交的(Haykin，S.，1996，Adaptive filter theory：Prentice-Hall，Inc)。然而，实际地震信号很难满足这种假设。

假设y(n)是地震观测信号，而预测步长为q，q≥1。常规预测反褶积算法通过下式得到q步预测值，即

其中，a(i)，i＝1，…，m是滤波器系数。

预测信号

与观测信号y(n)的误差ε(n)表示如下：

$\mathrm{\ϵ}\left(n\right)=y\left(n\right)-\hat{y}\left(n\right)=y\left(n\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}a\left(i\right)y(n-q-i+1)---\left(2\right)$

实际上，误差信号ε(n)就是采用常规预测反褶积算法恢复的一次波信号。将常规预测反褶积算法简记为TPD。TPD算法通过最小化下列目标函数来获得滤波器系数(Robinson，E.，1957，Predictive decomposition of seismic trace：Geophysics，22，767-779)：

$E=\underset{n}{\mathrm{\Σ}}{\mathrm{\ϵ}}^2\left(n\right),---\left(3\right)$

而它的最小二乘解是：

a＝[Q^TQ]^-1Q^Tb，               (4)

其中

a＝[a(l)，…，a(m)]^T，        (6)

b＝[r(q)，…，r(q+m-1)]^T，    (7)

并且

r(i)＝E{y(n)y(n+i)}＝r(-i).   (8)

当求得a(i)，i＝1…m后，将其代入式(2)就可以得到估计的一次波信号。

然而，常规的预测反褶积算法需要假设地震信号中的一次波和多次波具有正交性，而实际地震数据通常都不会满足这个假设，这在很大程度上限制了预测反褶积算法的应用效果。1985年Walden用大量资料数据验证了实际地震信号具有非高斯的特性，更确切地说是具有超高斯性(Walden，A.，1985，Non-Gaussianreflectivity，entropy，and deconvolution：Geophysics，50，2862-2888)。利用地震数据的非高斯特性，人们实现了最小熵反褶积(Wiggins，R.，1978，Minimum entropy deconvolution：Geoexplor.，16，21-35)，基于独立成分分析的多次波自适应相减(Verschuur，D.，Berkhout，A.，Wapenaar，C.，1992，Adaptive surface-related multiple elimination：Geophysics，57，1166-1177Lu，W.，Luo，Y.，Zhao，B.，2004，Adaptive multiple subtraction usingindependent component analysis：Chinese J.Geophys.(in Chinese)，47，886-891Lu，W.，Mao，F.，2005，Adaptive Multiple Subtraction UsingIndependent Component Analysis：The Leading edge，24，282-284)，取得了较好的应用效果。

发明内容

本发明的目的在于提供一种预测地震波中的多次波和一次波信号的方法。

本发明提供的预测地震波中的多次波信号的方法，包括以下步骤：

1)用如下式(9)的目标函数优化滤波器系数a(i)；

max J(a(i))＝[E{G(ε₀(n))}-E{G(v(n))}]²，i＝1，2，…，m；    (9)

2)将步骤1)中得到的优化滤波器系数a(i)代入式(1)，得到多次波信号

式(9)和式(1)中，ε₀(n)是误差信号ε(n)的归一化信号，v(n)是标准高斯信号，G(·)是非线性函数，q是地震观测信号的预测步长，i是滤波器系数的序号，n是地震观测信号的序号。

所述G(·)为式(9)中的一个非线性函数，具体可为：

G(x)＝-exp(-0.5x²)    (10)

当所述G(·)为式(10)中的函数时，所述用式(9)的目标函数优化滤波器系数a(i)具体可为迭代求解如下方程组；

$R(i,s)=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y(n-s-q+1)\right],i=1,...,m,s=1,...,m,---\left(12\right)$

$b\left(i\right)=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y\left(n\right)\right],i=1,...,m.---\left(13\right)$

式(11)、式(12)和式(13)中，m为滤波器的阶数。

本发明还提供了一种预测地震波中的一次波信号的方法，是将多次波信号

代入式(2)得到一次波信号ε(n)；

$\mathrm{\ϵ}\left(n\right)=y\left(n\right)-\hat{y}\left(n\right)=y\left(n\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}a\left(i\right)y(n-q-i+1);---\left(2\right)$

式(2)中，所述y(n)为地震观测信号；所述

是通过上述方法计算得到的多次波。

本发明提出并实现了基于非高斯性最大化的预测地震波中的多次波和一次波信号的方法，特点是：

1)利用地震信号的非高斯特性。利用地震信号的非高斯性(Walden，A.，1985，Non-Gaussian reflectivity，entropy，and deconvolution：Geophysics，50，2862-2888)建立合适的目标函数，从而克服传统的预测反褶积算法需要假设一次波和多次波是正交的限制。

2)建立了新的基于非高斯性最大化的优化目标函数。

3)根据目标函数构造迭代的非线性方程组，通过多次迭代更新预测滤波器的系数，输出估计的一次波信号。

4)克服传统的预测反褶积算法需要假设一次波和多次波是正交的限制。

以下的实施例便于更好地理解本发明，但并不限定本发明。

附图说明

图1为本发明方法的计算机程序流程图。

图2为人工合成数据多次波压制：a)一次波；b)多次波+一次波；c)采用TPD算法得到的结果；d)采用PDMNG算法得到的结果。

图3为ISNR随参数λ变化的曲线。

具体实施方式

当G(x)＝-exp(-0.5x²)时，为了优化式(9)的目标函数，对J(a(i))关于a(i)，i＝1，…，m求偏导如下：

$\frac{\∂J\left(a\left(i\right)\right)}{\∂a\left(i\right)}=0.---\left(14\right)$

从而得到：

$2\left[E\right\{G\left({\mathrm{\ϵ}}_0\left(n\right)\right)\}-E\{G\left(v\left(n\right)\right)\left\}\right][\frac{\∂}{\∂a\left(i\right)}E\{G\left({\mathrm{\ϵ}}_0\left(n\right)\right)\}-\frac{\∂}{\∂a\left(i\right)}E\{G\left(v\left(n\right)\right)\left\}\right]=0.---\left(15\right)$

由于地震信号是非高斯分布的，同时E{G(v(n))}是一个常量，从而E{G(ε₀(n))}-E{G(v(n))}≠0，并且 $\frac{\∂}{\∂a\left(i\right)}E\left\{G\left(v\left(n\right)\right)\right\}=0.$ 式(15)可以改写成：

$\frac{\∂}{\∂a\left(i\right)}E\left\{G\left({\mathrm{\ϵ}}_0\left(n\right)\right)\right\}=0.---\left(16\right)$

按如下步骤改写式(16)：

$\frac{\∂}{\∂a\left(i\right)}\left[\frac{1}{N}\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}G\left({\mathrm{\ϵ}}_0\left(n\right)\right)\right]$

$=\frac{1}{N}\left[\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\frac{\∂G\left({\mathrm{\ϵ}}_0\left(n\right)\right)}{\∂a\left(i\right)}\right]$

$=\frac{1}{N}\left[\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left(\frac{\mathrm{dG}\left({\mathrm{\ϵ}}_0\left(n\right)\right)}{d{\mathrm{\ϵ}}_0\left(n\right)}\right)\frac{d{\mathrm{\ϵ}}_0\left(n\right)}{\mathrm{d\ϵ}\left(n\right)}\frac{\∂\mathrm{\ϵ}\left(n\right)}{\∂a\left(i\right)}\right].---\left(17\right)$

$=\frac{1}{N}\left[\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}g\left({\mathrm{\ϵ}}_0\left(n\right)\right)\frac{d{\mathrm{\ϵ}}_0\left(n\right)}{\mathrm{d\ϵ}\left(n\right)}\frac{\∂\mathrm{\ϵ}\left(n\right)}{\∂a\left(i\right)}\right]$

$=0$

因为

$\frac{\∂\mathrm{\ϵ}\left(n\right)}{\∂a\left(i\right)}=\frac{\∂}{\∂a\left(i\right)}[y\left(n\right)-\underset{t=1}{\overset{m}{\mathrm{\Σ}}}y(n-t-q+1)a\left(t\right)]$

$=-\underset{t=1}{\overset{M}{\mathrm{\Σ}}}\frac{\∂}{\∂a\left(i\right)}\left[y(n-t-q+1)a\left(t\right)\right],---\left(18\right)$

$=-y(n-i-q+1)$

同时

$\frac{d{\mathrm{\ϵ}}_0\left(n\right)}{\mathrm{d\ϵ}\left(n\right)}=\frac{d}{\mathrm{d\ϵ}\left(n\right)}\left(\frac{1}{\sqrt{\underset{z=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ϵ}}^2\left(z\right)}}\mathrm{\ϵ}\left(n\right)\right),---\left(19\right)$

$=\frac{1}{\sqrt{\underset{z=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ϵ}}^2\left(z\right)}}$

可以按如下步骤简化式    (17)

$\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[{\mathrm{\ϵ}}_0\left(n\right)\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)\right]=0,---\left(20\right)$

$\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[(y\left(n\right)-\underset{s=1}{\overset{m}{\mathrm{\Σ}}}y(n-s-q+1)a\left(s\right))\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)\right]=0,---\left(21\right)$

$\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\left(\underset{s=1}{\overset{m}{\mathrm{\Σ}}}y(n-s-q+1)a\left(s\right)\right)\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)\right]$

$=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[y\left(n\right)\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)\right]---\left(22\right)$

$\underset{s=1}{\overset{m}{\mathrm{\Σ}}}\left[a\left(s\right)\left(\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left(\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y(n-s-q+1)\right)\right)\right]---\left(23\right)$

$=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[y\left(n\right)\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)\right]$

重写式(23)成矩阵向量的形式：

其中

$R(i,s)=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y(n-s-q+1)\right],i=1,...,m,s=1,...,m,---\left(12\right)$

$b\left(i\right)=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y\left(n\right)\right],i=1,...,m.---\left(13\right)$

通过以上推导，在以下实施例中，采用如下方法预测地震波中的一次波信号：

假设y(n)是地震观测信号，而预测步长为q，q≥1。

1)向计算机输入以下设定量，并初始化迭代计数器：

零均值化的地震观测信号y(n)，n＝1，…，N，设当前迭代次数l＝1，最大迭代次数L；

2)采用零均值随机数作为滤波器系数a(i)；i＝1，...，m的初值；

3)计算机按以下步骤计算归一化的误差信号ε₀(n)：

①计算机按下式计算误差信号ε(n)，

$\mathrm{\ϵ}\left(n\right)=y\left(n\right)-\hat{y}\left(n\right)=y\left(n\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}a\left(i\right)y(n-q-i+1)$

②对误差信号ε(n)进行归一化，

${\mathrm{\ϵ}}_0\left(n\right)=\frac{\mathrm{\ϵ}\left(n\right)}{\left|\right|\mathrm{\ϵ}\left(n\right)\left|\right|};$

4)按以下步骤构造求取滤波器系数a(i)；i＝1，...，m的迭代方程组：

①计算机按下式计算迭代方程组的元素R(i，s)和b(i)，

$R(i,s)=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y(n-s-q+1)\right],i=1,...,m,s=1,...,m,$

$b\left(i\right)=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y\left(n\right)\right],i=1,...,m.$

②构造迭代方程组，

5)按以下步骤迭代更新滤波器系数：

①将滤波器系数a(i)；i＝1，...，m的初值代入步骤4)的①计算迭代方程组的元素R(i，s)和b(i)；

②利用最小二乘法求解步骤4)的②中的方程组来更新系数a(i)；i＝1，...，m；

③令l＝l+1，若l≤L，将更新后的滤波器系数代入步骤3)，重复进行步骤3)至5)的操作；

④若l＞L，停止迭代，得到滤波器系数a(i)；i＝1，...，m；

6)输出估计一次波信号 $\mathrm{\ϵ}\left(n\right)=y\left(n\right)-\hat{y}\left(n\right)=y\left(n\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}a\left(i\right)y(n-q-i+1).$

下述实施例中的实验方法，如无特殊说明，均为常规方法。

实施例1、仿真试验

在dell WS690工作站上进行了下面的试验，该工作站是双颗Intel四核至强CPU，主频2.33GHZ，内存8G。

首先将地震反射系数和Ricker子波相卷积生成一次波信号x(n)，反射系数序列μ(n)是利用Bernoulli-Gaussian序列(John，J.K.，and Jerry，M.M.，1982，Maximum likelihood detection and estimation of Bernoulli-Gaussian processes：IEEE Transactions on Information Theory，28，482-488)生成，定义如下

μ(n)＝r(n)q(n)    (24)

其中，r(n)是方差为C的高斯白噪声序列，而q(n)是Bernoulli序列，并满足

$\mathrm{Pr}\left(q\left(n\right)\right)=\left\{\begin{array}{c}1-\mathrm{\λ},q\left(n\right)=0\\ \mathrm{\λ},q\left(n\right)=1\end{array}\right.---\left(25\right)$

式中，λ是可以用来控制序列μ(n)的稀疏性的参数。

图2a显示了一次实验中的一次波信号x(n)，λ＝0.1。模拟的多次波包括1-5阶多次波，每一阶的多次波是由一次波x(n)与同一个地震子波卷积得到，鸣震周期为100ms，地震信号长度为750个采样点，时间采样率为2ms。将生成的多次波信号和一次波信号相加就得到了观测信号y(n)，可以表示为

y(n)＝IFFT{X(ω)[1-F(ω)exp(-jqω)]⁵}    (26)

其中，IFFT{·}表示反傅利叶变换，X(ω)是一次波信号x(n)的频谱，F(ω)是用来生成多次波信号的频域滤波器。图2b显示了观测信号y(n)。

二、实施本发明(PDMNG算法)

为了得到预测滤波器系数a(i)；i＝1，...，m，我们采用如下迭代的方法求解：

1)向计算机输入以下设定量，并初始化迭代计数器：

地震信号长度N＝750采样点，滤波器长度m＝20，预测周期q为50个采样点，设当前迭代次数l＝1，最大迭代次数L＝3，零均值化的地震观测信号y(n)，n＝1，…，N；

2)采用零均值随机数作为滤波器系数a(i)；i＝1，...，m的初值；

3)计算机按以下步骤计算归一化的误差信号ε₀(n)：

①计算机按下式计算误差信号ε(n)，

$\mathrm{\ϵ}\left(n\right)=y\left(n\right)-\hat{y}\left(n\right)=y\left(n\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}a\left(i\right)y(n-q-i+1);$

②对误差信号ε(n)进行归一化，

${\mathrm{\ϵ}}_0\left(n\right)=\frac{\mathrm{\ϵ}\left(n\right)}{\left|\right|\mathrm{\ϵ}\left(n\right)\left|\right|};$

4)按以下步骤构造求取滤波器系数a(i)；i＝1，...，m的迭代方程组：

①计算机按下式计算迭代方程组的元素R(i，s)和b(i)，

$R(i,s)=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y(n-s-q+1)\right],i=1,...,m,s=1,...,m,$

$b\left(i\right)=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y\left(n\right)\right],i=1,...,m.$

②构造迭代方程组，

5)按以下步骤迭代更新滤波器系数：

①将滤波器系数a(i)；i＝1，...，m的初值代入步骤4)的①计算迭代方程组的元素R(i，s)和b(i)；

②利用最小二乘法求解步骤4)的②中的方程组来更新系数a(i)；i＝1，...，m；

③令l＝l+1，若l≤L，将更新后的滤波器系数代入步骤3)，重复进行步骤3)至5)的操作；

④若l＞L，停止迭代，得到滤波器系数a(i)；i＝1，...，m；

6)输出估计一次波信号 $\mathrm{\ϵ}\left(n\right)=y\left(n\right)-\hat{y}\left(n\right)=y\left(n\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}a\left(i\right)y(n-q-i+1).$

得到的结果显示在图2d。

三、对比试验

用TPD算法对信号y(n)进行多次波压制，结果显示在图2c。

分别用常规预测反褶积算法(TPD)和基于非高斯性最大化的预测反褶积算法(PDMNG)对信号y(n)进行多次波压制，可以看出，PDMNG更好的恢复了一次波。

实施例2、PDMNG算法处理结果的信噪比

为了定量评价PDMNG算法在处理不同非高斯性地震信号时的性能，在下面的实验中，设定参数λ的取值范围为[0.1，1]，取值步长为0.1。对每一种参数λ的取值，独立实验300次，计算PDMNG算法处理结果相对于TPD算法处理结果所提高信噪比的平均值ISNR。图3显示了ISNR随参数λ的变化曲线，可以看出，当λ越小(即信号的非高斯性越大)，相对于TPD算法，采用PDMNG算法压制多次波提高的信噪比越大。

Claims (3)

1.一种预测地震波中的多次波信号的方法，包括以下步骤：

1)用如下式(9)的目标函数优化滤波器系数a(i)：

max J(a(i))＝[E{G(ε₀(n))}-E{G(v(n))}]²，i＝1，2，…，m    (9)；

2)将步骤1)中得到的优化滤波器系数a(i)代入下式(1)，得到多次波信号

式(9)和式(1)中，ε₀(n)是误差信号ε(n)的归一化信号，v(n)是标准高斯信号，G(·)是非线性函数，J(a(i))是用来衡量ε₀(n)的非高斯性程度的变量，q是地震观测信号的预测步长，i是滤波器系数的序号，n是地震观测信号的序号，m是滤波器的阶数，y(n-q-i+1)是地震观测信号；

所述G(·)为下式(10)中的函数：

G(x)＝-exp(-0.5x²)    (10)。

2.如权利要求1所述的方法，其特征在于：所述用式(9)的目标函数优化滤波器系数a(i)的方法为迭代求解如下方程组：

$R(i,s)=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y(n-s-q+1)\right],$ i＝1，…，m，s＝1，…，m(12)；

$b\left(i\right)=\underset{n=q+m}{\overset{N}{\mathrm{\Σ}}}\left[\exp (-0.5{\mathrm{\ϵ}}_0^2\left(n\right))y(n-i-q+1)y\left(n\right)\right],$ i＝1，…，m    (13)；

式(11)、式(12)和式(13)中，m为滤波器的阶数，N为地震观测信号y的长度，q为地震观测信号的预测步长。

3.一种预测地震波中的一次波信号的方法，是将多次波信号

代入下式(2)得到一次波信号ε(n)；

$\mathrm{\ϵ}\left(n\right)=y\left(n\right)-\hat{y}\left(n\right)=y\left(n\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}a\left(i\right)y(n-q-i+1)---\left(2\right);$

式(2)中，所述y(n)为地震观测信号；所述

是通过权利要求1或2所述的方法计算得到的。

Images