1 引言

受各种记录设备所固有的物理限制, 大多数真实世界中获取到的图像都包含了一定程度的噪声和失真^[1], 下面的模型可以描述期望图像退化至观测图像的过程:

$ y = A{\textit{z}} + e $

(1)

其中, $ y \in {R^N} $ 是由观测图像的像素组成的 $ N $ 维列向量, $ {\textit{z}} \in {R^N} $ 是对应的期望图像, $ e \in {R^N} $ 假定是均值为0, 方差为 $ {\sigma ^2} $ 的加性高斯白噪声, 噪声主要发生在图像采集和传输过程中的模拟电路中, 是实际图像中最常见的噪声. 另外, $ A \in {R^{N \times N}} $ 是由模糊核产生的降级矩阵, 在我们所关注的去噪问题中, $ A = I $ .

一类常见的去噪算法基于某些特定变换域, 例如傅里叶变换、快速傅里叶变换、离散余弦变换、小波变换, 这类算法主要利用稀疏性, 使用较少的非零系数来表示信号. 确定变换域后还需要确定阈值, 阈值函数一般有硬阈值和软阈值两种选择^[2]. 硬阈值方法将小波系数全部切除, 其余的全部保留, 软阈值方法为了避免频域中的阶跃, 将保留的系数都减去了阈值.

另一类方法基于空域, 直接在原图像上进行处理, 这类方法还可分为局部和非局部两种. 常见的局部方法包括高斯滤波器、最小均方滤波器、双边滤波器、维纳滤波器, 局部算法的基本思想是利用像素之间的相关性, 而像素的相关性会由于噪声水平变高而严重受损, 因此局部方法在高噪声条件下表现不佳. 非局部方法中最常用的是非局部均值滤波(NLM), 它利用了图像的非局部相似性, 是一种典型的逐点去噪算法^[3].

图信号处理的方法为图像去噪提供了一种新思路. 非局部图变换算法通过图信号处理在图像块中构造出数据自适应的变换域用于去噪, 并在深度图像去噪问题中表现优秀^[4]. Meyer等人利用图拉普拉斯矩阵的低维特征向量进行去噪^[5]. Yan等人提出了GNNLG, 在非局部方法的基础上, 结合图拉普拉斯先验项进行去噪^[6].

目前基于图信号的去噪的研究中, 主要方法都是通过图像块匹配后, 使用图信号的方法对每一个图像块分别进行去噪, 关于图像块的去噪比较类似于传统方法中基于变换域的方法, 本质上是利用图傅里叶变换得到一个新的变换域. 而本文所提出的算法虽然同样基于图信号处理, 但更类似于基于空域方法, 直接针对图像本身设计滤波器得到去噪图像.

在本文的方法中, 首先选择适当的相似性度量, 计算出图像对应的权重矩阵以及拉普拉斯矩阵, 根据所得矩阵的滤波效果, 提出了一个基于图信号的目标函数, 这个目标函数使用权重矩阵来拟合数据保真度项, 并同时加入拉普拉斯先验项以确保恢复信号的平滑性. 我们随后说明了所提出的目标函数的合理性, 并最终使用迭代的方式对优化问题进行快速求解.

本文的其余安排如下: 首先在第2节中介绍准备工作, 然后在第3节中提出并详细描述了我们的算法, 并分别在第4节和第5节中介绍实验和结论.

2 准备工作 2.1 图信号处理的相关概念

根据图信号处理的相关理论, 可以将图像建模为无向加权图 $ G\left( {V, E, W} \right) $ 上的信号^{[7, 8]}. 通过将图像中的每个像素都看作一个顶点, 可以得到顶点集合 $ V $ ; 边缘集合 $ E $ 中的每一条边连接了一对顶点; 每条边 $ \left( {i, j} \right) \in E $ 对应一个非负边缘权重 $ {w_{i, j}} $ , 权重越大代表边缘两端的顶点 $ i, j \in V $ 的相关性越强. 由于 $ G $ 为无向加权图, 因此有 $ {w_{i, j}} = {w_{j, i}} $ , 显然 $ W $ 是一个对称矩阵. 一个 $ 3 \times 3 $ 图像对应的图结构与相应的权重矩阵如图1所示. 图像的像素值被认为是图 $ G $ 顶点上的图信号, $y = {\left[ {{y_1}, {y_2}, \cdots, {y_N}} \right]^{\rm{T}}}$ . 于是, 对于每一个顶点 $ i \in V $ , 本文定义了一个向量来描述位于该顶点的位置和像素强度信息, 即:

$ {r_i} = \left[ {{x_{i, 1}}, {x_{i, 2}}, {y_i}} \right] $

(2)

其中, $ {x_{i, 1}}, {x_{i, 2}} $ 分别代表顶点 $ i $ 在二维平面的空间坐标, $ {y_i} $ 代表顶点 $ i $ 的像素强度.

确定图权重的一种常用方法是使用距离的指数衰减函数, 即:

$ {w_{i, j}} = \exp \left( { - dist\left( {{r_i}, {r_j}} \right)/{h^2}} \right) $

(3)

其中, $ dist\left( {{r_i}, {r_j}} \right) $ 代表顶点 $ {r_i} $ 和 $ {r_j} $ 之间的某种距离度量, $ {h^2} $ 是给定的参数, 这样的选择在直觉上也是合理的, 因为 $ {r_i} $ 与 $ {r_j} $ 距离越接近, 权重越接近1; 反之, 距离越疏远, 权重也就越小. 本文认为较大权重 $ {w_{i, j}} $ 所连接的顶点之间应该具有相似的图信号值.

关于常用的距离度量在表1中给出, 所列举距离度量的直观图形解释如图2所示. 由于经典回归滤波器对感兴趣信号的底层结构缺乏更强的适应性, 双边滤波器在低信噪比情况下无法提供有效性能, 因此在后文中均选择非局部均值滤波器作为权重矩阵中的距离度量 $ dist\left( {{r_i}, {r_j}} \right) $ . 在非局部均值滤波器中, $ {\widehat{y}}_{i} $ 和 $ {\widehat{y}}_{j} $ 分别指以顶点 $ {r}_{i} $ 和 $ {r}_{j} $ 为中心的图像块.

3 基于图滤波的图像去噪算法 3.1 图信号的谱分析

图信号处理中的一个重要算子是图拉普拉斯算子 $ L = D - W $ , 其中, $ D $ 是一个对角矩阵, 它的第 $ i $ 个对角元素等于 $ W $ 第 $ i $ 行的所有元素之和.

对拉普拉斯矩阵进行特征分解, 有 $ Lu = \lambda u $ , $ \lambda $ 为特征值, $ u $ 为相应的特征向量. 由于 $ L $ 的每一行与常向量的内积为0, 假设 $ u $ 为一个常向量, 则 $ Lu = 0u = \vec 0 $ , 另外对于无向图结构, $ L $ 是半正定的实对称矩阵, 于是可以写作 $L = U\Lambda {U^{\rm{T}}}$ , 其中 $ \Lambda = diag\left( {{\lambda _1}, {\lambda _2}, \cdots, {\lambda _N}} \right) $ 是一个由 $ L $ 的特征值所组成的对角矩阵, $ U $ 由特征向量构成, 每一列代表一个特征向量, 并且满足 ${U^{ - 1}} = {U^{\rm{T}}}$ . 根据特征向量的定义, 拉普拉斯矩阵 $ L $ 的每一个特征向量 $ {u_k} $ 都满足 $ L{u_k} = {\lambda _k}{u_k} $ .

接下来说明, 对于图信号 $ {\textit{z}} $ , ${{\textit{z}}^{\rm{T}}}L{\textit{z}}$ 表示信号 $ {\textit{z}} $ 在 $ L $ 所延拓的空间内的平滑度.

由 $ L{u_k} = {\lambda _k}{u_k} $ , 左乘 ${u_k}^{\rm{T}}$ , 得 ${u_k}^{\rm{T}}L{u_k} = {u_k}^{\rm{T}}{\lambda _k}{u_k}$ , 特征向量可以进行常数倍的放缩, 因此经过适当的约定, 令 ${u_k}^{\rm{T}}{u_k} = 1$ , 可以得到:

$ {u_k}^{\rm{T}}{\lambda _k}{u_k} = {\lambda _k}, k = 1, 2, \cdots, N $

(4)

由于拉普拉斯矩阵 $ L $ 的特征向量 $ {u_1}, {u_2}, \cdots, {u_N} $ 可以构成一个完备正交集, 因此可以将图信号 $ {\textit{z}} $ 分解为不同特征向量的系数和. 分解后 $ {\textit{z}} $ 对应的低频分量越多, ${{\textit{z}}^{\rm{T}}}L{\textit{z}}$ 的值越小, 对应的高频分量越多, ${{\textit{z}}^{\rm{T}}}L{\textit{z}}$ 的值也就越大. 在极端情况下, 如果 ${\textit{z}} = {u_1} = \left[ {1, 1, \cdots, 1} \right]$ , 那么此时 ${{\textit{z}}^{\rm{T}}}L{\textit{z}} = 0$ , 因此如果希望优化后的信号更加平滑, 可以在目标函数中加入 ${{\textit{z}}^{\rm{T}}}L{\textit{z}}$ 作为先验项, 这样倾向于得到平滑度更高的 $ {\textit{z}} $ 作为最优解.

换一个角度, 如果将 $ L $ 看作滤波矩阵, 那么它作用在图信号 $ {\textit{z}} $ 时得到 $ L{\textit{z}} $ , $ L $ 的效果相当于一个自适应的高通滤波器, 假如 $ {\textit{z}} $ 为直流向量, 那么滤波后的结果为 $ \vec 0 $ . 基于这样的性质, 通过计算图像所对应的拉普拉斯矩阵, 就可以过滤图像所对应的图信号, 形成自适应滤波.

同理, 由于 $ L = D - W $ , $ W $ 和 $ L $ 的特征向量是相同的, 对应的特征值逆序, 因此 $ W $ 也可以认为是一个自适应低通滤波器.

3.2 问题建模与优化

本文提出下面的优化问题:

$ \hat {\textit{z}} = \arg \mathop {\min }\limits_{\textit{z}} {\left( {y - {\textit{z}}} \right)^{\rm{T}}}W\left( {y - {\textit{z}}} \right) + \eta {{\textit{z}}^{\rm{T}}}L{\textit{z}} $

(5)

其中, 第一项为数据保真度项, 用来衡量观测输入 $ y $ 和期望输出 $ {\textit{z}} $ 之间的差异; 第二项是数据自适应差分项, 有助于获得更平滑的输出 $ {\textit{z}} $ .

对于目标函数 $E\left( {\textit{z}} \right) = {\left( {y - {\textit{z}}} \right)^{\rm{T}}}W\left( {y - {\textit{z}}} \right) + \eta {{\textit{z}}^{\rm{T}}}L{\textit{z}}$ , 关于 $ E\left( {\textit{z}} \right) $ 的优化可以通过计算函数梯度直接得到其最优解的数学表达式.

首先计算 $ E\left( {\textit{z}} \right) $ 的梯度:

$ \nabla E\left( {\textit{z}} \right) = - W\left( {y - {\textit{z}}} \right) + \eta L{\textit{z}} $

(6)

令 $ \nabla E\left( {\hat {\textit{z}}} \right) = - W\left( {y - \hat {\textit{z}}} \right) + \eta L\hat {\textit{z}} = 0 $ , 解得:

$ \hat {\textit{z}} = {\left( {W + \eta L} \right)^{ - 1}}Wy $

(7)

考虑到权重矩阵 $ W \in {R^{N \times N}} $ , 矩阵求逆的时间复杂度过高, 不便于计算, 因此采用梯度下降法对问题进行求解.

由式(6), 迭代更新表达式为:

$ \begin{split} {{\textit{z}}_k} &= {{\textit{z}}_{k - 1}} - {\mu _k}\nabla E\left( {{{\textit{z}}_{k - 1}}} \right) \\ & = {{\textit{z}}_{k - 1}} + {\mu _k}\left[ {W\left( {y - {{\textit{z}}_{k - 1}}} \right) - \eta L{{\textit{z}}_{k - 1}}} \right] \\ \end{split} $

(8)

其中, $ \;{\mu _k} $ 代表梯度下降法每一轮更新的步长.

如果每一轮都要求出最优的 $ {\mu _k} $ , 就会转换为对 $ \;{\mu _k} $ 的优化问题, 显然会使计算量会变大, 于是可以采用Armijo线搜索确定每一轮的步长.

3.3 算法流程

总结归纳后得到算法1.

算法1. 基于图结构滤波的图像去噪算法

1) 非局部均值核计算输入图像 $\scriptstyle y $ 的权重矩阵 $\scriptstyle W $ ;

2) 根据权重矩阵 $\scriptstyle W $ 计算得到拉普拉斯矩阵 $\scriptstyle L $ ;

3) 循环执行下列迭代直到收敛:

$\scriptstyle {{\textit{z}}_k} = {{\textit{z}}_{k - 1}} - {\mu _k}\nabla E\left( {{{\textit{z}}_{k - 1}}} \right)$

4) 得到输出图像近似等于 $\scriptstyle \widehat{\textit{z}}={\left(W+\eta L\right)}^{-1}Wy$ .

算法1中, 步长 $ \;{\mu _k} $ 通过Armijo线搜索确定.

4 实验分析

为了验证本文所提出的方法的性能, 使用图3中的6张灰度图像作为测试图像, 如果要推广至彩色图像去噪, 将算法分别应用于各个通道即可. 在这些图像上添加加性高斯白噪声, 均值为0, 噪声方差从0.2到400, 从而得到噪声图像. 对照算法分别选用NLM^[3], BM3D^[9], WNNM^[10], 这3种算法都是去噪效果较好且比较常用的算法. NLM的基本思想是基于图像块构建图像的逐点估计. WNNM基于加权核范数进行去噪. BM3D基于图像块进行三维滤波. 所有实验都是基于Matlab, 运行在LENOVO 82DN环境下, 所使用的处理器是Intel Core i5-10210U.

其中, $ n $ 为每像素的比特数, 对于灰度图像, $ n=8 $ , ${\textit{MSE}}$ 表示当前图像和参考图像之间的均方误差. ${\textit{PSNR}}$ 的单位是 $ \mathrm{d}\mathrm{B} $ , 该值越大, 代表失真越小.

对于图像恢复效果的衡量, 使用峰值信噪比( ${\textit{PSNR}}$ )作为性能指标^[11]. 定义为:

$ {\textit{PSNR}} = 10{\lg}\frac{{{{\left( {{2^n} - 1} \right)}^2}}}{\textit{MSE}} $

(9)

表2中给出了不同方法的 ${\textit{PSNR}}$ 指标. 噪声方差处于0.2和400之间. 可以看到, 与其他的去噪算法相比, 本文的方法在所有图像中的 ${\textit{PSNR}}$ 指标都是最佳的. BM3D和NLM在各个条件下呈现出相似的去噪水平, 而我们的方法在噪声方差为0.2时比BM3D的表现平均提升了26.8 dB. 可以观察到本文的方法在低噪声水平下具有更明显的优势, 这是由于在本文的目标函数中的数据保真度项通过权重矩阵的拟合, 可以更加自适应地还原出原图像的细节信息.

视觉上的主观比较如图4, 以Goldhill叠加方差为400的噪声为例. 使用本文的方法所生成的图像纹理更加清晰. 相比WNNM, 本文的方法没有造成过度模糊, 相比NLM, 本文的方法显然能更好地去除噪声, 相比BM3D, 本文的方法消除了其中的振铃效应. 因此, 使用本文的方法可以获得更好的细节, 极大抑制了噪声且避免了振铃, 同时纹理信息清晰可见.

5 结论与展望

本文提出了一种基于图信号的图像去噪通用方法. 这个方法基于以下两点分析: (1)将基于NLM滤波计算得到的图结构用于描述图像的内在框架. (2) 基于图像计算的权重矩阵和拉普拉斯矩阵分别可以看作自适应低通和高通滤波器, 可用于表征图信号的特征. 我们创新性地提出了一种基于图信号的目标函数用于去噪, 对于所提出的目标函数, 给出了理论的最优解和迭代逼近算法. 实验表明, 本文所提出的方法在数值指标和视觉效果上都优于BM3D和WNNM等前沿去噪算法.