多次反射波常出现于海洋地震勘探资料中，而且通常能量很强，从而对一次波的处理和成像造成强烈干扰，所以对当前广泛应用的成像方法而言，多次波压制预处理尤为重要(Jakubowicz，1998；李翔和胡天跃，2009；王保丽等，2014).在众多的多次波压制方法中，基于数据驱动的3D SRME(Surface-related multiple elimination)方法被认为是最有可能成功预测和压制表面多次波的方法(Dragoset et al.，2010).基于Berkhout(1980)提出的数据矩阵的 思想，该方法最初由Verschuur等提出(1992).依据 惠更斯原理(Berkhout and Verschuur，1997; Verschuur and Berkhout，1997; Sun，1999)，SRME方法分为两步来实现：一是多次波预测；二是多次波自适应相减.而多次波预测过程又可分为两步：一是计算多次波的贡献道集(Kabir and Verschuur，1995；van Dedem and Verschuur，2001；Hokstad and Sollie，2003)；二是对贡献道集进行叠加，形成预测的多次波道集(van Dedem and Verschuur，2005).SRME方法虽基于数据驱动，可处理复杂构造地震资料，但有效预测多次波的一个重要假设是输入全地震数据(Verschuur，1991；Dragoset et al.，2006)，也就是说每个检波点位置都有炮点(炮点间隔等于检波点间隔).常规观测系统采集的海洋二维地震资料，在进行简单的数据规则化后，上述假设条件就基本得到满足，所以在二维地震资料多次波压制处理中，SRME方法得到成功应用.

从原理上来讲，2D SRME方法可直接扩展到3D，而且三维 SRME方法避免了二维算法的所有假设(Matson and Abma，2005；Dragoset et al.，2010).然而，当前海洋三维地震资料采集的纵测线稀少，因此横测线方向的采样间隔较大，这就存在计算成本高、存储以及空间采样不足等问题，致使全三维SRME技术一直未能得到有效应用.

将SRME方法应用于三维地震资料的预处理，目前常规的处理方法有两种：一是在多次波压制前，对三维数据做规则化处理，重建缺失的炮记录和道记录，进而依据SRME基本理论，实现全三维多次波预测.主要缺点是在叠前会形成大量的地震数据，占用大量磁盘空间，预测计算时，也会耗费大量计算 机时，增大地震数据的处理周期(Bisley and MacNeill，2008). 而且叠前数据的规则化处理，也存在较大的误差(Dragoset et al.，2010)；二是依据二维理论，通过诸多的假设、振幅匹配以及乘经验因子等方法将三维数据利用二维算法近似预测多次波，其有悖于地下介质的真实情况，最终使得多次波压制效果不佳(Dragoset and Jeričević，1998；Kelamis and Verschuur，2000).

Verschuur(1991)和Ross等(1999)等将三维地震资料按照二维方法进行多次波压制预处理，也就是说对不同的检波线依次进行处理，在完成简单的数据规则化后，可完成依据波动理论的多次波压制处理.很显然，该方法不能充分考虑三维介质的倾角和方位角等信息，对于地层构造简单的地区可以应用，但是对于稍微复杂的地层，该方法在多次波预测的走时、相位和振幅等诸多方面都存在很大误差，难以取得理想的多次波预处理结果(van Dedem and Verschuur，2005).因此，进行全三维SRME实用方法研究是地震勘探的迫切需求.van Dedem和Verschuur(2001)、Hokstad和Sollie(2003)提出贡献道集概念，为3D SRME的有效实现探索新的思路.本文借鉴贡献道集的思想，利用基于稀疏反演方法实现3D SRME表面多次波压制.

依据实际海洋地震资料采集的特点计算三维地震资料表面多次波的贡献道集.先对纵测线进行道集记录的简单外推和规则化后，形成可按照二维方法求取贡献道集的地震数据，横测线方向数据很稀疏，但是横测线方向地震同相轴可视为双曲线或抛物线，因此可将该同相轴沿给定的曲线路径进行积分，进而可进行反演计算，以代替横测线方向的求和过程，降低3D SRME对运算数据的要求，避免数据重建的处理，减少贡献道集菲涅尔带求和过程中的空间假频.由于抛物Radon变换求解矩阵小，频率域解耦和计算效率高等优点，本文针对三维地震数据横测线稀疏的情况，并依据其抛物线型同相轴的假设，进行稀疏反演三维表面多次波压制方法研究，并对理论和实际模型进行了试算，显示了多次波压制方法的效果和实用性.

Verschuur(1992)等提出了压制表面多次波的2D SRME级数展开法，即：

其中，P(ω)为只含有效波的地震数据，D(ω)为含多次波的地震数据，A(ω)为与表面反射系数r₀和震源子波 S(ω)相关的表面算子，公式为

其中，r₀为表面反射系数，S ^－1(ω)为逆源子波.

从式(1)和(2)可知，2D SRME方法依赖于原始地震数据、表面反射系数以及震源子波信息.但是实际上震源子波未知，Verschuur等(1992)假设多次波压制后地震数据能量最小，利用最小二乘方法实现了多次波的自适应衰减.

理论上讲，级数展开法要求预测所有阶表面多次波，而且要求数据矩阵的多次相乘计算.基于多次波形成的物理原理，Berkhout和Verschuur(1997)等给出了多次波压制的SRME迭代计算法，迭代形式如下：

其中，i代表迭代次数.一般来说，首次迭代时，取 P ⁽¹⁾(ω)= D(ω)，类似于级数展开法，通过预测多次波与总波场的最小二乘匹配可估计出震源子波，进而实现多次波压制处理.

在公式(3)中，定义 IR(ω)代表第i次迭代的多次波压制结果，即 IR(ω)= P ⁽ⁱ⁾(ω)，定义 M(ω)为预测多次波.则多次波预测公式可表述为

对任一频率，依据式(4)预测多次波需要计算矩阵乘法，也就是对矩阵 IR(ω)第i行和矩阵 D(ω)的第j列进行内积计算.

由公式(4)，单道二维地震数据的多次波预测可进一步表示为

将二维多次波预测方法直接扩展到三维，得到3D SRME多次波预测公式为

其中，x_r、y_r、x_s、y_s分别表示检波点与炮点的空间坐标，x_k、y_k分别表示求和空间范围内震源和检波点位置的坐标变量.

依据贡献道集的概念，式(6)中的求和可进一步用贡献道集作以说明，若令 M _xy为多次波贡献道集，则有

上式也可写成褶积表达式为

因此，应用贡献道集的思想，地震数据多次波预测的数学描述为

其中，m _xy为多次波贡献道集的时空域表示.不难 发现，SRME方法压制多次波可分为三个步骤：即贡献道集计算、贡献道集叠加和多次波自适应相减.

设计如图 1a所示的三维地质模型，测试分析二维和三维SRME方法的多次波预测和压制效果.模型最上层为海水，下面为含断层、陡倾角和尖灭等复杂构造的地层.深度域速度模型三个空间方向的网格数均为400，网格间距为5 m.设计Inline方向和Crossline方向均100个检波器接收的观测系统，也就是道间距和炮间距均为20 m，震源子波为主频25 Hz的零相位Ricker子波，时间方向采样间隔为4 ms，记录道长为2 s.按传统三维SRME多次波预测对地震数据的要求，在每个检波点都布设了炮点，共模拟10000炮地震记录.图 1b为正演的三维某单炮地震记录，可见能量较强的表面多次波.

图 1

Fig. 1

图 1 三维复杂模型与正演模拟炮记录 (a)三维速度模型；(b)正演模拟三维单炮记录. Fig. 1 hree dimensional complex geological model and modeling shot record (a)3D velocity model;(b)The 3D shot record modeled by wave-equation method.

文中应用公式(7)对图 1所示的复杂模型数据进行三维多次波的预测与压制试算，同时也应用公式(5)对抽取三维数据中一条二维测线(Inline50)进行二维 SRME的计算.2D和3D SRME算法预测和压制多次波效果如图 2所示.图 2a为正演模拟的含多次波Inline50线共零偏移距剖面，从图中可以看出含有能量较强的多次波.图 2b和2c分别为2D SRME和3D SRME算法的多次波预测结果，对比分析可知，3D SRME算法预测结果的同相轴更加连续，特别需要说明的是图中箭头所示位置处，适应三维复杂介质的全三维 SRME算法能够清晰预测出多次波.如前所述，当地层倾角较大或者较为复杂时，2D SRME对多次波预测的各种假设已经不符合实际地质条件，也就是说对于上述地震数据体，2D SRME算法预测的多次波出现较为明显的走时误差.3D SRME应用全三维的预测算法，其预测的多次波位置符合实际数据情况，因此预测复杂构造数据多次波时，3D SRME算法精度要明显高于2D SRME算法.图 2d和2e分别为采用相同自适应相减方法(石颖等，2013)得到的2D SRME和3D SRME多次波压制结果，由上两图中箭头所示位置可知，利用2D SRME方法多次波压制不彻底，多次波能量有较大残余，相比而言，3D SRME方法压制多次波的效果较为理想.

图 2

Fig. 2

图 2 2D和3D SRME预测及压制多次波对比 (a)含多次波的原始共偏移距剖面；(b)2D SRME方法预测的多次波；(c)3D SRME方法预测的多次波；(d)基于2D SRME方法预测的多次波压制结果；(e)基于3D SRME方法预测的多次波压制结果. Fig. 2 Comparison of multiple prediction and suppression using 2D and 3D SRME (a)The modeled common offset section;(b)Predicted surface-related multiple by 2D SRME;(c)Predicted surface-related multiple by 3D SRME; (d)Multiple suppression based on 2D SRME multiple prediction;(e)Multiple suppression based on 3D SRME multiple prediction.

对图 1所示的理论模型，设计观测系统的炮间隔与检波点间隔相等，因此模拟数据为全波场数据，符合SRME方法预测表面多次波的数据假设条件.对实际地震数据来说，其Inline方向采样较密，稍做预处理后，可满足多次波预测计算的要求，但是由于观测系统的限制，其Crossline方向通常较为稀疏，如果直接按照公式(9)所述方法预测多次波，将会产生严重的假频，通常采用数据插值和规则化的方法，获得全波场数据，再利用SRME方法计算，但是数据插值与规则化后，将产生更大的地震数据体，因此计算效率和存储成本是制约SRME方法有效应用的两个重要问题.不论是二维还是三维SRME方法，多次波预测的计算效率可通过应用CPU/GPU异构平台加速计算来解决(石颖等，2013)，对海量地震数据的规则化预处理而言，波场存储则具有很大的挑战性.

由公式(9)可知，多次波预测由贡献道集计算与叠加两步组成.对于实际的三维地震资料而言，预测单道多次波的贡献道集为三维数据体，可将贡献道集计算分为两个步骤.为此，定义如下数据向量：

公式(10)表示沿主线方向进行叠加求和后得到的数据向量.理论上，再按照Crossline方向进行叠加求和就得到最终预测的多次波.实际上，Crossline方向数据通常较为稀疏，无法进行求和运算.van Dedem和Verschuur(2001)证明公式(10)表示的数据的时距曲线可近似表示为如下双曲线形式：

其中，τ为截距时间，y_a为双曲线顶点位置坐标，v为速度，满足v=v_nmo/2，v_nmo为正常时差校正速度.

将式(11)进行Taylor展开，可得到抛物线方程为

其中，q=1/2τv²为曲率参数，公式(12)的假设条件 是勘探目标层的埋藏深度远大于Crossline方向的宽度，对多数实际三维地震数据，上述假设是可以满足的.

对于频率域的数据，公式(12)的抛物Radon变换形式可表示为

其中，q_i(i=1，2，…，N_q)是曲率参数，y_aj(j=1，2，…，N_ya)为抛物线顶点位置，上述方程用矩阵形式表示为

其中，L 为算子矩阵，上角标H代表复共轭转置，D 为数据矩阵.其反变换可表示为

公式(14)模型空间的估计可通过最小二乘反演问题的计算来实现.假设模型空间是稀疏的，在模型求解时应用Cauchy概率密度函数和Gauss-Cauchy准则，可得到稀疏的高分辨率的解(Sacchi and Ulrych，1995; Trad et al.，2003; Hokstad and Sollie，2003). 模型空间的最大后验估计，也就是其解的形式为

公式(16)为用于求解超定方程的加权最小二乘解(Wang，2003).但本文求解的是数据空间点数小于模型空间的情况，也就是所谓的欠定问题，根据Sacchi和Ulrych(1995)以及Hokstad和Sollie(2003)等给出的矩阵形式，可将公式(16)写为

公式(17)可用于欠定方程的求解.其中，I 为单位矩阵，λ=σ_n²/σ_m²，σ_n²表征地震数据中的噪声水平，σ_m²为控制解的稀疏性的参数. Q 为对角阵，其主对角元素为

式中，m_j为模型空间矩阵的第j个分量.由式(18)可知，式(17)所表示的为模型空间的非线性稀疏解，该方程可应用迭代法进行求解(Sacchi and Ulrych，1995; Hokstad and Sollie，2003).

应用非线性抛物稀疏反演替代多次波贡献道集稀疏的Crossline方向的叠加，可大幅减少三维数据的叠前规则化预处理.同时这里需要指出的是在得到预测的多次波数据时，不需进行模型空间数据的反变换和求和运算，而直接对稀疏的模型空间数据叠加即得到预测的多次波，但是这种三维多次波预 测结果存在相位误差.依据稳相近似理论(Bleistein，1984; Wapenaar，1992)，可应用如下公式进行相位校正：

模型空间数据求和后，应用公式(19)进行校正，就得到最终三维表面多次波预测数据为

上述基于稀疏反演的三维表面多次波预测主要包括四个步骤，即：沿Inline方向贡献道集求和，Crossline方向贡献道集数据的稀疏反演，反演模型空间数据相位校正以及校正后叠加，得到如式(20)所示的三维表面多次波预测结果.为了对比分析预测效果，文中也进行了多次波的自适应相减计算(Berkhout and Verschuur，1997; Abma et al.，2005; 石颖等，2013). 不同于常规的SRME方法，基于稀疏反演的表面多次波预测的方法，可避免Crossline方向的数据规则化处理，有效缓解数据存储的压力，降低三维SRME方法对地震采集数据的要求.

对所形成的基于稀疏反演三维表面多次波预测方法，文中也应用CPU/GPU并行加速技术提高三维数据体多次波预测和压制的计算效率.为验证稀疏反演3D SRME方法，设计了含四层水平层状介质的理论模型，正演模拟时采用全3D地震数据采集形式，观测系统包括56条炮线，在x方向和y方向分别放置56个检波器，每条炮线等间隔放置56个震源点，每炮3136道接收，道间距和炮间距均为15 m，时间采样间隔为4 ms.抽取地震数据中的5炮进行显示和分析(图 3)，图 3a为原始含多次波模型数据，图 3b为采用传统2D SRME方法预测的多次波，图 3c为利用传统3D SRME方法对全波场数据预测的多次波，图 3d为抽取全波场数据中四条测线(线间距为75 m)，利用稀疏反演方法预测的多次 波.分析可知，相比于2D SRME方法，传统3D SRME和基于稀疏反演的3D SRME方法改善了预测效果.此外，抛物稀疏反演的3D SRME方法避免了Crossline方向的大规模数据规则化，极大地降低了计算成本，节约了大部分存储空间.

图 3

Fig. 3

图 3 表面多次波预测与对比 (a)原始炮数据；(b)2D SRME方法预测；(c)传统3D SRME方法预测；(d)抛物稀疏反演3D SRME方法预测. Fig. 3 Surface-related multiple prediction comparison for different prediction approaches (a)Shot records with multiple；(b)Predicted multiple shot gathers using 2D SRME；(c)Predicted multiple shot gathers using traditional 3D SRME；(d)Predicted multiple using the proposed sparse inversion 3D SRME prediction algorithm.

为更好理解本文给出的稀疏反演表面多次波预测原理，图 4给出了三维地震数据中单道多次波预测的对比分析.该方法可有效避免数据规则化，特别是稀疏Crossline方向数据的大规模重建处理.基于图 3抽出的四条横测线数据(线间距为75 m)，在完成沿主线方向的贡献道集计算与叠加后，图 4a给出了某一接收点处 Crossline贡献道集剖面，在稀疏反演三维表面多次波预测计算中，该道集就是稀疏反演的输入数据，根据公式(17)计算，稀疏反演的模型空间数据如图 4b所示.图 4c为单道多次波预测数据的对比，重点对比密集Crossline贡献道集叠加(传统3D SRME方法)和模型道集叠加(稀疏反演方法)多次波预测结果.其中第1道为正演模拟的含多次波的原始道数据；基于全波场数据，利用传统3D SRME方法预测的多次波如图 4c中第2道所示；第3道为输入的数据Inline网格不变，但只有4条主线(横测线方向稀疏)，采用本文稀疏反演三维SRME方法预测的多次波.与第2道对比可知，稀疏反演方法预测多次波在能量级别与走时上和传统3D SRME方法基本一致，但是相位出现了反转现象，需要应用公式(20)进行稳态相位校正，图 4c第4道为相位校正后的结果，和第2道对比可知，稀疏反演多次波预测结果经相位校正后和传统全数据SRME方法预测得到的多次波无明显差别，相位校正后的数据体就是多次波自适应减的输入数据.

图 4

Fig. 4

图 4 稀疏反演三维SRME方法预测多次波单道数据分析 (a)Crossline多次波贡献道集(作为稀疏反演的输入数据)；(b)稀疏反演的模型空间域剖面；(c)原始单道数据(1)、密集道集叠加(2)、模型道叠加(3)与相位校正道(4)对比. Fig. 4 Single trace analysis of multiple prediction by three-dimensional SRME using sparse inversion (a)The Crossline multiple contributions used as input for the sparse inverse processing;(b)Sparse inversion data in model space domain;(c)The single trace of original data(trace1)，stack of dense gathers result(trace2)，stack of model gathers result(trace3)，and the phase correction trace(trace4).

在理论模型验证的基础上，本文也将研究形成的稀疏反演三维表面多次波预测与压制算法应用于中国某海洋探区的实际地震资料.该数据有6条Inline线(拖缆)，每条拖缆有144个检波器，线间距为100 m，道间距为12.5 m，炮间距为50 m，炮线和检波线平行，时间方向采样间隔为4 ms.根据稀疏反演多次波预测对地震数据的要求，对Inline线方向数据进行了规则化处理，规则化后炮间距也为12.5 m，计算所用炮数为204炮.规则化后数据Inline测线方向道间距为12.5 m，而Crossline方向的道间距为100 m，由于横测线方向间距过大，传统的3D SRME方法是无法预测准确的多次波数据，而且会带来严重的假频干扰.

应用稀疏反演三维表面多次波压制方法对主线规则化后的海洋地震数据进行了多次波预测和压制处理，为对比分析计算效果，文中也应用二维SRME方法对该实际数据进行了处理.图 5a为三维海洋数据体中的某条Inline线的含多次波的原始共偏移距剖面，可以看出该数据含有较强的海底一阶、二阶甚至高阶的表面多次波.图 5b为2D SRME方法预测的表面多次波，图 5c为基于稀疏反演三维方法预测的表面多次波.对比可知，三维稀疏反演方法预测的表面多次波的同相轴连续性更好，与原始多次波数据在能量对比关系上具有更好的一致性，也就是说可更好的与原始剖面的多次波同相轴进行拟合，为多次波自适应相减的有效处理提供了数据基础.

图 5

Fig. 5

图 5 三维实际海洋数据表面多次波压制 (a)原始地震数据共偏移距剖面；(b)2D SRME方法预测多次波；(c)稀疏反演3D SRME方法预测多次波； (d)2D SRME方法多次波压制；(e)稀疏反演3D SRME方法多次波压制. Fig. 5 3D marine dataset for multiple suppression (a)The common offset section of original seismic data;(b)The multiple predicted by 2D SRME;(c)The multiple predicted by sparse inversion 3D SRME;(d)The demultiple section by 2D SRME;(e)The demultiple section by sparse inversion 3D SRME.

针对两种方法的预测结果，文中采用前述理论模型测试的多次波压制算法进行处理，2D SRME及稀疏反演3D SRME多次波压制结果分别如图 5d和5e所示.对图中所示的双程走时约为0.55 s位置处箭头所示同相轴分析表明，2D SRME方法压制多次波后，其海底一阶表面多次波同相轴的能量仍很强，而3D稀疏反演的SRME方法中多次波能量得到很好的压制，同样，在双程走时约为1.5 s和1.8 s附近的黑色箭头也显示3D稀疏反演SRME方法多次波压制效果明显优于2D算法.双程走时约为0.8 s处箭头所示同相轴为一有效波同相轴，这一点从原始共偏移距剖面和预测的多次波剖面对比可以判断出，由于2D SRME预测的多次波存在振幅和相位的误差，使得自适应多次波匹配相减无法有效实现，有效波的振幅受到一定程度的损伤，与2D SRME方法相比，基于全三维假设的3D稀疏反演的SRME方法可较为精确的预测出多次波走时，进而可与原始数据中表面多次波有效拟合，可实现保幅的表面多次波压制处理.

相比于二维方法压制三维地震数据的多次波，三维SRME方法考虑了地震波在地下介质中传播机理和实际情况，更为适合复杂构造海洋地震资料的多次波压制处理，但是该算法对数据的采集要求很高.实际采集地震数据的横测线方向通常较为稀疏，为提高计算效率、避免大数据体的规则化处理以及规则化数据的存储，本文基于贡献道集的概念，将三维表面多次波贡献道集求取分两步实现，也就是在完成纵测线方向的多次波贡献积分后，应用高分辨率抛物Radon变换稀疏反演代替横测线的多次波道集求和，进而完成三维表面多次波的预测.应用CPU/GPU并行加速技术提高三维地震资料多次波预测和压制的计算效率，可大幅缩短海洋地震资料预处理周期.对于稀疏反演的三维多次波预测数据和全三维全数据预测的多次波存在的相位差异，文中应用稳相原理近似完成预测多次波的相位校正.理论模拟和实际海洋勘探三维数据的试算表明，稀疏反演三维表面多次波压制方法具有计算精度高、适用范围广和易于实现的特点.