一种矿用微震监测信号的处理方法与流程

本发明涉及矿山微震监测技术领域，尤其涉及一种矿用微震监测信号的处理方法。

背景技术：

在矿山开采中，人为开采挖掘、爆破以及车辆机械振动会破坏矿体内原有应力平衡状态，导致围岩变形、移动和破坏，引起岩层塌陷，严重威胁着矿山井下工作人员的安全。因此，地压监测在深部开采中占有极其重要的地位，矿山岩体在变形破坏的整个过程中几乎都伴随着裂纹的产生、扩展、摩擦，积聚的能量在释放的过程中会产生量级较小的震动，即微震事件，微震监测是地压监测中的一种重要的监测手段，它是通过对物体震动进行监测，得到微震事件发生的时间、位置及震级，从而判断岩体承受压力大小、受到破坏程度等情况，来对其安全程度进行评价的方法，为矿山进行安全预警、控制地压灾害提供了有效依据。其中，由于微震事件震级一般都很微小，而现场由于施工、爆破等产生持续而较大的震动事件，使采集环境信噪比降低，对采集系统对微震事件波形进行采集非常困难。因此，滤波处理是微震监测系统的重要模块。在微震系统中，波形质量的好坏直接影响着微震数据处理结果，包括震源定位、震级等信息，而这些震源信息需要为进一步的安全风险评估提供可靠依据。

矿用微震监测信号滤波方法是矿用微震监测系统中的一个重要环节，数据质量好坏(数据信噪比的高低)直接影响着微震波形的分析结果。现有技术中常用的带通滤波、高通滤波、低通滤波等基础滤波方法，都是基于频域来实现不同频率噪声的处理，这些方法实现思想均是先对信号进行傅里叶变换，根据需要选取一个窗函数来圈出我们所需要的频率范围，接着在频谱上仅保留频域范围内的所有信号，从而滤除频带以外的所有信号，再用反傅里叶变换方法恢复成时域信号，即可得到指定频带的滤波信号。但上述类滤波方法对具有与微震信号不同频域的噪声抑制十分有效，但频域相似的噪声仍会产生较大干扰，使得微震信号质量不高；且该类滤波无法自动搜索有效信号所在的实际频率范围，需要人为参与指定频带来进行滤波，因此并没有完全过滤不同频带的噪声；而且也无法有效分析有效震动信号实际发生的时间，对于频率相似而发生时间不同的噪声也难以剔除，所以得到的滤波信号中仍然掺杂许多复杂未知噪声。

技术实现要素：

本发明的目的是提供一种矿用微震监测信号的处理方法，该方法能够分析得到有效信号实际发生的时间段，并可根据其时间段更精细地处理信号发生前后的噪声，提高信号整体信噪比。

本发明的目的是通过以下技术方案实现的：

一种矿用微震监测信号的处理方法，所述方法包括：

步骤1、从微震监测数据采集系统中接收现场采集的原始信号，通过信号升维法绘制原始信号在时间-频率二维空间中的能量分布图，得到清晰的时间-频率-能量谱；

步骤2、基于步骤1所得到的时间-频率-能量谱，利用无监督式聚类算法搜索并提取出有效震动信号所在的频带和时间段；

步骤3、基于所述时间-频率-能量谱，以及步骤2中提取出的有效震动信号所在的频带和时间段，建立基于全局噪声能量的滤波阈值，以及k个基于有效震动信号频带噪声能量的滤波阈值；

步骤4、根据步骤3中得到的基于全局噪声能量的滤波阈值，以及k个基于有效震动信号频带噪声能量的滤波阈值，从时域上完成信号滤波处理，得到滤波处理后的有效震动信号。

由上述本发明提供的技术方案可以看出，上述方法能够分析得到有效信号实际发生的时间段，并可根据其时间段更精细地处理信号发生前后的噪声，提高信号整体信噪比。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的矿用微震监测信号的处理方法流程示意图；

图2为按本发明实施例所述方法进行信号处理后的波形对比示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

下面将结合附图对本发明实施例作进一步地详细描述，如图1所示为本发明实施例提供的矿用微震监测信号的处理方法流程示意图，所述方法包括：

步骤1、从微震监测数据采集系统中接收现场采集的原始信号，通过信号升维法绘制原始信号在时间-频率二维空间中的能量分布图，得到清晰的时间-频率-能量谱；

在该步骤中，原始信号是指由传感器接收到的带有未知噪声和有效震动信号的混合信号。通过信号升维分析方法绘制原始信号在时间-频率二维空间中的能量分布图的方法具体如下：

采用小波变换将采集到的带有噪声的原始信号升维到时间-频率-能量的三维空间中，再利用其瞬时频率对频带并进行能量压缩，得到清晰的时间-频率-能量谱，其过程具体为：

采集的原始信号f(t)表示为：

式中，ak(t)为第k个谐波分量的瞬时振幅；θk(t)为第k个谐波分量的瞬时相位；e(t)为噪声或误差；k为可分解分量的个数；

对f(t)进行连续小波变换(cwt)，得到的小波系数wf(a,b)，表示为：

式中，a为尺度因子；b为平移因子；ψ^*为共轭小波函数；

其在频率域等价变换为：利用小波系数wf(a,b)得到的相位来计算其瞬时频率ωf(a,b)：

利用计算得到的瞬时频率，可以建立(a,b)→(ωf(a,b),b)之间的映射关系，经过同步压缩变换对时间—尺度平面的能量进行重新分配，将其转化为时间—频率平面。离散情况下，尺度坐标为(δa)i＝ai-ai-1，频率坐标为δω＝ωl-ωl-1，tf(ωl,b)在ωl的中心位置，其中ωl∈[ωl-δω/2,ωl+δω/2]。压缩变换tf(ωl,b)公式为：

步骤2、基于步骤1所得到的时间-频率-能量谱，利用无监督式聚类算法搜索并提取出有效震动信号所在的频带和时间段；

从时间-频率-能量谱中可以看出，有效震动信号的特征为能量分布集中，能量值较大，并且该信号仅发生在一定的时间段内；噪声信号的特征多为能量分布分散、能量值普遍较小，部分噪声为连续信号，遍布整个时域。

本实施例中采用k-means聚类算法对震动信号出现部分进行无监督式聚类学习，能够精准地找出有效震动信号所在的频带和时间段，其原理为：

首先，根据时间-频率-能量谱中信号能量的密集程度找出震动信号能量分布中心，确定中心个数k，即为该原始信号中包含k个震动信号；其次，对时间-频率-能量谱分别从时域的角度和频域的角度对能量进行叠加，得到其时域能量分布图和频域能量分布图，基于能量图，获取震动信号的k个初始聚类中心点的时间和频率范围；接着，采用k-means聚类算法，依据最近邻原则，将离各聚类中心最近的点分别划分到k个类中，并排除不符合时间和频率范围规定的点，计算平均值得到新的k个聚类中心，如此循环往复，依次迭代，使得评价准则函数j达到最优时，迭代结束得到k个震动信号能量分布区域；最后，分别获取这k个震动信号能量分布区域的时间和频率范围。

基于上述原理，本实例中在微震信号频段10～600hz的时间-频率-能量谱中集中搜索微震信号所在区域，以提取出所述有效震动信号所在的频带和时间段的具体过程为：

基于所得到的时间-频率-能量谱，首先基于时域得到该信号的时域能量分布图，能量集中的地方会形成k个波峰，即k个信号能量集群c＝{c1,c2,...,ck}；

结合微震监测数据采集系统，这里根据数据采集系统中对有效震动信号的触发及去触发设置即可知道一条原始数据中有效震动信号的最大个数kmax，则k≤kmax。这里可采用sta/lta长短时窗方法来进行触发和去触发操作，找到波峰个数k及各个波峰的宽度，即可得到k个聚类中心的大致时间范围；

接着基于所得到的k个聚类的大致时间范围，分别获得其各自时间宽度上的频域能量曲线，找出波峰位置，确定这k个聚类中心各自的频率范围；

则这k个聚类中心的时间及频率约束条件为：

其中，x代表时间，y代表频率，λ代表第λ个集群；这里，当波形宽度小于0.05秒，则认为是突发性噪声引起的振动，不是有效震动信号；上述x是数据采样点对应的时间，但不是数据所对应的时间，例如1秒内采样点有10个，则第10个采样点代表1秒；

接着计算能量值大于噪声水平的数据点到初始聚类中心的距离dij,λ：

再依据最近邻原则将离各个聚类中心较近的点划分到这k个聚类中心，找出有效震动信号发生的区域，即确定所述有效震动信号所在的频带和时间段。

步骤3、基于所述时间-频率-能量谱，以及步骤2中提取出的有效震动信号所在的频带和时间段，建立基于全局噪声能量的滤波阈值，以及k个基于有效震动信号频带噪声能量的滤波阈值。

在该步骤中，首先将全频带信号进行能量叠加，形成全局时域能量曲线s0，对s0中的点按能量大小重新排序，分为噪声区域、尾波区域和信号区域，这里针对噪声区域部分进行噪声水平评估；其中，该重排序列的第一个突变点之前部分即为背景噪声数据，通过对背景噪声数据分析即可测量出背景噪声水平，其具体操作为：

先根据现场大量数据得到的经验值来设定环境噪声最低信号阈值，对重排序列中最低信号阈值以下的部分，从重排序列第二个开始逐个计算两点之间的差值，对该差值序列进行排列，差值大于之前差值的平均值2倍，则可确定该点为突变点ny，则该点之前的部分为噪声信号区域，该原始信号全局噪声能量水平naverage为：

naverage＝(value(loc(ny)*20％)-value(loc(ny)*10％))*9

其中，naverage即为全局噪声能量水平，loc代表该点在能量曲线中的位置，value代表该点所在的能量值；

滤波时，基于全局噪声能量的滤波阈值即为该全局噪声能量水平naverage；

接着针对步骤2中所得到的k个有效震动信号，分别对其所在的频带基于时域进行压缩，分别得到时域能量曲线si,i＝1,2,...k，再按照上述操作步骤得到其各自的噪声能量阈值ni,i＝1,2,...,k，即k个基于有效震动信号频带噪声能量的滤波阈值。

步骤4、根据步骤3中得到的基于全局噪声能量的滤波阈值，以及k个基于有效震动信号频带噪声能量的滤波阈值，从时域上完成信号滤波处理，得到滤波处理后的有效震动信号。

在该步骤中，所述从时域上完成信号滤波处理，得到滤波处理后的有效震动信号的过程具体包括：

首先利用基于全局噪声能量的滤波阈值naverage对整个波形进行压制，从而剔除整体背景噪声干扰，得到：

其中，s0为全局时域能量曲线；

再利用得到的k个基于有效震动信号频带噪声能量的滤波阈值ni,i＝1,2,...,k，逐个对波形进行噪声压制：

对全局滤波信号f1(t)与k个滤波阈值ni,i＝1,2,...,k逐个进行对比，对能量小于ni的信号，对能量大于ni的信号保留，进行k次滤波后，最终得到滤波处理后的有效震动信号fk+1(t)为：

其中，si,i＝1,2,...k是k个频段的时域能量曲线；

如图2所示为按本发明实施例所述方法进行信号处理后的波形对比示意图。

值得注意的是，本发明实施例中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

综上所述，本发明实施例所述方法具有如下优点：

1)能够将时域信号分解到三维空间中，得到其时间-频率-能量谱，能够清晰地分辨出信号在时频域内的能量分布，从而辨别出有效信号和噪声的区别，进而解决频域相似噪声对有效信号的影响，提高信噪比；

2)能够自动搜索有效信号所在的实际频率范围，解决预先设定有效信号频带不准确导致的噪声干扰；

3)能够分析得到有效信号实际发生的时间段，并可根据其时间段更精细地处理信号发生前后的噪声，提高信号整体信噪比；

4)用时域滤波方法取代波形重构方法，解决了因波形重构引起的运算速度慢、信号丢失等问题，并极大程度地保留原始信号特征，快速有效地提高信号信噪比。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。