CN111913460B - 一种基于序列相关局部保持投影算法的故障监测方法 - Google Patents

Abstract

本发明公开一种基于序列相关局部保持投影算法的故障监测方法，旨在推理出一种全新的数据特征挖掘算法，并基于此算法实施故障监测。具体来讲，本发明方法首先将时序相关性最大化与局部近邻结构保持最小化合并成一个目标函数；其次，在求解的过程进一步保证各投影变换向量之间的相互正交特性；最后，利用提取的潜在特征以及模型误差实施故障监测。与传统方法相比：本发明方法中涉及的序列相关局部保持投影算法是一种全新的特征提取算法，它在投影变换的过程中同时考虑了自相关性特征与局部近邻特征，并且保证了投影变换向量的正交特性，能够更全面地挖掘训练数据中潜藏的有用信息。因此，本发明方法是一种更为优选的故障监测方法。

Description

一种基于序列相关局部保持投影算法的故障监测方法

技术领域

本发明涉及一种数据驱动的故障监测方法，尤其涉及一种基于序列相关局部保持投影算法的故障监测方法。

背景技术

对工业过程对象的运行状态进行实时监测是保证生产安全运行与维持产品质量稳定性的直接保障，工业过程对象综合自动化体系的设计离不开故障监测系统。当前，工业发展已进入以“大数据”为代表的信息化建设阶段，利用采样数据来监测过程运行状态是否发生故障测，早已成为工业自动化领域的研究热门之一。通常来讲，数据驱动的故障监测方法的核心思想在于：如何对过程正常数据进行有效地挖掘以提取能反应过程运行状态的潜在有用信息。由于现代工业过程规模的复杂化趋势，挖掘采样数据的单一特征往往无法体现出的现代工业过程对象的复杂特性。可以说，如何全面挖掘出过程数据中潜藏的有用信息，并建立更适于监测现代工业过程对象的故障监测模型，一直以来都是该研究领域所面临的主要问题。

在各式各样的特征提取算法中，提取训练数据的时序相关动态特征与近邻局部结构特征是较为常见的。其中，时序相关动态特征的提取一般是通过自回归模型或增广矩阵的思路来实现。最为代表性的动态建模方法是：动态主元分析(Dynamic PrincipalComponent Analysis，缩写：DPCA)。此外，局部保持投影(Locality PreservingProjections，缩写：LPP)算法是最经典的近邻局部结构特征方法。LPP在投影变换的过程中尽量保留数据点在空间距离上的近邻分布特征。虽然，通过增广矩阵引入延时测量值可将数据的时序自相关性考虑进来，可将LPP直接拓展成动态LPP方法。但是，训练数据的时序相关特征却没有得到充分的挖掘。

在时序相关特征挖掘方面，可通过最大化潜在特征成分的时序相关性为目标，设计新的特征提取算法。近年来，有学者曾提出动态潜变量建模方法与动态内部PCA方法，考虑了潜在特征成分的时序相关性。然而，由于局部近邻特征与自相关特征都是数据中潜藏的有用特征，在实施故障监测时，这两个特征的提取都应予以充分考虑。

发明内容

本发明所要解决的主要技术问题是：设计出一种能同时提取时序相关性与近邻局部特征的投影变换算法，并基于此算法实施故障监测。具体来讲，本发明方法首先将时序相关性最大化与局部近邻结构保持最小化合并成一个目标函数；其次，在求解的过程进一步保证各投影变换向量之间的相互正交特性；最后，利用提取的潜在特征以及模型误差实施故障监测。

本发明解决上述技术问题所采用的技术方案为：一种基于序列相关局部保持投影算法的故障监测方法，包括以下步骤：

步骤(1)：在生产过程正常运行状态下，按照采样时间先后依次采集n个样本数据组成矩阵X∈R^n×m，并计算矩阵X中各行向量的均值向量μ与标准差向量δ，其中m为测量变量的总个数、R为实数集，R^n×m表示n×m维的实数矩阵。

步骤(2)：利用均值向量μ与标准差向量δ对矩阵X中各行向量实施标准化处理，从而得到矩阵

其中

为标准化处理后的数据向量，i＝1，2，…，n，上标号T表示矩阵或向量的转置。

步骤(3)：设置自相关阶数为d，按照如下所示公式①分别构造Future矩阵X_F与Past矩阵X_P：

从公式①可以看出，Future矩阵X_F与Past矩阵X_P中相同行的行向量体现了采样数据的时序性。为了使潜在特征的时序相关性最大化，特设计如下所示的目标函数：

上式中，p∈R^m×1表示投影变换向量，β＝[β₁，β₂，…，β_d]^T∈R^d×1为序列相关系数向量，d为自相关阶数(一般可取d＝2)，符号

表示Kronecker乘法，

的计算结果如下所示：

从上式②中可以看出，经投影变换向量p∈R^m×1转换后的潜在特征需满足时序自相关性最大的要求。

此外，LPP算法旨在投影变换后依旧能保持各个样本数据的局部近邻结构特征，其实现的目标函数如下所示：

上式中，W∈R^n×n表示毗邻矩阵，对角矩阵D∈R^n×n对角线上的元素为毗邻矩阵W中各列向量之和，即

D_ii表示对角矩阵D中第i行第i列元素。其中毗邻矩阵W中第i行第j列元素W_ij的具体计算方式如下所示：

若想同时考虑序列相关特征与局部近邻结构特征，需要在求解多个投影变换向量p₁，p₂，…，p_A的同时，满足公式②与公式④中的目标函数。由于J₁是最大化问题而J₂是最小化问题，因此可设计如下所示的目标函数：

其中，矩阵L＝D-W。为不失一般性，可设

此外，为保证A个投影变换向量之间相互正交，还需增加正交约束条件：p_a ^Tp₁＝p_a ^Tp₂＝…＝p_a ^Tp_a-1＝0，其中a＝1，2，…，A表示求解的第a个投影变换向量。上式⑥最终转变成如下所示的优化问题：

再通过拉格朗日乘子法即可求解上式⑦：构造如下所示的拉格朗日函数L：

计算函数L相对于p_a与β_a的偏微分：

上式中，I_m与I_d分别表示m×m维与d×d维的单位矩阵。根据极值求解思路，上式⑨与⑩中的偏微分都等于零。因此，可推理出：

上两式中，

在公式

中等号两边同时左乘p_a ^T，即可得到：

所以，拉格朗日乘子λ即等价于公式⑥中的目标函数值。

在公式

中等号两边同时依次左乘

则可得到a-1个等式关系：

上式

则可以等价的写成如下所示的矩阵形式：

上式

中，γ^(a-1)＝[γ₁，γ₂，…，γ_a-1]^T，A^(a-1)＝[p₁，p₂，…，p_a-1]，

将公式

等号两边同时左乘

后再减去公式

即可得到：

因此，公式

定义了一个常规的特征值问题，投影变换向量p_a即是最大特征值所对应的特征向量。此外，公式

同样定义了一个常规的特征值问题，序列相关系数向量β_a即为最大特征值所对应的特征向量。

由于矩阵G的计算涉及到序列相关系数向量β_a，而β_a的求解需要在p_a已知的前提下根据公式

计算得到。所以，p_a与β_a的求解过程相互耦合，可通过如下所示的迭代收敛方式予以同时求解。

步骤(A)：初始化a＝1与初始化β_a为任意一个d×1维的实数向量。

步骤(B)：根据公式

计算矩阵G后，判断是否满足条件：a＜2？若是，则求解特征值问题：

最大特征值λ所对应的特征向量p_a；若否，则求解特征值问题：

最大特征值λ所对应的特征向量p_a。

步骤(C)：对p_a实施归一化处理

并求解特征值问题Hβ_a＝ηβ_a最大特征值对应的特征向量β_a，其中

步骤(D)：对β_a实施归一化处理

后，判断β_a是否收敛？若否，则返回步骤(B)；若是，则得到第a个投影变换向量p_a，并执行步骤(E)。

步骤(E)：判断是否满足条件：a＜A？若是，则置a＝a+1后，初始化β_a为任意一个d×1维的实数向量，再返回步骤(B)；若否，则得到A个投影变换向量p₁，p₂，…，p_A。

步骤(4)：设置投影变换向量的个数为A，根据上述步骤(A)至步骤(E)求解得到A个投影变换向量p₁，p₂，…，p_A后，即可建立模型：

其中A＜m，得分矩阵

载荷矩阵P＝[p₁，p₂，…，p_A]，误差矩阵

步骤(5)：根据公式Λ＝S^TS/(n-1)计算协方差矩阵Λ后，利用主元分析算法对误差矩阵E实施分解：

具体的实施过程如下所示：

步骤(5.1)：根据公式C＝E^TE/(n-1)计算得到协方差矩阵C后，求解特征值问题：Cv＝ζv所有特征值ζ₁≥ζ₂≥…≥ζ_m所对应的特征向量v₁，v₂，…，v_m，这里要求特征向量都为单位长度。

步骤(5.2)：设置主元个数θ为满足条件：(ζ₁+ζ₂+…+ζ_θ)/(ζ₁+ζ₂+…+ζ_m)≥0.85的最小值后，将v₁，v₂，…，v_θ组成主元载荷矩阵V＝[v₁，v₂，…，v_θ]，并将θ个特征值ζ₁，ζ₂，…，ζ_θ组成对角矩阵Θ。

步骤(5.3)：分别根据公式

计算主元得分矩阵U与主元误差矩阵

步骤(6)：根据公式ψ＝diag{SΛ^-1S^T}、ξ＝diag{UΘ^-1UT}、和

分别计算监测指标向量ψ、ξ、和Q，其中diag{}表示将矩阵对角线上的元素变成向量的操作。

步骤(7)：分别将监测指标向量ψ、ξ、和Q中第n/100个最大元素记做监测指标上限ψ_lim、ξ_lim、和Q_lim。

以上步骤(1)至步骤(7)为本发明方法的离线建模阶段，离线建模完成后，即可实施对工业过程对象的在线监测。

步骤(8)：采集最新采样时刻的样本数据x_new∈R^m×1，并利用均值向量μ与标准差向量δ对x_new实施标准化处理得到向量

步骤(9)：根据公式

u_new＝e_newV、和

分贝计算得分向量s_new、误差向量e_new、主元得分向量u_new、和主元误差向量

步骤(10)：根据公式ψ_new＝s_newΛ^-1s_new ^T、ξ_new＝u_newΘ^-1u_new ^T、

分别计算出监测指标ψ_new、ξ_new、和Q_new的具体数值。

步骤(11)：判断是否满足条件：ψ_new≤ψ_lim且ξ_new≤ξ_lim且Q_new≤Q_lim？若是，则当前采样时刻过程正常运行，返回步骤(8)继续实施对下一采样时刻数据的监测；若否，则当前采样时刻过程进入异常工况状态。

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

首先，本发明方法中涉及的序列相关局部保持投影算法是一种全新的特征提取算法，它在投影变换的过程中同时考虑了自相关性特征与局部近邻特征，并且保证了投影变换向量的正交特性，能够更全面地挖掘训练数据中潜藏的有用信息。其次，在具体实施案例中，通过TE过程上的监测详情对比，验证了本发明方法相比于传统故障监测方法在故障监测上的优越性。因此，本发明方法是一种更为优选的故障监测方法。

附图说明

图1为本发明方法的实施流程图。

图2为序列相关局部保持投影算法的实施流程图。

图3本发明方法与传统方法的故障监测效果对比图。

具体实施方式

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

如图1所示，本发明公开一种基于序列相关局部保持投影算法的故障监测方法。现结合一个具体的实施案例来阐述本发明方法的具体实施方式。

所测试的过程对象为TE过程，该过程原型是伊斯曼化工生产车间的一个实际工艺流程。目前，TE过程因其流程的复杂性，已作为一个标准实验平台被广泛用于故障检测研究。整个TE过程包括22个测量变量、12个操作变量、和19个成分测量变量。所采集的数据分为22组，其中包括1组正常工况下的数据集与21组故障数据。而在这些故障数据中，有16个是已知故障类型，如冷却水入口温度或进料成分的变化、阀门粘滞、反应动力学漂移等，还有5个故障类型是未知的。为了对该过程进行监测，选取如表1所示的33个过程变量，接下来结合该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过程正常运行状态下，按照采样时间先后依次采集n＝960个样本数据x₁，x₂，…，x_n组成矩阵X＝[x₁，x₂，…，x_n]^T∈R^960×33，并计算矩阵X中各行向量的均值向量μ与标准差向量δ。

步骤(2)：利用均值向量μ与标准差向量δ对矩阵X中各行向量实施标准化处理，从而得到矩阵

其中

为标准化处理后的数据向量。

步骤(3)：设置自相关阶数为d，按照上述公式①分别构造Future矩阵X_F与Past矩阵X_P。

步骤(4)：设置投影变换向量个数为A＝10与近邻阶数为k＝5，按照如图2所示的实施流程求解得到A个投影变换向量p₁，p₂，…，p_A后，即可建立模型：

步骤(5)：根据公式Λ＝S^TS/(n-1)计算协方差矩阵Λ后，利用主元分析算法对误差矩阵E实施分解：

步骤(6)：根据公式ψ＝diag{SΛ^-1S^T}、ξ＝diag{UΘ^-1U^T}、和

分别计算监测指标向量ψ、ξ、和Q，其中diag{}表示将矩阵对角线上的元素组成向量的操作。

步骤(7)：分别将监测指标向量ψ、ξ、和Q中第n/100≈10个最大元素记做监测指标上限ψ_lim、ξ_lim、和Q_lim。

离线建模阶段完成后，即可实施在线过程监测。采集TE过程在第19种故障工况条件下的测试数据，其中前160个样本数据采集自正常运行状态，故障样本数据从第161个采样时刻引入。

步骤(8)：采集最新采样时刻的样本数据x_new∈R^33×1，并利用均值向量μ与标准差向量δ对x_new实施标准化处理得到向量

步骤(9)：根据公式

u_new＝e_newV、和

分贝计算得分向量s_new、误差向量e_new、主元得分向量u_new、和主元误差向量

步骤(10)：根据公式ψ_new＝s_newΛ^-1s_new ^T、ξ_new＝u_newΘ^-1u_new ^T、

分别计算出监测指标ψ_new、ξ_new、和Q_new的具体数值。

步骤(11)：判断是否满足条件：ψ_new≤ψ_lim且ξ_new≤ξ_lim且Q_new≤Q_lim？若是，则当前采样时刻过程正常运行，返回步骤(8)继续实施对下一采样时刻数据的监测；若否，则当前采样时刻过程进入异常工况状态。

将本发明方法与其他传统过程监测方法(包括LPP与动态PCA)对TE过程该故障的监测详情对比于图3中。从图3中可以明显地发现，本发明方法的故障漏报情况更少，监测效果得到了有效的改善。

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

Claims (1)

1.一种基于序列相关局部保持投影算法的故障监测方法，其特征在于，包括以下步骤：

首先，离线建模阶段包括如下所示步骤(1)至步骤(7)；

步骤(1)：在生产过程正常运行状态下，按照采样时间先后依次采集n个样本数据组成矩阵X∈R^n×m，并计算矩阵X中各行向量的均值向量μ与标准差向量δ，其中m为测量变量的总个数、R为实数集，R^n×m表示n×m维的实数矩阵；

步骤(2)：利用均值向量μ与标准差向量δ对矩阵X中各行向量实施标准化处理，从而得到矩阵

其中

为标准化处理后的数据向量，i＝1，2，…，n，上标号T表示矩阵或向量的转置；

步骤(3)：设置自相关阶数为d，按照如下所示公式①分别构造Future矩阵X_F与Past矩阵X_P：

步骤(4)：设置投影变换向量个数为A与近邻阶数为k，根据如下所示步骤(4.1)至步骤(4.5)求解得到A个投影变换向量p₁，p₂，…，p_A后，即可建立模型：

其中得分矩阵

载荷矩阵P＝[p₁，p₂，…，p_A]、误差矩阵

步骤(4.1)：根据如下所示公式②计算毗邻矩阵W中第i行第j列元素W_ij：

上式中，i＝1，2，…，n、j＝1，2，…，n；

步骤(4.2)：计算矩阵L＝D-W，其中对角矩阵D∈R^n×n对角线上的元素为毗邻矩阵W中各列向量之和，并初始化a＝1与初始化β_a为任意一个d×1维的实数向量；

步骤(4.3)：根据公式

计算矩阵G后，判断是否满足条件：a＜2；若是，则求解特征值问题：

最大特征值λ所对应的特征向量p_a；若否，则求解特征值问题：

最大特征值所对应的特征向量p_a；其中，I_m表示m×m维的单位矩阵、符号

表示Kronecker乘法、

A^(a-1)＝[p₁，p₂，…，p_a-1]、

步骤(4.4)：对p_a实施归一化处理

并求解特征值问题Hβ_a＝ηβ_a最大特征值η对应的特征向量β_a，其中

I_d表示d×d维的单位矩阵；

步骤(4.5)：对β_a实施归一化处理

后，判断β_a是否收敛；若否，则返回步骤(4.3)；若是，则得到第a个投影变换向量p_a，并执行步骤(4.6)；

步骤(4.6)：判断是否满足条件：a＜A；若是，则置a＝a+1后，初始化β_a为任意一个d×1维的实数向量，再返回步骤(4.3)；若否，则得到A个投影变换向量p₁，p₂，…，p_A；

步骤(5)：根据公式Λ＝S^TS/(n-1)计算协方差矩阵Λ后，利用主元分析算法对误差矩阵E实施分解：

其中U表示主元得分矩阵、V表示主元载荷矩阵、

表示主元误差矩阵；

步骤(6)：根据公式ψ＝diag{SΛ^-1S^T}、ξ＝diag{UΘ^-1U^T}、和

分别计算监测指标向量ψ、ξ、和Q，其中diag{}表示将矩阵对角线上的元素变成向量的操作；

步骤(7)：分别将监测指标向量ψ、ξ、和Q中第n/100个最大元素记做监测指标上限ψ_lim、ξ_lim、和Q_lim；

其次，在线故障监测阶段包括如下所示步骤(8)至步骤(11)；

步骤(8)：采集最新采样时刻的样本数据x_new∈R^m×1，并利用均值向量μ与标准差向量δ对x_new实施标准化处理得到向量

步骤(9)：根据公式

u_new＝e_newV、和

分别计算得分向量s_new、误差向量e_new、主元得分向量u_new、和主元误差向量

步骤(10)：根据公式ψ_new＝s_newΛ^-1s_new ^T、ξ_new＝u_newΘ^-1u_new ^T、

分别计算出监测指标ψ_new、ξ_new、和Q_new的具体数值；

步骤(11)：判断是否满足条件：ψ_new≤ψ_lim且ξ_new≤ξ_lim且Q_new≤Q_lim；若是，则当前采样时刻过程正常运行，返回步骤(8)继续实施对下一采样时刻数据的监测；若否，则当前采样时刻过程进入异常工况状态。

Images