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

孟令博,耿修瑞

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

0 引 言

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

1 主峭度分析算法簡介

(1)

(2)

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

算法 1 PKA算法輸入?yún)?shù):讀入圖像數(shù)據(jù)X={x1,x2,…,xN};∥首先對輸入的數(shù)據(jù)進行白化∥①計算均值向量μ=mean(X);②按照F=ED-1/2ET計算圖像的白化矩陣,其中,E表示圖像協(xié)方差矩陣的特征矩陣,D=diag(λ1,λ2,…,λL)表示對應(yīng)的圖像協(xié)方差矩陣的特征值;③得到白化后的圖像矩陣^X=FT(X-μ),^X記為白化后的圖像數(shù)據(jù)∥計算協(xié)峭度張量矩陣∥④計算協(xié)偏度張量K=1N∑Ni=1^xi ^xi ^xi ^xi∥搜索投影方向∥⑤ 令P=Ip×p⑥ for i=1:p 初始化向量ui(0),并單位化 令迭代次數(shù)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)的表示方法。相應(yīng)地,數(shù)據(jù)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)成立。

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

(26)

式中,λ≠0,函數(shù)兩邊分別對u和λ求導(dǎo),有

(27)

令

(28)

得到如下方程組,即

(29)

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

(30)

(31)

X1=P1X={P1x1,P1x2,P1xN}

(32)

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

(33)

相應(yīng)的拉格朗日函數(shù)為

(34)

式中,λ≠0,函數(shù)兩邊分別對u和λ求導(dǎo),有

(35)

令

(36)

得到如下方程組

(37)

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

(38)

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

(39)

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

改進后PKA算法Matlab風(fēng)格的計算流程如下。

算法 2 改進的PKA算法輸入?yún)?shù):讀入圖像數(shù)據(jù)X={x1,x2,…,xN};∥首先對輸入的數(shù)據(jù)進行白化∥①計算均值向量μ=mean(X);②按照F=ED-1/2ET計算圖像的白化矩陣,其中,E表示圖像協(xié)方差矩陣的特征向量矩陣,D=diag(λ1,λ2,…,λL)表示對應(yīng)的圖像協(xié)方差矩陣的特征值;③得到白化后的圖像矩陣^X=FT(X-μ),^X記為白化后的圖像數(shù)據(jù)∥計算協(xié)峭度張量矩陣∥④計算協(xié)偏度張量,K=1N∑Ni=1^xi ^xi ^xi ^xi∥搜索投影方向∥⑤ 令P=Ip×p⑥ for i=1:p 初始化向量ui(0),并單位化 令迭代次數(shù)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 仿真實驗

仿真實驗部分分為兩部分:一部分為基于模擬數(shù)據(jù)的仿真實驗,一部分為基于真實數(shù)據(jù)的仿真實驗。模擬數(shù)據(jù)實驗部分用于比較改進后的PKA算法和原始的PKA算法,主要比較兩者的收斂速度和計算時間。基于真實數(shù)據(jù)的仿真實驗主要比較本文提出的特征提取算法對小目標探測的有效性。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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