CN108170648A - 一种基于已知数据回归的非高斯过程监测方法 - Google Patents

Abstract

本发明公开一种基于已知数据回归的非高斯过程监测方法，旨在将非高斯独立元成分转换为高斯分布的误差信息，以实现对非高斯过程对象实施精准的过程监测。具体来讲，本发明方法首先利用独立元分析(ICA)算法挖掘非高斯过程对象的独立元成分，然后在逐个假设变量数据缺失的条件下，利用已知数据回归(KDR)估计出独立元。最后，利用独立元估计误差实施基于主元分析模型的过程监测。一般而言，ICA算法能揭露出原始数据的本质，以ICA模型为基础通过KDR得到的估计误差是服从高斯分布的。而受益于误差的高斯分布特性，本发明方法所描述的正常区域更为精确，能显著提升传统ICA模型用于非高斯过程监测的故障检测能力。

Description

一种基于已知数据回归的非高斯过程监测方法

技术领域

本发明涉及一种数据驱动的过程监测方法，尤其涉及一种基于已知数据回归的非高斯过程监测方法。

背景技术

实施在线故障检测是保证生产安全与维持产品质量稳定的基本手段，对故障检测的研究伴随着整个生产工业的发展历程。当前的故障检测方法可以大致分为两类，其一是基于机理模型的故障检测方法，其二是基于数据的故障检测方法。基于机理模型的故障检测方法依赖于过程某些变量货参数的实际值与根据模型推理出的估计值之间的误差来实施故障检测。也就是说，如何生成误差是设计基于机理模型的故障检测方法的核心所在。基于数据的故障检测方法主要依赖于数据，不需要过程对象的机理模型，因此很适合于现代工业过程监测系统的实施与建立。与基于机理模型的故障检测方法生成误差的理念不同，基于数据的故障检测方法通过对过程数据进行挖掘，提取出潜藏的有用信息实施故障检测。

近年来，各种不同的数据挖掘算法都在故障检测领域找到了用武之地。其中，当以主元分析(Principal Component Analysis，PCA)与独立元分析(Independent ComponentAnalysis，ICA)两种算法最为常见。PCA与ICA算法用于故障检测的实施流程大同小异，最主要的不同之处在于两者提取潜藏成分的出发点不一样。具体来讲，PCA挖掘数据变量间的相关性特征，并使提取的潜藏成分(即主元)最大化的保留原始数据的方差信息。而ICA算法实在高阶统计量的指引下，挖掘出数据中潜藏的独立元成分。该独立元信息具备非高斯性，能更好地揭露原始数据的本质。由于现代工业过程采样数据一般不会满足高斯分布假设，因此ICA方法相比于PCA方法具有更广泛的应用性。而且，通过理论与实际的研究叶都发现，ICA方法能一般取得优越于PCA方法的故障检测效果。然而，ICA算法用于故障检测时通常是依赖距离型的统计量作为监测指标。例如，监测独立元变化情况常采用平方马氏距离，监测模型残差的变化一般使用平方欧式距离。而ICA提取的独立元成分是按照非高斯最大化原则来的，肯定不会服从高斯分布。从几何空间的角度来看，马氏距离与欧氏距离分别定义的是一个超椭球体与超球体。只有在数据服从高斯分布的前提查下，正常数据的可能变化范围能够完全填充超椭球或超球体内部空间。ICA方法中用马氏距离监测独立元的变化情况，在独立元不满足高斯分布时，该超椭球体内部空间有可能存在较多的“空洞”。某些故障样本数据转变成独立元后若刚好处于这些“空洞”位置，ICA是无法将其甄别出来的。由此可见，高斯分布对于PCA与ICA方法的重要性程度。

针对此问题，可行的解决方法是采用数据分布描述方法，如多变量的核密度估计法或支持向量描述。虽然，这两种方法能定义出一个精确描述正常波动范围的边界线。但是，这两类方法都需要提前设定好相应的模型参数。若模型参数设置不当，所定义的正常变化范围要么太紧凑导致较多的误报情况发生，要么太松散导致漏报率过高。此外，在没有充分先验知识的前提下，如何为这两类方法确定模型参数一直都是一个公开而未得到很好解决的问题。另一种解决思路，可以是在不丢失过程数据特征的前提下，通过某种方式将不明分布情况的数据转换成服从高斯分布的数据。现有文献中也存在几种高斯变换方法，但都很难直接应用于故障检测。若是按照基于机理模型误差生成的方式，通过数据模型生成某些变量的估计值，那么相应的估计误差一般而言是服从高斯分布的。这种方法可行但是其难点在于如何通过数据模型产生估计值。

发明内容

本发明所要解决的主要技术问题是：如何在ICA模型的基础上，将非高斯独立元成分转换为高斯分布的误差成分信息，以实现对非高斯过程对象实施精准的过程监测。本发明提供一种基于已知数据回归的非高斯过程监测方法，该方法通过逐个假设测量变量数据缺失，然后利用已知数据回归(KnownDataRegression，KDR)估计出相应的独立元成分，最后利用独立元估计误差实施过程监测。

本发明解决上述技术问题所采用的技术方案为：一种基于已知数据回归的非高斯过程监测方法，包括以下步骤：

(1)收集生产过程正常运行状态下的数据样本，组成训练数据集X∈R^n×m，并对每个变量进行标准化处理，得到均值为0，标准差为1的新数据矩阵其中，n为训练样本数，m为过程测量变量数，R为实数集，R^n×m表示n×m维的实数矩阵。

(2)利用ICA算法为建立相应的ICA模型：初始化变量下标号i＝1，为d个独立成分列向量组成的矩阵，W∈R^m×d为分离矩阵，A∈R^m×d为混合矩阵，E∈R^n×m表示模型误差，上标号T表示矩阵或向量的转置。利用ICA算法为建立ICA模型的具体实施过程如下所示：

①计算的协方差矩阵其中C∈R^m×m；

②计算矩阵C的所有特征值和特征向量，并剔除小于0.0001的特征值及其对应的特征向量，得到特征向量矩阵P＝[p₁，p₂，…，P_M]∈R^m×M以及特征值对角矩阵D＝diag(λ₁，λ₂，…，λ_M)∈R^M×M；

值得注意的是，这里求解得到的特征向量p₁，p₂，…，p_M都必须是单位长度的向量。

③根据公式对进行白化处理，得到Z∈R^n×M，并初始化i＝1；

④取列向量c_i为M×M维单位矩阵中的第i列，

⑤按照如下所示公式更新c_i，即：

c_i←E{Zg(c_i ^TZ)}-E{h(c_i ^TZ)}c_i (1)

上式(1)中，E{}表示求取期望值(即向量的平均值)，函数g和h的具体形式如下所示：

g(u)＝tanh(u) (2)

h(u)＝[sech(u)]² (3)

上式(2)与(3)中，u为函数自变量，在这里指代c_i ^TZ中的元素。

⑥对更新后的向量c_i依次按照下式进行正交标准化处理：

c_i←c_i/||c_i|| (5)

⑦重复步骤⑤～⑥直至向量c_i收敛，并保存向量c_i；

⑧判断i＜M？若是，置i＝i+1后，重复步骤④-⑧；若否，执行步骤⑨；

⑨将得到的所有M个向量c₁，c₂，…，c_M组成矩阵C＝[c₁，c₂，…，c_M]∈R^M×M，并按照如下所示公式计算分离矩阵W₀∈R^m×M与混合矩阵A₀∈R^m×M：

A₀＝PD^1/2C (6)

W₀＝PD-^1/2C (7)

⑩计算A₀中每一列向量的长度，分别记为L₁，L₂，…，L_M，并将L₁，L₂，…，L_M按照数值大小进行降序排列得到l₁，l₂，…，l_M，那么保留的独立成分个数d为满足下列条件的最小值：

将A₀中列向量长度最大的d个列向量组成新的混合矩阵A∈R^m×d，同时从W₀中取出与A对应的列向量组成新分离矩阵W∈R^m×d；

最后得到的ICA模型为

(3)假设矩阵中第i列数据缺失，为不失一般性，可将矩阵表述成其中，为假设缺失的数据(实为矩阵中第i列)，由矩阵中剩余的列组成，为已知数据。

(4)利用最小二乘的思路构建已知数据与独立元成分矩阵S之间的回归模型，即：

上式中，回归矩阵E_i∈R^n×d为独立元估计误差矩阵。

值得指出的是，独立元估计误差E_i的秩rank(E_i)＝1，也就是说E_i中存在较多冗余信息。关于证明rank(E_i)＝1的具体思路如下：

估计误差矩阵E_i可以按照如下公式进行推算：

那么，秩rank(E_i)的运算就满足如下所示公式：

由于实为矩阵中第i列，则又因为误差矩阵E_i一般不为零矩阵，因此有rank(E_i)＝1。

(5)对估计误差实施奇异值分解，即：

E_i＝U_iΛ_iV_i ^T (12)

其中，U_i与V_i为酉矩阵，对角矩阵Λ_i实际上只包含了一个非零奇异值，这是因为rank(E_i)＝1。因此，从误差E_i到向量U_i之间的变换矩阵为Θ_i＝V_iΛ_i ^-1。

(6)根据公式U_i＝E_iΘ_i计算出剔除冗余信息后的误差向量U_i，并判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(3)；若否，则将得到的误差向量组成矩阵U＝[U_i，U₂，…，U_m]后继续执行下一步骤。

(7)利用PCA算法为包含独立元估计误差的矩阵U建立相应的PCA故障检测模型，保留模型参数集其中H∈R^m×k为载荷矩阵，A∈R^k×k为对角矩阵，与Q_c分别为监测统计量的控制上限。具体的实施过程如下所示：

①计算U的协方差矩阵Z＝U^TU/(n-1)；

②求解Z所有特征值γ₁≥γ₂≥…≥γ_m所对应的特征向量h₁，h₂…，h_m；

③设置保留的主成分个数k为满足如下所示条件的最小值，并将对应的k个特征向量组成载荷矩阵H＝[h₁，h₂…，h_k]；

④得到对角矩阵A＝diag{γ₁，γ₂，…，γ_k}；

⑤根据如下所示公式分别确定监测统计量T²与Q对应的控制上限T_c ²与Q_c：

上两式中，置信水平α＝99％，F_α(k，n-k)表示自由度为k与n-k的F分布，表示权重为g＝v/2b，自由度为h＝2a²/b的χ²分布，a与b分别是U对应的Q统计量的估计均值和估计方差。

(8)收集新采样时刻的数据样本x∈R^1×m，对其实施与步骤(1)中相同的标准化处理得到新数据向量后，初始化i＝1。

(9)假设向量中第i个数据缺失，同理，可表示成其中，x_i ^#为第i个缺失的数据，由向量中除缺失数据以外的元素组成。

(10)利用如下所示公式计算出对应于向量在缺失第i个数据的前提下相应的独立元估计误差e_i：

上式中，独立元实际值独立元估计值

(11)利用公式u_i＝e_iΘ_i计算消除冗余信息后的误差u_i后，判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(9)；若否，则将得到的误差组成向量u＝[u₁，u₂，…，u_m]并继续执行下一步骤。

(12)调用PCA故障检测模型参数集Φ，并根据如下所示公式计算统计监测指标T²与Q的具体数值：

T²＝uHA^-1H^Tu^T (17)

Q＝u(I-HH^T)u^T (18)

判断T²与Q的具体数值是否大于对应控制上限与Q_c？若否，则当前样本为正常工况采样；若是，则当前采样数据有可能来自故障工况，理应继续监测接下来的3～6个新样本，若都超限，则说明当前工况已出现故障，若都没超限，则说明当前工况仍旧处于正常状态。

与传统方法相比，本发明方法的优势在于：

本发明方法在传统ICA模型的基础上，通过逐个假设缺失数据并利用KDR推算出相应的独立元计估计值，从而将非高斯分布的独立元成分转换成高斯分布的估计误差，并以误差作为监测对象实施基于PCA的在线故障监测。一般而言，ICA算法能揭露出原始数据的本质，以ICA模型为基础通过KDR得到的估计误差通常是服从高斯分布的。而受益于误差的高斯分布特性，本发明方法所描述的正常区域更为精确，能显著提升传统ICA模型用于非高斯过程监测的故障检测能力，是一种更为优选的非高斯过程监测方法。

附图说明

图1为本发明方法的基本原理示意图。

图2为原始数据的非高斯分布情况散点图。

图3为独立元估计误差的高斯分布检验图。

具体实施方式

下面结合附图与具体的实施案例对本发明方法进行详细的说明。

如图1所示，本发明公开一种基于已知数据回归的非高斯过程监测方法。下面结合一个具体的工业过程的例子来说明本发明方法的具体实施过程，以及相对于现有方法的优越性。

按照如下所示公式设计一个5个测量变量的非高斯过程：

其中，源信号s＝[s₁，s₂]是根据如下所示公式生成的：

上两式(19)与(20)中，测量噪声v服从均值为0，标准差为0.2的高斯分布，t₁与t₂都为在区间[0，1]上均匀分布的随机数。根据上两个公式模拟仿真出1000个采样样本，以便按照传统ICA、ICA-PCA、和本发明方法建立相应的故障检测模型。可以很明显的发现，这1000个训练数据采样于非高斯过程对象，相应的散点图如图2所示。

步骤1)：对训练数据进行标准化处理，得到新数据矩阵

步骤2)：对建立ICA模型并初始化变量下标号i＝1。

步骤3)：假设矩阵中第i列数据缺失，为不失一般性，可将矩阵表述成其中，为假设缺失的数据(实为矩阵中第i列)，由矩阵中剩余的列组成，为已知数据；

步骤4)：利用最小二乘的思路构建已知数据与独立元成分矩阵S之间的回归模型，保留回归矩阵

步骤5)：对估计误差实施奇异值分解，得到误差E_i到向量U_i之间的变换矩阵为Θ_i＝V_iΛ_i ^-1。

步骤6)：根据公式U_i＝E_iΘ_i计算出剔除冗余信息后的误差向量U_i，并判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(3)；若否，则将得到的误差向量组成矩阵U＝[U₁，U₂，…，U₅]后继续执行下一步骤。

为了展示本发明方法得到的误差是服从高斯分布的，特将U₁，U₂，…，U₅的高斯分布检验结果显示于图3中。显而易见的是，高斯分布检验结果都近乎成一条直线，相应的数据是服从高斯分布的。

步骤7)：利用PCA算法为包含独立元估计误差的矩阵U建立相应的PCA故障检测模型，保留模型参数集

根据表1中所列的2中故障工况，分别对应仿真出相应的测试数据集，每类测试数据集各包含1000个样本，且故障工况于201个采样点时引入。相应的在线故障监测实施过程如下所示：

步骤8)：收集新采样时刻的数据样本x∈R^1×m，对其实施与步骤(1)中相同的标准化处理得到新数据向量∈R^1×m后，初始化i＝1。

步骤9)：假设向量中第i个数据缺失，同理，可表示成

步骤10)：利用公式计算出对应于向量在缺失第i个数据的前提下相应的独立元估计误差e_i。

步骤11)：利用公式u_i＝e_iΘ_i计算消除冗余信息后的误差u_i后，判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(9)；若否，则将得到的误差组成向量u＝[u₁，u₂，…，u_m]并继续执行下一步骤。

步骤12)：调用PCA故障检测模型参数集Φ，并计算统计监测指标T²与Q的具体数值，实施在线故障检测。

针对表1中的两类故障工况，三种不同方法所取得故障检测率同样列于表1中，取得最大检测率的数值已经用粗体标出。本发明方法针对这两类故障的监测结果都是最优的，而且过程监测性能的提升很显著。

表1：非高斯过程故障工况及其相应的故障检测结果(％)。

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

Claims (3)

1.一种基于已知数据回归的非高斯过程监测方法，其特征在于，包括以下步骤：

离线建模阶段的实施过程如下所示：

(1)收集生产过程正常运行状态下的数据样本，组成训练数据集X∈R^n×m，并对每个变量进行标准化处理，得到均值为0，标准差为1的新数据矩阵其中，n为训练样本数，m为过程测量变量数，R为实数集，R^n×m表示n×m维的实数矩阵；

(2)利用独立元分析(ICA)算法为建立相应的ICA模型：初始化变量下标号i＝1，为d个独立成分列向量组成的矩阵，W∈R^m×d为分离矩阵，A∈R^m×d为混合矩阵，E∈R^n×m表示模型误差，上标号T表示矩阵或向量的转置；

(3)假设矩阵中第i列数据缺失，为不失一般性，可将矩阵表述成其中，为假设缺失的数据(实为矩阵中第i列)，由矩阵中剩余的列组成，为已知数据；

(4)利用最小二乘的思路构建已知数据与独立元成分矩阵S之间的回归模型，即：

上式中，回归矩阵E_i∈R^n×d为独立元估计误差矩阵；

(5)对估计误差实施奇异值分解，即：

E_i＝U_iΛ_iV_i ^T (2)

其中，U_i与V_i为酉矩阵，对角矩阵Λ_i实际上只包含了一个非零奇异值，这是因为rank(E_i)＝1；因此，从误差E_i到向量U_i之间的变换矩阵为Θ_i＝V_iΛ_i ^-1；

(6)根据公式U_i＝E_iΘ_i计算出剔除冗余信息后的误差向量U_i，并判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(3)；若否，则将得到的误差向量组成矩阵U＝[U₁，U₂，…，U_m]后继续执行下一步骤；

(7)利用PCA算法为包含独立元估计误差的矩阵U建立相应的PCA故障检测模型，保留模型参数集其中H∈R^m×k为载荷矩阵，A∈R^k×k为对角矩阵，与Q_c分别为监测统计量的控制上限；

在线故障检测的实施过程如下所示：

(8)收集新采样时刻的数据样本x∈R^1×m，对其实施与步骤(1)中相同的标准化处理得到新数据向量后，初始化i＝1；

(9)假设向量中第i个数据缺失，同理，可表示成其中，x_i ^#为第i个缺失的数据，x_i ^*由向量中除缺失数据以外的元素组成；

(10)利用如下所示公式计算出向量在缺失第i个数据的前提下相应的独立元估计误差e_i：

上式中，独立元实际值独立元估计值

(11)利用公式u_i＝e_iΘ_i计算消除冗余信息后的误差u_i后，判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(9)；若否，则将得到的误差组成向量u＝[u₁，u₂，…，u_m]并继续执行下一步骤；

(12)调用PCA故障检测模型参数集Φ，并根据如下所示公式计算统计监测指标T²与Q的具体数值：

T²＝uHA^-1H^Tu^T (4)

Q＝u(I-HH^T)u^T (5)

判断T²与Q的具体数值是否大于对应控制上限与Q_c？若否，则当前样本为正常工况采样；若是，则当前采样数据有可能来自故障工况，理应继续监测接下来的3～6个新样本，若都超限，则说明当前工况已出现故障，若都没超限，则说明当前工况仍旧处于正常状态。

2.根据权利要求1所述的一种基于已知数据回归的非高斯过程监测方法，其特征在于，所述步骤(2)中为建立相应的ICA模型的具体实施过程如下所示：

①计算的协方差矩阵其中C∈R^m×m；

②计算矩阵C的所有特征值和特征向量，并剔除小于0.0001的特征值及其对应的特征向量，得到特征向量矩阵P＝[p₁，p₂，…，p_M]∈R^m×M以及特征值对角矩阵D＝(λ₁，λ₂，…，λ_M)∈R^M×M；值得注意的是，这里求解得到的特征向量p₁，p₂，…，p_M都必须是单位长度的向量；

③根据公式对进行白化处理，得到Z∈R^n×M，并初始化i＝1；

④取列向量c_i为M×M维单位矩阵中的第i列，

⑤按照如下所示公式更新c_i，即：

c_i←E{Zg(c_i ^TZ)}-E{h(c_i ^TZ)}c_i (6)

上式(6)中，E{ }表示求取期望值(即向量的平均值)，函数g和h的具体形式如下所示：

g(u)＝tanh(u) (7)

h(u)＝[sech(u)]² (8)

上式(7)与(8)中，u为函数自变量，在这里指代c_i ^TZ中的元素；

⑥对更新后的向量c_i依次按照下式进行正交标准化处理：

c_i←c_i/||c_i|| (10)

⑦重复步骤⑤～⑥直至向量c_i收敛，并保存向量c_i；

⑧判断i＜M？若是，置i＝i+1后，重复步骤④～⑧；若否，执行步骤⑨；

⑨将得到的所有M个向量c₁，c₂，…，c_M组成矩阵C＝[c₁，c₂，…，c_M]∈R^M×M，并按照如下所示公式计算分离矩阵W₀∈R^m×M与混合矩阵A₀∈R^m×M：

A₀＝PD^1/2C (11)

W₀＝PD^-1/2C (12)

⑩计算A₀中每一列向量的长度，分别记为L₁，L₂，…，L_M，并将L₁，L₂，…，L_M按照数值大小进行降序排列得到l₁，l₂，…，l_M，那么保留的独立成分个数d为满足下列条件的最小值：

将A₀中列向量长度最大的d个列向量组成新的混合矩阵A∈R^m×d，同时从W₀中取出与A对应的列向量组成新分离矩阵W∈R^m×d；

最后得到的ICA模型为

3.根据权利要求1所述的一种基于已知数据回归的非高斯过程监测方法，其特征在于，所述步骤(7)中为误差U建立相应的PCA故障检测模型的具体实施过程如下所示：

①计算U的协方差矩阵Ξ＝U^TU/(n-1)；

②求解Ξ所有特征值γ₁≥γ₂≥…≥γ_m所对应的特征向量h₁，h₂…，h_m；

③设置保留的主成分个数k为满足如下所示条件的最小值，并将对应的k个特征向量组成载荷矩阵H＝[h₁，h₂…，h_k]；

④得到对角矩阵A＝diag{γ₁，γ₂，…，γ_k}；

⑤根据如下所示公式分别确定监测统计量T²与Q对应的控制上限：

上两式中，置信水平α＝99％，F_α(k，n-k)表示自由度为k与n-k的F分布，表示权重为g＝v/2b，自由度为h＝2a²/b的χ²分布，a与b分别是U对应的Q统计量的估计均值和估计方差。