基于小波變換算法的普朗銅礦爆破振動信號能量分布研究

劉明武，黃永輝，李爭榮，賈皓琦，甘等俊，李洪超，張智宇

（1.云南迪慶有色金屬有限責(zé)任公司，云南 迪慶 674400；2.昆明理工大學(xué)電力工程學(xué)院，云南 昆明 650500；3.昆明理工大學(xué)國土資源工程學(xué)院，云南 昆明 650093；4.昆明理工大學(xué)城市學(xué)院，云南 昆明 650051）

0 引言

地下礦山爆破面臨十分復(fù)雜的環(huán)境[1-2]，井下爆破產(chǎn)生的地震波、沖擊波、飛石等對井下巷道、支護(hù)和工業(yè)設(shè)施等會造成嚴(yán)重的破壞。炸藥爆炸能量除一部分對爆區(qū)主體產(chǎn)生影響外，還有相當(dāng)一部分轉(zhuǎn)化為振動、噪聲和沖擊波。而在大部分的爆破振動監(jiān)測中，只將最大振動速度作為評判安全依據(jù)，忽略了爆破振動能量在頻域上分布是不均勻的。即使整體爆破振動速度沒有達(dá)到臨界值，但是在某頻段的能量集中也可能導(dǎo)致事故的發(fā)生[3]，因此，有必要對爆破振動的時域特征進(jìn)行進(jìn)一步的討論。龍翼等[4]利用決策樹分類算法能有效識別巖體破裂信號，消除噪音信號的影響。魯超等[5]利用小波與EMD 算法得出測點(diǎn)信號能量熵隨頻率降低呈先增后降的趨勢。汪平等[6]通過HHT 分析了小間距巷道爆破時，相鄰的巷道爆炸洞壁處振動能量大部分集中于40 Hz。何理等[7]對開挖巷道的上方邊坡進(jìn)行爆破振動監(jiān)測并進(jìn)行小波分析，發(fā)現(xiàn)隨著距離增加，能量衰減呈現(xiàn)先急速后平緩的趨勢。張聲輝等[8]通過EMD 算法對爆破振動信號主成分進(jìn)行了提取，發(fā)現(xiàn)EMD 算法分解出來的IMF 分量呈現(xiàn)出一定尺度范圍的模態(tài)。韋嘯等[9]通過EEMD 算法對信號進(jìn)行分析，得出隧道爆破能量一般集中在5～80 Hz，且主要分布在50 Hz 以下低頻段。部分學(xué)者[10-13]通過EMD 算法、EEMD 算法、HHT 算法等對爆破振動信號進(jìn)行能量研究，但是EMD 算法存在明顯的缺點(diǎn)，比如端點(diǎn)混疊效應(yīng)、模態(tài)混疊效應(yīng)等并不能很好地對爆破振動信號進(jìn)行研究分析，而且在大部分的爆破振動監(jiān)測中，只將振動幅值最大速度作為評判安全依據(jù)，忽視了爆破振動能量在頻域上的分布情況。

為解決上述EMD 算法等出現(xiàn)端點(diǎn)混疊效應(yīng)以及不能精確分析爆破振動能量在頻域上的分布情況等問題，可采用小波變換算法進(jìn)行分析。目前有學(xué)者[14-19]采用小波變換算法、改進(jìn)的小波算法以及分形理論等方法對爆破振動信號的頻率、能量分布特征進(jìn)行研究，利用小波分析可以將爆破振動信號按照不同的主振頻帶對信號進(jìn)行分解及重構(gòu)，從而能夠更加精準(zhǔn)地評價爆破振動危害。

本文以普朗銅礦為背景，通過數(shù)據(jù)處理軟件對自然崩落法鉆爆施工過程中采集到的爆破振動數(shù)據(jù)進(jìn)行分析，借助小波變換算法、HHT 算法、STFT 算法研究普朗銅礦爆破振動信號能量衰減規(guī)律，為普朗銅礦自然崩落法鉆爆施工過程中底部結(jié)構(gòu)在爆破振動荷載下的穩(wěn)定分析提供依據(jù)。

1 小波變換算法和HHT 分析法

1.1 小波變換算法

小波變換算法是將小波函數(shù)處于不同分解水平上，窗口尺寸隨分解水平變化的情況下，進(jìn)行信號時域局部分析，計算見式（1）。

式中：ψ（t）為一個基本小波函數(shù)或母小波函數(shù)；L2（R）為能量限量的信號空域；t為時間，s。

將 ψ（t）滑移和伸展得以取得1 組小波。對于接連不斷的情形，計算見式（3）。

式中：ɑ為尺寸因素；b為伸展因素。對于失散情形，計算見式（4）。

由式（3）和式（4）可知，在爆破振動領(lǐng)域處理振動信號使用最多的小波基為db8，因此，本文選用小波基為db8。

子頻帶能量分布：將被剖析信號拆分到第n層，假定與Sn，j呼應(yīng)的能效為En，j，計算見式（5）。

式中：Xj，k為重組信號Sn，j（t）的散落點(diǎn)的波譜幅值；j=0,1,2,…,2n-1；k=0,1,2,…,m；m為信號散落取樣基數(shù)。

設(shè)被分析信號的總能量為E0，計算見式（6）。

每個頻譜的能量占被解析信號總能量的比率，計算見式（7）。

由式（5）～式（7）可獲得信號經(jīng)過小波包分解之后不同子頻帶內(nèi)的能量情況，進(jìn)而可發(fā)現(xiàn)爆破振動信號傳播時能量變換的規(guī)律。

1.2 HHT 分析法

希爾伯特變換（HHT）是對非穩(wěn)態(tài)爆破振動信號進(jìn)行時頻分析的一種方法，突破了傅里葉變換分析方法（FT 變換）的局限性，計算見式（8）和式（9）。

式中：PV為柯西主分量；xi（t）為爆破振動的分量；ci（t）為振動信號的第一個分量。

將振幅H（ω，t）對時間T積分，得到Hilbert 的邊際譜h（ω），計算見式（10）。

式中：ω為頻率；H（ω，t）為爆破振動信號的振幅。

Hilbert 的瞬時能量IE（t）是將振幅H（ω，t）的平方對頻率積分得到，計算見式（11）。

2 爆破振動監(jiān)測及分析

2.1 工程概況

普朗銅礦所采用的自然崩落法在國內(nèi)及世界范圍內(nèi)均處于先進(jìn)水平。礦山位于香格里拉，其巖體以塊狀石英二長斑巖、石英閃長玢巖為主，巖石堅硬，中等穩(wěn)固，圍巖穩(wěn)定性較好。本次拉底爆破設(shè)計孔數(shù)為144 個；孔網(wǎng)參數(shù)最小抵抗線為2.0 m；最大孔底距為1.9～2.4 m；孔網(wǎng)密集系數(shù)為0.95～1.20；設(shè)計裝藥總長度為772.4 m；裝藥系數(shù)為70.2%；炸藥總用量為2 295 kg；每米深孔崩落礦量為5.45 t/m。

2.2 測點(diǎn)布置及監(jiān)測成果

選定南部區(qū)域S12N 巷道爆區(qū)為試驗(yàn)對象，爆破總藥量為2 295 kg，最大單響裝藥量為579.4 kg。依次把爆破測振儀放置在S13N 巷道、S13S 巷道、S14N 巷道、S14S 巷道，傳感器X軸方向?qū)?zhǔn)爆區(qū)，整體布置情況如圖1 所示，爆破振動測試儀離爆區(qū)最近距離約為50 m，離爆區(qū)最遠(yuǎn)距離約為200 m，具體爆破振動監(jiān)測成果見表1。

表1 爆破振動監(jiān)測成果Table 1 Results of blasting vibration monitoring

圖1 測點(diǎn)布置圖Fig.1 Layout of measuring points

2.3 振動波能量分布衰減情況

為了更好地分析出振動波能量分布衰減情況，對S13N 巷道、S13S 巷道、S14N 巷道、S14S 巷道進(jìn)行監(jiān)測，一號機(jī)、二號機(jī)爆心距分別為50 m、100 m，距離爆區(qū)較近，將其歸類為爆區(qū)近區(qū)；三號機(jī)、四號機(jī)爆心距分別為150 m、200 m，距離爆區(qū)較遠(yuǎn)，將其歸類為爆區(qū)遠(yuǎn)區(qū)，具體振動信號波形如圖2 所示。由圖2 可知，四臺爆破振動測試儀測出的振動波形時間均在2 000 ms 以內(nèi)，其波形緊密程度不一。通過小波變換算法對任意一個振動波形進(jìn)行分解，得到6 個模態(tài)分量IMF，具體情況如圖3 所示。

圖2 爆破振動波形Fig.2 Blasting vibration waveform

圖3 IMF 分量Fig.3 IMF component

由圖3 可知，模態(tài)分量IMF 在局部頻帶區(qū)間上呈由高到低的有序排列，同時表現(xiàn)出高頻、低幅值和波長短等特征，這些成分是原始信號在復(fù)雜環(huán)境監(jiān)測條件下受噪聲干擾的重要組成部分，充分證明小波變換算法能夠有效反映信號的局部性特征。

當(dāng)爆破振動信號頻率超出爆破振動測試儀的最小工作頻率后，則可能引起信號失真。因此，當(dāng)采用小波變換算法對爆破振動信號進(jìn)行分解時，其分解層數(shù)要根據(jù)爆破振動測試儀的工作頻帶來確定[20-21]。此次爆破振動測試儀工作頻率為16 kHz，為了進(jìn)一步開展爆破振動能量衰減的研究，運(yùn)用小波變換算法對爆破振動信號進(jìn)行九層的小波分解，選用db8 作為本次分析的小波基，根據(jù)Nyquist 定理，其頻率為8 kHz，小波為九層分解，則有512 個小波包，最低頻帶為0～15.62 Hz。根據(jù)小波變換算法，分別把4 個振動信號分解成8 個頻率段，每個頻率段的能量占比具體情況如圖4 所示。

圖4 0～125 Hz 能量占比Fig.4 Percentage of 0-125 Hz energy

由圖4 可知，振動波傳播到50 m、100 m 時，其能量分布情況為第2 頻率段能量占比上升達(dá)到最高而后減?。徽駝硬▊鞑サ?50 m、200 m 時，其能量分布情況為第4 頻率段能量占比上升達(dá)到最高而后減小。將每個頻率段能量占比數(shù)值化，見表2。

表2 能量分布情況Table 2 Distribution of energy

由表2 可知，一號機(jī)、二號機(jī)最高子頻帶能量占比位于第2 頻率段，能量占比分別為24.4%、25.3%，三號機(jī)、四號機(jī)最高子頻帶能量占比位于第4 頻率段，能量占比分別為22.9%、21.9%。在15.62～＜31.2 Hz頻率段，一號機(jī)、二號機(jī)能量占比24.4%、25.3%，三號機(jī)、四號機(jī)能量占比10.8%、10.1%，爆區(qū)近區(qū)比爆區(qū)遠(yuǎn)區(qū)能量占比平均多出14.1%；在46.87～＜93.75 Hz頻率段，三號機(jī)、四號機(jī)能量占比61.7%、58.7%，一號機(jī)、二號機(jī)能量占比37.8%、45.3%，爆區(qū)遠(yuǎn)區(qū)比爆區(qū)近區(qū)能量占比平均多出18.65%。結(jié)果表明，爆區(qū)遠(yuǎn)區(qū)較爆區(qū)近區(qū)，在46.87～＜93.75 Hz 頻率段能量較高，能量分布相對集中，第4 頻率段、第5 頻率段、第6 頻率段能量占比均高于爆區(qū)近區(qū)。

通過把4 個爆破振動波形的每條子頻帶能量占比畫成折線圖，可以看出不同爆心距的爆破振動波形能量運(yùn)動趨勢，具體如圖5 所示。

圖5 能量運(yùn)動趨勢Fig.5 Trend of energy movement

由圖5 可知，一號機(jī)、二號機(jī)能量運(yùn)動趨勢為在第2 頻率段能量占比達(dá)到最大，在后續(xù)的頻率段中能量占比迅速減小；三號機(jī)、四號機(jī)能量運(yùn)動趨勢為在第1 頻率段、第2 頻率段、第3 頻率段能量占比緩慢增加，直到第4 頻率段能量占比達(dá)到最大，在第4頻率段、第5 頻率段、第6 頻率段能量占比相較一號機(jī)、二號機(jī)減少速度較為緩慢，三號機(jī)、四號機(jī)較一號機(jī)、二號機(jī)能量分布更為寬廣。

結(jié)果表明，爆破振動波在傳播距離較近時，能量會在短時間內(nèi)迅速上升而后會隨著時間增加能量占比迅速減小，且能量分布不密集；振動波在傳播距離較遠(yuǎn)時，能量上升速度緩慢并且能量分布相對較集中。

2.4 振動波HHT 頻譜分析

HHT 算法在信號處理和故障診斷中被廣泛應(yīng)用，該方法由小波變換算法和希爾伯特-黃變換兩部分組成。小波變換算法可以將信號分解成多個子帶，從而獲取模態(tài)分量IMF。IMF 是原始信號的本征模態(tài)，可以反映出信號的局部特征。希爾伯特-黃變換可以將IMF 轉(zhuǎn)換成瞬時頻率和能量值，進(jìn)一步分析其時變性和頻域特征，具體效果如圖6 所示。

圖6 邊際譜能量分布Fig.6 Energy distribution of marginal spectrum

由圖6 可知，爆破振動能量無論傳播距離的遠(yuǎn)近，其能量主要分布在0～200 Hz 之間。由圖6（a）和圖6（b）可知，離爆區(qū)較近的區(qū)域，最高能量幅值均超過2×10-4dB，集中在0～50 Hz；由圖6（c）和圖6（d）可知，離爆區(qū)較遠(yuǎn)的區(qū)域，其能量幅值均在1.8×10-4dB，遠(yuǎn)遠(yuǎn)不如爆區(qū)近區(qū)的能量幅值，最高能量幅值分布在50～100 Hz 之間。通過HHT 算法分別將爆區(qū)近區(qū)、爆區(qū)遠(yuǎn)區(qū)能量三維化，計算得出能量三維圖，具體效果如圖7 所示。

圖7 能量三維圖Fig.7 Three-dimensional diagram of energy

由圖7 可知，爆區(qū)近區(qū)、爆區(qū)遠(yuǎn)區(qū)能量三維圖運(yùn)動趨勢對應(yīng)了圖6 邊際譜能量分布規(guī)律，爆區(qū)近區(qū)大部分能量分布在0～50 Hz，爆區(qū)遠(yuǎn)區(qū)大部分能量分布在40～80 Hz。爆區(qū)遠(yuǎn)區(qū)0～50 Hz 之間的部分能量小于爆區(qū)近區(qū)的能量，爆區(qū)遠(yuǎn)區(qū)50～80 Hz 之間部分能量大于爆區(qū)近區(qū)能量。爆區(qū)近區(qū)能量主要分布在0～50 Hz，且整體能量分布趨勢相對較窄，爆區(qū)遠(yuǎn)區(qū)大部分能量分布在40～80 Hz，且整體能量分布趨勢相對較寬。爆區(qū)近區(qū)能量在達(dá)到最高點(diǎn)以后迅速衰減，爆區(qū)遠(yuǎn)區(qū)能量在達(dá)到最高點(diǎn)以后衰減緩慢，且能量分布相對集中。

2.5 振動波短時傅里葉變換分析

短時傅里葉變換（STFT）是時頻分析方法之一，該方法把信號劃分為時間域內(nèi)的幾個短時段，并在每一個短時段內(nèi)對窗函數(shù)加窗，然后進(jìn)行傅里葉變換以獲得各個時刻的頻率分量。STFT 算法比離散傅里葉變換及連續(xù)傅里葉變換有更高的時間分辨率及頻率分辨率。為進(jìn)一步研究爆破振動波能量衰減規(guī)律，采用STFT 算法對爆區(qū)近區(qū)、爆區(qū)遠(yuǎn)區(qū)振動波形所含能量進(jìn)行STFT 變換運(yùn)算，如圖8所示。

由圖8 可知，圖中Z軸代表能量功率譜幅值，爆區(qū)近區(qū)能量高功率幅值均在8.59～14.23 dB/Hz，時間分布在0.13～0.23 s，爆區(qū)遠(yuǎn)區(qū)能量高功率幅值均在8.35～11.98 dB/Hz，時間分布在0.15～0.35 s。爆區(qū)近區(qū)能量高功率幅值普遍大于爆區(qū)遠(yuǎn)區(qū)能量高功率幅值，是由爆破振動傳播距離不同而產(chǎn)生的差異。爆區(qū)近區(qū)、爆區(qū)遠(yuǎn)區(qū)功率譜密度時頻三維圖后續(xù)走向也存在差異，隨著爆破振動傳播的增加，爆區(qū)遠(yuǎn)區(qū)的能量在達(dá)到最高點(diǎn)以后較爆區(qū)近區(qū)能量衰減緩慢，并且時間越長，后續(xù)能量差越大，進(jìn)一步證明爆區(qū)遠(yuǎn)區(qū)能量衰減較為緩慢，能量分布相對集中。

為了使振動波形中的能量功率趨勢表現(xiàn)得更清晰，使用數(shù)據(jù)處理軟件繪制出功率趨勢圖，具體效果如圖9 所示。

圖9 功率趨勢圖Fig.9 Trend chart of power

由圖9 可知，爆區(qū)近區(qū)著時間增加功率迅速減小，且能量分布不集中，而爆區(qū)遠(yuǎn)區(qū)功率隨時間增加，功率較爆區(qū)近區(qū)衰減速度更為緩慢。

3 結(jié)論

1）普朗銅礦爆區(qū)近區(qū)主振頻率主要分布在15.62～＜31.20 Hz，爆區(qū)遠(yuǎn)區(qū)主振頻率主要分布在46.87～＜93.75 Hz；15.62～＜31.2 Hz 頻率段爆區(qū)近區(qū)比爆區(qū)遠(yuǎn)區(qū)能量占比平均高出14.1%，46.87～＜93.75 Hz 頻率段爆區(qū)遠(yuǎn)區(qū)要比爆區(qū)近區(qū)能量占比平均高出18.65%。

2）爆破振動波在傳播距離較近時，能量會在短時間內(nèi)迅速上升，而后隨主振頻率增加能量占比迅速減小，且能量分布不集中；振動波在傳播距離較遠(yuǎn)時，能量上升速度緩慢且能量達(dá)到最高后其后續(xù)衰減速度較為緩慢，且能量分布更為集中。

3）實(shí)測振動信號能量特征表明，基于小波包算法能量分析方法可與爆破振動非平穩(wěn)信號匹配，不僅準(zhǔn)確了解爆破振動能量-頻譜分布，還給出了不同頻帶上的振動能量分布和衰減規(guī)律，為指導(dǎo)井下構(gòu)筑物結(jié)構(gòu)抗震設(shè)計和工程爆破監(jiān)測提供分析依據(jù)。