基于代數余子式的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