基于Gram行列式的快速端元提取方法

孫 康，帥 通，陳金勇

(1.中國(guó)電子科技集團(tuán)公司第五十四研究所，河北 石家莊 050081；2.武漢大學(xué) 測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079)

基于Gram行列式的快速端元提取方法

孫 康1，2，帥 通1，陳金勇1

(1.中國(guó)電子科技集團(tuán)公司第五十四研究所，河北 石家莊 050081；2.武漢大學(xué) 測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079)

高光譜圖像端元提取往往涉及到高維空間中單形體體積的計(jì)算，使用無(wú)需降維的體積公式能夠避免信息損失，但卻具有極大的計(jì)算復(fù)雜度。針對(duì)這一缺點(diǎn)進(jìn)行了研究，提出了基于Gram行列式快速的端元提取算法。該算法不需要計(jì)算單形體體積，而是利用了體積公式的遞推關(guān)系，大大降低了計(jì)算復(fù)雜度。模擬和真實(shí)數(shù)據(jù)試驗(yàn)表明，該算法在保證高精度端元提取的同時(shí)，具有極快的端元提取速度。

高光譜圖像；端元提取；單形體；線性混合模型；Gram行列式

0 引言

混合像元分析是高光譜遙感中信息提取的重要手段之一，端元提取(Endmember Extraction)是混合像元分析至關(guān)重要的步驟。現(xiàn)有的端元提取方法可以分為端元選擇和端元生成2大類[1]。常見(jiàn)的端元選擇算法包括N-FINDR[2]、VCA[3]和GEM[4]等。這類算法的主要優(yōu)點(diǎn)是計(jì)算量小，但是需要純像元假設(shè)。與端元生成方法(如MVC-NMF[5]和GOM[6])相比，端元選擇往往具有較低的計(jì)算復(fù)雜度，并且獲得的端元具有實(shí)際物理意義，因此得到了廣泛研究和應(yīng)用。

端元選擇的最重要的依據(jù)之一是端元構(gòu)成的單行體體積最大。基于最大體積的原則，發(fā)展了大量的算法如N-FINDR、SGA[7]、OBA[8]和MVHT[9]。然而對(duì)于這些基于體積的算法而言，有些需要降維，如N-FINDR、SGA以及MVHT，這會(huì)造成一定的信息損失。

耿修瑞等提出了一種不需要降維的體積公式[10]，使用該體積公式進(jìn)行端元選擇得到的結(jié)果優(yōu)于需要降維的算法，本文中，該算法稱為New Volume Formula for Simplex(NVFS)。由于不需要降維，這個(gè)算法能夠有效地提取低概率目標(biāo)的端元，但是一個(gè)嚴(yán)重的缺點(diǎn)是計(jì)算太慢。本文提出了一種基于Gram行列式的快速端元選擇算法——Gram Determinant Based Algorithm(GDA)，使得端元提取的計(jì)算復(fù)雜度大大降低，同時(shí)具有更優(yōu)的端元提取結(jié)果。

1 端元提取背景

1.1 線性混合模型原理

1.2 高維空間單形體體積

對(duì)于p個(gè)端元，記為ei(i=1,…,p)，為了計(jì)算端元所構(gòu)成的體積，一個(gè)經(jīng)常使用的體積公式如下：

(1)

耿修瑞在文獻(xiàn)[10]中給出了一個(gè)不需要降維的體積公式：

(2)

從式(2)可以看出，由于ATpAp始終都是方陣，因此不需要降維也可計(jì)算體積。Ap中所有端元都減去e1，表示將數(shù)據(jù)移至以e1為原點(diǎn)。NVFS使用的就是如式(2)所示的體積公式，在進(jìn)行端元提取時(shí)，NVFS采用了跟N-FINDR相同的策略，先隨機(jī)初始化端元，然后逐像元計(jì)算體積，取得最大體積的像元作為新入選端元，直至體積不再增大。NVFS對(duì)于每個(gè)像元都要反復(fù)計(jì)算其與其他端元構(gòu)成的體積，加上該體積表達(dá)無(wú)需降維，因此計(jì)算量極大，限制了其在實(shí)際中的應(yīng)用。同時(shí)受隨機(jī)初始化的影響，計(jì)算結(jié)果并不穩(wěn)定[11]。

2 基于Gram行列式的端元提取

2.1 端元初始化策略改進(jìn)

在進(jìn)行端元提取之前，需要確定端元的個(gè)數(shù)，本文中使用虛擬維數(shù)VD[12]估計(jì)。NVFS先隨機(jī)產(chǎn)生所有端元，然后再逐個(gè)替換，直至不能替換為止，一般要運(yùn)行3～5次才收斂[13]。由于式(2)計(jì)算體積與維數(shù)無(wú)關(guān)，因此可以采取類似SGA的策略，一個(gè)一個(gè)地按順次求取端元，這樣程序只需要運(yùn)行一次即可，而且結(jié)果具有可重復(fù)性。

在進(jìn)行端元初始化的時(shí)候，第1個(gè)端元的選取是基于這樣的事實(shí)：距離原點(diǎn)最遠(yuǎn)的像元必是端元[4]，因此本文選擇亮度最大的像元作為第1個(gè)端元。然后再按照式(2)計(jì)算體積，體積最大的就是第2個(gè)端元，依次進(jìn)行下去，直至所有的端元全部提取為止，這樣只進(jìn)行一次外部循環(huán)即可，同時(shí)結(jié)果具有魯棒性。

2.2 端元得分指標(biāo)設(shè)計(jì)

在使用式(2)作為體積公式提取第k個(gè)端元時(shí)，令Bk=ATkAk，當(dāng)k≥2時(shí)，可以發(fā)現(xiàn)Bk和Bk+1之間有如下遞推的關(guān)系：

(3)

式中，

dk=[(ek+1-e1)T(e2-e1)，(ek+1-e1)T(e3-e1)，…, (ek+1-e1)T(ek-e1)]T;

ak=(ek+1-e1)T(ek+1-e1)。

依據(jù)分塊矩陣行列式計(jì)算規(guī)則[14]矩陣，可以推導(dǎo)出Bk和Bk+1的行列式之間的關(guān)系為：

(4)

2.3 方法計(jì)算復(fù)雜度分析

令p為端元數(shù)，N為像元數(shù)，L為波段數(shù)，則NVFS需要對(duì)所有像元依次計(jì)算體積det(Bp)，而計(jì)算p階方陣的行列式計(jì)算復(fù)雜度為pη，2.3≤η≤2.9[14]，計(jì)算Bp需要(p-1)L次乘法(L為波段數(shù))，因此替換一個(gè)端元所需的計(jì)算量為(pη+pL-L)N(N為像元數(shù)目)。該算法需要替換p個(gè)端元，且一般需要替換3～5次才會(huì)收斂，因此提取p個(gè)端元，NVFS中所需的計(jì)算復(fù)雜度為σ(pη+1+p2L-pL)N，其中2.3≤η≤2.9，3≤σ≤5。

在同樣大小數(shù)據(jù)下，對(duì)于GDA，每次循環(huán)只需計(jì)算每個(gè)像素的端元得分ak-dTkB-1kdk，在一次循環(huán)中，所有像元“共享”ak和B-1k，因此只需計(jì)算一次，然后計(jì)算各個(gè)像元的端元得分指標(biāo)，端元得分指標(biāo)最大的像元就是新的端元。dk的計(jì)算需要L次乘法，N個(gè)像元的第k+1個(gè)端元得分指標(biāo)ak-dTkB-1kdk需要的乘法次數(shù)為(k2+L)N。因此使用GDA提取p個(gè)端元，需要p-1次循環(huán)，而初始化第一個(gè)端元的計(jì)算量為L(zhǎng)*N，故GDA所需的總計(jì)算量為：

可以看出，GDA使得計(jì)算量大大減少。值得一提的是端元得分指標(biāo)的表達(dá)形式為二次型，對(duì)其規(guī)范化以后，可以使用快速的矩陣乘積代替逐像元的循環(huán)計(jì)算。

3 方法性能分析與驗(yàn)證

本文將對(duì)N-FINDR、SGA、MVHT、NVFS、OBA和本文提出的算法這6種算法進(jìn)行比較。試驗(yàn)分為模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)2部分進(jìn)行，分別比較這6種算法在不同信噪比下提取端元的準(zhǔn)確度以及相應(yīng)的計(jì)算時(shí)間，其中端元提取的準(zhǔn)確度采用rmsSAE和rmsSID[3]作為度量指標(biāo)。所有試驗(yàn)均在同一配置的機(jī)器上進(jìn)行，硬件配置為CPU Intel (R) Core(TM) i5-2400@3.10 GHz，內(nèi)存2 GB，系統(tǒng)環(huán)境Windows XP Professional Service Pack 3，使用Matlab 7.0作為軟件編譯環(huán)境。

3.1 模擬數(shù)據(jù)試驗(yàn)

3.1.1 試驗(yàn)1：端元精確度比較

使用上述方法產(chǎn)生了10 000個(gè)像元，并且加入一系列強(qiáng)度的高斯白噪聲，SNR分別為10 dB、15 dB、20 dB、25 dB和30 dB。各個(gè)方法提取端元的rmsSAE和rmsSID結(jié)果如表1所示。

表1 不同方法的端元選擇結(jié)果精度比較

從表1中可以看出，GDA與OBA的端元提取結(jié)果是一樣的，這是因?yàn)槎呤褂玫捏w積公式是相同的，而且選擇初始端元的方法是一樣的。

表1的對(duì)比結(jié)果表明，所有算法的性能都隨著SNR的降低而降低(N-FINDR有幾次例外情況，但多次平均之后也符合這個(gè)規(guī)律)。在SNR較高時(shí)(SNR為30 dB和25 dB時(shí))，使用式(2)的算法，端元提取性能較好，這是因?yàn)樗惴ú恍枰稻S，能夠全面地利用所有光譜信息，因而可以獲得更為準(zhǔn)確的結(jié)果。而需要降維的算法，降維之后信息會(huì)有所損失。值得注意的是，本文改進(jìn)了NVFS的初始化方式及搜尋策略之后，算法性能也有所提高。

3.1.2 試驗(yàn)2：端元得分指標(biāo)隨端元數(shù)的變化

此外，本文還發(fā)現(xiàn)端元指標(biāo)的一個(gè)重要特性，即端元指標(biāo)隨著端元數(shù)目的增加單調(diào)下降。本試驗(yàn)分別用4、5、6、7和8個(gè)實(shí)際端元產(chǎn)生模擬數(shù)據(jù)，并控制噪比為30 dB，然后分別比較了估計(jì)端元數(shù)從1～15時(shí)，端元得分隨端元數(shù)的變化趨勢(shì)，結(jié)果如圖1所示。

圖1 端元得分隨端元數(shù)的增加而單調(diào)下降

從圖1中可以看出，端元得分隨著估計(jì)端元個(gè)數(shù)的增加是單調(diào)下降的。當(dāng)估計(jì)端元數(shù)目小于實(shí)際端元數(shù)目時(shí)，下降的速度很快；當(dāng)估計(jì)端元數(shù)目等于實(shí)際端元數(shù)目時(shí)，該得分急速下降，對(duì)于絕大多數(shù)情況，此時(shí)下降幅度最大；而當(dāng)估計(jì)端元數(shù)目大于實(shí)際端元數(shù)目的時(shí)候，該得分幾乎不變。因此，端元得分可以看作是端元的重要程度，本算法中，最重要的端元最先被找出來(lái)，之后提取的端元重要程度是依次下降的。當(dāng)提取到某個(gè)端元重要程度急速下降時(shí)，則意味著提取的端元已經(jīng)比實(shí)際端元多了一個(gè)，如果此后端元得分指標(biāo)變化異常平緩，則更加說(shuō)明已經(jīng)提取了足夠的端元。端元得分的這個(gè)性質(zhì)可以為估計(jì)端元數(shù)目提供一種輔助性的手段。

3.1.3 試驗(yàn)3：程序運(yùn)行時(shí)間比較

本文還比較了在固定端元個(gè)數(shù)的情況下，計(jì)算時(shí)間隨像元個(gè)數(shù)的變化情況，考慮到實(shí)際端元情況，試驗(yàn)選取p=6，結(jié)果如圖2所示。

從圖2可以看出，在端元個(gè)數(shù)固定為6個(gè)時(shí)，所有算法的計(jì)算時(shí)間隨像元數(shù)都成相同的變化趨勢(shì)，這是因?yàn)檫@些算法的時(shí)間對(duì)于像元個(gè)數(shù)N是同樣的一次函數(shù)。GDA是計(jì)算速度最快的算法，它在提取600*600大小的圖像中的6個(gè)端元用時(shí)僅為0.25 s。

3.2 真實(shí)高光譜數(shù)據(jù)試驗(yàn)

這部分試驗(yàn)中，本文使用實(shí)際高光譜數(shù)據(jù)檢驗(yàn)上述幾種算法的性能。數(shù)據(jù)是由AVIRIS獲取的Nevada Cuprite礦區(qū)圖像。圖像大小為190*250，包含從370～2 500 nm的224個(gè)波段，光譜分辨率大約10 nm，去除了水汽嚴(yán)重吸收以及信噪比較低的1～2、104～113、148～167以及221～224的36個(gè)波段，剩余188個(gè)波段。在使用VD[12]進(jìn)行維數(shù)估計(jì)發(fā)現(xiàn)，端元數(shù)目隨著不同的虛警率在16～20之間。參照其他的試驗(yàn)，本試驗(yàn)令端元數(shù)目p=18。

各個(gè)方法的計(jì)算時(shí)間以及獲得的體積如表2所示。從表2中可以看出，GDA能夠獲得最大的端元體積，由于各個(gè)算法的目標(biāo)函數(shù)是獲得最大體積，因此GDA的性能是最優(yōu)的。此外，GDA的計(jì)算時(shí)間僅為0.607 3 s，具有最高的計(jì)算效率。

表2 不同算法的計(jì)算時(shí)間以及獲得的體積

為進(jìn)一步對(duì)各個(gè)端元精度進(jìn)行比較，本文將各方法提取的典型端元與USGS光譜中對(duì)應(yīng)的光譜進(jìn)行比較，比較的指標(biāo)為光譜角度，結(jié)果如表3所示，其中“-”表示沒(méi)有提取到相應(yīng)端元。

表3 不同算法提取的端元與USGS光譜庫(kù)中光譜的角度(°)

表3中的數(shù)據(jù)表明，使用式(2)的端元提取方法(NVFS、OBA和GDA)獲得的端元提取精度高于使用式(1)的算法(NFINDR、SGA和MVHT)。與NVFS相比，GDA使用順次提取端元的方式提高了端元提取精度。盡管GDA在提取某些端元時(shí)，如(Kaolinite等)，精度比一些算法低了一些，但是GDA提取的端元的精度整體優(yōu)于其他算法。

4 結(jié)束語(yǔ)

端元提取使用無(wú)需降維的體積公式可以避免信息損失，因而能更準(zhǔn)確地選擇所需的端元，然而該算法具有極高的計(jì)算復(fù)雜度。本文提出了一種快速端元選擇算法GDA，在保證提取端元精度的同時(shí)，使速度有了極大的提升。相對(duì)于其他基于體積的端元選擇算法，GDA具有不需要降維、計(jì)算量小、可以提供一個(gè)端元數(shù)據(jù)參考指標(biāo)以及較高魯棒性等優(yōu)點(diǎn)。

[1] BIOUCAS-DIAS J M，PLAZA A，DOBIGEON N，et al.Hyperspectral Unmixing Overview:Geometrical，Statistical，and Sparse Regression-based Approaches[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2012，5(2):354-379.

[2] WINTER M E.N-FINDR:an Algorithm for Fast Autonomous Spectral End-member Determination in Hyperspectral Data[C]∥SPIE’s International Symposium on Optical Science，Engineering，and Instrumentation，1999:266-275.

[3] NASCIMENTO J M，BIOUCAS-DIAS J M.Vertex Component Analysis:a Fast Algorithm to Unmix Hyperspectral Data[J].IEEE Transactions on Geoscience and Remote Sensing，2005，43(4):898-910.

[4] GENG X，XIAO Z，JI L，et al.A Gaussian Elimination Based Fast Endmember Extraction Algorithm for Hyperspectral Imagery[J].ISPRS Journal of Photogrammetry and Remote Sensing，2013(79):211-218.

[5] MIAO L，QI H.Endmember Extraction from Highly Mixed Data Using Minimum Volume Constrained Nonnegative Matrix Factorization[J].IEEE Transactions on Geoscience and Remote Sensing，2007，45(3):765-777.

[6] GENG X，JI L，ZHAO Y，et al.A New Endmember Generation Algorithm Based on a Geometric Optimization Model for Hyperspectral Images[J].IEEE Geoscience and Remote Sensing Letters，2013，10(4):811-815.

[7] CHANG C I，WU C C，LIU W M，et al.A New Growing Method for Simplex-based Endmember Extraction Algorithm[J].IEEE Transactions on Geoscience and Remote Sensing，2006，44(10):2 804-2 819.

[8] TAO X，WANG B，ZHANG L，et al.A New Endmember Extraction Algorithm Based on Orthogonal Bases of Subspace Formed By Endmembers [C]∥Geoscience and Remote Sensing Symposium，2007:2 006-2 009.

[9] LIU J，ZHANG J.A New Maximum Simplex Volume Method Based on Householder Transformation for Endmember Extraction[J].IEEE Transactions on Geoscience and Remote Sensing，2012，50(1):104-118.

[10]GENG X，ZHAO Y，WANG F，et al.A New Volume Formula for a Simplex and Its Application to Endmember Extraction for Hyperspectral Image Analysis[J].International Journal of Remote Sensing，2010，31(4):1 027-1 035.

[11]PLAZA A，CHANG C I.Impact of Initialization on Design of Endmember Extraction Algorithms[J].IEEE Transactions on Geoscience and Remote Sensing，2006，44(11):3 397-3 407.

[12]CHANG C I，DU Q.Estimation of Number of Spectrally Distinct Signal Sources in Hyperspectral Imagery[J].IEEE Transactions on Geoscience and Remote Sensing，2004，42(3):608-619.

[13]PLAZA A，MARTINEZ P，PEREZ R，et al.A Quantitative and Comparative Analysis of Endmember Extraction Algorithms from Hyperspectral Data[J].IEEE Transactions on Geoscience and Remote Sensing，2004，42(3):650-663.

[14]BAUR W，STRASSEN V.The Complexity of Partial Derivatives[J].Theoretical Computer Science，1983(22):317-330.

孫 康 男，(1988—)，博士，工程師。主要研究方向：遙感圖像處理與高性能計(jì)算的研究。

陳金勇 男，(1970—)，研究員。主要研究方向：航天信息處理與航天地面應(yīng)用系統(tǒng)研究。

A New Endmember Extraction Method Based on Gram Determinant

SUN Kang1，2，SHUAI Tong1，CHEN Jin-yong1

(1.The54thResearchInstituteofCETC，ShijiazhuangHebei050081，China;2.StateKeyLaboratoryofInformationEngineeringinSurveying，MappingandRemoteSensing，WuhanUniversity，WuhanHubei430079,China)

The endmember extraction methods for hyperspectral imagery generally involve the computation of simplex volume in high-dimension space.The employment of volume without dimensionality reduction though avoids the information loss，suffers from high computational complexity.This paper puts forward a new endmember extraction method based on Gram determinant which is expected to relieve the efficiency issue.The algorithm does not involve the computation of the volumes，instead，it uses the recurrence relationship of volumes which greatly reduces the computational complexity.The experiments on simulated and real datasets verify the accuracy and efficiency of the proposed method.

hyperspectral imagery;endmember extraction;simplex;LMM;Gram determinant

10.3969/j.issn.1003-3106.2016.11.07

孫 康,帥 通,陳金勇.基于Gram行列式的快速端元提取方法[J].無(wú)線電工程,2016,46(11):26-29，46.

2016-08-22

中國(guó)博士后科學(xué)基金資助項(xiàng)目(2015M580217);河北省博士后科學(xué)基金資助項(xiàng)目(B2015005003)。

TP751

A

1003-3106(2016)11-0026-04