吳彥博，張 達(dá)，冀 虎，戴 銳，石雅倩

(1.礦冶科技集團有限公司，北京 102628；2.金屬礦山智能開采技術(shù)北京市重點實驗室，北京 102628；3.國家金屬礦綠色開采國際聯(lián)合研究中心，北京 102628)

我國經(jīng)濟快速發(fā)展導(dǎo)致金屬、非金屬等資源需求量日益增大，礦石開采量不斷增加，逐步形成了大量的采空區(qū)。部分采空區(qū)能夠依靠有效的充填、注漿等方法進(jìn)行及時治理，然而仍然有大批采空區(qū)不能夠被及時處理，地下采空區(qū)面積越來越大，從而導(dǎo)致地面塌陷、井下發(fā)生巖爆、冒頂、塌方等地壓災(zāi)害，不僅對井下安全開采造成威脅，也對周邊生態(tài)環(huán)境產(chǎn)生不可估量的破壞。微震監(jiān)測技術(shù)能夠利用井下巖體在擠壓變形、破壞以及微裂隙產(chǎn)生過程中所發(fā)出的震動波來進(jìn)行巖體分析和穩(wěn)定性評估，是關(guān)鍵的地壓監(jiān)測技術(shù)之一。但在微震監(jiān)測系統(tǒng)實際監(jiān)測過程中，系統(tǒng)采集到的數(shù)據(jù)不僅包含巖體破裂的微震信號，還包含開采爆破、機械作業(yè)、人員活動、電干擾、鑿巖信號等復(fù)雜震動信號，而且在傳播過程中這些震動信號會在時域和頻域上產(chǎn)生疊加，導(dǎo)致真正用于巖體穩(wěn)定性評估的微震信號難以分辨，微震分析中的數(shù)據(jù)不準(zhǔn)確，分析結(jié)果亦不可靠。

為了解決上述問題，首先要實現(xiàn)對震動信號開展時域及頻域的分析，只有分離出有效震動信號才能夠開展有效的波形識別，從而辨別出微震信號?，F(xiàn)已有大量文獻(xiàn)采用傅里葉頻譜分析或小波分析等手段對礦山微震信號進(jìn)行分析與識別，但仍存在通用性不足、計算方法復(fù)雜、計算量大等問題。本文研究了一種基于譜峰識別和頻域濾波技術(shù)的礦山多源震動信號分離與提取有效方法，能夠降低噪聲干擾對礦山微震信號的影響，為進(jìn)一步實現(xiàn)礦震信號的處理及分析提供新的思路。

1 基礎(chǔ)理論

1.1 傅立葉變換及其反變換

傅里葉變換是目前信號頻域分析最常用的方法之一，其能夠?qū)⑿盘枏臅r域轉(zhuǎn)換到頻域，進(jìn)而從頻域上得到很多在時域上不易發(fā)現(xiàn)的信號特征?？焖俑盗⑷~變換(fast fourier transform，F(xiàn)FT)是傅立葉變換中一種的高效算法，能夠?qū)r域信號轉(zhuǎn)到復(fù)數(shù)域信號來快速完成頻域轉(zhuǎn)換。其中，復(fù)數(shù)域信號中的實數(shù)序列可以看做虛部為零的復(fù)數(shù)，例如某實數(shù)信號x(n)可認(rèn)為是將實部加上數(shù)值為零的虛部，變成復(fù)信號y(n)=x(n)+j0，從而進(jìn)行FFT計算。MATLAB提供的FFT工具函數(shù)可以實現(xiàn)任意點數(shù)的混合基運算，本文在實際應(yīng)用研究中用C語言實現(xiàn)了混合基FFT算法，作為本文數(shù)據(jù)處理的軟件工具。

1.2 窗函數(shù)

FFT對有限長的信號做頻譜分析時，會產(chǎn)生頻譜泄露現(xiàn)象，原本該頻率的能量會泄漏到其他頻率上，導(dǎo)致頻譜不準(zhǔn)確。當(dāng)輸入信號的頻率不是FFT分辨率的整數(shù)倍時，信號的能量就會向整個頻域擴散，此時原本幅度比較小頻率就會疊加屬于其他頻率的泄露能量，使得原本小幅度的頻率能量無法在頻譜中顯示其真實的信息。因此，需要建立窗函數(shù)與信號疊加，可以在一定程度上防止或減少能譜外泄效應(yīng)。目前常見的窗函數(shù)有矩形窗、漢寧窗、海明窗、三角窗高斯窗等，具體根據(jù)實際數(shù)據(jù)和所要解決的問題而定。

由于窗函數(shù)同樣會造成信號產(chǎn)生不同程度的失真，所以在選擇和應(yīng)用窗函數(shù)時，一方面要考慮窗譜本身的特性，另一方面還要考慮信號本身的特點。如果對頻率分辨率要求較高，而譜幅值精度要求不高，可直接選用主瓣窄的矩形窗。在實際應(yīng)用中也可以選擇多種窗函數(shù)進(jìn)行比較和選擇。對于采集觸發(fā)得到的微震信號，一般能量比較集中，而波形特征又不統(tǒng)一，所以很難選擇一種固定的窗函數(shù)形式，為保證主頻帶能量盡量保留、其他頻帶能夠快速壓制，因此選用一種參數(shù)可設(shè)的凱澤窗(圖1)，通過調(diào)整控制參數(shù)來修改窗的形態(tài)，從而達(dá)到不同的阻帶衰減要求。其中，乘數(shù)因子可采用零階第一類修正貝塞爾函數(shù)。

圖1 凱澤窗Fig.1 Kaiser window

2 FFT頻域濾波

通過FFT可以知道信號序列中包含的頻率成分，以及各頻率成分的振幅大?。环粗ㄟ^快速傅立葉逆變換(inverse fast fourier transform，IFFT)，又可以將頻域的信號變換到時域，從而得到與原來信號長度相同的時間序列。

2.1 FFT計算

輸入實數(shù)序列信號數(shù)據(jù)如圖2所示，經(jīng)過FFT計算的結(jié)果為復(fù)數(shù)序列對稱共軛形式，其對稱實數(shù)域和共軛虛數(shù)域如圖3所示。

圖2 原始微震信號及能量頻譜Fig.2 Primary seismic signal and energy spectrum

圖3 對稱實數(shù)域與共軛虛數(shù)域Fig.3 Symmetric real number domain andconjugate imaginary domain

為保證對頻譜復(fù)數(shù)序列進(jìn)行FFT反變換之后能夠得到實數(shù)序列，需要保留頻譜復(fù)數(shù)序列的共軛對稱特性。以如圖4所示的長度位為9的奇數(shù)頻譜序列X1(0∶8)，和長度為10的偶數(shù)頻譜序列X2(0∶9)為例。X1(0)和X2(0)部分代表直流分量，可直接保留不做處理。序列X1(1∶8)分作X1(1∶4)和X1(5∶8)兩個部分，即這兩個復(fù)數(shù)序列為共軛關(guān)系，實部相對稱。序列X2(1∶9)分作X2(1∶4)和X2(6∶9)兩個部分，這兩個復(fù)數(shù)序列滿足共軛關(guān)系，其中實部對稱，X2(5)為對稱中心直接保留即可。對滿足以上條件的頻譜序列進(jìn)行FFT反變換，直接取其實部，即可得到濾波后的實數(shù)序列。

圖4 奇數(shù)和偶數(shù)長度序列Fig.4 Odd & even length sequence

2.2 頻域濾波

通過FFT可以在頻域中對信號的部分頻段做濾波處理，再通過IFFT得到時間序列信號。值得注意的是，在進(jìn)行IFFT時需要將奈奎斯特頻率之后對應(yīng)的頻點位置通過共軛對稱的方法來補全，以得到正確的時域信號。因此整個濾波過程可以分解為五個步驟。

第一步：通過局部極值法獲取譜峰頻點，見式(1)。

S{xi}={xi|f(xi+1)≥

f(xi)}∩{xi|f(xi-1) ≤f(xi)}

(1)

第二步：過濾選取有效的頻譜點。

方法一：在局部極值集合的基礎(chǔ)上多次迭代，見式(2)。

S′{xi}={xi|S(xi+1)≥

S(xi)}∩{xi|S(xi-1)≤S(xi)}?

S″{xi}?…Sr{xi)}

(2)

通過多次迭代計算局部極值，得到最終需要的譜峰頻點集合。迭代次數(shù)根據(jù)頻譜特征來確定，迭代次數(shù)越多，獲得的譜峰越少。如果信號頻譜比較集中可以適當(dāng)減少迭代次數(shù)，反之如果信號頻譜比較復(fù)雜，則可適當(dāng)提高迭代次數(shù)，目的是獲得占主要成分的頻帶分量，并且考慮頻域覆蓋度。

方法二：以能量權(quán)重優(yōu)先為原則，保留能量最大的頻譜分量，見式(3)。

Sr{xi}={xi|f(xi)≥λ}

(3)

式中，λ為可調(diào)閾值，取值可以參考局部極值頻點集合的3/4百分位數(shù)。方法二計算過程較為簡單，但會損失一部分能量較小的頻域分量，因此該方法更適用于頻譜分布比較集中的信號。

第三步：選取合適的頻譜分量帶寬，取值標(biāo)準(zhǔn)可參考1/2峰值處帶寬的2倍。在減少譜峰串?dāng)_的基礎(chǔ)上帶寬選取可以寬一些，以降低信號的能量損失。

第四步：獲取分量并進(jìn)行濾波處理，利用平頂窗函數(shù)加權(quán)計算。加權(quán)系數(shù)S(n)見式(4)。

(4)

式中：n為序號；N為信號長度，0≤n≤N-1。該過程可以提高信號的信噪比，同時會改善信號分量譜峰間串?dāng)_問題。

第五步：頻域數(shù)據(jù)經(jīng)過IFFT得到時域信號分量，完成信號的頻域濾波過程完成。

3 信號濾波與信源提取實例分析

3.1 信號分量粗提取

以某礦山的微震和爆破兩種實際信號為例進(jìn)行說明。在對信號分量進(jìn)行濾波提取之前，首先進(jìn)行頻段的劃分和提取。圖5和圖6中的兩個例子先不考慮信號的實際頻譜特征，大致選取幾個頻段進(jìn)行分量提取。圖5為微震信號實例，分別選取100～400 Hz、200～400 Hz、400～600 Hz三個頻段。其中，圖5(a)為原始微震信號波形，圖5(b)為0～200 Hz頻段分量，圖5(c)為200～400 Hz頻段分量，圖5(d)為400～600 Hz頻段分量。圖6為爆破信號實例，分別選取0～400 Hz和1 000～1 200 Hz兩個頻段。其中，圖6(a)為原始爆破信號波形；圖6(b)為100～400 Hz頻段分量；圖6(c)為1 000～1 200 Hz頻段分量。

圖5 微震信號實例Fig.5 Example for microseismic signal

圖6 爆破信號實例Fig.6 Example for blasting signal

信號能夠?qū)崿F(xiàn)多源分離的關(guān)鍵因素是譜峰識別和帶寬提取，識別有效的譜峰可以提高信號分量的識別能力，這也是信號多源分離的前提；合理提取帶寬能夠提高分量信號的還原度。分析表明，微震信號更偏向低頻，能量時間分布周期較長，信源成分較單一和集中。

對比各信號分量與原始信號可以發(fā)現(xiàn)，圖5(b)分量信號與原始信號尾波十分接近，圖5(c)分量信號與微震主能量時段相吻合，而圖5(d)分量則無明顯的波形特征；圖6(b)與爆破信號形態(tài)接近，圖6(c)為爆破信號的高頻部分，與原始信號相比無明顯特征。通過以上兩個信號分量粗提取的例子可以看到，提取的信號分量特征不明顯，分量信號波形形態(tài)也比較隨機，不能得到統(tǒng)一的分量數(shù)據(jù)，無法有效開展進(jìn)一步的統(tǒng)計分析。

3.2 頻譜細(xì)分

礦山現(xiàn)場的微震或爆破信號比較復(fù)雜，所以很難獲得規(guī)律的一致性特征。在得到信號的頻譜之后識別有效譜峰，認(rèn)為每一個譜峰是一種信源分量，然后按照信號頻譜的譜峰來進(jìn)行頻帶劃分，針對每個頻帶進(jìn)行窄帶濾波，從而提取出信號中的多種信源分量。頻帶劃分見圖7。

圖7 微震信號細(xì)分頻譜Fig.7 Subdivision spectrum for microseismic signal

3.3 信源分量提取

信源分量提取應(yīng)遵守以下原則：①兩個相鄰的抖動譜峰需要合并為一個來提取；②譜峰帶寬需要保留一部分延拓重疊數(shù)據(jù)；③使用合適的窗函數(shù)生成濾波系數(shù)，以降低信號分量失真度。如圖8所示，由以上分解到的時域分量數(shù)據(jù)可以分析得到：①各個分量時域信號觸發(fā)特征明顯；②觸發(fā)到時數(shù)據(jù)一致性；③低頻分量的振幅相對高頻偏大；④能量頻譜分布和時域分量趨勢一致，振幅、能量較大的頻段信號量對應(yīng)的能量譜值也比較大；⑤各個分量的觸發(fā)到時可以為原始完整信號的觸發(fā)到時提供評估基礎(chǔ)。

圖8 多頻帶細(xì)分信號分量Fig.8 Multiband subdivision of signal components

4 結(jié) 語

本文通過FFT頻譜分析和頻域濾波，提出了一種礦山微震多源信號的實際可行的分析與分解方法。本方法以多源信號的譜峰分布為依據(jù)，進(jìn)行信源成分識別，利用FFT窄帶濾波有效實現(xiàn)了信源分量的分離與提取。通過實際數(shù)據(jù)分析實驗結(jié)果表明，在信號源成分未知的情況下能夠靈活有效實現(xiàn)了信號的多源分離，得到信源分量波形。信源分量的觸發(fā)到時與原始信號的一致性，也從另一個角度印證了該方法的有效性。該方法為微震信號的分類識別、觸發(fā)到時評估、能量震級的評估以及進(jìn)一步的信號特征關(guān)聯(lián)分析和樣本庫統(tǒng)計提供了合理有效的數(shù)據(jù)基礎(chǔ)和技術(shù)支持，對礦山微震監(jiān)測領(lǐng)域的信號分析提供了新的思路。