侯增福， 劉镕源， 閆柏琨， 譚琨

(1.中國礦業(yè)大學(xué)國土環(huán)境與災(zāi)害監(jiān)測國家測繪地理信息局重點實驗室，徐州 221116； 2.中國自然資源航空物探遙感中心，北京 100083)

0 引言

高光譜遙感影像不同于全色和多光譜遙感影像，具有光譜分辨率高、圖譜合一的特點，在地物目標(biāo)探測領(lǐng)域具有獨特的優(yōu)勢，可廣泛應(yīng)用于環(huán)境監(jiān)測和軍事偵察等領(lǐng)域，然而在實際應(yīng)用中研究者往往很難獲得足夠的先驗知識來表征目標(biāo)類別的統(tǒng)計信息，因此在沒有可用先驗信息輔助的情況下完成異常目標(biāo)的探測，成為了近年來高光譜遙感影像目標(biāo)探測領(lǐng)域的研究重點[1]。

在高光譜影像中，異常像元的光譜往往不同于周圍背景像元的光譜信息，這就為異常像元能被探測出來創(chuàng)造了條件。由Reed和Yu在1990年發(fā)展起來的RXD(Reed-X detector)算法[2]，通過計算待探測像元與背景的馬氏距離來完成異常探測，該算法選取整幅影像作為背景信息，故又稱為全局RX(global RX，GRX)，由于使用全圖均值和協(xié)方差來估計背景均值與協(xié)方差矩陣會影響探測精度，故對此改進(jìn)的使用局部計算代替全局計算的RX又稱之為局部RX(local RX，LRX)[3-5]。然而在真實高光譜影像中，背景信息復(fù)雜，使用估計的協(xié)方差與均值向量來表示背景信息并不準(zhǔn)確?；诖颂岢龅囊恍└倪M(jìn)算法，如權(quán)重RXD(weighted-RXD，W-RXD)算法[6]和基于線性濾波的RXD(linear filter-based RXD，LF-RXD)算法[6]，這2種算法均旨在通過提高背景信息的估計來提高影像中異常被探測出的概率。一些基于核理論的探測算法，如較為經(jīng)典的非線性核心RX探測(Kernel RX)算法[7]，相比于傳統(tǒng)的探測算法[8]在異常探測中獲得了較好的探測效果。

近年來，基于信號稀疏表示的算法也被應(yīng)用于高光譜圖像目標(biāo)探測問題上[9]。然而這種算法僅僅考慮了影像的光譜信息，并沒有顧及空間信息，故將其應(yīng)用于異常探測中，往往難以取得令人滿意的效果[1]。一種基于協(xié)同表示的異常探測算法[10](collaborative-representation-based detector，CRD)認(rèn)為每一個背景像元都可以被其空間臨域像元近似表示，而異常像元則不能，并在應(yīng)用中取得了不錯的探測效果。不同于信號稀疏表示的算法，趙銳等[11]通過在異常探測器的背景信息構(gòu)建中引入魯棒性分析方法，提出了一種在核特征空間中具有魯棒性的異常探測算法； 張樂飛等[12]基于張量數(shù)據(jù)模型和張量代數(shù)運算，針對遙感數(shù)據(jù)多維或高維的特點提出了一種基于張量學(xué)習(xí)機(jī)的遙感影像目標(biāo)探測算法； 彭波等[13]基于Cholesky分解，將高維矩陣的求逆運算轉(zhuǎn)換為求解下三角線性系統(tǒng)，提出了基于Cholesky分解的逐像元實時高光譜異常探測算法。

目前，一些關(guān)于低秩分解的算法也被應(yīng)用于高光譜異常探測中，如較為經(jīng)典的魯棒性主成分分析(robust principal component analysis，RPCA)算法[14]被應(yīng)用于高光譜圖像的異常探測中[15]，其中影像部分僅僅為單子空間表示，并沒有考慮到高光譜影像中較為復(fù)雜的背景地物。針對該情況提出的低秩表示(low-rank representation，LRR)模型[16]，將低秩矩陣表示為多個子空間的線性組合。然而這種算法在使用時將自身作為字典，對應(yīng)不同高光譜影像，最優(yōu)參數(shù)往往不同，這是一個非常明顯的缺陷。Xu等[9]首次將LRR模型引入到高光譜影像的異常探測中，提出了基于低秩和稀疏表示(low-rank and sparse representation，LRaSR)的異常探測算法。另一種基于低秩表示與學(xué)習(xí)字典(low-rank representation and learned dictionary，LRRaLD)的算法[17]在LRR模型的基礎(chǔ)之上引入了僅包含背景光譜信息的學(xué)習(xí)字典，實現(xiàn)了高光譜背景與異常的有效分離，從而提高了算法的魯棒性。

然而，由于高光譜本身數(shù)據(jù)的冗余性，使用上述算法進(jìn)行異常探測時，往往需要較大的計算代價，如何在保留最大有用信息的同時，減少波段數(shù)量，從而達(dá)到減小計算代價的目的也就成為了研究的熱點。基于此產(chǎn)生的數(shù)據(jù)降維算法可概括為2類： 特征提取[18-20]和波段選擇。近些年來，波段選擇算法廣泛應(yīng)用于遙感影像的分類研究中，并取得了不錯的分類效果，如： 基于聚類分析的自組織特征映射神經(jīng)網(wǎng)絡(luò)(self-organizing feature map，SOM)[21-22]、流形學(xué)習(xí)應(yīng)用于高光譜遙感影像[23]和最佳分形波段選擇模型[24]。目前較為流行的蟻群優(yōu)化算法也已經(jīng)被應(yīng)用于高光譜圖像的降維中[25-28]。

在考慮到計算的復(fù)雜度和時間效率等綜合因素后，本文引入了一種基于波段相似尺度的線性預(yù)測(linear prediction and band similarity metric，LPaBS)算法[29]，對原始影像進(jìn)行預(yù)處理，即在原始波段特征空間進(jìn)行選擇，找到波段差異性最大的波段，從而形成原空間的一個子集，在最大程度上保留了波段的原始信息,同時降低了維度； 然后對選擇的數(shù)據(jù)子集進(jìn)行低秩表示與字典學(xué)習(xí)，并使用傳統(tǒng)經(jīng)典RXD算法進(jìn)行異常探測，旨在減少計算代價的同時提高探測精度，較好地實現(xiàn)高光譜影像的異常探測。

1 算法研究

1.1 波段選擇

高光譜遙感影像所具有的大量光譜波段為更加精細(xì)的地物分類與異常探測提供了極其豐富的信息，隨著波段數(shù)的增多，其光譜特征組合方式更是以指數(shù)形式增長，導(dǎo)致了信息的冗余和數(shù)據(jù)處理復(fù)雜性的提高。分類器和探測器的性能在很大程度上依賴于數(shù)據(jù)降維的特征提取結(jié)果，依賴于這些特征是否能夠精確地描述對象的特征[30]。本文所引入的LPaBS算法，通過在原始波段特征空間進(jìn)行選擇，找到波段差異性最大的波段，從而形成原空間的一個子集。

為了在高光譜影像中選出最具代表性的波段子集，需要某種尺度來衡量波段間的相似程度，常用方法有JM距離和空間相關(guān)性等，本文基于LPaBS算法選出差異性最大的波段A1和A2作為初始子集?，并通過?線性表示出其他波段，再繼續(xù)通過基于LPaBS算法不斷更新波段子集?，直到達(dá)到所要求的波段數(shù)。該方法中有2個初始參數(shù)，分別為初始波段對?和波段數(shù)。假設(shè)原始影像集共用N個波段，基于LPaBS法則[29]為： ①尋找初始波段對B1和B2，初始子集?={B1,B2}； ②根據(jù)評判標(biāo)準(zhǔn)選出波段B3，并通過?=?∪Bi升級子集?； ③繼續(xù)步驟②，當(dāng)子集?滿足要求時停止。初始波段對的選擇算法流程為： ①隨機(jī)選出波段A1，并將其余N-1個波段投影到A1的正交子空間上，找出最大投影波段記為A2； ②將其他的N-1個波段投影到A2的正交子空間上，找到該投影空間上最大投影波段并記為A3； ③若A3=A1，則認(rèn)為A1與A2是包含最多信息的波段對，停止循環(huán)，并將其作為初始波段對?，否則進(jìn)入步驟④； ④對于波段Ai，繼續(xù)步驟②和③，直到Ai-1=Ai+1，則將Ai-1與Ai作為初始波段對。關(guān)于評判波段相似性標(biāo)準(zhǔn)，假設(shè)當(dāng)前波段子集?={B1,B2}，通過波段B1與B2線性預(yù)測波段B′，即

B′=a0+a1B1+a2B2，

(1)

a=(a0,a1,a2)T=(XTX)-1XTy，

(2)

式中：a為參數(shù)向量，可以最小化線性預(yù)測誤差；X為一個N×3維矩陣，第1列元素全為1，第2列包含波段B1中所有像元，第3列包含波段B2中所有像元；y為一個包含波段B中所有像元的N×1維列向量； 波段B′為對波段B的線性預(yù)測值，通過最小二乘算法求解，即

(3)

式中E為預(yù)測波段誤差。通過計算所有波段與B′之間的E，并找出最大誤差對應(yīng)的波段B3，B3即為所求波段。

1.2 低秩與學(xué)習(xí)字典

高光譜數(shù)據(jù)存在一個低維線性子空間，通過尋找該空間來實現(xiàn)對高光譜數(shù)據(jù)的降維處理，如經(jīng)典的主成分分析(principal component analysis，PCA)算法，然而當(dāng)高光譜數(shù)據(jù)中存在較大噪聲或異常時，則不能取得理想效果[31]。有學(xué)者提出了RPCA[14]的算法，基于影像矩陣源于一個子空間的假設(shè)，將圖像數(shù)據(jù)分解為低秩部分與稀疏部分[16]，即

Y=L+S，

(4)

式中：Y，L，S∈Rb×p，p為像元個數(shù)，b為波段數(shù)；L表示低秩矩陣；S表示稀疏矩陣。

不同于RPCA，假設(shè)LRR是基于高光譜數(shù)據(jù)矩陣由多個子空間構(gòu)成，即

Y=DZ+S，

(5)

式中：D∈Rb×m表示字典，m為字典原子個數(shù)；Z∈

Rm×p為系數(shù)矩陣；S為包含異常值的稀疏矩陣。由于rank(DZ)≤rank(Z)故公式(5)的求解等價于

(6)

(7)

(8)

公式(8)的增廣拉格朗日函數(shù)為[16-17]

(9)

在解決問題的過程中[16-17]，字典起到了很關(guān)鍵的作用，初始字典選取的好壞決定了字典收斂程度與收斂速度，在以往的研究中有些學(xué)者提出了使用數(shù)據(jù)本身作為字典的算法，在這種情況下平衡參數(shù)β起到了決定性的作用。若參數(shù)過小，探測率不高； 若參數(shù)過大，虛警率提高，針對這些問題有學(xué)者提出使用學(xué)習(xí)字典算法，這樣在很好地解決平衡參數(shù)問題的同時提高了探測率，字典學(xué)習(xí)[17]的過程如下：

1)輸入： 數(shù)據(jù)矩陣Y和字典原子個數(shù)m。

2)初始化：m=200，γ=0.01，μ=10，ε=10-6，D為歸一化隨機(jī)正值。

步驟1： 從高光譜影像中隨機(jī)選擇m個像元。

步驟2： 進(jìn)行稀疏編碼，其公式為

(10)

步驟3： 升級字典，其公式為

(11)

步驟4： 字典D歸一化處理。

步驟5：μ→0.998μ。

3)輸出： 學(xué)習(xí)字典D。

1.3 算法流程

提出的異常探測算法流程如圖1所示，其主要過程為： ①使用線性預(yù)測法則對高光譜遙感影像進(jìn)行波段選擇，獲得最終的高光譜波段子集； ②將高光譜波段子集轉(zhuǎn)換為二維圖像數(shù)據(jù)； ③利用字典學(xué)習(xí)過程進(jìn)行字典學(xué)習(xí)，獲得D； ④利用公式(9)將二維數(shù)據(jù)矩陣分解為L和E； ⑤利用RXD算法對E進(jìn)行異常探測，獲得最終探測結(jié)果。

圖1 算法流程Fig.1 Framework of the proposed method

2 實驗結(jié)果

為了對該算法進(jìn)行驗證，使用了4幅高光譜影像進(jìn)行驗證，其中1幅基于HyMap數(shù)據(jù)的模擬數(shù)據(jù)和3幅Hyperion，HYDICE，Hyspex真實高光譜遙感影像數(shù)據(jù)。

2.1 模擬數(shù)據(jù)

為了更好地驗證本文算法效果，首先使用了1幅基于HyMap機(jī)載高光譜成像儀的模擬數(shù)據(jù)，該影像為2006年6月拍攝于美國馬薩諸塞州區(qū)域，影像大小為280像元×800像元，如圖2(a)。

(a) B14(R)，B8(G)，B1(B)真彩色合成影像

(b) 實驗數(shù)據(jù)區(qū)域(c) 真實異常地物

圖2HyMap數(shù)據(jù)集

Fig.2HyMapdataset

該影像含有126個波段，去除水汽吸收波段后余121個波段，實驗中截取影像左側(cè)中間大小150像元×150像元區(qū)域合成模擬數(shù)據(jù)，如圖2(b)，選擇紅色棉布光譜作為異常光譜。基于線性混合模型，使用異常點埋入的方法生成模擬數(shù)據(jù)，其表達(dá)式為

z=ft+(1-f)b，

(12)

式中：z為合成異常數(shù)據(jù)；f為豐度分?jǐn)?shù)；t為異常光譜；b為背景光譜。采用埋點的方法隨機(jī)生成25個異常值，其中豐富分?jǐn)?shù)為從0.05～1之間以等差數(shù)列形式生成25個豐度值。

為了評估本文提出的基于波段相似性尺度線性預(yù)測的低秩表示與學(xué)習(xí)字典(linear prediction and band similarity metric and low-rank representation and learned dictionary， LPaBS-LRRaLD)算法的優(yōu)越性，分別與GRX，LRX，基于馬氏距離的非監(jiān)督最近鄰規(guī)則子空間(unsupervised nearest regularized subspace with Mahalanobis distance，UNRS-MD)[32]與LRRaLD等算法進(jìn)行對比分析。

首先通過基于HyMap數(shù)據(jù)合成的模擬數(shù)據(jù)來證明提出算法的可行性，模擬數(shù)據(jù)為121個波段，通過LPaBS 算法選出具有代表性的80個波段。各種算法異常探測結(jié)果及對應(yīng)的接收者操作特性曲線(receiver operating characteristic，ROC)如圖3所示，表1列出了每種算法的ROC曲線下面積(area under ROC curve，AUC)與運行時間2個定量評價指標(biāo)。

(a) GRX(b) LRX(c) UNRS-MD

(d) LRRaLD(e) LPaBS-LRRaLD(f) ROC曲線

圖3HyMap模擬數(shù)據(jù)探測結(jié)果及ROC曲線
Fig.3DetectionresultsandROCcurvesofHyMapsimulationdataset

表1 HyMap模擬數(shù)據(jù)AUC與耗時性比較Tab.1 Comparison of AUC and execution time using HyMap simulation data

由圖3與表1可以看出，使用LRRaLD和LPaBS-LRRaLD算法能夠獲得較好的探測效果，但是LPaBS-LRRaLD算法在探測精度與運行時間上都優(yōu)于LRRaLD算法。

2.2 真實數(shù)據(jù)

2.2.1 HYDICE 數(shù)據(jù)

Urban數(shù)據(jù)是由HYDICE機(jī)載高光譜成像儀于城市上空拍攝而得到的空間分辨率近1 m的高光譜遙感影像，整幅影像大小為307像元×307像元，包含210個波段，去除低信噪比與水汽吸收波段后剩余174個波段，截取整幅影像右上角80像元×100像元的子塊與其對應(yīng)的真實異常地物如圖4所示。

(a) B49(R),B36(G),B18(B)真彩色合成影像(b) 實驗數(shù)據(jù)區(qū)域(c) 真實異常地物

圖4HYDICE數(shù)據(jù)集

Fig.4HYDICEdataset

為了利用HYDICE高光譜成像儀獲取的Urban數(shù)據(jù)證明本文提出算法的可行性，首先使用 LPaBS 算法選出具有代表性的100個波段，再采用LRRaLD算法進(jìn)行異常探測。各種不同算法的異常探測結(jié)果及其對應(yīng)的ROC曲線如圖5所示，表2列出了各種算法的AUC和運行時間。

(a) GRX(b) LRX(c) UNRS-MD

(d) LRRaLD(e) LPaBS-LRRaLD(f) ROC曲線

圖5HYDICE數(shù)據(jù)集探測結(jié)果及ROC曲線
Fig.5DetectionresultsandROCcurvesofHYDICEdataset

表2 HYDICE數(shù)據(jù)集AUC與耗時性比較Tab.2 Comparison of AUC and execution time using HYDICE data set

通過圖5可以看出，在各種算法中，使用LRRaLD和LPaBS-LRRaLD算法能夠獲得較好的探測效果。由表2數(shù)據(jù)可知，GRX算法運行時間最短，但是探測精度明顯低于LRRaLD和LPaBS-LRRaLD算法，LPaBS-LRRaLD算法在探測精度與運行時間上明顯優(yōu)于其他算法。

2.2.2 Hyperion 數(shù)據(jù)

Hyperion遙感影像數(shù)據(jù)含有242個波段，光譜分辨率為10 nm ，波長范圍為357～2 576 nm。實驗使用的影像數(shù)據(jù)采集于2008年，影像中主要包括美國印第安納州的農(nóng)業(yè)區(qū)。去除低信噪比與未定標(biāo)波段后余149個波段，在實驗中使用含有真實異常數(shù)據(jù)的150像元×150像元大小的子區(qū)域(如圖6所示)完成實驗，影像中異常值主要為儲物倉庫和屋頂。同樣，為了驗證本文算法的有效性，首先基于LPaBS 算法選出具有代表性的100個波段，然后使用學(xué)習(xí)字典并利用增廣拉格朗日公式求解進(jìn)行異常探測，各種不同算法的探測結(jié)果及其對應(yīng)的ROC曲線如圖7所示，每種算法的AUC與運行時間如表3所示。

(a) B29(R)，B23(G)，B16(B)假彩色合成影像(b) 真實異常地物

圖6Hyperion數(shù)據(jù)集

Fig.6Hyperiondataset

(a) GRX(b) LRX(c) UNRS-MD

(d) LRRaLD(e) LPaBS-LRRaLD(f) ROC曲線

圖7Hyperion數(shù)據(jù)集探測結(jié)果及ROC曲線
Fig.7DetectionresultsandROCcurvesofHyperiondataset

表3 Hyperion數(shù)據(jù)集AUC與耗時性比較Tab.3 Comparison of AUC and execution time using Hyperion data set

通過圖7可以看出，使用Hyperion數(shù)據(jù)進(jìn)行異常探測，除LRX算法外，其他探測算法均具有較高的探測精度，其中LPaBS-LRRaLD算法探測精度最高。由表3數(shù)據(jù)可知，GRX算法運行時間最短，但是探測精度均低于UNRS-MD，LRRaLD和LPaBS-LRRaLD算法，其中LPaBS-LRRaLD算法在探測精度與運行時間上明顯優(yōu)于其他算法。

2.2.3 Hyspex 數(shù)據(jù)

該數(shù)據(jù)是由機(jī)載Hyspex高光譜成像儀于2014年11月在徐州泉山區(qū)附近拍攝的條帶影像，包含160個可見光波段與288個短波紅外波段，光譜范圍為415 ～2 508 nm，空間分辨率近1 m，實驗中截取條帶中340像元×260像元區(qū)域，其中異常值為彩鋼房屋，影像與真實異常地物如圖8所示。

(a) B63(R)，B38(G)，B13(B)真彩色合成影像(b) 真實異常地物

圖8Hyspex數(shù)據(jù)集

Fig.8Hyspexdataset

在Hyspex數(shù)據(jù)集中，首先基于LPaBS 算法選出具有代表性的45個波段，然后使用學(xué)習(xí)字典算法進(jìn)行異常探測。

各種探測算法實驗結(jié)果及各種算法的ROC曲線如圖9所示，每種算法的AUC與運行時間如表4所示。

(a) GRX(b) LRX(c) UNRS-MD

(d) LRRaLD(e) LPaBS-LRRaLD(f) ROC曲線

圖9Hyspex數(shù)據(jù)集探測結(jié)果及ROC曲線
Fig.9DetectionresultsandROCcurvesofHyspexdataset

表4 Hyspex數(shù)據(jù)集AUC與耗時性比較Tab.4 Comparison of AUC and execution time using Hyspex data set

通過圖9可以看出，使用Hyspex數(shù)據(jù)進(jìn)行異常探測，其中的 LPaBS-LRRaLD算法探測精度明顯高于其他算法。由表4數(shù)據(jù)可知，除GRX算法運行時間最短外，LPaBS-LRRaLD算法在探測精度與運行時間上明顯優(yōu)于其他算法。

2.3 實驗小結(jié)

通過模擬數(shù)據(jù)與真實數(shù)據(jù)實驗，可以看出本文提出算法的可行性，對照實驗中所使用的RXD算法與最近提出的UNRS-MD和LRRaLD算法都是基于原始數(shù)據(jù)進(jìn)行的異常探測，而本文提出算法在降維的同時探測精度得到了提高，尤其是Hyspex數(shù)據(jù)，波段數(shù)從448個減少到45個，很大程度上去除了信息冗余，在減少計算代價的同時探測精度也得到了提高。

3 結(jié)論

針對全局高光譜異常，將低秩分解的算法引入到高光譜異常探測中，并通過低秩分解將圖像表征為背景低秩矩陣與稀疏矩陣，在求解背景低秩矩陣過程中采用學(xué)習(xí)字典來提高背景字典的準(zhǔn)確性與魯棒性，同時顧及維數(shù)災(zāi)難對高光譜影響異常探測的影響。首先，采用基于波段相似性線性預(yù)測的算法進(jìn)行降維，在保持原有波段信息不變性的同時有效地去除數(shù)據(jù)冗余； 然后，結(jié)合學(xué)習(xí)字典算法在低秩分解過程中提高背景與異常信息可分性的同時，更好地挖掘數(shù)據(jù)本身的低秩特性，從而達(dá)到快速收斂； 最后，使用傳統(tǒng)的RXD算法對稀疏矩陣進(jìn)行異常探測。

實驗表明，本文算法與同類算法相比，在高光譜影像異常探測中，在進(jìn)一步降低計算代價的同時，提高了異常探測率，因此該算法更具有實際應(yīng)用意義。由于學(xué)習(xí)字典的隨機(jī)性，會使得背景字典中存在異常的小概率事件發(fā)生，針對這種情況，如何找到完全不存在異常的背景字典來表征背景矩陣，從而使得背景與異常更加有效地分離將是需要進(jìn)一步研究的問題； 同時也建議嘗試其他探測算法對稀疏矩陣進(jìn)行異常探測，以達(dá)到更高精度的探測率。