多轴复频移近似完全匹配层弹性波模拟

罗玉钦 刘 财

(吉林大学地球探测科学与技术学院，吉林长春 130026)

0 引言

为满足地震波正、反演的需要，边界条件的研究引起人们的极大关注。目前应用最广泛的边界条件是由Berenger[1-2]在研究麦克斯韦方程时提出的完全匹配层(Perfectly matched layers，PML)。Chew等[3]和Collino等[4]将其解释为复坐标伸展变换的结果。当入射波角度很大时，衰减系数会变得很小，因此常规PML无法吸收大角度入射波，且在极低频时产生奇异值。复频移(Complex frequency shifted，CFS)技术对复坐标伸展变换(CCS)做进一步处理，原来为1的项被尺度因子替换并在大角度入射波的吸收中起作用，在分母中加上频移因子以消除低频奇异值[5]。Roden等[6]提出基于CFS卷积的PML(即C-PML)； Komatitsch等[7]进一步将其推广到弹性介质中。

此外，近似PML(Nearly PML，NPML)[8]不同于上述常规PML，它是直接对波场进行变换而得到，不改变波场形式，也不需做卷积处理。该方法[8]自提出后很快被推广试用于声波介质、弹性介质以及双相介质中[9-11]。但是，由于该方法在大入射角时吸收能力较弱且在复杂介质的模拟中存在不稳定现象，因此，在地震正、反演中现今主流方法仍旧是C-PML。

对于边界条件的稳定性，Festa等[12]指出PML边界与面波之间会降低吸收效果并产生不稳定性。Mezafajardo等[13]认为常规PML在各向同性和各向异性中都不满足严格的渐进稳定性，并引入多轴技术，即在多个正交方向上同时引入衰减因子。同时，衰减因子对边界吸收效果有很大影响。Collino等[14]提出内截断边界与PML层的距离呈指数关系的衰减函数； 此后Groby等[15]对其进行了改进和整理。近年来，中国学者也对PML做了许多研究[16-28]。田坤等[29]将多轴技术应用于C-PML中； 陈可洋[30]提出了正弦型和余弦型衰减函数； 罗玉钦等[31]对衰减因子进行了系统研究，提出的衰减因子将吸收效果提升了20%～60%。

虽然NPML拥有其自身优势，但它仍存在稳定性与吸收效果上的问题。本文基于NPML并引入频移因子、尺度因子及稳定性因子，提出了多轴复频移NPML(MCFS-NPML)吸收边界，旨在提升NPML对大角度入射波的吸收且增强其稳定性； 在此基础上采用文献[31]提出的衰减函数，进一步提升MCFS-NPML的吸收能力。本文详细分析MCFS-NPML的稳定性和吸收效果，探讨各因子在改善吸收效果及稳定性方面所起的作用，证实了MCFS-NPML拥有更好的吸收效果和稳定性； 同时，将该衰减函数与广泛运用的指数型衰减函数进行比对，结果表明在该衰减函数下吸收效果得到进一步改善。

1 理论研究

1.1 近似完全匹配层(NPML)吸收边界

本文采用交错网格有限差分法，因此采用一阶速度—应力形式方程组。将方程组变换到频率域

(1)

式中：ω是角频率；λ、μ是拉梅常数；ρ是密度；Vx、Vz为频率域弹性波场速度分量；Txx、Tzz、Txz为频率域弹性波场应力分量。

采用复坐标伸展变换的形式为

(2)

(3)

将式(3)代入式(1)，并直接作用在速度分量和应力分量上，得到复坐标伸展变换后的方程组

(4)

(5)

(6)

同理，采用相同方法也将其他变量变换到时间域，有

(7)

且有

(8)

ξ=τxx,τzx,τzz,vx,vzm=x,z

式中：τxx、τzx、τzz是时间域应力；vx、vz为时间域速度参数。分别与频率域中的Txx、Tzx、Tzz和Vx、Vz相对应。

从式(7)可见主控方程的形式并未发生变化。这里需引入8个辅助变量替换原主控方程中的变量。在程序编写时不需改动主代码，易于实现。式(8)中的原始控制方程中的辅助变量是原变量经伸展变换得到，是将原波场直接变换后再代入相同形式的控制方程求解，因此易于将NPML推广到其他介质。在计算区域中衰减因子值为0，故式(8)中原变量与新变量是相同的，因此在编程时可只在PML区域中进行波场变换以减少存储空间。

1.2 复频移近似完全匹配层(CFS-NPML)吸收边界

先分析NPML的平面波解，经过复坐标变换之后方程具有新的平面波解

(9)

修改式(3)，并引入尺度因子βx(x)和频移因子ηx(x)，可得

(10)

得到新的平面波解的衰减项为

(11)

尺度因子βx(x)使大角度入射波的传播方向向法线方向弯曲，导致衰减系数增加； 频移因子ηx(x)对入射角的影响不大，但可有效地避免奇异值的出现。当入射角较大时，在常规PML中波只能在很浅的介质中传播，无法达到较好的吸收效果，尺度因子的引入增强了边界对这一类波的吸收[7]。衰减函数、尺度因子与频移因子的表达式为

(12)

(13)

式中：L是PML的层厚；R是理论边界反射系数；l为PML区计算点与内部边界的距离；η0的经验值范围是0～6；β0的经验值范围是1～10； 指数qβ、qη与n相同，通常取1～3。

以式(10)替换式(3)，如Txx

(14)

(15)

=[ηx(x)+iω]Txx

(16)

=ηx(x)Txx+iωTxx

(17)

转换到时间域，有

(18)

同理，求得其余变量的微分方程且替换式(8)，则可简写为

(19)

ξ=τxx,τzx,τzz,vx,vzm=x,z

采用交错网格有限差分数值模拟时，将式(19)离散，得到如下离散格式

(20)

1.3 多轴复频移NPML(MCFS-NPML)吸收边界

(21)

式中αx的取值不再是0，而与原衰减函数成比例，比例系数P(x/z)即是稳定性因子 。

图1 衰减因子PML区域示意图

将双衰减剖面应用于NPML，如区域2，有

(22)

对于 CFS-NPML，在区域1和区域2中的处理方法与上面不一样。同样以区域1为例，变换函数为

(23)

将式(23)代入式(14)中，得到

(24)

这样就避免了单剖面衰减引起的平面波解呈指数增长。采用多轴衰减剖面引入了稳定性因子，使系数矩阵特征值向负半轴移动，消除了不稳定性。稳定性因子应小于1，且随介质的复杂程度而改变，不宜太大，否则会导致虚假反射增强。

2 数值模拟

在地震波模拟中分别实现NPML、M-NPML、CFS-NPML及MCFS-NPML。选取的时间步长为1ms，加载震源是主频为25Hz的雷克子波。构建模型的尺寸是5200m×1200m(包含边界厚度)，x方向和z方向网格步长都是10m。纵波速度为3000m/s，横波速度为1400m/s，密度为2000kg/m3，PML层数L设为10，R取10～6。将震源加载于x=2610m、z=110m处，贴近于PML边界，往两侧传播的波近似平行于上边界传播，截取上述几种近似完全匹配层下的地震波波场快照。

图2a和图2b分别显示上述四种PML在0.90s和0.95s时刻地震波传播情况。从上到下依次呈现稳定性因子、尺度因子、频移因子对匹配层的影响及三个因子共同作用下的波场； 图2a第一排标示的①、②和③分别对应直达波及在下、上边界形成的反射波； 图2b第一、第四排的第一列中黑色椭圆即匹配层中残余能量团，该能量团是因波在边界处衰减不完全产生的损耗波，且入射角度越大，能量团越明显。该能量团的存在会导致虚假反射回传到计算区域，且能量团的积累使计算不稳定。

通过图2a第一排可见增加稳定性因子时反射波能量被增强； 从图2b第一排中可看出稳定性因子值越大，边界区域越干净，消除了大部分的残余能量团。对于CFS-NPML，选取多组参数分析频移因子和尺度因子对近似完全匹配层的影响。通过图2a对比NPML与η0=0、β0=2 时的CFS-NPML，可见残余能量团稍有减弱，且下边界处的虚假反射也得到削减，但幅度较小。

对比η0、β0分别取1、2和5这几组波场快照可知：改变η0值会显著影响吸收效果，且η0取值越大吸收效果越好，但在边界中的波“拖尾”现象更严重(图2b第三排)； 但β0取值过大，反而会影响吸收效果。对于MCFS-NPML、CFS-NPML中边界处未能完全吸收的残余能量被再次削弱，MCFS-NPML能进一步清除PML层中的能量以保持匹配层的稳定且再次压制了因大角度入射所带来的虚假反射(图2b第四排)。

图3中0时刻介质中的能量迅速增加，在约0.30s时因波的下侧已到达边界，波的能量逐渐减弱。图中实线和虚线分别对应内区间(不包括PML边界)和整体区域的能量曲线。在0.83s时刻波抵达左右两侧，经过PML层作用能量迅速降低，在约2.50s时残余的反射波再次抵达左、右两边界，能量再次减少。图3a中NPML下虚线与实线的差异较大，表明有大量能量残留在边界中。该残留能量会使PML层变得不稳定，同时过多的能量积累会影响边界的吸收性能。M-NPML曲线差随稳定性因子取值增大而逐渐减小。当稳定性因子取0.1时，虚线与实线几乎重合，证明边界中的大部分能量均被消除。对于CFS-NPML，从β0取2、η0取0、1、2、3、5这几组能量曲线(图3c)看出，η0值越大，在约0.83s时能量曲线最低。

正演模拟中，期望地震波在第一次达到边界处时就被尽可能地吸收，防止地震波反射回计算区域。频移因子改善了PML吸收性能，但随着η0增大，左右两侧初次未被吸收的微弱反射波再次到达边界时，整体区间的能量曲线扰动剧烈，是不稳定的。β0取值过大也会有这样的扰动，且CFS-NPML中曲线差依旧很大； 而在引入稳定性因子之后(MCFS-NPML)，能量扰动和曲线差过大将被压制。

图2 三种因子不同组合下的波场快照

图4是模型中x=110m、z=110m点处0～3.20s时段的x方向速度分量的变化曲线，可清晰地看到上述三个因子对常规入射波吸收效果的影响。稳定性因子取0.1时，来自底界的虚假反射波最强。从反射波振幅看，CFS-NPML比NPML和M-NPML低很多，来自另一侧的反射波随η0增大而减小； 对于MCFS-NPML，不仅反射波能量微弱，来自另一侧的反射波衰减也特别明显。证实了M-NPML会增强虚假反射，CFS-NPML和MCFS-NPML却有利于入射波的吸收。

为进一步测试边界的稳定性与吸收性能，采用了复杂Marmousi模型。在该模型中不稳定性被放大，震源靠近边界，大角度的入射波进入上边界并截取0.5s和4.0s时的波场快照(图5)。

在NPML中，上边界存在强烈的地震波残留(图5a，红色椭圆所示)，CFS-NPML也是如此。从4.0s波场快照可见NPML下整个计算区域已因误差积累而发散。至于M-NPML(图5b、图5c)，随着稳定性因子的增加，波场残留减弱。β0取2、η0取0时的CFS-NPML(图5d)与NPML对比，虽然产生累积误差，但污染区域相对较小，表明β0有微弱压制作用； 当η0取2时(图5e)也未见波场被污染，表明频移因子对稳定性有一定作用。而MCFS-NPML(图5f)因为稳定性因子的引入，上边界比CFS-NPML更干净。

图3 能量衰减曲线

图4 反射波振动曲线

图5 0.5s (左)和4.0s(右)时刻Marmousi模型下的波场快照

简而言之，在提高稳定性方面，三个因子都有贡献，但机理不同。稳定性因子沿平行匹配层方向增加一个衰减剖面，吸收沿该方向传播的耗散波； 频移因子避免了低频奇异值的产生； 尺度因子使波向法线方向弯曲，增强对匹配层中波的吸收能力。Marmousi模型模拟中的不稳定随频移因子的引入而消除，故该不稳定是低频奇异值所致。虽然CFS-NPML下计算保持稳定，但上界残余能量团十分明显，因此需同时引入三种因子使NPML具有更好的稳定性和吸收效能。

3 衰减函数

尝试通过修改衰减函数，进一步提升边界吸收能力。用网格算法模拟地震波传播，衰减函数不再是一条连续的曲线，匹配层之间的衰减函数值差异较大，引起虚假反射。通过增加层数缩小该差异，但会增加存储量，降低计算效率。探讨在不增加层数的基础上进一步削弱离散差异带来的虚假反射。设计下式

(25)

(26)

波场模拟采用MCFS-NPML，且β0=2、η0=2及P=0.005，模型尺寸为2000m×2000m。当层数取8时，γ=0.065、δ=0.094这一组系数下波反射最弱。图6b中蓝线和黑线对应传统指数型衰减函数和式(26)给出的函数下的反射情况，可看出新函数下的虚假反射更弱。新函数取8层达到的吸收效果与原函数取10层相当。

图6 吸收效果图(a)与波形曲线(b)

4 结束语

为将NPML更好地应用于正反演，本文推导出MCFS-NPML吸收边界条件。该过程中引入了频移因子、尺度因子和稳定性因子，这三种因子保证了波场模拟稳定有效地进行。将它们同时引入弥补了仅引入某种单一因子的不足，相互弥补了各自固有的不利影响。通过分析衰减函数的引入过程以及分布方式，改进了衰减函数，最终达到进一步提升吸收效果的目的。