一种基于非凸低秩矩阵分解的高光谱图像去噪方法与流程

本发明隶属于计算机图像、信息处理技术、遥感技术等领域，尤其涉及一种基于非凸低秩矩阵分解的高光谱图像去噪方法。

背景技术：

随着数字遥感技术的快速发展，高光谱图像广泛应用于军事侦察、环境科学、地质勘探、生物医学成像等领域。由于传感器的机械故障，图像传输故障等等各种因素的影响，高光谱图像在采集与传输过程中不可避免地会受到各种噪声，例如高斯噪声、脉冲噪声、条纹等的污染，严重制约了高光谱图像的进一步应用。同时，高光谱图像维数的急剧增长，导致了“维数灾难”。因此，设计一种既可以去除各种混合噪声，同时可以解决“维数灾难”的高光谱图像去噪方法成为了遥感图像应用中关键技术之一。

传统的高光谱去噪方法往往有以下缺陷：1)仅仅考虑高斯白噪声，但是在高光谱图像中往往含有多种混合噪声；2)仅仅用到单波段图像信息，没有考虑高光谱图像中丰富的地物空间特性和光谱特性。近年来，随着压缩感知理论的蓬勃发展，基于低秩矩阵分解的高光谱图像去噪方法越来越受到广大研究者的关注。此类算法多以核范数凸近似矩阵秩函数。然而，基于核范数的低秩矩阵分解算法在高光谱图像去噪中往往出现以下几种问题：

1、充分利用了高光谱图像丰富的光谱信息，但是只能有效地去除高斯噪声和部分稀疏噪声；

2、基于核范数的凸模型在实际高光谱图像去噪中实际效果差，并伴随着高光谱图像维数的增长，计算时间显著延长；

3、基于核范数的凸模型往往会导致矩阵秩估计过大，导致条带噪声无法去除；

4、基于双边随机投影的低秩矩阵分解算法，往往需要提前给出矩阵的秩等先验信息。

参考文献：

[1]Dabov K,Foi A,Katkovnik V,et al.Image denoising by sparse 3-D transform-domain collaborative filtering[J].IEEE Transactions on Image Processing,2007,16(8):2080-95.

[2]Zhang H,He W,Zhang L,et al.Hyperspectral Image Restoration Using Low-Rank Matrix Recovery[J].IEEE Transactions on Geoscience and Remote Sensing,2014,52(8):4729-4743.

[3]He W,Zhang H,Zhang L,et al.Hyperspectral Image Denoising via Noise-Adjusted Iterative Low-Rank Matrix Approximation[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2015,8(6):1-12.

[4]Xie Y,Qu Y,Tao D,et al.Hyperspectral Image Restoration via Iteratively Regularized Weighted Schatten,p-Norm Minimization[J].I EEE Transactions on Geoscience and Remote Sensing,2016,54(8):4642–4659.

技术实现要素：

发明目的：本发明旨在解决遥感图像领域中高光谱图像去噪的技术问题，提供一种基于非凸低秩矩阵分解的高光谱图像去噪方法，从而有效地去除高光谱图像存在的混合噪声的同时保存图像中光谱和空间特征。

技术方案：本发明公开了一种基于非凸低秩矩阵分解的高光谱图像去噪方法，其核心在于利用非凸低秩矩阵分解技术线性化分离高光谱图像，包括以下步骤：

步骤1，输入含有混合噪声的原始高光谱图像，高光谱图像是由几十乃至数百个连续波段图像组成的三维图像数据。初始化子块的尺寸，扫描的步长，然后把原始图像分割成重叠的子块，以保留高光谱图像的局部细节；

步骤2，生成矩阵：对每个子块的第k个波段向量化，然后对所有向量化的波段根据字典顺序排列成一个矩阵；

步骤3，非凸低秩矩阵分解：通过非凸低秩矩阵分解处理，将步骤2得到的含有混合噪声的矩阵，线性化分解成一个低秩矩阵，一个稀疏矩阵，以及一个无结构化矩阵；此过程中，以矩阵秩函数和稀疏矩阵l_2,1范数的线性和为目标函数，以含有混合噪声矩阵的线性化分解为约束条件，然后利用增广拉格朗日方法求解此问题；

步骤4，还原子块：将低秩矩阵的每个子列合成一个矩阵，将得到的所有矩阵按光谱维排列成三维子块，并对重叠部分图像像素值均值化处理；

步骤5，还原高光谱图像：将均值化处理的三维子块排列成和原始高光谱图像相同的三维图像数据，最后得到去噪后的高光谱图像。

本发明中，步骤1包括以下步骤：

步骤1-1，输入原始含有混合噪声的高光谱图像，此图像是由几十乃至数百个连续波段图像组成的三维图像数据d，其中d_h,d_w,d_s分别代表高光谱图像的高度，宽度和波段数，表示实数域中一个维数分别为d_h,d_w,d_s的三阶张量；

步骤1-2，初始化子块的尺寸为p，扫描步长为s，其中1≤s<p。将原始高光谱图像分割成重叠的子块，共计个。

本发明中，步骤2包括以下步骤：

步骤2-1，对于选取的第i个子块d_i，表示实数域中一个维数分别为p,p,d_s的三阶张量，将第k个波段R^p×p向量化为表示实数域中一个维数为p²的向量，对所有向量化的波段根据字典顺序排列成矩阵D_i，表示实数域中一个维数分别为p²,d_s的矩阵；

步骤2-2，重复步骤2-1，共计次。

本发明中，步骤3包括以下步骤：

步骤3-1，由于在对目标的空间特征成像的同时，对每个空间像元经过色散形成几十个乃至几百个窄波段以进行连续的光谱覆盖，因此干净的高光谱图像在光谱维上表现出低秩的特性。同时考虑到高光谱图像中含有的噪声，例如：高斯噪声，脉冲噪声，条纹，离群点等等，因此，本发明考虑以下的高光谱图像去噪模型：

D_i＝L_i+S_i+N_i，

其中L_i,S_i,N_i分别为低秩矩阵、稀疏矩阵和无结构化矩阵，

步骤3-2，通过下式进行鲁棒性主成分分析：

s.t.D_i＝L_i+S_i+N_i，

其中，s.t.表示满足条件，min表示最小化函数，||L_i||*表示低秩矩阵L_i的核范数，||S_i||₁表示稀疏矩阵S_i的1范数，λ是正则化参数。然而由于核范数平等地极小化所有的奇异值，导致不能很好地近似秩函数。本发明考虑如下非凸矩阵秩近似：

其中，||L_i||_γ表示低秩矩阵L_i的非凸矩阵秩近似，Σ代表求和运算，σ_j(L_i)表示低秩矩阵L_i的第j个奇异值，γ是一个大于0的常数，代表非凸稀疏罚函数，本发明中采用以下两个非凸稀疏罚函数：

或者

同时考虑到高光谱图像中多含有纵向条纹，因此采用稀疏矩阵S_i的2,1范数即以去除条带噪声，其中(S_i)_mn表示变量S_i的第m行第n列元素，表示开平方；

步骤3-3，采用非凸低秩矩阵分解模型：

s.t.D_i＝L_i+S_i+N_i，

利用增广拉格朗日方法求解上述模型。其中此模型的增广拉格朗日函数如下：

其中Λ_i是拉格朗日乘子，ρ是罚参数，||·||_F是矩阵的F范数，<·>表示矩阵的内积，然后依次迭代更新变量(L_i,S_i,N_i,Λ_i)，第k+1次迭代格式如下：

和分别表示变量L_i经过k+1次迭代后的值、变量S_i经过k+1次迭代后的值、变量N_i经过k+1次迭代后的值和变量Λ_i经过k+1次迭代后的值，argmin表示子问题最小化点；

步骤3-4，通过下式求解变量L_i经过k+1次迭代后的值

其中由于||L_i||_γ是非凸的，容易证明是连续的、凹的、光滑的、可微的并且在区间[0,+∞)是非减的。记m＝min{p²,d_s}，然后对||L_i||_γ在点处一阶泰勒展开，则上式变为：

其中，表示变量L_i在第k次迭代的第j个奇异值，是非凸稀疏罚函数的导数。然后得到的封闭解：

其中T_i^k＝U∑V^T是矩阵T_i^k的奇异值分解，U,V分别表示矩阵T_i^k的左酉矩阵和右酉矩阵，V^T表示右酉矩阵的转置，且max表示最大化，diag表示将向量变成一个对角矩阵；

步骤3-5，通过下式求解变量S_i经过k+1次迭代后的值

其中则的封闭解的第j列是：

其中||·||₂表示向量的2范数；

步骤3-6，通过下式求解变量N_i经过k+1次迭代后的值

此问题是一个最小二乘问题，的封闭式解如下：

步骤3-7，通过下式更新变量Λ_i经过k+1次迭代后的值

步骤3-8，重复步骤3-4至步骤3-7，共计次。

本发明中，步骤4包括以下步骤：

步骤4-1，对低秩矩阵选取矩阵L_i的第a列，合成一个p×p矩阵，将所有低秩矩阵按光谱维排列成三维子块

步骤4-2，重复步骤4-1，共计次；

步骤4-3，计算重叠部分图像像素值，本发明均值化处理重叠部分图像像素值。

本发明中，步骤5包括：

将均值化处理的子块排列成和原始高光谱图像相同的三维d_h×d_w×d_s图像数据，最后得到去噪后的高光谱图像。

有益效果：

1)本发明中的高光谱去噪方法去噪速度快，更加鲁棒性。对原始高光谱图像分割成重叠的子块后，可以并行化实现非凸低秩矩阵分解，显著地减少方法的运行时间，并将原始图像线性化分割成干净的高光谱图像，稀疏噪声和高斯随机噪声。

2)本发明中的高光谱去噪方法充分利用了高光谱图像丰富的光谱信息，同时考虑到不同类型的噪声，更加有效地去除混合噪声。由于将原始高光谱图像分割成重叠的子块，并对重叠部分像素值均值化处理，有效地保持高光谱图像中的局部细节信息。

3)本发明中的高光谱图像去噪方法具有较强的可扩展性和灵活性。理论上本发明中的方法对任意大小的高光谱图像都可以进行去噪。

附图说明

图1为本发明高光谱图像去噪方法的基本流程图。

图2为HYDICE urban dataset数据集去噪效果图，轻度噪声污染的第87波段去噪效果。

图3为HYDICE urban dataset数据集去噪效果图，中度噪声污染的第104波段去噪效果。

图4为HYDICE urban dataset数据集去噪效果图，轻度噪声污染的第108波段去噪效果。

图5为EO-1Hyperion Australia dataset数据集去噪效果图，轻度噪声污染的第61波段去噪效果。

图6为EO-1Hyperion Australia dataset数据集去噪效果图，轻度噪声污染的第123波段去噪效果。

图7为EO-1Hyperion Australia dataset数据集去噪效果图，中度噪声污染的第88波段去噪效果。

图8为EO-1Hyperion Australia dataset数据集去噪效果图，重度噪声污染的第51波段去噪效果。

具体实施方式

下面结合附图和具体实施方式对本发明做更进一步的具体说明，但本发明的应用范围不限于此：

本方法的流程图如图1所示，可主要分为四大过程：首先将原始高光谱图像分割成重叠的子块，并将各个子块矩阵化；其次进行并行化的非凸低秩矩阵分解；再次将低秩矩阵还原成三维数据块；最后将重叠部分像素值均值化，并还原成去噪后的高光谱图像。

具体地说，如图1所示，本发明公开了一种基于非凸低秩矩阵分解的高光谱图像去噪方法，主要包括以下几个步骤：

步骤1，分割高光谱图像：高光谱图像是由几十乃至数百个连续波段图像组成的三维图像数据。初始化子块的尺寸，扫描的步长，然后把原始图像分割成重叠的子块，以保留高光谱图像的局部细节；

步骤2，生成矩阵：对每个子块的第k个波段向量化，然后对所有向量化的波段根据字典顺序排列成一个矩阵；

步骤3，非凸低秩矩阵分解：经过非凸低秩矩阵分解处理，将步骤2得到的含有混合噪声矩阵，线性化分解成一个低秩矩阵，一个稀疏矩阵，一个无结构化矩阵。此过程中，以矩阵秩函数和稀疏矩阵l_2,1范数的线性和为目标函数，以含有混合噪声矩阵的线性化分解为约束条件，然后利用增广拉格朗日方法求解此问题；

步骤4，还原子块：将低秩矩阵的每个子列合成一个矩阵，然后将所有矩阵按光谱维排列成三维子块，并对重叠部分图像像素值均值化处理；

步骤5，还原高光谱图像：将均值化处理的子块排列成和原始高光谱图像相同的三维图像数据，最后得到去噪后的高光谱图像。

对于步骤1，分割高光谱图像的具体实施细节如下步骤：

步骤1-1，输入原始含有混合噪声的高光谱图像，此图像是由几十乃至数百个连续波段图像组成的三维图像数据，记为其中d_h,d_w,d_s分别代表高光谱图像的高度，宽度和波段数；

步骤1-2，初始化子块的尺寸p，扫描步长s，其中1≤s<p。将原始高光谱图像分割成重叠的子块，共计个。

对于步骤2，生成矩阵的具体实施细节如下步骤：

步骤2-1，对于选取的子块将第k个波段R^p×p向量化为(k＝1,…,d_s)，然后对所有向量化的波段根据字典顺序排列成矩阵

步骤2-2，重复步骤2-1，共计次。

对于步骤3，非凸低秩矩阵分解的具体实施细节如下步骤：

步骤3-1，由于在对目标的空间特征成像的同时，对每个空间像元经过色散形成几十个乃至几百个窄波段以进行连续的光谱覆盖，因此干净的高光谱图像在光谱维上表现出低秩的特性。同时考虑到高光谱图像中含有的噪声，例如：高斯噪声，脉冲噪声，条纹，离群点等等，因此，本发明考虑以下的高光谱图像去噪模型：

D_i＝L_i+S_i+N_i

其中分别为低秩矩阵、稀疏矩阵和无结构化矩阵；

步骤3-2，传统的鲁棒性主成分分析可以写成：

s.t.D_i＝L_i+S_i+N_i

其中，s.t.表示满足条件，min表示最小化函数，||L_i||_*表示矩阵L_i的核范数，||S_i||₁表示矩阵S_i的1范数，λ是正则化参数。然而由于核范数平等地极小化所有的奇异值，导致不能很好地近似秩函数。本发明考虑如下非凸矩阵秩近似：

其中，||L_i||_γ表示低秩矩阵L_i的非凸矩阵秩近似，∑代表求和运算，σ_j(L_i)表述矩阵L_i的第j个奇异值，γ是一个大于0的常数，代表非凸稀疏罚函数，本发明中采用以下两个非凸稀疏罚函数：

同时考虑到高光谱图像中多含有纵向条纹，因此采用稀疏矩阵S_i的2,1范数即以去除条带噪声，其中(S_i)_mn表示变量S_i的第m行第n列元素，表示开平方；

步骤3-3，采用非凸低秩矩阵分解模型：

s.t.D_i＝L_i+S_i+N_i

利用增广拉格朗日方法求解上述模型。其中此模型的增广拉格朗日函数是：

其中Λ_i是拉格朗日乘子，ρ是罚参数，||·||_F是矩阵的F范数，<·>表示矩阵的内积。然后依次迭代更新变量(L_i,S_i,N_i,Λ_i)，第k+1次迭代格式如下：

步骤3-4，关于变量L_i子问题求解，

其中由于||L_i||_γ是非凸的。

容易证明是连续的、凹的、光滑的、可微的并且在区间[0,+∞)是非减的。记m＝min{p²,d_s}，然后对||L_i||_γ在点处一阶泰勒展开，那么子问题可以写成：

其中是非凸稀疏罚函数的导数。然后得到此子问题的封闭解：

其中T_i^k＝U∑V^T是矩阵T_i^k的奇异值分解，U,V分别表示矩阵T_i^k的左酉矩阵和右酉矩阵，V^T表示右酉矩阵的转置，且max表示最大化，diag表示将向量变成一个对角矩阵；

步骤3-5，关于变量S_i子问题求解，

其中此子问题的封闭解的第j列是：

步骤3-6，关于变量N_i子问题求解，

此子问题是一个最小二乘问题，其封闭式解是：

步骤3-7，更新拉格朗日乘子，

步骤3-8，重复步骤3-4至步骤3-7，共计次。

对于步骤4，还原子块的具体实施细节如下步骤：

步骤4-1，对步骤3生成的低秩矩阵选取矩阵L_i的第a列，合成一个p×p矩阵，然后将所有低秩矩阵按光谱维排列成三维子块

步骤4-2，重复步骤4-1，共计次；

步骤4-3，计算重叠部分图像像素值，本发明均值化处理重叠部分图像像素值。

对于步骤5，还原图像的具体实施细节如下步骤：

步骤5-1，将均值化处理的子块排列成和原始高光谱图像相同的三维d_h×d_w×d_s图像数据，最后得到去噪后的高光谱图像。

实施例

本实施的实验硬件环境是：Intel-Core4i767003.4GHz，16G内存，显卡NVIDIAGeForce GTX950。软件环境是MATLAB 2015b。测试图像是来源于网络上公开的高光谱图像数据集：Hyperspectral Digital Imagery Collection Experiment(HYDICE)urban dataset,the Earth Observing-1(EO-1)Hyperion Australia dataset,and the Airborne Visible/InfraredImaging Spectrometer(AVIRIS)Indian Pines dataset。为了验证方法的有效性，本发明选取了近几年提出的先进去噪方法：视频块匹配三维过滤方法(Video Block Matching 3-D filtering，VBM3D^[1])，低秩矩阵恢复方法(Low-Rank Matrix Recovery，LRMR^[2])，噪声自适应的迭代低秩矩阵近似方法(Noise-Adjusted Iterative Low-Rank Matrix Approximation，NAILRMA^[3])，基于加权p-范数的低秩矩阵近似方法(Weighted Schattenp-normIterativeLow-Rank Matrix Approximation，WSN-LRMA^[4])。本发明的方法：非凸低秩矩阵分解方法(Nonconvex Low-Rank Matrix Approximation，NonLRMA)。

实例一：本实例选择遥感领域典型数据集—HYDICEurban数据集，它是一种被高斯噪声、脉冲噪声、条纹、大气、水文等污染的影像数据，是评价高光谱去噪方法的理想数据。整幅图像的尺寸是307×307，共包含210个波段。容易看出本发明中的方法和其他比较方法都可以去除轻度噪声污染波段中的噪声，同时还保存高光谱图像的局部细节信息(如图2所示，其中图2中(a)表示原始含有噪声图像，(b)-(f)分别是由算法VBM3D，LRMR，NAILRMA，WSN-LRMA，NonLRMA去噪之后的图像，剩余的图像都是如此操作)；对于中度噪声污染的波段，部分去噪算法不能有效地去除噪声，相反，本方法不仅去除各种噪声，也很好地保存图像的局部细节信息，如图3所示；对于重度噪声污染的波段，除了本方法可以在去除噪声，同时可以恢复主要细节信息，如图4所示。

实例二：本实例选择遥感领域典型数据集—EO-1Hyperion Australia数据集，此数据集是于2012年12月4日拍摄，原始图像尺寸是256×3858，共包含242个波段。由于篇幅限制，本方法选取其中的200×400的子图像块。此数据集中多含有纵向条纹，如图5-图8所示。容易看出只有本方法可以有效地去除条带噪声，保存图像局部细节信息。