基于頻域濾波的礦山微震信號分離方法研究

吳彥博，張 達，冀 虎，戴 銳，石雅倩

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

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

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

1 基礎理論

1.1 傅立葉變換及其反變換

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

1.2 窗函數(shù)

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

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

為保證對頻譜復數(shù)序列進行FFT反變換之后能夠得到實數(shù)序列，需要保留頻譜復數(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)兩個部分，即這兩個復數(shù)序列為共軛關系，實部相對稱。序列X2(1∶9)分作X2(1∶4)和X2(6∶9)兩個部分，這兩個復數(shù)序列滿足共軛關系，其中實部對稱，X2(5)為對稱中心直接保留即可。對滿足以上條件的頻譜序列進行FFT反變換，直接取其實部，即可得到濾波后的實數(shù)序列。

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

2.2 頻域濾波

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

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

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

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

(1)

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

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

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

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

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

(2)

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

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

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

(3)

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

第三步：選取合適的頻譜分量帶寬，取值標準可參考1/2峰值處帶寬的2倍。在減少譜峰串擾的基礎上帶寬選取可以寬一些，以降低信號的能量損失。

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

(4)

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

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

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

3.1 信號分量粗提取

以某礦山的微震和爆破兩種實際信號為例進行說明。在對信號分量進行濾波提取之前，首先進行頻段的劃分和提取。圖5和圖6中的兩個例子先不考慮信號的實際頻譜特征，大致選取幾個頻段進行分量提取。圖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

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

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

3.2 頻譜細分

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

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

3.3 信源分量提取

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

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

4 結 語

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