一种基于相关性分析的矿山微震信号初至波时刻提取方法

摘要

一种基于相关性分析的矿山微震信号初至波时刻提取方法，属于矿山微震信号的分析处理方法。其功能为智能识别两信号时间差、实现自动移位对齐和判断矿山微震信号相关性，该方法主要分为五个步骤：一是读取矿山微震信号；二是求解两矿山微震信号的相关函数；三是根据求解出的相关函数，求出最大相关时的时间差；四是根据计算出的时间差对其中一个信号进行移位；五是对对齐后的两矿山微震信号求解其相关系数，判断其相关性。本发明较好的解决了智能识别矿山微震信号到达时间差及判断矿山微震信号是否来自同一震源的问题。

说明书

技术领域

本发明涉及矿山微震信号的分析处理方法，特别涉及于一种基于相关性分析的矿山微 震信号初至波时刻提取方法。

背景技术

在自然地震、工程爆破、爆炸等震动事件中，其中有一部分能量必然转化为震动波的 形式，它会以震源为中心向周围传播。矿山微震信号初至波拾取是震动研究中非常关键和 重要的问题。在震源定位中，准确快捷的拾取到初至波的时刻是进行准确震源定位或信号 分析的基础。对采集到的矿山微震信号进行分析和处理，并以此来确定矿山微震信号初至 波到达时刻的技术称为初至波拾取技术。初至波拾取技术在军事、民用和工业工程领域都 具有非常广泛的应用研究。

而矿震是采矿活动诱发的矿井岩体突然失稳破坏的动力现象，它严重威胁着矿井生产 与矿工的生命安全，对其进行实时监测和预警具有重要的理论和现实意义。矿震信号相对 于地震信号而言，具有震级小、震源浅和影响范围有限等特点，所以可以称之为矿山微震 信号。由于可以通过分析矿山微震信号的分布特点及发生频率大致标注出危险地带，实现 有效的监测与预警，且因矿山微震信号本身能量释放较小，震波传输距离有限，从而要求 矿山微震信号的定位精度相比普通地震信号更高。而定位精度的提高又需要初至波时刻更 为快速准确的提取。

目前一般采用的矿山微震信号初至波自动拾取技术并无专门针对矿山微震信号的特性 和要求的。初至波自动拾取目的是要确定其信号中纯噪声信号和有效信号之间的分界时 刻，通常都是根据矿山微震信号的振幅、频率和相位的变化来确定这一时刻。传统的初至 拾取方法主要分为两大类：一类是基于地震记录瞬时特征的方法，如极值法(峰值检测)、 差分法；这类方法对噪声比较敏感，当地震记录的噪声较严重时，难以准确拾取初至。另 一类方法是基于地震记录整体特征的方法，如相关法；这类方法虽然对噪声有较好的抑制 作用，但受到地震道之间相似性等因素的影响，对于复杂地震记录，初至拾取的精度也会 受到影响。

到目前为止，已经提出了许多初至拾取的方法，如人工拾取法、相关法、能量比法、 最大振幅法、分形维法及神经网络法等。

人工拾取法简单易行，但受人为因素和主观因素影响较大，容易引入人为误差，会直 接导致结果误差增大。

Gelchinsky和Shtivelman提出了一种用相邻道进行互相关的方法，它假设各道的脉 冲形状不发生变化。相关法受初至的续至影响较大，并且与子波的选择关系十分密切，同 时要求选择合适的时窗范围，这些在实际工程中都有一定的难度。

Hatherrly提出了线性最小平方预测技术与拐点校正相结合的方法，他提出首先识别 第一个峰值和拐点，然后估算二者的统计差值。黄成之等采用统计方法将地震初至波记录 分成信号和噪音两个部分，并使这两部分统计特征之间的差别为最大。统计特征法收到地 震波形相似性影响较大，会对精度造成一定影响。

能量比值法使用周期内的信号能量与总时窗能量的比值，对初至比较敏感，续至波衰 减比较快，所以将比值的最大值点作为初至的近似值并作适当的时移，即为初至时刻。 Coppens提出了在不同大小的时窗内进行能量比较的方法。江玉乐等提出同极性能量比值 法，即改进的能量比法。由于能量比值法的抗干扰能力还不够好，所以对于初至波形发生 明显变化的地区拾取的初至时刻不够准确。

时间域分形维方法拾取初至的过程必须插值，且结果强烈的依靠插值的准确性。Fabio  Boschetti等提出了一种基于分形维的初至检测算法，该方法是基于地震道随着信号出现其 分形维值发生变化的特征来确定地震道初至。但其对时窗和步长的选取十分敏感，稍有不 慎就会严重的影响其结果。

神经网络法利用多参数特征进行模式识别，充分利用地震记录的瞬时特征和整体特 征。由于神经网络不仅具有并行处理、自组织自学习能力，而且具有高度鲁棒性、容错性 和高度的映射、计算和分类能力。庄东海等采用将地震记录初至拾取看作一个模式识别过 程，充分利用地震记录的瞬时特征和整体特征，用人工神经网络方法进行地震记录初至拾 取，其能获得较好的实际效果。但其主要缺点是算法复杂性太高，搜索需要一定时间。

地震信号通过小波多分辨分解，可以有效地分离、消除噪声，有利于分形维与神经网 络法提高拾取初至的精度。罗光提出了改进型基于小波变换的初至波拾取方法，杨俊峰提 出了基于小波变换的三分向震相识别法、能量因子法。但有时其仍需要人工选取数据段来 进行拾取以减少拾取时间。

除此之外，有一些方法还依赖该道与其近道之间的对比，虽然这一类方法对噪声有一 定的压制作用，但受到地震道之间相似性等因素的影响，对于复杂地震记录，初至拾取的 精度也会受到影响。

综上所述，地震波信号是一个一维的时间序列，由于该序列仅仅表明了时间和幅度之 间的关系，且有噪声信号的干扰，对直接进行地震波初至的拾取造成了相当大的困难。解 决此类问题可采用如上文所述的智能算法，但由于智能算法普遍具有较高的算法复杂度， 不适合用于需要快速、并有较高定位精度的矿山微震环境，而本发明较好的解决了对初至 波时间的拾取这一问题，并具有较低的算法复杂度。

发明内容

本发明的目的是针对已有技术中存在的问题，提供一种基于相关性分析的矿山微震信 号初至波时刻提取方法，解决矿山微震信号处理中初至波时间智能识别、矿山微震信号对 齐和判断矿山微震信号是否来自同一震源的问题。

实现本发明目的的技术方案：本发明的矿山微震信号初至波时刻提取方法包括五个步 骤：一是读取矿山微震信号；二是求解两矿山微震信号的相关函数；三是根据求解出的相 关函数，求出最大相关时的时间差；四是根据计算出的时间差对其中一个信号进行移位； 五是对对齐后的两矿山微震信号求解其相关系数，判断其相关性；具体方法步骤如下：

(一)读取矿山微震信号：将两矿山微震信号读取进系统，定义x_n为一矿山微震信号 序列，y_n为另一矿山微震信号序列，要求两序列必须等长，设定其序列长度为L，且两信 号采样频率均为f_s；

(二)求解两矿山微震信号的相关函数：根据下式求解矿山微震信号相关函数，在相 关函数取得最大值时，根据采样点及采样率计算出矿山微震信号间的时间差； 根据离散信号序列的相关函数公式计算相关函数R_xy为：

$R_{\mathrm{xy}}\left(m\right)=\underset{n=-\infty }{\overset{\infty }{\Sigma }}x\left(n\right)y(n+m)$

其中变量m取值范围为0至L；

(三)根据求解出的相关函数，求出最大相关时的时间差：根据计算出的时间差对矿 山微震信号进行对齐；根据相关函数、采样点与采样率计算最大相关时的时间差，取得相 关函数最大值R_max为：

R_max＝max[|R_xy(m)|]

(四)根据计算出的时间差对其中一个信号进行移位：对对齐后的矿山微震信号进行 相关性分析，并设定阈值以判断矿山微震信号是否来自于同一震源； 取得初始时刻至相关函数绝对值最大时刻的采样点总数，记作N，则时间偏差T_offset：

T_offset＝N/f_s

(五)对对齐后的两矿山微震信号求解其相关系数，判断其相关性：

两矿山微震信号的相关系数矩阵R_corr为：

$R_{\mathrm{corr}}=\left(\begin{array}{cc}{\rho }_{11}& {\rho }_{12}\\ {\rho }_{21}& {\rho }_{22}\end{array}\right)$

其中ρ_ij＝E((X_i-E(X_i))·(Y_j-E(Y_j)))，其中变量i、变量j取值范围均为1至2，且X为x_n序列的总体，Y为y_n序列的总体，E为数学期望；求出的相关系数矩阵主对角线表示自相 关性，其始终为1；副对角线是两信号的相关性，绝对值越大，相关性越好，其取值范围 为[0，1]；分界值可根据实际情况不同设定为0到1之间的任何值；最终判断程序是否需 要结束，若需结束则结束方法，若需继续运行则返回读取数据步骤，读取新数据进行移位 对齐和相关性比较。

有益效果，由于采用了上述方案，智能识别两信号时间差、实现自动移位对齐和判断 矿山微震信号相关性，实现智能识别矿山微震信号到达时间差；实现矿山微震信号的智能 对齐，以便于存储、分析和处理；通过相关性比较两矿山微震信号是否来自于同一震源； 解决了智能识别矿山微震信号到达时间差及判断矿山微震信号是否来自同一震源的问题。

所述矿山微震信号初至波时刻提取方法由于采用相关性分析的方法，算法简便，相比 背景技术中的模糊识别方法、时间域分形维方法、神经网络方法、小波变换方法等方法具 有较小的程序复杂度，更适合应用于存在大量震动数据、对及时性要求较高场合的使用， 相比人工识别方法、能量比值法等方法减少了主观因素产生的误差。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的 示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

图1为本发明方法的流程图。

图2为本发明实施例的系统结构示意图。

图3为本发明实施例的a、b、c、d信号波形图。

图4为本发明方法实施例中a信号与b信号相关函数及移位对齐效果仿真图。

图5为本发明方法实施例中a信号与d信号相关函数及移位对齐效果仿真图。

具体实施方式

下面结合附图中实施例对本发明作进一步说明：

矿山微震信号初至波时刻提取方法包括五个步骤：一是读取矿山微震信号；二是求解 两矿山微震信号的相关函数；三是根据求解出的相关函数，求出最大相关时的时间差；四 是根据计算出的时间差对其中一个信号进行移位；五是对对齐后的两矿山微震信号求解其 相关系数，判断其相关性；具体方法步骤如下：

(一)读取矿山微震信号：将两矿山微震信号读取进系统，定义x_n为一矿山微震信号 序列，y_n为另一矿山微震信号序列，要求两序列必须等长，设定其序列长度为L，且两信 号采样频率均为f_s；

(二)求解两矿山微震信号的相关函数：根据下式求解矿山微震信号相关函数，在相 关函数取得最大值时，根据采样点及采样率计算出矿山微震信号间的时间差； 根据离散信号序列的相关函数公式计算相关函数R_xy为：

$R_{\mathrm{xy}}\left(m\right)=\underset{n=-\infty }{\overset{\infty }{\Sigma }}x\left(n\right)y(n+m)$

其中变量m取值范围为0至L；

(三)根据求解出的相关函数，求出最大相关时的时间差：根据计算出的时间差对矿 山微震信号进行对齐；根据相关函数、采样点与采样率计算最大相关时的时间差，取得相 关函数最大值R_max为：

R_max＝max[|R_xy(m)|]

(四)根据计算出的时间差对其中一个信号进行移位：对对齐后的矿山微震信号进行 相关性分析，并设定阈值以判断矿山微震信号是否来自于同一震源； 取得初始时刻至相关函数绝对值最大时刻的采样点总数，记作N，则时间偏差T_offset：

T_offset＝N/f_s

(五)对对齐后的两矿山微震信号求解其相关系数，判断其相关性：

两矿山微震信号的相关系数矩阵R_corr为：

$R_{\mathrm{corr}}=\left(\begin{array}{cc}{\rho }_{11}& {\rho }_{12}\\ {\rho }_{21}& {\rho }_{22}\end{array}\right)$

其中ρ_ij＝E((X_i-E(X_i))·(Y_j-E(Y_j)))，其中变量i、变量j取值范围均为1至2，且X 为x_n序列的总体，Y为y_n序列的总体，E为数学期望；求出的相关系数矩阵主对角线表示 自相关性，其始终为1；副对角线是两信号的相关性，绝对值越大，相关性越好，其取值 范围为[0，1]；分界值可根据实际情况不同设定为0到1之间的任何值。最终判断程序是 否需要结束，若需结束则结束方法，若需继续运行则返回读取数据步骤，读取新数据进行 移位对齐和相关性比较。

本发明利用相关函数对信号进行相关性分析，从而得到初至波时刻、进行信号自动对 齐和判断矿山微震信号是否来自同一震源。

附图1为本发明方法的流程图，所述基于相关性分析的矿山微震信号初至波时刻提取 方法，按照附图1的流程完成时刻提取、信号对齐及判断信号是否来自同一震源，其基本 步骤如下：

读取矿山微震信号x_n，矿山微震信号y_n，其采样频率均为f_s；

根据离散信号序列的相关函数公式计算相关函数：

$R_{\mathrm{xy}}\left(m\right)=\underset{n=-\infty }{\overset{\infty }{\Sigma }}x\left(n\right)y(n+m)---\left(1\right)$

根据相关函数、采样点与采样率计算最大相关时的时间差

取得相关函数最大值：

R_max＝max[|R_xy(m)|]   (2)

取得初始时刻至相关函数绝对值最大时刻的采样点总数，记作N

则时间偏差：

T_offset＝N/f_s   (3)

依据上一步骤计算出的T_offset将两信号移位对齐，以便下一步骤计算相关系数；

根据公式求两个信号的相关系数矩阵：

$R_{\mathrm{corr}}=\left(\begin{array}{cc}{\rho }_{11}& {\rho }_{12}\\ {\rho }_{21}& {\rho }_{22}\end{array}\right)---\left(4\right)$

其中ρ_ij＝E((X_i-E(X_i))·(Y_j-E(Y_j)))，且X为x_n序列的总体，Y为y_n序列的总体，E 为数学期望。求出的相关系数矩阵主对角线表示自相关性，其始终为1；副对角线是两信 号的相关性，绝对值越大，相关性越好。一般情况下两信号相关系数绝对值在0-0.3之间 为微相关；0.3-0.5之间为实相关；0.5-0.8之间为显著相关；0.8-1.0之间为高度相关。该 方法通过设定分界值来判定两信号是否具有高度相关性，从而判断其是否来自同一震源， 该分界值可以根据实际需要进行更改。

所述基于相关性分析的矿山微震信号初至波时刻提取方法，可应用于网络化分布式矿 山震源定位系统中，附图2为该系统结构示意图，下面结合附图2对所述方法的具体实施 方式做一详细说明。附图2中，圆点为传感器节点，包含微震传感器，可以监测矿山微震 信号。当震源发生震动时，临近的传感器节点可以采集到矿山微震数据。附图3为收到矿 山微震信号的四个不同震动传感器记录的波形，其中a、b、c三个信号来自于同一震源， d来自于不同于其它三个信号的另一震源。在本实施例中，根据大量样本及相关经验设定 判定是否来自同一震源的分界值为0.5。

首先对a与b信号进行相关性分析。附图4为本发明方法实施例中a信号与b信号相 关函数及移位对齐效果仿真图。据式(1)计算a与b信号的相关函数并作图，求得的相 关函数具体分布见附图4。据式(2)取得相关函数最大值R_max＝250.8038。据式(3)计算 a与b信号两信号时间差，T_offset＝-0.0393(s)。据式(4)计算两信号相关系数矩阵。在 相关系数矩阵中主对角线始终为1，副对角线即两信号相关性，绝对值越大则相关性越好。 求得的a与b信号相关系数矩阵为：

$R_{\mathrm{ab}}=\left(\begin{array}{cc}1.0000& -0.5462\\ -0.5462& 1.0000\end{array}\right)$

将其相关系数矩阵副对角线取绝对值，即得a信号与b信号相关性为0.5462，大于设 定的分界值0.4，则所述方法判定a信号与b信号来自不同震源。

其次对a与d信号进行相关性分析。附图5为本发明方法实施例中a信号与d信号相 关函数及移位对齐效果仿真图。据式(1)计算a与d信号的相关函数并作图，求得的相 关函数具体分布见附图5。据式(2)取得相关函数最大值R_max＝31.5554。据式(3)计算 a与b信号两信号时间差，T_offset＝0.0533(s)。据式(4)计算两信号相关系数矩阵。求得 的a与d信号相关系数矩阵为：

$R_{\mathrm{ad}}=\left(\begin{array}{cc}1.0000& -0.0566\\ -0.0566& 1.0000\end{array}\right)$

将其相关系数矩阵副对角线取绝对值，即得a信号与d信号相关性为0.0566，小于设 定的分界值0.4，则所述方法判定a信号与d信号来自不同震源。

以上结合附图对本发明的具体实施方式作了说明，但这些说明不能被理解为限制了本 发明的范围，本发明的保护范围由随附的权利要求书限定，任何在本发明权利要求基础上 的改动都是本发明的保护范围。