Hyperspectral Mixed Noise Image Restoration Based on Non-Convex Low-Rank Tensor Decomposition and Group Sparse Total Variation
-
摘要:
高光谱图像（Hyperspectral Image，HSI）在采集的过程中会被大量混合噪声污染，会影响遥感图像后续应用的性能，因此从混合噪声中恢复干净的HSI成为了重要的预处理过程。在本文中，提出了一种基于非凸低秩张量分解和群稀疏总变分正则化的高光谱混合噪声图像恢复模型；一方面，采用对数张量核范数来逼近HSI的低秩特性，可以利用高光谱数据固有的张量结构，同时减少对较大奇异值的收缩以保留图像更多细节特征；另一方面，采用群稀疏总变分正则化来增强HSI的空间稀疏性和相邻光谱间的相关性。并采用ADMM（Alternating Direction Multiplier Method）算法求解，实验证明该算法易于收敛。在模拟和真实的高光谱图像实验中，与其他方法相比，该方法在去除HSI混合噪声方面具有更好的性能。
Abstract:Hyperspectral images (HSIs) are polluted by a large amount of mixed noise during the acquisition process, which affects the performance of subsequent applications of remote sensing images. Therefore, restoring clean HSI from the mixed noise is an important preprocessing step. In this study, a hyperspectral mixed noise image restoration model based on nonconvex low-rank tensor decomposition and group-sparse total variational regularization is proposed. On the one hand, by using logarithmic tensor nuclear norm to approximate the low-rank characteristics of the HSI, the inherent tensor structure of hyperspectral data can be utilized, and the shrinkage of larger singular values can be reduced to preserve more detailed features of the image. On the other hand, the group sparse total variational regularization can be used to enhance the spatial sparsity of the HSI and correlation between adjacent spectra. ADMM algorithm is used to solve the problem, and an experiment shows that the algorithm converges easily. In simulated and real hyperspectral image experiments, this method performs better in removing HSI mixed noise when compared to other methods.
-
0. 引言
利用同一场景的光谱仪来获取不同光谱下的图像被称为高光谱成像。它包含了比普通图像更丰富的信息，大大提高了地面覆盖识别和特征分析的准确性和可靠性。因此，它有多种应用，包括环境研究、农业、军事、地理等[1]。同时由于高光谱图像提供了丰富的光谱信息，在分类[2]、超分辨率[3]和解混[4]等方面也发挥着重要作用。然而由于高光谱传感器在采集数据时容易受到外界环境的影响，数据不可避免地会受到高斯噪声、椒盐噪声、条带和死线等各种混合噪声污染，从而降低了HSI的质量，限制了其后续的应用[5]。因此，从已被噪声污染的HSI中恢复干净的HSI是很重要且有意义的研究。
在过去的几十年里，许多不同的HSI恢复方法被提出来提高其质量。其中较简单的技术是使用传统的二维灰度图像和一维信号方法来实现HSI逐带[6]或逐像素[7]去除噪声。然而，这些去噪方法没有考虑到HSI光谱域的低秩先验特性。利用这一特性，Candes等人[8]提出了一个最小化核范数的鲁棒主成分分析（Robust Principal Component Analysis，RPCA）模型，以有效地获得清晰的图像。Gu等人[9]提出了一种加权核范数极小化（Weighted Kernel Norm Minimization，WNNM）模型，通过分配不同的权重到奇异值。He[10]等人建立了针对不同波段不同噪声强度的低秩矩阵逼近方法，有效地提高了图像恢复的信噪比。然而，张量矩阵化破坏了张量的高阶结构，这对去噪后的局部细节和重要信息产生了影响。
因此，许多算法引入基于低秩张量的方法来描述HSI的低秩特性。包括Tucker分解[11]和张量奇异值分解（tensor Singular Value Decomposition，t-SVD）[12]。然而，如文献[13]所示，低秩正则化并不足以描述HSI的空间先验信息。总变分正则化是图像处理中保持局部空间分段平滑性常用的工具，因此，许多人提出了基于低秩矩阵/张量分解框架和各种类型的总变分正则化模型，以同时探索HSI的空间和光谱先验。例如，He[13]等人将HTV正则化引入低秩矩阵分解框架，以提高恢复结果。在文献[14]中，将低Tucker秩模型和SSTV（Spatial Spectrum Total Variation）正则化项组合分别利用全局空间谱和光谱相关性，增强空间信息；如范等[15]将SSTV正则化加到低阶张量分解框架中（Spatial Spectrum Total Variation-Low Rank Tensor Decomposition，SSTV-LRTF），使用张量核范数（Tensor Kernel Norm，TNN）近似HSI的低秩属性并同时利用空间和光谱域之间的分段光滑性来完成HSI去噪。Chen[16]等人提出加权群稀疏正则化低秩张量分解方法，将空间差分图像的群稀疏性正则化和Tucker分解结合起来，有效探索了不同光谱波段空间差图像的共享群稀疏模式。近期，许多基于深度学习的算法被提出用于HSI图像去噪，例如3DADNet[17]，SSCAN[18]，然而这些方法都只能去除单一的高斯噪声，而高光谱图像往往都会被混合噪声污染，这会使其在实际运用中达不到理想的效果。总之，低秩张量框架和总变分正则化的结合在高光谱混合噪声图像恢复时可以达到最优的结果。
然而，SSTV-LRTF的方法中采用张量核范数作为其低秩的凸松弛，它测量非零奇异值的L1范数，这不是张量低秩的一个很好的近似[19]。其次，它平等地对待每个奇异值，因此可能无法很好地保存主要信息。这是因为较大的奇异值通常对应于主要信息，如轮廓、锐利的边缘和光滑的区域，因此收缩应小于较小的奇异值[20-24]。因此，为解决上述问题，我们采用对数张量核范数作为其低秩的非凸松弛。
同时由于SSTV[14]中在空间维度做差分时采用简单的L1范数，其并不能充分表征张量梯度域的空间稀疏结构，最终导致去噪后的图像过于平滑。Chen等[16]设计了一个加权L2, 1正则化器来探索张量空间维度的群稀疏性，取得更好的结果。为弥补SSTV正则项不足，提出了新的群稀疏总变分正则项，来保持HSI局部空间光谱相关性的同时探索其空间维度的群稀疏性。
为了去除高光谱数据的混合噪声，将对数张量核范数和群稀疏总变分正则项结合起来。在模拟和真实HSI数据实验中与其他去噪方法相比，取得了更好的结果。
主要贡献有以下3点：
1）在非凸低秩张量分解模型（Non-convex low-rank tensor factorization，NCLRTF）中，采用对数张量核范数来作为HSI的低秩非凸松弛，可以更好地近似其低秩并保留图像更多的特征信息。
2）提出了新的群稀疏总变分正则项，其考虑HSI空间光谱维度的局部相关性的同时可以探索HSI共享的群稀疏模式，因此能够保留图像更多的细节特征。并将其与非凸低秩张量分解的模型结合起来，提高了图像恢复的结果。
3）为了求解该模型，我们采用了ADMM（交替方向乘子法）算法[25]来获得全局最优解，实验结果表明，该算法易于收敛。
1. 相关工作
1.1 HSI退化模型
基于观测到的HSI总会受到混合噪声污染，假设被混合噪声污染的HSI数据用Y∈RM×N×p表示，其中M×N为空间域，p为谱带数。它们可以被建模为：
$$ Y=X+S+N $$ (1) 式中：X∈RM×N×p为干净的高光谱数据；S∈RM×N×p和N∈RM×N×p分别为稀疏噪声和高斯噪声。HSI恢复的重点是如何从观察到的噪声图像中恢复干净图像。
1.2 HSI低秩张量分解恢复模型
在文献[15]中，将高光谱数据每个波段的二维图像重塑为大小为M×1×N的横向切片，然后将这些横向切片排列，得到大小为M×p×N的三维张量$ \hat X $。
可由以下操作实现：Map(RM×N×p)→(RM×p×N)，其逆算子为：iMap(RM×p×N)→(RM×N×p)。对于三向张量$ \hat X \in {R^{M \times p \times N}} $，可以找到两个张量A∈RM×k×N和B∈Rk×p×N满足：
$$ \hat X = A*B = \sum\limits_{j = 1}^k {A\left( {:,j,:} \right)*B\left( {j,:,:} \right)} $$ (2) 因此使用低秩张量分解重构HSI退化模型：
$$ Y={\rm{iMap}}(A*B)+S+N $$ (3) 对于三向张量$ \hat X \in {R^{M*p*N}} $，假设k是它的管状张量秩的上界[26]，可得：
$$ {\left\| {\hat X} \right\|_{{\text{TNN}}}} = \mathop {\inf }\limits_{A,B} \left\{ {\frac{1}{2}\left\| A \right\|_{\text{F}}^2 + \left\| B \right\|_{\text{F}}^2:\hat X = A*B} \right\} $$ (4) 式中：张量A∈RM×k×N和B∈Rk×p×N。
因此可以得到以下低秩张量分解的去噪模型[15]：
$$ \begin{array}{l} \mathop {\min }\limits_{X,S} \lambda {\left\| S \right\|_1} + {\left\| {{\text{Map}}\left( X \right)} \right\|_{{\text{TNN}}}}\;{\text{s}}{\text{.t}}. \hfill \\ \left\| {Y - X - S} \right\|_{\text{F}}^2 \leqslant \varepsilon \quad {\text{rank}}\left( {{\text{Map}}\left( X \right)} \right) \leqslant k \hfill \\ \end{array} $$ (5) 该模型等同于受秩约束的TRPCA[27]的拉格朗日形式。
2. 提出的模型
式(5)中采用张量核范数作为其低秩的凸松弛，它测量非零奇异值的L1范数，平等地对待每个奇异值，因此可能无法很好地保存主要信息。这是因为较大的奇异值通常对应于主要信息，如轮廓、边缘和光滑的区域，因此较大的奇异值收缩应小于较小的奇异值[28]。为此采用一个对数张量核范数（Logarithmic tensor kernel norm，LogTNN）来作为低秩的非凸松弛，其可以对较大的奇异值减少收缩以保持主要信息，对较小的奇异值增加收缩以抑制噪声[19]。
张量X∈RM×N×p的基于对数的张量核范数定义为[19]：
$$ {\text{LogTNN}}\left( {X,\varepsilon } \right) = \sum\limits_{i{\text{ = 1}}}^p {{\text{LogMNN}}\left( {{{\bar X}^{\left( i \right)}},\varepsilon } \right)} $$ (6) 式中：$ {\bar X^{\left( i \right)}} $表示$ \bar X = {\text{fft}}\left( {X,\left[ {\ } \right],k} \right) $的第i个切片，其中：
$$ {\text{LogMNN}}\left( {X,\varepsilon } \right) = \sum\limits_{i = 1}^m {\log \left( {{\sigma _i}\left( \boldsymbol X \right),\varepsilon } \right)} $$ (7) 这里σi(X)是X的第i个奇异值，ε＞0是一个常数。因此基于LogTNN的HSI去噪模型表述为：
$$ \begin{array}{l} \mathop {\min }\limits_{X,S} \lambda {\left\| S \right\|_1} + {\text{LogTNN}}\left( {{\text{Map}}\left( X \right),\varepsilon } \right)\quad {\text{s}}{\text{.t}}. \hfill \\ \left\| {Y - X - S} \right\|_{\text{F}}^2 \leqslant \varepsilon \quad {\text{rank}}\left( {{\text{Map}}\left( X \right)} \right) \leqslant k \hfill \\ \end{array} $$ (8) 该模型可以充分利用HSI的全局低秩特性，并减小对较大奇异值的收缩来保留图像更多的细节特征。但该模型中缺乏HSI空间维度的先验信息，而总变分正则化是保持图像局部空间平滑的有效工具。
2.1 SSTV正则化
如文献[29-30]中所述，每个HSI波段都可以看作是一个灰度图像，因此它在空间维数上具有局部分段平滑性。空间光谱总变分（SSTV）可以充分保持HSI空间光谱的局部平滑性，其定义为：
$$ {\left\| X \right\|_{{\text{SSTV}}}} = {\left\| {{D_x}X} \right\|_1} + {\left\| {{D_y}X} \right\|_1} + {\left\| {{D_z}X} \right\|_1} $$ (9) 式中：Dx表示水平方向上的差分算子；Dy是垂直方向上的差分算子；Dz是在每个像素的光谱特征上的一维有限差分算子。SSTV的模型为了充分利用空间谱差图像的稀疏特性，一般使用凸L1范数来描述其稀疏先验。虽然L1范数是促进各频带分段光滑结构的有效约束条件，但是其只描述了非零元素的数量，而忽略了非零元素的局部群结构。因此，Chen等人提出了群稀疏正则化[16]，其可以表示为：
$$\begin{array}{l} {\left\| {\boldsymbol W \odot DX} \right\|_{2,1}} = \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {W\left( {i,j} \right)} } {\left\| {{D_x}X\left( {i,j,;} \right)} \right\|_2} + \hfill \\ \quad \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {W\left( {i,j} \right)} } {\left\| {{D_y}X\left( {i,j,;} \right)} \right\|_2} \hfill \\ \end{array} $$ (10) 式中：W是权重矩阵。采用L2, 1范数的群稀疏性正则化表示在HSI的每个波段边界区域的梯度值明显较大，而平滑区域的梯度值较小。与普通的TV正则项相比，提高了HSI在空间维度的稀疏特性。
为保持HSI空间光谱的局部平滑性的同时可以探索HSI空间维度的群稀疏特性，新的群稀疏TV正则化可以表示为：
$$ ||X|{|_{{\text{GSTV}}}} = ||\boldsymbol W \odot DX|{|_{2,1}} + ||{D_z}X|{|_1} $$ (11) 该GSTV（Group Sparse Total Variation）正则化保证HSI空间维度的群稀疏特性的同时可以增强相邻光谱间的相关性，来保持HSI空间光谱的局部平滑性，以保留图像更多的细节特征，提高了去噪的性能。
2.2 NCLRGSTV模型
基于对数核范数去噪模型是利用高光谱图像的全局低秩特性，来抑制噪声，而GSTV正则项是用来探索HSI空间维度的群稀疏性和增强相邻光谱间的相关性，因此可以结合二者的优势，来恢复含有混合噪声的高光谱图像。为此提出NCLRGSTV模型（Non-Convex Low-Rank Tensor Factorization Group Sparse Total Variation）：
$$ \begin{array}{l} \mathop {\min }\limits_{X,S} {\text{LogTNN}}\left( {{\text{Map}}\left( X \right),\varepsilon } \right) + {\lambda _1}{\left\| S \right\|_1} + {\lambda _2}{\left\| X \right\|_{{\text{GSIV}}}} \hfill \\ {\text{s}}{\text{.t}}{\text{.}}\quad \left\| {Y - X - S} \right\|_{\text{F}}^2 \leqslant \varepsilon \quad {\text{rank}}\left( {{\text{Map}}\left( X \right)} \right) \leqslant k \hfill \\ \end{array} $$ (12) 式中：λ2是用来控制对数张量核范数和GSTV之间权衡的参数。
2.3 优化NCLRGSTV模型
通过ADMM算法引入辅助变量F，式(12)可以重写为：
$$ \begin{array}{l} \mathop {\min }\limits_{X,F,S} {\text{LogTNN}}\left( {\hat X,\varepsilon } \right) + {\lambda _1}{\left\| S \right\|_1} + {\lambda _2}{\left\| F \right\|_{{\text{GSTV}}}} \hfill \\ {\text{s}}{\text{.t}}{\text{.}}\quad \left\| {Y - X - S} \right\|_{\text{F}}^2 \leqslant \varepsilon \quad {\text{rank}}\left( {\hat X} \right) \leqslant k\quad F = X \hfill \\ \end{array} $$ (13) 式中：$ \hat X = {\text{Map}}\left( X \right) $。在ADMM框架内X, F, S交替更新为：
$$ \begin{array}{*{20}{c}} {{X^{k + 1}} = \mathop {\arg \min }\limits_{rank\left( {\hat X} \right) \leqslant k} L\left( {X,{F^t},{S^t},P_1^t,P_2^t} \right)} \\ {{F^{k + 1}} = \mathop {\arg \min }\limits_F L\left( {{X^{t + 1}},F,{S^t},P_1^t,P_2^t} \right)} \\ {{S^{k + 1}} = \mathop {\arg \min }\limits_F L\left( {{X^{t + 1}},{F^{t + 1}},S,P_1^t,P_2^t} \right)} \end{array} $$ (14) 1）更新Xt+1，Xt+1的子问题为：
$$ \begin{array}{l} {X^{t + 1}} = \mathop {\arg \min }\limits_{\hat X} \tau {\text{LogTNN}}\left( {\hat X,\varepsilon } \right) + \left\| {\hat X - \hat Q} \right\|_F^2 \hfill \\ {\text{s}}{\text{.t}}{\text{.}}\quad {\text{rank}}\left( {\hat X} \right) \leqslant k \hfill \\ \end{array} $$ (15) 式中：$ \hat Q = {\text{Map}}\left( Q \right),Q = \frac{1}{2}\left( {Y + X - S + \left( {P_1^t + P_2^t} \right)} \right) $，上式求解可以参考文献[19, 31]。
$$ {X^{t + 1}} = D_{It}^{\tau ,\varepsilon }\left( {\hat Q} \right) = \boldsymbol U*\boldsymbol S_{It}^{\tau ,\varepsilon }*{\boldsymbol V^{\text{T}}} $$ (16) 式中：$ \hat Q = \boldsymbol U*\boldsymbol S*{\boldsymbol V^T} $，$ S_{It}^{\tau ,\varepsilon } = {\text{ifft}}\left( {\bar S_{It}^{\tau ,\varepsilon },\left[ {\ } \right],3} \right) $。
2）更新Ft+1，Ft+1的子问题为：
$$ {F^{t + 1}} = \mathop {{\text{arg}}\min }\limits_{_F} {\lambda _2}||F|{|_{{\text{GSTV}}}} + \frac{\beta }{2}||F + \frac{{P_1^t}}{\beta } - {X^{t + 1}}||_F^2 $$ (17) 上式可由ADMM算法求解。F的子问题是一个最小二乘问题，可以采用高效快速傅里叶变换（Fast Fourier transform，FFT）求解：
$$ F = {\text{ifftn}}\left( {\frac{\boldsymbol G}{{\beta 1 + \beta {{\left| {{\text{fftn}}\left( D \right)} \right|}^2}}}} \right) $$ (18) 式中：$ \boldsymbol G = {\text{fftn}}\left( {X + \left( {\frac{{{P_1}}}{\beta }} \right)} \right) + {\boldsymbol D^{\text{T}}}\left( {R + \frac{{{P_3}}}{\beta },C + \frac{{{P_4}}}{\beta }} \right) $，fftn和ifftn是快速的三维傅里叶变换和它的反变换。
3）更新St+1，St+1的子问题为：
$$ \mathop {\min }\limits_S {\lambda _1}{\left\| S \right\|_1} + \frac{\beta }{2}\left\| {Y - {X^{t + 1}} - S + \frac{{P_1^t}}{\beta }} \right\|_{\text{F}}^2 $$ (19) 通过应用软阈值收缩算子，St+1子问题可以用以下方法精确求解：
$$ {S^{t + 1}} = {\text{shrink}}\left( {Y - {X^{t + 1}} + \frac{{P_1^t}}{\beta },\frac{{{\lambda _1}}}{\beta }} \right) $$ (20) 4）更新拉格朗日乘子：
$$ \begin{array}{*{20}{l}} {{P_1} = P_1^t + \beta \left( {{X^{t + 1}} - {F^{t + 1}}} \right)} \\ {{P_2} = P_2^t + \beta \left( {Y - {X^{t + 1}} - {S^{t + 1}}} \right)} \\ {{P_3} = P_3^t + \beta \left( {R - DF} \right)} \\ {{P_4} = P_4^t + \beta \left( {C - {D_z}F} \right)} \end{array} $$ (21) 总结步骤(1)~(4)的过程，可以得到所提HSI恢复模型的NCLRGSTV最优解，如算法1所示。
算法1：NCLRGSTV求解输入的优化过程：
输入：噪声图像Y，参数λ1，λ2，k，ε，β，t=0
1:初始化：X=F=S=P1=P2=P3=P4=0
2:当不收敛时：
3:通过(16)更新X
4:通过(18)更新F
5:通过(20)更新S
6:通过(21)更新得到P1, P2, P3, P4
7: t=t+1
8:检查收敛条件：$ {\left\| {{X^{k + 1}} - {X^k}} \right\|_{\text{F}}}/{\left\| {{X^k}} \right\|_{\text{F}}} < {10^{ - 3}} $
如果收敛，输出恢复HIS。
2.4 计算复杂度分析
设输入的HSI大小为M×N×p，式(13)中的复杂优化问题被分为几个子问题。更新X子问题的计算复杂度为O(MNp(log(MNp)+p))；更新F子问题采用fft进行优化，需要O(MNp(log(MNp)))的计算；更新R, S, C是软阈值收缩操作，其计算复杂度为O(MNp)；因此整个过程的计算复杂度总计为O(MNp(log(MNp)+p))。
3. 实验结果与讨论
为了验证优化后的NCLRGSTV模型在去除混合噪声中的性能，将NCLRGSTV模型分别应用于模拟数据实验和真实数据实验。此外，为了证明NCLRGSTV算法的去噪有效性，将其与相同条件下的LRMR[10]、LRTV[13]、LRTDTV[14]，LRTDGS[16]、FRCTR-PnP[32]5种高效去噪方法进行了比较，在进行实验之前，所使用的HSI的每个波段的灰度值归一化处理。
3.1 模拟数据实验
模拟实验中，采用由反射光学系统成像光谱仪（rose-03）收集的Pavia城市中心数据集，其原始规模为1096×1096×102。由于Pavia城市中心数据集的一些光谱波段被噪声污染严重，不能作为去噪的参考[30]。因此，去掉该数据的前几个波段，选取子图像的大小为200×200×80。由于模拟实验给出了地面真实度的HSI，我们采用了5个定量图像指标进行比较，包括各个波段平均峰值信噪比（Mean Peak Signal Noise Ratio，MPSNR）、平均结构相似度（Mean Structural Similarity，MSSIM）、平均特征相似度（Mean Feature Similarity，MFSIM）、平均光谱角映射（Mean Spectral Angle Mapping，MSAM）和相对全局合成维度误差（Erreur Relative Globale Adimensionnelle de Synthese，ERGAS）。MPSNR、MSSIM和MFSIM越大，MSAM和ERGAS越小，表示算法的去噪性能越好。
Case 1：对HSI每个波段加入不同强度的高斯噪声，零均值高斯噪声的方差在[0, 0.2]间选择，同时每个波段加入不同强度的椒盐噪声，其百分比在[0, 0.2]间选择。
Case 2：在Case1的基础上，在Pavia City Center数据集的第50~第60波段添加条带噪声。条带噪声的数目在10~30之间随机变化。
Case 3：在Case2基础上，在Pavia City Center数据集的54~64增加死线。上述区域的死线数量在[3, 10]之间随机变化，死线宽度在[1, 3]之间随机变化。
1）定量分析
从表 1中可以看到随着噪声种类的增加，各种去噪算法的去噪性能随之下降，但是NCLRGSTV算法的各项指标除MSAM外都为最优，这充分说明了对数张量核范数和群稀疏总变分的结合能够充分抑制图像中的高斯噪声、椒盐噪声、条带和死线等稀疏噪声。同时可以发现，FRCTR-PnP方法在去除每个波段的噪声强度不同时，效果并不理想且很耗时。虽然LRTDGS中群稀疏正则化探索了空间的群稀疏性，提升了去噪性能，但其忽略了局部空间光谱的相关性，会使MSAM值增大，而提出群稀疏总变分考虑到了这一点，可以看到，NCLRGSTV算法的MSAM值与LRTDGS算法相比下降0.01。与此同时，采用对数张量核范数来表示低秩的非凸松弛也表明了其有效性，与采用核范数或Tucker分解的算法相比，其中MPSNR值提升了1 dB左右，MSSIM和MFSIM提升了0.01左右，ERGAS下降了10左右，这也说明了该算法去除混合噪声的优势。NCLRGSTV运行时间虽然不是最快的，但耗时并不长。
表 1 Pavia city center数据集的不同去噪方法的定量评价结果Table 1. Quantitative evaluation results of different denoising methods in Pavia city center data setsCase Indexes Noise LRMR LRTV LRTDTV LRTDGS FRCTR-PnP NCLRGSTV Case 1 MPSNR 14.144 33.336 34.356 34.743 35.380 34.557 36.369 MSSIM 0.2143 0.9341 0.9444 0.9457 0.9506 0.9370 0.9637 MFSIM 0.5985 0.9590 0.9626 0.9646 0.9647 0.9630 0.9761 MSAM 0.6676 0.0833 0.0545 0.0495 0.0637 0.1331 0.0514 ERGAS 707.54 74.698 65.280 70.351 61.441 109.32 51.975 Time/s - 43.046 23.234 61.463 47.482 371.04 71.641 Case 2 MPSNR 14.118 33.175 34.291 34.710 35.294 34.251 36.232 MSSIM 0.2142 0.9332 0.9439 0.9457 0.9496 0.9348 0.9632 MFSIM 0.5976 0.9588 0.9627 0.9643 0.9710 0.9608 0.9757 MSAM 0.6687 0.0846 0.0547 0.0494 0.0625 0.1304 0.0519 ERGAS 707.93 75.787 65.678 61.582 59.506 108.85 52.485 Time/s - 43.294 22.994 61.906 44.786 397.41 72.997 Case 3 MPSNR 14.092 33.083 34.193 34.652 35.220 34.338 35.969 MSSIM 0.2114 0.9330 0.9437 0.9454 0.9491 0.9356 0.9619 MFSIM 0.5955 0.9587 0.9624 0.9641 0.9707 0.9618 0.9746 MSAM 0.6720 0.0855 0.0553 0.0493 0.0641 0.1207 0.0538 ERGAS 709.14 76.431 66.452 61.936 60.275 100.18 54.261 Time/s - 43.680 22.733 61.790 45.851 373.44 77.096 2）视觉质量比较
为了更好地观察这几个算法的去噪效果，将Case1中各种算法去噪后的第20波段的图像进行对比，如图 1所示。从(b)中可以看到图像已经完全被噪声污染，基本观察不到任何图像特征。经过各种算法去噪后，可以看到图像的基本特征，然而LRMR并不能完全去除噪声，如图(c)所示，图像中还有残留的噪声；LRTV和LRTDTV虽然能够去除噪声，但是会使图(d)和图(e)中的图像细节特征模糊；FRCTR-PnP也能取得较好的视觉效果，LRTDGS和NCLRGSTV能够取得最好的视觉效果，这是由于提出的NCLRGSTV算法可以对较大的奇异值减少收缩来保留图像更多的细节特征，增大对小的奇异值的收缩来去除噪声。如图 2所示，在增加了条带噪声后，其他算法仍然有上述问题，而NCLRGSTV算法中结合了群稀疏总变分正则项可以保持空间光谱的局部平滑性，从而抑制条带等稀疏噪声。
![]()
图 1 Case1中各种算法去噪后第20波段对比Figure 1. Comparison of the 20th band after denoising of various algorithms in Case1![]()
图 2 Case2中各种算法去噪后第58波段对比Figure 2. Comparison of the 58th band after denoising of various algorithms in Case23）定性分析
图 3、图 4和图 5为实验中各种算法在每个频段的PSNR和SSIM指标曲线，可以看到，在图 3和图 4中提出的NCLRGSTV算法在Pavia City Center数据集每个波段都达到最高值，这是由于采用对数张量核范数能够减少对较大奇异值的收缩，可以保留图像更多的细节特征，群稀疏总变分也可以在探索空间群稀疏性的同时增强相邻光谱间的相关性。虽然在图 5中该方法并没有在所有波段都达到最高值，但是在所有波段的平均值为最优。
3.2 真实数据实验
由于该模型对HSI中含有高斯噪声、椒盐噪声、条带和死线噪声等混合噪声的去除具有很好的效果，因此，实验中我们采用了具有类似真实噪声的HSI数据集HYDICE Urban进行性能评估[16]。与模拟实验一样，在测试恢复实验之前，将每个图像波段的灰度值归一化。
为了方便对比，用其中一个典型的波段来表示恢复结果。图 6表示第109波段的恢复结果，从图 6(a)中可以看到，图像已经完全被上述4种噪声污染。经过各种方法恢复后，大量噪声被去除，可以看图像的基本特征；但从图(b)、(c)和(f)中的放大图可以看到，仍有一些条纹未能被去除，而图(d)中的图像过于平滑，使图像的细节变得模糊；而提出的NCLRGSTV方法消除了更多的噪声，图像较为清晰，保留了图像很多的细节特征。
![]()
图 6 HYDICE Urban数据第109波段恢复图像比较Figure 6. Comparison of restored images in the 109th band of HYDICE Urban data为了更清楚地表示各种方法的恢复结果，如图 7所示，给出了恢复后的第109波段的垂直平均剖面图。纵轴表示每行的平均数字数值，横轴表示行数。如图 7(a)所示，受到噪声的影响曲线有强烈的波动，经过恢复后，曲线的波动明显减小。图 7(b)、(c)和(f)还有部分噪声未能去除，而(d)中有些部分过于平滑。相比之下，图(e)和(h)都可以获得更合理的平均轮廓结果。这也说明该模型去噪的有效性。
![]()
图 7 HYDICE Urban数据第109波段水平平均剖面比较Figure 7. Comparison of the 109th band horizontal mean profile in HYDICE Urban data3.3 参数分析与讨论
NCLRGSTV模型受到一些参数的影响，比如惩罚参数β和常数ε，期望的秩k，正则化参数λ1和λ2。像其他算法一样，将惩罚参数β=10－2设为初始值，并在每次迭代中更新为β=min(1.5β, 106β)，这种方法便于算法收敛，在实际中已被证明是有效的。常数ε的设置可以参考文献[17]，将其设置过大或过小都不能获得较好的结果，因此将其设置为中等大小，令ε=70。
图 8(a)为根据不同秩k得出的客观评价指标MPSNR。如图所示，该算法对期望的k比较敏感，当k=3时获得较大的值。因此，对期望秩k的估计越准确，恢复的结果就越好。图 8(b)是正则化参数λ1和λ2对PSNR的影响，如图所示，随着λ1和λ2增大，PSNR值也随之减小，而当λ1=0.01和λ2=0.02时，该算法可以取得较大值。因此，在实验中将其设置为k=3、λ1=0.01和λ2=0.02。
在文献[33]中，保证了NCLRGSTV模型理论上的收敛性。此外，图 9显示了所提出的NCLRGSTV方法的收敛曲线。超过20次迭代后，数值不再出现明显的变化，这表明该算法具有较高的收敛速度。
4. 结论
混合噪声的复杂性给高光谱图像处理和分析带来了巨大的挑战，因此，提出了一种非凸低秩张量分解和群稀疏总变分的高光谱图像恢复模型，该模型采用对数张量核范数利用高光谱图像的低秩特性时，减少对较大奇异值的收缩，保留更多的细节特征；同时结合群稀疏总变分正则化，增强HSI空间稀疏性和相邻光谱的相关性，达到去除高斯噪声、椒盐噪声、条带和死线等稀疏噪声。利用ADMM算法进行求解，实验中也易于收敛。与流行的恢复算法进行对比试验发现，该算法去除混合噪声具有一定的优势。
然而，该算法在处理HSI混合噪声中含有大量死线噪声时，仍有少量死线噪声残留难以被去除。除此之外，该算法还有需要改进的地方，如模型参数的自适应设置问题。未来可以利用即插即用框架嵌入一些深度去噪先验[34-35]，来增强模型去除噪声的能力。
-
![]()
图 1 Case1中各种算法去噪后第20波段对比
Figure 1. Comparison of the 20th band after denoising of various algorithms in Case1
![]()
图 2 Case2中各种算法去噪后第58波段对比
Figure 2. Comparison of the 58th band after denoising of various algorithms in Case2
![]()
图 6 HYDICE Urban数据第109波段恢复图像比较
Figure 6. Comparison of restored images in the 109th band of HYDICE Urban data
![]()
图 7 HYDICE Urban数据第109波段水平平均剖面比较
Figure 7. Comparison of the 109th band horizontal mean profile in HYDICE Urban data
表 1 Pavia city center数据集的不同去噪方法的定量评价结果
Table 1 Quantitative evaluation results of different denoising methods in Pavia city center data sets
Case Indexes Noise LRMR LRTV LRTDTV LRTDGS FRCTR-PnP NCLRGSTV Case 1 MPSNR 14.144 33.336 34.356 34.743 35.380 34.557 36.369 MSSIM 0.2143 0.9341 0.9444 0.9457 0.9506 0.9370 0.9637 MFSIM 0.5985 0.9590 0.9626 0.9646 0.9647 0.9630 0.9761 MSAM 0.6676 0.0833 0.0545 0.0495 0.0637 0.1331 0.0514 ERGAS 707.54 74.698 65.280 70.351 61.441 109.32 51.975 Time/s - 43.046 23.234 61.463 47.482 371.04 71.641 Case 2 MPSNR 14.118 33.175 34.291 34.710 35.294 34.251 36.232 MSSIM 0.2142 0.9332 0.9439 0.9457 0.9496 0.9348 0.9632 MFSIM 0.5976 0.9588 0.9627 0.9643 0.9710 0.9608 0.9757 MSAM 0.6687 0.0846 0.0547 0.0494 0.0625 0.1304 0.0519 ERGAS 707.93 75.787 65.678 61.582 59.506 108.85 52.485 Time/s - 43.294 22.994 61.906 44.786 397.41 72.997 Case 3 MPSNR 14.092 33.083 34.193 34.652 35.220 34.338 35.969 MSSIM 0.2114 0.9330 0.9437 0.9454 0.9491 0.9356 0.9619 MFSIM 0.5955 0.9587 0.9624 0.9641 0.9707 0.9618 0.9746 MSAM 0.6720 0.0855 0.0553 0.0493 0.0641 0.1207 0.0538 ERGAS 709.14 76.431 66.452 61.936 60.275 100.18 54.261 Time/s - 43.680 22.733 61.790 45.851 373.44 77.096 
下载: 导出CSV
-
[1] Bioucas Dias J M, Plaza A, Camps Valls G, et al. Hyperspectral remote sensing data analysis and future challenges[J]. IEEE Trans. Geosci. Remote Sens. , 2013, 1: 6-36.
[2] LI Hengchao, WANG Weiye, PAN Lei, et al. Robust capsule network based on maximum correntropy criterion for hyperspectral image classification[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2020, 13: 738-751. DOI: 10.1109/JSTARS.2020.2968930
[3] DIAN R, LI S, FANG L. Learning a low tensor-train rank representation for hyperspectral image super-resolution[J]. IEEE Trans. Neural Netw. Learn. , Syst. , 2019, 30(9): 2672-2683. DOI: 10.1109/TNNLS.2018.2885616
[4] Bioucas Dias J M, Plaza A, Dobigeon N, et al. Hyperspectral unmixing overview: geometrical, statistical, and sparse regression-based approaches[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2012, 5(2): 354-379. DOI: 10.1109/JSTARS.2012.2194696
[5] ZHANG Pengdan, NING Jifeng. Hyperspectral image denoising via group sparsity regularized hybrid spatio-spectral total variation[J]. Remote. Sens., 2022, 14: 2348. DOI: 10.3390/rs14102348
[6] Elad M, Aharon M. Image denoising via sparse and redundant representations over learned dictionaries[J]. IEEE Trans. Image Process., 2006, 15(12): 3736-3745. DOI: 10.1109/TIP.2006.881969
[7] Green A A, Berman M, Switzer P, et al. A transformation for ordering multispectral data in terms of image quality with implications for noise removal[J]. IEEE Transactions on Geoscience and Remote Sensing, 1988, 26(1): 65-74. DOI: 10.1109/36.3001
[8] Candès E J, LI X, MA Y, et al. Robust principal component analysis?[J]. J. ACM, 2011, 58: 1-37.
[9] GU S, ZHANG L, ZUO W, et al. Weighted nuclear norm minimization with application to image denoising[C]//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2014: 2862–2869.
[10] HE W, ZHANG H, ZHANG L, et al. Hyperspectral image denoising via noise-adjusted iterative low-rank matrix approximation[J]. IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens, . 2015, 8: 3050-3061.
[11] Karami A, Yazdi M, Asli A Z. Noise reduction of hyperspectral images using kernel non-negative tucker decomposition[J]. IEEE J. Sel. Topics Signal Process. , 2011, 5(3): 487-493. DOI: 10.1109/JSTSP.2011.2132692
[12] FAN H, CHEN Y, GUO Y, et al. Hyperspectral image restoration using low-rank tensor recovery[J]. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. , 2017, 10(10): 4589-4604. DOI: 10.1109/JSTARS.2017.2714338
[13] HE W, ZHANG H, ZHANG L, et al. Total-variation-regularized low-rank matrix factorization for hyperspectral image restoration[J]. IEEE Trans. Geosci. Remote Sens. , 2016, 54(1): 178-188. DOI: 10.1109/TGRS.2015.2452812
[14] WANG Y, PENG J, ZHAO Q, et al. Hyperspectral image restoration via total variation regularized low-rank tensor decomposition[J]. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., 2018, 11(4): 1227-1243. DOI: 10.1109/JSTARS.2017.2779539
[15] FAN H, LI C, GUO Y, et al. Spatial–spectral total variation regularized low-rank tensor decomposition for hyperspectral image denoising[J]. IEEE Trans. Geosci. Remote Sens. , 2018, 56(10): 6196-6213. DOI: 10.1109/TGRS.2018.2833473
[16] CHEN Y, HE W, Yokoya N, et al. Hyperspectral image restoration using weighted group sparsity-regularized low-rank tensor decomposition[J]. IEEE Transactions on Cybernetics, 2020, 50(8): 3556-3570. DOI: 10.1109/TCYB.2019.2936042
[17] SHI Q, TANG X, YANG T, et al. Hyperspectral image denoising using a 3-D attention denoising network[J]. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59(12): 10348-10363. DOI: 10.1109/TGRS.2020.3045273
[18] WANG Z, SHAO Z, HUANG X, et al. SSCAN: a spatial–spectral cross attention network for hyperspectral image denoising[J]. IEEE Geoscience and Remote Sensing Letters, 2022, 19: 1-5.
[19] ZHENG Y B, HUANG T Z, ZHAO X L, et al. Mixed noise removal in hyperspectral image via low-fibered-rank regularization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 58(1): 734-749. DOI: 10.1109/TGRS.2019.2940534
[20] Fazel M, Hindi H, Boyd S P. Log-det heuristic for matrix rank minimization with applications to Hankel and Euclidean distance matrices[C]//Proc. IEEE Amer. Control Conf., 2003, 3: 2156-2162.
[21] Mohan K, Fazel M. Iterative reweighted algorithms for matrix rank minimization[J]. J. Mach. Learn. Res., 2012, 13(1): 3441-3473.
[22] DENG Y, DAI Q, LIU R, et al. Low-rank structure learning via nonconvex heuristic recovery[J]. IEEE Trans. Neural Netw. Learn. Syst., 2013, 24(3): 383-396. DOI: 10.1109/TNNLS.2012.2235082
[23] DONG W, SHI G, HU X, et al. Nonlocal sparse and low-rank regularization for optical flow estimation[J]. IEEE Trans. Image Process. , 2014, 23(10): 4527-4538. DOI: 10.1109/TIP.2014.2352497
[24] JI T Y, HUANG T Z, ZHAO X L, et al. A non-convex tensor rank approximation for tensor completion[J]. Applied Mathematical Modelling, 2017, 48: 410-422. DOI: 10.1016/j.apm.2017.04.002
[25] Boyd S, Parikh N, CHU E, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers[J]. Found. Trends Mach. Learn., 2011, 3(1): 1-122.
[26] ZHANG Z, LIU D, Aeron S, et al. An online tensor robust PCA algorithm for sequential 2D data[C]//Proc. Int. Conf. Acoust., Speech, Signal Process., 2016: 2434-2438.
[27] LU C, FENG J, CHEN Y, et al. Tensor robust principal component analysis: exact recovery of corrupted low-rank tensors via convex optimization[C]//IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016: 5249-5257.
[28] GU S, ZHANG L, ZUO W, et al. Weighted nuclear norm minimization with application to image denoising[C]//Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR), 2014: 2862-2869.
[29] Aggarwal H K, Majumdar A. Hyperspectral image denoising using spatio-spectral total variation[J]. IEEE Geosci. Remote Sens. Lett. , 2016 13(3): 442-446.
[30] HE W, ZHANG H, SHEN H, et al. Hyperspectral image denoising using local low-rank matrix recovery and global spatial–spectral total variation[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2018, 11(3): 713-729. DOI: 10.1109/JSTARS.2018.2800701
[31] XIE Q, ZHAO Q, MENG D, et al., Multispectral images denoising by intrinsic tensor sparsity regularization[C]//Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR), 2016: 1692-1700.
[32] LIU Y Y, ZHAO X L, ZHENG Y B, et al. Hyperspectral image restoration by tensor fibered rank constrained optimization and plug-and-play regularization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 1-17.
[33] Eckstein J, YAO W. Understanding the convergence of the alternating direction method of multipliers: theoretical and computational perspectives[J]. Pacific Journal of Optimization, 2015, 11(4): 619-644.
[34] LUO Y S, ZHAO X L, JIANG T X, et al. Hyperspectral mixed noise removal via spatial-spectral constrained unsupervised deep image prior[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14: 9435-9449. DOI: 10.1109/JSTARS.2021.3111404
[35] LAI Z, WEI K, FU Y. Deep plug-and-play prior for hyperspectral image restoration[J]. Neurocomputing, 2022, 481: 281-293. DOI: 10.1016/j.neucom.2022.01.057
计量
- 文章访问数: 0
- HTML全文浏览量: 0
- PDF下载量: 0
