CN105760843A - 一种滚动轴承早期故障特征提取方法 - Google Patents

Abstract

一种滚动轴承早期故障特征提取方法，本发明提供了一种基于稀疏优化的滚动轴承早期故障特征方法，包括如下步骤：步骤(1)采集滚动轴承加速度信号作为待分析信号；步骤(2)建立待分析信号的稀疏优化函数并求解得到待分析信号中的周期冲击成分；步骤(3)对周期冲击成分进行包络解调分析得到故障特征频率。本发明利用滚动轴承故障信号本身具有稀疏性的先验知识，而不需基于故障信号在字典变换下系数的稀疏性，避免了选择字典不合适而导致诊断失误的问题。

Description

一种滚动轴承早期故障特征提取方法

技术领域

本发明涉及滚动轴承故障诊断方法，尤其设计一种基于l₁范数最小化算法的滚动轴承早期故障特征提取方法。

背景技术

滚动轴承是机械设备中应用最广泛的零件之一，也是旋转机械易损件之一。滚动轴承在使用过程中总会经历正常、早期微弱故障、严重故障到失效的过程。严重故障阶段意味着滚动轴承的故障已经发展到中晚期，故障特征明显且容易提取；早期微弱故障阶段的特征提取相对来说比较困难，因为早期阶段故障特征微弱，且其他运动部件的信息以及环境干扰也会被引入到轴承系统形成背景噪声，从而使得轴承的早期故障难以监测和诊断。

针对滚动轴承发生故障时所具有的非平稳特征，国内外学者对滚动轴承早期故障特征提取方法进行了大量的研究。目前主要有小波变换、经验模态分解和Hilbert变换等非平稳信号分析方法在滚动轴承故障诊断中有广泛应用。

国内涉及信号稀疏表示的故障诊断的发明专利有“一种基于稀疏分解的风力发电机组故障特征提取方法”(201310471280.X)，是基于形态成分分析原理用不同的稀疏表达字典将信号分解为谐波、冲击以及噪声三个成分，其中采用离散余弦变换字典提取谐波成分，采用离散小波变换字典提取冲击成分，使得原本难以发现的故障特征凸显出来。发明专利“一种基于复合Q因子基算法的轴承故障诊断方法”(201210515071.6)是利用复合Q因子的小波变换产生相应的高Q因子基及低Q因子基，将原始信号在复合Q因子基上进行分解，利用相应的Q因子基提取出故障冲击信号成分。目前已有的基于稀疏表示的故障诊断方法，都是基于故障信号在给定的字典变换下系数是稀疏的先验条件，所以诊断结果的准确性很大程度上依赖于字典选择是否合适。若字典选择不合适将直接影响稀疏表示的结果，从而导致故障诊断结果的不准确。

发明内容

本发明基于上述技术问题，提出了一种滚动轴承早期故障特征提取方法，利用滚动轴承故障信号本身具有稀疏性的先验知识，而不需基于故障信号在字典变换下系数的稀疏性，构造了l₁范数最小化问题，通过求解该最小化问题获得分析信号中的周期冲击成分。

本发明的技术方案如下：

一种滚动轴承早期故障特征提取方法，包括如下步骤：

步骤(1)采集滚动轴承加速度信号作为待分析信号；

步骤(2)利用待分析信号中的周期冲击成分具有的稀疏性先验知识，建立其l₁范数最小化问题并求解得到待分析信号中的周期冲击成分；

步骤(3)对周期冲击成分进行包络解调分析得到故障特征频率。

进一步地，所述步骤(2)建立待分析信号的l₁范数最小化问题如下：

$x=\arg \underset{x}{min}\{F\left(x\right)=\frac{1}{2}|\left|H(y-x)\right|{|}_2^2+\mathrm{\λ}\left|\right|Wx\left|{|}_1\right\}$

其中，x为待分析信号中的周期冲击成；y为待分析信号；H为高通滤波器；λ为正则化参数；W为加权系数。

所述高通滤波器H为零相位非因果二阶巴特沃兹高通滤波器；具体设计如下：

差分方程描述为：

a₁y(n+1)+a₀y(n)+a₁y(n-1)＝-x(n+1)+2x(n)-x(n-1)

高通滤波器H如下形式：

H＝BA^-1

式中，矩阵B为：

$B=\left[\begin{array}{ccccccc}-1& 2& -1& & & & \\ & -1& 2& -1& & & \\ & & -1& 2& -1& & \\ & & & -1& 2& -1& \\ & & & & -1& 2& -1\end{array}\right]$

矩阵B的大小是(N-2)×N，N是待分析信号x的长度；矩阵A为：

$A=\left[\begin{array}{ccccc}{a}_0& a_1& & & \\ a_1& a_0& a_1& & \\ & a_1& a_0& a_1& \\ & & a_1& a_0& a_1\\ & & & a_1& a_0\end{array}\right]$

矩阵A的大小为(N-2)×(N-2)；差分方程中系数a₀和a₁需要满足a₀-2a₁＝4；设置在截止频率ω_c处滤波器的增益为0.5，由频响函数得到：

$a_0=\frac{4}{1+{\mathrm{cos\ω}}_c}$

将a₀和a₁带入到A中，得到高通滤波器H。

进一步地，求解所述待分析信号的l₁范数最小化问题采用最小化优化算法；具体如下：

步骤(21)设置迭代次数l＝0，设置加权系数W_i ^(l)＝1,i＝1,...,N；

步骤(22)求解如下优化问题：

$x^{\left(l\right)}=\arg \underset{x}{min}\{F\left(x\right)=\frac{1}{2}|\left|H(y-x)\right|{|}_2^2+\mathrm{\λ}\left|\right|W^{\left(l\right)}x\left|{|}_1\right\}$

将优化问题目标函数等效成光滑凸函数形式，将优化问题转化如下形式：

$G\left(x\right)=\frac{1}{2}\left|\right|{\mathrm{BA}}^{-1}(y-x)|{|}_2^2+\frac{\mathrm{\λ}}{2}{x}^T{W}^T\mathrm{\Λ}Wx+c$

矩阵Λ是对角矩阵形式：

${\mathrm{\Λ}}^{\left(l\right)}=\mathrm{\λ}\frac{{\mathrm{\φ}}^{\′}\left(x^{\left(l\right)}\right)}{x^{\left(l\right)}}=\mathrm{\λ}/\mathrm{\ψ}\left(x^{\left(l\right)}\right)$

式中，φ(x)＝||x||₁；

根据最小化优化更新等式为：

$x^{(l+1)}=\arg \underset{x}{min}G\left(x^{\left(l\right)}\right)$

求导获得G(x)得极小值，即可得到：

x^(l+1)＝A(B^TB+A^T(W^TΛ^(l)W)A)^-1B^TBA^-1y

步骤(23)对每个i＝1,...,N更新权重系数：

$W_i^{(l+1)}=\frac{1}{\left|x_i^{\left(l\right)}\right|+\mathrm{\ϵ}};$

步骤(24)直到收敛或迭代次数l达到最大迭代次数l_max，若不符合终止条件则跳到步骤(22)。

本发明的有益效果：

1)本发明利用滚动轴承故障信号本身具有稀疏性的先验知识，而不需基于故障信号在字典变换下系数的稀疏性，避免了选择字典不合适而导致诊断失误的问题。

2)在优化目标函数中包含了恢复冲击特征的保真项，保真项中的高通滤波器H是一个零相位非因果二阶巴特沃兹滤波器，滤波器中A和B是带状矩阵，这样做的优点在于有效地提高了算法的计算效率。

3)为了使得目标函数中正则化项l₁范数项可以无限趋近于l₀范数项，在l₁范数项中加入加权系数W，有利于提高信号的复原特性。

4)采用最小化优化算法对优化目标函数中l₁范数项进行等效转化，使得目标函数转化成光滑凸函数，则可通过直接求导的方式获得稀疏成分x。这种直接对目标函数求导获得故障信号的方式，可以有效的提高计算效率。

5)通过最小化优化方法和滤波算法，将滚动轴承中的混合成分进行分离，可以获得周期性冲击成分、谐波成分和噪声成分。l₁范数最小化算法不仅可以实现混合信号的分解，同时实现了降噪过程。

附图说明

图1为本发明的滚动轴承故障微弱故障特征提取的流程图。

图2所示为采集到的滚动轴承外圈早期故障的时域波形。

图3所示式对时域波形进行傅里叶变换获得的频谱图。

图4(a)(b)(c)分别表示周期冲击成分，谐波成分和噪声成分三种信号成分示意图。

图5为对周期冲击信号进行包络解调分析获得的包络解调谱。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

图1为本发明基于l₁范数最小化算法的滚动轴承故障微弱故障特征提取的流程图。下面结合流程图对l₁范数最小化算法原理进行详细说明。

利用加速度传感器对滚动轴承进行采集，获得振动加速度信号作为待分析信号y，图2所示为采集到的滚动轴承外圈早期故障的时域波形，图3所示是对时域波形进行傅里叶变换获得的频谱图，从图2和图3中均不能得到故障特征。采集到故障轴承的振动信号y中包含机械转频成分f、周期冲击成分x和噪声成分w。其中周期冲击信号主要是由轴承故障产生的，谐波信号主要是机械转频信号，而噪声主要是背景噪声。故障信号表现出周期冲击特征，故其具有稀疏性。同时，假设振动信号中的噪声方差为σ²，信号长度为N，则数据保真度约束为这样就形成如下的约束优化问题：

$\begin{array}{ccc}\arg \underset{x}{min}\left|\right|x|{|}_0& suchthat& \left|\right|y-x-f|{|}_2^2\≤{\mathrm{N\σ}}^2\end{array}---\left(1\right)$

通过推导不难得到H(y-x)≈y-x-f，其中H为高通滤波器。选择合适的λ，则可将式(1)的约束优化转化为如下形式的无约束优化问题：

$\hat{x}=\arg \underset{x}{min}\{F\left(x\right)=\frac{1}{2}|\left|H(y-x)\right|{|}_2^2+\mathrm{\λ}\left|\right|x\left|{|}_0\right\}---\left(2\right)$

式(2)的优化目标函数融合了恢复冲击特征的保真项和利用冲击特征稀疏性先验知识建立的正则化项。其中正则化项中l₀范数是非凸的，而且是NP难问题，用l₁范数替换式(2)中的l₀范数，则可形成一个可解的凸优化问题，即：

$\hat{x}=\arg \underset{x}{min}\{F\left(x\right)=\frac{1}{2}|\left|H(y-x)\right|{|}_2^2+\mathrm{\λ}\left|\right|x\left|{|}_1\right\}---\left(3\right)$

在式(2)的惩罚项是l₀范数项，其本质特征是对大的系数和小的系数的惩罚是一致的；为了解决l₁范数中对系数惩罚的不一致，在l₁范数项设置一个权重因子，使得对非零系数的惩罚更趋向于一致。则将式(3)转化为如下形式的优化函数：

$x=\arg \underset{x}{min}\{F\left(x\right)=\frac{1}{2}|\left|H(y-x)\right|{|}_2^2+\mathrm{\λ}\left|\right|Wx\left|{|}_1\right\}---\left(4\right)$

(1)为了求解式(1)中稀疏成分必须合理构造一个高通滤波器H，为了避免滤波后信号不必要的失真，构造一个零相位非因果巴特沃兹滤波器，差分方程描述为：

a₁y(n+1)+a₀y(n)+a₁y(n-1)＝-x(n+1)+2x(n)-x(n-1)(5)

高通滤波器H可以写成如下形式：

H＝BA^-1(6)

式中，A和B可以写成带状矩阵形式，带状矩阵是矩阵中大部分元素为零，而非零元素都集中在以主对角线为中心的带状区域中，在计算存储时只保留非零元素，故在运算过程中可以有效地提高算法的计算效率。矩阵B为：

$B=\left[\begin{array}{ccccccc}-1& 2& -1& & & & \\ & -1& 2& -1& & & \\ & & -1& 2& -1& & \\ & & & -1& 2& -1& \\ & & & & -1& 2& -1\end{array}\right]$

矩阵B的大小是(N-2)×N，N是输入信号x的长度。矩阵A为：

$A=\left[\begin{array}{ccccc}{a}_0& a_1& & & \\ a_1& a_0& a_1& & \\ & a_1& a_0& a_1& \\ & & a_1& a_0& a_1\\ & & & a_1& a_0\end{array}\right]$

矩阵A的大小为(N-2)×(N-2)。式(4)中系数a₀和a₁需要满足a₀-2a₁＝4。设置在截止频率ω_c处滤波器的增益为0.5，可由频响函数得到：

$a_0=\frac{4}{1+{\mathrm{cos\ω}}_c}$

设置高通滤波的截止频率f_c＝0.05Hz，则可以得到ω_c＝2πf_c＝0.1π，可以求得a₀＝2.050，a₁＝-0.975。将a₀和a₁代入到A中，即可通过式(6)求得高通滤波器H。

(3)由于优化目标函数中包含l₁范数项，是非光滑凸函数，无法通过对目标函数求导的方式获得稀疏成分x。故采用最小化优化算法对目标函进行等效处理，算法过程如下：

a.设置迭代次数l＝0，设置加权系数W_i ^(l)＝1,i＝1,...,N；

b.求解如下优化问题：

$x^{\left(l\right)}=\arg \underset{x}{min}\{F\left(x\right)=\frac{1}{2}|\left|H(y-x)\right|{|}_2^2+\mathrm{\λ}\left|\right|W^{\left(l\right)}x\left|{|}_1\right\}$

采用最小化优化算法将目标函数等效成一系列简单最小问题，简单的最小问题是光滑凸函数形式，将式(4)的优化问题转化如下形式：

$G\left(x\right)=\frac{1}{2}\left|\right|{\mathrm{BA}}^{-1}(y-x)|{|}_2^2+\frac{\mathrm{\λ}}{2}{x}^T{W}^T\mathrm{\Λ}Wx+c---\left(7\right)$

矩阵Λ是对角矩阵形式：

${\mathrm{\Λ}}^{\left(l\right)}=\mathrm{\λ}\frac{{\mathrm{\φ}}^{\′}\left(x^{\left(l\right)}\right)}{x^{\left(l\right)}}=\mathrm{\λ}/\mathrm{\ψ}\left(x^{\left(l\right)}\right)---\left(8\right)$

式中，φ(x)＝||x||₁。根据最小化优化更新等式为：

$x^{(l+1)}=\arg \underset{x}{min}G\left(x^{\left(l\right)}\right)---\left(9\right)$

对式(7)求导获得G(x)得极小值，即可得到：

x^(l+1)＝A(B^TB+A^T(W^TΛ^(l)W)A)^-1B^TBA^-1y(10)

c.对每个i＝1,...,N更新权重系数：

$W_i^{(l+1)}=\frac{1}{\left|x_i^{\left(l\right)}\right|+\mathrm{\ϵ}};$

d.直到收敛或迭代次数l达到最大迭代次数l_max，若不符合终止条件则跳到b。

根据式(10)可以获得周期冲击特征信号则谐波成分为：

$\hat{f}=LPF(y-\stackrel{\→}{x})---\left(11\right)$

式中，LPF为低通滤波器，y为加速度传感器采集到的振动信号，低通滤波器LPF为：

LPF＝I-H(12)

式中，I为单位矩阵，将式(12)带入到(11)中，即可获得谐波成分

噪声成分为：

$w\≈H(y-\stackrel{\→}{x})---\left(13\right)$

通过求解式(10)、(11)和(13)即可获得周期冲击成分，谐波成分和噪声成分，实现滚动轴承振动信号中多成分分解。三种信号成分分别对应于图4(a)(b)(c)所示。

对周期冲击信号进行包络解调分析，获得图5所示的包络解调谱，从图中即可看到故障特征频率及其倍频。

Claims (4)

1.一种滚动轴承早期故障特征提取方法，其特征在于：包括如下步骤：

步骤(1)采集滚动轴承加速度信号作为待分析信号；

步骤(2)利用待分析信号中的周期冲击成分具有的稀疏性先验知识，建立其l₁范数最小化问题并求解得到待分析信号中的周期冲击成分；

步骤(3)对周期冲击成分进行包络解调分析得到故障特征频率。

2.根据权利要求1所述的滚动轴承早期故障特征提取方法，其特征在于：所述步骤(2)建立待分析信号中的周期冲击成分的l₁范数最小化问题如下：

$x=\arg \underset{x}{min}\{F\left(x\right)=\frac{1}{2}|\left|H(y-x)\right|{|}_2^2+\mathrm{\λ}\left|\right|Wx\left|{|}_1\right\}$

其中，x为待分析信号中的周期冲击成；y为待分析信号；H为高通滤波器；λ为正则化参数；W为加权系数。

3.根据权利要求2所述的滚动轴承早期故障特征提取方法，其特征在于：所述高通滤波器H为零相位非因果二阶巴特沃兹高通滤波器；具体设计如下：

差分方程描述为：

a₁y(n+1)+a₀y(n)+a₁y(n-1)＝-x(n+1)+2x(n)-x(n-1)

高通滤波器H如下形式：

H＝BA^-1

式中，矩阵B为：

$B=\left[\begin{array}{ccccccc}-1& 2& -1& & & & \\ & -1& 2& -1& & & \\ & & -1& 2& -1& & \\ & & & -1& 2& -1& \\ & & & & -1& 2& -1\end{array}\right]$

矩阵B的大小是(N-2)×N，N是待分析信号x的长度；矩阵A为：

$A=\left[\begin{array}{ccccc}{a}_0& a_1& & & \\ a_1& a_0& a_1& & \\ & a_1& a_0& a_1& \\ & & a_1& a_0& a_1\\ & & & a_1& a_0\end{array}\right]$

矩阵A的大小为(N-2)×(N-2)；差分方程中系数a₀和a₁需要满足a₀-2a₁＝4；设置在截止频率ω_c处滤波器的增益为0.5，由频响函数得到：

$a_0=\frac{4}{1+{\mathrm{cos\ω}}_c}$

将a₀和a₁带入到A中，得到高通滤波器H。

4.根据权利要求2所述的滚动轴承早期故障特征提取方法，其特征在于：求解所述待分析信号的l₁范数最小化问题采用最小化优化算法；具体如下：

步骤(21)设置迭代次数l＝0，设置加权系数W_i ^(l)＝1,i＝1,...,N；

步骤(22)求解如下优化问题：

$x^{\left(l\right)}=\arg \underset{x}{min}\{F\left(x\right)=\frac{1}{2}|\left|H(y-x)\right|{|}_2^2+\mathrm{\λ}\left|\right|W^{\left(l\right)}x\left|{|}_1\right\}$

将优化问题目标函数等效成光滑凸函数形式，将优化问题转化如下形式：

$G\left(x\right)=\frac{1}{2}\left|\right|{\mathrm{BA}}^{-1}(y-x)|{|}_2^2+\frac{\mathrm{\λ}}{2}{x}^T{W}^T\mathrm{\Λ}Wx+c$

矩阵Λ是对角矩阵形式：

${\mathrm{\Λ}}^{\left(l\right)}=\mathrm{\λ}\frac{{\mathrm{\φ}}^{\′}\left(x^{\left(l\right)}\right)}{x^{\left(l\right)}}=\mathrm{\λ}/\mathrm{\ψ}\left(x^{\left(l\right)}\right)$

式中，φ(x)＝||x||₁；

根据最小化优化更新等式为：

$x^{(l+1)}=\arg \underset{x}{min}G\left(x^{\left(l\right)}\right)$

求导获得G(x)得极小值，即可得到：

x^(l+1)＝A(B^TB+A^T(W^TΛ^(l)W)A)^-1B^TBA^-1y

步骤(23)对每个i＝1,...,N更新权重系数：

$W_i^{(l+1)}=\frac{1}{\left|x_i^{\left(l\right)}\right|+\mathrm{\ϵ}};$

步骤(24)直到收敛或迭代次数l达到最大迭代次数l_max，若不符合终止条件则跳到步骤(22)。