1 引言

拓宽地震勘探信号的频带，不仅能够增强地震成像效果，提高分辨率，而且能够提高波阻抗反演和全波形反演的精度(ten Kroode et al., 2013；Soubaras et al., 2013).而海上拖缆地震采集的记录受海水面的虚反射(也称鬼波)影响，存在陷波特性，使得地震记录的频带变窄，不仅降低了地震剖面的分辨率，而且缺失了地震反演所需要的低频分量(Kragh et al., 2010；李套山等，1997).因此，近几年来，国内外开展了海水宽频带地震采集新方法试验(吴志强等，2013；吴志强，2014; 余本善和孙乃达，2015)，旨在充分利用地震采集数据的特点和有针对性的数据处理手段压制鬼波.这些试验包括海上双传感器地震拖缆采集(Carlson et al., 2007)、上下拖缆采集技术(Hill et al., 2006；刘春成等，2013；赵仁永等，2011)，以及变深度拖缆(VDS)或倾斜缆采集技术(Soubaras and Dowle, 2010；Sablon et al., 2011)，其中倾斜缆采集技术因其充分利用同一条拖缆中陷波频率的多样性，通过在F-XY域、T-X域及τ-p域运用专门的宽频处理算法有效压制鬼波(Soubaras，2010；Soubaras，2012；Lin et al., 2011；Wang and Peng, 2012; Wang et al., 2013；Song et al., 2015)，能较好地拓宽地震记录频带，因此，该技术在国内得到了研究和推广应用，并取得了较好的应用效果(谢玉洪等，2012；张振波和许自强，2014；许自强等，2015；王冲等，2016).现有的技术在推导鬼波滤波算子时仅假定海水面反射系数为-1，没有考虑到崎岖海底和目的层深度使得鬼波和一次反射波的振幅差异系数不再是一个恒定值而是随偏移距的变化而变化，进而影响了压制鬼波的效果，尤其在远偏移距时此问题更严重，为此，本文推导出频率慢度域中鬼波滤波算子以及最小二乘反演算法求解上行波算法，该鬼波滤波算子与不同水平慢度对应的鬼波延迟时间和海水面反射系数有关.并针对实际缆深误差引起的鬼波延迟时间估计误差以及鬼波与一次反射波振幅差异随偏移距和目的层深度的变化而难以给定一个固定值的问题，基于计算出的下行波与实际下行波之间的平方误差最小理论实现自适应反演迭代计算最优鬼波与一次反射波振幅差异系数和鬼波延迟时间.合成及某海上采集的倾斜缆数据处理结果表明，该方法能较好地压制变深度缆由海水面产生的鬼波，达到了拓宽地震记录频带的目的.

2 频率慢度域鬼波滤波算子

如图 1所示，粗实线显示的倾斜拖缆在A点处的沉放深度为z_i，偏移距为x_i，i=1, …N，N为单炮记录总道数.不妨假设某一出射角为θ_j的上行平面波到达A点处的波前面为AB，该波前面的截距时间为τ，波场为u^U(p_j, τ)，p_j为水平方向(沿着x方向)的慢度，且有

(1)

该平面波在到达海水面C点处，海水面反射下行达到A点处的下行波为u^D(p_j, τ)，在A点处接收到的总波场u^S(p_j, τ)是上行波场和下行波场之和，也即

(2)

由于下行波场u^D(p_j, τ)是上行波场在海水面反射向海底方向传播的下行波场，对于复杂构造介质尤其崎岖海底，地下同一反射界面上不同反射点反射的振幅是不同的，而在同一接收点上接收到的上行波是深层反射界面上某点处的一次反射波，同一接收点上接收的鬼波是深层反射界面上另外一点处的一次反射波上行至海水面并经海水面反射的下行波，而且这两个上行波在海底出射点也不同，导致在同一接收点上接收到的一次反射波和海水面反射的鬼波除了符号相反外，还存在振幅差异，该差异与构造复杂程度和目的层深度有关，也即鬼波和一次反射波的振幅差异随偏移距的变化而变化，因此，不妨假定鬼波与一次反射波振幅差异用系数R_j=R(p_j)表示，其值大小与出射角有关，根据平面波传播理论以及图 1中的几何关系，可以得到

(3)

(4)

令垂直方向(沿着z方向)的慢度为q_j，则有

(5)

式(5)代入式(4)，则式(3)可以写成

(6)

这里，

(7)

式(6)代入式(2)得

(8)

则对于M个水平慢度p_j(j=1, …M)在接收点A处的总波场d(x_i, t)可以表示为

(9)

基于时间域离散线性τ-p正反变换可知，式(9)中的u^U(p_j, τ)就是A点处的上行波场u(x, t)的τ-p正变换u^U(p, τ)的第p_j分量，且有τ=t-p_jx_i.

对式(9)做关于时间t的傅里叶变换，得

(10)

这里，G(x_i, p_j, ω)为鬼波滤波算子在频率慢度域表示形式，即

(11)

3 自适应最小二乘反演迭代去鬼波算法和步骤

对式(3)做关于时间t的傅里叶变换，得

(12)

对于倾斜缆N道而言，对每一个频率分量ω，式(10)可表示为矩阵形式：

(13)

这里，D (x, ω)为N个元素组成列向量，即

(14)

U (p, ω)为M个元素组成列向量，即

(15)

G (x, p, ω)是由N×M个元素组成的矩阵，其元素可表示为

(16)

在利用式(7)计算延迟时间以及给定海水面反射系数前提下，就可采用最小二乘反演方法求解线性方程组(13)，可得(Yilmaz，1989；Song et al., 2015)

(17)

式中，ε为阻尼系数，I为单位矩阵.

由于拖缆沉放深度通常有误差，导致由式(7)计算的鬼波延迟时间τ_j有误差；其次，鬼波与一次反射波的振幅差异系数与平面波出射角有关，难以用一个恒定的值表示，也即式(16)难以准确计算，为此，本文采用自适应迭代计算最优振幅差异系数R_j和鬼波延迟时间τ_j，然后采用式(17)求解上行波，也即压制鬼波后的波场.

一方面，从式(12)可以看出，对某个慢度p_j传播的上行波(一次反射波)施加一个算子R_je^-jωτ_j便得到了该慢度传播的下行波(鬼波)；另一方面，由于倾斜缆接收到的总波场是上行波场和下行波场之和，因而如果上行波已知的话，就可以从下式计算下行波场，即

(18)

式中，是倾斜缆接收到的总波场频率域Radon正变换.因此，分别基于式(12)和式(18)计算出的下行波应该相等，但是，正如前面所述，鬼波延迟时间τ_j有误差，鬼波与一次反射波的振幅差异系数难以用一个恒定的值表示，为此，这里采用最小二乘反演迭代算法，找到反射波与鬼波之间的最优平均鬼波延迟时间和鬼波与一次反射波的最佳振幅差异系数，也即利用式(12)计算出的下行波与式(18)计算的下行波之间的平方误差最小求得最优平均鬼波延迟时间和最佳振幅差异系数，依次迭代，直至误差小于某一阈值为止.最优平均鬼波延迟时间和最佳振幅差异系数迭代计算方法和步骤如下：

设第k次迭代时，对于某个水平慢度p_j，可以基于式(17)计算出所有频率分量ω对应的上行波U^(k)(p_j, ω)，则下行波可通过下式求得，即

(19)

则第k次迭代时，与某个水平慢度p_j对应的最优平均鬼波延迟时间和最佳振幅差异系数可通过使得式(19)得到的实际下行波U^{D (k)}(p_j, ω)与式(12)计算的理论下行波之间的平方误差最小求得，即

(20)

根据式(20)，采用最小二乘反演迭代，求得第k次迭代时的最优平均鬼波延迟时间τ_j^(k)和最佳振幅差异系数R_j^(k)，再利用式(16)也就求得了第k+1次迭代时的最优鬼波滤波算子矩阵G ^(k)(x, p, ω)中的元素G^(k)(x_i, p_j, ω)，对于所有的水平慢度p_j，重复使用式(19)-(20)，就可以得到矩阵G ^(k)(x, p, ω)中的所有元素，再基于式(17)，可得到第k+1次迭代时的压制鬼波的结果，即上行波向量为U ^(k+1)(x, p, ω).从而反复使用式(18)-(20)不断地反演迭代，直至误差E小于某一阈值为止.

为了最小二乘求解式(20)，令

(21)

化简得

(22)

由于鬼波与一次波之间的延迟时间随着偏移距的增大而减少，因此，最大延迟时间为2倍缆深最大值除以海水层速度，即τ_max=2z_max/V_w，z_max为拖缆最大沉放深度，假定因缆深沉放深度误差引起的延迟时间误差小于半个周期，因此，延迟时间τ_j∈(0, τ_max+T/2)，T为子波周期.而振幅差异系数R_j∈[-1, 0)，因此，可以用牛顿法求解非线性方程式(22)得到第k次迭代时的最优平均鬼波延迟时间τ_j^(k)和最佳振幅差异系数R_j^(k).

频率域τ-p变换自适应最小二乘迭代反演压制鬼波方法的具体步骤可概括如下：

(1)将倾斜缆炮集记录d(x, t)进行1D傅里叶变换得到向量D (x, ω)，再进行频率域Radon正变换，得到；考虑到平面波分解是以倾斜电缆沉放深度处的波场值进行计算，相对于水平电缆接收而言，倾斜电缆接收相当于给地下原本水平界面增加了一个倾角，同时考虑到实际地下界面的起伏，这里Randon变换中的水平视慢度的取值为-1/V_min~1/V_min，V_min为地层最小层速度；

(2)对于所有的频率分量ω、所有的道和水平慢度p_j，先假定振幅差异系数R_j为-1，并用式(7)计算鬼波延迟时间τ_j，利用式(11)计算得到G(x_i, p_j, ω)，进而可以得到矩阵G (x, p, ω)；

(3)对于所有的频率分量ω和水平慢度p_j，利用式(16)求解得到初始的上行波场向量U (p, ω)；

(4)基于最优平均鬼波延迟时间和最佳振幅差异系数迭代计算步骤求得每一个水平慢度对应的最优平均鬼波延迟时间向量和最佳振幅差异系数；

(5)对于所有的频率分量ω和水平慢度p_j，利用式(11)计算得到矩阵G (x, p, ω)，进而基于式(17)求解得到最终的上行波场向量U (p, ω)；

(6)进行频率域Radon反变换，得到；

(7)将进行反傅里叶变换，得到u(x, t)，即得到去除鬼波后的对应倾斜缆的炮集记录.

4 合成倾斜缆炮集记录试算

为了验证变深度缆最小二乘反演迭代去鬼波技术对由海水自由反射界面所引起的鬼波的压制效果，对一个复杂盐丘模型(图 2a)的合成变深度缆记录进行去鬼波试算，正演模拟时，采用主频为30 Hz的雷克子波，记录道数为300道，道间距为10 m，采样间隔为2 ms，记录长度为3.0 s，偏移距为0 m，单边放炮.斜缆采集是不同接收点处的拖缆沉放深度不同，斜缆的形态有多种，通常采用抛物线形态(Amundsen and Landr∅, 2014)，其数学表达式如下：

(23)

这里可以看出，抛物线形态斜缆由拖缆沉放最小深度或起始深度z_min、最大沉放深度z_max以及拖缆倾斜段水平距离x_c等三个参数控制，深度z是沿着Z轴测量的，炮检距x是沿着X轴测量的，第一式提供斜缆沉放深度数据，而第二式提供水平缆沉放深度.起始深度z_min通常为几米到十多米，x_c也即拖缆倾斜段起始点到终止点之间的水平长度，其值从几百米到几千米.如图 2b所示的拖缆三个参数分别为6 m、50 m和1500 m.

图 3a为基于有限差分法波动方程正演模拟的变深度缆单炮记录，炮点位置在400 m，明显看出，在第一道时间在1.1 s处有一组海底界面上的反射波同相轴，在1.5 s、1.8 s、2.1 s、2.5 s等处分别有一组地层界面上的反射波同相轴，可以看出每组反射波同相轴具有两个强相位特征，第一个强相位对应一次反射波，而续至的第二个强相位对应海水面的鬼波，而且随着炮检距的增大，一次波和鬼波逐渐明显分开又逐渐靠近，也可以从图 4a和4c显示的局部放大结果看出这个特征.图 3d为τ-p域炮集记录，也可以看出这些特征，随着慢度的增大，一次波和鬼波逐渐明显分开又逐渐靠近.

基于自适应最小二乘反演迭代去鬼波算法，对原始单炮记录进行去鬼波处理.首先对图 3a显示的t-x域炮集记录d(x, t)对时间t进行傅里叶变换，得到各道的频谱D(x, ω)，其振幅谱如图 5a所示，再进行频率域Radon正变换，得到，其振幅谱如图 5c所示，从图上可以看出，鬼波存在所产生的陷波特征，如图中箭头所示，这些陷波频率随着道号(偏移距)或慢度的变化而明显变化，这也是倾斜缆采集陷波频率多样性特征.通过自适应最小二乘迭代反演压制鬼波的步骤，求解得到最终的上行波场U(p, ω)，其振幅谱如图 5d所示，再进行频率域Radon反变换，得到，其振幅谱如图 5b所示，从图上可以看出，由于鬼波得到压制，陷波能量得到补偿，如图中箭头所示，这也是倾斜缆采集能够得到宽频的原因之一.

通过对进行反傅里叶变换，得到u(x, t)，即得到去除鬼波后的对应倾斜缆的炮集记录，为了比较，这里将去鬼波前后的t-x域和τ-p域倾斜缆单炮记录及其残差显示在一起，图 3b、3e分别为去鬼波后的t-x域和τ-p域倾斜缆单炮记录，图 3c、3f分别为t-x域和τ-p域残差或所去除的鬼波数据，图 4分别为鬼波压制前后的单炮记录近偏移距及远偏移距局部放大显示结果，对比明显可以看出，去鬼波后的单炮记录上的反射波主要为一个正相位和一个负相位特征，远近偏移距道上的鬼波都能得到较好的压制，说明了自适应最小二乘反演迭代去鬼波算法可以较好地消除倾斜缆处对应的海水面的鬼波.

图 6为原始和去鬼波后的t-x域单炮记录的频率波数谱(FK谱)，明显看出，原始炮集记录的FK谱(图 6a)存在明显的陷波特征(箭头所示)，而去鬼波后的炮集记录的FK谱(图 6b)上的陷波得到了消除，能量得到补偿，频谱连续，频带得到了拓宽，更进一步验证了自适应最小二乘反演迭代去鬼波算法可以有效压制大偏移距接收缆鬼波，弥补了由于鬼波存在引起的频谱缺失.

5 实际倾斜缆采集数据试算

实际海上斜缆采集观测系统参数道间距为12.5 m，炮间距为25 m，接收长度6.25 s，采样率为2 ms，震源深度为5 m，斜缆三个参数分别为5 m、50 m和3000 m.由于海况影响，实际采集时电缆沉放深度与设计的理论深度之间存在一定的误差，实际采集时海水面也是波浪起伏，导致鬼波实际到达接收点的时间与理论时间存在一个差值，鬼波与一次波的振幅也不同.

图 7a和图 7d分别显示了涌浪干扰滤波后的t-x域、τ-p域倾斜缆单炮记录，明显可以看出每组反射波同相轴具有两个强相位特征，第一个强相位对应一次反射波，而续至的第二个强相位对应海水面的鬼波，而且随着炮检距的增大，一次波和鬼波逐渐明显分开又逐渐靠近.

基于自适应最小二乘反演迭代去鬼波算法，对原始单炮记录进行去鬼波处理，图 8显示了去鬼波前后的不同域记录的振幅谱，与倾斜缆合成记录振幅谱特征类似，对比图中箭头位置处的振幅谱，可以看出去鬼波前的记录由于鬼波存在产生陷波特征，而在去鬼波后的记录频谱上，陷波能量得到较好的补偿，说明了自适应最小二乘反演迭代去鬼波算法可以较好地消除倾斜缆处对应的海水面的鬼波.

图 7b、7e分别为实际斜缆采集的去鬼波后的t-x域和τ-p域倾斜缆单炮记录，图 7c、7f分别为t-x域和τ-p域残差或所去除的鬼波数据，对比明显可以看出，去鬼波后的单炮记录上的反射波主要为一个正相位和一个负相位特征，远近偏移距道上的鬼波都能得到较好的压制，进一步说明了自适应最小二乘反演迭代去鬼波算法可以较好地消除倾斜缆处对应的海水面的鬼波.

图 9为原始和去鬼波后的t-x域单炮记录的频率波数谱(FK谱)，明显看出，整炮原始记录的FK谱(图 9a)存在明显的陷波特征(椭圆圈所示)，而去鬼波后的整炮记录的FK谱(图 9b)上的陷波得到了消除，能量得到补偿，频谱连续，频带得到了拓宽.从图 9c和图 9d显示的炮集记录中对应海底反射的FK谱更明显看出这些特征，进一步验证了自适应最小二乘反演迭代去鬼波算法可以有效压制大偏移距接收缆鬼波，弥补了由于鬼波存在引起的频谱缺失.

图 10显示出了偏移距为2800 m的原始及去鬼波后的共偏移距剖面，可以明显看出，原始共偏移距剖面上反射波出现了两条同相轴，使得分辨率降低，尤其对深层反射同相轴影响更为严重；而去鬼波结果中压制了其中的一条同相轴，分辨率得到了提高.这也可以从图 11显示出的整条地震记录剖面的一维平均叠加频谱分析结果看出，在原始记录频谱图(细实线)上有多个陷波频率，而在去鬼波压制后的频谱图(粗实线)上，陷波消去了，说明本文提出的去鬼波算法能充分利用变深度缆陷波频率发散特征，有效弥补陷波，扩宽频带.

还可以从压制鬼波前后的共中心点道集记录及其自相关谱(图 12)分析去鬼波的效果，首先采用波场延拓法对变深度缆采集的原始和去鬼波炮集数据进行校正到水平缆，然后进行抽道集和动校正，得到动校正后的共中心点道集记录.对比图 12a和图 12b，从动校正后的原始共中心点道集记录，可以明显看出存在于一次反射波后的鬼波，且存在明显的剩余时差，而动校正后的去鬼波后的道集记录上，只有一次反射波同相轴，鬼波得到较好的压制.图 12c和图 12d分别对应去鬼波前后的CMP道集上黑色粗框窗口内各道的自相关图.对比可以看出，在图 12c上，由于鬼波的存在，使得每一条自相关曲线上都呈现出多个极值，而在图 12d上，由于鬼波得到压制，使得每一条自相关曲线上只有一个极值，这些图件进一步说明了频率慢度域自适应最小二乘反演迭代去鬼波算法在压制变深度缆鬼波的应用效果.

图 13显示了去震源和接收缆鬼波前后的偏移叠加剖面，对比可以看出，原始偏移叠加剖面上海底及目的层由正负峰组成的多条同相轴在去鬼波后的最终处理剖面上变成了单峰组成的同相轴，海水面引起的鬼波得到明显压制，地震记录的分辨率得到明显提高，达到了宽频处理的目的.

6 结论

通过模型变深度缆模拟记录及实际采集的变深度缆去鬼波处理结果分析，可以得出如下结论：

(1)频率慢度域导出的鬼波算子中同一水平慢度对应不同深度反射层的鬼波延迟时间以及鬼波与一次波的振幅差异系数相同，克服了在时空域中鬼波算子中在同一炮检距下难以给定一个固定的鬼波延迟时间以及鬼波振幅差异系数的问题.

(2)在频率慢度域可以通过理论下行波与实际下行波之间的平方误差最小实现自适应反演迭代计算得到鬼波最优延迟时间以及振幅差异系数，进而通过频率慢度域阻尼最小二乘反演算法计算得到压制鬼波后的上行波，合成记录验证了该算法的可靠性.

(3)将频率慢度域自适应最小二乘反演迭代算法应用于某海上变深度缆去鬼波技术后的炮集记录，其振幅谱频带明显比原始记录的振幅谱频带宽，在最终处理的偏移叠加剖面上，中高频谱振幅能量得到提高，使得频带变得更宽，中深层反射的分辨率得到明显提高，信噪比也得到提高，断点更清晰.