孟令博，耿修瑞，楊煒暾

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

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

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

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

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

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

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

(1)

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

(2)

(3)

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

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

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

(4)

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

假設(shè)觀測信號X是源信號S線性混合而成，經(jīng)過向量w可以分離出一個源信號分量wTX。一般情況下，在處理觀測信號之前會進行數(shù)據(jù)白化，因此，不妨假設(shè)，觀測信號X的均值向量為0，協(xié)方差矩陣為單位陣，令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，基于峭度指標(biāo)的FastICA算法的目標(biāo)函數(shù)可以最終表示為

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

(6)

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

(7)

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

綜合上文，基于峭度指標(biāo)的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次迭代甚至幾次迭代內(nèi)就會得到最優(yōu)解。但它的運算量依然很大，而且循環(huán)迭代求解過程需要全部數(shù)據(jù)的參與，因此下面提出一種改進的FastICA算法以減少運算量，新算法引入?yún)f(xié)峭度張量的概念使得循環(huán)迭代過程中不需要所有像元的參與。

2 改進的FastICA算法

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

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

(9)

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

(10)

2.2 改進的FastICA算法原理

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

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

(11)

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

(12)

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

改進的FastICA算法偽代碼1)首先,計算圖像的協(xié)峭度張量K=1n∑ni=1xixixixi; 2)令P=I,P為正交投影矩陣;3)for i=1∶p,p為需要提取的獨立成分數(shù);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即為特征提取后的結(jié)果。

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

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

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

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

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

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

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

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

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

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

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

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

本節(jié)使用真實的多光譜圖像數(shù)據(jù)比較FastICA算法及改進FastICA算法的計算效率。該多光譜數(shù)據(jù)為Landsat5 TM反射率數(shù)據(jù)，獲取時間為2001年6月30日，該數(shù)據(jù)為江西北部的鄱陽湖部分區(qū)域數(shù)據(jù)，如圖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算法提取結(jié)果Fig.5 The feature extraction results of theoriginal FastICA and improved FastICA

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

5 總結(jié)

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

附錄：定理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)