孫 苗,楊鈞凱,吳 立

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

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

以煙臺(tái)某地下洞室爆破地震波監(jiān)測信號為研究對象,通過進(jìn)行STFT、CWT、DWT、HHT 和改進(jìn)HHT的爆破地震波監(jiān)測信號時(shí)頻對比分析,最終發(fā)現(xiàn)改進(jìn)HHT算法是一種更具自適應(yīng)的爆破地震波信號處理算法,能有效抑制HHT在處理含噪地震波信號出現(xiàn)的模態(tài)混淆和端點(diǎn)發(fā)散的現(xiàn)象。改進(jìn)HHT算法能有效、準(zhǔn)確地對爆破地震波監(jiān)測信號進(jìn)行時(shí)頻分析,提取爆破地震波信號蘊(yùn)含的時(shí)頻特征參數(shù),為研究爆破振動(dòng)有害效應(yīng)提供了重要的依據(jù)。

1 多種時(shí)頻分析方法原理分析

1.1 短時(shí)傅里葉變換

短時(shí)傅里葉變換(STFT)[7]是在傳統(tǒng)傅里葉變換的基礎(chǔ)上通過加窗處理,實(shí)現(xiàn)可局部反應(yīng)信號時(shí)頻域信息,詳見式(1)。

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

(1)

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

1.2 連續(xù)小波變換

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

(2)

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

1.3 離散小波變換

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

(3)

1.4 HHT算法

HHT包含兩部分,第一部分是Huang N E提出的經(jīng)驗(yàn)?zāi)B(tài)分解(EMD);第二部分是Hilbert變換。HHT將時(shí)間信號通過EMD得到一組固有模態(tài)函數(shù)(IMF),再對IMF進(jìn)行Hilbert變換,得到Hilbert譜,即可描繪出信號的時(shí)頻譜和幅值譜。

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

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

圖1 EMD運(yùn)算流程Fig.1 EMD operation flow

1.4.2 Hilbert變換

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

(4)

(5)

1.5 改進(jìn)HHT算法

1.5.1 EMD模態(tài)混淆抑制研究

考慮到爆破地震波監(jiān)測信號多為含噪信號,噪聲的存在導(dǎo)致EMD得到的IMF產(chǎn)生嚴(yán)重的模態(tài)混淆,因此解決EMD模態(tài)混淆最根本的途徑就是對爆破地震波進(jìn)行降噪處理。

對EMD進(jìn)行改進(jìn)得到補(bǔ)充集合經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMD)[11],CEEMD是在原始含噪監(jiān)測信號S(t)中添加2個(gè)方向相反的噪聲信號,并分別(S1(t)=S(t)+正方向隨機(jī)噪聲,S2(t)=S(t)+反方向隨機(jī)噪聲)進(jìn)行EMD,即成對呈相反方向增加隨機(jī)噪聲。

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

1.5.2 端點(diǎn)處理

由EMD運(yùn)算原理可知,IMF產(chǎn)生的原理是不斷求均值的過程,均值來自于極大值和極小值包絡(luò)之差,由于端點(diǎn)不可能同時(shí)處于極大值和極小值,導(dǎo)致EMD在端點(diǎn)處發(fā)散,這種發(fā)散會(huì)從端點(diǎn)處向信號中間擴(kuò)散,數(shù)據(jù)越短,影響越大。

端點(diǎn)處理最直接的方法就是將原信號端點(diǎn)進(jìn)行延拓,使得原端點(diǎn)向中間移動(dòng),進(jìn)而使EMD免受端點(diǎn)效應(yīng)的影響。具體實(shí)現(xiàn)如下:

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

情況1:tmax1

(6)

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

(7)

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

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

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

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

2.1 端點(diǎn)處理爆破地震波信號STFT 時(shí)頻分析

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

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

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

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

2.2 爆破地震波信號CWT時(shí)頻分析

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

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

2.3 爆破地震波信號DWT時(shí)頻分析

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

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

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

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

2.4 爆破地震波信號HHT時(shí)頻分析

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

圖8 實(shí)測爆破振動(dòng)信號EMD結(jié)果Fig.8 EMD results of measured blasting vibration signal

不難發(fā)現(xiàn)由于噪聲的混入,導(dǎo)致IMF4和IMF5在采樣點(diǎn)500～700區(qū)間內(nèi)出現(xiàn)了向低頻發(fā)展的趨勢;IMF6在采樣點(diǎn)750～1 150區(qū)間內(nèi)出現(xiàn)了向低頻發(fā)展的趨勢,即中高頻有向中低頻混淆的趨勢,這種變化趨勢對時(shí)頻分析的準(zhǔn)確性影響很大。

進(jìn)一步分析發(fā)現(xiàn)IMF4、IMF5和IMF6在左端點(diǎn)出現(xiàn)了向低頻發(fā)散的現(xiàn)象;IMF8和IMF9在信號的左右端點(diǎn)都存在發(fā)散現(xiàn)象,端點(diǎn)發(fā)散也會(huì)對IMF真實(shí)性產(chǎn)生不利影響,進(jìn)而影響HHT時(shí)頻分析精度。

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

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

圖9直觀展示了噪聲信號以及端點(diǎn)效應(yīng)對信號整體時(shí)頻分布的影響,不僅降低了時(shí)頻分辨率,同時(shí)導(dǎo)致實(shí)際的時(shí)頻信息真實(shí)性受損,時(shí)頻分布欠穩(wěn)定。

2.5 爆破地震波信號改進(jìn)HHT時(shí)頻分析

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

圖10 基于改進(jìn)HHT的時(shí)頻能量三維圖Fig.10 Three dimensional diagram of time-frequency energy based on improved HHT

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

3 結(jié)論

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

2)CWT中a和τ的變化是連續(xù)的,導(dǎo)致計(jì)算量大、計(jì)算難度大,研究結(jié)果引入了大量的冗余成分,導(dǎo)致小波逆變換重構(gòu)不唯一;

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

4)EMD分解具有唯一性,噪聲的存在及端點(diǎn)效應(yīng)均會(huì)影響IMF的穩(wěn)定性和真實(shí)性,導(dǎo)致Hilbert變換后得到的時(shí)頻分析結(jié)果精度不高。改進(jìn)HHT算法通過抑制模態(tài)混淆和端點(diǎn)效應(yīng),可得到高分辨率時(shí)頻分析結(jié)果,對爆破振動(dòng)特性研究和爆破危害控制具有重要的現(xiàn)實(shí)意義。