Degradation Modeling and Remaining Useful Life Prediction with Bivariate Time Scale
-
摘要: 基于退化建模的剩余寿命预测（Remaining useful life，RUL）是当前可靠性领域研究的热点.现有的退化模型都是针对单个时间尺度下的退化设备，缺少对设备性能变化与多个时间尺度相关的退化建模与剩余寿命预测方法.鉴于此，本文基于Wiener过程提出了一种双时间尺度随机退化建模与剩余寿命预测方法，用随机比例系数描述不同时间尺度之间的不确定关系，推导出丫首达时间意义下设备的双时间尺度剩余寿命分布，讨论了其与基于单时间尺度退化模型得到的剩余寿命分布之间的关系，并给出了基于历史退化数据的未知参数极大似然估计方法.最后，将所提方法应用到惯性平台关键器件陀螺仪的退化建模与剩余寿命预测中，验证了方法的有效性.Abstract: Degradation modeling based remaining useful life (RUL) prediction is currently a research hotspot in the field of reliability. However, almost all the existing models consider the degrading equipment with one time scale, and no particular approach has been developed for deteriorating equipment whose degradation process involves more than one time scale. To fill this gap, a bivariate-time-scale degradation model, in which the uncertain relationship between the two time scales is depicted through a stochastic proportional coefficient, is presented based on Wiener process. Under the concept of the first hitting time, a bivariate-time-scale RUL distribution is derived and its relationship with the RUL distribution under single time scale is discussed. Besides, a maximum likelihood estimation method is introduced to estimate unknown parameters using historical degradation data. A case study of degradation data from gyroscopes in an inertial platform demonstrates the effectiveness of the proposed method.
-
受工作环境和自身材料性能恶化等因素综合作用的影响, 设备的性能不可避免地发生退化, 并逐步累积成失效.对于航空航天、高铁列车、运载火箭、钻井平台等大型复杂系统, 一些关键设备的退化失效将带来难以估量的经济损失、人员伤亡和环境破坏[1].如果能够在关键设备失效前, 准确地对其进行剩余寿命预测(Remaining useful life, RUL), 就能够及时主动地采取预防措施, 避免灾难的发生.传统的基于失效数据的方法需要同类设备的失效数据, 耗时长、成本高, 不能适应新形势下设备寿命预测的需求.基于设备退化测量数据建立退化模型, 进而对失效设备RUL预测的方法是一条经济可行的途径[2-4].因此, 设备的退化建模与剩余预测逐渐得到学者和工程技术人员的关注、研究和应用, 成为可靠性领域的一个研究热点[5-9].
设备的退化轨迹往往具有随机性, 故基于随机过程的方法成为退化建模中的主流方法.文献中已有的相关研究可分为两类:针对硬失效的退化建模方法和针对软失效的方法.针对硬失效的建模方法, 通常将设备的随机退化过程作为协变量融入到比例风险模型的失效率函数中, 刻画退化过程对设备失效的影响, 进而实现设备的RUL预测[10-11]; 针对软失效的方法, 一般将设备随机退化过程首次到达某个失效阈值的时间定义为设备的寿命, 通过求解该首达时间的分布实现设备的寿命和剩余寿命预测[4-5].由于第二类方法具有物理意义直观明确、计算简单便捷、便于融入随机效用等优点, 且工程实际中大量随机退化系统的失效模式属于软失效, 因此本文基于第二类方法展开研究.
时间尺度是退化建模和剩余寿命预测的基本要素之一, 也是一个常被现有方法所忽视的问题.事实上, 目前文献中绝大多数的退化模型都只考虑了单一的时间尺度, 包括退化模型中的自然时间、Li-ion电池退化模型中的充放电次数、金属试样疲劳试验中的循环次数等, 使用的都是单个时间尺度.在对实际系统的退化建模和RUL预测中, 也会面临双时间尺度、甚至多时间尺度的问题.例如, 在评估汽车轮胎保修期之前, 需要确定一个时间和行驶路程双时间尺度下的寿命分布[12]; 在对飞机起落架进行检修决策时, 需要考虑起落架服役的总时长、处于工作状态的时长以及经受的起落次数三个不同时间尺度下的寿命[13-14].再如, 研究表明通电测试对贮存中的惯性仪表的性能退化有明显的加速作用, 而维护和备件订购等决策都基于自然时间[15].因此, 有必要考虑双时间尺度下惯性仪表的退化建模和剩余寿命预测问题.此外, 不同时间尺度通常也不是相互独立的, 由于设备运行过程中的随机性, 不同时间尺度之间的关系也难以通过确定性的函数描述.这些问题虽有一些学者展开了研究, 但现有的双时间尺度或多时间尺度失效模型都是基于传统的失效数据[13-14, 16], 尚未发现关于双时间尺度下设备随机退化过程建模和RUL预测的文献报道.
鉴于以上讨论, 本文提出了一种基于Wiener过程的双时间尺度随机退化建模和剩余寿命预测方法, 利用随机比例系数描述两个时间尺度之间的关系, 实现双时间尺度下设备的退化过程建模.在此基础上, 推导了双时间尺度下设备寿命和RUL分布的解析解, 并讨论了其与基于单时间尺度退化模型得到结果之间的关系.同时, 给出了基于历史退化数据的未知参数极大似然估计方法.最后, 通过惯性平台关键器件陀螺仪的退化建模和RUL估计验证了所提方法的有效性.
1. 随机退化建模
基于Wiener过程的退化模型能够对非单调的退化设备进行建模, 且具有解析形式的首达时间分布, 已经广泛应用于激光发生器、LED灯、金属材料等的退化建模和剩余寿命预测[17-19].因此, 本文基于Wiener过程进行双时间尺度下设备的退化建模.下面, 首先给出单时间尺度下基于Wiener过程的退化模型.
1.1 单时间尺度随机退化建模
令$t$时刻设备的退化状态为$X(t)$, 基于Wiener过程来刻画设备的退化, 即
$ X(t) = {\lambda _0} + {\lambda _1}t + {\sigma _B}B(t) $
(1) 其中, $\lambda_0$表示设备的初始退化水平, 漂移系数$\lambda_1$表示设备的退化速率, 扩散系数$\sigma_B$和标准布朗运动$B(t)$一起, 描述系统健康状态的动态不确定性. $\lambda_0$, $\lambda_1$, $\sigma_B$及$B(t)$相互独立.这里假设设备的退化表现为指标$X(t)$的增加, 对于性能指标逐渐降低的设备, 其退化过程可通过简单的数学处理使其转化为逐渐增加的退化过程.
显然, 式(1) 所描述的退化模型适合对线性的退化过程进行建模.当设备的退化轨迹呈现非线性特性时, 可通过适当的变换, 将退化轨迹线性化之后进行建模.这种变换可以通过对退化水平本身实现, 例如将指数型的退化轨迹对数化[20], 也可以通过时间尺度变换实现[21].下文中的方法可以扩展到可线性化的非线性退化建模中与RUL预测中.
通常当设备的状态超过某个预定的失效阈值$w$时, 不再满足实际使用需求, 认为设备失效.故将性能指标超过阈值$w$的首达时间定义为设备的寿命$T$, 即
$ T = \inf \{ t:X(t) \geq w|t \geq 0\} $
(2) 相应地, $t_k$时刻系统的剩余寿命$L_k$可表示为
$ {L_k} = \inf \{ {l_k}:X({t_k} + {l_k}) \geq w|X({t_k}) \leq w\} $
(3) 1.2 双时间尺度随机退化建模
当考虑两个时间尺度$t$和$\tau$时, 设备的性能退化状态$X(t, \tau)$可以通过式(4) 描述
$ X(t, \tau ) = {\lambda _0} + {\lambda _1}t + {\sigma _B}B(t) + {\lambda _2}\tau + {\sigma _W}W(\tau ) $
(4) 其中, $W(\tau)$是时间尺度$\tau$下的一个标准布朗运动, 与扩散系数$\sigma_W$一起表征由时间尺度$\tau$引入的动态不确定性, $\lambda_2$表示与时间尺度$\tau$相关的退化速率.当然, 也可以借鉴双时间尺度失效模型的常用的方法, 解算出设备的使用率, 进而借鉴加速退化模型的方法把使用率结合到退化模型中, 把问题简化成一个维度, 并考虑两个时间尺度之间的关系[22].为突出双时间尺度与单时间尺度退化模型的联系与区别, 本文将设备寿命定义为
$ [T, \Lambda]\! = \inf \{ t:X(t, \tau )\! \geq w, \Lambda (t) = \tau |t \geq 0, \tau \! \geq 0\} $
(5) 其中, $\Lambda(t)$用于描述时间尺度$t$与$\tau$之间的关系.
模型(1) 中假设漂移系数$\lambda_1$与扩散系数$\sigma_B$相互独立, 模型(4) 中同样假设漂移系数$\lambda_1$, $\lambda_1$与扩散系数$\sigma_W$, $\sigma_B$相互独立.这种假设被大多数基于Wiener过程的退化模型所采用.然而, Ye等最近考虑了漂移系数与扩散系数之间的关系, 假设其之间满足比例关系[23].当考虑漂移系数与扩散系数之间的关系时, 同样可基于文献[23]中的寿命分布结果, 构建双时间尺度下的退化模型, 推导过程类似.
模型(4) 与现有的多退化量模型差异显著.多退化量指的是表征设备性能的指标是多元的, 即单维的退化状态$X(t)$应该表示成多维向量的形式$\boldsymbol{X}(t)$.而本文考虑的双时间尺度退化模型中, 刻画设备性能退化状态的随机变量仍是一维的, 但是时间尺度由$t$变成了$[t, \tau]$.可以认为, 多退化量模型针对的是函数向量表示的退化过程, 而本文考虑的是用两个时间尺度的二元标量函数刻画的退化过程.
如前所述, 时间尺度之间通常是相互关联的.本文也主要考虑两个时间尺度相互关联时的情形.当两个时间尺度$t$和$\tau$相互不相关时, 寿命分布的推导及参数估计等都将更加简单.为简单起见, 两个时间尺度$t$和$\tau$之间的关系可以用如下函数表述
$ \tau = \Lambda (t) = \gamma t $
(6) 其中, $\gamma$表示时间尺度$\tau$所占时间尺度$t$的比例.例如, 若$t$和$\tau$分别表示自然时间与工作时间, 则$\gamma$为工作率.为了描述设备运行过程中的不确定性, 这里假设$\gamma$为服从正态分布的随机变量.当然, 根据设备的具体运行情况, 也可以选择服从其他类型分布的随机变量刻画不同时间尺度之间的关系.此外, 这种关联并不要求其中一个时间尺度为另外一个时间尺度的一部分, 例如实例分析中通电时间与储存时间的关系.对于其他不同形式的双时间尺度, 例如轮胎保养中涉及的行驶里程与使用时间[23], 虽互不包含, 但仍可用式(6) 进行描述.
类似地, 设备在某给定时刻$t_k$的剩余寿命$[L_k, \Psi_k]$可定义为
$ \begin{align}\label{EQ7} [{L_k}, {\Psi _k}] \approx&\ \inf \{ {l_k}:X({t_k} + {l_k}, {\tau _k} + {\eta _k}) \geq w, \nonumber \\ &\ \Lambda ({t_k} + {l_k}) - \Lambda ({t_k}) = {\eta _k}|X({t_k}, {\tau _k}) < w \} \end{align} $
(7) 本文的目的是求解$[T, \Lambda]$和$[L_k, \Psi_k]$的概率密度函数(Probability density function, PDF)及累积概率分布函数(Cumulative distribution function, CDF).为叙述方便, 下文分别将单时间尺度和双时间尺度退化模型记为$M_1$和$M_2$.
对比单时间尺度退化模型(1) 和双时间尺度退化模型(4) 可知, 本文提出的双时间尺度模型是基于单时间尺度提出来的, 能够将单时间尺度模型包含为特例, 当$\lambda_2=0$且$\sigma^2_W=0$时, 双时间尺度模型即为单时间尺度模型; 同时, 推导双时间尺度下的寿命分布也需要基于单时间尺度退化模型的相关结果.然而, 双时间尺度退化模型并不仅仅是两个单时间尺度模型的简单叠加, 还需要考虑了两个时间尺度之间的相互关系.此外, 单时间尺度退化模型结构简单, 需要估计的参数少, 但往往忽略了另一个时间尺度对退化过程的影响; 双时间尺度模型虽然需要确定的参数更多, 但能够提供更多的寿命分布信息, 同时也为考虑不同因素对退化过程的影响提供了一条新的途径.
2. 寿命分布和剩余寿命分布推导
为了计算$M_1$下设备寿命的PDF, 需要用到单时间尺度下的一些寿命分布的一些结论以及条件概率公式.具体地
$ \begin{align}\label{EQ8} & {f_{T, \Lambda |\boldsymbol{\theta}, {M_2}}}(t, \tau |\boldsymbol{\theta}, {M_2}) =\nonumber \\ & \qquad {f_{T|\Lambda, \boldsymbol{\theta}, {M_1}}}(t|\tau , \boldsymbol{\theta}, {M_1}) {f_{\Lambda |\boldsymbol{\theta} , {M_1}}}(\tau |\boldsymbol{\theta}, {M_1}) \end{align} $
(8) 其中, ${\boldsymbol{\theta}^{\rm{T}}} = [{\lambda _0}, {\lambda _1}, {\lambda _2}]$, ${f_{T|\Lambda, \boldsymbol{\theta}, {M_1}}}(t|\tau, \boldsymbol{\theta}, {M_1})$表示$t$时刻$\tau = \Lambda (t)$条件下的$T$的PDF, ${f_{\Lambda |\boldsymbol{\theta}, {M_1}}}(\tau |\boldsymbol{\theta} $, ${M_1})$用于描述时间尺度$t$与$\tau$之间的关系.为计算${f_{T|\Lambda, \boldsymbol{\theta}, {M_1}}}(t|\tau, \boldsymbol{\theta}, {M_1})$, 可将$X(t, \tau)$超过首达固定阈值$w$的时间分布, 转化成$X(t)$首达随机失效阈值$\tilde{w}$的时间分布.这里的随机阈值$\tilde{w}$可以表示为$\tilde{w} =w- {\lambda _2}\tau - {\sigma _W}W(\tau )$, 根据Wiener过程的性质, 随机阈值$\tilde{w}$服从正态分布, 即$ \tilde{w}\sim{\rm N}(w -$ ${\lambda _2}\tau, \sigma _W^2\tau )$.设备在$[{t_k}, {\tau _k}]$时刻RUL的分布${f_{{L_k}, {\Psi _k}|\boldsymbol{\theta}, {M_2}}}( {l_k}, {\eta _k} |\boldsymbol{\theta}, {M_2})$可以采用类似的方法进行计算, 即
$ \begin{align}\label{EQ9} &{f_{{L_k}, {\Psi _k}|\boldsymbol{\theta}, {M_2}}}({l_k}, {\eta _k}|\boldsymbol{\theta}, {M_2})=\nonumber\\ &\qquad {f_{{L_k}|{\Psi _k}, \boldsymbol{\theta}, {M_2}}}({l_k}|{\eta _k}, \boldsymbol{\theta}, {M_1}) {f_{{\Psi _k}|\boldsymbol{\theta} , {M_2}}}({\eta _k}|\boldsymbol{\theta}, {M_1}) \end{align} $
(9) 因此, 首先给出单时间尺度$t$, 即$M_1$下设备的寿命分布.根据文献[18]中的结论, 对于式(1) 描述的随机退化设备, 给定$\lambda_0$和$\lambda_1$时, 其寿命$T$和RUL$L_k$都有解析形式的PDF, 即
$ \begin{align}\label{EQ10} &{f_{T|\boldsymbol{\lambda}, {M_1}}} (t|\boldsymbol{\lambda}, {M_1}) = \frac{{{w_0}}}{{\sqrt {2\pi \sigma _B^2{t^3}} }} \exp \left[{-\frac{{{{({w_0}-{\lambda _1}t)}^2}}}{{2t\sigma _B^2}}} \right]\end{align} $
(10) $ \begin{align}&{f_{{L_k}|\boldsymbol{\lambda}, {M_1}}} ({l_k}|\boldsymbol{\lambda}, {M_1}) =\\ &\qquad \qquad \frac{{{w_k}}}{{\sqrt {2\pi \sigma _B^2{l_k^3}} }} \exp \left[{-\frac{{{{({w_k}-{\lambda _1}t)}^2}}}{{2l_k\sigma _B^2}}} \right] \end{align} $
(11) 其中, ${\boldsymbol{\lambda}^{\rm{T}}} = [{\lambda _0}, {\lambda _1}]$, ${w_0} = w - {\lambda _0}$, ${w_k} = w - {\lambda _k}$.
相应地, 如式(12) 和式(13) 所示, 设备的寿命$T$和RUL$L_k$的CDF可据逆高斯分布的性质得到, 已被许多文献引用与验证, 证明过程可参考文献[24].
$ \begin{align}\label{EQ12} {F_{T|\boldsymbol{\lambda}, {M_1}}}&(\left. t \right|\boldsymbol{\lambda}, {M_1}) = 1 - \Phi \left( {\frac{{{w_0} - {\lambda _1}t}}{{{\sigma _B}\sqrt t }}} \right) +\nonumber\\ & \exp \left[{\frac{{2{\lambda _1}{w_0}}}{{\sigma _B^2}}} \right]\Phi \left( {\frac{{ - {w_0} - {\lambda _1}t}}{{{\sigma _B}\sqrt t }}} \right) \end{align} $
(12) $ \begin{align}\label{EQ13} {F_{\left. {{L_k}} \right|\boldsymbol{\lambda}, {M_1}}}&(\left. {{l_k}} \right|\boldsymbol{\lambda}, {M_1}) = 1 - \Phi \left( {\frac{{{w_k} - {\lambda _1}{l_k}}}{{{\sigma _B}\sqrt {{l_k}} }}} \right) + \nonumber\\ &\exp \left[{\frac{{2{\lambda _1}{w_k}}}{{\sigma _B^2}}} \right]\Phi \left( {\frac{{ - {w_k} - {\lambda _1}{l_k}}}{{{\sigma _B}\sqrt {{l_k}} }}} \right) \end{align} $
(13) 根据式(10) 和式(11), 利用文献[18]中的引理1, 可以推导出单时间尺度下设备寿命和RUL的条件PDF, 即${f_{T|\Lambda, \boldsymbol{\theta}, {M_1}}}(t|\tau, \boldsymbol{\theta}, {M_1})$和${f_{{L_k}, {\Psi _k}|\boldsymbol{\theta}, {M_1}}}(\left. {{l_k}, {\eta _k}} \right|\boldsymbol{\theta}, {M_1})$的表达式, 分别如式(14) 和式(15) 所示.
$ \begin{align}\label{EQ14} &{f_{T|\Lambda, \boldsymbol{\theta}, {M_2}}}(\left. t \right|\tau, \boldsymbol{\theta}, {M_2}) =\notag\\[2mm] &\qquad \sqrt {\frac{{\sigma _B^2}}{{2\pi {t^2}(\sigma _B^2t + \sigma _W^2\tau )}}} \times\notag\\[2mm] &\qquad \frac{{{w_0}\sigma _B^2t + (\sigma _W^2\tau - \sigma _B^2){\lambda _2}t \tau }}{{\sigma _B^2t + \sigma _W^2\tau }} \times\nonumber\\[2mm] & \qquad \exp \left[{-\frac{{{{({w_0}-{\lambda _1}t-{\lambda _2}\tau )}^2}}}{{2(t\sigma _B^2 + \sigma _W^2\tau )}}} \right] \end{align} $
(14) $ \begin{align} &{f_{{L_k}|{\Psi _k}, \boldsymbol{\theta}, {M_2}}}(\left. {{l_k}} \right|{\eta _k}, \boldsymbol{\theta}, {M_2}) =\notag\\[2mm] &\qquad \sqrt {\frac{{\sigma _B^2}}{{2\pi l_k^2(\sigma _B^2{l_k} + \sigma _W^2{\eta _k})}}} \times\nonumber\\[2mm] & \qquad \frac{{{w_k}\sigma _B^2{l_k} + (\sigma _W^2{\eta _k} - \sigma _B^2){\lambda _2}{l_k}{\eta _k}}}{{\sigma _B^2{l_k} + \sigma _W^2{\eta _k}}} \times\nonumber \\[2mm] &\qquad\exp \left[{-\frac{{{{({w_k}-{\lambda _1}{l_k}- {\lambda _2}{\eta _k})}^2}}}{{2(\sigma _B^2{l_k} + \sigma _W^2{\eta _k})}}} \right] \end{align} $
(15) 利用文献[19]中的相关结论, 其对应的CDF分别如式(16) 和式(17) 所示.
$ \begin{align}\label{EQ16} &{F_{T|\Lambda, \boldsymbol{\theta}, {M_2}}}(\left. t \right|\tau, \boldsymbol{\theta}, {M_2}) = \nonumber \\ &\qquad 1 - \Phi \left( {\frac{{{w_0} - {\lambda _1}t - {\lambda _2} \tau }}{{\sqrt {\sigma _B^2t + \sigma _W^2\tau } }}} \right) +\nonumber \\ &\qquad \exp \left[{\frac{{2{\lambda _1}({w_0}-{\lambda _2}\tau )}}{{\sigma _B^2}} + \frac{{2\lambda _1^2\sigma _W^2\tau }}{{\sigma _B^4}}} \right] \times\nonumber \\ &\qquad \Phi \left( {\frac{{ - {w_0} - {\lambda _1}t - {\lambda _2}\tau - \dfrac{2{\lambda _1}\sigma _W^2\tau }{\sigma _B^2}}}{{\sqrt {\sigma _B^2t + \sigma _W^2\tau } }}} \right) \end{align} $
(16) $ \begin{align} &{F_{{L_k}|{\Psi _k}, \boldsymbol{\theta}, {M_2}}}(\left. {{l_k}} \right|{\eta _k}, \boldsymbol{\theta}, {M_2}) = \nonumber \\ &\qquad 1 - \Phi \left( {\frac{{{w_k} - {\lambda _1}{l_k} - {\lambda _2}{\eta _k}}} {{\sqrt {\sigma _B^2{l_k} + \sigma _W^2\tau } }}} \right) +\nonumber \\ &\qquad \exp \left[{\frac{{2{\lambda _1}({w_k}-{\lambda _2}{\eta _k})}}{{\sigma _B^2}} + \dfrac{{2\lambda _1^2\sigma _W^2{\eta _k}}}{{\sigma _B^4}}} \right] \times\nonumber \\[1mm] &\qquad \Phi \left( {\frac{{ - {w_0} - {\lambda _1}{l_k} - {\lambda _2}{\eta _k} - \dfrac{2{\lambda _1}\sigma _W^2{\eta _k}}{\sigma _B^2}}}{{\sqrt {\sigma _B^2{l_k} + \sigma _W^2{\eta _k}} }}} \right) \end{align} $
(17) 根据式(6) 中的相关假设, ${f_{\Lambda |\boldsymbol{\theta}, {M_1}}}(\tau |\boldsymbol{\theta}, {M_1})$和${f_{{\Psi _k}|\boldsymbol{\theta}, {M_1}}}({\eta _k}|\boldsymbol{\theta}, {M_1})$为正态分布, 表达式为
$ \begin{align}\label{EQ18} &{f_{\Lambda |\boldsymbol{\theta}, {M_1}}}(\tau |\boldsymbol{\theta}, {M_1}) = \frac{1}{{\sqrt {2\pi \sigma _\gamma ^2{t^2}} }} \times\nonumber \\ &\qquad \exp \left[{-\frac{{{{(\tau-{\mu _\gamma }t)}^2}}}{{2\sigma _\gamma ^2{t^2}}}} \right] \\ \end{align} $
(18) $ \begin{align}\label{EQ19}& {f_{{\Psi _k}|\boldsymbol{\theta}, {M_1}}}({\eta _k}|\boldsymbol{\theta}, {M_1}) = \frac{1}{{\sqrt {2\pi \sigma _\gamma ^2\eta _k^2} }} \times\nonumber \\&\qquad \exp \left[{-\frac{{{{({\eta _k}-{\mu _\gamma }{l_k})}^2}}}{{2\sigma _\gamma ^2\eta _k^2}}} \right] \end{align} $
(19) 将式(14)、(18) 及式(15)、(19) 代入式(8) 和式(9), 可得到双时间尺度下设备寿命与RUL的分布.
进一步, 可以基于以上双时间尺度寿命分布的结果得到一些有用的结论.具体地, 通过求解式(8) 和式(9) 关于时间尺度$t$和$\tau$的边缘分布, 可得到单个时间尺度下的寿命和RUL分布, 如式(20) 和式(21) 所示.若$t$为设备的储存寿命, $\tau$为设备的累计工作时间, 可通过计算$T$和$L_k$的边缘分布, 求解考虑测试影响的系统的贮存寿命和剩余贮存寿命分布.
$ \begin{align}\label{EQ20} &{f_{T|\boldsymbol{\theta}, {M_2}}}(\left. {t, \tau } \right|\boldsymbol{\theta}, {M_2})=\\ &\qquad \int {{f_{T, \Lambda |\boldsymbol{\theta}, {M_1}}}(t, \tau |\boldsymbol{\theta}, {M_1})} {\rm d}\tau = \nonumber\\ &\qquad \int {f_{T|\Lambda, \boldsymbol{\theta}, {M_1}}}(t|\tau , \boldsymbol{\theta}, {M_1}) \times \nonumber\\ &\qquad {f_{\Lambda |\boldsymbol{\theta}, {M_1}}}(\tau |\boldsymbol{\theta}, {M_1}){\rm d}\tau\end{align} $
(20) $ \begin{align} & {f_{{L_k}|\boldsymbol{\theta}, {M_2}}}(\left. {{l_k}} \right|v, {M_2})= \nonumber \\ &\qquad \int {{f_{{L_k}, {\Psi _k}|\boldsymbol{\theta}, {M_1}}}(\left. {{l_k}, {\eta _k}} \right|\theta, {M_1})} {\rm d}{\eta _k} =\nonumber\\ &\qquad \int {f_{{L_k}|{\Psi _k}, \boldsymbol{\theta} , {M_2}}}({l_k}|{\eta _k}, \boldsymbol{\theta}, {M_1}) \times \nonumber\\ &\qquad f_{{\Psi _k}|\boldsymbol{\theta}, {M_2}}({\eta _k}|\boldsymbol{\theta}, {M_1}) {\rm d}{\eta _k} \end{align} $
(21) 同理, 可通过计算$\Lambda $和${\Psi _k}$的边缘分布, 求解考虑贮存期系统性能退化影响的间歇工作系统的工作寿命和剩余工作寿命, 具体过程与式(20) 和式(21) 类似.需要说明的是, 通过求解双时间尺度下寿命和RUL分布的对其中一个时间尺度的边缘分布得到的其中一个时间尺度下设备的寿命和RUL分布, 将考虑另外一个时间尺度对设备性能退化的影响, 以及两个时间尺度之间的关系, 将会使寿命和RUL估计更准确.这一点将在第4节陀螺仪的RUL预测结果中得到进一步验证.
本节的结论都是在模型参数已知的前提下得到的, 而工程实际中往往需要根据设备的历史退化数据确定模型的参数.因此, 下一节将考虑基于设备历史观测数据的模型参数估计方法.
3. 模型参数估计
本节根据Wiener过程的性质, 基于设备的历史性能观测数据, 利用极大似然估计来对模型中的未知参数进行估计.假设有$N$个系统, 其中第$n$个系统有$m_n$次测量数据, 其第$j$次测量得到的结果和对应的时间分别为$x_{n, j}$, $t_{n, j}$和$\tau_{n, j}$.进一步, 令$\boldsymbol{x}_n^{\rm{T}} =$ $[{x_{n, 1}}$, ${x_{n, 2}}, \cdots, {x_{n, {m_n}}}]$, $\boldsymbol{x} = \{ \boldsymbol{x}_1^{\rm{T}}, \boldsymbol{x}_2^{\rm{T}}, \cdots, \boldsymbol{x}_N^{\rm{T}}\}$分别表示第$n$个设备的退化数据和所有设备的退化数据, 未知参数记为${\boldsymbol{\varphi} ^{\rm{T}}} = [{\lambda _0}, {\lambda _1}, {\lambda _2}, \sigma _B^2, \sigma _W^2]$. $\boldsymbol{x}_n$服从多元正态分布${\rm MVN}({{\boldsymbol\mu} _n}, {{\boldsymbol\Sigma} _n})$, 故其对数似然函数为
$ \begin{align}\label{EQ22} \ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x}) =& - \frac{{{m_n}}}{2}\sum\limits_{n = 1}^N {\ln (2\pi )} - \frac{1}{2}\sum\limits_{n = 1}^N {\ln \left| {{\boldsymbol{\Sigma} _n}} \right|} - \nonumber\\ &\ \frac{1}{2}\sum\limits_{n = 1}^N {{{({x_n} - {\boldsymbol{\mu} _n})}^{\rm{T}}}\boldsymbol{\Sigma}_n^{-1}({x_n} - {\boldsymbol{\mu} _n})} \end{align} $
(22) 其中, ${\boldsymbol{\mu}_n} = {{\Gamma}_n}\boldsymbol{\theta}$, ${\boldsymbol{\Sigma}_n} = {\sigma}_B^2{T}_n + \sigma _W^2{{\Lambda}_n}$,
$ \begin{align} &{{\Gamma} _n} = \left[{\begin{array}{*{20}{c}} 1&{{t_{n, 1}}}&{{\tau _{n, 1}}} \\ 1&{{t_{n, 2}}}&{{\tau _{n, 2}}} \\ \vdots&\vdots&\vdots \\ 1&{{t_{n, {m_n}}}}&{{\tau _{n, {m_n}}}} \end{array}} \right]\nonumber \\ & {{T}_n} = \left[{\begin{array}{*{20}{c}} {{t_{n, 1}}}&{{t_{n, 1}}}& \cdots &{{t_{n, 1}}} \\ {{t_{n, 1}}}&{{t_{n, 2}}}& \cdots &{{t_{n, 2}}} \\ \vdots&\vdots&\ddots&\vdots \\ {{t_{n, 1}}}&{{t_{n, 2}}}& \cdots &{{t_{n, {m_n}}}} \end{array}} \right]\nonumber \\[2mm] & {{\Lambda} _n} = \left[{\begin{array}{*{20}{c}} {{\tau _{n, 1}}}&{{\tau _{n, 1}}}& \cdots &{{\tau _{n, 1}}} \\ {{\tau _{n, 1}}}&{{\tau _{n, 2}}}& \cdots &{{\tau _{n, 2}}} \\ \vdots&\vdots&\ddots&\vdots \\ {{\tau _{n, 1}}}&{{\tau _{n, 2}}}& \cdots &{{\tau _{n, {m_n}}}} \end{array}} \right]\nonumber \end{align} $
求解使得$\ell (\left. \boldsymbol{\varphi} \right|x)$最大的参数$\boldsymbol{\varphi}$的值, 作为参数的估计值$\hat{\boldsymbol{ \varphi}}$.对其求偏导数可得
$ \begin{align}\label{EQ23} \frac{{\partial \ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}}{{\partial \boldsymbol{\theta} }} =& \sum\limits_{n = 1}^N {{\Gamma} _n^{\rm{T}}\boldsymbol{\Sigma} _n^{ - 1}({\boldsymbol{x}_n} - {{\Gamma} _n}\boldsymbol{\theta} )} \end{align} $
(23) $ \begin{align} \label{EQ24} \frac{{\partial \ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}}{{\partial \sigma _B^2}} =&\ \frac{1}{2}\sum\limits_{n = 1}^N {{({\boldsymbol{x}_n} - {\boldsymbol{\mu} _n})}^{\rm{T}}} \times\nonumber \\ &\ \boldsymbol{\Sigma} _n^{ - 1} {{T}_n}\boldsymbol{\Sigma} _n^{ - 1}({\boldsymbol{x}_n} - {\boldsymbol{\mu} _n}) -\nonumber \\ & \sum\limits_{n = 1}^N {\rm{tr}}\left(\boldsymbol{\Sigma} _n^{ - 1}{{T}_n}\right) \end{align} $
(24) $ \begin{align}\label{EQ25} \frac{{\partial \ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}}{{\partial \sigma _W^2}} =&\sum\limits_{n = 1}^N {{({\boldsymbol{x}_n} - {\boldsymbol{\mu} _n})}^{\rm{T}}} \times\nonumber \\ &\ \boldsymbol{\Sigma} _n^{ - 1} {{\Lambda} _n}\boldsymbol{\Sigma} _n^{ - 1} ({\boldsymbol{x}_n} - {\boldsymbol{\mu} _n}) -\nonumber\\ & \sum\limits_{n = 1}^N {\rm tr}\left(\boldsymbol{\Sigma} _n^{ - 1}{{\Lambda} _n}\right) \end{align} $
(25) 其中, $\rm{tr}(\cdot)$表示矩阵的迹运算.除式(23) 外, 令偏导数为零并不能得到解析的参数估计值, 因此只能采用优化算法优化似然函数$\ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})$得到参数的点估计值$ \hat{\boldsymbol{{\varphi}}}$.在本文的实例研究中, 采用的优化算法是单纯形法, 即MATLAB中的多维优化函数fminsearch.
为了得到未知参数$\boldsymbol{\varphi}$的区间估计, $\ell (\left. \boldsymbol{\varphi}\right|\boldsymbol{x})$对未知参数$\boldsymbol{\varphi}$求二阶导数, 如式(26)~(31) 所示.
$ \begin{align}\label{EQ26} &\frac{{{\partial ^2}\ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}}{{\partial \boldsymbol{\theta} \partial {\boldsymbol{\theta} ^{\rm{T}}}}} = - \sum\limits_{n = 1}^N {{\Gamma} _n^{\rm{T}} \boldsymbol{\Sigma} _n^{ - 1}{{\Gamma}_n}}\end{align} $
(26) $ \begin{align} &\frac{{{\partial ^2}\ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}}{{\partial \sigma _B^2\partial \boldsymbol{\theta} }} = \frac{{{\partial ^2}\ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}}{{\partial {\boldsymbol{\theta} ^{\rm{T}}}\partial \sigma _B^2}}= \nonumber\\ &\qquad -\sum\limits_{n = 1}^N {{\Gamma} _n^{\rm{T}} \boldsymbol{\Sigma} _n^{ - 1}{{\rm {T}}_n} \boldsymbol{\Sigma} _n^{ - 1}({\boldsymbol{x}_n} - {\boldsymbol{\mu} _n})} \\ \end{align} $
(27) $ \begin{align} &\frac{{{\partial ^2}\ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}}{{\partial \sigma _W^2\partial \boldsymbol{\theta} }} = \frac{{{\partial ^2}\ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}} {{\partial {\boldsymbol{\theta} ^{\rm{T}}}\partial \sigma _W^2}}= \nonumber\\ &\qquad - \sum\limits_{n = 1}^N {{\Gamma} _n^{\rm{T}} \boldsymbol{\Sigma} _n^{ - 1}{{\Lambda} _n} \boldsymbol{\Sigma} _n^{ - 1}({\boldsymbol{x}_n} - {\boldsymbol{\mu} _n})} \end{align} $
(28) $ \begin{align} &\frac{\partial^2\ell(\left. \boldsymbol{\varphi} \right|\boldsymbol{x})} {\partial \sigma _B^2\partial \sigma_B^2} =\frac{1}{2}\sum\limits_{n=1}^N {\rm tr}\left(\boldsymbol{\Sigma}_n^{-1} {T}_n\boldsymbol{\Sigma}_n^{-1} {T}_n \right)- \nonumber\\ &\qquad\sum\limits_{n=1}^N (\boldsymbol{x}_n- \boldsymbol{\mu} _n)^{\rm{T}}\left(\boldsymbol{\Sigma}_n^{- 1} {T}_n\boldsymbol{\Sigma}_n^{-1} {T}_n\boldsymbol{\Sigma}_n^{-1}\right)\times\nonumber\\ &\qquad(\boldsymbol{x}_n- \boldsymbol{\mu}_n) \end{align} $
(29) $ \begin{align} &\frac{{{\partial ^2}\ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}}{{\partial \sigma _W^2\partial \sigma _W^2}}=\frac{1}{2}\sum\limits_{n = 1}^N {{\rm{tr}} \left( {\boldsymbol{\Sigma} _n^{ - 1}{{\Lambda} _n} \boldsymbol{\Sigma} _n^{ - 1}{{\Lambda} _n}} \right)} - \nonumber\\ &\qquad\sum\limits_{n = 1}^N {{({x_n} - {\boldsymbol{\mu} _n})}^{\rm{T}}} \left( {\boldsymbol{\Sigma} _n^{ - 1}{{\Lambda} _n} \boldsymbol{\Sigma} _n^{ - 1}{{\Lambda} _n} \boldsymbol{\Sigma} _n^{ - 1}} \right)\times\nonumber\\ &\qquad({\boldsymbol{x}_n} - {\boldsymbol{\mu}_n}) \end{align} $
(30) $ \begin{align} &\frac{{{\partial ^2}\ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}}{{\partial \sigma _W^2\partial \sigma _B^2}}{\rm{= }}\frac{{{\partial ^2}\ell (\left. \boldsymbol{\varphi} \right|\boldsymbol{x})}}{{\partial \sigma _B^2\partial \sigma _W^2}}= \nonumber\\ &\qquad\frac{1}{2} \sum\limits_{n = 1}^N {{\rm{tr}}\left( {\boldsymbol{\Sigma} _n^{ - 1} {{\boldsymbol{ T}}_n}\boldsymbol{\Sigma} _n^{ - 1}{\boldsymbol {\Lambda} _n}} \right)} - \nonumber \\ &\qquad \frac{1}{2}\sum\limits_{n = 1}^N {{({\boldsymbol{x}_n} - {\boldsymbol{\mu} _n})}^{\rm{T}}} \times\nonumber\\ &\qquad\boldsymbol{\Sigma} _n^{ - 1} \left( {{{{T}}_n}\boldsymbol{\Sigma} _n^{ - 1} {{\Lambda} _n} + {{\Lambda} _n} \boldsymbol{\Sigma} _n^{ - 1}{{{T}}_n}} \right) \times\nonumber\\ &\qquad\boldsymbol{\Sigma} _n^{ - 1}({\boldsymbol{x}_n} - {\boldsymbol{\mu} _n}) \end{align} $
(31) 基于以上二阶微分的结果, 可知参数$\boldsymbol{\varphi}$的信息阵$\boldsymbol{\mathcal{I}}\left( \varphi \right)$为
$ \begin{align}\label{EQ32} \boldsymbol{\mathcal{I}}\left( \boldsymbol{\varphi} \right) = {\rm{diag}}\{ {{\boldsymbol{\mathcal{I}}_1} \left( \boldsymbol{\varphi} \right), {\boldsymbol{\mathcal{I}}_2}\left( \boldsymbol{\varphi} \right)} \} \end{align} $
(32) 其中, ${\rm{diag}}\{ \cdot, \cdot \} $为分块对角矩阵, $\boldsymbol{\mathcal{I}}_1(\boldsymbol{\varphi})$与$\boldsymbol{\mathcal{I}}_2(\boldsymbol{\varphi})$分别为
$ \boldsymbol{\mathcal{I}}_1\left( \boldsymbol{\varphi} \right) = \sum\limits_{n = 1}^N {{\Gamma} _n^{\rm{T}}\boldsymbol{\Sigma} _n^{ - 1}{{\Gamma} _n}} $
(33) $ {\boldsymbol{\mathcal{I}}_2}\left( \boldsymbol{\varphi} \right) = \left[{\begin{array}{*{20}{c}} \mathcal{I}_2^{11} &\mathcal{I}_2^{12} \\ \mathcal{I}_2^{11} &\mathcal{I}_2^{22} \end{array}} \right] $
(34) 其中,
$ \begin{align} \mathcal{I}_2^{11}(\boldsymbol{\varphi})&=\frac{1}{2}{\sum\limits_{n = 1}^N {{\rm{tr}}\left( {\boldsymbol{\Sigma} _n^{ - 1}{{{T}}_n}\boldsymbol{\Sigma} _n^{ - 1}{{{ T}}_n}} \right)} }\nonumber\\ \mathcal{I}_2^{12}(\boldsymbol{\varphi})&=\frac{1}{2}{\sum\limits_{n = 1}^N {{\rm{tr}}\left( {\boldsymbol{\Sigma}_n^{ - 1}{{{T}}_n}\boldsymbol{\Sigma} _n^{ - 1}{{\Lambda} _n}} \right)} }\\ \mathcal{I}_2^{21}(\boldsymbol{\varphi})&=\frac{1}{2}{\sum\limits_{n = 1}^N {{\rm{tr}}\left( {\boldsymbol{\Sigma} _n^{ - 1}{{{T}}_n}\boldsymbol{\Sigma}_n^{ - 1}{{\Lambda} _n}} \right)} }\nonumber\\ \mathcal{I}_2^{22}(\boldsymbol{\varphi})&=\frac{1}{2}{\sum\limits_{n = 1}^N {{\rm{tr}}\left( {\boldsymbol{\Sigma} _n^{ - 1}{{\Lambda} _n}\boldsymbol{\Sigma} _n^{ - 1}{{\Lambda} _n}} \right)}} \end{align} $
(35) 根据极大似然估计的性质, 当样本足够大时, 可将点估计的结果$\hat{\boldsymbol{ \varphi}}$代入Fisher信息阵中, 即$\boldsymbol{\mathcal{I}}\left( \boldsymbol{\varphi}\right)$, 进而得到参数的区间估计.具体地, 未知参数$\boldsymbol{\varphi}$的置信度为$a$的区间估计为
$ \begin{align}\label{EQ35} \left[{{{\hat \varphi }_k}-{z_\frac{a}{2}}\sqrt {\widehat {\operatorname{var} }\left( {{{\hat \varphi }_k}} \right)}, {{\hat \varphi }_k} + {z_\frac{a}{2}}\sqrt {\widehat {\operatorname{var} }\left( {{{\hat \varphi }_k}} \right)} } \right] \end{align} $
(36) 其中, $k = 1, \cdots, 5$, ${z_a} = {\Phi ^{ - 1}}(1 - a)$表示标准正态分布的$1-a$分位数, $\widehat {\operatorname{var} }\left( {{{\hat \varphi }_k}} \right)$为$\operatorname{var} \left( {{{\hat \varphi }_k}} \right)$的估计值, 而$\operatorname{var} \left( {{{\hat \varphi }_k}} \right)$可基于Fisher信息阵和Delta方法得到, 即
$ \begin{align}\label{EQ36} \operatorname{var} \left( {{{\hat \varphi }_k}} \right) = \mathcal{I}_{kk}^{-1}\left( \boldsymbol{\varphi} \right) \end{align} $
(37) 其中, $\mathcal{I}_{kk}^{-1}\left(\boldsymbol{\varphi} \right)$为$\boldsymbol{\mathcal{I}}^{-1}\left(\boldsymbol{\varphi} \right)$对角线上的第$k$个元素.
当得到参数的点估计后, 可用图形化的方法表征双时间尺度退化模型的拟合优度(Goodness of fit).具体地, 令$y_{n, j} = x_{n, j}-x_{n, j-1}$为第$n$条退化轨迹的在第$j-1$个状态监测时刻的退化增量, 易知${y_{n, j}}$服从均值和方差分别为${\mu _{n, j}} = {\lambda _1}\Delta {t_{n, j}} + {\lambda _1}\Delta {\Lambda _{n, j}}$与$\sigma _{n, j}^2 = \sigma _B^2\Delta {t_{n, j}} + \sigma _W^2\Delta {\Lambda _{n, j}}$的正态分布, 其中, $\Delta {t_{n, j}} = {t_{n, j}} - {t_{n, j - 1}}$, $\Delta {\Lambda _{n, j}} = {\Lambda _{n, j}} - {\Lambda _{n, j - 1}}$.故有如下统计量
$ {\xi _{n, j}} = \frac{1}{\sigma _{n, j}^2}\left(y_{n, j}-\mu _{n, j} \right)^2 $
(38) 服从自由度为1的$\chi _1^2$分布.因此, 可根据观测量$\boldsymbol{x}$的增量, 利用$\chi _1^2$分布的分位数--分位数图(Quantile-quantile plot, Q-Q plot)来定性地表示模型的拟合优度.单时间尺度假设下退化增量服从正态分布, 可用正态分布的Q-Q图定性地表示模型的拟合优度.当定性的模型选择标准不能满足模型选择的要求时，可以选择用模型的极大对数似然函数值$\ell (\left. \hat{\boldsymbol{\varphi}} \right|\boldsymbol{x})$或者AIC (Akaike information criterion)准则定量地比较和选择退化模型. AIC准则根据下式计算
$ \rm{AIC}= - 2\ell(\left. \hat{\boldsymbol{\varphi}} \right| \boldsymbol{x}) + 2 \it{d} (\boldsymbol{\varphi}) $
(39) 其中, $d(\boldsymbol{\varphi})$表示模型参数的个数.一般而言, 极大对数似然函数$\ell (\left.\hat{\boldsymbol{\varphi}}\right|\boldsymbol{x})$的值越大, AIC的值越小, 对应的退化模型对退化数据的拟合越好.
将本节参数估计的结果代入上节的寿命和剩余寿命分布中, 即可得到设备的寿命或者剩余寿命分布.至此, 本文实现了基于历史监测数据的设备退化建模与剩余寿命预测.下一节将本文方法应用于陀螺仪漂移系数退化数据的分析中, 对所提方法进行说明和验证.
4. 实例研究
陀螺仪是惯性平台的核心器件, 常用于火箭、导弹等长期储存一次使用的系统中.在陀螺仪的长期储存过程中, 其性能不可避免地发生退化, 具体体现为其漂移系数的增加.当陀螺仪的漂移系数超过一定阈值时, 不能再满足使用需求, 故认为陀螺仪失效.工程实际中超过90 %的陀螺仪失效为陀螺漂移系数超差[25].为了保证贮存陀螺仪的性能能够满足使用的要求, 必须定期或不定期地对陀螺仪进行通电检测.而通电检测将会加速陀螺仪的性能退化, 故累计通电时间也是衡量陀螺仪寿命的一项重要指标.因此, 在对陀螺仪进行剩余寿命预测时, 必须综合考虑储存时间$t$和检测(通电)时间$\tau$两个时间尺度.
下面利用文献[15]提供的陀螺仪漂移率历史测试数据对本文方法进行验证.该组数据包括了6组陀螺仪的漂移率数据、储存时间数据以及每次检测时长数据.如图 1所示, 陀螺仪的检测是不定期进行的, 且每次测试的时间为随机变量.基于陀螺仪测试的时间数据以及每次测试的时长数据, 利用极大似然估计即可以得到式(6) 中参数$\gamma$所服从的正态分布为N$(0.2916, {0.0564^2})$.
受篇幅限制, 本文仅给出#3陀螺仪退化轨迹的相关结果.首先, 利用所提极大似然估计方法计算模型中的未知参数, 所得结果如表 1所示.时间尺度$t$下未知参数的似然函数明显大于时间尺度$\tau$下的似然函数, 这是因为陀螺仪的性能退化主要是由于通电测试引起的.进一步对比表 1中的数据可知, 双时间尺度下退化模型中未知参数的对数似然函数更大, 而时间尺度$\tau$下AIC准则只略大于双时间尺度下的AIC.因此, 综合两个指标, 双时间尺度模型能够更好地拟合陀螺仪的退化轨迹.
基于表 1中的参数估计, 可计算双时间尺度下该陀螺仪的寿命分布, 所得结果如图 2所示.需要说明的是, 为了得到该陀螺真实的失效时间, 这里选择的是该陀螺最后一次测试的漂移值作为失效阈值, 即$w=0.2516$, 而不是文献[15]所设定的$w=0.5$.这并不影响实验的效果. 图 2中寿命的PDF考虑了两个时间尺度之间的随机关系, 故相较于单时间尺度下一维的寿命PDF, 图 2能够提供更多的寿命分布信息.
表 1 #3陀螺仪的模型参数估计结果Table 1 Estimated parameters based on degradation path #3时间尺度 $\lambda_0$ $\lambda_1$ $\lambda_2$ $\sigma_B$ $\sigma_W$ $\ell(\hat{\boldsymbol{\varphi}}|\boldsymbol{x})$ AIC $t$ -2.8959E-3 3.0662E-3 - 1.0845E-2 - 18.8689 -31.7378 $\tau$ -3.5932E-3 - 1.1095E-2 - 1.1251E-2 23.4324 -40.8648 $[t, \tau]$ 3.6148E-3 9.5126E-5 1.1440E-2 1.8146E-5 5.5160E-3 24.9530 -39.9060 此外, 利用式(20) 和式(21) 计算该陀螺仪RUL在单个时间尺度$t$下的PDF, 并利用类似的方法计算单个时间尺度$\tau$下陀螺仪的剩余寿命的分布.同时, 为了量化比较双时间尺度和单时间尺度下剩余寿命预测的准确性, 引入均方误差(Mean squared error, MSE)指标.该指标既考虑了RUL预测值的准确性, 又考虑了RUL分布的不确定性, 是RUL预测中最常用的指标.其定义为
$ \begin{align}\label{EQ39} {\rm{MSE}}_k=\int_0^\infty {{{\left({l_k} - {{\tilde l}_k}\right)}^2}} {f_{{L_k}|\boldsymbol{\theta}, {M_2}}}(\left. {{l_k}} \right|\boldsymbol{\theta}, {M_2}){\rm d}{l_k}\nonumber \end{align} $
其中, ${\tilde l_k}$为设备$t_k$时刻的真实RUL, 当RUL紧致地分布在设备真实的RUL周围时, ${\rm MSE}_k$的取值较小, RUL估计的结果也就更准确.
图 3(a)和图 3(b)分别给出了时间尺度$t$下通过双时间尺度和单时间尺度得到的设备RUL的PDF, 即剩余贮存寿命的PDF.容易看出, 通过双时间尺度得到的结果, 即考虑贮存过程中测试对性能退化的影响时, 得到的RUL的PDF更紧致, 说明预测的不确定性更小, 准确度更高.从图 4中的MSE结果, 能看出通过双时间尺度得到RUL估计的MSE明显小于通过单时间尺度得到的结果, 进一步证实了所提方法的准确性.
图 4 储存时间尺度下$M_1$和$M_2$下剩余寿命预测的MSEFig. 4 MSE comparison of $M_1$ and $M_2$ under time scale $t$在时间尺度$\tau$下, 两种方法得到陀螺仪RUL的PDF分别如图 5(a)和图 5(b)所示.整体而言, 图 5中的PDF曲线相比于时间尺度下的结果更紧致.因为陀螺仪的退化主要是由于通电测试引起的, 通电时间与退化过程的相关程度更高. 图 6中MSE远小于图 4中的结果.对比图 5(a)与图 5(b)可知, 在时间尺度$\tau$下, 通过双时间尺度得到的RUL估计更准确, 由双时间尺度得到的RUL对应的MSE明显小于通过单个时间尺度得到结果.
图 6 检测时间尺度下$M_1$和$M_2$下剩余寿命预测的MSEFig. 6 MSE comparison of $M_1$ and $M_2$ under time scale $\tau$综合以上分析结果, 基于双时间尺度退化建模的RUL预测方法有两个明显优点. 1) 该方法考虑了不同时间尺度之间的随机关系, 能提供更全面的设备寿命信息; 2) 预测结果能较方便地将转化到单时间尺度下, 且转化后的结果比基于单时间尺度退化模型的预测方法精度更高.
5. 结论
针对工程实际中一类性能变化过程与两个时间尺度相关联的随机退化设备, 本文提出了一种基于Wiener过程的双时间尺度随机退化建模与RUL估计方法.通过随机系数, 刻画时间尺度之间的不确定关系.在首达时间意义下, 推导了双时间尺度下设备寿命和RUL的分布, 给出了基于设备历史观测数据的模型参数极大估计方法.陀螺仪RUL预测的实例应用结果表明, 所提方法不仅能提供更全面的RUL分布信息, 而且预测结果能方便地转换到单时间尺度下, 并有效提高单时间尺度下RUL预测的准确性.
虽然所提出方法较之前的方法有一定的改进, 但仍有一些问题值得进一步研究. 1) 基于Wiener过程的模型一般适合对非单调退化过程进行建模, 而工程实际中也存在大量的非单调的随机退化过程, 因此，可以采用与本文类似的思路, 构建基于Gamma过程或者逆高斯过程的双时间尺度退化模型; 2) 虽然目前文献中提出了一些将非线性退化过程转化为线性化模型进行分析的方法, 但工程实际中存在一些难以线性化的退化过程, 因此, 有必要在非线性化框架下将所提的双时间尺度退化模型进行扩展; 3) 工程实际中可能既有离散的时间尺度(例如充放电次数、循环次数等), 又有连续的时间尺度(例如自然时间、储存时间等), 如何描述不同时间尺度之间的关系, 并将其融入到双时间尺度退化模型中值得进一步研究; 4) 双时间尺度的失效数据分析模型是存在的, 现实中有一些系统同时有失效数据和退化数据, 如何对在双时间尺度下实现对失效数据和退化数据的联合分析, 有待进一步考虑.
-
图 4 储存时间尺度下$M_1$和$M_2$下剩余寿命预测的MSE
Fig. 4 MSE comparison of $M_1$ and $M_2$ under time scale $t$
图 6 检测时间尺度下$M_1$和$M_2$下剩余寿命预测的MSE
Fig. 6 MSE comparison of $M_1$ and $M_2$ under time scale $\tau$
表 1 #3陀螺仪的模型参数估计结果
Table 1 Estimated parameters based on degradation path #3
时间尺度 $\lambda_0$ $\lambda_1$ $\lambda_2$ $\sigma_B$ $\sigma_W$ $\ell(\hat{\boldsymbol{\varphi}}|\boldsymbol{x})$ AIC $t$ -2.8959E-3 3.0662E-3 - 1.0845E-2 - 18.8689 -31.7378 $\tau$ -3.5932E-3 - 1.1095E-2 - 1.1251E-2 23.4324 -40.8648 $[t, \tau]$ 3.6148E-3 9.5126E-5 1.1440E-2 1.8146E-5 5.5160E-3 24.9530 -39.9060 
下载: 导出CSV
-
[1] 周东华, 魏慕恒, 司小胜.工业过程异常检测、寿命预测与维修决策的研究进展.自动化学报, 2013, 39(6):711-722 http://www.aas.net.cn/CN/abstract/abstract18097.shtmlZhou Dong-Hua, Wei Mu-Heng, Si Xiao-Sheng. A survey on anomaly detection, life prediction and maintenance decision for industrial processes. Acta Automatica Sinica, 2013, 39(6):711-722 http://www.aas.net.cn/CN/abstract/abstract18097.shtml [2] 司小胜, 胡昌华, 周东华.带测量误差的非线性退化过程建模与剩余寿命估计.自动化学报, 2013, 39(5):530-541 http://www.aas.net.cn/CN/abstract/abstract17879.shtmlSi Xiao-Sheng, Hu Chang-Hua, Zhou Dong-Hua. Nonlinear degradation process modeling and remaining useful life estimation subject to measurement error. Acta Automatica Sinica, 2013, 39(5):530-541 http://www.aas.net.cn/CN/abstract/abstract17879.shtml [3] 李鑫, 吕琛, 王自力, 陶小创.考虑退化模式动态转移的健康状态自适应预测.自动化学报, 2014, 40(9):1889-1895 http://www.aas.net.cn/CN/abstract/abstract18458.shtmlLi Xin, Lv Chen, Wang Zi-Li, Tao Xiao-Chuang. Selfadaptive health condition prediction considering dynamic transfer of degradation mode. Acta Automatica Sinica, 2014, 40(9):1889-1895 http://www.aas.net.cn/CN/abstract/abstract18458.shtml [4] Pecht M G. Prognostics and Health Management of Electronics. New Jersey:John Wiley, 2008. [5] Ye Z S, Xie M. Stochastic modelling and analysis of degradation for highly reliable products. Applied Stochastic Models in Business and Industry, 2015, 31(1):16-32 doi: 10.1002/asmb.2063 [6] Ye Z S, Chen N. The inverse Gaussian process as a degradation model. Technometrics, 2014, 56(3):302-311 doi: 10.1080/00401706.2013.830074 [7] Le Son K, Fouladirad M, Barros A. Remaining useful lifetime estimation and noisy Gamma deterioration process. Reliability Engineering and System Safety, 2016, 149:76-87 doi: 10.1016/j.ress.2015.12.016 [8] Ye Z S, Wang Y, Tsui K L, Pecht M. Degradation data analysis using Wiener processes with measurement errors. IEEE Transactions on Reliability, 2013, 62(4):772-780 doi: 10.1109/TR.2013.2284733 [9] Li N P, Lei Y G, Lin J, Ding S X. An improved exponential model for predicting remaining useful life of rolling element bearings. IEEE Transactions on Industrial Electronics, 2015, 62(12):7762-7773 doi: 10.1109/TIE.2015.2455055 [10] Zhou Q, Son J B, Zhou S Y, Mao X F, Salman M. Remaining useful life prediction of individual units subject to hard failure. IIE Transactions, 2014, 46(10):1017-1030 doi: 10.1080/0740817X.2013.876126 [11] Son J B, Zhou Q, Zhou S Y, Mao X F, Salman M. Evaluation and comparison of mixed effects model based prognosis for hard failure. IEEE Transactions on Reliability, 2013, 62(2):379-394 doi: 10.1109/TR.2013.2259205 [12] Lawless J, Hu J, Cao J. Methods for the estimation of failure distributions and rates from automobile warranty data. Lifetime Data Analysis, 1995, 1(3):227-240 doi: 10.1007/BF00985758 [13] Lawless J F, Crowder M J, Lee K A. Analysis of reliability and warranty claims in products with age and usage scales. Technometrics, 2009, 51(1):14-24 doi: 10.1198/TECH.2009.0002 [14] Singpurwalla N D, Wilson S P. Failure models indexed by two scales. Advances in Applied Probability, 1998, 30(4):1058-1072 doi: 10.1017/S000186780000879X [15] Wang Z Q, Hu C H, Wang W B, Zhou Z J, Si X S. A case study of remaining storage life prediction using stochastic filtering with the influence of condition monitoring. Reliability Engineering and System Safety, 2014, 132:186-195 doi: 10.1016/j.ress.2014.07.015 [16] Duchesne T, Lawless J. Alternative time scales and failure time models. Lifetime Data Analysis, 2000, 6(2):157-179 doi: 10.1023/A:1009616111968 [17] 王小林, 郭波, 程志君.融合多源信息的维纳过程性能退化产品的可靠性评估.电子学报, 2012, 40(5):977-982. http://www.cnki.com.cn/Article/CJFDTOTAL-DZXU201205018.htmWang Xiao-Lin, Guo Bo, Cheng Zhi-Jun. Reliability assessment of products with Wiener process degradation by fusing multiple information. Acta Electronica Sinica, 2012, 40(5):977-982 http://www.cnki.com.cn/Article/CJFDTOTAL-DZXU201205018.htm [18] 彭宝华, 周经伦, 冯静, 刘学敏.金属化膜脉冲电容器剩余寿命预测力法研究.电子学报, 2011, 39(11):2674-2679 http://www.cnki.com.cn/Article/CJFDTOTAL-DZXU201111033.htmPeng Bao-Hua, Zhou Jing-Lun, Feng Jing, Liu Xue-Min. Residual lifetime prediction of metallized film pulse capacitors. Acta Electronica Sinica, 2011, 39(11):2674-2679 http://www.cnki.com.cn/Article/CJFDTOTAL-DZXU201111033.htm [19] Si X S, Zhou D H. A generalized result for degradation model-based reliability estimation. IEEE Transactions on Automation Science and Engineering, 2014, 11(2):632-637 doi: 10.1109/TASE.2013.2260740 [20] Gebraeel N. Sensory-updated residual life distributions for components with exponential degradation patterns. IEEE Transactions on Automation Science and Engineering, 2006, 3(4):382-393 doi: 10.1109/TASE.2006.876609 [21] Whitmore G A, Schenkelberg F. Modelling accelerated degradation data using Wiener diffusion with a time scale transformation. Lifetime Data Analysis, 1997, 3(1):27-45 doi: 10.1023/A:1009664101413 [22] Ye Z S, Murthy D N P, Xie M, Tang L C. Optimal burn-in for repairable products sold with a two-dimensional warranty. IIE Transactions, 2013, 45(2):164-176 doi: 10.1080/0740817X.2012.677573 [23] Ye Z S, Chen N, Shen Y. A new class of Wiener process models for degradation analysis. Reliability Engineering and System Safety, 2015, 139:58-67 doi: 10.1016/j.ress.2015.02.005 [24] Folks J L, Chhikara R S. The inverse Gaussian distribution and its statistical application-a review. Journal of the Royal Statistical Society Series B, 1978, 40(3):263-289. http://www.jstor.org/stable/2984691 [25] 杜党波, 张伟, 胡昌华, 周志杰, 司小胜, 张建勋.含缺失数据的小波-卡尔曼滤波故障预测方法.自动化学报, 2014, 40(10):2115-2125 http://www.aas.net.cn/CN/abstract/abstract18486.shtmlDu Dang-Bo, Zhang Wei, Hu Chang-Hua, Zhou Zhi-Jie, Si Xiao-Sheng, Zhang Jian-Xun. A failure prognosis method based on wavelet-Kalman filtering with missing data. Acta Automatica Sinica, 2014, 40(10):2115-2125 http://www.aas.net.cn/CN/abstract/abstract18486.shtml 期刊类型引用(5)
1. 王康，齐金平. 基于超椭球Markov的列车控制中心剩余使用寿命预测. 铁路计算机应用. 2024(02): 67-73 . 
百度学术2. 孟海宁，童新宇，谢国，张贝贝，黑新宏. 基于RUL和SVs-GFF的云服务器老化预测方法. 自动化学报. 2024(10): 2036-2048 . 本站查看3. 裴洪，司小胜，胡昌华，郑建飞，张建勋，董青. 零寿命标签下退化设备剩余寿命预测方法. 电子学报. 2023(07): 1939-1948 . 百度学术4. 康守强，周月，王玉静，谢金宝，MIKULOVICH Vladimir Ivanovich. 基于改进SAE和双向LSTM的滚动轴承RUL预测方法. 自动化学报. 2022(09): 2327-2336 . 本站查看5. 马健，滕召胜，邱伟，马丽莎，刘颉. 基于IGA-BP神经网络的智能电能计量设备退化趋势研究. 仪器仪表学报. 2021(11): 123-131 . 百度学术其他类型引用(21)
-
下载:
计量
- 文章访问数: 1880
- HTML全文浏览量: 461
- PDF下载量: 567
- 被引次数: 26
