希爾伯特—黃變換(HHT)在EH-4數據去噪處理中的應用

黃澤佼，徐子東，羅晗，黃遠生

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

0 引言

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

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

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

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

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

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

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

2 HHT在工頻噪聲去除中的應用

2.1 經驗模態(tài)分解(EMD)分解過程

經驗模態(tài)分解就是將信號x(t)分解成一系列滿足條件的IMF的過程，具體實現步驟[3,12-14](圖1)為：

圖1 一維經驗模態(tài)分解的篩分流程Fig.1 Sieve graph divided by empirical mode decomposition(EMD)

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

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

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

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

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

2.2 經驗模態(tài)分解(EMD)的模態(tài)混疊問題

當信號存在跳躍性變化或是間斷時，分解過程中一些時間尺度就會丟失，信號本身存在的極值點分布不均勻，這時就會造成分解的混亂，出現模態(tài)混疊的情況，這時候所分解得到的各階IMF分量就不存在所謂的物理意義。模態(tài)混疊包括兩個方面：①一個本征模態(tài)函數IMF中同時包含了不屬于同一頻率內的兩個頻率；②同一個本征模態(tài)函數包含了尺度差異較大的信號分量[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方法對待分解的仿真信號進行經驗模態(tài)分解，得到如圖4。從原理上分析，由于待分解的信號中存在間斷性跳躍信號，會造成包絡線的突變，EMD分解時會出現模態(tài)混疊現象。事實上，從圖4中可知，各階的IMF分量，即不同的頻率之間互相干擾嚴重，EMD并不能夠將待分解的3個信號分量完全分解出來。

圖4 EMD分解結果Fig.4 The result of EMD decomposition

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

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

2.3 聚合經驗模態(tài)分解(EEMD)抑制模態(tài)混疊

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

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

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

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

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

圖6 EEMD分解結果Fig.6 The result of EEMD decomposition

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

圖7 實測信號EEMD分解及各階本征模態(tài)函數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ū)域內的能量很強，且分布不均勻。在0～200 Hz內可以看到，50 Hz左右有一紅色的、能量很強且頻率固定的水平條帶；另外在150 Hz左右也存在類似的高能量水平條帶。可以判斷這2個能量高的條帶是由50 Hz及其諧波所引起的工頻噪聲。另外，從希爾伯特譜也可以看到，圖中的能量點越多就表示被分解的信號的能量就越強，越有利于數據的處理與解釋。所以，在實際數據處理中，可以通過時間序列的希爾伯特譜，將能量弱的疊加去除，保留能量強的疊加，這樣在一定程度上可以增加信號的信噪比。

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

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

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

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

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

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

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

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

圖11為某勘查區(qū)某測點的EH-4時間序列，該點旁有輸電線，采集受到嚴重的干擾。由圖11可知，測點高頻段的電道未見明顯的工頻干擾，高頻段的磁道、中頻段的電道磁道及低頻段的電道磁道均受到工頻干擾較為明顯。本文采用的噪聲組數為200組，均方差0.001，根據上述EEMD分析方法，去除工頻噪聲并將信號進行重構，將重構后的時間序列寫成二進制Y文件，導入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模式，中、低頻信號在未去除工頻噪聲之前，有效信號被嚴重壓制，視電阻率曲線出現了病態(tài)或發(fā)散的現象，曲線在中、低頻段存在大量的缺點、間斷，走勢形態(tài)完全無法判斷。而在中、低頻信號去除工頻噪聲之后，視電阻率曲線的中低頻段曲線缺點及間斷現象得到較好地改善，曲線的走勢形態(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 結論

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

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

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