本章把前面介绍的有限差分算子用于二维非均匀弹性介质中的地震波传播计算。一个完整的数值正演问题不仅要离散控制方程，还必须同时处理介质参数、变量的网格位置、时间推进、稳定性、自由表面、人工边界和震源。

二维地球介质模型

在二维直角坐标系 $(x,z)$ 中，通常令 $x$ 为水平方向，$z$ 为竖直方向。各物理量只随 $x,z,t$ 变化，并假定在垂直于计算平面的方向上不变。

各向同性弹性介质由密度和两个独立弹性参数描述：

$$\rho=\rho(x,z),\qquad \lambda=\lambda(x,z),\qquad \mu=\mu(x,z).$$

也可以使用 P 波速度 $\alpha$、S 波速度 $\beta$ 和密度表示：

$$\alpha=\sqrt{\frac{\lambda+2\mu}{\rho}}, \qquad \beta=\sqrt{\frac{\mu}{\rho}},$$

$$\mu=\rho\beta^2, \qquad \lambda=\rho\left(\alpha^2-2\beta^2\right).$$

二维弹性波可以分解为相互独立的两类：

• P-SV 波：位移或速度位于 $x-z$ 平面内，包括 P 波和垂直极化的 SV 波；
• SH 波：位移垂直于 $x-z$ 平面，只包含水平极化的 S 波。

本章主要讨论 P-SV 波。

二维弹性波动方程

动量守恒方程

设位移为

$$\mathbf u=(u_x,u_z)^T,$$

应力分量为 $\sigma_{xx},\sigma_{zz},\sigma_{xz}$，单位体积体力为 $f_x,f_z$。动量守恒方程为

$$\rho\frac{\partial^2u_x}{\partial t^2} =\frac{\partial\sigma_{xx}}{\partial x} +\frac{\partial\sigma_{xz}}{\partial z}+f_x,$$

$$\rho\frac{\partial^2u_z}{\partial t^2} =\frac{\partial\sigma_{xz}}{\partial x} +\frac{\partial\sigma_{zz}}{\partial z}+f_z.$$

本构关系

小变形条件下，应变为

$$\varepsilon_{xx}=\frac{\partial u_x}{\partial x}, \qquad \varepsilon_{zz}=\frac{\partial u_z}{\partial z},$$

$$\varepsilon_{xz}=\frac12 \left( \frac{\partial u_x}{\partial z} +\frac{\partial u_z}{\partial x} \right).$$

各向同性线弹性本构关系为

$$\sigma_{xx} =(\lambda+2\mu)\frac{\partial u_x}{\partial x} +\lambda\frac{\partial u_z}{\partial z},$$

$$\sigma_{zz} =\lambda\frac{\partial u_x}{\partial x} +(\lambda+2\mu)\frac{\partial u_z}{\partial z},$$

$$\sigma_{xz} =\mu\left( \frac{\partial u_x}{\partial z} +\frac{\partial u_z}{\partial x} \right).$$

将本构关系代入动量方程，可以得到只含位移的二阶方程；保留应力变量则得到更适合交错网格离散的一阶或混合方程组。

三种常用方程形式

传统网格位移方程

矢量形式为

$$\rho\frac{\partial^2\mathbf u}{\partial t^2} =(\lambda+\mu)\nabla(\nabla\cdot\mathbf u) +\mu\nabla^2\mathbf u+\mathbf f$$

在均匀介质中成立。非均匀介质中弹性参数位于空间导数内部，不能把 $\lambda$ 和 $\mu$ 简单移到导数外部。位移格式只需保存位移分量，但对介质界面和混合导数的离散较复杂。

位移-应力方程

保留位移和应力，可以写成

$$\rho\frac{\partial^2u_x}{\partial t^2} =\partial_x\sigma_{xx}+\partial_z\sigma_{xz}+f_x,$$

$$\rho\frac{\partial^2u_z}{\partial t^2} =\partial_x\sigma_{xz}+\partial_z\sigma_{zz}+f_z,$$

再用位移空间导数计算应力。该形式时间上仍为二阶。

速度-应力一阶方程

定义速度

$$v_x=\frac{\partial u_x}{\partial t}, \qquad v_z=\frac{\partial u_z}{\partial t}.$$

则 P-SV 波方程写为

$$\rho\frac{\partial v_x}{\partial t} =\frac{\partial\sigma_{xx}}{\partial x} +\frac{\partial\sigma_{xz}}{\partial z}+f_x,$$

$$\rho\frac{\partial v_z}{\partial t} =\frac{\partial\sigma_{xz}}{\partial x} +\frac{\partial\sigma_{zz}}{\partial z}+f_z,$$

$$\frac{\partial\sigma_{xx}}{\partial t} =(\lambda+2\mu)\frac{\partial v_x}{\partial x} +\lambda\frac{\partial v_z}{\partial z},$$

$$\frac{\partial\sigma_{zz}}{\partial t} =\lambda\frac{\partial v_x}{\partial x} +(\lambda+2\mu)\frac{\partial v_z}{\partial z},$$

$$\frac{\partial\sigma_{xz}}{\partial t} =\mu\left( \frac{\partial v_x}{\partial z} +\frac{\partial v_z}{\partial x} \right).$$

速度-应力形式所有方程对时间均为一阶，不含位移方程中的混合二阶导数，适合交错网格和显式时间推进。

网格离散与变量配置

设

$$x_i=x_0+i\Delta x, \qquad z_k=z_0+k\Delta z, \qquad t_n=n\Delta t.$$

传统网格

传统网格把 $u_x,u_z$ 及介质参数定义在同一节点。优点是数据结构直观，缺点是某些一阶导数、混合导数和介质界面参数需要额外平均，且可能产生网格解耦。

交错网格

标准二维速度-应力交错网格的一种典型配置为：

• $\sigma_{xx},\sigma_{zz},\lambda,mu$ 位于整数节点 $(i,k)$；
• $v_x$ 位于 $(i+\tfrac12,k)$；
• $v_z$ 位于 $(i,k+\tfrac12)$；
• $\sigma_{xz}$ 位于 $(i+\tfrac12,k+\tfrac12)$。

时间上也可以交错：应力位于整数时间层 $t_n$，速度位于半整数时间层 $t_{n+1/2}$。这样每个空间导数自然落在被更新变量所在的位置。

介质非均匀时，密度和剪切模量需要插值到速度或剪切应力节点。例如速度节点处常使用密度倒数的适当平均；$\sigma_{xz}$ 节点处的 $\mu$ 常采用调和平均，以更合理地处理介质界面。

二阶精度速度-应力交错网格格式

用

$$\delta_xf_{i+1/2,k} =\frac{f_{i+1,k}-f_{i,k}}{\Delta x},$$

$$\delta_zf_{i,k+1/2} =\frac{f_{i,k+1}-f_{i,k}}{\Delta z}$$

表示二阶精度中心差分。速度更新为

$$v_{x,i+1/2,k}^{n+1/2} =v_{x,i+1/2,k}^{n-1/2} +\frac{\Delta t}{\rho_{i+1/2,k}} \left[ \frac{\sigma_{xx,i+1,k}^n-\sigma_{xx,i,k}^n}{\Delta x} +\frac{\sigma_{xz,i+1/2,k+1/2}^n-\sigma_{xz,i+1/2,k-1/2}^n}{\Delta z} +f_{x,i+1/2,k}^n \right],$$

$$v_{z,i,k+1/2}^{n+1/2} =v_{z,i,k+1/2}^{n-1/2} +\frac{\Delta t}{\rho_{i,k+1/2}} \left[ \frac{\sigma_{xz,i+1/2,k+1/2}^n-\sigma_{xz,i-1/2,k+1/2}^n}{\Delta x} +\frac{\sigma_{zz,i,k+1}^n-\sigma_{zz,i,k}^n}{\Delta z} +f_{z,i,k+1/2}^n \right].$$

应力更新为

$$\sigma_{xx,i,k}^{n+1} =\sigma_{xx,i,k}^{n} +\Delta t\left[ (\lambda+2\mu)_{i,k} \frac{v_{x,i+1/2,k}^{n+1/2}-v_{x,i-1/2,k}^{n+1/2}}{\Delta x} +\lambda_{i,k} \frac{v_{z,i,k+1/2}^{n+1/2}-v_{z,i,k-1/2}^{n+1/2}}{\Delta z} \right],$$

$$\sigma_{zz,i,k}^{n+1} =\sigma_{zz,i,k}^{n} +\Delta t\left[ \lambda_{i,k} \frac{v_{x,i+1/2,k}^{n+1/2}-v_{x,i-1/2,k}^{n+1/2}}{\Delta x} +(\lambda+2\mu)_{i,k} \frac{v_{z,i,k+1/2}^{n+1/2}-v_{z,i,k-1/2}^{n+1/2}}{\Delta z} \right],$$

$$\sigma_{xz,i+1/2,k+1/2}^{n+1} =\sigma_{xz,i+1/2,k+1/2}^{n} +\Delta t\,\mu_{i+1/2,k+1/2} \left[ \frac{v_{x,i+1/2,k+1}^{n+1/2}-v_{x,i+1/2,k}^{n+1/2}}{\Delta z} +\frac{v_{z,i+1,k+1/2}^{n+1/2}-v_{z,i,k+1/2}^{n+1/2}}{\Delta x} \right].$$

上述格式按“由旧应力更新速度，再由新速度更新应力”的顺序循环推进。

高阶空间差分

为了减小数值频散，可以把二阶交错网格导数替换为高阶算子。四阶精度一阶导数为

$$\left.\frac{\partial f}{\partial x}\right|_{i+1/2} \approx\frac{1}{\Delta x} \left[ \frac98(f_{i+1}-f_i) -\frac1{24}(f_{i+2}-f_{i-1}) \right].$$

一般 $2M$ 阶算子写为

$$\left.\frac{\partial f}{\partial x}\right|_{i+1/2} \approx\frac1{\Delta x} \sum_{m=1}^{M}a_m \left(f_{i+m}-f_{i-m+1}\right).$$

提高空间阶数可以扩大有效波数范围，但会增加每次更新的运算量、边界模板宽度和并行通信宽度。时间上若仍采用二阶 leapfrog 格式，整体时间精度不会随空间阶数同步提高。

稳定性条件

对均匀介质和平面波进行 Von Neumann 分析，二阶时间、二阶空间标量波动格式满足近似 CFL 条件

$$c_{\max}^2\Delta t^2 \left(\frac1{\Delta x^2}+\frac1{\Delta z^2}\right) \leq1.$$

因此

$$\Delta t\leq \frac{1}{c_{\max} \sqrt{\Delta x^{-2}+\Delta z^{-2}}}.$$

当 $\Delta x=\Delta z=h$ 时，

$$\Delta t\leq\frac{h}{\sqrt2c_{\max}}.$$

弹性波计算中取最大 P 波速度 $c_{\max}=\alpha_{\max}$。高阶算子的稳定上限还取决于差分系数和方程形式，实际程序应采用相应理论条件并保留安全系数。

注意

CFL 条件只保证数值解不发散。为了得到可靠波形，时间步长通常还要满足时间采样精度要求，空间步长也要满足最短波长的采样要求。

数值频散与空间采样

对二维二阶差分标量波动方程，离散频散关系为

$$\sin^2\left(\frac{\omega\Delta t}{2}\right) =c^2\Delta t^2 \left[ \frac{\sin^2(k_x\Delta x/2)}{\Delta x^2} +\frac{\sin^2(k_z\Delta z/2)}{\Delta z^2} \right].$$

数值频率与解析关系 $\omega=c\sqrt{k_x^2+k_z^2}$ 不再完全一致，因而数值相速度依赖波长和传播方向。短波长分量误差最大，并会使波形随传播距离逐渐失真。

设震源有效最高频率为 $f_{\max}$，最小 S 波速度为 $\beta_{\min}$，模型中的最短弹性波长通常由 S 波控制：

$$\lambda_{\min}=\frac{\beta_{\min}}{f_{\max}}.$$

若要求每个最短波长至少有 $N_g$ 个网格点，则

$$h\leq\frac{\lambda_{\min}}{N_g}.$$

$N_g$ 取值取决于差分阶数和允许误差。二阶格式需要较多节点，高阶格式可以减少每波长节点数，但不能低于采样极限。

非均匀网格

浅部低速结构往往需要更小网格，而深部高速区域可用较大网格。非均匀网格可以降低总节点数，但普通等距中心差分公式不再适用。

在节点 $x_i$ 左右间距分别为 $h_{i-1}=x_i-x_{i-1}$ 和 $h_i=x_{i+1}-x_i$ 时，一阶和二阶导数权重需由非均匀节点上的 Taylor 展开重新求解。另一种做法是引入计算坐标 $\xi$，用坐标映射 $x=x(\xi)$ 把物理非均匀网格转为计算域等距网格，并根据

$$\frac{\partial}{\partial x} =\frac{1}{x_\xi}\frac{\partial}{\partial\xi}$$

修正微分算子。网格尺度突变会产生人为反射，因此网格间距应平滑变化。

边界条件

模型边界常见四类条件：

1. 自由表面条件；
2. 固定或刚性边界条件；
3. 周期边界条件；
4. 人工吸收边界条件。

自由表面

地表与空气接触时，表面牵引为零。若地表为水平面 $z=0$，法向量沿 $z$ 方向，则

$$\sigma_{zz}=0, \qquad \sigma_{xz}=0.$$

自由表面会产生反射以及 P-SV 转换。在交错网格中，可以令相应应力节点恰好位于自由表面并直接设置为零；对于缺失的表面外节点，可利用应力的奇偶延拓或虚拟节点构造差分值。边界处理必须与内部差分阶数和变量布局一致。

人工边界

有限计算区域的侧边界和底边界并不是真实反射界面。如果简单截断，外传波会发生强烈反射并污染计算结果。常用处理包括：

• 阻尼或海绵吸收层；
• 单程波近似吸收边界；
• 完全匹配层（PML）；
• 扩展模型并在外层逐渐衰减波场。

简单阻尼层可以在每一时间步把波场乘以位置相关衰减函数 $d(x,z)$：

$$q^{n+1}\leftarrow d(x,z)q^{n+1}, \qquad 0<d\leq1.$$

衰减应从物理区域向外平滑增强，突变的阻尼本身也会造成反射。

震源处理

体力源与力偶源

体力源直接加入动量方程的 $f_x$ 或 $f_z$。单力源具有方向性；一对方向相反、距离趋于零的力构成力偶。地震断层源通常用震源矩张量表示，其作用可等效为应力方程中的源项或空间导数形式的体力。

爆破源

各向同性爆破源在二维情况下可同时加入两个正应力分量：

$$\sigma_{xx}\leftarrow\sigma_{xx}+S(t)G(x,z),$$

$$\sigma_{zz}\leftarrow\sigma_{zz}+S(t)G(x,z).$$

它主要激发 P 波；双力偶源则同时产生具有辐射方向性的 P 波和 S 波。

震源时间函数

Ricker 子波常写为

$$S(t)= \left[1-2\pi^2f_0^2(t-t_0)^2\right] \exp\left[-\pi^2f_0^2(t-t_0)^2\right].$$

$f_0$ 是主频，$t_0$ 用于把主要能量移到计算开始之后。为了设计网格，常按震源谱中仍有显著能量的最高频率估计 $f_{\max}$，而不能只使用主频。

空间离散

理想点源包含无限高波数，直接放在单个节点上容易产生网格噪声。可用 Gaussian 核

$$G(x,z)=\exp\left[-\frac{(x-x_s)^2+(z-z_s)^2}{r_s^2}\right]$$

将震源平滑分配到若干节点。若震源不在网格点上，也可以根据插值权重分配给邻近节点，并保持总力或总矩守恒。

完整数值计算流程

1. 读入 $\alpha,\beta,\rho$，计算 $\lambda,\mu$；
2. 根据最高频率和最小波速确定空间步长；
3. 根据最大 P 波速度和差分阶数确定稳定时间步长；
4. 分配速度、应力和介质参数数组，并按交错位置进行参数平均；
5. 设置震源、接收器、自由表面和吸收层；
6. 初始化全部波场变量为零；
7. 由应力更新速度；
8. 加入速度方程中的体力源并处理速度边界；
9. 由速度更新应力；
10. 加入应力型震源并施加自由表面条件；
11. 记录接收点波形和指定时刻波场快照；
12. 滚动时间层，直至达到总计算时间。

结果检查与常见问题

• 到时检查：均匀模型中首波到时应接近距离除以理论波速；
• 对称性检查：对称模型和对称震源应产生相应对称波场；
• 能量检查：无吸收、无震源时能量不应无故增长；
• 网格收敛检查：减小 $h$ 和 $\Delta t$ 后，主要波形应趋于稳定；
• 边界检查：人工边界反射应明显弱于目标波场；
• 界面检查：参数平均错误会在界面产生非物理散射；
• 变量位置检查：交错网格数组下标错半个网格是最常见的程序错误之一；
• 震源检查：过窄的空间源或过高的频率会激发无法解析的短波。

版权所有

版权归属：Staaaaaaaaar

许可证：署名-非商业性-相同方式共享 4.0 国际 (CC-BY-NC-SA-4.0)

编辑此页

贡献者: Staaaaaaaaar

上一页切比雪夫伪谱方法

下一页有限单元方法