简介：一套用于电磁波走时层析成像的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 是观测走时向量，长度是 m；s 是待求慢度向量，长度是 n；A 是灵敏度矩阵（也叫射线路径矩阵），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(m+1)) / 2; ymid = (path_y(m) + path_y(m+1)) / 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(m+1)-path_x(m))^2 + (path_y(m+1)-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, col+1) = 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, col+nx) = 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 的写法是有讲究的，因为 L'L 的量纲与 A'A 不同，把阻尼系数平方后再乘规则化矩阵，可以让 λ 的数值在不同网格剖分下保持一个相对直观的范围。我一般先试 λ = 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 倍，得到一套分辨率随尺度下降的数据。把这些结果和现场最终成像放一起看，每个异常体到底能不能信、可信到什么程度，心里就有底了。实际操作里，我还会把棋盘格测试用的走时数据里加进一条明显的粗差，看看反演程序在脏数据下如何表现——一个稳定的程序应该让粗差影响只局限在局部区域，而不是污染整个模型。
这几年处理过好几套不同来源的走时数据，越来越觉得走时层析返演这东西，程序写得完整不难，真正难的是知道程序每个环节在什么条件下会失效。宁可多花半小时做一次棋盘格测试，也别拿着一次成像就去下结论。希望这些思路能帮你在自己的反演数据上少走弯路，把更多时间花在对结果的解释上。
本文还有配套的精品资源，点击获取
