超像素低秩高光譜稀疏解混

李 璠

(南昌工程學(xué)院 江西省水信息協(xié)同感知與智能處理重點(diǎn)實(shí)驗(yàn)室，江西 南昌 330099)

受遙感器空間分辨力限制，高光譜圖像中存在大量的混合像元，阻礙了高光譜數(shù)據(jù)的精細(xì)化解譯[1-2]。高光譜解混旨在提取高光譜圖像中的純凈光譜(端元提取)并估計(jì)每種端元在混合像元中的占比(豐度估計(jì))[3]。假設(shè)所有端元光譜都存在于一個(gè)大的可用光譜庫(kù)中，則解混過(guò)程可等同于在光譜庫(kù)中尋找可以有效表達(dá)混合像元的最佳光譜特征子集[4]。由于一幅圖像中所含端元數(shù)量通常遠(yuǎn)少于光譜庫(kù)中的端元個(gè)數(shù)，豐度系數(shù)自然表現(xiàn)出稀疏特性，因此，稀疏線性回歸算法常用于解決這一優(yōu)化問(wèn)題，稱為稀疏解混[5]。例如，基于變量分裂和增廣拉格朗日的稀疏解混算法(Sparse Unmixing algorithm via variable Splitting and Augmented Lagrangian，SUnSAL)采用l1范數(shù)約束豐度系數(shù)的稀疏性[4]，協(xié)同SUnSAL算法(Collaborative SUnSAL，CLSUnSAL)采用l2，1范數(shù)約束豐度矩陣的行稀疏性[6]。此外，迭代加權(quán)稀疏策略通過(guò)懲罰豐度系數(shù)中的非零值來(lái)增強(qiáng)解的稀疏性，取得了良好效果，如雙重加權(quán)稀疏解混算法(Double Reweighted Sparse Unmixing，DRSU)在l1范數(shù)正則項(xiàng)中引入兩種加權(quán)因子促進(jìn)豐度矩陣在光譜域和空間域的稀疏性[7]。這些算法只考慮豐度系數(shù)的稀疏特征，而忽視了像元的空間結(jié)構(gòu)。

利用高光譜圖像豐富的空間信息可以有效約束豐度估計(jì)優(yōu)化問(wèn)題的解空間，因此一系列空間先驗(yàn)稀疏解混方法被提出[8]。例如，基于變量分裂增廣拉格朗日和全變差稀疏解混算法(Sparse Unmixing via variable Splitting Augmented Lagrangian and Total Variation，SUnSAL-TV)通過(guò)各向異性全變差正則項(xiàng)刻畫相鄰像元之間的平滑性[9]，譜—空加權(quán)稀疏解混算法(Spectral-Spatial Weighted Sparse Unmixing，S2WSU)利用包含像元鄰域信息的空間加權(quán)因子促進(jìn)圖像的連續(xù)性[10]，非局部稀疏解混算法(Non-Local Sparse Unmixing，NLSU)采取空間結(jié)構(gòu)自相似正則化方法挖掘圖像的空間高階結(jié)構(gòu)信息[11]。融入空間信息的稀疏解混方法在解混性能方面表現(xiàn)出巨大的潛力，但采用固定形狀和大小的窗口探索像元之間的相關(guān)性限制了對(duì)圖像空間結(jié)構(gòu)信息的精確描述。

高光譜圖像中相似像元的豐度向量具有相關(guān)性，反之，對(duì)豐度矩陣施加低秩約束可以保持圖像的空間一致性[12]。例如，交替方向稀疏低秩解混算法(Alternating Direction Sparse and Low-Rank Unmixing，ADSpLRU)和聯(lián)合稀疏塊低秩解混算法(Joint-Sparse-Blocks and Low-Rank Unmixing，JSpBLRU)分別利用3×3滑動(dòng)窗口和全局的豐度加權(quán)核范數(shù)實(shí)現(xiàn)豐度系數(shù)的低秩逼近，以保留高光譜圖像的低維空間結(jié)構(gòu)[13-14]。事實(shí)上，相對(duì)于預(yù)定義窗口或整幅圖像，局部同質(zhì)區(qū)域內(nèi)的像元具有更高的相似度。超像素分割技術(shù)以像元的光譜特征相似度和空間距離為度量，將像元聚類為不同的空間組，自適應(yīng)生成互不交疊的不規(guī)則同質(zhì)區(qū)域，每個(gè)區(qū)域稱為一個(gè)超像素[15]。超像素分割已廣泛應(yīng)用于高光譜圖像處理，如分類[16]、去噪[17]、超分辨率[18]等。

綜上，受低秩表示理論和迭代加權(quán)稀疏思想的啟發(fā)，結(jié)合超像素分割技術(shù)，本文提出一種超像素低秩稀疏解混算法(SupLRSU)，該算法將基于超像素的低秩約束和光譜加權(quán)稀疏約束聯(lián)合施加在豐度矩陣上，既保持了圖像內(nèi)在的局部低維空間結(jié)構(gòu)，促進(jìn)圖像的空間一致性，又誘導(dǎo)豐度矩陣的行稀疏性，促使圖像中所有像元的豐度向量聯(lián)合稀疏。

1 相關(guān)工作

1.1 線性光譜混合模型

線性光譜混合模型假設(shè)每個(gè)像元光譜都可以表示為端元光譜的線性組合[3]，其矩陣形式記為

Y=MX+N， s.t.X≥0，1Tx=1.

(1)

其中，Y=[y1，y2，…yn]∈d×n表示高光譜圖像，d為波段數(shù)，n為像元數(shù)，M=[m1，m2，…，ml]∈d×l表示端元矩陣，l為端元數(shù)，X=[x1，x2，…xn]∈l×n表示豐度矩陣，N∈d×n表示噪聲或模型誤差，X≥0為豐度“非負(fù)性”約束(Abundance Nonnegative Constraint，ANC)，1Tx=1為豐度“和為一”約束(Abundance Sum-to-one Constraint，ASC)。

1.2 高光譜稀疏解混

稀疏解混用已知光譜庫(kù)A∈d×m代替端元矩陣M∈d×l。同一場(chǎng)景中存在的地物種類一般少于20種，遠(yuǎn)小于光譜庫(kù)中光譜特征的數(shù)量(l?m)，因此豐度矩陣X∈m×n表現(xiàn)出稀疏特性[4]。稀疏解混模型可表示為

(2)

式(2)中l(wèi)0范數(shù)的求解是典型的NP(Non-deterministic Polynomial)難問(wèn)題，在滿足有限等距性質(zhì)(Restricted Isometry Property，RIP)的條件下，l0與l1優(yōu)化問(wèn)題有相同的最優(yōu)解，式(2)可轉(zhuǎn)化為

(3)

1.3 基于低秩先驗(yàn)的高光譜稀疏解混

高光譜圖像鄰域像元通常由相同的幾種端元按照相同或相近的比例組合而成，因而豐度矩陣表現(xiàn)出低秩性。為捕獲這種低秩稀疏的數(shù)據(jù)結(jié)構(gòu)，在式(3)中加入豐度秩約束項(xiàng)，得到豐度估計(jì)的稀疏減秩回歸問(wèn)題[13]，表示如下：

(4)

其中，rank(X)表示X的秩，τ≥0是正則化參數(shù)。式(4)中rank(X)不易求解，常用其核范數(shù)‖X‖*替代，記為

(5)

其中，σi(X)表示X的第i個(gè)奇異值，r=rank(X)。加權(quán)核范數(shù)最小化方法(Weighted Nuclear Norm Minimization，WNNM)通過(guò)為每個(gè)奇異值分配不同的權(quán)值實(shí)現(xiàn)奇異值收縮[20]，其計(jì)算方式如下：

(6)

綜上，基于低秩先驗(yàn)的高光譜稀疏解混模型表示如下：

(7)

模型(7)通過(guò)最小化加權(quán)核范數(shù)實(shí)現(xiàn)豐度系數(shù)的低秩矩陣逼近，促使豐度向量線性相關(guān)，雖提高了解混精度，但忽略了圖像的空間細(xì)節(jié)信息。

2 本文方法

2.1 超像素分割

為更準(zhǔn)確地描述高光譜圖像的空間相關(guān)性，采用超像素分割技術(shù)將圖像劃分為若干個(gè)互不交疊的同質(zhì)區(qū)域。簡(jiǎn)單線性迭代聚類算法(Simple Linear Iterative Clustering，SLIC)是一種結(jié)構(gòu)簡(jiǎn)單、快速有效的超像素生成算法[21]。SLIC算法基本思想是局部區(qū)域k-均值(k-means)聚類，具體流程為：初始化p個(gè)種子點(diǎn)，并將搜索空間限定在以種子點(diǎn)為中心、與超像素大小成比例的區(qū)域，計(jì)算搜索區(qū)域內(nèi)各個(gè)像元到種子點(diǎn)的加權(quán)距離，將每個(gè)像元?dú)w類至距離最近的種子點(diǎn)。每輪歸類結(jié)束，以同一類像元的均值為新的種子點(diǎn)重復(fù)上述步驟，迭代至收斂。SLIC算法一方面限定了k-均值的搜索空間，減少了優(yōu)化過(guò)程中距離的計(jì)算次數(shù)，降低了算法的復(fù)雜度，另一方面以像元的顏色接近度與空間位置關(guān)系為指標(biāo)設(shè)置加權(quán)距離，靈活控制超像素的緊湊性和大小。

原始的SLIC算法主要面向CIELAB顏色空間的彩色圖像，需要修改加權(quán)距離的計(jì)算方式以適應(yīng)三維的高光譜數(shù)據(jù)[15]。為保持一個(gè)超像素區(qū)域內(nèi)像元的相似度，改進(jìn)的SLIC算法同時(shí)考慮像元的光譜距離和空間距離，定義中心像元(xi，yi)與搜索區(qū)域內(nèi)任一像元(xj，yj)的加權(quán)距離如下：

(8)

其中，zi和zj分別表示像元(xi，yi)和(xj，yj)的觀測(cè)光譜向量，dspe表示像元的光譜距離，采用光譜角距離(Spectral Angle Distance，SAD)度量，dspa表示像元的空間距離，采用歐式距離度量，常數(shù)ω>0表示空間距離和光譜距離之間的權(quán)重，D表示像元(xi，yi)和(xj，yj)的加權(quán)距離。

采用改進(jìn)的SLIC算法對(duì)高光譜圖像進(jìn)行空間預(yù)處理，將整幅圖像標(biāo)記為一組不規(guī)則的空間異構(gòu)區(qū)域，而每個(gè)區(qū)域內(nèi)像元的觀測(cè)光譜高度相似。

2.2 超像素低秩高光譜稀疏解混模型

為更有效地利用圖像的空間上下文信息，用每個(gè)超像素區(qū)域的豐度低秩約束取代全局性或固定窗口的低秩約束，僅使每個(gè)同質(zhì)區(qū)域內(nèi)的像元共享相同的稀疏結(jié)構(gòu)。此外，為促使圖像中所有像元共享光譜庫(kù)相同的活動(dòng)端元子集，我們還在稀疏正則項(xiàng)中引入光譜加權(quán)因子誘導(dǎo)豐度矩陣的行稀疏性?；谏鲜隹紤]，提出超像素低秩高光譜稀疏解混模型(SupLRSU)，表示如下：

(9)

其中，X(s)表示圖像中第s個(gè)超像素，p為超像素個(gè)數(shù)，H∈m×m是一個(gè)對(duì)角矩陣，表示光譜權(quán)重，H中第i行的值計(jì)算如下：

(10)

其中，X(i，：)表示X的第i行向量，‖·‖2表示l2范數(shù)。光譜權(quán)重H迭代更新，其每一行的權(quán)值與豐度矩陣相應(yīng)行向量的l2范數(shù)成反比，因此，每次迭代將更大程度地懲罰豐度系數(shù)中具有較小值的行，引導(dǎo)豐度矩陣只保留少數(shù)的非零行，提高豐度矩陣的行稀疏性。

2.3 模型求解

受文獻(xiàn)[10]啟發(fā)，我們采用內(nèi)外雙循環(huán)迭代方案求解式(9)模型，在外循環(huán)中根據(jù)式(10)更新光譜權(quán)重H，在內(nèi)循環(huán)中利用交替方向乘子法(Alternating Direction Method of Multipliers，ADMM)更新豐度矩陣X，具體過(guò)程如下。

引入5個(gè)輔助變量U，V1，V2，V3，V4，將式(9)改寫為如下形式：

(11)

(12)

其增廣拉格朗日函數(shù)為

(13)

其中，μ>0為罰參數(shù)，D=(D1，D2，D3，D4)為約束條件GU+BV=0相關(guān)的拉格朗日乘子。

根據(jù)ADMM方法，對(duì)式(13)中U和V分別迭代求解下列優(yōu)化問(wèn)題，更新拉格朗日乘子D：

(14)

計(jì)算式(14)得到各個(gè)變量的迭代更新公式。SupLRSU算法流程如下所示。其中，sof r(α，β)≡sign(α)max(|α|-β，0)表示軟閾值函數(shù)，svtb，β(X)≡U(σ(X)-βdiag((b))+VT表示加權(quán)奇異值閾值算子。

算法1:SupLRSU算法1.Initialization:set t,k,s=0,λ,τ,μ,ε>0,U(0),V(0)1,V(0)2,V(0)3,V(0)4,D(0)1,D(0)2,D(0)3,D(0)42.Repeat:3.Update the spectral weighting facor:h(t+1)=1‖(U(t)-D(t)2)(i,:)‖2+ε,i=1,…,m4.Repeat:5.U(k+1)←(ATA+3I)-1(AT(V(k)1+D(k)1)+V(k)2+D(k)2+V(k)3+D(k)3+V(k)4+D(k)4)6.V(k+1)1←11+μ[Y+μ(AU(k+1)-D(k)1)]7.V(k+1)2←softU(k+1)-D(k)2,λμH(t)()8.for s=1:p9.(V(k+1)3)(s)←svtb,τμ((U(k+1)-D(k)3)(s))10.end for11.V(k+1)4←max(U(k+1)-D(k)4,0)12.Update the Lagrange multipliers:13.D(k+1)1←D(k)1-AU(k+1)1+V(k+1)114.D(k+1)2←D(k)2-U(k+1)2+V(k+1)215.D(k+1)3←D(k)3-U(k+1)3+V(k+1)316.D(k+1)4←D(k)4-U(k+1)4+V(k+1)417.Update iteration:k←k+118.Until stopping condition19.U(t+1)←U(k)20.D(t+1)←D(k)221.Update iteration:t←t+122.Until stopping condition

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

為了驗(yàn)證本文所提超像素低秩高光譜稀疏解混算法(SupLRSU)的性能，本節(jié)采用模擬和真實(shí)高光譜數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，與SUnSAL，CLSUnSAL，DRSU，ADSpLRU，JSpBLRU算法進(jìn)行對(duì)比，分析實(shí)驗(yàn)結(jié)果。模擬實(shí)驗(yàn)采用信號(hào)與重建誤差比(Signal Reconstruction Error，SRE)和豐度重構(gòu)正確率(Probability of Success，Ps)進(jìn)行定量評(píng)價(jià)。SRE(dB)定義如下：

(15)

(16)

3.1 模擬數(shù)據(jù)實(shí)驗(yàn)

從美國(guó)地質(zhì)勘探局(United States Geological Survey，USGS)發(fā)行的splib06波譜數(shù)據(jù)庫(kù)中隨機(jī)選擇240條端元光譜組成光譜庫(kù)A∈224×240(光譜波段數(shù)為224)。從光譜庫(kù)A中隨機(jī)選取9種端元，按照M.D.Iordache等制作的豐度分布圖線性組合生成模擬數(shù)據(jù)。隨后加入信噪比(SNR)為30dB、40dB和50dB的高斯噪聲，形成三組模擬數(shù)據(jù)集。實(shí)驗(yàn)涉及算法的正則化參數(shù)均已調(diào)優(yōu)，最大迭代次數(shù)均設(shè)為500次，殘差閾值均設(shè)為1e-5。9種端元的真實(shí)豐度分布如圖1所示。

圖1 模擬數(shù)據(jù)中9種端元的真實(shí)豐度圖

表1展示了不同解混算法對(duì)三組模擬數(shù)據(jù)集進(jìn)行豐度估計(jì)所獲得的SRE(dB)和Ps值。從表1可以看出，在低SNR情況下，所提SupLRSU算法的SRE(dB)值和Ps值均優(yōu)于其他算法，并展現(xiàn)出顯著優(yōu)勢(shì)，證明SupLRSU算法具有良好的抗噪性能。在高SNR情況下，所有算法的Ps值均大幅提升，SupLRSU算法的SRE(dB)值略高于其他算法或與其他算法基本持平。

表1 不同解混算法豐度估計(jì)結(jié)果的SRE(dB)和Ps值(30次平均)及相應(yīng)正則化參數(shù)設(shè)置

圖2展示了模擬數(shù)據(jù)真實(shí)豐度圖和SNR=30dB時(shí)不同解混算法的估計(jì)豐度圖。為使視覺(jué)效果更加清晰，從結(jié)果中隨機(jī)選擇500個(gè)像元進(jìn)行展示。此外，從表1可知SUnSAL和CLSUnSAL算法的豐度估計(jì)結(jié)果較差，因而這兩種算法的估計(jì)豐度圖在本文中不做展示。從圖2可以看出，每種算法基本都能從光譜庫(kù)中識(shí)別出9種端元，但所提SupLRSU算法對(duì)端元的識(shí)別更加準(zhǔn)確。SUnSAL-TV估計(jì)豐度圖中有明顯的代表錯(cuò)誤端元的異常豐度線，DRSU估計(jì)豐度圖中有較多異常點(diǎn)，ADSpLRU和JSpBLRU的豐度估計(jì)圖包含少量具有較小值的異常豐度線。

圖2 真實(shí)豐度圖及不同解混算法獲得的估計(jì)豐度圖(SNR=30dB)

為直觀比較豐度估計(jì)的精確性，圖3展示了SNR=30dB時(shí)不同解混算法獲得的端元9的豐度分布圖，其他端元豐度圖的表現(xiàn)與端元9基本一致，因此不做展示。從估計(jì)豐度與真實(shí)豐度的差值圖很容易看出，所提SupLRSU算法的豐度圖與端元9的真實(shí)豐度圖最接近。SUnSAL-TV得到的豐度圖出現(xiàn)了較明顯的過(guò)平滑現(xiàn)象，DRSU的豐度圖含有大量噪聲點(diǎn)，與同樣引入低秩約束的ADSpLRU和JSpBLRU算法相比，SupLRSU算法的估計(jì)豐度圖在保證分段平滑的同時(shí)保持了空間異質(zhì)性結(jié)構(gòu)，并且保留了更多的細(xì)節(jié)信息。

圖3 不同解混算法獲得的端元9的豐度圖以及與真實(shí)豐度的差值圖(SNR=30dB)

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

真實(shí)數(shù)據(jù)實(shí)驗(yàn)采用機(jī)載可見(jiàn)紅外成像光譜儀(Airborne Visible Infrared Imaging Spectrometer，AVIRIS)采集的Cuprite 礦區(qū)數(shù)據(jù)。該數(shù)據(jù)集光譜波段數(shù)為224，波長(zhǎng)范圍為0.4～2.5μm，光譜分辨率為10nm，大小為250×191像素。剔除數(shù)據(jù)中低信噪比和低吸水率波段(1～2、105～115、150～170和223～224)，實(shí)驗(yàn)數(shù)據(jù)的實(shí)際波段數(shù)為188。實(shí)驗(yàn)所用的端元光譜庫(kù)來(lái)源與模擬實(shí)驗(yàn)一致，刪除相應(yīng)的低信噪比和低吸水率波段后，保留188個(gè)波段。實(shí)驗(yàn)參考圖4所示Tricorder軟件[22]制作的礦物分類圖，對(duì)不同算法的豐度估計(jì)結(jié)果進(jìn)行定性分析。

圖4 USGS獲得的內(nèi)華達(dá)州Cuprite礦區(qū)數(shù)據(jù)

圖5以明礬石(Alunite)、水銨長(zhǎng)石(Buddingtonite)和玉髓(Chalcedony)三種礦物為例，展示不同算法對(duì)Cuprite數(shù)據(jù)集的解混結(jié)果?？梢钥闯鯯UnSAL-TV算法得到的豐度圖均不同程度地存在過(guò)平滑現(xiàn)象，DRSU算法的豐度圖空間一致性較差(如玉髓)。所提SupLRSU算法的估計(jì)豐度圖噪聲更少(如水銨長(zhǎng)石)，更接近Tricorder參照?qǐng)D，表明該算法能有效處理真實(shí)數(shù)據(jù)。

圖5 不同算法對(duì)Cuprite數(shù)據(jù)集解混得到的豐度圖

4 結(jié)論

本文在空間先驗(yàn)稀疏解混方法的基礎(chǔ)上，一方面利用改進(jìn)的SLIC算法對(duì)高光譜圖像進(jìn)行分割，生成像元光譜高度相似的不規(guī)則同質(zhì)區(qū)域，對(duì)同質(zhì)區(qū)域的豐度向量組施加低秩約束，在促進(jìn)圖像空間一致性的同時(shí)，保護(hù)圖像的空間異質(zhì)性結(jié)構(gòu)。另一方面，在l1范數(shù)正則項(xiàng)中引入光譜加權(quán)因子，提高豐度矩陣的行稀疏性，增強(qiáng)算法從光譜庫(kù)中識(shí)別實(shí)際端元的能力。模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)實(shí)驗(yàn)證實(shí)，與SUnSAL，CLSUnSAL，DRSU，ADSpLRU，JSpBLRU算法相比，在大多數(shù)情況下，本文提出的SupLRSU算法在豐度估計(jì)方面能夠獲得更高的精度，特別是在低信噪比情況下，具有明顯優(yōu)勢(shì)，同時(shí)能更準(zhǔn)確地從光譜庫(kù)中識(shí)別端元。