CN106529563A - 基于双图稀疏非负矩阵分解的高光谱波段选择方法 - Google Patents

Abstract

本发明公开了一种基于双图稀疏非负矩阵分解的高光谱波段选择方法。其处理步骤包括：输入需处理的高光谱数据；分别计算高光谱数据数据空间和特征空间的相似度矩阵；分别计算高光谱数据数据空间和特征空间的相似度对角矩阵；初始化高光谱数据重构矩阵；设置迭代次数；更新重构矩阵；判断是否达到最大迭代次数，若达到则得到数据空间的非负矩阵分解因子，否则返回更新重构矩阵继续迭代，直至达到最大迭代次数；构造并输出高光谱数据的波段选择矩阵。本发明解决了原始高光谱图像中存在许多冗余波段的问题，剔除了冗余信息，降低了数据维数，选择出具有更高判别性的波段，提高了高光谱图像分类的准确性。用于高光谱图像分类之前的波段选择处理。

Description

基于双图稀疏非负矩阵分解的高光谱波段选择方法

技术领域

本发明属于图像处理技术领域，涉及高光谱图像分类，具体是一种基于双图稀疏非负矩阵分解的高光谱波段选择方法，主要用于高光谱图像分类之前的波段选择处理。

背景技术

为了探测宽波段遥感中不可探测的物质，高光谱遥感技术得到了快速的发展，与此同时为了剔除冗余信息，降低高光谱图像数据维数，很多高光谱波段选择方法相继被提出。

2013年，V Kumar，J Hahn和AM Zoubir在他们发表的论文“BAND SELECTION FORHYPERSPECTRAL IMAGES BASED ON SELF-TUNING SPECTRAL CLUSTERING”中提出了一种基于谱聚类的自调整高光谱波段选择方法。该方法首先采用谱聚类将所有波段聚类为K类，然后用主成份分析法得到每一类中波段的协方差矩阵，并计算该协方差矩阵的特征值以及所对应的特征向量；按照一定的比例在各类中选取特征值以及所对应的特征向量作为类的特征向量基，各类的特征向量基组成一个变换矩阵，然后通过变换矩阵将高维的高光谱图像数据变换为低维的数据矩阵，达到降低数据维度的目的。这种方法将相似的波段分成一类，分块处理，很好的去除了冗余。然而谱聚类和变换矩阵是相互影响的，而该方法将谱聚类与变换矩阵分步进行计算，所以该方法没有用到高光谱图像数据的局部几何信息而且缺乏学习策略，不能很好地选择出具有判别性的波段。从数据原始的物理意义上来看，该方法属于特征提取方法，对原始数据采用一个变换矩阵进行低维映射来实现降维的目的，缺乏一定的物理意义。

西安电子科技大学张向荣等人申请的专利“基于低秩表示的高光谱图像波段选择方法”(专利申请号：CN201510411250.9，公开号：CN105046276A)提出了一种低秩表示聚类的高光谱图像波段选择方法。该方法对高光谱图像进行归一化处理并进行低秩表示，再使用增强拉格朗日乘子法来求解低秩表示系数，然后基于低秩表示系数对所有波段进行聚类，最后从每个聚类中选择出最具代表性的波段作为最终选择的波段，并对其进行聚类。该方法可以去除波段中的冗余信息，选择具代表性的波段用于分类，降低了高光谱图像数据的维度，提高了高光谱图像的分类准确度。然而在进行低秩表示系数求解学习时，没有利用高光谱图像数据的局部几何信息，而且没有用到特征的自表示信息，所以导致得到的低秩表示系数不够准确，最终选择的波段也缺乏一定的代表性。

现有的研究方法中大多只采用了数据空间的局部几何信息，并不能充分挖掘高光谱图像的潜在信息，另外没有加入稀疏约束，无法保证矩阵的稀疏性，也没有利用特征自表示信息这一重要的特征属性信息。

发明内容

本发明的目的在于解决上述现有技术存在的问题，提出一种更具波段判别性、高光谱图像分类准确性更高的基于双图稀疏非负矩阵分解的高光谱波段选择方法。

本发明是一种基于双图稀疏非负矩阵分解的高光谱波段选择方法，其特征在于，包括如下步骤：

(1)输入需要处理的高光谱图像数据，并预设所需波段选择数目；

(2)分别计算高光谱图像数据的数据空间和特征空间的相似度矩阵：

采用权重测度算法，分别计算高光谱图像数据的数据空间和特征空间中所有波段的波段相似度，得到高光谱图像数据的数据空间和特征空间中所有波段的波段相似度矩阵；

(3)分别计算高光谱图像数据的数据空间和特征空间的相似度对角矩阵：

对高光谱图像数据的数据空间和特征空间中所有波段的波段相似度矩阵分别进行对角化处理，得到高光谱图像数据的数据空间和特征空间的相似度对角矩阵；

(4)初始化高光谱图像数据的重构矩阵：

采用随机矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子、高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，采用单位矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵三个重构矩阵；

(5)设置循环迭代次数：

将初始迭代次数设置为0，最大迭代次数设置为10～30；

(6)更新高光谱图像数据的重构矩阵：

(6a)利用高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的更新公式，得到当前迭代次数下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子；

(6b)利用高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子的更新公式，得到当前迭代次数下更新的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子；

(6c)利用高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵的更新公式，得到当前迭代次数下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵；

(7)判断是否达到最大迭代次数，若满足则执行步骤(8)，否则将当前循环迭代次数加1并返回执行步骤(6)直至达到最大迭代次数；

(8)得到高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子；

(9)构造并输出高光谱图像数据的波段选择矩阵：

(9a)采用波段评价向量公式，基于高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，计算高光谱图像数据的波段评价向量；

(9b)将高光谱图像数据的波段评价向量中的元素值降序排序，并按预设的波段选择数目从中选出元素值最大的波段构造成高光谱图像数据的波段选择矩阵，输出高光谱图像数据的波段选择矩阵。

本发明将双图(数据图和特征图)模型与非负矩阵分解方法相结合引入到高光谱图像的波段选择上，并对数据空间的非负矩阵分解因子加入稀疏约束，即稀疏自表示项。数据空间和特征空间的非负矩阵分解因子可以在交替迭代的过程中相互作用更新，充分发挥双图模型的效果。依据计算得到的波段评价向量从原始高光谱图像中选择出具有判别性的波段，保留了原始高光谱图像的物理意义，同时提高了高光谱图像分类的准确性。

本发明与现有技术相比具有以下优点：

第一，本发明计算波段相似度矩阵充分利用了数据空间和特征空间的局部几何信息，尤其是引入了特征空间的局部几何信息，解决了现有方法不能充分挖掘高光谱图像的潜在信息的问题，提高了高光谱图像分类的准确性。

第二，本发明将双图模型与非负矩阵分解方法相结合，并引入到高光谱图像的波段选择上，利用高光谱图像数据的重构矩阵的更新公式，使得数据空间的含特征空间局部几何信息的非负矩阵分解因子和特征空间的含数据空间局部几何信息的非负矩阵分解因子在交替迭代更新的过程中相互作用更新，可以充分发挥双图模型的效果，选择出更具判别性的波段，提高了高光谱图像分类的准确性。

第三，本发明对数据空间的含特征空间局部几何信息的非负矩阵分解因子进行稀疏约束，即体现特征自表示信息的稀疏自表示项，不仅反映了在特征属性信息中的重要性，而且保证了数据空间中的含特征空间局部几何信息的非负矩阵分解因子的稀疏性。

附图说明

图1是本发明的流程框图；

图2是本发明与现有技术的仿真实验结果对比曲线图；

图3是本发明与现有技术的仿真效果彩图，其中图3(a)为原始的Indian Pines高光谱图像的真值彩图，图3(b)为波段全选方法(ALL)的KNN分类效果彩图，图3(c)为基于谱聚类的自调整波段选择方法(SC)选择50个波段时的KNN分类效果彩图，图3(d)为本发明(DSNMF)选择50个波段时的KNN分类效果彩图。

图4是本发明与现有技术的仿真效果灰度图，其中图4(a)为原始的IndianPines高光谱图像的真值灰度图，图4(b)为波段全选方法的KNN分类效果灰度图，图4(c)为基于谱聚类的自调整波段选择方法选择50个波段时的KNN分类效果灰度图，图4(d)为本发明选择50个波段时的KNN分类效果灰度图。

具体实施方式

下面结合附图对本发明做详细描述。

与遥感图像相比，高光谱图像不仅具有更丰富的信息，而且在处理技术上，为光谱数据进行更有效的分析处理提供了可能，所以高光谱图像不仅引起了遥感界的兴趣，同时受到了其它领域(如农学、医学等)的极大关注。高光谱图像有着非常广阔的应用前景，尤其是在高光谱图像的分类方面，而高光谱图像分类前的波段选择对于高光谱图像分类起着举足轻重的作用，因此很多高光谱波段选择方法相继被提出。然而现有的波段选择方法大多存在只采用了数据空间的局部几何信息，没有利用特征空间的局部几何信息等问题，不能选择出更具判别性的波段，因此高光谱图像分类也不够准确。为此本发明经过进一步的分析和创新，提出了一种基于双图稀疏非负矩阵分解的高光谱波段选择方法。

实施例1

本发明是一种基于双图稀疏非负矩阵分解的高光谱波段选择方法，参见附图1，包括如下步骤：

(1)输入需要处理的高光谱图像数据，并预设所需波段选择数目：

本发明中用户可以根据自身需求选择相应的波段选择数目，比如参见附图2，可从中选择5-50个波段选择数目，一般分类准确度随着波段选择数目的增加而增加。

(2)分别计算高光谱图像数据的数据空间和特征空间的相似度矩阵：

采用权重测度算法，本发明中采用计算高光谱图像数据的数据空间中所有波段的波段之间的欧式距离这一重要指标来评价波段相似度，得到高光谱图像数据的数据空间中所有波段的波段相似度矩阵；采用权重测度算法，本发明中采用计算高光谱图像数据的特征空间中所有波段的波段之间的欧式距离这一重要指标来评价波段相似度，得到高光谱图像数据的特征空间中所有波段的波段相似度矩阵。

(3)分别计算高光谱图像数据的数据空间和特征空间的相似度对角矩阵：

对步骤(2)中得到的高光谱图像数据的数据空间中所有波段的波段相似度矩阵进行对角化处理，得到高光谱图像数据的数据空间的相似度对角矩阵；对步骤(2)中对高光谱图像数据特征空间中所有波段的波段相似度矩阵进行对角化处理，得到高光谱图像数据的特征空间的相似度对角矩阵。

本发明在步骤(2)、(3)中都采用了一种双图模型，同时引入了数据空间和特征空间的局部几何信息，有利于充分挖掘高光谱图像数据的潜在信息，选择出具有更高判别性的波段，提高高光谱图像分类的准确性。

(4)初始化高光谱图像数据的重构矩阵：

采用随机矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子、高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，本发明在Matlab R2016a运行环境下采用rand函数进行初始化；采用单位矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵，在Matlab R2016a运行环境下采用eye函数进行初始化。

(5)设置循环迭代次数：

将初始迭代次数设置为0，最大迭代次数设置为10～30。当循环迭代次数较小时，运算时间较短，但高光谱图像分类的准确性有限；当循环迭代次数较大时，运算时间较长，但高光谱图像分类的准确性较高。综合考虑以上因素，本例中设置最大迭代次数为20。

(6)更新高光谱图像数据的重构矩阵：

(6a)利用高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的更新公式，得到当前迭代次数下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子。

(6b)利用高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子的更新公式，得到当前迭代次数下更新的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子。

(6c)利用高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵的更新公式，得到当前迭代次数下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵。

本发明在数据空间的非负矩阵分解因子的更新公式和特征空间的非负矩阵分解因子的更新公式中，不仅利用到上次循环迭代的数据空间和特征空间的局部几何信息，而且可以相互更新作用，充分发挥双图模型的优势，选择出具有更高判别性的波段，提高高光谱图像分类的准确性。

(7)判断是否达到最大迭代次数，若满足则执行步骤(8)，否则将当前循环迭代次数加1并返回执行步骤(6)，直至达到最大迭代次数。

(8)得到高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子。

(9)构造并输出高光谱图像数据的波段选择矩阵：

(9a)采用波段评价向量公式，基于高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，计算高光谱图像数据的波段评价向量。

(9b)将高光谱图像数据的波段评价向量中的元素值降序排序，并按预设的波段选择数目从中选出元素值最大的波段构造成高光谱图像数据的波段选择矩阵，并输出高光谱图像数据的波段选择矩阵。

本发明利用高光谱图像数据分别计算高光谱图像数据的数据空间和特征空间的相似度矩阵和相似度对角矩阵。利用高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的更新公式、高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子的更新公式、高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵的更新公式相互作用。也就是说在更新高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的过程中，用到了特征空间的局部几何信息，在更新高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子的过程中，用到了数据空间的局部几何信息。由此不断地进行迭代更新，充分发挥双图模型的优势，有利于挖掘高光谱图像数据的潜在信息，提高高光谱图像分类的准确性。

实施例2

基于双图稀疏非负矩阵分解的高光谱波段选择方法同实施例1，本例中设置最大迭代次数为10，其中步骤(2)中所述的权重测度算法如下：

2.1：根据下面两个式子，分别计算高光谱图像数据的数据空间和特征空间的所有波段之间的欧式距离：

数据空间的所有波段之间的欧式距离：

特征空间的所有波段之间的欧式距离：

其中，O^P和O^S分别表示高光谱图像数据的数据空间和特征空间的所有波段之间的欧式距离，表示开方操作，X表示高光谱图像数据，*表示Hadamard矩阵乘积操作，1_n×d表示n×d的全1矩阵，1_d×n表示d×n的全1矩阵，n表示高光谱图像数据的像素点总数，d表示高光谱图像数据的波段总数，^T表示转置操作。

2.2：根据下面两个式子，分别计算高光谱图像数据的数据空间和特征空间的相似度矩阵；

数据空间的相似度矩阵：

特征空间的相似度矩阵：

其中，W^P和W^S分别表示高光谱图像数据的数据空间和特征空间的相似度矩阵，exp(·)表示指数操作，O^P和O^S分别表示高光谱图像数据的数据空间和特征空间的所有波段之间的欧式距离，σ表示高斯尺度参数。

现有的权重测度算法有很多，本发明中利用的欧式距离是一种常用的评价波段相似度的测度算法，此方法物理意义明确而且计算简单。

实施例3

基于双图稀疏非负矩阵分解的高光谱波段选择方法同实施例1-2，本例中设置最大迭代次数为15，其中步骤(3)中所述的对角化处理是按照下面两个式子求得：

数据空间的相似度对角矩阵：

特征空间的相似度对角矩阵：

其中，D^P和D^S分别表示高光谱图像数据的数据空间和特征空间的相似度对角矩阵，diag(·)表示生成对角矩阵操作，∑表示叠加操作，[W^P]_i表示高光谱图像数据的数据空间的相似度矩阵的第i列，d表示高光谱图像数据的波段总数，[W^S]_j表示高光谱图像数据的特征空间的相似度矩阵的第j列，n表示高光谱图像数据的像素点总数。

相对于现有的波段选择方法，本发明中采用的双图模型，同时利用了数据空间和特征空间的局部几何信息，所以在计算相似度矩阵以及相似度对角矩阵时都在数据空间和特征空间这两个空间进行计算，即使用了双图模型。

实施例4

基于双图稀疏非负矩阵分解的高光谱波段选择方法同实施例1-3，本例中设置最大迭代次数为25，其中步骤(4)中所述的对重构矩阵初始化包括如下步骤：

4.1：采用随机矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子为d×k的随机矩阵，k表示所选取的高光谱图像数据的波段选择数目目，d表示高光谱图像数据的波段总数；

4.2：采用随机矩阵法初始化高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子为n×k的随机矩阵，k表示所选取的高光谱图像数据的波段选择数目目，n表示高光谱图像数据的像素点总数；

4.3：采用单位矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵为d×d的单位矩阵，d表示高光谱图像数据的波段总数。

本发明前两个重构矩阵即步骤4.1和4.2中所完成的矩阵是采用非负矩阵分解(NMF)方法分解的两个非负矩阵分解因子，第三个重构矩阵即步骤4.3中所完成的矩阵则主要体现了特征的自表示信息。

实施例5

基于双图稀疏非负矩阵分解的高光谱波段选择方法同实施例1-4，本例中设置最大迭代次数为30，其中步骤(6a)中所述的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的更新公式如下：

其中，P^(t+1)表示第t+1次迭代下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，P^(t)表示第t次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，S^(t)表示第t次迭代下的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，U^(t)表示第t次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵，X表示高光谱图像数据，W^P表示高光谱图像数据的数据空间的相似度矩阵，D^P表示高光谱图像数据的数据空间的相似度对角矩阵，^T表示转置操作，α为数据空间的正则参数，取值范围为{10^-3,10^-2,10^-1,10⁰,10¹,10²,10³}，θ为稀疏参数，取值范围为{300,800,2000,4000,6000,8000}。在实际的实验中，用户可根据具体的高光谱图像自行调整数据空间的正则参数和稀疏参数设置。

本发明的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的更新公式中体现了双图模型和非负矩阵分解的结合，在迭代过程中，不仅利用到了数据空间的局部几何信息，也利用到了特征空间的局部几何信息。

实施例6

基于双图稀疏非负矩阵分解的高光谱波段选择方法同实施例1-5，其中步骤(6b)中所述的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子的更新公式如下：

其中，S^(t+1)表示第t+1次迭代下更新的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，S^(t)表示第t次迭代下的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，P^(t)表示第t次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，X表示高光谱图像数据，W^S表示高光谱图像数据的特征空间的相似度矩阵，D^S表示高光谱图像数据的特征空间的相似度对角矩阵，β为特征空间的正则参数，取值范围为{10^-3,10^-2,10^-1,10⁰,10¹,10²,10³}，^T表示转置操作。在实际的实验中，用户可根据具体的高光谱图像自行调整特征空间的正则参数设置。

与实施例5相似，本发明中采用的双图模型和非负矩阵分解相结合的方法，在高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子的迭代公式中可以看出，不仅用到特征空间的局部几何信息，而且利用到了数据空间的含特征空间局部几何信息的非负矩阵分解因子。因此，本发明充分发挥了双图模型的优势，提高高光谱图像分类的准确性。

实施例7

基于双图稀疏非负矩阵分解的高光谱波段选择方法同实施例1-6，其中步骤(9a)中所述的波段评价向量公式如下：

其中，G表示高光谱图像数据的波段评价向量，P表示高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，*表示Hadamard矩阵乘积操作，[P*P]_i表示[P*P]的第i列，k表示所选取的高光谱图像数据的波段选择数目目，∑表示叠加操作。

根据求出的波段评价向量，将高光谱图像数据的波段评价向量中的元素值降序排序，并按预设所需波段选择数目选择出元素值最大的波段构造成高光谱图像数据的波段选择矩阵，并输出高光谱图像数据的波段选择矩阵。最后依据高光谱图像数据的波段选择矩阵对高光谱图像数据进行分类，不仅可以解决原始高光谱图像中存在许多冗余波段的问题，剔除冗余信息，降低数据维数，而且可以选择出具有更高判别性的波段，提高高光谱图像分类的准确性。

下面给出一个完整的例子，结合附图对本发明做进一步详细描述。

实施例8

基于双图稀疏非负矩阵分解的高光谱波段选择方法同实施例1-7，参见附图1，具体步骤如下。

步骤1，输入需要处理的高光谱图像数据，并预设所需波段选择数目。

在本发明实验中，输入的高光谱图像数据是根据高光谱图像Indian Pines得到的，原始图像参见附图4(a)。通常情况下，高光谱图像分类的准确性随着波段选择数目的增加而增大。假设需要与现有技术在不同波段选择数目下进行分类准确度对比时，可以改变波段选择数目，比如针对高光谱图像Indian Pines可以选择5-50个波段。假设需要与现有技术进行高光谱图像分类效果图像的对比时，可以固定波段选择数目进行分类效果图像的对比，比如针对高光谱图像Indian Pines固定选择50个波段。

步骤2，分别计算高光谱图像数据的数据空间和特征空间的相似度矩阵。

利用权重测度算法，分别计算高光谱图像数据的数据空间和特征空间的波段相似度，得到高光谱图像数据的数据空间和特征空间的相似度矩阵。

权重测度算法如下：

2.1：根据下面两个式子，分别计算高光谱图像数据的数据空间和特征空间的所有波段之间的欧式距离：

数据空间的相似度矩阵：

特征空间的相似度矩阵：

其中，O^P和O^S分别表示高光谱图像数据的数据空间和特征空间的所有波段之间的欧式距离，表示开方操作，X表示高光谱图像数据，*表示Hadamard矩阵乘积操作，1_n×d表示n×d的全1矩阵，1_d×n表示d×n的全1矩阵，n表示高光谱图像数据的像素点总数，d表示高光谱图像数据的波段总数，^T表示转置操作；

2.2：根据下面两个式子，分别计算高光谱图像数据的数据空间和特征空间的相似度矩阵；

数据空间的相似度矩阵：

特征空间的相似度矩阵：

其中，W^P和W^S分别表示高光谱图像数据的数据空间和特征空间的相似度矩阵，exp(·)表示指数操作，O^P和O^S分别表示高光谱图像数据的数据空间和特征空间的所有波段之间的欧式距离，σ表示高斯尺度参数。

步骤3，分别计算高光谱图像数据的数据空间和特征空间的相似度对角矩阵。

分别对高光谱图像数据的数据空间和特征空间的相似度矩阵进行对角化处理，得到高光谱图像数据的数据空间和特征空间的相似度对角矩阵。

根据下面两个式子进行对角化处理：

数据空间的相似度对角矩阵：

特征空间的相似度对角矩阵：

其中，D^P和D^S分别表示高光谱图像数据的数据空间和特征空间的相似度对角矩阵，diag(·)表示生成对角矩阵操作，∑表示叠加操作，[W^P]_i表示高光谱图像数据的数据空间的相似度矩阵的第i列，d表示高光谱图像数据的波段总数，[W^S]_j表示高光谱图像数据的数据空间的相似度矩阵的第j列，n表示高光谱图像数据的像素点总数。

本发明中采用的双图模型，对数据空间和特征空间的相似度矩阵以及相似度对角矩阵同时进行计算，充分挖据了高光谱图像数据的局部几何信息，提高了分类的准确性，而且可以进行并行处理，提高运算效率。

步骤4，初始化高光谱图像数据的重构矩阵。

对高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子、高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子以及高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵三个重构矩阵进行初始化。

对重构矩阵初始化的具体步骤如下：

4.1：采用随机矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子为d×k的随机矩阵，k表示所选取的高光谱图像数据的波段选择数目目，d表示高光谱图像数据的波段总数；

4.2：采用随机矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子为n×k的随机矩阵，k表示所选取的高光谱图像数据的波段选择数目目，n表示高光谱图像数据的像素点总数；

4.3：采用单位矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵为d×d的单位矩阵，d表示高光谱图像数据的波段总数。

具体实施中，采用Matlab R2016a软件中的rand函数对高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子和高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子进行初始化，采用Matlab R2016a软件中的eye函数对高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵进行初始化。

步骤5，设置循环迭代次数。

将初始迭代次数设置为0，最大迭代次数设置为20。

步骤6，更新高光谱图像数据的重构矩阵。

6.1：利用高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的更新公式，得到当前迭代次数下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子。

高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的更新公式如下：

其中，P^(t+1)表示第t+1次迭代下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，P^(t)表示第t次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，S^(t)表示第t次迭代下的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，U^(t)表示第t次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵，X表示高光谱图像数据，W^P表示高光谱图像数据的数据空间的相似度矩阵，D^P表示高光谱图像数据的数据空间的相似度对角矩阵，α为数据空间的正则参数，取值范围为{10^-3,10^-2,10^-1,10⁰,10¹,10²,10³}，θ为稀疏参数，取值范围为{300,800,2000,4000,6000,8000}，^T表示转置操作。

6.2：利用高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子的更新公式，得到当前迭代次数下更新的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子。

高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子的更新公式如下：

其中，S^(t+1)表示第t+1次迭代下更新的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，S^(t)表示第t次迭代下的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，P^(t)表示第t次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，X表示高光谱图像数据，W^S表示高光谱图像数据的特征空间的相似度矩阵，D^S表示高光谱图像数据的特征空间的相似度对角矩阵，β为特征空间的正则参数，其取值范围为{10^-3,10^-2,10^-1,10⁰,10¹,10²,10³}，^T表示转置操作。

6.3：利用高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵的更新公式，得到当前迭代次数下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵。

高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵的更新公式如下：

其中，U^(t+1)表示第t+1次迭代下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵，P^(t+1)表示第t+1次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，*表示Hadamard矩阵乘积操作，[P^(t+1)*P^(t+1)]_i表示[P^(t+1)*P^(t+1)]的第i列，∑表示叠加操作，表示开方操作，diag(·)表示生成对角矩阵操作，k表示所选取的高光谱图像数据的波段选择数目目。

本发明对数据空间的含特征空间局部几何信息的非负矩阵分解因子加入了体现特征自表示信息的稀疏约束，即稀疏自表示项，这样不仅体现了在特征属性信息中的各特征的重要性程度，而且保证了数据空间中的含特征空间局部几何信息的非负矩阵分解因子的稀疏性。

步骤7，判断是否达到最大迭代次数，若满足则执行步骤8，否则将当前循环迭代次数加1并执行步骤6。

步骤8，得到高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子。

步骤9，构造并输出高光谱图像数据的波段选择矩阵。

9.1：采用波段评价向量公式，基于高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，计算高光谱图像数据的波段评价向量。

波段评价向量公式如下：

其中，G表示高光谱图像数据的波段评价向量，P表示高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，*表示Hadamard矩阵乘积操作，[P*P]_i表示[P*P]的第i列，k表示所选取的高光谱图像数据的波段选择数目目，∑表示叠加操作。

9.2：将高光谱图像数据的波段评价向量中的元素值降序排序，并按预设所需波段选择数目选择出元素值最大的波段构造成高光谱图像数据的波段选择矩阵，并输出高光谱图像数据的波段选择矩阵。

参见附图2可知，采用高光谱图像Indian Pines，固定选择50个波段，用本发明处理的波段选择矩阵通过KNN分类器后，分类准确度可以达到0.7030，也就是说在同样的KNN分类器下，采用本发明很明显要优于其他两种方法。采用本发明处理后的分类效果图像参见附图3(d)和附图4(d)，其中错分点比其他方法少，其他两种方法的分类效果图像参见附图3(b)、附图4(b)以及附图3(c)、附图4(c)。尤其参见上述图像的右下角区域，在附图3(b)、附图4(b)的右下角区域分布有数十处大大小小的错分点，有参见附图3(c)、附图4(c)的右下角区域同样也分布有数十处大大小小的错分点。而本发明处理后的分类效果较好，参见附图3(d)、附图4(d)，本发明在该区域的错分点只有两三处。纵观整体图像，采用本发明处理后的分类效果图像的错分点数目也明显少于其他两种方法处理后的分类效果图像。因此可知，本发明可以选择出更具判别性的波段，为高效的高光谱图像分类奠定了基础，是一种合理有效的基于双图稀疏非负矩阵分解的高光谱波段选择方法。

下面结合仿真实验对本发明再做进一步详细描述。

实施例9

基于双图稀疏非负矩阵分解的高光谱波段选择方法同实施例1-8。

仿真实验条件：

本发明的仿真实验采用的硬件测试平台是：处理器为Inter Core i5，主频为2.50GHz，内存8GB；软件平台为：Windows 7旗舰版64位操作系统、Matlab R2016a进行仿真测试。

仿真实验内容：

对通过光谱仪AVIRIS获取的Indian Pines高光谱图像进行波段选择和分类算法仿真。

本发明的仿真实验采用的Indian Pines高光谱图像，图中包含16类主要植被，每类植被含有220个波段。实验中，主要针对16类主要植被进行研究，并且要除去水吸收波段，因此实验中只用到10366个像素点和200个波段，即得到200×10366的高光谱图像数据。

用本发明对图像数据矩阵进行波段选择，并将选择后的结果用KNN工具箱中的KNN分类器进行分类，用分类准确率来验证波段选择效果，实验中随机选取7％的像素点作为训练样本，其余的作为测试样本，独立运行10次取平均值作为分类结果。

仿真结果分析：

图2是本发明与现有技术的实验结果对比图，图2显示了KNN分类器对波段全选方法(ALL)、基于谱聚类的自调整波段选择方法(SC)以及本发明(DSNMF)的选择结果进行分类所得到的分类准确度对比。图2中的横坐标表示波段选择数目k，纵坐标表示分类准确度(AC)。图2中的以正方形标记的曲线表示波段全选方法仿真实验的结果，以星号标记的曲线代表基于谱聚类的自调整波段选择方法仿真实验的结果，以圆圈标记的曲线代表本发明仿真实验的结果。

由图2可以看出，本发明在大部分波段选择数目下的分类结果都优于基于谱聚类的自调整波段选择方法，尤其是在波段选择数目较大时，而且随着波段选择数目的增加，分类准确度一直处于增加的状态。在波段选择数目为45和50时，分类准确度可以达到0.6937和0.7030，高于波段全选方法，这表示本发明不仅可以实现对数据的降维，而且还可以提高分类准确度。

本发明针对原始高光谱图像中存在许多冗余波段的问题，提出的高光谱波段选择方法，剔除了冗余信息，降低了数据维数，同时选择出具有更高判别性的波段，提高了高光谱图像分类的准确性。

实施例10

基于双图稀疏非负矩阵分解的高光谱波段选择方法同实施例1-8，仿真实验条件和仿真内容同实施例9：

图3和图4分别是本发明与现有技术的仿真效果彩图和灰度图，图3(a)和图4(a)为原始的Indian Pines高光谱图像的真值彩图和灰度图，图3(b)和图4(b)为波段全选方法的KNN分类效果彩图和灰度图，图3(c)和图4(c)为基于谱聚类的自调整波段选择方法选择50个波段时的KNN分类效果彩图和灰度图，图3(d)和图4(d)为本发明选择50个波段时的KNN分类效果彩图和灰度图。由图3和图4可以看出，本发明得到的分类效果图更接近原始的Indian Pines高光谱图像的真值图。就区域间的划分界限来说，采用本发明处理后的分类效果图像明显清晰于其他两种方法。就分类的错分点来说，采用本发明处理后的分类效果图像的错分点数目也明显少于采用其他两种方法处理后的分类效果图像，尤其是右下角的那块区域，错分点很少，分类准确度很高，充分说明了本发明可以得到很好的波段选择结果。

综上所述，本发明提出的一种基于双图稀疏非负矩阵分解的高光谱波段选择方法，具体处理步骤如下：(1)输入需要处理的高光谱图像数据；(2)分别计算高光谱图像数据的数据空间和特征空间的相似度矩阵；(3)分别计算高光谱图像数据的数据空间和特征空间的相似度对角矩阵；(4)初始化高光谱图像数据的重构矩阵；(5)设置循环迭代次数；(6)更新高光谱图像数据的重构矩阵；(7)判断是否达到最大迭代次数，若达到最大迭代次数则执行步骤(8)，否则将当前循环迭代次数加1并执行步骤(6)；(8)得到高光谱图像数据的数据空间的非负矩阵分解因子(9)构造并输出高光谱图像数据的波段选择矩阵。本发明中不仅采用了一种新的双图模型并引入了稀疏自表示项高光谱波段选择，而且将非负矩阵分解的方法和双图模型相结合，在更新数据空间的非负矩阵分解因子时利用到了特征空间的局部几何信息，在更新特征空间的非负矩阵分解因子时利用到了数据空间的局部几何信息，因此它们在交替迭代的过程中可以相互作用更新，充分发挥双图模型的效果。本发明解决了原始高光谱图像中存在许多冗余波段的问题，剔除了冗余信息，降低了数据维数，选择出具有更高判别性的波段，提高了高光谱图像分类的准确性。

Claims (8)

1.一种基于双图稀疏非负矩阵分解的高光谱波段选择方法，其特征在于，包括如下步骤：

(1)输入需要处理的高光谱图像数据，并预设所需波段选择数目；

(2)分别计算高光谱图像数据的数据空间和特征空间的相似度矩阵：

采用权重测度算法，分别计算高光谱图像数据的数据空间和特征空间中所有波段的波段相似度，得到高光谱图像数据的数据空间和特征空间中所有波段的波段相似度矩阵；

(3)分别计算高光谱图像数据的数据空间和特征空间的相似度对角矩阵：

对高光谱图像数据的数据空间和特征空间中所有波段的波段相似度矩阵分别进行对角化处理，得到高光谱图像数据的数据空间和特征空间的相似度对角矩阵；

(4)初始化高光谱图像数据的重构矩阵：

采用随机矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子、高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，采用单位矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵三个重构矩阵；

(5)设置循环迭代次数：

将初始迭代次数设置为0，最大迭代次数设置为10～30；

(6)更新高光谱图像数据的重构矩阵：

(6a)利用高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的更新公式，得到当前迭代次数下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子；

(6b)利用高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子的更新公式，得到当前迭代次数下更新的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子；

(6c)利用高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵的更新公式，得到当前迭代次数下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵；

(7)判断是否达到最大迭代次数，若满足则执行步骤(8)，否则将当前循环迭代次数加1并返回执行步骤(6)直至达到最大迭代次数；

(8)得到高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子；

(9)构造并输出高光谱图像数据的波段选择矩阵：

(9a)采用波段评价向量公式，基于高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，计算高光谱图像数据的波段评价向量；

(9b)将高光谱图像数据的波段评价向量中的元素值降序排序，并按预设的波段选择数目从中选出元素值最大的波段构造成高光谱图像数据的波段选择矩阵，并输出高光谱图像数据的波段选择矩阵。

2.根据权利要求1所述的基于双图稀疏非负矩阵分解的高光谱波段选择方法，其特征在于，步骤(2)中所述的权重测度算法如下：

2.1：分别计算高光谱图像数据的数据空间和特征空间的所有波段之间的欧式距离：

数据空间的所有波段之间的欧式距离：

$O^P=\sqrt{(X*X)1_{n\×d}+1_{d\×n}(X^T*X^T)-2{\mathrm{XX}}^T}$

特征空间的所有波段之间的欧式距离：

$O^S=\sqrt{(X^T*X^T)1_{d\×n}+1_{n\×d}(X*X)-2X^{T}X}$

其中，O^P和O^S分别表示高光谱图像数据的数据空间和特征空间的所有波段之间的欧式距离，表示开方操作，X表示高光谱图像数据，*表示Hadamard矩阵乘积操作，1_n×d表示n×d的全1矩阵，1_d×n表示d×n的全1矩阵，n表示高光谱图像数据的像素点总数，d表示高光谱图像数据的波段总数，^T表示转置操作；

2.2：分别计算高光谱图像数据的数据空间和特征空间的相似度矩阵；

数据空间的相似度矩阵：

$W^P=\exp (-\frac{O^P}{2{\mathrm{\σ}}^2})$

特征空间的相似度矩阵：

$W^S=\exp (-\frac{O^S}{2{\mathrm{\σ}}^2})$

其中，W^P和W^S分别表示高光谱图像数据的数据空间和特征空间的相似度矩阵，exp(.)表示指数操作，O^P和O^S分别表示高光谱图像数据的数据空间和特征空间的所有波段之间的欧式距离，σ表示高斯尺度参数。

3.根据权利要求1所述的基于双图稀疏非负矩阵分解的高光谱波段选择方法，其特征在于，步骤(3)中所述的对角化处理是指按照下面两个式子求得：

数据空间的相似度对角矩阵：

$D^P=diag\left(\underset{i=1}{\overset{d}{\Σ}}{\[W^P\]}_i\right)$

特征空间的相似度对角矩阵：

$D^S=diag\left(\underset{j=1}{\overset{n}{\Σ}}{\[W^S\]}_j\right)$

其中，D^P和D^S分别表示高光谱图像数据的数据空间和特征空间的相似度对角矩阵，diag(.)表示生成对角矩阵操作，∑表示叠加操作，[W^P]_i表示高光谱图像数据的数据空间的相似度矩阵的第i列，d表示高光谱图像数据的波段总数，[W^S]_j表示高光谱图像数据的特征空间的相似度矩阵的第j列，n表示高光谱图像数据的像素点总数。

4.根据权利要求1所述的基于双图稀疏非负矩阵分解的高光谱波段选择方法，其特征在于，步骤(4)中所述的对重构矩阵初始化包括如下步骤：

4.1：采用随机矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子为d×k的随机矩阵，k表示所选取的高光谱图像数据的波段数目，d表示高光谱图像数据的波段总数；

4.2：采用随机矩阵法初始化高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子为n×k的随机矩阵，k表示所选取的高光谱图像数据的波段数目，n表示高光谱图像数据的像素点总数；

4.3：采用单位矩阵法初始化高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵为d×d的单位矩阵，d表示高光谱图像数据的波段总数。

5.根据权利要求1所述的基于双图稀疏非负矩阵分解的高光谱波段选择方法，其特征在于，步骤(6a)中所述的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的更新公式如下：

$P^{(t+1)}=P^{\left(t\right)}\frac{{\mathrm{XS}}^{\left(t\right)}+{\mathrm{\αW}}^P{P}^{\left(t\right)}}{P^{\left(t\right)}{\left(S^{\left(t\right)}\right)}^T{S}^{\left(t\right)}+{\mathrm{\αD}}^P{P}^{\left(t\right)}+{\mathrm{\θU}}^{\left(t\right)}{P}^{\left(t\right)}}$

其中，P^(t+1)表示第t+1次迭代下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，P^(t)表示第t次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，S^(t)表示第t次迭代下的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，U^(t)表示第t次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵，X表示高光谱图像数据，W^P表示高光谱图像数据的数据空间的相似度矩阵，D^P表示高光谱图像数据的数据空间的相似度对角矩阵，^T表示转置操作，α为数据空间的正则参数，取值范围为{10^-3,10^-2,10^-1,10⁰,10¹,10²,10³}，θ为稀疏参数，取值范围为{300,800,2000,4000,6000,8000}。

6.根据权利要求1所述的基于双图稀疏非负矩阵分解的高光谱波段选择方法，其特征在于，步骤(6b)中所述的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子的更新公式如下：

$S^{(t+1)}=S^{\left(t\right)}=\frac{X^T{P}^{\left(t\right)}+{\mathrm{\βW}}^S{S}^{\left(t\right)}}{S{\left(P^{\left(t\right)}\right)}^T{P}^{\left(t\right)}+{\mathrm{\βD}}^S{S}^{\left(t\right)}}$

其中，S^(t+1)表示第t+1次迭代下更新的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，S^(t)表示第t次迭代下的高光谱图像数据的特征空间的含数据空间局部几何信息的非负矩阵分解因子，P^(t)表示第t次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，X表示高光谱图像数据，W^S表示高光谱图像数据的特征空间的相似度矩阵，D^S表示高光谱图像数据的特征空间的相似度对角矩阵，^T表示转置操作，β为特征空间的正则参数，取值范围为{10^-3,10^-2,10^-1,10⁰,10¹,10²,10³}。

7.根据权利要求1所述的基于双图稀疏非负矩阵分解的高光谱波段选择方法，其特征在于，步骤(6c)中所述的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵的更新公式如下：

$U^{(t+1)}=diag\left(\frac{1}{2\sqrt{\underset{i=1}{\overset{k}{\Σ}}{\[P^{(t+1)}*P^{(t+1)}\]}_i}}\right)$

其中，U^(t+1)表示第t+1次迭代下更新的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子的对角矩阵，P^(t+1)表示第t+1次迭代下的高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，*表示Hadamard矩阵乘积操作，[P^(t+1)*P^(t+1)]_i表示[P^(t+1)*P^(t+1)]的第i列，∑表示叠加操作，表示开方操作，diag(.)表示生成对角矩阵操作，k表示所选取的高光谱图像数据的波段数目。

8.根据权利要求1所述的基于双图稀疏非负矩阵分解的高光谱波段选择方法，其特征在于，步骤(9a)中所述的波段评价向量公式如下：

$G=\underset{i=1}{\overset{k}{\Σ}}{\[P*P\]}_i$

其中，G表示高光谱图像数据的波段评价向量，P表示高光谱图像数据的数据空间的含特征空间局部几何信息的非负矩阵分解因子，*表示Hadamard矩阵乘积操作，[P*P]_i表示[P*P]的第i列，k表示所选取的高光谱图像数据的波段数目，∑表示叠加操作。