趙春暉，崔士玲，趙艮平，鐘衛(wèi)

（哈爾濱工程大學(xué)信息與通信工程學(xué)院，黑龍江哈爾濱150001）

一種改進(jìn)迭代誤差分析端元提取算法

趙春暉，崔士玲，趙艮平，鐘衛(wèi)

（哈爾濱工程大學(xué)信息與通信工程學(xué)院，黑龍江哈爾濱150001）

迭代誤差分析（IEA）算法是應(yīng)用比較廣泛的端元提取算法之一，針對(duì)IEA端元提取算法計(jì)算量大的缺點(diǎn)，從減少參與迭代過(guò)程中的像元數(shù)目進(jìn)行改進(jìn)。根據(jù)凸面幾何理論，混合像元位于其端元構(gòu)成子空間內(nèi)部，這部分像元到其端元正交子空間（OSP）投影值理論上為零，在迭代求下一個(gè)端元過(guò)程中，可以將這部分無(wú)用像元去除，從而減少每步迭代過(guò)程像元數(shù)目。采用模擬數(shù)據(jù)和真實(shí)高光譜數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，證明改進(jìn)算法與原算法提取端元精度相同，隨著端元提取個(gè)數(shù)的增多，參與迭代過(guò)程的像元數(shù)目逐次減少，比原始的IEA端元提取算法減少了計(jì)算時(shí)間。

端元提取；迭代誤差分析；正交子空間投影；最小二乘算法

近年來(lái)，高光譜遙感技術(shù)快速發(fā)展［1］。目前，高光譜遙感已經(jīng)在農(nóng)業(yè)生產(chǎn)、環(huán)境監(jiān)測(cè)、地質(zhì)勘探、軍事對(duì)抗等多個(gè)領(lǐng)域得到成功應(yīng)用［2］。由于地面的復(fù)雜多樣性及傳感器空間分辨率的限制，導(dǎo)致混合像元在高光譜圖像上廣泛存在［3］。因此，如何有效地解決混合像元是高光譜圖像分析的難題之一。

在混合像元分解中，端元提取是關(guān)鍵。自動(dòng)端元提取是高光譜圖像端元提取面臨的一大挑戰(zhàn)，在過(guò)去的幾十年里已經(jīng)提出了很多端元提取算法。有代表性的方法包括像素純度指數(shù)（PPI）［4］、光學(xué)實(shí)時(shí)自適應(yīng)光譜識(shí)別系統(tǒng)（ORASIS）［5］、迭代誤差分析（IEA）［6］、內(nèi)部最大體積（N?FINDR）［7］和自動(dòng)形態(tài)學(xué)端元提?。ˋMEE）［8］，等。IEA是一種不需要對(duì)原始數(shù)據(jù)進(jìn)行降維或者去冗余而直接對(duì)數(shù)據(jù)處理的端元提取算法，以誤差大小作為判斷端元的標(biāo)準(zhǔn)［3］。已經(jīng)證明IEA端元提取算法是具有較好魯棒性的廣泛應(yīng)用的端元提取算法之一。由于IEA是順序迭代求取每個(gè)端元，隨著端元數(shù)目的增大，它的計(jì)算量也隨之增大。在迭代的過(guò)程中，并不是所有的混合像元都對(duì)端元的提取有幫助。針對(duì)此問(wèn)題，本文在深入分析最小二乘算法與正交子空間投影之間關(guān)系的基礎(chǔ)上，發(fā)現(xiàn)由已經(jīng)提取端元構(gòu)成凸面單形體內(nèi)部的像元到已提取端元構(gòu)成的正交子空間的投影理論上為零，這部分像元對(duì)下一個(gè)端元提取無(wú)用，可以剔除，以此降低每步迭代計(jì)算量。

1 迭代誤差分析算法

迭代誤差分析（iterative error analysis，IEA）以線性光譜混合模型的代數(shù)學(xué)描述為基礎(chǔ)，是一種不需要對(duì)原始數(shù)據(jù)進(jìn)行降維或者去冗余而直接對(duì)數(shù)據(jù)處理的端元提取算法，以誤差大小作為判斷端元的標(biāo)準(zhǔn)［3］。線性光譜混合模型［9］可以表示為

式中：y為混合像元，E是端元，a是混合系數(shù)，M是端元的數(shù)量，N表示噪聲影響。

該模型有2個(gè)約束條件：

1）aj≥0非負(fù)約束；

迭代誤差分析端元［3］提取算法首先給定一個(gè)初始向量，然后采用全約束最小二乘方法對(duì)混合像元y進(jìn)行解混，每次解混操作得到式（1）的一個(gè)估計(jì)a＾，利用已有的端元E和a＾可得到y(tǒng)的一個(gè)估計(jì)y＾，那么誤差可表示為

對(duì)應(yīng)的y'為誤差向量，產(chǎn)生最長(zhǎng)誤差向量的像元記為新的端元，加入端元集合并開(kāi)始下一次迭代，直到找到M個(gè)端元或誤差向量長(zhǎng)度足夠小。需要注意的是，初始向量只參與第1個(gè)端元的選擇，得到第1個(gè)端元后，初始向量不再參加后面端元提取的運(yùn)算。

有時(shí)，為了減小噪聲和誤差，在每次迭代中允許選擇誤差向量最長(zhǎng)且光譜角足夠小的若干像元取平均值，得到新的端元。

2 改進(jìn)的IEA端元提取算法

由第1部分IEA端元提取算法可以看出，最耗時(shí)的部分在于每次迭代過(guò)程誤差向量模值的計(jì)算，如果能夠減少每次迭代過(guò)程中的像元數(shù)目，那么計(jì)算量就可以降下來(lái)。誤差計(jì)算過(guò)程與正交子空間投影之間存在如下聯(lián)系：

式中：I為與PE同階的單位矩陣，為端元E的正交子空間。

對(duì)于傳統(tǒng)線性光譜混合模型，y可以表示成不感興趣部分（已經(jīng)提出端元以及由它們構(gòu)成的混合像元）與感興趣部分（剔除不感興趣剩余的部分）的混合：

式中：d為感興趣部分，αp為d所對(duì)應(yīng)的豐度，U為不感興趣部分，γ為U所對(duì)應(yīng)的豐度，N為噪聲。

通過(guò)均方根誤差的計(jì)算式（3）可以分為2部分計(jì)算，先計(jì)算前一部分yTP⊥E，然后再和混合像元y相乘，求最大值得到端元。假如已經(jīng)提取出m個(gè)端元EM＝(e1，e2，…eM)，構(gòu)造它們的正交子空間：

由式（7）可以看出，經(jīng)過(guò)將混合像元投影到已提取端元正交子空間上，背景信號(hào)（由已提取端元構(gòu)成凸面單形體內(nèi)部的混合像元）被消除，原始噪聲也得到了投影壓縮。

為了更直觀的說(shuō)明去除像元的過(guò)程，用提取出3個(gè)端元的幾何關(guān)系為例進(jìn)行說(shuō)明，如圖1所示，3個(gè)實(shí)黑點(diǎn)代表已提取的3個(gè)端元e1、e2和e3，其中圓圈代表由e1、e2和e3構(gòu)成單形體內(nèi)部混合像元，方形代表其他像元，圓圈代表的像元到端元e1、e2和e3正交子空間的投影為零，可以在下一步迭代過(guò)程中去掉，新的端元來(lái)源于方形代表的像元。

圖1 3個(gè)端元混合像元去除示意圖Fig.1 Mixture of three endmembers illustrating the pixels removal

改進(jìn)IEA端元提取算法步驟：

1）給定一個(gè)初始向量（第1個(gè)初始端元選擇的是距離原點(diǎn)最遠(yuǎn)的若干像元的均值），并且確定端元個(gè)數(shù)；

2）用得到的E計(jì)算均方根誤差，對(duì)圖像數(shù)據(jù)中每個(gè)向量通過(guò)式（3）求取均方根誤差，得到誤差矩陣；

3）選擇誤差矩陣中模值最大的點(diǎn)，為了減少噪聲和誤差，在每次迭代的過(guò)程中選擇誤差向量最長(zhǎng)且光譜角足夠小的若干個(gè)像元取平均值，將得到的平均值作為第1個(gè)端元e1；

4）求端元e1的正交子空間，將所有的像元投影到e1的正交子空間上，將投影值小于閾值ε的像元去掉，得到新的觀測(cè)值為y1；

5）用e1對(duì)y1中每個(gè)向量通過(guò)式（3）求取均方根誤差，為了提高端元提取的準(zhǔn)確性，可以選擇全約束最小二乘算法對(duì)豐度進(jìn)行估計(jì)，將誤差最大且光譜角足夠小的若干個(gè)像元取平均值作為第2個(gè)端元e2；

重復(fù)上面步驟，直至找到規(guī)定數(shù)量的端元為止。

3 仿真實(shí)驗(yàn)

式中：〈x·y〉是2個(gè)向量x和y的內(nèi)積。

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

仿真高光譜數(shù)據(jù)，從美國(guó)地質(zhì)調(diào)查局（USGS）礦物光譜庫(kù)選取8種地物光譜，分別為明礬石、銨長(zhǎng)石、方解石、高嶺石、白云母、榍石、黃鉀鐵礬和綠脫石，用這8種地物合成40 000個(gè)具有224波段的高光譜數(shù)據(jù)，每個(gè)像元由8種純端元中的若干種隨機(jī)混合而成，并添加8種地物的純光譜，共40 008個(gè)像元，為模擬真

IEA端元提取算法是經(jīng)典的端元提取算法之一，在文獻(xiàn)［10?11］已經(jīng)和其他端元提取算法進(jìn)行了比較，本文對(duì)改進(jìn)的算法與原始的IEA算法進(jìn)行對(duì)比。為了驗(yàn)證改進(jìn)算法的有效性，分別用合成模擬數(shù)據(jù)和真實(shí)的高光譜數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。計(jì)算機(jī)的硬件配置為：處理器為Intel Core i3，主頻為2.13 GHz，內(nèi)存為DDR32 GB，仿真軟件為Matlab 7.11版本。

為了評(píng)價(jià)端元提取的精度，采用光譜角距離作為評(píng)價(jià)指標(biāo)，光譜角距離［12］定義為實(shí)的高光譜影像，附加信噪比為25 dB的高斯噪聲。

為了驗(yàn)證算法的有效性，做2組實(shí)驗(yàn)，第1組實(shí)驗(yàn)驗(yàn)證改進(jìn)算法提取端元光譜與原始算法提取端元相同，并且精度相同，說(shuō)明去除部分無(wú)用像元并不影響端元提取的精度；第2組實(shí)驗(yàn)驗(yàn)證在保證端元提取精度前提下，改進(jìn)算法隨著端元提取個(gè)數(shù)的增多，在每次迭代過(guò)程中減少的像元數(shù)目和減少像元百分比與投影閾值之間的關(guān)系以及改進(jìn)算法比IEA算法減少的時(shí)間。

改進(jìn)算法和原始算法都能夠準(zhǔn)確的提取出8種端元，為了更客觀的評(píng)價(jià)提取出端元的準(zhǔn)確度，表1給出改進(jìn)算法與原始算法端元提取精度對(duì)比，可以看到，光譜角都很小，提取出的效果很好，并且2種算法提取出端元的精度相同，說(shuō)明改進(jìn)的算法并不影響端元提取的精度。

表2是在保證端元提取精度基礎(chǔ)上，不同閾值下的每步迭代減少的像元個(gè)數(shù)及占總像元的百分比。橫向來(lái)看，隨著提取端元個(gè)數(shù)的增多，每步減少的像元數(shù)目逐次增多，在閾值設(shè)置從0.016～0.02時(shí)，第7次迭代減少的像元數(shù)目達(dá)到90%以上，證明該算法能夠大大減少運(yùn)算量?？v向來(lái)看，隨著投影閾值的變化，每步迭代減少的像元數(shù)目是不同的，當(dāng)閾值為0.020時(shí)，達(dá)到理想的效果，計(jì)算時(shí)間比較小。合理設(shè)置投影閾值對(duì)計(jì)算量降低十分關(guān)鍵。

表1 改進(jìn)算法與原始IEA提取端元的精度對(duì)比Table1 Comparing the precision of improved endmember extraction algorithm with the original IEA

表2 不同閾值下的每步迭代減少的像元個(gè)數(shù)及占總像元的百分比Table 2 The number of removal pixels and the percent of total pixels per iteration under different thresholds

表3是改進(jìn)算法在不同閾值下與原始IEA算法運(yùn)行20次平均時(shí)間對(duì)比，改進(jìn)算法隨著投影閾值的增大，計(jì)算時(shí)間逐步降低，當(dāng)閾值為0.02時(shí)，比原始算法計(jì)算時(shí)間降低了約50%，大大減小了計(jì)算時(shí)間，但是當(dāng)閾值降低到0.012時(shí)，計(jì)算速度反而比原來(lái)還要慢，原因在于每次迭代減少的像元數(shù)目急速減少，而閾值的判定也需要花費(fèi)時(shí)間，因此合理有效地設(shè)置閾值可以比原始IEA大大降低計(jì)算時(shí)間。

表3 不同閾值下改進(jìn)算法與原始IEA算法運(yùn)行20次平均時(shí)間對(duì)比Table3 Average time comparison of improved algorithm with the original IEA algorithm running 20 times under different thresholds s

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

真實(shí)的高光譜數(shù)據(jù)采用Cuprite礦物數(shù)據(jù)，該影像由AVIRIS傳感器獲得，而且含有該地區(qū)對(duì)應(yīng)的地質(zhì)調(diào)查圖。實(shí)驗(yàn)所用數(shù)據(jù)大小為250×191像元，光譜范圍為0.369～2.480 μm，共224波段，去除因水汽吸收和傳感器噪聲等低信噪比波段（1～2、104～113、148～167、221～224），實(shí)驗(yàn)用的數(shù)據(jù)共188波段。

從Cuprite礦物數(shù)據(jù)中提取出6種端元，分別為沙漠地表、玉髓、綠脫石、明礬石、榍石和蒙脫石，提取出的端元在文獻(xiàn)［13］中可以查到。由模擬實(shí)驗(yàn)得出：設(shè)置合理的閾值可以在保證提取精度的基礎(chǔ)上，最大限度減少計(jì)算量。在真實(shí)數(shù)據(jù)中設(shè)置的閾值為0.002，在此閾值下可以獲得良好的效果，提取6種端元，總共迭代5步，在保證端元提取精度的前提下，每次迭代中減少的端元數(shù)目與所占的百分比如表4所示，可以看出，隨著端元提取個(gè)數(shù)的增多，每次迭代減少的像元數(shù)逐次下降。改進(jìn)算法與原始算法提取出的端元光譜相同，將提取出的端元光譜與USGS光譜庫(kù)進(jìn)行對(duì)比，如圖2所示。

表4 每步迭代減少的像元個(gè)數(shù)及占總像元的百分比Table4The number of removal pixels and the per?cent of total pixels per iteration

圖2 提取端元光譜與真實(shí)光譜對(duì)比Fig.2 Extraction endmember spectral compared with the real spectral

從圖中可以看出，提取出的端元光譜和真實(shí)光譜基本吻合。為了更客觀地評(píng)價(jià)提取端元結(jié)果的精確性，采用光譜角距離進(jìn)行評(píng)價(jià)，光譜角越小，精度越高，如表5所示，可以看出改進(jìn)算法與原始IEA端元提取算法提取端元精度相同，但是改進(jìn)算法所用的時(shí)間明顯減少，表5最后一行給出的是運(yùn)行20次的平均時(shí)間，改進(jìn)算法所用時(shí)間比原始算法大約減少了25%。

表5 改進(jìn)算法與原始IEA提取端元的精度及時(shí)間對(duì)比Table5 Comparing the precision and time of im?provedendmemberextractionalgorithm with the original IEA

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

本文提出了一種改進(jìn)的IEA端元提取算法，通過(guò)將像元到已提取端元正交子空間投影值，去除已提取端元構(gòu)成單形體內(nèi)部的像元，從而減少每次迭代過(guò)程中像元的個(gè)數(shù)，提高了算法的效率。實(shí)驗(yàn)結(jié)果證明，改進(jìn)算法在保證端元提取精度的基礎(chǔ)上，大大減少了參與每次迭代的像元數(shù)目，整體計(jì)算時(shí)間降低。由于噪聲等的影響，投影閾值的設(shè)定需要根據(jù)實(shí)際情況設(shè)定，如何自動(dòng)的選擇閾值是下一步研究的主要問(wèn)題。

［1］王立國(guó)，魏芳潔.結(jié)合遺傳算法和蟻群算法的高光譜圖像波段選擇［J］.中國(guó)圖象圖形學(xué)報(bào)，2013，18（2）：235?242.

WANG Liguo，WEI Fangjie.Band selection for hyperspectral imagery based on combination of genetic algorithm and ant colony algorithm［J］.Journal of Image and Graphics，2013，18（2）：235?242.

［2］王立國(guó)，趙春暉.高光譜圖像處理技術(shù)［M］.北京：國(guó)防工業(yè)出版社，2013：18?29.

［3］張兵，高連如.高光譜圖像分類與目標(biāo)探測(cè)［M］.北京：科學(xué)出版社，2011：30?55.

［4］BOARDMAN J W.Automating spectral unmixing of AVIRIS data using convex geometry concepts［J］.Summaries of the Fourth Annual JPL Airborne Geoscience Workshop，JPL Publication，1993，1：11?14.

［5］BOWLES J，ANTONIADES J，BAUMBACK M，et al.Real time analysis of hyperspectral data sets using NRL's ORASIS algorithm［J］.Proceedings of the SPIE，1997，3118：38?45.［6］NEVILLE R A，STAENZ K，SZEREDI T，et al.Automatic endmember extraction from hyperspectral data for mineral ex?ploration［C］／／Proceedings 21st Canadian Symposium Re?mote Sensing，Ottawa，Canada，1999：21?24.

［7］趙春暉，成寶芝，楊偉超.利用約束非負(fù)矩陣分解的高光譜解混算法［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2012，33（3）：377?382.

ZHAO Chunhui，CHENG Baozhi，YANG Weichao.Algo?rithm for hyperspectral unmixing using constrained nonnega?tive matrix factorization［J］.Journal of Harbin Engineering University，2012，33（3）：377?382.

［8］PLAZA A，MARTINEZ P，PEREZ R，et al.Spatial／spec?tral endmember extraction by multidimensional morphological operations［J］.IEEE Transactions on Geoscience and Re?mote Sensing，2002，40：2025?2041.

［9］王玉磊，趙春暉，齊濱.基于光譜相似度量的高光譜圖像異常檢測(cè)算法［J］.吉林大學(xué)學(xué)報(bào)：工學(xué)版，2013，43：148?153.

WANG Yulei，ZHAO Chunhui，QI Bin.Hyperspectral anom?aly detection algorithm based on spectral similarity scale［J］.Journal of Jilin University：Engineering and Technology Edi?tion，2013，43：148?153.

［10］PLAZA A，MARTINEZ P，PEREZ R，et al.A quantitative and comparative analysis of endmember extraction algo?rithms from hyperspectral data［J］.IEEE Transactions on Geoscience and Remote Sensing，2004，42：650?663.

［11］WINTER M E，WINTER E M.Comparison of approaches for determining end?members in hyperspectral data［J］.Proceedings IEEE Aerospace Conference.Big Sky，USA，2000，3：18?25.

［12］丁海勇，史文中.利用卡方分布改進(jìn)N?FINDR端元提取算法［J］.遙感學(xué)報(bào)，2013，17（1）：130?137.

DING Haiyong，SHI Wenzhong.Fast N?FINDR algorithm for endmember extraction based on Chi?square distribution［J］.Journal of Remote Sensing，2013，17（1）：130?137.

［13］JOSé M P，JOSé M N，BIOUCAS D.Vertex component a?nalysis：a fast algorithm to unmix hyperspectral data［J］.IEEE Transactions on Geoscience and Remote Sensing，2005，43：898?910.

Endmember extraction algorithm based on improved iterative error analysis

ZHAO Chunhui，CUI Shiling，ZHAO Genping，ZHONG Wei
（College of Information and Communications Engineering，Harbin Engineering University，Harbin 150001，China）

Iterative error analysis（IEA）algorithm is one of the widely used endmember extraction algorithms.As for the compute-intensive shortcoming，this paper improves it from reducing the number of pixels participated in the iterative process.According to the convex geometry theory，mixed pixel locates on its internal subspace constituted by its endmembers.The projection value of these cells to their endmember orthogonal subspace（OSP）is zero，which can be removed during the iterative process for deriving the next endmember，thereby reducing the number of pixels in every iteration step.The experiments with simulated data and real hyperspectral data were done，proving that the improved algorithm has the same precision of the extracted endmember with original algorithm，and with the increasing of the number of extracted endmembers，the number of pixels participated in the iteration decreases grad?ually，which reduces the computation time than the endmember extraction of the original IEA algorithm.

endmember extraction；iterative error analysis；orthogonal subspace projection；least?squares algorithm

10.3969／j.issn.1006?7043.201310094

http：／／www.cnki.net／kcms／doi／10.3969／j.issn.1006?7043.201310094.html

TP751.1

A

1006?7043（2015）02?0257?05

2013?10?31.網(wǎng)絡(luò)出版時(shí)間：2014?11?27.

國(guó)家自然科學(xué)基金資助項(xiàng)目（61405041，61275010），黑龍江省自然科學(xué)基金重點(diǎn)資助項(xiàng)目（ZD201216），哈爾濱市優(yōu)秀學(xué)科帶頭人基金資助項(xiàng)目（RC2013XK009003）.

趙春暉（1965?），男，教授，博士生導(dǎo)師.

趙春暉，E?mail：zhaochunhui＠hrbeu.edu.cn