一种层间多次波压制方法与流程

本发明涉及海洋和陆地地球物理勘探技术领域，尤其涉及一种层间多次波压制方法。

背景技术：

在海洋和陆地地震勘探中，由于地下强反射界面的存在，地震波在强反射界面间多次反射形成层间多次波，层间多次波与一次波相互叠加干扰，严重降低了地震资料的分辨率，加大了对有效波识别的难度，影响了地震成像质量和地震解释的真实性和可靠性。因此，衰减或消除层间多次波是地震资料处理中一个重要环节。

为了消除层间多次波的干扰，提高资料分辨率，物探界提出了两类多次波压制方法：

(1)基于一次波与多次波之间特性差异的滤波法。滤波法包括预测反褶积法、f-k滤波法、radon变换法、聚束滤波法等。当假设条件较好地满足时，滤波法能有效地衰减或消除多次波，并且计算量小，易于实现，效率高，但是滤波法需要较多的地下假设信息，当一次波和多次波之间的特征差异很小或没有时，很难获得理想的效果，甚至会严重损伤一次波。

(2)基于波动理论的预测相减法。预测相减法避免了滤波法的局限性，无需先验信息，是多次波压制方法的主要发展趋势。预测相减法主要包括反馈迭代法和逆散射级数法。针对层间多次波压制，反馈迭代法需要人工干预，通过一层一层估计多次波产生算子来预测层间多次波，而逆散射级数法完全数据驱动，无需人工干预，通过算法自身进行预测，一次可以预测出所有的层间多次波，是目前最先进的层间多次波压制方法。

现有的逆散射级数法的不足在于，由于预测层间多次波过程中至少有三个子波进行褶积，所以预测的层间多次波波形和振幅发生了很大变化，同时使同相轴变“胖”并使能量差异变大，因此子波的存在对层间多次波预测的准确性有很大的影响。

技术实现要素：

本发明针对上述现有技术的不足，提供了一种层间多次波压制方法，包括以下步骤：

对当前采集到的地震数据进行预处理，并对预处理后的地震数据进行子波去除处理，得到去除子波后的地震数据；

对所述去除子波后的地震数据进行抽道集，得到全波场数据；

对所述全波场数据进行傅里叶变换，得到频率-波数域数据；

对所述频率-波数域数据进行常数背景速度偏移，得到伪深度域数据；

根据所述伪深度域数据，利用逆散射级数法进行层间多次波预测，得到去除子波后的层间多次波数据；

对所述去除子波后的层间多次波数据进行反傅里叶变换，得到空间-时间域数据；

对所述空间-时间域数据进行子波补偿处理，得到补偿子波后的层间多次波数据；

将所述预处理后的地震数据与所述补偿子波后的层间多次波数据进行匹配相减，得到层间多次波压制后的地震数据。

在一个实施例中，利用自适应相减法从所述预处理后的地震数据中减去所述补偿子波后的层间多次波数据，得到层间多次波压制后的地震数据。

在一个实施例中，所述对当前采集到的地震数据进行预处理，包括对海洋地震数据进行预处理和对陆地地震数据进行预处理。

在一个实施例中，对当前采集到的海洋地震数据进行包括低切滤波、去气泡、压制涌浪噪音、去除直达波、压制鬼波和压制自由表面多次波的预处理。

在一个实施例中，对当前采集到的陆地地震数据进行包括静校正、异常振幅相干噪音压制、地表一致性处理、剩余振幅补偿和剩余静校正的预处理。

在一个实施例中，所述频率-波数域数据为去除子波后的单频平面波场数据，其通过以下表达式表示：

b1(kg，ks，qg+qs)＝-2iqsa(ω)^-1d(kg，ks，ω)

其中，b1(kg，ks，qg+qs)表示去除子波后的单频平面波场数据，i表示虚数单位，ks和kg分别表示震源和检波器的水平方向波数，qs和qg分别表示相应的震源和检波器的垂向波数，a(ω)表示子波，ω表示圆周频率，d(kg，ks，ω)表示地震数据。

在一个实施例中，qs满足qg满足c0为背景介质速度。

在一个实施例中，通过以下表达式进行层间多次波预测，得到去除子波后的层间多次波数据；

其中，b3(kg，ks，qg+qs)表示去除子波后的层间多次波数据，k1和k2分别表示相应震源和检波器的水平方向波数积分变量，q1和q2分别表示相应震源和检波器的垂向波数变量，εs和εg分别表示震源和检波器的深度，ε表示使约束关系z1＞z2和z3＞z2严格成立的常量，zj＝1，2，3表示常数背景速度域成像的伪深度，b1(k2，ks，zj)表示伪深度域数据。

在一个实施例中，ε的值与子波的长度有关。

在一个实施例中，所述补偿子波后的层间多次波数据通过以下表达式表示：

m(kg，ks，ω)＝(-2iqs)^-1a(ω)b3(kg，ks，qg+qs)

其中，m(kg，ks，ω)表示补偿子波后的层间多次波数据。

与现有技术相比，本发明的一个或多个实施例可以具有如下优点：

本发明的层间多次波压制方法，通过增加子波去除和子波补偿步骤，可以有效提高层间多次波预测的准确度，从而压制层间多次波，提高地震资料分辨率，并改善成像品质。

本发明的其它特征和优点将在随后的说明书中阐述，并且部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例共同用于解释本发明，并不构成对本发明的限制。在附图中：

图1为现有的正演反射波和层间多次波散射模型示意图；

图2为现有的反演中三个散射点不同位置关系组成的波示意图；

图3为本发明实施例的层间多次波压制方法的流程图；

图4为本发明实施例的模拟数据与层间多次波预测结果对比分析示意图；

图5为本发明实施例的模拟数据中的层间多次波与预测结果对比示意图；

图6为本发明实施例的海洋数据层间多次波压制前后对比示意图；

图7为本发明实施例的海洋数据层间多次波压制前后放大图；

图8为本发明实施例的陆地数据层间多次波压制前后炮集图；

图9为本发明实施例的陆地数据层间多次波压制前后叠加剖面图。

具体实施方式

以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。

为了更清楚地理解本发明，下面首先对现有的逆散射级数预测层间多次波的物理机制进行分析。

散射理论是关于散射波场(实际波场与背景波场之间的差异)和扰动介质(实际介质与背景介质之间的差异)之间关系的基本理论。实际波场和背景波场分别满足两个基本波动方程：

lg＝δ(1)

l0g0＝δ(2)

其中，l和l0分别为实际介质和背景介质中的微分算子，g和g0分别为实际介质和背景介质中的波场，δ为狄拉克算子。

定义扰动算子v≡l0-l和散射波场ψs≡g-g0，散射方程可以表达为：

g＝g0+g0vg(3)

把方程(3)代入到自身中，得到正演散射级数

g＝g0+g0vg0+g0vg0vg0+g0vg0vg0vg0+…(4)

已知介质信息(实际介质与背景介质之间的差异v)和背景波场g0，通过公式(4)即可求得实际波场g，同时也可获得散射波场ψs＝g-g0，即

这里(ψs)n是v的n阶方程式，为散射波场ψs的一部分。公式(4)和(5)在任意位置成立，当散射波场在测量表面时称之为地震数据d＝(ψs)ms。扰动算子v可以展开成级数形式，得到

v＝v1+v2+v3+…(6)

把公式(6)代入到公式(5)中，在测量表面令方程两端同阶相等得到逆散射级数

d＝g0v1g0(7)

g0v2g0＝-g0v1g0v1g0(8)

由公式(7)-(9)可知，理论上可通过已知的实际波场、背景波场和背景介质一步一步求得v1，v2，v3等，并最终获得扰动算子v和实际介质信息，但是逆散射级数的总体收敛性并没有得到验证，因此可以选取一些子级数来处理一部分针对性问题，例如去除表面多次波子级数、去除层间多次波子级数、成像子级数和反演子级数等。

本发明实施例只涉及去除层间多次波子级数，因此先从正演(公式(5))的角度分析点散射模型反射波和层间多次波产生的不同物理机制。层间多次波的定义是在自由表面以下的任何一个界面，至少发生一次下行反射的波，下行反射的次数称为层间多次波的阶数。

图1为现有的正演反射波和层间多次波散射模型示意图。如图1所示，一阶层间多次波发生了一次下行反射，二阶层间多次波发生了两次下行反射。此外，下行反射的散射点深度要低于相邻两个散射点的深度，而且层间多次波同散射点的个数有关，至少需要三个散射点才能形成层间多次波。

通过正演可知，一阶和二阶散射级数无法形成层间多次波。同理，在反演中，一阶和二阶的逆散射级数也不可能形成层间多次波。将公式(9)简化如下：

其中，v31＝v1g0v2，v32＝v2g0v1，v33＝v1g0v1g0v1。由公式(10)可知，层间多次波从三阶逆散射级数中的第三项才开始产生，即至少有三个散射点。

图2为现有的反演中三个散射点不同位置关系组成的波示意图。图2示出了三阶逆散射级数第三项的不同组合，与正演类似，只有满足空间位置条件(z1＞z2和z3＞z2)和阶数限制条件的逆散射级数才能形成层间多次波。

发明人在分析了逆散射级数预测层间多次波的物理机制之后，从数学上推导了现有的逆散射级数预测层间多次波的算法，具体如下。

无限空间常速度背景介质中的格林函数g0，可以表示为

其中，i表示虚数单位，q表示垂向波数，q′表示垂向波数积分变量，z表示深度，z′表示另一个深度，ω表示圆周频率。

将右端积分部分分成一个主值和在奇点q′＝±q处的函数值，即

其中，pv表示主值，δ表示狄拉克算子。

只考虑在奇点处的贡献，把公式(12)右边第二项代入到逆散射级数中的三阶方程公式(10)中右边第三项可得(省略变量ω)

式中ks和kg分别表示震源和检波器的水平方向波数，qs和qg分别表示相应的震源和检波器的垂向波数，k和k′表示相应震源和检波器的水平方向波数积分变量，q和q′分别表示相应震源和检波器的垂向波数变量。此处只考虑奇点处的贡献，而不考虑主值的贡献，其物理意义在于后者在逆散射级数中起着异常速度扰动的作用，同速度模型有关，前者在逆散射级数中只同背景速度模型有关。

定义单频平面波场数据b1为

根据公式(13)和(14)，引入z1＞z2且z3＞z2可构造得出层间多次波的预测算法，

其中，εs和εg为震源和检波器的深度，qs和qg为震源和检波器的垂向波数，qs满足qg满足k1和k2分别表示相应震源和检波器的水平方向波数积分变量，q1和q2分别表示相应震源和检波器的垂向波数变量，c0为背景介质速度。zj＝1，2，3为常数背景速度域成像的伪深度。ε表示使约束关系z1＞z2和z3＞z2严格成立的常量，在实际数据处理中它的值与子波的长度有关。输入数据中如果三个一次波满足“低-高-低”约束关系将产生一阶层间多次波，如果一阶层间多次与两个一次波满足“低-高-低”约束关系就形成二阶层间多次波，同理高阶层间多次波也可以通过低阶层间多次波加以合成。由于输入数据中含有所有信号包括一次波和各阶层间多次波，所以该方法可以一次预测出所有的层间多次波。

公式(15)是现有的逆散射级数层间多次波预测算法，但是由于预测层间多次波过程中至少有三个子波进行褶积，所以预测的层间多次波波形和振幅发生了很大变化，同时使同相轴变“胖”并使能量差异变大，因此子波的存在对层间多次波预测的准确性有很大的影响。

针对上述现有的逆散射级数层间多次波预测算法的不足，本发明实施例从散射理论出发，引入了子波的影响，提出了一种改进的逆散射级数层间多次波压制方法。

图3为本发明实施例的层间多次波压制方法的流程图。如图3所示，可以包括以下步骤s310至s380。

在步骤s310中，对当前采集到的地震数据进行预处理，并对预处理后的地震数据进行子波去除处理，得到去除子波后的地震数据。

具体地，对当前采集到的地震数据进行预处理，包括对当前采集到的海洋地震数据进行预处理和对当前采集到的陆地地震数据进行预处理。对海洋数据的预处理包括低切滤波、去气泡、压制涌浪噪音、去除直达波、压制鬼波和压制自由表面多次波等，对陆地数据的预处理包括静校正、异常振幅相干噪音压制、地表一致性处理、剩余振幅补偿和剩余静校正等。当前采集到的地震数据通过预处理后，首先进行子波去除处理，得到去除子波后的地震数据。

在步骤s320中，对去除子波后的地震数据进行抽道集，得到全波场数据。

在步骤s330中，对全波场数据进行傅里叶变换，得到频率-波数域数据。

频率-波数域数据为去除子波后的单频平面波场数据，其通过以下表达式表示：

b1(kg，ks，qg+qs)＝-2iqsa(ω)^-1d(kg，ks，ω)(16)

其中，b1(kg，ks，qg+qs)表示去除子波后的单频平面波场数据，i表示虚数单位，ks和kg分别表示震源和检波器的水平方向波数，qs和qg分别表示相应的震源和检波器的垂向波数，a(ω)表示子波，ω表示圆周频率，d(kg，ks，ω)表示地震数据。优选地，qs满足qg满足c0为背景介质速度。

与表达式(14)相比，本发明实施例一中的单频平面波场数据去除了子波，得到了去除子波后的频率-波数域数据。

在步骤s340中，对频率-波数域数据进行常数背景速度偏移，得到伪深度域数据。

在步骤s350中，根据伪深度域数据，利用逆散射级数法进行层间多次波预测，得到去除子波后的层间多次波数据。

通过以下表达式进行层间多次波预测，得到去除子波后的层间多次波数据；

其中，b3(kg，ks，qg+qs)表示去除子波后的层间多次波数据，k1和k2分别表示相应震源和检波器的水平方向波数积分变量，q1和q2分别表示相应震源和检波器的垂向波数变量，εs和εg分别表示震源和检波器的深度，ε表示使约束关系z1＞z2和z3＞z2严格成立的常量，zj＝1，2，3表示常数背景速度域成像的伪深度，b1(k2，ks，zj)表示伪深度域数据。其中，ε的值与子波的长度有关。

在步骤s360中，对去除子波后的层间多次波数据进行反傅里叶变换，得到空间-时间域数据。

在步骤s370中，对空间-时间域数据进行子波补偿处理，得到补偿子波后的层间多次波数据。

补偿子波后的层间多次波数据通过以下表达式表示：

m(kg，ks，ω)＝(-2iqs)^-1a(ω)b3(kg，ks，qg+qs)(18)

其中，m(kg，ks，ω)表示补偿子波后的层间多次波数据。

在步骤s380中，将预处理后的地震数据与补偿子波后的层间多次波数据进行匹配相减，得到层间多次波压制后的地震数据。

优选地，利用自适应相减法从预处理后的地震数据中减去补偿子波后的层间多次波数据，得到层间多次波压制后的地震数据。

下面通过对模拟数据和实际地震资料的处理，来对本实施例的层间多次波压制方法进行验证。

图4为本发明实施例的模拟数据与层间多次波预测结果对比分析示意图。其中图4包括图4(a)-(i)。

首先，采用简单层状模型产生模拟数据，该模型有三层，层速度分别为2000m/s、3200m/s和6100m/s，第一层深度为300m，第二层深度为700m。震源和检波器在表面，震源采用主频25hz的雷克子波，共2001道，道间距为5m，采样间隔为0.005s，记录时间为1.5s。图4(a)为有限差分法产生的模拟数据，该模拟数据含有两个一次波和一阶、二阶层间多次波。

其次，用改进的逆散射级数法预测层间多次波，也即采用图3所示的层间多次波压制方法预测层间多次波。具体地，通过地震数据预处理-子波去除-抽道集-傅里叶变换-常数背景速度偏移-层间多次波预测-反傅里叶变换-子波补偿，预测得到如图4(b)所示的层间多次波。用现有的逆散射级数法预测层间多次波。具体地，通过地震数据预处理-抽道集-傅里叶变换-常数背景速度偏移-层间多次波预测反傅里叶变换，预测得到如图4(c)所示的层间多次波。

可以看出，改进的逆散射级数法预测的层间多次波与模拟数据中的层间多次波能量相近，而现有的逆散射级数法预测的层间多次波波形变“胖”且与模拟数据中的层间多次波相比能量差距较大。

为了进一步对比波形和振幅，从图4(a)-(c)中选取近偏移(图4(d)-(f))和远偏移距(图4(g)-(i))数据进行比较。图4(d)为原始数据中近偏移距的层间多次波波，图4(e)和图4(f)分别为用改进的和现有的逆散射级数法预测层间多次波的近偏移距道数据。

图5为本发明实施例的模拟数据中的层间多次波与预测结果对比示意图。其中，图5(a)和图5(b)分别为模拟数据中的层间多次波与用改进的逆散射级数法预测的层间多次波在近偏移距和远偏移距的对比。

可以看出，改进的逆散射级数法预测的层间多次波与模拟数据中的层间多次波波形和能量基本一致，如图5(a)，而现有的逆散射级数法预测的层间多次波波形差异较大，如图4(f)。同样，在远偏移距数据(图4(g)-(i))对比中也可得出一样的结论，如图5(b)。

因此，通过模型数据测试验证了改进的逆散射级数法提高了层间多次波预测的准确性。

下面通过选取一套深海海洋数据和一套陆地数据对逆散射级数层间多次波压制方法进行分析和验证。

海洋数据震源深度为8.5m，检波器深度9.5m，采样间隔0.004s，记录长度9s，每炮239道，共1751炮，道间距25m，炮间距50m，最小偏移距160m。

图6为本发明实施例的海洋数据层间多次波压制前后对比示意图。其中，图6(a)为经过预处理和表面多次波压制后的数据。由于层间多次波算法的计算量特别大，选取我们关心的岩底结构即方框内的部分进行层间多次波压制。图6(b)为层间多次波压制后的结果。

图7为本发明实施例的海洋数据层间多次波压制前后放大图。为了进一步比较，放大方框内的数据可得层间多次波压制前后的结果如图7(a)和图7(b)，可以看出，在此海洋数据的复杂结构下，无需任何地下信息，逆散射级数法可以有效地压制层间多次波而不损伤有效信号，提高了成像的品质。

陆地数据采样间隔0.002s，记录长度6.5s，每炮2880道(240道*12线，道间距50m，炮间距50m，接收线距400m。

图8为本发明实施例的陆地数据层间多次波压制前后炮集图。其中，图8(a)为陆地数据经过预处理的炮集，可以看出中间道一次波与多次波相互叠加干扰，部分位置无法分辨出一次波例如箭头处，经过逆散射级数压制层间多次波后，一次波得以分辨，并且同相轴更清晰如图8(b)。

图9为本发明实施例的陆地数据层间多次波压制前后叠加剖面图。其中，图9(a)为层间多次波压制前的叠加剖面，可见在4.2s的箭头处有一个层间多次波形成的同相轴，此层间多次波由箭头上面的两个强反射层产生，经过逆散射级数法压制后，该层间多次波同相轴得以消除，降低了后续解释的难度。

总之，通过实际海洋和陆地数据的处理验证了逆散射级数法对层间多次波压制的有效性，提高地震资料的分辨率，改善了成像品质，提高了解释的可靠性。

通过模拟数据和实际数据的处理进行验证，证明了在上述改进的层间多次波压制方法步骤中，增加了子波去除步骤和子波补偿步骤，通过引入子波的影响，提高了层间多次波预测的准确度，降低了自适应相减的依赖程度。而且，该方法完全数据驱动，无需人工干预，无需已知地下结构，适用于复杂的地形和地质情况。利用该方法进行层间多次波压制后，可以有效地提高资料分辨率，并最终改善成像品质。

虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的保护范围，仍须以所附的权利要求书所界定的范围为准。