0 引言

地震波数值模拟是了解复杂地质条件下地震波传播规律的有效工具，可以用来指导地震观测系统设计、数据解释和检验地震数据处理方法。然而，受计算内存和计算时间的限制，在对波动方程进行数值求解时，我们总是在一个有限的空间区域进行计算，这样自然就引入了人为边界，从而产生人为边界反射。这些人为边界反射严重干扰了由反射界面和自由界面反射的有效波。为此，已经发展了多种消除人为边界反射的方法。Clayton等^[1]在傍轴近似理论基础上导出了二维声波方程吸收边界条件，该吸收边界条件对特定的入射角和频率范围内的边界反射具有较好的吸收效果。Cerjan等^[2]通过在计算边界附近引入损耗介质来衰减向外传播的波；由于在两个具有不同吸收系数损耗介质的界面上同样会产生反射，因此需要有较多的吸收层才能达到很好吸收效果。Berenger^[3]针对电磁场有限域波动问题给出了一种高效的完全匹配层(PML)吸收边界条件，并证明该方法可以完全吸收来自各个方向、各种频率的电磁波，而不产生任何反射。目前人们已将PML思想引入到地震波数值模拟中，并发展了多种规则网格下的PML吸收边界条件^[4-6]。

然而在实际地震勘探中，起伏地表和复杂构造是普遍存在的。为此，地球物理学家们提出了多种模拟起伏地表下地震波场的数值方法。其中有限差分法因其计算效率高、对内存需求较小而得到广泛应用。Jih等^[7]为了研究高山和峡谷地区的地震波传播特点，对Ilan^[8]提出的一种处理二维折线状起伏地表的有限差分法进行了改进；但是，该方法需要做复杂的坐标变换，且很难推广到三维情况下。Oprsal等^[9]将真空法自由边界条件引入到阶梯状近似的起伏地表上；该方法的缺点是精度低，并且容易出现不稳定现象。Frankel等^[10]提出了一种密度渐变方法，该方法实际上是对真空法的一种改进，同样使用阶梯网格对起伏地表进行近似；但是这种方法精度较低，而且会在阶梯状网格角点处引起虚假散射波。为了彻底消除阶梯状网格引起的虚假散射波，人们发展了曲线网格法^[11-14]。Appelo等^[15]提出了一种非常有效的曲线坐标系离散自由边界条件，并证明了该方法的有效性和稳定性。后来，Lan等^[16]把该方法扩展到了三维各向异性介质的地震波数值模拟中。相比于矩形网格有限差分法，曲线网格法不仅可以有效消除阶梯状网格引起的虚假散射波，而且可以减少离散化时的网格点数，从而节省内存、提高计算效率。然而，到目前为止，有关曲线坐标系下PML吸收边界条件的文章还比较少。Festa等^[17]首次将曲线坐标系下的PML吸收边界条件应用到谱元法地震波数值模拟中。Gao等^[18]利用三角网格计算地震波场时提出了一种非规则网格PML吸收边界条件。

为了有效消除起伏地表下地震波数值模拟引起的边界反射，本文提出了一种适用于曲线坐标系下二阶弹性波方程的PML吸收边界条件，并将该边界条件应用到正交曲线坐标系下的地震波数值模拟中。文中采用孙建国等^[19]提出的正交贴体网格对起伏地表模型进行网格剖分。由于边界条件在正交曲线坐标系下无需进行坐标旋转和插值运算，从而简化了计算过程，提高了计算精度。

1 方法原理

笛卡尔坐标系下二维非均匀各向同性弹性波方程为

式中：u和v分别表示x、z方向上的位移；ρ为密度；λ和μ为lame系数。

对于起伏地表下的地震波数值模拟，为了压制阶梯状网格引起的人为虚假散射波，网格线与起伏地表需要保持一致。这种网格被称为“贴体网格^[20]”，已经有许多研究者在地震波数值模拟中采用了这种网格^{[15, 21-22]}。贴体网格可以通过由计算空间到物理空间的坐标变换来获得(图 1)。通过这种变换，将曲线坐标系(ξ, η)映射成物理域上的笛卡尔直角坐标系。本文采用孙建国等^[18]提出的方法生成正交贴体网格，具体的网格生成过程不做赘述。

通过生成贴体网格，笛卡尔坐标系下的每个网格点都由曲线坐标系唯一确定。因此，我们可以通过链式法则，由笛卡尔坐标系(x, z)下的波动方程得到曲线坐标系(ξ, η)下的波动方程。

根据链式法则，有

式中：ξ_x、ξ_z、η_x、η_z分别表示∂ξ(x, z)/∂x、∂ξ(x, z)/∂z、∂η(x, z)/∂x、∂η(x, z)/∂z；∂_x、∂_ξ、∂_z、∂_η分别表示变量对x、ξ、z、η的偏导数。利用式(1)和式(2)可以得到曲线坐标系下的弹性波方程：

为了得到地震波在PML里的波动方程，把波动方程(3)、(4)变换到频率域，并对频率域里的波动方程做如下坐标拉伸变换^[23]：

式中：为复数坐标扩展函数；s为积分变量；d_s(x)为衰减系数；，ω为圆频率。由式(5)可以得到复数坐标系中空间导数的表达式：

式中，d_s≥0，是使地震波振幅在PML里指数衰减的系数。把式(6)带入到由式(3)和式(4)得到的复杂拉伸坐标系下的波动方程，得：

式中：U和V分别表示u和v对t的傅里叶变换；d_ξ=－(n+1)lg(R)(v_P/(2δ))(ξ/δ)^2[23]，是ξ方向上的衰减系数；d_η=－(n+1)lg(R)(v_P/(2δ))·(η/δ)²是η方向上的衰减系数；v_P为纵波速度，δ为PML的厚度，n为常数，用来控制衰减系数的变化率，通常取2~3，R是理论反射系数，通常取10^-2~10^-6。

利用波场分解技术^[24]，把位移水平分量分解成3个部分：U=U₁+U₂+U₃，带入到方程(7)的左边，并把带入后的方程分解成如下3个方程：

对方程(9)、(10)、(11)直接做反傅里叶变换，可以得到时间域PML方程，但是这样会引入卷积运算。为此我们用(iω+d_ξ)²、(iω+d_ξ)(iω+d_η)和(iω+d_η)²分别乘上方程(9)、(10)和(11)的两边，并引入两个中间变量R₁和R₂(见附录)，再作反傅里叶变换，可以得到时间域PML吸收边界条件：

式中：r₁和r₂分别表示R₁和R₂的反傅里叶变换；u=u₁+u₂+u₃，u₁、u₂和u₃分别表示U₁、U₂和U₃的反傅里叶变换；∂_t表示变量对时间t求导；d_ξ′=∂d_ξ/∂ξ；d_η′=∂d_η/∂η。

同理，把位移水平分量分解成3个部分：V=V₁+V₂+V₃, 并引入两个中间变量R₃和R₄(见附录)，可根据式(8)得到如下时间域PML方程：

式中：r₃和r₄分别表示R₃和R₄的反傅里叶变换；v=v₁+v₂+v₃，v₁、v₂和v₃分别表示V₁、V₂和V₃的反傅里叶变换。

由方程(12)和(13)可以看出，通过在频率域中引入中间变量R₁、R₂、R₃和R₄，可以避免将频率域里的PML方程转化到时间域时引起的卷积运算。

在计算域中利用PML吸收边界的基本方法如图 2所示。在常规区域，令d_ξ=d_η=0，方程(12)和(13)退化成标准的波动方程(3)、(4)。在区域1中，令d_ξ≠0，d_η=0；在区域2中，令d_ξ=0，d_η≠0；在区域3中，令d_ξ≠0, d_η≠0。这样在计算区域的左、右和底边界都有PML吸收介质，波由区域内通过边界传播到完全匹配层时不会产生任何反射。

在顶边界我们采用自由界面条件。由于采用了正交贴体网格，因此可以直接实施笛卡尔坐标系下的自由边界条件，无需做复杂的坐标转换^[25]。

2 有限差分格式

由于波在PML区域内传播是逐步衰减的，故可以采用低精度的差分进行离散，因此方程(12)的差分格式可以写成：

式中：(U₁)_{i, j}^k、(U₂)_{i, j}^k、(U₃)_{i, j}^k、U_{i, j}^k、V_{i, j}^k、(R₁)_{i, j}^k、(R₂)_{i, j}^k分别表示u₁、u₂、u₃、u、v、r₁、r₂在k时刻的离散值；h为空间网格步长；Δt为时间采样间隔。计算上述方程时，先根据U^k和V^k计算R₁^k和R₂^k，再计算U₁^k+1、U₂^k+1、U₃^k+1，最后根据U^k+1=U₁^k+1+U₂^k+1+U₃^k+1计算U^k+1。同理可以求取方程(13)，具体过程不再赘述。

3 模型试算

为了验证本文方法的有效性，我们设计了一个正弦函数起伏地表模型(图 3a)。该模型的地表可用正弦函数表示为，即模型地表的起伏周期为1 200 m，相对最大高差为400 m；模型大小为5 000 m×5 000 m。模型介质纵、横波速度及密度分别为v_P=3 000 m/s、v_S=2 500 m/s、ρ=2 000 kg/m³。震源采用主频为20 Hz的Ricker子波，位于(x, z)=(2 500 m, 50 m)处。利用正交贴体网格对模型进行网格剖分，网格点数为501×501，平均网格间距为10 m，吸收边界宽度为20个网格点(图 3b)。由图 3b可见，本文算法生成的网格具有很好的正交性和光滑性。图 4和图 5分别为t=1.6 s的水平分量和垂直分量波场快照。从图 4、图 5中可以看出，本文提出的PML吸收边界在左、右边界以及左上角、右上角的吸收效果均非常理想，即使入射角较大，边界反射也吸收得非常干净。图 6和图 7分别为t=2.4 s的水平分量和垂直分量波场快照。从图 6、7中可以看出，本文方法在左下角、右下角和下边界同样能达到理想的吸收效果。

图 8和图 9是模拟地震记录的水平分量和垂直分量。由图 8、9可见：不加吸收边界时，波在边界上产生明显的边界反射；而采用PML吸收边界时，波在边界上几乎没有产生边界反射，而且对反射纵波和横波都有较好的吸收效果。

4 结束语

本文提出了一种曲线坐标系下的PML吸收边界条件，并结合有限差分法模拟了起伏地表条件下的地震波场。该方法通过在频率域中引入4个中间变量，避免了复杂的卷积运算。数值算例表明，本文提出的PML吸收边界对入射到人工边界的各种类型波都具有很好的吸收效果，保证了有效波场不受虚假反射波的干扰。

附录

为了避免将频率域PML吸收方程变换到时间域时引起的卷积运算，引入4个中间变量：