简介这套基于 MATLAB 的二维射线追踪与地震声波正演源码包面向地球物理、地震勘探专业的初学者与研究者用于模拟地震波在地层中的传播路径与接收信号。程序涵盖射线理论基础、几何扩散法、速度模型构建、源项与接收器设置、数值求解如龙格-库塔法及结果可视化等核心模块既适合理解正演原理也可作为二次开发的基础。资源共 29 个文件以 m 脚本为主另含 1 个 mat 速度模型数据压缩包大小 342KBm 文件分工明确包括主演示程序、射线追踪函数、速度模型加载与绘图工具等便于按模块调用。已有 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(k1) x(k) (dt/6)*(k1x 2*k2x 2*k3x k4x); z(k1) 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.mP波P波斯奈尔角相等常规纵波偏移traceray_ps.mP波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顺序检查所有事件道比人工抽查高效得多。本文还有配套的精品资源点击获取