随着科学技术的不断进步与发展, 在实际工程中, 设备内部结构的复杂性、运行载荷和外部环境的多变性也随之不断增加, 导致设备的退化过程通常具有随机性、阶段性、多样性、相关性和混合性等特征^[1-2], 这给设备的健康管理带来了一定的难度和挑战. 为了避免由于退化失效而造成人员与财产损失, 有必要研究该类设备的退化建模与剩余寿命(Remaining useful life, RUL)预测问题, 进而保障其安全、可靠运行. 但是, 如何准确预测此类随机退化设备的剩余寿命仍是目前研究的难点与热点问题.

在实际工程中, 常常会遇到混合退化过程的设备^[1-2], 但现有文献处理方法是将设备的退化过程简化成单一的线性或非线性Wiener过程^[3-9], 可以比较容易得到设备RUL的解析解. 例如对于具有线性扩散项的维纳退化过程, Si等^[3]研究了三层不确定因素影响下线性随机退化系统的RUL预测方法; Huang等^[4]在文献[3]基础上, 提出一种带自适应漂移系数的隐线性Wiener模型, 用以描述设备的退化轨迹, 引入自适应漂移参数及测量噪声, 同时考虑设备历史数据及测量噪声问题, 进一步提高了设备RUL的预测精度; 王玺等^[5] 基于线性Wiener过程, 针对新研发光电产品提出一种剩余寿命自适应预测方法, 克服了现有RUL预测方法中, 当前时刻估计的随机参数与上一时刻随机参数的后验估计完全相等的潜在假设, 进而提高了RUL的预测准确性.

对于具有非线性扩散项的Wiener退化过程, Si等^[6]使用时空变换方法, 获得了RUL的近似解析解. 司小胜等^[7]考虑了退化过程中的测量不确定性, 然后根据维纳过程的统计性质, 建立一个状态空间模型预测RUL, 但是RUL的分布中未考虑测量不确定性. 郑建飞等^[8] 在文献[7]基础上, 还考虑了个体差异的不确定性, 并在RUL分布的推导中考虑了测量不确定性和个体差异性, 进一步提高了RUL的预测精度. Cai等^[9]考虑了退化量与监测量之间的非线性关系, 推导得到了其RUL的概率密度函数(Probability density function, PDF), 且通过激光器数据验证了该方法能够明显提高预测精度并降低预测结果的不确定性.

但是, 这些退化过程都有两个潜在的假设: 1)假设设备的退化过程是一种单一的线性或非线性退化形式, 忽略了各种退化形式在整个退化过程所占的比重是不同的, 影响效果是有差别的; 2)假设设备的退化过程是一种无记忆效应马尔科夫过程, 忽略了监测数据之间可能存在的长期依赖性和相关性对RUL的预测具有一定的影响. 例如在发动机引擎性能退化数据^[10]、高炉的性能退化数据^[11-12]、锂电池的性能退化数据^[13]等数据中发现了记忆效应, 即未来退化状态受到当前状态与历史状态的共同影响. 文献[14]在分数布朗运动(Fractional Brownian motion, FBM)的框架下, 建立了设备的退化模型, 使用复杂弱收敛定理将FBM近似为标准布朗运动(Brownian motion, BM), 然后推导了RUL的PDF. 文献[15]建立了考虑退化状态之间存在相关性的新型退化模型, 并将其应用于高炉与发动机性能退化数据中, 通过求取RUL数值PDF的方式验证了此方法的有效性和优越性. 文献[16]在文献[15]基础上, 使用更为简单的弱收敛理论推导了RUL的近似PDF, 避免大量的数值积分, 同时考虑随机效应的影响, 并通过锂电池性能退化数据进行验证, 验证结果表明, 该方法能够有效提高RUL的预测精度. 但是, 现有基于FBM过程模型的设备RUL预测方法最大的不足是均未实时更新模型中的相关参数, 仅利用同批设备的历史退化数据对模型中的未知参数进行估计. 此外, 文献[15]从结构组成和退化机理两个方面分析了惯性导航陀螺仪的随机退化过程具有混合性特征. 相比于传统单一退化形式的Wiener过程对设备的随机退化过程建模, 通过建立混合随机退化模型, 能够进一步提高RUL预测精度. 文献[15]详细论述和证明了考虑退化过程的混合性, 能够提高设备RUL的预测精度.

鉴于此, 本文将重点研究在记忆效应影响下同时包含线性退化过程和非线性退化过程的混合随机退化设备建模与RUL预测问题. 主要关注以下3个问题: 1)如何建立存在记忆效应的混合退化过程模型; 2)如何在首达时间(First hitting time, FHT)条件下, 推导设备RUL的解析表达形式; 3)如何根据监测数据实时更新退化模型参数, 并实现RUL分布的自适应更新. 针对以上3个问题, 本文首先建立了一种基于分数布朗运动的混合退化模型, 考虑了记忆效应对未来退化过程的影响; 进一步, 在弱收敛性理论和FHT的概念下, 推导了混合退化设备RUL的近似解析表达式; 然后, 利用共性参数离线估计和随机参数自适应更新的策略, 实现RUL的自适应预测; 最后, 将本文方法应用于数值仿真例子和陀螺仪的监测数据上进行验证.

1.   基于FBM的混合随机退化模型

首先, 令$ X({{t}}) $表示混合退化设备在$ t $时刻的退化量. 基于上文的具体分析, 受记忆效应影响的混合退化设备模型由以下3部分组成^[15]: 1)线性退化元件引起的线性退化过程; 2)非线性退化元件引起的非线性退化过程; 3)含有记忆效应的随机波动. 因此, 在文献[11]基础上, 在$ t $时刻, 基于FBM混合退化模型可以表示为:

式中, $ X(0)$表示混合退化设备的初始退化状态, 为了不失一般性, 假设$ X(0)= 0 $ (在实际中, 若$ X(0)\ne 0 $, 可以通过平移手段将其转化为零^[16-17]). $ \lambda t $代表混合退化设备的线性趋势, $ \int_0^t {\eta(\gamma ;\beta )} {\rm d}\gamma $代表非线性退化趋势; 令$ \lambda $和$ \alpha $为随机变量, 用于描述由结构差异、外部环境等差异引起的同批设备中不同个体差异性. $ \beta $和$ \sigma $是共性参数, 用于反映同批设备的共同特性^[18-19]. 为了刻画线性退化部分与非线性退化部分之间的相关性, 本文假设$ \lambda $和$ \alpha $服从二维正态分布. 另外, $ {B_H}(t)$为分数布朗运动, 描述了退化过程中带有记忆效应的随机波动性.

定义 1^{[13, 20]}. 赫斯特指数$ H $的取值区间满足$ 0 < H $$ < 1 $. 赫斯特指数为$ H $的分数布朗运动, $ {B_H}\left( t \right)$是一个定义在概率空间$ \left( {\Omega ,F,P} \right)$上, 且协方差函数为:

的中心化高斯过程^[21]:

式中, $ {K_H}(t - s)$定义为:

$ \Gamma(\cdot)$为伽马函数, 具体形式为:

由式(3)可以看出, FBM是BM增量的非线性移动平均值^[20], 这进一步引入了记忆效应, 其中:

为了简化$ {\sigma _B} $的计算, 假设式(1)中的$ {B_H}(t)$为标准FBM^[13]. 标准FBM满足^[21]:

1) $ {B_H}\left( 0 \right)= 0 $且${\rm{E}}\left( {{B_H}\left( 0 \right)} \right)= 0$;

2) $\forall t \ge 0,\; {\rm{E}}( {{B_H}{{\left( t \right)}^2}} )= {t^{2H}}$;

3)分数布朗运动具有自相似性和平稳增量;

4)分数布朗运动的样本轨道是连续但几乎不可微的;

5)当$ {1}/{2} < H < 1 $时, $ {B_H}(t)$具有长期依赖性、相关性, 即若令$r\left( n \right)= {\rm{cov}}( {B_H}\left( 1 \right),$${B_H}\left( {n + 1} \right)\;- {B_H}\left( n \right))$则有$\sum\nolimits_{n = 1}^{ \infty } {r\left( n \right)= \infty }$.

因此对于任意的$ \lambda \in{\bf{R}} $, $ \alpha \in {\bf{R}} $, $ \sigma \in{\bf{R}} $当且仅当$ {1}/{2} < H < 1 $时, $ X\left( t \right)$的增量具有长期相关性.

注1. FBM是一种具有长期依赖性、自相关性的连续非马尔科夫过程^[20], 其增量是固定且相关的, 并且引入了长程相关的结构^[22-23]. 赫斯特指数$ H $可以测量整个退化轨迹之间的长期依赖性. 根据$ H $的不同, FBM可以分为三种类型, 当$ 0 < H < 0.5 $时, 退化轨迹是遵循均值回归的规则; 当$H = 0.5$时, 退化轨迹是无记忆效应的随机游走过程; 当$ 0.5 < H < 1 $时, 未来的退化轨迹会沿当前的退化趋势发展.

注2. 令$\{ {{S_t},t \in \bf{R}} \}$是一个有平稳增量的随机过程且$\{ {r( n ),n \in \bf{N}{^{*}}} \}$是由$\forall n \in \bf{N}{^{*}}$, $r( n )= {\rm{E}}( {{S_{n + 1}}{S_n}} )$定义的序列, 当且仅当$\sum\nolimits_{n \in \bf{N}{^{*}}} r( n )=$$\infty$时, 就可称$\{ {{S_t},t \in \bf{R}} \}$具有长期相关性^[24].

2.   基于混合随机退化模型的RUL预测

基于FHT的定义^[25], 当失效阈值为$ \omega $时, 目标设备在任意时间$ {t_k} $的剩余寿命$ {L_k} $定义如下:

基于式(1)、式(7)和弱收敛理论^[26] 推导得到目标设备RUL的近似PDF, 如定理1所述. 定理1考虑了关于FBM的更简单的弱收敛方案, 提出了基于FBM过程的混合退化模型的RUL分布.

定理1. 对于混合退化模型(1), 若第$ m $个混合退化设备在$ {t_{m,k}} $时刻的退化状态为$ {x_m}({t_k})$, 且随机参数$ \lambda $和$ \alpha $的联合分布满足${\pi _0}\left( {\lambda ,\alpha } \right)\sim {\rm{N}}( {\mu _\lambda },{\mu _\alpha }, \sigma _\lambda ^2, \sigma _\alpha ^2,\rho )$, 在FHT的概念和弱收敛定理下, 其在 $ {t_{m,k}} $时刻估计的剩余寿命${l_{m,k}}$的近似PDF如下:

式中,

式中, $ {\omega _{m,k}} = {\omega _m} - {x_m}\left( {{t_k}} \right)$, $ \tilde h = h\left( {{t_k} + {l_k}} \right)- h\left( {{t_k}} \right)$, $ \bar h = h\left( {{t_k} + {l_k} + \Delta l} \right)- h\left( {{t_k} + {l_k}} \right)$.

定理1的证明见附录A.

下面对混合退化模型中的未知参数进行实时估计更新, 从而实现混合退化设备RUL的自适应预测.

3.   随机混合退化模型的参数估计

3.1   共性参数与随机参数的离线估计

令$\{ {{X_m}( {{t_{m,n}}}),\;1 \le m \le M,\;1 \le n \le {N_{{m}}}}\}$为通过状态监测的手段获得的同一批次随机退化设备的历史退化数据, 其中$ M $为随机设备的总个数, $ {N_m} $为第$ m $个随机退化设备监测数据的个数, $ {t_{m,n}} $为第$ m $个随机退化设备的第 $ n $个观测时刻. 本文将${X_m}( {{t_{m,n}}})$记作$ {X_{{{m}},n}} $. 采用极大似然估计算法, 利用同批次随机退化设备的历史退化监测数据, 对随机退化模型的共性参数和随机参数先验分布中的超参数进行离线估计. 令$\Theta = [ {\mu _{\lambda 0}},{\mu _{\alpha 0}},\sigma _{\lambda 0}^2,\sigma _{\alpha 0}^2, {\rho _0}, \beta , {\sigma ^2},H]$为带估计的模型参数, 其中$ {\mu _{\lambda 0}} $和$ \sigma _{\lambda 0}^2 $为随机参数$ \lambda $先验分布中的超参数, $ {\mu _{\alpha 0}} $和$ \sigma _{\alpha 0}^2 $为随机参数$ \alpha $先验分布中的超参数.

根据式(1), 第$ m $个随机退化设备在$ {t_{m,n}} $时刻的退化状态可表示为:

式中, $ {\lambda _0} $与$ {\alpha _0} $分别表示随机参数$ \lambda $和$ \alpha $的先验值, 将其联合先验分布记作${\pi _0} \sim {\rm{N}}( {\mu _{\lambda 0}},$${\mu _{\alpha 0}},\sigma _{\lambda 0}^2, \sigma _{\alpha 0}^2, {\rho _0})$.

为了方便推导, 令:

假设同批次不同设备个体之间的历史退化数据都是互不相关的. 因此, 对第$ m $个随机退化设备的所有历史监测数据可以进一步写为:

定理2. $ {X_m} $的联合分布为多元正态分布, 其均值与方差如下:

定理2证明见附录B.

根据定理2可以得到模型参数$ \Theta $的似然函数如下:

式(19)对$ {\mu _{\lambda 0}} $和$ {\mu _{\alpha 0}} $分别求一阶偏导, 可得:

令式(20)、式(21)分别为零, 可得:

其中

将式(22)、式(23)代入式(19), 可得:

可以看出, 似然函数(24)具有高维的特征, 直接把式(24)极大似然化很难得到其余参数的极大估计值. 本文首先利用Matlab中的fminsearch函数, 求取极大似然估计值(该函数基于Nelder-Mead单纯形法对最小化执行同步多维搜索), 进而得到参数$ \sigma _{\lambda 0}^2 $、$ \sigma _{\alpha 0}^2 $、$ {\rho _0} $、$ \beta $、$ {\sigma ^2} $和$ H $的估计值; 然后, 将其代入式(22)、式(23), 得到相应$ {\mu _{\lambda 0}} $、$ {\mu _{\alpha 0}} $的极大似然估计值; 最后, 采用贝叶斯推理方法, 利用设备的实时退化监测数据对随机参数进行实时更新.

3.2   随机参数的实时更新

第3.1节通过参数离线估计的方法得到了随机参数$ \lambda $和$ \alpha $的联合先验分布, 下面对随机参数进行实时更新. 假设第$ m $个设备在前$ {t_k} $时刻一共获得$ k $个退化监测数据$ X_k^* $, 记作:

式中, $ {t_{{i}}}\left( {i = 1,2, \cdots ,k} \right)$表示对应的退化监测时刻.

基于贝叶斯推理理论, 利用获得的退化数据$ X_k^* $更新$ \lambda $和$ \alpha $的联合后验分布, 即:

式中, $ {\pi _0}\left( {\lambda ,\alpha } \right)$是$ \lambda $和$ \alpha $的联合先验分布, $ X_k^* $服从多元正态分布即$X_k^* \sim {\rm{N}}\left( {{\mu _k},{\varSigma _{k}^{-1}}} \right)$, $ \mu _k $和$\varSigma _{k}^{-1}$的具体形式参考式(14)、式(15). 因此, 在已知$ \lambda $和$ \alpha $情况下, 退化监测数据$ X_k^* $的联合分布为:

将估计得到的联合分布$ {\pi _0}\left( {\lambda ,\alpha } \right)$、式(27)代入式(26), 对$ \lambda $和$ \alpha $的联合后验分布进行更新, 具体更新结果如下:

其中

$ Z $为与$ \lambda $和$ \alpha $无关的常数项. 式(28)中, $X_k^{*\mathrm{T}}\;\times \varSigma _{k}^{-1}{T_A}$ 和 $T_A^\mathrm{T}\varSigma _{k}^{-1}X_k^*$、$X_k^{*\mathrm{T}}\varSigma _{k}^{-1}{T_B}$ 和 $T_B^\mathrm{T}\;\times \varSigma _{k}^{-1} X_k^*$、$T_A^\mathrm{T}\varSigma _{k}^{-1}{T_B}$和$T_B^\mathrm{T}\varSigma _{k}^{-1}{T_A}$都是常数值, 并且互为转置矩阵, 由矩阵的性质可以得到$X_k^{*\mathrm{T}} \varSigma _{k}^{-1}{T_A}= T_A^\mathrm{T}\;\times \varSigma _{k}^{-1}X_k^*$、$X_k^{*\mathrm{T}}\varSigma _{k}^{-1}{T_B}= T_B^\mathrm{T}\varSigma _{k}^{-1}\times {X_k^*}$、$T_A^\mathrm{T}\varSigma _{k}^{-1}{T_B}= T_B^\mathrm{T}\varSigma _{k}^{-1}{T_A}$.

因此, 式(28)可以进一步简化为:

式中, $ Z_1 $是与$ \lambda $和$ \alpha $无关的常数项. 本文假设$ \lambda $和$ \alpha $的联合先验分布是服从二元正态分布的, 根据正态分布的特性, $ \lambda $和$ \alpha $的联合后验分布也是服从二元正态分布的, 满足$\left\{ {\lambda ,\alpha |X_k^*} \right\} \sim {\rm{N}}( {\mu _{\lambda k}},{\mu _{\alpha k}},\sigma _{\lambda k}^2, \sigma _{\alpha k}^2,{\rho _k} )$, 因此$ p\left( {\lambda ,\alpha |X_k^*} \right)$的具体形式为:

基于式(29)、式(30), $ \lambda $和$ \alpha $在$ {t_k} $时刻联合后验分布中的相关参数可以由下式得到:

式中, $F = T_A^\mathrm{T}\varSigma _{k}^{-1}{T_A}$, $G = T_B^\mathrm{T}\varSigma _{k}^{-1}{T_B}$, $H = T_A^\mathrm{T}\;\times $ $\varSigma _{k}^{-1}{T_B}$, $I = - T_A^\mathrm{T}\varSigma _{k}^{-1}X_k^*$, $J = - T_B^\mathrm{T}\varSigma _{k}^{-1}X_k^*$.

当在$ t_k $时刻获取到最新的观测数据$ X_k^*\left( {{t_k}} \right)$后, 可以基于式(31) ~ (35), 对$ \lambda $和$ \alpha $的超参数进行实时更新.

将第3.1节参数的离线估计与第3.2节参数的实时更新归结为算法1. 至此, 完成了模型里面随机参数的实时更新. 下面将本文方法应用到数值仿真和陀螺仪实际退化监测数据中.

　 算法1. 参数的离线估计与实时更新算法

1)参数离线估计阶段

输入. ${\tilde X_{1:M}}$, $1 \le m \le M$, ${t_{m,n}}$, $1 \le {{n}} \le {N_{{m}}}$.

输出. ${t_k}$时刻${\mu _{\lambda k}}$、${\mu _{\alpha k}}$、$\sigma _{\lambda k}^2$、$\sigma _{\mu k}^2$、${\rho _k}$的实时更新值.

a)初始化模型参数$\Theta$;

b)最大化式(24), 得到参数$\sigma _{\lambda 0}^2$、$\sigma _{\alpha 0}^2$、${\rho _0}$、$\beta$、${\sigma ^2}$、$H$ 的极大似然估计值;

c)将b)中得到的参数极大似然估计值代入式(22)、式(23), 得到${\mu _{\lambda 0}}$和${\mu _{\alpha 0}}$的估计值;

d)输出模型参数$\Theta$的估计值.

2)随机参数的实时更新

输入. $X_k^* = {\left[ {X_k^*\left( {{t_1}} \right),X_k^*\left( {{t_2}} \right), \cdots ,X_k^*\left( {{t_k}} \right)} \right]^\mathrm{T}}$, ${t_{{i}}}\;( i = 1, 2, \cdots, k)$和参数离线估计阶段获得的离线估计值$\Theta$.

a)在任一时刻${t_k}$, 计算出A、B、C、D、E、F、G、H、I 、J;

b)将a)得到的A、B、C、D、E、F、G、H、I 、J值代入式(31) ~ (35), 可以得到${\mu _{\lambda k}}$、${\mu _{\alpha k}}$、$\sigma _{\lambda k}^2$、$\sigma _{\mu k}^2$、${\rho _k}$在${t_k}$时刻的更新值.

4.   数值仿真与实例验证

本节将本文所提方法应用到数值例子与实际例子中, 验证其是否有效. 采用线性退化模型、非线性退化模型和本文提出的模型来拟合退化数据, 并且比较在三种模型下设备RUL的预测结果. 三种模型为: 1)本文设计的模型; 2)模型1. 基于带线性漂移 的Wiener过程模型^[5]; 3)模型2. 基于带非线性漂移的Wiener过程模型^{[8, 27]}; 4)模型3. 基于带非线性漂移的FBM过程模型(单一退化形式)^[13-14] (未考虑随机系数的实时更新).

4.1   数值仿真

1)仿真数据

首先, 根据混合随机退化模型$X({{t}})= \lambda t \;+ \alpha \int_0^t {\eta(\gamma ;\beta )} {\rm{d}}\gamma + \sigma {B_H}(t)$产生所需的退化数据, 其中假设$ \eta(\gamma ;\beta )= \beta {\gamma ^{\beta - 1}} $. 由于使用幂函数来刻画设备非线性退化部分被广泛应用于退化建模领域, 则该模型可简化为$ X({{t}})= \lambda t + \alpha {t^\beta } + \sigma {B_H}(t)$. 根据模型假设$(\lambda ,\alpha )\sim {\rm{N}}\left( {{\mu _\lambda },{\mu _\alpha },\sigma _\lambda ^2,\sigma _\alpha ^2,\rho } \right)$, 设置仿真参数: $ {\mu _\lambda }=1 $, $ {\mu _\alpha }=0.265 $, $ \sigma _\lambda ^2=4.5 \times {10^{{\rm{ - }}4}} $, $\sigma _\alpha ^2= 3.35 \;\times {10^{{\rm{ - }}4}}$, $ \rho = - 0.468 $, $ \sigma = 0.05 $, $ \beta = 0.2 $, $H = 0.87$. 同时, 设置仿真步长$ \tau = 0.1 $. 使用小波合成的方法^[23]得到标准分数布朗运动$ {B_H}(t)$. 最后, 利用Matlab生成30组退化数据作为同批次设备的历史退化数据. 仿真历史退化数据如图1所示.

图 1  30组仿真历史退化数据

Fig. 1  30 sets of simulated history degradation data

下载: 全尺寸图片 幻灯片

图1是30组仿真退化轨迹, 仿真监测时间为0 ~ 15 s, 采样间隔为0.1 s, 每条轨迹共150个采样点. 图2为其中任意一条退化轨迹. 作为带预测设备的实时监测数据, 为了简单化, 将该轨迹的最后一个监测数据$ X(15)= 13.1 $作为失效阈值, 即$ \omega = 13.1 $. 则该设备的剩余寿命可近似为14.8 s.

图 2  待预测设备实时监测退化数据

Fig. 2  Real-time monitoring degradation data of the equipment for prediction

下载: 全尺寸图片 幻灯片

下面利用图2的仿真退化数据对本文所提模型和RUL预测的有效性进行验证.

2)结果对比

为了比较的公平性, 模型2、模型3中的非线性部分也采用幂函数的形式. 首先, 基于历史退化数据, 采用极大似然估计算法得到共性参数与随机参数的先验估计值. 为了对比度量四种模型之间的拟合的精准度与估计的准确度, 使用赤池信息准则(Akashi information criterion, AIC)^[28]、贝叶斯信息准则(Bayesian information criterion, BIC)^[29]测量预测模型的拟合程度. AIC、BIC的值越小, 则拟合程度越高; 反之, 则拟合越差. AIC的具体公式为:

式中, $ L\left( \Theta \right)$表示似然函数值, $ p $是未知参数的总个数.

BIC引入了改进惩罚项, 惩罚项大于AIC的, 可以有效避免大样本的过拟合问题, 其具体表达式为:

式中, $ L\left( \Theta \right)$表示似然函数值, $ p $是未知参数的总个数, $ n $为样本数据量.

表1为四种模型参数(本文方法、模型1、模型2、模型3)的先验估计值.

表 1  四种模型参数的先验估计值

Table 1  The parameters＇ prior estimates ofthe four models

参数	本文方法	模型1	模型2	模型3
${\mu _\lambda }$	1.0754	2.0251	—	—
${\mu _\alpha }$	0.2547	—	3.1547	3.3643
$\sigma _\lambda ^2$	$4.56\times10^{-4}$	0.1618	—	—
$\sigma _\alpha ^2$	$3.47\times10^{-4}$	—	0.0176	0.0103
$\rho$	−0.4793	—	—	—
${\sigma ^2}$	0.53362	0.51362	0.58362	0.54112
$\beta$	0.1547	—	0.6035	0.6101
$H$	0.8472	—	—	0.8130
$\ln L\left( \Theta \right)$	57.462	46.328	48.633	51.278
${\rm{AIC}}$	−98.924	−86.656	−89.266	−92.556
${\rm{BIC}}$	−87.724	−82.456	−83.666	−85.556

下载: 导出CSV 

| 显示表格

由表1可知, 本文方法的AIC与BIC值最小, 并且其中的3个方差参数总体也比较小, 这说明本文模型参数估计的准确性更高. 主要原因是本文所提模型既考虑退化过程的混合性, 还考虑了退化状态之间可能存在的相关性.

得到参数的离线估计值后, 基于待预测设备的实时监测退化数据, 采用贝叶斯更新的方法实现$ \lambda $和$ \alpha $的实时更新, 进而实现RUL的实时预测. 图3是本文模型中$ \lambda $和$ \alpha $的实时更新过程.

图 3  随机参数的实时更新过程

Fig. 3  Real-time updating process of random parameters

下载: 全尺寸图片 幻灯片

由图3可以看出: 1)随着仿真退化数据的不断获取, 随机参数不断实时更新; 2) $ \rho $的实时更新值为负, 表明$ \lambda $和$ \alpha $是负相关, 进一步表明混合随机退化设备的线性退化部分与非线性退化部分是相互抑制的.

下面对比在四种方法下得到的RUL实时预测结果. 预测时间点共取10个, 从第5个时间点(第50个采样点)到第14个时间点(第140个采样点), 间隔为10个采样点, 具体RUL对比见图4. 同时选取第5、8、11、13个时间点的RUL预测结果, 通过二维平面图形式做进一步对比, 对比图如图5所示.

图 4  4种方法在各个时间点处RUL预测的对比

Fig. 4  Comparison of RUL prediction by four methods at each time

下载: 全尺寸图片 幻灯片

图 5  四种方法在第5、8、11、13个时间点RUL预测结果的对比

Fig. 5  Comparison of RUL prediction by four methods at 5, 8, 11, 13th time

下载: 全尺寸图片 幻灯片

在图4和图5中, 蓝色曲线为本文方法RUL的PDF, 红色曲线为模型1得到的PDF, 黑色曲线为模型2得到的PDF, 绿色曲线为模型3得到的PDF. 可以看出: 1) 四种模型随着获取的退化数据增多, RUL预测的精度越来越高; 2)本文方法得到的RUL的PDF能够更好地覆盖真实的RUL, 其预测均值更接近真实的RUL, 较模型1、模型2、模型3得到的PDF精度更高; 3)本文方法预测得到的RUL的PDF更为尖锐和紧凑, 这说明本文方法预测的不确定性比其他三种方法预测的不确定性更小. 由图4、图5可以直观地看出, 本文方法优于其他三种方法.

下面从定量的角度分析四种方法的优劣. 使用均方误差(Mean squared error, MSE)指标来评价四种方法预测RUL的精度, 其既可以表征RUL预测结果的准确性, 还可以表征RUL预测的不确定性, 是常用的误差评定指标, 可以较好地对比不同预测方法的优劣. 待测设备在$ {t_k} $的RUL的MSE可以表示为:

所有预测时间点处的均方误差之和为总体均方误差(Total mean squared error, TMSE). 很明显, TMSE值越小, 则该方法的RUL预测精度就越高. 通过计算得到四种方法在第10个时间点处的TMSE分别为1.0051$ \times 10^{2} $、5.1308$ \times 10^{3} $、$1.1007\times 10^{3} $、1.5041$ \times 10^{2} $. 通过对比可知, 本文方法的RUL预测精度高于传统模型1、模型2、模型3的RUL预测精度; 相比于模型1和模型2, 精度约提高了一个数量级. 定量分析结果与图4、图5直观得到的结论是一致的.

4.2   实例验证

陀螺仪是惯性导航系统、导弹制导与控制系统的关键设备, 决定着导弹导航与命中的精度. 但是, 随着陀螺仪工作年限的增加, 在外部复杂环境与内部随机应力的影响下, 陀螺仪的性能可能会随之发生退化, 主要表现在其漂移系数的不断增大, 如果漂移系数值增大到一定程度, 陀螺仪就无法正常工作, 即陀螺仪发生失效.

本文获取的5组某型号陀螺仪的退化监测数据^[30], 每组73个数据, 采样间隔为2.5 h, 陀螺仪失效阈值设定为0.37 (°/h), 陀螺仪寿命约为180.5 h. 首先, 基于极大似然估计算法, 利用其中4组退化数据^[30]对模型参数进行离线估计, 得到其先验估计值; 其次, 利用剩下一组数据为待预测的陀螺仪的实时监测数据, 采用贝叶斯更新的方法, 对模型随机参数进行实时更新; 最后, 进行设备RUL的实时预测, 并且将模型1、模型2、模型3得到的RUL进行比较. 图6为待预测陀螺仪的实时监测数据.

图 6  待预测陀螺仪漂移数据

Fig. 6  Drift data of the gyroscope for prediction

下载: 全尺寸图片 幻灯片

为了比较公平, 模型2、模型3中的非线性形式和本文模型的非线性部分都采用幂函数形式. 调用算法1中的离线估计阶段得到模型参数的先验估计值, 如表2所示.

表 2  陀螺仪退化模型参数的先验估计值

Table 2  A parameters＇ prior estimate of the gyroscope degradation model

参数	本文方法	模型1	模型2	模型3
${\mu _\lambda }$	0.0151	0.1348	—	—
${\mu _\alpha }$	0.0373	—	0.1503	0.1671
$\sigma _\lambda ^2$	$2.34\times 10^{-6}$	0.0028	—	—
$\sigma _\alpha ^2$	$7.82 \times 10^{-6}$	—	$5.76\times 10^{-4}$	$4.76\times 10^{-4}$
$\rho$	$-0.4503$	—	—	—
${\sigma ^2}$	0.00192	0.00582	0.00422	0.00272
$\beta$	0.2114	—	0.1342	0.1742
$H$	0.8011	—	—	0.7903
$\ln L\left( \Theta \right)$	41.871	31.013	34.337	38.699
${\rm{AIC}}$	−67.742	−56.026	−60.674	−67.398
${\rm{BIC}}$	−70.942	−57.226	−62.274	−69.398

下载: 导出CSV 

| 显示表格

下面基于陀螺仪模型参数的先验估计值, 调用算法1中的参数实时更新阶段, 在每一个时间点对模型随机参数进行实时更新, 模型随机参数的实时更新过程如图7所示. 由图7可知, 随着陀螺仪监测数据获取的增加, 模型中随机参数不断实时更新, 并趋于稳定状态, 这表明随机参数的估计值越来越接近真实的模型参数值.

图 7  陀螺仪退化模型随机参数的实时更新过程

Fig. 7  Real-time updating process of random parameters of gyroscope degradation model

下载: 全尺寸图片 幻灯片

下面选择第152.5 h ~ 175 h之间的10个时间点为例, 比较分析所使用四种方法对该型号陀螺仪RUL预测的精准程度. 在所选的10个时间点, 四种方法预测该型号陀螺仪RUL的PDF如图8所示. 由图8可知: 1)四种方法预测得到RUL的PDF在10个时间点都能够比较好地覆盖真实的RUL, 说明了四种方法能够有效地预测该型号陀螺仪的RUL. 2)随着监测数据的增加, 参数不断实时更新, 四种方法预测RUL的PDF随着时间变得越来越尖、越来越窄, 说明四种方法对RUL预测的不确定度随着监测数据的累积而不断减小. 3)本文方法得到RUL的PDF与其他三种方法相比较, 其RUL的PDF明显更高、更紧致, 说明本文方法预测得到的RUL更准确, 并且预测的不确定度更小. 其原因是, 相比于模型1和模型2单一的马尔科夫建模方式, 本文方法不仅考虑了陀螺仪内部不同类型的随机退化过程, 还考虑了退化监测数据之间的长期相关性、记忆性, 克服了模型1和模型2中的强马尔科夫限制. 相比于模型3, 本文方法考虑了设备的退化混合性, 更符合设备真实的退化过程.

图 8  在不同时间点处四种方法预测RUL的PDF对比

Fig. 8  Comparison of RUL＇s PDFs by four prediction methods at each time

下载: 全尺寸图片 幻灯片

进一步通过计算, 得到四种方法在10个时间点的TMSE分别为${1.8154} \times 10^{2}$、$2.6321 \times {10}^2 $ 、${3.6543}\times $${10^2} $、${2.1031}\times 10^{2}$. 可知, 本文方法的RUL预测精度高于其他三种方法.

5.   结束语

本文针对随机退化设备, 提出一种考虑随机退化过程中退化状态之间可能存在长期相关性、记忆性特点的混合随机退化模型, 克服了目前研究中潜在的假设: 1)假设设备的随机退化过程为单一的线性或非线性形式, 忽略了设备内部的多种类型退化过程分别对RUL预测的影响; 2)将设备的随机退化过程简化成无记忆效应的马尔科夫过程, 受到强马尔科夫性的限制. 最后, 通过本文方法实现了随机退化设备的RUL寿命自适应预测, 且预测精度高于同等条件下传统模型的预测精度. 主要结论有一下3点:

1)基于FBM的混合随机退化模型, 能够更为准确、合理地描述设备随机退化过程中可能含有的记忆效应. 通过实例验证和与传统的建模方法比较, 其RUL预测结果更加精确、不确定性更小;

2)本文采用的参数估计方法有效地利用了设备的历史退化信息和实时监测信息, 从而更加准确地实现了设备RUL的实时预测.

综上所述, 本文提出的基于FBM的混合随机退化模型对随机退化设备建模更为合理与准确, 并且RUL预测结果优于传统方法, 具有一定的工程实用价值. 下一步的研究方向为如何对预测模型中共性参数进行自适应更新.

附录A 定理1的证明

证明. 当$ n \to \infty $时, 缩放过程$ \left\{ {{z_{n,{H_p}}}(t),t \ge 0} \right\} $弱收敛于$ \sqrt {{L_p}} {B_{{H_p}}}(t)$^[26], 其中:

$ {z_{n,{H_p}}}(t)$的具体形式为:

式中, $\left\lfloor nt \right\rfloor$表示不超过的最大整数, 即$\left\lfloor nt \right\rfloor$为向下取整符号. $ \left\{ {{\xi _i},i \in Z} \right\} $是独立同分布的随机序列, 满足${\rm{E}}\left( {{\xi _0}} \right)= 0$和$\rm{E}\left( {\xi _0^2} \right)= 1$. 当$ k > 0 $时, $ {a_k} = p{k^{p - 1}} $; 否则, $ {a_k} = 0 $.

因此, 只需要解决以下形式的退化过程的RUL分布:

其中

基于弱收敛理论, 可以将模型(1)转换成式(37). 文献[27]指出, 时间重新缩放的FBM保持零均值高斯过程, 然后在适当的假设下给出了相应RUL分布的近似表达式^[13]. 具体来说, 考虑以下退化模型:

假设失效阈值为$ \omega $, 则第$ m $个设备在$ {t_k} $时刻RUL的PDF如下:

式中, $\Delta \Lambda ( {{t_k} \;+\; {l_k};\theta } )\;= \;\Lambda ( {{t_k}\; + \;{l_k};\theta } )- \Lambda ( {{t_k};\theta } )$, $\Delta \Psi ( {t_k} + {l_k};\vartheta )\,=\, \Psi ( {{t_k}\; +\; {l_k};\vartheta } )- \Psi ( {{t_k};\vartheta } )$, $\Delta \tau ( {t_k}\; +$ $ {l_k};\varepsilon )=\tau ( {t_k} + {l_k};\varepsilon )- \tau ( {{t_k};\varepsilon } )$, $ {\omega _k} = \omega - x( {{t_k}} )$.

令$\Lambda ( {{t_k};\theta } )= {t_k},$$\Psi ( {{t_k};\vartheta } )= \int_0^{{t_k}} {\eta ( {\gamma ;\beta } )} {\rm{d}}\gamma,$ $\Delta \Lambda\;\times$ $ ( {t_k} + {l_{m,k}}; \theta )= {l_{m,k}}$, ${\omega _{m,k}} \;= \;{\omega _m}\; -\; {x_m}( {{t_k}} )$, $\Delta \Psi\;\times$ $( {{t_k} + {l_{m,k}};\vartheta } ) = \int_{{t_k}}^{{t_k} + {l_{m,k}}} {\eta ( {\gamma ;\beta } )} {\rm{d}}\gamma$, $\Delta \tau ( {t_k} + {l_{m,k}};\varepsilon )= h( {{t_k}+ {l_{m,k}}} )- h( {{t_k}} )= \Delta h( {{t_k} + {l_{m,k}}} ),$ 则:

令${\xi _{m,k}} = \left( {{l_{m,k}}\sigma _\lambda ^2 + \Delta \Psi \left( {{t_k} + {l_{m,k}};\vartheta } \right)\sigma _\alpha ^2} \right){\omega _k} \;+$ $( {\mu _\lambda } +{\mu _\alpha } ){\sigma ^2}\Delta h\left( {{t_k} + {l_{m,k}}} \right),$ $ 则\;{g_{m,k}} \;可以简化为:$

定理1的具体使用方法见文献[13]. □

附录B 定理2的证明

证明. 基于式(14)和FBM的特性, 可以比较容易得到${\mu _{{m}}}= {\mu _{\lambda 0}}{T_{m,A}} + {\mu _{\alpha 0}}{T_{m,B}}$. 下面对$\varSigma _{m}^{-1}$进行推导.

对于监测数据${X_m} = [ {X_{m,1}},{X_{m,2}}, \cdots ,X_{m,{N_{{m}}}} ]^{\rm{T}}$, 其协方差阵可以记为如下:

根据式(14), 可得:

令$ {t_{m,1}} \le {t_{m,{n_i}}} \le {t_{m,{n_j}}} \le {t_{m,n}} $. 根据协方差定义, 可得:

首先, 推导${\rm{E}}\left( {{X_{m,{n_i}}},{X_{m,{n_j}}}} \right)$:

其次, 推导${\rm{E}}\left( {{X_{m,{n_i}}}} \right){\rm{E}}\left( {{X_{m,{n_j}}}} \right)$:

将式(B4)、式(B5)代入式(B3), 可得:

最后, 将式(B6)、式(B2)代入式(30), 可得${\sum _{m=1}^M}$. □