CN104793606A - 基于改进的kpca和隐马尔科夫模型的工业故障诊断方法 - Google Patents

Abstract

本发明公开了一种基于改进的KPCA和隐马尔科夫模型的工业故障诊断方法，属于工业过程监控与诊断技术领域。本发明通过相似性分析方法的引入大大提高了KPCA在大样本情况下的计算效率，并且利用隐马尔科夫模型极强的动态过程时间序列建模能力和时序模式分类能力，对工业过程的故障进行分类。因此与其它现有的方法相比，由于充分考虑了工业数据的非线性特性以及海量数据的特点，本发明方法不仅可以降低计算的复杂度，而且可以更有效处理过程的非线性特征，从而对于非线性工业故障诊断具有更高的准确率。

Description

基于改进的KPCA和隐马尔科夫模型的工业故障诊断方法

技术领域

本发明属于工业过程监控与故障诊断领域，特别涉及一种基于改进的KPCA和隐马尔科夫模型的工业故障诊断方法。

背景技术

随着工业过程复杂性的增长，工业过程监测和诊断的有效性对于保障生产过程安全、维持产品质量和优化产品利益变得日益重要。

对于过程监控和故障诊断问题，传统的方法大多采用多元统计过程监控技术(Multivariable Statistical Process Monitoring,MSPM)，其中以主元分析(Principal Component Analysis,PCA)、偏最小二乘(Partial Least Squares,PLS)和独立成分分析(Independent Component Analysis,ICA)为代表等方法已在工业过程监控中得到了成功的应用。传统的PCA、ICA等方法均假设过程变量间的关系是线性的，但是实际中测量变量难以满足这个假设条件，常呈现强烈的非线性特性。尽管核方法的引入，如核ICA，核PCA(Kernel PCA)等方法被提出用于解决变量间的非线性，但是，上述方法具有以下缺点，原始数据从输入空间映射到高维特征空间后变得冗余，而且核矩阵是一个以样本数目为大小的方阵。随着样本数目的增多，计算量不断增大，而在工业过程，样本数目往往是巨大的，因此用原始的KPCA来提取变量间的非线性关系是困难的，由此，可能会引起故障诊断发生错误。

发明内容

本发明的目的在针对现有技术的不足，提供一种基于改进的KPCA和隐马尔科夫模型的工业故障诊断方法，通过相似性分析方法的引入大大提高了KPCA在大样本情况下的计算效率，并且利用隐马尔科夫模型极强的动态过程时间序列建模能力和时序模式分类能力。因此，本方法在工业过程具有强非线性以及存在海量数据的情况下，可以实现有效地故障诊断。

一种基于改进的KPCA和隐马尔科夫模型的工业故障诊断方法，该方法的步骤如下：

步骤一：离线建模，对工业过程采集的离线数据，在原空间进行相似性分析，去除相似样本，得到数据集Ω₁，引入核函数φ将数据集Ω₁投影到高维特征空间，同样在特征空间进行数据的相似性分析，去除相似样本，得到数据集Ω₂，根据数据集Ω₂中数据计算相应的核矩阵以及核主元，将得到的核主元作为隐马尔科夫模型(HMM)的观测序列，训练隐马尔科夫模型；

步骤二：在线诊断，对在线采集的数据在原空间进行相似性分析，去除相似样本，同样在特征空间进行数据的相似性分析，去除相似样本，计算相应的核矩阵以及核主元，即得到相应的观测序列后，选择最匹配的HMM模型，从而判断故障类型。

步骤一所述的离线建模过程如下：

1)工业过程采集的第i个时刻监测数据 $N_i={[{\tilde{x}}_1,{\tilde{x}}_2,...,{\tilde{x}}_{n_i}]}^T,i=1,...,t_1(n_i<t_1),$ 其中n_i表示样本个数；

2)在原输入空间进行相似性分析，去除相似样本，得到数据集Ω₁，

2.1)初始时只有2个数据点，即新数据点根据以下提出的相似性指数公式进行相似性的判断，

$S_i=\frac{1}{2}{e}^{\left|\right|x_{\mathrm{new}}-{\tilde{x}}_i\left|\right|}+\frac{1}{2}{\left(\frac{{(x_{\mathrm{new}}-{\tilde{x}}_i)}^T(x_{\mathrm{new}}-{\tilde{x}}_{i-1})}{{\left|\right|x_{\mathrm{new}}-{\tilde{x}}_i\left|\right|}_2\left|\right|x_{\mathrm{new}}-{\tilde{x}}_{i-1}\left|\right|}\right)}^2---\left(1\right)$

其中x_new表示新数据点；

2.2)如果新数据点被引入，否则该样本点被舍弃，即N_i＝N_i-1，

其中γ₀是预先设定好的一个足够小的值，满足

2.3)对原始输入空间进行相似性分析之后，被保留的数据点个数为n₁，所以我们获得了数据集Ω₁；

3)引入核函数φ将数据集Ω₁投影到高维特征空间，在特征空间进行相似性分析，去除相似样本，得到数据集Ω₂，

3.1)特征空间中相似性指数公式定义如下：

$S_f=\frac{1}{2}{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right)\left|\right|}_2+\frac{1}{2}{\left(\frac{{(\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right))}^T(\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{i-1}\right))}{{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right)\left|\right|}_2{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{i-1}\right)\left|\right|}_2}\right)}^2---\left(2\right),$

3.2)如果新数据点被引入，否则该样本点被舍弃，

其中γ₁是预先设定好的一个足够小的值，满足

3.3)对特征空间进行相似性分析之后，被保留的数据点个数为n₂，所以我们获得了数据集Ω₂；

4)对数据集Ω₂，计算核矩阵K以及核主元，

4.1)计算核矩阵

K_ij＝<Φ(x_i),Φ(x_j)>＝K(x_i,x_j)   (3)

其中使用径向基核函数σ＝rm，r为常数；

4.2)对核矩阵K进行中心化处理，

$\tilde{K}=K-1_{n_i}K-K{1}_{n_i}+1_{n_i}K{1}_{n_i}---\left(4\right)$

其中

4.3)计算主成分tk

$t_k=<v_k,\mathrm{\&Phi;}\left(x\right)>=\underset{i=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^k<\mathrm{\&Phi;}\left(x_i\right),\mathrm{\&Phi;}\left(x\right)>=\underset{i=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^{k}K(x_i,x)---\left(5\right)$

5)将得到的核主元作为隐马尔科夫模型的观测序列O，训练隐马尔科夫模型，得到其参数λ＝(A,B,π,M,N)；

其中，A是隐含状态转移概率矩阵，描述了HMM模型中各个状态之间的转移概率，见下式(6)～(7)

$\begin{array}{c}{a}_{S_i,S_j}=a_{i,j}=P[q_{t+1}=S_j|q_t=S_i],1\&le;i,j\&le;N\\ \underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{i,j}=\mathrm{1,1}\&le;i\&le;N\end{array}---\left(7\right)$

式中，简记为a_i,j，表示在t时刻，状态为S_i的条件下，在t+1时刻状态是S_j的概率；

B是观测值概率密度矩阵，见下式

$\begin{array}{c}B=\{b_j\left(O\right),j=\mathrm{1,2},...N\}\\ b_j\left(O\right)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}H[O,{\mathrm{\&mu;}}_{\mathrm{jm}},U_{\mathrm{jm}}],1\&le;j\&le;N\end{array}---\left(9\right)$

其中O是观测向量，C_jm是隐含状态j的第m个混合成分的混合系数，H是对数凹或椭圆对称密度，μ_jm是隐含状态S_j的第m个混合成分的均值向量，U_jm是隐含状态S_j的第m个混合成分的协方差矩阵，C_jm满足随机约束：

$\begin{array}{c}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}=\mathrm{1,1}\&le;j\&le;N\\ C_{\mathrm{jm}}\&GreaterEqual;\mathrm{0,1}\&le;j\&le;N,1\&le;m\&le;M\end{array}---\left(10\right)$

C_jm、μ_jm、U_jm的重估是这样的：C_jm的重估值是HMM在隐含状态S_j带有第k个混合成分的次数的期望除以HMM处于隐含状态S_j的次数的期望，μ_jm的重估值是由第k个混合成分说明的观测矢量的部分期望值，U_jm是由第k个混合成分说明的观测矢量的部分协方差,见公式(11)～(13)，

$\begin{array}{c}{\stackrel{\&OverBar;}{C}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}\\ {\mathrm{\&gamma;}}_t(j,k)=\left[\frac{{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}\right]\left[\frac{C_{\mathrm{jk}}H[O_t,{\mathrm{\&mu;}}_{\mathrm{jk}},U_{\mathrm{jk}}]}{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}H[O_t,{\mathrm{\&mu;}}_{\mathrm{jm}},U_{\mathrm{jm}}]}\right]\end{array}---\left(11\right)$

${\stackrel{\&OverBar;}{\mathrm{\&mu;}}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)\&CenterDot;O_t}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(12\right)$

${\stackrel{\&OverBar;}{U}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)\&CenterDot;(O_t-{\mathrm{\&mu;}}_{\mathrm{jk}}){(O_t-{\mathrm{\&mu;}}_{\mathrm{jk}})}^{\&prime;}}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(13\right)$

π是初始状态概率矩阵，N是隐含状态的数目，M是每个隐含状态下的高斯混合成分的数目。

步骤二所述的在线诊断程如下：

a)在线采集工业过程的数据；

b)在原空间进行相似性分析，去除相似样本，得到数据集Ω₁，引入核函数φ将数据集Ω₁投影到高维特征空间，同样在高维特征空间进行数据的相似性分析，去除相似样本，得到数据集Ω₂，根据数据集Ω₂中数据计算相应的核矩阵以及核主元；

c)以核主元作为观测序列，选择最匹配的HMM模型，从而判断故障类型。

所述的工业过程具有非线性。

所述的工业故障为高炉冶炼过程故障。

一种所述的方法用于高炉冶炼过程故障诊断。

本发明有以下优势：

1.本发明首次提出一种基于改进的KPCA和隐马尔科夫模型的工业故障诊断方法，实现对工业过程的故障诊断；

2.本发明能够解决过程数据存在的非线性，从而能够更有效地诊断。

3.本发明能够提高核矩阵的计算速度，从而能够更有效地诊断。

附图说明

图1是本发明方法的流程框图。

具体实施方式

本发明提出的一种基于改进的KPCA和隐马尔科夫模型的工业故障诊断方法，如图1所示，包括以下各步骤：

步骤一所述的离线建模过程如下：

1)工业过程采集的第i个时刻监测数据 $N_i={[{\tilde{x}}_1,{\tilde{x}}_2,...,{\tilde{x}}_{n_i}]}^T,i=1,...,t_1(n_i<t_1),$ 其中n_i表示样本个数；

2)在原输入空间进行相似性分析，去除相似样本，得到数据集Ω₁

2.1)初始时只有2个数据点，即新数据点根据以下提出的相似性指数公式进行相似性的判断

$S_i=\frac{1}{2}{e}^{\left|\right|x_{\mathrm{new}}-{\tilde{x}}_i\left|\right|}+\frac{1}{2}{\left(\frac{{(x_{\mathrm{new}}-{\tilde{x}}_i)}^T(x_{\mathrm{new}}-{\tilde{x}}_{i-1})}{{\left|\right|x_{\mathrm{new}}-{\tilde{x}}_i\left|\right|}_2\left|\right|x_{\mathrm{new}}-{\tilde{x}}_{i-1}\left|\right|}\right)}^2---\left(1\right)$

其中x_new表示新数据点。

2.2)如果新数据点被引入，否则该样本点被舍弃，即N_i＝N_i-1。其中γ₀是预先设定好的一个足够小的值，满足

2.3)对原始输入空间进行相似性分析之后，被保留的数据点个数为n₁，所以我们获得了数据集Ω₁。

3)引入核函数φ将数据集Ω₁投影到高维特征空间，在特征空间进行相似性分析，去除相似样本，得到数据集Ω₂

3.1)特征空间中相似性指数公式定义如下：

$S_f=\frac{1}{2}{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right)\left|\right|}_2+\frac{1}{2}{\left(\frac{{(\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right))}^T(\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{i-1}\right))}{{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right)\left|\right|}_2{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{i-1}\right)\left|\right|}_2}\right)}^2---\left(2\right)$

3.2)如果新数据点被引入，否则该样本点被舍弃。

其中γ₁是预先设定好的一个足够小的值，满足

3.3)对特征空间进行相似性分析之后，被保留的数据点个数为n₂，所以我们获得了数据集Ω₂。

4)对数据集Ω₂，计算核矩阵K以及核主元

4.1)计算核矩阵

K_ij＝<Φ(x_i),Φ(x_j)>＝K(x_i,x_j)   (3)

其中使用径向基核函数σ＝rm，r为常数；

4.2)对核矩阵K进行中心化处理

$\tilde{K}=K-1_{n_i}K-K{1}_{n_i}+1_{n_i}K{1}_{n_i}---\left(4\right)$

其中

4.3)计算主成分_tk

$t_k=<v_k,\mathrm{\&Phi;}\left(x\right)>=\underset{i=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^k<\mathrm{\&Phi;}\left(x_i\right),\mathrm{\&Phi;}\left(x\right)>=\underset{i=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^{k}K(x_i,x)---\left(5\right)$

5)将得到的核主元作为隐马尔科夫模型的观测序列O，训练隐马尔科夫模型，得到其参数λ＝(A,B,π,M,N)；

其中，A是隐含状态转移概率矩阵，描述了HMM模型中各个状态之间的转移概率，见下式(6)～(7)

$\begin{array}{c}{a}_{S_i,S_j}=a_{i,j}=P[q_{t+1}=S_j|q_t=S_i],1\&le;i,j\&le;N\\ \underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{i,j}=\mathrm{1,1}\&le;i\&le;N\end{array}---\left(7\right)$

式中，简记为a_i,j，表示在t时刻，状态为S_i的条件下，在t+1时刻状态是S_j的概率；

B是观测值概率密度矩阵，见下式

$\begin{array}{c}B=\{b_j\left(O\right),j=\mathrm{1,2},...N\}\\ b_j\left(O\right)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}H[O,{\mathrm{\&mu;}}_{\mathrm{jm}},U_{\mathrm{jm}}],1\&le;j\&le;N\end{array}---\left(9\right)$

其中O是观测向量，C_jm是隐含状态j的第m个混合成分的混合系数，H是对数凹或椭圆对称密度，本发明采用混合高斯概率密度。μ_jm是隐含状态S_j的第m个混合成分的均值向量，U_jm是隐含状态S_j的第m个混合成分的协方差矩阵。C_jm满足随机约束：

$\begin{array}{c}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}=\mathrm{1,1}\&le;j\&le;N\\ C_{\mathrm{jm}}\&GreaterEqual;\mathrm{0,1}\&le;j\&le;N,1\&le;m\&le;M\end{array}---\left(10\right)$

C_jm、μ_jm、U_jm的重估是这样的：C_jm的重估值是HMM在隐含状态S_j带有第k个混合成分的次数的期望除以HMM处于隐含状态Sj的次数的期望。μ_jm的重估值是由第k个混合成分说明的观测矢量的部分期望值，U_jm是由第k个混合成分说明的观测矢量的部分协方差,见公式(11)～(13)。

$\begin{array}{c}{\stackrel{\&OverBar;}{C}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}\\ {\mathrm{\&gamma;}}_t(j,k)=\left[\frac{{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}\right]\left[\frac{C_{\mathrm{jk}}H[O_t,{\mathrm{\&mu;}}_{\mathrm{jk}},U_{\mathrm{jk}}]}{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}H[O_t,{\mathrm{\&mu;}}_{\mathrm{jm}},U_{\mathrm{jm}}]}\right]\end{array}---\left(11\right)$

${\stackrel{\&OverBar;}{\mathrm{\&mu;}}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)\&CenterDot;O_t}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(12\right)$

${\stackrel{\&OverBar;}{U}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)\&CenterDot;(O_t-{\mathrm{\&mu;}}_{\mathrm{jk}}){(O_t-{\mathrm{\&mu;}}_{\mathrm{jk}})}^{\&prime;}}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(13\right)$

π是初始状态概率矩阵，N是隐含状态的数目，M是每个隐含状态下的高斯混合成分的数目。

步骤二所述的在线诊断过程如下：

1)在线采集工业过程的数据；

2)在原空间进行相似性分析，去除相似样本，得到数据集Ω₁，引入核函数φ将数据集Ω₁投影到高维特征空间，同样在高维特征空间进行数据的相似性分析，去除相似样本，得到数据集Ω₂，根据数据集Ω₂中数据计算相应的核矩阵以及核主元；

3)以核主元作为观测序列，选择最匹配的HMM模型，从而判断故障类型；

实施例

钢铁冶炼作为国民经济中最重要的基础产业之一，是衡量一个国家的经济水平和综合国力的重要指标。而高炉炼铁是钢铁工业生产流程中最重要的环节，所以对大型高炉非正常工况诊断与安全运行方法进行研究具有重要意义。

高炉是一个巨大的密闭反应容器，其内部冶炼过程是在高温、高压条件下，经过一系列复杂的物理化学和传热反应，是一个典型的“黑箱”操作。正是由于高炉内部的复杂性，使得其监测过程具有非线性、非高斯性以及多模态等特性。因此，我们提出的方法对高炉故障监测具有适应性。下面结合柳钢2号高炉来说明本发明方法的有效性。

成立于1958年的柳钢炼铁厂，是一个有着56年辉煌历史的设备先进、装备水平较高的大型冶炼企业，主要产品为生铁，副产品有炉尘、炉渣、高炉煤气等。它拥有7座现代化高炉，高炉整体有效容积为11750立方米，其中2号高炉有效容积为2000立方米，是目前广西最大的高炉。新高炉投产后，炼铁厂将具备年产生铁1000万吨以上的综合能力。

接下来结合该具体过程对本发明的实施步骤进行详细地阐述：

步骤一：离线建模

1)假设传感器采集5种工况下的监测数据，包括1种正常工况，以及4种故障，悬料、塌料、管道行程以及炉凉。对于每种情况采集第i个时刻监测数据 $N_i={[{\tilde{x}}_1,{\tilde{x}}_2,...,{\tilde{x}}_{n_i}]}^T,i=1,...,t_1(n_i<t_1),$ 其中n_i表示样本个数；

2)对于5种情况都分别在原空间进行相似性分析，去除相似样本，得到数据集Ω₁

2.1)在初始时只有2个数据点，即新数据点根据以下提出的相似性指数公式进行相似性的判断

$S_i=\frac{1}{2}{e}^{\left|\right|x_{\mathrm{new}}-{\tilde{x}}_i\left|\right|}+\frac{1}{2}{\left(\frac{{(x_{\mathrm{new}}-{\tilde{x}}_i)}^T(x_{\mathrm{new}}-{\tilde{x}}_{i-1})}{{\left|\right|x_{\mathrm{new}}-{\tilde{x}}_i\left|\right|}_2\left|\right|x_{\mathrm{new}}-{\tilde{x}}_{i-1}\left|\right|}\right)}^2$

其中x_new表示新数据点。

2.2)如果新数据点被引入，否则该样本点被舍弃，即N_i＝N_i-1。其中γ₀是预先设定好的一个足够小的值，满足

2.3)对原始输入空间进行相似性分析之后，被保留的数据点个数为n₁，所以我们获得了数据集Ω₁。

3)引入核函数φ将数据集Ω₁投影到高维特征空间，在特征空间进行相似性分析，去除相似样本，得到数据集Ω₂

3.1)特征空间中相似性指数公式定义如下：

$S_f=\frac{1}{2}{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right)\left|\right|}_2+\frac{1}{2}{\left(\frac{{(\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right))}^T(\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{i-1}\right))}{{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right)\left|\right|}_2{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{i-1}\right)\left|\right|}_2}\right)}^2---\left(2\right)$

3.2)如果新数据点被引入，否则该样本点被舍弃。

其中γ₁是预先设定好的一个足够小的值，满足

3.3)对特征空间进行相似性分析之后，被保留的数据点个数为n₂，所以我们获得了数据集Ω₂。

4)对数据集Ω₂，计算核矩阵K以及核主元

4.1)计算核矩阵

K_ij＝<Φ(x_i),Φ(x_j)>＝K(x_i,x_j)   (3)

其中使用径向基核函数σ＝rm，r为常数；

4.2)对核矩阵K进行中心化处理

$\tilde{K}=K-1_{n_i}K-K{1}_{n_i}+1_{n_i}K{1}_{n_i}---\left(4\right)$

其中

4.3)计算主成分tk

$t_k=<v_k,\mathrm{\&Phi;}\left(x\right)>=\underset{i=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^k<\mathrm{\&Phi;}\left(x_i\right),\mathrm{\&Phi;}\left(x\right)>=\underset{i=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^{k}K(x_i,x)---\left(5\right)$

5)将得到的核主元作为隐马尔科夫模型的观测序列O，训练隐马尔科夫模型，得到其参数λ＝(A,B,π,M,N)；

其中，A是隐含状态转移概率矩阵，描述了HMM模型中各个状态之间的转移概率，见下式(6)～(7)

$\begin{array}{c}{a}_{S_i,S_j}=a_{i,j}=P[q_{t+1}=S_j|q_t=S_i],1\&le;i,j\&le;N\\ \underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{i,j}=\mathrm{1,1}\&le;i\&le;N\end{array}---\left(7\right)$

式中，简记为a_i,j，表示在t时刻，状态为S_i的条件下，在t+1时刻状态是S_j的概率；

B是观测值概率密度矩阵，见下式

$\begin{array}{c}B=\{b_j\left(O\right),j=\mathrm{1,2},...N\}\\ b_j\left(O\right)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}H[O,{\mathrm{\&mu;}}_{\mathrm{jm}},U_{\mathrm{jm}}],1\&le;j\&le;N\end{array}---\left(9\right)$

其中O是观测向量，C_jm是隐含状态j的第m个混合成分的混合系数，H是对数凹或椭圆对称密度，本发明采用混合高斯概率密度。μ_jm是隐含状态S_j的第m个混合成分的均值向量，U_jm是隐含状态S_j的第m个混合成分的协方差矩阵。C_jm满足随机约束：

$\begin{array}{c}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}=\mathrm{1,1}\&le;j\&le;N\\ C_{\mathrm{jm}}\&GreaterEqual;\mathrm{0,1}\&le;j\&le;N,1\&le;m\&le;M\end{array}---\left(10\right)$

C_jm、μ_jm、U_jm的重估是这样的：C_jm的重估值是HMM在隐含状态S_j带有第k个混合成分的次数的期望除以HMM处于隐含状态S_j的次数的期望。μ_jm的重估值是由第k个混合成分说明的观测矢量的部分期望值，U_jm是由第k个混合成分说明的观测矢量的部分协方差,见公式(11)～(13)。

$\begin{array}{c}{\stackrel{\&OverBar;}{C}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}\\ {\mathrm{\&gamma;}}_t(j,k)=\left[\frac{{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}\right]\left[\frac{C_{\mathrm{jk}}H[O_t,{\mathrm{\&mu;}}_{\mathrm{jk}},U_{\mathrm{jk}}]}{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}H[O_t,{\mathrm{\&mu;}}_{\mathrm{jm}},U_{\mathrm{jm}}]}\right]\end{array}---\left(11\right)$

${\stackrel{\&OverBar;}{\mathrm{\&mu;}}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)\&CenterDot;O_t}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(12\right)$

${\stackrel{\&OverBar;}{U}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)\&CenterDot;(O_t-{\mathrm{\&mu;}}_{\mathrm{jk}}){(O_t-{\mathrm{\&mu;}}_{\mathrm{jk}})}^{\&prime;}}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(13\right)$

π是初始状态概率矩阵，N是隐含状态的数目，M是每个隐含状态下的高斯混合成分的数目。

步骤二：在线诊断过程

a)在线采集工业过程的数据；

b)在原空间进行相似性分析，去除相似样本，得到数据集Ω₁，引入核函数φ将数据集Ω₁投影到高维特征空间，同样在高维特征空间进行数据的相似性分析，去除相似样本，得到数据集Ω₂，根据数据集Ω₂中数据计算相应的核矩阵以及核主元；

c)以核主元作为观测序列，选择最匹配的HMM模型，从而判断故障类型；

上述实施例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明做出的任何修改和改变，都落入本发明的保护范围。

Claims (6)

1.一种基于改进的KPCA和隐马尔科夫模型的工业故障诊断方法，其特征在于，该方法的步骤如下：

步骤一：离线建模，对工业过程采集的离线数据，在原空间进行相似性分析，去除相似样本，得到数据集Ω₁，引入核函数φ将数据集Ω₁投影到高维特征空间，同样在特征空间进行数据的相似性分析，去除相似样本，得到数据集Ω₂，根据数据集Ω₂中数据计算相应的核矩阵以及核主元，将得到的核主元作为隐马尔科夫模型(HMM)的观测序列，训练隐马尔科夫模型；

步骤二：在线诊断，对在线采集的数据在原空间进行相似性分析，去除相似样本，同样在特征空间进行数据的相似性分析，去除相似样本，计算相应的核矩阵以及核主元，即得到相应的观测序列后，选择最匹配的HMM模型，从而判断故障类型。

2.根据权利要求1所述的方法，其特征在于，步骤一所述的离线建模过程如下：

1)工业过程采集的第i个时刻监测数据其中n_i表示样本个数；

2)在原输入空间进行相似性分析，去除相似样本，得到数据集Ω₁，

2.1)初始时只有2个数据点，即新数据点根据以下提出的相似性指数公式进行相似性的判断，

$S_i=\frac{1}{2}{e}^{\left|\right|x_{\mathrm{new}}-{\tilde{x}}_i\left|\right|}+\frac{1}{2}{\left(\frac{{(x_{\mathrm{new}}-{\tilde{x}}_i)}^T(x_{\mathrm{new}}-{\tilde{x}}_{i-1})}{{\left|\right|x_{\mathrm{new}}-{\tilde{x}}_i\left|\right|}_2\left|\right|x_{\mathrm{new}}-{\tilde{x}}_{i-1}\left|\right|}\right)}^2---\left(1\right)$

其中x_new表示新数据点；

2.2)如果新数据点被引入，否则该样本点被舍弃，即N_i＝N_i-1，其中γ₀是预先设定好的一个足够小的值，满足

2.3)对原始输入空间进行相似性分析之后，被保留的数据点个数为n₁，所以我们获得了数据集Ω₁；

3)引入核函数φ将数据集Ω₁投影到高维特征空间，在特征空间进行相似性分析，去除相似样本，得到数据集Ω₂，

3.1)特征空间中相似性指数公式定义如下：

$S_f=\frac{1}{2}{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right)\left|\right|}_2+\frac{1}{2}{\left(\frac{{(\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right))}^T(\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{i-1}\right))}{{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_i\right)\left|\right|}_2{\left|\right|\mathrm{\&Phi;}\left(x_{\mathrm{new}}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{i-1}\right)\left|\right|}_2}\right)}^2---\left(2\right),$

3.2)如果新数据点被引入，否则该样本点被舍弃，其中γ₁是预先设定好的一个足够小的值，满足

3.3)对特征空间进行相似性分析之后，被保留的数据点个数为n₂，所以我们获得了数据集Ω₂；

4)对数据集Ω₂，计算核矩阵K以及核主元，

4.1)计算核矩阵

K_ij＝<Φ(x_i),Φ(x_j)>＝K(x_i,x_j)    (3)

其中使用径向基核函数 $K(x,y)=\exp (-\frac{{\left|\right|x-y\left|\right|}^2}{\mathrm{\&sigma;}}),\mathrm{\&sigma;}=\mathrm{rm},$ r为常数；

4.2)对核矩阵K进行中心化处理，

$\tilde{K}=k-1_{n_i}K-{K1}_{n_i}+1_{n_i}K{1}_{n_i}---\left(4\right)$

其中

4.3)计算主成分t_k

$t_k=<v_k,\mathrm{\&Phi;}\left(x\right)>=\underset{i=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^k<\mathrm{\&Phi;}\left(x_i\right),\mathrm{\&Phi;}\left(x\right)>=\underset{i=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^{k}K(x_i,x)---\left(5\right)$

5)将得到的核主元作为隐马尔科夫模型的观测序列O，训练隐马尔科夫模型，得到其参数λ＝(A,B,π,M,N)；

其中，A是隐含状态转移概率矩阵，描述了HMM模型中各个状态之间的转移概率，见下式(6)～(7)

$\begin{array}{c}{a}_{S_i,S_j}=a_{i,j}=P[q_{t+1}=S_j|q_t=S_i],1\&le;i,j\&le;N\\ \underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{i,j}=\mathrm{1,1}\&le;i\&le;N\end{array}---\left(7\right)$

式中，简记为a_i,j，表示在t时刻，状态为S_i的条件下，在t+1时刻状态是S_j的概率；

B是观测值概率密度矩阵，见下式

$\begin{array}{c}B=\{b_j\left(O\right),j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;N\}\\ b_j\left(O\right)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}H[O,{\mathrm{\&mu;}}_{\mathrm{jm}},U_{\mathrm{jm}}],1\&le;j\&le;N\end{array}---\left(9\right)$

其中O是观测向量，C_jm是隐含状态j的第m个混合成分的混合系数，H是对数凹或椭圆对称密度，μ_jm是隐含状态S_j的第m个混合成分的均值向量，U_jm是隐含状态S_j的第m个混合成分的协方差矩阵，C_jm满足随机约束：

$\begin{array}{c}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}=\mathrm{1,1}\&le;j\&le;N\\ C_{\mathrm{jm}}\&GreaterEqual;\mathrm{0,1}\&le;j\&le;N,1\&le;m\&le;M\end{array}---\left(10\right)$

C_jm、μ_jm、U_jm的重估是这样的：C_jm的重估值是HMM在隐含状态S_j带有第k个混合成分的次数的期望除以HMM处于隐含状态S_j的次数的期望，μ_jm的重估值是由第k个混合成分说明的观测矢量的部分期望值，U_jm是由第k个混合成分说明的观测矢量的部分协方差,见公式(11)～(13)，

$\begin{array}{c}{\stackrel{\&OverBar;}{C}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}\\ {\mathrm{\&gamma;}}_t(j,k)=\left[\frac{{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}\right][\frac{C_{\mathrm{jk}}H[O_t,{\mathrm{\&mu;}}_{\mathrm{jk}},U_{\mathrm{jk}}]}{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jm}}H[O_t,{\mathrm{\&mu;}}_{\mathrm{jm}},U_{\mathrm{jm}}]}\end{array}---\left(11\right)$

${\stackrel{\&OverBar;}{\mathrm{\&mu;}}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)\&CenterDot;O_t}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(12\right)$

${\stackrel{\&OverBar;}{U}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)\&CenterDot;(O_t-{\mathrm{\&mu;}}_{\mathrm{jk}}){(O_t-{\mathrm{\&mu;}}_{\mathrm{jk}})}^{\&prime;}}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(13\right)$

π是初始状态概率矩阵，N是隐含状态的数目，M是每个隐含状态下的高斯混合成分的数目。

3.根据权利要求2所述的方法，其特征在于，步骤二所述的在线诊断程如下：

a)在线采集工业过程的数据；

b)在原空间进行相似性分析，去除相似样本，得到数据集Ω₁，引入核函数φ将数据集Ω₁投影到高维特征空间，同样在高维特征空间进行数据的相似性分析，去除相似样本，得到数据集Ω₂，根据数据集Ω₂中数据计算相应的核矩阵以及核主元；

c)以核主元作为观测序列，选择最匹配的HMM模型，从而判断故障类型。

4.根据权利要求1所述的方法，其特征在于，所述的工业过程具有非线性。

5.根据权利要求1所述的方法，其特征在于，所述的工业故障为高炉冶炼过程故障。

6.一种根据权利要求1-4任一项所述的方法用于高炉冶炼过程故障诊断。