一种改进的顶点成分分析端元提取算法

袁博，张杰林

（1.中国矿业大学（北京），地球科学与测绘工程学院，北京100083；2.核工业北京地质研究院，遥感信息与图像分析技术国家级重点实验室，北京100029）

一种改进的顶点成分分析端元提取算法

袁博1，张杰林2

（1.中国矿业大学（北京），地球科学与测绘工程学院，北京100083；2.核工业北京地质研究院，遥感信息与图像分析技术国家级重点实验室，北京100029）

顶点成分分析算法需要预先提供端元数目，端元数目正确与否对结果会产生较大影响，其算法在实际应用中多次运行的结果不稳定。针对上述缺点提出了一种改进的顶点成分分析端元提取算法。该方法在n维光谱空间中生成n个彼此正交的单位向量，在此基础上生成与之具有一定夹角的单位向量，将光谱空间中的像元点分别投影在单位向量上以获取端元。结果表明改进的顶点成分分析端元算法提高了端元提取结果的稳定性。

端元提取；顶点成分分析；纯净像元指数；正交向量；蚀变信息

高光谱遥感技术是20世纪80年代发展起来的首次实现图谱合一的对地观测技术，其光谱分辨率可以达到5～10 nm［1］，可为每个像元提供数十个或数百个波段的光谱信息，可以产生一条完整并且连续的光谱曲线，使得许多在多光谱遥感中不能分辨的物质得以识别。地球表面并不是由单一的物质组成的，当遥感图像中的一个像元包含了不同波谱特性的物质时就会产生混合像元，在高光谱遥感中混合像元是广泛存在的，高光谱遥感图像中每个像元的光谱信息是其对应地表物质光谱信息的综合，不同的物质具有不同的光谱响应特征，而一个像元仅用一条波谱曲线来代表这些特征［2］，因此对高光谱遥感数据进行分析就必须进行混合像元的解混。当像元中只有一种地物时，该像元就被称为端元，端元光谱就是该区域的特征光谱，因此端元提取是解混的前提。

像元混合分为线性混合和非线性混合。如果像元内物质混合的尺度大，那么混合像元被视为线性混合；如果像元内物质混合尺度小，混合紧密，则混合像元被视为非线性混合［3］。因为不同时间地点获得的高光谱数据具体的非线性混合模式复杂多变，通常很难找到与实际情况相吻合的非线性混合模型，而线性混合模型具有简单易于操作并且精度满足需要的优点，因此目前线性混合模型应用比较广泛。顶点成分分析（Vertex Component Analysis，VCA）是在线性混合模型基础上提出的端元提取算法，该算法需要提前确定端元数目，并且多次提取端元的结果不稳定，针对VCA的这一缺点笔者提出了一种改进的顶点成分分析端元提取算法。

1 端元提取算法研究现状

目前比较常用的端元提取算法包括纯净像元指数法（Pure Pixel Index，PPI）、内部体积最大法（N-FINDR）、单形体体积增长算法（Simplex Growing Algorithm，SGA）、VCA算法等［4］。

PPI算法是Boardman（1995）提出的在多光谱和高光谱影像中寻找波谱最纯净的像元，从而提取端元的方法［2］。该算法可以从高光谱数据中分离出光谱较纯净的像元，减少参与端元选取的像元数量。之后利用纯净像元生成n维散点图（n指纯净像元个数），进行交互式可视化操作选出图像端元［5］。PPI方法使用随机初始化向量进行投影，具有随意性，不可避免地产生冗余向量，不能有效区分各个端元的投影，并且该方法用到大量的投影变换，计算效率较低，消耗时间较长。

Winter（1999）提出了内部体积最大法（NFINDR），通过寻找具有最大体积的单形体自动地逐步获取影像中的端元［6］。该方法首先随机选择一定数量的像元作为端元，并计算以这些端元为顶点构成的单形体的体积，再将新的像元依次代替端元，重新计算单形体的体积，如果体积增大，则该像元替代之前端元，直到没有像元能使单形体的体积增大时可以得到的所有的端元［6］。但是当数据存在噪声时，该算法使得一些点将不可避免地出现在单形体之外，因此数据含有噪声时结果准确性将降低。

SGA算法是对N-FINDR算法的改进，通过顺序搜索方法来寻找构成最大单形体体积的端元［7］。该方法需要重复的体积计算，计算复杂度比较高，因此效率较低。

VCA算法理论基础为线性混合模型，一个像元矢量的光谱信号用线性混合模型描述为［8］:

式（1）:B—L×1维的观测波谱，其中L—高光谱图像的波段数；M＝［m1，m2，…，me］—一个L×e的端元集矩阵，其中mi（i＝1，2…e）—端元的光谱，其中e—端元个数；k—形态学比例因子；A＝［a1，a2，…，ae］T，—e×1维的丰度值向量；n—存在于观测波谱B各谱段中的L×1维高斯白噪声。

在线性混合模型下，观测向量集合组成一个凸锥C＝{B∈RL:R＝kMA，1TA＝1，A≧0，k≧0}。凸锥C投影到超平面RTu＝1形成一个单形体S＝{y∈RL:y＝R/（RTu），R∈CR}。VCA首先随机生成一个初始向量，找到单形体S中在初始向量上具有最大投影长度的像元作为第一个端元，然后计算与该端元正交的向量作为第二个投影向量，将数据投影找到第二个端元。然后依次生成新的与已提取端元所张成的空间相正交投影向量，提取相应的端元，当端元数目满足预先确定的数目e时算法结束［9］。以图1中两个端元为例，第一次迭代中，单形体数据被投影到第一个方向f1，最大的投影向量对应着端元m1。接下来的迭代中，将数据投影到与端元m1正交的方向f2，此时最大的投影向量便对应着端元m2。

VCA算法需要事先确定要提取的端元数目，端元数目正确与否对丰度值求解具有较大影响，但在实际很多情况下不能预先确定研究区域内端元个数，此时VCA方法不具有适用性［8］。另外VCA方法在没有误差和噪声的理论前提下经过n次投影确定出n个端元，但是高光谱数据中不可避免地含有噪声，实际投影所得结果可能是噪声，因此只凭一次投影就确定一个端元会存在很大的误差。通过VCA方法多次确定出的端元集差异性较大，其多次运行结果不稳定，该方法的准确率较低。

何晓宁等在2013年对VCA算法进行了改进。该改进算法利用图像空间信息排除噪声，确定候选端元，再通过病态矩阵分析候选端元矩阵的有效性，经过迭代后获得端元［8］。该方法需要在计算前提供端元数目，但是在实际应用中并不总能确定端元数目，并且该算法具有计算量大的缺点。

图1 VCA原理示意图Fig.1Schematic diagram of VCA

2 对顶点成分分析（VCA）的改进算法

高光谱数据可以被视为在n维空间中的点，点的坐标是各个波段的反射率值。改进的VCA端元提取算法首先对数据进行最小噪声分离（Minimum Noise Fraction，MNF）变换，去除数据中的噪声，确定含有有用信息的最小波段数n，然后在变换后的n维空间中生成一个不与原始数据集正交的单位投影向量a1，再计算出与向量a1正交的单位向量a2作为第二投影向量，此时单位向量a1和单位向量a2组成一个二维向量空间。生成第三个单位向量a3，使之与向量a1和a2所张成的空间正交。以此类推直到第n个单位向量an，使an与a1、a2、…、an-1构成的n-1维向量空间正交，这里n指波段数，此时得到n个彼此之间相互正交的单位向量，即找到了n维空间中的一组基向量，n维空间中的任何像元点都可以有效地投影到单位向量上。

为了保证有效的提取端元，需要适当增加投影向量的数量。计算与已知n个向量中每两个向量ai和aj均成π/4夹角的中间单位向量aij，i＝1…n，j＝i+1…n，经计算可得到n（n-1）/2个中间单位向量aij，于是最终得到n＋n（n－1）/2个投影向量，n为选取的波段数。最后将n维空间中的所有像元点投影到单位向量上，记录到每个点投影后为最大投影的次数，然后设定阈值，选择最大投影次数大于阈值的像元为端元。

图2 改进的VCA原理示意图Fig.2Schematic diagram of improved VCA

改进的VCA端元提取算法原理如图2所示，ai为生成的单位向量，aj为与ai正交的单位向量，aij为与ai、aj向量夹角均为π/4的单位向量，图中的点1、点2、点3代表3个不同像元在n维波谱空间的位置，虚线表示投影方向。记录像元点在各个向量上投影后具有最大投影的次数，点1被记录为最大投影的次数为1，点2被记录为最大投影的次数为2，点3被记录为最大投影的次数为3。

在该算法计算之前，高光谱数据首先经过MNF处理，选取特征值较大数十个波段，总的投影向量个数控制在数百个。该方法原理简单，确保了n维空间中不同位置顶点有同等机率投影到端点位置，最终选取数百个投影向量，相比VCA方法，适当增加了投影向量，避免了投影向量产生的随意性，具有一定的抗噪声能力，使提取的端元信息更准确，保证了端元提取结果的稳定性。

3 实验分析

3.1 数据来源

在IDL8.0环境下实现改进的VCA端元提取算法，并采用核工业北京地质研究院采集的SASI航空高光谱遥感测量数据验证算法的实用性。实验区域的面积约为200 km2。

SASI是加拿大生产的机载短波红外（SWIR）高光谱成像系统，其主要技术指标见表1。

表1 SASI航空成像光谱仪主要技术指标Table 1Main technical performances of SASI airborne imaging spectroradiometer

因为SASI航空高光谱数据中1 370～1 415 nm及1 835～1 925 nm波谱范围内共有11个波段受大气水汽影响严重，故将其去除。最终选择950～1 355 nm、1 430～1 820 nm、1 940～2 450 nm共90个波段参与计算。

3.2 数据处理与结果

首先对数据进行预处理，采用统计学模型中的经验线性法（Empirical line），利用已测得的白板数据作为已知点的地面光谱对高光谱影像进行反射率的计算。之后对数据进行MNF操作去除噪声，得到特征值分布（图3）。

图3 MNF特征值分布曲线Fig.3Curve of the eigenvalues of MNF

根据图3中的特征值曲线并观察MNF图像可知处于20右侧的波段趋于一致为噪声波段，因此选择前20个波段参与之后的运算。将选取的波段带入改进的VCA提取端元算法，经过210次迭代运算，程序自动选择出端元，生成基于空间位置的端元分布图，为了验证本算法的稳定性，经过多次计算生成的端元空间位置分布见图4。

图4显示改进的VCA算法多次提取结果的端元分布图中端元位置保持一致，验证了改进算法的结果具有稳定性。而VCA算法多次提取端元的结果不一致，并且不能将图像中的端元全部提取出来。

通过改进的VCA端元提取算法可以获得端元波谱曲线，再利用ENVI软件中自带的波谱分析功能将端元波谱与USGS波谱库中数据进行对比分析可以识别出各个端元。因为选取的SASI数据为90个波段，光谱分辨率为15 nm，所以该步骤需要事先将USGS波谱库数据重采样到同样波段数和光谱分辨率才可以进行。对比USGS波谱库中光谱曲线，经分析可以得到图像中的端元包括白云母、片沸石、黄钾铁矾和深绿辉石4种蚀变矿物。

白云母（K2Al4（Si6Al2O20）（OH，F）4）具有440 nm、900 nm Fe2＋吸收谱带；1 400 nm-OH倍频、合频谱带；在2 200 nm、2 350 nm、2 450 nm存在吸收特征［10］，其诊断性光谱特征是Al-OH键在2 180～2 228 nm之间的尖锐而深的单一吸收特征，以及在2 340 nm和2 440 nm附近的较弱Al-OH特征［11］；片沸石（（Ca，Na）2～3Al3（Al，Si）2Si13O36·12H2O），具有965 nm、1 160 nm、1 430 nm、1 940 nm水引起的吸收特征；黄钾铁矾（KFe3（OH）6（SO4）2），具有1 400 nm、1 900 nm、2 500 nm-OH振动谱带，900 nm Fe3＋吸收波谱特征［10］；深绿辉石（Ca（Mg，Fe3＋，Al）（Si，Al）2O6），具有1 050 nm Fe2＋谱带和850 nm、450 nm Fe3＋谱带［10］。

由改进算法提取的端元波谱与VCA提取的端元波谱及波谱库波谱对比（图5）可见改进算法比VCA算法提取的波谱更接近于波谱库。

以所提取的4种端元的光谱曲线为基础利用混合调制匹配滤波方法进行蚀变信息的填图，经运算可得蚀变信息提取结果（图6）。结果表明以改进算法提取的端元光谱为依据能更准确更详细地提取出蚀变信息。

图5 光谱曲线对比Fig.5Comparison of spectral curve

图6 蚀变信息提取结果Fig.6Extracting results of alteration information

（，Continued on page 44）（，Continued from page 38）

4 结论

提出了改进的VCA端元提取算法，计算n维空间的一组基向量并选择与之成特定夹角的单位向量共同作为投影向量，在选出的数百个单位向量上进行像元投影，确定阈值选出端元。改进的VCA算法无需预先确定图像的端元数目，通过适量投影自动选择端元，算法同时改进了VCA算法不具有抗噪声能力的缺点，在保证计算效率的同时提高了端元提取结果的稳定性。

［1］芶盛.高光谱遥感图像光谱特征提取与匹配技术研究［D］.成都：成都理工大学，2011.

［2］林娜，杨武年，刘汉湖.基于高光谱遥感的岩矿端元识别及信息提取研究［J］.遥感信息，2011（5）：114-117.

［3］Hassan N，Hashim M.Decomposition of mixed pixels of ASTER satellite data for mapping Chengal（Neobalanocarpusheimiisp）tree［C］//Control System，Computing and Engineering（ICCSCE），2011 IEEE International Conference on.Penang：IEEE，2011：74-79.

［4］夏威.高光谱遥感图像解混和波段选择方法研究［D］.上海:复旦大学，2013.

［5］Boardman J W，Kruse F.Analysis of Imaging Spectrometer Data Using N-Dimensional Geometry and a Mixture-Tuned Matched Filtering Approach［J］.GeoscienceandRemoteSensing，IEEE Transactions on，2011，49（11）：4 138-4 152.

［6］钱进，邓喀中，刘冬.基于优化N-FINDR算法的高光谱遥感影像矿物识别［J］.金属矿山，2012，08：84-87，91.

［7］王丽姣.基于单形体体积增长的高光谱图像端元提取及快速实现［D］.杭州:浙江大学，2015.

［8］何晓宁，曹建农，高坡，等.基于顶点成分分析法的端元提取改进算法［J］.测绘通报，2013（7）：30-34.

［9］Nascimento J M P，Dias J M B.Vertex component analysis：A fast algorithm to unmix hyperspectral data［J］.IEEE Transactions on Geoscience and Re mote Sensing，2005，43（4）：898-910.

［10］燕守勋，张兵，赵永超，等.矿物与岩石的可见-近红外光谱特性综述［J］.遥感技术与应用，2003，18（4）：191-201.

［11］刘圣伟，甘甫平，闫柏琨，等.成像光谱技术在典型蚀变矿物识别和填图中的应用［J］.中国地质，2006，33（1）：178-186.

An improved algorithm of end member extraction based on Vertex Component Analysis

YUAN Bo1，ZHANG Jielin2
（1.China University of Mining&Technology（Beijing），College of Geoscience and Surveying Engineering，Beijing 100083，China；2.National Key Laboratory of Remote Sensing Information and Image Analysis Technology，Beijing Research Institute of Uranium Geology，Beijing 100029，China）

Vertex Component Analysis algorithm needs to determine the number of end members at the beginning，therefore the determined number has great influence on the extract results，and the algorithm was found not stable in practical applications.To overcome the shortcomings of Vertex Component Analysis，this paper proposed an improved algorithm.This improved algorithm first find orthogonal unit vectors in the n-dimensional spectral space，and generates a number of unit vectors which have certain angle with the original unit vectors，then pixel points in the spectral space are projected on all of the unit vectors to extract end members.The result shows that this improved algorithm improves stability of the extraction result of end members.

extraction of endmembers；Vertex Component Analysis；Pure Pixel Index；orthogonal vectors；alteration information

TP79；P57

A

1672-0636-（2016）01-0033-06

10．3969/j．issn．1672-0636．2016．01．006

2015-08-20

袁博（1989—），男，河北保定人，中国矿业大学（北京）硕士研究生；研究方向:高光谱遥感。E-mail:boziaaaaa@163.com