基于HDP-CHMM的机械设备性能退化评估

王 恒， 季 云， 朱龙彪， 刘 肖

(南通大学机械工程学院 南通，226019)

引 言

为了保证机械设备长期安全运行，工业现场需要对机械设备进行状态监测和性能评估，实现预知维护和管理。隐马尔可夫模型为一双重随机过程：一个是Markov链，用来描述隐藏状态间的转移概率，此为基本随机过程；另一个是描述每个隐藏状态下产生观测值的一般随机过程。HMM这一双重随机特性，能够很好地描述机械设备运行过程中观测到的退化征兆信号(如振动、转速和位移等)与隐藏的衰退状态之间的随机关系，在机械设备性能退化评估与预测中得到广泛应用[1-2]。

传统HMM的定义和估计过程存在严重不足，例如，HMM参数训练学习中需要预先确定模型的数学结构，不能解决在模型训练过程中出现的过适应或欠适应问题。刘新民等[3]建立了基于支持向量机(support vector machine, 简称SVM)与HMM串联结构的故障诊断模型，实现了对直升机减速箱的状态识别和诊断。张继军等[4]利用HMM自身的状态识别和转移回溯能力，结合多智能体遗传算法(multi-agent genetic algorithm, 简称MAGA)实现了温控放大器的状态分类。上述文献中HMM的初始状态数都是依据传统经验设定的，缺乏科学性和通用性。曾庆虎[5]提出了基于最小描述长度准则(minimum description length, 简称MDL)学习算法，通过MDL优化调整状态数，但模型状态数还是需要预先确定，计算过程较为复杂。腾红智等[6]基于CHMM对齿轮箱全寿命过程的退化状态识别进行了研究，提出了基于K均值算法和交叉验证相结合的状态数优化方法，但K均值聚类算法仍需要预先确定状态数，且交叉验证对不同的退化状态数都要进行训练并检验分类错误率，计算时间长，效率低。张星辉等[7]建立了一组聚类方法评价指标，利用K均值聚类算法对状态特征进行聚类，通过指标评定结果从中选取模型的最优状态数。但是设备运行过程中，随着监测数据的更新，退化状态数也需要随之不断更新，关于如何有效地确定最优退化状态数还需要进一步深入研究。

针对HMM参数的优化问题，提出了一种基于分层狄利克雷过程(hierarchical Dirichlet process, 简称HDP)和HMM相结合的机械零件性能退化评估方法，将HDP引入HMM中，结合HMM良好的分析和建模能力，实现运行过程中的退化状态识别。轴承退化性能评估结果表明了该方法的有效性和适用性。

1 HDP模型

1.1 DP定义

狄利克雷过程定义为关于一组分布或者随机测度的分布，假设参数服从一类样本空间上的宽先验分布，参数的后验分布通过采样推断，该狄利克雷模型及其扩展模型则具有良好的聚类特性。近几年，狄利克雷模型已经在机器学习、生物信息学、文本聚类和图像分割等方面有较好的应用。Ferguson首次提出狄利克雷过程的定义，G0为测度空间Θ上的随机概率测度，参数α为正实数。对于测度空间Θ的任意有限划分A1,,Ar，如果存在如下关系

(G(A1),，G(Ar))～Dir(αG0(A1)，,αG0(Ar))

(1)

则G服从由基分布G0和参数α组成的DP，即

G～DP(α,G0)

(2)

1.2 HDP过程

狄利克雷过程可以实现一组数据的聚类，但在研究多组数据的聚类问题时，单纯利用狄利克雷混合模型无法进行建模分析，这大大限制了其应用。针对这个问题，引入分层狄利克雷过程[8-9]。分层狄利克雷过程将基础分布G0扩展为一个服从狄利克雷过程的随机概率测度G，即G～DP(γ,H)，多个数据源Gj～DP(α,G)将共享基础分布G的离散原子，实现多组关联数据的聚类。

HDP无法实现对其过程的采样，在实际应用中往往采用不同形式的构造实现HDP过程的应用。HDP的Stick-breaking构造可以分为两层进行[10]。第1层Dirichlet过程为

(3)

其中：H为基础分布；β为第1层分布的权重；k为任意初始聚类数；φk为抽样的原子。

第2层每组Gj的Stick-breaking构造方式为

(4)

其中：Gj为扩展分布；aj为第2层分布中每组的权重。

在HDP模型中，狄利克雷过程作为参数的先验分布存在，假设模型中的观测数据xji表示第j组第i个观测数据，其分布形式定义如式(5)，两层狄利克雷过程分别用上述的Stick-breaking构造[10]。

(5)

HDP模型参数的含义如下：xji为观测数据；k为类别标签的初始聚类数；γ，α为HDP模型的超参数；β，aj为HDP模型每层的权重系数；λ为基础分布的采样参数；φk为随机抽样的离散原子。

2 基于HDP-CHMM的设备性能退化评估算法

2.1 基于HDP的退化隐状态数确定

隐状态数K的确定是HMM模型训练和测试的关键，但目前HMM隐状态数大多根据经验人为设定，或者通过训练一个最可能的隐状态数来达到近似的目的，这些方法很难将各种类别都考虑到，且新的数据中也可能有未知类型出现，依靠训练样本得到的固定模型结构对新的观测数据的适用性和涵盖性不强。HDP算法不依赖于训练样本，而且随着数据的变化，模型结构能够实现自适应调整，实现动态聚类。基于HDP的隐状态数确定算法流程图如图1所示，主要步骤如下。

图1 基于HDP模型的隐状态数确定Fig.1 The number of hidden states of HDP

1) 初始化HDP模型的参数，任意给定隐状态数K，迭代次数M，超参数α，γ。

2) 假设输入的多组观测数据xji服从多项式分布，分布密度函数为f(·|θ)；观测数据分布的参数θji服从其共轭分布DP分布，分布密度函数为h(·)；HDP的权重系数β和aj分别通过式(3)和式(4)中的Stick-breaking构造先验分布，观测数据xji的条件分布为

(6)

其中：fkxji(xji)为观测数据x中除了xji外的条件概率分布。

3) 基于马尔科夫链蒙特卡罗(Markov chain Monte Carlo, 简称MCMC)统计模拟方法简化条件概率的计算，并实现模型后验参数的更新，利用Augmented representation后验采样算法[11]实现对参数β的更新(ajk类似)，根据超参数γ和α对β进行采样

(7)

其中：m.1,,m.k利用直接分配后验采样算法[12]进行采样。

(8)

其中：mjk为第j组的聚类数目；Mjk为除了第j组的第k类之外的聚类数；nj.k为第j组观测数据中属于k类的数据总数；s(n,m)为Stirling数；γ为权重的第K+1个值。

通过对β和m的交叉采样实现参数更新，并在Stick-breaking构造中运用直接分配后验采样算法对观测数据的指示因子Zji进行采样，每次只更新一个数据的聚类属性，当某个类簇中元素个数为0时，K减1，继续迭代，待聚类数目稳定时停止迭代，获得最终的状态数目K。

4) 超参数γ和α决定DP过程的离散程度，初始值的选取对聚类结果影响很大，大大影响了HDP模型的通用性[13]。因此，在进行后验参数的更新时，同时对超参数进行更新。设α～Γ(a,b)，对任意k=1,2，，n，P(k丨α)服从Beta分布，则

(9)

其中：β′为Beta分布。

从β′(α+1,n)中抽样参数η，则参数α和η的联合分布为

p(α,η|k)∝p(α)αk-1(α+n)ηα(1-η)n-1

(10)

参数α的后验分布为

p(α|η,k)∝αa+k-2(α+n)e-α(b-log(η))

(11)

则

(α|η,k)～aηG(a+k,b-log(η))+

(1-aη)G(a+k-1,b-log(η))

(12)

其中：aη为权重。

(α|η，k)服从两个Gamma分布之和，且分布系数分别为aη和1-aη，超参数α和γ分别从后验分布aη/(1-aη)中不断迭代采样实现更新。

2.2 基于CHMM的机械设备退化状态评估

基于CHMM的机械设备退化状态评估算法的流程图如图2所示，主要步骤如下。

图2 基于CHMM的设备退化状态识别Fig.2 Identification of equipment degradation status based on CHMM

1) 根据步骤1～4确定的状态数K作为CHMM模型的输入参数，在一定约束条件下，随机初始化模型的其他参数(初始状态Q，转移矩阵A，观测值矩阵B，混合高斯数N)。

2) 运用CHMM模型进行建模，将观测数据作为模型的输入，通过Baum-Welch算法对观察数据进行训练，利用EM算法求概率参数模型的最大似然估计，重估权值w、均值μ和方差ξ，通过不断迭代估计模型参数λ=(Q,A,B)。

3) 根据训练得到的CHMM模型，利用Viterbi算法计算已知模型参数λ最可能的隐藏状态序列，即P(O|λ)的最大值，通过路径回溯求得每个观察序列最可能的状态，即得到退化状态转移曲线。

3 应 用

应用研究采用美国USFI/UCR的智能维护中心提供的轴承全寿命数据[14]，图3为实验装置。4个ZA-2115双列轴承并列安装在同一轴上，由恒定转速2 kr/min的直流电机驱动，在轴承和轴上加载约26 671 N的弹性径向载荷，采用加速度传感器采集振动信号，采样频率为20 kHz。轴承1在连续运转约163 h外圈出现严重损伤，共采集982组数据，采用轴承1的全寿命数据进行建模与分析。

图3 轴承加速寿命实验装置示意图Fig.3 Schematic diagram of the bearing life device

DP模型的扩展模型HDP利用分层构造DP原理，能够实现多组特征值的聚类。笔者基于HDP模型对轴承偏度指标(skewness)和峭度指标(kurtosis)进行训练，以获取CHMM模型中的最佳退化状态数。以轴承峭度指标为例，分析HDP模型中超参数α的不同初始值对聚类结果的影响，由图4可知，当α分别取2，20，2 000时，该模型聚类结果均能收敛到相同的值，可见，HDP模型的超参数初始值选取对聚类结果不敏感。

图4 不同α值影响分析Fig.4 Effect analysis of different valuesα

选取峭度指标和偏度指标作为HDP模型的观测数据。假设观测数据服从Multinomial分布，观测数据分布的参数服从DP分布。设定初始聚类数目K=50，聚集参数α=20，γ=1，迭代次数M=300，通过Stick-breaking构造HDP，利用不断采样更新参数，实现基于HDP的自动聚类，结果如图5所示，聚类结果趋近于5，且聚类结果比图4中单一特征量输入更稳定。因此，将K=5作为CHMM模型隐状态数进行模型的训练和测试。

利用混合高斯模型来拟合各状态下的观测值概率密度函数，在构造连续隐马尔科夫模型时，高斯分量数目N取3，基于Baum-Welch算法重估参数(Q,A,B)，获得经过多次迭代的重估模型。利用Viterbi算法计算P(O丨λ)的最大值，即得到退化状态转移曲线。

图5 基于HDP模型的多组观测数据混合聚类Fig.5 Multi group observation data mixed clustering of HDP

基于HDP-CHMM获得的轴承退化状态转移曲线如图6所示。轴承从正常状态到失效状态的全寿命历程中一共出现了5次不同状态，分别为正常状态、早期退化状态1、中度退化状态2、严重退化状态3和失效状态。通过该模型可以找出轴承运行时的早期故障点在576号文件处，严重故障点在930号文件处，并能识别轴承在运行过程中的一系列退化状态，为轴承的早期维护和维修提供理论指导。

图6 基于HDP-CHMM的轴承退化状态识别Fig.6 Identification of bearing degradation status based on HDP-CHMM

为了验证HDP-CHMM算法的有效性，与文献[15]中基于K-S检验的轴承退化状态识别方法进行了对比。由图6，7可知，轴承1全寿命数据在基于K-S检验的退化评估中，早期故障点在第533序列处，严重故障点在962处，与基于HDP-CHMM算法识别的早期故障和严重故障点大体一致。在轴承性能退化阶段，由于振动波动异常剧烈，基于距离的K-S检验法将轴承振幅跳变剧烈处的状态划分过多，在工作过程中不利于轴承实际状态的判定。

图7 基于K-S检验的轴承退化状态识别Fig.7 Identification of bearing degradation status based on K-S test

4 结 论

1) 提出了一种基于HDP-CHMM的机械设备性能退化评估方法，利用HDP的分层聚类特性，将轴承多个特征指标作为模型的输入，使聚类结果更加准确，解决了传统的CHMM状态数必须预先设定的不足，实现了模型结构的自适应调整。

2) 滚动轴承全寿命数据建模结果表明，HDP-CHMM可以有效地找出轴承运行时的早期故障点和严重故障点，便于设置预警机制，并能够识别轴承在运行过程中的一系列退化状态，为基于状态的设备退化评估提供了一种新的方法。

3) HDP-CHMM算法和K-S检验算法对比可知，HDP-CHMM算法可以对轴承实际运行过程中的状态转移进行建模，有效识别出轴承运行中的不同退化状态，比K-S检验算法基于距离的诊断方法更符合设备实际退化过程。