一种基于KPCA混合模型的多工况过程监控方法

摘要

本发明公开了一种基于KPCA混合模型的多工况过程监控方法，属于工业过程监控与诊断技术领域。本发明结合高斯混合模型和核主元分析模型，利用核主元分析在处理工业过程非线性、降低数据维度方面的优势以及高斯混合模型在处理非高斯、多工况等问题的优越性能，对工业过程进行监控。相比于现有的其他方法，由于充分考虑了工业过程的非高斯、非线性、多工况的复杂特性，本发明方法能更精确地估计各工况的统计特性，从而更准确及时地检测出多工况过程的各种故障。

说明书

技术领域

本发明属于工业过程监控与故障诊断领域，特别涉及一种基于KPCA混合模 型的多工况过程监控方法。

背景技术

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

对于过程监控和故障诊断问题，传统的方法大多采用多元统计过程监控技 术(Multivariable Statistical Process Monitoring,MSPM)，其中以主元分析 (Principal Component Analysis,PCA)和偏最小二乘(Partial Least Squares,PLS) 为代表等方法已在工业过程监控中得到了成功的应用。传统的MSPM方法均假 设过程数据服从高斯分布，变量之间是线性关系且数据来自单一的操作工况下， 但是实际中测量数据难以满足这些假设条件，常呈现非高斯、非线性和多工况 等特性。尽管，一些改进的方法，如针对非高斯ICA，核PCA(KPCA)等也被 提出。但是，当上述非高斯、非线性和多工况等特性同时存在的时，这些方法 仍然无法很好地解决。

近年来，基于PCA混合模型的多工况过程监控方法被提出用于解决上述问 题。将混合高斯模型和PCA相结合,用EM算法估计模型的工况数以及工况的 分布参数和主元数,并对每个子模型构建T²、SPE统计量实现时多工况过程的监 控。此种方法对每一个高斯成分模型建立了一个PCA模型，但是传统PCA只能 处理变量间的线性关系，并不能提取变量间的非线性信息，而在工业过程，变 量间的非线性关系是普遍存在的。由此，在工况存在强非线性情况下，可能会 引起故障检测发生错误。

发明内容

本发明的目的在针对现有技术的不足，提供一种基于KPCA混合模型的多 工况过程监控方法，利用核主元分析在处理工业过程非线性、降低数据维度方 面的优势以及高斯混合模型在处理非高斯、多工况等问题的优越性能，对工业 过程进行监控，从而更准确及时地检测出多工况过程的各种故障。

一种基于KPCA混合模型的多工况过程监控方法，该方法的步骤如下：

步骤一：离线建模，收集多工况过程正常运行的数据，构建高斯混合模型， 在高斯混合模型基础上，对每个高斯元空间进行KPCA变换和降维，建立了基 于KPCA的混合模型。采用EM算法来估计模型参数，并按照传统的KPCA方 法构建各工况的控制限；

步骤二：在线检测，采集在线运行数据，利用步骤一所述的混合模型将监 控样本按其后验概率大小分类到相应的工况中，并按照传统的KPCA方法计算 统计量，如果该统计量超出步骤一所建立的相对应的控制限，则判断故障发生。

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

1)利用工业过程采集的多工况监测数据构成X＝[x₁,x₂,…,x_n]^T∈R^n×m，其中 m表示监控变量的个数，n表示样本个数，x_i∈R^m,i＝1,…,n表示第i个样本；

2)它在有限GMM模型下的概率密度函数表示为其中 K表示GMM中混合高斯成分数目，w_i表示第i个单高斯成分的混合系数，且满 足和θ＝{θ₁,…,θ_K}分别表示局部和全局高斯模型参数集， 即均值向量μ_i和协方差矩阵Σ_i。相应的第i个分量的多元高斯密度函数可表示为 $g\left(x\right|\theta )=\frac{1}{{\left(2\pi \right)}^{m/2}{\left|{\Sigma }_i\right|}^{1/2}}\exp [-\frac{1}{2}{(x-{\mu }_i)}^T{\Sigma }_i^{-1}(x-{\mu }_i)].$

3)第i个单高斯成分样本子集表示为θ_i＝{μ_i,Σ_i}。 对数据集X_i进行KPCA投影。

3.1)引入核函数φ将数据集X_i投影到高维特征空间F，表示为

Φ:R^m→F        (1)

3.2)计算核矩阵K

K_ij＝〈Φ(x_i),Φ(x_j)〉＝K(x_i,x_j)      (2)

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

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

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

其中

3.4)计算主成分t_k

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

4)采用EM算法来估计模型参数

4.1)E-step

$p^{\left(s\right)}\left(C_i\right|x_j)=\frac{w_i^{\left(s\right)}g\left(x_i\right|{\mu }_i^{\left(s\right)},{\Sigma }_i^{\left(s\right)})}{\underset{k=1}{\overset{K}{\Sigma }}{w}_k^{\left(s\right)}g\left(x_j\right|{\mu }_k^{\left(s\right)},{\Sigma }_k^{\left(s\right)})}---\left(5\right)$

其中p^(s)(C_i|x_j)表示第s次迭代后第j个训练样本属于第i个高斯成分的后 验概率。

4.2)M-step

${\mu }_i^{(s+1)}=\frac{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)x_j}{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)}---\left(6\right)$

${\Sigma }_i^{(s+1)}=\frac{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)(x_j-{\mu }_i^{(s+1)}){(x_j-{\mu }_i^{(s+1)})}^T}{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)}---\left(7\right)$

$w_i^{(s+1)}=\frac{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)}{N}---\left(8\right)$

${\lambda }_{i,j}^{s+1}{\alpha }_{i,j}={\tilde{K}}_i^{(s+1)}{\alpha }_{i,j}---\left(9\right)$

其中，和α_i,j分别表示第(s+1)次迭代后，第i个高斯 成分的均值，协方差，先验概率和KPCA特征向量。

5)确定T²和SPE统计限

${T^2}_{\lim }=\frac{p(n_i-1)}{n_i(n_i-p)}{F}_{\alpha }(p,n_i-p)---\left(10\right)$

${\mathrm{SPE}}_{\lim }=g{\chi }_h^2---\left(11\right)$

步骤二所述的在线检测过程如下：

1)在线采集测量数据y∈R^m；

2)判断采集样本属于的工况类别。计算则x_t所属工 况为$c=\arg \underset{i=\mathrm{1,2},...,K}{\max }p(x_t\in C_i);$

3)计算采集样本的T²和SPE统计量。

T²＝[t₁,…,t_p]Λ^-1[t₁,…,t_p]^T        (12)

$\mathrm{SPE}={\left|\right|\phi \left(x\right)-{\phi }_p\left(x\right)\left|\right|}^2=\underset{i=1}{\overset{N}{\Sigma }}{t}_i^2-\underset{i=1}{\overset{P}{\Sigma }}{t}_i^2---\left(13\right)$

4)比较T²、SPE与式(10)、(11)中所建立的检测控制限T²_lim、SPE_lim之间的 大小，如果统计量超出控制限，则判断故障发生；如果统计量低于控制限，则 说明过程正常运行。

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

本发明有以下优势：

1.本发明首次提出一种基于KPCA混合模型的多工况过程监控方法，实现 对复杂过程的监测；

2.本发明能够解决过程数据存在的非高斯、非线性和多模态等问题，从而能 够更有效地监控。

具体实施方式

本发明提出的一种基于KPCA混合模型的多工况过程监控方法，包括以下 各步骤：

步骤一：离线建模

1)利用工业过程采集的多工况监测数据构成X＝[x₁,x₂,…,x_n]^T∈R^n×m，其中 m表示监控变量的个数，n表示样本个数，x_i∈R^m,i＝1,…,n表示第i个样本；

2)它在有限GMM模型下的概率密度函数表示为其中 K表示GMM中混合高斯成分数目，w_i表示第i个单高斯成分的混合系数，且满 足和θ＝{θ₁,…,θ_K}分别表示局部和全局高斯模型参数集， 即均值向量μ_i和协方差矩阵Σ_i。相应的第i个分量的多元高斯密度函数可表示为 $g\left(x\right|\theta )=\frac{1}{{\left(2\pi \right)}^{m/2}{\left|{\Sigma }_i\right|}^{1/2}}\exp [-\frac{1}{2}{(x-{\mu }_i)}^T{\Sigma }_i^{-1}(x-{\mu }_i)].$

3)第i个单高斯成分样本子集表示为θ_i＝{μ_i,Σ_i}。 对数据集X_i进行KPCA投影。

3.1)引入核函数φ将数据集X_i投影到高维特征空间F，表示为

Φ:R^m→F       (1)

3.2)计算核矩阵K

K_ij＝〈Φ(x_i),Φ(x_j)〉＝K(x_i,x_j)          (2)

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

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

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

其中

3.4)计算主成分t_k

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

4)采用EM算法来估计模型参数

4.1)E-step

$p^{\left(s\right)}\left(C_i\right|x_j)=\frac{w_i^{\left(s\right)}g\left(x_i\right|{\mu }_i^{\left(s\right)},{\Sigma }_i^{\left(s\right)})}{\underset{k=1}{\overset{K}{\Sigma }}{w}_k^{\left(s\right)}g\left(x_j\right|{\mu }_k^{\left(s\right)},{\Sigma }_k^{\left(s\right)})}---\left(5\right)$

其中p^(s)(C_i|x_j)表示第s次迭代后第j个训练样本属于第i个高斯成分的后 验概率。

4.2)M-step

${\mu }_i^{(s+1)}=\frac{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)x_j}{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)}---\left(6\right)$

${\Sigma }_i^{(s+1)}=\frac{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)(x_j-{\mu }_i^{(s+1)}){(x_j-{\mu }_i^{(s+1)})}^T}{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)}---\left(7\right)$

$w_i^{(s+1)}=\frac{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)}{N}---\left(8\right)$

${\lambda }_{i,j}^{s+1}{\alpha }_{i,j}={\tilde{K}}_i^{(s+1)}{\alpha }_{i,j}---\left(9\right)$

其中，和α_i,j分别表示第(s+1)次迭代后，第i个高斯 成分的均值，协方差，先验概率和KPCA特征向量。

5)确定T²和SPE统计限

${T^2}_{\lim }=\frac{p(n_i-1)}{n_i(n_i-p)}{F}_{\alpha }(p,n_i-p)---\left(10\right)$

${\mathrm{SPE}}_{\lim }=g{\chi }_h^2---\left(11\right)$

步骤二：在线监测

1)在线采集测量数据y∈R^m；

2)判断采集样本属于的工况类别。计算则x_t所 属工况为$c=\arg \underset{i=\mathrm{1,2},...,K}{\max }p(x_t\in C_i);$

3)计算采集样本的T²和SPE统计量。

T²＝[t₁,…,t_p]Λ^-1[t₁,…,t_p]^T           (12)

$\mathrm{SPE}={\left|\right|\phi \left(x\right)-{\phi }_p\left(x\right)\left|\right|}^2=\underset{i=1}{\overset{N}{\Sigma }}{t}_i^2-\underset{i=1}{\overset{P}{\Sigma }}{t}_i^2---\left(13\right)$

4)比较T²、SPE与式(10)、(11)中所建立的检测控制限T²_lim、SPE_lim之间的 大小，如果统计量超出控制限，则判断故障发生；如果统计量低于控制限，则 说明过程正常运行。

实施例

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

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

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

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

步骤一：离线建模

1)假设传感器采集高炉上料系统、喷吹系统、热风炉送风系统、回收煤气 与除尘、上料系统和渣水处理系统6种监测数据，且每种监测数据构成 X＝[x₁,x₂,…,x_n]^T∈R^n×m，其中m表示监控变量的个数，n表示样本个数 x_i∈R^m,i＝1,…,n表示第i个样本；

2)它在有限GMM模型下的概率密度函数表示为其中 K表示GMM中混合高斯成分数目，w_i表示第i个单高斯成分的混合系数，且满

足θ_i＝{μ_i,Σ_i}和θ＝{θ₁,…,θ_K}分别表示局部和全局高斯模型参数集， 即均值向量μ_i和协方差矩阵Σ_i。相应的第i个分量的多元高斯密度函数可表示为 $g\left(x\right|\theta )=\frac{1}{{\left(2\pi \right)}^{m/2}{\left|{\Sigma }_i\right|}^{1/2}}\exp [-\frac{1}{2}{(x-{\mu }_i)}^T{\Sigma }_i^{-1}(x-{\mu }_i)].$

3)第i个单高斯成分样本子集表示为θ_i＝{μ_i,Σ_i}。 对数据集X_i进行KPCA投影。

3.1)引入核函数φ将数据集X_i投影到高维特征空间F，表示为

Φ:R^m→F        (1)

3.2)计算核矩阵K

K_ij＝〈Φ(x_i),Φ(x_j)〉＝K(x_i,x_j)       (2)

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

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

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

其中

3.4)计算主成分t_k

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

4)采用EM算法来估计模型参数

4.1)E-step

$p^{\left(s\right)}\left(C_i\right|x_j)=\frac{w_i^{\left(s\right)}g\left(x_i\right|{\mu }_i^{\left(s\right)},{\Sigma }_i^{\left(s\right)})}{\underset{k=1}{\overset{K}{\Sigma }}{w}_k^{\left(s\right)}g\left(x_j\right|{\mu }_k^{\left(s\right)},{\Sigma }_k^{\left(s\right)})}---\left(5\right)$

其中p^(s)(C_i|x_j)表示第s次迭代后第j个训练样本属于第i个高斯成分的后 验概率。

4.2)M-step

${\mu }_i^{(s+1)}=\frac{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)x_j}{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)}---\left(6\right)$

${\Sigma }_i^{(s+1)}=\frac{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)(x_j-{\mu }_i^{(s+1)}){(x_j-{\mu }_i^{(s+1)})}^T}{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)}---\left(7\right)$

$w_i^{(s+1)}=\frac{\underset{j=1}{\overset{n_i}{\Sigma }}{p}^{\left(s\right)}\left(C_i\right|x_j)}{N}---\left(8\right)$

${\lambda }_{i,j}^{s+1}{\alpha }_{i,j}={\tilde{K}}_i^{(s+1)}{\alpha }_{i,j}---\left(9\right)$

其中，和α_i,j分别表示第(s+1)次迭代后，第i个高斯 成分的均值，协方差，先验概率和KPCA特征向量。

5)确定T²和SPE统计限

${T^2}_{\lim }=\frac{p(n_i-1)}{n_i(n_i-p)}{F}_{\alpha }(p,n_i-p)---\left(10\right)$

${\mathrm{SPE}}_{\lim }=g{\chi }_h^2---\left(11\right)$

步骤二：在线监测

1)在线采集测量数据y∈R^m；

2)判断采集样本属于的工况类别。计算则x_t所 属工况为$c=\arg \underset{i=\mathrm{1,2},...,K}{\max }p(x_t\in C_i);$

3)计算采集样本的T²和SPE统计量。

T²＝[t₁,…,t_p]Λ^-1[t₁,…,t_p]^T      (12)

$\mathrm{SPE}={\left|\right|\phi \left(x\right)-{\phi }_p\left(x\right)\left|\right|}^2=\underset{i=1}{\overset{N}{\Sigma }}{t}_i^2-\underset{i=1}{\overset{P}{\Sigma }}{t}_i^2---\left(13\right)$

4)比较T²、SPE与式(10)、(11)中所建立的检测控制限T²_lim、SPE_lim之间的 大小，如果统计量超出控制限，则判断故障发生；如果统计量低于控制限，则 说明过程正常运行。

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