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

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

0 引言

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

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

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

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

1 大地電磁信號噪聲干擾分類及工頻噪聲的特點

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

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

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

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

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

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

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

2)通過上下包絡(luò)線曲線方程計算得到平均值曲線，得到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)，即得到第一個IMF。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2)對加有高斯白噪聲的信號進行EMD分解，得到K個IMF和一個剩余分量Res；

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

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

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

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

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

2.4 HHT去噪過程分析

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

圖8 EH-4實測信號時間—頻率—能量希爾伯特譜Fig.8 Time-frequency-energy spectrum of EH-4 signal

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

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

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

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

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

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

3 EH-4去工頻噪聲實例分析

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

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

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

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

(a)TE模式 (b)TM模式圖13 TE模式(a)和TM模式(b)去工頻噪聲前后視電阻率曲線對比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ù)處理的方法，以工頻噪聲作為主要研究對象，通過分析得出以下結(jié)論:

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

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