弹性波数值模拟中的高斯型混合吸收边界条件及其GPU并行

王绍文 宋 鹏*②③ 谭 军②③ 解 闯 毛士博 王倩倩

(①中国海洋大学海洋地球科学学院，山东青岛 266100；②青岛国家海洋科学与技术实验室，山东青岛 266100；③中国海洋大学海底科学与探测技术教育部重点实验室，山东青岛 266100)

0 引言

地震波正演模拟技术是研究地震波传播规律的有效手段，并在地震勘探的采集、处理以及反演等各个环节中都有着广泛的应用。地震波正演模拟必须引入人工边界界定计算区域，而人工边界的处理不当，往往会产生虚假反射污染中心波场，从而降低波场模拟的精度，因而人工边界的处理一直是波动方程数值计算的重要研究内容。

当前的人工边界处理方法主要有三类。第一类是衰减边界方法[1-3]。该类方法是在靠近边界处引入一个衰减区，在衰减区内地震波传播遵循阻尼波动方程，能量按指数衰减。由于该类方法存在阻尼因子难以确定、计算量大和内存消耗较大等问题，因此应用较少。第二类是完全匹配层(PML)方法。该类方法由Berenger[4]在电磁波数值模拟中首先提出，通过在中心波场计算区域外加上一系列的吸收层、层内引入衰减因子以达到消除边界反射的目的。王守东[5]、Hastings等[6]、Collino等[7]将PML法成功应用于声波和弹性波正演模拟。随后， Komatitsch等[8-9]、熊章强等[10]实现了非分裂PML法；罗玉钦等[11]、杨茜娜等[12]发展了复频移及多轴复频移PML法；陈可洋[13]和Wang等[14]提出了基于余弦型和基于高斯型衰减函数的PML方法。PML法现已在地震波数值模拟中获得广泛应用，然而为达到更好的吸收效果，PML法往往需要设计几十甚至上百个匹配层，意味着巨大的计算量和内存消耗，严重限制了该方法在大规模波场模拟中的应用。第三类为基于单程波近似的吸收边界条件。该类方法由Enguist等[15-16]首先提出；Higdon[17-19]提出了一种与其等价的渐进吸收边界条件，可以设计吸收任意角度的入射波，显著提升了吸收精度和适应性。然而由于单程波的高阶近似会导致边界条件方程中出现高阶的混合时间、空间偏导数项，给数值求解带来困难，因此常规方法多采用二阶边界条件，但对大角度入射波的吸收效果较差，其整体吸收效果也难以满足高精度数值模拟的要求。为进一步提升基于单程波近似的吸收边界条件的效果，宋鹏等[20-21]发展了优化系数的四阶吸收边界条件算法。Baffet等[22]、Hagstrom等[23]以及Potter等[24]也通过引入辅助变量实现了高阶吸收边界条件，但这些算法同样存在计算过程复杂且计算效率低等问题。

为提升常规基于单程波近似的吸收边界条件的效果，在不提高阶数的前提下，Liu等[25]提出了一种混合吸收边界条件。该方法在中心波场与边界之间增加一个过渡区域，通过线性加权使过渡区域内双程波波场与单程波波场平滑过渡，最终实现了人工边界的高效吸收。随后，Liu等[26-27]、任志明等[28]以及Liu等[29]实现了三维声波及弹性波正演模拟中的混合吸收边界，Ren等[30]将混合吸收边界推广到频率域数值模拟中，取得了较好的吸收效果。刘学义等[31]将混合吸收边界用于各向异性介质，Wang等[32]、薛浩等[33]分别将混合吸收边界应用于逆时偏移和最小二乘逆时偏移。

当前混合吸收边界已在地震波场正反演延拓中得到广泛应用，然而，该方法的加权系数仍以线性和指数型[34]为主，未能兼顾内边界和外边界的综合吸收效果，在吸收效率上并没有达到最优；此外，当前大规模的波场模拟往往借助于GPU并行提高计算效率，但混合吸收边界差分格式的GPU加速效率远远低于中心波场，严重影响了整个波场模拟的计算效率。本文提出一种基于高斯型加权系数的混合吸收边界条件方法，其更能兼顾内、外吸收边界的吸收效果，具有更高效的整体边界吸收效率；同时本文还推导了一种更适用于GPU加速的一阶Higdon混合吸收边界条件差分格式，可显著提升边界计算的GPU并行效率。

1 基于高斯型加权系数的混合吸收边界条件

1.1 弹性波数值模拟中混合吸收边界原理

二维空间中的一阶速度—应力弹性波方程[35]可以表示为

(1)

式中：vx、vz分别为水平方向和垂直方向波场速度分量；σxx、σzz、σxz为三个应力分量；ρ为密度，λ和μ为拉梅系数，与纵波速度VP和横波速度VS关系为

(2)

(3)

式(1)的交错网格有限差分格式为[36-37]

(4)

(5)

(6)

(7)

(8)

式中：U和W分别表示速度分量vx、vz的离散形式；P、Q和R分别表示应力分量σxx、σzz、σxz的离散形式；Δt、Δx、Δz分别为离散的时间步长、空间水平和垂直网格间隔；i、j、k为空间和时间离散序号；2N为差分阶次；am(1≤m≤N)为差分系数，可通过求解下式方程获取[38]

(9)

应用混合吸收边界方法进行数值模拟时，整个计算区域被划分为中心波场区域Ⅰ、边界过渡区域 Ⅱ 以及边界区Ⅲ三部分(图1)，其波场计算流程如下：

图1 二维混合吸收边界条件示意图

(3)对于M层区域Ⅱ，对双程波场值和单程波场值进行加权求和得到混合波场

ul=(1-wl)ut+wlu°l=1，2，…，M

(10)

式中wl为加权系数。

而对于区域Ⅱ和区域Ⅲ的单程波场计算，由于二阶Higdon吸收边界条件稳定性较差，通常采用改进的一阶Higdon吸收边界条件[28]，其左边界方程为

(11)

(12)

则可将一阶Higdon左边界算子表示为[19]

(13)

其差分格式为

(14)

式中与空间位置和网格剖分相关的常量为[19]

(15)

(16)

(17)

其中b为常数，通常取值为0.5。

混合边界的角点区域需做特殊处理，以左上角速度水平分量为例，其角点差分格式为

(18)

1.2 基于高斯型加权系数的混合吸收边界

混合边界通过引入过渡区实现波场混叠压制边界反射，过渡区与中心波场相邻的边界可称为内边界，过渡区的最外层边界可称为外边界，混合边界的残余反射通常来自于内边界和外边界反射。内、外边界的残余反射能量可通过混合加权系数的选取进行调整，常用的线性加权系数和指数型加权系数[34]分别为

(19)

(20)

线性加权系数(图2黑线所示)在内、外边界处梯度变化相同，而由于内边界更靠近内部波场，因此应用线性加权系数的混合边界其通常会产生较强的内边界反射；指数型加权系数(图2绿线所示)在内边界处变化更加缓慢，使内边界与中心波场区域得到更好的耦合，因此其内边界残余反射得到有效压制，但另一方面，由于其在内边界处的变化过于缓慢导致外边界处的变化过于剧烈，从而易产生较强的外边界反射。因此一个更为理想的混合边界条件必须全面考虑地震波在内、外边界处的综合边界反射压制效果。

鉴于以上分析，本文提出了一种高斯型混合加权系数

(21)

式中α为高斯型权系数的阶数，控制内、外边界处的变化率。大量的实验证明，当α=2.5时，高斯型混合边界既能保持稳定性，又可达到良好的边界吸收效果。

高斯型加权系数(图2红线)在内边界处梯度变化虽然比指数型加权系数稍显剧烈，但其较线性加权系数已更为平缓；而外边界处，高斯型加权系数又比指数型加权系数更平缓，较线性加权系数稍剧烈，因此理论上其应有更优的内、外边界残余反射综合压制效果。

当然，采用混合吸收边界时，整个边界的吸收效果还与过渡区层数有关，Liu等[27]证明了当过渡区厚度较小时，混合边界易不稳定，并且考虑到为追求更高精度的波场模拟效果，当前数值模拟时通常大多采用20层以上的混合吸收边界。因此本文将以20层和30层为例讨论高斯型混合加权系数的吸收效果。

图2 线性、指数型、高斯型加权系数变化曲线

1.3 常速模型数值模拟实验

为验证本文提出的高斯型混合吸收边界条件的吸收效果，本文设计3200m×3200m、VP=4000m/s、VS=2000m/s、ρ= 2100kg/m3的均匀介质模型。纵、横向空间网格间距均为8m，时间步长为1ms，差分格式精度为空间10阶。震源位于模型中心，采用主频为25Hz的Ricker子波加载于速度垂直分量上，记录长度为2s。

应用20层混合吸收边界得到1s时刻的vx、vz分量波前快照如图3所示，波前快照的水平和垂直切片(位置见图3红色虚线)如图4所示。应用30层混合吸收边界得到1s时刻的vx、vz分量波前快照如图5所示，波前快照的水平和垂直切片(位置见图5红色虚线)如图6所示。

由图3～图6可见，无论是20层还是30层吸收边界，对于vx、vz分量，线性加权系数其波前快照中仍存在较强的内边界残余反射(图3左、图5左箭头所示)，指数型加权系数其波前快照中则存在着明显的外边界残余反射(图3中、图5中箭头所示)，而高斯型加权系数的波前快照中无明显边界残余反射波。

图3 20层混合吸收边界不同加权系数的1s时刻vx分量(a)、vz分量(b)波场模拟快照对比 左：线性；中：指数型；右：高斯型

图4 20层混合吸收边界不同加权系数的1s时刻vx分量(a)、vz分量(b)波场模拟切片对比

图5 30层混合吸收边界不同加权系数的1s时刻vx分量(a)、vz分量(b)波场模拟快照对比 左：线性；中：指数型；右：高斯型

图6 30层混合吸收边界不同加权系数的1s时刻vx分量(a)、vz分量(b)波场模拟切片对比

2 适用于GPU加速的一阶Higdon混合吸收边界条件差分格式

当前大规模的波场模拟往往借助于GPU并行提高计算效率。但应用于大规模波场模拟的混合边界差分格式的GPU加速比远远低于中心波场计算的GPU加速比，严重影响了整个波场模拟的计算效率。本文推导了一种新的混合边界差分计算格式，可显著提高边界计算的GPU加速比。

2.1 适用于GPU并行的混合吸收边界差分格式

由式(14)可知，基于常规混合吸收边界差分格式进行GPU并行时，由于每层边界上k+1时刻的波场值计算都要用到该层内侧点k+1时刻的波场值，所以应用GPU计算边界时无法同时启动大量的线程实现边界各点值的同时计算，无法实现高效的GPU并行。为了解决常规边界差分格式中需要由内向外依次计算、无法实现GPU高效并行的问题，本文推导了一种新的边界条件差分格式，具体推导过程如下(以左边界和左上角点为例)。

为避免边界区域计算时，差分格式的右侧不出现k+1时刻项，将时间偏导数和水平方向空间偏导数分别表示为

(22)

(23)

则可将一阶Higdon边界算子表示为

(24)

(25)

对于左上角点，可通过上边界和左边界的加权平均得到

(26)

对比式(14)与式(25)可知，当计算k+1时刻波场值时，本文推导的差分格式只需用到已计算出的k时刻波场值，因此同一时间切片内所有网格点可同时计算，能实现GPU的高效并行加速。

在常规混合吸收边界算法中，由于第l层边界计算与第l+1层边界计算存在计算先后顺序关系，过渡区Ⅱ中的混合波场计算需要由最内侧开始逐层向外进行，因此同一时间层内，前者需多次启动边界计算核函数(与边界层数相当)，且每次核函数只能计算本吸收层内的网格点数，严重降低了GPU的加速效率。基于本文差分格式的混合吸收边界算法，只需一次启动和计算边界计算核函数，可实现所有吸收层内的所有网格点的同时计算，更适合于大规模波场模拟的GPU高效并行计算。

2.2 常速模型数值实验

应用上述均匀模型和模拟参数，验证本文推导的混合吸收边界差分格式的正确性。20、30层常规差分格式与本文差分格式vx、vz分量波前快照如图7、图8所示。

由图7、图8可以看出，在边界层数相同的条件下，两种差分格式的混合吸收边界均能有效压制边界反射。

图7 不同边界差分格式的20层混合吸收边界数值模拟的0.6(左)、1.0(中)和1.4s(右)时刻波前快照对比(a)常规差分格式vx分量； (b)常规差分格式vz分量； (c)本文差分格式vx分量； (d)本文差分格式vz分量

为评估两种差分格式的计算效率，本文基于RTX-2080Ti-11GB的GPU应用上述均匀模型进行测试，结果如表1所示。

由表1可见，若采用20层边界，本文边界差分格式计算效率为常规差分格式3.25倍；若采用30层边界，本文边界差分格式计算效率为常规差分格式4.34倍，因此本文差分格式更适于大规模波场模拟的GPU高效并行计算。

表1 两种边界差分格式正演模拟计算效率对比

图8 不同边界差分格式30层混合吸收边界数值模拟的0.6(左)、1.0(中)和1.4s(右)时刻波前快照对比(a)常规差分格式vx分量； (b)常规差分格式vz分量； (c)本文差分格式vx分量； (d)本文差分格式vz分量

3 复杂模型实验

本文截取了Marmousi模型的一部分(图9)进行正演模拟测试。该模型尺寸为3000m×3000m，纵横波速比为1.73，常密度为2100kg/m3。数值模拟中，横向和纵向网格间隔分别为5m和4m，时间步长为0.4ms。震源采用加载于vz分量上、主频为25Hz的Ricker子波，位于(1500m，4m)处。

未加吸收边界、常规差分格式的高斯型混合吸收边界(30层)和本文差分格式的高斯型混合吸收边界(30层)模拟的1.4s时vx、vz分量的波前快照如图10所示。由图可见，未加吸收边界的波场中有明显的强边界反射，而无论采用常规边界差分格式还是本文差分格式的30层混合吸收边界，均可有效吸收边界反射，表明基于高斯型混合加权吸收边界能够适应于复杂构造模型的边界处理。

基于单卡GPU(型号：RTX-2080Ti-11GB)，常规差分格式的高斯型混合吸收边界的部分Mar-mousi模型计算时间为16.6335s，本文差分格式计算时间为5.3510s，可见本文差分格式的计算效率是常规格式的3倍以上，说明本文的边界差分格式更适合大规模数值模拟的GPU高效并行。

图9 部分Marmousi纵波速度模型

图10 Marmousi模型不同边界模拟的1.5s时刻的vx分量(左)、vz分量(右)波场快照(a)不加边界； (b)常规差分格式的高斯型混合吸收边界； (c)本文差分格式的高斯型混合吸收边界

4 结论

本文在系统分析线性混合边界和指数型混合边界残余反射形成机制的基础上，提出了一种能兼顾内外边界吸收效果的高斯型混合吸收边界条件。在此基础上，推导了更适用于GPU加速的一阶Higdon混合吸收边界条件差分格式及相应的角点计算差分格式。模型实验结果表明，高斯型混合吸收边界能更有效地压制内边界和外边界反射，比线性和指数型混合吸收边界具有更高的吸收效率，适用于复杂构造模型的边界处理；本文推导的边界及其角点的边界差分格式，GPU并行的计算效率更高，更适用于大规模弹性波场数值模拟。

将本文提出的混合吸收边界算法推广至三维弹性波数值模拟及逆时偏移和全波形反演，是今后的研究方向。

十分感谢中国海洋大学地球探测软件技术研发团队提供了资金和软、硬件支持，以及团队成员提供了建设性意见。