基于端元優化的非線性高光譜分解算法

唐曉燕，高 昆，劉 瑩，倪國強

(1．北京理工大學光電學院，光電成像技術與系統教育部重點實驗室，北京100081;2．南陽理工學院電子與電氣工程學院，河南南陽473004)

1 引言

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

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

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

2 基于分層貝葉斯模型參數估計的雙線性混合模型

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

2．1 廣義雙線性模型

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

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

2．2 似然估計

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

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

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

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

2．3 參數先驗分布

1)豐度先驗分布:豐度的和為1約束可以表示為一個端元的豐度αk*是其他端元豐度的函數，如式(3)所示。由于沒有該參數的其他信息，αk*的先驗選擇的是Sk*統一分布。

其中，ζ1和 ζ2為兩個超參數。設ζ1=2 ，ζ2=ζ，可得到一個可調的超參數ζ。

4)超參數先驗分布:超參數ζ對分解的精度影響很大，其先驗分布利用Jefffery無信息先驗原則確定為:

2．4 θ的后驗分布

參數向量θ的聯合后驗概率分布可用下面公式來計算:

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

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

2．5 基于MCMC的參數估計

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

為了獲得對參數 θ =(αT，γT，σ2)T的估計，從待估參數的后驗分布來抽樣估計。根據Gibbs采樣，很容易產生 NMC個采樣，Xθ={σ2(t)，α(t)k*，γ(t)}t=1，…，NMC。這些采樣逐漸接近式(7)中所示的聯合后驗分布，形成的馬爾科夫鏈的統計分布為f。結果，這些參數可經驗的采用最后Nr=NMC－Nbi個采樣值的均值來計算。

其中，x為感興趣的參數。

3 基于端元優化的非線性光譜分解算法

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

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

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

3)對誤差較大的點yi，對豐度向量α中的元素由大到小進行排序，確定排序后的豐度向量為α珘=［，…］，該排序結果對應于端元對該混合像元的貢獻大小;

5)對 αs重新計算 REs，如果 REs＜RE，就接受αs;否則保留原來的α。

4 實驗結果

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

4．1 模擬數據實驗

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

表1 模擬圖像的產生模型和參數

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

表2 解混算法的性能比較

表3給出了用實際端元進行解混時兩幅模擬圖像所需的計算時間，時間的單位為秒。實驗采用的計算機硬件環境為Inter Core i3雙核CPU 3．07GHz，內存2GB，軟件環境為 Microsoft Windows 7、MATLAB 2012(b)。根據表3所示，FCLS解混算法所需時間最少，本文提出的OE－GBM算法所需時間與GBM算法解混所需時間相當。

表3 用實際端元進行解混的計算時間(s)

4．2 真實數據實驗

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

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

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

圖2 圣地亞哥機場真實值解混的重建誤差(RE)

5 結論及討論

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

［1］ TONG Qingxi，ZHANG Bing，ZHENG Lanfen．Hyperspectral remote sensing［M］．Beijing:Higher Education Press，2006．(in Chinese)童慶禧，張兵，鄭蘭芬．高光譜遙感——原理、技術與應用［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］．遙感學報，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的混合α穩定分布參數貝葉斯推理［J］．西安理工大學學報，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．