唐曉燕，高 昆，劉 瑩，倪國(guó)強(qiáng)

(1．北京理工大學(xué)光電學(xué)院，光電成像技術(shù)與系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室，北京100081;2．南陽(yáng)理工學(xué)院電子與電氣工程學(xué)院，河南南陽(yáng)473004)

1 引言

由于地面的復(fù)雜多樣性及傳感器空間分辨率的限制，高光譜圖像上存在大量混合像元，它不僅影響了基于高光譜圖像的地物識(shí)別精度，而且已經(jīng)成為高光譜遙感向定量化方向深入發(fā)展的主要障礙［1－2］。

高光譜圖像可以利用混合像元光譜分解技術(shù)來(lái)估計(jì)組成混合像元的端元比例(豐度)。混合像元光譜分解模型可分為兩類(lèi):線性光譜混合模型和非線性光譜混合模型。線性光譜混合模型假設(shè)一個(gè)光子只看到一種物質(zhì)，像元的光譜是各個(gè)端元的線性組合［3］。雙線性光譜混合模型在線性模型的基礎(chǔ)上做了改進(jìn)，將兩種物質(zhì)之間的散射作為乘積項(xiàng)加入線性模型［4－5］，提高了模型的精度。

對(duì)于非線性混合模型的光譜分解是個(gè)很具挑戰(zhàn)性的問(wèn)題。幾乎所有的基于非線性模型的解混算法都是采用最小二乘估計(jì)［4，6］。近來(lái)一些學(xué)者提出了基于支撐向量回歸［7］和神經(jīng)網(wǎng)絡(luò)［8］的非線性分解方法。Halimi等［4］采用分層貝葉斯算法來(lái)估計(jì)廣義雙線性模型(Generalized Bilinear Model，GBM)的豐度系數(shù)和噪聲變量，取得了較好的結(jié)果，但是該算法是利用所有的端元來(lái)進(jìn)行計(jì)算的。自然界地物分布的復(fù)雜性導(dǎo)致了高光譜圖像通常是非均質(zhì)的，例如背景不單一、地物空間分布不均勻，較高的地物存在的陰影影響等。本文針對(duì)非均質(zhì)背景下，混合像元中包含的端元種類(lèi)一般是不相同的［9］，且當(dāng)端元中包含高度較高的地物類(lèi)型時(shí)，高光譜圖像中不可避免受到陰影的影響，導(dǎo)致光譜分解精度下降，利用端元的可變性，對(duì)基于分層貝葉斯模型非線性光譜分解算法進(jìn)行改進(jìn)，加入陰影端元對(duì)混合像元端元集進(jìn)行優(yōu)化，用優(yōu)化端元子集采用基于分層貝葉斯模型的非線性光譜分解算法來(lái)獲取相應(yīng)端元的豐度。

2 基于分層貝葉斯模型參數(shù)估計(jì)的雙線性混合模型

基于貝葉斯推理方法進(jìn)行參數(shù)估計(jì)的基本思路［10］是基于觀測(cè)樣本、參數(shù)的先驗(yàn)分布與似然函數(shù)，對(duì)所研究的對(duì)象或問(wèn)題建立參數(shù)后驗(yàn)分布的概率模型(Bayesian模型)。

2．1 廣義雙線性模型

廣義雙線性混合模型［4］如式(1)所示:

其中，y∈RL為每個(gè)像元的光譜列向量，L為譜段數(shù)。端元光譜 mk=k=1，2，…，R 。α =是每個(gè)端元在像元點(diǎn)y中所占的豐度分?jǐn)?shù)向量，n= [n1，n2，…，nR]T是附加的白噪聲，通常假設(shè)其為獨(dú)立的，且服從0均值方差為σ2的高斯分布，n ～ Ν(0L，σ2IL)，IL是L×L矩陣。

2．2 似然估計(jì)

式(1)的觀測(cè)模型和噪聲n的高斯參數(shù)滿足:

2)非線性參數(shù)先驗(yàn)分布:由于豐度之間的相互作用通常是小于單個(gè)豐度值，參數(shù)γi，j被假設(shè)是非負(fù)的且小于1。γi，j沒(méi)有其他信息，因此被假設(shè)滿足［0，1］之間的均勻分布［11］。

3)噪聲方差的先驗(yàn)分布:通常假設(shè)σ2的先驗(yàn)服從共軛逆 Gamma分布［12］:

||·||指的是標(biāo)準(zhǔn)l2范數(shù)。y為高光譜圖像中一個(gè)像元的光譜向量，是已知的數(shù)據(jù)，未知參數(shù)向量θ=(αT，γT，σ2)T。下面介紹根據(jù)分層貝葉斯模型估計(jì)GBM 的未知參數(shù)向量 θ =(αT，γT，σ2)T。

2．3 參數(shù)先驗(yàn)分布

1)豐度先驗(yàn)分布:豐度的和為1約束可以表示為一個(gè)端元的豐度αk*是其他端元豐度的函數(shù)，如式(3)所示。由于沒(méi)有該參數(shù)的其他信息，αk*的先驗(yàn)選擇的是Sk*統(tǒng)一分布。

其中，ζ1和 ζ2為兩個(gè)超參數(shù)。設(shè)ζ1=2 ，ζ2=ζ，可得到一個(gè)可調(diào)的超參數(shù)ζ。

4)超參數(shù)先驗(yàn)分布:超參數(shù)ζ對(duì)分解的精度影響很大，其先驗(yàn)分布利用Jefffery無(wú)信息先驗(yàn)原則確定為:

2．4 θ的后驗(yàn)分布

參數(shù)向量θ的聯(lián)合后驗(yàn)概率分布可用下面公式來(lái)計(jì)算:

其中，∝表示“成正比例”，f(θζ)是式(2)定義的似然函數(shù)。假設(shè)所有參數(shù)是先驗(yàn)獨(dú)立的，則f(θζ)

因?yàn)閒(α)需滿足和為1和非負(fù)約束，用最小均方根誤差和最大先驗(yàn)估計(jì)結(jié)合的后驗(yàn)分布并不容易確定。

2．5 基于MCMC的參數(shù)估計(jì)

MCMC的基本思想是構(gòu)建一條馬爾科夫鏈，使其平衡分布為P，然后對(duì)這條馬爾科夫鏈模擬并對(duì)其平衡分布采樣。而GibbS抽樣方法是最簡(jiǎn)單的MCMC方法，適用于條件概率分布容易計(jì)算的情況。，代入式(6)可得到 θy的后驗(yàn)分布:它依據(jù)所有其他變量的當(dāng)前值，對(duì)其中每一個(gè)變量進(jìn)行迭代抽樣。經(jīng)過(guò)一段時(shí)間的迭代后，可以認(rèn)為時(shí)刻x(t+k)的邊緣分布趨向平穩(wěn)，此時(shí)它收斂，而在收斂出現(xiàn)前的Nbi次迭代中，由于初值的影響，各狀態(tài)的邊際分布還不能認(rèn)為是P(x1，x2，…，xl)。因此在估計(jì)期望時(shí)，應(yīng)該把前Nbi次迭代舍去。

為了獲得對(duì)參數(shù) θ =(αT，γT，σ2)T的估計(jì)，從待估參數(shù)的后驗(yàn)分布來(lái)抽樣估計(jì)。根據(jù)Gibbs采樣，很容易產(chǎn)生 NMC個(gè)采樣，Xθ={σ2(t)，α(t)k*，γ(t)}t=1，…，NMC。這些采樣逐漸接近式(7)中所示的聯(lián)合后驗(yàn)分布，形成的馬爾科夫鏈的統(tǒng)計(jì)分布為f。結(jié)果，這些參數(shù)可經(jīng)驗(yàn)的采用最后Nr=NMC－Nbi個(gè)采樣值的均值來(lái)計(jì)算。

其中，x為感興趣的參數(shù)。

3 基于端元優(yōu)化的非線性光譜分解算法

已有的最小二乘算法和基于分層貝葉斯模型的非線性光譜解混算法，均假設(shè)圖像中所有像元內(nèi)包含的端元種類(lèi)是相同的，即利用所有的端元組成的端元矩陣對(duì)影像中混合像元進(jìn)行光譜分解。然而，混合像元中包含的端元種類(lèi)一般是不相同的。實(shí)際上，只有在多種地物的交界處的像元才可能包含多數(shù)或者全部的端元種類(lèi)。而且，當(dāng)端元中包含高度較高的地物類(lèi)型時(shí)，高光譜圖像中不可避免地會(huì)出現(xiàn)陰影。本文利用端元的可變性，加入陰影端元對(duì)混合像元端元集進(jìn)行優(yōu)化，并用對(duì)優(yōu)化的端元子集利用分層貝葉斯模型進(jìn)行豐度估計(jì)。算法的具體步驟如下:

1)對(duì)于圖像 Y= ［y1，y2，…，yN］，根據(jù)初始端元集 M= ［m1，m2，…，mp］，利用分層貝葉斯模型計(jì)算每個(gè)像元的豐度向量α和非線性系數(shù)γ;

2)計(jì)算每個(gè)像元的重建誤差(Reconstruction Error，RE)，對(duì) RE 利用 Ostu 算法［13］進(jìn)行閾值分割，獲得誤差較大的像元點(diǎn)集;

3)對(duì)誤差較大的點(diǎn)yi，對(duì)豐度向量α中的元素由大到小進(jìn)行排序，確定排序后的豐度向量為α珘=［，…］，該排序結(jié)果對(duì)應(yīng)于端元對(duì)該混合像元的貢獻(xiàn)大小;

5)對(duì) αs重新計(jì)算 REs，如果 REs＜RE，就接受αs;否則保留原來(lái)的α。

4 實(shí)驗(yàn)結(jié)果

本文采用模擬數(shù)據(jù)和圣地亞哥機(jī)場(chǎng)真實(shí)數(shù)據(jù)來(lái)對(duì)提出的算法進(jìn)行驗(yàn)證。為了驗(yàn)證算法的有效性，采用全約束的最小二乘算法(Fully Constrained Least Square，F(xiàn)CLS)和基于分層貝葉斯估計(jì)GBM的算法與本文提出優(yōu)化端元的GBM(Optimal Endmembers GBM，OE－GBM)進(jìn)行對(duì)比。試驗(yàn)結(jié)果采用估計(jì)豐度和實(shí)際豐度的豐度均方根誤差［14］(A-bundance RootMean Square Error，ARMSE)，像元光譜估計(jì)值和實(shí)際值的重建誤差和光譜角分布［15］(Spectral Angle Map，SAM)來(lái)對(duì)各算法進(jìn)行評(píng)價(jià)，具體計(jì)算公式如式(9)～(11)。這三個(gè)參數(shù)越小代表算法效果越好。

4．1 模擬數(shù)據(jù)實(shí)驗(yàn)

在ENVI軟件自帶的光譜庫(kù)中抽取三種純物質(zhì):瀝青、路、樹(shù)，根據(jù)AVIRIS高光譜圖像首先進(jìn)行重采樣得到220個(gè)譜段，然后去除水汽吸收和高噪聲譜段，剩余162個(gè)譜段用于生成兩幅模擬數(shù)據(jù)。圖像 I1采用 LMM 模型，α1，α2，α3滿足［0，1］直接的均勻分布，且滿足和為1約束，從中隨機(jī)選取10%的點(diǎn)，用陰影端元代替某個(gè)端元。圖像I2采用GBM 模型，參數(shù) α1，α2，α3滿足［0，1］之間的均勻分布，且滿足和為 1 約束，非線性系數(shù) γ1，2，γ1，3，γ2，3滿足［0，1］之間的均勻分布，從中隨機(jī)選取10%的點(diǎn)，用陰影端元代替某個(gè)端元。圖像I1和I2都加入方差σ2=3．0×10－4的高斯噪聲。模擬圖像的產(chǎn)生模型和參數(shù)如表1所示。

表1 模擬圖像的產(chǎn)生模型和參數(shù)

預(yù)熱的迭代次數(shù)Nbi=300，迭代次數(shù)Nr=700，用式計(jì)算參數(shù)估計(jì)值。對(duì)比FCLS、GBM和本文提出的OE－GBM算法的效果。表2給出了不同算法分解兩幅圖像的重建誤差(RE)和光譜角分布(SAM)和豐度均方根誤差(ARMSE)。對(duì)于圖像I1，OE－GBM算法得到了最小RE、SAM和ARMSE，效果最好。由于圖像I1是采用線性混合，因此FCLS的解混效果優(yōu)于GBM。對(duì)于圖像I2，OE－GBM算法的RE、SAM和ARMSE也優(yōu)于其它兩種FCLS和GBM算法。每幅圖像比較合適的算法是采用相應(yīng)模型解混的算法。但是，也可看出，采用本文提出OE－GBM算法對(duì)于各種模型產(chǎn)生的圖像，均有較好的解混效果。

表2 解混算法的性能比較

表3給出了用實(shí)際端元進(jìn)行解混時(shí)兩幅模擬圖像所需的計(jì)算時(shí)間，時(shí)間的單位為秒。實(shí)驗(yàn)采用的計(jì)算機(jī)硬件環(huán)境為Inter Core i3雙核CPU 3．07GHz，內(nèi)存2GB，軟件環(huán)境為 Microsoft Windows 7、MATLAB 2012(b)。根據(jù)表3所示，F(xiàn)CLS解混算法所需時(shí)間最少，本文提出的OE－GBM算法所需時(shí)間與GBM算法解混所需時(shí)間相當(dāng)。

表3 用實(shí)際端元進(jìn)行解混的計(jì)算時(shí)間(s)

4．2 真實(shí)數(shù)據(jù)實(shí)驗(yàn)

真實(shí)數(shù)據(jù)實(shí)驗(yàn)采用美國(guó)圣地亞哥機(jī)場(chǎng)的AVIRIS數(shù)據(jù)。波長(zhǎng)為0．389～2．467μm，除去低信噪比和水吸收譜段(譜段1 －2，104－113，148－167，221－224)，后剩余188個(gè)譜段數(shù)據(jù)。為了減少計(jì)算的時(shí)間復(fù)雜度，從原始圖像中截取大小為50×50的子圖(如圖1所示)。從圖1可看出，端元數(shù)目為4，分別是:飛機(jī)、混凝土1、混凝土2、硬土。由于相同地物光譜的可變性，4種端元的參考光譜本文從圖中人工獲得，例如飛機(jī)的參考光譜，從飛機(jī)機(jī)身的中心位置抽取4個(gè)像元，求其光譜的均值作為飛機(jī)的參考光譜。

圖1 圣地亞哥機(jī)場(chǎng)的偽彩色圖像(R:譜段90;G:譜段50;B:譜段27)

用OE－GBM、FCLS、GBM三種算法進(jìn)行解混得到的RE圖像如圖2所示。從圖2中可看出，GBM解混算法的RE優(yōu)于FCLS，但是由于光照的不均勻和陰影的影響，在飛機(jī)區(qū)域的誤差較大。OE－GBM算法考慮了陰影端元，且對(duì)參與分解的端元集做了優(yōu)化，因此RE小于FCLS和GBM算法。

圖2 圣地亞哥機(jī)場(chǎng)真實(shí)值解混的重建誤差(RE)

5 結(jié)論及討論

本文針對(duì)非均質(zhì)背景下，混合像元中包含的端元種類(lèi)一般是不相同的，且當(dāng)端元中包含高度較高的地物類(lèi)型時(shí)，高光譜圖像中不可避免會(huì)受到陰影的影響，導(dǎo)致光譜分解精度下降，利用端元的可變性對(duì)基于分層貝葉斯模型的雙線性光譜解混算法進(jìn)行改進(jìn)。加入陰影端元對(duì)混合像元的端元集進(jìn)行優(yōu)化，并用優(yōu)化端元子集采用基于分層貝葉斯模型的雙線性光譜解混算法進(jìn)行光譜解混，實(shí)驗(yàn)結(jié)果表明，提出的OE－GBM算法能很好的解決高光譜圖像中存在的陰影，光譜分解效果優(yōu)于 FCLS和 GBM算法。

［1］ TONG Qingxi，ZHANG Bing，ZHENG Lanfen．Hyperspectral remote sensing［M］．Beijing:Higher Education Press，2006．(in Chinese)童慶禧，張兵，鄭蘭芬．高光譜遙感——原理、技術(shù)與應(yīng)用［M］．北京:高等教育出版社，2006．

［2］ GAO Xiaojian，GUO Baofeng，YU Ping．Classification of hyperspectral remote sensing image based on spatialspectral integration［J］．Laser ＆ Infrared，2013，43(11):1296 －1300．(in Chinese)高曉健，郭寶峰，于平．高光譜空譜一體化圖像分類(lèi)研究［J］．激光與紅外，2013，43(11):1296 －1300．

［3］ Dobigeon N，Tourneret J，Richard C，et al．Nonlinear unmixing of hyperspectral images:Models and algorithms［J］．IEEE Signal Processing Magazine，2014，31(1):82－94．

［4］ FANWenyi，HU Baoxin，Miller J，etal．Comparative study between a new nonlinearmodel and common linearmodel for analysing laboratory simulated－forest hyperspectral data［J］．International Journal of Remote Sensing，2009，30(11):2951－2962．

［5］ Halimi A，Altmann Y，Dobigeon N，et al．Nonlinear unmixing of hyperspectral images using a generalized bilinearmodel［J］．IEEE Transactions on Geoscience and Remote Sensing，2011，49(11):4153 －4162．

［6］ Combe JP，Launeau P，Carrère V，et al．Mapping microphytobenthos biomass by non－linear inversion of visible－ infrared hyperspectral images［J］．Remote Sensing of Environment，2005，98(4):371 －387．

［7］ WU Bo，ZHANG Liangpei，LIPingxiang．Unmixing hyperspectral imagery based on support vector nonlinear approximating regression［J］．Journal of Remote Sensing，2006，10(3):312 －318．(in Chinese)吳波，張良培，李平湘．基于支撐向量回歸的高光譜混合像元非線性分解［J］．遙感學(xué)報(bào)，2006，10(3):312 －318．

［8］ Liu W G，Wu E Y．Comparison of non －linear mixture models:sub － pixel classification［J］．Remote Sensing of Environment，2005，94(2):145 －154．

［9］ Rogge D M，Rivard B，Zhang J，et al．Iterative spectral unmixing for optimizing per － pixel endmember sets［J］．IEEE Transactions on Geoscience and Remote Sensing，2006，44(12):3725 －3736．

［10］ CHEN Yajun，LIU Ding，LIANG Junli．Bayesian inference on parameters formixtures ofα－stable distributions based on markov chainmonte carlo［J］．Journal of Xi＇an University of Technology，2012，28(4):385 －391．(in Chinese)陳亞軍，劉丁，梁軍利．基于MCMC的混合α穩(wěn)定分布參數(shù)貝葉斯推理［J］．西安理工大學(xué)學(xué)報(bào)，2012，28(4):385－391．

［11］ Somers B，Cools K，Delalieux S，et al．Nonlinear hyperspectralmixture analysis for tree cover estimates in orchards［J］．Remote Sensing of Environment，2009，113(6):1183－1193．

［12］ Dobigeon N，Tourneret JY，Chang C I．Semi－ supervised linear spectral unmixing using a hierarchical Bayesian model for hyperspectral imagery［J］．IEEE Transactions on Signal Processing，2008，56(7):2684 －2695．

［13］ Otsu N．A threshold selectionmethod from gray－levelhistograms［J］．Automatica，1975，11(285 －296):23 －27．

［14］ Plaza J，Plaza A，Pérez R，et al．Joint linear/nonlinear spectral unmixing of hyperspectral image data［C］//IEEEGeoscience and Remote Sensing Symposium，2007:4037－4040．

［15］ Keshava N，Mustard JF．Spectral unmixing［J］．IEEE Signal Processing Magazine，2002，19(1):44 －57．