海洋地震是目前最为有效的海洋油气勘探技术，然而在海洋地震勘探中，由于海水表面以及海底或海底之下的地质界面均为高波阻抗界面，会产生强多次波，严重影响地震数据的信噪比，为后续地震高精度成像带来困难^[1]，因此近几十年来，多次波的预测与剔除始终是海洋地震勘探的研究热点之一。

近年来，以波动方程为基础的方法在多次波压制方面取得了显著进展，该类方法主要包括波场外推法^[2]、自由界面多次波衰减方法(SRME)^[3-6]和逆散射级数法^[5]等，目前已成为实际应用中的主要多次波压制手段。波动方程类多次波消除方法其多次波压制通常可以分为两个主要步骤：多次波预测和匹配衰减，即首先，通过波动理论对地震数据中的多次波进行预测，随后将预测出的多次波与原始数据进行匹配并衰减，从而去除多次波成分。在此过程中，多次波匹配衰减算法是决定波动方程类多次波消除方法处理效果的关键步骤^[7]。

当前，多次波匹配衰减方法主要分为两类：一类是多道维纳滤波法，另一类方法则是基于数学变换域的自适应匹配衰减。在多道维纳滤波类方法方面，Verschuur等^[5]通过最小平方准则成功压制了预测的自由界面多次波，使SRME方法得以应用于实际地震数据处理；Wang^[8]在此基础上，提出了扩展多道匹配滤波和伪多道匹配滤波，显著提高了多次波的压制效果；Li等^[9-10]开发了均衡拟多道匹配滤波方法，以在一次波与多次波同相轴交叉时提升自适应衰减性能；Li等^[11]进一步优化了伪多道匹配滤波，增强了保幅的能力；另一类方法则是基于数学变换域的自适应匹配衰减，Zhou等^[12]提出在抛物线拉东域中应用的匹配衰减法，通过抛物线拉东变换处理地震数据和多次波记录，再利用蒙版滤波去除原始记录中的多次波成分；Wang^[13]利用拉东域中的模板数据在时域中实现自适应相减；Shi等^[14]通过在拉东域投影构建非线性滤波器进行自适应相减；Herrmann等^[15]提出了一种基于曲波域的新方法来分离一次波与多次波，而Dong等^[16]基于此方法开发了复曲波变换的多次波自适应衰减技术。

相对来讲，多道维纳滤波法在处理各种多次波时表现良好，但当一次波和多次波的同相轴交叉或重叠时，可能会削弱有效波的完整性；而变换域的匹配衰减方法在处理相交的多次波时表现更优，但对多次波预测的精确度要求较高。然而基于波动理论的多次波预测算法，由于多道加权求和、克希霍夫积分孔径的限制^[17]、二维侧面反射效应、观测误差和自由界面因子难以精确获取等多种因素，多次波预测的精确度难以得到保证，为此本文将时域变换道引入曲波变换域，实现了一种曲波域多道最小平方滤波的多次波匹配衰减方法，该方法结合时域变换对多次波预测误差的适应性和曲波变换对一次波与多次波的分离优势，有效解决了复杂海域中多次波预测不准确导致的多次波匹配衰减精度低的问题。

1 方法原理 1.1 曲波变换

Candès等人分别提出了第一代和第二代曲波变换(Curvelet)^[18-20]。第一代曲波变换是通过脊波变换间接实现的，尽管其在提取图像边缘信息方面表现出色，但由于其结构复杂、冗余度高以及计算量庞大，限制了该方法的进一步发展与应用。直到Candès等人摒弃了传统的脊波变换路径，直接构建了新的曲波基函数，从而提出了全新的曲波构造方法，这也标志着第二代曲波变换的诞生^[19-20]。与第一代相比，第二代曲波变换在结构上更加简洁，计算效率显著提升，相应的变换公式为

$ \boldsymbol{c}(j, l, \boldsymbol{k})=\iint \boldsymbol{d}(x, t) \overline{\boldsymbol{\varphi}_{j, l, k}(x, t)} \mathrm{d} x \mathrm{~d} t $

(1)

式中：$\boldsymbol{d}(x, t)$为地震记录；$\boldsymbol{\varphi}_{j, l, \boldsymbol{k}}(x, t)$为曲波基函数；$\boldsymbol{c}(j, l, \boldsymbol{k})$为曲波系数，j、l、k为曲波的尺度、方向和位置，对于二维曲波变换k包含两个维度，分别记作k_x与k_y；x、t分别为空间位置与旅行时。

Candès等^[21]提出了用于快速离散曲波变换的卷绕(Wrapping)和USFFT算法。由于地震记录可视为沿地震道方向与时间方向均离散采样的二维矩阵，因此可采用卷绕算法，公式(1)经离散化可得

$ \boldsymbol{c}(j, l, \boldsymbol{k})=\sum\limits_{x=0}^{x_{\max }} \sum\limits_{t=0}^{t_{\max }} \boldsymbol{d}(x, t) \overline{\boldsymbol{\varphi}_{j, l, \boldsymbol{k}}(x, t)} 。$

(2)

式中x、t均为离散采样，且0≤x≤x_max、0≤t≤t_max。

曲波在时空域中呈现中间粗、两边细的“针状”形态，并且在某固定方向震荡、其它方向是平滑的。正是由于曲波存在方向性强的特点，当沿着地震同相轴方向时相关性好因而能获得极大的曲波系数；而当方向与地震同相轴方向逐渐偏离时，对应的曲波系数值迅速衰减；在与地震同相轴垂直时曲波系数达到最小，其数值接近于零，因此曲波变换能够较好地提取地震数据的方向特征，实现地震数据的稀疏表示。由于多次波与一次波同相轴在时空域普遍存在相交甚至部分重叠，因此时空域中难以在有效保持一次波的前提下消除多次波干扰，而曲波变换则可根据相交同相轴方向上的差异达到显著分离二者的目的。

1.2 多次波记录的时域变换

曲波变换属于多尺度的数学变换，其中，小尺度参数提取的为数据的低频成分，用以描述数据的概貌；而大尺度参数提取的为数据的高频特征，反映地震同相轴的细节特征，因此若要在曲波域有效去除多次波，则预测多次波记录与原始数据应具有相近的尺度特征，即二者的振幅谱要基本一致。但是，利用波动方程法进行多次波预测时，由于输入记录中多道加权求和的计算过程，使得多次波记录的频带宽度明显变窄，且主频由高频端向低频端移动，从而影响了多次波的压制效果。

为了更好的匹配原始记录中的多次波成分，本文应用傅里叶变换实现基于预测的多次波记录构建主频具有明显差异的多次波数据，然后选择与原始记录中多次波匹配最佳的记录进行自适应衰减。

输入时空域的多次波记录m(x, t)，其傅里叶正变换与反变换可如下表示：

$ \left\{\begin{array}{l} M(x, \omega)=\int\limits_{-\infty}^{+\infty} m(x, t) e^{-j \omega t} \mathrm{~d} t \\ m(x, t)=\frac{1}{2 \pi} \int\limits_{-\infty}^{+\infty} M(x, \omega) e^{j \omega t} \mathrm{~d} \omega \end{array}\right.。$

(3)

式中：ω为角频率；M(x, ω)表示频率域的多次波记录。

根据傅里叶变换的性质，对m(x, t)求时域积分等价于频率域多次波记录除以jω，即

$ \int\limits_{-\infty}^t m(x, t) \mathrm{d} t \leftrightarrow \frac{1}{j \omega} M(x, \omega)。$

(4)

多次波记录m(x, t)的n次时域积分结果可表示为

$ m_n^I(x, t)=\frac{1}{2 \pi} \int\limits_{-\infty}^{+\infty}(j \omega)^{-n} M(x, \omega) e^{j \omega t} \mathrm{~d} \omega 。$

(5)

式中要求n>0。

同理，对m(x, t)求时域微分等价于频率域多次波记录M(x, ω)乘以jω，即

$ \frac{\mathrm{d}[m(x, t)]}{\mathrm{d} t} \leftrightarrow j \omega M(x, \omega)。$

(6)

因此，多次波记录的n次时域微分可表示为

$ m_n^D(x, t)=\frac{1}{2 \pi} \int\limits_{-\infty}^{+\infty}(j \omega)^n M(x, \omega) e^{j \omega t} \mathrm{~d} \omega。$

(7)

式中要求n>0。

综合公式(5)与公式(7)可以统一多次波记录的时域微积分过程：

$ \tilde{m}^q(x, t)=\frac{1}{2 \pi} \int\limits_{-\infty}^{+\infty}(j \omega)^q M(x, \omega) e^{j \omega t} \mathrm{~d} \omega。$

(8)

式中$ \tilde{m}^q(x, t) $表示多次波记录的时域积分或时域微分结果。

在公式中，当q < 0时，$ \tilde{m}^q(x, t) $表示时域积分记录，其主频向低频端移动；而q>0时，$ \tilde{m}^q(x, t) $为时域微分记录，主频向高频端移动；当q=0时，$ \tilde{m}^q(x, t) $表示多次波记录本身。

1.3 基于时域变换道曲波域滤波的多次波自适应相减

为了使存在预测误差的多次波干扰有效匹配衰减，设置时域变换阶数q的范围q∈[q₁, q₂]，利用预测的多次波记录m(x, t)创建系列时域积分与时域微分记录$ \tilde{m}^q(x, t) $。多次波的自适应相减是在曲波域而非时域中完成的，故而要根据公式(2)将原始地震记录d(x, t)以及系列多次波记录转换至曲波域，获得原始地震记录和多次波记录在曲波域中的系数C_d与$ \widetilde{\boldsymbol{C}}_m^q $，以便在该域中进行多次波的自适应相减处理。由于位置变量k包含两个维度，对于尺度为j与方向为l的曲波系数，其属于水平方向为k_x、垂直方向为k_y的二维复矩阵。为保证滤波过程的稳定性，针对任意位置k(k_x, k_y)的曲波系数值，以目标点为中心，构建一个边长为n的方形窗口，并从中提取数据以形成数据块C_d与$ \widetilde{\boldsymbol{C}}_m^q $，各数据块均包含n²个样点。

$ \boldsymbol{C}_d=\left|\begin{array}{cccc} c_{11} & c_{12} & \cdots & c_{1 n} \\ c_{21} & c_{22} & \cdots & c_{2 n} \\ \vdots & \vdots & \vdots & \vdots \\ c_{n 1} & c_{n 2} & \cdots & c_{n n} \end{array}\right|, \widetilde{\boldsymbol{C}}_m^q=\left|\begin{array}{cccc} \tilde{c}_{11}^q & \tilde{c}_{12}^q & \cdots & \tilde{c}_{1 n}^q \\ \tilde{c}_{21}^q & \tilde{c}_{22}^q & \cdots & \tilde{c}_{2 n}^q \\ \vdots & \vdots & \vdots & \vdots \\ \tilde{c}_{n 1}^q & \tilde{c}_{n 2}^q & \cdots & \tilde{c}_{n n}^q \end{array}\right|。$

(9)

式中，c_ij与_ij^q分别所截取记录块C_d与$ \widetilde{\boldsymbol{C}}_m^q $中的样点，i、j表示样点序号(1≤i≤n; 1≤j≤n)。

通过相关分析的方法确定与C_d匹配最佳的多次波分量，令$ \widetilde{\boldsymbol{C}}_m^q $与C_d的互相关函数为

$ \boldsymbol{w}_q=\sum\limits_{i=1}^n \sum\limits_{j=1}^n \boldsymbol{c}_{i j} \tilde{\boldsymbol{c}}_{i j}^q。$

(10)

式中，w_q为$ \widetilde{\boldsymbol{C}}_m^q $与C_d的互相关，而q∈[q₁, q₂]。

根据|w_q|数值的大小进行排序，据此确定最大值|w_q|_max时的q_max，然后利用相应的多次波数据块为$ \widetilde{\boldsymbol{C}}_m^{q_{\max }} $进行多次波自适应相减。滤波因子f可根据最小平方准则获得，即

$ \boldsymbol{f}=\underset{f}{\operatorname{arg}}{\rm{min}}\left\|\boldsymbol{C}_d-\widetilde{\boldsymbol{C}}_m^{q_{\max }} \boldsymbol{f}\right\|_{2 }。$

(11)

式中f为滤波因子向量$\left[f_1, f_2, f_3, \cdots f_q\right]^{\mathrm{T}}$。

由于利用最佳匹配的多次波记录块$ \widetilde{\boldsymbol{C}}_m^{q_{\max }} $进行自适应相减，滤波算子f的长度为1。对于基于L₂范数的最小平方滤波需使下式中的误差能量平方和e值为最小，即

$ e=\left\|\boldsymbol{C}_d-f \cdot \widetilde{\boldsymbol{C}}_m^q\right\|_2 \text { 。} $

(12)

式中f表示长度为1的滤波算子f。

由于滤波算子f的长度为1，从而避免线性方程组的求解而直接计算出f，即

$ f=\frac{\left(\widetilde{\boldsymbol{C}}_m^q\right)^{\mathrm{T}} \boldsymbol{C}_d}{\left(\widetilde{\boldsymbol{C}}_m^q\right)^{\mathrm{T}} \widetilde{\boldsymbol{C}}_m^q+\varepsilon}。$

(13)

式中：$ \left(\widetilde{\boldsymbol{C}}_m^q\right)^{\mathrm{T}} $表示$ \widetilde{\boldsymbol{C}}_m^q $的转置矩阵；ε表示阻尼系数(ε∈(0, 1)，一般取0.001)。

求出滤波因子之后，即可实现针对原始数据曲波系数C_d的多次波匹配相减处理，即

$ \boldsymbol{C}_p=\boldsymbol{C}_d-f \cdot \widetilde{\boldsymbol{C}}_m^q 。$

(14)

式中C_p表示曲波域数据块的自适应相减结果。

针对原始地震记录所有尺度、方向的曲波系数矩阵C_d，对其中每个样点重复公式(9)—(14)的处理步骤，即可获得消除多次波的曲波系数矩阵；将其反变换回时空域可得到压制多次波后的地震记录。

2 理论模型实验

建立多次波同相轴与一次波同相轴具有较小时差的理论模型，然后采用主频为35 Hz的雷克子波作为震源，通过弹性波方程的有限差分模拟得到含有多次波的原始地震记录。图 1(a)展示了原始记录，每炮200道，最小偏移距为0 m，道间距为10 m。应用基于反馈环理论的自由界面多次波预测方法，可得到仅包含多次波的炮集记录如图 1(b)所示，箭头所示位置即为多次波同相轴，其与有效信号时差较小。

使用图 1(b)中的多次波记录，分别对图 1(a)的炮记录采用常规曲波变换和本文方法进行多次波剔除，处理结果分别如图 2(a)和图 3(a)中所示，这两种方法剔除掉的多次波分别显示在图 2(b)和图 3(b)中。对比这两种方法的多次波衰减结果能够看出，在一次波与多次波同相轴交叉位置(虚线圈内)，由于同相轴过于接近，同时由于预测多次波与实际多次波记录存在波形相位上的差异，常规曲波变换多次波剔除结果中存在明显的多次波残余(箭头指向的近道位置)；而本文基于时域变换道曲波域滤波的多次波匹配衰减方法，即使在多次波预测不完全准确的情况下，仍能有效地压制多次波干扰，并保持有效波同相轴的连续性。对比剔除的多次波记录可知，本文提出的方法在多次波压制方面具有更好精度。

3 野外资料处理实验

本文选用的二维测线ML_A位于海底硬度高、地形复杂的区域，水深在1 100~1 650 m之间，该测线采集数据包括2 362炮，每炮包含180道，且道间距和炮间距均为26.6 m，最小偏移距为212.8 m，测线原始炮集记录中存在强烈的自由界面多次波和其他复杂的多次波。

在预测自由界面多次波之前，对ML_A测线的原始炮集进行预处理，包括去除直达波、滤除低于3 Hz的低频信号以及球面扩散补偿等。预处理后的第760炮记录如图 4(a)所示，预测得到的仅包含多次波的第760炮记录如图 4(b)所示。

根据图 4(a)中的原始炮集记录与图 4(b)中的预测多次波记录，分别采用常规曲波变换法和时域变换道的曲波域滤波法进行多次波自适应衰减，两种方法使用相同的窗口大小(10×10)，后者的时域变换参数q范围为-1≤q≤2。图 5(a)和图 5(b)显示了使用两种方法消除多次波后的炮集记录，去除的多次波干扰分别如图 6(a)与图 6(b)所示。对比可见，本文方法处理后期海底及下层波阻抗界面的强多次波已大部分消除，而常规曲波变换的处理结果中仍存在明显的多次波残留(矩形区域所示)。

分别基于原始炮集记录、常规曲波变换多次波剔除后的炮集记录以及应用本文方法多次波剔除后的炮集记录进行克希霍夫积分叠前时间偏移处理，得到的剖面分别见图 7—图 9。

对比图 7—图 9可知，在一次波偏移速度场下进行成像时，尽管多次波已经得到了部分压制，但在图 7的虚线框区域，仍然可以看到明显的多次波同相轴。通过对比图 8和图 9中两种多次波压制方法的结果，可以发现，基于时域变换道曲波域滤波的多次波匹配衰减方法(即本文方法)所得的剖面中，虚线框内的强多次波同相轴几乎完全消除，而常规曲波变换处理后的剖面中仍然存在明显的多次波残余，由此进一步验证，本文提出的方法在复杂崎岖海底环境中，相较于传统的曲波变换方法，不仅能够更好地保留一次波信息，还显著提升了多次波的压制效果，从而可获得更加清晰和精确的地震成像结果。

4 结语

本文提出一种时域变换道曲波域滤波的多次波匹配衰减方法，此方法充分利用了时域变换对波形和相位的良好适应性以及曲波变换在有效分离一次波和多次波方面的优势。理论模型测试以及实际数据处理结果都表明，该方法在多次波剔除效果上优于传统的曲波变换方法，其能够在复杂海底地形(如崎岖海底区域)中有效地压制多次波，同时尽量避免对一次波的损伤，更适合复杂地质构造区域的多次波剔除需求。