0 引言

在地震波动正反问题数值模拟中计入实际介质阻尼可有效提升实测波形记录的幅值、相位(到时)的模拟精度(Savage et al., 2010)；提升基于实测记录反演地下波速结构的精度及可靠性(Kurzmann et al., 2013; Groos et al., 2014).地震学常基于线性滞弹性阻尼理论(Aki and Richards, 2002)建立阻尼介质本构，其频域本构与弹性介质频域本构形式上完全一致，区别在于前者中弹性模量为与频率相关的复数，后者弹性模量为常实数.计入复弹性模量在频域波动数值模拟中易于实现(Stekl and Pratt, 1998)，但在时域数值模拟中却十分困难.因复弹性模量与频域应变乘积的傅里叶反变换由卷积给出，且就一般形式的复弹性模量而言，离散卷积计算需存储以往所有时刻的应变值.基于广义标准线性体模型建立阻尼介质本构可规避这一问题(Blanc et al., 2016).广义标准线性体频域本构中复弹性模量M(ω)为若干单极点有理分式之和：M(ω)= ，其中M_l、a_l、b_l为大于零的正实数，其时域本构易于在时域波动模拟中实现，因的傅里叶反变换为h(t)+(a_l-b_l)[e^-b_ltH(t)*h(t)](H(t)为单位阶跃函数)，且指数函数与任意函数的卷积可利用迭代计算得到(Luebbers and Hunsberger, 1992; Ramadan, 2003).在任意有限宽频带范围，广义标准线性体模型能准确拟合介质Q值、波速随频率的变化，为此在滞弹性时域波动数值模拟中得到广泛应用(Fichtner, 2011).

利用时域波动数值模拟方法：有限差分、有限体积、(间断)有限元，求解近场波动问题需引入人工边界(Artificial boundary)截断无限域以获取有限计算区，并在截断处设置无反射边界条件(Non-reflecting Boundary Condition: NRBC)使得从有限计算区向无限域行进的外行波能无反射地穿过人工边界并不再返回至计算区(廖振鹏，2002).由于滞弹性近场波动数值模拟研究起步较晚，NRBC的研究多局限于完全弹性介质(Givoli, 1991, 2004).对弱阻尼介质或对NRBC精度要求不高的情形，可直接采用弹性NRBC(王振宇和刘晶波，2004；Askan et al., 2007; Bielak et al., 2011)，但对强阻尼介质或对NRBC精度要求较高情形，将高精度弹性NRBC未加修正地应用于滞弹性情形会导致数值失稳或模拟精度降低等问题.

应用高精度NRBC可有效提升计算效率.Liu最早开展了滞弹性近场波动问题高精度NRBC—PML的研究，分别就导电媒介电磁波方程以及黏性介质声波方程(阻尼介质本构基于单个开尔文体建立)构建了对应的PML(Liu, 1997; Liu and Tao, 1997)，其研究工作表明PML能高效衰减从有限计算区辐射出来的外行波动且数值稳定性较好，清晰揭示了PML的应用前景.

PML是当前时域滞弹性近场波动数值模拟中应用最广泛的高精度NRBC，其等价于无限域波动方程的解析延拓(Chew et al., 1997; Collino and Monk, 1998)，可通过对无限域波动方程进行复坐标延拓得到.对速度-应力形式波动方程，滞弹性PML的相关研究工作包括：有限差分(Liu and Tao, 1997; Ely et al., 2008; Martin and Komatitsch, 2009); 有限元(Groby and Tsogka, 2006；Ma and Liu, 2006)；对位移-应力(应变)形式波动方程，其中应力(应变)作为独立存储变量，相关研究包括：有限元(Basu and Chopra 2004；Assimaki et al., 2012)；对位移形式波动方程，相关研究包括：有限差分(Moczo et al., 2014)；有限元(Ping et al., 2016).直接基于位移形式波动方程建立PML及相应的数值离散格式，能显著降低计算存储，提高计算效率，避免额外引入变量的边界条件难于设置等问题(Kreiss et al., 2002；Moczo et al., 2014)，但相比前两类PML而言，第三类PML中多极点有理分式与频域函数乘积的傅里叶反变换计算问题更为突出，需对不同重复极点个数进行分类处理，计算对应的傅里叶反变换.然而，当重复极点个数n大于等于3时，难以建立反变换中t^n-1*h(t)的稳定数值离散格式.这一问题在采用高阶PML以提高宽频带近场波动数值模拟精度的研究中十分突出(Connolly et al., 2014; Feng et al., 2017).为此，研究人员提出通过对b_l进行微调以避免多重极点(Ping et al., 2016；Xie et al., 2016)，但均没有给出具体准则.截至目前，还没有普适且简便的多极点有理分式与频域函数乘积的傅里叶反变换计算方案，多重极点的分类反变换极大增加了编程工作量及计算量，同时可能导致数值失稳.

针对上述问题，考虑位移形式波动方程及其高阶集中质量有限元离散，本文构建了适用于滞弹性近场波动勒让德谱元模拟的完美匹配层.首先，采用作者提出的基于复坐标延拓技术变换弱形式无限域波动方程以构建PML的方法(谢志南和章旭斌，2017)，建立了弱形式时域滞弹性PML，规避以往基于复坐标延拓技术独立延拓无限域内波动方程及界面条件可能导致的：延拓后所得场方程和界面条件匹配不合理而引发数值失稳、计算精度低下等问题(Duru and Kreiss, 2014).其次，结合PML计算精度对微调复坐标延拓函数中的参数不敏感的特点，给出了明确的参数微调准则，并建立了一种普适、简便的多极点有理分式与频域函数乘积傅里叶反变换的计算方法，极大简化了PML计算.基于该方法可实现任意高阶PML(Feng et al., 2017)以及多轴PML(Meza-Fajardo and Papageorgiou, 2008).最后，基于新构建滞弹性完美匹配层，结合高阶勒让德谱元建立了滞弹性近场波动问题离散方案.

文中第一节给出了无限域内滞弹性波动方程及基于广义标准线性体建立的阻尼介质时域本构.第二节详细阐述了弱形式滞弹性PML的构建.第三节给出了一种PML参数的微调方案及多极点卷积简便计算方法.在此基础上，将弱形式滞弹性PML与勒让德谱元结合建立滞弹性近场波动谱元模拟方案.第四节基于数值算例阐明滞弹性PML的精度、计算效率及新建方案的稳定性.最后给出本文工作小结.

1 无限域滞弹性波动方程

就近场波动问题而言，常假定无限域为层状构造(图 1)，在小应变前提假定下，无限域内滞弹性介质波动满足平衡方程，即：

(1)

ρ为介质密度，为x_i方向加速度.∂_j为对x_j的偏导算子，σ_ij为应力张量分量.介质分界面处满足位移与应力连续条件，即：

(2)

自由地表处应力为零，即：

(3)

n_j为界面外法线单位矢量分量.无穷远处位移、应力为零且满足无穷远辐射条件(廖振鹏，2002).

与弹性介质不同，滞弹性介质应力(应变)与应变(应力)历史相关，应力-应变关系满足Boltzmann叠加原理.就各向同性滞弹性介质而言，其频域本构可表示为

(4)

时域本构为

(5)

(6)

为σ_ij、ε_kl的傅里叶变换，ε_kl为应变张量分量，δ_ij为狄拉克张量分量.M^v(ω)为的傅里叶变换，v=κ, μ.ψ^κ、ψ^μ分别为膨胀应力与纯剪切应力松弛函数，满足因果律：

(7)

与松弛函数满足因果律对应，复模量M^v(ω)的实部与虚部满足Kramers-Krönig关系(Waters et al., 2000).将(5)、(6)式代入(1)式得到位移形式波动方程，即方程中只有位移为独立变量.

基于广义标准线性体模型(图 2)建立阻尼介质本构，可得：

(8)

N为标准线性体个数，下标l为标准线性体的编号，l=1, …, N.τ_σl^v、τ_εl^v为第l个SLS对应的应力松弛时间和应变蠕变时间.τ_σl^v=η_l^v/δM_l^v，τ_εl^v=τ_σl^v(M_Rl^v+δM_l^v)/M_Rl^v，M_Rl^v、δM_l^v为第l个标准线性体中的弹簧元件的弹性系数，η_l^v为第l个标准线性体中黏性阻尼元件的阻尼系数.M^v(ω)为

(9)

M^v(0)=M_R^v为广义标准线性体的松弛弹性模量，为未松弛弹性模量，.基于(9)式可定义广义标准线性体的Q值为(O′Connell and Budiansky, 1978)：

(10)

结合(10)式, 利用最小误差的非线性优化方法拟合广义标准线性体模型Q值与实际介质的Q值可获取τ_εl^v、τ_σl^v取值.M_R取值通过令广义标准线性体与目标介质在某一选定参考频率ω_ref点处的波传播的相速度值(参考波速)相等得到(Blanc et al., 2016).图 3给出了在地震波动数值模拟常用频带[0.1 Hz, 10 Hz]内，基于广义标准线性体模型对两类阻尼介质(介质Q值为常数以及Q值随频率变化)的Q值及相速度拟合结果.这里令阻尼介质在参考频率ω_ref=1 Hz处的相速度为1.由图 3可知，增加标准线性体的个数，在给定频带范围内，可准确拟合介质Q值及波传播的相速度.另外，由于广义标准线性体在选定频率点处的相速度计算较为繁琐(谢志南等，2018)，且广义标准线性体在ω_ref=0 Hz及ω_ref=∞ Hz退化为弹性体，为此常将介质弹性波速作为参考频率ω_ref=0 Hz(Carcione et al., 1988; Martin and Komatitsch, 2009)或ω_ref=∞ Hz处的参考波速(Carcione，2010).值得注意的是，若采用介质弹性波速作为0 Hz或∞Hz处的参考波速，在人为选取模拟频率范围不同的情形，将导致不同频带范围内拟合得到的相速度频散曲线存在较大差异，这一差异可通过给定滞弹性介质在某一频率点处的参考波速加以消除(Emmerich and Korn, 1987；Kristek and Moczo, 2003；袁士川等，2018).但这一差异并不影响PML的精度和稳定性，本文不做深入讨论.

2 弱形式时域滞弹性PML

弱形式时域滞弹性PML可通过对无限域内弱形式滞弹性波动方程进行复坐标延拓得到(谢志南和章旭斌，2017).这一过程可分解为：(1)利用虚功(虚位移)原理建立弱形式无限域波动方程；(2)将这一方程变换至频域，并对其做复坐标延拓得到弱形式频域PML；(3)对弱形式频域PML进行坐标变换将其转换为实坐标系下弱形式频域PML；进而基于傅里叶反变换得到弱形式时域PML.下面以外部无限域A区PML构建为例阐述这一过程.

2.1 利用虚位移原理建立弱形式无限域波动方程

与波动问题的强形式表达(1)—(3)式不同，基于虚位移原理可建立有限区域的弱形式平衡方程即：弱形式波动方程.这一方程耦合了波动方程(1)以及相应的边界条件：(2)、(3)式和无穷远辐射条件.下面阐明图 4中A区弱形式波动方程的建立过程.引入边界条件允许的x_i方向虚位移w_i，将其与波动方程(1)点乘，然后做区域积分，利用高斯定理得到：

(11)

Ω_A表征A区域面积或体积，Γ_MI^Ω_A为介质A、B的分界面，Γ_FB^Ω_A为A区自由地表，Γ_AB^Ω_A为A区人工边界.n_j为外法线矢量分量，n_jdΓ=dΩ/dx_j.二维的dΩ为面积微元，三维的dΩ为体积微元.式(11)隐含假定了A区右侧无穷远处位移、应力为零的无穷远辐射条件.因自由表面处应力为0，即：

(12)

由于Γ_MI^Ω_A、Γ_MI^ΩB外法线单位向量方向相反且介质分界面处应力连续，即：

(13)

由此得到A区弱形式波动方程为

(14)

2.2 将弱形式波动方程变换至频域，并对其做复坐标延拓得到弱形式频域PML

基于傅里叶变换将(14)式变换至频域，即：

(15)

分别为位移和应力张量分量的傅里叶变换.假定坐标原点位于A区左下角(图 4)，引入复坐标延拓，即：

(16)

s_i为解析的复坐标延拓函数.常用的复坐标延拓函数为单极点复频移坐标延拓函数，即：

(17)

其中d_i(x_i)为n_i外法线方向行波衰减因子，κ_i(x_i)为坐标缩放系数，α_i为频移参数.当κ_i(x_i)=1, α_i(x_i)=0时，复频移延拓函数(17)式退化为传统延拓函数，即：

(18)

基于复频移延拓函数构建PML能提高PML对快衰波以及大角度入射平面波的吸收精度，但与基于传统延拓函数构建PML相比，在一定程度上降低了对低频(ω≤α_i(x_i))波动的吸收能力(Festa and Vilotte, 2005).为此，研究人员提出采用多极点复坐标延拓函数构建高阶PML以解决这一问题(Correia and Jin, 2005).公式为

(19)

Feng等(2017)将高阶PML应用于宽频带地震波动数值模拟，结果表明高阶PML精度优于基于单极点复坐标延拓函数构建所得PML.

将(15)式中实坐标替换为复坐标(16)式得到：

(20)

(21)

(22)

为对的偏导算子.二维情形：，三维情形：.复坐标延拓不会改变波动方程解的形式(Teixeira and Chew, 2000).考虑原平衡方程(14)的平面谐波解e^{iωt-ik_xx-ik_yy-ik_zz}，与之对应，仅对x坐标进行复坐标延拓后所得方程的解为.考虑，可得e^{iωt-ik_xx-ik_yy-ik_zz}e^-k_xxd/ω.从这一表达式可知，在实数坐标系下考察复坐标延拓方程的解，其解会沿延拓方向进行空间衰减，且在延拓界面复坐标延拓后方程的解与原方程的解完全一致，为此：

(23)

进而得到A区弱形式频域PML的表达式为

(24)

2.3 对弱形式频域PML进行坐标变换将其转换为实坐标系下弱形式频域PML

由于坐标变换不会改变方程的解，因此在这一过程中我们无须考虑方程中与时间、空间变量之间的关系.实-复坐标变换可通过微分计算链式法则实现，即：

(25)

(26)

(27)

二维情形：，三维情形：.将(25)—(27)式代入(24)式整理得到：

(28)

(29)

利用傅里叶逆反变换将(28)式转换至时域，即：

(30)

式中，J(t)=F^-1[-ω²s_Ω]，，F^-1表示傅里叶反变换.附录A给出了二维情形下(24)式与(30)式的详尽表达式.

3 任意多极点有理分式与频域函数乘积傅里叶反变换的简化计算及滞弹性PML时空离散

由(30)式可知，即便基于单极点复坐标延拓函数构建PML，亦需处理多极点有理分式与频域函数乘积的傅里叶反变换计算问题.下面以式(31)为例，即：

(31)

其中h(t)为任意函数，阐述这一傅里叶反变换普适、简化计算方法.该方法的关键是找出重复极点，通过微调将多重极点转化为单重极点，然后对单重极点有理分式乘积进行依次分解.我们仅调整PML的参数，保持1/τ_σl^v不变.PML参数微调的依据在于k_i(x_i)及α_i(x_i)的变动不会改变PML的完美匹配特性(Collino and Tsogka, 2001; Feng et al., 2017)，且PML精度对较大范围的k_i(x_i)及α_i(x_i)取值并不敏感(Zhang and Shen, 2010; Feng et al., 2017).

考虑到常用延拓函数参数取值方案，即：

(32)

(33)

(34)

其中f₀为模拟频带主频，V^P为PML相邻计算区的P波波速，|x_i|/L为PML层内点距人工边界距离与PML层厚的比值，反射系数R一般取为0.001.依据Zhang和Shen(2010)以及Collino和Tsogka(2001)，a的取值可在1.0~2.0之间，对不同模型κ_max的优化取值区间分布在1.0~13.0之间.但上述区间范围内PML吸收精度变化不大.这里将调整系数取为aπf₀Δx₀/(8L)，其中Δx₀为PML延拓方向相邻节点之间最小间距.这一调整系数为PML中α(x_i)参数最小取值的八分之一.通过微调参数，我们保证每个极点之间的差值大于这一调整系数，且保证调整后的极点仍为大于0的正实数，以κ₁为例，若|α₁-α₂|＜aπf₀Δx₀/(8L)，则将max(α₁, α₂)调整为max(α₁, α₂)+aπf₀Δx₀/(8L).进一步，对非重复极点，基于有理分式展开，即：

(35)

(36)

(37)

(38)

按照上述方法可依次将分解为若干个单极点有理分式的叠加.进一步，由于(a_l+ 的傅里叶反变换为h(t)+(a_l-b_l)[e^-b_ltH(t)*h(t)]，可通过迭代计算得到e^-b_ltH(t)*h(t)离散值, g(t)的计算流程为

(1) 对多极点函数进行逐次降极有理因式分解，即：

(39)

两个单极点有理因式乘积的分解由(35)—(38)式给出.

(2) 将分离得到的单极点函数与h(t)分别作卷积，然后逐个叠加得到g(t)，即：

(40)

在此基础上，采用固定边界截断PML(Duru and Kreiss, 2014)，与建立弹性PML及有限计算区内弹性波动方程的勒让德时空离散一致(Komatitsch and Tromp, 1999; Savage et al., 2010; 谢志南和章旭斌, 2017)，可建立滞弹性PML以有限计算区内的滞弹性波动方程勒让德谱元时空离散方案.详细表达式参看附录1.

4 数值算例

本节基于算例检验滞弹性PML的精度、计算效率及滞弹性近场谱元模拟方案的长持时稳定性.首先，基于全空间均匀介质波动问题论证滞弹性PML对体波的吸收精度.进一步，阐明即便对于介质阻尼较大情形，PML的计算效率优于“完美匹配”物理阻尼吸收层(吸收层内介质与相接计算区阻尼介质完全一致).其次，基于半空间均匀介质波动问题论证滞弹性PML对瑞利面波的吸收精度.最后，基于层状半空间波动问题论证滞弹性PML对复杂波动的吸收精度及新构建滞弹性近场谱元模拟方案的长持时稳定性，同时揭示不加修正地将“弹性介质”PML应用于滞弹性近场波动数值模拟可导致数值失稳.

4.1 全空间均匀介质波动模拟

给定介质参数：密度ρ=2000 kg·m^-3，考虑强衰减阻尼介质Q^κ=12，Q^μ=5，基于一级近似得到的介质弹性P波，S波波速分别为V^P=900 m·s^-1，V^S=500 m·s^-1.

在坐标原点(0 m, 0 m)y方向上施加垂直向下集中力荷载，荷载时间函数为主频f₀=1 Hz的Ricker子波，即：

(41)

其中A为无量纲化幅值系数.通过调整这一系数使得接收点处竖向位移幅值为1.在x、y=±3000 m处设置人工边界截断无限域.

在0.1~10 Hz频带范围基于标准线性体建立阻尼介质时域本构，由于模型尺寸与模拟主波长的比值较小，依据谢志南等(2018)可知，仅需采用包含2个标准线性体的广义标准线性体即可满足模拟精度要求.将弹性波速作为ω=∞ Hz处参考波速，采用最小误差的非线性优化方法可拟合得到广义标准线性体系数.

基于单位插值函数为4阶，单元节点数为25的正方形勒让德谱元对有限计算区进行空间离散.依据谱元模拟精度要求，考虑阻尼介质本构拟合频带最高频率10 Hz，将单元长度设为50 m，大致为一个波长范围内包含一个单元.时间离散采用Newmark-Beta预估校正法(β=0, γ=1/2)，时间离散步长取为8×10^-4s.

在人工边界外设置两类高精度吸收层形式人工边界条件.(1)滞弹性PML，PML层厚为600 m(三层单元)，令a=1, κ_max=1, R=0.001，V^P=900 m·s^-1，复坐标延拓函数中参数取值由(32)—(34)式给出.(2)“完美匹配”物理阻尼吸收层.利用滞弹性介质自身的衰减特性，直接在人工边界外设置同一阻尼介质的吸收层亦能保证计算区域辐射的波动能无反射地穿过人工边界且在吸收层内逐渐衰减.阻尼吸收层厚取为4500 m、9000 m以及13500 m，分别对应为5λ_P、10λ_P、15λ_P，其中λ_P为主频f₀=1 Hz处P波波长，λ_P=V^P/f₀=900 m.

图 5给出了采用滞弹性PML与不同层厚“完美匹配”物理阻尼吸收层情形下，(2800 m, 2800 m)处观测点的水平和竖向位移模拟波形及对应的大区域解.由图 5可知，3层单元滞弹性PML能有效吸收从计算区域向外辐射的体波，模拟波形与大区域解基本完全吻合，且优于“完美匹配”物理阻尼吸收层.考虑到标准线性体个数为2情形，单个滞弹性PML单元的计算量为单个“完美匹配”物理阻尼吸收层内单元计算量的6倍，本算例中PML层内单元个数为396.5λ_P、10λ_P、15λ_P厚，阻尼吸收层内单元个数分别为4725、13500、26325，由此可知即便对强衰减阻尼介质，滞弹性PML的计算效率亦远优于“完美匹配”物理阻尼吸收层.对弱衰减介质，由于介质对波的衰减速度变慢，滞弹性PML的计算效率更为突出.图 6给出了采用滞弹性PML情形，计算区内总能量E随时间的衰减曲线，公式为

(42)

其中为全局速度矢量，第一项为动能，第二项为势能.从图 6可知，新构建滞弹性近场谱元模拟方案是长持时稳定的.

4.2 半空间均匀介质波动模拟

考虑半空间均匀介质波动模拟，介质参数同上一算例.自由面设置在y=0 m处.在(500 m, 0 m)y方向上施加垂直向下方向集中力荷载，荷载时间函数为主频f₀=1 Hz的Ricker子波.在x=0 m、5000 m以及y=-2500 m处设置人工边界截断无限域.

在0.1~10 Hz频带范围基于标准线性体建立阻尼介质时域本构，由于模型尺寸与模拟主波长的比值较小，采用包含3个标准线性体的广义标准线性体即可满足精度要求.将弹性波速作为ω=∞ Hz处参考波速，采用最小误差的非线性优化方法可拟合得到广义标准线性体系数.

基于单位插值函数为4阶，单元节点数为25的正方形勒让德谱元对有限计算区进行空间离散.单元长度设为50 m.时间离散采用Newmark-Beta预估校正法(β=0, γ=1/2)，时间离散步长取为8×10^-4s.

在人工边界外层厚为150 m的滞弹性PML，复坐标延拓函数内参数取值同算例1.图 7给出不同时刻位移波场快照图.图 8给出(3500 m, 0 m)处接收点水平及竖向位移时程，且为对比分析，图 8亦给出了人工边界条件设为黏性边界情形下接收点处位移时程(Komatitsch and Tromp, 1999).由图 7、8可知，滞弹性PML对强衰减的瑞利面波同样具有极高的吸收精度.同时图 8表明，在滞弹性介质中，瑞利面波幅值并不一定总是远大于P波振幅，在本算例中瑞利面波对应的位移幅值与P波位移振幅相接近.

4.3 层状半空间波动模拟

给定两层介质参数，上覆软基岩介质参数：密度ρ₁=2000 kg·m^-3，P波波速V₁^P=1800 m·s^-1，S波波速V₁^S=1040 m·s^-1，Q₁^κ=120，Q₁^μ=50，覆盖层厚度为500 m.下卧基岩半空间介质参数为：ρ₂=2550 kg·m^-3，P波波速V₂^P=4500 m·s^-1，S波波速V₂^S=2630 m·s^-1，介质品质因子Q₂^κ=680，Q₂^μ=263.自由表面位于y=0 m处.在(5000 m, -6000 m)处施加位错震源，矩震级M_W=3.3，地震矩张量为.震源时间函数为主频f₀=1 Hz的Ricker子波.

在x=0 m、20000 m、y=-10500 m处设置人工边界截断无限域.采用包含2个标准线性体的广义标准线性体在[0.1 Hz~10 Hz]频带范围近似滞弹性介质本构.在采用最小误差的非线性优化方法拟合广义标准线性体系数过程中，将弹性波速作为ω=∞ Hz处参考波速.基于单位插值函数为4阶，单元节点数为25的正方形勒让德谱元对有限计算区进行空间离散.单元长度设为100 m，时间离散采用Newmark-Beta预估校正法(β=0, γ=1/2)，时间离散步长取为8×10^-4s.

在人工边界处设置三类人工边界条件：(1)弹性PML，层厚为300 m(三层单元)，令a=1, κ_max=1, R=0.001复坐标延拓函数中参数值由(32)至(34)给出，V^P=max(V₁^P, V₂^P)=4500 m·s^-1；(2)滞弹性PML，层厚为300 m(三层单元)，其中延拓函数参数取值同与弹性PML；(3)黏性边界.总计算时长为80 s，即迭代步数100000步.图 9中分别采用不同人工边界条件计算区域内的能量衰减曲线.由图 9可知，将弹性PML直接应用于滞弹性波动数值模拟将导致数值失稳.基于滞弹性PML和黏性边界建立的滞弹性近场波动数值模拟方案长持时稳定.进一步，图 10给出接收点处(3700 m, 0 m)和(15700 m, 0 m)处位移时程，如图 10所示，滞弹性PML能有效吸收由于介质分界面存在而导致的复杂外行波动，吸收精度远优于黏性边界.计算所得接收点位移时程与大区域解几乎完全一致.本文算例所采用的均为四边形单元，但文中结果可直接推广至Meng和Fu(2017)建立的三角形有限元单元波动模拟优化方案中.

5 结论

本文基于复坐标延拓技术变换弱形式无限域滞弹性波动方程建立了弱形式时域滞弹性PML.给出了明确的PML参数微调准则，建立了多极点有理分式与频域函数乘积傅里叶反变换的普适、简便计算方法，极大地简化了PML计算.基于该方法可实现任意高阶PML.在此基础结合谱元建立了滞弹性近场波动谱元模拟方案.基于算例检验了PML的精度、计算效率.结果表明：PML对体波，面波均有极好的吸收精度.新构建滞弹性近场谱元模拟方案长持时稳定.同时阐明了不加修正地将“弹性介质”PML应用于滞弹性近场波动数值模拟可能导致数值失稳这一问题.新方案可直接推广应用至基于任意高阶滞弹性PML的宽频带地震波动正反问题研究(Komatitsch et al., 2016; Feng et al., 2017).

附录A 二维弱形式频域PML的详尽表达式

基于分量形式表述弱形式频域滞弹性PML(28)式得到：

(A1)

(A2)

(A3)

(A4)

(A5)

利用多极点的分离方法，并作傅里叶逆变换可得：

(A6)

(A7)

(A8)

(A9)

其中：

(A10)

(A11)

式中各系数为

(A12)

(A13)

(A14)

(A15)

(A16)

(A17)

(A18)

(A19)

等式(A1)左侧第一项的时域表达式涉及卷积积分项，可通过两极点的分离得到，即：

(A20)

式中各系数为

(A21)

(A22)

(A23)

(A24)

(A25)

将(A20)式及(A6)—(A7)式代入(30)式整理得到：

(A26)

采用勒让德谱元离散(A26)可得局部节点系运动方程为

(A27)

其中：

(A28)

(A29)

(A30)

(A31)

(A32)

(A33)

(A34)

(A35)

(A36)

(A37)

u_q^e为单元q节点上的位移矢量，N_q表示单元q节点上插值函数，N_{q, x}、N_{q, y}单元q节点上插值函数对x、y的偏导数.经组装可得到近场波动勒让德谱元模拟所需全部节点上的运动方程组为

(A38)

全局矢量、P、R分别由、P_p^e、R_p^e组装得到.M为全局对角质量阵.考虑计算计算区内部荷载，须在式(A37)右端补充相应的全局荷载向量.

致谢  感谢审稿专家在论文评审过程中提出的宝贵意见和建议，感谢廖振鹏研究员一直以来对作者工作的支持，感谢马新生提供的波场快照图绘制程序初稿，感谢陈少林教授在本文撰写过程中提出的修改意见.