田东升，王云专，李义鹏，石 颖，柯 璇，李婷婷，刘淑芬

（1.东北石油大学 地球科学学院，黑龙江 大庆 163318； 2.北华航天工业学院 电子工程系，河北 廊坊 065000）

单程和双程波动方程叠前深度偏移方法

田东升1，王云专1，李义鹏2，石 颖1，柯 璇1，李婷婷1，刘淑芬1

（1.东北石油大学 地球科学学院，黑龙江 大庆 163318； 2.北华航天工业学院 电子工程系，河北 廊坊 065000）

叠前深度偏移是获得地下构造映像的有效手段，而基于波动方程的叠前深度偏移方法对速度横向变化剧烈的地层有更好的适应性.分析基于单程波方程的相移法、相移加插值法、频率—空间域有限差分法、傅里叶有限差分法和基于双程波方程的逆时偏移方法，借助于地堑模型与盐丘模型，测试5种逆时偏移方法成像复杂构造的精度和适应性.结果表明，基于波动方程的叠前深度偏移方法可实现横向变速地下构造成像，相比于基于双程波方程的逆时偏移方法，单程波方程方法对垂直断层等高陡倾角构造成像有局限性；逆时偏移方法对垂直断层、盐丘下边界等复杂构造可以清晰成像，辅以精确的地层速度，逆时偏移方法在地震资料成像领域中有广阔的发展和应用前景.

叠前深度偏移；相移法；相移加插值；频率—空间域有限差分；傅里叶有限差分；逆时偏移

0 引言

随着油气勘探目标日趋复杂，具体表现在断块小、倾角陡、纵横向速度变化剧烈等方面，常规偏移方法很难达到地震资料处理要求，加之巨大的数据处理量，亟需研究高精度、高效率的偏移算法，为精细地质构造解释及储层识别提供重要依据.叠前深度偏移是获得精确地下构造的有效途径.在数学解法上，叠前深度偏移分为两类：基于波动方程积分解的射线偏移和基于波动方程微分解的波动方程叠前深度偏移.波动方程叠前深度偏移较射线偏移在处理横向变速问题上具有更强的适应性，不存在射线偏移法中成像点的多值走时问题.其中波动方程叠前深度偏移分为单程波法和双程波法.

20世纪70年代，Claerbout J F[1]提出应用波动方程进行偏移，采用有限差分法求解得到单程波动方程15°近似公式，该方法在主传播方向小范围内具有较好的效果，对宽方位地震波传播模拟并不理想，尤其对于地震数据中包含的水平和陡倾角反射信息的偏移效果不明显.Gazdag J最初提出的频率—波数域相移法[2]仅适用于对地下横向速度不发生变化的介质进行成像；进而提出的相移加插值方法[3]可以适应存在横向速度变化的介质，但需要频繁计算参考波场，效率较低[4].随后，人们将有限差分法与傅里叶法结合研究，提出混合偏移方法[5-6]，对横向速度变化强烈的介质成像效果良好.王玉学等[7]推导上行波方程的两种高阶近似表达式.冯凤萍等[8]研究加吸收层的三维45°上行波方程的隐式差分格式.然而，所有单程波偏移方法存在局限性，即单程波算子在成像大角度传播的波时将发生相位改变和振幅减弱的现象，无法对陡倾角进行成像；另外，单程波法也无法成像回转波.

采用双程波动方程进行偏移能很好地适应剧烈的横向速度变化，可以有效解决复杂地质体成像问题，最典型的方法就是逆时偏移方法.该方法最早由Whitmore D N等[9]提出，最初应用于处理叠后资料.逆时偏移方法简单、易于实现，不对波动方程做任何近似，从而不存在倾角的限制，可以对透射波、多次波、绕射波等进行成像.近些年，随着计算机硬件技术的迅猛发展及勘探要求的日益提高，逆时偏移方法的研究也从叠后走向叠前[10-11]，从二维走向三维.GPU加速计算技术的引入[12]和噪音压制策略的研究[13]推动逆时偏移技术的发展.同时，如全波形反演的精准速度建模方法的研究[14]也加速波动方程叠前深度偏移的研究进程.

笔者首先阐述单程和双程波动方程叠前深度偏移方法的基本原理，分析不同方法优缺点，通过地堑模型和二维盐丘模型进行成像测试，分析不同方法成像复杂构造的精度和适应性，为高精度叠前深度偏移方法的工业化应用提供依据.

1 单程波方程叠前深度偏移方法

1.1 相移和相移加插值法

在纵波勘探中，假设地下介质为均匀各向同性介质，并且介质密度恒定，可以用声波方程描述地震波的传播，二维形式为

式中：v为地下介质速度；W 为t时刻（x，z）位置处的波场值；t为时间；x，z为空间方向.式（1）在频率—波数域的解析解为

式中：kx为水平波数；kz为垂直波数；d z为深度延拓步长；w为频率.

通过标量波动方程的频散关系得到上、下行波分解的形式表达式，将它代入式（2）并求偏导数，得到频率—波数域单程波方程为

式（3）为相移法的基本原理.相移法的优点在于稳定性，对网格间距没有要求，但是仅能对横向速度无变化的介质进行偏移；因此也决定它无法对复杂地质体准确成像.

为了能更好地适应横向速度变化，引入相移加插值法，将式（2）分解成2个独立的公式：

式中：v′为v 的近似，v′≠v（x，z）.

利用式（4）对每道进行时移；再根据式（5）使用v（x，z）的最小值和最大值进行相移；最后对两个相移波场的模数和相位角进行线性插值，得到最终的波场信息.相移加插值法较相移法有很大的改进，可对横向速度变化缓慢的介质进行成像[4].

1.2 频率—空间域有限差分法

对式（3）二阶近似后分离成两部分，可得频率—空间域有限差分法表达式：

式（6）为绕射项方程，可以通过有限差分求解，实现对绕射波的收敛.式（7）可看作折射项方程，在频率—空间域表示为

式中：Dxx为x方向二阶导数；Dz为z方向一阶导数.

通过有限差分法或相移法计算式（8），可以完成横向速度变化的时差校正[15].对偏移倾角不大的地质构造，该方法成像效果较理想.

1.3 混合偏移法

混合偏移法是将有限差分法与傅里叶法相结合，主要包括裂步傅里叶法和傅里叶有限差分法.裂步傅里叶法基于扰动理论，将速度场分离成背景速度项和扰动项.该方法首先在频率—波数域以背景速度延拓，而后在频率—空间域利用表示速度横向变化的扰动项进行相移；在成像横向速度变化强烈的介质时，成像效果不理想.

傅里叶有限差分法以裂步傅里叶法为基础，在波场延拓过程中，引入自适应有限差分算子.由于在选取偏移速度时，采用介质速度v和参考速度vr偏移得到的结果不一致，根据二者偏移误差d，可得傅里叶有限差分法表达式：

第Ⅰ部分是在频率—波数域进行的相位移算子；第Ⅱ部分是裂步傅里叶法在频率—空间域进行的相位移算子；第Ⅲ部分为有限差分算子.当地下介质为层状介质时，式（9）只保留第Ⅰ部分，该算法变为相移法；当地下介质速度横向变化剧烈时，该算法变为有限差分法[16].因此，傅里叶有限差分法是裂步傅里叶法和频率—空间域有限差分法的混合偏移方法，兼顾二者各自的优点.

2 双程波方程叠前深度偏移方法

相对于单程波方程偏移算法，通过求解双程波方程得到波场传播信息的典型方法是叠前深度逆时偏移方法，实现步骤主要包括：（1）震源处波场沿时间轴正向延拓，保存各个时刻波场信息；（2）检波点处波场沿时间轴反向延拓；（3）对同一时刻的两个波场进行成像，完成单炮逆时深度偏移.对各炮的成像结果进行叠加，得到叠前逆时偏移结果.

在逆时偏移计算中，对式（1）进行高阶有限差分，空间方向2N阶差分格式[17]为

时间方向差分格式为

逆时偏移的成像条件主要有激发时刻成像、互相关成像和除法成像.目前最常用的成像条件是互相关成像条件，即

式中：s（x，z，t）、r（x，z，t）分别为震源波场与检波点波场.

采用Robert G Clapp提出的随机边界条件[18]，无需存储波场正传过程中所有时刻的波场，节省大量的存储空间；利用GPU/CPU协同并行技术加速逆时偏移计算[19]，提高算法的计算效率；由于逆时偏移成像结果中含有大量低频噪音，通过拉普拉斯算子法对噪音进行压制[20-21].不同波动方程叠前深度偏移方法优缺点见表1.

表1 波动方程叠前深度偏移方法优缺点Table 1 Comparison of different wave equation prestack depth migration approaches

3 计算精度分析

为测试单程和双程波动方程叠前深度偏移方法对复杂构造的偏移效果，分别成像地堑模型和二维盐丘模型，分析不同算法的成像精度和适应性.

3.1 地堑模型

地堑模型网格点数为200×100，网格大小为5 m×5 m，上层介质速度为2 000 m/s，下层速度为3 000 m/s，激发震源采用Ricker子波，主频为40 Hz，震源与检波器始于地面最左端，炮间距20 m，共50炮，每炮200道接收，道间距为5 m.

地堑速度模型见图1（a）；采用频率—空间域有限差分算法得到的偏移结果见图1（b）、（c），其中图1（c）将正演中的棱柱波切除.由图1（b）、（c）可知，棱柱波对成像结果存在影响（图中方框位置）.另外，成像剖面中上层存在较大频散，是由算法本身造成的（图中箭头指示位置）.采用傅里叶有限差分算法得到的偏移结果见图1（d）、（e），其中图1（e）将正演中的棱柱波切除.由图1（d）、（e）可知，傅里叶有限差分法对棱柱波不能准确归位，因此对垂直断层成像效果不理想.采用裂步傅里叶法得到的偏移结果见图1（f）.由图1（f）可知，下部同相轴不清晰，边界无法识别.基于单程波方程的3种算法对水平分层可有效成像，但无法成像垂直断层.

采用逆时偏移方法经拉普拉斯算子法压制低频噪音后得到的偏移结果见图1（g）.由图1（g）可知，成像剖面清晰，水平分层明显，垂直断层得到很好归位（图中椭圆位置）.相对于单程波算法，逆时偏移方法可以对高陡倾角准确成像，即对横向速度变化剧烈的介质成像效果较好.

图1 不同叠前深度偏移算法成像地堑模型Fig.1 Graben model imaging result by different prestack depth migration approaches

3.2 二维盐丘模型

二维盐丘模型网格点数为150×649，网格大小为24.384 m×24.384 m；震源采用Ricker子波，主频为18 Hz；共325炮，每炮176道接收，炮间距为48.768 m.

二维盐丘速度模型见图2（a），采用有限差分算法得到的偏移结果见图2（b）.由图2（b）可知，盐丘模型整体成像不清晰，上部断层不明显，上边界较模糊，下边界几乎无法识别.采用相移加插值法得到的偏移结果见图2（c），采用傅里叶有限差分算法得到的偏移结果见图2（d）.由图2（c）、（d）可知，两者对横向变速的地质体成像有明显改善，上部断层成像不清晰，下部边界不准确.总体上，在单程波法偏移中，相移加插值算法成像精度较高，但耗时最长.采用叠前逆时偏移方法，经拉普拉斯算子法压制低频噪音后得到的成像结果见图2（e）.由图2（e）可知，与单程波偏移结果相比，逆时偏移得到的二维盐丘剖面成像清晰，上部断层得到很好归位（图中方框位置），盐丘下边界同相轴明显（图中椭圆位置），下部陡倾角得到清晰成像，成像剖面质量良好.

图2 不同叠前深度偏移算法成像结果Fig.2 Imaging result by different prestack depth migration approaches

4 结论

（1）基于波动方程的叠前深度偏移方法不受高频近似和多值走时的影响，对速度横向变化剧烈的地层有很好的适应性.

（2）单程波方程方法在成像大角度传播的地震波场时，存在相位改变和振幅削弱的问题；另外，因无法成像回转波，导致难以对陡倾角构造成像.

（3）基于双程波的逆时偏移法可有效成像各种地震波，利用精准的地层速度，对横向变速地层和高陡构造地层均能高精度成像，但是存储尤其是3D数据存储问题还有待于解决.

[1] Claerbout J F.Toward a unified theory of reflector mapping[J].Geophysics，1971，36（3）：467-481.

[2] Gazdag J.Wave equation migration with the phase-shift method[J].Geophysics，1978，43（7）：1342-1351.

[3] Gazdag J，Sguazzero P.Migration of seismic data by phase shift plus interpolation[J].Geophysics，1984，49（2）：124-131.

[4] 张钋，李幼铭，刘洪.几类叠前深度偏移方法的研究现状[J].地球物理学进展，2000，15（2）：30-39.

Zhang Po，Li Youming，Liu Hong.The situation of several prestack depth migration methods[J].Progress in Geophysics，2000，15（2）：30-39.

[5] Stoffa P L，Fokkema J T.Split-step Fourier migration[J].Geophysics，1990，55：410-421.

[6] Ristow D，Ruhl T.Fourier finite-difference migration[J].Geophysics，1994，59（12）：1882-1893.

[7] 王玉学，周瑞芬，张廷全.三维上行波的高阶近似[J].大庆石油学院学报，2003，27（2）：120-122.

Wang Yuxue，Zhou Ruifen，Zhang Tingquan.High order approximation of 3D one way wave equation[J].Journal of Daqing Petroleum Institute，2003，27（2）：120-122.

[8] 冯凤萍，周瑞芬，刘畅.加吸收层的三维45°上行波方程的隐式差分格式[J].大庆石油学院学报，2005，29（4）：98-100.

Feng Fengping，Zhou Ruifen，Liu Chang.Implicit difference scheme of 45 degree up-going wave equation with an absorbing layer[J].Journal of Daqing Petroleum Institute，2005，29（4）：98-100.

[9] Whitmore D N.Iterative depth imaging by back time propagation[C].53rd Annual International Meeting，SEG Expanded Abstracts，1983：382-385.

[10] Yoon K，Marfurt K J，Starr W.Challenges in reverse-time migration[C].74th Annual International Meeting，SEG Expanded Abstracts，2004，23（1）：1057-1060.

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

Liu Hongwei，Li Bo，Liu Hong，et al.The algorithm of high order finite difference prestack reverse time migration and GPU implementation[J].Chinese Journal of Geophysics，2010，53（7）：1725-1733.

[12] 李博，刘红伟，刘国峰，等.地震叠前逆时偏移算法的CPU/GPU实施对策[J].地球物理学报，2010，53（12）：2938-2943.

Li Bo，Liu Hongwei，Liu Guofeng，et al.Computational strategy of seismic pre-stack reverse time migration on CPU/GPU[J].Chinese Journal of Geophysics，2010，53（12）：2938-2943.

[13] Zhang Yu，Sun J.Practical issues of reverse time migration：True amplitude gathers，noise removal and harmonic-source encoding[J].First Break，2009，26，29-35.

[14] 杨午阳，王西文，雍学善，等.地震全波形反演方法综述[J].地球物理学进展，2013，28（2）：0766-0776.

Yang Wuyang，Wang Xiwen，Yong Xueshan，et al.The review of seismic full waveform inversion method[J].Progress in Geophysics，2013，28（2）：0766-0776.

[15] 马淑芳，李振春.波动方程叠前深度偏移方法综述[J].勘探地球物理进展，2007，30（7）：153-161.

Ma Shufang，Li Zhenchun.Review of wave equation prestack depth migration methods[J].Progress in Exploration Geophysics，2007，30（7）：153-161.

[16] 郑晓.起伏地表单程波叠前深度偏移方法研究[D].北京：中国地质大学，2011，34-43.

Zheng Xiao.Method study of rugged topography prestack one-way equation depth migration[D].Beijing：China Unversity of Geosciences，2011，34-43.

[17] 石颖，柯璇，田东升，等.复杂构造地震数据叠前逆时偏移方法[J].数学的实践与认识，2013，43（10）：206-213.

Shi Ying，Ke Xuan，Tian Dongsheng，et al.Seismic data prestack reverse time migraion approach on complicated construction[J].Mathematics in Practice and Theory，2013，43（10）：206-213.

[18] Robert G Clapp.Reverse time migration with random boundaries[C].79th Annual International Meeting，SEG Expanded Abstracts，2009：2809-2813.

[19] 石颖，陆加敏，柯璇，等.基于GPU并行加速的叠前逆时偏移方法[J].东北石油大学学报，2012，36（4）：111-115.

Shi Ying，Lu Jiamin，Ke Xuan，et al.Prestack reverse time migration based on GPU parallel accelerating algorithm[J].Journal of Northeast Petroleum University，2012，36（4）：111-115.

[20] 刘红伟，刘洪，邹振.地震叠前逆时偏移中的去噪与存储[J].地球物理学报，2010，53（9）：2171-2180.

Liu Hongwei，Liu Hong，Zou Zhen.The problems of denoise and storage in seismic reverse time migration[J].Chinese Journal of Geophysics，2010，53（9）：2171-2180.

[21] 石颖，王建民，田东升，等.应用于叠前逆时偏移的拉普拉斯算子去噪法[J].大庆石油地质与开发，2012，31（5）：154-157.

Shi Ying，Wang Jianmin，Tian Dongsheng，et al.Denoise approach by Laplace operator applied in prestack reverse time migration[J].Petroleum Geology and Oilfield Development in Daqing，2012，31（5）：154-157.

TE132.1

A

2095-4107（2014）04-0039-06

2013-02-12；

任志平

黑龙江省教育厅科学技术研究项目（12511025）

田东升（1989-），男，硕士研究生，主要从事地震资料处理方面的研究.通讯作者：石 颖，E-mail：shiyingdqpi@163.com

DOI 10.3969/j.issn.2095-4107.2014.04.006