基于特征优化的高光谱遥感影像降维算法

王小宁，宋伟东

(辽宁工程技术大学 测绘与地理科学学院，辽宁 阜新 123000)

0 引言

高光谱遥感处在当今遥感技术发展的前沿领域，相比于传统的遥感技术，高光谱影像具有更为丰富的地物信息，对于准确实现目标识别和地物分类等至关重要[1]。但Hughes现象严重影响高光谱数据的应用效率，采用数据降维的方式提取高光谱影像中有利于区分不同地物目标的信息是解决Hughes现象的重要手段。

主成分分析(principal component analysis，PCA)变换是一种最常用的降维算法[2]，通过在低维空间中找出几个互不相关的主成分来表示高维数据，实现数据的降维。该算法可以滤除影像波段间的二阶相关性，将有效信息集中于前几个主成分。然而高光谱影像各波段间相关性差异较大，对于相关性较弱的不同波段，基于全局统计特性的PCA变换降维效果不理想。此外，因该算法忽略波段间的重要性差异，不能有效处理那些对于特征值贡献较小，却对分类结果至关重要的波段，导致无法充分利用高光谱数据的差异性信息。为此，Chen等[3]提出一种改进算法，即加权主成分分析(weighted principal component analysis,WPCA)。其基本思想是根据不同波段特征在后续分类过程中所起作用的程度，定义相应权值，以达到在影像降维过程中保留重要信息的目的。

PCA变换对噪声和异常值较为敏感，若影像中包含过多噪声或异常值，则严重影响影像降维效果。为了克服以上缺陷，Green等[4]提出最小噪声分离(minimum noise fraction，MNF)变换。经MNF变换后的数据按照信噪比大小进行排列，在一定程度上能够克服噪声对影像质量的影响。然而高光谱影像混合像元现象较为明显，MNF变换仅基于空间邻域关系计算噪声分数，存在较大误差且计算结果不具稳定性，因而算法效果提升并不明显。Yang等[5]在PCA基础上提出二维主成分分析(complete two dimensional principal component analysis,2DPCA)，该算法基于二维影像构造协方差矩阵，并据此计算最优投影方向，获得影像特征表示。但本质上来说，2DPCA算法是一种局部化的PCA算法，其协方差阵所含信息量明显少于PCA，且该算法难以有效去除高光谱影像波段间的相关性。Scholkopf等[6]将核方法引入PCA，提出核主成分分析(kernel principal component analysis,KPCA)，将影像映射至高维核空间，并利用PCA变换提取数据主要信息。Zhang等[7]在此基础上提出超像素核主成分分析(SuperKPCA)，采用影像分割技术获得同质区域，再利用KPCA降低每个同质区域的维度，得到更为精确的高光谱影像的非线性低维特征，但核方法的本质是将低维空间映射至高维空间，增加算法运行时间复杂度且在如何选择核及配置核参数方面仍有许多理论问题需要研究。此外，上述算法均未充分考虑波段间的相关性，应用于弱相关性波段降维效果不理想。

为解决上述算法对相关性较弱波段降维效果不显著且易忽略能够有效区分不同地物类别特征波段的问题，本文提出一种特征优化降维算法。该算法首先借助影像波段间相关系数矩阵确定初始子空间[8]。然后基于相似性度量准则，利用K-means聚类更新子空间，以保证各子空间内波段具有较强的相关性。最后，利用PCA变换提取第一主成分以实现子空间降维。由于各子空间并集未覆盖全部波段，对剩余波段采用基于矩阵模式的波段选择(band selection method based on matrix mode，BSMM)算法进行处理[9]，利用离散系数和光谱相关系数，定义矩阵因子，并基于该因子计算量化指标，从而选取合适数量的波段组合。将提取出的各第一主成分和波段选择结果叠加作为最终降维结果。较之上述算法，本文算法能够在有效降低影像维度的同时极大限度地保留原始影像信息、增强目标间的可分性，使得降维结果更适用于后续目标识别与分类。

1 算法描述

高光谱影像由于传感器波段范围的重叠、地形因素和地物的反射特性，导致相近波段间具有较强的相关性，波段间信息冗余度较高。使用聚类算法，能够将全部波段划分为几个组内相关性较强、组间差异性明显的波段组集合，每个集合称为一个子空间。鉴于K-means算法处理大数据集的高效性，其在高光谱数据子空间划分中应用较为广泛[10]。但直接采用K-means算法对高光谱影像各波段进行聚类时，存在以下缺陷[11]：①传统聚类算法对初始聚类中心选取较为敏感，若初始聚类中心相距较近，则会将同一聚类划分至不同组别；②聚类个数需要提前给定，无法实现自动判别；③单纯使用一种相似性度量准则作为距离判定标准，鲁棒性较差。为了有效克服上述问题，提出改进的K-means算法对高光谱影像各波段进行聚类分析，用以确定子空间，并对其进行PCA变换，获得降维结果。对于相关性较弱的剩余波段采用BSMM算法进行降维处理，叠加2次降维结果实现高光谱影像的最终降维。

1.1 相似性度量准则

假定X={xi|xi=(xi1，xi2，…，xiN)，i=1，2，…，l}是一幅待处理高光谱影像。其中l为波段总数；N为像素总数；xi=(xi1，xi2， …，xiN)、xi’=(xi’1，xi’2，…，xi’N)分别为第i和i’个波段的N维向量。本文分别选取基于欧式距离、相关系数和信息散度的光谱相似性度量准则作为距离的评价指标[12]，衡量波段与聚类间相似性。

欧氏距离计算如式(1)所示[13]。

(1)

式中：k为像素索引；D(xi，xi ′)表征波段间的欧式距离。若D(xi，xi ′)值越大，表征波段间相关性就越弱；反之，则表征波段间相关性越强。

以相关系数作为统计变量衡量光谱间相似性，其计算公式见式(2)。

(2)

光谱信息散度是典型的基于信息测度的相似性度量准则。设波段xi=(xi1，xi2，…，xiN)，xi’=(xi’1，xi’2，…，xi’N)，分别对应的概率向量为pxi=(pxi1，pxi2，…，pxiN)和pxi ′=(pxi ′1，pxi ′2，…，pxi ′N)，其中，

(3)

光谱信息散度计算见式(4)。

SID(xi，xi ′)=DKL(xi‖xi ′)+DKL(xi ′‖xi)

(4)

其中，

(5)

式中：DKL为相对熵。SID值越高，则波段间相关性越弱；反之，波段间相关性越强。当SID为零时，表征2个波段像素强度完全一致。

1.2 子空间划分及降维

计算影像波段间的相关系数矩阵，设定相关系数阈值T1。从矩阵的第一行开始遍历所有列中相关系数大于T1的波段，若遍历至第li列相关系数均大于T1，划分第1个波段组为1～li波段。而后从li+1行开始，遍历li+1列至最后1列，按上述方法继续搜索相关系数大于T1的波段，构成第2个波段组，以此类推，直至所有波段遍历完毕，结束操作。经大量实验验证，T1存在一个最佳取值范围0.88～0.92，在这个范围内取值基本可以满足实验要求。若T1≥0.92，则会因条件过于苛刻，导致波段组内波段数较少；若T1≤0.88，则会将相关性较弱的波段划分至同一波段组，影响后续实验效果。因此，在最佳取值范围内取均值0.9作为阈值，较为合理。受噪声和异常值等因素的影响，某些波段组内波段数可能较少，选取波段总数的10%作为阈值T2，用以衡量各波段组内的波段数。剔除波段数小于T2的波段组，只保留波段数大于T2的波段组作为初始子空间。经实验表明，T2的设定能够较为准确地判断初始子空间个数。若波段组个数较多，可适当降低T2取值。统计波段数大于T2的组别，计为K组。

为获取相关性较强的子空间，依据上述3种相似性度量准则，应用改进K-means算法分别进行聚类。设置聚类个数为K，并在各初始子空间中随机选取一个波段作为初始聚类中心。假定选取的初始聚类中心为C={cj|cj=(cj1，cj2，…，cjN)，j=1，2，…，K}。其中，j为聚类索引；cj表征第j个聚类中心的N维向量。依次遍历所有波段与初始聚类中心的距离，遵循距离最近原则将波段划分至相应的子空间，若划分结果未发生变化即子空间趋于稳定，停止迭代；否则，重新确定聚类中心，新的聚类中心为各子空间内所有波段的均值，表达见式(6)。

(6)

式中：nj为第j个聚类的波段总数；cj为更新后的聚类中心。取3种度量准则下聚类结果交集，得到最终子空间划分结果，表示为S={sj|j=1，2，…，K}。此时可保证各子空间内波段间相关性较强，对各子空间分别进行PCA变换提取第一主成分，首先计算子空间内波段均值，表达如式(7)所示。

(7)

计算子空间内波段的协方差矩阵，见式(8)。

(8)

tj=ej1(sj-μj)

(9)

式中：tj表示第j个子空间的第一主成分，则各子空间降维结果为T=(t1，t2， …，tK)。

1.3 基于BSMM的波段选择

由于各子空间并集未覆盖全部波段，对于未被子空间覆盖的剩余波段，应用BSMM算法进行波段选择，首先计算各波段离散系数，见式(10)。

(10)

式中：ζi表示第i个波段像元值的标准差；σi表示影像第i波段的离散系数。σi越大，表征第i波段影像所含信息量越大。

其中，

(11)

式中：xik表示第i波段第k个像元的像元值，结合离散系数和相关系数定义矩阵因子mi，i ′，表达见式(12)。

mi，i ′=(σi×σi ′)/SCM(xi，xi ′)

(12)

式中：mi，i ′越大，表征2个波段所含信息量越大，且波段间相关性越小。为避免波段选择过程中，由于谱间高度相关性造成的相关系数未发挥其自身作用的情况，对ζi进行区间映射，映射见式(13)。

σi=(Fmax-Fmin)(σi-σmin)/(σmax-σmin)+Fmin

(13)

式中：Fmax、Fmin分别表征相关系数矩阵的最大、最小值；σmax、σmin分别表征影像各波段离散系数的最大、最小值，且σi取值范围为[0，1]。计算量化指标wi如式(14)所示。

(14)

式中：wi值越大，表征该波段所含信息量越大，且与其他波段间相关性越小。为进行波段选择，将wi值按照从大到小顺次排列，以剩余波段总数的10%为阈值T3，选择相应的最佳波段。将子空间降维结果与波段选择结果叠加作为最终降维结果。

1.4 算法流程

1)计算波段间相关系数矩阵，利用阈值优先约束划分初始子空间。

2)依据初始子空间确定聚类个数及聚类中心，选取上述相似性度量准则，应用K-means算法进行波段聚类，取不同聚类结果交集确定为最终子空间。

3)对各子空间进行PCA变换，提取第一主成分，作为子空间降维结果。

4)对于未被子空间覆盖的剩余波段，采用BSMM算法进行降维。

5)叠加2次降维结果，实现最终降维。

2 实验结果与评价

2.1 实验数据

为了验证算法的有效性，选取2幅影像(图1)进行降维实验，第1幅是由HYDICE光谱仪采集的华盛顿哥伦比亚特区高光谱影像。影像原始大小为1 280×307个像元，有210个波段。实验区截取部分影像为418×188个像元，去除噪声和水吸收波段，剩余191个波段，波长范围为0.4～2.4 μm，光谱分辨率为10 nm，图1(a)是利用27、60、17波段作为RGB分量合成假彩色影像，通过目视解译的方法选取水体、建筑物、树木、草地、裸地和道路6种地物进行分类。第2幅是由ROSIS传感器采集的意大利帕维亚大学地区高光谱影像。影像原始大小为610×340 个像元，有115个波段，去除噪声较多波段，剩余103个波段，波长范围为0.43～0.86 μm，空间分辨率为1.3 m。图1(b)是利用26、66、16波段作为RGB分量合成假彩色影像，包含道路、建筑物、裸地、草地、树木和金属板6种地物类别。实验环境是Intel(R)Core(TM)i5-5200U CPU@2.2 GHz/4 GB内存/Matlab2013b。

图1 实验影像

2.2 降维结果

图2所示分别对应2幅高光谱影像各波段间的相关系数矩阵灰度图[14]，图中亮度越高表征波段间相关性越强。通过目视观察不难发现，原始影像可划分为互不重叠的相关性较强的若干子空间。本文算法通过设定相关系数和最小波段数阈值，计算得到的子空间个数与目视判读一致，进而间接说明本文算法能够较为准确地确定子空间个数。

图2 相关系数矩阵灰度图

基于3种相似性度量准则，分别应用K-means算法进行波段聚类，取不同聚类结果的交集，得到最终子空间。2幅影像对应的子空间划分结果如图3所示。图3(a)对应的4个子空间分别为1～26波段、34～54波段、56～84波段以及103～191波段。

图3 子空间划分结果

图3(b)对应的4个子空间分别为2～24波段、26～41波段、44～69波段以及74～103波段。对各子空间进行PCA变换，分别提取第一主成分，作为子空间降维结果。

显然各子空间并集未覆盖全部波段，对图3(a)子空间划分结果共剩余26个波段，以T3为阈值选取量化指标排在前两位的波段分别为第55和92波段，将其与子空间提取的4个第一主成分叠加，最终将影像降至六维。图3(b)包含103个波段，根据子空间划分结果剩余8个波段，说明影像波段间相关性均较强，以T3为阈值算得波段数不足1，因此不进行波段选择，最终影像降至四维。

PCA和MNF因算法简单且能相对有效地降维，被广泛应用于高光谱影像降维处理；本文采用BSMM对未被子空间覆盖的剩余波段进行波段选择；SuperKPCA是最新提出的超像素核主成分分析算法，对每个超像素进行核主成分分析，用以解决PCA无法处理高光谱影像中普遍存在的非线性特征问题，因此，选取PCA、MNF、BSMM和SuperKPCA作为对比算法，将高光谱影像降维至本文算法相同维度，并评价各算法降维结果，验证本文算法的有效性。

2.3 降维结果评价

为评价各降维算法优劣，本文基于类别可分性、信息量和波段间相关性3种指标进行评价。理想情况下，降维结果信息量大、相关性小、可分性强，则说明降维效果较优。

其中，最有效的是可分性评价方法中的分类精度评价法。对利用PCA、MNF、BSMM、SuperKPCA和本文算法得到的降维结果在相同环境下选取同一样本执行支持向量机(support vector machine，SVM)分类[15]。图4(a)分别对应上述5种算法对华盛顿哥伦比亚特区降维后的分类结果，图4(b)则分别对应于帕维亚大学降维后的分类结果。观察图4(a)，不难看出图4(a1)和图4(a3)中均存在多处建筑物被误分为道路的现象，且建筑物分类结果轮廓较为模糊；图4(a2)中则出现裸地与建筑物误分的情况，同时植被区域分类结果不够完整；图4(a4)相较于前3种分类结果，草地与树木误分现象较严重；相比之下，图4(a5)中各地物类别分类结果较为完整，建筑物轮廓相对清晰，误分像素较少，分类结果较为理想。比较图4(b)，图4(b1)和图4(b2)中大量建筑物被误分为道路，且在类内方差较大的裸地区域存在较多的误分像素；图4(b3)和图4(b4)中右上角金属板被误分成建筑物，同时存在部分建筑物被误分为道路的现象；而图4(b5)中虽有个别建筑物误分成道路，但整体看来，误分像素明显少于图4(b1)～图4(b4)的分类结果。通过目视判读，不难发现本文算法的分类结果较为完整，优于其余对比算法。

对各算法分类结果做定量评价，分别计算其总体分类精度和Kappa系数[16-17]，如表1所示。总体分类精度为准确分类的像素个数与像素总数之比，Kappa系数用以表征算法分类结果与理想分类结果的一致性程度，通常当Kappa≥0.75时，表明分类精度较高。

观察表1不难发现，本文算法的总体分类精度和Kappa系数均高于对比算法。对于华盛顿哥伦比亚特区降维结果，本文算法的总体分类精度为85.9%，Kappa系数为0.802，其总体分类精度和Kappa系数较PCA算法提升7.5%和0.097，较MNF算法提升12.4%和0.108，较BSMM算法提升8.3%和0.101，较SuperKPCA算法提升1.2%和0.019；对于帕维亚大学降维结果，本文算法的总体分类精度为80.3%，Kappa系数为0.752，其总体分类精度和Kappa系数较PCA算法提升5.4%和0.077，较MNF算法提升9.9%和0.101，较BSMM算法提升5.2%和0.059，较SuperKPCA算法提升1.8%和0.02。通过对5种算法分类结果的定量评价，得出本文算法较其余算法具有更强的优越性。

图4 降维后分类结果

表1 总体分类精度和Kappa系数

除用分类结果来评价各降维结果的优劣，本文还采用基于联合熵反映信息量以及相关性这2个评价因子衡量降维结果，如表2所示[18-19]。观察表2发现，对于图1中华盛顿哥伦比亚特区和帕维亚大学2幅影像，本文算法得到的降维结果联合熵分别为3.915 2和2.837 0，高于对比算法，进而定量说明本文算法降维结果所含信息量较大。

表2 联合熵和相关性指标

相关性表征波段间信息冗余度[20]，通常来讲，降维得到的影像各波段间相关性越小，表明信息冗余度越低，降维结果越好，本文使用平均相关性来评价各波段间相关性。观察表中数据可知，本文算法平均相关性相对较低，说明降维结果中各维度间数据相关性不高，进而验证本文算法可行性。

3 结束语

本文利用PCA变换对相关性较强波段降维效果显著的优势，提出基于特征优化的高光谱影像降维算法。算法利用相关系数矩阵自动判别子空间个数，且在初始子空间中选取聚类中心，能够保证初始聚类中心覆盖各个聚类，以提高K-means算法聚类精度。依据3种相似性度量准则分别进行聚类，取聚类结果交集作为最终子空间，保证子空间具有高度相关性，进而采用PCA变换进行降维，使得降维结果更加可靠。对于未被子空间覆盖且信息量大的波段，为防止信息损失，应用BSMM算法进行波段选择。借助不同的降维评价指标，定性定量地分析各降维结果，验证本文算法对高光谱影像降维的有效性及优越性。然而在进行波段选择时，算法设置经验阈值判断波段选择个数，如何更为准确地对剩余波段进行建模，确定其降维维度，将是本文下一步研究方向。