第4章 经典变分模型：Tikhonov、TV、小波与非局部先验

本章核心：所有的变分图像处理问题本质上都是在寻找一个平衡点——在相信观测数据（Data Fidelity）与满足先验假设（Regularization）之间博弈。本章不讨论“怎么解”（那是后续章节算法的任务），而是专注于“怎么算好解”——即能量泛函（Energy Functional）的设计。

4.1 开篇：病态问题与 MAP 视角

在第1章中，我们建立了通用成像模型 $y = A(x) + n$。 如果我们试图直接求逆 $x_{est} = A^{-1}y$，哪怕 $A$ 是可逆的，结果通常也是灾难性的。这是因为成像逆问题大多是 病态的（Ill-posed）。

4.1.1 为什么直接求逆会失败？

大多数成像算子 $A$（如模糊、降采样）本质上是低通滤波器。它们极大地衰减了高频信息。

  • 当我们试图求逆（$A^{-1}$）时，本质上是在进行高通滤波。
  • 观测噪声 $n$ 通常包含大量高频分量。
  • $A^{-1}$ 会将噪声的频率成分放大几个数量级，导致恢复出的图像被杂乱的噪声淹没，原本的图像结构荡然无存。

4.1.2 贝叶斯视角（MAP Estimation）

正则化并非凭空捏造，它在统计学上有严格的解释。最大后验概率估计（MAP）告诉我们： $$ P(x|y) \propto P(y|x) \cdot P(x) $$ 取负对数后，我们得到变分问题的标准形式： $$ \min_x \underbrace{-\log P(y|x)}_{\text{数据保真项 } D(Ax, y)} + \underbrace{-\log P(x)}_{\text{正则项 } \lambda R(x)} $$

  • 数据项：描述噪声分布（如高斯噪声对应 $L_2$ 距离）。
  • 正则项：描述我们要恢复的图像 $x$ 的先验概率。例如，如果我们认为“自然图像应当是平滑的”，就对应了特定的 $R(x)$。

4.2 二次正则与 Tikhonov：平滑的代价

最经典、数学上最易处理的正则化是 Tikhonov 正则化（在统计学中称为岭回归 Ridge Regression）。

4.2.1 模型形式

$$ E(x) = \frac{1}{2}|Ax - y|_2^2 + \frac{\lambda}{2} |\nabla x|_2^2 $$ 或者更一般的形式 $|Lx|_2^2$，其中 $L$ 是 Tikhonov 矩阵（通常是单位阵或差分矩阵）。

4.2.2 物理与几何直觉

最小化 $|\nabla x|_2^2$ 等价于假设图像是布朗运动的结果，或者说图像平滑度服从高斯分布。 从欧拉-拉格朗日方程（Euler-Lagrange Equation）来看，其梯度流对应于热传导方程（Heat Equation）： $$ \frac{\partial x}{\partial t} = \Delta x \quad (\text{扩散过程}) $$ 这意味着 Tikhonov 正则化会像热量扩散一样，让图像的高亮和阴影区域相互中和。

4.2.3 优点与致命缺陷

  • 优点：
    • 解析解：如果 $A$ 是线性的，这只是一个最小二乘问题，可以通过求解线性方程组 $(A^TA + \lambda \nabla^T\nabla)x = A^Ty$ 直接得到答案（闭式解）。
    • 计算极快：利用 FFT（如果 $A$ 是卷积）可以瞬间求解。
  • 缺陷：
    • 各向同性扩散：它不管你是边缘还是平坦区域，一视同仁地进行平滑。
    • 边缘模糊：图像中最重要的高频信息（边缘）被惩罚得最重（因为边缘处梯度很大，平方后更大）。

Rule of Thumb：永远不要在追求视觉质量的自然图像恢复中单独使用 Tikhonov。它只适用于非常具体的科学成像（如某些对平滑度有物理要求的流体场恢复）。

4.3 全变分（TV）：保边利器

1992年，Rudin, Osher 和 Fatemi 提出了著名的 ROF 模型，引入了全变分（Total Variation）正则化。这是现代图像处理的基石。

4.3.1 ROF 模型

$$ \min_x \frac{1}{2}|x - y|_2^2 + \lambda |\nabla x|_1 $$ 关键变化：梯度的 $L_1$ 范数 替代了 $L_2$ 范数。

4.3.2 几何直觉：Coarea Formula 与等高线

为什么 $L_1$ 能保边？除了“稀疏性”解释外，几何解释更为深刻。 根据 Coarea Formula，一个函数的 TV 范数等于其所有等高线（Level Sets）的长度之和： $$ TV(x) = \int_{-\infty}^{+\infty} \text{Length}(\{z | x(z) = c\}) dc $$ 最小化 TV，本质上是在缩短图像中所有物体的边界长度。

  • 噪声通常表现为细碎、周长很长的小斑点。TV 会迅速抹除这些小斑点（因为这能大幅减少周长）。
  • 真实的边缘是长且连贯的。TV 会倾向于把弯曲的边缘拉直，但不会模糊它（不会像热扩散那样把陡坡变成缓坡）。
ASCII 效果对比：一维信号

原始 (Step)    含噪观测        Tikhonov (L2)        TV (L1)
   +-----+       . : .            /                 +-----+
   |     |     : .   : .         /                  |     |
   |     |       . :            /                   |     |
---+     +---  ...     ...   ---+      +---      ---+     +---
                               (边缘变斜坡)        (边缘保持陡峭)

4.3.3 各向同性 (Iso) vs 各向异性 (Aniso)

这是工程实现中最容易混淆的细节。

  1. 各向同性 TV (Iso-TV): $$ |\nabla x|_{1,2} = \sum_i \sqrt{(\nabla_h x)_i^2 + (\nabla_v x)_i^2} $$

  • 几何：惩罚梯度的欧几里得长度。具有旋转不变性。

    • 效果：倾向于产生圆滑的拐角。 2. 各向异性 TV (Aniso-TV): $$ |\nabla x|_{1,1} = \sum_i (|\nabla_h x|_i + |\nabla_v x|_i) $$

  • 几何：惩罚梯度的曼哈顿长度（L1 of L1）。不具备旋转不变性。

    • 效果：倾向于产生水平和垂直的边缘。如果你处理的是草图、二维码或建筑物，Aniso-TV 甚至可能比 Iso-TV 更好；但对于人脸或自然风景，它会产生块状伪影。

4.3.4 阶梯效应 (Staircasing Effect)

TV 的阿喀琉斯之踵。由于 TV 的解倾向于分段常数（Piecewise Constant），当它试图逼近一个平滑的渐变区域（如天空、皮肤）时，会产生阶梯状的色块。

  • 数学原因：对于线性斜坡数据，TV 正则化的解要么是该斜坡本身（$\lambda$ 很小），要么直接变成平坦（$\lambda$ 很大），中间状态往往不稳定并退化为阶梯。

4.4 改进 TV：高阶与结构化正则

为了保留 TV 的优点（保边）并去除缺点（阶梯），数学家们进行了大量修正。

4.4.1 Huber-TV：可微性的妥协

为了解决 TV 在 0 点不可微导致的数值困难，以及缓解阶梯效应，常用 Huber 函数近似 $L_1$： $$ H_\epsilon(z) = \begin{cases} \frac{1}{2\epsilon}z^2 & |z| \le \epsilon \\ |z| - \frac{\epsilon}{2} & |z| > \epsilon \end{cases} $$

  • 小梯度处（平坦区）：表现为 $L_2$（Tikhonov），避免阶梯，产生平滑渐变。
  • 大梯度处（边缘）：表现为 $L_1$（TV），保持边缘锐利。

    Gotcha：$\epsilon$ 参数不仅影响算法收敛速度，也直接决定了平滑与保边的阈值。

4.4.2 TGV (Total Generalized Variation)

TGV 是目前变分图像恢复中公认效果最好的凸正则项之一。 它假设图像 $x$ 近似于一个函数 $v$ 的积分。模型形式（二阶 TGV）： $$ \min_{x, v} \alpha_1 |\nabla x - v|_1 + \alpha_0 |\mathcal{E}(v)|_1 $$

  • 直觉：引入了一个辅助向量场 $v$ 来吸收图像的导数。
    • 如果图像是线性斜坡，$\nabla x$ 是常数。我们可以让 $v = \nabla x$，第一项为0。第二项惩罚 $v$ 的导数（即 $x$ 的二阶导），也是0。
    • 结论：TGV 对多项式曲面（如斜坡）的惩罚为0！它能完美重建斜坡而不产生阶梯，同时保留边缘。

4.5 小波与 Frame 正则：频域的稀疏

TV 是在梯度域（空间域的差分）找稀疏性。小波（Wavelet）则是在变换域找稀疏性。

4.5.1 紧框架（Tight Frame）与过完备性

现代图像处理很少使用正交小波基（Orthogonal Basis），因为它们缺乏平移不变性（平移图像会导致系数剧变，产生伪影）。 常用的是过完备字典/紧框架（Over-complete Frame/Dictionary），如 Curvelet, Contourlet, Undecimated Wavelets (UWT)。

  • 特点：系数比像素多得多（冗余），这种冗余带来了对平移和纹理更好的鲁棒性。

4.5.2 分析模型 vs 合成模型

这是稀疏建模中最重要的两个流派，设 $W$ 为变换算子（如小波变换）。

  1. 分析模型 (Analysis Model / Cosparse): $$ \min_x \frac{1}{2}|Ax - y|_2^2 + \lambda |Wx|_1 $$
  • 逻辑：把图像 $x$ 扔进变换域，强迫系数稀疏。
    • 几何：解位于相空间中多个超平面的交集。
    • 优势：结果通常更自然，伪影较少，适合恢复任务。
  1. 合成模型 (Synthesis Model / Sparse Coding): $$ \min_\alpha \frac{1}{2}|A W^T \alpha - y|_2^2 + \lambda |\alpha|_1, \quad \text{令 } x = W^T \alpha $$
  • 逻辑：寻找一组稀疏系数 $\alpha$，用它们组装出图像。
    • 优势：模型简单，可用 ISTA/OMP 等算法求解。
    • 劣势：如果字典 $W$ 不够好，重建出的图像会有明显的“基函数长相”（如小波的振铃伪影）。

Rule of Thumb：如果是去噪或去模糊，首选分析模型（或 Analysis-Synthesis 混合）。如果是做图像压缩或特征提取，选合成模型。

4.6 非局部正则（Non-local Means / NLTV）

TV 和小波都是局部的（只看邻域）。但自然图像具有自相似性（Self-similarity）：眼睛的纹理在左眼和右眼都出现，甚至头发的纹理在全图随处可见。

4.6.1 图拉普拉斯视角

非局部正则化通过构建一个全连接（或大范围连接）的加权图来利用这一特性。

  • 对于像素 $i$ 和 $j$，计算它们周围 Patch 的相似度权重 $w_{ij}$。
  • 变分 NLTV： $$ J_{NL}(x) = \sum_{i,j} w_{ij} |x_i - x_j| $$ 这实际上是定义在加权图上的 TV。如果 $i$ 和 $j$ 的 Patch 很像（$w_{ij}$ 大），我们就强迫 $x_i$ 和 $x_j$ 的值也要接近。

4.6.2 效果与代价

  • 效果：这是去噪领域的“核武器”。对于周期性纹理、重复图案，NLTV 的效果远超 TV 和小波。它能从乱码中恢复出整齐的栅格。
  • 代价：计算 $w_{ij}$ 需要全图（或大搜索窗）搜索，复杂度极高。通常作为后处理或精细调节步骤。

4.7 数据项的鲁棒化：应对非高斯噪声

正则项决定了图像长什么样，数据项（Data Term） 决定了我们如何容忍误差。

| 数据项形式 | 概率解释 | 物理场景 | 优化特性 |

数据项形式 概率解释 物理场景 优化特性
$L_2$: $\frac{1}{2}|Ax-y|_2^2$ 高斯噪声 $\mathcal{N}(0, \sigma^2)$ 通用传感器热噪声 光滑凸，梯度是线性的 $(Ax-y)$
$L_1$: $|Ax-y|_1$ 拉普拉斯/脉冲噪声 传输错误、椒盐噪声、死像素 非光滑凸，需用 Prox
Huber / Charbonnier 混合分布 既有高斯底噪又有偶尔的大离群值 光滑凸，鲁棒性好
Kullback-Leibler (KL) 泊松分布 (Poisson) 光子计数成像（天文、显微、低光照、医学PET/CT） 凸，但非Lipschitz梯度

关键直觉：$L_2$ 数据项对大误差极其敏感（因为是平方）。如果图像中有遮挡（Inpainting）或死像素，使用 $L_2$ 会导致这些错误被“平均”到周围区域，形成模糊的光晕。此时必须用 $L_1$ 数据项，它能自动“忽略”离群点。

4.8 本章小结：模型选择指南

| 模型 | 核心特征 | 优点 | 缺点 | 适用场景 |

模型 核心特征 优点 缺点 适用场景
Tikhonov ($L_2$) 惩罚梯度平方 极快，闭式解 边缘模糊 光流平滑项、曲面拟合
TV ($L_1$) 惩罚梯度模长 保边，对比度不变 阶梯效应，纹理丢失 工业检测、卡通图、医学分割预处理
TGV 一阶+二阶导联合 保边且保斜坡 计算稍慢 高质量自然图像恢复
Wavelet ($L_1$) 频域稀疏 保留纹理细节 振铃伪影 (Gibbs) 纹理丰富的图像去噪
NLTV 非局部 Patch 相似 纹理恢复极强 极慢 极低信噪比、重复纹理修复

4.9 练习题

基础题

  1. 梯度推导：考虑正则项 $R(x) = \frac{1}{2}|\nabla x|_2^2$。请写出其在离散域的梯度表达式。假设使用周期边界条件。

    Hint 回忆拉普拉斯算子 $\Delta$ 与 $\nabla$ 的关系。
    答案 $\nabla R(x) = -\Delta x$（负拉普拉斯算子）。在离散域对应卷积核 $\begin{bmatrix} 0 & -1 & 0 \\ -1 & 4 & -1 \\ 0 & -1 & 0 \end{bmatrix}$。

  2. 各向异性计算：给定 3x3 图像块 $\begin{bmatrix} 1 & 1 & 1 \\ 5 & 5 & 5 \\ 1 & 1 & 1 \end{bmatrix}$。计算中心像素的 Aniso-TV 和 Iso-TV贡献（仅考虑中心点的前向差分）。

    Hint $\nabla_h$: 右边减自己; $\nabla_v$: 下边减自己。

  3. L1 vs L2 数据项：给定一组观测值 $y = [1, 1, 1, 100, 1]$（明显有一个离群值 100）。 (a) 求 $x$ 最小化 $\sum (x - y_i)^2$。 (b) 求 $x$ 最小化 $\sum |x - y_i|$。 哪个结果更能代表数据的“真实”水平？

    答案 (a) $x = \text{mean}(y) = 20.8$。 (b) $x = \text{median}(y) = 1$。 显然 (b) 对离群值鲁棒。

挑战题 / 思考题

  1. TGV 的零空间：二阶 TGV 的能量泛函对什么样的图像 $x$ 值为 0？这解释了它为什么能避免阶梯效应？

    Hint 考虑 $x$ 为仿射函数（平面）$x(i,j) = ai + bj + c$ 时，$\nabla x$ 是什么？能否找到对应的 $v$ 使各项均为 0？

  2. 分析与合成的等价性：证明当变换矩阵 $W$ 是正交阵（Orthonormal）时，分析模型与合成模型是完全等价的。如果 $W$ 是冗余的（Over-complete），这种等价性还成立吗？

    Hint 正交阵满足 $W^T W = I$。代换变量 $x = W^T \alpha \Rightarrow \alpha = Wx$。对于冗余字典，行数 > 列数，逆变换不唯一，模型不再等价。

  3. 尺度不变性：经典的 TV 模型 $\min \frac{1}{2}|x-y|^2 + \lambda |\nabla x|_1$ 是对比度不变（Contrast Invariant）的吗？即如果 $y \to cy$，最优解 $x^*$ 也会变成 $cx^*$ 吗？如果不是，如何修改模型使其具有这一性质？

    Hint 检查数据项和正则项的齐次性（Homogeneity）。二次项是 2 次齐次，L1 是 1 次齐次。不匹配。这也是为什么 $\lambda$ 难调的原因。

4.10 常见陷阱与错误 (Gotchas)

  1. 离散化梯度的方向错位

    • 错误：实现 TV 时，梯度算子 $\nabla$ 用了前向差分，而散度算子（$\nabla$ 的转置）也用了前向差分，或者简单地用 imgradient 这种经过平滑的算子。
    • 后果：优化无法收敛，或者结果出现奇怪的相移。
    • 正解：必须严格保证伴随关系（Adjoint Property）。如果 $\nabla$ 是前向差分（$x_{i+1}-x_i$），那么 $\nabla^*$ 必须是后向差分（$x_{i}-x_{i-1}$）的负数。可以在 Chapter 17 找到代码验证方法。

  2. 边界条件（Boundary Conditions）

    • 现象：恢复出的图像四周有一圈黑框或奇怪的波纹。
    • 原因：FFT 实现卷积时默认是周期边界（Circular）。图像左边的边缘会和右边相互作用。
    • 对策：在进行计算前对图像进行 pad（如对称填充），或者使用 Neumann 边界条件（反射边界）推导差分矩阵。

  3. $\lambda$ 的数值敏感性

    • 陷阱：直接对 0-255 的图像进行优化。
    • 后果：数据项 $|x-y|^2$ 的值会非常大（$255^2 \approx 65000$），导致你需要设置 $\lambda$ 为 0.0001 甚至更小，不仅难调，还容易遇到浮点精度问题。
    • 建议：始终将图像归一化到 [0, 1] 区间处理。此时 $\lambda$ 通常在 $10^{-3}$ 到 $10^{-1}$ 之间。

  4. 把 Prox 算子当梯度用

    • 错误：在使用梯度下降法求解 TV 时，试图直接对 $|\nabla x|_1$ 求导。
    • 问题：$L_1$ 在 0 处不可导。虽然可以用次梯度，但直接梯度下降会震荡。
    • 对策：对于非光滑项（TV, L1），必须使用 近端梯度法（Proximal Gradient）、ADMM 或 Primal-Dual 算法，或者使用 Huber 平滑近似后再求导。不要强行求导。