2024年5月16日发(作者：褚振博)
声波方程逆时偏移中的无分裂PML吸收边界条件
王鹏飞;何兵寿
【摘 要】逆时偏移是当前地震资料处理的前沿技术,吸收边界条件是逆时偏移技术
的重要组成部分.目前,业界常用的基于波场分裂思路的完全匹配层PML(Perfectly
Matched Layer)吸收边界条件在改善边界吸收效果和改进偏移成像质量方面发挥
了巨大作用,但这种方法需要在边界处为分裂后的各个分量开辟额外的内存空间,且
需要分边角处理,增加了逆时偏移技术的内存负担和计算开销.为了减少内存负担和
提高计算效率,首先从双程声波方程出发,推导了声波方程的无分裂PML吸收边界
条件,然后给出该边界条件下波动方程逆时延拓的数值实现过程.理论分析和模型实
验结果得到:无分裂PML边界条件具有与分裂算法相同的边界吸收效果,并且基于
无分裂PML的逆时偏移算法效率更高,更便于程序代码的编写和GPU
(GraphicProcessing Unit)的并行实现.
【期刊名称】《工程地球物理学报》
【年(卷),期】2015(012)005
【总页数】8页(P583-590)
【关键词】逆时偏移;吸收边界条件;无分裂完全匹配层;计算效率
【作 者】王鹏飞;何兵寿
【作者单位】中国海洋大学海底科学与探测技术教育部重点实验室,山东青岛
266100;中国海洋大学海底科学与探测技术教育部重点实验室,山东青岛266100
【正文语种】中 文
【中图分类】P631.5
逆时偏移是当前地震资料成像处理领域的热点技术。目前，国内外在逆时偏移领域
的研究工作主要集中在以下几个方面：①波动方程逆时延拓算法[1,2]研究，包括
差分格式推导、边界伪反射压制和数值频散压制等；②成像准则和成像方法[3-8]
研究；③逆时偏移噪声压制方法[9-14]研究，主要研究波场延拓过程中由层间反射
导致的低频噪音的压制方法；④逆时偏移存储策略[15-17]研究，主要用于降低逆
时偏移的临时文件存储量和硬盘访问量，拓展逆时偏移处理的并行性并提高处理效
率；⑤偏移处理的并行算法[18-20]研究，包括CPU(Central Processing Unit)并
行、GPU并行以及CPU+GPU协同并行。
本文内容属于逆时延拓算法和逆时偏移并行处理范畴，对声波方程逆时偏移中常用
的PML吸收边界条件进行了改进，使之由分裂形式变为无分裂形式，避免了计算
区域的分割，降低了内存开销，提升了计算效率，且简化了代码编写工作。本文在
GPU计算平台CUDA(Compute Unified Device Architecture)上实现了无分裂
PML条件的逆时延拓。模型实验结果表明，无分裂PML边界条件具有与分裂算法
相同的边界吸收效果，且降低了内存消耗，减少了运算次数，更便于程序代码的编
写和GPU的并行实现。
逆时偏移可以通过以下三步实现：①炮点波场正向延拓；②接收点波场逆时延拓；
③炮检波场互相关成像。其中炮点波场的正向延拓方法与对应地震波方程的正演方
法相同，本文不赘述，只简要介绍接收点波场逆时延拓方法与波场成像方法。
2.1 接收点波场的逆时延拓
以二维情况为例，各向同性介质中的一阶速度-应力声波方程为
其中，vx、vz是质点在x和z方向上的震动速度(m/s)；p为位移(m)；ρ为密度
(kg/m3)；v为纵波速度(m/s)；x、z分别为空间坐标(m)；t为时间(s)。
(1)式在交错网格空间中的高阶差分格式[21]为
其中，n为时间离散序号；i、j分别为x和z方向的空间离散序号；Δt为时间步长；
Δx、Δz分别为空间x和z方向的网格大小；L为空间差分的阶数；Dm表示阶数
m=1,2,…，L时的差分系数，Dm的计算公式为[22]
m-1a=1La=m+1
式中，tmax为记录长度(ms)；Su(x,z,t)、Sd(x,z,t)、Sl(x,z,t)、Sr(x,z,t)分别为炮点
波场的上、下、左、右行波；Ru(x,z,t)、Rd(x,z,t)、Rl(x,z,t)、Rr(x,z,t)为分别为接
收点波场的上、下、左、右行波；S(x,z,t)为炮点波场；I(x,z)为对应网格点偏移成
像结果。
吸收边界条件是波场延拓的重要组成部分。目前，业界常用的PML吸收边界在消
除截断边界伪反射和改进成像质量方面发挥了巨大作用。但分裂式PML吸收边界
条件需要对边界区域进行分块处理，从而增加了逆时偏移的内存需求，并使软件开
发工作变得繁琐复杂。因此，为了克服上述缺陷，本文对逆时偏移的无分裂PML
吸收边界条件及其延拓算法进行研究。
3.1 分裂式PML吸收边界条件及其局限
依据PML方程的分裂思路[3]，(1)式中的第一个方程的吸收边界条件为：
p=px+pz
其中，px和pz分别为p分量在x与z方向的分裂算子；σx和σz分别为x和z
方向的衰减因子。当σx和σz为零时，式(6)退化成为式(1)。
逆时偏移时将计算空间划分成PML镶边区域和成像区域两个区域(图1)，具体的
实现方法有两种：①在整个计算空间都采用相同的分裂算法，但不同衰减区域的衰
减因子不同；②PML镶边区域采用式(6)计算，成像区域采用式(1)计算。其中，第
一种方法易于编程实现但使内存需求增加了一倍；第二种方法虽然没有增加成像区
域的内存需求，只增加了PML镶边区域内存需求，但对不同区域需要采用不同的
算法。两套算法耦合使编程变得繁琐，特别是在移植入GPU端时需对各个区域分
配不同的线程，这种分散区块的处理方式不利于提高其计算效率。
3.2 无分裂PML吸收边界
为了克服分裂式PML吸收边界条件在内存需求和GPU并行等方面的不足，本文
对此进行改进，并给出一种无分裂PML吸收边界条件。
式(1)在复数坐标系下做傅氏变换[26,27]，可得到：
其中为频率域中的位移(m)；为频率域中质点在x和z方向上的振动速度(m/s)；ρ
为密度(kg/m3)；v为纵波速度(m/s)。
在频率域复数坐标系下，分别对位移分量和速度分量做衰减，经推导可得到(7)式
在复数坐标系下的无分裂PML边界条件的控制方程[27]，然后将得到的控制方程
转换到时间域，可得到声波方程的无分裂吸收边界条件为
其中， Ovx、Ovz、Opx、Opz为中间辅助量，其控制方程[27]为
式(8)和式(9)即为本文给出的声波方程无分裂PML吸收边界条件。方程组(8)右端
为正常项与衰减项的加和，因此时间域内的求解可以通过以下两个步骤来完成：①
对包含PML区域的整个计算区域进行常规方程求解；②在第一步中减去相关PML
区域的衰减影响。对应图1，无分裂PML吸收边界条件可以将区域3并入区域1、
2，不需要对其进行单独处理。同时，该方法不需要对各个波场分量进行分裂，避
免了各分量分裂算子的额外存储。
3.3 无分裂PML吸收边界条件的实现步骤
以式(8)中第一个方程式为例，逆时偏移时间域的迭代求解过程如下：
1)p分量在x、z方向上的全局更新：
2)p分量在x、z方向上PML区域内减去衰减量：
3)衰减量在x、z方向上PML区域内的更新：
3.4 无分裂PML边界条件的吸收效果
式(6)和式(9)中的σx、σz是相同的衰减因子，理论上只要它们的各参数取值相同，
就具有相同的衰减效果。为验证无分裂PML边界吸收效果，对图2所示的水平层
状模型进行炮点的正向延拓和接收点的逆时延拓。其中，接收点波场记录采用如下
参数正演得到：纵波源激发，震源主频为35 Hz的雷克子波，中间放炮，炮点两
侧共300道接收，道间距为0.005 km，记录长度为1.25 s，PML层厚度为10层。
图3和图4分别为600 ms时刻的无分裂和分裂边界条件下的炮点波场快照；图5
和图6分别为600 ms时刻无分裂和分裂边界条件下的接收点波场快照。图7、图
8、图9、图10分别为600 ms时刻无分裂和分裂边界条件下炮点波场PML区域
放大10倍、100倍的快照。由图3～图10可知：两种边界对炮点波场和接收点
波场的吸收效果是相同的，它们均很好地吸收了入射到计算边界的外行波。
4.1 无分裂PML吸收边界条件在GPU并行逆时偏移中的实现
采用基于CUDA架构的GPU和CPU联合并行[17-19]计算技术实现声波方程的逆
时偏移，其具体流程为：①将地震记录和初始波场值由主机内存传至设备显存；②
将GPU的处理器分为若干个计算区块，每一个区块的线程对应数组中的一个单元；
③利用GPU进行炮点波场的正向延拓，其中边界条件采用本文提出的无分裂PML
吸收边界条件，延拓过程中依据有效边界理论[17,28]记录下边界处部分网格点的
波场值；④依据有效边界理论进行炮点波场重构[28]，同时进行接收点波场逆时延
拓，并在延拓过程中实现各成像点处炮检波场的互相关，在炮点波场正向延拓和接
收点波场逆时延拓过程中仍然采用本文提出的无分裂PML边界条件来吸收截断边
界处的外行波；⑤输出成像结果，并释放主机端和设备端的内存与显存空间；⑥重
复1～5步直至完成所有炮集的偏移成像；⑦从偏移结果中抽出共成像点道集并叠
加。
4.2 两种边界条件的GPU加速比分析
对Marmousi速度模型(图11)的合成单炮数据进行逆时偏移处理，单炮记录按以
下观测系统，通过双程声波方程有限差分正演获得：纵波源激发，震源主频为35
Hz的雷克子波，中间放炮，炮点两侧共2 301道接收，测线两端固定不动，只移
动炮点，道间距为0.005 km，炮间距为0.05 km，采样间隔为0.35 ms，记录长
度为3.5 s，共得200炮合成记录。分别基于分裂式PML吸收边界条件与无分裂
PML吸收边界条件对所有炮集进行逆时偏移处理并作共成像点叠加，可得到该模
型数据的逆时偏移叠加剖面(图12、图13)。对比图12和图13可以得到，两种边
界条件均可以得到准确的成像结果。
表1为采用分裂式PML吸收边界条件和无分裂PML吸收边界条件算法的运算时
间对比，其中程序运行所用的GPU型号为GTX560Ti，拥有448个流处理器(sp)、
2G显存，显存频率为4400MHz，核心频率为950MHz，显存位宽为256 bit。
可见，基于无分裂PML吸收边界条件的逆时延拓可将二维声波方程的GPU加速
比提升至分裂式PML吸收边界条件时的1.2倍。
1)本文推导的声波方程逆时偏移的无分裂PML吸收边界条件具有与常规分裂算法
相同的炮检波场吸收效果。
2)与分裂式PML吸收边界条件相比，本文算法无需在边界处对各个分量进行分裂
处理，避免了常规方法中各分量分裂算子的额外内存开销。
3)Marmousi模型算例表明：基于无分裂PML吸收边界条件的逆时延拓相比于分
裂式PML吸收边界条件，在GPU端具有1.2倍左右的加速比。
4)本文思路不仅适用于各向同性介质中的声波方程，还可用于各向异性和弹性波方
程的逆时偏移。
【相关文献】
[1]McMechan G A. Migration by extrapolation of time dependent boundary values[J].
Geophysical Prospecting, 1983, 31(2):413-420.
[2]Sun R, McMechan G A,Lee C S,et al. Prestack scalar reverse-time depth migration of 3D
elastic seismic data [J]. Geophysics, 2006, 71(5):S199-S207.
[3]何兵寿，张会星，魏修成，等.双程声波方程叠前逆时深度偏移的成像条件[J].石油地球物理勘探，
2010，45(2):237-243.
[4]Chattopadhyay S, McMechan G A. Imaging conditions for prestack reverse-time
migration [J]. Geophysics, 2008, 73(3):S81-S89.
[5]Loewenthal D, Hu L. Two methods for computing the imaging conditions for common
shot pre-stack migration [J].Geophysics, 1991, 56(3): 378-381.
[6]Steve T H. Reverse-time depth migration: Ipad2 and imaging condition[J]. Geophysics,
1987, 52(8): 1 060.
[7]Kaelin B, Guitton A. Imaging conditions for reverse time migration[C]// SEG Technical
Program Expanded Abstracts. USA: Society of Exploration Geophysics, 2006:2 594-2 598.
[8]Schleicher J,Novais A, Costa J. A comparison of imaging conditions for wave equation
shot profile migration [J]. Geophysics, 2008,73 (6):S219-S227.
[9]Loewenthal D, Stoffa P L, Faria E L. Suppresing the unwanted reflections of the full wave
equation[J]. Geophysics.1987,52(7):1 007-1 012.
[10]Gutton A, Kaelin B, Biondi B. Least-square attenuation of reverse-time migration
artifacts[C]// SEG Technical Program Expanded Abstracts. USA: Society of Exploration
Geophysics, 2006: 2 348-2 352 .
[11]杜启振，朱钇同，张明强，等.叠前逆时深度偏移低频噪声压制策略研究[J].地球物理学报，
2013，56(7):2 391-2 401.
[12]Kwangjin Yoon, Kurt J Marfurt. Reverse-time migrstion using the Poynting vector[J].
Exploration Geophysics, 2006, 37(9):102-106.
[13]郭鹏，何兵寿，郭敏.基于坡印廷矢量的声波方程叠前逆时偏移成像条件[J].煤炭学报，2011，
36(8):1 291-1 295.
[14]王京，屈念念，程飞.逆时偏移反射波压制方法及对比分析[J].工程地球物理学报，2013，
10(4):480-486.
[15]William W Symes. Reverse-time migration with optimal checkpointing[J]. Geophysics,
2007, 72(5):SM213-SM221.
[16]Eric Dussaud, William W Symes, Paul Williamson, et al. Computational strategies for
reverse-time migrations[C]// SEG Technical Program Expanded Abstracts. USA: Society of
Exploration Geophysics, 2008:2 267-2 271.
[17]Robert G Clapp. Reverse-time migration with randon boundaries[C]// SEG Technical
Program Expanded Abstracts. USA: Society of Exploration Geophysics, 2009:2 809-2 813.
[18]李博，刘国峰，刘洪.地震叠前时间偏移的一种图形处理器提速实现方法[J].地球物理学报，
2009，52(1):242-252.
[19]刘红伟，李博，刘洪，等.地震叠前逆时偏移高阶有限差分算法及GPU实现[J]. 地球物理学报，
2010，53(7):1 725-1 733.
[20]刘红伟，刘洪，邹振，等.地震叠前逆时偏移中去噪与存储[J].地球物理学报，2010，53(9):2
171-2 180.
[21]董良国，马在田，曹景忠，等.一阶弹性波方程交错网格高阶差分算法解法[J].地球物理学报，
2000，43(3):411-419.
[22]牟永光，裴正林.三维复杂介质地震数值模拟[M].北京：石油工业出版社，2005:34-34.
[23]Zeng Y, He J, Liu Q H. The application of the perfectly matched layer in numerical
modeling of wave propagation in poroelastic media[J].Geophysics,2001,66(4):1 258-1 266.
[24]Chen T, He B S. A normalized wave-field separation cross-correlation imaging
condition for reverse-time migration based on Poynting vector[J]. Applied Geophysics,
2014,11(2):158-166.
[25]Tang C, Wang D L. Reverse-time mifration with source wave-field reconstruction and
wave-field decomposition[J]. Global Geology, 2012,31(4):803-812.
[26]秦臻，任培罡，姚姚，等.弹性波正演模拟中PML吸收边界条件的改进[J].地球科学(中国地质
大学学报)，2009，34(4):659-663.
[27]潘威，成景旺，闫立爽.二维弹性波方程频率空间域有限差分正演研究[J].工程地球物理学报，
2013，10(6):858-862.
[28]王保利，高静怀，陈文超，等.地震叠前逆时偏移的有效边界存储策略[J].地球物理学报，2012，
55(7):2 413-2 421.
2024年5月16日发(作者：褚振博)
声波方程逆时偏移中的无分裂PML吸收边界条件
王鹏飞;何兵寿
【摘 要】逆时偏移是当前地震资料处理的前沿技术,吸收边界条件是逆时偏移技术
的重要组成部分.目前,业界常用的基于波场分裂思路的完全匹配层PML(Perfectly
Matched Layer)吸收边界条件在改善边界吸收效果和改进偏移成像质量方面发挥
了巨大作用,但这种方法需要在边界处为分裂后的各个分量开辟额外的内存空间,且
需要分边角处理,增加了逆时偏移技术的内存负担和计算开销.为了减少内存负担和
提高计算效率,首先从双程声波方程出发,推导了声波方程的无分裂PML吸收边界
条件,然后给出该边界条件下波动方程逆时延拓的数值实现过程.理论分析和模型实
验结果得到:无分裂PML边界条件具有与分裂算法相同的边界吸收效果,并且基于
无分裂PML的逆时偏移算法效率更高,更便于程序代码的编写和GPU
(GraphicProcessing Unit)的并行实现.
【期刊名称】《工程地球物理学报》
【年(卷),期】2015(012)005
【总页数】8页(P583-590)
【关键词】逆时偏移;吸收边界条件;无分裂完全匹配层;计算效率
【作 者】王鹏飞;何兵寿
【作者单位】中国海洋大学海底科学与探测技术教育部重点实验室,山东青岛
266100;中国海洋大学海底科学与探测技术教育部重点实验室,山东青岛266100
【正文语种】中 文
【中图分类】P631.5
逆时偏移是当前地震资料成像处理领域的热点技术。目前，国内外在逆时偏移领域
的研究工作主要集中在以下几个方面：①波动方程逆时延拓算法[1,2]研究，包括
差分格式推导、边界伪反射压制和数值频散压制等；②成像准则和成像方法[3-8]
研究；③逆时偏移噪声压制方法[9-14]研究，主要研究波场延拓过程中由层间反射
导致的低频噪音的压制方法；④逆时偏移存储策略[15-17]研究，主要用于降低逆
时偏移的临时文件存储量和硬盘访问量，拓展逆时偏移处理的并行性并提高处理效
率；⑤偏移处理的并行算法[18-20]研究，包括CPU(Central Processing Unit)并
行、GPU并行以及CPU+GPU协同并行。
本文内容属于逆时延拓算法和逆时偏移并行处理范畴，对声波方程逆时偏移中常用
的PML吸收边界条件进行了改进，使之由分裂形式变为无分裂形式，避免了计算
区域的分割，降低了内存开销，提升了计算效率，且简化了代码编写工作。本文在
GPU计算平台CUDA(Compute Unified Device Architecture)上实现了无分裂
PML条件的逆时延拓。模型实验结果表明，无分裂PML边界条件具有与分裂算法
相同的边界吸收效果，且降低了内存消耗，减少了运算次数，更便于程序代码的编
写和GPU的并行实现。
逆时偏移可以通过以下三步实现：①炮点波场正向延拓；②接收点波场逆时延拓；
③炮检波场互相关成像。其中炮点波场的正向延拓方法与对应地震波方程的正演方
法相同，本文不赘述，只简要介绍接收点波场逆时延拓方法与波场成像方法。
2.1 接收点波场的逆时延拓
以二维情况为例，各向同性介质中的一阶速度-应力声波方程为
其中，vx、vz是质点在x和z方向上的震动速度(m/s)；p为位移(m)；ρ为密度
(kg/m3)；v为纵波速度(m/s)；x、z分别为空间坐标(m)；t为时间(s)。
(1)式在交错网格空间中的高阶差分格式[21]为
其中，n为时间离散序号；i、j分别为x和z方向的空间离散序号；Δt为时间步长；
Δx、Δz分别为空间x和z方向的网格大小；L为空间差分的阶数；Dm表示阶数
m=1,2,…，L时的差分系数，Dm的计算公式为[22]
m-1a=1La=m+1
式中，tmax为记录长度(ms)；Su(x,z,t)、Sd(x,z,t)、Sl(x,z,t)、Sr(x,z,t)分别为炮点
波场的上、下、左、右行波；Ru(x,z,t)、Rd(x,z,t)、Rl(x,z,t)、Rr(x,z,t)为分别为接
收点波场的上、下、左、右行波；S(x,z,t)为炮点波场；I(x,z)为对应网格点偏移成
像结果。
吸收边界条件是波场延拓的重要组成部分。目前，业界常用的PML吸收边界在消
除截断边界伪反射和改进成像质量方面发挥了巨大作用。但分裂式PML吸收边界
条件需要对边界区域进行分块处理，从而增加了逆时偏移的内存需求，并使软件开
发工作变得繁琐复杂。因此，为了克服上述缺陷，本文对逆时偏移的无分裂PML
吸收边界条件及其延拓算法进行研究。
3.1 分裂式PML吸收边界条件及其局限
依据PML方程的分裂思路[3]，(1)式中的第一个方程的吸收边界条件为：
p=px+pz
其中，px和pz分别为p分量在x与z方向的分裂算子；σx和σz分别为x和z
方向的衰减因子。当σx和σz为零时，式(6)退化成为式(1)。
逆时偏移时将计算空间划分成PML镶边区域和成像区域两个区域(图1)，具体的
实现方法有两种：①在整个计算空间都采用相同的分裂算法，但不同衰减区域的衰
减因子不同；②PML镶边区域采用式(6)计算，成像区域采用式(1)计算。其中，第
一种方法易于编程实现但使内存需求增加了一倍；第二种方法虽然没有增加成像区
域的内存需求，只增加了PML镶边区域内存需求，但对不同区域需要采用不同的
算法。两套算法耦合使编程变得繁琐，特别是在移植入GPU端时需对各个区域分
配不同的线程，这种分散区块的处理方式不利于提高其计算效率。
3.2 无分裂PML吸收边界
为了克服分裂式PML吸收边界条件在内存需求和GPU并行等方面的不足，本文
对此进行改进，并给出一种无分裂PML吸收边界条件。
式(1)在复数坐标系下做傅氏变换[26,27]，可得到：
其中为频率域中的位移(m)；为频率域中质点在x和z方向上的振动速度(m/s)；ρ
为密度(kg/m3)；v为纵波速度(m/s)。
在频率域复数坐标系下，分别对位移分量和速度分量做衰减，经推导可得到(7)式
在复数坐标系下的无分裂PML边界条件的控制方程[27]，然后将得到的控制方程
转换到时间域，可得到声波方程的无分裂吸收边界条件为
其中， Ovx、Ovz、Opx、Opz为中间辅助量，其控制方程[27]为
式(8)和式(9)即为本文给出的声波方程无分裂PML吸收边界条件。方程组(8)右端
为正常项与衰减项的加和，因此时间域内的求解可以通过以下两个步骤来完成：①
对包含PML区域的整个计算区域进行常规方程求解；②在第一步中减去相关PML
区域的衰减影响。对应图1，无分裂PML吸收边界条件可以将区域3并入区域1、
2，不需要对其进行单独处理。同时，该方法不需要对各个波场分量进行分裂，避
免了各分量分裂算子的额外存储。
3.3 无分裂PML吸收边界条件的实现步骤
以式(8)中第一个方程式为例，逆时偏移时间域的迭代求解过程如下：
1)p分量在x、z方向上的全局更新：
2)p分量在x、z方向上PML区域内减去衰减量：
3)衰减量在x、z方向上PML区域内的更新：
3.4 无分裂PML边界条件的吸收效果
式(6)和式(9)中的σx、σz是相同的衰减因子，理论上只要它们的各参数取值相同，
就具有相同的衰减效果。为验证无分裂PML边界吸收效果，对图2所示的水平层
状模型进行炮点的正向延拓和接收点的逆时延拓。其中，接收点波场记录采用如下
参数正演得到：纵波源激发，震源主频为35 Hz的雷克子波，中间放炮，炮点两
侧共300道接收，道间距为0.005 km，记录长度为1.25 s，PML层厚度为10层。
图3和图4分别为600 ms时刻的无分裂和分裂边界条件下的炮点波场快照；图5
和图6分别为600 ms时刻无分裂和分裂边界条件下的接收点波场快照。图7、图
8、图9、图10分别为600 ms时刻无分裂和分裂边界条件下炮点波场PML区域
放大10倍、100倍的快照。由图3～图10可知：两种边界对炮点波场和接收点
波场的吸收效果是相同的，它们均很好地吸收了入射到计算边界的外行波。
4.1 无分裂PML吸收边界条件在GPU并行逆时偏移中的实现
采用基于CUDA架构的GPU和CPU联合并行[17-19]计算技术实现声波方程的逆
时偏移，其具体流程为：①将地震记录和初始波场值由主机内存传至设备显存；②
将GPU的处理器分为若干个计算区块，每一个区块的线程对应数组中的一个单元；
③利用GPU进行炮点波场的正向延拓，其中边界条件采用本文提出的无分裂PML
吸收边界条件，延拓过程中依据有效边界理论[17,28]记录下边界处部分网格点的
波场值；④依据有效边界理论进行炮点波场重构[28]，同时进行接收点波场逆时延
拓，并在延拓过程中实现各成像点处炮检波场的互相关，在炮点波场正向延拓和接
收点波场逆时延拓过程中仍然采用本文提出的无分裂PML边界条件来吸收截断边
界处的外行波；⑤输出成像结果，并释放主机端和设备端的内存与显存空间；⑥重
复1～5步直至完成所有炮集的偏移成像；⑦从偏移结果中抽出共成像点道集并叠
加。
4.2 两种边界条件的GPU加速比分析
对Marmousi速度模型(图11)的合成单炮数据进行逆时偏移处理，单炮记录按以
下观测系统，通过双程声波方程有限差分正演获得：纵波源激发，震源主频为35
Hz的雷克子波，中间放炮，炮点两侧共2 301道接收，测线两端固定不动，只移
动炮点，道间距为0.005 km，炮间距为0.05 km，采样间隔为0.35 ms，记录长
度为3.5 s，共得200炮合成记录。分别基于分裂式PML吸收边界条件与无分裂
PML吸收边界条件对所有炮集进行逆时偏移处理并作共成像点叠加，可得到该模
型数据的逆时偏移叠加剖面(图12、图13)。对比图12和图13可以得到，两种边
界条件均可以得到准确的成像结果。
表1为采用分裂式PML吸收边界条件和无分裂PML吸收边界条件算法的运算时
间对比，其中程序运行所用的GPU型号为GTX560Ti，拥有448个流处理器(sp)、
2G显存，显存频率为4400MHz，核心频率为950MHz，显存位宽为256 bit。
可见，基于无分裂PML吸收边界条件的逆时延拓可将二维声波方程的GPU加速
比提升至分裂式PML吸收边界条件时的1.2倍。
1)本文推导的声波方程逆时偏移的无分裂PML吸收边界条件具有与常规分裂算法
相同的炮检波场吸收效果。
2)与分裂式PML吸收边界条件相比，本文算法无需在边界处对各个分量进行分裂
处理，避免了常规方法中各分量分裂算子的额外内存开销。
3)Marmousi模型算例表明：基于无分裂PML吸收边界条件的逆时延拓相比于分
裂式PML吸收边界条件，在GPU端具有1.2倍左右的加速比。
4)本文思路不仅适用于各向同性介质中的声波方程，还可用于各向异性和弹性波方
程的逆时偏移。
【相关文献】
[1]McMechan G A. Migration by extrapolation of time dependent boundary values[J].
Geophysical Prospecting, 1983, 31(2):413-420.
[2]Sun R, McMechan G A,Lee C S,et al. Prestack scalar reverse-time depth migration of 3D
elastic seismic data [J]. Geophysics, 2006, 71(5):S199-S207.
[3]何兵寿，张会星，魏修成，等.双程声波方程叠前逆时深度偏移的成像条件[J].石油地球物理勘探，
2010，45(2):237-243.
[4]Chattopadhyay S, McMechan G A. Imaging conditions for prestack reverse-time
migration [J]. Geophysics, 2008, 73(3):S81-S89.
[5]Loewenthal D, Hu L. Two methods for computing the imaging conditions for common
shot pre-stack migration [J].Geophysics, 1991, 56(3): 378-381.
[6]Steve T H. Reverse-time depth migration: Ipad2 and imaging condition[J]. Geophysics,
1987, 52(8): 1 060.
[7]Kaelin B, Guitton A. Imaging conditions for reverse time migration[C]// SEG Technical
Program Expanded Abstracts. USA: Society of Exploration Geophysics, 2006:2 594-2 598.
[8]Schleicher J,Novais A, Costa J. A comparison of imaging conditions for wave equation
shot profile migration [J]. Geophysics, 2008,73 (6):S219-S227.
[9]Loewenthal D, Stoffa P L, Faria E L. Suppresing the unwanted reflections of the full wave
equation[J]. Geophysics.1987,52(7):1 007-1 012.
[10]Gutton A, Kaelin B, Biondi B. Least-square attenuation of reverse-time migration
artifacts[C]// SEG Technical Program Expanded Abstracts. USA: Society of Exploration
Geophysics, 2006: 2 348-2 352 .
[11]杜启振，朱钇同，张明强，等.叠前逆时深度偏移低频噪声压制策略研究[J].地球物理学报，
2013，56(7):2 391-2 401.
[12]Kwangjin Yoon, Kurt J Marfurt. Reverse-time migrstion using the Poynting vector[J].
Exploration Geophysics, 2006, 37(9):102-106.
[13]郭鹏，何兵寿，郭敏.基于坡印廷矢量的声波方程叠前逆时偏移成像条件[J].煤炭学报，2011，
36(8):1 291-1 295.
[14]王京，屈念念，程飞.逆时偏移反射波压制方法及对比分析[J].工程地球物理学报，2013，
10(4):480-486.
[15]William W Symes. Reverse-time migration with optimal checkpointing[J]. Geophysics,
2007, 72(5):SM213-SM221.
[16]Eric Dussaud, William W Symes, Paul Williamson, et al. Computational strategies for
reverse-time migrations[C]// SEG Technical Program Expanded Abstracts. USA: Society of
Exploration Geophysics, 2008:2 267-2 271.
[17]Robert G Clapp. Reverse-time migration with randon boundaries[C]// SEG Technical
Program Expanded Abstracts. USA: Society of Exploration Geophysics, 2009:2 809-2 813.
[18]李博，刘国峰，刘洪.地震叠前时间偏移的一种图形处理器提速实现方法[J].地球物理学报，
2009，52(1):242-252.
[19]刘红伟，李博，刘洪，等.地震叠前逆时偏移高阶有限差分算法及GPU实现[J]. 地球物理学报，
2010，53(7):1 725-1 733.
[20]刘红伟，刘洪，邹振，等.地震叠前逆时偏移中去噪与存储[J].地球物理学报，2010，53(9):2
171-2 180.
[21]董良国，马在田，曹景忠，等.一阶弹性波方程交错网格高阶差分算法解法[J].地球物理学报，
2000，43(3):411-419.
[22]牟永光，裴正林.三维复杂介质地震数值模拟[M].北京：石油工业出版社，2005:34-34.
[23]Zeng Y, He J, Liu Q H. The application of the perfectly matched layer in numerical
modeling of wave propagation in poroelastic media[J].Geophysics,2001,66(4):1 258-1 266.
[24]Chen T, He B S. A normalized wave-field separation cross-correlation imaging
condition for reverse-time migration based on Poynting vector[J]. Applied Geophysics,
2014,11(2):158-166.
[25]Tang C, Wang D L. Reverse-time mifration with source wave-field reconstruction and
wave-field decomposition[J]. Global Geology, 2012,31(4):803-812.
[26]秦臻，任培罡，姚姚，等.弹性波正演模拟中PML吸收边界条件的改进[J].地球科学(中国地质
大学学报)，2009，34(4):659-663.
[27]潘威，成景旺，闫立爽.二维弹性波方程频率空间域有限差分正演研究[J].工程地球物理学报，
2013，10(6):858-862.
[28]王保利，高静怀，陈文超，等.地震叠前逆时偏移的有效边界存储策略[J].地球物理学报，2012，
55(7):2 413-2 421.