曹小玲 劉開元 嚴良俊

(①長江大學油氣資源與勘探技術教育部重點實驗室，湖北武漢 430100； ②非常規(guī)油氣湖北省協(xié)同創(chuàng)新中心，湖北武漢 430100；③長江大學信息與數學學院，湖北荊州 434023； ④中國石油川慶鉆探工程有限公司地球物理勘探公司，四川成都 610213)

1 引言

大地電磁測深(MT)是一種被動源電磁探測技術，通過在地表測量正交電、磁場分量的擾動值，研究地下介質的地電結構[1]。MT探測深度的范圍從地表到上地幔，它既可以研究大尺度地質單元的電磁感應、深部地殼地幔電性結構及殼幔中的高導層，也可以開展油氣資源普查、地熱資源勘查、天然地震活動前兆監(jiān)測以及地下水資源勘查和環(huán)境保護等，因而得到廣泛應用[2]。但大地電磁測深觀測的天然電磁場具有信號弱、頻帶寬等特點，極易受到各種噪聲的干擾，因此大地電磁數據的去噪一直是研究熱點。地球物理學家們提出一系列技術，這些技術主要集中在四個方面：遠參考方法、穩(wěn)健阻抗估計、時頻轉換技術和時間域去噪[3-7]。

傳統(tǒng)的數據處理基于傅里葉變換—功率譜估計，利用傅里葉變換將信號通過一個低通或帶通濾波器，映射到整個頻域內加以分析處理，但無法表達信號的時域局部性質，而這恰恰是大地電磁信號作為一種非平穩(wěn)信號的最基本和最關鍵的特征。為此，小波分析被陸續(xù)應用于大地電磁的去噪處理[8-13]。王書明等[14,15]根據MT信號的雙相干系數、三相干系數以及在z平面上的零點分布，分析研究不同地區(qū)實測MT信號的特征，得出MT信號一般具有非高斯性、非線性和非最小相位特征,這個結論啟發(fā)人們用新的信號處理方法解決大地電磁的噪聲問題。近年來，隨著盲源分離技術的迅猛發(fā)展，獨立分量分析理論[16-19]被廣泛應用于生物醫(yī)學信號處理、語音信號處理、圖像處理、地震勘探等多個領域，在信號去噪中表現優(yōu)異。獨立分量分析要求源信號在統(tǒng)計特性上具有非高斯性、非依賴性、非白性[16-19]，大地電磁信號恰好符合這些要求；另一方面，獨立分量分析不需要知道信號的先驗知識，即它只要通過觀測信號就能推測甚至恢復原始信號，這對大地電磁測深來說意義重大，因為由于各種客觀因素的制約，導致在大地電磁測深中想要獲得相關的先驗知識一般來說是非常困難的。因此，使用獨立分量分析技術對大地電磁信號的各種信息成分進行分析和提取，尤其是對噪聲成分進行分析和提取，符合大地電磁測深實際應用的特點和需要。獨立分量分析處理的對象是一組相互統(tǒng)計獨立的信號源經線性組合而產生的混合信號,最終從混合信號中提取各獨立的信號分量。就大地電磁噪聲干擾問題,目標信號和噪聲干擾分別由不同的信號源產生，因此彼此相互獨立。由統(tǒng)計原理，相互獨立的隨機變量必然不相關； 但反過來卻不一定成立[20]，即不相關的隨機變量不一定相互獨立。獨立分量分析方法在統(tǒng)計獨立意義下對混合信號進行分離，以達到消除噪聲的目的，這顯然比基于不相關性度量的傳統(tǒng)自適應濾波方法的效果更好。因此，根據大地電磁測深數據的特點和其噪聲的來源，以及大地電磁資料的處理流程和解釋流程，研究基于獨立分量分析的大地電磁數據噪聲消除方法，并通過實際大地電磁數據驗證方法的有效性。

2 獨立分量分析

獨立分量分析(ICA)[16-19]是20世紀90年代發(fā)展起來的一種信號處理技術，它是從多維統(tǒng)計數據中找出隱含因子或分量的方法。從線性變換和線性空間角度來說，源信號為相互獨立的非高斯信號，可以看作線性空間的基信號，而觀測信號則為源信號的線性組合。ICA就是在源信號和線性變換均未知的情況下，從觀測的混合信號中估計出數據空間的基本結構或源信號。ICA的算法流程如圖1所示：在源信號s(t)中各分量相互獨立的假設下，由觀察信號x(t)通過解混系統(tǒng)B把他們分離開來，使輸出信號y(t)逼近s(t)。

圖1 ICA的算法流程圖

2.1 數據預處理

ICA通常需要進行零均值化和白化預處理。

一般情況下，觀測數據是具有相關性的，因此需要對零均值化后的數據進行白化處理，以去除相關性、簡化后續(xù)分量的提取過程，數據的白化處理還能有效增強算法的收斂性。

若均值為零的隨機向量Z=(Z1,…,ZM)T滿足E(ZZT)=I(I為單位矩陣，M為隨機向量中元素的個數，E(·)是期望值算子)，則稱Z為白化向量。白化的本質是去相關。對觀測信號X，需尋找一個變換，使X經變換后變成白化向量，即有：Z=QX,其中Q為白化矩陣。

零均值化和白化作為ICA的預處理可以有效地降低問題的復雜度。

2.2 FastICA算法

FastICA算法又稱固定點(Fixed-Point)算法，它采用批處理方式，速度比較快,是一種快速尋優(yōu)迭代算法。FastICA算法有多種形式，本文采用基于負熵最大的FastICA算法。它以負熵最大作為搜尋方向，可以實現順序提取獨立源；此外，還采用了定點迭代的優(yōu)化算法，使收斂更加快速、穩(wěn)健。

MT時間序列信號可視為隨機變量在某個時刻的采樣，因此可以將時間序列信號設為隨機變量。假定這里有m個觀測信號，第i個觀測信號xi是由n個相互獨立的未知源信號s1、s2、…、sn混合而成(xi、sj均為隨機變量)，即

xi=ai，1s1+ai，2s2+…+ai，nsn

(1)

上式中ai,j(j=1,2,…,n)為常數系數。用矢量X表示觀測信號變量(x1,x2,…,xm)T，矢量S表示源信號變量(s1,s2,…,sn)T(m≥n)，Am×n表示混合矩陣(ai,j)，則式(1)可表示為

X=AS

(2)

FastICA的目標就是尋求分離矩陣Wm×n，以使Y=WTX=WTAS具有最大的非高斯性。

設離散隨機變量y的熵定義為

(3)

式中：ki指y的可能取值；P(·)指求概率。式中對數一般取2為底，單位為比特。其負熵定義為

Ng(y)=H(yGauss)-H(y)

(4)

式中yGauss是一個與y具有相同方差的高斯隨機變量。負熵Ng(y)可以作為隨機變量y非高斯性的測度。由于計算離散隨機變量的熵需要知道隨機變量的概率分布，其不易獲取，因此常采用下式計算負熵

Ng(y)={E[g(y)]-E[g(yGauss)]}2

(5)

式中g(·)為特定非線性函數。

FastICA算法原理如下。

首先，通過對E[g(WTX)]進行優(yōu)化獲得y的負熵的最大近似值。根據Kuhn-Tucker條件，在E[(WTX)2]=‖W‖2=1的約束下，E[g(WTX)]的最優(yōu)值在滿足下式的點上獲得

E[Xg(WTX)]+βW=0

(6)

然后，利用牛頓迭代法求解上述方程。記P=E[Xg(WTX)]+βW，可得P的雅可比矩陣為

JP=E[XXTg′(WTX)]-βI

(7)

式中g′表示g的一次導數。由于數據被白化，有E(XXT)=I,故

E[XXTg′(WTX)] ≈E(XXT)·E[g′(WTX)]

=E[g′(WTX)]I

(8)

因而JP變成了對角陣，容易求逆。如此便得到近似牛頓迭代公式

(9)

式中：W*是W進行一次迭代后的取值；β=E[WTXg(WTX)]。該式簡化后就得到FastICA算法的迭代公式

(10)

MT觀測的是電場和磁場的時間序列，把獲得的時間序列看作一個隨機變量在特定時候的采樣，即將記錄的電場或磁場時間序列視為隨機變量x1，x2，…，xm,則可以利用上面所述的FastICA算法進行迭代處理，便可在獲得分離矩陣W后獲得具有最大的非高斯性的隨機變量矢量

Y=WTX=WTAS

(11)

Y既包含原始信號也包含噪聲，但因為可以將大地電磁資料中的原始電磁信號和噪聲看作是相互獨立的信號源，因此可以利用FastICA算法進行原始信號與噪聲的分離和提取。

2.3 基于小波分析的獨立分量分析

經過多次試驗后發(fā)現，小波變換對高信噪比的含噪信號去噪效果較好，對低信噪比的含噪信號去噪效果較差，而屬于盲源分離的獨立分量分析算法對低信噪比的含噪信號的去噪效果卻較好。

根據盲源分離理論[21]，在不知道源信號及混合矩陣的任何信息的情況下，只需假設源信號是相互統(tǒng)計獨立的，就能使用獨立分量分析將源信號從混合信號中分離出來。在大地電磁勘探的實際應用中，并不能通過有效方法證明源信號之間是嚴格獨立的，更不能證明源信號是嚴格線性瞬時混合的，但從各種源信號的不同存在方式來看，它們之間至少是獨立的，也就是說它們滿足相互統(tǒng)計獨立的要求。從另一方面來說，由于大地電磁有效信號與噪聲信號的來源不一致，因而它們的自身屬性也不相同，所以從屬性角度看它們也是相互獨立的。因此，可以把大地電磁資料中的有效信號與噪聲信號視為相互統(tǒng)計獨立的信號，這樣就可以使用獨立分量分析算法來對它們進行分離。

但獨立分量分析也有其弊端，那就是它的適用條件要求采集到的信號的數目要大于或者等于源信號的數目[16-19]，這就意味著必須同時獲得至少兩組在同時同地采集到的數據才能將所得數據分離為有效信號和噪聲信號兩組數據。而且，一般而言，使用的采集信號的組數越多，提供的噪聲特征的信息量越豐富，因此在許多情況下要求使用在同一時間、同一地點采集到的至少三組大地電磁信號來進行實驗和研究。而在實際野外作業(yè)中，由于地質條件、氣候條件、環(huán)境條件及人力物力條件的制約，要獲得這樣的三組及以上的數據是非常困難的，且成本非常高。受陳艷等[22]研究的啟發(fā)，首先對信號進行小波分解和去噪處理，再應用獨立分量分析技術進行盲源分離，以提取源信號和噪聲。為了減小觀測信號的信噪比對去噪效果的影響，本文引入動態(tài)的自適應因子?，根據觀測信號的幅值情況不同而取不同的?值；另外，由于噪聲一般包含在小波高頻分量中，所以重構信號由不含噪聲的低頻分量(近似信號)獲得。根據如下公式重構信號

Y=R+?*T

(12)

式中：R為小波近似信號；T為ICA分離后提取的有效源信號；Y為去噪后信號。

基于小波分析的ICA去噪的算法流程如圖2所示。

圖2 基于小波分析的ICA去噪算法流程

3 仿真試驗

為研究獨立分量分析算法對信號處理的效果，首先進行仿真試驗。分別采用周期信號和非周期信號進行實驗研究。

以t表示時間變量，周期信號取y1=20sin150t+5cos100t，非周期信號取leleccum信號，即y2=leleccum(t)，這里的leleccum信號為matlab所提供的一種非周期非規(guī)則信號。

分別運用七種方法進行去噪： ①啟發(fā)式軟閾值去噪法(5層分解)； ②啟發(fā)式硬閾值去噪法(5層分解)； ③SURE軟閾值去噪法(5層分解)； ④SURE硬閾值去噪法(5層分解)； ⑤自適應默認閾值去噪法(5層分解)； ⑥重構信號Y=R+∑T′去噪法(T′為ICA分離后未經處理的原始結果信號)(5層分解)； ⑦由式(12)所示的去噪法(5層分解)。得到的實驗結果見圖3。限于篇幅，這里只顯示加入噪聲(SNR=10)時y1的實驗結果。

圖3 原始信號和去噪結果 (a)原始信號和含噪信號； (b)含噪信號5層小波分解后的結果； (c)ICA處理后的結果； (d)七種方法分別處理后的結果

為了比較去噪效果，這里采用Pearson(皮爾遜)相關系數ρY,O度量去噪后信號Y與原始信號O(無噪信號)之間的相關度

(13)

顯然Y、O包含的元素個數相等。

表1和表2是分別運用上述七種去噪方法對信號y1和y2進行去噪處理后的結果。可以看出，無論哪一種信號、無論信噪比大小，本文所采用的去噪方法的去噪效果都比較穩(wěn)定。所采用的恢復信號既包含低頻小波系數，又含有獨立分量分析分離出的源信號，并且通過自適應因子動態(tài)改變獨立分量分析算法的分離效果，使觀測信號的信噪比對去噪結果的影響減小。

由表2可以發(fā)現，對信號y2而言，方法⑥明顯比方法⑦更好，這是因為信號y2含有的高頻成分較多，把ICA分離得到的原始信號進行直接疊加，減少了高頻成分的損失。也就是對變化較快的高頻信號而言，方法⑥會比方法⑦效果更好。但就MT信號來說，大多數情況下需要獲取的是低頻率信號、而不是高頻信號，因此應用于大地電磁測深數據時，方法⑦整體性能仍然優(yōu)于方法⑥。

圖4對比了方法①與方法⑦的穩(wěn)定性，可以發(fā)現方法⑦的穩(wěn)定性優(yōu)于方法①。其他方法(方法②～方法⑥)與方法⑦的穩(wěn)定性比較結果與之類似，這里不再贅述。

表1 信號y1七種方法去噪處理后的相關系數

表2 信號y2七種方法去噪處理后的相關系數

圖4 方法①與方法⑦的穩(wěn)定性對比圖 (a)數據y1計算結果； (b)數據y2計算結果

4 實際觀測資料去噪處理

選取含有噪聲的兩個測點(分別記為A和B)進行研究，其中測點A采自內蒙古境內，測點B采自河南省境內，因環(huán)境影響，測點B含噪聲較多。對這兩個測點數據運用本文所述的方法進行去噪處理，得到圖5和圖6所示的一系列結果(限于篇幅，僅展示1000個采樣點的結果)。

圖7和圖8分別是A、B兩點去噪前、后的視電阻率曲線和相位曲線?？梢钥闯觯\用本文方法進行去噪處理后，無論視電阻率曲線還是相位曲線，都

圖5 測點A去噪前(左)、后(右)的時間序列對比 (a)Ex； (b)Ey; (c)Hx; (d)Hy; (e)Hz

圖6 測點B去噪前(左)、后(右)的時間序列對比 (a)Ex； (b)Ey; (c)Hx; (d)Hy; (e)Hz

圖7 測點A(a)和測點B(b)去噪前(左)、后(右)視電阻率曲線對比

圖8 測點A(a)和測點B(b)去噪前(左)、后(右)相位曲線對比

變得更加光滑、流暢，即去噪效果都比較好。雖然仍有極少數頻點去噪效果并不理想，但這些頻點在數據預處理時可以作為“飛點”予以刪除，而并不影響整條曲線的處理結果。

5 結論

本文結合小波分析與獨立分量分析的優(yōu)勢，提出基于小波分析的獨立分量分析去噪方法并對大地電磁數據進行去噪。理論上證明了將獨立分量分析應用于大地電磁數據去噪中的可行性；模擬信號仿真實驗表明該方法的去噪穩(wěn)定性優(yōu)于傳統(tǒng)小波閾值去噪方法；實際大地電磁觀測資料去噪結果說明，此方法能比較有效地去除噪聲，特別是在非極低頻區(qū)域，去噪效果尤其明顯。經該方法處理后的視電阻率曲線和相位曲線變得光滑、流暢，為后續(xù)數據處理提供了高質量的數據。

[1] 謝成良.大地電磁測深資料綜合處理軟件系統(tǒng)研究[學位論文].北京:中國地質大學(北京),2013.

[2] 程正璞,胡祥云,李燁等.民豐凹陷大地電磁探測研究.石油地球物理勘探,2016,51(2):391-403. Cheng Zhengpu,Hu Xiangyun,Li Ye et al.A magnetotelluric survey in Minfeng Sag.OGP,2016,51(2):391-403.

[3] 嚴良俊,胡文寶,陳清禮等.遠參考MT方法及其在南方強干擾地區(qū)的應用.石油天然氣學報,1998，18(4):34-38. Yan Liangjun,Hu Wenbao,Chen Qingli et al.Application of remote reference MT to noisy area.Journal of Oil and Gas,1998，18(4):34-38

[4] 蔡劍華,龔玉蓉,王先春.基于Hilbert-Huang變換的大地電磁測深數據處理.石油地球物理勘探,2009,44(5):617-621,629. Cai Jianhua,Gong Yurong,Wang Xianchun.Magnetotelluric data processing based on Hilbert-Huang transform in oil and gas exploration.OGP,2009,44(5):617-621，629.

[5] 張全勝,王家映.大地電磁測深資料的去噪方法.石油地球物理勘探,2004，49(增刊1):17-23. Zhang Quansheng,Wang Jiaying.The method of de-noising for magnetotelluric data.OGP,2004，49(S1):17-23.

[6] 蔡劍華,王先春.EMD在大地電磁信號分析中的問題及解決方法.石油地球物理勘探,2016,51(1):204-210. Cai Jianhua,Wang Xianchun.Problems and solutions of EMD in magnetotelluric signal analysis.OGP,2016,51(1):204-210.

[7] 湯井田,李晉,肖曉等.數學形態(tài)濾波與大地電磁噪聲壓制.地球物理學報,2012,55(5):1784-1793. Tang Jingtian,Li Jin,Xiao Xiao et al.Mathematical morphology filtering and noise suppression of magnetotelluric sounding data.Chinese Journal of Geophy-sics,2012,55(5):1784-1793.

[8] Trad D O,Travassos J M.Wavelet filtering of magnetotelluric data.Geophysics,2000,65(2):482-491.

[9] Manoj C,Nagarajan N.The application of artificial neural networks to magnetotelluric time-series analysis.Geophysical Journal of the Royal Astronomical Society,2010,153(2):409-423.

[10] Weckmann U,Magunia A,Ritter O.Effective noise se-paration for magnetotelluric single site data processing using a frequency domain selection scheme.Geophysical Journal International,2005,161(3):635-652.

[11] 何蘭芳,王緒本,何展翔等.MT時間序列的小波去噪分析.地震地質,2001,23(2):222-226. He Lanfang,Wang Xuben,He Zhanxiang et al.Wavelet-based denoising of MT time series.Seismology & Geology,2001,23(2):222-226.

[12] 嚴家斌,劉貴忠.基于小波變換的脈沖類電磁噪聲處理.煤田地質與勘探,2007,35(5):61-65. Yan Jiabin,Liu Guizhong.Like-impulse electromagnetic noise processing based wavelet transform.Coal Geology & Exploration,2007,35(5):61-65.

[13] 蔡劍華,肖曉.基于小波自適應閾值去噪的MT信號處理方法.地球物理學進展,2015,30(6):2433-2439. Cai Jianhua,Xiao Xiao.Method of processing magnetotelluric signal based on the adaptive threshold wavelet.Progress in Geophysics,2015,30(6):2433-2439.

[14] 王書明,王家映.大地電磁信號統(tǒng)計特征分析.地震學報,2004,26(6):669-674. Wang Shuming,Wang Jiaying.Analysis on statistic characteristics of magnetotelluric signal.Acta Seismologica Sinica,2004,26(6):669-674.

[15] 王書明,王家映.關于大地電磁信號非最小相位性的討論.地球物理學進展,2004,19(2):216-221. Wang Shuming,Wang Jiaying.Discussion on the non-minimum phase of magnetotelluric signals.Progress in Geophysics,2004,19(2):216-221.

[16] Rinen A,Oja E.A fast fixed-point algorithm for independent component analysis.MIT Press,USA，1997,1483-1492.

[17] Hyvarinen A,Oja E.Independent component analysis:algorithms and applications.Neural Networks,2000,13(4):411-430.

[18] 楊福生,洪波.獨立分量分析的原理與應用.北京:清華大學出版社,2006,97-101.

[19] 海韋里恩,卡爾胡恩,奧亞.獨立成分分析.北京:電子工業(yè)出版社,2014,138-152.

[20] 盛驟,謝式千.概率論與數理統(tǒng)計及其應用.北京:高等教育出版社,2004,87.

[21] 余先川,胡丹.盲源分離理論與應用.北京:科學出版社,2011,35-42,51-104.

[22] 陳艷,何英,朱小會.基于小波變換的獨立分量分析及其在圖像分離中的應用.現代電子技術,2007,30(24):131-134. Chen Yan,He Ying,Zhu Xiaohui.ICA based wavelet transform and its application to image separation.Modern Electronics Technique,2007,30(24):131-134.

[23] 胡昌華.基于MATLAB的系統(tǒng)分析與設計—小波分析.陜西西安:西安電子科技大學出版社,1999,223-225.

[24] Donoho D L,Johnstone I M.Ideal spatial adaptation by wavelet shrinkage.Biometrika,1994,81(3):425-455.