蔡劍華,胡惟文, 覃業(yè)貴, 陳卿冶

?

MT信號(hào)處理中小波分析與Hilbert-Huang變換的比較

蔡劍華*1,胡惟文1, 覃業(yè)貴2, 陳卿冶1

(1. 湖南文理學(xué)院 物理與電子科學(xué)學(xué)院,湖南 常德, 415000; 2. 解放軍93062部隊(duì), 吉林 吉林, 132102)

以實(shí)測(cè)大地電磁信號(hào)為分析對(duì)象, 對(duì)比研究了小波分析和Hilbert-Huang變換在大地電磁(magneto- telluric, MT)信號(hào)分解和時(shí)頻分析上面的特點(diǎn). 給出了小波分析和Hilbert-Huang變換的數(shù)學(xué)原理, 分析了2種方法在MT信號(hào)分解中的分解方法及其物理意義, 討論了2種方法在MT信號(hào)時(shí)頻譜表征能力上的異同, 為大地電磁信號(hào)分析方法的選擇提供用了有價(jià)值的參考.

大地電磁信號(hào); Hilbert-Huang變換; 小波分析; 分解; 時(shí)頻分析

非線性信號(hào)處理方法是分析、處理非平穩(wěn)大地電磁信號(hào)的有力工具[1], 其中小波分析和Hilbert- Huang變換在地球物理領(lǐng)域應(yīng)用的較為廣泛, 他們都能將信號(hào)進(jìn)行自適應(yīng)的分解, 并得到時(shí)頻譜圖, 揭示信號(hào)的時(shí)變特征, 尤其在希望獲得頻率成分隨時(shí)間變化規(guī)律的地震信號(hào)處理和大地電磁法中應(yīng)用較多[2—4]. 本文以實(shí)測(cè)大地電磁信號(hào)為分析對(duì)象, 對(duì)比研究小波分析和Hilbert-Huang變換在大地電磁信號(hào)分解和時(shí)頻分析上面的優(yōu)缺點(diǎn), 試圖為大地電磁信號(hào)分析方法的選擇提供用有價(jià)值的參考.

1 小波變換

則()被稱為一個(gè)基本小波函數(shù), 將小波函數(shù)()進(jìn)行伸縮和平移變換就可以得到一組小波[5].

在連續(xù)情況下, 小波為:

通過適當(dāng)?shù)剡x擇尺度因子和平移因子, 可得到一個(gè)伸縮窗[6].

在離散情況下, 小波為:

與函數(shù)()相應(yīng)的離散小波變換系數(shù)及其重構(gòu)公式為:

式(7)中是一個(gè)與信號(hào)無關(guān)的常數(shù).

2 Hilbert-Huang變換

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

經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)是依據(jù)信號(hào)本身的時(shí)間尺度特征, 將信號(hào)分解為一組模態(tài)函數(shù)(intrinsic mode function, IMF)[11], IMF滿足以下2個(gè)定義條件[7]: ①對(duì)于一列數(shù)據(jù)極值點(diǎn)和過零點(diǎn)數(shù)目必須相等或至多相差一點(diǎn); ②在任意點(diǎn), 由局部極大點(diǎn)和極小點(diǎn)構(gòu)成的兩條包絡(luò)線的平均值為0. EMD分解方法如下[7]: ①找出信號(hào)()的所有極大點(diǎn)和極小點(diǎn), 將其用三次樣條函數(shù)分別擬合為原數(shù)據(jù)序列的上、下包絡(luò)線, 上、下包絡(luò)線的均值為平均包絡(luò)線1, 將原數(shù)據(jù)序列減去1可得到一個(gè)去掉低頻的新數(shù)據(jù)序列1. 一般1不是一個(gè)平穩(wěn)數(shù)據(jù)序列, 為此重復(fù)以上過程次, 使所得到的平均包絡(luò)線趨于0, 此時(shí)的1n就是第1個(gè)IMF(1), 它表示信號(hào)數(shù)據(jù)中的最高頻成分; ②用()減去1得到一個(gè)去掉高頻成分的新數(shù)據(jù)序列, 重復(fù)步驟(1), 得到一系列c和最后一個(gè)序列r,r代表()的均值或趨勢(shì)項(xiàng). 那么原序列()可表示為IMF分量和一個(gè)殘余項(xiàng)的和:

2.2 Hilbert變換與時(shí)頻譜

獲得IMF分量后, 就可以對(duì)每階IMF(設(shè)為())作Hilbert變換, 得到:

為柯西主值.()和()共同組合為一解析信號(hào)(), 采用極坐標(biāo)形式表示:

其中:

根據(jù)瞬時(shí)頻率的定義:

對(duì)每一階IMF作Hilbert變換, 并求出相應(yīng)解析函數(shù)的幅值譜和瞬時(shí)頻率, 從而原始信號(hào)可以表示為[8]:

可見, 信號(hào)的幅值能表示為時(shí)間、瞬時(shí)頻率的函數(shù)(,), 從而獲得信號(hào)幅值的時(shí)間和頻率分布即Hilbert譜.

3 MT信號(hào)的小波分解與時(shí)頻譜

圖1中EH-4儀器采用低頻工作模式, 采樣率為12 kHz. 由于數(shù)據(jù)量龐大, 此處僅取E分量的4 096個(gè)點(diǎn)進(jìn)行分析. 小波變換在不同尺度下分解信號(hào)是多分辨率的, 小波變換的尺度因子使得小波變換具有“變焦距特性”, 被譽(yù)為數(shù)學(xué)上的顯微鏡, 小波變換中可改變尺度因子得到一個(gè)伸縮窗. 圖2(a)所示為利用sym5小波基對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行分解的小波分解圖. 圖中, a6及d1—d6是小波分解后得到的小波分量, 其中a6為第6層分解的近似系數(shù)、d1—d6為前6層分解的高頻系數(shù), 各層系數(shù)對(duì)應(yīng)于從高頻到低頻的濾波過程.

圖1 原始的大地電磁信號(hào)

小波分解中小波基函數(shù)的選擇不是任意的, 也不是唯一的. 因此, 小波基的選擇是小波分析在實(shí)際應(yīng)用中的一個(gè)重要問題. 如圖2(b)所示為相同的圖1所示信號(hào)采用db5小波基分解后的分解圖. 可見, 其分解系數(shù)與圖2(a)所示是不盡相同的, 尤其是高頻部分的信息. 同一個(gè)問題小波基不同, 其結(jié)果也不同, 甚至相差很大.

圖3所示為利用sym5小波基分析得到的大地電磁信號(hào)時(shí)頻譜圖. 從圖3中可清楚地看出大地電磁信號(hào)能量隨時(shí)頻的分布情況, 譜圖表現(xiàn)出良好的局部化特征, 較精確地描述了大地電磁信號(hào)的時(shí)頻分布. 但也可以清楚地看到在大地電磁信號(hào)小波譜中, 能量呈分散特征, 分別率也不高, 存在混疊的現(xiàn)象, 且除主成分外還存在的諧波成分, 其對(duì)應(yīng)的頻率范圍也被展寬.

圖3 大地電磁信號(hào)的小波時(shí)頻圖

4 MT信號(hào)的EMD分解和時(shí)頻譜

EMD分解技術(shù)在時(shí)間域中表現(xiàn)為從小尺度到大尺度的層層濾波, 在頻率域中EMD分解過程表現(xiàn)為從高頻到低頻的層層分解, 無需選擇基函數(shù), 是完全自適應(yīng)于被分解的信號(hào)的. 采用EMD方法對(duì)圖1所示的大地電磁信號(hào)進(jìn)行分解, 圖4所示為大地電磁信號(hào)的EMD分解圖及其各階對(duì)應(yīng)的功率譜. 從圖4中可知, 信號(hào)被自適應(yīng)的分解為9層, 其中imf1為白噪聲, imf2—imf8為優(yōu)勢(shì)子頻帶, residual為趨勢(shì)項(xiàng), 每一階IMF分量都有不同的振幅和頻率, 這為大地電磁信號(hào)的去噪等后續(xù)工作帶來了方便. 整個(gè)頻段內(nèi)低頻成分的振幅較大, 這與進(jìn)行數(shù)據(jù)采集時(shí)使用儀器的低頻模式采樣的事實(shí)相吻合, 也說明了分解的正確性.

按照式(12)計(jì)算得到信號(hào)的Hilbert能量譜如圖5所示. 從圖5中可以得出: 對(duì)比圖3中小波分析方法得到的時(shí)頻譜圖, 圖3中時(shí)頻譜表現(xiàn)的優(yōu)勢(shì)頻率與Hilbert能量譜中表征的優(yōu)勢(shì)頻率是一致的, 說明信號(hào)的能量主要集中在低頻部分, 這與信號(hào)的實(shí)際吻合得很好. 但圖5所示所示的Hilbert能量譜清晰而詳細(xì)地顯示了大地電磁信號(hào)能量隨時(shí)頻的具體分布. 與圖3相比可以得到: Hilbert譜具有更好的局部化能力和更好的時(shí)頻聚集性, 大地電磁數(shù)據(jù)的突變點(diǎn)、持續(xù)時(shí)段和頻帶能量分布均能夠清晰地顯現(xiàn); HHT方法對(duì)大地電磁數(shù)據(jù)的動(dòng)態(tài)變化過程刻畫得比較清楚, 反應(yīng)了每時(shí)段都有各自的頻率特性、能量差異, 小波時(shí)頻譜圖難以揭示出這些細(xì)微性變化.

圖5 大地電磁信號(hào)的Hilbert 時(shí)頻譜

5 結(jié)論

小波變換是一種被廣泛應(yīng)用于大地電磁信號(hào)分析的方法, 它克服了短時(shí)傅里葉分析中窗口大小不隨頻率變化的缺點(diǎn), 是一種比較理想的信號(hào)處理數(shù)學(xué)工具. 但小波譜的能量在頻率范圍內(nèi)分布較寬, 并存在選擇小波基和確定分解層數(shù)的問題. HHT方法中其EMD分解是完全自適應(yīng)于信號(hào)的, 具有多分辨分析和完全重構(gòu)的特性; 且取消了窗函數(shù)的作用, 具有完全的局部時(shí)頻特性. Hilbert時(shí)頻譜可準(zhǔn)確描述大地電磁信號(hào)的時(shí)變特征, 能更好地揭示電磁場(chǎng)信號(hào)的分布規(guī)律, 是一種更精確的MT信號(hào)分析方法.

[1] 王書明, 王家映. 大地電磁信號(hào)統(tǒng)計(jì)特征分析[J]. 地震學(xué)報(bào), 2004, 26(6): 669—674.

[2] 嚴(yán)家斌, 劉貴忠, 柳建新. 小波變換在天然電磁場(chǎng)信號(hào)時(shí)間序列處理中的應(yīng)用[J]. 地質(zhì)與勘探, 2008, 44(3): 75—78.

[3] CAI Jianhua, TANG Jingtian, HUA Xirui. An analysis method for magnetotelluric data based on the Hilbert-Huang Transform [J]. Exploration Geophysics, 2009, 40: 197—205.

[4] 湯井田, 化希瑞, 曹哲民, 等. Hilbert-Huang變換與大地電磁噪聲壓制[J]. 地球物理學(xué)報(bào), 2008, 51(2): 603—610.

[5] Donoho D L. De-noising by soft-thresholding [J]. IEEE Transaction on Information Theory, 1995, 41(3):613-627.

[6] Donoho D L, Johnstone I M. Adapting to unknown smoothness via wavelet shrinkage [J]. J American Statist Assoc, 1995, 90: 1200—1224.

[7] Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis [J]. Proc R Soc Lond Ser A, 1998, 454: 93—95.

[8] Huang N E, Wu M L, Long S R, et al. A confidence limit for the empirical mode decomposition and Hilbert spectral analysis [J]. Proceeding of Royal Society London A, 2003, 459: 2317—2345.

Comparison of wavelet analysis and Hilbert-Huang transform in MT signal procession

CAI JianHua1,HU WeiWen1, QING YeGui2, CHEN QinYe1

(1. Department of Physics and Electronics, Hunan University of Arts and Science, Changde 415000, China; 2. Chinese People's Liberation Army 93062 Troops, Jilin 132102, China)

Measured magnetotelluric signal is taken as the analysis object, and comparative study of wavelet analysis and Hilbert-Huang transform is applied in magnetotelluric (MT) signal procession from signal decomposition to time-frequency analysis. The mathematical theory of Wavelet analysis and Hilbert-Huang transform are expressed. The decomposition method of two ways and their physical significance for MT signal are analyzed, and the similarities and differences to characterize the spectrum of MT signal are discussed between the two methods. These provide a valuable reference foe the choice of analysis method of magnetotelluric signal procession.

magnetotelluric signal; Hilbert-Huang transform; wavelet analysis; decomposition; time-frequency analysis

P 631

1672-6146(2014)02-0029-06

10.3969/j.issn.1672-6146.2014.02.007

通訊作者email: cjh1021cjh@163.com.

2014-05-04

國家自然科學(xué)基金項(xiàng)目(41304098), 湖南省自然科學(xué)基金項(xiàng)目(12JJ4034); 湖南省教育廳青年項(xiàng)目(13B076); 湖南文理學(xué)院重點(diǎn)學(xué)科建設(shè)項(xiàng)目資助-無線電物理; 湖南省重點(diǎn)實(shí)驗(yàn)室“光電信息集成與光學(xué)制造技術(shù)”; “湖南省光電信息技術(shù)校企聯(lián)合人才培養(yǎng)基地”; 常德市科技局博士創(chuàng)新創(chuàng)業(yè)項(xiàng)目(2013BS01); 常德市科技局項(xiàng)目(2011JC02).

(責(zé)任編校:劉剛毅)