本發(fā)明涉及微震監(jiān)測技術領域。具體涉及一種工程尺度下微震信號及p波初至自動識別算法，該方法可廣泛用于礦業(yè)工程、水利水電工程、石油工程、巖土工程以及地下工程。

背景技術：
微震源定位是微震監(jiān)測與災害預警的重要組成部分，而微震信號的識別及P波初至拾取是微震源定位的關鍵?，F(xiàn)有的微震信號識別方法大部分來源于地震領域，方法很多，主要有：根據(jù)在時間域能量和能量變化構建特征函數(shù)時間域的STA/LTA算法；根據(jù)地震信號到達前后地震波形數(shù)據(jù)統(tǒng)計性質的差別，如AIC算法；此外還有神經(jīng)網(wǎng)絡法、高階統(tǒng)計法、以及基于小波理論的小波變換法等。STA/LTA及其改進算法，由于算法簡單、計算效率高、適于實時處理，在地震領域被廣泛應用。這些方法用于地震信號的處理具有較好的適用性，但工程尺度下產(chǎn)生的巖石破裂信號與天然地震信號不同：1)地震信號頻率一般為幾十赫茲以下，而微震信號的頻率一般為幾十到幾百赫茲，有的高達幾千赫茲；2)地震信號持續(xù)時間較長，一般大于1秒，微震信號持續(xù)時間較短，一般不超過0.1秒；3)微震與天然地震相比，一般震級小，較小的微震事件不易被拾??；4)工程環(huán)境下由于機械施工、電氣干擾、車輛行駛等因素造成背景噪音復雜，使得監(jiān)測到的大多數(shù)微震信號信噪比較低。因此，將地震領域的信號識別及P波到時拾取方法用于工程尺度微震信號是不合適的，主要表現(xiàn)為：1)難以識別低信噪比的微弱信號；2)P波初至自動拾取誤差較大。為此，針對上述不足，發(fā)明一種可以有效識別工程尺度下微震信號及其p波初至的自動識別的方法，提高工程尺度下微震信號及其P波初至的拾取準確率，進而提高巖爆、礦震、塌方等地質災害預警的及時性與準確性，是非常必要的。

技術實現(xiàn)要素：
本發(fā)明的目的在于克服現(xiàn)有技術存在的問題，提供了一種工程尺度下微震信號及p波初至自動識別算法，提高工程尺度下微震信號及其P波初至的拾取準確率，進而提高巖爆、礦震、塌方等地質災害預警的及時性與準確性。一種工程尺度下微震信號及p波初至自動識別算法，包括以下步驟：步驟1、讀取振動傳感器實時監(jiān)測到的設定時窗內波形的采樣點數(shù)據(jù)，建立平面直角坐標系，其中X軸對應波形的采樣點編號，Y軸對應波形的幅值，在平面直角坐標系中對時窗內波形的采樣點進行對稱校準；步驟2、計算采樣點個數(shù)N，設置零交叉點個數(shù)上限M1，零交叉點個數(shù)下限M2；設置微震信號時間判斷閥值T；步驟3、設定第A個采樣點的加權因子K(A)和特征函數(shù)CF(A)；設定第A個采樣點短時平均值STA(A)和長時平均值LTA(A)；設定第A個采樣點的動態(tài)閾值r(A)；其中A為采樣點個數(shù)的編號；步驟4、以平面直角坐標系中編號最小的采樣點為當前采樣點，當前采樣點的編號為i；步驟5、設置算法參數(shù)變量M、S、L、t，M、S、L、t的初值均設置為0，其中M為零交叉點個數(shù)；S為計數(shù)器；L的計算公式為L＝3+M/3；t為拾取出的微震信號的持續(xù)時間；步驟6、比較當前采樣點的短時平均值STA(i)和動態(tài)閾值r(i)的大??；步驟7、若當前采樣點的短時平均值STA(i)小于當前采樣點的動態(tài)閥值r(i)，則當前采樣點不是P波初至點，此時將當前采樣點的下一個采樣點作為當前采樣點，返回到步驟(6)；若當前采樣點的短時平均值STA(i)大于等于當前點的動態(tài)閥值r(i)，則將當前采樣點的編號i賦值給k，計算當前采樣點的幅值P1，將當前采樣點的幅值大小P1賦值給緩存最大幅值P0，然后進行步驟(8)；步驟8、計算編號為k的采樣點的下一個采樣點k+1的幅值P2；步驟9、如果采樣點的編號為k的采樣點的下一個采樣點k+1不是零交叉點，將緩存最大幅值P0和采樣點k+1的幅值P2中較大值賦值給緩存最大幅值P0，將k+1賦值給k，當k的值小于波形采樣點個數(shù)N時，返回步驟8；當k的值大于等于波形采樣點個數(shù)N時，則將當前采樣點的下一個采樣點作為當前采樣點，返回步驟5；如果采樣點的編號為k的下一個采樣點k+1是零交叉點，則零交叉點個數(shù)M值加1，將k+1賦值給k，然后進行步驟10；步驟10、若零交叉點個數(shù)M值大于零交叉點個數(shù)上限M1值，則當前采樣點不為P波初至點，則將當前采樣點的下一個采樣點作為當前采樣點，返回步驟5；若零交叉點個數(shù)M值小于等于零交叉點個數(shù)上限M1值，則計算當前采樣點的判定參數(shù)δ＝LTA(i)×M，然后進行步驟11；步驟11、比較編號為k的采樣點的短時平均值STA(k)與當前采樣點的判定參數(shù)δ的大?。喝艟幪枮閗的采樣點的短時平均值STA(k)小于等于當前采樣點的判定參數(shù)δ，計數(shù)器S值歸零，將P2賦值給P0，此時將k加1，返回步驟8；若編號為k的采樣點的短時平均值STA(k)大于當前采樣點的判定參數(shù)δ，計數(shù)器S值加1，計算L＝3+M/3，進入步驟12；步驟12、比較計數(shù)器S值與L的大小：若計數(shù)器S值小于等于L，則將P2賦值給P0，將k加1，返回步驟8；若計數(shù)器S值大于L，計算從當前采樣點i到編號為k的采樣點的這段波形的持續(xù)時間t，t＝(k-i)/f，其中f為采樣頻率；步驟13、當t小于微震信號時間判斷閥值T或零交叉點個數(shù)M小于零交叉點個數(shù)下限M2時，則當前采樣點i與編號為k的采樣點之間的采樣點不是微震信號，將編號為k的采樣點的下一個采樣點作為當前采樣點，清零M、P0、P1、P2、S、L的取值，進行步驟5；當t大于等于時間閥值T且零交叉點個數(shù)M大于等于零交叉點個數(shù)下限M2時，則當前采樣點i與編號為k的采樣點之間的采樣點為微震信號，將當前采樣點的編號i賦值給T1，將k賦值給T2，其中編號為T1的采樣點為微震信號的P波初至點，編號為T2的采樣點為信號結束點，進入步驟14；步驟14、將編號為k的采樣點的下一個采樣點作為當前采樣點，清零M、P0、P1、P2、S、L的取值，進行步驟5，直至將設定時窗內波形的采樣點數(shù)據(jù)分析完。如上所述的對稱校準包括以下步驟：選取時窗中設定個數(shù)的連續(xù)的采樣點作為校正樣本，若校正樣本中幅值為正的采樣點個數(shù)與幅值為負的采樣點個數(shù)的比值不在設定比值范圍內，則設定時窗內波形的采樣點數(shù)據(jù)需對稱校準，計算校正樣本幅值的平均值，將設定時窗內波形的采樣點幅值與校正樣本幅值的平均值相加，完成對稱校準；若校正樣本中幅值為正的采樣點個數(shù)與幅值為負的采樣點個數(shù)的比值不在設定比值范圍內，則設定時窗內波形的采樣點數(shù)據(jù)不需對稱校準。如上所述的第A個采樣點的加權因子K(A)基于以下公式：所述的第A個采樣點的特征函數(shù)CF(A)基于以下公式：其中，A為對稱校準后波形中采樣點個數(shù)的編號，y(A)為對稱校準后波形的第A個采樣點的幅值，為y(A)一階差分；所述的第A個采樣點的短時平均值STA(A)基于以下公式：STA(A)＝STA(A-1)+C3×[CF(A)-STA(A-1)]所述的第A個采樣點的長時平均值LTA(A)基于以下公式：LTA(A)＝LTA(A-1)+C4×[CF(A)-LTA(A-1)]所述的第A個采樣點的動態(tài)閾值r(A)基于以下公式：r(A)＝C5×LTA(A)其中，C3為短時平均系數(shù)，C4為長時平均系數(shù)，C5為觸發(fā)閥值。當P波到來時，將引起微震信號振幅或者頻率等波形特征變化，若信號的信噪比越低，則P波初至越不明顯。因而提高拾取算法對該類變化的敏感性，將能有利于提高工程尺度下微震信號的P波初至拾取能力。而算法中的加權因子K、特征函數(shù)CF、動態(tài)長短時平均比值ε均能反應波形的特征變化，直接影響拾取準確率，動態(tài)長短時平均比值ε為長時平均值LTA(A)與短時平均值STA(A)的比值：分析本發(fā)明的加權因子K、特征函數(shù)CF、動態(tài)長短時平均比值ε對信號振幅或者頻率波形等特征變化的敏感性，分別如圖2、3、4。對比結果分析本發(fā)明三種變量在頻率或振幅變化處數(shù)值改變明顯，達到峰值的迅速，因此在理論上本發(fā)明對信噪比低的信號及其P波初至的拾取效果明顯。從錦屏引水隧洞實時監(jiān)測的數(shù)據(jù)中隨機選取555個微震信號為樣本，定義兩個拾取結果判斷準則：1)信號成功拾?。耗軐颖拘盘栕詣幼R別出來而不產(chǎn)生漏波與誤判。這樣將有利于人工對自動拾取出來的信號部分進行進一步處理。信號成功拾取個數(shù)越多，信號拾取率越高。2)P波初至拾取準確：將自動算法拾取P波初至結果與人工拾取結果相比，若拾取誤差不大于3個采樣點，則認為該算法自動拾取結果準確。本發(fā)明信號拾取率為94.23％，P波拾取準確率為73.51％，拾取率完全滿足微震領域目前定位誤差。分析2015年6月15日至2015年7月15日錦屏極深地下實驗室微震監(jiān)測數(shù)據(jù)，一共524個微震信號，109個事件，將人工拾取P波、本發(fā)明拾取結果用于定位，結果分別如圖6、7所示，其中7#實驗室是巖石破裂高風險區(qū)域，可以看出本發(fā)明定位結果微震事件的聚集度較高，微震事件空間分布規(guī)律與人工拾取結果接近，與工程實際情況相符。將人工拾取過程中信噪比非常低、P波不明顯、定位結果明顯異常的13個事件，如圖6中的三角形所示，用本發(fā)明拾取定位分析：1)本發(fā)明拾取的P波結果均能成功定位，如圖7中的三角形所示；2)本發(fā)明該13個事件拾取定位效果相對于人工拾取改進明顯，絕大部分在定位誤差允許的范圍內。有益效果：本發(fā)明針對工程尺度下微震信號特點，改進地震領域的識別算法，本發(fā)明能夠準確的自動拾取工程尺度下的微震信號特別是性噪比較低的弱信號，同時提高了自動拾取P波初至的準確率，進而提高巖爆、礦震、塌方等地質災害預警的及時性與準確性。附圖說明圖1為本發(fā)明流程圖；圖2為本發(fā)明加權因子K對振幅或頻率變化敏感性分析圖；如圖2(a)當振幅發(fā)生變化時，本發(fā)明的加權因子振幅約增大2倍；如圖2(b)當頻率發(fā)生變化時，本發(fā)明K值變化幅度較為明顯。圖3為本發(fā)明特征函數(shù)CF對振幅或頻率變化敏感性分析圖；圖3(a)模擬背景噪音中出現(xiàn)幅值較小的微震信號，即幅值由1突變到2時，本發(fā)明特征函數(shù)在突變點的幅值由1到16，突變前后波形對稱軸由Y＝1變?yōu)閅＝16；本發(fā)明在突變前后幅值變化和對稱軸偏移明顯，特征函數(shù)曲線突變迅速，圖3(b)模擬信號出現(xiàn)頻率變化時，本發(fā)明特征函數(shù)在突變點的幅值由1到15.6，幅值變化明顯，幅值大小增長迅速。圖4為本發(fā)明ε比值對振幅或頻率變化敏感性分析圖；圖4(a)所示：本發(fā)明ε變化曲線在波形振幅變化后3個采樣內即達到峰值，曲線變化明顯，；如圖4(b)當頻率變化時，人工能準確拾取出本發(fā)明頻率改變位置，曲線變化明顯。圖5為本發(fā)明實例1拾取過程S-L值變化圖；圖6人工拾取微震信號P波定位效果左視圖；圖7本發(fā)明自動拾取微震信號P波定位效果左視圖。具體實施方式為了使本發(fā)明的技術手段、創(chuàng)作特征、工作流程、使用方法達成目的與功效易于明白了解，下面結合具體實施例，進一步闡述本發(fā)明。本發(fā)明的保護范圍不受以下實例的限制。錦屏地下實驗室垂直巖石覆蓋達2400米，是目前世界巖石覆蓋最深的實驗室，工程區(qū)的最大主應力可達63MPa，屬典型的高應力區(qū)，其中7#、8#實驗室工程施工時已遭遇到強～極強巖爆，給現(xiàn)場人員和設備安全造成了嚴重威脅和損失?，F(xiàn)場微震監(jiān)測過程中發(fā)現(xiàn)，由于傳感器距離巖爆風險較高的7-8#實驗室施工段距離較遠，大概為400～500m，微震信號在傳播過程中衰減嚴重，且現(xiàn)場監(jiān)測背景噪音復雜，監(jiān)測到的微震信號大多數(shù)信噪比低，能量幅值小，P波初至不明顯，致使微震信號不易識別，P波不易拾取。利用該工程微震樣本信號，將本發(fā)明拾取和人工拾取結果對比分析，利用粒子群算法優(yōu)化出適于該工程的本發(fā)明最佳參數(shù)組：C3＝0.25，C4＝0.08，C5＝2，上限M1＝300，下限M2＝15，T＝0.01。本實例以分析該工程某時窗內實時監(jiān)測采樣點數(shù)據(jù)為例加以說明。實例1：一種工程尺度下微震信號及p波初至自動識別算法，流程如圖1所示，其步驟如下：(1)、讀取振動傳感器實時監(jiān)測到的設定時窗內波形的采樣點數(shù)據(jù)，建立平面直角坐標系X0Y，其中X軸對應波形的采樣點編號，從平面直角坐標系的原點到X軸正軸方向編號從1開始，按照編號從小到大排列采樣點，采樣點編號為正整數(shù)，Y軸對應波形的幅值，在平面直角坐標系中對時窗內波形的采樣點進行對稱校準，進入步驟(2)。波形對稱校準：選取時窗中設定個數(shù)(在本例中為1000個)連續(xù)的采樣點作為校正樣本，若校正樣本中幅值為正的采樣點個數(shù)與幅值為負的采樣點個數(shù)的比值不在設定比值范圍內(在本例中為小于0.95或大于1.05)，則選取的時窗內波形的采樣點數(shù)據(jù)需對稱校準，計算校正樣本幅值的平均值，將步驟(1)中設定時窗內波形的采樣點幅值與校正樣本幅值的平均值相加，完成對稱校準，這樣有利于準確統(tǒng)計波形中的零交叉點個數(shù)M，減小算法拾取誤差。若校正樣本中幅值為正的采樣點個數(shù)與幅值為負的采樣點個數(shù)的比值在設定比值范圍內(在本例中為大于等于0.95小于等于1.05)，則步驟(1)中設定時窗內波形的采樣點數(shù)據(jù)不需對稱校準。(2)計算對稱校準后平面直角坐標系中波形的采樣點個數(shù)N。設置零交叉點個數(shù)上限M1＝300，零交叉點個數(shù)下限M2＝15；設置微震信號時間判斷閥值T＝0.01，單位為秒。微震信號持續(xù)時間短，設置零交叉點個數(shù)上限M1和零交叉點個數(shù)下限M2，能有效避免持續(xù)時間較短的電氣尖峰信號干擾和對大部分持續(xù)時間較長的非微震信號誤判。(3)、利用對稱校準后平面直角坐標系中的波形的采樣點數(shù)據(jù)計算本發(fā)明的第A個采樣點的加權因子K(A)：和第A個采樣點的特征函數(shù)CF(A)：其中A為對稱校準后波形中采樣點個數(shù)的編號，y(A)為對稱校準后波形的第A個采樣點的幅值，為y(A)一階差分。計算短時平均值STA(A)和長時平均值LTA(A)：STA(A)＝STA(A-1)+C3×[CF(A)-STA(A-1)]LTA(A)＝LTA(A-1)+C4×[CF(A)-LTA(A-1)]其中STA(A)為第A個采樣點的短時平均值；C3為短時平均系數(shù)，實例1中C3取值為0.25；LTA(A)為第A個采樣點長時平均值；C4為長時平均系數(shù)，實例1中C4取值為0.08。計算r(A)＝C5×LTA(A)，r(A)為第A個采樣點的動態(tài)閥值，C5為觸發(fā)閥值，實例1中C5值為取2。(4)、以平面直角坐標系中編號最小的采樣點為當前采樣點，當前采樣點的編號為i。(5)、設置算法參數(shù)變量M、S、L、t，初值均設置為0，其中M為零交叉點個數(shù)；S為計數(shù)器，起計數(shù)作用；L的計算公式為L＝3+M/3；t為拾取出的微震信號的持續(xù)時間。(6)、比較當前采樣點的短時平均STA(i)值和動態(tài)閥值r(i)的大小。(7)、若當前采樣點的短時平均值STA(i)小于當前采樣點的動態(tài)閥值r(i)，則當前采樣點不是P波初至點，此時將當前采樣點的下一個采樣點作為當前采樣點，返回到步驟(6)；若當前采樣點的短時平均值STA(i)大于等于當前點的動態(tài)閥值r(i)，則當前采樣點可能為P波初至點，將當前采樣點的編號i賦值給k，即為k＝i，計算當前采樣點的幅值大小P1，為找出當前采樣點后第一個峰值，將當前采樣點的幅值大小P1賦值給緩存最大幅值P0暫存為最大峰值，然后進行步驟(8)。(8)、計算對稱校準后坐標系波形中采樣點的編號為k的下一個采樣點k+1的幅值P2。(9)、如果對稱校準后波形中采樣點的編號為k的下一個采樣點k+1不是零交叉點，比較緩存最大幅值P0、采樣點k+1的幅值P2的大小，將較大值賦值給緩存最大幅值P0，將k+1賦值給k，當k的值小于波形采樣點個數(shù)N時，返回步驟(8)；當k的值大于等于波形采樣點個數(shù)N時，則將當前采樣點的下一個采樣點作為當前采樣點，即i＝i+1，返回步驟(5)。如果對稱校準后坐標系波形中采樣點的編號為k的下一個采樣點k+1是零交叉點，則零交叉點個數(shù)M值加1，即M＝M+1，將k+1賦值給k，然后進行步驟(10)。(10)、判斷若零交叉點個數(shù)M值大于零交叉點個數(shù)上限M1值，則當前采樣點不為P波初至點，則將當前采樣點的下一個采樣點作為當前采樣點，即i＝i+1，返回步驟(5)。若零交叉點個數(shù)M值小于等于零交叉點個數(shù)上限M1值，則計算當前采樣點的判定參數(shù)δ＝LTA(i)×M，然后進行步驟(11)。(11)、比較編號為k的采樣點的短時平均值STA(k)與當前采樣點的判定參數(shù)δ的大?。喝艟幪枮閗的采樣點的短時平均值STA(k)小于等于當前采樣點的判定參數(shù)δ，計數(shù)器S值歸零，將P2賦值給P0，此時將k加1，即k＝k+1，返回步驟(8)；若編號為k的采樣點短時平均值STA(k)大于當前采樣點的判定參數(shù)δ，計數(shù)器S值加1，即S＝S+1，計算L＝3+M/3，進入步驟(12)。(12)、比較計數(shù)器S值與L的大?。喝粲嫈?shù)器S值小于等于L，則將P2賦值給P0，此時將k加1，即k＝k+1，返回步驟(8)；若計數(shù)器S值大于L，計算對稱校準后坐標系波形中從當前采樣點到編號為k的采樣點的這段波形的持續(xù)時間t，t＝(k-i)/f，其中f為采樣頻率。(13)、當t小于微震信號時間判斷閥值T，T為0.01，單位為秒，或零交叉點個數(shù)M小于零交叉點個數(shù)下限M2時，M2為15，則判斷當前采樣點i與編號為k的采樣點之間的采樣點不是微震信號，需對此時的采樣點k后的信號繼續(xù)拾取，此時將對稱校準后波形坐標系中編號為k的采樣點的下一個采樣點作為當前采樣點，清零M、P0、P1、P2、S、L的取值，進行步驟(5)；當t大于等于微震信號時間閥值T且零交叉點個數(shù)M大于等于零交叉點個數(shù)下限M2時，則判斷當前采樣點i與編號為k的采樣點之間的采樣點為微震信號，將i賦值給T1，k賦值給T2，其中編號為T1的采樣點為微震信號的P波初至點，編號為T2的采樣點為信號結束點，進入步驟(14)。(14)、繼續(xù)對下一個信號自動分析，對稱校準后波形坐標系中編號為k的采樣點的下一個采樣點作為當前采樣點，清零M、P0、P1、P2、S、L的取值，進行步驟(5)，直至將設定時窗內波形的采樣點數(shù)據(jù)分析完。以上顯示和描述了本發(fā)明的基本原理和主要特征和本發(fā)明的優(yōu)點。本行業(yè)的技術人員應該了解，本發(fā)明不受上述實施例的限制，上述實施例和說明書中描述的只是說明本發(fā)明的原理，在不脫離本發(fā)明精神和范圍的前提下，本發(fā)明還會有各種變化和改進，這些變化和改進都落入要求保護的本發(fā)明范圍內。本發(fā)明要求保護范圍由所附的權利要求書及其等效物界定。