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

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

0 引言

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

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

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

1 方法原理

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

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

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

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

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

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

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

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

圖1 EMD流程

(1)

最后可以得到

(2)

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

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

(3)

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

(4)

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

(5)

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

Hilbert 邊際譜：

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

(6)

(7)

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

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

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

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

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

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

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

3 實(shí)際資料處理

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

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

(10)

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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