ＦａｓｔＩＣＡ算法及其在地震信號(hào)去噪中的應(yīng)用

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

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

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

中圖分類(lèi)號(hào)：TP301.6文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：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 引言

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

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

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

1 ICA方法的基本原理

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

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

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

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

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

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

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

a)源信號(hào)S(t)各分量均為平穩(wěn)隨機(jī)過(guò)程，且相互獨(dú)立。

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

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

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

1.3 模型解的不確定性

1.3.1 幅度的不確定性

對(duì)式(1)，可以改寫(xiě)為

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

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

1.3.2 輸出順序的不確定性

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

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

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

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

2 基于負(fù)熵極大的FastICA算法

2.1 FastICA算法

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

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

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

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

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

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

預(yù)處理結(jié)束，就可以進(jìn)行非高斯最大化計(jì)算(即分離過(guò)程)，此時(shí)必須先確定非高斯性的度量標(biāo)準(zhǔn)。非高斯性的度量標(biāo)準(zhǔn)主要有峭度(kurtosis)和負(fù)熵(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ù)尋找獨(dú)立分量時(shí)，易受大幅度隨機(jī)脈沖干擾的影響，而負(fù)熵則顯得更為穩(wěn)健。

2.2 基于負(fù)熵極大的FastICA算法

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

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

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

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

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

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

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

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

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

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

對(duì)式(9)作進(jìn)一步近似：

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

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

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

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

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

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

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

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

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

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

b）k=k+1。

c）對(duì)w進(jìn)行調(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）如果算法收斂，則求出一個(gè)獨(dú)立分量y=wX；否則轉(zhuǎn)入c）繼續(xù)迭代。

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

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

3.1 去噪思想

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

3.2 仿真實(shí)驗(yàn)

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

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

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

利用FastICA算法進(jìn)行分離，得到地震信號(hào)和噪聲信號(hào)。

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

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

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

4 結(jié)束語(yǔ)

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

參考文獻(xiàn)：

［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).基于盲分離的肺音信號(hào)中心音干擾的去除[J].太原理工大學(xué)學(xué)報(bào)，2002，33(3):329-331.

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

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

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

[10]胡天躍.地震資料疊前去噪技術(shù)的現(xiàn)狀與未來(lái)[J].地球物理學(xué)進(jìn)展，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]陳香朋，曹思遠(yuǎn).第二代小波變換及其在地震信號(hào)去噪中的應(yīng)用[J].石油物探，2004（6）:51-54.

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

[15]康永強(qiáng)，田夢(mèng)君.一種快速I(mǎi)CA算法及其在腦電信號(hào)處理中的應(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]劉喜武，劉洪，李幼銘.快速獨(dú)立分量變換與去噪初探[J].中國(guó)科學(xué)院研究生院學(xué)報(bào)，2003，20(4):488-492.

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