考虑随机冲击影响的自适应Wiener过程剩余寿命预测方法

董青,郑建飞,*,胡昌华,余铜辉,牟含笑

1. 火箭军工程大学 导弹工程学院,西安 710025 2. 火箭军驻西安地区第一军事代表室,西安 710100

随着高新技术的日益发展,现代工业设备正朝着高度自动化、集成化发展。由于寿命周期中受到振动冲击、工况切换、机械磨损、化学侵蚀、负载变化以及能量消耗等内外部因素影响的综合作用,此类设备的性能与健康状态不可避免的出现退化现象,最终导致设备失效,甚至造成难以挽回的生命财产及环境损失。如果能在设备失效之前预测其寿命或剩余寿命,进而采用相应的维护策略,则可有效避免灾难性事故,从而达到降低成本和保护环境的目的。作为预测与健康管理(Prognostics and Health Management, PHM)的核心技术,剩余寿命预测(Remaining Useful Life, RUL)是对设备进行科学维修决策的前提和基础。因此,准确预测设备剩余寿命来降低突发失效风险具有重要的理论意义和工程应用价值。

经过几十年的发展,剩余寿命预测技术取得了丰硕的理论成果并得到广泛应用。Pecht在关于PHM的专著中,将剩余寿命预测方法主要分为基于机制建模和基于数据驱动建模的方法2类。柴天佑指出复杂的工业过程往往具有多变量、强耦合、强非线性、动态工况变化、难以用数学模型描述等复杂特性,数据驱动的方法主要包括人工智能方法和统计数据驱动的方法。Si等将数据驱动的方法分为基于直接监测数据的方法和基于间接监测数据的方法。其中直接监测数据的方法又可分为基于随机系数回归模型的方法、基于Wiener过程的方法、基于Gamma过程的方法、基于马尔科夫链的方法以及逆高斯过程的方法。相对于其他随机过程,Wiener过程是一类具有高斯独立分布增量的非单调退化过程。凭借其良好的数学特性,已被广泛应用于设备可靠性分析和寿命预测领域。

在工程实际中,设备除去正常工作外,还会受到单粒子效应、温度振动冲击、腐蚀等外部环境因素的影响,而这些额外因素可以看作广义的随机冲击,对设备退化过程造成一定影响。张延静等针对受到随机冲击影响的单部件系统,考虑冲击对退化量的影响,建立竞争失效可靠性模型,并运用到维修决策中。孙富强等针对存在冲击韧性的系统,考虑冲击对退化量和退化率2方面的影响,建立非线性Wiener过程的竞争失效模型进行可靠性评估。王浩伟等针对导弹武器系统,假定突发失效概率与整弹性能退化相关,采用Gamma过程和Weibull分布分别作为退化失效模型和突发失效模型进行可靠性评估和RUL预测。王华伟等利用航空发动机失效信息,分别对性能退化失效和突发失效建立RUL预测模型。通过建立混合Weibull可靠性模型,量化性能退化失效对突发失效的影响,实现突发失效退化情况下的航空发动机剩余寿命预测。白灿等提出一种考虑随机冲击影响的非线性退化设备RUL预测方法,在首达时间意义下给出寿命和剩余寿命的近似解表达式,极大缩短了计算时间,为后续开展设备健康维护奠定了良好基础。

基于上述分析,考虑冲击与退化之间的相互影响已经开展许多工作,但大多应用于设备可靠性评估领域,而在寿命预测领域仍需要进一步深入拓展。在现有基于Wiener过程进行退化建模和剩余寿命预测方法研究中,对于设备寿命周期中现实存在的同时具有随机冲击影响、监测间隔不均匀、与先验数据测量频率不一致的情况尚未考虑。同时,对于模型漂移系数的更新,仅局限于实时监测数据下的更新,而未考虑未来RUL预测中该模型自适应漂移的可变性。

综上,本文考虑随机冲击对退化过程的影响,提出一种基于自适应Wiener过程的非线性退化建模和剩余寿命预测方法。首先,通过自适应Wiener过程描述自身退化过程；其次,通过正态分布刻画单次冲击对设备退化水平的影响,并将随机冲击融入到自身退化建模中,在首达时间意义下,推导出设备随机冲击影响下剩余寿命分布解析式；再次,通过构建状态空间模型,基于Kalman滤波框架和期望最大化算法,利用退化数据实现参数和RUL的自适应更新；最后,通过数值仿真、惯导系统陀螺仪以及锂电池实例,验证所提方法进行剩余寿命预测的精确性和实用性。

1 考虑随机冲击影响的自适应Wiener过程退化建模

1.1 自适应Wiener过程退化模型

令()表示设备在运行中的退化过程,现有常用的Wiener过程模型为

()=++()

(1)

式中：、分别为退化过程的漂移系数和扩散系数；()为标准Brownian运动；为退化初值,通常假设=0。

由于上述Wiener过程模型存在测量间隔不均匀、测量频率不一致、忽略自适应漂移可变性3点不足,当不考虑随机冲击影响时,对于设备退化过程{(),≥0},采用自适应Wiener过程来描述

(2)

式中：()为一个符合Wiener过程且随时间变化的漂移系数;为初始漂移率,>0;为自适应漂移项的扩散系数;()为设备监测的退化量;为退化过程的扩散系数;()为1个独立于()的标准Brownian运动;为时间的1个参数；(,)为1个随时间增长的非线性函数,为了方便表示,后续可简写为()。当=0,=1时,设备呈线性退化,此时退化模型变为如式(1)所示的Wiener过程。

1.2 随机冲击模型

现有考虑退化与冲击相互影响的研究中,大多采用泊松过程来描述系统受到的随机冲击过程。相比于其他随机过程,泊松过程主要有以下3点优势：① 泊松过程是一种点过程,用来表征随机冲击这种单事件效应现象是合理的；② 泊松过程具有无记忆属性,即冲击是随机发生的；③ 泊 松密度()可为任意形式,能较好地描述随机冲击的出现频次。

(3)

式中：()为随机冲击引起的退化水平变化量之和；为单次冲击引起的退化水平突变量。

上述随机冲击模型已经被应用于考虑随机冲击影响的退化建模中,而本文基于此模型,与非线性自适应Wiener过程相结合,建立退化模型,具体方法在1.3节展开介绍。

1.3 考虑随机冲击影响的退化模型

将随机冲击模型融入设备退化模型中,如图1 所示。图1中:横坐标为退化时间,纵坐标为退化量(),为失效阈值,(=1,2,3,…)为第次状态监测时间,为设备剩余寿命分布。当()达到预先设定的失效阈值时,即认定设备失效。

图1 考虑随机冲击影响时设备退化过程Fig.1 Degradation process of equipment with random shook

用随机过程{(),>0}表示设备退化过程,考虑随机冲击影响,设备实际退化过程分为2部分：设备正常退化和随机冲击对退化水平的影响。此时,设备退化模型为

(4)

2 随机冲击影响下设备剩余寿命预测

2.1 设备剩余寿命预测

在随机退化建模框架下,根据首达时间定义设备的寿命为

=inf{：()≥|(0)<)

(5)

式中：()为设备退化状态；(0)为退化初值;为设定的失效阈值。

由于Brownian运动的随机性,寿命为服从逆高斯分布的随机变量,根据文献[15]自适应Wiener过程寿命分布的概率密度函数(Probability Density Function, PDF)为

(6)

若～(,),、都为常数,为正实数,则以下期望公式成立：

(7)

根据引理1,可通过全概率公式,得到考虑随机冲击影响下自适应Wiener过程的寿命分布PDF为

(8)

根据随机过程首达时间的概念,当{(),≥0}首次达到失效阈值时,即认为设备失效。因此,基于监测数据0:={,,…,},将设备在时刻的剩余寿命定义为

=inf{：(+)≥|()<)

(9)

由式(8)可得,考虑随机冲击影响时自适应Wiener过程的剩余寿命分布PDF为

(10)

2.2 设备剩余寿命预测在线更新

为实现同时考虑时变不确定性、个体差异性、随机冲击影响等因素时设备RUL预测的在线更新,需构建退化过程的离散模型。此外,当受到随机冲击时,设备的退化速率可能会受到不同程度的影响,并且同样的冲击对不同退化速率设备的影响也不尽相同。因此,在自适应Wiener过程中,表征设备退化速率的漂移系数()和随机冲击影响的退化量可能会相互影响,且该影响随着设备运行的时间不断变化。鉴于此,本节将式(4)中的退化模型改写为

(11)

(12)

2) 状态估计

(13)

3) 协方差更新

(14)

若～(,)的元正态分布,且、为已知的常实数,、为已知的维常值向量,则

(15)

根据引理2,通过全概率公式可求得同时考虑随机冲击和个体差异性时,首达时间意义下设备RUL分布的PDF为

(16)

基于泊松分布的概率密度函数,产生个随机变量表示随机冲击次数,结果记为(),(),…,()，有

(17)

计算步骤3得到RUL分布的均值,令其作为随机冲击影响下设备剩余寿命的PDF

(18)

3 参数估计

(19)

(20)

E步：

(21)

M步：

(22)

根据状态空间模型,基于贝叶斯定理和条件概率公式,未知参数的对数联合似然函数为

(23)

基于构建的状态空间模型,以上条件期望值可通过RTS(Rauch-Tung-Striebel)平滑滤波计算。具体步骤如下：

1) 平滑滤波

(24)

2) 协方差初始值

(25)

3) 后向迭代

(26)

(27)

交替执行E步和M步,直至估计的结果满足预设的算法终止条件,即最大迭代次数或参数收敛。当有新的状态监测数据时,将之前估计得到的参数作为EM算法的初值,并基于新监测数据执行EM算法,以获得相关参数的最新估计值。

4 实例验证

4.1 数值仿真

由表1可见,相比于白灿等所提方法、不考虑随机冲击影响的方法,本文所提方法的参数估计结果更为准确,说明了本文所提方法的合理性和有效性,为后续准确预测剩余寿命奠定基础。

表1 模型参数估计值Table 1 Estimated values of model parameters

此外,设备退化过程首次超过预先设定的阈值即为首达时间,经过大量模拟仿真可得到寿命分布直方图如图2所示,并与白灿等所提方法得到的寿命分布PDF、本文所提方法得到的寿命分布PDF进行比较。为了后续便于表示,令本文所提模型记为M0,白灿等提出的随机冲击退化模型记为M1,不考虑随机冲击影响的退化模型记为M2。

图2 不同模型寿命分布PDFFig.2 PDF of life distribution with different models

图3～图5的结果反映了随机冲击对寿命分布的影响,即： ① 增大,加速设备退化,设备寿命随之不断下降,且寿命预测的不确定性增高； ② 增大,寿命分布的均值几乎不变,但寿命预测的不确定性增高； ③ 减小,设备寿命下降且预测不确定性也随之下降。

图3 不同均值uI下寿命分布PDFFig.3 PDF of life distribution with different mean values

图4 不同方差下寿命分布PDFFig.4 PDF of life distribution with different variances

图5 不同泊松密度μ下寿命分布PDFFig.5 PDF of life distribution with different Poisson densities

4.2 陀螺仪实例分析

陀螺仪是导弹等运载体惯性导航系统的核心设备,其性能的优劣直接决定了导航精度和任务的成败,同时具有结构复杂、精度要求高、价格昂贵等特点。通常,陀螺仪在安装到导弹后,会随弹经历长时间储存、测试、拆装、运输等不同环境过程,易受到随机振动冲击等影响,对表征陀螺仪精度性能的陀螺仪误差模型系数的退化过程产生加速退化影响,最终影响陀螺仪服役时间。因此,为了准确预测陀螺仪剩余寿命,为后续维护管理提供科学依据,必须考虑随机冲击对陀螺仪退化的影响,某机械转子陀螺仪在寿命周期中漂移系数退化数据如图6所示,其中的数据突变是由运输振动冲击导致的退化加速。

在工程实际中,该陀螺仪漂移系数超过预先设定阈值=0.37 (°)/h,认为其不再满足使用需求,即认为发生失效故障。为了验证所提模型RUL预测的有效性,用前67个数据估计模型参数,并用剩余数据对陀螺仪的RUL更新。利用上一节提出的EM算法估计未知参数。

图6 陀螺仪退化数据Fig.6 Gyro degradation data

图7 参数更新过程Fig.7 Update process of parameter

图8 参数P0|0更新过程Fig.8 Update process of parameter P0|0

图9 陀螺仪漂移数据拟合效果Fig.9 Gyro drift data fitting effect

在此,给出3种模型各自的RUL预测结果,具体如图10、图11所示。结果表明,相比于模型M2,模型M0考虑了设备运行过程中随机冲击对退化量和退化率的影响,更加符合实际退化情况,因此能显著提高预测准确度；与模型M1相比,模型M0基于非线性自适应Wiener过程进行退化建模,能克服测量间隔不均匀、测量频率不一致的影响,同时该模型将漂移系数视为Wiener过程,在未来RUL预测时漂移参数自适应变化。此外,MO模型考虑到同一批次设备不同个体间退化率的差异性的影响,并将随机冲击对退化率的影响融入状态空间模型中,实现RUL在线更新。因此,本文所提模型M0剩余寿命分布的PDF窄而尖锐,预测结果的不确定性较小,能获得更好的预测结果。

图10 不同模型RUL预测结果Fig.10 RUL prediction results of different models

图11 不同监测点RULFig.11 RUL prediction at different monitoring points

为了对不同模型预测结果进行量化比较,进而验证所提模型M0的有效性。定义设备的均方根误差(Root Mean Squared Error, RMSE)为

(28)

图12表明,随着设备不断退化,预测结果的RMSE也随之减小。所提模型M0对应RMSE最小,进一步证明了M0较M1、M2具有更好的模型拟合性,能有效降低预测结果的不确定性。

图12 不同模型RUL预测的均方根误差Fig.12 Root mean square error of RUL prediction with different models

为了定量分析预测结果,在此引入绝对误差(Absolute Error, AE)和评分函数(Scoring Function)进行性能评价。评分函数的计算公式为

(29)

由表2可知,3种模型RUL预测结果的绝对误差在不断减小,且模型M0对应的绝对误差最小,即预测结果更准确。评分函数主要用于描述预测值与真实值高估与低估的关系。一般评分函数分值越小,预测结果越准确。由式(29)可以计算得到,3种模型对应结果分别为0.124 1、4.816 9、25.097 0,模型M0对应的评分函数分值最小。接下来,采用美国国家航空航天局(NASA)公开的锂电池数据集进一步验证本文所提方法的有效性和通用性。

表2 不同监测点绝对误差Table 2 Absolute error at different monitoring points

4.3 锂电池实例分析

锂电池的容量退化数据随着充放电循环而不断降低,在退化过程中,由于充放电循环的中止,电池在没有电子压力的情况下,其容量数据可能会有所恢复。而这种容量的恢复可认为一种负向的随机冲击施加到退化过程中。本文从公开的锂电池数据#5、#6、#7、#18中,选择#5电池进行验证,锂电池退化数据如图13所示。在工程实际中,一般认为电池容量衰减30%将失效。

图13 锂电池退化数据Fig.13 Lithium battery degradation data

通过初始电池容量值减去每次充放电循环后的电池容量,得到符合Wiener过程的锂电池真实退化数据和预测退化数据,具体如图14所示。图14结果表明,本文方法预测数据与真实数据拟合效果较好,能较好地反映锂电池退化情况。

图14 锂电池退化数据拟合效果Fig.14 Lithium battery degradation data fitting effect

根据参数估计结果代入模型RUL预测解析式中,可得到模型M0和模型M1的RUL预测结果,具体如图15、图16所示。

由图15、图16可见,与模型M1相比,本文所提模型M0的预测效果更为准确。尤其在退化初期,模型M0能动态调整漂移参数,预测精度较高。此外,模型M0同时考虑时变不确定性、同一批次不同设备间的个体差异性以及随机冲击影响,考虑情况更加全面,预测结果也更为准确。

图15 不同模型RUL预测结果Fig.15 RUL prediction results of different models

图16 不同监测点RULFig.16 RUL prediction at different monitoring points

同样地,为了定性和定量分析本文所提方法的有效性,给模型M0、M1的RUL预测结果的均方根误差和绝对误差,分别如图17、表3所示。

图17 不同模型RUL预测的均方根误差Fig.17 Root mean square error of RUL prediction with different models

表3 不同监测点绝对误差Table 3 Absolute error at different monitoring points

图17表明,随着设备运行时间的推进,模型M0和M1预测结果的RMSE不断减小。在模型M0和M1中,所提模型M0对应RMSE较小,表明模型M0能有效降低预测结果的不确定性。

由表3可见,模型M0的绝对误差始终小于模型M1,其预测结果更为准确。综上,本文所提模型M0预测准确度高、不确定性小,并且具有一定的适用性,可为后续设备的维护决策提供理论依据。

5 结 论

针对考虑随机冲击影响的随机退化设备,本文提出一种基于非线性自适应Wiener过程的剩余寿命预测方法。

1) 利用非线性自适应Wiener过程和正态分布分别描述设备正常退化和随机冲击对设备退化量的影响,建立了能有效融合随机冲击影响的自适应Wiener过程随机退化模型。

2) 在首达时间意义下推导出剩余寿命解析式,构建状态空间模型实现RUL在线更新。通过仿真实验以及惯性导航系统陀螺仪、锂电池数据验证所提方法的有效性。

3) 本文所提方法基于非线性自适应Wiener过程,能够克服退化设备测量间隔分布不均匀、监测数据的测量频率与历史数据频率不一致的情况,并且考虑将来退化过程中自适应漂移的可变性以及随机冲击对设备退化的影响,可为存在随机冲击影响的设备剩余寿命预测提供理论依据,且具有一定的潜在应用价值。