0 引言

地震波场数值模拟是复杂地区地震勘探常用手段^[1]，尤其在逆时偏移成像^[2]、全波形反演^[3]等领域。这些模拟方法可分为波动方程法和几何射线法两种主要类型。相较于几何射线法，波动方程法模拟的地震波场包含了更多的地震波传播信息，能为地震波传播机理研究及后期地质综合解释提供更多基础数据。因此，波动方程法在地震波场数值模拟中得到越来越多的应用。

波动方程数值解法主要有有限差分法^[4-5]、伪谱法^[6-7]、有限元法^[8-10]等。其中，基于快速傅里叶变换的伪谱法对于带限波场的正演可达到无限阶差分算子精度^[11]，被广泛应用于复杂介质的波场正演模拟与偏移成像。传统的波动方程求解方法往往采用规则的均匀网格对模型进行离散，以模型的最小速度决定网格尺寸，在高速区易产生“过采样”现象。而变网格方法能根据精度要求，在模型速度较低区域采用细网格、在速度较高区域采用粗网格，以便兼顾波场模拟的效率与精度。

Chen等^[12]提出基于坐标变换的梯形网格方法，使用梯形网格对介质进行剖分，并通过坐标变换在新坐标系下采用伪谱法对弹性波波场做数值模拟。该方法计算效率高，同时避免了一般变网格方法易造成的虚假反射。

高静怀等^[13]最近提出基于深度域均匀采样的梯形网格有限差分模拟法，应用横向采样间隔随深度增加而逐渐线性增大、纵向采样间隔恒定的梯形网格剖分模型，并采用有限差分法求解梯形坐标系波动方程，模拟的地震波场与实际情况相符。但由于在深度域进行均匀采样，对于含有高速夹层的模型，仍会出现一定程度过采样。

Wu等^[14]对上述梯形网格剖分方法^[13]做了改进，在深度方向进行非均匀采样，进一步减少了网格数量。通过Marmousi模型试算显示，梯形网格所占内存只有常规网格的25%，同时计算时间减少了66%，大大提高了计算效率。

本文基于Wu等^[14]的研究，提出一种梯形网格伪谱法地震波模拟方法。内容包括：定义与速度相关的梯形坐标变换；推导梯形坐标系下波动方程表达式，同时为消除人工截断边界反射波的影响，定义了相关的吸收边界；研究时间差分的稳定性条件；利用典型模型测试数值结果，并与常规网格伪谱法和有限差分法结果进行对比；最后给出所得结论。

1 基本理论 1.1 梯形坐标变换

根据Wu等^[14]的理论，可定义如下直角坐标系与梯形坐标系之间的变换关系

$ \left\{ {\begin{array}{*{20}{l}} {x = \frac{{{x_0} - \alpha }}{{1 + {\gamma _{{z_0}}}}}}\\ {{z_0} = g(z)} \end{array}} \right. $

(1)

式中：(x, z)为梯形坐标系坐标，(x₀, z₀)为对应的直角坐标系坐标；α为横向位置参数，一般与震源横向位置相同；γ为形态参数，与模型速度随深度变化的程度相关；g(z)为深度域映射函数。

由式(1)可看出，梯形坐标系中在矩形网格的横向采样间隔确定为Δx时，对应直角坐标系中梯形网格在深度z₀处横向采样间隔为

$ \Delta {x_0}({z_0}) = (1 + \gamma {z_0})\Delta x $

(2)

Δx的求取方法：首先由Δx_0max(z₀)=$\frac{v_{\min }\left(z_{0}\right)}{f_{0} N_{G x}}$确定直角坐标系中梯形网格在深度为z₀时可取的最大横向采样间隔，其中v_min(z₀)为模型在深度为z₀时的最小速度，f₀为震源的主频率，N_Gx为设定的一个波长内最小的横向采样点数。设直角坐标系梯形网格的最大深度的横向覆盖范围与给定模型相同，根据求得的直角坐标系梯形网格最大深度的横向采样间隔即可求出梯形坐标系矩形网格的横向采样点数。然后，通过一个深度相关的函数，可逐步求出直角坐标系梯形网格每一深度的横向采样间隔，取其中最小值即为梯形坐标系矩形网格的横向采样间隔。最后，由已得到的梯形坐标系矩形网格的横向采样点数和采样间隔判断此时梯形坐标系矩形网格的横向覆盖范围是否达到原始模型横向覆盖范围的50%，若不满足，则可通过增加梯形坐标系矩形网格的横向采样点数增大梯形坐标系矩形网格的横向覆盖范围。

得到Δx后，设梯形坐标系矩形网格的纵向采样间隔Δz=Δx，则对于深度采样函数g(z)，可用以下迭代方法得到

$ \left\{ {\begin{array}{*{20}{l}} {g(0) = 0}\\ {g(z + \Delta z) = g(z) + \frac{{{v_{{\rm{min}}}}[g(z)]}}{{{f_0}{N_{Gz}}}}} \end{array}} \right. $

(3)

式中N_Gz为一个波长内的最小纵向采样点数。

通过三次样条插值得到g(z)的一阶、二阶导数g′(z)与g″(z)。最后，可得到形态参数

$ \gamma = \frac{{\frac{{\Delta {x_0}({z_{{\rm{0}}{\kern 1pt} {\rm{max}}}})}}{{\Delta x}} - 1}}{{{z_{{\rm{0}}{\kern 1pt} {\rm{max}}}}}} $

(4)

式中z_{0 max}为给定模型的最大深度。

图 1为梯形网格剖分示意图。

1.2 带PML吸收边界的声波方程的梯形坐标系表征

由式(1)通过链式法则可求得两个坐标系之间一阶、二阶偏导的关系

$ {\frac{\partial }{{\partial {x_0}}} = \frac{1}{{1 + \gamma g(z)}}\frac{\partial }{{\partial x}}} $

(5)

$ {\frac{\partial }{{\partial {z_0}}} = - \frac{{\gamma x}}{{1 + \gamma g(z)}}\frac{\partial }{{\partial x}} + \frac{1}{{{g^\prime }(z)}}\frac{\partial }{{\partial z}}} $

(6)

$ {\frac{{{\partial ^2}}}{{\partial x_0^2}} = \frac{1}{{{{[1 + \gamma g(z)]}^2}}}\frac{{{\partial ^2}}}{{\partial {x^2}}}} $

(7)

$ \begin{array}{l} \frac{{{\partial ^2}}}{{\partial z_0^2}} = \frac{{2{\gamma ^2}x}}{{{{[1 + \gamma g(z)]}^2}}}\frac{\partial }{{\partial x}} + \frac{{{\gamma ^2}{x^2}}}{{{{[1 + \gamma g(z)]}^2}}}\frac{{{\partial ^2}}}{{\partial {x^2}}} - \\ \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{{2\gamma x}}{{[1 + \gamma g(z)]{g^\prime }(z)}}\frac{{{\partial ^2}}}{{\partial x\partial z}} - \frac{{{g^{\prime \prime }}(z)}}{{{{[{g^\prime }(z)]}^3}}}\frac{\partial }{{\partial z}} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\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}{{{{[{g^\prime }(z)]}^2}}}\frac{{{\partial ^2}}}{{\partial {z^2}}}} \end{array} \end{array} $

(8)

直接计算式(8)的空间混合偏导项$\frac{\partial^{2}}{\partial x \partial z}$会产生较大误差。本文通过坐标变换，将混合偏导项转化为非混合偏导项^[15]

$ \left[ {\begin{array}{*{20}{l}} x\\ z \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos \theta }&{ - \cos \theta }\\ {\sin \theta }&{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sin \theta } \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{x^\prime }}\\ {{z^\prime }} \end{array}} \right] $

(9)

式中θ为沿x方向的网格线与对角线的夹角。由于在梯形坐标系网格中Δx=Δz，所以θ=$\frac{\pi}{4}$。

将式(9)代入空间混合偏导项，即得

$ \frac{{{\partial ^2}}}{{\partial x\partial z}} = \frac{1}{{2\sin 2\theta }}\left( {\frac{{{\partial ^2}}}{{\partial {x^{\prime 2}}}} - \frac{{{\partial ^2}}}{{\partial {z^{\prime 2}}}}} \right) = \frac{1}{2}\left( {\frac{{{\partial ^2}}}{{\partial {x^{\prime 2}}}} - \frac{{{\partial ^2}}}{{\partial {z^{\prime 2}}}}} \right) $

(10)

根据完全匹配层吸收边界条件相关理论^[16-19]，可将波场u分裂为三项u₁、u₂、u₃，在计算区域外部引入完全匹配层。直角坐标系中带PML吸收边界的二阶声波方程分以下三种情形表示。

在内部计算区域

$ \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{{v^2}}}\frac{{{\partial ^2}{u_1}}}{{\partial {t^2}}} = \frac{{{\partial ^2}u}}{{\partial x_0^2}}}\\ {\frac{1}{{{v^2}}}\frac{{{\partial ^2}{u_2}}}{{\partial {t^2}}} = \delta ({x_0} - {x_{{\rm{0s}}}})\delta ({z_0} - {z_{{\rm{0s}}}})f(t)}\\ {\frac{1}{{{v^2}}}\frac{{{\partial ^2}{u_3}}}{{\partial {t^2}}} = \frac{{{\partial ^2}u}}{{\partial z_0^2}}} \end{array}} \right. $

(11)

在左侧和右侧PML区域

$ \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{{v^2}}}{{\left[ {\frac{\partial }{{\partial t}} + d({x_0})} \right]}^2}{u_1} = \frac{{{\partial ^2}u}}{{\partial x_0^2}}}\\ {\frac{1}{{{v^2}}}{{\left[ {\frac{\partial }{{\partial t}} + d({x_0})} \right]}^3}{u_2} = - {d^\prime }({x_0})\frac{{\partial u}}{{\partial {x_0}}}}\\ {\frac{1}{{{v^2}}}\frac{{{\partial ^2}{u_3}}}{{\partial {t^2}}} = \frac{{{\partial ^2}u}}{{\partial z_0^2}}} \end{array}} \right. $

(12)

在上侧和下侧PML区域

$ \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{{v^2}}}\frac{{{\partial ^2}{u_1}}}{{\partial {t^2}}} = \frac{{{\partial ^2}u}}{{\partial x_0^2}}}\\ {\frac{1}{{{v^2}}}{{\left[ {\frac{\partial }{{\partial t}} + d({z_0})} \right]}^3}{u_2} = - {d^\prime }({z_0})\frac{{\partial u}}{{\partial {z_0}}}}\\ {\frac{1}{{{v^2}}}{{\left[ {\frac{\partial }{{\partial t}} + d({z_0})} \right]}^2}{u_3} = \frac{{{\partial ^2}u}}{{\partial z_0^2}}} \end{array}} \right. $

(13)

式中：(x_0s, z_0s)为震源在直角坐标系中位置；衰减函数d(i)及其一阶导数d′(i)可表示为

$ {d(i) = \frac{{3v}}{{2L}}{{\left( {\frac{{{s_i}}}{L}} \right)}^2}{\rm{ln}}\frac{1}{R}} $

$ {{d^\prime }(i) = \frac{{3v}}{{{L^3}}}{\rm{ln}}\frac{1}{R}{s_i}} $

其中：i=x₀或z₀；L为PML吸收层的厚度；s_i为PML层内节点到计算区域边界的距离；R为衰减系数，一般取10^-3。

将式(5)~式(8)代入式(11)~式(13)，即得梯形坐标系中带PML吸收边界的三种情形的二阶声波方程表达式。

在内部计算区域

$ \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{{v^2}}}\frac{{{\partial ^2}{u_1}}}{{\partial {t^2}}} = \frac{1}{{{{[1 + \gamma g(z)]}^2}}}\frac{{{\partial ^2}u}}{{\partial {x^2}}}}\\ {\frac{1}{{{v^2}}}\frac{{{\partial ^2}{u_2}}}{{\partial {t^2}}} = \delta \left( {x - {x_{\rm{s}}}} \right)\delta \left( {z - {z_{\rm{s}}}} \right)f(t)}\\ {\frac{1}{{{v^2}}}\frac{{{\partial ^2}{u_3}}}{{\partial {t^2}}} = \frac{{2{\gamma ^2}x}}{{{{[1 + \gamma g(z)]}^2}}}\frac{{\partial u}}{{\partial x}} + \frac{{{\gamma ^2}{x^2}}}{{{{[1 + \gamma g(z)]}^2}}}\frac{{{\partial ^2}u}}{{\partial {x^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} {\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{{\gamma x}}{{[1 + \gamma g(z)]{g^\prime }(z)}}\left( {\frac{{{\partial ^2}u}}{{\partial {x^{\prime 2}}}} - \frac{{{\partial ^2}u}}{{\partial {z^{\prime 2}}}}} \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} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{g^{\prime \prime }}(z)}}{{{{\left[ {{g^\prime }(z)} \right]}^3}}}\frac{{\partial u}}{{\partial z}} + \frac{1}{{{{\left[ {{g^\prime }(z)} \right]}^2}}}\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \end{array}} \right. $

(14)

在左侧和右侧PML区域

$ \left\{ \begin{array}{l} \frac{1}{{{v^2}}}{\left[ {\frac{\partial }{{\partial t}} + d(x)} \right]^2}{u_1} = \frac{1}{{{{[1 + \gamma g(z)]}^2}}}\frac{{{\partial ^2}u}}{{\partial {x^2}}}\\ \frac{1}{{{v^2}}}{\left[ {\frac{\partial }{{\partial t}} + d(x)} \right]^3}{u_2} = - {d^\prime }(x)\frac{1}{{1 + \gamma g(z)}}\frac{{\partial u}}{{\partial x}}\\ \frac{1}{{{v^2}}}\frac{{{\partial ^2}{u_3}}}{{\partial {t^2}}} = \frac{{2{\gamma ^2}x}}{{{{[1 + \gamma g(z)]}^2}}}\frac{{\partial u}}{{\partial x}} + \frac{{{\gamma ^2}{x^2}}}{{{{[1 + \gamma g(z)]}^2}}}\frac{{{\partial ^2}u}}{{\partial {x^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} {\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{{\gamma x}}{{[1 + \gamma g(z)]{g^\prime }(z)}}\left( {\frac{{{\partial ^2}u}}{{\partial {x^{\prime 2}}}} - \frac{{{\partial ^2}u}}{{\partial {z^{\prime 2}}}}} \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} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{g^{\prime \prime }}(z)}}{{{{\left[ {{g^\prime }(z)} \right]}^3}}}\frac{{\partial u}}{{\partial z}} + \frac{1}{{{{\left[ {{g^\prime }(z)} \right]}^2}}}\frac{{{\partial ^2}u}}{{\partial {z^2}}} \end{array} \right. $

(15)

在上侧和下侧PML区域

$ \left\{ \begin{array}{l} \frac{1}{{{v^2}}}\frac{{{\partial ^2}{u_1}}}{{\partial {t^2}}} = \frac{1}{{{{[1 + \gamma g(z)]}^2}}}\frac{{{\partial ^2}u}}{{\partial {x^2}}}\\ \frac{1}{{{v^2}}}{\left[ {\frac{\partial }{{\partial t}} + d(z)} \right]^3}{u_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} = - {d^\prime }(z)\left[ {\frac{1}{{{g^\prime }(z)}}\frac{{\partial u}}{{\partial z}} - \frac{{\gamma x}}{{1 + \gamma g(z)}}\frac{{\partial u}}{{\partial x}}} \right]\\ \frac{1}{{{v^2}}}{\left[ {\frac{\partial }{{\partial t}} + d(z)} \right]^2}{u_3}\\ {\kern 1pt} {\kern 1pt} {\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{\gamma ^2}x}}{{{{[1 + \gamma g(z)]}^2}}}\frac{{\partial u}}{{\partial x}} + \frac{{{\gamma ^2}{x^2}}}{{{{[1 + \gamma g(z)]}^2}}}\frac{{{\partial ^2}u}}{{\partial {x^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} \frac{{\gamma x}}{{[1 + \gamma g(z)]{g^\prime }(z)}}\left( {\frac{{{\partial ^2}u}}{{\partial {x^{\prime 2}}}} - \frac{{{\partial ^2}u}}{{\partial {z^{\prime 2}}}}} \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} \frac{{{g^{\prime \prime }}(z)}}{{{{\left[ {{g^\prime }(z)} \right]}^3}}}\frac{{\partial u}}{{\partial z}} + \frac{1}{{{{\left[ {{g^\prime }(z)} \right]}^2}}}\frac{{{\partial ^2}u}}{{\partial {z^2}}} \end{array} \right. $

(16)

式中(x_s, z_s)为震源在梯形坐标系中位置。

此处的衰减函数d(i)及其一阶导数d′(i)的表达式在形式上与直角坐标系时相同，且i=x或z。

1.3 傅里叶伪谱法原理

由傅里叶变换的相关理论^[20]，若函数f(x)在Hilbert空间中可积，则有

$ \begin{array}{l} \frac{{{{\rm{d}}^m}}}{{{\rm{d}}{x^m}}}f(x) = \frac{1}{{2\pi }}\int_{ - \infty }^\infty {{{({\rm{i}}k)}^m}} F(k){{\rm{e}}^{{\rm{i}}kx}}{\rm{d}}k\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\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 \in \mathscr{R} \end{array} $

(17)

式中F(k)为f(x)的傅里叶变换。

以均匀网格离散f(x)，则式(17)变为

$ \begin{array}{*{20}{l}} {\frac{{{{\rm{d}}^m}}}{{{\rm{d}}{x^m}}}f(n\Delta x) = \frac{1}{{N\Delta x}}\sum\limits_{l = - (\frac{N}{2} - 1)}^{\frac{N}{2}} {{{({\rm{i}}l\Delta k)}^m}} F(l\Delta k){{\rm{e}}^{{\rm{i}}2\pi \frac{{nl}}{N}}}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} n = 0, \cdots ,N - 1} \end{array} $

(18)

式中$\Delta k=\frac{2 \pi}{N \Delta x}$。

通过上述方法就能得到f(x)在各个离散点的导数值。对于波场函数，只需改变相关参数即可计算其偏导值。在数值模拟中常采用快速傅里叶变换加速计算。

1.4 稳定性条件

与z轴夹角为θ的平面波的解析式为

$ u(x,z) = u_0^*{\kern 1pt} {{\rm{e}}^{{\rm{i}}{k_x}x + i{k_z}z - {\rm{i}}\omega t}} $

(19)

不带PML边界的梯形坐标系伪谱法求解声波方程的差分格式为

$ \begin{array}{l} \frac{{u_{m,n}^{j + 1} - 2u_{m,n}^j + u_{m,n}^{j - 1}}}{{{v^2}\Delta {t^2}}} = \frac{{{\gamma ^2}{x^2} + 1}}{{{{[1 + \gamma g(z)]}^2}}}\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}}} \right)_{m,n}^j + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{2{\gamma ^2}x}}{{{{[1 + \gamma g(z)]}^2}}}\left( {\frac{{\partial u}}{{\partial x}}} \right)_{m,n}^j - \frac{{{g^{\prime \prime }}(z)}}{{{{\left[ {{g^\prime }(z)} \right]}^3}}}\left( {\frac{{\partial u}}{{\partial z}}} \right)_{m,n}^j - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{2\gamma x}}{{[1 + \gamma g(z)]{g^\prime }(z)}}\left( {\frac{{{\partial ^2}u}}{{\partial x\partial z}}} \right)_{m,n}^j + \frac{1}{{{{\left[ {{g^\prime }(z)} \right]}^2}}}\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)_{m,n}^j \end{array} $

(20)

将式(19)代入式(20)，可得

$\frac{2}{v^{2} \Delta t^{2}}[1-\cos (\omega \Delta t)]=\frac{\gamma^{2} x^{2}+1}{[1+\gamma g(z)]^{2}} k_{x}^{2}-\\ \;\;\;\;\;\;\frac{2 \gamma^{2} x}{[1+\gamma g(z)]^{2}} {\rm{i}} k_{x}-\frac{2 \gamma x}{[1+\gamma g(z)] g^{\prime}(z)} k_{x} k_{z}+\\ \;\;\;\;\;\;\frac{g^{\prime \prime}(z)}{\left[g^{\prime}(z)\right]^{3}} \mathrm{i} k_{z}+\frac{1}{\left[g^{\prime}(z)\right]^{2}} k_{z}^{2} $

(21)

式中：ω为角频率；k=$\frac{\omega}{v}$为波数，且有k_x=ksinθ，k_z=kcosθ。

由于差分误差会随着波数的增大而增大，因此考虑最大波数的情形，即$k_{x}=\frac{\pi}{\Delta x}$，$k_{z}=\frac{\pi}{\Delta z}$。令

$ \begin{array}{l} {s_1} = \frac{{{\gamma ^2}{x^2} + 1}}{{{{[1 + \gamma g(z)]}^2}}}\frac{{{\pi ^2}}}{{\Delta {x^2}}} - {\rm{i}}\frac{{2{\gamma ^2}x}}{{{{[1 + \gamma g(z)]}^2}}}\frac{\pi }{{\Delta x}} - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\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\gamma x}}{{[1 + \gamma g(z)]{g^\prime }(z)}}\frac{{{\pi ^2}}}{{\Delta x\Delta z}} + {\rm{i}}\frac{{{g^{\prime \prime }}(z)}}{{{{\left[ {{g^\prime }(z)} \right]}^3}}}\frac{\pi }{{\Delta z}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\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}{{{{\left[ {{g^\prime }(z)} \right]}^2}}}\frac{{{\pi ^2}}}{{\Delta {z^2}}} \end{array} $

(22)

$ \cos (\omega \Delta t) \approx 1 - \frac{{{v^2}\Delta {t^2}}}{2}{s_1} $

(23)

再令s=1-$\frac{v^{2} \Delta t^{2}}{2}$s₁，则当|s|＜1时，该差分格式是稳定的，因此可得稳定性条件为

$ \Delta t < \frac{2}{{v\sqrt { {\rm{Re}} ({s_1})} }} $

(24)

2 数值算例

利用不同速度模型脉冲响应验证本文方法的正确性，并与常规网格伪谱法及常规网格有限差分法结果进行对比。测试均基于Intel(R) Core(TM) i5-6300HQ CPU及MATLAB编程。其中常规网格伪谱法与常规网格有限差分法均采用中心差分近似时间偏导，常规网格有限差分法采用八阶差分算子近似空间偏导^[21-22]。由于梯形坐标系伪谱法可在大尺度时间步长下保持良好稳定性，因此采用Erik等^[23-24]提出的时域离散变换消除大尺度时间步长带来的时间频散。

2.1 均匀介质模型

在一个速度为3000m/s、深度为2000m、宽度为1000m的梯形坐标系均匀介质模型中，设置横向、纵向采样间隔均为25m，形态参数γ=5×10^-4。在模型中心放置主频为10Hz的Ricker子波作为震源，震源时延为0.075s。使用时间步长为1ms的梯形坐标系伪谱法做波场模拟。得到t=0.4s时的波场快照(图 2a)。对该波场快照做坐标变换即可得到直角坐标系波场快照。但由于梯形网格剖分的空间采样点数较少，因而得到的直角坐标系波场快照会缺失部分波场信息，可通过对直角坐标系波场快照插值得到更多波场信息(图 2b)。图 3显示由梯形坐标系中位于(500m，0m)处的检波器接收到的信号插值得到的直角坐标系(1000m，0)处的波场时间曲线，可见梯形网格伪谱法求得的数值解与解析解^[25]高度吻合，验证了本文方法的正确性。

2.2 三层速度模型

测试选用图 4所示的2000m×2000m的三层速度模型。各层的速度、厚度(从上到下)依次为2000、4000、3000m/s，500、500、1000m。梯形网格伪谱法的实际模拟范围是图中红线包围的(梯形)区域。设置直角坐标系梯形网格的横向采样间隔由最顶层的10.5470m增至最底层的18.5185m，深度域映射函数如图 5所示，形态参数为γ=3.697×10^-4。

在模型顶层正中央放置主频为20Hz的Ricker子波作为震源，使用时间步长为1ms的梯形坐标系伪谱法做波场模拟，得到t₁=0.18s、t₂=0.42s、t₃=0.81s三个时刻梯形坐标系波场快照(图 6a~图 6c)。经过坐标变换和插值，得到相应的直角坐标系波场快照(图 6d~图 6f)。同时，在相同的三层速度模型中使用时间步长为1ms，横向、纵向采样间隔均为10m的常规网格伪谱法和常规网格有限差分法得到同样三个时刻的波场快照(图 6g~图 6i、图 6j~图 6l)。可见梯形网格因未覆盖左右两个角形区域，使用梯形网格伪谱法求得的波场快照会缺失部分浅层大角度信息，但在其他区域与常规网格方法结果一致。

2.3 Marmousi模型

选取的Marmousi模型(图 7a)中红色线圈定范围是梯形网格伪谱法的实际模拟区域。设置梯形网格的横向采样间隔由最浅层的8.7184m增至最底层的17.3493m，深度域映射函数如图 7c所示，可算出形态参数γ=3.3003×10^-4。在模型顶层正中间放置主频为15Hz的Ricker子波震源，震源时延为0.125s。使用时间步长为0.943ms的梯形网格伪谱法做波场模拟，得到t₁=1s、t₂=2s时的梯形坐标系波场快照(图 8a、图 9a)。

经过坐标变换和插值，得到t₁=1s、t₂=2s时的直角坐标系波场快照(图 8b、图 9b)。在Marmousi原始速度模型中使用时间步长为0.436ms，横向采样间隔为12.5m，纵向采样间隔为4m的常规网格伪谱法做波场模拟，得到t₁=1s、t₂=2s时的波场快照(图 8c、图 9c)。使用时间步长为0.537ms，横向采样间隔为12.5m，纵向采样间隔为4m的常规网格有限差分法进行波场模拟，得到的t₁=1s、t₂=2s时的波场快照(图 8d、图 9d)。可见三种方法在Marmousi模型中正演得到的同一时刻波场快照并无差别。为了对结果做更详细对比，抽取了三种方法得到的t₁=1s时波场快照的第269、第369及第469道的深度域波形(图 10)，可见三者波形高度一致，充分体现了梯形网格伪谱法的高精度。

表 1列出使用三种方法进行地震波数值模拟时各项参数。可看到相较于常规网格剖分，梯形网格剖分的网格数减少了69.26%；相较于常规网格伪谱法，梯形网格伪谱法的计算时间减少了58.16%；相较于常规网格有限差分法，梯形网格伪谱法的计算时间减少了60.33%。考虑到梯形网格伪谱法的实际模拟区域比常规方法减少约25%，因此对于同一模拟区域，梯形网格伪谱法的计算时间比常规方法只减少约45%。同时，三种方法的精度相当。这充分展现了深度非均匀采样梯形网格伪谱法进行地震波模拟时的高效率和高精度。

3 结论

本文提出了一种梯形网格伪谱法地震波模拟方法。相较于常规均匀网格方法，梯形网格更契合地震波传播速度由浅入深逐渐增大的物理现象，在低速区域采用细网格、高速区域采用粗网格，同时在深度域采用非均匀采样，这样可大幅度减少内存需求，提高计算效率。对均匀介质模型的数值测试表明，本文所提方法符合解析解，同时数值频散较小。三层速度模型的测试结果表明，本文方法正演得到的地震波场与常规方法正演得到的地震波场相比，除了部分浅层大角度信息缺失，无其他区别。针对Marmousi模型的测试表明，梯形网格伪谱法在复杂模型中同样可得很好模拟效果。相较于常规网格剖分，梯形网格剖分能减少约69%的采样点数；相较于常规网格伪谱法，梯形网格伪谱法计算时间能减少约58%；相较于常规网格有限差分法，梯形网格伪谱法的计算时间能减少约60%。上述测试结果表明梯形网格伪谱法地震波模拟方法是一种高效、高精度地震波模拟方法，具有较高应用价值。

当然，本文方法也存在一定的局限性：①梯形网格并未覆盖左右两个角形区域，因此使用梯形网格伪谱法做波场模拟会缺失此部分波场信息；②由于采用梯形网格剖分只拟合地震波传播速度由浅入深逐渐增加的整体趋势，因此对于浅层高速或深层低速模型仍会出现“过采样”。