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

一种基于nmf的微震弱信号识别方法 Download PDF

Info

Publication number
CN107992802A
CN107992802A CN201711107747.7A CN201711107747A CN107992802A CN 107992802 A CN107992802 A CN 107992802A CN 201711107747 A CN201711107747 A CN 201711107747A CN 107992802 A CN107992802 A CN 107992802A
Authority
CN
China
Prior art keywords
mrow
msub
msup
munderover
frequency
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201711107747.7A
Other languages
English (en)
Other versions
CN107992802B (zh
Inventor
张法全
王海飞
肖海林
毛学港
王国富
叶金才
王小红
贾小波
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guilin University of Electronic Technology
Original Assignee
Guilin University of Electronic Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Guilin University of Electronic Technology filed Critical Guilin University of Electronic Technology
Priority to CN201711107747.7A priority Critical patent/CN107992802B/zh
Publication of CN107992802A publication Critical patent/CN107992802A/zh
Application granted granted Critical
Publication of CN107992802B publication Critical patent/CN107992802B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/288—Event detection in seismic signals, e.g. microseismics
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2411—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on the proximity to a decision surface, e.g. support vector machines
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/12—Classification; Matching

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Biology (AREA)
  • Environmental & Geological Engineering (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Remote Sensing (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Emergency Management (AREA)
  • Business, Economics & Management (AREA)
  • Signal Processing (AREA)
  • Acoustics & Sound (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

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)中，STxb,f为信号xt时频谱，f为采样频率，t为采样时刻，b为时 间位移参数，i为虚数单位；
2)由于受到海森堡不确定性原则的限制，其时频分辨率有限。时频谱分布在f0附近的一 定范围内，有一定的虚假频带宽度，因此在S变换的基础上，引入时频谱的重排，以提高弱 振幅信号的时频分辨率，重排时频谱如下：
上述公式(2)中是信号x(t)的重排谱矩阵，fk是S变换的离散频率，为重 排重心，为重排带宽；重排的意义：重排得到的谱图，它在任和点处的频谱值等于该点附近谱图值的和，达到提高时频分辨率的目的；
3)对步骤2)得到的重排谱矩阵做非负矩阵分解，分解得到频域基向量WN×r和时域位置向量Hr×M；
4)对步骤3)得到的频域基向量WN×r和时域位置向量Hr×M，提取尖锐度、导数平方和、 信息熵和稀疏度构造微震信息的特征空间；
所述的尖锐度SH，能够反应信号的频域能量信息如下：
上述公式(3)中，为每一个频域基向量{wi}i＝1,...,r的傅里叶 变换(FT)，而
所述的导数平方和SD，能够直观地对向量中元素的均匀性进行度量，频域基向量{wi}i＝1,...,r的导数平方和定义如下：
上述公式(4)中，wi′n为{wi}i＝1,...,r中相邻元素之差即：
所述的信息熵EN，反映序列的平均信息量，分别采用下式计算频域基向量{wi}i＝1,...,r和 时域位置向量{hi}i＝1,...,r的信息熵EN：
稀疏度SP，反映元素之间偏移量，对{hi}i＝1,...,r计算稀疏度如下：
上述公式(7)中，M为hi的长度。当向量中所有元素相等时，其稀疏度SP＝1；
5)从每个wi提取SHwi,SDwi,ENwi三个特征向量，每个hi提取ENhi,SPhi两个特征参量，随着i取值的不同，从微震信号的时频谱矩阵中提取的特征集F可以表示如下：
F＝SHwi,SDwi,ENwi,ENhi,SPhi；... (8)
6)对步骤5)的微震信号进行分类，以公式(8)构造的特征集为输入，输出微震信号的标签，其中岩石破裂信号的标签为1，爆破信号的标签为2，机械噪声的标签为3；
7)对微震信号x(t)进行分类，分类器的最优分类面表达式如下：
s.t yiwφxi)+b≥1-ζi
上述公式(9)中，w和b分别为最优分类面的权向量和偏差；ζi为松弛变量，ζi≥0； C为惩罚系数；
φxi由式(9)可得到：
上述公式(10)中，r为控制高斯核宽度的参数，采用网格搜索算法确定联合参数c和r；
8)分别选取微震信号特征集F的1-100组数据作为模式识别的训练组，101-200组数据 作为预测组，经分类器输出训练组微震信号的类别标签；
经过上述步骤，完成微震信号的分类，正确率达到98％。
步骤3)中，具体步骤如下：
3.1)对做变量代换，具体如下：
3.2)求VT N×M的分解矩阵WN×r、Hr×M，具体步骤如下：
式中，WN×r为基向量矩阵，Hr×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＝x1,x2,...,xn T。
(2)对采集的数据，进行等间隔采样用matlab绘制了如图2所示其波形 图。
(3)对信号x(n)进行重排S变换。
(4)如果x的列大于行，则进入步骤(5)，否则进入步骤(6)。
(5)对x做转置，并赋给A本身x＝xT；保证x是一个列向量。
(6)如果输入参数的个数为1，进入步骤(7)，否则进入步骤(8)。
(7)令最小采样频率为fmin＝0，最大采样频率fmax＝fNyquist，采样间隔T＝1。
(8)显示参数：最小采样频率fmin、最大采样频率fmax、采样间隔T、抽样频率间隔Tf、时间序列的长度N。
(10)采样时间t＝Tf,2Tf,...,N-1Tf。
(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)
(24)如果输入参数的个数小于2。
(25)令N＝xrow。
(26)h＝G N/4。
(27)trow＝size t,2,tcol＝size t,1。
(28)如果xclo不等于1；则结束，否则进入(29)。
(29)如果trow不等于1,则结束，否则进入(30)。
(30)如果则结束，否则进入(30)。
(31)hrow＝size h,2,hcol＝size h,1,Lh＝hrow-1/2。
(32)如果hcol≠1||hrow/2≠0，则进入结束，否则进入(33)。
(33)如果tcol＝1；则Dt＝1；否则进入(34)。
(34)
(35)Dt＝Dmin。
(36)
(37)tfr＝|tfr|2
(38)绘制重排后的时频谱图，如图4所示，
微震信号特征提取。具体实现步骤如下：首先对重排谱矩阵做非负矩阵分解，得出分解矩阵，步骤如下：
(39)为了方便描述，做变量代换，如式(3)：
VT N×M＝rtft (3)
(40)求VT N×M的分解矩阵WN×r、Hr×M，具体步骤如下，如式(4)
式中，WN×r为基向量矩阵，Hr×M为系数矩阵。
(41)采用乘法迭代法不断更新W和H矩阵的元素，迭代式(5)：
式中a表示当前迭代次数，其中i≤m,j≤n,a≤min m,n)。
(42)分解矩阵WN×r、Hr×M进行特征提取，也即是对微震信号的特征提取。
(43)提取特征一：尖锐度SH，尖锐度(Sharpness)能够反应信号的频域能量信息如式 (6)：
其中是每一个频域基向量{wi}i＝1,...,r的傅里叶变换(FT)，而
(44)提取特征二：导数平方和SD。导数平方和能够直观地对向量中元素的均匀性进行度 量。频域基向量{wi}i＝1,...,r的导数平方和定义如式(7)：
式中，n为{wi}i＝1,...,r中相邻元素之差即
(45)提取特征三：信息熵EN。信息熵反映序列的平均信息量，分别采用式(8)和式(9)计算频域基向量{wi}i＝1,...,r和时域位置向量{hi}i＝1,...,r的信息熵EN。
(46)提取特征四：稀疏度SP。稀疏度主要反映元素之间偏移量，对{hi}i＝1,...,r计算稀疏 度如式(8)
式中，M为hi的长度。当向量中所有元素相等时，其稀疏度SP＝1。
综上所述，从每个wi提取了三个特征向量，而每个hi提取了两个特征参量。因此，随着i取值的不同，从微震信号的时频谱矩阵中提取的特征集F可以 表示为式(9):
(47)F＝SHwi,SDwi,ENwi,ENhi,SPhi；...
(48)选定训练集和测试集。
(49)数据预处理：对载入的矩阵做归一化处理，即：
其中x,y∈Rn；xmin＝min x；xmax＝max x，归一化的效果是原始数据被规整到0,1范围内，即yi∈0,1,i＝1,2...n；这种归一化方式称为[0,1]区间归一化。
(50)训练和预测：用训练集对SVM分类器进行训练，用得到的模型对测试集进行标签预测。预测结果如图4所示。

Claims (2)

1.一种基于NMF的微震弱信号识别方法，其特征在于，具体包括如下步骤：
1)对于采集的任意时变信号x(t)进行S变换，具体如下：
<mrow> <msub> <mi>ST</mi> <mi>x</mi> </msub> <mi>f</mi> <mo>,</mo> <mi>b</mi> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mi>f</mi> <mo>|</mo> </mrow> <msqrt> <mrow> <mn>2</mn> <mi>&amp;pi;</mi> </mrow> </msqrt> </mfrac> <munderover> <mo>&amp;Integral;</mo> <mrow> <mo>-</mo> <mi>&amp;infin;</mi> </mrow> <mi>&amp;infin;</mi> </munderover> <msup> <mi>xte</mi> <mfrac> <mrow> <mo>-</mo> <mi>t</mi> <mo>-</mo> <msup> <mi>b</mi> <mn>2</mn> </msup> <msup> <mi>f</mi> <mn>2</mn> </msup> </mrow> <mn>2</mn> </mfrac> </msup> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>f</mi> <mi>t</mi> </mrow> </msup> <mi>d</mi> <mi>t</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
上述公式(1)中，STxb,f为信号x t时频谱，f为采样频率，t为采样时刻，b为时间位移参数，i为虚数单位；
2)由于受到海森堡不确定性原则的限制，其时频分辨率有限，时频谱分布在f0附近的一定范围内，有一定的虚假频带宽度，因此在S变换的基础上，引入时频谱的重排，以提高弱振幅信号的时频分辨率，重排时频谱如下：
<mrow> <msup> <msub> <mi>ST</mi> <mi>x</mi> </msub> <mi>r</mi> </msup> <msub> <mover> <mi>f</mi> <mo>^</mo> </mover> <mi>l</mi> </msub> <mo>,</mo> <mi>b</mi> <mo>=</mo> <mi>&amp;Delta;</mi> <msup> <msub> <mover> <mi>f</mi> <mo>^</mo> </mover> <mi>l</mi> </msub> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>&amp;times;</mo> <munder> <mi>&amp;Sigma;</mi> <mrow> <msub> <mi>f</mi> <mi>k</mi> </msub> <mo>:</mo> <mrow> <mo>|</mo> <mrow> <msub> <mover> <mi>f</mi> <mo>^</mo> </mover> <mi>l</mi> </msub> <msub> <mi>f</mi> <mi>k</mi> </msub> <mo>,</mo> <mi>b</mi> <mo>-</mo> <msub> <mover> <mi>f</mi> <mo>^</mo> </mover> <mi>l</mi> </msub> </mrow> <mo>|</mo> </mrow> <mo>&amp;le;</mo> <mi>&amp;Delta;</mi> <msub> <mover> <mi>f</mi> <mo>^</mo> </mover> <mi>l</mi> </msub> <mo>/</mo> <mn>2</mn> </mrow> </munder> <mo>|</mo> <mrow> <msub> <mi>ST</mi> <mi>x</mi> </msub> <mi>f</mi> <mo>,</mo> <mi>b</mi> </mrow> <mo>|</mo> <msub> <mi>f</mi> <mi>k</mi> </msub> <msub> <mi>&amp;Delta;f</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
上述公式(2)中是信号x(t)的重排谱矩阵，fk是S变换的离散频率，为重排重心，为重排带宽；
3)对步骤2)得到的重排谱矩阵做非负矩阵分解，分解得到频域基向量WN×r和时域位置向量Hr×M；
4)对步骤3)得到的频域基向量WN×r和时域位置向量Hr×M，提取尖锐度、导数平方和、信息熵和稀疏度构造微震信息的特征空间；
所述的尖锐度SH，能够反应信号的频域能量信息如下：
<mrow> <msub> <mi>SH</mi> <mrow> <mi>w</mi> <mi>i</mi> </mrow> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <msub> <mi>n</mi> <mn>0</mn> </msub> </mrow> <mrow> <mi>N</mi> <mo>/</mo> <mn>4</mn> </mrow> </munderover> <mo>|</mo> <mrow> <msub> <mi>W</mi> <mi>i</mi> </msub> <mi>k</mi> </mrow> <mo>|</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
上述公式(3)中，为每一个频域基向量{wi}i＝1,...,r的傅里叶变换(FT)，而
所述的导数平方和SD，能够直观地对向量中元素的均匀性进行度量，频域基向量{wi}i＝1,...,r的导数平方和定义如下：
<mrow> <msub> <mi>SD</mi> <mrow> <mi>w</mi> <mi>i</mi> </mrow> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msup> <msub> <mi>w</mi> <mi>i</mi> </msub> <mo>&amp;prime;</mo> </msup> <msup> <mi>n</mi> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
上述公式(4)中，wi′n为{wi}i＝1,...,r中相邻元素之差即：wi′n＝win+1-winn＝1,...,N-1；
所述的信息熵EN，反映序列的平均信息量，分别采用下式计算频域基向量{wi}i＝1,...,r和时域位置向量{hi}i＝1,...,r的信息熵EN：
<mrow> <msub> <mi>EN</mi> <mrow> <mi>w</mi> <mi>i</mi> </mrow> </msub> <mo>=</mo> <mo>-</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msup> <msub> <mi>w</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> <mi>n</mi> <mi> </mi> <mi>log</mi> <mi> </mi> <msup> <msub> <mi>w</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> <mi>n</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <msub> <mi>EN</mi> <mrow> <mi>h</mi> <mi>i</mi> </mrow> </msub> <mo>=</mo> <mo>-</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msup> <msub> <mi>h</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> <mi>m</mi> <mi> </mi> <mi>log</mi> <mi> </mi> <msup> <msub> <mi>h</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> <mi>m</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
稀疏度SP，反映元素之间偏移量，对{hi}i＝1,...,r计算稀疏度如下：
<mrow> <msub> <mi>SP</mi> <mrow> <mi>h</mi> <mi>i</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <msqrt> <mi>M</mi> </msqrt> <mo>-</mo> <mrow> <mo>(</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>h</mi> <mi>i</mi> </msub> <mi>m</mi> <mo>)</mo> </mrow> <mo>/</mo> <msqrt> <mrow> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msup> <msub> <mi>h</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> <mi>m</mi> </mrow> </msqrt> </mrow> <mrow> <msqrt> <mi>M</mi> </msqrt> <mo>-</mo> <mn>1</mn> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
上述公式(7)中，M为hi的长度。当向量中所有元素相等时，其稀疏度SP＝1；
5)从每个wi提取SHwi,SDwi,ENwi三个特征向量，每个hi提取ENhi,SPhi两个特征参量，随着i取值的不同，从微震信号的时频谱矩阵中提取的特征集F可以表示如下：
F＝SHwi,SDwi,ENwi,ENhi,SPhi；... (8)
6)对步骤5)的微震信号进行分类，以公式(8)构造的特征集为输入，输出微震信号的标签，其中岩石破裂信号的标签为1，爆破信号的标签为2，机械噪声的标签为3；
7)对微震信号x(t)进行分类，分类器的最优分类面表达式如下：
<mrow> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mrow> <mi>w</mi> <mo>,</mo> <mi>b</mi> <mo>,</mo> <mi>&amp;zeta;</mi> <mi>i</mi> </mrow> </munder> <mo>:</mo> <mo>|</mo> <mo>|</mo> <mi>w</mi> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>/</mo> <mn>2</mn> <mo>+</mo> <mi>C</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>&amp;zeta;</mi> <mi>i</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
s.t yiwφxi+b≥1-ζi
上述公式(9)中，w和b分别为最优分类面的权向量和偏差；ζi为松弛变量，ζi≥0；C为惩罚系数；
φxi由式(9)可得到：
<mrow> <msub> <mi>Kx</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>&amp;equiv;</mo> <msub> <mi>&amp;phi;x</mi> <mi>i</mi> </msub> <msub> <mi>&amp;phi;x</mi> <mi>j</mi> </msub> <mo>=</mo> <msup> <mi>e</mi> <mrow> <mi>r</mi> <mo>|</mo> <mo>|</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>
上述公式(10)中，r为控制高斯核宽度的参数，采用网格搜索算法确定联合参数c和r；
8)分别选取微震信号特征集F的1-100组数据作为模式识别的训练组，101-200组数据作为预测组，经分类器输出训练组微震信号的类别标签；
经过上述步骤，完成微震信号的分类，正确率达到98％。
2.根据权利要求1所述的一种基于NMF的微震弱信号识别方法，其特征在于，步骤3)中，具体步骤如下：
3.1)对做变量代换，具体如下：
<mrow> <msub> <msup> <mi>V</mi> <mi>T</mi> </msup> <mrow> <mi>N</mi> <mo>&amp;times;</mo> <mi>M</mi> </mrow> </msub> <mo>=</mo> <msup> <msub> <mi>ST</mi> <mi>x</mi> </msub> <mi>r</mi> </msup> <msub> <mover> <mi>f</mi> <mo>^</mo> </mover> <mi>l</mi> </msub> <mo>,</mo> <mi>b</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>
3.2)求VT N×M的分解矩阵WN×r、Hr×M，具体步骤如下：
<mrow> <msub> <msup> <mi>V</mi> <mi>T</mi> </msup> <mrow> <mi>N</mi> <mo>&amp;times;</mo> <mi>M</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>W</mi> <mrow> <mi>M</mi> <mo>&amp;times;</mo> <mi>r</mi> </mrow> </msub> <msub> <mi>H</mi> <mrow> <mi>r</mi> <mo>&amp;times;</mo> <mi>N</mi> </mrow> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>r</mi> </munderover> <msub> <mi>w</mi> <mi>i</mi> </msub> <msub> <mi>h</mi> <mi>i</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>
式中，WN×r为基向量矩阵，Hr×M为系数矩阵。
3.3)采用乘法迭代法不断更新W和H矩阵的元素，迭代式(13)：
<mrow> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>a</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>a</mi> </mrow> </msub> <mfrac> <mrow> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <mfrac> <msub> <mi>V</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mrow> <msub> <mi>WH</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> </mrow> </mfrac> <msub> <mi>H</mi> <mrow> <mi>a</mi> <mi>j</mi> </mrow> </msub> </mrow> <mrow> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <mi>H</mi> <mrow> <mi>a</mi> <mi>j</mi> </mrow> </msub> </mrow> </mfrac> <mo>,</mo> <msub> <mi>H</mi> <mrow> <mi>a</mi> <mi>j</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>H</mi> <mrow> <mi>a</mi> <mi>j</mi> </mrow> </msub> <mfrac> <mrow> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <mfrac> <msub> <mi>V</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mrow> <msub> <mi>WH</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> </mrow> </mfrac> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>a</mi> </mrow> </msub> </mrow> <mrow> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>a</mi> </mrow> </msub> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>
上述公式(13)中，a表示当前迭代次数，其中i≤m,j≤n,a≤min m,n。
CN201711107747.7A 2017-11-10 2017-11-10 一种基于nmf的微震弱信号识别方法 Active CN107992802B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711107747.7A CN107992802B (zh) 2017-11-10 2017-11-10 一种基于nmf的微震弱信号识别方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711107747.7A CN107992802B (zh) 2017-11-10 2017-11-10 一种基于nmf的微震弱信号识别方法

Publications (2)

Publication Number Publication Date
CN107992802A true CN107992802A (zh) 2018-05-04
CN107992802B CN107992802B (zh) 2021-12-21

Family

ID=62030750

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711107747.7A Active CN107992802B (zh) 2017-11-10 2017-11-10 一种基于nmf的微震弱信号识别方法

Country Status (1)

Country Link
CN (1) CN107992802B (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109214469A (zh) * 2018-10-24 2019-01-15 西安交通大学 一种基于非负张量分解的多源信号分离方法
CN109307886A (zh) * 2018-08-24 2019-02-05 中国石油天然气股份有限公司 多次波自适应相减方法及装置
CN110428369A (zh) * 2019-06-20 2019-11-08 中国地质大学（武汉） 基于信息熵的chnmf遥感图像解混算法
CN110458760A (zh) * 2019-06-20 2019-11-15 中国地质大学（武汉） 基于信息熵的hnmf遥感图像解混算法
CN110471104A (zh) * 2019-08-26 2019-11-19 电子科技大学 基于智能特征学习的叠后地震反射模式识别方法
CN112505481A (zh) * 2020-11-20 2021-03-16 云南电网有限责任公司普洱供电局 基于近邻传播聚类的35kV电力线路故障行波提取方法
CN112968741A (zh) * 2021-02-01 2021-06-15 中国民航大学 基于最小二乘向量机的自适应宽带压缩频谱感知算法
CN113049252A (zh) * 2021-03-25 2021-06-29 成都天佑路航轨道交通科技有限公司 一种列车轴承箱的故障检测方法
CN114155870A (zh) * 2021-12-02 2022-03-08 桂林电子科技大学 低信噪比下基于spp和nmf的环境音噪声抑制方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4649524A (en) * 1983-07-05 1987-03-10 Potash Corporation Of Saskatchewan Mining Limited Integrated acoustic network
US20030067843A1 (en) * 2001-10-05 2003-04-10 Jean-Francois Therond Method intended for detection and automatic classification, according to various selection criteria, of seismic events in an underground formation
CN103064111A (zh) * 2012-12-12 2013-04-24 中国石油天然气集团公司 一种基于形态滤波的微地震事件识别方法
CN103786069A (zh) * 2014-01-24 2014-05-14 华中科技大学 一种机械加工设备的颤振在线监测方法
CN105527650A (zh) * 2016-02-17 2016-04-27 中国科学院武汉岩土力学研究所 一种工程尺度下微震信号及p波初至自动识别算法
CN105974469A (zh) * 2016-06-30 2016-09-28 马克 岩质边坡开挖扰动作用下的微震监测预警分析系统及方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4649524A (en) * 1983-07-05 1987-03-10 Potash Corporation Of Saskatchewan Mining Limited Integrated acoustic network
US20030067843A1 (en) * 2001-10-05 2003-04-10 Jean-Francois Therond Method intended for detection and automatic classification, according to various selection criteria, of seismic events in an underground formation
CN103064111A (zh) * 2012-12-12 2013-04-24 中国石油天然气集团公司 一种基于形态滤波的微地震事件识别方法
CN103786069A (zh) * 2014-01-24 2014-05-14 华中科技大学 一种机械加工设备的颤振在线监测方法
CN105527650A (zh) * 2016-02-17 2016-04-27 中国科学院武汉岩土力学研究所 一种工程尺度下微震信号及p波初至自动识别算法
CN105974469A (zh) * 2016-06-30 2016-09-28 马克 岩质边坡开挖扰动作用下的微震监测预警分析系统及方法

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109307886A (zh) * 2018-08-24 2019-02-05 中国石油天然气股份有限公司 多次波自适应相减方法及装置
CN109214469A (zh) * 2018-10-24 2019-01-15 西安交通大学 一种基于非负张量分解的多源信号分离方法
CN110428369A (zh) * 2019-06-20 2019-11-08 中国地质大学（武汉） 基于信息熵的chnmf遥感图像解混算法
CN110458760A (zh) * 2019-06-20 2019-11-15 中国地质大学（武汉） 基于信息熵的hnmf遥感图像解混算法
CN110428369B (zh) * 2019-06-20 2021-10-08 中国地质大学（武汉） 基于信息熵的chnmf遥感图像解混方法
CN110458760B (zh) * 2019-06-20 2021-11-05 中国地质大学（武汉） 基于信息熵的hnmf遥感图像解混方法
CN110471104A (zh) * 2019-08-26 2019-11-19 电子科技大学 基于智能特征学习的叠后地震反射模式识别方法
CN112505481A (zh) * 2020-11-20 2021-03-16 云南电网有限责任公司普洱供电局 基于近邻传播聚类的35kV电力线路故障行波提取方法
CN112968741A (zh) * 2021-02-01 2021-06-15 中国民航大学 基于最小二乘向量机的自适应宽带压缩频谱感知算法
CN113049252A (zh) * 2021-03-25 2021-06-29 成都天佑路航轨道交通科技有限公司 一种列车轴承箱的故障检测方法
CN114155870A (zh) * 2021-12-02 2022-03-08 桂林电子科技大学 低信噪比下基于spp和nmf的环境音噪声抑制方法
CN114155870B (zh) * 2021-12-02 2024-08-27 桂林电子科技大学 低信噪比下基于spp和nmf的环境音噪声抑制方法

Also Published As

Publication number Publication date
CN107992802B (zh) 2021-12-21

Similar Documents

Publication Publication Date Title
CN107992802A (zh) 一种基于nmf的微震弱信号识别方法
CN105956526B (zh) 基于多尺度排列熵的低信噪比微震事件辨识方法
Liu et al. Deep learning and recognition of radar jamming based on CNN
CN109359523B (zh) 一种基于svm多分类算法的卫星导航干扰类型识别方法
CN106682615B (zh) 一种水下弱小目标检测方法
CN107301381A (zh) 基于深度学习和多任务学习策略的雷达辐射源识别方法
US10621506B2 (en) Apparatus and method for activity detection and classification from sensor data
CN103487513B (zh) 一种空间碎片撞击毁伤声发射信号类型识别方法
CN103412557A (zh) 一种适于非线性过程在线监控的工业故障检测与诊断方法
CN103743980A (zh) 基于pso优化的svm的电能质量扰动识别与分类方法
Li et al. Life grade recognition method based on supervised uncorrelated orthogonal locality preserving projection and K-nearest neighbor classifier
CN105447464A (zh) 一种基于pso的电能质量扰动识别与分类方法
CN104330258A (zh) 一种基于特征参量的滚动轴承故障灰色关联度辨识方法
CN103413134A (zh) 基于稀疏分解的地面移动目标微动信号特征提取
CN104732970A (zh) 一种基于综合特征的舰船辐射噪声识别方法
CN107392123A (zh) 一种基于相参积累消噪的射频指纹特征提取和识别方法
CN103994820B (zh) 一种基于微孔径麦克风阵列的运动目标识别方法
CN114037020A (zh) 微地震事件识别分类方法、装置、设备及存储介质
Omar et al. Fourier Domain Kernel Density Estimation-based Approach for Hail Sound classification
CN114519372B (zh) 基于支持向量机的一维距离像目标识别方法
CN110222386A (zh) 一种行星齿轮退化状态识别方法
CN111275003B (zh) 基于类最优高斯核多分类支持向量机的微震信号识别方法
Kamal et al. Novel class detection of underwater targets using self-organizing neural networks
CN102789780B (zh) 基于谱时幅度分级向量辨识环境声音事件的方法
Liu et al. Bearing performance degradation assessment using linear discriminant analysis and coupled HMM

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant