高光譜影像子空間分析孤立森林異常目標探測方法

黃遠程，薛園園，李朋飛

西安科技大學測繪科學與技術學院，西安 710054

遙感中的異常信息通常以低概率的形式出現在地表環境中，且與周圍環境存在光譜或空間屬性上明顯差異[1-2]。異常目標檢測被認為是一個二分類問題，它將待測像素分為異常類或背景類[3]，它的特點是在背景和異常目標無先驗知識的情況下識別出感興趣目標[4-5]。而高光譜圖像具有數百個光譜波段，為每個像素提供了豐富的光譜信息，這一優勢使其成為異常目標檢測最有價值的遙感手段[6-7]，比如利用高光譜影像檢測海洋環境背景下的船舶溢油[8]和機場的飛機，發現地表的礦物[9]等。

異常目標檢測有兩類代表性的方法：①基于統計模型的異常檢測，這類方法假設背景服從特定的概率分布，而背離該概率分布的樣本則為異常。其中，應用最為廣泛的是RX算子[10-11]，它建立在背景服從高斯分布假設的基礎上，通過計算待探測像素與背景的馬氏距離來判斷像素點是否屬于異常目標。但由于背景并不一定簡單地服從單一的高斯分布，故而產生較高的虛警率。局部RX算子則在背景的協方差與均值估計時容易受潛在的異常信號影響。②基于誤差重建的異常檢測，這類方法不對背景的概率分布做一般假設，而是采用構建背景字典來對待檢測像元進行重構[12-13]，以重構誤差來反映異常。但該方法多以線性稀疏表達為主，無法有效地處理線性不可分問題；此外字典構建是關鍵，全局字典依賴于非監督學習算法獲得，局部字典則由局部鄰域像素構成，局部構建的字典也易受可能的異常像元影響[14]。

與上述方法不同，孤立森林(isolation forest,iForest)異常檢測定義異常樣本為分布稀疏且易被孤立的離群點[15]，如圖1(b)點x0。iForest在子采樣的基礎上使用隨機超平面切割特征空間來構建孤立樹(isolation tree，iTree)，如圖1(a)中正常的樣本點xi需多次切割到達更深的樹葉節點上才能被孤立；圖1(b)中異常點x0只需較少的切割次數就能夠被孤立，它靠近孤立樹的根部，有較短的路徑長度。t個iTrees組成iForest，然后用生成的iForest來評估測試數據，iForest中具有短的平均路徑長度的點被識別為異常，而這個異常度值是通過將森林里的所有樹用集成的方法計算得到的。它為解決異常檢測問題提供了一個新的視角，即不需要對背景概率分布做出假設，且具有集成多種不同特征的優勢。但孤立樹的構建將決定著iForest模型異常檢測的性能。特別的，若將iForest直接應用于高維高光譜遙感影像進行異常檢測，則容易出現兩個主要的問題：①遙感影像中地物復雜，對輸入特征進行隨機閾值剖分時不易有效的區分背景與目標；②仍有重要的光譜特征易在建完樹后沒被使用[16-17]。這是因為高光譜影像波段眾多且相鄰波段之間相關性大，iTrees的每一個結點對應的特征均采用隨機的方式定義，無法有效地保證具有區別能力的關鍵特征被選中。

圖1 孤立森林對特征空間中數據點切割Fig.1 Schematic diagram of data point cutting in feature space

為了利用iForest求解異常目標問題的優勢，同時克服其直接應用于高光譜影像分析的困難，本文提出了一種子空間分析孤立森林異常目標探測方法，其總體流程如圖2所示。針對問題①，采用正交子空間分析[18]技術增強輸入特征影像中潛在異常目標與背景之間的對比度，以減少背景的干擾；針對問題②，對增強的目標特征影像利用主成分分析方法將其投影到一個低維子空間進行降維[19]，從而降低維數過高給iForest算法帶來的不確定性；在此基礎上，采用在整幅圖像中隨機采樣的方式來構建孤立森林獲得初步檢測結果，稱之為全局iForest異常檢測；最后利用滑動窗口進行局部的iForest異常檢測來細化初始異常檢測結果，得到最終的異常分布。

圖2 基于子空間分析的iForest異常目標檢測流程Fig.2 iForest anomaly detection flowchart with subspace analysis and global-local processing

1 本文方法

1.1 正交子空間抑制背景類信息

假設r∈RL×1為L維高光譜圖像的任意像元，根據正交子空間分析理論，假設目標子空間為S∈RL×c，背景子空間U∈RL×k，它們構成了影像的特征空間，則圖像的每個像元可由目標子空間S和背景子空間U進行線性表達，其中c和k分別是目標子空間和背景子空間的大小。如式(1)所示

r=Sα+Uγ+ε

(1)

(2)

然后將其乘以每一個像元矢量如式(3)

(3)

即像元向量投影于背景的正交子空間中，這樣消除了式(1)的第2項，使得與U相關的背景信息被抑制后僅剩下目標相關的信息。這一思想稱為高光譜影像正交子空間投影(OSP)的分析[20]。

本文采用兩種方法定義式(1)的U，分別為主成分特征子空間法和典型類判別特征子空間法。這兩類方法定義的背景子空間和目標相關子空間符合式(1)的線性模型。

(1) 主成分特征子空間法：假設[U1,U2,…,UL]為高光譜影像的主成分特征向量，其對應的特征值由大到小排列[21]。由于異常目標在影像中分布稀少，其信息主要表現在較小特征值對應的特征中，故可定義背景子空間為最大的k個特征值對應的特征向量組成，即U=[U1,U2,…,Uk]。

(2) 典型類判別特征子空間法：首先對高光譜圖像進行k-means聚類[22]，統計聚類后的每一類像元數的占比，并定義預估異常像元占比δ，將像元數量占比大于這個值的類視為正常類或稱背景類。接著采用Fisher判別分析[23]，最小化類內散度和最大化類間散度，對所有的k個背景類樣本進行判別分析計算，得到的判別特征向量組成背景子空間U=[U1,U2,…,Uk-1]。

1.2 背景正交子空間增強變換的孤立森林異常檢測

結合正交子空間分析的特點和經典的iForest算法提出如下計算方法。首先采用主成分特征子空間法或典型類判別特征子空間法估計得到的背景子空間對輸入圖像進行正交子空間投影，從而抑制背景，增強目標，使主要的信息與異常相關；為降低因高光譜影像維數過高給算法帶來的風險，接著對輸入的特征繼續采用主成分分析進行降維處理，使主要的異常信息集中在幾個主成分中。在此基礎上利用iForest對主成分特征影像進行訓練和對待檢測樣本進行測試。具體的訓練和測試方法如下。

(1) 訓練階段是建立孤立森林。假設孤立樹iTree的數目為t，由t棵iTree構成孤立森林iForest。iTree的建立是通過對訓練集的遞歸分割來建立的，遞歸分割結束的標志是所有的樣本被孤立，或者樹達到了指定的高度。假設樹的高度限制為l，l與子樣本數量n的關系為l=ceiling(log2(n))，近似等于樹的平均高度。以下為樹的構造過程。

Step1：從被背景子空間抑制后的高光譜影像中隨機選擇n個樣本作為子樣本集，X={x1,x2,…,xn}，樣本的特征維數d是主成分分析降維后的特征維數。

Step2：為了構建一棵孤立樹，需要從d個特征中隨機選擇一個特征q及其分割值p，p值介于該特征的最小值和最大值之間，與經典的iForest不同在于p對輸入特征的隨機剖分能更容易區分背景與目標像元。

Step3：假設T是孤立樹的一個節點，它要么是沒有子節點的葉子節點，要么是只有兩個子節點(Tl,Tr)的內部節點，每一步分割，都包含特征q和分割值p，對特征q的任一樣本qi，若其滿足qi

Step4：子節點遞歸的重復Step2和Step3，直到樹的高度達到了限制值l或節點上的樣本都相同。

(2) 測試階段計算待檢測樣本的異常值。每一個測試樣本x(x∈Rd×1)的異常值由期望路徑長度E(h(x))表達得到，樹的路徑長度h(x)表示樣本x從iTree的根節點到葉節點過程中經過的邊的數目。為了便于評價異常值，對E(h(x))采用歸一化處理(4)來定義最終的異常分值

S(x)=2-E(h(x))/C(n)

(4)

式中，C(n)是n個樣本構建二叉樹的平均路徑長度，定義如式(5)

C(n)=2H(n-1)-(2(n-1)/n)

(5)

式中，H(·)為調和函數，用ln(·)+0.577 215 664 9(歐拉常數)估算。對高光譜影像中的每一個像元均進行上述步驟，得到異常檢測結果圖D∈Rh×w，h和w分別為影像的行數和列數。

異常值0

1.3 局部化iForest細化處理

影像中亮度值比較大的像元具有離群的特點，易被誤檢為異常。這類地物與感興趣的異常目標不同在于它們往往成面狀或線狀。為了消除其對異常的影響，采用二次iForest計算，稱之為局部iForest細化處理，其實質是在局部范圍內再進行一次iForest異常檢測，并將1.2節的處理結果稱為初始異常檢測結果。具體步驟如下。

(1) 二值化：將初始的異常檢測結果圖D轉換為二值化圖像B∈Rh×w

(6)

式中，閾值s由Otsu[24]算法確定；Bi和Di分別為B和D對應的像元。

(2) 基于分塊分析的iForest異常檢測：對二值圖像進行形狀分析，找出二值圖像B的所有連通分量，并計算連通分量內的像元數。接著對圖像進行分塊，相鄰塊之間設置一定的重疊度。評估每一個塊內最大的連通分量的像元數與塊內像元總數的比值θ，如果該值大于給定的閾值，則在塊內構造一個局部iForest來重新定義該窗口中像素的異常分數，并以此值代替全局值。不失一般性，本文設置的圖像塊的大小為20×20，塊間重疊度為4個像素，比值θ的閾值設置為0.3。若相鄰塊均有局部操作，則重疊區域的異常分數取相鄰塊的異常分數的平均值。

2 試驗驗證與分析

2.1 試驗數據集說明

利用5組高光譜數據集來驗證本文方法的性能，如圖3(a)、(b)、(c)和(d)為4幅分別命名為airport1、airport2、beach、harbor的AVIRIS高光譜傳感器獲取的影像[25]，4幅影像進行目標檢測前去掉了水吸收波段和信噪聲比低的波段；圖3(e)是用具有46個波段的CRI高光譜成像儀[26]獲取的真實高光譜數據影像(grassland)，圖3中的二值圖像為影像中對應的目標分布的參考圖。各數據的場景特點和目標空間分布相關信息描述如表1，比如圖3(a)為機場環境下有13架飛機目標，異常像元個數為144個。

2.2 試驗的實現與對比

根據采用的背景子空間估計方式的差異，以及是否進行主成分分析降維處理，對異常檢測方法有不同形式的實現。主成分背景子空間變換iForest異常檢測，記為(PCA subspaceiForest detector,PSF)，即背景子空間由主成分特征子空間法構造，影像經正交子空間背景抑制后，各像元的異常值的計算由全局iForest計算得到；局部細化的PSF，記為(locally refined PSF,LPSF)，即利用局部iForest對PSF得到的結果檢測結果圖再進行細化；聚類判別背景子空間變換iForest異常檢測，記為(cluster discriminant subspaceiForest detector,CDSF)，與PSF不同的在于背景子空間由典型類判別特征子空間法得到；局部細化的CDSF，記為(locally refined CDSF,LCDSF)；主成分降維的局部細化主成分背景子空間變換iForest異常檢測，記為(dimension reduced LPSF,DLPSF)，即對高光譜影像進行正交子空間背景抑制處理后用PCA主成分分析降維，之后各像元的異常值的計算與LPSF的計算方式相同。

圖3 高光譜影像及其對應的異常目標參考Fig.3 Hyperspectral images and their target reference images

表1 試驗數據描述及其試驗參數設置

為便于分析，將本文方法處理結果與其他全局化異常檢測方法的結果對比。包括不做任何預處理的全局iForest，全局RX(Global RX-GRX)，和經正交子空間背景抑制后的GRX，分別是主成分背景子空間變換GRX異常檢測(PCA background subspace GRX,Ps-GRX)和聚類判別背景子空間變換GRX異常檢測(cluster discriminant background subspace GRX,Cs-GRX)。

2.2.1 試驗的參數設置

本文所有基于iForest的算法中樹的數目t統一設置為100，樣本隨機采樣大小n為256。背景子空間參數設置方面。

主成分特征向量背景子空間參數k則通過AUC(area under the curve)[27]評價方法生成每一個試驗數據用不同值時PSF異常目標檢測結果。圖4結果展示的點線圖反映5幅影像采用的主成分個數與AUC值之間的關系，根據圖4試驗結果設置每一個試驗數據的最優參數見表1中“PCA子空間數目k”。CDSF方法中背景子空間的確定涉及兩個參數，分別為聚類數目和異常像元占影像面積的比例[28]。根據本文所用數據的特點，并進行多次試驗后選擇了效果最好的參數，各場景的上述兩參數設置見表1最后兩行。

圖4 PCA子空間參數k的確定Fig.4 Determination of PCA subspace parameter k

2.2.2 試驗對比分析

2.2.2.1 背景正交子空間增強結果對比

基于表1給定的參數，圖5展示了9種不同異常檢測方法在5個高光譜數據上進行處理后異常目標響應強度的灰度圖。其中基于iForest的方法異常值根據式(4)的定義將小于0.5的像元全部設為0。

圖5 異常目標檢測結果Fig.5 Detection maps for anomaly targets

對照目標參考圖計算目標與背景分離的箱狀表達圖[29](圖6)，箱狀圖表現了異常檢測特征圖像中背景和目標響應值的分散情況，紅色和藍色分別代表異常目標與背景檢測值的主體分布，距離越遠，則說明算法的性能越好。在計算GRX的箱狀圖之前先將其檢測得到的特征圖進行歸一化。由圖6中可以看出iForest相對于GRX在harbor和grassland兩幅圖像的目標響應值的數字范圍總體上好于GRX，但與GRX相同的是目標的響應最大值比背景的小；在airport2中GRX結果明顯優于iForest；從Ps-GRX和Cs-GRX的分離圖中可以看出正交子空間投影對GRX的數值動態變化范圍影響甚微，甚至結果不如GRX；兩種類型的子空間增強方法中PSF在各試驗數據的背景與目標的分離性比CDSF好，也比其他方法的分離性好。

圖6 目標與背景分離Fig.6 Separability maps of target and background

表2計算了圖5各結果的AUC(area under curve)值，AUC越接近1.0，檢測方法可靠性越高。采用子空間分析的iForest的AUC值比GRX，以及同樣采用子空間分析的Ps-GRX和Cs-GRX的都高，也比經典iForest高。以airport1為例，PSF的AUC值相對于經典iForest，GRX、Ps-GRX和Cs-GRX分別提高了9.38%、12.8%、14.52%和13.22%；而CDSF則分別提高了0.65%、3.8%、5.4%和4.2%。精度提高非常明顯說明經子空間背景抑制處理后異常目標檢測變得更容易。同時也注意到子空間分析并沒有提高經典GRX的結果。

表2 算法的AUC分數

2.2.2.2 局部iForest對異常檢測結果的影響

由圖5檢測圖可以明顯地看到全局與局部相結合的異常檢測算法的誤檢減少更為明顯。結合表2，以airport1為例，LPSF相對于PSF的AUC值得到了提高1.43%；LCDSF相對于CDSF的AUC值得到了提高2.82%。這說明局部處理進一步提高了精度。

2.2.2.3 主成分分析特征降維對iForest異常檢測算法的影響

用不同維數的特征對不同數據進行異常檢測的AUC統計值如圖7所示，airport1、airport2、beach、harbor和grassland 5個場景分別在保留2、4、6、2、1個特征時取得最大值。可以看出背景抑制后的影像在有限個主成分特征上獲得了更好的檢測結果，目標的形狀更為完整。其對應的AUC值見表2，以airport1為例，DLPSF相對于LPSF的AUC值提高了0.38%。

2.2.2.4 時間效率對比

在CPU為3.60 GHz、內存為8 GB配置的計算機上比較了異常檢測算法的運行時間。所有算法都在Matlab軟件上實現，表3給出每種方法在每個數據上的運行時間，可以看出iForest需要耗費比GRX更長的時間，時間消耗的主要原因是孤立樹建立和測試過程；正交子空間投影簡化了建樹過程，因而其消耗時間比經典iForest少，主成分分析定義的比聚類判別分析的耗費時間更短。局部細化計算需要額外的時間，降維能使耗費時間進一步減少。

圖7 AUC分數隨保留特征維數的變化Fig.7 Changes of AUC scores with retained feature dimensions

3 結 論

不同于經典高光譜圖像異常目標檢測方法偏向于正常背景樣本的描述，本文采用iForest異常檢測思想著重于孤立異常點。但并不是直接將iForest算法應用高光譜圖像，而是提出了在基于特征子空間正交分析的背景抑制和主成分降維處理后，將圖像全局iForest異常值計算初步結果和局部分塊計算結果結合。綜合試驗分析與討論得出如下結論：①以主成分特征子空間分析為基礎的LPSF的異常特征值與背景具有更好的可分性，其能夠獲得更高目標計算精度和計算效率；②全局和局部結合的iForest異常目標檢測更靈活地檢測出了感興趣的目標，并獲得了良好的檢測精度，豐富了算法的內涵；③主成分特征降維使得iForest算法較好地適用于高光譜圖像分析。但從圖6中可以看到檢測結果圖中存在一定的噪聲，這是論文未來需要解決的問題。

表3 試驗數據集的算法運行時間