滚动轴承集成期望最大化和粒子滤波的寿命预测模型

摘要

滚动轴承集成期望最大化和粒子滤波的寿命预测模型，首先采用峭度指标对轴承健康状态进行实时监测，确定寿命预测起始时刻；当满足预测起始条件后，采用有效值对轴承剩余寿命进行预测；在预测阶段，采用期望最大化方法对模型参数进行评估，同时采用粒子滤波方法对轴承状态进行评估，通过对模型参数和轴承状态的准确评估，提高剩余寿命预测精度，本发明能够实现对模型参数和滚动轴承状态的准确评估，并且在滚动轴承寿命预测中表现出了比传统指数模型更好的预测效果。

说明书

技术领域

本发明涉及滚动轴承剩余寿命预测技术领域，具体涉及滚动轴承集成期 望最大化和粒子滤波的寿命预测模型。

背景技术

滚动轴承广泛应用于机械设备中，由于复杂多变的工况，滚动轴承故障 时有发生。为了保证设备的安全运行，传统定期维修策略需要投入大量人力 物力对滚动轴承进行定期检修，对存在故障或安全隐患的滚动轴承进行更换。 滚动轴承从发生故障到完全失效往往经历很长的衰退期，对存在早期故障的 滚动轴承进行更换势必会大大缩短滚动轴承的有效服役周期，造成资源的浪 费。采用预防性维修策略，对滚动轴承的剩余寿命进行有效预测，可以延长 滚动轴承的有效服役周期，提高经济效益。因此如何对滚动轴承进行有效剩 余寿命预测正受到国内外学者的广泛关注。

基于模型的滚动轴承剩余寿命预测方法，试图采用物理或数学模型对滚 动轴承的衰退趋势进行描述，并根据观测数据对模型参数进行评估，以预测 滚动轴承未来时刻健康衰退趋势和剩余寿命。美国普渡大学Gebraeel等人提 出的指数预测模型与滚动轴承加速衰退趋势相适应，因此在滚动轴承剩余寿 命预测中得到广泛应用。在指数预测模型中，模型参数和状态评估的准确性 是影响模型预测精度的两个关键因素。Gebraeel最先提出的指数预测模型采 用贝叶斯框架对模型参数进行评估。为了提高参数评估的准确性，清华大学 司小胜等人对指数模型进行改进，采用期望最大化和贝叶斯更新相结合的方 法对模型参数进行评估，得到了较好的参数评估效果。但是以上研究工作都 存在一个弊端，即直接将滚动轴承观测值作为状态值进行寿命预测，没有对 滚动轴承的真实健康状态进行评估。实际上，观测值与滚动轴承真实状态之 间存在一定差异，而且观测值中有大量随机噪声的干扰。因此，传统指数预 测模型只对模型参数进行评估而忽略了对滚动轴承状态的评估，导致模型预 测精度降低。

发明内容

为了克服上述现有技术存在的缺点，本发明的目的在于提供滚动轴承集 成期望最大化和粒子滤波的寿命预测模型，采用期望最大化方法对模型参数 进行评估，采用粒子滤波对滚动轴承状态进行评估，以提高指数模型的预测 精度。

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

滚动轴承集成期望最大化和粒子滤波的寿命预测模型，包括以下步骤：

1)实时监测并采集滚动轴承振动信号，计算滚动轴承峭度指标和有效值；

2)计算滚动轴承健康时刻峭度指标的均值μ和标准差σ，以确定健康状 态下峭度指标的3σ区间[μ-3σ,μ+3σ]；

3)判断滚动轴承连续l+1个时刻的峭度指标{m_p+i}_i＝0:l是否满足 {|m_p+i-μ|>3σ}_i＝0:l，如果满足该条件，则确定m_p所对应的时刻t_p为寿命预测的 起始时刻；

4)从寿命预测起始时刻开始，将滚动轴承振动信号有效值带入衰退模型：

其中，x_i为t_i时刻状态值，是已知常数，θ′，β′和σ为三个未知参数， σB(t_i)～N(0,σ²t_i)服从布朗运动，对上式求对数得到以下变形形式：

其中，θ＝ln(θ′)服从分布β＝β′-σ²/2服从分布且规 定以简化计算；

对模型参数进行初始化，得初始参数μ₀＝μ_θ,0，μ₁＝μ_β,0和${\sigma }_1^2={\sigma }_{\beta ,0}^2;$

5)从概率密度函数中进行随机采样，得到N_s个初始粒 子粒子权值为其中Δt＝t_j-t_j-1为时间间隔；

6)t_k时刻得到观测序列S_1:k＝{s₁,s₂,...,s_k}，在参数已知条件下，观测序列 的条件概率密度为：

$p\left(S_{1:k}\right|\theta ,\beta )={\left(\frac{1}{\sqrt{2\pi {\sigma }^2\mathrm{\Delta t}}}\right)}^k\exp [-\frac{{(s_1-\theta -\beta t_1)}^2}{2{\sigma }^2{t}_1}-\underset{j=2}{\overset{k}{\Sigma }}\frac{{(s_j-s_{j-1}-\mathrm{\beta \Delta t})}^2}{2{\sigma }^2\mathrm{\Delta t}}];---\left(3\right)$

根据贝叶斯理论得到参数θ和β的联合概率密度函数为：

$\left(\begin{array}{c}p(\theta ,\beta |S_{1:k})\propto p\left(S_{1:k}\right|\theta ,\beta )p(\theta ,\beta )\\ \propto \exp [-\frac{{(s_1-\theta -\beta t_1)}^2}{2{\sigma }^2{t}_1}-\underset{j=2}{\overset{k}{\Sigma }}\frac{{(s_j-s_{j-1}-\mathrm{\beta \Delta t})}^2}{2{\sigma }^2\mathrm{\Delta t}}]\exp [-\frac{{(\theta -{\mu }_0)}^2}{2{\sigma }_0^2}]\exp [-\frac{{(\beta -{\mu }_1)}^2}{2{\sigma }_1^2}]\\ \propto \frac{1}{2\pi {\sigma }_{\theta ,k}{\sigma }_{\beta ,k}\sqrt{1-{\rho }_k^2}}\exp [-\frac{1}{2(1-{\rho }_k^2)}(\frac{{(\theta -{\mu }_{\theta ,k})}^2}{{\sigma }_{\theta ,k}^2}-2{\rho }_k\frac{(\theta -{\mu }_{\theta ,k})(\beta -{\mu }_{\beta ,k})}{{\sigma }_{\theta ,k}{\sigma }_{\beta ,k}}\\ +\frac{{(\beta -{\mu }_{\beta ,k})}^2}{{\sigma }_{\beta ,k}^2}\left)\right]\end{array}\right)---\left(4\right)$

由此得到参数更新公式如下：

$\left(\begin{array}{c}{\mu }_{\theta ,k}=\frac{(s_1{\sigma }_0^2+{\mu }_0{\sigma }^2{t}_1)({\sigma }^2+{\sigma }_1^2{t}_k)-{\sigma }_0^2{t}_1(s_k{\sigma }_1^2+{\mu }_1{\sigma }^2-0.5{\sigma }^4)}{({\sigma }_0^2+{\sigma }^2{t}_1)({\sigma }_1^2{t}_k+{\sigma }^2)-{\sigma }_0^2{\sigma }_1^2{t}_1}\\ {\sigma }_{\theta ,k}^2=\frac{{\sigma }_0^2{\sigma }^2{t}_1({\sigma }^2+{\sigma }_1^2{t}_k)}{({\sigma }_0^2+{\sigma }^2{t}_1)({\sigma }_1^2{t}_k+{\sigma }^2)-{\sigma }_0^2{\sigma }_1^2{t}_1}\\ {\mu }_{\beta ,k}=\frac{(s_k{\sigma }_1^2+{\mu }_1{\sigma }^2-0.5{\sigma }^4)({\sigma }_0^2+{\sigma }^2{t}_1)-{\sigma }_1^2(s_1{\sigma }_0^2+{\mu }_0{\sigma }^2{t}_1)}{({\sigma }_0^2+{\sigma }^2{t}_1)({\sigma }_1^2{t}_k+{\sigma }^2)-{\sigma }_0^2{\sigma }_1^2{t}_1}\\ {\sigma }_{\beta ,k}^2=\frac{{\sigma }_1^2{\sigma }^2{t}_1({\sigma }_0^2+{\sigma }^2{t}_1)}{({\sigma }_0^2+{\sigma }^2{t}_1)({\sigma }_1^2{t}_k+{\sigma }^2)-{\sigma }_0^2{\sigma }_1^2{t}_1}\\ {\rho }_k=\frac{-{\sigma }_0{\sigma }_1\sqrt{t_1}}{\sqrt{({\sigma }_0^2+{\sigma }^2{t}_1)({\sigma }_1^2{t}_k+{\sigma }^2)}}\end{array}\right)---\left(5\right)$

采用以上公式对参数μ_θ,k,μ_β,k和进行更新；

7)为待评估的模型参数，计算似然函数：

$\left(\begin{array}{c}\ln \left[p(S_{1:k},\theta ,\beta |{\Theta }_k)\right]=\ln \left[p\left(S_{1:k}\right|\theta ,\beta ,{\Theta }_k)\right]+\ln \left[p(\theta ,\beta |{\Theta }_k)\right]\\ =-\frac{k}{2}\ln \left(\mathrm{\Delta t}\right)-\frac{k+2}{2}\ln \left(2\pi \right)-\frac{k}{2}\ln \left({\sigma }_k^2\right)-\frac{{(s_1-\theta -\beta t_1)}^2}{2{\sigma }_k^2{t}_1}-\underset{j=2}{\overset{k}{\Sigma }}\frac{{(s_j-s_{j-1}-\mathrm{\beta \Delta t})}^2}{2{\sigma }_k^2\mathrm{\Delta t}}\\ -\frac{1}{2}\ln \left({\sigma }_{0,k}^2\right)-\frac{1}{2}\ln \left({\sigma }_{1,k}^2\right)-\frac{{(\theta -{\mu }_{0,k})}^2}{2{\sigma }_{0,k}^2}-\frac{{(\beta -{\mu }_{1,k})}^2}{2{\sigma }_{1,k}^2}\end{array}\right)---\left(6\right)$

似然函数的期望值为：

其中为第i次评估的结果，计算满足条件 的参数为第i+1次评估的结果，

$\left(\begin{array}{c}{\hat{\sigma }}_k^{2(i+1)}=\frac{1}{k}\left(\frac{s_1^2-2s_1({\mu }_{\theta ,k}+{\mu }_{\beta ,k}{t}_1)+{\sigma }_{\theta ,k}^2+{\sigma }_{\beta ,k}^2+2t_1({\rho }_k{\sigma }_{\theta ,k}{\sigma }_{\beta ,k}+{\mu }_{\theta ,k}{\mu }_{\beta ,k})+t_1^2({\mu }_{\beta ,k}^2+{\sigma }_{\beta ,k}^2)}{t_1}\right)\\ +\underset{j=2}{\overset{k}{\Sigma }}\frac{{(s_j-s_{j-1})}^2-(s_j-s_{j-1})\mathrm{\Delta t}{\mu }_{\beta ,k}+{\left(\mathrm{\Delta t}\right)}^2({\mu }_{\beta ,k}^2+{\sigma }_{\beta ,k}^2)}{\mathrm{\Delta t}}\\ {\hat{\mu }}_{0,k}^{(i+1)}={\mu }_{\theta ,k},{\hat{\sigma }}_{0,k}^{2(i+1)}={\sigma }_{\theta ,k}^2,\\ {\hat{\mu }}_{1,k}^{(i+1)}={\mu }_{\beta ,k},{\hat{\sigma }}_{1,k}^{2(i+1)}={\sigma }_{\beta ,k}^2;\end{array}\right)---\left(8\right)$

8)将衰退模型改写为以下形式：

s_k＝s_k-1+βΔt+σW(Δt)，   (9)

其中W(Δt)＝B(t_k)-B(t_k-1)，由此得到重要密度函数为：

$p\left(s_k\right|s_{k-1})=\frac{1}{\sqrt{2\mathrm{\pi \Delta t}\left(\sqrt{{\sigma }_{\beta ,k}^2\mathrm{\Delta t}+{\sigma }_k^2}\right)}}\exp [-\frac{{(s_k-s_{k-1}-{\mu }_{\beta ,k}\mathrm{\Delta t})}^2}{2\mathrm{\Delta t}({\sigma }_{\beta ,k}^2\mathrm{\Delta t}+{\sigma }_k^2)}],---\left(10\right)$

从以上重要密度函数中进行重要性采样，得到粒子集

9)采用t_k时刻的观测值S_k对粒子权值进行更新，

$w_k^i=w_{k-1}^i\frac{1}{\sqrt{2\pi {\sigma }_k^2{t}_k}}\exp [-\frac{{(s_k-s_k^i)}^2}{2{\sigma }_k^2{t}_k}],w_k^i=\frac{w_k^i}{\underset{i=1}{\overset{N_s}{\Sigma }}{w}_k^i},---\left(11\right)$

采用下式计算有效粒子数，

${\tilde{N}}_{\mathrm{eff}}=\frac{1}{\underset{i=1}{\overset{N_s}{\Sigma }}{\left(w_k^i\right)}^2},---\left(12\right)$

如果有效粒子数小于阈值N_T，需要根据粒子权值大小进行重采样，得到新的 粒子集${\left\{{\alpha }_k^{i^{*}}\right\}}_{i=1:N_s}$使其满足$P({\alpha }_k^{i^{*}}={\alpha }_k^i)=w_k^i,$粒子权值重置为$w_k^i=1/N_s;$

10)采用粒子集对滚动轴承状态进行评估，

${\hat{s}}_k=\underset{i=1}{\overset{N_s}{\Sigma }}\left(w_k^i{s}_k^i\right),---\left(13\right)$

然后将状态评估结果带入下式对滚动轴承剩余寿命的概率密度函数进行 预测，

$f_{L_k|S_{1:k}}\left(l_k\right|S_{1:k})=\frac{\gamma -{\hat{s}}_k}{\sqrt{2\pi l_k^3({\sigma }_{\beta ,k}^2{l}_k+{\sigma }^2)}}\exp [-\frac{{(\gamma -{\hat{s}}_k-{\mu }_{\beta ,k}{l}_k)}^2}{2l_k({\sigma }_{\beta ,k}^2{l}_k+{\sigma }^2)}],l_k\ge 0,---\left(14\right)$

其中l_k为t_k时刻的剩余寿命，γ为滚动轴承失效阈值。

本发明利用期望最大化方法在参数评估方面的优势和粒子滤波方法在状 态评估方面的优势，将两种方法集成到指数预测模型中，对滚动轴承的剩余 寿命进行预测，克服了原有指数模型只注重参数评估而忽视了状态评估的弊 端，采用粒子滤波方法对滚动轴承的健康状态进评估，状态评估结果能更好 地反映滚动轴承的真实衰退趋势，通过同时对模型参数和滚动轴承状态的准 确评估，提高了剩余寿命的预测精度。

附图说明

图1为滚动轴承集成期望最大化和粒子滤波的寿命预测模型流程图。

图2为PRONOSTIA实验台整体结构图。

图3为两个测试滚动轴承全寿命期内的振动信号。

图4为两个测试滚动轴承峭度指标和预测起始时刻选择结果图。

图5为两个测试滚动轴承有效值和对应的预测起始时刻。

图6为采用改进指数模型对两个测试滚动轴承有效值从预测起始时刻进 行评估的结果。

图7为采用传统指数模型和改进指数模型对两个测试滚动轴承剩余寿命 从预测起始时刻开始进行预测的结果对比图。

具体实施方式

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

如图1所示，滚动轴承集成期望最大化和粒子滤波的寿命预测模型，包 括以下步骤：

1)实时监测并采集滚动轴承振动信号，计算滚动轴承峭度指标和有效值；

2)计算滚动轴承健康时刻峭度指标的均值μ和标准差σ，以确定健康状 态下峭度指标的3σ区间[μ-3σ,μ+3σ]；

3)判断滚动轴承连续l+1个时刻的峭度指标{m_p+i}_i＝0:l是否满足 {|m_p+i-μ|>3σ}_i＝0:l，如果满足该条件，则确定m_p所对应的时刻t_p为寿命预测的 起始时刻；

4)从寿命预测起始时刻开始，将滚动轴承振动信号有效值带入衰退模型：

其中，x_i为t_i时刻状态值，是已知常数，θ′，β′和σ为三个未知参数， σB(t_i)～N(0,σ²t_i)服从布朗运动，对上式求对数得到以下变形形式：

其中，θ＝ln(θ′)服从分布β＝β′-σ²/2服从分布且规 定以简化计算；

对模型参数进行初始化，得初始参数μ₀＝μ_θ,0，μ₁＝μ_β,0和${\sigma }_1^2={\sigma }_{\beta ,0}^2;$

5)从概率密度函数中进行随机采样，得到N_s个初始粒 子粒子权值为其中Δt＝t_j-t_j-1为时间间隔；

6)t_k时刻得到观测序列S_1:k＝{s₁,s₂,...,s_k}，在参数已知条件下，观测序列 的条件概率密度为：

$p\left(S_{1:k}\right|\theta ,\beta )={\left(\frac{1}{\sqrt{2\pi {\sigma }^2\mathrm{\Delta t}}}\right)}^k\exp [-\frac{{(s_1-\theta -\beta t_1)}^2}{2{\sigma }^2{t}_1}-\underset{j=2}{\overset{k}{\Sigma }}\frac{{(s_j-s_{j-1}-\mathrm{\beta \Delta t})}^2}{2{\sigma }^2\mathrm{\Delta t}}];---\left(3\right)$

根据贝叶斯理论得到参数θ和β的联合概率密度函数为：

$\left(\begin{array}{c}p(\theta ,\beta |S_{1:k})\propto p\left(S_{1:k}\right|\theta ,\beta )p(\theta ,\beta )\\ \propto \exp [-\frac{{(s_1-\theta -\beta t_1)}^2}{2{\sigma }^2{t}_1}-\underset{j=2}{\overset{k}{\Sigma }}\frac{{(s_j-s_{j-1}-\mathrm{\beta \Delta t})}^2}{2{\sigma }^2\mathrm{\Delta t}}]\exp [-\frac{{(\theta -{\mu }_0)}^2}{2{\sigma }_0^2}]\exp [-\frac{{(\beta -{\mu }_1)}^2}{2{\sigma }_1^2}]\\ \propto \frac{1}{2\pi {\sigma }_{\theta ,k}{\sigma }_{\beta ,k}\sqrt{1-{\rho }_k^2}}\exp [-\frac{1}{2(1-{\rho }_k^2)}(\frac{{(\theta -{\mu }_{\theta ,k})}^2}{{\sigma }_{\theta ,k}^2}-2{\rho }_k\frac{(\theta -{\mu }_{\theta ,k})(\beta -{\mu }_{\beta ,k})}{{\sigma }_{\theta ,k}{\sigma }_{\beta ,k}}\\ +\frac{{(\beta -{\mu }_{\beta ,k})}^2}{{\sigma }_{\beta ,k}^2}\left)\right]\end{array}\right)---\left(4\right)$

由此得到参数更新公式如下：

$\left(\begin{array}{c}{\mu }_{\theta ,k}=\frac{(s_1{\sigma }_0^2+{\mu }_0{\sigma }^2{t}_1)({\sigma }^2+{\sigma }_1^2{t}_k)-{\sigma }_0^2{t}_1(s_k{\sigma }_1^2+{\mu }_1{\sigma }^2-0.5{\sigma }^4)}{({\sigma }_0^2+{\sigma }^2{t}_1)({\sigma }_1^2{t}_k+{\sigma }^2)-{\sigma }_0^2{\sigma }_1^2{t}_1}\\ {\sigma }_{\theta ,k}^2=\frac{{\sigma }_0^2{\sigma }^2{t}_1({\sigma }^2+{\sigma }_1^2{t}_k)}{({\sigma }_0^2+{\sigma }^2{t}_1)({\sigma }_1^2{t}_k+{\sigma }^2)-{\sigma }_0^2{\sigma }_1^2{t}_1}\\ {\mu }_{\beta ,k}=\frac{(s_k{\sigma }_1^2+{\mu }_1{\sigma }^2-0.5{\sigma }^4)({\sigma }_0^2+{\sigma }^2{t}_1)-{\sigma }_1^2(s_1{\sigma }_0^2+{\mu }_0{\sigma }^2{t}_1)}{({\sigma }_0^2+{\sigma }^2{t}_1)({\sigma }_1^2{t}_k+{\sigma }^2)-{\sigma }_0^2{\sigma }_1^2{t}_1}\\ {\sigma }_{\beta ,k}^2=\frac{{\sigma }_1^2{\sigma }^2{t}_1({\sigma }_0^2+{\sigma }^2{t}_1)}{({\sigma }_0^2+{\sigma }^2{t}_1)({\sigma }_1^2{t}_k+{\sigma }^2)-{\sigma }_0^2{\sigma }_1^2{t}_1}\\ {\rho }_k=\frac{-{\sigma }_0{\sigma }_1\sqrt{t_1}}{\sqrt{({\sigma }_0^2+{\sigma }^2{t}_1)({\sigma }_1^2{t}_k+{\sigma }^2)}}\end{array}\right)---\left(5\right)$

采用以上公式对参数μ_θ,k,μ_β,k和进行更新；

7)为待评估的模型参数，计算似然函数：

$\left(\begin{array}{c}\ln \left[p(S_{1:k},\theta ,\beta |{\Theta }_k)\right]=\ln \left[p\left(S_{1:k}\right|\theta ,\beta ,{\Theta }_k)\right]+\ln \left[p(\theta ,\beta |{\Theta }_k)\right]\\ =-\frac{k}{2}\ln \left(\mathrm{\Delta t}\right)-\frac{k+2}{2}\ln \left(2\pi \right)-\frac{k}{2}\ln \left({\sigma }_k^2\right)-\frac{{(s_1-\theta -\beta t_1)}^2}{2{\sigma }_k^2{t}_1}-\underset{j=2}{\overset{k}{\Sigma }}\frac{{(s_j-s_{j-1}-\mathrm{\beta \Delta t})}^2}{2{\sigma }_k^2\mathrm{\Delta t}}\\ -\frac{1}{2}\ln \left({\sigma }_{0,k}^2\right)-\frac{1}{2}\ln \left({\sigma }_{1,k}^2\right)-\frac{{(\theta -{\mu }_{0,k})}^2}{2{\sigma }_{0,k}^2}-\frac{{(\beta -{\mu }_{1,k})}^2}{2{\sigma }_{1,k}^2}\end{array}\right)---\left(6\right)$

似然函数的期望值为：

其中为第i次评估的结果，计算满足条件 的参数为第i+1次评估的结果，

$\left(\begin{array}{c}{\hat{\sigma }}_k^{2(i+1)}=\frac{1}{k}\left(\frac{s_1^2-2s_1({\mu }_{\theta ,k}+{\mu }_{\beta ,k}{t}_1)+{\sigma }_{\theta ,k}^2+{\sigma }_{\beta ,k}^2+2t_1({\rho }_k{\sigma }_{\theta ,k}{\sigma }_{\beta ,k}+{\mu }_{\theta ,k}{\mu }_{\beta ,k})+t_1^2({\mu }_{\beta ,k}^2+{\sigma }_{\beta ,k}^2)}{t_1}\right)\\ +\underset{j=2}{\overset{k}{\Sigma }}\frac{{(s_j-s_{j-1})}^2-(s_j-s_{j-1})\mathrm{\Delta t}{\mu }_{\beta ,k}+{\left(\mathrm{\Delta t}\right)}^2({\mu }_{\beta ,k}^2+{\sigma }_{\beta ,k}^2)}{\mathrm{\Delta t}}\\ {\hat{\mu }}_{0,k}^{(i+1)}={\mu }_{\theta ,k},{\hat{\sigma }}_{0,k}^{2(i+1)}={\sigma }_{\theta ,k}^2,\\ {\hat{\mu }}_{1,k}^{(i+1)}={\mu }_{\beta ,k},{\hat{\sigma }}_{1,k}^{2(i+1)}={\sigma }_{\beta ,k}^2;\end{array}\right)---\left(8\right)$

8)将衰退模型改写为以下形式：

s_k＝s_k-1+βΔt+σW(Δt)，   (9)

其中W(Δt)＝B(t_k)-B(t_k-1)，由此得到重要密度函数为：

$p\left(s_k\right|s_{k-1})=\frac{1}{\sqrt{2\mathrm{\pi \Delta t}\left(\sqrt{{\sigma }_{\beta ,k}^2\mathrm{\Delta t}+{\sigma }_k^2}\right)}}\exp [-\frac{{(s_k-s_{k-1}-{\mu }_{\beta ,k}\mathrm{\Delta t})}^2}{2\mathrm{\Delta t}({\sigma }_{\beta ,k}^2\mathrm{\Delta t}+{\sigma }_k^2)}],---\left(10\right)$

从以上重要密度函数中进行重要性采样，得到粒子集

9)采用t_k时刻的观测值S_k对粒子权值进行更新，

$w_k^i=w_{k-1}^i\frac{1}{\sqrt{2\pi {\sigma }_k^2{t}_k}}\exp [-\frac{{(s_k-s_k^i)}^2}{2{\sigma }_k^2{t}_k}],w_k^i=\frac{w_k^i}{\underset{i=1}{\overset{N_s}{\Sigma }}{w}_k^i},---\left(11\right)$

采用下式计算有效粒子数，

${\tilde{N}}_{\mathrm{eff}}=\frac{1}{\underset{i=1}{\overset{N_s}{\Sigma }}{\left(w_k^i\right)}^2},---\left(12\right)$

如果有效粒子数小于阈值N_T，需要根据粒子权值大小进行重采样，得到新的 粒子集${\left\{{\alpha }_k^{i^{*}}\right\}}_{i=1:N_s}$使其满足$P({\alpha }_k^{i^{*}}={\alpha }_k^i)=w_k^i,$粒子权值重置为$w_k^i=1/N_s;$

10)采用粒子集对滚动轴承状态进行评估，

${\hat{s}}_k=\underset{i=1}{\overset{N_s}{\Sigma }}\left(w_k^i{s}_k^i\right),---\left(13\right)$

然后将状态评估结果带入下式对滚动轴承剩余寿命的概率密度函数进行 预测，

$f_{L_k|S_{1:k}}\left(l_k\right|S_{1:k})=\frac{\gamma -{\hat{s}}_k}{\sqrt{2\pi l_k^3({\sigma }_{\beta ,k}^2{l}_k+{\sigma }^2)}}\exp [-\frac{{(\gamma -{\hat{s}}_k-{\mu }_{\beta ,k}{l}_k)}^2}{2l_k({\sigma }_{\beta ,k}^2{l}_k+{\sigma }^2)}],l_k\ge 0,---\left(14\right)$

其中l_k为t_k时刻的剩余寿命，γ为滚动轴承失效阈值。

下面结合实施例对本发明做详细描述。

实施例：采用PRONOSTIA实验台上采集的滚动轴承加速寿命实验数据 对本发明进行验证。

PRONOSTIA实验台如图2所示，该实验台通过对滚动轴承气压加载， 使滚动轴承在高负荷条件下工作，可以在数小时内实现滚动轴承从正常状态 退化到完全失效。实验过程中，滚动轴承转速为1800rpm，负载为4kN。采 用加速度传感器对滚动轴承振动信号进行采样，采样频率为25.6kHz，数据 长度为2560，每次采样持续时间为0.1s，采样间隔为10s。当振动幅值超过 20m/s²时，滚动轴承完全失效。两组实验滚动轴承的全寿命振动信号如图3 所示。

从振动信号中提取峭度指标，对寿命预测起始时刻进行选择，选择过程 如图4所示，由健康历史数据确定滚动轴承峭度指标的3σ区间，然后对滚动 轴承峭度指标进行实时监测，当峭度指标m_p超出3σ区间即|m_p-μ|>3σ，说 明此时滚动轴承状态异常，滚动轴承异常状态可能由滚动轴承故障引起，也 可能由随机噪声引起，为了排除随机噪声对峭度指标带来的干扰，连续观测t_p时刻之后的l个时刻的峭度指标{m_p+i}_i＝1:l是否也超出了3σ区间，即判断 {|m_p+i-μ|>3σ}_i＝1:l，如果满足以上条件，说明滚动轴承的异常状态确实由故障 引起。因此选择t_p时刻为寿命预测的起始时刻。

确定预测起始时刻之后，从该时刻开始，将滚动轴承有效值带入衰退模 型，对滚动轴承剩余寿命进行预测，两滚动轴承的有效值如图5所示，可以 看出，所选的时刻与滚动轴承加速衰退的起始时刻相近，作为滚动轴承的预 测起始时刻较为合适，本发明采用粒子滤波方法对滚动轴承有效值进行评估， 评估结果如图6中所示，从图中可以看出，由于随机噪声的影响，有效值局 部波动性较强，而粒子滤波方法评估结果可以有效抑制随机噪声的干扰，反 映出滚动轴承真实的加速衰退趋势，因此本发明通过采用粒子滤波方法，实 现了对滚动轴承健康状态的有效评估。

图7显示了滚动轴承剩余寿命预测结果，从图中可以看出，在寿命预测 前期，传统指数模型和改进指数模型的预测效果都不理想，随着时间的推移， 两种方法都逐渐收敛到真实值，但是改进指数模型的收敛速度更快，预测精 度更高，其原因为，在刚开始预测时，由于观测数据点数较少，模型参数无 法得到有效评估，因此引起较大的预测误差，当观测点数逐渐增多，模型参 数获得较准确的评估结果，此时状态评估结果便成为影响寿命预测结果的主 要因素，由于改进指数模型采用粒子滤波方法实现了对滚动轴承健康状态的 准确评估，因此其预测效果明显好于传统指数模型。

通过采用滚动轴承加速寿命实验数据验证了本发明在滚动轴承剩余寿命 预测中的优势，由于传统指数模型只是对模型参数进行评估，而滚动轴承观 测值中包含大量的随机噪声，影响了剩余预测结果的精度，本发明对传统指 数模型进行改进，采用粒子滤波方法对滚动轴承健康状态进行有效评估，从 而提高了指数模型对滚动轴承剩余寿命的预测精度。

本发明所提出的滚动轴承集成期望最大化和粒子滤波的寿命预测模型， 并不只局限于滚动轴承的剩余寿命预测，还可以应用于其他机械电子产品的 剩余寿命预测问题。大量研究工作证明，本方法适用于各类具有指数衰退形 式的机电产品的剩余寿命预测。实施者只需对本方法相应步骤进行适当调整， 以适应不同产品的应用需求。另外，本发明提供了一种同时对模型参数和产 品健康状态进行准确评估的思路，应当指出，在不脱离本发明构思的前提下， 所做的调整和变形，也应视为本发明的保护范围。