结合局部空谱信息的高光谱图像多端元提取

杨华东，郝永平

(沈阳理工大学 1.信息科学与工程学院；2.机械工程学院，沈阳 110159)

高光谱遥感成像技术同时获取观测目标区域的二维几何空间信息与一维光谱信息，具有光谱分辨率高，成像波段多且连续的特点，使人们能够从遥感图像中获取的信息有了质的转变[1]。然而，由于高光谱传感器空间分辨率偏低及地物复杂多样等因素，图像中多数像元的光谱信号都是多种地物光谱信号的综合，即混合像元现象。混合像元现象给地物准确识别和精细分类带来极大困难，严重阻碍了高光谱遥感技术定量化应用。

高光谱解混技术将高光谱图像中的混合像元光谱分解为一组组分光谱，称为端元，与相应组分所占的比例，称为丰度，是解决混合像元的有效方法[2]。提取混合像元中基本组成成分特征光谱的过程称为端元提取，而计算各个端元在混合像元中所占比例的过程称为丰度反演。端元提取是混合像元分解的必要前提，端元提取的好坏是混合像元分解效果的关键。

当前研究中，基于光谱数据在其特征空间中呈凸面单体的几何特性，从特征空间搜索单体顶点是端元提取最成功、最广泛的方法。典型方法包括：像元纯度指数(PPI)[3]、自动目标生成(ATGP)[4]、顶点成分分析(VCA)[5]、NFINDR[6]、单体增长算法(SGA)[7]、正交基方法(OBA)[8]、最小体积约束非负矩阵分解(MVC-NMF)[9]、最小体积单体分析(MVSA)[10]等。还有一些学者针对图像中的空间信息进行研究，发展了融合空间信息的端元提取方法，如自动形态学端元提取(AMEE)[11]、空间预处理(SPP)[12]端元提取算法等。这些方法分别在不同技术发展阶段获得了较好的端元提取效果，但这些方法都有一个共同的缺陷：端元提取未考虑端元变异性，只为每类地物提取出一个端元，端元光谱代表性差，很难准确刻画一个地物类别，使得混合像元分解精度难以提高。

针对端元变异问题，本文提出一种结合局部空间光谱信息的多端元提取方法(MEELSI)，通过将原始高光谱图像划分为若干不重叠的图像块，然后进行分块端元提取，并利用局部空间信息精选优化端元光谱，最后对端元集聚类分组实现地物光谱多样化表达。与传统仅基于光谱信息的单端元提取方法相比，MEELSI在端元提取过程中充分利用了图像的光谱信息和局部空间约束，提取的端元能更好地刻画地物类别，有效提高混合像元分解精度。另一方面，MEELSI可以为影像中的每个地物类别提取多个端元光谱，有效解决遥感影像中的光谱变异问题。

1 线性光谱混合模型

(1)

基于线性混合模型，L维像元向量r可以看作是高光谱数据特征空间中的一个点，在满足非负和全加性约束条件的情况下，全部光谱数据在该空间扩展成一个凸面单形体(也称单体)，而端元恰好位于这个单形体的顶点处，因此，端元提取问题就转化为寻求凸面单形体的问题。

2 算法描述

本文提出的结合局部空谱信息的MEELSI算法包含四个关键步骤：(1)图像子空间划分；(2)基于ATGP的端元自动提取；(3)空间信息约束下的端元优化选择；(4)基于光谱相似性的端元聚类。

2.1 图像子空间划分

遥感图像中，受地势、光照、环境等因素影响，通常会发生光谱变异现象，即同一地物类型在不同的区域所呈现的特征光谱曲线不尽相同，即“同物异谱”现象，图像空间跨度较大时更是如此。因此，在不同的图像空间区域分别提取端元光谱，为每类地物获取多个端元光谱更能刻画地物类型属性，同时也可以有效解决端元光谱变异问题，提高端元提取精度。

MEELSI算法第一步，将原始高光谱图像划分为不重叠的图像子区域(图像块)。本文采用最简单的方法对图像进行分割，即以固定的行列数目对图像进行划分，并限制分割后的图像块大小为20×20。

2.2 基于ATGP的端元自动提取

在利用ATGP进行高光谱图像端元提取之前，需要确定端元个数。Chang等提出了虚拟维数(virtual dimensionality,VD)[13]概念来解决这个问题，并给出了相应的确定VD噪声白化的HFC(noise-hitened Harsanyi-Farrand-Chang,NWHFC)算法。本文中高光谱图像的端元个数q采用NWHFC算法确定，然后利用ATGP算法分别在每个图像子区域提取候选端元。

ATGP自动提取端元过程简述如下。设由NWHFC方法确定的端元个数为q，像元向量为r。选择图像中具有最大长度的像元向量作为初始端元向量t0，选择公式如式(2)所示。

(2)

(3)

式中UT和U-1分别是U的转置矩阵和逆矩阵。

2.3 空间信息约束下的端元优化选择

ATGP算法易受噪声和异常像元的影响。文献[11-12]中指出，图像中纯像元一般位于空间同质区，而混合像元通常位于地物类型过渡区或异质区。因此，在候选端元邻域内(以候选像元为中心大小为d×d的窗口内，d=3,5,7)搜寻与候选端元光谱近似的像元光谱，以其均值光谱作为新的端元光谱，可以在一定程度上抵消噪声影响。另一方面，异常像元与其邻域像元的光谱差异性较大，所以通过比较候选端元与其邻域像元的光谱相似性，可以有效剔除异常像元。

本文采用光谱角距离(Spectral Angle Distance，SAD)衡量光谱相似性，像元x和y的SAD定义如下。

(4)

SAD的值越小表明两个像元的光谱相似性越高。通过合理设定SAD阈值区间，可以甄别异常像元和同质区像元。

异常像元光谱与其邻域像元光谱差异显著，光谱角的值较大，因此，判定异常像元的SAD阈值一般较大，本文实验中设为10(度)。而同质区内包含相同地物类型的纯像元光谱差异较小，这类光谱差异主要由噪声引起。因此，SAD阈值一般设为2～5(度)。实验中SAD阈值的大小与图像空间分辨率、噪声水平以及地物分布等因素都有关系，应根据实际情况作相应调整，如噪声水平低或空间分辨率高时，阈值应适当减小；反之亦然。

2.4 基于光谱相似性的端元聚类

在对所有图像子区域进行空间信息约束下的端元优化后，利用聚类技术对所有端元进行聚类分组，本文采用K-means聚类算法[14]对端元集分类。算法中聚类个数q由NWHFC方法确定，q个初始聚类中心通过在端元集上执行ATGP算法获得。在迭代更新聚类中心时，采用SAD作为距离度量。

2.5 MEELSI算法流程描述

MEELSI方法流程如图1所示。首先对高光谱图像进行图像空间区域划分，然后在所有图像块上利用ATGP算法提取候选端元，得到整幅图像的候选端元集；为提高端元提取精度，利用光谱相似性测度和局部空间信息分析对所有候选端元进行优化，剔除候选端元集中的异常像元和纯度较低的候选端元；最后利用K-means聚类算法对优化后的端元集进行分组，获得地物的多端元光谱集。

图1 MEELSI算法流程

3 实验结果与分析

为验证结合局部空谱信息的高光谱图像多端元提取算法的正确性和有效性，这一部分将选择ATGP和NFINDR算法与本文提出的MEELSI算法在仿真高光谱影像和真实高光谱影像数据上进行实验验证。实验采用的计算机硬件环境为Inter(R) Core(TM) i7-4700MQ CPU 2.4GHz，内存8GB，软件环境为Microsoft Windows 8.1专业版，Matlab R2016b。

为评价端元提取和丰度估计精度，分别采用SAD和均方根误差(Root Mean Square Error,RMSE)两个指标定量比较。SAD计算公式见式(4)，RMSE计算公式如式(5)所示。

(5)

3.1 仿真高光谱数据实验

从美国地质调查局(USGS)光谱库中选取3类地物6种端元光谱，分别为明矾石(Alunite-1和Alunite-2，简记为A1和A2)、水铵长石(Budding-tonite-1和Buddingtonite-2，简记为B1和B2)和高岭石(Kaolinite-1和Kaolinite-2，简记为K1和K2)。利用上述6种光谱合成一幅具有200个波段120×90的高光谱图像，示意图见图2，生成方式如下。

图2 合成图像示意图

将图像分为12个大小为30×30的图像块，每个图像块内的像元由两种端元光谱混合而成，混合比例为1:1，从左至右，自上而下端元光谱组成情况依次为(A1,B1)、(A1,B2)、(A2,B1)、(A2,B2)、(A1,K1)、(A1,K2)、(A2,K1)、(A2,K2)、(B1,K1)、(B1,K2)、(B2,K1)、(B2,K2)。为了模拟纯像元，分别用大小为10×10的纯像元块替换上述12个图像块的中心区域(图中阴影区域)，纯像元块的组成依次为K1、K2、K1、K2、B1、B2、B1、B2、A1、A2、A1、A2。

为模拟真实的高光谱影像，对影像先进行模板大小为9×9的均值滤波处理，然后附加信噪比(SNR)为30dB的高斯白噪声。

表1是MEELSI算法与ATGP、NFINDR算法获得的端元光谱与真实端元光谱的SAD(以弧度值表示)对比结果。

表1 提取的端元光谱与真实端元光谱的SAD

由于ATGP和NFINDR都是单端元提取方法，因此只能为每类地物提取一条光谱，表1中“*”号表示该方法未提取到这种端元光谱。总体来看，MEELSI算法可以为每类地物提取多个端元代表，而且MEELSI算法的SAD值最小，精度最高。

为进一步验证端元提取结果对光谱解混精度的影响，本文分别利用ATGP、NFINDR和MEELSI算法提取的端元进行丰度估计。对于单端元提取算法ATGP和NFINDR，采用经典的全约束最小二乘估计法(Fully Constrained Least Squares,FCLS)[15]估计丰度。对于本文提出的多端元提取算法MEELSI，则利用经典的多端元光谱混合分析(Multiple Endmember Spectral Mixture Analysis,MESMA)[16]估计丰度。表2是丰度估计的RMSE对比结果。

表2 真实丰度与估计丰度的RMSE

从表2中数据可以看出，对仿真高光谱图像进行多端元提取和多端元混合像元分解，其光谱解混精度要高于传统单端元光谱解混。

3.2 真实高光谱数据实验

真实高光谱影像是由AVIRIS成像仪于1997年获取的位于美国内华达州南部Cuprite矿区的高光谱数据。图像包含224个波段，空间分辨率为20m，光谱分辨率为10nm，光谱范围为0.4～2.5μm。实验中选取该影像的200×200子区域作为真实实验数据，如图3所示。删除水吸收和低信噪比波段后，剩余188个波段的数据供实验使用。

图3 真实高光谱影像

利用NWHFC方法估计的端元数目为14，参考Clark等[17]在该地区的实地调查和矿物填图结果，本文在实验时，选取该地区内分布较多且发生光谱变异较多的5种典型地物类型进行研究。ATGP、NFINDR和MEELSI算法提取的端元光谱与USGS光谱库中相应矿物光谱的SAD结果见表3所示。

表3 Cuprite数据实验结果

从表3中可以看出，MEELSI提取的端元光谱和光谱库中的矿物光谱总体上最为相似。

4 结论

提出了结合局部空间和光谱信息的高光谱图像多端元提取方法，提取的端元能更好地刻画地物类别，可以有效解决遥感影像中的光谱变异问题，提高了端元提取和光谱解混精度。在仿真数据集上，发生光谱变异的端元提取成功率达到100%，真实高光谱图像数据集上的端元提取结果也更接近光谱库中的参考光谱。未来将进一步改进本文的方法，使其适用于数据混合度较高的高光谱图像，同时设计更有效的空间信息分析策略，进一步提高算法精度。