第1章 导论：从成像模型到优化问题

1. 开篇：看不见的世界与数学的魔法

在计算机视觉与图像处理的浩瀚领域中，我们通常面临两类截然不同的任务：

  1. 高层视觉（High-level Vision / Analysis）：教计算机“看懂”图片。例如：这张图里有一只猫吗？它的位置在哪里？这通常属于判别模型（Discriminative Models）的范畴。
  2. 底层视觉（Low-level Vision / Restoration）：教计算机“看清”图片。例如：这张图充满了噪点，能把噪点去掉吗？这张图模糊了，能复原清晰吗？这属于生成模型或逆问题（Inverse Problems）的范畴。

本书专注于底层视觉。在这个领域，我们不仅是把图像当作像素矩阵，更是将其视为物理信号的离散采样。我们的目标是通过数学建模，逆转物理退化过程。

如果你想理解如何把模糊的监控画面变清晰，或者想知道 MRI（核磁共振）是如何从一堆毫无意义的频率数据中重建出大脑图像的，那么欢迎来到变分与优化的世界。在这里，我们将物理定律翻译成目标函数，将图像质量的追求转化为数学约束。

本章学习目标：

  • 深刻理解成像的前向模型 $y = A(x) + n$ 及其物理对应。
  • 掌握病态性（Ill-posedness）的本质，理解为什么“直接求逆”是死路一条。
  • 建立贝叶斯推断与变分正则化之间的桥梁。
  • 熟悉常用的先验假设（稀疏性、低秩性、平滑性）及其对应的数学范数。
  • 获得优化算法的宏观地图，知道何时该用梯度下降，何时该用近端算子。

2. 逆问题全景图：$y = A(x) + n$

几乎所有的图像复原问题，都可以被浓缩进一个极其简洁的线性方程组中：

$$ y = A(x) + n $$ 这个公式虽然简单，但它是连接现实世界与数学世界的桥梁。让我们逐一解剖其中的变量。

2.1 变量的解剖

  • $x \in \mathbb{R}^N$（真实图像）：这是我们要寻找的“圣杯”。在数学推导中，我们通常把一张 $H \times W$ 的二维图像拉直（Vectorize）为一个长度 $N = H \cdot W$ 的列向量。
    • 注意：虽然数学上写成向量，但在工程实现时，我们尽量保持其 $H \times W$ 的矩阵结构以利用空间局部性。
  • $y \in \mathbb{R}^M$（观测数据）：这是传感器实际吐给我们的数据。它可能是模糊的照片、缺失像素的图像，甚至是经过傅里叶变换后的频谱（如 MRI）。注意 $M$ 可能等于、小于甚至远小于 $N$。
  • $n \in \mathbb{R}^M$（噪声）：不可避免的随机扰动。它来源于光子散粒噪声、电路热噪声或传输误差。
  • $A \in \mathbb{R}^{M \times N}$（线性算子/成像矩阵）：这是物理退化过程的数学描述。

2.2 算子 $A$ 的七十二变

根据算子 $A$ 的不同构造，我们定义了不同的图像处理任务。理解 $A$ 的结构是建模的第一步。

| 任务 (Task) | 算子 $A$ 的性质 | 物理意义 | 矩阵形态 (直观) |

任务 (Task) 算子 $A$ 的性质 物理意义 矩阵形态 (直观)
去噪 (Denoising) 单位矩阵 $I$ 理想成像，仅受噪声干扰 对角线为1
去模糊 (Deblurring) 卷积矩阵 $H$ 镜头失焦、物体运动、大气湍流 托普利兹 (Toeplitz) 或循环矩阵
超分辨率 (Super-Resolution) 降采样 $D$ $\times$ 模糊 $H$ 传感器像素密度不足 扁矩阵 (Fat Matrix)，行数少
修复/补全 (Inpainting) 掩膜矩阵 $M$ 遮挡、坏点、文字水印覆盖 从 $I$ 中抽掉若干行
压缩感知 (Compressed Sensing) 随机投影 $\Phi$ 单像素相机、加速成像 随机高斯/伯努利矩阵
断层重建 (CT/MRI) 变换域采样 $\mathcal{F}$ 测量的是投影或频率 傅里叶/Radon 变换矩阵

Rule of Thumb (工程法则) 在实际代码中，永远不要显式生成矩阵 $A$！ 对于一张 $1024 \times 1024$ 的图像，$N \approx 10^6$。矩阵 $A$ 的大小是 $10^6 \times 10^6$，存储它需要数 TB 的内存。 做法：我们实现两个函数（或类方法）：

  1. A(x)：模拟前向过程。
  2. At(y)：模拟 $A^T$（伴随算子/转置矩阵）的过程。 所有的优化算法只需要这两个函数即可运行。

2.3 ASCII 可视化：从 $x$ 到 $y$

      [ 真实世界 x ]                 [ 算子 A 的作用 ]                 [ 观测 y ]

      +---------+                   (Convolution)                  +---------+
      |  Clear  |                     Shift & Sum                  |  Blurry |
      |  Image  |    -------->      Kernel  K        -------->     |  Image  |
      +---------+                                                  +---------+
          |                                                             ^

          |                                                             ^
          |                                                             |
          v                                                             |

      (Vectorized)                                                  (Add Noise)
      [ x_1 ]                                                           |
      [ x_2 ]       x      [      A (Huge Matrix)      ]      =     [ y_clean ]  +  n
      [ ... ]              [ (Implicitly structured)   ]
      [ x_N ]

3. 病态性（Ill-posedness）：为什么不能直接求 $x = A^{-1}y$？

既然 $y = Ax$，那为什么不直接计算 $x = A^{-1}y$ 呢？这引出了逆问题的核心困难——病态性。 法国数学家 Hadamard 定义了“适定问题”（Well-posed problem）必须满足三个条件：

  1. 存在性 (Existence)：解存在。
  2. 唯一性 (Uniqueness)：解是唯一的。
  3. 稳定性 (Stability)：解连续依赖于数据（数据的微小变化只会导致解的微小变化）。

图像逆问题通常违反这三条中的至少一条，尤其是第三条。

3.1 稳定性的灾难

假设 $A$ 可逆，但它通常是一个“平滑”算子（如模糊）。这意味着 $A$ 的奇异值（Singular Values）衰减非常快，包含很多接近于 0 的微小值。 根据 SVD 分析，$A^{-1}$ 的奇异值是 $A$ 的倒数。 $$ \hat{x} = A^{-1}(y_{true} + n) = x_{true} + \underbrace{A^{-1}n}_{\text{Explosion!}} $$ 如果 $A$ 有一个特征值 $\sigma_i = 10^{-6}$，那么对应的噪声分量将被放大 $10^6$ 倍！恢复出来的图像将被巨大的噪声淹没，完全不可用。

3.2 唯一性的缺失

在超分辨率或修复问题中，$A$ 是不可逆的“扁”矩阵（欠定系统，Under-determined）。这意味着有无穷多个 $x$ 都能满足 $Ax = y$。 例如，填补缺失的像素：我们可以填黑色，也可以填白色，方程都成立。但只有一种填法看起来像“真实的图像”。

4. 变分建模：贝叶斯视角的救赎

为了从无数可能的解中通过“选美”挑出最合理的一个，或者为了压制噪声的爆炸，我们需要引入正则化（Regularization）。

4.1 贝叶斯公式推导

我们利用贝叶斯定理： $$ P(x|y) = \frac{P(y|x) P(x)}{P(y)} $$ 我们想要最大化后验概率（MAP）： $$ \hat{x}_{MAP} = \arg\max_x P(y|x) P(x) $$ 对两边取负对数（Negative Log-Likelihood），问题转化为最小化能量函数： $$ \min_x \underbrace{-\log P(y|x)}_{\text{数据保真项 } D(Ax,y)} + \underbrace{-\log P(x)}_{\text{正则项 } \lambda R(x)} $$ 这就构成了本书通用的变分模型（Variational Model）： $$ \min_x \quad D(Ax, y) + \lambda R(x) $$ 其中 $\lambda$ 是平衡参数，控制了我们多大程度上相信观测数据，多大程度上相信先验知识。

5. 数据保真项 $D(Ax, y)$：噪声决定形式

这一项由噪声的物理统计特性决定。

5.1 高斯噪声 (Gaussian Noise) $\to$ $\ell_2$ 范数

如果 $n \sim \mathcal{N}(0, \sigma^2I)$，则 $P(y|x) \propto \exp(-|Ax-y|_2^2 / 2\sigma^2)$。 $$ D(Ax, y) = \frac{1}{2} |Ax - y|_2^2 $$

  • 特点：处处可导（光滑），计算梯度容易（$A^T(Ax-y)$）。
  • 缺点：对离群点（Outliers）极度敏感（平方惩罚）。如果图像中有坏点，$\ell_2$ 会导致整个图像模糊以“平均”这个错误。

5.2 脉冲/拉普拉斯噪声 (Impulse/Laplace) $\to$ $\ell_1$ 范数

如果噪声服从拉普拉斯分布（尾巴比高斯长），或者包含椒盐噪声。 $$ D(Ax, y) = |Ax - y|_1 = \sum_i |(Ax)_i - y_i| $$

  • 特点：鲁棒（Robust），能容忍极端错误。
  • 缺点：在 0 处不可导，优化难度大（需要近端梯度法或 ADMM）。

5.3 泊松噪声 (Poisson) $\to$ Kullback-Leibler 散度

常见于低光照成像、天文摄影、显微成像（光子计数过程）。 $$ D(Ax, y) = \langle Ax, 1 \rangle - \langle y, \log(Ax) \rangle $$

  • 特点：非对称。需要强制约束 $x \ge 0$。

6. 正则项 $R(x)$：先验决定画质

这一项决定了我们恢复出的图像长什么样。这是过去几十年图像处理研究的核心。

6.1 Tikhonov 正则化 ($\ell_2$)

$$ R(x) = | \Gamma x |_2^2 $$ 通常 $\Gamma$ 是梯度算子。

  • 假设：图像是光滑的，能量有限。
  • 效果：虽然抑制了噪声，但会平滑掉图像的边缘，看起来像失焦了一样。现在很少单独使用。

6.2 全变分 (Total Variation, TV)

$$ R(x) = | \nabla x |_1 = \sum_{i,j} \sqrt{(\partial_x x)_{ij}^2 + (\partial_y x)_{ij}^2} $$

  • 假设：图像由分段常数区域组成（Piecewise Constant）。梯度是稀疏的（大部分地方梯度为0，边缘处梯度大）。
  • 效果：这是图像处理的里程碑。它能完美地保留边缘的同时去除噪声。
  • 副作用：阶梯效应（Staircasing），平滑的渐变区域会被变成像梯田一样的色块。

6.3 稀疏表示 (Sparsity)

$$ R(x) = | \Psi x |_1 $$ 其中 $\Psi$ 是某种变换（小波 Wavelet、DCT、Curvelet）。

  • 假设：自然图像在变换域是稀疏的（大部分系数为 0）。
  • 效果：比 TV 更精细，能保留纹理。

6.4 低秩 (Low-Rank)

$$ R(X) = | X |_* = \sum \sigma_i(X) $$ 针对矩阵形式的图像或视频数据。

  • 假设：数据内部有极强的相关性（如视频背景静止、图像包含重复纹理）。
  • 应用：去背景、去阴影、矩阵补全。

7. 优化算法谱系：你将学到什么？

一旦模型确定（即选定了 $D$ 和 $R$），剩下的就是纯粹的数学计算问题：求最小值。本书后续章节将详细介绍以下几类算法：

[ 优化方法族谱 ]

1. 光滑无约束 (Smooth Unconstrained)
   |-- 梯度下降 (Gradient Descent)  -> 也就是现在的深度学习训练基础
   |-- 牛顿法 / 拟牛顿法 (L-BFGS)   -> 收敛快，但内存消耗大
   |-- 共轭梯度 (Conjugate Gradient)-> 求解线性方程组 (Ax=b) 的神器

2. 近端方法 (Proximal Methods) -> 处理 R(x) 不可导的情况 (如 L1, TV)
   |-- ISTA (Iterative Shrinkage-Thresholding)
   |-- FISTA (Fast ISTA) -> 加速版，图像处理标配

3. 算子分裂 / 增广拉格朗日 (Splitting / ALM) -> 处理复杂约束
   |-- ADMM (交替方向乘子法) -> 极其灵活，万能胶水
   |-- Split Bregman -> 专门针对 L1/TV 问题的高效解法
   |-- Primal-Dual (Chambolle-Pock) -> 数学优美，处理鞍点问题

4. 非凸优化 (Non-convex) -> 更好的先验，更难的求解
   |-- IRLS (迭代重加权最小二乘)
   |-- MM (Majorization-Minimization)

选择指南：

  • 如果 $D$ 和 $R$ 都可导 $\to$ 梯度下降 / L-BFGS。
  • 如果 $R$ 是简单的 $\ell_1$ (Lasso) $\to$ FISTA。
  • 如果 $R$ 很复杂 (TV)，或者有硬约束 (像素值在0到1之间) $\to$ ADMM 或 Primal-Dual。

8. 本章小结

  1. 物理模型：$y=Ax+n$ 是一切的起点。$A$ 的性质决定了任务类型。
  2. 病态性：由于 $A$ 通常有零空间或衰减的奇异值，直接求逆会放大噪声。
  3. 变分法：通过 $\min D(Ax,y) + \lambda R(x)$ 求解。
    • $D$ 拟合数据（拒接噪声）。
    • $R$ 引入先验（补充信息）。
  4. 权衡：$\lambda$ 越大，结果越符合先验（规则、平滑）；$\lambda$ 越小，结果越贴近观测（可能包含噪声）。

9. 练习题

基础题

Q1.1 (算子的理解): 对于“去模糊”任务，如果模糊核是 $3 \times 3$ 的均匀模糊（每个元素都是 $1/9$），对于 $100 \times 100$ 的图像。

  1. 矩阵 $A$ 的维度是多少？
  2. $A$ 的每一行有多少个非零元素？
  3. 这说明 $A$ 具有什么性质？
Hint

每个输出像素只与周围的输入像素有关。

Answer
  1. $A$ 的维度是 $10000 \times 10000$ ($N=100\times100$)。
  2. 每一行对应输出图像的一个像素，它由输入图像的 $3 \times 3$ 区域加权和得到。所以除了边界行，大部分行有 9 个非零元素。
  3. $A$ 是稀疏矩阵 (Sparse Matrix)。这再次印证了我们不能存储 $A$，但可以快速计算 $Ax$（卷积）。

Q1.2 (L2 vs L1): 考虑一维信号恢复问题。观测值 $y=3$，模型 $x$。 情况 A: $D(x) = (x-3)^2$。 情况 B: $D(x) = |x-3|$。 如果现在增加一个离群观测 $y_{outlier} = 100$，两个观测一起考虑。 A 的目标: $\min_x (x-3)^2 + (x-100)^2$ B 的目标: $\min_x |x-3| + |x-100|$ 请分别求出最优解，并说明哪个对离群值更鲁棒。

Hint

A求导置0（即求平均值）。B寻找中位数。

Answer

A: $2(x-3) + 2(x-100) = 0 \implies x = 51.5$。结果被 100 严重拉偏。 B: 在 $[3, 100]$ 之间的任何 $x$ 都能使目标函数最小（值为 97）。通常取中位数。如果我们有三个点 $3, 3, 100$，中位数是 $3$。 L1 范数（B）更鲁棒。

Q1.3 (正则化参数): 在超分辨率任务中，如果 $\lambda \to \infty$，恢复出来的图像 $x$ 会变成什么样？（假设 $R(x) = |\nabla x|_2^2$）。

Hint

$\lambda \to \infty$ 意味着我们只在乎 $R(x)$ 最小，完全不在乎 $y$。

Answer

我们会得到一张常数图像（全灰或全黑），因为常数图像的梯度为 0，使 $R(x)$ 达到最小。这意味着我们完全忽略了观测数据。

挑战题

Q1.4 (贝叶斯与 MAP): 证明：如果先验 $P(x) \propto e^{-\lambda |x|_1}$（拉普拉斯分布），似然 $P(y|x) \propto e^{-\frac{1}{2\sigma^2}|Ax-y|_2^2}$（高斯分布）。推导出的 MAP 估计等价于 Lasso 问题（$\ell_2$ 数据项 + $\ell_1$ 正则）。并给出 $\lambda_{opt}$ 与噪声方差 $\sigma^2$ 以及拉普拉斯参数 $b$ 的关系。

Hint

写出 $-\log P(x|y)$，观察两项系数的比值。

Answer

$$ -\log P(y|x) \propto \frac{1}{2\sigma^2}|Ax-y|_2^2 $$ $$ -\log P(x) \propto \lambda_{laplace} |x|_1 $$ 总目标: $\frac{1}{2\sigma^2}|Ax-y|_2^2 + \lambda_{laplace} |x|_1$。 乘以 $2\sigma^2$ 后: $|Ax-y|_2^2 + (2\sigma^2 \lambda_{laplace}) |x|_1$。 所以正则化参数 $\lambda$ 与噪声方差 $\sigma^2$ 成正比。噪声越大，正则化强度应该越大。

Q1.5 (逆犯罪 Inverse Crime): 在学术论文的仿真实验中，如果你用模型 $y = A_{sim}x + n$ 生成数据，然后用完全相同的 $A_{sim}$ 去进行恢复 $\hat{x} = \arg\min |A_{sim}x - y|^2 + R(x)$，你的结果通常会好得出奇。 但在真实世界中，真实的退化过程 $A_{real}$ 永远无法被完美的数学模型 $A_{sim}$ 描述（例如离散化误差、未建模的光学畸变）。 这种现象被称为“逆犯罪”。如何设计实验避免这个问题？

Hint

生成数据时应该使用更精细的模型或不同的网格。

Answer

生成数据时，应该在更高分辨率或连续域模拟物理过程，然后降采样得到 $y$。 恢复时，只能在离散网格上构建 $A$。 例如：用 $2048 \times 2048$ 的图模拟模糊和噪声，然后下采样到 $512 \times 512$ 作为观测 $y$。恢复时假设 $x$ 是 $512 \times 512$ 的。这能模拟模型误差（Modeling Error）。

10. 常见陷阱与错误 (Gotchas)

  1. 混淆图像坐标与矩阵坐标：
    • 在数学公式中，图像是向量 $x$。在代码中，是矩阵 img[row, col]。
    • 陷阱：图像坐标通常是 $(x, y) = (col, row)$。在定义梯度算子或运动模糊方向时，搞反 $x$ 和 $y$ 是最常见的 Bug。
  2. $\lambda$ 的单位依赖性：
    • 许多初学者问：“$\lambda$ 应该设为 0.1 还是 1000？”
    • 这完全取决于你的图像数值范围！如果你的像素是 $[0, 1]$，数据项误差也是 $0.0x$ 量级。如果像素是 $[0, 255]$，误差可能是 $10000$ 量级。
    • Rule of Thumb：始终将图像归一化到 $[0, 1]$ (float) 进行计算。这样你的 $\lambda$ 经验值可以在不同图片间迁移。
  3. 忽略边界条件：
    • 当你计算 $\nabla x$ 或 $H * x$ 时，边缘像素怎么算？
    • FFT 卷积默认是周期性的（图像右边连着左边）。如果不做处理，边缘会出现强烈的十字伪影。
    • 对策：在去模糊前使用 edgetaper (Matlab) 或在变分模型中引入边界处理。
  4. 把优化当魔法：
    • 如果 $A$ 丢失的信息太多（如严重的模糊 + 强噪声），没有任何优化算法能找回细节。
    • 不要指望变分方法能像深度生成模型（GAN/Diffusion）那样“脑补”出原本不存在的纹理。变分法只能利用先验约束解空间，而不能创造信息。