基于特征融合和粒子滤波的滚动轴承剩余寿命预测方法

摘要

基于特征融合和粒子滤波的滚动轴承剩余寿命预测方法，在指标计算部分，首先从轴承振动信号中提取原始特征，采用相关性聚类方法对提取的原始特征进行聚类，然后从各类中选取一个典型特征构成最优特征集，最后采用加权融合方法将特征集融合为最终的衰退指标；在寿命预测部分，先对衰退指标进行平滑处理和重采样，将时间间隔调整为期望值，采用最小二乘拟合计算状态空间模型初始参数，然后根据新的观测数据对模型参数进行实时更新，最后预测轴承剩余寿命，本发明寿命预测结果与真实值差异小，应用效果好。

说明书

技术领域

本发明属于滚动轴承剩余寿命预测技术领域，具体涉及一种基于特征融合和粒子滤波的滚动轴承剩余寿命预测方法。 

背景技术

滚动轴承在旋转机械中广泛应用，其健康状况直接关系到机械设备的安全运行。由于滚动轴承经常在高速重载的恶劣环境下工作，磨损、疲劳点蚀等故障时有发生。一旦轴承出现故障，势必会对设备的安全运行造成严重威胁，轻则引起设备停机的生产事故，重则导致机毁人亡的重大灾难。由于各滚动轴承的有效寿命差异极大，传统的定期维修策略不仅费时费力，而且可靠性不高，无法满足工程的实际需求。预防性维修策略可以对滚动轴承的健康状况进行准确评估，并对其剩余寿命进行有效预测，保证在最小的维修成本下实现最大化的生产效率，因此越来越受到人们的关注。对滚动轴承的剩余寿命进行有效预测可为维修决策提供关键依据，是预防性维修的基础和前提。 

寿命预测的两个关键技术环节就是指标提取和趋势预测。 

1）指标提取就是从轴承的振动信号中提取出包含故障信息的特征指标，对轴承的健康状况进行评估。常用的特征指标有RMS、峰值、峭度、波形指标以及小波包能量指标等。这些特征指标往往对故障的某个阶段较为敏感，例如RMS只在故障严重阶段才会呈现出明显的变化趋势，而峭度和波形指标在故障初期波动较大，随着故障的发展会逐渐趋于平缓。这些特征只能片面地反映故障的局部信息，无法对轴承的健康状况进行全面描述。因此，要想实现对轴承健康状况的准确评估，就必须从振动信号中提取出能够全面反 应轴承故障信息的特征指标。 

2）趋势预测就是采用已知数据序列训练预测模型，对其参数进行迭代更新，使其能够包含指标的趋势信息，然后对未来发展趋势进行预测。粒子滤波方法以其在处理非线性非高斯问题方面的优势，已被广泛应用于机电设备的趋势预测领域。粒子滤波方法中，系统在t_k=k·△t时刻的状态用以下状态空间模型描述： 

$\left(\begin{array}{c}{x}_k=f_k(x_{k-1},{\theta }_k,{\omega }_k)\\ y_k=g_k(x_k,v_k)\end{array}\right)$

其中，x_k是系统在t_k时刻的状态； 

f_k：是状态传递方程； 

g_k：是观测方程； 

θ_k是由模型参数组成的参数向量； 

y_k是由振动信号计算得到的特征指标； 

ω_k，ν_k分别代表过程噪声和观测噪声引起的噪声误差。 

假设初始参数已知，首先进行单步预测得到先验概率：p(x_k|y_1:k-1)=∫p(x_k|x_k-1)p(x_k-1|y_1:k-1)dx_k-1。然后代入新的观测数据对上式进行更新得到后验概率：$p\left(x_k\right|y_{1:k})=\frac{p\left(y_k\right|x_k)p\left(x_k\right|y_{1:k-1})}{p\left(y_k\right|y_{1:k-1})}.$粒子滤波采用蒙特卡罗模拟方法，将状态空间描述为一系列具有相同概率分布的粒子以及相应权值序列：$p\left(x_k\right|y_{1:k})\approx \underset{i=1}{\overset{N}{\Sigma }}{\tilde{w}}_k\left(x_{0:k}^i\right)·\delta (x_{0:k}-x_{0:k}^i).$采用序贯重要性采样算法对式中的粒子权值进行更新：$w_k=w_{k-1}\frac{p\left(y_k\right|x_k)p\left(x_k\right|x_{k-1})}{q\left(x_k\right|x_{0:k-1},y_{1:k})}.$为了进一步简化计算过程，提高计算效率，式中q(x_k|x_0:k-1,y_1:k)选为p(x_k|x_k-1)，上式简化为：w_k=w_k-1·p(y_k|x_k)。 

粒子滤波方法能够充分利用观测数据对状态空间模型的参数进行迭代更 新，不断调整采样粒子及其权值，并通过状态传递方程对未来时刻的系统状态进行预测，达到寿命预测的目的。但是，在状态空间模型初始参数的选择上，目前大都采用根据已有的历史数据进行人为设定的方法，这就要求不同样本之间模型的参数必须相差不大。然而，滚动轴承作为机械传动系统中经常使用的零部件，通常在载荷、转速大幅度波动的恶劣环境下工作。复杂的运行工况，多变的外部环境使得各轴承状态模型的初始参数变化极大，人为设定的模型参数往往无法有效描述其衰退趋势，使得寿命预测结果与真实值差异较大，应用效果不佳，限制了粒子滤波方法在滚动轴承剩余寿命预测中的应用。亟需寻找有效的模型初始参数选择方法，对该问题进行解决。 

发明内容

为了克服上述现有技术的缺陷，本发明提供了一种基于特征融合和粒子滤波的滚动轴承剩余寿命预测方法，寿命预测结果与真实值差异小，应用效果好。 

为了达到上述目的，本发明采取的技术方案为： 

基于特征融合和粒子滤波的滚动轴承剩余寿命预测方法，包含以下步骤： 

第1步，从轴承振动信号中提取M个原始特征，根据各特征之间的相关性系数大小，采用相关性聚类方法对提取的原始特征进行聚类，聚类原则为：类内特征相关性最大，类间特征相关性最小，聚类方法具体过程如下： 

1.1计算M个原始特征的相关系数矩阵，并初始化聚类数K； 

1.2选择相关系数最小的两个特征为第一类和第二类的中心，然后选择与已有类中心平均相关系数最小的特征为下一类的中心，直到选出K个类中心为止； 

1.3将剩余的M-K个特征依次归入与其平均相关系数最大的类中； 

第2步，防止同类特征间的信息冗余，从各类中选取一个典型特征构成最优特征集，采用加权融合方法将特征集融合为最终的衰退指标，过程如下： 

2.1计算各特征的趋势性指标，计算公式为 

$T=\frac{N\underset{t=1}{\overset{N}{\Sigma }}(t·F_t)-\underset{t=1}{\overset{N}{\Sigma }}t·\underset{t=1}{\overset{N}{\Sigma }}{F}_t}{\sqrt{[N\underset{t=1}{\overset{N}{\Sigma }}{t}^2-{\left(\underset{t=1}{\overset{N}{\Sigma }}\right)]}^2[N\underset{t=1}{\overset{N}{\Sigma }}{F}_t^2-{\left(\underset{t=1}{\overset{N}{\Sigma }}{F}_t\right)}^2]}}$

其中，N是原始数据长度，F_t是第t次采样的特征指标，趋势性指标取值范围为：-1≤T≤+1； 

2.2从各类特征中选择趋势性指标最大的特征组成最优特征集； 

2.3归一化特征值，计算正常状态空间的特征均值矢量： 

$V_i=\frac{1}{p}\underset{t=1}{\overset{p}{\Sigma }}{\bar{F}}_{i,t}(i=1,...,K)$

其中，p为正常状态空间维数； 

2.4归一化最优特征集的趋势性指标； 

$W_i=\frac{\left|T_i\right|}{\underset{i=1}{\overset{K}{\Sigma }}\left|T_i\right|}(i=1,...,K)$

2.5计算归一化特征集与正常状态空间均值V的曼哈顿距离，并用归一化趋势性指标进行加权，得最终的衰退指标： 

$D_t=\underset{i=1}{\overset{K}{\Sigma }}(W_i·{|\bar{F}}_{i,t}-V_i|)(t=p+1,...,N);$

第3步，对衰退指标进行平滑处理和重采样，以消除噪声影响，将时间间隔调整为期望值；再采用最小二乘拟合对已知数据序列进行曲线拟合，得到状态空间模型初始参数；然后根据新的观测数据对模型参数进行实时更新；最后对未来时刻的轴承状态进行预测，统计各粒子达到失效阈值的时间，计算剩余寿命概率分布情况。 

相关性聚类方法根据特征值之间的相关性对原始特征进行聚类，采用“类间特征相关性最小”原则选择各类中心，采用“类内特征相关性最大”原则进行特征归类，保证了聚类结果的类间距最大，类内距最小。通过该方法聚类之后，具有相同性质的特征指标被聚为同一类，不同类的特征之间分别反 映了轴承不同的故障信息，既保证了故障信息的全面性，又有效地防止了信息冗余。 

特征融合方法采用趋势性指标对各原始特征进行评估，并选择出趋势性最强的特征指标组成最优特征集，然后采用加权融合算法，得到最后的衰退指标。计算得到的衰退指标，既包含了所有原始特征的故障信息，又具有较强的趋势性，克服了原始特征无法全面描述轴承健康状况的问题，是较为理想的轴承寿命预测指标。 

本发明在寿命预测模块采用最小二乘拟合方法首先对轴承状态空间模型参数进行初始化，克服了粒子滤波方法模型初始参数人为设定，无法满足工程应用需求的问题。通过采用IEEE PHM2012挑战赛数据对本方法进行验证，证明本方法可以自适应确定模型初始参数，对轴承衰退趋势进行合理描述，将预测结果与自适应神经模糊推理系统（ANFIS）进行比较，发现该方法预测效果较好。 

附图说明

图1为本发明的滚动轴承剩余寿命预测方法流程图。 

图2为相关性聚类算法流程图。 

图3为PRONOSTIA实验台结构图。 

图4为训练轴承和测试轴承水平、竖直两个方向的振动信号时域波形图。 

图5为从测试轴承原始特征中选择出来的最优特征集。 

图6为最终计算得到的训练轴承和测试轴承的衰退指标。 

图7为模型参数的更新过程图。 

图8为采用基于特征融合和粒子滤波方法和基于自适应神经模糊推理系统方法对测试轴承进行寿命预测的结果对比图。 

具体实施方式

下面结合附图和实施例对本发明做进一步详细描述： 

如图1所示，基于特征融合和粒子滤波的滚动轴承剩余寿命预测方法，包括以下步骤： 

第1步，从轴承振动信号中提取M个原始特征。为了能够充分挖掘振动信号中隐含的故障信息，原始特征可以从时域、频谱和时频域中同时提取，轴承振动信号可以是来自不同测点的振动数据，以保证对设备健康状况的全面监测。然后采用相关性聚类方法对提取的原始特征进行聚类，聚类原则为：类内特征相关性最大，类间特征相关性最小，如图2所示，聚类方法具体过程如下： 

1.1计算M个原始特征的相关系数矩阵，并初始化聚类数K； 

1.2选择相关系数最小的两个特征为第一类和第二类的中心，然后选择与已有类中心平均相关系数最小的特征为下一类的中心，直到选出K个类中心为止； 

1.3将剩余的M-K个特征依次归入与其平均相关系数最大的类中； 

第2步，防止同类特征间的信息冗余，从各类中选取一个典型特征构成最优特征集，采用加权融合方法将特征集融合为最终的衰退指标，过程如下： 

2.1计算各特征的趋势性指标，计算公式为 

$T=\frac{N\underset{t=1}{\overset{N}{\Sigma }}(t·F_t)-\underset{t=1}{\overset{N}{\Sigma }}t·\underset{t=1}{\overset{N}{\Sigma }}{F}_t}{\sqrt{[N\underset{t=1}{\overset{N}{\Sigma }}{t}^2-{\left(\underset{t=1}{\overset{N}{\Sigma }}\right)]}^2[N\underset{t=1}{\overset{N}{\Sigma }}{F}_t^2-{\left(\underset{t=1}{\overset{N}{\Sigma }}{F}_t\right)}^2]}}$

其中，N是原始数据长度，F_t是第t次采样的特征指标，趋势性指标取值范围为：-1≤T≤+1； 

2.2从各类特征中选择趋势性指标最大的特征组成最优特征集； 

2.3归一化特征值，计算正常状态空间的特征均值矢量： 

$V_i=\frac{1}{p}\underset{t=1}{\overset{p}{\Sigma }}{\bar{F}}_{i,t}(i=1,...,K)$

其中，p为正常状态空间维数； 

2.4归一化最优特征集的趋势性指标； 

$W_i=\frac{\left|T_i\right|}{\underset{i=1}{\overset{K}{\Sigma }}\left|T_i\right|}(i=1,...,K)$

2.5计算归一化特征集与正常状态空间均值V的曼哈顿距离，并用归一化趋势性指标进行加权，得最终的衰退指标： 

$D_t=\underset{i=1}{\overset{K}{\Sigma }}(W_i·{|\bar{F}}_{i,t}-V_i|)(t=p+1,...,N);$

第3步，对衰退指标进行平滑处理和重采样，以消除噪声影响，将时间间隔调整为期望值；再采用最小二乘拟合计算状态空间模型初始参数；然后根据新的观测数据对模型参数进行实时更新；最后对轴承剩余寿命进行预测。具体为： 

3.1采用“loess”滤波器对衰退指标进行平滑处理，以消除噪声误差带来的影响，然后按照时间间隔△t的要求，对指标进行重采样； 

3.2采用Paris-Erdogan模型对轴承状态进行描述： 

$\frac{\mathrm{dx}}{\mathrm{dN}}=C{\left(\mathrm{\Delta k}\right)}^n$

其中，x表示轴承损伤程度，如裂纹长度、剥落面积等，N表示应力循环次数，C和n是与材料特性有关的常数，△k表示应力强度因数，计算公式为根据粒子滤波状态空间模型的描述方法对以上模型进行变形得到以下形式： 

$\left(\begin{array}{c}{x}_k=x_{k-1}+A·x_{k-1}^m·\mathrm{\Delta t}\\ y_k=x_k+v_k\end{array}\right)$

其中，A和m是需要计算的模型参数，ω_k由于模型参数的不确定性而被忽略，假设y_k与x_k之间呈线性关系，是由于运行环境和设备本身所引起的测量噪声，得到轴承状态空间模型之后，将已知序列代入状态传递方程，采用最小二乘拟合方法对参数A和m进行求解； 

3.3由状态空间模型表达式，计算权值更新系数为： 

$p\left(y_k\right|x_k)=\frac{1}{\sqrt{2\pi }{\sigma }_v}{e}^{-\frac{1}{2}{\left(\frac{y_k-x_{k-1}-A·x_{k-1}^m·\mathrm{\Delta t}}{{\sigma }_v}\right)}^2}$

根据新的观测数据对粒子权值进行实时更新w_k=w_k-1·p(y_k|x_k)，权值更新之后对粒子进行重采样得到新的粒子集合表示k时刻第i个粒子，通过粒子和权值的更新，不断修正模型参数的概率分布； 

3.4参数更新完毕之后，保持模型参数不变，采用状态传递方程对未来时刻的轴承状态进行预测，统计各粒子达到失效阈值的时间，计算剩余寿命概率分布情况。 

实施例：采用PRONOSTIA实验台上采集的滚动轴承加速寿命实验数据验证本发明方法的正确性。 

PRONOSTIA实验台如图3所示，该实验台专门设计用于滚动轴承故障诊断、趋势预测方法的验证。实验台有三部分组成：传动机构、加载部分和数采系统。实验轴承转速为1800rpm，负载为4000N。采样频率为25.6kHz，数据长度为2560，每次采样持续时间为0.1s，采样间隔为10s，安装在轴承座水平和竖直两个方向的加速度传感器同时采样。轴承从正常开始，直到完全失效为止。实验数据包括两组轴承，一个训练轴承和一个测试轴承，其振动信号波形如图4所示，其中（a）和（b）为训练轴承的水平和竖直振动信号，（c）和（d）为测试轴承的水平和竖直振动信号。当振动幅值超过20m/s²时，轴承完全失效。两个轴承的实际使用寿命为27600s和17790s。 

为了充分挖掘振动信号的故障信息，分别从各测点振动信号中提取2个双曲函数指标，19个时域统计特征和16个小波包能量指标。2个双曲函数指标的计算公式如下： 

$\sigma \left(\log \right[x_i+{(x_i^2-1)}^{1/2}\left]\right)$

$\sigma \left(\log \right[x_i+{(x_i^2+1)}^{1/2}\left]\right)$

19个时域统计特征如表1所示。16个小波包能量指标是对原始信号进行三层小波包分解，所对应的8个频段的节点能量和能量比指标。每组振动信号中提取37个特征，总共提取74个特征，用于后续的特征聚类和融合。 

表1时域统计特征 

采用相关性聚类算法对测试轴承的74组特征指标进行聚类，聚类数设为8，然后从中提取出趋势性最强的特征组成最优特征集。测试轴承所对应的最优特征集如图5所示。由图可以看出，8组特征分别具有不同的特性，包含了轴承不同的故障信息。水平偏斜度、水平峰峰值和竖直1频段能量比3个指标对故障初期不敏感，当达到严重失效阶段时才呈现出急速上升趋势。另外，三者之间也存在差异，水平峰峰值较水平偏斜度在前期有明显的噪声波动，而竖直1频段能量比在后期有明显的噪声波动。水平1频段能量比对故障初期变化较为敏感，对故障严重期反而敏感性不强。竖直均值在故障前期一直呈现平稳波动，直到严重失效期波动幅值明显增大，而水平均值在整个运行阶段一直保持平稳波动。竖直8频段能量比在故障严重期有明显的下降趋势，而水平3频段能量比在故障前期存在幅值下降趋势，到故障严重期幅值反而有所升高。从以上各特征趋势的分析可知，各个特征在故障的不同阶段呈现出不同的变化趋势，有些特征对故障初期较为敏感，有些特征对故障 严重期较为敏感，有些特征甚至不包含任何故障信息。 

下面采用加权融合算法对最优特征集进行特征融合，训练轴承和测试轴承融合后的衰退指标分别如图6（a）、（b）所示。由图6（b）可以看出，测试轴承衰退趋势大致可以分为三个阶段：正常阶段、故障发展期和严重阶段。在正常阶段衰退指标平稳波动。故障发展期轴承衰退指标缓慢升高，说明在这段时期滚动轴承故障逐渐扩展。到严重阶段，幅值突然增大，并呈现大幅度波动趋势，说明此时轴承故障急剧扩展，由故障引起的振动噪声增加。对比图5和图6（b）可以发现，融合后的衰退指标不仅包含原始特征所有的故障信息，而且趋势性也比原始特征要好，滚动轴承故障演变趋势能够通过衰退指标明显地反映出来，是较为理想的趋势预测指标。 

得到轴承衰退指标之后，对指标进行平滑处理，然后将时间间隔设定为△t=1min，对指标进行重采样。假设已知数据时间段为0～170min，采用已知数据序列对模型初始参数进行求解得A=0.0100，m=0.9282。随机产生5000个粒子，使用观测数据对粒子分布和权值进行不断更新，从而实现对模型参数的更新，模型参数均值最终更新为A=0.0187，m=1.2710。参数更新过程如图7所示。 

要想实现对滚动轴承的剩余寿命预测，必须设定一个失效阈值。本实例测试轴承的失效阈值根据训练轴承衰退指标的最大幅值设定。根据图6（a），训练轴承的衰退指标最大值约为0.6，因此测试轴承失效阈值同样设定为0.6。保持参数A和m不变，对测试轴承的衰退趋势进行预测，对各粒子失效时间进行统计，从而得到轴承剩余寿命分布。为了验证粒子滤波（PF）方法的有效性，采用自适应神经模糊推理系统（ANFIS）对轴承剩余寿命进行预测，将两种方法的预测结果进行对比，对比结果如图8所示。测试轴承的真实使用寿命为17790s，时间间隔设定为△t=1min之后，使用寿命换算为296min。在t=170min时刻，测试轴承的真实剩余寿命为126min。粒子滤波方法的预测结果为116min，相对误差为-8.62%。自适应神经模糊推理系统的预测结果 为96min，相对误差为-23.81%。粒子滤波方法预测结果要好于自适应神经模糊推理系统的预测结果。 

通过以上实施例具体处理过程分析以及实验结果对比可以发现，本发明提出的基于特征融合和粒子滤波的滚动轴承剩余寿命预测方法，可以充分利用轴承振动信号，提取有效趋势预测指标，对模型初始参数进行准确计算，并能实现对轴承剩余寿命的准确评估。该方法与自适应神经模糊推理系统相比，预测精度较高。 

本发明所提出的基于特征融合和粒子滤波的寿命预测方法，并不只局限于滚动轴承的剩余寿命预测，还可以应用于其他机械电子产品的剩余寿命预测问题，实施者只需对本方法相应步骤进行适当调整，以适应不同产品的应用需求。应当指出，在不脱离本发明构思的前提下，所做的调整和变形，也应视为本发明的保护范围。