CN109102477B - 一种基于非凸低秩稀疏约束的高光谱遥感图像恢复方法 - Google Patents

Abstract

基于非凸低秩稀疏约束的高光谱遥感图像恢复方法，属于遥感图像处理中的高光谱遥感图像处理领域，为解决现有的高光谱遥感图像恢复技术不能有效去除噪声的问题，提高图像恢复质量，包括以下步骤：输入高光谱遥感图像；初始化权重系数矩阵、迭代次数和收敛阈值,初始化子图像尺寸和扫描步长，划分子块；建立图像恢复模型；引入辅助变量和正则项系数，利用最大值最小值算法迭代求解；判断恢复结果是否满足收敛条件；达到迭代次数获得符合要求的高光谱恢复图像，否则返回相应步骤继续迭代操作；计算权重系数矩阵，给每个子块分配合适的权重；还原高光谱遥感图像，得到最终恢复的高光谱遥感图像。去噪效果明显，并且保留图像细节。

Description

一种基于非凸低秩稀疏约束的高光谱遥感图像恢复方法

技术领域

本发明属于遥感图像处理中的高光谱图像处理领域，具体涉及一种基于非凸低秩稀疏约束模型的高光谱遥感图像恢复方法。

发明背景

高光谱遥感成像技术结合光谱分析和光学成像技术，探测目标的二维几何空间及一维光谱信息，获取高分辨率的连续、窄波段的图像数据。目前，高光谱成像技术发展迅速，因其蕴含丰富的地物光谱信息而被广泛应用于环境研究、地质勘探和军事监视等领域。但是由于传感器的物理缺陷、光子效应、传输损失和校准误差等原因，实际中获得的高光谱遥感图像常常会出现各种类型的退化，比如：噪音污染、条纹影响和数据缺失等。这不仅严重降低了图像的质量，还在很大程度上影响了分类、目标检测和解混等高光谱遥感图像的进一步应用。因此，高光谱遥感图像的恢复是一个重要的亟待解决的问题。

目前，多种不同高光谱遥感图像的恢复方法被应用于高光谱遥感图像的预处理工作中。近几年有学者提出利用高光谱遥感图像数据的相关性构建稀疏和低秩约束优化方程恢复高光谱遥感图像的方法。比如利用利用高光谱遥感图像低秩特性的低秩矩阵恢复(Low-Rank Matrix Recovery，LRMR)方法、噪声自适应迭代低秩矩阵分解(Noise-AdjustedIterative Low-Rank Matrix Approximation，NAILRMA)以及利用高光谱遥感图像稀疏特性的全变分正则化低秩矩阵分解(Total-Variation-Regularized Low-Rank MatrixApproximation，LRTV)方法等。但是此类算法多以低秩约束或稀疏约束的凸松弛l₁范数进行约束优化模型的构建，但是l₁范数由于收缩效应有时存在估计偏差，影响图像恢复质量；此外，当图像数据矩阵维度增加时，收敛速度常不能达到要求。

发明内容

本发明的目的旨在针对上述已有技术的不足，提出一种基于非凸低秩和非凸稀疏约束的高光谱遥感图像恢复方法，在去除高光谱遥感图像中噪声的同时保留图像细节，以提高高光谱遥感图像的恢复质量，满足后续应用要求。

为实现上述目的，本发明所采取的技术方案是：

一种基于非凸低秩稀疏约束的高光谱遥感图像恢复方法，所述方法包含以下步骤：

(1)输入含有混合噪声的高光谱遥感图像Y，Y的大小为M×N×B，其中M和N分别代表各波段高光谱遥感图像的长和宽，B代表高光谱遥感图像波段数；

定义恢复后的图像X^*，初始化变量X⁽⁰⁾＝Y，X^*是X⁽⁰⁾经过一定迭代次数后的恢复图像，初始化X^*的权重系数矩阵Weight，初始化迭代次数和收敛阈值；

定义子块的尺寸p和扫描步长s，并将高光谱遥感图像Y按照预设的尺寸p和扫描步长s划分为若干个部分重叠的子块；

(2)对于选取的第i个子块Y_i，将受噪声污染的高光谱遥感图像分解为

其中

为没有噪声污染的干净的高光谱遥感图像，J_i为包括高斯噪声、条带噪声、脉冲噪声的混合噪声；根据

求解逆问题并引入非凸低秩约束正则项LR(X_i)和非凸稀疏约束正则项S(X_i)，构建高光谱遥感图像的每个子块的恢复模型：

式(2)是利用高光谱遥感图像的低秩性和稀疏性的先验信息建立的图像恢复模型，其中，||·||表示矩阵二范数，α₁和α₂是用于平衡数据一致性项||X_i-Y_i||²和正则项LR(X_i)和S(X_i)的两个正则项系数；X_i表示第i个子块的当前迭代次数的恢复图像；

(3)分别在LR(X_i)和S(X_i)中引入辅助变量R、T_g，利用交替迭代最小化算法分别迭代求解恢复图像X^*和辅助变量R、T_g，过程如下：

首先利用含有辅助变量R、T_g的X^*的二次函数最大化惩罚项，构造恢复模型的代价函数，然后利用迭代交替最小化方法最小化代价函数，求解得到X^*、R、T_g；

(4)判断当前恢复图像是否满足收敛条件，若满足，进入步骤(5)，否则进入步骤(3)，继续利用交替迭代最小化算法迭代求解恢复图像X^*和辅助变量R、T_g；

(5)判断是否满足迭代次数，若满足，则获得最终恢复的高光谱遥感图像子块

并执行步骤(6)，否则返回步骤(3)继续用交替迭代最小化算法迭代求解恢复图像

和辅助变量R、T_g；

(6)计算X^*的权重系数矩阵Weight；

(7)将恢复好的子块

图像按照原位置排列，乘上步骤(6)中得到的权重系数矩阵Weight，还原高光谱遥感图像，得到最终恢复的高光谱遥感图像X^*。

本发明的有益效果是：

本发明的技术要点为：输入高光谱遥感图像；初始化权重系数矩阵、迭代次数和收敛阈值,初始化子图像尺寸和扫描步长，划分子块；建立图像恢复模型；引入辅助变量和正则项系数，利用最大值最小值算法迭代求解；判断恢复结果是否满足收敛条件；达到迭代次数获得符合要求的高光谱恢复图像，否则返回相应步骤继续迭代操作；计算权重系数矩阵，给每个子块分配合适的权重；还原高光谱遥感图像，得到最终恢复的高光谱遥感图像。本发明解决了现有的高光谱遥感图像恢复技术不能有效去除噪声的问题，提高了图像恢复质量。

利用本发明提出的基于非凸低秩稀疏约束的方法恢复高光谱遥感图像，更加充分地挖掘了高光谱遥感图像潜在的低秩性和稀疏性，图像的恢复质量比较高，去噪效果明显，并且有效地保留了图像的边缘细节，能够更好的满足高光谱遥感图像分类、目标检测等后续应用的需求。实验结果表明，与低秩矩阵恢复方法(LRMR)，全变分正则化低秩矩阵分解方法(LRTV)，噪声自适应迭代低秩矩阵逼近方法(NAILRMA)相比，本发明方法能在相同的噪声水平下获得更高的图像恢复质量，去噪效果明显，并且保留图像细节，有效实现了高光谱遥感图像恢复的目的。

附图说明

图1为本发明方法流程图(基于非凸低秩稀疏约束的高光谱遥感图像恢复方法流程图)；

图2为EO-1Hyperion Australia数据集恢复效果图，第52波段恢复效果。

图3为HYDICE Washington DC Mall数据集恢复效果图，第1波段恢复效果。

具体实施方法

下面结合附图和实施实例对本发明进行详细说明。

如图1所示，本发明的具体实施步骤如下：

一种基于非凸低秩稀疏约束的高光谱遥感图像恢复方法，所述方法包含以下步骤：

(1)输入含有混合噪声的高光谱遥感图像Y，Y的大小为M×N×B，其中M和N分别代表各波段高光谱遥感图像的长和宽，B代表高光谱遥感图像波段数；

定义恢复后的图像X^*，初始化变量X⁽⁰⁾＝Y，X^*是X⁽⁰⁾经过一定迭代次数后的恢复图像，初始化X^*的权重系数矩阵Weight，初始化迭代次数和收敛阈值；

定义子块的尺寸p和扫描步长s，并将高光谱遥感图像Y按照预设的尺寸p和扫描步长s划分为若干个部分重叠的子块；

(2)对于选取的第i个子块Y_i，将受噪声污染的高光谱遥感图像分解为

其中

为没有噪声污染的干净的高光谱遥感图像，J_i为包括高斯噪声、条带噪声、脉冲噪声的混合噪声；根据

求解逆问题并引入非凸低秩约束正则项LR(X_i)和非凸稀疏约束正则项S(X_i)，构建高光谱遥感图像的每个子块的恢复模型：

式(2)是利用高光谱遥感图像的低秩性和稀疏性的先验信息建立的图像恢复模型，其中，||·||表示矩阵二范数，α₁和α₂是用于平衡数据一致性项||X_i-Y_i||²和正则项LR(X_i)和S(X_i)的两个正则项系数；X_i表示第i个子块的当前迭代次数的恢复图像；

(3)分别在LR(X_i)和S(X_i)中引入辅助变量R、T_g，利用交替迭代最小化算法分别迭代求解恢复图像X^*和辅助变量R、T_g，过程如下：

首先利用含有辅助变量R、T_g的X^*的二次函数最大化惩罚项，构造恢复模型的代价函数，然后利用迭代交替最小化方法最小化代价函数，求解得到X^*、R、T_g；

(4)判断当前恢复图像是否满足收敛条件，若满足，进入步骤(5)，否则进入步骤(3)，继续利用交替迭代最小化算法迭代求解恢复图像X^*和辅助变量R、T_g；

(5)判断是否满足迭代次数，若满足，则获得最终恢复的高光谱遥感图像子块

并执行步骤(6)，否则返回步骤(3)继续用交替迭代最小化算法迭代求解恢复图像

和辅助变量R、T_g；

(6)计算X^*的权重系数矩阵Weight；

(7)将恢复好的子块

图像按照原位置排列，乘上步骤(6)中得到的权重系数矩阵Weight，还原高光谱遥感图像，得到最终恢复的高光谱遥感图像X^*。

在步骤一中，初始化X^*的权重系数矩阵Weight为大小是M×N×B的零矩阵；定义子块的尺寸为p×p×B，扫描步长为s，其中1≤s≤p；按照预设的子块尺寸和扫描步长把高光谱遥感图像划分为若干个部分重叠的子块，共计

个；p×p×B表示子块的长度×宽度×波段数。

步骤(2)中包含以下操作步骤：

步骤2-1，选取第i个子块Y_i，，i＝1,2,…,

分解受噪声污染的高光谱遥感图像，建立如下退化模型：

其中，Y_i表示含有混合噪声的高光谱遥感图像，

表示恢复后没有噪声影响的高光谱遥感图像，J_i表示包含高斯噪声、脉冲噪声、条带噪声的混合噪声，Y_i、

J_i的大小均为p×p×B；

步骤2-2，基于秩函数的非凸近似函数的Schatten-p范数给出LR(X_i)的表达式，具体形式如下：

其中，0<p₁<1，σ_j是X_i的第j个奇异值；

S(X_i)为X_i的非凸稀疏约束正则项，表示为X_i的TV表达式的l_p2范数，0<p₂<1，具体形式如下：

其中，TV(X_i)为X_i的全变分(Total Variation，TV)表达式，具体形式如下：

其中，D_x，D_y，D_z分别表示X_i在x，y，z轴方向上的差分算子，x，y，z轴所在坐标系以图像X_i的第一个波段的第一个像素点为坐标原点，x和y表示高光谱遥感图像的空间轴，z轴表示光谱轴。

步骤(3)中包含以下操作步骤：

步骤3-1，采用最大值最小值的方法来求解式(2)，引入两个辅助变量R和T_g，利用二次函数最大化惩罚项，得到下式：

其中，||·||_F表示弗罗贝尼乌斯范数；0<p₁<1，T_g中g＝1,2,3时分别表示在x，y，z轴方向上T的值，β₁、β₂分别是LR(X_i)和S(X_i)的收敛系数；

将式(6)、(7)代入式(2)，得到需要求解的优化问题如下：

式中，C为恢复优化问题的代价函数，当C取到最小值时，即得到变量(X_i,R,{T_g})的最优解，C具体表示为：

步骤3-2，求解上述优化问题，依次迭代更新变量(X_i,R,T_g)，第n+1次迭代格式如下：

其中，α₁,α₂为固定大小的系数，β₁,β₂为随着迭代不断增大的系数，T_g,n+1、R_n+1、X_i,n+1分别表示变量T_g经过n+1次迭代后的值、变量R_n经过n+1次迭代后的值和变量X_i,n经过n+1次迭代后的值；

步骤3-3，求解辅助变量T_g,n+1：

首先求中间变量P：

T_g,n+1可以通过下式求解：

步骤3-4，求解辅助变量R_n+1：

首先对X_i,n奇异值分解：

X_i,n＝UΣV^T(15)

其中，U是奇异值分解的右奇异矩阵，V是左奇异矩阵，Σ是X_i,n的奇异值构成的对角矩阵，V^T表示左奇异矩阵V的转置；

R_n+1可以通过下式求解：

其中，(·)₊表示取矩阵大于0的部分；

步骤3-5，利用共轭梯度法来求解变量X_i经过n+1次迭代后的值X_i,n+1：

式(13)的欧拉-拉格朗日方程：

其中，I表示单位矩阵(恒等运算符)；

式(17)的右边表达式的值在每次迭代的过程中都是常数，式(17)可视为一个线性方程组，可用共轭梯度法进行求解X_i,n+1。

步骤(4)中，所述收敛条件为：

其中，Threshold为收敛阈值，Cⁿ⁺¹和Cⁿ分别代表代价函数C经过n+1次和n次迭代求解后的值。

在步骤(6)中，X^*的权重系数矩阵Weight的求解过程：

高光谱遥感图像上像素点每参与迭代求解一次，即在权重系数矩阵Weight相应位置加1，迭代完成后即得到最终的权重系数矩阵。

本发明将通过以下实验进行有效性的说明

1.实验环境：

实验室用台式机参数：CPU为Inter(R)Core(TM)i5-4590，主频为3.30GHz，内存为8G，操作系统为Win10 64位系统，实验软件为MatlabR2014a。

2.实验结果及结果分析

本发明实验所用的图像数据为网络上公开的高光谱遥感图像数据集：Hyperspectral Digital Imagery Collection Experiment(HYDICE)Washington DCMall，the Earth Observing-1(EO-1)Hyperion Australia dataset。为了验证方法的有效性，本发明与近几年提出的高光谱遥感图像恢复方法进行对比：低秩矩阵恢复方法(Low-Rank Matrix Recovery,LRMR)，全变分正则化低秩矩阵分解方法(Total-Variation-Regularized Low-Rank Matrix Approximation，LRTV)，噪声自适应迭代低秩矩阵逼近方法(Noise-Adjusted Iterative Low-Rank Matrix Approximation，NAILRMA)。本发明的方法：非凸低秩和非凸稀疏高光谱恢复方法(Nonconvex Low-Rank and Nonconvex Total-Variation-Regularized Hyperspectral Image Restoration，NonLRTV)。

实例一：本实例选择遥感领域典型数据集——EO-1Hyperion Australia数据集，原始图像的尺寸是3858×256×242，这里因为篇幅限制，且去除可见光近红外和短波红外范围内的重叠波段后，选取尺寸为400×200×150，将其每个像素值减去最小像素值后，再除以最大与最小像素值之差，得到归一化到[0，1]后的数据。使用本发明提出的方法进行图像恢复处理，并用上述提到的LRMR、LRTV、NAILRMA进行对比实验。实验结果如图2所示，其中图2(a)为第52波段原始的高光谱遥感图像，图2(b)为LRMR方法效果，图2(c)为LRTV方法效果，图2(d)为NAILRMA方法效果，图2(e)为本发明方法效果。可以看出，图2(a)中含有大量噪声，而图2(b)、图2(c)和图2(d)都有明显的噪声，虽然去除部分噪声，但是去除效果有限。图2(e)中本发明的去噪效果明显，不仅很好地去除了噪声，而且很大程度上保留了丰富的细节。

实例二：本实例选择遥感领域典型数据集——Hyperspectral Digital ImageryCollection Experiment(HYDICE)Washington DC Mall数据集，原始图像的尺寸是1208×307×191，因为篇幅限制，选取图像尺寸为256×256×11。与实例一相同，对原始图像数据进行归一化处理，人为地对所有波段随机添加均值为25dB的高斯噪声。然后用本发明提出的方法进行恢复处理，并与其他方法进行对比。实验结果如图3所示，其中图3(a)为其中一张清晰的高光谱遥感图像数据，图3(b)为添加高斯噪声和条带噪声后的图像，图3(c)为LRMR方法效果(SNR＝10.8786dB，MSSIM＝0.66377)，图3(d)为LRTV方法效果(SNR＝10.6603dB，MSSIM＝0.74259)，图3(e)为NAILRMA方法效果(SNR＝15.5521dB，MSSIM＝0.83878)，图3(f)为本发明方法效果(SNR＝17.087dB，MSSIM＝0.89488)。可以看出，图3(c)和图3(e)仍有明显噪声，图3(d)虽然没有明显噪声，但是过于模糊，丢失了许多细节。图3(f)中本发明的去噪效果明显，能够在有效去除大部分噪声的同时，保留绝大部分的细节，恢复效果最好，并且评价指标SNR、MSSIM的值较前三种方法均有较大程度的提高。

在相同的噪声水平下，利用本发明中提出的方法恢复出的高光谱遥感图像具有较高的信噪比，并且与清晰的未加噪声的高光谱遥感图像的平均结构相似度最高。可见，相比于其他三种现已提出的方法，本发明所提方法可以获得更好的恢复结果，与为加噪声的原数据图像最为接近，恢复效果最好。

Claims (3)

1.一种基于非凸低秩稀疏约束的高光谱遥感图像恢复方法，其特征在于，所述方法包含以下步骤：

(1)输入含有混合噪声的高光谱遥感图像Y，Y的大小为M×N×B，其中M和N分别代表各波段高光谱遥感图像的长和宽，B代表高光谱遥感图像波段数；

定义恢复后的图像X^*，初始化变量X⁽⁰⁾＝Y，X^*是X⁽⁰⁾经过一定迭代次数后的恢复图像，初始化X^*的权重系数矩阵Weight，初始化迭代次数和收敛阈值；

定义子块的尺寸p和扫描步长s，并将高光谱遥感图像Y按照预设的尺寸p和扫描步长s划分为若干个部分重叠的子块；

(2)对于选取的第i个子块Y_i，将受噪声污染的高光谱遥感图像分解为

其中

为没有噪声污染的干净的高光谱遥感图像，J_i为包括高斯噪声、条带噪声、脉冲噪声的混合噪声；根据

求解逆问题并引入非凸低秩约束正则项LR(X_i)和非凸稀疏约束正则项S(X_i)，构建高光谱遥感图像的每个子块的恢复模型：

式(2)是利用高光谱遥感图像的低秩性和稀疏性的先验信息建立的图像恢复模型，其中，||·||表示矩阵二范数，α₁和α₂是用于平衡数据一致性项||X_i-Y_i||²和正则项LR(X_i)和S(X_i)的两个正则项系数；X_i表示第i个子块的当前迭代次数的恢复图像；

(3)分别在LR(X_i)和S(X_i)中引入辅助变量R、T_g，利用交替迭代最小化算法分别迭代求解恢复图像X^*和辅助变量R、T_g，过程如下：

首先利用含有辅助变量R、T_g的X^*的二次函数最大化惩罚项，构造恢复模型的代价函数，然后利用迭代交替最小化方法最小化代价函数，求解得到X^*、R、T_g；

(4)判断当前恢复图像是否满足收敛条件，若满足，进入步骤(5)，否则进入步骤(3)，继续利用交替迭代最小化算法迭代求解恢复图像X^*和辅助变量R、T_g；

(5)判断是否满足迭代次数，若满足，则获得最终恢复的高光谱遥感图像子块

并执行步骤(6)，否则返回步骤(3)继续用交替迭代最小化算法迭代求解恢复图像

和辅助变量R、T_g；

(6)计算X^*的权重系数矩阵Weight；

(7)将恢复好的子块

图像按照原位置排列，乘上步骤(6)中得到的权重系数矩阵Weight，还原高光谱遥感图像，得到最终恢复的高光谱遥感图像X^*；

步骤(1)中，初始化X^*的权重系数矩阵Weight为大小是M×N×B的零矩阵；定义子块的尺寸为p×p×B，扫描步长为s，其中1≤s≤p；按照预设的子块尺寸和扫描步长把高光谱遥感图像划分为若干个部分重叠的子块，共计

个；p×p×B表示子块的长度×宽度×波段数；

步骤(2)中包含以下操作步骤：

步骤2-1，选取第i个子块Y_i，

分解受噪声污染的高光谱遥感图像，建立如下退化模型：

其中，Y_i表示含有混合噪声的高光谱遥感图像，

表示恢复后没有噪声影响的高光谱遥感图像，J_i表示包含高斯噪声、脉冲噪声、条带噪声的混合噪声，Y_i、

J_i的大小均为p×p×B；

步骤2-2，基于秩函数的非凸近似函数的Schatten-p范数给出LR(X_i)的表达式，具体形式如下：

其中，0<p₁<1，σ_j是X_i的第j个奇异值；

S(X_i)为X_i的非凸稀疏约束正则项，表示为X_i的TV表达式的

范数，0<p₂<1，具体形式如下：

其中，TV(X_i)为X_i的全变分(Total Variation，TV)表达式，具体形式如下：

其中，D_x，D_y，D_z分别表示X_i在x，y，z轴方向上的差分算子，x，y，z轴所在坐标系以图像X_i的第一个波段的第一个像素点为坐标原点，x和y表示高光谱遥感图像的空间轴，z轴表示光谱轴；

步骤(3)中包含以下操作步骤：

步骤3-1，采用最大值最小值的方法来求解式(2)，引入两个辅助变量R和T_g，利用二次函数最大化惩罚项，得到下式：

其中，||·||_F表示弗罗贝尼乌斯范数；0<p₁<1，T_g中g＝1,2,3时分别表示在x，y，z轴方向上T的值，β₁、β₂分别是LR(X_i)和S(X_i)的收敛系数；

将式(6)、(7)代入式(2)，得到需要求解的优化问题如下：

式中，C为恢复优化问题的代价函数，当C取到最小值时，即得到变量(X_i,R,{T_g})的最优解，C具体表示为：

步骤3-2，求解上述优化问题，依次迭代更新变量(X_i,R,T_g)，第n+1次迭代格式如下：

其中，β₁,β₂为随着迭代不断增大的系数，T_g,n+1、R_n+1、X_i,n+1分别表示变量T_g经过n+1次迭代后的值、变量R_n经过n+1次迭代后的值和变量X_i,n经过n+1次迭代后的值；

步骤3-3，求解辅助变量T_g,n+1：

首先求中间变量P：

T_g,n+1通过下式求解：

步骤3-4，求解辅助变量R_n+1：

首先对X_i,n奇异值分解：

X_i,n＝UΣV^T (15)

其中，U是奇异值分解的右奇异矩阵，V是左奇异矩阵，Σ是X_i,n的奇异值构成的对角矩阵，V^T表示左奇异矩阵V的转置；

R_n+1通过下式求解：

其中，(·)₊表示取矩阵大于0的部分；

步骤3-5，利用共轭梯度法来求解变量X_i经过n+1次迭代后的值X_i,n+1：

式(13)的欧拉-拉格朗日方程：

其中，I表示单位矩阵；

式(17)的右边表达式的值在每次迭代的过程中都是常数，式(17)视为一个线性方程组，用共轭梯度法进行求解X_i,n+1。

2.根据权利要求1所述的一种基于非凸低秩稀疏约束的高光谱遥感图像恢复方法，其特征在于，

步骤(4)中，所述收敛条件为：

其中，Threshold为收敛阈值，Cⁿ⁺¹和Cⁿ分别代表代价函数C经过n+1次和n次迭代求解后的值。

3.根据权利要求2所述的一种基于非凸低秩稀疏约束的高光谱遥感图像恢复方法，其特征在于，

在步骤(6)中，X^*的权重系数矩阵Weight的求解过程：

高光谱遥感图像上像素点每参与迭代求解一次，即在权重系数矩阵Weight相应位置加1，迭代完成后即得到最终的权重系数矩阵。

Images