一种用于振动信号特征提取的多重超阶分析方法与流程
本发明涉及旋转机械振动信号的处理技术领域，尤其涉及一种用于振动信号特征提取的多重超阶分析方法。
背景技术：
故障诊断方法主要关注系统的运行状态，能够及时发现故障并指导维修，对于提高系统可靠性有重要作用。一般情况下，当旋转部件出现故障时，损伤点在载荷区域和其他元件接触会产生冲击作用，加剧部件的振动，因此振动信号携带了重要的诊断信息，是设备状态识别的重要依据。同时振动在信号上通常表现出非平稳性与多重分形特征。因此从复杂的故障信号中提取能表征故障的特征信息成为旋转机械故障诊断的关键问题之一。
传统的信号处理方法如快速傅里叶变换(fft)等在处理平稳信号时有着很好的效果，但并不适用于非平稳信号。小波变换由于本质上仍然是一种窗口可调傅里叶变换，受限于小波基的长度，缺乏自适应性。为了分析复杂时间序列信息，去趋势波动分析(detrendedfluctuationanalysis,dfa)被提出用来量化非平稳时间的标度指数。然而，dfa方法仅依靠一个2阶波动函数来对时间序列的单个标度进行分析，只适合于处理一维的单重分形时间序列，难以分析具有多重分形特征的多标度时间序列。多重分形去趋势波动分析(multifractaldetrendedfluctuationanalysis，mf-dfa)可以克服dfa的缺陷，有效分析非平稳时间序列的多重分形特征。但是对于类型相近的故障，mf-dfa方法存在提取参数接近，交叉以及不同程度的状态混叠的问题，影响故障特征提取精度。
技术实现要素：
本发明所要解决的技术问题在于提供一种故障分类精确的多重超阶分析方法。
本发明是通过以下技术方案解决上述技术问题的：
一种用于振动信号提取的多重超阶分析方法，包括以下步骤：
步骤1：对待处理的原始时间序列x(t),(t＝1,2,3,…,n)取极值
步骤2：计算原序列的极值增量序列
步骤3：构造步骤1获得的时间序列的去均值序列y(i)；
步骤4：改变去均值序列的时间尺度，将得到的去均值序列y(i)按时间尺度长度s划分为ns＝int(n/s)个独立的区间，int()表示向下取整；从数据的反方向再次以相同的长度分段，共得到2ns个区间；
步骤5：用最小二乘法对每个区间进行高阶多项式趋势拟合；
步骤6：计算每个区间数据的均方误差；
步骤7：对所有区间计算不同分形阶数的波动函数；
步骤8：改变步骤4中的时间尺度，多次重复步骤4-7，获取广义赫斯特指数；
步骤9：计算时间序列的奇异指数与多重分形奇异谱，提取特征点作为信号特征。
优选地，步骤1中的极值满足以下条件：
其中，1≤j≤n-2。
优选地，步骤3构造的去均值序列为：
其中，
优选地，步骤5中拟合得到的高阶多项式趋势为：
yv(i)＝a1ik+a2ik-1+…+aki+ak+1(i＝1,2,…,s；k＝1,2,…)
其中，yv(i)代表第v段数据上第i点所对应的趋势，v＝1,2,…,2ns。
优选地，步骤6中计算均方误差的方法如下：
优选地，步骤7中分形阶数的波动函数如下：
其中，q为分形阶数。
优选地，步骤8中的广义赫斯特指数满足以下条件：
fq(s)∝sh(q)
其中，h(q)为广义赫斯特指数。
优选地，步骤9中计算奇异指数和多重分形奇异谱的方法如下：
f(α)＝q[α-h(q)]+1
其中，α为奇异指数，f(α)为多重分形奇异谱。
优选地，步骤1和2中的极值增量序列与原始时间序列之间的波动特征关系如下：对于原始时间序列xt,(t＝1,2,3,…,n)，其分配函数如下：
zq(l)＝[|xt+l-xt|q]
其中，[·]为期望；q为矩指数即分形阶数；
若原始时间序列存在长相关性，则分配函数满足幂律特征
zq(l)～lτ(q)
式中～表示等价关系，τ(q)为标度指数；
其增量序列δxt的两点相关函数c(l)如下：
c(l)＝[δxt·δxt+l]
若原始时间序列为长程相关的平稳高斯时间序列，则c(l)满足：
c(l)～l-γ
h(q)＝(1-γ)/2
其中，0＜γ＜1为相关指数。
优选地，步骤9中的多重分形奇异谱的计算方法如下：
广义赫斯特指数h(q)与标度指数τ(q)存在以下关系：
τ(q)＝qh(q)-1
经过legendre变换，得到
α＝dτ(q)/dq
f(α)＝qα-τ(q)
在多重分形奇异谱上提取作为信号分类依据的特征点分别为：
左端点：(αmin,f(αmin))
右端点：(αmin,f(αmax))
极值点：(α0,fmax)
其中，fmax＝f(α0),α0∈[αmin,αmax]。
本发明提供的用于振动信号特征提取的多重超阶分析方法的优点在于：提取了原始信号的极值增量序列，突出了原有信号中的波动与冲击特征，更有利于旋转机械的故障特征提取；结合了多重去趋势波动分析，可以滤除非平稳趋势对时间序列的影响，消除或削弱信号中的噪声信号的分形特征，有效分析非平稳时间序列的多重分形特征。
附图说明
图1是本发明的实施例所提供的用于振动信号特征提取的多重超阶分析方法的流程图；
图2是本发明的实施例所提供的用于振动信号特征提取的多重超阶分析方法获取的广义赫斯特指数曲线图；
图3是本发明的实施例所提供的用于振动信号特征提取的多重超阶分析方法获取的多重分形奇异谱示意图；
图4是本发明的实施例所提供的用于振动信号特征提取的多重超阶分析方法与mf-dfa方法分类结果的对比图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。
本实施例选用美国凯斯西储大学电气工程实验室的轴承数据中心获取的轴承数据作为本发明方法分析的振动信号，试验轴承选用skf6205-2rsjemskf深沟球轴承，工作转速为1750rpm，采样频率为12khz；为分析本发明提出方法对于不同故障状态振动信号的特征提取能力，分别选取正常、内圈故障、外圈故障和滚子故障4种故障类型，故障点蚀直径均为0.18mm，另外选取点蚀直径为0.54mm的内圈故障来检验本发明提供的多重超阶分析方法对于故障损伤程度的识别效果，一共分为正常、滚子故障、内圈轻度故障、内圈重度故障、外圈故障5种信号类型，每种状态信号取20组。
如图1所示，一种用于振动信号特征提取的多重超阶分析方法，包括以下步骤：
步骤1：对待处理的原始时间序列x(t),(t＝1,2,3,…,n)取极值
其中，1≤j≤n-2。
步骤2：计算原序列的极值增量序列
通过提取原始时间序列的极值增量序列，突出了原始时间序列中的波动和冲击特征。
步骤3：构造步骤1获得的时间序列的去均值序列y(i)，
其中，
步骤4：改变去均值序列的时间尺度，将得到的去均值序列y(i)按时间尺度长度s划分为ns＝int(n/s)个独立的区间，int()表示向下取整；从数据的反方向再次以相同的长度分段，共得到2ns个区间；本实施例中时间尺度优选为
步骤5：用最小二乘法对每个区间进行高阶多项式趋势拟合，得到
yv(i)＝a1ik+a2ik-1+…+aki+ak+1(i＝1,2,…,s；k＝1,2,…)
其中，yv(i)代表第v段数据上第i点所对应的趋势，v＝1,2,…,2ns。
步骤6：计算每个区间数据的均方误差；
步骤7：对2ns个区间计算不同分形阶数的波动函数；
其中，q为分形阶数，优选实施例中q∈[-5,5]。
步骤8：改变步骤4中的时间尺度，多次重复步骤4-7，获取广义赫斯特指数，所述广义赫斯特指数满足以下条件：
fq(s)∝sh(q)
其中，h(q)为广义赫斯特指数。
步骤9：计算时间序列的奇异指数与多重分形奇异谱：
f(α)＝q[α-h(q)]+1
其中，α为奇异指数，f(α)为多重分形奇异谱，然后在多重分形奇异谱上提取特征点作为信号特征。
步骤1和2中的极值增量序列与原始时间序列之间的波动特征关系如下：
对于原始时间序列xt,(t＝1,2,3,…,n)，其分配函数如下：
zq(l)＝[|xt+l-xt|q]
其中，[·]为期望；q为矩指数即分形阶数；
若原始时间序列存在长相关性，则分配函数满足幂律特征
zq(l)～lτ(q)
式中～表示等价关系，τ(q)为标度指数；
其增量序列δxt的两点相关函数c(l)如下：
c(l)＝[δxt·δxt+l]
若原始时间序列为长程相关的平稳高斯时间序列，则c(l)满足：
c(l)～l-γ
h(q)＝(1-γ)/2
其中，0＜γ＜1为相关指数。
步骤9中的多重分形奇异谱的计算方法如下：
广义赫斯特指数h(q)与标度指数τ(q)存在以下关系：
τ(q)＝qh(q)-1
经过legendre变换，得到
α＝dτ(q)/dq
f(α)＝qα-τ(q)
本实施例中得到的广义赫斯特指数和多重分形奇异谱分别如图2和图3所示；在图3中选取以下特征点作为作为信号分类依据：
左端点：(αmin,f(αmin))
右端点：(αmin,f(αmax))
极值点：(α0,fmax)
其中，fmax＝f(α0),α0∈[αmin,αmax]。
如图4所示，本实施例提供的方法的分类结果与mf-dfa的结果对比可以发现，本发明提供的方法对于信号的不规则程度更加敏感，可有效突出信号中的波动与冲击特征，从而得到更准确的分类结果。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，在不脱离本发明的精神和原则的前提下，本领域普通技术人员对本发明所做的任何修改、等同替换、改进等，均应落入本发明权利要求书确定的保护范围之内。
- 还没有人留言评论。精彩留言会获得点赞！