CN104504296A - 高斯混合隐马尔可夫模型和回归分析的剩余寿命预测方法 - Google Patents

Abstract

本发明公开了一种高斯混合隐马尔可夫模型和回归分析的剩余寿命预测方法，具体步骤如下：(1)使用小波包算法分解进行特征矢量提取；(2)建立不同故障模式的高斯混合隐马尔科夫模型库；(3)建立不同故障模式的零部件失效时间模型；(4)零部件故障模式识别和失效状态评估；(5)利用回归分析预测剩余使用寿命。本发明融合数据驱动和概率统计的剩余寿命预测方法，充分利用各自的优势：隐马尔科夫模型对剩余使用寿命进行预测，具有随机性和有实时性；将零部件失效过程划分为多个阶段，采用概率统计和回归分析的手段对当前使用时间进行修正，提高了剩余寿命预测的精度。本发明具有预测精度高、运算速度快、实时性强、成本低等优点。

Description

高斯混合隐马尔可夫模型和回归分析的剩余寿命预测方法

技术领域

本发明属于故障诊断与寿命预测领域，具体涉及一种基于高斯混合隐马尔可夫模型和回归分析的剩余寿命预测方法。

背景技术

在机电产品竞争力的需求中，除了需要不断提升产品的实用性、可靠性、安全性，同时也要求不断降低产品的维护和保养成本。目前机电产品大多采用事后维修和定期维修，这两种维修方法虽具有可操作性，但由于没有考虑到设备实际的运行及老化状况，这些维护方式的效果其实并不理想，导致维护不足或维护过度的现象。

基于状态的维护(CBM)作为一种可行的解决方案，有助于降低了传统的定期维护、周期性保养和事后维护的成本。目前关于剩余寿命预测方法主要可分为以下三类：

(1)基于物理模型的寿命预测方法

该方法是基于分析模型来描述系统的行为，其优点在于提供了更高的精度。但是，现实中的系统往往是非线性，而且退化机制通常都是随机的，因而很难以分析模型的形式获得。

(2)基于可靠性的寿命预测方法

该方法使用的数据主要聚集在一个重要时期的经验反馈(维护和运行数据，故障时间等)，而且后者常用这些来估计剩余寿命时间。其优点在于当有相当数量的经验反馈数据是可用时，可靠性预测的方法是一种最简单实用的方法，然而其预测的精度是不如基于模型和数据驱动的方法的，尤其是当操作条件是可变的，或者是在一个缺乏经验数据的新系统的情况下时，这种问题尤为突出。

(3)基于数据驱动的寿命预测方法

该方法是将原始监测数据转变为相关的信息和系统的行为模型的，包括退化，而不是建立在全面的系统的物理和人类的知识基础之上的模型。该方法使用人工智能(AI)工具(神经元网络，贝叶斯网络，)或统计模型学习退化模型，来预测未来的健康状况和相关系统的剩余使用寿命。这种方法的原则主要是从监测的数据中抽象出行为模型，然后利用抽象出来的行为模型来评估当前系统的健康状态。

还有将上述方法的两者或三者融合在一起，产生的交叉方法，仍存在成本高、预测周期长、实时性不佳等问题。

发明内容

针对现有技术存在的上述技术问题，本发明的目的是提供一种预测精度高、运算速度快、实时性强、成本低的基于高斯混合隐马尔可夫模型和回归分析的剩余寿命预测方法。

为实现上述目的，本发明的一种高斯混合隐马尔可夫模型和回归分析的剩余寿命预测方法，具体步骤如下：

(1)使用小波包算法分解进行特征矢量提取；

小波包算法分解公式为：

$c_{j-1,k}=\underset{m=-\∞}{\overset{+\∞}{\mathrm{\Σ}}}{\hat{h}}_{m-2k}{c}_{j,m}$

$d_{j-1,k}=\underset{m=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{\hat{g}}_{m-2k}{c}_{j,m}$

其中，c_j,k为尺度系数，d_j,k为小波系数，分别是正交滤波器共轭系数；j,k是连续小波变换中的尺度参数和平移参数离散化取值，且j,k∈z；j表示小波包分解层数，即分辨率为2^j的小波包分解，j＝1,2,L,n；k,m是小波包分解在2^j子空间的平移量，表示分解到j层对应的节点数的平移量。

小波包算法重构公式为：

$c_{j,k}=\underset{m=-\∞}{\overset{+\∞}{\mathrm{\Σ}}}{h}_{k-2m}{c}_{j-1,m}+\underset{m=-\∞}{\overset{+\∞}{\mathrm{\Σ}}}{g}_{k-2m}{d}_{j-1,m}$

对于离散信号，利用小波包分解重构算法，得到某一频段或几个频段的信号成分，进而提取原始信号的特征矢量；

(2)建立不同故障模式的高斯混合隐马尔科夫模型库；

使用高斯混合模型来拟合各状态下的观测值概率密度函数，即

$b_j\left(o_t\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{w}_{\mathrm{jm}}{b}_{\mathrm{jm}}\left(o_t\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{w}_{\mathrm{jm}}N(o_t,{\mathrm{\μ}}_{\mathrm{jm}},C_{\mathrm{jm}}),1\≤j\≤N,1\≤m\≤M$

式中，M是状态s_j的高斯分量数目，相当于离散HMM中每个状态对应的观测值数目；w_jm是状态s_j的第m个高斯分布的权值，μ_jm和C_jm分别是状态s_j的第m个高斯分布的均值向量和协方差矩阵，b_jm(o_t)是状态s_j的第m个高斯分布，此分布是一个多维正态随机变量概率密度函数，其维数D就是特征向量的维数；

则多维正态随机变量概率密度函数为：

$\begin{array}{c}{b}_{\mathrm{jm}}\left(o_t\right)=N(o_t,{\mathrm{\μ}}_{\mathrm{jm}},C_{\mathrm{jm}})\\ =\frac{1}{\sqrt{{\left(2\mathrm{\π}\right)}^D\left|C_{\mathrm{jm}}\right|}}\exp (-\frac{1}{2}{(o_t-{\mathrm{\μ}}_{\mathrm{jm}})}^{\′}{C}_{\mathrm{jm}}^{-1}(o_t-{\mathrm{\μ}}_{\mathrm{jm}}))\end{array}$

γ_t(j,m)为在给定模型参数λ和观测值序列O的条件下，t时刻模型处于状态s_j并且对应该状态下的第m个高斯分布的联合概率，即

γ_t(j,m)＝P(q_t＝s_j,x_j,t＝X_j,m|O,λ)

式中，x_j,t表示t时刻状态s_j的高斯分布，X_j,m表示状态s_j的第m个高斯分布，γ_t(j,m)使用前向变量和后向变量求得，

${\mathrm{\γ}}_t(j,m)=\frac{{\mathrm{\α}}_t\left(j\right){\mathrm{\β}}_t\left(j\right)}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_t\left(i\right){\mathrm{\β}}_t\left(i\right)}\×\frac{w_{j,m}N(o_t,{\mathrm{\μ}}_{j,m},C_{j,m})}{\underset{n=1}{\overset{M_j}{\mathrm{\Σ}}}{w}_{j,n}N(o_t,{\mathrm{\μ}}_{j,n},C_{j,n})}$

根据变量γ_t(j,m)，重估权值w_j,m、均值向量μ_j,m和协方差矩阵C_j,m，公式如下：

${\hat{w}}_{j,m}=\frac{\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{\γ}}_t(j,m)}{\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\underset{n=1}{\overset{M_j}{\mathrm{\Σ}}}{\mathrm{\γ}}_t(j,n)},{\widehat{\mathrm{\μ}}}_{j,m}=\frac{\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{\γ}}_t(j,m)o_t}{\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{\γ}}_t(j,m)},{\mathrm{\Σ}}_{j,m}=\frac{\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{\γ}}_t(j,m)(o_t-{\mathrm{\μ}}_{j,m}){(o_t-{\mathrm{\μ}}_{j,m})}^T}{\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{\γ}}_t(j,m)};$

即可得出高斯混合隐马尔科夫模型参数重估均值向量μ_j,m、协方差矩阵C_j,m和权值w_j,m参数；

(3)建立不同故障模式的零部件失效时间模型；

齿轮磨损按四种损伤状态分为分为初期磨损、正常磨损、中度磨损、剧烈磨损4个时间段，每个阶段时间分布密度函数均满足正态分布，其概率密度函数为

$f({\mathrm{\&mu;}}_i,{\mathrm{\&sigma;}}_i,t)=\frac{1}{{\mathrm{\&sigma;}}_i\sqrt{2\mathrm{\&pi;}}}\exp [-\frac{{(t-{\mathrm{\&mu;}}_i)}^2}{2{\mathrm{\&sigma;}}_i^2}],0\&le;t<\&infin;$

式中μ_i是齿轮第i种磨损状态经历时间的期望值，σ_i为齿轮第i种磨损状态标准差，表示齿轮磨损状态的分散程度，μ_i,σ_i,t为时间单位；

(4)零部件故障模式识别和失效状态评估；

采用小波包分析提取当前的设备状态监测数据的特征向量，这些特征向量作为高斯混合隐马尔科夫模型的观测值向量，采用前向算法计算当前观测值向量在隐马尔科夫模型库的概率P(O|λ)，并根据Viterbi算法确定观测序列出现的最优路径，依次确定最大概率来判断零部件所处的故障模式，评估其失效状态和损伤程度；

(5)利用回归分析预测剩余使用寿命；

概率P(O|λ)记为p，经小波包和预测处理后可判断出齿轮磨损处于第i损伤状态，其概率密度函数参数可以确定为μ_i,σ_i；

记齿轮磨损最终损伤状态的概率密度函数参数为μ_z,σ_z，则存在下列关系

$\left\{\begin{array}{c}\mathrm{\&Delta;\&mu;}={\mathrm{\&mu;}}_z-{\mathrm{\&mu;}}_i\\ \mathrm{\&Delta;\&sigma;}={\mathrm{\&sigma;}}_z-{\mathrm{\&sigma;}}_i\end{array}\right.$

记齿轮磨损失效寿命为T，则剩余使用寿命T_RUL可表示为

T_RUL＝T-t_n+c×Δμ+d×Δσ+e×lnp

式中，c,d,e分别代表回归模型的回归系数。

本发明融合数据驱动和概率统计的剩余寿命预测方法，可以充分利用各自的优势，一方面，隐马尔科夫模型具有严谨的数据结构和可靠计算性能，它是在实时状态监测的基础上，建立描述故障状态转移的马氏链和反映故障状态与状态信息关系的随机过程，对零部件故障伤状态进行评估，对剩余使用寿命进行预测，具有随机性和有实时性；另一方面，考虑机电设备的载荷不确定对剩余寿命预测的影响，将零部件失效过程划分为多个阶段，采用概率统计和回归分析的手段对当前使用时间进行修正，提高了剩余寿命预测的精度。因此本发明具有预测精度高、运算速度快、实时性强、成本低等优点。

附图说明

图1为本发明预测方法的流程图；

图2为齿轮磨损剩余使用寿命测试平台示意图；

图3为本发明齿轮磨损故障原始信号的波形图；

图4为本发明齿轮磨损故障第三层小波包分解节点0-3的波形图；

图5为本发明齿轮磨损故障第三层小波包分解节点4-7的波形图；

图6为本发明齿轮磨损故障第三层小波分解与重构各频段节点能量图；

图7为本发明齿轮初期磨损HMM训练迭代过程图；

图8为本发明齿轮正常磨损HMM训练迭代过程图；

图9为本发明齿轮中度磨损HMM训练迭代过程图；

图10为本发明齿轮剧烈磨损HMM训练迭代过程图；

图11为本发明齿轮不同磨损损伤状态的概率密度函数图。

具体实施方式

如图1所示，本发明的高斯混合隐马尔可夫模型和回归分析的剩余寿命预测方法。该评估方法包括：

(1)使用小波包算法分解进行特征矢量提取。使用小波包分解对设备状态监测信号进行处理，提取零件不同故障模式的特征信息，作为高斯混合隐马尔可夫模型(GMM-HMM)的观察值。

(2)建立不同故障模式的高斯混合隐马尔科夫模型库。利用不同故障模式的历史数据进行样本训练，用来生成关键零部件不同状态下的隐马尔科夫模型库，属于离线学习阶段。

(3)建立不同故障模式的零部件失效时间模型。对不同故障模式下的时间历程进行概率统计，建立不同故障模式的失效时间模型，得到不同故障模式下零部件失效时间分布函数。

(4)零部件故障模式识别和失效状态评估。计算当前的设备状态监测数据在零部件故障HMM模型库的概率P(O|λ)，并根据出现的最大概率来确定零部件所处的故障模式及损伤程度。

(5)利用回归分析预测剩余使用寿命。其特征在于：根据步骤(3)不同故障模式下零部件失效时间正态分布函数参数，根据步骤(4)确定零部件的目前损伤状况及其相关的概率值，采用回归分析预测零部件剩余使用寿命。

接下来对各步骤的具体操作过程进行详细介绍。

步骤(1)的具体过程为：

关键零部件性能劣化状态和依据对象的工作性质来划分，如齿轮磨损状态(劣化状态)可分为：初期磨损、正常磨损、中度磨损、剧烈磨损。它们对应于隐马尔科夫模型(HMM)的状态空间，是一个不可逆过程。

在图2齿轮磨损剩余使用寿命测试平台采集的试验数据，如图3所示。利用小波分解和小波重构手段对实验数据进行处理，可以较好地用于非高斯噪声的消除。在小波包分解过程中，使用“Daubechies”小波族将从原始信号(振动信号)的第三个层次中的节点能量提取出来，对应齿轮不同磨损状态的特征向量，作为HMM的参数估计和状态预测的输入向量，如图4所示。

小波包分析Mallet算法的基本思想：假定已经计算出某函数或信号f(t)∈L²(R)在分辨率2^j下的离散逼近A_jf(t)，则f(t)在分辨率2^j-1的离散逼近A_j-1f(t)可通过用离散低通滤波器H对A_jf(t)滤波获得。根据Mallat算法的分解思想，A_jf(t)分解为粗糙像(低频)A_j-1f(t)与细节D_j-1f(t)之和，即有

A_jf(t)＝A_j-1f(t)+D_j-1f(t)    (1)

一旦多分辨分析的尺度函数φ(t)被确定，小波函数ψ(t)也可以构造出来。当构造的尺度函数是正交的，则此多分辨分析为正交多分辨分析。Mallat算法正是基于正交多分辨分析进行的。函数f(t)的离散逼近A_jf(t)和细节部分D_jf(t)可分别表示为：

$A_{j}f\left(t\right)=\underset{-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{c}_{j,k}{\mathrm{\&phi;}}_{j,k}\left(t\right)---\left(2\right)$

$D_{j}f\left(t\right)=\underset{-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{d}_{j,k}{\mathrm{\&psi;}}_{j,k}\left(t\right)---\left(3\right)$

式中c_j,k和d_j,k分别为2^j分辨率下的尺度(或粗糙像)系数和小波(或细节)系数。j,k是连续小波变换中的尺度参数和平移参数离散化取值，且j,k∈z；j表示小波包分解层数，即分辨率为2^j的小波包分解，j＝1,2,L,n；k是小波包分解在2^j子空间的平移量，表示分解到j层对应的节点数的平移量。

由多分辨分析的伸缩性及包容性，可得到尺度函数和小波函数的双尺度方程。

$\mathrm{\&phi;}\left(t\right)=\sqrt{2}\underset{k=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{h}_k\mathrm{\&phi;}(2t-k)---\left(4\right)$

$\mathrm{\&psi;}\left(t\right)=\sqrt{2}\underset{k=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{g}_k\mathrm{\&phi;}(2t-k)---\left(5\right)$

式(4)中，h_k为尺度函数的双尺度方程展开系数，也称低通滤波器系数。式(5)中，g_k为小波双尺度方程展开系数，也称高通滤波器系数。g_k与h_k二者关系可表示为

g_k＝(-1)^k-1h_1-k,k∈Z    (6)

根据多分辨分析通过一系列计算，小波包Mallat算法分解公式可表示为：

$c_{j-1,k}=\underset{m=-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{\hat{h}}_{m-2k}{c}_{j,m}---\left(7\right)$

$d_{j-1,k}=\underset{m=-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{\hat{g}}_{m-2k}{c}_{j,m}---\left(8\right)$

式(7)、式(8)中分别是正交滤波器共轭系数，m与k表示意义相同。

小波包Mallat算法重构公式可表示为：

$c_{j,k}=\underset{m=-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{h_{k-2m}c}_{j-1,m}+\underset{m=-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{g}_{k-2m}{d}_{j-1,m}---\left(9\right)$

对于离散信号，利用小波包分解重构算法，得到某一频段或几个频段的信号成分，进而提取原始信号的特征矢量。事实上，对应于每一个给定零部件的特定状态的原始数据，通过采用小波包分解转化为一个能量矩阵，矩阵每个列向量对应的原始信号的映像。由图3、图4、图5、图6可以看出，通过小波包分解重构得到的各层高频系数的信号相比原始更为缓和，达到消噪的目的，这也更直观的展现信号小波分解与重构在信号消噪中的重要作用。

步骤(2)的具体过程为：

HMM是双重随机过程，首先描述状态之间转移的Markov链是随机过程，另外描述状态和观测变量之间的统计对应关系也是随机过程。站在观察者的角度，只能看到观察值，不能直接看到状态，而是通过一个随机过程去感知状态的存在及其特性。隐马尔科夫模型(HMM)的预测方法，与人工神经网络、模糊逻辑、专家系统等黑箱模型相比，HMM有一个比较明确清晰的预测过程，并且运算速度快，预测精度高。

通常情况下，HMM模型认为某种状态下的观测值都是离散固定值信号，并且使用离散概率密度的方法对转移概率和观察概率进行建模。然而在实际的状态监测和故障诊断过程中，作为观测值的特征信号通常是连续变化的，并非几个确定的值。虽然连续信号可以通过各种离散方法(如矢量量化)进行离散化编码，但是该过程可能会造成有效信息的大量损失。因此，在状态监测和故障诊断过程中，使用带连续观测值概率密度函数的HMM(简称CHMM)比观测值为离散符号的HMM(简称HMM)将更有优势。在实际的故障诊断信号处理和特征提取过程中，故障特征信号往往是一个多维特征向量。这些连续变化的多维特征信号可以用几个概率密度函数混合来表征，由于高斯混合模型(Gaussian Mixture Model,GMM)可以无限逼近任意分布，因此，通常使用GMM来拟合各状态下的观测值概率密度函数，即

$b_j\left(o_t\right)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{w}_{\mathrm{jm}}{b}_{\mathrm{jm}}\left(o_t\right)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{w}_{\mathrm{jm}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jm}},C_{\mathrm{jm}}),1\&le;j\&le;N,1\&le;m\&le;M---\left(10\right)$

式中，M是状态s_j的高斯分量数目，相当于离散HMM中每个状态对应的观测值数目；w_jm是状态s_j的第m个高斯分布的权值，μ_jm和C_jm分别是状态s_j的第m个高斯分布的均值向量和协方差矩阵，b_jm(o_t)是状态s_j的第m个高斯分布，此分布是一个多维正态随机变量概率密度函数，其维数D就是特征向量的维数。根据前面对多维特征向量的分析处理，多维正态随机变量概率密度函数为

$\begin{array}{c}{b}_{\mathrm{jm}}\left(o_t\right)=N(o_t,{\mathrm{\&mu;}}_{\mathrm{jm}},C_{\mathrm{jm}})\\ =\frac{1}{\sqrt{{\left(2\mathrm{\&pi;}\right)}^D\left|C_{\mathrm{jm}}\right|}}\exp (-\frac{1}{2}{(o_t-{\mathrm{\&mu;}}_{\mathrm{jm}})}^{\&prime;}{C}_{\mathrm{jm}}^{-1}(o_t-{\mathrm{\&mu;}}_{\mathrm{jm}}))\end{array}---\left(11\right)$

从上面分析可知，对于GMM-HMM的观察值矩阵B是利用多维高斯密度函数进行建模的，即用均值向量μ_jm和协方差矩阵C_jm和混合系数(权值)w_jm来表征b_j(o_t)。一个HMM模型通常采用一个三元组λ＝(π,A,B)来表示，那么一个GMM-HMM也可以用类似的参数来表示。GMM-HMM的前向-后向算法、Viterbi算法和Baum-Welch算法与离散HMM基本一致，只是后者在计算和重估观察值概率分布矩阵B时更为复杂一些。

定义变量γ_t(j,m)为在给定模型参数λ和观测值序列O的条件下，t时刻模型处于状态s_j并且对应该状态下的第m个高斯分布的联合概率。

γ_t(j,m)＝P(q_t＝s_j,x_j,t＝X_j,m|O,λ)    (12)

式中，x_j,t表示t时刻状态s_j的高斯分布，X_j,m表示状态s_j的第m个高斯分布。γ_t(j,m)可以使用前向变量和后向变量求得。

${\mathrm{\&gamma;}}_t(j,m)=\frac{{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}\&times;\frac{w_{j,m}N(o_t,{\mathrm{\&mu;}}_{j,m},C_{j,m})}{\underset{n=1}{\overset{M_j}{\mathrm{\&Sigma;}}}{w}_{j,n}N(o_t,{\mathrm{\&mu;}}_{j,n},C_{j,n})}---\left(13\right)$

然后，根据变量γ_t(j,m)，可以重估权值w_j,m、均值向量μ_j,m和协方差矩阵C_j,m，公式如下：

${\hat{w}}_{j,m}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,m)}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{M_j}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,n)}---\left(14\right)$

${\widehat{\mathrm{\&mu;}}}_{j,m}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,m)o_t}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,m)}---\left(15\right)$

${\mathrm{\&Sigma;}}_{j,m}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,m)(o_t-{\mathrm{\&mu;}}_{j,m}){(o_t-{\mathrm{\&mu;}}_{j,m})}^T}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,m)}---\left(16\right)$

在故障诊断中，通常采集到的是连续信号的离散值，使用连续HMM将更有优势。由传感器记录的原始数据信号经过预处理，使用WPD提取出每个节点的能量；然后节点的能量被用来估计参数(π，A，B)和GMM-HMM的时间参数(留在每个状态的持续时间)。采用高斯混合模型(Gaussian Mixture Model,GMM)拟合各种典型故障状态下特征向量的观测值概率密度函数，研究基于GMM-HMM关键零部件故障诊断与状态评估算法，利用Baum-Welch重估算法，观察值序列在最大似然估计(EM算法)下学习训练，得到不同故障模式下HMM模型参数，建立各种状态下HMM模型库。此阶段属于离线学习阶段。齿轮磨损四种损伤状态(初期磨损、正常磨损、中度磨损、剧烈磨损)HMM训练过程如图7、图8、图9、图10所示，其中齿轮中度磨损得到的HMM模型参数如下：

$\mathrm{\&pi;}=\left(\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right),A=\left(\begin{array}{cccc}0.9687& 0.0313& 0& 0\\ 0& 0.7627& 0.2373& 0\\ 0& 0& 0.9959& 0.0041\\ 0& 0& 0& 1.0000\end{array}\right)$

高斯混合隐马尔科夫模型参数重估均值向量μ_j,m和协方差矩阵C_j,m、权值w_j,m参数为下表1所示。

表1高斯混合隐马尔科夫模型参数重估表

步骤(3)的具体过程为：

现场收集设备运行故障数据记录、试验数据或利用故障仿真数据，根据现场、试验故障数据分析产品的故障分布规律，对不同故障模式下的时间历程进行概率统计，采用正态分布、威布尔分布等建立不同故障模式的失效时间模型，得到不同故障模式下零部件失效时间分布函数。假设齿轮磨损按四种损伤状态分为4个时间段，每个阶段时间分布密度函数均满足正态分布，其概率密度函数为

$f({\mathrm{\&mu;}}_i,{\mathrm{\&sigma;}}_i,t)=\frac{1}{{\mathrm{\&sigma;}}_i\sqrt{2\mathrm{\&pi;}}}\exp [-\frac{{(t-{\mathrm{\&mu;}}_i)}^2}{2{\mathrm{\&sigma;}}_i^2}],0\&le;t<\&infin;---\left(17\right)$

式中μ_i是齿轮第i种磨损状态经历时间的期望值，σ_i为齿轮第i种磨损状态标准差，表示齿轮磨损状态的分散程度，μ_i,σ_i,t可以是分钟、小时、天等时间单位。齿轮磨损四种损伤状态的正态分布函数f(μ_i,σ_i,t)，如图7、图8、图9、图10所示。从图7、图8、图9、图10可看出齿轮磨损四种损伤的均值基本是先后时间关系，也说明齿轮磨损是不可逆转的时间历程。

步骤(4)的具体过程为：

采用小波包分析提取当前的设备状态监测数据的特征向量，这些特征向量作为GMM-HMM的观测值向量，采用前向算法计算当前观测值向量在HMM模型库的概率P(O|λ)，并根据Viterbi算法确定观测序列出现的最优路径，依次确定最大概率来判断零部件所处的故障模式，评估其失效状态和损伤程度。Viterbi算法解码寻找状态最优路径如图7所示。利用Viterbi算法，计算在HMM模型中某个退化状态产生某些观察序列的概率，确定观测序列最可能的隐藏状态。由传感器检测当前的齿轮磨损状态监测数据，经小波包分解重构处理，提取三层小波包节点能量信号，作为观测矢量输入到HMM模型库，计算观测矢量在HMM库出现的概率，并根据出现的最大概率来评估齿轮磨损的损伤状态。此阶段属于在线预测阶段。

步骤(5)的具体过程为：

建立基于回归分析的剩余寿命预测模型，利用历史数据对回归参数进行辨识，回归参数包括与零部件使用时间、设备状态监测数据在HMM模型库的概率P(O|λ)、不同故障模式下零部件失效时间分布函数(μ_i,σ_i)等有关的系数。根据第3、4步结果，利用不同故障模式下零部件失效时间分布函数，当前时刻的设备状态监测数据在HMM模型库的概率P(O|λ)，利用回归分析模型预测零部件剩余使用寿命。

齿轮磨损损伤状态可分为初期磨损、正常磨损、中度磨损、剧烈磨损，每个损伤状态时间概率密度函数如式(17)。若当前时刻t_n的状态监测数据在HMM模型库的概率P(O|λ)设为p，经小波包和预测处理后可判断出齿轮磨损处于第i损伤状态，其概率密度函数参数可以确定为μ_i,σ_i。

假设齿轮磨损最终损伤状态的概率密度函数参数为μ_z,σ_z，则存在下列关系

$\left\{\begin{array}{c}\mathrm{\&Delta;\&mu;}={\mathrm{\&mu;}}_z-{\mathrm{\&mu;}}_i\\ \mathrm{\&Delta;\&sigma;}={\mathrm{\&sigma;}}_z-{\mathrm{\&sigma;}}_i\end{array}\right.---\left(18\right)$

又假设齿轮磨损失效寿命为T，则剩余使用寿命T_RUL可表示为

T_RUL＝T-t_n+c×Δμ+d×Δσ+e×ln p    (19)

式(19)中，c,d,e分别代表回归模型的回归系数。

采用基于高斯混合隐马尔可夫模型和回归分析的方法对齿轮磨损失效过程中剩余寿命进行预测，误差精度可达到9-10％。

本发明提出的基于高斯混合隐马尔可夫模型和回归分析的剩余寿命预测方法，融合数据驱动和概率统计的剩余寿命预测方法，可以充分利用各自的优势。一方面，隐马尔科夫模型具有严谨的数据结构和可靠计算性能，它是在实时状态监测的基础上，建立描述故障状态转移的马氏链和反映故障状态与状态信息关系的随机过程，对零部件故障伤状态进行评估，对剩余使用寿命进行预测，具有随机性和有实时性。另一方面，考虑机电设备的载荷不确定对剩余寿命预测的影响，将零部件失效过程划分为多个阶段，采用概率统计和回归分析的手段对当前使用时间进行修正，提高了剩余寿命预测的精度。

Claims (1)

1.一种高斯混合隐马尔可夫模型和回归分析的剩余寿命预测方法，其特征在于，具体步骤如下：

(1)使用小波包算法分解进行特征矢量提取；

小波包算法分解公式为：

$c_{j-1,k}=\underset{m=-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{\hat{h}}_{m-2k}{c}_{j,m}$

$d_{j-1,k}=\underset{m=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{\hat{g}}_{m-2k}{c}_{j,m}$

其中，c_j,k为尺度系数，d_j,k为小波系数，分别是正交滤波器共轭系数；j,k是连续小波变换中的尺度参数和平移参数离散化取值，且j,k∈z；j表示小波包分解层数，即分辨率为2^j的小波包分解，j＝1,2,L,n；k,m是小波包分解在2^j子空间的平移量，表示分解到j层对应的节点数的平移量。

小波包算法重构公式为：

$c_{j,k}=\underset{m=-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{h}_{k-2m}{c}_{j-1,m}+\underset{m=-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{g}_{k-2m}{d}_{j-1,m}$

对于离散信号，利用小波包分解重构算法，得到某一频段或几个频段的信号成分，进而提取原始信号的特征矢量；

(2)建立不同故障模式的高斯混合隐马尔科夫模型库；

使用高斯混合模型来拟合各状态下的观测值概率密度函数，即

$b_j\left(o_t\right)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{w}_{\mathrm{jm}}{b}_{\mathrm{jm}}\left(o_t\right)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{w}_{\mathrm{jm}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jm}},C_{\mathrm{jm}}),1\&le;j\&le;N,1\&le;m\&le;M$

式中，M是状态s_j的高斯分量数目，相当于离散HMM中每个状态对应的观测值数目；w_jm是状态s_j的第m个高斯分布的权值，μ_jm和C_jm分别是状态s_j的第m个高斯分布的均值向量和协方差矩阵，b_jm(o_t)是状态s_j的第m个高斯分布，此分布是一个多维正态随机变量概率密度函数，其维数D就是特征向量的维数；

则多维正态随机变量概率密度函数为：

$\begin{array}{c}{b}_{\mathrm{jm}}\left(o_t\right)=N(o_t,{\mathrm{\&mu;}}_{\mathrm{jm}},C_{\mathrm{jm}})\\ =\frac{1}{\sqrt{{\left(2\mathrm{\&pi;}\right)}^D\left|C_{\mathrm{jm}}\right|}}\exp (-\frac{1}{2}{(o_t-{\mathrm{\&mu;}}_{\mathrm{jm}})}^{\&prime;}{C}_{\mathrm{jm}}^{-1}(o_t-{\mathrm{\&mu;}}_{\mathrm{jm}}))\end{array}$

γ_t(j,m)为在给定模型参数λ和观测值序列O的条件下，t时刻模型处于状态s_j并且对应该状态下的第m个高斯分布的联合概率，即

γ_t(j,m)＝P(q_t＝s_j,x_j,t＝X_j,m|O,λ)

式中，x_j,t表示t时刻状态s_j的高斯分布，X_j,m表示状态s_j的第m个高斯分布，γ_t(j,m)使用前向变量和后向变量求得，

${\mathrm{\&gamma;}}_t(j,m)=\frac{{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}\&times;\frac{w_{j,m}N(o_t,{\mathrm{\&mu;}}_{j,m},C_{j,m})}{\underset{n=1}{\overset{M_j}{\mathrm{\&Sigma;}}}{w}_{j,n}N(o_t,{\mathrm{\&mu;}}_{j,n},C_{j,n})}$

根据变量γ_t(j,m)，重估权值w_j,m、均值向量μ_j,m和协方差矩阵C_j,m，公式如下：

${\hat{w}}_{j,m}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,m)}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{M_j}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,n)},{\widehat{\mathrm{\&mu;}}}_{j,m}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,m)o_t}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,m)},{\mathrm{\&Sigma;}}_{j,m}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,m)(o_t-{\mathrm{\&mu;}}_{j,m}){(o_t-{\mathrm{\&mu;}}_{j,m})}^T}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,m)};$

即可得出高斯混合隐马尔科夫模型参数重估均值向量μ_j,m、协方差矩阵C_j,m和权值w_j,m参数；

(3)建立不同故障模式的零部件失效时间模型；

齿轮磨损按四种损伤状态分为分为初期磨损、正常磨损、中度磨损、剧烈磨损4个时间段，每个阶段时间分布密度函数均满足正态分布，其概率密度函数为

$f({\mathrm{\&mu;}}_i,{\mathrm{\&sigma;}}_i,t)=\frac{1}{{\mathrm{\&sigma;}}_i\sqrt{2\mathrm{\&pi;}}}\exp [-\frac{{(t-{\mathrm{\&mu;}}_i)}^2}{2{\mathrm{\&sigma;}}_i^2}],0\&le;t<\&infin;$

式中μ_i是齿轮第i种磨损状态经历时间的期望值，σ_i为齿轮第i种磨损状态标准差，表示齿轮磨损状态的分散程度，μ_i,σ_i,t为时间单位；

(4)零部件故障模式识别和失效状态评估；

采用小波包分析提取当前的设备状态监测数据的特征向量，这些特征向量作为高斯混合隐马尔科夫模型的观测值向量，采用前向算法计算当前观测值向量在隐马尔科夫模型库的概率P(O|λ)，并根据Viterbi算法确定观测序列出现的最优路径，依次确定最大概率来判断零部件所处的故障模式，评估其失效状态和损伤程度；

(5)利用回归分析预测剩余使用寿命；

概率P(O|λ)记为p，经小波包和预测处理后可判断出齿轮磨损处于第i损伤状态，其概率密度函数参数可以确定为μ_i,σ_i；

记齿轮磨损最终损伤状态的概率密度函数参数为μ_z,σ_z，则存在下列关系

$\left\{\begin{array}{c}\mathrm{\&Delta;\&mu;}={\mathrm{\&mu;}}_z-{\mathrm{\&mu;}}_i\\ \mathrm{\&Delta;\&sigma;}={\mathrm{\&sigma;}}_z-{\mathrm{\&sigma;}}_i\end{array}\right.$

记齿轮磨损失效寿命为T，则剩余使用寿命T_RUL可表示为

T_RUL＝T-t_n+c×Δμ+d×Δσ+e×ln p

式中，c,d,e分别代表回归模型的回归系数。