基于分數群稀疏混合范式和空間正則化的高光譜解混

王 傘, 王立國, 王麗鳳

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

0 引 言

高光譜遙感技術廣泛應用于環境監測、精準農業、林業監測和礦產勘探等領域[1]。由于成像傳感器的空間分辨率有限,觀察到的像元常由多種類型的地物組成,這種像元稱為混合像元。高光譜解混合旨在估計混合像元中的純光譜信號(端元)和相應比例(豐度)[2]。在高光譜解混研究中應用最廣泛的光譜混合模型為線性混合模型(Linear sectral mixing model, LSMM),其中混合像元的反射率等于其相應的豐度和端元線性加權組合。利用LSMM,基于幾何、統計和稀疏的高光譜解混方法逐步發展起來[3]。隨著美國地質調查局(United States Geological Survey, USGS)光譜庫的廣泛使用,基于稀疏的高光譜解混方法已成為研究熱點。高光譜稀疏解混方法旨在找到一個光譜庫的最佳子集對高光譜圖像的混合像元進行建模[4]。文獻[5]首次采用l1范式來誘導豐度的稀疏性,提出了基于變量分裂和增廣拉格朗日的稀疏解混算法(Sparse unmixing by variable splitting and augmented Lagrangian, SUnSAL)。為了結合空間信息,文獻[6]將全變分(Total variation, TV)正則項引入SUnSAL,提出了SUnSAL-TV算法。TV正則化揭示了空間均質性,即兩個相鄰像元中含有同類端元的豐度相似。文獻[7]通過將TV正則化和加權非局部低秩張量正則化引入到協作稀疏變分增廣拉格朗日解混方法[8](Collaborative SUnSAL, CLSUnSAL),提出了一種加權非局部低秩張量稀疏解混方法(Weighted nonlocal low-rank tensor decomposition method for sparse unmixing, WNLTDSU),同時利用了局部空間平滑性、全局稀疏性和非局部低階性。SUnSAL-TV、CLSUnSAL和WNLTDSU假設一個端元代表一種類型的材料,忽略了端元光譜可變性。解決這個問題的一個更為理想的方法是引入群稀疏性;即光譜庫被組織成多個端元群,每個端元群代表一種類型的地物。解混時,混合像元中只有少數端元處于活躍狀態;因此,估計的豐度是稀疏的,并且具有群結構。文獻[9]提出了三種群稀疏混合范式lG,2,1、lG,1,2和lG,1,q。lG,2,1范式稱為Group Lasso(G-Lasso),它誘導少數豐度群是活躍的,而這些活躍豐度群中的大部分豐度都是活躍的。lG,1,2范式稱為Elitist Lasso(E-Lasso),它誘導多數豐度群活躍,而這些活躍的豐度群體中少數豐度是活躍的[9]。lG,1,q范式稱為Fractional Lasso(F-Lasso),它誘導少數豐度群和每個群中的少數豐度活躍[9]。當端元和端元群的數量較大時(即場景中有多種類型的地物),lG,1,q范式的性能優于lG,2,1和lG,1,2范式。群稀疏混合范式針對具有群結構的光譜庫進行豐度稀疏誘導,獲得的豐度具有群結構,即每個端元群包含的多端元均參與解混并具有對應豐度,群稀疏混合范式能較好地將端元光譜可變性納入解混模型;但群稀疏混合范式沒有考慮豐度空間均質性特征對豐度的約束,導致估計的豐度欠平滑,因而將TV正則化引入到群稀疏解混,有助于提高豐度估計精度。

本文將TV空間正則化引入到分數混合范式lG,1,q中,提出了一種基于圖像的光譜庫提取、分數群稀疏混合范數和TV空間正則化的高光譜解混算法(Hyperspectral unmixing algorithm based on image spectral library extraction, fractional group sparse mixing norm, and TV space regularization, HU-ISLE-FTV)。采用一種新的基于圖像的光譜庫提取方法,稱為基于頂點分量分析(Vertex component analysis, VCA)和概率輸出支持向量機(Probabilistic outputs support vector machine, PSVM)的自動光譜庫提取算法(Automatic spectral library extraction algorithm based on VCA and PSVM, ASLE-VPSVM)。傳統稀疏解混方法利用通用光譜庫解混,容易產生光譜庫與待解混數據失配問題,同時,通用光譜庫中通常利用單個端元光譜代表一類地物,未能很好地體現端元光譜可變性。ASLE-VPSVM針對待解混數據提取端元,克服了稀疏解混算法對通用光譜庫的依賴性,提高了它們的適用性。同時,ASLE-VPSVM提取多個端元代表一類地物,充分反映了客觀存在端元光譜可變性,能夠提高稀疏解混方法的解混性能。此外,ASLE-VPSVM可以用作無監督或半監督方法。傳統稀疏解混算法注重豐度稀疏性誘導,忽略了豐度空間均質性,即兩個相鄰像元中含有的同類端元的豐度是相似的。HU-ISLE-FTV算法首先應用ASLE-VPSVM提取多端元作為光譜庫,然后將TV空間正則化引入到F-Lasso混合范式中,在強化稀疏誘導的同時,施加豐度空間均質性約束,有效地提高豐度估計精度。同時,為HU-ISLE-FTV算法設計了一種高效簡潔的交替方向乘子迭代法(Alternating direction method of multipliers, ADMM)[10],提高了算法運算效率。

1 相關工作

1.1 ASLE-VPSVM算法

以提取3類端元為例表示ASLE-VPSVM算法原理,如圖1所示。ASLE-VPSVM算法首先從原始高光譜圖像中隨機選擇像元構成隨機圖像子集,然后在該子集運用VCA算法[11]提取端元m1、m2和m3,并利用PSVM[12]判斷提取端元(概率值為1)的鄰域像元的概率估計值,該估計值體現了端元鄰域像元歸屬于端元對應所屬類別的可能性,概率估計值越大,則端元鄰域像元為同類端元的可能性越大。ASLE-VPSVM方法在針對每個隨機圖像子集提取端元時,包含2次提取過程,過程1利用VCA提取端元,過程2利用PSVM判斷過程1中已提取端元的鄰域像元是否為同類端元,如符合條件,則提取為端元?？梢?ASLE-VPSVM既充分利用了光譜特征空間,即采用VCA方法,同時充分利用空間信息,即利用PSVM判斷端元空間鄰域像元為端元的可能性。

圖1 ASLE-VPSVM算法原理圖

ASLE-VPSVM算法具體實現過程如下:

步驟1: 從原始高光譜圖像中隨機選擇一定比例的像元構成圖像子集。在該圖像子集中,運用VCA算法提取P個端元及其位置索引。P為原始高光譜圖像中存在的端元類別數量,端元的位置索引是指其在原始圖像中的位置坐標。

步驟2: 根據步驟1中提取的端元位置索引找到提取的P個端元的空間鄰域像元(采用4鄰域像元)作為候選端元。將已提取端元用作訓練集,使用PSVM算法[12]預測候選端元概率。當候選端元的概率值大于閾值時,將其提取為真實端成員并添加到已提取端元。已提取端元是由步驟1和步驟2中提取的端元構成。

步驟3: 重復步驟1和2若干次后終止迭代。然后,從提取端元集合中移除冗余樣本,最終獲得光譜庫。ASLE-VPSVM可用作無監督方法或半監督方法。當ASLE-VPSVM用作無監督方法時,在步驟1中,應首先應用HySime算法[13]估計端元種類數量P,并通過VCA提取端元。然后將已提取端元的平均向量用作參考向量,通過計算光譜角度距離區分新提取的P端元屬于哪一類端元。無監督ASLE-VPSVM可以區分端員類別,但不知道每種端元類別的具體屬性。在通常情況下,P和參考向量已知,ASLE-VPSVM則可用作半監督方法。

1.2 HU-ISLE-FTV算法

令Y∈L×N代表一個具有L個光譜波段且包含N個像元的高光譜圖像數據,M∈L×n代表一個具有群結構且包含n個端元的端元光譜庫;A∈n×N代表豐度矩陣,其中每一列對應Y中一個像元的豐度向量a;則線性光譜混合模型可表述為:

Y=MA

(1)

當M由ASLE-VPSVM算法獲得后,利用F-LASSO建立稀疏解混優化問題如下:

(2)

lG,1,q可以被視為l1和lq的一種“組合”。如圖2(a)所示,群結構被定義為豐度向量的P個子向量, 圖2(a)中每個黑框對應一個豐度子向量,每個子向量內的豐度構成一個豐度群,群中綠色、藍色、紅色和黃色代表數值不為0的豐度,其余淡色代表數值為0的豐度。lG,1,q首先計算每個子向量的l1范式以獲得P維向量,然后計算該P維向量的lq范式以獲得最終結果。在豐度非負的約束條件下,計算每個子向量的l1范式相當于計算子向量中所有豐度的總和。如圖2(b)所示,lG,1,q范式可以很容易地擴展到豐度矩陣,按列計算l1范式并對結果計算lq范式。不失一般性,假設a為某一像元對應的豐度向量,則:

(a) lG,1,q作用于豐度向量 (b) lG,1,q作用于豐度矩陣

a=[a1,1,a1,2,…,a1,m1,a2,1,a2,2,…,a2,m2,…,aP,1,aP,2,…,aP,mP]T

(3)

由于a內每一個豐度都非負,那么求解a的lG,1,q范式相當于a的每個群內豐度相加得到一個P維向量,然后再求這個P維向量的lq范式。若定義G∈P×n為:

(4)

式中:第g行包含mg個連續的1,mg為第g個群包含豐度個數,那么式(3)等價為:

(5)

將TV正則化約束引入到F-LASSO解混模型,得到基于分數混合范式和總變分約束的解混模型如下:

(6)

定義Hh為一個線性算子,Hh:P×N表示計算當前像元與其空間垂直鄰域像元的豐度向量之差,如Hh(GA)=[d1,d2,…,dN],式中di=ai-aih,i和ih分別表示某一像元及其垂直鄰域。同樣,定義Hν為一個線性算子,Hν:P×N表示計算當前像元與其空間水平鄰域像元的豐度向量之差,如Hν(GA)=[d1,d2,…,dN],式中di=ai-aiν,i和iν分別表示某一像元及其水平鄰域。

根據上述兩個線性算子進一步定義:

(7)

根據分數混合范式和TV正則定義,式(6)的等效優化問題如下:

(8)

利用變量分裂的方法給出式(8)的約束等價公式:

(9)

式(9)的增廣拉格朗日方程為:

(10)

利用ADMM將式(10)的優化過程分解為若干個簡化的子問題進行求解。

(1) 變量A的優化子問題如下:

(11)

令式(11)中的優化目標函數對A求導并等于0,可得:

A←(MTM+2μGTG+μI)-1 (MTY+μGT(V1+D1+V2+D2)+μ(V4+D4))

(12)

(2) 變量V1的優化子問題如下:

(13)

參考文獻[9]解式(13),可得:

V1←Sq,τ(GA-D1)

(14)

(3) 變量V2的優化子問題如下:

(15)

令式(15)中的優化目標函數對V2求導并等于0,可得:

V2←(HTH+I)-1(GA-D2+HT(V3+D3))

(16)

(4) 變量V3的優化子問題如下:

(17)

應用矢量軟閾值[9]求解式(17),可得:

V3←Softτ(HV2-D3)

(18)

(5) 變量V4的優化子問題如下:

(19)

根據文獻[16]解式(19),可得:

V4←projΔ(A-D4)

(20)

(6) 更新拉格朗日乘子D:

(21)

綜上分析,總結HU-ISLE-FTV算法的ADMM如下:

步驟1: 輸入高光譜圖像數據Y,應用ASLE-VPSVM提取端元光譜庫M;

步驟2: 利用全約束最小二乘解混方法(Fully constrained least square method, FCLSM)獲得初始化豐度矩陣A,并進一步利用式(9)中約束條件計算初始化V;同時初始化D=0;

步驟3: 利用式(12)更新A;

步驟4: 利用式(14)、式(16)、式(18)和式(20)分別更新V1、V2、V3和V4;

步驟5: 用式(21)更新D1、D2、D3和D4;

步驟6: 循環執行步驟3至步驟5,直到滿足迭代停止條件。

2 實驗結果與分析

2.1 實驗設置

2.1.1 實驗數據

采用一個人工模擬數據集(DC1)和一個真實高光譜數據集Houston進行實驗。模擬數據集DC1采用線性混合模型,由3類真實光譜和3個模擬豐度圖人工線性合成而得到。3類真實光譜選自于USGS數字光譜庫,分別為赤鐵礦(Hematite)、白云母(Muscovite)和橄欖石(Olivine),上述光譜包含224個波段,均勻分布在0.4～2.5 μm。三類光譜分別包含12、13和17個端元,這樣選取端元的目的是真實模擬高光譜數據的端元光譜可變性。模擬豐度圖依據文獻[6]中的方法生成,如圖3第一列所示。每次從3類真實光譜中各隨機抽取一個端元并結合豐度圖中每個像元對應豐度,線性合成混合像元,最終得到63×63×224模擬數據集DC1。真實高光譜數據集Houston為2012年6月在美國德克薩斯州休斯頓大學校園拍攝的圖像的子集,該圖像包含區間0.380～1.050 μm的144個波段,具有2.5 m的空間分辨率。截取原數據集的105×128×144子集命名為Houston,該數據集包含4類不同地物:植被(Vegetation)、紅色屋頂(Red roofs)、混凝土(Concrete)和瀝青(Asphalt)。

圖3 7種算法估計DC1數據中3類端元豐度圖

2.1.2 實驗算法及參數設置方法

參與實驗的算法有FCLSU、SUnSAL-TV、G-LASSO、E-LASSO、F-LASSO、WNLTDSU和HU-ISLE-FTV等7種算法。其中,G-LASSO和E-LASSO需要設置的參數為μ和λ;SUnSAL-TV和WNLTDSU需要設置的參數為μ、λ和λTV;F-LASSO需要設置的參數為μ、λ和q;GSME-FTV需要設置的參數為μ、λ、λTV和q。μ設置為固定值0.05,λ和λTV在以下數值中進行取值:{10-5,5×10-4,10-3,5×10-3,0.01,0.05,0.1,0.3,0.5,1.00}。q在3個數值范圍內取值,分別為0.001～0.01、0.01～0.1、0.1～1;每個數值區間包含10個等間距數值。上述參數遍歷所有組合,可尋找最優參數。

2.1.3 實驗結果評估指標

實驗采用均方根誤差(Root mean square error, RMSE)、光譜角距離(Spectral angular distanceroot, SAD)和重構誤差(Reconstruction error, RE)作為量化指標評估解混結果,具體計算方式如下:

(22)

2.2 DC1數據實驗

參與比較的7種算法利用ASLE-VPSVM所提取的相同光譜庫進行解混。針對DC1數據,本實驗運行ASLE-VPSVM 60次迭代,獲得8個赤鐵礦端元、9個白云母端元和7個橄欖石端元。圖3顯示了7種算法解混獲得的豐度圖,各種方法的量化評估結果如表1所示。圖3中每一行對應一類端元,由上至下分別對應赤鐵礦、白云母和橄欖石;每一列對應一種算法,由左至右分別對應真實豐度、FCLS、SUnSAL-TV、G-LASSO、E-LASSO、F-LASSO、WNLTDSU和HU-ISLE-FTV。由圖3可以看出,與真實豐度圖相比,7種方法估計的豐度圖均不同程度包含噪聲點,這是因為DC1中每個像元由3個隨機選擇的端元線性合成;端元光譜可變性很大,這影響了各種方法的解混性能。然而,7種算法均獲得了較為清晰的豐度圖,這是因為實驗使用的ASLE-VPSVM提取的光譜庫具有群結構,即每種類型的材料都由多個端元代表,可以更好地反映端元光譜可變性,從而驗證了所提出的ASLE-VPSVM的有效性。圖3顯示,FCLSU和WNLTDSU估算的豐度圖比其他五種方法估算的豐度圖更為模糊,尤其是赤鐵礦和橄欖石豐度圖。SUnSAL-TV、G-Lasso、E-Lasso和F-Lasso估算的豐度圖相似。與其他方法估計的豐度圖相比,HU-ISLE-FTV估算的豐度圖與真實豐度圖更為吻合。表1表明,HU-ILSE-FTV在RMSE、RE和運行時間方面均優于其他算法,但SAM略大于其它大部分算法,HU-ISLE-FTV的運行時間明顯短于其他稀疏解混算法,證實HU-ISLE-FTV設計的ADMM效率比較高。綜上所述,HU-ISLE-FTV解混合成數據的性能優于其他參與比較的算法。

表1 7種算法解混DC1數據量化評估結果

2.3 Hunston數據實驗

針對Hunston數據,運行了ASLE-VPSVM 60次迭代,獲得15個植被端元、23個紅色屋頂端元、21個混凝土端元和19個瀝青端元。7種算法估計的豐度圖如圖4所示,表2列出了各種方法的量化評估結果。圖中每一行對應一類端元,由上至下分別對應植被、紅色屋頂、混凝土和瀝青;每一列對應一種算法,由左至右分別對應FCLSM、SUnSAL-TV、G-LASSO、E-LASSO、F-LASSO、WNLTDSU和HU-ISLE-FTV。從視覺上看,幾乎所有豐度圖都清楚地顯示了各類端元的分布。這驗證了ASLE-VPSVM提取的光譜庫能夠很好地反映端元光譜可變性,驗證了HU-ISLE-FTV算法的有效性。定量評估結果如表2所示,對于真實高光譜數據,由于無法獲得其真實的豐度和端元,因此使用RE和運行時間作為評估指標。在所有方法中,HU-ISLE-FTV的RE最小,HU-ISLE-FTV的運行時間略高于FCLSU,遠低于其他方法。綜上所述,HU-ISLE-FTV解混真實數據的性能優于其他參與比較的算法。

表2 7種算法解混Hunston數據量化評估結果

3 結 論

提出了HU-ISLE-FTV算法,該算法設計了一種基于VCA和PSVM的自動光譜庫提取方式,稱為ASLE-VPSVM。ASLE-VPSVM提取多個端元代表一類地物,較好地解決了端元光譜可變性對高光譜解混的影響,并使稀疏解混算法不再依賴于通用光譜庫。HU-ISLE-FTV將TV空間正則化引入分數群稀疏混合范數,在強化豐度稀疏誘導的同時施加空間均質性約束,有效地提升了豐度估計精度。對人工模擬數據集和真實數據集進行實驗,結果表明HU-ISLE-FTV性能優于現有的6種代表性稀疏解混算法。未來的研究,將集中于將所提出的方法應用于非線性解混模型。