CN107505652B - 一种基于能量分布特征的矿山微震信号辨识方法 - Google Patents

Abstract

本发明公开了一种基于能量分布特征的矿山微震信号辨识方法，属于信号分析及识别领域，包括如下步骤：读取待辨识微震信号x(t)；对x(t)进行VMD分解，得到K个按照频率从高到低顺序排列的变分模态分量；计算出各模态分量的频带能量，提取各模态分量占原信号的能量百分比值构成能量分布向量P；以能量分布向量P为基础计算出能量分布X轴重心系数cx；根据辨识阈值T识别矿山微震信号，若cx>T为矿山煤岩体破裂微震信号，若cx≤T为爆破震动信号；最后对辨识阈值T的值进行自适应更新。本发明能有效区分煤岩体破裂微震信号和爆破震动信号，具有自适应性强、准确性高等特点。

Description

一种基于能量分布特征的矿山微震信号辨识方法

技术领域

本发明属于信号分析及识别领域，具体涉及一种基于能量分布特征的矿山微震信号辨识方法。

背景技术

微震监测是近年来发展起来的先进且行之有效的煤岩动力灾害监测手段，它能够对煤岩体微震活动实时、连续、在线监测，形成微震监测数据。由于矿山环境复杂，存在现场背景噪声、爆破震动等大量干扰信号，使得微震监测系统无法准确识别并记录有效微震事件，后期需要依靠技术人员人工识别出有效微震事件，严重影响了微震监测系统的识别效率。由于煤矿爆破作业经常发生，而煤岩体微震和爆破震动波形又极为相似，采用人工识别方式，经常出现误处理，识别难度大。

目前，针对矿山微震信号波形识别的常用时频分析法包括傅立叶变换、小波变换、小波包变换、频率切片小波变换和EMD等，传统傅立叶变换主要用于分析周期性平稳信号，对包含有尖峰和突变的随机性、非平稳性微震信号分析效果欠佳；小波分析能同时进行时频分析，但需要选择合适的小波基才能达到较好的分解效果；EMD能较好地处理随机非平稳信号，但EMD方法存在边界效应及模态混叠现象，导致EMD具有不稳定性和不唯一性。这些方法用于信号分析时均存在一定程度的弊端，为信号辨识增加了难度，误判率高。

发明内容

针对现有技术中存在的上述问题，本发明提出了一种基于能量分布特征的矿山微震信号辨识方法，采用变分模态分解(VMD)对信号进行时频分析。VMD是一种新的信号分解方法，相比于其它模态分解技术，它具有坚实的理论基础，消除了模态混叠问题，克服了现有技术的不足，具有良好的频域自适应分解效果。

为了实现上述目的，本发明采用如下技术方案：

一种基于能量分布特征的矿山微震信号辨识方法，包括如下步骤：

步骤1：读取待辨识微震信号x(t)，其中，t＝1，2，…，N，N为微震信号的采样点个数；

步骤2：对待辨识微震信号x(t)进行VMD分解，得到K个按照频率从高到低顺序排列的变分模态分量{u₁,…u_k,…,u_K}：

对待辨识微震信号x(t)采用VMD分解为K个变分模态分量，约束条件为使各个模态的估计带宽之和最小，且各模态之和等于待辨识微震信号x(t)，约束变分模型描述为式(1)和式(2)：

s.t.∑ku_k＝x(t) (2)；

其中，x(t)表示待辨识微震信号，{u_k}:＝{u₁,…,u_K}代表分解得到的K个有限带宽的变分模态分量，{ω_k}:＝{ω₁,…,ω_K}表示各分量的频率中心，δ(t)为狄拉克(Dirac)函数，*表示卷积，j²＝-1；式(2)中，为对所有的变分模态求和；

为求解式(1)和式(2)的最优解，引入扩展的Lagrange将约束变分问题变为非约束变分问题，其表达式为式(3)：

其中，α为二次惩罚因子，λ(t)为拉格朗日乘法算子；求解该变分问题的具体步骤如下：

步骤2.1：定义变分模态分量个数K值与惩罚因子α的值；

步骤2.2：初始化n＝0；

步骤2.3：令n＝n+1，执行整个循环；

步骤2.4：执行内层第一个循环，根据式(4)更新u_k：

其中，为待辨识微震信号x(t)的傅立叶变换，

步骤2.5：令k＝k+1，重复步骤2.4，直到k＝K，结束内层第一个循环；

步骤2.6：执行内层第二个循环，根据式(5)更新ω_k：

步骤2.7：令k＝k+1，重复步骤2.6，直到k＝K，结束内层第二个循环；

步骤2.8：执行外层循环，根据式(6)更新λ：

其中，τ为拉格朗日乘法算子λ(t)的更新步长参数；

步骤2.9：重复步骤2.3至步骤2.8，直到满足迭代停止条件如式(7)所示，结束整个循环，输出结果，得到K个变分模态分量：

其中，ε为求解精度；

步骤3：计算待辨识微震信号x(t)的能量分布向量P：

根据公式(8)计算各模态分量u_k对应的能量E_k：

其中，x_ik(t)(i＝1,2,…N；k＝1,2,…,K；N为采样点个数，K为变分模态个数)表示模态分量u_k时序序列的离散点幅值；

根据每个模态分量u_k的能量以及待辨识微震信号x(t)的总能量，可以得到模态分量u_k的能量百分比值从而得到该微震信号的能量分布向量P，即P＝[P(1),…,P(k),…,P(K)]；

步骤4：计算待辨识微震信号x(t)的能量分布X轴的重心系数cx：

根据公式(9)计算能量分布X轴重心系数cx：

步骤5：根据辨识阈值T识别待辨识微震信号x(t)，若cx>T为矿山煤岩体破裂微震信号，cx≤T为爆破震动信号；

步骤6：自适应更新辨识阈值T的值：

根据方程组(10)更新辨识阈值T：

其中，W₁为训练集中煤岩体破裂微震信号的cx值集合，W₂为训练集中爆破震动信号的cx值集合。

本发明原理如下：

为实现煤岩体破裂微震信号和爆破震动信号的有效分类辨识，本发明利用两种微震信号能量分布差异显著的特点，首先读取待辨识微震信号并进行VMD分解，得到K个按照频率从高到低顺序排列的变分模态分量；计算出各模态分量的频带能量，提取各模态分量占原信号的能量百分比值构成能量分布向量P；以能量分布向量P为基础计算出能量分布X轴重心系数cx；根据辨识阈值T识别矿山微震信号，若cx>T时，检测微震信号为矿山煤岩体破裂微震信号，若cx≤T时，检测微震信号为爆破震动信号。该方法可以实现对煤岩体破裂微震信号和爆破震动信号的有效辨识。

本发明采用以上技术方案，与现有技术现比，具有以下优点：

本发明依据VMD良好频谱分解特征对矿山微震信号进行自适合剖分，依据煤岩体破裂微震信号和爆破震动信号在不同频段上能量分布的显著差异，通过计算微震信号能量分布X轴重心系数，实现对两类矿山微震信号的有效辨识，该方法具有算法简单、自适应性和实时性强的特点，具有很好的技术价值和应用前景。

附图说明

图1为本发明一种基于能量分布特征的矿山微震信号辨识方法的流程图。

图2为待辨识微震信号x(t)的示意图及其时频图。

图3为待辨识微震信号x(t)经VMD分解后得到的6个变分模态分量示意图及其时频图。

图4为待辨识微震信号x(t)的能量分布直方图。

图5为15组煤岩体破裂微震测试信号的能量向量、重心系数及其辨识结果图。

图6为15组爆破震动测试信号的能量向量、重心系数及其辨识结果图。

图7为测试组微震信号分类识别结果。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

一种基于能量分布特征的矿山微震信号辨识方法，其流程如图1所示，具体包括如下步骤：

步骤1：读取待辨识微震信号x(t)，其中，t＝1，2，…，N，N为微震信号的采样点个数；

步骤2：对待辨识微震信号x(t)进行VMD分解，得到K个按照频率从高到低顺序排列的变分模态分量{u₁,…u_k,…,u_K}：

对待辨识微震信号x(t)采用VMD分解为K个变分模态分量，约束条件为使各个模态的估计带宽之和最小，且各模态之和等于待辨识微震信号x(t)，约束变分模型描述为式(1)和式(2)：

s.t.∑_ku_k＝x(t)(2)；

其中，x(t)表示待辨识微震信号，{u_k}:＝{u₁,…,u_K}代表分解得到的K个有限带宽的变分模态分量，{ω_k}:＝{ω₁,…,ω_K}表示各分量的频率中心，δ(t)为狄拉克(Dirac)函数，*表示卷积，j²＝-1；式(2)中，为对所有的变分模态求和；

为求解式(1)和式(2)的最优解，引入扩展的Lagrange将约束变分问题变为非约束变分问题，其表达式为式(3)：

其中，α为二次惩罚因子，λ(t)为拉格朗日乘法算子；求解该变分问题的具体步骤如下：

步骤2.1：定义变分模态分量个数K值与惩罚因子α的值；

步骤2.2：初始化n＝0；

步骤2.3：令n＝n+1，执行整个循环；

步骤2.4：执行内层第一个循环，根据式(4)更新u_k：

其中，为微震信号x的时序序列x(t)的傅立叶变换，

步骤2.5：令k＝k+1，重复步骤2.4，直到k＝K，结束内层第一个循环；

步骤2.6：执行内层第二个循环，根据式(5)更新ω_k：

步骤2.7：令k＝k+1，重复步骤2.6，直到k＝K，结束内层第二个循环；

步骤2.8：执行外层循环，根据式(6)更新λ：

其中，τ为拉格朗日乘法算子λ(t)的更新步长参数；

步骤2.9：重复步骤2.3至步骤2.8，直到满足迭代停止条件如式(7)所示，结束整个循环，输出结果，得到K个变分模态分量：

其中，ε为求解精度；

步骤3：计算待辨识微震信号x(t)的能量分布向量P：

根据公式(8)计算各模态分量u_k对应的能量E_k：

其中，x_ik(t)(i＝1,2,…N；k＝1,2,…,K；N为采样点个数，K为变分模态个数)表示模态分量u_k时序序列的离散点幅值；

根据每个模态分量u_k的能量以及待辨识微震信号x(t)的总能量，可以得到模态分量u_k的能量百分比值从而得到该微震信号的能量分布向量P，即P＝[P(1),…,P(k),…,P(K)]；

步骤4：计算待辨识微震信号x(t)的能量分布X轴的重心系数cx：

根据公式(9)计算能量分布X轴重心系数cx：

步骤5：根据辨识阈值T识别待辨识微震信号x(t)，若cx>T为矿山煤岩体破裂微震信号，cx≤T为爆破震动信号；

步骤6：自适应更新辨识阈值T的值：

根据方程组(10)更新辨识阈值T：

其中，W₁为训练集中煤岩体破裂微震信号的cx值集合，W₂为训练集中爆破震动信号的cx值集合。

如图2所示，步骤1获取以时间(s)为横轴，振幅为纵轴，采样频率fs＝1000Hz的待辨识微震信号x(t)，t＝1*1/fs,2*1/fs,…,5000*1/fs，微震信号采样点数据见表1。

表1监测信号采样点数据(可以存储于Excel表中)

序号	采样点(N)	振幅
			1*1/fs	1	4.34E-08
2*1/fs	2	1.69E-07
			3*1/fs	3	1.41E-07
4*1/fs	4	-2.43E-07
			5*1/fs	5	-6.50E-07
…	…
			4999*1/fs	4999	-6.10E-06
5000*1/fs	5000	-6.25E-06

按照步骤2的VMD算法对待辨识微震信号x(t)进行变分模态分解，取K＝6，二次惩罚因子α＝2000，分解后的6个变分模态分量及其时频谱如图3所示。

按照步骤3的方法，计算出各模态分量的频带能量，提取各模态分量占原信号的能量百分比值，得到该微震信号的能量分布向量P，即P＝[0.02，0.09，3.01，1.45，13.41，82.02]。

按照步骤4的方法，应用公式(9)计算待辨识微震信号x(t)的能量分布X轴重心系数cx，得到cx＝0.957，图4为该微震信号的能量分布直方图，图中黑色实心圆为该微震信号的能量分布平面重心位置。

按照步骤5的方法，根据辨识阈值T＝0.56以及待检测微震信号x(t)的cx＝0.957对该待测微震信号进行辨识，因为cx>T，所以，该待检测微震信号x(t)辨识为煤岩体破裂微震信号。

测试组中分别给出了15组煤岩体破裂微震信号和15组爆破震动微震信号，煤岩体破裂微震测试信号的能量向量、重心系数及其辨识结果如图5所示；爆破震动微震测试信号的能量向量、重心系数及其辨识结果如图6所示。根据检测结果，该测试组共30组微震信号，其中29组辨识正确，1组辨识错误，辨识正确率为96.67％。测试组微震信号分类识别结果如图7所示。

按照步骤6的方法，将测试组数据加入已有训练组中，对辨识阈值自适应更新为T＝0.61，以继续提高辨识正确率。

微震信号是非平稳随机信号，其频率分布较为分散，不同种类的微震信号在不同频带能量分布差异显著，因此可根据这一特点，提取微震信号的能量分布向量，通过计算能量分布重心系数并与辨识阈值进行对比，即可实现对待测微震信号的分类辨识。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (1)

1.一种基于能量分布特征的矿山微震信号辨识方法，其特征在于：包括如下步骤：

步骤1：读取待辨识微震信号x(t)，其中，t＝1，2，…，N，N为微震信号的采样点个数；

步骤2：对待辨识微震信号x(t)进行VMD分解，得到K个按照频率从高到低顺序排列的变分模态分量{u₁,…u_k,…,u_K}：

对待辨识微震信号x(t)采用VMD分解为K个变分模态分量，约束条件为使各个模态的估计带宽之和最小，且各模态之和等于待辨识微震信号x(t)，约束变分模型描述为式(1)和式(2)：

s.t.∑_ku_k＝x(t) (2)；

其中，x(t)表示待辨识微震信号，{u_k}:＝{u₁,…,u_K}代表分解得到的K个有限带宽的变分模态分量，{ω_k}:＝{ω₁,…,ω_K}表示各分量的频率中心，δ(t)为狄拉克(Dirac)函数，*表示卷积，j²＝-1；式(2)中，为对所有的变分模态求和；

为求解式(1)和式(2)的最优解，引入扩展的Lagrange将约束变分问题变为非约束变分问题，其表达式为式(3)：

其中，α为二次惩罚因子，λ(t)为拉格朗日乘法算子；

求解该变分问题的具体步骤如下：

步骤2.1：定义变分模态分量个数K值与惩罚因子α的值；

步骤2.2：初始化n＝0；

步骤2.3：令n＝n+1，执行整个循环；

步骤2.4：执行内层第一个循环，根据式(4)更新u_k：

其中，为待辨识微震信号x(t)的傅立叶变换，

步骤2.5：令k＝k+1，重复步骤2.4，直到k＝K，结束内层第一个循环；

步骤2.6：执行内层第二个循环，根据式(5)更新ω_k：

步骤2.7：令k＝k+1，重复步骤2.6，直到k＝K，结束内层第二个循环；

步骤2.8：执行外层循环，根据式(6)更新λ：

其中，τ为拉格朗日乘法算子λ(t)的更新步长参数；

步骤2.9：重复步骤2.3至步骤2.8，直到满足迭代停止条件如式(7)所示，结束整个循环，输出结果，得到K个变分模态分量：

其中，ε为求解精度；

步骤3：计算待辨识微震信号x(t)的能量分布向量P：

根据公式(8)计算各模态分量u_k对应的能量E_k：

其中，x_ik(t)(i＝1,2,…N；k＝1,2,…,K；N为采样点个数，K为变分模态个数)表示模态分量u_k时序序列的离散点幅值；

根据每个模态分量u_k的能量以及待辨识微震信号x(t)的总能量，可以得到模态分量u_k的能量百分比值从而得到该微震信号的能量分布向量P，即P＝[P(1),…,P(k),…,P(K)]；

步骤4：计算待辨识微震信号x(t)的能量分布X轴的重心系数cx：

根据公式(9)计算能量分布X轴重心系数cx：

步骤5：根据辨识阈值T识别待辨识微震信号x(t)，若cx>T为矿山煤岩体破裂微震信号，cx≤T为爆破震动信号；

步骤6：自适应更新辨识阈值T的值：

根据方程组(10)更新辨识阈值T：

其中，W₁为训练集中煤岩体破裂微震信号的cx值集合，W₂为训练集中爆破震动信号的cx值集合。