MATLAB地震射线追踪工具包：9函数实现建模、正演与速度反演

简介：本资源是一套面向地球物理、电子信息工程及应用数学等专业学习者的地震波传播仿真教学工具包，聚焦地震射线追踪算法实现与分层介质模型构建，适用于课程设计、毕业设计及科研入门阶段的原理验证与代码实践。压缩包共116个文件，以107个MATLAB函数（.m）为核心，涵盖主控脚本、射线路径求解、速度模型生成、时间场反演、横波可视化等关键模块；辅以5个Markdown说明文档和4个文本配置文件，清晰呈现算法逻辑、参数设置与运行流程；整体体积仅330KB，轻量易部署。目前已有224人学习下载，资源提供完整可运行源码、实测地层数据集及详细使用说明，特别包含射线管理器（tool_tx_manager）、时间组设定（fun_set_timegroup）、双参数最小二乘反演（fun_dmplstsqr2）等实用工具函数，便于读者理解正演建模思路、调试追踪过程并拓展自定义模型。

1. 地震波走时不是直线，但射线追踪能算出它真实路径——这套 MATLAB 工具包把复杂地层建模、初至波路径求解、速度场反演全打包进 9 个核心函数里

地质勘探中，地震波在地下并非匀速直线传播：遇到不同岩性界面会折射、反射，穿过低速带会弯曲，绕过高速体则加速。传统手工画射线或查表法误差大、不可复现；商业软件又黑箱难调试、参数不透明。这套基于 MATLAB 实现的地震射线追踪与地层模型仿真资源，用 9 个结构清晰的.m文件（main.m为入口，fun_calmod.m构建速度模型，fun_dmplstsqr2.m执行最小二乘反演），把从二维/三维分层介质建模、射线路径数值积分（采用四阶 Runge-Kutta 求解 Hamilton 方程）、初至时间正演计算，到基于走时残差的速度模型迭代更新，全部拆解为可读、可调、可验证的代码逻辑。它不依赖任何工具箱（仅需基础 MATLAB + Optimization Toolbox 中的fmincon或lsqnonlin，若无该工具箱可用fun_dmplstsqr.m的纯矩阵最小二乘替代），适合地球物理方向本科生做课程设计、研究生复现经典射线追踪算法、工程师快速验证某区块速度结构合理性。你不需要懂变分法推导，但得能看懂dx/dt = p_x / m,dp_x/dt = -∂H/∂x这类哈密顿系统离散化过程——这正是tool_r_in_editor.m和tool_tx_manager.m封装的底层逻辑。

2. 从分层速度模型构建到射线路径数值积分：fun_calmod.m与main.m协同完成正演全流程

2.1 分层介质建模：fun_calmod.m支持三种典型地层结构定义方式

fun_calmod.m是整个正演仿真的地基模块，它不直接生成网格，而是输出一个可被fun_txin_maker.m调用的结构体mod，其中包含z,vp,vs,rho四个关键字段（深度、纵波速度、横波速度、密度）。其输入支持三种常见工程场景：

• 分段常速层：传入向量z_layer = [0, 500, 1200, 2500]（单位：米）和对应vp_layer = [1800, 2400, 3100, 4200]，函数自动线性插值生成 100 层精细剖面；
• 平滑变化层：指定z_min,z_max,vp_min,vp_max,alpha（指数衰减系数），生成vp(z) = vp_min + (vp_max - vp_min) * (1 - exp(-alpha*(z-z_min)))类型的速度梯度；
• 用户自定义函数句柄：如@(z) 2000 + 1.5*z + 0.0002*z.^2，直接参与后续射线追踪的局部速度查表。

提示：fun_calmod.m默认输出mod.z为等间距深度向量（步长由dz参数控制，默认 10 米），但mod.vp是逐点计算值，不进行样条平滑——这是为避免射线积分时因插值引入虚假速度跳变，导致路径发散。若需更精细控制，可修改第 67 行z = z_min:dz:z_max;为非均匀采样。

2.2 射线路径求解：main.m调用fun_txin_maker.m启动四阶 Runge-Kutta 积分器

main.m是主控脚本，其核心流程是：加载模型 → 设置炮检对（source-receiver pairs）→ 调用fun_txin_maker.m计算每一对的射线路径与走时。fun_txin_maker.m内部实现的是基于哈密顿力学的射线追踪：

% fun_txin_maker.m 关键片段（简化示意） function [t_arr, x_ray, z_ray, p_x, p_z] = fun_txin_maker(mod, src, rec, opts) % src = [x_s, z_s], rec = [x_r, z_r] % 初始化射线参数：位置(x,z)与动量(p_x,p_z)，满足 |p|^2 = 1/v^2 x0 = src(1); z0 = src(2); p0 = ray_direction_initial(src, rec, mod); % 初值估计：直线方向+速度加权 % 四阶 Runge-Kutta 积分，步长 h 自适应（由 opts.h_init 控制） t = 0; x = x0; z = z0; p = p0; x_ray = x0; z_ray = z0; t_arr = t; p_x = p(1); p_z = p(2); while z < rec(2) && norm([x-rec(1), z-rec(2)]) > opts.tol_pos k1 = hamilton_deriv(x, z, p, mod); k2 = hamilton_deriv(x+h/2, z+h/2, p+h/2*k1, mod); k3 = hamilton_deriv(x+h/2, z+h/2, p+h/2*k2, mod); k4 = hamilton_deriv(x+h, z+h, p+h*k3, mod); p = p + h/6*(k1 + 2*k2 + 2*k3 + k4); % ... 更新 x,z,t 并检查是否越界 end end

其中hamilton_deriv计算哈密顿方程右端项：
dx/dt = p_x / (p_x^2 + p_z^2) * v^2，
dz/dt = p_z / (p_x^2 + p_z^2) * v^2，
dp_x/dt = -p_x * p_z * ∂v/∂z + p_x^2 * ∂v/∂x（忽略高阶项），
dp_z/dt = p_z^2 * ∂v/∂z - p_x * p_z * ∂v/∂x。

注意：∂v/∂x和∂v/∂z由mod结构体中的速度场通过gradient函数有限差分近似，不使用符号微分——这保证了任意复杂vp(z)函数均可接入，且计算稳定。若你的模型含陡峭界面（如断层），建议将opts.h_init设为 1~5 米（默认 10 米），否则 RK4 步长过大易跳过速度突变点。

2.3 走时与路径可视化：fun_vin_Swave_plot.m直接绘制横波射线与速度剖面叠加图

fun_vin_Swave_plot.m不仅画图，更承担结果验证功能。它接收fun_txin_maker.m输出的x_ray,z_ray,t_arr，并叠加mod.z,mod.vs曲线：

% 绘制横波射线（蓝色虚线）与速度剖面（红色实线） figure; hold on; plot(x_ray, z_ray, 'b--', 'LineWidth', 1.5); % 射线路径 plot(mod.vp, mod.z, 'r-', 'LineWidth', 2); % 纵波速度（常作为参考） xlabel('Horizontal Distance (m)'); ylabel('Depth (m)'); title(sprintf('S-wave Ray Tracing: %d points, Total Time %.3f s', ... length(x_ray), t_arr(end))); % 添加走时等时线（可选） if ~isempty(opts.time_levels) contour(X_grid, Z_grid, T_grid, opts.time_levels, 'k:', 'LineWidth', 0.8); end

该函数强制 y 轴反转（set(gca,'YDir','reverse')），符合地质剖面惯例；且自动缩放坐标轴使射线完全可见。若发现射线在某深度突然“折断”或发散，大概率是mod.vs在该处未定义（NaN）或fun_calmod.m生成的z向量未覆盖该深度——此时应检查mod.z范围是否 ≥rec(2)。

3. 从走时残差到速度模型更新：fun_dmplstsqr2.m实现带约束的最小二乘反演

3.1 反演问题建模：将速度模型参数化为分段线性节点

fun_dmplstsqr2.m解决的是经典地震层析反演问题：已知多组炮检对的观测走时t_obs，如何调整地下速度模型v(z)使正演走时t_calc与之匹配？它不采用网格单元参数化（易病态），而是将v(z)表示为N个控制节点(z_i, v_i)的分段线性函数。z_i固定（取mod.z的等间距子集，如每 200 米一个节点），待优化变量仅为v_i。目标函数为：

$$\min_{\mathbf{v}} \sum_{j=1}^{M} \left( t_{\text{calc}}^{(j)}(\mathbf{v}) - t_{\text{obs}}^{(j)} \right)^2 + \lambda \sum_{i=2}^{N-1} \left( v_i - \frac{v_{i-1}+v_{i+1}}{2} \right)^2$$

第二项是平滑约束（λ 由opts.lambda控制），抑制高频噪声。

3.2 雅可比矩阵解析计算：避免数值差分带来的精度损失

反演成败关键在于雅可比矩阵J_ij = ∂t_calc^{(j)}/∂v_i的准确性。fun_dmplstsqr2.m采用伴随状态法（Adjoint State Method）解析求导，而非对每个v_i做一次正演（计算量爆炸）。其核心是：对每条射线，沿积分路径反向求解伴随方程，得到∂t/∂v_i与射线在z_i附近停留时间成正比。代码中体现为：

% 在 fun_dmplstsqr2.m 内部，对第 j 条射线： % 1. 正向积分得路径 (x_k, z_k) 和时间 t_k % 2. 反向求解伴随变量 (q_x,k, q_z,k) 满足： % dq_x/ds = - (∂²H/∂x∂p_x)*q_x - (∂²H/∂x∂p_z)*q_z % dq_z/ds = - (∂²H/∂z∂p_x)*q_x - (∂²H/∂z∂p_z)*q_z % 3. 则 ∂t/∂v_i ≈ sum_{k where z_k near z_i} (q_z,k * ∂v/∂v_i) * Δs_k % 其中 ∂v/∂v_i 是分段线性基函数的导数（三角形函数）

提示：若你的观测走时含明显异常值（如某道信噪比极低），应在调用前用isoutlier(t_obs)标记并剔除，fun_dmplstsqr2.m不内置鲁棒加权机制。若必须保留，可将opts.weight设为向量，其元素为1./abs(t_obs - t_init)的归一化值。

3.3 约束优化执行：fmincon与lsqnonlin双模式切换

fun_dmplstsqr2.m默认调用fmincon（需 Optimization Toolbox），因其支持上下界约束v_min ≤ v_i ≤ v_max（如opts.v_bounds = [1500, 6000]）。若无该工具箱，则回退至lsqnonlin，此时需手动添加平滑项到残差向量：

% 伪代码：当 opts.solver == 'lsqnonlin' residual = [t_calc - t_obs; ... % 数据拟合残差 sqrt(lambda)*diff(v,2)]; % 平滑残差（二阶差分） [v_opt, resnorm] = lsqnonlin(@(v) my_residual(v), v0, [], [], opts);

两种模式下，v0均由fun_calmod.m的初始模型提供。反演后，新模型通过fun_set_timegroup.m更新到mod结构体中，供下一轮正演使用。

4. 多源协同管理与时间组设置：tool_tx_manager.m与fun_set_timegroup.m构建灵活实验框架

4.1tool_tx_manager.m：统一管理炮点、检波点、时间组的元数据容器

地震实验常涉及多炮多道（如 12 炮 × 48 道），手动维护坐标易错。tool_tx_manager.m将所有空间与时间信息封装为结构体tx_mgr：

tx_mgr = struct(... 'src', [x1,z1; x2,z2; ...], ... % Ns×2 矩阵，每行一个炮点 'rec', [x1,z1; x2,z2; ...], ... % Nr×2 矩阵，每行一个检波点 't_group', {g1,g2,...}, ... % 元胞数组，g1={1,5,12} 表示第1组含道1/5/12 't_obs', {t1,t2,...} ... % 对应每组的观测走时向量 );

fun_maintain_toolbox.m提供辅助函数：add_source(tx_mgr, [x,z])动态添加炮点；split_by_depth(tx_mgr, 1000)按检波点深度分组；export_to_segy(tx_mgr, 'data.sgy')导出标准 SEG-Y 格式（需额外 SEGY 工具箱）。关键设计是t_group字段——它允许将物理上连续的检波点划分为逻辑时间组，例如“浅层反射组”（0–800m）、“中层折射组”（800–2000m），每组独立反演，避免深部慢速带污染浅部分辨率。

4.2fun_set_timegroup.m：动态绑定模型与时间组，支持迭代反演中的模型热替换

在多轮反演中，需将新反演的速度模型mod_new快速应用到特定时间组的正演中。fun_set_timegroup.m完成此映射：

% 将 mod_new 应用于 tx_mgr 的第 2 组（索引为 2） tx_mgr = fun_set_timegroup(tx_mgr, mod_new, 2); % 内部操作： % 1. 提取该组对应的 src/rec 索引：idx_src = tx_mgr.src_idx{2}; idx_rec = tx_mgr.rec_idx{2}; % 2. 用 mod_new 重新计算这些炮检对的走时：t_calc = fun_txin_maker(mod_new, ...); % 3. 更新 tx_mgr.t_calc{2} = t_calc;

此函数确保tx_mgr始终持有当前最优模型的正演结果，无需重复调用main.m。若某次反演后tx_mgr.t_calc{2}与tx_mgr.t_obs{2}残差仍大，可单独对该组调用fun_dmplstsqr2.m，聚焦优化局部速度结构。

4.3 实战技巧：用tool_r_in_editor.m快速调试单条射线的初始动量

当某条射线积分失败（如t_arr为空或z_ray超出范围），最高效调试方式是可视化初始动量p0的敏感性。tool_r_in_editor.m提供交互式界面：

% 在命令行运行： tool_r_in_editor(mod, [100,0], [800,1500]); % 弹出窗口：滑块调节 p_x, p_z，实时显示射线路径与终点深度

窗口内，p_x滑块范围设为[-0.1, 0.1]（对应水平动量），p_z为[0.01, 0.05]（垂直动量）。观察当p_z过小时，射线无法到达接收深度；过大时则提前上翘。找到使z_ray(end)最接近1500的p0，将其设为fun_txin_maker.m的opts.p0_init，可大幅提升收敛成功率。此技巧比盲目减小步长更本质——它直击射线追踪的初值问题（Two-Point Boundary Value Problem）。

5. 验证与边界处理：用fun_dmplstsqr.m替代反演、用fun_vin_Swave_plot.m诊断速度跳变

5.1 无优化工具箱时的稳健替代：fun_dmplstsqr.m的纯矩阵最小二乘解法

当缺乏 Optimization Toolbox 时，fun_dmplstsqr.m提供降维方案：将速度模型固定为K个常速层（K通常取 5~10），每层速度v_k为未知数。正演走时t_j是v_k的显式函数（通过射线穿越各层厚度Δz_k计算：t_j = Σ Δz_k / v_k）。于是问题转化为线性系统A*v = t_obs，其中A_jk = Δz_k^{(j)}（第j条射线在第k层的厚度）。fun_dmplstsqr.m用v = A\t_obs求解，并添加 Tikhonov 正则化：

% A 是 M×K 矩阵（M 条射线，K 层），L 是 K×K 光滑算子（如二阶差分矩阵） v_opt = (A'*A + lambda*L'*L) \ (A'*t_obs);

此方法计算快、稳定性好，但牺牲了速度场的连续性表达能力。适用于教学演示或快速初值估计。

5.2 诊断速度不连续：fun_vin_Swave_plot.m的双纵轴模式识别界面反射

fun_vin_Swave_plot.m隐藏功能是启用双纵轴，同时显示速度剖面与射线路径曲率：

% 在调用时添加选项： fun_vin_Swave_plot(mod, x_ray, z_ray, 'curvature_on', true); % 自动计算射线曲率 κ = |d²z/dx²| / (1 + (dz/dx)²)^(3/2) % 并在右侧纵轴绘制 κ(z)，峰值位置即为强速度梯度区（如层界面）

若κ(z)在z=1200m处出现尖峰，而mod.vp在此处是平滑过渡，则说明模型未能刻画该界面——需在fun_calmod.m中增加一个控制节点或改用分段常速定义。反之，若κ平缓但mod.vp有陡坎，则可能是模型过度参数化，应增大fun_dmplstsqr2.m的opts.lambda。

5.3 关键边界条件处理：main.m中的opts.boundary控制射线终止行为

地下介质存在天然边界（地表z=0、莫霍面z=30km）。main.m通过opts.boundary字段定义射线如何响应：

启用'reflect'时，fun_txin_maker.m会在z=0处翻转p_z符号，并继续积分；启用'refract'则需计算临界角θ_c = asin(v_upper/v_lower)，并在满足条件时重置p_x,p_z。这些选项让同一套代码可模拟从浅层折射勘探到深部面波研究的多种物探模式。

本文还有配套的精品资源，点击获取