基于改进变分模态分解的滚动轴承故障诊断方法

2020-06-13孟宗,吕蒙,殷娜,李晶

计量学报 2020年6期

孟宗, 吕蒙, 殷娜, 李晶

(燕山大学河北省测试计量技术及仪器重点实验室, 河北秦皇岛 066004)

1 引言

滚动轴承在各种旋转机械中应用广泛,其运行状态往往影响整台机器的性能,对滚动轴承进行状态监测以及故障诊断具有重要意义。当轴承发生损伤时,其振动信号表现为非线性、非平稳、非高斯特征,且信号微弱、调制性强、背景噪声大,使故障特征淹没在强背景噪声中,难以提取[1,2]。目前常见的非线性非平稳处理方法有:经验模态分解(empirical mode decomposition,EMD)、局部均值分解(local mean decomposition,LMD)等。以上方法在轴承故障信号处理中都得到了良好的效果[3],但也有各自的局限性。EMD分解缺乏严格的数学基础,容易产生模态混叠、端点效应[4],并且存在算法效率低等问题;LMD分解存在迭代计算量大、由解调引起的信号突变等问题[5]。变分模态分解(variational mode decomposition,VMD)是最近提出的一种信号分解方法,该方法具有较好的抗噪性和更高的频域分辨力,能够提取出较微弱的信号成分[6]。然而VMD分解必须提前设定分解层数,层数设置不当会引起模态混叠或者大量虚假分量的产生[7]。文献[8]使用蛙跳算法搜索VMD算法的最优参数,文献[9]通过粒子群算法搜索VMD分解的最佳参数,以上方法都以包络熵作为判断标准。这些方法不但计算繁琐,其效果在很大程度上受到优化算法的参数设定的影响,并且没有足够的理论依据。

本文利用改进的奇异值分解(singular value decomposition, SVD)对故障信号降噪[10～12],通过奇异值平均下降速度法确定有效奇异值个数;并通过能量占比确定VMD分解的最佳分解层数,将根据相关系数判断为分解自同一调频部分的分量信号作为一个整体处理,最后对分解得到的IMF分量进行包络解调提取故障特征,判断故障类型。

2 改进变分模态分解方法

2.1 变分模态分解原理

VMD算法通过不断地迭代搜索变分模态最优解来确定每个分量的中心频率和带宽。将原信号非递归地分解为一系列限带宽的固有模态函数,最终得到的分量信号在频域中具有一定的稀疏特性。进行K阶分解时(K是分解层数),产生的约束变分问题为:

(1)

式中:∑kuk=f为约束条件;{uk}={u1,…,uK}和{ωk}={ω1,…,ωK}为所有的分量和中心频率,k=1,…,K。引入增广拉格朗日函数来求解以上约束变分问题[13],即:

L({uk},{ωk},λ)=

(2)

利用交替方向乘数法(ADMM)求上述拉格朗日函数的鞍点,详细步骤如下。

2)令n=n+1,进行K次循环,先对uk更新,随即更新相应的ωk,直到K个分量全部更新完毕;

3)采用对偶上升法更新λ:

(3)

判断是否满足收敛条件:

(4)

若不满足则返回步骤(2),满足则循环停止,得到最终的K个分量。

2.2 能量占比确定分解层数

VMD分解需要提前设定分解层数,理想情况下的分解层数恰好等于信号中包含的窄带分量个数。分解层数与衰减因子α共同作用,对分解效果产生影响[14]:

1)层数太少导致分解不彻底,某些频带内的有效信息包含在其他分量中从而产生模态混叠(α偏小)或作为“噪声”丢弃(α偏大);

2)太多的分解层数会将噪声作为有效分量分离出来(α偏小),或者使得若干分量的中心频率重合(α偏大)。

VMD分解得到的全部分量的能量之和与原信号能量接近程度,称作能量占比,反映了VMD分解效果。能量占比的定义如下:

(5)

式中:uk为第k个分量信号;f为原信号。由于VMD本质上是一种维纳滤波,远离分量中心频率的频段会产生衰减,因此总有能量占比P<1。

能量占比随分解层数的变化情况如下:

1)当分解层数K设置偏小时,无论α的值如何,都会导致能量的较大幅度丢失,此时的能量占比小于理想情况,分解层数越接近理想情况,则能量占比越大;

2)当K设置偏大,α也偏大时,也会导致原信号能量的损失,能量占比再次下降;

3)而当α偏小时,部分噪声会作为有效分量被分解出来,但由于噪声信号在频带上能量分散,因此分解自噪声的分量一般能量很小,此时能量占比相对最佳分解层数时会略微增大。

因此,在分解前对信号进行去噪,并合理设置α的前提下,能量占比首次达到极值时的分解层数就是最佳分解层数。

2.3 相关系数识别调频成分

实际采集到的信号大多为调制信号,存在调幅调频成分,相对于单一频率的信号,调频信号的频谱较宽,由于VMD算法只能保证分量的准正交性,因此VMD分解时有可能将信号中的调频成分分解到多个分量之中[15]。当α较大时,导致IMF分量的频带较宽,更容易识别调频成分,但也容易产生模态混叠;α较小时,模态混叠较少,但调频成分更容易被分解到多个分量之中。

由于实际信号中既有频率接近并互相独立的单一频率成分,又有占频带较宽的调频成分,因此不存在一个α的最优值能够兼顾频域分辨力和保留调频成分的能力。

本文通过将α设置为归一化采样频率的2倍,保留了VMD的频域分辨力,对于调频分量是否被分解到不同分量中,通过中心频率相邻的分量间的相关系数进行判断并作出修正,相关系数为:

(6)

式中:ui、ui+1为中心频率相邻的2个分量信号,其中Cov(ui,ui+1)为ui与ui+1的协方差,D(ui)为ui的方差。相关系数越大,说明两者在频域上有更多的重叠,ρ>0.1时,认为两分量分解自同一调频部分,相加作为一个分量处理。

3 基于SVD和改进VMD的故障诊断方法

本文使用奇异值分解中的奇异值平均下降速度法对原信号进行降噪处理,对于以N个点组成的采样时间序列f构造的Hankel矩阵Am×n,其中m为最接近N/2的整数[16],必然存在两正交矩阵U∈Rm×m和V∈Rn×n,其中Rm×m表示m×m维的欧氏空间,使得:

A=UDVT

(7)

式中: 方阵D=(diag(σ1,σ2,…,σq),0)或其转置;0表示零矩阵,q=min(m,n),σ一般按降序排列。

为合理描述奇异值的变化,本文引入奇异值平均下降速度的概念,前k个奇异值的平均下降速度vk为:

(8)

实际下降速度的累计量为:

dk=vk×(k-1)=σ1-σk

(9)

总体平均下降速度累计量定义为:

(10)

实际信号分解得到的奇异值的下降速度整体上是逐渐减小的,并且有关系式:

c=dk-d0k>0, 1

(11)

在定义域(1,N)内存在最大值cmax,以cmax对应的横坐标kc为界,取前kc-1个奇异值重新构造矩阵,得到降噪后的信号。

用SVD对信号进行降噪,提高信号的信噪比,然后将降噪后的信号用VMD进行分解。针对VMD分解过程中的分解层数难以确定的问题,通过能量占比来评价分解效果,从而选取最优分解层数,并利用分量信号间的相关系数判断分量是否属于同一调频分量。其诊断流程图如图1所示。

图1 故障诊断框图

4 仿真实验

构造仿真信号f0(t)=f1(t)+f2(t),其中f1(t)和f2(t)的表达式为:

f1(t)=sin(4 π t)×sin(240 π t)

(12)

f2(t)=(1+cos(7 π t))cos(290 π t+2cos(10 π t))

(13)

该仿真信号由一个调频调幅部分和一个调幅部分组成,n(t)为加性高斯白噪声,服从正态分布:

n(t)～N(0,σ2)

(14)

其中σ为n(t)的标准差,这里令σ的值等于f0(t)的标准差,构造带噪声信号f(t):

f(t)=f0(t)+n(t)

(15)

对f(t)进行SVD降噪,用奇异值平均下降速度法确定有效奇异值个数,并重构信号。降噪前后的信号时域波形如图2所示。

图2 降噪前后仿真信号的时域波形

从图2(a)可以看到,信号的上下包络并不关于x轴对称,说明仿真信号f(t)中除了一个调频调幅成分和一个调幅成分外,还掺杂有随机噪声。通过奇异值平均下降速度法对f(t)降噪,降噪后的信号频谱如图2(b)所示,对比图2(a)和图2(b),经过降噪有效去除了信号中的随机噪声和脉冲干扰。

不同的分解层数下,分解得到的能量占比如表1所示。发现分解层数K=3时能量占比最高,分解层数K>3时能量占比开始下降,这是由于在分解层数设置过高时,导致若干分量的中心频率向同一个值收敛,使能量占比下降。

表1 不同K值下仿真信号分解结果的能量占比

根据表1,设置分解层数为3层对仿真信号进行VMD分解,然后计算中心频率相邻的分量信号的相关系数。如表2所示,分量的序号按其中心频率从小到大确定。

表2 相邻分量的相关系数

根据表2,ρ23>0.1,判断2号分量和3号分量分解自同一调频部分,即由f2(t)分解而来。1号分量分解自相对独立的f1(t)。为验证以上判断,求取u1与f1的瞬时频率以及u23与f2的瞬时频率分别进行比较,如图3所示。

图3 分量与原信号成分瞬时频率

从图3可以看出,分量信号与原信号中对应的频率成分的瞬时频率基本重合,图3(b)验证了根据相关系数得出的第2和第3分量分解自原信号中同一调频成分的判断。通过分量信号的瞬时频率,不仅可以判断出原信号由一个调幅调频部分和一个调幅部分组成,甚至可以进一步估算两者的解析式。而图中存在的几处分量与原信号成分的瞬时频率曲线不重合的地方,其在很短的时间内瞬时频率变化很大,可以判断为冲激噪声残留,并不影响最终的结果。

5 实测信号分析

本实验所用数据采自美国SQI公司生产的机械故障综合模拟试验台,实验平台如图4所示。

图4 故障诊断实验平台

5.1 内圈故障分析

滚动轴承转速为1 200 r/min, 人为在轴承内圈上设置故障点, 经计算得到其故障特征频率为98.3 Hz, 采样频率12.8 KHz,采样点数4 096,用1号传感器采集振动信号,其时域波形如图5所示。

图5 内圈故障信号时域波形

首先对原始信号使用改进的SVD降噪处理,根据奇异值平均下降速度法,选取2 048个奇异值的前347个进行重构。经过降噪处理后,信号充分滤除了杂波,然后对降噪后的信号进行VMD分解,首先确定分解层数。经过试验,当分解层数K=5时,分解结果的能量占比基本不再上升,然后开始缓慢下降。选择5层作为分解层数,对故障信号进行VMD分解,将各分量按中心频率从小到大排列,并求取中心频率相邻的分量信号间的相关系数,能量占比如表3所示,相关系数如表4所示。

表3 不同K值下仿真信号分解结果的能量占比

表4 相邻分量的相关系数

根据表4可知,分量1与分量2、分量3与分量4分解自同一调频部分,分别合并这些分量,然后用Hilbert变换求取各个有效分量的包络谱,由于分量信号能量相差很大,取包含能量最多的3个IMF分量信号的包络谱进行分析。如图6所示。

图6 改进方法下内圈故障数据主要分量包络谱

然后使用EMD分解算法对降噪后的信号进行分解,与本文方法进行对比。将分解得到的IMF分量按能量从大到小排列,选取包含能量最多的3个IMF分量求包络谱,如图7所示。

在图7(b)中,特征频率处不存在峰值,说明这是一个虚假分量;而图6中,根据本文方法分解得到的全部3个分量的包络谱中都出现了特征频率峰值。这是因为在分解过程中,部分频率成分被分解到了多个分量之中,本文方法通过合并这部分相邻分量,有效避免了虚假分量的产生。

5.2 外圈故障分析

外圈故障的时域波形如图8所示。

对该信号使用文中方法进行处理,首先使用SVD降噪,选取前311个奇异值重构信号,然后根据能量占比,确定分解层数K=6,包含能量最多的3个分量的包络谱如图9。滚动轴承内圈旋转频率为20 Hz,经计算得知滚动轴承内圈故障特征频率为60.96 Hz。在图9(c)分量的包络谱中,60.6 Hz处有明显的峰值,即为故障信号的特征频率。