理论反射系数对于完全匹配层边界的吸收效果影响分析❋

张晓波， 宋　鹏,2,3❋❋， 李金山,2,3， 谭　军,2,3， 刘照仑， 夏冬明,2,3， 姜秀萍,2,， 赵　波,2

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



理论反射系数对于完全匹配层边界的吸收效果影响分析❋

张晓波1， 宋鹏1,2,3❋❋， 李金山1,2,3， 谭军1,2,3， 刘照仑1， 夏冬明1,2,3， 姜秀萍1,2,， 赵波1,2

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

摘要：本文基于平面波理论和空间二阶精度差分格式推导出声波方程数值模拟中的完全匹配层边界条件反射系数公式，并基于该公式系统分析了理论反射系数对于完全匹配层吸收效果的影响。实验结果表明，对于小入射角度(小于30°)的入射波，当理论反射系数R取10-4时，完全匹配层有最优的吸收效果；对于入射角度为40°、50°和60°的入射波，当理论反射系数R分别取10-5、10-6和10-8时完全匹配层边界反射系数最小；而对更大角度的入射波，需要更小的理论反射系数R才能使得完全匹配层达到最优的吸收效果。

关键词：声波方程；有限差分数值模拟；完全匹配层；理论反射系数

引用格式：张晓波，宋鹏，李金山，等.理论反射系数对于完全匹配层边界的吸收效果影响分析[J]. 中国海洋大学学报(自然科学版)， 2016,46(6)： 84-89.

ZHANG Xiao-Bo,SONG Peng,LI Jin-Shan,et al. Analysis of the impactof the theoretical reflection coefficient on the perfectly matched layer absorbing boundary[J]. Periodical of Ocean University of China, 2016，46(6)： 84-89.

数值模拟是研究地震波传播规律的有效手段，在地震勘探的采集、处理以及反演等环节中有着广泛应用。地震波数值模拟中，必须用有限的计算空间来模拟实际的无限空间，这就引入了人工边界问题。如果人工边界处理不当，不仅会产生污染中心波场的虚假反射，而且会使地震波的运动学和动力学特征产生畸变，影响模拟精度，因此，人工边界的正确处理是地震波高精度数值模拟的前提。

目前常用的人工边界处理方法主要有吸收边界条件(ABC)方法和完全匹配层(PML)边界条件两种。ABC方法[1-5]主要通过在边界处应用单程外行波方程来消除边界反射，由于单程波方程往往具有小角度入射的假设，这种方法很难对大角度入射波取得满意的吸收效果。PML方法主要通过在中心波场计算区域外加入吸收层，并在吸收层内引入衰减因子来达到消除边界反射的目的，理论上PML方法不存在对地震波入射角和频率的假设，其具有更好的吸收效果，因此PML方法已成为目前应用最为广泛的边界处理手段[6-11]。

对于PML边界条件来说，当吸收层数相同时，其吸收效果取决于PML中的吸收衰减因子和理论反射系数，衰减因子的设定及其特性已有大量文献发表[12-16]，本文不再赘述。在PML理论反射系数的研究方面，相关文献或讨论很少，目前较为常用的方法或思路是：在边界处理过程中直接给定一个相对较小的常数作为理论反射系数，然后依据波场快照和合成记录判断PML吸收效果。事实上理论反射系数的设置对于PML的吸收效果有重要影响，该系数的设定需要理论依据作为支撑，而这恰好是当前国内外在PML研究领域的盲区。

本文基于平面波理论并采用二阶差分精度格式，推导出了离散情况下频率域PML的边界反射系数的计算公式，并以声波方程数值模拟为例，分析了理论反射系数对于PML吸收效果的影响规律。

1一阶速度-应力声波方程及其PML边界处理

二维各向同性介质中的声波方程为[13]：

(1)

数值求解式(1)时，常将模拟空间划分为3个区域：中心波场区域、PML区域和角点区域(见图1)。

PML区域的控制方程为[13]：

(2)

其中：Px和Pz分别为P在x和z方向上的分量；dx和dz分别为x和z方向上的衰减因子；其常用表达式为：

其中：m为PML层内的计算点到中心波场边界的距离；R为理论反射系数；L是PML层厚度。实际计算时，在中心波场区域，dx=0，dz=0。在PML区域1，dx=0，dz=d(m)。在PML区域2，dx=d(m)，dz=0；在角点区域dx=d(m)，dz=d(m)。本文主要以式(3)为基础讨论理论反射系数R对于PML吸收效果的影响。

2离散条件下频率域PML的边界反射系数计算公式推导

式(1)在均匀介质中有如下形式的解[12]：

(4)

其中：c是纵波波速；θ为入射角；u0为波动幅值；t为时间；ω为频率；x和y表示空间坐标。以x方向的吸收衰减为例，PML层中声波方程的解[12]为

(5)

(6)

其中：i为虚数单位；x0为x方向上与PML层邻接的中心波场边界的位置。将式(6)代入式(5)中，然后结合式(4)可得PML层与中心波场的振幅比为

(7)

当波在PML层中传播并反射回来之后，可以很容易地得到经PML层吸收衰减之后的反射系数为

(8)

将吸收衰减因子d(m)(以(3)式为例)代入式(8)得：

(9)

从式(9)中可以看出，在介质连续的情况下，反射系数仅与入射角度和理论反射系数的值有关，当垂直入射时，反射系数Rm与理论反射系数R相等，并且对于任意角度的入射波，R的取值越小，PML边界的反射系数越小，即PML边界的吸收效果最好。但在实际模拟时，由于有限差分需要对空间进行离散化采样，空间离散化会对反射系数产生影响，因此当理论反射系数不断减小时，实际反射系数并不随之减小[17]。若要深入分析实际模拟时理论反射系数对于PML边界的吸收效果影响，需推导出离散条件下的PML边界反射系数公式。

因为PML吸收边界条件在各个方向上的吸收原理具有一致性，因此本文以x方向的吸收衰减为例推导离散条件下的PML边界反射系数公式。令式(2)中的衰减因子dx=d(m)，dz=0，则PML控制方程写为

(10)

(11)

(12)

为了构造一种非分裂的PML表达式，这里引入Bérenger变量[18]

(13)

将其代入式(12)中得

(14)

(15)

令Δx=Δz=h，将上式写成空间二阶精度有限差分的格式

(16)

(17)

为了方便计算，这里将Pl,j写作Pl,j=Ple-i(kzjh)，其中kz=ksinθ，这里k表示波数，θ表示入射角。然后将其代入式(17)中可得

(18)

(19)

(20)

1)当l=-1时，γ-1=γ-2=1，则有

(21)

2)当l=0时，γ-1=1，则有

(22)

3)当l=nl-2时，则有

(23)

(24)

式中

利用上述方程组，代入式(25)所示的反射系数的表达式[18]求出反射系数，从而得到网格离散情况下理论反射系数与实际反射系数的关系。

(25)

其中：kx=kcosθ，这里k表示波数；θ表示入射角。

3理论反射系数对于PML边界吸收效果影响分析

以均匀介质模型(速度为2 500m/s)为例研究理论反射系数与实际反射系数的关系，模型网格大小为2 m×2 m，计算时频率取为15 Hz，PML层数由10层并以10层为间隔依次增至50层，并且将理论反射系数由10-1并以10-1为比例因子依次降至10-15，得到入射角分别为0°～70°(间隔为10°)情况下的理论反射系数与PML边界反射系数的关系图线(关系曲线如图2所示，其中红色实线、黑色实线、蓝色实线、红色虚线和蓝色虚线分别表示PML层数为10、20、30、40和50层时的理论反射系数与PML边界反射系数关系曲线，绿色虚线指示最小实际反射系数所对应的理论反射系数值)。

由图2可以看出，当PML层数为10、20、30、40或50层时，由于空间的离散化采样，并不是理论反射系数R的取值越小，PML边界吸收效果越好，而是存在一个最优的R值使得PML边界的反射系数达到极小值，并且不同的入射角度其最优的R值不同，对于小入射角度(小于30°)的入射波，当理论反射系数R取10-4时有最优的吸收效果，对于入射角度为40°、50°和60°的入射波，当理论反射系数R分别取10-5、10-6和10-8时其边界反射系数最小，而对于更大角度的入射波，需要更小的理论反射系数才能达到最优的边界吸收效果。

4结论

本文在详细论述PML基本原理的基础上，基于平面波理论和空间二阶精度差分格式推导出声波方程数值模拟中离散条件下的完全匹配层边界反射系数公式，并基于该公式系统分析了理论反射系数对于PML边界吸收效果的影响。理论分析和实验结果得出结论如下：

(1)基于声波方程实际数值模拟时，无论PML层数为10、20、30、40或50层，由于空间的离散化采样，并不是理论反射系数R的取值越小，PML边界吸收效果越好，而是存在一个最优的R值可以使PML边界的实际反射系数值最小。

(2)对不同的入射角度其最优的理论反射系数R值不同，对于小入射角度(小于30°)的入射波，当理论反射系数R取10-4时有最优的吸收效果，对于入射角度为40°、50°和60°的入射波，当理论反射系数R分别取10-5、10-6和10-8时其边界反射系数最小，而对更大角度的入射波，需要更小的理论反射系数R才能达到最优的边界吸收效果。

参考文献：

[1]Engquist B， Majda A． Absorbing boundary conditions for numerical simulation of waves [J]． Proceedings of the National Academy of Sciences， 1977， 74(5)： 1765-1766．

[2]Engquist B， Majda A． Radiation boundary conditions for acoustic and elastic wave calculations [J]． Communications on Pure and Applied Mathematics， 1979， 32(3)： 313-357．

[3]Higdon R L． Absorbing boundary conditions for difference approximations to the multidimensional wave equation [J]． Mathematics of Computation， 1986， 47(176)： 437-459．

[4]HigdonR L. Numerical absorbing boundary conditions for the wave equation [J]. Mathematics of Computation, 1987, 49(179): 65-90.

[5]宋鹏， 王修田． 优化系数的四阶吸收边界条件 [J]． 中国海洋大学学报(自然科学版)， 2008， 38(2)： 251-258．

Song Peng， Wang Xiutian． The fourth-order absorbing boundary conditions with optimal parameters [J]． Periodical of Ocean University of China， 2008， 38(2)： 251-258．

[6]Bérenger J P． A perfectly matched layer for the absorption of electromagnetic waves [J]． Journal of Computational Physics， 1994， 114(1)： 185-200．

[7]Chew W C， Liu Q H． Perfectly matched layers for elastodynamics: A new absorbing boundary condition [J]． Journal of Computational Acoustics， 1996， 4(4)： 341-359．

[8]Wang H X． A Study of damping factors in perfectly matched layers for the numerical simulation of seismic waves [J]． Applied Geophysics， 2013， 10(1)： 63-77．

[9]韩令贺， 何兵寿， 张会星． 含直立裂缝粘弹性介质地震波场正演模拟 [J]． 山东科技大学学报(自然科学版)， 2010， 29(1)： 17-25．

Han Linghe， He Bingshou， Zhang Huixing． Simulation of seismic wave field in viscoelastic media with vertical fractures [J]． Journal of Shandong University of Science and Technology， 2010， 29(1)： 17-25．

[10]Hastings F D， Schneider J B， Broschat S L． Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation[J]. The Journal of the Acoustical Society of America， 1996， 100： 3061-3069．

[11]何兵寿， 张会星， 范国苗． 双程声波方程逆时深度偏移 [J]． 物探化探计算技术， 2010， 32(2)： 113-119， 107．

He Bingshou， Zhang Huixing，Fan Guomiao． Reverse-time depth migration of two way acoustic wave equations [J]． Computing Techniques for Geophysical and Geochemical Exploration，2010，32(2)： 113-119，107．

[12]Collino F and Tsogka C． Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media [J]． Geophysics，2001，66(1)： 294-307．

[13]王守东． 声波方程完全匹配层吸收边界 [J]． 石油地球物理勘探，2003，38(1)： 31-34．

Wang Shoudong． Absorbing boundary condition for acoustic waveequation by perfectly matched layer [J]． Oil Geophysical Prospecting，2003，38(1)： 31-34．

[14]陈可洋． 完全匹配层吸收边界条件研究 [J]． 石油物探， 2010， 49(5)： 473-477．

Chen Keyang． Study on perfectly matched layer absorbing boundary condition [J]． Geophysical Prospecting for Petroleum， 2010，49(5)： 473-477．

[15]李飞，程久龙，杨思通． 基于Cerjan衰减函数的PML吸收边界条件 [J]． 物探与化探， 2012， 36(6)： 996-1000．

Li Fei，Cheng Jiulong，Yang Sitong． A study of PML absorbing boundary based on Cerjian attenuation function [J]． Geophysical & Geochemical Exploration， 2012， 36(6)： 996-1000．

[16]李宁． 完美匹配层理论及其在地震波模拟中的应用 [D]． 北京： 中国地震局工程力学研究所，2006．

Li Ning． The theory and the application of PerfectlyMatched Layer in the seismic wave simulation [D]． Beijing： Institute of Engineering Mechanics，China Earthquake Administration，2006．

[17]Wu Z，Fang J． Numerical implementation and performance of perfectly matched layer boundary condition for waveguide structures [J]． IEEE Transactions on Microwave Theory and Techniques，1995， 43(12)： 2676-2683．

[18]Collino F，Monk P B． Optimizing the perfectly matched layer [J]． Computer Methods in Applied Mechanics and Engineering，1998，164(1)： 157-171．

责任编辑徐环

Analysis of the Impact of the Theoretical Reflection Coefficient on the Perfectly Matched Layer Absorbing Boundary

ZHANG Xiao-Bo1, SONG Peng1,2,3, LI Jin-Shan1,2,3, TAN Jun1,2,3,LIU Zhao-Lun1, XIA Dong-Ming1,2,3, JIANG Xiu-Ping1,2, ZHAO Bo1,2

(1.College of Marine Geo-science, Ocean University of China, Qingdao 266100, China; 2.Key Lab of Submarine Geosciences and Prospecting Techniques Ministry of Education, Qingdao 266100, China; 3.Qingdao National Laboratory for Marine Science and Technology, Qingdao 266100, China)

Abstract:This paper derives the reflection coefficient formula of the perfectly matched layers (denoted as “PML”)for the numerical simulation of the acoustic wave equation, based on the plane wave theory and the second-order accuracy spatial difference scheme. Then, the impact of the theoretical reflection coefficientRon the absorption effect of the PML is analyzed. The experimental results demonstrate that, for the incident wave with a small incident angle (less than 30°), whenRis equal to 10-4, the PML absorbing boundary condition has the optimal absorption effect.And for the incident wave with the incident angles of 40°, 50° and 60°, when theRis equal to 10-5, 10-6and 10-8respectively, the PML absorbing boundary condition has the smallestactual reflection coefficients.While for the incident wave with a large incident angle, in order to achieve the optimal absorption effect, it is necessary to take a smallerR.

Key words:acoustic wave equation; finite difference numerical simulation; perfectly matched layers; theoretical reflection coefficient

基金项目：❋ 国家自然科学基金项目(41574105)资助

收稿日期：2015-05-26；修改日期：2015-06-14

作者简介：张晓波(1990-)，男，博士生，主要从事地震波正演模拟和属性分析研究。E-mail：xiaobozhang@outlook.com ❋❋通讯作者：E-mail：pengs@ouc.edu.cn

中图法分类号：P631.4

文献标志码：A

文章编号：1672-5174(2006)06-084-06

DOI：10.16441/j.cnki.hdxb.20150190

Supported by the National Natural Science Foundation of China (41574105)