一种基于扩展卡尔曼滤波算法的轴承故障诊断及预测方法

摘要

一种基于扩展卡尔曼滤波算法的轴承故障诊断及预测方法，包括以下步骤：1)采集轴承全寿命周期振动信号；2)利用振动信号构建AR模型，对振动信号进行滤波分析，突显与故障相关的信号；3)运用小波包变换，提取与小波包系数相关的能量信息，构建特征参数；4)进行马氏距离的计算，基于此构建健康指数，对值为非负、非高斯分布的健康指数，运用Box-Cox变换，转换成高斯分布的数据，利用高斯分布的性质和逆Box-Cox变换，确定相关异常阈值范围；5)对耗损期的健康指数数据拟合分析，构建退化模型和状态空间模型，利用当前数据和扩展卡尔曼滤波算法更新模型参数并预测轴承剩余寿命。本发明预测精度较高、耗时较短。

说明书

技术领域

本发明属于轴承故障诊断及预测领域，尤其涉及一种基于扩展卡尔曼滤波算法的轴承故障诊断及预测方法。

背景技术

轴承是旋转机械中不可缺少的零部件，在电力、石化、冶金、机械、航空航天以及一些军事工业部门中广泛使用，是保证精密机床、高速铁路、风力发电机等重要装备设施精度、性能、寿命和可靠性的核心零部件，但同时也是最容易损伤的部件之一。据统计，旋转机械的许多故障都是由于轴承损坏引起的。轴承发生故障，轻则降低或失去系统的某些功能，重则造成严重的甚至是灾难性的后果。因此轴承状态监测、故障诊断以及故障预测一直是近年来的研究重点。振动信号广泛地应用于轴承的状态监测以及故障预测中，一般采用选择时域统计量作为寿命特征参数的方法进行处理，但是由于振动信号具有非线性、模糊性等特点，且轴承在耗损期故障形式差别较大，存在不同轴承样本间的实际工作寿命差异明显等问题，应用已有的方法通常无法解决这样的问题。此外，其他方法，如神经网络、人工智能等，在预测精度、耗时等方面性能表现不尽如人意。

发明内容

为了克服现有的轴承故障诊断及预测方法的预测精度较低、耗时较长的不足，本发明提供了一种预测精度较高、耗时较短的基于扩展卡尔曼滤波算法的轴承故障诊断及预测方法。

为了解决上述技术问题提供的技术方案为：

一种基于扩展卡尔曼滤波算法的轴承故障诊断及预测方法，所述方法包括以下步骤：

S1.采集轴承的全寿命周期振动信号；

S2.利用轴承健康时的振动信号构建AR模型，运用该AR模型对采集到的振动信号进行滤波，突出振动信号中与故障相关的信息，考虑一个时间序列x[1],x[2],…,x[m]，其中，m为时间序列的总数，p阶AR模型公式如下：

$>x\left[t\right]=\underset{i=1}{\overset{p}{\Sigma }}a\left[i\right]x[t-i]+e\left[t\right]---\left(1\right)$>

上式中，x[t]为第t个时间序列，1≤t≤m，a[i]是第i个模型系数，1≤i≤p，通过使用Levinson-Durbin递归方法估计得到，e[t]是均值为零，方差为σ的白噪声，p是AR模型的最优阶数，通过使用赤池信息量准则AIC，选取AIC为最小的AR模型阶数，其定义为：

AIC(j)＝M·ln(σ²)+2·j(2)

上式中，AIC(j)为j阶AR模型所对应的AIC值，M为所取时间序列样本的数量，利用该AR模型滤波分析轴承的振动信号，当轴承健康时，得到的残差信号是白噪声，当轴承发生故障时，得到的残差信号是与当前故障密切相关的故障信号；

S3.运用小波包变换对AR模型滤波后的残差信号进行分析，构建小波包系数对应的能量特征用于马氏距离计算；

S4.进行马氏距离的计算，基于马氏距离构建出表示轴承健康状况的指数，所构建的健康指数在整个轴承寿命周期上，呈现为两头高，中间低的曲线，对轴承健康时非高斯分布的健康指数，运用Box-Cox变换，转换成高斯分布的数据，利用高斯分布的性质和逆Box-Cox变换，确定轴承发生异常时的健康指数的阈值；

S5.拟合分析轴承耗损期的健康指数数据，构建退化模型并建立状态空间模型，利用当前观测到的健康指数数据和扩展卡尔曼滤波算法更新模型参数，并预测剩余寿命，过程如下：

对耗损期的健康指数数据，拟合分析构建如下的退化模型：

HI(k)＝a·exp(b·k)+c·exp(d·k)(3)

上式中，HI(k)为轴承在k时刻的健康指数，k为时间参数，a,b,c,d为耗损期轴承健康指数数据拟合出的退化模型参数，基于该退化模型构建状态方程：

$>a_k=a_{k-1}+w_{k-1}^a---\left(4\right)$>

$>b_k=b_{k-1}+w_{k-1}^b---\left(5\right)$>

$>c_k=c_{k-1}+w_{k-1}^c---\left(6\right)$>

$>d_k=d_{k-1}+w_{k-1}^d---\left(7\right)$>

上式中，a_k,b_k,c_k,d_k和a_k-1,b_k-1,c_k-1,d_k-1为分别在k时刻和k-1时刻的状态变量a,b,c,d的值，为在k-1时刻，独立的且分别对应状态变量a,b,c,d的噪声；

同时构建测量方程

HI_k＝a_k·exp(b_k·k)+c_k·exp(d_k·k)+v_k(8)

上式中，HI_k为在k时刻健康指数的测量值，v_k为在k时刻的测量噪声；

利用扩展卡尔曼滤波算法更新状态方程和测量方程参数至k时刻，按公式(3)计算k+l时刻的健康指数HI(k+l)

HI(k+l)＝a_k·exp(b_k·(k+l))+c_k·exp(d_k·(k+l))(9)

上式中，l＝1,2,…,∞；计算使得不等式(10)成立的l的值，并记录l的最小值为在k时刻预测的轴承剩余寿命；

HI(k+l)＞故障阀值(10)。

进一步，所述S3中，对AR模型滤波后的残差信号进行小波包变换分析，将每个小波包系数所对应的能量信息提取后，作为特征参数反映轴承当前的健康状态：

E_s＝∫|x(t)|²dt(11)

上式中，E_s是小波包系数重构信号x(t)的能量。

再进一步，所述S4中，对S3所得的特征参数，基于这些特征参数计算马氏距离；

设第i组数据标准化的特征向量为x_i，其包含着n个特征参数，则该组数据对应的马氏距离MD(i)为：

$>MD\left(i\right)=\frac{1}{n}{x}_i{C}^{-1}{x}_i^T---\left(12\right)$>

上式中，x_i＝[x_i1,x_i2,…,x_in]，x_i1表示特征向量x_i中第1个元素，x_i2表示特征向量x_i中第2个元素，x_in表示特征向量x_i中第n个元素，x_i^T是x_i的转置矩阵，C为轴承健康时特征参数间的协方差矩阵，C^-1是C的逆矩阵；

基于马氏距离计算构建出表示轴承健康状况的指数，对轴承健康时的非高斯分布的健康指数，运用Box-Cox变换，转换成高斯分布的数据，利用高斯分布的性质和逆Box-Cox变换，确定轴承发生异常时健康指数的阈值；

健康指数包括水平和垂直两个方向的信息，其由下式计算得到：

$>HI\left(k\right)={({\mathrm{MD}}_H{\left(k\right)}^{1/2n_1}+{\mathrm{MD}}_V{\left(k\right)}^{1/2n_2})}^{1/2}---\left(13\right)$>

上式中，n₁和n₂分别为在水平振动信号和垂直振动信号上提取的特征个数，MD_H和MD_V分别为在水平振动信号和垂直振动信号上计算得到的马氏距离；

设健康指数的变量x＝[x₁,x₂,…,x_N]，N为健康指数样本的个数，变量x中的第i个数据x_i经Box-Cox变换后得到的数据通过下式计算出来：

$>y_i=\left(\begin{array}{cc}\frac{x_i^{\lambda }-1}{\lambda },& \lambda \ne 0\\ \ln \left(x_i\right),& \lambda =0\end{array}\right)---\left(14\right)$>

上式中，y_i为对应于x_i经Box-Cox变换后得到的值，1≤i≤N，λ表示为变换参数，通过下式的最大似然函数LLF(y,λ)估计得到：

$>LLF(y,\lambda )=-\frac{N}{2}·\ln \left(\underset{i=1}{\overset{N}{\Sigma }}\frac{{(y_i-\bar{y})}^2}{N-1}\right)+(\lambda -1)·\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{lnx}}_i---\left(15\right)$>

上式中，y＝[y₁,y₂,…,y_N]，

本发明的技术构思为：通过采集轴承振动信号，利用AR模型和小波包变换，突出轴承故障信息并构建特征参数，基于马氏距离提出健康指数，利用Box-Cox变换确定轴承发生异常时的健康指数的阈值，利用该阀值把轴承整个寿命周期区分为如下三个阶段：磨合期，有效寿命期和耗损期。拟合分析轴承耗损期的健康指数数据，构建退化模型并建立状态空间模型，利用当前观测到的健康指数数据和扩展卡尔曼滤波算法更新模型参数，并预测剩余寿命。

本发明的有益效果为：预测精度较高、耗时较短。

附图说明

图1为基于扩展卡尔曼滤波算法的轴承故障诊断及预测方法流程图；

图2为轴承全寿命周期健康指数示意图；

图3为轴承信号的频谱演变图；

图4为轴承信号经AR模型滤波后残差信号的频谱演变图；

图5为轴承健康时健康指数数据的正态概率图；

图6为针对轴承耗损期健康指数数据分别利用指数函数和两个指数函数之和拟合出的曲线；

图7为轴承健康时健康指数数据经过Box-Cox变换后的直方图；

图8为轴承寿命实验得到的7个轴承在耗损期的健康指数数据；

图9为轴承在耗损期阶段的剩余寿命预测。

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1～图9，一种基于扩展卡尔曼滤波算法的轴承故障诊断及预测方法，所述方法包括以下步骤：

S1.采集轴承的全寿命周期振动信号；

S2.利用轴承健康时的振动信号构建AR模型，运用该AR模型对采集到的振动信号进行滤波，突出振动信号与故障相关的信息，从而有利于后续健康指数构建，故障诊断及预测工作；

S3.运用小波包变换对AR模型滤波后的残差信号进行分析，构建小波包系数对应的能量特征用于马氏距离计算；

S4.进行马氏距离的计算，基于马氏距离构建出表示轴承健康状况的指数，对轴承健康时非高斯分布的健康指数，运用Box-Cox变换，转换成高斯分布的数据，利用高斯分布的性质和逆Box-Cox变换，确定轴承发生异常时，健康指数的阈值，利用该阀值把轴承整个寿命周期分为三个阶段：磨合期，有效寿命期和耗损期；

S5.拟合分析轴承耗损期的健康指数数据，构建退化模型并建立状态空间模型，利用当前观测到的健康指数数据和扩展卡尔曼滤波算法更新模型参数，并预测剩余寿命。

所述S1中，如附图2所示，轴承的全寿命周期可以分为三个阶段：磨合期，有效寿命期和耗损期。

所述S2中，利用轴承健康时的振动信号构建AR模型，运用该AR模型对采集到的振动信号进行滤波，突出振动信号中与故障相关的信息，考虑一个时间序列x[1],x[2],…,x[m]，其中，m为时间序列的总数，p阶AR模型公式如下：

$>x\left[t\right]=\underset{i=1}{\overset{p}{\Sigma }}a\left[i\right]x[t-i]+e\left[t\right]---\left(1\right)$>

上式中，x[t]为第t个时间序列，1≤t≤m，a[i]是第i个模型系数，1≤i≤p，通过使用Levinson-Durbin递归方法估计得到，e[t]是均值为零，方差为σ的白噪声，p是AR模型的最优阶数，通过使用赤池信息量准则AIC，选取AIC为最小的AR模型阶数，其定义为：

AIC(j)＝M·ln(σ²)+2·j(2)

上式中，AIC(j)为j阶AR模型所对应的AIC值，M为所取时间序列样本的数量，利用该AR模型滤波分析轴承的振动信号，当轴承健康时，得到的残差信号是白噪声，当轴承发生故障时，得到的残差信号是与当前故障密切相关的故障信号；

所述S3中，对AR模型滤波后的残差信号进行小波包变换分析，将每个小波包系数所对应的能量信息提取后，作为特征参数反映轴承当前的健康状态：

E_s＝∫|x(t)|²dt(11)

上式中，E_s是小波包系数重构信号x(t)的能量。

所述S4中，对S3所得的特征参数，基于这些特征参数计算马氏距离；

设第i组数据标准化的特征向量为x_i，其包含着n个特征参数，则该组数据对应的马氏距离MD(i)为：

$>MD\left(i\right)=\frac{1}{n}{x}_i{C}^{-1}{x}_i^T---\left(12\right)$>

上式中，x_i＝[x_i1,x_i2,…,x_in]，x_i1表示特征向量x_i中第1个元素，x_i2表示特征向量x_i中第2个元素，x_in表示特征向量x_i中第n个元素，x_i^T是x_i的转置矩阵，C为轴承健康时特征参数间的协方差矩阵，C^-1是C的逆矩阵；

基于马氏距离计算构建出表示轴承健康状况的指数，对轴承健康时的非高斯分布的健康指数，运用Box-Cox变换，转换成高斯分布的数据，利用高斯分布的性质和逆Box-Cox变换，确定轴承发生异常时健康指数的阈值；

健康指数包括水平和垂直两个方向的信息，其由下式计算得到：

$>HI\left(k\right)={({\mathrm{MD}}_H{\left(k\right)}^{1/2n_1}+{\mathrm{MD}}_V{\left(k\right)}^{1/2n_2})}^{1/2}---\left(13\right)$>

上式中，n₁和n₂分别为在水平振动信号和垂直振动信号上提取的特征个数，MD_H和MD_V分别为在水平振动信号和垂直振动信号上计算得到的马氏距离；

设健康指数的变量x＝[x₁,x₂,…,x_N]，N为健康指数样本的个数，变量x中的第i个数据x_i经Box-Cox变换后得到的数据通过下式计算出来：

$>y_i=\left(\begin{array}{cc}\frac{x_i^{\lambda }-1}{\lambda },& \lambda \ne 0\\ \ln \left(x_i\right),& \lambda =0\end{array}\right)---\left(14\right)$>

上式中，y_i为对应于x_i经Box-Cox变换后得到的值，1≤i≤N，λ表示为变换参数，通过下式的最大似然函数LLF(y,λ)估计得到：

$>LLF(y,\lambda )=-\frac{N}{2}·\ln \left(\underset{i=1}{\overset{N}{\Sigma }}\frac{{(y_i-\bar{y})}^2}{N-1}\right)+(\lambda -1)·\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{lnx}}_i---\left(15\right)$>

上式中，y＝[y₁,y₂,…,y_N]，

所述S5中，拟合分析轴承耗损期的健康指数数据，构建退化模型并建立状态空间模型，利用当前观测到的健康指数和扩展卡尔曼滤波算法更新模型参数，并预测剩余寿命，过程如下：

对耗损期的健康指数数据，拟合分析构建如下的退化模型：

HI(k)＝a·exp(b·k)+c·exp(d·k)(3)

上式中，HI(k)为轴承在k时刻的健康指数，k为时间参数，a,b,c,d为耗损期轴承健康指数数据拟合出的退化模型参数，基于该退化模型构建状态方程：

$>a_k=a_{k-1}+w_{k-1}^a---\left(4\right)$>

$>b_k=b_{k-1}+w_{k-1}^b---\left(5\right)$>

$>c_k=c_{k-1}+w_{k-1}^c---\left(6\right)$>

$>d_k=d_{k-1}+w_{k-1}^d---\left(7\right)$>

上式中，a_k,b_k,c_k,d_k和a_k-1,b_k-1,c_k-1,d_k-1为分别在k时刻和k-1时刻的状态变量a,b,c,d的值，为在k-1时刻，独立的且分别对应状态变量a,b,c,d的噪声；

同时构建测量方程

HI_k＝a_k·exp(b_k·k)+c_k·exp(d_k·k)+v_k(8)

上式中，HI_k为在k时刻健康指数的测量值，v_k为在k时刻的测量噪声；

利用扩展卡尔曼滤波算法更新状态方程和测量方程参数至k时刻，按公式(3)计算k+l时刻的健康指数HI(k+l)

HI(k+l)＝a_k·exp(b_k·(k+l))+c_k·exp(d_k·(k+l))(9)

上式中，l＝1,2,…,∞；计算使得不等式(10)成立的l的值，并记录l的最小值为在k时刻预测的轴承剩余寿命；

HI(k+l)＞故障阀值(10)。

本实施例利用PRONOSTIA平台轴承全周期寿命数据对基于扩展卡尔曼滤波的轴承故障诊断及故障预测方法进行验证。具体过程如下：

(1)采集轴承的振动信号。通过加速度传感器采集水平方向和垂直方向的振动信号，信号每10s采集一次，每一次采集时长为0.1s。数据采样频率为25.6kHz；

(2)利用轴承健康时的振动信号构建AR模型，运用该AR模型对采集到的振动信号进行滤波，突出振动信号与故障相关的信息，从而有利于后续健康指数构建，故障诊断及预测工作；

(3)运用小波包变换对AR模型滤波后的残差信号进行分析，构建小波包系数对应的能量特征用于马氏距离计算；

(4)进行马氏距离的计算，基于此构建健康指数，对值为非负、非高斯分布的健康指数，运用Box-Cox变换，转换成高斯分布的数据，利用高斯分布的性质和逆Box-Cox变换，确定轴承发生异常时所对应的健康指数阈值。由附图5可以得知轴承的健康指数数据没有服从高斯分布，因此利用Box-Cox变换。如附图7所示，变换后的数据服从平均值为0.26，标准差为0.15的高斯分布，由此得到轴承发生异常时所对应的健康指数阈值为2.235；

(5)对耗损期的健康指数数据拟合分析，附图6为轴承耗损期数据用指数函数和两个指数函数之和的拟合曲线，从图中可以看出，利用双指数函数可以较好的拟合耗损期轴承数据。构建退化模型并建立状态空间模型，利用当前数据和扩展卡尔曼滤波算法更新模型参数并预测轴承剩余寿命。附图8显示了7组轴承数据在耗损期内的健康指数与时间的关系。通过计算，故障阈值为19.41，即健康指数超过该值时认定轴承是失效的。利用扩展卡尔曼滤波算法更新模型参数和预测剩余寿命，建立剩余寿命预测模型为：

HI(k+l)＝a_k·exp(b_k·(k+l))+c_k·exp(d_k·(k+l))(9)

上式中，l＝1,2,…,∞；计算使得不等式(10)成立的l的值，并记录l的最小值为在k时刻预测的轴承剩余寿命；

HI(k+l)＞故障阀值(10)。

附图9表示轴承数据的预测曲线，从曲线中可以看出，一开始由于数据不足，预测曲线与实际寿命曲线的误差较大，随着数据量的不断增大，最终的预测曲线与实际曲线相吻合。有效的验证了扩展卡尔曼滤波算法在轴承故障检测和故障预测中的可行性。