0 引言

波动方程数值模拟是勘探地震学重要的基础研究内容。在地震数据采集阶段，可用于优化野外地震观测系统；在地震数据处理阶段，可用于检验处理方法的合理性，优化处理流程；在地震资料解释阶段，可用于检验和论证解释结果的正确性。同时，波动方程数值模拟还是逆时偏移^[1-6]和全波形反演^[7-10]的基础及核心内容。

波动方程数值模拟的方法主要有伪谱法^[11-14]、有限元法^[15-17]和有限差分法^[18-21]。伪谱法采用正反傅里叶变换实现空间偏导数的全局求解，虽然其模拟精度较高，但不适用于物性参数变化剧烈的模型，并且过多的傅里叶变换导致计算量大且计算效率低。有限元法基于区域分割和变分原理对等效积分弱形式的波动方程进行求解，进而模拟地震波的传播，该方法的模拟精度较高，并且对复杂模型具有良好的适用性，但由于其在迭代的每个时刻都需要求解大型矩阵的逆，因此占用内存较多、计算量巨大。有限差分法的本质是用差分代替微分从而近似求解波动方程，该方法具有适用广泛、占用内存少、计算效率高及编程实现简单等优点，因此，有限差分法是目前应用最普遍的一种波动方程数值模拟方法。

Alford等^[19]在1974年研究了声波方程有限差分数值模拟的精度问题，指出有限差分法存在数值频散(网格频散)现象，严重影响数值模拟的精度。此后，压制数值频散、提高模拟精度成为波动方程有限差分数值模拟的重要研究内容。Dablain^[22] 1986年指出数值频散可以分为时间数值频散(时间频散)和空间数值频散(空间频散)。时间频散使得相速度偏大，模拟波场中会出现“相位超前”的现象；空间频散使得相速度偏小，模拟波场中会出现“相位滞后”的现象。同时，Dablain^[22]还指出高阶有限差分方案能有效减小数值频散、提高模拟精度，但时间高阶有限差分方案会显著增加计算量，并使稳定性降低。因此，时间高阶有限差分方案很少被使用，发展空间高阶有限差分方案是提高模拟精度的有效途径。

Fornberg ^[23]和刘洋等^[24]1998年均给出了任意偶数阶有限差分法及其差分系数计算方法，被称为传统2 M 阶有限差分法(T2 M-FD)。T2 M-FD的基本思想是将波动方程中的Laplace算子表示为常规直角坐标系中坐标轴上的网格点差分近似的M 个Laplace算子的加权平均，是应用最普遍的一种高精度有限差分数值模拟方法。为了减小数值频散，提高数值模拟精度，交错网格有限差分法、隐式有限差分法和低秩有限差分法也相继被提出。交错网格有限差分法主要针对一阶应力速度形式的波动方程，相比规则网格有限差分法，它具有更高的模拟精度和更强的稳定性^[25-26]。隐式有限差分法需要求解大型线性方程组，因此占用内存多、计算效率低^[27-29]。低秩有限差分法需要作大量的空间傅里叶变换，计算效率也很低^[30-31]。

Liu等^[32]2009年指出波动方程是在时间空间域同时求解的，而T2 M-FD仅基于空间域频散关系计算差分系数，存在不合理性，本质上只具有二阶模拟精度。在保持差分格式不变的情况下，Liu等^[32]提出基于时间空间域频散关系和平面波理论计算差分系数，这种改进的差分方法被称为时空域2 M阶有限差分法(TS2 M-FD)。在二维模拟中，TS2 M-FD沿8个传播方向能够达到2 M 阶差分精度；在三维模拟中，TS2 M-FD沿48个方向能够达到2 M 阶差分精度；但无论是二维或是三维模拟中，除这些传播方向以外，TS2 M-FD沿其他方向传播时，仍然只具有2阶模拟精度。相比T2 M-FD，TS2 M-FD的模拟精度明显提高，并且稳定性也得到增强。Liu等^[33]在2013年进一步改进了时间空间域2 M 阶有限差分法，提出了菱形网格有限差分法。这种差分方法沿任意传播方向都能够达到2 M 阶模拟精度，使得模拟精度进一步提高；然而，增大菱形网格阶数的同时，计算量也会急剧增加，并且这种菱形网格有限差分法很难推广到三维波动方程的数值模拟。

在频率空间域，普遍采用混合网格有限差分格式求解波动方程^[34-40]。混合网格有限差分格式的基本出发点是充分利用旋转直角坐标系中距离差分中心点更近的网格点，对波动方程进行差分离散，将波动方程中的Laplace算子表示为常规直角坐标系和旋转直角坐标系中的网格点差分近似出的Laplace算子的加权平均。Jo等^[34]1996年提出的混合9点有限差分和Shin等^[35]1998年提出的混合25点有限差分均是较具代表性的2种混合网格有限差分格式。混合9点有限差分格式将Laplace算子表示为常规直角坐标系中差分近似出的1个Laplace算子和45°旋转直角坐标系中差分近似出的1个Laplace算子的加权平均。混合25点有限差分格式将Laplace算子表示为常规直角坐标系中差分近似出的2个Laplace算子、45°旋转直角坐标系中差分近似出的2个Laplace算子、26.6°旋转直角坐标系中差分近似出的1个Laplace算子以及63.4°旋转直角坐标系中差分近似出的1个Laplace算子，合计6个Laplace算子的加权平均。混合网格有限差分格式在兼顾计算效率的同时有效地提高了数值模拟的精度。

在T2 M-FD的基础上，将频率空间域混合网格有限差分格式的基本思想引入到时间空间域，即将波动方程中的Laplace算子表示为常规直角坐标系中坐标轴上的网格点差分近似出的M 个Laplace算子和不同角度旋转直角坐标系中的网格点差分近似出的N 个Laplace算子的加权平均，构建出时间空间域波动方程混合网格有限差分格式，并称之为混合网格2 M+N 型有限差分格式(M2 M+N-FD)。借鉴Liu等^[32-33]提出的基于时间空间域频散关系和平面波理论的差分系数计算方法，推导出M2 M+N-FD的差分系数计算方法，进行频散及稳定性分析，并与T2 M-FD及TS2 M-FD的相关结果进行对比。利用均匀介质模型和层状介质模型进行数值模拟试验，以期证实M2 M+N-FD能够更有效地压制数值频散、提高模拟精度。最后，将M2 M+N-FD应用于Marmousi模型的逆时偏移成像，以验证该差分方法的有效性和普遍适用性。

1 混合网格有限差分格式构建

差分格式构建是波动方程混合网格有限差分数值模拟的基础。在阐述波动方程传统2 M 阶有限差分法基本原理的基础上，引出构建混合网格2 M+N 型有限差分格式的基本思想及方法原理。

1.1 传统2 M 阶有限差分格式

地震波在二维常密度声学介质中传播时，波动方程可以表示为

$ \frac{{{\partial ^2}P}}{{\partial {x^2}}} + \frac{{{\partial ^2}P}}{{\partial {z^2}}} = \frac{1}{{{v^2}}}\frac{{{\partial ^2}P}}{{\partial {t^2}}} $

(1)

式中：P = P(x, z, t)为标量声波波场；v = v(x, z)为地震波在介质中的传播速度；${\nabla ^2}P = \frac{{{\partial ^2}P}}{{\partial \;{x^2}}} + \frac{{{\partial ^2}P}}{{\partial \;{z^2}}} $为Laplace算子；t为时间。

T2 M-FD采用的是时间二阶、空间2 M 阶的有限差分方案，对式(1)中的二阶时间偏导数进行二阶有限差分离散近似，可以得到

$ \frac{{{\partial ^2}P}}{{{\partial ^2}t}} \approx \frac{1}{{{\tau ^2}}}\left( {P_{0,0}^1 - 2P_{0,0}^0 + P_{0,0}^{ - 1}} \right) $

(2)

式中：P_{m, n}^j = P(x + mh, z + nh, t + jτ)，τ为时间采样间隔，h 为空间采样间隔；P_{0, 0}⁰ = P(x, z, t)为任意参考时刻t和参考位置(x, z)处的波场值。

构建差分格式过程中，P_{0, 0}⁰表示差分中心点(图 1)。

T2 M-FD采用图 1所示的差分方案对式(1)中的Laplace算子进行差分离散近似，可以得到M 种Laplace算子的差分近似表示方法

$ \begin{array}{l} {\nabla ^2}P = \frac{{{\partial ^2}P}}{{\partial {x^2}}} + \frac{{{\partial ^2}P}}{{\partial {z^2}}}\\ \;\;\;\;\;\;\; \approx \frac{1}{{{m^2}{h^2}}}\left( {P_{m,0}^0 + P_{0,m}^0 - 4P_{0,0}^0 + P_{0, - m}^0 + P_{ - m,0}^0} \right)\\ \;\;\;\;\;\;\;\;\;\;\left( {m = 1,2, \cdots ,M} \right) \end{array} $

(3)

式中：h 为空间采样间隔。

将这M 个Laplace差分近似表达式进行加权平均，可以得到

$ \begin{array}{l} {\nabla ^2}P = \frac{{{\partial ^2}P}}{{\partial {x^2}}} + \frac{{{\partial ^2}P}}{{\partial {z^2}}}\\ \;\;\;\;\;\;\; \approx \sum\limits_{m = 1}^M {\frac{{{c_m}}}{{{m^2}{h^2}}}\left( {P_{m,0}^0 + P_{0,m}^0 - 4P_{0,0}^0 + P_{0, - m}^0 + P_{ - m,0}^0} \right)} \end{array} $

(4)

式中：c_m 为权系数。

将式(2)和式(4)代入式(1)可以得到

$ \begin{array}{*{20}{c}} {\frac{1}{{{h^2}}}\sum\limits_{m = 1}^M {{a_m}\left( {P_{m,0}^0 + P_{0,m}^0 - 4P_{0,0}^0 + P_{0, - m}^0 + P_{ - m,0}^0} \right)} \approx }\\ {\frac{1}{{{v^2}{\tau ^2}}}\left( {P_{0,0}^1 - P_{0,0}^0 + P_{0,0}^{ - 1}} \right)} \end{array} $

(5)

式中：a_m = c_m/m²，为差分系数。

式(5)是T2 M-FD对式(1)的有限差分离散形式。TS2 M-FD和T2 M-FD的差分格式完全相同，因此它们相应的波动方程有限差分离散形式也完全相同，但它们的差分系数计算方法不同，T2 MFD是基于空间域频散关系计算的差分系数，TS2 M-FD是基于时间空间域频散关系和平面波理论计算的差分系数^[28]。

1.2 混合网格2 M+N 型有限差分格式

T2 M-FD和TS2 M-FD均采用图 1所示的差分方案对波动方程进行差分离散，它们仅利用常规直角坐标系中坐标轴上的网格点来差分近似波动方程中的Laplace算子，主要通过增大M 的取值来提高数值模拟的精度，随着M 取值的增大，新增加的网格点距离差分中心点越来越远，从而对提高模拟精度的贡献越来越小。

受频率空间域混合网格有限差分格式的启发，旋转直角坐标系中的网格点也可以用来差分近似Laplace算子。利用图 2(a)所示45°旋转直角坐标系中的网格点差分近似Laplace算子，可以得到

$ {\nabla ^2}P = \frac{{{\partial ^2}P}}{{\partial {x^2}}} + \frac{{{\partial ^2}P}}{{\partial {z^2}}} \approx \frac{1}{{2{h^2}}}\left( {P_{1, - 1}^0 + P_{1,1}^0 - 4P_{0,0}^0 + P_{ - 1, - 1}^0 + P_{ - 1,1}^0} \right) $

(6)

具体推导过程见附录A。

利用图 2(b)所示26.6°和63.4°旋转直角坐标系中的网格点差分近似Laplace算子，可以得到

$ \begin{array}{l} {\nabla ^2}P \approx \frac{1}{{10{h^2}}}\left[ {\left( {P_{2, - 1}^0 + P_{1,2}^0 - 4P_{0,0}^0 + P_{ - 1, - 2}^0 + P_{ - 2,1}^0} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left( {P_{1, - 2}^0 + P_{2,1}^0 - 4P_{0,0}^0 + P_{ - 2, - 1}^0 + P_{ - 1,2}^0} \right)} \right] \end{array} $

(7)

具体推导过程见附录B。

同样，还可以推导出图 2(c)~(d)所示旋转直角坐标系中的网格点差分近似Laplace算子的表达式。按照相同方法，还能构造出更为丰富的旋转网格有限差分格式。

在图 1所示T2 M-FD的基础上，依次引入图 2各个旋转直角坐标系中的网格点，可以构造出如图 3所示的混合网格有限差分方案。即：将图 1和图 2(a)中的差分格式进行组合，可以得到图 3(a)所示的M2 M+N-FD (N= 1)；将图 1和图 2(a)~(b)中的差分格式进行组合，可得到图 3(b)所示的M2 M+ N-FD(N = 2)；将图 1和图 2(a)~(c)中的差分格式进行组合，可以得到图 3(c)所示的M2 M+N-FD (N = 3)；将图 1和图 2(a)~(d)中的差分格式进行组合，可得到图 3(d)所示的M2 M+N-FD(N = 4)。若对图 1继续组合更多的旋转网格有限差分格式、增大N 的取值，便可以构建出更为丰富的混合网格有限差分格式。

混合网格有限差分格式可以充分利用位于旋转直角坐标系、但距离差分中心点更近的网格点来差分近似波动方程中的Laplace算子，这在理论上更为合理。以M2 M+N-FD (N= 1)为例，进一步阐述混合网格有限差分的基本方法和原理。

图 3(a)所示的M2 M+N-FD(N= 1)由图 1和图 2(a)中的差分格式组合而成，将式(4)和式(6)加权求和，可以得到M2 M+N-FD(N = 1)差分近似出的Laplace算子表达式

$ \begin{array}{*{20}{c}} {{\nabla ^2}P \approx \sum\limits_{m = 1}^M {\frac{{{c_m}}}{{{m^2}{h^2}}}\left( {P_{m,0}^0 + P_{0,m}^0 - 4P_{0,0}^0 + P_{0, - m}^0 + P_{ - m,0}^0} \right)} + }\\ {\frac{{{c_{1,1}}}}{{2{h^2}}}\left( {P_{1, - 1}^0 + P_{1,1}^0 - 4P_{0,0}^0 + P_{ - 1, - 1}^0 + P_{ - 1,1}^0} \right)} \end{array} $

(8)

式中：c_m (m= 1，2，…，M)，c_1，1为权系数。

式(8)表明，M2 M+N-FD(N = 1)是将Laplace算子表示为常规直角坐标系中差分近似的M 个Laplace算子和45°旋转直角坐标系中差分近似的1个Laplace算子的加权平均。同理，M2 M+N-FD本质上是将Laplace算子表示为常规直角坐标系中差分近似的M 个Laplace算子和不同角度旋转直角坐标系中差分近似的N 个Laplace算子的加权平均，这是构建混合网格有限差分格式的基本方法。

将式(2)和式(8)代入式(1)可以得到

$ \begin{array}{*{20}{c}} {\frac{1}{{{h^2}}}\sum\limits_{m = 1}^M {{a_m}\left( {P_{m,0}^0 + P_{0,m}^0 - 4P_{0,0}^0 + P_{0, - m}^0 + P_{ - m,0}^0} \right)} + }\\ {\frac{{{a_{1,1}}}}{{{h^2}}}\left( {P_{1, - 1}^0 + P_{1,1}^0 - 4P_{0,0}^0 + P_{ - 1, - 1}^0 + P_{ - 1,1}^0} \right) \approx }\\ {\frac{1}{{{v^2}{\tau ^2}}}\left( {P_{0,0}^1 - P_{0,0}^0 + P_{0,0}^{ - 1}} \right)} \end{array} $

(9)

式中：a_m = c_m/m²；a_1，1 = c_1，1/2；a_m(m = 1，2，…，M)，a_1，1为差分系数。

式(9)是M2 M+N-FD(N = 1)对式(1)的有限差分离散形式。利用同样的方法，可以推导出当N 取不同值时，其他M2 M+N-FD相应的波动方程的有限差分离散形式。

2 差分系数计算

差分系数计算是有限差分数值模拟重要的研究内容，计算方法的优劣将会直接影响数值模拟的精度和稳定性。本节继续以M2 M+N-FD (N= 1)为例，来阐述M2 M+N-FD差分系数计算的原理和方法。

式(1)在均匀介质中存在平面波解，其离散形式为

$ P_{m,n}^j = {{\rm{e}}^{i\left[ {{k_x}\left( {x + mh} \right) + {k_z}\left( {z + nh} \right) - \omega \left( {t + j\tau } \right)} \right]}} $

(10)

$ {k_x} = k\cos \theta ,{k_z} = k\sin \theta $

(11)

式中：k 为波数；k_x ，k_z 分别为水平波数和垂直波数；ω 为圆频率；θ为平面波传播方向与x 轴正半轴的夹角；i 为虚单位，即i² = -1。

将式(10)代入(9)式，利用Euler公式e^iθ = cos θ + i sin θ，可以得到

$ \begin{array}{l} \frac{1}{{{h^2}}}\sum\limits_{m = 1}^M {{a_m}\left[ {\cos \left( {m{k_x}h} \right) + \cos \left( {m{k_z}h} \right) - 2} \right]} + \\ \;\;\;\;\;\;\;\;\frac{{{a_{1,1}}}}{{{h^2}}}\left[ {\cos \left( {{k_x}h - {k_z}h} \right) + \cos \left( {{k_x}h + {k_z}h} \right)} \right] \approx \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{{{v^2}{\tau ^2}}}\left[ {\cos \left( {\omega \tau } \right) - 1} \right] \end{array} $

(12)

通过泰勒公式，展开式(12)中的余弦函数，可得到

$ \begin{array}{*{20}{c}} {\sum\limits_{n = 1}^\infty {\sum\limits_{m = 1}^M {{m^{2n}}{a_m}\frac{{{{\left( { - 1} \right)}^n}k_x^{2n}{h^{2n - 2}}}}{{\left( {2n} \right)\mathit{!}}}} } + \sum\limits_{n = 1}^\infty {\sum\limits_{m = 1}^M {{m^{2n}}{a_m}\frac{{{{\left( { - 1} \right)}^n}k_z^{2n}{h^{2n - 2}}}}{{\left( {2n} \right)\mathit{!}}}} } + }\\ {{a_{1,1}}\left\{ {\sum\limits_{m = 1}^M {\frac{{{{\left( { - 1} \right)}^n}{{\left( {{k_x} - {k_z}} \right)}^{2n}}{h^{2n - 2}}}}{{\left( {2n} \right)\mathit{!}}}} + \sum\limits_{m = 1}^M {\frac{{{{\left( { - 1} \right)}^n}{{\left( {{k_x} + {k_z}} \right)}^{2n}}{h^{2n - 2}}}}{{\left( {2n} \right)\mathit{!}}}} } \right\} \approx }\\ {\sum\limits_{n = 1}^\infty {{{\left( { - 1} \right)}^n}\frac{{{r^{2n - 2}}{k^{2n}}{h^{2n - 2}}}}{{\left( {2n} \right)\mathit{!}}}} } \end{array} $

(13)

式中：r = vτ/h，为Courant稳定性条件。

使式(13)中左右两边k_x²k_z²的系数对应相等，保留二项式展开过程中二项式系数组合数的表示形式，可以得到

$ C_4^2\left[ {{1^2} \times {{\left( { - 1} \right)}^2} + {1^2} \times {1^2}} \right]{a_{1,1}} = C_2^1{r^2} $

(14)

式中：C₄²和C₂¹为组合数。

采用组合数的表示形式，有助于更便利地将这种计算差分系数的方法推广至其他M2 M+N-FD (N ≠ 1)。

使式(13)中左右两边k_x²ⁿ或k_z²ⁿ的系数对应相等，可以得到

$ \sum\limits_{m = 1}^M {{m^{2n}}{a_m}} + {a_{1,1}}\left[ {{1^2} + {1^2}} \right] = {r^{2n - 2}}\;\;\;\;\left( {n = 1,2, \cdots ,M} \right) $

(15)

将式(15)写成矩阵形式，可以得到

$ \left[ {\begin{array}{*{20}{c}} \begin{array}{l} {1^2}\\ {1^4}\\ \vdots \\ {1^{2M}} \end{array}&\begin{array}{l} {2^2}\\ {2^4}\\ \vdots \\ {2^{2M}} \end{array}&\begin{array}{l} \cdots \\ \cdots \\ \\ \cdots \end{array}&\begin{array}{l} {M^2}\\ {M^4}\\ \vdots \\ {M^{2M}} \end{array} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{a_1}}\\ {{a_2}}\\ \vdots \\ {{a_3}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {1 - 2{a_{1,1}}}\\ {{r^2} - 2{a_{1,1}}}\\ \vdots \\ {{r^{2M - 2}} - 2{a_{1,1}}} \end{array}} \right] $

(16)

求解式(14)可以得到差分系数a_1，1= r²/6，再将a_1，1的解代入式(16)，便可以求解出差分系数a₁，a₂，…，a_M；这样M2 M+N-FD (N= 1)的M+1个差分系数a₁，a₂，…，a_M 及a_1，1就全部求解出来了。利用相同的方法，可以求解当N 取不同值时其他M2 M+N-FD的差分系数。

M2 M+N-FD的差分系数与M，N 及r 的取值有关，TS2 M-FD的差分系数与M 及r 的取值有关，T2 M-FD的差分系数仅与M 的取值有关。

表 1给出了T2 M-FD (M= 6)，TS2 M-FD (M= 6)和M2 M+N-FD(M = 5；N = 1，2，3，4)共6种差分格式的差分系数，计算差分系数时采用的速度及时间、空间采样间隔参数分别为：v = 3 000 m/s，τ= 0.001 s，h= 10 m。

3 频散分析

地震波在非耗散介质中传播时不会产生频散现象，即不同频率成分的地震波传播速度相同，但利用有限差分法求解波动方程时，对波动方程进行有限差分离散，会导致不同频率成分的地震波传播速度不再相同，即出现数值频散现象。时间数值频散(时间频散)使相速度变大而出现“相位超前”，空间数值频散(空间频散)使相速度变小而出现“相位滞后”。数值频散是利用有限差分求解波动方程时的固有特性，只能减小，无法消除，且数值频散的大小直接反映波动方程有限差分数值模拟的精度。

通常使用归一化相速度δ来描述有限差分格式的数值频散

$ \delta = \frac{{{v_{{\rm{ph}}}}}}{v} $

(17)

式中：v_ph为相速度。

结合相速度的定义(v_ph = ω/k)，通过式(12)和式(17)可给出M2 M+N-FD (N= 1)的数值频散关系式

$ \delta = \frac{{{v_{{\rm{ph}}}}}}{v} = \frac{\omega }{{vk}} = \frac{G}{{2{\rm{ \mathsf{ π} }}r}}\arccos \left( {{r^2}c + 1} \right) $

(18)

$ \begin{array}{l} c = \sum\limits_{m = 1}^M {{a_m}\left[ {\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}m}}{G}\cos \theta } \right) + \cos \left( {\frac{{2{\rm{ \mathsf{ π} }}m}}{G}\sin \theta } \right) - 2} \right]} + \\ \;\;\;\;\;{a_{1,1}}\left[ {\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}}}{G}\cos \theta - \frac{{2{\rm{ \mathsf{ π} }}}}{G}\sin \theta } \right) + } \right.\\ \;\;\;\;\;\left. {\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}}}{G}\cos \theta + \frac{{2{\rm{ \mathsf{ π} }}}}{G}\sin \theta } \right) - 2} \right] \end{array} $

(19)

式中：G = λ/h；λ为波长；G 为单位波长内网格点的个数。

δ值与1越接近，数值频散误差越小。δ > 1时，表明相速度偏大，存在时间数值频散，模拟波场中会出现“相位超前”的现象；δ < 1时，表明相速度偏小，存在空间数值频散，模拟波场中会出现“相位滞后”的现象。利用相同方法，可以推导出T2 M-FD，TS2 M-FD和其他M2 M+N-FD的频散关系。

图 4为T2 M-FD(M = 3，6，10，15)，TS2 M-FD (M = 3，6，10，15)和M2 M+N-FD(M = 3，6，10，15；N = 1)的频散曲线。计算频散曲线时采用的速度和时间、空间采样间隔参数分别为：v= 3 000 m/s，τ= 0.001 s，h = 10 m。对比分析图 4中的频散曲线，可以看出：

(1) T2 M-FD (M= 3，6，10，15)的相速度与真实速度之比约等于1，对应1/G 的区间很小，大致为0~0.075。增大M 取值，这一区间也并没有随之增大。这说明，增大M 取值不能增大T2 M-FD的差分精度阶数。

(2) TS2 M-FD(M = 3，6，10，15)的相速度与真实速度之比约等于1，对应1/G 的区间较小，大致为0~0.125。增大M 取值，这一区间同样也没有随之增大。因此，增大M 取值，也不能增大TS2 M-FD的差分精度阶数。与T2 M-FD(M = 3，6，10，15)相比，TS2 M-FD(M = 3，6，10，15)的相速度与真实速度之比约等于1时，对应1/G 的区间明显增大，说明TS2 M-FD能更有效地压制数值频散，提高模拟精度。

(3) T2 M-FD和TS2 M-FD在M 取值较小时(M= 3左右)，增大M 取值，虽然不能增大相速度与真实速度之比约等于1时对应1/G 的区间，即不能增大差分精度的阶数，但能够减小数值频散的幅值；当M 取值较大时(M ≥ 6)，继续增大M 取值，既不能增大相速度与真实速度之比约等于1时对应1 G 的区间，也基本不能进一步减小数值频散的幅值。这说明，T2 M-FD和TS2 M-FD不能仅依靠增大M 取值来持续提高数值模拟精度，而应该寻找新的差分格式，这正是本次研究提出M2 M+N-FD的基本出发点。

(4) M2 M+N-FD (M= 3；N= 1)的相速度与真实速度之比约等于1，对应1/G 的区间大致为0~0.175；当M ≥ 6时，这一区间增大为0~0.25，是TS2 M-FD的2倍、T2 M-FD的3.33倍。因此，M2 M+N-FD能够更有效的压制数值频散，提高模拟精度。

图 5为M2 M+N-FD (M= 3，10，15；N= 1，2，3，4)的频散曲线，计算频散曲线时采用的速度和时间、空间采样间隔参数分别为：v = 3 000 m/s，τ= 0.001 s，h= 10 m。需要注意，图 5(a)~(d)、图 5(e)~(h)和图 5(i)~(l)这3组频散曲线子图的纵轴刻度不相同。对比分析这3组频散曲线，可以看出：

(1) 对比图 5(a)~(d)中M2 M+N-FD(M = 6；N = 1，2，3，4)的各频散曲线，特征基本相同，数值模拟精度也基本相当。由于增大N 取值的同时会增大数值模拟的计算量、降低计算效率，因此，当M 的取值不是足够大，如M = 6时，建议采用M2 M+N-FD(N = 1)进行数值模拟。这样既可以保证模拟精度，又可以兼顾计算效率。

(2) 对比图 5(e)~(h)中M2 M+N-FD(M = 10；N= 1，2，3，4)的各频散曲线，当N 取值由1增大到2时，数值频散的误差明显减小，继续将N 取值增大至3或4时，频散误差却不再减小，因此，当M 的取值较大，如M= 10时，可以采用M2 M+N-FD (N= 2)进行数值模拟，而没有必要采用M2 M+N-FD(N = 3，4)。

(3) 对比图 5(i)~(l)中M2 M+N-FD(M = 15；N= 1，2，3，4)的各频散曲线，随着N 取值的增大，数值频散的误差明显减小，因此，当M 的取值很大，如M = 15时，可以采用M2 M+N-FD(M = 4)进行数值模拟。

(4) 将图 5(b)和图 5(c)、图 5(f)和图 5(g)、图 5(j)和图 5(k)两两进行对比，发现频散曲线特征的差别微小，这是由于按照传统差分精度阶数分析方法，M2 M+N-FD(N = 2)和M2 M+N-FD(N = 3)的差分精度阶数相同，均为6阶差分精度。这里只给出结论，不再给出差分精度阶数的相关公式推导。

建议多数情况下采用M2 M+N-FD(N = 1)，且M 取值在6左右，进行数值模拟；当对模拟精度有非常苛刻的要求并且传播距离很长时，可以考虑采用M2 M+N-FD(N = 2，4)，且M 取大于10以上的值，进行数值模拟；N 取值大于4的M2 M+N-FD基本不会被采用。

前面提到Liu等^[33]2013年提出的菱形网格有限差分格式是M2 M+N-FD的特殊形式(M 和N 取特定值时)，四阶菱形网格有限差分格式等价于M2 M+N-FD(M = 4；N = 4)。在菱形网格有限差分格式中，增大M 的取值时，N 的取值也随之急剧增大。前文分析表明：过大的N 值通常不能有效地进一步减小数值频散、提高模拟精度，反而会增大计算量、降低计算效率；M2 M+N-FD更具一般性，能够合理地自由选择M 和N 的取值，从而有效兼顾模拟精度和计算效率。

图 6为T2 M-FD(M = 6，7)，TS2 M-FD(M = 6，7)和M2 M+N-FD (M= 5，6；N= 1)的频散曲线，计算频散曲线时采用的速度和时间与空间采样间隔参数分别为：v= 3 000 m/s，τ = 0.001 s，h= 10 m。

由于T2 M-FD(M = 6)，TS2 M-FD(M = 6)和M2 M+N-FD(M = 5；N = 1)这3种有限差分格式的Laplace算子长度均为25点，因此三者具有基本相同的计算量。对比分析这3种有限差分格式的频散曲线[图 6(a)~(c)]，可以看出：

(1) T2 M-FD(M = 6)具有严重的数值频散现象，主要为时间频散(δ > 1)。

(2) 与T2 M-FD (M = 6)相比，TS2 M-FD (M = 6)的频散误差幅值明显减小，但频散曲线仍然发散，具有一定的时间频散和空间频散。TS2 M-FD (M=6)的频散特性与地震波传播的方向角θ有关：θ = 0°时，主要为空间频散(δ < 1)；θ = 45°时，主要为时间频散(δ > 1)；θ从0°增大至45°的过程中，会依次出现空间频散减小、消失及时间频散增强等现象。

(3) M2 M+N-FD(M = 5；N = 1)的频散曲线收敛性较好，数值频散误差在这3种有限差分格式中最小。

综上所述，在计算量相同的情况下，T2 M-FD的模拟精度最低，模拟结果(波场快照或炮集)中会出现严重的时间频散(δ > 1，相位超前)现象；TS2 MFD的模拟精度中等，模拟结果中会同时出现一定的空间频散(δ < 1，相位滞后)和时间频散(δ > 1，相位超前)现象；M2 M+N-FD的模拟精度最高，数值频散最小。

T2 M-FD(M = 7)，TS2 M-FD(M = 7)和M2 M+ N-FD(M = 6；N = 1)有限差分格式的Laplace算子长度均为29点，数值模拟的计算量也基本相同。对比分析这3种频散曲线[图 6(d)~(f)]，结论和图 6(a)~(c)的对比结果一致。

4 稳定性分析

时间空间域有限差分数值模拟是通过利用差分代替微分迭代来求解偏微分波动方程的，只有迭代过程稳定收敛，这种代替才算合理。稳定性条件描述地震波在介质中传播的速度、空间及时间采样间隔与差分阶数及差分系数之间必须满足的关系，也称为Courant稳定性条件或CFL条件(CourantFriedrichs-Lewy Condition)。

根据式(12)可以得到

$ \begin{array}{l} \frac{1}{{{r^2}}}\left[ {1 - \cos \left( {\omega \tau } \right)} \right] \approx \sum\limits_{m = 1}^M {{a_m}\left[ {2 - \cos \left( {m{k_x}h} \right) - \cos \left( {m{k_z}h} \right)} \right]} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{a_{1,1}}\left[ {2 - \cos \left( {{k_x}h - {k_z}h} \right) - \cos \left( {{k_x}h + {k_z}h} \right)} \right] \end{array} $

(20)

取最大空间波数(尼奎斯特波数)k_x = k_z = π/h ，并考虑0 ≤ 1-cos (ω τ) ≤ 2，得到M2 M+N-FD (N= 1)的稳定性条件

$ r \le S = \frac{1}{{\sqrt {2\sum\limits_{m = 1}^{{M_1}} {{a_{2m - 1}}} } }} $

(21)

式中：S 为稳定性因子；$ {M_1} = {\rm{int}}\left[{\frac{{M + 1}}{2}} \right]$，int表示取整。

式(21)表明M2 M+N-FD(N = 1)的稳定性与M的取值及差分系数有关，同时，差分系数又与r 的取值相关。

利用相同方法，可以推导出T2 M-FD，TS2 MFD和其他M2 M+N-FD的稳定性条件及稳定性因子的表达式。

图 7为T2 M-FD，TS2 M-FD和M2 M+N-FD (N = 1，2，3，4)的稳定性因子随r 增大而变化的曲线。根据式(21)，稳定性条件要求r ≤ S，因此，稳定性因子S 随r 增大而变化的曲线与直线r 的交点即为该差分格式所允许的最大r 取值。所允许的r 取值越大，表明该差分格式的稳定性越强。

图 8为T2 M-FD，TS2 M-FD和M2 M+N-FD (N= 1，2，3，4)共6种差分格式满足稳定性条件r ≤ S的最大r 取值随M 取值的变化曲线。可看以出：T2 M-FD的稳定最弱，TS2 M-FD的稳定性较T2 M-FD有明显增强，M2 M+N-FD (N= 1)的稳定性与TS2 M-FD基本相当，M2 M+N-FD(N = 2，3，4)的稳定性比TS2 M-FD略强。因此，增大M2 M+N-FD中N的取值，会适当增强其稳定性，但同时会增大数值模拟的计算量。

5 数值模拟

频散分析结果表明：M2 M+N-FD能更有效地压制数值频散、提高模拟精度。为了验证这一结论，本节将采用T2 M-FD，TS2 M-FD和M2 M+N-FD分别在均匀介质模型和层状介质模型上进行数值模拟试验，并对模拟结果进行对比分析，最后将M2 M+N-FD推广应用于Marmousi模型的逆时偏移。

5.1 均匀介质模型数值模拟

设计1个简单的均匀介质模型，大小为nx×nz= 4 001×4 001，速度v = 3 000 m/s，空间采样间隔Δx= Δz = h = 10 m，震源位于模型中心(2 001，2 001)，震源子波为30 Hz雷克子波，时间采样间隔τ = 0.001 s。分别利用T2 M-FD (M= 7)，TS2 M-FD (M= 7)和M2 M+N-FD(M = 6；N = 1)对该均匀介质模型进行数值模拟，这3种有限差分格式的Laplace算子长度均为29点，具有基本相同的计算量。

考虑到模型的对称性，绘制波场快照时，只画出了右上角的1/4部分。图 9(a)~(c)分别为T2 M-FD(M = 7)，TS2 M-FD(M = 7)和M2 M+N-FD(M = 6；N = 1)在6.6 s时刻的波场快照，图 9(d)~(f)分别为图 9(a)~(c)中蓝色方框部分的放大显示。图 9(d)~(f)中，每一个波场快照自下而上对应的地震波传播方向角θ的变化范围均为0°~45°。可以看出：

(1) T2 M-FD (M= 7)正演模拟的波场快照中出现严重的“相位超前”现象，即存在严重的时间频散，在3种差分格式中模拟精度最低。这与图 6(a)中频散曲线分析给出的结论一致。

(2) 随着传播方向角θ由0°变化至45°，TS2 MFD(M = 7)正演模拟的波场快照中出现由“相位滞后”逐渐过渡为“相位超前”的现象，即由空间频散过渡到时间频散；相比T2 M-FD(M = 7)，TS2 M-FD (M= 7)的数值频散程度明显减弱，模拟精度明显提高。这与图 6(b)中频散曲线分析给出的结论完全一致。

(3) M2 M+N-FD(M = 6；N = 1)正演模拟的波场快照中无明显的数值频散，在3种差分格式中模拟精度最高。这与图 6(c)中频散曲线分析得到的结论也一致。

均匀介质模型数值模拟结果表明：计算量基本相等时，与T2 M-FD和TS2 M-FD相比，M2 M+N-FD能更有效地压制数值频散、提高模拟精度，进一步验证了频散分析给出的结论(参见图 6)。

5.2 层状介质模型数值模拟

设计1个简单的双层模型，模型参数为h₁= 8 km，v₁ = 2 600 m/s，h₂= 7 km，v₂ = 3 200 m/s，模型宽度为20 km，模型网格化参数为Δx =Δz = h = 10 m。震源位于模型中部，深度为10 m，震源子波为30 Hz雷克子波，时间采样间隔τ= 0.001 s。分别利用T2 M-FD (M = 7)，TS2 M-FD(M = 7)和M2 M+N-FD(M = 6；N = 1)对该双层模型进行数值模拟，这3种有限差分格式的Laplace算子长度均为29点，具有基本相同的计算量。

图 10(a)~(c)分别为T2 M-FD (M= 7)，TS2 MFD(M = 7)和M2 M+N-FD(M = 6；N = 1)在4.6 s时刻的波场快照，图 10(d)~(f)分别为图 10(a)~(c)中黑色方框部分的放大显示。可以看出：T2 M-FD (M= 7)数值模拟的波场快照中的直达波、反射波和透射波均存在严重的数值频散；TS2 M-FD(M = 7)数值模拟的波场快照中的直达波、反射波和透射波也存在一定的数值频散；M2 M+N-FD(M = 6；N = 1)数值模拟的波场快照中无明显的数值频散。

图 11(a)~(c)分别为采用T2 M-FD(M = 7)，TS2 M-FD(M = 7)和M2 M+N-FD(M = 6；N = 1)这3种差分方法进行层状介质模型数值模拟得到的炮集记录，图 11(d)~(f)分别为图 11(a)~(c)中黑色方框部分的放大显示，图 11(g)~(i)分别为图 11(a)~(c)中蓝色方框部分的放大显示。可看出：T2 M-FD (M= 7)存在严重的时间频散[图 11(d)和图 11(g)中出现严重的“相位超前”现象]，波形发生严重畸变；TS2 M-FD(M = 7)存在一定的时间频散[图 11(e)中出现“相位超前”现象]和空间频散[图 11(h)中出现“相位超前”现象]，波形发生一定程度的畸变；M2 M+N-FD(M = 6；N = 1)无明显数值频散，波形未发生畸变。

层状介质模型数值模拟试验再次表明：在计算量基本相等的情况下，T2 M-FD具有严重的时间频散，模拟精度最低；TS2 M-FD存在一定的时间和空间频散，模拟精度中等；M2 M+N-FD无明显数值频散，模拟精度最高。

5.3 Marmousi模型数值模拟及逆时偏移

为进一步验证混合网格有限差分法对复杂模型的适用性，利用M2 M+N-FD(M = 6；N = 1)对Marmousi模型进行正演模拟和逆时偏移(图 12)。

图 12(a)为Marmousi速度模型。模型网格nz× nx = 751×2 301，空间采样间隔Δx =Δz = h = 4 m。震源为20 Hz的雷克子波，时间采样间隔τ= 0.000 3 s。共正演115炮，炮点深度均为12 m，炮间距为80 m (即间隔20个网格点)，每炮1 150道接收，检波点深度也均为12 m，道间距为8 m。

图 12(b)和图 12(c)分别为第28炮和第58炮的正演模拟炮集，图中均无明显的数值频散、模拟精度较高。图 12(d)为逆时偏移成像结果，图中可见，成像精度较高。

Marmousi模型正演模拟和逆时偏移结果进一步证明了混合网格有限差分法具有普遍适用性并易于推广应用于逆时偏移成像处理。

6 结论

(1) 利用常规直角坐标系和旋转直角坐标系中的网格点一起差分近似波动方程中的Laplace算子，构建出一种新型的混合网格有限差分格式，即M2 M+N-FD；并推导出该混合网格有限差分格式基于时间空间域频散关系和平面波理论的差分系数计算方法。

(2) T2 M-FD和TS2 M-FD均是仅利用常规直角坐标系中坐标轴上的网格点差分近似波动方程中Laplace算子，当M 取值较小时，若增大M 的取值，能够有效减小数值频散的幅值、提高模拟精度；当M 取值较大时，若继续增大M 的取值，新增加的网格点距离中心点的距离会越来越远，而对提高模拟精度的贡献越来越小。本次研究提出的M2 M+ N-FD，能充分利用常规直角坐标系和旋转直角坐标系中网格点差分近似波动方程中的Laplace算子，与T2 M-FD和TS2 M-FD相比，M2 M+N-FD能更有效地减小数值频散、进一步提高模拟精度。通过合理选择M 和N 的取值，可以使M2 M+N-FD在获得较高数值模拟精度的同时兼顾计算效率。在一般的数值模拟中，选择M2 M+N-FD (N= 1)即可。

(3) 在计算量基本相等时，T2 M-FD具有严重的时间频散，模拟精度最低；TS2 M-FD具有一定的时间频散和空间频散，模拟精度中等；M2 M+N-FD的数值频散最小，模拟精度最高。稳定性分析表明，T2 M-FD的稳定性最差；TS2 M-FD的稳定性较T2 M-FD具有明显增强；M2 M+ N-FD(N = 1)与TS2 M-FD的稳定性基本相当；增大N 的取值时，会进一步增强M2 M+N-FD的稳定性，同时也会增加数值模拟的计算量。

(4) 均匀介质模型和层状介质模型的数值模拟试验进一步验证了频散分析得出的结论，即当计算量基本相等时，M2 M+N-FD能更有效地压制数值频散且模拟精度更高。将M2 M+N-FD应用于Marmousi模型的数值模拟和逆时偏移，高精度的数值模拟和偏移成像结果充分证明了M2 M+N-FD的有效性和普遍适用性。

(5) M2 M+N-FD只适用于空间采样间隔Δx = Δz 的情况，而Δx ≠Δz 的情况有待进一步研究。M2 M+N-FD是基于声波方程推导出的，目前尚未对各项异性和弹性波方程的适用性进行深入研究。

附录A

图 2(a) 45°旋转直角坐标系中的网格点表示Laplace算子。将利用二元函数的泰勒展开推导出利用45°旋转直角坐标系[参见图 2(a)]中的网格点差分近似Laplace算子的表达式。

对45°旋转直角坐标系[图 2(a)]中的4个绿色网格点利用二元函数的泰勒展开，省略三阶及以上的偏导项，可以得到

$ \begin{array}{l} P_{1,1}^0 = P\left( {x + \Delta x,z + \Delta z,t} \right) \approx {\rm{P}}\left( {x,z,t} \right) + {P_x}\left( {x,z,t} \right)\Delta x + \\ \;\;\;\;\;\;\;\;{P_z}\left( {x,z,t} \right)\Delta z + \frac{1}{2}{P_{xx}}\left( {x,z,t} \right)\Delta {x^2} + \frac{1}{2}{P_{zz}}\left( {x,z,t} \right)\Delta {z^2} + \\ \;\;\;\;\;\;\;\;{P_{xz}}\left( {x,z,t} \right)\Delta x\Delta z \end{array} $

(A-1)

利用相同的方法还可以得到

$ \begin{array}{l} P_{ - 1,1}^0 = P\left( {x - \Delta x,z + \Delta z,t} \right) \approx P\left( {x,z,t} \right) - {P_x}\left( {x,z,t} \right)\Delta x + \\ \;\;\;\;\;\;\;\;{P_z}\left( {x,z,t} \right)\Delta z + \frac{1}{2}{P_{xx}}\left( {x,z,t} \right)\Delta {x^2} + \frac{1}{2}{P_{zz}}\left( {x,z,t} \right)\Delta {z^2} - \\ \;\;\;\;\;\;\;\;{P_{xz}}\left( {x,z,t} \right)\Delta x\Delta z \end{array} $

(A-2)

$ \begin{array}{l} P_{ - 1, - 1}^0 = P\left( {x - \Delta x,z - \Delta z,t} \right) \approx {\rm{P}}\left( {x,z,t} \right) - {P_x}\left( {x,z,t} \right)\Delta x - \\ \;\;\;\;\;\;\;\;{P_z}\left( {x,z,t} \right)\Delta z + \frac{1}{2}{P_{xx}}\left( {x,z,t} \right)\Delta {x^2} + \frac{1}{2}{P_{zz}}\left( {x,z,t} \right)\Delta {z^2} + \\ \;\;\;\;\;\;\;\;{P_{xz}}\left( {x,z,t} \right)\Delta x\Delta z \end{array} $

(A-3)

$ \begin{array}{l} P_{1, - 1}^0 = P\left( {x + \Delta x,z - \Delta z,t} \right) \approx P\left( {x,z,t} \right) + {P_x}\left( {x,z,t} \right)\Delta x - \\ \;\;\;\;\;\;\;\;{P_z}\left( {x,z,t} \right)\Delta z + \frac{1}{2}{P_{xx}}\left( {x,z,t} \right)\Delta {x^2} + \frac{1}{2}{P_{zz}}\left( {x,z,t} \right)\Delta {z^2} - \\ \;\;\;\;\;\;\;\;{P_{xz}}\left( {x,z,t} \right)\Delta x\Delta z \end{array} $

(A-4)

将式(A-1) ~ (A-4)相加，可以得到

$ \begin{array}{l} P_{1,1}^0 + P_{ - 1,1}^0 + P_{ - 1, - 1}^0 + P_{1, - 1}^0 \approx 4P\left( {x,z,t} \right) + 2{P_{xx}}\left( {x,z,t} \right)\Delta {x^2} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;2{P_{zz}}\left( {x,z,t} \right)\Delta {z^2} \end{array} $

(A-5)

当Δx = Δz = h时，可以得到

$ {\nabla ^2}P = {P_{xx}} + {P_{zz}} \approx \frac{1}{{2{h^2}}}\left( {P_{1,1}^0 + P_{ - 1,1}^0 - 4P_{0,0}^0 + P_{1, - 1}^0 + P_{ - 1, - 1}^0} \right) $

(A-6)

式(A-6)为利用45°旋转直角坐标系[图 2(a)]中的网格点差分近似Laplace算子的表达式。

注意式(A-6)得出的前提条件Δx = Δz ，这是本次研究提出的混合网格有限差分格式只适用于纵横空间采样间隔相同的情况的原因。

附录B

图 2(b)中26.6°和63.4°旋转直角坐标系中的网格点表示Laplace算子。将利用二元函数的泰勒展开推导出利用26.6°和63.4°旋转直角坐标系[参见图 2(b)]中的网格点差分近似Laplace算子的表达式。

对图 2(b)中26.6°旋转直角坐标系中的4个绿色网格点($P_{2, 1}^0, P_{-1, 2}^0, P_{-2, -1}^0, P_{1, - 2}^0 $)利用二元函数的泰勒展开，省略三阶及以上的偏导项，可以得到

$ \begin{array}{l} P_{2,1}^0 = P\left( {x + 2\Delta x,z + \Delta z,t} \right) \approx P\left( {x,z,t} \right) + 2{P_x}\left( {x,z,t} \right)\Delta x + \\ \;\;\;\;\;\;\;\;{P_z}\left( {x,z,t} \right)\Delta z + 2{P_{xx}}\left( {x,z,t} \right)\Delta {x^2} + \frac{1}{2}{P_{zz}}\left( {x,z,t} \right)\Delta {z^2} + \\ \;\;\;\;\;\;\;\;2{P_{xz}}\left( {x,z,t} \right)\Delta x\Delta z \end{array} $

(B-7)

$ \begin{array}{l} P_{ - 1,2}^0 = P\left( {x - \Delta x,z + 2\Delta z,t} \right) \approx {\rm{P}}\left( {x,z,t} \right) - {P_x}\left( {x,z,t} \right)\Delta x + \\ \;\;\;\;\;\;\;\;2{P_z}\left( {x,z,t} \right)\Delta z + \frac{1}{2}{P_{xx}}\left( {x,z,t} \right)\Delta {x^2} + \\ \;\;\;\;\;\;\;\;2{P_{zz}}\left( {x,z,t} \right)\Delta {z^2} - 2{P_{xz}}\left( {x,z,t} \right)\Delta x\Delta z \end{array} $

(B-8)

$ \begin{array}{l} P_{ - 2, - 1}^0 = P\left( {x - 2\Delta x,z - \Delta z,t} \right) \approx P\left( {x,z,t} \right) - \\ \;\;\;\;\;\;\;\;2{P_x}\left( {x,z,t} \right)\Delta x - {P_z}\left( {x,z,t} \right)\Delta z + 2{P_{xx}}\left( {x,z,t} \right)\Delta {x^2} + \\ \;\;\;\;\;\;\;\;\frac{1}{2}{P_{zz}}\left( {x,z,t} \right)\Delta {z^2} + 2{P_{xz}}\left( {x,z,t} \right)\Delta x\Delta z \end{array} $

(B-9)

$ \begin{array}{l} P_{1, - 2}^0 = P\left( {x\Delta x,z - 2\Delta z,t} \right) \approx P\left( {x,z,t} \right) + {P_x}\left( {x,z,t} \right)\Delta x - \\ \;\;\;\;\;\;\;\;2{P_z}\left( {x,z,t} \right)\Delta z + \frac{1}{2}{P_{xx}}\left( {x,z,t} \right)\Delta {x^2} + \\ \;\;\;\;\;\;\;\;2{P_{zz}}\left( {x,z,t} \right)\Delta {z^2} - 2{P_{xz}}\left( {x,z,t} \right)\Delta x\Delta z \end{array} $

(B-10)

将式（B-7）~（B-10）相加，可以得到

$ \begin{array}{l} P_{2,1}^0 + P_{ - 1,2}^0 + P_{ - 2, - 1}^0 + P_{1, - 2}^0 \approx 4P\left( {x,z,t} \right) + 5{P_{xx}}\left( {x,z,t} \right)\Delta {x^2} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;5{P_{zz}}\left( {x,z,t} \right)\Delta {z^2} \end{array} $

(B-11)

当Δx = Δz = h时，可以得到

$ {\nabla ^2}P = {P_{xx}} + {P_{zz}} \approx \frac{1}{{5{h^2}}}\left( {P_{2,1}^0 + P_{1, - 2}^0 - 4P_{0,0}^0 + P_{ - 1,2}^0 + P_{ - 2, - 1}^0} \right) $

(B-12)

式(B-12)是利用图 2(b)中26.6°旋转直角坐标系中的4个绿色网格点($P_{2, 1}^0, P_{-1, 2}^0, P_{-2, -1}^0, P_{1, - 2}^0 $)和中心点P_{0, 0}⁰差分近似Laplace算子的表达式。

利用相同的方法，还可以推导出利用图 2(b)中63.4°旋转直角坐标系中的4个绿色网格点($ P_{1, 2}^0, P_{-2, 1}^0, P_{-1, -2}^0, P_{2, -1}^0$)和中心点P_{0, 0}⁰差分近似Laplace算子的表达式

$ {\nabla ^2}P = {P_{xx}} + {P_{zz}} \approx \frac{1}{{5{h^2}}}\left( {P_{1,2}^0 + P_{2, - 1}^0 - 4P_{0,0}^0 + P_{ - 2,1}^0 + P_{ - 1, - 2}^0} \right) $

(B-13)

将式(B-12)和式(B-13)相加，可以得到

$ \begin{array}{*{20}{c}} {{\nabla ^2}P \approx \frac{1}{{10{h^2}}}\left[ {\left( {P_{2,1}^0 + P_{1, - 2}^0 - 4P_{0,0}^0 + P_{ - 1,2}^0 + P_{ - 2, - 1}^0} \right) + } \right.}\\ {\left. {\left( {P_{1,2}^0 + P_{2, - 1}^0 - 4P_{0,0}^0 + P_{ - 2,1}^0 + P_{ - 1, - 2}^0} \right)} \right]} \end{array} $

(B-14)

图 2(b)由26.6°和63.4°等2套旋转直角坐标系组成。因此，式(B-14)为利用图 2(b)中26.6°和63.4°旋转直角坐标系中的网格点差分近似Laplace算子的表达式。