简介面向地球物理专业的正则化反演MATLAB程序包适合新手及有一定数值计算经验的开发人员可用于解决反演问题中的不适定性与数据噪声干扰也适用于课堂演示、课题实验与工程反演任务。压缩包内共61个文件主体为57个.m源码文件另附PDF说明文档、变更记录与运行日志等辅助内容整体大小约1.11MB按正则化方法、测试问题、参数选择与求解器等模块划分便于按需调用。资源已有915人浏览学习涵盖Tikhonov正则化、TSVD、截断完全最小二乘、广义交叉验证与L曲线选参等常用方法并内置多个标准测试问题如热传导、图像模糊等典型场景可帮助读者快速理解反演算法原理、对比不同正则化策略的效果。代码采用模块化函数设计支持修改与扩展既能从基础示例入手逐步演练也可作为科研与工程数值反演的参考工具。1. 正则化反演不是多算几步而是给物探数值反演装上“避雷针”做重磁电震反演的人都有过这种经历实测数据拿来直接求解 Axb得到的剖面像心电图模型值上下乱跳幅度比真实物性差几个数量级网格加密一点解就完全换一张脸。问题不在数据质量而在反演本身是离散不适定问题——系数矩阵的奇异值可以跨十几个数量级观测噪声在求逆过程中被成倍放大。正则化反演就是在这种数学结构里加约束、做截断、选参数让解在拟合观测与保持合理之间取折中。这套 MATLAB 源码覆盖 Tikhonov、TSVD、MTSVD、迭代正则化、L 曲线与 GCV 参数选择等全套工具。新手能照着跑通标准流程有经验的人也可以直接替换正演算子用到重力、磁法、电阻率和地震走时反演的实际数值反演任务里。2. 先看奇异值谱再动手反演离散病态性与 Picard 条件的判断2.1 为什么正则化反演中“精确解”反而是最差解地下物性分布与地面观测之间通常用第一类 Fredholm 积分方程描述。离散之后得到线性系统 AxbA 的条件数用 cond(A)s_max/s_min 来度量。实测反演里这个比值轻松到 10^10 以上。直接最小二乘解里每个分量都被 1/s_i 缩放小奇异值对应的分量把观测噪声放得巨大解偏离真实模型很远。这是数学结构使然与迭代次数、求解器精度无关。正则化工具箱里所有方法本质都是在滤波因子 f_i 上做文章。Tikhonov 的 f_is_i^2/(s_i^2λ^2)TSVD 的 f_i 在 i≤k 时取 1、ik 取 0。参数 λ 或 k 一旦确定滤波因子就把大奇异值方向保留、小奇异值方向压制。后面的章节都围绕这个视角展开。2.2 csvd.m 与 picard.m诊断数据里还有没有“可反演信息”先说一个我拿到新数据必做的事在加任何正则化手段之前先对 A 做一次 SVD。工具箱里这个动作由 csvd.m 完成它对矩形矩阵也适用返回 [U,s,V]s 是奇异值向量。接着用 picard.m 画三条曲线|u_i^T b|、s_i、|u_i^T b|/s_i。这里 u_i^T b 是右端项在第 i 个左奇异向量上的投影代表观测数据在该 SVD 分量上的“信息含量”。% 用 shaw 测试问题生成一个标准反演场景 n 128; [A, b, x] shaw(n); % 生成离散积分算子、右端项和真实模型 rng(42); % 固定随机种子保证可复现 e randn(size(b)); % 高斯白噪声 noise 1e-3 * norm(b) / norm(e); bn b noise * e; % 带噪右端项 [U, s, V] csvd(A); % 紧凑 SVDA 不必是方阵 picard(U, s, bn); % 诊断离散 Picard 条件跑完这段后看 picard 图的对数坐标。只要看到 |u_i^T b| 随 i 前几十个点缓慢下降而 s_i 下降更快那么 |u_i^T b|/s_i 那条线必然在第几十个分量后掉头向上。向上的位置就标出了噪声主导区的起点。比值曲线先降后升说明数据满足离散 Picard 条件噪声主导区在尾部如果比值从头到尾都在上升说明观测数据本身信噪比过低任何正则化手段都救不回来这时该回去检查观测系统而不是调参数。csvd 的 s 输出按降序排列长度是 min(m,n)和 U 的列数严格对应所以 picard 里三个量能画在同一组横坐标上。提示picard.m 输出的是一个窗口图需要和 csvd.m 返回的 s 配套使用。不要拿 picard 去诊断原始 A 的条件数它诊断的是 A、b 联合决定的离散 Picard 条件。2.3 网格剖分对病态性的推波助澜地球物理反演里网格一加密病态性指数上升。原因很直观观测点数量有限地下网格越多方程组越欠定相邻网格对地表观测的响应差异又随深度增大而减小。A 的条件数随网格数的平方量级增长SVD 谱末端的奇异值会被推向机器精度。这时候用满秩 SVD 已经没有意义数值上也都是零空间噪声。所以后面所有正则化方法都基于“保留前 k 个奇异值方向”或“用 λ 压缩尾部”两种思路而不是去精确求逆。下面这张表是我看 Picard 图时常用的判读指标。曲线正常形态异常形态s_i从大到小近似几何衰减尾部水平趋近机器精度|u_i^T b|缓慢下降高频处仍不见收敛|u_i^T b| / s_i先降后升持续上升或全程平缓如果第三条曲线全程平缓说明噪声被平均摊到所有分量里反演结果对 λ 极其敏感需要额外加先验约束。3. Tikhonov 正则化的两种参数选择L 曲线与 GCV 实测对比3.1 从滤波因子看 λ 的语义Tikhonov 的解可以写成 SVD 域内的显式形式x_λ Σ (s_i^2 / (s_i^2 λ^2)) · (u_i^T b / s_i) · v_is_i^2 / (s_i^2 λ^2) 是 Tikhonov 的滤波因子。λ 越大越小奇异值方向被压得越多解越光滑λ 越小越接近最小二乘解。关键是 λ 的量纲跟数据幅度、奇异值尺度强相关不存在一个统一默认值。这就是为什么每次反演都必须单独做参数选择而不是从别的工区复制一个 λ 过来。3.2 l_curve.m 找拐角残差范数与解范数的平衡L 曲线把 log‖Ax-b‖ 作为横轴、log‖x‖ 作为纵轴扫过一串 λ 画成一条 L 形曲线。拐角处对应一个折中点再增加 λ解范数降得慢但残差涨得快再减小 λ残差降得慢但解范数涨得快。工具箱的 l_curve.m 不仅画图还返回数值的拐角位置 reg_corner。% 对 2.2 节生成的带噪右端项做 Tikhonov 正则化 [U, s, V] csvd(A); % 用 L 曲线自动选择正则化参数Tikh 表示按 Tikhonov 扫描 [reg_corner, rho, eta, reg_param] l_curve(U, s, bn, Tikh); % 用拐角处的 lambda 解出模型 x_lambda tikhonov(U, s, V, bn, reg_corner); % 对比恢复结果与真实模型 figure; plot(x, k-, LineWidth, 1.5); hold on; plot(x_lambda, r--, LineWidth, 1.2); legend(真实模型, Tikhonov解); title([L曲线选参: \lambda , num2str(reg_corner)]);l_curve 的第四个入参 Tikh 告诉它用 Tikhonov 的残差-范数对如果改成 tsvd 或 mtsvd它会按截断类方法重算曲线。rho 是残差二范数eta 是解的二范数reg_param 是对应扫过的 λ 序列。reg_corner 通常落在噪声水平对应的奇异值附近这符合一般认知数据噪声决定你可信的分量个数。我习惯把 reg_param 里拐角前后各 5 个点对应的解都画出来看解对 λ 的敏感度这在展示反演可靠性时很有说服力。3.3 gcv.m 统计视角不依赖噪声范数GCV广义交叉验证的思想是把某条观测暂时剔除用剩余数据预测它选取使预测误差最小的 λ。gcv.m 的调用方式% 用 GCV 选择 Tikhonov 正则化参数 [reg_gcv, G, reg_param_gcv] gcv(U, s, bn, Tikh); % 用 GCV 给出的 lambda 求解 x_gcv tikhonov(U, s, V, bn, reg_gcv);G 是每个候选 λ 的 GCV 函数值画出来是一条平滑单谷曲线谷底就是 reg_gcv。和 L 曲线比GCV 不需要估计噪声范数实现简单但实际数据上它偏保守——选出的 λ 常常偏小解比 L 曲线更毛糙。原因是 GCV 在数值上对过拟合的惩罚不够强。我用两个指标一起看L 曲线给一个 λGCV 给一个 λ两个解对比残差和解范数差异在 10% 以内说明选参可靠差太多就要怀疑数据里还有系统误差。两个方法的侧重点对比如下。选参方式需要的输入偏差倾向L 曲线SVD 分解拐角不明显时估值偏大GCVSVD 分解无需噪声范数噪声大时估值偏小3.4 参数选择前的两个预处理陷阱第一b 必须做与 A 一致的单位归一化。l_curve 的横纵轴双对数尺度对量纲很敏感如果 b 整体放大了 1000 倍reg_corner 会移动但解不变——这种位移容易让人误判。第二不要把 reg_param 里的 λ 直接拿去给下一个数据集用。λ 对 A 的奇异值谱非常敏感换网格、换观测布局都要重算。我见过不少反演程序把 λ 写成常量换一个工区就出问题这需要在参数选择函数里把 λ 设计成可计算的量而不是写死一个数。4. TSVD 与 MTSVD用截断代替光滑惩罚4.1 截断法在正则化反演中的位置TSVD 把 SVD 展开中高于 k 的分量直接丢弃x_k Σ_{i1}^{k} (u_i^T b / s_i) · v_ik 是截断参数等价于一个硬门槛滤波因子。Tikhonov 是软压缩TSVD 是硬切除。对大多数物探反演问题两者解的形状接近差别主要出现在尾部奇异值对应的模型细节Tikhonov 保留了一部分被压缩的高频细节TSVD 则完全看不到。如果你的目标模型本身就简单平滑TSVD 往往更干净。4.2 tsvd.m 与 mtsvd.m 的差异点TSVD 的调用极简用第 2 节得到的 SVD 分解和 L 曲线给出的截断指标就能算% 用 L 曲线按 TSVD 方式找最优截断点 [reg_corner, rho, eta, reg_param] l_curve(U, s, bn, tsvd); k_tsvd round(reg_corner); % 截断点取整 % 执行截断奇异值分解 [x_tsvd, rho_tsvd, eta_tsvd] tsvd(U, s, V, bn, k_tsvd);l_curve 返回的 reg_corner 对 TSVD 扫描的是整数 k这里 round 一下即可。剩余分量全部归零所以解范数比 Tikhonov 小光滑度都由前 k 个分量决定。MTSVDmtsvd.m是修正思路它同样取前 k 个主分量但剩下的分量不是直接丢掉而是通过一个模型约束空间 L 找一个在约束下残差最小的解。正则化工具箱里 mtsvd 的入参除了 U、s、V、b、k还接受一个 L。如果 L 取单位阵MTSVD 退化成带阻尼的 TSVDL 取二阶差分算子解会被约束得更平滑。这对物探解释有意义重力反演的密度分界可以把 L 设计成界面深度的平滑算子磁化率反演则更常取单位阵。方法滤波因子参数物探适用场景Tikhonovs_i^2 / (s_i^2 λ^2)λ 连续大多数重磁反演默认起点TSVDi≤k 为 1其余为 0k 整数对解范数敏感、需要紧凑模型MTSVD截断后按 L 约束修正k 与 L 矩阵带先验几何约束的物性分界4.3 截断点 k 的物理含义k 的物理含义值得多说一句它大致对应奇异值衰减到噪声平台的位置。shaw(128) 这个例子里奇异值从 10^0 到 10^-16 基本呈几何衰减加 10^-3 噪声时噪声主导起点大约在 s_i ≈ 10^-3 附近对应 k 通常在 20 到 30 之间。取大了放进噪声分量取小了丢失深部信息——深部网格在 SVD 谱里恰好在靠后的分量上这就是物探反演里“深部分辨率与稳定性不可兼得”的数学表达。5. 网格规模上来之后迭代正则化与 Lanczos 降维5.1 为什么二维三维网格下直接 SVD 会失去意义前几章示例都是 n128 的一维问题csvd 很快。一旦进入二维网格比如 100×100 的密度模型A 就是 10000×10000 的量级直接 SVD 的 O(n^3) 完全扛不住。而且正演算子通常以函数或稀疏矩阵存在根本不会显式构造全体 A。迭代法在这时替代 SVD 成为主选。工具箱里 cgls.m、pcgls.m、plsqr_b.m、bidiag.m、lanc_b.m 都是这组场景的工具它们不计算完整 SVD而是在 Krylov 子空间里逐步逼近。5.2 迭代次数就是正则化参数cgls 与 plsqr_b 实战CGLS共轭梯度最小二乘的迭代过程从零模型开始每步在一个 Krylov 子空间内降低残差。迭代步数少时解偏向光滑低频继续迭代会一步步把对应小奇异值的高频成分引进来。所以 CGLS 里迭代次数就等价于前面的 k 或 λ。下面是标准用法% Afun 是正演算子函数句柄接收模型向量 x返回观测值 % b_obs 是实测异常数据 k_iter 30; % 30 步共轭梯度最小二乘 [X_cgls, rho_cgls, eta_cgls] cgls(Afun, b_obs, k_iter);如果 A 是显式矩阵直接传 A 也可以Afun 返回 x 序列时则传函数句柄。rho_cgls 是每一步的残差范数eta_cgls 是解范数用这两个向量就能画出迭代版的 L 曲线拐角对应的迭代次数就是合适的截止点。preconditioned 版本是 pcgls.mP 需要你提供一个逼近 A 的预条件子把谱聚拢减少迭代步数。plsqr_b.m 则是带预条件的 LSQR 变体对大网格投影更稳% plsqr_b 三步配置50 次迭代返回解列矩阵与残差、范数序列 [X, rho, eta] plsqr_b(Afun, b_obs, 50);这里的 X 是每一迭代步的解拼成的列矩阵rho 和 eta 是每一步的残差范数和解范数。和 cgls 一样从 rho 的变化能看出迭代是否进入噪声拟合阶段。5.3 混合策略Lanczos 双对角化加小规模 Tikhonov实际项目里我常用混合方法先用 lanc_b.m 对 A 做 k 步 Lanczos 双对角化得到一个小型的双对角矩阵 B再在 B 上做正则化。B 的规模只有 k×k工具箱里所有参数选择函数都能直接上等于把大问题压成小问题再选 λ避免了大矩阵 SVD 的不可行。lanc_b 返回的 U、B、V 把原问题投影到了 k 维基空间投影后的右端项也参与后续的 l_curve 或 gcv。这样处理后即便 A 的维度到了 10^5 量级整个正则化计算也只花在 60×60 的小矩阵上。如果还想更快plsqr_b.m 可以在双对角化过程中同时迭代求解省掉显式回转。如果你后续打算做深度学习方向这套迭代生成的模型样本也可以当合成数据来用生成不同信噪比的“观测到模型”训练对。6. 测试问题集与验收流程上真实数据前先跑这一套6.1 标准测试问题的物理背景对照工具箱内置的 shaw、phillips、baart、foxgood、heat、blur、ilaplace、ursell 这组测试问题覆盖了图像重建、核密度估计、热传导逆问题、拉普拉斯逆变换等典型不适定场景。用它们不是因为“官方推荐”而是它们有解析真解可以定量评估正则化参数选得好不好。我在新工区反演前都会先在这组问题上跑一遍准备用的算法验证代码链路没有 bug再替换真实算子。下面是常用问题的背景和病态程度。函数原始物理背景病态程度shaw一维图像重建积分方程严重phillips核密度估计中等baart含对数核的第一类方程严重foxgood强振荡核积分方程非常严重heat热传导时间反演严重blur图像去模糊中等ilaplace拉普拉斯逆变换严重ursell无界核积分方程极度病态6.2 regudemo.m 一键对比各方法表现工具箱自带的 regudemo.m 无参数运行会把 shaw、phillips、baart、foxgood、heat 这几个问题依次用 Tikhonov、TSVD 各做一遍并画出对应的 L 曲线、GCV 曲线和解对比图。第一次打开这个工具箱的人我建议先跑它。regudemo;它能在半小时内让你建立对“λ 和 k 在什么量级上是合理的”的直觉判断也能顺便验证你的 MATLAB 环境里所有 m 文件路径是否都配好了。6.3 把测试问题结论迁移到实测数据的三个检查点换真实数据时不要只换 A 和 b还要做三个检查。第一对实测 b 跑 picard 诊断确认数据中存在可反演分量而不是纯噪声第二同时跑 l_curve 与 gcv 各给一个 λ对比两个解的残差范数与解范数若差异在可接受范围说明选参可靠若两者相差很大优先查系统误差而非继续调参数第三把反演结果正演回代计算回代残差与数据噪声水平。残差比噪声水平高出一个量级说明模型过于光滑欠拟合残差接近机器精度说明过拟合了噪声——合理目标是把残差控制在噪声水平的 1.1 到 1.3 倍。本文还有配套的精品资源点击获取