基于改進(jìn)MP-WVD算法的核電廠建設(shè)爆破振動(dòng)信號(hào)處理方法*

梅 比，汪旭光，楊仁樹

（1.中國(guó)礦業(yè)大學(xué)（北京）力學(xué)與建筑工程學(xué)院，北京 100083；2.北京礦冶研究總院，北京 100160）

我國(guó)核電建設(shè)進(jìn)入加速期，由于核電站規(guī)模較大，無法多機(jī)組同時(shí)建設(shè)，因此在后期機(jī)組施工時(shí)就要考慮到已建機(jī)組的安全穩(wěn)定，尤其是爆破施工產(chǎn)生的振動(dòng)影響[1]。研究表明，建（構(gòu)）筑物在地震波作用下產(chǎn)生的動(dòng)力響應(yīng)特征與波的時(shí)頻能量特性密切相關(guān)。目前時(shí)-頻分析算法有各自優(yōu)缺點(diǎn)：STFT算法(short-time Fourier transform)最初為傅里葉變換應(yīng)用于時(shí)頻分析提供了可能，但是郭濤等[2]通過對(duì)傳統(tǒng)的STFT算法與小波算法進(jìn)行對(duì)比分析得到了傳統(tǒng)的STFT以三角函數(shù)作為基函數(shù)，無法準(zhǔn)確與爆破振動(dòng)這種非穩(wěn)態(tài)隨機(jī)過程相匹配的結(jié)論。馬華原等[3]將小波變化方法運(yùn)用于核電施工爆破振動(dòng)分析取得了良好的效果，但也發(fā)現(xiàn)小波變換雖然引入了更加適配的小波基，但時(shí)頻分辨率受測(cè)不準(zhǔn)原理限制。李夕兵等[4]、魏新江等[5]改進(jìn)了HHT算法(Hilbert-Huang transform)并將其運(yùn)用于爆炸振動(dòng)分析中，獲得了具有比小波變換分辨率更高的時(shí)頻分布圖譜。但隨著研究的深入，HHT變換提出的自適應(yīng)基底導(dǎo)致算法邊界效應(yīng)較大，這一點(diǎn)目前仍是個(gè)難題。另一方面，周輝等[6]對(duì)傳統(tǒng)的匹配追蹤算法進(jìn)行了改進(jìn)，提供了消除經(jīng)典WVD算法(Wigner-Ville distribution)交叉干擾項(xiàng)的思路，在此基礎(chǔ)上，本文將改進(jìn)型MP算法(matching pursuit algorithm)與傳統(tǒng)WVD算法相結(jié)合，成功解決了交叉項(xiàng)干擾的問題，同時(shí)很好地控制了算法的計(jì)算復(fù)雜度。結(jié)合核電爆破信號(hào)實(shí)例，取得了較好的分析效果。

1 算法基本原理概述

1.1 WVD分布

WVD算法最初是由Ville將其應(yīng)用到信號(hào)的時(shí)-頻分析領(lǐng)域的[7]，針對(duì)一個(gè)時(shí)間序列W(t,ω)根據(jù)特征函數(shù)方法推導(dǎo)出WVD表達(dá)式：

式中，t為時(shí)間，ω為角頻率。由式(1)可知，在信號(hào)WVD表達(dá)中，不存在任何形式的窗函數(shù)。因此，WVD分相對(duì)STFT、小波變換而言有著更高的時(shí)間和頻率分辨率靈活性。

1.2 MP算法

MP算法是由Mallat于1993年提出的一種信號(hào)分解算法[8]。其核心原理是將信號(hào)以字典原子為基，進(jìn)行分解。過程如圖1所示：選取與信號(hào)Xn匹配程度最高的原子Ψn，并求出投影值an和差值信號(hào)Xn+1。此時(shí)得到殘差信號(hào)Xn+1，重復(fù)進(jìn)行原子匹配過程，將其投影到與其最相近的原子Ψn+1上，得到差值信號(hào)Xn+2。以此類推，直到殘余信號(hào)的能量小于設(shè)定的閾值為止。

MP算法是一種貪婪算法，該方法與統(tǒng)計(jì)學(xué)中使用的投影追蹤算法和波形增益矢量量化有密切關(guān)系。信號(hào)的分解需要在超完備子波庫(kù)中進(jìn)行，其中，超完備意為信號(hào)的分解目標(biāo)在信號(hào)所組成的空間中足夠密集，這也就是最終無法以一組正交基進(jìn)行描述。對(duì)于任意一個(gè)有限維 Hilbert 空間H，D為此空間內(nèi)的一個(gè)超完備詞典，設(shè)信號(hào)為f∈H，長(zhǎng)度為N，D中的元素滿足：

圖1 分解過程示意圖Fig.1 Schematic diagram of decomposition process

式中：gγ為子波分解算子，Γ為伽馬函數(shù)。

MP算法就是將信號(hào)f垂直投影到子波庫(kù)D的子波上。設(shè)則f可以表示為：

式中：Rf的意義為原始信號(hào)f通過子波匹配進(jìn)行一次分解的殘余[9]。為了達(dá)到最好的分解效果，就必須使殘余信號(hào)盡可能小，因此就必然要使內(nèi)積項(xiàng)取最大。很顯然與Rf是正交的，因此：

設(shè)，經(jīng)過了n次迭代(n>>0)得到殘余信號(hào)=f此時(shí)再選擇一個(gè)匹配子波使其匹配即：

就是進(jìn)行了n+1次迭代得到的差值。因此MP算法即是利用式（5）描述的一個(gè)重復(fù)迭代過程，若迭代m次，則f可以表示為：

2 MP-WVD算法的缺陷及改進(jìn)

WVD分布為非線性時(shí)頻分布，也就是說其不具備線性分布的可加性，即兩信號(hào)和的 WVD并不等于每一個(gè)信號(hào)的 WVD之和[10]。

令x(t)=x1(t)+x2(t)，則：

式中：2Re[Wx1+x2(t,ω)]是x1(t)和x2(t)的交叉項(xiàng)。t為時(shí)間，ω為角頻率。

由式（7）可以看出，交叉項(xiàng)的存在對(duì)信號(hào)的時(shí)-頻分布產(chǎn)生了很大的干擾。因此，將MP算法與WVD結(jié)合起來，將原始信號(hào)細(xì)分成基本原子，然后再對(duì)每個(gè)信號(hào)做WVD，將得到的結(jié)果疊加，在消除干擾項(xiàng)的基礎(chǔ)上得到更清晰的時(shí)頻分布圖譜。具體示例如下。

示例信號(hào)x(t)，如圖2所示。對(duì)示例信號(hào)x直接進(jìn)行WVD運(yùn)算，得到的時(shí)頻分布情況如圖3所示?？梢钥闯鰣D中不僅有4個(gè)真實(shí)分量的時(shí)頻能量分布，還在任意2個(gè)分量的時(shí)頻中心連線的中點(diǎn)處出現(xiàn)了交叉項(xiàng)，嚴(yán)重干擾了對(duì)信號(hào)真實(shí)時(shí)頻分布的判讀。

圖2 示例信號(hào)x(t)Fig.2 Example signal x(t)

圖3 WVD時(shí)頻分布圖Fig.3 WVD time-frequency distribution

現(xiàn)將示例信號(hào)x(t)代入MP算法，分解為x1、x2、x3、x4共 4個(gè)子項(xiàng)，如圖4所示。

圖4 MP算法分解結(jié)果Fig.4 Decomposition results of MP algorithm

再分別對(duì)4個(gè)子波進(jìn)行WVD運(yùn)算，得到其各自的WVD時(shí)頻分布結(jié)果，如圖5所示；將各子波的時(shí)頻分布合并，即得到了原示例信號(hào)的時(shí)頻分布情況，如圖6所示。

可見，WVD變換在任意兩個(gè)有效信號(hào)之間產(chǎn)生一個(gè)交叉項(xiàng)。圖3中信號(hào)有4個(gè)分量共產(chǎn)生6個(gè)交叉項(xiàng)（其中交叉項(xiàng)5和6重疊），對(duì)N個(gè)分量則會(huì)產(chǎn)生N×(N-1)/2個(gè)交叉項(xiàng) 。而圖4、5、6中采用MP子波分解方法，成功剔除了交叉項(xiàng)干擾。

另一方面，MP算法也存在著缺陷。其計(jì)算量過大，需耗費(fèi)很長(zhǎng)機(jī)時(shí)來完成一次信號(hào)分解[11]。而通常時(shí)頻分析數(shù)據(jù)量都較大，因此研究者希望通過改進(jìn)得到一種更高效的算法。

算法選用Gabor子項(xiàng)的控制參量有：振幅、頻率、中心時(shí)間和相位[12]。對(duì)過完備子波庫(kù)的掃描過程利用窮舉法對(duì)這上述參量進(jìn)行優(yōu)化選擇。引入HHT算法，通過EMD(empirical mode decomposition)分解得到信號(hào)的固有模態(tài)IMF(intrinsic mode function)分量。對(duì)IMF分量進(jìn)行Hilbert變換，得到信號(hào)的瞬時(shí)優(yōu)勢(shì)頻率以及相位，并將其代入MP算法中，即可將4參數(shù)掃描運(yùn)算降維至2參數(shù)掃描，大大減少了程序循環(huán)步數(shù)。程序流程如圖7所示。

與傳統(tǒng)算法進(jìn)行效率對(duì)比，分別用常規(guī)MP算法與改進(jìn)后得MP算法對(duì)合成信號(hào)進(jìn)行處理：

由表1中對(duì)比分析可以看出，經(jīng)過改進(jìn)后的算法效率有了明顯提升，對(duì)同一信號(hào)分解所用機(jī)時(shí)有著50倍左右的差距。

圖5 各個(gè)子波的WVD分布Fig.5 WVD distribution of each wavelet

圖7 程序框圖Fig.7 Program diagram

表1 兩種方法所用機(jī)時(shí)對(duì)比Table 1 Comparison of machine time between two methods

3 基于改進(jìn)MP-WVD算法的核電爆破振動(dòng)信號(hào)分析

3.1 工程概況

結(jié)合福建漳州核電廠一期工程場(chǎng)地平整土石方工程施工實(shí)例中監(jiān)測(cè)到的爆破振動(dòng)信號(hào)進(jìn)行分析。本工程施工區(qū)域以5#、6#機(jī)組及其廠區(qū)西側(cè)為界，功能設(shè)施為5#機(jī)組與6#機(jī)組主生產(chǎn)區(qū)和廠前區(qū)[13-14]。開挖后場(chǎng)平標(biāo)高，13.5 m；土方量，603 200m3；石方量，5 878 300m3；回填方量，3 573 300 m3；挖溝槽土方，648 m3；石方，1 513 m3；邊坡預(yù)裂面積，39 600 m2。

爆破振動(dòng)監(jiān)測(cè)點(diǎn)設(shè)置3個(gè)測(cè)點(diǎn)：人家村村委會(huì)東側(cè)民房設(shè)置測(cè)點(diǎn)1，距離爆破施工位置859 m；南山村地標(biāo)設(shè)置測(cè)點(diǎn)2，距離爆破施工位置678 m；水壩閘門設(shè)置測(cè)點(diǎn)3，距離爆破施工位置291 m。本文中選取三標(biāo)段爆破實(shí)例，深孔爆破，采用乳化90炸藥，其爆破技術(shù)參數(shù)如表2所示。

表2 爆破參數(shù)表Table 2 Blasting parameters

3.2 數(shù)據(jù)處理

以水壩閘門測(cè)點(diǎn)的振動(dòng)信號(hào)的鉛垂分量（Z軸）數(shù)據(jù)為例，利用改進(jìn)MP-WVD算法進(jìn)行分析。本次爆破振動(dòng)測(cè)試采用了TC-4850型爆破測(cè)振儀，采樣頻率為10 000 Hz，采集到的典型爆破振動(dòng)信號(hào)如圖8所示。第一步做MP分解，所得到的子波集合如圖9所示。分解得到的子波均具有良好緊支性[15]，可良好地呈現(xiàn)出信號(hào)的細(xì)節(jié)。

圖8 原始信號(hào)Fig.8 Original signal

圖9 子波集合Fig.9 Collection of wavelets

數(shù)據(jù)處理過程中發(fā)現(xiàn)，MP算法用于爆破振動(dòng)信號(hào)高頻去燥取得了十分理想的效果，與其他濾波算法進(jìn)行對(duì)比研究，結(jié)果如圖10所示。

從圖10可看出，MP算法得到的重構(gòu)信號(hào)保留了幾乎所有的原始信號(hào)振動(dòng)細(xì)節(jié)，同時(shí)剔除了高頻噪聲。小波變換得到的重構(gòu)信號(hào)也能保持與原始信號(hào)的高度吻合，但有較大噪聲殘留，信噪比不如MP算法。EMD算法由于采用了自適應(yīng)基底[16]，其算法速度非常快，但自適應(yīng)基底導(dǎo)致濾波效果不穩(wěn)定，容易將有用信息一同剔除；其邊界效應(yīng)也會(huì)引入多余信號(hào)，重構(gòu)信號(hào)和原始信號(hào)差值較大。

將所得原子矩陣代入WVD算法，分別計(jì)算每個(gè)原子的時(shí)頻分布，再逐一疊加，得到信號(hào)的總體時(shí)頻分布如圖11所示，基于小波變換的時(shí)頻分析結(jié)果如圖12所示。

圖10 濾波效果對(duì)比Fig.10 Comparison of filtering effects

由圖12可知，小波變換結(jié)果不僅頻率分辨率較低，且時(shí)間分辨率也不如MP-WVD算法得到結(jié)果，且存在一定的邊界溢出問題。MP-WVD算法不僅有著較高的頻率分辨率，同時(shí)較好地反映了地震波的頻率分布隨時(shí)間的變化而發(fā)生改變。從圖11中可以清晰地分辨出此振動(dòng)信號(hào)的特征：振動(dòng)峰值在0.13 s左右到達(dá)，頻率中心為13 Hz。此時(shí)刻地震波頻率分布最寬，在10～40 Hz的范圍內(nèi)都有分布。最后高頻成分迅速衰減，只剩余13 Hz分量持續(xù)了較長(zhǎng)時(shí)間。

將該測(cè)點(diǎn)的水平徑向分量（X軸）以及水平切向分量（Y軸）信號(hào)分別代入MP-WVD算法，得到其時(shí)頻分布情況，分別如圖13、14所示。

從圖13、14中可以看出，水平切向振動(dòng)信號(hào)的頻率中心為21 Hz，水平徑向振動(dòng)信號(hào)在21 Hz和13 Hz處都有能量聚集。3個(gè)方向的振動(dòng)信號(hào)均具有初始頻帶寬、高頻衰減快的特點(diǎn)。不同點(diǎn)在于：水平切向的信號(hào)沒有出現(xiàn)持續(xù)時(shí)間較長(zhǎng)的13 Hz分量，且水平切向的能量分布明顯比其他2個(gè)方向的更集中。

圖12 小波變換結(jié)果Fig.12 Result of wavelet transform

圖13 水平徑向分量（X軸）Fig.13 Horizontal radial component

圖14 水平切向分量（Y軸）Fig.14 Horizontal tangential component

對(duì)信號(hào)的時(shí)頻分布函數(shù)分別進(jìn)行時(shí)間積分和頻率積分，得到信號(hào)的瞬時(shí)能量譜如圖15所示，頻率邊際譜如圖16所示。

圖15 瞬時(shí)能量譜Fig.15 Instantaneous energy spectrum

圖16 頻率邊際譜Fig.16 Marginal spectrum

頻率邊際譜中縱軸E·f1表示能量在頻率軸的分布密度, 量綱為J/Hz。

不同方向的傳感器收集的地震波類型不同，其中瑞利波的質(zhì)點(diǎn)運(yùn)動(dòng)方向?qū)?yīng)X軸和Z軸，勒夫波的質(zhì)點(diǎn)運(yùn)動(dòng)方向?qū)?yīng)Y軸和X軸。邊際譜縱坐標(biāo)E·f-1表示能量在頻率軸上的分布密度。由瞬時(shí)能量譜可以看出，水平切向（Y軸）的振動(dòng)峰值率先到達(dá)，鉛垂方向（Z軸）的振動(dòng)峰值有較大的時(shí)延，而水平徑向（X軸）的振動(dòng)有多個(gè)峰值。4個(gè)方向的瞬時(shí)能量曲線反映出了在此次爆破引起的振動(dòng)中，勒夫波率先到達(dá)測(cè)點(diǎn)，瑞利波相較勒夫波有約0.05 s的時(shí)延。

頻率邊際譜反映出信號(hào)的各個(gè)頻率成分在時(shí)間全局上的累加，與傅氏幅頻譜意義不同的是，傅氏譜只能反映信號(hào)的頻率成分存在的概率而不能描述每個(gè)分量在整個(gè)振動(dòng)過程中的能量份額。由圖16可以看出，水平切向和鉛垂方向振動(dòng)頻率較為集中，其頻率中心分別為13 Hz和21 Hz。水平徑向的振動(dòng)由于混雜了勒夫波和瑞利波兩個(gè)成分，所以具有多個(gè)頻率中心。尤其需注意的是在36 Hz左右，X軸信號(hào)的邊際譜有較強(qiáng)峰值出現(xiàn)，但在時(shí)頻譜上并未見較強(qiáng)的時(shí)頻能量集中，表明能量在此頻率處集中但在時(shí)間上并不集中，邊際譜的峰值是由于全局時(shí)間累積產(chǎn)生的，這種峰值小但加載時(shí)間長(zhǎng)的特殊能量加載形式也需高度重視，防止產(chǎn)生損傷累積效應(yīng)。

將南山村測(cè)點(diǎn)以及人家村測(cè)點(diǎn)所采集到的數(shù)據(jù)代入算法進(jìn)行處理，得到的結(jié)果如圖17、18所示。

圖17 南山村測(cè)點(diǎn)數(shù)據(jù)分析結(jié)果Fig.17 Data analysis result of Nanshan village

圖18 人家村測(cè)點(diǎn)數(shù)據(jù)分析結(jié)果Fig.18 Data analysis result of Renjia village

從圖17、18中可以看出，由于測(cè)試點(diǎn)距離起爆點(diǎn)較遠(yuǎn)，瑞利波在傳遞過程中的色散效應(yīng)產(chǎn)生了到達(dá)時(shí)間差，導(dǎo)致振動(dòng)能量在時(shí)域上的分布更分散。

在頻域上，能量分布更加向低頻集中，絕大部分能量分布在10～20 Hz頻帶，20 Hz以上頻帶幾乎沒有分布?？梢钥闯鲈诒频卣鸩ㄖ校哳l成分衰減較快，而低頻成分傳播較遠(yuǎn)，爆破遠(yuǎn)區(qū)的抗震設(shè)計(jì)當(dāng)以應(yīng)對(duì)10～20 Hz頻段的低頻振動(dòng)為主。

4 結(jié) 論

（1）將MP算法與WVD分布結(jié)合，有效地消除了交叉項(xiàng)的干擾，進(jìn)而發(fā)揮出WVD分布對(duì)瞬態(tài)信號(hào)敏感以及高分辨率的特點(diǎn)。引入HHT算法對(duì)信號(hào)的瞬時(shí)頻率以及瞬時(shí)相位先行確定，可有效降低計(jì)算復(fù)雜度，機(jī)時(shí)平均可縮短至之前的約2%。（2）改進(jìn)型MP算法用于處理地震波數(shù)據(jù)有著很好的去燥效果，能在盡量保留信號(hào)局部特征的基礎(chǔ)上剔除噪聲信號(hào)，相比小波濾波以及EMD濾波性能較為突出。（3）MP-WVD算法對(duì)于核電爆破施工振動(dòng)信號(hào)的分析結(jié)果較之小波變換方法具有更高的分辨率以及細(xì)節(jié)刻畫能力。（4）由WVD時(shí)頻譜衍生出的瞬時(shí)能量譜能夠清晰地反映出信號(hào)的能量峰值到達(dá)時(shí)間以及加載次數(shù)。頻率邊際譜能夠描述任意頻率成分在時(shí)間全局上的累積，這些特性在研究結(jié)構(gòu)動(dòng)態(tài)響應(yīng)過程中均可納入?yún)⒖肌?/p>