一种改进的设备剩余寿命估计方法

一种改进的设备剩余寿命估计方法

1退化数据驱动的规制方法作为一种新兴技术，规划和健康管理（phm）得到了工程实践的证实，可以降低维护成本，提高设备的可靠性和安全性，减少失败事件的风险。phm的主要问题是状态监测数据，并根据设备的冗余生活时间（ul）进行评估，并确定最合适的维护时间。传统的RUL估计方法以寿命数据为基础,通过对寿命数据的统计分析确定设备的寿命分布,基于此得到RUL的分布.但在实际中,短期内通常难以获取足够多的寿命数据,或者对于一些代价昂贵的设备,获取寿命数据的成本过高.由此,以寿命数据为基础的RUL估计方法难以取得满意的结果.考虑到设备的寿命总是有限的,随着设备的老化或环境、负载的作用,设备的性能会随之退化,当累积到一定程度时,导致退化失效.因此,通过对设备退化数据的监测,建立其退化轨迹的模型,进而估计设备的RUL是一条经济可行的途径.目前,退化数据驱动的方法已发展成为进行设备RUL估计的主要方法,文献对这类方法进行详细系统的综述.有关退化建模的研究,最早可追溯到Gertsbackh和Kordonskiy利用退化数据评定设备的可靠性.Nelson对20世纪90年代以前的退化建模方法进行综述.Lu和Meeker提出一种一般的随机系数回归模型来建模退化测量数据.该文献在退化建模领域具有重要的影响,在该工作基础上出现许多扩展和变形.对于单个服役设备的退化过程来讲,采用该方法就意味着确定的退化轨迹,即所有的模型系数都是确定的.因此,这类方法只适用于描述一批同类设备的退化过程,其涉及的数据分析方法主要是统计分析中的纵向数据分析方法,难以对单个服役设备的退化特征进行具体描述.为了实现对服役设备的退化建模,Gebraeel等提出一类指数退化模型用以描述轴承的磨损退化过程.在这类方法的建模过程中,在获取退化监测数据后,通过Bayesian更新,对模型的随机参数进行后验估计.但是,对于退化模型中非随机的未知参数(如高斯误差项的方差参数、随机参数的先验分布中的参数)没有提出相应的估计方法,而是假设存在多个同类型设备的历史退化数据,通过对这些离线历史数据的统计分析,经验确定模型的非随机参数.实验研究表明,对这些非随机参数选择的恰当与否会较大地影响退化建模及RUL估计的结果,因此有必要研究这些非随机参数的估计方法.再者,这类RUL估计方法仍需要多个同类设备退化的历史数据,但对于具体的设备使用者而言,工程实际中通常难以得到多个设备的数据,因此现有方法具有一定的局限性.考虑到以上的问题,本文提出一种Bayesian更新与期望最大化(ExpectationMaximization,EM)算法协作下退化数据驱动的RUL估计方法,以求实现对单个服役设备退化过程的建模以及RUL估计.具体地,首先利用文献中的指数退化模型来描述设备的退化过程,对模型的随机参数进行Bayesian更新,进而得到RUL的概率分布函数和点估计.其次,对退化模型中的非随机未知参数,利用运行设备到当前时刻的所有监测数据,基于EM算法给出具体的参数估计方法,并证明参数迭代估计中每步得到的结果是唯一的最优解.最后通过数值仿真和实际数据应用研究,表明本文方法可对单个设备退化过程进行建模,有效估计退化模型中的未知参数,进而得到满意的RUL估计.2退出模型与剩余寿命的估计2.1同类型设备中各参数的耦合模型指数类的退化模型是一类典型的描述退化累积过程的模型.到目前为止,这类模型已在设备腐蚀、轴承磨损、液晶显示器等设备的退化建模中得到广泛应用[13,14,15,17,18,19].因此本文主要考虑采用这类模型来表述退化过程的情况.令X(t)表示t时刻的退化量,那么在离散时间监测点t1,t2,…,采用指数退化模型,可将tk时刻的退化量描述为X(tk)=φ+θexpβtk+ε(tk)-σ22),其中φ为已知的常数.考虑到在一批同类型设备中,每个设备在运行过程中的情况都可能有差异,从而显示不同的退化轨迹.因此,在退化模型中,有必要考虑不同设备之间存在的个体差异.本文将θ和β作为随机参数,用以描述设备的个体差异,ε(tk)是随机误差项,这里将θ和β作为随机参数的主要原因在于θ和β对于退化轨迹的期望值有很大影响,而σ2只刻画退化过程的不确定性.为便于后面的建模,在本文中,引入如下文献中广泛采用的模型假设[5,9,12,13,14,15,18]:1)θ和β是高斯分布的,且满足lnθ～N(μ0,σ20)和β～N(μ1,σ21);2)ε(tk)是高斯分布的,即ε(tk)～N(0,σ2),且有ε(0)=0;3)ε(t1),ε(t2),…是独立同分布的随机变量;4)θ,β和ε(tk)是相互独立的.对于指数类的退化模型,通常采用对数变换,得到变换后的形式,以简化处理的程序.因此,本文将tk时刻对数变换后的退化量记为S(tk).由此对应式(1),变换后的退化量可描述为S(tk)=ln(X(tk)-φ)=lnθ-σ22+βtk+ε(tk)=θ′+βtk+ε(tk),其中θ′=lnθ-σ22,由此θ′～N(μ′0,σ20),μ′0=μ0-σ22.基于以上的模型和假设,一旦得到设备退化的监测数据,就可通过Bayesian理论得到θ′和β的后验估计.在后文中,不失一般性,仅考虑φ=0的情况,并且为了符号的简化,记Sk=S(tk)=ln(X(tk)).2.2,令S1:k={S1,S2,…,Sk},其中Sk=ln(X(tk)-φ),那么由于ε(t1),ε(t2),…,ε(tk)是独立同分布的,在给定θ′和β的条件下,采样数据S1:k的联合分布是多变量高斯分布的,即p(S1:kθ′,β)=1k∏j=1√2πσ2exp[-k∑j=1(Sj-θ′-βtj)22σ2].由于θ′和β的先验分布是高斯分布,与采样分布p(S1:kθ′,β)是共轭的,根据Bayesian理论可知:θ′和β关于S1:k的联合后验分布仍然是高斯分布,也就是说θ′,βS1:k～N(μθ′,k,σ2θ′,k,μβ,k,σ2β,k,ρk).因此有p(θ′,βS1:k)∝p(S1:kθ′,β)p(θ′,β)∝12πσθ′,kσβ,k√1-ρ2kexp-12(1-ρ2k)·(θ′-μθ′,k)2σ2θ′,k-2ρk(θ′-μθ′,k)(β-μβ,k)σθ′,kσβ,k+(β-μβ,k)2σ2β,k)].经过一些代数的运算,上式对应的参数μθ′,k,σ2θ′,k,μβ,k,σ2β,k,ρk可由以下命题给出.命题1给定到tk时刻的所有退化观测数据S1:k,θ′和β的联合后验分布为双变量高斯分布,即θ′,βS1:k～N(μθ′,k,σ2θ′,k,μβ,k,σ2β,k,ρk),其中μθ′,k=(k∑i=1Siσ20+μ′0σ2)(k∑i=1t2iσ21+σ2)-(k∑i=1tiσ20)(k∑i=1Sitiσ21+μ1σ2)(kσ20+σ2)(k∑i=1t2iσ21+σ2)-(k∑i=1tiσ21)(k∑i=1tiσ20),μβ,k=(kσ20+σ2)(k∑i=1Sitiσ21+μ1σ2)-(k∑i=1tiσ21)(k∑i=1Siσ20+μ′0σ2)(kσ20+σ2)(k∑i=1t2iσ21+σ2)-(k∑i=1tiσ21)(k∑i=1tiσ20),σ2θ′,k=ˉσ2σ21k∑i=1t2iσ21+σ2(kσ20+σ2)(k∑i=1t2iσ21+σ2)-(k∑i=1ti)2σ20σ21‚σ2β,k=ˉσ2σ20kσ20+σ2(kσ20+σ2)(k∑i=1t2iσ21+σ2)-(k∑i=1ti)2σ20σ21‚ρk=-σ0σ1k∑i=1ti√kσ20+σ2√k∑i=1t2iσ21+σ2,ˉσ2=σ2σ20σ21,μθ′,k、μβ,k是θ′、β的后验估计的均值,σ2θ′,k、σ2β,k为θ′、β后验估计的方差,ρk是对应的相关系数.2.3类目数下t+tkt时培训条件下的规训参数估计基于Bayesian更新的思想,在2.2节得到θ′和β的后验估计后,要实现RUL的估计,首先需要基数据S1:k预测设备未来的退化趋势.基于退化趋势的预测结果,本文采取以下方法进行RUL估计:当退化量超过给定的失效阈值w>0时,认为设备发生失效,进而寿命终止(例如,ISO2372和ISO10816经常被用来确定振动阈值的大小).因此估计RUL的问题就转化为估计退化量达到失效阈值w的时间问题.由此,对于将来t+tk时刻的退化量S(t+tk),在给定S1:k的条件下,有以下命题成立.命题2给定到tk时刻的所有退化观测数据S1:k,则预测的t+tk时刻对应的退化量S(t+tk)为高斯分布,即S(t+tk)S1:k∼Ν(˜μ(t+tk),˜σ2(t+tk)),其中˜μ(t+tk)=μθ′,k+μβ,k(t+tk)-σ22,˜σ2(t+tk)=σ2θ′,k+σ2β,k(t+tk)2+σ2+2ρk(t+tk)σθ′,kσβ,k.(1)令T表示tk时刻的RUL,即从当前时刻tk到失效时间的间隔.由上面给出的退化失效的定义及对退化模型采取的对数变换特征,可得T应满足S(T+tk)=lnw.由此给定S1:k时,RUL的条件累积分布函数FΤ|S1:k(t):(Ζ≥lnw-˜μ(t+tk)√˜σ2(t+tk))=Φ(g(t)),其中g(t)=˜μ(t+tk)-lnw√˜σ2(t+tk)‚Z为标准正态分布的随机变量,ϕ(⋅)为标准正态分布随机变量的累积分布函数.考虑到T表示tk时刻的RUL(到失效的时间间隔),因此T是非负的随机变量.为保证这一事实,将得到的分布FΤ|S1:k(t)在T≥0的条件下进行截尾,即得到如下的估计结果:FΤ|S1:k,Τ≥0(t)=Ρr(Τ≤t|S1:k,Τ≥0)=Ρr(0≤Τ≤t|S1:k)Ρr(Τ≥0|S1:k)=ϕ(g(t))-ϕ(g(0))1-ϕ(g(0)).(2)对应的RUL的条件概率密度函数如下:fΤ|S1:k,Τ≥0(t)=dFΤ|S1:k,Τ≥0(t)dt=ϕ(g(t))1-Φ(g(0))g′(t),(3)其中,ϕ(⋅)为标准正态分布随机变量的概率密度函数.式(2)和式(3)分别得到RUL的条件分布函数和条件概率密度函数,直接通过式(2)和式(3)计算tk时刻RUL的点估计(如期望)较困难,因此在本文中采用如下方法得到RUL的点估计:根据退化失效的定义有S(T+tk)=lnw,于是通过预测的退化量的均值˜μ(t+tk)替代S(T+tk),使得˜μ(t+tk)=lnw,从而得到tk时刻RUL的点估计RULk如下:RULk=lnw-μθ′,k+σ22μβ,k-tk,其中,RULk可作为RUL的期望的近似.在以上的退化建模和RUL估计过程中,参数σ2,μ′0,μ1,σ20,σ21都是未知的.文献中利用同类设备的历史退化数据来估计σ2,μ′0,μ1,σ20,σ21,而且一旦由历史数据确定这些参数后,就不再更新.然而,对于实际服役设备而言,需要通过实时监测的数据不断更新模型参数,使得估计的RUL能更好地反映设备的当前实际运行状态.第3节将主要讨论利用到当前时刻tk的监测数据S1:k,基于EM算法的未知参数估计问题.3为了估计未知参数Θ=[σ2,μ′0,μ1,σ20,σ21],本文采用极大似然估计的方法,在tk时刻的监测数据S1:k得到后,计算关于S1:k的对数似然函数如下:lk(Θ)=ln[p(S1:k|Θ)]=k∑j=2ln[p(Sj|S1:j-1,Θ)],其中,p(S1:kΘ)表示退化数据S1:k的联合概率密度函数.因此在tk时刻Θ的极大似然估计ˆΘk:Θ^k=argmaxΘlk(Θ),其中,argmaxΘlk(Θ)表示lk(Θ)关于Θ取最大值时对应的参数向量Θ的取值,即Θ^k.然而在退化模型式(2)中,参数θ′和β是随机的,且不能直接观测得到,因此采用上式难以直接优化得到极大似然估计Θ^k.EM算法为解决这一极大似然估计问题提供可行的解决思路,其主要思想是基于Bayesian规则,通过p(S1:kΘ)和p(S1:k,θ′,β|Θ)之间的关系,使得估计Θ可通过迭代以下两步实现:E-step和M-step.1)E-step.计算l(ΘΘ^k(i))=Eθ′,β|S1:k,Θ^k(i){lnp(S1:k,θ′,βΘ)},(4)其中,Θ^k(i)表示基于退化数据S1:k进行估计的过程中,第i次迭代得到的结果.2)M-step.计算Θ^k(i+1)=argmaxΘl(Θ|Θ^k(i)).(5)上面的迭代过程将得到一系列参数估计{Θ^(0)k,Θ^(1)k,Θ^(2)k,…}.通过对EM算法的理论分析可得到随着迭代次数的增加,似然函数lk(Θ)也将不断增加,因此在最大化似然函数的意义下,将得到越来越好的估计.在实际应用EM算法进行参数估计的过程中,通常在Θ^(i)k和Θ^(i+1)k的差异小于一个较小的数值(阈值)时,终止迭代计算的过程,将最后一次估计的结果作为tk时刻的参数估计结果.有关EM算法的理论性质和收敛性分析可参考文献、,限于篇幅,这里不再赘述.下面给出利用上述EM算法对未知参数Θ进行估计的具体过程.为了表示估计的参数依赖于到当前时刻tk的所有退化监测数据S1:k,此后,将基于S1:k估计的参数表示为Θk=[σk2,μ′0,k,μ1,k,σ0,k2,σ1,k2].令EM算法中第i次迭代得到的估计为Θ^ki=[σ^k2(i),μ^′0,k(i),μ^1,k(i),μ^i,k2(i)],那么完整的对数似然函数可表示为lnp(S1:k,θ′,βΘk)=lnp(S1:kθ′,β,Θk)+lnp(θ′,βΘk)=-k+22ln2π-k2lnσk2-∑kj=1(Sj-θ′-βtj)22σk2-12lnσ0,k2-12lnσ1,k2-(θ′-μ′0,k)22σ0,k2-(β-μ1,k)22σ1,k2.(6)因此,从推论、式(4)和式(6)可得l(ΘkΘ^k(i))=Eθ′,β|S1:k,Θ^k(i){lnp(S1:k,θ′,βΘk)}=-k+22ln2π-k2lnσk2-∑kj=1Sj2-2Sj(μθ′,k+μβ,ktj)+μθ′,k2+σθ′,k2+2tj(ρkσθ′,kσβ,k+μθ′,kμβ,k)+tj2(μβ,k2+σβ,k2)2σk2-12lnσ0,k2-12lnσ1,k2-μθ′,k2+σθ′,k2-2μθ′,kμ′0,k+μ′0,k22σ0,k2-μβ,k2+σβ,k2-2μβ,kμ1,k+μ1,k22σ1,k2.(7)令∂l(Θk|Θ^k(i))∂Θk=0,可得第i+1步的参数估计Θ^k(i+1)如下:σ^k2(i+1)=1k∑kj=1(Sj2-2Sj(μθ′,k+μβ,ktj)+μθ′,k2+σθ′,k2+2tj(ρkσθ′,kσβ,k+μθ′,kμβ,k)+tj2(σβ,k2+μβ,k2)),(8)μ^′0,k(i+1)=μθ′,k,σ^0,k2(x+1)=σθ′,k2,μ^1,k(i+1)=μβ,k,σ^1,k2(i+1)=σβ,k2.(9)命题3由式(8)和式(9)得到的估计Θ^k(i+1)是l(ΘkΘ^k(i))的唯一最大值点.证明由式(7)可知,从式(8)和式(9)得到的Θ^k(i+1)是满足∂l(Θk|Θ^k(i))∂Θk=0的唯一解.因此,为证明Θ^k(i+1)是l(ΘkΘ^k(i))的唯一最大值点,计算∂2l(Θk|Θ^k(i))∂Θk∂ΘkΤ如下:∂2l(Θk|Θ^k(i))∂Θk∂ΘkΤ=[k2σk4-φσk600000-1σ0,k20μ0,k´-μθ´,kσ0,k4000-1σ1,k20μ1,k-μβ,kσ1,k40μ0,k´-μθ´,kσ0,k4012σ0,k4-ψ1σ0,k6000μ1,k-μβ,kσ1,k4012σ1,k4-ψ2σ1,k6],(10)其中ψ1=μθ′,k2+σθ′,k2-2μθ′,kμ0,k+μ0,k2,ψ2=μβ,k2+σβ,k2-2μβ,kμ1,k+μ1,k2,φ=∑kj=1(Sj2-2Sj(μθ′,k+μβ,ktj)+μθ′,k2+σθ′,k2+2tj(ρkσθ′,kσβ,k+μθ′,kμβ,k)+tj2(σβ,k2+μβ,k2)).可证明在由式(8)和式(9)得到的估计Θk=Θ^k(i+1)时,式(10)中矩阵对应的顺序主子式满足Δ1Θk=Θ^k(i+1)＜0,Δ2Θk=Θ^k(i+1)＞0,Δ3Θk=Θ^k(i+1)＜0,Δ4Θk=Θ^k(i+1)＞0,Δ5Θk=Θ^k(i+1)＜0.由此表明,矩阵式(10)在Θk=Θ^k(i+1)时是负定的.再由Θ^k(i+1)是满足∂l(Θk|Θ^k(i))∂Θk=0的唯一解,可知式(8)和式(9)得到的估计Θ^k(i+1)是l(Θk|Θ^k(i))的唯一最大值点.上述的参数估计方法是基于到tk时刻的所有退化数据S1:k估计对参数的结果,由于S1:k对监测时刻的依赖关系,在任一监测时刻的数据得到后都可采用上面方法进行参数估计.此外,通过命题3可看出,在本文提出的基于EM算法的参数估计方法中,M步只需要单次的计算,而且每次计算的结果都为式(5)的最优解,由此使得参数估计方法具有快速且非常简单的计算程序,这对于工程实现具有重要的现实意义.4数值模拟和示例研究4.1不同方法在不同退化量时的仿真结果为了验证本文提出的退化建模和参数估计方法,以及对应的剩余寿命估计方法,在本小节将通过一个数值例子实现本文方法,并与文献方法对比.为产生仿真的退化数据,本文利用文献中建立的退化模型来产生一组用于验证本文方法对单个服役设备退化过程建模和剩余寿命估计的模拟数据,即使用X(t)=φ+θexp(βt+ε(t)-σ22).值得注意的是,在文献中需要仿真多组数据以实现模型参数的估计.相比之下,本文方法更适合于单个服役设备.在数据产生的过程中,令模型的参数满足以下条件:φ=0,lnθ-σ2/2～N(0.016,2×10-8),β～N(0.005,10-8),σ2=0.008.但给定采样间隔为5,失效阈值w=12.4时,得到一仿真的退化轨迹如图1所示,共得到99个监测数据,对应的超过失效阈值的时间约为497.5.为了与文献方法对比,令随机参数θ和β的先验分布满足lnθ～N(0.05,0.02)和β～N(0.1,0.01),且令σ2的初始值为0.4.分别采用本文方法和文献的方法,采用式(6)中第一个方程得到的退化轨迹的估计,如图1所示.本文方法和文献的方法估计退化量时得到的结果与仿真数据之间的均方误差分别为0.0425和0.6457.明显可见,本文方法由于引入EM算法对随机参数先验分布中的参数和σ2进行估计,可明显提高退化建模的精度.相比之下,由于文献中没有考虑对θ和β的先验分布中的参数及σ2进行实时估计,因此退化估计的结果很多程度上依赖于初始参数的选择,而本文方法可有效克服这一问题.下面以先验参数μ′0=μ0-σ2/2和μ1为例,采用本文方法,随着退化数据的累计,对μ′0=μ0-σ2/2和μ1估计的结果如图2所示.从图2估计的结果可看出,对μ′0和μ1最终的估计结果将随数据的积累不断逼近仿真数据产生时所用参数的真实值,其它参数的估计结果类似,限于篇幅不再列出.以上结果表明本文方法的有效性.对于本文方法和文献方法在剩余寿命的点估计,可通过式(6)中的第2个方程计算得到.对应于各个采样点的估计结果对比如图3所示.从图中可看出,本文方法较文献方法在剩余寿命的点估计的性能方面,有明显优势,可快速精确地估计剩余寿命.具体地,本文方法和文献方法得到的剩余寿命的点估计与仿真得到的实际剩余寿命之间的平均相对误差分别为0.1014和0.5705,表明本文方法可明显提高剩余寿命估计的精度.为了对比剩余寿命的分布,以最后一个采样点得到的结果为例,进一步进行对比.对应于本文方法得到的剩余寿命的条件概率密度函数与文献得到的密度函数如图4.图4进一步表明本文方法在估计剩余寿命时也能产生更好结果,文献方法估计得到剩余寿命明显与实际情况存在很多差异.另外,采用本文方法得到的剩余寿命分布更窄,表明剩余寿命估计的不确定性更小.这一点在PHM中的管理决策中尤其重要,具体的讨论可参考文献.通过上述这些仿真对比,充分说明本文方法的优越性.4.2剩余寿命估计模型本小节基于惯性导航平台的退化监测数据,将提出的剩余寿命估计方法应用到惯性导航平台的剩余寿命预测中.惯性导航平台在军事设备、航空航天系统等领域有着重要应用,其健康状态的主要评价标准是平台系统的精度测试,然而平台系统的精度主要受陀螺仪漂移系数的影响.在实际中,平台系统陀螺转子高速旋转势必造成转轴的磨损,随着工作时间的累积,引起漂移系数的增大和系统性能的退化,最终导致失效.其次,随着导航系统工作环境的变化(如温度和湿度),陀螺仪的性能也会逐步发生退化,进而影响这个导航系统的精度.在这种情况下,及时有效估计导航平台的