(1.中國地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院， 武漢 430074；2.廣州航海高等?？茖W(xué)校 計算機與信息工程系， 廣州 510725；3.新疆維吾爾自治區(qū)地震局， 烏魯木齊 830011）

摘 要：ICA算法是求解盲源分離問題的有效算法。建立了ICA算法的數(shù)學(xué)模型，對模型的求解條件及多解性進行了分析。給出一種基于負熵極大的FastICA算法，討論該算法在地震信號去噪中的應(yīng)用。仿真實驗驗證了該算法的有效性。

關(guān)鍵詞：獨立分量分析算法；地震信號；去噪

中圖分類號：TP301.6文獻標(biāo)志碼：A

文章編號：1001-3695(2009)04-1432-03

FastICA algorithm and its application in seismic signal noise elimination

ZHANG Nian1，2，LIU Tian-you1，LI Jie1，3

（1.Institute of Geophysics Geomatics， China University of Geosciences， Wuhan 430074， China;2. Dept. of Computer Science Information Technology， Guangzhou Maritime College， Guangzhou 510725， China;3.Earthquake Administration of Xinjiang Uygur Autonomous Region， Urumchi 830011， China）

Abstract:ICA algorithm is a effective algorithm for blind source separation(BSS) problem.This paper established themathematical model of ICA algorithm ， and analyzed its solving condition and multiplicity solutions.Gave a Fast ICA algorithm based on maximum negentropy， discussed its application in seismic signal noise elimination. A simulation experiment verifies the effectiveness ofthe algorithm.

Key words:ICA(independent component analysis) algorithm; seismic signal; noise elimination

0 引言

獨立分量分析(independent component analysis，ICA)是伴隨BSS(blind source separation)問題而發(fā)展起來的一種信號分離算法。自Comon于1994年首次闡述ICA的概念并基于高階統(tǒng)計量構(gòu)造代價函數(shù)以來[1]，許多學(xué)者對該方法進行了研究，提出了一系列算法[2~5]。與傳統(tǒng)的盲源分離算法相比，ICA算法基于數(shù)據(jù)的高階統(tǒng)計量，分離得到的各分量不僅互不相關(guān)(二階統(tǒng)計獨立)，而且還盡可能高階統(tǒng)計獨立，故更能反映數(shù)據(jù)的本質(zhì)特征。因此，ICA算法在信號處理領(lǐng)域受到廣泛的關(guān)注，逐漸應(yīng)用于許多工程領(lǐng)域[6~9]。

地震勘探中存在著大量噪聲，隨著勘探工作的逐步深入，勘探條件更加復(fù)雜，勘探難度越來越大，目標(biāo)越來越隱蔽，采集條件越來越惡劣，使得采集的樣本信號中包含大量干擾信息，嚴(yán)重影響了地震剖面質(zhì)量，所以在使用地震資料前，必須進行去噪。

目前采用的地震信號去噪方法，如帶通濾波、F-K域消除、小波變換分頻去噪、徑向預(yù)測濾波、矢量分解去噪等[10～13]，在各類地震信號去噪中起到了很大作用，但也存在一些缺點，如需要犧牲部分有效信號，而且必須滿足一定的信噪比才可使用等。本文針對許多地震記錄強噪聲、弱信號的特點進行研究，嘗試采用ICA算法對地震信號去噪，提取有用信號。

1 ICA方法的基本原理

1.1 ICA方法的數(shù)學(xué)模型

圖1是ICA方法的簡化數(shù)學(xué)模型。

圖中S(t)=[s1(t)，s2(t)，…， sn(t)]T是由n個源信號構(gòu)成的n維矢量；X(t)=[x1(t)，x2(t)， …，xm(t)]T表示通過傳感器輸出的m個觀測信號。S(t)與X(t)的關(guān)系可用如下方程描述[14]：

X(t)=AS(t)(1)

其中：A為m×n維混合矩陣，含義為n個源信號通過線性組合得到m個觀測值。ICA方法要解決的問題是：在矩陣A未知、S(t)除知道其分量獨立外，無其他先驗知識的情況下，求解混合矩陣W，使得處理結(jié)果Y(t)=WX(t)中各分量盡可能互相獨立，且可以作為S(t)的有效估計值。

1.2 模型假設(shè)條件

在上述模型中，由于缺少混合矩陣A的結(jié)構(gòu)信息，會導(dǎo)致方程的多解。在應(yīng)用上述模型時，需作如下假設(shè)[15]：

a)源信號S(t)各分量均為平穩(wěn)隨機過程，且相互獨立。

b)源信號S(t)各分量中最多只能有一個服從高斯分布。

c)m≥n，即觀測信號的個數(shù)不少于源信號的個數(shù)。為便于計算，一般假設(shè)m=n。

d)混合矩陣A是列滿秩的。

1.3 模型解的不確定性

1.3.1 幅度的不確定性

對式(1)，可以改寫為

X(t)=AS(t)=AB-1BS(t)= A′S′(t)(2)

其中：B為一對角矩陣，A′=AB-1，S′(t)= BS(t)。在ICA算法中，惟一可以利用的信息是X(t)，很難確定由X(t)恢復(fù)出來的信號到底是對S(t)的估計，還是對S′(t)的估計，這導(dǎo)致了分離信號幅度的不確定性。

1.3.2 輸出順序的不確定性

將式(1)寫成如下形式：

X(t)=AS(t)=AD-1DS(t)= A″S″(t)(3)

其中:D為置換矩陣，即D的每一行每一列都只有一個元素為1，其余元素均為0；A″= AD-1，S″(t)= DS(t)，則無法確定恢復(fù)出來的信號是對S(t)的估計，還是對S″(t)的估計，這使得解的順序不確定。

考慮到上述不確定性，筆者在使用ICA算法求解BSS問題時，主要關(guān)心信號的波形，而對幅度和順序則不考慮。

2 基于負熵極大的FastICA算法

2.1 FastICA算法

FastICA算法是芬蘭赫爾辛基工業(yè)大學(xué)的Hyvarinen等人[3]最先提出來的。該算法基于非高斯性最大化原理，使用固定點(fixed-point)迭代理論尋找Wx的非高斯性最大值，故有時又稱為固定點算法。該算法采用牛頓迭代算法對觀測變量x的大量采樣點進行批處理，每次從觀測信號中分離一個獨立分量，是獨立分量分析的一種快速算法。為減少算法需要估計的參數(shù)，在運行FastICA算法之前，需要對數(shù)據(jù)進行預(yù)處理，即去均值和白化。

去均值又稱中心化，一般希望處理后的信號為零均值。其算法為

X=X-E{X}(4)

白化是對去均值后的觀測信號向量X施加一個線性變換V，使得新向量Z的各個分量間互不相關(guān)。白化過程將混合矩陣A變換成一個新的正交矩陣B，如下式：

E{ZZT}=E{BSSTB}=BE{SST}BT=BBT=I(5)

新的觀測向量Z由具有單位方差且互不相關(guān)的各分量zi組成，白化過程并沒有改變ICA問題的性質(zhì)。

預(yù)處理結(jié)束，就可以進行非高斯最大化計算(即分離過程)，此時必須先確定非高斯性的度量標(biāo)準(zhǔn)。非高斯性的度量標(biāo)準(zhǔn)主要有峭度(kurtosis)和負熵(negentropy)兩種，分別定義如下：

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

J(x)=∫p(x) log p(x)/[pG(x)]dx=HG(x)-H(x)(7)

由于峭度存在不穩(wěn)定性，用峭度作為目標(biāo)函數(shù)尋找獨立分量時，易受大幅度隨機脈沖干擾的影響，而負熵則顯得更為穩(wěn)健。

2.2 基于負熵極大的FastICA算法

2.2.1 目標(biāo)函數(shù)

負熵是一種有效度量隨機變量非高斯性的目標(biāo)函數(shù)，可以用來衡量隨機變量的獨立性。式(7)給出了負熵的一般定義?？梢钥闯?，負熵的值總是非負的，只有當(dāng)x為高斯分布時，負熵才為零。

負熵還有一個重要性質(zhì)：對所有可逆線性變換，負熵的值保持不變。故從統(tǒng)計學(xué)角度來講，負熵是度量隨機變量非高斯性的最優(yōu)工具之一。但是在負熵計算中涉及到對隨機變量概率密度的估計，其有效性依賴于對估計參數(shù)的正確選擇，且計算量很大，需要尋找更簡單的計算方法，對負熵進行合理而有效的近似。

一種傳統(tǒng)的對負熵的估計方法來源于對概率密度函數(shù)的Gram-Charlier的近似展開：

J(y)≈1/12k3(y)2+1/48 k4(y)2=1/12(E{y3})2+1/48(kurt(y))2(8)

該方法由于使用了峭度而使得估計的魯棒性不好。為解決此問題，很多時候采用基于最大熵原理的負熵近似：

J(y)≈ni=1ki[E{Gi(y)}-E{Gi(v)}]2(9)

其中：ki為一正常量；v為零均值、單位方差的高斯變量；函數(shù)Gi()為非二次函數(shù)。如下是被實踐證明性能較好的G()函數(shù)：

G1(u)=(1/a1) log cos（a1u）；1≤a1≤2

G2(u)=-(1/a2) exp(-a2u2/2)；a2≈1

對式(9)作進一步近似：

J(y)∝[E{G(y)}-E{G(v)}]2(10)

由y=wTX(y為其中一個獨立分量， w為分離矩陣W的一行，X為混合矩陣)，負熵的近似函數(shù)[16]也可表示為

G(W)∝[E{G(wTX)}-E{G(v)}]2(11)

問題轉(zhuǎn)換為求分離矩陣W，使分離出的估計信號y=wTX能使函數(shù)JG(W)達到最大。由Kuhn-Tucker 條件，轉(zhuǎn)換為無限制條件的優(yōu)化問題，得到目標(biāo)函數(shù)為

F(w)=E{G(wTX) }+c(‖w‖2-1)(12)

用牛頓迭代法對該目標(biāo)函數(shù)求解，得到FastICA 算法的迭代式：

w+=E{Xg(wTX)}-E{g′(wTX)}w，w*=w+/‖w+‖

2.2.2 算法的實現(xiàn)步驟

在計算機上實現(xiàn)上述FastICA算法，可以按照如下的迭代步驟進行：

a)令k=0，初始化權(quán)值向量w(0)。

b）k=k+1。

c）對w進行調(diào)整。

w(k+1)=E{Xg(wT(k)X)}-E{g′(wT(k)X)}w(k)

d)歸一化處理w(k+1)= w(k+1)/‖w(k+1)‖。

e）如果算法收斂，則求出一個獨立分量y=wX；否則轉(zhuǎn)入c）繼續(xù)迭代。

采用上述方法，依次迭代出W的權(quán)值向量w1T，w2T，…，wkT，構(gòu)成ICA算法的分離矩陣W的行向量。

3 FastICA算法在地震信號去噪中的應(yīng)用

3.1 去噪思想

地震勘探信號中包括有效信號s1和隨機干擾信號s2。如果s1與s2服從非高斯分布，且統(tǒng)計獨立，則可以按照ICA算法進行去噪[17，18]。圖2是ICA算法的去噪模型。其基本思想是：取一道記錄的兩次觀測結(jié)果，或近似取鄰近的兩道記錄（x1、x2），x1和x2都是s1和s2的線性組合；利用FastICA算法對觀測信號進行分離，得到y(tǒng)1和y2，它們分別為s1和s2的估計值。

3.2 仿真實驗

1）模擬地震信號(s1) 圖3(a)是由一個Ricker子波(頻率為50 Hz)和一均勻分布的隨機序列通過卷積生成的地震記錄。

2）模擬噪聲信號(s2) 圖3(b)是一均勻分布的隨機噪聲信號，均值為0.462 5，方差為0.247 2。圖中以橫軸為采樣點序號(n)，縱軸為波幅。

3）觀測信號(x1、x2) 分別如圖4(a)(b)所示，均為s1與s2的混合值。

利用FastICA算法進行分離，得到地震信號和噪聲信號。

4）分離出的地震信號(y1) 如圖5(a)所示。

5）分離出的噪聲信號(y2) 如圖5(b)所示。

將分離出的地震信號圖5(a)與原來的模擬地震信號圖3(a)進行對比，可以發(fā)現(xiàn)，除幅度、相位和正負號不同外，信號的其余特征非常相近，說明分離出的地震信號能夠有效反映原地震信號的特征，用該算法進行去噪是有效的。

4 結(jié)束語

隨著地震勘探的逐步深入，勘探區(qū)域已經(jīng)從平地延伸到山區(qū)、沙漠、灘海等區(qū)域，信號采集條件越來越惡劣，包含的干擾信息也越來越多，研究ICA算法在地震信號去噪中的應(yīng)用具有重要的現(xiàn)實意義。另外，本文采用的FastICA算法，不僅可以用于地震信號去噪，也可考慮用于重、磁等信號去噪，具有廣泛的應(yīng)用前景。

參考文獻：

［1］

COMON P.Independent component analysis，a new concept[J].Signal Processing，1994，36(3):287- 314.

[2]BELL A J，SEJNOWSKI T J.An information-maximization approach to blind separation and blind deconvolution[J].Neural Computation，1995，7(6):1004-1034.

[3]HYVARINEN A，OJA E.A fast fixed-point algorithm for independent component analysis [J].Neural Computation，1997，9(7):1483-1492.

[4]LEE T W，GIROLAMI M，SEJNOWSKI T J，et al.Independent component analysis using an extended infomax algorithm for mixed sub-Gaussian and super-Gaussian sourees[J].Neural Computation，1999，11(2):417-441.

［5］ HYVARINEN A，OJA E.Independent component analysis:algorithms and applications[J].Neural Networks，2000，13(4-5):411-430.

[6]百培瑞，牛海軍，師衛(wèi).基于盲分離的肺音信號中心音干擾的去除[J].太原理工大學(xué)學(xué)報，2002，33(3):329-331.

[7]夏文靜，傅行軍.基于ICA在強背景噪聲振動信號中的去噪研究[J].汽輪機技術(shù)，2006，48(2):121-123.

[8]鐘飛，譚中軍，史鐵林，等.基于ICA和小波變換的軸承故障特征提取[J].微計算機信息，2007，23(10):154-155.

[9]曾生根，夏德深.獨立分量分析在多光譜遙感圖像分類中的應(yīng)用[J].計算機工程與應(yīng)用，2004，40(21):108-110.

[10]胡天躍.地震資料疊前去噪技術(shù)的現(xiàn)狀與未來[J].地球物理學(xué)進展，2002，2(17):218-223.

[11]AL-YAHYA K M.Application of the partial Karhunen-Loeve transform to suppress random noise in seismic sections[J].Geophysical Prospecting，1993，39(1):77-93.

[12]DAUBECHIESL S W.Factoring wavelet transform into lifting steps[J].Journal of Fourier Analysis Application，1994，4(3):247-269.

[13]陳香朋，曹思遠.第二代小波變換及其在地震信號去噪中的應(yīng)用[J].石油物探，2004（6）:51-54.

[14]楊福生，洪波.獨立分量分析的原理與應(yīng)用[M].北京:清華大學(xué)出版社，2006:7-46.

[15]康永強，田夢君.一種快速ICA算法及其在腦電信號處理中的應(yīng)用[J].科學(xué)技術(shù)與工程，2006，6(3):529-533.

[16]SHI Zhen-wei，TANG Huan-wen，TANG Yi-yuan.A new fixed-point algorithm for independent component analysis[J].Neurocomputing，2004，56:467-473.

[17]劉喜武，劉洪，李幼銘.快速獨立分量變換與去噪初探[J].中國科學(xué)院研究生院學(xué)報，2003，20(4):488-492.

[18]焦衛(wèi)東，楊世錫，吳昭同.基于獨立分量分析的噪聲消除技術(shù)研究[J].浙江大學(xué)學(xué)報:工學(xué)版，2004，38(7):872-876.