基于改進MP-WVD算法的核電廠建設爆破振動信號處理方法*

梅 比，汪旭光，楊仁樹

（1.中國礦業大學（北京）力學與建筑工程學院，北京 100083；2.北京礦冶研究總院，北京 100160）

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

1 算法基本原理概述

1.1 WVD分布

WVD算法最初是由Ville將其應用到信號的時-頻分析領域的[7]，針對一個時間序列W(t,ω)根據特征函數方法推導出WVD表達式：

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

1.2 MP算法

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

MP算法是一種貪婪算法，該方法與統計學中使用的投影追蹤算法和波形增益矢量量化有密切關系。信號的分解需要在超完備子波庫中進行，其中，超完備意為信號的分解目標在信號所組成的空間中足夠密集，這也就是最終無法以一組正交基進行描述。對于任意一個有限維 Hilbert 空間H，D為此空間內的一個超完備詞典，設信號為f∈H，長度為N，D中的元素滿足：

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

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

MP算法就是將信號f垂直投影到子波庫D的子波上。設則f可以表示為：

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

設，經過了n次迭代(n>>0)得到殘余信號=f此時再選擇一個匹配子波使其匹配即：

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

2 MP-WVD算法的缺陷及改進

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

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

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

由式（7）可以看出，交叉項的存在對信號的時-頻分布產生了很大的干擾。因此，將MP算法與WVD結合起來，將原始信號細分成基本原子，然后再對每個信號做WVD，將得到的結果疊加，在消除干擾項的基礎上得到更清晰的時頻分布圖譜。具體示例如下。

示例信號x(t)，如圖2所示。對示例信號x直接進行WVD運算，得到的時頻分布情況如圖3所示。可以看出圖中不僅有4個真實分量的時頻能量分布，還在任意2個分量的時頻中心連線的中點處出現了交叉項，嚴重干擾了對信號真實時頻分布的判讀。

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

圖3 WVD時頻分布圖Fig.3 WVD time-frequency distribution

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

圖4 MP算法分解結果Fig.4 Decomposition results of MP algorithm

再分別對4個子波進行WVD運算，得到其各自的WVD時頻分布結果，如圖5所示；將各子波的時頻分布合并，即得到了原示例信號的時頻分布情況，如圖6所示。

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

另一方面，MP算法也存在著缺陷。其計算量過大，需耗費很長機時來完成一次信號分解[11]。而通常時頻分析數據量都較大，因此研究者希望通過改進得到一種更高效的算法。

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

與傳統算法進行效率對比，分別用常規MP算法與改進后得MP算法對合成信號進行處理：

由表1中對比分析可以看出，經過改進后的算法效率有了明顯提升，對同一信號分解所用機時有著50倍左右的差距。

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

圖7 程序框圖Fig.7 Program diagram

表1 兩種方法所用機時對比Table 1 Comparison of machine time between two methods

3 基于改進MP-WVD算法的核電爆破振動信號分析

3.1 工程概況

結合福建漳州核電廠一期工程場地平整土石方工程施工實例中監測到的爆破振動信號進行分析。本工程施工區域以5#、6#機組及其廠區西側為界，功能設施為5#機組與6#機組主生產區和廠前區[13-14]。開挖后場平標高，13.5 m；土方量，603 200m3；石方量，5 878 300m3；回填方量，3 573 300 m3；挖溝槽土方，648 m3；石方，1 513 m3；邊坡預裂面積，39 600 m2。

爆破振動監測點設置3個測點：人家村村委會東側民房設置測點1，距離爆破施工位置859 m；南山村地標設置測點2，距離爆破施工位置678 m；水壩閘門設置測點3，距離爆破施工位置291 m。本文中選取三標段爆破實例，深孔爆破，采用乳化90炸藥，其爆破技術參數如表2所示。

表2 爆破參數表Table 2 Blasting parameters

3.2 數據處理

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

圖8 原始信號Fig.8 Original signal

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

數據處理過程中發現，MP算法用于爆破振動信號高頻去燥取得了十分理想的效果，與其他濾波算法進行對比研究，結果如圖10所示。

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

將所得原子矩陣代入WVD算法，分別計算每個原子的時頻分布，再逐一疊加，得到信號的總體時頻分布如圖11所示，基于小波變換的時頻分析結果如圖12所示。

圖10 濾波效果對比Fig.10 Comparison of filtering effects

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

將該測點的水平徑向分量（X軸）以及水平切向分量（Y軸）信號分別代入MP-WVD算法，得到其時頻分布情況，分別如圖13、14所示。

從圖13、14中可以看出，水平切向振動信號的頻率中心為21 Hz，水平徑向振動信號在21 Hz和13 Hz處都有能量聚集。3個方向的振動信號均具有初始頻帶寬、高頻衰減快的特點。不同點在于：水平切向的信號沒有出現持續時間較長的13 Hz分量，且水平切向的能量分布明顯比其他2個方向的更集中。

圖12 小波變換結果Fig.12 Result of wavelet transform

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

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

對信號的時頻分布函數分別進行時間積分和頻率積分，得到信號的瞬時能量譜如圖15所示，頻率邊際譜如圖16所示。

圖15 瞬時能量譜Fig.15 Instantaneous energy spectrum

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

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

不同方向的傳感器收集的地震波類型不同，其中瑞利波的質點運動方向對應X軸和Z軸，勒夫波的質點運動方向對應Y軸和X軸。邊際譜縱坐標E·f-1表示能量在頻率軸上的分布密度。由瞬時能量譜可以看出，水平切向（Y軸）的振動峰值率先到達，鉛垂方向（Z軸）的振動峰值有較大的時延，而水平徑向（X軸）的振動有多個峰值。4個方向的瞬時能量曲線反映出了在此次爆破引起的振動中，勒夫波率先到達測點，瑞利波相較勒夫波有約0.05 s的時延。

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

將南山村測點以及人家村測點所采集到的數據代入算法進行處理，得到的結果如圖17、18所示。

圖17 南山村測點數據分析結果Fig.17 Data analysis result of Nanshan village

圖18 人家村測點數據分析結果Fig.18 Data analysis result of Renjia village

從圖17、18中可以看出，由于測試點距離起爆點較遠，瑞利波在傳遞過程中的色散效應產生了到達時間差，導致振動能量在時域上的分布更分散。

在頻域上，能量分布更加向低頻集中，絕大部分能量分布在10～20 Hz頻帶，20 Hz以上頻帶幾乎沒有分布。可以看出在爆破地震波中，高頻成分衰減較快，而低頻成分傳播較遠，爆破遠區的抗震設計當以應對10～20 Hz頻段的低頻振動為主。

4 結 論

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