基于CEEMDAN的塌落觸地振動信號最優(yōu)降噪光滑模型算法*

趙明生，茍倩倩，張光雄，池恩安，胡娟新,姚穎康

(1.保利新聯(lián)爆破工程集團(tuán)有限公司，貴陽 550002；2.貴州大學(xué) 礦業(yè)學(xué)院，貴陽 550025； 3.保利民爆哈密有限公司，哈密 839000；4.武漢爆破有限公司，武漢 430000)

隨著安全環(huán)保意識的提高，對城鎮(zhèn)爆破拆除作業(yè)所帶來的塌落觸地振動危害控制越來越嚴(yán)格。在建(構(gòu))筑物倒塌過程中，炸藥釋放的能量使得部分構(gòu)件解體，造成上部結(jié)構(gòu)失去穩(wěn)定性從而塌落，觸地沖擊地面瞬間會產(chǎn)生很大的能量，其中，部分能量促使地面介質(zhì)產(chǎn)生振動，以地震波的形式向外傳播，并在地表產(chǎn)生震動[1]。由于受爆破地震波監(jiān)測點的隨機(jī)性、儀器誤差、磁場干擾、傳播介質(zhì)等因素影響，使爆破振動信號攜帶一定的噪聲，難以精確識別其模態(tài)信息，從而影響振動危害評價的準(zhǔn)確度。

目前，爆破振動信號降噪主要采取二代小波變換、集合經(jīng)驗?zāi)B(tài)分解及集合經(jīng)驗?zāi)B(tài)分解和小波閾值共同作用等方法[2-4]，小波變換分析法有效性很大程度上依賴小波基函數(shù)的選取，存在一定程度的主觀因素影響[5]。集成經(jīng)驗?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition,EEMD)是在原始信號中多次添加不同白噪聲，分別進(jìn)行EMD分解，然后對多次 EMD 分解的 IMF 進(jìn)行平均而得到最終的實際分量，它能有效地改善EMD方法所存在的模態(tài)混疊現(xiàn)象[6]。但是EEMD在信號重構(gòu)時存在完備性不足，且計算耗時長的問題[7-9]。近年來，由Torres提出的具有自適應(yīng)噪聲的完整集成經(jīng)驗?zāi)B(tài)分解[10](Complete Ensemble Empirical Mode Decomposition with Adaptive Noise,CEEMDAN)，通過在每個階段添加有限次的自適應(yīng)白噪聲，計算唯一余項信號獲取IMF，實現(xiàn)在較少的平均次數(shù)下，幾乎完全重構(gòu)分解后的固有模態(tài)，實現(xiàn)精確的信號特征提取。該方法能有效解決EMD分解的模態(tài)混疊問題，同時還能克服EEMD分解的完備性較差及多次增加集成次數(shù)來降低重構(gòu)誤差而導(dǎo)致的計算效率低的問題[11-13]。

基于此，提出一種基于CEEMDAN的爆破振動信號最優(yōu)光滑降噪算法，以設(shè)計仿真信號為研究對象，將其分別進(jìn)行EMD、EEMD及CEEMDAN的最優(yōu)光滑降噪算法效果對比，并將該算法應(yīng)用于實測水塔爆破觸地振動信號降噪分析。

1 建立最優(yōu)降噪光滑模型

1.1 設(shè)計濾波器

將原始含噪信號進(jìn)行CEEMDAN分解后的IMF分量進(jìn)行高通、低通、帶通濾波器的組合設(shè)計[14]。其中，高通濾波常用來于增強(qiáng)高音成分以驅(qū)動揚聲器的高音單元，低通濾波可以增強(qiáng)中音成分，而帶通濾波可以對特定頻率范圍內(nèi)的信號進(jìn)行篩選，達(dá)到選擇性傳輸?shù)哪康摹?/p>

高通濾波器可設(shè)計為

(1)

低通濾波器可設(shè)計為

(2)

帶通濾波器可設(shè)計為

(3)

式中：a、q為濾波器的上截斷參數(shù)；b、p為濾波器的下截斷參數(shù)，其取值應(yīng)根據(jù)不同信號的濾波要求來優(yōu)化確定。

假設(shè)由CEEMDAN分解出的固有模態(tài)分量IMF有N個，則具體的濾波器算法表達(dá)式為

A1=s(n)-IMF1

A2=s(n)-IMF1-IMF2

A3=s(n)-IMF1-IMF2-IMF3

?

An=s(n)-IMF1-IMF2-IMF3-…-IMFn

(4)

1.2 定義降噪偏均方誤差

(5)

由定義可知，降噪均方誤差越小，降噪后的曲線越接近原始信號，即兩個曲線相似度越高。

1.3 定義曲線曲折度

設(shè)有兩條曲線P(t)、q(t)(0≤t≤1)，如果曲線P在點P(1)處的二階倒數(shù)等于曲線Q在點Q(0)的二階倒數(shù)，則曲線光滑。

(6)

將點P(1)、Q(0)處的二階倒數(shù)按照離散公式展開得

(7)

(8)

由于曲線光滑，則點P(1)、Q(0)一階導(dǎo)數(shù)相等，即可得到在x=x0點處的光滑指標(biāo)

SNx=x0=f(x0+2h)-f(x0-2h)-

2[f(x0+h)-f(x0-h)]

(9)

式中：h為步長，取值1；當(dāng)x0趨于零時，x0附近的點對應(yīng)曲線越光滑，即曲線的曲折度越小，信號越光滑。

1.4 基于CEEMDAN最優(yōu)降噪光滑模型判定準(zhǔn)則

在模型的建立中，降噪偏差均方差MSEf的值越小，降噪后的信號與原始信號相似度越高，同時濾波算法的光滑性指標(biāo)SMSEf的值越小，曲線越光滑，即目標(biāo)函數(shù)值越小降噪效果越好?；诖?，建立約束條件

(10)

則爆破振動信號濾波算法目標(biāo)函數(shù)可以寫成

min{f}=min{αSMSEf+(1-α)MSEf}

(11)

式中，α為曲線降噪均方誤差權(quán)重因子，1-α為曲線曲折度權(quán)重因子，其取值按照信號處理過程中具體要求而定?；贑EEMDAN的最優(yōu)光滑降噪流程如圖1所示。

2 仿真信號降噪分析

設(shè)計仿真信號對上述降噪算法進(jìn)行驗證，仿真信號Y(t)分別由低頻正弦信號x1，調(diào)頻信號x2及功率為0.01的高斯白噪聲x3合成。

(12)

式中:t=[0,4]，采樣頻率為2000 Hz，各組成信號及混合仿真信號時域波形如圖2所示。

對仿真信號分別進(jìn)行EMD、EEMD及CEEMDAN分解。其中，EMD分解得到10個固有模態(tài)分量一個余項，EEMD分解得到12個固有模態(tài)分量及一個余項，CEEMDAN分解后得到13個固有模態(tài)分量及余項。按照式(4)EMD分解將會產(chǎn)生11個帶通濾波算法，同理，EEMD及CEEMDAN分解分別產(chǎn)生13和14個帶通濾波算法，將仿真信號分別用三種方法分解后組成的算法進(jìn)行降噪處理，求出每個算法的降噪均方誤差及光滑度，按照約束條件確定目標(biāo)函數(shù)值，并進(jìn)行歸一化處理，當(dāng)曲線降噪均方誤差權(quán)重因子取0.5時，三種方法分解后組成的算法的降噪誤差曲線如圖2～圖4所示。

結(jié)合圖2至圖4可以看出：當(dāng)降噪均方誤差權(quán)重因子取0.5時，基于EMD分解的11個算法中，A2算法目標(biāo)函數(shù)值最小，其值為0.1653；基于EEMD分解的13個算法中，A2算法目標(biāo)函數(shù)值最小，其值為0.3179；基于CEEMDAN分解的14個算法中，A4算法目標(biāo)函數(shù)值最小，其值為0.1353，即基于CEEMDAN分解的A4算法的目標(biāo)函數(shù)值最小，降噪效果最好。圖2至圖4中目標(biāo)函數(shù)曲線在算法A10之后幾乎都趨于平衡，這是由于從A10算法之后減去的固有模態(tài)分量都是低頻，而高頻噪聲在算法A10之前幾乎被完全過濾掉。

圖5為基于CEEMDAN分解的算法A4、A5、A6、A7降噪的仿真信號與原始仿真信號，由圖可以看出：經(jīng)算法A4降噪的信號在低頻部分明顯比原仿真信號光滑，并且重構(gòu)后的信號幾乎與原仿真信號重合，算法A5僅次于算法A4，而經(jīng)算法A6、A7降噪的信號在1 s之后頻率明顯較低，出現(xiàn)了過度降噪的情況，與原始仿真信號差別較大。

圖6可以直觀的看出CEEMDAN分解的A4算法降噪重構(gòu)的仿真信號比原始信號光滑且相似度高，降噪信號與未加噪信號x3相關(guān)度高達(dá)0.9969；并且重構(gòu)信號與降噪信號之差與原始(x1+x2)噪聲信號極其相似，其相關(guān)度為0.9134，降噪效果較好。

表1計算出三種方法分解下組成的最優(yōu)算法的含噪信號與降噪誤差比，CEEMDA分解的A4算法的值最小，該值越小，降噪誤差越大，降噪信號越偏離含噪信號，則降噪效果越好。

實際上，當(dāng)降噪均方誤差權(quán)重分別取0.1～0.3時，基于CEEMDAN分解的算法A5都最優(yōu)，其次是算法A4；當(dāng)降噪均方誤差權(quán)重分別取0.4～0.9時，算法A4最優(yōu)。圖7為降噪均方誤差權(quán)重因子分別取0.1～0.9時，CEEMDAN分解的14個降噪算法對應(yīng)的目標(biāo)函數(shù)，由圖可以看出：在14個算法中，算法A4的平均目標(biāo)函數(shù)值最小。

表1 含噪信號與降噪誤差比Table 1 Noise signal to noise reduction error ratio

3 塌落觸地振動信號時頻分析

為了驗證算法的實用性，以實測水塔爆破觸地振動信號為依據(jù)，對其進(jìn)行基于CEEMDAN的最優(yōu)降噪光滑模型算法降噪后進(jìn)行希爾伯特變換對比不同測點的頻率及能量分布。

3.1 工程概況

待拆除水塔高38.1 m，其中，塔身高28.6 m，倒錐傘形水箱高9.5 m。水塔下部支筒為圓筒形，筒體內(nèi)徑為2.0 m，壁厚0.18 m，筒內(nèi)設(shè)有5個設(shè)備檢修平臺，上部水箱為倒錐形結(jié)構(gòu)，傘形水箱最外徑為10.64 m。水塔筒體及水箱配有主筋和箍筋，為研究高聳構(gòu)筑物觸地振動規(guī)律，采用原地坍塌總體爆破方案，爆破切口布設(shè)在水塔底部以及各層檢修平臺處，共計6個爆破切口，其結(jié)構(gòu)示意圖及切口布置如圖8所示，起爆網(wǎng)路延期時間如表2所示。利用TC-4850N爆破測振儀實測其速度(儀器采樣頻率設(shè)置2000 Hz，采樣時長6 s)，振動測點布置如圖9所示，振動監(jiān)測數(shù)據(jù)如表3所示，其中三個測點信號三向合速度依次記為S1、S2、S3，降噪后的信號記為S1-1、S2-2、S3-3。

3.2 信號降噪

針對實測的信號S1、S2、S3進(jìn)行CEEMDAN分解后分別得到13個固有模態(tài)分量IMF及余項r(t)，同樣產(chǎn)生14個帶通濾波算法，考慮到觸地振動頻率一般在1～10 Hz，屬于低頻率[15]，濾波之后光滑度權(quán)重大于相似度，所以均方降噪誤差權(quán)重因子取值0.3或0.4較為合理。當(dāng)α=0.3時，14個帶通算法的降噪均方誤差、折線曲折度、目標(biāo)函數(shù)如圖10所示(限于文章篇幅只列取S1信號降噪效果圖)。

由圖10可看出，當(dāng)降噪均方誤差權(quán)重因子取0.3時，算法A4的目標(biāo)函數(shù)值最小，即算法A4的降噪效果最好；由圖11可以看出，經(jīng)算法A4降噪重構(gòu)后的信號與原始爆破振動信號相似度較高，在2.5～4 s信號幅值有所減小，這可能是水塔塌落過程中結(jié)構(gòu)碰撞或觸地引起的噪聲；圖12截取S1及S1-1信號0～1 s數(shù)據(jù)，可以看出，基于CEEMDAN分解的A4算法降噪重構(gòu)后信號比原始信號光滑且相似度高，與仿真信號降噪效果一致。

表2 起爆網(wǎng)路延期時間(單位：ms)Table 2 Detonation network delay time(unit:ms)

表3 水塔爆破觸地振動數(shù)據(jù)監(jiān)測結(jié)果Table 3 Water tower blasting ground vibration data monitoring results

為了進(jìn)一步說明實測信號在三種方法分解下的降噪效果，表4表計算出三個信號在三種方法分解下組成最優(yōu)光滑降噪算法的信噪比，相比EMD分解A2算法和EEMD分解A2算法，由CEEMDAN分解的A4算法的信噪比值最小，說明降噪效果最好，與仿真信號降噪效果一致。

3.3 信號邊際能量譜分析

對振動信號的邊際譜幅值的平方在全局時間上進(jìn)行積分得到振動信號的邊際能量譜，信號邊際能量譜上的曲線函數(shù)代表在不同頻率時信號的總能量幅值?，F(xiàn)對降噪后的S1-1、S2-2、S3-3水平信號進(jìn)行希爾伯特-黃變換，得到其邊際能量譜，如圖13所示。

表4 含噪信號與降噪誤差比Table 4 Noise signal to noise reduction error ratio

由圖13可以看出：信號S1-1、S2-2、S3-3邊際譜在頻率為5 Hz左右處能量出現(xiàn)最大值，說明其在頻率5Hz附近能量最集中；在能量達(dá)到最大值后緊接著又出現(xiàn)了兩個幅值較大的曲線，這主要是由于切口之間爆破時間間隔小，前后切口爆破倒塌觸地產(chǎn)生的沖擊波累加導(dǎo)致頻率有所增大。信號S1-1,S2-2,S3-3的邊際能量譜幅值逐漸的減小，也間接的說明振動波的能量在傳播過程中在不斷地衰減，并且隨著傳播距離的增加衰減速度逐漸變慢。

3.4 信號三維時頻能量譜分析

為直觀的看到降噪后信號的時頻、能量分布，截取前4 s的S1-1、S2-2、S3-3水平信號再次經(jīng)過CEEMDAN分解后，進(jìn)一步進(jìn)行希爾伯特黃變換得到三維時頻能量譜，如圖14所示。

由圖13可以看出，經(jīng)塌落觸地振動信號的主要頻率集中在頻率0～20 Hz的頻帶內(nèi)，大于20 Hz的頻帶內(nèi)攜帶的能量較少；偶爾在高頻帶段內(nèi)有時會突然出現(xiàn)局域能量峰值，通過峰值以后迅速衰減，這一現(xiàn)象表明塌落觸地振動產(chǎn)生的地震波在高頻部分衰減速度比在低頻部分的衰減速度要快；隨著測點距離增加，能量向低頻帶集中，并且3 s后各頻帶所攜帶能量明顯減少；塌落觸地振動波在傳播過程中不同頻帶能量的持續(xù)時間不同，低頻部分振動波持續(xù)時間更長，幾乎貫穿地震波傳播的整個過程，這表明振動波主要以低頻波攜帶的能量為主。

4 結(jié)論

實測的非平穩(wěn)爆破振動信號往往夾雜較大的噪聲，因此，信號降噪在原始信號特征提取時顯得尤為重要，本文將仿真信號分別進(jìn)行EMD、EEMD、CEEMDAN分解，通過設(shè)計合理的濾波器，得到不同條件下的濾波算法，根據(jù)降噪均方誤差指標(biāo)、曲線曲折度指標(biāo)建立最優(yōu)光滑降噪算法的判斷準(zhǔn)則，最終確定基于CEEMDAN分解得到的算法A4降噪效果最好，其信噪比值最小。

將實測塌落觸地振動信號進(jìn)行降噪處理，對比分析了三種分解方法下信號S1的信噪比，結(jié)果表明，基于CEEMDAN分解的最優(yōu)降噪算法信噪比值最小，驗證了該算法模型具有良好的降噪能力以及工程實用性，為實現(xiàn)精細(xì)爆破及提高振動危害評價的準(zhǔn)確度都具有一定的應(yīng)用價值。

對降噪后的觸地振動信號進(jìn)行時頻分析，結(jié)果表明，塌落觸地振動信號的主要頻率集中在頻率0～20 Hz的頻帶，且隨著距離的增加能量向低頻帶集中。