简介：这套基于 MATLAB 的二维射线追踪与地震声波正演源码包，面向地球物理、地震勘探专业的初学者与研究者，用于模拟地震波在地层中的传播路径与接收信号。程序涵盖射线理论基础、几何扩散法、速度模型构建、源项与接收器设置、数值求解（如龙格-库塔法）及结果可视化等核心模块，既适合理解正演原理，也可作为二次开发的基础。资源共 29 个文件，以 m 脚本为主，另含 1 个 mat 速度模型数据，压缩包大小 342KB；m 文件分工明确，包括主演示程序、射线追踪函数、速度模型加载与绘图工具等，便于按模块调用。已有 1286 人学习使用。对于希望快速上手 MATLAB 地震波模拟、并掌握射线追踪算法实现细节的读者，这套源码能提供完整可运行的工程骨架、典型示例（如 Marmousi 模型测试）以及调试排错参考，是理论结合实践的有益资料。
1. 为什么地震声波正演还值得自己写一套MATLAB射线追踪
在拿到Marmousi模型或任意二维速度网格时，跑通一套射线追踪正演，并不只是为了画几张漂亮的波路径图。射线追踪真正解决的是地震勘探里最朴素的问题：从震源出发的波，走哪条路、花多长时间、以多大振幅到达检波器。这个问题搞不定，后续的速度分析、偏移成像、观测系统设计全都缺少一个可验证的前向算子。MATLAB里虽然有众多工具箱，但正演射线追踪很少直接给现成函数，因为射线路径对网格、界面和速度梯度的依赖太具体，必须按自己的模型手工实现。
这套二维射线追踪源程序的价值在于：它把地震声波正演从教科书公式变成了可运行的代码，文件里同时包含了打靶法（shoot-ray）和弯曲法（bend-ray）两类策略，还带上了Marmousi标准模型和基础偏移程序。适合正在做速度建模、初至拾取验证或观测系统覆盖次数分析的人作为骨架工程，也适合课程设计阶段需要跑通完整正演流程的学生。需要提醒的是，它不是一个点开就能出漂亮结果的成品应用，而是一套需要理解射线几何和数值传播逻辑的源程序，运行、调试、改造的过程才是真正有价值的部分。
2. 源程序拆解：从文件清单看清二维射线追踪的骨架
拿到压缩包后不要急着运行某个demo，先按文件命名把模块边界划清楚。这套程序的文件名非常直白，shootray系列负责从震源向外发射射线，traceray系列负责给定两点后反算路径，rayfan是扇形射线簇，rayvxz系列是速度模型上的走时计算，demoprep系列是数据准备，raymig系列是偏移。先读Contents.m，它相当于程序自带的地图，能帮你确认当前版本收录了哪些接口，也能避免漏掉某个mat数据文件导致运行中断。
2.1 主控流程与数据准备
raytrace_demo.m是入口，它做的事很典型：加载速度模型、设置震源和检波器、调用射线追踪函数、绘制结果。自己写demo时可以完全复用这个骨架。
% raytrace_demo.m 主控流程示意 clear; close all; load marmousi_mod.mat; % 载入Marmousi标准速度模型 vmodel = marmousi_mod; % 二维速度矩阵，单位 m/s src = [2000 20]; % 震源坐标：距离2000m，深度20m rec = [5200 800]; % 检波器坐标：距离5200m，深度800m % 调用射线追踪，返回路径和时间场 [tpath, ttime] = traceray(vmodel, src, rec, 'method', 'rk4'); % 可视化：速度模型上叠加射线路径 imagesc(vmodel'); colormap(jet); hold on; plot(tpath(:,1), tpath(:,2), 'w-', 'LineWidth', 1.5);这里marmousi_mod是程序自带的二维速度模型，行对应深度方向，列对应水平距离方向，坐标单位是米，速度单位是米/秒。src和rec是长度为2的行向量，分别代表水平坐标和深度坐标。traceray的method参数指定数值积分方式，默认可选'euler'或'rk4'，后者的精度更高，但计算时间约增加一倍。
demoprep.m和demoprep2.m这两份准备脚本经常被人跳过，但它们做的事情很关键：一是把测井或层速度数据插值成规则网格；二是预计算球面扩散补偿系数。sphdiv.m就是几何扩散补偿模块，它根据射线走过的路径长度计算振幅衰减因子，这个因子在后续用射线做正演波形合成时是必须乘上的。
2.2 速度模型的网格语义与梯度连续性
rayvelmod.m是速度模型接口，它负责根据传入的坐标从速度网格中插值出该点的速度值。写这类插值函数时有一个容易忽略的细节：MATLAB矩阵的索引从1开始，而物理坐标从0开始，必须先做坐标偏移再取网格索引。
% rayvelmod.m 核心逻辑示意 function [c, dcx, dcz] = rayvelmod(vmodel, dx, dz, x, z) nx = size(vmodel, 2); nz = size(vmodel, 1); ix = floor(x / dx) + 1; % 水平网格索引 iz = floor(z / dz) + 1; % 深度网格索引 ix = max(1, min(nx, ix)); % 边界截断 iz = max(1, min(nz, iz)); c = vmodel(iz, ix); % 线性插值速度 enddcx和dcz是速度沿水平方向和深度方向的梯度，射线追踪的核心方程里需要它们来计算慢度矢量的变化率。这个文件决定了程序是使用解析速度梯度还是数值差分梯度。drayvec.m和drayveclin.m的差别也在这里：前者假设速度场是连续变化的，用解析方式计算梯度；后者把速度层当作线性层，层内速度按深度线性变化，层间允许不连续。如果模型是Marmousi这种含断层和透镜体的复杂结构，我更倾向于用drayveclin，因为解析梯度在速度突变面上会给出误导性的射线偏转方向。
2.2.1 速度间断处的射线处理
eventraymod.m和eventraymig.m这两个文件专门处理速度界面上的事件。射线到达速度界面时，需要按斯奈尔定律计算透射和反射方向，这要求程序能检测射线是否跨过了速度不连续面。检测方法并不复杂，比较射线当前点的速度和前一步的速度，如果两者之差超过设定的相对阈值，就认为射线撞上了界面。阈值一般设为1%到2%，太低会把正常梯度误判为界面，太高会漏掉低速层的顶面。
2.3 “不成功”目录名不代表程序跑不起来
压缩包目录名是“二维射线追踪程序地震声波正演源程序不成功”，很多新手误以为这个程序有Bug。这里的“不成功”更合理的理解是：射线追踪算法本身存在大量路径不收敛的情况。例如射线在低速透镜体里反复折射、在临界角附近发生全反射、或者进入阴影区后根本无法到达期望检波器，这些都是物理上的“不成功路径”，而不是代码崩溃。理解了这一点，调试时就不会一看到空结果就怀疑程序，而是去检查震源出射角和速度模型的合理性。
提示：调试这类射线追踪程序，建议先跑
sphdiv.m确认几何扩散曲线是否平滑连续。如果振幅曲线出现跳变，说明射线路径在某个界面上发生了意料之外的折射，先修路径再谈波形。
3. 打靶法与扇形射线簇：地震波如何从震源一条条射出
射线追踪的核心是求解程函方程在射线路径上的投影，通俗说就是：给定初始位置和初始出射方向，一步步推出射线接下来的位置。shootray.m实现的就是标准的初值打靶法，它把射线追踪当作常微分方程初值问题来解，只要给定震源坐标和出射角，就能沿射线传播方向推进整条路径。
3.1 打靶法为什么能覆盖地下结构
打靶法的本质是对出射角进行离散采样。rayfan.m生成扇形射线簇，就是从震源位置发出一系列等角度间隔的射线，覆盖一个扇形角度范围。这些射线在非均匀介质中会自动弯曲，从而把地下的速度异常信息“画”出来。实现上，rayfan.m会反复调用shootray，每次只改变初始出射角。
% rayfan.m 扇形射线簇生成示意 angles = -60:5:60; % 出射角范围：-60度到60度，间隔5度 for i = 1:length(angles) theta0 = angles(i) * pi / 180; [xpath, zpath] = shootray(vmodel, src, theta0, 'rk4', maxstep); drawray(xpath, zpath); % 绘制单条射线 end出射角间隔的选择需要在计算效率和覆盖密度之间权衡。间隔5度适合观察整体结构，但会漏掉窄小的低速异常体；间隔1度能识别更细的结构，但计算量接近五分钟量级（取决于网格大小和射线长度）。shootrayvxz.m和shootrayvxz_g.m是打靶法的变体，前者在速度-深度坐标系下逐层推进，后者加入了梯度修正，适合处理速度垂向变化占主导的沉积层模型。
3.2 四阶Runge-Kutta推进射线路径
射线传播的常微分方程组基于慢度矢量p = grad(t)构成，二维情况下需要同时更新水平慢度px和垂直慢度pz，再据此更新射线位置。标准的四阶Runge-Kutta方法在MATLAB里实现非常干净。
% shootray.m 核心传播循环（RK4示意） dt = 0.001; % 时间步长，单位秒 for k = 1:nsteps [k1x, k1z] = rayvelmod(vmodel, x(k), z(k)); [k2x, k2z] = rayvelmod(vmodel, x(k)+0.5*dt*k1x, z(k)+0.5*dt*k1z); [k3x, k3z] = rayvelmod(vmodel, x(k)+0.5*dt*k2x, z(k)+0.5*dt*k2z); [k4x, k4z] = rayvelmod(vmodel, x(k)+dt*k3x, z(k)+dt*k3z); x(k+1) = x(k) + (dt/6)*(k1x + 2*k2x + 2*k3x + k4x); z(k+1) = z(k) + (dt/6)*(k1z + 2*k2z + 2*k3z + k4z); end这段代码的关键在于rayvelmod返回的并不是速度，而是慢度分量变化率。这里的x(k)和z(k)是当前射线位置，k1到k4是经典的四阶Runge-Kutta斜率，本质上是射线方向在四个子步上的偏转量。时间步长dt的选择需要匹配网格尺寸，一般让单步传播距离不超过网格间距的1/5，否则射线可能直接跨过一个速度为异常值的网格点。
3.3 射线簇与球面扩散的配合检验
rayfan_a.m在rayfan.m基础上增加了振幅属性，它会按几何扩散规律对每条射线分配振幅权重。这里有一个常见的认知误区：射线簇里相邻射线之间的距离增大，并不直接等同于球面扩散衰减，而是聚焦和散焦效应的叠加。真正的球面扩散补偿在sphdiv.m里，它按射线路径累计长度计算衰减因子。
| 函数名 | 输入要点 | 输出内容 | 适用场景 |
|---|---|---|---|
shootray.m | 震源坐标、出射角、数值方法 | 射线路径坐标序列 | 单条射线追踪 |
shootrayvxz.m | 震源坐标、出射角、速度-深度模型 | 逐层路径与走时 | 沉积层追踪 |
shootrayvxz_g.m | 含梯度修正的模型参数 | 带梯度偏转的路径 | 速度梯度显著的模型 |
rayfan.m | 震源坐标、角度范围、角间隔 | 扇形射线簇 | 覆盖范围观察 |
sphdiv.m | 射线路径坐标 | 几何扩散振幅因子 | 波形正演前的补偿 |
sphdiv.m的工作方式比较直观——对射线路径上每两个相邻采样点求距离增量，累加得到路径长度，再计算1 / sqrt(path_length)作为振幅衰减因子。控制台里看到振幅急剧下降的拐点时，不必急着调参数，先用plot检查该区域的射线是否发生了聚焦或散焦，这是几何传播的必然结果。
4. 从打靶到弯曲：traceray与PP/PS转换波的边值策略
打靶法的致命缺点是效率低：为了找到一条从震源到固定检波器的射线，需要打出几十条射线再插值逼近。traceray.m的思路完全不同，它直接以两端固定点（震源和检波器）作为边界条件，通过迭代修正路径形状来满足最小走时原理，这种两点射线追踪方法即弯曲法。
4.1 弯曲法的路径参数化与走时约束
弯曲法先把初始猜测路径离散成若干控制点，然后反复调整控制点位置，直到整条路径的走时达到局部极小。初始猜测通常使用直线连接震源和检波器，如果速度模型复杂，直线路径穿过的网格误差太大，迭代可能陷入局部极小。testray.m就是用来测试这种收敛情况的，它会输出每次迭代的走时残差，帮助你判断当前路径是否真的收敛到了全局最小走时。
% traceray.m 弯曲法迭代核心示意 path = linspace(src, rec, npoints); % 初始直线路径 TT = calcTT(vmodel, path); % 计算初始走时 for iter = 1:maxiter grad = calcGrad(vmodel, path); % 沿路径计算走时梯度 path = path - alpha * grad; % 负梯度方向修正路径 TT_new = calcTT(vmodel, path); if abs(TT_new - TT) < tol break; % 走时变化小于容差，收敛 end TT = TT_new; endalpha是步长因子，设置过大会导致路径在目标位置附近振荡，设置过小会需要几百上千次迭代才能收敛，经验取值范围是0.1到0.5之间。tol是走时容差，当地震数据采样率为1毫秒时，tol设为1e-6秒足够。梯度计算采用有限差分，对每个控制点分别在水平和深度方向加一个小扰动，观察走时变化量除以扰动量即可得到梯度近似值。
4.2 PP波与PS转换波的路径拆分
traceray_pp.m和traceray_ps.m处理的是两种不同的地震波类型。PP波是纵波从震源下行、经反射点再以纵波上行回到检波器；PS波则是下行纵波、反射后转换成横波上行。两者的射线路径在反射点处不对称，必须分别追踪下行段和上行段。
| 函数名 | 下行段波型 | 上行段波型 | 反射点条件 | 主要用途 |
|---|---|---|---|---|
traceray_pp.m | P波 | P波 | 斯奈尔角相等 | 常规纵波偏移 |
traceray_ps.m | P波 | S波 | 纵波入射角与横波反射角满足速度比关系 | 转换波成像 |
shootraytosurf.m | 任意角度 | 自由表面出射 | 出射点落在地表 | 地表接收记录 |
traceray_ps.m的难点在于反射点处需要同时满足两种波型的斯奈尔定律，而P波和S波速度比未知。常见的处理方式是先假设一个纵横波速度比（沉积岩中约为1.7到2.0），确定反射点初始位置，再通过迭代修正该比例。
4.2.1 事件记录与出射角校验
eventraymod.m负责把追踪成功的射线整理成事件记录，包含射线编号、出射角、到达时间、振幅衰减因子。这个模块还会做一个重要的合理性检查：检波器处的出射角必须在自由表面的接收范围内（通常是垂直方向±15度），如果超出这个范围，说明射线路径虽然数学上收敛，但物理上不可能被地表检波器记录到，程序应跳出该事件而不是直接保存。
4.3 初至走时与速度模型的双向验证
rayvxz_demo.m演示了一个很有用的验证思路——把走时计算结果和速度模型对照，确认走的路径确实经过目标地质体。具体操作是在速度模型上用imagesc显示速度分布，并叠加射线路径和等走时线。如果射线路径明显绕过了高速异常体，说明数值迭代过程中梯度计算方向有误，优先检查rayvelmod里速度梯度的正负号是否为水平方向深度方向的正确对应。
drayvec.m和drayveclin.m在这里容易被混淆。前者是向量化实现，一次计算多条射线在某个位置的慢度方向导数，适合批量处理扇形射线簇；后者是线性层解析解，用于速度在层内线性变化、层间不连续的情况。复杂模型用drayvec更稳妥，因为drayveclin对分层的依赖太强，遇到Marmousi这种断层错断的模型容易给出错误的层间归属。
5. 用射线走时做偏移和振幅保真校验
把正演射线追踪往前推一步，就是射线偏移。raymig.m和normraymig.m把观测走时沿射线路径反向投回到地下反射点，normray.m负责把幅值归一化到反射界面处，避免浅层强振幅压制深层弱反射。
在运行raymig.m之前，需要先用demoprep2.m把观测数据和速度模型整理到同一个网格结构中。一个快速验证流程是：加载marmousi_mod.mat，对每个炮集追踪射线簇，然后用sphdiv.m做几何扩散校正，最后把校正后的走时交给normraymig进行归一化偏移。
% 射线偏移执行示意 load marmousi_mod.mat; vmodel = marmousi_mod; [raypaths, traveltimes] = raytrace_demo(vmodel); image_profile = normraymig(vmodel, raypaths, traveltimes); imagesc(image_profile'); colorbar; title('射线偏移成像剖面');偏移结果的质量可以从两个角度检验：一是检查单条射线是否在反射点处满足反射定律，二是统计叠加剖面的信噪比。实际使用中，如果偏移剖面出现明显横向不连续条带，多数原因是射线覆盖不均匀，可以在rayfan里减小角度间隔重新生成射线簇。
clearrays.m值得多说一句，它的作用是从射线集合中剔除那些走时残差过大的病态射线。判断标准通常是：单条射线的走时残差超过该炮集中值的三倍以上，并且连续多次迭代都无法减小。这些射线大多穿过了速度间断面的临界角区域，强行保留会污染偏移剖面。
最后分享一个具体技巧：在traceray的迭代收敛之后，把残差走时和观测初至做差分，如果差值小于半个采样间隔（即1毫秒以下），说明该射线路径完全可靠。这个过程可以用testray.m顺序检查所有事件道，比人工抽查高效得多。
本文还有配套的精品资源，点击获取
