一種基於相關性分析的礦山微震信號初至波時刻提取方法
2023-05-11 22:13:46
一種基於相關性分析的礦山微震信號初至波時刻提取方法
【專利摘要】一種基於相關性分析的礦山微震信號初至波時刻提取方法，屬於礦山微震信號的分析處理方法。其功能為智能識別兩信號時間差、實現自動移位對齊和判斷礦山微震信號相關性，該方法主要分為五個步驟：一是讀取礦山微震信號；二是求解兩礦山微震信號的相關函數；三是根據求解出的相關函數，求出最大相關時的時間差；四是根據計算出的時間差對其中一個信號進行移位；五是對對齊後的兩礦山微震信號求解其相關係數，判斷其相關性。本發明較好的解決了智能識別礦山微震信號到達時間差及判斷礦山微震信號是否來自同一震源的問題。
【專利說明】—種基於相關性分析的礦山微震信號初至波時刻提取方法
【技術領域】
[0001]本發明涉及礦山微震信號的分析處理方法，特別涉及於一種基於相關性分析的礦山微震信號初至波時刻提取方法。
【背景技術】
[0002]在自然地震、工程爆破、爆炸等震動事件中，其中有一部分能量必然轉化為震動波的形式，它會以震源為中心向周圍傳播。礦山微震信號初至波拾取是震動研究中非常關鍵和重要的問題。在震源定位中，準確快捷的拾取到初至波的時刻是進行準確震源定位或信號分析的基礎。對採集到的礦山微震信號進行分析和處理，並以此來確定礦山微震信號初至波到達時刻的技術稱為初至波拾取技術。初至波拾取技術在軍事、民用和工業工程領域都具有非常廣泛的應用研究。
[0003]而礦震是採礦活動誘發的礦井巖體突然失穩破壞的動力現象，它嚴重威脅著礦井生產與礦工的生命安全，對其進行實時監測和預警具有重要的理論和現實意義。礦震信號相對於地震信號而言，具有震級小、震源淺和影響範圍有限等特點，所以可以稱之為礦山微震信號。由於可以通過分析礦山微震信號的分布特點及發生頻率大致標註出危險地帶，實現有效的監測與預警，且因礦山微震信號本身能量釋放較小，震波傳輸距離有限，從而要求礦山微震信號的定位精度相比普通地震信號更高。而定位精度的提高又需要初至波時刻更為快速準確的提取。
[0004]目前一般採用的礦山微震信號初至波自動拾取技術並無專門針對礦山微震信號的特性和要求的。初至波自動拾取目的是要確定其信號中純噪聲信號和有效信號之間的分界時刻，通常都是根據礦山微震信號的振幅、頻率和相位的變化來確定這一時刻。傳統的初至拾取方法主要分為兩大類:一類是基於地震記錄瞬時特徵的方法，如極值法(峰值檢測)、差分法；這類方法對噪聲比較敏感，當地震記錄的噪聲較嚴重時，難以準確拾取初至。另一類方法是基於地震記錄整體特徵的方法，如相關法；這類方法雖然對噪聲有較好的抑制作用，但受到地震道之間相似性等因素的影響，對於複雜地震記錄，初至拾取的精度也會受到影響。
[0005]到目前為止，已經提出了許多初至拾取的方法，如人工拾取法、相關法、能量比法、最大振幅法、分形維法及神經網絡法等。
[0006]人工拾取法簡單易行，但受人為因素和主觀因素影響較大，容易引入人為誤差，會直接導致結果誤差增大。
[0007]Gelchinsky和Shtivelman提出了一種用相鄰道進行互相關的方法,它假設各道的脈衝形狀不發生變化。相關法受初至的續至影響較大，並且與子波的選擇關係十分密切，同時要求選擇合適的時窗範圍，這些在實際工程中都有一定的難度。
[0008]Hatherrly提出了線性最小平方預測技術與拐點校正相結合的方法，他提出首先識別第一個峰值和拐點，然後估算二者的統計差值。黃成之等採用統計方法將地震初至波記錄分成信號和噪音兩個部分，並使這兩部分統計特徵之間的差別為最大。統計特徵法收到地震波形相似性影響較大，會對精度造成一定影響。
[0009]能量比值法使用周期內的信號能量與總時窗能量的比值，對初至比較敏感，續至波衰減比較快，所以將比值的最大值點作為初至的近似值並作適當的時移，即為初至時刻。Coppens提出了在不同大小的時窗內進行能量比較的方法。江玉樂等提出同極性能量比值法，即改進的能量比法。由於能量比值法的抗幹擾能力還不夠好，所以對於初至波形發生明顯變化的地區拾取的初至時刻不夠準確。
[0010]時間域分形維方法拾取初至的過程必須插值，且結果強烈的依靠插值的準確性。Fab1 Boschetti等提出了一種基於分形維的初至檢測算法，該方法是基於地震道隨著信號出現其分形維值發生變化的特徵來確定地震道初至。但其對時窗和步長的選取十分敏感，稍有不慎就會嚴重的影響其結果。
[0011]神經網絡法利用多參數特徵進行模式識別，充分利用地震記錄的瞬時特徵和整體特徵。由於神經網絡不僅具有並行處理、自組織自學習能力，而且具有高度魯棒性、容錯性和高度的映射、計算和分類能力。莊東海等採用將地震記錄初至拾取看作一個模式識別過程，充分利用地震記錄的瞬時特徵和整體特徵，用人工神經網絡方法進行地震記錄初至拾取，其能獲得較好的實際效果。但其主要缺點是算法複雜性太高，搜索需要一定時間。
[0012]地震信號通過小波多分辨分解，可以有效地分離、消除噪聲，有利於分形維與神經網絡法提高拾取初至的精度。羅光提出了改進型基於小波變換的初至波拾取方法，楊俊峰提出了基於小波變換的三分向震相識別法、能量因子法。但有時其仍需要人工選取數據段來進行拾取以減少拾取時間。
[0013]除此之外，有一些方法還依賴該道與其近道之間的對比，雖然這一類方法對噪聲有一定的壓製作用，但受到地震道之間相似性等因素的影響，對於複雜地震記錄，初至拾取的精度也會受到影響。
[0014]綜上所述，地震波信號是一個一維的時間序列，由於該序列僅僅表明了時間和幅度之間的關係，且有噪聲信號的幹擾，對直接進行地震波初至的拾取造成了相當大的困難。解決此類問題可採用如上文所述的智能算法，但由於智能算法普遍具有較高的算法複雜度，不適合用於需要快速、並有較高定位精度的礦山微震環境，而本發明較好的解決了對初至波時間的拾取這一問題，並具有較低的算法複雜度。
【發明內容】
[0015]本發明的目的是針對已有技術中存在的問題，提供一種基於相關性分析的礦山微震信號初至波時刻提取方法，解決礦山微震信號處理中初至波時間智能識別、礦山微震信號對齊和判斷礦山微震信號是否來自同一震源的問題。
[0016]實現本發明目的的技術方案:本發明的礦山微震信號初至波時刻提取方法包括五個步驟:一是讀取礦山微震信號；二是求解兩礦山微震信號的相關函數；三是根據求解出的相關函數，求出最大相關時的時間差；四是根據計算出的時間差對其中一個信號進行移位；五是對對齊後的兩礦山微震信號求解其相關係數，判斷其相關性；具體方法步驟如下:
[0017](一)讀取礦山微震信號:將兩礦山微震信號讀取進系統，定義\為一礦山微震信號序列，yn為另一礦山微震信號序列，要求兩序列必須等長，設定其序列長度為L，且兩信號採樣頻率均為fs ；
[0018](二)求解兩礦山微震信號的相關函數:根據下式求解礦山微震信號相關函數，在相關函數取得最大值時，根據採樣點及採樣率計算出礦山微震信號間的時間差；根據離散信號序列的相關函數公式計算相關函數Rxy為:
CO
[0019]尺'.r(m)=Yj x(n) y (11 + m)
?=-00
[0020]其中變量m取值範圍為O至L ;
[0021](三)根據求解出的相關函數，求出最大相關時的時間差:根據計算出的時間差對礦山微震信號進行對齊；根據相關函數、採樣點與採樣率計算最大相關時的時間差，取得相關函數最大值Rmax為:
[0022]Rmax = max [ | Rxy (m) | ]
[0023](四)根據計算出的時間差對其中一個信號進行移位:對對齊後的礦山微震信號進行相關性分析，並設定閾值以判斷礦山微震信號是否來自於同一震源；取得初始時刻至相關函數絕對值最大時刻的採樣點總數，記作N，則時間偏差Irffsrt:
[0024]Toffset = N/fs
[0025](五)對對齊後的兩礦山微震信號求解其相關係數，判斷其相關性:
[0026]兩礦山微震信號的相關係數矩陣Rc^為:
P rP\ I P\ 2 η
[0027]R,?n = YJ
Al Pu
[0028]其中p u = E ((X1-E (Xi)).(Yj-E (Yj)))，其中變量1、變量j取值範圍均為I至2，且X為χη序列的總體，Y為yn序列的總體，E為數學期望；求出的相關係數矩陣主對角線表示自相關性，其始終為I ;副對角線是兩信號的相關性，絕對值越大，相關性越好，其取值範圍為[0，I];分界值可根據實際情況不同設定為O到I之間的任何值；最終判斷程序是否需要結束，若需結束則結束方法，若需繼續運行則返回讀取數據步驟，讀取新數據進行移位對齊和相關性比較。
[0029]有益效果，由於採用了上述方案，智能識別兩信號時間差、實現自動移位對齊和判斷礦山微震信號相關性，實現智能識別礦山微震信號到達時間差；實現礦山微震信號的智能對齊，以便於存儲、分析和處理；通過相關性比較兩礦山微震信號是否來自於同一震源；解決了智能識別礦山微震信號到達時間差及判斷礦山微震信號是否來自同一震源的問題。
[0030]所述礦山微震信號初至波時刻提取方法由於採用相關性分析的方法，算法簡便，相比【背景技術】中的模糊識別方法、時間域分形維方法、神經網絡方法、小波變換方法等方法具有較小的程序複雜度，更適合應用於存在大量震動數據、對及時性要求較高場合的使用，相比人工識別方法、能量比值法等方法減少了主觀因素產生的誤差。
【專利附圖】
【附圖說明】
[0031]此處所說明的附圖用來提供對本發明的進一步理解，構成本申請的一部分，本發明的示意性實施例及其說明用於解釋本發明，並不構成對本發明的不當限定。
[0032]圖1為本發明方法的流程圖。
[0033]圖2為本發明實施例的系統結構示意圖。
[0034]圖3為本發明實施例的a、b、C、d信號波形圖。
[0035]圖4為本發明方法實施例中a信號與b信號相關函數及移位對齊效果仿真圖。
[0036]圖5為本發明方法實施例中a信號與d信號相關函數及移位對齊效果仿真圖。
【具體實施方式】
[0037]下面結合附圖中實施例對本發明作進一步說明:
[0038]礦山微震信號初至波時刻提取方法包括五個步驟:一是讀取礦山微震信號；二是求解兩礦山微震信號的相關函數；三是根據求解出的相關函數，求出最大相關時的時間差；四是根據計算出的時間差對其中一個信號進行移位；五是對對齊後的兩礦山微震信號求解其相關係數，判斷其相關性；具體方法步驟如下:
[0039](一)讀取礦山微震信號:將兩礦山微震信號讀取進系統，定義\為一礦山微震信號序列，yn為另一礦山微震信號序列，要求兩序列必須等長，設定其序列長度為L，且兩信號採樣頻率均為fs ；
[0040](二)求解兩礦山微震信號的相關函數:根據下式求解礦山微震信號相關函數，在相關函數取得最大值時，根據採樣點及採樣率計算出礦山微震信號間的時間差；根據離散信號序列的相關函數公式計算相關函數Rxy為:
OO
[0041 ] R'' (/? ) = Σ-ν(/?) ν(/； + ///)
?=-οο
[0042]其中變量m取值範圍為O至L ;
[0043](三)根據求解出的相關函數，求出最大相關時的時間差:根據計算出的時間差對礦山微震信號進行對齊；根據相關函數、採樣點與採樣率計算最大相關時的時間差，取得相關函數最大值Rmax為:
[0044]Rmax = max [ | Rxy (m) | ]
[0045](四)根據計算出的時間差對其中一個信號進行移位:對對齊後的礦山微震信號進行相關性分析，並設定閾值以判斷礦山微震信號是否來自於同一震源；取得初始時刻至相關函數絕對值最大時刻的採樣點總數，記作N，則時間偏差Irffsrt:
[0046]Toffset = N/fs
[0047](五)對對齊後的兩礦山微震信號求解其相關係數，判斷其相關性:
[0048]兩礦山微震信號的相關係數矩陣Rc^為:
[0049]Korr = V]
P1、Pll
[0050]其中p u = E ((X1-E (Xi)).(Yj-E (Yj)))，其中變量1、變量j取值範圍均為I至2，且X為χη序列的總體，Y為yn序列的總體，E為數學期望；求出的相關係數矩陣主對角線表示自相關性，其始終為I ;副對角線是兩信號的相關性，絕對值越大，相關性越好，其取值範圍為[0，1];分界值可根據實際情況不同設定為O到I之間的任何值。最終判斷程序是否需要結束，若需結束則結束方法，若需繼續運行則返回讀取數據步驟，讀取新數據進行移位對齊和相關性比較。
[0051]本發明利用相關函數對信號進行相關性分析，從而得到初至波時刻、進行信號自動對齊和判斷礦山微震信號是否來自同一震源。
[0052]附圖1為本發明方法的流程圖，所述基於相關性分析的礦山微震信號初至波時刻提取方法，按照附圖1的流程完成時刻提取、信號對齊及判斷信號是否來自同一震源，其基本步驟如下:
[0053]讀取礦山微震信號Xn，礦山微震信號yn，其採樣頻率均為fs;
[0054]根據離散信號序列的相關函數公式計算相關函數:
[0055]Rxy(m)= ^-v(//)v(// + /；7)(1.)
?=—CO
[0056]根據相關函數、採樣點與採樣率計算最大相關時的時間差
[0057]取得相關函數最大值:
[0058]Rmax = max [ | Rxy (m) | ] (2)
[0059]取得初始時刻至相關函數絕對值最大時刻的採樣點總數，記作N
[0060]則時間偏差:
[0061]Toffset = N/fs (3)
[0062]依據上一步驟計算出的Irffset將兩信號移位對齊，以便下一步驟計算相關係數；
[0063]根據公式求兩個信號的相關係數矩陣:
[0064]R = [；" ^ '"I(4)
Pl I P22
[0065]其中Pij = E((X1-E(Xi)).(Yj-E(Yj)))，且X為xn序列的總體，Y為yn序列的總體，E為數學期望。求出的相關係數矩陣主對角線表示自相關性，其始終為I ;副對角線是兩信號的相關性，絕對值越大，相關性越好。一般情況下兩信號相關係數絕對值在0-0.3之間為微相關；0.3-0.5之間為實相關；0.5-0.8之間為顯著相關；0.8-1.0之間為高度相關。該方法通過設定分界值來判定兩信號是否具有高度相關性，從而判斷其是否來自同一震源，該分界值可以根據實際需要進行更改。
[0066]所述基於相關性分析的礦山微震信號初至波時刻提取方法，可應用於網絡化分布式礦山震源定位系統中，附圖2為該系統結構示意圖，下面結合附圖2對所述方法的【具體實施方式】做一詳細說明。附圖2中，圓點為傳感器節點，包含微震傳感器，可以監測礦山微震信號。當震源發生震動時，臨近的傳感器節點可以採集到礦山微震數據。附圖3為收到礦山微震信號的四個不同震動傳感器記錄的波形，其中a、b、c三個信號來自於同一震源，d來自於不同於其它三個信號的另一震源。在本實施例中，根據大量樣本及相關經驗設定判定是否來自同一震源的分界值為0.5。
[0067]首先對a與b信號進行相關性分析。附圖4為本發明方法實施例中a信號與b信號相關函數及移位對齊效果仿真圖。據式(I)計算a與b信號的相關函數並作圖，求得的相關函數具體分布見附圖4。據式(2)取得相關函數最大值Rmax = 250.8038。據式(3)計算a與b信號兩信號時間差，Irffset = -0.0393 (S)。據式(4)計算兩信號相關係數矩陣。在相關係數矩陣中主對角線始終為1，副對角線即兩信號相關性，絕對值越大則相關性越好。求得的a與b信號相關係數矩陣為:
廠 1. -0.5462'
[0068]Rab =
ab -0.5462 1,00
[0069]將其相關係數矩陣副對角線取絕對值，即得a信號與b信號相關性為0.5462，大於設定的分界值0.4，則所述方法判定a信號與b信號來自不同震源。
[0070]其次對a與d信號進行相關性分析。附圖5為本發明方法實施例中a信號與d信號相關函數及移位對齊效果仿真圖。據式(I)計算a與d信號的相關函數並作圖，求得的相關函數具體分布見附圖5。據式(2)取得相關函數最大值Rmax = 31.5554。據式(3)計算a與b信號兩信號時間差，Irffsrt = 0.0533 (s)。據式(4)計算兩信號相關係數矩陣。求得的a與d信號相關係數矩陣為:
Γ 1.0 -0.0566'
[0071]RaiJ - _ 0 0566 丨 0_
[0072]將其相關係數矩陣副對角線取絕對值，即得a信號與d信號相關性為0.0566，小於設定的分界值0.4，則所述方法判定a信號與d信號來自不同震源。
[0073]以上結合附圖對本發明的【具體實施方式】作了說明，但這些說明不能被理解為限制了本發明的範圍，本發明的保護範圍由隨附的權利要求書限定，任何在本發明權利要求基礎上的改動都是本發明的保護範圍。
【權利要求】
1.一種基於相關性分析的礦山微震信號初至波時刻提取方法，其特徵是:礦山微震信號初至波時刻提取方法包括五個步驟:一是讀取礦山微震信號；二是求解兩礦山微震信號的相關函數；三是根據求解出的相關函數，求出最大相關時的時間差；四是根據計算出的時間差對其中一個信號進行移位；五是對對齊後的兩礦山微震信號求解其相關係數，判斷其相關性；具體方法步驟如下: (一)讀取礦山微震信號:將兩礦山微震信號讀取進系統，定義Xn為一礦山微震信號序列，yn為另一礦山微震信號序列，要求兩序列必須等長，設定其序列長度為L，且兩信號採樣頻率均為fs ； (二)求解兩礦山微震信號的相關函數:根據下式求解礦山微震信號相關函數，在相關函數取得最大值時，根據採樣點及採樣率計算出礦山微震信號間的時間差；根據離散信號序列的相關函數公式計算相關函數Rxy為:
C0 RX] (m) =+ "?)
?=-co 其中變量m取值範圍為Ο至L; (三)根據求解出的相關函數，求出最大相關時的時間差:根據計算出的時間差對礦山微震信號進行對齊；根據相關函數、採樣點與採樣率計算最大相關時的時間差，取得相關函數最大值Rmax為:
Rmax = max [ I Rxy (m) | ] (四)根據計算出的時間差對其中一個信號進行移位:對對齊後的礦山微震信號進行相關性分析，並設定閾值以判斷礦山微震信號是否來自於同一震源；取得初始時刻至相關函數絕對值最大時刻的採樣點總數，記作N，則時間偏差TfM:
Toffset = N/fs (五)對對齊後的兩礦山微震信號求解其相關係數，判斷其相關性: 兩礦山微震信號的相關係數矩陣艮。?為:
R = [Pl 1 Pl2] corr LI
Ρ:.ι P12 其中= E((X1-E(Xi)).(Υ」-Ε(Υ」)))，其中變量1、變量j取值範圍均為1至2，且X為xn序列的總體，Y為yn序列的總體，E為數學期望；求出的相關係數矩陣主對角線表示自相關性，其始終為1 ;副對角線是兩信號的相關性，絕對值越大，相關性越好，其取值範圍為[0，1];分界值可根據實際情況不同設定為0到1之間的任何值；最終判斷程序是否需要結束，若需結束則結束方法，若需繼續運行則返回讀取數據步驟，讀取新數據進行移位對齊和相關性比較。
【文檔編號】G01N1/30GK104266894SQ201410453967
【公開日】2015年1月7日 申請日期:2014年9月5日 優先權日:2014年9月5日
【發明者】張申, 張然, 程婷婷 申請人:中國礦業大學