简介一套用于电磁波走时层析成像的MATLAB反演程序面向地球物理、地质勘探及科研人员帮助解决地下速度模型构建与走时数据反演问题。压缩包共5个文件全部为.m脚本整体仅4KB代码简洁涵盖走时正演、层析反演迭代、模型更新等核心过程各脚本功能独立分别承担数据处理、正演模拟、反演计算与结果输出等任务便于阅读和二次开发。已有290人学习下载适合作为学习走时层析成像算法及MATLAB实现的入门样例。借助这套程序可以完成从正演模拟到反演计算的完整流程理解通过最小化观测走时与预测走时误差逼近真实速度结构的原理结合当前速度模型逐层修正、迭代收敛最终获得地下介质的速度分布为地质勘探、工程探测等实际应用提供可运行的工具脚本。1. 走时层析成像的MATLAB反演先看清这个程序解决什么问题做过工程物探或岩土检测的人多半遇到过这样的场景现场数据采回来一堆电磁波或弹性波走时甲方只看你最终能不能给出一张漂亮的内部速度分布图。市面上商业软件能出图但参数封在黑匣子里换一个观测系统、换一种介质分布结果就变得不可控。走时层析成像travel time tomography恰恰是用MATLAB就能从头实现的一类反演算法——通过初至走时反推介质内部的速度结构或衰减参数在堤坝隐患检测、混凝土结构评估、既有道路病害筛查里都有人用。这个标题里的diancibo显然指向电磁波类数据做法上与地震波走时层析几乎同构核心都是射线追踪加迭代求解。这套程序适合谁适合手头有走时数据且想自己控制反演全流程、能接受自己调参数踩坑的开发者不适合只想要一键出图的人。2. 走时层析反演的基本原理从射线追踪到线性方程组2.1 走时的物理含义与离散化方式电磁波或弹性波在介质中传播从发射点到接收点会经历一个时间开销这就是走时。不同介质中传播速度不同如果目标区域内部存在低速异常体波经过它时走时就会比正常路径更长。走时层析反演做的事情就是利用大量发射-接收对的走时数据反过来求解这个慢度场速度的倒数 s(x) 1/v(x)的空间分布。把连续介质离散成一个个网格单元这是数值反演的第一步。我在实际处理中通常会把目标区域剖分成规则矩形网格每个网格内认为慢度恒定这一步直接决定了未知数的个数。假设有 m 条射线剖分成 n 个网格每条射线在各个网格里穿过的路径长度是不同的于是走时和慢度之间构成一个线性关系每条射线的总走时等于它穿过每个网格的长度乘以该网格慢度之和。这样的处理带来的直接好处是反演问题被转化成一个线性方程组 d A·s其中 d 是观测走时向量长度是 ms 是待求慢度向量长度是 nA 是灵敏度矩阵也叫射线路径矩阵A(i,j) 表示第 i 条射线在第 j 个网格内的传播路径长度。当 m n 且射线覆盖充分时方程组是超定的可以用最小二乘思想来求解当 m n 时则欠定必须引入正则化约束才能得到稳定解。实际观测数据往往有噪声所以严格拟合每一条走时并非上策。我一般会在求解时写上一个阻尼项允许解在满足数据的同时不要过于跳变——这就是后面要讲到的阻尼最小二乘。搞清楚这个线性关系是理解整个程序的关键所在后续所有 MATLAB 代码都是围绕这个矩阵的构建与求解展开的。2.2 射线追踪与灵敏度矩阵的构建灵敏度矩阵的构建是走时层析里最费工夫的一环。常见做法有两个一是最短路径射线追踪如 Dijkstra 算法二是打靶法shooting或弯曲法bending。前者把介质切成细网格在所有网格点之间搜出一条真实意义上的最短传播路径后两者从震源出发按某种方向推进再逐步修正射线的出射角使射线终点恰好落在接收点位置。我在 MATLAB 中实现时最常用的是最短路径类的粗糙网络射线追踪先在剖分好的网格上叠加一层更细的辅助点相邻点之间的距离由局部慢度决定然后调用图搜索算法求最短路径。因为相邻点之间的距离其实就是走时的近似求到的最短路径就是正演的射线。这一步的好处是稳定基本不会出现射线发散或者打不到接收点的情况代价是要额外维护一个较大规模的稀疏图。构建灵敏度矩阵的代码可以参考下面这段写法核心是记录每条射线经过了哪些网格以及各自的路径长度。先把网格剖分参数准备好再对每一对发射-接收点调用射线追踪最后把结果填入稀疏矩阵 A 中。% 灵敏度矩阵构建射线追踪 路径长度积分 % nx, ny: 网格列数与行数; dx, dy: 网格尺寸 % rays: 结构体数组, 每个元素含 sr(cell) 坐标, rc(rx) 坐标, travel 观测走时 A spalloc(nrays, nx*ny, nrays * max_path_nodes); for k 1:nrays % 1. 对第 k 条射线做最短路径射线追踪 [path_x, path_y] shortest_path(rays(k).sx, rays(k).sy, ... rays(k).rx, rays(k).ry, slowness_grid); % 2. 把连续路径离散化到网格单元, 累加每个网格内的路径长度 for m 1:length(path_x)-1 % 取路径中点判断落在哪个网格 xmid (path_x(m) path_x(m1)) / 2; ymid (path_y(m) path_y(m1)) / 2; ix min(nx, max(1, floor((xmid - x0) / dx) 1)); iy min(ny, max(1, floor((ymid - y0) / dy) 1)); seg_len sqrt((path_x(m1)-path_x(m))^2 (path_y(m1)-path_y(m))^2); col (iy-1)*nx ix; A(k, col) A(k, col) seg_len; end end这段代码里shortest_path 函数需要自己实现一般用多阶 Dijkstra 或基于费马原理的快速推进法FMM替代。判断路径中点落在哪个网格时要特别注意坐标系对齐x0 和 y0 是网格区域左下角坐标如果没对齐会导致路径被安放到相邻网格反演结果整体错位。A 矩阵使用稀疏存储是必须的因为真实数据里射线数量往往几千条网格数上万密集矩阵根本存不下。2.3 为什么用初至走时而不是全部波形题目里的「走时」两个字不是随便写的。实际电磁波或者声波数据里包含的信息远不止初至时间还有振幅、频率、相位等但这些信息对噪声的敏感度和正演模拟复杂度完全不同。走时信息最大的优势是数据预处理代价低、抗噪能力强、正演又快又稳一条射线一趟正演不过是解一个最短路径问题。这个特性在工程尺度问题里极其重要因为场地数据往往脏波形信噪比低硬去拟合完整波形的反演全波形反演极其容易翻车。当然这也带来一个代价走时层析的分辨率受限于波长和射线覆盖密度无法像全波形反演那样恢复细节。对工程检测来说找到异常体的位置和大致边界往往就够了这不值得为了细部特征去承担一个高得多的计算量。另一层考虑是初至拾取的误差会直接进到反演里后面章节会专门讨论粗差与噪声的影响。如果拾取的走时本身错得离谱再精细的算法也没法还原真实介质。3. 在MATLAB里跑通反演主流程从数据准备到成像输出3.1 数据组织与读取走时文件如何规范化拿到的原始数据通常是发射点坐标、接收点坐标和对应的初至走时有可能几百条也有可能好几千条。动手写反演之前我一般先把数据整成统一的 MAT 结构体。数据格式之中的坑在坐标单位不少现场记录把坐标写成厘米或者工程坐标而走时写成毫秒、微秒混用之后慢度单位就会乱掉。统一用国际单位制米、秒是省心的选择。% 统一数据结构的示例 % src: 发射点坐标 [x1 y1; x2 y2; ...] 单位米 % rec: 接收点坐标 [x1 y1; x2 y2; ...] 单位米 % tt: 走时数值, 单位秒, 与观测记录一一对应 data struct(); data.src src; % 或用 load(src.txt) data.rec rec; data.tt tt; % 检查坐标范围, 若发现坐标量级异常则做缩放 fprintf(坐标范围: x [%.1f, %.1f], y [%.1f, %.1f]\n, ... min(src(:,1)), max(src(:,1)), min(src(:,2)), max(src(:,2))); fprintf(走时范围: [%.6f, %.6f] s\n, min(tt), max(tt));走时数据里经常混进去个别异常大的值多半是初至拾错了或者触发延迟反演前必须先做直方图检查。我在很多场合下会先看走时随距离的关系如果同一个距离上走时分布特别离散说明数据质量堪忧需要回到拾取环节复查。反演程序本身不会替你判断数据好坏这一步省了后面全是白算。3.2 用阻尼最小二乘求解反演方程灵敏度矩阵构建好了数据也准备好了反演数学模型就是求解d A·s 噪声。这里的 s 是个大向量直接用反斜杠符号 s A\d 解超定方程组是可以的但结果往往惨不忍睹因为 A 对慢度分布变化的响应极不均匀——没有射线穿过的网格在方程组里对应全零列解出来就是无意义的噪声值。这时候需要引入规则化常见手法是阻尼最小二乘s (AᵀA λLᵀL)⁻¹ Aᵀ d其中 L 可以是单位矩阵 I也可以是平滑算子如拉普拉斯算子。λ 是阻尼系数控制数据拟合与模型平滑之间的权重。L 选 I 最简单相当于给慢度加了一个模长限制L 选拉普拉斯算子则强制相邻网格慢度不要突变成像更平滑。工程反演里我一般会用二阶差分算子做 L成像结果更接近真实介质慢变假设。% 阻尼最小二乘求解核心代码 % A 为灵敏度矩阵, d 为观测走时向量, lam 为阻尼系数 % L 为平滑算子, 这里构造一个简单的横向二阶差分 Lx sparse(nx*ny, nx*ny); for iy 1:ny for ix 2:nx-1 col (iy-1)*nx ix; Lx(col, col) -2; Lx(col, col-1) 1; Lx(col, col1) 1; end end % 纵向二阶差分 Ly sparse(nx*ny, nx*ny); for ix 1:nx for iy 2:ny-1 col (iy-1)*nx ix; Ly(col, col) -2; Ly(col, col-nx) 1; Ly(col, colnx) 1; end end L [Lx; Ly]; % 组装法方程并求解 ATA A * A; ATd A * d; model (ATA lam^2 * (L*L)) \ ATd; % 还原为二维慢度/速度场 slow reshape(model, nx, ny); vel 1 ./ slow;这里 lam^2 的写法是有讲究的因为 LL 的量纲与 AA 不同把阻尼系数平方后再乘规则化矩阵可以让 λ 的数值在不同网格剖分下保持一个相对直观的范围。我一般先试 λ 0.1 再上下调两个数量级观察成像稳定度具体判定标准会在下一章细说。3.3 迭代式反演与残差监控走时层析前两节描述的是线性化反演。真实现场里射线路径本身取决于慢度分布而慢度分布恰恰是待求量这就产生了非线性。通常做法是走「正演-反演-更新模型-重新正演」的迭代流程先用猜测的初始慢度模型比如均匀背景算出射线路径和合成走时再求观测走时与合成走时的残差利用残差做一次线性反演得到慢度修正量更新模型后回到正演环节。这样反复迭代 5 到 10 次残差会逐步收敛。% 迭代反演主循环 % vel_init: 初始速度模型(通常用均匀介质) % max_iter: 最大迭代次数 model 1 ./ vel_init(:); for iter 1:max_iter % 基于当前模型重新计算灵敏度矩阵(射线追踪) A build_sensitivity_matrix(model, data); % 计算合成走时与残差 t_syn A * model; residual data.tt - t_syn; % 用阻尼最小二乘求修正量 dmodel solve_damped_ls(A, residual, lam, L); % 更新模型 model model dmodel; % 输出当前迭代残差RMS, 用于判断收敛 rms(iter) sqrt(mean(residual.^2)); fprintf(iter %d, RMS residual %.6f s\n, iter, rms(iter)); end每轮迭代重新做射线追踪是必要的也是必要的代价。如果只在一开始固定射线路径、后面迭代不再更新路径就成了线性反演遇到速度异常较强的介质会计算出明显偏差——射线走的路其实绕过了异常体线性近似不再成立。判断迭代是否该停看残差 RMS 是否进入平台期一般前 3 次迭代下降很快后面下降幅度很小再迭代下去只会放大噪声。实践中我习惯把 max_iter 设在 8 到 12既保证路径更新充分又不至于让过拟合吃掉反演质量。4. 反演参数怎么设规则化、网格剖分与迭代终止4.1 阻尼系数与规则化强度选法阻尼系数是走时层析里最需要靠经验打磨的参数。λ 太小第二项平滑约束的约束力趋近于零反演完全贴合数据噪声会直接被映射成高速/低速的假异常λ 太大第二项占主导解得出来的模型几乎是一条平坦直线什么异常都看不到。合适的 λ 在两者之间取一个折中——数据残差能接受模型也不至于出现过大的速度跳变。判断 λ 是否合适有一个很实际的办法做一个 λ 扫描比如 0.01 到 100 按对数间隔取 10 个值每个 λ 下跑一次单步线性反演画出残差和模型粗糙度随 λ 的变化曲线取残差开始快速上升之前那段平滑区域里的拐点值。不需要每次都做完整迭代单步反演就够选 λ。我常用的量化指标是模型粗糙度定义为 ||L·s||²它衡量相邻网格之间的差异总量。把 λ 从大到小推进你会看到残差一路下降、模型粗糙度不断上升。工程上一般取残差落在观测走时噪声水平附近的 λ例如拾取的走时误差在 0.5 微秒量级就找让拟合残差 RMS 大约等于这个水平的 λ 值既拟合到数据噪声下限又保持模型尽量平滑。4.2 网格剖分尺寸的选择逻辑网格尺寸直接决定未知数个数和射线穿过网格的路径长度记录精度。剖得太粗一个网格内包含多种介质反演结果的横向边界被严重平滑异常体边界模糊成一片剖得太细未知数增多灵敏度矩阵稀疏度加大而且射线覆盖密度不足的网格会变成反演盲区。网格大小的选择不是单纯跟着异常体最小尺寸走而是要结合射线密度来定每个网格内至少要有几条射线穿过否则解在那个位置不可靠。实际工程中边界网格往往射线很少我一般会在网格剖分时预留一圈虚拟零慢度外扩区域或者把目标区域整体向外扩半个网格让边缘网格多少有点射线覆盖。网格尺寸还直接影响阻尼系数 λ 的量纲前面代码里把 λ 做了平方处理就是为了让 λ 对不同网格尺寸不至于变化过大。经验上说网格尺寸约等于预估异常体尺寸的四分之一到六分之一比较合适太小了除了让成像变花哨没有实际意义。另外一个要注意的是不同方向网格大小不一定一致。有的观测系统比如地表发射、井中接收射线密集分布在横向纵向射线稀疏这时可以适当拉大纵向网格尺寸、减少纵向未知数避免大量欠定网格白白拉高求解难度。MATLAB 里 reshape 以后的方向顺序要跟灵敏度矩阵里的列排序严格一致这一步搞反了图像会转 90 度排查时非常坑。4.3 迭代终止条件别让反演过头迭代次数不是一个需要精调的参数但需要一个终止标准。我一般准备两个条件哪个先满足就停残差下降率低于某个阈值如上一步残差的 2%或者迭代次数到达上限。后一个只是安全网真正决定质量的是第一个。实际操作里初始模型如果给得不准前几次迭代残差下降非常明显到 5 次以后基本不降这时停下来的反演结果是最干净的。还有一种常见情况是残差降到一定程度后反而开始回升这通常意味着模型被过度修正某些区域的慢度调整过头导致重新正演后的合成走时更差了。遇到这种情况检查一下每次迭代后的模型最大修正量如果某一步修正量突然超过正常值的几倍就要考虑这次反演的发散风险。处理办法是给每一步的修正量加一个压缩因子比如限制单步修正不超过当前模型最大值的 10%虽然会牺牲一些收敛速度但不会让反演彻底失控。5. 走时层析反演的常见问题与避坑记录5.1 网格过细导致棋盘格噪声现象反演出来的速度剖面出现规则分布的小块高低速交替视觉上很像棋盘格与真实地质分布明显不符。 原因网格剖得太细未知数数量相对数据量偏多在没有足够射线覆盖的网格位置上数据中的噪声被当作真实信号拟合进模型形成自发性的假异常。 解决降低网格密度或增大阻尼系数或同时给模型加入更重的平滑约束。我在处理一个混凝土结构检测数据时把网格从 0.05 米加密到 0.02 米后棋盘格就出现了退回 0.05 米再加大 λ 后成像干净得多。另一个关键检查点是射线覆盖图把每条射线经过的网格计数画出来成像异常只出现在覆盖数高的区域才能放心低覆盖区域的数值基本不可信。5.2 反演发散或残差震荡现象迭代过程中残差不降反升或者在一个值附近来回跳动速度模型出现离谱的大值。 原因初始模型与实际介质差异过大射线追踪在错误模型下计算出错误的路径进而反演出错误的修正量更新后模型更差另一个常见原因是阻尼系数设置过小或单步修正量没有限幅。 解决先改用均匀初始模型让射线路径基于平均速度计算通常比拍脑袋给一个分层的初始模型更稳更直接的办法是对修正量施加限幅dmodel(dmodel cap) cap把每步最大速度改变限制在 15% 以内。这会让收敛慢一点但保证过程可控。我在检查残差震荡时还会看一个指标数据中是否存在个别与周围走时差异极大的野值。一只野鼠咬断了触发线、现场设备漏触发都能让一条走时记录差出几个数量级这一条射线变成反演中的杠杆点把整个解拉偏。定位这种野值的方法是做正演后算出每条射线的残差残差超过整体 RMS 3 倍以上的射线标记出来逐个核对。5.3 射线覆盖不均匀造成的成像盲区现象某些区域反演后完全无变化始终保持初始模型仔细看射线图这些区域根本没有射线穿过。 原因观测系统几何布局上的局限——发射点和接收点都在同一侧时另一侧的射线密度天然为零。 解决最好的办法是在数据采集阶段就做一次覆盖仿真用均匀介质的射线追踪画出覆盖图再决定测线布置对于已经采完的数据反演时要在覆盖数为零的网格上打掩码避免程序在那类网格上硬解出无意义的数值。成像时在这些区域加上透明掩膜不让空区域干扰判断。5.4 走时拾取误差与反演结果偏移现象反演出的异常体位置与真实位置在水平方向或深度方向有系统性偏移。 原因初至拾取存在系统性延迟比如触发信号晚到、仪器零时没有校零整条测线的走时都偏大或偏小。均匀的零时偏移通常只会影响速度绝对值但不同测段偏移不一致时就会变成路径长度上的假异常。 解决在反演中显式估计一个零时校正量每炮或每条射线一个 t0 偏移与慢度一起联合反演或在反演前先做一致性校正——取同一对收发路径在不同时刻重复观测的走时差标定出仪器延迟。这个坑在电磁波数据里特别常见因为电子设备的触发链路延迟不像地震仪那样容易校零很多现场人员把这当成噪声忽略了。6. 用棋盘格恢复测试验证反演程序的可靠性程序写完了、成像也出来了但怎么证明反演结果可信到这个规模最直接的办法是棋盘格恢复测试checkerboard recovery test。构造一个速度模型在其中交替放置高低速方块像棋盘一样排列用同一套观测系统在这个棋盘模型上正演出合成走时再拿这套走时去做反演。如果反演能恢复出棋盘图案说明观测系统配合这套参数在相应尺度上有分辨能力如果图案糊掉或者完全恢复不出来那就要调整观测参数或反演策略了。% 棋盘格恢复测试 % 注意: 这里用的正演就是反演程序里的射线追踪走时积分 % 因此该测试检验的是反演参数与观测系统, 不检验正演本身 nx 40; ny 40; dx 0.05; base_slow ones(nx, ny) / 1500; % 背景速度 1500 m/s checker base_slow; for iy 1:ny for ix 1:nx if mod(floor((ix-1)/5) floor((iy-1)/5), 2) 0 checker(ix, iy) 1 / 1800; % 高速异常 1800 m/s else checker(ix, iy) 1 / 1300; % 低速异常 1300 m/s end end end % 用棋盘模型正演合成走时, 然后加入轻微高斯噪声 syn_tt A_dense * checker(:); obs_tt syn_tt 0.02 * max(abs(syn_tt)) * randn(size(syn_tt)); % 反演并评估 [inv_model, ~] run_tomography(obs_tt, params); figure; subplot(1,2,1); imagesc(checker); title(真实棋盘); subplot(1,2,2); imagesc(reshape(inv_model, nx, ny)); title(恢复棋盘);棋盘格测试真正检验的是这套反演参数在多大空间尺度上是可靠的。测试步骤分为四步先在同样的收发布置下用棋盘模型正演得到合成走时其次在合成走时中加入与现场噪声水平接近的高斯噪声再按原参数跑完整迭代反演最后对比恢复出的棋盘格边界是否清晰、位置是否准确。如果恢复效果勉强甚至失败需要缩小网格尺寸、增大阻尼或者重新审视射线覆盖密度是否满足小尺度分辨要求。跑测试的时候我一般会连续检查三个尺度的棋盘格尺寸比如异常体尺寸分别取网格边长的 1 倍、2 倍、4 倍得到一套分辨率随尺度下降的数据。把这些结果和现场最终成像放一起看每个异常体到底能不能信、可信到什么程度心里就有底了。实际操作里我还会把棋盘格测试用的走时数据里加进一条明显的粗差看看反演程序在脏数据下如何表现——一个稳定的程序应该让粗差影响只局限在局部区域而不是污染整个模型。这几年处理过好几套不同来源的走时数据越来越觉得走时层析返演这东西程序写得完整不难真正难的是知道程序每个环节在什么条件下会失效。宁可多花半小时做一次棋盘格测试也别拿着一次成像就去下结论。希望这些思路能帮你在自己的反演数据上少走弯路把更多时间花在对结果的解释上。本文还有配套的精品资源点击获取
