趙春暉, 崔士玲, 趙艮平

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

改進(jìn)的多端元高光譜解混算法

趙春暉, 崔士玲, 趙艮平

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

針對(duì)經(jīng)典多端元光譜混合模型(MESMA)存在著計(jì)算量大，端元預(yù)選繁瑣以及過(guò)擬合等缺點(diǎn)，提出了一種改進(jìn)的多端元解混算法。該算法根據(jù)正交子空間投影具有分離感興趣信號(hào)與不感興趣信號(hào)的特點(diǎn)，將像元投影到全部地物端元(每類地物選擇一條類內(nèi)光譜)構(gòu)成的正交子空間上，按照投影值確定構(gòu)成混合像元每類地物的類內(nèi)光譜，在下一步迭代求解的過(guò)程中，分離出已確定地物類內(nèi)光譜的像元，降低計(jì)算量，然后根據(jù)重構(gòu)誤差變化量確定最優(yōu)端元個(gè)數(shù)，避免過(guò)擬合。實(shí)驗(yàn)結(jié)果表明，改進(jìn)的算法反演豐度誤差和解混時(shí)間都比原有算法降低很多。

高光譜；多端元；光譜解混；正交子空間投影；誤差變化量；類內(nèi)光譜變化

近年來(lái)，高光譜遙感技術(shù)快速發(fā)展，已經(jīng)在農(nóng)業(yè)生產(chǎn)、環(huán)境監(jiān)測(cè)、地質(zhì)勘探、軍事對(duì)抗等多個(gè)領(lǐng)域得到成功應(yīng)用。由于地面的復(fù)雜多樣性及傳感器空間分辨率的限制，使得高光譜圖像普遍存在著混合像元的情況，混合像元的存在嚴(yán)重影響后續(xù)圖像處理的精度[1-2]。

傳統(tǒng)光譜混合模型一個(gè)類別用一個(gè)固定的光譜端元來(lái)代表，由于高光譜圖像空間幅度大導(dǎo)致每類地物內(nèi)的光譜變化一般很大，在這種條件下固定端元光譜很難準(zhǔn)確刻畫一個(gè)類別，導(dǎo)致解混精度不高[3]。Roberts D 等提出多端元光譜混合模型(MESMA)[4]，Asner G 等應(yīng)用蒙特卡羅來(lái)估計(jì)豐度值及其信度區(qū)間來(lái)克服光譜變化(AutoMCU)[5]， B. Somers 根據(jù)類內(nèi)類間光譜變化，提出穩(wěn)定區(qū)域解混(SZU)[6]算法，解混效果較好，但是參數(shù)的變化對(duì)解混精度影響很大。歸一化光譜混合模型 (NSMA)以及基于導(dǎo)數(shù)的光譜解混算法(DSU)，這兩者物理意義不夠明確。其中 MESMA 算法是至今應(yīng)用最為廣泛的一種多端元光譜解混算法[7]，但是其計(jì)算量大，端元預(yù)選繁瑣，而且采用重構(gòu)誤差最小判斷最優(yōu)端元組合，會(huì)導(dǎo)致過(guò)擬合現(xiàn)象，解混時(shí)的采用的端元個(gè)數(shù)要比真實(shí)構(gòu)成混合像元端元個(gè)數(shù)大。針對(duì)上述問(wèn)題，本文提出了一種改進(jìn)的多端元解混算法，有效地解決了這一難題。

1 傳統(tǒng)線性光譜混合模型

傳統(tǒng)線性光譜混合模型[8]可以表示為：

(1)

式中：x為觀測(cè)向量，ej是端元，aj是混合系數(shù)，M是端元的數(shù)量，N代表噪聲影響。

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

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

由最小二乘算法反演豐度，如式(2)，每個(gè)像元由所有的端元E反演，然后重構(gòu)像元如式(3)。

(2)

(3)

在傳統(tǒng)線性光譜混合模型中，一個(gè)類別用一個(gè)固定的光譜端元E來(lái)代表。然而，高光譜圖像空間幅度大導(dǎo)致了類內(nèi)光譜變化一般很大，在這種條件下一個(gè)端元光譜很難準(zhǔn)確刻畫一個(gè)類別，導(dǎo)致光譜解混精度不高，改進(jìn)多端元光譜混合模型算法可以很好地解決這個(gè)問(wèn)題。

2 正交子空間投影算法

對(duì)于傳統(tǒng)線性光譜混合模型，觀測(cè)向量x可以表示成感興趣部分與不感興趣部分的混合[9]，再疊加上噪聲，如下：

x=dαp+Uγ+N

(4)

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

(5)

(6)

由式(6)可以看出，經(jīng)過(guò)將觀測(cè)信號(hào)投影到非感興趣正交子空間上，背景信號(hào)被消除，原始噪聲也得到了投影壓縮。

3 改進(jìn)MESMA算法

傳統(tǒng)線性光譜混合模型端元矩陣E是固定不變的，然而，實(shí)際上不同像元包含的端元類別與端元個(gè)數(shù)是不同的，如果采用固定的端元解混，反演的豐度是不準(zhǔn)確的，顯然多端元光譜解混模型更加合理。

線性多端元光譜混合模型如下：

(7)

為了解決這個(gè)問(wèn)題，采用均方根誤差變化量 △RMSE來(lái)判斷，如式(7)所示，如果 △RMSE小于設(shè)定的閾值，說(shuō)明多出的這個(gè)端元對(duì)混合像元解混重構(gòu)誤差影響微弱，該端元不構(gòu)成該混合像元，于是選擇個(gè)數(shù)較小的端元組合解混相關(guān)的像元，反之選擇較大個(gè)數(shù)端元組合模型：

(8)

(9)

1)構(gòu)造 3-端元組合為：

(10)

投影值為零的像元U1是由端元E1構(gòu)成的，分離出像元U1，以此類推，分別計(jì)算剩下的像元到端元E2、E3、E4的正交子空間投影得到投影值小于設(shè)定閾值的像元，得到結(jié)果表示為U2、U3、U4，這樣通過(guò)正交子空間投影逐步將像元U1、U2、U3、U4分離。

3)確定構(gòu)成像元U1、U2、U3、U4最優(yōu)端元個(gè)數(shù)。

圖1 解混流程圖Fig. 1 Unmixing flow chart

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

為了驗(yàn)證改進(jìn)算法的有效性，采用真實(shí)數(shù)據(jù)和模擬數(shù)據(jù)進(jìn)行仿真實(shí)驗(yàn)，并將改進(jìn)的算法與另外3種算法進(jìn)行對(duì)比。計(jì)算機(jī)的硬件配置為：處理器為IntelCorei3，主頻為 2.13GHz，內(nèi)存為DDR3 4GB，仿真軟件為MATLAB7.11 版本。

為了客觀地評(píng)價(jià)混合像元分解的精度，采用定性和定量?jī)煞N評(píng)價(jià)準(zhǔn)標(biāo)[10]。定性評(píng)價(jià)采用光譜解混分量圖：它是某類別成分在圖像各像元中所占比例的二維灰度顯示，每個(gè)類別對(duì)應(yīng)一幅解混分量圖。定量分析對(duì)于已知真實(shí)數(shù)據(jù)的精度評(píng)價(jià)采用反演豐度與參考值之間的均方根誤差RMSE為：

(11)

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

從美國(guó)地質(zhì)調(diào)查局(USGS)礦物光譜庫(kù)選取 6 種礦物，分別為明礬石、高嶺石、蒙脫石、綠脫石、白云母和高嶺土，每類礦物選擇 5 條類內(nèi)光譜，用這 30 條光譜合成具有 5 000 個(gè)像元的具有 224 波段的高光譜數(shù)據(jù)，每個(gè)像元由 6 類礦物(每類礦物選擇一條光譜)中的若干種隨機(jī)混合而成，滿足豐度非負(fù)、和為一條件限制，為模擬真實(shí)的高光譜圖像，附加信噪比為 35dB的高斯白噪聲。

表1 解混豐度誤差對(duì)比

表1 是 6 種礦物反演豐度誤差對(duì)比，表中顯示改進(jìn)算法反演豐度誤差均比其他3種算法低，并且前三者算法誤差都低于傳統(tǒng)固定端元解混算法，對(duì)于克服類內(nèi)光譜變化均起到了積極作用，同時(shí)說(shuō)明改進(jìn)算法的有效性。

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

實(shí)驗(yàn)采用了美國(guó)航空可見(jiàn)光/紅外成像光譜儀所測(cè)美國(guó)印第安納州實(shí)驗(yàn)田高光譜數(shù)據(jù)。圖像大小為 144×144，220 波段，將原始的 220 波段中受噪聲影響較大的一些波段去除后選取 200 波段作為仿真研究對(duì)象。從中選擇 3 類地物，干草，喬木，玉米，分別有 489、1294、834 樣本點(diǎn)，平均化后獲得3類代表光譜如圖 2 所示，其空間分布標(biāo)號(hào)圖如圖 3(a) 所示。

圖2 干草、喬木、玉米光譜曲線Fig. 2 Spectral curves of hong-windrowed woods and corn

(a) 真實(shí)地物分布標(biāo)號(hào)圖

(b) OSPMESMA算法解混分量圖

(c) MESMA算法解混分量圖

(d) SZU算法解混分量圖

(e) LSMA算法解混分量圖圖3 3類地物解混分量圖對(duì)比Fig. 3 Unmixng component figures of three kinds of feature

表 2 是4種算法解混豐度均方根誤差對(duì)比，其中 OSPMESMA 和 MESMA 算法誤差是在 11 種端元組合條件下得到的結(jié)果。從中可以看到，OSPMESMA 算法誤差相對(duì)其他3種算法大大降低。在解混分量圖中，灰度越亮代表該類別在所對(duì)應(yīng)的混合像元中的豐度值越大，反之亦然。圖 3 是更直觀的表示，OSPMESMA 算法3類地物得到很好的體現(xiàn)，解混結(jié)果與真實(shí)地物分布十分接近，而其他3種算法解混結(jié)果較模糊，不能有效區(qū)分不同地物。其中 MESMA 和 SZU 算法不論解混誤差上，還是分量圖上，效果都比傳統(tǒng)固定端元解混效果優(yōu)良，它們?cè)诳朔悆?nèi)光譜變化上均起到了積極作用，也說(shuō)明多端元模型更加合理。

表2 解混豐度誤差對(duì)比

表 3 和表 4 是在類內(nèi)變化光譜增多的情況下，OSPMESMA 和 MESMA 算法解混豐度誤差和解混時(shí)間對(duì)比，由于 SZU 和 LSMA 算法解混用到的解混端元是固定的，解混結(jié)果和類內(nèi)變化光譜個(gè)數(shù)沒(méi)有關(guān)系，所以表中只列出了這兩者的結(jié)果?？梢钥吹?，隨著端元組合數(shù)的增大，OSPMESMA 算法解混誤差始終小于 MESMA 算法，時(shí)間也比后者短，精度和效率上都較優(yōu)良。

表3 解混豐度誤差隨端元組合數(shù)變化結(jié)果

Table 3 The change of abundance error with the number of endmember combinations

Endmembercombinations711151923OSPMESMA0.06360.06350.06320.06260.0621MESMA0.14220.13650.14140.13520.1342

表4 解混時(shí)間隨端元組合數(shù)變化結(jié)果

Table 4 The change of time with the number of endmember combinations

Endmembercombinations711151923OSPMESMA0.3640.5480.7080.9151.36MESMA2.385.368.2310.714.9

圖 4 和圖 5 是4種算法隨著端元組合數(shù)增大，解混豐度誤差和時(shí)間變化趨勢(shì)。SZU 和 LSMA 算法豐度誤差和時(shí)間始終保持不變，OSPMESMA 和 MESMA 算法豐度誤差呈下降趨勢(shì)，而解混時(shí)間逐次遞增，其中 OSPMESMA 算法兩者變化都比較緩慢，而 MESMA 算法解混時(shí)間相對(duì)于解混豐度誤差變化劇烈，解混時(shí)間急速上升。

圖4 解混時(shí)間Fig. 4 Unmixing time

圖5 解混豐度誤差Fig. 5 Abundance error of unmxing

圖 6 是4種算法3種地物豐度誤差分布直方圖，直觀地看到，OSPMESMA 算法豐度誤差聚集在 0 值附近，誤差幾乎都落在 -0.2～0.2 之間，而其他3種算法誤差分布散落，落在 -0.2 和 0.2 以外的像元較多，可見(jiàn)改進(jìn)的算法相對(duì)較好。

圖6 豐度誤差直方圖

4.3 解混結(jié)果分析

通過(guò)上述仿真結(jié)果的對(duì)比，可以得到OSPMESMA 算法解混豐度誤差和時(shí)間都優(yōu)于其他3種算法，采用正交子空間投影確定像元所屬M(fèi)-端元組合，避免了不必要端元組合重復(fù)迭代求解，節(jié)約計(jì)算時(shí)間，同時(shí)采用誤差變化量評(píng)價(jià)標(biāo)準(zhǔn)代替最小均方根誤差，克服了過(guò)擬合現(xiàn)象。MESMA 算法需要窮盡所有端元組合，為每個(gè)像元尋找到最優(yōu)端元組合，計(jì)算時(shí)間相當(dāng)長(zhǎng)，通過(guò)上述實(shí)驗(yàn)可以看到，隨著端元組合數(shù)的增大，豐度誤差降低很小，而時(shí)間卻急速上升，由于尋優(yōu)的過(guò)程中，采用最小重構(gòu)均方根誤差為評(píng)價(jià)標(biāo)準(zhǔn)，會(huì)導(dǎo)致過(guò)擬合現(xiàn)象，解混時(shí)用的端元個(gè)數(shù)要大于實(shí)際端元個(gè)數(shù)，重構(gòu)誤差雖然降低了，但是豐度誤差卻增大了。SZU 算法雖然去除了類內(nèi)光譜變化較大的波段，但是仍然采用固定端元解混，不能有效區(qū)分三類地物。LSMA 算法是傳統(tǒng)固定端元解混，效果最差。

5 結(jié)論

本文提出了一種改進(jìn)的多端元高光譜解混算法，避免了隨著地物類別以及類內(nèi)光譜個(gè)數(shù)的增加，端元組合數(shù)的增大導(dǎo)致計(jì)算量的增大，降低了計(jì)算復(fù)雜度。同時(shí)采用誤差變化代替最小重構(gòu)均方根誤差作為評(píng)價(jià)標(biāo)準(zhǔn)，避免了過(guò)擬合現(xiàn)象。采用模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，結(jié)果表明 OSPMESMA 算法解混豐度誤差和時(shí)間都優(yōu)于其他3種算法，證明提出的算法可以克服類間光譜變化，多端元解混模型更加合理。在正交子空間投影確定像元 M-端元?dú)w屬時(shí)，根據(jù)投影值小于某個(gè)很小的閾值來(lái)確定，閾值是在多次試驗(yàn)的基礎(chǔ)上選擇的，如何自適應(yīng)地選擇閾值是下一步工作的重點(diǎn)。

[1]趙春暉, 成寶芝, 楊偉超. 利用約束非負(fù)矩陣分解的高光譜解混算法[J]. 哈爾濱工程大學(xué)學(xué)報(bào), 2012, 33(3): 377-382. ZHAO Chunhui, CHENG Baozhi, YANG Weichao. Algorithm for hyperspectral unmixing using constrained nonnegative matrix factorization[J]. Journal of Harbin Engineering University, 2012, 33(3): 377-382.

[2]王立國(guó), 張晶. 基于線性光譜混合模型的光譜解混改進(jìn)模型[J]. 光電子·激光, 2010, 21(8): 1222-1226. WANG Liguo, ZHANG Jing. An improved spectral unmixing modeling based on linear spectral mixing modeling[J]. Journal of Optoelectronics·Laser, 2010, 21(8): 1222-1226.

[3]王立國(guó), 趙春暉. 高光譜圖像處理技術(shù)[M]. 北京: 國(guó)防工業(yè)出版社, 2013：1-20. WANG Liguo, ZHAO Chunhui. Processing techniques of hyperspectral imagery[M]. Beijing: National Defense Industry Press, 2013：1-20.

[4]ROBERTS D A, GARDNER M, CHURCH R, et al. Mapping chaparral in the Santa Monica Mountains using multiple endmember spectral mixture models[J]. Remote Sensing of Environment, 1998, 65(3): 267-279.

[5]ASNER G P, LOBELL D B. A biogeophysical approach for automated SWIR unmixing of soils and vegetation[J]. Remote Sensing of Environment, 2000, 74(1): 99-112.

[6]SOMERS B, DELALIEUX S, VERSTRAETEN W W, et al. An automated waveband selection technique for optimized hyperspectral mixture analysis[J]. International Journal of Remote Sensing, 2010, 31(20): 5549-5568.

[7]SOMERS B, ASNER G P, TITS L, et al. Endmember variability in spectral mixture analysis: a review[J]. Remote Sensing of Environment, 2011, 115(7): 1603-1616.

[8]宋梅萍, 張甬榮, 安居白, 等. 基于有效端元集的雙線性解混模型[J]. 光譜學(xué)與光譜分析, 2014, 34(1): 196-200. SONG Meiping, ZHANG Yongrong, AN Jubai, et al. Effective endmembers based bilinear unmixing model[J]. Spectroscopy and Spectral Analysis, 2014, 34(1): 196-200.

[9]吳波, 張良培, 李平湘. 非監(jiān)督正交子空間投影的高光譜混合像元自動(dòng)分解[J]. 中國(guó)圖象圖形學(xué)報(bào), 2004, 9(11): 1392-1396. WU Bo, ZHANG Liangpei, LI Pingxiang. Unsupervised orthogonal subspace projection approach to unmix hyperspectral imagery automatically[J]. Journal of Image and Graphics, 2004, 9(11): 1392-1396.

[10]DEMARCHI L, CANTERS F, CHAN J C, et al. Mapping sealed surfaces from CHRIS/Proba data: a multiple endmember unmixing approach[C]// Proceedings of the 2nd Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS). Reykjavik, 2010: 1-4.

An improved multi-endmember hyperspectral unmixing algorithm

ZHAO Chunhui，CUI Shiling，ZHAO Genping

(College of Information and Communication Engineering, Harbin Engineering University, Harbin 150001，China)

The classical multi-endmember spectral mixture analysis model has shortcomings in computation intensity, cumbersome endmember preselection and over-fitting. To overcome these shortcomings, an improved multi-endmember unmixing algorithm is proposed here. Using the characteristics of orthogonal subspace projection that can distinguish signals of interest, it projects pixels onto the orthogonal subspace composed of all of endmembers of the entire surface feature class. Each class selects only one intra-class spectrum. Then it determines the intra-class spectrum of every feature class to which pixels belong according to their projection values. These pixels are isolated in the next iteration in order to reduce computation. Then the optimal number of endmember combinations can be determined according to the reconstruction error variation, which avoids over-fitting. Experiment results show that the inversion abundance error and unmixing time of the improved algorithm are reduced compared to the original algorithm.

hyperspectral unmixing;multi-endmember; unmixing algorithm; orthogonal subspace projection; error variation; intra-class spectral variability

2014-05-17.

時(shí)間：2015-07-27.

國(guó)家自然科學(xué)基金資助項(xiàng)目(61571145,61405041)；中國(guó)博士后科學(xué)基金資助項(xiàng)目(2014M551221)；黑龍江省自然科學(xué)基金重點(diǎn)資助項(xiàng)目(ZD201216)；哈爾濱市優(yōu)秀學(xué)科帶頭人基金資助項(xiàng)目(RC2013XK009003).

趙春暉(1965-), 男, 教授, 博士生導(dǎo)師.

趙春暉，E-mail：zhaochunhui@hrbeu.edu.cn.

10.3969/jheu.201405042

TP751.1

A

1006-7043(2015)09-1281-06

網(wǎng)絡(luò)出版地址：http://www.cnki.net/kcms/detail/23.1390.U.20150727.1254.002.html