Hyperspectral image denoising method based on depth image prior

利用遥感技术获取的高光谱图像，作为一个数据立方体，结合了用于表征物质反射特性的光谱信息以及用于刻画物质2维空间分布的图像信息^[1]。凭借其强大的空间与光谱信息的表达能力，高光谱图像被广泛应用于物质分类^[2]、环境监测、军事目标的检测与分析等^[3]领域。高光谱图像在采集与处理过程中容易受到各种噪声的干扰^[4]，这对后续的高光谱图像的应用，如高光谱图像的分类^[5]、目标检测^[6]以及目标分割^[7]等，都有一定的影响。因此，在有效去除噪声的基础上，恢复重建高分辨率的图像成为遥感图像领域研究的热点。

基于先验的去噪方法可以精确地表达出高光谱图像的空谱信息分布特点。依据空谱特征描述方式的不同，先验方法可以更为精确地划分为基于平滑的去噪方法^[8]、基于全局自相似的去噪方法^[8]、基于低秩的去噪方法^[9]以及基于深度学习的去噪方法。本文作者^[10]曾通过构建空谱平滑图正则的方式, 实现基于塔克分解的低秩模型算法的扩展，但现有的基于平滑的方法容易陷入过平滑造成的细节信息丢失。ZENG等人^[11]通过张量分解(block term decomposition，BTD) 的手段深入探究相似立体补丁组的光谱局部联系与空间的全局相关性，但基于全局自相似的方法采用的块匹配方法往往涉及较大计算量，会造成运算成本加剧。WANG等人^[12]构建了一个融合了空谱全变分的低秩塔克分解模型，可以充分利用高光谱图像的空谱全局与局部相关性并恢复出高质量的去噪图像。然而，基于矩阵模型的算法面临着立体结构信息被破坏的问题。XIONG等人^[9]尝试将全局注意力机制模块引入到传统的U型网络(U-Net)，可以实现同时捕获空间的全局与局部信息。上述基于监督学习的方法通常涉及大量的数据训练以及超参数，会导致计算负担的增加。

变换域的出现为高光谱图像的去噪研究开辟了一条全新的道路。图像经变换后, 主要信息会集中到低频部分，细节信息则在高频部分^[13]。针对变换域的特点从张量秩^[14]、平滑^[15]、深度学习^[16]等不同角度对去噪方法进行了探索。KONG等人^[14]研究了基于紧框架的张量纤维秩以及对应的凸松弛来表达三模张量核范数，可以更好地刻画高光谱图像不同维度的低秩性。但总体来看，基于变换域的张量核范数方法由于涉及到大数据的奇异值分解计算，面临着计算复杂度过高的问题。XU等人^[15]尝试通过非负矩阵分解得到丰度矩阵在曲线波域的稀疏特性，但基于变换域平滑的方法不能充分利用不同频段的图像特征信息。考虑到小波变换可以利用不同频段，实现高斯噪声与条带噪声的有效分离，PAN等人^[16]设计了多层3维卷积低频网络和结合渐进式空谱混合卷积块的高频卷积神经网络，可以较好地学习变换域下图像不同频段的空谱潜在特征，然而基于变换域的深度学习模型仍存在着可解释性差以及网络复杂度高等问题。

为了克服上述问题，本文作者提出将基于变换域的正则与深度学习模型相结合，构建无监督学习的噪声去除模型。首先，对无监督去噪模型展开探索，规避了监督去噪模型缺乏大量高光谱训练数据的缺陷，同时通过增加传统高光谱图像的先验优化改善了基于深度图像先验涉及到的半拟合现象；然后通过构建端元矩阵与丰度矩阵在变换域中的稀疏正则实现空谱局部平滑；最后, 设计迭代算法求解模型得到高质量的恢复图像，同时加速网络学习并获得更好的性能。

假设F ∈ R_m×n×p表示高光谱立体图像，R表示实数集，利用矩阵分解可将高光谱立体数据F重塑为矩阵形式Y =[y₁, …, y_i, …, y_c]，其中y_i∈ R_p被定义为第i个像素对应的光谱向量。由于纯高光谱数据可以利用线性混合模型分解为端元矩阵与丰度矩阵的乘积，则噪声扰动的退化模型可以定义为：

式中: A ∈ R_p×r表示包含r个端元和p个光谱波段的端元矩阵; X ∈ R_r×c是丰度矩阵; c=m×n代表空间像素数目，其每一列表示对应像素包含的r个端元的分布比例；S ∈ R_p×c表示稀疏噪声；N ∈ R_p×c代表高斯噪声或误差。

在高光谱图像退化模型式(1)的基础上，引入各变量的正则约束定义基于高斯噪声最小化的去噪模型：

式中: ‖·‖_F为Frobenius范数; λ_Aθ_A(A)、λ_Xθ_X(X)和λ_Sθ_S(S)分别代表着端元矩阵、丰度矩阵以及稀疏噪声的正则项。利用深度图像先验分别学习端元矩阵与丰度矩阵的特征，得到空谱深度图像先验去噪模型：

式中: f_ξ1(·)与f_ξ2(·)分别表示端元矩阵与丰度矩阵的深度图像先验神经网络；ξ₁和ξ₂代表对应的可学习网络权重；Z₁∈ R_p×r及Z₂∈ R_r×mn用于生成空谱矩阵的随机输入。

无监督深度图像先验是一个类似于U-Net的带有跳跃连接的沙漏类型卷积神经网络，将随机生成的张量视作网络输入，通过不断调整网络参数值实现图像特征信息的学习。它的特点是: 不需要训练网络用于训练，可以规避监督去噪模型缺乏大量高光谱训练数据的缺陷。基于此，本文作者采用深度图像先验分别学习丰度矩阵与端元矩阵的特征信息，具体的网络结构如图 1所示。图中，upsample是上采样；BatchNorm2D是深度学习中的一种层，用于进行批归一化操作；MaxPooling是最大池化。

Figure 1.  Overall network flowchart of deep image prior

设计的深度图像先验网络结构分别包括5层编码层、解码层以及跳跃链接。编码层的每层设计了两组相同的卷积、批标准化与激活模块，可以更好地提取图像的特征信息，其中2维卷积(convolution 2-D，Conv2D)采用大小为3、步长为1的卷积核，带泄露修正线性单元(leaky rectified linear unit，LeakyReLU)作为激活模块可以加速网络的学习过程；解码层通过一系列的上采样、卷积、批标准化与激活等操作，用于解析学习到的图像特征信息；跳跃连接采用大小为1、步长为1的Conv2D模块以及批标准化与激活模块建立相同层级间的编解码特征之间的联系，可以减少下采样操作所引起的信息损失。

对于去噪模型式(3)，如果仅使用空间与光谱深度图像先验去学习空谱特征，则会导致半拟合现象，造成去噪性能的下降。因此本文中对稀疏噪声施加L₁范数正则用于增强其稀疏特性，对端元矩阵与丰度矩阵的局部相关性构建紧框架域平滑正则。综上所述，将去噪模型式重新定义为：

式中: W表示紧框架变换矩阵，对于多网络结构优化问题，可以交替迭代ξ₁和ξ₂网络参数以及稀疏噪声S直到收敛。

用于生成端元矩阵与丰度矩阵的ξ₁和ξ₂神经网络优化子问题如下：

其可以利用动量梯度下降算法求解。

式中: ∇_ξ1(·)与∇_ξ2(·)分别表示运用反向传播算法计算得到的ξ₁和ξ₂梯度；α对应着通过自适应矩估计算法确定的步长；k是当前迭代的次数。由此，在更新了网络参数ξ₁和ξ₂的基础上，可以得到对应的网络输出f_ξ1(Z₁)^(k+1)与f_ξ2(Z₂)^(k+1)。

针对稀疏噪声S的优化子问题如下：

其闭式解可以利用软阈值算子得到：

本文作者所提出的基于空谱深度图像先验结合紧框架平滑与噪声稀疏先验(spatial spectral depth tight sparse，SSDTS)算法的迭代流程描述如下：

(a) 输入Y，λ_A，λ_X，λ_S；

(b) 随机初始化A⁽⁰⁾和X⁽⁰⁾，S⁽⁰⁾=0，k=0；

(c) 当不满足迭代停止条件时，借助式(6)更新ξ₁^(k+1), 借助式(7)更新ξ₂^(k+1), 借助式(9)更新S^(k+1), k=k+1；

(d) 判断收敛条件为‖ξ₁^(k+1)ξ₂^(k+1)-ξ₁^(k)ξ₂^(k)‖_F²/ ‖ξ₁^(k)ξ₂^(k)‖_F²≤0.001；

(e) 结束循环；

(f) 输出L =f_ξ1(Z₁)f_ξ2(Z₂)。

本文中采用模拟数据集和真实数据集测试和评估SSDTS算法对高光谱图像的去噪能力。在此基础上，选取4种广泛使用的算法作为对比算法，包括低秩矩阵恢复(low-rank matrix recovery，LRMR)^[17]、全变分正则化的低秩矩阵分解(the TV-regularized low-rank matrix factorization，LRTV)^[18]、利用子空间和深度图像先验的高光谱图像混合噪声去除方法(hyperspectral image mixed noise removal method using subspace representation and deep CNN image prior，HySuDeep)^[19]和2维深度图像先验(deep image prior 2D，DIP2D)^[20]。与此同时，为了准确可靠地评价图像重建性能，采用4个定量评价指标，分别为平均峰值信噪比(mean peak signal-to-noise ratio，MPSNR)、平均结构相似度(mean structural similarity index measure，MSSIM)、光谱角制图(spectral angle mapper，SAM)和相对无量纲全局误差(erreur relative globale adimensionnelle de synthèse，ERGAS)。另外，作为实验的预处理，将高光谱图像的灰度值归一化到[0, 1]区间。上述所有深度学习的实验结果均在GeForce GTX 2080 Ti处理器上使用Python环境的pytorch库产生。

对于模拟实验，本文作者设计了两个高光谱图像数据集作为参考。一个是Washington DC Mall的子图像，尺寸为120×256×256个像素; 另一个是Salinas数据集，包含192×192个空间像素和130个光谱波段。

方案Ⅰ(高斯噪声与椒盐噪声)：采用分布比例为0.1的高斯噪声与椒盐噪声污染高光谱图像的所有波段。

方案Ⅱ(高斯噪声与条带噪声)：对于受到与方案Ⅰ相同的高斯噪声污染的高光谱图像，随机选择波段添加条带噪声，其波段数量占所有波段数目的0.5，每个波段中条带的数量从5~15间随机变化。

方案Ⅲ(高斯噪声与截止噪声)：包含与方案Ⅰ相同的高斯噪声和随机选择的波段上的截止噪声，其波段数量占所有波段数目的0.5。与此同时，截止噪声的宽度从1~3随机变化，每个波段存在的截止噪声总数在^{[6, 10]}之间。

方案Ⅳ(高斯噪声与混合稀疏噪声)：在每个高光谱图像波段中，添加与方案Ⅰ、方案Ⅱ和方案Ⅲ相同的高斯噪声、椒盐噪声、条带噪声和截止噪声。

如图 2所示，Washington DC Mall图像遭受到严重的高斯噪声与椒盐噪声的污染，LRTV和DIP2D算法生成的去噪图像仍残存着明显的噪声，LRMR和HySuDeep算法去除噪声效果较好，但存在着图像细节信息的丢失，本文中提出的SSDTS算法可以最大程度地消除噪声，并准确地描述图像的细节。对比图 3可以发现，本文中所提出的SSDTS算法还具有较为理想的去噪结果，LRMR、LRTV和DIP2D算法对于高斯噪声与多种稀疏噪声的混合具有较差的抑制性，HySuDeep算法对应的重构图像面临着过度平滑的现象。综上所述，SSDTS算法在各种复杂、严重的噪声干扰下均表现出较好的视觉恢复性能。

Figure 2.  The denoised results of the 100th band in Washington DC Mall dataset under case Ⅰ

Figure 3.  The denoised results of the 100th band in Salinas dataset under case Ⅳ

表 1和表 2中对比了SSDTS算法与4种算法在Washington DC Mall和Salinas数据集上产生的评价指标值。相对于其它算法，SSDTS算法产生的平均峰值信噪比值提高了1 dB以上，对应的ERGAS值至少降低了10%。这说明本文中提出的方法在严重的高斯噪声或稀疏噪声的污染下，依然能获得更为理想的重构图像。从整体结果来看，SSDTS算法在几乎所有评价指标上都取得了令人满意的性能。

图 4和图 5中绘制了SSDTS算法与所有对比算法在受到不同噪声方案污染的Washington DC Mall和Salinas数据集上产生的峰值信噪比(peak signal-to-noise ratio，PSNR)和ERGAS曲线。值得注意的是，本文中提出的SSDTS算法产生的PSNR曲线优于由对比算法生成的其它曲线；而SSDTS算法对应的ERGAS曲线在绝大多数波段得到最低值。上述性能曲线变化情况表明，SSDTS去噪算法在大部分波段获得了最理想的结果。同时，也从侧面证明，不同于仅采用深度图像先验的DIP2D算法，本文中所提出的空谱深度图像先验结合紧框架域平滑正则在图像去噪领域中可以产生更加满意的结果。

Figure 4.  PSNR and ERGAS values of all bands in Washington DC Mall dataset under case Ⅰ and Ⅲ

Figure 5.  PSNR and ERGAS values of all bands in Salinas dataset under case Ⅱ and Ⅳ

为了验证本文作者所提出的算法对于真实噪声的去除仍具有有效性，选取由传感器直接采集的WHU HongHu数据集，包含150个光谱波段和192×192个空间像素。图 6展示出经过所有去噪方法的优化处理后的WHU HongHu数据集在第2个波段的重构结果。可以明显观察到，HySuDeep与DIP2D算法对真实噪声的抑制效果有限。其它两个对比算法可以实现较好地去噪，但会造成图像局部特征信息的模糊。相比之下，由于深度图像先验模型与传统的变换域平滑先验的结合，本文作者提出的算法能够在最大限度地消除噪声的同时, 精确地保留图像细节信息。

Figure 6.  The denoised results of the 2nd band in WHU HongHu dataset

由于缺乏干净图像作为参考，采用水平均值剖面图作为图像质量评价指标，通过遥感影像平均像元亮度值(digital number，DN)记录地物的灰度值。仔细对比图 7中各个剖面图的变化趋势可知，由LRMR、LRTV、HySuDeep与DIP2D这4种对比算法产生的剖面图在局部位置仍存在明显的波动，这说明对应的去噪图像中存在着噪声残留。而本文作者提出的SSDTS算法对应的水平均值剖面图平滑效果最好。

Figure 7.  Horizontal mean profile of the 2nd band in WHU HongHu dataset

图 8反映了方案Ⅱ下的模拟数据集的MPSNR与λ_A、λ_X、λ_S和端元数目r的关系。由图可知，当λ_A的值设为0.1、λ_X的值设为0.2、λ_S值设为0.1、r选择为5时，在Washington DC Mall与Salinas数据集上产生的MPSNR值均达到了最优性能的近似值。

Figure 8.  Relationship between parameter values and MPSNR

从图 9中可以明显观察到，随着迭代次数的增加，MPSNR对应的曲线均逐渐趋向于平稳的理想结果，这表明本文作者提出的SSDTS算法具有理想的收敛性，同时也从侧面证明SSDTS算法中采用的噪声稀疏与变换域平滑正则项可以有效改善深度图像网络存在的半拟合问题，在充分学习图像特征达到最佳评价指标结果后, 没有产生继续学习噪声特征导致性能下降的现象。

Figure 9.  Relationship between iteration number and MPSNR in different cases

为了验证本文作者提出的将深度学习与传统先验结合的有效性，设计了空谱深度图像先验(spatial spectral depth image prior，SSDIP)算法、空谱深度图像先验结合稀疏噪声正则(spatial spectral depth sparse noise，SSDSN) 算法、空谱深度图像先验结合紧框架平滑先验正则(spatial spectral depth frame smoothing，SSDFS)算法和提出的SSDTS算法的对比实验。从表 3可知，在受到不同程度混合噪声干扰下，相对于其它3种算法，本文作者提出的SSDTS算法均呈现出更加理想的去噪表现，这证明了本文中将空谱深度图像先验与紧框架平滑以及噪声稀疏先验相结合的合理性。

从表 4中的图像去噪运行时间来看，LRMR算法恢复图像最耗时，HySuDeen和本文作者所提出的算法运行时间很接近，两者相差非常小，这证明了本文作者提出的算法不止在客观评价指标和视觉恢复效果上有良好的表现，在运行时间上也具有一定的优势。

为解决深度图像的半拟合问题，改善矩阵分解对立体特征的破坏，本文作者利用特定的深度图像先验网络分别学习矩阵分解得到的更加明显的端元与丰度矩阵上特征信息；同时，构建端元与丰度矩阵上的紧框架平滑正则可以有效改善网络的半拟合现象并加速网络训练; 最后，设计了高效迭代算法来实现神经网络与优化问题相结合模型的高效求解。多组模拟高光谱图像数据实验结果表明，所提出的SSDTS算法可以实现MPSNR性能1 dB以上的提升，ERGAS性能有至少10%的降低，进一步验证了本文中算法在混合噪声去除方面具有优异的表现。