8.2 8.2 不确定性量化与反问题

8.2 不确定性量化与反问题

本节摘要：正问题"由因求果"通常良态，反问题"由果溯因"往往病态——观测的微小噪声在反演中被放大成参数的巨大摆动。本节用一个热传导反演工单演示病态的来源，用 Tikhonov 正则化治标（给出稳定的点估计），再用贝叶斯反演治本（给出参数的后验分布），MCMC 采样落地实现。

本节能力目标

阅读完本节，你应当能够：

1. 解释反问题病态的几何原因（信息被正向过程压缩）；
2. 实现 Tikhonov 正则化并用 L 曲线选正则参数；
3. 写出贝叶斯反演的框架代码，用 MCMC 采后验。

工单：从边界温度反推热源

一根一维导热杆，内部藏着一个热源，我们只能在杆端测到若干温度读数，要反推热源的强度分布。正问题（给定热源算边界温度）是第 5 章的常规操作；反问题把箭头反过来，麻烦从此开始。

几何直觉理解病态：正问题像一个"压缩打包"过程——多种不同的热源分布可能产生几乎相同的边界温度（信息在正向传播中被平均掉了）。反演就是"解压缩"：要把被压掉的细节还原回来，而观测噪声让还原过程极度过敏。第 2 章的 SVD 语言在此直接复用：正问题算子的小奇异值方向就是"几乎不影响观测"的参数方向，反演时除以小奇异值，噪声被放大到天文数字。

import numpy as np

rng = np.random.default_rng(21)
n = 60

# 正问题算子：把热源分布映射为边界温度（平滑平均型算子）
x = np.linspace(0, 1, n)
A = np.exp(-((x[:, None] - x[None, :])**2) / 0.05) / np.sqrt(0.05*np.pi)

source_true = np.sin(2*np.pi*x)                    # 真实热源
obs_clean = A @ source_true
obs_noisy = obs_clean + 0.001 * rng.standard_normal(n)

# 朴素反演：直接除奇异值
U, s, Vt = np.linalg.svd(A, full_matrices=False)
coeff = (U.T @ obs_noisy) / s
naive = Vt.T @ coeff
print("奇异值跨度:", f"1e0 到 1e{s.mean():0f}（量级 {np.log10(s[0]/s[-1]):.0f} 个数量级）")
print(f"朴素反演最大偏差 {np.max(np.abs(naive - source_true)):.2f}")
print(f"观测噪声水平仅 1e-3，被放大了 {np.max(np.abs(naive - source_true))/0.001:.0f} 倍")
​

千分之一的观测噪声被放大成参数上的满幅振荡——病态不是数值bug，是问题本身缺信息。

Tikhonov 正则化：承认缺信息

治法与第 2 章的岭回归同构：在目标里加参数范数的平方惩罚，等价于把小奇异值方向的分量软性压扁而不是除以它。正则参数拉姆达控制"信观测多少"：零则回到病态，无穷则参数归零。选拉姆达的经典工具是 L 曲线：横轴解的范数、纵轴残差范数，对数坐标下常呈 L 形，拐点（曲率最大处）是"残差与解规模"的最佳折中。

import numpy as np

rng = np.random.default_rng(21)
n = 60
x = np.linspace(0, 1, n)
A = np.exp(-((x[:, None] - x[None, :])**2) / 0.05) / np.sqrt(0.05*np.pi)
source_true = np.sin(2*np.pi*x)
obs_noisy = A @ source_true + 0.001 * rng.standard_normal(n)
U, s, Vt = np.linalg.svd(A, full_matrices=False)

def tikhonov(lam):
    """谱域滤波形式的 Tikhonov 解"""
    f = s / (s**2 + lam**2)                       # 滤波因子
    return Vt.T @ (f * (U.T @ obs_noisy))

for lam in [1e-6, 1e-4, 1e-2, 1.0]:
    sol = tikhonov(lam)
    resid = np.linalg.norm(A @ sol - obs_noisy)
    size = np.linalg.norm(sol)
    err = np.max(np.abs(sol - source_true))
    print(f"lam={lam:.0e}: 残差 {resid:.4f}, 解范数 {size:.2f}, "
          f"真误差 {err:.3f}")
​

拉姆达扫过四个量级，误差先降后升——中间某处是甜点区。Tikhonov 给的是一个稳定的点估计，但它回避了一个问题："这个估计有多可信？"回答这个问题要把话事权交给贝叶斯。

贝叶斯反演：给参数发"可信度分布"

贝叶斯观点把未知参数当随机变量：先验分布编码"观测之前我们相信什么"（比如"热源不会太离谱"的高斯先验），似然函数编码"给定参数时观测出现的可能性"（噪声模型），后验分布由两者按贝叶斯公式合成。后验不是一个数，是一整张分布图——它的宽度就是不确定性的诚实度量。

后验通常没有解析形式，MCMC（马尔可夫链蒙特卡洛）是通用采样器：构造一条以后验为平稳分布的随机游走（梅特罗波利斯规则：提议新点，按后验概率之比决定去留），走足够多步后，链条上的样本分布近似后验。

import numpy as np

rng = np.random.default_rng(88)

def log_prior(source, tau=5.0):
    """高斯先验：热源幅度不太大"""
    return -0.5 * np.sum(source**2) / tau**2

def log_likelihood(source, sigma=0.001):
    """观测噪声模型：高斯似然"""
    r = A @ source - obs_noisy
    return -0.5 * np.sum((r / sigma)**2)

def log_posterior(source):
    return log_prior(source) + log_likelihood(source)

# 梅特罗波利斯采样（模式坐标批量提议, 演示用低维参数化）
def metropolis(n_samples=8000, step=0.02):
    theta = np.array([1.0, 0.5])              # 两参数的正弦形状系数
    basis = np.vstack([np.sin(2*np.pi*x), np.sin(4*np.pi*x)])
    def to_source(th): return basis.T @ th
    lp = log_posterior(to_source(theta))
    chain = [theta.copy()]
    accept = 0
    for _ in range(n_samples):
        prop = theta + step * rng.standard_normal(2)
        lp_new = log_posterior(to_source(prop))
        if np.log(rng.random()) < lp_new - lp:
            theta, lp = prop, lp_new
            accept += 1
        chain.append(theta.copy())
    return np.array(chain), accept / n_samples

chain, rate = metropolis()
burn = chain[len(chain)//2:]                    # 丢弃暂态半程
mean_est = burn.mean(axis=0)
ci = np.quantile(burn, [0.05, 0.95], axis=0)
print(f"接受率 {rate:.2f}（健康区间约 0.2-0.5）")
print(f"后验均值系数 {np.round(mean_est, 3)}（真值 [1.0, 0.0]）")
print(f"90% 可信区间: 第一个参数 [{ci[0][0]:.3f}, {ci[1][0]:.3f}]")
​

后验均值贴近真值，且第一个系数带着一个宽度有限的区间——**"热源主波强度是 0.98，九成把握在 0.94 到 1.02 之间"**这样的表述才是反演工单的合格交付物。MCMC 的工程注意点：暂态要丢弃（收敛前样本不代表后验）、接受率调到健康区间（步长太大几乎全拒、太小原地打转）、多链起点交叉验证收敛。这些诊断与第 4 章的验证纪律一脉相承。

⚠️ 常见坑：把 MCMC 的前几百个样本当后验用。链条需要"热身"，收敛诊断（多链对比、自相关时间）没过之前的样本会系统性偏向前期起点，结论整体带偏。

💡 关键直觉：正则化与贝叶斯在高斯场合数学等价——Tikhonov 的拉姆达对应先验精度，选拉姆达的 L 曲线对应证据最大化。两条路线殊途同归，但贝叶斯语言把"不确定性"从副产品升级为交付主体，这正是"可信建模"时代要的语法。

反问题治疗路线

本节要点回顾

• 反问题天然病态：正问题压缩信息，小奇异值方向反演放大噪声；
• Tikhonov 压扁小奇异方向，L 曲线拐点定正则强度，交付点估计；
• 贝叶斯把不确定性升级为交付主体，后验宽度即可信度；
• MCMC 三纪律：丢暂态、调接受率、多链验收敛；
• 正则化与贝叶斯在高斯场合等价，语言不同，药效相通。