基于G-AHSMM的设备剩余寿命预测研究

张青山， 张思岩， 肖 萌， 徐 伟

(沈阳工业大学 管理学院， 沈阳 110870)

在中国制造2025的背景下，以物联网、大数据、云计算为基础，制造业企业的生产车间逐渐向自动化、数字化、智能化转变。因此，设备的维修管理工作成为制造业企业迫切关注的问题，而制定合理的维修策略可以有效降低设备的维修成本、提高设备可用度，减少不必要的损失[1]。而作为设备故障预测与健康管理技术(PHM)的核心内容，设备剩余寿命(RUL)预测的研究热度只增不减。

一、设备寿命预测方法比较分析

设备剩余寿命预测在设备维修领域占有重要地位，其中，服役设备的寿命预测尤其重要，它不仅是保障企业生产系统正常运行的前提，也是企业查定计算生产能力、编制生产作业计划、制定设备维修和更新改造计划、作出设备投资预算等的重要基础。特别是面对生产制造单元内各道生产工序相互关联的复杂设备，如何精准预测设备剩余寿命，是学术界关注的热点，也是备受工业制造业企业关注而又亟待解决的现实难题。

对于RUL预测方法的研究，经历了一百多年的历程，大致可分为三个阶段：物理模型驱动方法、知识驱动方法和数据驱动方法。早在一百多年以前，基于物理模型驱动的预测方法就已初现雏形，随着研究的不断深入，模型逐渐扩展完善。如景所立等结合软件分析了电气化列车某部件疲劳裂纹的扩展特性，并对其扩展寿命进行预测[2]。由于专家知识系统不易建立，故基于知识驱动的预测方法在设备寿命方面应用较少，多用于设备故障的诊断或其他领域。基于数据驱动的方法是随着设备监测系统的发展而不断完善的，也是目前较为先进的预测方法。这类方法是通过分析处理监测系统收集到的设备运行数据而进行预测的，具体包括神经网络、支持向量机、深度学习、隐马尔可夫模型、贝叶斯理论等。表1展示了各个阶段常用的或有代表性的预测方法使用的情景条件及其优缺点，并对各个方法的应用文献进行举例。

表1 设备寿命预测方法对比

在解决设备剩余寿命预测的实际问题中，尤其是对于复杂设备来说，以上很多方法中的数学模型不能很好地描述设备退化过程，模型解释性较差，而隐马尔可夫模型(HMM)在具有强大理论知识的基础上，还具有丰富的数学结构[3]，从而对模型的解释更加具有说服性。对于HMM没有显示时间结构的这一缺点，可以通过由隐马尔可夫模型扩展到具有描述设备驻留时间的隐半马尔可夫模型(HSMM)解决。传统的HMM和HSMM要求各观测值之间相互独立，这一假设违背了现实性，而本文将设备历史监测数据引入观测概率矩阵的计算，克服了这一局限性。

二、隐马尔可夫系列模型与拓展

1. 隐马尔可夫系列模型

对于隐马尔可夫模型(HMM)的研究应用，学者们发表了大量文献进行探讨。如Fan等将显式状态驻留时间的分布参数引入HMM中，发现HSMM具有更大灵活性，且克服了HMM的局限性[13-14]。虽然HSMM比HMM更加贴合实际，但也存在不足之处。为了克服传统HSMM的局限性，学者们提出了很多基于传统HSMM进行改进的模型，如基于自适应的AHSMM[15]、埃尔朗分布和隐半马尔可夫相结合的寿命预测模型E-HSMM[16]、基于动态粒子群算法的PSO-HSMM[17]、基于EM算法的集成分段隐半马尔可夫模型即EM-SHSMM[18]等，均在不同方向、程度上对HSMM进行了改进。

2. 隐马尔可夫模型的拓展

综合以上文献可以看出，起初是通过HMM对设备进行预测性维护，但由于该模型中缺乏对于状态驻留时间的表述，进而提出HSMM用于改进设备状态的驻留时间问题。对于状态驻留时间的描述，学者们进行了大量研究，通过应用指数分布、高斯分布、埃尔朗分布等来表示这一过程，但都存在一定局限性。而有研究显示：当形状参数β取不同值时，伽马分布可转化成指数分布(β=1)和k阶埃尔朗分布(β=k且为正整数)，且伽马分布更加适合描述设备的退化过程。因此，基于前人的启发，本文提出采用伽马分布来估计设备状态驻留时间参数，使其具有更高的现实拟合性。而在前文提到的文献中，虽有将设备历史运行数据融入HSMM转移概率矩阵中的尝试，但没有考虑设备自相关观测量对观测概率矩阵参数估计的影响。本文在伽马(Gamma)分布的基础上，将这一因素考虑在内，计算设备的剩余有效寿命，使其预测结果更加精准。

三、设备剩余寿命预测模型构建

1. 构建思路

只有在明确设备退化过程的前提下，才能精准预测设备的剩余寿命，因此可将设备的整个生命周期大致分为几个退化状态。通过监测系统收集到的历史数据，对各个状态下的模型即模型库进行训练，用以明确设备运行到当下时刻所处的退化状态，进而根据对应状态下的模型，计算其状态转移概率和状态驻留时间及相关参数，最后得到该设备的剩余有效寿命。

2. G-AHSMM

HSMM是HMM的扩展形式，将一个描述状态驻留时间概率矩阵的参数考虑在内，能够更好地描述状态转移关系，使其处理实际问题具有更高的准确性[19]。设备剩余寿命预测模型(G-AHSMM)是以HSMM为基础，引入观测值的自相关观测量，再以伽马分布进行参数估计而建立的，其基本元素参数包括[20]：

(1) 隐藏状态集：S={s1，s2，…，sN}。N为隐藏状态数，其对应的状态序列Q={q1，q2，…，qT}，T为序列长度。

(2) 观测状态集：V={v1，v2，…，vM}。M为可观测的状态数，其对应的观测序列O={o1，o2，…，oT}，ot为t时刻所处的观测状态，且ot∈{v1，v2，…，vM}，t∈[1，T]。

(4) 状态转移概率分布矩阵：A={aij(d)}。aij(d)=P(qt+1=sj|qt=si，dt(i)=d)，i、j∈[1，N]，d∈[1，Di]。aij(d)表示设备从t时刻的隐藏状态si转移到状态sj的概率，且在t时刻停留了时间d；Di为设备在si状态驻留的最大时间。

(5) 观测值概率分布矩阵：B={bi(ot|ot-1)}。bi(ot|ot-1)=P(ot=vk|ot-1=vk-1，qt=si)，t∈[1，T]，k∈[1，M]。bi(ot|ot-1)表示设备在t-1时刻观测值基础上t时刻的观测条件概率。

(6) 状态驻留时间概率分布：D={pi(d)}。pi(d)=P(di=d|qt=si)，其中i∈[1，N]，d∈[1，Di]，表示设备当下状态为si且还要在此状态停留时间d的概率，d为设备在状态si下的驻留时间。

因此，一个标准的G-AHSMM需要由A、B、π、N、M、D几个主要元素组成，且它们之间存在一定的关系，即当A、B确定后，N与M便已知，故可将G-AHSMM记为λ=(π，A，B，D)，其具体结构如图1所示。

图1 G-AHSMM基本框架

标准HSMM假设系统的观测值是相互独立的，即在t时刻出现的观测量ot只与当前系统所处状态qt(qt∈S)有关，而与t时刻之前的观测量{o1，o2，…，ot-1}和状态{q1，q2，…，qt-1}无关，很显然，这并不符合实际。因此，本文将设备的自相关观测量整合到状态观测概率矩阵的参数估计中，建立基于自相关观测量的隐半马尔可夫模型，再以伽马分布进行参数估计，建立G-AHSMM用于设备寿命预测。

3. 模型基本问题

不论是HMM还是HSMM，亦或基于HSMM改进的G-AHSMM，均有三个基本问题需要解决[21]：

(1) 评估问题，即估计观测序列O的概率。求观测序列O在模型λ下出现的概率P(O|λ)时，常用前向后向(forward backward)算法[22]解决，相对其他两个问题较为简单。

(2) 解码问题，也称预测问题。在给定观测序列O的前提下，找出最可能出现的与之相对应的隐藏状态序列S。常用维特比(Viterbi)算法[23]解决该问题，复杂程度适中。

4. 算法改进

前向算法和后向算法均采用动态规划的原理，是解决隐半马尔可夫系列模型基本问题常用的算法。本文对陈永刚等提出的前向后向算法[25]进行了改进，得出以下中间变量。

已知G-AHSMM中λ=(π，A，B，D)，观测序列O={o1，o2，…，oT}，定义以下变量：

(1) 前向概率。定义设备在t时刻的隐藏状态为sj且在此状态驻留了时间d的情况下，观测到序列为o1，o2，…，ot的概率为前向概率αt(j，d)，则

αt(j，d)=P(o1，o2，…，ot，sj=qt，dt(j)=d|λ)=

(t∈[1，T-1])

(1)

前向概率的计算以t-d时刻隐藏状态下的前向概率αt-d(i)为基础，再与对应的状态转移概率αij(d)相乘求和，得到t时刻隐藏状态为sj的概率，在此基础上引入观测变量及其自相关观测量，即可得到t时刻观测到序列为o1，o2，…，ot且隐藏状态为sj的概率，其递推过程如图2所示。

图2 前向算法递推过程

进而可以得出对于给定模型λ，观测序列O={o1，o2，…，oT}出现的概率公式为

(2)

(2) 后向概率。定义设备在t时刻的隐藏状态为si且在此状态驻留了时间d，观测到序列为ot+1，ot+2，…，oT的概率为后向概率βt(i，d)，那么

βt(i，d)=P(o1，o2，…，ot，sj=qt，dt(i)=d|λ)=

(t=T-1，T-2，…，1)

(3)

后向概率的递推过程与前向概率恰好相反，以t+d时刻隐藏状态的后向概率βt+d(j)为基础，加入相对应的状态转移概率、观测概率等相关参数后，即可得到设备在t时刻的后向概率βt(i，d)，其递推过程如图3所示。

图3 后向算法逆向递推过程

(3) 定义在观测到序列O={o1，o2，…，oT}的前提下，设备在t时刻处于状态si，持续了时间d后转移到状态sj的概率为ξt(i，j，d)，在得到该状态下的前向概率、后向概率、状态转移概率、状态驻留时间概率、观测概率的前提下，可对其进行计算，表达式为

ξt(i，j，d)=P(qt=si，qt+d=sj，dt(i)=d|O，λ)=

(4)

(4) 定义设备在观测到o1，o2，…，oT的前提下，在t时刻处于状态si的概率为γt(i，d)，通过前向、后向概率计算可得，其表达式为

γt(i，d)=P(qt=si，dt(i)=d|O，λ)=

(5)

(5) 定义设备在状态si持续了时间d后观测到o1，o2，…，oT的概率为φt(i，d)，φt(i，d)在ξt(i，j，d)的基础上将设备在t时刻处于状态si这一条件作为前提，使得在计算过程中后向概率的代入有所差异，其表示形式为

φt(i，d)=P(dt(i)=d|qt=si，O，λ)=

(6)

5. 参数估计

(1) 基于自相关观测量的观测概率

(7)

bi(ot|ot-1)=P(ot=vk|ot-1=vk-1，qt=si)=

(8)

(2) 基于伽马分布的参数估计

假设设备的状态驻留时间服从伽马分布，即d～Γ(r)，其概率密度函数表达式为

(9)

式中：r为形状参数；θ为尺度参数；Γ(r)为伽马函数，满足均值为r/θ，方差为r/θ2。

若r>0且为整数，则伽马函数可定义为

(r-1)!

(10)

将式(10)代入式(9)得

(11)

综上，得到模型λ=(π，A，B，D)各参数估计值的表达式分别如式(12)～(16)所示：

(12)

(13)

(14)

(15)

(16)

(17)

(18)

6. 剩余寿命的时间构成

(19)

若T为设备的寿命周期，则ρ的表达式为

(20)

定义设备处于状态sk时，其剩余寿命为RULk，表达式为

RULk=ak，k(Dk+RULk+1)+ak，k+1RULk+1

(k∈[1，N-1])

(21)

当k=N时，表示设备到达生命周期的最后一个阶段，即故障发生前的状态，则此状态的剩余有效寿命只取决于当前状态的驻留时间，表达式为

RULN=aN，NDN

(22)

基于式(22)对式(21)不断迭代，即可计算出设备在各个状态时的剩余有效寿命。

四、G-AHSMM的模型求解

基于G-AHSMM的设备剩余有效寿命预测流程如下：

(1) 根据G-AHSMM设计分类器，形成模型库。若某一设备从正常运行到出现故障共经历m个退化状态，就需要建立m个分类器，可记为λ1，λ2，…，λm。

(2) 通过获得的数据对模型进行训练，运用式(12)～(16)对模型的参数不断修正，当找出使得观测序列出现概率P(O|λm)最大的各个参数时，则终止迭代，模型训练结束。

(3) 将当前要预测剩余寿命的设备观测序列代入各分类器中，计算观测序列条件概率P(O|λm)，必有一模型λi(i=1，2，…，m)使得P(O|λm)最大，那么该模型状态就是设备当前所处状态。

(4) 根据式(17)～(19)计算设备在各个健康状态下的持续驻留时间，再将得到的结果代入式(20)，计算设备有效剩余寿命。

模型求解流程如图4所示。

图4 G-AHSMM的求解流程

表2 不同健康状态之间的转移概率

五、例证分析

本文以美国宇航局艾姆斯预测数据存储库的涡轮风扇发动机为例，在G-AHSMM中将涡轮风扇发动机的整个寿命周期分为4个阶段，分别为正常状态N、退化状态D1、退化状态D2、故障状态F。共选取整个生命周期的60组历史监测数据进行分析，并将其划分为两个数据集——训练数据集(含30组数据)和测试数据集(含余下的30组数据)，分别用于G-AHSMM的训练和测试。在模型训练中，通过设置使得模型共迭代50次，收敛误差限定为e=0.000 001，得出G-AHSMM的参数训练曲线如图5所示。可以看出，在训练步数接近20的情况下达到限定误差，证明G-AHSMM的数据处理能力较好。

图5 G-AHSMM训练曲线

利用30组涡轮风扇发动机历史数据训练得到4个健康状态之间的转移概率，结果如表2所示。再通过式(17)、(18)计算发动机在不同状态下的均值和方差，结果如表3所示。

表3 不同状态下驻留时间的均值和方差

基于表2的数据和式(19)、(20)计算发动机在4个健康状态下的持续驻留时间，结果如表4所示。

表4 不同状态下的持续驻留时间

通过式(21)及表2、3中数据可计算发动机处于不同状态下的剩余有效寿命。假设发动机当前的健康状态为退化状态D2，则由表2～4可知，a33=0.927 7，a34=0.072 3，D3=22.760 7，RUL4=8.780 8。则RUL3=a33(D3+RUL4)+a34RUL4=29.895 9。而此健康状态下的实际剩余有效寿命为28，误差为6.767 8%(误差计算方式参考王宁等的研究[26])。运行结果说明，模型具有一定的可靠性。

将余下30组测试数据分别代入模型库中的4个G-AHSMM中，计算似然概率对数值。其中概率数值最大的G-AHSMM所对应的状态即为该测试数据形成的观测序列所处的健康状态，测试结果如表5所示。

表5 G-AHSMM 4种状态下测试数据运行结果

由表5可知，基于G-AHSMM的状态识别率为97.50%。由于用于测试模型的数据只有30组，故识别率没有十分接近100%，若测试组数更多，识别率定会高于97.50%。

再将这30组数据代入传统HSMM，测试结果如表6所示。对比两个模型的识别率可以发现，G-AHSMM具有更高的预测精度。HSMM训练过程与G-AHSMM相类似，此处不再赘述。

表6 传统HSMM运行结果

每个健康状态下选取2组样本数据分析G-AHSMM预测RUL的精准度，并对比传统HSMM的误差率，结果如表7所示。表7数据说明，基于设备状态监测数据驱动的剩余寿命预测模型的预测结果更加准确。

表7 寿命预测结果对比

六、结 论

本文提出了基于G-AHSMM的设备剩余寿命预测方法。该方法是在对已有设备寿命预测方法比较分析的基础上，将基于自相关观测量的隐半马尔可夫模型加以拓展，并与伽马分布相结合，构建了设备状态监测数据驱动的剩余寿命预测模型G-AHSMM。该模型规避了以往“状态观测值之间相互独立”的不实假设，并通过例证分析证实，其比传统HSMM具有更高的现实拟合性、模型可靠性和预测精准性，这也是大数据驱动下设备寿命预测研究的新进展。

本文提出的G-AHSMM是在设备运行数据完备的情况下对其剩余有效寿命的预测，没有考虑数据缺失问题。因此，下一步的研究方向是通过相关算法对缺失的数据进行填补，在降低数据分析难度的同时，使得预测结果更加精确。