黃澤佼，徐子?xùn)|，羅晗，黃遠(yuǎn)生

(1.海南省海洋地質(zhì)資源與環(huán)境重點(diǎn)實(shí)驗(yàn)室，海南 ?？?570206；2.海南水文地質(zhì)工程地質(zhì)勘察院，海南 ?？?571100)

0 引言

大地電磁測(cè)深法是一種以天然交變電磁場(chǎng)為場(chǎng)源的電磁勘探方法，其特點(diǎn)為：信噪比弱、頻帶分布較寬、極易受到各種噪聲的干擾，是典型非線性、非平穩(wěn)信號(hào)[1-2]，加之人文干擾噪聲日益加劇，嚴(yán)重影響了電磁勘探的處理和解釋精度。自大地電磁法提出以來，隨著解釋需求的提高，去噪方法一直都是學(xué)者們所關(guān)心的話題。

傅里葉變換、小波變換等方法常用來分析復(fù)雜的信號(hào)，但其具有較大的局限性[3]：傅里葉變換只適合分析頻率固定的平穩(wěn)信號(hào)，并不適合分析非平穩(wěn)、非線性的大地電磁信號(hào)；而小波變換本質(zhì)上是一種窗口可調(diào)的傅里葉變換，其窗口內(nèi)的信號(hào)必須是平穩(wěn)的，并沒有擺脫傅里葉變換的束縛，分解效果很大程度上依賴于小波基函數(shù)的選取。另外，小波變換是非適應(yīng)性的，小波基一經(jīng)選定，在整個(gè)信號(hào)分析過程中是固定不變的。

1998年，美國(guó)國(guó)家航空航天局(NASA)首席專家Norden E.Huang(黃鍔)院士在Proceeding of the Royal Society of London(英國(guó)皇家學(xué)會(huì)會(huì)刊)發(fā)表了一篇經(jīng)典文章，一個(gè)全新的時(shí)頻分析方法——希爾伯特—黃變換(Hilbert-Huang transformation, HHT)由此產(chǎn)生[4-5]。通過HHT對(duì)信號(hào)進(jìn)行一維經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)，將信號(hào)分解為各階本征模態(tài)函數(shù)(intrinsic mode function，IMF)和一個(gè)剩余分量(residue, res)，各階的IMF頻率由高到低變化，且每一階的IMF分量有其自身的物理意義。再對(duì)IMF進(jìn)行Hilbert變換，可得到包含時(shí)間—能量—頻率三維離散時(shí)頻譜的分布特征。希爾伯特—黃變換不僅具有多分辨率的特性，而且具有自適應(yīng)性，該方法的提出很好地解決了非線性、非平穩(wěn)信號(hào)的時(shí)頻分析問題[6-7]?；Ｈ稹锏萚8-9]引入HHT方法對(duì)EH-4原始時(shí)間序列進(jìn)行時(shí)頻分析，根據(jù)EMD將復(fù)雜信號(hào)分解為隨頻率變換的各階本征模態(tài)函數(shù)，對(duì)受干擾的原始時(shí)間序列進(jìn)行處理，得到了很好的效果。但是其文章沒有對(duì)EMD分解存在的一些問題，如“模態(tài)混疊”、“端點(diǎn)效應(yīng)”等進(jìn)行討論及處理。

本文主要針對(duì)EH-4數(shù)據(jù)中常見的干擾噪聲——工頻噪聲，采用希爾伯特—黃變換(HHT)進(jìn)行去噪處理。對(duì)于工頻噪聲，由于我國(guó)工業(yè)電流的基頻并不是穩(wěn)定的50 Hz，而是在其左右波動(dòng)，應(yīng)用其他的去噪方法效果不佳，尤其是陷波法。本文采用希爾伯特—黃變換(HHT)對(duì)實(shí)測(cè)的時(shí)間序列進(jìn)行EMD分解，結(jié)合各階本征模態(tài)分量(IMF)的希爾伯特譜，在一定程度上能將工頻噪聲的影響去除，但從重構(gòu)后的信號(hào)上看，EMD分解過程中會(huì)產(chǎn)生嚴(yán)重的“模態(tài)混疊”[10]及“端點(diǎn)效應(yīng)”。針對(duì)此問題，本文引入聚合經(jīng)驗(yàn)?zāi)B(tài)(EEMD)法，很好地抑制了“模態(tài)混疊”及“端點(diǎn)效應(yīng)”。通過實(shí)測(cè)數(shù)據(jù)表明，基于聚合經(jīng)驗(yàn)?zāi)B(tài)(EEMD)的希爾伯特—黃變換(HHT)是一種有效的去噪手段，為大地電磁信號(hào)的去噪提供了一條有效的路徑。

1 大地電磁信號(hào)噪聲干擾分類及工頻噪聲的特點(diǎn)

嚴(yán)家斌[11]將噪聲歸類分為：人文噪聲、隨機(jī)噪聲、場(chǎng)源噪聲和地質(zhì)噪聲。有些干擾信號(hào)在時(shí)間序列上具有明顯特征，如工頻噪聲等；有的干擾信號(hào)在頻率域上有很強(qiáng)的特征，如地磁噪聲等；而有的干擾信號(hào)在時(shí)間域上和頻率域上均無任何明顯特征，如地質(zhì)噪聲等。所以，噪聲類型不同，其表現(xiàn)出的電磁特征也是不一樣，因此，在物探工作過程中，研究噪聲的形成機(jī)制、分布特征及規(guī)律以及壓制與消除噪聲的方法是至關(guān)重要的。

本文主要考慮的是人文噪聲中尤為常見的工頻噪聲。人文噪聲源于社會(huì)生產(chǎn)活動(dòng)中產(chǎn)生的電磁噪聲，主要有高壓電線、無線電通訊、鐵路公路等。人文噪聲的信號(hào)強(qiáng)度比天然電磁信號(hào)強(qiáng)幾十甚至幾百倍，常會(huì)造成視電阻率曲線病態(tài)或發(fā)散，嚴(yán)重影響電磁勘探的效果，因此在采集數(shù)據(jù)時(shí)盡可能要偏離供電線，最好是能協(xié)調(diào)當(dāng)?shù)毓╇姴块T停電。但是有時(shí)不可避免的會(huì)在有供電環(huán)境的干擾下進(jìn)行信號(hào)的觀測(cè)采集，當(dāng)外在環(huán)境影響無法剔除時(shí)，這就迫切需要數(shù)據(jù)處理技術(shù)的提高，所以，如何在強(qiáng)工頻干擾的環(huán)境下提取弱的有效信號(hào)就顯得尤為重要了。

2 HHT在工頻噪聲去除中的應(yīng)用

2.1 經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)分解過程

經(jīng)驗(yàn)?zāi)B(tài)分解就是將信號(hào)x(t)分解成一系列滿足條件的IMF的過程，具體實(shí)現(xiàn)步驟[3,12-14](圖1)為：

圖1 一維經(jīng)驗(yàn)?zāi)B(tài)分解的篩分流程Fig.1 Sieve graph divided by empirical mode decomposition(EMD)

1)獲得信號(hào)x(t)全部極值點(diǎn)，并通過插值方法構(gòu)建上下包絡(luò)線的線性方程μ0(t)、d0(t)，一般常用的插值方法為三次樣條插值。

2)通過上下包絡(luò)線曲線方程計(jì)算得到平均值曲線，得到m0(t)=[μ0(t)+d0(t)]/2；記h1(t)=x(t)-m0(t)。

3)判斷h1(t)是否滿足上面的IMF條件，若不滿足，令x(t)=h1(t)，重復(fù)步驟1)、2)，直到hn(t)滿足IMF條件，記c1=hn(t)，即得到第一個(gè)IMF。

4)將原始信號(hào)減去c1(t)，得到分解剩余項(xiàng)，即Res1(t)=x(t)，對(duì)Res1(t)重復(fù)前4個(gè)步驟的分解，即可獲得一系列的單分量信號(hào)ci(t)(i=1,2,…,n)。

5)事先給定特定的值，當(dāng)剩余項(xiàng)小于該值時(shí)，則分解結(jié)束，可得：

2.2 經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)的模態(tài)混疊問題

當(dāng)信號(hào)存在跳躍性變化或是間斷時(shí)，分解過程中一些時(shí)間尺度就會(huì)丟失，信號(hào)本身存在的極值點(diǎn)分布不均勻，這時(shí)就會(huì)造成分解的混亂，出現(xiàn)模態(tài)混疊的情況，這時(shí)候所分解得到的各階IMF分量就不存在所謂的物理意義。模態(tài)混疊包括兩個(gè)方面：①一個(gè)本征模態(tài)函數(shù)IMF中同時(shí)包含了不屬于同一頻率內(nèi)的兩個(gè)頻率；②同一個(gè)本征模態(tài)函數(shù)包含了尺度差異較大的信號(hào)分量[15-17]。

建立如圖2待分解的仿真信號(hào)，其中圖3為模擬信號(hào)的分量。由圖3可知，待分解的仿真信號(hào)由3個(gè)信號(hào)分量組成，其中有1個(gè)信號(hào)分量為間斷跳躍性信號(hào)、1個(gè)正弦信號(hào)及1個(gè)線性信號(hào)。

圖2 待分解的仿真電壓信號(hào)Fig.2 The simulation voltage signal to be decomposed

a—間斷性跳躍信號(hào)；b—正弦信號(hào)；c—線性信號(hào)a—intermittent skip signal; b—sinusoidal signal; c—linear signal圖3 模擬信號(hào)的分量Fig.3 The components of simulation signal

運(yùn)用EMD方法對(duì)待分解的仿真信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，得到如圖4。從原理上分析，由于待分解的信號(hào)中存在間斷性跳躍信號(hào)，會(huì)造成包絡(luò)線的突變，EMD分解時(shí)會(huì)出現(xiàn)模態(tài)混疊現(xiàn)象。事實(shí)上，從圖4中可知，各階的IMF分量，即不同的頻率之間互相干擾嚴(yán)重，EMD并不能夠?qū)⒋纸獾?個(gè)信號(hào)分量完全分解出來。

圖4 EMD分解結(jié)果Fig.4 The result of EMD decomposition

圖5為某工區(qū)的實(shí)測(cè)數(shù)據(jù)，由圖5原始信號(hào)可以明顯地看到，數(shù)據(jù)受到了很強(qiáng)的工頻干擾，完全無法觀察到天然電磁場(chǎng)的特性。將實(shí)際信號(hào)通過EMD自適應(yīng)地分解得到18階固有本征模態(tài)函數(shù)(IMF)，如圖5(b)～(s)所示。通過分析可知，IMF5～I(xiàn)MF10分量包含了不屬于同一頻率段的信息，每一階的IMF分量并沒有表現(xiàn)出某一尺度范圍的模態(tài)，即同一階的IMF分量表現(xiàn)出兩個(gè)甚至多個(gè)頻率的信息，各階IMF分量彼此間出現(xiàn)了嚴(yán)重的模態(tài)混疊現(xiàn)象。

圖5 實(shí)測(cè)EH-4信號(hào)EMD分解及各階本征模態(tài)函數(shù)Fig.5 EMD decomposition and levels of intrinsic mode function based on measured EH-4 signal

2.3 聚合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)抑制模態(tài)混疊

EEMD的目的在于使白噪聲相互抵消，其原理[18]為：高斯白噪聲的頻率分布很均勻，當(dāng)信號(hào)混入高斯白噪聲之后，信號(hào)的極值點(diǎn)就會(huì)發(fā)生改變，使極值點(diǎn)分布更為均勻，包絡(luò)線不存在畸變的現(xiàn)象，從而抑制了模態(tài)混疊現(xiàn)象。另外，高斯白噪聲具有隨機(jī)性，每次加入的噪聲信號(hào)都是獨(dú)立的，所以，只有加入了足夠的噪聲組，利用噪聲的不相關(guān)性使其均值為零，才能消除高斯白噪聲對(duì)真實(shí)信號(hào)的影響。EEMD的分解步驟[15-18]為：

1)在原始信號(hào)x(t)的N次分解時(shí)(N>1)加入均值為零的高斯白噪聲；

2)對(duì)加有高斯白噪聲的信號(hào)進(jìn)行EMD分解，得到K個(gè)IMF和一個(gè)剩余分量Res；

3)由于高斯白噪聲屬于隨機(jī)噪聲，具有不相關(guān)的特性，其均值為零。將各階相對(duì)應(yīng)的IMF分量進(jìn)行求和取總平均，以消除多次加入高斯白噪聲對(duì)真實(shí)IMF和Res的影響，得到最終EEMD分解結(jié)果。

對(duì)上一節(jié)提到的仿真信號(hào)采用EEMD進(jìn)行分解：由圖6可知，通過加入的白噪聲，EEMD分解很好地分離了3個(gè)模擬信號(hào)分量，模態(tài)混疊現(xiàn)象得到了很好地抑制，使得到的各階IMF分量有實(shí)質(zhì)的物理意義。

圖6 EEMD分解結(jié)果Fig.6 The result of EEMD decomposition

上小節(jié)提到，對(duì)實(shí)測(cè)數(shù)據(jù)的原始時(shí)間序列進(jìn)行EMD分解，由于模態(tài)混疊現(xiàn)象，在信號(hào)重構(gòu)的過程中就有可能導(dǎo)致有效信號(hào)也被剔除，各階IMF分量也失去了分解的意義。對(duì)此，本文采用EEMD方法對(duì)原始數(shù)據(jù)進(jìn)行分解處理。本次EEMD分解采用的噪聲方差為0.01，組數(shù)為400的白噪聲，經(jīng)分解得到各階本征模態(tài)函數(shù)，見圖7。由圖7可知，各階的IMF分量只包含其自身頻率相近的信號(hào)成分，EMD分解中的模態(tài)混疊現(xiàn)象已經(jīng)被EEMD完全抑制。

圖7 實(shí)測(cè)信號(hào)EEMD分解及各階本征模態(tài)函數(shù)Fig.7 EEMD decomposition and each intrinsic mode function based on measured EH-4 signal

2.4 HHT去噪過程分析

圖8為信號(hào)通過HHT變換后得到的時(shí)間—能量—頻率三維離散時(shí)頻譜。EH-4低頻段的頻率為10～1 000 Hz，從希爾伯特譜上可以看出，這個(gè)區(qū)域內(nèi)的能量很強(qiáng)，且分布不均勻。在0～200 Hz內(nèi)可以看到，50 Hz左右有一紅色的、能量很強(qiáng)且頻率固定的水平條帶；另外在150 Hz左右也存在類似的高能量水平條帶?？梢耘袛噙@2個(gè)能量高的條帶是由50 Hz及其諧波所引起的工頻噪聲。另外，從希爾伯特譜也可以看到，圖中的能量點(diǎn)越多就表示被分解的信號(hào)的能量就越強(qiáng)，越有利于數(shù)據(jù)的處理與解釋。所以，在實(shí)際數(shù)據(jù)處理中，可以通過時(shí)間序列的希爾伯特譜，將能量弱的疊加去除，保留能量強(qiáng)的疊加，這樣在一定程度上可以增加信號(hào)的信噪比。

圖8 EH-4實(shí)測(cè)信號(hào)時(shí)間—頻率—能量希爾伯特譜Fig.8 Time-frequency-energy spectrum of EH-4 signal

將EEMD分解的各階本征模函數(shù)做希爾伯特變換，進(jìn)行頻譜譜分析，得到圖9。分析找出哪階IMF是由于工頻電引起的，將該階IMF置零，然后將信號(hào)進(jìn)行重構(gòu)，重構(gòu)信號(hào)在一定程度上就能消除工頻電的干擾。從圖9可知，在頻率域中，EEMD分解過程表現(xiàn)為從高頻到低頻的濾波過程。具體表現(xiàn)為：IMF1～I(xiàn)MF4分量能夠很好地體現(xiàn)原始電磁信號(hào)的高頻細(xì)節(jié)信息；IMF5分量在150 Hz左右幅值表現(xiàn)為突然性的“尖窄”跳躍且變化較大，同樣IMF6分量及IMF7分量在50 Hz左右幅值也存在同一規(guī)律的異常跳躍；IMF8～I(xiàn)MF10能夠很好地體現(xiàn)原始電磁信號(hào)的低頻細(xì)節(jié)信息。綜上所述，從各階本征模態(tài)分量的頻譜中可以看出IMF5～I(xiàn)MF7分量中包含了異常信息且與工頻噪聲的特征吻合較好，應(yīng)作為噪聲源去掉。

圖9 EEMD分解的各階IMF分量的頻譜Fig.9 The frequency spectrum of IMF decomposed by EEMD

分別將EMD與EEMD分析處理后的信號(hào)進(jìn)行重構(gòu)，可得圖10。由EMD 分解重構(gòu)后的信號(hào)可知，在重構(gòu)后信號(hào)的首部和尾部出現(xiàn)了信號(hào)的變形(如圖中紅框內(nèi)所示)，這種現(xiàn)象稱之為“端點(diǎn)效應(yīng)”[5]。

a—原始信號(hào);b—EEMD分解重構(gòu)后的信號(hào);c—EMD分解重構(gòu)后的信號(hào);d—信號(hào)重構(gòu)后的誤差a—original signal ;b—signal reconfiguration decomposed by EEMD;c—signal reconfiguration decomposed by EMD; d—the error of signal reconstruction圖10 去噪后信號(hào)重構(gòu)結(jié)果Fig.10 The result of signal reconfiguration after denoising

這種現(xiàn)象引起的原因?yàn)?，?dāng)信號(hào)的邊界端點(diǎn)不是極值點(diǎn)時(shí)，這就導(dǎo)致構(gòu)成上下包絡(luò)線的三次樣條曲線在數(shù)據(jù)序列的兩端出現(xiàn)發(fā)散。解決的方式有兩種:其一，對(duì)于短數(shù)據(jù)序列，可以將原始信號(hào)進(jìn)行擴(kuò)邊，即對(duì)原始數(shù)據(jù)的首尾各加一定量的數(shù)值，在運(yùn)算結(jié)束后進(jìn)行裁邊，恢復(fù)序列的原始長(zhǎng)度。另一種方法為“掐頭去尾”[19-20]。另外，其重構(gòu)后的信號(hào)局部有非正常的起跳(如圖中綠框內(nèi)所示)，預(yù)測(cè)其引起的原因?yàn)橛捎谀B(tài)混疊的原因?qū)⑦^多的IMF分量刪除，使相應(yīng)的有效信號(hào)丟失所造成。

由EEMD 分解重構(gòu)后的信號(hào)可知，EEMD在一定的程度上解決了邊界問題所引起的“端點(diǎn)效應(yīng)”，且重構(gòu)的信號(hào)其間沒有不正常的波動(dòng)，可見模態(tài)混疊現(xiàn)象得到明顯的抑制，信號(hào)恢復(fù)了天然信號(hào)所具備的特點(diǎn)，信號(hào)分解重構(gòu)后的誤差正是所去除的工頻噪聲，信號(hào)整體上平穩(wěn)變化且頻率是固定的。

3 EH-4去工頻噪聲實(shí)例分析

圖11為某勘查區(qū)某測(cè)點(diǎn)的EH-4時(shí)間序列，該點(diǎn)旁有輸電線，采集受到嚴(yán)重的干擾。由圖11可知，測(cè)點(diǎn)高頻段的電道未見明顯的工頻干擾，高頻段的磁道、中頻段的電道磁道及低頻段的電道磁道均受到工頻干擾較為明顯。本文采用的噪聲組數(shù)為200組，均方差0.001，根據(jù)上述EEMD分析方法，去除工頻噪聲并將信號(hào)進(jìn)行重構(gòu)，將重構(gòu)后的時(shí)間序列寫成二進(jìn)制Y文件，導(dǎo)入EH-4自帶的處理系統(tǒng)IMAGEM。去噪后的時(shí)間序列如圖12所示，各頻段的噪聲得到了很好地壓制。

(a) 高頻段 (b)中頻段 (c)低頻段圖11 某測(cè)點(diǎn)的原始時(shí)間序列示意Fig.11 The original time series of measuring point

(a)高頻段 (b)中頻段 (c)低頻段圖12 去工頻噪聲后的時(shí)間序列示意Fig.12 The time series after de-noising power frequency

通過IMAGEM軟件處理所得到的阻抗視電阻率曲線如圖13所示，圖中紅色框線顯示：無論是TE還是TM模式，中、低頻信號(hào)在未去除工頻噪聲之前，有效信號(hào)被嚴(yán)重壓制，視電阻率曲線出現(xiàn)了病態(tài)或發(fā)散的現(xiàn)象，曲線在中、低頻段存在大量的缺點(diǎn)、間斷，走勢(shì)形態(tài)完全無法判斷。而在中、低頻信號(hào)去除工頻噪聲之后，視電阻率曲線的中低頻段曲線缺點(diǎn)及間斷現(xiàn)象得到較好地改善，曲線的走勢(shì)形態(tài)趨于明顯，尤其是TM模式，中、低頻的電阻率曲線更加的平滑?？梢?，去噪的效果很明顯。

(a)TE模式 (b)TM模式圖13 TE模式(a)和TM模式(b)去工頻噪聲前后視電阻率曲線對(duì)比Fig.13 Contrast of apparent resistivity curve of TE mode(a) and TM mode(b) before and after de-noising power frequency

4 結(jié)論

本文研究了EH-4數(shù)據(jù)處理的方法，以工頻噪聲作為主要研究對(duì)象，通過分析得出以下結(jié)論:

1) 通過對(duì)野外實(shí)測(cè)大地電磁信號(hào)進(jìn)行希爾伯特—黃(HHT)變換、二維Hilbert時(shí)頻分析，結(jié)果顯示HHT變換能充分體現(xiàn)大地電磁信號(hào)的特征，該方法在一定程度能去除工頻干擾噪聲。本文應(yīng)用聚合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)法，通過實(shí)際數(shù)據(jù)的處理，較好地抑制了EMD分解所引起的模態(tài)混疊現(xiàn)象，在去除信號(hào)的同時(shí)保留了更多有效信號(hào)，在一定程度上也很好地解決了“端點(diǎn)效應(yīng)”的問題。

2) 采用基于EEMD的HHT變換對(duì)實(shí)際數(shù)據(jù)進(jìn)行去噪，在中、低頻信號(hào)去除工頻噪聲之后，視電阻率曲線的中、低頻段曲線間斷的現(xiàn)象得到明顯的改善，曲線的走勢(shì)形態(tài)趨于明顯，充分說明基于EEMD的HHT變換在去工頻噪聲中的可行性。