0 引言

逆时偏移的思想起源于数学上的利用有限差分求解波动方程^[1]，后被广泛应用到地震勘探领域^[2-4]。逆时偏移不受倾角限制，能对强横向变速地下结构进行成像。传统地震资料只记录弹性波场的垂直分量，忽略了大部分甚至全部的S波信息^[5]，相比之下，多分量地震资料的弹性波逆时偏移可以提供更多的地下构造信息。如PS波成像能够补偿P波能量弱的区域的信息，且通过比较PP与PS成像剖面可以验证亮点反射^[6]。目前海上多分量地震资料越来越多^[7]，随着地震采集技术的发展和计算机性能的提高，弹性波逆时偏移的应用也会越来越广泛。

实现逆时偏移成像一般需要三个步骤：震源波场的获取；接收波场的重构；利用成像条件提取成像值^[8-9]。因此，地震波场的精确求解是影响成像剖面质量的直接因素^[10-11]。由于通常采用有限差分求解波动方程，正演波场会存在数值频散^[12-14]，对后期成像和反演皆不利。压制数值频散的途径主要有三类：增加差分算子长度^[15-17]、通量校正传输^[18]、优化差分系数^[19-23]。前两类方法都能在一定程度上压制数值频散，但会引入额外的计算量，且不能保证算法的稳定性^[24]。优化差分系数的方法在不增加计算量、不破坏算法的稳定性和收敛性的前提下，从数值频散的本质出发，通过合理布设各个网格的权值分配，使数值解更接近于真实解，是压制数值频散的有效方法。传统上，常采用泰勒近似求取差分系数^[25]，由于直接忽略高阶项，导致其模拟结果数值频散严重。Liu等^[26]利用平面波原理推导了一种隐式的有限差分交错网格策略，提高了波场模拟的精度，但隐式格式需要求解更多的方程，难以用于逆时偏移；Yang等^[27]在空间—波数域构建了新的频散关系，在最小二乘意义下，求得一阶交错网格差分系数的最优解，有效地压制了空间频散，但该方法比传统方法的时间频散更严重；雍鹏等^[28]提出一种基于声波方程的等效交错网格方法，能同时压制空间频散和时间频散。

目前有限差分算法主要基于两类网格：规则网格和交错网格。对于弹性波数值模拟，交错网格往往具有更高的模拟精度^[29-30]，但目前大多数交错网格都基于一阶速度—应力方程，该方程含多个中间变量，在成像和反演过程中未被使用，因此浪费了计算机内存。等效交错网格^[31]使用两次一阶交错网格离散二阶空间偏导数，其模拟精度等价于常规的一阶交错网格。因此，使用等效交错网格策略，就可以利用二阶位移方程的交错网格实现弹性波正演、成像和反演，与传统一阶速度—应力方程的交错网格相比，有效节约了计算内存。

另一方面，实现逆时偏移，需要保存全部时刻的震源波场，尤其是三维介质的情况下，需要消耗大量的存储，因此超额的存储量是制约逆时偏移用于巨量地震资料处理的瓶颈^[32-33]。近年来许多学者对逆时偏移的存储问题进行了研究，并取得了一定的成效：Symes^[34]利用少量的最优检验点控制波场的存储方式，检验点之间的波场通过重新正演获取，实现了计算量和存储量的折中；Clapp^[35]提出一种边界存储策略，该方法只需保存边界区域的震源波场，反向传播时，利用边界区域的波场作为边界条件重构之前时刻的波场，该方法虽然增加了计算量，但大幅度降低了逆时偏移的存储需求。

本文基于优化差分系数的思想，提出了一种平面波优化差分系数的等效交错网格算法，主要目的在于利用低阶的优化差分算子替代传统高阶差分算子，在保证精度情况下，减少计算量。研究表明，本文的6、10阶空间差分算子分别能达到传统交错网格法10、18阶的精度，有效降低了计算量。在内存需求方面，将有效边界存储的策略引入弹性波逆时偏移，利用边界附近区域的部分波场值控制逆时偏移的波场存储方式，同时采用二阶位移弹性波方程，避免了中间变量，有效降低了内存占用。模型测试证实了本文方法相比于传统交错网格法在计算量和内存占用方面具有优势。

1 基于平面波优化的等效交错网格有限差分法 1.1 等效交错网格差分格式

在均匀各向同性弹性介质中，用位移表示的二维弹性波动方程为

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}u}}{{\partial {t^2}}} = v_{\rm{P}}^2\frac{{{\partial ^2}u}}{{\partial {x^2}}} + v_{\rm{S}}^2\frac{{{\partial ^2}u}}{{\partial {z^2}}} + (v_{\rm{P}}^2 - v_{\rm{S}}^2)\frac{{{\partial ^2}w}}{{\partial x\partial z}}}\\ {\frac{{{\partial ^2}w}}{{\partial {t^2}}} = v_{\rm{P}}^2\frac{{{\partial ^2}w}}{{\partial {z^2}}} + v_{\rm{S}}^2\frac{{{\partial ^2}w}}{{\partial {x^2}}} + (v_{\rm{P}}^2 - v_{\rm{S}}^2)\frac{{{\partial ^2}u}}{{\partial x\partial z}}} \end{array}} \right. $

(1)

式中：u、w分别为在x、z方向的位移分量；t表示波传播的时刻；v_P、v_S分别为P波和S波速度。

以式(1)为例，说明等效交错网格有限差分法近似空间偏导的基本原理。二阶空间偏导数的离散格式为

$ \begin{array}{*{20}{l}} {D_{xx}^u = \frac{{{\partial ^2}u_{0,0}^0}}{{\partial {x^2}}} = \frac{1}{{\Delta {x^2}}}\sum\limits_{l = 1}^M {\sum\limits_{m = 1}^M {{a_l}} } {a_m}[(u_{m + l - 1,0}^0 - u_{ - m + l,0}^0) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (u_{m - l,0}^0 - u_{ - m - l + 1,0}^0)]} \end{array} $

(2)

$ \begin{array}{*{20}{l}} {D_{zz}^u = \frac{{{\partial ^2}u_{0,0}^0}}{{\partial {z^2}}} = \frac{1}{{\Delta {z^2}}}\sum\limits_{l = 1}^M {\sum\limits_{m = 1}^M {{b_l}} } {b_m}[(u_{0,m + l - 1}^0 - u_{0, - m + l}^0) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (u_{0,m - l}^0 - u_{0, - m - l + 1}^0)]} \end{array} $

(3)

$ \begin{array}{*{20}{l}} {D_{xz}^w = \frac{{{\partial ^2}w_{0,0}^0}}{{\partial x\partial z}} = \frac{1}{{\Delta x\Delta z}}\sum\limits_{l = 1}^M {\sum\limits_{m = 1}^M {{c_l}} } {c_m}[({w_{m - \frac{1}{2},l - \frac{1}{2}}} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} w_{ - m + \frac{1}{2},l - \frac{1}{2}}^0 - (w_{m - \frac{1}{2},l + \frac{1}{2}}^0 - {w_{ - m + \frac{1}{2}, - l + \frac{1}{2}}})]} \end{array} $

(4)

式中：M为空间差分算子长度的一半；u_{l, m}⁰为波场u在空间网格点(x=lΔx，z=mΔz)处0时刻的值，Δx、Δz为空间采样间隔；a=(a₁, a₂, …, a_M)、b=(b₁, b₂, …, b_M)、c=(c₁, c₂, …, c_M)为差分系数。

时间二阶导数的离散格式为

$ D_{tt}^u = \frac{{{\partial ^2}u_{0,0}^0}}{{\partial {t^2}}} = \frac{{u_{0,0}^1 - 2u_{0,0}^0 + u_{0,0}^{ - 1}}}{{\Delta {t^2}}} $

(5)

式中Δt为时间采样间隔。同理，可得其他时间和空间二阶导数D_tt^w、D_zz^w、D_xx^w、D_xz^u。值得注意的是，w分量定义在错开的半个网格点上。式(1)的等效交错网格有限差分离散格式为

$ \left\{ {\begin{array}{*{20}{l}} {D_{tt}^u = v_{\rm{P}}^2D_{xx}^u + v_{\rm{S}}^2D_{zz}^u + (v_{\rm{P}}^2 - v_{\rm{S}}^2)D_{xz}^w}\\ {D_{tt}^w = v_{\rm{P}}^2D_{zz}^w + v_{\rm{S}}^2D_{xx}^w + (v_{\rm{P}}^2 - v_{\rm{S}}^2)D_{xz}^u} \end{array}} \right. $

(6)

1.2 基于平面波优化的差分算子求取

若采用泰勒近似法求解差分系数，则a=b=c，均为传统一阶交错网格的差分系数，可由下式求解

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{1^0}}&{{3^0}}& \cdots &{{{(2M - 1)}^0}}\\ {{1^2}}&{{3^2}}& \cdots &{{{(2M - 1)}^2}}\\ {}&{}& \vdots &{}\\ {{1^{2M - 2}}}&{{3^{2M - 2}}}& \cdots &{{{(2M - 1)}^{2M - 2}}} \end{array}} \right] \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{c}} {1{a_1}}\\ {3{a_2}}\\ \vdots \\ {(2M - 1){a_M}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1\\ 0\\ \vdots \\ 0 \end{array}} \right] \end{array} $

(7)

实际上，泰勒近似的方法所确定的差分系数并非针对方程本身，而只是近似空间偏导数，然后代入波动方程求解。目前有不少的优化方法是针对方程本身构建频散关系，从而优化差分系数，然而它们大多基于声波方程^[36]，然后直接过渡到弹性波方程^[37-38]，造成了不必要的精度损失。

由式(1)可以看出，对于弹性波方程，空间偏导数与速度项耦合在一起，与声波方程不同的是，这些速度项不能被单独分离，因此用有限差分去近似空间偏导数时，P、S波速度影响了网格点权值的分配。为解决这一问题，本文定义a≠b≠c，通过构建三种时空频散关系求解这三组差分系数。此外，声波方程的平面波解并不适用于弹性波方程，要想获取更为准确的差分系数，必须从弹性波方程本身出发，通过矩阵分解的方式求取式(1)的平面波解，然后构建新的时空频散关系，从而获得优化的差分算子。

1.2.1 弹性波的平面波解

二阶位移弹性波方程在波数域可表示为

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}\tilde u}}{{\partial {t^2}}} = - [v_{\rm{P}}^2k_x^2\tilde u + v_{\rm{S}}^2k_z^2\tilde u + (v_{\rm{P}}^2 - v_{\rm{S}}^2){k_x}{k_z}\tilde w]}\\ {\frac{{{\partial ^2}\tilde w}}{{\partial {t^2}}} = - [v_{\rm{P}}^2k_z^2\tilde w + v_{\rm{S}}^2k_x^2\tilde w + (v_{\rm{P}}^2 - v_{\rm{S}}^2){k_x}{k_z}\tilde u]} \end{array}} \right. $

(8)

式中(k_x, k_z)=(kcosθ, ksinθ)=k为波数矢量，其中θ为波的传播方向与水平面的夹角。将上式写成矩阵形式

$ \left[ {\begin{array}{*{20}{c}} {\frac{{{\partial ^2}\tilde u}}{{\partial {t^2}}}}\\ {\frac{{{\partial ^2}\tilde w}}{{\partial {t^2}}}} \end{array}} \right] = - \left[ {\begin{array}{*{20}{c}} {v_{\rm{P}}^2k_x^2 + v_{\rm{S}}^2k_z^2}&{(v_{\rm{P}}^2 - v_{\rm{S}}^2){k_x}{k_z}}\\ {(v_{\rm{P}}^2 - v_{\rm{S}}^2){k_x}{k_z}}&{v_{\rm{S}}^2k_x^2 + v_{\rm{P}}^2k_z^2} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\tilde u}\\ {\tilde W} \end{array}} \right] $

(9)

上式右端系数为实对称矩阵，根据实对称矩阵的性质可写为

$ \left[ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}\tilde u}}{{\partial {t^2}}}}\\ {\frac{{{\partial ^2}\tilde w}}{{\partial {t^2}}}} \end{array}} \right] = \mathit{\boldsymbol{P}}\left[ {\begin{array}{*{20}{l}} { - v_{\rm{P}}^2k}&{}\\ {}&{ - v_{\rm{S}}^2k} \end{array}} \right]{\mathit{\boldsymbol{P}}^{ - 1}}\left[ {\begin{array}{*{20}{l}} {\tilde u}\\ {\tilde \varpi } \end{array}} \right] $

(10)

式中

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{P}} = \left[ {\begin{array}{*{20}{l}} {\frac{{{k_x}}}{k}}&{ - \frac{{{k_z}}}{k}}\\ {\frac{{{k_z}}}{k}}&{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{k_x}}}{k}} \end{array}} \right]}\\ {{\mathit{\boldsymbol{P}}^{ - 1}} = {\mathit{\boldsymbol{P}}^{\rm{T}}}} \end{array}} \right. $

(11)

式(10)的通解可以表示为

$ \begin{array}{l} \left[ {\begin{array}{*{20}{l}} {\tilde u(\mathit{\boldsymbol{k}},t)}\\ {\tilde w(\mathit{\boldsymbol{k}},t)} \end{array}} \right] = \mathit{\boldsymbol{P}}\left[ {\begin{array}{*{20}{l}} {{\rm{exp}}( \pm {\rm{i}}{v_{\rm{P}}}kt)}&{}\\ {}&{{\rm{exp}}( \pm {\rm{i}}{v_{\rm{S}}}kt)} \end{array}} \right] \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{P}}^{ - 1}}\left[ {\begin{array}{*{20}{l}} {\tilde u(\mathit{\boldsymbol{k}},0)}\\ {\tilde w(\mathit{\boldsymbol{k}},0)} \end{array}} \right] \end{array} $

(12)

式中$\tilde u$(k, 0)、$\tilde w$(k, 0)为初始值。将上式展开，可得弹性波动方程的平面波解

$ \left\{ \begin{array}{l} \tilde u(\mathit{\boldsymbol{k}},t)\\ \begin{array}{*{20}{c}} { = \frac{1}{{{k^2}}}[k_x^2{\rm{exp}}( \pm {\rm{i}}{v_{\rm{P}}}kt) + k_z^2{\rm{exp}}( \pm {\rm{i}}{v_{\rm{S}}}kt)]\tilde u(\mathit{\boldsymbol{k}},0) + }\\ {\frac{1}{{{k^2}}}[{\rm{exp}}( \pm {\rm{i}}{v_{\rm{P}}}kt) - {\rm{exp}}( \pm {\rm{i}}{v_{\rm{S}}}kt)]{k_x}{k_z}\tilde w(\mathit{\boldsymbol{k}},0)} \end{array}\\ \tilde w(\mathit{\boldsymbol{k}},t)\\ \begin{array}{*{20}{c}} { = \frac{1}{{{k^2}}}[k_z^2{\rm{exp}}( \pm {\rm{i}}{v_{\rm{P}}}kt) + k_x^2{\rm{exp}}( \pm {\rm{i}}{v_{\rm{S}}}kt)]\tilde w(\mathit{\boldsymbol{k}},0) + }\\ {\frac{1}{{{k^2}}}[{\rm{exp}}( \pm {\rm{i}}{v_{\rm{P}}}kt) - {\rm{exp}}( \pm {\rm{i}}{v_{\rm{S}}}kt)]{k_z}{k_x}\tilde u(\mathit{\boldsymbol{k}},0)} \end{array} \end{array} \right. $

(13)

1.2.2 时空频散关系

对式(6)左端项做时间Fourier变换，然后将式(13)代入，可得

$ \begin{array}{*{20}{l}} {{\rm{FT}}[D_{tt}^u] = \frac{2}{{{k^2}\Delta {t^2}}}\{ k_x^2[{\rm{cos}}({v_{\rm{P}}}k\Delta t) - 1] + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k_z^2[{\rm{cos}}({v_{\rm{S}}}k\Delta t) - 1]\} \tilde u(\mathit{\boldsymbol{k}},t) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{2{k_x}{k_z}}}{{{k^2}\Delta {t^2}}}[{\rm{cos}}({v_{\rm{P}}}k\Delta t) - {\rm{cos}}({v_{\rm{S}}}k\Delta t)]\tilde w(\mathit{\boldsymbol{k}},t)} \end{array} $

(14)

$ \begin{array}{*{20}{l}} {{\rm{FT}}[D_{tt}^w] = \frac{2}{{{k^2}\Delta {t^2}}}\{ k_z^2[{\rm{cos}}({v_{\rm{P}}}k\Delta t) - 1] + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k_x^2[{\rm{cos}}({v_{\rm{S}}}k\Delta t) - 1]\} \tilde w(\mathit{\boldsymbol{k}},t) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{2{k_x}{k_z}}}{{{k^2}\Delta {t^2}}}[{\rm{cos}}({v_{\rm{P}}}k\Delta t) - {\rm{cos}}({v_{\rm{S}}}k\Delta t)]\tilde u(\mathit{\boldsymbol{k}},t)} \end{array} $

(15)

同理，对式(6)右端项做空间Fourier变换，可得

$ {\rm{FT}}[D_{xx}^u] = - {\left\{ {\sum\limits_{m = 1}^M 2 {a_m}{\rm{sin}}[(m - 0.5){k_x}h]} \right\}^2}\tilde u(\mathit{\boldsymbol{k}},t) $

(16)

$ {\rm{FT}}[D_{zz}^u] = - {\left\{ {\sum\limits_{m = 1}^M 2 {b_m}{\rm{sin}}[(m - 0.5){k_z}h]} \right\}^2}\tilde u(\mathit{\boldsymbol{k}},t) $

(17)

$ \begin{array}{*{20}{l}} {{\rm{FT}}[D_{xz}^w] = - \left\{ {\sum\limits_{m = 1}^M 2 {c_m}{\rm{sin}}[(m - 0.5){k_z}h]} \right\} \times }\\ {\quad {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\sum\limits_{m = 1}^M 2 {c_m}{\rm{sin}}[(m - 0.5){k_z}h]} \right\}\tilde w(\mathit{\boldsymbol{k}},t)} \end{array} $

(18)

$ {\rm{FT}}[D_{zz}^w] = - {\left\{ {\sum\limits_{m = 1}^M 2 {a_m}{\rm{sin}}[(m - 0.5){k_z}h]} \right\}^2}\tilde w(\mathit{\boldsymbol{k}},t) $

(19)

$ FT[D_{xx}^w] = - \left\{ {\sum\limits_{m = 1}^M 2 {b_m}{\rm{sin}}{{[(m - 0.5){k_x}h]}^2}} \right\}\tilde w(\mathit{\boldsymbol{k}},t) $

(20)

$ \begin{array}{*{20}{l}} {{\rm{FT}}[D_{xz}^u] = - \left\{ {\sum\limits_{m = 1}^M 2 {c_m}{\rm{sin}}[(m - 0.5){k_x}h]} \right\} \times }\\ {\quad {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\sum\limits_{m = 1}^M 2 {c_m}{\rm{sin}}[(m - 0.5){k_z}h]} \right\}\tilde u(\mathit{\boldsymbol{k}},t)} \end{array} $

(21)

式中h=Δx=Δz为空间采样间隔。

联立式(14)~式(21)，化简可得P、S和PS波时空频散关系

$ \begin{array}{l} {\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)\\ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{{v_{\rm{P}}^2\Delta {t^2}}}{{{h^2}}}\left\{ {{{\left[ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}((m - 0.5){k_x}h)} \right]}^2} + } \right.}\\ {\left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left[ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}((m - 0.5){k_z}h)} \right]}^2}} \right\}} \end{array} \end{array} $

(22)

$ \begin{array}{l} {\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{S}}}k\Delta t)\\ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{{v_{\rm{S}}^2\Delta {t^2}}}{{{h^2}}}\left\{ {{{\left[ {\sum\limits_{m = 1}^M {{b_m}} {\rm{sin}}((m - 0.5){k_x}h)} \right]}^2} + } \right.}\\ {\left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left[ {\sum\limits_{m = 1}^M {{b_m}} {\rm{sin}}((m - 0.5){k_z}h)} \right]}^2}} \right\}} \end{array} \end{array} $

(23)

$ \begin{array}{l} \frac{{{k_x}{k_z}}}{{{k^2}}}[{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t) - {\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{S}}}k\Delta t)]\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{{(v_{\rm{P}}^2 - v_{\rm{S}}^2)\Delta {t^2}}}{{{h^2}}}\left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}[(m - 0.5){k_x}h]} \right\} \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}[(m - 0.5){k_z}h]} \right\} \end{array} $

(24)

对于式(22)~式(24)，本文采用牛顿迭代法使其误差最小化，以求解差分系数a、b、c。

1.2.3 牛顿迭代法

式(22)~式(24)是非线性的，无法直接求解，只能依靠迭代逼近方法求解。首先以式(22)为例，利用二阶范数构建目标泛函

$ \varPhi (\mathit{\boldsymbol{a}}) = \frac{1}{2}\sum\limits_{k = 0}^{{k_e}} {\sum\limits_{\theta = 0}^{\pi /4} {\left\| {\frac{{r_1^2q(\mathit{\boldsymbol{a}})}}{{{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}} - 1} \right\|_2^2} } $

(25)

式中：k_e=π/h为最大优化波数；r₁=v_PΔt/h；

$ \begin{array}{*{20}{l}} {q(\mathit{\boldsymbol{a}}) = {{\left\{ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}[(m - 0.5){k_x}h]} \right\}}^2} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left\{ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}[(m - 0.5){k_z}h]} \right\}}^2}} \end{array} $

(26)

式(25)的牛顿迭代法解可写为

$ {\mathit{\boldsymbol{a}}^{i + 1}} = {\mathit{\boldsymbol{a}}^i} - {({\mathit{\boldsymbol{H}}^i})^{ - 1}}{\mathit{\boldsymbol{g}}^i} $

(27)

式中：i为迭代次数；H为Hessian矩阵；g为梯度向量。g和H可通过对式(25)分别求取一阶和二阶导数获得

$ \begin{array}{l} \frac{{\partial \varPhi (\mathit{\boldsymbol{a}})}}{{\partial {a_m}}} = \sum\limits_{k = 0}^{{k_{\rm{e}}}} {\mathop \sum \limits_{\theta = 0}^{\pi /4} } \left[ {\frac{{r_1^2q(\mathit{\boldsymbol{a}})}}{{{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}} - 1} \right] \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{r_1^2}}{{{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}}\frac{{\partial q(\mathit{\boldsymbol{a}})}}{{\partial {a_m}}} \end{array} $

(28)

$ \begin{array}{l} \frac{{\partial \varPhi (\mathit{\boldsymbol{a}})}}{{\partial {a_m}\partial {a_n}}} = \sum\limits_{k = 0}^{{k_e}} {\sum\limits_{\theta = 0}^{\pi /4} {\left\{ {\left[ {\frac{{r_1^2q(a)}}{{{{\sin }^2}(0.5{v_P}k\Delta t)}} - 1} \right] \times } \right.} } \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {\frac{{r_1^2}}{{{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}}\frac{{{\partial ^2}q(\mathit{\boldsymbol{a}})}}{{\partial {a_m}\partial {a_n}}} + }\\ {\left. {{{\left[ {\frac{{r_1^2}}{{{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}}} \right]}^2}\frac{{\partial q(\mathit{\boldsymbol{a}})}}{{\partial {a_n}}}\frac{{\partial q(\mathit{\boldsymbol{a}})}}{{\partial {a_m}}}} \right\}} \end{array} \end{array} $

(29)

为了加快收敛速度，将迭代初始值设为基于泰勒展开法的一阶交错网格差分系数，误差阈值设为1%，由于牛顿法的快速收敛性，通常只需要迭代几次即可收敛。因此，就可利用式(27)获得P波项的差分系数a。按照同样的方法，将v_P换成v_S即可获得差分系数b。

当求解频散方程式(24)时，如果直接将左端的sin²(0.5v_PkΔt)-sin²(0.5v_SkΔt)除到方程的右端，则可能会使分母为零，进而导致迭代算法不收敛。为了解决这一问题，可将sin²(0.5v_PkΔt)-sin²(0.5v_SkΔt)拆成两项，然后分别与P波速度和S波速度匹配。构建的目标函数为

$ \begin{array}{*{20}{l}} {\varPhi (\mathit{\boldsymbol{c}}) = \frac{1}{2}\sum\limits_{k = 0}^{{k_e}} {\sum\limits_{\theta = 0}^{\pi /4} {\left\| {\frac{{v_{\rm{P}}^2{k^2}\Delta {t^2}q(\mathit{\boldsymbol{c}})}}{{{h^2}{k_x}{k_z}{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{P}}}k\Delta t)}} - 1} \right\|_2^2 + } } }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}\sum\limits_{k = 0}^{{k_e}} {\sum\limits_{\theta = 0}^{\pi /4} {\left\| {\frac{{v_{\rm{S}}^2{k^2}\Delta {t^2}q(\mathit{\boldsymbol{c}})}}{{{h^2}{k_x}{k_z}{\rm{si}}{{\rm{n}}^2}(0.5{v_{\rm{S}}}k\Delta t)}} - 1} \right\|_2^2} } } \end{array} $

(30)

式中

$ \begin{array}{*{20}{l}} {q(\mathit{\boldsymbol{c}}) = \left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}[(m - 0.5){k_x}h]} \right\} \times }\\ {\quad {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}[(m - 0.5){k_z}h]} \right\}} \end{array} $

(31)

按照求解差分系数a的方法，以式(30)为目标泛函，可求得差分系数c。

此外，通过特征值分析方法^[39]，可求得本文方法的稳定性条件为

$ \begin{array}{l} {\rm{max}}\left[ {{v_{{\rm{max}}}}\Delta t\sqrt {\frac{1}{{\Delta {x^2}}}\sum\limits_{m = 1}^M | {a_m}{|^2} + \frac{1}{{\Delta {z^2}}}\sum\limits_{m = 1}^M | {b_m}{|^2}} ,} \right.\\ {\kern 1pt} \left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {v_{\max }}\Delta t\sqrt {\left( {\frac{1}{{\Delta {x^2}}} + \frac{1}{{\Delta {z^2}}}} \right)\sum\limits_{m = 1}^M | {c_m}{|^2}} } \right] \le 1 \end{array} $

(32)

1.3 频散分析

定义P波、S波和PS波相速度误差分别为

$ {E_{\rm{P}}} = \frac{{v_{\rm{P}}^2\Delta {t^2}s(\mathit{\boldsymbol{a}})}}{{{h^2}{\rm{si}}{{\rm{n}}^2}(\pi f\Delta t)}} $

(33)

$ {{E_{\rm{S}}} = \frac{{v_{\rm{S}}^2\Delta {t^2}s(\mathit{\boldsymbol{b}})}}{{{h^2}{\rm{si}}{{\rm{n}}^2}(\pi f\Delta t)}}} $

(34)

$ {{E_{{\rm{PS}}}} = \frac{{v_{\rm{P}}^2\Delta {t^2}{s_1}(c) + v_{\rm{S}}^2\Delta {t^2}{s_2}(\mathit{\boldsymbol{c}})}}{{{h^2}{\rm{sin}}2\theta {\kern 1pt} {\rm{si}}{{\rm{n}}^2}(\pi f\Delta t)}}} $

(35)

式中

$ \begin{array}{*{20}{l}} {s(\mathit{\boldsymbol{a}}) = {{\left\{ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{P}}}}}{\rm{cos}}\theta } \right]} \right\}}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left\{ {\sum\limits_{m = 1}^M {{a_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{P}}}}}{\rm{sin}}\theta } \right]} \right\}}^2}} \end{array} $

(36)

$ \begin{array}{*{20}{l}} {s(\mathit{\boldsymbol{b}}) = {{\left\{ {\sum\limits_{m = 1}^M {{b_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{S}}}}}{\rm{cos}}\theta } \right]} \right\}}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left\{ {\sum\limits_{m = 1}^M {{b_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{S}}}}}{\rm{sin}}\theta } \right]} \right\}}^2}} \end{array} $

(37)

$ \begin{array}{*{20}{l}} {{s_1}(\mathit{\boldsymbol{c}}) = {{\left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{P}}}}}{\rm{cos}}\theta } \right]} \right\}}^2} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{P}}}}}{\rm{sin}}\theta } \right]} \right\}}^2}} \end{array} $

(38)

$ \begin{array}{*{20}{l}} {{s_2}(\mathit{\boldsymbol{c}}) = {{\left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{S}}}}}{\rm{cos}}\theta } \right]} \right\}}^2} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left\{ {\sum\limits_{m = 1}^M {{c_m}} {\rm{sin}}\left[ {\frac{{2\pi (m - 0.5)hf}}{{{v_{\rm{S}}}}}{\rm{sin}}\theta } \right]} \right\}}^2}} \end{array} $

(39)

式(33)~式(35)中利用关系k=2πf/v将自变量设为频率f，根据采样定理可知f∈[0，v_S/(2h)]。

图 1为均匀低速介质中传统方法和本文方法的三种波的频散误差曲线，其中模拟参数为：v_P=2500m/s、v_S=1443m/s、h=15m、M=5、Δt=0.5ms。可以发现，三种波的频散误差曲线形态不同，间接地证明了对于弹性波方程，有必要采取三组不同的差分系数单独地分析这三种波的频散。由图 1a可以看出：随着频率的增大，传统方法的频散误差值逐渐大于1，说明出现了时间频散，而本文的方法能够在一定程度上压制时间频散，但由于平均效应，会损失一定的空间精度；随着传播角度的增大，这种空间精度的损失会越来越小。由图 1b可以看出，S波产生了严重的空间频散，本文方法对S波的频散有很好的压制效果。图 1c为PS波频散误差曲线，其频散误差介于P波与S波之间，与传统方法相比，本文方法也能有效压制PS波频散。

为了验证本文方法在高速介质中的优势，将速度增大为：v_P=4000m/s、v_S=2309m/s，其余参数不变，计算三种波的频散曲线(图 2)。可以看出，在高速介质中，三种波都出现了很严重的时间频散，P波时间频散最为严重，随着频率的增大，频散误差越来越大。与传统方法相比，本文方法能够有效压制时间频散，尤其是在高频段。另外，对比图 2与图 1可以看出，图 2的有效频带范围要宽于图 1(红圈所示)，说明随着速度的增大，有效频带范围变宽。

2 弹性波逆时偏移成像策略 2.1 矢量波分离和成像条件

由于弹性波中的P波和S波是耦合在一起的，在成像的过程会存在串扰假象，严重影响成像剖面的质量，因此在成像之前需要对弹性波进行分离。为了便于推导，将式(1)写成矢量形式

$ \frac{{{\partial ^2}\mathit{\boldsymbol{u}}}}{{\partial {t^2}}} = v_{\rm{P}}^2\nabla (\nabla \cdot \mathit{\boldsymbol{u}}) + v_{\rm{S}}^2\nabla \times \nabla \times \mathit{\boldsymbol{u}} $

(40)

式中u=(u, w)表示位移矢量。P波与S波的耦合关系可以表示为

$ \mathit{\boldsymbol{u}} = \mathit{\boldsymbol{p}} + \mathit{\boldsymbol{s}} $

(41)

式中：p表示P波分量；s为S波分量。将式(41)代入式(40)，有

$ \frac{{{\partial ^2}\mathit{\boldsymbol{p}}}}{{\partial {t^2}}} + \frac{{{\partial ^2}\mathit{\boldsymbol{s}}}}{{\partial {t^2}}} = v_{\rm{P}}^2\nabla (\nabla \cdot \mathit{\boldsymbol{u}}) + v_{\rm{S}}^2\nabla \times \nabla \times \mathit{\boldsymbol{u}} $

(42)

由于P波是无旋场，S波是无散场，因此对式(42)两边同时求散度和旋度，可得

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}\mathit{\boldsymbol{p}}}}{{\partial {t^2}}} = v_{\rm{P}}^2\nabla (\nabla \cdot \mathit{\boldsymbol{u}})}\\ {\frac{{{\partial ^2}\mathit{\boldsymbol{s}}}}{{\partial {t^2}}} = v_{\rm{S}}^2\nabla \times \nabla \times \mathit{\boldsymbol{u}}} \end{array}} \right. $

(43)

令u=p_x+s_x，w=p_z+s_z，将式(43)展开为笛卡尔坐标(x，z)的形式，有

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}{p_x}}}{{\partial {t^2}}} = v_{\rm{P}}^2\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}w}}{{\partial x\partial z}}} \right)}\\ {\frac{{{\partial ^2}{p_z}}}{{\partial {t^2}}} = v_{\rm{P}}^2\left( {\frac{{{\partial ^2}w}}{{\partial {z^2}}} + \frac{{{\partial ^2}u}}{{\partial x\partial z}}} \right)}\\ {\frac{{{\partial ^2}{s_x}}}{{\partial {t^2}}} = v_{\rm{S}}^2\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}} - \frac{{{\partial ^2}w}}{{\partial x\partial z}}} \right)}\\ {\frac{{{\partial ^2}{s_z}}}{{\partial {t^2}}} = v_{\rm{S}}^2\left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}} - \frac{{{\partial ^2}u}}{{\partial x\partial z}}} \right)} \end{array}} \right. $

(44)

基于上式，就可以利用有限差分法外推P波和S波两个方向的波场，再利用内积成像公式提取PP波、PS波成像值

$ \left\{ \begin{array}{l} {I_{{\rm{PP}}}}({\rm{ }}\mathit{\boldsymbol{x}}{\rm{ }}){\rm{ }} = \int {[{p_x}({\mathit{\boldsymbol{x}}_{\rm{s}}},t){p_x}({\mathit{\boldsymbol{x}}_{\rm{r}}},t) + {p_z}({\mathit{\boldsymbol{x}}_{\rm{s}}},t){p_z}({\mathit{\boldsymbol{x}}_{\rm{r}}},t)]{\rm{d}}t} \\ {I_{{\rm{PS}}}}({\rm{ }}\mathit{\boldsymbol{x}}{\rm{ }}){\rm{ }} = \int {[{p_x}({\mathit{\boldsymbol{x}}_{\rm{s}}},t){s_x}({\mathit{\boldsymbol{x}}_{\rm{r}}},t) + {p_z}({\mathit{\boldsymbol{x}}_{\rm{s}}},t){s_z}({\mathit{\boldsymbol{x}}_{\rm{r}}},t)]{\rm{d}}t} {\rm{ }} \end{array} \right. $

(45)

式中：x_s为炮点坐标；x_r为检波点坐标。由于地震勘探主要使用P波震源，本文没有给出SP波和SS波成像过程。

2.2 有效边界存储策略

逆时偏移技术在求解波动方程的时候需要存储全部时刻的震源波场和接收波场，需大量的内存而难以处理巨量数据。将平衡计算量和存储量的有效边界存储策略引入到弹性波逆时偏移中，可有效节约存储空间，如图 3所示。值得注意的是，为了保证空间差分算子长度的边界区域波场值的正确性，需要进行两次波场反传，因此会增加部分计算量。

边界存储的基本做法是：先进行一次震源波场的正向延拓，但只存储边界区域波场，由于接收波场的重构是震源波场的反过程，因此在反传时，可以将边界区域的波场作为边界条件重构之前时刻的波场，这样会增加一些计算量，但大大地降低了存储量，而且边界区域的波场无需重新计算，直接取出替换即可，有效缓解了计算负担。

研究发现，要重构反传波场，无需存储全部边界区域的波场，只需要保证部分边界区域的波场值是正确的，就可以进行波场反向延拓，即有效边界存储策略，该区域的厚度与有限差分空间阶数有关。如采用10阶空间有限差分，根据传统有限差分法的延拓公式，只需要保证边界区域厚度为10以及计算区域边界厚度为10的区域波场即可，如图 3黄色区域所示。

如果不采用边界存储策略，弹性波逆时偏移需要存储量为

$ {S_1} = ({N_x} + 2 \times {N_{\rm{b}}}) \times ({N_z} + 2 \times {N_{\rm{b}}}) \times {N_t} \times 4{\rm{B}} $

(46)

如果采用边界存储策略，需要的存储量为

$ {S_2} = 2 \times {N_{\rm{b}}} \times ({N_x} + {N_z} + 2 \times {N_{\rm{b}}}) \times {N_t} \times 4{\rm{B}} $

(47)

如采用有效边界存储策略，需要的存储量为

$ {S_3} = 4 \times N \times ({N_x} + {N_z}) \times {N_t} \times 4{\rm{B}} $

(48)

式中：N_x、N_z分别为横纵向网格点数；N_t为时间采样点数；N_b为边界层厚度；N=2M为有限差分空间阶数。假设模型尺寸为201×201，时间样点数为1000，边界层厚度为30，有限差分阶次为10，则三种方法所需的存储量分别为0.254、0.103和0.060GB。当模型更大时，有效边界存储的优势会更加明显。

3 模型试算 3.1 洼陷模型

首先选择比较简单的洼陷模型进行测试，详细讨论本文方法在计算量、内存需求和模拟精度三个方面的优劣。模型尺寸为301×301，时间采样间隔为0.5ms，空间采样间隔为15m，总接收时长为3.0s，采用全接收方式，布置14炮，炮间距为300m，道间距为15m，纵横波速度比为1.73，P波速度模型如图 4所示。选用PML边界条件，厚度为30个空间采样间隔。采用主频为30Hz的Ricker子波进行波场模拟。

图 5展示了传统一阶交错网格法(简称传统方法)与本文平面波优化的等效交错网格法(简称本文方法)水平分量单炮地震记录，有限差分空间阶数均为10，可见本文方法能有效压制数值频散(图 5红框所示)。

图 6a、图 6b分别为传统方法和本文方法的PP波成像剖面，可以看出，本文方法PP波成像剖面的偏移噪声能量更弱(黑色箭头所示)，且对同相轴的刻画更精细(红框所示)。图 6c、图 6d为两种方法的PS波成像剖面，实施了与图 6a、图 6b相同的增益和滤波处理，与PP波成像结果相比，PS波成像剖面存在更为严重的噪声和假象，这与频散分析结论一致，但PS波成像分辨率要高于PP波成像。对比图 6c和图 6d可见，本文方法的PS成像剖面质量要优于传统方法(黑色箭头和红圈所示)。

为了验证本文方法利用低阶的空间差分能实现传统方法高阶差分的成像效果，抽取两种方法、不同差分阶次的PP波(图 7)、PS波(图 8)成像单道进行对比，并计算与参考道的残差

$ E = \sum\limits_{n = 1}^{{N_z}} {(|\mathit{\boldsymbol{I}}_n^{{\rm{cal}}} - \mathit{\boldsymbol{I}}_n^{{\rm{ref}}}|)} $

(49)

式中：I^cal为PP波或PS波成像结果；I^ref为参考道，是传统方法差分算子长度为20时的逆时偏移结果(接近无频散)。

由图 7和图 8可见，本文方法6阶空间差分能达到传统方法10阶空间差分的成像效果，而10阶空间差分则能达到传统方法18阶的成像效果，因此在相同的精度要求下，本文方法使用的差分阶次低，能有效节约计算成本。

以有效边界存储策略为例，表 1统计了本文方法和传统方法采用不同阶次有限差分进行洼陷模型逆时偏移的用时，可见：差分算子阶次越高，计算效率越低；如利用6阶本文方法代替10阶传统方法，计算效率可提升约37%；如利用10阶本文方法代替18阶传统方法，计算效率可提升约31%。对比洼陷模型不同阶次、不同方法、采用不同存储策略的内存占用(表 1)，可见有效边界存储策略能节省大量内存，且空间差分阶数越低节省得越多。

3.2 大庆模型

选取大庆油田A区块速度模型(图 9a)验证本文方法对复杂模型的适用性。模型尺寸为1250×248，Δx=12m，Δz=8m，纵横波速度比为1.73。震源选取主频为30Hz的Ricker子波，Δt=0.5ms，有限差分算子长度为10。观测系统采用全接收的方式布置51炮，炮间距为240m，道间距为12m，接收时长为4.0s。

图 9b和图 9c分别为传统方法和本文方法模拟的单炮记录，图 10a和图 10b分别为图 9b和图 9c红框中的局部放大。对比图 10a和图 10b红色箭头处可以看出，传统方法产生了严重的数值频散，本文方法利用更优的差分算子，能有效压制数值频散。

图 11为大庆模型本文方法弹性波逆时偏移结果。可以看出，本文方法能较好地展示成像剖面上的构造层位，对孔洞的几何边界刻画也较为精确，说明了本文方法对复杂模型的有效性。

抽取6、10阶本文方法成像结果与10、18阶传统方法进行单道对比，并计算其残差(图 12、图 13)。可见，6和10阶本文方法成像精度分别与10、18阶传统方法相当，意味着在相同的精度要求下，利用本文方法进行波场外推所需的差分算子更短，因此能有效降低节约计算量。

4 结论

(1) 本文通过推导二阶弹性波方程的平面波解，构建了新的时空频散关系，利用牛顿迭代法，获得了优化的差分系数。平面波优化的差分算子在低速介质中能压制空间频散，在高速介质中能有效压制时间频散。与一阶速度—应力方程相比，避免了中间变量的出现，因此能减少正演和偏移的内存占用。

(2) 将本文提出的平面波优化的差分算子应用于有效边界存储弹性波逆时偏移，有效地降低了存储需求，减小了计算量。

(3) 洼陷模型和复杂的大庆模型弹性波逆时偏移结果表明，6、10阶本文方法能达到10、18阶传统方法的成像精度，因此能有效减少计算量。