CN113537102B - 一种微震信号的特征提取方法 - Google Patents

一种微震信号的特征提取方法 Download PDF

Info

Publication number
CN113537102B
CN113537102B CN202110832860.1A CN202110832860A CN113537102B CN 113537102 B CN113537102 B CN 113537102B CN 202110832860 A CN202110832860 A CN 202110832860A CN 113537102 B CN113537102 B CN 113537102B
Authority
CN
China
Prior art keywords
frequency
segmentation
points
frequency band
microseismic
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.)
Active
Application number
CN202110832860.1A
Other languages
English (en)
Other versions
CN113537102A (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.)
Shenzhen Smart Microelectronics Technology Co ltd
Original Assignee
Spl Electronic Technology Co ltd
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 Spl Electronic Technology Co ltd filed Critical Spl Electronic Technology Co ltd
Priority to CN202110832860.1A priority Critical patent/CN113537102B/zh
Publication of CN113537102A publication Critical patent/CN113537102A/zh
Application granted granted Critical
Publication of CN113537102B publication Critical patent/CN113537102B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. analysis, for interpretation, for correction
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/21Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02Preprocessing
    • G06F2218/04Denoising
    • G06F2218/06Denoising by applying a scale-space analysis, e.g. using wavelet analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08Feature extraction
    • G06F2218/10Feature extraction by analysing the shape of a waveform, e.g. extracting parameters relating to peaks
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30Assessment of water resources

Abstract

本发明涉及一种微震信号的特征提取方法,属于信号处理技术领域。方法包括:获取去趋势项后的微震信号;计算去趋势项后的微震信号的功率谱,找出功率谱中的局部极大值点;将局部极大值点左右相邻的极小值点作为预分割频点,计算属于不同的局部极大值点的相邻预分割频点之间的频率差值;判断频率差值的大小,若某两相邻预分割频点之间的频率差值小于等于分割阈值,则根据这两个相邻预分割频点确定最终的分割频点,若某两相邻预分割频点之间的频率差值大于分割阈值,则这两个相邻预分割频点即为最终的分割频点;按照最终的分割频点对功率谱进行频段的划分,进而提取特征值。本发明避免了过分割而形成无效分量,提高了对微震信号分割的准确性。

Description

一种微震信号的特征提取方法
技术领域
本发明涉及一种微震信号的特征提取方法,属于信号处理技术领域。
背景技术
微震监测系统通过分析生产活动所引起的微震事件来监测岩体的稳定状态。但矿山环境复杂,干扰因素众多,将有效的微震事件从大量的数据样本中识别出来是发挥微震监测系统预测预警的基础性问题,而特征提取是对微震信号进行准确分析和精确识别的关键步骤,具有重要的研究意义。
目前常用的特征提取流程是通过时频分析方法将信号分解成多个特征分量,再利用SVD或者信息熵理论得到特征分量的差异性特征值。典型的时频分析方法主要有傅里叶变换(FT),经验模态分解(EMD)、小波变换(WT)和小波包变换(WPT)。这些方法在一定程度上得到了信号的特征分量,但还存在一些不足。FT适用于平稳信号的分析,但面对具有非平稳和突发瞬态的微震信号时难以进行有效分解;WPT和WT虽然适用于非平稳信号处理,但对信号的分解效果很大程度依赖于分解尺度和基函数的选择,自适应性不够好;EMD自适应性好,但存在模态混叠、端点效应的问题,而且EMD缺乏严格理论支撑。
为此,2013年Gills提出一种新的时频分析方法——经验小波变换(EWT),它是在WT理论框架下提出的,所以理论基础完备,并且借助于小波分解快速算法使得自身具有较高的计算效率。EWT首先利用傅里叶变换得到信号的频域表示,再利用局部极大值的中点作为分界频点对信号频谱进行划分,根据不同的频段构建经验小波,并获得多个窄带频率分量,具有EMD的自适应性。近年来EWT在不同领域的信号分析中得到应用,如电机振动信号,心电信号、变压器振动信号等。但随着研究的深入,EWT存在的一些问题也在显现,EWT方法的有效性分解主要取决于频谱的有效分割,但微震信号频域特性复杂,传统的依靠局部极大值的频谱分割方法会将同一频带的成分划分到不同的频带,使得同一分量被划分成两个分量,造成过分割的现象,另外,传统的EWT方法在分析非平稳、局部噪声信号时,非平稳分量以及噪声分量产生的一些局部极大值可能出现并错误的保持在峰值序列中,而有效信号的极大值可能不保持在峰值序列中,导致不当的分割,进而使得微震信号的特征提取准确性低。
发明内容
本申请的目的在于提供一种微震信号的特征提取方法,用以解决现有特征提取方法准确性的问题。
为实现上述目的,本申请提出了一种微震信号的特征提取方法的技术方案,包括以下步骤:
1)获取待特征提取的微震信号,对微震信号进行去趋势项处理,得到去趋势项后的微震信号;
2)计算去趋势项后的微震信号的功率谱,找出功率谱中的局部极大值点;
3)将局部极大值点左右相邻的极小值点作为预分割频点,计算属于不同的局部极大值点的相邻预分割频点之间的频率差值;
4)判断频率差值的大小,若某两相邻预分割频点之间的频率差值小于等于分割阈值,则根据这两个相邻预分割频点确定最终的分割频点,若某两相邻预分割频点之间的频率差值大于分割阈值,则这两个相邻预分割频点即为最终的分割频点;
5)按照最终的分割频点对功率谱进行频段的划分,进而提取特征值。
本发明的微震信号的特征提取方法的技术方案的有益效果是:本发明在对功率谱进行划分时,首先将局部极大值点左右相邻的极小值点作为预分割频点,接着判断这些预分割频点是否为最终的分割频点:根据不同的局部极大值点的相邻预分割频点之间的频率差值进行判断,若频率差值小于等于分割阈值,则还需要根据这两个相邻预分割频点确定最终的分割频点,若频率差值大于分割阈值,则这两个相邻预分割频点即为最终的分割频点。本发明通过极大极小值频点进行初步频谱分割,进而利用根据采样频率自适应调控的阈值对初步划分的频谱进一步进行频点优化,避免了过分割或者不当分割而形成无效分量,提高了对微震信号分割的准确性。
进一步的,为了减小功率谱的畸变,使得功率谱更加接近与真实谱,所述步骤2)中,计算的功率谱为最大熵功率谱。
进一步的,为了提高分割频点的计算效率,所述步骤4)中,若某两相邻预分割频点之间的频率差值小于等于分割阈值,则这两个相邻预分割频点的频率平均值对应的频点为最终的分割频点。
进一步的,所述步骤5)中,提取各频段的特征值的步骤包括:
构建各频段内的小波函数和尺度函数,利用傅里叶反变换得到各频段分量的时域表示;
根据各频段分量的时域表示构造矩阵;
对矩阵进行SVD分解,进而提取特征值。
进一步的,为了减小无效频段分量的影响,得到各频段分量的时域表示后,还包括利用相关性计算对时域表示的频段进行筛选的步骤:计算各频段分量的时域表示与去趋势项后的微震信号的相关系数,若相关系数小于相关系数阈值,则剔除该频段分量的时域表示。
进一步的,为了充分提取每个频段分量的特征信息,更加准确的得到每个频段分量的特征值,所构造的矩阵包括每个频段分量对应的Hankel矩阵,对每个Hankel矩阵进行SVD分解,得到的每个Hankel矩阵的最大奇异值即为对应频段分量的特征值。
进一步的,为了更好的进行去趋势处理,步骤1)中利用最小二乘法进行去趋势项处理,趋势项为:
Figure BDA0003176163450000031
其中,Dt为趋势项;fs为微震信号的采样率;t为采样点序号;ai(i=0,1,2,...,K)为待定系数;K为趋势项的阶数。
进一步的,最大熵功率谱的计算过程为:
Figure BDA0003176163450000032
其中,S(ω)为最大熵功率;ω为频率;
Figure BDA0003176163450000033
为第m阶的预报误差方差估计;/>
Figure BDA0003176163450000034
为第m阶第j次的自回归系数;fs为微震信号的采样率;/>
Figure BDA0003176163450000035
为虚数符号。
进一步的,通过Burg递推算法计算预报误差方差估计和自回归系数。
进一步的,相关系数的计算过程为:
Figure BDA0003176163450000036
其中:co(k)为第k个频段分量的相关系数;x'(t)为去趋势项后的微震信号,
Figure BDA0003176163450000037
为x'(t)的平均值;xk(t)为第k个频段分量的时域表示;/>
Figure BDA0003176163450000038
为各频段分量的时域表示的平均值;N为频段分量的数量;t为采样点序号;M为微震信号的采样点数量。
附图说明
图1是本发明微震信号的特征提取方法的流程图;
图2是本发明小波函数和尺度函数的滤波器组图。
具体实施方式
微震信号的特征提取方法实施例:
本发明的主要构思在于,基于微震信号的特点,为了避免过分割的现象,本发明在得到最大熵功率谱中的局部极大值点后,以局部极大值点左右相邻的极小值点作为预分割频点,对预分割频点进行判断:计算属于不同的局部极大值点的相邻预分割频点之间的距离,若距离小于等于分割阈值,则根据这两个相邻预分割频点再次确定最终的分割频点,若距离大于分割阈值,则这两个相邻预分割频点即为最终的分割频点,根据最终的分割频点对最大熵功率谱进行频段划分,进而得到各频段的特征值。
具体的,微震信号的特征提取方法如图1所示,包括以下步骤:
1)获取检测到的微震信号,对微震信号进行模数转换后进行M点截断,得到M点采样值序列的离散信号x(t),(t=1,2,...M),将离散信号x(t)作为待特征提取的微震信号。
2)利用最小二乘法对步骤1)的得到的离散信号x(t)进行去趋势项处理,得到去趋势项的离散信号x'(t)。
趋势项是由于传感器的零点漂移,基础运动等引起的信号波形偏移,本步骤中,利用最小二乘法进行去趋势项处理的过程如下:
利用K阶多项式(也即估计值)来拟合离散信号x(t)中的趋势项Dt
Figure BDA0003176163450000041
其中,fs为微震信号的采样率;ai(i=0,1,2,...,K)为待定系数;t为采样点序号;K为趋势项的阶数。待定系数ai的确定依据是使得K阶多项式Dt与x(t)之间的误差平方和最小,有如下公式:
Figure BDA0003176163450000042
令E对ai的偏导数为零,可产生K+1元线性方程组,第j元线性方程为:
Figure BDA0003176163450000043
通过求解(3)式,可得到K+1个待定系数ai,将求出的待定系数ai代入估计值Dt中,即可得到估计值Dt
去趋势项的离散信号x'(t)为:x'(t)=x(t)-Dt (4)
当然作为其他实施方式,还可以采用插值法进行去趋势项处理,但插值法需严格的通过数据,当数据带有误差时,就会使预测值偏离真实数据,产生过拟合现象。
3)根据步骤2)中得到的去趋势项的离散信号x'(t),结合伯格(Burg)算法计算出自回归系数
Figure BDA0003176163450000051
和预报误差方差估计值/>
Figure BDA0003176163450000052
进而将两个参数带入最大熵功率谱S(ω)的计算公式中计算出最大熵功率谱S(ω),并将频率轴归一化到[0,π]。
本步骤中,最大熵功率谱S(ω)的计算公式如下:
Figure BDA0003176163450000053
其中,
Figure BDA0003176163450000054
为第m阶的预报误差方差估计;/>
Figure BDA0003176163450000055
为第m阶第j次的自回归系数;/>
Figure BDA0003176163450000056
为虚数符号。参数/>
Figure BDA0003176163450000057
和/>
Figure BDA0003176163450000058
的确定需要利用Burg递推算法,计算过程如下:
(1)计算预报误差方差估计的初始值
Figure BDA0003176163450000059
Figure BDA00031761634500000510
Figure BDA00031761634500000511
式中,
Figure BDA00031761634500000512
为前向预测误差的初始值;/>
Figure BDA00031761634500000513
为后向预测误差的初始值;
(2)依据如下公式计算第m阶的反射系数Km
Figure BDA00031761634500000514
式中,
Figure BDA00031761634500000515
为第m阶的前向预测误差;/>
Figure BDA00031761634500000516
为第m阶的后向预测误差;p是总阶数量;
(3)计算第m阶第j次的自回归系数
Figure BDA00031761634500000517
Figure BDA00031761634500000518
(4)计算第m阶的预报误差方差估计
Figure BDA00031761634500000519
Figure BDA00031761634500000520
(5)计算预测滤波器输出:
Figure BDA00031761634500000521
Figure BDA00031761634500000522
式中,
Figure BDA00031761634500000523
为Km的复共轭。
(6)令m=m+1,重复步骤(2)~(5),直到m=p。
将求得的参数代入公式(5)中可求得最大熵功率谱。
当然,作为其他实施方式,还存在自相关法和改进协方差法算出自回归系数和预报误差方差估计值,但改进协方差法计算繁琐,编程较困难,自相关法计算最简单,谱估计的分辨率相对较差,Burg法是较为通用的方法,计算不太复杂,谱估计质量较好。
4)找出功率谱中的局部极大值点,将局部极大值点左右相邻的极小值点作为预分割频点,计算属于不同的局部极大值点的相邻预分割频点之间的距离d(这里的距离是指相邻预分割频点之间的频率差值),若距离d小于等于分割阈值α,则根据这两个相邻预分割频点确定最终的分割频点,若距离d大于分割阈值α,则这两个相邻预分割频点即为最终的分割频点,根据最终的分割频点对最大熵功率谱S(ω)进行频段划分。
本步骤中的局部极大值点的确定过程为:假设需要分割的频段数为N个,用ωn(n=1,2,...,N+1)表示频段分割点。找出最大熵功率谱S(ω)中所有的局部极大值,从大到小排列,保留前N个局部极大值对应的频段分割点为局部极大值点。
本步骤中,若属于不同的局部极大值点的相邻预分割频点之间的距离d小于分割阈值α,则将这两个相邻预分割频点的频率平均值作为最终的分割频点。例如:相邻局部极大值点分别为A1和A2,A1位于A2的左侧,也即A1<A2;局部极大值点A1左右相邻的极小值点分别为ω1,ω2,局部极大值点A2左右相邻的极小值点分别为ω3,ω4,也即ω1,ω2,ω3,ω4均为预分割频点;
属于A1和A2的相邻预分割频点为ω2,ω3,计算ω2,ω3之间的距离d,若d≤α(经过多次实验验证,α一般取10/fs效果较好),则剔除ω2,ω3这两个预分割频点,将
Figure BDA0003176163450000061
作为最终的分割频点;若d>α,则确定ω2,ω3为最终的分割频点。
5)在划分的各个频段内构建对应的小波函数ψk(ω)和尺度函数
Figure BDA0003176163450000062
这里的k表示第k个频段,k=0,1,2,…,N-1。
本步骤中,第k个频段的小波函数ψk(ω)和尺度函数
Figure BDA0003176163450000063
的频率响应如图2所示,计算公式如下:
Figure BDA0003176163450000064
Figure BDA0003176163450000071
其中,ω为归一化频率,ωn∈[0,π]为分割频点,tn=γωn(0<γ<1)为经验小波滤波器过渡带宽度,其中
Figure BDA0003176163450000072
β(x)=x4(35-85x+70x2-20x3),这里的x代表自变量。
6)根据各个频段内的小波函数和尺度函数,利用傅里叶反变换得到尺度系数和细节系数,进而得到各频段分量的时域表示。
应用傅里叶反变换计算S(ω)×ψk(ω)和
Figure BDA00031761634500000713
从而得到各个分量x0(t),x1(t),...xN-1(t)的时域表示。
具体的,本步骤中通过信号x'(t)与小波函数和尺度函数的内积运算获取信号分解的尺度系数
Figure BDA0003176163450000073
k=0,也即/>
Figure BDA0003176163450000074
和细节系数/>
Figure BDA0003176163450000075
计算公式如下:
Figure BDA0003176163450000076
其中,
Figure BDA0003176163450000077
和/>
Figure BDA0003176163450000078
分别为经验尺度函数/>
Figure BDA0003176163450000079
和经验小波函数ψk(t)的复共轭;
Figure BDA00031761634500000710
和ψk(ω)分别为/>
Figure BDA00031761634500000711
和ψk(t)的Fourier变换;F-1[·]为Fourier反变换(即傅里叶反变换)。
至此,可计算各个频段分量的时域表示:
Figure BDA00031761634500000712
式中:“*”代表卷积。
7)利用相关性分析对x0(t),x1(t),...xN-1(t)共N个频段分量的时域表示进行筛选。
本步骤中,对得到的N个频段分量的时域表示进行筛选,筛选是依据各个频段分量与去趋势项后的微震信号相关系数值的大小,当频段分量的时域表示与去趋势项后的微震信号的相关系数值小于0.03(相关系数阈值)时,即可认为该频段分量为无效分量,进行剔除。相关性计算公式为:
Figure BDA0003176163450000081
其中:co(k)为第k个频段分量的相关系数;x'(t)为去趋势项后的微震信号,
Figure BDA0003176163450000082
为x'(t)的平均值;xk(t)为第k个频段分量的时域表示;/>
Figure BDA0003176163450000083
为各频段分量的时域表示的平均值;N为频段分量的数量;t为采样点序号;M为微震信号的采样点数量。
8)利用筛选后的每个频段分量的时域表示构造对应的Hankel矩阵,并对每个矩阵进行SVD分解,取每个矩阵最大的奇异值作为表征该分量的特征值。
本步骤中,Hankel矩阵Hα的构造形式如下:
Figure BDA0003176163450000084
式中:xk,n+1为第k个频段分量第n+1个采样点的时域表示;且1<n<M;M是截断的采样点长度,其中矩阵Hα行列数按以下规则选取:
(1)若M为偶数,Hankel矩阵的行数为m=M/2+1,列数为n=M/2;
(2)若M为奇数,Hankel矩阵的行数和列数都为m=(M+1)/2;
分别对生成的m×n阶矩阵Hα进行SVD,SVD分解形式如下:
Figure BDA0003176163450000085
式中:r(r≤min(m,n))为矩阵Hα的秩,U∈Rm×m和V∈Rn×n为正交矩阵,ui和vi为矩阵U和V的第i列向量,其中矩阵
Figure BDA0003176163450000086
并且σi为第i列向量的奇异值;Hi为第i列向量对应的矩阵,σ=diag(σ12,...,σp)(p=min(m,n)),将最大奇异值σ作为表征对应矩阵Hα的特征值。
上述实施例中,关于功率谱的具体实现,采用最大熵功率谱,作为其他实施方式,也可以为AR功率谱、或者music功率谱等,本发明对此不做限制。
上述实施例中,为了提高特征提取的效率,在某两相邻预分割频点之间的频率差值小于等于分割阈值时,将两个相邻预分割频点的频率平均值对应的频点为最终的分割频点,作为其他实施方式,也可以通过设置不同权值的加权方式确定最终的分割频点,例如最终的分割频点为βω2+(1-β)ω3,β为ω2的权重,0<β<1可依据A1和A2的幅值进行β的确定,若A1>A2,相应的β值大于1-β。
上述实施例中,采用SVD分解方法提取分段后的特征值,并且为了提高特征值提取的准确性,对每个频段分量的时域表示构造对应的Hankel矩阵,并对每个矩阵进行SVD分解,取每个矩阵最大的奇异值作为表征该分量的特征值,作为其他实施方式,也可以采用现有的常规的SVD分解方式进行特征值的提取。奇异值分解(SVD)是一种正交变换,通过将原矩阵转化为一个对角矩阵,得到的奇异值可以有效反映原矩阵的一些特征,目前常规的SVD是将时频分析得到的各个分量按行排列成复合矩阵进行分解,最后得到从大到小顺序排列的奇异值,然而实际中各个分量对应的奇异值并不是从大到小排列的,也就是说得到的奇异值并不能反映对应分量的特征信息,因此采用常规的分解方式可能对各个分量的特征信息产生丢失。或者也可以采用信息熵理论进行特征值的提取,本发明对特征值提取的具体实施方式不做限制。
上述实施例中,为了减小无效频段分量的影响,得到各频段分量的时域表示后,还包括利用相关性计算对时域表示的频段进行筛选的步骤,作为其他实施方式,保证所有频段分量为有效频段分量的情况下,也可以不进行筛选。
本发明提出的特征提取方法不仅适用于微震信号,还适用于过分割的其他信号,适用范围广。
本发明在对功率谱进行划分时,首先将局部极大值点左右相邻的极小值点作为预分割频点,接着判断这些预分割频点是否为最终的分割频点:根据不同的局部极大值点的相邻预分割频点之间的频率差值进行判断,若频率差值小于等于分割阈值,则还需要根据这两个相邻预分割频点确定最终的分割频点,若频率差值大于分割阈值,则这两个相邻预分割频点即为最终的分割频点。本发明避免了过分割而形成无效分量,提高了对微震信号分割的准确性。

Claims (10)

1.一种微震信号的特征提取方法,其特征在于,包括以下步骤:
1)获取待特征提取的微震信号,对微震信号进行去趋势项处理,得到去趋势项后的微震信号;
2)计算去趋势项后的微震信号的功率谱,找出功率谱中的局部极大值点;局部极大值点的确定过程为:需要分割的频段数为N个,找出功率谱中所有的局部极大值从大到小排列,保留前N个局部极大值对应的频段分割点为局部极大值点;
3)将局部极大值点左右相邻的极小值点作为预分割频点,计算属于不同的局部极大值点的相邻预分割频点之间的频率差值;
4)判断频率差值的大小,若某两相邻预分割频点之间的频率差值小于等于分割阈值,则根据这两个相邻预分割频点确定最终的分割频点,若某两相邻预分割频点之间的频率差值大于分割阈值,则这两个相邻预分割频点即为最终的分割频点;
5)按照最终的分割频点对功率谱进行频段的划分,进而提取特征值。
2.根据权利要求1所述的微震信号的特征提取方法,其特征在于,所述步骤2)中,计算的功率谱为最大熵功率谱。
3.根据权利要求1所述的微震信号的特征提取方法,其特征在于,所述步骤4)中,若某两相邻预分割频点之间的频率差值小于等于分割阈值,则这两个相邻预分割频点的频率平均值对应的频点为最终的分割频点。
4.根据权利要求1所述的微震信号的特征提取方法,其特征在于,所述步骤5)中,提取各频段的特征值的步骤包括:
构建各频段内的小波函数和尺度函数,利用傅里叶反变换得到各频段分量的时域表示;
根据各频段分量的时域表示构造矩阵;
对矩阵进行SVD分解,进而提取特征值。
5.根据权利要求4所述的微震信号的特征提取方法,其特征在于,得到各频段分量的时域表示后,还包括利用相关性计算对时域表示的频段进行筛选的步骤:计算各频段分量的时域表示与去趋势项后的微震信号的相关系数,若相关系数小于相关系数阈值,则剔除该频段分量的时域表示。
6.根据权利要求4所述的微震信号的特征提取方法,其特征在于,所构造的矩阵包括每个频段分量对应的Hankel矩阵,对每个Hankel矩阵进行SVD分解,得到的每个Hankel矩阵的最大奇异值即为对应频段分量的特征值。
7.根据权利要求1所述的微震信号的特征提取方法,其特征在于,步骤1)中利用最小二乘法进行去趋势项处理,趋势项为:
Figure QLYQS_1
其中,Dt为趋势项;fs为微震信号的采样率;t为采样点序号时刻;ai为待定系数,i=0,1,2,…,K;K为趋势项的阶数。
8.根据权利要求2所述的微震信号的特征提取方法,其特征在于,最大熵功率谱的计算过程为:
Figure QLYQS_2
其中,S(ω)为最大熵功率;ω为频率;
Figure QLYQS_3
为第m阶的预报误差方差估计;/>
Figure QLYQS_4
为第m阶第j次的自回归系数;fs为微震信号的采样率;/>
Figure QLYQS_5
为虚数符号。
9.根据权利要求8所述的微震信号的特征提取方法,其特征在于,通过Burg递推算法计算预报误差方差估计和自回归系数。
10.根据权利要求5所述的微震信号的特征提取方法,其特征在于,相关系数的计算过程为:
Figure QLYQS_6
其中:co(k)为第k个频段分量的相关系数;x'(t)为去趋势项后的微震信号,
Figure QLYQS_7
为x'(t)的平均值;xk(t)为第k个频段分量的时域表示;/>
Figure QLYQS_8
为各频段分量的时域表示的平均值;k=0,2,…,N-1;N为频段分量的数量;t为采样点序号;M为微震信号的采样点数量。
CN202110832860.1A 2021-07-22 2021-07-22 一种微震信号的特征提取方法 Active CN113537102B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110832860.1A CN113537102B (zh) 2021-07-22 2021-07-22 一种微震信号的特征提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110832860.1A CN113537102B (zh) 2021-07-22 2021-07-22 一种微震信号的特征提取方法

Publications (2)

Publication Number Publication Date
CN113537102A CN113537102A (zh) 2021-10-22
CN113537102B true CN113537102B (zh) 2023-07-07

Family

ID=78120537

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110832860.1A Active CN113537102B (zh) 2021-07-22 2021-07-22 一种微震信号的特征提取方法

Country Status (1)

Country Link
CN (1) CN113537102B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114265113B (zh) * 2021-12-23 2022-08-12 中国矿业大学(北京) 微地震事件识别方法、装置和电子设备
CN114383959B (zh) * 2022-01-13 2024-02-23 安徽惠洲地质安全研究院股份有限公司 用于岩体性质评定的随钻检测装置及方法
CN116203646B (zh) * 2023-05-04 2023-07-14 山东省地质测绘院 一种确定地质资源量的勘探数据处理系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104089624A (zh) * 2014-07-18 2014-10-08 赵佳 计步方法及终端设备
CN109765055A (zh) * 2019-01-31 2019-05-17 杭州安脉盛智能技术有限公司 基于ewt、谱有效值和knn的滚动轴承故障检测方法及系统
CN111275244A (zh) * 2020-01-10 2020-06-12 南京航空航天大学 一种车速时间序列分频预测方法
CN111652031A (zh) * 2019-12-02 2020-09-11 昆明理工大学 一种基于改进经验小波变换的滚动轴承故障诊断方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE102004049457B3 (de) * 2004-10-11 2006-07-06 Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. Verfahren und Vorrichtung zur Extraktion einer einem Audiosignal zu Grunde liegenden Melodie
CN108375472A (zh) * 2018-02-12 2018-08-07 武汉科技大学 基于改进经验小波变换的轴承故障诊断方法及系统装置
CN109211568B (zh) * 2018-09-19 2019-11-15 四川大学 基于条件经验小波变换的滚动轴承故障诊断方法
CN109525215B (zh) * 2018-09-29 2023-02-28 长安大学 一种采用峭度谱确定子频带边界的经验小波变换方法
CN109799090B (zh) * 2019-01-08 2020-09-18 长安大学 采用频带3分区的经验小波变换的轴承特征频率提取方法
CN110671613B (zh) * 2019-10-15 2021-03-16 重庆邮电大学 基于改进经验小波变换的流体管道泄漏信号时延估计方法
CN111623982B (zh) * 2020-06-15 2021-03-26 大连理工大学 一种基于apewt和imomeda的行星齿轮箱早期故障诊断方法
CN112411324B (zh) * 2020-09-30 2022-03-29 武汉光谷卓越科技股份有限公司 一种线结构光路面跳车检测方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104089624A (zh) * 2014-07-18 2014-10-08 赵佳 计步方法及终端设备
CN109765055A (zh) * 2019-01-31 2019-05-17 杭州安脉盛智能技术有限公司 基于ewt、谱有效值和knn的滚动轴承故障检测方法及系统
CN111652031A (zh) * 2019-12-02 2020-09-11 昆明理工大学 一种基于改进经验小波变换的滚动轴承故障诊断方法
CN111275244A (zh) * 2020-01-10 2020-06-12 南京航空航天大学 一种车速时间序列分频预测方法

Also Published As

Publication number Publication date
CN113537102A (zh) 2021-10-22

Similar Documents

Publication Publication Date Title
CN113537102B (zh) 一种微震信号的特征提取方法
Kanjilal et al. On multiple pattern extraction using singular value decomposition
CN108845250B (zh) 基于振动信号特征提取的有载分接开关故障识别方法
CN109633431B (zh) 基于振动信号特征提取的有载分接开关故障识别方法
Wang et al. Blind source extraction of acoustic emission signals for rail cracks based on ensemble empirical mode decomposition and constrained independent component analysis
CN108303740B (zh) 一种航空电磁数据噪声压制方法及装置
CN108645620B (zh) 一种基于信息熵和多尺度形态学的滚动轴承早期故障诊断方法
CN111985426A (zh) 一种基于变分模态分解的海杂波混合去噪算法
Barat et al. INTELLIGENT AE SIGNAL FILTERING METHODS.
CN113780055A (zh) 一种momeda与压缩感知的滚动轴承故障诊断方法
CN114330459A (zh) 基于多层分解的振源响应信噪分离方法及系统
CN108108666B (zh) 一种基于小波分析和时频单源检测的混合矩阵估计方法
CN110865375B (zh) 一种水中目标检测方法
CN112817056B (zh) 一种大地电磁信号去噪方法及系统
CN113657268B (zh) 一种应用于风电机组齿轮箱故障诊断的信号自动分解方法
CN112285793B (zh) 一种大地电磁去噪方法及系统
CN115310477A (zh) 基于分形特征和捕食者算法的泵机设备故障声音检测方法及其系统
CN115563480A (zh) 基于峭度比系数筛选辛几何模态分解的齿轮故障辨识方法
CN110688981B (zh) 一种振动信号去噪的模态混叠消除方法
CN112082793A (zh) 一种基于SCA和FastICA的旋转机械耦合故障诊断方法
CN112082792A (zh) 一种基于mf-jade的旋转机械故障诊断方法
CN110765881B (zh) 一种基于主成分分析的小波基选择方法
CN117007313A (zh) 基于多尺度符号动力熵高密度小波的机械故障诊断方法
CN114397010A (zh) 基于小波分解的瞬态信号声成像方法
CN112329667A (zh) 基于经验小波变换和多尺度熵的水质数据清洗方法

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
CP01 Change in the name or title of a patent holder

Address after: Room 4005, block a, block 8, area C, Wanke Yuncheng phase III, Liuxin 4th Street, Xili community, Xili street, Nanshan District, Shenzhen, Guangdong 518000

Patentee after: Shenzhen Smart Microelectronics Technology Co.,Ltd.

Address before: Room 4005, block a, block 8, area C, Wanke Yuncheng phase III, Liuxin 4th Street, Xili community, Xili street, Nanshan District, Shenzhen, Guangdong 518000

Patentee before: SPL ELECTRONIC TECHNOLOGY CO.,LTD.

CP01 Change in the name or title of a patent holder