改進的主峭度分析算法及其在高光譜圖像小目標檢測中的應用

孟令博,耿修瑞

(1.中國科學院空間信息處理與應用系統技術重點實驗室,北京 100190;2.中國科學院電子學研究所,北京 100190; 3.中國科學院大學,北京 100049)

0 引 言

主成分分析(principal component analysis,PCA)算法是一種廣泛的應用于數據壓縮和降維的方法[1]。它以數據的方差為指標進行數據降維,數據的協方差矩陣比較大的特征值對應的特征方向即為數據的投影方向。但小目標由于含有較小的信息量(方差)容易被忽略掉。所以,這種降維方法可能會不利于對高光譜圖像的分析和處理[2]。同樣的,最大噪聲成分(maximum noise fraction,MNF)和基于噪聲調整的主成分分析(noise-adjusted principal component,NAPC)算法,也存在同樣的問題[3-5]。關于這一問題的典型例子是對于AVIRIS傳感器獲取的Nevada Cuprite礦區的高光譜圖像的處理和分析,這可能是由于該礦區含有的礦物種類較多,且每種礦物的樣本較少,難以得到可靠的統計信息,從而可能會被基于二階統計信息的特征提取方法所忽略[5]?；诟唠A統計的分析方法在信號處理和分析領域的應用,為我們提供了新的思路。雖然這些小目標的存在對于整個數據的方差影響不大,但是這些小目標的存在對于數據的峭度值影響較大?；讵毩⒊煞址治?independent component analysis,ICA)的降維方法,即是一種利用數據的高階統計特性進行降維的技術,并且在高光譜圖像處理領域得到了越來越多的關注和研究[6-9]。在ICA的眾多解法中,基于固定點迭代的算法(FastICA),以收斂速度快、分離效果好而被廣泛的使用,且FastICA算法使用多種目標函數來使其達到性能最佳化[10-15]。近年來,文獻[16]提出了一種基于峭度指標的特征提取算法——主峭度分析(principal kurtosis analysis,PKA)算法。該算法通過引入協峭度張量來計算數據在某個方向上的峭度值,在迭代求解峭度極值的過程中不需要所有像元的參與,且也不需要選擇步長參數。它將求解偏度極值方向的問題轉化為求解協峭度張量的特征值和特征向量的問題。算法的收斂性能較好,且具有較好的穩定性。和基于峭度指標的FastICA算法相比,PKA算法可以降低算法的計算復雜度、減少計算時間[16]。但該方法使用歸一化的四階中心矩來計算峭度,算法在某些情況下的收斂性能不穩定,且收斂速度較慢,因而本文提出了一種新的基于協峭度張量的峭度計算方法,該方法具有更快的收斂速度,對于不同數據的適應性更好。為了測試使用MPKA降維后的小目標探測效果,仿真實驗部分采用約束最小能量算子(constrained energy minimization,CEM)[17-18]對降維后的數據進行小目標探測。

1 主峭度分析算法簡介

(1)

(2)

該算法采用固定點迭代法求解峭度極值方向,算法Matlab風格的計算流程如下。

算法 1 PKA算法輸入參數:讀入圖像數據X={x1,x2,…,xN};∥首先對輸入的數據進行白化∥①計算均值向量μ=mean(X);②按照F=ED-1/2ET計算圖像的白化矩陣,其中,E表示圖像協方差矩陣的特征矩陣,D=diag(λ1,λ2,…,λL)表示對應的圖像協方差矩陣的特征值;③得到白化后的圖像矩陣^X=FT(X-μ),^X記為白化后的圖像數據∥計算協峭度張量矩陣∥④計算協偏度張量K=1N∑Ni=1^xi ^xi ^xi ^xi∥搜索投影方向∥⑤ 令P=Ip×p⑥ for i=1:p 初始化向量ui(0),并單位化 令迭代次數k=1⑦ while(不滿足停止條件)dou(k+1)i=P(K×1(Pu(k)i)×2(Pu(k)i)×3(Pu(k)i))u(k+1)i=u(k+1)i‖u(k+1)i‖k=k+1 end⑧ 令ui=u(k+1)i,U=[u1 u2 … ui];⑨ P=I-UUT end

2 改進的主峭度分析算法

主峭度分析算法采用如下的歸一化的四階中心矩來計算峭度,即

(3)

而另一種對于零均值的隨機變量的峭度計算方法更為推薦和使用[8,13,19-21],其表示為

kurt(y)=E{y4}-3(E{y2})

(4)

改進的主峭度分析算法使用式(4)的表示方法。相應地,數據X在u方向上的峭度值可表示為

(5)

可以證明

E{(uTX)4}-3(E{(uTX)2})

(6)

證明過程如下。

(7)

(1)首先計算E{(uTX)4}-3(E{(uTX)2})

令y=uTX,則

(8)

那么

(9)

對于的方差

D(y)=E{y2}=E{uTXXTu}=1

(10)

令

(11)

式(11)中第j個元素為zj(1≤j≤N),有

(12)

令

(13)

則

(14)

所以

(15)

綜上有

E{(uTX)4}-3(E{(uTX)2})=

(16)

(17)

則S的任一元素為

(18)

令

(19)

則D的任一元素為

(20)

(21)

又

(22)

則

(23)

那么

(24)

又u為單位向量

所以

(25)

比較式(16)和式(25)可知,式(6)成立。

采用拉格朗日乘數法求解式(5)的極值優化問題,建立拉格朗日函數方程為

(26)

式中,λ≠0,函數兩邊分別對u和λ求導,有

(27)

令

(28)

得到如下方程組,即

(29)

采用固定點迭代法求解上述方程組,得到迭代公式為

(30)

(31)

X1=P1X={P1x1,P1x2,P1xN}

(32)

利用式(5),圖像數據在u方向上的峭度值為

(33)

相應的拉格朗日函數為

(34)

式中,λ≠0,函數兩邊分別對u和λ求導,有

(35)

令

(36)

得到如下方程組

(37)

采用固定點迭代法求解上述方程組,得到迭代公式為

(38)

根據上述討論,可以推導出搜索第(k+1)個向量時的迭代公式為

(39)

式中,Pk=I-UUT,U由已經搜索到的方向向量組成,U=[u1u2…uk]。

改進后PKA算法Matlab風格的計算流程如下。

算法 2 改進的PKA算法輸入參數:讀入圖像數據X={x1,x2,…,xN};∥首先對輸入的數據進行白化∥①計算均值向量μ=mean(X);②按照F=ED-1/2ET計算圖像的白化矩陣,其中,E表示圖像協方差矩陣的特征向量矩陣,D=diag(λ1,λ2,…,λL)表示對應的圖像協方差矩陣的特征值;③得到白化后的圖像矩陣^X=FT(X-μ),^X記為白化后的圖像數據∥計算協峭度張量矩陣∥④計算協偏度張量,K=1N∑Ni=1^xi ^xi ^xi ^xi∥搜索投影方向∥⑤ 令P=Ip×p⑥ for i=1:p 初始化向量ui(0),并單位化 令迭代次數k=1⑦ while(不滿足停止條件)dou(k+1)i=P(K×1(Pu(k)i)×2(Pu(k)i)×3(Pu(k)i)-3Pu(k)i)u(k+1)i=u(k+1)i‖u(k+1)i‖k=k+1 end⑧ 令ui=u(k+1)i,U=[u1 u2 … ui];⑨ P=I-UUT end

3 仿真實驗

仿真實驗部分分為兩部分:一部分為基于模擬數據的仿真實驗,一部分為基于真實數據的仿真實驗。模擬數據實驗部分用于比較改進后的PKA算法和原始的PKA算法,主要比較兩者的收斂速度和計算時間?；谡鎸崝祿姆抡鎸嶒炛饕容^本文提出的特征提取算法對小目標探測的有效性。

3.1 模擬數據

本小節設計了一組基于模擬數據的仿真實驗來比較改進后的PKA算法和原始的PKA算法的收斂速度。收斂速度的比較分為兩個部分:一是算法的計算時間,由于兩種算法僅是迭代過程不同,因此計算時間的比較只統計迭代過程需要的計算時間;二是算法迭代過程實際的迭代次數。首先隨機生成了一組圖像尺寸大小為500×500,波段數分別為5,10,15,20,25,30,35,40,45,50的模擬圖像數據用于特征提取,分別記錄兩種算法的計算時間及相應的迭代次數。然后生成了一組圖像波段數為30,圖像尺寸大小分別為100×100,200×200,…,1 000×1 000的模擬數據用于測試算法,同樣的記錄兩種算法的計算時間及相應的迭代次數。這兩種算法的迭代次數如表1所示。

表1 PKA算法和改進的PKA算法迭代次數比較Table 1 Comparison of iterations between PKA algorithm and improved PKA algorithm

從表1可以看出,隨著波段數的增加,PKA算法和改進的PKA算法所需要的迭代次數均增加,但是PKA算法的增加速度明顯快于改進后的PKA算法;但是算法迭代次數并不一定隨著像元數目的增加而增多,這是因為這兩種算法的迭代過程不需要所有像元的參與,僅需要協峭度張量的參與,而協峭度張量的大小僅僅和圖像的波段數有關。表1顯示,改進后的PKA算法的迭代次數遠少于原始PKA算法的迭代次數,從而改進后的PKA算法的收斂速度確實快于PKA算法的收斂速度。此外,實驗過程中還發現,PKA算法在求解同一個數據的不同極值方向時,迭代次數變化比較大,而改進后的PKA算法則比較穩定。這兩種算法相應的計算時間曲線如圖1所示。

圖1 兩種算法的計算時間變化趨勢Fig.1 Calculating time trend of two algorithms

由圖1可知,在相同的圖像尺寸及波段數目條件下,改進的PKA算法的計算時間遠少于原始的PKA算法的時間。在圖像尺寸不變的情況下,隨著波段數目的增加,兩種算法的計算時間均增加;但是在圖像波段數目一定的情況下,圖像尺寸的變化對計算時間的影響不大。綜上所述,改進的PKA算法相比較于原始的PKA算法具有更快的收斂,節省計算時間。

3.2 真實數據

高光譜遙感技術的重要優勢之一,即是對于低概率出露目標的有效探測?；诮y計分析的角度出發,這類低概率出露目標在圖像中的出露面積較小,甚至以亞像元的形式存在。本文提出的基于協峭度張量的特征提取技術特別有利于這類目標的探測。仿真實驗使用envi自帶的Cuprite數據用于測試特征提取后的結果對小目標的探測效果。該數據中礦物Buddingtonite具有較低的分布概率,實驗選擇該礦物作為感興趣小目標。Buddingtonite礦物的光譜可以從envi自帶的光譜庫中獲得,該小目標的光譜如圖2(a)所示。根據光譜庫中的光譜,利用光譜匹配算法使用全波段數據進行光譜匹配,制作了Buddingtonite的地面真值圖,如圖2(b)所示。

圖2 實驗中待檢測的小目標信息Fig.2 Small target information to be detected in experiment

實驗將PKA算法與本文提出的MPKA算法進行比較,利用這2種算法進行特征提取,分別選擇了前5、10、15個特征成分,然后將提取的特征成分使用CEM方法進行目標檢測,基于各個方法提取的特征成分用于CEM目標檢測得到的結果如圖3示。

圖3 MPKA和PKA算法提取的特征成分用于CEM檢測的結果Fig.3 Results of CEM detection with feature components extracted by the MPKA and PKA algorithm

由圖3可知,隨著選擇的特征成分數目增多,得到的CEM檢測結果都在變好。由于PKA和MPKA都是基于峭度指標的特征提取算法,因而它們的異常檢測結果相當。但是,MPKA算法在計算時間上占有優勢。PKA和MPKA算法搜索峭度極值方向時的迭代次數如表2所示。

表2 PKA和MPKA算法搜索峭度極值方向時的迭代次數Table 2 Number of iterations in searching the direction of kurtosis extremum between PKA and MPKA algorithm

由表2可知,MPKA算法迭代次數較少且相對穩定。算法提取15個特征成分花費的時間分別為:PKA算法計算時間約為25.947 2 s,MPKA算法的計算時間約為4.306 7 s,可見改進后的PKA算法計算時間更少,提高了計算效率。

4 結 論

本文提出了改進的主峭度分析算法——MPKA算法,并在此基礎上提出了基于主峭度分析算法的小目標檢測算法。MPKA算法使用了更為常用的峭度計算方式,從而使算法的收斂過程更加穩定而快速。實驗結果表明,相比PKA算法,MPKA算法需要的迭代次數更少、的計算時間更短,且利用該方法進行降維后的數據能夠有效地檢測出圖像中的小目標。