基于代数余子式的N-FINDR快速端元提取算法

李 琳孟令博孙 康赵永超

①（中国科学院电子学研究所 北京 100190）

②（中国科学院大学 北京 100049）

基于代数余子式的N-FINDR快速端元提取算法

李 琳*①②孟令博①②孙 康①②赵永超①

①（中国科学院电子学研究所 北京 100190）

②（中国科学院大学 北京 100049）

基于高光谱图像特征空间几何分布的端元提取方法通常可分为投影类算法和单形体体积最大类算法，通常前者精度不好，后者计算复杂度较高。该文提出一种基于代数余子式的快速N-FINDR端元提取算法（FCA），该算法融合了投影类算法速度快和单形体体积最大类算法精度高的优势，利用像元投影到端元矩阵元素的代数余子式构成的向量上的方法，寻找最大体积的单形体。此外，该算法在端元搜索方面较为灵活，每次迭代都可用纯度更高的像元代替已有端元，因此能保证用该端元确定的单形体，可以将特征空间中全部像元包含在内。仿真和实际高光谱数据实验结果表明，该文算法在精准提取出端元的同时，收敛速度非常快。

图像处理；高光谱；端元提取；单形体；体积最大；代数余子式；投影

1 引言

高光谱遥感图像空间分辨率通常相对较低，再加上自然界地物的复杂性、多样性，故混合像元普遍存在于高光谱遥感图像中。而混合像元的存在正是像元级遥感分类和面积测量精度难以达到使用要求的主要原因。因此，光谱解混是分析高光谱图像的关键所在［1］。在光谱解混过程中，端元提取算法可以提取出图像中只包含一种特征地物的纯像元。

常见的端元提取算法通常是在线性混合模型成立条件下［2，3］，基于高光谱图像特征空间几何分布的。该几何分布是指图像中的像元分布在以端元为顶点的单形体结构中［4］。目前，绝大多数端元提取算法都致力于寻找单形体的顶点，具体可分为基于特征空间的投影算法和单形体体积最大算法。其中投影算法有：纯像元指数算法（Pixel Purity Index， PPI）［5］，顶点成分分析算法（Vertex Component Analysis，VCA）［6］，正交子空间投影算法（Orthogonal Subspace Projection， OSP）［7］等。这类算法通常将像元投影到一组向量上，位于线段端点的像元是纯像元的可能性最大。其优点是复杂度较低，缺点是投影向量的选择具有一定随机性，端元提取的结果很大程度上依赖于向量的选择，因此精度不高。而另一类基于单形体体积最大的算法则精度较好。如N-FINDR算法［8］，通过寻找具有最大体积的单形体来自动获取图像中的所有端元，但是其算法复杂度很高，因而应用受到一定限制。故有很多算法对其改进，如单形体增长算法（Simplex Grow ing A lgorithm， SGA）［9］和基于Householder变换的体积最大算法（Maximum Volume based on Householder Transformation， MVHT）［10］等。尽管SGA和MVHT算法在计算复杂度上低于N-FINDR，但是其寻找端元所用的是贪婪算法，即从第1个初始端元出发，逐步一个个地确定后面的端元，因而这种搜索方法总是做出在当前看来最好的选择，从局部而非整体考虑结果。此方法缺点在于一旦端元入选后就不能将其删除，故容易陷入局部极值。而体积最大算法（New Volum e Form u la for Sim p lex， NVFS）［11］虽然在端元搜索方面和N-FINDR一样，即：已寻找到的端元仍可被其他像元替代，但是其计算复杂度比N-FINDR还高，因此应用也受到了一定限制。

由此可知，投影算法虽然复杂度较低，但精度不高；单形体体积最大算法虽然精度较好，但算法复杂度较高。所以在此背景下，本文提出一种基于代数余子式的N-FINDR快速端元提取算法（A Fast endm em ber ex traction m ethod based on Cofactor of a determ inant A lgorithm， FCA）。该算法是对N-FINDR的改进，有如下几点创新和优势：首先，FCA本质思想是寻找最大体积的单形体，故克服了投影算法精度低的缺点；其次，FCA可将图像全部像元投影到一个端元矩阵元素产生的代数余子式的向量上，故具有投影算法复杂度低的优势；最后，在端元搜索方面，FCA采用随机搜索算法，所选择的端元在计算过程中可随时被其他纯度更高的像元替代，故克服了SGA和MVHT采用的贪婪搜索算法的缺点，减少了提取出异常极值的情况。综上，FCA兼备投影算法和单形体体积最大算法的优点，可以高精度高速度地提取出纯像元。

2 线性混合模型和N-FINDR算法

2.1 线性混合模型

线性光谱混合模型（Liner M ixture M odel，LMM）假设：太阳入射辐射只与一种地物表面发生作用，物体间没有相互作用。故图像中任意一个像元光谱r可以表示为

其中ie为端元光谱，p为图像中端元个数，ic为像元光谱r中每个端元所占的丰度，n为误差和噪声。当图像满足线性混合模型，且不考虑噪声影响时，图像中的所有像元都被包含在一个以端元为顶点的单形体中［12］。如图1所示，2个波段3个端元的图像像元分布在特征空间中以端元为顶点的三角形中。故可通过寻找构成特征空间中最大单形体体积的顶点来提取端元。N-FINDR即是基于这种思想的端元提取算法。

2.2 N-FINDR算法

N-FINDR通过计算式（3）和式（4）来提取一组端元。

3 FCA算法描述

FCA算法是建立在线性混合模型成立条件下，高光谱图像特征空间中单形体体积最大的一类方法。本质上和N-FINDR原理一样，但N-FINDR需要大量运算，尤其是很多矩阵运算需要大量迭代，时间复杂度高，不利于工程实际中对高光谱图像的大量处理。本文通过将图像所有像元投影到以端元矩阵元素的代数余子式构成的向量上的方法对N-FINDR第2种实现方式进行改进，使算法在保持端元提取精度的同时，时间复杂度大大降低。

图1 2维特征空间下含有3个端元的图像像元三角形结构（图中圆圈“o”代表像元）

3.1 算法原理及步骤

应用元素与其代数余子式乘积之和求行列式det（E）的方法可以降低算法复杂度。即

由式（6）可知，元素ei，1的代数余子式的值与ei，1的具体值无关，只与ei，1在方阵中的位置有关。因此，用不同像元代替端元e1，相同位置的元素对应的代数余子式相同。故可以得到端元e1的p个元素所在位置对应的代数余子式向量为

由式（6）可知，M e1对于端元e1，即第1列的所有元素而言是常量，故可将图像的N个像元投影到这一向量上。

图2 图像像元在1M e上的投影（图中圆圈“o”代表像元）

此外，FCA的终止条件和N-FINDR相同，即：当两次迭代得到的最大体积相同时，算法结束。这种终止条件保证了在前一次迭代中确定的端元，在下一次迭代过程中都有可能被其他纯度更高的像元替代，避免了如SGA和MVHT算法在确定一个端元后不可替换的缺点，保证了一直从整体角度考虑问题。

3.2 FCA算法步骤

步骤1 根据虚拟维度（Virtual Dim ensionality，VD）［13］计算出图像中的端元个数p。

步骤2 初始化方阵E

（1）用PCA等降维方法将原图像光谱维降到p-1维；

（2）任取降维后图像中p个像元的光谱作为e1， e2，…， ep，代入式（3）中矩阵E的第2行至第p行，同时矩阵E的第1行置1，此时矩阵E为p阶方阵；

（3）设Vm ax= 0；迭代次数it= 0。

步骤3 提取端元

（1）方阵E去掉第j列，即第j个端元光谱向量，得到矩阵M a， j=1，2，…， p；

（2）去掉矩阵M a的第i行，即可得到原方阵E第i行第j列的元素所对应的代数余子式的方阵，设该方阵为M b，阶数为p-1， i=1，2，…， p；

（3）计算第i行第j列元素所对应的代数余子式：（-1）（i+j）×det（M b）；

（4）分别计算出第j列每个元素所对应的代数余子式的值，得到一个含有p个代数余子式值的行向量；

（5）用该代数余子式的行向量乘以像元矩阵，即可得到N个p阶方阵行列式的值。N为图像像元个数；

（6）用N个行列式中绝对值最大的值对应的像元向量代替矩阵E的第j列，直到j= p时结束；

（7）根据式（4）计算出此时的体积V（ E）；

（8）如果V（ E）＞Vm ax ，则令Vmax=V（ E），it=it+ 1；否则退出循环，算法结束。

4 实验结果与讨论

4.1 仿真数据实验

仿真数据由已知特定的5种端元光谱和相应的丰度矩阵构成。端元光谱均从USGS光谱库中得到，分别为：A luniteGDS83， CalciteWS272， Desert_ VarnishGDS141， KaoliniteCM 9， Nontronite GDS41，其波段数为210。丰度矩阵由Dirichlet［5］分布构成，满足丰度非负与和为1约束［14］。将端元矩阵与丰度矩阵相乘，并加上不同强度的高斯白噪声，即得到了本实验所用的仿真数据。

试验中参与比较的5种算法为：N-FINDR，SGA， MVHT， NVFS， FCA，均为基于单形体体积最大的端元提取算法。在端元提取的基础上，利用全约束最小二乘法（Full Constraint Least Square，FCLS）［15］进行丰度估计，最后，通过分析算法执行时间和端元提取以及丰度估计的精度综合评价5种算法的性能。

实验1 算法复杂度分析

如图3（a）和图3（b），当像元数N和端元数p分别增加时，NVFS复杂度最高，FCA比N-FINDR复杂度低了很多，特别是当p逐渐增加时，N-FINDR和FCA的算法复杂度差距会进一步加大。此外，在p和N增加的两种情况下，FCA的算法复杂度几乎和MVHT相同，都是最低的。但是实验2和实验3将证明，在实际运算中FCA所用时间要短于MVHT，同时由于MVHT采用的贪婪搜索算法更容易陷入局部极值，而FCA采用的是随机搜索算法，因而FCA往往可以获得较高的端元提取精度。

实验2 算法时间复杂度实验

该组实验分别改变了端元个数和像元个数，以测试分析5种算法运算时间。时间测试的硬件环境为：Intel（R）Core（TM）i5-3470， CPU3.20 GHz，W indow s 7， MATLAB R 2010a。

（1）改变仿真数据端元个数，以测试不同算法的时间复杂度与端元个数的关系。端元个数p从2增加到10，均是从USGS矿物光谱库中得到，像元个数N为5002， SNR为20 dB。

（2）改变仿真数据像元个数，以测试不同算法时间复杂度与像元个数的关系。像元个数N为1002，2002， 3002， 4002， 5002， 6002， 7002，端元个数p为8，均是从USGS矿物光谱库中得到，SNR为20 dB。

从图4（a）和图4（b）中可看出随着端元个数和像元个数增加，5种算法的运行时间都在增加，其中FCA时间最短，几乎是N-FINDR的百分之一。可见，应用投影思想的FCA的时间复杂度要大大优于其他4种算法。值得注意的是，实验1中MVHT与FCA复杂度几乎相同，而本实验中FCA运行时间却短于MVHT，这主要是因为应用MATLAB实现这两种算法时， MVHT需要用到较多的循环语句实现，而FCA的投影过程使用矩阵运算就可实现。在MATLAB软件中，循环语句的解释时间要长于矩阵运算，故FCA在实际中运算速度快于MVHT。

图3 各种算法的运算复杂度

图4 不同情况下算法运行时间

图4 不同情况下算法运行时间

实验3 算法抗噪声实验

本实验选取光谱信息散度（Spectral Information D ivergence， SID）［16］和均方根误差（Root M ean Square Error， RMSE）［3］两个指标来衡量不同噪声强度下，各算法提取端元的性能。SID用以衡量提取的端元和原端元的相似程度，取其均值作为性能评价的标准，即。RMSE用以衡量原图像和反混图像的相似程度，取其均值作为性能评价标准，即。SID和RMSE的值越低，所比较的两个对象越相似。

实验中，高斯白噪声的强度SNR分别为：10 dB，15 dB， 20 dB， 25 dB， 35 dB， 40 dB。端元个数p为5，像元个数N为2002。为避免误差，对仿真数据进行50次实验，取其均值作为最终结果。

表1 不同噪声强度下算法RMSE（×10-2）性能比较

表2 不同噪声强度下算法SID（×10-3）性能比较

如表1和表2所示，随着SNR降低，SID和RMSE逐渐增加，即5种算法的精度都在下降。同时，在不同信噪比条件下，FCA的SID和RMSE值均低于MVHT和SGA， 与N-FINDR相近。说明SGA和MVHT所采用的逐一确定端元，并且确定后不再更改的方法不能很好地从整体把握单形体体积最大的思想，所选择的点有可能为异常极值点。而以FCA为代表的算法，所选择的端元在运算过程中，随时可被能构成更大体积的像元代替，更能从整体考虑特征空间中像元的分布情况，更具有灵活性，因此算法精度也更好。

4.2 实际数据实验

本节使用的数据是美国Nevada州南部沙漠地区的Cuprite数据。它是由机载可见光及红外成像光谱仪（A irborne V isib le/ In frared Im aging Spectrom eter， AV IRIS）拍摄于1997年6月19日。该数据共有224个波段，波长范围为370 nm到2500 nm，光谱分辨率为10 nm。截取原始数据中一块大小为250× 191比较典型的数据用于算法实现［6］。其中波段1～2， 104～113， 148～167， 221～224共36个波段因为信噪比太低或为水吸收波段而被移除，余下的188个波段用于算法验证。虚拟维度VD估计的端元个数及相应的虚警率pf如表3所示，虚警率pf越低，被错判为纯像元的可能性越低。但依据该地区地物的实际分布［17］，和相关文献资料［10，18］，以及基于不遗漏端元的考虑，本实验将端元个数设定为18。分别用SGA， MVHT， NVFS， N-FINDR，FCA这5种算法提取出该图像的18个端元进行丰度反演并重新混合图像，图5为用FCA提取的该地区几种典型矿物端元的丰度反演结果。利用丰度反演的结果初步确定相应地物，并与USGS矿物光谱库比较，结果如表4。表中数值表示光谱相似度，?是光谱角度匹配（Spectral Angle M apper， SAM）和光谱特征拟合（Spectral Feature Fitting， SFF）的加权平均值，权重分别为0.6和0.4。此数值越接近1，说明两个光谱曲线越相似。表4的最后给出了5种算法求得的光谱相似度的均方根。

图5 FCA提取的几种典型矿物端元的丰度反演结果

表3 虚拟维度VD估计的端元个数及相应的虚警率

表4 不同算法提取出的端元光谱与USGS光谱库中光谱的相似度比较

表5 算法运行时间（s）

表4中的数据表明用FCA算法提取出的端元光谱与USGS光谱库中光谱较为相似，其结果优于其他几种算法。同时表5列出的算法运行时间表明，FCA算法运行速度最快，是N-FINDR的百倍多。综上，实际数据表明，本论文提出的FCA算法相对于其他相似算法而言，具有较好的端元提取精度和较快的算法运行时间，可以快速准确地实现高光谱图像端元提取。

5 结束语

本文提出一种基于代数余子式的N-FINDR快速端元提取算法，简称FCA。该算法通过寻找在特征空间中最大体积的单形体来提取端元。在计算过程中，FCA可得到一个代数余子式构成的向量，并能将所有像元投影到这一向量上，从而将复杂度极高的行列式计算转化为复杂度较低的内积（投影）运算，进而得到具有最大体积的单形体。由于FCA将循环迭代计算体积的方法改进为计算投影，故FCA具有投影类端元提取算法计算复杂度低的优势。同时，在端元搜索方面，不同于SGA和MVHT所采用的逐一确定端元并且端元不可更改的方法，FCA所确定的端元都可在下一次迭代过程中被纯度更高的像元替换，故降低了提取出异常极值点的概率，并且保证所确定的端元都是特征空间中由全部像元构成的最大体积的单形体顶点。这也是FCA算法保持良好精度的原因。综上，经过仿真和真实实验验证，FCA算法可以在计算复杂度很低的情况下保持很好的端元提取精度，这对于大数据量的高光谱图像端元提取的快速处理很有意义。

［1］ Sun Kang， Geng Xiu-rui， Wang Pan-shi， et al.. A fast endmember extraction algorithm based on gram determinant［J］. IEEE Geoscience and Remote Sensing Letters，2014， 11（6）: 1124-1128.

［2］ Ji Lu-yan， Geng Xiu-rui， Yu Kai， et al.. A new non-negative matrix factorization method based on barycentriccoordinates for endmember extraction in hyperspectral remote sensing［J］. International Journal of Rem ote Sensing，2013， 34（19）: 6577-6586.

［3］ Deng Cheng-bin and Wu Chang-shan. A spatially adaptive spectral m ixture analysis for mapping subpiexl urban im pervious surface distribution［J］. Rem ote Sensing of Environment， 2013， 133（1）: 62-70.

［4］ Beauchem in M. A m ethod based on spatial and spectral information to reduce the solution space in endmember extraction algorithm s［J］. Rem ote Sensing Letters， 2014， 5（5）: 471-480.

［5］ Boardm an J， K ruse F， and G reen R. M app ing target signatures via partial unmixing of AVIRIS data［C］. Proceedings of the Summaries of the VI Jet Propulsion Laboratory Airborne Earth Science Workshop， Pasadena，USA， 1995: 23-26.

［6］ Nascimento P and Bioucas D. Vertex com ponent analysis: a fast algorithm to unm ix hyperspectral data［J］. IEEE Transactions on Geoscience and Remote Sensing， 2005， 43（4）: 898-910.

［7］ Chang Chein-I. Hyperspectral Imaging: Techniques for Spectral Detection and Classification［M］. New York: K luwer Academ ic/ Plenum Pub lishers， 2003: 73-88.

［8］ W inter M. N-FINDR: an algorithm for fast autonom ous spectral endmember determ ination in hyperspectral data［C］. Proceedings of the SPIE Imaging Spectrometry V， San Diego，USA， 1999: 266-275.

［9］ Chang Chein-I， Wu Chao-cheng， Liu W ei-m in， et al.. A new grow ing method for sim plex-based endmember extraction algorithm［J］. IEEE Transactions on Geoscience and Rem ote Sensing， 2006， 44（10）: 2804-2819.

［10］ Liu Jun-m in and Zhang Jiang-she. A new m aximum sim p lex volume method based on householder transformation for endmember extraction［J］. IEEE Transactions on Geoscience and Rem ote Sensing， 2012， 50（1）: 104-118.

［11］ Geng X iu-rui， Zhao Yong-chao， W ang Fu-xiang， et al.. A new volume formula for a sim plex and its application to endmember extraction for hyperspectral image analysis［J］. Internation Journal of Rem ote Sensing， 2010， 31（4）: 1027-1035.

［12］ Shi Chen and Wang Le. Incorporation spatial in formation in spectral unm ixing: a review［J］. Rem ote Sensing of Environment， 2014， 149（2）: 70-87.

［13］ Chang Chein-I and Du Qian. Estimation of number of spectrally distinct signal sources in hyperspectral imagery［J］. IEEE Transactions on Geoscience and Rem ote Sensing， 2004，42（3）: 608-619.

［14］ Zhang Xian-da. Matrix Analysis and Applications［M］. Beijing: Tsinghua University Press， 2004: 70-81.

［15］ Heinz D and Chang Chein-I. Fully constrained least squares linear spectral m ixture analysis method for material quantification in hyperspectral imagery［J］. IEEE Transactions on Geoscience and Rem ote Sensing， 2001， 39（3）: 529-545.

［16］ Chang Chein-I. An in formation-theoretic approach to spectral variability， sim ilarity， and discrim ination for hyperspectral image analysis［J］. IEEE Transactions on Information Theory， 2000， 46（5）: 1927-1932.

［17］ Roger N， Clark R， Gregg A， et al.. Evolution in imaging spectroscopy analysis and sensor signal-to-noise: an exam ination of how far we have come［C］. Proceedings of the Summ aries of the 6th Annual JPL A irborne Geoscience Workshop， Pasadena， USA， 1996: 4-8.

［18］ Geng Xiu-rui， Xiao Zheng-qing， Ji Lu-yan， et al.. A Gaussian elim ination based fast endmember extraction algorithm for hyperspectral im agery［J］. ISPRS Journal of Photogrammetry and Rem ote Sensing， 2013， 7（9）: 211-218.

李 琳： 女，1989年生，博士生，研究方向为高/多光谱数据优化和混合像元分析.

孟令博： 女，1989年生，博士生，研究方向为高光谱遥感图像处理和水光谱学.

孙 康： 男，1988年生，博士生，研究方向为高光谱特征提取与高性能计算.

赵永超： 男，1971年生，研究员，研究方向为高光谱遥感图像去噪、大气纠正和辐射定标.

A Fast N-FINDR Algorithm Based on Cofactor of a Determ inant

Li Lin①②Meng Ling-bo①②Sun Kang①②Zhao Yong-chao①①（Institute of Electronics， Chinese Academ y of Sciences， Beijing 100190， China）
②（University of Chinese Academy of Sciences， Beijing 100049， China）

Endmember extraction methods based on geometric distribution of hyperspectral images usually divide into pro jection algorithm and the maximum volume formula for sim plex， which the former has lower com putational com p lexity and the latter has better precision. A Fast endm ember extraction m ethod based on Cofactor of a determ inant A lgorithm （FCA） is p roposed. The algorithm combines the two kinds of algorithm s， and which means it has a high speed and accuracy performance for endmember extraction. FCA finds the max volume of simp lex by making pixels p roject to vectors， which are com posed of the cofactors of elements in endmember determ inant. Besides， FCA is flexible in endmember search， for it can use higher purity pixels to rep lace the endm embers extracted in the last iteration， which ensures that all the endmembers extracted by FCA are the vertices of simp lex. The theoretical analysis and experiments on both simulated and real hyperspectral data demonstrate that the proposed algorithm is a fast and accu rate algorithm for endm em ber extraction.

Image p rocessing； Hyperspectral； Endmember extraction； Simp lex； Maximum volume； Cofactor；Pro jection

TP751.1

： A

：1009-5896（2015）05-1128-07

10.11999/JEIT140923

2014-07-15 收到，2014-11-28改回

*通信作者：李琳 linli_iecas8@163.com