地震波场数值模拟是研究地震波传播规律的主要手段, 也是基于波动方程进行地震反演和成像的基础。有限差分方法具有实现简单、应用范围广、内存需求小、易于并行的优点, 被广泛应用于理论地震学和勘探地震研究中。目前的有限差分算法可以根据变量在计算网格上的分布不同而分成同位网格和交错网格两种网格系统。交错网格布局又可以根据波场分量在计算网格上的相对位置关系被大致分为三种格式:标准交错网格^[1-2]、完全交错网格^[3-4]和旋转交错网格^[5]。交错网格有限差分方法通常采用中心差分的方式求解一阶速度-应力方程组, 可以避免同位网格简单中心差分格式的奇偶失联问题^[6], 是目前求解波动方程数值解的主流方法之一。

有限差分波动方程数值解存在不同类型的误差, 其中与地震波频率有关的色散/耗散误差可以通过使用高阶差分算子或者优化格式来压制^[7-10]。而包含内部速度间断面的复杂速度模型采用有限差分网格离散时, 由于网格与速度界面不重叠, 同一物理模型采用不同的介质离散方法, 因此会导致不同的网格离散模型, 进而引起正演结果的不同。ALTERMAN等^[11]采用虚拟网格点和显式实施内部边界条件的方式处理二阶位移方程中的水平间断面问题。VIRIEUX^{[1, 12]}发展了基于一阶速度-应力方程的交错网格有限差分方法, 通过网格点离散取值表征速度模型, 并且展示了这种方法在大泊松比情况下的稳定性。GRAVES^[13]给出了交错网格有限差分算法处理三维弹性介质参数的格点平均方法。MUIR等^[14]首次将等效介质参数理论^[15]应用到有限差分方法中。MOCZO等^[16]提出了体积分平均法, 用来处理交错网格有限差分算法中的介质参数间断面, 解决了物理界面在网格格点之间变化时的表征问题。MOCZO等^[17]、KRISTEK等^[18]将离散后界面附近的介质近似为正交各向异性介质来实施间断面的边界条件, 从而有效地压制了界面误差。LISITSA等^[19]通过平面波方程分析了交错网格有限差分方法处理界面的理论误差。VISHNEVSKY等^[20]对标准交错网格、完全交错网格和旋转交错网格在层状界面下的误差收敛情况进行了总结。但是, 上述这些研究都没有定量评估和给出在任意复杂模型中准确计算(满足给定误差要求)反射/透射波的网格大小要求, 尤其是没有考虑准确计算P-S转换波所要求的网格大小, 所以目前进行有限差分地震模拟通常仅考虑格式本身的色散和耗散误差对网格大小的限制, 而没有考虑介质离散方法不同导致的与介质界面有关的反射/透射波计算准确性的影响。

本文首先介绍了交错网格有限差分方法原理, 然后总结了目前有限差分地震模拟中所使用的介质离散方法, 最后定量分析了不同介质离散方法对不同地震波场模拟准确性的影响。

1 交错网格有限差分方法

在二维笛卡尔坐标系下, 波动方程的速度-应力形式表示为:

$ \rho \frac{{\partial {v_x}}}{{\partial t}} = \frac{{\partial {\tau _{xx}}}}{{\partial x}} + \frac{{\partial {\tau _{xz}}}}{{\partial z}} $

(1a)

$ \rho \frac{{\partial {v_z}}}{{\partial t}} = \frac{{\partial {\tau _{zz}}}}{{\partial z}} + \frac{{\partial {\tau _{xz}}}}{{\partial x}} $

(1b)

$ \frac{{\partial {\tau _{xx}}}}{{\partial t}} = {C_{11}}\frac{{\partial {v_x}}}{{\partial x}} + {C_{12}}\frac{{\partial {v_z}}}{{\partial z}} $

(2a)

$ \frac{{\partial {\tau _{zz}}}}{{\partial t}} = {C_{21}}\frac{{\partial {v_z}}}{{\partial z}} + {C_{22}}\frac{{\partial {v_x}}}{{\partial x}} $

(2b)

$ \frac{{\partial {\tau _{xz}}}}{{\partial t}} = {C_{33}}\left( {\frac{{\partial {v_z}}}{{\partial x}} + \frac{{\partial {v_x}}}{{\partial z}}} \right) $

(2c)

式中:v_x与v_z代表速度分量; τ_xx, τ_zz与τ_xz代表应力分量; ρ表示介质密度; t表示时间; C_ij(i=1, 2, 3;j=1, 2, 3)为弹性刚度矩阵, 各向同性介质情况下, C_ij由两个独立的分量构成, 即拉梅常数λ与μ。

为了分析问题方便, 假定计算网格是矩形, 显式M(M=2, 4, 6, 8, …)阶交错网格差分算子D_x定义为:

$ \begin{array}{l} \frac{{\partial f}}{{\partial x}} = {D_x}\left[ f \right] = \frac{1}{h}\sum\limits_{m = 1}^M {{\alpha _m}\left[ {f\left( {x + mh - 0.5h} \right) - } \right.} \\ \;\;\;\;\;\;\;\;\left. {f\left( {x - mh + 0.5h} \right)} \right] \end{array} $

(3)

式中:h表示网格步长, α_m表示差分系数。则交错网格格式如下:

$ {{\hat \rho }_{i + 1/2,j}}{D_t}\left[ {{v_x}} \right]_{i + 1/2,j}^{n - 1/2} = {D_x}\left[ {{\tau _{xx}}} \right]_{i + 1/2,j}^{n - 1/2} + {D_z}\left[ {{\tau _{xz}}} \right]_{i + 1/2,j}^{n - 1/2} $

$ {{\hat \rho }_{i,j + 1/2}}{D_t}\left[ {{v_z}} \right]_{i,j + 1/2}^{n - 1/2} = {D_x}\left[ {{\tau _{xz}}} \right]_{i,j + 1/2}^{n - 1/2} + {D_z}\left[ {{\tau _{zz}}} \right]_{i,j + 1/2}^{n - 1/2} $

$ {D_t}\left( {{\tau _{xx}}} \right)_{i,j}^n = {{\hat C}_{11i,j}}{D_x}\left[ {{v_x}} \right]_{i,j}^n + {{\hat C}_{12i,j}}{D_z}\left[ {{v_z}} \right]_{i,j}^n $

$ {D_t}\left( {{\tau _{zz}}} \right)_{i,j}^n = {{\hat C}_{21i,j}}{D_x}\left[ {{v_x}} \right]_{i,j}^n + {{\hat C}_{22i,j}}{D_z}\left[ {{v_z}} \right]_{i,j}^n $

$ \begin{array}{*{20}{c}} {{D_t}\left[ {{\tau _{xz}}} \right]_{i + 1/2,j + 1/2}^n = {{\hat C}_{33i + 1/2,j + 1/2}}\left( {{D_x}\left[ {{v_z}} \right]_{i + 1/2,j + 1/2}^n + } \right.}\\ {\left. {{D_z}\left[ {{v_x}} \right]_{i + 1/2,j + 1/2}^n} \right)} \end{array} $

(4)

式中:n, i, j取正整数; $\hat \rho $ 与 ${\mathit{\boldsymbol{\hat C}}_{ij}}$ 分别代表计算网格格点上的密度与介质弹性参数。由于有限差分计算需要将实际物理模型离散为计算域中具体的网格点, 会出现实际界面位置与计算网格错位或不一致的情况, 因此需要对真实模型上的密度ρ与刚度系数矩阵C_ij进行特殊修正, 以逼近真实界面边界条件。

2 介质参数离散方法

对于介质参数不连续的内部界面, 其边界条件应该满足位移连续与法向牵引力连续的条件。目前较为常用的介质参数离散方法有直接取值法、格点平均法^[13]、体积分平均法^[16]与正交各向异性等效介质法^[17-18]。本文以一个倾斜界面模型为例, 说明不同介质离散方法的差异(图 1)。如图 1a所示, 不同介质的分界面倾斜穿过网格, 倾角为30°, 上层介质参数为v_P1=2 000m/s, v_S1=1 410m/s, ρ₁=1 160kg/m³; 下层介质参数为v_P2=4 000m/s, v_S2=2 310m/s, ρ₂=2 050kg/m³。图 1b是采用广义反射-透射系数方法^[21-23]按图 1a中的观测系统计算的接收波形(通过水平层计算后旋转得到), 震源为炸药震源。其中第一个事件(第一个检波器0~1.0s之间)是直达P波, 第二个事件(第一个检波器1.5~2.0s之间)是反射P波, 第三个事件是P-S转换波(第一个检波器2.3~2.6s之间)。以下介绍不同的介质离散方法、图 1a所示模型用不同介质离散方法得到的不同离散模型, 以及与图 1b对应的不同事件计算结果的差异(图 2, 图 3)。

2.1 直接取值法

直接取值法即直接将物理模型根据计算域划分的网格予以离散, 每个格点取其与物理模型相对应点的介质参数。界面倾斜穿过网格时, 该方法会形成阶梯状的介质参数间断界面, 见图 2a。根据惠更斯原理, 此时阶梯上的每一个点相当于一个点震源, 从而导致界面上产生人工的虚假波。图 3a展示了采用直接取值法离散介质所产生的散射波形, 可以看到在每道地震记录的P-S转换波后都伴随着明显的波形扭曲(第一个检波器2.6s之后)。因此直接取值法不能精确模拟界面的实际位置, 误差较为明显。

2.2 格点平均法

地震波穿过锯齿状弹性系数间断面时, 会在计算波形上产生明显的虚假散射波。为此GRAVES^[13]提出了格点平均法, 半网格点的介质参数通过整网格点的介质参数平均实现一定程度的平滑, 降低了虚假散射波的程度。该方法假定在计算域离散实际模型时, 整网格点{(x, y)|x∈Z且z∈Z}上的介质参数由物理模型直接给出, 半网格点{(x, y)|x∉Z或z∉Z}上的介质参数则需要通过其周围整网格点上的对应参数取平均获得。整网格点上的拉梅常数λ和μ可以由物理模型直接给出; 而剪应力对应位置的密度参数 $\hat \rho $ 与拉梅常数 $\hat \mu $ 由于分布在半网格点上, 故需要通过对其周围整网格点上的对应参数取平均来获得, 即:

$ \hat \rho = {\left( {\frac{2}{{\rho \left[ {i + 1,j} \right]}} + \frac{2}{{\rho \left[ {i,j + 1} \right]}}} \right)^{ - 1}} $

(5)

$ \begin{array}{l} \hat \mu = \left( {\frac{4}{{\mu \left[ {i,j} \right]}} + \frac{4}{{\mu \left[ {i + 1,j} \right]}} + \frac{4}{{\mu \left[ {i,j + 1} \right]}} + } \right.\\ \;\;\;\;\;\;\;{\left. {\frac{4}{{\mu \left[ {i + 1,j + 1} \right]}}} \right)^{ - 1}} \end{array} $

(6)

格点平均法对界面附近的介质做了简单的平均化处理, 避免了介质参数的剧烈间断(图 2b), 但这种方法本身无法对单位网格之内的界面变化进行细致的描述, 因此当界面位置与网格位置没有耦合的时候, 其模拟精度较差(图 3b)。

2.3 体积分平均法

当界面在单位网格间变化的时候, 简单的格点平均法无法准确反映界面位置。MOCZO等^[16]提出的体积分平均法可以更好地描述界面位置, 提高反射/透射波的模拟精度。在体积分平均法中, 某个点在计算域的介质参数不仅仅由这个点本身所在位置的介质参数决定, 还受到以它为中心的单位网格区域内的介质参数的影响。

设界面位于x=0处, 在一维情况下对波动方程在单位网格内做空间积分:

$ \int\limits_{ - h/2}^{h/2} {\rho \dot v{\rm{d}}x} = \int\limits_{ - h/2}^{h/2} {{\tau _{,x}}{\rm{d}}x} $

(7)

对(7)式应用均值定理, 得到:

$ \dot v\left( 0 \right)\int\limits_{ - h/2}^{h/2} {\rho {\rm{d}}x} = \frac{1}{2}\left[ {{\tau _{,x}}\left( {{0^ - }} \right) + {\tau _{,x}}\left( {{0^ + }} \right)} \right]\int\limits_{ - h/2}^{h/2} {{\rm{d}}x} $

(8)

则可以定义密度分量ρ的积分平均为:

$ {\rho ^{\rm{A}}} = \frac{1}{h}\int\limits_{ - h/2}^{h/2} {\rho \left( x \right){\rm{d}}x} $

(9)

ρ^A反映了这个被积分单元的平均响应。将其推广到二维情况下, 同理有:

$ {\rho ^{\text{A}}} = \frac{1}{{{h^2}}}\iint {\rho {\text{d}}x{\text{d}}z} $

(10)

$ {\lambda ^{\text{H}}} = {\left[ {\frac{1}{{{h^2}}}\iint {\frac{1}{\lambda }{\text{d}}x{\text{d}}z}} \right]^{ - 1}} $

(11)

$ {\mu ^{\text{H}}} = {\left[ {\frac{1}{{{h^2}}}\iint {\frac{1}{\mu }{\text{d}}x{\text{d}}z}} \right]^{ - 1}} $

(12)

如公式(10)至公式(12)所示, 对于网格上的某一点, 应用体积分平均法对其周围单位空间的介质参数进行积分然后取平均, 即可得到修正后的弹性介质参数。图 2c展示了体积分平均法处理后的介质参数, 与直接取值法和格点平均法相比, 这种方法能够反映界面在单位网格间的具体变化。同时, 体积分平均法对倾斜界面做了合适的处理, 与直接取值法相比有效压制了因倾斜界面介质离散不准确而产生的散射波(图 3c)。

2.4 正交各向异性等效介质法

层状介质研究表明, 由五个弹性系数描述的均匀横向各向同性介质是层状各向同性介质叠加物的长波等效^{[15, 24]}。其等效介质的弹性参数可以通过对层状介质拉梅常数的合理平均来近似获得。正交各向异性等效介质法基于这一等效理论, 通过将介质界面附近的弹性参数等效为一种正交各向异性介质弹性参数, 模拟界面在单位网格内的不同位置。

二维各向同性介质情况下:

$ \left[ {\begin{array}{*{20}{c}} {{\tau _{xx}}}\\ {{\tau _{zz}}}\\ {{\tau _{xz}}} \end{array}} \right] = \mathit{\boldsymbol{C}}\left[ {\begin{array}{*{20}{c}} {{\varepsilon _{xx}}}\\ {{\varepsilon _{zz}}}\\ {{\varepsilon _{xz}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\lambda + 2\mu }&\lambda &0\\ \lambda &{\lambda + 2\mu }&0\\ 0&0&{2\mu } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\varepsilon _{xx}}}\\ {{\varepsilon _{zz}}}\\ {{\varepsilon _{xz}}} \end{array}} \right] $

(13)

式中:ε_xx, ε_xz与ε_zz表示应变分量。将(13)式应力-应变关系应用于一个水平介质界面, τ_zz, τ_xz与ε_xx在界面两端连续, τ_xx, ε_xz与ε_zz在界面两端不连续, 用上标“+”表示在界面两端具有连续性的分量, 上标“-”表示在界面两端间断的分量, 则(13)式可以改写为:

(14)

其中,

$ {{\mathit{\boldsymbol{\hat C}}}_1} = \left[ {\begin{array}{*{20}{c}} {\lambda + 2\mu }&0\\ 0&{2\mu } \end{array}} \right]\\ {{\mathit{\boldsymbol{\hat C}}}_2} = \left[ {\begin{array}{*{20}{c}} \lambda \\ 0 \end{array}} \right]\\ {{\mathit{\boldsymbol{\hat C}}}_3} = \left[ {\lambda + 2\mu } \right] $

(15)

(14)式中弹性刚度矩阵被横线与竖线分为4个部分, 将每一个部分看作一个整体变量, 解出界面两端间断的分量τ^-与ε^-, 同时对等式两端取平均可得:

$ \begin{array}{l} \left\langle {{\tau ^ - }} \right\rangle = \left\langle {{{\mathit{\boldsymbol{\hat C}}}_3} - \mathit{\boldsymbol{\hat C}}_2^{\rm{T}}\mathit{\boldsymbol{\hat C}}_1^{ - 1}{{\mathit{\boldsymbol{\hat C}}}_2}} \right\rangle \varepsilon _{i,j}^ + + \left\langle {\mathit{\boldsymbol{\hat C}}_2^{\rm{T}}\mathit{\boldsymbol{\hat C}}_1^{ - 1}} \right\rangle \tau _{i,j}^ + \\ \left\langle {{\varepsilon ^ - }} \right\rangle = - \left\langle {\mathit{\boldsymbol{\hat C}}_1^{ - 1}{{\mathit{\boldsymbol{\hat C}}}_2}} \right\rangle \varepsilon _{i,j}^ + + \left\langle {\mathit{\boldsymbol{\hat C}}_1^{ - 1}} \right\rangle \tau _{i,j}^ + \end{array} $

(16)

式中, “〈〉”表示对单位网格求积分平均。(16)式反映了位于界面上的微小单元的应力-应变关系, 由于〈τ⁺〉=τ_{i, j}, 〈ε⁺〉=ε_{i, j}, 因此可改写为:

$ \left\langle \tau \right\rangle = \mathit{\boldsymbol{\hat C}}\left\langle \varepsilon \right\rangle $

(17)

方程(17)可以看作是根据界面边界条件对应力-应变关系进行了修正, 使非界面上的点与界面上的点同样满足该应力-应变关系。修正后的弹性系数刚度矩阵 $\mathit{\boldsymbol{\hat C}}$ 为:

$ \mathit{\boldsymbol{\hat C}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\Pi }_x}}&{{\lambda _{zx}}}&0\\ {{\lambda _{zx}}}&{{\mathit{\Pi }_z}}&0\\ 0&0&{2{\mu _{zx}}} \end{array}} \right] $

(18)

设M=λ+2μ, 则:

$ \begin{array}{l} {\mathit{\Pi }_x} = {\left\langle {{{\left\langle {M - \frac{{{\lambda ^2}}}{M}} \right\rangle }^z} + {{\left[ {{{\left\langle {\frac{\lambda }{M}} \right\rangle }^z}} \right]}^2}{{\left\langle M \right\rangle }^{{\rm{H}}z}}} \right\rangle ^{{\rm{H}}x}}\\ {\mathit{\Pi }_z} = {\left\langle {{{\left\langle {M - \frac{{{\lambda ^2}}}{M}} \right\rangle }^x} + {{\left[ {{{\left\langle {\frac{\lambda }{M}} \right\rangle }^z}} \right]}^2}{{\left\langle M \right\rangle }^{{\rm{H}}x}}} \right\rangle ^{{\rm{H}}z}}\\ {\lambda _{zx}} = {\left\langle M \right\rangle ^{{\rm{H}}zx}}{\left\langle {\frac{\lambda }{M}} \right\rangle ^{zx}}\\ {\mu _{zx}} = {\left\langle \mu \right\rangle ^{{\rm{H}}zx}} \end{array} $

(19)

式中:〈〉^H表示调和平均, 通过这样的平均化处理将界面附近的介质视为等效的正交各向异性介质, 将其弹性刚度矩阵代入应力-应变关系以模拟界面上的边界条件。与体积分平均法一样, 正交各向异性等效介质法能够准确地描述界面在单位网格之间的实际位置(图 2d), 从而能够获得更加准确的地震记录(图 3d)。

3 不同介质离散方法对不同事件模拟准确性的影响

本文通过改变计算域网格的每波长网格数(points per wavelength, PPW)来观察界面误差的变化情况, 分析了当每波长网格数逐渐增大时, 由界面计算产生的不同反射/透射波的相对误差。为了定量评估交错网格有限差分算法对界面计算的准确性, 首先定义相对误差为:

$ \sigma = \frac{{{{\left\| {{u^{{\rm{fd}}}} - {u^{{\rm{ref}}}}} \right\|}_2}}}{{{{\left\| {{u^{{\rm{ref}}}}} \right\|}_2}}} $

(20)

式中:‖u‖₂表示u的二范数, u^fd表示有限差分解, u^ref表示参考解。本文通过截取相应事件窗口的波形计算相对误差, 使用PPW来评估网格的分辨率。对于Ricker子波, 其PPW定义为:

$ {\rm{PPW}} \sim \frac{{{v_{\min }}}}{{2{f_{\rm{c}}}\Delta h}} $

(21)

式中:f_c为Ricker子波中心频率, Δh为空间网格步长, v_min表示P波或S波的最小速度。

在具体探索各种介质参数处理方法在不同波形下的相对误差PPW约束之前, 首先观察不同反射/透射角度对相对误差的影响。如图 4所示设置一个两层模型, 上层介质参数为v_P1=2 000m/s, v_S1=1 000m/s, ρ₁=1 200kg/m³, 下层介质参数为v_P2=4 000m/s, v_S2=2 000m/s, ρ₂=2 100kg/m³, 震源位于(1 000m, 1 500m)处。在距离界面上方400m与下方300m深度, 分别沿着图 4中的虚线布置两排检波器, 相对误差随不同反射/透射角的变化如图 5所示。由于某些位置反射/透射波与P-S转换波发生了耦合, 因此在展示相对误差随反射/透射角度的变化时, 使用整体波形计算其相对误差。

3.1 透射波和反射波相对误差随网格的变化

在图 5中相对误差较大的位置附近, 同时在时间上能够区分反射/透射波与P-S转换波的前提下, 将反射波接收点置于(1 930m, 2 060m)处, 透射波接收点置于(2 260m, 3 160m)处。通过保持震源、检波器与界面相对位置不变, 逐渐缩小网格空间步长以得到不同的PPW。为了使获得的结论更加一般化, 将实际物理界面位置置于两个网格之间(不与任何网格格点发生耦合), 使用主频为10Hz的Ricker子波作为震源子波, 分别对比不同反射/透射波在使用不同介质离散方法后的相对误差。为了避免色散误差的影响, 本文采用基于泰勒展开的10阶交错网格格式, 该格式满足格式的色散误差要求的PPW为5左右, 而下面测试的PPW最小是6, 因此所有测试都满足格式的色散误差要求。

图 6分别展示了反射波与透射波相对误差随PPW的变化, 可以看到不同介质离散方法对P波影响很大。直接取值法具有极大的不稳定性, 其对于波形的误差显著大于其它方法。就不同波形而言, 相同情况下透射波的准确性普遍高于反射波。图 7展示了反射波与透射波分别在相对误差约为10%与5%时的波形对比。

3.2 P-S转换波误差随网格的变化

图 8展示了P-S转换波相对误差随PPW的变化。由于界面上发生了P-S波的转换, 因此转换波(反射)的精度显著低于其它波形。由于直接取值法与格点平均法不能精确描述界面在网格间的具体位置, 因此它们在粗网格情况下误差较大, 而通过增大PPW提高模拟精度的计算代价高昂。P-S转换波在相对误差约为10%与5%时的波形对比见图 9。

3.3 多次反射波误差随网格的变化

为了观察多次反射波相对误差随网格的变化, 设置一个平层模型, 模型大小为n_x×n_z=2 500m×1 500m, 包括3层介质:顶层与底层纵波速度为4 000m/s, 横波速度为2 000m/s, 密度2 050kg/m³; 中间低速层纵波速度为2 000m/s, 横波速度为1 000m/s, 密度1 410kg/m³。震源(700m, 800m)与检波点(1 330m, 770m)均位于模型中间的低速层。图 10展示了多次反射波在不同PPW情况下, 使用不同介质离散方法得到的相对误差。图 11和图 12展示了误差约10%和5%时的波形对比(实际误差分别为9.61%与4.79%), 可以看到在相对误差约为5%时, 波形的拟合已经相对较好。

表 1与表 2分别展示了4种介质离散方法达到10%与5%的相对误差时, 不同类型反射/透射波各自的PPW要求。

4 结束语

本文针对介质的弹性系数处理对于交错网格波场模拟准确性造成的影响进行了分析, 给出了不同介质离散方式下各个反射/透射波相对误差的PPW约束参考, 得到以下结论：

1) 在使用交错网格有限差分进行正演计算时, 不仅要考虑格式本身的色散误差, 还需要考虑介质离散方式对模拟精度的影响, 特别是在复杂模型中, 即使在满足色散误差的情况下, 仍然存在较为严重的界面误差, 这类误差对整体波形有着不可忽视的影响。

2) 界面误差对PPW的要求普遍大于4阶格式的色散误差, 因此, 在选取交错网格的网格步长时, 为了保证误差在一定范围之内, 应该以横波最小速度对应的PPW为约束, 在可以接受的误差范围内首先确定网格的最大空间步长。

在实际交错网格计算中, 首先应该确定计算精度要求, 根据所采取的介质参数处理算法, 参考界面误差随PPW的变化关系, 以最小横波速度为基准, 确定网格的最大空间步长, 再根据柯朗-弗里德里希斯-列维(CFL)稳定性条件确定时间步长, 从而在确保计算准确性的前提下, 减少计算资源需求。