CN102289673A - 基于偏最小二乘法的高光谱遥感图像波段选择方法 - Google Patents

Abstract

本发明属于高光谱遥感图像处理技术领域，具体为一种基于偏最小二乘法的高光谱遥感图像波段选择的方法。本发明利用偏最小二乘法提取成分保留高光谱图像变异信息且与分类信息相关程度高的特点，将光谱矩阵与隶属度矩阵乘积的能量作为选择波段的标准，通过迭代求取已选择波段递归的残差来选择下一组波段，实现波段选择的过程。本发明能有效克服传统多光谱图像波段选择方法计算复杂度高，需要去除相关波段的缺点。高光谱遥感图像分类实验结果表明，高光谱遥感图像使用本发明进行波段后具有良好的分类效果。本发明对于高效利用高光谱图像的信息资源有着重要的应用价值。

Description

基于偏最小二乘法的高光谱遥感图像波段选择方法

技术领域

本发明属于遥感图像处理技术领域，具体涉及一种基于偏最小二乘法的高光谱遥感图像波段选择方法。

背景技术

高光谱遥感在电磁波谱的紫外、可见光、近红外和中红外区域中的数十至数百个非常窄且连续的光谱段上同时获取信息，从而得到每个像素完整连续的光谱曲线^[1]。高光谱遥感技术的发展为人们提供了详尽精确的地物光谱信息，但是大量高度相关的波段信息也为高光谱图像的进一步处理带来了困难。高维数据不仅会造成冗余，提高计算复杂度，而且由于Hughes现象^[2]（即训练样本固定时，分类精度随特征维数的增加而下降），使得高光谱图像的分类效果受到严重影响。因此，如何有效选取特征对于高光谱图像处理是非常重要的。一般来说，选择有效特征主要采取特征提取或波段选择这两种方法。

相对于特征提取的方法，波段选择不改变原光谱数据，保留原始波段的物理含义以及地物的光谱特性，在相关领域有着重要应用价值。当前对高光谱图像波段选择方法的研究，主要针对其在分类中的应用，具体可分为有监督和无监督两种。两者的区别在于，无监督的方法不需要任何先验信息，而有监督的方法需要事先为每类地物在遥感图像上采集样本数据之后，通过学习训练再进行波段选择^[3]。相对而言，有监督的波段选择方法往往能够取得更高的分类精度，而且在与有监督的分类方法结合使用时，不需要专门为波段选择寻找训练数据。本发明是一种用于分类的高光谱图像有监督波段选择方法。

用于分类的有监督传统波段选择的方法一般有以下两种思路：一、基于特征子集搜索的波段选择，将高光谱数据中选择的波段作为特征子集，以某种准则进行全局或局部搜索。常见的准则有：各种距离（巴式距离、JM距离），信息论的指标（散度、变换散度）等^[4]。由于最优的穷举法搜索运算量极大，实际中常采用次优的特征子集搜索算法（如前向搜索法）得到有效波段。这类方法物理意义明确，且易于实施，但计算相对较为耗时。二、基于特征值分解的一系列方法，以最小误分类典型分析为例，该方法通过求解费舍尔分离矩阵的特征值和特征向量来计算各波段对分离矩阵的能量贡献，再经过对能量贡献排序和去除相关波段得到有效波段^[5]。

偏最小二乘法作为一种新型的多元统计数据分析方法，是一种多因变量对多自变量的回归建模方法，可以实现多种数据分析方法的综合应用，尤其适合于多重线性相关问题，在分析光谱数据中有着重要价值^[6]。在高光谱遥感图像分类的研究中，偏最小二乘法提取的潜在向量不但与分类信息的相关度高，而且尽可能多地保留高光谱数据的变异信息。因此，选择偏最小二乘法通过特征值分解计算各个波段能量贡献在保留有效信息方面有其优势。此外，偏最小二乘法中迭代求取潜在向量回归残差的步骤可以有效避免选择波段冗余的情况。

发明内容

本发明的目的在于提供一种计算复杂度低、算法简单而有效的高光谱遥感图像波段选择方法。

本发明提出的高光谱遥感图像像波段选择方法，是一种基于特征值分解，进行波段排序的波段选择方法。与传统方法直接对所有波段排序不同，本发明通过迭代求取已选择波段递归的残差，再通过排序得到当前这迭代最优波段，以此逐个得到所选择的波段。该方法将偏最小二乘法目标函数中各个波段能量贡献作为选择波段的标准，通过迭代求取已选择波段递归的残差逐个得到选择的波段。此外，本发明由于避免了特征子集搜索和大矩阵特征值分解的运算，具有更快的运算速度，为高光谱遥感图像波段选择提供了一种新的有效的快速算法。

先介绍有关偏最小二乘法的概念：

偏最小二乘法最早是由Wold在1983年提出的一系列用于线性回归、特征提取的算法。偏最小二乘法不仅关注自变量的信息，而且关注自变量与因变量之间的相关关系，因此偏最小二乘法所提取的特征信息对于回归分析而言往往更为有用。偏最小二乘法之所以能够做到这一点，与它的目标函数选择有关。偏最小二乘法的核心思想是要找到一组自变量 的潜在向量和因变量的潜在向量，要求、在保留数据信息的同时，能够很好地解释。、能分别很好代表与中的数据变异信息，根据主成分分析原理，应该要求、的方差尽可能大；另一方面，由于回归建模的需要，又要求对有很强的解释能力，有典型相关分析的思路，与的相关度应达到最大值。综合起来，偏最小二乘法要求与的协方差达到最大。作如下定义，、分别是和的潜在向量，、是对应的单位轴，偏最小二乘法的目标函数如下：

(1)

式(1)中分别为自变量和因变量的投影轴，为自变量和因变量的协方差。如果采用拉格朗日算法，和分别是去均值后的自变量和因变量，则有

(2)

求解目标函数取最大值，可以归结为以下特征值问题：

(3)

由式(3)可知，

是矩阵

的特征向量，对应的特征值为

。

为目标函数值，要求取最大值。所以，

是对应于

矩阵最大特征值的单位特征向量。其中，具有上标T的向量或矩阵表示该向量或矩阵的转置，下同。

根据和可得到成分,。偏最小二乘回归分别实施对的回归以及对 的回归。如果回归方程已经达到满意的精度，则算法终止；否则,将利用被解释后的残余信息以及被解释后的残余信息进行第二轮的成分提取：

	(4)
			(5)

如此往复，通过迭代求取潜在向量回归的残存矩阵就能够求出所需要的潜在向量组。在实际计算中，对

求残差的过程可以被省略。可以证明的，直接用

进行下一次迭代的计算结果是完全等价的。

本发明提供的高光谱遥感图像波段选择方法基于偏最小二乘法的高光谱遥感图像波段选择方法描述如下。

基于偏最小二乘法的高光谱遥感图像波段选择

由于偏最小二乘法在保留自变量（遥感图像）变异信息的同时反映其与因变量（地物分类信息）的联系，因此基于偏最小二乘法的波段选择与单独考虑信息量或类间分离性的方法相比，在保留有效信息方面更有优势。

将训练数据集作为自变量，隶属度矩阵作为因变量，由偏最小二乘法计算其潜在向量的投影轴。偏最小二乘法得到的潜在向量保留了大量有利于分类的信息，分析利用潜在向量可以帮助有效波段的选择。由偏最小二乘法的性质可知，潜在向量的投影轴是样本空间的一组标准正交基。根据文献[5]中载荷因子的定义，可得：

(6)

式(6)中，

是投影权重向量，

是第个投影轴的第

个元素，

是其权重，

是投影轴个数，

代表了第

个波段的贡献。在文献[5]中，Chang提出了一系列基于特征值分解的波段选择算法，以特定矩阵(如协方差矩阵，费舍尔判别矩阵)的特征值作为对应投影轴的权重，根据式(6)判断各个波段的重要程度来进行波段选择。该方法利用各个波段的重要程度进行排序，只考虑单个波段的作用，而不考虑所选波段集合中波段之间的关系，最后仅仅简单地去除候选集合中高度相关的波段，未必能够使得所选波段达到最优。同时由于偏最小二乘法所求得的每一个潜在向量投影轴并非由同一个矩阵特征值分解得到，而是通过迭代得到残差矩阵最大特征值所对应的特征向量，因此基于偏最小二乘法的波段选择中

的构造不能直接像文献[5]一样取对应特征向量的特征值。

考察偏最小二乘法中迭代求取残差的过程，在偏最小二乘法中，根据当前自变量计算一个潜在向量，再由这个潜在向量计算当前自变量被它解释后的残余信息，然后由这个自变量的残差计算下一个潜在向量，直到得到足够多的潜在向量。在波段选择中也可以利用这一求残差过程尽可能剔除已选择波段的影响，将偏最小二乘法中迭代求取潜在向量的过程与波段选择的过程相统一。将当前所选择的波段归一化后作为当前的潜在向量,计算自变量被它解释后的残差，再由此选择下一个波段。

此时在每一次迭代中，就可以使用式(6)的构造,

取当前计算得到的所有特征值。将基于特征值分解的波段选择算法和偏最小二乘法中的迭代算法结合，得到一个最优的波段，再计算已选择波段回归的残差，根据残差选择下一个波段，如此反复，直到得到所有需要的波段。

最后，将本发明中的基于基于偏最小二乘法的高光谱遥感图像波段选择的算法流程归结如下：

步骤 1、将高光谱数据转换为矩阵形式，根据训练数据的地物信息构造隶属度矩阵

，其中N为训练样本数，L为波段数，M为地物类别数；

步骤 2、计算乘积矩阵的所有特征向量，特征值，并由式（7）计算投影权重向量：

,

(7)

根据式(6)定义，计算载荷因子：

,

(8)

其中

为

的第个元素；

步骤 3、从L个载荷因子中选择最大值，将

记为选择的波段；

步骤 4、根据已选择波段对做递归并求残差：

(9)

其中

为

的第

个列向量；

步骤 5、当所选择的波段数已满足要求，则结束循环；否则，返回步骤 2，继续迭代。

本发明的优点

本发明为一种新的基于偏最小二乘法的高光谱遥感图像波段选择方法。本发明利用偏最小二乘法提取成分保留高光谱图像变异信息且与分类信息相关程度高的特点，将光谱矩阵与隶属度矩阵乘积的能量作为选择波段的标准，通过迭代求取已选择波段递归的残差来选择下一组波段，实现波段选择的过程。本发明能有效克服传统多光谱图像波段选择方法计算复杂度高，需要去除相关波段的缺点。与传统的高光谱波段选择方法相比，采用偏最小二乘法的目标函数计算各波段能量贡献更有效地保留了有用信息，通过迭代求已选择波段回归的残差，再由此逐个选取波段，避免了冗余波段的重复选择。同时，由于避免了特征子集搜索和大矩阵特征值分解的运算，新算法具有更快的运算速度。高光谱数据实验表明，与传统方法相比，在选取相同波段数的情况下可取得较高的分类精度，证实了算法的有效性和对于不同数据的适用性。此外由于避免了特征子集搜索和大矩阵特征值分解的运算，新算法具有更快的运算速度。

附图说明

图1. Indian Pine 高光谱遥感数据。其中，(a)伪彩色图，(b)地面真实。

图2. Washington DC 广场高光谱遥感数据。

图3. Indian Pine AVIRIS 高光谱遥感数据六种波段选择方法平均分类精度比较。其中，(a)(b)是采用KNN进行分类，(c)(d)是采用SVM进行分类，(e)(f)是采用LDA进行分类，其中(a)(c)(e)是训练集与测试集取1:2进行分类，(b)(d)(f)则是训练集与测试集取1:4进行分类。

图4. Washington DC HYDICE 高光谱遥感数据六种波段选择方法平均分类精度比较。其中，(a)(b)是采用KNN进行分类，(c)(d)是采用SVM进行分类，(e)(f)是采用LDA进行分类，其中(a)(c)(e)是训练集与测试集取1:2进行分类，(b)(d)(f)则是训练集与测试集取1:8进行分类。

具体实施方式

在高光谱遥感数据实验中，本方法与以下五种波段选择方法进行比较：JM方法以任意两类之间的JM距离作为评判指标进行特征子集搜索，实验采用运算量较小且在波段选择中表现较好的前向搜索算法SFS进行子集搜索^[4]。最大方差主成分分析(Maximum Variance Principal Components Analysis, MVPCA)和最小误分类典型分析(Minimum Misclassification Canonical Analysis, MMCA)是Chang提出的基于特征值分解的方法^[5]，这两种方法分别以线性判别矩阵和协方差矩阵特征值分解求出各个波段的重要性，进行排序并去除相关波段。最小估计丰度方差(Minimum Estimated Abundance Covariance, MEAC)的方法以丰度估计矩阵的协方差最小作为选择波段集合的指标进行特征子集搜索^[7]。特征权重(Feature Weighting，FW)的方法以散度作为各类主成分分析的权重确定各个波段在各类别中的重要性排序，再通过投票确定各个波段的重要性，进行排序并去除相关波段^[8]。

通过比较各种波段选择方法所选择的波段在最近邻(K-Nearest Neighbour, KNN)^[9]、线性判别分析(Linear Discriminant Analysis, LDA)^[10]、支撑向量机(Support Vector Machine，SVM)^[11]三种分类器下的分类精度，以说明本方法的适用性。

下面，分别用两组实际遥感图像数据为例说明具体的实施方式。

实验例 1. AVIRIS 高光谱遥感数据

表1. Indian Pine地物真实对照

类别	对应地物	分类样本个数
			C1	Corn-notill	1256
C2	Corn-min	726
			C3	Grass/Pasture	431
C4	Grass/Tress	626
			C5	Hay-windrowed	443
C6	Soybeans-notil	828
			C7	Soybeans-min	2284
C8	Soybeans-clean	503
			C9	Woods	1198

实验使用1992年由机载可见光及红外成像光谱仪(Airborne Visible/Infrared Imaging Spectrometer，AVIRIS)采集的美国印第安纳州Pine测试点高光谱影像（https://engineering.purdue.edu/~biehl/MultiSpec/提供该数据的免费下载）。影像块大小为145×145个像素，总共220个光谱波段。波长范围从0.4~2.5μm，光谱分辨率为10nm，空间分辨率为17m。该数据已被广泛地用于遥感图像的分类研究。

第1-4，103-113以及148-166波段由于信噪比太低或为水吸收波段而被移除，剩下186个波段。取第70、86、136波段分别作为R、G、B分量合成伪彩色图，如图1(a)所示。选择{Corn-notill, Corn-min, Grass/Pasture, Grass/Tress, Hay-windrowed, Soybeans-notill, Soybeans-min, Soybeans-clean, Woods}这9种主要地物进行分类，其余区域作为空白区域不用于分类。其对应的地物真实如图1(b)所示。这9个类各自所对应的地物和其分类样本个数如表1所示。

由于Indian Pine测试点的高光谱图像数据存在混合像元，几种相近的农作物之间易混淆，因此分类的难度较大。为保证分类效果，考虑给予充分的训练。使用KNN、LDA和SVM分类，分别以训练集与测试集取1:2和1:4进行分类实验，结果如图3所示。其中横坐标均为波段个数，而纵坐标为分类精度。图中的虚直线表示的是使用全波段所得到的分类精度。

下面我们讨论各种分类方法的情况。KNN在有充足训练的情况下，分类效果优于SVM，但是各种波段选择方法在不同选择波段数下的平均分类精度趋势大致相近。由于LDA本身就要进行降维，因此进行波段选择的效果并不明显。但是正因为LDA分类方法中有特征提取的步骤，从LDA的分类精度可以看出各种波段选择方法丢失可用信息的情况。

不同的分类方法和训练集的大小对分类精度的高低有直接的影响，但对各种波段选择方法的分类精度随所选波段数变化趋势的影响并不明显。

比较各种波段选择方法，在训练数据足够多时，JM表现稳定，随着选择波段数的增加有着不错的分类效果；MMCA的分类效果略逊于JM；MEAC与FW的表现相近。MVPCA作为无监督的方法参与比较，主要是为了比较去除相关冗余波段对于波段选择是否有着决定性的影响，因此MVPCA的效果远逊于其他方法也从侧面说明了有监督的方法的必要性。本文的方法在选择波段较少时，所选择波段的分类精度要明显优于其它的波段选择方法。随着波段数增加，其分类效果也略优于其他方法。

实验例 2. HYDICE 高光谱遥感数据

实验使用由(Hyperspectral Digital Imagery Collection Experiment, HYDICE)传感器采集的Washington DC高光谱影像。影像块大小为1280×307个像素，总共210个光谱波段。该数据的地面分辨率高达2.8 m，混合像元少，各种地物的差别较为明显。

去除第103-106，138-148，207-210这些噪声和水吸收波段，剩余191个光谱波段。取第53、32、20波段分别作为R、G、B分量合成伪彩色图，如图2所示。

根据文献[12]中的地物信息选择Water, Grass, Shadow, Street, Path, Trees, Roof这7种地物进行分类，其余区域作为空白区域不用于分类，选取的各类样本数如表2所示。

表2. Washington DC地物真实对照

对应地物	分类样本个数
		Water	1224
Grass	175
		Shadow	97
Street	416
		Path	405
Tree	1928
		Roof	3834

Washington DC Mall数据由于图像本身分辨率高，地物中以纯像元为主，分类精度相对较高,因此相应地降低训练集的比例。使用KNN、LDA和SVM分类，分别以训练集与测试集取1:2和1:8进行分类实验，结果如图2所示。其中横坐标均为波段个数，而纵坐标为分类精度。图中的虚直线表示的是使用全波段所得到的分类精度。

由图4中可以看出，只要求保留光谱信息的MVPCA方法也取得了不错的分类效果，而强调类间距离的方法，如MMCA、JM的表现就不如在实验例1中那么突出。FW方法在这个数据的分类中有着不错的表现。本文的方法由于兼顾了分类信息和光谱变异信息而始终有着不错的表现。由此可以看出，本文提出的方法对于这一类地物类别差异明显的遥感图像同样有着良好的效果。

计算时间比较

表3比较了实验例1和实验例2中各种波段选择方法在选取15个波段时的运算时间。JM和MEAC这两种方法需要进行特征子集搜索，即使采用运算量相对较小的次优搜索算法运算，依然相对耗时，且运算耗时与选择波段数近似成正比。MMCA、MVPCA、FW以及本文提出的波段选择方法或多或少都需要做奇异值分解，但是由于本文的方法只需要对高光谱数据

和隶属度矩阵

的协方差阵做奇异值分解，而

每行的维数（即地物类别数目）要远少于高光谱的波段数。因此也就降低了运算量，减少了时间消耗。本文的方法既不需要做特征子集搜索，也不需要对大矩阵求特征值分解，其耗时要远少于JM距离等特征子集搜索的算法，且略优于同类的基于特征值分解的方法。

表3. 两次实验中选取15个波段时的耗时比较 （单位：秒）

参考文献：

[1] C.-I. Chang. Hyperspectral imaging: techniques for spectral detection and classification[M]. New York: Plenum, 2003.

[2] G.F. Hughes. On the mean accuracy of statistical pattern recognizers[J]. IEEE Trans. Inf. Theory, 1986, 14(1):56-63.

[3] S. Serpico and L. Bruzzone. A new search algorithm for feature selection in hyperspectral remote sensing images[J]. IEEE Trans. Geosci. Remote Sens., 2001, 39(7):1360-1367.

[4] J. A. Richards and X. Jia. Remote sensing digital image analysis: an introduction[M]. 4th ed. Berlin, Germany: Springer-Verlag, 2006.

[5] C.-I. Chang, Q. Du, T.-L. Sun, and M.L.G. Althouse. A joint band prioritization and band-decorrelation approach to band selection for hyperspectral image classification[J]. IEEE Trans. Geosci. Remote Sens., 1999, 37(6):2631-2641.

[6] R. Rosipal, N. Krämer. Overview and recent advances in partial least squares[J]. Lecture Notes in Computer Science, 2006, 3940: 34-51.

[7] H. Yang, Q. Du, H. Su, and Y. Sheng. An efficient method for supervised hyperspectral band selection[J]. IEEE Geosci. Remote Sens. Letters, 2011, 8(1):138-142.

[8] R. Huang and M. He. Band selection based on feature weighting for classification of hyperspectral data[J]. IEEE Geosci. Remote Sens. Letters, 2005, 2(2):156-159.

[9] T. Cover and P. Hart. Nearest neighbor pattern classification[J]. IEEE Trans. Inf. Theory, 1967, 13(1):21-27.

[10] Q. Du and C-I Chang. A linear constrained distance-based discriminant analysis for hyperspectral image classification[J]. Pattern Recognition, 2001, 34: 361-373.

[11] F. Melgani and L. Bruzzone. Classification of hyperspectral remote sensing images with support vector machines[J]. IEEE Trans. Geosci. Remote Sens., 2004, 42(8):1778-1790.

J.A. Benediktsson, J.A. Palmason and J.R. Sveinsson. Classification of hyperspectral data from urban areas based on extended morphological profiles[J]. IEEE Trans. Geosci. Remote Sens., 2005, 43(3):480-491。

Claims (1)

1.一种基于偏最小二乘法的高光谱遥感图像波段选择的方法，其特征在于通过迭代求取已选择波段递归的残差来选择波段，具体步骤如下：

步骤 1、将高光谱数据转换为矩阵形式

，根据训练数据的地物信息构造隶属度矩阵

，其中N为训练样本数，L为波段数，M为地物类别数；

步骤 2、计算乘积矩阵的所有特征向量

，特征值

，根据式（1）计算投影权重向量

：

,

(1)

带上标T的表示

的转置；根据式(2)计算载荷因子：

,

(2)

其中

为

的第

个元素；

步骤 3、从L个载荷因子中选择最大值

，将

记为选择的波段；

步骤 4、根据已选择波段对做递归并求残差：

(3)

其中

为的第

个列向量； 带上标T的

表示

向量的转置；

步骤 5、当所选择的波段数已满足要求，则结束循环；否则，返回步骤 2，继续迭代。

Images