孟令博,耿修瑞

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

0 引 言

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

1 主峭度分析算法簡介

(1)

(2)

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

算法 1 PKA算法輸入?yún)?shù):讀入圖像數(shù)據(jù)X={x1,x2,…,xN};∥首先對輸入的數(shù)據(jù)進(jìn)行白化∥①計(jì)算均值向量μ=mean(X);②按照F=ED-1/2ET計(jì)算圖像的白化矩陣,其中,E表示圖像協(xié)方差矩陣的特征矩陣,D=diag(λ1,λ2,…,λL)表示對應(yīng)的圖像協(xié)方差矩陣的特征值;③得到白化后的圖像矩陣^X=FT(X-μ),^X記為白化后的圖像數(shù)據(jù)∥計(jì)算協(xié)峭度張量矩陣∥④計(jì)算協(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 改進(jìn)的主峭度分析算法

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

(3)

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

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

(4)

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

(5)

可以證明

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

(6)

證明過程如下。

(7)

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

令y=uTX,則

(8)

那么

(9)

對于的方差

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

(10)

令

(11)

式(11)中第j個(gè)元素為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)

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

(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)

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

(38)

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

(39)

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

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

算法 2 改進(jìn)的PKA算法輸入?yún)?shù):讀入圖像數(shù)據(jù)X={x1,x2,…,xN};∥首先對輸入的數(shù)據(jù)進(jìn)行白化∥①計(jì)算均值向量μ=mean(X);②按照F=ED-1/2ET計(jì)算圖像的白化矩陣,其中,E表示圖像協(xié)方差矩陣的特征向量矩陣,D=diag(λ1,λ2,…,λL)表示對應(yīng)的圖像協(xié)方差矩陣的特征值;③得到白化后的圖像矩陣^X=FT(X-μ),^X記為白化后的圖像數(shù)據(jù)∥計(jì)算協(xié)峭度張量矩陣∥④計(jì)算協(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í)驗(yàn)

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

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

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

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

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

圖1 兩種算法的計(jì)算時(shí)間變化趨勢Fig.1 Calculating time trend of two algorithms

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

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

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

圖2 實(shí)驗(yàn)中待檢測的小目標(biāo)信息Fig.2 Small target information to be detected in experiment

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

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

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

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

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

4 結(jié) 論

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