装备剩余寿命预测平行仿真中模型动态演化方法*

葛承垄，朱元昌，邸彦强，崔浩浩

(陆军工程大学石家庄校区， 河北 石家庄 050003)

装备平行仿真首次完整提出是在2016年亚仿/秋季仿真大会上[1-2]。装备平行仿真旨在将实际装备和仿真系统结合在一起，仿真系统利用传感器采集的装备信息演化仿真模型，实际装备受益于仿真系统的仿真结果，从而提高装备的运用效能。以这种模式运行的仿真系统称之为平行仿真系统。

装备平行仿真中的一个重要概念是实时数据驱动的模型演化，即在装备实时数据驱动下，仿真模型形态或参数进行适应性调整的过程，包括模型适宜性选择和模型参数演化两个方面的内涵，这被认为是其区别于以往仿真技术的典型特征。以往仿真技术中，仿真模型侧重于一次性构建，仿真运行后模型形态和模型参数不再发生变化，并不存在模型演化过程。在装备平行仿真中，模型演化的目的是为了满足动态变化的仿真需求，通过基于实时装备数据的动态仿真，动态提高仿真逼真度和仿真预测准确度，为仿真预测和辅助决策提供数据支持。仿真模型及其演化属于装备平行仿真的模型理论范畴，也是装备平行仿真研究中的基础问题。装备平行仿真还具有虚实共生、平行运行、数据驱动等技术特征。目前，装备平行仿真主要包括两个研究分支：一是面向战场指挥决策领域的平行仿真，代表性的研究学者包括邱晓刚[3]、周芳[4]、窦林涛[5]等，此类平行仿真已初步建立演化建模框架，包括模型动态匹配、模型类别修正、模型参数校准和实体模型行为意图更新等，相关演化方法的研究正在逐步展开；二是装备维修保障领域中面向装备剩余寿命(Remaining Useful Life，RUL)预测的平行仿真，已完成模型构建机理研究和模型演化方法设计。本文关注的是面向装备RUL预测的平行仿真中模型演化方法研究。

在装备维修保障领域中，由于外部冲击、磨损、疲劳、腐蚀等原因，装备的性能将不可避免地发生退化，最终引起装备故障甚至造成严重事故，因此需要视情对装备进行健康维护。近年来，预测与健康管理技术是实现装备健康维护的主要方法[6]，它通过预测装备的RUL，对装备进行合理的维修与管理，保证设备运行的安全性、可靠性与经济性，其关键内容包括RUL预测与健康管理两个方面。其中，RUL预测是预测与健康管理技术的核心内容[7]。RUL预测主要是指解算RUL分布，即概率密度函数(Probability Density Function，PDF)，它表征了寿命预测的不确定性，是维修决策的重要依据。

当前，RUL预测方法可分为失效物理模型法、基于统计的方法和人工智能法[8]。对于复杂装备来说，其失效机理很难获得，因此后两种方法得到更多关注。基于统计的方法和人工智能法分别通过统计模型、机器学习进行数据拟合，从而预测RUL。然而基于统计的方法和人工智能法仍存在以下典型不足：第一，在预测过程中通常假定装备退化模式在整个预测周期中是固定不变的，利用单个模型进行预测，然而退化模式却可能呈现多种退化模式混合的情况，其中连续退化与未知离散冲击混合是一种典型情形；第二，模型参数往往通过极大似然估计等离线方法获得，不能随着退化数据的积累而在线更新；第三，难以得到RUL概率密度函数的解析表达式，制约预测方法的实时性，尤其是人工智能法，它还需要大量训练样本，这在实际中往往无法满足。因此，针对现实中普遍存在的连续退化与未知离散冲击混合的退化过程，迫切需要在RUL预测中同时考虑模型适宜性选择和参数在线演化。依据其技术原理和典型特征，装备平行仿真为解决此类RUL预测问题提供了新思路。本文研究涉及面向混合退化装备RUL预测的平行仿真中仿真模型构建、模型动态演化以及基于平行仿真的RUL实时预测等问题。

1 平行仿真建模

1.1 建模分析

文献[1,9]指出，构建装备性能退化状态空间模型(State Space Model，SSM)是面向装备RUL预测平行仿真的建模方向。退化状态估计是RUL预测的前提，装备性能退化SSM兼顾了退化过程的动态性和时变性，易于进行退化状态估计，有利于RUL预测。特别地，考虑到退化状态是经平行仿真得到的，故可称之为仿真退化状态。

由于基于统计的方法更容易得到RUL概率密度函数的解析表达式，本文将随机过程方法与SSM结合起来构建平行仿真模型。随机过程适宜描述退化过程的随机性和不确定性，Wiener过程和Gamma过程是最常用的随机过程，但Gamma过程适合于描述具有严格单调特征的退化过程，适用条件过于苛刻，而Wiener过程具有适用范围广、首达时(First Hitting Time，FHT)分布明确等建模优势。因此，在平行仿真建模中将其与SSM建模法结合起来，构建Wiener状态空间模型(Wiener SSM，WSSM)是一种合理选择[9]。特别地，为描述带未知离散冲击的混合退化过程，宜建立多态WSSM作为平行仿真模型，它包括连续退化模型和含未知离散冲击的退化模型两种形态。

1.2 多态WSSM构建

多态WSSM构建涉及装备混合退化状态方程和观测方程。装备混合退化状态方程包括两种形态，一种是连续退化状态方程，另一种是带未知离散冲击的退化状态方程。首先，利用Wiener过程构建连续退化状态方程，有

x(t)=x(0)+ηt+σB(t)

(1)

式中：{x(t),t≥0}是由标准Brownian运动{B(t),t≥0}驱动的连续退化过程，且有B(t)～N(0,t)；x(0)为初始退化状态；η、σ分别是标准Brownian运动的漂移系数和扩散系数。对式(1)进行Euler离散化，得到在离散时间点tk(k=1,2,…)上的不考虑离散冲击的状态方程

(2)

随机离散冲击的未知特性是指冲击时刻未知，其对装备性能造成的损伤可用一个冲击变量表示。根据冲击的离散特性，将冲击造成的损伤融入式(2)中，得到带未知离散冲击的退化状态方程，有

(3)

式中，D为冲击造成的损伤。假设冲击的到达服从Poisson过程[10]，即有

(4)

式中，ρ为冲击到达率，M(tk)表示从初始时刻至k时刻冲击出现的总数量，n的取值为0或1[11]。为简便起见，这里假定Poisson过程与Brownian运动相互独立。

装备退化观测数据y(t)与x(t)的随机关系可由观测方程描述，即

y(t)=x(t)+π(t)

(5)

式中，π(t)～N(0,φ2)，φ2为测量噪声的方差，并假设π(t)与Brownian运动B(t)相互独立。对式(5)进行Euler离散化，得到在离散时间点tk(k=1,2,…)上的观测方程为

yk=xk+φζk

(6)

根据式(1)～(6)，可以得到在离散时间点tk(k=1,2,…)上的多态WSSM为

(7)

2 多态WSSM演化

2.1 模型演化分析

在SSM建模方式下，仿真退化状态可通过实际退化观测值驱动多态WSSM运行而获得，仿真退化数据与实际退化观测值构成完全退化数据。在SSM建模框架下，仿真退化状态的实质是一种隐含退化状态。为实现多态WSSM动态演化，提出如下演化机制：一是基于交互多模型(Interactive Multiple Model，IMM)滤波[12]的模型软切换，即利用观测数据修正仿真预测结果，并动态计算不同仿真模型的概率，实现仿真输出和实际观测数据的同化，得到仿真退化状态估计；二是在实时退化观测数据驱动下在线估计模型参数。模型软切换和参数在线估计不是孤立执行的，而是相互迭代，通过二者的不断迭代，动态校正模型预测输出，从而动态提高仿真预测准确度。

在多态WSSM中，冲击的未知特性使得无法通过观测数据得知是否有冲击到达，这就导致无法得知仿真模型切换的判断条件，因此需要利用IMM滤波计算不同模型形态的概率。由于装备性能退化过程受离散冲击影响，具有退化轨迹突变、非平稳的特点，考虑在IMM滤波中采用强跟踪滤波器[13]。基于IMM强跟踪滤波(IMM Strong Tracking Filter，IMMSTF)的模型软切换包括5个阶段，即模型输入交互、模型预测、滤波、模型概率计算和模型输出交互。多态WSSM是一种含有隐含状态的状态空间模型，故考虑利用期望最大化(Expectation Maximum，EM)算法[14]实现模型参数演化，EM算法能有效估计此类SSM的模型参数。综合以上分析，装备平行仿真中模型动态演化示意如图1所示。

图1 装备平行仿真中模型动态演化Fig.1 Model dynamic evolution of equipment parallel simulation

2.2 基于IMMSTF的模型软切换

(8)

(9)

(10)

(11)

(12)

(13)

(14)

则仿真模型u的概率为

(15)

(16)

(17)

2.3 模型参数演化

在多态WSSM中，模型参数θ包括μ0、Σ0、η、σ、D和φ，其中μ0、Σ0表示初始退化状态x0的均值和协方差。根据EM算法，在监测时刻k、EM算法第j步时，模型参数θ={μ0,Σ0,η,σ,D,φ}的估计可由式(18)获得。

(18)

(19)

2.3.1 计算联合对数似然函数的数学期望

(20)

根据多态WSSM可知

p(yi|xi,θ)=N(yi;xi,φ2)

p(x0|θ)=N(x0;μ0,Σ0)

(21)

(22)

IMMBS算法包括后向滤波和模型状态融合两个阶段。后向滤波是指从最新的测量值开始后向滤波，其运算流程与IMMSTF类似，但是二者也存在明显差异。IMMSTF先执行输入交互再执行一步预测，而后向滤波则是先执行一步预测再执行输入交互。后向滤波可以分为后向一步预测、后向输入交互、后向滤波更新、后向模型概率计算和后向输出融合5个步骤。特别地，后向滤波的后向一步预测方程为：

(23)

(24)

(25)

(26)

(27)

(28)

2.3.2 最大化联合对数似然函数的数学期望

时刻k处、期望最大化算法第j步时，θ的估计可利用对联合对数似然函数的数学期望取偏导数并令偏导数为0求得，即

(29)

求解可得参数θ的在线估计值，即

(30)

3 剩余寿命分布实时解算

根据FHT的概念，装备RULT定义为x(t)首次通过失效阈值w的时间[17]，即

T(w)=inf{t:x(t)≥w|x(0)

(31)

鉴于多态维纳状态空间模型在某特定时刻可能存在式(2)、式(3)两种仿真退化状态方程，无法直接得到寿命分布的PDF，需对x(t)进行转化。根据文献[11]可知，维纳过程可改写为

x(t)=xk+η(t-tk)+σB(t-tk)

(32)

其中，B代表布朗运动。将式(2)、式(3)进行融合，即将n次Poisson冲击造成的损伤融入式(32)中，有

x(t)=xk+nD+η(t-tk)+σB(t-tk)

(33)

根据式(31)中RUL的定义和维纳过程的FHT分布性质可知，多态WSSM中以n次Poisson冲击、xk、θ和Yk为条件的RUL概率密度函数服从逆高斯分布，即

(34)

其中，Tk为监测时刻k处的剩余寿命，其概率密度函数可以写为

f(Tk|n,xk,θ,Yk)

(35)

引理1 设Ω～N(γ,ξ2)，A、B、C为常数，则有

(36)

引理的证明见文献[18]。

定理1对于混合退化过程{x(t),t≥0}，时刻k处装备剩余寿命概率密度函数为

(37)

fT(Tk|θ,Yk)

□

根据数学期望的定义，平行仿真系统利用数值积分计算RUL的数学期望值，即

(38)

根据k时刻RUL概率密度函数和数学期望的解析表达式，平行仿真系统能够以在线、实时的方式计算RUL的概率密度函数和期望值，为装备维修决策提供数据支撑。

4 实例研究

4.1 数据介绍

机械装备中轴承的性能退化是典型的含有冲击特性的混合退化过程，本文采用某轴承全寿命试验数据进行模型演化方法验证。该数据由法国FEMTO-ST研究所提供[19]，全寿命试验在PRONOSTIA平台上进行，近年来被广泛应用于可靠性领域的方法验证。试验过程工况条件以及试验数据相关信息详见文献[11]，本文仍以轴承1_3为例进行实例研究。试验过程中采集的是轴承1_3的振动加速度信号，该信号的均方根(Root Mean Square，RMS)值是常用的退化特征量，计算公式为

(39)

式中，N是采样点数，这里N=2560，xi为第i个采样点对应的振动加速度信号。轴承1_3的RMS如图2所示。根据图2可知，轴承1_3的性能退化过程冲击特性明显，适宜用于验证本文方法。根据文献[11]可知，失效阈值设为4.714 5，即w=4.714 5。

图2 轴承1_3的均方根值Fig.2 RMS of bearing 1_3

4.2 多态WSSM动态演化与RUL实时预测

多态WSSM参数初始设置为x0=0.2、η=0.02、σ=0.5、D=0.02、ρ=0.5、τ=1、φ=0.1，Markov模型转移概率矩阵P=[0.50.5;0.60.4]，模型初始概率分别为0.6和0.4，即μ0=[0.60.4]T。得到的退化轨迹对比如图3所示，仿真退化轨迹和实际退化轨迹差异较小，表明在实时退化数据驱动下，仿真退化轨迹能有效逼近实际退化轨迹。利用均方根误差(Root Mean Square Error，RMSE)来量化退化轨迹对比结果，其计算公式为

(40)

其中，r=842为监测时间点数目。经计算，两种退化轨迹的RMSE仅为3.497%，充分说明本文提出的模型动态演化方法能有效实现轴承1_3性能退化过程的建模与仿真。

图3 退化轨迹对比Fig.3 Comparison of the degraded traces

模型概率如图4所示。由图4可知，在t1500至t1765监测时间段内，由于Poisson冲击特性并不显著，此时线性退化特性较为明显，线性退化模型(模型1)的概率明显高于带未知离散Poisson冲击退化模型(模型2)的概率，线性退化模型的概率保持在0.74左右，带未知离散Poisson冲击退化模型的概率保持在0.26左右，此时处于主导地位的模型是线性退化模型。但随着时间的推移，Poisson冲击特性越发显著，突出表现在t1766、t1827、t1877、t2130、t2234等时刻，带未知离散Poisson冲击退化模型的概率总体呈现动态上升的趋势，反之，线性退化模型的概率呈现动态下降的趋势。在退化后期，带未知离散Poisson冲击退化模型的概率超过线性退化模型的概率，说明带离散Poisson冲击退化模型更适宜描述当前的退化过程。从以上分析可以看出，本文方法能有效实现仿真模型的“软”切换，使之满足寿命预测对模型适宜性的需求。值得注意的是，由于模型库中仅存在两种模型，在模型软切换条件下，二者模型概率之和为1，因此两个模型概率曲线关于概率μ=0.5对称。

图4 模型概率Fig.4 Model probability

在轴承退化数据的动态驱动下，模型参数θ的在线估计结果如图5所示。由图5可知，漂移系数η在0.004附近动态波动，波动区间为[0,0.012]，尤其是在t1766、t1827、t1877等冲击特性明显的时刻，η的数值波动较大，反映出轴承1_3的退化速率加快；扩散系数σ能较快收敛，在冲击特性明显的时刻，扩散系数发生较大波动，并达到新的收敛状态，扩散系数的收敛有利于获得稳定的剩余寿命概率密度函数，同时扩散系数的收敛值较小，有利于获得更为狭窄的剩余寿命概率密度函数，提高剩余寿命预测的准确性；Poisson冲击造成的损伤D在初始时刻附近波动较大，后在区间[0.03,0.07]内动态波动，将损伤D考虑到剩余寿命预测中，能有效减少甚至避免“欠维修”的发生；测量误差的标准差φ收敛速度较快，收敛值约为0.01，反映出测量误差的波动逐渐稳定。

(a) 参数η在线演化(a) Online evolution of parameter η

(b) 参数σ在线演化(b) Online evolution of parameter σ

(c) 参数D在线演化(c) Online evolution of parameter D

(d) 参数φ在线演化(d) Online evolution of parameter φ图5 多态WSSM参数在线演化Fig.5 Online evolution of polymorphic WSSM parameters

在每一个监测点，随着模型演化的执行，平行仿真系统利用式(37)和式(38)预测轴承1_3的剩余寿命，包括剩余寿命概率密度函数及其数学期望。以100个监测时间点长度为预测间隔，可得到轴承1_3在t1600，t1700，…，t2300处共计8个监测点的剩余寿命预测结果，如图6所示。

图6 轴承1_3在不同监测点处的剩余寿命预测结果Fig.6 RUL prognostic results of bearing 1_3 at different monitoring time

根据图6可知，在t1600，t1700，…，t2300等监测时间点，RUL概率密度曲线均能有效覆盖实际RUL，并且随着轴承1_3退化数据的动态注入，RUL概率密度函数逐渐收窄，右偏特性减弱，正态特性越发明显，说明通过模型软切换和参数在线估计，模型匹配度逐渐提高，RUL预测的不确定性逐渐减小。此外，RUL期望值与真实RUL之间的误差较小，RUL期望值接近PDF峰值对应的RUL，表明PDF的不确定性较小，有利于辅助维修方案的制订。

4.3 比较研究

图7 不同方法剩余寿命预测对比结果Fig.7 Comparative RUL prediction results of different methods

图8 第1900个监测时间点处对比结果Fig.8 Comparative results at the 1900th monitoring time

为量化预测结果，给出两个评价指标：平均相对精度(Mean Relative Accuracy, MRA)和总均方差(Total Mean Square Error, TMSE)，它们的定义分别为：

(41)

(42)

经计算得到的方法评价结果如表1所示。根据表1可知，较之不考虑模型软切换的演化方法，本文方法的MRA明显较大，说明经过模型动态演化后剩余寿命预测的相对预测精度更高；本文方法的TMSE明显较小，说明剩余寿命预测的误差更小。

表1 方法评价结果

5 结论

以连续退化与未知离散冲击混合的混合退化装备剩余寿命预测为背景，提出一种装备平行仿真中模型演化方法。该方法构建了多态Wiener状态空间模型，利用交互多模型强跟踪滤波动态计算模型概率，实现模型软切换，利用期望最大化算法在线更新模型参数。通过模型软切换和参数在线估计的迭代，动态提高了仿真逼真度，得到了首达时意义下的剩余寿命分布解析表达式。通过模型动态演化，有效提高了剩余寿命预测准确度，减小了预测不确定性。该方法还具有可扩展性，可为解决更为复杂条件下平行仿真在寿命预测领域中的应用问题提供模型演化框架。