地球物理学进展

基于泰勒展开和采样逼近法的各向异性介质弹性波模拟时空高阶旋转交错网格法

PDF下载

李皖来 ¹ , 赵前进 ¹^,^* , 韦慧 ² , 王弈轩 ¹

地球物理学进展 | 应用地球物理学（油气、金属矿、海洋、环境、工程、仪器等） 2026,41(1): 289-310

收起

地球物理学进展 |应用地球物理学（油气、金属矿、海洋、环境、工程、仪器等） 2026 , 41 (1) : 289 -310

基于泰勒展开和采样逼近法的各向异性介质弹性波模拟时空高阶旋转交错网格法

全屏

李皖来¹, 赵前进¹^,^*, 韦慧², 王弈轩¹

作者信息

¹安徽理工大学数学与大数据学院, 淮南 232000

²桂林电子科技大学数学与计算科学学院, 桂林 541004

通讯作者:

赵前进, 男, 1967年生, 教授, 博士, 研究方向为有理插值与逼近、数字图像处理.E-mail: qjzhao@aust.edu.cn

作者简介:

李皖来, 女, 2000年生, 硕士研究生, 主要从事偏微分方程研究.E-mail: 3037104153@qq.com

Temporal high-order rotated-staggered-grid finite difference method of anisotropic media elastic wave equations

WanLai LI¹, QianJin ZHAO¹^,^*, Hui WEI², YiXuan WANG¹

Affiliations

¹School of Mathematics and Big Data, Anhui University of Science and Technology, Huainan 232000, China

²School of Mathematics and Computing Science, Guilin University of Electronic Technology, Guilin 541004, China

出版时间: 2026-02-20 doi: 10.6038/pg2026JJ0228

文章导航

摘要

收起

旋转交错网格有限差分法对各向异性介质弹性波波场进行数值模拟时, 通常存在数值频散现象.为此, 给出了基于泰勒展开与采样逼近法耦合的时域高阶差分系数计算模型, 构建了改进的时空域高阶旋转交错网格格式.为了验证改进算法的性能, 将改进算法与一般交错网格以及其他旋转交错网格进行比较, 分析各算法的频散程度以及稳定性.通过数值实验结果发现, 基于泰勒展开与采样逼近法的差分算子系数优化策略改进的算法在各向同性介质、VTI/TTI各向异性介质及修正BP模型中均能显著抑制弹性波场的频散.

关键词

弹性波方程 / 有限差分法 / 旋转交错网格 / 频散分析

Abstract

收起

When numerically simulating the elastic wave field in anisotropic media using the rotated staggered grid finite difference method, numerical dispersion phenomena usually occur. To solve this problem, a time-domain high-order difference coefficient calculation model based on the coupling of Taylor expansion and sampling approximation method is proposed, and an improved space-time domain high-order rotating staggered grid scheme is constructed. To verify the performance of the modified algorithm, the improved algorithm is compared with the general staggered grid and other rotated staggered grids, and the dispersion degree and stability of each algorithm are analyzed. Through the results of numerical experiments, it is found that the algorithm improved by the optimization strategy of the difference operator coefficients based on the Taylor expansion and the sampling approximation method can significantly suppress the dispersion of the elastic wave field in isotropic media, VTI/TTI anisotropic media, and the modified BP model.

Key words

Elastic wave equation / Finite-difference method / Rotated-staggered-grid / Dispersion analysis

引用本文

李皖来, 赵前进, 韦慧, 王弈轩. 基于泰勒展开和采样逼近法的各向异性介质弹性波模拟时空高阶旋转交错网格法. 地球物理学进展, 2026 , 41 (1) : 289 -310 . DOI: 10.6038/pg2026JJ0228

WanLai LI, QianJin ZHAO, Hui WEI, YiXuan WANG. Temporal high-order rotated-staggered-grid finite difference method of anisotropic media elastic wave equations[J]. Progress in Geophysics, 2026 , 41 (1) : 289 -310 . DOI: 10.6038/pg2026JJ0228

正文

收起

0　引言

收起

各向异性波动方程是地震勘探领域的一个重要研究内容，为研究复杂介质中地震波的传播机理提供了一种有效手段.目前，求解各向异性波动方程数值解中常用的方法，如有限差分法(Doostaki et al., 2023; Yan and Xu, 2024)、伪谱法(Lin et al., 2024)、有限元法(Zhao et al., 2023)、边界元法(Li et al., 2023)、无网格法(钟霖等, 2021)等.有限差分法的计算复杂度较低、计算效率高且便于实现，因而被广泛应用于地震波方程的数值模拟(Almquist and Dunham, 2021; Behera, 2022; 王啟等, 2025; 芦永明等, 2025; Zhou et al., 2023).然而，时空网格离散的疏密性往往会引起数值频散现象(Fei and Larner, 1995).因此，为了解决频散问题，研究者通过改进差分离散格式逐步优化算法性能.交错网格有限差分法(Staggered-Grid Finite-Difference，SFD)通过将应力与速度分量分离至不同空间网格节点，提升变量采样密度降低频散误差，显著提高了各向同性介质中的模拟精度而被广泛应用，然而，该方法在处理各向异性介质时，因应力-速度耦合关系与网格对齐方向失配，容易出现计算不稳定的问题.为克服此问题，旋转交错网格(Rotated-Staggered-Grid Finite-Difference, RSFD)通过旋转网格坐标系使所有物理量仅在两个不同的位置定义，在弹性模量变化较大时可以避免空间插值操作，从而在各向异性介质的模拟中获得更稳定和令人信服的解决方案(Chen et al., 2006).

Gold与Shapiro(Gold et al., 1997)于1997年提出旋转交错网格，通过动态调整局部网格坐标系与介质各向异性主轴方向的一致性，从而抑制由各向异性参数空间异质性引起的数值不稳定现象.Saenger等(Saenger et al., 2000)进一步系统化完善旋转交错网格理论框架，通过旋转法向坐标系实现应力与位移分量离散，并验证了其处理各向异性介质弹性波方程的稳定性与频散特性优势.旋转交错网格方法因强鲁棒性成为各向异性弹性波模拟的研究热点(魏皓等, 2024; 徐云贵等, 2023; Zhang and Schmitt, 2025)，但其离散格式仍存在显著数值频散效应，直接影响波场保真度(Yang et al., 2015).

数值频散现象是任何有限差分方法都不可避免的，传统基于泰勒级数展开的高阶交错网格格式在中、低频范围内可以部分抑制数值频散，但计算的差分系数不能满足计算精度的要求(Yang et al., 2015).针对这一关键问题，许多学者研究了不同空间差分系数的优化策略，致力于提升算法精度(王超群等, 2025; Liu and Sen, 2011; Yang et al., 2017; Zhang and Yao, 2013).Liu和Sen(2011)基于时空域色散关系与平面波理论，构建了一维标量波动方程的交错网格求解方案；Zhang和Yao(2013)采用模拟退火算法在预设误差容限内优化波数覆盖范围，实现了差分算子常数系数的最优配置；Yang等(2017)通过极小化极大算法改进泰勒展开法，提出了具有优化特性的隐式交错网格格式.上述研究虽提升了标准交错网格的数值精度，但针对旋转交错网格的频散抑制缺乏系统性解决方案.在旋转交错网格优化方面，Liu(2014)基于全局优化方法构建了RSFD最佳空间差分系数；Yang等(2015)针对横向各向同性(TTI)介质，分别基于采样逼近法(Sampling Approximation, SA)与最小二乘法，提出了两种高阶空间差分优化方案.值得注意的是，Liu和Yang等提出的方法显著地提高了各向异性介质弹性波方程空间偏导数的求解精度，但其时间偏导数项只具有二阶精度.为了进一步提高旋转交错网格时间算子的差分精度，Gao和Huang(2017)通过三阶泰勒展开将时间导数转换为空间导数，并基于Khan和Ohba (2003)的高阶导数差分近似方法计算差分算子系数，使得各向异性介质弹性波方程具有四阶时间精度.然而，该方案很难直接拓展使得时间精度具有更高的阶数.最新研究进展中，针对TTI介质弹性波方程，Xie和He(2024)结合Lax-Wendroff校正法与网格加密技术，通过最小二乘优化实现了时空精度的协同提升，但该方法在强各向异性介质中的频散抑制作用仍需进一步优化提升.

本文提出一种融合泰勒级数展开法(Taylor Series Expansion, TE)与采样逼近法(SA)的时空域高阶RSFD优化方法，旨在解决各向异性二维三分量弹性波方程数值模拟中的时空耦合频散抑制难题.为了方便起见，下面简称该方法为SA+TE-RSFD.首先，系统阐述SFD、RSFD及Xie-RSFD方法的基本原理，推导SA+TE-RSFD的时域高阶差分系数计算公式.其次，对各算法进行频散与稳定性分析.最后，针对各向同性介质、VTI/TTI各向异性介质及改进BP模型进行数值模拟，通过对比一般交错网格以及其他旋转交错网格算法的波场快照，讨论各算法对数值频散的抑制性.

1　算法设计

收起

以二维各向异性介质的三分量一阶速度-应力弹性波方程为例，研究旋转交错网格的基本格式，并结合采样逼近和泰勒展开公式对差分系数优化.控制方程如下：

$\left\{\begin{array}{l}\rho \frac{\partial v_x}{\partial t}=\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \sigma_{x z}}{\partial z}, \rho \frac{\partial v_y}{\partial t}=\frac{\partial \sigma_{x y}}{\partial x}+\frac{\partial \sigma_{y z}}{\partial z}, \rho \frac{\partial v_z}{\partial t}=\frac{\partial \sigma_{x z}}{\partial x}+\frac{\partial \sigma_{z z}}{\partial z} \\\frac{\partial \sigma_{x x}}{\partial t}=D_{11} \frac{\partial v_x}{\partial x}+D_{13} \frac{\partial v_z}{\partial z}+D_{14} \frac{\partial v_y}{\partial z}+D_{15}\left(\frac{\partial v_z}{\partial x}+\frac{\partial v_x}{\partial z}\right)+D_{16} \frac{\partial v_y}{\partial x} \\\frac{\partial \sigma_{z z}}{\partial t}=D_{13} \frac{\partial v_x}{\partial x}+D_{33} \frac{\partial v_z}{\partial z}+D_{34} \frac{\partial v_y}{\partial z}+D_{35}\left(\frac{\partial v_z}{\partial x}+\frac{\partial v_x}{\partial z}\right)+D_{36} \frac{\partial v_y}{\partial x} \\\frac{\partial \sigma_{y z}}{\partial t}=D_{14} \frac{\partial v_x}{\partial x}+D_{34} \frac{\partial v_z}{\partial z}+D_{44} \frac{\partial v_y}{\partial z}+D_{45}\left(\frac{\partial v_z}{\partial x}+\frac{\partial v_x}{\partial z}\right)+D_{46} \frac{\partial v_y}{\partial x} \\\frac{\partial \sigma_{x z}}{\partial t}=D_{15} \frac{\partial v_x}{\partial x}+D_{35} \frac{\partial v_z}{\partial z}+D_{45} \frac{\partial v_y}{\partial z}+D_{55}\left(\frac{\partial v_z}{\partial x}+\frac{\partial v_x}{\partial z}\right)+D_{56} \frac{\partial v_y}{\partial x} \\\frac{\partial \sigma_{x y}}{\partial t}=D_{16} \frac{\partial v_x}{\partial x}+D_{36} \frac{\partial v_z}{\partial z}+D_{46} \frac{\partial v_y}{\partial z}+D_{56}\left(\frac{\partial v_z}{\partial x}+\frac{\partial v_x}{\partial z}\right)+D_{66} \frac{\partial v_y}{\partial x}\end{array}, \right.$

(1)

式中：σ=σ(x, z)表示应力；v_x=v_x(x, z)、v_y=v_y(x, z)、v_z=v_z(x, z)分别为质点振动速度场在三个正交方向x、y和z上的分量；ρ=ρ(x, z)为介质密度；D₁₁~D₆₆为刚度系数，由Thomsen参数或介质弹性参数组成.

下面介绍旋转交错网格、时空高阶旋转交错网格的基本理论及差分系数计算方法，进而给出一种基于系数优化的高阶时空域旋转交错网格有限差分法.

1.1　旋转交错网格法

针对各向异性介质中弹性波动方程，采用传统交错网格有限差分法(Staggered-Grid Finite-Difference, SFD)离散时，部分空间导数需通过空间插值逼近，这将产生误差，从而导致数值解的精度降低(Chen et al., 2006).2000年，Saenger等(Saenger et al., 2000)通过旋转法向坐标系，使所有应力分量和位移分量都在同一位置，从而避免空间插值提高算法精度，即旋转交错网格(Rotated-Staggered-Grid Finite-Difference, RSFD).该算法中，旋转后的坐标系x^R和z^R方向与法向坐标系x和z方向上的空间导数之间的关系为：

$\frac{\partial}{\partial x}=\frac{\Delta r}{2 \Delta x}\left(\frac{\partial}{\partial x^{\mathrm{R}}}-\frac{\partial}{\partial z^{\mathrm{R}}}\right), \frac{\partial}{\partial z}=\frac{\Delta r}{2 \Delta z}\left(\frac{\partial}{\partial x^{\mathrm{R}}}+\frac{\partial}{\partial z^{\mathrm{R}}}\right), $

其中Δx、Δz分别是法向坐标x和z轴网格间距，

$\Delta r=\sqrt{\Delta x^2+\Delta z^2}$

是旋转坐标网格间距.

相较一般的交错网格，Saenger设计的旋转交错网格法其空间偏导数的计算精度得到提高，其时间偏导数的离散依旧采用二阶差分，低精度的差分格式往往出现时间频散问题(Xie and He, 2024).因此，提高时间离散精度也备受学者关注(Gao and Huang, 2017; Xie and He, 2024).为表述方便，统一将Saenger旋转交错网格法简称为S-RSFD.采用S-RSFD算法计算弹性波方程时，每个空间算子参与计算的网格点数量为.在S-RSFD方案基础上，Xie和He(2024)在空间算子中增加离轴网格点来减少数值时间频散，提出一种时空高阶旋转交错网格格式，简称Xie-RSFD.在计算弹性波方程时，因增加了4^Q－1个计算网格点(Q为空间算子离轴网格点的长度)，Xie-RSFD每个空间算子所需网格点数量提高为4P+4^Q－1.Xie-RSFD算法采用泰勒级数展开来计算差分算子系数，并基于最小二乘法对系数进行优化，能有效模拟不同介质的非均匀性并提高了时间偏导数的计算精度.在旋转交错网格中，选择恰当的方案对差分算子系数进行计算可以有效提升算法性能，因此，本文对时空高阶差分算子系数进行优化，进一步优化提升旋转交错网格在强各向异性介质中的频散抑制作用.

1.2　差分算子系数计算

泰勒级数展开法(Taylor Series Expansion, TE)是时空高阶旋转交错网格差分算子系数计算常用方法之一，其将标准交错网格系数直接引入各向异性介质模拟框架，理论严谨性存疑(Xie and He, 2024).旋转交错网格差分系数的计算通常需要充分考虑时间偏导数的多向性特征，而单纯依赖泰勒级数展开的系数求解方法难以准确满足各向异性介质数值模拟的精度要求.为此，学者们提出了多种差分算子系数的优化方法，如时空域色散关系的方法(Liu and Sen, 2011)、模拟退火法(Zhang and Yao, 2013)、最小二乘法(Least Squares, LS)(Liu, 2014)以及采样逼近法(Sampling Approximation, SA)(Yang et al., 2015)等.

采样逼近法是通过选取特定波数范围内的采样点，构建空间偏导数近似误差最小化模型，获得差分算子系数的近似解的方法.该方法基于空间偏导数计算精度主要取决于k_x和k_z(波数的x和z分量)叠加效应的理论特性，能够在相同计算精度前提下有效降低旋转交错网格的频散误差(Yang et al., 2015).值得注意的是，Yang等(2015)系统评估了TE、LS、SA三种差分算子系数计算方法在旋转交错网格中的性能差异，其实验结果表明在有效波数范围内，采样逼近法与最小二乘法的优化格式均能展现出显著的频散抑制优势，采样逼近法在低波数区域的表现尤为突出，这为改善时空域高阶旋转交错网格的频散特性提供了新的优化方向.

因此，本文将选择泰勒展开法和采样逼近法作为时空高阶精度空间差分算子系数的优化方法，建立适用于各向异性介质速度-应力弹性波方程的优化格式.该优化算法记为TE+SA-RSFD.TE+SA-RSFD算法的核心思想是通过理论约束与数值优化结合计算差分系数.具体分三步实现：(1)基于旋转交错网格离散方程，推导时空频散关系；(2)泰勒展开建立系数理论约束；(3)在波数域采样频散关系构建方程组，数值求解最优系数.下面以TTI介质弹性波方程为例，推导TE+SA-RSFD差分系数的计算过程.

TTI介质中，TE+SA-RSFD的时间偏导数离散格式延续Saenger-RSFD的时间偏导数离散格式，计算空间偏导数时采用Xie-RSFD的空间算子模板(Xie and He, 2024)，为表述方便，统一用u表示各应力和速度分量，即：

$\left\{\begin{align*}\frac{\partial u_{m, n}^{l}}{\partial t} \approx & \boldsymbol{D}_{t} u_{m, n}^{l}=\frac{1}{\Delta t}\left(u_{m, n}^{l+\frac{1}{2}}-u_{m, n}^{l-\frac{1}{2}}\right) \\\frac{\partial u_{m, n}^{l}}{\partial x} \approx & \tilde{\boldsymbol{D}}_{x} R u_{m, n}^{l}=\frac{1}{2 \Delta x}\left[\sum\limits_{p=1}^{P} a_{p}\left(u_{m+\Delta P, n+\Delta P}^{l+\frac{1}{2}}-u_{m-\Delta P, n-\Delta P}^{l+\frac{1}{2}}\right)-\sum\limits_{p=1}^{P} a_{\mathrm{p}}\left(u_{m-\Delta P, n+\Delta P}^{l+\frac{1}{2}}-u_{m+\Delta P, n-\Delta P}^{l+\frac{1}{2}}\right)+\right. \\& \sum\limits_{p=1}^{Q-1} \sum\limits_{q=1}^{Q-p} a_{p, q}\left(u_{m+\Delta Q_{1}, n+\Delta Q_{2}}^{l+\frac{1}{2}}-u_{m-\Delta Q_{2}, n-\Delta Q_{1}}^{l+\frac{1}{2}}+u_{m+\Delta Q_{2}, n+\Delta Q_{1}}^{l+\frac{1}{2}}-u_{m-\Delta Q_{1}, n-\Delta Q_{2}}^{l+\frac{1}{2}}\right)- \\& \left.\sum\limits_{p=1}^{Q-1} \sum\limits_{q=1}^{Q-p} a_{p, q}\left(u_{m-\Delta Q_{1}, n+\Delta Q_{2}}^{l+\frac{1}{2}}-u_{m+\Delta Q_{2}, n-\Delta Q_{1}}^{l+\frac{1}{2}}+u_{m-\Delta Q_{2}, n+\Delta Q_{1}}^{l+\frac{1}{2}}-u_{m+\Delta Q_{1}, n-\Delta Q_{2}}^{l+\frac{1}{2}}\right)\right] \\\frac{\partial u_{m, n}^{l}}{\partial z} \approx & \tilde{\boldsymbol{D}}_{z^{R}} u_{m, n}^{l}=\frac{1}{2 \Delta z}\left[\sum\limits_{p=1}^{P} a_{\mathrm{p}}\left(u_{m+\Delta P, n+\Delta P}^{l+\frac{1}{2}}-u_{m-\Delta P, n-\Delta P}^{l+\frac{1}{2}}\right)+\sum\limits_{p=1}^{P} a_{\mathrm{p}}\left(u_{m-\Delta P, n+\Delta P}^{l+\frac{1}{2}}-u_{m+\Delta P, n-\Delta P}^{l+\frac{1}{2}}\right)+\right. \\& \sum\limits_{p=1}^{Q-1} \sum\limits_{q=1}^{Q-p} a_{p, q}\left(u_{m+\Delta Q_{1}, n+\Delta Q_{2}}^{l+\frac{1}{2}}-u_{m-\Delta Q_{2}, n-\Delta Q_{1}}^{l+\frac{1}{2}}+u_{m+\Delta Q_{2}, n+\Delta Q_{1}}^{l+\frac{1}{2}}-u_{m-\Delta Q_{1}, n-\Delta Q_{2}}^{l+\frac{1}{2}}\right)+ \\& \left.\sum\limits_{p=1}^{Q-1} \sum\limits_{q=1}^{Q-p} a_{p, q}\left(u_{m-\Delta Q_{1}, n+\Delta Q_{2}}^{l+\frac{1}{2}}-u_{m+\Delta Q_{2}, n-\Delta Q_{1}}^{l+\frac{1}{2}}+u_{m-\Delta Q_{2}, n+\Delta Q_{1}}^{l+\frac{1}{2}}-u_{m+\Delta Q_{1}, n-\Delta Q_{2}}^{l+\frac{1}{2}}\right)\right]\end{align*},\right.$

(2)

其中：

$\Delta P=p-\frac{1}{2}, \Delta Q_{1}=p-q-\frac{1}{2} ; \Delta Q_{2}=p+q-\frac{1}{2} ; u_{m, n}^{l}=u(m, n, l) ; ~ \tilde{\boldsymbol{D}}_{x^{\mathrm{R}}} 、\tilde{\boldsymbol{D}}_{z^{\mathrm{R}}}$

为时域高阶旋转交错网格的空间偏导数算子；

$a_{p} 、a_{p, q}$

为空间算子系数.

基于旋转交错网格离散方程和平面波理论对频散关系的推导过程见附录A，频散关系式为

$\left[-2 \sin \left(\frac{\omega \Delta t}{2}\right)\right]^{2}=r^{2}\left(A^{2}+B^{2}\right).$

(3)

下面介绍采用泰勒展开法求解TE+SA-RSFD空间模板中离轴网格点的差分系数a_{p, q}.波数与差分系数关系式的推导过程见附录A，比较波数与差分系数关系式(A5)两边(k_xh)^2l和(k_zh)^2l的系数，得到方程组：

$\begin{array}{l}\sum\limits_{m=1}^{l} d_{\mathrm{m}} d_{l-m+1}= & \sum\limits_{m=1}^{l} \frac{(-1)^{l-1} r^{2 l-2}}{2^{2 l}(2 m-1)!(2 l-2 m+1)!} \\& (1 \leqslant l \leqslant P), \end{array}$

(4)

对方程组(4)求解可得dm.与此同时，比较关系式(A5)两侧k^2l－2ξxk^2ξzh^2l的系数，得到方程组：

$\begin{array}{l}\sum\limits_{i, j, w, v=1,l=i+j+v+w-1}^{Q-1} e_{i, j} e_{w, v}+2 \sum\limits_{i, j, m=1, l=i+j+m-1}^{Q-1} e_{i, j} d_{\mathrm{m}} \\\;\;\;\;\;\;\;\;\;\;= 2 \sum\limits_{m=1}^{l} \frac{(-1)^{l-1} r^{2 l-2}}{2^{2 l}(2 m-1)!(2 l-2 m+1)!} \\\;\;\;\;\;\;\;\;\;\; (1 \leqslant l \leqslant Q-1),\end{array}$

(5)

对方程组(5)求解可得e_{i, j}.最后，再将e_{i, j}代入式(A7)可求得a_{p, q}.

对于TE+SA-RSFD空间模板中其余点的差分系数ap，采用采样逼近法求得.例如，式(A2)中的k_xh=A，k_xh是空间偏导数化简剩余项，A是旋转交错网格离散格式化简的剩余项，为实现方程右端项对左端项的渐近逼近，通过采用采样逼近法控制k_xh∈[0, m]内的误差.具体计算步骤以沿x分量的偏导数为例，在[0, m]内等间距地取P个采样点，用β(1), β(2), …, β(P)表示.当k_z=0时，式(3)中仅与k_x相关的频散关系表示为

$\begin{align*}\sin \left(\frac{\omega \Delta t}{2}\right)= & r \sum\limits_{p=1}^{P} a_{p} \sin \left[\left(p-\frac{1}{2}\right) k_{x} h\right]+2 r \sum\limits_{p=1}^{Q-1} \sum\limits_{q=1}^{Q-P} a_{p, q} \\& \sin \left[\left(p-\frac{1}{2}\right) k_{x} h\right] \cos \left(q k_{x} h\right) .\end{align*}$

(6)

通过在频散关系(6)中用k_z替换k_x，可以变为与k_z有关的频散关系.根据仅与k_x相关的频散关系建立以下误差函数：

$\max\limits_{k_x h \in[0, m]}\left|\sum\limits_{p=1}^p a_{\mathrm{p}} r \frac{\sin \left[(p-1 / 2) k_x h\right]}{\sin \left[0.5 k_x h r\right]}+2 \sum\limits_{p=1}^{Q-1} \sum\limits_{q=1}^{Q-p} a_{p, q} r \frac{\sin \left[( p - 1 / 2 ) k _ { x } h \left[\cos \left[q k_x h\right]\right.\right.}{\sin \left[0.5 k_x h r\right]}-1\right| \leqslant \eta,$

其中，η是最大允许误差，u是k_xh的上限.在最大允许误差趋于零的情况下：

$\begin{gathered}1-2 r \sum\limits_{p=1}^{Q-1} \sum\limits_{q=1}^{Q-p} a_{p, q} \frac{\sin \left[(p - 1 / 2) k _ { x } h \left[\cos \left[q k_{x} h\right]\right.\right.}{\sin \left[0.5 k_{x} h r\right]} \\\approx r \sum\limits_{p=2}^{P} a_{\mathrm{p}} \frac{\sin \left[(p-1 / 2) k_{x} h\right]}{\sin \left[0.5 k_{x} h r\right]},\end{gathered}$

相应地，对于

$k_{x} h=\beta(n), n=1, \cdots, P$

，使用这些采样点满足

$P$

个线性方程:

$\begin{array}{rl}1-2 & r \sum\limits_{p=1}^{Q-1} \sum\limits_{q=1}^{Q-p} a_{p, q} \frac{\sin [(p-1 / 2) \beta(n)[\cos [q \beta(n)]}{\sin [0.5 \beta(n) r]} \\\approx & r \sum\limits_{p=2}^{P} a_{p} \frac{\sin [(p-1 / 2) \beta(n)]}{\sin [0.5 \beta(n) r]}, n=1, \cdots, P ,\end{array}$

(7)

将由泰勒展开法求得的

$a_{p, q}$

代入线性方程组(7)求解，可得基于采样逼近法的旋转交错网格系数

$a_{1}$

，

$a_{2}, \cdots, a_{\mathrm{P}}$

.计算弹性波方程时，因TE+SA－Xie－RSFD采用Xie－RSFD的空间算子模板，每个空间算子所需网格点数为

$4^{P}+4^{Q-1}$

2　频散及稳定性分析

收起

一般地，基于平面简谐波理论，可利用傅里叶分析和相速度相对误差推导各算法的频散关系式与稳定性条件. 为评估算法的频散特性及验证其数值稳定性，基于Madariaga的交错网格(M-SFD)、Saenger的旋转交错网格(S-RSFD)、采用不同系数计算方法的时空域高阶RSFD(TE-Xie-SFD、TE+LS-Xie-SFD)和本文改进优化算法TE+SA-RSFD的频散关系式和稳定性准则进行讨论及建模稳定性比较.

2.1　频散分析

在基于网格离散化的弹性波场数值模拟中，交错网格通过对连续弹性波方程进行离散近似，不可避免地引入计算误差.该误差会导致波场传播过程中群速度与相速度产生频率依赖性，由此产生的非物理频散现象称为数值频散.数值频散作为有限差分法的固有属性，无法被彻底消除，但可通过优化算法设计实现有效压制(Yang et al., 2015).分析数值频散效应通常基于频率-波数域的频散关系式，该关系式可揭示相速度相对误差

$\frac{V_{\mathrm{fd}}}{V}$

的分布特征及其频散强度.交错网格的频散关系式为

$\delta_{\mathrm{fd}}=\frac{V_{\mathrm{fd}}}{V}=\frac{2}{k h r} \arcsin \left[r \sqrt{\left(\sum\limits_{p=1}^{P} c_{\mathrm{p}} \sin \left[\left(p-\frac{1}{2}\right) k_{\mathrm{x}} h\right]\right)^{2}+\left(\sum\limits_{p=1}^{P} c_{\mathrm{p}} \sin \left[\left(p-\frac{1}{2}\right) k_{\mathrm{z}} h\right]\right)^{2}}\right], $

(8)

相比之下，旋转交错网格通过优化算子结构，其频散关系式展现出更优的波场保真特性.其频散关系式为

$\left\{\begin{array}{l}\delta_{\mathrm{fd}}=\frac{V_{\mathrm{fd}}}{V}=\frac{1}{k h r} \arccos \left[1-\frac{r^{2}\left(I^{2}+J^{2}\right)}{2}\right] \\I=\sum\nolimits_{p=1}^{P} b_{\mathrm{p}}\left\{\sin \left[\left(p-\frac{1}{2}\right) k_{\mathrm{x}} h+\left(p-\frac{1}{2}\right) k_{\mathrm{z}} h\right]+\sin \left[\left(p-\frac{1}{2}\right) k_{\mathrm{x}} h-\left(p-\frac{1}{2}\right) k_{\mathrm{z}} h\right]\right\} \\J=\sum\nolimits_{p=1}^{P} b_{\mathrm{p}}\left\{\sin \left[\left(p-\frac{1}{2}\right) k_{\mathrm{x}} h+\left(p-\frac{1}{2}\right) k_{\mathrm{z}} h\right]-\sin \left[\left(p-\frac{1}{2}\right) k_{\mathrm{x}} h-\left(p-\frac{1}{2}\right) k_{\mathrm{z}} h\right]\right\}\end{array}\right.,$

(9)

特别地，时空高阶旋转交错网格通过引入时间-空间耦合差分格式，进一步提升了数值频散压制效果，其频散关系式为

$\delta_{\mathrm{fd}}=\frac{V_{\mathrm{fd}}}{V}=\frac{1}{k h r} \arccos \left[1-\frac{r^{2}\left(A^{2}+B^{2}\right)}{2}\right], $

(10)

式(8)~(10)频散关系式中，

$c_{\mathrm{p}}$

和

$b_{\mathrm{p}}$

分别表示交错网格与旋转交错网格差分算子系数.

$\delta_{\mathrm{fd}}$

表示相速度的相对误差值，若

$\delta_{\mathrm{fd}}>1$

则表示出现时间频散，相位提前；若

$\delta_{\mathrm{fd}}<1$

，则表示出现空间频散，相位延迟(Zhou et al.，2022).为了直观呈现各算法的频散程度，根据频散关系式(8)~(10)绘制交错网格M-SFD、旋转交错网格S-RSFD、TE-Xie-RSFD、TE+LS-Xie-SFD和本文改进算法TE+SA-RSFD的频散曲线图.

为研究空间算子长度对旋转交错网格算法S-RSFD、TE-Xie-RSFD和TE+SA-RSFD频散程度的影响，讨论下面两种情形下的各算法的频散程度：

(1) 取定参数r=0.4，传播角

$\theta=\frac{{\rm{\mathsf{π}}}}{8}$

，讨论各旋转交错网格算法在空间算子长度取2~8情形下，相速度的相对误差随kh变化情况，其结果如图 1所示.图中红色虚线表示相速度相对误差δ_fd的基准线(δ_fd=1)，频散曲线与基准线距离越接近，表示频散程度越弱.

(2) 取定r=0.4，kh=2.2，S-RSFD、TE-Xie-RSFD和TE+SA-RSFD在空间算子长度取2~8情形下，各算法相速度的相对误差随传播角θ变化的频散情况如图 2所示.

由图 1可知，随着kh∈[1, 3]的增加，S-RSFD、TE-Xie-RSFD和TE+SA-RSFD的频散曲线均随空间算子长度P的增大越来越接近基准线.由此可知，随空间精度阶数O(h^2P)的提高，各算法的频散程度逐渐减小，但当空间算子长度增加到一定的长度后，各算法的频散程度没有显著变化，即随着空间算子精度的提高频散并没有得到很好地抑制.由图 1a、b、c，可以看出，在kh∈[1, 2.1]之间时，S-RSFD算法存在时间频散现象，随着空间长度P越大时间频散程度越小而空间频散程度越大；TE-Xie-RSFD通过修改空间差分算子提高了交错网格的时间计算精度，时间频散程度较低，但是空间频散较大；改进算法TE+SA-RSFD在低空间算子长度(P < 4)情形下，其频散曲线较S-RSFD时间频散程度较小，较TE-Xie-RSFD空间频散明显变小，但与S-RSFD一样存在时间频散.

由图 2可知，随着传播角

$\theta \in\left[\frac{{\rm{\mathsf{π}}}}{60}, \frac{{\rm{\mathsf{π}}}}{5}\right]$

的增加，S-RSFD、TE-Xie-RSFD和TE+SA-RSFD的频散曲线均随空间算子长度P的增大趋近于基准线，这表明提高空间差分精度可有效降低算法频散误差，但当P≥7后，继续增加算子长度无显著频散抑制效果.由图 2a—c，可以看出，当θ∈ π 60, π 12之间时，S-RSFD算法存在时间频散现象，随P值增大时间频散程度减小，空间频散程度呈现反向增大趋势；TE-Xie-RSFD在时间频散控制方面表现优异，但空间频散程度较大；改进算法TE+SA-RSFD的空间频散程度较TE-Xie-RSFD显著减小，在

$\theta \in\left[\frac{{\rm{\mathsf{π}}}}{60}, \frac{{\rm{\mathsf{π}}}}{12}\right]$

之间时，时间频散程度较S-RSFD减小.

下面讨论同一空间算子长度P=5情形下，取定参数r=5.33，讨论不同算法的相速度相对误差随kh值变化情况，这时传播角θ取为

$\frac{{\rm{\mathsf{π}}}}{3}$

，其结果如图 3a所示；讨论不同算法的相速度相对误差随传播角θ变化的频散曲线如图 3b所示，这时kh值取为2.2.

由图 3a可知，交错网格法M-SFD的时间频散程度较大，随着kh的增加会出现空间频散；S-RSFD与M-SFD的频散变换较为一致，但S-RSFD在kh取1.8左右开始出现空间频散；TE-Xie-RSFD的空间频散程度较大，在研究范围内未出现时间频散；在kh较小时，TE+LS-Xie-RSFD时间频散程度较小，随着kh的增大其存在较明显的空间频散；在kh较小时，改进算法TE+SA-RSFD的频散程度较小，当kh>1后，其空间频散程度增大，但较TE-Xie-RSFD、TE+LS-Xie-RSFD两种算法，对频散具有更好的抑制作用.

由图 3b可知，随着

$\theta$

的增加，M-SFD无空间频散现象，而时间频散程度较大；S-RSFD在

$\theta \in\left[\frac{{\rm{\mathsf{π}}}}{60}, \frac{{\rm{\mathsf{π}}}}{10}\right]$

之间时间频散程度较大，在

$\theta$

取

$\frac{{\rm{\mathsf{π}}}}{10}$

左右开始出现空间频散；TE-Xie-RSFD的空间频散程度随着θ增大存在较明显的空间频散，在研究范围内未出现时间频散；改进算法TE+SA-RSFD的频散曲线和TE+LS-Xie-RSFD的频散曲线在θ∈[

$\frac{{\rm{\mathsf{π}}}}{60}$

$\frac{{\rm{\mathsf{π}}}}{10}$

]之间基本重合，当θ>

$\frac{{\rm{\mathsf{π}}}}{10}$

后，TE+SA-RSFD较TE+LS-Xie-RSFD空间频散程度略大.

综上所述，在同一空间算子长度情形下，TE+SA-RSFD方案时间频散程度较M-SFD、S-RSFD和TE+LS-Xie-RSFD小，空间频散程度较TE-Xie-RSFD小.

取定

$r=5.33$

，波数范围[0.3，2.3]、传播角

$\left[\frac{{\rm{\mathsf{π}}}}{60}, \frac{{\rm{\mathsf{π}}}}{3}\right]$

，三个不同差分系数计算方法下的各算法的归一化相速度误差随波数与传播角度变化的等高线分布热图，如图 4所示.由图 4可知: TE-Xie-RSFD方案归一化相速度误差值低于0.05的区域占比小于

$50 \%$

，误差值小于0.01的区域占比约

$25 \%$

；TE+LS-Xie-RSFD方案在相同参数区间内，误差值低于0.03的区域占比提升至约

$50 \%$

，但误差值小于0.01的区域占比与TE-Xie-RSFD基本持平，整体精度呈微弱优势；TE+SA-RSFD方案表现显著优化，误差值低于0.05的区域占比约

$70 \%$

，误差值小于0.01的区域占比提升至

$40 \%$

，可将归一化相速度误差有效控制在

$10 \%$

以下.

综上分析，在波数和入射角参数变化域内，差分算子系数优化方案TE+SA-RSFD较TE-Xie-RSFD、TE+LS-Xie-RSFD具有更优的频散抑制能力，其归一化相速度误差小于0.01的区域范围最大.

2.2　稳定性分析

在弹性波方程有限差分数值离散体系中，稳定性条件表征了地震波传播速度、时空采样步长、差分算子长度及差分系数间的约束关系.为构建改进算法的稳定性理论框架，基于离散化弹性波方程的频散关系式(8)~(10)，根据

$0 \leqslant \sin ^{2}\left(\frac{r k h}{2}\right) \leqslant 1$

，且设定

$k_{x}=k_{z}=\frac{{\rm{\mathsf{π}}}}{2 h}$

，经推导获得稳定性判据不等式.该判据描述了差分格式在数值计算过程中保持收敛性所需满足的时空参数阈值.值得注意的是，不同网格离散策略对应的稳定性准则具有差异.对于交错网格有限差分法，其稳定性条件可表述为

$\Delta t v_{p} \frac{1}{h}\left(\sqrt{2} \sum\nolimits_{p=1}^{P}\left|c_{p}\right|\right) \leqslant 1 .$

(11)

而采用旋转交错网格架构的旋转交错网格方法，通过优化应力-速度分量耦合关系，其稳定性准则改进为

$\Delta t v_{p} \frac{1}{h}\left(\sum\nolimits_{p=1}^{P}\left|b_{p}\right|\right) \leqslant 1 .$

(12)

进一步地，时空高阶旋转交错网格在扩展差分算子维度的同时，通过构建时空耦合加权系数矩阵，发展出具有更高容许度的稳定性判据，其稳定性准则为

$\Delta t v_{p} \frac{1}{h}\left(\sum\nolimits_{p=1}^{P}\left|a_{p}\right|+2 \sum\nolimits_{p=1}^{Q-1} \sum\nolimits_{q=1}^{Q-p}\left|a_{p, q}\right|\right) \leqslant 1 . $

(13)

为了对比各算法符合稳定条件的时空间隔区间，根据稳定性准则(11)~(13)绘制M-SFD、S-RSFD、TE-Xie-RSFD、TE+LS-Xie-SFD和本文改进算法TE+SA-RSFD的稳定曲线图.

为研究空间算子长度对交错网格算法M-SFD、旋转交错网格算法S-RSFD和改进算法TE+SA-RSFD符合稳定条件的时空间隔区间的影响，讨论在空间算子长度取2~8时，下面两种情形下的各算法的稳定曲线：

(1) 取定空间间隔h为15 m，V_p=4000 m/s，讨论各算法稳定值随时间间隔Δt的变化情况，其结果如图 5所示.图中红色虚线为稳定值阙值1，各算法稳定曲线位于虚线下方时，即稳定值低于1时满足稳定性条件，反之则失稳.

(2) 取定时间间隔Δt为0.1 ms，V_p=4000 m/s，讨论各算法稳定值随空间间隔h的变化情况如图 6所示.

由图 5可知，M-SFD、S-RSFD和TE+SA-RSFD的稳定性条件均与时间间隔呈正相关，随着空间算子长度P的增大，各算法的稳定性越差，由此可知，高精度离散将加剧算法时间间隔的稳定性约束，但当空间算子长度增加到一定的长度后，再度提高空间算子精度，满足稳定条件的时间间隔区间长度无显著变化.由图 5a、b、c，可以看出，当Δt大于0.19 ms时，M-SFD算法满足稳定性条件；S-RSFD算法通过修改空间差分算子增大了符合稳定条件的时间间隔区间长度，在空间算子长度较大(P=8)情形下，在Δt=0.27 ms处突破稳定阙值；改进算法TE+SA-RSFD在Δt < 0.27 ms之间满足稳定性条件.由此可知，改进算法TE+SA-RSFD较M-SFD和S-RSFD可以在更宽泛的时间间隔变化范围内维持算法稳定性.

由图 6可知，M-SFD、S-RSFD和TE+SA-RSFD的稳定性条件均与空间间隔呈负相关，且随空间算子长度P的增大，各算法的稳定性越差，当P≥7后，继续增加算子长度，各算法的稳定性空间间隔约束区间均无明显变化.由图 6a—c可以看出，当在h∈[0.3, 0.75]之间，M-SFD算法失稳；h∈[0.3, 0.55]时，S-RSFD大于稳定阙值；h∈[0.3, 0.5]时，改进算法TE+SA-RSFD不满足稳定性条件.由此可知，改进算法TE+SA-RSFD较M-SFD和S-RSFD可以在更宽泛的空间间隔变化范围内维持算法稳定性.

同一空间算子长度P=5情形下，取定V_p=4000 m/s，讨论相同空间步长h=15 m情形下，各算法取不同时间步长的稳定性变化情况，其结果如图 7a所示.同时，讨论相同时间步长Δt=0.1 ms情形下，各算法取不同空间间隔大小的稳定性变化，如图 7b所示.

由图 7a可知，固定空间步长情形下，时间步长越小稳定性越好.当时间间隔Δt大于0.2 ms，交错网格法M-SFD不稳定，M-SFD稳定性准则(11)左端项存在系数

$\sqrt{2}$

的放大效应，其速度-空间步长比值更易突破临界阈值，导致稳定性条件失效；当时间间隔Δt小于0.28 ms时，S-RSFD符合稳定性条件；TE-Xie-RSFD可取更大的时间间隔，仍满足稳定性条件；当时间间隔Δt小于0.29 ms时，TE+LS-Xie-RSFD、TE+SA-RSFD满足稳定性条件.由此可知，考虑计算效率时，可适当增大时间步长以提高计算效率，TE-Xie-RSFD的最优，TE+LS-Xie-RSFD与本文改进算法TE+SA-RSFD次之.

由图 7b可知，固定时间步长Δt=0.1 ms情形下，空间步长越大则稳定性越好.M-SFD在空间间隔h∈[0.3, 0.8]失稳，失稳区间长度显著大于其他算法；当空间间隔h∈[0.3, 0.55]之间，S-RSFD算法失稳；h∈[0.3, 0.4]时TE-Xie-RSFD大于稳定阙值；h∈[0.3, 0.52]时，TE+LS-Xie-RSFD失稳；在h≤0.52时，改进算法TE+SA-RSFD不满足稳定性条件.由此可知，考虑计算精度时，尽可能取较小的空间步长，此时TE-Xie-RSFD的最优，TE+LS-Xie-RSFD与本文改进算法TE+SA-RSFD次之.

综上频散以及稳定性分析，改进算法TE+SA-RSFD通过系数优化，能够在更大幅度的时空步长调整中保持稳定运行，相较于传统交错网格(M-SFD)和传统旋转网格方法(S-RSFD)有明显优势.与此同时，相较于稳定性条件更强的TE-Xie-RSFD算法，TE+SA-RSFD具有比其更低的频散程度.因此，TE+SA-RSFD算法既具有较优的频散抑制能力，又能适应更灵活的步长设置，实现了频散与稳定性之间的优化平衡.

3　数值模拟实验

收起

为验证改进算法TE+SA-RSFD的数值性能，选取基于Madariaga交错网格法(M-SFD)、Saenger旋转交错网格法(S-RSFD)以及采用不同系数计算策略的时空域高阶RSFD方案(TE-Xie-RSFD、TE+LS-Xie-RSFD)进行对比，对问题(1)进行数值模拟.实验模型为：(1)各向同性介质模型；(2)各向异性介质模型：各向异性VTI介质模型、各向异性TTI介质模型；(3)修正BP模型.数值实验中，在各类模型几何中心设置雷克子波震源进行波场模拟实验，观测在介质复杂度递进场景下波场快照的波形特征.为了消除边界反射干扰现象，三类实验均采用PML边界条件，衰减函数为

$d_i=\frac{-3 \log (s) v_{\mathrm{p}} L_i^2}{2 L^2}$

，其中，d_i表示x轴方向或者z轴方向的衰减函数；s=0.000001；L_i表示计算区域到匹配层边界的距离；L表示匹配层的厚度.

3.1　各向同性介质模型

设计一个均匀各向同性介质模型，模型参数：纵、横波速度和密度分别为2500 m/s，1470 m/s，2.5 g/cm³，模型尺寸nx×nz=400×400，模型空间采样间隔h=Δx=Δz=3 m，时间采样间隔Δt=0.1 ms，边界取模型四周300 m.在模型中心(x, z)=(600 m, 600 m)处释放60 Hz主频的雷克子波作为震源.为保证各算法离散精度一致，TE-Xie-RSFD、TE+LS-Xie-RSFD与TE+SA-RSFD的空间算子离轴网格点长度Q均取3.图 8和图 9分别给出了四种算法((M-SFD、S-RSFD、TE-Xie-RSFD、TE+LS-Xie-RSFD)及改进算法(TE+SA-RSFD)在空间离散精度4阶、18阶(算子长度P=2和P=9)条件下模拟得到的压力场的x分量波场快照.

通过图 8和图 9中五种交错网格波场快照的对比分析，发现：

(1) 当空间离散精度O(h^2P)由O(h⁴)提升至O(h¹⁸)时，M-SFD、S-RSFD、TE-Xie-RSFD、TE+LS-Xie-RSFD和TE+SA-RSFD算法产生的波场快照数值频散程度均存在衰减现象，验证了高精度格式的频散抑制机理.

(2) 对比图 8a、b和图 9a、b，由于M-SFD空间算子中的计算网格点数量少于S-RSFD，部分网格点不参与计算，导致M-SFD的波场快照较S-RSFD波形间断区域更大，存在波前断裂现象.

(3) 如图 9b、c所示，在高精度机制下(P=9)，相较于S-RSFD算法产生的波场快照，TE-Xie-RSFD通过扩展空间差分算子模板增加了参与计算的网格点数，使其波场快照中的波前断裂区域减少，与此同时，TE-Xie-RSFD波场中的频散现象较S-RSFD增强.

(4) 由图 9c、d可知，TE+LS-Xie-RSFD通过LS法优化系数使其波场快照的频散程度相较TE-Xie-RSFD降低，但受最小二乘约束产生了波形几何畸变(波形矩形化)，同时，波形间断区域增加.

(5) 如图 9c—e所示，改进算法TE+SA-RSFD综合了平衡频散程度与波形保真度，其波场快照的频散现象较TE-Xie-RSFD减弱，波形分辨率与TE+LS-Xie-RSFD相当.

3.2　各向异性介质模型

为了验证优化算法应用于各向异性弹性波方程中的模拟效果，分别设计一个各向异性VTI介质模型和一个复杂各向异性TTI介质模型对比各算法在二维三分量弹性波方程求解的数值频散程度.VTI介质模型参数：纵、横波速度和密度分别为2500 m/s，1470 m/s，2.5 g/cm³，Thomsen参数ε=0.3333，δ=0.885，γ=0.25；TTI介质模型参数：纵、横波速度和密度分别为2500 m/s，1470 m/s，2.5 g/cm³，Thomsen参数ε=0.3333，δ=0.885，γ=0.25，TTI介质模型是由VTI介质模型通过逆时针绕z轴旋转60度产生.两个模型尺寸都设置为：nx×nz=400×400，空间采样间隔h=Δx=Δz=3 m，时间采样间隔Δt=0.1 ms，边界取模型四周300 m.在两个模型中心(x, z)=(600 m, 600 m)处释放60 Hz主频的雷克子波作为震源.为保证各算法离散精度一致，五种算法的空间算子长度P均取9(空间离散精度为18阶)，TE-Xie-RSFD、TE+LS-Xie-RSFD与TE+SA-RSFD的空间算子离轴网格点长度Q均取3.

图 10是五种交错网格M-SFD、S-RSFD、TE-Xie-RSFD、TE+LS-Xie-RSFD和TE+SA-RSFD基于VTI模型模拟得到的压力场x分量的波场快照，通过对比分析发现：

(1) 如图 10a、b所示，在基于VTI介质模型的数值实验中，M-SFD算法产生的波场快照呈现单一纵波椭圆波前特征，而S-RSFD方法表现出更复杂的各向异性响应：除纵波椭圆波前外，可观测到具有横向椭圆特征的快横波与慢横波分裂现象.

(2) 由图 10b、c可知，TE-Xie-RSFD方法相较于S-RSFD方法，其波场快照的外围波形分辨率提升，并能有效揭示介质内部各向异性结构特征，同时，该方法伴随更强烈的数值频散现象，表明了计算精度与数值稳定性之间存在权衡关系.

(3) 如图 10c、d所示，TE+LS-Xie-RSFD较TE-Xie-RSFD算法在验中有效抑制了数值频散，分辨率进一步提升，但导致其快、慢横波波形连续性降低，波前断裂现象增加，推测此现象与最小二乘约束导致的相位畸变相关.

(4) 如图 10e所示，改进算法TE+SA-RSFD产生的波场快照完整保留了介质的各向异性特征，分辨率与TE+LS-Xie-RSFD产生的波场快照相当.

图 11是各算法基于TTI模型的波场快照，如图 11d、e可知，改进算法TE+SA-RSFD产生的波场快照在TTI介质模型实验中完整保留了介质的各向异性特征，分辨率较TE+LS-Xie-RSFD相当，另外，TE+SA-RSFD的波场快照较TE+LS-Xie-RSFD数值频散现象显著减弱，且未出现明显波前断裂现象，验证了改进算法TE+SA-RSFD在各向异性介质数值实验中的频散抑制优势.

为了获得更加直观的观察效果，图 12是以传统交错网格M-SFD(图中红线)为参照，各算法基于各向异性介质模型下在坐标点z=270处的地震波形记录.由图 12a可知，在VTI介质模型实验中，改进算法TE+SA-RSFD的波形与TE+LS-Xie-RSFD几乎一致，波后幅度较均TE+LS-Xie-RSFD略大，存在一定程度的空间频散误差；由图 12b可知，在TTI介质模型实验中，TE+SA-RSFD不管波前还是波后幅度较TE+LS-Xie-RSFD和S-RSFD都更接近M-SFD，空间和时间频散误差较小.综合来说，在各向异性介质模型数值实验中，TE+SA-RSFD具有频散抑制优势.

3.3　修正BP模型

针对复杂各向异性介质中弹性波传播的数值模拟精度问题，基于修正BP模型开展对比实验，通过分析各算法各向异性弹性波方程求解过程中不同算法的数值频散特性，系统验证优化算法的计算优势.图 13是修正后的BP网络模型参数示意图，模型尺寸为nx×nz=1575×900, 模型空间采样间隔h=Δx=Δz=5 m，时间采样间隔Δt=0.1 ms，在模型中心采用40 Hz主频的雷克子波作为震源.为保证各算法离散精度一致，五种算法的空间算子长度P均取9，TE-Xie-RSFD、TE+LS-Xie-RSFD与TE+SA-RSFD中的空间算子额外点的长度Q取3.

图 14是五种交错网格M-SFD、S-RSFD、TE-Xie-RSFD、TE+LS-Xie-RSFD和TE+SA-RSFD分别基于修正BP模型在t=0.25 s、t=0.55 s时刻弹性波数值模拟的波场快照，由图 14可知，

(1) 交错网格与旋转交错网格算法在复杂介质中的波场重构能力呈现显著差异：如图 14a所示，交错网格算法M-SFD波场快照在t=0.25 s时刻因受限于计算域内部分网格节点未参与运算，局部网格节点激活不足导致波前重构缺失；当模拟时间延拓至0.55 s时，其波场快照中无法观测到波形特征，此现象揭示了M-SFD在长时模拟中存在的稳定性缺陷.

(2) 各旋转交错网格算法通过改进空间离散结构，在弹性参数梯度突变区域表现出分辨率优势：如图 14b—e所示，S-RSFD、TE-Xie-RSFD、TE+LS-Xie-RSFD和TE+SA-RSFD在t=0.25 s时刻波场快照中均可辨识典型的不对称畸变特征，且维持了波场传播的时空连续性；随模拟时间延拓至0.55 s时，旋转交错网格算法均可在长时模拟中维持稳定性，同时，如图 14c—e所示，时空高阶旋转交错网格算法TE-Xie-RSFD、TE+LS-Xie-RSFD与TE+SA-RSFD均可清晰分离准纵波(qP)、快横波(qS1)与慢横波(qS2)分量，较S-RSFD算法分辨率提升.

为比较改进算法和其他旋转交错网格算法的数值性能，图 15是S-RSFD、TE-Xie-RSFD、TE+LS-Xie-RSFD和TE+SA-RSFD在t=0.55 s时刻基于修正BP模型波场快照的细节放大图，由图 15可知：

(1) 比较图 15a、b可知，TE-Xie-RSFD通过增加网格点优化计算，其波场快照中箭头处不对称畸变特征较S-RSFD更清晰，不足的是，TE-Xie-RSFD波场快照中伴随强烈数值频散现象，同时虚线圈内的伪影反映了有限差分法中计算精度与数值稳定性之间的固有矛盾.

(2) 比较图 15b、c可知，TE+LS-Xie-RSFD引入最小二乘约束优化技术，伪影在很大程度上被压制，同时，其波场快照中红色箭头处仍存在数值频散现象.

(3) 如图 15e所示，改进算法TE+SA-RSFD波场快照中各向异性波形内部结构分辨率较S-RSFD显著提高，伪影较TE+LS-Xie-RSFD被进一步被压制，且TE+SA-RSFD通过时空耦合离散策略展现出更精细的波型解耦能力，红色箭头处数值频散现象较TE-Xie-RSFD和TE+LS-Xie-RSFD减少，维持了波场传播的时空连续性，验证了采样逼近约束策略对各向异性特征的保护机制.

图 16是各算法基于修正BP模型在坐标点z=370处的地震波形记录.以传统交错网格M-SFD(图中红线)为参照，观察波形曲线可以看到：由图 16可知，在t=0.25 s和t=0.55 s时刻，改进算法TE+SA-RSFD的波前还是波后幅度较TE-Xie-RSFD和TE+LS-Xie-RSFD均更接近M-SFD，也就是说，在修正BP模型数值实验中，TE+SA-RSFD对时空频散的抑制较TE-Xie-RSFD和TE+LS-Xie-RSFD更具优势.

以上实验证明，改进算法TE+SA-RSFD在各向同性介质、VTI介质、TTI介质及复杂各向异性介质中均具有较优波形分辨率，保持各向同性或各向异性波场特征完整性的同时，求解过程中能够很好地压制各种介质弹性波方程数值模拟中的数值频散现象.

4　结论与展望

收起

基于泰勒展开法与采样逼近法协同优化差分算子系数，本文给出了一种时空高阶旋转交错网格有限差分法.该方法通过融合泰勒展开法的理论完备性与采样逼近法的数值适应性，优化空间导数算子系数，显著提升了各向异性弹性波方程的数值求解精度.根据频散分析和稳定性分析结果，改进算法的稳定性得到保证，且能降低数值频散程度，相速度相对误差在10%以内.层状介质、VTI/TTI各向异性模型及修正BP模型的数值实验证实，该方法在各向同性及各向异性介质中均展现出更优的数值频散压制能力，有效解决了传统方法精度与稳定性难以兼顾的瓶颈，为裂缝介质地震波传播机理研究提供了高效计算工具.

交错网格与旋转交错网格的差分系数仅依赖于偏导数，与介质弹性参数无关，而本文的时空高阶旋转交错网格改进方案需将弹性参数与偏微分方程耦合求解，导致时空域高阶差分系数需依据弹性参数和偏导数共同确定.这一特性使得各向异性介质模拟前需逐点计算差分系数，显著增加了内存消耗.未来研究可通过设计混合差分格式，在保证精度的前提下减少存储需求，进一步推动改进方案在地震波场模拟中的实际应用.

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

Almquist M, Dunham E M. 2021. Elastic wave propagation in anisotropic solids using energy-stable finite differences with weakly enforced boundary and interface conditions. Journal of Computational Physics, 424: 109842, doi: 10.1016/j.jcp.2020.109842.

Behera L. 2022. Elastic anisotropic finite-difference full-wave modeling and imaging of 2D tilted transversely isotropic (TTI) media. Journal of Applied Geophysics, 207: 104837, doi: 10.1016/j.jappgeo.2022.104837.

Chen H, Wang X M, Zhao H B. 2006. A rotated staggered grid finite-difference with the absorbing boundary condition of a perfectly matched layer. Chinese Science Bulletin, 51(19): 2304-2314, doi: 10.1007/s11434-006-2127-8.

Doostaki R, Hosseini M M, Salemi A. 2023. A new simultaneously compact finite difference scheme for high-dimensional time-dependent PDEs. Mathematics and Computers in Simulation, 212: 504-523, doi: 10.1016/j.matcom.2023.05.008.

Fei T, Larner K. 1995. Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport. Geophysics, 60(6): 1830-1842, doi: 10.1190/1.1443915.

Gao K, Huang L J. 2017. An improved rotated staggered-grid finite-difference method with fourth-order temporal accuracy for elastic-wave modeling in anisotropic media. Journal of Computational Physics, 350: 361-386, doi: 10.1016/j.jcp.2017.08.053.

Gold N, Shapiro S A, Burr E. 1997. Modelling of high contrasts in elastic media using a modified finite difference scheme. //SEG Technical Program Expanded Abstracts. SEG, 1850-1853.

Khan I R, Ohba R. 2003. Taylor series based finite difference approximations of higher-degree derivatives. Journal of Computational and Applied Mathematics, 154(1): 115-124, doi: 10.1016/S0377-0427(02)00816-6.

Li K P, Yang T, Jiang W, et al. 2023. Isogeometric boundary element method for isotropic damage elastic mechanical problems. Theoretical and Applied Fracture Mechanics, 124: 103802, doi: 10.1016/j.tafmec.2023.103802.

Lin B, Ma Y, Wang C S. 2024. A Lawson-time-splitting extended Fourier pseudospectral method for the Gross-Pitaevskii equation with time-dependent low regularity potential. Journal of Computational Physics, 512: 113133, doi: 10.1016/j.jcp.2024.113133.

Liu Y, Sen M K. 2011. Scalar wave equation modeling with time-space domain dispersion-relation-based staggered-grid finite-difference schemes. Bulletin of The Seismological Society of America, 101(1): 141-159, doi: 10.1785/0120100041.

Liu Y. 2014. Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling. Geophysical Journal International, 197(2): 1033-1047, doi: 10.1093/gji/ggu032.

Lu Y M, Yang J X, Zhang J F, et al. 2025. Traveltime calculation for qP and qSV waves based on fast fitting algorithm in 2D TTI media. Chinese Journal of Geophysics (in Chinese), 68(5): 1911-1921, doi: 10.6038/cjg2024R0795.

Saenger E H, Gold N, Shapiro S A. 2000. Modeling the propagation of elastic waves using a modified finite-difference grid. Wave Motion, 31(1): 77-92, doi: 10.1016/S0165-2125(99)00023-2.

Wang C Q, Wu G C, Yang L Y. 2025. Finite-difference numerical simulation method for two-dimensional scalar waves in a 25-point Laplace-Fourier domain. Progress in Geophysics (in Chinese), 40(6): 2698-2710, doi: 10.6038/pg2025HH0496.

Wang Q, Cheng B J, Xu T J, et al. 2025. Two-dimensional elastic wave numerical simulation in dual-phase HTI porous media. Progress in Geophysics (in Chinese): 1-12, doi: 11.2982.p.20250627.1131.023.

Wei H, Tong S Y, Xin C Q, et al. 2024. Numerical simulation of 2D three-component surface waves in TI media based on the rotated staggered grid. Oil Geophysical Prospecting (in Chinese), 59(5): 1022-1036, doi: 10.13810/j.cnki.issn.1000-7210.2024.05.010.

Xie J R, He B S. 2024. Spatial-temporal high-order rotated-staggered-grid finite-difference scheme of elastic wave equations for TTI medium. Journal of Computational Physics, 498: 112684, doi: 10.1016/j.jcp.2023.112684.

Xu Y G, Liao J P, Zhou L, et al. 2023. Comparing the adaptability of three different staggered-grids for the simulation of 3D anisotropic fractured media and improving the Lebedev method. Chinese Journal of Geophysics (in Chinese), 66(3): 1164-1179, doi: 10.6038/cjg2022Q0089.

Yan H Y, Xu T. 2024. Suppressing spatial dispersion of seismic finite-difference modeling with the improved pix2pix algorithm. Journal of Computational Physics, 511: 113129, doi: 10.1016/j.jcp.2024.113129.

Yang L, Yan H Y, Liu H. 2015. Optimal rotated staggered-grid finite-difference schemes for elastic wave modeling in TTI media. Journal of Applied Geophysics, 122: 40-52, doi: 10.1016/j.jappgeo.2015.08.007.

Yang L, Yan H Y, Liu H. 2017. An optimal implicit staggered-grid finite-difference scheme based on the modified Taylor-series expansion with minimax approximation method for elastic modeling. Journal of Applied Geophysics, 138: 161-171, doi: 10.1016/j.jappgeo.2017.01.020.

Zhang J H, Yao Z X. 2013. Optimized finite-difference operator for broadband seismic wave modeling. Geophysics, 78(1): A13-A18, doi: 10.1190/geo2012-0277.1.

Zhang O M, Schmitt D R. 2025. An optimized 2D/3D finite-difference seismic wave propagator using rotated staggered grid for complex elastic anisotropic structures. Computers & Geosciences, 196: 105850, doi: 10.1016/j.cageo.2024.105850.

Zhao M L, Shi J Y, Wang L Q. 2023. A three-dimensional Petrov-Galerkin finite element interface method for solving inhomogeneous anisotropic Maxwell's equations in irregular regions. Computers & Mathematics with Applications, 152: 364-377, doi: 10.1016/j.camwa.2023.10.035.

Zhong L, Ma S F, Joka D L. 2021. Meshless method for delay partial differential equations with piecewise continuous arguments. Chinese Journal of Computational Mechanics (in Chinese), 38(2): 160-165, doi: 10.7511/jslx20191217001.

Zhou H Y, Liu Y, Wang J. 2022. Time-space domain scalar wave modeling by a novel hybrid staggered-grid finite-difference method with high temporal and spatial accuracies. Journal of Computational Physics, 455: 111004, doi: 10.1016/j.jcp.2022.111004.

Zhou X H, Huo S D, Liang Y, et al. 2023. Finite-difference method for modeling the surface wave propagation with surface topography in anisotropic-viscoelastic media. Journal of Applied Geophysics, 217: 105161, doi: 10.1016/j.jappgeo.2023.105161.

芦永明, 杨建勋, 张剑锋, 等. 2025. 基于快速拟合算法的二维TTI介质qP波和qSV波走时计算. 地球物理学报, 68(5): 1911-1921, doi: 10.6038/cjg2024R0795.

王超群, 吴国忱, 杨凌云. 2025.25点Laplace-Fourier域二维标量波有限差分数值模拟方法. 地球物理学进展, 40(6): 2698-2710, doi: 10.6038/pg2025HH0496.

王啟, 程冰洁, 徐天吉, 等. 2025. 双相-HTI多孔介质二维弹性波数值模拟. 地球物理学进展. 1-12, doi: 11.2982.p.20250627.1131.023.

魏皓, 童思友, 辛成庆, 等. 2024. 基于旋转交错网格的TI介质二维三分量面波模拟. 石油地球物理勘探, 59(5): 1022-1036, doi: 10.13810/j.cnki.issn.1000-7210.2024.05.010.

徐云贵, 廖建平, 周林, 等. 2023. 三维各向异性裂缝介质正演模拟的三种交错网格适应性比较及Lebedev方法的改进. 地球物理学报, 66(3): 1164-1179, doi: 10.6038/cjg2022Q0089.

钟霖, 马淑芳, 莱蒙. 2021. 无网格法求解一类分段连续型延迟偏微分方程. 计算力学学报, 38(2): 160-165, doi: 10.7511/jslx20191217001.

2026年第41卷第1期

PDF下载

189

引用本文

BibTeX

文章信息

doi: 10.6038/pg2026JJ0228

接收时间：2025-06-26
首发时间：2026-09-23
出版时间：2026-02-20

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2025-06-26
修回日期：2025-09-10

基金

作者信息

¹安徽理工大学数学与大数据学院, 淮南 232000

²桂林电子科技大学数学与计算科学学院, 桂林 541004

通讯作者:

赵前进, 男, 1967年生, 教授, 博士, 研究方向为有理插值与逼近、数字图像处理.E-mail: qjzhao@aust.edu.cn

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/dqwlxjz/CN/10.6038/pg2026JJ0228

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT