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

摘要

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

主权项

滚动轴承集成期望最大化和粒子滤波的寿命预测模型，其特征在于，包括以下步骤：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)$根据贝叶斯理论得到参数θ和β的联合概率密度函数为：$\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 -{\mathrm{\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})]\end{array}---\left(4\right)$由此得到参数更新公式如下：${\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}---\left(5\right)$${\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)}}$采用以上公式对参数μ_θ,k,μ_β,k和进行更新；7)为待评估的模型参数，计算似然函数：$\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 -{\mathrm{\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}---\left(6\right)$似然函数的期望值为：$\begin{array}{c}l\left({\Theta }_k\right|{\hat{\Theta }}_k^{\left(i\right)})=E_{\theta ,\beta |S_{1:k},{\hat{\Theta }}_k^{\left(i\right)}}\left\{\ln \right[p(S_{1:k},\theta ,\beta |{\Theta }_k)\left]\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{1}{2}\ln \left({\sigma }_{0,k}^2\right)-\frac{1}{2}\ln \left({\sigma }_{1,k}^2\right)\\ -\frac{s_1^2-2s_1({\mu }_{\theta ,k}+{\mu }_{B,k}{t}_1)+{\mu }_{\theta ,k}^2+{\sigma }_{\theta ,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)}{2{\sigma }_k^2{t}_1}\\ -\underset{j=2}{\overset{k}{\Sigma }}\frac{{(s_j-s_{j-1})}^2-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)}{2{\sigma }_k^2\mathrm{\Delta t}}\\ -\frac{{\mu }_{\theta ,k}^2+{\sigma }_{\theta ,k}^2-2{\mu }_{\theta ,k}{\mu }_{0,k}+{\mu }_{0,k}^2}{2{\sigma }_{0,k}^2}-\frac{{\mu }_{\beta ,k}^2+{\sigma }_{\beta ,k}^2-2{\mu }_{1,k}{\mu }_{1,k}+{\mu }_{1,k}^2}{2{\sigma }_{1,k}^2}\end{array}---\left(7\right)$其中为第i次评估的结果，计算满足条件$\partial \left[l\left({\Theta }_k\right|{\hat{\Theta }}_k^{\left(i\right)})\right]/\partial {\Theta }_k=0$的参数为第i+1次评估的结果，$\begin{array}{c}{\hat{\sigma }}_k^{2(i+1)}=\frac{1}{k}(\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}\\ +\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}}\end{array}---\left(8\right)$${\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;$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}({\sigma }_{\beta ,k}^2\mathrm{\Delta t}+{\sigma }_k^2)}}\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)$采用下式计算有效粒子数，${\hat{N}}_{\mathrm{eff}}=\frac{1}{\underset{i=1}{\overset{N_s}{\Sigma }}{\left(w_k^i\right)}^2},---\left(12\right)$如果有效粒子数小于阈值NT，需要根据粒子权值大小进行重采样，得到新的粒子集使其满足$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时刻的剩余寿命，γ为滚动轴承失效阈值。