基于濾波張量分析的高光譜目標檢測方法*

楊帥，計璐艷，耿修瑞

(中國科學院大學， 北京 100049；中國科學院空天信息創(chuàng)新研究院 中國科學院空間信息處理與應用系統(tǒng)技術重點實驗室， 北京 100094)

高光譜遙感指具有高光譜分辨率的遙感數(shù)據(jù)獲取、處理、分析和應用的科學與技術，通常采用覆蓋一定波譜范圍的成像光譜儀和非成像光譜儀兩種傳感器獲取數(shù)據(jù)，利用大量窄波段電磁波獲取感興趣目標的理化信息[1]。隨著高光譜遙感技術的發(fā)展，現(xiàn)在的傳感器已經(jīng)可以獲取地物從可見到紅外的多個甚至數(shù)百個波段[2]。這些豐富的光譜信息使得高光譜遙感數(shù)據(jù)能夠提供更為精細的地物信息，從而提高地物的分類及目標識別的精度。現(xiàn)在高光譜手段已經(jīng)在軍事偵察、環(huán)境監(jiān)測、地質調查及精準農(nóng)業(yè)等方面得到極大應用[3]，例如遙感技術以其能快速、連續(xù)、大面積地獲取地面信息的優(yōu)勢，迅速成為環(huán)境監(jiān)測的重要手段[4]。

在眾多的圖像應用中，目標檢測一直是高光譜遙感領域研究的熱點內容。目標檢測大概有以下幾種方式：1)基于光譜的吸收特征。比如植被的歸一化植被指數(shù)[5]，歸一化差分水指數(shù)[6]等。這些方法都是基于光譜吸收特征，利用目標反射率計算，需要知道目標的光譜吸收特征，但是要求背景地物不具有這種特征。2)簡單的光譜匹配，比如基于光譜間的最小距離匹配算法[7]，光譜角度匹配法(spectral angle mapping，SAM)[8]等，光譜匹配檢測是通過與目標光譜先驗特征相關或匹配來尋找目標，感興趣目標的光譜特征可從光譜庫中或是從同一觀測場景中的已知目標像素處獲取[9]。但一方面由于“異物同譜”現(xiàn)象的存在，目標和背景的光譜可能只有細微的差別；另一方面簡單的光譜匹配方法沒有利用任何背景地物信息，難以獲得良好檢測效果。一個好的目標檢測算法還需要考慮背景和目標在特征空間中的分布結構[10]，從而提高檢測效果。3)基于單形體幾何的目標檢測，Ren和Chang構造了一系列基于信號子空間投影和最小二乘原理的匹配檢測方法，包括基于正交子空間投影(orthogonal subspace projection，OSP)的方法[11]，在此基礎上發(fā)展了推廣正交子空間投影(generalized orthogonal subspace projection，GOSP)[12]；利用信號雜波噪聲模型的干擾子空間投影(interference subspace projection，ISP)算法[13]；利用噪聲子空間的噪聲子空間投影算法(noise subspace projection，NSP)[14]等。這種方法需要把所有背景地物的端元都提出來，然后通過正交投影或者正交子空間投影壓制背景。4)基于統(tǒng)計的目標檢測，這種方法很好地利用了圖像的統(tǒng)計信息(主要是二階統(tǒng)計信息)，是最常用的一系列目標檢測算法。尋麗娜等提出約束能量最小(constrained energy minimization，CEM)檢測器[15]。Manolakis和Shan提出MF(matched filter)算法[16]，Kraut和Scharf提出自適應余弦估計(adaptive cosine estimate，ACE)算法[17]，近些年還有Geng等提出的慧眼(clever eye，CE)算法[18]。5)非線性目標檢測，如基于核函數(shù)的約束能量最小化(kernel-based constrained energy minimization，KCEM)[19]、基于核函數(shù)的匹配濾波(kernel-based matched filter，KMF)算法[20]，基于張量表示的高光譜圖像目標檢測[21]，分段線性策略(piecewise linear strategy，PLS)[22]。這類方法有很強的非線性檢測能力，但是也容易把目標像元歸成背景。

以上這些方法都是基于單時相的遙感圖像，雖然取得了一定效果并被廣泛應用，但是沒有有效利用遙感圖像時間維度信息。由于水文、植被、土壤狀況等具有季節(jié)性變化的特征，待檢測目標及周邊環(huán)境在不同季節(jié)有不同的光譜，因此多時相數(shù)據(jù)對于地物目標分類與檢測具有重要意義。如使用多時相的SPOT-4植被傳感器數(shù)據(jù)的中國東北森林的分類[23]和基于多時相的MODIS圖像的中國南部水稻田分布檢測[24]等。在多時相遙感圖像處理領域，Geng等提出基于多重線性函數(shù)的多時相目標檢測方法(filter tensor analysis，F(xiàn)TA)[25]，將不同時相的遙感圖像作為輸入來檢測目標，充分利用多時相遙感圖像中豐富的地物時間維度信息，提高檢測精度。Xi等在此基礎上提出多時相多目標濾波張量算法[26]，能夠實現(xiàn)在多時相遙感數(shù)據(jù)上同時檢測多個目標。

FTA對多時相遙感數(shù)據(jù)進行目標檢測的物理基礎在于，數(shù)據(jù)一般在同一個時相內光譜變化較小，而在不同的時相間光譜差異較大。類似地，對于高光譜數(shù)據(jù)，在光譜維通?？煞譃榭梢姽?、近紅外、短波紅外等波段范圍，在很多情況下，地物在不同譜段間的光譜差異往往大于同一譜段內的光譜差異。而現(xiàn)有的高光譜目標檢測算法往往對各個波段區(qū)間不加區(qū)分、統(tǒng)一對待，無法充分利用圖像波段的物理信息。因此為了最優(yōu)化地利用波段物理信息，本文將數(shù)據(jù)的光譜維與時間維對應，即在波段維對數(shù)據(jù)進行分段，利用多時相目標檢測算法FTA對單時相數(shù)據(jù)進行目標檢測，形成一種新的目標檢測非線性算法，稱之為分波段FTA(band-divided FTA，bd-FTA)算法。并與現(xiàn)有的目標檢測算法進行比較，以驗證本文方法的有效性。

1 數(shù)據(jù)及背景

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

為了測試分波段FTA在目標檢測中的性能，選擇兩景AVIRIS高光譜圖像：Salinas和Indian Pines數(shù)據(jù)(數(shù)據(jù)來源：http:∥www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes)。圖像的具體描述如下：

1)Salinas高光譜數(shù)據(jù)。該數(shù)據(jù)的獲取時間是1998年，地點是美國加州的薩利納斯谷，由AVIRIS高光譜傳感器獲取，原始數(shù)據(jù)包含220個波段，波段覆蓋范圍為400～2 500 nm，去除水汽吸收嚴重的第 108～112，154～167，和第 224 個波段后，還剩下204個波段。圖像大小為512像素×217像素，其空間分辨率較高，可達到3.7 m。從圖1(a)可以看到，地物空間分布相對均勻，共有16種地物覆蓋類型；其中第13種植物lettuce romaine 6wk 被選為待檢測目標。

2)Indian Pines高光譜數(shù)據(jù)。該數(shù)據(jù)由美國國家航空航天局在1992年利用AVIRIS采集于美國印第安納州，其波長覆蓋范圍為400～2 500 nm，原始數(shù)據(jù)為220個波段，剔除第104～108，150～163和第220個水吸收波段后，剩余200個有效波段。圖像大小為145像素×145像素，空間分辨率約為20 m，共包含16種地物覆蓋類別；其中第9種植物Oats 被選為待檢測目標，如圖1(b)所示。

圖1 2個數(shù)據(jù)集的假彩色圖像、各類別真值圖及目標真值圖Fig.1 The false color images, groundtruth classes， and the groundtruth images of the target for Salinas and Indian Pines data

1.2 CEM算法

CEM在保持目標輸出為常數(shù)的情況下，從能量的角度對背景進行壓制，該算法的設計思想與數(shù)字信號處理領域的線性約束最小方差波束形成器類似，其在保持特定信號輸出不變的同時，將其他可能存在的干擾信號進行衰減壓制。CEM算子及其濾波輸出結果[15]為

(1)

(2)

其中：R為數(shù)據(jù)的自相關矩陣，r為圖像像元光譜，d為目標光譜向量，wCEM為CEM算子，yCEM為輸出結果。

1.3 FTA算法

通常，光譜目標在不同的時間階段具有不同的光譜特征。因此，對于M時相遙感數(shù)據(jù)集，可以使用一個M階張量濾波器進行目標檢測，它對應于M重線性函數(shù)，在保持目標輸出值不變的同時，最小化檢測器輸出能量，這個算法被稱為FTA。

(3)

其中

(4)

(5)

(6)

其中：?表示克羅內克積。可以看到FTA算子在形式上與CEM是一致的。

其中L=L1×L2×…×LM，因此總共的時間復雜度是

O((L1×L2×…×LM)2(N+1)+

(L1×L2×…×LM)3).

可以明顯看出，F(xiàn)TA的時間復雜度隨著N和L1,L2,…,LM增加而增加，尤其是當時相數(shù)M或者相應的波段數(shù)L1,L2,…,LM較大時，計算復雜度將會變得十分巨大。

2 分波段FTA算法

在多時相目標檢測中，F(xiàn)TA能夠最優(yōu)化地利用各個時相的光譜信息，從而取得比CEM和MF更加令人滿意的結果。在單時相高光譜目標檢測中，F(xiàn)TA對時相信息的優(yōu)化策略可以被借鑒，以實現(xiàn)對波段信息的更優(yōu)利用。分波段FTA將多時相數(shù)據(jù)的時相與高光譜數(shù)據(jù)的波段相對應，將高光譜數(shù)據(jù)在波段維進行波段分區(qū)，實現(xiàn)對單時相高光譜數(shù)據(jù)的目標檢測。算法流程圖如圖2所示，主要包括以下3個步驟：

1)波段分區(qū)：將單時相高光譜圖像劃分為幾個波段區(qū)間，每個波段范圍對應多時相數(shù)據(jù)的每個時相，相當于將單時相高光譜圖像變換為多時相高光譜圖像，即將波段信息轉換為時相信息，以充分利用波段的物理信息。

圖2 分波段FTA 算法流程圖Fig.2 Flow chart of band-divided FTA

2)數(shù)據(jù)降維：對于高光譜數(shù)據(jù)，由于波段數(shù)多，即使分段后，各個波段區(qū)間的波段數(shù)仍然較多，而FTA算法計算及空間復雜度較高，因此需對各個波段區(qū)間采用降維處理，本文選擇性能較優(yōu)的最小噪聲分離(minimum noise fraction，MNF)法[27]，獲得有限量波段。

3)FTA目標檢測：將分波段降維后的數(shù)據(jù)利用FTA目標檢測算法，實現(xiàn)對單時相高光譜圖像目標進行檢測。

下面對各個部分進行分別介紹。

2.1 波段區(qū)間劃分

一方面，高光譜傳感器在不同波段范圍的成像機理往往不同，比如可見光(400～750 nm)、近紅外(750～1 000 nm)、短波紅外(1 000～2 500 nm)等。另一方面，在很多情況下，地物在不同譜段間的光譜差異往往大于同一譜段內的光譜差異。本文實驗將數(shù)據(jù)分為2個波段區(qū)間。首先從傳感器角度而言，AVIRIS可見光波段探測器類型為Si探測器，近紅外波段為InGaAr探測器，而短波紅外探測器為InSb探測器。所以本文從傳感器角度而言，首先將高光譜數(shù)據(jù)分為可見光、近紅外及短波紅外波段3個區(qū)間，這種波段區(qū)間劃分稱為分段方法1。

從另一個角度，光譜特征曲線而言，本文以Salinas數(shù)據(jù)為例，如圖3所示為Salinas高光譜圖像隨機獲取100個像素點的光譜特征曲線。從中能夠看出，由于水吸收影響，光譜特征曲線在820，950，1 150，1 400和1 900 nm附近反射率出現(xiàn)了谷值，而可見光部分則無該情況出現(xiàn)。因此，從光譜特征曲線角度，本文將高光譜數(shù)據(jù)分為可見光波段(400～750 nm)，近紅外和短波紅外波段(750～2 500 nm)2個區(qū)間，這種波段區(qū)間劃分稱為分段方法2。

2.2 分波段降維

FTA算法的時間和空間復雜度均很高，特別是當數(shù)據(jù)波段數(shù)較大時，對計算機的內存要求會大大增加，且計算時間將遠遠無法滿足實際應用需求。通常，在高光譜圖像目標檢測中如果計算復雜度較高，則需要進行降維處理。本文采用MNF，該方法可以判定圖像數(shù)據(jù)內在的維數(shù)(即波段數(shù))，分離數(shù)據(jù)中的噪聲，以減少后續(xù)處理中的計算需求量。MNF本質上是兩次層疊的主成分變換。第1次變換(基于估計的噪聲協(xié)方差矩陣)用于分離和重新調節(jié)數(shù)據(jù)中的噪聲，該操作使變換后的噪聲數(shù)據(jù)只有最小的方差且沒有波段間的相關。第2步是對噪聲白化數(shù)據(jù)的標準主成分變換。為了進一步進行波譜處理，通過檢查最終特征值和相關圖像來判定數(shù)據(jù)的內在維數(shù)。數(shù)據(jù)空間可分為兩部分：一部分與較大特征值和相對應的特征圖像相關，其余部分與近似相同的特征值以及噪聲占主導地位的圖像相關。具體地，

1)噪聲估計，假設數(shù)據(jù)中噪聲的協(xié)方差矩陣CN，將其對角化為矩陣DN，即

DN=UTCNU,

(7)

式中：DN為CN的特征值按照降序排列的對角矩

圖3 Salinas數(shù)據(jù)隨機選取100個像素光譜曲線Fig.3 Spectral characteristic curves of 100 randomly selected pixels in Salinas dataset

陣；U為由特征向量組成的正交矩陣，進一步變換得

I=OTCNO,

(8)

(9)

其中：I為單位矩陣，O為變換矩陣。

2)對噪聲數(shù)據(jù)進行標準主成分變換，即

G=OTCXO,

(10)

式中：CX為數(shù)據(jù)X的協(xié)方差矩陣；G為經(jīng)過O變換后的矩陣，進一步將其對角化為矩陣F

F=VTGV,

(11)

式中：F為G的特征值按照降序排列的對角矩陣；因此，對數(shù)據(jù)X可通過以下表達式進行MNF變換

Y=OVX.

(12)

2.3 分波段FTA

(13)

此處，運算符×j(j=1,2,…,P)代表模j乘法運算符。式(13)中的張量向量積W×1y(1)×2y(2)×…×Py(P)是個標量，可以用下式表示：

W×1yi(1)×2yi(2)×…×Pyi(P)=

(yi(P)?yi(P-1)?…?yi(1))Tvec(W).

(14)

其中:?表示克羅內克積,vec()是一個將L1×L2×…×LP維的張量轉化為L1×L2×…×LP維的向量的運算符，同理可得

W×1d(1)×2d(2)×…×Pd(P)=

(d(P)?d(P-1)?…?d(1))Tvec(W).

(15)

(16)

可得分波段FTA濾波算子的表達式

(17)

類似FTA，可以得到分波段FTA的計算復雜度為O((L1×L2×…×LP)2(N+1)+(L1×L2×…×LP)3)，其中L1,L2,…,LP為經(jīng)過MNF降維后的波段數(shù)，如實驗中經(jīng)過MNF降維后，波段數(shù)由幾十個降低至幾個，復雜度大大降低。FTA復雜度較高，很難直接應用于高光譜數(shù)據(jù)。本文由于對各個波段區(qū)間使用了MNF變換，大大降低了各個波段區(qū)間的波段數(shù)，因此相比于原始FTA，分段FTA的計算和空間復雜度均得到極大的降低，可應用于高光譜目標檢測。

3 實驗與分析

3.1 實驗技術指標

在本實驗中，將本文所提出的分波段FTA算法檢測結果與前文所提的CEM，KCEM，KMF，CE算法相比較。為了評估不同方法的目標檢測結果，選用接收者操作特性曲線 (receiver operating characteristic curve，ROC)[28]描述目標檢測算法的性能。ROC曲線是指在特定刺激條件下，以被試在不同判斷標準下所得的虛報概率為橫坐標，以擊中概率為縱坐標，畫得的各點的連線。通常來講，如果某條ROC曲線一直在其他ROC曲線的左上方，表明該曲線對應的檢測算法的性能是最優(yōu)的。但是對于性能相近的方法而言，其ROC曲線往往會交疊在一起，僅僅從視覺的角度上，很難分辨其優(yōu)劣性，為此需要一些具體的數(shù)值指標對分類器的性能進行輔助判定。一種較為常用的性能指標為曲線下面積(area under curve，AUC)表示ROC曲線下各部分的面積求和。AUC越大表明檢測效果越好[29]。

3.2 Salinas數(shù)據(jù)實驗

在Salinas高光譜數(shù)據(jù)實驗中，第13種植物lettuce romaine 6wk 被選為待檢測目標。針對2種不同的分段方法進行實驗。

3.2.1 分段方法1

首先204個波段被分為3個波段區(qū)間——可見光(400～750 nm)、近紅外(400～1 000 nm)及短波紅外波段(1 000～2 500 nm)，可見光有波段1～43共43個，近紅外有波段44～68共25個，短波紅外有波段69～204共136個。每個波段區(qū)間用MNF降維后得到6個波段，利用分波段FTA處理得到檢測結果。為保證波段數(shù)一致，204個波段用MNF降維獲取18個波段，利用CEM、KCEM、KMF、CE處理得到檢測結果。5種算法的結果和ROC曲線如圖4所示，從目視判讀角度可以看到，分波段FTA對于lettuce romaine 6wk的檢測更加準確，對于背景噪聲的抑制均優(yōu)于其他算法，特別是目標鄰近區(qū)域。5種檢測算法的AUC結果如表1所示，分波段FTA的AUC值(加粗，下同)為最大，效果最優(yōu)。

圖4 Salinas實驗結果圖及ROC曲線圖(分段方法1)Fig.4 The experiment results of Salinas and ROC curve (band division method one)

3.2.2 分段方法2

類似地，根據(jù)分段方法2，將數(shù)據(jù)分成可見光波段(400～750 nm)，近紅外和短波紅外波段(750～2 500 nm)2個波段區(qū)間，然后利用MNF將2個波段區(qū)間均降為6個波段，輸入CEM、KCEM、KMF、CE，分波段FTA進行目標檢測，如圖5、表1所示，可以看到，與分段方法1結果類似，分波段FTA性能最優(yōu)。

圖5 Salinas實驗結果圖及ROC曲線圖(分段方法2)Fig.5 The experiment results of Salinas and ROC curve (band division method two)

表1 Salinas實驗2種分段方法各檢測算法AUC值Table 1 AUC of each detection algorithm in Salinas experiment with two band division methods

3.3 Indian Pines數(shù)據(jù)實驗

在Indian Pines高光譜數(shù)據(jù)實驗中，第9種植物Oats 被選為待檢測目標。與3.2節(jié)相同，針對2種不同的分段方法進行實驗。

3.3.1 分段方法1

5種算法的結果和ROC曲線圖如圖6所示，從目視判讀角度可以看到，分波段FTA對于Oats的檢測更加準確，圖中藍色框目標區(qū)域的亮度更高。5種檢測算法的AUC結果如表2所示，分波段FTA的AUC值為最大，效果最優(yōu)。

3.3.2 分段方法2

5種算法的結果和ROC曲線圖如圖7所示，5種檢測算法的AUC結果如表2所示，結果與分段方法1類似，分波段FTA藍色框目標區(qū)域更加明顯，對于周圍背景的抑制效果更好，AUC值再次證明了分波段FTA的有效性。

表2 Indian Pines實驗2種分段方法各檢測算法AUC值Table 2 AUC results of each detection algorithm in Indian Pines experiment with two band division methods

4 總結

隨著遙感技術的快速發(fā)展，遙感圖像的獲取變得極為方便，遙感圖像的質量也愈來愈好，空間、時間以及光譜分辨率均得到很大的提升。但是作為遙感圖像處理領域的一個很重要的分支，常用的遙感圖像目標檢測算法對于波段的物理信息利用不夠充分。本研究將高光譜圖像的波段維信息與時間維信息相對應，對最新的多時相目標檢測算法FTA進行改進，提出分波段FTA算法，實現(xiàn)對高光譜數(shù)據(jù)波段信息的最佳利用。經(jīng)Salinas和Indian Pines高光譜圖像實驗驗證，在波段區(qū)間數(shù)為2和3的2種不同波段區(qū)間劃分方法中，分波段FTA目標檢測算法的 AUC值均為最大，性能優(yōu)于CEM、KCEM、KMF和CE。其中，波段區(qū)間數(shù)為2的AUC值優(yōu)于波段區(qū)間數(shù)為3的AUC值，波段區(qū)間數(shù)增加時似乎分波段FTA算法的性能下降。

在本次試驗中，基于FTA的高光譜圖像目標檢測算法相對于算法CEM及KCEM等能夠獲得更優(yōu)效果，算法仍可以從以下方面改進：1)波段區(qū)間的劃分對于基于FTA的高光譜圖像目標檢測性能存在一定的影響，需進一步研究分波段區(qū)間數(shù)目對檢測結果的影響，分波段依據(jù)何種準則對波段進行分類可取得最優(yōu)效果，以充分利用目標波段的物理信息。2)由于FTA算法復雜度較高，而高光譜圖像波段數(shù)多，因此針對單時相高光譜圖像的目標檢測需要數(shù)據(jù)降維處理，本實驗以特征提取中的MNF為例，其他數(shù)據(jù)降維方法如波段選擇對分波段FTA算法的影響有待研究。下一步的工作目標是從上述兩個角度進一步優(yōu)化算法，以期提高目標檢測算法的準確性和穩(wěn)定性。