基于Hilbert-Huang 變換的大地電磁去噪研究

陳鈞,嚴(yán)良俊,周磊

(1.油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(長江大學(xué))，湖北 武漢 430100；2.非常規(guī)油氣省部共建協(xié)同創(chuàng)新中心，湖北 武漢 430100)

0 引言

大地電磁在電磁法勘探領(lǐng)域有著廣泛的應(yīng)用，但是其信號弱、易受干擾問題一直難以解決。隨著研究的深入發(fā)展，國內(nèi)外大地電磁測深前處理已經(jīng)取得許多成果，研究者提出的去噪方法種類繁多：胡家華等[1]對大地電磁噪聲源進(jìn)行了分析并認(rèn)為人文噪聲是主要干擾源；王書明等[2-3]比較了大地電磁功率譜和高功率譜，同時對大地電磁信號的統(tǒng)計特征進(jìn)行了分析。研究者主要利用大地電磁信號和噪聲的不同特征對他們加以分離或者根據(jù)干擾特征壓制噪聲：王輝等[4]使用改進(jìn)的遠(yuǎn)參考方法對礦區(qū)近場源噪聲進(jìn)行壓制并且分析了近場源噪聲的特點(diǎn)；劉俊峰等[5]引入了一種多結(jié)構(gòu)元素組成的濾波器用于消除大地電磁信號噪聲；李晉等[6]運(yùn)用遞歸分析的方法進(jìn)行了大地電磁信號的信號和噪聲識別，然后對噪聲進(jìn)行壓制。本文采用的HHT方法是根據(jù)噪聲和信號的時頻譜特征將它們加以分離。

隨著人文活動的加劇，許多以往的數(shù)據(jù)處理方式已不能有效壓制近場源效應(yīng)，所獲數(shù)據(jù)與后續(xù)工作的精度要求不符。如何有效壓制近場效應(yīng)已經(jīng)成了大地電磁測深技術(shù)前處理的一個重要研究方向。在近場效應(yīng)下，視電阻率曲線表現(xiàn)為一條斜率 45°的線段，相位曲線則為 0°(對數(shù)坐標(biāo))[4]。近年來，Hilbert-Huang 變換在處理數(shù)據(jù)的實(shí)踐中得到廣泛應(yīng)用：Huang等[7]將HHT應(yīng)用于地球物理并且闡述了HHT相比其他處理方法在處理非線性信號上的優(yōu)越性；白大為等[8]將HHT方法應(yīng)用于帶天然氣水合物勘探中提取反射信號；楊洋等[9]將HHT解析包絡(luò)的思想結(jié)合小波變換提出一種新的CSEM信號噪聲評價方法；蔡劍華、湯井田等[10-11]也先后將HHT應(yīng)用于大地電磁噪聲壓制。在驗(yàn)證去噪效果方面，U.Weckmann等[12]提出根據(jù)電磁場強(qiáng)極化和響應(yīng)函數(shù)的誤差分布來檢驗(yàn)去噪效果的方法，他提出正常含噪聲較少信號的分布函數(shù)應(yīng)該是趨近于高斯分布的，并且使用大量的實(shí)際數(shù)據(jù)驗(yàn)證了該觀點(diǎn)。

本文把HHT與加拿大鳳凰公司的SSMT2000系統(tǒng)相結(jié)合，應(yīng)用到某處實(shí)地采集的大地電磁測深數(shù)據(jù)中，并研究了HHT對近場干擾的壓制效果。在介紹HHT原理的基礎(chǔ)上，先數(shù)值模擬了經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)，確認(rèn)了其對復(fù)雜信號的處理能力之后，進(jìn)一步對采集到的大地電磁數(shù)據(jù)電磁各道進(jìn)行分解，根據(jù)希爾伯特時頻譜對信號篩選去噪，最后再結(jié)合SSMT2000的計算結(jié)果和去噪前后的電磁場極化方向的響應(yīng)函數(shù)分布情況進(jìn)行驗(yàn)證。

1 方法原理

Huang等人提出的HHT變換擴(kuò)大了Hilbert Transform的適用范圍，它分為兩個部分：經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)和Hilbert spectrum分析。Huang利用信號的頻率特征，根據(jù)波峰波谷的包絡(luò)線在時間尺度上對信號進(jìn)行分解，指出1個IMF必須滿足以下2個條件：①在整個IMF數(shù)據(jù)序列中，極值點(diǎn)和零點(diǎn)的數(shù)量差應(yīng)當(dāng)不超過1個；②在任何時間點(diǎn)上，由信號序列局部最大值定義的包絡(luò)和局部極小值定義的包絡(luò)，兩者的均值必須為零。滿足以上2個條件的IMF 已被實(shí)踐證明適用于 Hilbert 變換[7]。

用程序進(jìn)行EMD分解流程大致如下。

1)初始信號X(t)0先用插值(一般用三次樣條插值)后數(shù)據(jù)找到波峰值、波谷值，保存波峰、波谷橫坐標(biāo)和波峰、波谷值。

2)對波峰、波谷進(jìn)行插值(一般用三次樣條插值)，因此產(chǎn)生上下包絡(luò)線并保存。

3)用循環(huán)語句將插值后的每個點(diǎn)的波峰、波谷信號求和再除以2，即得到均值信號(記為mik)，并保存在一個新數(shù)組中。

4)計算原始信號和mik的差值，記為hik：hik=X(t)-mik。如果hik符合基本面模式分量的條件，就作為第i個IMF輸出保存：sj=hik；如果不滿足條件，就對hik繼續(xù)前三步。

5)令此時的原始信號減去hik，得到剩余信號ri：ri=X(t)i-sj。將ri作為初始信號繼續(xù)進(jìn)行上面步驟，得到不同的IMF和剩余信號ri，直到剩余信號ri不再有波峰或波谷時，分解結(jié)束。

圖1給出了程序流程。流程中第四步有一個篩選判斷，一般情況下設(shè)定一個停止準(zhǔn)則(本文使用毛煒等提出的準(zhǔn)則3)[13]：

圖1 EMD流程

(1)

最后可以得到

(2)

式中：sj(t)為各個IMF，rn(t)為殘余分解量。

篩選后對剩余的各個IMF分量進(jìn)行Hilbert變換：

(3)

式(3)中：ai(t)、φi(t)是IMF分量的時變幅度和相位，rn(t)一般不進(jìn)行Hilbert變換。在傳統(tǒng)的傅里葉分析中，為了定義局部頻率值起碼需要一個完整的周期。對于非線性非平穩(wěn)信號傅氏變換不夠準(zhǔn)確，于是Huang等人根據(jù)Hilbert-Huang變換提出了新的瞬時頻率定義方法：

(4)

式中N為IMF個數(shù)。式(4)可以表示每一個IMF的瞬時頻率。根據(jù)前文給出的瞬時頻率定義，可以把原始信號(式(3))表示為

(5)

式(5)中實(shí)部是Hilbert譜，一般記為H(ω,t)。

Hilbert 邊際譜：

H(ω,t)=Re[x(t)]，

(6)

(7)

式中：T是總信號跨度；h(ω)為Hilbert邊際譜，即H(ω,t)做時間積分可得到Hilbert邊際譜。

2 經(jīng)驗(yàn)?zāi)B(tài)分解數(shù)值模擬：

以模擬信號y=x+1.5sin(5x)+sin(20x)為例進(jìn)行EMD分解，該信號的頻率由15 s內(nèi)的0.8 Hz和3.2 Hz這2個正弦電場信號和1個線性電場信號組成，對它進(jìn)行分解。圖2a是原始電場信號，圖2b—d中的藍(lán)色信號是分解得到的固有模態(tài)函數(shù)IMF1-2和剩余信號，橙色線為構(gòu)成原始電場信號的正確信號分量。顯然，圖中IMF1對應(yīng)3.2 Hz正弦信號，IMF2對應(yīng) 0.8 Hz正弦信號，而剩余分量則是對應(yīng)線性信號。由于端點(diǎn)效應(yīng)數(shù)據(jù)在端點(diǎn)處少量數(shù)據(jù)點(diǎn)不擬合，從模擬中也能發(fā)現(xiàn)經(jīng)驗(yàn)?zāi)B(tài)分解對信號的分解和還原都很成功。

圖2 數(shù)值模擬模擬經(jīng)驗(yàn)?zāi)B(tài)分解

在對模擬信號分解重構(gòu)成功的基礎(chǔ)上再對該信號添加一個鋸齒波干擾模擬噪聲(圖3)，繼續(xù)驗(yàn)證EMD分解。圖4顯示分離出的IMF5和添加的鋸齒波干擾十分接近，而剩余信號res是分解剩余的線性信號(見圖3)，可見EMD分解成功，壓制了鋸齒波噪聲。

a—加入噪聲后的電場信號;b～f—IMF1～I(xiàn)MF5;g—剩余信號

圖4 實(shí)際噪聲(noise)信號(a)和分離的噪聲(IMF5)信號(b)

3 實(shí)際資料處理

為了進(jìn)一步研究HHT對大地電磁信號近場干擾的壓制效果，以加拿大鳳凰公司V8大地電磁測深系統(tǒng)實(shí)測采集的大地電磁測深數(shù)據(jù)為例進(jìn)行驗(yàn)證。該儀器采集數(shù)據(jù)時記錄電道和磁道里面Ex、Ey、Hx、Hy、Hz5個分量，并將數(shù)據(jù)以二進(jìn)制形式儲存。使用長江大學(xué)嚴(yán)良俊老師課題組開發(fā)的ts view軟件將采集到的tsn文件轉(zhuǎn)化為csv文件(從二進(jìn)制編碼裝換到ASCII碼)，用程序?qū)ζ溥M(jìn)行去噪處理后再將文件替換回tsn文件。

已知對于任意1組正交的電磁場水平分量Ex、Hy和Ey、Hx，有波阻抗關(guān)系：

(10)

根據(jù)阻抗將卡尼亞視電阻率求出。本文使用了甘肅某地采集的大地電磁測深數(shù)據(jù)，同時使用加拿大鳳凰公司的V8測深系統(tǒng)進(jìn)行視電阻率和功率譜的自動篩選，由于其計算過程在SSMT2000軟件里面進(jìn)行，故在此直接展示計算結(jié)果。

采集地區(qū)地勢較平坦空曠。觀察視電阻率曲線(圖5)，發(fā)現(xiàn)TM信號在1～0.1 Hz時曲線呈現(xiàn)出45°形態(tài)，同時這一頻段的相位曲線趨近于0°(相位散點(diǎn)圖中負(fù)的相位加了180°使TE、TM模式下的相位都處于0～180°內(nèi))，TM曲線的視電阻率和相位表現(xiàn)出明顯的近場干擾效應(yīng)；觀測電磁場極化方向響應(yīng)函數(shù)的分布(圖6)，發(fā)現(xiàn)分布函數(shù)的高斯性被嚴(yán)重破壞，電場磁場散點(diǎn)高度聚集，說明信號受到很大的噪聲干擾污染。這些現(xiàn)象表明該數(shù)據(jù)的TM曲線受到嚴(yán)重的近場干擾。

圖5 未去噪時間序列對應(yīng)的死頻帶視電阻率和相位曲線

圖6 未去噪時的極化方向響應(yīng)函數(shù)

將V8測深系統(tǒng)的3個不同采用頻率ts3-ts5(2 400、150、15 Hz)各5道信號(Ex、Hy、Ey、Hx、Hz)分別導(dǎo)出，針對連續(xù)采樣的ts5以及含較多噪聲頻段的ts4加以處理。得到各階IMF之后，利用式(4)對各個IMF進(jìn)行希爾伯特變換求其瞬時頻率，然后根據(jù)Hilbert時頻譜篩選IMF重構(gòu)信號(圖7)。希爾伯特時頻譜刻畫的是瞬時頻率，噪聲頻率含量較少的IMF散點(diǎn)會表現(xiàn)出比較連續(xù)而不是離散的形態(tài)。

圖7 噪聲頻率含量較多IMF的希爾伯特時頻譜

完成噪聲壓制工作后將重構(gòu)的信號寫回二進(jìn)制。使用SSMT2000計算視電阻率和相位(圖8、圖9)并與圖5、圖6進(jìn)行了對比。圖8顯示近場干擾得到了有效壓制，其中相位散點(diǎn)圖中負(fù)的相位加了180°，使得TE、TM模式下的相位都處于0°～180°內(nèi)，經(jīng)過HHT方法壓制后視電阻率曲線和相位曲線不再表現(xiàn)出近場干擾特征；圖9顯示散點(diǎn)趨勢趨近于高斯分布，說明經(jīng)過HHT處理后電場、磁場的極化響應(yīng)函數(shù)高斯性已經(jīng)得到較大修復(fù)。以上結(jié)果表明Hilbert-Huang變換可以對強(qiáng)近場干擾有效壓制。

圖8 去噪之后重新計算的視電阻率和相位

圖9 去噪后的極化方向響應(yīng)函數(shù)

4 結(jié)論

文中數(shù)值模擬了電磁信號的分解重構(gòu)，針對添加了鋸齒波干擾的模擬信號分解重構(gòu)，成功壓制了大部分鋸齒波干擾，同時保證了重構(gòu)波形的準(zhǔn)確性，并總結(jié)探究了經(jīng)驗(yàn)?zāi)B(tài)分解的適用性。

研究發(fā)現(xiàn)：對于非線性復(fù)雜混合信號，經(jīng)驗(yàn)?zāi)B(tài)分解可以將信號分解成簡單的固有模態(tài)信號，根據(jù)希爾伯特邊際譜和時頻譜很容易區(qū)分噪聲(IMF1-2)和其余部分的有用信號(其余信號和剩余信號等)，而且相比于傳統(tǒng)的傅氏變換，經(jīng)驗(yàn)?zāi)B(tài)分解對非穩(wěn)態(tài)非線性信號有著更好的表現(xiàn)。此外，傳統(tǒng)遠(yuǎn)參考方法難以壓制1 Hz和0.1 Hz附近的近場干擾噪聲，而使用V8測深系統(tǒng)自帶的Robust系統(tǒng)結(jié)合Hilbert-Huang變換可以有效壓制近場噪聲。

在處理信號過程中，發(fā)現(xiàn)處理后的信號有部分低頻數(shù)據(jù)丟失，視電阻率和相位的功率譜不夠平滑連續(xù)，人工篩選信號重構(gòu)也可能造成有效信號丟失并且工作效率不高。如何有針對性地精確去噪、減少有效信號的丟失，是下一步研究的方向。