一种基于空间结构信息约束的高光谱图像解混方法

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


背景技术：

2.高光谱遥感是一种重要的综合对地观测技术。高光谱遥感图像中包含了丰富的对地观测信息，即：它同时包含了反映特征辐射的光谱信心和反映特征二维空间的图像信息。高光谱图像的每个像素是视场中物质在多个连续波段的电磁波信号下的综合反映，因而高光谱图像还提供了对地观测的特征属性探测信息。通过对高光谱图像的解混处理，可提供对地观测物质构成的探测方法，而且结合地理信息还能够提供更为丰富的环境与灾害的监测方法等；通过对此问题的研究，还能为高光谱遥感图像提供更为有效的分析方法，进而能够更好地应用于高光谱遥感图像的分类与检测等任务。然而，由于高光谱图像的空间分辨率较低，这导致了高光谱图像的每个像素中同时记录了多种物质的光谱信息，这严重限制了对高光谱遥感图像的进一步分析与应用。因而，需要对高光谱遥感图像进行有效的解混处理。
3.为了实现对高光谱图像进行有效的解混，通过在解混过程中引入高光谱图像的像素间的空间结构信息是一种常见的方法。空间信息的引入需要依据合理的先验知识，这些先验信息需要满足实际的高光谱图像的结构特点，以及更符合实际的物理意义。不同的空间信息引入所使用的先验假设和方法也有所不同，在目前现有的一些同类型的解混方法中，通过假设相邻像素间的谱信息是相似的先验知识，然后手动构建能够刻画相邻像素间空间结构关系的空间矩阵，该矩阵强制性的刻画了相邻像素间的结构关系。然而，手动构建的空间矩阵所刻画的空间结构信息往往并不符合高光谱图像中像素间的真实空间结构。通过引入合理的空间结构信息，可在一定程度上提升高光谱图像解混的性能。因此，这需要一种更为合理的空间信息被用于高光谱图像的解混任务中。


技术实现要素：

4.本发明基于高光谱遥感图像相邻像素间存在着一定的空间结构信息，提供了一种基于空间结构信息约束的高光谱图像解混方法。本发明在高光谱图像的像素间的空间信息刻画上引入了一种新的方法，该方法基于观测的高光谱图像可利用图学习的方式学习得到能够刻画像素间空间结构信息的空间矩阵，其能够自适应不同情况的高光谱图像，更加鲁棒和准确，可为相关的专业人员对高光谱图像的进一步分析与应用提供一定的理论分析依据。
5.本发明的目的是通过以下技术方案实现的：
6.一种基于空间结构信息约束的高光谱图像解混方法，包括如下步骤：
7.步骤1：高光谱图像的预处理
8.首先将3-d的高光谱图像y分割为由相邻多个波段子图像所构成的图像集合
其中，k为图像集合的数量，每个图像集合包含了原高光谱图像的相邻的多个波段的子图像，即：其中，nk为与第k个波段相邻的子图像所构成的集合，yj为原高光谱图像y的第j个波段的子图像；
9.然后，每个图像集合yi由下式进行融合处理：
[0010][0011]
其中，|nk|表示集合的大小，包含了每个图像集合yk中的平均空间信息；
[0012]
步骤2：空间矩阵的图学习
[0013]
步骤2.1、计算对角权重矩阵w
x
和wy，其对角元素按下式计算：
[0014][0015][0016]
其中，d
x
和dy分别为水平和垂直两个方向上的离散差分算子，g
σ
是标准差为σ的高斯滤波器，ξ为图学习过程设定的超参数，“*”为卷积算子，为基于元素运算的乘积算子；
[0017]
步骤2.2、计算图laplacian矩阵
[0018][0019]
其中，d
x
和dy分别为水平和垂直两个方向上的离散差分算子，上角标
t
为矩阵的转置运算；
[0020]
步骤2.3、计算空间矩阵h
x
和h
y`
按下式计算：
[0021][0022][0023]
其中，w
1/2
表示对矩阵w的所有元素取平方根运算；
[0024]
步骤2.4、在图学习的第t次迭代过程中，变量按下式进行更新：
[0025][0026]
其中，in表示维数为n的单位阵，上角标-1
表示矩阵求逆运算；
[0027]
步骤2.5、图学习过程的终止条件：
[0028]
若迭代次数t＝q
max
，则终止迭代并按下式输出空间矩阵h
x
和hy：
[0029][0030][0031]
否则图学习过程转到步骤2.1继续执行；
[0032]
步骤3：高光谱图像的解混
[0033]
步骤3.1、解混过程的外层循环：
[0034]
解混过程外层循环的权重矩阵w1和w2按下式计算：
[0035][0036]
其中，ε为解混过程设定的超参数，x
(l)
(i,:)表示在解混过程外层循环的第l次迭代时，估计矩阵x的第i行的所有列的元素所构成的行向量，表示解混过程内层循环的第l次迭代时，估计矩阵x的所有元素的绝对值；
[0037]
步骤3.2、变量更新：
[0038]
更新内层循环迭代次数：k＝k+1；
[0039]
将内层循环求解的x赋值给外层循环：x
(l+1)
＝x
(k+1)
；
[0040]
更新外层循环迭代次数：l＝l+1；
[0041]
步骤3.3、终止条件：
[0042]
若解混过程外层循环的迭代次数l＝l
max
，则终止迭代并输出矩阵x
(k+1)
；否则转到步骤3.1继续执行；当依据步骤3.3停止时，输出的x
(k+1)
就是从高光谱图像y中估计出的丰度矩阵。
[0043]
相比于现有技术，本发明具有如下优点：
[0044]
1、本发明与目前广泛用于高光谱遥感图像解混的其他稀疏解混方法相比，所提出的解混方法具有更好的解混性能和相对较低的计算代价。由于所提出的方法在高光谱图像的空间结构信息刻画方面更加准确，所以在很大程度上提升了高光谱图像中像素间的空间结构关系，可更精确的识别出实际场景中的物质在图像中的位置和所占的比例，能够使光谱信息相同或相似的像素能够被更好的划分到同一个同质区域内，这使得该方法在实际应用问题中能够更准确的识别出场景中实际存在的不同物质或材料，以及这些物质或材料在图像中相对更精确的位置和所占的比例。因此，本发明的方法具有精度相对较高、可靠性好的特点，适合推广使用。
[0045]
2、与现有的高光谱图像解混方法中引入空间信息所不同的是，本发明利用图学习的方式来引入空间信息，通过图学习方式所得到的空间矩阵与手动构建的空间矩阵相比，可自适应不同情况的图像且更加鲁棒和准确。
[0046]
3、本发明基于两个模拟高光谱图像数据集dc1、dc2和一个真实的铜矿石(cuprite)数据集进行实验，这些数据集在高光谱图像的解混任务中被广泛的用于定量分
析解混方法的性能。基于该数据集的实验表明，本发明中所提出的解混方法远好于当前的其他同类型的稀疏解混方法。
附图说明
[0047]
图1为采用现有的sunsal、sunsal-tv、s2wsu、mua-slic及sbglsu五种稀疏解混方法与本发明中所述的解混方法，在输入信噪比为10db时，对高光谱图像模拟数据集dc1进行解混时的性能表现图；
[0048]
图2为采用现有的sunsal、sunsal-tv、s2wsu、mua-slic及sbglsu五种稀疏解混方法与本发明中所述的解混方法，在输入信噪比为10db时，对高光谱图像模拟数据集dc2进行解混时的性能表现图；
[0049]
图3为采用现有的sunsal、sunsal-tv、s2wsu、mua-slic及sbglsu五种稀疏解混方法与本发明中所述的解混方法，对真实的铜矿石(cuprite)高光谱图像的数据集进行解混时的性能表现图；
[0050]
图中，sunsal、sunsal-tv、s2wsu、mua-slic以及sbglsu为五种稀疏解混对比方法所生成的丰度图，glbdwsu为本发明所提出的方法所生成的丰度图，ground truth为用于衡量解混性能的真实参考丰度图，tricorder maps为用于定量分析解混性能的参考丰度图。
具体实施方式
[0051]
下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。
[0052]
本发明提供了一种基于空间结构信息约束的高光谱图像解混方法，所述方法包括如下步骤：
[0053]
步骤1：高光谱图像的预处理
[0054]
为了同时减少高光谱图像谱域中的冗余信息和计算成本，首先将3-d的高光谱图像y分割为由相邻多个波段子图像所构成的图像集合其中，k为图像集合的数量，每个图像集合包含了原高光谱图像的相邻的多个波段的子图像，即：其中，nk为与第k个波段相邻的子图像所构成的集合，yj为原高光谱图像y的第j个波段的子图像；然后，每个图像集合yi由下式进行融合处理：
[0055][0056]
其中，|nk|表示集合的大小，包含了每个图像集合yk中的平均空间信息。
[0057]
高光谱图像预处理过程中，不同数据集的参数设置如下：
[0058]
模拟数据集dc1中图像集合数量k：k＝5。
[0059]
模拟数据集dc2中图像集合数量k：k＝5。
[0060]
真实的集铜矿石数据集中图像集合数量k：k＝20。
[0061]
步骤2：空间矩阵的图学习
[0062]
基于融合的图像所提出的空间矩阵的图学习算法的相关参数及输入量设置如下：
[0063]
所有高光谱图像数据集图学习过程的初始化设置：
[0064]
模拟数据集dc1的相关参数设置如下：高斯滤波器的标准差：σ＝1.2；图学习的正则参数：λ
l
＝0.003；超参数：ξ＝0.001；图学习过程的最大迭代次数：q
max
＝5。
[0065]
模拟数据集dc2的相关参数设置如下：高斯滤波器的标准差：σ＝0.5；图学习的正则参数：λ
l
＝0.008；超参数：ξ＝0.01；图学习过程的最大迭代次数：q
max
＝5。
[0066]
铜矿石数据集的相关参数设置如下：高斯滤波器的标准差：σ＝0.1；图学习的正则参数：λ
l
＝0.05；超参数：ξ＝15；图学习过程的最大迭代次数：q
max
＝10。
[0067]
具体步骤如下：
[0068]
步骤2.1、计算对角权重矩阵w
x
和wy，其对角元素按下式计算：
[0069][0070][0071]
其中，d
x
和dy分别为水平和垂直两个方向上的离散差分算子，g
σ
是标准差为σ的高斯滤波器，ξ为图学习过程设定的超参数，“*”为卷积算子，为基于元素运算的乘积算子，以及除法和绝对值都是基于元素的运算；
[0072]
步骤2.2、计算图laplacian矩阵
[0073][0074]
其中，d
x
和dy分别为水平和垂直两个方向上的离散差分算子，上角标
t
为矩阵的转置运算；
[0075]
步骤2.3、计算空间矩阵h
x
和h
y`
按下式计算：
[0076][0077][0078]
其中，w
1/2
表示对矩阵w的所有元素取平方根运算；
[0079]
步骤2.4、在图学习的第t次迭代过程中，变量按下式进行更新：
[0080][0081]
其中，in表示维数为n的单位阵，上角标-1
表示矩阵求逆运算；
[0082]
步骤2.5、图学习过程的终止条件：
[0083]
若迭代次数t＝q
max
，则终止迭代并按下式输出空间矩阵h
x
和hy：
[0084][0085][0086]
否则图学习过程转到步骤2.1继续执行。
[0087]
步骤3：高光谱图像的解混
[0088]
所有高光谱图像数据集解混过程的相关参数及输入量设置如下：
[0089]
模拟数据集dc1的相关参数设置：正则参数：λ1＝0.03，λ2＝0.0006；罚参数：μ＝0.2；超参数：ε＝0.0001；内外层循环的最大迭代次数分别设定为：k
max
＝10，l
max
＝10。
[0090]
模拟数据集dc2的相关参数设置：正则参数：λ1＝0.0004，λ2＝0.0004；罚参数：μ＝0.2；超参数：ε＝0.0001；内外层循环的最大迭代次数分别设定为：k
max
＝10，l
max
＝10。
[0091]
真实铜矿石数据集的相关参数设置：正则参数：λ1＝0.0003，λ2＝0.00003；罚参数：μ＝0.1；超参数：ε＝0.0001；内外层循环的最大迭代次数分别设定为：k
max
＝11，l
max
＝4。
[0092]
算法所要估计的所有数据集丰度矩阵的初始化x
(0)
设置如下：
[0093][0094]
式(8)中的变量说明：“svd”表示奇异值分解运算，上角标
t
表示矩阵的转置运算，a为已知的谱库矩阵，y为观测的高光谱图像，“diag(
·
)”表示构建一个对角阵，或以向量的形式返回一个矩阵对角线上的元素，μ为解混过程设定的罚参数。
[0095]
解混过程的辅助变量u
(0)
的初适化设置：u
(0)
＝x
(0)
；
[0096]
解混过程的辅助变量v
(0)
的初适化设置：v
(0)
＝x
(0)
；
[0097]
解混过程的辅助变量的初适化设置：
[0098]
解混过程的辅助变量的初适化设置：
[0099]
解混过程的lagrangian乘子初始化设置：为全零矩阵。
[0100]
具体步骤如下：
[0101]
步骤3.1、解混过程的外层循环：
[0102]
解混过程外层循环的权重矩阵w1和w2按下式计算：
[0103][0104]
其中，ε为解混过程设定的超参数，x
(l)
(i,:)表示在解混过程外层循环的第l次迭
代时，估计矩阵x的第i行的所有列的元素所构成的行向量，表示解混过程内层循环的第l次迭代时，估计矩阵x的所有元素的绝对值；
[0105]
具体步骤如下：
[0106]
步骤3.1.1、解混过程的内层循环：
[0107]
解混过程内层循环的变量x按下式进行计算：
[0108][0109]
其中，im表示维数为m的单位阵，上角标
t
表示矩阵的转置运算；
[0110]
步骤3.1.2、解混过程的内层循环：
[0111]
解混过程内层循环的变量u按下式进行计算：
[0112][0113]
其中，为基于元素运算的乘积算子，soft表示软阈值算子，定义为
[0114]
soft(u,θ)＝sign(u)
·
max(|u|-θ,0)
ꢀꢀꢀ
(14)；
[0115]
其中，sign为符号函数；
[0116]
步骤3.1.3、解混过程的内层循环：
[0117]
解混过程内层循环的变量v按下式进行计算：
[0118][0119]
其中，[v]
+
＝max(v,0)，in是维数为n的单位阵，上角标
t
表示矩阵的转置运算；
[0120]
步骤3.1.4、解混过程的内层循环：
[0121]
解混过程内层循环的变量z
x
和zy按下式进行计算：
[0122][0123]
其中，soft的定义方式与式(14)相同；
[0124]
步骤3.1.5、解混过程的内层循环：
[0125]
解混过程内层循环的lagrangian乘子按下式进行更新：
[0126][0127]
步骤3.2、必要的变量更新：
[0128]
更新内层循环迭代次数：k＝k+1；
[0129]
将内层循环求解的x赋值给外层循环：x
(l+1)
＝x
(k+1)
；
[0130]
更新外层循环迭代次数：l＝l+1。
[0131]
步骤3.3、算法的终止条件：
[0132]
若解混过程外层循环的迭代次数l＝l
max
，则终止迭代并输出矩阵x
(k+1)
；否则算法转到步骤3.1继续执行；当算法依据步骤3.3停止时，算法输出的x
(k+1)
就是从高光谱图像y中估计出的丰度矩阵。
[0133]
为验证本发明中所提出的解混方法的有效性，我们分别采用sunsal、sunsal-tv、s2wsu、mua-slic以及sbglsu五种稀疏解混方法和本发明中所述的解混方法(glbdwsu)，使用了两个合成的模拟数据集dc1、dc2和一个真实的铜矿石(cuprite)数据集进行解混处理。模拟数据集dc1的高光谱图像包含了224个波段，每个波段的图像由75
×
75的像素构成，包含25个以网格方式空间排列的正方形，其中每个正方形是一个同质区域，其像素以相同的丰度生成。模拟数据集dc2含有224个波段，每个波段的图像包含了100
×
100个像素，其通过使用以高斯随机场为中心的狄利克雷分布的相应丰度分数组成。真实的铜矿石(cuprite)数据集包含了188个波段，每个波段的图像由250
×
191的像素构成。
[0134]
我们在实验中，对于数据集dc1选择了第1端元、第3端元和第5端元的丰度图作为代表；对数据集dc2选择了第1端元、第3端元和第8端元的丰度图作为代表；对于真实的铜矿石数据集选择了明矾石(alunite)、布丁英石(buddingtonite)以及玉髓(chalcedony)这三种矿物作为代表，所提出的方法在每个数据集上的解混性能分别如图1、图2和图3所示。从图中可以看出，本发明所提出的解混方法的性能表现都远好于同类型的其他稀疏解混方法，此实验充分证实了所提出方法的有效性。
[0135]
此外，为了更进一步的说明本发明所提出方法的有效性，我们利用信号重构误差(signal-reconstruction-error,sre)来度量所提出解混方法的性能，该指标被广泛的应用于高光谱解混方法性能的度量中。信号重构误差sre的具体定义方式如下：
[0136][0137]
其中，x和分别为真实的丰度图和估计的丰度图，通常情况，sre值越大说明解混方法的性能越好。
[0138]
在实验中，我们分别对模拟数据集dc1和dc2加入了信噪比为10db的高斯噪声，以用于模拟含有噪声的高光谱图像的数据集。我们分别采用sunsal、sunsal-tv、s2wsu、mua-slic以及sbglsu五种稀疏解混方法和本发明中所述的解混方法，对加入10db高斯噪声的模拟数据集dc1和dc2进行了解混处理，所提出的方法的sre值如表1所示。从表1可以看出，本发明所提出的解混方法的性能表现都远好于同类型的其他稀疏解混方法，此实验更进一步的证实了所提出方法的有效性。
[0139]
表1各种方法在模拟数据集dc1、dc2上输入信噪比为10db时的解混sre值和计算时间
[0140][0141]
注：由于铜矿石数据是真实的高光谱图像数据，不存在真实的丰度值作为参考，因此没有列出关于铜矿石数据集的性能值的表格。