0 引言

气枪是海洋地震勘探中的一种常用震源，其性能稳定、可靠，在一定的工作条件下，产生的子波一致性很好。气枪激发时，将高压空气释放到水中，迅速形成球形气泡，由于气泡内的压力远远大于周围海水的静压，导致气泡迅速扩张，瞬间产生第一个正的波形，即地震子波。在浮出水面之前，气泡在水中不断膨胀、压缩，形成多个气泡脉冲^[1]，每一次气泡胀缩都产生地震波。由于气泡不断胀缩导致气枪子波多次振荡，使地震子波长度变长、波形复杂，降低了地震资料的信噪比和分辨率。因此，压制气泡效应一直是海洋地震勘探中的一项重要研究课题。

为了提高震源能量、抑制气泡效应，人们提出利用气枪阵列替代单枪震源采集海洋地震数据。如：Ziolkowski等^[2]提出气枪阵列远场子波数值模拟理论和方法；倪成洲等^[3]、陈浩林等^[4]、杨怀春等^[5]先后进行了气枪阵列参数优选和子波模拟及应用；李绪宣等^[6]研究了立体气枪阵列。在地震数据处理方面：Wood等^[7]提出基于维纳滤波器的气泡效应压制方法；李高林^[8]、Chen等^[9]利用反演方法压制气泡效应；Davison等^[10]、任婷等^[11]利用深水区直达波提取子波压制气泡效应。上述方法都集中于理论研究层面，缺少实际应用及影响因素、局限性等分析。如，由于在浅水区直达波与一次波混叠在一起，文献[10]、文献[11]的方法无法取得好效果。

目前在实际地震数据处理中，考虑到计算效率、可操作性等问题，通常采用基于远场子波的预测反褶积方法压制气泡效应^[12]。在绝大部分情况下，远场子波是理论子波，与实际采集的气枪子波差别较大，因此基于远场子波的气泡效应压制结果存在明显的气泡残余。为此，本文提出了一种数据驱动的气泡效应压制方法，即利用数据的自相关提取子波，基于提取的子波和不含气泡效应的期望子波求取气泡压制算子，进而压制数据中的气泡效应。同时，深入分析了气泡压制算子长度、白噪系数等影响气泡效应压制效果的因素。

1 基于子波提取的气泡效应压制原理 1.1 基于远场子波压制气泡效应的缺陷

目前生产中最常用的压制气泡效应方法是基于远场子波的预测反褶积方法。通过直接测量才能获得真正的远场子波，但是在实际地震勘探中获得远场子波的采集条件非常苛刻，操作成本很高。因此，在实际生产中很少采用直接测量法获得远场子波，绝大部分远场子波是利用商业软件的数值模拟方法得到的。由于模拟过程中理论推导与计算过程中的近似条件的误差、对气泡振荡过程中的阻尼机制研究不充分、数值算法本身的缺陷，使模拟的远场子波与真正的气枪子波差别较大，因此基于远场子波的气泡效应压制结果存在明显的气泡残余。

1.2 从地震数据中提取最小相位子波

地震子波x(t)的振幅谱为

$ A(\omega ) = \left| {\int_{ - \infty }^\infty x (t){{\rm{e}}^{ - {\rm{j}}\omega t}}{\rm{d}}t} \right| $

(1)

式中：t为时间；ω为角频率。根据信号处理原理，有

$ A(\omega ) = \frac{1}{{2\pi }}\left| {\int_{ - \infty }^\infty {{R_{xx}}} (\tau ){{\rm{e}}^{ - {\rm{j}}\omega \tau }}{\rm{d}}\tau } \right| $

(2)

式中R_xx(τ)为地震子波的自相关，τ为时移量。地震数据处理假设地下介质的反射系数是随机的^[13]，基于此基本假设，有

$ {R_{xx}}(\tau ) = {R_{dd}}(\tau ) $

(3)

式中R_dd(τ)为地震数据的自相关。根据式(1)~式(3)，由R_dd(τ)可以求得A(ω)^[14]。

一般认为海上气枪震源激发的子波是最小相位子波，其相位为^[15]

$ \varphi (\omega ) = - {\mathop{\rm HT}\nolimits} [A(\omega )] = - \ln A(\omega ) \cdot \frac{1}{{\pi \omega }} $

(4)

则地震子波为

$ x(t) = \frac{1}{{2\pi }}\int_{ - \infty }^\infty {\left[ {A(\omega ){{\rm{e}}^{{\rm{j}}\varphi (\omega )}}} \right]} {{\rm{e}}^{{\rm{j}}\omega t}}{\rm{d}}\omega $

(5)

式中HT(·)为Hilbert变换。

采用F工区的实际拖缆数据进行试验。图 1为归一化远场子波及从数据中提取的气枪子波。由图可见，从数据中提取的气枪子波(图 1b)和远场子波(图 1a)整体形态类似，但是两者的振幅、气泡的起始时间及起始振幅差别很大。这是由于远场子波是理论计算子波(图 1a)，海面反射系数、激发与接收条件等都呈理想状态，没有考虑实际施工条件对子波的影响。图 2为归一化远场子波及从数据中提取的气枪子波的频谱。由图可见：由于远场子波采用的海面反射系数为-1.0，因此远场子波(图 2a)的检波点(红色箭头)与震源(蓝色箭头)虚反射的陷波点效应均明显强于提取子波(图 2b)，与实际数据不一致；提取子波的气泡效应(图 2b红色圆圈处)更明显，类似于实际数据。

1.3 基于提取的子波压制气泡效应

从地震数据中提取气枪子波x(t)识别气泡效应，利用

$ x\left( t \right){\rm{ }}H\left( t \right) = y\left( t \right) $

(6)

得到不含气泡效应的期望子波y(t)^[16]。式中H(t)为与气泡起始时间有关的窗函数。

以x(t)为输入、y(t)为期望，则有

$ x\left( t \right) * f\left( t \right) = y\left( t \right) $

(7)

利用最小二乘准则，从式(7)中求得滤波算子f(t)(气泡压制算子)，则压制气泡效应后的地震记录为

$ {d^\prime }(t) = d(t) * f(t) $

(8)

式中d(t)为压制气泡效应前的地震记录。

基于提取子波的气泡效应压制方法可分区域或者分道集处理，以应对气枪漏气、天气变化等采集条件剧烈变化时引起的气枪子波变化问题，从而改善气泡效应压制效果。

1.4 应用因素分析

上述基于提取子波的气泡效应压制方法的实际应用效果受很多因素影响，主要有气泡压制算子长度和白噪系数。

1.4.1 气泡压制算子长度

气泡压制算子长度由数据中的气泡延续时间决定^[17]。如果气泡压制算子长度过短，则气泡效应压制结果存在明显的气泡残余。

图 3为A区、B区海上地震数据。由图可见：A区气泡延续时间非常长，在双程反射时间1300ms附近还存在明显的气泡能量(图 3a)；B区气泡延续时间很短，在双程反射时间大于350ms后几乎看不到气泡能量(图 3b)。因此应采用不同的气泡压制算子长度压制A区、B区数据的气泡效应。

1.4.2 白噪系数

利用式(7)求解气泡压制算子f(t)时，为使求解过程稳定，需要预白化处理。当加入白噪后，使压制气泡效应后数据的能量降低，白噪系数越大，能量降低程度越大。但是由于海洋地震数据的特殊性，其存在固有虚反射(鬼波)，选择不同的白噪系数会影响虚反射的陷波效应，进而影响虚反射压制效果。

海洋地震数据的虚反射会在地震子波和数据上引起明显的陷波效应^[18]，若虚反射的延迟时为T，则虚反射算子g(t)在频率域可以表示为

$ G(\omega ) = 1 + r{{\rm{e}}^{ - {\rm{j}}\omega T}} $

(9)

式中r为海面的反射系数。当r=-1.0时，G(ω)振幅谱为

$ Z(\omega ) = 2|\sin (\omega T)| $

(10)

在实际地震数据采集中，由于海浪的存在，海面不可能为镜面反射，此时海面反射系数为接近-1.0的常数。图 4为虚反射算子的振幅谱。由图可见，在0、150Hz附近存在明显的陷波效应。

考虑虚反射的陷波效应，x(t)和y(t)可分别表示为

$ {x(t) = w(t) * g(t)} $

(11)

$ {y(t) = {w^\prime }(t) * g(t)} $

(12)

式中w(t)和w′(t)分别为压制气泡效应前、后的子波，两者的区别集中在低频部分，则式(7)变为

$ w(t) * g(t) * f(t) = {w^\prime }(t) * g(t) $

(13)

将式(13)变换到频率域

$ {W(\omega )G(\omega )F(\omega ) = {W^\prime }(\omega )G(\omega )} $

(14)

则

$ {G(\omega )F(\omega ) = \frac{{{W^\prime }(\omega )}}{{W(\omega )}}G(\omega )} $

(15)

设$G'\left( \omega \right) = \frac{{W'\left( \omega \right)}}{{W\left( \omega \right)}}G\left( \omega \right)$，则式(15)变为

$ {G(\omega )F(\omega ) = {G^\prime }(\omega )} $

(16)

因此

$ {F(\omega ) = \frac{{{G^\prime }(\omega )}}{{G(\omega )}} = \frac{{{G^\prime }(\omega )\overline {G(\omega )} }}{{|G(\omega ){|^2}}}} $

(17)

根据最小平方反褶积的预白化处理，通过式(16)求取滤波算子F(ω)时，需要加入白噪系数为λ的白噪^[19-20]，即

$ {F^\prime }(\omega ) = \frac{{{G^\prime }(\omega )\overline {G(\omega )} }}{{|G(\omega ){|^2} + \lambda {r_{gg}}(0)}} $

(18)

式中r_gg(0)是延迟时为零的g(t)自相关。则加入白噪前、后的滤波算子的关系为

$ \begin{array}{*{20}{l}} {\frac{{{F^\prime }(\omega )}}{{F(\omega )}} = \frac{{{G^\prime }(\omega )\overline {G(\omega )} }}{{|G(\omega ){|^2} + \lambda {r_{gg}}(0)}}\frac{{|G(\omega ){|^2}}}{{{G^\prime }(\omega )\overline {G(\omega )} }}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{{|G(\omega ){|^2}}}{{|G(\omega ){|^2} + \lambda {r_{gg}}(0)}}} \end{array} $

(19)

对于原始地震数据S(ω)，分别利用未加白噪的滤波算子(式(17))和加入白噪的滤波算子(式(18))压制气泡效应，分别得到数据D(ω)和D′(ω)

$ {D(\omega ) = F(\omega )S(\omega )} $

(20)

$ {{D^\prime }(\omega ) = {F^\prime }(\omega )S(\omega )} $

(21)

则由式(19)~式(21)得

$ {D^\prime }(\omega ) = \frac{{|G(\omega ){|^2}}}{{|G(\omega ){|^2} + \lambda {r_{gg}}(0)}}D(\omega ) $

(22)

考虑虚反射的陷波效应，则

$ {D(\omega ) = {S_0}(\omega )G(\omega )} $

(23)

$ {{D^\prime }(\omega ) = S_0^\prime (\omega )G(\omega )} $

(24)

式中S₀ (ω)和S′₀(ω)分别为S(ω)压制气泡效应前、后未考虑虚反射效应的分量。

设加入白噪后的虚反射算子G″(ω)满足

$ {{S_0}(\omega ){G^{\prime \prime }}(\omega ) = S_0^\prime (\omega )G(\omega )} $

则

$ {{D^\prime }(\omega ) = {S_0}(\omega ){G^{\prime \prime }}(\omega )} $

(25)

将式(25)、式(23)代入式(22)，得

$ {G^{\prime \prime }}(\omega ) = \frac{{|G(\omega ){|^2}}}{{|G(\omega ){|^2} + \lambda {r_{gg}}(0)}}G(\omega ) $

(26)

由信号分析原理，得

$ {r_{gg}}(0) = \int_{ - \infty }^\infty | G(\omega ){|^2}{\rm{d}}\omega $

(27)

当r=-1.0时，有

$ {r_{gg}}(0) = 2\pi $

(28)

将式(28)代入式(26)，得

$ {G^{\prime \prime }}(\omega ) = \frac{{|G(\omega ){|^2}}}{{|G(\omega ){|^2} + 2\lambda \pi }}G(\omega ) $

则G″(ω)的振幅谱为

$ U(\omega ) = \frac{{2|\sin (\pi \omega \tau ){|^3}}}{{|\sin (\pi \omega \tau ){|^2} + 2\lambda \pi }} $

(29)

利用理论数据说明白噪系数对虚反射振幅的影响。图 5为海面反射系数分别为0.80、-0.96时，加入不同白噪的虚反射算子的振幅谱及其局部放大。由图可见：①当加入不同的白噪后，虚反射算子振幅谱的形态基本一致，陷波点对应的频率也没有改变，但是能量减小，并且随着白噪系数增加，陷波点处的振幅谱数值更小。上述现象说明，按照实际的拖缆深度和海面反射系数压制虚反射时，无法弥补这种人为增强的陷波效应，因此难以取得较好的虚反射压制效果甚至无效。②随着海面反射系数增大(图 5a~图 5b)，人为增强的陷波效应的幅度增大。因此在实际压制气泡效应时，选择白噪系数为0.01%~1.00%为宜。

2 实际资料应用

图 6展示了利用不同方法对拖缆数据A的压制气泡效应的效果。由图可见：基于远场子波的气泡效应压制方法能够去除部分气泡效应，但是在叠加剖面中还存在明显的剩余气泡能量(图 6b箭头、圆圈和方框处)；当气泡压制算子长度明显小于气泡延续时间时，所提方法的气泡效应压制结果中仍存在气泡残余(图 6c箭头、圆圈和方框处)；当气泡压制算子长度不小于气泡延续时间时，所提方法的气泡效应压制效果更好(图 6d箭头、圆圈和方框处)。图 7为气泡效应压制前、后叠加数据的频谱及其低频部分放大。由图可见：基于远场子波的气泡效应压制结果存在很明显的气泡残余(图 7b红色箭头处)，并且损害了高频信号(图 7a红色方框处)；基于提取子波的气泡效应压制方法更好地压制了气泡效应(图 7b红色箭头处)，几乎不损伤其余频带的信号(图 7a红色方框处)，且当气泡压制算子长度不小于气泡延续时间时，气泡残余更少，频谱更平滑(图 7b黑色圆圈中的绿线)。

图 8为基于提取子波的方法压制拖缆数据B气泡效应前、后叠加剖面。由图可见，压制气泡效应后，明显削弱了低频、虚假同相轴(图 8b)。图 9为气泡效应压制前、后的叠加数据频谱。由图可见：①压制气泡效应后，消除了低频异常振幅，振幅谱更光滑、自然，证明了所提方法的有效性；②当白噪系数分别为0.001%和0.100%时，在压制气泡能量的同时，振幅谱变化很小，尤其是陷波点处的能量几乎无变化；③随着白噪系数增加，能量减小，陷波点处能量减小的趋势更明显，即人为增强了陷波效应，这与理论数据的试验结果(图 5)类似。因此，为了不人为增加陷波效应，在实际压制气泡效应时宜采用尽可能小的白噪系数。

3 结论

本文提出了一种基于提取子波压制气泡效应的方法，该法从数据中提取气枪子波，然后根据提取的子波和不含气泡效应的期望子波求取气泡压制算子压制地震数据气泡效应。这种方法是一种数据驱动的气泡效应压制方法，克服了基于远场子波压制气泡效应方法中的模拟远场子波没有考虑实际采集条件对子波影响的问题，能更好地压制气泡效应，提高数据信噪比和分辨率，同时能更好地保持数据保真度。

文中还分析了影响基于提取子波压制气泡效应方法效果的两个因素：①气泡压制算子长度，影响气泡压制效果，由数据中的气泡延续时间决定；②用于求取气泡压制算子的白噪系数，影响虚反射陷波效应，且随着白噪系数增加，会人为增强虚反射的陷波效应，影响虚反射压制效果。因此，求取气泡压制算子时宜采用尽可能小的白噪系数。上述认识对其他基于子波的气泡效应压制方法具有借鉴意义。