在海洋和陆地地震勘探中, 由于地下强反射界面的存在, 地震波经多次反射形成层间多次波, 并与一次波叠加干扰, 严重降低了地震资料的分辨率, 加大了有效波识别的难度, 影响了地震成像品质和地震解释成果的可靠性, 因此, 衰减或去除层间多次波是地震资料处理的重要环节。为去除层间多次波的干扰, 提高地震资料分辨率, WEGLEIN^[1]提出了两类多次波压制方法:一类是基于一次波与多次波之间特性差异的滤波法; 另一类是基于波动理论的预测相减法。滤波法包括预测反褶积法、Radon变换法、f-k滤波法和聚束滤波法等, 此类方法在满足假设条件时能有效衰减或去除多次波, 且计算量小, 计算效率高, 但需要假设地下信息。当一次波和多次波之间的特征差异很小或不存在时, 很难获得理想的压制效果, 甚至会严重损伤一次波。预测相减法避免了滤波法需要假设地下信息的局限, 无需先验信息, 代表了多次波压制的发展趋势, 主要包括反馈迭代法和逆散射级数法。反馈迭代法需人工干预, 利用逐层指定多次波产生算子预测层间多次波。逆散射级数法^[2]完全由数据驱动, 无需人工干预, 可一次预测出所有的层间多次波, 是目前最先进的层间多次波压制方法。MATSON等^[3]将逆散射级数法应用于海洋地震勘探数据处理; FU等^[4]将该方法应用于陆上地震勘探资料处理并取得了较好的效果; LIANG等^[5]分析了逆散射级数法预测时出现的虚假层间多次波问题, 并利用高阶项加以去除; ZOU等^[6]将逆散射级数法的应用从衰减层间多次波拓展到了完全去除层间多次波; WU等^[7]将逆散射级数法应用于含衰减因子Q的地震资料的层间多次波压制, 验证了方法的有效性; YANG等^[8]分析了震源子波和辐射方向模式对层间多次波预测的影响, 改进了逆散射级数法预测层间多次波的精度; JIN等^[9]提出了改进的1.5维逆散射级数法预测层间多次波。

本文从散射理论出发, 首先详细推导了逆散射级数层间多次波压制方法, 然后改进了常规层间多次波压制处理流程, 在逆散射级数法预测层间多次波前、后去除和补偿子波来提高层间多次波预测精度, 从而降低对自适应相减的依赖度, 最后将该方法应用于模型数据和实际数据的处理, 结果证明了方法的有效性和适用性。

1 方法原理 1.1 逆散射级数预测层间多次波的物理机制

散射理论是关于散射波场和扰动介质关系的基本理论。在散射理论^[10]中, 实际介质包括参考介质(或称背景介质)和扰动介质, 相应地, 实际波场也包括参考波场(背景波场)和扰动波场(或称散射波场)。在参考介质中传播的波场称为参考波场, 在实际介质中传播的波场称为实际波场, 扰动波场为实际波场和参考波场的差^[11]。实际波场和背景波场分别满足以下波动方程:

$ LG = \delta $

(1)

$ {L_0}{G_0} = \delta $

(2)

式中:L, L₀分别为实际介质和背景介质的微分算子; G, G₀分别为实际介质和背景介质的波场; δ为狄拉克算子。我们定义扰动算子V≡L₀－L, 散射波场ψ_S≡G－G₀, 散射方程(Lippmann-Schwinger)可以表示为:

$ G = {G_0} + {G_0}VG $

(3)

将(3)式中G代入到(3)式右侧, 得到正演散射级数:

$ \begin{array}{l} G = {G_0} + {G_0}V{G_0} + {G_0}V{G_0}V{G_0} + {G_0}V{G_0}V{G_0}V{G_0} + \\ \;\;\;\;\;\; \cdots \end{array} $

(4)

已知扰动算子V(与介质信息有关)和背景波场G₀, 散射波场为:

$ \begin{array}{l} {\psi _S} = {G_0}V{G_0} + {G_0}V{G_0}V{G_0} + {G_0}V{G_0}V{G_0}V{G_0} + \cdots \\ \;\;\;\;\; = {\left( {{\psi _S}} \right)_1} + {\left( {{\psi _S}} \right)_2} + {\left( {{\psi _S}} \right)_3} + \cdots \end{array} $

(5)

式中:(ψ_S)_n为V的n阶方程式, 是散射波场ψ_S的一部分。(4)式和(5)式适用于任意位置, 在地表测得的散射波场称为地震数据D=(ψ_S)_ms。根据D的阶数展开扰动算子V, 可得:

$ V = {V_1} + {V_2} + {V_3} + \cdots $

(6)

式中:V_n为数据D的n阶方程式。将(6)式代入到(5)式中, 可得以下一阶、二阶和三阶逆散射级数:

$ D = {G_0}{V_1}{G_0} $

(7)

$ {G_0}{V_2}{G_0} = - {G_0}{V_1}{G_0}{V_1}{G_0} $

(8)

$ \begin{array}{l} {G_0}{V_3}{G_0} = - {G_0}{V_1}{G_0}{V_2}{G_0} - {G_0}{V_2}{G_0}{V_1}{G_0} - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{G_0}{V_1}{G_0}{V_1}{G_0}{V_1}{G_0} \end{array} $

(9)

利用(7)式、(8)式和(9)式可得V₁, V₂, V₃等, 并最终获得扰动算子V和实际介质信息。但逆散射级数的总体收敛性并未验证, 故可选取部分子级数进行针对性处理, 如去除表面多次波子级数、去除层间多次波子级数、成像子级数和反演子级数等。

本文选取去除层间多次波子级数, 从正演的角度分析点散射模型反射波和层间多次波产生的物理机制。层间多次波为在自由表面以下的任何界面发生多次反射的波, 发生下行反射的次数为层间多次波的阶数。如图 1所示, 一阶层间多次波发生了1次下行反射, 二阶层间多次波发生了两次下行反射。此外, 下行反射波的散射点深度必须小于相邻两个散射点的深度, 且层间多次波的形成与散射点的个数有关。由正演模拟可知, 至少要有3个散射点才能形成层间多次波, 故一阶和二阶散射级数无法形成层间多次波。同理, 反演时一阶和二阶的逆散射级数也不能形成层间多次波。将公式(9)简化如下:

$ {G_0}{V_3}{G_0} = - {G_0}{V_{31}}{G_0} - {G_0}{V_{32}}{G_0} - {G_0}{V_{33}}{G_0} $

(10)

式中:V₃₁=V₁G₀V₂; V₃₂=V₂G₀V₁; V₃₃=V₁G₀·V₁G₀V₁。由(10)式可知, 层间多次波从三阶逆散射级数中的第3项才开始产生(至少包含3个散射点)。

图 2给出了三阶逆散射级数第3项中散射点的不同位置组合, z₁, z₂, z₃分别为从左到右3个散射点的深度, 与正演模拟类似, 只有满足z₁＞z₂且z₃＞z₂和阶数限制条件的逆散射级数才能形成层间多次波。

1.2 逆散射级数预测层间多次波的算法

无限空间中常速度背景介质在两点z和z′之间的格林函数, 可表示为:

$ G\left( {z,z',q} \right) = - \frac{{{{\rm{e}}^{{\rm{i}}q\left| {z - z'} \right|}}}}{{2{\rm{i}}q}} = - \frac{1}{{2{\rm{ \mathsf{ π} }}}}\int_{ - \infty }^\infty {\frac{{{{\rm{e}}^{{\rm{i}}q\left( {z - z'} \right)}}}}{{{q^2} - {{q'}^2}}}{\rm{d}}q'} $

(11)

式中:q为垂向波数。将公式(11)右端积分项分为主值和奇点q′=±q处的函数值, 即:

$ \begin{array}{*{20}{c}} {\frac{1}{{{q^2} - {{q'}^2}}} = PV\left( {\frac{1}{{{q^2} - {{q'}^2}}}} \right) + \frac{1}{{2{\rm{ \mathsf{ π} }}}}\left\{ {{\rm{i}}\pi \frac{1}{{2\left| q \right|}}\left[ {\delta \left( {q' - } \right.} \right.} \right.}\\ {\left. {\left. {\left. q \right) + \delta \left( {q' + q} \right)} \right]{{\rm{e}}^{{\rm{i}}q'\left( {z - z'} \right)}}} \right\}} \end{array} $

(12)

只考虑格林函数G₀在奇点处的贡献, 将公式(12)右边第2项代入到逆散射级数的三阶公式(10)中右边第3项可得(省略变量ω) :

$ \begin{matrix} {{V}_{33}}\left( {{k}_{\text{g}}},-{{q}_{\text{g}}},{{k}_{\text{s}}},{{q}_{\text{s}}} \right)=\int{\int_{-\infty }^{+\infty }{\cdot \left[ {{V}_{1}}\left( {{k}_{\text{g}}},-{{q}_{\text{g}}},k,q \right)\cdot \right.}} \\ \left. {{V}_{1}}\left( k,q,{k}',-{q}' \right){{V}_{1}}\left( {k}',-{q}',{{k}_{\text{s}}},{{q}_{\text{s}}} \right) \right]/\left[ \left( 2\text{i}q \right)\cdot \right. \\ \left. \left( 2\text{i}{q}' \right) \right]\text{d}k\text{d}{k}'+\int{\int_{-\infty }^{+\infty }{\left[ {{V}_{1}}\left( {{k}_{\text{g}}},-{{q}_{\text{g}}},k,q \right)\cdot \right.}} \\ \left. {{V}_{1}}\left( k,q,{k}',{q}' \right){{V}_{1}}\left( {k}',{q}',{{k}_{\text{s}}},{{q}_{\text{s}}} \right) \right]/\left[ \left( 2\text{i}q \right)\left( 2\text{i}{q}' \right) \right]\cdot \\ \text{d}k\text{d}{k}'+\int{\int_{-\infty }^{+\infty }{\left[ {{V}_{1}}\left( {{k}_{\text{g}}},-{{q}_{\text{g}}},k,-q \right){{V}_{1}}\left( k,-q, \right. \right.}} \\ \left. \left. {k}',-{q}' \right){{V}_{1}}\left( {k}',-{q}',{{k}_{\text{s}}},{{q}_{\text{s}}} \right) \right]/\left[ \left( 2\text{i}q \right)\left( 2\text{i}{q}' \right) \right]\cdot \\ \text{d}k\text{d}{k}'+\int{\int_{-\infty }^{+\infty }{\left[ {{V}_{1}}\left( {{k}_{\text{g}}},-{{q}_{\text{g}}},k,-q \right){{V}_{1}}\left( k,-q, \right. \right.}} \\ \left. \left. {k}',-{q}' \right){{V}_{1}}\left( {k}',{q}',{{k}_{\text{s}}},{{q}_{\text{s}}} \right) \right]/\left[ \left( 2\text{i}q \right)\left( 2\text{i}{q}' \right) \right]\text{d}k\text{d}{k}' \\ \end{matrix} $

(13)

式中:k_s, k_g分别表示震源和检波器的水平方向波数; q_s, q_g分别表示相应的震源和检波器的垂向波数。此处只考虑奇点处的贡献, 而不考虑主值的贡献, 其物理意义在于奇点处的贡献在逆散射级数中起着异常速度扰动的作用, 与速度模型有关, 主值的贡献在逆散射级数中只与背景速度模型有关。

单频平面波场b₁为:

$ \begin{matrix} {{b}_{1}}\left( {{k}_{\text{g}}},{{k}_{\text{s}}},{{q}_{\text{g}}}+{{q}_{\text{s}}} \right)=-2\text{i}{{q}_{\text{s}}}D\left( {{k}_{\text{g}}},{{k}_{\text{s}}},\omega \right) \\ =\frac{{{V}_{1}}\left( {{k}_{\text{g}}},-{{q}_{\text{g}}},{{k}_{\text{s}}},{{q}_{\text{s}}} \right)}{-2\text{i}{{q}_{g}}} \\ \end{matrix} $

(14)

根据公式(13)和(14), 引入z₁＞z₂且z₃＞z₂约束条件, 可构造得出层间多次波场的预测算法^[10]:

$ \begin{array}{l} {b_3}\left( {{k_{\rm{g}}},{k_{\rm{s}}},{q_{\rm{g}}} + {q_{\rm{s}}}} \right) = \frac{1}{{{{\left( {2{\rm{ \mathsf{ π} }}} \right)}^2}}}\int_{ - \infty }^\infty {{\rm{d}}{k_1}} \int_{ - \infty }^\infty {{\rm{d}}{k_2}{{\rm{e}}^{ - {\rm{i}}{q_1}\left( {{\varepsilon _{\rm{g}}} - {\varepsilon _{\rm{s}}}} \right)}} \cdot } \\ \;\;\;{{\rm{e}}^{{\rm{i}}{q_2}\left( {{\varepsilon _{\rm{g}}} - {\varepsilon _{\rm{s}}}} \right)}} \times \int_{ - \infty }^\infty {{\rm{d}}{z_1}{b_1}\left( {{k_{\rm{g}}},{k_1},{z_1}} \right){{\rm{e}}^{{\rm{i}}\left( {{q_\text{g}} + {q_{\rm{l}}}} \right){z_1}}}} \cdot \\ \;\;\;\int_{ - \infty }^{{z_1} - \varepsilon } {{\rm{d}}{z_2}{b_1}\left( {{k_{\rm{1}}},{k_2},{z_2}} \right){{\rm{e}}^{ - {\rm{i}}\left( {{q_1} + {q_{\rm{2}}}} \right){z_2}}}} \times \int_{{z_2} + \varepsilon }^\infty {{\rm{d}}{z_3}} \cdot \\ \;\;\;{b_1}\left( {{k_2},{k_s},{z_3}} \right){{\rm{e}}^{{\rm{i}}\left( {{q_2} + {q_s}} \right){z_3}}} \end{array} $

(15)

式中:ε_s, ε_g分别为震源和检波器的深度; q_s, q_g分别为震源和检波器的垂向波数; k₁, k₂分别为震源和检波器的水平波数积分变量; k_s, k_g分别为震源和检波器的水平波数, 满足$ {q_i} = {\rm{sgn}}\left( \omega \right)\sqrt {{{(\omega /{c_0})}^2} - k_i^2} $, 其中, i为震源s或检波器g, c₀为背景介质速度; ε为使z₁＞z₂和z₃＞z₂(“低-高-低”)严格成立的参数, 在实际数据处理中ε与子波波长相关。如果3个一次波满足“低-高-低”约束关系将产生一阶层间多次波, 如果一阶层间多次波与两个一次波满足“低-高-低”约束关系就形成了二阶层间多次波, 同理高阶层间多次波也可由低阶层间多次波合成。输入数据中含有的信号包括一次波和各阶层间多次波, 因此利用该算法可一次性预测全部层间多次波。

利用逆散射级数法((15)式)预测层间多次波时至少需对3个子波进行褶积处理, 预测的层间多次波波形和振幅发生了很大变化, 使得同相轴变“胖”, 能量差异变大, 因此子波的存在影响了层间多次波预测的准确性。本文改进了层间多次波处理流程, 在逆散射级数层间多次波的预测步骤前、后分别增加了子波去除和子波补偿的步骤以提高预测的准确性。改进后的单频平面波场为:

$ {b_1}\left( {{k_{\rm{g}}},{k_{\rm{s}}},{q_{\rm{g}}} + {q_{\rm{s}}}} \right) = - 2{\rm{i}}{q_{\rm{s}}}A{\left( \omega \right)^{ - 1}}D\left( {{k_{\rm{g}}},{k_{\rm{s}}},\omega } \right) $

(16)

式中:A(ω)为子波。由震源波场的远场近似或利用格林理论^[10]可求得海洋地震资料的远场子波, 进而衰减子波; 利用盲反褶积可衰减陆地地震资料的子波。由于难以准确提取子波, 故不易完全去除子波的影响, 但是近似子波也很大地提高了多次波预测的准确性。多次波预测后, 经子波补偿后的层间多次波M为:

$ M\left( {{k_{\rm{g}}},{k_{\rm{s}}},\omega } \right) = {\left( { - 2{\rm{i}}{q_{\rm{s}}}} \right)^{ - 1}}A\left( \omega \right){b_3}\left( {{k_{\rm{g}}},{k_{\rm{s}}},{q_{\rm{g}}} + {q_{\rm{s}}}} \right) $

(17)

逆散射级数层间多次波压制方法的优点如下:①完全由数据驱动, 无需人工干预; ②无需已知地下介质信息和速度模型, 适用于各种复杂的地形和地质情况; ③可一次性预测所有层间多次波。

在实际地震资料处理时, 首先进行预处理, 海洋地震资料的预处理包括低切滤波、去气泡、压制涌浪噪声、去除直达波、压制鬼波和压制自由表面多次波等, 陆地地震资料的预处理包括静校正、异常振幅相干噪声压制、地表一致性处理、剩余振幅补偿和剩余静校正等; 然后对预处理后的地震资料进行去除子波、构建全波场数据、傅里叶变换和常速度偏移等处理; 最后采用逆散射级数层间多次波压制方法预测并压制层间多次波。图 3给出了常规和改进的逆散射级数层间多次波压制处理流程。

2 实例分析 2.1 模拟数据应用

利用两个简单层状模型数据验证和分析逆散射级数层间多次波压制方法。模型数据1包括两个反射层, 第1层速度为1 500 m/s, 第2层速度为2 000 m/s, 震源为主频25 Hz的雷克子波, 采样间隔0.002 s, 记录长度为2.0 s。图 4a为模型数据1, 其中前2个信号为一次波, 随后的4个信号为一阶、二阶、三阶和四阶层间多次波。图 4b和图 4c分别为采用常规逆散射级数法预测的层间多次波和压制多次波后的结果。可以看出, 常规逆散射级数法预测了层间多次波准确的双程旅行时, 但是振幅和波形与实际差异明显, 自适应相减后存在部分多次波残留。图 4d为去除子波后的模拟数据, 图 4e为采用改进的逆散射级数法预测的层间多次波, 图 4f为采用改进的逆散射级数法压制层间多次波后的结果。该结果表明应用改进的逆散射级数法可准确预测层间多次波的双程旅行时, 且波形和振幅与模型数据图 4a中的层间多次波基本一致, 因此自适应相减法可更易去除层间多次波。图 4g, 图 4h和图 4i分别为模拟数据1及应用常规和改进的逆散射级数法压制层间多次波后结果的灰度显示。

模拟数据2同样采用主频25 Hz的雷克子波, 共2 001道, 道间距5 m, 记录长度为1.5 s。图 5a为模拟数据2, 包含两个一次波和一阶、二阶层间多次波。图 5b和图 5c分别为利用改进和常规逆散射级数法预测的层间多次波。可以看出, 改进的逆散射级数法预测的层间多次波与模拟数据2中的层间多次波能量相近, 而常规逆散射级数法预测的层间多次波波形变“胖”且与模拟数据2中的层间多次波的能量差异较大。为了进一步对比波形和振幅, 从图 5a、图 5b和图 5c中选取近偏移距(图 5d, 图 5e和图 5f)和远偏移距(图 5g, 图 5h和图 5i)数据进行比较, 图 5d为模型数据中近偏移距的层间多次波, 图 5e和图 5f分别为采用改进和常规的逆散射级数法预测的层间多次波的近偏移距数据; 相应地, 图 5g为模型数据中远偏移距的层间多次波, 图 5h和图 5i分别为采用改进和常规的逆散射级数法预测层间多次波的远偏移距数据。图 6a和图 6b分别为模拟数据2中的层间多次波和应用改进的逆散射级数法预测的层间多次波在近偏移距和远偏移距的对比, 可以看出, 该方法预测的层间多次波与模拟数据中的层间多次波波形和能量基本一致。两个模型数据测试均证明改进的逆散射级数法提高了层间多次波预测的准确性。

2.2 实际数据应用

利用逆散射级数层间多次波压制方法分别对1套深海地震数据和1套陆地地震数据进行多次波压制处理。

深海地震数据采集参数为:震源和检波器深度分别为8.5 m和9.5 m, 采样间隔0.004 s, 记录长度9 s, 每炮239道, 共1 751炮, 炮间距50 m, 道间距25 m, 最小偏移距160 m。图 7a为预处理和压制表面多次波后得到的地震数据, 由于层间多次波算法的计算量大, 故只选取具有岩底结构的部分数据(图 7a蓝线内)进行层间多次波压制, 图 7b为压制层间多次波后结果。为便于比较, 放大图 7蓝线内压制层间多次波前、后的结果分别如图 8a和图 8b所示, 可以看出, 逆散射级数法可以有效地压制层间多次波而无损有效信号, 提高了成像品质。

陆地地震数据采样间隔0.002 s, 记录长度6.5 s, 每炮2 880道(240道×12线), 道间距50 m, 炮间距50 m, 接收线距400 m。图 9a为预处理后的炮集, 可以看出中间道附近一次波与多次波相互叠加干扰, 部分位置(如箭头处)处无法分辨出一次波, 采用改进的逆散射级数法压制层间多次波后, 被掩盖的一次波得到恢复, 且同相轴更清晰, 如图 9b所示。图 10a为压制层间多次波前的叠加剖面, 箭头处有1个层间多次波形成的同相轴, 此层间多次波由箭头上方的两个强反射层产生。采用改进的逆散射级数法压制层多次波后, 去除了该多次波的同相轴, 降低了后续解释的难度(图 10b)。

上述深海和陆地数据的处理结果验证了逆散射级数法对层间多次波压制的有效性, 提高了地震资料的分辨率和成像品质。

3 结论

本文从散射理论出发, 提出了一种改进的逆散射级数层间多次波压制方法, 并给出了改进的逆散射级数法层间多次波处理流程。该方法通过在逆散射级数法预测层间多次波前、后去除和补偿子波, 提高了层间多次波预测的准确度, 降低了对自适应相减的依赖程度。该方法具有以下优点:①完全数据驱动; ②无需人工干预; ③无需已知地下信息和速度模型; ④适用于复杂的地形和地质情况。模型数据测试证明该方法能有效提高层间多次波预测的准确性并提高压制效果。实际数据处理结果表明, 该方法可以有效压制地震资料的层间多次波, 提高资料分辨率, 改善成像品质, 适用于复杂的海上和陆地地震数据处理, 尤其对深层和超深层地质目标的勘探具有广阔的应用前景。

致谢: 衷心感谢休斯敦大学WEGLEIN教授在本文论证期间给予的指导和帮助。