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

張君生，劉瑞安，王 磊

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

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

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

實際應用中，有用信號與噪聲信號往往是獨立的，本研究將較弱的有效地震信號和較強的工頻干擾看成2 個獨立的信源，接收到的信號是兩者混合，采用FastICA 算法消除地震信號中的工頻干擾.

1 ICA 算法

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

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

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

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

2 FastICA 算法

FastICA 算法又稱Fixed-Point 算法，是一種通過牛頓迭代法使目標函數最優化的迭代算法，具有收斂速度快、穩健性高等優點.與梯度算法相比，FastICA 算法不要求學習步長，可通過確定數據塊的學習方式調整分離矩陣系數.為減少算法需要估計的參數，簡化算法的計算需要對數據進行預處理.FastICA 算法的預處理步驟是對觀測信號進行白化處理，使白化后的分量x 為非相關的，且滿足E{xxT}=I.

在非高斯性最大化方法中，FastICA 算法基于負熵來進行估計.為求出y=wTx 的投影方向，首先利用式（3）進行最大化

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

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

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

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

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

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

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

基于式（9），FastICA 算法即固定點迭代算法的基本結構為：

（1）權向量w 的初始化選擇，k=0.

（2）對wk+1進行調整.

3 實驗結果

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

為了滿足FastICA 的應用條件，構造觀測信號x2（t）=a21s（t）+a22n（t）.實際情況下，s（t）未知，工頻干擾信號滿足n（t）=A0sin（2πf0t+θ0）.在不考慮幅值的情況下，構造觀測信號x2（t）=0×s（t）+a22（sin 2πf0t+θ0）以符合ICA 的要求.觀測信號可進一步表示為

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

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

式（12）中：X 和S 分別為觀測信號向量和源信號向量；A 為混合矩陣.

在地震信號中選取3 000 個采模點進行降噪實驗，FastICA 算法對混有工頻干擾的地震信號的干擾消除結果如圖2 所示.

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

信號的相關系數描述了信號之間的線性相關性，表1 給出經過20 次實驗所得的分離信號和源信號的相關系數絕對值.由表1 中數據可以看出，分離信號與源信號的相關系數絕對值均高于0.998 8，說明分離效果較為理想.

表1 分離信號與源信號的相關系數絕對值Tab.1 Absolute values of correlation of separative signal and source signal

4 對比實驗

4.1 不同方法的比較

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

FIR 陷波濾波器對混有工頻干擾的地震信號濾波的結果如圖3 所示.由圖3 可以看出，輸出信號（圖3c）較輸入信號（圖3a）存在明顯延遲，嚴重影響了地震信號分離的時效性.

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

表2 給出FIR 陷波濾波器、LMS 自適應濾波器和FastICA 算法輸出結果相關系數絕對值的比較結果.由表2可以明顯看出，采用FastICA 算法得到的恢復信號與源信號的相關系數要高于其他2 種典型的傳統濾波方法，這說明FastICA 算法消除地震信號中的工頻干擾能力優于傳統方法.

表2 3 種方法相關系數絕對值的對比Tab.2 Comparison of absolute values of correlation of 3 methods

4.2 采樣點對結果的影響

考慮到ICA 本身的統計特性，減少各信號的采樣點至1 000 個時，實驗結果如圖5 所示.

由圖5 可知，ICA 算法的準確性取決于采樣點的數量.采樣點越少，即先驗知識越少，恢復信號與原信號的相關系數就越??；采樣點越多，得到的相關系數越大，但隨之而來的算法的時間度亦將變大，從而給算法的時效性帶來了影響.

表3 給出1 000 個采樣點時FastICA 算法消除工頻干擾后分離信號與源信號相關系數的絕對值.同樣進行20 次實驗，1 000 個采樣點相關系數絕對值的均值為0.971 8，顯然小于3 000 個采樣點的均值0.998 8，說明采樣點減少后，恢復出來的信號失真率較高.進一步說明采樣點對FastICA 算法的影響，也說明ICA 算法這種基于信號統計特性的算法存在的不足之處.

表3 選取1 000 個采樣點時，分離信號與源信號的相關系數絕對值Tab.3 Results of absolute values of correlation of separative signal and source signal in 1 000 sampling points

5 結論

根據ICA 的數學模型，分析了不同ICA 算法以及它們之間的聯系.將FastICA 算法應用于消除地震信號中的工頻干擾，并對輸出信號的相關系數絕對值進行對比.結果表明：經FastICA 算法處理后，信號的相關系數絕對值優于傳統方法，說明在消除微信號工頻干擾的問題中，與傳統工頻干擾消除技術相比，FastICA 算法可以更準確地從混合信號中分離出實際所需要的源信號，證明FastICA 算法在消除地震信號中的工頻干擾方面具有良好的應用前景.

[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]王世海，陳向東，畢雪，等.多分辨率子帶分解的獨立分量分析算法在紅外圖像去噪上的應用[J].電子技術應用，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]白琳，陳豪.一種奇異混合陣的盲信號提取算法[J].空間電子技術，2012（1）：6—10，22.

[11]李遠清，張麗清，劉永清.一種奇異混合信號盲分離的神經網絡模型[J].控制理論與應用，2001，18（1）：105—108.

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

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

[14]楊福生，洪波，唐慶玉.獨立分量分析及其在生物醫學工程中的應用[J].國外醫學：生物醫學工程分冊，2000，23（3）：129—134.

[15]郭松.獨立分量分析及其應用研究[D].武漢：武漢大學，2004：13—20.

[16]劉喜武，劉洪，李幼銘.獨立分量分析及其在地震信息處理中的應用初探[J].地球物理學進展，2003，18（1）：90—96.

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