一维波动方程全波形反演与数值模拟系统

资 源 简 介

本项目旨在构建一个基于MATLAB的完整计算平台，用于一维标量波动方程的数值正演与参数反演。项目首先利用高阶有限差分法（Finite Difference Method）实现波动方程的正演模拟，能够精确计算波在非均匀介质中的传播过程并生成合成地震记录。核心反演模块基于全波形反演（FWI）理论或最小二乘优化框架，通过构建目标函数（观测数据与计算数据的L2范数误差），利用伴随状态法高效计算梯度，并采用共轭梯度法或L-BFGS等优化算法迭代更新介质波速参数，从而从边界观测数据中重构内部介质的速度分布。此外，项目包含丰富的数据可视化功能，能够动态展示波场传播快照、实时绘制反演过程中模型参数的收敛轨迹、显示目标函数残差下降曲线，并提供观测波形与反演波形的拟合对比图，以直观评估反演效果和精度。

详 情 说 明

一维波动方程参数反演与可视化模拟系统

项目简介

本项目是一个基于 MATLAB 开发的计算地球物理模拟平台，专注于一维标量波动方程的数值模拟与参数反演。系统利用高阶有限差分法（Finite Difference Method）进行波场正演，并基于全波形反演（Full Waveform Inversion, FWI）理论，采用伴随状态法（Adjoint State Method）从观测数据中重构地下介质的速度分布。

该代码结构清晰，集成了物理建模、数值计算、反演优化及动态可视化模块，旨在演示如何通过边界或地表观测到的地震波形数据，反推内部介质的物理属性。

功能特性

• 高精度正演模拟：实现了基于有限差分算法的一维波动方程数值解，支持Courant稳定性条件自动计算时间步长。
• 非均匀介质建模：支持构建包含背景速度和高斯型异常体的复杂速度模型，模拟波在非均匀介质中的传播、反射与透射。
• 全波形反演（FWI）内核：实现了基于L2范数目标函数的反演流程，利用伴随状态法高效计算梯度，能够反演介质波速。
• 梯度优化策略：内置梯度计算、时间反转（Time Reversal）、梯度平滑处理及源点伪影压制（Source Masking）算法。
• 动态可视化系统：在反演迭代过程中实时展示速度模型更新轨迹、目标函数收敛曲线、灵敏度核（梯度）分布以及波形拟合情况。

系统要求

• MATLAB R2016b 或更高版本
• 无需额外工具箱，使用标准数学与绘图库

使用方法

1. 确保所有脚本文件位于 MATLAB 当前路径下。
2. 直接运行主程序入口脚本。
3. 程序将自动初始化物理参数、生成真实模型与初始模型、进行正演模拟生成观测数据，随后进入反演循环。
4. 运行过程中会弹出一个动态监视窗口，实时显示反演进度；结束后将生成最终波场对比图。

技术实现细节

本项目核心逻辑封装在主程序中，具体实现流程分析如下：

1. 基础参数与模型构建

程序首先定义了模拟区域的几何尺寸（2000m）和模拟时长（1.5s），根据最大估计速度自动计算满足稳定性的时间步长。

• 真实模型：构建了一个背景速度为 2000m/s 的均匀介质，并在中心位置叠加了一个高斯型高速异常体（幅度+500m/s），用于模拟地下构造。
• 初始模型：作为反演的起始点，采用常速度模型（2000m/s），假设未知异常体存在。

2. 观测系统设置

• 震源：使用 10Hz 主频的雷克子波（Ricker Wavelet），设置在网格左侧（索引5），低频选择有助于缓解FWI中的周波跳跃问题。
• 检波器：在全计算区域内每隔5个网格点布设一个接收器，用于记录全波场数据。

3. 正演模拟流程

通过有限差分算子求解波动方程，计算波在介质中的传播。程序首先利用真实速度模型生成“观测数据”（Observed Data），这代表了野外实际采集到的地震记录，作为反演的目标。

4. 全波形反演（FWI）主循环

反演过程采用迭代优化策略，最大迭代次数设为50次，主要步骤包括：

• 残差计算：在当前速度模型下进行正演，计算预测数据与观测数据的差值，并计算L2范数作为目标函数（Misfit）。
• 伴随源构造：将数据残差在时间轴上进行翻转（Time Reversal），作为新的震源项加载到检波器位置。
• 伴随场模拟：求解伴随方程，实际上是利用时间反转后的残差进行一次“反传”模拟。
• 梯度计算：

* 计算正演波场的二阶时间导数（加速度场）。 * 将伴随场翻转回物理时间轴。 * 利用互相关成像条件（Zero-lag cross-correlation）计算梯度：对正演加速度场与伴随场在时间域进行积分。 * 应用系数 $2/v^3$ 得到关于速度的梯度。

• 梯度预处理：

* 平滑：对原始梯度应用滑动平均滤波器，消除高频数值噪声。 * 源端压制：应用空间Mask窗口，强制将震源附近的梯度置零，防止源点强能量导致的数值伪影。

• 模型更新：

* 采用最速下降法（Steepest Descent）方向。 * 使用归一化梯度与固定步长相结合的策略更新速度模型。 * 物理约束：强制将更新后的速度限制在 [1500, 4000] m/s 的合理物理范围内。

5. 可视化分析

程序在迭代过程中实时维护一个四分屏图形窗口：

• 左上：展示真实模型、初始模型、当前迭代模型及更新后模型的速度剖面对比。
• 右上：绘制目标函数随迭代次数下降的收敛曲线（对数坐标）。
• 左下：显示当前计算得到的灵敏度核（梯度）的空间分布。
• 右下：选取特定检波器位置，对比观测波形与当前反演波形的拟合程度。

6. 最终结果验证

迭代结束后，程序利用最终反演得到的速度模型重新进行一次正演，并生成 X-T（空间-时间）域的波场快照，直观展示波场传播特征与真实情况的吻合度。