利用粒子滤波与谱峭度的滚动轴承故障诊断方法

摘要

本发明公开了一种利用粒子滤波与谱峭度的滚动轴承故障诊断方法，涉及粒子滤波降噪处理以及谱峭度计算。该方法提出了在快速谱峭度的基础上，利用粒子滤波降噪处理提高信噪比，解决了快速谱峭度在低信噪比情况下可行性差的问题。首先建立振动信号的状态方程；然后提取背景噪声，将其与状态方程之和作为观测方程；最后联立状态方程与观测方程建立状态空间模型；采用粒子滤波算法对信号重新估计，得到新的序列即是降噪之后的信号；最后用快速谱峭度的方法，获取最佳分析频带，得出故障频率。本发明的滚动轴承故障诊断方法，降低故障信号中的噪声干扰，提高信噪比，实现了滚动轴承早期微弱故障的诊断。

说明书

技术领域

本发明涉及一种利用粒子滤波与谱峭度的滚动轴承故障诊断方法，涉及粒子滤波降噪处 理以及谱峭度计算。

背景技术

滚动轴承是机械设备尤其是旋转机械中重要的零部件之一，但是其寿命的随机性较大， 容易损坏，目前还无法准确预测其寿命的长短。近年来，科技和工业生产都在快速的发展， 机械设备逐渐向高速化、大型化和自动化发展，这在提高生产力的同时对设备的安全维护也 提出了更高的要求，在某一环节出现一个微小的问题都将造成不可挽回的巨额损失。鉴于滚 动轴承在机械设备中的地位，可以得出滚动轴承能否正常运行工作关系着机械设备乃至生产 线的正常运转，因而掌握滚动轴承运行的工作状态以及故障的形成和发展，是目前机械故障 诊断领域中所研究的重要课题之一。

根据故障诊断技术机理的不同，诊断滚动轴承故障的常用方法有振动分析法、噪声分析 法、油样分析法、温度分析法、油膜电阻法、声发射诊断等。这些分析方法中，振动分析法 是运用最广泛的一种方法。振动信号包含丰富的设备异常或故障的信息，通过对振动信号的 实时在线监测和分析，可以判断机械系统及其部件是否正常运行，并以此确定机械故障产生 的原因，判断故障产生的具体部位和损伤程度。较常用的振动诊断方法有倒频谱分析、特征 参数分析法、冲击脉冲法、包络分析法、小波分析等。虽然诊断方法很多，但是对早期的微 弱故障诊断仍然是一个难题。因为在故障早期，所需频段信号能量非常微弱、特征信号成分 极易被噪声和其他频段信号俺没，未经任何处理的原始频段微弱信号中很难提取出有用的特 征信号成分，如何提取被淹没在噪声无用频段内的微弱故障特征信号是目前国内外研究焦点。

发明内容

振动信号是典型的非平稳、非线性的时序信号，粒子滤波在解决非线性滤波问题时，随 机量不会受到高斯分布的制约，比高斯模型具有更广泛的分布，对非线性特性有很强的建模 能力。本发明提供了一种利用粒子滤波与谱峭度的滚动轴承故障诊断方法，解决了普通的快 速谱峭度方法难以实现在低信噪比的情况下诊断出轴承故障的问题。

一种利用粒子滤波与谱峭度的滚动轴承故障诊断方法，具体过程如下：

粒子滤波是以状态空间方程为主对信号滤波，所以首先要建立轴承振动信号的状态空间 方程。对于一个非线性、非高斯系统，状态空间模型的一般表达式为：

$\left(\begin{array}{c}{x}_k=f(x_{k-1},v_k)\\ z_k=h(x_k,n_k)\end{array}\right)---\left(1\right)$

其中，v_k和n_k分别是过程噪声和观测噪声，并且都是相互独立的，协方差分别是Q_k和R_k的零均值噪声序列，f和h为状态转移方程和观测方程。

众所周知，轴承振动信号是典型的非线性、非高斯时间序列。所以用AR时序模型作为 状态方程，模型的残差方差即为过程噪声。提取原始信号的背景噪声作为观测噪声，将其与 状态方程加载一起作为观测方程。综上所述，振动信号的状态空间表达式可以写为：

$\left(\begin{array}{c}{x}_t=\underset{t=1}{\overset{p}{\Sigma }}{a}_t\left(t\right)x_{t-1}+{\sigma }_{e_t}{e}_t\\ y_t=x_t+{\sigma }_{n_t}{n}_t\end{array}\right)---\left(2\right)$

式中，p为AR模型阶数，a_t(t)是模型参数序列，为模型残差方差，即过程噪声，为信号背景噪声，具体获得方法为：用小波分解函数对采集的振动信号进行多尺度分解,得 到小波变换系数。再设定合适的阈值,将大于阈值的小波系数置为零,小于阈值的小波系数保 留下来。对经过处理的小波系数进行小波重构,就得到所需的噪声信号，即观测噪声。e_t和n_t都是服从N～(0,1)的标准正太分布。模型参数估计方法为Yule-Walker法，模型定界准则是 AIC准则。

粒子滤波的本质是利用蒙特卡罗积分法求解贝叶斯积分运算，它的关键是核心粒子的权 值大小要能不断被修正和调整，最后获得精确的估计值，下文将标准粒子滤波的具体流程。

首先用一组随机粒子和这些粒子的权值大小来近似表达后验概率密度函数 p(x_k|z_1:k)：$p\left(x_k\right|z_{1:k})\approx \underset{i=1}{\overset{N_s}{\Sigma }}{\omega }_k^i\times \delta (x_k-x_k^i)---\left(3\right)$

式中，δ(·)是狄拉克函数，z_1:k是观测量。如果这些粒子是通过真实概率密度函数 p(x_k|z_1:k)获得，那么式(3)可写为：

$p\left(x_k\right|z_{1:k})\approx \underset{i=1}{\overset{N_s}{\Sigma }}\frac{\delta (x_k-x_k^i)}{N}---\left(4\right)$

然而，在很多情况下获得真实概率密度函数粒子是非常困难的，所以就利用重要性概率 函数q(x_k|z_1:k)获取粒子。假设p(x_k|z_1:k)与一个分析函数π(x_k|z_1:k)是线性关系，则是(3)中 的权值大小可以通过下式计算：

${\omega }_k^i\propto \frac{\pi \left(x_k^i\right|z_{1:k})}{q\left(x_k^i\right|z_{1:k})}\propto \frac{p\left(x_k^i\right|z_{1:k})}{q\left(x_k^i\right|z_{1:k})},i=\mathrm{1,2,3},...N_s---\left(5\right)$

假设重要性函数可以分解为：

$q(x_k^i,x_{k-1}^i|z_{1:k})=q\left(x_k^i\right|x_{k-1}^i,z_{1:k})\times p(x_{k-1}^i,z_{1:k-1}),i=\mathrm{1,2,3},...N_s---\left(6\right)$

将式(5)和时(6)代入到式(4)中更新粒子权值大小

$\left(\begin{array}{c}{\omega }_k^i\propto \frac{p\left(x_k^i\right|z_{1:k})}{q\left(x_k^i\right|z_{1:k})}\propto \frac{p\left(z_k\right|x_k^i)\times p\left(x_k^i\right|z_{1:k-1})\times p\left(x_{k-1}^i\right|z_{1:k-1})}{q\left(x_k^i\right|x_{k-1}^i,z_{1:k})\times q\left(x_{k-1}^i\right|z_{1:k-1})}\\ \propto {\omega }_{k-1}^i\times \frac{p\left(z_k\right|x_k^i)\times p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,z_{1:k})},i=\mathrm{1,2,3},...N_s\end{array}\right)---\left(7\right)$

当满足条件:的时候，式(7)可以写为：

${\omega }_k^i={\omega }_{k-1}^i\times \frac{p\left(z_k\right|x_k^i)\times p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,z_k)}/(\underset{i=1}{\overset{N_s}{\Sigma }}{\omega }_{k-1}^i\times \frac{p\left(z_k\right|x_k^i)\times p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,z_k)}),i=\mathrm{1,2,3},...N_s---\left(8\right)$

又因为先验概率密度函数通常选择为：所以式(8)可以写为：

${\omega }_k^i={\omega }_{k-1}^i\times p\left(z_k\right|x_k^i)/(\underset{i=1}{\overset{N_s}{\Sigma }}{\omega }_{k-1}^i\times p\left(z_k\right|x_k^i)),i=\mathrm{1,2,3},...N_s---\left(9\right)$

粒子滤波结束后结合快速谱峭度法提取特征故障频率。谱峭度的定义式为：

$K_x\left(f\right)=\frac{S_4\left(f\right)}{{\left(S_2\left(f\right)\right)}^2}-2,f\ne 0---\left(10\right)$

式中，为信号的n阶谱矩，其中，E<·>——均值运算，E|·|取模 运算。X(f,t)是信号x(t)在频率f处的复包络。

快速谱峭度算法采用的是在全频带范围内1/3和二叉树结构设计的带通滤波器分析信号。 算法如下：

(1)构造一个截止频率f_c＝1/8+ε的低通滤波器h(n)，其中f_s＝1,ε＞0。低通滤波器 的截止频率增大ε是为了消除过度带的影响。并基于h(n)构造准解析低通滤波器h₀(n)和准解 析高通滤波器h₁(n)，其分析频带分别为[0,1/4]、[1/4,1/2]，即：

$\left(\begin{array}{c}{h}_0\left(n\right)=h\left(n\right)e^{\mathrm{j\pi n}/4}\\ h_1\left(n\right)=h\left(n\right)e^{j3\mathrm{\pi n}/4}\end{array}\right)---\left(11\right)$

(2)分别以h₀(n)和h₁(n)低通和高通滤波，对滤波结果做2倍降采样，如图1所示。如此 迭代进行，可以获得如图2所示的滤波器树及相应的滤波结果。其中表示第k层第_i个 滤波器产生的滤波结果，其中i取值为0到2^k-1；同时也是信号在中心频带 f_i＝(i+2^-1)2^k-1和带宽(Δf)_k＝2^k-1的复包络，2倍降采样的目的是保证滤波器每一层中的数据 长度与原始数据相同。

(3)对于每一个滤波器结果，按照式(10)计算峭度：

$K=\frac{E<{\left|c_k^i\left(n\right)\right|}^4>}{{\{E<{\left|c_k^i\left(n\right)\right|}^2>\}}^2}-2---\left(11\right)$

(4)将所有的谱峭度汇总，从而得到“快速谱峭度图”。

根据峭度图选择峭度最大值频带区间解调分析，计算功率谱，对功率谱做FFT变换得到 频谱图，根据频谱信息做出故障类型的诊断。

本发明提供了一种利用粒子滤波与谱峭度的滚动轴承故障诊断方法，降低故障信号中的 噪声干扰，提高信噪比，实现了滚动轴承早期微弱故障的诊断。

附图说明

图1为低通/高通滤波器分解图。

图2为二叉树滤波器组及相应的滤波结果。

图3内环故障时域振动信号。

图4外环故障时域振动信号。

图5内环振动信号的AIC曲线图。

图6外环振动信号的AIC曲线图。

图7滤波之后的内环时域信号。

图8滤波之后的外环时域信号。

图9为内环信号降噪后频谱图。

图10为外环信号降噪后频谱图。

图11内环信号滤波之前的频谱图。

图12外环信号滤波之前的频谱图。

具体实施方式

以下结合技术方案和附图详细叙述本发明的实施例。

实施例的数据来自美国西储大学，数据为电机驱动端轴承故障数据,轴承型号为 6205-2RS JEM SKF，深沟球轴承。轴承参数如下表所示：

表1 轴承参数表

内圈直径 外圈直径 厚度 滚动体直径D_d节经d_m接触角α 滚珠数Z 25mm 52mm 15mm 7.94mm 139.04mm 60° 9

实验电机转速1797rpm/m，故障尺寸为0.1778mm，内外环的故障特征频率计算公式如下： 内环$\mathrm{BPFI}=\frac{{\mathrm{Zf}}_r}{2}(1+\frac{D_d}{d_m}\cos \alpha )=159.7\mathrm{Hz}$

外环$\mathrm{BPFO}=\frac{{\mathrm{Zf}}_r}{2}(1-\frac{D_d}{d_m}\cos \alpha )=107.4\mathrm{Hz}$

观察图3和图4外所示的内环、外环时域振动信号，无法看出冲击信号，所以不能判断 是否存在故障。要得到故障特征频率只有作进一步的分析处理。

图5和图6分别是内环和外环的AIC曲线图。AIC准则认为：AIC值越小，模型精度越 高。但是这样模型的阶数就会变大，计算量就会上升。在上升到一定阶段时，模型的精确度 度提高非常缓慢。综合考虑两者因素，我们选定当AIC值下降缓慢时候的阶数为模型阶数。 内环的AIC曲线在阶数为20时，值已经趋于平缓，下降幅度很小；而外环在第8阶时下降 幅度最小。所以模型阶数分别定位20阶、8阶。

其次采用小波降噪方法得到信号的背景噪声，即为观察噪声，观察噪声与AR模型之和 即为观察方程；最后联立AR模型与观察方程便是粒子滤波所需的状态空间方程。利用粒子 滤波方法对信号重新估计滤波，得到滤波后的时域信号如图7、图8所示，对比观察滤波前 后的信号，可以发现滤波之后信号幅值下降的明显。结合谱峭度方法计算峭度值，根据峭度 图选取峭度值最大子带为最优解调分析频带。提取特征频率如图9、图10所示。

在图9、图10频谱图中，很明显的观察到故障特征频率以及它们的倍频信号，可以说明 本方法能成功地实现轴承故障诊断。为了进一步体现所提出方法的优越性，对比分析了粒子 滤波之前的处理结果，如图11、图12所示。

图11是内环滤波之前的频域图，图中只能观察到轴频信号(因为实验转速为1797rpm/m， 所以轴频是29.95Hz)及其倍频，没有故障频率。说明滤波之前谱峭度无法得到内环的故障 信息。图12是外环振动信号滤波之前的频域图，频谱图中虽然可以看到故障频率107.7Hz， 但是存在着边频带谱线，干扰信号太多，可视效果很差。再反观滤波之后频谱如图10所示， 边频带频率被成功滤除，没有干扰信号，结果清晰一目了然。

上述对比分析再一次证明了所提出的方法的优越性。