地震信號時(shí)頻分析中的希爾伯特黃變換研究

周竹生， 羅勇濤

(中南大學(xué)　地球科學(xué)與信息物理學(xué)院，長沙 410083)

?

地震信號時(shí)頻分析中的希爾伯特黃變換研究

周竹生， 羅勇濤*

(中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083)

摘要：隨著時(shí)頻分析方法的發(fā)展，生產(chǎn)研究上對復(fù)雜信號的時(shí)頻分析有了更高的要求。這里簡要介紹了希爾伯特—黃變換的原理和實(shí)現(xiàn)步驟，然后對合成信號進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)和希爾伯特譜分析，進(jìn)而將HHT方法運(yùn)用到實(shí)際地震信號的時(shí)頻分析中。在此基礎(chǔ)上，運(yùn)用經(jīng)驗(yàn)?zāi)B(tài)分解對合成地震信號進(jìn)行了閥值去噪，證明了該去噪方法的有效性。

關(guān)鍵詞：希爾伯特—黃變換； 經(jīng)驗(yàn)?zāi)B(tài)分解； 地震信號； 時(shí)頻分析； 時(shí)頻濾波

0引言

在地震勘探中，由于地下介質(zhì)的復(fù)雜性以及采集過程中的各種因素的干擾，獲得的地震信號具有非線性、非平穩(wěn)等特征，是一種典型的非平穩(wěn)隨機(jī)信號[1-2]。在以往的地震資料分析處理等研究中，傅立葉變換和功率譜估計(jì)都是基于平穩(wěn)信號分析處理理論建立起來的，只能得到整個(gè)信號的頻率成分，不能同時(shí)得到信號的時(shí)間和頻率分辨率。短時(shí)傅立葉變換、小波變換和Cohen類等時(shí)頻分析方法，雖然可以得到信號在不同時(shí)間、不同頻率處的能量密度和強(qiáng)度，但由于這類方法還是基于傅立葉分析理論，因而也受到傅立葉分析不足的制約。1998年Norden. E. Huang等[3]發(fā)表了希爾伯特黃變換(Hilbert-Huang Transform，HHT)的相關(guān)論文。該方法主要分兩步：①自適應(yīng)的經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)，得到有限個(gè)固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)；②進(jìn)行希爾伯特變換，得到信號的瞬時(shí)頻率和振幅，進(jìn)而得到希爾伯特譜。

地震信號作為一種復(fù)雜的非線性信號，其時(shí)頻分析方法自然也可以采用HHT方法，這里將HHT方法應(yīng)用到實(shí)際地震信號的分析中，討論其應(yīng)用效果。在此基礎(chǔ)上將經(jīng)驗(yàn)?zāi)B(tài)分解閥值去噪方法，運(yùn)用到合成地震道集的去噪中，驗(yàn)證了該去噪方法的有效性。

1HHT方法的原理

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

在進(jìn)行HHT方法的時(shí)頻分析之前，需要定義固有模態(tài)函數(shù)(IMF)。它必須滿足兩個(gè)方面的要求，①信號曲線過零點(diǎn)的數(shù)值和曲線的極值點(diǎn)數(shù)不能相差超過一個(gè)；②在任意一個(gè)局部范圍內(nèi)，極大值的包絡(luò)線和極小值的包絡(luò)線是關(guān)于時(shí)間軸近似對稱的[3]。

經(jīng)驗(yàn)?zāi)B(tài)分解就是將原信號分解為有限個(gè)不同尺度的模態(tài)函數(shù)，具體流程是：

1)找出原始信號x(t)所有的極值點(diǎn)，并將其用三次樣條函數(shù)插值成為原數(shù)據(jù)序列的上包絡(luò)線bmax(t)和下包絡(luò)線bmin(t)，并求出其上包絡(luò)線和下包絡(luò)線的平均m1(t)。

(1)

2)將原數(shù)據(jù)序列減去平均包絡(luò)后，得到一個(gè)新的數(shù)據(jù)序列。

h1(t)=x(t)-m1(t)

(2)

判斷h1(t)是否為IMF函數(shù)，若不滿足IMF條件，將h1(t)看作新的x(t)，重復(fù)上述處理過程，直到h1(t)滿足條件時(shí)記c1(t)為IMF(1)：

c1(t)=h1(t)

(3)

然后令：

r(t)=x(t)-c1(t)

(4)

3)將看作新的x(t)，重復(fù)以上步驟1)和步驟2)，即可依次得到固有模態(tài)函數(shù)的各階分量IMF(k)，直到滿足給定的終止條件時(shí)篩選結(jié)束。終止條件為：

(5)

根據(jù)大量實(shí)踐表明，迭代閾值SD設(shè)為0.2～0.3比較合適[3]。

4)原始的數(shù)據(jù)序列即可由這些IMF分量以及一個(gè)均值或趨勢項(xiàng)r(t)表示。

(6)

也就是說分解得到的各分量可以通過簡單相加得到原信號，所以EMD分解得到的IMF分量都是原信號固有的模態(tài)特征。

1.2瞬時(shí)分量和希爾伯特時(shí)頻譜

將分解得到的每一階IMF分量進(jìn)行希爾伯特變換：

(7)

式中cpv代表柯西主值( Cauchy Principal Value)。

構(gòu)造解析信號:

(8)

于是得到幅值函數(shù)：

(9)

相位函數(shù)：

(10)

進(jìn)一步可以求出瞬時(shí)頻率：

(11)

根據(jù)幅值函數(shù)和瞬時(shí)頻率，定義原始信號為：

(12)

這種將信號幅度表示在時(shí)間和頻率的平面上的分布就稱為Hilbert時(shí)頻譜，簡稱為Hilbert譜，用表示為：

(13)

式中ωj(t)=ω,Hilbert邊際譜的定義為：

(14)

邊際譜表示每個(gè)頻率在時(shí)間上的累計(jì)幅值，在統(tǒng)計(jì)意義上就是將所有的數(shù)據(jù)進(jìn)行時(shí)間上的累加。

2EMD分解閥值去噪

常規(guī)的去噪方法(如F-K濾波、小波變換、廣義S變換等)，其理論都受傅立葉變換的限制，濾波效果都會(huì)受到濾波的通放帶的影響。希爾伯特-黃變換是以瞬時(shí)頻率為信號變化的基本量、特征模態(tài)為基本時(shí)域信號的一種新方法，不會(huì)像上述方法受傅立葉變換理論的約束，可以在全時(shí)段兼顧到頻率分辨率和時(shí)間分辨率，具有很好的局部特性和自適應(yīng)性，適用于非平穩(wěn)信號的濾波和去噪。

在實(shí)際運(yùn)用中，通常閥值是固定的：

(15)

式中：σ2為噪聲的方差；N為采樣的長度。但是通常獲得精確的噪聲方差是很困難的，所以在計(jì)算中，常常采用估算的方法對噪聲方差進(jìn)行計(jì)算。

(16)

式中：mae表示平均絕對誤差；ch是IMF1的高頻系數(shù)。這是因?yàn)閏h是分解得到的最精細(xì)的系數(shù)，主要就是噪聲系數(shù)。

3信號的EMD分解與時(shí)頻分析

3.1合成信號的EMD分解

參考文獻(xiàn)[6]中合成信號由一個(gè)調(diào)幅調(diào)頻信號和一個(gè)正弦波疊加而成，其中調(diào)幅調(diào)頻信號的基頻 為50 Hz，其調(diào)制頻率為20 Hz，正弦波的頻率為120 Hz。

圖1(a)所示為合成信號和分解得到的6階IMF時(shí)間域波形曲線。第一行表示原始合成信號；IMF1表示一階序列(c1(t))，是第一個(gè)分解出來的高頻信號，對應(yīng)于原始信號中的120 Hz正弦波；IMF2表示一階序列(c2(t))，是第二個(gè)分解出來的信號，對應(yīng)于原始信號中的50 Hz的調(diào)頻調(diào)幅信號；IMF3、IMF4、IMF5、IMF6和r6分別代表三階、四階、五階、六階序列(c3(t)、c4(t)、c5(t)、c6(t))和剩余殘量。根據(jù)EMD的分解原理，高頻信號應(yīng)該首先被分解出來，即120 Hz的正弦波應(yīng)該是IMF1;其次應(yīng)該是50 Hz的調(diào)幅調(diào)頻信號，后面得到的IMF分量就是信號的趨勢分量和殘余信號了。這也證明了該方法的正確性和有效性。

3.2HHT時(shí)頻變換和常規(guī)時(shí)頻分析方法

對信號進(jìn)行EMD分解后，將得到的所有IMF分量進(jìn)行希爾伯特變換，再綜合各階瞬時(shí)頻譜，就得到信號的HHT時(shí)頻譜(圖2)。

從圖2(a)可看到兩條紅色曲線，分別代表了50 Hz為基頻的調(diào)幅調(diào)頻信號和120 Hz的正弦波。根據(jù)時(shí)頻譜圖可以看到，能量幅值隨時(shí)間和頻率是動(dòng)態(tài)變化的，而且不同頻率的細(xì)微不同也能區(qū)分出來。圖2(b)所示為對應(yīng)的HHT 邊際譜，它是不同頻率在整個(gè)信號時(shí)間內(nèi)的能量積累。

圖3(a)和圖3(b)分別是短時(shí)傅立葉變換和平滑偽Wigner-Ville分布，可以看到，雖然這兩種方法都可以分辨出50Hz和120Hz，但不能刻畫出相應(yīng)頻率的波形變化特征。這也正是由于這兩種方法都是基于傅立葉變換原理，不能進(jìn)行局部的時(shí)頻分析。

根據(jù)上面合成信號的時(shí)頻分析，可以看到HHT方法相比傳統(tǒng)方法具有很好的優(yōu)越性，它不需要人為選擇基函數(shù)，可以自適應(yīng)地進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，得到的固有模態(tài)函數(shù)是原信號的固有信息，而且可以根據(jù)需要改變幅值，從而能夠描述信號的非平穩(wěn)程度。

3.3實(shí)測地震信號的EMD分解和時(shí)頻分析

圖1　合成信號的EMD分解及其頻譜Fig．1　EMD decomposition and frequency spectrum of signal and IMF components(a)信號及6階IMF時(shí)間域波形；(b)信號及6階IMF的頻譜

圖2　合成信號的HHT時(shí)頻譜Fig．2　HHT time-frequency spectrum of signal(a) 合成信號的二維HHT時(shí)頻譜； (b) 合成信號的 HHT邊際譜

圖3　合成信號的時(shí)頻分析Fig．3　Time-frequency analysis of composite signal(a) 合成信號的短時(shí)傅立葉變換；(b) 合成信號的平滑偽Wigner-Ville分布

圖4　實(shí)測地震信號Fig．4　Original seismic signal

實(shí)際地震信號產(chǎn)生的復(fù)雜性以及采集過程中的各種干擾原因，地震檢波器獲得的信號往往是非線性、非平穩(wěn)的。圖4是一個(gè)原始地震波形，采集時(shí)間間隔為0.001 s。通過HHT方法對該地震信號進(jìn)行時(shí)頻分析，以說明HHT方法的有效性和刻畫信號非線性的能力。

圖 5(a)所示為實(shí)測地震信號及經(jīng)過EMD分解后的8階固有模態(tài)函數(shù)(IMF)。由EMD的分解原理，可知各分量是按頻率從高到低排列的，其中高頻成分震動(dòng)明顯，可以用于初動(dòng)檢測；低頻分量可以看作是簡單的非線性變化，可以用來消除趨勢干擾和進(jìn)行相關(guān)的校正。圖 5(b)所示為信號與各階IMF的頻譜，可見隨著分解次數(shù)的增加，得到的IMF分量能量也逐漸減弱。由圖5可見，整個(gè)記錄的能量主要集中在30 Hz以內(nèi)的低頻成分，這是由于檢波器在地表附近，而震源在地下幾千米，經(jīng)過地下介質(zhì)的吸收衰減，只有低頻信號能被檢波器接收到。

圖5　地震信號的EMD分解及其頻譜Fig．5　EMD decomposition and frequency spectrum of seismic signal and IMF components(a) 地震信號和8階IMF時(shí)間域波形；(b)地震信號和8階IMF的頻譜

圖7　地震數(shù)據(jù)的去噪Fig．7　The seismic data denoising(a)原始單炮記錄；(b)F-K濾波結(jié)果；(c)EMD分解閥值去噪結(jié)果

圖6是對圖5(a)所示的所有IMF分量進(jìn)行Hilbert變換得到的時(shí)頻譜。圖6(a)中的點(diǎn)、線表示能量變換特征的時(shí)頻演化過程，線表示能量的聚集性較好，且顏色越明顯，能量越大，反之，能量越??；點(diǎn)表示能量相對分散。由此可見，希爾伯特譜能很好地刻畫出信號幅值在時(shí)間和頻率上的變化情況，地震信號的非線性、非平穩(wěn)特征都能很好的體現(xiàn)出來，在細(xì)微的變化上，能準(zhǔn)確地表現(xiàn)出突變點(diǎn)的位置。圖6(b)為地震信號的HHT邊際譜，對應(yīng)圖6(a)中每個(gè)頻率在信號持續(xù)時(shí)間內(nèi)的能量積累，顯然能量集中在0 Hz～30 Hz頻段內(nèi)，這與實(shí)際的采集情況是相吻合的。

由上述分析，說明HHT方法具有完全的自適應(yīng)性，避免了人為選擇基函數(shù)和窗函數(shù)的影響，而且其時(shí)頻分辨率不受Heisenberg測不準(zhǔn)原理的限制，具有完全的局部時(shí)頻分析能力[8-10]。

4EMD分解閥值去噪效果測試

在地震資料的去噪中，需要對每一道地震信號進(jìn)行EMD分解閥值去噪，然后將處理后的地震道進(jìn)行疊加輸出，就可以得到多道的單炮記錄或者疊加剖面。圖7(a)為實(shí)際的單炮記錄，可以看到地震信號受面波干擾嚴(yán)重，其反射波與面波混疊，反射波的同相軸分辨困難，地震信號的信噪比較低。

相較于F-K濾波，EMD閥值去噪不僅壓制了面波，加強(qiáng)了淺層信號的強(qiáng)度，使同相軸更加清晰，而且各個(gè)時(shí)頻段的能量均沒有大的損失，在單炮記錄上每一道數(shù)據(jù)均能很好地顯示。

5結(jié)論

1)運(yùn)用EMD方法可以對信號進(jìn)行自適應(yīng)的分解，而且分解得到的每一個(gè)IMF分量都代表了原始信號的固有特征，因此可以選擇不同的IMF進(jìn)行時(shí)頻分析，具有很好的靈活性。

2)HHT時(shí)頻譜在地震信號的時(shí)頻分析中，能夠反映地震信號隨時(shí)間和時(shí)頻動(dòng)態(tài)變化的主要特征，且時(shí)頻分辨率不受Heisenberg測不準(zhǔn)原理的限制，這為地震信號的識別、分析提供了有利條件。

3)HHT 時(shí)頻譜的能量主要集中在一定的時(shí)間和頻率范圍內(nèi)，而不是整個(gè)時(shí)頻空間，這也真實(shí)反映了非線性、非平穩(wěn)序列的本質(zhì)特征。

4)HHT變換能夠很好地反映出信號的局部特性，也就能夠很好地反映出噪聲的頻率成分和范圍，因此運(yùn)用經(jīng)驗(yàn)?zāi)B(tài)分解閥值能夠很好地對信號進(jìn)行時(shí)頻濾波，去噪效果也比較理想。

[1]鄭治真，胡勁波.瞬時(shí)頻譜分析及其在地震趨勢估計(jì)中的應(yīng)用[J].地震學(xué)報(bào)，1993，15(l)：68-75.

ZHENG Z Z，HU J B.Instaneous spectrum analysis and its application in estimation of earthquake tendency[J].Acta Seismologica Sinica,1993，15(l)：68-75.(In Chinese)

[2]劉希強(qiáng)，周蕙蘭，李 紅. 基于小波包變換的地震數(shù)據(jù)時(shí)頻分析方法[J].西北地震學(xué)報(bào)，2000，22( 2)：143-147.

LIU X Q，ZHOU H L，LI H.The time-frequency analysis method about seismic data based on wavelet packet transform[J].North-Western Seismological Journal,2000，22( 2)：143-147.(In Chinese)

[3]HUANG N E, SHEN Z, LONG S R, et al.The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proc.R.Soc.London,1998,454:903-995

[4]WU Z，HUANG N E．A study of the characteristics of white noise using the empirical mode decomposition method[J]．Proceedings of the Royal Society A，2004，460(2046)：1597-1611．

[5]Wu Z H，HUANG N E．Ensemble empirical mode decomposition：a noise-assisted data analysis method[J]Advances in Adaptive Data Analysis,2009，1(1)：1-41．

[6]武安緒,吳培稚,蘭從欣,等.Hilbert-Huang 變換與地震信號的時(shí)頻分析研究[J].中國地震，2005，21(2)：207-215.

WU A X, WU P Z, LAN C X, et al. Hilbert-Huang transform and time—frequency analysis of seismic signal[J].Earthquake Research in China, 2005, 21(2): 207-216.(In Chinese)

[7]湯井田,蔡劍華,任政勇,等.Hilbert-Huang 變換與大地電磁信號的時(shí)頻分析研[J].中南大學(xué)學(xué)報(bào)，2009，40(5)：1399-1406.

TANG J T，CAI J H，REN Z Y，et al.Hilbert-Huang transform and time-frequency analysis of magnetotelluric signal[J].Journal of Central South University.2009，40(5)：1399-1406.(In Chinese)

[8]陳斌.時(shí)頻分析及在地震信號分析中的應(yīng)用研究[D].成都：成都理工大學(xué),2007.

CHENG B.Time-frequency analysis and its implication in seismic signal[D].Chengdu: Chengdu University of Technology,2007.(In Chinese)

[9]楊培杰,印興耀,張廣智.希爾伯特黃變換地震信號時(shí)頻分析與屬性提取[J].地球物理學(xué)進(jìn)展,2007,22(5):1585-1590.

YANG P J，YIN X Y，ZHANG G Z.Seismic signal time-frequency analysis and attributes extraction based on HHT[J].Progress in geophysics,2007,22(5):1585-1590.(In Chinese)

[10]黃光南.希爾伯特黃變換及其在地震資料分析處理中的應(yīng)用[D].青島：中國海洋大學(xué),2009.

HUANG G N.Hilbert Huang transform and applications of it in seismic data analysis and processing[D].Qingdao: Ocean University of China,2009.(In Chinese)

Research on Hilbert Huang transform in time-frequency analysis of seismic signals

ZHOU Zhu-sheng， LUO Yong-tao*

(School of Geoscience and Info-Physics,Central South University,Changsha410083,China)

Abstract:With the development of time-frequency analysis, higher feasibility and accuracy on production has been required for complex signal. Firstly, the time-frequency analysis method of Hilbert-Huang transform and conception of transient frequency are introduced. The empirical mode decomposition and time?frequency distribution of a compound signal is carried on and then the actual seismic data were analyzed and processed based on Hilbert-Huang transform (HHT) spectrum. A time-frequency filtering method based on the transform is used by tasting on a synthetic seismic records example. The result proves its effectiveness.

Key words:Hilbert-Huang transform; empirical mode decomposition; seismic signals; time-frequency analysis； time-frequency filtering

中圖分類號：P 631.4

文獻(xiàn)標(biāo)志碼：A

DOI：10.3969/j.issn.1001-1749.2016.01.09

文章編號：1001-1749(2016)01-0059-08

作者簡介：周竹生(1965-)，男，博士，教授，主要從事資源勘察、工程地球物理勘探、應(yīng)用軟件研制及數(shù)據(jù)庫開發(fā)等方面的教學(xué)和研究工作，E-mail：672390136@qq.com。*通信作者：羅勇濤(1990-)，男，碩士，主要從事工程地震勘探的研究，E-mail：1206387954@qq.com。

收稿日期：2015-01-16改回日期：2015-04-13