高光譜解混方法研究

嚴 陽，華文深，劉 恂，崔子浩

(陸軍工程大學石家莊校區 電子與光學工程系， 石家莊 050003)

引 言

高光譜圖像通過光譜儀采集，能同時獲得光譜信息和圖像信息，具有圖譜合一的優點[1]，使得其在軍事目標檢測、農作物分類、礦物探測等多個領域都得到廣泛應用。高光譜圖像的光譜分辨率在不斷提高，但是空間分辨率仍舊較低。由于高光譜遙感圖像在采集圖像時，是以像元為單位來獲取地面物體的光譜信息，高光譜圖像中的每一個像元都對應著具有一定面積的地表區域，而區域的大小由光譜儀的空間分辨率決定。因此，當空間分辨率較低時，圖像中將會出現大量混合像元，導致目標的分類精度降低。若一個像元里僅僅包含一種物體，則該像元是純凈像元，包含純凈的光譜信號的像元稱之為端元;當光譜儀空間分辨率較低時，一個像元里含有多種物質混合，則包含混合光譜信號的像元稱作混合像元[2]。高光譜像元解混主要分為線性和非線性解混，線性光譜解混相對簡單易行，非線性解混考慮了物體的二次散射效應，更符合實際光譜采集情況，但情況相對復雜，相關因素眾多，導致解混難度大，仍處于研究初期[3-4]。本文中主要介紹線性光譜解混，經過線性光譜解混,得到混合像元中各個端元在每個像元中所占的比例，該比例稱為豐度，該過程稱為豐度反演。

高光譜線性解混主要包括基于幾何學和基于統計學光譜解混的兩類方法[2]，基于幾何學方法,計算簡單、計算量小，先提取端元，然后求解豐度；而基于統計學方法，通過迭代求解端元和豐度，計算量較大。

基于幾何學光譜解混的典型方法有：像元純凈指數(pixel purity index,PPI)[5]、N-FINDR[6]、頂點成分分析(vertex component analysis,VCA)[7]、自動目標生成(automatic target generation process,ATGP)[8]、正交子空間投影(orthogonal subspace projection,OSP)[9]和連續最大角凸錐(sequential maximum angle convex cone,SMACC)[10]算法等。

基于統計學光譜解混的方法是一種混合信號分離的方法，即把一個矩陣分解成兩個矩陣的乘積，主要有:非負矩陣分解(nonnegative matrix factorization,NMF)[11]、獨立成分分析(independent component analysis,ICA)[12]和迭代誤差分析(iterative error analysis,IEA)[13]等。

近年來，國內外眾多學者對混合像元分解方法進行了深入的研究，同時對現有的方法進行了改進和優化，還有結合圖像的空間信息進行端元提取。2002年,PLAZA提出了采用自動形態學端元提取(automatic morphological endmember extraction,AMEE)[14]的方法，充分利用了高光譜圖像的空間信息和光譜信息，提高了解混精度。

1 解混方法

假設高光譜數據表達為X={x1,x2,…,xN}，其包含N個像元、L個波段，其中xi={y1,y2,…,yL},xi表示第i個像元的光譜矢量,yL為第L個波段的光譜值，i=1,2,…,N，根據線性解混模型，高光譜數據又可以表示為P個端元的光譜矩陣及對應的豐度矩陣的組合[2]：

X=f(E,S)+n

(1)

(2)

式中,f(E,S)表示端元矩陣和豐度矩陣的函數關系式，E=(e1,e2,…,eP)∈RL×P代表包含P個端元的光譜矩陣，eP為第P個端元的光譜；S=(s1,s2,…,sP)T∈RP×N是像元對應的豐度矩陣，sz代表第z個端元的豐度值；n為實際存在的誤差及噪聲項;(2)式即為豐度約束條件。

1.1 基于幾何學的解混方法

基于幾何學的解混方法有一個前提：圖像數據中各物質都存在純凈像元(即端元)。幾何學方法是根據像元在幾何空間的分布特性來提取端元，然后再對圖像進行解混。理想情況下認為高光譜數據集的所有像元的空間分布恰好位于一個凸面單形體中[15]。單形體包含所有的數據點，數據點構成的單形體空間中，單形體的頂點就是端元。因此，端元提取是通過尋找對應單形體的頂點來獲得端元[16]。在2維空間數據中，單形體看作一個三角形，3維空間數據中單形體是一個棱錐，分別如圖1和圖2所示。在多維空間中單形體是個多面體。

圖1 2維空間中單形體示意圖

圖2 3維空間中單形體示意圖

高光譜像元線性解混通常分為端元提取和豐度反演兩個部分。基于幾何學的端元提取方法中常規方法有PPI,VCA,N-FINDR方法等。

1.1.1 像元純凈指數PPI 像元純凈指數PPI[5]是采用投影的方法來尋找端元。首先，使用最小噪聲分離(minimum noise fraction,MNF)[17]降維，然后將所有像元向穿過單形體空間中的隨機方向投影，將被投影到兩端的點作為潛在的端元。統計每個像元在所有投影方向作為端元點的總次數，次數越多，說明屬于端元的概率越大。通過迭代投影的方法找到集合中相對純凈的像元，設定閾值，提取端元。

PPI算法思想簡單易行、計算量較小，但是也存在著較多的缺點。PPI的方法屬于半監督式方法，端元提取過程中需要人為進行選擇，缺乏自動性。同時，由于高光譜維數高，需要先用MNF對數據進行降維，即將高維高光譜數據降為低維數據。有多個投影方向，導致測試耗時量大。同時投影方向隨機導致每次提取的結果有差異，算法的魯棒性差。對此，ZENG[18]引入Givens旋轉矩陣選擇了0°,45°,90°,135°這4個基本方向作為投影方向,大幅減少原算法中的投影方向，而且選取的方向又涵蓋了各個基本方向，改進了PPI算法，減少了計算，使得算法魯棒性更好。

1.1.2 頂點成分分析VCA 頂點成分分析(VCA)[7]改進了投影方式。首先計算高光譜圖像中的端元數目，假定待提取的端元數目為P,通過P次迭代計算就可以找到單形體的頂點，即提取出端元。首先，找出圖像數據中具有最大投影長度的像元，即模值最大的像元點，作為第1個端元。然后找到一個與該端元光譜方向正交的空間向量，作為第2次迭代的投影方向，將數據集往該方向上投影，找出投影中最大投影值對應的像元，將此像元作為第2個端元。然后將這兩個端元組成端元集確定下次的投影方向，找到下一個端元，以此類推，計算出P個端元。

VCA算法計算簡單、計算量小，并且結果穩定，提取結果常常被作為其他算法的初始化端元。

1.1.3 N-FINDR 不同于PPI和VCA算法，N-FINDR算法是通過計算、比較單形體的體積大小，找到構成最大體積對應的像元，作為端元，端元集構成最大單形體體積[6]。首先，對高光譜數據用MNF進行降維，數據維度降成P-1，然后隨機選擇P個初始的端元點作為端元集。計算由端元集組成的單形體體積。體積計算公式為：

(3)

(4)

式中,eP代表端元向量，E為端元集構成的矩陣，V為矩陣的體積，|E|為E的行列式，abs為取絕對值。依次更換一個像元，重新計算體積，反復進行迭代、比較，遍歷所有像元點，直至找到構成最大體積的端元集。

N-FINDR算法計算量較大，實際使用時往往設置迭代次數。因此初始端元的選擇影響著最后的端元提取結果。

對于N-FINDR算法的改進有很多方法。GENG等人[19]提出一種新的體積計算公式，通過改進體積計算方法，使計算擺脫維數限制。假設單形體的P個頂點為e1,e2,…,eP;AP-1=[e2-e1,e3-e1,…,eP-e1]，對行列式進行變換。則單形體體積計算公式為:

(5)

(6)

式中，det表示行列式的值，T表示轉置。由于(AP-1)TAP-1是方陣，所以此公式能應用于任何維度的高光譜數據，則采用N-FINDR算法時不需要先進行降維處理。

采用上述幾何學方法提取端元之后，再通過豐度反演求解出各個端元對應的豐度。目前已經有大量算法來進行豐度反演。如最小二乘法、凸面幾何學分析方法、獨立成分分析法、濾波向量法、投影尋蹤法、端元投影向量法、正交子空間投影法等[20]。其中線性光譜解混中應用最廣泛的方法是基于最小二乘法的反演算法。

1.1.4 最小二乘法線性光譜解混 依據豐度非負以及和為一這兩個約束條件，可以將最小二乘法分為無約束最小二乘法(unconstrained least squares,UCLS)、非負約束最小二乘法(nonnegative constrained least squares,NCLS)、和為一約束最小二乘法((sum-to-one constrained least squares,SCLS)、全約束最小二乘法(fully constrained least squares,FCLS)[21]4種。

在不考慮任何約束時，用無約束最小二乘法(UCLS)求解(1)式，豐度可以表示為：

SUCLS=(ETE)-1ETX

(7)

當加入非負約束(NCLS)時，由于存在不等式約束，模型沒有解析解，利用迭代方法獲得最優解，迭代求解公式為：

(8)

式中，λ為拉格朗日乘子：

(9)

當僅加入和為一約束(SCLS)時，利用拉格朗日乘子法求解，可得豐度向量為：

λ(1-1TX)

(10)

(ETE)-11(1T(ETE)-11)(1-1TX)

(11)

式中,F(S)為拉格朗日函數，1為含有P個1的列向量。

FCLS考慮了所有約束條件，更符合實際情況，本文中采用FCLS進行豐度反演。HEINZ等人[21]對高光譜數據矩陣X和端元矩陣E引入影響因子δ，對兩個矩陣做了變形求解：

(12)

利用這種變形，FCLS的解可以根據UCLS算法求得。

1.1.5 濾波向量法線性光譜解混 像元x可以看作信號光譜T和背景光譜B兩部分組成，表示為：

xi=T+B+n

(13)

式中，n為噪聲。將端元矩陣分為信號光譜和背景光譜兩部分，用背景光譜的正交子空間設計濾波器，使混合像元進行濾波處理，去除背景光譜，剩下信號光譜，求解信號端元的豐度。依次更改端元矩陣中的信號光譜與背景光譜，得到混合像元中每個端元對應的豐度值[22]。

綜上所述，幾何學方法基于幾何理論，首先提取端元，然后再進行豐度反演。但是統計學方法能夠同時得到端元和豐度矩陣。

1.2 基于統計學方法

高光譜線性解混方法中基于統計學的方法，主要包括非負矩陣分解、獨立成分分析和迭代誤差分析等。

1.2.1 非負矩陣分解 非負矩陣分解(NMF)是將一個矩陣分解成兩個非負矩陣乘積的過程。解混時使得端元矩陣和豐度矩陣的重構的圖像數據，相對原圖像數據誤差最小化。LEE等人提出用歐氏距離來表示這一過程，常采用歐氏距離的平方[23]：

(14)

(15)

式中,‖·‖F代表F范數。NMF采用梯度更新迭代公式。通常計算量較大，一般設置最大迭代次數和最大誤差得到端元矩陣和豐度矩陣。在實際應用中難以得到全局最優，經常陷入局部最優。為此通過引入約束條件來緩解局部最優的問題，例如，MIAO等人[24]通過把最小體積約束附加到NMF中，使得最終的端元集盡可能準確。另外還有把平滑性和稀疏性,以及以最小距離作為約束條件，都得到了更好的結果[25-26]。

在NMF的初始端元設置時，可以將基于幾何學方法提取的端元作為初始端元，通過改進初始端元的選取，避免最終結果為局部最優解，同時可以減少迭代計算時間。

1.2.2 獨立成分分析 獨立成分分析(ICA)屬于一種非監督盲源信號分離的方法。通常假設信號源是獨立，且數據是非高斯分布[12]。高光譜數據作為混合信號，將端元光譜或者豐度作為源信號，應用ICA進行盲解混。BAYLISS等人[27]以端元光譜作為源信號，使用ICA進行解混。后來經過其他學者的研究發現，以豐度作為信號源解混效果更好，具有更多的統計信息[28-29]。

ICA方法對數據所做的非高斯分布獨立性假設與實際情況不相吻合，真實數據大體上符合高斯分布，這是ICA 解混方法的主要問題。

基于統計學的方法能夠同時提取出端元矩陣和豐度矩陣，但是求解過程相對復雜，計算量大。

1.2.3 迭代誤差分析 迭代誤差分析(IEA)不需要對數據進行降維，通過計算誤差的大小判斷端元位置。首先將數據均值作為初始向量，然后由最小二乘法解混，求解一個估計豐度，根據已有端元和豐度重構圖像，找到誤差最大的像元作為新的端元；再由更新的端元集迭代，再次解混，直至找到P個端元或誤差達到設定值。

IEA通過不斷更新端元集，對每一個像元迭代求解，能同時求解端元和豐度，解混精度相對較高，但是每次都要計算各像元的模值，總體計算量較大。對此，ZHAO等人提出一種改進的迭代誤差分析方法[13]，求端元集的正交子空間，把所有像元投影到該子空間中，去掉投影值小于閾值的像元，計算剩余像元的均方根誤差。通過投影的方法減少了冗余像元的計算。

1.3 空譜信息結合的方法

高光譜不僅含有光譜信息，同時還包含普通圖像的空間信息。2002年,PLAZA將數學形態學方法應用于多光譜圖像像元解混中，實現了多光譜圖像結合空間信息的端元自動提取(automatic morphological endmember extraction,AMEE)[14], 融合空間信息提高了解混精度。同時，使得AMEE運用于高光譜端元提取成為現實。

AMEE算法采用數學形態學中膨脹和腐蝕兩種運算方法，將空間和光譜信息結合再對高光譜數據進行端元提取,該方法不需要提前進行降維。首先,設置最小和最大空間窗，稱為結構元素。圖像數據在結構元素內進行膨脹和腐蝕基本操作，依次在每個鄰域空間中得到最純光譜像元和混合度最高的像元，用形態學離心指數(morphological eccentricity index,MEI)判斷結構元素內像元純度，其中MEI值由光譜角余弦值計算。依次增大結構元素大小直至預設的最大值，求取所有的MEI值的平均值作為閾值，若像元的MEI值大于該閾值，則該像元屬于端元[14]。此后，又有學者提出了一些空間預處理的方法[30-32]，結合空間信息,得到了比較好的解混效果。

高光譜線性解混方法及其改進優化方法有很多，各種方法也各有特點。目前這些方法中，主要存在的問題是對圖像中異常點的處理，異常點對端元提取的影響十分嚴重，主要的解決思路是根據物體的聚類特性，結合空間信息消除異常點。

2 算法對比

本文中用1997年AVIRIS獲取的美國內達華州Cuprite礦物區真實高光譜數據進行端元提取，然后利用FCLS解混，比較端元矩陣和豐度矩陣的重構圖像與原數據的均方根誤差(root mean square error,RMSE)[16]大小，綜合判斷端元提取方法的特點。對Cuprite數據裁剪得到大小為100×100的高光譜數據，圖3中為第40、第21、第11波段合成的RGB圖像，數據波長范圍為0.37μm～2.51μm，低信噪比波段第98～第128,第148～第170被移除，本實驗中總共采用170個波段。

實驗中采用PPI,N-FINDR,ATGP,SMACC和IEA這5種端元提取方法和FCLS進行豐度反演。5種端元提取及豐度反演之后的RMSE分別如圖3b～圖3f所示，圖3g為RMSE的圖例。均方根誤差圖像中，每個像元的誤差值由圖例所示顏色來表示，數值越大，代表算法的誤差越大。從圖中可以看出，在本實驗中，IEA和N-FINDR算法的解混精度相對其它算法的解混精度更高，這兩種算法對噪聲的抑制效果更佳。表1為5種方法的RMSE數值。通過各個算法總體的RMSE數值也可以看出，IEA和N-FINDR算法的均方根誤差相對更小，而ATGP和SMACC這兩種算法精度最低，因為其通過尋找最亮的像元點，及與該像元點差別最大的點，算法受誤差及噪聲影響較大。根據實驗結果及過程將各方法的特點進行了歸納對比，如表2所示。

圖3 5種端元提取方法的RMSE

a—RGB b—PPI c—N-FINDR d—ATGP e—SMACC f—IEA g—legend

表1 端元提取算法的均方根誤差

表2 常用端元提取算法的特點

從這幾種典型的端元提取方法可以看出，基于幾何學的解混方法復雜度較低，基于統計學盲源信號分離的方法復雜度較高，但解混精度相對更高。大部分算法屬于全自動，只有PPI算法需要手動劃分感興趣區，屬于半自動。同時可以看到，大部分算法中考慮空間信息的算法相對較少。

3 結束語

對高光譜解混方法中的線性解混方法進行了簡要的介紹與總結。總體而言，高光譜線性解混方法是目前研究的重點內容。基于幾何學理論的解混方法相對較多，其思路簡單易行。基于統計學的方法，算法復雜度相對更高。

筆者認為今后應該主要從以下幾個方面對高光譜解混作進一步研究：

(1)將高光譜圖像中的空間信息與光譜信息相結合。高光譜圖像不僅僅包含光譜信息，而且還具有普通圖像的空間位置信息。將空間信息與光譜信息進行融合，能夠減少圖像中異常點的影響。使用張量表達的方法融合高光譜圖像的多種特性，提高解混精度。

(2)充分結合高光譜圖像數據自身獨特的特點，結合稀疏性和空間中相同物質的聚類特性來進行解混。

(3)借鑒線性解混方法，依據非線性解混方法的特點，研究采用更符合實際情況的非線性解混方法來解混。非線性解混方法考慮了二次散射，更加符合實際物體間光譜的作用情況。

(4)改善和優化非線性解混算法的核函數形式，簡化計算。非線性解混算法相對復雜，現有的非線性解混大部分利用高斯徑向基核函數。考慮如何針對高光譜數據特性，改進核函數，改善其性能，提高解混效率。

隨著高光譜遙感技術的發展，圖像中包含的信息將會更加豐富，數據量更大，高光譜解混將會更加重要，在單個像元以及其子像元方面的研究會更加深入。