辅助微分方程完全匹配层在声波方程数值模拟中的应用

赵建国，史瑞其，陈竞一，赵维俊，王宏斌，潘建国

1.中国石油大学（北京）地球物理与信息工程学院/油气资源与探测国家重点实验室/CNPC物探重点实验室，北京 102249

2.中国海洋石油研究总院，北京 100027

3.塔尔萨大学地球科学学院，美国 塔尔萨 74104

4.沈阳地质矿产研究所，沈阳 110034

5.中石油勘探开发研究院西北分院，兰州 730000

0 引言

在地震波数值模拟中，为了在有限计算区域模拟无限大区域中地震波的传播，需要对计算区域截断边界给出吸收边界条件以消除边界处产生的虚假数值反射。目前，广泛使用的主要有2种构造吸收边界条件的方式：一种是通过波动方程因式分解获得单行波方程的旁轴近似边界条件，如Clayton－Engquist吸收边界条件［1］、廖氏吸收边界条件［2］、Higdon吸收边界条件［3］等；另一种是在边界上引入吸收介质，使得地震波在无反射地进入吸收介质后被衰减掉，如完全匹配层（pefectly matched layer，PML）［4］。PML作为一种高效的吸收边界，在地震波正演中被广泛使用。PML最初采用对波场变量进行非物理场分裂的方式来实现，使得计算内域与PML吸收边界区域的数值计算格式完全不同，实现方式较为繁琐。后来，复拉伸坐标［5］概念的引入使得PML可以通过非分裂场方式实现。

在使用非分裂场方式实现PML的算法中，需要引入辅助变量并在时间域迭代中对其进行更新。通过辅助变量更新方式的不同，可以划分出一系列不同的PML算法，如通过递归积分［6］和递归卷积［7］算法。以上算法使得辅助变量的时间计算格式局限于时间二阶精度。但是，在波动方程数值模拟中，为了获得更高的模拟精度，常常需要使用高阶时间递推格式；因此，上述非分裂PML在这种情况下的应用受到限制。针对该问题，作者介绍了一种通过辅助微分方程求解辅助变量的 ADE－PML（auxiliary differential equation PML）［8］算法，将其用于高阶时间精度（四阶Runge－Kutta）的声波方程数值模拟中。通过数值模拟，考察ADE－PML在高阶精度时间递推格式声波方程模拟中的应用效果并进行理论分析。

1 方法原理

在二维直角坐标系下，声波方程的一阶压力－速度形式在时间域如下：

其中：p为声压；体积模量K＝ρc2，ρ为质量密度，c为声波在介质中的传播速度；vx和vz分别为沿x和z方向的速度场分量。注意，这里vx和vz是指介质中质点振动速度，c为波在介质中传播的群速度。

本节将介绍一种使用辅助微分方程在时间域更新辅助变量的非分裂完全匹配层算法，这种方法可以适应高阶精度的时间递推格式。本文以四阶Runge－Kutta时间递推格式为例说明ADE－PML在高阶时间递推格式中的应用。为构造PML，需要对空间偏导数做复拉伸坐标变换。为简单起见，首先以x方向为例，只讨论x正半坐标轴的形式，并假设PML区域从x＝0开始。然后将x坐标替换为复拉伸坐标：

其中：sx是复拉伸函数，决定着PML的吸收效果。则关于x方向的偏导可以改写为

z方向同理。因此，只需将关于x，z方向的偏导∂x，∂z替换为复拉伸坐标下的∂˜x和∂˜z后即可。由于PML需在频率域导出，以式（1－1）的频率域方程为例：

其中：i为虚数单位；ω为角频率。当确定复拉伸函数sx，sz的具体形式后，即可将式（4）变换回时间域。传统的复拉伸函数为

其中：dx≥0是使地震波振幅在PML内呈指数级下降的衰减系数。由于式（5）只对行波有效，对平行于PML和计算内域交界传播的隐失波和极低频波无效［9－10］，因此，只能增加 PML厚度或将震源放置于远离PML的位置，用以增加内存和计算时间。为解决这个问题，Kuzuoglu等［11］将复拉伸函数的概念扩展，提出了复频移拉伸函数：

其中：κx和αx为辅助衰减系数；κx≥1，作用是将倾斜入射波矫正为垂直入射波；αx≥0使复平面的极点从实轴移动到虚负半平面上，作用是消除切入射波引起的低频波和隐失波。为避免数值反射，各衰减系数需要渐进变化。当取κx＝1和αx＝0时，式（6）即退化为传统复拉伸函数（5）。笔者主要采用指数函数［12－14］来定义 PML 的衰减系数分布：

其中：L是PML厚度；pd，pα和pκ为常数，通常取pd＝2，pα＝1，pκ＝2；dmax，αmax和κmax的最优值依赖于计算模型参数或所选数值模拟方法。根据文献［12］，衰减系数dmax通过关于垂直入射平面波的目标反射系数R计算：

其中：cmax为计算模型的最大声波速度。根据文献［13］和［14］，辅助衰减系数αmax和κmax由以下关系确定：

其中：cmin为计算模型的最小声波速度；W0为当前所用空间离散格式每个波长所需的最小采样点数；Δh为最大网格间距；f0为震源中心频率。

ADE－PML的主要思想为将1/sx分解为2项之和，并使得iω只存在于其中一项的分母中：

则式（4）转化为

其中，为与空间偏导数∂xvx相关的辅助变量，定义为

将式（12）改写为

将式（11）代入式（5），由频率域变换回时间域后，得到非分裂形式的PML方程（同理处理项）：

可见要得到时域PML方程，只需将空间偏导项∂xA替换为∂xA/κx＋ψx即可。将式（13）变换回时间域，得到关于辅助变量的微分方程：

同理可得关于的辅助微分方程。ADE－PML的特点在于通过直接在时间域离散微分方程式（15）来求解辅助变量，因此适用于任意时间精度的递推格式的数值模拟。下面以四阶Runge－Kutta（以下简称RK4）格式为例说明ADE－PML在高阶时间精度声波方程数值模拟中的应用。

在时间域直接对式（15）进行RK4时间离散，辅助变量满足其中：Δt为离散时间步长；下标i对应每个RK4循环中第i步迭代；当θ取0，1/2和1时，分别对应求取辅助变量的显式、半隐式和隐式格式。对式（16）化简可得

其中：

下面给出求解一阶速度－压力声波方程式（1）的RK4时间递推算法。设w＝（vx，vz，p），则由n时刻的wn更新n＋1时刻的wn＋1的过程可表达为

其中：算子L（·）表示式（1）中等式的右端项；RK4系数α2＝1/2，α3＝1/2，α4＝1；β1＝1/6，β2＝2/6，β3＝2/6，β4＝1/6。

2 数值模拟

为验证ADE－PML的吸收效果，使用四阶精度旋转交错网格［15－16］进行空间离散。模型计算区域为一5 200m×1 600m的长方形区域，由651×201个网格点构成，有限差分元胞边长为Δx＝Δz＝8 m。介质密度ρ＝2 000kg/m3，声波速度cp＝3 000 m/s。时间上采用RK4高阶时间更新格式，这里取Δt＝1.2ms，模拟持续时长3.6s，共计3 000时间步。采用点声压震源激发，震源子波为一导数高斯函数，其中a＝（πf0）2，震源中心频率f0＝15Hz，子波时延t0＝1.2/f0＝0.06s，激发点位于（176m，496m）处；设置2个检波点，检波点1和检波点2分别位于（1 456m，496 m）和（176m，4 896m）处；计算区域四周由PML吸收边界包围，厚度为80m，占10个有限差分元胞。震源和检波点2靠近左侧PML，检波点1靠近右侧PML。

为检验复频移拉伸算子的效果，分别给出使用传统拉伸算子（5）的ADE－PML（这里称为非复频移ADE－PML）和复频移拉伸算子（6）的 ADE－PML条件下的模拟计算结果。

该算例中，由于震源靠近左侧吸收边界，形成了高角度入射PML的情况。从波场快照图1a可见，非复频移ADE－PML对高角度入射波衰减不完全，入射波在左侧PML内形成损耗波，这些波产生了低频虚假反射，且虚假反射能量随着入射角度和传播距离的增大而增强：从1.68s的波场可见，在远偏移距低频虚假反射明显影响了求解区域的波场。而图1b中复频移ADE－PML边界由于使用复频移拉伸算子，对高角度入射波吸收很好，在左侧PML内没有产生损耗波和低频反射，在图1b中1.68s波场中，即使在远偏移距也没有可见的虚假反射。

在地震记录图2a中，随着偏移距的增大，高角度入射产生的虚假反射振幅逐渐变强，在直达波同相轴之后出现了一些“拖尾”现象，这是由于波场快照图1中在左侧吸收边界附近产生的低频虚假反射所致，从图1中1.68s的波场可见，在远偏移距这些低频反射尤为强烈。而图2中复频移ADE－PML由于使用了复频移拉伸算子，对高角度入射波吸收很好，没有出现图1中的低频虚假反射；在图2中1.68s波场中，即使在远偏移距也没有产生可见的虚假反射，在地震记录图3b中，同非复频移ADEPML相比（图3a），复频移ADE－PML条件下仅可见直达波同相轴而没有产生图3a中的“拖尾”现象。

从检波点1地震记录（图3a）可见，在检波点1处，非复频移ADE－PML和复频移ADE－PML的计算结果都与精确解（即解析解）吻合良好，说明在垂直或小角度入射情况下，使用非复频移和复频移拉伸算子都可以高效吸收外行波场能量。因此，复频移拉伸算子对垂直入射波的吸收无明显改善。

而从检波点2的地震记录（图3b）中可见：非复频移ADE－PML条件下的地震记录上产生了明显的虚假反射抖动，与参考解吻合不佳，这是因为检波器2记录到了高角度入射波产生的低频虚假反射，这些虚假能量随着偏移距增大而增强；而在检波器2处，复频移ADE－PML情况下的记录与参考解吻合良好。这说明，复频移拉伸算子对切入射波的吸收有明显改善。

图1 非复频移（a）和复频移（b）ADE－PML边界条件下的声压场p的波场快照Fig.1 Snapshots of the pressure field pof acoustic wave in the non－CFS（a）and CFS（b）ADE－PML condition

图2 非复频移（a）和复频移（b）ADE－PML条件下在直线x＝xs上接收到的地震记录Fig.2 Seismograms recorded at line x＝xsin the non－CFS（a）and CFS（b）ADE－PML condition

图3 非复频移和复频移ADE－PML边界条件下得到的检波点1（a）和检波点2（b）的声波方程地震记录同精确解的对比Fig.3 Seismograms calculated in non－CFS ADE－PML and CFS ADE－PML conditions at receiver 1（a）and receiver 2（b）

图4 显式、隐式和半隐式更新ADE－PML辅助变量条件下在检波器1（a）和检波器2（b）处得到的地震记录Fig.4 Seismograms calculated using the explicit，implicit and semi－implicit update scheme for the ADE－PML auxiliary variables at receiver 1（a）and receiver 2（b）

图5 显式、隐式和半隐式更新ADE－PML辅助变量条件下在检波器1（a）和检波器2（b）处得到的测试解同精确解的差异Fig.5 Differences between the test and reference solution calculated using the explicit，implicit and semi－implicit update scheme for the ADE－PML auxiliary variables at receiver 1（a）and receiver 2（b）

图6 波场总能量衰减Fig.6 Total energy decay of the wavefield

如前文式（15）中所示，对于ADE－PML辅助变量的计算可使用显式、隐式、和半隐式更新格式，这里分别对3种辅助变量的更新格式下的ADE－PML吸收效果进行了检验。如前述复频移拉伸算子可以改善吸收效果，因此在复频移ADE－PML边界条件下计算。从图4中可见，3种更新格式得到的地震记录均无虚假反射振动，因为复频移拉伸算子的作用，在检波点2处也与精确解吻合良好。图5中绘出了3种更新格式同精确解的差异曲线，用于细微比较三者的差异。从图中差异曲线的摆动幅度可见，由于数值频散效应在远偏移矩的积累，图5b的差异曲线比图5a中的摆动大，无论是垂直入射的检波点1还是高角度入射的检波点2处，显式格式比隐式和半隐式格式同精确解的吻合程度好。

图6为3种更新格式在长时间数值模拟中（这里进行了2×105时间步，共240s的模拟计算）的波场总能量变化情况，从波场能量的变化情况可检验 ADE－PML吸收边界的稳定性［7，13］。求解区域Ω中（不包括吸收边界区域）的总能量E通过下式计算：

从图6a可见：当数值模拟开始时（约0.1s左右），震源将能量注入介质使得总能量急剧上升；之后由于吸收边界作用，波场能量开始逐渐下降。在约1.7 s时能量出现一个陡降，因为这时声波完全离开求解区域，理论上讲该时刻之后求解区域中的能量为虚假能量。由图6b可见，在数值模拟后期（6s之后），波场总能量持续下降至10－20J数量级，说明显式、隐式和半隐式记忆变量更新格式都具有长时间稳定性，并且仔细观察可发现3种更新格式的能量衰减速度存在差异，隐式情况下能量衰减最快。

3 结论和讨论

数值模拟表明，使用复频移拉伸坐标算子的ADE－PML比使用非复频移拉伸算子（即传统拉伸坐标算子）的ADE－PML更能有效消除高角度入射情况下的虚假反射。对ADE－PML辅助变量的3种更新格式的检验说明，3种格式在吸收效果上不存在明显差异，但显式条件下可以获得相对较好的吸收效果；长时间的能量衰减计算表明，3种更新格式条件下的ADE－PML都稳定至2×105时间步，说明ADE－PML在数值模拟中具有长时间稳定性。由于ADE－PML技术具有非分裂场格式，降低了编程难度，且由于适应高阶精度时间格式，可以进一步推广到波动方程高精度偏移成像中。

（Reference）：

［1］Clayton R，Engquist B.Absorbing Boundary Conditions for Acoustic and Elastic Wave Equations［J］.Bulletin of the Seismological Society of America，1977，67（6）：1529－1540.

［2］Liao Z P，Wong H L，Yang B P，et al.A Transmitting Boundary for Transient Wave Analysis［J］.Scientia Sinica，1984，27（10）：1063－1076.

［3］Higdon R.Absorbing Boundary Conditions for Difference Approximations to the Multidimensional Wave Equation［J］.Mathematics of Computation，1986，47（176）：437－459.

［4］Bérenger J P.A Perfectly Matched Layer for the Absorption of Electromagnetic Waves［J］.Journal of Computational Physics，1994，114（2）：185－200.

［5］Chew W C，Weedon W H.A 3－D Perfectly Matched Medium from Modified Maxwell’s Equations with Stretched Coordinates［J］.Microwave and Optical Technology Letters，1994，7（13）：599－604.

［6］Drossaert F H，Giannopoulos A.A Nonsplit Complex Frequency Shifted PML Based on Recursive Integration for FDTD Modeling of Elastic Waves［J］.Geophysics，2007，72（2）：9－17.

［7］Komatitsch D，Martin R.An Unsplit Convolutional Perfectly Matched Layer Improved at Grazing Incidence for the Seismic Wave Equation［J］.Geophysics，2007，72（5）：155－167.

［8］Ramadan O.Auxiliary Differential Equation Formulation：An Efficient Implementation of the Perfectly Matched Layer［J］.IEEE Microwave and Wireless Components Letters，2003，13（2）：69－71.

［9］Moerloose J D，Stuchly M A.Behavior of Berenger’s ABC for Evanescent Waves［J］.IEEE Microwave and Wireless Components Letters，1995，5（10）：344－346.

［10］Moerloose J D，Stuchly M A.Reflection Analysis of PML ABC’s for Low－Frequency Applications［J］.IEEE Microwave and Guided Wave Letters，1996，6（4）：177－179.

［11］Kuzuoglu M，Mittra R.Frequency Dependence of the Constitutive Parameters of Causal Perfectly Matched Anisotropic Absorbers［J］.IEEE Microwave and Guided Wave Letters，1996，6（12）：447－449.

［12］Collino F，Tsogka C.Application of the PML Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media［J］.Geophysics，2001，66（1）：294－307.

［13］Festa G，Vilotte J P.The Newmark Scheme as Velocity－Stress Time－Staggering：An Efficient PML Implementation for Spectral－Element Simulations of Elastodynamics［J］.Geophysical Journal International，2005，161（3）：789－812.

［14］Zhang W，Shen Y.Unsplit Complex Frequency－Shifted PML Implementation Using Auxiliary Differential Equations for Seismic Wave Modeling［J］.Geophysics，2010，75（4）：T141－T154.

［15］周辉，谢春临，王尚旭，等.复杂地表地震资料叠前深度偏移成像［J］.吉林大学学报：地球科学版，2012，42（1）：262－268.Zhou Hui，Xie Chunlin，Wang Shangxu，et al.Prestack Depth Migration for Geological Structures with Complicated Surface［J］.Journal of Jilin University：Earth Science Edition，2012，42（1）：262－268.

［16］孟庆生，樊玉清，张珂，等.高阶有限差分法管波传播数值模拟［J］.吉林大学学报：地球科学版，2011，41（1）：292－298.Meng Qingsheng，Fan Yuqing，Zhang Ke，et al.Tube Wave Propagation Numerical Simulation Based on High Order Finite－Difference Method［J］.Journal of Jilin University：Earth Science Edition，2011，41（1）：292－298.