一種改進的基于峭度指標的FastICA算法*

孟令博，耿修瑞，楊煒暾

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

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

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

因此，本文提出一種基于峭度指標的FastICA算法的等價算法。通過引入協峭度張量，將FastICA的固定點迭代問題轉化為代數形式的張量計算，在每次迭代中避免所有像元的參與，因而大大降低計算復雜度。

1 基于峭度指標的FastICA算法

ICA算法是一種從盲源分離技術發展而來的算法[15]，它能夠將多維觀測信號分解為盡量相互統計獨立的獨立分量，分離出的分量是源信號的一種近似，因而在圖像處理、語音識別、遠程通信等領域得到廣泛應用。

為了更為方便和簡潔地說明ICA，假設觀測信號是源信號的線性變換。假設觀測信號為

(1)

它由相互獨立且均值為零的源信號

(2)

(3)

ICA可以描述為在源信號和混合矩陣都未知的情況下，根據觀測矩陣X找到一個分離矩陣W能夠盡可能地分離出源信號S。FastICA是ICA算法中最為流行，應用最廣泛的一種算法。相比于其他算法，FastICA算法計算效率高，運算量小，相對容易實現。

對于零均值的隨機變量x，其峭度可以表示為

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

(4)

顯然，對于高斯信號，它的峭度值為0。

假設觀測信號X是源信號S線性混合而成，經過向量w可以分離出一個源信號分量wTX。一般情況下，在處理觀測信號之前會進行數據白化，因此，不妨假設，觀測信號X的均值向量為0，協方差矩陣為單位陣，令y=wTX，則它的峭度計算可以進行如下簡化[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滿足約束條件‖w‖=1，基于峭度指標的FastICA算法的目標函數可以最終表示為

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

(6)

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

(7)

式中：“⊙”符號代表矩陣點乘，表示2個矩陣對應位置的每個元素相乘。下同。

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

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

2)將其代入迭代公式

w(k+1)=

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

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

(8)

上述算法具有較快的收斂速度，一般情況下，它在10～100次迭代甚至幾次迭代內就會得到最優解。但它的運算量依然很大，而且循環迭代求解過程需要全部數據的參與，因此下面提出一種改進的FastICA算法以減少運算量，新算法引入協峭度張量的概念使得循環迭代過程中不需要所有像元的參與。

2 改進的FastICA算法

改進的FastICA算法是對原始FastICA算法的迭代公式的改進，通過引入協峭度張量，在降低算法計算復雜度的同時避免所有像元參與算法的迭代過程。下面，首先介紹協峭度張量的概念。

2.1 協峭度張量的計算

(9)

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

(10)

2.2 改進的FastICA算法原理

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

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

(11)

則式(7)的迭代公式變為

(12)

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

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

2.3 算法復雜度分析

比較基于峭度指標的FastICA算法及其改進算法的迭代公式可知，改進的FastICA算法的迭代公式更加簡潔，且不需要全部像元點的參與，因此會節省更多的計算時間。為了更加準確地比較這兩種算法，需要對它們進行詳細的計算復雜度分析。

3 基于模擬數據的實驗分析

原始的FastICA算法和改進后的FastICA算法雖然都使用固定點迭代算法，但是二者有本質上的不同。原始的FastICA算法在迭代過程中需要全部像元的參與，改進后的FastICA算法則依賴預先計算好的協峭度張量。為了更加客觀地比較兩種算法，FastICA算法和改進的FastICA算法使用相同的初始向量。所有的仿真實驗均在同一電腦上使用相同的操作系統完成。本文中的算法全部使用MATLAB7.1軟件實現。

3.1 計算時間和數據維數之間的關系

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

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

圖2 不同數據維數條件下計算協峭度張量的時間占總的計算時間的比重Fig.2 Ratio of the calculation time for cokurtosistensor to the total calculation time

3.2 計算時間和像元數目之間的關系

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

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

4 基于真實數據的實驗分析

本節使用真實的多光譜圖像數據比較FastICA算法及改進FastICA算法的計算效率。該多光譜數據為Landsat5 TM反射率數據，獲取時間為2001年6月30日，該數據為江西北部的鄱陽湖部分區域數據，如圖4所示。該圖像分辨率為30 m，尺寸為800×800，共包含6個波段，分別為485，569，660，840，1 676和2 223 nm。

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

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

將原始的基于峭度指標的FastICA算法和改進后的FastICA算法應用于該多光譜圖像，將特征提取結果進行詳細對比分析發現，原始的基于峭度指標的FastICA算法和改進后的FastICA算法提取結果一致，特征提取結果如圖5所示。同時基于峭度指標的FastICA算法需要的時間為8.529 2 s，而改進后的FastICA算法僅僅需要2.275 5 s。因此，改進后的FastICA算法在保持原有結果不變的同時計算速度遠快于原始的FastICA算法。

5 總結

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

附錄：定理2.1的證明

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

令

(A1)

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

(A2)

則

(A3)

即

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

(A4)

為超對稱張量.

(A5)

則S的任一元素為

(A6)

(A7)

(A8)

(A9)

那么，

(A10)