吳 靜，吳 立，孫 苗，路亞妮，韓燕華

(1.湖北工程學(xué)院 土木工程學(xué)院，孝感 432000；2.中國地質(zhì)大學(xué)(武漢)工程學(xué)院，武漢 430074；3.湖北國土資源職業(yè)學(xué)院 環(huán)境與工程學(xué)院，武漢 430090)

經(jīng)驗(yàn)?zāi)B(tài)分解[1](Empirical Mode Decomposition，EMD)是一種根據(jù)爆破地震波信號本身特性進(jìn)行分解的自適應(yīng)算法，在爆破地震波信號時(shí)頻分析中得到了廣泛的應(yīng)用。解決爆破地震波信號在分解時(shí)產(chǎn)生的端點(diǎn)效應(yīng)問題，可有效提高EMD的分解精度[2]，得到反映真實(shí)爆破特性的固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)[3]。

爆破地震波信號通過EMD得到IMF需要進(jìn)行多次篩選，篩選的本質(zhì)是根據(jù)信號全部極大值點(diǎn)確定的上包絡(luò)線，和全部極小值點(diǎn)確定的下包絡(luò)線計(jì)算信號的局部均值。但信號的端點(diǎn)不可能同時(shí)處于極大值或極小值，且存在信號的端點(diǎn)不一定是極值點(diǎn)的情況，因此上、下包絡(luò)線在端點(diǎn)處會產(chǎn)生發(fā)散的現(xiàn)象，使得EMD的結(jié)果失真[4-6]。一維信號存在開始和結(jié)束節(jié)點(diǎn)；二維圖像存在邊界，無論是一維信號的開始和結(jié)束節(jié)點(diǎn)還是二維圖形的邊界，在進(jìn)行處理的時(shí)候都會受到影響，因此可以說絕大多數(shù)算法都存在端點(diǎn)效應(yīng)[7]。為了緩解EMD固有的端點(diǎn)效應(yīng)問題，有學(xué)者們提出了極值延拓法[8]、多項(xiàng)式擬合法以及端點(diǎn)鏡像法等方法[9,10]，上述方法多關(guān)注信號端點(diǎn)處的局部變化趨勢，欠缺對信號全局進(jìn)行考慮；也有學(xué)者提出考慮信號全局關(guān)聯(lián)性而弱化信號端點(diǎn)變化趨勢的方法[11]。

針對以上研究現(xiàn)狀，本文提出了一種基于邊界局部特征尺度自適應(yīng)匹配延拓(BLCS-AME)的EMD端點(diǎn)效應(yīng)抑制法，該方法可同時(shí)考慮信號在端點(diǎn)處的局部變化趨勢和原始信號本身的全局特性。該方法不僅能有效地抑制端點(diǎn)效應(yīng)，同時(shí)可實(shí)現(xiàn)信號特征參數(shù)的精確提取。

1 邊界局部特征尺度自適應(yīng)匹配延拓法

BLCS-AME方法由兩步組成：第一步，進(jìn)行邊界局部特征尺度延拓，即考慮信號端點(diǎn)幅值變化趨勢，同時(shí)考慮信號全局時(shí)間和端點(diǎn)時(shí)間間隔的內(nèi)在聯(lián)系對信號端點(diǎn)進(jìn)行延拓；第二步，進(jìn)行自適應(yīng)匹配，即將延拓好的局部特征尺度在信號全局找到與之匹配度最高的時(shí)間序列。

1.1 邊界局部特征尺度延拓

信號的變化趨勢在端點(diǎn)處和在信號內(nèi)部均有體現(xiàn)，信號全局時(shí)間和端點(diǎn)時(shí)間間隔之間存在內(nèi)在聯(lián)系。根據(jù)這種聯(lián)系，在信號左右端點(diǎn)各求出一個(gè)極大值點(diǎn)和極小值點(diǎn)對應(yīng)的時(shí)間參數(shù)，將時(shí)間參數(shù)導(dǎo)入到由端點(diǎn)極大(小)值點(diǎn)建立的多項(xiàng)式中，可求極大(小)值點(diǎn)對應(yīng)的幅值參數(shù)。

(1)邊界局部特征尺度延拓時(shí)間參數(shù)

以信號左端點(diǎn)為例，找到信號所有的極大值點(diǎn)發(fā)生的時(shí)刻，記為tmax 1,tmax 2,…,tmax i(i=1,2,3,…,M)，同理找到信號所有的極小值點(diǎn)發(fā)生的時(shí)刻，記為tmin 1,tmin 2,…,tmin i(i=1,2,3,…,N)。需要延拓的極大值點(diǎn)和極小值點(diǎn)的時(shí)刻分別記為tmax 0和tmin 0，其計(jì)算分以下4種情況：

情況1：tmax 1

(1)

(2)

情況2：tmax 1tmin N，M=N+1，tmax 0求解同式(2)，tmin 0求解見式(3)。

(3)

情況3：tmax 1>tmin 1∩tmax M>tmin N，M=N，tmax 0和tmin 0求解分別見式(4)和式(5)。

(4)

(5)

情況4：tmax 1>tmin 1∩tmax M

(6)

(2)邊界局部特征尺度延拓幅值參數(shù)

以信號左端點(diǎn)為例，找到全局信號中所有的極大值點(diǎn)發(fā)生時(shí)刻所對應(yīng)的幅值，記為xmax 1,xmax 2,…,xmax i(i=1,2,3…M)，同理找到全局信號中所有的極小值點(diǎn)發(fā)生時(shí)刻所對應(yīng)的幅值，記為xmin 1,xmin 2,…,xmin i(i=1,2,3…N)，根據(jù)端點(diǎn)附近幅值變化趨勢延拓出左端點(diǎn)附近的一個(gè)極值大值點(diǎn)和一個(gè)極小值點(diǎn)。

邊界局部特征尺度延拓具體步驟如下：

第一步：在信號的左端點(diǎn)附近，取最靠近左端點(diǎn)的a個(gè)極值點(diǎn)，即xmax 1,xmax 2,xmax 3,...,xmax a，a的取值和樣本容量有關(guān)；

第二步：對所取出的xmax 1,xmax 2,xmax 3,...,xmax a擬合出多項(xiàng)式，帶入tmax 0求出xmax 0，實(shí)現(xiàn)左端點(diǎn)極大值延拓，即(tmax 0,xmax 0)。

同上，可以實(shí)現(xiàn)左端點(diǎn)極小值點(diǎn)延拓，即(tmin 0,xmin 0)。

1.2 自適應(yīng)匹配延拓

信號的變化趨勢在端點(diǎn)和信號內(nèi)部均有體現(xiàn)，因此在信號內(nèi)部尋找一段和邊界局部特征尺度延拓結(jié)果匹配度最高的曲線是能夠?qū)崿F(xiàn)的。

自適應(yīng)匹配延拓具體步驟如下：

第一步：取1.1節(jié)得到左端點(diǎn)延拓結(jié)果(tmax 0,xmax 0)和(tmin 0,xmin 0)，以及與其最靠近的兩對極大(小)值點(diǎn)，即(tmax 1,xmax 1)、(tmin 1,xmin 1)和(tmax 2,xmax 2)和(tmin 2,xmin 2)，一共包括6個(gè)極值點(diǎn)的一段時(shí)間序列，記做l0，其包含k個(gè)采樣點(diǎn)；

第二步：將原始信號分成n個(gè)采樣點(diǎn)數(shù)為k的時(shí)間序列，記做l1,l2,l3,...,ln；

第三步：計(jì)算n個(gè)時(shí)間序列的自適應(yīng)匹配系數(shù)(Adaptive matching，ξam)，ξam計(jì)算見式(7)。其中，A=max{xmax 0,xmax 1,xmax 2}，B=max{xmin 0,xmin 1,xmin 2}。

(7)

第四步：篩選出一個(gè)和l0匹配度最高的li(i=0,1,2,…,n)。求自適應(yīng)匹配系數(shù)ξam的最小值ξam min=min{ξam1,ξam2,ξam3,…,ξamn}。ξam min所對應(yīng)的li即為所求。

第五步：將和l0匹配度最高的時(shí)間序列l(wèi)i平移到l0的位置，進(jìn)行EMD。

2 仿真信號多種端點(diǎn)效應(yīng)抑制方法對比研究

采用仿真信號進(jìn)行多種EMD端點(diǎn)效應(yīng)抑制對比研究，再結(jié)合EMD得到的IMF和原始信號之間的相關(guān)系數(shù)和誤差標(biāo)準(zhǔn)差進(jìn)行EMD端點(diǎn)效益抑制評價(jià)。

仿真信號S(t)由頻率為10 Hz、30 Hz、60 Hz的3個(gè)正弦信號組成，即S(t)=x1(t)+x2(t)+x3(t)，其中：x1(t)=sin(2×π×10×t)；x2(t)=sin(2×π×30×t)；x3(t)=sin(2×π×60×t)。采樣點(diǎn)N=200，采樣時(shí)間t分別取0,π/100,2π/100,3π/100,…,π，即在[0,2π]上，中間以π/100為間距取點(diǎn)，仿真信號如圖1(a)所示。

對S(t)進(jìn)行未采取端點(diǎn)效應(yīng)抑制措施的EMD，得到的IMF如圖1(b)，可發(fā)現(xiàn)IMF1、IMF2和IMF3分量的兩端都出現(xiàn)了不同程度的發(fā)散，隨著分解的進(jìn)行IMF3的發(fā)散最嚴(yán)重，并有向數(shù)據(jù)內(nèi)部發(fā)展的趨勢。采用邊界局部特征尺度自適應(yīng)匹配延拓法、極值延拓法、多項(xiàng)式擬合法得到的分解結(jié)果分別見圖1(c)、圖1(d)、圖1(e)，考慮到圖片的分辨率以及清晰度問題，這里僅展示不同端點(diǎn)效應(yīng)抑制方法得到的IMF在左端點(diǎn)處的變化趨勢。

圖 1 仿真信號分解結(jié)果

進(jìn)一步計(jì)算分量和對應(yīng)正弦信號之間的相關(guān)系數(shù)(rxy)和誤差標(biāo)準(zhǔn)差(Dsde)，計(jì)算公式見式(8)和式(9)。式中，N為采樣點(diǎn)數(shù)，x對應(yīng)x1(t)，x2(t)和x3(t)，y對應(yīng)IMF1、IMF2和IMF3。

(8)

(9)

邊界局部特征尺度自適應(yīng)匹配延拓法、極值延拓法、多項(xiàng)式擬合法的相關(guān)系數(shù)(rxy)和誤差標(biāo)準(zhǔn)差(Dsde)的計(jì)算結(jié)果見表1。

表 1 端點(diǎn)效應(yīng)抑制方法的評價(jià)指標(biāo)

由圖1和表1可知，邊界局部特征尺度自適應(yīng)匹配延拓法得到的3個(gè)IMF分量反映了仿真信號所包含的3個(gè)正弦信號，與對應(yīng)正弦信號之間相關(guān)性最高、誤差最小。即邊界局部特征尺度自適應(yīng)匹配延拓法得到的IMF精度最高，對EMD端點(diǎn)效應(yīng)的抑制效果最好；極值延拓法的效果略優(yōu)于多項(xiàng)式擬合法，尤其是對低頻分量的抑制。

綜上得出結(jié)論：邊界局部特征尺度自適應(yīng)匹配延拓法(BLCS-AME)對EMD端點(diǎn)效應(yīng)具有較好的抑制能力，能夠得到高精度的IMF。

3 爆破地震波信號EMD端點(diǎn)效應(yīng)抑制應(yīng)用

3.1 工程背景

工程爆破的發(fā)展極大提高了工作效率，給我國基礎(chǔ)設(shè)施工程建設(shè)帶來巨大的便利，但其產(chǎn)生的地震效應(yīng)對周邊建筑物、構(gòu)筑物及周圍環(huán)境的影響也日益突出[12-14]，特別是在城市進(jìn)行的爆破施工，其影響更為顯著。以湖北某水下鉆孔爆破工程為對象，研究BLCS-AME方法在爆破地震波信號EMD端點(diǎn)效應(yīng)抑制中的應(yīng)用，得到反映爆破地震波真實(shí)屬性的IMF，分析爆破產(chǎn)生的地震效應(yīng)對爆破區(qū)周邊建筑物的影響。

某爆破區(qū)施工環(huán)境圖如圖2所示，爆破區(qū)位于水下10～15 m，位于某城市道路交叉口附近，其中，主路距離爆破區(qū)150～180 m，輔路距離爆破區(qū)240～280 m，在爆破區(qū)與輔路之間距離爆破區(qū)100～120 m處有一棟地上7層地下1層花園洋房。為了對控制爆破地震波的有害效應(yīng)提供依據(jù)，開展了現(xiàn)場爆破振動(dòng)監(jiān)測，將花園洋房作為本次爆破施工重點(diǎn)監(jiān)控對象。

圖 2 爆破施工環(huán)境圖(單位：m)

3.2 監(jiān)測結(jié)果及分析

現(xiàn)場爆破采用TC-4850型爆破測振儀對該樓棟進(jìn)行監(jiān)測，設(shè)置的信號采樣頻率為4000 sps。選取一條典型爆破振動(dòng)信號曲線作為研究對象，如圖3所示，測試時(shí)間為1.8 s。根據(jù)奈奎斯特采樣定理，實(shí)測爆破地震波信號的奈奎斯特頻率值為2000 Hz，共包括4096個(gè)采樣點(diǎn)。

對圖3所示信號進(jìn)行基于邊界局部特征尺度自適應(yīng)匹配延拓法的EMD分解，分解結(jié)果見圖4(a)，由圖4(a)可知，僅IMF4的右端點(diǎn)處出現(xiàn)了輕微的發(fā)散現(xiàn)象，其余分量的端點(diǎn)效應(yīng)都得到了很好的控制。

基于邊界局部特征尺度自適應(yīng)匹配延拓法得到的IMF，經(jīng)過Hilbert變換，進(jìn)一步分析求出每個(gè)IMF的邊際譜，如圖4(b)所示，由圖4(b)可知，基于邊界局部特征尺度自適應(yīng)匹配延拓法的EMD分解得到的IMF從高頻到低頻依次排列，每個(gè)IMF都攜帶爆破地震波信號一組特定的頻率信號。

圖 3 地震波監(jiān)測信號

圖 4 地震波監(jiān)測信號分解圖

進(jìn)一步計(jì)算，得到如圖5所示的信號總邊際譜，其縱坐標(biāo)ESD為能量譜密度(Energy spectral density，ESD)，由圖5可知水下鉆孔爆破地震波能量主要集中在0～50 Hz，這與文獻(xiàn)[15,16]得出的結(jié)論一致。信號的主頻為最大能量密度對應(yīng)的頻率[17]，因此本信號主頻為14.286 Hz。

圖 5 信號的總邊際譜

3.3 數(shù)值模擬分析

通過有限元軟件對該洋房進(jìn)行建模分析，得到洋房的三維模型圖如圖6所示，并計(jì)算該洋房的自振頻率。該洋房前8階陣型對應(yīng)的自振頻率見表2。

圖 6 洋房三維模型圖

表 2 洋房前8階陣型對應(yīng)的自振頻率(單位：Hz)

由表2可見該洋房的第2階陣型與本次爆破信號主頻十分接近。該結(jié)果表明，本工程若繼續(xù)采用目前的爆破參數(shù)，產(chǎn)生的爆破地震波極有可能引起該洋房共振，實(shí)際施工中必須采取相應(yīng)控制措施才能保證該洋房的安全。

上述案例分析結(jié)果證明：本文提出的邊界局部特征尺度自適應(yīng)匹配延拓法，不僅有助于抑制EMD產(chǎn)生的端點(diǎn)效應(yīng)，得到蘊(yùn)含爆破地震波特征參數(shù)的高精度IMF分量，而且有助于進(jìn)一步識別爆破地震波效應(yīng)，實(shí)現(xiàn)爆破地震波危害控制。

4 結(jié)論

(1)邊界局部特征尺度自適應(yīng)匹配延拓法(BLCS-AME)不僅考慮了信號在端點(diǎn)處的局部變化趨勢，同時(shí)通過波形匹配保留了信號原始波形的整體特征，是一種更具自適應(yīng)性的EMD端點(diǎn)效應(yīng)抑制算法。

(2)對比分析邊界局部特征尺度自適應(yīng)匹配延拓法(BLCS-AME)和其他常規(guī)端點(diǎn)效應(yīng)抑制方法對仿真信號的分解結(jié)果，發(fā)現(xiàn)BLCS-AME能夠在有效抑制EMD端點(diǎn)效應(yīng)的基礎(chǔ)上，實(shí)現(xiàn)仿真信號細(xì)節(jié)特征參數(shù)的提取，得到精度更高的IMF。

(3)基于邊界局部特征尺度自適應(yīng)匹配延拓法(BLCS-AME)得到的IMF經(jīng)過Hilbert變換得到的邊際譜，可有效提取爆破地震波所蘊(yùn)含的頻率-能量信息，該信息有助于爆破振動(dòng)特征的識別，為科學(xué)制定抗震措施提供依據(jù)。