Spectral element-finite difference hybrid methods for seismic wave simulation
-
摘要:
为了充分发挥SEM和FDM各自在数值模拟复杂介质中地震波传播的优点，同时避免各自缺点，本文提出了模拟地震波传播的谱元-有限差分（SEM-FDM）混合方法．该混合方法在起伏地形附近采用谱元法计算地表附近的地震波传播；在远离起伏地形的区域之下，采用有限差分法计算地震波传播；通过构建数据交换区实现两种方法的耦合．结果表明，本文提出的谱元-有限差分混合方法能有效模拟任意非均匀介质中地震波传播，具有对复杂地形的刻画能力，能处理自由地表边界条件，且能适应模型内速度间断面．通过模型试算，并与谱元法进行对比，验证了该混合方法用于地震波模拟的有效性和高精度性．
Abstract:Efficient methods for seismic wave simulations play important roles in seismic waveform inversion. Currently, the spectral element method （SEM） and finite difference method （FDM） are mostly frequently used for simulating seismic waveforms in waveform inversion due to their efficiency. Compared with other methods, such as finite element method and pseudo spectral method, the computational costs of the SEM and FDM are relatively low, which is greatly important for seismic wave modeling in large-scale models and waveform inversion. Each method, the SEM or FDM, has its advantages and disadvantages when simulating seismic waveforms in complex medium. For example, the SEM can simulate seismic waves in arbitrary heterogeneous medium with velocity discontinuities due to its high-accuracy. In addition, the free-surface boundary can be automatically satisfied. However, designing a mesh with good quality is generally time-consuming. In some cases, it is difficult to construct an appropriate mesh for the SEM. The FDM generally uses regular meshes for seismic wave modeling. Therefore, it is extremely convenient to design a mesh for the FDM. However, free-surface boundary condition cannot be automatically satisfied and special treatments are required to incorporate the free-surface boundary condition. Although some special treatments are proposed in the past decades, algorithms for the approximation of free-surface boundary condition commonly have very low accuracy. Thus, it is difficult to accurately simulate surface waves using the FDM. To retain the advantages and avoid the disadvantages of the SEM and FDM, we develop a hybrid method, which is call SEM-FDM. For this method, the propagation of seismic waves near the free-surface boundary is simulated by the SEM, and in the area far away from the free-surface boundary is modeled by the FDM. A layer for data exchanging is used for coupling the two methods. The SEM-FDM hybrid method is efficient for the simulation of seismic waves in arbitrary heterogeneous media. The hybrid method has the properties for handling undulating topography, free-surface boundary condition and velocity discontinuities. Based on a plane wave analysis, we derive the stability conditions of the hybrid method, and present the stability conditions for spatial accuracy with different orders. In order to demonstrate the validity of the SEM-FDM hybrid method in seismic wave modeling, we construct four models for numerical tests. The first model is used to show the accuracy and efficiency of the SEM-FDM hybrid method. By a comparison with analytical solutions and results from the FDM and SEM, the SEM-FDM is indeed very suitable for seismic wave modeling. The second and third models are used to display the ability of the SEM-FDM hybrid method for modeling seismic wave in complex medium. Although relatively sample grids are adopted, the SEM-FDM can still accurately model seismic waves. The last is a really geological model, which is used to exhibit the usefulness of the proposed method in real cases. By a comparison with the results from the SEM, we demonstrate that the SEM-FDM hybrid method is a useful tool in real cases. Although the areas near the free-surface area are discretized by irregular elements, the process for generate irregular elements are not difficult. In our future work, we will use Fourier series to approximate the free surface and automatically generate a mesh based on the Fourier series and control parameters.
-
引言
地震波正演模拟有助于研究地震波传播规律，是研究地球内部结构的基础，同时在资源勘探和地震灾害评估中发挥着重要作用（Tromp et al，2005；Tong et al，2014；Liu et al，2015a；Shen et al，2022；祝贺君等，2023）．就地震波运动方程的直接求解而言，研究人员已发展出多种数值模拟方法，例如：有限差分法（Zhang，Yao，2013）、有限元法（finite difference method，缩写为FDM）（Marfurt et al，1984）、谱元法（spectral-element method，缩写为SEM）（Komatitsch，Vilotte，1998；刘少林等，2021）等，每种方法各有其优缺点，其中有限差分法和谱元法是目前关于地震波运动方程的数值求解使用最广泛的方法（桑莹泉等，2021）．
有限差分法（FDM）采用泰勒（Taylor）级数近似原理，利用网格点上的波场值直接对方程包含的空间微分项离散（Zhou et al，2021），该方法具有数值格式简单、程序易于实现、计算速度快等优点（Zhang，Yao，2013），因此，FDM在不同尺度的地震波场模拟问题中得到了广泛应用（Moczo et al，2007；Zhang et al，2019）．然而，FDM在模拟地震波时存在一些不足，其不足表现为网格划分不灵活、自由边界条件处理困难以及数值频散现象明显等（Liu，2014；Zhang，Yao，2013）．前人提出了多种策略克服FDM存在的不足，例如：曲线网格方法（Rao，Wang，2013）、优化系数方法（Zhang，Yao，2013）、自由地表边界条件处理方法（Lan，Zhang，2011）．虽然这些改进策略对FDM某方面的不足具有一定的改进效果，但其它方面的不足依然存在。例如：曲线网格的FDM难以以高阶差分格式应用于地震波数值模拟；FDM对自由地表边界条件的近似精度低（苏波等，2019）。
谱元法（SEM）基于变分原理求解弱形式的地震波运动方程，其结合了谱方法的高精度性和有限元法的网格灵活性（Komatitsch，Vilotte，1998），且能自然满足自由地表边界条件（刘少林等，2022）。当采用勒让德（Legendre）型正交插值多项式和高斯-勒让德-洛巴托（Gauss-Legendre-Lobatto，缩写为GLL）数值求积公式时，由于数值积分点与插值点重合的特点（Komatitsch，Tromp，2002；刘少林等，2021），SEM得到的质量矩阵为对角矩阵，而对角矩阵的求逆运算计算量较小，因此SEM用于地震波场模拟时具有较高的计算效率（刘少林等，2021）。在波形伴随层析成像之中，越来越多的研究利用SEM计算地震波形（Chen et al，2015；Lei et al，2020）。
虽然谱元法得到了较为广泛的应用，但该方法存在两方面不足。其一，相比于规则网格的FDM，SEM的计算量明显要大（刘少林等，2021）。在物理坐标向局部坐标转换的过程中，任意方向的空间微分算子要转换为局部坐标所有方向微分算子的加权（Komatitsch，Tromp，2002），该过程增加了计算量；其二，当模型内界面复杂时，网格剖分存在一定困难。谱元法采用六面体单元（三维情形）刻画地下介质结构，当模型内部界面复杂时，难以简单快捷构建高质量网格（Shen et al，2022）。
为了实现复杂含起伏地表模型中的地震波高效数值模拟，前人发展了有限元-有限差分混合方法（Galis et al，2008）。在该混合方法中，利用有限元法模拟起伏地表附近地震波场，采用FDM模拟其他区域中的地震波场。由于有限元法计算量较大，为了提高计算效率，有限元法所计算的区域仅占模型的一小部分（Galis et al，2008）。特别需要指出的是，在该混合方法中FEM难以使用高阶格式得到高精度地震波场。
为了实现复杂模型中地震波场的高精度、高效数值模拟，本文发展了SEM-FDM混合方法。该方法在地表附近采用SEM模拟地震波传播，在远离地表区域利用FDM计算地震波场。混合方法吸收了SEM自由地表边界条件自然满足的优点，具有对复杂地形刻画能力；同时，采用针对任意非均匀介质的有限差分格式，模拟地震波在远离地表区域中的传播，在保证算法计算效率的同时避免了生成复杂网格的难点问题。通过模型试算，并与SEM和FDM等方法进行对比，以证实SEM-FDM混合方法是一种高效的地震波模拟方法。
1. SEM-FDM混合方法原理
以2D非均匀各向同性介质中地震波运动方程为例介绍SEM-FDM混合方法，对于3D或各向异性介质的SEM-FDM混合方法与之类似。首先简要介绍SEM，更详细的推导请见刘少林等（2021）。地震波运动方程可写为：
$$ \begin{array}{l}\left\{\begin{array}{l}\rho \ddot{\boldsymbol{u}}＝\nabla ·{\boldsymbol{T}} ＋ {\boldsymbol{f}} \text{，} \\ {\boldsymbol{T}}＝{\boldsymbol{C}}:\nabla {\boldsymbol{u}} \text{，} \end{array}\right.\end{array} $$ (1) 式中$ \boldsymbol{u}＝\left(u\text{，} v\right) $为位移向量，T为应力张量，$ \boldsymbol{f}＝\left({f}_{x}\text{，} {f}_{{\textit{z}}}\right) $，$ \nabla ＝\left(\frac{ {\text{∂}} }{ {\text{∂}} x}\text{，} \frac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\right) $，ρ为介质质量密度，$ \boldsymbol{C} $为弹性张量，$ \boldsymbol{u} $上两点表示时间二阶微分。利用任意测试向量w乘以式（1），并利用分部积分和格林公式，式（1）可变为：
$$ \begin{array}{c}{\int }_{\mathrm{\Omega }}\mathrm{\rho }\boldsymbol{w}· \ddot{\boldsymbol{u}}{\mathrm{d}}^{2}{\boldsymbol{x}}=-{\int }_{\mathrm{\Omega }}\nabla \boldsymbol{w}:\boldsymbol{T}{\mathrm{d}}^{2}{\boldsymbol{x}} ＋ {\int }_{ {\text{∂}} \mathrm{\Omega }}\boldsymbol{w}· \boldsymbol{T}· \mathbf{n}\mathrm{d}{\mathrm{s}} ＋ {\int }_{\mathrm{\Omega }}\boldsymbol{w}· \boldsymbol{f}{\mathrm{d}}^{2}{\boldsymbol{x}} \text{，} \end{array} $$ (2) 式中Ω为研究区域，n为边界外法向。根据自由地表边界条件和自然边界条件，式（2）中关于边界的积分项为零。将模型剖分为Ne个非重叠单元，即有$ \mathrm{\mathit{\Omega}}＝\bigcup_{\mathrm{i}＝1}^{N_{\mathrm{e}}}\mathrm{\mathit{\Omega}}_{\mathrm{i}} $。在每个单元上将物理坐标（x，z）变换成局部坐标（ξ，η），其中ξ和η的取值范围为$ ［ $−1，1$ ］ $。对于每个局部坐标，根据方程$ （ 1-{\mathrm{\alpha }}^{2} ） {\mathrm{L}}_{\mathrm{N}}' （ \mathrm{\alpha } ） ＝0 $（其中LN（α）为$ N $阶勒让德多项式）的零点确定插值点。确定插值点以后，每个单元上利用插值多项式近似$ \boldsymbol{u} $和$ \boldsymbol{w} $，其近似公式为：
$$ \begin{array}{c}\left\{\begin{array}{c}{\boldsymbol{u}}^{\mathrm{e}} （ \xi \text{，} \eta ） {\text{≈}} \sum\limits _{i＝0}^{N}\sum\limits _{j＝0}^{N}{l}_{i} （ \xi ） {l}_{j} （ \eta ） {\boldsymbol{u}}_{ij}^{\mathrm{e}} \text{，} \\ {\boldsymbol{w}}^{\mathrm{e}} （ \xi \text{，} \eta ） {\text{≈}} \sum\limits _{i＝0}^{N}\sum\limits _{j＝0}^{N}{l}_{i} （ \xi ） {l}_{j} （ \eta ） {\boldsymbol{w}}_{ij}^{\mathrm{e}} \text{，} \end{array}\right.\end{array} $$ (3) 式中l为插值多项式，上标e表示单元编号，下标i和j表示插值节点索引。将式（3）代入式（2），在每个单元上采用GLL数值求积公式，式（2）对应的离散形式为：
$$ \begin{array}{c}\sum\limits_{\mathit{\mathit{\mathrm{\mathit{e}}}}＝1}^{N_{\mathrm{e}}}\mathbf{\mathit{\boldsymbol{M}}}^{\mathrm{e}}\ddot{\mathbf{\mathit{\boldsymbol{U}}}}^{\mathrm{e}}＝-\sum\limits_{\mathrm{e}＝1}^{N_e}\mathbf{\mathit{\boldsymbol{K}}}^{\mathrm{e}}\mathbf{\mathit{\mathbf{U}}}^{\mathrm{e}}＋\sum\limits_{\mathrm{e}＝1}^{N_e}\mathbf{\boldsymbol{F}}^{\mathrm{e}}\text{，} \end{array} $$ (4) 式中Ue为离散位移向量，Me为单位质量矩阵，Ke为单元刚度矩阵，Fe为单元震源向量。
为了便于介绍FDM，将式（1）写成如下等价形式：
$$ \begin{array}{c}\left\{\begin{array}{c}\rho \dfrac{{ {\text{∂}} }^{2}u}{ {\text{∂}} {t}^{2}}＝\dfrac{ {\text{∂}} }{ {\text{∂}} x}\left[ （ \lambda ＋ 2\mu ） \dfrac{ {\text{∂}} u}{ {\text{∂}} x} ＋ \lambda \dfrac{ {\text{∂}} v}{ {\text{∂}} {\textit{z}}}\right] ＋ \dfrac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\left[\mu \left(\dfrac{ {\text{∂}} v}{ {\text{∂}} x} ＋ \dfrac{ {\text{∂}} u}{ {\text{∂}} {\textit{z}}}\right)\right] ＋ {f}_{x} \text{，} \\ \rho \dfrac{{ {\text{∂}} }^{2}v}{ {\text{∂}} {t}^{2}}＝\dfrac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\left[\lambda \dfrac{ {\text{∂}} u}{ {\text{∂}} x} ＋ （ \lambda ＋ 2\mu ） \dfrac{ {\text{∂}} v}{ {\text{∂}} {\textit{z}}}\right] ＋ \dfrac{ {\text{∂}} }{ {\text{∂}} x}\left[\mu \left(\dfrac{ {\text{∂}} v}{ {\text{∂}} x} ＋ \dfrac{ {\text{∂}} u}{ {\text{∂}} {\textit{z}}}\right)\right] ＋ {f}_{{\textit{z}}} \text{，} \end{array}\right.\end{array} $$ (5) 式中λ和μ为各项同性介质的拉梅常数。在式（5）的右端项中，对于只包含关于x或z的偏导数项的FDM离散格式，其离散格式具有相似性。以$ \frac{ {\text{∂}} }{ {\text{∂}} x}\left[\left(\lambda ＋ 2\mu \right)\frac{ {\text{∂}} u}{ {\text{∂}} x}\right] $为例介绍，2N阶有限差分近似格式为：
$$ \begin{split} \frac{ {\text{∂}} }{ {\text{∂}} x}{\left[ （ \lambda ＋ 2\mu ） \frac{ {\text{∂}} u}{ {\text{∂}} x}\right]}_{{i}_{0}\text{，} {j}_{0}} {\text{≈}} & \frac{1}{\Delta x}\sum _{i=0}^{N-1}{b}_{i}\left[{ （ \lambda ＋ 2\mu ） }_{{i}_{0} ＋ i ＋ \frac{1}{2}\text{，} {j}_{0}}{\left(\frac{ {\text{∂}} u}{ {\text{∂}} x}\right)}_{{i}_{0} ＋ i ＋ \frac{1}{2}\text{，} {j}_{0}}-{ （ \lambda ＋ 2\mu ） }_{{i}_{0}-i-\frac{1}{2}\text{，} {j}_{0}}{\left(\frac{ {\text{∂}} u}{ {\text{∂}} x}\right)}_{{i}_{0}-i-\frac{1}{2}\text{，} {j}_{0}}\right] =\\ & \frac{1}{{\Delta x}^{2}}\sum _{i=0}^{N-1}{b}_{i}\left[{ （ \lambda ＋ 2\mu ） }_{{i}_{0} ＋ i ＋ \frac{1}{2}\text{，} {j}_{0}}\sum _{j=0}^{N-1}{b}_{j}\left({u}_{{i}_{0} ＋ i ＋ j ＋ 1\text{，} {j}_{0}}-{u}_{{i}_{0} ＋ i-j\text{，} {j}_{0}}\right) \right. -\\ & \left. { （ \lambda ＋ 2\mu ） }_{{i}_{0}-i-\frac{1}{2}\text{，} {j}_{0}}\sum _{j=0}^{N-1}{b}_{j}\left({u}_{{i}_{0}-i ＋ j\text{，} {j}_{0}}-{u}_{{i}_{0}-i ＋ j-1\text{，} {j}_{0}}\right)\right]\text{\text{，}} \end{split} $$ (6) 式中下标i0和j0分别为x和z方向上的有限差分网格点编号，bi为有限差分系数（Liu，2014），Δx为x方向的网格间距。应该说明的是，式（6）涉及半网格点上的拉梅系数值，可利用整数网格点上的值和调和平均公式（Ma et al，2018）计算。以$ \frac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\left[\lambda \frac{ {\text{∂}} v}{ {\text{∂}} x}\right] $为例介绍式（5）中的混合导数有限差分离散格式，$ \frac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\left[\lambda \frac{ {\text{∂}} v}{ {\text{∂}} x}\right] $的近似公式如下（Ma et al，2018）：
$$ \begin{split} {\left[\frac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\left(\lambda \frac{ {\text{∂}} v}{ {\text{∂}} x}\right)\right]}_{{i}_{0} \text{，} {j}_{0}}＝& \sum _{j＝1}^{N}{\omega }_{j}\frac{\left[{\lambda }_{{i}_{0} \text{，} {j}_{0} ＋ j}{\left(\frac{ {\text{∂}} v}{ {\text{∂}} x}\right)}_{{i}_{0} \text{，} {j}_{0} ＋ j}-{\lambda }_{{i}_{0} \text{，} {j}_{0}-j}{\left(\frac{ {\text{∂}} v}{ {\text{∂}} x}\right)}_{{i}_{0} \text{，} {j}_{0-j}}\right]}{\mathrm{\Delta }{\textit{z}}}＝\\& \sum _{j＝1}^{N}\frac{{\omega }_{j}}{\mathrm{\Delta }{\textit{z}}}\left[{\lambda }_{{i}_{0} \text{，} {j}_{0} ＋ j}\sum _{{i}_{0}＝1}^{N}{\omega }_{{i}_{0}}\frac{\left[{v}_{{i}_{0} ＋ i \text{，} {j}_{0} ＋ j}-{v}_{{i}_{0}-i \text{，} {j}_{0} ＋ j}\right]}{\mathrm{\Delta }x}-{\lambda }_{{i}_{0} \text{，} {j}_{0}-j}\sum _{{i}_{0}＝1}^{N}{\omega }_{{i}_{0}}\frac{\left[{v}_{{i}_{0} ＋ i \text{，} {j}_{0}-j}-{v}_{{i}_{0}-i \text{，} {j}_{0}-j}\right]}{\mathrm{\Delta }x}\right]＝\\&\frac{1}{\mathrm{\Delta }x\mathrm{\Delta }{\textit{z}}}\sum _{{j}_{0}＝1}^{N}\frac{{\omega }_{{j}_{0}}}{2{j}_{0}}\left[{\lambda }_{{i}_{0} \text{，} {j}_{0} ＋ j}\sum _{{i}_{0}＝1}^{N}\frac{{\omega }_{{i}_{0}}}{2{i}_{0}}\left({v}_{i ＋ {i}_{0} \text{，} j ＋ {j}_{0}}-{v}_{i-{i}_{0} \text{，} j ＋ {j}_{0}}\right)\right.-\\&\left. {\lambda }_{{i}_{0} \text{，} {j}_{0}-j}\sum _{{i}_{0}＝1}^{N}\frac{{\omega }_{{i}_{0}}}{2{i}_{0}}\left({v}_{i ＋ {i}_{0} \text{，} j-{j}_{0}}-{v}_{i-{i}_{0} \text{，} j-{j}_{0}}\right)\right]\text{\text{，}} \end{split} $$ (7) 式中Δz为z方向网格间距，ωi为有限差分系数。应该指出的是，虽然bi和ωi都为有限差分系数，但其值不同。
使用SEM模拟起伏地表附近地震波传播，利用FDM模拟远离地表处的地震波场，即为SEM-FDM混合方法，该方法的优势在于网格建模只需关注起伏地表地形，而无需考虑远离地表处的非规则结构，这种处理方式在保持SEM对起伏地表模型处理能力的同时，极大地简化了SEM网格建模（图1）。SEM-FDM混合方法的核心在于在SEM网格的底边界和FDM的上边界设置数据交换区，设置数据交换区是为了在SEM网格的底边界和FDM网格的上边界满足位移连续性边界条件。将计算区域Ω分解为Ω1和Ω2，即Ω=Ω1＋Ω2，Ω1为模型浅部区域，包含起伏地表，Ω2为模型深部区域，地震波运动方程的求解可表示为（假设震源位于Ω1）：
![]()
图 1 SEM-FDM混合方法的网格模型示意图Figure 1. Schematic of the grid model in the SEM-FDM hybrid method$$ \left\{\begin{array}{l}\ddot{\boldsymbol{U}}＝\left(\genfrac{}{}{0pt}{}{{\boldsymbol{U}}_{1}}{{\boldsymbol{U}}_{2}}\right)\text{，}\boldsymbol{x}\in {\varOmega } \text{，} \\ {\ddot{\boldsymbol{U}}}_{1}＝{K}^{\mathrm{S}\mathrm{E}\mathrm{M}}\left({\boldsymbol{U}}_{1} ＋ \boldsymbol{f}\right)\text{，}\boldsymbol{x}\in {{\varOmega }}_{1} \text{，} \\ {\ddot{\boldsymbol{U}}}_{2}＝{K}^{\mathrm{F}\mathrm{D}\mathrm{M}}\left({\boldsymbol{U}}_{2}\right)\text{，}\boldsymbol{x}\in {{\varOmega }}_{2} \text{，} \end{array}\right. $$ (8) 式中U为离散网格点上的位移向量，其下标1和2分别代表Ω1和Ω2研究区域，KSEM表示SEM的空间离散算子（Liu et al，2017），KFDM表示FDM的空间离散算子（Liu et al，2015a）。当有限差分的阶数为2N时，数据交换区的厚度为$N · \Delta {\textit{z}} $，即Ω1和Ω2重叠区域的厚度为$N · \Delta {\textit{z}} $。此时通过以下三步操作保证SEM网格的底边界和FDM网格的上边界满足位移连续边界条件。第一，在t时刻，由式（7）计算Ω1区域GLL节点上的加速度值，以及Ω2区域有限差分节点上的加速度值，并根据时间中心差分公式（刘少林等，2022），得到t＋Δt （其中Δt为时间离散步长）时刻GLL节点和有限差分节点上的波场值。第二，在t＋Δt时刻，数据交换区域第一至第N层有限差分网格点上的波场可通过邻近GLL点上的波场值插值近似（图2a），插值近似多项式为：
![]()
图 2 数据交换区有限差分（a）和GLL （b）网格点上网格点上的波场数据插值示意图（a） 数据交换区的位移插值；（b） 数据交换区底边界GLL点上的位移插值Figure 2. Schematic of data interpolations for the wavefields on the FD （a） and GLL （b） grid nodes in the data exchange domain$$ \left\{\genfrac{}{}{0pt}{}{{{\boldsymbol{U}}}_{\mathrm{F}\mathrm{D}}\left(x,{\textit{z}}\right) {\text{≈}} \sum _{i＝0}^{M}\sum _{j＝0}^{M}{L}_{ij}{{\boldsymbol{U}}}_{\mathrm{G}\mathrm{L}\mathrm{L}}\left({x}_{i},{{\textit{z}}}_{j}\right),}{{L}_{ij}\left(x,{\textit{z}}\right)＝\prod _{k＝0 \text{，} k\ne i}^{M}\dfrac{x-{x}_{k}}{{x}_{i}-{x}_{k}}\prod _{k＝0 \text{，} k\ne j}^{M}\dfrac{{\textit{z}}-{{\textit{z}}}_{k}}{{{\textit{z}}}_{j}-{{\textit{z}}}_{k}},}\right. $$ (9) 式中下标FD表示有限差分网格节点，下标GLL表示SEM的GLL网格节点，L为拉格朗日插值多项式，M为插值多项式阶数。实际计算表明，插值阶数为3时，可保证数据交换的精度较高；进一步提高插值阶数，对于数据交换精度提升贡献较小。第三，在t＋Δt时刻，数据交换区底边界GLL点的位移值可通过临近的有限差分节点上的波场值插值近似（图2b），近似公式如下：
$$ \boldsymbol{U}_{\mathrm{G}\mathrm{L}\mathrm{L}}\left(x,\text{z}\right)\text{≈}\sum_{i＝0}^{\mathrm{\mathit{M}}}\sum_{j＝0}^{\mathrm{\mathit{M}}}L_{ij}\boldsymbol{U}_{\mathrm{F}\mathrm{D}}\left(x_i,\text{z}_j\right). $$ (10) SEM-FDM混合方法模拟地震波场的流程可由图3描述。相比于SEM和FDM，SEM-FDM混合方法通过设置数据交换区的方式满足研究区内部位移连续性边界条件。应该指出的是，在t＋Δt时刻，在近似数据交换区第一至第N层有限差分网格点上的位移时，避免使用数据交换区底边界上GLL点上的位移值用于插值近似。此外，在数据交换区，除了底边界GLL点上位移值需要插值近似得到t＋Δt时刻值，其它GLL点上的位移值已经满足连续性边界条件而无需再插值近似。
![]()
图 3 SEM-FDM混合方法模拟地震波传播示意图Figure 3. Schematic of wavefield simulation using the SEM-FDM hybrid method数值算法精度主要由网格尺度和算法数值阶数来决定。研究表明网格尺寸（或一个波长内的采样点数）是影响数值算法的最主要因素（Liu et al，2015b，2017）。在使用SEM-FDM混合方法模拟地震波场时，为了避免因SEM和FDM精度不匹配而在数据交换区产生虚假散射，在网格设计时，SEM的GLL点平均间距应与FDM的网格间距相近。SEM和FDM的高阶格式能有效较弱数值频散（De Basabe，Sen，2007；Liu et al，2015a），但数值格式阶数过高会导致计算效率下降，SEM的常用阶数为4—5 （Komatitsch，Tromp，2002），FDM的常用差分阶数为4—6 （Tan，Huang，2014）。为了在数值精度与计算效率之间达到平衡，建议一个波长内采样点数为8—10，SEM-FDM混合方法中SEM的阶数选择5，FDM的阶数选择8—10。
2. 稳定性分析
根据Ma等（2014）和Liu等（2015b，2017）的研究，地震波运动方程的稳定性条件可表示为：
$$ \frac{\Delta t\mathrm{\mathit{v}}_{\mathrm{P}}}{h}\sqrt{1＋\frac{a}{r^2}}\text{≤}b\text{，} $$ (11) 式中vP为P波速度，γ＝vP/vS为P波与S波速度比，h为网格间距（或GLL点的平均间距），a和b为待定常数。De Basabe等（2007）和Liu等（2017）指出，当SEM的空间插值阶数为1时，SEM与标准网格2阶的FDM等价（Liu et al，2017），此时a和b的值都为1；当SEM的空间插值阶数大于1时，a为0，b值列于表1。在对式（5）中只包含x或z方向的导数项近似时，由于利用了半网格点处的空间一阶导数，以至于本文SEM-FDM混合方法中的FDM差分格式与标准网格FDM有一定差别，且稳定性条件与标准网格的FDM之间具有一定差别。根据平面波分析原理（Liu et al，2015b）可得到SEM-FDM混合方法中FDM的稳定性条件，a和b的取值见表2。
表 1 空间不同插值阶数SEM的稳定性参数a和b值Table 1. The stability parameters a and b for the SEMs with different interpolation orders空间插值阶数 a b 1 1 1 2 0 0.66 3 0 0.56 4 0 0.48 5 0 0.41 6 0 0.35 7 0 0.31 8 0 0.28 9 0 0.25 表 2 空间不同差分阶数FDM的稳定性参数a和b值Table 2. The stability parameters a and b for the different order FDMs空间插值阶数 a b 2 1.08 1.00 4 1.07 0.80 6 0.93 0.73 8 0.73 0.69 10 0.70 0.67 12 0.55 0.65 14 0.36 0.63 16 0.31 0.62 18 0.21 0.61 在获得了SEM和FDM各自的稳定性条件以后，可根据SEM和FDM稳定性的下限得到SEM-FDM混合方法的稳定性条件，例如SEM-FDM混合方法中SEM采用5阶插值和FDM采用8阶空间差分，其稳定性条件为：
$$ \begin{split} \\[-8pt] \frac{{\Delta t\mathrm{V}}_{\mathrm{P}}}{h} {\text{≤}} 0.41. \end{split} $$ (12) 3. 数值算例
本节将利用四个数值算例检验SEM-FDM混合方法的计算精度与效率。数值算例中SEM-FDM混合方法的SEM阶数选为5，FDM空间差分阶数选为8。为了减弱边界截断而引起的虚假反射，采用二阶位移形式的完美匹配层（perfectly matched layer，缩写为PML）吸收边界条件（刘有山等，2012）吸收虚假反射。采有较小的时间步长，以减弱时间离散误差对结果的干扰（衡量时间离散大小的标准是式（10）中的库朗数），3.1—3.3节中的时间离散间隔为0.1 ms，3.4节中时间离散间隔为1 ms。
3.1 Lamb问题
均匀介质模型如图4所示，P波速度为3.2 km/s，S波速度为1.8 km/s，密度ρ为2.6×103 kg/m3，一个单位集中力位于（
1000 m，100 m）处，震源时间函数是频率为15 Hz的雷克（Ricker）子波。两个接收器R1和R2分别位于（1500 m，800 m）和（1700 m，0 m）。模型0—1 km深度部分用SEM计算，正方形单元的边长为20 m；模型1—2 km深度部分用FDM计算，网格间隔为5 m。数据交换区位于z＝1000 m附近。![]()
图 4 用于检验SEM-FDM混合方法计算效率的均匀介质模型模型0—1 km深度部分用SEM计算；1—2 km部分用FDM计算。黑色五角星代表震源；两个黑色三角形对应两个接收器Figure 4. A homogeneous model used to verify the efficiency of the SEM-FDM hybrid methodThe wavefields in upper （0−1 km depths） and lower （1−2 km depths） parts of the model are calculated by the SEM and FDM，respectively. A star represents a source；two triangles denote two receivers R1 and R2上述模型的波场模拟结果如图5所示。为了对比需要，图5也展示了5阶SEM和8阶FDM的波场模拟结果。从t＝0.5 s的快照可看出SEM-FDM混合方法、SEM和FDM波场模拟结果大致相同（图5）。波场快照显示除了明显的直达波，还可看到地表的反射波和转换波；由于震源离地表较近，可在地表附近观测到明显的瑞雷面波。从图5a和5b可看到地震波场具有较好的连续性，在z＝1 000 m附近的数据交换区未见明显的虚假数值散射现象。
![]()
图 5 均匀介质模型中t＝0.5 s时刻的波场水平分量（左）和垂直分量（右）快照（a，b） SEM-FDM方法计算；（c，d） SEM方法计算；（e，f） FDM方法计算Figure 5. Snapshots of displacement wavefield of horizontal （left） and vertical （right） components in a homogeneous model at t＝0.5 s generated by SEM-FDM （a，b），SEM （c，d） and FDM （e，f），respectively为了说明SEM-FDM混合方法的波场模拟精度，图6展示了合成波形模拟结果，并将数值模拟结果与解析解（De Hoop，1959）对比。可以看出：在模型内部的接收器R1处SEM-FDM混合方法、SEM和FDM的模拟结果与解析解均较接近（图6a和6b），最大相对误差约为3%（图6c和6d）；在地表接收器R2处SEM-FDM混合方法和SEM模拟精度较高（图6e和6f），FDM模拟结果存在明显误差（图6e和6f），最大相对误差约15%，其误差明显大于SEM-FDM混合方法和SEM，这是缘于FDM对自由地表边界条件近似精度减低（苏波等，2019）。由于SEM-FDM混合方法能自动满足自由地表边界条件，因此能高精度模拟瑞雷面波（图6e和6f）。
![]()
图 6 由解析解得到的合成地震记录及其与三种数值方法得到的地震记录对比（a，b） R1处位移水平和垂直分量；（c，d） R1处合成波形水平和垂直分量的误差；（e，f） R2处位移水平和垂直分量；（g，h） R2处合成波形水平和垂直分量的误差. 黑线、蓝线、红线和青色为解析解、SEM-FDM、SEM和FDMFigure 6. Comparison of synthetic seismograms generated by analytical solution and three numerical methods（a，b） The horizontal and vertical components of displacements recorded at the R1；（c，d） the errors of the horizontal and vertical components of the displacements at the R1；（e，f） the horizontal and vertical components of displacements at the R2，respectively；（g，h） the errors of the horizontal and vertical components of the displacements at the R2表3对比了SEM-FDM混合方法、SEM和FDM三种方法的计算内存和时间。从表3可看出SEM-FDM混合方法的计算内存和时间介于SEM与FDM之间。SEM除了需要存储GLL点上的波场值以外，还需要存储GLL节点上雅克比（Jacobi）矩阵行列式的值（刘少林等，2021，2022），因此SEM的内存开销要大于FDM，但SEM和FDM的内存消耗仍处于同一数量级（表3；刘少林等，2021）。由于SEM-FDM混合方法为SEM和FDM两种方法的综合，因此其内存消耗介于SEM与FDM之间。当使用SEM方法计算空间任意方向的微分时，需要将空间坐标转化为三个参考坐标的加权（刘少林等，2021），而规则网格的FDM方法不存在坐标转化过程（Liu，2014），因此SEM的计算量略大于FDM，但SEM和FDM的计算量仍处于统一数量级。因此对于SEM-FDM混合方法而言，其计算量介于SEM与FDM之间。
表 3 三种地震波正演模拟方法运行内存和计算时间对比Table 3. Computer memory and computation time of three forward modeling methods正演模拟方法 运行内存/kb 计算时间/s SEM-FDM 68157 912.12 SEM 80740 961.31 FDM 28835 821.25 3.2 层状模型
为了说明SEM-FDM混合方法对模型内界面的处理能力，设计如图7所示模型，该模型具有分层结构，且内部含倾斜界面。在SEM-FDM混合方法中，深度500 m之上采用SEM计算，SEM的网格单元平均尺寸为20 m，深度500 m之下利用FDM计算，FDM网格间隔为5 m，数据交换区位于z＝500 m附近。模型上层介质P波速度为2.9 km/s，S波速度为1.7 km/s，密度为2.3×103 kg/m3；下层介质P波速度为3.2 km/s，S波速度为1.8 km/s，密度为2.6×103 kg/m3。一个爆炸源位于（1 000 m，100 m），其时间函数为15 Hz的雷克子波。两个接收器R1和R2分别位于（500 m，800 m）和（500 m，
1400 m）。![]()
图 7 用于地震波场模拟的层状模型模型内部存在一倾斜界面（黑色实线）. 五角星代表震源；黑色三角形为两个接收器R1和R2Figure 7. A layered model used for simulation of seismic wavefieldsAn inclined interface （solid black line） is in the model. The star located at （1000 m，100 m） represents a source；black triangles R1 and R2 denotes two receiversSEM-FDM混合方法地震波场模拟结果如图8所示。为了对比，也给出了5阶SEM波场模拟结果。需要指出的是，SEM网格考虑了模型内部界面，即网格剖分时沿着倾斜界面剖分。图8显示SEM-FDM混合方法与SEM方法所得的波场快照模拟结果具有很好的一致性，从波场快照可看到穿过倾斜界面的透射P波和pP波，以及倾斜界面的反射Pp波、Ps波和pPp波。接收器R1和R2处的合成地震记录如图9所示，可以看出SEM-FDM混合方法模拟结果与参考波形（6阶SEM计算得到）高度一致，相对误差在1%左右，这说明SEM-FDM混合方法能较好处理含速度间断面的模型。
![]()
图 8 层状模型中t＝0.5 s时刻的波场水平分量（左）和垂直分量（右）快照（a，b） SEM-FDM计算所得；（c，d） SEM计算所得.Figure 8. Snapshots of displacement wavefields of horizontal component （left） and horizontal component （right） in a layered model at t＝0.5 s generated by the SEM-FDM （a，b） and SEM （c，d）The left plots show the displacement wavefields of horizontal component；left plots display the displacement wavefields of vertical component. The upper and lower plots are generated by the SEM-FDM and SEM，respectively![]()
图 9 混合方法得到的合成地震记录与参考解对比（a，b） R1处位移水平和垂直分量；（c，d） R1处合成波形水平和垂直分量的误差；（e，f） R2处位移水平和垂直分量；（g，h） R2处合成波形水平和垂直分量的误差Figure 9. Comparison of synthetic seismograms generated by the SEM-FDM hybrid method and reference solution（a，b） The horizontal and vertical components of displacements recorded at the R1，respectively；（c，d） Errors of the horizontal and vertical components of the displacements at the R1；（e，f） the horizontal and vertical components of displacements at the R2；（g，h） the errors of the horizontal and vertical components of the displacements at the R23.3 起伏地表模型
为了检验SEM-FDM混合方法对起伏地表模型的处理能力，设计了如图10所示的起伏地表模型。模型分为两层，每层介质参数与3.2节中模型相同。一个爆炸震源位于（1 000 m，100 m）处，震源时间函数是频率为15 Hz的雷克子波。200个接收器均匀放置于地表，其平均水平间隔为16 m。
![]()
图 10 含起伏地表的层状模型中地震波场模拟地表曲线代表起伏地形; 模型内部实线表示倾斜界面. 五角星表示震源，其位置为（1 000 m，100 m）. 三角形表示接收器，其平均水平间隔为16 m，为了显示需要，接收器每隔四个显示一个Figure 10. Simulation of seismic wavefields in a layered model with topographyThe topography is marked by the black line. The star located at （1 000 m，100 m） represents a source. The triangles denote receivers，the average horizontal distance of which is about 16 m. For visualization purposes，only one-fifth stations are displayed单炮地震记录模拟结果如图11所示。图11a和11d由SEM-FDM混合方法计算得到；为了对比需要，图11b和11e展示了5阶SEM计算结果，可见SEM-FDM混合方法和SEM计算结果具有很好的一致性；图11c和11f展示了两种方法模拟结果的相对差，相对差在0.8%以内。通过对比可见，SEM-FDM混合方法保持了SEM对起伏地表的处理能力，能高精度模拟起伏地表模型中的地震波传播。
![]()
图 11 SEM-FDM混合方法计算的单炮地震记录与SEM计算结果对比图（a）和（b）分别为SEM-FDM和SEM方法计算的位移水平分量，图（c）为其相对差；图（d）和（e）分别为SEM-FDM和SEM方法计算的位移垂直分量，图（f）为其相对差Figure 11. A comparison of single shot profiles generated by the SEM-FDM hybrid method and SEMPlots （a） and （b） show displacement wavefields of horizontal component by the SEM-FDM and SEM，respectively，and plot （c） shows the relative differences between the wavefields shown in plots （a） and （b）. Plots （d） and （e） display displacement wavefields of vertical component by the SEM-FDM and SEM，respectively，and plot （f） displays the relative differences between the wavefields shown in plots （d） and （e）3.4 龙门山断裂带地壳模型
为了检验SEM-FDM混合方法在实际地质模型中的地震波模拟能力，选取龙门山断裂带地壳模型（图12）。图12显示龙门山断裂带东西两侧的地形相差较大，东侧的高程较低，西侧的高程较高。模型P波（图12b）和S波（图12c）速度结构来源于川滇地区速度模型2.0 （Liu等，2023）。选取二维垂直剖面（图12a中红色线表示其位置）用于正演合成测试。一个爆炸震源位于（110 km，13 km），震源时间函数是主频为1 Hz的雷克子波。两个接收器R1和R2分别位于（77.4 km，−3.4 km）和（136.4 km，−0.6 km）（其中z坐标的负值表示海平面之上）。SEM-FDM混合方法在深度z＝5 km之上利用SEM计算，网格单元的平均尺寸为800 m，在z＝5 km之下采用FDM计算，网格间隔为320 m，数据交换区位于z＝5 km附近。
![]()
图 12 龙门山断裂带地壳模型F1：龙门山断裂带（LMSF）；F2：岷江断裂（MJF）；F3：虎牙断裂（HYF）；F4：龙泉山断裂（LQSF）；F5：安宁河断裂（ANHF）；SGB: 松潘-甘孜块体; SCB: 四川盆地. （a） 地质构造图. 红色线为AB垂直剖面所在位置，用于二维数值测试；（b） P波速度结构；（c） S波速度结构Figure 12. The crustal model of the Longmenshan fault zoneF1：Longmenshan fault zone （LMSF）；F2：Minjiang fault （MJF）；F3：Huya fault （HYF）；F4：Longquanshan fault （LQSF）；F5：Anninghe fault （ANHF）. （a） The topography around the Longmenshan fault zone. The black lines denote faults. The red line AB represent a vertical profile，that is used for a 2D numerical test.（b） P-wave velocity model；（c） S-wave velocity modelSEM-FDM混合方法的地震波数值模拟结果如图13和14所示。为了对比需要，图13和14也展示了5阶SEM数值模拟结果。可见：SEM-FDM混合方法得到的t＝13 s时刻波场快照显示了清晰的P，pP和pS震相，波场模拟结果（图13a和13c）与SEM模拟结果一致（图13b和图13d）。R1和R2处的合成波形模拟结果（图14）表明SEM-FDM混合方法与SEM得到的合成波形具有很好的一致性。以上模拟结果说明SEM-FDM能模拟复杂地壳模型中的地震波场。
![]()
图 13 地壳模型中t＝13 s时刻的波场快照. 左图为位移水平分量；右图为位移垂直分量（a）和（b）为SEM-FDM混合方法计算得到；（c）和（d）为SEM计算得到Figure 13. Snapshots of wavefields in the crustal model at t＝13 s. The left and right plots show horizontal and vertical components of displacement wavefieldsPlots （a） and （b） are generated by the SEM-FDM hybrid method；plots （c） and （d） are computed by the SEM![]()
图 14 地壳模型中合成地震记录左图为位移水平分量；右图为位移垂直分量。图（a）和（b）为R1处合成地震记录；图（c）和（d）为R2处合成地震记录. 黑线为参考地震图; 蓝线为SEM-FDM混合方法计算得到Figure 14. Synthetic seismograms in a crustal modelLeft and right plots show synthetic seismograms of horizontal and vertical components，respectively；Plots （a） and （b） present seismograms recorded by the R1；plots （a） and （b） display seismograms recorded by the R2. The black lines are the reference seismograms. The blue lines are the synthetic seismograms computed by the SEM-FDM hybrid method4. 讨论与结论
本文提出了模拟地震波场的SEM-FDM混合方法，该方法在接近地表区域利用SEM模拟地震波传播，在远离地表区域采用FDM计算地震波场，通过设置数据交换区将SEM和FDM耦合，使得模型内部的位移连续。发展SEM-FDM混合方法是为了简化网格建模，即在网格建模时只用考虑起伏地表附近模型的非规则性，在远离地表区域利用规则有限差分网格刻画模型内部界面。SEM-FDM混合方法综合了SEM和FDM的各自优势，可简单、便捷地模拟复杂模型中的地震波场。通过四组数值算例证实了SEM-FDM混合方法模拟地震波场的有效性。本文设计的地下介质模型相对简单，主要目的是为了简化参考、对比对象SEM的网格建模，便于SEM-FDM混合方法计算结果与SEM对比。SEM-FDM混合方法可适用于任意模型，并不受模型的限制。
在利用SEM-FDM混合方法模拟地震波场过程中，利用商业软件Hypermesh （https://www.altair.com）建立地表附近SEM网格模型，为了进一步提高网格建模效率，可构造高阶多项式拟合地表地形，在此基础之上自动生成SEM网格，作者将另文介绍SEM网格自动生成技术。基于网格自动生成技术，SEM-FDM混合方法将成为一种高效、便捷的正演模拟工具，可模拟任意复杂介质中的地震波场。
在波形伴随成像过程中，随着模型迭代更新，反演得到的速度模型与初始模型之间可能存在较大差别。目前，基于SEM的波形伴随成像在模型迭代过程中难以根据反演模型重新剖分网格。由于SEM网格未能有效地刻画反演后的速度结构异常体，这必然造成合成波形精度较低，进而影响成像结果。本文的SEM-FDM混合方法利用较密的有限差分网格刻画速度结构的高波数成分，因此无须在反演迭代过程中再重新划分网格。作者将另文介绍基于SEM-FDM混合方法的波形伴随成像方法与软件。
-
![]()
图 1 SEM-FDM混合方法的网格模型示意图
Figure 1. Schematic of the grid model in the SEM-FDM hybrid method
![]()
图 2 数据交换区有限差分（a）和GLL （b）网格点上网格点上的波场数据插值示意图
（a） 数据交换区的位移插值；（b） 数据交换区底边界GLL点上的位移插值
Figure 2. Schematic of data interpolations for the wavefields on the FD （a） and GLL （b） grid nodes in the data exchange domain
![]()
图 3 SEM-FDM混合方法模拟地震波传播示意图
Figure 3. Schematic of wavefield simulation using the SEM-FDM hybrid method
![]()
图 4 用于检验SEM-FDM混合方法计算效率的均匀介质模型
模型0—1 km深度部分用SEM计算；1—2 km部分用FDM计算。黑色五角星代表震源；两个黑色三角形对应两个接收器
Figure 4. A homogeneous model used to verify the efficiency of the SEM-FDM hybrid method
The wavefields in upper （0−1 km depths） and lower （1−2 km depths） parts of the model are calculated by the SEM and FDM，respectively. A star represents a source；two triangles denote two receivers R1 and R2
![]()
图 5 均匀介质模型中t＝0.5 s时刻的波场水平分量（左）和垂直分量（右）快照
（a，b） SEM-FDM方法计算；（c，d） SEM方法计算；（e，f） FDM方法计算
Figure 5. Snapshots of displacement wavefield of horizontal （left） and vertical （right） components in a homogeneous model at t＝0.5 s generated by SEM-FDM （a，b），SEM （c，d） and FDM （e，f），respectively
![]()
图 6 由解析解得到的合成地震记录及其与三种数值方法得到的地震记录对比
（a，b） R1处位移水平和垂直分量；（c，d） R1处合成波形水平和垂直分量的误差；（e，f） R2处位移水平和垂直分量；（g，h） R2处合成波形水平和垂直分量的误差. 黑线、蓝线、红线和青色为解析解、SEM-FDM、SEM和FDM
Figure 6. Comparison of synthetic seismograms generated by analytical solution and three numerical methods
（a，b） The horizontal and vertical components of displacements recorded at the R1；（c，d） the errors of the horizontal and vertical components of the displacements at the R1；（e，f） the horizontal and vertical components of displacements at the R2，respectively；（g，h） the errors of the horizontal and vertical components of the displacements at the R2
![]()
图 7 用于地震波场模拟的层状模型
模型内部存在一倾斜界面（黑色实线）. 五角星代表震源；黑色三角形为两个接收器R1和R2
Figure 7. A layered model used for simulation of seismic wavefields
An inclined interface （solid black line） is in the model. The star located at （1000 m，100 m） represents a source；black triangles R1 and R2 denotes two receivers
![]()
图 8 层状模型中t＝0.5 s时刻的波场水平分量（左）和垂直分量（右）快照
（a，b） SEM-FDM计算所得；（c，d） SEM计算所得.
Figure 8. Snapshots of displacement wavefields of horizontal component （left） and horizontal component （right） in a layered model at t＝0.5 s generated by the SEM-FDM （a，b） and SEM （c，d）
The left plots show the displacement wavefields of horizontal component；left plots display the displacement wavefields of vertical component. The upper and lower plots are generated by the SEM-FDM and SEM，respectively
![]()
图 9 混合方法得到的合成地震记录与参考解对比
（a，b） R1处位移水平和垂直分量；（c，d） R1处合成波形水平和垂直分量的误差；（e，f） R2处位移水平和垂直分量；（g，h） R2处合成波形水平和垂直分量的误差
Figure 9. Comparison of synthetic seismograms generated by the SEM-FDM hybrid method and reference solution
（a，b） The horizontal and vertical components of displacements recorded at the R1，respectively；（c，d） Errors of the horizontal and vertical components of the displacements at the R1；（e，f） the horizontal and vertical components of displacements at the R2；（g，h） the errors of the horizontal and vertical components of the displacements at the R2
![]()
图 10 含起伏地表的层状模型中地震波场模拟
地表曲线代表起伏地形; 模型内部实线表示倾斜界面. 五角星表示震源，其位置为（1 000 m，100 m）. 三角形表示接收器，其平均水平间隔为16 m，为了显示需要，接收器每隔四个显示一个
Figure 10. Simulation of seismic wavefields in a layered model with topography
The topography is marked by the black line. The star located at （1 000 m，100 m） represents a source. The triangles denote receivers，the average horizontal distance of which is about 16 m. For visualization purposes，only one-fifth stations are displayed
![]()
图 11 SEM-FDM混合方法计算的单炮地震记录与SEM计算结果对比
图（a）和（b）分别为SEM-FDM和SEM方法计算的位移水平分量，图（c）为其相对差；图（d）和（e）分别为SEM-FDM和SEM方法计算的位移垂直分量，图（f）为其相对差
Figure 11. A comparison of single shot profiles generated by the SEM-FDM hybrid method and SEM
Plots （a） and （b） show displacement wavefields of horizontal component by the SEM-FDM and SEM，respectively，and plot （c） shows the relative differences between the wavefields shown in plots （a） and （b）. Plots （d） and （e） display displacement wavefields of vertical component by the SEM-FDM and SEM，respectively，and plot （f） displays the relative differences between the wavefields shown in plots （d） and （e）
![]()
图 12 龙门山断裂带地壳模型
F1：龙门山断裂带（LMSF）；F2：岷江断裂（MJF）；F3：虎牙断裂（HYF）；F4：龙泉山断裂（LQSF）；F5：安宁河断裂（ANHF）；SGB: 松潘-甘孜块体; SCB: 四川盆地. （a） 地质构造图. 红色线为AB垂直剖面所在位置，用于二维数值测试；（b） P波速度结构；（c） S波速度结构
Figure 12. The crustal model of the Longmenshan fault zone
F1：Longmenshan fault zone （LMSF）；F2：Minjiang fault （MJF）；F3：Huya fault （HYF）；F4：Longquanshan fault （LQSF）；F5：Anninghe fault （ANHF）. （a） The topography around the Longmenshan fault zone. The black lines denote faults. The red line AB represent a vertical profile，that is used for a 2D numerical test.（b） P-wave velocity model；（c） S-wave velocity model
![]()
图 13 地壳模型中t＝13 s时刻的波场快照. 左图为位移水平分量；右图为位移垂直分量
（a）和（b）为SEM-FDM混合方法计算得到；（c）和（d）为SEM计算得到
Figure 13. Snapshots of wavefields in the crustal model at t＝13 s. The left and right plots show horizontal and vertical components of displacement wavefields
Plots （a） and （b） are generated by the SEM-FDM hybrid method；plots （c） and （d） are computed by the SEM
![]()
图 14 地壳模型中合成地震记录
左图为位移水平分量；右图为位移垂直分量。图（a）和（b）为R1处合成地震记录；图（c）和（d）为R2处合成地震记录. 黑线为参考地震图; 蓝线为SEM-FDM混合方法计算得到
Figure 14. Synthetic seismograms in a crustal model
Left and right plots show synthetic seismograms of horizontal and vertical components，respectively；Plots （a） and （b） present seismograms recorded by the R1；plots （a） and （b） display seismograms recorded by the R2. The black lines are the reference seismograms. The blue lines are the synthetic seismograms computed by the SEM-FDM hybrid method
表 1 空间不同插值阶数SEM的稳定性参数a和b值
Table 1 The stability parameters a and b for the SEMs with different interpolation orders
空间插值阶数 a b 1 1 1 2 0 0.66 3 0 0.56 4 0 0.48 5 0 0.41 6 0 0.35 7 0 0.31 8 0 0.28 9 0 0.25 
下载: 导出CSV
表 2 空间不同差分阶数FDM的稳定性参数a和b值
Table 2 The stability parameters a and b for the different order FDMs
空间插值阶数 a b 2 1.08 1.00 4 1.07 0.80 6 0.93 0.73 8 0.73 0.69 10 0.70 0.67 12 0.55 0.65 14 0.36 0.63 16 0.31 0.62 18 0.21 0.61
下载: 导出CSV
表 3 三种地震波正演模拟方法运行内存和计算时间对比
Table 3 Computer memory and computation time of three forward modeling methods
正演模拟方法 运行内存/kb 计算时间/s SEM-FDM 68157 912.12 SEM 80740 961.31 FDM 28835 821.25
下载: 导出CSV
-
刘少林,杨顶辉,徐锡伟,李小凡,申文豪,刘有山. 2021. 模拟地震波传播的三维逐元并行谱元法[J]. 地球物理学报,64(3):993–1005. Liu S L,Yang D H,Xu X W,Li X F,Shen W H,Liu Y S. 2021. Three-dimensional element-by-element parallel spectral-element method for seismic wave modeling[J]. Chinese Journal of Geophysics,64(3):993–1005 (in Chinese).
刘少林,杨顶辉,孟雪莉,汪文帅,徐锡伟,李小凡. 2022. 模拟地震波传播的优化质量矩阵Legendre谱元法[J]. 地球物理学报,65(12):4802–4815. Liu S L,Yang D H,Meng X L,Wang W S,Xu X W,Li X F. 2022. A Legendre spectral element method with optimal mass matrix for seismic wave modeling[J]. Chinese Journal of Geophysics,65(12):4802–4815 (in Chinese).
刘有山,刘少林,张美根,马德堂. 2012. 一种改进的二阶弹性波动方程的最佳匹配层吸收边界条件[J]. 地球物理学进展,27(5):2113–2122. Liu Y S,Liu S L,Zhang M G,Ma D T. 2012. An improved perfectly matched layer absorbing boundary condition for second order elastic wave equation[J]. Progress in Geophysics,27(5):2113–2122 (in Chinese).
桑莹泉,刘有山,徐涛,白志明,解桐桐. 2021. 远震波场正演模拟方法及应用[J]. 地球与行星物理论评,52(6):569–586. Sang Y Q,Liu Y S,Xu T,Bai Z M,Xie T T. 2021. Forward modeling method and application of teleseismic wavefield[J]. Reviews of Geophysics and Planetary Physics,52(6):569–586(in Chinese).
苏波,李怀良,刘少林,杨顶辉. 2019. 修正辛格式有限元法的地震波场模拟[J]. 地球物理学报,62(4):1440–1452. Su B,Li H L,Liu S L,Yang D H. 2019. Modified symplectic scheme with finite element method for seismic wavefield modeling[J]. Chinese Journal of Geophysics,62(4):1440–1452 (in Chinese).
祝贺君,刘沁雅,杨继东. 2023. 地震学全波形反演进展[J]. 地球与行星物理论评(中英文),53(3):287–317. Zhu H J,Liu Q Y,Yang J D. 2023. Recent progress on full waveform inversion[J]. Reviews of Geophysics and Planetary Physics,54(3):287–317 (in Chinese).
Galis M,Moczo P,Kristek J. 2008. A 3-D hybrid finite-difference-finite-element viscoelastic modelling of seismic wave motion[J]. Geophys J Int,175(1):153–184. doi: 10.1111/j.1365-246X.2008.03866.x
Chen M,Niu F,Liu Q,Tromp J,Zheng X. 2015. Multiparameter adjoint tomography of the crust and upper mantle beneath East Asia:1. Model construction and comparisons[J]. J Geophys Res-Sol Ea,120:1762–1786. doi: 10.1002/2014JB011638
De Basabe J D,Sen M K,2007. Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equation[J]. Geophysics,72(6):T81–T95.
De Hoop A T,1959. The surface line source problem[J]. Appl Sci Res,8(4):349–356.
Komatitsch D,Tromp J. 2002. Spectral-element simulations of global seismic wave propagation-Ⅰ. Validation[J]. Geophys J Int,149(2):390–412. doi: 10.1046/j.1365-246X.2002.01653.x
Komatitsch D,Vilotte J P. 1998. The spectral element method:an efficient tool to simulate the seismic response of 2D and 3D geological structures[J]. B Seismol Soc Am,88(2):368–392. doi: 10.1785/BSSA0880020368
Lan H Q,Zhang Z J. 2011. Comparative study of the free-surface boundary condition in two-dimensional finite-difference elastic wave field simulation[J]. J Geophys Eng. 8:275–286.
Lei W,Ruan Y Y,Bozdağ E,Peter D,Lefebvre M,Komatitsch D,Tromp J,Hill J,Podhorszki N,Pugmire D. 2020. Global adjoint tomography-model GLAD-M25[J]. Geophys J Int,233(1):1–21.
Liu Y,2014. Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling[J]. Geophys J Int,197(2):1033–1047.
Liu S L,Li X F,Wang W S,Xu L,Li B F. 2015a. A modified symplectic scheme for seismic wave modeling[J]. J Appl Geophys,116:110–120. doi: 10.1016/j.jappgeo.2015.03.007
Liu S L,Li X F,Wang W S,Zhu T. 2015b. Source wavefield reconstruction using a linear combination of the boundary wavefield in reverse time migration[J]. Geophysics,80(6):S203–S212. doi: 10.1190/geo2015-0109.1
Liu S L,Yang D H,Ma J. 2017. A modified symplectic PRK scheme for seismic wave modeling[J]. Comput Geosci-UK,99:28–36. doi: 10.1016/j.cageo.2016.11.001
Liu Y,Yu Z,Zhang Z Q,Yao H J,Wang W T,Zhang H J,Fang H J,Fang L H. 2023. The high-resolution community velocity model V2.0 of southwest China,constructed by joint body and surface wave tomography of data recorded at temporary dense arrays[J]. Sci China Earth Sci,66(10):2368–2385. doi: 10.1007/s11430-022-1161-7
Ma J,Yang D H,Tong P,Ma X. 2018. TSOS and TSOS-FK hybrid methods for modelling the propagation of seismic waves[J]. Geophys J Int,214(3):1665–1682. doi: 10.1093/gji/ggy215
Ma X,Yang D H,Song G J,Wang M X. 2014. A low-dispersion symplecitic partitioned Rung–Kutta method for solving seismic-wave equations:I. Scheme and theoretical analysis[J]. Bull Seismol Soc Am,104(5):2206–2225. doi: 10.1785/0120120210
Marfurt K J. 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[J]. Geophysics,49(5):533–549. doi: 10.1190/1.1441689
Moczo P,Robertsson J,Eisner L. 2007. The finite-difference time-domain method for modeling of seismic wave propagation[J]. Adv Geophys,48:421–516.
Rao Y,Wang Y. 2013. Seismic waveform simulation with pseudo-orthogonal grids for irregular topographic models[J]. Geophys J Int,194(3):1778–1788. doi: 10.1093/gji/ggt190
Shen W H,Yang D H,Xu X W,Yang S X,Liu S L. 2022. 3D simulation of ground motion for the 2015 MW7.8 Gorkha earthquake,Nepal,based on the spectral element method[J]. Nat Hazards,112:2853–2871. doi: 10.1007/s11069-022-05291-1
Tan P,Huang L J. 2014. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation[J]. Geophys J Int,197:1250–1267. doi: 10.1093/gji/ggu077
Tong P,Zhao D,Yang D,Chen J,Liu Q. 2014. Wave-equation-based travel-time seismic tomography-Part 1:Method[J]. Solid Earth,5(2):1151–1168. doi: 10.5194/se-5-1151-2014
Tromp J,Tape C,Liu Q Y. 2005. Seismic tomography,adjoint methods,time reversal and banana-doughnut kernels[J]. Geophys J Int,160(1):195–216.
Zhang J H,Yao Z X. 2013. Optimized finite-difference operator for broadband seismic wave modeling[J]. Geophysics,78(1):A13–A18. doi: 10.1190/geo2012-0277.1
Zhang W Q,Zhang Z G,Fu H H,Li Z B,Chen X F. 2019. Importance of Spatial Resolution in Ground Motion Simulations With 3-D Basins:An Example Using the Tangshan Earthquake[J]. Geophys Res Lett,46(21):11915–11924. doi: 10.1029/2019GL084815
Zhou H,Liu Y,Wang J. 2021. Elastic wave modeling with high-order temporal and spatial accuracies by a selectively modified and linearly optimized staggered-grid finite-difference scheme[J]. IEEE T Geosci Remote,60(60):1–22.
计量
- 文章访问数: 196
- HTML全文浏览量: 25
- PDF下载量: 67
