CN109669412B - 一种基于新型动态独立成分分析的非高斯过程监测方法 - Google Patents

Abstract

本发明公开一种基于新型动态独立成分分析的非高斯过程监测方法，旨在结合可处理自相关动态性数据的动态内部主成分分析模型与可处理非高斯数据的独立成分分析模型的优势。具体来讲，本发明方法首先利用动态内部主成分分析算法分别提取出自相关的动态特征成分与交叉相关的静态特征成分。其次，在对特征成分进行白化处理后，利用合并后白化特征成分作为初始独立成分并迭代求取动态独立成分变量模型。最后，基于此种动态模型实施动态非高斯过程监测。可以说，本发明方法利用了动态内部主成分分析算法分开提取动态成分与静态成分的能力，再进一步结合能提取非高斯特征成分的独立成分分析算法。因此，本发明方法是一种可行的动态非高斯过程监测方法。

Description

一种基于新型动态独立成分分析的非高斯过程监测方法

技术领域

本发明涉及一种数据驱动的过程监测方法，尤其涉及一种基于新型动态独立成分分析的非高斯过程监测方法。

背景技术

随着先进测量技术与计算技术的飞速发展，现代工业过程逐步走向数字化管理，工业“大数据”热潮正是在这个背景下产生的。由于生产过程对象可以离线存储与在线测量海量的数据，这些数据蕴含着能体现生产过程运行状态的信息，利用采样数据实施过程运行状态的监测于是乎得到了较多学者们的青睐。事实上，无论是学术界还是工业界，都投入了大量的人力与物力研究以故障检测与诊断为核心任务的过程监测方法。在数据驱动的过程监测研究领域，统计过程监测是被研究得最多的方法，其中当以主成分分析(PrincipalComponent Analysis，PCA)与独立成分分析(Independent Component Analysis，ICA)为最主流的实施技术手段。一般而言，由于ICA算法能够挖掘出数据中潜藏的非高斯成分信息，更能揭示对象的本质，因此更适合于非高斯过程对象的监测。

建立数据驱动的过程监测模型时，需要充分挖掘正常工况下采样数据的特征。由于采样时间间隔较短，采样数据不可避免的存在时间序列上的自相关性。因此，数据的自相关性这种动态特征是必须考虑的一个问题。针对动态过程监测问题的研究，最常见的思路就是使用增广矩阵，将数据的自相关性与交叉相关性混淆在一起实施特征提取。此外，也有学者在在PCA算法优化目标的基础上，对目标函数进行改进从而有效地挖掘出数据的自相关性特征。这种算法被称为动态内部主成分分析(Dynamic Inner Principal ComponentAnalysis，DiPCA)算法，在自相关特征挖掘之后通过自回归模型过滤掉自相关性特征，进而使用传统PCA算法建模。DiPCA算法虽能挖掘自相关性特征，但是未能保留自相关特征，同时也未能更进一步的挖掘出自相关的非高斯独立成分信息。因此，DiPCA算法的非高斯过程监测性能还有待改进。

发明内容

本发明所要解决的主要技术问题是：提出一种新型动态独立成分分析建模算法，并利用该算法建立故障检测模型实施非高斯过程监测。

本发明解决上述技术问题所采用的技术方案为：一种基于新型动态独立成分分析的非高斯过程监测方法，包括以下所示步骤：

(1)采集生产过程正常运行状态下的样本，组成训练数据矩阵X∈R^n×m，并计算矩阵X中各列向量的均值μ₁，μ₂，…，μ_m以及标准差δ₁，δ₂，…，δ_m，对应组成均值向量μ＝[μ₁，μ₂，…，μ_m]^T与标准差向量δ＝[δ₁，δ₂，…，δ_m]，其中，n为训练样本数，m为过程测量变量数，R为实数集，R^n×m表示n×m维的实数矩阵，上标号T表示矩阵或向量的转置。

(2)根据如下所示公式对矩阵X实施标准化处理得到矩阵

上式(1)中，U∈R^n×m是由n个相同的均值向量μ组成的矩阵，即U＝[μ，μ，…，μ]^T，对角矩阵Ф中对角线上的元素由标准差向量δ组成。

(3)以

为训练数据，利用动态内部主成分分析算法求解模型载荷矩阵P与动态得分矩阵

其中，x_i∈R^m×1为标准化处理后的第i个样本数据，i＝1，2，…，n，Q为分解矩阵。

在给出动态内部主成分分析的详细实施过程之前，首先介绍下该算法的基本原理。设矩阵X_d＝[x_d，x_d+1，…，x_n-D+d-1]^T，其中下标号d＝1，2，…，D，D为自相关样本的个数。DiPCA算法的优化目标函数如下所示：

上式(2)中，X_D+1＝[x_D+1，x_D+2，…，x_n]^T，Z＝[X₁，X₂，…，X_D]，向量β＝[β₁，β₂，…，β_D]^T为权重系数向量、列向量p为载荷向量，

表示Kronecker乘积，具体的计算方式如下所示：

通过上述公式(2)可以发现，DiPCA算法的目标在于使提取出的特征成分的序列自相关性最大，但同时又保证了PCA各特征成分方差最大化的特点。通过拉格朗日乘子法即可解决公式(2)中定义的优化问题。

由于需要优化得到多个载荷向量p₁，p₂，…，P_M，可采取迭代求解的方式得到相应的结果，具体的实施过程如下所示：

①设置自相关样本数据个数D，并初始化a＝1与初始化

②初始化向量w_a∈R^m×1为m×m维单位矩阵中的第a列，根据公式t＝Yw_a计算得分向量t。

③根据公式β＝[t₁，t₂，…，t_D]^Tt_D+1计算向量β后，并根据公式β＝β/||β||单位化处理向量β，其中列向量t_j＝t(j：n-D+j-1)，下标号j＝1，2，…，D+1，t(j：n-D+j-1)表示取向量t中第j至第n-D+j-1个元素组成列向量。

④根据如下所示公式更新向量w_a：

上式(4)中，β_d为向量β中的第d个元素。

⑤根据公式w_a＝w_a/||w_a||对向量w_a实施单位化处理后，判断向量w_a是否收敛？若否，则计算t＝Yw_a后返回步骤③；若是，则得到第a个向量w_a后执行步骤⑥。

⑥根据公式p_a＝Y^Tt/(t^Tt)计算第a个载荷向量p_a，并根据公式Y＝Y-tp_a ^T更新数据矩阵Y。

⑦判断是否满足条件a＜m？若是，则置a＝a+1后返回步骤②；若否，则将m个载荷向量p₁，p₂，…，p_m以及向量w₁，w₂，…，w_m分别组成初始载荷矩阵P＝[p₁，p₂，…，p_m]与矩阵W＝[w₁，w₂，…，w_m]后执行步骤⑧。

⑧根据公式Q＝W(P^TW)计算分离矩阵Q，并根据公式

计算动态得分矩阵T。

(4)计算矩阵T中各列向量的方差λ₁，λ₂，…，λ_m，将方差小于10^-4的列向量从矩阵T中剔除从而得到矩阵T₀∈R^n×M，并对应地剔除载荷矩阵P中相应的列向量得到P₀，对应地剔除矩阵Q中相应的列向量得到Q₀，其中M为保留的动态得分向量的个数。

(5)根据公式Z₀＝T₀Λ^-1/2对动态得分矩阵T₀实施白化处理，其中Λ＝T₀ ^TT₀/(n-1)是一个对角矩阵。

(6)以Z₀作为独立成分的初始估计值，调用独立成分分析迭代算法求解得到转换矩阵C，从而使Z₀进一步转换成动态独立成分矩阵S＝Z₀C，具体的实施过程如下所示：

①初始化k＝1。

②设置向量c_k为M×M维单位矩阵中的第k列。

③按照公式c_k←E{Z₀g(c_k ^TZ₀)}-E{h(c_k ^TZ₀)}c_k更新向量c_k，其中函数g(u)＝tanh(u)、函数h(u)＝[sech(u)]²、u＝c_k ^TZ₀表示函数自变量、E{}表示计算均值。

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

c_k←c_k/||c_k|| (6)

⑤重复步骤③～④直至向量c_k收敛，并保存向量c_k。

⑥判断是否满足条件k＜M？若是，则置k＝k+1后返回步骤②；若否，则将得到的向量c₁，c₂，…，c_M组成初始转换矩阵C₀＝[c₁，c₂，…，c_M]。

⑦根据公式S₀＝Z₀C₀计算得到动态潜独立成分矩阵S₀，S₀中各列向量即为各独立成分向量，需按照非高斯性大小做进一步排序处理。

⑧按照均值为0方差为1的高斯分布随机生成一个n×1维的列向量ε，并按照如下所示公式计算S₀中各独立成分向量的非高斯性大小J₁，J₂，…，J_M：

J_k＝|E{log cosh(s_k)}-E{log cosh(ε)}| (7)

上式(7)中，s_k为矩阵S₀中第k列。

⑨根据J₁，J₂，…，J_M的数值大小进行降序排序，并对应地变更矩阵S₀中各列向量的先后排序得到S，以及矩阵C₀中各列向量的先后排序得到转换矩阵C。

(7)设置保留的动态独立成分个数η，将转换矩阵C中第η+1行至第M行的向量删除，从而得到新转换矩阵C_η，并根据公式S_η＝Z₀C_η计算新动态独立成分矩阵S_η。

(8)根据公式I²＝diag(S_ηS_η ^T)与SPE＝diag(E_ηE_η ^T)计算训练数据对应的监测统计指标I²与SPE，其中diag()表示将矩阵对角线的元素单独作为列向量的操作运算，模型误差矩阵

(9)计算I²指标的均值

与标准差

并计算SPE指标的均值μ_SPE与标准差δ_SPE，那么监测统计量I²与SPE的控制上限分别为

与SPE_lim＝μ_SPE+3δ_SPE。

(10)保留如下所示模型参数以备在线监测时调用：步骤(1)中的均值向量μ、步骤(2)中的对角矩阵Ф、步骤(4)中的载荷矩阵P₀与分离矩阵Q₀、步骤(5)中的对角矩阵Λ、步骤(7)中的转换矩阵C_η、以及步骤(9)中的控制上限

与SPE_lim。

(11)收集新采样时刻的数据样本x∈R^m×1，并根据公式

对x实施标准化处理得到

(12)根据如下所示公式计算在线样本数据的动态独立成分s与模型误差e：

(13)根据公式I²＝ss^T与SPE＝ee^T计算得到在线样本数据的监测统计指标I²与SPE。

(14)判断是否满足条件

且SPE≤SPE_lim？若是，则当前样本采集自正常工况，返回步骤(11)继续监测下一时刻的样本数据；若否，则当前监测样本为故障样本。

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

首先，本发明方法将动态内部主成分分析算法作为独立成分分析算法的白化过程，发挥了DiPCA方法能挖掘动态自相关特征的能力。其次，本发明方法对动态特征成分做进一步转换，使得到的动态独立成分信息挖掘出数据中潜藏的非高斯特征成分。因此，本发明方法是一种全新的非高斯过程监测方法。此外，具体实施案例中将会验证本发明方法的优越性，从而说明本发明方法是一种更为优选的非高斯过程监测方法。

附图说明

图1为本发明方法离线建模阶段的实施流程图。

图2为本发明方法在线监测的实施流程图。

图3为TE过程故障监测详情对比图。

具体实施方式

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

本发明公开一种基于新型动态独立成分分析的非高斯过程监测方法，下面结合一个具体的工业过程的例子来说明本发明方法的具体实施过程，以及相对于现有方法的优越性。

应用对象是来自于美国田纳西-伊斯曼(TE)化工过程实验，原型是伊斯曼化工生产车间的一个实际工艺流程。目前，TE过程因其流程的复杂性，已作为一个标准实验平台被广泛用于故障检测研究。整个TE过程包括22个测量变量、12个操作变量、和19个成分测量变量。该TE过程对象可以模拟仿真多种不同的故障类型，如物料进口温度阶跃变化、冷却水故障变化等等。为了对该过程进行监测，选取如表1所示的33个过程变量。由于采样间隔时间较短，TE过程采样数据不可避免的存在序列自相关性，接下来结合该TE过程对本发明具体实施步骤进行详细的阐述。

表1：TE过程监测变量。

序号	变量描述	序号	变量描述	序号	变量描述
						1	物料A流量	12	分离器液位	23	D进料阀门位置
2	物料D流量	13	分离器压力	24	E进料阀门位置
						3	物料E流量	14	分离器塔底流量	25	A进料阀门位置
4	总进料流量	15	汽提塔等级	26	A和C进料阀门位置
						5	循环流量	16	汽提塔压力	27	压缩机循环阀门位置
6	反应器进料	17	汽提塔底部流量	28	排空阀门位置
						7	反应器压力	18	汽提塔温度	29	分离器液相阀门位置
8	反应器等级	19	汽提塔上部蒸汽	30	汽提塔液相阀门位置
						9	反应器温度	20	压缩机功率	31	汽提塔蒸汽阀门位置
10	排空速率	21	反应器冷却水出口温度	32	反应器冷凝水流量
						11	分离器温度	22	分离器冷却水出口温度	33	冷凝器冷却水流量

首先，离线建模的实施流程如图1所示，利用TE过程正常工况下采样的960个样本数据建立非高斯过程监测模型，包括以下步骤：

(1)采集生产过程正常运行状态下的样本，组成训练数据矩阵X∈R^960×33，并计算矩阵X的均值向量μ＝[μ₁，μ₂，…，μ₃₃]^T与标准差向量δ＝[δ₁，δ₂，…，δ₃₃]。

(2)根据公式

对矩阵X实施标准化处理得到

其中U∈R^n×m是由n个相同的均值向量μ组成的矩阵，对角矩阵Ф中对角线上的元素由标准差向量δ组成。

(3)以

为训练数据，利用动态内部主成分分析算法求解模型载荷矩阵P与动态得分矩阵

(4)计算矩阵T中各列向量的方差λ₁，λ₂，…，λ₃₃，将方差小于10^-4的列向量从矩阵T中剔除从而得到矩阵T₀∈R^960×31，并对应地剔除载荷矩阵P中相应的列向量得到P₀，对应地剔除矩阵Q中相应的列向量得到Q₀。

(5)根据公式Z₀＝T₀Λ^-1/2对动态得分矩阵T₀实施白化处理，其中Λ＝T₀ ^TT₀/(n-1)是一个对角矩阵。

(6)以Z₀作为独立成分的初始估计值，调用独立成分分析迭代算法求解得到转换矩阵C，从而使Z₀进一步转换成动态独立成分矩阵S＝Z₀C。

(7)设置保留的动态独立成分个数η＝9，将转换矩阵C中第η+1行至第M行的向量删除，从而得到新转换矩阵C_η，并根据公式S_η＝Z₀C_η计算新动态独立成分矩阵S_η。

(8)根据公式I²＝diag(S_ηS_η ^T)与SPE＝diag(E_ηE_η ^T)计算训练数据对应的监测统计指标I²与SPE。

(9)计算I²指标的均值

与标准差

并计算SPE指标的均值μ_SPE与标准差δ_SPE，那么监测统计量I²与SPE的控制上限分别为

与SPE_lim＝μ_SPE+3δ_SPE。

(10)保留如下所示模型参数以备在线监测时调用：步骤(1)中的均值向量μ、步骤(2)中的对角矩阵Ф、步骤(4)中的载荷矩阵P₀与分离矩阵Q₀、步骤(5)中的对角矩阵Λ、步骤(7)中的转换矩阵C_η、以及步骤(9)中的控制上限

与SPE_lim。

实施在线监测的流程如图2所示，利用TE过程故障条件下的960个采样数据实施在线故障检测，其中前160个样本数据为正常，后800个数据才采集自故障工况。

(11)收集新采样时刻的数据样本x∈R^33×1，并根据公式

对x实施标准化处理得到

(12)计算在线样本数据的动态独立成分s与模型误差e。

(13)根据公式I²＝ss^T与SPE＝ee^T计算得到在线样本数据的监测统计指标I²与SPE。

(14)判断是否满足条件

且SPE≤SPE_lim？若是，则当前样本采集自正常工况，返回步骤(11)继续监测下一时刻的样本数据；若否，则当前监测样本为故障样本。

如图3所示，本发明方法与基于ICA的传统非高斯过程监测方法在监测故障数据时的监测详情。从图3中可以很明显地发现，本发明方法在故障发生后漏报情况明显优越于基于ICA的非高斯过程监测方法。因此，可以说本发明方法具有更可靠的过程监测性能。

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

Claims (1)

1.一种基于新型动态独立成分分析的非高斯过程监测方法，其特征在于，包括以下步骤：

步骤(1)：采集生产过程正常运行状态下的样本，组成训练数据矩阵X∈R^n×m，并计算矩阵X中各列向量的均值μ₁，μ₂，…，μ_m以及标准差δ₁，δ₂，…，δ_m，对应组成均值向量μ＝[μ₁，μ₂，…，μ_m]^T与标准差向量δ＝[δ₁，δ₂，…，δ_m]，其中，n为训练样本数，m为过程测量变量数，R为实数集，R^n×m表示n×m维的实数矩阵，上标号T表示矩阵或向量的转置；

步骤(2)：根据如下所示公式对矩阵X实施标准化处理得到矩阵

上式(1)中，U∈R^n×m是由n个相同的均值向量μ组成的矩阵，即U＝[μ，μ，…，μ]^T，对角矩阵Φ中对角线上的元素由标准差向量δ组成；

步骤(3)：以

为训练数据矩阵，利用动态内部主成分分析算法求解模型载荷矩阵P与动态得分矩阵

具体的实施过程如下所示；其中，x_i∈R^m×1为标准化处理后的第i个样本数据，i＝1，2，…，n，Q为分解矩阵；

①设置自相关样本数据个数D，并初始化a＝1与初始化

②初始化向量w_a∈R^m×1为m×m维单位矩阵中的第a列，根据公式t＝Yw_a计算得分向量t；

③根据公式β＝[t₁，t₂，…，t_D]^Tt_D+1计算向量β后，并根据公式β＝β/||β||单位化处理向量β，其中列向量t_j＝t(j：n-D+j-1)，t(j：n-D+j-1)表示取向量t中第j至第n-D+j-1个元素组成列向量，j＝1，2，…，D+1；

④根据如下所示公式更新向量w_a：

上式中，β_d为向量β中的第d个元素；

⑤根据公式w_a＝w_a/||w_a||对向量w_a实施单位化处理后，判断向量w_a是否收敛；若否，则计算t＝Yw_a后返回步骤③；若是，则得到第a个向量w_a后执行步骤⑥；

⑥根据公式p_a＝Y^Tt/(t^Tt)计算第a个载荷向量p_a，并根据公式Y＝Y-tp_a ^T更新数据矩阵Y；

⑦判断是否满足条件a＜m；若是，则置a＝a+1后返回步骤②；若否，则将m个载荷向量p₁，p₂，…，p_m以及向量w₁，w₂，…，w_m分别组成初始载荷矩阵P＝[p₁，p₂，…，p_m]与矩阵W＝[w₁，w₂，…，w_m]后执行步骤⑧；

⑧根据公式Q＝W(P^TW)计算分离矩阵Q，并根据公式

计算动态得分矩阵T；

步骤(4)：计算动态得分矩阵T中各列向量的方差λ₁，λ₂，…，λ_m，将方差小于10^-4的列向量从T中剔除从而得到矩阵T₀∈R^n×M，并对应地剔除载荷矩阵P中相应的列向量得到P₀，对应地剔除矩阵Q中相应的列向量得到Q₀，其中M为保留的动态得分向量的个数；

步骤(5)：根据公式Z₀＝T₀Λ^-1/2对矩阵T₀实施白化处理，其中，Λ＝T₀ ^TT₀/(n-1)为对角矩阵；

步骤(6)：以Z₀作为独立成分的初始估计值，调用独立成分分析迭代算法求解得到转换矩阵C，从而使Z₀进一步转换成动态独立成分矩阵S＝Z₀C，具体的实施过程如下所示；

①初始化k＝1；

②设置向量c_k为M×M维单位矩阵中的第k列；

③按照公式c_k←E{Z₀g(c_k ^TZ₀)}-E{h(c_k ^TZ₀)}c_k更新向量c_k，其中函数g(u)＝tanh(u)、函数h(u)＝[sech(u)]²、u＝c_k ^TZ₀表示函数自变量、E{}表示计算均值；

④对更新后的向量c_k依次按照如下所示公式进行正交标准化处理：

c_k←c_k/||c_k|| (4)

⑤重复步骤③～④直至向量c_k收敛，并保存向量c_k；

⑥判断是否满足条件k＜M；若是，则置k＝k+1后返回步骤②；若否，则将得到的向量c₁，c₂，…，c_M组成初始转换矩阵C₀＝[c₁，c₂，…，c_M]；

⑦根据公式S₀＝Z₀C₀计算得到动态潜独立成分矩阵S₀，S₀中各列向量即为各独立成分向量；

⑧按照均值为0方差为1的高斯分布随机生成一个n×1维的列向量ε，并按照如下所示公式计算S₀中各独立成分向量的非高斯性大小J₁，J₂，…，J_M：

J_k＝|E{log cosh(s_k)}-E{log cosh(ε)}| (5)

上式中，s_k为矩阵S₀中的第k列向量；

⑨根据J₁，J₂，…，J_M的数值大小进行降序排序，并对应地变更矩阵S₀中各列向量的先后排序得到S，以及矩阵C₀中各列向量的先后排序得到转换矩阵C；

步骤(7)：设置保留的动态独立成分个数η，将转换矩阵C中第η+1行至第M行的向量删除，从而得到新转换矩阵C_η，并根据公式S_η＝Z₀C_η计算新动态独立成分矩阵S_η；

步骤(8)：根据公式I²＝diag(S_ηS_η ^T)与SPE＝diag(E_ηE_η ^T)计算训练数据对应的监测统计指标I²与SPE，其中diag()表示将矩阵对角线的元素单独作为列向量的操作，模型误差矩阵

步骤(9)：计算I²指标的均值

与标准差

并计算SPE指标的均值μ_SPE与标准差δ_SPE，那么监测统计量I²与SPE的控制上限分别为

与SPE_lim＝μ_SPE+3δ_SPE；

步骤(10)：保留如下所示模型参数以备在线监测时调用：步骤(1)中的均值向量μ、步骤(2)中的对角矩阵Φ、步骤(4)中的载荷矩阵P₀与分离矩阵Q₀、步骤(5)中的对角矩阵Λ、步骤(7)中的转换矩阵C_η、以及步骤(9)中的控制上限

与SPE_lim；

步骤(11)：收集新采样时刻的数据样本x∈R^m×1，并根据公式

对x实施标准化处理得到

步骤(12)：根据如下所示公式计算在线样本数据的动态独立成分s与模型误差e：

步骤(13)：根据公式I²＝ss^T与SPE＝ee^T计算得到在线样本数据的监测统计指标I²与SPE；

步骤(14)：判断是否满足条件

且SPE≤SPE_lim；若是，则当前样本采集自正常工况，返回步骤(11)继续监测下一时刻的样本数据；若否，则当前监测样本为故障样本。

Images