图像去噪的目的是在消除噪声的同时最大程度地保留图像细节信息。图像在采集、传输和存储过程中经常受到噪声污染，导致图像质量降低，甚至影响图像增强(李健等，2021)、识别(贺敏雪等，2021)、分类(尹红等，2019)和分割(江宗康等，2020)等后续处理工作。因此，图像去噪是一项非常重要的图像预处理工作。

滤除高斯噪声，常见的空间域算法是基于变分法的去噪模型，可以通过最小化能量函数达到平滑图像的目的。Rudin等人(1992)提出一种非线性滤波算法，即全变分算法(total variation，TV)，可以通过减小噪声图像的绝对梯度的积分来滤除噪声。Goldstein和Osher(2009)利用分裂Bregman快速算法提出了各向异性全变分和各项同性全变分。Adam和Paramesran(2019)为了更好地平衡消除阶梯效应和保留边缘信息之间的关系，提出一种结合重叠群稀疏和混合非凸高阶全变分模型(hybrid non-convex higher order TV with overlapping group sparsity，HNHOTV-OGS)。Parisott等人(2020)设计了一种高阶方向性全变分(total direction variation，TDV)，在高阶梯度中加入各项异性方向信息。然而，由于全变分完全依赖于图像的局部特征，在边缘部分或梯度较大处容易受到噪声干扰而丢失细节信息，从而产生阶梯效应。

针对局部滤波这一缺点，Buades等人(2005)提出了非局部均值算法(non-local means，NLM)，主要思想是将图像块作为一个处理单元，计算中心图像块与所有邻域图像块之间的相似性权重，最后得到的去噪图像块即邻域图像块的加权平均值。在NLM算法取得显著成效后，人们对基于图像块的非局部去噪算法进行改进。例如，将全变分与非局部相似性相结合，Gilboa和Osher(2007)提出了非局部全变分(non-local total variation，NLTV)，Liu等人(2014)设计了非局部广义相对全变分(non-local version of the generalized relative TV，NLGRTV)，Li等人(2017)设计了正则化非局部全变分(regularized non-local total variation，RNLTV)，Wang等人(2019)提出了基于结构相似性的非局部全变分。不同于局部全变分模型仅通过图像像素灰度值的变化来除噪，这类非局部全变分算法利用了邻域内图像块结构的相似性，可以更好地保留图像的细节信息。

值得注意的是，现有许多去噪算法原本是针对灰度图像提出的，在处理彩色图像时只是简单地对3个颜色分量进行单独去噪。此外，大部分去噪算法忽略了真实噪声在图像不同区域内的分布差异，而是简单假定整幅图像的噪声强度是均匀的，并设置同样的参数，导致这些算法在每个颜色分量和每个图像块上，甚至每次迭代中的去噪强度都是一样的，从而造成过平滑和欠平滑共存现象。

为了解决以上问题，本文提出一种自适应非局部3维全变分(adaptive non-local 3D total variation，ANL3DTV)去噪模型，利用一种非局部3维全变分正则项来捕获图像块内和块间的相关性信息，同时引入一个自适应权重矩阵，根据每次迭代后去噪图像中残留噪声的特征，调整算法在不同颜色分量和不同图像块上的去噪强度。本文从数学上对ANL3DTV模型进行推导求解，分析该模型的收敛性，最后通过仿真去噪实验验证了算法的有效性。

给定噪声图像${\boldsymbol{Y}}={\boldsymbol{\tilde X }}+{\boldsymbol{N}}$，其中${\boldsymbol{\tilde X }}$是对应的原始图像或未加噪图像，${\boldsymbol{N}}$是高斯噪声，去噪图像${\boldsymbol{X}}$为

$\boldsymbol{X}_i=\sum\limits_{j \in S_i} \omega_{i j} \boldsymbol{Y}_j $

(1)

$ \omega_{i j}=\frac{1}{W} \cdot \exp \left(-\frac{\left\|\boldsymbol{Y}_i-\boldsymbol{Y}_j\right\|_{\mathrm{F}}^2}{h^2}\right) $

(2)

式中，归一化参数$W=\sum\limits_{j \in \boldsymbol{S}_i} \exp \left(-\left\|\boldsymbol{Y}_i-\boldsymbol{Y}_j\right\|_{\mathrm{F}}^2 / h^2\right)$。${\boldsymbol{X}}_{i}$表示去噪图像中的第$i$个图像块，${\boldsymbol{Y}}_{i}$和${\boldsymbol{Y}}_{j}$分别表示噪声图像中的第$i$个和第$j$个图像块，$ω_{ij}$表示图像块${\boldsymbol{Y}}_{i}$和${\boldsymbol{Y}}_{j}$间的相似性权重，${\boldsymbol{S}}_{i}$是中心像素位于$i$的搜索窗口，$h$是平滑因子，$||·||_\text{F}$表示弗罗贝尼乌斯范数。

全变分去噪模型(Rudin等，1992)可以表示为

$\min\limits _{\boldsymbol{X}} \frac{1}{2}\|\boldsymbol{Y}-\boldsymbol{X}\|_{\mathrm{F}}^2+\lambda\|\nabla \boldsymbol{X}\|_{\ell} $

(3)

式中，${\boldsymbol{Y}}$和${\boldsymbol{X}}$分别为加噪信号和去噪信号，$\nabla$是梯度算子，$λ>0$是正则化参数。$||\nabla {\boldsymbol{X}}||_{1}$称为各向异性全变分(Yan和Lu，2015)，$||\nabla {\boldsymbol{X}}||_{2}$称为各向同性全变分(Rudin等，1992)。两种全变分计算为

$ \|\nabla \boldsymbol{X}\|_1=\left\|\nabla_x \boldsymbol{X}\right\|_1+\left\|\nabla_y \boldsymbol{X}\right\|_1 $

(4)

$ \|\nabla \boldsymbol{X}\|_2=\sqrt{\left\|\nabla_x \boldsymbol{X}\right\|_2^2+\left\|\nabla_y \boldsymbol{X}\right\|_2^2} $

(5)

式中，$\nabla _{x}$和$\nabla _{y}$分别用来计算2维信号在$X$方向和$Y$方向上的梯度，$||·||_{1}$和$||·||_{2}$分别表示1范数和2范数。

已知噪声图像${\boldsymbol{Y}}∈{\mathbf{R}}^{M×N}$，即${\boldsymbol{Y}}$是大小为$M×N$的实数矩阵，非局部全变分的去噪模型可以表示(Gilboa和Osher，2009；Lou等，2010)为

$ \min\limits _{\boldsymbol{X}} \frac{1}{2}\|\boldsymbol{Y}-\boldsymbol{X}\|_{\mathrm{F}}^2+\lambda\left\|\nabla_\omega \boldsymbol{X}\right\|_{\mathrm{F}} $

(6)

$ \left\|\nabla_\omega \boldsymbol{X}\right\|_{\mathrm{F}}=\sum\limits_{i \in \boldsymbol{\varOmega}} \sqrt{\sum\limits_{j \in \boldsymbol{S}_i}\left(x_i-x_j\right)^2 \omega_{i j}} $

(7)

$ \omega_{i j}=\frac{1}{W} \cdot \exp \left(-\frac{\left\|\boldsymbol{X}_i-\boldsymbol{X}_j\right\|_{\mathrm{F}}^2}{h^2}\right) $

(8)

$ W=\sum\limits_{j \in \boldsymbol{S}_i} \exp \left(-\left\|\boldsymbol{X}_i-\boldsymbol{X}_j\right\|_{\mathrm{F}}^2 / h^2\right) $

(9)

式中，$x_{i}$和$x_{j}$分别表示${\boldsymbol{X}}$的第$i$个和第$j$个像素值，$\boldsymbol{\varOmega}=\{1, 2, \cdots, M\} \times\{1, 2, \cdots, N\}$。

ANL3DTV算法的主要思想是通过计算彩色图像块内的局部梯度和图像块间的非局部梯度，并利用不同图像块和颜色分量上的噪声强度差异，最终实现自适应去噪。ANL3DTV去噪算法流程如图 1所示，具体如下：

1) 将给定的噪声图像${\boldsymbol{Y}}$分成$K$个相互重叠的图像块，每个图像块大小为$p×p×3$，然后依次以每个图像块为中心图像块，在对应的搜索窗内找到前$m$个与中心图像块最相似的邻域图像块，并构成图像块组。

2) 将图像块组中的每个彩色图像块拉直成长度为3$p^{2}$列向量，得到一个2维矩阵${\boldsymbol{Y}}_{i}∈{\mathbf{R}}^{3p^{2}×m}$。

3) 对${\boldsymbol{Y}}_{i}$进行去噪，得到去噪后的2维矩阵${\boldsymbol{X}}_{i}∈{\mathbf{R}}^{3p^{2}×m}$。

4) 将所有去噪后的${\boldsymbol{X}}_{i}$返回至原始位置，并对重叠区域进行加权平均，得到最终的去噪图像${\boldsymbol{X}}$。

步骤3)中，对${\boldsymbol{Y}}_{i}$进行去噪的ANL3DTV算法的去噪模型为

$ \min\limits _{\dot{\boldsymbol{X}}_i, \boldsymbol{W}_i} \frac{1}{2}\left\|\boldsymbol{W}_i \circ\left(\boldsymbol{Y}_i-\boldsymbol{X}_i\right)\right\|_{\mathrm{F}}^2+\lambda\left\|\nabla \boldsymbol{X}_i\right\|_1 $

(10)

式中，$λ>0$为正则化参数，$ \circ $表示矩阵点乘。$ \nabla {\boldsymbol{X}} _{i}$=[${\boldsymbol{D}}_{x}{\boldsymbol{X}}_{i}, {\boldsymbol{D}}_{y}{\boldsymbol{X}}_{i}, {\boldsymbol{X}}_{i}{\boldsymbol{D}}_{z}]∈{\mathbf{R}}^{3p^{2}×3 m}$是${\boldsymbol{X}}_{i}$的梯度矩阵，${\boldsymbol{D}}_{x}∈{\mathbf{R}}^{3p^{2}×3p^{2}}$、${\boldsymbol{D}}_{y}∈{\mathbf{R}}^{3p^{2}×3p^{2}}$和${\boldsymbol{D}}_{z}∈{\mathbf{R}}^{m×m}$是梯度算子。其中，${\boldsymbol{D}}_{x}{\boldsymbol{X}}_{i}$和${\boldsymbol{D}}_{y}{\boldsymbol{X}}_{i}$用来计算图像块内的水平和垂直梯度，即局部像素间的相似性；${\boldsymbol{X}}_{i}{\boldsymbol{D}}_{z}$用来计算相似图像块间的差值，即非局部像素间的相似性。此外，${\boldsymbol{W}}_{i}∈{\mathbf{R}}^{3p^{2}×m}$是一个自适应权重矩阵，计算为

$\boldsymbol{w}_{i j}=\left[w_{i j, R} \boldsymbol{I}_{p^2 \times 1} ; w_{i j, G} \boldsymbol{I}_{p^2 \times 1} ; w_{i j, B} \boldsymbol{I}_{p^2 \times 1}\right] $

(11)

$ w_{i j, c}=\frac{1}{\sigma_{i j, c}+\varepsilon}, c \in\{\mathrm{R}, \mathrm{G}, \mathrm{B}\} $

(12)

$\sigma_{i j, c}=\sqrt{\max \left(0, \sigma_c^2-\left\|\boldsymbol{y}_{i j, c}-\boldsymbol{x}_{i j, c}\right\|_2^2 / p^2\right)} $

(13)

式中，${\boldsymbol{y}}_{ij}$=[${\boldsymbol{y}}_{ij, R}, {\boldsymbol{y}}_{ij, G}, {\boldsymbol{y}}_{ij, B}$]、${\boldsymbol{x}}_{ij}$=[${\boldsymbol{x}}_{ij, R}, {\boldsymbol{x}}_{ij, G}, {\boldsymbol{x}}_{ij, B}$]和${\boldsymbol{w}}_{ij}$分别为${\boldsymbol{Y}}_{i}$、${\boldsymbol{X}}_{i}$和${\boldsymbol{W}}_{i}$的第$j$个列向量，${\boldsymbol{y}}_{ij, c}, {\boldsymbol{x}}_{ij, c}, {\boldsymbol{I}}_{p^{2}×1}∈{\mathbf{R}}^{p^{2}×1}$，${\boldsymbol{I}}_{p^{2}×1}$是一个全1列向量，[$σ_\text{R}$, $σ_\text{G}$, $σ_\text{B}$]是彩色图像的3个颜色分量上的噪声强度，这里为噪声标准差。

为了求解ANL3DTV去噪模型，引入辅助变量${\boldsymbol{H}}$，得到式(10)的增广拉格朗日函数，具体为

$\min\limits _{\varTheta=\left\{\boldsymbol{X}_i, \boldsymbol{H}\right\}, \varPhi, \beta} \mathfrak{J}(\boldsymbol{\varTheta} ; \boldsymbol{\varPhi}, \boldsymbol{\beta}) $

(14)

式中

$\begin{gathered} \mathfrak{J}(\boldsymbol{\varTheta} ; \boldsymbol{\varPhi}, \beta)=\frac{1}{2}\left\|\boldsymbol{W}_i \circ\left(\boldsymbol{Y}_i-\boldsymbol{X}_i\right)\right\|_{\mathrm{F}}^2+\lambda\|\boldsymbol{H}\|_1+ \\ \left\langle\boldsymbol{\varPhi}, \boldsymbol{H}-\nabla \boldsymbol{X}_i\right\rangle+\frac{\beta}{2}\left\|\boldsymbol{H}-\nabla \boldsymbol{X}_i\right\|_{\mathrm{F}}^2 \end{gathered} $

(15)

式中，${\boldsymbol{\varPhi }}$是增广拉格朗日乘子，$β$是惩罚因子。然后，利用交替方向乘子法求解式(12)，即先优化一个变量，同时固定其他变量，交替进行，直至收敛。初始化相关变量，第$t+1$次迭代计算过程如下：

1) 求解${\boldsymbol{X}}_{i}$。

关于${\boldsymbol{X}}_{i}$的优化模型为

$\min \limits_{\boldsymbol{x}_i}\left\|\hat{\boldsymbol{W}}_i\left(\boldsymbol{y}_i-\boldsymbol{x}_i\right)\right\|_{\mathrm{F}}^2+\left\|\boldsymbol{h}^{(t)}-\boldsymbol{D} \boldsymbol{x}_i+\boldsymbol{\phi}^{(t)} / \beta^{(t)}\right\|_{\mathrm{F}}^2 $

(16)

$ \begin{aligned} & \boldsymbol{x}_i^{(t+1)}=\left[\hat{\boldsymbol{W}}_i^2+\beta^{(t)} \boldsymbol{D}^{\mathrm{T}} \boldsymbol{D}\right]^{-1} \cdot \\ & {\left[\hat{\boldsymbol{W}}_i^2 \boldsymbol{y}_i+\boldsymbol{D}^{\mathrm{T}}\left(\boldsymbol{\beta}^{(t)} \boldsymbol{h}^{(t)}+\boldsymbol{\phi}^{(t)}\right)\right]} \end{aligned} $

(17)

式中，${\boldsymbol{y}}_{i}$、${\boldsymbol{x}}_{i}$、${\boldsymbol{h}}$和${\boldsymbol{\phi }}$分别是矩阵${\boldsymbol{Y}}_{i}$、${\boldsymbol{X}}_{i}$、${\boldsymbol{H}}$和${\boldsymbol{\varPhi }}$的向量形式，$\boldsymbol{D}=\left[\begin{array}{l} \boldsymbol{I}_m \otimes \boldsymbol{D}_x \\ \boldsymbol{I}_m \otimes \boldsymbol{D}_y \\ \boldsymbol{D}_z \otimes \boldsymbol{I}_{p^2} \end{array}\right] \in \mathbf{R}^{9 m p^2 \times 3 m p^2}$是梯度算子，${\boldsymbol{I}}_{m}∈{\mathbf{R}}^{m×m}$和${\boldsymbol{I}}_{p^{2}}∈{\mathbf{R}}^{p^{2}×p^{2}}$是单位矩阵，$⊗$是克罗内克乘积算子，$\hat{\boldsymbol{W}}_{i}=\text{diag}(vec({\boldsymbol{W}}_{i}))∈{\mathbf{R}}^{3mp^{2}×3mp^{2}}$，$\text{diag}(·)$用来构造对角矩阵，$vec(·)$是向量化算子。

2) 求解${\boldsymbol{H}}$。

关于${\boldsymbol{H}}$的优化模型为

$\min\limits _{\boldsymbol{H}} \frac{\lambda}{\beta^{(t)}}\|\boldsymbol{H}\|_1+\frac{1}{2} \| \boldsymbol{H}-\nabla \boldsymbol{X}_i^{(t+1)}+\boldsymbol{\varPhi}^{(t)} /\left.\beta^{(t)}\right|_{\mathrm{F}} ^2 $

(18)

计算为

$\boldsymbol{H}^{(t+1)}=S_{\lambda / \beta^{(t)}}\left(\beta^{(t)} \nabla \boldsymbol{X}_i^{(t+1)}-\boldsymbol{\varPhi}^{(t)}\right) $

(19)

式中，$S_{α}(b)=\text{sign}(b)·\max(|b|-α, 0)$, sign为符号函数。

3) 求解${\boldsymbol{\varPhi }}$。

$ \boldsymbol{\varPhi}^{(t+1)}=\boldsymbol{\varPhi}^{(t)}+\beta^{(t)}\left(\boldsymbol{H}^{(t+1)}-\nabla \boldsymbol{X}_i^{(t+1)}\right) $

(20)

4) 求解$β$。

$ \beta^{(t+1)}=\tau \beta^{(t)} $

(21)

定理1  式(12)中，$\{{\boldsymbol{X}}^{(t)}_{i}\}$和$\{{\boldsymbol{H}}^{(t)}\}$是收敛序列，且满足

$\lim\limits _{t \rightarrow \infty}\left\|\boldsymbol{H}^{(t+1)}-\nabla \boldsymbol{X}_i^{(t+1)}\right\|_{\mathrm{F}}=0 $

(22)

证明  首先，利用式(18)和式(19)得到

$ \begin{gathered} \left\|\boldsymbol{\varPhi}^{(t+1)}\right\|_{\mathrm{F}}=\left\|\boldsymbol{\varPhi}^{(t)}+\beta^{(t)}\left(\boldsymbol{H}^{(t+1)}-\nabla \boldsymbol{X}_i^{(t+1)}\right)\right\|_{\mathrm{F}}= \\ \beta^{(t)} \| S_{\lambda / \beta^{(t)}}\left(\beta^{(t)} \nabla \boldsymbol{X}_i^{(t+1)}-\boldsymbol{\varPhi}^{(t)}\right)- \\ \left(\beta^{(t)} \nabla \boldsymbol{X}_i^{(t+1)}-\boldsymbol{\varPhi}^{(t)}\right) \|_{\mathrm{F}} \leqslant \\ \beta^{(t)} \cdot \lambda / \beta^{(t)} \cdot \sqrt{9 m p^2}=3 p \lambda \sqrt{m} \end{gathered} $

(23)

所以，$\{{\boldsymbol{\varPhi }}^{(t)}\}$是有上界的，从而可以推出$\{{\boldsymbol{\varPhi }}^{(t+1)}-{\boldsymbol{\varPhi }}^{(t)}\}$也是有上界的。

因为${\boldsymbol{\varTheta}}^{(t+1)}$是函数$\mathfrak{J}({\boldsymbol{\varTheta}}$; ${\boldsymbol{\varPhi }}^{(t)}, β^{(t)})$的全局最优解，所以有$\mathfrak{J}({\boldsymbol{\varTheta}}^{(t+1)}$; ${\boldsymbol{\varPhi }}^{(t)}, β^{(t)}$)≤$\mathfrak{J}({\boldsymbol{\varTheta}}^{(t)}$; ${\boldsymbol{\varPhi }}^{(t)}, β^{(t)})$。

令$a$为序列$\{{\boldsymbol{\varPhi }}^{(t+1)}$－${\boldsymbol{\varPhi }}^{(t)}\}$的上界，可得

$ \begin{gathered} \mathfrak{I}\left(\boldsymbol{\varTheta}^{(t+2)} ; \boldsymbol{\varPhi}^{(t+1)}, \beta^{(t+1)}\right) \leqslant \\ \mathfrak{I}\left(\boldsymbol{\varTheta}^{(t+1)} ; \boldsymbol{\varPhi}^{(t+1)}, \beta^{(t+1)}\right)=\mathfrak{I}\left(\boldsymbol{\varTheta}^{(t+1)} ; \boldsymbol{\varPhi}^{(t)}, \beta^{(t)}\right)+ \\ \frac{\beta^{(t+1)}+\beta^{(t)}}{2\left(\beta^{(t)}\right)^2}\left\|\boldsymbol{\varPhi}^{(t+1)}-\boldsymbol{\varPhi}^{(t)}\right\|_{\mathrm{F}}^2 \leqslant\\ \mathfrak{J}\left(\boldsymbol{\varTheta}^{(t+1)} ; \boldsymbol{\varPhi}^{(t)}, \beta^{(t)}\right)+\frac{\beta^{(t+1)}+\beta^{(t)}}{2\left(\beta^{(t)}\right)^2} a^2 \leqslant \\ \Im\left(\boldsymbol{\varTheta}^{(1)} ; \boldsymbol{\varPhi}^{(0)}, \beta^{(0)}\right)+\frac{a^2}{\beta^{(0)}} \sum\limits_{t=0}^{\infty} \frac{\tau+1}{2 \tau^t}<+\infty \end{gathered} $

(24)

所以，得到的函数序列$\{\mathfrak{J}({\boldsymbol{\varTheta}}^{(t+1)}$; ${\boldsymbol{\varPhi }}^{(t)}, {\boldsymbol{ \varPsi }}^{(t)}, β^{(t)})\}$也是有上界的。因此，$\{{\boldsymbol{X}}^{(t)}_{i}\}$和$\{{\boldsymbol{H}}^{(t)}\}$都是有界序列。

然后，利用强凸函数的定义和性质(Boyd和Vandenberghe，2013)可得，存在$a_{i}>0(i=1, 2)$使以下不等式成立。即

$\begin{gathered} \frac{a_1}{2}\left\|\boldsymbol{X}_i^{(t+1)}-\boldsymbol{X}_i^{(t)}\right\|_{\mathrm{F}}^2+\frac{a_2}{2}\left\|\boldsymbol{H}^{(t+1)}-\boldsymbol{H}^{(t)}\right\|_{\mathrm{F}}^2 \leqslant \\ \mathfrak{I}\left(\boldsymbol{\varTheta}^{(t)} ; \boldsymbol{\varPhi}^{(t)}, \beta^{(t)}\right)-\mathfrak{I}\left(\boldsymbol{\varTheta}^{(t+1)} ; \boldsymbol{\varPhi}^{(t+1)}, \beta^{(t+1)}\right)+ \\ \frac{\beta^{(t)}+\beta^{(t+1)}}{2\left(\beta^{(t)}\right)^2}\left\|\boldsymbol{\varPhi}^{(t+1)}-\boldsymbol{\varPhi}^{(t)}\right\|_{\mathrm{F}}^2 \end{gathered} $

(25)

由于$\mathfrak{J}({\boldsymbol{\varTheta}}$; ${\boldsymbol{\varPhi }}, β)≥0$，并且序列$\{{\boldsymbol{\varPhi }}^{(t)}\}$是有上界的，可得

$\begin{array}{r} \sum\limits_{t=0}^{\infty}\left\{\frac{a_1}{2}\left\|\boldsymbol{X}_i^{(t+1)}-\boldsymbol{X}_i^{(t)}\right\|_{\mathrm{F}}^2+\right. \\ \left.\frac{a_2}{2}\left\|\boldsymbol{H}^{(t+1)}-\boldsymbol{H}^{(t)}\right\|_{\mathrm{F}}^2\right\}<+\infty \end{array} $

(26)

因此，序列$\{{\boldsymbol{X}}^{(t)}_{i}\}$和$\{{\boldsymbol{H}}^{(t)}\}$是收敛的，以及

$ \begin{gathered} \lim _{t \rightarrow \infty}\left\|\boldsymbol{H}^{(t+1)}-\nabla \boldsymbol{X}_i^{(t+1)}\right\|_{\mathrm{F}}= \\ \lim _{t \rightarrow \infty} \frac{1}{\beta^{(t)}}\left\|\boldsymbol{\varPhi}^{(t+1)}-\boldsymbol{\varPhi}^{(t)}\right\|_{\mathrm{F}}=0 \end{gathered} $

(27)

为了测试ANL3DTV算法的去噪效果，与TV、TDV、NLTV、RNLTV、NLGRTV和HNHOTV-OGS等算法进行对比。软件平台采用MATLAB R2016b，硬件平台采用Intel Core (TM) i7-7700U CPU处理器，16 GB内存，Windows 10操作系统。测试数据为从TID2013(tampere image database 2013)数据库(Ponomarenko等，2013)和一些常用的测试图像(例如Lena、Peppers、Baboon等)中随机选取的20幅尺寸为341 × 256像素和256 × 256像素的彩色图像，并用两种方式给这些图像添加高斯噪声。第1种噪声情况是给整幅图像添加同一强度的高斯噪声，$σ$为10、30和50；第2种噪声情况是给图像的3个颜色通道添加不一样强度的高斯噪声，[$σ_\text{R}, σ_\text{G}, σ_\text{B}$]为[5, 15, 10]、[40, 50, 30]、[5, 40, 15]和[40, 5, 25]。定量评价指标采用峰值信噪比(peak signal-to-noise ratio，PSNR)和结构相似性(structural similarity，SSIM)。ANL3DTV算法中图像块边长$p$为7，相似块数量$m$为60，参数$τ$为1.5。对比算法按相关文献进行参数设置。

表 1列出了第1种噪声情况下ANL3DTV和对比算法在测试数据集上的平均PSNR和平均SSIM结果。由表 1可知，当$σ$为10、30和50时，ANL3DTV的PSNR分别为32.33 dB、26.92 dB和24.57 dB，SSIM分别为92.99%、81.68%和73.57%。ANL3DTV在3种噪声强度下的平均PSNR和平均SSIM分别比TV高出1.22 dB和3.45%，比TDV高出0.40 dB和0.92%，比NLTV高出0.33 dB和0.95%，比RNLTV高出0.61 dB和1.14%，比NLGRTV高出1.06 dB和2.39%，比HNHOTV-OGS高出0.16 dB和0.12%。

表 2列出了第2种噪声情况下ANL3DTV和对比算法在测试数据集上的平均PSNR和平均SSIM结果。由表 2可知，当[$σ_\text{R}, σ_\text{G}, σ_\text{B}$]为[5, 15, 10]、[40, 50, 30]、[5, 40, 15]和[40, 5, 25]时，ANL3DTV的PSNR结果分别为31.62 dB、24.49 dB、27.47 dB和26.81 dB，SSIM分别为92.88%、73.02%、85.94%和81.00%。ANL3DTV$_\text{w}$为式(10)中的保真项内去掉权重矩阵${\boldsymbol{W}}_{i}$时得到的算法。可以看出，添加权重矩阵${\boldsymbol{W}}_{i}$后，ANL3DTV的平均PSNR和平均SSIM分别比ANL3DTV$_\text{w}$提升了0.89 dB和3.01%。ANL3DTV在4种噪声情况下的平均PSNR和平均SSIM结果分别比TV高出1.32 dB和3.02%，比TDV高出1.76 dB和6.13%，比NLTV高出0.29 dB和0.19%, 比HNHOTV-OGS高出0.55 dB和0.71%。

图 2展示了噪声强度为$σ$=50时ANL3DTV与其他对比算法的去噪视觉效果。从放大的细节图可以看出，TV、TDV、RNLTV和NLGRTV过度平滑了图像的结构，尤其是NLGRTV，几乎完全丢失了百叶窗的结构信息。虽然NLTV和HNHOTV-OGS保留了相对较多的细节信息，但是NLTV引入了大量明显的伪影，HNHOTV-OGS也受到强噪声的影响，误将噪声当做有用的边缘信息保留了下来。相比之下，ANL3DTV仅丢失了很少的细节，而且更有效地滤除了噪声。

图 3和图 4分别展示了噪声强度[$σ_\text{R}, σ_\text{G}, σ_\text{B}$]为[5, 15, 10]和[40, 5, 25]时ANL3DTV与对比算法的去噪视觉效果。可以看出，TV、TDV、HNHOTV-OGS和NLTV丢失了较多的图像细节。例如水面上的波纹、白墙上的文字和屋顶的瓦片等。在图 4(a)中，可以看到黄色墙面上有3条很细的电线，在ANL3DTV的去噪图像中还能隐约看出，而在对比算法的去噪图像中均完全消失。

为了进一步说明去噪算法保持图像纹理的能力，图 5展示了噪声强度$σ$=10和[$σ_\text{R}, σ_\text{G}, σ_\text{B}$]=[5, 40, 15]时，噪声图像和未加噪图像在红、绿和蓝3个颜色通道上的梯度分布。其中，横坐标为$|\nabla x_{i}|=|\nabla _{x}x_{i}|+|\nabla _{y}x_{i}|$，即灰度图像${\boldsymbol{X}}$在像素$i$处的梯度的绝对值。由图 5(a)可知，由于噪声的影响，图像中接近零梯度的像素数目下降较为明显。在$|\nabla x_{i}|>25$处的梯度分布变化不大，此时3个颜色通道的梯度分布非常相似。由图 5中绿色通道对应的梯度分布可以看出，随着绿色通道的噪声强度$σ$从10增大到40，梯度分布曲线逐渐向右偏移，尤其是$σ$=40时，梯度$|\nabla x_{i}|>25$对应的像素数目明显增加。

图 6显示了噪声强度$σ$=10和$σ$=[5, 40, 15]时NLTV、HNHOTV-OGS和ANL3DTV等3种算法获得的去噪图像的梯度分布图。可以看出，随着噪声强度$σ$的增加，NLTV和HNHOTV-OGS过度估计接近零梯度的像素数目，丢失了大量边缘信息。尤其在[$σ_\text{R}, σ_\text{G}, σ_\text{B}$]=[5, 40, 15]时，NLTV和HNHOTV-OGS的梯度分布曲线明显向左偏移。相比之下，ANL3DTV在两种噪声强度下的梯度分布都更接近未加噪图像。实验表明，ANL3DTV在保持图像纹理方面表现更佳，并且对噪声更鲁棒。

ANL3DTV中需要调试的参数为正则化参数$λ$和惩罚因子$β$。图 7显示了噪声强度$σ$= 50时不同参数值对去噪后的PSNR和SSIM的影响。由图 7可知，当$λ$∈(0.07, 0.1]，$β$∈(0.07, 0.1]时，ANL3DTV得到最优的PSNR和SSIM值。对于其他噪声强度，使用类似方法选择得到参数设置。在第1种噪声情况下，当$σ$为10、30和50时，$λ$的取值范围分别为(0.01, 0.07]、(0.04, 0.1]和(0.07, 0.1]，$β$的取值范围分别为(0.5, 1]、(0.07, 0.1]和(0.07, 0.1]。在第2种噪声情况下，当[$σ_\text{R}, σ_\text{G}, σ_\text{B}$]为[5, 15, 10]、[40, 50, 30]、[5, 40, 15]和[40, 5, 25]时，$λ$的取值范围分别为(0.01, 0.07]、(0.001, 0.007]、(0.004, 0.007]和(0.007, 0.01]，$β$的取值范围分别为(0.01, 0.04]、(0.07, 0.1]、(0.04, 0.1]和(0.04, 0.1]。

图 8展示了噪声强度$σ$=50时，ANL3DTV的PSNR随迭代次数$t$增多的变化曲线。可以看出，当$t$＞40时，PSNR趋于稳定。实验结果表明，ANL3DTV具有良好的收敛性。

ANL3DTV和对比算法在341×256像素图像上的平均运行时间如表 3所示。可以看出，TV所需时间最短，RNLTV和ANL3DTV所需时间相对较长。主要是因为ANL3DTV需要为每一个图像块搜索相应的非局部相似图像块，并计算3维全变分，而TV只需处理像素的局部梯度信息。

针对现有许多彩色图像去噪算法未充分利用彩色图像的非局部和局部相似性信息，并且忽略真实噪声在不同图像块和颜色分量间的分布差异，提出了一种自适应非局部3维全变分(ANL3DTV)去噪算法。ANL3DTV利用一个非局部3维全变分正则项捕获彩色图像块内和块间的梯度信息，同时引入一个自适应权重矩阵，可以根据噪声强度的变化来调节算法在不同图像块和不同颜色分量上的去噪强度。实验从定量评价和视觉效果两个方面比较ANL3DTV与6个基于全变分的算法在不同高斯噪声强度下的去噪性能。结果表明，相较于对比算法，ANL3DTV在PSNR和SSIM指标上分别增加了0.16~1.76 dB和0.12%~6.13%。同时，ANL3DTV提升了去噪视觉效果，对噪声变化具有更强的鲁棒性。

由于ANL3DTV在搜索相似性图像块时耗费时间相对较长，因此不利于算法拓展至实时应用。另外，虽然ANL3DTV在去除高斯噪声上取得了令人满意的结果，但对于非高斯噪声的处理存在局限性。在后续工作中，针对彩色图像混合噪声的分布特点，将考虑彩色图像空间的几何结构特征，结合一些较先进的图像处理技术，对相似性测量方式和去噪模型进行优化，进一步提升算法的综合性能。