CN103376466B - 一种多次波压制方法 - Google Patents

Abstract

本发明提供了一种多次波压制方法，属于油气及煤层气地震勘探与开发领域。本发明方法包括：(1)选择地震道；(2)截取数据；(3)确定多次波的周期；(4)，定义目标多次波时窗；(5)计算目标多次波得滤波算子；(6)褶积与滤波；(7)自相关检测滤波结果；(8)迭代处理；(9)将原始地震数据中的其它所有地震道的地震道数据依次作为选出的地震道数据x(t)，对每个选出的地震道数据x(t)重复步骤(4)至步骤(8)。本发明方法针对每个多次波调节衰减系数的值，获得了满意的滤波效果，有效地压制了多次波；本发明方法适用于二维地震数据和三维地震数据。

Description

一种多次波压制方法

技术领域

本发明属于油气及煤层气地震勘探与开发领域，具体涉及一种多次波压制方法。

背景技术

压制多次波的方法可分为两大类：一类是基于有效波和多次波之间差异的滤波方法，简称为滤波方法，比较典型且经常被应用的方法包括预测反褶积、Radon变换、二维或三维傅里叶变换和聚束滤波等；另一类是基于波动方程的预测减去法，简称为波动方程预测减去方法，这种方法通过波动方程模拟实际地震波场或反演地震数据来预测多次波，之后把预测出来的多次波数据从原始地震数据中减去，这类方法中典型的方法包括波场外推法、反馈环法、反散射级数法和SRME方法等。

其中，滤波方法目前利用的是多次波和有效波之间的主要差异特征，即多次波的周期性和有效波与多次波之间的可分离性。这类方法是利用多次波和有效波之间在地震波波传播上的运动学方面(如传播时间、传播速度等)上的差异，通过各种变换手段，把时间空间域(称之为“旧域”)含有多次波的地震数据映射到其他特殊的域(称之为“新域”)。映射之后，有效波与多次波在新域会呈现出较之旧域更为明显的差异性，可以通过各种变换手段，把有效波和多次波分离开，之后就可以在地震数据中滤除多次波，最后再将“新域”中的数据进行一次反变换，就得到了常用的时间空间域中多次波被滤除后的地震数据体。

目前，主要有三种不同的波动方程多次波压制方法：波场外推法，即利用波场外推来模拟多次波；反馈环和反散射级数法，则是通过叠前反演来预测多次波。可以认为：波场外推法是模型驱动的，而反馈环和反散射级数法是数据驱动的。

常规预测反褶积方法是：根据地震数据自相关中多次波的衰减关系求出一个平均衰减系数，计算一个滤波算子，然后用此衰减系数和该滤波算子对所有多次波进行预测反褶积滤波处理。由于单个衰减系数并不符合每个多次反射波的实际衰减情况，结果对任何一个多次波都没有得到最好的压制效果。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种多次波压制方法，在计算出滤波算子后，针对每个多次波调节衰减系数的值，最后获得满意的滤波效果，有效地压制多次波。

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

一种多次波压制方法，所述方法包括以下步骤：

(1)，选择地震道：

从原始地震数据中选择一道含有有效波和多次波的地震道数据x(t)；所述原始地震数据包括多个地震道的数据，一个地震道的数据就是地震道数据；所述地震道数据x(t)是指所选地震道的全部样点值；所述t表示所选地震道的全部样点数；

(2)，截取数据：

对该地震道数据x(t)进行自相关计算得到x(t)的自相关结果；所述x(t)的自相关结果中的第一部分代表地震子波，后面周期性出现的孤立震动即为多次波；得到x(t)的自相关结果即完成了截取数据的工作；

(3)，确定多次波的周期：

将从步骤(2)得到的x(t)的自相关结果的零时刻开始到第一个孤立震动的最大值之间的时延设定为预测步长α，该预测步长α就是多次波的重复周期；

(4)，定义目标多次波时窗：

在该地震道数据x(t)上，将第一个多次波设为目标多次波，设定第一个多次波的时窗的起始点是有效波的第一个点，时窗的长度设为L，时窗的长度L大于1个多次波的重复周期α，但小于2个多次波的重复周期α；

(5)，计算目标多次波得滤波算子：

计算位于时窗内的地震道数据x(n)的自相关得到x(n)的自相关结果r(n)，将所述x(n)自相关结果r(n)转换成对称矩阵放在方程(1)的左边，而右边列矩阵为第α个元素开始到时窗结束这一段内的地震道数据的自相关结果；这时方程(1)中的未知部分即为滤波算子a(n)，即(a₀，a₁，a₂，…，a_n-1)；然后采用高斯主元素消去法求得滤波算子a(n)；所述n表示选出的地震道上位于时窗内的样点数；方程(1)中的点表示省略号；

$\left[\begin{array}{ccccc}{r}_0& r_1& r_2& ...& r_{n-1}\\ r_1& r_0& r_1& ...& r_{n-2}\\ r_2& r_1& r_0& ...& r_{n-3}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ r_{n-1}& r_{n-2}& r_{n-3}& ...& r_0\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ .\\ .\\ .\\ a_{n-1}\end{array}\right..=\left[\begin{array}{c}{r}_{\mathrm{\α}}\\ r_{\mathrm{\α}+1}\\ r_{\mathrm{\α}+2}\\ .\\ .\\ .\\ r_{\mathrm{\α}+n-1}\end{array}\right]---\left(1\right)$

(6)，褶积与滤波：

将时窗中的地震道数据x(n)与a(n)进行褶积，得到的褶积结果A即为该时窗中的目标多次波的预测值数组；从褶积结果A找到褶积结果A在该地震道数据中的起始位置和终了位置，所述起始位置是在时窗中的第α个样点处，所述终了位置是在时窗的结束位置处；将位于所述起始位置和终了位置之间的该地震道上的每个样点的样点值减去该样点对应的褶积结果A的值，其它位置处的样点值保持不变，就得到了消除了目标多次波干扰的新的地震道数据x2(t)；

(7)，自相关检测滤波结果：

计算x2(t)的自相关得到其自相关结果，从该自相关结果上观察目标多次波处的自相关幅值，所述自相关幅值的大小反映滤波效果，自相关幅值越低，滤波效果越好；

设衰减系数k，将衰减系数k与a(n)相乘后的值作为新的a(n)，重复步骤(6)和步骤(7)直到找到能获得最佳滤波效果的衰减系数k，此时完成对目标多次波的压制；

(8)，迭代处理

保持时窗长度L不变，将所述时窗下移一个预测步长α，将下一个多次波作为目标多次波，然后重复第(5)步到(8)步，直到处理完该地震道数据上的所有多次波，就得到经过多次波压制后的该地震道的新的地震道数据。

所述方法在步骤(8)后进一步包括步骤(9)：将原始地震数据中的其它所有地震道的地震道数据依次作为选出的地震道数据x(t)，对每个选出的地震道数据x(t)重复步骤(4)至步骤(8)，这样就完成对所述原始地震数据的多次波压制，得到新的地震数据；

所有地震道采用的预测步长均与步骤(3)中得到的预测步长α相同；所有地震道采用的时窗的长度L都是相同的。

所述步骤(1)中的样点值是指地震仪器接收人工地震时所接收到的数据，其是一个浮点数。

所述步骤(5)中，对于在地震道数据x(t)上所开的最后一个时窗，如果该时窗中的数据不能达到时窗的长度L，为了计算的统一性及计算机算法实现的可行性，将数据按时窗的长度L补齐，补齐的部分填充0值，即将方程(1)中的3个矩阵的尾部空缺元素都写为0。

所述步骤(7)中，衰减系数k为浮点数，其有效范围在0～5之间。

所述方法能够用于三维地震数据体，所述三维地震数据体是由多个线数据组成，每个所述线数据就相当于一个所述的原始地震数据，依次采用所述方法完成对每个线数据的处理后，就完成了整个三维地震数据体的多次波压制。

与现有技术相比，本发明的有益效果是：实际地震数据中的多次波振幅衰减规律并不符合理想的一维理论，常规预测反褶积计算一次滤波算子，从地震道的自相关结果中提取平均衰减系数，并将该衰减系数应用于滤波算子处理多次波组中的每个多次波，使得滤波效果不理想。而利用本发明提出的目标反褶积法计算出滤波算子后，针对每个多次波调节系数值，最后获得了满意的滤波效果，有效地压制了多次波。

附图说明

图1是本发明多次波压制方法的步骤框图。

图2是本发明多次波压制方法的实施例中的含多次波的原始地震数据。

图3是本发明多次波压制方法的实施例中的单道多次波地震数据。

图4是本发明多次波压制方法的实施例中的单道多次波自相关数据。

图5是本发明多次波压制方法的实施例中的预测步长的确定。

图6是本发明多次波压制方法的实施例中的第一个目标多次波的时窗范围。

图7是本发明多次波压制方法的实施例中的开窗次数的示意图。

图8是利用本发明多次波压制方法进行的第一个试验中的消除第一个多次波后得到的地震数据。

图9是利用本发明多次波压制方法进行的第一个试验中的消除第二个多次波后得到的地震数据。

图10是利用本发明多次波压制方法进行的第一个试验中的消除第三个多次波后得到的地震数据。

图11是利用本发明多次波压制方法进行的第一个试验中的消除第四个多次波后得到的地震数据。

图12是利用本发明多次波压制方法进行的第二个试验中的将预测步长设置成55ms后得到的地震数据。

图13是利用本发明多次波压制方法进行的第二个试验中的将预测步长设置成60ms后得到的地震数据。

图14是利用本发明多次波压制方法进行的第二个试验中的将预测步长设置成65ms后得到的地震数据。

图15是利用本发明多次波压制方法进行的第二个试验中的将预测步长设置成70ms后得到的地震数据。

图16是利用本发明多次波压制方法进行的第二个试验中的将时窗长度设置成120ms后得到的地震数据。

图17是利用本发明多次波压制方法进行的第二个试验中的将时窗长度设置成100ms后得到的地震数据。

图18是利用本发明多次波压制方法进行的第二个试验中的在模型中增加有效反射或不相关的噪音后得到的地震数据。

图19本发明中对于一个原始地震数据中的所有地震道数据使用相同的参数的示意图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

本发明一种多次波压制方法是一种改进的方法，创造性地提出了目标反褶积法，在设计滤波算子时考虑多次反射之间的真实振幅比，而不去理会所谓的一维理论。

本发明方法的基本原理是：在计算出滤波算子后，不是按从地震数据自相关中得到的平均衰减系数对每个多次波进行处理，而是在(确定)开预测反褶积时窗(后面提到的时窗都是所述预测反褶积时窗的简称)时一次只处理一个多次波(又称为多次反射波)(从图7看，在某些时窗中会包含到两个多次波，但这对处理不会有任何影响，因为本发明方法确立了多次波周期长度，也就是预测步长α，所以处理时不会干扰到后一个多次波的数据)，然后给求出的滤波算子乘一个合适的衰减系数k，再用它和时窗内的地震数据褶积得到该多次波的预测值。从地震道中减去该预测值后再做自相关检查滤波效果。本发明方法在一定的范围内试验若干个衰减系数k，当从自相关结果观察到的多次波剩余分量最小时，说明找到了最合适的衰减系数k。然后在消除了第一个多次波的地震数据上将时窗下移包含下一个多次波，并重复上述处理步骤，挨个地处理后面的多次波，直到所有多次波都处理完毕即完成了该地震数据的多次波压制。

如图1所示，所述方法包括以下步骤：

1，选择地震道：

从原始地震数据(如图2所示)中选择一道含有有效波和多次波干扰的地震道数据x(t)，所述原始地震数据包括多个地震道的数据，一个地震道对应的数据就是地震道数据。具体是通过现有的地震数据显示软件，能很容易的发现某些地震道有多次波，将这些地震道数据选出即可。这里的选出指得到此地震道的全部样点值。图2是含多次波地震道数据的多道显示，该例子中有四道，图中的第一个方框表示所有道的有效波，后面每个方框表示所有道的多次波，四个方框表示有四个多次波，图3是从图2中选出的某一道的单道显示(地震数据包括多个地震道数据，每一道对应的数据就是地震道数据)，图2和图3中的横、纵坐标是一样的，横坐标都是样点数，纵坐标都是计算的样点值(即地震仪器接收人工地震时所收到的数据，每个点叫“采样点”，它是一个浮点数。)大小，图2中的数值大小没有标上，图2的纵坐标表示地震道的道号。

2，截取数据

选择的地震道数据x(t)(即图3中的单道多次波地震数据)所示区域内的全部数据，其包含一次波和所有多次波。对这个地震道数据进行自相关计算，在数学上有相应的自相关函数可以使用，自相关(Autocorrelation)函数是描述随机信号一个时刻与另一个时刻的依赖关系，即研究t时刻与t+t′时刻两个随机变量的相关性。两个有界、可积函数f₁(x，y)和的乘积积分， $f_{11}(x,y)=\∫\underset{-\∞}{\overset{\∞}{\∫}}{f}_1(\mathrm{\ξ},\mathrm{\η})f_1^{*}(\mathrm{\ξ}-x,\mathrm{\η}-y)\mathrm{d\ξd\η}$ 称为函数f₁(x，y)的自相关函数。简称自相关。x(t)的自相关结果显示如图4所示，x(t)的自相关结果更能反映多次波的重复周期，因为地震道的自相关具有子波自相关的特征。x(t)的自相关结果的第一部分代表地震子波，后面周期性出现的孤立震动一般是多次波。

3，确定多次波周期

从x(t)的自相关结果上找出多次波的重复周期。将从x(t)的自相关结果的零时刻开始到第一个孤立震动的最大值之间的时延(如图5中的方框所示)设定为预测步长α，预测步长α也就是多次波的重复周期，如图5所示。软件设计在操作界面上提供预测步长参数的输入窗口。

4，选择目标多次波时窗

在图6所示的地震道上针对目标多次波(每个多次波依次作为目标多次波)开时窗。从单道中获得周期同样用于其它3个地震道，也就是说，在多道处理时用同样的周期参数进行多次波的压制。图6中设定了第一个多次波的时窗范围，其长度从一次反射的起点(第1次开窗的反射起点是有效波的起点，即有效波的第一个点)开始，时窗长度＞1个多次波周期，但＜2个多次波周期。后面所有多次波的时窗都与第一个多次波的时窗大小一样。

5，计算目标多次波的滤波算子

计算第一个时窗中的地震道数据x(n)(x(t)是整个地震道数据，x(n)是时窗中的地震道数据，即整个地震道数据的一部分。)的自相关得到自相关结果r(n)(一个地震道有很多样点组成，t表示整道地震数据的样点数，n表示所开时窗内地震数据的样点数)。将自相关结果r(n)(步骤2中是用整道地震数据自相关来确定多次波的周期，这里是将时窗中的地震数据自相关来压制多次波)转换成对称矩阵放在方程(1)的左边，而右边列矩阵为第α个元素(此处的α就是前面所说的多次波的周期)开始到第一个时窗结束这一段内的地震道数据的自相关结果，实际上等于原始数据与预测输出的互相关，原始数据(公式(1)的最左边)与最右边互相关的结果就是a(n)。

这时方程(1)中的未知部分即为滤波算子a(n)，即(a₀，a₁，a₂，…，a_n-1)。采用高斯主元素消去法来求出滤波算子a(n)。对于在地震道数据上所开的最后一个时窗(如图7中的第5次开窗形成的时窗)，有可能其中的数据不能达到时窗的长度，为了计算的统一性及计算机算法实现的可行性，一般都将数据按时窗长度补齐，补齐的部分可以填充0值，即方程(1)中的3个矩阵的尾部空缺元素都写为0(因为公式(1)是针对所有时窗的，而只有时窗中的数据长度不够时才补充0，所以在公式(1)中没有体现出这些0来。)，这样并不会影响整体数据的正确性。

高斯主元素消去法(请参考《用高斯主元素消去法解线性方程组》、《达州职业技术学院学报》、2004年1期)是顺序消去法的一种改进，它的基本思想是在逐次消元时总是选绝对值最大的元素(称之为主元)做除数，按顺序消去法的步骤消元，式中对称矩阵为已知，r_α也为已知，用高斯主元素消去法即可求出a(n)。

$\left[\begin{array}{ccccc}{r}_0& r_1& r_2& ...& r_{n-1}\\ r_1& r_0& r_1& ...& r_{n-2}\\ r_2& r_1& r_0& ...& r_{n-3}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ r_{n-1}& r_{n-2}& r_{n-3}& ...& r_0\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ .\\ .\\ .\\ a_{n-1}\end{array}\right..=\left[\begin{array}{c}{r}_{\mathrm{\α}}\\ r_{\mathrm{\α}+1}\\ r_{\mathrm{\α}+2}\\ .\\ .\\ .\\ r_{\mathrm{\α}+n-1}\end{array}\right]---\left(1\right)$

该步骤中的逐次计算出每个时窗算子就可以消除现有多次波压制技术用通用性函数对整道数据做多次波处理时带来的误差，也就是说在用本发明方法处理多次波时已考虑了多次波之间的振幅差异了。

6，褶积与滤波

将输入x(n)与a(n)褶积，得到第一个多次波的预测值。从输入x(t)中减去褶积结果，就得到消除了第一个多次波干扰的地震道数据x2(t)。x(t)表示原始地震道数据，x(n)与a(n)褶积结果A就是某一段时窗中得到的一个多次波的预测值数组，从褶积结果A可以得到这段A在此地震道中的起始位置和终了位置，将这段原始地震道数据的每个样点值减去A的每个样点值即可。从第一个多次波的开始处，即第1次开窗中的第α个样点减起。从图7中可以看出第一个多次波的开始处就是第1次开窗中的第α个样点处，也就是第2次开窗的起始处，而第二个多次波的开始处就是第2次开窗中的第α个样点处，也就是第3次开窗的起始处，以此类推。

7，自相关检测滤波结果

重新进行x2(t)的自相关，从自相关的输出上观察第一个多次波处的自相关幅值，幅值的大小反映滤波效果。自相关输出幅值越低，滤波效果越好。可以试验不同的衰减系数k(为浮点数，在0～5之间为有效范围)与a(n)相乘后再与x(n)褶积，看对结果的影响。再重复本步骤的计算，直到找到获得最佳滤波效果的衰减系数k，此时完成对第一个多次波的压制。软件设计提供衰减系数k的输入窗口。

8，迭代处理

对第一个多次波获得满意的压制效果后，将时窗下移一个预测步长α，针对第二个多次波，再重复第5到7步的处理，直到整道多次波组结束。如图7所示，第1次开窗至第5次开窗所开的时窗的大小都是一样的，不同的就是每个相邻两个时窗相差一个预测步长α。因此，开取时窗要从上一个多次波的起点开始，时窗大小与上一个时窗一样，具体如图7所示，每次所开时窗的长度不变，但后一个时窗的起始位置并不是前一个时窗的末端，前后相邻两个时窗是有部分重叠的，目的是为了保证数据的连续性。而后一个时窗的起始位置最好是在上一个多次波波形起伏的第一个点，偏差不能多。因为每次是下移一个预测步长α，所以减去褶积结果时，要从本次包含的多次波开始处减起，即每个时窗中的第α个样点处减起。

9，结果数据

依次对所选择的地震道数据中的所有多次波做完以上步骤的处理后，即得到了压制多次波后的该地震道的地震道数据，对于原始地震数据中的其它所有道的地震道数据再利用以上步骤所得到的参数(如图19所示。)，进行循环计算即可，这样就完成对所述原始地震数据的多次波压制，得到新的地震数据。

本发明方法适用于二维地震数据也同样适用于三维地震数据，对于三维地震数据体实际上可以采用二维地震数据的方式来完成，三维地震数据体分为若干的线数据，线数据就相当于一个二维地震剖面，依次用本发明的方法完成这些线后，即相当于完成了整个三维数据体。计算方式与二维完全一样，就不再附图说明了。

下面通过几个试验来说明本发明的效果：下面的试验用到的软件是针对本发明方法编写的验证软件，不是任何商业软件或其他软件，以图8为例说明该验证软件的界面，软件界面中最上面的窗口显示的是原始地震数据(是一道地震数据)，在“窗口开取长度”处可设置所开时窗的大小，“道截取时间长度”表示图中显示的原始地震数据的长度；界面中的第二个窗口显示的是原始地震数据自相关后的结果，在“多次波周期起始位置”处可设置多次波的起始样点位置，在“多次波周期结束位置”处可设置单个多次波的结束样点位置；界面中的第三个窗口显示的是每一步针对目标多次波执行预测反滤波后的自相关结果，用户可以从此窗口评估旋转算子参数的效果；界面中的第4个窗口显示的是消除每个多次波后的地震道数据，因为软件设计成可以针对每个多次波分布执行预测反滤波，所以这里设计有该窗口；界面中在第4个窗口下方还为用户提供了操作控制旋转，软件允许执行单步调试，也可以执行连续处理，如果进行单步调试，确认“启动单步调试”的选择框，选择单步调试后对话框显示需要进行处理的“总步数”吗，如果进行连续处理，则确认“连续处理”的选择框。总步数是根据截取的地震道长度自动确定的，所以在截取的地震道数据时一般应包含整个多次波波组；在“系统设置”参数输入区输入衰减系数k并按“执行本步”按钮，程序将对当前时窗执行目标反褶积处理；自相关窗口和滤波后数据窗口将显示处理后的结果；如果对处理结果不满意，可以输入新的衰减系数k后再按“执行本步”按钮，如此循环直至获得满意的滤波效果为止。

1，利用本发明方法对非典型(即不按一维规律衰减的)多次波模型数据的试验：

用本发明方法对不按一维规律衰减的多次波模型数据进行处理。采用单步执行(单步执行就是指每次压制一个多次波)的方式，分别设置衰减系数k为1.0、0.94、1.1和0.92。图8、图9、图10和图11是执行每一步的显示结果。从这些图中可以看到，每一个多次波都得到了很好的压制(又称为消除或滤波)效果。

2，预测步长、时窗长度和噪音对本发明方法的影响

分别改变预测步长、时窗长度、和在模型中添加噪音(或有效反射)对本发明方法进行了测试。

其中，图12，图13，图14和图15显示的是依次完成的多次波压制效果，分别将预测步长(即软件界面中的“多次波周期结束位置”)设定为55ms，60ms，64ms和70ms。从测试结果看到，预测步长应该严格地等于多次波的重复周期(即第一个测试中的所用的64ms)。否则将得不到好的滤波效果。

图16，图17是对同样模型数据分别将时窗长度(即软件界面中的“窗口开取长度”)设置成120ms和100ms处理的结果。可以看出，只要时窗长度大于一个多次波周期但小于2个多次波周期，都能取得满意的多次波滤波效果。

最后图18的地震数据中，增加了3个与多次波周期不相干的同相轴，它们可以是有效反射波或噪音。结果表明本发明方法只能滤除多次波，对这些同相轴不起作用。图18是加入噪音后压制多次波的效果示意，最后留下来的不是多次波，而是噪音。图18表明本方法只压制多次波，而不会影响其它不相关的数据。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (5)

1.一种多次波压制方法，其特征在于：所述方法包括以下步骤：

(1)，选择地震道：

从原始地震数据中选择一道含有有效波和多次波的地震道数据x(t)；所述原始地震数据包括多个地震道的数据，一个地震道的数据就是地震道数据；所述地震道数据x(t)是指所选地震道的全部样点值；所述t表示所选地震道的全部样点数；

(2)，截取数据：

对该地震道数据x(t)进行自相关计算得到x(t)的自相关结果；所述x(t)的自相关结果中的第一部分代表地震子波，后面周期性出现的孤立震动即为多次波；

(3)，确定多次波的周期：

将从步骤(2)得到的x(t)的自相关结果的零时刻开始到第一个孤立震动的最大值之间的时延设定为预测步长α，该预测步长α就是多次波的重复周期；

(4)，定义目标多次波时窗：

在该地震道数据x(t)上，将第一个多次波设为目标多次波，设定第一个多次波的时窗的起始点是有效波的第一个点，时窗的长度设为L，时窗的长度L大于1个多次波的重复周期α，但小于2个多次波的重复周期α；

(5)，计算目标多次波的滤波算子：

计算位于时窗内的地震道数据x(n)的自相关得到x(n)的自相关结果r(n)，将所述x(n)自相关结果r(n)转换成对称矩阵放在方程(1)的左边，而右边列矩阵为第α个元素开始到时窗结束这一段内的地震道数据的自相关结果；这时方程(1)中的未知部分即为滤波算子a(n)，即(a₀,a₁,a₂,…,a_n-1)；然后采用高斯主元素消去法求得滤波算子a(n)；所述n表示选出的地震道上位于时窗内的样点数；

$\left[\begin{array}{ccccc}{r}_0& r_1& r_2& ...& r_{n-1}\\ r_1& r_0& r_1& ...& r_{n-2}\\ r_2& r_1& r_0& ...& r_{n-3}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ r_{n-1}& r_{n-2}& r_{n-3}& ...& r_0\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ .\\ .\\ .\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}{r}_{\mathrm{\α}}\\ r_{\mathrm{\α}+1}\\ r_{\mathrm{\α}+2}\\ .\\ .\\ .\\ r_{\mathrm{\α}+n-1}\end{array}\right]---\left(1\right)$

(6)，褶积与滤波：

将时窗中的地震道数据x(n)与a(n)进行褶积，得到的褶积结果A即为该时窗中的目标多次波的预测值数组；从褶积结果A找到褶积结果A在该地震道数据中的起始位置和终了位置，所述起始位置是在时窗中的第α个样点处，所述终了位置是在时窗的结束位置处；将位于所述起始位置和终了位置之间的该地震道上的每个样点的样点值减去该样点对应的褶积结果A的值，其它位置处的样点值保持不变，就得到了消除了目标多次波干扰的新的地震道数据x2(t)；

(7)，自相关检测滤波结果：

计算x2(t)的自相关得到其自相关结果，从该自相关结果上观察目标多次波处的自相关幅值，所述自相关幅值的大小反映滤波效果，自相关幅值越低，滤波效果越好；

设衰减系数k，将衰减系数k与a(n)相乘后的值作为新的a(n)，重复步骤(6)和步骤(7)直到找到能获得最佳滤波效果的衰减系数k，此时完成对目标多次波的压制；

其中，衰减系数k为浮点数，在0～5之间为有效范围；

(8)，迭代处理：

保持时窗长度L不变，将所述时窗下移一个预测步长α，将下一个多次波作为目标多次波，然后重复第(5)步到(8)步，直到处理完该地震道数据上的所有多次波，就得到经过多次波压制后的该地震道的新的地震道数据。

2.根据权利要求1所述的多次波压制方法，其特征在于：所述方法在步骤(8)后进一步包括步骤(9)：将原始地震数据中的其它所有地震道的地震道数据依次作为选出的地震道数据x(t)，对每个选出的地震道数据x(t)重复步骤(4)至步骤(8)，这样就完成对所述原始地震数据的多次波压制，得到新的地震数据；

所有地震道采用的预测步长均与步骤(3)中得到的预测步长α相同；所有地震道采用的时窗的长度L都是相同的。

3.根据权利要求1或2所述的多次波压制方法，其特征在于：所述步骤(1)中的样点值是指地震仪器接收人工地震时所接收到的数据，其是一个浮点数。

4.根据权利要求1或2所述的多次波压制方法，其特征在于：所述步骤(5)中，对于在地震道数据x(t)上所开的最后一个时窗，如果该时窗中的数据不能达到时窗的长度L，为了计算的统一性及计算机算法实现的可行性，将数据按时窗的长度L补齐，补齐的部分填充0值，即将方程(1)中的3个矩阵的尾部空缺元素都写为0。

5.根据权利要求2所述的多次波压制方法，其特征在于：所述方法能够用于三维地震数据体，所述三维地震数据体是由多个线数据组成，每个所述线数据就相当于一个所述的原始地震数据，依次采用所述方法完成对每个线数据的处理后，就完成了整个三维地震数据体的多次波压制。