大地電磁的小波變換—獨立分量分析去噪

曹小玲 劉開元 嚴(yán)良俊

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

1 引言

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

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

2 獨立分量分析

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

圖1 ICA的算法流程圖

2.1 數(shù)據(jù)預(yù)處理

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

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

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

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

2.2 FastICA算法

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

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

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

(1)

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

X=AS

(2)

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

設(shè)離散隨機變量y的熵定義為

(3)

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

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

(4)

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

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

(5)

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

FastICA算法原理如下。

首先，通過對E[g(WTX)]進行優(yōu)化獲得y的負(fù)熵的最大近似值。根據(jù)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(shù)′表示g的一次導(dǎo)數(shù)。由于數(shù)據(jù)被白化，有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 基于小波分析的獨立分量分析

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

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

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

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層分解)； ⑤自適應(yīng)默認(rèn)閾值去噪法(5層分解)； ⑥重構(gòu)信號Y=R+∑T′去噪法(T′為ICA分離后未經(jīng)處理的原始結(jié)果信號)(5層分解)； ⑦由式(12)所示的去噪法(5層分解)。得到的實驗結(jié)果見圖3。限于篇幅，這里只顯示加入噪聲(SNR=10)時y1的實驗結(jié)果。

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

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

(13)

顯然Y、O包含的元素個數(shù)相等。

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

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

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

表1 信號y1七種方法去噪處理后的相關(guān)系數(shù)

表2 信號y2七種方法去噪處理后的相關(guān)系數(shù)

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

4 實際觀測資料去噪處理

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

圖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)去噪前(左)、后(右)相位曲線對比

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

5 結(jié)論

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

[1] 謝成良.大地電磁測深資料綜合處理軟件系統(tǒng)研究[學(xué)位論文].北京:中國地質(zhì)大學(xué)(北京),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] 嚴(yán)良俊,胡文寶,陳清禮等.遠(yuǎn)參考MT方法及其在南方強干擾地區(qū)的應(yīng)用.石油天然氣學(xué)報,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變換的大地電磁測深數(shù)據(jù)處理.石油地球物理勘探,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] 湯井田,李晉,肖曉等.數(shù)學(xué)形態(tài)濾波與大地電磁噪聲壓制.地球物理學(xué)報,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時間序列的小波去噪分析.地震地質(zhì),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] 嚴(yán)家斌,劉貴忠.基于小波變換的脈沖類電磁噪聲處理.煤田地質(zhì)與勘探,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] 蔡劍華,肖曉.基于小波自適應(yīng)閾值去噪的MT信號處理方法.地球物理學(xué)進展,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)計特征分析.地震學(xué)報,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] 王書明,王家映.關(guān)于大地電磁信號非最小相位性的討論.地球物理學(xué)進展,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] 楊福生,洪波.獨立分量分析的原理與應(yīng)用.北京:清華大學(xué)出版社,2006,97-101.

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

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

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

[22] 陳艷,何英,朱小會.基于小波變換的獨立分量分析及其在圖像分離中的應(yīng)用.現(xiàn)代電子技術(shù),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)分析與設(shè)計—小波分析.陜西西安:西安電子科技大學(xué)出版社,1999,223-225.

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