﻿ 时间域广义2.5D一阶波动方程曲网格有限差分法数值模拟
 石油地球物理勘探  2021, Vol. 56 Issue (6): 1262-1278  DOI: 10.13810/j.cnki.issn.1000-7210.2021.06.008 0
 文章快速检索 高级检索

### 引用本文

YANG ShangBei, BAI ChaoYing, ZHOU Bing. Curvilinear-grid finite-difference numerical simulation method for generalized first-order 2.5D time-domain wave equation. Oil Geophysical Prospecting, 2021, 56(6): 1262-1278. DOI: 10.13810/j.cnki.issn.1000-7210.2021.06.008.

### 文章历史

① 中国地质调查局西安地质调查中心, 陕西西安 710054;
② 长安大学地质工程与测绘学院地球物理系, 陕西西安 710054;
③ 哈利法科学与技术大学地球科学系, 阿布扎比 2533

Curvilinear-grid finite-difference numerical simulation method for generalized first-order 2.5D time-domain wave equation
YANG ShangBei , BAI ChaoYing , ZHOU Bing
① Xi'an Centre of Geological Survey, China Geological Survey, Xi'an, Shaanxi 710054, China;
② Institute of Geophysics, School of Geological Engineering and Geomatics, Chang'an University, Xi'an, Shaanxi 710054, China;
③ Department of Earth Sciences, Khalifa University of Science and Technology, Abu Dhabi 2533, UAE
Abstract: The 2.5D seismic wavefield numerical stimulation employs the point source in 2D geological models to calculate 3D seismic wavefields. In this paper, we present a generalized 2.5D first-order time-domain wave equation that can be applied to different media (acoustic isotropic, elastic isotropic, and elastic anisotropic) and various boundary conditions (acoustic free-surface, solid free-surface, and solid-liquid boundary). The wave equation is solved by a curvilinear-grid finite-difference method. A comparison of 2.5D numerical solutions, 3D analytic solutions, and 3D numerical solutions in different homogeneous medium models (acoustic isotropic, elastic isotropic, and elastic anisotropic) verifies the correctness of the derived equation and the numerical solution method. It also demonstrates that compared with the 3D nume-rical method, the 2.5D numerical method has great advantages in calculation efficiency and memory footprint. The 2D numerical solutions cannot be applied directly in that they suffer significant amplitude distortion and phase shifts due to an artificial line source applied in this method. The results of numerical experiments show that the proposed 2.5D numerical simulation method can be applied to geological models with different boundary conditions (acoustic free-surface, solid free-surface, and solid-liquid boundary). In addition, unlike the 2D wavefield numerical simulation method, the 2.5D method can be directly employed to process actual point source observation data such as 2.5D reverse-time migration.
Keywords: wavefield simulation    generalized first-order wave equation    2.5D    finite-difference method    curvilinear-grid
0 引言

1 时间域2.5D广义一阶速度—应力波动方程

1.1 弹性各向异性介质波动方程

 $\rho \ddot{u}_{j}=\sigma_{i j, i}+F_{j}$ (1)

 $\mathit{\boldsymbol{\sigma }} = \mathit{\boldsymbol{cu}}$ (2)

 $\left\{\begin{array}{l} \dot{\widetilde{v}}_{x}=\frac{\tilde{\sigma}_{x x, x}+\tilde{\sigma}_{x z, z}+\mathrm{i} k_{y} \tilde{\sigma}_{x y}+\tilde{s}_{x}}{\rho} \\ \dot{\tilde{v}}_{y}=\frac{\tilde{\sigma}_{x y, x}+\tilde{\sigma}_{y z, z}+\mathrm{i} k_{y} \tilde{\sigma}_{y y}+\tilde{s}_{y}}{\rho} \\ \dot{\widetilde{v}}_{z}=\frac{{\tilde{\sigma}_{x z, x}}{+\tilde{\sigma}_{z z, z}}+\mathrm{i} k_{y} \tilde{\sigma}_{y z}+\tilde{s}_{z}} {\rho} \end{array}\right.$ (3)

 \left\{\begin{array}{l} \begin{aligned} \dot{\tilde{\sigma}}_{x x}=&c_{11} \tilde{v}_{x, x}+c_{15} \widetilde{v}_{x, z}+c_{16} \tilde{v}_{y, x}+c_{14} \widetilde{v}_{y, z}+c_{15} \tilde{v}_{z, x}+ \\ & c_{13} \tilde{v}_{z, z}+\mathrm{i} k_{y}\left(c_{16} \tilde{v}_{x}+c_{12} \widetilde{v}_{y}+c_{14} \tilde{v}_{z}\right) \end{aligned}\\ \begin{aligned} \dot{\tilde{\sigma}}_{x y}=&c_{16} \widetilde{v}_{x, x}+c_{56} \widetilde{v}_{x, z}+c_{66} \widetilde{v}_{y, x}+c_{46} \widetilde{v}_{y, z}+c_{56} \widetilde{v}_{z, x}+ \\ &c_{36} \widetilde{v}_{z, z}+\mathrm{i} k_{y}\left(c_{66} \widetilde{v}_{x}+c_{26} \widetilde{v}_{y}+c_{46} \widetilde{v}_{z}\right) \end{aligned}\\ \begin{aligned} \dot{\tilde{\sigma}}_{y y}=&c_{12} \tilde{v}_{x, x}+c_{25} \tilde{v}_{x, z}+c_{26} \tilde{v}_{y, x}+c_{24} \widetilde{v}_{y, z}+c_{25} \tilde{v}_{z, x}+ \\ &c_{23} \tilde{v}_{z, z}+\mathrm{i} k_{y}\left(c_{26} \tilde{v}_{x}+c_{22} \tilde{v}_{y}+c_{24} \tilde{v}_{z}\right) \end{aligned}\\ \begin{aligned} \dot{\widetilde{\sigma}}_{x z}=&c_{15} \widetilde{v}_{x, x}+c_{55} \widetilde{v}_{x, z}+c_{56} \widetilde{v}_{y, x}+c_{45} \widetilde{v}_{y, z}+c_{55} \widetilde{v}_{z, x}+\\ &c_{35} \widetilde{v}_{z, z}+\mathrm{i} k_{y}\left(c_{56} \widetilde{v}_{x}+c_{25} \widetilde{v}_{y}+c_{45} \widetilde{v}_{z}\right)\\ \dot{\widetilde{\sigma}}_{y z}=&c_{14} \widetilde{v}_{x, x}+c_{45} \widetilde{v}_{x, z}+c_{46} \widetilde{v}_{y, x}+c_{44} \widetilde{v}_{y, z}+c_{45} \widetilde{v}_{z, x}+\\ &c_{34} \tilde{v}_{z, z}+\mathrm{i} k_{y}\left(c_{46} \tilde{v}_{x}+c_{24} \tilde{v}_{y}+c_{44} \widetilde{v}_{z}\right)\\ \dot{\widetilde{\sigma}}_{z z}=&c_{13} \widetilde{v}_{x, x}+c_{35} \widetilde{v}_{x, z}+c_{36} \widetilde{v}_{y, x}+c_{34} \widetilde{v}_{y, z}+c_{35} \widetilde{v}_{z, x}+\\ &c_{33} \tilde{v}_{z, z}+\mathrm{i} k_{y}\left(c_{36} \widetilde{v}_{x}+c_{23} \tilde{v}_{y}+c_{34} \widetilde{v}_{z}\right) \end{aligned} \end{array}\right. (4)

 $\left\{\begin{array}{l} \tilde{\boldsymbol{v}}=\left(\widetilde{v}_{x}, \widetilde{v}_{y}, \widetilde{v}_{z}\right)^{\mathrm{T}} \\ \tilde{\boldsymbol{s}}=\frac{1}{\rho}\left(\tilde{s}_{x}, \tilde{s}_{y}, \tilde{s}_{z}\right)^{\mathrm{T}} \\ \tilde{\boldsymbol{\sigma}}_{1}=\left(\tilde{\sigma}_{x x}, \tilde{\sigma}_{x y}, \tilde{\sigma}_{y y}\right)^{\mathrm{T}} \\ \tilde{\boldsymbol{\sigma}}_{2}=\left(\tilde{\sigma}_{x z}, \tilde{\sigma}_{y z}, \tilde{\sigma}_{z z}\right)^{\mathrm{T}} \\ \tilde{\boldsymbol{Y}}^{(\mathrm{S})}=\left(\tilde{\boldsymbol{v}}^{\mathrm{T}}, \tilde{\boldsymbol{\sigma}}_{1}^{\mathrm{T}}, \tilde{\boldsymbol{\sigma}}_{2}^{\mathrm{T}}\right)^{\mathrm{T}} \\ \tilde{\boldsymbol{s}}{}^{(\mathrm{S})}=(\tilde{\boldsymbol{s}}, 0,0)^{\mathrm{T}} \end{array}\right.$ (5)

 $\dot{\tilde{\boldsymbol{Y}}}{}^{(\mathrm{S})}=\boldsymbol{A}^{(\mathrm{S})} \partial_{x} \tilde{\boldsymbol{Y}}{}^{(\mathrm{S})}+\boldsymbol{B}^{(\mathrm{S})} \partial_{z} \tilde{\boldsymbol{Y}}{}^{(\mathrm{S})}+\boldsymbol{C}^{(\mathrm{S})} \tilde{\boldsymbol{Y}}{}^{(\mathrm{S})}+\tilde{\boldsymbol{s}}^{(\mathrm{S})}$ (6)

1.2 声波介质波动方程

 $\sigma_{i j}=-\delta_{i j} P$ (7)

P可表示为

 $P=-K u_{j, j}+s(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right)$ (8)

 $\dot{v}_{j}=-\frac{P_{, j}}{\rho}$ (9)

 $\dot{P}=-K v_{j, j}+\dot{s}(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right)$ (10)

 $\dot{P}=-K\left(\partial_{x} v_{x}+\partial_{y} v_{y}+\partial_{z} v_{z}\right)+\dot{s}(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right)$ (11)

 $\left\{\begin{array}{l} \dot{\tilde{v}}_{x}=-\frac{\partial_{x} \widetilde{P}}{\rho} \\ \dot{\tilde{v}}_{y}=-\frac{\mathrm{i} k_{y} \widetilde{P}}{\rho} \\ \dot{\tilde{v}}_{z}=-\frac{\partial_{z} \widetilde{P}}{\rho} \end{array}\right.$ (12)
 $\dot{\widetilde{P}}=-K\left(\partial_{x} {\widetilde{v}}_{x}+\mathrm{i} k_{y} {\tilde{v}}_{y}+\partial_{z} \widetilde{v}_{z}\right)+\dot{\widetilde{s}}(t) \delta\left(x-x_{\mathrm{s}}\right)$ (13)

 $\left\{\begin{array}{l} \tilde{\boldsymbol{Y}}{}^{(\mathrm{W})}=\left(\widetilde{v}_{x}, \widetilde{v}_{y}, \tilde{v}_{z}, \widetilde{P}\right)^{\mathrm{T}} \\ \tilde{\boldsymbol{s}}{}^{(\mathrm{W})}=\left[0,0,0, \dot{\widetilde{s}}(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right)\right]^{\mathrm{T}} \end{array}\right.$ (14)

 $\dot{\tilde{\boldsymbol{Y}}}{}^{(\mathrm{W})}=\boldsymbol{A}^{(\mathrm{W})} \partial_{x} \tilde{\boldsymbol{Y}}{}^{(\mathrm{W})}+\boldsymbol{B}^{(\mathrm{W})} \partial_{z} \tilde{\boldsymbol{Y}}{}^{(\mathrm{W})}+\boldsymbol{C}^{(\mathrm{W})} \tilde{\boldsymbol{Y}}{}^{(\mathrm{W})}+\tilde{\boldsymbol{s}}{}^{(\mathrm{W})}$ (15)

1.3 声波自由地表

 $\left\{\begin{array}{l} \dot{\tilde{v}}_{x}=-\frac{\partial_{x} \widetilde{P}}{\rho} \\ \dot{\tilde{v}}_{y}=0 \\ \dot{\tilde{v}}_{z}=-\frac{-\partial_{z} \widetilde{P}}{\rho} \\ \dot{\tilde{P}}=0 \end{array}\right.$ (16)

{\mathit{\boldsymbol{\tilde Y}}^{({\rm{AW}})}}{\mathit{\boldsymbol{\tilde s}}^{({\rm{AW}})}}定义与式(14)相同，则可以将上式写成矩阵形式  \begin{aligned} \dot{\tilde{\boldsymbol{Y}}}{}^{(\mathrm{AW})}=& \boldsymbol{A}^{(\mathrm{AW})} \partial_{x} \tilde{\boldsymbol{Y}}{}^{(\mathrm{AW})}+\boldsymbol{B}^{(\mathrm{AW})} \partial_{z} \tilde{\boldsymbol{Y}}{}^{(\mathrm{AW})}+\\ & \boldsymbol{C}^{(\mathrm{AW})} \tilde{\boldsymbol{Y}}{}^{(\mathrm{AW})}+\tilde{\boldsymbol{s}}{}^{(\mathrm{AW})} \end{aligned} (17) 式中：上标“(AW)”表示声波自由地表；系数矩阵A(AW)B(AW)C(AW)分别为  \boldsymbol{A}^{(\mathrm{AW})}=-\rho^{-1}\left(\begin{array}{llll} 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array}\right) (18)  \boldsymbol{B}^{(\text {AW })}=-\rho^{-1}\left(\begin{array}{llll} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{array}\right) (19)  \boldsymbol{C}^{(\text {AW })}=\bf{0} (20) 对比式(18)~式(20)与附录B，可以发现式(17)是式(15)的特例。 1.4 固体自由地表 固体自由地表边界条件需要满足地表处法向应力为0，即  \boldsymbol{\sigma} \cdot \boldsymbol{n}=0 (21) 展开为  \left\{\begin{array}{l} \sigma_{x x} n_{1}+\sigma_{x y} n_{2}+\sigma_{x z} n_{3}=0 \\ \sigma_{x y} n_{1}+\sigma_{y y} n_{2}+\sigma_{y z} n_{3}=0 \\ \sigma_{x z} n_{1}+\sigma_{y z} n_{2}+\sigma_{z z} n_{3}=0 \end{array}\right. (22) 式中n =(n1n2n3)表示界面的法向向量。将式(22)写成矩阵形式  \left(\begin{array}{cccccc} n_{1} & n_{2} & 0 & n_{3} & 0 & 0 \\ 0 & n_{1} & n_{2} & 0 & n_{3} & 0 \\ 0 & 0 & 0 & n_{1} & n_{2} & n_{3} \end{array}\right)\left(\begin{array}{l} \sigma_{x x} \\ \sigma_{x y} \\ \sigma_{y y} \\ \sigma_{x z} \\ \sigma_{y z} \\ \sigma_{z z} \end{array}\right)=0 (23) 将上式写为子矩阵格式  \boldsymbol{\varGamma}_{2} \boldsymbol{\sigma}_{2}=-\boldsymbol{\varGamma}_{1} \boldsymbol{\sigma}_{1} (24) 式中  \left\{\begin{array}{l} \boldsymbol{\sigma}_{1}=\left(\sigma_{x x}, \sigma_{x y}, \sigma_{y y}\right) \\ \boldsymbol{\sigma}_{2}=\left(\sigma_{x z}, \sigma_{y z}, \sigma_{z z}\right) \end{array}\right. (25)  \left\{\begin{array}{l} \boldsymbol{\varGamma}_{1}=\left(\begin{array}{lll} n_{1} & n_{2} & 0 \\ 0 & n_{1} & n_{2} \\ 0 & 0 & 0 \end{array}\right) \\ \boldsymbol{\varGamma}_{2}=\left(\begin{array}{ccc} n_{3} & 0 & 0 \\ 0 & n_{3} & 0 \\ n_{1} & n_{2} & n_{3} \end{array}\right) \end{array}\right. (26) 定义  \begin{aligned} \boldsymbol{S} &=-\boldsymbol{\varGamma}_{2}^{-1} \boldsymbol{\varGamma}_{1} \\ &=\left(\begin{array}{ccc} n_{3}^{-1} & 0 & 0 \\ 0 & n_{3}^{-1} & 0 \\ -n_{1} n_{3}^{-2} & -n_{2} n_{3}^{-2} & n_{3}^{-1} \end{array}\right)\left(\begin{array}{ccc} n_{1} & n_{2} & 0 \\ 0 & n_{1} & n_{2} \\ 0 & 0 & 0 \end{array}\right) \\ &=\left(\begin{array}{ccc} -n_{1} n_{3}^{-1} & -n_{2} n_{3}^{-1} & 0 \\ 0 & -n_{1} n_{3}^{-1} & -n_{2} n_{3}^{-1} \\ n_{1}^{2} n_{3}^{-2} & 2 n_{1} n_{2} n_{3}^{-2} & n_{2}^{2} n_{3}^{-2} \end{array}\right) \end{aligned} (27) 式(23)可改写为  \boldsymbol{\sigma}_{2}=\boldsymbol{S} \boldsymbol{\sigma}_{1} (28) 由上式可知，应力分量σ2可由σ1和表面法向分量n得到，只有应力分量σ1未知。假设自由地表由z=z0(x)给出，其坡度为tanα=z0(x)，法向分量为n =(－sinα, 0, cosα)，则式(27)变为  \boldsymbol{S}=\left(\begin{array}{ccc} z_{0}^{\prime}(x) & 0 & 0 \\ 0 & z_{0}^{\prime}(x) & 0 \\ {\left[z_{0}^{\prime}(x)\right]^{2}} & 0 & 0 \end{array}\right) (29) 根据式(28)和式(29)可得  \left\{\begin{array}{l} \dot{\boldsymbol{\sigma}}_{2}=\boldsymbol{S}\dot{\boldsymbol{\sigma}}_{1} \\ \partial_{x} \boldsymbol{\sigma}_{2}=\partial_{x} \boldsymbol{S} \boldsymbol{\sigma}_{1}+\boldsymbol{S} \partial_{x} \boldsymbol{\sigma}_{1} \\ \partial_{z} \boldsymbol{\sigma}_{2}=\boldsymbol{S} \partial_{z} \boldsymbol{\sigma}_{1} \end{array}\right. (30) 将上式代入式(3)，可得  \left\{\begin{aligned} \dot{\tilde{v}}_{x}=&\left[\tilde{\sigma}_{x x, x}+z_{0}^{\prime}(x) \tilde{\sigma}_{x x, z}+\mathrm{i} k_{y} \tilde{\sigma}_{x y}+\tilde{s}_{x}\right] / \rho \\ \dot{\tilde{v}}_{y}=&\left[\tilde{\sigma}_{x y, x}+z_{0}^{\prime}(x) \tilde{\sigma}_{x y, z}+\mathrm{i} k_{y} \tilde{\sigma}_{y y}+\tilde{s}_{y}\right] / \rho \\ \dot{\tilde{v}}_{z}=&\left\{z_{0}^{\prime}(x) \tilde{\sigma}_{x x, x}+\left[z_{0}^{\prime}(x)\right]^{2} \tilde{\sigma}_{x x, z}+\right.\\ &\left.\mathrm{i} k_{y} z_{0}^{\prime}(x) \tilde{\sigma}_{x y}+z_{0}^{\prime \prime}(x) \tilde{\sigma}_{x x}+\tilde{s}_{z}\right\} / \rho \end{aligned}\right. (31) 由式(4)可得  \left\{\begin{array}{l} \dot{\widetilde{\sigma}}_{x x}=c_{11} \tilde{v}_{x, x}+c_{15} \tilde{v}_{x, z}+c_{16} \widetilde{v}_{y, x}+c_{14} \widetilde{v}_{y, z}+c_{15} \widetilde{v}_{z, x}+c_{13} \widetilde{v}_{z, z}+\mathrm{i} k_{y}\left(c_{16} \widetilde{v}_{x}+c_{12} \widetilde{v}_{y}+c_{14} \widetilde{v}_{z}\right) \\ \dot{\widetilde{\sigma}}_{x y}=c_{16} \widetilde{v}_{x, x}+c_{56} \widetilde{v}_{x, z}+c_{66} \widetilde{v}_{y, x}+c_{46} \widetilde{v}_{y, z}+c_{56} \widetilde{v}_{z, x}+c_{36} \widetilde{v}_{z, z}+\mathrm{i} k_{y}\left(c_{66} \widetilde{v}_{x}+c_{26} \widetilde{v}_{y}+c_{46} \widetilde{v}_{z}\right) \\ \dot{\widetilde{\sigma}}_{y y}=c_{12} \widetilde{v}_{x, x}+c_{25} \widetilde{v}_{x, z}+c_{26} \widetilde{v}_{y, x}+c_{24} \widetilde{v}_{y, z}+c_{25} \widetilde{v}_{z, x}+c_{23} \widetilde{v}_{z, z}+\mathrm{i} k_{y}\left(c_{26} \widetilde{v}_{x}+c_{22} \widetilde{v}_{y}+c_{24} \widetilde{v}_{z}\right) \end{array}\right. (32) 定义  \left\{\begin{array}{l} \tilde{\boldsymbol{Y}}{}^{(\mathrm{AS})}=\left(\widetilde{v}_{x}, \widetilde{v}_{y}, \widetilde{v}_{z}, \tilde{\sigma}_{x x}, \tilde{\sigma}_{x y}, \tilde{\sigma}_{y y}\right)^{\mathrm{T}} \\ \tilde{\boldsymbol{s}}{}^{(\mathrm{AS})}=\rho^{-1}\left(\tilde{s}_{x}, \tilde{s}_{y}, \tilde{s}_{z}, 0,0,0\right)^{\mathrm{T}} \end{array}\right. (33) 则式(30)和式(31)可简写为  \begin{gathered} \dot{\tilde{\boldsymbol{Y}}}{}^{(\mathrm{AS})}=\boldsymbol{A}^{(\mathrm{AS})} \partial_{x} \tilde{\boldsymbol{Y}}{}^{(\mathrm{AS})}+\boldsymbol{B}^{(\mathrm{AS})} \partial_{z} \tilde{\boldsymbol{Y}}{}^{(\mathrm{AS})}+ \\ \boldsymbol{C}^{(\mathrm{AS})} \tilde{\boldsymbol{Y}}{}^{(\mathrm{AS})}+\tilde{\boldsymbol{s}}{}^{(\mathrm{AS})} \end{gathered} (34) 式中：上标“(AS)”表示固体自由地表边界；系数矩阵A(W)B(W)C(W)的具体表达式见附录C。式(34)为固体自由地表处的2.5D地震波控制方程。与式(5)、式(6)相比，式(33)、式(34) 根据固体自由地表边界条件消除了应力分量{\mathit{\boldsymbol{\tilde \sigma }}_2}，只需求解应力分量 {\mathit{\boldsymbol{\tilde \sigma }}_1}和速度向量\mathit{\boldsymbol{\tilde v }} 1.5 固—液边界 在固—液边界上，例如由z=z0(x)定义的二维海底界面，边界条件需要满足法向应力和法向速度分量连续，并且切向应力分量为0，即  \tilde{\boldsymbol{\sigma}}\left[t, x, z_{0}(x)\right] \cdot \boldsymbol{n}=-\widetilde{P}\left[t, x, z_{0}(x)\right] \cdot \boldsymbol{n} (35)  \tilde{\boldsymbol{v}}{}^{(\mathrm{S})}\left[t, x, z_{0}(x)\right] \cdot \boldsymbol{n}=\tilde{\boldsymbol{v}}{}^{(\mathrm{W})}\left[t, x, z_{0}(x)\right] \cdot \boldsymbol{n} (36)  \left\{\begin{array}{l} \tilde{\boldsymbol{\sigma}}\left[t, x, z_{0}(x)\right] \cdot \boldsymbol{\tau}^{(1)}=0 \\ \tilde{\boldsymbol{\sigma}}\left[t, x, z_{0}(x)\right] \cdot \boldsymbol{\tau}^{(2)}=0 \end{array}\right. (37) 式中：\mathit{\boldsymbol{\tilde \sigma }}是固体介质中的应力；{\mathit{\boldsymbol{\tilde v}}^{({\rm{S}})}}{\mathit{\boldsymbol{\tilde v}}^{(W)}}分别为界面两侧固体和液体介质的振动速度矢量。根据界面处的坡度tanα=z0(x)，界面处的切向矢量为

 $\left\{\begin{array}{l} \boldsymbol{\tau}^{(1)}=(0,1,0)^{\mathrm{T}} \\ \boldsymbol{\tau}^{(2)}=(\cos \alpha, 0, \sin \alpha)^{\mathrm{T}} \end{array}\right.$ (38)

 $\left\{\begin{array}{l} \tilde{\sigma}_{12}=\tilde{\sigma}_{22}=\tilde{\sigma}_{23}=0 \\ \tilde{\sigma}_{11}=-z^{\prime}{}_{0}(x) \tilde{\sigma}_{13} \\ \tilde{\sigma}_{13}=-z^{\prime}{}_{0}(x) \tilde{\sigma}_{33} \end{array}\right.$ (39)

 $\tilde{\sigma}_{33}=\frac{-1}{1+\left[z_{0}^{\prime}(x)\right]^{2}} \widetilde{P}$ (40)

 $\dot{\tilde{P}}=-\left\{1+\left[z_{0}^{\prime}(x)\right]^{2}\right\} \dot{\widetilde{\sigma}}_{33}$ (41)

 \begin{aligned} &\partial_{x} \tilde{v}_{x}^{(\mathrm{W})}+\mathrm{i} k_{y} \tilde{v}_{y}^{(\mathrm{W})}+\partial_{z} \tilde{v}_{z}^{(\mathrm{W})} \\ &=\frac{1+\left[z_{0}^{\prime}(x)\right]^{2}}{K}\left[\left(c_{13}, c_{36}, c_{35}\right) \cdot \partial_{x} \tilde{\boldsymbol{v}}^{(\mathrm{S})}+\right. \\ &\left.\left(c_{35}, c_{34}, c_{33}\right) \cdot \partial_{x} \tilde{\boldsymbol{v}}^{(\mathrm{S})}+\mathrm{i} k_{y}\left(c_{36}, c_{23}, c_{34}\right) \cdot \tilde{\boldsymbol{v}}^{(\mathrm{S})}\right] \end{aligned} (42)

 $\tilde{v}_{y}^{(\mathrm{W})}=\frac{1+\left[z_{0}^{\prime}(x)\right]^{2}}{K}\left[c_{36} \widetilde{v}_{x}^{(\mathrm{S})}+c_{23} \widetilde{v}_{y}^{(\mathrm{S})}+c_{34} \widetilde{v}_{z}^{(\mathrm{S})}\right]$ (43)

 $\widetilde{v}_{z}^{(\mathrm{W})}=z_{0}^{\prime}(x)\left[\widetilde{v}_{x}^{(\mathrm{W})}-\widetilde{v}_{x}^{(\mathrm{S})}\right]+\widetilde{v}_{z}^{(\mathrm{S})}$ (44)

 \begin{aligned} &\dot{\tilde{P}}=-\left\{1+\left[z_{0}^{\prime}(x)\right]^{2}\right\}\left[\left(c_{13}, c_{36}, c_{35}\right) \cdot \partial_{x} \tilde{\boldsymbol{v}}^{(\mathrm{S})}+\right. \\ &\left.\left(c_{35}, c_{34}, c_{33}\right) \cdot \partial_{z} \tilde{\boldsymbol{v}}^{(\mathrm{S})}+\mathrm{i} k_{y}\left(c_{36}, c_{23}, c_{34}\right) \cdot \tilde{\boldsymbol{v}}^{(\mathrm{S})}\right] \end{aligned} (45)

 $\left\{\begin{array}{l} \dot{\widetilde{v}}_{x}^{(\mathrm{S})}=\rho_{\mathrm{S}}^{-1} \frac{z_{0}^{\prime}(x)}{\left[z_{0}^{\prime}(x)\right]^{2}+1}\left\{-z_{0}^{\prime}(x) \partial_{x} \widetilde{P}+\right. \\ \ \ \ \ \ \ \ \ \ \ \ \left.\partial_{z} \widetilde{P}+\frac{-2 z_{0}^{\prime \prime}(x)}{\left[z_{0}^{\prime}(x)\right]^{2}+1} \widetilde{P}\right\} \\ \dot{v}_{y}^{(\mathrm{S})}=\rho_{\mathrm{S}}^{-1}\left[-z_{0}^{\prime}(x) \partial_{z} \tilde{\sigma}_{12}+\partial_{z} \tilde{\sigma}_{23}\right] \\ \dot{v}_{z}^{(\mathrm{S})}=\rho_{\mathrm{S}}^{-1} \frac{1}{\left[z_{0}^{\prime}(x)\right]^{2}+1}\left\{-z_{0}^{\prime}(x) \partial_{x} \widetilde{P}-\right. \\ \ \ \ \ \ \ \ \ \ \ \ \left.\partial_{z} \widetilde{P}+\frac{-z_{0}^{\prime \prime}(x)\left[1-\left(z_{0}^{\prime}(x)\right)^{2}\right]}{\left[z_{0}^{\prime}(x)\right]^{2}+1}\widetilde{P}\right\} \end{array}\right.$ (46)

 $\dot{\tilde{v}}{}_{x}^{(\mathrm{W})}=\rho_{\mathrm{W}}^{-1}\left(-\partial_{x} \widetilde{P}\right)$ (47)

 $\left\{\begin{array}{l} \tilde{\boldsymbol{Y}}{}^{(\mathrm{WS})}=\left[\tilde{\boldsymbol{v}}^{(\mathrm{WS})}, \tilde{\mathrm{P}}, \tilde{\sigma}_{12}, \tilde{\sigma}_{23}\right]^{\mathrm{T}} \\ \tilde{\boldsymbol{v}}{}^{(\mathrm{WS})}=\left[\widetilde{v}_{x}^{(\mathrm{W})}, \widetilde{v}_{x}^{(\mathrm{S})}, \widetilde{v}_{y}^{(\mathrm{S})}, \widetilde{v}_{z}^{(\mathrm{S})}\right] \\ \tilde{\boldsymbol{s}}{}^{(\mathrm{WS})}=\left\{\begin{array}{l} {\left[0,0,0,0, \dot{\widetilde{s}}(t) \delta\left(\boldsymbol{x-x}_{s}\right), 0,0\right]^{\mathrm{T}}(\text { 液体中 })} \\ \left(0, \tilde{s}_{x}, \tilde{s}_{y}, \tilde{s}_{z}, 0,0,0\right)^{\mathrm{T}} \text { (固体中) } \end{array}\right. \end{array}\right.$ (48)

 $\begin{gathered} \dot{\tilde{\boldsymbol{Y}}}{}^{(\mathrm{WS})}=\boldsymbol{A}^{(\mathrm{WS})} \partial_{x} \tilde{\boldsymbol{Y}}{}^{(\mathrm{WS})}+\boldsymbol{B}^{(\mathrm{WS})} \partial_{z} \tilde{\boldsymbol{Y}}{}^{(\mathrm{WS})}+ \\ \boldsymbol{C}^{(\mathrm{WS})} \tilde{\boldsymbol{Y}}{}^{(\mathrm{WS})}+\tilde{\boldsymbol{s}}{}^{(\mathrm{WS})} \end{gathered}$ (49)

 $\dot{\tilde{\boldsymbol{Y}}}=\boldsymbol{A} \partial_{x} \tilde{\boldsymbol{Y}}+\boldsymbol{B} \partial_{z} \tilde{\boldsymbol{Y}}+{\boldsymbol{C}}\tilde{\boldsymbol{Y}}+\tilde{\boldsymbol{s}}$ (50)

2 曲线网格有限差分算法 2.1 曲线坐标系与直角坐标系的转换关系

 图 1 曲线网格和计算网格的变换关系示意图
 $\left\{\begin{array}{l} x=x(\xi, \eta) \\ z=z(\xi, \eta) \end{array}\right.$ (51)

 $\left\{\begin{array}{l} 1=x_{, \xi} \xi_{, x}+x_{, \eta} \eta_{, x} \\ 0=x_{, \xi} \xi_{, z}+x_{, \eta} \eta_{, z} \\ 0=z_{, \xi} \xi_{, x}+z_{, \eta} \eta_{, x} \\ 1=z_{, \xi} \xi_{, z}+z_{, \eta} \eta_{, z} \end{array}\right.$ (52)

 $\left\{\begin{array}{l} \xi_{, x}=\frac{1}{J} z_{, \eta} \\ \xi_{, z}=-\frac{1}{J} x_{, \eta} \\ \eta_{, x}=-\frac{1}{J} z_{, \xi} \\ \eta_{, z}=\frac{1}{J} x_{, \xi} \end{array}\right.$ (53)

2.2 曲线坐标系中的广义一阶速度—应力波动方程

 \begin{aligned} \dot{\tilde{\boldsymbol{Y}}}=&\left(\partial_{x} \xi \boldsymbol{A}+\partial_{z} \xi \boldsymbol{B}\right) \partial_{\xi} \tilde{\boldsymbol{Y}}+\\ &\left(\partial_{x} \eta \boldsymbol{A}+\partial_{z} \eta \boldsymbol{B}\right) \partial_{\eta} \widetilde{\boldsymbol{Y}}+\boldsymbol{C} \widetilde{\boldsymbol{Y}}+\widetilde{\boldsymbol{s}} \end{aligned} (54)

Hixon等[24]采用Tam等[25]提出的DRP (Dispersion Relation Preserving)方法对同位MacCormack格式的系数进行了优化，得到了低频散、低耗散的DRP/opt MacCormack格式，本文采用该格式有限差分算子离散方程。

 $L_{\xi}^{\mathrm{F}}(\tilde{\boldsymbol{Y}})_{i}=\frac{1}{\Delta \xi} \sum\limits_{j=-1}^{3} a_{j} \tilde{\boldsymbol{Y}}_{i+j}$ (55)
 $L_{\mathtt{ξ}}^{\mathrm{B}}(\tilde{\boldsymbol{Y}})_{i}=\frac{1}{\Delta \xi} \sum\limits_{j=-1}^{3}-a_{j} \tilde{\boldsymbol{Y}}_{i-j}$ (56)

 $\left\{\begin{array}{l} a_{-1}=-0.30874 \\ a_{0}=-0.63260 \\ a_{1}=1.23300 \\ a_{2}=-0.33340 \\ a_{3}=0.04168 \end{array}\right.$ (57)

 $\left\{\begin{array}{l} \hat{L}{}^{\mathrm{FF}}(\tilde{\boldsymbol{Y}})=\left(\partial_{x} \xi \boldsymbol{A}+\partial_{z} \xi \boldsymbol{B}\right) L_{\xi}^{\mathrm{F}}(\tilde{\boldsymbol{Y}})+\left(\partial_{x} \eta \boldsymbol{A}+\partial_{z} \eta \boldsymbol{B}\right) L_{\eta}^{\mathrm{F}}(\tilde{\boldsymbol{Y}}) \\ \hat{L}{}^{\mathrm{BB}}(\tilde{\boldsymbol{Y}})=\left(\partial_{x} \xi \boldsymbol{A}+\partial_{z} \xi \boldsymbol{B}\right) L_{\xi}^{\mathrm{B}}(\tilde{\boldsymbol{Y}})+\left(\partial_{x} \eta \boldsymbol{A}+\partial_{z} \eta \boldsymbol{B}\right) L_{\eta}^{\mathrm{B}}(\tilde{\boldsymbol{Y}}) \\ \hat{L}{}^{\mathrm{FB}}(\tilde{\boldsymbol{Y}})=\left(\partial_{x} \xi \boldsymbol{A}+\partial_{z} \xi \boldsymbol{B}\right) L_{\xi}^{\mathrm{F}}(\tilde{\boldsymbol{Y}})+\left(\partial_{x} \eta \boldsymbol{A}+\partial_{z} \eta \boldsymbol{B}\right) L_{\eta}^{\mathrm{B}}(\tilde{\boldsymbol{Y}}) \\ \hat{L}{}^{\mathrm{BF}}(\tilde{\boldsymbol{Y}})=\left(\partial_{x} \xi \boldsymbol{A}+\partial_{z} \xi \boldsymbol{B}\right) L_{\xi}^{\mathrm{B}}(\tilde{\boldsymbol{Y}})+\left(\partial_{x} \eta \bf{A}+\partial_{z} \eta \boldsymbol{B}\right) L_{\eta}^{\mathrm{F}}(\tilde{\boldsymbol{Y}}) \end{array}\right.$ (58)

 $\left\{\begin{array}{l} \boldsymbol{h}^{(1)}=\Delta t \hat{L}{}^{\mathrm{FF}}\left(\tilde{\boldsymbol{Y}}{}^{N}\right) \\ \boldsymbol{h}^{(2)}=\Delta t \hat{L}{}^{\mathrm{BB}}\left[\tilde{\boldsymbol{Y}}{}^{N}+\alpha_{1} \boldsymbol{h}^{(1)}\right] \\ \boldsymbol{h}^{(3)}=\Delta t \hat{L}{}^{\mathrm{FF}}\left[\tilde{\boldsymbol{Y}}{}^{N}+\alpha_{2} \boldsymbol{h}^{(2)}\right] \\ \boldsymbol{h}^{(4)}=\Delta t \hat{L}{}^{\mathrm{BB}}\left[\tilde{\boldsymbol{Y}}{}^{N}+\alpha_{3} \boldsymbol{h}^{(3)}\right] \\ \tilde{\boldsymbol{Y}}{}^{N+1}=\tilde{\boldsymbol{Y}}{}^{N}+\beta_{1} \boldsymbol{h}^{(1)}+\beta_{2} \boldsymbol{h}^{(2)}+\beta_{3} \boldsymbol{h}^{(3)}+\beta_{4} \boldsymbol{h}^{(4)} \end{array}\right.$ (59)

 $\left\{\begin{array}{l} \boldsymbol{h}^{(1)}=\Delta t \hat{L}{}^{\mathrm{FB}}\left(\tilde{\boldsymbol{Y}}{}^{N+1}\right) \\ \boldsymbol{h}^{(2)}=\Delta t \hat{L}{}^{\mathrm{BF}}\left[\tilde{\boldsymbol{Y}}{}^{N+1}+\alpha_{1} \boldsymbol{h}^{(1)}\right] \\ \boldsymbol{h}^{(3)}=\Delta t \hat{L}{}^{\mathrm{FB}}\left[\tilde{\boldsymbol{Y}}{}^{N+1}+\alpha_{2} \boldsymbol{h}^{(2)}\right] \\ \boldsymbol{h}^{(4)}=\Delta t \hat{L}{}^{\mathrm{BF}}\left[\tilde{\boldsymbol{Y}}{}^{N+1}+\alpha_{3} \boldsymbol{h}^{(3)}\right] \\ \tilde{\boldsymbol{Y}}{}^{N+2}=\tilde{\boldsymbol{Y}}{}^{N+1}+\beta_{1} \boldsymbol{h}^{(1)}+\beta_{2} \boldsymbol{h}^{(2)}+\beta_{3} \boldsymbol{h}^{(3)}+\beta_{4} \boldsymbol{h}^{(4)} \end{array}\right.$ (60)

 $\left\{\begin{array}{l} \boldsymbol{h}^{(1)}=\Delta t \hat{L}{}^{\mathrm{BB}}\left(\tilde{\boldsymbol{Y}}{}^{N+2}\right) \\ \boldsymbol{h}^{(2)}=\Delta t \hat{L}{}^{\mathrm{FF}}\left[\tilde{\boldsymbol{Y}}{}^{N+2}+\alpha_{1} \boldsymbol{h}^{(1)}\right] \\ \boldsymbol{h}^{(3)}=\Delta t \hat{L}{}^{\mathrm{BB}}\left[\tilde{\boldsymbol{Y}}{}^{N+2}+\alpha_{2} \boldsymbol{h}^{(2)}\right] \\ \boldsymbol{h}^{(4)}=\Delta t \hat{L}{}^{\mathrm{FF}}\left[\tilde{\boldsymbol{Y}}{}^{N+2}+\alpha_{3} \boldsymbol{h}^{(3)}\right] \\ \tilde{\boldsymbol{Y}}{}^{N+3}=\tilde{\boldsymbol{Y}}{}^{N+2}+\beta_{1} \boldsymbol{h}^{(1)}+\beta_{2} \boldsymbol{h}^{(2)}+\beta_{3} \boldsymbol{h}^{(3)}+\beta_{4} \boldsymbol{h}^{(4)} \end{array}\right.$ (61)

 $\left\{\begin{array}{l} \boldsymbol{h}^{(1)}=\Delta t \hat{L}{}^{\mathrm{BF}}\left(\tilde{\boldsymbol{Y}}{}^{{N}+3}\right) \\ \boldsymbol{h}^{(2)}=\Delta t \hat{L}{}^{\mathrm{FB}}\left[\tilde{\boldsymbol{Y}}{}^{{N}+3}+\alpha_{1} \boldsymbol{h}^{(1)}\right] \\ \boldsymbol{h}^{(3)}=\Delta t \hat{L}{}^{\mathrm{BF}}\left[\tilde{\boldsymbol{Y}}{}^{{N}+3}+\alpha_{2} \boldsymbol{h}^{(2)}\right] \\ \boldsymbol{h}^{(4)}=\Delta t \hat{L}{}^{\mathrm{FB}}\left[\tilde{\boldsymbol{Y}}{}^{{N}+3}+\alpha_{3} \boldsymbol{h}^{(3)}\right] \\ \tilde{\boldsymbol{Y}}{}^{{N}+4}=\tilde{\boldsymbol{Y}}{}^{{N}+3}+\beta_{1} \boldsymbol{h}^{(1)}+\beta_{2} \boldsymbol{h}^{(2)}+\beta_{3} \boldsymbol{h}^{(3)}+\beta_{4} \boldsymbol{h}^{(4)} \end{array}\right.$ (62)

3 震源、人工边界和波数域采样 3.1 震源和吸收边界

 $s(t)=\left[1-2 {\rm{ \mathit{ π} }}^{2} f_{\mathrm{c}}^{2}\left(t-t_{0}\right)^{2}\right] \exp \left[-{\rm{ \mathit{ π} }}^{2} f_{\mathrm{c}}^{2}\left(t-t_{0}\right)^{2}\right]$ (63)

3.2 波数域采样方法

 \begin{aligned} \boldsymbol{Y}\left(t_{{N}}, x, y, z\right) \approx & \frac{2 k_{y}^{\mathrm{c}}}{{\rm{ \mathit{ π} }}\left(M_{k y}-1\right)} \sum\limits_{m=1}^{M_{k y}} \tilde{\boldsymbol{Y}}\left(t_{{N}}, x, k_{y}^{m}, z\right) \mathrm{e}^{\mathrm{i} k_{y}^{m} y} \\ & k_{y}^{m} \in\left[-k_{y}^{\mathrm{c}}, k_{y}^{\mathrm{c}}\right] \end{aligned} (64)

 $v_{i j}\left(t_{N}, x, y, z\right)= \begin{cases}\frac{k_{y}^{\mathrm{c}}}{{\rm{ \mathit{ π} }}\left(N_{k y}-1\right)} \sum\limits_{m=1}^{M_{k y}} \tilde{v}_{i j}\left(t_{N}, x, k_{y}^{m}, z\right) \cos \left(k_{y}^{m} y\right) & i+j \text { 为偶数 } \\ \frac{k_{y}^{\mathrm{c}}}{{\rm{ \mathit{ π} }}\left(N_{k y}-1\right)} \sum\limits_{m=1}^{M_{k y}} {\rm{i}}\tilde{v}_{i j}\left(t_{N}, x, k_{y}^{m}, z\right) \sin \left(k_{y}^{m} y\right) & i+j \text { 为奇数 }\end{cases}$ (65)

 $\begin{gathered} P\left(t_{N}, x, y, z\right)=\frac{k_{y}^{\mathrm{c}}}{{\rm{ \mathit{ π} }}\left(M_{k y}-1\right)} \times \\ \sum\limits_{m=1}^{M_{k y}} \widetilde{P}\left(t_{N}, x, k_{y}^{m}, z\right) \cos \left(k_{y}^{m} y\right) \end{gathered}$ (66)
4 2.5D与3D数值模拟的精度和效率分析

 图 2 地表水平均匀介质模型2.5D模拟的波场快照 (a)声波介质的P分量，0.15s；(b)、(c)分别为各向同性介质的vx和vz分量，0.09s；(d)、(e)分别为VTI-1介质的中vx和vz分量，0.08s

 图 3 地表水平均匀介质模型2.5D数值解与3D数值解、3D解析解的对比 (a)声波介质中P分量；(b)、(c)分别为各向同性介质的vx和vz分量；(d)、(e)分别为VTI-1介质的vx和vz分量

2.5D与3D数值模拟计算机内存占用和CPU运行时间的对比如表 2所示，可以看出，2.5D数值模拟所需CPU时间和计算机内存都远小于3D，特别是计算机内存。所以本文的2.5D波场数值模拟方法则具有较高的实际应用价值，可以直接应用于地震数据反演方法。

5 复杂模型的波场数值模拟 5.1 起伏地表均匀模型

 图 4 模型网格剖分示意图 红星和蓝三角分别为震源和检波器

 图 5 起伏地表声波均匀介质模型2.5D数值解vx(a)和vz(b)分量在0.4s的波场快照

 图 6 起伏地表各向同性均匀弹性介质模型2.5D数值解vx(a)和vz(b)分量在0.32s的波场快照

 图 7 起伏地表均匀VTI-1介质模型2.5D数值解vx(a)和vz(b)分量在0.32s的波场快照

 图 8 起伏地表声波均匀介质模型2.5D与2D数值解对比 (a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

 图 9 起伏地表均匀各向同性弹性介质模型2.5D与2D数值解的对比 (a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

 图 10 起伏地表均匀VTI-1介质模型2.5D与2D数值解的对比 (a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示
5.2 固—液耦合模型

 图 11 固—液耦合模型2.5D数值解vx(a)和vz(b)分量在0.28s的波场快照 图中起伏曲线为正弦起伏界面

 图 12 固—液耦合模型2.5D与2D数值解波形对比 (a)vx分量；(b)vz分量
5.3 三层混合模型

 图 13 三层模型示意图 星和三角分别表示震源和接收器点位置

 图 14 三层模型2.5D数值解vx分量在不同时刻的波场快照 (a)0.12s；(b)0.20s；(c)0.28s；(d)0.36s；(e)0.44s；(f)0.52s

 图 15 三层模型第一个检波器的2.5D与2D数值模拟记录对比 (a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

 图 16 三层模型第二个检波器的2.5D与2D数值模拟记录对比 (a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

 图 17 三层模型第三个检波器的2.5D与2D数值模拟记录对比 (a)vx分量；(b)vz分量；(c)图a的归一化显示；(d)图b的归一化显示。下图为上图方框的放大显示

6 结论

 $\boldsymbol{A}^{(\mathrm{S})}=\left(\begin{array}{ccc} 0 & \boldsymbol{U}_{1} & \boldsymbol{U}_{2} \\ \boldsymbol{W}_{1} & 0 & 0 \\ \boldsymbol{W}_{4} & 0 & 0 \end{array}\right)$ (A-1)
 $\boldsymbol{B}^{(\mathrm{S})}=\left(\begin{array}{ccc} 0 & 0 & \boldsymbol{U}_{3} \\ \boldsymbol{W}_{2} & 0 & 0 \\ \boldsymbol{W}_{5} & 0 & 0 \end{array}\right)$ (A-2)
 $\boldsymbol{C}^{(\mathrm{S})}=\left(\begin{array}{ccc} 0 & \boldsymbol{U}_{4} & \boldsymbol{U}_{5} \\ \boldsymbol{W}_{3} & 0 & 0 \\ \boldsymbol{W}_{6} & 0 & 0 \end{array}\right)$ (A-3)

 $\boldsymbol{U}_{1}=\rho^{-1}\left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{array}\right)$ (A-4)
 $\boldsymbol{U}_{2}=\rho^{-1}\left(\begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{array}\right)$ (A-5)
 $\boldsymbol{U}_{3}=\rho^{-1}\left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right)$ (A-6)
 $\boldsymbol{U}_{4}={\rm{i}} k_{y} \rho^{-1}\left(\begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{array}\right)$ (A-7)
 $\boldsymbol{U}_{5}=\mathrm{i} k_{y} \rho^{-1}\left(\begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 1 & 0 \end{array}\right)$ (A-8)
 $\boldsymbol{W}_{1}=\left(\begin{array}{ccc} c_{11} & c_{16} & c_{15} \\ c_{16} & c_{66} & c_{56} \\ c_{12} & c_{26} & c_{25} \end{array}\right)$ (A-9)
 $\boldsymbol{W}_{2}=\left(\begin{array}{lll} c_{15} & c_{14} & c_{13} \\ c_{56} & c_{46} & c_{36} \\ c_{25} & c_{24} & c_{23} \end{array}\right)$ (A-10)
 $\boldsymbol{W}_{3}=\mathrm{i} k_{y}\left(\begin{array}{ccc} c_{16} & c_{12} & c_{14} \\ c_{66} & c_{26} & c_{46} \\ c_{26} & c_{22} & c_{24} \end{array}\right)$ (A-11)
 $\boldsymbol{W}_{4}=\left(\begin{array}{lll} c_{15} & c_{56} & c_{55} \\ c_{14} & c_{46} & c_{45} \\ c_{13} & c_{36} & c_{35} \end{array}\right)$ (A-12)
 $\boldsymbol{W}_{5}=\left(\begin{array}{lll} c_{55} & c_{45} & c_{35} \\ c_{45} & c_{44} & c_{34} \\ c_{35} & c_{34} & c_{33} \end{array}\right)$ (A-13)
 $\boldsymbol{W}_{6}=\mathrm{i} k_{y}\left(\begin{array}{ccc} c_{56} & c_{25} & c_{45} \\ c_{46} & c_{24} & c_{44} \\ c_{36} & c_{23} & c_{34} \end{array}\right)$ (A-14)

 $\boldsymbol{A}^{(\mathrm{W})}=\left(\begin{array}{cccc} 0 & 0 & 0 & -\rho^{-1} \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ -K & 0 & 0 & 0 \end{array}\right)$ (B-1)
 $\boldsymbol{B}^{(\mathrm{W})}=\left(\begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\rho^{-1} \\ 0 & 0 & -K & 0 \end{array}\right)$ (B-2)
 $\boldsymbol{C}^{(\mathrm{W})}=-\mathrm{i} k_{y}\left(\begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \rho^{-1} \\ 0 & 0 & 0 & 0 \\ 0 & K & 0 & 0 \end{array}\right)$ (B-3)

 $\boldsymbol{A}^{(\mathrm{AS})} =\left(\begin{array}{cc} 0 & \boldsymbol{A}_{1}^{(\mathrm{AS})} \\ \boldsymbol{A}_{2}^{(\mathrm{AS})} & 0 \end{array}\right)$ (C-1)
 $\boldsymbol{B}^{(\mathrm{AS})} =\left(\begin{array}{cc} 0 & \boldsymbol{B}_{1}^{(\mathrm{AS})} \\ \boldsymbol{B}_{2}^{(\mathrm{AS})} & 0 \end{array}\right)$ (C-2)
 $\boldsymbol{C}^{(\mathrm{AS})} =\left(\begin{array}{cc} 0 & \boldsymbol{C}_{1}^{(\mathrm{AS})} \\ \boldsymbol{C}_{2}^{(\mathrm{AS})} & 0 \end{array}\right)$ (C-3)

 $\boldsymbol{A}_{1}^{(\mathrm{AS})}=\rho^{-1}\left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \partial_{x} z_{0} & 0 & 0 \end{array}\right)$ (C-4)
 $\boldsymbol{A}_{2}^{(\mathrm{AS})}=\left(\begin{array}{ccc} c_{11} & c_{16} & c_{15} \\ c_{16} & c_{66} & c_{56} \\ c_{12} & c_{26} & c_{26} \end{array}\right)$ (C-5)
 $\boldsymbol{B}_{1}^{(\mathrm{AS})}=\rho^{-1}\left(\begin{array}{ccc} \partial_{x} z_{0} & 0 & 0 \\ 0 & \partial_{x} z_{0} & 0 \\ \left(\partial_{x} z_{0}\right)^{2} & 0 & 0 \end{array}\right)$ (C-6)
 $\boldsymbol{B}_{2}^{(\mathrm{AS})}=\left(\begin{array}{ccc} c_{15} & c_{14} & c_{13} \\ c_{56} & c_{46} & c_{36} \\ c_{25} & c_{24} & c_{23} \end{array}\right)$ (C-7)
 $\boldsymbol{C}_{1}^{(\mathrm{AS})}=\rho^{-1}\left(\begin{array}{ccc} 0 & \mathrm{i} k_{y} & 0 \\ 0 & 0 & \mathrm{i} k_{y} \\ \partial_{x x} z_{0} & \mathrm{i} k_{y} \partial_{x} z_{0} & 0 \end{array}\right)$ (C-8)
 $\boldsymbol{C}_{2}^{(\mathrm{AS})}=\mathrm{i} k_{y}\left(\begin{array}{lll} c_{16} & c_{12} & c_{14} \\ c_{66} & c_{26} & c_{46} \\ c_{26} & c_{22} & c_{24} \end{array}\right)$ (C-9)

 $\boldsymbol{A}^{(\mathrm{WS})} =\left(\begin{array}{cc} \bf{0} & \boldsymbol{A}_{1}^{(\mathrm{WS})} \\ \boldsymbol{A}_{2}^{(\mathrm{WS})} & \bf{0} \end{array}\right)$ (D-1)
 $\boldsymbol{B}^{(\mathrm{WS})} =\left(\begin{array}{cc} \bf{0} & \boldsymbol{B}_{1}^{(\mathrm{WS})} \\ \boldsymbol{B}_{2}^{(\mathrm{WS})} & \bf{0} \end{array}\right)$ (D-2)
 $\boldsymbol{C}^{(\mathrm{WS})} =\left(\begin{array}{cc} \bf{0} & \boldsymbol{C}_{1}^{(\mathrm{WS})} \\ \boldsymbol{C}_{2}^{(\mathrm{WS})} & \bf{0} \end{array}\right)$ (D-3)

 $\boldsymbol{A}_{1}^{(\mathrm{WS})}=\left(\begin{array}{ccc} -\rho_{\mathrm{w}}^{-1} & 0 & 0 \\ \frac{-\rho_{\mathrm{s}}^{-1}\left(\partial_{x} z_{0}\right)^{2}}{\left(\partial_{x} z_{0}\right)^{2}+1} & 0 & 0 \\ 0 & 0 & 0 \\ \frac{-\rho_{\mathrm{s}}^{-1} \partial_{x} z_{0}}{\left(\partial_{x} z_{0}\right)^{2}+1} & 0 & 0 \end{array}\right)$ (D-4)
 $\boldsymbol{A}_{2}^{(\mathrm{WS})}=-K\left(\begin{array}{llll} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array}\right)$ (D-5)
 $\boldsymbol{B}_{1}^{(\mathrm{WS})}=\rho_{\mathrm{s}}^{-1}\left(\begin{array}{ccc} 0 & 0 & 0 \\ \frac{\partial_{x} z_{0}}{\left(\partial_{x} z_{0}\right)^{2}+1} & 0 & 0 \\ 0 & -\partial_{x} z_{0} & 1 \\ \frac{-1}{\left(\partial_{x} z_{0}\right)^{2}+1} & 0 & 0 \end{array}\right)$ (D-6)
 $\boldsymbol{B}_{2}^{(\mathrm{WS})}=K\left(\begin{array}{cccc} -\partial_{x} z_{0} & \partial_{x} z_{0} & 0 & -1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array}\right)$ (D-7)
 $\boldsymbol{C}_{1}^{(\mathrm{WS})}=\frac{-\rho_{\mathrm{s}}^{-1}}{\left[\left(\partial_{x} z_{0}\right)^{2}+1\right]^{2}}\left(\begin{array}{ccc} 0 & 0 & 0 \\ -2 \partial_{x} z_{0} \partial_{x x} z_{0} & 0 & 0 \\ 0 & 0 & 0 \\ \partial_{x x} z_{0}\left[1-\left(\partial_{x} z_{0}\right)^{2}\right] & 0 & 0 \end{array}\right)$ (D-8)
 $\boldsymbol{C}_{2}^{(\mathrm{WS})}=-\mathrm{i} k_{y}\left[\left(\partial_{x} z_{0}\right)^{2}+1\right]\left(\begin{array}{cccc} 0 & c_{36} & c_{23} & c_{34} \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array}\right)$ (D-9)

 [1] 李青阳, 吴国忱, 段沛然. 准规则网格高阶有限差分法非均质弹性波波场模拟[J]. 石油地球物理勘探, 2019, 54(3): 539-550. LI Qingyang, WU Guochen, DUAN Peiran. Elastic wavefield forward modeling in heterogeneous media based on the quasi-regular grid high-order finite difference[J]. Oil Geophysical Prospecting, 2019, 54(3): 539-550. [2] Zhu J, Dorman J. Two-dimensional, three-component wave propagation in a transversely isotropic medium with arbitrary-orientation-finite-element modelling[J]. Geophysics, 2000, 65(3): 934-942. DOI:10.1190/1.1444789 [3] Ge Z, Chen X. An efficient approach for simulating wave propagation with the boundary element method in multilayered media with irregular interfaces[J]. Bulletin of the Seismological Society of America, 2008, 98(6): 3007-3016. DOI:10.1785/0120080920 [4] 谭文卓, 吴帮玉, 李博, 等. 梯形网格伪谱法地震波场模拟[J]. 石油地球物理勘探, 2020, 55(6): 1282-1291. TAN Wenzhuo, WU Bangyu, LI Bo, et al. Seismic wave simulation using a trapezoid grid pseudo-spectral method[J]. Oil Geophysical Prospecting, 2020, 55(6): 1282-1291. [5] 张庆朝, 朱国维, 周俊杰, 等. TTI介质qP波伪谱法正演模拟[J]. 石油地球物理勘探, 2019, 54(2): 302-311. ZHANG Qingchao, ZHU Guowei, ZHOU Junjie, et al. qP-wave numerical simulation in TTI media with pseudo-spectral method[J]. Oil Geophysical Prospecting, 2019, 54(2): 302-311. [6] Zhou B, Greenhalgh S A. 3-D frequency-domain seismic wave modelling in heterogeneous, anisotropic media using a Gaussian quadrature grid approach[J]. Geophysical Journal International, 2011, 184(1): 507-526. DOI:10.1111/j.1365-246X.2010.04859.x [7] Dormy E, Tarantola A. Numerical simulation of elastic wave propagation using a finite volume method[J]. Journal of Geophysical Research: Solid Earth, 1995, 100(B2): 2123-2133. DOI:10.1029/94JB02648 [8] Zhang Y, Duan L, Xie Y. A stable and practical implementation of least-squares reverse time migration[J]. Geophysics, 2015, 80(1): V23-V31. DOI:10.1190/geo2013-0461.1 [9] Li Y E, Demanet L. Full-waveform inversion with extrapolated low-frequency data[J]. Geophysics, 2016, 81(6): R339-R348. DOI:10.1190/geo2016-0038.1 [10] Zhou B, Greenhalgh S A, Hansruedi M. 2.5-D frequency-domain seismic wave modelling in heteroge-neous, anisotropic media using a Gaussian quadrature grid technique[J]. Computer and Geosciences, 2012, 39(1): 18-33. [11] Auer L, Nuber A M, Greenhalgh S A, et al. A critical appraisal of asymptotic 3D-to-2D data transformation in full-waveform seismic crosshole tomography[J]. Geophysics, 2013, 78(6): R235-R247. DOI:10.1190/geo2012-0382.1 [12] Williamson P R, Pratt R G. A critical review of the acoustic wave modelling procedure in 2.5 dimensions[J]. Geophysics, 1995, 60(2): 591-595. DOI:10.1190/1.1443798 [13] Song Z, Williamson P R. Frequency-domain acoustic-wave modeling and inversion of crosshole data: Part I, 2.5-D modeling method[J]. Geophysics, 1995, 60(3): 784-795. DOI:10.1190/1.1443817 [14] Zhou B, Greenhalgh S A. An adaptive wavenumber sampling strategy for 2.5D seismic-wave modelling in the frequency-domain[J]. Pure and Applied Geophysics, 2006, 163(7): 1399-1416. DOI:10.1007/s00024-006-0081-7 [15] Zhou B, Greenhalgh S A. A damping method for the computation of the 2.5-D Green's function for arbitrary acoustic media[J]. Geophysical Journal International, 1998, 133(1): 111-120. DOI:10.1046/j.1365-246X.1998.1331474.x [16] Maokun L, Vladimir D, Aria A, et al. A 2.5D finite-difference algorithm for elastic wave modeling using near-optimal quadratures[J]. Geophysics, 2016, 81(4): T155-T162. DOI:10.1190/geo2015-0550.1 [17] Takenaka H, Kennett B L N. A 2.5-D time-domain elastodynamic equation for plane-wave incidence[J]. Geophysical Journal International, 1996, 125(2): F5-F9. DOI:10.1111/j.1365-246X.1996.tb00001.x [18] Novais A, Santos L T. 2.5D finite-difference solution of the acoustic wave equation[J]. Geophysical Prospecting, 2005, 53(4): 523-531. DOI:10.1111/j.1365-2478.2005.00488.x [19] Xiong J L, Lin Y, Abubakar A, et al. 2.5-D forward and inverse modelling of full-waveform elastic seismic survey[J]. Geophysical Journal International, 2013, 193(2): 938-948. DOI:10.1093/gji/ggt013 [20] Zhang W, Chen X F. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[J]. Geophysical Journal International, 2006, 167(1): 337-353. DOI:10.1111/j.1365-246X.2006.03113.x [21] 李芳, 曹思远, 孙建国, 等. 2.5D傅氏变换法声波方程数值模拟及精度分析[J]. 石油地球物理勘探, 2009, 44(3): 276-281. LI Fang, CAO Siyuan, SUN Jianguo, et al. Numeric simulation and precision analysis of acoustic equation by 2.5D Fourier transform[J]. Oil Geophysical Prospecting, 2009, 44(3): 276-281. [22] Aki K, Richards P G. Quantitative Seismolog(2nd Edition)[M]. Sausalito, CA: University Science Books, 2002. [23] 蒋丽丽, 孙建国. 基于Poisson方程的曲网格生成技术[J]. 世界地质, 2008, 27(3): 298-305. JIANG Lili, SUN Jianguo. Source terms of elliptic system in grid generation[J]. Global Geology, 2008, 27(3): 298-305. DOI:10.3969/j.issn.1004-5589.2008.03.011 [24] Hixon R, Turkel E. Compact implicit MacCormack-type schemes with high accuracy[J]. Journal of Computational Physics, 2000, 158(1): 51-70. DOI:10.1006/jcph.1999.6406 [25] Tam C K, Webb J C. Dispersion-relation-preserving finite difference schemes for computational acoustics[J]. Journal of Computational Physics, 1993, 107(2): 262-281. DOI:10.1006/jcph.1993.1142 [26] Hixon R. Evaluation of a high-accuracy MacCormack-type scheme using benchmark problems[J]. Journal of Computational Acoustics, 1997, 6(3): 291-305. [27] Zhou B, Moosoo W, Greenhalgh S, et al. Generalized stiffness reduction method to remove the artificial edge-effects for seismic wave modelling in elastic anisotropic media[J]. Geophysical Journal Internatio-nal, 2020, 220(2): 1394-1408. DOI:10.1093/gji/ggz529 [28] Vavryuk V. Asymptotic green's function in homoge-neous anisotropic viscoelastic media[J]. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2007, 463(2): 2689-2707.