zehua

菜单

Administrator
发布于 2026-09-09 / 0 阅读
0
0

Python TP4：图像反卷积优化

#实验报告#M2 ATSI

1. 病态逆问题

在数据处理中，感兴趣的参数通常是通过间接方式观测的。测量结果来自一个称为"正向模型"的系统，该系统将未知参数 x 与观测数据 y 联系起来。这个模型可以是仪器模型或物理模型，并依赖于额外的参数 \theta。在实际中，数据永远不会被模型完美地再现：这种偏差通常被归因于"噪声"。

1.1. 问题描述

这里的目标是消除图像的模糊，即从模糊图像中恢复清晰图像。这类问题属于反演问题，因为实际测量系统具有有限的分辨率、动态范围和响应时间。因此，测量值是感兴趣物理量的退化版本。

如果这些退化有时可以忽略不计，在某些情况下可能需要进行特定处理。例如，使用不完美对焦拍摄的图像是模糊的：一个点变成一个斑点，观测图像是这些斑点的叠加。

已知仪器的脉冲响应和模糊图像，问题就变成了重建清晰图像。

1.2. 数学建模

描述这种变换的最简单模型是线性时不变滤波器，即卷积。这样的系统在图1中以示意图形式表示。

x \longrightarrow \boxed{H} \longrightarrow \oplus \longrightarrow y
\uparrow
b

图1 — 正向模型。

*T. Rodet, F. Orieux

这里 x_{n,m} 表示真实图像，y_{n,m} 表示可用数据。为了考虑测量误差和建模误差，我们引入噪声 b_{n,m}。在缺乏关于该噪声结构的特定信息的情况下，我们假设它是（空间上的）白噪声。相应的观测方程形式如下

y_{n,m} = \sum_{p=0}^{P-1} \sum_{q=0}^{Q-1} h_{n-p,m-q} x_{p,q} + b_{n,m} \tag{1}

要完成的工作是根据测量值 y_{n,m} 构建 x_{n,m} 的估计量 \hat{x}_{n,m}。

需要特别注意的是，这里涉及的是低通滤波器。大多数测量设备都属于这种类型：它们无法"跟踪过高的频率"。因此很明显，系统会消除图像的高频分量，这就是为什么这个问题特别困难：必须在高频似乎不存在于数据中的情况下重建它们。

2. 一维反卷积

这样提出的图像反卷积问题是在估计理论的框架下，借助最小二乘方法和二次正则化来处理的。主要结果将在下面回顾，特别是所涉及的准则以及其最小值构成所寻求的图像估计量。

在实际中，为了在普通计算机上以合理的计算时间获得这些估计量，本题目建议在后续部分进行一些近似。这些近似的细节将在下一部分以一维形式呈现，即信号反卷积而非图像反卷积。这种呈现方式使我们能够以比二维更清晰的方式研究这些近似及其优势。向二维的过渡只是一维情况的扩展。尽管如此，在PYTHON下要实现的编程部分只涉及二维情况。

2.1. 一维建模

观测模型方程 (1) 在一维情况下变为

y_n = \sum_{p=0}^{P-1} h_p x_{n-p} + b_n.

如果我们有 N 个 n 的值，并且有 N 个观测值，可以将得到的方程以矩阵形式写成

\boldsymbol{y} = \boldsymbol{H}\boldsymbol{x} + \boldsymbol{b}.

向量 \boldsymbol{y} 汇集了所有观测数据（二维情况下为模糊图像），向量 \boldsymbol{x} 汇集了待恢复的信号（二维情况下为清晰图像），\boldsymbol{b} 汇集了噪声样本。矩阵 \boldsymbol{H}，称为卷积矩阵，具有以下经典的结构

\boldsymbol{H} = \begin{bmatrix} h_P & \cdots & h_1 & 0 & & \\ 0 & h_P & \cdots & h_1 & 0 & \\ & 0 & \ddots & \cdots & \ddots & 0 \\ & & 0 & h_P & \cdots & h_1 \\ & & & 0 & \ddots & \vdots \\ & & & & 0 & h_P \end{bmatrix}.

这是一个 N 维方阵（如果有 N 个观测值），具有Toeplitz结构。因此，现在的问题是根据数据向量 \boldsymbol{y} 和矩阵 \boldsymbol{H} 来估计向量 \boldsymbol{x}。

2.2. 最小二乘法

解决这个问题最简单的技术是最小化最小二乘准则

Q_{\text{MC}}(\boldsymbol{x}) = (\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x})^t(\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x}).

最小二乘意义下的估计对应于表达式的最小值

\hat{\boldsymbol{x}} = (\boldsymbol{H}^t\boldsymbol{H})^{-1}\boldsymbol{H}^t\boldsymbol{y}. \tag{2}

需要求逆的矩阵 \boldsymbol{H}^t\boldsymbol{H} 的大小为 N \times N，如果 N 很大（对于二维图像就是这种情况），直接求逆在计算成本方面是困难的。存在多种快速方法来求逆该矩阵并计算 \hat{\boldsymbol{x}}。

2.2.1. 循环近似求解

本实验中将讨论的第一种方法是基于循环矩阵的性质（参见附录）。这些性质允许用对角矩阵"替换"方程 (2) 中出现的矩阵。矩阵乘积和矩阵求逆在计算时间方面变得非常便宜。此外，这些矩阵向对角矩阵的"转换"是通过FFT实现的，然而为此必须接受卷积矩阵的循环近似。即用以下矩阵替换 \boldsymbol{H}

\tilde{\boldsymbol{H}} = \begin{bmatrix} h_P & \cdots & h_1 & 0 & & \\ 0 & h_P & \cdots & h_1 & 0 & \\ 0 & \ddots & \cdots & \ddots & 0 \\ & 0 & h_P & \cdots & h_1 \\ \ddots & & 0 & \ddots & \vdots \\ \cdots & h_1 & & & 0 & h_P \end{bmatrix}. \tag{3}

这种近似在于在矩阵 \boldsymbol{H} 的左下角添加一些项使其成为循环矩阵。需要注意的是，当 N（信号或图像的大小）远大于 P（脉冲响应的长度）时，这种近似效果更好。在实际中，对于 N = 256 个点的信号和 P = 10 的脉冲响应，这种近似的影响并不大。

由于矩阵 \tilde{\boldsymbol{H}} 是循环的，可以在傅里叶基下将其对角化（参见附录）

\tilde{\boldsymbol{H}} = \boldsymbol{F}^\dagger \boldsymbol{\Lambda}_h \boldsymbol{F}. \tag{4}

矩阵 \boldsymbol{\Lambda}_h 是对角矩阵，其对角元素是 \boldsymbol{H} 的特征值

\lambda_h^1, \lambda_h^2, \ldots, \lambda_h^N

这些特征值可以通过对 \boldsymbol{H} 的第一行进行FFT非常容易地获得，即对计算了 N 个点的脉冲响应进行FFT。换句话说

\boldsymbol{\lambda} = \begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_N \end{bmatrix} = \boldsymbol{F} \begin{bmatrix} h_1 \\ \vdots \\ h_P \\ 0 \\ \vdots \\ 0 \end{bmatrix}

从方程 (2) 出发，利用对角化关系式 (4) 以及课程中关于矩阵 \boldsymbol{F} 的关系，可以得到

\hat{\boldsymbol{x}} = \left(\tilde{\boldsymbol{H}}^t\tilde{\boldsymbol{H}}\right)^{-1}\tilde{\boldsymbol{H}}^t\boldsymbol{y} \tag{5}
= \left(\boldsymbol{F}^\dagger \boldsymbol{\Lambda}_h^\dagger \boldsymbol{F} \boldsymbol{F}^\dagger \boldsymbol{\Lambda}_h \boldsymbol{F}\right)^{-1} \boldsymbol{F}^\dagger \boldsymbol{\Lambda}_h^\dagger \boldsymbol{F} \boldsymbol{y}
= \boldsymbol{F}^\dagger \boldsymbol{\Lambda}_h^{-1} \boldsymbol{F} \boldsymbol{y}. \tag{6}

由此，左乘 \boldsymbol{F} 可得

\boldsymbol{F}\hat{\boldsymbol{x}} = \boldsymbol{\Lambda}_h^{-1}\boldsymbol{F}\boldsymbol{y}.

最后令 \hat{\boldsymbol{X}} = \boldsymbol{F}\hat{\boldsymbol{x}} 和 \boldsymbol{Y} = \boldsymbol{F}\boldsymbol{y} 为傅里叶变换后的 \hat{\boldsymbol{x}} 和 \boldsymbol{y}。注意 \hat{\boldsymbol{X}} 和 \boldsymbol{Y} 也可以通过对 \hat{\boldsymbol{x}} 和 \boldsymbol{y} 进行FFT获得。在傅里叶域中，问题变为

\hat{\boldsymbol{X}} = \boldsymbol{\Lambda}_h^{-1}\boldsymbol{Y}

最后，构造向量 \boldsymbol{g}_{\text{MC}}，其元素为 \boldsymbol{\Lambda}_h^{-1} 的对角线值

g_{\text{MC}}^n = 1/\lambda_n \quad \text{pour} \quad n = 1, 2, \ldots N. \tag{7}

然后通过两个向量 \boldsymbol{g}_{\text{MC}} 和 \boldsymbol{Y} 的"逐元素乘积"（记作 \cdot，称为Hadamard乘积）得到向量 \hat{\boldsymbol{X}}

\hat{\boldsymbol{X}} = \boldsymbol{g}_{\text{MC}} \cdot \boldsymbol{Y}.

注意 \boldsymbol{Y}（\boldsymbol{y} 的傅里叶变换）和 \hat{\boldsymbol{X}}（\boldsymbol{x} 的傅里叶变换）之间的这个关系是一个"傅里叶域滤波"关系，其"离散频率传递函数"是 \boldsymbol{g}_{\text{MC}}。总之，要遵循的步骤是

  1. 构造脉冲响应在 N 个点上的FFT \boldsymbol{\lambda}_h。
  2. 根据方程 (7) 构造频率传递函数向量 \boldsymbol{g}_{\text{MC}}。
  3. 构造数据的FFT \boldsymbol{Y}。
  4. 通过 \boldsymbol{Y} 乘以 \boldsymbol{g}_{\text{MC}} 计算 \hat{\boldsymbol{X}}。
  5. 通过对 \hat{\boldsymbol{X}} 进行逆FFT回到空间域，得到 \hat{\boldsymbol{x}}。

2.2.2. 梯度下降算法

本实验中将讨论的第二种方法是梯度下降法。由于准则是二次的，它是凸的并且只有一个全局最小值，参见方程 (2)。为了找到准则的最小值，寻找最陡下降的方向（这就是所谓的梯度方向），然后沿着这个方向移动。每次重复都会接近最小值。

首先回顾准则 Q_{\text{MC}} 的表达式

Q_{\text{MC}}(\boldsymbol{x}) = (\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x})^t(\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x})
= \boldsymbol{x}^t\boldsymbol{H}^t\boldsymbol{H}\boldsymbol{x} - 2\boldsymbol{x}^t\boldsymbol{H}^t\boldsymbol{y} + \boldsymbol{y}^t\boldsymbol{y}.

为了计算准则的梯度，对向量 \boldsymbol{x} 求导

\frac{\partial Q_{\text{MC}}(\boldsymbol{x})}{\partial \boldsymbol{x}} = 2\boldsymbol{H}^t\boldsymbol{H}\boldsymbol{x} - 2\boldsymbol{H}\boldsymbol{y}
= 2\boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x} - \boldsymbol{y}).

记 \boldsymbol{e}_k = (\boldsymbol{H}\boldsymbol{x} - \boldsymbol{y}) 为模型输出与数据之间的误差。为了实现算法，取任意的 \boldsymbol{x}^0，然后应用以下迭代

\boldsymbol{x}^{k+1} = \boldsymbol{x}^k - \alpha\boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x}^k - \boldsymbol{y}). \tag{8}

如果

\alpha < \frac{2}{\|\boldsymbol{H}^t\boldsymbol{H}\|_M},

则算法收敛到准则的最小值，其中 \|\cdot\|_M 是矩阵范数，等于最大奇异值的模的平方。总之，要遵循的步骤是

  1. 将 \boldsymbol{x}^0 初始化为任意值。

  2. 通过将 \boldsymbol{x}^k 与脉冲响应卷积来计算 \boldsymbol{H}\boldsymbol{x}^k。

  3. 计算误差 \boldsymbol{e}_k，(\boldsymbol{H}\boldsymbol{x}^k - \boldsymbol{y})。

  4. 通过将上一步的结果与翻转的脉冲响应卷积来计算 \boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x}^k - \boldsymbol{y})。

  5. 将上一步的结果乘以 \alpha 并从 \boldsymbol{x}^k 中减去。

  6. 返回步骤2，直到算法收敛。

你会发现结果在收敛时并不令人满意。唯一的解决方案是进行早停（early stopping）。

2.3. 正则化或惩罚

当你使用最小二乘方法时，你会发现它并不总是令人满意的。在某些情况下，得到的图像仍然非常嘈杂，非常不稳定。你需要解释为什么会这样，这实际上是第4.2.3节中向你提出的问题之一。

无论如何，为了纠正这个缺陷，可以修改最小二乘准则以确保得到的图像具有一定的规则性。最简单的想法是惩罚待恢复信号相邻样本之间的强烈差异。准则于是采用以下正则化形式

Q_{\text{MCR}}(\boldsymbol{x}) = (\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x})^t(\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x}) + \mu\sum_{n=1}^{N-1}(x_{n+1} - x_n)^2
= (\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x})^t(\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x}) + \mu\boldsymbol{x}^t\boldsymbol{D}^t\boldsymbol{D}\boldsymbol{x}

其中 \boldsymbol{D} 是一阶差分矩阵，大小为 (N-1) \times N，定义为

\boldsymbol{D} = \begin{bmatrix} 1 & -1 & 0 & 0 & \ldots & \ldots & 0 \\ 0 & 1 & -1 & 0 & \ldots & \ldots & 0 \\ \vdots & & \ddots & \ddots & & & \vdots \\ \vdots & & & \ddots & \ddots & & \vdots \\ 0 & \ldots & \ldots & 0 & 1 & -1 & 0 \\ 0 & \ldots & \ldots & 0 & 0 & 1 & -1 \end{bmatrix}.

这个准则和简单最小二乘准则一样，是关于 \boldsymbol{x} 的二次函数，其最小值的表达式已知

\hat{\boldsymbol{x}} = (\boldsymbol{H}^t\boldsymbol{H} + \mu\boldsymbol{D}^t\boldsymbol{D})^{-1}\boldsymbol{H}^t\boldsymbol{y}.

2.3.1. 循环近似

出于与简单最小二乘情况相同的原因，我们接受对 \boldsymbol{H} 和 \boldsymbol{D} 的循环近似，以便用对角矩阵上的计算来"替换"任意矩阵上的计算。关于 \boldsymbol{H}，我们保留方程 (3) 给出的近似 \tilde{\boldsymbol{H}}。关于 \boldsymbol{D}，只需添加最后一行即可得到

循环近似 \tilde{\boldsymbol{D}}，其维度为 N \times N

\tilde{\boldsymbol{D}} = \begin{bmatrix} 1 & -1 & 0 & 0 & \ldots & \ldots & 0 \\ 0 & 1 & -1 & 0 & \ldots & \ldots & 0 \\ \vdots & & \ddots & \ddots & & & \vdots \\ \vdots & & & \ddots & \ddots & & \vdots \\ 0 & \ldots & \ldots & 0 & 1 & -1 & 0 \\ 0 & \ldots & \ldots & 0 & 0 & 1 & -1 \\ -1 & \ldots & \ldots & 0 & 0 & 0 & 1 \end{bmatrix}.

近似的最小值于是采用以下形式

\hat{\boldsymbol{x}} = (\tilde{\boldsymbol{H}}^t\tilde{\boldsymbol{H}} + \mu\tilde{\boldsymbol{D}}^t\tilde{\boldsymbol{D}})^{-1}\tilde{\boldsymbol{H}}^t\boldsymbol{y}

从这个表达式出发，利用附录中关于矩阵 \boldsymbol{F} 的关系，并参考简单最小二乘情况下方程 (5) 到 (6) 的计算，证明

\hat{\boldsymbol{x}} = \boldsymbol{F}^\dagger(\boldsymbol{\Lambda}_h^\dagger\boldsymbol{\Lambda}_h + \boldsymbol{\Lambda}_d^\dagger\boldsymbol{\Lambda}_d)^{-1}\boldsymbol{\Lambda}_h^\dagger\boldsymbol{F}\boldsymbol{y}

其中 \boldsymbol{\Lambda}_h 和 \boldsymbol{\Lambda}_d 分别是 \tilde{\boldsymbol{H}} 和 \tilde{\boldsymbol{D}} 的特征值对角矩阵。我们继续记 \boldsymbol{\lambda}_h = \left[\lambda_h^1, \lambda_h^2, \ldots, \lambda_h^N\right] 为 \tilde{\boldsymbol{H}} 的特征值向量，记 \boldsymbol{\lambda}_d = \left[\lambda_d^1, \lambda_d^2, \ldots, \lambda_d^N\right] 为 \tilde{\boldsymbol{D}} 的特征值向量。在你的报告中说明如何获得这些特征值。

如前所述，令 \hat{\boldsymbol{X}} = \boldsymbol{F}\hat{\boldsymbol{x}} 和 \boldsymbol{Y} = \boldsymbol{F}\boldsymbol{y} 为 \hat{\boldsymbol{x}} 和 \boldsymbol{y} 的傅里叶变换，并推导出

\hat{\boldsymbol{X}} = (\boldsymbol{\Lambda}_h^\dagger\boldsymbol{\Lambda}_h + \mu\boldsymbol{\Lambda}_d^\dagger\boldsymbol{\Lambda}_d)^{-1}\boldsymbol{\Lambda}_h^\dagger\boldsymbol{Y}.

最后，构造向量 \boldsymbol{g}_{\text{MCR}}，其元素为矩阵 (\boldsymbol{\Lambda}_h^\dagger\boldsymbol{\Lambda}_h + \mu\boldsymbol{\Lambda}_d^\dagger\boldsymbol{\Lambda}_d)^{-1}\boldsymbol{\Lambda}_h^\dagger 的对角线值

g_{\text{MCR}}^n = \frac{(\lambda_h^n)^*}{|\lambda_h^n|^2 + \mu|\lambda_d^n|^2} \quad \text{pour} \quad n = 1, 2, \ldots, N. \tag{9}

然后解释为什么向量 \hat{\boldsymbol{X}} 可以通过两个向量 \boldsymbol{g}_{\text{MCR}} 和 \boldsymbol{Y} 的"逐元素乘积"获得

\hat{\boldsymbol{X}} = \boldsymbol{g}_{\text{MCR}} \cdot \boldsymbol{Y}

注意这里仍然是一个"傅里叶域滤波"关系，其"离散传递函数"是 \boldsymbol{g}_{\text{MCR}}。总之，要遵循的步骤是

  1. 构造脉冲响应在 N 个点上的FFT \boldsymbol{\lambda}_h。
  2. 构造 [1, -1] 在 N 个点上的FFT \boldsymbol{\lambda}_d。
  3. 根据方程 (9) 构造频率传递函数向量 \boldsymbol{g}_{\text{MCR}}。
  4. 构造数据的FFT \boldsymbol{Y}。
  5. 计算 \hat{\boldsymbol{X}} 作为频率传递函数 \boldsymbol{g}_{\text{MCR}} 乘以 \boldsymbol{Y} 的乘积。
  6. 通过对 \hat{\boldsymbol{X}} 进行逆FFT回到空间域，得到 \hat{\boldsymbol{x}}。

2.3.2. 梯度下降算法

用新准则 Q_{\text{MCR}} 重写第2.2.2部分的方程

Q_{\text{MCR}}(\boldsymbol{x}) = (\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x})^t(\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x}) + \mu\boldsymbol{x}^t\boldsymbol{D}^t\boldsymbol{D}\boldsymbol{x}
= \boldsymbol{x}^t\boldsymbol{H}^t\boldsymbol{H}\boldsymbol{x} - 2\boldsymbol{x}^t\boldsymbol{H}^t\boldsymbol{y} + \boldsymbol{y}^t\boldsymbol{y} + \mu\boldsymbol{x}^t\boldsymbol{D}^t\boldsymbol{D}\boldsymbol{x}.

为了计算准则的梯度，对向量 \boldsymbol{x} 求导

\frac{\partial Q_{\text{MCR}}(\boldsymbol{x})}{\partial \boldsymbol{x}} = 2\boldsymbol{H}^t\boldsymbol{H}\boldsymbol{x} - 2\boldsymbol{H}\boldsymbol{y} + 2\mu\boldsymbol{D}^t\boldsymbol{D}\boldsymbol{x} \tag{10}
= 2\boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x} - \boldsymbol{y}) + 2\mu\boldsymbol{D}^t\boldsymbol{D}\boldsymbol{x}. \tag{11}

为了实现算法，取任意的 \boldsymbol{x}^0，然后应用以下迭代

\boldsymbol{x}^{k+1} = \boldsymbol{x}^k - \alpha\left(\boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x}^k - \boldsymbol{y}) + \mu\boldsymbol{D}^t\boldsymbol{D}\boldsymbol{x}^k\right) \tag{12}

其中最优步长为

\alpha_{opt} = \frac{\boldsymbol{g}_k^t\boldsymbol{g}_k}{\boldsymbol{g}_k^t(\boldsymbol{H}^t\boldsymbol{H} + \mu\boldsymbol{D}^t\boldsymbol{D})\boldsymbol{g}_k} \tag{13}

其中 \boldsymbol{g}_k 是在 \boldsymbol{x}_k 处的梯度

\boldsymbol{g}_k = \left(\boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x}^k - \boldsymbol{y}) + \mu\boldsymbol{D}^t\boldsymbol{D}\boldsymbol{x}^k\right). \tag{14}

总之，要遵循的步骤是

  1. 将 \boldsymbol{x}^0 初始化为任意值。
  2. 通过将 \boldsymbol{x}^k 与脉冲响应卷积来计算 \boldsymbol{H}\boldsymbol{x}^k。
  3. 计算 \boldsymbol{H}\boldsymbol{x}^k 与数据 \boldsymbol{y} 之间的差 \boldsymbol{e}_k。
  4. 通过将上一步的结果与翻转的脉冲响应卷积来计算 \boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x}^k - \boldsymbol{y})。
  5. 通过将 \boldsymbol{x}^k 与微分滤波器 [-1, 1] 卷积来计算 \boldsymbol{D}\boldsymbol{x}^k。
  6. 通过将 \boldsymbol{D}\boldsymbol{x}^k 与翻转的微分滤波器 [1, -1] 卷积来计算 \boldsymbol{D}^t\boldsymbol{D}\boldsymbol{x}^k。
  7. 使用方程 (12) 更新 \boldsymbol{x}^k。
  8. 返回步骤2，直到算法收敛。

2.4. 扩展到二维

2.4.1. 循环近似方法

可以进行与一维情况类似的严格计算来建立二维情况的结果。然而，这些计算要繁琐得多，因为矩阵具有更复杂的分块Toeplitz结构，每个块本身也具有Toeplitz结构。此外，在两个方向上的循环近似更难以处理。这就是为什么我们不将二维情况作为一维情况的扩展来讨论。但有几点说明。

— 图像、脉冲响应、平滑项现在都是二维的，这就是为什么FFT被FFT2D取代。

— 更准确地说，如果待处理的图像有 N 行 N 列，则所有FFT2D都应在 N 行 N 列上进行。
— 频率传递函数也是二维的：每个空间频率都有一个维度，对应两个空间频率。
— 算子 \boldsymbol{D} 在二维中分解为两个算子 \boldsymbol{D}_1 和 \boldsymbol{D}_2，它们分别在图像的行和列上进行一维导数运算。

从一维情况推导出的步骤结构如下所示，首先是简单最小二乘情况，然后是正则化最小二乘情况。

  1. 构造脉冲响应在 N \times N 个点上的FFT2D \boldsymbol{\lambda}_h。

  2. 根据方程 (7) 构造频率传递函数矩阵 \boldsymbol{g}_{\text{MC}}。

  3. 构造数据的FFT2D \boldsymbol{Y}。

  4. 计算 \hat{\boldsymbol{X}} 为 \boldsymbol{g}_{\text{MC}} 乘以 \boldsymbol{Y} 的逐元素乘积。

  5. 通过对 \hat{\boldsymbol{X}} 进行逆FFT2D回到空间域，得到 \hat{\boldsymbol{x}}。

  6. 构造脉冲响应在 N \times N 个点上的FFT2D \boldsymbol{\lambda}_h。

  7. 构造 (0 \; 1 \; \text{-}1) 在 N \times N 个点上的FFT2D \boldsymbol{\lambda}_{d1}。

  8. 构造

\begin{pmatrix} 0 \\ 1 \\ \text{-}1 \end{pmatrix}

在 N \times N 个点上的FFT2D \boldsymbol{\lambda}_{d2}。

  1. 根据方程 (9) 构造频率传递函数矩阵 \boldsymbol{g}_{\text{MCR}}。
  2. 构造数据的FFT2D \boldsymbol{Y}。
  3. 计算 \hat{\boldsymbol{X}} 为 \boldsymbol{g}_{\text{MC}} 乘以 \boldsymbol{Y} 的逐元素乘积。
  4. 通过对 \hat{\boldsymbol{X}} 进行逆FFT2D回到空间域，得到 \hat{\boldsymbol{x}}。

二维梯度下降方法的扩展是通过用二维卷积运算替换卷积运算来实现的，借助函数 scipy.signal.convolve2d

注1. 函数 scipy.ndimage.convolve 对图像边界的处理接受不同的参数。在我们的情况下
— 计算正向，即 \boldsymbol{H}，使用选项 'valid'。
— 计算伴随，即 \boldsymbol{H}^t，使用选项 'full'。

注2. 如果要用函数 scipy.ndimage.convolve 计算 \boldsymbol{H}，你需要使用二维核 d

d = np.random.standard_normal((5, 5))

则 \boldsymbol{H}^t 的计算是用翻转的核进行的，即

d_transpose = np.fliplr(np.flipud(d))

3. 图像反卷积

3.1. 待处理图像

在PYTHON中加载文件 data1.npz（用于循环近似算法）和 data2.npz（用于梯度下降算法）中提供的数据。这些文件包含模糊图像、完美图像（scipy的ascent）、脉冲响应 ir 和滤波器的频率响应 fr。

绘制脉冲响应 h_{n,m} 及其对应的频率传递函数 H(\nu_x, \nu_y)。花时间识别两个频率轴、零频率、最大频率等。在报告中清楚地呈现这些结果。滤波器的性质是什么？

3.2. 反卷积的实现

要实现的PYTHON程序必须通过最小化两个准则来实现反卷积：简单最小二乘和正则化最小二乘，均为二维。

3.2.1. 循环近似反演

你来到了本次实验最重要的部分：结果的分析和评论。

注意，当然要看获得图像的整体外观，特别观察边界的清晰度，例如摄影师的外套和天空之间的边界。也要关注面部或相机设备的细节。

使用逆滤波器会发生什么？为什么？使用正则化方法，定性地解释获得的结果。参数 \mu 的影响是什么？将获得的结果与考虑到正则化项形式的预期结果进行比较。

3.2.2. 使用梯度下降算法的反演

在准则 Q_{MC} 上用不同的 \alpha 值测试你的算法： \alpha = 1.5，\alpha = 0.015，\alpha = 2.03。你注意到什么？在准则 Q_{MCR} 上用 \mu = 0.005 和 0 < \alpha < 0.2 测试你的算法。将获得的结果与之前的结果进行比较，特别是在图像边界方面。也比较两种方法的计算时间。编写计算每一步最优alpha的函数。

4. 凸优化

在前面的部分中，我们观察到考虑关于所寻找对象的某种平滑度的先验信息可以获得令人满意的解。然而，这个先验假设图像中没有轮廓（灰度级的剧烈变化）。

因此，我们将引入一个允许强烈惩罚小变化的先验信息

像素之间的差异，这对应于噪声，同时对相邻像素之间的重要差异惩罚较轻，这对应于图像的轮廓。

4.1. 理论研究

4.1.1. 准则的构造

我们希望像之前的实验一样，有一个数据保真项和一个项来惩罚连续像素之间的差异，借助一个稍后定义的函数 \varphi：

Q_{\text{MCR}_2}(\boldsymbol{x}) = (\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x})^t(\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x}) + \mu\sum_{n=1}^{N-1}\varphi(x_{n+1} - x_n) \tag{15}
= (\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x})^t(\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x}) + \mu\boldsymbol{\Phi}(\boldsymbol{D}\boldsymbol{x}) \tag{16}

4.1.2. 关于函数 \varphi 作用的讨论

函数 \varphi 用于惩罚相邻像素之间的差异。如果取 \varphi(\tau) = \tau^2，则得到与前一次实验相同的正则化最小二乘准则。我们在重建中观察到图像的轮廓被平滑了，因为像素之间的重要差异被强烈惩罚。为了更好地保留轮廓，同时惩罚由噪声引起的图像波动，我们将使用Huber函数：

\varphi(\tau) = \begin{cases} \tau^2 & \text{si } |\tau| \leq s, \\ 2s|\tau| - s^2 & \text{si } |\tau| \geq s. \end{cases} \tag{17}

图2 — 凸函数 \varphi 的例子：虚线为 x^2，实线为阈值 s = 0.5 时的Huber函数。

从图2中可以观察到Huber函数始终是凸函数，并且该函数在阈值 s 以下具有二次行为，在阈值 s 以上具有线性行为。这个函数对像素之间对应于轮廓的高差异惩罚较轻。

4.1.3. 准则的最小化

方程 (15) 的准则是凸准则，因为数据保真项是二次的，而 \boldsymbol{x} 方差的惩罚项是凸的。但它不再是二次的，这意味着我们不再能显式地知道准则的最小值。因此不可能直接使用基于循环矩阵求逆的快速算法。

但由于准则是严格凸的，我们可以通过下降算法找到最小值，例如梯度方向下降算法。

梯度计算 首先展开方程 (15) 的准则

Q_{\text{MCR}_2}(\boldsymbol{x}) = (\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x})^t(\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x}) + \mu\boldsymbol{\Phi}(\boldsymbol{D}\boldsymbol{x})
= \boldsymbol{x}^t\boldsymbol{H}^t\boldsymbol{H}\boldsymbol{x} - 2\boldsymbol{x}^t\boldsymbol{H}^t\boldsymbol{y} + \boldsymbol{y}^t\boldsymbol{y} + \mu\boldsymbol{\Phi}(\boldsymbol{D}\boldsymbol{x})

现在对准则关于 \boldsymbol{x} 求导来计算梯度

\frac{\partial Q_{\text{MCR}_2}(\boldsymbol{x})}{\partial \boldsymbol{x}} = \boldsymbol{g} = 2\boldsymbol{H}^t\boldsymbol{H}\boldsymbol{x} - 2\boldsymbol{H}^t\boldsymbol{y} + \mu\boldsymbol{D}^t\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}\boldsymbol{x}) \tag{18}
= 2\boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x} - \boldsymbol{y}) + \mu\boldsymbol{D}^t\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}\boldsymbol{x}) \tag{19}

其中 \boldsymbol{\Phi} 的梯度写为

\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{x}) = \begin{pmatrix} \varphi'(x_1) \\ \vdots \\ \varphi'(x_n) \end{pmatrix} \tag{20}

其中

\varphi'(\tau) = \begin{cases} -2s & \text{si } \tau \leq -s, \\ 2\tau & \text{si } |\tau| \leq s, \\ 2s & \text{si } \tau \geq s. \end{cases} \tag{21}

梯度方向下降算法 用 \boldsymbol{x}^0 的任意值初始化，然后沿梯度的反方向下降来最小化准则

\boldsymbol{x}^{(k+1)} = \boldsymbol{x}^{(k)} - \alpha\left(\boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x}^{(k)} - \boldsymbol{y}) + \mu\boldsymbol{D}^t\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}\boldsymbol{x}^{(k)})\right). \tag{22}

二次近似下的最优步长

\alpha_{opt} = \frac{\boldsymbol{g}_k^t\boldsymbol{g}_k}{\boldsymbol{g}_k^t(\boldsymbol{H}^t\boldsymbol{H} + \mu\boldsymbol{D}^t\frac{\partial^2\boldsymbol{\Phi}}{\partial\boldsymbol{x}\partial\boldsymbol{x}^t}(\boldsymbol{D}\boldsymbol{x}^{(k)})\boldsymbol{D})\boldsymbol{g}_k} \tag{23}

其中Hessian矩阵为

\frac{\partial^2\boldsymbol{\Phi}}{\partial\boldsymbol{x}\partial\boldsymbol{x}^t}(\boldsymbol{x}) = \begin{pmatrix} \varphi''(x_1) & 0 & \ldots & 0 \\ 0 & \varphi''(x_2) & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \ldots & 0 & \varphi''(x_n) \end{pmatrix} \tag{24}

其中

\varphi''(\tau) = \begin{cases} 2 & \text{si } |\tau| \leq s, \\ 0 & \text{sinon.} \end{cases} \tag{25}
  1. 将 \boldsymbol{x}^0 初始化为任意值。
  2. 通过将 \boldsymbol{x}^{(k)} 与脉冲响应卷积来计算 \boldsymbol{H}\boldsymbol{x}^{(k)}。
  3. 计算 \boldsymbol{H}\boldsymbol{x}^{(k)} 与数据 \boldsymbol{y} 之间的差 e_k。
  4. 通过将上一步的结果与翻转的RI卷积来计算 \boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x}^{(k)} - \boldsymbol{y})（使用 rotate 命令）。
  5. 通过将 \boldsymbol{x}^{(k)} 与微分滤波器 [-1 \; 1] 卷积来计算 \boldsymbol{D}\boldsymbol{x}^{(k)}。
  6. 借助方程 (20) 计算 \frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}\boldsymbol{x}^{(k)})。
  7. 通过将 \frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}\boldsymbol{x}^{(k)}) 与翻转的微分滤波器 [1 \; -1] 卷积来计算 \boldsymbol{D}^t\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}\boldsymbol{x}^{(k)})。
  8. 借助方程 (23) 计算 \alpha。
  9. 使用方程 (22) 更新 \boldsymbol{x}^{(k)}。
  10. 返回步骤2，直到算法收敛。（返回步骤2，直到梯度的范数大于某个我们任意设定的阈值）。

4.1.4. 扩展到二维

可以进行与一维情况类似的严格计算来建立二维情况的结果。然而，这些计算要繁琐得多。这就是为什么我们不将二维情况作为一维情况的扩展来讨论。但有几点说明。

— 图像、脉冲响应、平滑项现在都是二维的，这就是为什么卷积运算将是二维的。
— 算子 \boldsymbol{D} 在二维中分解为两个算子 \boldsymbol{D}_1 和 \boldsymbol{D}_2，它们分别在图像的行和列上进行一维导数运算。
— 微分算子 \boldsymbol{D}_1 和 \boldsymbol{D}_2 分别使用 convolution 命令计算。

方程 (22) 在二维中变为

\boldsymbol{x}^{(k+1)} = \boldsymbol{x}^{(k)} - \alpha\Big(\boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x}^{(k)} - \boldsymbol{y}) + \mu\Big(\boldsymbol{D}_1^t\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}_1\boldsymbol{x}^{(k)}) + \boldsymbol{D}_2^t\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}_2\boldsymbol{x}^{(k)})\Big)\Big). \tag{26}

从一维情况推导出的步骤结构如下所示。

  1. 将 \boldsymbol{x}^0 初始化为任意值，例如零。
  2. 通过使用 convolution 命令将 \boldsymbol{x}^{(k)} 与脉冲响应卷积来计算 \boldsymbol{H}\boldsymbol{x}^{(k)}。
  3. 计算 \boldsymbol{H}\boldsymbol{x}^{(k)} 与数据 \boldsymbol{y} 之间的差 e_k。
  4. 通过将上一步的结果与翻转的RI卷积来计算 \boldsymbol{H}^t(\boldsymbol{H}\boldsymbol{x}^{(k)} - \boldsymbol{y})（使用 rotate 命令）。
  5. 通过使用 convolution 将 \boldsymbol{x}^{(k)} 与微分滤波器 (0 \; 1 \; \text{-}1) 卷积来计算 \boldsymbol{D}_1\boldsymbol{x}^{(k)}。
  6. 通过微分滤波器
\begin{pmatrix} 0 \\ 1 \\ -1 \end{pmatrix}

借助 convolution 将 \boldsymbol{x}^{(k)} 卷积来计算 \boldsymbol{D}_2\boldsymbol{x}^{(k)}。

  1. 借助方程 (20) 计算 \frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}_1\boldsymbol{x}^{(k)}) 和 \frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}_2\boldsymbol{x}^{(k)})。
  2. 通过将 \frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}_1\boldsymbol{x}^{(k)}) 与翻转的微分滤波器卷积来计算 \boldsymbol{D}_1^t\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}_1\boldsymbol{x}^{(k)})
\begin{pmatrix} -1 & 1 & 0 \end{pmatrix}.
  1. 通过将 \frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}_2\boldsymbol{x}^{(k)}) 与翻转的微分滤波器卷积来计算 \boldsymbol{D}_2^t\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{x}}(\boldsymbol{D}_2\boldsymbol{x}^{(k)})
\begin{pmatrix} -1 \\ 1 \\ 0 \end{pmatrix}.
  1. 借助函数 AlphaOptConv 计算 \alpha。
  2. 使用方程 (26) 更新 \boldsymbol{x}^{(k)}。
  3. 返回步骤2，直到算法收敛。（返回步骤2，直到梯度的范数大于某个我们任意设定的阈值）。

4.2. 提供图像的反卷积

4.2.1. 待处理图像

你可以在前一部分已经处理过的数据上测试算法。

4.2.2. 二维反演的实现

要实现的程序必须实现反卷积方法，通过带有凸惩罚的正则化最小二乘对二维图像进行处理。

4.2.3. 评论

你来到了本次实验最重要的部分：结果的分析和评论。

注意，当然要看获得图像的整体外观，特别观察边界的清晰度。也要关注细节。

改变Huber函数的阈值 s 的值，并解释结果。

将获得的结果与前一次实验中的预期结果进行比较。分析不同类型正则化的作用。

A. 循环矩阵的对角化

本附录汇集了关于在傅里叶基下循环矩阵对角化的线性代数结果。这些结果使得所提方法的快速实现成为可能。

A.1. 循环矩阵

向量的循环移位（或旋转）

为了对向量进行循环移位（或旋转），只需将其元素向右移动一个位置，并将最后一个元素放到第一个位置，如下所示

[c_1, c_2, c_3, \ldots, c_{N-1}, c_N] \quad \longrightarrow \quad [c_N, c_1, c_2, c_3, \ldots, c_{N-1}].

循环矩阵

如果矩阵的每一行都是通过前一行的循环移位得到的，且第一行是通过最后一行的循环移位得到的，则称该矩阵是循环矩阵。因此，循环矩阵的形式为

\boldsymbol{C} = \begin{bmatrix} c_1 & c_2 & c_3 & \ldots & & c_N \\ c_N & c_1 & c_2 & & & c_{N-1} \\ c_{N-1} & c_N & c_1 & \ddots & & c_{N-2} \\ \vdots & & & \ddots & & \vdots \\ c_3 & c_4 & c_5 & & & c_2 \\ c_2 & c_3 & c_4 & \ldots & c_N & c_1 \end{bmatrix}.

首先，我们注意到这个矩阵具有Toeplitz结构，即其元素在对角线上是常数：c_1 在主对角线上，c_2 在第一条上对角线上，c_N 在第一条下对角线上，等等，反之亦然。我们还注意到每一列都是通过前一列的循环移位得到的，第一列是通过最后一列的循环移位得到的。因此，可以通过对行或列进行移位来构造循环矩阵。

本实验中提出的方法的快速实现利用了以下结果：循环矩阵在傅里叶基下是可对角化的。这个结果允许用对角矩阵的乘积来替换任意矩阵的乘积（更加经济），代价是进行一些FFT计算。

A.2. 傅里叶矩阵及其一些性质

记 \boldsymbol{F} 为傅里叶矩阵，它是大小为 N 的方阵，即大小为 N 的向量在 N 个点上的"FFT矩阵"。它由以下关系式表征：\boldsymbol{X} = \boldsymbol{F}\boldsymbol{x}，其中

— \boldsymbol{x} 是大小为 N 的任意向量，且
— \boldsymbol{X} 是大小为 N 的向量，包含其傅里叶变换在 0 到 1 - 1/N 之间均匀分布的 N 个频率点上的 N 个样本。

这个矩阵具有一些特殊的性质，对于第2.1节的计算很有用

\boldsymbol{F}^t = \bar{\boldsymbol{F}}
\boldsymbol{F}\boldsymbol{F}^\dagger = \boldsymbol{I}_N
\boldsymbol{F}^{-1} = \boldsymbol{F}^\dagger = \boldsymbol{F}^*.

A.3. 对角化

可以证明 \boldsymbol{C} 在傅里叶基 \boldsymbol{F} 下是可对角化的，即

\boldsymbol{C} = \boldsymbol{F}^\dagger\boldsymbol{\Lambda}_C\boldsymbol{F}

其中 \boldsymbol{\Lambda}_C 是 \boldsymbol{C} 的特征值对角矩阵。还可以证明，这些特征值可以通过矩阵 \boldsymbol{C} 的第一行 \boldsymbol{c} = [c_1, c_2, \ldots, c_N]^t 的傅里叶变换获得

\boldsymbol{\lambda}_C = \boldsymbol{F}\boldsymbol{c}

其中 \boldsymbol{\lambda}_C 是包含 \boldsymbol{C} 的 N 个特征值的向量，即 \boldsymbol{\Lambda}_C 的对角元素。向量 \boldsymbol{\lambda}_C 显然可以通过向量 \boldsymbol{c} 的FFT获得。

Python TP3：NumPy ...
数学工具（一）：概率论与矩阵代数（...

评论