基于交错网格的有限差分地震波场正演模拟程序

资 源 简 介

本程序实现了基于交错网格有限差分（SGFD）方法的地震波场数值模拟系统。其核心功能是利用交错网格技术对声波或弹性波运动方程进行离散化求解，通过在不同的空间网格点上交错定义速度分量与应力分量，极大地提高了数值模拟的稳定性并有效抑制了常规网格中常见的数值色散现象。程序支持对复杂非均匀地质模型进行精确模拟，允许用户自定义地质结构的纵波速度、横波速度及密度分布。系统内置了高阶空间差分算子，能够灵活处理不同精度要求的波动方程计算。该程序是开展全波形反演（FWI）研究不可或缺的正演计算引擎，同时也广泛应用于地震学基础

详 情 说 明

基于交错网格有限差分（SGFD）的地震波场正演模拟系统

项目介绍

本项目是一款基于交错网格有限差分方法（Staggered-Grid Finite Difference, SGFD）开发的二维弹性波场数值模拟工具。该程序旨在通过数值手段模拟地震波在地下复杂介质中的传播过程，包含波的反射、折射及衍射等物理现象。系统采用了高阶空间差分算子与完全匹配层（PML）吸收边界条件，能够输出高精度的波场快照和地震记录，为全波形反演（FWI）、地震成像研究及油气勘探理论分析提供核心计算引擎。

功能特性

• 二维弹性波模拟：支持同时计算纵波（P）和横波（S）的运动过程。
• 高阶空间精度：内置四阶空间有限差分算子，有效抑制数值色散，提高计算效率。
• 交错网格技术：通过在不同的空间位置交错定义速度分量（vx, vz）与应力分量（txx, tzz, txz），保证了数值模拟的稳定性。
• PML吸收边界：实现了分裂格式的完全匹配层边界处理，有效吸收模型边缘的虚假反射。
• 复杂介质支持：允许用户自定义非均匀速度模型和密度模型，内置双层介质示例。
• 实时可视化分析：程序执行过程中可动态显示垂直速度分量波场快照和实时地震记录（Common Shot Gather）。
• 全方位结果总结：模拟结束后自动汇总模型速度结构、地震记录图、震源子波形态及单道波形数据。

系统要求

• 运行环境：MATLAB R2016b 或更高版本。
• 硬件要求：建议 8GB 以上内存以满足大规模网格计算需求。
• 依赖模块：无需第三方工具箱，基于 MATLAB 原生数学运算库实现。

实现逻辑与功能结构

本程序的实现逻辑严格遵循地震波动力学方程的离散化求解过程，具体步骤如下：

1. 参数初始化与模型构建 程序首先定义探测区域的物理尺寸（网格点数与空间步长）及模拟持续时间。随后构建纵波速度（Vp）、横波速度（Vs）和密度（rho）的二维矩阵，并通过逻辑索引定义了具有速度差异的双层地质模型。

2. 媒质参数转换 利用获取的速度和密度数据，预先计算拉梅常数（$lambda$ 和 $mu$）。这些常数是后续更新应力分量的核心物理参数。

3. 离散化算子定义 系统定义了四阶交错网格差分系数（c1=9/8, c2=-1/24）。其核心逻辑在于利用周围四个相邻点的值来估计中心点的导数，从而在较低的网格密度下获得较高的计算精度。

4. PML 边界层预计算 为了消除边界反射，程序在模型四周设置了厚度为 npml 的吸收层。通过计算二次方衰减廓线，生成 x 方向和 z 方向的衰减因子（ax, az）和补偿算子（bx, bz）。此部分逻辑决定了 PML 区域内波场的衰减速度和吸收效果。

5. 震源激发 采用 Ricker 子波作为震源，通过预设的中心频率计算波形。在时间迭代过程中，将震源振幅加载到应力分量（txx 与 tzz）上，模拟一个膨胀点源的激发效果。

6. 核心时间迭代循环 程序执行循环迭代，每个时间步内包含以下计算逻辑：

• 速度场更新：根据应力的空间梯度更新 vx 和 vz 分量。在更新过程中，使用了 PML 分裂波场技术（vx_x, vx_z 等），分别处理不同方向的导数项。
• 应力场更新：在更新后的速度场基础上，计算速度空间梯度，进而更新三个应力分量（txx, tzz, txz）。更新过程考虑了介质的非均匀性，对剪切模量 $mu$ 进行了交错点处的调和平均处理。
• 记录与显示：实时提取地表位置的 vz 分量并存入地震记录矩阵，同时每隔固定步数刷新可视化图像。

关键算法与实现细节分析

• 交错网格差分方案：不同于常规网格，交错网格将速度分量放置在网格线段中点，将应力分量放置在网格中心或顶点。这种空间布局天然适配一阶偏微分方程组的求解，避免了网格解耦问题。
• 分裂格式 PML (Split-field PML)：程序将每一个波场分量拆分为两个子项（如 vx 拆分为 vx_x 和 vx_z），分别对应 x 方向和 z 方向的偏导数贡献。通过对每个子项应用独立的指数衰减，实现了在有限厚度内耗散波场能量的目的。
• 高阶差分系数应用：在代码中，c1 和 c2 直接参与梯度的数值计算。四阶精度不仅提高了计算性能，还允许在大步长、高频率的情况下保持较低的数值相位误差。
• 介质平均技术：在处理交错点上的物理参数时（如 rho_half 和 mu_half），程序采用了算术平均或调和平均的思路，确保了在速度不连续界面处弹性波传播的物理准确性。

使用方法

1. 启动程序：在 MATLAB 命令行窗口直接运行 main 函数即可开启模拟。
2. 调整模型：用户可通过修改“参数设置”部分中的 Vp、Vs 矩阵来设计自己的地下构造（如添加褶皱、断层或透镜体）。
3. 配置震源：修改 f0 可调整主频，修改 src_x 和 src_z 可改变震源位置。
4. 观测点设置：通过调整 rec_z 和 rec_x 的定义，可以模拟地表接收、井中观测或 VSP 观测模式。
5. 输出分析：模拟运行结束后，程序会弹出一个综合报告窗口，显示速度模型图、最终地震记录、震源形态以及指定观测点的单道波形，便于用户进行时域分析。