简介这是一套基于MATLAB的射线反演与SIRT迭代重建实现面向地球物理勘探、医学成像及科学计算领域的研究者与学习者。在地球物理中常利用地震波走时数据反演地下速度结构该资源即围绕这一常见问题以二维射线追踪为核心通过模拟射线路径、计算旅行时间并利用SIRT算法反复修正模型参数可处理非均匀介质中的走时反演与速度结构重建任务。压缩包共10个文件全部为.m脚本整体仅6KB内含SIRT迭代核心、旅行时计算、射线搜索、误差分析以及初始模型构建等模块并配备扇形投影与格林函数计算工具函数划分清晰且便于二次开发。目前已有197人学习下载。该资源覆盖从初始模型构建、射线寻径、正演旅行时到SIRT反演及误差分析的完整流程可直接在MATLAB中运行调试通过修改初始模型与观测参数即可快速对比反演效果适合需要上手射线类反演算法、研究SIRT收敛特性或开展相关课程实验的读者。1. GPU 与 SIRT 相遇时射线反演最慢的那一步不在反演本身先给结论用 Matlab 写 2D 射线反演真正卡时间的不是 SIRT 迭代更新而是正演的射线追踪路径计算。SIRTSimultaneous Iterative Reconstruction Technique同时迭代重建每一次迭代都要把每条射线的走时残差沿路径摊回网格路径不重新追踪、只重算走时的情况下一次迭代其实是一次稀疏矩阵乘。可网格一旦细化到几百乘几百、射线数量上千路径计算和稀疏矩阵装配就会把单次迭代拖到秒级几十轮迭代后训练等不起。本文围绕的是一条常见技术链路RayTrace2D 形式的二维射线追踪做正演配合 SIRT 做走时反演层析全程落在 Matlab 环境里再用 GPU 把最贵的计算段搬上去。看这篇东西的人多半在做超声 CT、跨孔地震、混凝土空洞探测或者课程论文里刚好接了一个「走时反演」的任务。这类问题网格不大、数据不海量但迭代次数多、路径耦合明显是一个很适合用 GPU 加速但又不能无脑把代码扔给gpuArray的典型场景。下面从正演方程开始逐步把这条路走通。2. 2D 射线反演的走时方程与灵敏度矩阵构造2.1 慢度域里的走时方程与 SIRT 的更新式射线走时反演最常见的控制方程是t L·s其中 t 是走时向量每条射线一个值s 是慢度向量速度的倒数单位 s/mL 是灵敏度矩阵或者说路径长度矩阵。L 的行对应一条射线列对应一个网格单元元素值为该射线在这个单元内穿过的长度。注意这里用的是慢度而不是速度原因有两个一是走时对慢度是线性的方程写成线性系统后能直接用代数重建二是反演更新时对慢度做加减法不会像对速度做加法那样产生分母为零的风险。SIRT 的更新公式常见有两种写法。一种是像素驱动pixel-driven视角对每个网格单元把所有穿过它的射线的残差加权求和后更新。另一种是射线驱动ray-driven视角把每条射线残差均匀分摊到路径经过的网格里。实际代码里多半混合使用核心更新表达式是s_new(j) s_old(j) λ · (Σ_i [ L(i,j) · Δt_i / P_i ]) / (Σ_i L(i,j) ε)这里 Δt_i t_obs(i) − t_cal(i) 是第 i 条射线的走时残差P_i 是第 i 条射线的路径总长度。λ 是松弛因子ε 是防止除零的小量。这个形式的好处是当某条射线在网格内路径较短但残差很大时它对每个网格的贡献会被 P_i 归一化避免长路径射线主导更新。实际使用中要特别注意所有网格被射线经过的次数差异巨大角点单元可能只有一条射线穿过而中心区域有几十条直接累加会导致边界单元慢度被抬得很夸张。2.2 用 RayTrace2D 思路给网格建模与射线路径编码2.2.1 规则网格下的索引设计与路径采样二维走时反演最常用的建模方案还是规则网格把探测区域切成 nx × ny 个矩形单元每个单元内慢度恒定射线在各单元交界处发生折射如果要考虑折射或者直接按直线段穿过衰减层析、初代超声 CT 常用。网格建好后最关键的是把「一条射线穿过哪些单元、各穿多长」这个信息编码好。我一般用一个 N×3 的数组保存一条射线路径第 1、2 列是穿过的网格索引第 3 列是对应长度。全部射线存成一个 cell避免用稠密矩阵浪费空间。下面是二维均匀网格下、给定首尾坐标后按直线段采样路径的 Matlab 实现function path straightRayPath(p1, p2, dx, dy, nx, ny) % p1, p2: [x, y] 起终点坐标物理单位 % dx, dy: 网格单元尺寸nx, ny: 网格列数、行数 % path: 结构体包含 idx 与 len % 把射线离散采样每个采样点对应一个网格 Nsample round(norm(p2 - p1) / min(dx, dy)) * 4; xs linspace(p1(1), p2(1), Nsample); ys linspace(p1(2), p2(2), Nsample); ix floor(xs / dx) 1; iy floor(ys / dy) 1; ix min(max(ix, 1), nx); % 边界截断 iy min(max(iy, 1), ny); % 合并相邻重复点同一网格内多次采样的长度要合并 idx_sub sub2ind([ny, nx], iy, ix); [uniq_idx, ~, ic] unique(idx_sub, stable); len accumarray(ic, 1) * min(dx, dy) / 4; path.idx uniq_idx; path.len len; end这段代码把射线按空间步长采点采样步长取网格最小边长的四分之一这样在穿过单元斜边时长度误差能控制在可接受范围。unique加accumarray的作用是把同一个单元内的多个采样点合并成一段路径长度这一步不做的话累加的路径长度会偏大反演出来的慢度整体偏小。实际项目中如果介质的慢度差异超过 20%直线路径假设就会引入系统性误差需要升级成最短路径法或 eikonal 方程求解。RayTrace2D 路线里常见做法是先跑一遍 fermat 原理下的最短路径正演把弯曲路径算出来再用本节的路径编码结构去承载。2.3 Matlab 里先验证正演走时一个最小可运行片段写反演之前务必先验证正演走时算得对不对。方法是给一个慢度已知的均匀模型用直线路径算理论走时与解析解t 距离 / 速度对比。% 均匀模型速度 1500 m/s慢度 1/1500 dx 0.05; dy 0.05; nx 40; ny 20; slow ones(ny, nx) / 1500; % 发射点在左边界中点接收点在右边界中点 p1 [0, 0.5]; p2 [2.0, 0.5]; path straightRayPath(p1, p2, dx, dy, nx, ny); t_calc sum(slow(path.idx) .* path.len); t_ana norm(p2 - p1) * (1 / 1500); fprintf(计算走时 %.6f ms解析走时 %.6f ms误差 %.3f%%\n, ... t_calc*1000, t_ana*1000, abs(t_calc - t_ana)/t_ana*100);注意这里path.idx是按列优先索引的对应sub2ind([ny, nx], iy, ix)在把慢度向量slow转成列向量后可以直接做线性索引。验证通过后把多条射线的路径组装成一个大 cell正演部分就算闭环了。这个最小片段跑出来的误差通常小于 1%如果你发现误差达到 2% 以上先查采样步长是不是太稀再查边界截断逻辑是否正确。3. Matlab 实现 SIRT松弛迭代、阻尼与残差收敛曲线3.1 为什么要用代数重建而不用直接求逆走时层析里灵敏度矩阵 L 的性质很差行数和列数都是几千的量级但每行只有几个非零元素条件数大得离谱直接求s L \ t在 Matlab 里虽然能算但结果会被噪声放大到不可用。SIRT 本质上是 Landweber 迭代的一种变体它不要求矩阵满秩也不要求观测数据完备只要初始模型合理、松弛因子选择得当就能在几十轮迭代内收敛到接近最小二乘解。更重要的是SIRT 很容易加先验约束慢度上下限、平滑项、阻尼项都可以直接塞进迭代式中。对射线反演场景SIRT 还有一个隐式优点射线的覆盖密度天然不均匀SIRT 的「残差按路径长度归一化」和「除以射线经过次数」这两步能有效緩解覆盖稀疏区域的慢度震荡。直接解正规方程时覆盖稀疏区域的列几乎为零逆矩阵会把不可观测量放大到离谱。3.2 一个可直接套用的 sirt 层析实现下面的函数实现了带慢度上下限约束的 SIRT 迭代输入是初始慢度向量、路径 cell 列表以及观测走时function [s, hist] sirt2d(s0, rays, tObs, nIter, lambda, sLow, sHigh) % s0: 初始慢度向量ny*nx, 1 % rays: 结构体数组含 idx 和 len 字段 % tObs: 观测走时向量 % nIter: 最大迭代次数 % lambda: 松弛因子典型范围 0.01 ~ 0.5 % sLow, sHigh: 慢度下限和上限用于物理约束 s s0(:); nRays numel(rays); nCells numel(s); hist zeros(nIter, 1); for k 1:nIter ds zeros(nCells, 1); cnt zeros(nCells, 1); for i 1:nRays idx rays(i).idx; L rays(i).len; tCal sum(s(idx) .* L); res tObs(i) - tCal; pathLen sum(L); if pathLen 0 ds(idx) ds(idx) L(:) * (res / pathLen); cnt(idx) cnt(idx) 1; end end upd lambda * ds ./ max(cnt, 1); upd(~isfinite(upd)) 0; s s upd; s min(max(s, sLow), sHigh); % 投影到物理可行域 hist(k) sqrt(mean((tObs - predictTimes(s, rays)).^2)); end end function t predictTimes(s, rays) nRays numel(rays); t zeros(nRays, 1); for i 1:nRays t(i) sum(s(rays(i).idx) .* rays(i).len); end end这个实现里有三个细节要说明。第一残差是用res / pathLen归一化而不是直接乘L这会让每条射线对整体慢度更新的贡献与自身长度无关只与相对误差有关。第二max(cnt, 1)的目的是避免无射线覆盖的单元被除以零但这同时意味着无覆盖区域的慢度完全不更新它们会保留初始模型的设定值。第三min(max(...))的投影操作非常关键实测中如果不做慢度上下限截断边缘区域的慢度会在二十轮迭代后漂移到离谱值。3.3 松弛因子与阻尼的参数选择松弛因子 λ 的取值直接决定迭代曲线形态。λ 太小收敛慢适合噪声大的数据λ 太大前几轮残差下降快但随后出现震荡高空间频率的慢度波动会越迭代越强。下面是不同 λ 在中等噪声数据下的行为对照λ 取值前 10 轮残差下降50 轮后残差模型平滑度建议场景0.01慢线性下降高最平滑信噪比低、数据量少0.050.1理想低较好通用首选0.20.5前 5 轮快低但震荡出现椒盐噪声数据干净、覆盖密集1.0 以上发散残差反弹严重噪化不推荐实际调参时我一般先用 λ 0.1 跑 20 轮看残差下降率如果每轮下降不足 0.2% 就提升到 0.2如果出现震荡就降到 0.05。还有一种自适应做法每轮记录残差均方根若连续三轮上升就把 λ 减半这个策略能省去手动干预。阻尼方面给更新量加一个随迭代次数衰减的系数也常见本质上是让早期迭代大胆、晚期迭代小心。4. GPU 加速路径选择gpuArray、parfor、还是 mexCUDA4.1 先算清楚瓶颈在哪射线追踪 vs 反演更新把一段 Matlab SIRT 代码直接改成gpuArray之前先用tic/toc把耗时拆开射线追踪路径构造、每条射线的走时计算、残差累加回填、松弛更新。我实测过多次数据点规模在「网格 100×80、射线 800 条」量级时各段耗时大致比例是这样的直线路径构造如果每次迭代都要做占总耗时 50%70%走时预测与残差累加占 20%30%松弛更新与约束投影占 5% 以下关键结论是如果算法设计成「每轮迭代都重新追踪路径」那么 GPU 要加速的核心其实是射线追踪或路径采样而不是矩阵运算。反过来如果路径只在反演前算一次走时预测变成稀疏矩阵乘这时gpuArray才有发挥空间。标题里的 RayTrace2D 如果实现对每条射线做逐点步进采样那这是一段循环嵌套的代码在 GPU 上反而不好写因为每个采样点之间存在数据依赖。备选方案有两个用parfor把不同射线并行起来或者改写成矩阵化批次路径计算一次算出所有射线的采样坐标。4.2 小模型用 gpuArray 改写矩阵更新段固定路径不重追的条件下SIRT 主循环可以写成稀疏矩阵运算这时 GPU 介入成本最低。% 预先装配稀疏灵敏度矩阵 A尺寸 nRays x nCells % A(i, j) 第 i 条射线在 j 单元内的路径长度 A sparse(nRays, nCells); for i 1:nRays A(i, rays(i).idx) rays(i).len; end % 转移到 GPU Ag gpuArray(A); s gpuArray(s0); tObsg gpuArray(tObs); for k 1:nIter tCal Ag * s; res tObsg - tCal; pathLen sum(Ag, 2); % 每条射线总长 resNorm res ./ pathLen; % 归一化残差 ds Ag * resNorm; % 残差按长度回填到网格 cnt sum(Ag, 1); % 每个网格被选中的次数 s s (lambda * ds) ./ max(cnt, 1); s min(max(s, sLow), sHigh); end s gather(s);这段代码把所有主要计算都搬到了 GPU 上。注意其中Ag * resNorm和sum(Ag, 2)是稀疏矩阵在 GPU 上的高效操作Matlab 对 gpuArray 稀疏矩阵乘法支持得不错。但是有两处必须提醒第一gpuArray(A)的前天是 A 已经是 sparse 类型不要把稀疏矩阵转 full 再上传那会直接耗尽显存第二在网格只有几千个单元、射线只有几百条的小规模问题上GPU 版本往往比 CPU 慢因为数据要上传到显存、计算完还要传回来传输开销盖过了计算收益。经验阈值是单元数超过 2 万、射线超过 2000 条时GPU 的优势才明显。4.3 mexCUDA 全 GPU 实现路线与数据结构如果 RayTrace2D 的路径每次迭代都要重算稀疏矩阵方式就失灵了这时靠谱的路子是写 CUDA 内核通过 mex 在 Matlab 里调用。不需要把整个反演搬到 CUDA只需要搬「路径采样 残差回填」这一段。// pathKernel.cu 核心逻辑示意 // 每个 CUDA 线程负责一条射线的路径采样与残差分配 __global__ void pathSirtKernel( const float* p1, const float* p2, // 射线端点 const float* s, // 当前慢度 const float* tObs, // 观测走时 float* ds, // 慢度累积更新量 float* cnt, // 单元覆盖计数 int nRays, int nx, int ny, float dx, float dy, float lambda, int nSamplePerRay) { int i blockIdx.x * blockDim.x threadIdx.x; if (i nRays) return; // 直线采样每个线程独立算整条路径互不依赖 // 残差累加时注意多个射线可能经过同一单元 // 使用 atomicAdd(ds[idx], len * res / pathLen); }这里的关键设计是原子加操作atomicAdd。由于多条射线会映射到同一个网格单元直接写ds[idx] ...会发生数据竞争必须用原子操作累加。代价是累加冲突的线程会串行排队不过对每单元被覆盖次数几十量级的场景开销可忽略。这个路线的加速收益比gpuArray版本大得多尤其是路径需要反复重追的反演问题。但开发成本高需要熟悉 CUDA 编程、会写 mex 接口、还要处理单双精度问题。Tesla 系列这类计算卡跑单精度速度快双精度会明显变慢如果正演的数据量不大单精度累计误差对最终慢度结果的影响通常在千分之一量级可以接受。反过来消费级游戏卡的双精度性能很弱遇到高精度需求时要提前验证。对于不想碰 CUDA 的团队折中方案是parforparfor i 1:nRays % 每条射线的路径独立计算并行化友好 path straightRayPath(p1(i,:), p2(i,:), dx, dy, nx, ny); tCal(i) sum(s(path.idx) .* path.len); endparfor的适用场景是路径计算耗时远大于调度开销且单机多核可用。要注意parfor里不能直接修改共享的ds数组需要用先归约后合并的模式否则 Matlab 会报变量冲突。就收益而言parfor通常能达到 24 倍加速和 GPU 的 1030 倍有量级差距但开发成本几乎为零。5. 收敛判据的实用技巧用相对残差和数值梯度双确认SIRT 迭代什么时候该停最粗糙的做法是「跑满 50 轮」这不靠谱不同模型和数据质量下SIRT 的收敛速度差异极大。先用相对残差下降率做判据记录当前轮的走时均方根残差 rmse_k计算下降率(rmse_{k-1} - rmse_k) / rmse_{k-1}当连续三到五轮下降率小于 1% 时停止迭代。这个判据的优点是鲁棒缺点是可能停在「残差高原」上即下降率很小但模型还没完全收敛。解法是再同时记录噪声水平估计值 sigma当 rmse 降到与 sigma 接近时立刻停因为继续迭代只会拟合数据噪声。第二个非常值得养成的习惯是数值梯度检查。SIRT 本质上是梯度类算法的变体它的更新方向应该接近走时残差对慢度的梯度。检查方法如下随机选一个单元给慢度加一个小扰动 delta_s例如 1e-6重算所有射线的走时用有限差分估算残差变化量把这个变化量与 SIRT 更新公式里该单元的更新量做方向对比。如果两者同号且幅度比例合理说明路径装配和索引逻辑没错如果出现方向相反一定是 idx 或 len 的数据对应错了。这个检查在模型较小时几分钟能跑完能拦下大量「反演看似收敛但结果是错的」的隐蔽 bug。最后一层验证是交叉验证式的怀疑把观测射线随机抽掉 10% 作为验证集用剩下的 90% 反演然后看验证集的走时残差是否与训练集相当。如果验证集残差远大于训练集说明迭代过拟合这时需要调小 λ 或提前停止。这套做法和深度学习里划分验证集的逻辑一脉相承在层析反演里价值很高因为反演结果往往是医生或工程判据的直接输入单纯追求训练集残差为零没有实际意义。本文还有配套的精品资源点击获取 GPU加速SIRT射线走时反演：从Matlab到CUDA的二维层析实践
简介这是一套基于MATLAB的射线反演与SIRT迭代重建实现面向地球物理勘探、医学成像及科学计算领域的研究者与学习者。在地球物理中常利用地震波走时数据反演地下速度结构该资源即围绕这一常见问题以二维射线追踪为核心通过模拟射线路径、计算旅行时间并利用SIRT算法反复修正模型参数可处理非均匀介质中的走时反演与速度结构重建任务。压缩包共10个文件全部为.m脚本整体仅6KB内含SIRT迭代核心、旅行时计算、射线搜索、误差分析以及初始模型构建等模块并配备扇形投影与格林函数计算工具函数划分清晰且便于二次开发。目前已有197人学习下载。该资源覆盖从初始模型构建、射线寻径、正演旅行时到SIRT反演及误差分析的完整流程可直接在MATLAB中运行调试通过修改初始模型与观测参数即可快速对比反演效果适合需要上手射线类反演算法、研究SIRT收敛特性或开展相关课程实验的读者。1. GPU 与 SIRT 相遇时射线反演最慢的那一步不在反演本身先给结论用 Matlab 写 2D 射线反演真正卡时间的不是 SIRT 迭代更新而是正演的射线追踪路径计算。SIRTSimultaneous Iterative Reconstruction Technique同时迭代重建每一次迭代都要把每条射线的走时残差沿路径摊回网格路径不重新追踪、只重算走时的情况下一次迭代其实是一次稀疏矩阵乘。可网格一旦细化到几百乘几百、射线数量上千路径计算和稀疏矩阵装配就会把单次迭代拖到秒级几十轮迭代后训练等不起。本文围绕的是一条常见技术链路RayTrace2D 形式的二维射线追踪做正演配合 SIRT 做走时反演层析全程落在 Matlab 环境里再用 GPU 把最贵的计算段搬上去。看这篇东西的人多半在做超声 CT、跨孔地震、混凝土空洞探测或者课程论文里刚好接了一个「走时反演」的任务。这类问题网格不大、数据不海量但迭代次数多、路径耦合明显是一个很适合用 GPU 加速但又不能无脑把代码扔给gpuArray的典型场景。下面从正演方程开始逐步把这条路走通。2. 2D 射线反演的走时方程与灵敏度矩阵构造2.1 慢度域里的走时方程与 SIRT 的更新式射线走时反演最常见的控制方程是t L·s其中 t 是走时向量每条射线一个值s 是慢度向量速度的倒数单位 s/mL 是灵敏度矩阵或者说路径长度矩阵。L 的行对应一条射线列对应一个网格单元元素值为该射线在这个单元内穿过的长度。注意这里用的是慢度而不是速度原因有两个一是走时对慢度是线性的方程写成线性系统后能直接用代数重建二是反演更新时对慢度做加减法不会像对速度做加法那样产生分母为零的风险。SIRT 的更新公式常见有两种写法。一种是像素驱动pixel-driven视角对每个网格单元把所有穿过它的射线的残差加权求和后更新。另一种是射线驱动ray-driven视角把每条射线残差均匀分摊到路径经过的网格里。实际代码里多半混合使用核心更新表达式是s_new(j) s_old(j) λ · (Σ_i [ L(i,j) · Δt_i / P_i ]) / (Σ_i L(i,j) ε)这里 Δt_i t_obs(i) − t_cal(i) 是第 i 条射线的走时残差P_i 是第 i 条射线的路径总长度。λ 是松弛因子ε 是防止除零的小量。这个形式的好处是当某条射线在网格内路径较短但残差很大时它对每个网格的贡献会被 P_i 归一化避免长路径射线主导更新。实际使用中要特别注意所有网格被射线经过的次数差异巨大角点单元可能只有一条射线穿过而中心区域有几十条直接累加会导致边界单元慢度被抬得很夸张。2.2 用 RayTrace2D 思路给网格建模与射线路径编码2.2.1 规则网格下的索引设计与路径采样二维走时反演最常用的建模方案还是规则网格把探测区域切成 nx × ny 个矩形单元每个单元内慢度恒定射线在各单元交界处发生折射如果要考虑折射或者直接按直线段穿过衰减层析、初代超声 CT 常用。网格建好后最关键的是把「一条射线穿过哪些单元、各穿多长」这个信息编码好。我一般用一个 N×3 的数组保存一条射线路径第 1、2 列是穿过的网格索引第 3 列是对应长度。全部射线存成一个 cell避免用稠密矩阵浪费空间。下面是二维均匀网格下、给定首尾坐标后按直线段采样路径的 Matlab 实现function path straightRayPath(p1, p2, dx, dy, nx, ny) % p1, p2: [x, y] 起终点坐标物理单位 % dx, dy: 网格单元尺寸nx, ny: 网格列数、行数 % path: 结构体包含 idx 与 len % 把射线离散采样每个采样点对应一个网格 Nsample round(norm(p2 - p1) / min(dx, dy)) * 4; xs linspace(p1(1), p2(1), Nsample); ys linspace(p1(2), p2(2), Nsample); ix floor(xs / dx) 1; iy floor(ys / dy) 1; ix min(max(ix, 1), nx); % 边界截断 iy min(max(iy, 1), ny); % 合并相邻重复点同一网格内多次采样的长度要合并 idx_sub sub2ind([ny, nx], iy, ix); [uniq_idx, ~, ic] unique(idx_sub, stable); len accumarray(ic, 1) * min(dx, dy) / 4; path.idx uniq_idx; path.len len; end这段代码把射线按空间步长采点采样步长取网格最小边长的四分之一这样在穿过单元斜边时长度误差能控制在可接受范围。unique加accumarray的作用是把同一个单元内的多个采样点合并成一段路径长度这一步不做的话累加的路径长度会偏大反演出来的慢度整体偏小。实际项目中如果介质的慢度差异超过 20%直线路径假设就会引入系统性误差需要升级成最短路径法或 eikonal 方程求解。RayTrace2D 路线里常见做法是先跑一遍 fermat 原理下的最短路径正演把弯曲路径算出来再用本节的路径编码结构去承载。2.3 Matlab 里先验证正演走时一个最小可运行片段写反演之前务必先验证正演走时算得对不对。方法是给一个慢度已知的均匀模型用直线路径算理论走时与解析解t 距离 / 速度对比。% 均匀模型速度 1500 m/s慢度 1/1500 dx 0.05; dy 0.05; nx 40; ny 20; slow ones(ny, nx) / 1500; % 发射点在左边界中点接收点在右边界中点 p1 [0, 0.5]; p2 [2.0, 0.5]; path straightRayPath(p1, p2, dx, dy, nx, ny); t_calc sum(slow(path.idx) .* path.len); t_ana norm(p2 - p1) * (1 / 1500); fprintf(计算走时 %.6f ms解析走时 %.6f ms误差 %.3f%%\n, ... t_calc*1000, t_ana*1000, abs(t_calc - t_ana)/t_ana*100);注意这里path.idx是按列优先索引的对应sub2ind([ny, nx], iy, ix)在把慢度向量slow转成列向量后可以直接做线性索引。验证通过后把多条射线的路径组装成一个大 cell正演部分就算闭环了。这个最小片段跑出来的误差通常小于 1%如果你发现误差达到 2% 以上先查采样步长是不是太稀再查边界截断逻辑是否正确。3. Matlab 实现 SIRT松弛迭代、阻尼与残差收敛曲线3.1 为什么要用代数重建而不用直接求逆走时层析里灵敏度矩阵 L 的性质很差行数和列数都是几千的量级但每行只有几个非零元素条件数大得离谱直接求s L \ t在 Matlab 里虽然能算但结果会被噪声放大到不可用。SIRT 本质上是 Landweber 迭代的一种变体它不要求矩阵满秩也不要求观测数据完备只要初始模型合理、松弛因子选择得当就能在几十轮迭代内收敛到接近最小二乘解。更重要的是SIRT 很容易加先验约束慢度上下限、平滑项、阻尼项都可以直接塞进迭代式中。对射线反演场景SIRT 还有一个隐式优点射线的覆盖密度天然不均匀SIRT 的「残差按路径长度归一化」和「除以射线经过次数」这两步能有效緩解覆盖稀疏区域的慢度震荡。直接解正规方程时覆盖稀疏区域的列几乎为零逆矩阵会把不可观测量放大到离谱。3.2 一个可直接套用的 sirt 层析实现下面的函数实现了带慢度上下限约束的 SIRT 迭代输入是初始慢度向量、路径 cell 列表以及观测走时function [s, hist] sirt2d(s0, rays, tObs, nIter, lambda, sLow, sHigh) % s0: 初始慢度向量ny*nx, 1 % rays: 结构体数组含 idx 和 len 字段 % tObs: 观测走时向量 % nIter: 最大迭代次数 % lambda: 松弛因子典型范围 0.01 ~ 0.5 % sLow, sHigh: 慢度下限和上限用于物理约束 s s0(:); nRays numel(rays); nCells numel(s); hist zeros(nIter, 1); for k 1:nIter ds zeros(nCells, 1); cnt zeros(nCells, 1); for i 1:nRays idx rays(i).idx; L rays(i).len; tCal sum(s(idx) .* L); res tObs(i) - tCal; pathLen sum(L); if pathLen 0 ds(idx) ds(idx) L(:) * (res / pathLen); cnt(idx) cnt(idx) 1; end end upd lambda * ds ./ max(cnt, 1); upd(~isfinite(upd)) 0; s s upd; s min(max(s, sLow), sHigh); % 投影到物理可行域 hist(k) sqrt(mean((tObs - predictTimes(s, rays)).^2)); end end function t predictTimes(s, rays) nRays numel(rays); t zeros(nRays, 1); for i 1:nRays t(i) sum(s(rays(i).idx) .* rays(i).len); end end这个实现里有三个细节要说明。第一残差是用res / pathLen归一化而不是直接乘L这会让每条射线对整体慢度更新的贡献与自身长度无关只与相对误差有关。第二max(cnt, 1)的目的是避免无射线覆盖的单元被除以零但这同时意味着无覆盖区域的慢度完全不更新它们会保留初始模型的设定值。第三min(max(...))的投影操作非常关键实测中如果不做慢度上下限截断边缘区域的慢度会在二十轮迭代后漂移到离谱值。3.3 松弛因子与阻尼的参数选择松弛因子 λ 的取值直接决定迭代曲线形态。λ 太小收敛慢适合噪声大的数据λ 太大前几轮残差下降快但随后出现震荡高空间频率的慢度波动会越迭代越强。下面是不同 λ 在中等噪声数据下的行为对照λ 取值前 10 轮残差下降50 轮后残差模型平滑度建议场景0.01慢线性下降高最平滑信噪比低、数据量少0.050.1理想低较好通用首选0.20.5前 5 轮快低但震荡出现椒盐噪声数据干净、覆盖密集1.0 以上发散残差反弹严重噪化不推荐实际调参时我一般先用 λ 0.1 跑 20 轮看残差下降率如果每轮下降不足 0.2% 就提升到 0.2如果出现震荡就降到 0.05。还有一种自适应做法每轮记录残差均方根若连续三轮上升就把 λ 减半这个策略能省去手动干预。阻尼方面给更新量加一个随迭代次数衰减的系数也常见本质上是让早期迭代大胆、晚期迭代小心。4. GPU 加速路径选择gpuArray、parfor、还是 mexCUDA4.1 先算清楚瓶颈在哪射线追踪 vs 反演更新把一段 Matlab SIRT 代码直接改成gpuArray之前先用tic/toc把耗时拆开射线追踪路径构造、每条射线的走时计算、残差累加回填、松弛更新。我实测过多次数据点规模在「网格 100×80、射线 800 条」量级时各段耗时大致比例是这样的直线路径构造如果每次迭代都要做占总耗时 50%70%走时预测与残差累加占 20%30%松弛更新与约束投影占 5% 以下关键结论是如果算法设计成「每轮迭代都重新追踪路径」那么 GPU 要加速的核心其实是射线追踪或路径采样而不是矩阵运算。反过来如果路径只在反演前算一次走时预测变成稀疏矩阵乘这时gpuArray才有发挥空间。标题里的 RayTrace2D 如果实现对每条射线做逐点步进采样那这是一段循环嵌套的代码在 GPU 上反而不好写因为每个采样点之间存在数据依赖。备选方案有两个用parfor把不同射线并行起来或者改写成矩阵化批次路径计算一次算出所有射线的采样坐标。4.2 小模型用 gpuArray 改写矩阵更新段固定路径不重追的条件下SIRT 主循环可以写成稀疏矩阵运算这时 GPU 介入成本最低。% 预先装配稀疏灵敏度矩阵 A尺寸 nRays x nCells % A(i, j) 第 i 条射线在 j 单元内的路径长度 A sparse(nRays, nCells); for i 1:nRays A(i, rays(i).idx) rays(i).len; end % 转移到 GPU Ag gpuArray(A); s gpuArray(s0); tObsg gpuArray(tObs); for k 1:nIter tCal Ag * s; res tObsg - tCal; pathLen sum(Ag, 2); % 每条射线总长 resNorm res ./ pathLen; % 归一化残差 ds Ag * resNorm; % 残差按长度回填到网格 cnt sum(Ag, 1); % 每个网格被选中的次数 s s (lambda * ds) ./ max(cnt, 1); s min(max(s, sLow), sHigh); end s gather(s);这段代码把所有主要计算都搬到了 GPU 上。注意其中Ag * resNorm和sum(Ag, 2)是稀疏矩阵在 GPU 上的高效操作Matlab 对 gpuArray 稀疏矩阵乘法支持得不错。但是有两处必须提醒第一gpuArray(A)的前天是 A 已经是 sparse 类型不要把稀疏矩阵转 full 再上传那会直接耗尽显存第二在网格只有几千个单元、射线只有几百条的小规模问题上GPU 版本往往比 CPU 慢因为数据要上传到显存、计算完还要传回来传输开销盖过了计算收益。经验阈值是单元数超过 2 万、射线超过 2000 条时GPU 的优势才明显。4.3 mexCUDA 全 GPU 实现路线与数据结构如果 RayTrace2D 的路径每次迭代都要重算稀疏矩阵方式就失灵了这时靠谱的路子是写 CUDA 内核通过 mex 在 Matlab 里调用。不需要把整个反演搬到 CUDA只需要搬「路径采样 残差回填」这一段。// pathKernel.cu 核心逻辑示意 // 每个 CUDA 线程负责一条射线的路径采样与残差分配 __global__ void pathSirtKernel( const float* p1, const float* p2, // 射线端点 const float* s, // 当前慢度 const float* tObs, // 观测走时 float* ds, // 慢度累积更新量 float* cnt, // 单元覆盖计数 int nRays, int nx, int ny, float dx, float dy, float lambda, int nSamplePerRay) { int i blockIdx.x * blockDim.x threadIdx.x; if (i nRays) return; // 直线采样每个线程独立算整条路径互不依赖 // 残差累加时注意多个射线可能经过同一单元 // 使用 atomicAdd(ds[idx], len * res / pathLen); }这里的关键设计是原子加操作atomicAdd。由于多条射线会映射到同一个网格单元直接写ds[idx] ...会发生数据竞争必须用原子操作累加。代价是累加冲突的线程会串行排队不过对每单元被覆盖次数几十量级的场景开销可忽略。这个路线的加速收益比gpuArray版本大得多尤其是路径需要反复重追的反演问题。但开发成本高需要熟悉 CUDA 编程、会写 mex 接口、还要处理单双精度问题。Tesla 系列这类计算卡跑单精度速度快双精度会明显变慢如果正演的数据量不大单精度累计误差对最终慢度结果的影响通常在千分之一量级可以接受。反过来消费级游戏卡的双精度性能很弱遇到高精度需求时要提前验证。对于不想碰 CUDA 的团队折中方案是parforparfor i 1:nRays % 每条射线的路径独立计算并行化友好 path straightRayPath(p1(i,:), p2(i,:), dx, dy, nx, ny); tCal(i) sum(s(path.idx) .* path.len); endparfor的适用场景是路径计算耗时远大于调度开销且单机多核可用。要注意parfor里不能直接修改共享的ds数组需要用先归约后合并的模式否则 Matlab 会报变量冲突。就收益而言parfor通常能达到 24 倍加速和 GPU 的 1030 倍有量级差距但开发成本几乎为零。5. 收敛判据的实用技巧用相对残差和数值梯度双确认SIRT 迭代什么时候该停最粗糙的做法是「跑满 50 轮」这不靠谱不同模型和数据质量下SIRT 的收敛速度差异极大。先用相对残差下降率做判据记录当前轮的走时均方根残差 rmse_k计算下降率(rmse_{k-1} - rmse_k) / rmse_{k-1}当连续三到五轮下降率小于 1% 时停止迭代。这个判据的优点是鲁棒缺点是可能停在「残差高原」上即下降率很小但模型还没完全收敛。解法是再同时记录噪声水平估计值 sigma当 rmse 降到与 sigma 接近时立刻停因为继续迭代只会拟合数据噪声。第二个非常值得养成的习惯是数值梯度检查。SIRT 本质上是梯度类算法的变体它的更新方向应该接近走时残差对慢度的梯度。检查方法如下随机选一个单元给慢度加一个小扰动 delta_s例如 1e-6重算所有射线的走时用有限差分估算残差变化量把这个变化量与 SIRT 更新公式里该单元的更新量做方向对比。如果两者同号且幅度比例合理说明路径装配和索引逻辑没错如果出现方向相反一定是 idx 或 len 的数据对应错了。这个检查在模型较小时几分钟能跑完能拦下大量「反演看似收敛但结果是错的」的隐蔽 bug。最后一层验证是交叉验证式的怀疑把观测射线随机抽掉 10% 作为验证集用剩下的 90% 反演然后看验证集的走时残差是否与训练集相当。如果验证集残差远大于训练集说明迭代过拟合这时需要调小 λ 或提前停止。这套做法和深度学习里划分验证集的逻辑一脉相承在层析反演里价值很高因为反演结果往往是医生或工程判据的直接输入单纯追求训练集残差为零没有实际意义。本文还有配套的精品资源点击获取