CN114841875A - 一种基于图学习和降噪的高光谱图像解混方法 - Google Patents

Abstract

本发明公开了一种基于图学习和降噪的高光谱图像解混方法，所述方法包括如下步骤：步骤1：从观测高光谱图像中图学习；步骤2：高光谱图像的降噪处理；步骤3：基于降噪后的高光谱图像Y解混；步骤4：终止条件。本发明利用图学习的方法从含有噪声的观测的高光谱图像中学习得到一个空间Laplacian矩阵，通过空间Laplacian矩阵可自适应不同情况的高光谱图像。该Laplacian矩阵保留了实际的高光谱图像的像素间的真实空间结构关系的同时还能够用于去除高光谱图像中存在的噪声，然后再基于去噪后的高光谱图像进行解混处理，这可在一定程度上提升高光谱图像解混的性能，能够自适应不同情况的高光谱图像，更加鲁棒和准确。

Description

一种基于图学习和降噪的高光谱图像解混方法

技术领域

本发明属于遥感技术的研究领域，涉及一种高光谱图像解混方法，具体涉及一种对降噪后的高光谱图像的像素进行解混处理并得到相应地表物质丰度表示的方法。

背景技术

在高光谱遥感中，高光谱图像记录了场景中的地表物质的光谱特征，它的每个像素的光谱包含数百个连续波段丰富的光谱信息。高光谱图像已广泛用于各个领域的许多应用，包括矿产勘探、环境监测、农业和军事应用等。在实践中，由于高光谱传感器的低空间分辨率和地表物质分布的高度复杂性，导致了高光谱图像包含高度的混合像素，其每个像素的光谱是来自多种不同材料的光谱特征的混合。混合像素的存在严重限制了高光谱图像的应用，例如，它会给高光谱图像分析在聚类和分类任务中带来巨大挑战。对高光谱的每个像素进行解混是解决这一问题的有效的方法，其已经成为了高光谱图像分析与进一步开发的必要过程。因而，需要对高光谱图像进行有效的解混处理。

在实际场景的高光谱图像中，由于光照条件、大气效应以及设备自身精度等因素的干扰，导致了实际的高光谱图像中存在大量的干扰信息，这些干扰信息可视作是高光谱图像中的噪声。目前已有的稀疏解混方法都是基于含有噪声的观测的高光谱图像进行解混，但由于观测图像受到了噪声的干扰，这会在很大程度上影响对高光谱图像的混合像素进行解混处理。

发明内容

本发明基于观测的高光谱遥感图像中存在噪声干扰这一动机，提供了一种基于图学习和降噪的高光谱图像解混方法。该方法基于降噪后的高光谱图像进行解混，去除了噪声对解混过程的干扰，可进一步提升高光谱图像的解混性能，可为相关的专业人员对高光谱图像的进一步分析与应用提供一定的理论分析依据。

本发明的目的是通过以下技术方案实现的：

一种基于图学习和降噪的高光谱图像解混方法，包括如下步骤：

步骤1：从观测高光谱图像中图学习

首先将3-D的观测高光谱图像S进行分割与融合，然后将预处理后的图像用于空间Laplacian矩阵的图学习；

步骤2：高光谱图像的降噪处理

基于步骤1图学习过程所学得的空间Laplacian矩阵L，对观测的含噪高光谱图像S按下式进行降噪处理：

Y＝S(I_n+λ_lL)^-1；

其中，I_n表示维数为n的单位阵，λ_l为正则参数，上角标^-1表示矩阵的求逆运算，Y为降噪后的高光谱图像；

步骤3：基于降噪后的高光谱图像Y解混

步骤3.1、解混过程中变量X按下式进行计算：

X^(k+1)＝[Λ₁Λ₂]₊；

其中，[z]₊＝max(z,0)，Λ₁＝(A^TA+μ(W^(k))^-2)^-1，Λ₂＝A^TY+μ(W^(k))^-1(V^(k)-U^(k))，上角标^-2表示对矩阵所有元素进行平方后的求逆运算，上角标^-1表示矩阵的求逆运算，上角标^T表示矩阵的转置运算；

步骤3.2、解混过程内层循环的变量V按下式进行计算：

其中，H＝W^(k)X^(k+1)+U^(k)，H(i,:)表示矩阵H的第i行；

步骤3.3、解混过程中Lagrangian乘子按下式进行更新：

U^(k+1)＝U^(k)+(W^(k))^-1X^(k+1)-V^(k+1)；

步骤3.4、解混过程中促进丰度稀疏性的权重矩阵W按下式计算：

其中，ε为解混过程设定的超参数，X^(k+1)(i,:)表示在解混过程的第k+1次迭代时，估计矩阵X的第i行；

步骤4：终止条件

若解混过程中的迭代次数k＝Q_max，则终止迭代并输出矩阵X^(k+1)；否则转到步骤3继续执行；当依据步骤4停止时，输出的X^(k+1)就是基于降噪后的高光谱图像Y所估计出的丰度矩阵。

相比于现有技术，本发明具有如下优点：

1、本发明与目前广泛用于高光谱遥感图像解混的其他稀疏解混方法相比，所提出的解混方法具有更好的解混性能和很低的计算代价。

2、本发明基于降噪后的高光谱图像进行解混，通过使用观测的高光谱图像利用图学习的方法从含有噪声的观测的高光谱图像中学习得到一个空间Laplacian矩阵，通过学习得到的空间Laplacian矩阵可自适应不同情况的高光谱图像。该Laplacian矩阵保留了实际的高光谱图像的像素间的真实空间结构关系的同时还能够用于去除高光谱图像中存在的噪声，然后再基于去噪后的高光谱图像进行解混处理，这可在一定程度上提升高光谱图像解混的性能，能够自适应不同情况的高光谱图像，更加鲁棒和准确。因此，本发明的方法具有计算代价低、可靠性好以及精度相对较高的特点，适合推广使用。

3、本发明基于国际公开的光谱数据库，使用了被广泛应用于高光谱图像解混任务中的两个合成数据集DC1、DC2和一个真实的铜矿石(Cuprite)数据集进行实验，这些数据集可被用于定量分析高光谱解混方法的性能。基于这些数据集的实验表明，本发明中所提出的解混方法远好于当前其他的稀疏解混方法。

附图说明

图1为采用现有的SUnSAL、SUnSAL-TV、S²WSU、MUA-SLIC及SBGLSU五种稀疏解混方法与本发明中所述的解混方法，在输入信噪比为10dB时，对高光谱图像模拟数据集DC1进行解混时的性能表现图；

图2为采用现有的SUnSAL、SUnSAL-TV、S²WSU、MUA-SLIC及SBGLSU五种稀疏解混方法与本发明中所述的解混方法，在输入信噪比为10dB时，对高光谱图像模拟数据集DC2进行解混时的性能表现图；

图3为采用现有的SUnSAL、SUnSAL-TV、S²WSU、MUA-SLIC及SBGLSU五种稀疏解混方法与本发明中所述的解混方法，对真实的铜矿石(Cuprite)高光谱图像的数据集进行解混时的性能表现图；

图中，Ground Truth为用于衡量解混性能的真实参考丰度图，SUnSAL、SUnSAL-TV、S²WSU、MUA-SLIC以及SBGLSU为五种稀疏解混对比方法所生成的丰度图，GLDWSU为本发明所提出的方法所生成的丰度图。

具体实施方式

下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。

本发明提供了一种基于图学习和降噪的高光谱图像解混方法，所述方法包括如下步骤：

步骤1：从观测高光谱图像中图学习

为了同时减少高光谱图像谱域中的冗余信息和计算成本，首先将3-D的观测高光谱图像S进行分割与融合，然后将预处理后的图像更好的用于空间Laplacian矩阵的图学习，具体步骤如下：

步骤1.1：首先对观测图像S进行分割处理，将图像S的每个波段的图像视为一个子图像，则这些子图像可用集合表示为

其中，s_i表示3-D观测高光谱图像的第i个波段的子图像，d为总的波段数；然后将d个子图像划分为由相邻波段所构成的多个子集

其中，K是子集的数量，即：

其中，N_k是由与第k个波段相邻的多个波段的子图像所构成的索引集合；最后，每个子集合S_i由下式进行融合：

其中，|N_k|表示集合的大小，

包含了每个子集S_k中的平均空间信息。

在预处理过程中，不同的数据集的参数设置如下：

模拟数据集DC1中子集数量K：K＝5。

模拟数据集DC2中子集数量K：K＝5。

真实的铜矿石数据集中子集数量K：K＝20。

步骤1.2：基于融合的图像

来学习空间Laplacian矩阵，在Laplacian矩阵的学习过程中，相关参数及输入量设置如下：

所有数据集在图学习过程中的初始化设置：

其中，k＝1,…,K。

模拟数据集DC1的相关参数设置：高斯滤波器的标准差：σ＝1.2；图学习的正则参数：λ_l＝0.003；超参数：ξ＝0.001；图学习过程的最大迭代次数：T_max＝5。

模拟数据集DC2的相关参数设置：高斯滤波器的标准差：σ＝0.5；图学习的正则参数：λ_l＝0.008；超参数：ξ＝0.01；图学习过程的最大迭代次数：T_max＝5。

铜矿石数据集的相关参数设置：高斯滤波器的标准差：σ＝0.1；图学习的正则参数：λ_l＝0.05；超参数：ξ＝15；图学习过程的最大迭代次数：T_max＝10。

具体步骤如下：

步骤1.2.1、计算对角权重矩阵W_r和W_c按下式计算：

其中，D_r和D_c分别为水平和垂直两个方向上的离散差分算子，G_σ是标准差为σ的高斯滤波器，ξ为图学习过程设定的超参数，“*”为卷积算子，

为基于元素运算的乘积算子，以及除法和绝对值都是基于元素的运算；

步骤1.2.2、计算空间Laplacian矩阵

其中，D_r和D_c分别为水平和垂直两个方向上的离散差分算子，上角标^T为矩阵的转置运算；

步骤1.2.3、在图学习的第t次迭代过程中，变量

按下式进行更新：

其中，I_n表示维数为n的单位阵，上角标^-1表示矩阵求逆运算；

步骤1.2.4、图学习过程的终止条件：

若迭代次数t＝T_max，则终止迭代并按下式输出空间Laplacian矩阵L：

否则空间Laplacian矩阵的图学习转到步骤1.2.1继续执行；

步骤2：高光谱图像的降噪处理

基于步骤1图学习过程所学得的空间Laplacian矩阵L，对观测的含噪高光谱图像S按下式进行降噪处理：

Y＝S(I_n+λ_lL)^-1 (7)；

其中，I_n表示维数为n的单位阵，λ_l为步骤1中所设定的正则参数，上角标^-1表示矩阵的求逆运算，Y为降噪后的高光谱图像。

步骤3：基于降噪后的高光谱图像Y解混

所有降噪后的高光谱图像数据集的解混过程的相关参数及输入量设置如下：

模拟数据集DC1、DC2的相关参数设置：正则参数：λ_u＝0.5；罚参数：μ＝0.01；超参数：ε＝0.0001；解混过程的最大迭代次数设置：Q_max＝100。

真实铜矿石数据集的相关参数设置：正则参数：λ_u＝0.1；罚参数：μ＝0.6；超参数：ε＝0.0001；解混过程的最大迭代次数设置：Q_max＝200。

算法所要估计的丰度矩阵的初始化X⁽⁰⁾设置如下：

[UF,SF]＝svd(A^TA)

X⁽⁰⁾＝IF*A^T*Y

式(8)中的变量说明：“svd”表示奇异值分解运算，“*”为矩阵乘积运算，上角标T表示矩阵的转置运算，A为已知的谱库矩阵，Y为降噪的高光谱图像，“diag(·)”表示构建一个对角阵，或以向量形式返回一个矩阵对角线上的元素，μ为解混过程设定的罚参数。

解混过程的辅助变量V⁽⁰⁾的初适化设置：V⁽⁰⁾＝Y；

解混过程的Lagrangian乘子初始化设置：U⁽⁰⁾为全零矩阵。

具体步骤如下：

步骤3.1、解混过程中变量X按下式进行计算：

X^(k+1)＝[Λ₁Λ₂]₊ (9)；

其中，[z]₊＝max(z,0)，Λ₁＝(A^TA+μ(W^(k))^-2)^-1，Λ₂＝A^TY+μ(W^(k))^-1(V^(k)-U^(k))，上角标^-2表示对矩阵所有元素进行平方后的求逆运算，上角标^-1表示矩阵的求逆运算，上角标T表示矩阵的转置运算；

步骤3.2、解混过程内层循环的变量V按下式进行计算：

其中，H＝W^(k)X^(k+1)+U^(k)，H(i,:)表示矩阵H的第i行；

步骤3.3、解混过程中Lagrangian乘子按下式进行更新：

U^(k+1)＝U^(k)+(W^(k))^-1X^(k+1)-V^(k+1) (11)；

步骤3.4、解混过程中促进丰度稀疏性的权重矩阵W按下式计算：

其中，ε为解混过程设定的超参数，X^(k+1)(i,:)表示在解混过程的第k+1次迭代时，估计矩阵X的第i行，i＝1,…,m。

步骤4：算法的终止条件

若解混过程中的迭代次数k＝Q_max，则终止迭代并输出矩阵X^(k+1)；否则算法转到步骤3继续执行；当算法依据步骤4停止时，算法输出的X^(k+1)就是基于降噪后的高光谱图像Y所估计出的丰度矩阵。

为验证本发明中所提出的解混方法的有效性，我们分别采用SUnSAL、SUnSAL-TV、S²WSU、MUA-SLIC以及SBGLSU五种稀疏解混方法和本发明中所述的解混方法(GLDWSU)，使用了两个合成的模拟数据集DC1、DC2和一个真实的铜矿石(Cuprite)数据集进行解混处理。模拟数据集DC1的高光谱图像包含了224个波段，每个波段的图像由75×75的像素构成，包含25个以网格方式空间排列的正方形，其中每个正方形是一个同质区域，其像素以相同的丰度生成。模拟数据集DC2含有224个波段，每个波段的图像包含了100×100个像素，其通过使用以高斯随机场为中心的狄利克雷分布的相应丰度分数组成。真实的铜矿石(Cuprite)数据集包含了188个波段，每个波段的图像由250×191的像素构成。

我们在实验中，对于数据集DC1选择了第3端元、第4端元和第5端元的丰度图作为代表；对于数据集DC2选择了第3端元、第8端元和第9端元的丰度图作为代表；对于真实的铜矿石数据集选择了明矾石(Alunite)、布丁英石(Buddingtonite)以及玉髓(Chalcedony)这三种矿物作为代表，所提出的方法在每个数据集的解混性能分别如图1、图2和图3所示。从图中可以看出，本发明所提出的解混方法的性能表现都远好于同类型的其他稀疏解混方法，此实验充分证实了所提出方法的有效性。

此外，为了更进一步的说明本发明所提出方法的有效性，我们利用信号的重构误差(signal-reconstruction-error,SRE)来度量所提出解混方法的性能，该指标被广泛的应用于高光谱解混方法性能的度量中。信号重构误差SRE的具体定义方式如下：

其中，X和

分别为真实的丰度图和估计的丰度图，通常情况，SRE值越大说明解混方法的性能越好。

在实验中，我们分别对模拟数据集DC1和DC2加入了信噪比为10dB的高斯噪声，以用于模拟含有噪声的高光谱图像的数据集。我们分别采用SUnSAL、SUnSAL-TV、S²WSU、MUA-SLIC以及SBGLSU五种稀疏解混方法和本发明中所述的解混方法，对加入10dB高斯噪声的模拟数据集DC1和DC2进行了解混处理，所提出的方法的SRE值和计算时间如表1和表2所示。从表1和表2可以看出，本发明所提出的解混方法的性能表现都远好于同类型的其他稀疏解混方法，同时还具有很低的计算代价。此实验更进一步的证实了本发明所提出方法的有效性。

表1各种方法在模拟数据集DC1、DC2上输入信噪比为10dB时的解混SRE值

Data	SREValue	SUnSAL	SUnSAL-TV	S<sup>2</sup>WSU	MUA-SLIC	SBGLSU	GLDWSU
								DC1	SRE(dB)	0.07	4.67	3.65	7.89	7.62	16.81
DC2	SRE(dB)	1.37	2.10	3.46	9.10	8.73	14.63

表2各种方法在所有高光谱数据集上的计算时间

Data

SUnSAL

SUnSAL-TV

S<sup>2</sup>WSU

MUA-SLIC

SBGLSU

GLDWSU

DC1

2.46s

57.53s

13.56s

1.94s

26.32s

2.78s

DC2

3.93s

88.92s

21.98s

2.51s

58.93s

4.60s

RealCupritedata

52.77s

693.44s

166.71s

47.15s

777.32s

86.56s

注：由于铜矿石数据是真实的高光谱图像数据，不存在真实的丰度值作为参考，因此没有列出关于铜矿石数据集的性能值的表格。

Claims (5)

1.一种基于图学习和降噪的高光谱图像解混方法，其特征在于所述方法包括如下步骤：

步骤1：从观测高光谱图像中图学习

首先将3-D的观测高光谱图像S进行分割与融合，然后将预处理后的图像用于空间Laplacian矩阵的图学习；

步骤2：高光谱图像的降噪处理

基于步骤1图学习过程所学得的空间Laplacian矩阵L，对观测的含噪高光谱图像S按下式进行降噪处理：

Y＝S(I_n+λ_lL)^-1；

其中，I_n表示维数为n的单位阵，λ_l为正则参数，上角标-1表示矩阵的求逆运算，Y为降噪后的高光谱图像；

步骤3：基于降噪后的高光谱图像Y解混

步骤3.1、解混过程中变量X按下式进行计算：

X^(k+1)＝[Λ₁Λ₂]₊；

其中，[z]₊＝max(z,0)，Λ₁＝(A^TA+μ(W^(k))^-2)^-1，Λ₂＝A^TY+μ(W^(k))^-1(V^(k)-U^(k))，上角标^-2表示对矩阵所有元素进行平方后的求逆运算，上角标^-1表示矩阵的求逆运算，上角标^T表示矩阵的转置运算；

步骤3.2、解混过程内层循环的变量V按下式进行计算：

其中，H＝W^(k)X^(k+1)+U^(k)，H(i,:)表示矩阵H的第i行；

步骤3.3、解混过程中Lagrangian乘子按下式进行更新：

U^(k+1)＝U^(k)+(W^(k))^-1X^(k+1)-V^(k+1)；

步骤3.4、解混过程中促进丰度稀疏性的权重矩阵W按下式计算：

其中，ε为解混过程设定的超参数，X^(k+1)(i,:)表示在解混过程的第k+1次迭代时，估计矩阵X的第i行；

步骤4：终止条件

若解混过程中的迭代次数k＝Q_max，则终止迭代并输出矩阵X^(k+1)；否则转到步骤3继续执行；当依据步骤4停止时，输出的X^(k+1)就是基于降噪后的高光谱图像Y所估计出的丰度矩阵。

2.根据权利要求1所述的基于图学习和降噪的高光谱图像解混方法，其特征在于所述步骤1的具体步骤如下：

步骤1.1：首先对观测图像S进行分割处理，将图像S的每个波段的图像视为一个子图像，则这些子图像可用集合表示为

其中，s_i表示3-D观测高光谱图像的第i个波段的子图像，d为总的波段数；然后将d个子图像划分为由相邻波段所构成的多个子集

其中，K是子集的数量，即：

其中，N_k是由与第k个波段相邻的多个波段的子图像所构成的索引集合；最后，每个子集合S_i由下式进行融合：

其中，|N_k|表示集合的大小，

包含了每个子集S_k中的平均空间信息；

步骤1.2：基于融合的图像

来学习空间Laplacian矩阵。

3.根据权利要求2所述的基于图学习和降噪的高光谱图像解混方法，其特征在于所述步骤1.2的具体步骤如下：

步骤1.2.1、计算对角权重矩阵W_r和W_c按下式计算：

其中，D_r和D_c分别为水平和垂直两个方向上的离散差分算子，G_σ是标准差为σ的高斯滤波器，ξ为图学习过程设定的超参数，“*”为卷积算子，

为基于元素运算的乘积算子；

步骤1.2.2、计算空间Laplacian矩阵

其中，D_r和D_c分别为水平和垂直两个方向上的离散差分算子，上角标^T为矩阵的转置运算；

步骤1.2.3、在图学习的第t次迭代过程中，变量

按下式进行更新：

其中，I_n表示维数为n的单位阵，上角标^-1表示矩阵求逆运算；

步骤1.2.4、图学习过程的终止条件：

若迭代次数t＝T_max，则终止迭代并按下式输出空间Laplacian矩阵L：

否则空间Laplacian矩阵的图学习转到步骤1.2.1继续执行。

4.根据权利要求2或3所述的基于图学习和降噪的高光谱图像解混方法，其特征在于所述步骤1.2中，在Laplacian矩阵的学习过程中，所有数据集在图学习过程中的初始化设置：

其中，k＝1,…,K。

5.根据权利要求1所述的基于图学习和降噪的高光谱图像解混方法，其特征在于所述步骤3中，丰度矩阵的初始化X⁽⁰⁾设置如下：

其中，“svd”表示奇异值分解运算，“*”为矩阵乘积运算，上角标^T表示矩阵的转置运算，A为已知的谱库矩阵，Y为降噪的高光谱图像，“diag(·)”表示构建一个对角阵，或以向量形式返回一个矩阵对角线上的元素，μ为解混过程设定的罚参数；

解混过程的辅助变量V⁽⁰⁾的初适化设置：V⁽⁰⁾＝Y；

解混过程的Lagrangian乘子初始化设置：U⁽⁰⁾为全零矩阵。

Images