CN106886793B - 基于判别信息和流形信息的高光谱图像波段选择方法 - Google Patents

Abstract

本发明提出了一个基于判别信息和流形信息的高光谱图像波段选择方法，主要解决现有波段选择方法未能有效利用数据流形信息，丢失数据原始物理意义及存在误差积累的问题。其实现步骤是：1)把高光谱图像转化为二维矩阵Y；2)对Y做归一化处理得到二维矩阵X；3)根据X分别构建核函数K和线性判别表达式；4)根据3)的结果计算图正则矩阵G；5)结合G迭代计算系数表示矩阵W；6)判断迭代次数k是否大于等于最大迭代次数，若是，输出最终的W，执行7)，否则，令k＝k+1，返回5)；7)根据W计算所有波段评分构成新的矩阵X^*，完成波段选择。本发明减小了累积误差，选出的波段更有代表性，可用于高光谱图像的预处理。

Description

基于判别信息和流形信息的高光谱图像波段选择方法

技术领域

本发明属于图像处理领域，具体涉及的是一种对高光谱图像进行波段选择的方法，可用于对农业、地质、大气和水文等领域中涉及的高光谱图像的预处理。

背景技术

受益于高光谱成像技术的飞速发展和成像质量的不断提升，高光谱图像在农业、地质、大气和水文等领域得到了越来越广泛的应用，图像中所包含的信息也越来越丰富。另一方面，丰富的图像信息也带来了数据量大、波段冗余等问题，制约了高光谱图像的存储、传输和处理。目前，国内外已经提出了不少高光谱波段选择方法来解决上述问题。

V.Kumar等人在其发表的论文“Band selection for hyperspectral imagesbased on self-tuning spectral clustering”(EUSIPCO,2013:1-50)中提出了一种基于自调整谱聚类的技术用于高光谱的波段选择。该方法可以大致总结为两个步骤：第一步，使用谱聚类方法把相似的波段归为一类，总共得到若干个类簇；第二步，利用主成分分析法分别处理每个类簇，计算得到每个类簇对应的特征向量基，该特征向量基就代表了对应类簇的低维表示。当合并所有类簇的特征向量基后，就得到了原始数据集低维下的表示矩阵，从而达到降维的目的。此方法的一个优点是算法简洁易于实现，能够较好去除冗余；另一个优点是分块处理不同类簇，易于并行计算，提高运行效率。但是此方法的缺点是第二步的主成分分析受到第一步谱聚类的影响，计算误差会累积，对最后的波段选择结果造成了不利影响；此外，由于第二步的主成分分析生成了全新的特征，丢失了原始数据的物理意义，可解释性变差。

西安电子科技大学在其申请的专利“基于低秩表示的高光谱图像波段选择方法”(专利申请号：CN201510411250.9，公开号：CN105046276A)中提出了一种基于低秩表示技术的高光谱图像波段选择方法。该方法首先对高光谱波段进行低秩表示，求解得到低秩表示系数，然后根据低秩表示系数矩阵对波段进行聚类得到若干个类簇，最后从每个类簇中选出离聚类中心最近的波段作为该类簇最具代表性的波段。最终得到的降维后的图像就是各类簇中代表性波段合并后的结果。该方法的优点是充分利用了低秩表示的优点增强了对噪声的抗干扰能力和对冗余的去除能力；再者，基于聚类技术的优势，该方法充分使用了波段特征的相关性，提高了后续分类的精度。但是，该方法的不足之处在于未能利用数据内在的流形信息，导致其学习到的低秩表示系数不够准确，影响对波段代表性的评判，即可能错误地判定了不同波段的重要性大小，导致最终选出的波段组合成的新图像的质量得不到保证。

发明内容

本发明的目的在于针对上述现有技术的不足，提出一种基于判别信息和流形信息的高光谱图像波段选择方法，以加强数据的可解释性，消除分步计算的累积误差，提高低秩表示系数的计算准确性，最终选择出最有代表性的波段组合成新的高质量图像。

本发明的技术方案是：把线性判别分析与核技巧结合起来，把数据的全局判别信息保留到判别矩阵中；再把该矩阵统一到图嵌入的框架内，计算得到一个包含判别信息的图正则矩阵，该矩阵同时保留了全局数据流形信息和判别信息；然后，通过回归学习的方法迭代更新得到原始数据的低维表示系数矩阵，利用该系数矩阵的范数约束得到波段特征的评分结果，从原始数据中选择评分高的波段作为代表性波段，组合成新的高光谱图像数据矩阵。其具体步骤包括如下：

(1)输入高光谱图像I并将其转化为二维数据矩阵Y；

(2)对二维数据矩阵Y做归一化处理，得到原始空间中归一化的二维数据矩阵X；

(3)使用映射函数φ把原始空间中的二维数据矩阵X表示为高维线性可分空间中对应的形式，并使用核技巧方法构建能表达该形式的核函数K；

(4)计算高维空间下数据的总散度矩阵

和类内散度矩阵

并将其散度矩阵代入标准线性判别模型中，得到高维空间下的线性判别表达式

其中：μ是不等于0的判别参数，I_n是单位矩阵，n是数据集中的样本数量，tr(·)是矩阵的秩，(·)^-1是矩阵的逆；

(5)计算一个同时包含数据判别信息和流形信息的图正则矩阵G：

(5a)通过一个常数表达式tr(H^TC_nH)对步骤(4)得到的线性判别表达式做等价变换，得到线性判别表达式的初次等价变换式：

其中：C_n是置中心矩阵，H是伪标签矩阵，矩阵(·)^T是(·)的转置矩阵；

(5b)使用核函数K对初次等价变换式再进一步等价变换，得到线性判别表达式的最终等价变换式：tr(H^T(C_n-C_n ^T(C_n+μK^-1)^-1C_n)H)；

(5c)根据流形学习理论，计算出等价变换后的线性判别表达式中同时带有判别信息和流形信息的图正则矩阵：G＝C_n-C_n ^T(C_n+μK^-1)^-1C_n；

(6)利用图正则矩阵G,使用更新迭代的方法计算系数表示矩阵W：

(6a)初始化系数表示矩阵W和伪标签矩阵H为元素值符合高斯分布的随机矩阵，初始化与W有关联的中间变量U为单位矩阵；

(6b)设置最大更新迭代次数为150，并设置当前迭代次数为k＝0；

(6c)利用中间变量U的更新公式

计算当前迭代k下的中间变量矩阵U^(k),其中：

是W的行向量的l₂范数，u_ii是U中第i行第i列元素；

(6d)利用伪标签矩阵变量H的更新公式计算当前迭代下的伪标签矩阵H^(k)，其中：h_ij是H中第i行第j列元素，G⁺＝(|G|+G)/2，G^-＝(|G|-G)/2，λ称为自由参数，取值范围是10⁸～10¹²，β称为约束项平衡参数，取值范围是10^-2～10³；

(6e)利用系数表示矩阵变量W的更新公式

计算当前迭代下的系数表示矩阵W^(k)，其中：w_ij是W中第i行第j列元素，α称为回归项平衡参数，取值范围是0.1～3；

(7)判断当前迭代次数k是否大于等于最大迭代次数，若是，输出最终的系数表示矩阵W，执行(8)，否则，令k＝k+1，返回步骤(6c)；

(8)完成整个波段选择过程：

(8a)计算每个波段的评分结果p_i；

(8b)计算所有波段的评分集合||W||_2,1；

(8c)选择评分靠前的一批波段组合成新的数据矩阵X^*，完成波段选择。

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

第一，本发明使用回归学习框架，通过更新迭代的方法得到关键的表示系数矩阵，克服了现有技术缺乏统一框架导致误差递进累积的缺点，能选出更具有代表性的波段。

第二，本发明利用核判别分析和图谱理论，计算得到了同时包含数据判别信息和流形信息的图正则矩阵，克服了现有技术未能充分利用数据全局流形信息的问题，提高了波段选择的准确性。

附图说明

图1是本发明的实现流程图；

图2是本发明与现有波段选择方法的分类精度对比图；

图3是本发明与现有波段选择方法的分类成像效果对比图。

具体实施方式

下面结合附图对本发明实施例和效果做进一步的描述。

参照图1，本发明的实现步骤依次包括：高光谱数据的二维归一化处理、构建核函数、构建高维空间下的线性判别表达式、计算图正则矩阵、计算系数表示矩阵、计算波段评分、生成新的高光谱图像。下面对这些步骤做具体介绍：

步骤1，输入高光谱图像并转化为二维数据矩阵。

在本发明实施例中，输入的高光谱图像为经典的Indian Pines。该图像为一个三维矩阵：I∈R^p×q×m，其中p×q表示该图像的像素点数目，m表示波段数目；为了便于后续处理，把p×q个像素点的矩阵重新排列成n行的像素点向量，从而把三维矩阵转化为二维矩阵Y∈R^n×m，这里n＝p×q表示二维数据矩阵的样本数目，R表示实数域；该高光谱图像的每一个像素点代表一个样本，整幅图像包括c个类别。

步骤2，归一化二维数据矩阵。

如前所述，二维矩阵Y∈R^n×m表示该高光谱图像包含m个样本和n个波段，将该二维矩阵按列归一化，也就是分别把每列中的元素归一化为0-1之间的值，得到归一化后的二维数据矩阵X∈R^n×m，避免数值较小的波段特征在后续处理中被忽略。

步骤3，构建核函数。

由于原始数据可能分布在一个非线性的流形中，而后续方法只能处理线性数据，需要把非线性数据线性化。通常的做法是先求数据的局部线性流形，再把所有局部线性流形做积分得到数据的整体流形。除该做法外，基于“非线性数据在其高维映射中线性可分”这个理论，也可以使用核技巧方法建立原始数据的波段特征空间与其高维映射间的联系来计算非线性流形。考虑到核技巧方法的准确度更高，本发明中使用并不限于用核技巧方法处理非线性流形，其处理步骤如下：

3a)引入映射函数φ(·)，把二维数据X从当前实数域空间R映射到高维线性可分空间Г，用φ(X)表示，φ(X)中的任一样本可表示为φ(x)；

3b)根据核技巧方法，本发明中选取高斯函数K作为核函数，则原始数据X与其对应的高维空间中的映射φ(X)间的关系可以表示为：

其中：x_i和x_j分别表示数据集中的第i个样本和第j个样本，||·||²是范数的平方，exp(·)表示核技巧方法。

3c)依据莫塞尔定理，原始空间中两个样本的高斯形式等同于其高维映射空间中的内积，即φ(x_i)^T·φ(x_j)，矩阵(·)^T是(·)的转置矩阵，该定理的数学表示为：

K(x_i，x_j)＝φ(x_i)^T·φ(x_j)，

3d)结合上述两个公式，可建立低维空间和高维空间中的数据关系为：

通过该数据关系实现把非线性数据转化为线性数据。

步骤4，构建高维线性可分空间Г中的线性判别表达式。

4a)定义高维空间下对应的总散度矩阵

和类间散度矩阵

的求解公式：

由于标准线性判别表达式的建立是依赖于总散度矩阵S_t和类间散度矩阵S_b，且使用的线性判别表达式要在高维空间Г下计算，因此本实例分别把S_t和S_b转换为其高维空间下对应的形式

和

其求解公式分别如下：

其中，C_n＝I_n-(1/n)1_n1_n ^T是置中心矩阵，其目的是使φ(X)中的同类别数据尽可能集中；I_n是n行n列的单位矩阵；1_n是n维的值全为1的向量；H是被松弛为连续形式的伪标签矩阵；

4b)结合

和

的求解公式，构建出高维空间Г中的线性判别表达式为：

其中：μ为不等于0的判别参数。

步骤5，计算同时包含判别信息和流形信息的图正则矩阵G。

5a)利用tr(H^TC_nH)＝c-1这个常数表达式，对步骤4得到的线性判别表达式做初次等价变换，得到初次等价的线性判别表达式：

5b)使用步骤3b)得到的核函数K对初次等价变换式再做进一步等价变换，得到最终的等价变换式：

5c)令上述最终等价变换式中的C_n-C_n ^T(C_n+μK^-1)^-1C_n＝G，基于“高维空间中距离相近的样本在其对应的低维空间中距离仍然相近”的流形学习思想，得出G中必然同时包含了判别信息和流形信息结论，故把G命名为图正则矩阵。

步骤6，计算系数表示矩阵W。

本发明采用交替更新迭代的方法计算系数表示矩阵，计算过程中包括三个相互联系的矩阵变量，分别是：系数表示矩阵W∈R^m×c、伪标签矩阵H∈R^n×c和中间变量矩阵U∈R^m×m，其计算步骤如下：

6a)在交替更新迭代之前先初始化上述三个矩阵变量：即利用MATLAB程序中的随机生成函数rand()，分别设置初代的系数表示矩阵W和伪标签矩阵H为元素值符合0～1间高斯分布的随机矩阵；利用MATLAB程序的单位矩阵函数eyes()，设置初代的中间变量U为m行m列的单位矩阵；

6b)设置更新迭代次数：依据经验，设置最大更新迭代次数为150次，同时，设置当前迭代次数为k＝1；

6c)利用中间变量矩阵U的更新迭代公式，按元素计算当前迭代次数k下中间变量矩阵U中的第i行第i列元素

其中，

表示系数表示矩阵W中第i个行向量的l₂范数；

6d)利用伪标签矩阵变量H的更新迭代公式，按元素计算当前迭代次数k下伪标签矩阵H中第i行第j列元素h^(k) _ij：

由于图像数据矩阵中的元素值都是正数，因此这里把图正则矩阵G分解为两个正矩阵G⁺＝(|G|+G)/2和G^-＝(|G|-G)/2，以保证更新公式的非负性，式中β为约束项平衡参数，λ称为自由参数，经验设置β的调参区间为[10^-2,10^-1,1,10¹,10²,10³]，λ的调参区间为[10⁸,10⁹,10¹⁰,10¹¹,10¹²]。

6e)利用系数表示矩阵W的更新迭代公式，按元素更新计算当前迭代次数k下新的系数表示矩阵W中第i行第j列元素w^(k) _ij：

其中，α为回归项平衡参数，根据经验设置α的调参区间为[0.1～3]。

步骤7，判断当前迭代次数k是否大于等于最大迭代次数，若是，输出最终的系数表示矩阵W，执行(8)，否则，令k＝k+1，返回步骤(6c)。

步骤8，完成整个波段选择过程。

波段选择是通过对系数表示矩阵W施加范数约束得到的。通常，范数约束方法包括l₁范数，l₂范数和l_2,1范数。l₁范数只能对波段特征做二元评价，即要么判定为重要，要么判定为不重要，不能给出各个波段特征的相对重要性。l₂范数克服了l₁范数的二元评价缺陷，对各个波段特征的相对重要性给出了准确数值，但是其选择能力弱，会遗漏一些重要特征。l_2,1范数结合了l₁范数和l₂范数各自的优点，既能给出各个波段特征的相对重要性，又不会遗漏重要特征，因此本实例中使用但不限于用l_2,1范数计算各个特征的重要性评分，计算步骤如下：

8a)通过系数表示矩阵W的行向量

的元素之和来度量对应的波段特征的重要性，计算对应波段的评分结果p_i：

其中：w_ij表示系数表示矩阵W第i行第j列元素，c表示二维数据矩阵X中的类别个数，p_i越高，表示对应波段越重要；

8b)对所有波段的评分进行排序，得到所有波段评分从高到低排序的集合||W||_2,1：

8c)从评分排序集合中选择排序靠前的一批波段组合成新的数据矩阵X^*，即新的高光谱图像，具体选择的波段个数根据实际情况而定。

本发明的效果可以通过以下仿真实验进一步说明：

1、仿真条件

本发明的仿真实验在Inter Core i5(TM)2Duo CPU、主频2.50GHz，内存6GB，windows7 64位系统平台上的Matlab 2014a上完成；

仿真使用的数据集为经典的Indian Pines高光谱图像数据，该图像总共包括20251个像素点和200个波段；

仿真使用的分类器为KNN分类器，该分类器的设置方式为：随机选择上述数据集中的10％数据样本做KNN分类器的训练数据集，把剩余的90％样本做KNN分类器的测试数据集；设置KNN分类器的近邻参数p＝5；令KNN分类器使用十折交叉验证方法计算分类精度；

仿真使用的对比算法为自调整谱聚类波段选择方法(SC)。

2、仿真内容

仿真1，用本发明和SC方法在KNN分类器上以Indian Pines为输入数据，进行分类精度实验，分类精度结果如图2所示，其中直角坐标系，横坐标表示波段选择数目s，纵坐标表示分类精度OA，波段数目s的取值为：5、10、15、20、25、30、35、40、45、50，分类精度OA的取值在0～1之间。图2中3条折线中，圆圈标注的折线是本发明的分类精度，雪花形标注的折线是SC方法的分类精度，直线是未进行波段选择直接分类的分类精度，即基于当前分类器的最高分类精度ALL。

从图2可见，只有波段选择个数少于20个时，SC方法的精度才好于本发明，而这主要是因为当波段个数太少时生成的新数据集缺乏判别信息和流形信息，导致本发明所用技术得不到体现。随着波段选择个数的增加，本发明的分类精度不断提高，其表现如下：

首先，当波段个数大于20个后，本发明的分类精度超过SC方法，表明此时本发明的技术优势得到了体现。

此外，当波段选择个数超过20个后，本发明的分类精度提升速度开始放缓，这表明本发明选择的前20个波段是最具有代表性的波段。

再次，当选择的波段个数达到50个时，本发明的分类精度已经接近当前分类器下的最大分类精度，这表明Indian Pines数据集中其余未被选择的150个波段中有许多冗余波段，而本发明在选择波段时很好的避开了这些冗余波段。

最后，对比SC方法的折线图，发现当SC方法选择的波段个数大于35个时，SC方法的分类精度反而有所下降，这说明SC方法选出的波段中存在一些噪声波段，而这一现象并未出现在本发明的折线结果中，表明本发明选择波段时不仅避开了冗余波段，而且避开了噪声波段。

仿真2，用本发明与SC方法在KNN分类器以Indian Pines为输入数据，进行分类成像效果实验，成像效果如图3所示，其中图3a是Indian Pines高光谱图像的原始真值图；图3b是KNN分类器下的最好成像效果图，即未进行波段选择得到的分类成像效果图；图3c是对SC方法选择50个波段生成的新图像进行KNN分类后得到的成像效果图；图3d是本发明选择50个波段进行KNN分类后得到的成像效果图。

由图3b和图3d可以看出，本发明的成像与KNN分类器的成像中的错分区域比较一致，这说明很多错分点来自KNN分类器本身；

由图3c和图3d可以看出，相比SC方法本发明的错分点更少，表明本发明比SC方法选出的波段更有代表性，成像效果更好。

Claims (5)

1.一种基于判别信息和流形信息的高光谱图像波段选择方法，包括如下步骤：

(1)输入高光谱图像I并将其转化为二维数据矩阵Y；

(2)对二维数据矩阵Y做归一化处理，得到原始空间中归一化的二维数据矩阵X；

(3)使用映射函数φ把原始空间中的二维数据矩阵X表示为高维线性可分空间中对应的形式，并使用核技巧方法构建能表达该形式的核函数K；

(4)计算高维空间下数据的总散度矩阵

和类内散度矩阵并将其散度矩阵代入标准线性判别模型中，得到高维空间下的线性判别表达式

其中：μ是不等于0的判别参数，I_n是单位矩阵，n是数据集中的样本数量，tr(·)是矩阵的秩，(·)^-1是矩阵的逆；

(5)计算一个同时包含数据判别信息和流形信息的图正则矩阵G：

(5a)通过一个常数表达式tr(H^TC_nH)对步骤(4)得到的线性判别表达式做等价变换，得到线性判别表达式的初次等价变换式：

其中：C_n是置中心矩阵，H是伪标签矩阵，矩阵(·)^T是(·)的转置矩阵；

(5b)使用核函数K对初次等价变换式再进一步等价变换，得到线性判别表达式的最终等价变换式：tr(H^T(C_n-C_n ^T(C_n+μK^-1)^-1C_n)H)；

(5c)根据流形学习理论，计算出等价变换后的线性判别表达式中同时带有判别信息和流形信息的图正则矩阵：G＝C_n-C_n ^T(C_n+μK^-1)^-1C_n；

(6)利用图正则矩阵G,使用更新迭代的方法计算系数表示矩阵W：

(6a)初始化系数表示矩阵W和伪标签矩阵H为元素值符合高斯分布的随机矩阵，初始化与W有关联的中间变量U为单位矩阵；

(6b)设置最大更新迭代次数为150，并设置当前迭代次数为k＝0；

(6c)利用中间变量U的更新公式

计算当前迭代k下的中间变量矩阵U^(k),其中：是W的行向量的l₂范数，u_ii是U中第i行第i列元素；

(6d)利用伪标签矩阵变量H的更新公式

计算当前迭代下的伪标签矩阵H^(k)，其中：h_ij是H中第i行第j列元素，G⁺＝(|G|+G)/2，G^-＝(|G|-G)/2，λ称为自由参数，取值范围是10⁸～10¹²，β称为约束项平衡参数，取值范围是10^-2～10³；

(6e)利用系数表示矩阵变量W的更新公式

计算当前迭代下的系数表示矩阵W^(k)，其中：w_ij是W中第i行第j列元素，α称为回归项平衡参数，取值范围是0.1～3；

(7)判断当前迭代次数k是否大于等于最大迭代次数，若是，输出最终的系数表示矩阵W，执行(8)，否则，令k＝k+1，返回步骤(6c)；

(8)完成整个波段选择过程：

(8a)计算每个波段的评分结果p_i；

通过W的行向量

的元素之和来度量对应的波段特征的重要性，计算对应波段的评分结果p_i：

其中：w_ij表示系数表示矩阵W第i行第j列元素，c表示二维数据矩阵X中的类别个数，p_i越高，表示对应波段越重要；

(8b)计算所有波段的评分集合||W||_2,1；

将所有波段的评分排序，得到所有波段评分排序的集合||W||_2,1：

(8c)选择评分靠前的一批波段组合成新的数据矩阵X^*，完成波段选择。

2.根据权利要求1所述的方法，其特征在于，步骤(1)中将输入的高光谱图像I转化为二维数据矩阵Y，是将输入的三维矩阵图像I∈R^p×q×m中的p×q个像素点矩阵重新排列成n行的像素点向量，从而把三维图像I转化为二维矩阵Y∈R^n×m，其中p×q表示该图像的像素点个数，m表示波段数目，n＝p×q表示二维数据矩阵的像素点数目，R表示实数域。

3.根据权利要求1所述的方法，其特征在于，步骤(2)对二维数据矩阵Y做归一化处理，是将二维数据矩阵Y按列归一化，即分别把每列中的元素归一化为0-1之间的值，得到归一化后的二维数据矩阵X∈R^n×m。

4.根据权利要求1所述的方法，其特征在于，步骤(3)中构建的核函数K，表示如下：

其中，exp(·)表示核技巧方法，x_i和x_j是二维数据矩阵X中的任意两个样本，||·||²是范数的平方，φ(x_i)和φ(x_j)分别表述x_i和x_j在高维空间中的映射，φ(x_i)^T·φ(x_j)是φ(x_i)与φ(x_j)的内积。

5.根据权利要求1所述的方法，其特征在于，步骤(4)中的总散度矩阵

和类内散度矩阵

分别表示如下：

其中：C_n为置中心矩阵，C_n＝I_n-(1/n)1_n1_n ^T，1_n表示元素值都是1的n维列向量，φ(X)是二维数据矩阵X在高维空间中的映射。

Images