论文笔记13：Hyperspectral Image Denoising Employing a Spectral–Spatial Adaptive Total Variation Model

引言

论文地址

过去的方法：小波，PCA，张量等。

TV算法具有很好的边缘信息保持能力，本文提出了一种光谱-空间自适应高光谱全变分（SSAHTV）去噪算法，该算法既考虑了不同波段间噪声强度的差异，又考虑了不同像素点之间的空间特性差异。

MAP高光谱图像去噪模型

高光谱图像的噪声退化模型可以写成$f=u+n$，其中$u=\left[u_1,u_2,\dots ,u_j\dots ,u_B\right]$为干净图像，大小为$M\times N\times B$，$f=$ $\left[f_1,f_2,\dots ,f_j,\dots ,f_B\right]$为噪声退化图像，$n=\left[n_1,n_2,\dots ,n_j,\dots ,n_B\right]$为加性噪声。

基于MAP估计理论，高光谱图像的去噪模型可以表示为以下约束最小二乘问题$$\hat{u}=\arg \min \left\{\sum _{j=1}^B{\Vert u_j-f_j\Vert }_2^2+\lambda R(u)\right\}$$

其中${\sum }_{j=1}^B{\Vert u_j-f_j\Vert }_2^2$为数据保真度项，表示观测噪声图像与原始清晰图像之间的保真度，$R(u)$是正则化项，它给出了原始清晰高光谱图像的先验模型。λ是正则化参数，它控制数据保真度和正则化项之间的折衷(tradeoff).

SSAHTV模型

TV Model

对于灰度图像$u$，TV模型定义如下：$$TV(u)=\sum _i\sqrt{{\left({\nabla }_i^{h}u\right)}^2+{\left({\nabla }_i^{v}u\right)}^2}$$

其中${\nabla }_i^h$和${\nabla }_i^v$分别对应于像素$i$处的水平和垂直一阶差分的线性算子，${\nabla }_i^{h}u=u_i-u_{r(i)}$和${\nabla }_i^{v}u=u_i-u_{b(i)}$，其中$r(i)$和$b(i)$表示像素右侧和下方最近邻。

Spectral Adaptive Hyperspectral TV Model

简单的逐波段高光谱TV模型定义如下：$$HTV(u{)}_1=\sum _{j=1}^{B}TV\left(u_j\right)$$

于是去噪模型可以写成$$\hat{u}=\arg \min \left\{\sum _{j=1}^B{\Vert u_j-f_j\Vert }_2^2+\lambda \sum _{j=1}^{B}TV\left(u_j\right)\right\}$$

对于上式，Euler–Lagrange方程[Ref: Deblurring of color images corrupted by impulsive noise]如下：$$\left(u_j-f_j\right)-\lambda \nabla \cdot \frac{\nabla u_j}{\left|\nabla u_j\right|}=0(1)$$

这意味着每个波段都被单波段TV模型单独去噪，在上式中，如果我们对所有波段使用相同的正则化参数λ，这意味着每个波段的正则化强度相等，将导致以下结果：1）如果使用较大的λ很好地对高噪声强度波段去噪，则低噪声强度波段将被过度平滑；2）相反，当用较小的λ对低强度波段进行良好的去噪处理时，高噪声强度波段中的噪声不会得到很好的抑制。虽然可以在不同的波段使用不同的λ来克服这一缺点，但是手动调整不同波段的λ非常耗时。

我们扩展了[Fast dual minimization of the vectorial total variation norm and applications to color image processing]中提出的CTV模型，定义了高光谱TV模型：$$\begin{array}{rl}HTV(u{)}_2& =\sum _{i=1}^{MN}\sqrt{\sum _{j=1}^B{\left({\nabla }_{ij}u\right)}^2}\\ {\left({\nabla }_{ij}u\right)}^2& ={\left({\nabla }_{ij}^{h}u\right)}^2+{\left({\nabla }_{ij}^{v}u\right)}^2\end{array}$$

于是去噪模型可以写成$$\hat{u}=\arg \min \left\{\sum _{j=1}^B{\Vert u_j-f_j\Vert }_2^2+\lambda \sum _{i=1}^{MN}\sqrt{\sum _{j=1}^B{\left({\nabla }_{ij}u\right)}^2}\right\}$$

如果取$u_j$的导数，上式的Euler–Lagrange方程可以写成$$\left(u_j-f_j\right)-\lambda \nabla \cdot \frac{\nabla u_j}{\sqrt{{\sum }_{j=1}^B{\left(\nabla u_j\right)}^2}}=0$$

$$\left(u_j-f_j\right)-\lambda \nabla \cdot \frac{\left|\nabla u_j\right|}{\sqrt{{\sum }_{j=1}^B{\left|\nabla u_j\right|}^2}}\cdot \frac{\nabla u_j}{\left|\nabla u_j\right|}=0$$

与(1)式对比，我们可以看到在上式中增加了调整参数$\frac{|\nabla u_j|}{\sqrt{{\sum }_{j=1}^B{|\nabla u_j|}^2}}$，以自动调整各波段的去噪强度。对于高噪声强度带，如$\frac{|\nabla u_j|}{\sqrt{{\sum }_{j=1}^B{|\nabla u_j|}^2}}$值较大，这些波段的去噪强度将很强。相反，对于低强度噪声的波段，如$\frac{|\nabla u_j|}{\sqrt{{\sum }_{j=1}^B{|\nabla u_j|}^2}}$值较小，对其进行弱去噪处理。

SSAHTV Model

在分析了高光谱TV模型的光谱自适应特性后，另一个重要的问题是在去噪过程中如何实现空间自适应，即如何根据空间结构分布，在同一波段内不同像素位置调整去噪强度。

对于高光谱图像$u$，我们首先使用以下公式计算每个波段的梯度信息：
$$\nabla u_j=\sqrt{{\left({\nabla }^h{u}_j\right)}^2+{\left({\nabla }^v{u}_j\right)}^2}$$

其中${\nabla }^h{u}_j$和${\nabla }^v{u}_j$分别代表$u_j$的水平和垂直方向的一阶梯度，$({\nabla }^h{u}_j{)}^2$代表${\nabla }^h{u}_j$每个元素的平方。然后，将每个波段的梯度信息相加$$G=\sqrt{\sum _{j=1}^B{\left(\nabla u_j\right)}^2}$$

设$G_i$为向量$G$的第$i$个元素$(i=1,\cdots ,MN)$，控制波段间(interband)去噪强度的权重参数$W_i$定义如下：$${\tau }_i=\frac{1}{1+\mu G_i}$$

$$W_i=\frac{{\tau }_i}{\bar{\tau }}\bar{\tau }=\frac{{\sum }_{i=1}^{MN}{\tau }_i}{MN}$$

其中$\mu$是常数参数，$\tau \in [0,1]$，$\bar{\tau }$为${\tau }_i$的均值。

为了使去噪过程在空间上具有自适应性，在$HTV(u{)}_2$中加入参数$W_i$(高光谱TV模型中第$i$个像素的空间权重)，将SSAHTV模型定义为
$$SSAHTV(u)=\sum _{i=1}^{MN}{W}_i\sqrt{\sum _{j=1}^B{\left({\nabla }_{ij}u\right)}^2}$$

对于平滑区域中的像素，梯度信息$G_i$的值将很小，并且空间权重$W_i$将具有较高的值。因此，我们将对这些像素点使用强大的去噪强度，从而更好地抑制平滑区域中的噪声。相反，对于边缘和纹理区域中的像素，$G_i$的值将较大，而空间参数$W_i$的值将较小。因此，将使用弱去噪强度对其进行处理，并保留边缘和细节信息。

本文使用的最终MAP去噪模型可以写成：
$$\hat{u}=\arg \min \left\{\sum _{j=1}^B{\Vert u_j-f_j\Vert }_2^2+\lambda \sum _{i=1}^{MN}{W}_i\sqrt{\sum _{j=1}^B{\left({\nabla }_{ij}u\right)}^2}\right\}(2)$$

分裂Bregman优化

[37] L. Bregman, “The relaxation method of finding the common points of convex sets and its application to the solution of problems in convex optimization,” USSR Comp. Math. Math., vol. 7, pp. 200–217, 1967.
[38] T. Goldstein and S. Osher, “The split bregman method for L1 regularized problems,” SIAM J. Imag. Sci., vol. 2, no. 2, pp. 323–343, May 2009.

首先，在优化过程中引入辅助变量$d$，然后按以下方式加至(2)中
$$\begin{array}{r}\hat{u}=\arg \min \left\{\Vert u-f{\Vert }_2^2+\lambda {\sum }_{i=1}^{MN}{W}_i\sqrt{{\sum }_{j=1}^B{\left(d_{ij}\right)}^2}\right\}\\ \text{ subject to }:d=\nabla u\end{array}$$

上式中的约束问题可用Bregman迭代法转化为无约束问题：
$$\hat{u}=\arg \min \left\{\Vert u-f{\Vert }_2^2+\lambda \sum _{i=1}^{MN}{W}_i\sqrt{\sum _{j=1}^B{\left(d_{ij}\right)}^2}+\beta \Vert d-\nabla u-b{\Vert }_2^2\right\}$$

在上式中，变量$b$也是一个辅助变量来加速迭代。上式的最小化可以通过以下两个子问题交替进行
subproblem $u:$
$$\hat{u}=\arg \min \left\{\Vert u-f{\Vert }_2^2+\beta \Vert d-\nabla u-b{\Vert }_2^2\right\}$$

subproblem $d:$
$$\hat{d}=\arg \min \left\{\lambda \sum _{i=1}^{MN}{W}_i\sqrt{\sum _{j=1}^B{\left(d_{ij}\right)}^2}+\beta \Vert d-\nabla u-b{\Vert }_2^2\right\}$$

在上述方程中，为了求解$u$子问题，必须求解以下方程$$(I-\beta \Delta )u^{k+1}=f+\beta {\nabla }^T(d-b)$$

$\Delta$为Laplace算子，对于上式中的线性函数，由于系统是严格对角的，最有效的方法之一是使用Gauss–Seidel迭代算法。

$d$子问题方程可以通过收缩算子求解：$$\begin{array}{rlr}d& =\mathrm{shrink}\left(\sqrt{\sum _{j=1}^B{\left({\nabla }_{ij}{u}^{k+1}+b^k\right)}^2},\frac{\lambda W}{\beta }\right)& \\ & =\max \left(0,\sqrt{\sum _{j=1}^B{\left({\nabla }_{ij}{u}^{k+1}+b^k\right)}^2}-\frac{\lambda W}{\beta }\right)\times \frac{{\nabla }_{ij}{u}^{k+1}+b}{\sqrt{{\sum }_{j=1}^B{\left({\nabla }_{ij}{u}^{k+1}+b\right)}^2}}(23)& \end{array}$$

收缩算子为软阈值法( s o f t ( u , a ) = s i g n ( u ) × m a x { ∣ u ∣ − a , 0 } soft(u,a)=sign(u)×max\{|u|-a,0\} soft(u,a)=sign(u)×max{∣u∣−a,0}). 最后，对于参数$b$，在每次迭代中需要按照以下方式进行更新$$b^{k+1}=b^k+\left(\nabla u^{k+1}-d^{k+1}\right)$$

Split Bregman方法的优点是将(2)中比较困难的优化问题分解为上述两个子问题，非常容易优化。利用Split Bregman优化算法，(2)中SSAHTV去噪模型的优化过程总结如下：

实验结果与讨论

模拟数据的实验：
数据：Washington DC Mall
情况一：对于不同波段，噪声强度是相等的：所有波段均加入相同的零均值高斯噪声分布，方差为0.02。平均信噪比(SNR)为22.39。
情况二：不同波段的噪声强度不同，在各波段加入不同方差的零均值高斯噪声，方差值在0 ~ 0.05之间随机选取。平均信噪比为22.28。
对比方法：局部自适应维纳滤波器(MATLAB中的wiener2函数)、小波硬阈值去噪、逐波段TV模型和高光谱TV模型。

评估指标：
其中${\mu }_u$表示原始干净图像的平均灰度值，${\sigma }_u$表示原始干净图像的方差，${\sigma }_{u\hat{u}}$表示协方差，$C_1$和$C_2$为常数，防止当${\mu }_u^2+{\mu }_{\hat{u}}^2$或${\sigma }_u^2+{\sigma }_{\hat{u}}^2$非常接近0时出现不稳定结果。

从图4、图5和图7、图8(见原论文)可以看出，在SSAHTV去噪结果中，不仅能更彻底地抑制噪声，而且很好地保留了边缘信息和细节信息。然而，局部自适应维纳滤波的去噪结果被过度平滑，大量的细节信息丢失。在小波硬阈值去噪的结果中，图像中产生了一些“伪影”，边缘信息也没有得到很好的保存。使用逐波段TV模型和高光谱TV模型去噪的结果也不如使用SSAHTV模型去噪的图像清晰，平滑区域仍然存在一些噪声。

在图6和图9(见原论文)中可以清楚地看到，使用SSAHTV模型的大部分波段的PSNR和SSIM值都高于其他四种方法的PSNR和SSIM值。在表I和表II的整体定量评估中，本文提出的SSAHTV模型在所有算法中MPSNR和MSSIM值最高。

真实数据的实验：数据：Indian Pines

3、110、204波段的去噪结果分别如图11-13所示。图14为三个波段的组合结果。为了进一步验证所提出的去噪方法的有效性，给出了去噪前后高光谱图像的分类结果，以供比较。我们首先从总样本中随机选取少量样本作为训练样本，然后使用支持向量机(SVM)分类方法对图像进行监督分类。最后，用剩余样本测试分类准确率。

SVM对整个高光谱图像的分类结果如图15(a) - (e)所示。在图16(a) - (e)中，我们也展示了PCA后SVM的分类结果。采用OA和kappa系数的分类精度评价结果见表IV和表V。

第二组实验数据：HYDICE urban image(包含相当多的高强度条纹和混合噪声波段)

结果见原论文。

参数$\mu$的鲁棒性讨论见原论文。

虽然所提出的SSAHTV模型在高光谱图像去噪问题上表现良好，但在某些方面还可以进一步改进。我们未来的工作将集中在模型构建过程中考虑光谱维度的梯度，以开发一个真实的三维TV模型。此外，可以使用三维分割或聚类结果从区域角度来约束去噪过程，而不是本文中的像素角度。