CN116339275A - 基于全结构动态自回归隐变量模型的多尺度过程故障检测方法 - Google Patents

Abstract

本发明公开了一种基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，在线收集过程的多尺度数据样本，利用多尺度全结构自回归动态隐变量模型得到当前时刻测试样本的动态T²统计量、每个测量尺度下的静态T²统计量和SPE统计量，并与预先获得的检测控制限对比，得出过程的在线检测结果。该方法在完全利用过程完整多尺度数据信息的同时，又能利用卡尔曼滤波和贝叶斯方法分别考虑数据的动态和静态特性，并对动态和静态潜隐变量实现精确的估计，通过降维得到的动态和静态潜隐变量从不同的数据子空间对故障产生响应。在故障检测精度和应用范围上实现提升。

Description

基于全结构动态自回归隐变量模型的多尺度过程故障检测 方法

技术领域

本发明涉及一种控制方法，具体是涉及一种基于全结构动态自回归隐变量模型的多尺度过程故障检测方法。

背景技术

随着现代流程工业的发展，工业生产规模和复杂程度逐渐增加，及时对大规模工业过程可能产生的故障进行检测受到广泛的重视。随着集散控制系统(DCS)在工业领域的广泛应用，大量生产过程变量以较高的测量尺度被保存在系统数据库中，一些调度层的中间变量和需要离线化验的关键质量变量的测量尺度则较低。这导致在复杂工业过程数据存在多尺度特性。随着多元统计过程监测技术的发展，海量过程数据通过降维、重构以及可视化实现对过程的实时监测，并广泛应用于化工、有色冶金、治污等领域。传统的静态过程监测模型如主成分分析(PCA)和偏最小二乘(PLS)在面向高维时序相关过程数据时监测效果不佳，而多数动态过程监测模型如动态PCA(DPCA)、规范变量分析(CVA)以及线性高斯状态空间模型(LGSSM)等模型描述了过程的时序相关特性但并未考虑到过程数据的多尺度特性。多尺度动态潜隐变量模型能够完全利用过程数据的多尺度信息，并利用期望最大化(EM)算法对模型参数进行估计。然而上述动态模型假设动态潜隐变量同样包含了过程的静态特性，因此在面对强耦合的复杂工业过程故障检测时效果受到影响。因此需要提出一种既能利用过程多尺度数据，同时又能分别考虑到过程动态特性和静态特性的工业过程故障检测方法。

发明内容

本发明的目的在于针对现有技术的不足，提供一种基于全结构动态自回归隐变量模型的多尺度工业过程故障检测方法。

一种基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，包括：在线收集化工过程的多尺度数据样本，待测试样本集，对测试样本集进行标准化，利用预先构建的多尺度全结构自回归动态隐变量模型(TMsARDLV模型)得到当前时刻测试样本的动态T²统计量、每个测量尺度下的静态T²统计量和SPE统计量，并与预先获得的检测控制限对比，得出过程的在线检测结果；所述多尺度全结构自回归动态隐变量模型中，多尺度数据样本与动态潜隐变量和各个测量尺度下的静态潜隐变量呈线性关系。

所述多尺度全结构自回归动态隐变量模型训练时，收集对应过程(比如造纸废水处理过程)正常运行状态下多种不同测量尺度的变量样本，组成建模用的训练样本集，对训练样本集进行标准化，然后用于所述多尺度全结构自回归动态隐变量模型的构建。

所述不同测量尺度是指采用不同的采样率进行采样。本发明中，采集的多个变量样本对应的采样率可以完全不同，也可以是部分相同。

作为一种实施方式，所述多尺度全结构自回归动态隐变量模型结构如下：

其中：x(k)表示k时刻模型的动态潜隐变量；z(k-1)包含了过去L个时刻该模型的动态潜隐变量，L为滞后时间；A为模型动态潜隐变量的状态转移矩阵；v(k)为第k时刻模型的动态噪声；y^ξ(k)为k时刻采集得到的变量；C^ξ(k)为第k时刻变量和动态潜隐变量之间的动态发散矩阵；t^ξ(k)为第k时刻的静态潜隐变量；Ψ^ξ(k)为第k时刻变量和静态潜隐变量之间的静态发散矩阵；w^ξ(k)为第k时刻变量的测量噪声；ξ为表示当前时刻变量被采集到的样本属性。

为便于处理多尺度数据以及后续公式的直观表达，模型训练和实际使用时引入采样系数λ，其表现形式如下：

所述多尺度全结构自回归动态隐变量模型结构中：y^ξ(k)、C^ξ(k)、t^ξ(k)、Ψ^ξ(k)、w^ξ(k)均取决于采样系数λ，具体讲：

y^ξ(k)的值源于{y₁,y₂,...y_m,...y_M}，并由采样系数决定其构成，y_m表示第m个测量尺度下的变量样本集；

C^ξ(k)的值源自{C₁,C₂,...,C_m,...,C_M}，并由采样系数决定其构成，C_m为第m个测量尺度下变量和动态潜隐变量之间的动态发散矩阵；

t^ξ(k)的值源于{t₁,t₂,...t_m,...t_M}，并由采样系数决定其构成，t_m为第m个尺度下的变量独有的静态潜隐变量；

Ψ^ξ(k)的值源于{Ψ₁,Ψ₂,...Ψ_m,...Ψ_M}，并由采样系数决定其构成，Ψ_m为第m个尺度的静态发散矩阵；

w^ξ(k)的值源于{w₁,w₂,...w_m,...w_M}，并由采样系数决定其构成，w_m为采用第m种测量尺度的变量的测量噪声，其服从高斯分布w_m～N(0,R_m)。

某时刻y^ξ(k)＝C^ξ(k)x(k)+Ψ^ξ(k)t^ξ(k)+w^ξ(k)的具体表现形式，由该时刻的采样系数决定。假设k时刻，E种测量尺度下的变量被采集，则该E种测量尺度下样本对应的采样系数λ＝1，E种测量尺度分别第m₁,m₂,...,m_i...,m_E测量尺度，来源于M个采样尺度，并由采样系数决定，此时，模型中变量的y^ξ(k)＝C^ξ(k)x(k)+Ψ^ξ(k)t^ξ(k)+w^ξ(k)表达式如下：

其中：

为k时刻采集到的第m_i个测量尺度对应的变量；/>

为k时刻采集到的第m_i个测量尺度下变量和动态潜隐变量之间的动态发散矩阵；/>

为k时刻采集到的第m_i个尺度下的变量独有的静态潜隐变量；/>

为k时刻采集到的第m_i个尺度的静态发散矩阵；/>

为k时刻采集到的第m_i个尺度的变量的测量噪声。

可以选择利用期望最大化算法来更新模型参数。即所述多尺度全结构自回归动态隐变量模型采用期望最大化算法优化得到。

进一步讲，采用期望最大化算法优化过程中(或利用期望最大化算法来更新模型参数时)，在E步利用卡尔曼滤波算法结合当前模型参数估计动态潜隐变量的后验概率，利用贝叶斯方法结合当前模型参数估计静态潜隐变量的后验概率；在M步中，则通过极大化似然函数的方式更新多尺度全结构动态自回归隐变量模型参数；最后，反复迭代E步和M步直至达到模型收敛条件。

可以选择首先利用多尺度全结构动态自回归隐变量模型训练后的参数得到训练样本的动态隐变量的后验期望和不同测量尺度下的静态潜隐变量的后验期望，利用动态潜隐变量的期望值和方差构建

统计量(动态T²统计量)，利用静态潜隐变量的期望值和方差构建T_s ²统计量(静态T²统计量)。基于模型的重构误差构建SPE统计量。

的控制限/>

由χ²分布得到，估计方法为：

其中，D为动态隐变量维度。

T_s ²的控制限

估计方法为：/>

即T_s ²服从χ²分布，其中

gh＝mean(T_s ²)

2g²h＝var(T_s ²)

其中mean()为求均值，var()为求方差，

代表卡方分布，g和h分别为卡方分布的系数和自由度；通过上式可以求得g和h，进而求得T_s ²统计量的控制限。

SPE的控制限SPE_lim估计方法为：

即SPE服从χ²分布，其中：

gh＝mean(SPE)

2g²h＝var(SPE)

其中mean()为求均值，var()为求方差，

代表卡方分布，g和h分别为卡方分布的系数和自由度；通过上式可以求得g和h，进而求得SPE统计量的控制限。

作为进一步优选的方案，一种基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，包括：

(I)针对多尺度过程，利用多尺度训练样本集进行训练得到多尺度全结构动态自回归隐变量模型以及检测控制限；

(II)在线收集新的过程中与训练样本集过程变量和关键质量变量对应的多尺度过程样本数据，得到测试样本集；

(III)对得到的测试样本集进行相同的标准化；

(IV)针对标准化后的测试样本集，利用得到的多尺度全结构动态自回归隐变量模型得到当前时刻测试样本的

(动态T²统计量)、T_s ²(各测量尺度下的静态T²统计量)和SPE统计量，通过与步骤(4)得到的检测控制限对比，得出过程的在线检测结果。

作为更具体的方案，一种基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，包括：

(1)采集过程中正常运行状态下的不同测量尺度的各种变量样本，组成建模用的训练样本集；

(2)对得到的训练样本集进行标准化，使得各个标准化后的变量值与潜隐变量间存在线性相关关系；

(3)基于预处理后的训练样本集构建多尺度全结构动态自回归隐变量模型；

(4)根据建立的多尺度全结构动态自回归隐变量模型，得到训练样本的

T_s ²和SPE统计量相应的检测控制限；

(5)在线收集新的过程中与训练样本集过程变量和关键质量变量对应的多尺度过程样本数据，得到测试样本集；

(6)对得到的测试样本集进行步骤(2)的标准化；

(7)采用步骤(4)得到的多尺度全结构动态自回归隐变量模型得到当前时刻测试样本的

T_s ²和SPE统计量，通过与步骤(4)得到的检测控制限对比，得出过程的在线检测结果。

建模和监测用的变量样本包括但不限于以下数据中的一种或多种，1)可以利用集散控制系统检测得到的设备层面变量(比如表征设备状态的电流、电压、功率、位移、转速等等，采样频率通常是毫秒级到秒级)、过程层面变量(比如表征过程运行状态的温度、压力、流量、液位、pH值等等，采样频率通常是分钟到小时级)，这类变量又称过程变量；(2)可以利用化验手段获取正常工作状态下比较难测的质量层数据(包括表征产品产量和质量变量，比如目标物质浓度等，采样频率通常为小时级到天级)，这种数据又被称为关键质量变量；(3)以及需要进行统计计算的指标层数据(表示运行经济指标、能耗指标等指标的数据，采样频率可能是周级到月级)。

所述过程变量样本采用集散控制系统收集；所述关键质量变量采用化验手段收集。本发明中，所述的化验手段包括但是不限于化学滴定、试纸检测、纯度检测(比如借助HPLC、LC-Ms等等进行的检测)、核磁检测等等测试方法。

本发明中，过程变量样本一般指那些可以由现有的传感器检测得到的变量，可以方便地通过集散控制系统收集，比如可以是温度、压力、流量等。所述的关键质量变量一般不能、很难或者是不宜采用现有传感器直接检测检测的变量，比如某一中间体或者原料的浓度等。作为一种实施方式本发明主要采用程变量样本或者程变量样本和关键质量变量进行模型的构建等。

假设工业过程共有M种不同的测量尺度，共收集到K个历史数据样本用于模型训练。设定采样时间段内，采集到的M种不同测量尺度的正常变量样本集Y，Y＝{Y₁；Y₂；...；Y_m；...；Y_M}；其中，第m种采样率的变量样本为Y_m；M种不同测量尺度样本的样本数量分别为K₁,K₂,...,K_M：

其中，R表示实数集；在M种多尺度过程变量和质量变量数据中，G_m为采用第m种采样率的过程变量或质量变量的变量数，K_m为第m种采样率下的每个过程变量质量变量的样本数；y_m(i)为第m种采样率下过程变量或质量变量对应的第i个样本。将这些数据存入历史数据库，组成建模用的训练样本集模型。

训练过程中以及实际获得在线数据后，都需要进行所述的标准化操作，经过所述标准化，使得每个过程变量或者每个关键质量变量中的每个元素值围绕0上下波动，大于0表明高于平均水平，小于0表明低于平均水平，且与潜隐变量间存在线性相关关系。进一步讲，所述标准化方法为：在某一测量尺度下，针对该尺度下的每个过程变量或者每个关键质量变量，先使得其中的每个元素减去其对应的每个过程变量或者每个关键质量变量的平均值，然后再除以该变量对应样本集的整体标准差。

建模过程中和实际检测中，对于数据的标准化方法相同。

在建模过程中，可以首先构建动态多尺度模型；然后将所构建的动态多尺度模型，拓展为全结构动态多尺度模型；最后基于得到的全结构动态多尺度模型，利用预处理后的训练样本集构建多尺度全结构动态自回归隐变量模型。

具体讲，由于经过标准化后的训练数据集与潜隐变量间存在线性相关关系，可得到如下的多尺度动态隐变量模型：

其中，x∈R^D为维度为D的该模型动态潜隐变量，x(k)表示k时刻该模型的动态潜隐变量，z(k-1)＝[x^T(k-1)x^T(k-2)...x^T(k-L)]^T∈R^DL×1包含了过去L个时刻该模型的动态潜隐变量，L为滞后时间。A∈R^D×DL为模型的动态潜隐变量的状态转移矩阵。v(k)为第k时刻模型的动态噪声且服从高斯分布v(k)～N(0,Q)。训练样本集包含了总共M个测量尺度，第m个测量尺度下变量和动态潜隐变量之间的线性关系由动态发散矩阵C_m表示，所有动态发散矩阵集为{C₁,C₂,...C_m,...C_M}；同理，第m个测量尺度下的过程变量由y_m表示且所有过程变量的合集为{y₁,y₂,...y_m,...y_M}；第m个测量尺度下测量噪声为u_m，其服从高斯分布u_m～N(0,Ω_m)，其中Ω_m为第m个测量尺度下测量噪声的方差，噪声方差集为{Ω₁,Ω₂,...Ω_m,...Ω_M}。

C^ξ(k)为模型在第k时刻的动态发散矩阵，具体可以表述为：

其中，E表示了在k时刻数据集Y所含数据对应的测量尺度总数；C_k1为第一个被采集到的过程变量数据或关键质量变量对应的采样率下的模型参数。因此C^ξ(k)的值源于{C₁,C₂,...C_m,...C_M}且取决于该时刻其对应的样本是否被采集到。同理y^ξ(k)源于{y₁,y₂,...y_m,...y_M}。u^ξ(k)为模型第k时刻变量的输入噪声矩阵，其服从如下的高斯分布：

同理，此处的Ω^ξ(k)源于{Ω₁,Ω₂,...Ω_m,...Ω_M}。

同时，基于以上所构建的动态多尺度隐变量模型，通过将每个测量尺度下的噪声矩阵u^ξ(k)进一步分解，可以得到如下多尺度全结构自回归动态隐变量模型：

其中，

为采用第m个测量尺度变量独有的静态潜隐变量，S_m是第m个测量尺度下静态潜隐变量的维度，且所有尺度静态潜隐变量的合集为{t₁,t₂,...t_m,...t_M}；

表示该模型在k时刻多尺度样本对应的静态潜隐变量，t^ξ(k)源于{t₁,t₂,...t_m,...t_M}且取决于该时刻其是否被采集。同理，Ψ_m为采用第m种测量尺度的静态发散矩阵且其合集为{Ψ₁,Ψ₂,...Ψ_m,...Ψ_M}；/>

源于{Ψ₁,Ψ₂,...Ψ_m,...Ψ_M}。w_m为采用第m种测量尺度的变量的测量噪声，其服从高斯分布w_m～N(0,R_m)，R_m为测量噪声的方差且所有测量尺度方差合集为{R₁,R₂,...R_m,...R_M}，同理：

w^ξ(k)～N(0,R^ξ(k))

R^ξ(k)＝diag{R_k1；R_k2；...；R_kE}

其中R^ξ(k)源于{R₁,R₂,...R_m,...R_M}。

在线收集新的造纸废水处理过程中与训练样本集过程变量对应的多尺度过程样本数据，采样时间间隔由采样率最高的过程变量的采样率确定(可以等于或者大于最高的采样率，当然也可以小于)，得到测试样本集，并进行与训练阶段相同的标准化。

在得到测试样本后，首先，采用得到的多尺度全结构线性动态自回归隐变量模型以及模型训练得到模型参数对测试样本进行处理，利用卡尔曼滤波算法计算测试样本的当前时刻的动态潜隐变量的期望值，利用贝叶斯方法计算测试样本当前时刻静态潜隐变量的期望值。利用动态潜隐变量和静态潜隐变量的期望值计算测试样本的

T_s ²和SPE统计量。判断其是否超过控制限，得出化工生产过程的在线故障检测结果。

本发明中所述检测控制限，动态T²统计量检测控制限由卡方分布得到。各测量尺度下的静态T²统计量检测控制限和SPE统计量检测控制限由训练样本所得对应统计量获得。

可以使用本发明的过程包括但不限于化工废水处理过程(比如造纸废水处理过程)、化工原料或者中间体制备过程、化工原料或者中间体的精制过程、农药/医药/医药中间体的制备过程、精制过程等。

本发明基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，以过程不同测量尺度下的过程变量和质量变量为建模样本，充分考虑了多尺度数据的自相关性和不同测量尺度数据内部的互相关性，分别提取全局动态潜隐变量和各尺度下局部静态潜隐变量，并在此模型的基础上建立了相应的故障检测方法。以解决复杂多尺度过程故障检测难题。该方法在完全利用过程完整多尺度数据信息的同时，又能利用卡尔曼滤波和贝叶斯方法分别考虑数据的动态和静态特性，并对动态和静态潜隐变量实现精确的估计，通过降维得到的动态和静态潜隐变量从不同的数据子空间对故障产生响应。在故障检测精度和应用范围上实现提升。

与现有技术相比，本发明的有益效果体现在：

本发明的基于全结构动态隐变量模型的多尺度造纸废水处理过程的故障检测方法，以过程不同测量尺度的过程变量以及质量变量作为建模样本，在充分利用完整多尺度数据的同时，分别考虑了数据的自相关性和互相关性并设计全局动态隐变量和每个尺度独有的静态隐变量，模型参数的估计通过期望最大化(EM)算法和卡尔曼滤波算法以及贝叶斯方法实现，并在此模型的基础上建立了相应的故障检测方法，以解决多尺度复杂耦合过程故障检测的难题。该方法在实现多采样率信息处理的同时，既可以完整利用数据信息，又能利用卡尔曼滤波充分考虑数据的动态特性，并对动态潜隐变量实现精确的估计，使降维后的少数动态潜隐变量实现对难以直接测量的多采样关键质量变量更好的估计与描述，从而在故障检测精度与应用范围上实现提升。

附图说明

图1R2S过程TMsARDLV对数似然函数随模型迭代训练次数的变化。

图2TMsARDLV模型R2S过程故障1检测结果：(a)动态T²和SPE统计量故障检测结果(b)静态T²统计量故障检测结果。

图3TMsARDLV模型R2S过程故障2检测结果：(a)动态T²和SPE统计量故障检测结果(b)静态T²统计量故障检测结果。

图4MLGSSM模型R2S过程故障1检测结果。

图5MLGSSM模型R2S过程故障2检测结果。

图6MFA模型R2S过程故障1检测结果。

图7MFA模型R2S过程故障2检测结果。

图8MR-PCA模型R2S过程故障1检测结果。

图9MR-PCA模型R2S过程故障2检测结果。

图10FA₁模型R2S过程故障1检测结果。

图11FA₁模型R2S过程故障2检测结果。

图12FA₂模型R2S过程故障1检测结果。

图13FA₂模型R2S过程故障2检测结果。

具体实施方式

以造纸废水处理过程为例，对本发明做进一步说明：

一种基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，该方法针对造纸废水处理过程中故障检测问题，首先收集正常工作正常工作状态下的各种变量，这些变量包括但不限于如下变量：(1)可以利用集散控制系统检测得到的设备层面变量(比如表征设备状态的电流、电压、功率、位移、转速等等，采样频率通常是毫秒级到秒级)、过程层面变量(比如表征过程运行状态的温度、压力、流量、液位等等，采样频率通常是分钟级到小时级(过程变量)；(2)可以利用化验手段获取正常工作状态下比较难测的质量层数据(包括表征产品产量和质量变量，采样频率通常为小时级到天级)(关键质量变量)；(3)以及需要进行统计计算的指标层数据(表示运行经济指标、能耗指标等指标的数据，采样频率可能是周级到月级)，利用这些数据作为模型训练集数据，建立多尺度全结构自回归动态隐变量模型。该模型结构由期望最大化算法估计得到。同时获取故障检测统计量控制限。在此基础上，对在线造纸废水处理过程进行采样，获得多尺度测试样本，然后可利用已有的模型结构估计测试样本的动态潜隐变量以及静态潜隐变量，构建故障检测统计量，并根据所得统计量是否在对应控制限内，对造纸废水处理过程中可能发生的故障进行检测，实现最终的故障检测效果。

本发明是一种基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，包括以下步骤：

第一步：利用集散控制系统收集造纸废水处理过程中正常运行的不同测量尺度的过程变量数据，同时利用化验手段获取比较难测的关键质量变量的多尺度正常数据，还可以根据需要，获取设备层面变量数据或/和指标层变量数据，并组成建模用的训练样本集，假设设定时间段内，采集到了M种不同尺度的变量样本集Y，Y＝{Y₁；Y₂；...；Y_m；...；Y_M}；其中，第m种测量尺度(采样频率)的变量样本为Y_m；M种不同尺度样本的样本数量分别为K₁,K₂,...,K_M：

其中，R表示实数集；在M种多尺度变量数据中，G_m为采用第m种尺度下变量的变量个数，K_m为第m种测量尺度下的每个变量采集的样本数；y_m(i)为第m种测量尺度下对应的第i个样本。将这些数据存入历史数据库，组成建模用的训练样本集。

第二步：对样本集Y进行标准化，即每个样本集中的元素先减去其对应的变量样本的平均值，然后再除以其对应变量样本集的整体标准差，使得各个标准化后的元素值围绕0上下波动，大于0表明高于平均水平，小于0表明低于平均水平，且与潜隐变量间存在线性相关关系，可得到如下的多尺度动态隐变量模型：

训练样本集总共包含M个测量尺度，

表示k时刻采集得到的变量，E表示该时刻被采集到的变量对应的测量尺度类型的个数，第m个测量尺度下的变量样本集由y_m表示；同理，/>

表示k时刻对应的变量和动态潜隐变量之间的动态发散矩阵，C_m为第m个测量尺度下变量和动态潜隐变量之间的动态发散矩阵；

表示k时刻对应的测量噪声，第m个测量尺度下测量噪声为u_m，其服从高斯分布u_m～N(0,Ω_m),m＝1,2,...,M，其中Ω_m为第m个测量尺度下测量噪声的方差。

为便于处理多尺度数据以及后续公式的直观表达，此处引入一种采样系数λ，其表现形式如下所示：

因此，C^ξ(k)的值都源自{C₁,C₂,...,C_M}并由k时刻对应的采样系数λ决定其构成，也就是说取决于该测量尺度下的变量是否被采集。同理，u^ξ(k)的值源于{u₁,u₂,...,u_M}并由k时刻对应的采样系数λ决定其构成。例如，在一个三采样率的系统中，M为3，在k时刻，假设仅第1种测量尺度的正常过程变量样本值被采集，则λ_k1＝1，λ_k2＝0，λ_k3＝0，C^ξ(k)＝C₁(k),u^ξ(k)＝u₁(k)，模型在k时刻结构为：

而当第1种和第3种测量尺度的正常变量样本值被采集时，第2种测量尺度对应的正常变量样本值未被采集，λ_k1＝λ_k3＝1且λ_k2＝0，C^ξ(k)＝[C₁(k)；C₃(k)],u^ξ(k)＝[u₁(k)；u₃(k)]，模型在k时刻结构为：

x∈R^D为维度为D的该模型动态潜隐变量，x(k)表示k时刻该模型的动态潜隐变量，z(k-1)＝[x^T(k-1)x^T(k-2)...x^T(k-L)]^T∈R^DL×1包含了过去L个时刻该模型的动态潜隐变量，L为滞后时间。其初值z(L)＝[x^T(L)x^T(L-1)...x^T(1)]^T且z(L)～N(μ_L,V_L)，μ_L为z(L)的期望，V_L为z(L)的方差；A∈R^D×DL为模型为动态潜隐变量的状态转移矩阵。v(k)为第k时刻模型的动态噪声，其方差为Q，具体分布为v(k)～N(0,Q)；

基于以上所构建的多尺度动态隐变量模型，考虑到系统数据的自相关性和互相关性，对模型进行拓展，可以得到如下多尺度全结构动态自回归隐变量模型：

此模型中，x(k)、A、z(k-1)、v(k)、y^ξ(k)、C^ξ(k)定义同上。所有测量尺度静态潜隐变量的合集为{t₁,t₂,...t_m,...t_M}，其中

为第m个尺度下的变量独有的静态潜隐变量，S_m为第m个测量尺度下静态潜隐变量的维度；/>

表示该模型在k时刻多尺度样本对应的静态潜隐变量，E表示k时刻所采集的数据对应的测量尺度类型的个数，t^ξ(k)源于{t₁,t₂,...t_m,...t_M}且取决于该时刻其对应的变量{Y₁；Y₂；...；Y_m；...；Y_M}是否被采集。同样引入采样系数λ，静态发散矩阵的合集为{Ψ₁,Ψ₂,...Ψ_m,...Ψ_M}，其中Ψ_m为第m个尺度的静态发散矩阵；/>

源于{Ψ₁,Ψ₂,...Ψ_m,...Ψ_M}且取决于采样系数λ。w^ξ(k)为过程变量的测量噪声，/>

w^ξ(k)同样来源于{w₁,w₂,...w_m,...w_M}，其中w_m为采用第m种测量尺度的变量的测量噪声，其服从高斯分布w_m～N(0,R_m)，R_m为测量噪声的方差，且所有测量尺度方差合集为{R₁,R₂,...R_m,...R_M}，同理：

w^ξ(k)～N(0,R^ξ(k))

R^ξ(k)＝diag{R₁；R₂；...；R_E}

其中R^ξ(k)源于{R₁,R₂,...R_m,...R_M}且取决于采样系数λ。C^ξ(k)同理。

第三步：利用期望最大化(EM)算法来更新模型参数，在E步利用卡尔曼滤波算法和贝叶斯方法结合当前模型参数分别估计动态潜隐变量和静态潜隐变量的后验概率(后验期望)；在M步中，则通过极大化似然函数的方式更新多尺度全结构动态自回归隐变量模型参数。最后，反复迭代E步和M步直至达到模型收敛条件。

首先，对模型参数{C_m,Ψ_m,R_m，A，Q，μ_L，V_L}，(m＝1,2,...,M)随机进行初始化；由于卡尔曼滤波的输入数据格式的需要，对模型做出适当变换，其变换结果如下：

其中，k时刻的扩展的动态潜隐变量为z(k)，其具体表达式为z(k)＝[x^T(k)x^T(k-1)...x^T(k-L)]^T∈R^DL×1，包含了过去L个值，L为滞后时间；

为变换后的动态矩阵，/>

表示变换后的k时刻的动态噪声，其方差为/>

具体分布为/>

模型参数的具体变换如下所示：

其中，I为单位阵；加粗的0为零矩阵；

模型完整对数极大似然估计公式如下：

其中，Θ表示极大似然函数值，constant表示任意常数。p(·)表示概率密度的计算，p(x(k)|z(k-1))表示求解x(k)关于z(k-1)的概率密度函数，p(t₁(k))表示第一个测量尺度下的静态隐变量的概率密度函数，以此类推。p(y^ξ(k)|x(k),t^ξ(k))表示求解在k时刻对应的y^ξ(k)关于{x(k),t^ξ(k)}的概率密度函数。

在模型参数估计的E步，根据当前的模型参数的初始值，利用卡尔曼滤波算法获取模型动态潜隐变量估计的更新值，其主要公式为：

其中，

表示利用k-1时刻的预测结果对k时刻训练样本的扩展的动态潜隐变量进行估计，/>

表示k-1时刻的训练样本的最优估计；/>

表示k时刻训练样本的扩展的动态潜隐变量的最优估计；/>

表示/>

对应的协方差，/>

表

对应的协方差；/>

表示/>

对应的协方差；K_k表示k时刻的卡尔曼增益。为了方便表示，将/>

表示为/>

将/>

表示为/>

同样在E步，根据当前模型初始值，利用贝叶斯方法对每个尺度下的静态隐变量后验期望进行估计，首先需要得到第m个测量尺度下变量在k时刻观测变量和根据动态潜隐变量得到的重构值之间的残差：

根据得到的残差计算得到第m个尺度下的静态潜隐变量后验期望，主要公式如下：

E(t_m(k)t_m(k)^T)＝(M_m)^-1+E(t_m(k))E(t_m(k))^T

其中，E(·)表示期望的计算，因此E(t_m(k))表示第m个尺度下独有的静态隐变量t_m(k)的后验期望，E(t_m(k)t_m(k)^T)表示E(t_m(k))对应的协方差。其中(M_m)^-1＝[Ψ_mΨ_m ^T+R_m ^-1]^-1；

比较新模型参数所对应的极大似然值Θ_new与其原模型参数对应的极大似然值Θ_old的差异，如果||Θ_new-Θ_old||²＜ε，则进入到第四步，否则，继续进行EM算法的迭代，其中，ε为模型收敛的阈值。

在M步，根据E步的结果，对极大似然函数求取关于模型参数的偏导，获取模型参数{C_m,Ψ_m,R_m，A，Q，μ_L，V_L}，(m＝1,2,...,M)的更新值

分别如下：

/>

其中，E(·)表示期望的计算，

表示动态隐变量的后验期望，/>

可以通过/>

得到，表示为：/>

表示

对应的协方差。/>

表示/>

和/>

的协方差/>

表示

对应的协方差。

K表示这个统一的样本空间的样本总数，K_m表示第m种测量尺度的实际样本数。

第四步，在模型参数估计完成收敛之后，根据EM算法步骤得到的训练集动态潜隐变量后验期望

根据训练集不同动态潜隐变量的所有样本值计算得到每个变量的协方差并构建协方差矩阵/>

接着计算得到数据的动态T²统计量记为/>

计算方式如下：

后期实际监测过程中，采用上述公式，实现对在先样本的动态T²统计量的计算；

统计量的控制限由/>

分布估计，表示为：

其中D为动态潜隐变量的维度。

根据EM算法步骤得到训练集各个尺度下的静态潜隐变量后验期望E(t_m(k)),m＝1,2...M；k＝1,2,...,K_m，根据训练集不同静态隐变量的所有样本值计算得到每个变量的协方差并构建协方差矩阵

接着计算得到数据在每个测量尺度下的静态T²统计量记为/>

计算方式如下：

T_s ²的控制限

估计方法为：/>

即T_s ²服从χ²分布，其中

gh＝mean(T_s ²)

2g²h＝var(T_s ²)

其中mean()为求均值，var()为求方差，

代表卡方分布，g和h分别为卡方分布的系数和自由度；通过上式可以求得g和h，进而求得T_s ²统计量的控制限。/>

通过模型对训练集样本的重构误差可以构建SPE统计量。首先计算模型对样本的重构值和方差以及重构误差e(k)：

因此SPE统计量可以计算得到，表示为：

SPE(k)＝e^T(k)Φ(k)^-1e(k)

SPE的控制限SPE_lim估计方法为：

即SPE服从χ²分布，其中：

gh＝mean(SPE)

2g²h＝var(SPE)

其中mean()为求均值，var()为求方差，

代表卡方分布，g和h分别为卡方分布的系数和自由度；通过上式可以求得g和h，进而求得SPE统计量的控制限。

第五步：在线收集新的造纸废水处理过程中与训练样本集过程变量对应的多采样率过程样本数据，采样时间间隔由采样率最高的变量的采样率确定(可以等于或者大于最高的采样率，当然也可以小于，本实施例中，选择等于)，得到测试样本集，并进行标准化，该步骤的标准化与第二步可以采用相同的标准化方法，其样本数量分别为K′₁,K′₂,...,K′_m：

第六步：采用多尺度全结构动态自回归隐变量模型以及模型训练得到模型参数{C_m,Ψ_m,R_m，A，Q，μ_L，V_L}，(m＝1,2,...,M)对测试样本进行处理，利用卡尔曼滤波算法计算测试样本的第k时刻的动态潜隐变量的期望值

其中，

表示k-1时刻的测试样本的最优预测结果；

表示利用k-1时刻的预测结果对k时刻测试样本的扩展的动态潜隐变量进行估计；/>

表示k时刻测试样本的扩展的动态潜隐变量的最优估计；/>

表示

对应的协方差；/>

表示/>

对应的协方差；/>

表示

对应的协方差；K_test(k)表示k时刻测试样本的卡尔曼增益；为了后续简化表示，

表示为/>

表示为/>

利用在k时刻卡尔曼滤波算法得到的测试样本的动态潜隐变量的期望值/>

计算出k时刻样本的/>

统计量：/>

其中

为训练集动态隐变量的协方差。

利用贝叶斯方法计算测试样本在第m个测量尺度下的第k时刻的动态潜隐变量的期望值t_m,test(k)：

M_m ^-1＝[Ψ_mΨ_m ^T+R_m ^-1]^-1

利用第m个测量尺度下的第k时刻的动态潜隐变量的期望值

计算出k时刻在第m个测量尺度下样本的T_s ²统计量：

其中cov(t_m)为训练集每个尺度下静态隐变量对应的协方差。

基于模型重构误差，构建SPE统计量，重构值和方差表示为：

基于此可以计算得到模型的重构误差以及SPE统计量：

第七步，判断其是否超过第四步得到的控制限，得出造纸废水处理过程的在线故障检测结果：如果超过所述的控制限，则判断为发生故障；否则判断生产过程正常。

应用例：

我们利用造纸污水处理厂真实的R2S厌氧反应器数据，进一步测试了所提出的TMsARDLV模型的故障检测结果。造纸废水处理过程是一个复杂过程。在此过程中，污水处理厂将废水转化为另一种废水，这种废水可以在对环境影响最小的情况下返回水循环，或直接重复利用。污水处理厂通常包括混合池、澄清池、调节容器、厌氧反应器、循环竖管、缺氧池、沼气池等。造纸废水处理系统首先将杂质转移到非水相，并去除大的悬浮物、沉淀物和浮油。然后通过微生物厌氧好氧反应对废水进行生物降解，最后通过二沉池、砂滤、超滤、反渗透循环利用废水。通过添加适合于改善水质的微生物产品，显著降低污泥产量和化学需氧量浓度。由于操作条件的复杂性和原料变化的多样性，多种故障时有发生。

除此之外，造纸废水过程同样是一个多尺度过程，不同变量的测量尺度不同。其中，过程变量有两个测量尺度，分别是一小时和两小时。而和质量相关的变量以及输出变量由于只在每天进行一次实验室化验分析得到数据，因此其测量尺度为每24小时检测一次。表1给出了R2S过程所有变量描述。由于输出型变量在本次实验并未使用，因此没有给出。

表1R2S过程所有输入型变量

实验结果和分析

在本节实验中，总共选取了两个R2S厌氧反应器中18个典型的过程和质量变量进行多尺度过程建模和监测。变量的详细说明见表2。

表2实验选取变量名

可以看出，过程有三种测量尺度，M＝3。其中变量(1-3)每小时采样一次，变量(4-11)每2小时被测量一次。这两类变量均采用在线硬件传感器采集。此外，变量(12-18)都是与质量相关的变量，每天在实验室检测并记录一次，即其测量尺度为24h/次。利用变量(1-18)在78天产生的数据构建了上文提出的TMsARDLV。整个过程共收集了1872个变量(1-3)样本，936个变量(4-11)样本以及78个变量(12-18)样本。

在TMsARDLV中，模型动态主成分维度选取为9，每个尺度下稳态隐变量维度设置为3，利用的历史数据长度L设置为3，最大迭代步数设置为200。TMsARDLV模型每次迭代似然函数如图1所示。可以看到在200步之前似然函数不断增大且逐渐收敛。接下来将对实验设置进行详细阐述。

为了和TMsARDLV模型进行对比，本节选取了多采样率线性高斯状态空间模型(MLGSSM)以及多采样率因子分析模型(MR-FA)、多采样率主成分分析模型(MR-PCA)、因子分析模型1(FA₁)、因子分析模型1(FA₂)进行比较。其中，FA₁为使用降采样方法的因子分析模型(FA)，在方法建模过程中所有样本都被降到最低测量尺度；FA₂则是只使用最快测量尺度下的样本来建模的FA模型。MLGSSM模型隐层动态主元个数同样设置为9；而MR-FA和MR-PCA主元个数设置为3；FA₁、FA₂的主元个数分别设置为3和2。上述主成分个数的选择都是基于交叉验证法。为了进行故障检测，本实验对从实际运行过程中采集的两种故障进行了测试。经过适当的预处理，提取到两类故障数据。这两类故障都包含了54天的测试数据，其中故障发生在第20天。在显著性水平为0.99的情况下，表3列出了不同模型对造纸废水故障数据的误报率和漏检率，其中故障0为测试正常数据，对应结果为其误报率。

R2S过程数据的实验结果验证了该方法的有效性和先进性，在取得较低故障误报率的情况下，TMsARDLV方法取得了卓越的结果，优于大多数对比模型。模型的故障检测结果如图2～13所示。

根据实验结果可以看出，本发明提出的TMsARDLV模型在造纸废水处理过程故障检测任务中展现了良好的结果。在保证低误报率的情况下，相比于对比方法，基于TMsARDLV模型的故障检测方法的故障漏检率最低，对于大多数复杂故障都能做到及时响应。相比于其他对比方法，所提出故障检测方法针对造纸废水处理过程中的故障特点设计了针对不同隐变量空间的故障检测指标，对于复杂耦合的过程能够从不同层面实现对故障的检测。此外，所提出的方法为高阶动态模型，对于高度时序相关的过程数据有更好的鲁棒性和更好的故障检测准确性。

表3展示了造纸废水处理过程数据实验结果。其中，F0表示当在线监测数据为正常样本时，各个模型将正常样本误检测为故障的比率，称之为误报率。由于模型结构的不同，不同的模型存在不同数量的故障检测指标，但是都基本分为T²和SPE两类故障检测指标。从表中可以看出，其中MLGSSM、FA₁、FA₂都只设计了单个T²和SPE统计量。而MR-FA和MR-PCA由于是多采样率(多尺度)模型，模型在每个尺度下的残差空间分别设计了SPE统计量，因此这两个模型有一个T²统计量和多个SPE统计量。而本发明提出的TMsARDLV模型针对动态隐空间和静态隐空间分别设计了T²统计量，因此专利方法具有两个T²统计量和一个SPE统计量用于过程故障检测。F1和F2分别表示当在线监测数据为故障1和故障2样本时，不同模型的故障检测漏检率，称为漏报率。针对两类不同故障的故障检测结果可以评估不同模型的故障检测表现。

表3-1造纸废水处理过程数据实验结果1

表3-2造纸废水处理过程数据实验结果2

表3-3造纸废水处理过程数据实验结果3

从表中可以看出，提出的专利方法TMsARDLV在不同类型故障下的不同类统计指标中均取得最优的结果。故障1中T²统计量比次最优方法(MLGSSM)的漏检率下降了50％。SPE统计量的漏检率也相应有所下降，对于故障2，专利方法在两类故障检测指标上的表现相比于对比方法都取得了故障检测效果上的大幅进步。同一故障同一指标下的最优结果用红色加粗标记。

Claims (10)

1.一种基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，其特征在于，包括：在线收集化工过程的多尺度数据样本，得测试样本集，对测试样本集进行标准化，利用多尺度全结构自回归动态隐变量模型得到当前时刻测试样本的动态T²统计量、每个测量尺度下的静态T²统计量和SPE统计量，并与预先获得的检测控制限对比，得出化工过程的在线检测结果；所述多尺度全结构自回归动态隐变量模型中，多尺度数据样本与动态潜隐变量和各个测量尺度下的静态潜隐变量呈线性关系。

2.根据权利要求1所述的基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，其特征在于，包括：

(I)针对多尺度过程，利用获得的多尺度训练样本集进行模型训练，得到多尺度全结构动态自回归隐变量模型以及动态T²统计量、每个测量尺度下的静态T²统计量和SPE统计量的检测控制限；

(II)在线收集新的多尺度过程中与训练样本集过程变量和关键质量变量对应的多尺度过程样本数据，得到测试样本集；

(III)对得到的测试样本集进行相同的标准化；

(IV)针对标准化后的测试样本集，利用得到的多尺度全结构动态自回归隐变量模型得到当前时刻测试样本的动态T²统计量、各测量尺度下的静态T²统计量和SPE统计量，通过与得到的检测控制限对比，得出多尺度过程的在线检测结果。

3.根据权利要求1所述的基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，其特征在于，收集对应化工过程正常运行状态下多种不同测量尺度的变量，组成建模用的训练样本集，对训练样本集进行标准化，然后用于所述多尺度全结构自回归动态隐变量模型的构建。

4.根据权利要求1所述的基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，其特征在于，所述多尺度全结构自回归动态隐变量模型结构如下：

其中：x(k)表示k时刻模型的动态潜隐变量；z(k-1)包含了过去L个时刻该模型的动态潜隐变量，L为滞后时间；A为模型动态潜隐变量的状态转移矩阵；v(k)为第k时刻模型的动态噪声；y^ξ(k)为k时刻采集得到的变量；C^ξ(k)为第k时刻采集到的变量和动态潜隐变量之间的动态发散矩阵；t^ξ(k)为第k时刻的采集到的变量对应的静态潜隐变量，Ψ^ξ(k)为第k时刻采集到的变量和对应静态潜隐变量之间的静态发散矩阵；w^ξ(k)为第k时刻采集到的变量的测量噪声；ξ为表示当前时刻变量被采集到的样本属性。

5.根据权利要求4所述的基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，其特征在于，模型训练和实际使用时引入采样系数，其表现形式如下：

所述多尺度全结构自回归动态隐变量模型结构中：

y^ξ(k)的值源于{y₁,y₂,...y_m,...y_M}，并由采样系数决定其构成，y_m表示第m个测量尺度下的变量样本集；

C^ξ(k)的值源自{C₁,C₂,...,C_m,...,C_M}，并由采样系数决定其构成，C_m为第m个测量尺度下变量和动态潜隐变量之间的动态发散矩阵；

t^ξ(k)的值源于{t₁,t₂,...t_m,...t_M}，并由采样系数决定其构成，t_m为第m个尺度下的变量独有的静态潜隐变量；

Ψ^ξ(k)的值源于{Ψ₁,Ψ₂,...Ψ_m,...Ψ_M}，并由采样系数决定其构成，Ψ_m为第m个尺度的静态发散矩阵；

w^ξ(k)的值源于{w₁,w₂,...w_m,...w_M}，并由采样系数决定其构成，w_m为采用第m种测量尺度的变量的测量噪声，其服从高斯分布w_m～N(0,R_m)。

6.根据权利要求1所述的基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，其特征在于，多尺度全结构自回归动态隐变量模型采用期望最大化算法优化得到。

7.根据权利要求6所述的基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，其特征在于，采用期望最大化算法优化过程中，在E步利用卡尔曼滤波算法和贝叶斯方法结合当前模型参数分别估计动态潜隐变量和静态潜隐变量的后验概率；在M步中，则通过极大化似然函数的方式更新多尺度全结构动态自回归隐变量模型参数；反复迭代E步和M步直至达到模型收敛条件。

8.根据权利要求1所述的基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，其特征在于，所述检测控制限由卡方分布直接得到或者利用训练样本所得对应统计量通过卡方分布获得，或者采用上述两种方法的结合获得。

9.根据权利要求8所述的基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，其特征在于，所述检测控制限分别由如下方法得到：

动态T²统计量的控制限

根据动态潜隐变量的维度由/>

分布估计得到；

每个测量尺度下的静态T²统计量的控制限

由训练集的每个测量尺度下的静态T²统计量利用χ²分布得到，其中每个测量尺度下的静态T²统计量由各个尺度下的静态潜隐变量后验期望、以及训练集不同静态隐变量的协方差矩阵计算得到；

SPE的控制限SPE_lim：由训练集的SPE统计量经过χ²分布得到，其中SPE统计量由训练集样本的重构误差得到。

10.根据权利要求1所述的基于全结构动态自回归隐变量模型的多尺度过程故障检测方法，其特征在于，所述多尺度过程为造纸废水处理过程。

Images