一种岩体破裂信号与爆破振动信号的非线性识别方法与流程

本发明涉及一种岩体破裂信号与爆破振动信号的识别方法，尤其是涉及一种非线性的岩体破裂信号与爆破振动信号识别方法。

背景技术：
微震监测作为一种有效的地压监测手段在国内外得到了广泛应用，而岩体破裂信号与爆破振动信号的识别对微震监测具有重要意义。但岩体破裂信号与爆破振动信号相似性较大，且受众多噪音信号干扰，自动识别较为困难。目前主要采用人工识别岩体微震信号，但人工识别易受个人因素影响，且识别数量有限，限制了微震监测的实时分析。目前针对岩体破裂信号与爆破振动信号自动识别的方法可分为：多参数统计法、机器学习法和波形时频法，这些方法通常包括特征提取和特征识别两个过程。常用的特征提取包括震级、能量、视应力、视体积、静态应力降、动态应力降和时域波形特征参数(振幅、频率)等，特征识别方法包括fisher分类法、Logistic回归法、随机森林法、神经网络法、支持向量机法和贝叶斯分类法等。常用的特征参数大多难以自动得到，波形时域特征参数是在单一尺度上分析得到的，其信息量较少，限制了岩体微震信号的自动识别。波形时频法能够在多尺度上获取岩体微震信号的信息，在信号分析领域中得到了广泛的应用。波形时频法的多尺度分析主要是借助小波分析、小波包分析和频率切片小波分析。唐守锋等(2011)根据采样定理和Mallat算法确定了小波分解最大尺度，并提出采用小波特征能谱系数作为煤岩破裂微震信号识别的定量表征；朱权洁等(2012a)运用小波包分析对微震信号进行了5层多尺度分解，并对岩石破裂信号和爆破震动信号的能量分布进行了对比分析；朱权洁等(2012b)结合小波分析与分形理论对爆破振动、岩石破裂及电磁干扰3类信号进行了5层小波包分解，并以筛选后的23个小波分形盒维数作为支持向量机识别的特征向量；赵国彦等(2015)采用频率切片小波变换对矿山岩体破裂信号和爆破振动信号特定频带能量比和相关系数进行了研究。小波分析、小波包分析和频率切片小波分析具有较好的自适应性，但很容易受到信号中相邻谐波成分的交叠影响，致使不同频带信号存在混叠。可见现有的岩体破裂信号与爆破振动信号识别方法存在较大的局限，需要研究一种适用性强、准确性高的自动识别方法。

技术实现要素：
本发明所要解决的技术问题是提供一种岩体破裂信号与爆破振动信号的非线性识别方法，该岩体微震信号非线性识别方法适用性强、准确性高。发明的技术解决方案如下：一种岩体破裂信号与爆破振动信号的非线性识别方法，包括以下步骤：步骤1：导入原始岩体微震信号时间序列x(n)导入原始岩体微震信号的时间序列x(n),n＝1,2,…,N，其中N为岩体微震信号的总采样点数，取N＝4000～7000，岩体微震信号采样频率f＝4000～7000Hz；步骤2：EMD分解归一化后的岩体微震信号；2.1采用公式x*(n)＝x(n)/|xmax(n)|归一化处理原始岩体微震信号，其中x*(n)为归一化后的岩体微震信号，|xmax(n)|为最大振幅的原始岩体微震信号的绝对值；2.2采用EMD分解(经验模态分解)公式(1)对归一化后的岩体微震信号进行分解，得到本征模态分量IMFj；式中，x*(n)为归一化后的岩体微震信号，IMFj为EMD分解得到的第j个本征模态分量，IMFj(n)为IMFj第n点的值，rm(n)为EMD分解得到的残余项，m为本征模态分量的个数；步骤3：SVD分解得到特征矩阵奇异值σi(i＝1,2,…,r)3.1采用式(2)计算各个本征模态分量IMFj与归一化后的岩体微震信号的相关系数co(j)，进而根据相关系数co(j)的大小筛选得到主要本征模态分量；其中，3.2采用SVD分解公式(3)(奇异值分解)计算筛选后得到的本征模态分量构成矩阵X＝[c1c2…cr]T的奇异值σi(i＝1,2,…,r)，其中，c1,c2,…,cr为筛选后得到的本征模态分量，r为筛选后得到的本征模态分量个数，T为矩阵的转置；式中，U、V分别为r×r和N×N阶的正交矩阵；S为r×N阶的斜对角线矩阵，即其中，σ＝diag(σ1,σ2,…,σr)，是对角线元素为σ1,σ2,…,σr的对角矩阵；σi(i＝1,2,…,r)为矩阵X的奇异值，且σ1≥σ2≥…≥σr；步骤4：使用Logistic回归模型计算信号为爆破振动信号的概率p(Z)；分别选取岩体破裂信号与爆破振动信号的M组数据作为Logistic回归模型的训练数据，将步骤3.2求得的σ1,σ2,…,σr作为岩体微震信号识别的自变量；Z作为岩体微震信号识别的因变量，设定Z＝1表示信号为爆破振动信号，Z＝0表示信号为岩体破裂信号，使用最大似然估计法求得Logistic回归模型的参数β0,β1,β2,…,βr；Logistic回归模型为：p(Z)＝1/(1+exp(-Z))＝1/(1+exp(-(β0+β1·σ1+β2·σ2+…+βr·σr)))(4)其中，β0为常数项，β1,β2,…,βr为与自变量σ1,σ2,…,σr有关的参数；步骤5：根据p(Z)的大小识别岩体微震信号所述步骤3中，筛选相关系数co(j)≥0.03的本征模态分量作为主要本征模态分量。所述步骤3中，取p(Z)＝0.5作为Logistic回归模型识别岩体微震信号的分界值，p(Z)>0.5识别为爆破振动信号；p(Z)≤0.5识别为岩体破裂信号。所述步骤1中，岩体微震信号的总采样点数N取为5000，采样频率f取为6000Hz。所述步骤4中，M取值为100。有益效果：本发明的一种岩体破裂信号与爆破振动信号的非线性识别方法，包括如下步骤：导入岩体微震信号的时间序列x(n),n＝1,2,…,N。设定岩体破裂信号标识类别为0，爆破振动信号标识类别为1；EMD分解得到归一化岩体微震信号的本征模态分量，并筛选得到主要本征模态分量c1,c2,…,cr。其中，r为筛选后本征模态分量的个数；SVD分解矩阵[c1c2…cr]T，得到其奇异值σi(i＝1,2,…,r)；Logistic模型计算爆破振动信号的概率p(Z)。其中，p(Z)＝1/(1+exp(-Z))＝1/(1+exp(-(β0+β1·σ1+β2·σ2+…+βr·σr)))；识别岩体微震信号：p(Z)>0.5爆破振动信号；p(Z)≤0.5岩体破裂信号。本发明借助EMD自适应分解得到岩体微震信号的本征模态分量，使岩体微震信号具有了多尺度信息，再利用SVD对筛选后本征模态分量构成的矩阵进行分解，得到其奇异值，实现数据降维和特征提取，有效解决了单一时间尺度分析信息量较少的问题。经验模态分解从信号自身局部特征出发进行自适应的时频分解，是一种更有效的时频局部化分析方法，且适用于非线性非平稳信号的分析。此外，Logistic模型回归能有效解决多参数间的非线性关系，使用简便，计算量小，使得分类识别具有较高的准确率。此方法具有适用性强、准确性高等特点。附图说明图1是本发明所述方法流程图。图2是岩体微震信号EMD分解过程图。其中，(a)为原始岩体微震信号，(b)为归一化岩体微震信号，(c)为EMD分解归一化岩体微震信号得到的本征模态分量。图3是各本征模态分量与归一化岩体微震信号的相关系数图。图4是SVD分解筛选后得到的本征模态分量构成矩阵得到的奇异值图。图5是EMD_SVD工程应用得到的奇异值箱型图。图6是Logistic模型的训练和预测效果图具体实施方式下面将结合附图1～6，对本发明提出的一种岩体破裂信号与爆破振动信号的非线性识别方法作进一步说明。本发明算法思想的描述如下：本发明借助EMD自适应分解得到岩体微震信号的本征模态分量，使岩体微震信号具有了多尺度信息，再利用SVD对筛选后本征模态分量构成的矩阵进行分解，得到其奇异值，实现数据降维和特征提取，有效解决了单一时间尺度分析信息量较少的问题。此外，Logistic模型回归能有效解决多参数间的非线性关系，使得分类识别具有较高的准确率。步骤1：导入岩体微震信号时间序列x(n)导入岩体微震信号的时间序列x(n),n＝1,2,…,N，其中N为岩体微震信号的总采样点数，取N＝4000～7000，岩体微震信号采样频率f＝4000～7000Hz；步骤2：EMD分解归一化岩体微震信号x(n)采用公式x*(n)＝x(n)/|xmax(n)|归一化处理岩体微震信号，其中x*(n)为归一化后岩体微震信号，|xmax(j)|为原信号最大振幅的绝对值。再采用式(1)经验模态分解(EMD)对归一化岩体微震信号进行分解，得到本征模态分量IMFj；式中x*(n)为归一化后的岩体微震信号，IMFj为EMD分解得到的第j个本征模态分量，IMFj(n)为IMFj第n点的值，rm(n)为EMD分解得到的残余项，m为本征模态分量的个数。步骤3：SVD分解得到特征矩阵奇异值σi(i＝1,2,…,r)采用式(2)计算本征模态分量IMFj与归一化信号的相关系数co(j)，进而根据相关系数的大小筛选得到主要本征模态分量。其中，采用式(3)奇异值分解(SVD)计算筛选后本征模态分量构成矩阵[c1c2…cr]T的奇异值σi(i＝1,2,…,r)。其中，c1,c2,…,cr为筛选后得到的本征模态分量，r为筛选后得到的本征模态分量个数，T为矩阵的转置。式中X＝[c1c2…cr]T；U、V分别为r×r和N×N阶的正交矩阵；S为r×N阶的斜对角线矩阵，即其中，σ＝diag(σ1,σ2,…,σr)，σi(i＝1,2,…,r)为矩阵的奇异值，且σ1≥σ2≥L≥σr。步骤4：Logistic模型计算爆破振动信号的概率p(Z)分别选取岩体破裂信号与爆破振动信号的M组数据作为Logistic回归模型的训练数据，将步骤3.2求得的σ1,σ2,…,σr作为岩体微震信号识别的自变量；Z作为岩体微震信号识别的因变量，设定Z＝1表示信号为爆破振动信号，Z＝0表示信号为岩体破裂信号，使用最大似然估计法求得Logistic回归模型的参数β0,β1,β2,…,βr；Logistic回归模型为：p(Z)＝1/(1+exp(-Z))＝1/(1+exp(-(β0+β1·σ1+β2·σ2+…+βr·σr)))(4)其中，β0为常数项，β1,β2,…,βr为与自变量σ1,σ2,…,σr有关的参数；步骤5：根据p(Z)识别岩体微震信号识别岩体微震信号：p(Z)>0.5为爆破振动信号；p(Z)≤0.5为岩体破裂信号。实施例1：图2是岩体微震信号EMD分解过程图。其中，(a)为原始岩体微震信号，(b)为归一化岩体微震信号，(c)为EMD分解归一化岩体微震信号得到的本征模态分量。图中岩体微震信号总采样点数N＝5000，采样频率f＝6000Hz。由图2(c)知归一化岩体微震信号EMD分解得到了7个本征模态分量IMFj(j＝1,2,…,7)，扩大了信号的分析尺度，但增加了数据量，且可能存在虚假分量。图3是各本征模态分量与归一化岩体微震信号的相关系数图。由图3知IMF1～IMF3与归一化信号的相关系数较大，IMF4～IMF7与归一化信号的相关系数较小，且IMF7与归一化信号的相关系数仅为0.02。由徐锋等(2014)的研究成果-相关系数小于0.03的IMF分量即为虚假分量，可知选择IMF1～IMF6作为主要本征模态分量是较为合理的。图4是SVD分解筛选后得到的本征模态分量构成矩阵得到的奇异值图。由图4知SVD分解降低了矩阵的维数，并实现了矩阵的特征提取。图5是EMD_SVD工程应用得到的奇异值箱型图。岩体破裂信号类别标识为0，爆破振动信号类别标识为1。工程应用的200个岩体破裂信号和200个爆破振动信号为开阳磷矿用沙坝矿IMS微震监测系统中随机抽取得到的每个微震事件中触发最早的信号(每个微震事件包含多个微震信号)。每个微震信号总采样点数N＝5000，采样频率f＝6000Hz。由图5知岩体破裂信号与爆破振动信号的奇异值σ1、σ2和σ3差异明显，而σ4、σ5和σ6存在一定差异，但差异不明显。鉴于很难找到一个或多个奇异值分界值识别岩体破裂信号与爆破振动信号，本专利采用一种非线性方法展开识别—Logistic回归法。图6是Logistic模型的训练和预测效果图。分别选取岩体破裂信号与爆破振动信号的1～100组数据作为Logistic模型的训练数据，再分别选取101～200组数据作为Logistic模型的检验数据。Logistic模型回归得到爆破振动信号的概率p(Z)＝1/(1+exp(-Z))＝1/(1+exp(-(8.459-0.777·σ1-0.198·σ2-0.359·σ3+0.027·σ4-0.186·σ5+2.171·σ6)))，进而得到岩体微震信号的训练和检验效果(图6和表1)。为更直观的理解本专利，将25个爆破振动信号与岩体破裂信号的奇异值原始数据和预测结果分别汇于表2和表3。由图6和表1知：Logistic模型的训练和预测效果均非常好，训练组准确率达到了92.5％，检验组准确率达到了86.5％，整体准确率达到了89.5％。可见，基于EMD_SVD和Logistic模型的特征提取和识别方法是一种非常有效的岩体微震识别方法。表1Logistic模型的训练和预测效果统计表表225组爆破振动信号的奇异值原始数据和预测结果表325组岩体破裂信号的奇异值原始数据和预测结果以上所述仅为本发明的实施例而已，并不用以限制本发明，凡在本发明精神和原则之内，所作任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。