简介基于MATLAB实现的二维射线追踪程序是一套面向地震声波正演模拟的源代码包适用于地球物理、地震勘探、声波传播等方向的教学演示与科研复现。压缩包共30个文件包含28个M脚本、1个MAT数据文件和1个Markdown说明文档整体大小仅411KBM脚本除入口主程序外还集成射线追踪、射线扇、射线迁移、走时计算、非均匀介质模型正演等多个功能模块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(k1) xr(k) ds * sinth; % 在 ds 小步长内近似直线传播 zr(k1) zr(k) ds * costh; tr(k1) 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 就足够。下面的表总结了常用推进函数的差异。函数名推进方式适用速度模型主要风险drayvec.m数值递推任意二维连续模型速度突变处发散drayveclin.m线性梯度解析解速度平滑变化强间断处精度下降shootray.m角度扫描调用推进函数初至波覆盖计算临界角射线密集traceray.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 的岩层成立页岩或含气层会明显偏离。正演模式核心函数输入要求输出声波初至正演shootray.m drayvec.mvp 模型直达波、折射波走时声波反射正演traceray_pp.mvp 模型 反射界面反射波路径与走时扇形波前显示rayfan_a.mvp 模型 震源位置波前快照图转换波正演traceray_ps.mvp vs 模型PS 转换波路径检波面射线shootraytosurf.m速度模型 接收点数组地面记录道数据若要运行转换波需要先构造 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 和平滑核大小。如果误差集中在某一深度的强反射界面那要重点检查是否出现了射线路径跳跃尤其是临界角附近的射线。用这段插值走时场的算法你可以在几分钟内完成一次正演质量体检再带着可靠的走时数据去做偏移或反演。本文还有配套的精品资源点击获取 MATLAB地震射线追踪正演：从程函方程到Marmousi模型实践
简介基于MATLAB实现的二维射线追踪程序是一套面向地震声波正演模拟的源代码包适用于地球物理、地震勘探、声波传播等方向的教学演示与科研复现。压缩包共30个文件包含28个M脚本、1个MAT数据文件和1个Markdown说明文档整体大小仅411KBM脚本除入口主程序外还集成射线追踪、射线扇、射线迁移、走时计算、非均匀介质模型正演等多个功能模块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(k1) xr(k) ds * sinth; % 在 ds 小步长内近似直线传播 zr(k1) zr(k) ds * costh; tr(k1) 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 就足够。下面的表总结了常用推进函数的差异。函数名推进方式适用速度模型主要风险drayvec.m数值递推任意二维连续模型速度突变处发散drayveclin.m线性梯度解析解速度平滑变化强间断处精度下降shootray.m角度扫描调用推进函数初至波覆盖计算临界角射线密集traceray.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 的岩层成立页岩或含气层会明显偏离。正演模式核心函数输入要求输出声波初至正演shootray.m drayvec.mvp 模型直达波、折射波走时声波反射正演traceray_pp.mvp 模型 反射界面反射波路径与走时扇形波前显示rayfan_a.mvp 模型 震源位置波前快照图转换波正演traceray_ps.mvp vs 模型PS 转换波路径检波面射线shootraytosurf.m速度模型 接收点数组地面记录道数据若要运行转换波需要先构造 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 和平滑核大小。如果误差集中在某一深度的强反射界面那要重点检查是否出现了射线路径跳跃尤其是临界角附近的射线。用这段插值走时场的算法你可以在几分钟内完成一次正演质量体检再带着可靠的走时数据去做偏移或反演。本文还有配套的精品资源点击获取