一種改進(jìn)的基于峭度指標(biāo)的FastICA算法*

孟令博，耿修瑞，楊煒暾

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

多/高光譜遙感圖像通常包含多個(gè)光譜波段[1]，可提供關(guān)于地物的豐富的光譜信息，但同時(shí)也會(huì)造成一定的數(shù)據(jù)冗余及數(shù)據(jù)處理過(guò)程的計(jì)算負(fù)擔(dān)[2]。因此，在很多情況下需要引入數(shù)據(jù)降維或者特征提取的手段對(duì)數(shù)據(jù)進(jìn)行預(yù)處理[3]。主成分分析[4-5](principal component analysis，PCA)是一種最常用的特征提取方法，并常被應(yīng)用于遙感圖像的數(shù)據(jù)降維。PCA以圖像方差為指標(biāo)進(jìn)行特征提取，而沒(méi)有考慮噪聲在光譜特征空間的分布，因此降維結(jié)果易受具有較大方差的噪聲影響。Green等[6]提出基于圖像信噪比的特征提取方法，即最大噪聲成分(maximum noise fraction，MNF)算法，以解決這一問(wèn)題。Lee等[7]提出噪聲調(diào)整的主成分分析方法(NAPC)，也能夠減少噪聲對(duì)特征提取結(jié)果的影響，它基本與MNF等價(jià)。此外，還有考慮了圖像空間特征的最大/最小自相關(guān)因子分析算法，它基于圖像數(shù)據(jù)的自相關(guān)性而不是方差，因此不受向量幅值的影響。

以上這些特征提取算法都是基于數(shù)據(jù)的二階統(tǒng)計(jì)特性，它一般假設(shè)圖像數(shù)據(jù)服從高斯分布。但實(shí)際的遙感圖像數(shù)據(jù)往往并不滿(mǎn)足這一條件，所以基于二階統(tǒng)計(jì)特性的特征提取方法很多情況下并不能有效反映數(shù)據(jù)中的感興趣結(jié)構(gòu)。比如，當(dāng)感興趣目標(biāo)為小目標(biāo)時(shí)，基于二階統(tǒng)計(jì)特性的特征提取方法往往將其作為不重要的信息舍棄。而基于高階統(tǒng)計(jì)特性的特征提取方法就能很好地解決這類(lèi)問(wèn)題[8-9]。高階統(tǒng)計(jì)特性在遙感圖像處理中有很多應(yīng)用，例如目標(biāo)自動(dòng)識(shí)別、異常檢測(cè)、光譜解混、特征提取等[10-12]?；诟唠A統(tǒng)計(jì)特性的特征提取方法絕大多數(shù)是基于獨(dú)立成分分析(independent component analysis，ICA)的方法及其改進(jìn)算法[13-15]，其中FastICA算法以其快速的收斂速度和較高的計(jì)算效率得到廣泛應(yīng)用。FastICA算法有4種優(yōu)化指標(biāo)，其中基于峭度指標(biāo)的FastICA算法具有明確的物理意義，且相比基于偏度指標(biāo)的FastICA算法有更好的收斂性能和更快的收斂速度，因此應(yīng)用最為廣泛[16-20]。但是，在求解各個(gè)獨(dú)立成分時(shí)每一步迭代過(guò)程均需所有像元的參與，所以當(dāng)圖像尺寸較大時(shí)往往計(jì)算量很大且耗時(shí)較長(zhǎng)。然而，遙感圖像尺寸一般較大，此時(shí)FastICA算法的計(jì)算優(yōu)勢(shì)大打折扣。

因此，本文提出一種基于峭度指標(biāo)的FastICA算法的等價(jià)算法。通過(guò)引入?yún)f(xié)峭度張量，將FastICA的固定點(diǎn)迭代問(wèn)題轉(zhuǎn)化為代數(shù)形式的張量計(jì)算，在每次迭代中避免所有像元的參與，因而大大降低計(jì)算復(fù)雜度。

1 基于峭度指標(biāo)的FastICA算法

ICA算法是一種從盲源分離技術(shù)發(fā)展而來(lái)的算法[15]，它能夠?qū)⒍嗑S觀測(cè)信號(hào)分解為盡量相互統(tǒng)計(jì)獨(dú)立的獨(dú)立分量，分離出的分量是源信號(hào)的一種近似，因而在圖像處理、語(yǔ)音識(shí)別、遠(yuǎn)程通信等領(lǐng)域得到廣泛應(yīng)用。

為了更為方便和簡(jiǎn)潔地說(shuō)明ICA，假設(shè)觀測(cè)信號(hào)是源信號(hào)的線(xiàn)性變換。假設(shè)觀測(cè)信號(hào)為

(1)

它由相互獨(dú)立且均值為零的源信號(hào)

(2)

(3)

ICA可以描述為在源信號(hào)和混合矩陣都未知的情況下，根據(jù)觀測(cè)矩陣X找到一個(gè)分離矩陣W能夠盡可能地分離出源信號(hào)S。FastICA是ICA算法中最為流行，應(yīng)用最廣泛的一種算法。相比于其他算法，F(xiàn)astICA算法計(jì)算效率高，運(yùn)算量小，相對(duì)容易實(shí)現(xiàn)。

對(duì)于零均值的隨機(jī)變量x，其峭度可以表示為

kurt(x)=E{x4}-3(E{x2})2，

(4)

顯然，對(duì)于高斯信號(hào)，它的峭度值為0。

假設(shè)觀測(cè)信號(hào)X是源信號(hào)S線(xiàn)性混合而成，經(jīng)過(guò)向量w可以分離出一個(gè)源信號(hào)分量wTX。一般情況下，在處理觀測(cè)信號(hào)之前會(huì)進(jìn)行數(shù)據(jù)白化，因此，不妨假設(shè)，觀測(cè)信號(hào)X的均值向量為0，協(xié)方差矩陣為單位陣，令y=wTX，則它的峭度計(jì)算可以進(jìn)行如下簡(jiǎn)化[20]

kurt(y)=E{(wTX)4}-3(E{(wTX)2})2

=E{(wTX)4}-3(wTXXTw)2

=E{(wTX)4}-3‖w‖4，

(5)

式中:w滿(mǎn)足約束條件‖w‖=1，基于峭度指標(biāo)的FastICA算法的目標(biāo)函數(shù)可以最終表示為

J(w)=E{(wTX)4}-3‖w‖4+F(‖w‖2)，

(6)

式中:F是關(guān)于約束條件的懲罰函數(shù)。固定點(diǎn)算法求解式(6)可以得到迭代公式

(7)

式中：“⊙”符號(hào)代表矩陣點(diǎn)乘，表示2個(gè)矩陣對(duì)應(yīng)位置的每個(gè)元素相乘。下同。

綜合上文，基于峭度指標(biāo)的FastICA算法步驟表示如下

1)首先隨機(jī)生成一個(gè)模為1的向量w(0)，并令k=0；

2)將其代入迭代公式

w(k+1)=

4)如果|(w(k+1))Tw(k)|的值不足夠接近1，則令k=k+1，跳到步驟2)；否則輸出向量w(k+1)。

上述過(guò)程得到的向量w(k+1)為分離矩陣W的一個(gè)列向量，(w(k+1))TX可以分離出一個(gè)源信號(hào)的近似，即一個(gè)獨(dú)立分量。為了分離出p個(gè)獨(dú)立分量，需要將上述步驟運(yùn)行p次。同時(shí)為保證每次分離出的獨(dú)立分量是不同的，要對(duì)w(k+1)進(jìn)行正交投影，使其與已求得的其他向量均正交，可通過(guò)下式進(jìn)行正交化

(8)

上述算法具有較快的收斂速度，一般情況下，它在10～100次迭代甚至幾次迭代內(nèi)就會(huì)得到最優(yōu)解。但它的運(yùn)算量依然很大，而且循環(huán)迭代求解過(guò)程需要全部數(shù)據(jù)的參與，因此下面提出一種改進(jìn)的FastICA算法以減少運(yùn)算量，新算法引入?yún)f(xié)峭度張量的概念使得循環(huán)迭代過(guò)程中不需要所有像元的參與。

2 改進(jìn)的FastICA算法

改進(jìn)的FastICA算法是對(duì)原始FastICA算法的迭代公式的改進(jìn)，通過(guò)引入?yún)f(xié)峭度張量，在降低算法計(jì)算復(fù)雜度的同時(shí)避免所有像元參與算法的迭代過(guò)程。下面，首先介紹協(xié)峭度張量的概念。

2.1 協(xié)峭度張量的計(jì)算

(9)

式中:“°”代表外積，xi°xi°xi°xi可以生成一個(gè)p維的4階張量。根據(jù)式(9)可知，

(10)

2.2 改進(jìn)的FastICA算法原理

定理2.1對(duì)于任意向量w=[w1,…,wp]T，有下面的結(jié)論成立：

X[(XTw)⊙(XTw)⊙(XTw)]/n

(11)

則式(7)的迭代公式變?yōu)?/p>

(12)

式中：×1×2×3代表模乘，下同。

改進(jìn)的FastICA算法偽代碼1)首先,計(jì)算圖像的協(xié)峭度張量K=1n∑ni=1xixixixi; 2)令P=I,P為正交投影矩陣;3)for i=1∶p,p為需要提取的獨(dú)立成分?jǐn)?shù);4)令k=0,并生成隨機(jī)單位向量w(k)i;5)將w(k)i進(jìn)行正交投影w(k)i=Pw(k)i;6)將w(k)i進(jìn)行迭代更新w(k+1)i=K×1w(k)i×2w(k)i×3w(k)i-3w(k)iw(k+1)i=w(k+1)i‖w(k+1)i‖ì?í????; 7)k=k+1;8)如果滿(mǎn)足停止條件,則輸出w(k+1)i,否則跳到5);9)令P=I-WWT,W為已求得的向量wi的集合W=[w1,…,wi];10)end;Y=WTX即為特征提取后的結(jié)果。

2.3 算法復(fù)雜度分析

比較基于峭度指標(biāo)的FastICA算法及其改進(jìn)算法的迭代公式可知，改進(jìn)的FastICA算法的迭代公式更加簡(jiǎn)潔，且不需要全部像元點(diǎn)的參與，因此會(huì)節(jié)省更多的計(jì)算時(shí)間。為了更加準(zhǔn)確地比較這兩種算法，需要對(duì)它們進(jìn)行詳細(xì)的計(jì)算復(fù)雜度分析。

3 基于模擬數(shù)據(jù)的實(shí)驗(yàn)分析

原始的FastICA算法和改進(jìn)后的FastICA算法雖然都使用固定點(diǎn)迭代算法，但是二者有本質(zhì)上的不同。原始的FastICA算法在迭代過(guò)程中需要全部像元的參與，改進(jìn)后的FastICA算法則依賴(lài)預(yù)先計(jì)算好的協(xié)峭度張量。為了更加客觀地比較兩種算法，F(xiàn)astICA算法和改進(jìn)的FastICA算法使用相同的初始向量。所有的仿真實(shí)驗(yàn)均在同一電腦上使用相同的操作系統(tǒng)完成。本文中的算法全部使用MATLAB7.1軟件實(shí)現(xiàn)。

3.1 計(jì)算時(shí)間和數(shù)據(jù)維數(shù)之間的關(guān)系

在本節(jié)中，設(shè)計(jì)一組實(shí)驗(yàn)考察數(shù)據(jù)的維數(shù)對(duì)兩種算法計(jì)算時(shí)間的影響。首先用MATLAB軟件中的rand函數(shù)產(chǎn)生10組模擬數(shù)據(jù)，每組數(shù)據(jù)尺寸為800×800，維數(shù)從2增加到15?？紤]到初始值的影響，對(duì)兩種算法設(shè)置相同的初始值，各運(yùn)行100次取平均計(jì)算時(shí)間進(jìn)行比較。圖1給出兩種算法的計(jì)算時(shí)間。從圖中可以看出隨著數(shù)據(jù)維數(shù)的增加，改進(jìn)的FastICA需要的計(jì)算時(shí)間急劇增加，當(dāng)數(shù)據(jù)維數(shù)超過(guò)11時(shí)，改進(jìn)的FastICA算法反而需要更多的計(jì)算時(shí)間。這也和上文分析的結(jié)果一致。圖2給出計(jì)算協(xié)峭度張量的時(shí)間占總的計(jì)算時(shí)間的比重，可以看出隨著數(shù)據(jù)維數(shù)的增加它所占的比重越來(lái)越大，甚至占用絕大部分的計(jì)算時(shí)間。因此，如果能夠準(zhǔn)確快速地估計(jì)數(shù)據(jù)的協(xié)峭度張量將極大地提高改進(jìn)的FastICA算法的計(jì)算效率。

圖1 原始FastICA算法和改進(jìn)的FastICA算法的計(jì)算時(shí)間隨波段數(shù)變化的趨勢(shì)Fig.1 Variations in the computing time of the original andimproved FastICA algorithms with the band number

圖2 不同數(shù)據(jù)維數(shù)條件下計(jì)算協(xié)峭度張量的時(shí)間占總的計(jì)算時(shí)間的比重Fig.2 Ratio of the calculation time for cokurtosistensor to the total calculation time

3.2 計(jì)算時(shí)間和像元數(shù)目之間的關(guān)系

為考察像元總的數(shù)目對(duì)兩種算法的影響，生成一組維數(shù)為6，像元數(shù)目從100×100到900×900的隨機(jī)模擬數(shù)據(jù)。同樣地，對(duì)兩種算法設(shè)置相同的初始值，統(tǒng)計(jì)其所需要的計(jì)算時(shí)間(100次平均)，如圖3所示。從圖中可以看出隨著像元數(shù)目的增多，兩種算法所需的計(jì)算時(shí)間都在逐漸增加，但改進(jìn)后的FastICA算法增加得更慢，且它需要的計(jì)算時(shí)間始終少于原始FastICA算法需要的計(jì)算時(shí)間。

圖3 原始FastICA算法和改進(jìn)的FastICA算法在不同像元數(shù)目條件下的計(jì)算時(shí)間比較Fig.3 Comparison of the computing time between the originalFastICA and improved FastICA at different pixel numbers

4 基于真實(shí)數(shù)據(jù)的實(shí)驗(yàn)分析

本節(jié)使用真實(shí)的多光譜圖像數(shù)據(jù)比較FastICA算法及改進(jìn)FastICA算法的計(jì)算效率。該多光譜數(shù)據(jù)為L(zhǎng)andsat5 TM反射率數(shù)據(jù)，獲取時(shí)間為2001年6月30日，該數(shù)據(jù)為江西北部的鄱陽(yáng)湖部分區(qū)域數(shù)據(jù)，如圖4所示。該圖像分辨率為30 m，尺寸為800×800，共包含6個(gè)波段，分別為485，569，660，840，1 676和2 223 nm。

圖4 鄱陽(yáng)湖第一波段圖像Fig.4 The first-band image of Poyang Lake

圖5 兩種FastICA算法提取結(jié)果Fig.5 The feature extraction results of theoriginal FastICA and improved FastICA

將原始的基于峭度指標(biāo)的FastICA算法和改進(jìn)后的FastICA算法應(yīng)用于該多光譜圖像，將特征提取結(jié)果進(jìn)行詳細(xì)對(duì)比分析發(fā)現(xiàn)，原始的基于峭度指標(biāo)的FastICA算法和改進(jìn)后的FastICA算法提取結(jié)果一致，特征提取結(jié)果如圖5所示。同時(shí)基于峭度指標(biāo)的FastICA算法需要的時(shí)間為8.529 2 s，而改進(jìn)后的FastICA算法僅僅需要2.275 5 s。因此，改進(jìn)后的FastICA算法在保持原有結(jié)果不變的同時(shí)計(jì)算速度遠(yuǎn)快于原始的FastICA算法。

5 總結(jié)

本文提出一種改進(jìn)的FastICA算法以降低FastICA的計(jì)算復(fù)雜度，它是原始FasICA算法的代數(shù)等價(jià)算法，能夠產(chǎn)生與之完全相同的特征提取結(jié)果。改進(jìn)的FastICA算法引入?yún)f(xié)峭度張量的概念，在迭代過(guò)程中不需要全部像元的參與，而僅僅需要協(xié)峭度張量的參與。當(dāng)數(shù)據(jù)維數(shù)不是特別高時(shí)，改進(jìn)的FastICA算法相比于原始的FastICA算法計(jì)算復(fù)雜度更低，能節(jié)省更多的計(jì)算時(shí)間。但當(dāng)數(shù)據(jù)維數(shù)很高時(shí)，改進(jìn)的FastICA算法實(shí)際上并不能節(jié)省計(jì)算時(shí)間。因此該方法更加適用于多光譜圖像的特征提取。改進(jìn)的FastICA算法中計(jì)算協(xié)峭度張量花費(fèi)絕大部分時(shí)間，因此尋求一種快速計(jì)算協(xié)峭度張量的方法可以極大地提高算法的計(jì)算效率，縮短計(jì)算時(shí)間。

附錄：定理2.1的證明

首先計(jì)算X[(XTw)⊙(XTw)⊙(XTw)]/n.

令

(A1)

式(A1)中第k個(gè)元素表示為

(A2)

則

(A3)

即

X[(XTw)⊙(XTw)⊙(XTw)]/n=

(A4)

為超對(duì)稱(chēng)張量.

(A5)

則S的任一元素為

(A6)

(A7)

(A8)

(A9)

那么，

(A10)