2. 反问题数学基础

反问题建模的第一步,不是急着写代码或套优化器,而是把“未知量—观测—映射”三者的数学关系摆清楚。本节的目标很明确:建立一套能支撑后续不适定性分析的线性代数语言,核心工具是向量空间、矩阵分解、奇异值分解(SVD)与条件数。理解这些工具,才能在看到一条反演结果时判断:它究竟是数据本身给出的,还是算法在放大噪声。

从算子方程看反问题的统一形式

正问题通常写成“已知模型参数 m,计算观测数据 d”,即 d = F(m)。反问题则是给定 d 求 m,记作 F(m) = d。当 F 为线性算子时,离散化后得到矩阵方程:

A x = b,   A ∈ R^(m×n),  x ∈ R^n,  b ∈ R^m

这里 x 是待辨识的参数向量(如介质电导率、源强分布、弹性参数),b 是观测数据向量,A 是由物理模型和离散格式决定的敏感度矩阵(雅可比矩阵的线性化形式)。反问题的不适定性,几乎全部藏在这个 A 的性质里:它是否满秩、奇异值衰减多快、条件数多大。因此后续所有分析,本质上都是对 A 的“体检”。

向量空间:给未知量和观测一个几何框架

把 x 看作参数空间 R^n 中的向量,b 看作数据空间 R^m 中的向量,A 就是从参数空间到数据空间的线性映射。这个视角带来三个直接有用的概念:

  • 值域 Range(A):所有可能被模型产生的数据集合。若真实数据 b 含有模型无法表达的成分(如测量噪声、未建模物理),这部分就落在值域之外,反演只能求其投影。
  • 零空间 Null(A):满足 A x = 0 的参数扰动方向。沿零空间改变参数,观测数据完全不变——这意味着这些参数分量根本不可辨识,任何反演结果在这些方向上都是任意的。
  • 正交投影:最小二乘解 x̂ = argmin‖Ax−b‖² 的几何含义,是把 b 投影到 Range(A) 上,残差 r = b − Ax̂ 与值域正交。

实操中判断零空间是否“作祟”,可对 A 做数值秩估计:用 numpy.linalg.matrix_rank(A, tol=...) 或直接看 SVD 的奇异值谱。若最小奇异值比最大奇异值低十几个数量级,说明数值上已接近奇异,零空间方向开始主导解的不稳定性。

矩阵分解:把耦合的参数关系拆开

直接求解 A x = b 往往因参数间强耦合而失败。矩阵分解的作用是把 A 拆成若干结构简单的因子,使参数关系显式化。常用三种:

  • LU 分解:A = LU,适合方阵、良态问题的直接求解,计算量约 O(n³)。它不揭示不适定性,只负责高效解方程。
  • QR 分解:A = QR,Q 正交、R 上三角,是求解最小二乘的稳定手段,数值上比正规方程 AᵀA x = Aᵀb 更可靠,因为后者会把条件数平方。
  • 特征分解:仅对对称正定或可对角化方阵适用,A = VΛV⁻¹。在反问题中常用于正则化项(如 Tikhonov 的 LᵀL)的分析。

需要强调:LU、QR 解决的是“算得稳”,而 SVD 解决的是“看得清”。前者是计算工具,后者是诊断工具,二者在反问题流程中角色不同。

奇异值分解:不适定性分析的显微镜

对任意实矩阵 A ∈ R^(m×n),SVD 给出:

A = U Σ Vᵀ
U ∈ R^(m×m) 正交,  V ∈ R^(n×n) 正交
Σ = diag(σ₁ ≥ σ₂ ≥ … ≥ σ_r > 0),  r = rank(A)

把它代入 A x = b,得到解的谱展开形式:

x = Σᵢ (uᵢᵀ b / σᵢ) vᵢ

这个式子信息量极大。uᵢ 是数据空间的正交基,vᵢ 是参数空间的正交基,σᵢ 是第 i 个“模式”的增益。反演过程就是把数据投影到 uᵢ 上得到系数 uᵢᵀb,再除以 σᵢ 放大回参数空间。问题在于:当 σᵢ 很小时,数据中哪怕极微小的噪声分量(量级 ε)也会被放大成 ε/σᵢ 的参数误差。这就是不适定性的谱解释——小奇异值方向是噪声放大器。

一个可操作的诊断流程:对 A 做 SVD,画出 σᵢ 随 i 的衰减曲线(对数纵轴)。若曲线在某个 i 之后陡降数个数量级,说明有效秩远小于名义秩,必须截断或正则化。工程上常取“奇异值比 σᵢ/σ₁ 低于 10⁻⁶~10⁻⁸”作为数值零的判据,具体阈值依数据噪声水平而定。

条件数:一个数概括稳定性

条件数定义为 κ(A) = σ₁/σ_r(2-范数下)。它衡量解的相对误差对数据相对误差的放大倍数:

‖δx‖/‖x‖ ≤ κ(A) · ‖δb‖/‖b‖

经验判断:κ 在 10¹~10² 属良态;10³~10⁵ 需警惕;超过 10⁶ 基本可判定为病态,直接求解不可信。注意条件数与尺度有关——若参数单位不统一(如一个参数是 10⁻³、另一个是 10⁶),A 的列尺度差异会人为抬高 κ。因此建模时应先做列归一化(每列除以自身范数),再评估条件数,否则得到的是“单位病态”而非“物理病态”。

条件数只反映最坏情况放大,不反映噪声在哪些方向被放大。真正指导正则化设计的是完整奇异值谱,而非单一 κ 值。二者应配合使用:κ 做快速筛查,SVD 谱做精细诊断。

把工具串成一条分析链

下面这张图把本节的核心关系组织成一条从矩阵到诊断结论的分析链,重点表现 SVD 如何把“数据—参数”的耦合拆成独立模式,以及小奇异值如何成为噪声放大通道。

反问题线性代数分析链:从矩阵到不适定性诊断 算子方程 A x = b A: 敏感度矩阵 奇异值分解 A = U Σ Vᵀ σ₁ ≥ σ₂ ≥ … ≥ σᵣ > 0 r = rank(A) 解的谱展开 x = Σᵢ (uᵢᵀb / σᵢ) vᵢ 数据投影 → 除以 σᵢ → 参数分量 小 σᵢ 方向 = 噪声放大通道 零空间 Null(A) A x = 0 的参数扰动 该方向参数不可辨识 → 需先验约束或降维 奇异值谱 σᵢ 对数纵轴看衰减速度 陡降 → 有效秩 < 名义秩 → 截断 SVD / 正则化 条件数 κ = σ₁/σᵣ κ > 10⁶ 判为病态 先列归一化再评估 → 快速筛查稳定性 诊断结论 → 正则化设计依据 κ 做快速筛查,SVD 谱定截断阈值,零空间定不可辨识方向

一个可复现的数值示例

设某一维参数辨识问题离散后有 6 个参数、8 个观测,敏感度矩阵 A 由相邻参数对观测的加权叠加构成。以下为示例代码,用于观察奇异值谱与条件数(数值为合理假设,非实测):

import numpy as np
A = np.array([[1.0, 0.8, 0.0, 0.0, 0.0, 0.0],
              [0.8, 1.0, 0.8, 0.0, 0.0, 0.0],
              [0.0, 0.8, 1.0, 0.8, 0.0, 0.0],
              [0.0, 0.0, 0.8, 1.0, 0.8, 0.0],
              [0.0, 0.0, 0.0, 0.8, 1.0, 0.8],
              [0.0, 0.0, 0.0, 0.0, 0.8, 1.0],
              [0.5, 0.0, 0.0, 0.0, 0.0, 0.5],
              [0.0, 0.5, 0.0, 0.0, 0.5, 0.0]])
U, s, Vt = np.linalg.svd(A, full_matrices=False)
print("奇异值:", np.round(s, 4))
print("条件数:", round(s[0]/s[-1], 1))

这类相邻耦合矩阵的典型表现是:前几个奇异值接近 2,随后快速衰减,最小奇异值可能低至 10⁻² 量级,条件数达到数百。若在数据 b 上叠加 1% 的相对噪声,则沿最小奇异值方向的参数误差可能被放大到 10% 以上。此时若直接求逆,结果会出现沿参数序列的高频振荡——这正是“噪声被小奇异值方向放大”的可视化表现。处理办法不是换更快的求解器,而是截断小奇异值(TSVD)或加入 Tikhonov 正则项,把放大倍数压回可接受范围。

常见误区:看到残差 ‖Ax̂−b‖ 很小就认为反演成功。在病态问题中,过拟合噪声同样能产生极小残差,但参数已严重偏离真值。判断标准应是“残差合理 + 参数在物理约束内 + 对噪声扰动稳定”,三者缺一不可。

本节建立的这套语言——向量空间定位可辨识性、矩阵分解保证计算稳定、SVD 揭示噪声放大通道、条件数给出快速筛查指标——是后续所有正则化方法与不适定性分析的共同基础。拿到任何一个反问题,先对敏感度矩阵做一次 SVD 和条件数体检,再决定用什么求解策略,应当成为固定习惯。