1 引 言

地震波在地下介质中传播时，不仅受介质的各向异性影响，还受地层介质内在的黏滞性影响^[1].因此，为了准确地描述地震波在地下介质中的传播特征和更精确地指导地震数据采集､处理和解释，应选取能够同时模拟地层各向异性特征和黏滞性特征的黏弹各向异性介质模型来进行波场数值模拟与分析.

目前，黏弹各向异性介质的地震波场数值模拟研究已经取得了一些进展.Carcione^[1-3]基于广义线性固体模型引入新的黏弹性各向异性本构关系，建立了黏弹各向异性介质中随时间而变化的松弛矩阵，而且推导出了黏弹VTI介质的波动方程，并模拟了纵波､横波的各向异性和衰减特征.张忠杰等^[4-5]用近似解析分析法模拟了黏弹EDA 介质中的波场.杜启振等^[6-8]在Carcione的基础上给出了以各向异性主轴方位为参数的各向异性黏弹性介质波动方程，并先后采用有限元法､伪谱法和有限差分法对方位各向异性黏弹性介质中的波场进行了模拟与分析.郭智奇等^[9]利用伪谱法模拟了黏弹VTI介质中的波场.谭未一等^[10]采用规则网格有限差分法模拟了黏弹VTI介质中的波场.李军等^[11]利用交错网格高阶有限差分对黏弹性VTI介质中的波场进行了数值模拟.然而，这些研究对黏弹各向异性介质的波场模拟主要针对两种特殊情况:黏弹VTI介质或黏弹HTI介质.而具有任意倾斜对称轴横向各向同性(TTI)介质在实际地层各向异性的表现中更具有广泛性^[12-14].因此，如果能给出黏弹TTI介质波动方程，并针对其进行波场数值模拟与分析，则更具有一般性.

黏弹TTI介质可以看成是黏弹VTI介质经过旋转一定的角度后得到的，其对称轴与在地面建立的观测坐标系存有一定的夹角，如果利用普通交错网格法求解二维三分量的黏弹TTI介质波动方程，为计算应力而求质点振动速度的空间偏导数时，有些偏导不能直接求解，需要进行漂移或借助辅助网格.而且在普通交错网格中，要求弹性模量定义在所有的半网格点和整网格点上，而实际中通常只定义在半网格点上，这就需要对弹性模量进行插值^[15].因此，当弹性模量变化很大时，普通交错网格法就容易出现计算不稳定的问题^[16].Saengere等^[17-18]首先提出了旋转交错网格有限差分算法，并指出其可以克服普通交错网格差分的上述缺陷.而后，一些学者利用旋转交错网格有限差分法分别对孔隙介质､各向异性介质等的波场进行了数值模拟^[19-21].

旋转交错网格有限差分的所有物理量只定义在两个不同的位置.其中，应力(或质点振动速度)定义在离散网格单元的中心，质点振动速度(或应力)定义在离散网格单元的顶点，故可将所有的弹性模量定义在相同的位置，且可与应力所在的位置相同^{[17， 19]}.因此，采用旋转交错网格有限差分法，在模拟非均匀介质时，弹性模量不需要进行平均;在模拟TTI介质时，计算应力的所有物理量的空间偏导数可以直接求解，且弹性模量也不再需要进行插值处理^[19].所以旋转交错网格差分法更适合黏弹TTI介质的波场模拟，并且相对于普通交错网格差分法，旋转交错网格差分法可以得到更加稳定､更可信的解.

本文在Carcione^[1-2]提出的黏弹各向异性本构关系基础上，给出黏弹TTI介质的二维三分量一阶速度-应力方程，然后采用旋转交错网格任意偶数阶精度有限差分格式求解该方程，并将PML 引入到黏弹TTI介质波动方程的旋转交错网格差分数值模拟中，实现了该类介质波场的旋转交错网格高阶有限差分数值模拟，并且根据模拟结果具体分析了地震波在传播过程中所表现出的各向异性和黏弹性特征.

2 黏弹TTI介质波动方程 2.1 观测坐标系内的本构关系

Carcione的黏弹各向异性理论认为介质对地震波的吸收衰减作用是由大量耗散机制引起的，且各种耗散机制可以用黏弹性的本构关系进行描述，采用其建立的应力应变关系^[1]，即

(1)

(2)

其中，x_v(0)(v= 1，2)表示0 时刻的弛豫模量，当v=1时，对应于纵波，当v=2时，对应于横波;T^T =(T₁，T₂，T₃，T₄，T₅，T₆)= (σ_xx，σ_yy，σ_zz，σ_yz，σ_zx，σ_xy)表示应力向量;S^T = (S₁，S₂，S₃，S₄，S₅，S₆)= (ε_xx，ε_yy，ε_zz，ε_yz，ε_zx，ε_xy)表示应变向量;A_LK和A^(v)_LK(L，K=1，2，…，6;v=1，2)分别表示非弛豫空间函数和弛豫空间函数，并且其与弹性参数有关;L_v(v=1，2)表示纵波或横波弛豫机制的总数;e(v)Jl为记忆变量，e^(v)_Jl·是记忆变量的时间导数;φ_vl(0)(l=1，2，…，6;v= 1，2)表示第l个弛豫机制的响应函数;τ^(v)_σl(l=1，2，…，6;v=1，2)表示第l个机制的应力弛豫时间.上述变量及函数的具体表达式可以参见参考文献^[1].

根据黏弹各向异性介质的本构关系，坐标旋转以后，可以得到观测坐标系下以任意方位角和倾角为参数的应力应变关系方程^[1]，即

(3)

式中，M为旋转变换矩阵^[12]，且

(4)

其中，θ和φ 分别为在观测坐标系下的方位角和倾角.

将式(3)变形可得

(5)

(6)

2.2 一阶速度-应力方程

由黏弹性各向异性介质的本构关系，给出黏弹TTI介质的二维三分量一阶速度-应力方程^{[1， 2]}.

运动方程:

(7a)

(7b)

(7c)

其中，ρ为介质的密度，v_x､v_y､v_z分别是质点速度x､y､z的分量，σ_xx､σ_zz为正应力，σ_xz､σ_xy､σ_yz为切应力.

应力-速度关系:

(8a)

(8b)

(8c)

(8d)

(8e)

其中，c^′_ij为介质的弹性系数，ζ_xx､ζ_zz､ζ_yz､ζ_xz､ζ_xy是记忆变量的导数，并且

记忆变量方程:

(9a)

(9b)

(9c)

(9d)

(9e)

其中，

v=1､2分别对应两种耗散机制;τ_σ^(v) ､τ_ε^(v) 分别是两种耗散机制的应力和应变的松弛时间;Q₁､Q₂ 分别为纵波､横波的品质因子;$\tilde \omega $为中心频率的倒数.

由式(7)､(8)､(9)构成黏弹TTI介质中波动方程组，以此为基础对黏弹TTI介质做二维三分量波场数值模拟.

3 旋转交错网格有限差分格式

图 1解释了旋转交错网格有限差分算法的二维空间交错方法.在旋转交错网格中，所有的速度分量都被定义在同一位置，而所有的应力分量位于另外的同一位置(黏弹TTI介质波动方程中的记忆变量的导数都定义在与应力分量相同的位置)，因此，也可以将介质密度和弹性模量定义在相同的位置^[19]，本文中将它们定义在与应力分量相同的位置上.

旋转交错网格有限差分法是先沿网格对角线方向计算相关物理量的空间导数，然后再将这些对角线方向的计算结果进行线性叠加来得到沿坐标轴方向的空间导数^[19].在二维空间的情形下，旋转交错网格有限差分沿坐标轴方向的差分算子的计算表达式为^[17]

(10)

(11)

(12)

(13)

式中$\tilde x$和$\tilde z$分别代表两个对角线的方向，Δx和Δz分别是沿坐标轴x､z方向的空间步长，Δr是对角线的长度，且∂/∂$\tilde x$和∂/∂$\tilde z$的离散差分算子的计算方法与普通交错网格有限差分中∂/∂$\tilde x$和∂/∂$\tilde z$的离散差分算子计算方法相同^[19].

假定${D_{\tilde x}}$和${D_{\tilde z}}$分别表示∂/∂$\tilde x$和∂/∂$\tilde z$在$\tilde x$和$\tilde z$方向的任意偶数阶偏导数，设函数u(x，z，t)连续，且具有2N+1阶导数，则

(14)

(15)

式中，N为差分阶数的一半;m=1，2，…，N;am为差分系数，其值详细推导与普通交错网格的差分系数的推导是一致的，详见参考文献^[22].

由式(12)-(15)可以推导出在x和z方向的偏导数，即

(16)

(17)

利用式(14)-(17)直接求解黏弹TTI介质的波动方程组，可以得到方程的旋转交错网格任意偶数阶精度差分格式.

结合各向异性介质中弹性波方程差分算法的稳定性条件和旋转交错网格的特性，给出黏弹TTI介质中二维一阶速度-应力方程的时间二阶精度､空间2N阶精度的旋转交错网格差分稳定性条件近似表达式^[23]，即

(18)

式中，

4 PML 吸收边界条件

Berenger于1994 年首次提出了PML 吸收边界条件，并将其运用来消除电磁波场数值模拟中的人为截断边界所产生的反射波^[24].随后许多学者将此方法引入到了地震波场的数值模拟中，并取得了比较好的效果^[25].本文将PML 引入到用旋转交错网格高阶有限差分求解的黏弹TTI介质波动方程中.在二维空间情形下，旋转交错网格有限差分法中实现PML 的思路，与普通交错网格有限差分中PML 算法的基本思路是相同的.根据PML 吸收边界条件的基本原理，采用时间域变量分裂的方法，对一阶速度-应力黏弹TTI介质波动方程进行波场变量分离，可以得到应用于黏弹TTI介质波动方程的PML 吸收边界条件.在二维空间的条件下，对于式(7)､(8)､(9)，每一个波场变量可分为两部分，即

(19)

(20)

(21)

式中上标x和z代表该项只与相应的空间导数有关.同时，可以得到分裂后的PML 吸收边界中带有衰减因子的波动方程(下面只以σxx为例，逐步推导给出其吸收层边界方程的具体形式)，即

(22)

(23)

其中，函数d(x)和d(z)是与x和z方向有关的衰减因子^[25].

在数值模拟过程中PML 吸收层被分成了不同的区域，如图 2所示，不同区域的衰减因子的取值是不同的.在b₁ 和b₃ 代表的吸收层区域，d(x)=0，d(z)>0;在b₂ 和b₄ 代表的吸收层区域，d(x)>0，d(z)=0;在a₁､a₂､a₃ 和a₄ 代表的PML 吸收层区域，d(x)>0，d(z)>0.

限于篇幅，仅以式(20)中σ_xx分量和式(22)为例进行差分离散，给出其二阶时间差分精度､高阶空间差分的旋转交错网格有限差分格式，即

(24)

(25)

${D_{\tilde x}}$和${D_{\tilde z}}$分别表示∂/∂$\tilde x$和∂/∂$\tilde z$在$\tilde x$和$\tilde z$方向的任意偶数阶偏导数式，故将式(14)和式(15)代入式(25)中，其余分裂变量采用同样的方式推导，即可求得PML边界条件的任意偶数阶旋转交错网格有限差分格式.

5 数值模拟与分析 5.1 PML吸收边界条件效果验证

为了验证PML 吸收边界的效果，设计了一个大小为2000m×2000m､方位角为45°､对称轴倾角为60°的均匀黏弹TTI介质模型，空间采样间隔Δx=Δz=5m, 时间步长为1ms.模型介质在自然坐标系下的弹性参数如表 1，模型的纵波Q值为80，横波Q值为60.震源子波采用Ricker子波，震源位置(1000m, 1000m)，主频为20Hz, 震源为x方向的横波源.

在数值模拟过程中采用二十阶空间差分精度､二阶时间差分精度的旋转交错网格有限差分格式对上述模型做二维三分量波场数值模拟.图 3(a, b, c)分别是考虑吸收边界条件时模拟得到的x､y､z三个分量的波场快照，在每个分量，四个波场快照从左到右对应的时间依次为600ms､840ms､1200ms和1800ms.从图 3可以看出，x､y､z三个分量入射到截断边界处的各类外行波都被完全吸收，这表明本文在用旋转交错网格差分法模拟黏弹TTI介质中波场时所采用的PML 边界条件对各类外行波都具有良好的吸收性能.

5.2 各向异性特征分析

设计了一个大小为2000m×2000m 均匀黏弹TI介质模型，空间采样间隔Δx=Δz=5 m, 时间步长为1ms.模型的纵波Q值为80，横波Q值为60，模型介质在自然坐标系下的弹性参数如表 1.震源子波采用Ricker子波，震源位置(1000m, 1000m)，主频为20Hz, 震源为x方向的横波源.数值模拟差分精度为:十八阶空间差分精度､二阶时间差分精度.图 4是黏弹TTI介质在方位角和对称轴倾角都为0时模拟得到的波场快照和单炮记录.图 5 是黏弹TTI介质在方位角为45°､对称轴倾角为60°时模拟得到的波场快照和单炮记录.

从图 4中的各个分量可以看出:当黏弹TTI介质的对称轴倾角和方位角都等于0时，均匀黏弹TI介中的地震波记录包括拟纵波qP 波波形和拟横波qS波波形;波场快照中qS 波波前面的形状与单炮传播记录中qS波波形初至所形成的三叉区很相似，qS波波形及同相轴三叉区的出现导致波场变得更加复杂化;且单炮记录和波场快照中y分量均为0.

从图 5中的各个分量可以看到:当黏弹TTI介质的对称轴倾角不为0时，均匀黏弹TTI介质中地震波记录的三个分量中都有qP波波形､快横波qS1波形和慢横波qS2波形，即存在横波分裂现象;从单炮传播记录中也可以看到，各类波的同相轴变得不对称，且不一定是双曲线形态，并且随着偏移距的增大，在其对称轴上倾的方向上，qP 波形的初至时间要逐渐小于下倾向上qP 波形的初至时间，而两类qS波的情况则相反^[27].

为了简明地分析黏弹TTI介质中地震波的反射特征，设计了一个大小为5000 m×3000 m 两层介质模型.模型上层为黏弹TTI介质，其方位角为45°､对称轴倾角为60°;下层为黏弹各向同性介质，空间采样间隔Δx=Δz=5 m, 时间步长为1ms.模型介质参数如表 2.震源子波采用Ricker子波，震源位置(2500m, 1200m)，主频为25 Hz, 震源为胀缩源.数值模拟差分精度为二十阶空间差分精度､二阶时间差分精度.

图 6是采用表 2中介质参数模拟得到的波场快照和单炮记录.从图 6 中可以看到，当上层黏弹TI介质的方位角和对称轴倾角都不为0时:(1)经过黏弹TTI介质与黏弹各向同性介质分界面的反射和透射后，其中地震反射波包括:拟纵波qP 波､转换的快横波qPS1 波､转换的慢横波qPS2 波､拟快横波qS1波､拟慢横波qS2波､转换的拟纵波qSP 波，且可看到明显的横波分裂现象;而在下层黏弹各向同性介质中的透射波只存在P 波､S 波､转换波SP波和转换波PS波;(2)在黏弹TTI介质中，qS波的波形初至会形成三叉区^[26]，这种特殊现象使地震波场变得更为复杂化，并特别容易引起误解;(3)黏弹TTI介质中，单炮记录上各类反射波的同相轴不再是双曲线型，且qS反射波的同相轴偏离最为明显.

5.3 黏弹性特征分析

保持介质的弹性参数不变，但改变其纵波和横波的Q值，设计四种黏弹TTI介质模型以及一个TTI弹性介质模型(Q值无穷大)，具体的介质参数如表 3.设计模型的大小为2000m×2000m､方位角为45°､对称轴倾角为60°，空间采样间隔Δx=Δz=5m, 时间步长为1ms.震源子波采用Ricker子波，震源位置(1000m, 1000m)，主频为20 Hz, 震源为x方向的横波源.利用二十阶空间差分精度､二阶时间差分精度对四个模型分别做正演模拟，取(1500m, 500m)处质点的x､y､z三个分量振动记录进行分析.图 7､图 8､图 9分别是对表 3中介质1､2､3､4模拟得到的质点(1500m, 500m)处的x､y､z分量记录.图 10是对表 3 中介质1､5､4 分别模拟得到的质点(1500m, 500 m)处的x分量记录对比图.而图 11是对表 3中介质1､5､3分别模拟得到的质点(1500m, 500m)处的x分量记录对比图.

从图 7､图 8､图 9中可以看到，无论是x､y､z任一分量，由于介质的黏滞性而使质点的振动振幅都明显衰减，并且qS 波比qP 波的衰减幅度要大，而且品质因子越小，振幅衰减的越多，接收到的子波的峰值时间越往后延迟，并且波形变宽､畸变也越严重.从图 10中可以看出，纵波品质因子Q_P 对应于膨胀滞弹性的形变，且Q_P 越小qP 波振幅衰减越严重.从图 11中可以看出，横波品质因子Q_S 对应于剪切滞弹性的形变，且Q_S 越小qS波振幅衰减越严重.

6 结 论

采用旋转交错网格高阶有限差分法模拟黏弹TTI介质中地震波场的传播过程，模拟精度较高，且易于实现.同时，将完全匹配层吸收边界条件应用到黏弹TTI介质旋转交错网格高阶有限差分模拟中，能有效地消除人工边界所产生的反射.

对模拟得到的波场进行分析表明:

(1) 黏弹TTI介质的各向异性主要影响波前面的形态，并且波前面形态复杂;qS 波波前面和同相轴的三分叉现象普遍，且其同相轴一般不再是双曲线型;当TI介质倾斜时，三个分量上均能够观测到横波分裂现象，而且各类波形的同相轴变得不对称;并且TI介质倾斜时，地震波通过界面反射后，其反射特征变得更为复杂，能够观测到反射横波的分裂现象.

(2) 黏弹TTI介质的黏弹性主要影响地震波的振幅衰减和波形畸变,对应于膨胀滞弹性形变的纵波品质因子主要影响qP 波的振幅衰减和波形畸变，对应于剪切滞弹性形变的横波品质因子主要影响qS波的振幅衰减和波形畸变;而且，当品质因子越小时，地震波振幅衰减的越多，接收到的子波的峰值时间越往后延迟，且其波形变宽､畸变也越严重.