基于波动方程三维表面多次波预测方法研究

石 颖，王维红，李 莹，井洪亮

1东北石油大学地球科学学院，大庆 163318

2大庆油田有限责任公司勘探开发研究院，大庆 163712

1 引 言

海洋地震勘探资料中表面多次波［1］的存在将严重降低地下构造的成像精度，因此，在海域反射地震数据处理中，表面多次波预测及压制［2］工作至关重要.随着勘探目标的日趋复杂化，复杂构造勘探对地球物理技术提出了新的挑战.在这些地区，受假设条件的限制，基于滤波的多次波压制技术［3－5］处理效果不明显，故而，多次波预测通常要利用数据驱动［6］的基于波动方程的方法［7－8］.近年来，基于波动方程的反馈迭代表面多次波压制（SRME）方法［9－11］凸显优势，该方法在无需地下介质任何信息的情况下，可同时预测所有阶表面多次波［12］，选择合理的自适应相减算法［13］，可有效压制复杂介质勘探地震资料的表面多次波.Berkhout和Verschuur首先基于反馈模型理论［14］提出2D级数展开压制表面多次波算法［15］，随后，两人又发展提升了该理论，在无需对多个预测项求和的情况下，基于反馈模型理论提出迭代压制表面多次波算法［16－17］，在某种程度上提高了计算的精度.基于这个思想，Kelamis，Wang，Van Groenestjin等［18－22］先后对基于波动方程的反馈迭代方法的理论和应用进行了有益的完善和改进，因此，其2D算法发展较快，在油气地震勘探资料处理，尤其是海洋地震资料处理中发挥着重要的作用.

基于波动方程的SRME技术虽可有效地压制多次波，但是，目前仍主要局限于2D［23］，或通过诸多的假设、振幅匹配以及乘经验因子等将三维数据利用二维算法近似计算.地下介质的真实结构为三维，因此不同地层之间的分界面具有三维结构，这表明在分界面上产生的一次波和多次波具有三维的传播和反射效应.因此，利用2D近似算法压制表面多次波，有悖于地下介质的真实情况，无法获得理想的压制效果.当考虑一阶表面多次波时，自由表面的反射点可能位于自由表面的任何位置，与产生多次波反射面的形状有关.考虑相对较为简单的情况，假定反射面在cross－line方向具有倾角，如图1所示，可观察到表面的反射点M并没有位于连接炮点和检波点的直线上，而是沿着倾角方向在横向上有移动.如果应用2D多次波预测方法，就是假定了所有的传播路径都发生在包含震源和检波点的垂直平面ABCD内，因此，也只能考虑in－line方向的震源和检波器组合.所以，要获得高精度的多次波预测和压制效果，需要在全3D意义上研究基于波动方程的SRME表面多次波压制方法.Van Dedem［24］分析了进行3D表面多次波压制的必要性，Levin［25］基于近似NMO的方法实现零偏移距的三维SRME，Van Dedem和Verschuur［26］利用稀疏反演的方法预测多次波，Kleemeyer，Lin，Matson，Moore也从不同角度对三维SRME方法进行过研究［27－30］，Baumstein［31］提出基于DMO实现三维SRME方法，利用DMO技术产生三维SRME的输入数据.Dragoset［32］提出，在所有的方法中，三维SRME方法是最有可能广泛应用于压制表面多次波的方法，理论上，3D算法避免了2D算法的所有假设，即使处理含噪音和不完整采样的地震数据，在噪音和缺失数据不是很严重的情况下，通常也会获得较2D算法更精确的结果.

图1 复杂介质中一阶表面多次波的射线路径Fig.1 Ray path of first order surface－related multiple from complicated medium

针对地震资料处理中2D表面多次波预测算法的局限性，文中提出全3D的表面多次波预测算法.首先基于3D多次波预测的基本理论，给出用于3D数据表面多次波预测的褶积矩阵表示方法；然后设计含表面多次波的3D模型，利用三维SRME方法进行测试；比较分析二维和三维方法的处理效果，阐明了三维SRME方法压制表面多次波的优越性.针对全三维SRME方法的计算量问题，文中引入了GPU和CPU协同并行加速计算技术［33－37］预测表面多次波，算法中涉及的密集型运算由GPU执行，主要是利用三维数据体形成的数据矩阵预测表面多次波.

2 3D多次波预测的基本理论

SRME为全数据驱动的表面多次波压制方法，其直接利用隐含于叠前地震数据中的信息预测多次波，而无需地下介质的任何其他信息，并且可同时预测出所有阶的表面多次波.已知与检波点xr有关的共接收点道集P（xr，xk，t）和与震源xs有关的共炮点道集P（xk，xs，t），利用它们的褶积和可预测多次波道m（xr，xs，t）.

在空间频率域（x，ω），基于波动方程的2D表面多次波预测公式可表示为

式中，P（xr，xk，ω）表示原始数据共接收点道集的1D傅里叶变换，P（xk，xs，ω）表示共炮点道集的1D傅里叶变换，xk表示下行反射点，X表示采集区域，MCG表示震源位于xs和检波点位于xr的多次波贡献道集，对所有的下行反射点xk进行计算.

实际上，地下介质是3D介质，而2D多次波预测算法忽略了地下空间的3D效应，因此，2D方法通常不能准确地预测表面多次波，而自适应匹配相减算法也仅能在有限的程度上改善这种预测的不匹配性，三维SRME方法考虑了地下介质的3D效应，是一种极为有效的多次波压制方法，能够压制大多数甚至是最复杂的3D表面多次波［38］.

全3D多次波预测算法描述了地震数据空间2D褶积过程.将2D多次波预测理论推广到全3D的情况，即

图2 观测系统示意图（a）共炮点观测系统；（b）共接收点观测系统.Fig.2 Geometry schematic diagram（a）Common shot geometry；（b）Common receiver geometry.

式中：M是预测的多次波，P是含多次波的总波场.利用地震数据与其本身在时空域的褶积，可实现多次波预测，由于SRME方法无需对地下介质做任何假设，是一种完全数据驱动的方法，因此该方法不仅要求密集采样的地下波场数据，也要求在x和y方向震源和检波器的采样密度相同，就是说在地表的每个位置都放置震源和检波器，使得炮间距等于道间距.

对密集采样的多次波贡献道集来说，执行式（2）等号右边的内部（第一个）求和时，意味着在每个cross－line位置，对所有的in－line位置数据求和，结果得到每个cross－line位置的一个多次波贡献道.执行式（2）等号右边的外部（第二个）求和，将所有cross－line位置数据相加，可获得最终的多次波预测道.

设计3D观测系统，在某震源位置放炮，在地表所有的接收点位置接收，产生共炮点道集，如图2a所示，其中圈示意检波点位置，星示意震源位置.移动震源的位置，可产生任意网格点的共炮点道集，见公式（2）中等式右边第二项.同样可确定共接收点道集，在每个检波点位置都放置震源，在与共炮点记录同样的网格上，选择一接收点位置，可模拟出共接收点道集，如图2b所示，移动接收点的位置，可产生任意网格点的共接收点道集.同时，密集网格上的每一点都将作为二次震源，激发能被所选择的接收点接收到的地震波，共接收点道集见公式（2）中等式右边第一项.将共炮点道集和共检波点道集逐道进行褶积运算，可得到密集采样的多次波贡献道集.多次波贡献道集可由一个体来表示，即地震道是in－line和cross－line坐标的函数.

假定3D采集是在x－y平面上理想的矩形网格区域上进行的，所有的网格点都是接收点，并且每一个网格点做震源位置一次.假定在x方向放置Nx个检波器，在y方向放置Ny个检波器，则每炮数据有NxNy个检波器接收，产生NxNy道地震记录，假定在此矩形区域的NxNy个位置都放炮，每个网格点做一次震源位置，每炮在NxNy个检波器上接收到的信号被记录下来，结果会得到（NxNy）2道数据.同2D方法一样，所有道变换到频域后，可在每个频率分量上构建数据矩阵，每个数据矩阵将包含单频分量上的全3D数据采集信息，如图3所示，在每个震源位置都做这样的重复，结果得到大小为NxNy的方阵.用这样的矩阵表示，矩阵的每列为3D共炮点道集，每行为3D共接收点道集，由此，可利用数据矩阵实现x和y方向2D空间褶积，进而预测3D地震勘探数据的表面多次波.

图3 3D数据预测表面多次波的矩阵表示Fig.3 3Ddata matrix representation for surface－related multiple prediction

3 理论模型应用实例

为了验证文中方法用于多次波预测的有效性，采用正演的3D模型数据进行试算.在较为理想的x和y方向都密集采样的3D数据体上测试前述的全3D多次波预测算法，3D地下介质的速度模型如图4所示，该模型包含两个水平地下反射界面.在地表，间距是15m的2D网格用于模拟全3D地震勘探，在x方向和y方向分别放置56个检波器，可得到3136个检波器的2D排列.在3D数据的模拟中，每个网格点用做震源位置一次，产生3136炮地震记录，每炮3136道地震数据.同2D算法，利用3D算法压制多次波仍需要两个过程，即多次波预测和自适应相减.图5a为理论模拟的含一阶表面多次波的地震炮记录，炮点位于第28条inline线，第28条crossline线的交点位置，共3136道记录，图5b为全3DSRME算法预测的多次波，图5c为借助于自适应相减算法计算的多次波压制结果，文中采用基于高频重建的自适应相减方法计算［39］.对三维数据，输入采样的总数等于所有维采样数的乘积，也因此增加了较大的计算成本［40］，本文采用GPU加速计算，以提高表面多次波预测的计算效率.

图4 3D水平层状介质表面多次波正演数据的速度模型Fig.4 Velocity model of forward modeling data with surface－related multiple for 3Dhorizontal layer medium

图5 某炮数据56条测线的地震记录（a）某点放炮模拟的地震数据；（b）3D算法预测的多次波；（c）多次波压制后的结果.Fig.5 Seismic record with 56lines from one shot（a）Simulated one shot record；（b）Predicted multiple based on the proposed 3Dpredictive approach；（c）Result after multiple suppression.

前已述及，基于波动方程的全3D表面多次波预测算法需要在每个频率分量上完成地震数据的空间褶积运算，由3D数据形成庞大的预测矩阵，若采用常规的CPU对数据矩阵执行串行预测计算，其较低的计算效率将在很大程度上阻碍该方法的发展和应用，文中引入GPU加速计算技术，用以完成大规模数据矩阵的相关运算.在矩阵运算中，由矩阵元素的乘法和加法所实现的矩阵乘法计算量巨大，因此，文中对矩阵进行棋盘划分，将预测矩阵划分为若干个子方阵，利用共享存储器访问优化实现矩阵乘法以及线程间通信.访问共享存储器的延迟远小于全局存储器，先将需要计算的数组数据读到共享存储器中，再利用共享存储器中的数据进行计算，具体实现时，首先从全局存储器将两个对应的子方阵载入共享存储器，一个线程负责加载一个元素，然后每个线程负责计算乘积中的一个元素.各线程将所有乘积结果汇总到寄存器中，再将结果写入全局存储器.这种棋盘划分方式充分利用了高速的共享存储器，同时节省了大量的全局存储器带宽.本文采用的计算平台为NVIDIA GeForce GTX 550Ti，计算能力为2.1，GPU主频为1.8GHz，而另一计算环境为主频3.30GHz的双核CPU，计算文中3D模型数据，CPU串行计算耗时158794s，GPU和CPU协同并行计算耗时963s，较常规的CPU串行计算，其效率可提高约165倍.

为更清楚地分析本文三维多次波预测与压制效果，图6给出了一条主线排列的单炮头记录（Inline1），该记录的炮点位于测线Inline1和测线Crossline26的交点.图6a为含一阶表面多次波的原始数据，图6b为应用本文给出的3D多次波预测算法预测的多次波炮记录，图6c为多次波自适应匹配相减的成果炮记录，分析可知，最小双程旅行时分别位于1.06s和1.3s附近的多次波被有效压制，单炮记录的信噪比明显提高.

图6 某条测线的地震记录（a）某条测线含多次波的地震数据；（b）3D方法预测的多次波；（c）多次波压制结果.Fig.6 Seismic record from one line（a）Seismic data with multiples from one line；（b）Predicted multiple based on the proposed 3Dapproach；（c）Result after multiple suppression.

4 3D和2D多次波预测方法对比分析

当前绝大多数海洋地震资料预处理中，经常利用2D算法预测表面多次波，结合自适应相减算法进行多次波压制.然而，在许多情况下，对地下介质和观测系统所做的2D假设在某种程度上是无效的，使得预测的多次波出现了很大的误差，难以在自适应相减中有效地减去，3D表面多次波预测方法理论上符合地震波在地下介质中的传播规律，与2D算法相比，利用全3D算法预测多次波可获得更高的精度，从而获取更为理想的多次波自适应压制结果.

为了更好地对比分析2D和3D多次波预测和压制效果，将上述3D理论模型数据的一条测线（inline1）上的56炮数据分别利用2D和全3D算法进行计算，其中12炮数据（每4炮抽取1炮显示）的多次波预测和压制效果如图7所示.3D数据inline1测线含表面多次波的12炮理论模型数据如图7a所示；图7b为2D算法预测的表面多次波；图7c是经2D算法预测，并与图7a所示数据进行自适应相减后的多次波压制结果；图7d为3D算法预测的表面多次波；图7e是经3D算法预测，并与图7a所示数据进行自适应相减后的多次波压制结果.值得提出的是，文中对2D和3D预测数据均采用同一自适应匹配相减算法.对比分析可知，2D算法用于3D数据的多次波预测和压制，在某种程度上可以减弱多次波的能量，提高地震资料的信噪比，但炮记录上仍可见到较强的多次波残余，而3D算法则能够高精度地预测多次波，因此也可获得较为理想的压制效果.

为进一步对比2D和3D多次波预测方法，现抽取图6a所示的第40道地震记录，2D和3D多次波预测和压制效果的单道记录分别如图8和9所示.图8a和图9a为原始地震道记录，图8b和图9b分别为利用2D和3D预测算法得到的多次波压制结果，图8b的1.06s和1.3s所对应的同相轴有较大的多次波残余，而图9b中多次波在对应位置压制效果明显.对比发现，3D算法的预测精度更高，压制效果也更为理想，二维算法只在某种程度上减弱了多次波的能量，但仍有残余存在.

5 结 论

本文研究给出了基于波动方程的全3D表面多次波预测算法，其考虑了地震波在地下介质中传播时所有可能的反射路径，对任一频率分量，空间方向的二维褶积可实现多次波预测，文中也对利用2D和3D算法预测和压制多次波的效果进行了对比分析.通过理论模型的试算分析发现，基于波动方程的全3D算法可高精度地预测表面多次波，其预测效果明显优于2D算法，借助于有效的自适应相减方法，可获得较为理想的表面多次波压制效果.文中也利用GPU加速技术提高3D算法的计算效率，重点用于利用数据矩阵预测多次波的计算.就文中模型来说，GPU加速计算效率可提高165倍左右.

除了计算效率问题，地震数据空间采样密度不足以及较大的存储空间要求也制约着该方法的发展，文中利用GPU加速技术解决了全3D算法计算成本问题，随着数据重建预处理技术以及计算机硬件技术的不断发展和完善，3DSRME方法将有更大的发展空间.

（

）

［1］ 郭书娟，李振春，仝兆岐等.表层多次波成像方法技术研究.地球物理学进展，2012，27（6）：2570－2576.

Guo S J，Li Z C，Tong Z Q，et al.Method and technique for imaging of surface－related multiples.ProgressinGeophysics（in Chinese），2012，27（6）：2570－2576.

［2］ 赵昌垒，叶月明，姚根顺等.线性拉东域预测反褶积在海洋多次波去除中的应用.地球物理学进展，2013，28（2）：1026－1032.

Zhao C L，Ye Y M，Yao G S，et al.Prediction deconvolution in linear radon domain on the application of ocean multiples attenuation.ProgressinGeophys.（in Chinese），2013，28（2）：1026－1032.

［3］ 胡天跃.地震资料叠前去噪技术的现状与未来.地球物理学进展，2002，17（2）：218－223.

Hu T Y.The current situation and future of seismic data prestack noise attenuation techniques.ProgressinGeophys.（in Chinese），2002，17（2）：218－223.

［4］ 王维红，崔宝文，刘洪.表面多次波衰减的研究现状与进展.地球物理学进展，2007，22（1）：156－164.

Wang W H，Cui B W，Liu H.Research progress in surfacerelated multiple attenuation.ProgressinGeophys.（in Chinese），2007，22（1）：156－164.

［5］ 石颖，邢小林.表面多次波压制的研究进展：回顾与展望.地球物理学进展，2011，26（6）：2046－2054.

Shi Y，Xing X L.Investigation progress on surface－related multiple suppression：review and outlook.Progressin Geophys.（in Chinese），2011，26（6）：2046－2054.

［6］ 吴迪，杨长春.数据驱动型多次波衰减方法的研究.地球物理学进展，2008，23（1）：98－103.

Wu D，Yang C C.Research on a fully－data driven multiple attenuation approach.ProgressinGeophys.（in Chinese），2008，23（1）：98－103.

［7］ 赵保宗，孙永清，李学聪.基于波动方程的多次波压制方法应用研究.地球物理学进展，2010，25（1）：272－281.

Zhao B，Sun Y Q，Li X C.The progress of the multiple attenuation method based on the wave－equation.Progressin Geophys.（in Chinese），2010，25（1）：272－281.

［8］ 李鹏，刘伊克，常旭等.多次波问题的研究进展.地球物理学进展，2006，21（3）：888－897.

Li P，Liu Y K，Chang X，et al.Progress on the multiple problems.ProgressinGeophys.（in Chinese），2006，21（3）：888－897.

［9］ 陈小宏，刘华锋.预测多次波的逆散射级数方法与SRME方法及比较.地球物理学进展，2012，27（3）：1040－1050.

Chen X H，Liu H F.Comparison between inverse scattering series method and SRME method in free surface related multiple prediction.ProgressinGeophys.（in Chinese），2012，27（3）：1040－1050.

［10］ 石颖，井洪亮，李莹.反馈迭代法压制表面多次波效果分析.地球物理学进展，2012，27（4）：1493－1500.

Shi Y，Jing H L，Li Y.Surface－related multiple suppression effect analysis by feedback iteration approach.Progressin Geophys.（in Chinese），2012，27（4）：1493－1500.

［11］ 肖二莲，陈瑜，万欢等.SRME多次波衰减方法在海洋地震资料中的应用.地球物理学进展，2010，25（3）：1057－1064.

Xiao E L，Chen Y，Wan H，et al.Surface－related multiple elimination on marine seismic data.ProgressinGeophys.（in Chinese），2010，25（3）：1057－1064.

［12］ Donno D，Chauris H，Noble M.Curvelet－based multiple prediction.Geophysics，2010，75（6）：WB255－WB263.

［13］ 石颖，王建民，井洪亮等.多道自适应匹配滤波方法压制表面多次波.地球物理学进展，2013，28（2）：785－792.

Shi Y，Wang J M，Jing H L，et al.Suppressing surfacerelated multiple by multi－trace adaptive matching filter approach.ProgressinGeophys.（in Chinese），2013，28（2）：785－792.

［14］ Berkhout A J.Seismic Migration，Practical aspects Pt.B：Imaging of Acoustic Energy by Wave Field Extrapolation.Amsterdam：Elsevier Science Publishing Company，Inc，1982.20－43.

［15］ Verschuur D J，Berkhout A J，Wapenaar C P A.Adaptive surface－related multiple elimination.Geophysics，1992，57（9）：1166－1177.

［16］ Berkhout A J，Verschuur D J.Estimation of multiple scattering by iterative inversion，partⅠ：Theoretical considerations.Geophysics，1997，62（5）：1586－1595.

［17］ Verschuur D J，Berkhout A J.Estimation of multiple scattering by iterative inversion，partⅡ：Practical aspects and examples.Geophysics，1997，62（5）：1596－1611.

［18］ Kelamis P G，Verschuur D J.Surface－related multiple elimination on land seismic data—Strategies via case studies.Geophysics，2000，65（3）：719－734.

［19］ Wang Y H.Multiple prediction through inversion：a fully data－driven concept for surface－related multiple attenuation.Geophysics，2004，69（2）：547－553.

［20］ Wang Y H.Multiple prediction through inversion：Theoretical advancements and real data application.Geophysics，2007，72（2）：V33－V39.

［21］ Van Groenestjin G J A，Verschuur D J.Estimating primaries by sparse inversion and application to near－offset data reconstruction.Geophysics，2009，74（3）：A23－A28.

［22］ Van Groenestjin G J A，Verschuur D J.Estimation of primaries and near－offset reconstruction by sparse inversion：Marine data applications.Geophysics，2009，74（6）：R119－R128.

［23］ Lin D C，Young J，Lin W J，et al.3DSRME prediction and subtraction practice for better imaging.75th Ann.Internat.Mtg.，Soc.of Expi.Geophys.，Expanded Abstracts，2005.2088－2091.

［24］ Van Dedem E J，Verschuur D J.Analysis of surface－related multiples in 3－D media.67th Annual International Meeting，SEG，Expanded Abstracts，1997.1180－1183.

［25］ Levin S A.Prestack poststack 3Dmultiple prediction.72nd Annual International Meeting，SEG，Expanded Abstracts，2002.2110－2113.

［26］ Van Dedem E J，Verschuur D J.3Dsurface－related multiple prediction using sparse inversion：Experience with field data.72nd Annual International Meeting，SEG，Expanded Abstracts，2002.2094－2097.

［27］ Kleemeyer G，Pettersson S E，Eppenga R，et al.It’s a magic：Industry first 3Dsurface multiple elimination and prestack depth migration on Ormen Lange.65th Annual International Meeting，EAGE，Extended Abstracts，2003.

［28］ Lin D C，Young J，Huang Y，et al.3DSRME application in the Gulf of Mexico.74th Annual International Meeting，SEG，Expanded Abstracts，2004.1257－1260.

［29］ Matson K H，Michell S，Abma R，et al.Advanced subsalt imaging and 3Dsurface multiple attenuation in Atlantis：A case study.74th Annual International Meeting，SEG，Expanded Abstracts，2004.1269－1272.

［30］ Moore I，Dragoset W H.Practical，3Dsurface－related multiple prediction（SMP）.74th Annual International Meeting，SEG，Expanded Abstracts，2004.1249－1252.

［31］ Baumstein A，Hadidi M T.3Dsurface－related multiple elimination：Data reconstruction and application to field data.Geophysics，2006，71（3）：E25－E33.

［32］ Dragoset B，Verschuur E，Moore L，et al.A perspective on 3Dsurface－related multiple elimination.Geophysics，2010，75（5）：75A245－75A261.

［33］ Zhang J H，Wang S Q，Yao Z X.Accelerating 3DFourier migration with graphics processing units.Geophysics，2009，74（6）：WCA129－WCA139.

［34］ Clapp R G，Fu H H，Lindtjorn O.Selecting the right hardware for reverse time migration.TheLeadingEdge，2010，29（1）：48－58.

［35］ 李博，刘国峰，刘洪.地震叠前时间偏移的一种图形处理器提速实现方法.地球物理学报，2009，52（1）：245－252.

Li B，Liu G F，Liu H.A method of using GPU to accelerate seismic pre－stack time migration.ChineseJ.Geophys.（in Chinese），2009，52（1）：245－252.

［36］ 刘国峰，刘洪，李博等.山地地震资料叠前时间偏移方法及其GPU实现.地球物理学报，2009，52（12）：3101－3108.

Liu G F，Liu H，Li B，et al.Method of prestack time migration of seismic data of mountainous regions and its GPU implementation.ChineseJ.Geophys.（in Chinese），2009，52（12）：3101－3108.

［37］ 刘红伟，李博，刘洪等.地震叠前逆时偏移高阶有限差分算法及GPU实现.地球物理学报，2010，53（7）：1725－1733.

Liu H W，Li B，Liu H，et al.The algorithm of high order finite difference pre－stack reverse time migration and GPU implemenation.ChineseJ.Geophys.（in Chinese），2010，53（7）：1725－1733.

［38］ Baumstein A，Hadidi M T，Hinkley D L，et al.A practical procedure for application of 3DSRME to conventional marine data.TheLeadingEdge，2005，24（3）：254－258.

［39］ 石颖，刘洪，邹振.基于波动方程表面多次波预测与自适应相减方法研究.地球物理学报，2010，53（7）：1716－1724.

Shi Y，Liu H，Zou Z.Surface－related multiples prediction based on wave equation and adaptive subtraction investigation.ChineseJ.Geophys.（in Chinese），2010，53（7）：1716－1724.

［40］ Xu S，Zhang Y，LambaréG.Antileakage Fourier transform for seismic data regularization in higher dimensions.Geophysics，2010，75（6）：WB113－WB120.