基于分?jǐn)?shù)群稀疏混合范式和空間正則化的高光譜解混

王 傘, 王立國, 王麗鳳

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

0 引 言

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

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

1 相關(guān)工作

1.1 ASLE-VPSVM算法

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

圖1 ASLE-VPSVM算法原理圖

ASLE-VPSVM算法具體實(shí)現(xiàn)過程如下:

步驟1: 從原始高光譜圖像中隨機(jī)選擇一定比例的像元構(gòu)成圖像子集。在該圖像子集中,運(yùn)用VCA算法提取P個(gè)端元及其位置索引。P為原始高光譜圖像中存在的端元類別數(shù)量,端元的位置索引是指其在原始圖像中的位置坐標(biāo)。

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

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

1.2 HU-ISLE-FTV算法

令Y∈L×N代表一個(gè)具有L個(gè)光譜波段且包含N個(gè)像元的高光譜圖像數(shù)據(jù),M∈L×n代表一個(gè)具有群結(jié)構(gòu)且包含n個(gè)端元的端元光譜庫;A∈n×N代表豐度矩陣,其中每一列對(duì)應(yīng)Y中一個(gè)像元的豐度向量a;則線性光譜混合模型可表述為:

Y=MA

(1)

當(dāng)M由ASLE-VPSVM算法獲得后,利用F-LASSO建立稀疏解混優(yōu)化問題如下:

(2)

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

(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內(nèi)每一個(gè)豐度都非負(fù),那么求解a的lG,1,q范式相當(dāng)于a的每個(gè)群內(nèi)豐度相加得到一個(gè)P維向量,然后再求這個(gè)P維向量的lq范式。若定義G∈P×n為:

(4)

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

(5)

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

(6)

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

根據(jù)上述兩個(gè)線性算子進(jìn)一步定義:

(7)

根據(jù)分?jǐn)?shù)混合范式和TV正則定義,式(6)的等效優(yōu)化問題如下:

(8)

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

(9)

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

(10)

利用ADMM將式(10)的優(yōu)化過程分解為若干個(gè)簡化的子問題進(jìn)行求解。

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

(11)

令式(11)中的優(yōu)化目標(biāo)函數(shù)對(duì)A求導(dǎo)并等于0,可得:

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

(12)

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

(13)

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

V1←Sq,τ(GA-D1)

(14)

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

(15)

令式(15)中的優(yōu)化目標(biāo)函數(shù)對(duì)V2求導(dǎo)并等于0,可得:

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

(16)

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

(17)

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

V3←Softτ(HV2-D3)

(18)

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

(19)

根據(jù)文獻(xiàn)[16]解式(19),可得:

V4←projΔ(A-D4)

(20)

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

(21)

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

步驟1: 輸入高光譜圖像數(shù)據(jù)Y,應(yīng)用ASLE-VPSVM提取端元光譜庫M;

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

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

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

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

步驟6: 循環(huán)執(zhí)行步驟3至步驟5,直到滿足迭代停止條件。

2 實(shí)驗(yàn)結(jié)果與分析

2.1 實(shí)驗(yàn)設(shè)置

2.1.1 實(shí)驗(yàn)數(shù)據(jù)

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

圖3 7種算法估計(jì)DC1數(shù)據(jù)中3類端元豐度圖

2.1.2 實(shí)驗(yàn)算法及參數(shù)設(shè)置方法

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

2.1.3 實(shí)驗(yàn)結(jié)果評(píng)估指標(biāo)

實(shí)驗(yàn)采用均方根誤差(Root mean square error, RMSE)、光譜角距離(Spectral angular distanceroot, SAD)和重構(gòu)誤差(Reconstruction error, RE)作為量化指標(biāo)評(píng)估解混結(jié)果,具體計(jì)算方式如下:

(22)

2.2 DC1數(shù)據(jù)實(shí)驗(yàn)

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

表1 7種算法解混DC1數(shù)據(jù)量化評(píng)估結(jié)果

2.3 Hunston數(shù)據(jù)實(shí)驗(yàn)

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

表2 7種算法解混Hunston數(shù)據(jù)量化評(píng)估結(jié)果

3 結(jié) 論

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