基于隐半马尔科夫模型的滚动轴承故障预测

1.隐半马尔科夫模型

隐半马尔科夫过程（Hidden Semi-Markov Model，HSMM）是在马尔科夫（Markov）的基础上发展起来的，马尔科夫链是状态和时间参数都离散的马尔科夫过程，可用于描述随机过程统计特性的概率模型。这里所说的随机过程，一般都是很长的随机序列，它是一个双重随机过程。其中之一就是马尔科夫链，描述状态之间的转移；另一个随机过程是描述状态和观测变量之间的统计对应关系。如图6-19所示为HSMM组成框架。

图6-19 HSMM组成框架

一个有N个状态的HSMM由以下六元组参数组成：

λ={N，M，π，A，B，P_i（d）} （6-37）

式中，N为模型的状态数；M为观测符号数，即每个状态可能包含的观测值数目；π为初始状态概率分布；A为状态转移概率矩阵；B为由输出概率密度函数组成的向量；P_j（d）为状态持续时间密度函数。具体地：

A={a_ij} （6-38）

式中，a_ij表示从宏观状态h_i转移到宏观状态h_j的概率，且有0≤a_ij≤1，1≤i，j≤N。

B={b_j（k）} （6-39）

式中，b_j（k）为输出概率密度函数，是状态S_j的观测符号概率，它是随机变量或随机矢量在各状态的观测概率空间中的分布，通常以混合高斯分布表示：

式中，O为观测向量序列O=[O（1），O（2），…，O（k）]；G为每个状态可能包含的高斯元个数；N代表多维高斯概率密度函数（正态分布）；ω_jg代表第j个状态第g个混合高斯权重；μ_jg代表第j个状态第g个混合高斯值的均值；U_jg代表第j个状态第g个混合高斯方差矩阵。

初始状态概率分布π={π₁，π₂，…，π_N}用于描述观测序列O在t=1时刻所处状态q属于模型中各状态的概率分布，即

π_i=p（q_i=S_i），i=1，2，…，N （6-41）

状态持续时间密度函数P_j（d）表示状态S_j持续d个时间单元的概率，用单高斯分布函数来表示：

P_j（d）=N（d|μ_j，σ_j），0≤d≤D （6-42）

式中，d为状态持续时间单元；D为最大状态持续时间单元；N代表多维高斯概率密度函数（正态分布）；μ_j代表第j个状态持续时间的均值；σ_j代表第j个状态持续时间的方差。

2.基于HSMM的滚动轴承的故障预测方法

如图6-20所示为基于HSMM的滚动轴承寿命预测流程。

首先，采用前文介绍的小波相关特征尺度熵的特征提取方法建立滚动轴承全寿命数据的小波相关特征尺度熵向量序列；然后，利用得到的全寿命预测特征向量训练一个具有所有退化状态特征的HSMM模型，得到各个退化状态对应的状态转移矩阵和每个退化状态的状态持续时间密度函数的均值和方差，进而求取各状态的驻留时间单元，并结合退化状态识别方法确定滚动轴承当前所处的退化状态；最后，通过后向递归算法计算滚动轴承处于当前状态的剩余寿命。具体步骤如下：

1）全寿命HSMM模型训练样本的建立。首先，按照前文所介绍的小波相关特征尺度熵的预测特征提取方法提取全寿命过程中每个观测数据的小波相关特征尺度熵值；进而对滚动轴承的全寿命过程进行描述，并通过小波相关特征尺度熵值的大小趋势，确定全寿命HSMM建模的状态数N。然后，按照设定的时间单元选择多组最终状态模式相同并且等时间间隔的全寿命数据进行预处理；进而提取小波相关特征尺度熵向量序列，并建立多组表征全寿命的观测向量序列，把这些序列作为训练全寿命HSMM模型的训练样本。

2）全寿命HSMM模型的训练。首先利用上步建立的观测样本训练一个全寿命的HSMM模型HSMM_d（λ_d），获得模型中的状态转移矩阵{a_ll}和全寿命过程中各个状态持续时间概率函数的均值μ（h_l）和方差σ²（h_l）。然后由式（6-43）求得各个状态的持续时间单元D（h_l）。

D（h_l）=μ（h_l）+ρσ²（h_l） （6-43）

式中，。

3）剩余使用寿命的预测。由以下后向递归方程得到滚动轴承剩余使用寿命（Remaining Useful Life，RUL）（假设部件当前健康状态为退化状态l，RUL_l表示部件处于健康状态l时的剩余使用寿命）：

滚动轴承当前健康状态为l-1时：

RUL_l-1=a_l-1，l-1[D（h_l-1）+RUL_l]+a_l-1，lRUL_l(www.xing528.com)

滚动轴承当前健康状态为l-2时：

RUL_l-2=a_l-2，l-2[D（h_l-2）+RUL_l-1]+a_l-2，l-1[RUL_l-1]

滚动轴承当前健康状态为l-3时：

RUL_l-3=a_l-3，l-3[D（h_l-3）+RUL_l-2]+a_l-3，l-2[RUL_l-2]

……

图6-20 基于HSMM的滚动轴承寿命预测流程

滚动轴承当前健康状态为l时：

RUL_l=a_l，l[D（h_l）+RUL_l+1]+a_l，l+1[RUL_l+1]

3.滚动轴承试验验证

轴承实验台整体结构如图6-21所示。其基本组成包括实验台基座、轴承座（轴承）、轴、飞轮、电动机、联轴器，此外还包括一些辅助器件，如轴端挡圈、键等。本实验选用深沟球轴承622进行轴承的实验研究，其基本参数为：内圈直径100mm，外圈直径180mm，宽度34mm，滚动体数（球数）10个，滚动体直径为25.4mm。

实验选择电动机和变频器进行轴的转速控制。电动机为YVP系列变频调速异步电动机YVP90L-A，它与变频调速器配套使用，构成交流变频无级调速系统，在5～100Hz范围内恒转矩5～100Hz范围、变频调速运行，在恒功率变频调速运行，且运行平稳、无转矩脉动现象，起动转矩大、起动电流小。其主要性能指标有：额定功率4.5kW，额定转矩9.5N·m，额定电压380V。变频器具有多种控制方式，可适应不同场合控制需求，最大功率可达220kW，高输出转矩，低速时可达150N·m，频率解析度高达0.01Hz，且具有良好的通信界面；通过通信界面进行程序编辑，设置变频器的频率，以控制电动机的转速。另外，分别安装加速度传感器，采用多分析仪（图6-21b）和配套的分析软件测量振动信号。为了避免在轴承上设置测点影响轴承的真实寿命，并且能同时采集两个方向的数据。选择下面两个测点进行数据采集，测点位置如图6-22所示。测点一：安装在轴承座顶部径向位置；测点二：安装在轴承座顶部轴向位置。

图6-21 轴承实验台整体结构

图6-22 测点位置

图6-23 小波能谱熵描述全寿命周期和演化趋势

采用1#滚动轴承第1通道（径向）总计100组数据来描述轴承的全寿命过程。对数据进行规范化预处理后，提取该100组数据的小波能谱熵描述全寿命的演化趋势，如图6-23所示。从图6-23可以看出，在记录数为22处小波能谱熵值有明显的增大，经检测该处为状态转移点。根据全寿命趋势中的各个突变点，近似的把全寿命过程分为正常状态、退化状态1、退化状态2、退化状态3、退化状态4等五个阶段。

首先进行滚动轴承退化状态的识别，为了增加模型的适用性和鲁棒性，需要每种状态对应的多个观测样本进行模型训练。由于1#～4#滚动轴承最终的故障模式相同，按照图6-23近似的把全寿命过程分为五个阶段，任意选取1#～4#滚动轴承对应每个状态的20组数据作为各个状态的训练数据。五种状态的HSMM模型训练完成后，就可以构建一个状态分类器。于是，对于当前状态的特征向量序列，可以根据改进Viterbi算法计算出该特征向量序列在HSMM模型库中各个模型下的输出对数似然概率。然后根据Bayes分类准则，找出输出最大概率的模型作为该状态的退化类型。分别采用5#～8#滚动轴承五种已知状态对应的各10组数据形成观测序列样本对分类器进行检验，检验的结果见表6-2，结果较为理想。

表6-2 五种已知状态对应的各10组数据形成观测序列样本对分类器进行检验的结果

建立滚动轴承的全寿命观测样本，并进行训练可以得到五种状态的持续时间单元的均值和方差，其输出结果见表6-3。

表6-3 五种状态的持续时间单元的均值和方差输出结果

结合上述信息，由HSMM的后向递归算法可以计算出各个状态的剩余寿命RUL。假如当前状态为正常状态，则该状态剩余寿命的均值为Mean__RUL1=58.3378（时间单元）。同理，该状态剩余寿命的方差Var__RUL1=4.9641（时间单元）。这就表示，如果当前的状态是正常状态，此时滚动轴承的剩余寿命为58.3378个时间单元，置信区间为4.9641个时间单元。现利用5#～8#实验数据对上述剩余寿命预测方法进行验证。发现5#、7#、8#轴承实际测得正常状态时剩余寿命时间处于区间[Mean__RUL1-Var__RUL1，Mean__RUL1+Var__RUL1]内，6#号轴承正常状态的实际剩余寿命在上述区间外，预测结果相对较为理想。