高光譜解混方法研究

嚴(yán) 陽，華文深，劉 恂，崔子浩

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

引 言

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

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

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

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

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

1 解混方法

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

X=f(E,S)+n

(1)

(2)

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

1.1 基于幾何學(xué)的解混方法

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

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

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

高光譜像元線性解混通常分為端元提取和豐度反演兩個部分?；趲缀螌W(xué)的端元提取方法中常規(guī)方法有PPI,VCA,N-FINDR方法等。

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

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

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

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

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

(3)

(4)

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

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

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

(5)

(6)

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

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

1.1.4 最小二乘法線性光譜解混 依據(jù)豐度非負(fù)以及和為一這兩個約束條件，可以將最小二乘法分為無約束最小二乘法(unconstrained least squares,UCLS)、非負(fù)約束最小二乘法(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)

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

(8)

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

(9)

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

λ(1-1TX)

(10)

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

(11)

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

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

(12)

利用這種變形，F(xiàn)CLS的解可以根據(jù)UCLS算法求得。

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

xi=T+B+n

(13)

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

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

1.2 基于統(tǒng)計學(xué)方法

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

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

(14)

(15)

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

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

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

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

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

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

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

1.3 空譜信息結(jié)合的方法

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

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

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

2 算法對比

本文中用1997年AVIRIS獲取的美國內(nèi)達華州Cuprite礦物區(qū)真實高光譜數(shù)據(jù)進行端元提取，然后利用FCLS解混，比較端元矩陣和豐度矩陣的重構(gòu)圖像與原數(shù)據(jù)的均方根誤差(root mean square error,RMSE)[16]大小，綜合判斷端元提取方法的特點。對Cuprite數(shù)據(jù)裁剪得到大小為100×100的高光譜數(shù)據(jù)，圖3中為第40、第21、第11波段合成的RGB圖像，數(shù)據(jù)波長范圍為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的圖例。均方根誤差圖像中，每個像元的誤差值由圖例所示顏色來表示，數(shù)值越大，代表算法的誤差越大。從圖中可以看出，在本實驗中，IEA和N-FINDR算法的解混精度相對其它算法的解混精度更高，這兩種算法對噪聲的抑制效果更佳。表1為5種方法的RMSE數(shù)值。通過各個算法總體的RMSE數(shù)值也可以看出，IEA和N-FINDR算法的均方根誤差相對更小，而ATGP和SMACC這兩種算法精度最低，因為其通過尋找最亮的像元點，及與該像元點差別最大的點，算法受誤差及噪聲影響較大。根據(jù)實驗結(jié)果及過程將各方法的特點進行了歸納對比，如表2所示。

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

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

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

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

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

3 結(jié)束語

對高光譜解混方法中的線性解混方法進行了簡要的介紹與總結(jié)?？傮w而言，高光譜線性解混方法是目前研究的重點內(nèi)容?；趲缀螌W(xué)理論的解混方法相對較多，其思路簡單易行?；诮y(tǒng)計學(xué)的方法，算法復(fù)雜度相對更高。

筆者認(rèn)為今后應(yīng)該主要從以下幾個方面對高光譜解混作進一步研究：

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

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

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

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

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