基于改進HHT的地下洞室爆破地震波信號時頻分析

孫 苗,楊鈞凱,吳 立

(1.湖北國土資源職業學院,武漢 430090;2.中國地質大學(武漢)巖土鉆掘與防護教育部工程研究中心, 武漢 430074;3.中國地質大學(武漢)工程學院, 武漢 430074)

爆破地震波監測信號表現為瞬時、突變和震蕩特征,屬于典型的非平穩信號[1]。時頻分析已成為處理這類信號的重要手段[2]。常用的時頻分析方法[3-5]有:短時傅里葉變換(STFT)、連續小波變換(CWT) 、離散小波變換(DWT)等。上述時頻分析方法建立基礎均是傅里葉變換,因而不可避免地會受到傅里葉變換分析非平穩信號缺陷的影響,如出現虛假頻率和多余分量等。針對這一問題,Huang N E[6]等提出了HHT算法,該算法依據數據本身的特性進行分解,從根本上突破了傅里葉變換的限制。但是實際施工中爆破地震波信號監測環境復雜,導致監測信號不可避免混有噪聲,噪聲的存在使得EMD產生嚴重的模態混淆。同時信號存在開始和結束節點,導致絕大多數算法在端點處無法避免會產生端點效應。模態混淆和端點效應都是影響HHT時頻分析精度的主要原因,因此為了獲得準確的爆破地震波特征參數,必須對傳統的HHT算法進行模態混淆抑制和端點處理。

以煙臺某地下洞室爆破地震波監測信號為研究對象,通過進行STFT、CWT、DWT、HHT 和改進HHT的爆破地震波監測信號時頻對比分析,最終發現改進HHT算法是一種更具自適應的爆破地震波信號處理算法,能有效抑制HHT在處理含噪地震波信號出現的模態混淆和端點發散的現象。改進HHT算法能有效、準確地對爆破地震波監測信號進行時頻分析,提取爆破地震波信號蘊含的時頻特征參數,為研究爆破振動有害效應提供了重要的依據。

1 多種時頻分析方法原理分析

1.1 短時傅里葉變換

短時傅里葉變換(STFT)[7]是在傳統傅里葉變換的基礎上通過加窗處理,實現可局部反應信號時頻域信息,詳見式(1)。

F(τ,ω)=STFT{f(t)}=FT{f(t)·w(t-τ)}

(1)

式中:f(t)為原始信號;w(t-τ)是一個以時刻τ為中心的窗函數,其作用是對τ附近的函數做傅里葉變換,得到τ附近的頻率信息。

1.2 連續小波變換

將任意L2(R)空間中的原始信號f(t)在小波基函數下展開,稱這種展開為f(t)的連續小波變換(CWT)[8],其表達式為

(2)

式中:a為伸縮因子;τ為平移因子;ψa,τ(t)是依賴參數a,τ的小波基函數。

1.3 離散小波變換

連續小波變換中a和τ的變化是連續的,離散小波變換(DWT)[9]是不連續的,DWT定義式為

(3)

1.4 HHT算法

HHT包含兩部分,第一部分是Huang N E提出的經驗模態分解(EMD);第二部分是Hilbert變換。HHT將時間信號通過EMD得到一組固有模態函數(IMF),再對IMF進行Hilbert變換,得到Hilbert譜,即可描繪出信號的時頻譜和幅值譜。

1.4.1 經驗模態分解

EMD算法實現過程如圖1所示, S(t)為爆破地震波監測信號;Smax(t)和Smin(t) 為上、下包絡線;M(t)為Smax(t)和Smin(t)的均值;D(t)為S(t)和M(t)的差值;Ri(t)為余項。

圖1 EMD運算流程Fig.1 EMD operation flow

1.4.2 Hilbert變換

Hilbert變換[10]的實質是對任意信號f(t)和h(t)做卷積,將f(t)換成IMF可實現任意IMF的Hilbert變換。

(4)

(5)

1.5 改進HHT算法

1.5.1 EMD模態混淆抑制研究

考慮到爆破地震波監測信號多為含噪信號,噪聲的存在導致EMD得到的IMF產生嚴重的模態混淆,因此解決EMD模態混淆最根本的途徑就是對爆破地震波進行降噪處理。

對EMD進行改進得到補充集合經驗模態分解(CEEMD)[11],CEEMD是在原始含噪監測信號S(t)中添加2個方向相反的噪聲信號,并分別(S1(t)=S(t)+正方向隨機噪聲,S2(t)=S(t)+反方向隨機噪聲)進行EMD,即成對呈相反方向增加隨機噪聲。

CEEMD的具體步驟如下:①成對地添加方向相反的隨機噪聲到監測信號中,得到S1(t)和S2(t);②對S1(t)和S2(t)分別進行EMD;③重復步驟①、②直到達到預設的添加次數為止;④將步驟②得到的結果進行總體平均。

1.5.2 端點處理

由EMD運算原理可知,IMF產生的原理是不斷求均值的過程,均值來自于極大值和極小值包絡之差,由于端點不可能同時處于極大值和極小值,導致EMD在端點處發散,這種發散會從端點處向信號中間擴散,數據越短,影響越大。

端點處理最直接的方法就是將原信號端點進行延拓,使得原端點向中間移動,進而使EMD免受端點效應的影響。具體實現如下:

找到信號所有的極大值點的坐標,即(tmax1,xmax1),(tmax2,xmax2)…(tmaxi,xmaxi) (i=1,2,3…M),同理找到信號所有的極小值點的坐標,記為(tmin1,xmin1),(tmin2,xmin2)…(tminj,xminj) (j=1,2,3…N),需要延拓的極大值點和極小值點時刻分別為tmax0和tmin0,其計算分以下2種情況。

情況1:tmax1

(6)

情況2:tmax1>tmin1,tmax0和tmin0求解見式(7):

(7)

對所有極大值點的縱坐標進行多項式擬合,并將計算得到的tmax0代入擬合式,便可計算出xmax0,xmin0的計算和xmax0一致。假設滿足情況1,(tmax1,xmax1),(tmin0,xmin0)和(tmax0,xmax0)組成的三角波即為延拓的結果,對此時的信號進行CEEMD即可實現EMD模態混淆和端點抑制處理。對此時得到的IMF進行Hilbert變換即可實現改進HHT爆破地震波信號時頻分析。

2 多種方法爆破地震波信號時頻分析對比研究

以煙臺某地下水封LPG洞庫爆破施工工程為依托,采用TC-4850型爆破測振儀進行監測,選取一條典型的中間主洞室爆破地震波監測信號為研究對象,該地震波監測信號如圖2所示。該信號采樣頻率為4 000 Hz,在0～0.899 s可采集3 600個數據點。

圖2 原始爆破地震波監測信號Fig.2 Original blasting seismic wave monitoring signal

2.1 端點處理爆破地震波信號STFT 時頻分析

采用STFT[Hamming(749)]對圖2爆破地震波監測信號進行頻譜分析,得到如圖3所示的時頻譜和能量譜密度圖,圖3b中的ESD為能量譜密度。從圖3可知,當選定了窗函數,即確定了信號的時間分辨率。窗函數越窄,時域特征越明顯,在窗內進行快速傅里葉變換會因數據點過少導致快速傅里葉變換精度降低,即頻率分辨率降低,導致其在應用過程中只能滿足一種分辨率需求。

圖3 STFT[Hamming(749)]時頻譜和能量譜密度Fig.3 Spectrum and energy spectral density at STFT[Hamming (749)]

對比圖3b和圖4b可發現,Hamming由749縮小到257,頻率分辨率降低,時頻譜精度降低,快速傅里葉變換分析結果真實性需要進一步研究。因此STFT適合頻率波動不大的平穩信號,不適合爆破地震波這種非平穩、非線性的信號。

圖4 STFT[Hamming(257)]時頻譜和能量譜密度Fig.4 Spectrum and energy spectral density at STFT[Hamming (257)]

2.2 爆破地震波信號CWT時頻分析

圖2爆破地震波監測信號CWT[scale=1:32,db5]處理后,可得到各個小波系數能量占比,如圖5b所示,同時得到小波系數在時間-尺度平面上的分布,如圖5c所示。觀察圖5可發現,CWT在分析的過程中雖然可計算出特定尺度-位移平面下的小波系數,但同時也引入了大量的冗余成分,該特點導致小波逆變換重構不唯一。

圖5 CWT[scale=1:32,db5]小波譜Fig.5 CWT[scale = 1:32, db5] wavelet spectrum

2.3 爆破地震波信號DWT時頻分析

通過DWT[db5(分解8層)]對圖2爆破地震波監測信號進行分解,得到圖6所示的結果。觀察圖6b可發現,信號的高頻蘊含在d1～d3中,低頻蘊含在d7、d8和a8中,可發現DWT能夠較好地描述信號的時頻特征,可用于非平穩信號時頻分析。

圖6 DWT[db5,8層]分量Fig.6 DWT[db5, 8th floor] component diagram

研究發現DWT雖能夠在一定程度描述信號的局部特征,但會因小波基選擇不同導致不同的分解重構結果,說明小波分量依賴于小波基的選擇。圖7a是不同小波基函數小波分量d1的對比圖,通過該圖可發現db5小波基和sym3小波基得到的d1分量之間存在明顯的差距,這也從側面反映,DWT分解并不唯一。進一步觀察圖6b[db5(分解8層)]得到的小波系數圖和圖7c[sym3(分解8層)]得到的小波系數圖,也可發現明顯的差異,因此在小波變換中,小波基函數的選擇顯得極其重要。

圖7 DWT[sym3,8層/db5,8層]小波譜Fig.7 DWT[sym3, 8th floor / db5, 8th floor] wavelet spectrum

2.4 爆破地震波信號HHT時頻分析

針對小波變換受限于小波基選擇的問題,采用EMD對圖2信號進行分解,得到如圖8的分解結果。從圖8可見:分量IMF1-IMF2表現為頻率高、幅值低、能量低。而即噪聲信號一般常表現為高頻低能,可初步確定IMF1-IMF2為在實際監測中混入的噪聲成分。

圖8 實測爆破振動信號EMD結果Fig.8 EMD results of measured blasting vibration signal

不難發現由于噪聲的混入,導致IMF4和IMF5在采樣點500～700區間內出現了向低頻發展的趨勢;IMF6在采樣點750～1 150區間內出現了向低頻發展的趨勢,即中高頻有向中低頻混淆的趨勢,這種變化趨勢對時頻分析的準確性影響很大。

進一步分析發現IMF4、IMF5和IMF6在左端點出現了向低頻發散的現象;IMF8和IMF9在信號的左右端點都存在發散現象,端點發散也會對IMF真實性產生不利影響,進而影響HHT時頻分析精度。

對圖8得到的IMF進行Hilbert變換得到信號的時間-頻率-能量譜密度三維圖,如圖9所示。

圖9 基于HHT的時頻能量三維圖Fig.9 Three dimensional diagram of time-frequency energy based on HHT

圖9直觀展示了噪聲信號以及端點效應對信號整體時頻分布的影響,不僅降低了時頻分辨率,同時導致實際的時頻信息真實性受損,時頻分布欠穩定。

2.5 爆破地震波信號改進HHT時頻分析

通過改進HHT算法求得圖2爆破地震波監測信號的時頻能量三維圖如圖10所示。對比圖9和圖10可發現,通過抑制模態混淆和端點效應得到信號視頻能量三維圖,在時間和頻率上都具有很高的分辨率。本次爆破能量主要集中在150～300 Hz這一頻率范圍內,這一研究結論和李洪濤[12]關于地下洞室爆破地震波信號頻率能量分布研究結果一致。從側面反映出,改進HHT算法是一種更具自適應的時頻分析算法,運算結果具有可靠性,運算更穩定。

圖10 基于改進HHT的時頻能量三維圖Fig.10 Three dimensional diagram of time-frequency energy based on improved HHT

爆破地震波時頻分析是為了讓爆破科研人員掌握爆破產生的危害效應,并制定對應的爆破危害控制措施[13]。改進HHT算法可實現根據信號本身特征進行解析,同時可實現傳統HHT算法模態混淆和端點效應抑制研究,得到高分辨率的時間-頻率-能量三維圖,實現爆破地震波時頻參數提取,有助于爆破危害分析和控制。對爆破振動特性研究和爆破地震波衰減規律學習具有重要的現實意義。

3 結論

1)STFT窗函數選定,便確定了信號的時間分辨率。因此只適合平穩、線性信號時頻分析;

2)CWT中a和τ的變化是連續的,導致計算量大、計算難度大,研究結果引入了大量的冗余成分,導致小波逆變換重構不唯一;

3)DWT中a和τ的變化不連續,能在一定程度描述信號的局部特征,但是過度依賴于小波基的選擇,分解不唯一;

4)EMD分解具有唯一性,噪聲的存在及端點效應均會影響IMF的穩定性和真實性,導致Hilbert變換后得到的時頻分析結果精度不高。改進HHT算法通過抑制模態混淆和端點效應,可得到高分辨率時頻分析結果,對爆破振動特性研究和爆破危害控制具有重要的現實意義。