0 引言

海洋地震勘探中，与表面相关多次波^[1-2]通常夹杂在一次反射波之间，难以去除，因此对与表面相关多次波的预测一直是国内外地震数据处理中的关键性问题^[3]。近些年发展起来的自由表面多次波压制(SRME)^[4-5]技术, 在多次波预测方面表现突出。这种方法利用数据驱动^[6]理念，在不需要任何地下介质先验性信息和假设的前提下，利用地震数据本身反馈迭代实现多次波预测，并可获得较好的效果。SRME算法最早由Berkhout和Verchuur提出，最初的算法是基于反馈模型理论^[7]的2D级数展开法压制多次波^[8]，之后两人又基于反馈模型理论提出了迭代法压制表面多次波算法^[9-10]，这种迭代的算法在很大程度上提高了计算精度，这就是最初的2D SRME算法。

虽然2D SRME方法在有些条件下可以有效压制表面多次波，但其算法未考虑多方向、多次波的贡献，在大多情况下其处理效果难以满足要求，这使得人们向3D SRME探索。1997年，Van Dedem等^[11]在SEG (society of exploration geophysicists)会议上阐述了3D算法压制多次波的必要性。2002年，Van Dedem等^[12]利用稀疏反演方法进行三维多次波预测，Levin等^[13]根据近似动校正(NMO)技术实现了零偏移距的3D SRME算法。2009年，Van Groenestjin等^[14-15]进行了稀疏反演方法多次波预测的研究，并实现了近偏移距数据的重构。2010年，Dragoset等^[16]指出，在所有多次波预测方法中，3D SRME方法最能广泛压制表面相关多次波。常规的3D SRME理论技术在国外已经发展成熟，近些年，更多的研究趋于将三维SRME技术应用到实际地震数据处理中。

然而，在国内，3D SRME理论技术刚刚成型，基本都处于理论模型试算阶段。以往理论中给出的观测系统都是理想的，观测网格为高密度采样地震数据，但是这样的数据与实际观测获得的地震数据相差甚远。海洋地震数据采集中，缆与缆之间的距离往往都在200 m以上，这样的采样导致地震数据在横测线方向过于稀疏，不能有效描述横测线方向的地下构造，难以有效应用3D SRME方法。

常规SRME处理过程中，采用最小平方匹配法进行多次波压制。这种方法起步较早，发展到现在已经基本成熟，其算法核心是通过分窗口进行多次波与原始数据的匹配相减。匹配相减原则是对窗口内数据求最小平方差，当数据满足输出要求时输出。虽然最小平方匹配发展得比较成熟，但是这种方法处理的结果使一次波的能量受到损失^[17-18]，影响地震数据处理的保真性。

针对上述问题，本文提出三维稀疏反演预测及曲波域匹配衰减技术。对于地震数据在横测线方向过于稀疏的问题，采用三维稀疏反演法多次波贡献道集重建，有效利用三维数据进行多次波预测。引入曲波域匹配相减方法^[19]，多次波在曲波域有良好的自适应性和光滑性，匹配过程中能够有效地避免数据的过匹配，从而避免一次波能量的损失。

1 稀疏反演多次波预测方法基本理论

SRME方法是利用地震数据自身褶积实现多次波预测，2D与3D方法的矩阵表示形式是类似的，具体区别在于矩阵的维数，如：

式中：表示震源点坐标为(x_s, y_s)、检波点坐标为(x_r, y_r)的预测多次波，它是由(x_k, y_k)点激发、(x_r, y_r)点接收的共检波点道集R(x_r, y_r, ω; x_k, y_k)与(x_s, y_s)点激发、(x_k, y_k)点接收的共炮点道集D(x_k, y_k, ω; x_s, y_s)地表一致性褶积获得的(时间域的褶积在频率域为乘积)；ω表示频率；x表示纵测线方向，y表示横测线方向。当输入的地震数据只包含x方向的数据时，即数据体为二维时，等式(1)为2D SRME算法；当输入的地震数据是包含x方向和y方向的三维数据体时，等式(1)为3D SRME算法。3D SRME将一个三维共检波点道集数据与一个三维共炮点道集数据进行多维褶积，得到的结果M_xy(x_s, y_s, x_k, y_k, x_r, y_r, ω)为三维多次波贡献道集，将三维多次波贡献道集再进行累加实现多次波预测。式(1)中的多次波贡献是分两步进行求取的。第一步，将多次波贡献道集沿纵测线方向进行累加，转换为横测线方向的多次波贡献道集：

第二步，沿横测线方向进行累加，得到三维预测多次波：

而三维海洋拖缆地震采集系统所获得的地震数据沿横测线方向是稀疏的，这种稀疏采样直接影响三维多次波预测算法第二步的累加结果。

图 1a、b分别展示了高密度采样与稀疏采样沿横测线方向的多次波贡献道集，图 1c给出两种道集累加预测的多次波与原始数据的对比。可以明显看出，稀疏采样所预测的多次波结果过于粗糙，而且引入了噪音干扰；这也是实际生产中遇到的难题。面对稀疏采样的地震数据，需要寻找一种合理的数据处理方法来有效地实现3D SRME。从图 1a、b可知，多次波贡献道集的时距曲线是双曲线型，而多次波的主要贡献来自于双曲线的顶点处，这个贡献集中的部分定义为多次波贡献的菲涅尔带，即图 1a、b中矩形框定的范围。在多次波预测过程中，菲涅尔带的数据决定了多次波的波至时间和振幅能量；菲涅尔带的多次波贡献对于准确的三维多次波预测是十分必要的。假设对稀疏采样的多次波贡献道集菲涅尔带进行了有效的信息拾取，然后对稀疏波场进行重建，就能够实现多次波预测。本文采用稀疏反演法对横测线方向上多次波贡献道集的菲涅尔带进行稀疏反演重构。

对于稀疏采样的横测线方向多次波贡献道集，稀疏反演方法采取一种稀疏变换来取代等式(3)中的累加过程。即将横测线方向多次波贡献道集进行参数化表示，目的是在稀疏的多次波贡献道集中提取有用的多次波信息，包括菲涅尔带顶点的走时(多次波贡献的走时)和位置(多次波贡献顶点的横测线方向的位置)，以及该双曲线的曲率(多次波贡献的曲率)。对于水平反射层和倾斜反射层，多次波贡献道集的时距曲线呈现双曲线型，都能够有效识别出菲涅尔带，所以这种转换是可以实现的。横测线方向多次波贡献道集参数化需要以下3个参数：多次波贡献道集顶点的旅行时τ、多次波贡献道集横测线顶点位置y₀以及多次波贡献道集双曲线顶端的曲率函数q。利用这些参数，可以将多次波贡献从数据空间d(y_k, t)转换到模型空间m(q, y₀, τ)：

该过程可以写成

式中：L为变换算子，L^T是L的共轭算子。同样，也可以建立从模型空间到数据空间的转换：

即

变换算子L的维度为数据点数与模型参数个数的乘积。

对于稀疏的地震数据，只有小部分横测线方向的多次波贡献可以被用作参数化输入(图 1b)。图 2是对图 1b数据进行初步稀疏变换(式(4))的结果。从图 2可见，输入数据的稀疏采样及模糊性导致了模型空间严重的混淆现象，从而降低了数据转换到模型空间的可靠性与保真性。因此，通过式(4)这种简单的变换所得到的模型空间数据远远达不到预期的效果。

等式(5)和(7)的变换是在高分辨率模型数据问题上提出的，所以在稀疏的地震数据上并不适用。解决上述模型空间模糊问题的办法是定义一个比低分辨率共轭算子L^T更好的反向模型估算算子L。为了这个目的，需要将数据空间到模型空间的转换重新定义为一个反问题。数据空间与模型空间之间包含噪声的转换关系如下：

式中，n代表数据中的噪声，假定为高斯随机噪声。解决这个自然不适定问题需要整合一个先验信息来调整这个反问题。

将柯西概率密度函数及高斯正则化整合到等式(8)中，通过最小二乘反问题解法得到了一个最好的后验估算m：

式中：λ=σ_n²/σ_m²，σ_m为模型稀疏参数，σ_n为噪声估计参数；Q是一个N_m×N_m的对角矩阵，其中的元素为

式(9)近似于式(8)反问题的阻尼最小二乘解，其中的差别是阻尼项λQ^-1。结合先验信息将线性反问题转变为非线性，这个过程需要一个非线性优化法来求解。

m的非线性反演解法通过两个重要参数控制：σ_n和σ_m。第一个参数σ_n对出现在数据中的(高斯)噪声进行估计。这个估计包含数据中的随机噪声以及由建模的不确定性和测量方法的错误假设引起的噪声。第二个参数σ_m控制解法中稀疏程度。当比值m_i²/σ_m²≫1时，阻尼项λQ^-1作用到模型参数m_i上的阻尼较小，趋于稳定；当m_i²/σ_m²∝0时，模型参数取值将接近零。这个机制只提高与数据本身相对应的模型参数，因此提高了模型空间的分辨率。这两个参数是通过λ连接到一起的，并将控制反演的稳定性(或是权衡稀疏性和数据之间的拟合)。反演过程中参数σ_n固定不变，参数σ_m用于调整反演精度。

如图 3所示，通过在模型空间参数上应用柯西概率密度函数，不仅调整了反演，而且提高了模型空间的分辨率，大多数离群的参数值都被消除，得到了拟合逼近更好、更加真实可靠的模型空间参数。另外，对模型空间引入更高级的高斯正规化，使得模型空间的能量得到了平滑。可见，在模型空间参数上应用柯西概率密度函数及高斯正则化，使得数据在模型空间模糊的问题得到了解决，达到了理想效果。

将模型空间中的多次波信息转换为数据空间的多次波数据，最简单的方法是用一个模型操作符将模型空间中的多次波信息转换为数据空间高密度采样的横测线方向多次波贡献轨迹，之后对其累加就可以得到预测的多次波。然而，通过利用固定相位近似法，可以直接从模型空间信息中获取多次波数据。

固定相位近似法应用于多次波预测是成功的，它对模型空间数据进行振幅和相位校正。这个校正是对模型空间参数q和τ进行的：

将校正的模型空间轨迹进行累加，就可以直接得到预测多次波数据：

式中，为利用固定相位近似法得到的预测多次波数据。

在稀疏反演多次波预测过程中，通过一个参数化的模型反演，对垂直测线方向多次波贡献道集中菲涅尔带的累加进行替换，使其转换到模型空间中，再对模型空间中多次波轨迹的振幅和相位做出校正，就达到了多次波准确预测的目的，至此三维稀疏反演预测多次波基本完成。

2 曲波域多次波匹配相减

常规的多次波处理过程分为两步。第一步，预测多次波。通过一些多次波预测方法将地震数据中的多次波预测出来。第二步，多次波的匹配减去。找到一个多次波与地震数据的匹配算子，通过这个匹配算子将多次波从输入数据中减去。

传统SRME方法中的多次波去除过程基于最小平方匹配相减法。该方法是一种自适应减去滤波法，对窗口内的预测多次波与原始数据求最小平方差，当数据满足输出要求时就输出，当不满足要求时就被压制掉, 如图 4a所示。这种方法虽然提出得相对较早，发展也趋于成熟，但是匹配过程中会对地震数据中的非多次波带来干扰，或引入新的噪声，或使有效数据发生缺失。这样将使得处理后的地震数据一次波部分能量受到损失，影响地震数据处理的保真性。为了更好地适应复杂地震数据的多次波匹配问题，引入曲波域匹配。曲波域在地震数据处理中有很好的表现，由于其自身具有较好的多尺度、多方向性及圆滑性，能够对地震数据进行更细致的划分，分尺度、分方向地对地震数据进行处理。此外，曲波变换对二维分段光滑函数具有更好的非线性逼近性，稀疏性更好，这样能够更好地贴近多次波数据的同相轴，同时不影响邻近的数据，如图 4b所示。所以，相较于最小平方匹配，曲波域匹配方法能够获得更好的多次波匹配结果。

曲波域匹配滤波基本流程为：首先，对原始地震数据做一个子波匹配；然后，将子波匹配的多次波转换到曲波域，在曲波域内乘以一个对角化的曲波系数；最后，进行曲波域反变换，就得到了曲波域匹配的多次波。

多次波匹配完成之后，采用软阀值法将匹配的多次波信息从地震数据中去除。软阀值的基本原理是：将小于阀值的系数赋予零值；将大于阀值的系数减去一个阀值大小的量。具体表达式为

式中：C_p为原始地震数据经曲波变换得到的Curvelet系数；T_m为经过曲波域匹配后的多次波Curvelet系数的模。

3 理论模型试验 3.1 三维稀疏反演预测多次波

为了验证本文论述的三维稀疏反演方法对于稀疏情况下多次波预测的有效性，采用三维地震反射记录进行试算。三维介质的速度模型如图 5所示。该模型包含3个地下反射层，每个反射界面都是多方向倾斜的。在地表观测系统中布置5条测线，测线长2 000 m，测线间距为200 m，测线上道间距为25 m、炮间距为25 m，获得一个2 048×81×81的三维数据体。

利用数据矩阵自身褶积实现多次波贡献道集的求取，对每条测线上的多次波贡献累加求取对应检波点横测线方向上的多次波贡献道集。选取0测线上1 500 m处数据点对应的横测线多次波贡献(图 6a)，通过稀疏反演将数据空间稀疏的多次波贡献转换到模型空间(图 6b)。图 6a中箭头所指处即为每条多次波贡献同相轴的菲涅尔带顶点，图 6b模型空间多次波能量出现的位置与此相同，验证了稀疏反演法能够有效识别出所需的多次波能量，并将它重构出来。

将2D SRME预测的多次波、原始数据及模型空间累加预测的多次波进行对比(图 6c)，可以明显看出两种方法预测的多次波在波至时间上的差别：2D SRME预测的多次波与原始数据存在时间差，出现了波至时间延后的情况，如2.5 s处；在稀疏反演法预测的结果中，多次波波峰出现的位置与原始数据一致，克服了常规2D SRME算法中的弊端，实现了多次波信息的准确预测。

选取0测线上1 000 m炮点处数据，截取多次波发育良好的部分，分别进行2D SRME和三维稀疏反演法多次波预测，结果如图 7所示。通过对比图 7a、b、c可以看出，2D SRME方法预测的多次波整体出现了时间延后的现象，而三维稀疏反演法预测的多次波保证了波至时间与原始数据一致；对比图 7中的黑色箭头处可以更明显地看出预测结果的差别。

为了更有效地展示预测结果的差别，将图 7中箭头处进行放大(图 8)。对比发现：2D SRME预测结果中出现了错误的交互关系，与原始数据不符；在三维稀疏反演法的预测结果中，多次波间的交互关系与原始数据完全一致，准确地预测了深层较弱的多次波，这将为有效压制多次波信息提供可靠的依据。

3.2 曲波域多次波匹配减去

为验证曲波域多次波匹配减去的有效性，对数据分别进行最小平方法匹配和曲波域多次波匹配。从图 9可以看出：经最小平方匹配后，一次波位置出现了噪音影响(图 9c)；相比于此，经曲波域匹配后的多次波在原本一次波的位置未引入其他干扰信息，保证了数据处理的可靠性，有效避免了地震数据受匹配过程影响而引入虚假信息(图 9d)。

在将多次波从原始数据中去除的过程中，采用3种不同的组合方法进行测验。第一种，将最小平方匹配的多次波利用最小平方相减法进行去除(图 10a)，处理结果中大部分多次波信息得到了有效压制，仅残留了少量的多次波痕迹；但是部分一次波信息受到了损伤，带来了过匹配问题。第二种，将经过曲波域匹配的多次波进行最小平方相减(图 10b)，此次处理没有损伤一次波信息，多次波的压制也取得了一定的效果；但还是残留了一小部分多次波的痕迹。第三种，将经过曲波域匹配的多次波采用曲波域阈值相减(图 10c)，处理结果在确保了一次波不受损伤的情况下，近乎完全地去除了多次波信息，相较于前两种匹配减去方法，曲波域匹配减去法在对于一次波能量的保护和多次波的压制方面，都取得了良好的效果。

4 结论

1) 本文提出了三维稀疏反演多次波预测方法。通过引入柯西概率密度函数和高斯准则法，对稀疏的横测线多次波贡献道集进行了参数化表示，避免了模型空间模糊现象，有效实现了多次波贡献道集在模型空间的重构，可以达到准确预测三维多次波的目的。

2) 采用理论数据对三维稀疏反演多次波预测方法进行测试，并与传统2D SRME方法处理效果进行对比的结果表明，三维稀疏反演多次波预测方法对稀疏的地震数据具有更好的处理效果，实现了多次波信息的准确估计。

3) 本文采用曲波域多次波匹配减去法对多次波进行去除。由于地震数据在曲波域具有多尺度、多方向以及圆滑的特性，在曲波域对多次波进行匹配，可以有效避免数据过匹配的问题。另外，曲波变换对二维分段光滑函数具有更好的非线性逼近性, 在曲波域进行多次波匹配并结合曲波阀值法进行多次波去除，可以使多次波去除得更彻底。

综上所述，本文提出的三维稀疏反演多次波预测方法及曲波域多次波匹配减去法具有良好的理论基础，避免了传统方法的弊端，可以对海洋地震数据进行准确有效的多次波压制。