0 引言

地震波场数值模拟在地震勘探多个技术环节中扮演着重要角色，采用不同数值算法离散求解各种地震波动方程是地震数值模拟的核心。在众多数值算法中，有限差分法因其简单易行、计算高效和灵活性高等优点而广受青睐^[1-5]。但是，采用有限差分算子近似离散波动方程中的时间和空间偏导数时，通常会在时间域或空间域产生较强的数值频散。因此，提高有限差分在波场延拓中的精度、效率和灵活性一直是重要的研究内容^[6-9]。

现有差分方法可以通过增加差分算子长度提高空间域模拟精度。在传统规则网格高阶差分基础上，学者们发展了多种改进有限差分法，如交错网格与旋转交错网格差分方法^{[2, 8-11]}、隐式差分方法^[12-16]、优化差分方法^[17-19]、基于时间—空间域频散关系求取差分系数的高阶差分方法^[20-21]、波场外推在频率—空间域实现的差分方法^[22-23]等。然而这些改进方法主要聚焦于提高空间偏导数的精度，但时间导数的精度仍然较低，容易引起较强的时间频散。

在有限差分求解时间导数方面，直接采用上述高阶格式对时间导数进行离散容易导致不稳定^[12]，但可以采用不同的空间导数组合实现基于Lax-Wendroff算法的声波方程时间高阶数值模拟^[24]。本质上，Lax-Wendroff方案属于空间域方法，求解精度和效率有待进一步提高。相比之下，基于组合差分模板的时间—空间高阶差分法能够在兼顾效率的前提下有效提高时间模拟精度。在构建改进差分模板方面，Liu等^[25]提出了菱形网格有限差分法求解二维拉普拉斯算子，提高了声波方程在时间—空间域同步任意偶数阶模拟精度。但是该菱形网格中时间和空间导数的差分算子必须保持一致，因此限制了灵活性和效率。Wang等^[26]结合菱形网格与传统正交差分模板，发展了时间高阶差分算子与空间高阶差分算子相互独立的声波方程数值模拟算法，有效提高了计算效率。上述两种方案主要是针对形式较为简单的规则网格差分方法。Tan等^[27]基于菱形模板的思想，将交错网格旁轴上的节点与轴向节点联合用于求解导数，提出了针对声波方程的时间四阶和六阶差分方案；胡自多等^[28]对混合网格差分方法进行了推广与应用，提出了混合交错网格逆时偏移，有效提高了成像精度。通过构建组合差分模板的思想，目前已经发展了多种类型的时间—空间高阶有限差分方法，包括时间高阶弹性波差分^[29]、优化时间高阶差分^[30]、隐式时间高阶差分^[31]、针对三维波动方程的时间高阶差分^[32-33]等。但是，上述有限差分方法在构建改进模板时，主要采用正方形和立方体网格单元进行计算区域离散，即沿空间不同方向的离散网格间距一致，因此限制了模板的灵活性及进一步应用。

为了提高二维和三维声波方程时间—空间高阶有限差分数值模拟的灵活性，本文提出基于矩形/长方体网格单元(沿不同空间方向网格步长不一致)的改进时间—空间高阶有限差分法，分别用于求解二维和三维声波方程。基于平面波理论推导了相应的时间—空间域频散关系，进而应用泰勒展开求取了高阶差分系数。数值分析和多个模型算例验证了本文方法的有效性。

1 理论与方法 1.1 基于长方体网格单元的三维波动方程离散

三维各向同性介质中常密度声波方程表达式为

$ \frac{{\partial }^{2}P}{\partial {t}^{2}}={v}^{2}{\nabla }^{2}P $

(1)

式中：P=P(x, y, z, t)为声波波场；v=v(x, y, z, t)为介质速度；t为时间坐标；x、y和z为空间坐标；$ {\nabla }^{2} $=$ {\partial }^{2}/\partial {x}^{2}+ $$ {\partial }^{2}/\partial {y}^{2}+{\partial }^{2}/\partial {z}^{2}\mathrm{为}\mathrm{空}\mathrm{间}\mathrm{拉}\mathrm{普}\mathrm{拉}\mathrm{斯}\mathrm{算}\mathrm{子}, $$ {\partial }^{2}/\partial {x}^{2} $、$ {\partial }^{2}/ $$ \partial {y}^{2}\mathrm{和}{\partial }^{2}/\partial {z}^{2} $为空间不同方向的二阶偏导数。

采用有限差分法求解式(1)时，通常选用高阶差分格式对空间导数进行离散，即

$ \begin{array}{l}{\nabla }^{2}{P}_{i,l,f}^{o}\approx \frac{1}{{h}_{x}^{2}}\left[{a}_{0}{P}_{i,l,f}^{o}+\sum\limits_{m=1}^{M}{a}_{m}\left({P}_{i+m,l,f}^{o}+{P}_{i-m,l,f}^{o}\right)\right]+\\ \frac{1}{{h}_{y}^{2}}\left[{a}_{0}{P}_{i,l,f}^{o}+\sum\limits_{m=1}^{M}{a}_{m}\left({P}_{i,l+m,f}^{o}+{P}_{i,l-m,f}^{o}\right)\right]+\\ \frac{1}{{h}_{z}^{2}}\left[{a}_{0}{P}_{i,l,f}^{o}+\sum\limits_{m=1}^{M}{a}_{m}\left({P}_{i,l,f+m}^{o}+{P}_{i,l,f-m}^{o}\right)\right]\end{array} $

(2)

式中：i、l、f表示空间坐标网格索引；o表示时间坐标网格索引；M为空间项差分算子长度的一半；a_m为空间域高阶差分系数^[20]；h_x、h_y和h_z分别表示不同方向的空间网格步长。定义

$ \left\{\begin{array}{l}{R}_{yx}=\frac{{h}_{y}}{{h}_{x}}\\ {R}_{zx}=\frac{{h}_{z}}{{h}_{x}}\end{array}\right.  $

(3)

首选二阶形式对时间导数进行离散

$ \frac{{\partial }^{2}{P}_{i,l,f}^{o}}{\partial {t}^{2}}\approx \frac{1}{{\left(\mathrm{\Delta }t\right)}^{2}}\left(-2{P}_{i,l,f}^{o}+{P}_{i,l,f}^{o-1}+{P}_{i,l,f}^{o+1}\right) $

(4)

式中Δt为时间步长。

采用传统方案(式(2)和式(4))对式(1)进行数值离散能够达到空间域高阶精度，但时间域精度较低。如图 1所示，将轴向坐标轴上的网格节点(三个坐标轴上的黑色节点)与非轴向网格节点(蓝色八面体所包含的空间节点)进行有效结合。在传统差分方案的基础上，构建组合差分模板，得到改进的时间导数离散差分格式

$ \begin{array}{l}\frac{{\partial }^{2}{P}_{i,l,f}^{o}}{\partial {t}^{2}}\approx \frac{1}{\left(\mathrm{\Delta }{t}^{2}\right)}\left(-2{P}_{i,l,f}^{o}+{P}_{i,l,f}^{o-1}+{P}_{i,l,f}^{o+1}\right)-\\ \;\;\;\;\;\;\frac{{v}^{2}}{{h}_{x}^{2}}\left(\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{A}_{1}+\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-q}\sum\limits_{q=1}^{m-n}{A}_{2}\right)\end{array} $

(5)

$ \begin{array}{l}{A}_{1}={d}_{m,n,0}\left({P}_{i+m,l+n,f}^{o}+{P}_{i-m,l+n,f}^{o}+{P}_{i+m,l-n,f}^{o}+{P}_{i-m,l-n,f}^{o}\right)+\\ {d}_{m,0,n}\left({P}_{i+m,l,f+n}^{o}+{P}_{i-m,l,f+n}^{o}+{P}_{i+m,l,f-n}^{o}+{P}_{i-m,l,f-n}^{o}\right)+\\ {d}_{0,m,n}\left({P}_{i,l+m,f+n}^{o}+{P}_{i,l-m,f+n}^{o}+{P}_{i,l+m,f-n}^{o}+{P}_{i,l-m,f-n}^{o}\right)\end{array} $

(6)

$ \begin{array}{l}{A}_{2}={d}_{m,n,q}\left({P}_{i+m,l+n,f+q}^{o}+{P}_{i-m,l+n,f+q}^{o}+{P}_{i+m,l-n,f+q}^{o}+\right.\\ {P}_{i-m,l-n,f+q}^{o}+{P}_{i+m,l+n,f-q}^{o}+{P}_{i-m,l+n,f-q}^{o}+\\ \left.{P}_{i+m,l-n,f-q}^{o}+{P}_{i-m,l-n,f-q}^{o}\right)\end{array} $

(7)

式中：N为时间项差分算子长度的一半；d_{m, n, 0}、d_{m, 0, n}、d_{0, m, n}和d_{m, n, q}均为改进时间项离散格式中的高阶差分系数。式(5)的右端第一项为传统时间离散格式，A₁项和A₂项对应空间八面体上的节点，用于补偿传统时间离散格式的精度。

相应地，空间导数的离散格式修改为

$ {\nabla }^{2}{P}_{i,l,f}^{o}\approx \frac{{v}^{2}}{{h}_{x}^{2}}\left({d}_{\mathrm{0,0},0}{P}_{i,l,f}^{o}+\sum\limits_{m=1}^{M}{A}_{3}\right) $

(8)

$ \begin{array}{l}{A}_{3}={d}_{m,\mathrm{0,0}}\left({P}_{i+m,l,f}^{o}+{P}_{i-m,l,f}^{o}\right)+\\ {d}_{0,m,0}\left({P}_{i,l+m,f}^{o}+{P}_{i,l-m,f}^{o}\right)+\\ {d}_{\mathrm{0,0},m}\left({P}_{i,l,f+m}^{o}+{P}_{i,l,f-m}^{o}\right)\end{array} $

(9)

式(9)右端对应图 1中正交坐标轴上的节点；d_{0, 0, 0}、d_{m, 0, 0}、d_{0, m, 0}和d_{0, 0, m}为改进离散格式中的高阶差分系数。

将式(5)与式(8)结合，即可在时间域和空间域同步实现对式(1)的任意偶数阶精度差分求解。与前人成果不同，本文提出的时间高阶离散格式可以在不同方向采用非等间距的网格步长($ {h}_{x}\ne {h}_{y}\ne {h}_{z} $)，有效增加了差分模板的灵活性。

采用泰勒展开法对改进差分格式中的高阶系数进行求解。首先，基于平面波理论，标量波场P可表示为

$ {P}_{i,l,j}^{o}={\mathrm{e}}^{i\left[{k}_{x}\left(x+i{h}_{x}\right)+{k}_{y}\left(y+l{h}_{y}\right)+{k}_{z}\left(z+f{h}_{z}\right)-\omega \left(t+o\mathrm{\Delta }t\right)\right]} $

(10)

式中：ω为角频率；k_x、k_y和k_z分别为x、y和z方向的波数，且与总波数k之间具有如下关系

$ \left\{\begin{array}{l}\omega =kv\\ {k}_{x}=k\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{c}\mathrm{o}\mathrm{s}\varphi \\ {k}_{y}=k\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\varphi \\ {k}_{z}=k\mathrm{s}\mathrm{i}\mathrm{n}\theta \end{array}\right. $

(11)

式中：θ为从水平面测量的相对于垂向的平面波传播角度；ϕ为平面波方位角。

基于式(10)平面波解，将式(5)和式(8)代入式(1)，可进一步化简整理为

$ \begin{array}{l}{d}_{\mathrm{0,0},0}+2\sum\limits_{m=1}^{M}{A}_{4}+4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{A}_{5}+8\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{A}_{6}\\ \begin{array}{cc}& \end{array}={r}^{-2}\left[-2+2\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega \mathrm{\Delta }t\right)\right]\end{array} $

(12)

$ \begin{array}{l}{A}_{4}={d}_{m,\mathrm{0,0}}\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)+{d}_{0,m,0}\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{y}{h}_{y}\right)+\\ {d}_{\mathrm{0,0},m}\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{z}{h}_{z}\right)\end{array} $

(13)

$ \begin{array}{l}{A}_{5}={d}_{m,n,0}\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(n{k}_{y}{h}_{y}\right)+\\ \begin{array}{cc}& \end{array}{d}_{m,0,n}\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(n{k}_{z}{h}_{z}\right)+\\ \begin{array}{cc}& \end{array}{d}_{0,m,n}\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{y}{h}_{y}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(n{k}_{z}{h}_{z}\right)\end{array} $

(14)

$ {A}_{6}={d}_{m,n,q}\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(n{k}_{y}{h}_{y}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(q{k}_{z}{h}_{z}\right) $

(15)

式中：$ r=v\mathrm{\Delta }t/{h}_{x} $为x方向的库朗数；2M为空间差分算子长度。

采用泰勒级数对上式中的余弦函数进行展开，可得

$ \begin{array}{l}{d}_{\mathrm{0,0},0}+2\sum\limits_{m=1}^{M}\sum\limits_{j=0}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{j}{m}^{2j}}{\left(2j\right)!}{A}_{7}+\\ 4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}\sum\limits_{ε =0}^{\mathrm{\infty }}\sum\limits_{\delta =0}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{ε +\delta }{m}^{2ε }{n}^{2\delta }}{\left(2ε \right)!\left(2\delta \right)!}{A}_{8}+\\ 8\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{d}_{m,n,q}\sum\limits_{ε =0}^{\mathrm{\infty }}\sum\limits_{\delta =0}^{\mathrm{\infty }}\sum\limits_{\gamma =0}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{ε +\delta +\gamma }{m}^{2ε }{n}^{2\delta }{q}^{2\gamma }}{\left(2ε \right)!\left(2\delta \right)!\left(2\gamma \right)!}{A}_{9}\\ \approx 2\sum\limits_{j=1}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{j}{r}^{2j-2}{\left(k{h}_{x}\right)}^{2j}}{\left(2j\right)!}\end{array} $

(16)

$ {A}_{7}={d}_{m,\mathrm{0,0}}{k}_{x}^{2j}{h}_{x}^{2j}+{d}_{0,m,0}{k}_{y}^{2j}{h}_{y}^{2j}+{d}_{\mathrm{0,0},m}{k}_{z}^{2j}{h}_{z}^{2j} $

(17)

$ \begin{array}{l}{A}_{8}={d}_{m,n,0}{k}_{x}^{2ε }{k}_{y}^{2\delta }{h}_{x}^{2ε }{h}_{y}^{2\delta }+{d}_{m,0,n}{k}_{x}^{2ε }{k}_{z}^{2\delta }{h}_{x}^{2ε }{h}_{z}^{2\delta }+\\ \begin{array}{cc}& \end{array}{d}_{0,m,n}{k}_{y}^{2ε }{k}_{z}^{2\delta }{h}_{y}^{2ε }{h}_{z}^{2\delta }\end{array} $

(18)

$ {A}_{9}={k}_{x}^{2ε }{k}_{y}^{2\delta }{k}_{z}^{2\gamma }{h}_{x}^{2ε }{h}_{y}^{2\delta }{h}_{z}^{2\gamma } $

(19)

结合二项式定理

$ \begin{array}{l}{k}^{2j}={\left({k}_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}\right)}^{j}\\ =\sum\limits_{ε =0}^{j}\sum\limits_{\delta =0}^{ε }\frac{j!{k}_{x}^{2\left(j-ε \right)}{k}_{y}^{2\left(ε -\delta \right)}{k}_{z}^{2\delta }}{\left(j-ε \right)!\left(ε -\delta \right)!\delta !}\end{array} $

(20)

进一步将式(16)整理为

$ \begin{array}{l}{d}_{\mathrm{0,0},0}+2\sum\limits_{m=1}^{M}\left({d}_{m,\mathrm{0,0}}+{d}_{0,m,0}+{d}_{\mathrm{0,0},m}\right)+\\ 4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}\left({d}_{m,n,0}+{d}_{m,0,n}+{d}_{0,m,n}\right)+\\ 8\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{d}_{m,n,q}+2\sum\limits_{m=1}^{M}\sum\limits_{j=1}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{j}{m}^{2j}}{\left(2j\right)!}{A}_{10}{h}_{x}^{2j}+\\ 4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{A}_{11}+8\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{s=1}^{m-n}{d}_{m,n,q}{A}_{12}\\ \approx 2\sum\limits_{j=1}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{j}{r}^{2j-2}}{\left(2j\right)!}{A}_{13}{h}_{x}^{2j}+\\ 2\sum\limits_{j=1}^{\mathrm{\infty }}\sum\limits_{ε =1}^{j-1}\frac{{\left(-1\right)}^{j}j!{r}^{2j-2}}{\left(2j\right)!\left(j-ε \right)!ε !}{A}_{14}{h}_{x}^{2j}+\\ 2\sum\limits_{j=1}^{\mathrm{\infty }}\sum\limits_{ε =1}^{j-1}\sum\limits_{\delta =1}^{ε -1}\frac{{\left(-1\right)}^{j}j!{r}^{2j-2}}{\left(2j\right)!\left(j-ε \right)!\left(ε -\delta \right)!\delta !}{A}_{15}{h}_{x}^{2j}\end{array} $

(21)

$ {A}_{10}={d}_{m,\mathrm{0,0}}{k}_{x}^{2j}+{d}_{0,m,0}{k}_{y}^{2j}{R}_{yx}^{2j-2}+{d}_{\mathrm{0,0},m}{k}_{z}^{2j}{R}_{zx}^{2j-2} $

(22)

$ \begin{array}{l}{A}_{11}={d}_{m,n,0}\sum\limits_{j=1}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{j}\left({n}^{2j}{k}_{y}^{2j}{R}_{yx}^{2j}+{m}^{2j}{k}_{x}^{2j}\right)}{\left(2j\right)!}{h}_{x}^{2j}+\\ {d}_{m,0,n}\sum\limits_{j=1}^{\infty }\frac{{\left(-1\right)}^{j}\left({n}^{2j}{k}_{z}^{2j}{R}_{zx}^{2j}+{m}^{2j}{k}_{x}^{2j}\right)}{\left(2j\right)!}{h}_{x}^{2j}+\\ {d}_{0,m,n}\sum\limits_{j=1}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{j}\left({n}^{2j}{k}_{z}^{2j}{R}_{zx}^{2j}+{m}^{2j}{k}_{y}^{2j}{R}_{yx}^{2j}\right)}{\left(2j\right)!}{h}_{x}^{2j}+\\ \sum\limits_{ε =1}^{\mathrm{\infty }}\sum\limits_{\delta =1}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{ε +\delta }{m}^{2ε }{n}^{2\delta }}{\left(2ε \right)!\left(2\delta \right)!}{A}_{16}{h}_{x}^{2ε +2\delta }\end{array} $

(23)

$ \begin{array}{l}{A}_{12}=\sum\limits_{j=1}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{j}}{\left(2j\right)!}{A}_{17}{h}_{x}^{2j}+\\ \sum\limits_{ε =1}^{\mathrm{\infty }}\sum\limits_{\delta =1}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{ε +\delta }}{\left(2ε \right)!\left(2\delta \right)!}{A}_{18}{h}_{x}^{2ε +2\delta }+\\ \sum\limits_{ε =1}^{\mathrm{\infty }}\sum\limits_{\delta =1}^{\mathrm{\infty }}\sum\limits_{\gamma =1}^{\mathrm{\infty }}\frac{{\left(-1\right)}^{ε +\delta +\gamma }{m}^{2ε }{n}^{2\delta }{q}^{2\gamma }}{\left(2ε \right)!\left(2\delta \right)!\left(2\gamma \right)!}{A}_{19}{h}_{x}^{2ε +2\delta +2\gamma }\end{array} $

(24)

$ {A}_{13}={k}_{x}^{2j}+{k}_{y}^{2j}+{k}_{z}^{2j} $

(25)

$ {A}_{14}={k}_{x}^{2j-2ε }{k}_{y}^{2ε }+{k}_{x}^{2j-2ε }{k}_{z}^{2ε }+{k}_{y}^{2j-2ε }{k}_{z}^{2ε } $

(26)

$ {A}_{15}={k}_{x}^{2j-2ε }{k}_{y}^{2ε -2\delta }{k}_{z}^{2\delta } $

(27)

$ \begin{array}{l}{A}_{16}={d}_{m,n,0}{k}_{x}^{2ε }{k}_{y}^{2\delta }{R}_{yx}^{2\delta }+{d}_{m,0,n}{k}_{x}^{2ε }{k}_{z}^{2\delta }{R}_{zx}^{2\delta }+\\ {d}_{0,m,n}{k}_{y}^{2ε }{k}_{z}^{2\delta }{R}_{yx}^{2ε }{R}_{zx}^{2\delta }\end{array} $

(28)

$ {A}_{17}={m}^{2j}{k}_{x}^{2j}+{n}^{2j}{k}_{y}^{2j}{R}_{yx}^{2j}+{q}^{2j}{k}_{z}^{2j}{R}_{zx}^{2j} $

(29)

$ \begin{array}{l}{A}_{18}={m}^{2ε }{n}^{2\delta }{k}_{x}^{2ε }{k}_{y}^{2\delta }{R}_{yx}^{2\delta }+{m}^{2ε }{q}^{2\delta }{k}_{x}^{2ε }{k}_{z}^{2\delta }{R}_{zx}^{2\delta }+\\ \begin{array}{cc}& \end{array}{n}^{2ε }{q}^{2\delta }{k}_{y}^{2ε }{k}_{z}^{2\delta }{R}_{yx}^{2ε }{R}_{zx}^{2\delta }\end{array} $

(30)

$ {A}_{19}={k}_{x}^{2ε }{k}_{y}^{2\delta }{k}_{z}^{2\gamma }{R}_{yx}^{2\delta }{R}_{zx}^{2\gamma } $

(31)

比较式(21)两端关于$ {h}_{x}^{0} $和$ {h}_{x}^{2j}\left(j > 0\right) $项的系数，可得

$ \begin{array}{l}{d}_{\mathrm{0,0},0}+2\sum\limits_{m=1}^{M}\left({d}_{m,\mathrm{0,0}}+{d}_{0,m,0}+{d}_{\mathrm{0,0},m}\right)+\\ \begin{array}{cc}& \end{array}4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}\left({d}_{m,n,0}+{d}_{m,0,n}+{d}_{0,m,n}\right)+\\ \begin{array}{cc}& \end{array}8\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{d}_{m,n,q}=0\end{array} $

(32)

$ \begin{array}{l}\sum\limits_{m=1}^{M}\frac{{\left(-1\right)}^{j}{m}^{2j}}{\left(2j\right)!}{A}_{20}+\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}\left[\frac{{\left(-1\right)}^{j}}{\left(2j\right)!}\left({d}_{m,n,0}{A}_{21}+\right.\right.\\ \left.\left.{d}_{m,0,n}{A}_{22}+{d}_{0,m,n}{A}_{23}\right)+\sum\limits_{ε =1}^{j-1}\frac{{\left(-1\right)}^{j}{m}^{2j-2ε }{n}^{2ε }}{\left(2j-2ε \right)!\left(2ε \right)!}{A}_{24}\right]+\\ 4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{d}_{m,n,q}\left[\frac{{\left(-1\right)}^{j}{A}_{25}}{\left(2j\right)!}+\sum\limits_{ε =1}^{j-1}\frac{{\left(-1\right)}^{j}{A}_{26}}{\left(2j-2ε \right)!\left(2ε \right)!}+\right.\\ \left.\sum\limits_{ε =1}^{j-1}\sum\limits_{\delta =1}^{ε -1}\frac{{\left(-1\right)}^{j}{m}^{2j-2ε }{n}^{2ε -2\delta }{q}^{2\delta }}{\left(2j-2ε \right)!\left(2ε -2\delta \right)!\left(2\delta \right)!}{A}_{27}\right]\\ \approx \frac{{\left(-1\right)}^{j}{r}^{2j-2}{A}_{28}}{\left(2j\right)!}+\sum\limits_{ε =1}^{j-1}\frac{{\left(-1\right)}^{j}j!{r}^{2j-2}{A}_{29}}{\left(2j\right)!\left(j-ε \right)!ε !}+\\ \sum\limits_{ε =1}^{j-1}\sum\limits_{\delta =1}^{ε -1}\frac{{\left(-1\right)}^{j}j!{r}^{2j-2}{A}_{30}}{\left(2j\right)!\left(j-ε \right)!\left(ε -\delta \right)!\delta !}\end{array} $

(33)

$ {A}_{20}={d}_{m,\mathrm{0,0}}{k}_{x}^{2j}+{d}_{0,m,0}{k}_{y}^{2j}{R}_{yx}^{2j-2}+{d}_{\mathrm{0,0},m}{k}_{z}^{2j}{R}_{zx}^{2j-2} $

(34)

$ {A}_{21}={n}^{2j}{k}_{y}^{2j}{R}_{yx}^{2j}+{m}^{2j}{k}_{x}^{2j} $

(35)

$ {A}_{22}={n}^{2j}{k}_{z}^{2j}{R}_{zx}^{2j}+{m}^{2j}{k}_{x}^{2j} $

(36)

$ {A}_{23}={n}^{2j}{k}_{z}^{2j}{R}_{zx}^{2j}+{m}^{2j}{k}_{y}^{2j}{R}_{yx}^{2j} $

(37)

$ \begin{array}{l}{A}_{24}={d}_{m,n,0}{k}_{x}^{2j-2ε }{k}_{y}^{2ε }{R}_{yx}^{2ε }+{d}_{m,0,n}{k}_{x}^{2j-2ε }{k}_{z}^{2ε }{R}_{zx}^{2ε }+\\ \begin{array}{cc}& \end{array}{d}_{0,m,n}{k}_{y}^{2j-2ε }{k}_{z}^{2ε }{R}_{yx}^{2j-2ε }{R}_{zx}^{2ε }\end{array} $

(38)

$ {A}_{25}={m}^{2j}{k}_{x}^{2j}+{n}^{2j}{k}_{y}^{2j}{R}_{yx}^{2j}+{q}^{2j}{k}_{z}^{2j}{R}_{zx}^{2j} $

(39)

$ \begin{array}{l}{A}_{26}={m}^{2j-2ε }{n}^{2ε }{k}_{x}^{2j-2ε }{k}_{y}^{2ε }{R}_{yx}^{2ε }+\\ {m}^{2j-2ε }{q}^{2ε }{k}_{x}^{2j-2ε }{k}_{z}^{2ε }{R}_{zx}^{2ε }+\\ {n}^{2j-2ε }{q}^{2ε }{k}_{y}^{2ε }{k}_{z}^{2ε }{R}_{yx}^{2j-2ε }{R}_{zx}^{2ε }\end{array} $

(40)

$ {A}_{27}={k}_{x}^{2j-2ε }{k}_{y}^{2ε -2\delta }{k}_{z}^{2\delta }{R}_{yx}^{2ε -2\delta }{R}_{zx}^{2\delta } $

(41)

$ {A}_{28}={k}_{x}^{2j}+{k}_{y}^{2j}+{k}_{z}^{2j} $

(42)

$ {A}_{29}={k}_{x}^{2j-2ε }{k}_{y}^{2ε }+{k}_{x}^{2j-2ε }{k}_{z}^{2ε }+{k}_{y}^{2j-2ε }{k}_{z}^{2ε } $

(43)

$ {A}_{30}={k}_{x}^{2j-2ε }{k}_{y}^{2ε -2\delta }{k}_{z}^{2\delta } $

(44)

比较式(33)两端关于$ {k}_{x}^{2j} $、$ {k}_{y}^{2j} $、$ {k}_{z}^{2j} $、$ {k}_{x}^{2j-2ε }{k}_{y}^{2ε } $、$ {k}_{x}^{2j-2ε }{k}_{z}^{2ε } $、$ {k}_{y}^{2j-2ε }{k}_{z}^{2ε } $和$ {k}_{x}^{2j-2ε }{k}_{y}^{2ε -2\delta }{k}_{z}^{2\delta } $项的系数，可得

$ \begin{array}{l}\sum\limits_{m=1}^{M}{m}^{2j}{d}_{m,\mathrm{0,0}}+2\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{m}^{2j}\left({d}_{m,n,0}+{d}_{m,0,n}\right)+\\ \begin{array}{cc}& \end{array}4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{m}^{2j}{d}_{m,n,q}={r}^{2j-2}\end{array} $

(45)

$ \begin{array}{l}\sum\limits_{m=1}^{M}{m}^{2j}{d}_{0,m,0}+2\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}\left({n}^{2j}{d}_{m,n,0}+{m}^{2j}{d}_{0,m,n}\right)+\\ \begin{array}{cc}& \end{array}4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{n}^{2j}{d}_{m,n,q}={R}_{yx}^{-2j}{r}^{2j-2}\end{array} $

(46)

$ \begin{array}{l}\sum\limits_{m=1}^{M}{m}^{2j}{d}_{\mathrm{0,0},m}+2\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{n}^{2j}\left({d}_{m,0,n}+{d}_{0,m,n}\right)+\\ \begin{array}{cc}& \end{array}4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{q}^{2j}{d}_{m,n,q}={R}_{zx}^{-2j}{r}^{2j-2}\end{array} $

(47)

$ \begin{array}{l}\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{m}^{2j-2ε }{n}^{2ε }{d}_{m,n,0}+2\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{m}^{2j-2ε }{n}^{2ε }{d}_{m,n,q}\\ =\frac{j!\left(2j-2ε \right)!\left(2ε \right)!{R}_{yx}^{-2ε }{r}^{2j-2}}{2\left(2j\right)!\left(j-ε \right)!ε !}\end{array} $

(48)

$ \begin{array}{l}\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{m}^{2j-2ε }{n}^{2ε }{d}_{m,0,n}+2\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{m}^{2j-2ε }{n}^{2ε }{d}_{m,n,q}\\ =\frac{j!\left(2j-2ε \right)!\left(2ε \right)!{R}_{zx}^{-2ε }{r}^{2j-2}}{2\left(2j\right)!\left(j-ε \right)!ε !}\end{array} $

(49)

$ \begin{array}{l}\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{m}^{2j-2ε }{n}^{2ε }{d}_{0,m,n}+2\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{m}^{2j-2ε }{n}^{2ε }{d}_{m,n,q}\\ =\frac{j!\left(2j-2ε \right)!\left(2ε \right)!{R}_{yx}^{2ε -2j}{R}_{zx}^{-2ε }{r}^{2j-2}}{2\left(2j\right)!\left(j-ε \right)!ε !}\end{array} $

(50)

$ \begin{array}{l}\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{m}^{2j-2ε }{n}^{2ε -2\delta }{q}^{2\delta }{d}_{m,n,q}\\ =\frac{j!\left(2j-2ε \right)!\left(2ε -2\delta \right)!\left(2\delta \right)!{R}_{yx}^{2\delta -2ε }{R}_{zx}^{-2\delta }{r}^{2j-2}}{4\left(2j\right)!\left(j-ε \right)!\left(ε -\delta \right)!\delta !}\end{array} $

(51)

式中参数j、ε和δ的取值范围为

$ \left\{\begin{array}{l}j\in \left[1,N\right]\\ j > ε > \delta \end{array}\right.      $

(52)

综上所述，联合求解式(32)、式(45)~式(51)，即可获得改进模板中的高阶差分系数d_{0, 0, 0}、d_{m, 0, 0}、d_{0, m, 0}、d_{0, 0, m}、d_{m, n, 0}、d_{m, 0, n}、d_{0, m, n}和d_{m, n, q}。表 1给出了改进差分模板中时间四阶差分系数，相关参数为M=6, N=2, r=0.36, R_yx=1.2, R_zx=1.5。将相应差分系数代入式(5)和式(8)，即可实现三维地震波场高精度外推。

1.2 基于矩形网格单元的二维波动方程离散

二维常密度声波方程形式与方程(1)类似, 其中拉普拉斯算子的形式变为$ {\nabla }^{2}={\partial }^{2}/\partial {x}^{2}+{\partial }^{2}/\partial {z}^{2} $。

基于图 2中所示差分模板，联合应用非轴向和轴向的网格节点二维声波方程中的时间导数和空间导数，得到改进的波场递推公式

$ \begin{array}{l}\frac{1}{\left(\mathrm{\Delta }{t}^{2}\right)}\left(-2{P}_{i,l}^{o}+{P}_{i,l}^{o-1}+{P}_{i,l}^{o+1}\right)-\\ \frac{{v}^{2}}{{h}_{x}^{2}}\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{c}_{m,n}\left({P}_{i+m,l+n}^{o}+{P}_{i-m,l+n}^{o}+\right.\\ \left.{P}_{i+m,l-n}^{o}+{P}_{i-m,l-n}^{o}\right)\\ \approx \frac{{v}^{2}}{{h}_{x}^{2}}\left\{{c}_{\mathrm{0,0}}{P}_{i,l}^{o}+\left[{c}_{m,0}\left({P}_{i+m,l}^{o}+{P}_{i-m,l}^{o}\right)+\right.\right.\\ \left.{c}_{0,m}\left.\left({P}_{i,l+m}^{o}+{P}_{i,l-m}^{o}\right)\right]\right\}\end{array} $

(53)

式中c_{0, 0}、c_{m, 0}、c_{0, m}和c_{m, n}均为二维改进模板中差分系数。基于前文所述，结合平面波理论和泰勒级数展开即可求解二维改进差分格式中的高阶差分系数，此处不再赘述。

2 数值精度分析 2.1 相速度频散分析

相速度相对误差是一种常用的评价有限差分数值精度的方法。通过计算不同差分方法的相速度相对误差值，对其精度进行对比、分析。为了简便起见，对相关方法进行了简写表示(表 2)。

首先，定义参数$ \kappa $为数值相速度$ {v}_{\mathrm{F}\mathrm{D}} $与精确相速度$ v $的比值

$ \kappa =\frac{{v}_{\mathrm{F}\mathrm{D}}}{v}=\frac{1}{rk{h}_{x}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{-1}\left({r}^{2}Q+1\right)-1 $

(54)

采用$ \left|\kappa -1\right| $描述相速度频散误差，越接近于0，频散误差越小。算子Q在不同差分方案中形式不同，对于TFD-2D (2M, 2)格式

$ Q=\sum\limits_{m=1}^{M}{a}_{m}\left\{\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)-1\right]+{R}_{zx}^{-2}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{z}{h}_{z}\right)-1\right]\right\} $

(55)

对于MRFD-2D (2M, 2N)格式

$ \begin{array}{l}Q=\sum\limits_{m=1}^{M}\left\{{c}_{m,0}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)-1\right]+{c}_{0,m}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{z}{h}_{z}\right)-1\right]\right\}+\\ 2\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{c}_{m,n}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(n{k}_{y}{h}_{y}\right)-1\right]\end{array} $

(56)

对于TFD-3D (2M, 2)格式

$ \begin{array}{l}Q=\sum\limits_{m=1}^{M}{a}_{m}\left\{\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)-1\right]+\right.\\ \left.\begin{array}{cc}& \end{array}{R}_{yx}^{-2}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{y}{h}_{y}\right)-1\right]+{R}_{zx}^{-2}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{z}{h}_{z}\right)-1\right]\right\}\end{array} $

(57)

对于MCFD-3D (2M, 2N)格式

$ \begin{array}{l}Q=\sum\limits_{m=1}^{M}\left\{{d}_{m,\mathrm{0,0}}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)-1\right]+\right.\\ \left.{d}_{0,m,0}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{y}{h}_{y}\right)-1\right]+{d}_{\mathrm{0,0},m}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{z}{h}_{z}\right)-1\right]\right\}+\\ 2\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}\left\{{d}_{m,n,0}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(n{k}_{y}{h}_{y}\right)-1\right]+\right.\\ \begin{array}{cc}& {d}_{m,0,n}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(n{k}_{z}{h}_{z}\right)-1\right]+\end{array}\\ \begin{array}{cc}& \left.{d}_{0,m,n}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{y}{h}_{y}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(n{k}_{z}{h}_{z}\right)-1\right]\right\}+\end{array}\\ 4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{q=1}^{m-n}{d}_{m,n,q}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(m{k}_{x}{h}_{x}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(n{k}_{y}{h}_{y}\right)\times \right.\\ \left.\mathrm{c}\mathrm{o}\mathrm{s}\left(q{k}_{z}{h}_{z}\right)-1\right]\end{array} $

(58)

采用均匀二维和三维模型(v=3000 m/s)测试不同方法的数值频散精度。图 3为R_zx=1.5时二维不同方法计算的相速度频散曲面，Δt=1.5 ms、h_x=10 m。图 4为R_yx=1.2和R_zx=1.5时三维两种方法计算的相速度频散曲面，Δt=1.2 ms、h_x= 10 m。由图可见，在相同的网格比和角度参数条件下，TFD-2D (16, 2)和TFD-3D (16, 2)计算的频散曲面远离参考平面，表明其存在较强频散误差。而MRFD-2D (16, 2N)和MCFD-3D (16, 2N)的计算结果与参考平面更接近，能够有效压制频散误差的干扰。此外，在改进差分方案中，随着算子N的增加，模拟精度得到了进一步提升。

2.2 稳定性分析

稳定性条件是衡量有限差分方法优劣的一个重要因素。特征值分析通常被用于求解有限差分的稳定性因子^{[20, 25]}。基于特征值分析得到不同差分方法中的库朗数r满足如下稳定性约束条件

$ r\le s $

(59)

式中s为稳定性因子，其在不同差分方法中具有不同的计算形式，对于TFD-2D (2M, 2)格式

$ s={\left\{\frac{\left(1+{R}_{zx}^{-2}\right)}{2}\sum\limits_{m=1}^{M}{a}_{m}\left[1+{\left(-1\right)}^{m-1}\right]\right\}}^{-\frac{1}{2}} $

(60)

对于MRFD-2D (2M, 2N)格式

$ \begin{array}{l}s=\left\{\frac{1}{2}\sum\limits_{m=1}^{M}\left({c}_{m,0}+{c}_{0,m}\right)\left[1+{\left(-1\right)}^{m-1}\right]+\right.\\ {\left.\begin{array}{cc}& \end{array}\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}{c}_{m,n}\left[{\left(-1\right)}^{m+n-1}+1\right]\right\}}^{-\frac{1}{2}}\end{array} $

(61)

对于TFD-3D (2M, 2)格式

$ s={\left\{\frac{\left(1+{R}_{yx}^{-2}+{R}_{zx}^{-2}\right)}{2}\sum\limits_{m=1}^{M}{a}_{m}\left[1+{\left(-1\right)}^{m-1}\right]\right\}}^{-\frac{1}{2}} $

(62)

对于MCFD-3D (2M, 2N)格式

$ \begin{array}{l}s=\left\{\frac{1}{2}\sum\limits_{m=1}^{M}\left({d}_{m,\mathrm{0,0}}+{d}_{0,m,0}+{d}_{\mathrm{0,0},m}\right)\left[{\left(-1\right)}^{m-1}+1\right]+\right.\\ 2\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{N-m}\left({d}_{m,n,0}+{d}_{m,0,n}+{d}_{0,m,n}\right)\left[{\left(-1\right)}^{m+n-1}+1\right]+\\ {\left.4\sum\limits_{m=1}^{N-1}\sum\limits_{n=1}^{m-1}\sum\limits_{s=1}^{m-n}{d}_{m,n,l}\left[{\left(-1\right)}^{m+n+l-1}+1\right]\right\}}^{-\frac{1}{2}}\end{array} $

(63)

图 5、图 6为不同差分方案计算的二维和三维稳定性曲线图。从稳定性曲线对比结果可知：随着差分阶数（M）的增加，稳定性条件逐渐变差。当差分算子N长度一致时，相比之下，TFD-2D (2M, 2)和TFD-3D (2M, 2)具有最小的稳定性因子，表明传统方法容易出现不稳定。与传统方法相比，MRFD-2D (2M, 2N)和MCFD-3D(2M, 2N)具有更高的稳定性条件，因此当空间步长确定时可以选用更大的时间步长，有效提高计算效率。此外，稳定性随着N的增大而呈递增趋势。需要说明的是，图中的稳定性曲线均基于x方向的空间步长$ {h}_{x} $计算，其余方向的稳定性因子可以采用式(3)换算。

3 模型算例 3.1 均匀模型

首先采用均匀模型对相关方法的数值精度进行测试。二维模型尺寸为4000 m×4800 m，矩形离散网格单元尺寸为10 m×12 m，介质速度为2800 m/s，震源采用为主频34 Hz的雷克子波。三维模型尺寸为2000 m×2400 m×2500 m，长方体离散网格单元尺寸为10.0 m×12.0 m×12.5 m，介质速度为2350 m/s，震源采用主频为25 Hz的雷克子波。震源均置于模型中央。

图 7和图 8分别为不同差分方法计算的均匀介质二维和三维波场切片。由图可见，传统时间二阶差分法(TFD-2D(16, 2)和TFD-3D(12, 2))计算的波场切片出现较强的时间频散误差干扰(如图中箭头所示)。相比之下，改进时间高阶差分法（MRFD-2D(16, 6)和MCFD-3D(12, 6))能够在保证空间模拟精度的前提下，显著降低时间频散误差。当模型参数一定时，改进差分方案在保证模拟精度的同时采用更大的时间步长，可以在一定程度上提高计算效率。

3.2 复杂模型

采用修改后的山麓模型对二维差分法的精度进行测试。如图 9所示，模型包含834×500个矩形网格单元，每个网格单元的尺寸为10 m×12 m。震源选用主频为28 Hz的雷克子波，置于(4170 m，120 m)。图 10a为传统方法采用较小时间步长计算出的参考地震记录局部显示，图 10b和图 10c分别为传统方法和改进方法采用大时间步长计算的局部地震记录，图 11为不同方法单道波形对比图。为了进一步凸显三维差分法在复杂模型正演中的优势，利用公式$ v=\sqrt[3]{v/1500}\times 2000 $对原始二维盐丘模型进行修改，得到改进的速度模型如图 12所示，将其沿y方向进行延拓，得到修改后的三维盐丘速度模型。模型尺寸为4000 m×4800 m×2125 m，长方体离散网格单元尺寸为10.0 m×12.0 m×12.5 m。震源采用主频为26 Hz的雷克子波，置于(2000.0 m，2400.0 m，187.5 m)。图 13和图 14分别为三维盐丘模型利用传统方法和改进方法计算的地震记录和单道波形对比。

对比复杂模型不同方法计算的地震记录与单道波形可知，传统差分法(TFD-2D (16, 2)和TFD-3D (12, 2))采用相对较小的时间步长时仍然可以观察到比较明显的数值频散(如图中箭头所示)，与参考波形存在较大相移误差。相比之下，当模型参数一致时，改进差分法(MRFD-2D (16, 6)和MCFD-3D (12, 4))采用相对较大的时间步长时单道记录与参考波形匹配较好，未出现明显数值频散。此外，在相同的计算平台上，改进方法能够采用更大的时间步长获得较高的模拟精度，既减少了迭代次数，又保证了计算效率。

综上所述，复杂模型算例表明，改进方案能够采用较大的时间步长和不一致的网格间距实现地震波场高精度数值模拟，证明了本文方法具有高效性和灵活性。

4 结论

通过结合差分模板中旁轴节点和轴向节点，本文提出了一种基于矩形网格和长方体网格离散单元的改进时间—空间高阶有限差分方法，分别用于求解二维和三维声波方程。结合平面波解和泰勒级数展开求取了改进差分方案中的高阶时空域差分系数。

数值分析表明本文方法具有更高的相速度频散精度和稳定性。多个模型算例证明当模型参数相同时，改进差分方案能够在保证空间模拟精度的前提下有效提高时间模拟精度。此外，与传统方法相比，本文方法能够采用较大的时间步长和不同的网格间距，具有更高的效率和灵活性。