基于FastICA 的工頻干擾消除算法

張君生，劉瑞安，王 磊

（天津師范大學(xué)物理與電子信息學(xué)院，天津 300387）

獨(dú)立成分分析以非高斯源信號(hào)為研究對(duì)象，在統(tǒng)計(jì)獨(dú)立的假設(shè)下，對(duì)多路觀測(cè)到的混合信號(hào)進(jìn)行盲分離，從而較完好地分離出隱含在混合信號(hào)中的源信號(hào).近年來，通過ICA 實(shí)現(xiàn)盲信源分離已經(jīng)逐漸成為信號(hào)處理領(lǐng)域的熱點(diǎn)問題之一.從1990 年代以來，研究人員從不同角度提出多種ICA 算法.其中，Nadal 等[1]提出了信息最大化（Infomax）基本原理，Bell 等[2]給出該方法的隨機(jī)梯度學(xué)習(xí)規(guī)則，并將其分別運(yùn)用于預(yù)測(cè)、時(shí)間序列分析和盲源分離等研究中.Hyv?rinen 等[3]提出并詳細(xì)論述了最大似然估計(jì)法（MLE）.

應(yīng)用過程中，ICA 算法所假設(shè)的條件，如信號(hào)之間的獨(dú)立性、線性、瞬時(shí)性以及不考慮混合時(shí)產(chǎn)生的噪聲等條件不一定能夠得到滿足，導(dǎo)致算法的應(yīng)用效果并不理想.為解決這些問題，很多改進(jìn)的ICA 方法應(yīng)運(yùn)而生，包括非線性ICA 法、約束ICA法和多分辨率子波段分解ICA 法等[4-9].除了對(duì)基本ICA 方法的改進(jìn)外，基于信號(hào)時(shí)間結(jié)構(gòu)和稀疏性的ICA 新算法不斷被提出，奇異混合[10-11]、時(shí)延混合[12]和噪聲ICA[13]等問題也吸引了眾多學(xué)者的關(guān)注.

實(shí)際應(yīng)用中，有用信號(hào)與噪聲信號(hào)往往是獨(dú)立的，本研究將較弱的有效地震信號(hào)和較強(qiáng)的工頻干擾看成2 個(gè)獨(dú)立的信源，接收到的信號(hào)是兩者混合，采用FastICA 算法消除地震信號(hào)中的工頻干擾.

1 ICA 算法

ICA 的基本模型[14]可以具體描述為：

式（1）中：x=（x1，x2，…，xn）T為n 維隨機(jī)觀測(cè)向量；s=（s1，s2，…，sm）T是m 維未知源信號(hào)，si稱為獨(dú)立分量；A 為滿秩的混合矩陣（為了簡化模型，假設(shè)A 為m=n 的方陣）；ai是混合矩陣的基向量.

為了尋找變換矩陣或稱解混矩陣W，對(duì)x 進(jìn)行線性變換，得到n 維輸出向量

在比例允許和順序不確定的前提下，y=s?是對(duì)si的一個(gè)估計(jì)量.以2 個(gè)源信號(hào)為例，ICA 線性組合模型的基本運(yùn)算流程如圖1 所示.

2 FastICA 算法

FastICA 算法又稱Fixed-Point 算法，是一種通過牛頓迭代法使目標(biāo)函數(shù)最優(yōu)化的迭代算法，具有收斂速度快、穩(wěn)健性高等優(yōu)點(diǎn).與梯度算法相比，F(xiàn)astICA 算法不要求學(xué)習(xí)步長，可通過確定數(shù)據(jù)塊的學(xué)習(xí)方式調(diào)整分離矩陣系數(shù).為減少算法需要估計(jì)的參數(shù)，簡化算法的計(jì)算需要對(duì)數(shù)據(jù)進(jìn)行預(yù)處理.FastICA 算法的預(yù)處理步驟是對(duì)觀測(cè)信號(hào)進(jìn)行白化處理，使白化后的分量x 為非相關(guān)的，且滿足E{xxT}=I.

在非高斯性最大化方法中，F(xiàn)astICA 算法基于負(fù)熵來進(jìn)行估計(jì).為求出y=wTx 的投影方向，首先利用式（3）進(jìn)行最大化

式（3）中：w 是m 維變量.在預(yù)處理后，式（3）的最大化問題可以轉(zhuǎn)化為E[G（wTx）]的優(yōu)化問題.在約束條件E[（wTx）2]的限制下，E[G（wTx）]的優(yōu)化可以通過式（4）獲得，

式（4）中：β 是常量，可通過β=E[w0Txg（w0Tx）]獲得；w0是w 的初始值，函數(shù)g（u）為

式（5）中：1≤a1≤2；a2≈1.利用牛頓迭代法解式（4），式（4）左邊部分記為F（w），其雅可比（Jacobian）矩陣

為了簡化求矩陣的逆，利用式（7）對(duì)式（6）中的E[xxTg′（wTx）]進(jìn)行近似，

自媒體又稱“個(gè)人媒體”或“公民媒體”，是普通大眾以現(xiàn)代化和電子化的手段，傳遞規(guī)范性及非規(guī)范性信息，提供與分享他們自身的事實(shí)、新聞的途徑。本研究通過微信和新浪微博等自媒體平臺(tái)，在搜索欄輸入“黃山民宿”，剔除純圖片，純營銷，純介紹性文字，得到最終有效游記和評(píng)價(jià)共198條。

亦可以得到式（7）的牛頓迭代法的近似解

為了提高算法的穩(wěn)定性，在迭代后利用wk+1=wk+1/‖wk+1‖歸一化w.在式（8）兩邊分別乘以β=E[g′（wTx）]，得到固定點(diǎn)算法的迭代公式

基于式（9），F(xiàn)astICA 算法即固定點(diǎn)迭代算法的基本結(jié)構(gòu)為：

（1）權(quán)向量w 的初始化選擇，k=0.

（2）對(duì)wk+1進(jìn)行調(diào)整.

3 實(shí)驗(yàn)結(jié)果

地震記錄中的有效信號(hào)可以看作是地震子波與地層反射系數(shù)序列的卷積，假設(shè)地層反射系數(shù)是非高斯分布的，即地震信號(hào)也非高斯分布，由于有效地震信號(hào)和50 Hz 工頻干擾分別由不同的信源產(chǎn)生，彼此相互獨(dú)立，因此地震信號(hào)和工頻干擾信號(hào)滿足ICA 的應(yīng)用限制條件[16].

為了滿足FastICA 的應(yīng)用條件，構(gòu)造觀測(cè)信號(hào)x2（t）=a21s（t）+a22n（t）.實(shí)際情況下，s（t）未知，工頻干擾信號(hào)滿足n（t）=A0sin（2πf0t+θ0）.在不考慮幅值的情況下，構(gòu)造觀測(cè)信號(hào)x2（t）=0×s（t）+a22（sin 2πf0t+θ0）以符合ICA 的要求.觀測(cè)信號(hào)可進(jìn)一步表示為

式（10）中：b=a12A0cos θ0，c=a12A0sin θ0.可以通過構(gòu)造2 組參考信號(hào)s2=sin（2πf0t）和s3=cos（2πf0t）與含有工頻干擾和地震信號(hào)的混合信號(hào)共同作為ICA 的輸入信號(hào)，分離出有效地震信號(hào)s（t），這樣就避免了對(duì)θ0的直接估計(jì)[17]，其數(shù)學(xué)模型表示為

利用FastICA 算法對(duì)式（11）中的混合矩陣A 和分離矩陣W=A-1進(jìn)行估計(jì)，并利用式（12）求出不含工頻干擾的地震信號(hào)s（t）.

式（12）中：X 和S 分別為觀測(cè)信號(hào)向量和源信號(hào)向量；A 為混合矩陣.

在地震信號(hào)中選取3 000 個(gè)采模點(diǎn)進(jìn)行降噪實(shí)驗(yàn)，F(xiàn)astICA 算法對(duì)混有工頻干擾的地震信號(hào)的干擾消除結(jié)果如圖2 所示.

圖2c 為混入50 Hz 工頻干擾的地震源信號(hào)，圖2d 和圖2e 為依據(jù)前文介紹算法構(gòu)造出的兩路參考信號(hào).圖2f～圖2h 為經(jīng)過FastICA 算法處理后三路信號(hào)（圖2c～圖2e）的恢復(fù)信號(hào)，可以看出，圖2g 所示的第2 路恢復(fù)信號(hào)為所需要的信號(hào).

信號(hào)的相關(guān)系數(shù)描述了信號(hào)之間的線性相關(guān)性，表1 給出經(jīng)過20 次實(shí)驗(yàn)所得的分離信號(hào)和源信號(hào)的相關(guān)系數(shù)絕對(duì)值.由表1 中數(shù)據(jù)可以看出，分離信號(hào)與源信號(hào)的相關(guān)系數(shù)絕對(duì)值均高于0.998 8，說明分離效果較為理想.

表1 分離信號(hào)與源信號(hào)的相關(guān)系數(shù)絕對(duì)值Tab.1 Absolute values of correlation of separative signal and source signal

4 對(duì)比實(shí)驗(yàn)

4.1 不同方法的比較

采用有限長單位沖激響應(yīng)（Finite Impulse Response，F(xiàn)IR）陷波濾波器和最小均方誤差（Least Mean Square，LMS）自適應(yīng)濾波器2 種傳統(tǒng)消除工頻干擾的方法與FastICA 算法進(jìn)行比較.實(shí)驗(yàn)中，陷波器的采樣頻率為fs=1 000 Hz，通過將邊界頻率轉(zhuǎn)換為歸一化頻率可減少其帶寬；LMS 自適應(yīng)濾波器的步長因子為0.001，自適應(yīng)濾波器可以自動(dòng)調(diào)整濾波器系數(shù)，以跟蹤輸入過程的變化，并實(shí)現(xiàn)工頻干擾的自適應(yīng)抵消.

FIR 陷波濾波器對(duì)混有工頻干擾的地震信號(hào)濾波的結(jié)果如圖3 所示.由圖3 可以看出，輸出信號(hào)（圖3c）較輸入信號(hào)（圖3a）存在明顯延遲，嚴(yán)重影響了地震信號(hào)分離的時(shí)效性.

利用LMS 自適應(yīng)濾波器消除地震信號(hào)中工頻干擾的結(jié)果如圖4 所示.加入正弦干擾（即工頻干擾）的原始信號(hào)（圖4b）在經(jīng)過LMS 自適應(yīng)濾波器的降噪處理后，得到的輸出信號(hào)如圖4c 所示.圖4d為誤差信號(hào)，對(duì)于地震信號(hào)而言，此結(jié)果的信號(hào)誤差較大.

表2 給出FIR 陷波濾波器、LMS 自適應(yīng)濾波器和FastICA 算法輸出結(jié)果相關(guān)系數(shù)絕對(duì)值的比較結(jié)果.由表2可以明顯看出，采用FastICA 算法得到的恢復(fù)信號(hào)與源信號(hào)的相關(guān)系數(shù)要高于其他2 種典型的傳統(tǒng)濾波方法，這說明FastICA 算法消除地震信號(hào)中的工頻干擾能力優(yōu)于傳統(tǒng)方法.

表2 3 種方法相關(guān)系數(shù)絕對(duì)值的對(duì)比Tab.2 Comparison of absolute values of correlation of 3 methods

4.2 采樣點(diǎn)對(duì)結(jié)果的影響

考慮到ICA 本身的統(tǒng)計(jì)特性，減少各信號(hào)的采樣點(diǎn)至1 000 個(gè)時(shí)，實(shí)驗(yàn)結(jié)果如圖5 所示.

由圖5 可知，ICA 算法的準(zhǔn)確性取決于采樣點(diǎn)的數(shù)量.采樣點(diǎn)越少，即先驗(yàn)知識(shí)越少，恢復(fù)信號(hào)與原信號(hào)的相關(guān)系數(shù)就越??；采樣點(diǎn)越多，得到的相關(guān)系數(shù)越大，但隨之而來的算法的時(shí)間度亦將變大，從而給算法的時(shí)效性帶來了影響.

表3 給出1 000 個(gè)采樣點(diǎn)時(shí)FastICA 算法消除工頻干擾后分離信號(hào)與源信號(hào)相關(guān)系數(shù)的絕對(duì)值.同樣進(jìn)行20 次實(shí)驗(yàn)，1 000 個(gè)采樣點(diǎn)相關(guān)系數(shù)絕對(duì)值的均值為0.971 8，顯然小于3 000 個(gè)采樣點(diǎn)的均值0.998 8，說明采樣點(diǎn)減少后，恢復(fù)出來的信號(hào)失真率較高.進(jìn)一步說明采樣點(diǎn)對(duì)FastICA 算法的影響，也說明ICA 算法這種基于信號(hào)統(tǒng)計(jì)特性的算法存在的不足之處.

表3 選取1 000 個(gè)采樣點(diǎn)時(shí)，分離信號(hào)與源信號(hào)的相關(guān)系數(shù)絕對(duì)值Tab.3 Results of absolute values of correlation of separative signal and source signal in 1 000 sampling points

5 結(jié)論

根據(jù)ICA 的數(shù)學(xué)模型，分析了不同ICA 算法以及它們之間的聯(lián)系.將FastICA 算法應(yīng)用于消除地震信號(hào)中的工頻干擾，并對(duì)輸出信號(hào)的相關(guān)系數(shù)絕對(duì)值進(jìn)行對(duì)比.結(jié)果表明：經(jīng)FastICA 算法處理后，信號(hào)的相關(guān)系數(shù)絕對(duì)值優(yōu)于傳統(tǒng)方法，說明在消除微信號(hào)工頻干擾的問題中，與傳統(tǒng)工頻干擾消除技術(shù)相比，F(xiàn)astICA 算法可以更準(zhǔn)確地從混合信號(hào)中分離出實(shí)際所需要的源信號(hào)，證明FastICA 算法在消除地震信號(hào)中的工頻干擾方面具有良好的應(yīng)用前景.

[1]NADAL J，PARGA N.Nonlinear neurons in low-noise limit：A factorial code maximizes information transfer[J].Network，1994，5（3）：565—581.

[2]BELL A，SEJNOWSKI T.An information maximization approach to blind separation and blind de-convolution[J].Neural Computation，1995，7（6）：1129—1159.

[3]HYV?RINEN A，OJA E.Independent component analysis：Algorithms and applications[J].Neural Networks，2000，13（4/5）：411—430.

[4]HYV?RINEN A，PAJUNEN P.Nonlinear independent component analysis：Existence and uniqueness results[J].Neural Networks，1999，12（3）：429—439.

[5]JUTTEN C，ZADEH M，HOSSEINI S.Three easy ways for separating nonlinear mixtures[J].Signal Processing，2004，84（2）：217—229.

[6]LIU W，RAJAPAKSE J.Approach and applications of constrained ICA[J].Trans Neural Networks，2005，16（1）：203—212.

[7]PLUMBLEY M.Algorithms for nonnegative component independent analysis[J].Trans Neural Networks，2003，14（3）：534—543.

[8]王世海，陳向東，畢雪，等.多分辨率子帶分解的獨(dú)立分量分析算法在紅外圖像去噪上的應(yīng)用[J].電子技術(shù)應(yīng)用，2007（6）：66—68.

[9]LEE T，LEWICKI M，GIROLAMI M，et al.Blind source separation of more sources than mixtures using over-complete representations[J].Signal Processing Letters，1999，4（5）：205—208.

[10]白琳，陳豪.一種奇異混合陣的盲信號(hào)提取算法[J].空間電子技術(shù)，2012（1）：6—10，22.

[11]李遠(yuǎn)清，張麗清，劉永清.一種奇異混合信號(hào)盲分離的神經(jīng)網(wǎng)絡(luò)模型[J].控制理論與應(yīng)用，2001，18（1）：105—108.

[12]肖俊，何為偉.一種基于信息最大化的時(shí)延混合盲分離算法[J].井岡山師范學(xué)院學(xué)報(bào)：自然科學(xué)版，2005，26（3）：12—14.

[13]HYV?RINEN A，KARHUNEN J，OJA C.Independent Component Analysis[M].New York：Wiley，2001：293—304.

[14]楊福生，洪波，唐慶玉.獨(dú)立分量分析及其在生物醫(yī)學(xué)工程中的應(yīng)用[J].國外醫(yī)學(xué)：生物醫(yī)學(xué)工程分冊(cè)，2000，23（3）：129—134.

[15]郭松.獨(dú)立分量分析及其應(yīng)用研究[D].武漢：武漢大學(xué)，2004：13—20.

[16]劉喜武，劉洪，李幼銘.獨(dú)立分量分析及其在地震信息處理中的應(yīng)用初探[J].地球物理學(xué)進(jìn)展，2003，18（1）：90—96.

[17]吳小培，詹長安，周荷琴，等.采用獨(dú)立分量分析的方法消除信號(hào)中的工頻干擾[J].中國科技大學(xué)學(xué)報(bào)，2000，30（6）：671—676.