CN107992802B - 一种基于nmf的微震弱信号识别方法 - Google Patents

Abstract

本发明公开了一种基于NMF的微震弱信号识别方法，首先采用S变换对微震信号进行时频分析，然后对时频谱在频率方向上进行重排。采用非负矩阵分解（NMF）对重排的时频矩阵分解得到频域基向量和时域位置向量，从中提取尖锐度、导数平方和、信息熵以及稀疏度等特征参量，构造微震信号的特征空间，最后采用最小二乘支持向量机（LSSVM）对其分类。该方法增强了低频弱信号，也提高了时频分辨率，具有很好的时域和频域局部化能力。

Description

一种基于NMF的微震弱信号识别方法

技术领域

本发明涉及信号处理技术领域，具体是一种基于NMF的微震弱信号识别方法。

背景技术

微震信号的监测对煤矿安全生产具有重要的意义。对煤矿坍塌预警越早，造成事故的几率就越小。煤矿塌方初期会伴随着一系列的微弱岩体破裂信号的，但这些微弱岩裂信号很容易被机械噪声所淹没，因此寻找到一种能提取到微弱信号特征的方式是保障煤矿安全生产的关键。微震信号属于典型的非线性、低信噪比信号，其包含岩石破裂微震信号(下文简称微震信号)、煤矿爆破微震信号(下文简称爆破信号)、机械震动和其他施工噪声微震信号。

目前针对微震信号的分析主要采用短时傅里叶变换(STFT)、连续小波变换(CWT)时频分析方法。其中短时傅里叶变换(STFT)是Gabor在1946年提出的，在频域具有较好的局部化能力，但在时域没有局部化能力，并且对时变的非平稳信号，很难找到一个合适的时间窗口来适应不同的时间段。20世纪80年代，Y.Meyer等人引入了连续小波变换(CWT)很好的弥补了快速傅里叶变换的不足，在时域和频域均有良好的局部化能力，但在同一个工程问题用不同的小波函数其分析结果相差甚远。1996年，Stockwell提出了S变换(ST)集成了短时傅里叶变换(STFT)和连续小波变换(CWT)的优点，但是受海森堡不确定性原理的制约，其分辨率有限。

1998年，Huang提出经验模态分解(EMD)分解，2012年，程军圣在Huang的基础上提出了局部特征尺度分解法(LCD)提高了EMD分解的效率，但两者均是通过提取一系列不同特征时间尺度的本征模态分量(IMF)，存在频率混淆、端点效应，很难全面的提取微震信号的特征信息。

发明内容

本发明的目的在于克服现有技术的不足，而提供一种基于NMF的微震弱信号识别方法，该方法可以增强低频弱信号，也可以提高时频分辨率，具有很好的时域和频域局部化能力。

实现本发明目的的技术方案是：

一种基于NMF的微震弱信号识别方法，具体包括如下步骤：

1)对于采集的任意时变信号x(t)进行s变换，具体如下：

上述公式(1)中，ST_x(b，f)为信号x(t)时频谱，f为采样频率，t为采样时刻，b为时间位移参数，i为虚数单位；

2)由于受到海森堡不确定性原则的限制，其时频分辨率有限。时频谱分布在f₀附近的一定范围内，有一定的虚假频带宽度，因此在S变换的基础上，引入时频谱的重排，以提高弱振幅信号的时频分辨率，重排时频谱如下：

上述公式(2)中

是信号x(t)的重排谱矩阵，f_k是S变换的离散频率，

为重排重心，

为重排带宽；重排的意义：重排得到的谱图，它在任和点

处的频谱值等于该点附近谱图值的和，达到提高时频分辨率的目的；

3)对步骤2)得到的重排谱矩阵

做非负矩阵分解，分解得到频域基向量W_N×r和时域位置向量H_r×M；

4)对步骤3)得到的频域基向量W_N×r和时域位置向量H_r×M，提取尖锐度、导数平方和、信息熵和稀疏度构造微震信息的特征空间；

所述的尖锐度SH，能够反应信号的频域能量信息如下：

上述公式(3)中，

为每一个频域基向量{w_i}_{i＝1，...，r}的傅里叶变换(FT)，而

所述的导数平方和SD，能够直观地对向量中元素的均匀性进行度量，频域基向量{w_i}_{i＝1，...，r}的导数平方和定义如下：

上述公式(4)中，w_i′(n)为{w_i}_{i＝1，...，r}中相邻元素之差即：

w_i′(n)＝w_i(n+1)-w_i(n)n＝1，...，N-1；

所述的信息熵EN，反映序列的平均信息量，分别采用下式计算频域基向量{w_i}_{i＝1，...，r}和时域位置向量{h_i}_{i＝1，...，r}的信息熵EN：

稀疏度SP，反映元素之间偏移量，对{h_i}_{i＝1，...，r}计算稀疏度如下：

上述公式(7)中，M为h_i的长度。当向量中所有元素相等时，其稀疏度SP＝1；

5)从每个w_i提取(SH_wi，SD_wi，EN_wi)三个特征向量，每个h_i提取(EN_hi，SP_hi)两个特征参量，随着i取值的不同，从微震信号的时频谱矩阵中提取的特征集F可以表示如下：

F＝(SH_wi，SD_wi，EN_wi，EN_hi，SP_hi；...) (8)

6)对步骤5)的微震信号进行分类，以公式(8)构造的特征集为输入，输出微震信号的标签，其中岩石破裂信号的标签为1，爆破信号的标签为2，机械噪声的标签为3；

7)对微震信号x(t)进行分类，分类器的最优分类面表达式如下：

s.t y_i(wφ(x_i)+b)≥1-ζ_i

上述公式(9)中，w和b分别为最优分类面的权向量和偏差；ζ_i为松弛变量，ζ_i≥0；C为惩罚系数；

φ(x_i)由式(9)可得到：

上述公式(10)中，r为控制高斯核宽度的参数，采用网格搜索算法确定联合参数c和r；

8)分别选取微震信号特征集F的1-100组数据作为模式识别的训练组，101-200组数据作为预测组，经分类器输出训练组微震信号的类别标签；

经过上述步骤，完成微震信号的分类，正确率达到98％。

步骤3)中，具体步骤如下：

3.1)对

做变量代换，具体如下：

3.2)求V^T _N×M的分解矩阵W_N×r、H_r×M，具体步骤如下：

式中，W_N×r为基向量矩阵，H_r×M为系数矩阵。

3.3)采用乘法迭代法不断更新W和H矩阵的元素，迭代式(13)：

上述公式(13)中，a表示当前迭代次数，其中i≤m，j≤n，a≤min(m，n)。

有益效果：(1)该算法在时频分析时能有效分辨出弱振幅信号。能很好的应用于微震预警，而不是传统的用来分析微震事故的参考。(2)该算法能全面提取微震信号的潜在信息，使微震信号特征集更全面，从而使分类结果更精准。(3)该算法运算简单、便捷实时性更高。

附图说明

图1为本发明的一种基于NMF的微震弱信号的辨别方法流程图；

图2为信号波形图，3种典型微震信号波形，(a)为岩石破裂信号波形，(b)为爆破信号波形(c)表示机械噪声波形；

图3为S变换和重排S变换时频图，(a)为S变换时频图，(b)为重排ST时频图；

图4为取35组信号进行仿真分析，实际分类和预测分类的结果如图所示，微震信号的标签为1，爆破信号的标签为2，非威震信号的标签为3。

具体实施方式

下面结合附图和实施例对本发明做进一步阐述，但不是对本发明的限定。

实施例：

一种基于NMF的微震弱信号的辨别方法，流程图如图1所示，具体包括如下步骤：

(1)在某煤矿采集的微震数据是x＝[x₁，x₂，...，x_n]^T。

(2)对采集的数据，进行等间隔采样x(n)＝x(kn)，用matlab绘制了如图2所示其波形图。

(3)对信号x(n)进行重排S变换。

(4)如果x的列大于行，则进入步骤(5)，否则进入步骤(6)。

(5)对x做转置，并赋给A本身x＝x^T；保证x是一个列向量。

(6)如果输入参数的个数为1，进入步骤(7)，否则进入步骤(8)。

(7)令最小采样频率为f_min＝0，最大采样频率f_max＝f_Nyquist，采样间隔T＝1。

(8)显示参数：最小采样频率f_min、最大采样频率f_max、采样间隔T、抽样频率间隔T_f、时间序列的长度N。

(10)采样时间t＝[T_f，2T_f，...，(N-1)T_f]。

(11)

取整，取不小于X的最小整数。

(11)

(12)显示频率的个数Spe_nelements；也就是采样个数。

(13)进入S变换。

(14)计算时间序列的长度n；输入为时间序列A。

(15)计算快速傅里叶变换FFT。

(16)

(17)计算均值

(18)计算采样点的S变换的值。

(19)

(20)画出时频图，如图3所示。

(21)由图3可以看出，由于受到噪声的影响，频谱分辨率较低，很多微弱信号被噪声淹没，这为信号特征的提取带来了极大地干扰，因此进行时频谱重排。

(22)如果输入的变量个数为0，则结束，否则进入(22)。

(23)x_row＝size(x，2)，x_col＝size(x，1)。

(24)如果输入参数的个数小于2。

(25)令N＝x_row。

(26)h＝G(N/4)。

(27)t_row＝size(t，2)，t_col＝size(t，1)。

(28)如果xclo不等于1；则结束，否则进入(29)。

(29)如果trow不等于1，则结束，否则进入(30)。

(30)如果2^nextpow2(N)≠N，则结束，否则进入(30)。

(31)h_row＝size(h，2)，h_col＝size(h，1)，Lh＝(h_row-1)/2。

(32)如果h_col≠1||h_row/2≠0，则进入结束，否则进入(33)。

(33)如果t_col＝1；则Dt＝1；否则进入(34)。

(34)D＝t(k)-t(k-1)。

(35)Dt＝Dmin。

(36)

(37)tfr＝|tfr|²

(38)绘制重排后的时频谱图，如图4所示，

微震信号特征提取。具体实现步骤如下：首先对重排谱矩阵

做非负矩阵分解，得出分解矩阵，步骤如下：

(39)为了方便描述，做变量代换，如式(3)：

V^T _N×M＝rtft (3)

(40)求V^T _N×M的分解矩阵W_N×r、H_r×M，具体步骤如下，如式(4)

式中，W_N×r为基向量矩阵，H_r×M为系数矩阵。

(41)采用乘法迭代法不断更新W和H矩阵的元素，迭代式(5)：

式中a表示当前迭代次数，其中i≤m，j≤n，a≤min(m，n)。

(42)分解矩阵W_N×r、H_r×M进行特征提取，也即是对微震信号的特征提取。

(43)提取特征一：尖锐度SH，尖锐度(Sharpness)能够反应信号的频域能量信息如式(6)：

其中

是每一个频域基向量{w_i}_{i＝1，...，r}的傅里叶变换(FT)，而

(44)提取特征二：导数平方和SD。导数平方和能够直观地对向量中元素的均匀性进行度量。频域基向量{w_i}_{i＝1，...，r}的导数平方和定义如式(7)：

式中，w_i′(n)为{w_i}_{i＝1，...，r}中相邻元素之差即w_i′(n)＝w_i(n+1)-w_i(n)n＝1，...，N-1。

(45)提取特征三：信息熵EN。信息熵反映序列的平均信息量，分别采用式(8)和式(9)计算频域基向量{w_i}_{i＝1，...，r}和时域位置向量{h_i}_{i＝1，...，r}的信息熵EN。

(46)提取特征四：稀疏度SP。稀疏度主要反映元素之间偏移量，对{h_i}_{i＝1，...，r}计算稀疏度如式(8)

式中，M为h_i的长度。当向量中所有元素相等时，其稀疏度SP＝1。

综上所述，从每个w_i提取了(SH_wi，SD_wi，EN_wi)三个特征向量，而每个h_i提取了(EN_hi，SP_hi)两个特征参量。因此，随着i取值的不同，从微震信号的时频谱矩阵中提取的特征集F可以表示为式(9)：

(47)F＝(SH_wi，SD_wi，EN_wi，EN_hi，SP_hi；...)

(48)选定训练集和测试集。

(49)数据预处理：对载入的矩阵做归一化处理，即：

其中x，y∈Rⁿ；x_min＝min(x)；x_max＝max(x)，归一化的效果是原始数据被规整到[0，1]范围内，即y_i∈[0，1]，i＝1，2...n；这种归一化方式称为[0，1]区间归一化。

(50)训练和预测：用训练集对SVM分类器进行训练，用得到的模型对测试集进行标签预测。预测结果如图4所示。

Claims (2)

1.一种基于NMF的微震弱信号识别方法，其特征在于，具体包括如下步骤：

1)对于采集的任意时变信号x(t)进行S变换，具体如下：

上述公式(1)中，ST_x(f,b)为信号x(t)时频谱，f为采样频率，t为采样时刻，b为时间位移参数，i为虚数单位；

2)由于受到海森堡不确定性原则的限制，其时频分辨率有限，时频谱分布在f₀附近的一定范围内，有一定的虚假频带宽度，因此在S变换的基础上，引入时频谱的重排，以提高弱振幅信号的时频分辨率，重排时频谱如下：

上述公式(2)中

是信号x(t)的重排谱矩阵，f_k是S变换的离散频率，

为重排重心，

为重排带宽；

3)对步骤2)得到的重排谱矩阵

做非负矩阵分解，分解得到频域基向量W_N×r和时域位置向量H_r×M；

4)对步骤3)得到的频域基向量W_N×r和时域位置向量H_r×M，提取尖锐度、导数平方和、信息熵和稀疏度构造微震信息的特征空间；

所述的尖锐度SH，能够反应信号的频域能量信息如下：

上述公式(3)中，

为每一个频域基向量{w_i}_i＝1,...,r的傅里叶变换，而

所述的导数平方和SD，能够直观地对向量中元素的均匀性进行度量，频域基向量{w_i}_i＝1,...,r的导数平方和定义如下：

上述公式(4)中，w_i′(n)为{w_i}_{i＝1，...，r}中相邻元素之差即：w_i′(n)＝w_i(n+1)-w_i(n)，n＝1，...，N-1；

所述的信息熵EN，反映序列的平均信息量，分别采用下式计算频域基向量{w_i}_{i＝1，...，r}和时域位置向量{h_i}_{i＝1，...，r}的信息熵EN：

稀疏度SP，反映元素之间偏移量，对{h_i}_{i＝1，...，r}计算稀疏度如下：

上述公式(7)中，M为h_i的长度，当向量中所有元素相等时，其稀疏度SP＝1；

5)从每个w_i提取(SH_wi，SD_wi，EN_wi)三个特征向量，每个h_i提取(EN_hi，SP_hi)两个特征参量，随着i取值的不同，从微震信号的时频谱矩阵中提取的特征集F表示如下：

F＝(SH_wi，SD_wi，EN_wi，EN_hi，SP_hj；...) (8)

6)对步骤5)的微震信号进行分类，以公式(8)构造的特征集为输入，输出微震信号的标签，其中岩石破裂信号的标签为1，爆破信号的标签为2，机械噪声的标签为3；

7)对微震信号x(t)进行分类，分类器的最优分类面表达式如下：

s.t y_i(wφ(x_i)+b)≥1-ζ_i

上述公式(9)中，w和b分别为最优分类面的权向量和偏差；ζ_i为松弛变量，ζ_i≥0；C为惩罚系数；

φ(x_i)由式(9)得到：

上述公式(10)中，r为控制高斯核宽度的参数，采用网格搜索算法确定联合参数c和r；

8)分别选取微震信号特征集F的1-100组数据作为模式识别的训练组，101-200组数据作为预测组，经分类器输出训练组微震信号的类别标签；

经过上述步骤，完成微震信号的分类，正确率达到98％。

2.根据权利要求1所述的一种基于NMF的微震弱信号识别方法，其特征在于，步骤3)中，具体步骤如下：

3.1)对

做变量代换，具体如下：

3.2)求V^T _N×M的分解矩阵W_N×r、H_r×M，具体步骤如下：

式中，W_N×r为基向量矩阵，H_r×M为系数矩阵；

3.3)采用乘法迭代法不断更新W和H矩阵的元素，迭代式(13)：

上述公式(13)中，a表示当前迭代次数，其中i≤m，j≤n，a≤min(m，n)。

Images