唐曉燕，高 昆，劉 瑩，倪國(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ù)來估計(jì)組成混合像元的端元比例(豐度)?；旌舷裨庾V分解模型可分為兩類:線性光譜混合模型和非線性光譜混合模型。線性光譜混合模型假設(shè)一個(gè)光子只看到一種物質(zhì)，像元的光譜是各個(gè)端元的線性組合［3］。雙線性光譜混合模型在線性模型的基礎(chǔ)上做了改進(jìn)，將兩種物質(zhì)之間的散射作為乘積項(xiàng)加入線性模型［4－5］，提高了模型的精度。

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

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

基于貝葉斯推理方法進(jìn)行參數(shù)估計(jì)的基本思路［10］是基于觀測(cè)樣本、參數(shù)的先驗(yàn)分布與似然函數(shù)，對(duì)所研究的對(duì)象或問題建立參數(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沒有其他信息，因此被假設(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)所示。由于沒有該參數(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)概率分布可用下面公式來計(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)過一段時(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)分布來抽樣估計(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è)采樣值的均值來計(jì)算。

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

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

已有的最小二乘算法和基于分層貝葉斯模型的非線性光譜解混算法，均假設(shè)圖像中所有像元內(nèi)包含的端元種類是相同的，即利用所有的端元組成的端元矩陣對(duì)影像中混合像元進(jìn)行光譜分解。然而，混合像元中包含的端元種類一般是不相同的。實(shí)際上，只有在多種地物的交界處的像元才可能包含多數(shù)或者全部的端元種類。而且，當(dāng)端元中包含高度較高的地物類型時(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;否則保留原來的α。

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

本文采用模擬數(shù)據(jù)和圣地亞哥機(jī)場(chǎng)真實(shí)數(shù)據(jù)來對(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)來對(duì)各算法進(jìn)行評(píng)價(jià)，具體計(jì)算公式如式(9)～(11)。這三個(gè)參數(shù)越小代表算法效果越好。

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

在ENVI軟件自帶的光譜庫(kù)中抽取三種純物質(zhì):瀝青、路、樹，根據(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ì)背景下，混合像元中包含的端元種類一般是不相同的，且當(dāng)端元中包含高度較高的地物類型時(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)高曉健，郭寶峰，于平．高光譜空譜一體化圖像分類研究［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．