VTI介质弹性波的有限差分正演模拟

Finite-Difference Forward Simulation of Elastic Waves in VTI Media

摘要: 本文以弹性参数导出了二阶位移波动方程，引入Thomsen参数使得弹性参数带有更直观的物理意义，进而以有限差分法对波动方程进行差分给出了波动方程二阶偏导数的时间二阶，空间二阶差分方程，讨论了均匀空间与三层模型中两种各项异性参数对波场特征的影响。

Abstract: This article derives the second-order displacement wave equation using elastic parameters, and introduces the Thomsen parameter to give the elastic parameters a more intuitive physical meaning. Then, the finite difference method is used to give the time second-order and spatial second-order difference equations of the second-order partial derivatives of the wave equation. The influence of two anisotropic parameters on the wave field characteristics in uniform space and three-layer models is discussed.

1. 引言

地震波数值模拟是通过人为设定所需要正演模型的各项参数以及设置震源，进行激发产生地震波，记录所产生的地震记录并进行分析，与反演相结合可以预测地下实际构造，并做出地下构造模型。地震波正演模拟研究的是地下介质中地震波的传播规律和岩石的性质。并且在研究地震资料方案、处理数据和模拟验证的研究等方面都起到了重大作用，由此可知，正演模拟是一个非常有效的手段。但是由于实际地质条件复杂、各向异性种类繁多，难以完全的符合现有的理论基础，对于研究存在着不可逾越的客观现实。在地震波正演模拟研究中，地震学家研究的主要还是横向各项同性(TI)介质中地震波传播的波场特征规律以及各种介质参数之间的关系。对地震勘探最有实际意义的就是横向各向同性介质，这是由于组成地球介质的岩石通常为层理结构。一般认为，在层理平面方向上，岩石的弹性性质相同；而在垂直于层理方向上，岩石的弹性性质是变化的。VTI介质就是具有垂直对称轴的横向各向同性介质。研究VTI介质中弹性波有限差分正演模拟具有重要的意义[1]-[3]。

TI介质示意图和正演流程图分别见图1和图2：

Figure 1. Schematic diagram of TI medium

图1. TI介质示意图

Figure 2. Forward simulation flowchart

图2. 正演模拟流程图

2. VTI介质中弹性波的波动方程

2.1. 本构方程(也称广义虎克定律)

虎克(Hooke)定律是表示应力与应变之间的关系，是地下介质中物体固有的性质，在弹性范围内，基于柯西弹性和格林弹性两大理论确定的，其表达为：

${\sigma }_{ij}={\displaystyle \sum _{k,l=1}^3{C}_{ijkl}}{e}_{kl}$ (1)

式中${\sigma }_{ij}$ 和$e_{ij}$ 分别为应力张量和应变张量；$C_{ijkl}$ 是一种比例系数，是由弹性介质性质所决定的弹性常数。其中$i,j,k,l$ 均取1到3的整数，代表$x,y,z$ 三个坐标。${\sigma }_{ij}$ 和$e_{ij}$ 各有九个分量，$C_{ijkl}$ 是一个具有81个分量的四阶张量C。

${\sigma }_{ij}$ 和$e_{ij}$ 都是对称张量，只有六个分量相互独立，所以弹性常数张量C变为具有36个分量的矩阵。为了简便本文采用了如下对应关系：

$xx\to 1,yy\to 2,zz\to 3,yz\left(zy\right)\to 4,zx\left(xz\right)\to 5,xy\left(yx\right)\to 6$

于是本构方程的数学表达式可展开成如下形式：

$\left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right]=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& C_{14}& C_{15}& C_{16}\\ C_{21}& C_{22}& C_{23}& C_{24}& C_{25}& C_{26}\\ C_{31}& C_{32}& C_{33}& C_{34}& C_{35}& C_{36}\\ C_{41}& C_{42}& C_{43}& C_{44}& C_{45}& C_{46}\\ C_{51}& C_{52}& C_{53}& C_{54}& C_{55}& C_{56}\\ C_{61}& C_{62}& C_{63}& C_{64}& C_{65}& C_{66}\end{array}\right]\left[\begin{array}{c}{e}_1\\ e_2\\ e_3\\ e_4\\ e_5\\ e_6\end{array}\right]$ (2)

弹性系数$C_{ij}$ 的$i,j$ 均取1到6的自然数，共36个，表示在该点上弹性体的性质。刚性系数$C_{ij}$ 必须等于$C_{ji}$ ，描述一个复杂的各向异性介质所需弹性系数总数减至21个。当某个弹性介质的性质具有对称性时，描述介质所需的弹性系数将进一步减少。地球物理学家已经证明实际地下结构由于地质年代沉积而表现出来的是横向各向同性介质(即TI介质)，即在某个方向表现为各向同性而在另一个方面表现为各向异性[1]。

针对具有垂直对称轴横向各向同性介质(VTI介质)进行研究，经近似后其弹性矩阵中只有5个不同的弹性参数，简化本构方程如下：

$\left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right]=\left[\begin{array}{cccccc}{C}_{11}& C_{11}-2C_{66}& C_{13}& 0& 0& 0\\ C_{11}-2C_{66}& C_{11}& C_{13}& 0& 0& 0\\ C_{13}& C_{13}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right]\left[\begin{array}{c}{e}_1\\ e_2\\ e_3\\ e_4\\ e_5\\ e_6\end{array}\right]$ (3)

其中，$C_{11}=C_{22},C_{12}=C_{21},C_{13}=C_{23}=C_{31}=C_{32},C_{44}=C_{55},C_{12}=C_{11}-2C_{66}$ 。

化简(3)式得：

$\{\begin{array}{l}{\sigma }_1=C_{11}{e}_1+\left(C_{11}-2C_{66}\right)e_2+C_{13}{e}_3\\ {\sigma }_2=\left(C_{11}-2C_{66}\right)e_1+C_{11}{e}_2+C_{13}{e}_3\\ {\sigma }_3=C_{13}{e}_1+C_{13}{e}_2+C_{33}{e}_3\\ {\sigma }_4=C_{44}{e}_4\\ {\sigma }_5=C_{44}{e}_5\\ {\sigma }_6=C_{66}{e}_6\end{array}$ (4)

2.2. 柯西方程(Cauchy方程)

设质点的位移为，$u=u\left(u_x,u_y,u_z\right)$ ，$u_x,u_y,u_z$ 分别为$x,y,z$ 方向上的位移分量，$e$ 为应变，两者的关系如下：

$e_1=\frac{\partial u_x}{\partial x},e_2=\frac{\partial u_y}{\partial y},e_3=\frac{\partial u_z}{\partial z},e_4=\frac{\partial u_y}{\partial z}+\frac{\partial u_z}{\partial y},e_5=\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x},e_6=\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}$ (5)

其中，$e_1,e_2,e_3$ 分别为$x,y,z$ 方向上的正应变，$e_{\text{4}},e_5,e_6$ 为切应变。

2.3. 运动微分方程(Navier方程)

设外力$F=F\left(f_x,f_y,f_z\right)$ 作用在弹性介质微小应变的情况下，其满足牛顿第二定律，运动微分方程如下：

$\{\begin{array}{l}\frac{\partial {\sigma }_1}{\partial x}+\frac{\partial {\sigma }_4}{\partial y}+\frac{\partial {\sigma }_5}{\partial z}+\rho F_x=\rho \frac{{\partial }^2{u}_x}{\partial t^2}\\ \frac{\partial {\sigma }_6}{\partial x}+\frac{\partial {\sigma }_2}{\partial y}+\frac{\partial {\sigma }_4}{\partial z}+\rho F_y=\rho \frac{{\partial }^2{u}_y}{\partial t^2}\\ \frac{\partial {\sigma }_5}{\partial x}+\frac{\partial {\sigma }_4}{\partial y}+\frac{\partial {\sigma }_3}{\partial z}+\rho F_z=\rho \frac{{\partial }^2{u}_z}{\partial t^2}\end{array}$ (6)

在忽略外力作用情况下，方程如下：

$\{\begin{array}{l}\frac{\partial {\sigma }_1}{\partial x}+\frac{\partial {\sigma }_4}{\partial y}+\frac{\partial {\sigma }_5}{\partial z}=\rho \frac{{\partial }^2{u}_x}{\partial t^2}\\ \frac{\partial {\sigma }_6}{\partial x}+\frac{\partial {\sigma }_2}{\partial y}+\frac{\partial {\sigma }_4}{\partial z}=\rho \frac{{\partial }^2{u}_y}{\partial t^2}\\ \frac{\partial {\sigma }_5}{\partial x}+\frac{\partial {\sigma }_4}{\partial y}+\frac{\partial {\sigma }_3}{\partial z}=\rho \frac{{\partial }^2{u}_z}{\partial t^2}\end{array}$ (7)

2.4. 二阶位移弹性波动方程

结合本构方程(4)、柯西方程(5)和运动微分方程(7)，推导二阶位移弹性波动方程，当只考虑x，z方向时，对y的偏导数都为0，整理化简化得：

$\{\begin{array}{l}{C}_{11}\frac{{\partial }^2{u}_x}{\partial x^2}+C_{44}\frac{{\partial }^2{u}_x}{\partial z^2}+\left(C_{13}+C_{44}\right)\frac{{\partial }^2{u}_z}{\partial x\partial z}=\rho \frac{{\partial }^2{u}_x}{\partial t^2}\\ C_{44}\frac{{\partial }^2{u}_z}{\partial x^2}+C_{33}\frac{{\partial }^2{u}_z}{\partial z^2}+\left(C_{13}+C_{44}\right)\frac{{\partial }^2{u}_x}{\partial x\partial z}=\rho \frac{{\partial }^2{u}_z}{\partial t^2}\end{array}$ (8)

这就是不考虑外力作用下的二阶位移弹性波动方程。

3. 波动方程的有限差分

3.1. 各向异性介质的Thomsen参数

3.1.1. VTI介质的Thomsen参数

VTI介质的弹性系数矩阵表达式如下：

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{11}-2C_{66}& C_{13}& 0& 0& 0\\ C_{11}-2C_{66}& C_{11}& C_{13}& 0& 0& 0\\ C_{13}& C_{13}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right]$ (9)

弹性矩阵参数没有直接的物理意义，由弹性矩阵推导出来的波动方程不能很直接的表现其各向异性性质，所以为了编程和研究带有直观物理意义的参数对地震波传播规律的影响，本文引用表示各向异性的Thomsen参数[4]。

$V_{P0}=\sqrt{\frac{C_{33}}{\rho }},V_{S0}=\sqrt{\frac{C_{44}}{\rho }},\gamma =\frac{C_{66}-C_{44}}{2C_{44}},\epsilon =\frac{C_{11}-C_{33}}{2C_{33}},\delta =\frac{{\left(C_{13}+C_{44}\right)}^2-{\left(C_{33}-C_{44}\right)}^2}{2C_{33}\left(C_{33}-C_{44}\right)}$ (10)

式中，$\rho$ 是介质密度；$V_{P0}$ 表示准P波的垂直速度；$V_{S0}$ 表示准S波的垂直速度；$\epsilon ,\gamma ,\delta$ 表示各向异性强度的三个无量纲因子。$\epsilon$ 是准P波各向异性参数，$\gamma$ 是准S波各向异性参数。$\delta$ 是连接准p波和准s波的一个参数，表明垂直方向上的各向异性性质变化。对于各向同性介质中$\epsilon ,\gamma ,\delta$ 为0。

将式(10)变形得到用Thomsen参数表示的VTI介质弹性系数的表达式

$\{\begin{array}{l}{C}_{11}=C_{22}=\rho \left(1+2\epsilon \right)V_{p0}^2\\ C_{33}=\rho V_{p0}^2\\ C_{44}=C_{55}=\rho V_{S0}^2\\ C_{66}=\rho \left(1+2\gamma \right)V_{S0}^2\\ C_{13}=C_{23}=\rho \sqrt{\left(V_{P0}^2-V_{S0}^2\right)\left(V_{P0}^2-V_{S0}^2+2\delta *V_{P0}^2\right)}-\rho V_{S0}^2\end{array}$ (11)

3.1.2. VTI介质的Thomsen参数的限制条件

通过(10)式我们知道了Thomsen参数中各分量的含义，然而在对各向异性介质进行资料处理时，Thomsen参数数值的选取必须符合物理学定律和实际地质情况，否则可能会因为参数选取不恰当而产生不符合实际情况乃至错误的结果。根据热力学定律和弹性常数的物理意义可以得出VTI介质能够稳定存在的充要条件是介质的弹性系数矩阵式为正定矩阵，即弹性系数矩阵(9)式的所有顺序主子式都大于零：

$\{\begin{array}{l}{C}_{11}>C_{66}>0\\ C_{44}>0\\ C_{33}\left(C_{11}-C_{66}\right)>C_{13}^2\end{array}$ (12)

理论上只要介质的弹性系数满足(12)，那么该介质就可能存在，但是根据实验室测量数据的，真实的限制条件比这个还要小，这里采用李磊[5] [6]根据弹性系数各物理分量提出来的限制条件。

$\{\begin{array}{l}0<2C_{66}<C_{11}\\ 0<C_{13}<\min \left(C_{11},C_{33}\right)\\ 0<C_{44}<\frac{3}{4}\min \left(C_{11},C_{33}\right)\end{array}$ (13)

由$C_{\text{11}}>0$ 得到

$\epsilon >-\text{0}.5$ (14)

各向异性参数$\delta$ 可以写成：

$\delta ={\delta }_L+\frac{{\delta }_L^2}{2f}=\frac{{\left({\delta }_L+f\right)}^2}{2f}-\frac{f}{2}$ (15)

式中${\delta }_L$ 是$\delta$ 的线性定义：

${\delta }_L=\frac{C_{13}+2C_{44}-C_{33}}{C_{33}}$ (16)

f的值介于0和1之间，实际勘探中波速比V_p0/V_s0范围大概在1.5到2.5之间。

$f=1-\frac{V_{S0}^2}{V_{P0}^2}$ (17)

由公式(15)可以看出$\delta \ge {\delta }_L$ ，但两者必定同号，且差别很小，通常可以忽略。联立公式(11)~(17)得到各向异性参数的限制条件如下：

$\{\begin{array}{l}\frac{1}{4}<f<1\\ \epsilon >-\frac{f}{2}\\ -\frac{1}{2}<\gamma <\frac{1+2\epsilon }{4\left(1-f\right)}-\frac{1}{2}\\ 1-2f<{\delta }_L<2\left(1-f\right)\\ \frac{1}{2f}-1<\delta <2\left(\frac{1}{f}-1\right)\end{array}$ (18)

从(18)式可以看出$\epsilon ,\gamma ,\delta$ 的取值主要受到准p波和准s波的参考速度比V_p0/V_s0的约束，其中$\gamma$ 的取值上限受$\epsilon$ 的限制，波速比V_p0/V_s0、$\gamma$ 和$\delta$ 都有上下限，而$\epsilon$ 没有上限。在大多数情概况下，$\epsilon$ 一般在−0.05到0.45之间，$\delta$ 一般在−0.2到0.3之间，$\gamma$ 一般在−0.05到0.2之间。

本文在参数对比环节只用到了$\epsilon$ 和$\delta$ ，通过固定其中一个值，改变另一个值，来判断其对准p波和准s波各向异性的影响，同时也可以同时增大两个值，进行对比，看在增大过程中，哪一个参数对各向异性的影响更明显[7]。

3.2. VTI介质网格法有限差分–时间差分

有限差分法是以差分原理为基础的一种数值计算方法。它用网格把点离散化，用各离散点上的函数的差商来近似替代该点的偏导数，把微分方程问题转化为一组相应的差分方程。然后解出差分方程在各离散点上的函数值，便可推出微分方程的数值解。本文采用差分离散方法是将二阶波动方程直接转化为二阶的时间空间中心差分进行离散[8]-[10]。

对于时间差分，设$u\left(x,z,t\right)$ ，当把空间分量看作不变时，讨论对时间的偏导数，这里只讨论$u_x$ 分量上的差分形式，$u_z$ 分量上可以进行类似形式推导。

$f\left(x\right)=f\left(x_0\right)+f^{\prime }\left(x_0\right)\left(x-x_0\right)+\frac{f^{″}\left(x_0\right)}{2!}{\left(x-x_0\right)}^2+\cdots +\frac{f^{\left(n\right)}\left(x_0\right)}{n!}{\left(x-x_0\right)}^n+\frac{f^{\left(n+1\right)}\left(\xi \right)}{\left(n+1\right)!}{\left(x-x_0\right)}^{\left(n+1\right)}$ (19)

上式是$f\left(x\right)$ 在$x_0$ 处的泰勒展开式，其中$\zeta$ 取值x到$x_0$ 之间，然后把$u_x\left(t+\Delta t\right)$ 和$u_x\left(t-\Delta t\right)$ 按泰勒公式展开。

$u_x\left(t+\Delta t\right)=u_x\left(t\right)+\frac{\partial u_x\left(t\right)}{\partial t}\Delta t+\frac{1}{2!}\frac{{\partial }^2{u}_x\left(t\right)}{\partial t^2}{\left(\Delta t\right)}^2+\frac{1}{3!}\frac{{\partial }^3{u}_x\left(t\right)}{\partial t^3}{\left(\Delta t\right)}^3+\frac{1}{4！}\frac{{\partial }^4{u}_x\left(t\right)}{\partial t^4}{\left(\Delta t\right)}^4+\cdots$ (20)

$u_x\left(t-\Delta t\right)=u_x\left(t\right)-\frac{\partial u_x\left(t\right)}{\partial t}\Delta t+\frac{1}{2!}\frac{{\partial }^2{u}_x\left(t\right)}{\partial t^2}{\left(\Delta t\right)}^2-\frac{1}{3!}\frac{{\partial }^3{u}_x\left(t\right)}{\partial t^3}{\left(\Delta t\right)}^3+\frac{1}{4！}\frac{{\partial }^4{u}_x\left(t\right)}{\partial t^4}{\left(\Delta t\right)}^4+\cdots$ (21)

联立式(20)和式(21)得到$u_x$ 在时间上的二阶微分：

$u_x\left(t+\Delta t\right)=2u_x\left(t\right)-u_x\left(t-\Delta t\right)+\frac{{\partial }^2{u}_x\left(t\right)}{{\partial }^2{t}^2}{\left(\Delta t\right)}^2+\frac{2}{4!}\frac{{\partial }^4{u}_x\left(t\right)}{{\partial }^{4}t}{\left(\Delta t\right)}^4+\cdots$ (22)

利用有限差分算法进行地震波正演模拟时，采用的时间层数越多，需要的计算机内存就会越大，在内存与时间允许的情况下，在正演模拟中可以使用四阶或者更高阶的有限差分进行近似计算，在综合考虑计算精度、计算时间和计算内存多方面的情况下，本文采用时间二阶差分进行正演模拟，格式如下：

$u_x\left(t+\Delta t\right)=2u_x\left(t\right)-u_x\left(t-\Delta t\right)+\frac{{\partial }^2{u}_x\left(t\right)}{\partial t^2}{\left(\Delta t\right)}^2+O{\left(\Delta t\right)}^4$ (23)

3.3. VTI介质网格法有限差分–空间差分

波动方程中存在空间二阶偏导数和混合偏导数，在对其差分时，稍微复杂一点，这里只需要先对x方向上的各个偏导数进行推导，然后可类比出z方向上的偏导数。

3.3.1. 空间二阶偏导数二阶差分近似推导

同样设$u\left(x,z,t\right)$ 将$u\left(x+\Delta x\right)$ 与$u\left(x-\Delta x\right)$ 在空间点按泰勒展开：

$u\left(x+\Delta x\right)=u\left(x\right)+\frac{\partial u}{\partial x}\Delta x+\frac{1}{2!}\frac{{\partial }^{2}u}{\partial x^2}{\left(\Delta x\right)}^2+\frac{1}{3!}\frac{{\partial }^{3}u}{\partial x^3}{\left(\Delta x\right)}^3+\frac{1}{4!}\frac{{\partial }^{4}u}{\partial x^4}{\left(\Delta x\right)}^4+\cdots$ (24)

$u\left(x-\Delta x\right)=u\left(x\right)-\frac{\partial u}{\partial x}\Delta x+\frac{1}{2!}\frac{{\partial }^{2}u}{\partial x^2}{\left(\Delta x\right)}^2-\frac{1}{3!}\frac{{\partial }^{3}u}{\partial x^3}{\left(\Delta x\right)}^3+\frac{1}{4!}\frac{{\partial }^{4}u}{\partial x^4}{\left(\Delta x\right)}^4+\cdots$ (25)

联立式(24)和式(25)得：

$u\left(x+\Delta x\right)=2u\left(x\right)-u\left(x-\Delta x\right)+\frac{{\partial }^{2}u\left(x\right)}{\partial x^2}{\left(\Delta x\right)}^2+\frac{2}{4！}\frac{{\partial }^{4}u}{\partial x^4}{\left(\Delta x\right)}^4+\cdots$ (26)

则x方向的二阶偏微分为：

$\frac{{\partial }^2{u}_x}{\partial x^2}=\frac{1}{{\left(\Delta x\right)}^2}\left\{\left[u_x\left(x+\Delta x\right)-2u_x\left(x\right)+u_x\left(x-\Delta x\right)\right]+O{\left(\Delta x\right)}^4+\cdots \right\}$ (27)

$\frac{{\partial }^2{u}_x}{\partial z^2}=\frac{1}{{\left(\Delta z\right)}^2}\left\{\left[u_x\left(z+\Delta z\right)-2u_x\left(z\right)+u_x\left(z-\Delta z\right)\right]+O{\left(\Delta x\right)}^4+\cdots \right\}$ (28)

同理可得z方向上的二阶偏微分为

$\frac{{\partial }^2{u}_z}{\partial x^2}=\frac{1}{{\left(\Delta x\right)}^2}\left\{\left[u_z\left(x+\Delta x\right)-2u_z\left(x\right)+u_z\left(x-\Delta x\right)\right]-O{\left(\Delta x\right)}^4+\cdots \right\}$ (29)

$\frac{{\partial }^2{u}_z}{\partial z^2}=\frac{1}{{\left(\Delta z\right)}^2}\left\{\left[u_z\left(z+\Delta z\right)-2u_z\left(z\right)+u_z\left(z-\Delta z\right)\right]-O{\left(\Delta z\right)}^4+\cdots \right\}$ (30)

3.3.2. 空间混合偏导数二阶差分近似推导

对于空间混合偏导数，其计算有些复杂，需要通过分两步分别求取一阶偏导数，先对x求偏导数，然后再对z求偏导数。空间一阶偏导数如下：

$\{\begin{array}{l}\frac{\partial u_x}{\partial z}=\frac{1}{2\Delta z}\left[u_x\left(z+\Delta z\right)-u_x\left(z-\Delta z\right)\right]+O{\left(\Delta z\right)}^2\\ \frac{\partial u_x}{\partial x}=\frac{1}{2\Delta x}\left[u_x\left(x+\Delta x\right)-u_x\left(x-\Delta x\right)\right]+O{\left(\Delta x\right)}^2\end{array}$ (31)

$\{\begin{array}{l}\frac{\partial u_z}{\partial x}=\frac{1}{2\Delta x}\left[u_z\left(x+\Delta x\right)-u_z\left(x-\Delta x\right)\right]+O{\left(\Delta x\right)}^2\\ \frac{\partial u_{\text{z}}}{\partial z}=\frac{1}{2\Delta z}\left[u_z\left(z+\Delta z\right)-u_{\text{z}}\left(z-\Delta z\right)\right]+O{\left(\Delta z\right)}^2\end{array}$ (32)

由式(31)和式(32)推导的空间二阶混合偏导数如下

$\frac{{\partial }^2{u}_x}{\partial x\partial z}=\frac{1}{4\Delta x\Delta z}\left[u_x\left(x+\Delta x,z+\Delta z\right)-u_x\left(x+\Delta x,z-\Delta z\right)-u_x\left(x-\Delta x,z+\Delta z\right)+u_x\left(x-\Delta x,z-\Delta z\right)\right]$ (33)

$\frac{{\partial }^2{u}_z}{\partial x\partial z}=\frac{1}{4\Delta x\Delta z}\left[u_z\left(x+\Delta x,z+\Delta z\right)-u_z\left(x+\Delta x,z-\Delta z\right)-u_z\left(x-\Delta x,z+\Delta z\right)+u_z\left(x-\Delta x,z-\Delta z\right)\right]$ (34)

3.4. 波动方程的有限差分方程

根据前面的波动方程存在空间二阶偏导数和混合偏导数，为了求解得到满足有限差分法的波动方程的差分形式，我们需要对波动方程进行差分，求得其有限差分形式，这里只需要对x方向上的各个偏导数进行推导，即可类比出z方向上的偏导数。推导的波动方程二阶偏导数的有限差分时间二阶，空间二阶差分格式如下：

$\begin{array}{c}{u}_x\left(x,z,t+\Delta t\right)=2u_x\left(x,z,t\right)-u_x\left(x,z,t-\Delta t\right)\\ +\frac{C_{11}{\left(\Delta t\right)}^2}{\rho {\left(\Delta x\right)}^2}\left[u_x\left(x+\Delta x,z,t\right)-2u_x\left(x,z,t\right)+u_x\left(x-\Delta x,z,t\right)\right]\\ +\frac{C_{44}{\left(\Delta t\right)}^2}{\rho {\left(\Delta z\right)}^2}\left[u_x\left(x,z+\Delta z,t\right)-2u_x\left(x,z,t\right)+u_x\left(x,z-\Delta z,t\right)\right]\\ +\frac{\left(C_{13}+C_{44}\right){\left(\Delta t\right)}^2}{4\rho \left(\Delta x\right)\left(\Delta z\right)}[u_z\left(x+\Delta x,z+\Delta z,t\right)-u_z\left(x+\Delta x,z-\Delta z,t\right)\\ -u_z\left(x-\Delta x,z+\Delta z,t\right)+u_z\left(x-\Delta x,z-\Delta z,t\right)]\end{array}$ (35)

$\begin{array}{c}{u}_z\left(x,z,t+\Delta t\right)=2u_z\left(x,z,t\right)-u_z\left(x,z,t-\Delta t\right)\\ +\frac{C_{11}{\left(\Delta t\right)}^2}{\rho {\left(\Delta x\right)}^2}\left[u_z\left(x+\Delta x,z,t\right)-2u_z\left(x,z,t\right)+u_z\left(x-\Delta x,z,t\right)\right]\\ +\frac{C_{44}{\left(\Delta t\right)}^2}{\rho {\left(\Delta z\right)}^2}\left[u_z\left(x,z+\Delta z,t\right)-2u_z\left(x,z,t\right)+u_z\left(x,z-\Delta z,t\right)\right]\\ +\frac{\left(C_{13}+C_{44}\right){\left(\Delta t\right)}^2}{4\rho \left(\Delta x\right)\left(\Delta z\right)}[u_x\left(x+\Delta x,z+\Delta z,t\right)-u_x\left(x+\Delta x,z-\Delta z,t\right)\\ -u_x\left(x-\Delta x,z+\Delta z,t\right)+u_x\left(x-\Delta x,z-\Delta z,t\right)]\end{array}$ (36)

4. 正演模拟

4.1. 震源——雷克子波

在地震勘探中常常使用炸药作为震源，而在数值模拟中，使用雷克子波来模拟震源，可以很好地达到同样的效果[11]。本文使用的雷克子波是带有空间与时间上的函数正好对应于炸药的传播过程。其表达式如下：

$f\left(t\right)=\frac{1}{4}\left[1-2{\pi }^2{f}_m^2{\left(t-\frac{1}{f_m}\right)}^2\right]{\text{e}}^{-{\pi }^2{f}_m^2{\left(t-\frac{1}{f_m}\right)}^2}$ (37)

上式中$f_m$ 为雷克子波主要频率，本文采用的雷克子波如图3所示，频率为25 Hz，时间域采样间隔为0.0001 s，采样点为2000个点，总时间长度为0.2 s。

Figure 3. Ricker wave

图3. 雷克子波

4.2. 边界问题

当进行正演数值模拟时，我们通常用的人为设置的模型是有边界的，当地震波传播到边界时，就会出现反射波，通常这种反射波，并不是想要的，于是就要建立边界条，设置吸收边界函数，来吸收处理掉这种干扰波。通过查阅文献，了解到前人提出了许多改进方法总的分为两类[12] [13]。1) 通过扩大原有的边界，来计算需要模拟的地方。2) 通过增加效果良好的吸收边界来消除反射波的影响。第一种方法是，通过扩大计算范围，来达到需要的模拟效果对计算量要求较高，第二种方法则是加上吸收函数，当地震波传播到这个地方了，边界就给吸收掉，更加适用于现在的需求，应用较广。本文并未加以边界处理。

4.3. 数值模拟

4.3.1. 均匀VTI介质弹性波波场分析

Table 1. Model parameters for a uniform VTI medium in one layer

表1. 一层均匀VTI介质模型参数

模型序号/参数	Vp0/(m/s)	Vs0/(m/s)	$\epsilon$	$\delta$	$\rho$ /(g/cm^3)
模型1	3000	1500	0.25	−0.1	1
模型2	3000	1500	0.25	0	1
模型3	3000	1500	0.25	0.1	1
模型4	3000	1500	−0.05	0.25	1
模型5	3000	1500	0	0.25	1
模型6	3000	1500	0.25	0.25	1
模型7	3000	1500	0	0	1

上表1模型大小均为512 m × 512 m，其空间采样间隔$\Delta x=\Delta z=2\text{m}$ ，其时间采样间隔为$\Delta t=0.1\text{ms}$ ，震源位置于(256 m, 256 m)处，震源采用25赫兹的雷克子波。图4即是一层均匀VTI介质的波场快照。

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 4. (a) Wavefield snapshot of VTI medium model 1 at 180 ms (The left represents the X component, and the right represents the Y component); (b) Wavefield snapshot of VTI medium model 2 at 180 ms (The left represents the X component, and the right represents the Y component); (c) Wavefield snapshot of VTI medium model 3 at 180 ms (The left represents the X component, and the right represents the Y component); (d) Wavefield snapshot of VTI medium model 4 at 180 ms (The left represents the X component, and the right represents the Y component); (e) Wavefield snapshot of VTI medium model 5 at 180 ms (The left represents the X component, and the right represents the Y component); (f) Wavefield snapshot of VTI medium model 6 at 180 ms (The left represents the X component, and the right represents the Y component); (g) Wavefield snapshot of VTI medium model 7 at 180 ms (The left represents the X component, and the right represents the Y component)

图4. (a) VTI介质模型1在180 ms时刻波场快照(左为X分量，右为Y分量)；(b) VTI介质模型2在180 ms时刻波场快照(左为X分量，右为Y分量)；(c) VTI介质模型3在180 ms时刻波场快照(左为X分量，右为Y分量)；(d) VTI介质模型4在180 ms时刻波场快照(左为X分量，右为Y分量)；(e) VTI介质模型5在180 ms时刻波场快照(左为X分量，右为Y分量)；(f) VTI介质模型6在180 ms时刻波场快照(左为X分量，右为Y分量)；(g) VTI介质模型7在180 ms时刻波场快照(左为X分量，右为Y分量)

模型1~模型3对比可以看出：当$\epsilon$ 不变，改变$\delta$ 时，当$\delta$ 值由负值向0逐渐增大时，准p波表现出由菱形向圆形的转变。准S波在水平分量上能量弱逐渐增强，垂直分量能量变化相对较小，但形状由竖向椭圆逐渐变为圆形，在两个分量上出现的三叉区由大变小；当$\delta$ 值由0增大时，准p波在水平和竖直方向上的能量和形状变化较小，而准s波在水平和竖直方向上的变化较大，由圆形向菱形转变，并出现三叉区。

从模型4~6对比可以看出，当固定$\delta$ 值，改变$\epsilon$ 值时，$\epsilon$ 值大于0时，准p波水平和垂直方向上呈现圆形向椭圆的转变。准s波在波场快照的两条对角线上出现三叉区，在形状上出现圆形到椭圆的转变。当$\epsilon$ 为负值时，准s波的波前面发生明显变化，出现横波分裂的三叉区现象。

模型7是$\epsilon =\delta =\gamma =0$ ，准p波和准s波都无各向异性，该介质是各向同性介质用来做对比参照。对比其他模型可以验证，纵波源在各向同性介质中只激发纵波，而在各向异性介质中即激发准s波也激发准p波。

综合模型1~7来看，在VTI介质中准s波和准p波的包络面发生改变，相对于各向同性介质，从圆形变成椭圆，说明各向异性介质中弹性波的速度值随波传播的方向变化而变化。

当$\delta$ 值为正时，$\epsilon$ 值对准p波影响较大，对准s波影响较小；$\delta$ 为负值，准s波会出现三叉区，$\delta$ 为正值则不会出现三叉区。当$\delta$ 值为正，且$\epsilon$ 值最小时，准s波和准p波的各向异性相对不明显。当$\epsilon$ 值最大，$\delta$ 值最小，准s波和准p波的各向异性最为明显[14]-[16]。

4.3.2. 均匀层状波场分析

Table 2. Medium parameters of three-layer horizontal layered model

表2. 三层水平层状模型介质参数

	层数	层厚/(m)	Vp0/(m/s)	Vs0/(m/s)	$\epsilon$	$\delta$	$\rho$ /(g/cm^3)
各向同性介质	第一层	200	2000	1500	0	0	2.0
	第二层	200	2500	1800			2.2
	第三层	200	3000	2200			2.5
VTI介质	第一层	200	2000	1500	0.2	0.15	2.0
	第二层	200	2500	1800			2.2
	第三层	200	3000	2200			2.5

表2中模型大小均为800 m × 600 m，其空间采样间隔$\Delta x=\Delta z=2\text{m}$ ，其时间采样间隔为$\Delta t=0.1\text{ms}$ ，震源位置于(400 m, 300 m)处，震源采用25赫兹的雷克子波。

Figure 5. Wavefield snapshot of isotropic media at 140 ms (The left represents the X component, and the right represents the Z component)

图5. 各项同性介质140 ms时刻波场快照(左为X分量，右为Z分量)

Figure 6. Wavefield snapshot of VTI medium at 140 ms (The left represents the X component, and the right represents the Z component)

图6. VTI介质140 ms时刻波场快照(左为X分量，右为Z分量)

从图5和图6中可以清晰地看出横波与纵波在遇到分界面时会发生透射与反射，在传播时的基本规律与均匀空间相符合。

5. 结论

综上所述，可以明显看出弹性波在VTI介质中各项参数对其的影响，通过控制变量法结合波长快照得到：准s波对参数$\delta$ 的较敏感，改变其值后准s波的各向异性较明显，准p波对参数$\epsilon$ 较敏感，改变其值后准p波的各向异性较明显。弹性波场是纵横波耦合的波场，在实际中采集的波场是非常复杂，弹性波的波场及波前面会随着各向异性参数变化而变化，各向异性介质中弹性波在遇到分界面时会产生反射和透射等现象，本文只是初步分析了均匀空间和三层模型的两个各项异性参数对波场特征的影响，为波场特征分析提供了参考，后面可对比纵横波的波速分别对波场特征的影响，并以加pml边界面优化。