MATLAB地震射线追踪正演：从程函方程到Marmousi模型实践

简介：基于MATLAB实现的二维射线追踪程序，是一套面向地震声波正演模拟的源代码包，适用于地球物理、地震勘探、声波传播等方向的教学演示与科研复现。压缩包共30个文件，包含28个M脚本、1个MAT数据文件和1个Markdown说明文档，整体大小仅411KB；M脚本除入口主程序外，还集成射线追踪、射线扇、射线迁移、走时计算、非均匀介质模型正演等多个功能模块，MAT文件提供Marmousi模型数据，MD文档则写明使用说明与运行步骤，结构清晰，检索方便。程序以主函数为入口，在MATLAB 2020b及以上环境可直接运行，输出射线路径与走时结果，替换数据即可适配自定义场景；源码结构清晰，便于课程设计、算法复现与二次开发。目前已有135人学习/下载，代码经测试可稳定运行，配合说明文档可快速上手，适合初学者对照学习，也可供科研人员参考。

1. 为什么地震正演里射线追踪还没被淘汰

地震正演不是只有有限差分和伪谱法才能做。射线追踪在计算机上只解一个高频近似的偏微分方程——程函方程，跳过波动场的每一点压强变化，直接给出能量传播的路径和走时。正因为省掉了大量网格迭代，它在速度建模、走时反演、初至拾取和偏移预处理这些环节里，至今仍然是工业界默认的“第一张图”。

这套程序包来自 CSDN 用户 IT狂飙 上传的二维射线追踪资源，压缩包里有 raytrace_demo.m、traceray_pp.m、shootray.m、rayfan.m、rayvxz_wave.m 等二十多个 m 文件，另外带一份使用说明文档和 Marmousi 速度模型数据。它解决的问题非常具体：在已知二维速度模型的前提下，模拟地震波从炮点出发，经过地下介质传播到检波点的射线路径和波前快照。替换速度模型数据后，双击运行 raytrace_demo.m 就能得到结果，适合刚接触地震声波正演的学生，也适合需要快速验证射线参数的反演工程师。射线追踪比全波形正演快一到两个数量级，但它的适用边界也很明显——高频近似下的射线路径忽略了低频绕射和复杂焦散效应，理解这一点才能用好这套代码。

2. 从程函方程到 MATLAB 射线推进：drayvec 与 drayveclin 的差异

射线追踪的数学根基是程函方程：二维介质中走时 T(x,z) 满足 |∇T| = 1/v(x,z)，其中 v 是声波速度。这个方程是一个非线性偏微分方程，不能像波动方程那样直接时间递推，需要用特征线法将其转化为射线路径上的常微分方程组。常见做法是引入射线参数 p，在水平分层介质中，p = sinθ / v 沿射线保持不变，θ 是射线与垂直方向的夹角。这样，射线路径的推进就变成从当前点出发，按步长 ds 走一小段，再根据局部速度更新角度。s是路径长度，最终走时是路径上一系列 ds/v 的累加。

压缩包里的 drayvec.m 做的就是这件事，它把位置、慢度和走时打包成向量，用递推或 Runge-Kutta 积分推进。更实用的实现是 drayveclin.m，它假设速度在局部线性变化，可以解析地计算射线圆弧，稳定性比纯数值递推好。程序包里这两个文件都保留了，说明作者在“通用性”和“鲁棒性”之间做了区分。下面给出一段符合 drayvec.m 行为的最小逻辑示意代码，便于理解参数含义。

function [xr, zr, tr] = trace_snell_step(x0, z0, ang0, vfunc, ds, nstep) % trace_snell_step 基于斯奈尔定律的射线推进示意 % x0, z0 震源位置 % ang0 射线初始入射角，与铅垂方向夹角，单位度 % vfunc 速度函数句柄，v = vfunc(x, z) % ds 每步弧长，建议取网格间距的 1/5 ~ 1/2 % nstep 最大推进步数 xr = zeros(nstep, 1); zr = zeros(nstep, 1); tr = zeros(nstep, 1); xr(1) = x0; zr(1) = z0; tr(1) = 0; p_snell = sind(ang0) / vfunc(x0, z0); % 水平分层假设下的射线参数 for k = 1:nstep-1 vc = vfunc(xr(k), zr(k)); % 当前点速度 sinth = p_snell * vc; % 当前点射线倾角正弦 if abs(sinth) >= 1 break; % 达到临界角或全反射，终止 end costh = sqrt(1 - sinth^2); xr(k+1) = xr(k) + ds * sinth; % 在 ds 小步长内近似直线传播 zr(k+1) = zr(k) + ds * costh; tr(k+1) = tr(k) + ds / vc; % 走时增量 = 弧长 / 局部速度 end xr = xr(1:k); zr = zr(1:k); tr = tr(1:k); % 截断未到达的填充零 end

代码逻辑并不复杂：先根据初始角度和速度求出射线参数 p_snell，然后每走一步都用当前速度更新射线方向。因为步长 ds 足够小，局部可以视为匀速直线，方向变化完全由速度梯度控制。这比直接求解完整射线方程更直观，也更容易调参。参数中 ds 是关键，它不能大于速度模型网格间距的一半，否则路径会呈现折线抖动；但也不宜小于网格间距的十分之一，否则计算量增加数倍，走时精度却没有本质提升。ang0 决定了射线能覆盖到的深度范围，角度太大会让射线停在临界角之前，导致深层无覆盖。

drayveclin.m 则会在速度梯度方向上用一种解析公式推进，它避免了因速度突变导致的角度跳变，在 Marmousi 这类强变化模型中更稳定。两个文件在包内的角色可以这样区分：drayvec.m 负责通用介质，drayveclin.m 负责速度梯度平滑的介质。你在替换自己的速度模型时，如果模型来自测井插值，建议优先用 drayveclin.m；如果模型是层状均匀，drayvec.m 就足够。下面的表总结了常用推进函数的差异。

shootray.m 在包内的作用是扫描一组入射角，密集地调用射线推进函数，从而生成从震源出发的射线族。traceray.m 则是给定震源和检波点，反向迭代调整初始角度，直到射线恰好穿过检波点。这两个文件一个做“扫描”，一个做“定位”，是地震正演炮集生成的两类基本入口。实际跑 demo 时，raymarmousi_demo.m 会自动调度它们，不需要手动逐个调用，但改参数时你必须清楚自己改的是哪个环节。

3. 用 Marmousi 模型跑通射线追踪正演：参数设置与单位陷阱

在 MATLAB 中打开这个资源包，最稳妥的路径是把所有 m 文件和 marmousi_mod.mat 放在同一个当前目录，然后运行 raytrace_demo.m。如果你用的是 MATLAB R2023b 或者更新的 R2026a，脚本通常不需要改动；但 R2023b 对 imagesc、plot 的默认坐标方向处理和旧版本略有不同，如果绘图时地层看起来上下颠倒，在绘图语句后补一句 set(gca,'YDir','reverse') 即可。接下来需要手动检查的是速度模型的单位，这一步比任何代码都容易出错。

Marmousi 模型在公开资料里有不同的存储单位。压缩包里的 marmousi_mod.mat 如果数值范围在 1500~5500 左右，单位是 m/s；如果数值范围在 1.5~5.5 左右，则是 km/s。不统一单位时，走时计算会差三个数量级，画图时射线路径看着正常，但等时线完全贴在一起或散到几十公里外。建议在加载后立即输出 min 和 max 确认，再做一次单位转换逻辑。

% 加载模型并做基础预处理 load('marmousi_mod.mat'); vp = marmousi_mod; % 或者: vp = marmousi, 取决于 mat 文件内的变量名 fprintf('速度范围: %.2f ~ %.2f\n', min(vp(:)), max(vp(:))); % 如果范围在 km/s，统一转为 m/s if max(vp(:)) < 10 vp = vp * 1000; end % 设定网格参数 dx = 10; dz = 10; % 每网格代表 10m [nz, nx] = size(vp); x = (0:nx-1) * dx; z = (0:nz-1) * dz; % 炮点 sx = max(x)/2; sz = z(3); % 深度约 20~30m % 角度扫描范围，避开垂直向下的大角度冗余 angles = linspace(-70, 70, 201); % 射线步长取网格间距的 1/3 ds = 0.3 * min(dx, dz); nstep = ceil(sqrt(max(x)^2 + max(z)^2) / ds) * 2;

这里把角度范围限制在 -70° 到 70°，有两个原因：第一，接近 90° 的出射角在浅层会大量折射到水平方向，导致射线长时间停留在低速带，覆盖图里只是一堆密集横线；第二，偏移距过大的射线实际到达时间太晚，对反射波成像贡献有限。201 条射线对于一张 200×300 的网格模型足够，若模型更复杂，需要提高到 301 或 401。ds 的计算综合考虑了网格间距和模型尺寸，后面再乘以 2 是为了给深层射线留足推进步数，避免路径在深部被截断。

运行 demoprep.m 之后，程序会把速度模型转为射线追踪需要的网格梯度场。这一步容易出现“Matrix dimensions must agree”的报错，多数是因为 marmousi_mod.mat 里的矩阵维度排列为 [nx, nz]，而代码默认是 [nz, nx]。处理方式很简单，在读取后判断行列关系，用 vp = vp' 转置即可。另外，模型里如果有 NaN 或者零速度点，射线推进会直接失败。常见做法是用 rayvelmod.m 对 vp 做一次高斯平滑，同时把异常值替换为邻域中值。

% 剔除异常速度点，防止射线计算发散 vp(vp <= 0 | isnan(vp)) = min(vp(vp > 0)); % 高斯平滑，核宽 5 vp_smooth = imgaussfilt(vp, 5);

运行完成后图形窗口会显示从炮点向两侧扇形散开的射线路径，底图是速度色标。射线层接近水平时是明显的弧形，因为 Marmousi 模型中有大量低速楔形体，射线沿低速通道会发生弯曲。如果你发现射线在某一深度突然全部密集弯向水平，那不是 bug，而是临界角效应——速度增加导致 sinθ 达到 1，发生了全反射。此时应该缩小角度扫描范围，把最大角度减小到 60° 左右，或者增大射线衰减权重。如果你需要的是直达波初至，可以在 demoprep 中关闭反射追踪开关，只保留 shootray.m 的输出。

4. 扇形射线波前与 P-S 转换波：rayfan 和 traceray_ps 的使用边界

射线追踪不只是画几条线。rayfan.m 和 rayfan_a.m 的作用是从一个震源点发出扇形射线束，同时计算每条射线的走时，把所有射线末端等走时点连起来，就得到波前。这种波前快照适合用来观察声波在非均匀介质中的传播形态，也是验证速度模型是否合理的最直观方式。rayfan.m 与 rayfan_a.m 的区别在于后者支持将波前结果直接以谱图形式输出，更方便叠加在速度场上。实际调用时，目标函数参数中需要指定起始角度、结束角度和射线数，起始角与结束角之间的间隔越小，波前越光滑。

% 生成扇形射线，并画波前 figure; imagesc(x, z, vp_smooth'); colormap(flipud(gray)); hold on; rayfan_a(vp_smooth, sx, sz, -80, 80, 161); set(gca, 'YDir', 'reverse'); xlabel('Distance (m)'); ylabel('Depth (m)'); title('Ray fan wavefront');

这段代码里 161 是射线数，角度间隔为 1°，这个密度下波前连线基本是连续曲线。如果你的数据是弹性介质而不是纯声学介质，那就不能只靠声波走时，还需要考虑转换波。程序包里的 traceray_ps.m 专门处理 P 波入射、S 波反射/透射的情况。它需要的输入不仅是 P 波速度模型，还有 S 波速度模型。在纯声波正演中，S 波速度可以按下图中常用的经验关系近似：vs = vp / 1.73，但要注意这个关系只对泊松比约 0.25 的岩层成立，页岩或含气层会明显偏离。

若要运行转换波，需要先构造 vs 模型并存储在单独变量中，例如 vs_model = vp_smooth / 1.73，然后在 demoprep2.m 中把 vs_model 作为第二个速度参数传入。traceray_ps.m 在追踪过程中会先走 P 波段，到达反射界面后分裂出 S 波段，再在检波点合成旅行时。要注意的是，当入射角大于临界角时，转换波会变成首波而非反射波，此时程序通常返回 NaN 走时，需要将其从结果中置零或剔除。处理方式如下：

% 剔除转换波中的 NaN 走时，避免绘制时断裂 travel_time_ps = travel_time_ps(:); travel_time_ps(~isfinite(travel_time_ps)) = NaN; % 用邻域插值填充可以挖的缺口（如果缺口不超过3个采样点） travel_time_ps = fillmissing(travel_time_ps, 'linear', 'MaxGap', 3);

重点提醒：射线追踪的“射线”本身不携带振幅，程序包里的 sphdiv.m 用于球面扩散补偿，normray.m 用于射线归一化。这两个文件的作用是给走时记录配一个相对振幅，否则你只能看到同相轴的形态，看不到能量衰减规律。实际输出到图形上的颜色条不是真振幅，它的单位是相对值，不宜与地震记录振幅混用。

5. 射线覆盖质量检查：从走时插值到视速度验证

拿到射线追踪结果后，第一件事不是看路径图，而是检查覆盖密度。射线在高速区会发散，在低速区会聚焦，覆盖密度越不均匀，后续偏移成像的“脚印”越重。最简单的检查办法是把每条射线路径累积到一个二维计数矩阵中，然后以对数色标显示。射线覆盖密度图的用法是：高密度区域是能量聚焦区，该处走时可信；低密度区域是照射盲区，偏移后会出现明显的噪声条带。

countMap = zeros(nz, nx); % rayPaths 为射线结构体数组，包含 xr, zr 字段 for k = 1:length(rayPaths) xi = round(rayPaths(k).x / dx) + 1; zi = round(rayPaths(k).z / dz) + 1; valid = xi >= 1 & xi <= nx & zi >= 1 & zi <= nz; idx = sub2ind([nz, nx], zi(valid), xi(valid)); countMap(idx) = countMap(idx) + 1; end imagesc(x, z, log10(countMap + 1)); colorbar; axis xy;

这段代码使用累加计数的方式统计射线覆盖次数。index 计算时先做边界检查能有效防止数组越界。覆盖密度图中的对数色标让稀疏区域也能看到纹理，避免被高速区的高密度值淹没。如果右侧盲区过大，建议加密该区域的射线角度而不是增大所有角度范围。

如果射线本身正常，但走时等时线画出来不够圆滑，多半是走时场没有转为规整网格。rayvxz_wave.m 输出的走时是稀疏散点，需要插值成网格才能进一步分析。使用 griddata 之后，可以用走时场的梯度反推视速度，这是验证正演质量的非常有效的一招——视速度应该与输入模型的速度一致。

% 使用散点走时插值成网格 [Xg, Zg] = meshgrid(x, z); Tgrid = griddata(ray_x, ray_z, ray_t, Xg, Zg, 'cubic'); Tgrid = fillmissing(Tgrid, 'linear'); % 插值漏洞填充 % 由走时梯度计算视速度 [~, Tx] = gradient(Tgrid, dx); [~, Tz] = gradient(Tgrid, dz); vp_est = 1 ./ sqrt(Tx.^2 + Tz.^2); % 与原始模型速度比较 ratio = abs(vp_est - vp_smooth) ./ vp_smooth; fprintf('相对误差小于5%%的网格占比: %.1f%%\n', mean(ratio(:) < 0.05) * 100);

注意 gradient(Tgrid, dx) 的输出有两个梯度分量，第一个是沿 z 方向，第二个才是沿 x 方向。这里的写法刻意把第一个返回参数忽略，只取 Tx。计算出的 vp_est 若在边界区域出现尖峰，是因为插值在边界处产生了过冲，处理方法是先对 Tgrid 做一次裁剪，去掉边界外 5 个像素再求梯度。当误差较大的网格占比超过 20%，就应该回来调整 ds 和平滑核大小。如果误差集中在某一深度的强反射界面，那要重点检查是否出现了射线路径跳跃，尤其是临界角附近的射线。用这段插值走时场的算法，你可以在几分钟内完成一次正演质量体检，再带着可靠的走时数据去做偏移或反演。

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