海洋拖缆地震数据采集时, 检波器接收到的记录包括来自地下反射界面的一次反射波和一次反射波经过海水表面反射到达检波器的虚反射(鬼波)^[1]。鬼波紧随一次波之后, 与一次波相互干涉, 降低了地震剖面的分辨率。同时, 鬼波的陷频特性使得地震资料频带变窄, 给地震资料反演和解释带来了困难^[2-3]。

为了获得宽频带的地震数据, 人们提出了各种各类的宽频采集方式^[4]。上下缆采集技术^[5-7]通过对上下缆数据进行合并, 达到去鬼波的目的。双检采集技术^[8]通过对压力和速度传感器采集的记录进行波场分离, 得到上行波记录。变深度缆采集技术^[9-10]利用常规偏移和镜像偏移联合反褶积方法来压制鬼波。但是, 这些宽频采集技术相对于常规水平拖缆采集技术而言, 采集成本较高, 并且需要特定的去鬼波方法。因此, 基于常规水平拖缆采集技术的鬼波压制方法仍然是研究的热点。JOVANOVICH等^[11]利用确定性反褶积压制鬼波, 其压制效果很大程度上取决于鬼波延迟时间估计的精确性, 但精确的鬼波延迟时间往往难以获取, 因此该方法并不能有效压制鬼波。WANG等^[12-13]基于bootstrap方法, 分别在f-xy域和τ-p域实现了去鬼波, 该方法需要通过反演得到最优鬼波延迟时间, 但为了保证反演的稳定性, 需要进行加窗处理, 增加了计算的复杂性。POOLE^[14]结合平面波传播理论, 对最小平方线性拉东方程进行改进, 在重建基准面的同时压制鬼波, 该方法利用最小平方反演直接分离出上行波场, 在提高计算效率的同时提升了鬼波的压制效果。在国内, 王芳芳等^[15]通过构造频率-波数域鬼波压制算子, 实现了数据驱动的鬼波压制; 管西竹等^[16]采用高效率的Fourier变换波场解析延拓方法, 有效解决了水平上下缆鬼波压制中有效信号畸变的问题; 李洪建等^[17]基于格林函数理论, 通过对压力波场和垂直速度波场的预测压制检波点鬼波; ZHANG等^[18]通过考虑实际粗糙海面反射系数对鬼波压制的影响, 利用高斯统计公式精确估计反射系数, 从而得到最优化鬼波压制算子。以上这些去鬼波算法均基于平坦海水表面的假设, 而实际的海水表面是起伏变化的, 忽略起伏海水表面的影响使得鬼波延迟时间和海水表面反射系数的估计存在误差, 导致鬼波压制不完全。

近年来, ROBERTSSON等^[19]利用压力波场信息和近似的压力波场梯度对海水表面进行成像, 实现起伏海水表面去鬼波, 但是该方法对检波器沉放深度和空间采样间隔有严格的限制。KRAGH等^[20]从压力检波器中提取低频信息实现起伏海水表面的成像, 成像质量很大程度上取决于低频信息提取的准确性。ORJI等^[21]将双检数据分解为上行波场和下行波场并分别延拓到海水表面, 再对海水表面进行成像, 该方法受压力和速度传感器各自不同的灵敏度影响较大^[22]。

本文基于POOLE^[14]提出的最小平方线性拉东反演去鬼波方法, 反演求解线性拉东方程得到海水表面处上行波波场, 将海水表面处上行波场正向和反向延拓至检波器处, 利用最小平方成像条件对海水表面进行成像, 以校正检波器实际沉放深度; 利用校正后的深度进行去鬼波处理, 并与POOLE提出的去鬼波方法进行对比。合成数据和实际海上水平拖缆采集数据处理结果表明, 基于起伏海水表面的鬼波压制方法能更好地压制鬼波, 达到拓宽地震记录频带的目的。

1 方法原理 1.1 海水表面处上行波波场的求取

根据频率域最小平方线性拉东变换原理^[23], 检波器接收到的总波场道集D(x, ω)的频率拉东反变换可表示为:

(1)

式中:V(p, ω)为D(x, ω)拉东正变换的结果; ω为角频率; p为慢度参数; x为偏移距。

对于每一个频率分量ω, 可将公式(1) 表示为矩阵形式:

(2)

式中:L=e^{-jωτ_{n, m}}; τ_{n, m}=p_m·x_n; x_n为炮集中第n道的偏移距; p_m为拉东域第m个慢度道所对应的慢度参数。模型V(p, ω)可通过最小平方反演求解线性方程(2) 得到。

对于海洋水平拖缆采集, 二维平面波传播如图 1所示, 第n道z_n处的检波器接收到总波场, 可以看作海水表面x_n+Δx_{n, m}处上行波场回传到检波器的波场与海水表面x_n-Δx_{n, m}处上行波场经海水表面反射到达检波器的下行波场之和, 可将(2) 式中延迟时间τ_{n, m}分解为τ_u和τ_g两项^[14], 其表达式如下:

(3)

(3) 式中Δx_{n, m}, Δτ_{n, m}的几何关系如下:

(4)

式中:v_w为海水速度; θ_m为出射角; z_n为第n道检波器的沉放深度。

结合海水表面处上行波场和检波器处总波场之间的关系, 通过对频率域最小平方线性拉东方程进行改造, 得到新的线性方程, 表达式如下:

(5)

(6)

式中:R为海水表面反射系数; L_ud为拉东反变换算子; L_u=e^-jωτ_u为上行波拉东反变换算子; L_d=e^-jωτ_g为鬼波拉东反变换算子; V₀(p_m)为海水表面处上行波波场拉东变换的结果。

POOLE先通过求解线性拉东方程(5), 得到海水表面处上行波波场V₀(p_m), 然后通过上行波拉东反变换算子L_u将海水表面处上行波波场V₀(p_m)延拓, 得到检波器处上行波波场, 即鬼波压制后的结果。但是, 该方法假设海水表面是平坦的, 忽略了起伏海水表面的影响, 使得海水表面反射系数和检波器沉放深度的估计存在误差。为了消除海水表面反射系数估计的误差, 本文假设海水表面反射系数与地震波频率f以及波场入射角θ有关, 即R(θ, f)=R₀exp(-2f²·σ²cos²θ/v_w²)。其中, R₀为平坦海水表面的反射系数(近似等于-1), σ为海水表面的均方根振幅(近似等于H/2.83, H为海水表面波高)^[11]。起伏海面对检波器实际沉放深度估计的影响通过起伏海水表面成像来消除, 从而实现最大化压制鬼波的目的。

1.2 起伏海水表面成像

对于水平拖缆采集数据, 通过求解公式(5) 得到拉东域海水表面上行波波场V₀(p, ω), 利用频率域拉东反变换得到时间-空间域海水表面上行波波场U₀(x, t)。

(7)

(8)

利用公式(7) 将海水表面处上行波波场反向向下延拓至沉放深度为z_R的检波器处, 得到检波器处的上行波波场; 利用公式(8) 将海水表面处上行波波场经海水表面反射正向向下延拓至沉放深度为z_R的检波器处, 得到检波器处的下行波场。将公式(7) 和公式(8) 延拓得到的检波器处上行波场和下行波场相加, 构造出一个新的检波器处总波场:

(9)

在海水表面平坦的情况下, 检波器理论沉放深度为z_R, 构造的检波器总波场理论上与检波器处接收到的实际总波场D(x, ω)之间差异最小。然而实际采集中由于海水表面起伏变化的影响, 导致检波器实际沉放深度发生变化, 只有在延拓深度等于检波器实际沉放深度的情况下, 两者差异最小。

因此, 假定实际起伏海水表面有效波高为z_sea, 选择向下延拓的深度z′_R=z_R+Δz, 其中Δz为海水表面起伏高度误差, Δz∈[-z_sea, z_sea], 每隔一定微小的深度间隔, 利用公式(7) 和公式(8) 将海水表面上行波场U₀(x, t)延拓至深度z′_R处; 对每个深度构造的总波场与实际接收的总波场利用最小平方成像条件^[24]

(10)

将每道最大成像值I(x, z′_R)对应的Δz作为该道海水表面起伏高度, 实现对起伏海水表面的成像，进而得到每道检波器实际沉放深度z′_R=z_R+Δz; 用z′_R校正公式(4) 中检波器理论沉放深度z_n, 再次利用公式(5)进行去鬼波处理, 得到校正后的海水表面上行波场; 最后, 利用上行波拉东反变换算子L_u将校正后的海水表面上行波场延拓, 得到校正后检波器处上行波波场, 即最终鬼波压制的结果。

在实际资料去鬼波处理中, 海水表面的起伏形态随时间变化, 可通过选定时窗, 使得时窗内海水表面影响一致或者相似, 准确估计海水表面反射系数和检波器实际沉放深度, 然后对时窗内的数据进行去鬼波处理。

2 数值模拟算例

为了验证基于起伏海水表面的鬼波压制方法的鬼波压制效果, 对层状模型(图 2)的水平拖缆单炮记录进行去鬼波处理。正演模拟时, 采用主频为80Hz的雷克子波, 记录道数为61, 道间距为4m, 采样间隔为0.5ms, 记录长度为0.3s, 拖缆沉放深度为15m。为了验证本文方法的抗噪性, 加入了信噪比为15的随机噪声。

模拟平坦海水表面和起伏海水表面下的单炮记录分别如图 3a和图 3b所示, 图 3c为图 3a和图 3b的差剖面。可见, 由于受到起伏海水表面的影响, 海水表面反射的鬼波同相轴不再是标准的双曲线, 呈现出与起伏海水表面波动吻合的弯曲(图 3中红色方框所示)。对比图 3a和图 3b中各道的旅行时差发现起伏海面引起的鬼波同相轴最大旅行时差为1.5ms, 与所设计的起伏海水表面完全吻合。相反, 一次反射波经过海底反射直接传播到检波器处, 并不受起伏海水表面的影响, 其同相轴符合标准的双曲线。

图 4a, 图 4b为利用POOLE提出的方法对平坦海水表面和起伏海水表面下模拟的单炮记录进行鬼波压制的结果。可见, POOLE提出的方法能完全压制平坦海水表面下模拟单炮记录中的鬼波; 但对于起伏海水表面, 由于该方法忽略了起伏海水表面的影响, 鬼波压制结果中存在残余的鬼波(图 4b中蓝色箭头所指处), 同时会引入较强的窜扰噪声(图 4b中红色箭头所指处)。

利用本文方法对起伏海水表面条件下模拟的单炮记录(图 3b)的海水表面形态进行成像, 结果如图 5中蓝色曲线所示, 黑色曲线表示平坦海水表面, 红色曲线表示模拟的起伏海水表面。可见, 成像的海水表面与模拟的起伏海水表面在形态上匹配良好, 部分差异主要是由主频80Hz子波的分辨率限制所引起的。

基于成像的海水表面对起伏海水表面下模拟的单炮记录进行鬼波压制, 得到的去鬼波结果如图 4c所示。与POOLE方法的去鬼波结果(图 4b)相比, 由于考虑了起伏海水表面对鬼波的影响, 残余鬼波得到了压制, 窜扰噪声也随之减少, 压制效果得到明显提升。但去鬼波结果中仍存在少量的窜扰噪声, 这是由海水表面成像误差所导致的。

对去鬼波前后记录进行频谱分析, 如图 6所示, POOLE方法(蓝色曲线)和本文方法(红色曲线)去鬼波后的振幅谱中陷波现象得到消除, 能量得到补偿, 频带得到拓宽。本文方法对于低频端能量提升更加明显(绿色方框所示)。

3 实际水平拖缆数据测试

利用海上某工区水平拖缆采集资料进一步验证本文方法的适用性(图 7), 其采样间隔2ms, 道间距12.5m, 拖缆沉放深度为25m, 采样总时间为10.572s, 震源沉放深度为15m。图 7a为经过直达波切除、滤波去噪、多次波压制等预处理后的原始单炮记录(选取2.6~3.4s, 10~190道)。首先, 利用本文方法对实际单炮记录的海水表面形态进行成像, 结果如图 8所示。然后, 基于成像的海水表面进行鬼波压制, 得到鬼波压制后的单炮记录, 如图 7c所示。与由POOLE方法压制鬼波后的单炮记录(图 7b)进行对比可见, 2种方法压制鬼波后的单炮记录中反射同相轴相位个数减少, 炮记录变得更加清晰, 分辨率得到提高。从两种方法差剖面(图 7d)中可以看出, 本文方法压制鬼波的记录中残余鬼波更少, 压制效果更佳。

从成像的海水表面形态上看, 30~80道和120~170道受海水表面起伏影响较大, 两种方法在相应区域的鬼波压制效果存在微小的差异(由图 7d差剖面可以明显看出)。为了验证本文方法的优越性, 选取去鬼波前后相应的区域进行频谱分析, 如图 9a和图 9b所示。在去鬼波前的炮记录中, 在30, 50, 60, 90Hz处存在明显的陷波现象, 经过鬼波压制, 2种方法都能有效消除检波器鬼波所引起的陷频特征, 拓宽地震记录的频带。其中50Hz为震源鬼波引起的陷波点, 因此, 能量变化不大。与POOLE方法相比, 本文方法无论是对低频能量还是高频能量的提升都具有明显优势(图 9黑色箭头所指处)。

4 结论

本文通过对起伏海水表面进行成像, 实现了基于起伏海水表面的拖缆鬼波压制。通过对模拟单炮数据和实际海上数据处理结果的分析, 可得以下结论和认识:

1) 基于起伏海水表面的拖缆鬼波压制方法能有效地压制拖缆采集数据中的鬼波, 恢复陷波点能量, 提高了地震资料的分辨率, 拓宽了地震资料的频带。

2) 相对于基于平坦海水表面的鬼波压制方法, 基于起伏海水表面的拖缆鬼波压制方法考虑了起伏海水表面对鬼波的影响, 鬼波压制效果更佳, 低频和高频能量的提升更加明显, 能够适应实际起伏海水表面情况。

3) 起伏海水表面的拖缆鬼波压制方法基于最小平方成像条件实现了对海水表面的成像, 但起伏海水表面的成像精度在一定程度上取决于地震数据的信噪比。