二维弹性波方程交错网格数值模拟系统

资 源 简 介

本项目旨在利用MATLAB编程环境，针对一阶速度-应力形式的弹性波方程进行高精度的数值求解与仿真。核心算法采用交错网格有限差分法（Staggered Grid Finite Difference Scheme），将速度分量和应力分量定义在空间网格的不同节点上，以提高差分计算的稳定性和精度。具体的差分格式实现为时间方向二阶精度、空间方向四阶精度（2T4S），这种高阶格式能够有效压制数值频散，适用于复杂波场的精细模拟。系统具备完整的模型构建功能，支持用户定义均匀介质或层状介质模型的物理参数（纵波速度、横波速度、密度）。在边界处理方面，项目集成了吸收边界条件（如PML完全匹配层或Cerjan衰减边界），以在有限的计算区域内模拟无限介质，消除人工截断边界产生的虚假反射波。最终，系统通过可视化模块，按照设定的时间步长实时绘制并输出弹性波在介质中传播的波场快照，直观展示纵波（P波）和横波（S波）的波前扩散、反射与透射现象，为地震勘探基础理论研究及波场特征分析提供可靠的数值实验平台。

详 情 说 明

基于交错网格有限差分法的二维一阶弹性波方程数值模拟系统

项目简介

本项目是一个基于MATLAB开发的高精度地球物理数值模拟系统。它利用交错网格有限差分法（Staggered Grid Finite Difference, SGFD）求解一阶速度-应力形式的二维弹性波方程。该系统专为地震勘探基础研究和波场特征分析设计，能够模拟非均匀介质（如层状模型）中的弹性波传播过程，并提供实时的波场快照可视化。

核心算法采用了时间方向二阶、空间方向四阶（2T4S）的高阶差分格式，相比传统低阶格式，能更有效地压制数值频散，提高模拟精度。

功能特性

• 高阶数值精度：实现了空间四阶精度、时间二阶精度的交错网格差分格式，确保波场模拟的保真度。
• 非均匀介质建模：支持定义复杂的介质模型。当前实现包含双层介质模型，通过设置不同的纵波速度（Vp）、横波速度（Vs）和密度（Rho），模拟波在波阻抗界面处的反射与透射。
• 吸收边界条件：集成了Cerjan型阻尼衰减边界条件（Sponge Boundary），通过在计算区域四周设置具有特定衰减系数的吸收层，有效消除人工截断边界产生的虚假反射。
• 震源模拟：内置雷克子波（Ricker Wavelet）生成器，模拟中心激发的爆炸源（P波源）。
• 实时可视化：提供双视图动态显示，同时展示垂直速度分量（Vz）和水平速度分量（Vx）的传播过程，并具备自动波场能量归一化显示功能。

系统要求与使用方法

系统要求

• MATLAB R2016a 或更高版本（需支持基本的矩阵运算与绘图功能）。
• 足够的内存以存储全波场网格数据（默认400x400网格）。

使用方法

1. 确保运行环境已配置好MATLAB。
2. 直接运行主程序脚本。
3. 程序将自动初始化模型、计算差分系数、并在绘图窗口中实时演示弹性波传播动画。

---

核心代码逻辑详解

本项目的主程序（main）严格按照有限差分数值模拟的流程编写，具体实现逻辑如下：

1. 模型初始化与参数设置

• 网格定义：定义了400x400的计算网格，空间步长为10米。
• 介质建模：构建了双层地质模型。上层为低速层（Vp=2500m/s），下层为高速层（Vp=3500m/s），分界面位于深度索引200处。
• 物理参数转换：程序不仅存储速度和密度，还根据弹性动力学公式预计算了拉梅常数（Lambda, Mu）以及浮力参数（密度的倒数 B），以减少时间循环中的重复计算。
• 稳定性控制：根据Courant-Friedrichs-Lewy (CFL) 条件，自动计算满足稳定性的时间步长 dt。

2. 吸收边界构建 (Cerjan Damping)

• 代码实现了一个名为 DAMP 的衰减系数矩阵。
• 在网格四周定义了厚度为30个网格点的吸收层。
• 在吸收层内，衰减系数呈二次方规律（抛物线型）从内部向边界增加，计算出指数衰减因子。该矩阵在每个时间步直接作用于全波场，模拟能量流出边界的效果。

3. 差分系数计算

• 定义了标准的交错网格四阶差分系数：$c_1 = 9/8$, $c_2 = -1/24$。
• 为了优化性能，程序预先计算了系数除以网格步长（c1/dx, c2/dx 等）的值。

4. 时域有限差分循环 (Time Stepping)

主循环是通过时间步的迭代来推进波场，核心逻辑遵循速度-应力交错更新机制：

• 速度分量更新 (Vx, Vz)：

* 通过自定义的四阶差分函数计算应力分量（Txx, Tzz, Txz）的空间导数。 * 利用运动方程的一阶形式，结合浮力参数 B 和应力导数更新速度场。 * 应用 DAMP 矩阵进行边界吸收。

• 震源加载：

* 采用“软震源”方式，在网格中心位置将雷克子波的时间幅值叠加到正应力分量（Txx, Tzz）上，模拟各向同性的爆炸震源。

• 应力分量更新 (Txx, Tzz, Txz)：

* 计算更新后速度分量（Vx, Vz）的空间导数。 * 利用胡克定律（本构方程），结合拉梅常数和速度导数更新应力场。 * 再次应用 DAMP 矩阵进行边界吸收。

• 可视化渲染：

* 每隔10个时间步刷新一次图像。 * 使用 imagesc 绘制垂直和水平速度分量。 * 包含动态色标（caxis）调整逻辑，根据波场当前的最大幅值自动缩放显示范围，确保在波幅随距离衰减后依然能清晰观察到波前细节。

---

关键算法与实现细节分析

空间四阶差分算子 (Vectorized Implementation)

代码中并未采用低效的双重 for 循环来计算空间导数，而是实现了一个高效的辅助函数 diff_x_4th（以及对应的z方向逻辑）。

• 矢量化切片：该算法利用MATLAB的矩阵切片功能。例如，计算 $x$ 方向导数时，利用 U(idx+1, :)、U(idx, :)、U(idx+2, :) 和 U(idx-1, :) 的线性组合来实现差分公式：

$$ frac{partial U}{partial x} approx frac{1}{Delta x} left[ c_1 (U_{i+1} - U_i) + c_2 (U_{i+2} - U_{i-1}) right] $$ 此处 $U_i$ 与 $U_{i+1}$ 的差分对应交错网格中 $i+1/2$ 位置的导数。

• 边界处理：差分计算仅在网格内部区域（索引 2 到 Nx-2）进行，避免了数组索引越界。边界附近的数值由 DAMP 边界条件自然压制，因此未进行特殊的高阶边界修正。

交错网格机制 (Staggered Grid)

代码严格遵循交错网格定义：

• 速度节点：定义在整数网格点或半整数网格点上（依赖于具体的网格配置，通常Vx与Vz在空间上错开半个步长）。
• 应力节点：定义在速度节点之间的位置。
• 这种交错定义使得在计算速度导数时自然对应应力节点位置，反之亦然，从而大大减小了数值频散。

性能优化

• 预计算：所有材料参数（如浮力 B、模量 Lambda/Mu）均预先计算为全网格矩阵，计算循环中仅涉及矩阵点乘（Hadamard product）。
• 双缓冲绘图：开启 DoubleBuffer 并使用 drawnow limitrate，在保证动画流畅度的同时最小化绘图对计算速度的拖累。