一种微震信号的特征提取方法与流程

1.本发明涉及一种微震信号的特征提取方法，属于信号处理技术领域。


背景技术：

2.微震监测系统通过分析生产活动所引起的微震事件来监测岩体的稳定状态。但矿山环境复杂，干扰因素众多，将有效的微震事件从大量的数据样本中识别出来是发挥微震监测系统预测预警的基础性问题，而特征提取是对微震信号进行准确分析和精确识别的关键步骤，具有重要的研究意义。
3.目前常用的特征提取流程是通过时频分析方法将信号分解成多个特征分量，再利用svd或者信息熵理论得到特征分量的差异性特征值。典型的时频分析方法主要有傅里叶变换(ft)，经验模态分解(emd)、小波变换(wt)和小波包变换(wpt)。这些方法在一定程度上得到了信号的特征分量，但还存在一些不足。ft适用于平稳信号的分析，但面对具有非平稳和突发瞬态的微震信号时难以进行有效分解；wpt和wt虽然适用于非平稳信号处理，但对信号的分解效果很大程度依赖于分解尺度和基函数的选择，自适应性不够好；emd自适应性好，但存在模态混叠、端点效应的问题，而且emd缺乏严格理论支撑。
4.为此，2013年gills提出一种新的时频分析方法——经验小波变换(ewt)，它是在wt理论框架下提出的，所以理论基础完备，并且借助于小波分解快速算法使得自身具有较高的计算效率。ewt首先利用傅里叶变换得到信号的频域表示，再利用局部极大值的中点作为分界频点对信号频谱进行划分，根据不同的频段构建经验小波，并获得多个窄带频率分量，具有emd的自适应性。近年来ewt在不同领域的信号分析中得到应用，如电机振动信号，心电信号、变压器振动信号等。但随着研究的深入，ewt存在的一些问题也在显现，ewt方法的有效性分解主要取决于频谱的有效分割，但微震信号频域特性复杂，传统的依靠局部极大值的频谱分割方法会将同一频带的成分划分到不同的频带，使得同一分量被划分成两个分量，造成过分割的现象，另外，传统的ewt方法在分析非平稳、局部噪声信号时，非平稳分量以及噪声分量产生的一些局部极大值可能出现并错误的保持在峰值序列中，而有效信号的极大值可能不保持在峰值序列中，导致不当的分割，进而使得微震信号的特征提取准确性低。


技术实现要素：

5.本技术的目的在于提供一种微震信号的特征提取方法，用以解决现有特征提取方法准确性的问题。
6.为实现上述目的，本技术提出了一种微震信号的特征提取方法的技术方案，包括以下步骤：
7.1)获取待特征提取的微震信号，对微震信号进行去趋势项处理，得到去趋势项后的微震信号；
8.2)计算去趋势项后的微震信号的功率谱，找出功率谱中的局部极大值点；
9.3)将局部极大值点左右相邻的极小值点作为预分割频点，计算属于不同的局部极大值点的相邻预分割频点之间的频率差值；
10.4)判断频率差值的大小，若某两相邻预分割频点之间的频率差值小于等于分割阈值，则根据这两个相邻预分割频点确定最终的分割频点，若某两相邻预分割频点之间的频率差值大于分割阈值，则这两个相邻预分割频点即为最终的分割频点；
11.5)按照最终的分割频点对功率谱进行频段的划分，进而提取特征值。
12.本发明的微震信号的特征提取方法的技术方案的有益效果是：本发明在对功率谱进行划分时，首先将局部极大值点左右相邻的极小值点作为预分割频点，接着判断这些预分割频点是否为最终的分割频点：根据不同的局部极大值点的相邻预分割频点之间的频率差值进行判断，若频率差值小于等于分割阈值，则还需要根据这两个相邻预分割频点确定最终的分割频点，若频率差值大于分割阈值，则这两个相邻预分割频点即为最终的分割频点。本发明通过极大极小值频点进行初步频谱分割，进而利用根据采样频率自适应调控的阈值对初步划分的频谱进一步进行频点优化，避免了过分割或者不当分割而形成无效分量，提高了对微震信号分割的准确性。
13.进一步的，为了减小功率谱的畸变，使得功率谱更加接近与真实谱，所述步骤2)中，计算的功率谱为最大熵功率谱。
14.进一步的，为了提高分割频点的计算效率，所述步骤4)中，若某两相邻预分割频点之间的频率差值小于等于分割阈值，则这两个相邻预分割频点的频率平均值对应的频点为最终的分割频点。
15.进一步的，所述步骤5)中，提取各频段的特征值的步骤包括：
16.构建各频段内的小波函数和尺度函数，利用傅里叶反变换得到各频段分量的时域表示；
17.根据各频段分量的时域表示构造矩阵；
18.对矩阵进行svd分解，进而提取特征值。
19.进一步的，为了减小无效频段分量的影响，得到各频段分量的时域表示后，还包括利用相关性计算对时域表示的频段进行筛选的步骤：计算各频段分量的时域表示与去趋势项后的微震信号的相关系数，若相关系数小于相关系数阈值，则剔除该频段分量的时域表示。
20.进一步的，为了充分提取每个频段分量的特征信息，更加准确的得到每个频段分量的特征值，所构造的矩阵包括每个频段分量对应的hankel矩阵，对每个hankel矩阵进行svd分解，得到的每个hankel矩阵的最大奇异值即为对应频段分量的特征值。
21.进一步的，为了更好的进行去趋势处理，步骤1)中利用最小二乘法进行去趋势项处理，趋势项为：
[0022][0023]
其中，d
t
为趋势项；f
s
为微震信号的采样率；t为采样点序号；a
i
(i＝0,1,2,...,k)为待定系数；k为趋势项的阶数。
[0024]
进一步的，最大熵功率谱的计算过程为：
[0025][0026]
其中，s(ω)为最大熵功率；ω为频率；为第m阶的预报误差方差估计；为第m阶第j次的自回归系数；f
s
为微震信号的采样率；为虚数符号。
[0027]
进一步的，通过burg递推算法计算预报误差方差估计和自回归系数。
[0028]
进一步的，相关系数的计算过程为：
[0029][0030]
其中：co(k)为第k个频段分量的相关系数；x'(t)为去趋势项后的微震信号，为x'(t)的平均值；x
k
(t)为第k个频段分量的时域表示；为各频段分量的时域表示的平均值；n为频段分量的数量；t为采样点序号；m为微震信号的采样点数量。
附图说明
[0031]
图1是本发明微震信号的特征提取方法的流程图；
[0032]
图2是本发明小波函数和尺度函数的滤波器组图。
具体实施方式
[0033]
微震信号的特征提取方法实施例：
[0034]
本发明的主要构思在于，基于微震信号的特点，为了避免过分割的现象，本发明在得到最大熵功率谱中的局部极大值点后，以局部极大值点左右相邻的极小值点作为预分割频点，对预分割频点进行判断：计算属于不同的局部极大值点的相邻预分割频点之间的距离，若距离小于等于分割阈值，则根据这两个相邻预分割频点再次确定最终的分割频点，若距离大于分割阈值，则这两个相邻预分割频点即为最终的分割频点，根据最终的分割频点对最大熵功率谱进行频段划分，进而得到各频段的特征值。
[0035]
具体的，微震信号的特征提取方法如图1所示，包括以下步骤：
[0036]
1)获取检测到的微震信号，对微震信号进行模数转换后进行m点截断，得到m点采样值序列的离散信号x(t)，(t＝1,2,...m)，将离散信号x(t)作为待特征提取的微震信号。
[0037]
2)利用最小二乘法对步骤1)的得到的离散信号x(t)进行去趋势项处理，得到去趋势项的离散信号x'(t)。
[0038]
趋势项是由于传感器的零点漂移，基础运动等引起的信号波形偏移，本步骤中，利用最小二乘法进行去趋势项处理的过程如下：
[0039]
利用k阶多项式(也即估计值)来拟合离散信号x(t)中的趋势项d
t
：
[0040]
[0041]
其中，f
s
为微震信号的采样率；a
i
(i＝0,1,2,...,k)为待定系数；t为采样点序号；k为趋势项的阶数。待定系数a
i
的确定依据是使得k阶多项式d
t
与x(t)之间的误差平方和最小，有如下公式：
[0042][0043]
令e对a
i
的偏导数为零，可产生k 1元线性方程组，第j元线性方程为：
[0044][0045]
通过求解(3)式，可得到k 1个待定系数a
i
，将求出的待定系数a
i
代入估计值d
t
中，即可得到估计值d
t
。
[0046]
去趋势项的离散信号x'(t)为：x'(t)＝x(t)
‑
d
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0047]
当然作为其他实施方式，还可以采用插值法进行去趋势项处理，但插值法需严格的通过数据，当数据带有误差时，就会使预测值偏离真实数据，产生过拟合现象。
[0048]
3)根据步骤2)中得到的去趋势项的离散信号x'(t)，结合伯格(burg)算法计算出自回归系数和预报误差方差估计值进而将两个参数带入最大熵功率谱s(ω)的计算公式中计算出最大熵功率谱s(ω)，并将频率轴归一化到[0,π]。
[0049]
本步骤中，最大熵功率谱s(ω)的计算公式如下：
[0050][0051]
其中，为第m阶的预报误差方差估计；为第m阶第j次的自回归系数；为虚数符号。参数和的确定需要利用burg递推算法，计算过程如下：
[0052]
(1)计算预报误差方差估计的初始值
[0053][0054][0055]
式中，为前向预测误差的初始值；为后向预测误差的初始值；
[0056]
(2)依据如下公式计算第m阶的反射系数k
m
：
[0057][0058]
式中，为第m阶的前向预测误差；为第m阶的后向预测误差；p是总阶数量；
[0059]
(3)计算第m阶第j次的自回归系数
[0060][0061]
(4)计算第m阶的预报误差方差估计
[0062][0063]
(5)计算预测滤波器输出：
[0064][0065][0066]
式中，为k
m
的复共轭。
[0067]
(6)令m＝m 1，重复步骤(2)～(5)，直到m＝p。
[0068]
将求得的参数代入公式(5)中可求得最大熵功率谱。
[0069]
当然，作为其他实施方式，还存在自相关法和改进协方差法算出自回归系数和预报误差方差估计值，但改进协方差法计算繁琐，编程较困难，自相关法计算最简单，谱估计的分辨率相对较差，burg法是较为通用的方法，计算不太复杂，谱估计质量较好。
[0070]
4)找出功率谱中的局部极大值点，将局部极大值点左右相邻的极小值点作为预分割频点，计算属于不同的局部极大值点的相邻预分割频点之间的距离d(这里的距离是指相邻预分割频点之间的频率差值)，若距离d小于等于分割阈值α，则根据这两个相邻预分割频点确定最终的分割频点，若距离d大于分割阈值α，则这两个相邻预分割频点即为最终的分割频点，根据最终的分割频点对最大熵功率谱s(ω)进行频段划分。
[0071]
本步骤中的局部极大值点的确定过程为：假设需要分割的频段数为n个，用ω
n
(n＝1,2,...,n 1)表示频段分割点。找出最大熵功率谱s(ω)中所有的局部极大值，从大到小排列，保留前n个局部极大值对应的频段分割点为局部极大值点。
[0072]
本步骤中，若属于不同的局部极大值点的相邻预分割频点之间的距离d小于分割阈值α，则将这两个相邻预分割频点的频率平均值作为最终的分割频点。例如：相邻局部极大值点分别为a1和a2，a1位于a2的左侧，也即a1＜a2；局部极大值点a1左右相邻的极小值点分别为ω1，ω2，局部极大值点a2左右相邻的极小值点分别为ω3，ω4，也即ω1，ω2，ω3，ω4均为预分割频点；
[0073]
属于a1和a2的相邻预分割频点为ω2，ω3，计算ω2，ω3之间的距离d，若d≤α(经过多次实验验证，α一般取10/f
s
效果较好)，则剔除ω2，ω3这两个预分割频点，将作为最终的分割频点；若d>α，则确定ω2，ω3为最终的分割频点。
[0074]
5)在划分的各个频段内构建对应的小波函数ψ
k
(ω)和尺度函数这里的k表示第k个频段，k＝0,1,2，
…
，n
‑
1。
[0075]
本步骤中，第k个频段的小波函数ψ
k
(ω)和尺度函数的频率响应如图2所示，计算公式如下：
[0076][0077][0078]
其中，ω为归一化频率，ω
n
∈[0,π]为分割频点，t
n
＝γω
n
(0<γ<1)为经验小波滤波器过渡带宽度，其中β(x)＝x4(35
‑
85x 70x2‑
20x3)，这里的x代表自变量。
[0079]
6)根据各个频段内的小波函数和尺度函数，利用傅里叶反变换得到尺度系数和细节系数，进而得到各频段分量的时域表示。
[0080]
应用傅里叶反变换计算s(ω)
×
ψ
k
(ω)和从而得到各个分量x0(t),x1(t),...x
n
‑1(t)的时域表示。
[0081]
具体的，本步骤中通过信号x'(t)与小波函数和尺度函数的内积运算获取信号分解的尺度系数k＝0，也即和细节系数计算公式如下：
[0082][0083]
其中，和分别为经验尺度函数和经验小波函数ψ
k
(t)的复共轭；和ψ
k
(ω)分别为和ψ
k
(t)的fourier变换；f
‑1[
·
]为fourier反变换(即傅里叶反变换)。
[0084]
至此，可计算各个频段分量的时域表示：
[0085][0086]
式中：“*”代表卷积。
[0087]
7)利用相关性分析对x0(t),x1(t),...x
n
‑1(t)共n个频段分量的时域表示进行筛选。
[0088]
本步骤中，对得到的n个频段分量的时域表示进行筛选，筛选是依据各个频段分量与去趋势项后的微震信号相关系数值的大小，当频段分量的时域表示与去趋势项后的微震信号的相关系数值小于0.03(相关系数阈值)时，即可认为该频段分量为无效分量，进行剔除。相关性计算公式为：
[0089][0090]
其中：co(k)为第k个频段分量的相关系数；x'(t)为去趋势项后的微震信号，为x'(t)的平均值；x
k
(t)为第k个频段分量的时域表示；为各频段分量的时域表示的平均值；n为频段分量的数量；t为采样点序号；m为微震信号的采样点数量。
[0091]
8)利用筛选后的每个频段分量的时域表示构造对应的hankel矩阵，并对每个矩阵进行svd分解，取每个矩阵最大的奇异值作为表征该分量的特征值。
[0092]
本步骤中，hankel矩阵h
α
的构造形式如下：
[0093][0094]
式中：x
k,n 1
为第k个频段分量第n 1个采样点的时域表示；且1<n<m；m是截断的采样点长度，其中矩阵h
α
行列数按以下规则选取：
[0095]
(1)若m为偶数，hankel矩阵的行数为m＝m/2 1，列数为n＝m/2；
[0096]
(2)若m为奇数，hankel矩阵的行数和列数都为m＝(m 1)/2；
[0097]
分别对生成的m
×
n阶矩阵h
α
进行svd，svd分解形式如下：
[0098][0099]
式中：r(r≤min(m,n))为矩阵h
α
的秩，u∈r
m
×
m
和v∈r
n
×
n
为正交矩阵，u
i
和v
i
为矩阵u和v的第i列向量，其中矩阵并且σ
i
为第i列向量的奇异值；h
i
为第i列向量对应的矩阵，σ＝diag(σ1,σ2,...,σ
p
)(p＝min(m,n))，将最大奇异值σ作为表征对应矩阵h
α
的特征值。
[0100]
上述实施例中，关于功率谱的具体实现，采用最大熵功率谱，作为其他实施方式，也可以为ar功率谱、或者music功率谱等，本发明对此不做限制。
[0101]
上述实施例中，为了提高特征提取的效率，在某两相邻预分割频点之间的频率差值小于等于分割阈值时，将两个相邻预分割频点的频率平均值对应的频点为最终的分割频点，作为其他实施方式，也可以通过设置不同权值的加权方式确定最终的分割频点，例如最终的分割频点为βω2 (1
‑
β)ω3,β为ω2的权重，0<β<1可依据a1和a2的幅值进行β的确定，若a1>a2,相应的β值大于1
‑
β。
[0102]
上述实施例中，采用svd分解方法提取分段后的特征值，并且为了提高特征值提取的准确性，对每个频段分量的时域表示构造对应的hankel矩阵，并对每个矩阵进行svd分解，取每个矩阵最大的奇异值作为表征该分量的特征值，作为其他实施方式，也可以采用现有的常规的svd分解方式进行特征值的提取。奇异值分解(svd)是一种正交变换，通过将原矩阵转化为一个对角矩阵，得到的奇异值可以有效反映原矩阵的一些特征，目前常规的svd
是将时频分析得到的各个分量按行排列成复合矩阵进行分解，最后得到从大到小顺序排列的奇异值，然而实际中各个分量对应的奇异值并不是从大到小排列的，也就是说得到的奇异值并不能反映对应分量的特征信息，因此采用常规的分解方式可能对各个分量的特征信息产生丢失。或者也可以采用信息熵理论进行特征值的提取，本发明对特征值提取的具体实施方式不做限制。
[0103]
上述实施例中，为了减小无效频段分量的影响，得到各频段分量的时域表示后，还包括利用相关性计算对时域表示的频段进行筛选的步骤，作为其他实施方式，保证所有频段分量为有效频段分量的情况下，也可以不进行筛选。
[0104]
本发明提出的特征提取方法不仅适用于微震信号，还适用于过分割的其他信号，适用范围广。
[0105]
本发明在对功率谱进行划分时，首先将局部极大值点左右相邻的极小值点作为预分割频点，接着判断这些预分割频点是否为最终的分割频点：根据不同的局部极大值点的相邻预分割频点之间的频率差值进行判断，若频率差值小于等于分割阈值，则还需要根据这两个相邻预分割频点确定最终的分割频点，若频率差值大于分割阈值，则这两个相邻预分割频点即为最终的分割频点。本发明避免了过分割而形成无效分量，提高了对微震信号分割的准确性。