0 引言

图像去噪作为低层视觉中的经典问题，已经被广泛地研究，但它仍然是一个活跃的课题。在一般情况下，图像去噪的目的是从它的噪声观测y=x+v中来恢复无噪声的图像x，其中v假设为噪声。各种图像去噪方法在过去的几十年中已经被提出，包括基于滤波的方法^[1]、基于混合的方法^[2]、基于全变分的方法^[3]、基于小波和曲波的方法^[4]、基于稀疏表达的方法^[5]和基于非局部自相似性 (Nonlocal Self-Similarity, NSS) 的方法^[6]等。

自然图像往往有许多重复的局部模块，一个局部图像块可以在整个图像找到许多类似的图像块, 这种被称为非局部自相似性先验是图像恢复中很成功的先验。为了利用这一先验知识，代表性的有基于稀疏表达的方案，该方案是把图像块看成从字典中选择几个元素的线性组合^[7]。其中字典可以从现成的词典选择 (例如小波和曲波)，也可以从自然图像块中学习得到。文献[8]提供了字典学习的方法，并成功地扩展到各种图像处理和计算机视觉应用。通过把图像块看作是多元变量向量的样本并考虑到自然图像是非高斯的，Yu等^[9]利用高斯混合模型对图像块进行建模，并分别取得了先进的去噪技术和图像复原技术的成果。非局部均值^[10]和非局部正规化^[11]的方法，与传统的基于局部自相似性的方法相比，极大地提高了图像去噪效果。

虽然非局部自相似性先验在图像去噪中发挥重要的作用，但是目前存在的方法仅仅利用噪声图像的非局部自相似性。如基于3D滤波的块匹配 (Block Matching with 3D filtering, BM3D) 算法^[12]收集噪声图像的非局部自相似块组成的立方体，并对其进行协同滤波的方式实现图像的去噪。Mairal等^[13]提出了同时稀疏编码学习 (Learned Simultaneous Sparse Coding, LSSC) 算法，通过组合稀疏编码来恢复NSS。文献[14]方法中认为NSS小块组成的矩阵奇异值服从拉普拉斯分布。Dong等^[5]把NSS和局部稀疏编码集成到一个框架中，这个框架显示出了强大的图像复原能力。通过假设矩阵的非局部相似的小块有一个低秩的结构，基于低秩最小化的方法^[15-16]也取得了非常有竞争力的去噪效果。然而这些模型都忽视了外部干净图像的非局部自相似先验，因此本文应用混合高斯模型 (Gaussian Mixed Model, GMM) 获取外部干净图像的非局部自相似先验，并将其应用到带权的核范数最小化 (Weighted Nuclear Norm Minimization, WNNM) 图像去噪模型中，提高模型的去噪性能。

1 带权的核范数最小化图像去噪算法

对于一幅图像y的局部分块y_j，本文可以通过像块匹配的方法在图像中找到它的非局部分块 (在实践中是一个足够大的局部窗口)。通过将这些非局部相似分块组成矩阵Y_j，假设X_j和N_j分别是原始图像分块矩阵和噪声图像分块矩阵时，可以得到Y_j=X_j+N_j。根据图像的先验，X_j应当是一个低秩矩阵，低秩矩阵逼近方式可以用于从Y_j估计X_j。聚集所有经过降噪的分块，可以估计整个图像。事实上，核范数最小化 (Nuclear Norm Minimization, NNM) 已经在文献[16]中被应用于视频降噪。文献[15]在图像降噪中将WNNM模型应用于Y_j来估计X_j。通过使用噪声方差σ_n²来归一化F范数的保真项‖Y_j-X_j‖_F², 可以得到以下函数：

${{\mathit{\boldsymbol{\hat X}}}_j} = \mathop {\arg \min }\limits_{{X_j}} \frac{1}{{\sigma _n^2}}\left\| {{\mathit{\boldsymbol{Y}}_j} - {\mathit{\boldsymbol{X}}_j}} \right\|_{\rm{F}}^2 + {\left\| {{\mathit{\boldsymbol{X}}_j}} \right\|_{w,*}}$

(1)

显然，该模型的主要问题是权向量w的确定。权重向量与对应的奇异值成反比^[17]，于是得到如下的权重公式：

${w_i} = c\sqrt n /\left( {{\sigma _i}({\mathit{\boldsymbol{X}}_j}) + \varepsilon } \right)$

(2)

这里σ_i(X_j) 是图像块组X_j对应的第i个奇异值，c是一个常数，n是Y_j中相似的分块的数量，ε=10^-16是为了避免分母为0。由于算法中σ_i(X_j) 的奇异值不能直接获得。

σ_i(X_j) 可以近似估计为:

$\hat \sigma ({\mathit{\boldsymbol{X}}_j}) = \sqrt {\max ({\sigma ^2}({\mathit{\boldsymbol{Y}}_j}) - n\sigma _n^2,0)} $

(3)

其中：σ_i(Y_j) 是Y_i的第i个奇异值。注意，当$\hat \sigma ({\mathit{\boldsymbol{X}}_j})$按非递增顺序排列，因此获得的权重w_i=1，2, …，n被确保为非递增序列。通过将上述步骤应用到每个分块并将所有分块组合，图像X就被复原了，具体细节如算法1所示。

算法1。

1)  对于噪声图像y，初始化${\mathit{\boldsymbol{\hat x}}}$⁽⁰⁾=y, y⁽⁰⁾=y，给定常量L, δ

2)  k=1:L, 执行

3)  $\mathit{\boldsymbol{y}}\left( k \right)\mathit{\boldsymbol{ = }}{{\mathit{\boldsymbol{\hat x}}}^{\left( {k - 1} \right)}} + \delta \left( {\mathit{\boldsymbol{y}} - {{\mathit{\boldsymbol{\hat y}}}^{\left( {k - 1} \right)}}} \right)$

4)  对y^(k)每一小块y_j执行

  {寻找相似分块集合Y_j

  奇异值分解[U, Σ, V]=SVD (Y_j)

  根据式 (2) 计算估计权向量w

  获得估计：${{\mathit{\boldsymbol{\hat X}}}_j}$=US_w(Σ)V^T}

5)  聚集${{\mathit{\boldsymbol{\hat X}}}_j}$来组成去噪图像${{\mathit{\boldsymbol{\hat x}}}^{\left( k \right)}}$，跳转到步骤2)

算法1中S_w(Σ)_ii被定义为:

${S_w}{\left( \mathit{\pmb{\Sigma}} \right)_{ii}} = \max \left( {{\mathit{\pmb{\Sigma}}_{ii}} - {w_i},0} \right)$

(4)

2 基于GMM和WNNM的混合去噪模型

自然图像通常有一些重复的局部模式，并且对于一个局部的图像块在整个图像中可以找到许多与之相似的图像块，这些图像块组成的矩阵通常是低秩的，这个被称为非局部自相似的先验是图像恢复最成功的先验之一。

带权核范数最小化图像去噪模型收集欧氏距离相近的图像块向量组成非局部自相似矩阵，这个非局部自相似矩阵的质量好坏直接影响去噪性能。当图像受到噪声污染，根据退化图像的欧氏距离来组成非局部自相似矩阵是存在问题的，因为噪声会影响非局部自相似矩阵的组成，进而影响该矩阵的低秩性，导致NSS先验失败。

在图像去噪中仅利用退化图像的非局部自相似性先验是不够的，可以利用外部的无噪声自然图像的非局部自相似性的先验。为此本文对无噪声的自然图像的非局部自相似性图像块进行训练产生外部NSS先验的GMM，然后用该GMM指导退化图像生成非局部自相似块矩阵，最后用WNNM算法进行图像的去噪。这种结合退化图像本身的NSS先验和外部干净图像NSS先验得到的非局部自相似块矩阵质量更高，有助于提升图像的去噪效果，这点也通过后面的仿真实验得以验证。

可以假设自然图像块源于一个K维子空间的潜在结构，所以一个给定的图像块x_i的概率可以被定义为K个高斯的加权和，如式 (5) 所示：

$p(\left. {{\mathit{\boldsymbol{x}}_i}} \right|\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}) = \sum\limits_{k = 1}^k {{\omega _k}{p_k}(\left. {{\mathit{\boldsymbol{x}}_i}} \right|{\mathit{\boldsymbol{\theta }}_k})} $

(5)

这里Θ=(ω₁, ω₂, …, ω_k, θ₁, θ₂, …, θ_K) 是一组参数且$\sum\limits_{k = 1}^K {{\omega _k} = 1} $。注意，每一个θ_k描述了平均值为μ_k、协方差矩阵为Σ_k的高斯密度函数p_k，且p_k(x_i|θ_k)=c*exp (0.5*(x_i-μ_k)^TΣ_k^-1(x_i-μ_k))。参数Θ可以由最大似然估计确定，通常利用的是期望最大化 (Expectation Maximization, EM) 算法。

引入类标签变量C∈{1, 2, …, K}。对于每一块x_i, 本文用GMM学习的概率密度函数 (Probability Density Function, PDF) 来判断x_i属于哪一类C_i，式 (6) 表示属于第k类的概率：

$p(k|{\mathit{\boldsymbol{x}}_i}{\rm{;}}\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}) = \frac{{{\omega _k}{p_k}({\mathit{\boldsymbol{x}}_i}{\rm{|}}{\theta _k})}}{{\sum\limits_{j = 1}^K {{\omega _j}} {p_j}({\mathit{\boldsymbol{x}}_i}|{\theta _j})}};k = 1,2, \cdots ,K$

(6)

根据式 (6) 找出其中的最大值赋值给C_i；然后，对于给定图像块x_i，本文根据式 (7) 的马氏距离来寻找与它相似的非局部分块，形成非局部相似矩阵Y_j，最后应用WNNM模型对含噪声图像进行去噪，详细步骤见算法2。

${D_{i,j}} = {\left( {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{x}}_j}} \right)^{\rm{T}}}\mathit{\pmb{\Sigma}}_{{C_i}}^{ - 1}\left( {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{x}}_j}} \right)$

(7)

算法2。

1) 选定K, 对自然无噪声的图像块用GMM方法进行训练，得到参数Θ

2) 对于噪声图像y，初始化${\mathit{\boldsymbol{\hat x}}}$⁽⁰⁾ =y, y⁽⁰⁾ =y, L, δ和噪声方差σ_n²

3) k=1:L, 执行

4) $\mathit{\boldsymbol{y}}\left( k \right)\mathit{\boldsymbol{ = }}{{\mathit{\boldsymbol{\hat x}}}^{\left( {k - 1} \right)}} + \delta \left( {\mathit{\boldsymbol{y}} - {{\mathit{\boldsymbol{\hat y}}}^{\left( {k - 1} \right)}}} \right)$

5) 对y^(k)每一小块y_j根据式 (6) 计算C_i

根据式 (7) 寻找相似分块集合Y_j，根据式 (2) 估计权向量w

奇异值分解[U, Σ, V]=SVD (Y_j)

获得估计：${{\mathit{\boldsymbol{\hat X}}}_j}$=US_w(Σ)V^T

6) 聚集${{\mathit{\boldsymbol{\hat X}}}_j}$来组成去噪图像${\mathit{\boldsymbol{\hat x}}}$^(k)返回步骤3)

3 模型的扩展及收敛性分析

上述方法解决的式 (1) 所示的去噪模型，该模型中的保真项，适用于高斯噪声的特点，因此对高斯噪声有很好的抑制效果，然而当对均匀噪声等其他噪声或者图像中还有粗差的效果不好，因此需要对式 (1) 所示的去噪模型的优化目标进行修正为：

$F\left( {{\mathit{\boldsymbol{X}}_j}} \right){\rm{ = }}h\left( {{\mathit{\boldsymbol{X}}_j}} \right) + {\left\| {{\mathit{\boldsymbol{X}}_j}} \right\|_{w,*}}$

(8)

其中h(·) 为在R^m×n→R⁺利普希茨连续可导函数，利普希茨常数为L(h)，h(·) 替换式 (1) 中的保真项更具一般性，对模型 (8) 求最小值可以利用下面的更新策略求解：

$\begin{array}{*{20}{l}} {{{\bf{X}}^{k + 1}}}&{ = \mathop {{\rm{arg}}{\mkern 1mu} {\rm{min}}}\limits_\mathit{\boldsymbol{X}} \sum\limits_{i = 1}^m {{{\left\| \mathit{\boldsymbol{X}} \right\|}_{w,*}}} + h\left( {{\mathit{\boldsymbol{X}}^k}} \right) + }\\ {}&{\left\langle {\nabla h\left( {{\mathit{\boldsymbol{X}}^k}} \right),\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}^k}} \right\rangle + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}^k}} \right\|_{\rm{F}}^2 = }\\ {}&{ = \mathop {{\rm{arg}}{\mkern 1mu} {\rm{min}}}\limits_\mathit{\boldsymbol{X}} \sum\limits_{i = 1}^m {{{\left\| \mathit{\boldsymbol{X}} \right\|}_{w,*}} + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}^k} + \frac{1}{\mu }\nabla h\left( {{\mathit{\boldsymbol{X}}^k}} \right)} \right\|_{\rm{F}}^2} } \end{array}$

(9)

其中μ>L(h)。

收敛性分析：

为了证明上述扩展模型的收敛性，需要证明当μ>L(h) 则数列由式 (9) 获得的序列{X^k}满足以下特性：

1) F{X^k}是单调递减的，即

$F({\mathit{\boldsymbol{X}}^k}) - F({\mathit{\boldsymbol{X}}^{k + 1}}) \ge \frac{{\mu - L(h)}}{2}\left\| {{\mathit{\boldsymbol{X}}^{\rm{k}}} - {\mathit{\boldsymbol{X}}^{k + 1}}} \right\|_{\rm{F}}^2 \ge 0$

2) $\mathop {\lim }\limits_{k \to \infty } \left( {{\mathit{\boldsymbol{X}}^k} - {\mathit{\boldsymbol{X}}^{k + 1}}} \right) = 0$

证明 因为X^k+1是式 (9) 的最优解，可以得到式 (10) :

$\begin{array}{l} \sum\limits_{i = 1}^m {{{\left\| {{\mathit{\boldsymbol{X}}^{k + 1}}} \right\|}_{w,*}} + h\left( {{\mathit{\boldsymbol{X}}^k}} \right) + \left\langle {\nabla h\left( {{\mathit{\boldsymbol{X}}^k}} \right),{\mathit{\boldsymbol{X}}^{k + 1}} - {\mathit{\boldsymbol{X}}^k}} \right\rangle + } \\ \quad \frac{\mu }{2}\left\| {{\mathit{\boldsymbol{X}}^{k + 1}} - {\mathit{\boldsymbol{X}}^k}} \right\|_{\rm{F}}^2 \le \sum\limits_{i = 1}^m {{{\left\| {{\mathit{\boldsymbol{X}}^k}} \right\|}_{w,*}} + h\left( {{\mathit{\boldsymbol{X}}^k}} \right)} + \\ \quad \left\langle {\nabla h\left( {{\mathit{\boldsymbol{X}}^k}} \right),{\mathit{\boldsymbol{X}}^k} - {\mathit{\boldsymbol{X}}^k}} \right\rangle + \frac{\mu }{2}\left\| {{\mathit{\boldsymbol{X}}^k} - {\mathit{\boldsymbol{X}}^k}} \right\|_{\rm{F}}^2 = \\ \quad \sum\limits_{i = 1}^m {{{\left\| {{\mathit{\boldsymbol{X}}^k}} \right\|}_{w,*}}} \end{array}$

(10)

另一方面，因为h在利普希茨连续可导，所以可以得到:

$\begin{array}{l} h\left( {{\mathit{\boldsymbol{X}}^{k + 1}}} \right) \le h\left( {{\mathit{\boldsymbol{X}}^k}} \right) + \left\langle {\nabla h\left( {{\mathit{\boldsymbol{X}}^k}} \right),{\mathit{\boldsymbol{X}}^k} - {\mathit{\boldsymbol{X}}^k}} \right\rangle + \\ \quad \frac{{L\left( h \right)}}{2}\left\| {{\mathit{\boldsymbol{X}}^{k + 1}} - {\mathit{\boldsymbol{X}}^k}} \right\|_{\rm{F}}^2 \end{array}$

(11)

根据式 (10) 和 (11) 可以得到:

$\begin{array}{l} F\left( {{\mathit{\boldsymbol{X}}^k}} \right) - F\left( {{\mathit{\boldsymbol{X}}^{k + 1}}} \right) = \sum\limits_{i = 1}^m {{{\left\| {{\mathit{\boldsymbol{X}}^k}} \right\|}_{w,*}}} + h\left( {{\mathit{\boldsymbol{X}}^k}} \right) - \\ \quad \sum\limits_{i = 1}^m {{{\left\| {{\mathit{\boldsymbol{X}}^{k + 1}}} \right\|}_{w,*}}} - h\left( {{\mathit{\boldsymbol{X}}^{k + 1}}} \right) \ge \\ \quad \frac{{\mu - L\left( h \right)}}{2}\left\| {{\mathit{\boldsymbol{X}}^{k + 1}} - {\mathit{\boldsymbol{X}}^k}} \right\|_{\rm{F}}^2 \end{array}$

(12)

由于μ>L(h)，所以F(X^k)≥F(X^k+1)。

将式 (12) 对k=1, 2, 3，…进行累加可得到：

$F({\mathit{\boldsymbol{X}}^1}) \ge \frac{{\mu - L(h)}}{2}\sum\limits_{k = 1}^{ + \infty } {\left\| {{\mathit{\boldsymbol{X}}^{k + 1}} - {\mathit{\boldsymbol{X}}^k}} \right\|_{\rm{F}}^2} $

(13)

进而得到：

$\mathop {\lim }\limits_{k \to \infty } {\rm{(}}{\mathit{\boldsymbol{X}}^k} - {\mathit{\boldsymbol{X}}^{k + 1}}) = 0$

(14)

从而得到模型的敛散性。

4 仿真实验

为了验证上面提出算法的有效性，本文将算法应用到实际的例子中，并将本文算法跟其他经典的算法 (包括BM3D^[12]、LSSC^[13]和WNNM^[15]) 作比较。本实验的计算机为普通PC，计算机的配置为CPU i7-4770k, 内存为12 GB。为了使本文的仿真实验更加逼近真实情况，文中添加的噪声为不同噪声类型的叠加。本文训练GMM的样本来自Berkeley数据集，从这些图像中提取了2×10⁶个图像块。图像块的尺寸和很多文献一样设置为7×7、8×8、9×9和10×10，分别针对方差在[0, 20]、[20, 40]、[40, 60]和[60, 100]四种情况。

下面本文将就算法性能 (文中采用峰值信噪比 (Peak Signal-to-Noise Ratio, PSNR) 作为评价指标) 和算法的运行效率两方面进行比较说明。

限于篇幅，图 1和图 2分别给出了C.man和Monarch这两张图片在噪声方差为30、50，本文算法与其他先进的算法去噪的直观效果。

表 1显示了在噪声方差为30、50、100的情况下PSNR性能的比较。从实验结果可以看出BM3D性能平均提高了0.49 dB, 相对于LSSC算法提高了0.38 dB, 相对于WNNM算法平均提高了0.11 dB。本文的算法相对于BM3D算法和LSSC算法有显著的提高，相对于WNNM算法也有明显的提升。这些提升也验证了本文用混合高斯模型学习外部无噪声的图像的非局部自相似性这一先验知识，并在这一先验知识的指导下采用马氏距离来形成非局部自相似图像块矩阵，代替其他算法采用欧氏距离形成非局部自相似图像块矩阵更加合理。

此外，需要注意的是，本文算法在大噪声的情况下相对于其他算法优势更明显，在噪声方差为100时，本文算法相对于BM3D算法性能平均提高了0.65 dB, 相对于LSSC算法提高了0.59 dB, 相对于WNNM算法平均提高了0.156 dB。本文方法随着噪声方差的增加去噪优势更加明显。

表 2显示了本文算法与其他先进算法处理256×256图像的时间运行时间的比较，算法的实验结果是在不同的噪声方差下，对不同图像运行100次的平均时间。从表 4可以看出，虽然本文算法运行时间与WNNM相当，高于BM3D算法，但相对于WNNM来说牺牲很少的运行效率却获得了更好的去噪效果。

5 结语

本文首先用无噪声的自然图像块训练混合高斯模型; 然后再用训练好的混合高斯模型指导采集非局部自相似块，用马氏距离代替欧氏距离来衡量图像块之间的相似性，并为了使模型更具一般性，提出了扩展模型; 最后，为了验证本文算法的有效性，与一些先进的算法进行仿真实验，仿真结果表明，所提算法在很多情况下具有优势。