地球物理射线追踪与偏移成像反演系统

资 源 简 介

本项目开发了一套完整的基于MATLAB的地球物理反演与成像软件框架，核心旨在通过数值模拟技术解析地下地质结构。功能涵盖三个主要方面：首先是高效的射线追踪系统，支持在复杂二维及三维速度模型（包括高斯平滑模型和块状模型）中进行波前延拓和射线路径计算，能够利用打靶法（Shooting Method）或弯曲法（Bending Method）精确获取地震波走时数据，为层析成像提供基础；其次是偏移成像处理，集成Kirchhoff叠前深度偏移（PSDM）算法，利用射线追踪计算的旅行时表将地震记录能量归位，重建地下界面的几何形态，显著提高复杂构造成像的分辨率；最后是反演优化模块，采用迭代线性化方法（如高斯-牛顿法）或全局优化算法，通过最小化观测走时与计算走时的残差目标函数来反演更新地下速度场分布。此外，程序包含交互式数据可视化工具，可动态展示射线传播过程、速度模型切片及最终的成像剖面，辅助研究人员评估反演质量。

详 情 说 明

项目：地球物理射线追踪与偏移成像反演系统

1. 项目简介

本项目实现了一套完整的轻量级地球物理数值模拟与成像软件框架。该系统完全基于MATLAB开发，集成在一个独立的脚本文件中，旨在通过数值计算手段解析地下地质结构。项目通过构建复杂的地下速度模型，模拟地震波传播（正演），并利用观测到的地震数据反推地下速度分布（层析成像）以及对地质界面进行成像（偏移处理）。

该代码主要用于教学演示及算法验证，直观展示了从地质建模、数据采集、速度反演到最终构造成像的地球物理全流程。

2. 功能特性

• 自定义地质建模：支持构建包含线性梯度背景和高斯型异常体的二维速度模型。
• 射线追踪系统：基于龙格-库塔（RK4）算法的试射法（Shooting Method）射线追踪。
• 合成地震记录生成：基于Born近似和雷克子波（Ricker Wavelet）卷积生成合成地震剖面。
• 走时层析成像：基于高斯-牛顿法的迭代线性化反演，支持稀疏矩阵构建与正则化求解。
• Kirchhoff叠前深度偏移：利用程函方程（Eikonal equation）计算旅行时场，实现深度域成像。
• 多视图可视化：交互式展示真实模型、射线路径、合成记录、反演过程及最终成像结果。

3. 系统要求

• MATLAB R2016b 或更高版本。
• 工具箱要求：Image Processing Toolbox (用于 imgaussfilt 函数)。
• 无需安装额外库，直接运行脚本即可。

4. 使用方法

1. 确保MATLAB当前路径包含 main.m 文件。
2. 直接运行 main 函数。
3. 程序将自动执行初始化、正演模拟、反演迭代和偏移成像，并在最后弹出综合结果窗口。

5. 代码实现细节与算法原理

本部分详细解析 main.m 代码中实际执行的逻辑流和关键算法。

5.1 模型构建与观测系统 (初始化阶段)

• 网格设置：定义了 $101 times 61$ 的正交网格，间距为 10m。
• 真实速度模型 ($V_{true}$)：由线性梯度背景 ($v = v_0 + kz$) 和一个位于中心的高速高斯异常体叠加而成，用于生成理论“观测”数据。
• 初始模型 ($V_{init}$)：仅包含线性梯度背景，作为反演的起始猜测值。
• 观测系统：在地表设置了 5 个炮点和 21 个检波器，覆盖主要成像区域。

5.2 正演模拟模块

该模块包含两部分核心计算：

1. 射线追踪 (Ray Tracing)：

* 采用 试射法 (Shooting Method)，从每个炮点以不同角度发射射线。 * 使用 RK4 (四阶龙格-库塔法) 对射线微分方程进行数值积分，计算射线路径。 * 数据捕获：程序判断射线末端是否回到地表并落在检波器附近。如果有，则记录走时（Travel Time）用于后续反演。

1. 合成地震记录 (Seismic Data Generation)：

* 并未求解全波波动方程，而是采用了高效的 Born近似 建模。 * 利用有限差分求解程函方程计算源点 ($T_s$) 和检波点 ($T_r$) 的旅行时场。 * 计算反射系数图 (Reflectivity = abs(del2(V_true)))，主要捕捉速度剧烈变化处。 * 在满足 $t = T_s + T_r$ 的时间点叠加雷克子波，生成含反射信息的地震记录。

5.3 速度模型反演 (层析成像)

程序执行 5 次迭代更新速度模型，核心逻辑如下：

• 正演重算：在当前反演模型中重新进行两点间射线追踪（或近似追踪）。
• 残差计算：计算观测走时与计算走时的差值 ($Delta t$)。
• Jacobian矩阵构建：

* 计算射线在每个网格单元内的穿行长度。 * 构建稀疏灵敏度矩阵 $G$，关联模型慢度变化与走时残差。

• 方程求解：

* 构建阻尼最小二乘目标函数：$(G^T G + lambda I) Delta m = G^T Delta d$。 * 求解慢度扰动量 ($Delta m$) 并更新模型。

• 物理约束与平滑：对更新后的速度模型进行数值范围限制和高斯平滑，保证迭代稳定性。

5.4 Kirchhoff叠前深度偏移 (PSDM)

利用反演得到的最终速度模型进行成像：

• 旅行时场计算：对每个炮点和检波点，再次调用程函方程求解器计算全场旅行时表。
• 成像条件：遍历所有地震道，将每个样点的值“涂抹”回地下满足 $t_{sample} = T_{source}(x,z) + T_{receiver}(x,z)$ 的等时面上。
• 滤波处理：对最终成像结果应用拉普拉斯算子 (del2)，去除低频成像噪声，提高分辨率。

5.5 结果可视化

代码最后生成一个包含六个子图的窗口：

1. True Model：显示真实速度场及部分采样射线路径。
2. Initial Model：显示反演前的初始梯度模型。
3. Seismogram：展示中间炮集（Common Shot Gather）的合成地震记录。
4. Inverted Model：显示经过迭代反演后的速度结构，展示对异常体的恢复情况。
5. Convergence：绘制RMS走时残差随迭代次数下降的曲线。
6. Migration Image：显示最终的Kirchhoff偏移剖面，重建地下反射界面的几何形态。

6. 注意事项

• 代码中的射线追踪使用了简化的捕获逻辑（判断射线末端与检波器距离），在稀疏观测系统下可能导致部分射线对丢失。
• 合成数据是通过Born近似生成的，主要用于演示偏移成像位置的准确性，而非精确的动力学特征。
• 代码片段尾部引用的辅助函数（如 V0_gradient_model, trace_ray_RK4, eikonal_solver 等）是 main 函数运行所必需的内部子函数。