趙春暉,崔曉辰

(哈爾濱工程大學 信息與通信工程學院,黑龍江 哈爾濱 150001)

基于全局和局部流形結(jié)構(gòu)的高光譜圖像特征提取算法

趙春暉,崔曉辰

(哈爾濱工程大學 信息與通信工程學院,黑龍江 哈爾濱 150001)

針對目前高光譜圖像基于流形學習的無監(jiān)督特征提取算法中只能夠單獨描述高維數(shù)據(jù)空間局部或者全局的幾何結(jié)構(gòu),并且沒有一種算法能夠同時保持高維數(shù)據(jù)全局和局部的幾何結(jié)構(gòu)的問題,提出了一種基于全局和局部流形結(jié)構(gòu)的無監(jiān)督特征提取算法(GLMS) 對高光譜圖像進行特征提取.算法基于流形學習基本理論,需要建立兩種保持流形結(jié)構(gòu)的近鄰圖,分別用來描述數(shù)據(jù)的全局和局部的流形結(jié)構(gòu),通過求解廣義特征值問題獲得重構(gòu)權(quán)值矩陣進而得到低維嵌入空間的最優(yōu)投影,以達到降維的目的.在AVIRIS高光譜圖像以及Indian Pine和Salina數(shù)據(jù)集上進行仿真對比實驗,結(jié)果表明,提出的算法在分類精度和計算效率上有較好的提高.

高光譜圖像; 無監(jiān)督特征提取; 全局和局部流形結(jié)構(gòu); 流形學習

高光譜圖像所具有的大量光譜維為地物信息提取提供了極其豐富的信息[1],有助于更加精細的地物分類和目標識別,然而,光譜維的擴展也會導致數(shù)據(jù)的冗余和處理數(shù)據(jù)復雜度的提高.因此,數(shù)據(jù)降維在高光譜圖像處理中起到了十分重要的作用,并且成為高光譜圖像分類、混合像元分解、異常目標檢測和地物反演等常用到的預處理技術(shù).高光譜圖像降維通常包括兩種方式:波段選擇和特征提取.與特征選擇不同的是,特征提取需要對各個光譜波段重新組合和優(yōu)化,通過數(shù)學變換對光譜波段進行壓縮,因此,從原始數(shù)據(jù)集中重建數(shù)據(jù)往往比直接選擇數(shù)據(jù)應(yīng)用到后續(xù)圖像處理會得到更好的效果[2].

近年來,流形學習[3](Manifold Learning)越來越受到學者們的廣泛關(guān)注,逐漸成為廣泛使用的非線性特征提取算法,該理論的基本假設(shè)是:高維空間的數(shù)據(jù)并不是均勻分布的[3],而是嵌入在低維流形上的.Bachmann等首次利用等距特征映射算法(Isomap)[3-4]學習高光譜圖像的低維嵌入流形曲面,Isomap算法是一種典型的全局保持流形結(jié)構(gòu)的無監(jiān)督特征提取算法.局部線性嵌入算法(LLE)[5]考慮到高維數(shù)據(jù)空間不存在全局的線性結(jié)構(gòu),但是在樣本點的鄰域范圍內(nèi),可以近似地看成符合局部線性歐氏空間.隨后提出的拉普拉斯特征映射算法[6](LE)和局部切空間排列算法[7](LTSA)也都是基于局部流形結(jié)構(gòu)這一假設(shè)條件.然而,LLE、LE、LTSA這些局部算法并不能給出確定的映射函數(shù),以至于不能將所提的算法應(yīng)用到新的數(shù)據(jù)集上,為此,He等對LE算法進行了改進,相繼提出了鄰域保持嵌入算法[8](NPE)和局部保持投影算法[9](LLP),解決了樣本外問題,即能夠求解出映射函數(shù).然而上述算法都只能夠?qū)W習到單一的流形結(jié)構(gòu),為此,王立志提出了一種線性局部與全局保持的多流行學習算法[10],旨在學習每類數(shù)據(jù)的幾何結(jié)構(gòu),該方法雖然解決了單一流行結(jié)構(gòu)問題,但需要利用到樣本的類別信息,為此本文提出了一種基于全局和局部流形結(jié)構(gòu)的無監(jiān)督特征提取算法(Global and Local Manifold Structure Feature Extraction Algorithm,GLMS),并且通過引入?yún)?shù)使得求解的重構(gòu)權(quán)值矩陣達到最優(yōu).無監(jiān)督特征提取算法通過構(gòu)造圖嵌入的模型來同時保持不同類別數(shù)據(jù)的幾何結(jié)構(gòu),包括全局和局部幾何結(jié)構(gòu),求解兩種近鄰圖的重構(gòu)權(quán)值矩陣,最后通過求解廣義特征值問題得到最佳投影.

1 基于全局和局部流形結(jié)構(gòu)的特征提取算法

1.1 問題描述

給定一組高光譜數(shù)據(jù)樣本點,首先要建立一個關(guān)于樣本數(shù)據(jù)的近鄰圖.由于需要考慮高維數(shù)據(jù)的全局和局部的幾何結(jié)構(gòu),因此需要建立兩種保持流形結(jié)構(gòu)的重建矩陣.其中一個用來描述局部的流形結(jié)構(gòu),利用近鄰圖中相互連接的數(shù)據(jù)樣本點的線性重構(gòu)權(quán)值矩陣;與此同時,另外一個描述全局的流形結(jié)構(gòu),利用最短路徑算法求解成對的數(shù)據(jù)樣本點,因為這些樣本點在近鄰圖上不是相互連接的.最后通過求解廣義特征值問題獲得最優(yōu)的嵌入投影函數(shù).這樣,不同的局部和全局流形結(jié)構(gòu)才能同時保持于低維的嵌入空間中.

令X={x1,x2,…,xN}(xi)∈(RD)表示高光譜圖像數(shù)據(jù)樣本集,Z={z1,z2,…,zN}(zi)∈(Rd)是其低維表示,其中,d

(1)

式中,投影矩陣W=[w1,w2,…,wd]∈RD×d.

1.2 保持局部幾何結(jié)構(gòu)模型

為了保持局部流形結(jié)構(gòu)[11],首先要建立K近鄰圖G,其中,從原始數(shù)據(jù)集中選擇m個樣本點xi(i=1,2,…,m),若樣本點i在樣本點j的鄰域范圍內(nèi),則將這兩個點連起來構(gòu)成連接.從某種意義來說,樣本點會分布在非線性的嵌入流形結(jié)構(gòu)上,但是可以假設(shè)每一個局部的近鄰域都是線性的,類似于局部線性嵌入算法中,可以將這些局部的幾何結(jié)構(gòu)通過線性的相關(guān)矩陣來重建每一個來自其近鄰點的樣本重構(gòu).重構(gòu)誤差函數(shù)如下式所示:

(2)

式中,NK(xi)表示樣本點的K近鄰.需要說明的是,式(2)的約束條件是:如果樣本點j不在樣本點i的K近鄰域內(nèi),那么可以通過求解最小化問題來得到變換矩陣.在低維嵌入空間中,為了保持局部的流形結(jié)構(gòu),需要對下式求取最小化:

trace(YMYT)=trace(ATXMXTA).

(3)

式中,Y=(y1,y2,…,ym),X=(x1,x2,…,xm),M=(I-W)T(I-W),I=diag(1,…,1).可以容易地知道,矩陣M是對稱且半正定矩陣.考慮到一維的情況,最優(yōu)的投影矩陣可以通過求解最小化問題得到解決,如下式所示:

(4)

式中,a是n×1的向量.

1.3 保持全局幾何結(jié)構(gòu)模型

在保持局部流形結(jié)構(gòu)[7]的同時,還要同時考慮高維數(shù)據(jù)的全局幾何結(jié)構(gòu).首先建立一個K近鄰圖G的互補圖Gc.在這個互補圖Gc中,如果樣本點i和j在近鄰圖G中不是互相連接的,那么就需要將這兩個樣本點相互連接起來.這樣,G+Gc就構(gòu)成了樣本點的完整近鄰圖.因此可以確定原始數(shù)據(jù)空間中的全局幾何結(jié)構(gòu),并且通過尋求一種嵌入映射f,可以保持測地線距離在低維嵌入的歐氏空間中,條件是兩個樣本點在互補圖Gc中是互相連接的,如下式所示:

(5)

式中,dM(xi,xj)表示流形結(jié)構(gòu)上的兩個樣本點間的測地線距離,d(f(xi),f(xj))表示在低維嵌入空間中的歐氏距離.在表達真實地物的高光譜圖像數(shù)據(jù)集中,基本的流形結(jié)構(gòu)通常是未知的,因此通過求解樣本點間的測地線距離的方法也是未知的.所以只能夠通過求解近鄰圖G上兩個樣本點間的最短路徑dG(xi,xj)來近似地估計出樣本點間的測地線距離dM(xi,xj).計算最短路徑問題的過程如下:初始化dG(xi,xj)=d(xi,xj),如果樣本點i,j是相互連接的,那么dG(xi,xj)=∞;然后對于數(shù)據(jù)集中的每一個樣本點重復上面的操作,通過求解如下最小化問題來代替求解最短路徑問題dG(xi,xj):

(6)

最終的矩陣DG={dG(xi,xj)}包含在近鄰圖G的所有樣本點的最短路徑,如下式所示:

(7)

1.4 算法描述

(8)

考慮線性的嵌入映射:f(x)=aTx.令yi=f(xi),Y=aTX,則

(9)

所以,最佳的投影矩陣可以通過求解如下最優(yōu)化問題得到解決:

(10)

式(10)的推導過程如下:

trace(aTXXTaaTXXTa).

(11)

(12)

經(jīng)過上述推導,最優(yōu)的投影向量可以表示為

aTXXTaaTXXTa)-λ2aTXMXTa]}.

(13)

(14)

向量ai可以將式(14)最大化,并且可以轉(zhuǎn)換為求解特征值和特征向量問題,如下式所示:

(15)

在實驗中,令λ1=常量,即能夠求解出在不受全局結(jié)構(gòu)影響下的局部最優(yōu)重構(gòu)權(quán)值矩陣;令λ2=常量,即能夠求解出鄰域內(nèi)每一個樣本點都用其鄰域內(nèi)的點線性表示時的全局最優(yōu)解.

令A=(a1,a2,…,ad),線性映射描述如下式所示:

(16)

式中,y為高光譜圖像數(shù)據(jù)的低維表示.

1.5 評價指標

總體的分類精度(Overall Classification Accuracy,OCA)是一個比值,分母代表所有原始的高光譜數(shù)據(jù)的樣本點的數(shù)量,分子代表經(jīng)過分類方法正確分類的數(shù)據(jù)個數(shù),如下式所示:

(17)

2 算法步驟及流程圖

2.1 算法步驟描述

輸入: 高光譜圖像數(shù)據(jù)樣本集X={x1,x2,…,xN}∈RD.

步驟1: 對輸入數(shù)據(jù)進行初始化.

步驟2: 根據(jù)式(4)求解出保持局部幾何結(jié)構(gòu)的重構(gòu)權(quán)值矩陣.

步驟3: 根據(jù)式(7)求解出保持全局幾何結(jié)構(gòu)的重構(gòu)權(quán)值矩陣.

步驟4: 根據(jù)式(14)求解出同時保持高維數(shù)據(jù)樣本的全局和局部幾何結(jié)構(gòu)的投影向量.

步驟5: 計算出重構(gòu)權(quán)值矩陣對應(yīng)的最大的d個特征值對應(yīng)的特征向量,構(gòu)成投影向量A.

步驟6: 求解高維數(shù)據(jù)的低維嵌入坐標,x→y=ATx.

輸出:特征提取算法處理后的低維數(shù)據(jù)Z={z1,z2,…,zN}∈Rd.

2.2 算法流程圖

圖1是求解保持局部幾何結(jié)構(gòu)的重構(gòu)權(quán)值矩陣的算法流程圖.

3 實驗結(jié)果及分析

3.1 實驗數(shù)據(jù)

實驗所用計算機搭載64-bits,主頻3.5 GHz Intel i7 4770 K CPU,內(nèi)存16 G,三級緩存8 MB.軟件環(huán)境:MATLAB 2011b.本文實驗采用AVIRIS高光譜遙感圖像.Indian Pine數(shù)據(jù)集取自1992年6月拍攝的美國印第安納州西北部印第安遙感試驗區(qū)的一部分,圖像大小為145×145像素,共包含220個波段,像素深度為16 bits,波長范圍為400～2 500 nm,包含16類確定類別的地物.AVIRIS Salina高光譜數(shù)據(jù)由AVIRIS成像光譜儀在加利福尼亞地區(qū)收集獲取,擁有工作波段224個,含有512×217個像素,空間分辨率為4 m,包含16種真實地物類別.圖2a為Indian Pine高光譜數(shù)據(jù)集的灰度圖,圖2b為Indian Pine高光譜數(shù)據(jù)集的類別標記圖.圖3a為Salina高光譜數(shù)據(jù)集的灰度圖,圖3b為Salina高光譜數(shù)據(jù)集的類別標記圖.

圖1 求解保持局部幾何結(jié)構(gòu)的重構(gòu)權(quán)值矩陣的算法流程圖Fig.1 The flow chart of reconstruction weight matrix of the local geometric structure

圖2 Indian Pine高光譜數(shù)據(jù)集的灰度圖及類別標記圖Fig.2 Gray scale map and labeled graph ofhyperspectral data set of Indian Pine

實驗的數(shù)據(jù)樣本X為無標簽數(shù)據(jù)樣本部分隨機地采樣于原始的高光譜數(shù)據(jù)集,其中,70%隨機地從原始數(shù)據(jù)集中選取作為訓練樣本,剩余的30%作為測試樣本.為了調(diào)查出訓練樣本大小對分類性能的影響,初始的訓練樣本個數(shù)的選取也應(yīng)不同.實驗中,分類算法選取支持向量機SVM算法,使用RBF Kernels(徑向基核函數(shù)).每次實驗做10次,取平均值,并且結(jié)果利用OCA和Kappa系數(shù)來評價特征提取算法的分類性能.選擇高斯核函數(shù),其寬為0.02,懲罰是10.基于流形學習的特征提取算法的近鄰參數(shù)選取不是固定的,在實驗中,對式(15)的求解通常是令其中一個參數(shù)常量化,即實驗中令λ2=1,保證在不受局部特性干擾的情況下求解出全局最優(yōu)解;并且將本文提出的算法和經(jīng)典的算法進行了比較.本次實驗分別使用PCA,Isomap,LLE,LE,NPE,LPP,LLTSA,SPP,GLMS等算法對性能進行驗證.

圖3 Salina高光譜數(shù)據(jù)集的灰度圖及類別標記圖Fig.3 Gray scale map and labeled graph ofhyperspectral data set of Salina

3.2 實驗結(jié)果

表1為不同特征提取算法的最高總體分類精度顯示.其中,黑體顯示的為每一種數(shù)據(jù)集中不同算法中最高的總體分類精度,括號里表達的是最高總體分類精度對應(yīng)的降維后的最優(yōu)波段數(shù).

表1 不同特征提取算法的最高總體分類精度Table1 The highest overall classification accuracy for different feature extraction

由表1可以總結(jié)出如下結(jié)論.

(1) 對于高光譜數(shù)據(jù)樣本集來說,不同的特征提取算法可以提高數(shù)據(jù)集的分類性能,大多數(shù)的信息都將能保留下來以便為后續(xù)的分類服務(wù).對于Indian Pine數(shù)據(jù)集,經(jīng)過PCA特征提取后的低維子空間經(jīng)過SVM后的總體分類精度并沒有原始數(shù)據(jù)集的精度高,因此可以推論:由于PCA是線性特征提取算法,其不能提取出原始數(shù)據(jù)集中由于大氣、水分等產(chǎn)生的非線性折射產(chǎn)生的數(shù)據(jù)間的相關(guān)性和非線性因素.

(2) 不同的算法在應(yīng)用到不同的數(shù)據(jù)集時,得出的總體分類精度也是有差別的.當標簽樣本的數(shù)量有限時,半監(jiān)督算法相對于其他不考慮正則化參數(shù)的算法來說得到的,總體分類精度就要減少許多.

(3) 與無監(jiān)督算法SPP比較,充分利用了高光譜數(shù)據(jù)少量先驗類別信息的算法,獲得了更高的總體分類精度.相比于其他基于流形學習的降維算法,由于不用考慮參數(shù)的影響,曲線平滑上升.而對參數(shù)敏感的流形算法,在一些特定的維數(shù)時,總體分類精度會不穩(wěn)定.

(4) 本文提出的算法并不是在所有的波段都表現(xiàn)良好,例如在維數(shù)為7～11維時,總體分類精度小于LLTSA算法.在維數(shù)大于14時,總體分類精度是所有算法中最高的,原因是考慮到了數(shù)據(jù)集局部的幾何關(guān)系,效果好于全局算法.而且由于算法不需要設(shè)置額外的參數(shù),因此在總體分類精度這項指標上比SPP算法要稍微弱一點.

圖4為Indian Pine高光譜數(shù)據(jù)集不同算法對應(yīng)的總體分類精度.由圖4可以看出,各算法降維到7波段后,總體分類精度趨于穩(wěn)定并且達到了最高點,由此得出:降維后的波段數(shù)過小,將損失高光譜數(shù)據(jù)樣本的原始信息,然而低維子空間的維數(shù)如果選取得過大,又會使得運算復雜度過高,失去了特征提取的意義,況且得到的低維數(shù)據(jù)與原始數(shù)據(jù)相比并不能體現(xiàn)出降維的優(yōu)勢.由于考慮了全局和局部的幾何結(jié)構(gòu),因此本文提出的算法在分類精度這一指標上要優(yōu)于其他保持單一流形結(jié)構(gòu)的特征提取算法.

圖4 Indian Pine高光譜數(shù)據(jù)集不同算法對應(yīng)的總體分類精度Fig.4 The overall classification accuracy ofdifferent algorithm graph of hyperspectraldata set of Indian Pine

圖5為Salina高光譜數(shù)據(jù)集不同算法對應(yīng)的分類圖.由圖5能夠直觀地看出各種數(shù)據(jù)源下的通過SVM算法分類后的效果圖,能夠得出結(jié)論:經(jīng)過GMLS算法降維后的低維數(shù)據(jù)應(yīng)用SVM支持向量機分類后得到的分類效果在幾種算法中最優(yōu),幾乎能夠識別出全部的不同的地物,而且可以看出其他算法都或多或少地將某些地物錯分.

圖5 Salina高光譜數(shù)據(jù)集不同算法對應(yīng)的分類圖Fig.5 The classification map of different algorithmgraph of hyperspectral data set of Salina

表2所示為Indian Pine數(shù)據(jù)集不同特征提取算法的運算時間.每種算法在相同的條件下,分別運算10次取均值.本文算法的計算復雜度主要包括近鄰圖搜索和廣義特征值求解,與樣本個數(shù)N和維數(shù)D有關(guān),而GLMS算法的計算復雜度主要包括尋找近鄰圖,計算復雜度為o(DN2).當N?D時,本文提出的算法雖然比以往傳統(tǒng)的算法耗費時間長,但是考慮到GLMS算法同時考慮了高維數(shù)據(jù)空間的不同幾何結(jié)構(gòu),而且本文提出的算法在后續(xù)的數(shù)據(jù)分類中顯示出較高的性能,計算時間也只是比SPP算法多了14.79 s,相比于優(yōu)點,這些時間幾乎可以忽略.

表2 Indian Pine數(shù)據(jù)集不同特征 提取算法的運算時間Table2 The operation time of different algorithm graph of hyperspectral data set of Indian Pine

4 結(jié) 語

考慮到半監(jiān)督算法需要事先知道原始數(shù)據(jù)集的先驗知識,而無監(jiān)督算法則不需要,然而,目前基于流形學習的無監(jiān)督算法只能夠單獨地描述局部或者全局的幾何結(jié)構(gòu),通過構(gòu)建兩種圖來實現(xiàn)不同流形結(jié)構(gòu)的統(tǒng)一.然而正是由于構(gòu)建近鄰圖要搜索全局結(jié)構(gòu)和局部結(jié)構(gòu)的圖構(gòu)造策略,造成了在算法復雜度上的提高,因此,算法的改進還要圍繞著如何在減少搜索近鄰圖的時間復雜度上來提高.

[1] Plaza A,Benediktsson J A,Boardman J W,et al.Recent Advances in Techniques for Hyperspectral Image Processing[J].Remote Sensing of Environment,2009,113(1):110-122.

[2] Green R O,Eastwood M L,Sarture C M,et al.Imaging Spectroscopy and the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS)[J].Remote Sensing of Environment,1998,65(3):227-248.

[3] Bachmann C M,Ainsworth T L.Bathymetric Retrieval from Manifold Coordinate Representations of Hyperspectral Imagery[C]∥IEEE International Conference on Geoscience & Remote Sensing Symposium.IEEE,2007:1548-1551.

[4] Bachmann C M,Ainsworth T L,Gillis D B,et al.Modeling Coastal Waters from Hyperspectral Imagery Using Manifold Coordinates[C]∥IEEE International Conference on Geoscience & Remote Sensing Symposium.IEEE,2006:356-359.

[5] Belkin M,Niyogi P.Laplacian Eigenmaps and Spectral Techniques for Embedding and Clustering[J].Advances in Neural Information Processing Systems,2002,1(14):585-591.

[6] 齊濱.高光譜圖像分類及端元提取方法研究[D].哈爾濱:哈爾濱工程大學,2012:13-34.

(Qi Bin.Research on Hyperspectral Image Classification and Endmember Extraction Methods[D].Harbin: Harbin Engineering University,2012:13-34.)

[7] Zhang T H,Yang J,Zhao D L,et al.Linear Local Tangent Space Alignment and Application to Face Recognition[J].Neurocomputing,2007,70(7):1547-1553.

[8] He X F,Cai D,Yan S C,et al.Neighborhood Preserving Embedding[J].IEEE International Conference on Computer Vision,2005,2:1208-1213.

[9] He X F,Niyogi P.Locality Preserving Projections[J].Advances in Neural Information Processing Systems,2004,2(16):153-160.

[10] 王立志.基于流形學習的高光譜圖像降維與分類研究[D].重慶:重慶大學,2012:34-56.

(Wang Lizhi.Research on Dimensionality Reduction and Classification of Hyperspectral Image Based on Manifold Learning[D].Chongqing: Chongqing University,2012:34-56.)

[11] 甘資先,周方俊,肖奕.多維尺度分析中的算法研究[J].清華大學學報:自然科學版,1991,31(6):20-27.

(Gan Zixian,Zhou Fangjun,Xiao Yi.Study on the Algorithm for Multidimensional Scaling[J].Journal of Tsinghua University: Natural Science,1991,31(6):20-27.)

【責任編輯: 李 艷】

Hyperspectral Image Feature Extraction Algorithm Based on Global and Local Manifold Structure

ZhaoChunhui,CuiXiaochen

(School of Information and Communication Engineering,Harbin Engineering University,Harbin 150001,China)

High spectrum image manifold learning unsupervised feature extraction algorithm can only separate description based on high dimensional data space whether it is local or global geometric structure or not,there is no algorithm can also maintain a high dimensional data of global and local geometric structure.A novel unsupervised feature extraction algorithm is proposed based on global and local manifold structure (GLMS) for feature extraction of hyperspectral image.Manifold learning algorithms are based on the basic theory,which establishes two maintain manifold structure neighbor graphs,and it is applied to describe the data of the global and local manifold structure,by solving the generalized eigenvalue problem to obtain the optimal projection reconstruction weight matrix and then to get the low dimensional embedding space to achieve the purpose of reducing dimension.In the AVIRIS hyperspectral image tests on Indian Pine and Salina data set,experimental results show that the proposed algorithm has better improvement in classification accuracy and computational efficiency.

hyperspectral image; unsupervised feature extraction; global and local manifold; manifold learning

2015-01-15

趙春暉(1965-),男,黑龍江湯原人,哈爾濱工程大學教授,博士生導師.

2095-5456(2015)04-0283-06

TN 911.2

A