基于自適應(yīng)輪動(dòng)歸類的微震震相識別方法

賈寶新，李峰

遼寧工程技術(shù)大學(xué) 土木工程學(xué)院，遼寧阜新 123000

0 引言

隨著礦山開采深度及規(guī)模的增加，開采所處地下環(huán)境的不確定性逐漸增強(qiáng)，表現(xiàn)為所處地層的原生地質(zhì)構(gòu)造易遭受大范圍擾動(dòng)甚至破壞，導(dǎo)致礦產(chǎn)資源在開采的同時(shí)也會(huì)伴隨著一些諸如沖擊地壓、采空區(qū)垮落、巖爆、巷道片幫冒頂?shù)拳h(huán)境與安全隱患(潘一山等，2007；崔峰等，2019；李金雨等，2019).礦產(chǎn)資源開采過程中普遍采用微震監(jiān)測技術(shù)對災(zāi)害源進(jìn)行監(jiān)測定位，其中微震震相識別方法是微震監(jiān)測的核心內(nèi)容，也是實(shí)施實(shí)時(shí)微震監(jiān)測以及間接提高災(zāi)害預(yù)測精確度的關(guān)鍵(陳澤等，2020；Yu et al.，2020).

微震震相識別方法主要是利用微震信號本身的頻率、功率密度、持續(xù)時(shí)間、振幅等屬性(賈寶新等，2021)，將其中包含的微震響應(yīng)與背景噪聲進(jìn)行區(qū)分，使得微震響應(yīng)可以被自動(dòng)識別并截取.針對微震響應(yīng)識別困難的問題，董隴軍等(2016)選取微震力矩對數(shù)、能量對數(shù)、傳感器觸發(fā)數(shù)量、首次峰值幅值對數(shù)等屬性作為Fisher分類法的特征參數(shù)，以此訓(xùn)練震相識別模型；朱權(quán)潔等(2013)提出了可利用到時(shí)、振幅等信息剔除異常波形的多通道微震聯(lián)合識別方法；姜福興等(2014)提出了一種將預(yù)處理與識別后再優(yōu)化相結(jié)合的微震單事件多通道聯(lián)合識別方法；Chen和Paul(2021)提出了利用量化角頻率、地震矩、應(yīng)力降等震源特性，對微震爆破事件實(shí)施自動(dòng)識別和鑒別的方法；趙明等(2019)利用真實(shí)余震事件和深度學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)(CNN)進(jìn)行模型訓(xùn)練，以此對地震波形進(jìn)行了分類與識別；Zhao等(2015)將微震波形重復(fù)性、尾波下降、信號主頻和發(fā)生時(shí)間作為Fisher分類法特征參數(shù)，以此訓(xùn)練震相識別模型；Dong等(2016)選取微震發(fā)生時(shí)間、地震矩、總輻射能、S波和P波能量比、拐角頻率和靜態(tài)應(yīng)力降作為微震識別的特征參數(shù)，并用Fisher分類法、樸素Bayes法和邏輯回歸進(jìn)行分類；Jiang等(2021)通過奇異值分解減少數(shù)據(jù)量后，利用微震信號時(shí)頻特性提出了基于隨機(jī)森林分類器的震相識別方法；Wu等(2016)提出了基于S變換、相位和頻率的相似性以及隨機(jī)組合分析P波波形的震相識別方法；Zhang等(2019)提出了一種集成經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)、奇異值分解(SVD)和極限學(xué)習(xí)機(jī)(ELM)相結(jié)合的震相識別方法.

以上研究成果在一定程度上提高了震相識別的準(zhǔn)確率，但也存在一些不足.例如改進(jìn)長短時(shí)窗法(short time-window average/long time-window average，簡稱STA/LTA)必須計(jì)算整段信號的長、短時(shí)窗比值，而整段信號中無效的背景噪聲占絕大部分，故其浪費(fèi)了大量的算力與時(shí)間；機(jī)器學(xué)習(xí)方法通常面臨樣本量不足的問題，導(dǎo)致其難以處理多種形式的微震信號；多參數(shù)分析法由于其特征參數(shù)復(fù)雜，導(dǎo)致參數(shù)匹配難以普適，方法的實(shí)際應(yīng)用較為困難.

針對以上問題，本文基于自適應(yīng)高通濾波、背景噪聲幅值上下界計(jì)算、輪動(dòng)圓半徑計(jì)算及輪動(dòng)迭代、超限點(diǎn)歸類分組4個(gè)信號處理流程提出了自適應(yīng)輪動(dòng)歸類法，并通過分別對模型試驗(yàn)與實(shí)際工程中所得的微震信號進(jìn)行震相識別，證明了自適應(yīng)輪動(dòng)歸類法相較改進(jìn)STA/LTA方法在識別準(zhǔn)確率、識別穩(wěn)定性、計(jì)算速度、計(jì)算穩(wěn)定性等方面均具有一定的優(yōu)勢.自適應(yīng)輪動(dòng)歸類法針對目前已有震相識別方法存在的問題提出了可行的解決方案：相較改進(jìn)STA/LTA方法，該方法只需對微震信號幅值超出限定范圍的超限點(diǎn)進(jìn)行歸類分組，分組結(jié)果即為震相識別結(jié)果，而不必遍歷整段信號并篩選出超出閾值的部分信號，以此節(jié)省了算力，又因其只關(guān)注數(shù)據(jù)點(diǎn)是否超出輪動(dòng)圓范圍，而不需考慮超出部分的具體數(shù)值，因此避免了長短時(shí)窗相比必須在閾值限定范圍內(nèi)的普適性不佳的問題；相較機(jī)器學(xué)習(xí)方法，該方法不需要大量且形式各異的樣本信號進(jìn)行模型訓(xùn)練，其自身的自適應(yīng)高通濾波步驟可將不同形式的微震信號濾波為微震響應(yīng)較為突出的信號形式，隨即以此為基礎(chǔ)并通過之后的處理步驟對微震響應(yīng)進(jìn)行識別；相較多參數(shù)分析法，該方法需要主動(dòng)調(diào)試的初始賦值參數(shù)僅為3個(gè)，分別為分組間隔Fint(s)、范圍擴(kuò)大系數(shù)Q以及采樣間隔Cint，并且其取值方法獨(dú)立且定義明確，不需要多次嘗試尋找最優(yōu)組合，因此該方法在保證識別質(zhì)量的前提下也可適應(yīng)廣泛的應(yīng)用場景.

1 自適應(yīng)輪動(dòng)歸類法基本原理

利用背景噪聲與微震響應(yīng)之間的關(guān)系，自適應(yīng)輪動(dòng)歸類法可以從微震信號中篩選出微震響應(yīng)，其原理是先計(jì)算可涵蓋絕大部分背景噪聲的特征值，再以該值的系數(shù)倍為界限找出大于該界限的所有微震信號數(shù)據(jù)點(diǎn)，即超限點(diǎn).最后以一定的標(biāo)準(zhǔn)對超限點(diǎn)進(jìn)行分類，所得結(jié)果即為震相識別結(jié)果.

1.1 自適應(yīng)輪動(dòng)歸類法運(yùn)行原理

自適應(yīng)輪動(dòng)歸類法包含以下4個(gè)主要步驟：

(1)自適應(yīng)高通濾波

為了降低背景噪聲對于微震響應(yīng)的干擾，并提高自適應(yīng)輪動(dòng)歸類法的識別準(zhǔn)確率與算法穩(wěn)定性，需對原始信號進(jìn)行高通濾波，而高通濾波所需的高通濾波下限需要人工輸入，故為了實(shí)現(xiàn)整體方法的自動(dòng)流程化運(yùn)行，需要一種可以處理多種形式微震信號的自適應(yīng)高通濾波處理方法.本步驟由自適應(yīng)高通濾波下限計(jì)算以及對原始信號進(jìn)行高通濾波兩部分組成，其中后者采用Matlab軟件中的巴特沃斯濾波器實(shí)現(xiàn)，以下主要說明由本文提出的自適應(yīng)高通濾波下限的計(jì)算方法.

背景噪聲通常由多種不同頻段、不同能量的噪聲信號組成，本步驟是以其中能量較大且頻率穩(wěn)定的噪聲為準(zhǔn)，利用高通濾波過濾該頻率及以下的信號，以此保留較高頻率的微震響應(yīng)信號以及一部分高頻率但能量低的噪聲信號，用以將原始信號中的微震響應(yīng)信號凸顯出來.從原始信號中隨機(jī)篩選出時(shí)長相同但起始時(shí)刻不同的多段截取信號，通過計(jì)算各截取信號的基礎(chǔ)過零率并將其中位數(shù)作為整段信號中背景噪聲主頻的估計(jì)值，以此獲得高通濾波下限并對原始信號進(jìn)行高通濾波.

針對高通濾波的自動(dòng)化運(yùn)行，本文提出了基礎(chǔ)過零率，利用周期函數(shù)判別周期的方式(單位時(shí)間內(nèi)自下而上穿過X軸次數(shù))來估算截取信號中能量較高且頻率穩(wěn)定的主要趨勢波形的頻率.此估算頻率對結(jié)果的精確度要求不高，因?yàn)殡S機(jī)截取多段信號中，其中可能包含微震響應(yīng)、施工設(shè)備響應(yīng)等干擾信號，并且背景噪聲在整段信號中的占比是最大的，所以為了減小誤差并找到目標(biāo)背景噪聲，需要截取多段信號并取其基礎(chǔ)過零率的中位數(shù)，以保證所得結(jié)果可以代表微震信號中能量較高且頻率穩(wěn)定的背景噪聲.

截取信號段數(shù)即分樣組數(shù)Ldz計(jì)算見后(5)式所示，基礎(chǔ)過零率計(jì)算方式如圖1所示，圖中10 s內(nèi)信號自下而上穿過X軸次數(shù)為9次，故基礎(chǔ)過零率為0.9 Hz.

圖1 基礎(chǔ)過零率計(jì)算示意圖

自適應(yīng)高通濾波前后的語譜圖對比如圖2所示.可見圖中288.018 Hz以下的低頻背景噪聲已被準(zhǔn)確過濾，圖中使用的微震原始信號及其自適應(yīng)高通濾波效果如圖3所示.

圖2 自適應(yīng)高通濾波前后的語譜圖對比

(2)背景噪聲幅值上下界計(jì)算

為了進(jìn)一步凸顯出微震響應(yīng)信號并初步確定輪動(dòng)圓半徑Rw的取值依據(jù)，需要找到可以基本涵蓋全部背景噪聲幅值的范圍，確保此范圍內(nèi)基本為背景噪聲，而范圍外基本為微震響應(yīng)，此范圍的上界與下界即為背景噪聲幅值上下界.

在原始信號經(jīng)過自適應(yīng)高通濾波提高信噪比后，引入正態(tài)分布中橫軸區(qū)間(μ-2σ,μ+2σ)的落點(diǎn)概率95%作為信號的無效置信區(qū)間，利用自Y=0逐步增加涵蓋范圍并判斷是否已包含整段信號中95%的數(shù)據(jù)點(diǎn)的方法，迭代比較出滿足上述條件的背景噪聲幅值上下界，確保所得結(jié)果可以涵蓋大部分背景噪聲及無效響應(yīng)并初步篩選出微震響應(yīng)信號.

為了更好的適應(yīng)不同信號特征變化，規(guī)定算法采集微震信號時(shí)間間隔Tcol(s)為

Tcol=|--[Fint/(1-95%)]/10?×10，

(1)

該公式限定了采集微震信號時(shí)間間隔Tcol為10的倍數(shù)，其最小值為10 s，并且可以根據(jù)具體的微震響應(yīng)持續(xù)時(shí)間以10的倍數(shù)進(jìn)行調(diào)整.這樣不僅保證了微震響應(yīng)時(shí)長較短情況下算法的適配性，也可適應(yīng)諸如礦震、地震等微震響應(yīng)時(shí)長較長的情況，方便了后續(xù)的算法計(jì)算.

自適應(yīng)高通濾波與背景噪聲幅值上下界計(jì)算示意圖如圖3所示.

圖3 自適應(yīng)高通濾波與背景噪聲幅值上下界計(jì)算示意圖

(3)輪動(dòng)圓半徑計(jì)算及輪動(dòng)迭代

為了篩選出少量可代表單次微震響應(yīng)的較高幅值數(shù)據(jù)點(diǎn)，為下一步超限點(diǎn)歸類分組提供基礎(chǔ)數(shù)據(jù)，使用輪動(dòng)圓半徑Rw組成的輪動(dòng)圓自整段信號的起點(diǎn)時(shí)刻至終點(diǎn)時(shí)刻，遍歷所有數(shù)據(jù)點(diǎn)并篩選出幅值大于輪動(dòng)圓半徑Rw的數(shù)據(jù)點(diǎn)，即超限點(diǎn).

少量且有較高幅值的數(shù)據(jù)點(diǎn)有利于歸類分組，而且可以較為突出的區(qū)別于背景噪聲，因?yàn)檩^高幅值數(shù)據(jù)點(diǎn)是由能量突變的微震響應(yīng)所造成的，其可區(qū)別于全程存在但能量較低的背景噪聲，并且由于能量的突變，微震響應(yīng)區(qū)域的幅值要高于背景噪聲區(qū)域且所含數(shù)據(jù)點(diǎn)要少于背景噪聲區(qū)域，所以篩選出少量且幅值較高的數(shù)據(jù)點(diǎn)更加符合微震響應(yīng)的特征.

輪動(dòng)圓圓心位于X軸上，其半徑取值為背景噪聲幅值上界Du與背景噪聲幅值下界Dd間的距離的范圍擴(kuò)大系數(shù)Q倍，通過遍歷整段信號所有數(shù)據(jù)點(diǎn)并與其進(jìn)行比較，篩選出大于輪動(dòng)圓半徑Rw的數(shù)據(jù)點(diǎn)，即可獲得下一步驟所需的所有超限點(diǎn).迭代篩選時(shí)采用整段信號中前后兩個(gè)元素絕對值的平均值來代表此處信號的幅值，這樣可以增加迭代比較的適用范圍，避免出現(xiàn)單個(gè)數(shù)據(jù)點(diǎn)激變導(dǎo)致的識別錯(cuò)誤.輪動(dòng)圓半徑Rw的計(jì)算公式如式(2)所示，超限點(diǎn)篩選條件如式(3)所示，輪動(dòng)圓篩選超限點(diǎn)示意圖如圖4所示.

圖4 輪動(dòng)圓篩選超限點(diǎn)示意圖

(4)超限點(diǎn)歸類分組

為了最終獲取微震信號中微震響應(yīng)的影響區(qū)域即對整段微震信號進(jìn)行震相識別，需要將輪動(dòng)圓篩選出的超限點(diǎn)進(jìn)行歸類分組，分為同一組代表一次震相識別結(jié)果，分為多組則代表多次震相識別結(jié)果.

超限點(diǎn)作為微震響應(yīng)的代表少量分布于整段微震信號中，由于不同次微震響應(yīng)之間存在一定的時(shí)間差，所以可以利用一個(gè)較大的時(shí)間間隔對超限點(diǎn)進(jìn)行歸類分組，若多個(gè)超限點(diǎn)之間的時(shí)間間隔不超過規(guī)定的時(shí)間間隔，則可分為一組，此組即為一次震相識別結(jié)果.輪動(dòng)圓半徑計(jì)算及輪動(dòng)迭代步驟中選出的全部超限點(diǎn)，將其按分組間隔Fint進(jìn)行歸類，若前后兩者之間的時(shí)間差大于分組間隔Fint則歸為一類，反之歸為下一類，如是選出所有類，即為所求震相識別結(jié)果.超限點(diǎn)分類效果(自適應(yīng)高通濾波后)即自適應(yīng)輪動(dòng)歸類法識別震相效果如圖5所示，圖中紅色點(diǎn)劃線矩形框?yàn)檎鹣嘧R別結(jié)果.

圖5 超限點(diǎn)歸類分組效果示意圖

1.2 自適應(yīng)輪動(dòng)歸類法的實(shí)現(xiàn)

(1)算法主線運(yùn)行流程

① 讀取微震信號數(shù)據(jù)S2(mm·s-1)與對應(yīng)時(shí)間數(shù)據(jù)S1(s)，并執(zhí)行②.

② 輸入以下數(shù)據(jù)：

采樣頻率Fs(Hz)，采樣間隔Cint，分組間隔Fint，范圍擴(kuò)大系數(shù)Q=2，背景噪聲幅值上界Du=0，背景噪聲幅值下界Dd=0，并執(zhí)行③.

③ 按照采樣間隔Cint對微震信號數(shù)據(jù)S2與對應(yīng)時(shí)間數(shù)據(jù)S1進(jìn)行數(shù)據(jù)規(guī)模縮減，并執(zhí)行④.

④ 計(jì)算自適應(yīng)高通濾波下限Hp(Hz)，隨后將所得結(jié)果作為高通濾波下限對微震信號數(shù)據(jù)S2進(jìn)行高通濾波，并執(zhí)行⑤.

⑤ 計(jì)算微震信號數(shù)據(jù)S2升序排列中前95%的元素個(gè)數(shù)并將其賦值給過濾總數(shù)Jnum，并執(zhí)行⑥.

⑥ 將0賦值給計(jì)數(shù)因子j，依次判斷微震信號數(shù)據(jù)S2內(nèi)元素是否同時(shí)滿足小于背景噪聲幅值上界Du且大于背景噪聲幅值下界Dd，若滿足將j+1賦值給計(jì)數(shù)因子j，否則執(zhí)行⑦.

⑦ 判斷計(jì)數(shù)因子j是否大于過濾總數(shù)Jnum，若滿足執(zhí)行⑨，否則執(zhí)行⑧.

⑧ 將Du+1賦值給背景噪聲幅值上界Du，將Dd-1賦值給背景噪聲幅值下界Dd，并執(zhí)行⑥.

⑨ 計(jì)算輪動(dòng)圓半徑Rw并依次判斷其是否滿足以下公式：

Rw=(Du+|Dd|)×Q,

(2)

Rw<(Ji+Ji+1)/2,

(3)

式中：i為索引序數(shù)，取值為1至微震信號數(shù)據(jù)S2的最大元素個(gè)數(shù)；Ji為迭代因子，表示微震信號數(shù)據(jù)S2中第i元素個(gè)元素的絕對值.

若滿足將迭代因子Ji的值賦值給待歸類矩陣P1的第i個(gè)元素，否則將0賦值給待歸類矩陣P1的第i個(gè)元素，并執(zhí)行⑩.

(2)自適應(yīng)高通濾波下限計(jì)算運(yùn)行步驟

① 讀取微震信號數(shù)據(jù)S2與對應(yīng)時(shí)間數(shù)據(jù)S1，輸入采樣頻率Fs，并執(zhí)行②.

② 計(jì)算分組間隔內(nèi)元素總數(shù)Lint，并執(zhí)行③.其中分組間隔內(nèi)元素總數(shù)Lint滿足以下公式：

Lint=|--Fint×Fs?,

(4)

③ 計(jì)算分樣組數(shù)Ldz以及微震信號數(shù)據(jù)S2內(nèi)元素總個(gè)數(shù)Gm，并執(zhí)行④.其中分樣組數(shù)Ldz滿足以下公式：

Ldz=?(Send-Sfirst)/(1/|--(Fs/Lint)?)_|,

(5)

式中，Sfirst為對應(yīng)時(shí)間數(shù)據(jù)S1中的首個(gè)元素的值(s)；Send為對應(yīng)時(shí)間數(shù)據(jù)S1中的末尾元素的值(s).

④ 計(jì)算出Ldz組時(shí)間索引數(shù)a1(s)與信號索引數(shù)a2(mm·s-1)，并將時(shí)間索引數(shù)a1的值賦值給分組矩陣P3第1行中對應(yīng)位置處的元素，將信號索引數(shù)a2的值賦值給分組矩陣P3第2行中對應(yīng)位置處的元素.其中a1與a2滿足以下公式：

a1=?Rand×Gm_|,

(6)

a2=a1-(Lint-1),

(7)

式中，Rand為0至1范圍內(nèi)的隨機(jī)數(shù).

⑤ 依次判斷分組矩陣P3每組數(shù)據(jù)中的信號索引數(shù)a2是否小于1，若滿足則執(zhí)行⑥，否則執(zhí)行⑦.

⑥ 按照式(6)與式(7)重新計(jì)算該組數(shù)據(jù)，直至信號索引數(shù)a2大于1，隨后將計(jì)算結(jié)果賦值給分組矩陣P3對應(yīng)位置處的元素，并執(zhí)行⑤.

⑦ 依次將分組矩陣P3中每組信號索引數(shù)a2與時(shí)間索引數(shù)a1索引微震信號數(shù)據(jù)S2與對應(yīng)時(shí)間數(shù)據(jù)S1中的各部分?jǐn)?shù)據(jù)記錄在數(shù)據(jù)記錄矩陣P4中，并執(zhí)行⑧.

⑧ 依次判斷數(shù)據(jù)記錄矩陣P4中每組數(shù)據(jù)中的微震信號數(shù)據(jù)的前后兩個(gè)元素是否滿足前者小于0且后者大于0，并統(tǒng)計(jì)出每組數(shù)據(jù)中滿足條件的元素的總個(gè)數(shù)，將其賦值給過零率記錄矩陣P5中對應(yīng)組號的元素，并執(zhí)行⑨.

⑨ 將過零率記錄矩陣P5中的每個(gè)元素除以其對應(yīng)時(shí)間數(shù)據(jù)的總時(shí)長，隨后計(jì)算出過零率記錄矩陣P5的中位數(shù)，即為所求的自適應(yīng)高通濾波下限.

自適應(yīng)輪動(dòng)歸類法流程圖如圖6所示.

圖6 自適應(yīng)輪動(dòng)歸類法流程圖

1.3 自適應(yīng)輪動(dòng)歸類法初始參數(shù)分析

本算法需要初始賦值的參數(shù)共有4個(gè)，分別是采樣頻率Fs、采樣間隔Cint、分組間隔Fint、范圍擴(kuò)大系數(shù)Q.

(1)采樣頻率

采樣頻率Fs據(jù)由監(jiān)測設(shè)備決定，后續(xù)數(shù)據(jù)處理操作不可對其進(jìn)行更改，根據(jù)已有研究成果(賈寶新等，2017；賈寶新和李峰，2022)，采樣頻率Fs越大，可捕獲的微震響應(yīng)細(xì)節(jié)越多.

(2)采樣間隔

采樣間隔Cint用于對信號進(jìn)行規(guī)?？s減，比如間隔10個(gè)數(shù)據(jù)點(diǎn)采樣一次，采樣間隔Cint即為10.這種縮減雖可通過減少計(jì)算量實(shí)現(xiàn)大幅提高計(jì)算速度，但卻會(huì)使原始信號的細(xì)節(jié)減少，由于微震響應(yīng)信號持續(xù)時(shí)間較短，這種采樣方式容易漏掉振幅等關(guān)鍵點(diǎn).故采樣間隔的選取應(yīng)當(dāng)根據(jù)震相識別整體的計(jì)算速度來進(jìn)行調(diào)整，應(yīng)當(dāng)在保證震相識別全過程計(jì)算時(shí)間小于截取信號時(shí)長的前提下進(jìn)行調(diào)整，這樣既能滿足實(shí)時(shí)監(jiān)測的要求，又能保證原始信號細(xì)節(jié)的不丟失.

(3)分組間隔

分組間隔Fint用于對超限點(diǎn)進(jìn)行分組，該值的選擇可根據(jù)微震響應(yīng)持續(xù)時(shí)間來進(jìn)行規(guī)定，比如鋼球撞擊室內(nèi)試驗(yàn)?zāi)Ｐ彤a(chǎn)生的微震響應(yīng)信號時(shí)長約為0.1 s(自子波振幅激變至恢復(fù)平穩(wěn)狀態(tài))，便可將分組間隔Fint設(shè)置為0.1 s.若是工程實(shí)際的微震響應(yīng)信號，可通過試測或者統(tǒng)計(jì)之前的微震持續(xù)時(shí)間，以此確定分組間隔Fint的選取，比如實(shí)測遼寧某煤礦礦震產(chǎn)生的微震信號持續(xù)時(shí)間約為10 s，此時(shí)分組間隔Fint便應(yīng)當(dāng)取值為10 s.

微震響應(yīng)信號中凡是幅值大于輪動(dòng)圓半徑Rw的數(shù)據(jù)點(diǎn)均會(huì)被篩選出來，所以在較高采樣頻率下，代表同一微震響應(yīng)的超限點(diǎn)數(shù)量是眾多的，它們之間的時(shí)間差遠(yuǎn)小于微震響應(yīng)整體的持續(xù)時(shí)間.故分組間隔Fint的選取其實(shí)并不敏感，只需確定不同工況下微震響應(yīng)大致的持續(xù)時(shí)間即可.

(4)范圍擴(kuò)大系數(shù)

范圍擴(kuò)大系數(shù)Q用于擴(kuò)大由背景噪聲幅值上下界產(chǎn)生的背景噪聲估計(jì)范圍.因?yàn)?5%置信區(qū)間的存在，可能有一部分非微震響應(yīng)信號穿過此范圍而被輪動(dòng)圓篩選為超限點(diǎn)，從而被算法錯(cuò)誤的認(rèn)為是微震響應(yīng)，故采用將背景噪聲幅值上下界產(chǎn)生的背景噪聲估計(jì)范圍乘以范圍擴(kuò)大系數(shù)Q來進(jìn)一步擴(kuò)大輪動(dòng)圓的篩選范圍，所得結(jié)果即為輪動(dòng)圓半徑Rw.范圍擴(kuò)大系數(shù)Q可根據(jù)現(xiàn)場的信號干擾源的情況進(jìn)行選取，該值的選取范圍通常在1.0～3.0之間，若在實(shí)驗(yàn)室這種比較安靜且背景噪聲影響小的環(huán)境下，該值可取1.0，若有像放炮或大型工程器械震動(dòng)干擾的情況下，該值可相應(yīng)增大.該值的選取并不需要十分精確，對現(xiàn)場收集的微震信號進(jìn)行試算識別后即可確認(rèn).

2 模型試驗(yàn)

震相識別的難點(diǎn)在于不同類型與不同組合的震相疊加導(dǎo)致難以采用明確的標(biāo)準(zhǔn)進(jìn)行識別(金標(biāo)，2020；賈寶新等，2020)，故現(xiàn)通過實(shí)驗(yàn)室搭建室內(nèi)模型并采用人工發(fā)震的方式，對自適應(yīng)輪動(dòng)歸類法的識別效果進(jìn)行驗(yàn)證.

2.1 試驗(yàn)?zāi)Ｐ蜆?gòu)建

利用石英砂、石膏、石灰、水以6∶3∶1∶2的比例混合后分層堆砌，構(gòu)建室內(nèi)模型，其尺寸為2.3 m×1 m×1.3 m.22個(gè)微震傳感器在按設(shè)計(jì)位置嵌入墻體，其三維空間坐標(biāo)如表1所示.圖7中黑色球體與對應(yīng)數(shù)字代表傳感器位置與標(biāo)號，黑色立方體與對應(yīng)字母為震源位置與標(biāo)號.

圖7 試驗(yàn)布局圖

表1 各檢波器三維空間坐標(biāo)

2.2 監(jiān)測設(shè)備

超高頻多通道構(gòu)造活動(dòng)監(jiān)測儀(Antenna-III)可監(jiān)測頻帶寬度0～50 kHz的微震信號，其采樣頻率最高可達(dá)100 kHz.設(shè)備參數(shù)見表2，超高頻多通道構(gòu)造活動(dòng)監(jiān)測儀各組成部分見圖8.

表2 超高頻多通道構(gòu)造活動(dòng)監(jiān)測儀設(shè)備參數(shù)

圖8 超高頻多通道構(gòu)造活動(dòng)監(jiān)測儀

2.3 試驗(yàn)過程

模型搭建完成后帶模具養(yǎng)護(hù)7天，每天定時(shí)灑水，以防止產(chǎn)生大的貫通裂隙.拆模后保持實(shí)驗(yàn)室通風(fēng)且有每天至少5個(gè)小時(shí)的日照時(shí)間，照此持續(xù)30天且確保模型已干透后，隨即可開始下一步試驗(yàn).因模型干透后具有一定的脆性，所以當(dāng)其被鋼球撞擊時(shí)會(huì)使表面產(chǎn)生一定的下陷與破裂，此現(xiàn)象可確保產(chǎn)生來自擠壓破壞的P波以及來自剪切破壞的S波.

本試驗(yàn)采用半徑1 cm的鋼球?qū)δＰ瓦M(jìn)行定長單擺撞擊作為發(fā)震源，單擺定長為30 cm，固定鐵架與鋼球之間采用棉質(zhì)線連接.鋼球撞擊位置分別為“震源A”至“震源F”，其三維空間坐標(biāo)如表3所示.

表3 各震源三維空間坐標(biāo)

試驗(yàn)時(shí)先將鋼球自然懸掛并剛好接觸模型表面的指定撞擊位置，隨后將鋼球拉至垂直于模型表面一定距離后放手，每間隔2 s左右進(jìn)行1次撞擊，以3次撞擊為1次數(shù)據(jù)記錄，形成1組包含由22個(gè)檢波器產(chǎn)生的22段微震信號的數(shù)據(jù)組.每個(gè)震源點(diǎn)重復(fù)進(jìn)行以上過程3次，產(chǎn)生3組共66段微震信號，試驗(yàn)現(xiàn)場照片如圖9所示.

圖9 試驗(yàn)現(xiàn)場

3 自適應(yīng)輪動(dòng)歸類法識別效果分析

通過整理模型試驗(yàn)所得微震信號數(shù)據(jù)，共得到6個(gè)震源與22個(gè)檢波器下不同信噪比的18組微震信號，即396段微震信號，其中每段信號包含3次微震響應(yīng).將1段信號作為一個(gè)整體對其包含的3次微震響應(yīng)進(jìn)行同時(shí)識別，若識別微震響應(yīng)次數(shù)為3次，則算作1次成功識別.選取震相識別中常用的STA/LTA方法作為對照，通過二者在模型試驗(yàn)微震信號數(shù)據(jù)下震相識別效果的對比，對自適應(yīng)輪動(dòng)歸類法的優(yōu)越性進(jìn)行驗(yàn)證.

3.1 定義微震信號信噪比

(1)計(jì)算原理

為區(qū)分信號質(zhì)量并保證信號的合理分類，定義了微震信號專用的微震信號信噪比(Signal-to-noise ratio of microseismic signal，簡稱SNRmss)(dB).該值利用的是微震信號中震相信號振幅的平均值與自適應(yīng)高通濾波過程中各截取信號振幅的最小值進(jìn)行比較得出，相較傳統(tǒng)信噪比定義該方法在背景噪聲大小的界定上更為準(zhǔn)確且適應(yīng)性強(qiáng).SNRmss由以下公式得出：

SNRmss=10×lg(Ssig/Snoi),

(8)

式中，Ssig為微震響應(yīng)特征值，取值為微震信號數(shù)據(jù)S2中的震相信號振幅的平均值(mm·s-1)；Snoi為背景噪聲特征值，取值為自適應(yīng)高通濾波過程中篩選出的Ldz段截取信號振幅的最小值(mm·s-1).

圖10a為各段截取信號的振幅大小，取其中位數(shù)即為Snoi取值，取微震信號數(shù)據(jù)S2中的微震響應(yīng)信號振幅的平均值即為Ssig的取值，前者如圖10b中黑色短劃線所示，后者如圖10b中綠色短劃線所示，對比圖中的原始微震信號可見二者能有效的代表背景噪聲與微震響應(yīng).

圖10 SNRmss原理示意圖

(2)模型試驗(yàn)下SNRmss的計(jì)算及分析

計(jì)算模型試驗(yàn)所得微震信號數(shù)據(jù)的SNRmss，并對結(jié)果進(jìn)行分類，所得箱型圖如圖11所示.圖中橫坐標(biāo)代表的微震信號序號，例如“A-1”表示的是震源A下第1組微震信號，其中包含22個(gè)檢波器采集的22段微震信號，每段微震信號包含3次微震響應(yīng).

圖11 模型試驗(yàn)下微震信號的SNRmss

自圖11可知，模型試驗(yàn)中不同震源下微震信號的SNRmss的差別較大，以其中位數(shù)為代表，變化范圍為0.099～7.818 dB，相同震源下微震信號的SNRmss的變化范圍為0.012～3.744 dB，說明發(fā)震方式與監(jiān)測環(huán)境對微震信號能否有效監(jiān)測的影響較大.同時(shí)，SNRmss差別較大的試驗(yàn)信號也為下一步對兩種方法震相識別效果的區(qū)分提供了有效條件.

3.2 改進(jìn)STA/LTA方法數(shù)據(jù)分析

(1)STA/LTA方法原理及改進(jìn)方法

作為一種判斷信號瞬時(shí)響應(yīng)的方法，STA/LTA方法常用于微震到時(shí)拾取及震相識別中.其中STA表征為微震信號因微震激勵(lì)到達(dá)而產(chǎn)生的信號激變趨勢；LTA表征為信號背景噪聲的變化趨勢(劉晗和張建中，2014).該方法的原理如下所示(劉曉明等，2017)：

(9)

(10)

(11)

式中，h為采樣時(shí)刻(s)；p為時(shí)窗內(nèi)迭代時(shí)刻(s)；Wl為LTA時(shí)窗(s)；Ws為STA時(shí)窗(s)；λ為響應(yīng)閾值；CF(p)為微震信號在p時(shí)刻對應(yīng)的特征函數(shù)值.

選擇特征函數(shù)為微震信號的振幅絕對值在時(shí)窗內(nèi)的疊加(賈寶新等，2021)，即：

(12)

(13)

(14)

式中，|A(p)|為p時(shí)刻微震信號振幅的絕對值.

因?yàn)镾TA/LTA方法只能篩選出信號幅值激變點(diǎn)而無法識別出震相，故將該方法與自適應(yīng)輪動(dòng)歸類法中的超限點(diǎn)歸類分組模塊進(jìn)行結(jié)合，組成改進(jìn)STA/LTA方法.

(2)計(jì)算速度與識別偏差分析

定義識別偏差為每組微震信號識別結(jié)果與微震響應(yīng)次數(shù)真值的差的絕對值和.

統(tǒng)計(jì)改進(jìn)STA/LTA方法對各段微震信號進(jìn)行震相識別所用的時(shí)間，隨后對結(jié)果進(jìn)行分析，統(tǒng)計(jì)結(jié)果如圖12所示.

由圖12可知，改進(jìn)STA/LTA方法的計(jì)算時(shí)間在0.630 s至1.532 s之間，平均計(jì)算時(shí)間為0.869 s，并且圖中每組數(shù)據(jù)之間的中位數(shù)與分布區(qū)間差別較大，表明該方法計(jì)算速度穩(wěn)定性較差，即對不同類型的微震信號所用的計(jì)算時(shí)間差別較大.

圖12 改進(jìn)STA/LTA方法計(jì)算時(shí)間

使用改進(jìn)STA/LTA方法對模型試驗(yàn)微震信號進(jìn)行震相識別，識別震相次數(shù)結(jié)果如圖13a所示.

由圖13a可知，改進(jìn)STA/LTA方法的震相識別結(jié)果波動(dòng)性相對較大，有106段微震信號的識別次數(shù)出現(xiàn)了識別偏差，并且識別偏差幅度最大為3次，識別準(zhǔn)確率為73.2%，識別偏差標(biāo)準(zhǔn)差為0.766，可見改進(jìn)STA/LTA方法在識別準(zhǔn)確率上欠佳，并且該方法的識別穩(wěn)定性較差且偏差幅度波動(dòng)大.

3.3 自適應(yīng)輪動(dòng)歸類法數(shù)據(jù)分析

使用自適應(yīng)輪動(dòng)歸類法對模型試驗(yàn)微震信號進(jìn)行震相識別，識別震相次數(shù)結(jié)果如圖13b所示.據(jù)圖可知，該方法的識別結(jié)果相對穩(wěn)定，有7段微震信號的識別次數(shù)出現(xiàn)了偏差，但識別偏差幅度最大僅為1次，識別準(zhǔn)確率為98.2%，識別偏差標(biāo)準(zhǔn)差為0.133，可見自適應(yīng)輪動(dòng)歸類法的識別穩(wěn)定性與準(zhǔn)確率均較為優(yōu)秀.

圖13 模型試驗(yàn)下兩種方法震相識別結(jié)果

統(tǒng)計(jì)自適應(yīng)輪動(dòng)歸類法對各段微震信號進(jìn)行震相識別所用的時(shí)間，利用箱型圖對結(jié)果進(jìn)行分析，統(tǒng)計(jì)結(jié)果如圖14所示.

圖14 自適應(yīng)輪動(dòng)歸類法計(jì)算時(shí)間

由圖14可知，自適應(yīng)輪動(dòng)歸類法的計(jì)算時(shí)間在0.069 s至0.789 s之間，平均計(jì)算時(shí)間為0.384 s，并且圖中每組數(shù)據(jù)之間的中位數(shù)與分布區(qū)間差別不大，表明該方法計(jì)算速度穩(wěn)定性較好，即對不同類型的微震信號所用的計(jì)算時(shí)間較為穩(wěn)定.

3.4 兩種方法的震相識別效果對比

改進(jìn)STA/LTA方法與自適應(yīng)輪動(dòng)歸類法對于震相識別的效果對比如圖13所示.試驗(yàn)中3次微震響應(yīng)為1段信號，若3次均識別成功，記作1次成果的震相識別，由圖13a可知，改進(jìn)STA/LTA方法成功識別的微震信號為290段，因?yàn)閿?shù)據(jù)總段數(shù)為396段，故其識別準(zhǔn)確率為73.2%，同理可得自適應(yīng)輪動(dòng)歸類法的識別準(zhǔn)確率為98.2%.統(tǒng)計(jì)兩種方法對模型試驗(yàn)數(shù)據(jù)的各項(xiàng)計(jì)算結(jié)果，所得二者在5項(xiàng)指標(biāo)下的識別效果對比如表4所示.

由表4可知，模型試驗(yàn)所得微震信號數(shù)據(jù)下，相較改進(jìn)STA/LTA方法，自適應(yīng)輪動(dòng)歸類法的識別準(zhǔn)確率提高了25.0%，識別偏差標(biāo)準(zhǔn)差為前者的17.4%，計(jì)算時(shí)間平均值為前者的44.2%，計(jì)算時(shí)間標(biāo)準(zhǔn)差為前者的67.4%.由此可見，自適應(yīng)輪動(dòng)歸類法在識別準(zhǔn)確率、識別穩(wěn)定性、計(jì)算速度、計(jì)算穩(wěn)定性等方面均優(yōu)于改進(jìn)STA/LTA方法.

表4 模型試驗(yàn)下兩種方法震相識別效果對比

3.5 影響因素相關(guān)性分析

計(jì)算兩種震相識別方法對各組微震信號的識別偏差，計(jì)算結(jié)果與SNRmss的關(guān)系如圖15所示.

圖15 SNRmss與識別偏差的關(guān)系

因監(jiān)測設(shè)備的采樣頻率是固定的100 kHz，所以將信號長度定義為微震信號文件的大小，該值與計(jì)算時(shí)間的關(guān)系如圖16所示.

圖16 信號長度與計(jì)算速度的關(guān)系

兩種方法的計(jì)算速度與識別偏差和兩項(xiàng)對比指標(biāo)的相關(guān)性系數(shù)如表5所示.

由表5可知，改進(jìn)STA/LTA方法的計(jì)算速度與識別準(zhǔn)確率高度依賴于微震信號長度，而自適應(yīng)輪動(dòng)歸類法則弱化了對該指標(biāo)的依賴，其兩項(xiàng)相關(guān)系數(shù)分別為改進(jìn)STA/LTA方法的86.1%及47.8%，說明本文方法對比改進(jìn)STA/LTA方法對信號長度的敏感度更低，可以保證對于長度不同信號的識別容錯(cuò)率與穩(wěn)定性更高.

表5 兩種方法的相關(guān)性分析結(jié)果

理論上信噪比越高，震相識別的準(zhǔn)確率應(yīng)當(dāng)越高，但震相識別方法對于SNRmss的相關(guān)性也不宜過低或過高，若過低會(huì)導(dǎo)致方法的識別準(zhǔn)確率忽高忽低，識別偏差難以掌握；過高會(huì)導(dǎo)致方法過于依賴高信噪比信號，使得面對因背景噪聲振幅大且變化復(fù)雜的微震信號的適應(yīng)性較差.然而實(shí)際工程中因?yàn)槠渌麢C(jī)械運(yùn)作等原因，微震信號往往信噪比較低，對SNRmss有過高依賴性的方法往往無法投入實(shí)際使用.本文統(tǒng)計(jì)的識別偏差為單組微震信號中各段信號的識別偏差絕對值的總和，所以識別偏差應(yīng)該與SNRmss呈負(fù)相關(guān)關(guān)系，故此本文以兩項(xiàng)指標(biāo)與SNRmss相關(guān)系數(shù)為-0.5為最佳，用以評判兩種震相識別方法的效果.

自表5中可知，自適應(yīng)輪動(dòng)歸類法的計(jì)算速度和識別偏差與SNRmss的相關(guān)系數(shù)較改進(jìn)STA/LTA方法更為接近-0.5，說明本文方法對SNRmss的依賴性適中，而改進(jìn)STA/LTA方法的兩項(xiàng)相關(guān)系數(shù)接近于0，說明該方法識別偏差不穩(wěn)定分布且后期識別準(zhǔn)確率提升空間有限.

綜上所述，自適應(yīng)輪動(dòng)歸類法在處理信號上的容錯(cuò)率與穩(wěn)定性較改進(jìn)STA/LTA方法更高，并且能夠更好的適應(yīng)且不過度依賴高信噪比的信號，后期可通過軟、硬件提高信號信噪比或者通過微震信號質(zhì)量控制方法篩選高信噪比信號，使其進(jìn)一步的提高震相識別的準(zhǔn)確率.

4 工程驗(yàn)證

為了進(jìn)一步驗(yàn)證自適應(yīng)輪動(dòng)歸類法對于微震信號震相識別的可靠性，選取遼寧某煤礦微震監(jiān)測工程中的微震信號進(jìn)行驗(yàn)證.

4.1 微震監(jiān)測工程介紹

遼寧某煤礦微震監(jiān)測工程的目的是為預(yù)防礦井沖擊地壓等動(dòng)力災(zāi)害，該項(xiàng)目共設(shè)置8個(gè)監(jiān)測臺陣，每個(gè)監(jiān)測臺陣布置1臺三分向微震檢波器，用以全天候監(jiān)測礦區(qū)范圍內(nèi)發(fā)生的微震事件.將7次不同時(shí)期與不同發(fā)震位置的微震事件所產(chǎn)生的微震信號進(jìn)行分類，將其震源分別命名為“震源I”至“震源O”，其中每次事件包含7個(gè)三分向檢波器的X軸、Y軸、Z軸3個(gè)分向的微震信號，每個(gè)分向的微震信號包含1次微震響應(yīng)，信號形式如圖17所示.

圖17 微震監(jiān)測工程中的微震信號示意圖

4.2 兩種方法的震相識別效果對比

下述分析中自適應(yīng)輪動(dòng)歸類法采用初始參數(shù)：采樣頻率Fs=5000 Hz，采樣間隔Cint=10，分組間隔Fint=10 s，范圍擴(kuò)大系數(shù)Q=2.

(1)識別結(jié)果分析

通過日常對礦區(qū)微震事件的監(jiān)測，共收集到7個(gè)震源、8個(gè)檢波器、3個(gè)分向的168段微震信號，將其用“震源-分向”的形式表示，例如“I-X”表示的是震源I下檢波器X軸分向監(jiān)測到的微震信號.利用改進(jìn)STA/LTA方法與自適應(yīng)輪動(dòng)歸類法對所有微震信號進(jìn)行震相識別，其識別結(jié)果如圖18所示，圖中每個(gè)長方體代表1段微震信號，其中每段微震信號的微震響應(yīng)次數(shù)真值為1次，縱軸識別震相次數(shù)代表的是對應(yīng)震相識別方法識別每段微震信號所得的微震響應(yīng)次數(shù)，識別震相次數(shù)與微震響應(yīng)次數(shù)真值之間的差的絕對值稱作識別偏差.

圖18 兩種方法的震相識別結(jié)果

由圖18a可知，改進(jìn)STA/LTA方法識別失敗次數(shù)為9次，識別偏差幅度最大為1次，識別失敗結(jié)果分布不規(guī)律.由圖18b可知，自適應(yīng)輪動(dòng)歸類法識別失敗次數(shù)為1次，識別偏差幅度最大為1次.

(2)計(jì)算時(shí)間分析

統(tǒng)計(jì)改進(jìn)STA/LTA方法與自適應(yīng)輪動(dòng)歸類法對上述168段微震信號進(jìn)行震相識別所用時(shí)間，所得結(jié)果如圖19所示.

圖19 兩種方法的計(jì)算時(shí)間

由圖19可知，改進(jìn)STA/LTA方法的計(jì)算時(shí)間總體要高于自適應(yīng)輪動(dòng)歸類法，其計(jì)算時(shí)間平均值為0.281 s，并且不同震源下的計(jì)算時(shí)間變化幅度較大，極差為0.303 s.相較改進(jìn)STA/LTA方法，自適應(yīng)輪動(dòng)歸類法方法的計(jì)算時(shí)間較短，計(jì)算時(shí)間平均值為0.115 s，并且計(jì)算時(shí)間變化范圍也較小，極差為0.152 s.

(3)兩種方法識別效果對比

通過對識別結(jié)果與計(jì)算時(shí)間的整理與分析，結(jié)果如表6所示.

表6 兩種方法震相識別效果對比

由表6可知，相較改進(jìn)STA/LTA方法，自適應(yīng)輪動(dòng)歸類法的識別準(zhǔn)確率提高了4.8%，識別偏差標(biāo)準(zhǔn)差為前者的34.2%，計(jì)算時(shí)間平均值為前者的40.9%，計(jì)算時(shí)間標(biāo)準(zhǔn)差為前者的55.9%，即該方法相較改進(jìn)STA/LTA方法在震相識別準(zhǔn)確率、穩(wěn)定性以及計(jì)算速度、穩(wěn)定性上均具有一定的優(yōu)勢.

5 結(jié)論

本文提出了用于微震信號震相識別的自適應(yīng)輪動(dòng)歸類法，通過該方法與改進(jìn)STA/LTA方法對模型試驗(yàn)所得微震信號震相識別效果的對比，分析了自適應(yīng)輪動(dòng)歸類法的優(yōu)越性，并在工程實(shí)踐中獲得驗(yàn)證，所得結(jié)論如下：

(1)模型試驗(yàn)中不同震源與相同震源下微震信號SNRmss的中位數(shù)變化范圍分別為0.099～7.818 dB及0.012～3.744 dB，說明發(fā)震方式與監(jiān)測環(huán)境對微震信號的有效監(jiān)測影響較大.

(2)自適應(yīng)輪動(dòng)歸類法的計(jì)算時(shí)間、識別偏差與信號長度的相關(guān)系數(shù)分別為改進(jìn)STA/LTA方法的86.1%及47.8%，并且計(jì)算時(shí)間、識別偏差與SNRmss的相關(guān)系數(shù)更為接近-0.5，說明本文方法可以保證對長度不同信號的識別容錯(cuò)率與穩(wěn)定性更高，并且能夠更好的適應(yīng)且不過度依賴高SNRmss的信號，保證后期有較大的識別準(zhǔn)確率提升空間.

(3)模型試驗(yàn)下，相較改進(jìn)STA/LTA方法，自適應(yīng)輪動(dòng)歸類法的識別準(zhǔn)確率提高了25.0%，識別偏差標(biāo)準(zhǔn)差為前者的17.4%，計(jì)算時(shí)間平均值、標(biāo)準(zhǔn)差分別為前者的44.2%、67.4%.實(shí)際工程下，本文方法的識別準(zhǔn)確率提高了4.8%，識別偏差標(biāo)準(zhǔn)差為前者的34.2%，計(jì)算時(shí)間平均值、標(biāo)準(zhǔn)差分別為前者的40.9%、55.9%.