一种基于相关性分析的矿山微震信号初至波时刻提取方法
【专利摘要】一种基于相关性分析的矿山微震信号初至波时刻提取方法，属于矿山微震信号的分析处理方法。其功能为智能识别两信号时间差、实现自动移位对齐和判断矿山微震信号相关性，该方法主要分为五个步骤：一是读取矿山微震信号；二是求解两矿山微震信号的相关函数；三是根据求解出的相关函数，求出最大相关时的时间差；四是根据计算出的时间差对其中一个信号进行移位；五是对对齐后的两矿山微震信号求解其相关系数，判断其相关性。本发明较好的解决了智能识别矿山微震信号到达时间差及判断矿山微震信号是否来自同一震源的问题。
【专利说明】—种基于相关性分析的矿山微震信号初至波时刻提取方法
【技术领域】
[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,0()()0
[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日
【发明者】张申, 张然, 程婷婷 申请人:中国矿业大学