频率—慢度域延拓关键参数优化的鬼波压制方法

青杰,顾汉明*,王建花

(1.中国地质大学(武汉)地球物理与空间信息学院，湖北武汉430074； 2.中海油研究总院有限责任公司，北京100028；3.海洋油气勘探国家工程研究中心，北京100028)

0 引言

海上勘探中，拖缆通常沉放于水面以下，受水面强反射的影响，海洋地震资料中地层反射波后会伴随反相位的鬼波。鬼波的陷波效应使得地震记录有效频带变窄，故压制鬼波是拓宽地震资料频带的重要途径。

为了得到分辨率和信噪比更高的宽频地震资料，上/下缆采集[1-3]、双检接收技术[4-5]、上/下源激发技术[6-7]等相继提出。此外，斜缆采集也被推出并应用，其通过变化检波器深度的采集方式使陷波特性多样化，压制鬼波后可较好地拓展地震资料的频带。Soubaras等[8-9]提出原始数据结合镜像数据联合反褶积压制鬼波；Sablon等[10]将多道联合反褶积技术应用于同步多级震源与二、三维斜缆数据；Wang等[11]提出了倾斜拖缆“bootstrap”处理方法，用于处理τ-p域中的斜缆数据；Song等[12]在高分辨率Radon变换基础上进行倾斜拖缆数据处理，拓宽了地震频带；张振波等[13]在珠江口盆地使用斜缆采集、处理方法获得了高品质的海上地震数据；王建花等[14]通过波场延拓生成镜像数据共轭迭代反演获得去鬼波地震数据；张威等[15]利用最小平方残差法反演求解海表面上行波，并延拓得到拖缆处的上行波波场；马继涛等[16]将基于波场外推和阈值截断的方法应用于鬼波压制；宋建国等[17]提出利用格林函数多次迭代压制鬼波的方法；张宝金等[18]研究了倍频视角下的宽频地震子波优化方法，用于鬼波压制后宽频地震数据子波整形；王兆旗等[19]对不同时刻的静态粗糙海面组合得到了动态海面，用以测试鬼波压制效果。在复杂采集环境中，假设条件下的鬼波参数与实际参数相比误差较大，为此，Grion等[20]提出了以相移法延拓波场来压制鬼波的算法，对澳大利亚近海起伏海面上的数据有着良好的应用效果；Matta等[21]提出了自适应去鬼波方法；许自强等[22]研究了抑制鬼波的优化联合反褶积算法；王冲等[23-24]分别在频率—波数域与频率—慢度域对斜缆进行迭代反演压制鬼波，扩展了地震频带；张威等[25]在平缆数据上对反演得到的海面波场正、反向延拓到拖缆处获得拖缆深度，消除了起伏海面影响；封强等[26]以负熵度量去鬼波结果的非高斯性来搜索得到更精确的海面反射系数和鬼波延迟时间，继而通过反褶积压制鬼波；王永强等[27]基于编码/解码思想，在Bayes反演框架下进行了鬼波预测与压制；陈胜红等[28]以自举法τ-p域稀疏反演鬼波压制等多种方法为核心，形成了宽频准三维处理集成技术。

鬼波滤波方法需考虑鬼波与一次波的延迟时间和振幅差异系数，但崎岖海底、起伏海面及目的层深度等因素会使鬼波参数受到一定的扰动，经常规方法处理后记录中会存在噪声。本文基于频率—慢度域鬼波压制方法构建了波场延拓算子，用原始记录分别加上、减去延拓一次的记录，以两个新记录乘积的绝对值之和作为度量参数，交替迭代寻找最优平均鬼波延迟时间和振幅差异系数，进而去除鬼波。另外，为降低鬼波参数时变的影响，利用滑动时间窗口将地震记录离散化，使得鬼波参数在窗口内近似相等，再分别对窗口内数据寻找最优鬼波参数以消除鬼波。斜缆合成数据和实际资料的试算结果表明，该方法能较好地消除鬼波，达到拓宽地震记录频带的目的。

1 频率—慢度域鬼波压制原理

在拖缆沉放至某一深度时，由于海水面起伏、海底崎岖以及复杂构造等因素，拖缆鬼波与一次波的振幅差异系数不再为-1，这里用R表示该差异系数。对于同一水平慢度pj的平面波，由图1中鬼波波场ug与一次波波场up的几何关系可得

图1 检波器接收记录示意图

(1)

式中：τ为波前面的截距时间；Δτ为鬼波与一次波的延迟时间；xi、zi分别为检波器炮检距离和深度，其中i=1,2,…,N，且N是记录的道数；θj为水平慢度pj对应平面波的出射角度,j=1,2,…,M，其中M为慢度个数；VW是海水波场传播速度。

检波器的波场u是一次波与鬼波相加的结果，即

u(pj,τ)=up(pj,τ)+uG(pj,τ)

=up(pj,τ)+Rup(pj,τ-Δτ)

(2)

该接收点总波场d(xj,t)为M个水平慢度之和，即

(3)

式中τ=t-pjxi，其中t为时间。对式(3)做关于时间t的傅里叶变换，得

(4)

式中ω为频率。对于N道数据而言，每一个频率ω有

D(x,ω)=G(x,p,ω)U(p,ω)

(5)

式中：D(x,ω)为有N个元素的列向量；U(p,ω)是有M个元素的列向量，即

D(x,ω)=[D(x1,ω),D(x2,ω),…,D(xN,ω)]T

(6)

U(p,ω)=[U(p1,ω),U(p2,ω),…,U(pM,ω)]T

(7)

G(x,p,ω)为N×M阶矩阵，即

G(xi,pj,ω)=(1+Re-iω Δ τ)e-iωpjxi

(8)

U=(GTG+μI)-1GTD

(9)

式中：μ为阻尼系数；I为单位矩阵。

2 频率—慢度域延拓原理

将式(8)中G(xi,pj,ω)分解可得

G(xi,pj,ω)=Gp(xi,pj,ω)-Gg(xi,pj,ω)

(10)

式中：Gp(xi,pj,ω)=e-iωpjxi；Gg(xi,pj,ω)=-Re-iω Δ τe-iωpjxi。由式(5)可得

D(x,ω)=[Gp(x,p,ω)-Gg(x,p,ω)]U(p,ω)

(11)

(12)

(13)

(14)

以此延拓类推，则有

(15)

3 鬼波参数扰动分析

由式(8)可见，鬼波算子与鬼波参数(振幅差异系数及延迟时间)有关。实际资料中，鬼波参数会因采集环境不同而有所变化，海底的复杂地形会使鬼波与一次波的路径相差较大，延迟时间因此不同；此外，由于检波器中接收的鬼波和上行波不是从同一个点反射，两者的振幅差异系数也会偏离假定的反射系数。海面是随时间起伏的，拖缆深度随之产生变化，鬼波延迟时间也随空间和时间变化；反射层深度会使得延迟时间时变，也会影响鬼波和上行波的振幅差异系数。当然影响鬼波延迟时间的因素还有海水中波场传播速度等。

在此以延迟时间的扰动简述崎岖海底、起伏海面和反射地层深度等三个环境因素对鬼波参数的影响。采用两层模型简单模拟各个因素引起的延迟时间扰动，基础模型如图2a所示。海底与海面相距500 m，选取的检波器位于炮检距为1000 m、深度为20 m的位置，海水中波场传播速度固定为1500 m/s，采用50 Hz的Ricker子波，各个模型参数因考察影响因素不同而改变。图2b是海底情况，海底反射情况复杂，此处仅以鬼波的海底反射点位置抬升简单演示。图2c对应于海面起伏情况，检波器的深度因此变化。图2d为目的层深度对延迟时间的影响，为操作方便对该情况进行简化，以海底深度变化查看目的层深度引起的延迟时间扰动。

图3a～图3c为图2b～图2d模型对应的数值模拟记录，为同一个检波器在影响因素改变时接收到的记录，用于观察延迟时间的变化，此处已将一次波校平，每个记录的第一道对应影响因素未改变的情况；图3d～图3f为图3a～图3c减去各自第一道记录的结果，这就使得各个因素对延迟时间的影响清晰可见。随着鬼波在海底的反射点位置抬升，延迟时间随之减小(图3a)，当上抬高度到达1.0 m后差值记录(图3d)能量就不能忽视；由图3b可见，当海面上升，鬼波与一次波的延迟时间随之增加，当海面上升0.5 m后，鬼波就与改变前有较大的相位差(图3h)；随着目的层深度的增加，对应的延迟时间也随之变化(图3c、图3f)，且前期变化快，后期渐缓，深度相差较大时，延迟时间差也对应增大。深度越大，反射时间越长，可见不同反射时间的延迟时间会有些许差距，但当目的层越深，鬼波与一次波的路径就越接近，反射地层深度造成的延迟时间扰动就越小。

图2 延迟时间影响因素

图3 海底模型及数值模拟记录

上述只是一个简单的模拟，实际可能的情况要复杂得多。如海面是随时间起伏的，则延迟时间也对应增大或减小，是时变的，在整个拖缆上海面的影响也是变化的，故而该因素对应的延迟时间扰动也是空变的；海底情形多变，并不仅是上抬，还有如倾斜、下陷、凸起等，影响复杂；反射地层深度造成的扰动变化快、慢难以统一度量，各因素交错作用。另外，炮检距也会对各个因素造成的延迟时间扰动产生影响。

以上情况说明，各因素作用会让鬼波参数产生随时间和空间变化的不容忽视的扰动，这使得鬼波参数在实际采集中难以估算。因此，在处理海洋地震记录时，需要相对地震记录求最优平均延迟时间和振幅差异系数，获得一个整体最优的处理结果。

4 参数寻优原理

4.1 寻优原理

设Gg+(xj,pj,ω)=-Rie-iω Δ τie-iωpjxi,GG-(xj,pj,ω)=-e-iω Δ τie-iωpjxi/Ri，则

(16)

(17)

(18)

(19)

然后，对Q+进行反傅里叶变换得时间—空间域记录q+，对Q-进行反傅里叶变换得时间—空间域记录q-。将寻优判定函数设置为q+与q-乘积绝对值之和，即

(20)

式中S(xi)为第i道的判定函数；q+相当于原记录加上延拓后记录；q-相当于原记录减去延拓后记录。当鬼波与一次波的延迟时间和振幅差异系数最优时，延拓后记录的一次波移动到原记录中鬼波的位置，q+中的原记录鬼波与延拓后记录的一次波相加抵消，得到如式(17)右边所表达的记录结果，即一次波和延拓后的鬼波；当延迟时间或振幅差异系数有误差时，原记录的鬼波与延拓后记录的一次波将不能完全抵消，会有残余噪声。q-与q+相对应，在q+中振幅因原记录鬼波与延拓记录一次波相加削弱的时间段，q-则因原记录鬼波与延拓后记录的一次波相减增强，鬼波参数最优时便如式(18)右边所示。

当q+与q-相乘并取绝对值相加时，原记录鬼波所在时间段的振幅变化是影响相乘结果的主要原因。若延迟时间和振幅差异系数最优，q+中原记录鬼波去除彻底，q-中增强后的原记录鬼波与q+中极小的原记录鬼波残留相乘并取绝对值相加后结果就较小；相对的若延迟时间和振幅差异系数有误差，则q+中原记录鬼波所在时间段鬼波残留较多，与q-中增强的鬼波进行相乘并取绝对值相加时结果就会偏大。由此，通过两个新构建记录的波形特征可对不同延迟时间和不同振幅差异系数的记录处理结果进行判定，找到最优处理结果对应的延迟时间和振幅差异系数。每道最小S对应的延迟时间和振幅差异系数即为所求结果。近炮检距记录可适当增加延拓次数以改善效果。

4.2 滑动时窗

从本文的正演地震记录(图8a)提取出一道记录，即图4a展示的第100道记录。图4a标注了该道6个主要由一次波和鬼波构成的反射波组，分别提取了这6个反射波组的鬼波延迟时间(图4b)，可以看到鬼波的延迟时间是随时间波动的，前期变化较快，后期平缓，可知鬼波参数是时变的。

图4 第100道记录点位图(a)及对应的延迟时间示意图(b)

图5为滑动时窗示意图，滑动时窗中每一个时窗长度对应一个最优延迟时间和振幅差异系数。分时窗时应该尽量将鬼波和上行波分在同一个时窗，前、后滑动窗口应该有重叠部分，以保证地震数据的全覆盖。分时窗时窗口不宜太大，太大则前、后鬼波参数相差较大，导致无法求得一个合适的延迟时间和振幅差异系数；窗口也不宜太小，太小则计算量增加，时间方向滑动窗口约为200～500 ms即可，在该时段采集环境对鬼波参数的影响相对一致，空间方向长度根据处理情况调整。开始计算时，可先选择一炮记录进行时窗长度试算，在尽量使鬼波与一次波在同一个窗口的前提下，以一定的间隔增加时窗长度，将每一个时窗长度在炮记录上进行处理测试，分别查看应用效果，选择其中效果较好的时窗长度中的高值，将之扩散应用到整个接收记录。在记录长度的中、后段可以适当增加时窗长度，减小计算量。具体的窗口可根据实际处理情况实时修改。

图5 滑动时窗示意图

4.3 鬼波压制步骤

(1)划分滑动时间窗口，保证覆盖整个地震记录；

(3)在所给范围内枚举延迟时间，对于每一个延迟时间Δτ，振幅差异系数不变，按式(16)～式(20)计算对应的S值；

(4)对每个延迟时间的S比较大小，最小值所对应延迟时间即为最优延迟时间；

(5)在振幅差异系数范围内枚举鬼波振幅差异系数，保持步骤(4)所得到的延迟时间不变，对于每一个振幅差异系数，由式(16)～式(20)计算对应的S值；

(6)比较每个振幅差异系数的S值，最小值所对应的振幅差异系数即为所求最优鬼波振幅差异系数；

(7)用所得的振幅差异系数，重复步骤(3)和步骤(4)，得到新的延迟时间，再由新延迟时间得到新的鬼波振幅差异系数，据此不断迭代，直至到达最大迭代次数；

(8)以所得到的延迟时间和振幅差异系数，由式(8)和式(9)计算去鬼波记录U(p,ω)，再进行频率域Radon反变换和傅里叶反变换，得到时间—空间域的去鬼波记录U(x,t)；

(9)对每个窗口内的数据重复步骤(2)～步骤(8)，最后将所有窗口数据整合成一个完整的记录。

利用延拓寻找最优平均鬼波延迟时间和振幅差异系数进行鬼波压制方法的流程图如图6所示。

图6 鬼波压制流程

5 合成数据试算

为验证本文方法去除鬼波的效果，利用图7a所示复杂盐丘模型的正演模拟记录进行试算。斜缆采用图7b抛物线式斜缆，炮点位于800 m处，检波器沉放深度为6～50 m，记录共250道，道间距为10 m，时间采样间隔为0.002 s，记录长度为2.9 s，采用主频为35 Hz的Ricker子波。

图7 复杂盐丘模型(a)与斜缆深度随炮检距变化曲线(b)

在原始炮集(图8a)中，因为海面的强反射，一次波后紧跟着一个反相的鬼波，鬼波与一次波的延迟时间随着炮检距的增大呈先增后减的状态。图8b是拖缆鬼波压制后的炮集。由图可见，经本文方法处理后，一次波后的反相同相轴基本消失，整个记录的相位变得相对单一。对参数寻优时需将数据转换到频率—慢度域进行延拓，再换算到时间—空间域进行处理结果判定。将原始炮集转换到时间—慢度域时(图8c)，可发现与图8a相似的特征，鬼波与一次波对应的同相轴先逐渐分开再靠拢，两者相位表现相反，各同相轴与时间—空间域的同相轴相对应。图8d是去鬼波炮集的时间—慢度域记录，其中与鬼波对应的同相轴得到有效压制，鬼波去除相对干净。

图8 不同域去鬼波前、后的合成记录

从图9a的近炮检距记录可以看到鬼波与一次波随着炮检距的增大逐渐分开的过程，经鬼波压制后(图9b)，鬼波同相轴的能量得到极大程度的衰减。可见本文方法对鬼波有着良好的去除效果。鬼波会与一次波相互干涉，使得频谱出现陷波点，陷波点处频谱能量会有较大衰弱(图10a)。在原始炮集的频率—空间域振幅谱(图10a)中可以看到陷波频率随检波器深度变化，显现出了斜缆陷波特征的多样性。在去除鬼波后，陷波点处的能量得到有效地恢复，频谱变得连续，拓宽了记录的有效带宽(图10b)，验证了本文方法可压制鬼波、弥补由鬼波引起的陷波点能量损失。

图9 时间—空间域去鬼波前(a)、后(b)合成记录局部图

图10 去鬼波前(a)、后(b)合成记录振幅谱

6 实际数据试算

为进一步验证本文方法的有效性，对实际海上斜缆采集数据进行试算，该数据震源深度为5 m，道间距为12.5 m，时间采样间隔为0.002 s，总计为300道，斜缆大致呈抛物线状，拖缆沉放深度为5～50 m。图11a是经过直达波切除、去噪等预处理的斜缆炮集记录。在实际资料中，因各种采集环境因素造成的扰动，鬼波与一次波的实际延迟时间与理论上有一定误差，鬼波与一次波的振幅差异系数也不再与假定的相同。

海上斜缆数据的反射波主要出现在2 s后(图11a)，鬼波紧随在上行波后，两者间的距离随着炮检距的增大呈先逐渐分开后渐渐靠近的变化，相位表现相反。经过鬼波参数的延拓寻优，扰动的鬼波参数得到修正，经过鬼波处理后，紧随上行波的反相同相轴得到极大压制(图11b)，记录更加清晰，有效信号更加突出，说明本文方法可较好地消除斜缆的鬼波。实际资料变换到时间—慢度域，同样出现鬼波与一次波成对出现的特征(图11c)，两者相位相反，鬼波与一次波同相轴随着慢度增加先逐渐分离再慢慢并拢，应用本文方法进行鬼波处理后，时间—慢度域中鬼波对应的同相轴基本消失(图11d)，表明鬼波得到有效去除。从图12的局部对比图也能看到一次波后的同相轴能量得到较大衰减，炮记录相位更加单一，鬼波得到较好地压制，进一步说明了基于延拓寻找最优参数的鬼波压制方法可较好地去除由海面虚反射引起的鬼波，提高资料的分辨率。

图12 图11去鬼波前(a)、后(b)实际资料局部放大显示

从图13a可以看到实际资料中鬼波引起的陷波特征以及斜缆采集数据中陷波频率多样性的特点。经过鬼波处理后(图13b)，陷波处的能量得到有效补偿，振幅谱连续性增加，使得地震记录的频带得到拓宽，较好地去除了鬼波带来的陷频影响，说明本文方法在实际数据处理中的有效性。在鬼波去除前后第100道至第150道记录振幅谱对比图中(图14)，可见到多个陷波频率以及陷波频率处明显的能量低值，经本文方法去除鬼波后，振幅谱变得更加连续，低频和高频的能量都有所增益，斜缆陷波频率多样性的特点得到较好利用，提高了资料的分辨率，表明本文方法对资料的宽频处理有着良好的作用。

图13 去鬼波前(a)、后(b)实际资料振幅谱

图14 鬼波压制前、后第100～150道振幅谱对比

7 结论

通过对斜缆合成数据以及海上实际斜缆资料的处理分析，可以得到如下结论：

(1)崎岖海底、起伏海面、地层反射构造及拖缆沉放深度误差等因素直接影响到鬼波相对于一次波的振幅及延迟时间的扰动，继而影响到传统方法鬼波压制效果；

(2)结合滑动时窗方法，基于频率—慢度域延拓算子，通过交替迭代求解，得到最优鬼波相对于一次波的振幅差异系数和延迟时间；

(3)基于频率—慢度域延拓的鬼波参数最优化拖缆鬼波压制方法应用于合成和实际海上拖缆地震炮集记录的鬼波压制，鬼波参数扰动的影响得到消除，鬼波压制效果得以提升，地震资料的频带明显拓宽。