邢浩洁 刘爱文 李小军 陈 苏 傅 磊

1） 中国北京 100081 中国地震局地球物理研究所

2） 中国北京 100124 北京工业大学城市建设学部

引言

人工边界条件是地震波动数值模拟技术的关键环节（廖振鹏，2002），是计算人工边界节点运动的各种方法的统称，其任务在于尽可能地“算准”外行波引起的边界点运动，以达到消除不合理边界虚假反射的目的.半个世纪以来，国内外研究人员提出了基于各种直接或间接方法适用于不同波动问题的人工边界条件，其中我国学者廖振鹏及其合作者所发展的多次透射公式（multi-transmitting formula，缩写为MTF），本文称之为廖氏透射人工边界，具有形式简单、易于实现、通用性好的优点（廖振鹏等，1984；Liao，Wong，1984）.该边界具有清晰明确的物理机制，即使用统一人工波速进行“误差波多次透射”，因而受到研究人员和工程师的青睐.

三十多年来，MTF 人工边界在场地地震反应和土-结相互作用数值计算中发挥了重要作用（李小军等，1995；Shiet al，2016；Huang，2018；陈少林等，2020），关于MTF 精度优化（廖振鹏，李小军，1995；Liao，1996）及数值稳定性问题（Liao，Liu，1992；李小军，廖振鹏，1996；廖振鹏等，2002；景立平，2004；章旭斌等，2015；Xie，Zhang，2017）的研究也随之不断深入.近几年，MTF 开始被用于基于谱元法的地震波动模拟（戴志军等，2015；邢浩洁，李鸿晶，2017a，b；于彦彦等，2017；Xinget al，2021a）.Xing 等（2021a）对二维SH 波动和P−SV 波动的研究结果表明MTF 边界的精度略低于完美匹配层（perfectly matched layer，缩写为PML）边界（李宁等，2007；Liuet al，2014；谢志南，章旭斌，2017），但要优于黏弹性边界（刘晶波，李彬，2005；杜修力等，2006）.最近，作者（邢浩洁等，2021；Xinget al，2021b）成功地将多个人工波速引入MTF 表达式当中，提出了一种多人工波速优化透射边界，记为caj-MTF，其中caj（j＝1，2，···，N）表示与边界阶次N相对应的多个人工波速.该边界所具有的多个人工波速配置能够与弹性波、两相介质波动、频散波动等问题中的多种物理波速一一对应（即各个人工波速参数可以分别取为各个物理波速），解决了传统廖氏透射边界中单一人工波速不便于选取的难题，使得复杂波动模拟中的边界性能得到显著提升.

与已有的廖氏透射边界优化形式（李小军，廖振鹏，1996；Liao，1996）相比，多人工波速透射边界与MTF 表达式更为接近，完全继承了MTF 简洁、精度可控和普遍适用性的优点.廖振鹏（2002）曾经提出一种含有多个人工波速的MTF 边界的推广形式，其解析表达式为

式中，u（s，t）为人工边界所在的指向内域的离散网格线上的波场，N为边界阶次，I为单位算子，B（caj）为表示时间-空间移动的离散算子，caj（j＝1，2，···，N）为一组可调的人工波速.离散算子B（caj）对波场u（s，t）的作用是将其转换为另一个时空节点（即参考点）处的波场，该点在时间上向前平移一个时刻Δt，空间上沿着人工边界节点所在的离散网格线向计算区域内部平移cajΔt距离；算子的乘积如B（ca1）B（ca2）表示时间-空间平移（2Δt，ca1Δt＋ca2Δt）距离，依此类推.当caj-MTF 边界中各个人工波速取相同数值时，即退化为MTF 边界的算子表达式［I−B（ca） ］Nu（s，t）＝0 （廖振鹏等，2002；景立平，2004）.邢浩洁等（2021）证明了这两个边界表达式具有等价性，虽然二者的表达形式和数值离散格式都不一样，但在波动数值模拟中的精度和稳定性几乎完全相同.不过，离散边界caj-MTF 比解析边界更易于与不同的内域离散方法相结合，前者只需将各个离散参考点的运动由临近的内域节点的运动进行插值表示即可.作者在前期工作中对caj-MTF 边界模拟标量波和弹性波的性能进行了理论分析，并且在有限差分和有限元波动数值模拟中予以验证，本文将进一步探讨该边界在谱元法地震波动模拟中的应用.

本文将首先介绍多人工波速优化透射边界caj-MTF 的表达式，阐明该边界提高边界性能的原理；然后给出在该边界谱元离散网格下的具体数值格式，并指出其与有限差分或有限元格式的差异；最后，采用简单雷克（Ricker）子波激励下的体波大角度透射和面波透射数值算例，研究该边界对不同外行波的反射率，并与传统MTF 边界、完美匹配层边界、黏弹性边界、一阶旁轴近似边界（Clayton，Engquist，1977）及自由边界模拟结果进行比较，以期获得关于谱元法地震波动模拟中不同人工边界基本性能的定量认识.

1 基于多人工波速的优化透射边界

多人工波速优化透射边界（caj-MTF）的具体原理详见Xing 等（2021b）和邢浩洁等（2021），这里仅简要介绍.常用的前三阶优化边界表达式为

图1 为caj-MTF 边界示意图.每个人工边界节点0 的运动都由一系列参考点（图中空心圆圈）的运动推算得到，这些参考点位于一条指向内域的离散网格线上（图中以局部坐标s表示）且分别处于不同时刻（t－Δt，t－2Δt，t－3Δt，···）.参考点的空间坐标ca1Δt，ca2Δt，ca1Δt＋ca2Δt，···由各个人工波速ca1，ca2，ca3和时间步长Δt共同决定.观察图1 并比较上述两种边界表达式可知，caj-MTF 边界对MTF 边界的优化是通过将各个时刻的单一参考点替换为多个参考点来实现的.邢浩洁等（2021）将这种改进解释为在保留“误差波多次透射”物理机制的同时对各阶误差波的表述加以优化.

图1 多人工波速优化透射边界（ caj-MTF）示意图（s0t 为定义人工边界条件的局部坐标系）Fig. 1 Schematic diagram of optimized transmitting boundary with multiple artificial wave velocities （ caj-MTF）（s0t is a local coordinate system that is used to define the artificial boundary condition）

2 优化透射边界的谱元计算格式

由于caj-MTF 边界具有与MTF 边界非常相似的形式，其数值实现方式与后者也基本一致（邢浩洁，李鸿晶，2017a）.如图2 所示，利用紧邻人工边界的谱单元节点来实现caj-MTF.对于任意人工边界节点0，要求该节点位于一条指向内域的离散网格线上，紧邻人工边界的谱单元在该离散网格线上的一组节点为0，1，···，M（M为单元插值函数的阶次）.在采用caj-MTF 计算（p＋1）Δt时刻边界节点0 的运动时，只需将式（3）中pΔt，（p－1）Δt，（p－2）Δt，···各个时刻的边界参考点（图中空心圆圈）的运动用谱单元节点0，1，···，M的运动插值表示即可.

图2 优化透射边界的谱元计算格式Fig. 2 Spectral-element computational scheme of the optimized transmitting boundary

根据上述插值方法得到优化透射边界caj-MTF 的谱元计算格式为

式中：u1，u2和u3表示单元节点运动组成的列向量；T1，T21，···，T33表示插值系数组成的行向量；sk（k＝0，1，···，M）表示从边界节点0 到编号为k的内域节点的距离；T31，T32和T33中的插值系数与T21，T22类似，为节省篇幅，不再具体列出.式（5）—（6）所表示的caj-MTF 边界的谱元计算格式与已有的MTF 边界的谱元计算格式（戴志军等，2015；邢浩洁，李鸿晶，2017a，b；于彦彦等，2017；Xinget al，2021a）非常相似，二者都简单地将各个边界参考点的运动用该方向上谱单元节点的运动插值表示.这种简单实用的施加方式为廖氏透射边界的重要特色，本文边界继承了这个优点.不过，本文边界也具有Xing 等（2021a）所指出的谱元法中MTF 数值格式的缺点，即它不能像传统MTF 有限差分或有限元格式那样由一阶边界格式的乘方来得到高阶边界格式，这使得谱元法中高阶透射边界的精度和稳定性受到数值离散误差的影响很大.因此，我们建议在谱元法地震波动模拟中使用MTF 或caj-MTF 边界时，边界阶次取N＝2.

3 波动模拟结果分析

3.1 计算模型和用于对比的几种经典人工边界

本文以谱元法模拟脉冲型地震动输入的波动模拟结果来检验多人工波速优化透射边界（caj-MTF）的实用性，并通过与几种经典人工边界以及自由边界（物理边界）进行比较分析，得到关于caj-MTF 边界误差量级及其相对优势的定量认识.采用如图3 所示的计算模型，分别考虑体波大角度透射问题和面波透射问题两种情形.

图3 人工边界性能分析模型（a） 体波大角度透射问题；（b） 面波透射问题Fig. 3 Analytical models of artificial boundary performance（a） The problem of body wave transmission over large incident angles；（b） Surface wave transmission problem

对于体波问题，波速设定为cP＝1 000 m/s，cS＝440 m/s，对应于泊松比ν＝0.38.模型尺寸为400 m×1 200 m，波源位于坐标（120 m，600 m），三个观察点位于左边界上，A（0，500 m），B（0，400 m）和C（0，200 m）.这样大的纵横波速比（或泊松比）和透射角度对边界性能提出了很高的要求.对于面波问题，波速cP＝1 000 m/s，cS＝577.4 m/s，对应的泊松比为ν＝0.25.模型尺寸为800 m×600 m，波源位于坐标（500 m，20 m），三个观察点位于右边界上，从地表向下依次为D（800 m，0），E（800 m，20 m）和F（800 m，80 m）.两个问题的波源均采用雷克子波，施加在节点的z方向位移uz之上，波源函数表达式为

式中，t0为脉冲起始时刻，f0为中心频率，与之对应的脉冲时间宽度为2/f0，对应的峰值频率fmax约为2.5f0.本文取t0＝0，f0＝10 Hz，则脉冲时间宽度为0.2 s，峰值频率约为25 Hz.输入波的时程及其傅里叶谱如图4 所示.

图4 输入波时程及其傅里叶幅值谱Fig. 4 Time history and Fourier spectrum of the incident wave

模型的空间离散采用25 节点正方形谱单元，单元尺寸为20 m×20 m.时间步长在满足内域时间积分稳定性的前提下，随着不同边界条件的影响进行一定调整.用于与本文caj-MTF 边界相比较的其它几种边界分别为：MTF 边界（Xinget al，2021a）、PML 边界（Zenget al，2011；Liuet al，2014）、黏弹性边界（杜修力等，2006）、一阶旁轴近似边界（Clayton，Engquist，1977）以及自由边界（物理边界，谱元法自动满足该条件）.采用扩展计算区域的数值解作为参考解，来衡量各个人工边界的反射误差.

3.2 体波大角度透射问题

不同边界条件模拟体波大角度透射问题所得到的z方向位移uz的波场快照如图5 所示，各子图的快照时刻分别为0.3 s，0.4 s，0.5 s，0.7 s.从图5 可以看出：自由边界这种物理边界会造成大量虚假反射，且边界反射存在波型转换（四种反射模式：P−P，P−S 和S−P，S−S）；所有人工边界都能够极大地降低不必要的虚假反射，其模拟结果远远优于自由边界，这体现出人工边界条件的重要性；PML 边界精度最高，已看不出边界反射；本文caj-MTF 边界的精度明显好于MTF、黏弹性边界和一阶旁轴近似边界，虽然不及PML 那样完美，但也达到了很好的模拟效果.由此可见，波场快照定性地反映了不同人工边界的性能差异.

图6 给出了观察点A，B和C处z方向位移uz的时程及其相对于参考解的误差曲线.为便于比较，图中未显示误差很大的自由边界结果以及与MTF 和黏弹性边界精度相当的一阶旁轴近似边界结果，后文表格将列出全部结果.通过观察点时程曲线可以看到，黏弹性边界结果总是小于参考解，表明该边界给计算模型提供的约束“偏刚”，另外几种边界结果则看不出明显规律.从误差曲线可知，除PML 之外的几种人工边界均产生了显著的误差.不过，caj-MTF 边界相对于MTF 边界和黏弹性边界具有明显优势，不仅其误差幅值大大小于后两者，还对大角度透射波（见观察点B和C的误差曲线）实现了较高精度的模拟.前文所指出的caj-MTF 边界能够较传统边界更好地适应含有不同物理波速或存在大角度透射（大角度对应于大的视传播速度）的复杂波动问题，在这里得到验证.

为定量地比较不同边界条件对体波模拟效果的差异，表1 详细列出了观察点A，B和C处波场分量ux和uz的数值解相对于参考解（扩展模型解）的误差曲线（uz的误差曲线见图6 右侧，ux的未画出）的峰值（绝对值），以及这些误差峰值相对于参考解时程峰值的百分比（对于每种人工边界条件，将各个观察点的结果进行平均）.

图6 观察点A （a），B （b）和C （c）处z 方向位移 uz的时程及误差曲线Fig. 6 Time histories and error curves of the displacement uz at the observation points A （a），B （b） and C （c）

由表1 可知，各个观察点的最大误差绝对值存在明显差异，从A点到B点再到C点，误差逐步降低，这体现了波传播过程中几何衰减所导致的幅值下降.不过，各个观察点的最大误差总体上随着不同边界条件呈现出一致的变化趋势，这个趋势最终体现在边界的峰值百分误差上.不同边界的峰值百分误差显示出边界精度由高到低依次为：PML 边界最高，caj-MTF 边界较高，MTF 边界、黏弹性边界和一阶旁轴近似边界中等，自由边界最低.各个峰值百分误差数值存在着数量级的差异，PML 边界低至3.2%，caj-MTF 边界为15.3%，MTF 边界、黏弹性边界和一阶旁轴近似边界分别为37.0%，62.1%，33.7%，自由边界则高达119.2%.由于峰值百分误差是对各个边界误差最大值的度量，而从图5 所示的波场快照来看，总体上人工边界对波动能量的反射比例要远小于这些峰值.因此可以得出结论，在对大纵横波速比和大角度透射的体波进行模拟时，除PML 边界具有近乎完美的性能之外，本文的caj-MTF 边界亦具有很高的精度.

表1 体波情形下A，B 和C 点误差曲线峰值（绝对值）以及每种边界条件下误差峰值与参考解峰值的百分比Table 1 The peaks （absolute value） of error curves of the points A，B，C，and the percentages between those error peaks and reference solution peaks on different boundary conditions：Body wave case

图5 不同边界条件模拟体波大角度透射问题的波场快照（时刻依次为0.3 s，0.4 s，0.5 s，0.7 s）（a） caj-MTF 边界；（b） MTF 边界；（c） PML 边界；（d） 黏弹性边界；（e） 一阶旁轴近似边界；（f） 自由边界（物理边界）Fig. 5 Wave field snapshots of body wave propagating over large incident angles obtained from different boundary conditions （time instants are 0.3 s，0.4 s，0.5 s，0.7 s）（a） caj-MTF boundary；（b） MTF boundary；（c） PML boundary；（d） Viscous-spring boundary；（e） The first-order paraxial-approximation boundary；（f） Free boundary （physical boundary）

3.3 面波透射问题

不同边界条件模拟面波透射问题所得到的z方向位移uz的波场快照如图7 所示，各子图的快照时刻分别为0.5 s，0.9 s，1.2 s.面波快照结果与上述体波情形相似，主要为：与自由边界存在的大量虚假反射相比，所有人工边界都能够消除绝大部分边界反射；PML 边界模拟效果近乎完美；本文caj-MTF 仅存在微量面波反射，精度很高；MTF、黏弹性边界和一阶旁轴近似边界则有着较为显著的面波反射，精度一般.

图7 不同边界条件模拟面波问题的波场快照（时刻依次为0.5 s，0.9 s，1.2 s）（a） caj -MTF 边界；（b） MTF 边界；（c） PML 边界；（d） 黏弹性边界；（e） 一阶旁轴近似边界；（f） 自由边界（物理边界）Fig. 7 Wave field snapshots of surface waves obtained by different boundary conditions（time instants are 0.5 s，0.9 s，1.2 s）（a） caj -MTF boundary；（b） MTF boundary；（c） PML boundary；（d） Viscous-spring boundary；（e） The first-order paraxial-approximation boundary；（f） Free boundary （physical boundary）

图8 为面波观察点D，E和F处uz的时程及其相对于参考解的误差曲线.同样地，为了便于比较，图中未给出误差很大的自由边界结果以及与MTF 和黏弹性边界精度相当的一阶旁轴近似边界结果，后文表格将全部列出.D点位于地表，E点和F点分别距离地表20 m和80 m.三个观察点处时程曲线的幅值逐步降低，表现出面波幅值从地表向下逐步衰减的规律；从误差曲线来看，三个观察点的结果均表明PML 边界误差最小，caj-MTF 有少量误差，MTF 边界和黏弹性边界的误差进一步增大.综上可见，面波问题的模拟结果再次证明caj-MTF 边界相对于MTF 和黏弹性边界具有非常明显的优势.

为定量比较不同边界条件对面波模拟效果的差异，表2 详细列出了观察点D，E和F处波场分量ux和uz的数值解相对于参考解（扩展模型解）的误差曲线（uz的误差曲线见图8 右侧，ux的未画出）的峰值（绝对值），以及这些误差峰值相对于参考解时程峰值的百分比（对于每种人工边界条件，将各个观察点的结果进行平均）.

图8 观察点D （a），E （b）和F （c）处z 方向位移 uz的时程及误差曲线Fig. 8 Time history and error curves of the displacement uz at the observation points D （a），E （b） and F （c）

由表2 可知，各个观察点的最大误差总体上随着不同边界条件呈现出一致的变化趋势，这种趋势最终体现在峰值百分误差上：PML 边界误差最小，caj-MTF 边界次之，MTF 边界、黏弹性边界、一阶旁轴近似边界误差中等，自由边界则误差极大.不同边界的面波反射误差存在数量级差异，PML 边界为4.8%，caj-MTF 边界为11.7%，MTF 边界、黏弹性边界、一阶旁轴近似边界分别为59.7%，30.4%，41.0%，自由边界高达425.4%.因此，面波模拟结果同样表明caj-MTF 边界是一种精度很高的人工边界，虽然未达到PML 边界那样近乎完美的精度，但是显著优于MTF 边界、黏弹性边界和一阶旁轴近似边界.

表2 面波情形下D，E 和F 点误差曲线峰值（绝对值）以及每种边界条件下误差峰值与参考解峰值的百分比Table 2 The peaks （absolute value） of error curves of the points D， E，F， and the percentages between those error peaks and reference solution peaks under different boundary conditions：Surface wave case

4 讨论与结论

在使用有限的数值离散模型研究无限域波动问题时，如何施加人工边界条件（即吸收边界条件）是尚未得到圆满解决的关键环节.本文将作者（邢浩洁等，2021；Xinget al，2021b）最近发展的多人工波速优化透射边界（caj-MTF）应用于高精度谱元法的地震波动模拟，并通过与现有的廖氏透射边界（MTF）、完美匹配层边界（PML）、黏弹性边界以及一阶旁轴近似边界进行比较分析，证明了这是一种便捷、高效的人工边界条件实现方法.

caj-MTF 边界对传统MTF 边界的计算波速参数进行优化，但保持了后者“误差波多次透射”的基本原理以及简单通用的表达形式，因此caj-MTF 边界像MTF 边界一样对于不同类型的波动问题和数值方法具有普适性.研究结果显示，在谱元法中使用caj-MTF 边界或传统MTF 边界时，二阶边界效果最好，而三阶及以上阶次的边界则会受到数值离散误差的影响反而降低精度，这是目前尚难以解决的问题.MTF 边界或传统MTF 边界时，二阶边界效果最好，而三阶及以上阶次的边界则会受到数值离散误差的影响反而降低精度，这是目前尚难以解决的问题.

与现有的PML 边界相比，caj-MTF 边界的精度相差不大，而在复杂程度、计算量以及通用性方面则有着巨大优势.与MTF 边界、黏弹性边界或一阶旁轴近似边界相比，caj-MTF 边界在精度方面具有显著优势，尤其是对于外行波以大角度射向人工边界和面波透射情形.

本文研究结果对于使用高精度谱元法进行地震波动模拟时人工边界条件的优化选取具有重要参考价值.虽然文中只给出了二维弹性波的数值算例，但caj-MTF 边界或MTF 边界的施加方式决定了它们可以无差别地应用于一维、二维及三维波动离散模型，而其它人工边界条件则不具备这样的简单通用性.下一步拟将该方法拓展应用到一些更具实际意义的复杂波动问题当中，如局部场地的地震波散射问题、饱和介质弹性波传播问题、低频震源问题等，以探讨提高边界精度的方法以及确保长持时计算稳定性所要采取的措施.