0 引言

波动方程数值模拟是研究复杂介质中波场传播规律的重要手段^[1-2]，是逆时偏移(RTM)^[3-4]和全波形反演(FWI)^[5-6]的重要基础。相比伪谱法^[7-8]和有限元法^[9-10]，有限差分法^[11-12]具有计算效率高、占用内存小、算法实现简单灵活等优点，已发展成为应用最普遍的一种波动方程数值模拟方法^[13-15]。然而，有限差分法利用差分算子近似波动方程中的微分算子，导致波场的传播速度与真实速度不相等，且不同频率的波场传播速度不同，即出现数值频散现象^[16]。固有的数值频散严重影响有限差分法的模拟精度^[17-18]，进而影响了RTM和FWI的精度^[19]，因此，压制数值频散、提高模拟精度是有限差分法的一项重要研究内容。

构建更合理的差分算子和改进差分系数计算方法是提高有限差分法模拟精度的两条重要途径。Dablain^[16]指出，高阶差分算子能够减小数值频散，提高模拟精度，但时间高阶差分算子会使内存占用和计算量显著增加。此后，学者普遍保持时间2阶差分算子不变，通过设计更合理的空间差分算子以提高模拟精度。Fornberg^[20]给出了基于泰勒级数展开的规则网格和交错网格2M阶(M表示空间差分算子中与差分中心点等距的坐标轴网格点的组数)空间差分算子的差分系数解析表达式，据此可以构建采用时间2阶差分算子、空间2M阶差分算子的规则网格和交错网格有限差分法，本文称之为常规规则网格有限差分法(C-FD)和常规交错网格有限差分法(C-SFD)。C-FD和C-SFD利用空间域频散关系和泰勒级数展开计算差分系数，尽管空间差分算子具有2M阶差分精度，但差分离散波动方程仍然仅具有2阶差分精度。波动方程有限差分数值模拟通过迭代求解差分离散波动方程实现，不应通过分开逼近时间差分算子或空间差分算子达到高阶差分精度，而应该使差分离散波动方程达到高阶差分精度，才能有效提高有限差分法的模拟精度。针对二阶标量波动方程和一阶应力—速度声波方程有限差分模拟，一些学者提出基于时空间域频散关系和泰勒级数展开计算差分系数，构建了时空域规则网格和交错网格有限差分法，这种时空域有限差分法使得相应的二维和三维差分离散波动方程分别沿8个和48个传播方向达到2M阶差分精度，但沿其他传播方向仍然仅具有2阶差分精度，呈现明显的数值各向异性^{[17, 21]}。除了上述基于泰勒展开的差分系数算法，基于最小二乘优化算法最小化频散关系误差，相速度误差和群速度误差求解差分系数也被广泛采用^[22-24]，这种基于最小二乘的差分系数算法能够有效提高中高频成分的模拟精度，但会在一定程度上损失低频成分的模拟精度，并且最小二乘算法通过迭代优化差分系数，计算量较大。Liu^[25-26]采用最小二乘优化相对空间域和时空域频散关系求解差分系数，将非线性最优化问题转化为线性最优化问题，求解过程不需要迭代，有效地提高了差分系数的计算效率。

仅通过改进差分系数算法提高模拟精度的效果有限，构建更合理的空间差分算子是有效提高模拟精度的另一重要途径。针对2阶标量波动方程，Liu等^[27]提出了一种菱形网格有限差分法，并基于时空域频散关系和泰勒展开计算差分系数，使得差分离散波动方程沿任意传播方向达到2M阶差分精度，有效地提高了二维标量波动方程的模拟精度和稳定性，但是其空间差分算子长度随M²急剧增大，导致计算量巨大，计算效率低。后来，Wang等^[28]又提出组合常规十字交叉型网格和菱形网格以构建空间差分算子，有效兼顾了计算效率和模拟精度。胡自多等^[29]借鉴频率域混合网格有限差分法^[30-31]的构建思路，将波动方程中的Laplace微分算子近似表示为直角坐标系中坐标轴网格点构建的Laplace差分算子和非坐标轴网格点构建的Laplace差分算子的加权平均，建立了一种混合网格有限差分法，它与Wang等^[28]提出的有限差分法相似，但构建思路不同。胡自多等^[32]导出了三维直角坐标系中非坐标轴网格点构建Laplace差分算子的方法，进一步构建了三维混合网格有限差分法，有效提高了三维标量波动方程的模拟精度和稳定性。针对速度—应力声波方程，Tan等^[33]提出联合利用坐标轴网格点和非坐标轴网格点构建空间差分算子近似波动方程中的1阶空间偏导数，建立了一种混合交错网格有限差分法，相应的差分离散声波方程可以达到4阶或6阶差分精度，有效提高了声波模拟精度。Ren等^[34-35]进一步发展了这种混合交错网格有限差分法，使得差分离散声波和弹性波方程最高可以达到8阶差分精度。Zhou等^[36-37]对混合交错网格有限差分法做进一步改进，使得差分离散声波和弹性波方程可以达到任意偶数阶差分精度，并采用两步线性优化方法简化了基于最小二乘的优化差分系数计算，进一步提高了声波和弹性波的模拟精度。然而，这类混合交错网格有限差分法不恰当地使用了非坐标轴网格点的对称性，通常将与差分中心点距离不相等的2组非坐标轴网格点赋予了相同的差分系数，导致难以推导差分系数的通解。

针对速度—应力声波方程数值模拟，本文发展了一种改进型的混合交错网格有限差分法(M-SFD)，利用坐标轴网格点和非坐标轴网格点构建空间差分算子时，确保与差分中心点距离相等的1组非坐标轴网格点赋予相同的差分系数，与差分中心点距离不等的任意2组非坐标轴网格点赋予不同的差分系数，使得所构建的M-SFD更为合理，并利用时空域频散关系和泰勒级数展开建立差分系数求解方程组，导出差分系数通解的解析表达式。

本文首先阐述了M-SFD中空间差分算子的构建方法，给出相应的差分离散声波方程，并导出差分系数通解；其次，开展频散分析和稳定性分析；最后，利用层状介质模型和中国塔里木盆地典型复杂构造模型开展数值模拟，并将M-SFD推广应用于RTM，利用塔里木复杂构造模型模拟数据开展RTM测试。测试表明，该方法能够有效消除由于数值频散造成的成像假象，从而提高深层成像精度和分辨率。

1 混合交错网格有限差分法的基本原理 1.1 差分离散声波方程的导出

二维速度—应力声波方程可表示为

$ \begin{array}{lll}\left\{\begin{array}{l}\frac{\partial P}{\partial t}+\kappa \left(\frac{\partial {\upsilon }_{x}}{\partial x}+\frac{\partial {\upsilon }_{z}}{\partial z}\right)=0\\ \frac{\partial {\upsilon }_{x}}{\partial t}+\frac{1}{\rho }\frac{\partial P}{\partial x}=0\\ \frac{\partial {\upsilon }_{z}}{\partial t}+\frac{1}{\rho }\frac{\partial P}{\partial z}=0\end{array}\right.& \begin{array}{l}\\ \end{array}& \end{array} $

(1)

式中：$ P=P\left(x, z, t\right) $为压力场；$ {\upsilon }_{x}={\upsilon }_{x}\left(x, z, t\right) $和$ {\upsilon }_{z}={\upsilon }_{z}\left(x, z, t\right) $分别为质点振动速度场的$ x $和$ z $分量；$ \kappa =\kappa \left(x, z\right) $为体积模量；$ \rho =\rho \left(x, z\right) $为介质密度。

交错网格有限差分法求解式(1)时，波场变量和弹性参数定义在交错的网格位置上，如图 1所示。波场变量$ P\left(x, z, t\right) $差分离散点定义在空间$ x $和$ z $方向的半网格点、时间$ t $方向的整网格点上(图 1红色圆圈)，其离散表达式为$ {P}_{m-1/2, n-1/2}^{j}=P\left[x+\left(m-1/2\right)h, z+\left(n-1/2\right)h, t+j\mathrm{\Delta }t\right] $，$ m, n, j $为整数，$ h $为空间采样间隔，$ \mathrm{\Delta }t $为时间采样间隔。

与C-SFD类似，M-SFD采用时间2阶差分算子近似式(1)中波场变量关于时间变量$ t $的1阶偏导数，$ \partial P/\partial t $、$ \partial {\upsilon }_{x}/\partial t $和$ \partial {\upsilon }_{z}/\partial t $可近似表示为

$ \left\{\begin{array}{l}{\left.\frac{\partial P}{\partial t}\right|}_{1/2, 1/2}^{1/2}\approx \frac{{P}_{1/2, 1/2}^{1}-{P}_{1/2, 1/2}^{0}}{\mathrm{\Delta }t}\\ {\left.\frac{\partial {\upsilon }_{x}}{\partial t}\right|}_{0, 1/2}^{0}\approx \frac{{\upsilon }_{x\left(0, 1/2\right)}^{1/2}-{\upsilon }_{x\left(0, 1/2\right)}^{-1/2}}{\mathrm{\Delta }t}\\ {\left.\frac{\partial {\upsilon }_{z}}{\partial t}\right|}_{1/2, 0}^{0}\approx \frac{{\upsilon }_{z\left(1/2, 0\right)}^{1/2}-{\upsilon }_{z\left(1/2, 0\right)}^{-1/2}}{\mathrm{\Delta }t}\end{array}\right. $

(2)

式中：$ {P}_{1/2, 1/2}^{0} $代表$ m=1, n=1, j=0 $时$ {P}_{m-1/2, n-1/2}^{j} $的值；$ {\upsilon }_{x\left(m, n-1/2\right)}^{j-1/2} $和$ {\upsilon }_{z\left(m-1/2, n\right)}^{j-1/2} $具有相似的表达式。

C-SFD仅利用坐标轴网格点构建空间2M阶差分算子(图 2a)近似波场变量关于空间的一阶偏导数，随着M取值的增大，新增网格点与差分中心点的距离也逐渐增大，对提高模拟精度的贡献逐渐减小。

本文所提M-SFD的基本思路是：联合利用坐标轴网格点和非坐标轴网格点构建一种混合型的空间差分算子近似波场变量关于空间的1阶偏导数。图 2b~图 2e为M-SFD(N=1、2、3、4)的空间差分算子示意图，N表示空间差分算子中与差分中心点等距的非坐标轴网格点的组数。相比C-SFD，M-SFD能有效利用距离差分中心点更近的非坐标轴网格点，理论上更合理。

本文提出的M-SFD与Tan等^[33]提出的时间和空间高阶交错网格有限差分法具有一定的相似性，但是他们不恰当地使用了非坐标轴网格点的对称性，将与差分中心点距离不相等的两组非坐标轴网格点赋予了相同的差分系数，例如将图 2d中标记为绿色②和③的两组非坐标轴网格点赋予了相同的差分系数，标记为绿色②的网格点与差分中心点的距离为$ \sqrt[]{13}h/2 $，而标记为绿色③的网格点与差分中心点的距离为$ \sqrt[]{17}h/2 $。这种不合理的赋值导致差分系数的解析解求解困难。本文提出的M-SFD给任意两组与差分中心点距离不等的网格点赋予不同的差分系数，理论上更合理，并且求解差分系数的解析解更容易。

M-SFD(N=1)与Tan等^[33]提出的时间4阶、空间2M阶交错网格有限差分法本质上完全相同。下面以M-SFD(N=2)为例阐述M-SFD的基本原理。利用图 2c中M-SFD(N=2)的空间差分算子近似方程(1)中波场变量关于空间变量$ x $和$ z $的一阶导数，$ \partial {\upsilon }_{x}/\partial x $、$ \partial {\upsilon }_{z}/\partial z $、$ \partial P/\partial x $和$ \partial P/\partial z $可以表示为

$ \left\{\begin{array}{l}{\left.\frac{\partial {\upsilon }_{x}}{\partial x}\right|}_{1/2, 1/2}^{1/2}\approx \frac{1}{h}\left\{\begin{array}{l}\sum\limits _{m=1}^{M}{a}_{m}\left[{\upsilon }_{x\left(m, 1/2\right)}^{1/2}-{\upsilon }_{x\left(-m+1, 1/2\right)}^{1/2}\right]+{b}_{1}\left[{\upsilon }_{x\left(1, 3/2\right)}^{1/2}-{\upsilon }_{x\left(0, 3/2\right)}^{1/2}+{\upsilon }_{x\left(1, -1/2\right)}^{1/2}-{\upsilon }_{x\left(0, -1/2\right)}^{1/2}\right]+\\ {b}_{2}\left[{\upsilon }_{x\left(2, 3/2\right)}^{1/2}-{\upsilon }_{x\left(-1, 3/2\right)}^{1/2}+{\upsilon }_{x\left(2, -1/2\right)}^{1/2}-{\upsilon }_{x\left(-1, -1/2\right)}^{1/2}\right]\end{array}\right\}\\ {\left.\frac{\partial {\upsilon }_{z}}{\partial z}\right|}_{1/2, 1/2}^{1/2}\approx \frac{1}{h}\left\{\begin{array}{l}\sum\limits _{m=1}^{M}{a}_{m}\left[{\upsilon }_{z\left(1/2, m\right)}^{1/2}-{\upsilon }_{z\left(1/2, -m+1\right)}^{1/2}\right]+{b}_{1}\left[{\upsilon }_{z\left(3/2, 1\right)}^{1/2}-{\upsilon }_{z\left(3/2, 0\right)}^{1/2}+{\upsilon }_{z\left(-1/2, 1\right)}^{1/2}-{\upsilon }_{z\left(-1/2, 0\right)}^{1/2}\right]+\\ {b}_{2}\left[{\upsilon }_{z\left(3/2, 2\right)}^{1/2}-{\upsilon }_{z\left(3/2, -1\right)}^{1/2}+{\upsilon }_{z\left(-1/2, 2\right)}^{1/2}-{\upsilon }_{z\left(-1/2, -1\right)}^{1/2}\right]\end{array}\right\}\\ {\left.\frac{\partial P}{\partial x}\right|}_{0, 1/2}^{0}\approx \frac{1}{h}\left\{\begin{array}{l}\sum\limits _{m=1}^{M}{a}_{m}\left[{P}_{m-1/2, 1/2}^{0}-{P}_{-m+1/2, 1/2}^{0}\right]+{b}_{1}\left[{P}_{1/2, 3/2}^{0}-{P}_{-1/2, 3/2}^{0}+{P}_{1/2, -1/2}^{0}-{P}_{-1/2, -1/2}^{0}\right]+\\ {b}_{2}\left[{P}_{3/2, 3/2}^{0}-{P}_{-3/2, 3/2}^{0}+{P}_{3/2, -1/2}^{0}-{P}_{-3/2, -1/2}^{0}\right]\end{array}\right\}\\ {\left.\frac{\partial P}{\partial z}\right|}_{1/2, 0}^{0}\approx \frac{1}{h}\left\{\begin{array}{l}\sum\limits _{m=1}^{M}{a}_{m}\left[{P}_{1/2, m-1/2}^{0}-{P}_{1/2, -m+1/2}^{0}\right]+{b}_{1}\left[{P}_{3/2, 1/2}^{0}-{P}_{3/2, -1/2}^{0}+{P}_{-1/2, 1/2}^{0}-{P}_{-1/2, -1/2}^{0}\right]+\\ {b}_{2}\left[{P}_{3/2, 3/2}^{0}-{P}_{3/2, -3/2}^{0}+{P}_{-1/2, 3/2}^{0}-{P}_{-1/2, -3/2}^{0}\right]\end{array}\right\}\end{array}\right. \begin{array}{l}\\ \\ \\ \end{array} $

(3)

式中$ {a}_{m} $(m=1，2，…，M)、$ {b}_{1}、{b}_{2} $均为差分系数。将式(2)和式(3)代入式(1)得到

$ \left\{\begin{array}{l}\frac{{P}_{1/2, 1/2}^{1}-{P}_{1/2, 1/2}^{0}}{\mathrm{\Delta }t}\approx -\frac{\kappa }{h}\left\{\begin{array}{l}\sum\limits _{m=1}^{M}{a}_{m}\left[{\upsilon }_{x\left(m, 1/2\right)}^{1/2}-{\upsilon }_{x\left(-m+1, 1/2\right)}^{1/2}\right]+{b}_{1}\left[{\upsilon }_{x\left(1, 3/2\right)}^{1/2}-{\upsilon }_{x\left(0, 3/2\right)}^{1/2}+{\upsilon }_{x\left(1, -1/2\right)}^{1/2}-{\upsilon }_{x\left(0, -1/2\right)}^{1/2}\right]+\\ {b}_{2}\left[{\upsilon }_{x\left(2, 3/2\right)}^{1/2}-{\upsilon }_{x\left(-1, 3/2\right)}^{1/2}+{\upsilon }_{x\left(2, -1/2\right)}^{1/2}-{\upsilon }_{x\left(-1, -1/2\right)}^{1/2}\right]\end{array}\right\}-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{\kappa }{h}\left\{\begin{array}{l}\sum\limits _{m=1}^{M}{a}_{m}\left[{\upsilon }_{z\left(1/2, m\right)}^{1/2}-{\upsilon }_{z\left(1/2, -m+1\right)}^{1/2}\right]+{b}_{1}\left[{\upsilon }_{z\left(3/2, 1\right)}^{1/2}-{\upsilon }_{z\left(3/2, 0\right)}^{1/2}+{\upsilon }_{z\left(-1/2, 1\right)}^{1/2}-{\upsilon }_{z\left(-1/2, 0\right)}^{1/2}\right]+\\ {b}_{2}\left[{\upsilon }_{z\left(3/2, 2\right)}^{1/2}-{\upsilon }_{z\left(3/2, -1\right)}^{1/2}+{\upsilon }_{z\left(-1/2, 2\right)}^{1/2}-{\upsilon }_{z\left(-1/2, -1\right)}^{1/2}\right]\end{array}\right\}\\ \frac{{\upsilon }_{x\left(0, 1/2\right)}^{1/2}-{\upsilon }_{x\left(0, 1/2\right)}^{-1/2}}{\mathrm{\Delta }t}\approx -\frac{1}{\rho h}\left\{\begin{array}{l}\sum\limits _{m=1}^{M}{a}_{m}\left({P}_{m-1/2, 1/2}^{0}-{P}_{-m+1/2, 1/2}^{0}\right)+{b}_{1}\left[{P}_{1/2, 3/2}^{0}-{P}_{-1/2, 3/2}^{0}+{P}_{1/2, -1/2}^{0}-{P}_{-1/2, -1/2}^{0}\right]+\\ {b}_{2}\left[{P}_{3/2, 3/2}^{0}-{P}_{-3/2, 3/2}^{0}+{P}_{3/2, -1/2}^{0}-{P}_{-3/2, -1/2}^{0}\right]\end{array}\right\}\\ \frac{{\upsilon }_{z\left(1/2, 0\right)}^{1/2}-{\upsilon }_{z\left(1/2, 0\right)}^{-1/2}}{\mathrm{\Delta }t}\approx -\frac{1}{\rho h}\left\{\begin{array}{l}\sum\limits _{m=1}^{M}{a}_{m}\left({P}_{1/2, m-1/2}^{0}-{P}_{1/2, -m+1/2}^{0}\right)+{b}_{1}\left[{P}_{3/2, 1/2}^{0}-{P}_{3/2, -1/2}^{0}+{P}_{-1/2, 1/2}^{0}-{P}_{-1/2, -1/2}^{0}\right]+\\ {b}_{2}\left[{P}_{3/2, 3/2}^{0}-{P}_{3/2, -3/2}^{0}+{P}_{-1/2, 3/2}^{0}-{P}_{-1/2, -3/2}^{0}\right]\end{array}\right\}\end{array}\right. $

(4)

式(4)为M-SFD(N=2)对式(1)的差分离散声波方程，同样还可以导出M-SFD(N=1、3、4)对式(1)的差分离散声波方程。

1.2 差分系数计算

合理的差分系数算法能有效提高交错网格有限差分法的模拟精度。C-SFD利用空间域频散关系和泰勒级数展开计算差分系数，使得空间差分算子达到2M阶差分精度，但是相应的差分离散声波方程仅具有2阶差分精度^[21]。交错网格有限差分法通过迭代求解差分离散声波方程实现声波数值模拟，为了提高模拟精度，应该设法提高差分离散声波方程的差分精度，而不是仅仅提高空间差分算子的差分精度。

M-SFD基于时空域频散关系和泰勒级数展开计算差分系数，将有助于提高差分离散声波方程的差分精度。下面以M-SFD(N=2)为例，阐述M-SFD的差分系数计算方法。均匀介质中，速度—应力声波方程(1)具有如下形式的离散平面波解

$ \left\{\begin{array}{l}{P}_{m-1/2, n-1/2}^{j}={A}_{P}{\mathrm{e}}^{i\left[{k}_{x}\left(x+\left(m-1/2\right)h\right)+{k}_{z}\left(z+\left(n-1/2\right)h\right)-\omega \left(t+j\mathrm{\Delta }t\right)\right]}\\ {\upsilon }_{x\left(m, n-1/2\right)}^{j-1/2}={A}_{{\upsilon }_{x}}{\mathrm{e}}^{i\left[{k}_{x}\left(x+mh\right)+{k}_{z}\left(z+\left(n-1/2\right)h\right)-\omega \left(t+\left(j-1/2\right)\mathrm{\Delta }t\right)\right]}\\ {\upsilon }_{z\left(m-1/2, n\right)}^{j-1/2}={A}_{{\upsilon }_{z}}{\mathrm{e}}^{i\left[{k}_{x}\left(x+\left(m-1/2\right)h\right)+{k}_{z}\left(z+nh\right)-\omega \left(t+\left(j-1/2\right)\mathrm{\Delta }t\right)\right]}\\ {k}_{x}=k\mathrm{c}\mathrm{o}\mathrm{s}\theta \\ {k}_{z}=k\mathrm{s}\mathrm{i}\mathrm{n}\theta \end{array}\right. $

(5)

式中：$ {A}_{P} $、$ {A}_{{\upsilon }_{x}} $和$ {A}_{{\upsilon }_{z}} $为平面波的振幅；$ k $为波数；$ \theta $为平面波传播方向与$ x $轴正向的夹角；$ \omega $为角频率。

将式(5)代入式(4)得到

$ \left\{\begin{array}{l}\frac{{A}_{P}}{\mathrm{\Delta }t}\mathrm{s}\mathrm{i}\mathrm{n}\frac{\omega \mathrm{\Delta }t}{2}\\ \approx -\frac{\kappa {A}_{{\upsilon }_{x}}}{h}\left\{\sum\limits _{m=1}^{M}{a}_{m}\mathrm{s}\mathrm{i}\mathrm{n}\left[\left(m-\frac{1}{2}\right){k}_{x}h\right]+\right.\\ \;\;\;\;\left.2\mathrm{c}\mathrm{o}\mathrm{s}\left({k}_{z}h\right)\left({b}_{1}\mathrm{s}\mathrm{i}\mathrm{n}\frac{{k}_{x}h}{2}+{b}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\frac{3{k}_{x}h}{2}\right)\right\}-\\ \;\;\;\;\frac{\kappa {A}_{{\upsilon }_{z}}}{h}\left\{\sum\limits _{m=1}^{M}{a}_{m}\mathrm{s}\mathrm{i}\mathrm{n}\left[\left(m-\frac{1}{2}\right){k}_{z}h\right]+\right.\\ \;\;\;\;2\mathrm{c}\mathrm{o}\mathrm{s}\left({k}_{x}h\right)\left.\left({b}_{1}\mathrm{s}\mathrm{i}\mathrm{n}\frac{{k}_{z}h}{2}+{b}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\frac{3{k}_{z}h}{2}\right)\right\}\\ \frac{{A}_{{\upsilon }_{x}}}{\mathrm{\Delta }t}\mathrm{s}\mathrm{i}\mathrm{n}\frac{\omega \mathrm{\Delta }t}{2}\\ \approx -\frac{{A}_{P}}{\rho h}\left\{\sum\limits _{m=1}^{M}{a}_{m}\mathrm{s}\mathrm{i}\mathrm{n}\left[\left(m-\frac{1}{2}\right){k}_{x}h\right]+\right.\\ \;\;\;\;2\mathrm{c}\mathrm{o}\mathrm{s}\left({k}_{z}h\right)\left.\left({b}_{1}\mathrm{s}\mathrm{i}\mathrm{n}\frac{{k}_{x}h}{2}+{b}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\frac{3{k}_{x}h}{2}\right)\right\}\\ \frac{{A}_{{\upsilon }_{z}}}{\mathrm{\Delta }t}\mathrm{s}\mathrm{i}\mathrm{n}\frac{\omega \mathrm{\Delta }t}{2}\\ \approx -\frac{{A}_{P}}{\rho h}\left\{\sum\limits _{m=1}^{M}{a}_{m}\mathrm{s}\mathrm{i}\mathrm{n}\left[\left(m-\frac{1}{2}\right){k}_{z}h\right]+\right.\\ \;\;\;\;\left.2\mathrm{c}\mathrm{o}\mathrm{s}\left({k}_{x}h\right)\left({b}_{1}\mathrm{s}\mathrm{i}\mathrm{n}\frac{{k}_{z}h}{2}+{b}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\frac{3{k}_{z}h}{2}\right)\right\}\end{array}\right. $

(6)

消去式(6)中的$ {A}_{P} $，$ {A}_{{\upsilon }_{x}} $和$ {A}_{{\upsilon }_{z}} $，并且考虑到$ \omega =vk $和$ \kappa =\rho {v}^{2} $，得到

$ \begin{array}{l}\frac{1}{{\left(v\mathrm{\Delta }t\right)}^{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\frac{rkh}{2}\approx \frac{1}{{h}^{2}}\left\{\sum\limits _{m=1}^{M}{a}_{m}\left[\mathrm{s}\mathrm{i}\mathrm{n}\left(m-\frac{1}{2}\right){k}_{x}h\right]+\right.\\ {\left.2\mathrm{c}\mathrm{o}\mathrm{s}\left({k}_{z}h\right)\left({b}_{1}\mathrm{s}\mathrm{i}\mathrm{n}\frac{{k}_{x}h}{2}+{b}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\frac{3{k}_{x}h}{2}\right)\right\}}^{2}+\\ \frac{1}{{h}^{2}}\left\{\sum\limits _{m=1}^{M}{a}_{m}\left[\mathrm{s}\mathrm{i}\mathrm{n}\left(m-\frac{1}{2}\right){k}_{z}h\right]+\right.\\ {\left.2\mathrm{c}\mathrm{o}\mathrm{s}\left({k}_{x}h\right)\left({b}_{1}\mathrm{s}\mathrm{i}\mathrm{n}\frac{{k}_{z}h}{2}+{b}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\frac{3{k}_{z}h}{2}\right)\right\}}^{2}\end{array} $

(7)

式中：$ v $为介质中声波的传播速度；$ r=v\mathrm{\Delta }t/h $为Courant条件数。

式(7)为M-SFD(N=2)给出的差分离散声波方程的频散关系，也称为时空域频散关系。泰勒展开其中的三角函数，得到

$ \begin{array}{l}{\left\{\sum\limits _{j=0}^{\mathrm{\infty }}{c}_{j}{\beta }_{j}{\left({k}_{x}/2\right)}^{2j+1}{h}^{2j}+2\left(\sum\limits _{j=1}^{\mathrm{\infty }}{\gamma }_{j}{k}_{z}^{2j}{h}^{2j}\right)\left[{b}_{1}\sum\limits _{j=0}^{\mathrm{\infty }}{\beta }_{j}{\left({k}_{x}/2\right)}^{2j+1}{h}^{2j}+{b}_{2}\sum\limits _{j=0}^{\mathrm{\infty }}{\beta }_{j}{\left(3{k}_{x}/2\right)}^{2j+1}{h}^{2j}\right]\right\}}^{2}+\\ {\left\{\sum\limits _{j=0}^{\mathrm{\infty }}{c}_{j}{\beta }_{j}{\left({k}_{z}/2\right)}^{2j+1}{h}^{2j}+2\left(\sum\limits _{j=1}^{\mathrm{\infty }}{\gamma }_{j}{k}_{x}^{2j}{h}^{2j}\right)\left[{b}_{1}\sum\limits _{j=0}^{\mathrm{\infty }}{\beta }_{j}{\left({k}_{z}/2\right)}^{2j+1}{h}^{2j}+{b}_{2}\sum\limits _{j=1}^{\mathrm{\infty }}{\beta }_{j}{\left(3{k}_{z}/2\right)}^{2j+1}{h}^{2j}\right]\right\}}^{2}\\ \approx {\left[\sum\limits _{j=0}^{\mathrm{\infty }}{r}^{2j}{\beta }_{j}{\left(k/2\right)}^{2j+1}{h}^{2j}\right]}^{2}\end{array} $

(8)

其中$ {c}_{j} $、$ {\beta }_{j} $和$ {\gamma }_{j} $的表达式为

$ \left\{\begin{array}{l}{c}_{j}=\sum\limits _{m=1}^{M}{\left(2m-1\right)}^{2j+1}{a}_{m}+2{b}_{1}+{3}^{2j+1}·2{b}_{2}\\ {\beta }_{j}=\frac{{\left(-1\right)}^{j}}{\left(2j+1\right)!}\\ {\gamma }_{j}=\frac{{\left(-1\right)}^{j}}{\left(2j\right)!}\end{array}\right. $

(9)

使式(8)左右两边$ {k}_{x}^{2}{k}_{z}^{2}{h}^{2} $的系数对应相等，得到

$ {c}_{0}\left({b}_{1}+3{b}_{2}\right)=\frac{{r}^{2}}{24} $

(10)

使式(8)左右两边$ {k}_{x}^{2}{k}_{z}^{4}{h}^{4} $(或$ {k}_{x}^{4}{k}_{z}^{2}{h}^{4} $)的系数对应相等，得到

$ \begin{array}{l}{c}_{0}\left(3{b}_{1}+33{b}_{2}\right)+{c}_{1}\left({b}_{1}+3{b}_{2}\right)+\\ \;\;\;\;\;\;12({b}_{1}+3{b}_{2}{)}^{2}=\frac{{r}^{4}}{10}\end{array} $

(11)

使式(8)左右两边$ {k}_{x}^{2j+2}{h}^{2j} $(或$ {k}_{z}^{2j+2}{h}^{2j} $)$ \left(j=\mathrm{0, 1}, 2, \dots , M-1\right) $的系数对应相等，得到

$ \left\{\begin{array}{cc}\begin{array}{l}\begin{array}{cc}{c}_{0}^{2}=1& j=0\end{array}\\ \sum\limits _{p=0}^{j}{c}_{p}{c}_{j-p}{\beta }_{p}{\beta }_{j-p}\\ =\sum\limits _{p=0}^{j}{\beta }_{p}{\beta }_{j-p}{r}^{2j}j=\mathrm{1, 2}, \cdots , M-1\end{array}& \end{array}\right. $

(12)

根据式(12)可得出$ {c}_{0}=\pm 1 $，当$ {c}_{0} $从$ 1 $变为$ -1 $，相应的差分系数为$ {a}_{1}, {a}_{2}, \cdots , {a}_{M};{b}_{1} $变为其相反数，对最终结果没有影响。这里取$ {c}_{0}=1 $，并进行计算和推导可以得到

$ \begin{array}{cc}{c}_{j}={r}^{2j}& j=\mathrm{0, 1}, \cdots , M-1\end{array} $

(13)

由式(13)可知$ {c}_{0}=1 $，$ {c}_{1}={r}^{2} $，将它们代入式(10)和式(11)并联立求解，可得到

$ \begin{array}{cc}\left\{\begin{array}{l}{b}_{1}=-\frac{3{r}^{4}}{640}+\frac{11{r}^{2}}{192}\\ {b}_{2}=\frac{{r}^{4}}{640}-\frac{{r}^{2}}{192}\end{array}\right.& \end{array} $

(14)

将式(13)代入式(9)得到

$ \left\{\begin{array}{l}\sum\limits _{m=1}^{M}{\left(2m-1\right)}^{2j+1}{a}_{m}+2{b}_{1}+{3}^{2j+1}2{b}_{2}={r}^{2j}\\ j=\mathrm{0, 1}, \cdots , M-1\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\left(15\right)\end{array}\right. $

(15)

式(15)可改写为如下矩阵方程

$ \left[\begin{array}{ccccc}{1}^{0}& {3}^{0}& {5}^{0}& {(2M-1)}^{0}& 1\\ {1}^{2}& {3}^{2}& {5}^{2}& \cdots & {\left(2M-1\right)}^{2}\\ {1}^{4}& {3}^{4}& {5}^{4}& \cdots & {\left(2M-1\right)}^{4}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {1}^{2M-2}& {3}^{2M-2}& {5}^{2M-2}& \cdots & {\left(2M-1\right)}^{2M-2}\end{array}\right]\left[\begin{array}{c}{a}_{1}+2{b}_{1}\\ 3\left({a}_{2}+2{b}_{2}\right)\\ 5{a}_{3}\\ ⋮\\ \left(2M-1\right){a}_{M}\end{array}\right]=\left[\begin{array}{c}1\\ {r}^{2}\\ {r}^{4}\\ ⋮\\ {r}^{2M-2}\end{array}\right] $

(16)

式(16)是一个范德蒙矩阵方程。

联合求解方程(14)和(16)得到

$ \left\{\begin{aligned} b_1= & -\frac{3 r^4}{640}+\frac{11 r^2}{192} \\ b_2= & \frac{r^4}{640}-\frac{r^2}{192} \\ a_1= & \prod\nolimits_{2 \leqslant k \leqslant M}\left[\frac{r^2-(2 k-1)^2}{1-(2 k-1)^2}\right]-2 b_1 \\ a_2= & \frac{1}{3} \prod\nolimits_{1 \leqslant k \leqslant M, k \neq 2}\left[\frac{r^2-(2 k-1)^2}{3^2-(2 k-1)^2}\right]-2 b_2 \\ a_m= & \frac{1}{2 m-1} \prod\nolimits_{1 \leqslant k \leqslant M, k \neq m} \frac{r^2-(2 k-1)^2}{(2 m-1)^2-(2 k-1)^2} \\ & m=3,4, \cdots, M \end{aligned}\right.$

(17)

式(17)为M-SFD(N=2)的差分系数通解，利用同样的方法，可以导出M-SFD(N=1、3、4)的差分系数通解，见附录A。

1.3 差分精度分析

可以将差分离散声波方程的频散关系误差函数关于空间采样间隔$ h $或时间采样间隔$ \mathrm{\Delta }t $的最小幂指数定义为差分离散声波方程的差分阶数。根据M-SFD(N=2)给出的差分离散声波方程的频散关系式(7)，定义误差函数$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=2\right)} $，表达式为

$ \begin{array}{l}{E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=2\right)}=\frac{1}{{h}^{2}}\left\{\sum\limits _{m=1}^{M}{a}_{m}\mathrm{s}\mathrm{i}\mathrm{n}\left[\left(m-\frac{1}{2}\right){k}_{x}h\right]+\right.\\ {\left.2\mathrm{c}\mathrm{o}\mathrm{s}\left({k}_{z}h\right)\left({b}_{1}\mathrm{s}\mathrm{i}\mathrm{n}\frac{{k}_{x}h}{2}+{b}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\frac{3{k}_{x}h}{2}\right)\right\}}^{2}+\\ \frac{1}{{h}^{2}}\left\{\sum\limits _{m=1}^{M}{a}_{m}\mathrm{s}\mathrm{i}\mathrm{n}\left[\left(m-\frac{1}{2}\right){k}_{z}h\right]+\right.\\ {\left.2\mathrm{c}\mathrm{o}\mathrm{s}\left({k}_{x}h\right)\left({b}_{1}\mathrm{s}\mathrm{i}\mathrm{n}\frac{{k}_{z}h}{2}+{b}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\frac{3{k}_{z}h}{2}\right)\right\}}^{2}-\\ \frac{1}{{\left(v\mathrm{\Delta }t\right)}^{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\frac{rkh}{2}\left(18\right)\end{array} $

(18)

结合式(8)~式(12)，$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=2\right)} $可化简为

$ \begin{array}{l}{E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=2\right)}=\sum\limits _{j=M}^{\mathrm{\infty }}\sum\limits _{p=0}^{j}\left({c}_{p}{c}_{j-p}-{r}^{2j}\right){\beta }_{p}{\beta }_{j-p}\frac{1}{{2}^{2j+2}}\left({k}_{x}^{2j+2}+{k}_{z}^{2j+2}\right){h}^{2j}+\\ \sum\limits _{j=3}^{\mathrm{\infty }}\sum\limits _{p=0}^{j-1}\left\{\frac{{\gamma }_{j-p}}{{2}^{2p+2}}\sum\limits _{q=0}^{p}2{c}_{q}\left[{b}_{1}+{3}^{2\left(p-q\right)+1}{b}_{2}\right]{\beta }_{q}{\beta }_{p-q}+\frac{{\gamma }_{p+1}}{{2}^{2\left(j-p\right)}}\sum\limits _{q=0}^{j-p-1}2{c}_{q}\left[{b}_{1}+{3}^{2\left(j-p-1-q\right)+1}{b}_{2}\right]{\beta }_{q}{\beta }_{j-p-1-q}-\right.\\ \left.\frac{{r}^{2j}{C}_{j+1}^{p+1}}{{2}^{2j+2}}\sum\limits _{q=0}^{j}{\beta }_{q}{\beta }_{j-q}\right\}{k}_{x}^{2p+2}{k}_{z}^{2\left(j-p\right)}{h}^{2j}\left(19\right)\end{array} $

(19)

其中$ {C}_{j+1}^{p+1} $为组合数；$ {c}_{j} $，$ {\beta }_{j} $和$ {\gamma }_{j} $的表达式见式(9)。

式(19)表明，$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=2\right)} $关于空间采样间隔$ h $的最小幂指数为6，因此M-SFD(N=2)给出的差分离散声波方程具有6阶差分精度。式(19)给出的$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=2\right)} $的表达式较为复杂，还可以从差分系数求解过程直接分析M-SFD(N=2)给出的差分离散声波方程的差分精度。根据M-SFD(N=2)的差分系数求解过程可以看出，$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=2\right)} $可以看成关于$ h $的多项式，差分系数会使$ {E}_{\mathrm{M}-\mathrm{S}\mathrm{F}\mathrm{D}\left(N=2\right)} $中$ {k}_{x}^{2j+2}{h}^{2j} $，$ {k}_{z}^{2j+2}{h}^{2j} $$ \left(j=\mathrm{0, 1}, 2, \dots , M-1\right) $、$ {k}_{x}^{2}{k}_{z}^{2}{h}^{2} $、$ {k}_{x}^{2}{k}_{z}^{4}{h}^{4} $和$ {k}_{x}^{4}{k}_{z}^{2}{h}^{4} $的系数均为零，$ {E}_{\mathrm{M}-\mathrm{S}\mathrm{F}\mathrm{D}\left(N=2\right)} $中系数不为零的最低次项为$ {k}_{x}^{2}{k}_{z}^{6}{h}^{6} $、$ {k}_{x}^{6}{k}_{z}^{2}{h}^{6} $和$ {k}_{x}^{4}{k}_{z}^{4}{h}^{6} $，因此，M-SFD(N=2)给出的差分离散声波方程具有6阶差分精度。

同样地，M-SFD(N=1)给出的差分离散声波方程的频散关系误差函数$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=1\right)} $可以看成关于$ h $的多项式，差分系数能够使得$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=1\right)} $中$ {k}_{x}^{2j+2}{h}^{2j} $、$ {k}_{z}^{2j+2}{h}^{2j} $$ \left(j=\mathrm{0, 1}, 2, \dots , M-1\right) $、和$ {k}_{x}^{2}{k}_{z}^{2}{h}^{2} $的系数为零，$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=1\right)} $中系数不为零的最低次项为$ {k}_{x}^{2}{k}_{z}^{4}{h}^{4} $和$ {k}_{x}^{4}{k}_{z}^{2}{h}^{4} $，因此M-SFD(N=1)给出的差分离散声波方程具有4阶差分精度；M-SFD(N=3)给出的差分离散声波方程的频散关系误差函数$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=3\right)} $可以看成关于$ h $的多项式，差分系数会使得$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=3\right)} $中$ {k}_{x}^{2j+2}{h}^{2j} $、$ {k}_{z}^{2j+2}{h}^{2j} $$ \left(j=\mathrm{0, 1}, 2, \dots , M-1\right) $、$ {k}_{x}^{2}{k}_{z}^{2}{h}^{2} $、$ {k}_{x}^{2}{k}_{z}^{4}{h}^{4} $、$ {k}_{x}^{4}{k}_{z}^{2}{h}^{4} $和$ {k}_{x}^{4}{k}_{z}^{4}{h}^{6} $的系数均为零，$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=3\right)} $中系数不为零的最低次项为$ {k}_{x}^{2}{k}_{z}^{6}{h}^{6} $和$ {k}_{x}^{6}{k}_{z}^{2}{h}^{6} $，因此M-SFD(N=3)给出的差分离散声波方程具有6阶差分精度；M-SFD(N=4)给出的差分离散声波方程的频散关系误差函数$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=4\right)} $可以看成关于$ h $的多项式，差分系数会使得$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=4\right)} $中$ {k}_{x}^{2j+2}{h}^{2j} $、$ {k}_{z}^{2j+2}{h}^{2j} $$ \left(j=\mathrm{0, 1}, 2, \dots , M-1\right) $，$ {k}_{x}^{2}{k}_{z}^{2}{h}^{2} $、$ {k}_{x}^{2}{k}_{z}^{4}{h}^{4} $、$ {k}_{x}^{4}{k}_{z}^{2}{h}^{4} $、$ {k}_{x}^{4}{k}_{z}^{4}{h}^{6} $、$ {k}_{x}^{2}{k}_{z}^{6}{h}^{6} $和$ {k}_{x}^{6}{k}_{z}^{2}{h}^{6} $的系数均为零，$ {E}_{\mathrm{M}\mathrm{-}\mathrm{S}\mathrm{F}\mathrm{D}\left(N=4\right)} $中系数不为零的关于$ h $的最低次项为$ {k}_{x}^{2}{k}_{z}^{8}{h}^{8} $、$ {k}_{x}^{8}{k}_{z}^{2}{h}^{8} $、$ {k}_{x}^{4}{k}_{z}^{6}{h}^{8} $和$ {k}_{x}^{6}{k}_{z}^{4}{h}^{8} $，因此M-SFD(N=4)给出的差分离散声波方程具有8阶差分精度。理论上，继续增大N的取值，M-SFD给出的差分离散声波方程可以达到任意偶数阶差分精度。C-SFD给出的差分离散声波方程仅具有2阶差分精度^[21]，因此，M-SFD能更有效地提高数值模拟精度。

表 1给出了当差分离散声波方程达到相同阶的差分精度时，本文的M-SFD与Tan等^[33]构建的时间和空间高阶交错网格差分格式需要的非坐标网格点数。对比可以看出，达到相同阶的差分精度，本文的M-SFD需要的非坐标轴网格点数更少，因而计算量更小，计算效率更高。

2 频散分析和稳定性分析 2.1 频散分析

数值频散能直接反映交错网格有限差分法的数值模拟精度。本文引入归一化相速度误差函数$ {\varepsilon }_{ph}\left(kh, \theta \right) $定量描述数值频散的大小。根据M-SFD(N=1)给出的差分离散声波方程的频散关系式(7)和相速度的定义$ {v}_{ph}=\omega /k $，可导出$ {\varepsilon }_{ph}\left(kh, \theta \right) $的表达式为

$ \left\{\begin{array}{l}{\varepsilon }_{ph}\left(kh, \theta \right)=\frac{{v}_{ph}}{v}-1=\frac{2}{rkh}\mathrm{s}\mathrm{i}{\mathrm{n}}^{-1}\left(r\sqrt[]{q}\right)-1\\ q=\left\{\sum\limits _{m=1}^{M}{a}_{m}\mathrm{s}\mathrm{i}\mathrm{n}\left[\left(m-\frac{1}{2}\right)kh\mathrm{c}\mathrm{o}\mathrm{s}\theta \right]\right.+\\ \;\;\;\;\;\;2{b}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\left(kh\mathrm{s}\mathrm{i}\mathrm{n}\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}{\frac{kh\mathrm{c}\mathrm{o}\mathrm{s}\theta }{2}}^{2}+\\\;\;\;\;\;\; \left\{\sum\limits _{m=1}^{M}{a}_{m}\mathrm{s}\mathrm{i}\mathrm{n}\left[\left(m-\frac{1}{2}\right)kh\mathrm{s}\mathrm{i}\mathrm{n}\theta \right]\right.+\\\;\;\;\;\;\;2{b}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\left(kh\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}{\frac{kh\mathrm{s}\mathrm{i}\mathrm{n}\theta }{2}}^{2}\end{array}\right. $

(20)

$ {\varepsilon }_{ph}\left(kh, \theta \right)=0 $时，相速度与真实速度相等，无数值频散；$ {\varepsilon }_{ph}\left(kh, \theta \right) > 0\mathrm{时} $，相速度偏大，呈现时间频散；$ {\varepsilon }_{ph}\left(kh, \theta \right) < 0 $时，相速度偏小，呈现空间频散。同理，可导出C-SFD和M-SFD(N=2、3、4)的归一化相速度误差$ {\varepsilon }_{ph}\left(kh, \theta \right) $的表达式。

图 3给出了不同M值时C-SFD和M-SFD的归一化相速度误差$ {\varepsilon }_{ph}\left(kh, \theta \right) $随$ kh $和$ \theta $变化的特征。可以看出：①M=2时，C-SFD呈现明显的空间频散，相速度误差达到-4.0%(图 3a左)；M增至5时，空间频散明显减小，相速度误差最大为-1.5%，但出现较明显的时间频散，相速度误差最大为2.0%(图 3b左)；M增至8，空间频散基本消失，时间频散进一步增大，相速度误差最大为2.0%(图 3c左)。②M=2时，M-SFD呈现明显的空间频散，相速度误差达到-4.0%(图 3a右)；M增至5时，空间频散明显减小，相速度误差最大为-3.0%(图 3b右)；M增至8时，空间频散进一步减小，相速度误差最大为-1.0%(图 3c右)。③M取值较小时(如M=2)，M-SFD和C-SFD均具有较严重的数值频散，模拟精度均较低；M取值较大时(如M=8)，M-SFD的数值频散明显小于C-SFD，因此M-SFD的模拟精度优于C-SFD。④如果将数值频散误差的绝对值小于1‰定义为高精度。M的取值从2增至8时，C-SFD的高精度模拟区域基本没有增大，而M-SFD的高精度模拟区域逐步增大。M取值大于5时，M-SFD的高精度模拟区域明显大于C-SFD。⑤M的取值从2增至8时，C-SFD和M-SFD的空间频散均逐步减小，说明增加空间差分算子中坐标轴网格点数有助于减小空间频散；M取值相同时，M-SFD的时间频散均小于C-SFD，说明增加空间差分算子中非坐标轴网格点数有助于减小时间频散。

图 4给出了不同M及N时，M-SFD的归一化相速度误差$ {\varepsilon }_{ph}\left(kh, \theta \right) $随$ kh $和$ \theta $变化的情况。需要注意M-SFD(M=10；N=1、2、3、4)和M-SFD(M=20；N=1、2、3、4)两组子图中等值线的刻度不同，色标代表的相速度误差范围也不同。可以看出：①M=10，N=1、2、3、4时，M-SFD的数值频散特征相似，基本能将相速度误差的绝对值控制在1‰以内(图 4左)。②M=20，N=1时，M-SFD主要表现出时间频散，相速度误差最大为3.0‰(图 4a右)；N增至2时，时间频散显著减小(图 4b右)；N增至4时，时间频散进一步减小(图 4d右)；M=20，N=4时，M-SFD基本能将相速度误差的绝对值控制在0.1‰以内(图 4d右)。③对于一般的高精度模拟，要求将相速度误差的绝对值控制在1‰以内，建议采用M-SFD(N=1)，M的取值为10左右；如果对模拟精度要求极高，要求将相速度误差的绝对值控制在0.1‰以内，可以采用M-SFD(N=4)，M的取值约为20。因此，M-SFD应根据模拟精度要求，合理选择M和N的取值以兼顾计算效率。④固定M的取值，N的取值为2和3时，M-SFD的相速度误差分布规律基本一致，这是因为此时M-SFD给出的差分离散声波方程均为6阶差分精度，见表 1.

2.2 稳定性分析

根据M-SFD(N=2)给出的差分离散声波方程的频散关系式(7)得到

$ \begin{aligned} & \frac{1}{r^2} \sin ^2 \frac{r k h}{2} \approx \frac{1}{h^2}\left\{\sum\limits_{m=1}^M a_m \sin \left[\left(m-\frac{1}{2}\right) k_x h\right]+\right. \\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left.2 \cos \left(k_z h\right)\left(b_1 \sin \frac{k_x h}{2}+b_2 \sin \frac{3 k_x h}{2}\right)\right\}+ \\ & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{h^2} \sum\limits_{m=1}^M a_m \sin \left[\left(m-\frac{1}{2}\right) k_z h\right]+ \\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left.2 \cos \left(k_x h\right)\left(b_1 \sin \frac{k_z h}{2}+b_2 \sin \frac{3 k_z h}{2}\right)\right\} \end{aligned} $

(21)

取空间波数$ {k}_{x}={k}_{z}=\mathrm{\pi }/h $，并且考虑到$ 0\le \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left(rkh/2\right)\le 1 $，得

$ r\le S=\frac{1}{\sqrt[]{2}\left|\sum\limits _{m=1}^{M}{\left(-1\right)}^{m-1}{a}_{m}-2{b}_{1}+2{b}_{2}\right|} $

(22)

式中$ S $为稳定因子。

式(22)为M-SFD(N=2)给出的差分离散声波方程的稳定性条件表达式。同样地，可以导出C-SFD和M-SFD(N=1、3、4)的稳定性条件表达式。图 5给出了稳定性条件约束下，$ r $的最大取值随M的变化曲线，称为稳定性曲线。可以看出：①C-SFD和M-SFD(N=1、2、3，4)的稳定性均随M取值的增大而降低；②M取值相同时，M-SFD(N=1、2、3，4)的稳定性强于C-SFD；③固定M取值，M-SFD的稳定性会随N取值的增大而增强，因此，增加空间差分算子中非坐标轴网格点数能增强稳定性；④N的取值为2和3时，M-SFD的稳定性完全相同，这是因为M-SFD(N=2、3)给出的差分离散声波方程均为6阶差分精度，见表 1。

3 数值模拟实例 3.1 层状介质模型

图 6给出了一个包含5个反射界面的层状介质速度模型，模型尺寸为$ 12\text{ }\mathrm{k}\mathrm{m}\times 12\text{ }\mathrm{k}\mathrm{m} $，模型空间采样间隔$ h=15\text{ }\mathrm{m} $，网格数为$ 801\times 801 $。震源采用主频为22 Hz的Ricker子波，位于(0.15 km，0.15 km)。C-SFD(M=12)采用时间采样间隔$ \mathrm{\Delta }t=0.5 $、$ 1.0\text{ }\mathrm{m}\mathrm{s} $；M-SFD(M=10；N=1)采用$ \mathrm{\Delta }t=1.5 $、$ 2.0\text{ }\mathrm{m}\mathrm{s} $；M-SFD(M=10；N=2，3，4)采用$ \mathrm{\Delta }t=2.0\text{ }\mathrm{m}\mathrm{s} $进行模拟。图 7和图 8分别给出了C-SFD和M-SFD(M=10；N=1)模拟得到的3.6 s时刻压力场$ P $和质点振动速度场$ x $分量$ {\upsilon }_{x} $的波场快照。图 9a给出了M-SFD(M=10；N=2)模拟生成的单炮记录。为对比方便，图 9b和图 9c分别为C-SFD(M=12)和M-SFD(M=10；N=1，2，3，4)采用不同时间采样间隔模拟生成的炮集记录中第781道不同时间区间范围内的波形对比。

表 2统计了C-SFD(M=12)和M-SFD(M=10；N=1，2，3，4)采用不同时间采样间隔进行模拟的耗时和加速比，模拟总时长为9.0 s，且均采用Inter Xeon Gold 5218 CPU处理器。

图 7和图 8中两组波场快照显示：C-SFD(M=12)采用$ \mathrm{\Delta }t=0.5\text{ }\mathrm{m}\mathrm{s} $进行模拟的波场快照中无明显数值频散，当时间采样间隔增至$ \mathrm{\Delta }t=1.0\text{ }\mathrm{m}\mathrm{s} $时，波场快照中出现明显的时间数值频散；M-SFD(M=10；N=1，2)采用$ \mathrm{\Delta }t=2.0\text{ }\mathrm{m}\mathrm{s} $进行模拟得到的波场快照均无明显数值频散。

图 9b和图 9c中的单道波形显示：C-SFD(M=12)采用$ \mathrm{\Delta }t=1.0\text{ }\mathrm{m}\mathrm{s} $进行模拟，波形畸变严重，采用$ \mathrm{\Delta }t=0.5\text{ }\mathrm{m}\mathrm{s} $进行模拟，波形仍存在一定程度的畸变；M-SFD(M=10；N=1)采用$ \mathrm{\Delta }t=2.0\text{ }\mathrm{m}\mathrm{s} $进行模拟，波形存在一定程度的畸变，采用$ \mathrm{\Delta }t=1.5\text{ }\mathrm{m}\mathrm{s} $进行模拟，波形仅存在微小畸变；M-SFD(M=10；N=2，3，4)采用$ \mathrm{\Delta }t=2.0\text{ }\mathrm{m}\mathrm{s} $进行模拟，波形畸变基本可以忽略，模拟波形与参考波形基本重合。

综合图 7~图 9给出的数值频散特征及表 2给出的计算效率统计结果，可以看出：M-SFD(M=10；N=1)采用$ \mathrm{\Delta }t=1.5\text{ }\mathrm{m}\mathrm{s} $比C-SFD(M=12)采用$ \mathrm{\Delta }t=0.5\text{ }\mathrm{m}\mathrm{s} $模拟的数值频散更小，精度更高，且计算效率提高了1.96倍；M-SFD(M=10；N=2)采用$ \mathrm{\Delta }t=2.0\text{ }\mathrm{m}\mathrm{s} $比M-SFD(M=10；N=1)采用$ \mathrm{\Delta }t=1.5\text{ }\mathrm{m}\mathrm{s} $模拟，数值频散更小，精度更高，且计算效率更高。

仔细对比图 9c中的单道波形会发现，采用相同的时间采样间隔$ \mathrm{\Delta }t=2.0\text{ }\mathrm{m}\mathrm{s} $，M-SFD(M=10；N=3)与M-SFD(M=10；N=2)的模拟精度基本一致，M-SFD(M=10；N=4)相比M-SFD(M=10；N=2，3)在模拟精度方面存在十分微弱的优势，需要在更长传播距离和传播时间的情况下，这种优势才能更好地体现出来。

3.2 复杂构造模型

图 10a为中国塔里木盆地典型复杂构造速度模型，模型尺寸为$ 18\text{ }\mathrm{k}\mathrm{m}\times 7.875\text{ }\mathrm{k}\mathrm{m} $，模型空间采样间隔$ h=15\text{ }\mathrm{m} $，网格数为$ 1201\times 526 $。震源采用主频为25 Hz的Ricker子波，位于(9 km，0.15 km)，C-SFD(M=10)和M-SFD(M=8；N=1)分别采用时间采样间隔$ \mathrm{\Delta }t=1.0\text{ }\mathrm{m}\mathrm{s} $和$ \mathrm{\Delta }t=1.5\text{ }\mathrm{m}\mathrm{s} $进行模拟。C-SFD(M=10)和M-SFD(M=8；N=1)的空间差分算子均包含20个网格点，单次时间迭代的计算量基本相同，模拟时长相同时，计算效率与时间采样间隔成正比。

图 10b给出了M-SFD(M=8；N=1)采用采$ \mathrm{\Delta }t=1.5\text{ }\mathrm{m}\mathrm{s} $模拟生成的单炮记录。为了对比方便，图 10c~图 10f给出了C-SFD(M=10)和M-SFD(M=8；N=1)不同时间采样间隔模拟生成单炮的局部放大显示。对比可以看出：C-SFD(M=10)采用$ \mathrm{\Delta }t=1.0\text{ }\mathrm{m}\mathrm{s} $和$ \mathrm{\Delta }t=1.5\text{ }\mathrm{m}\mathrm{s} $模拟生成的单炮记录中均存在明显的时间频散；M-SFD(M=8；N=1)采用两种$ \mathrm{\Delta }t $模拟生成的单炮记录中均无明显数值频散。

复杂构造模型模拟结果表明：计算效率基本相同时，M-SFD比C-SFD能更有效地压制数值频散，模拟精度更高；M-SFD可采用更大的时间采样间隔以提高计算效率，并保持更高的模拟精度。

4 混合交错网格有限差分法在逆时偏移中的应用

将M-SFD作为逆时偏移算法中的波场传播算子，以图 10a中的复杂构造模型进行逆时偏移测试。震源采用主频为25 Hz的Ricker子波。C-SFD(M=15)采用非常小的时间采样间隔$ \mathrm{\Delta }t=0.1\text{ }\mathrm{m}\mathrm{s} $，模拟生成150炮无数值频散的炮集资料作为逆时偏移的输入道集，每炮600道接收，炮间距120 m，道间距30 m。

分别以C-SFD(M=10)和M-SFD(M=8；N=1)作为逆时偏移算法中的波场传播算子，时间采样间隔$ \mathrm{\Delta }t=1.5\text{ }\mathrm{m}\mathrm{s} $，采用互相关成像条件，并对成像结果做Laplace滤波去除低频噪声。图 11给出了C-SFD(M=10)和M-SFD(M=8；N=1)的逆时偏移剖面。对比可以看出：C-SFD(M=10)的偏移成像结果中，深层同相轴存在大量由于数值频散造成的成像假象；M-SFD(M=8；N=1)的偏移成像结果中，由于数值频散造成的成像假象消失，深层同相轴能量更强，分辨率更高。因此，相比C-SFD，M-SFD作为逆时偏移的波场传播算子能有效提高深层构造的成像精度和分辨率。

5 结论

本文联合利用坐标轴网格点和非坐标轴网格点构建空间差分算子近似波场的一阶空间偏导数，建立了适用于速度—应力声波方程模拟的M-SFD，并基于时空域频散关系和泰勒级数展开导出了差分系数通解。通过差分精度分析、频散分析、稳定性分析、数值模拟和逆时偏移试验，得到以下结论：

(1) C-SFD给出的差分离散声波方程仅具有二阶差分精度，M-SFD(N=1、2、4)给出的差分离散声波方程可达到4、6、8阶差分精度，继续增大N的取值，理论上可以达到任意偶数阶差分精度；

(2) 计算效率基本相同时，M-SFD比C-SFD能更有效地压制数值频散，模拟精度更高；M-SFD可通过采用更大的时间采样间隔以提高计算效率，且模拟精度仍然高于C-SFD；

(3) M-SFD作为逆时偏移的波场传播算子能够更有效地消除数值频散造成的成像假象，进而提高深层的成像精度和分辨率。

附录A M-SFD(N=1，3，4)的差分离散声波方程及差分系数通解

与文中M-SFD(N=2)的差分系数求解过程类似，利用M-SFD(N=1)对声波方程(式(1))进行差分离散可以导出相应差分离散声波方程，再将离散平面波解代入差分离散声波方程导出相应的时空域频散关系，并利用泰勒公式对频散关系中的三角函数进行级数展开，然后，使得方程左右两边$ {k}_{x}^{2}{k}_{z}^{2}{h}^{2} $和$ {k}_{x}^{2j+2}{h}^{2j} $(或$ {k}_{z}^{2j+2}{h}^{2j} $)(j=0，1，2，…，M-1)的系数对应相等以构建关于差分系数的方程组，求解此方程组可以得到M-SFD(N=1)的差分系数a₁，a₂，…，a_M；b₁的通解为

$ \left\{\begin{array}{l} b_1=\frac{r^2}{24} \\ a_1=\prod\limits_{2 \leqslant k \leqslant M} \frac{r^2-(2 k-1)^2}{1-(2 k-1)^2}-\frac{r^2}{12} \\ a_m=\frac{1}{2 m-1} \prod\limits_{1 \leqslant k \leqslant M, k \neq m} \frac{r^2-(2 k-1)^2}{(2 m-1)^2-(2 k-1)^2} \\ \quad\quad\quad m=2,3, \cdots, M \end{array}\right. $

(A-1)

利用泰勒公式对M-SFD(N=3)给出的差分离散声波方程的时空域频散关系中的三角函数进行级数展开，然后，使得方程左右两边$ {k}_{x}^{2}{k}_{z}^{2}{h}^{2} $，$ {k}_{x}^{2}{k}_{z}^{4}{h}^{4} $(或$ {k}_{x}^{4}{k}_{z}^{2}{h}^{4} $)，$ {k}_{x}^{4}{k}_{z}^{4}{h}^{6} $和$ {k}_{x}^{2j+2}{h}^{2j} $(或$ {k}_{z}^{2j+2}{h}^{2j} $)$ \left(j=\mathrm{0, 1}, 2, \dots , M-1\right) $的系数对应相等以构建关于差分系数的方程组，求解此方程组可以得到M-SFD(N=3)的差分系数$ {a}_{1}, {a}_{2}, \dots , {a}_{M};{b}_{1}, {b}_{2}, {b}_{3} $的通解为

$ \left\{\begin{array}{l} b_1=\frac{r^6}{4480}-\frac{5 r^4}{128}+\frac{37 r^2}{576} \\ b_2=\frac{r^6}{4480}-\frac{r^4}{640}+\frac{r^2}{576} \\ b_3=-\frac{r^6}{4480}+\frac{r^4}{320}-\frac{r^2}{144} \\ a_1=\prod\limits_{2 \leqslant k \leqslant M} \frac{r^2-(2 k-1)^2}{1-(2 k-1)^2}-2 b_1-2 b_3 \\ a_2=\frac{1}{3} \prod\limits_{1 \leqslant k \leqslant M} \frac{r^2-(2 k-1)^2}{3^2-(2 k-1)^2}-2 b_2 \\ a_m=\frac{1}{2 m-1} \prod\limits_{\substack{1 \leqslant k \leqslant M \\ k \neq m}} \frac{r^2-(2 m-1)^2-(2 k-1)^2}{(2 m-1)^2} \\ \quad\quad m=3,4, \cdots, M \end{array}\right.$

(A-2)

利用泰勒公式对M-SFD(N=4)给出的差分离散声波方程的时空域频散关系中的三角函数进行级数展开，然后使得方程左右两边$ {k}_{x}^{2}{k}_{z}^{2}{h}^{2} $、$ {k}_{x}^{2}{k}_{z}^{4}{h}^{4} $(或$ {k}_{x}^{4}{k}_{z}^{2}{h}^{4} $)、$ {k}_{x}^{4}{k}_{z}^{4}{h}^{6} $、$ {k}_{x}^{2}{k}_{z}^{6}{h}^{6} $(或$ {k}_{x}^{6}{k}_{z}^{2}{h}^{6} $)和$ {k}_{x}^{2j+2}{h}^{2j} $(或$ {k}_{z}^{2j+2}{h}^{2j} $)(j=0，1，2，…，M-1)的系数对应相等以构建关于差分系数的方程组，求解此方程组可以得到差分系数$ {a}_{1}, {a}_{2}, \dots , {a}_{M};{b}_{1}, {b}_{2}, {b}_{3}, {b}_{4} $的通解为

$ \left\{\begin{array}{l} b_1=\frac{1}{240 \times 16}\left(-\frac{r^6}{7}-15 r^4+\frac{629}{3} r^2\right) \\ b_2=\frac{1}{240 \times 48}\left(-\frac{31 r^6}{7}+87 r^4-239 r^2\right) \\ b_3=\frac{1}{240 \times 64}\left(\frac{25 r^6}{7}-57 r^4+\frac{457 r^2}{3}\right) \\ b_4=\frac{1}{240 \times 24 \times 8}\left(r^6-15 r^4+37 r^2\right) \\ a_1=\prod\limits_{2 \leqslant k \leqslant M} \frac{r^2-(2 k-1)^2}{1-(2 k-1)^2}-2 b_1-2 b_3 \\ a_2=\frac{1}{3} \prod\limits_{\substack{1 \leqslant k \leqslant M \\ k \neq 2}} \frac{r^2-(2 k-1)^2}{3^2-(2 k-1)^2}-2 b_2-2 b_4 \\ a_m=\frac{1}{2 m-1} \prod\limits_{\substack{1 \leqslant k \leqslant M \\ k \neq m}} \frac{r^2-(2 m-1)^2-(2 k-1)^2}{(2 m-1} \\ \quad\quad\quad m=3,4, \cdots, M \end{array}\right. $

(A-3)