CN109522972A - 一种基于潜变量自回归模型的动态过程监测方法 - Google Patents

Abstract

本发明公开一种基于潜变量自回归模型的动态过程监测方法，旨在建立潜变量自回归模型，并以此为基础实施动态过程监测。具体来讲，本发明方法通过定义潜变量的自回归模型的最小二乘目标函数，推理出相应的特征挖掘算法后，再建立故障监测模型从而实施在线故障监测。由于本发明方法以建立潜变量自回归模型的目标，挖掘出了动态自相关的潜变量，并给出相应满足最小二乘条件的自回归模型。通过该潜变量自回归模型不仅可以挖掘出原训练数据中的自相关特征，而且还可以消除潜变量自相关性的影响。因此，本发明方法与传统动态过程监测方法是存在显著不同的，而且模型的可解释性更强。可以说，本发明方法是一种更为优选的动态过程监测方法。

Description

一种基于潜变量自回归模型的动态过程监测方法

技术领域

本发明涉及一种数据驱动的过程监测方法，尤其涉及一种基于潜变量自回归模型的动态过程监测方法。

背景技术

近年来，各行各业都开始兴起“大数据”的研究与应用热潮。与此同时，对数据的利用程度体现了工业过程对象的现代化程度。在工业界，尤其是流程工业生产车间，先进仪表技术与计算技术都得到了非常广泛的应用，生产过程对象可以离线存储与在线测量海量的数据。这些数据为工业过程研究与应用“大数据”方法提供了坚实的数据基础。以生产安全为例，这些采样数据中潜藏着能体现生产过程运行状态的信息，利用采样数据即可实施过程运行状态的监测。在最近的十几年内，工业界与学术界都投入了大量的人力与物力研究以故障监测核心任务的过程监测方法。在数据驱动的过程监测研究领域，统计过程监测是被研究得最多的方法，其中当以主元分析(Principal Component Analysis，PCA)算法为主流的实施技术手段。数据驱动的过程监测的核心本质在于：通过对正常工况下的采样数据实施特征挖掘，并利用该特征建立单分类的模型以实施监测。

如上所述，正常工况下采样数据的特征挖掘是建立数据驱动的过程监测模型时最关键的一步。考虑到工业过程对象的采样时间间隔较短，采样数据不可避免地存在时间序列上的自相关性。因此，数据的自相关性这种动态特征是必须考虑的一个问题。针对动态过程监测问题的研究，最常见的思路就是使用增广矩阵，将数据的自相关性与交叉相关性混淆在一起后，利用PCA算法实施特征提取。这就是科研文献中经典的动态PCA方法，在经典动态PCA方法的基础上，也有学者发明了多种新颖的建模思路，比如动态潜变量模型与动态内部PCA模型。通过定义不同的目标函数，可以有效挖掘出训练数据的动态自相关特征。此外，自回归模型也曾被用于动态过程监测，通过建立前后采样时刻的回归模型，以达到描述采样数据自相关特征的目的。

然而，潜在特征同样是存在自相关特征的。换句话说，潜在特征或潜变量同样需要自回归模型来描述动态自相关特征。但是，纵观已有的科研文献与专利资料中，没有涉及这方面的研究与应用。因此，推理出一种潜变量自回归模型很有必要。可以通过少数几个潜变量的自回国模型来挖掘训练数据中潜藏的有用特征，基于该模型即可实施新颖的动态过程监测。

发明内容

本发明所要解决的主要技术问题是：如何建立潜变量自回归模型，并以此为基础实施动态过程监测。为此，本发明公开一种基于潜变量自回归模型的动态过程监测方法。通过定义潜变量的自回归模型的最小二乘目标函数，推理出相应的特征挖掘算法后，再建立故障监测模型从而实施在线故障监测。

本发明解决上述技术问题所采用的技术方案为：一种基于潜变量自回归模型的动态过程监测方法，包括以下所示步骤：

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

(2)根据公式对矩阵X实施标准化处理得到矩阵其中，Ξ∈R^n×m是由n个相同的均值向量μ组成的矩阵，即Ξ＝[μ，μ，…，μ]^T。

(3)记矩阵设置自相关阶数为D后，根据X_d＝[x_d，x_d+1，…，x_n-D+d-1]^T依次构造矩阵X₁，X₂，…，X_D+1，其中下标号d＝1，2，…，D+1，x_i为标准化后的第i个数据样本，i＝1，2，…，n。

特征挖掘的关键在于找到一个投影向量w∈R^m×1，使投影变换后的潜变量的得分向量能捕获原数据矩阵的关键特征。在本发明方法中，这种关键特征为得分向量的自相关特征。根据自回归模型的思想，得分向量的自回归模型如下所示：

上式中，向量β_d为第d个自回归系数、e_d+1为自回归模型残差。一般来讲，自回归系数β₁，β₂，…，β_D的确定需要满足如下所示最小二乘目标函数：

由此可见，本发明方法所涉及的潜变量自回归模型需要确定两个未知项：投影向量w与自回归系数向量θ＝[β₁，β₂，…，β_D]^T，且相应的目标函数如下所示：

上式中，矩阵表示Kronecker内积，其具体的计算方式如下所示：

求解上述目标函数可以通过拉格朗日乘子法进行求解，首先利用拉格朗日乘子λ构造如下所示的拉格朗日函数J：

其次，计算函数J相对于向量w与向量θ的偏导数：

上式中，I_m与I_D分别表示m×m维与D×D维的单位矩阵，与的计算方式分别如公式(8)与公式(9)所示：

令公式(6)与公式(7)中的偏导数等于0，则可得到如下所示的等式：

令矩阵则公式(10)可简化成：这其实是一个广义特征值问题，投影向量w即为该广义特征值问题所对应的特征向量。对公式(10)中等号两边同时左乘向量w^T，即可得到如下所示等式关系：

因为所以公式(3)中的目标函数值等于特征值λ。由于公式(3)是求目标函数最小值，因此投影向量w为广义特征值问题：中最小正实数特征值所对应的特征向量，且需满足长度要求：

此外，公式(11)中定义了自回归系数向量θ与投影向量w之间的等式关系，即：

因此，在得到了投影向量w后即可根据公式(13)求解自回归系数向量θ。然而，投影向量w的计算需要已知自回归系数向量θ。由此可见，自回归系数向量θ与投影向量w的求解相互耦合，只能通过如下所示的往复迭代直至收敛的方式求得。

可先初始化自回归系数向量β＝[1，1，…1]^T；然后，根据广义特征值问题：中最小正实数特征值所对应的特征向量计算得到投影向量w；其次，根据公式(13)更新自回归系数向量θ；最后，根据更新后的自回归系数向量θ再次计算投影向量w。

一般而言，需要求解多个自回归系数向量θ与投影向量w。因此，在计算下一组自回归系数向量θ与投影向量w前，需要对矩阵实施处理：其中，

推理完潜变量自回归模型的算法后，即可建立潜变量自回归模型，如步骤(4)所示。

(4)设置潜变量的个数为K，根据如下所示步骤(4.1)至(4.6)求取得到K个潜变量的投影向量w₁，w₂，…，w_K、K个自回归系数向量θ₁，θ₂，…，θ_K、和K个载荷向量p₁，p₂，…，p_K。

(4.1)初始化k＝1与初始化向量θ_k＝[1，1，…，1]^T，即向量θ_k中所有元素都等于1。

(4.2)根据计算矩阵G后，再求解广义特征值问题：中最小正实数特征值对应的特征向量w。

(4.3)根据公式计算得到投影向量w_k后，在根据公式(13)更新自回归系数向量θ_k。

(4.4)判断向量θ_k是否收敛(判断的标准为：θ_k中的元素不再发生变化)？若否，则返回步骤(4.2)；若是，输出第k个自回归系数向量θ_k与投影向量w_k后执行步骤(4.5)。

(4.5)根据公式计算第k个潜变量的得分向量s_k后，根据公式计算载荷向量p_k。

(4.6)判断是否满足条件：k＜K？若是，则置k＝k+1后根据公式更新矩阵再返回步骤(4.2)；若否，则得到K个潜变量的投影向量w₁，w₂，…，w_K、K个自回归系数向量θ₁，θ₂，…，θ_K、以及K个载荷向量p₁，p₂，…，p_K。

(5)根据公式Θ＝P(W^TP)^-1计算投影变换矩阵Θ，其中矩阵W＝[w₁，w₂，…，w_K]，载荷矩阵P＝[p₁，p₂，…，p_K]。

(6)以更新后的矩阵为训练数据，利用主元分析算法对建立模型：其中U∈R^n×η、H∈R^m×η、和E∈R^n×m分别表示主元得分矩阵、主元载荷矩阵、和残差矩阵，η表示主元的个数。

(7)根据如下所示公式分别计算控制限ψ_lim、ξ_lim、和Q_lim

上式中，表示自由度为K、置信度为α的卡方分布所对应的值，F_{η，n-η，α}表示置信度为α、自由度分别为η与n-η的F分布所对应的值，表示自由度为h、置信度为α的卡方分布所对应的值，a和τ分别为Q统计量的估计均值和估计方差。

上述步骤(1)至步骤(7)为本发明方法的离线建模阶段，如下所示步骤(8)至步骤(13)为本发明方法的在线监测过程。

(8)采集最新采样时刻的样本数据x_t∈R^m×1，并找出其前D个采样时刻的样本数据x_t-1，x_t-2，…x_t-D，下标号t表示当前最新采样时刻。

(9)根据公式对x_t，x_t-1，x_t-2，…x_t-D分别实施标准化处理，对应得到其中j＝t，t-1，t-2，…，t-D。

(10)根据公式与公式分别计算潜变量的得分向量v_t，v_t-1，…，v_t-D与残差向量e。

(11)根据如下所示步骤(11.1)至步骤(11.3)消除潜变量的得分向量v_t中的自相关性。

(11.1)初始化k＝1。

(11.2)根据公式f_t，k＝v_t，k-V_kθ_k计算行向量f_t中的第k个元素f_t，k，其中v_t，k表示向量v_t中的第k个元素，V_k表示矩阵V＝[v_t-D ^T，v_t-D+1 ^T，…，v_t-1 ^T]中的第k行向量。

(11.3)判断是否满足条件：k＜K？若是，则置k＝k+1后返回步骤(11.2)；若否，则得到向量f_t。

(12)根据如下所示公式分别计算统计量ψ、ξ、和Q的具体数值：

上式中，矩阵Λ＝U^TU/(n-1)。

(13)判断是否满足条件ψ≤ψ_lim且ξ≤ξ_lim且Q≤Q_lim？若是，则当前样本采集自正常工况，返回步骤(8)继续监测下一时刻的样本数据；若否，则当前监测样本采集自故障工况。

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

本发明方法以建立潜变量自回归模型的目标，挖掘出动态自相关的潜变量，并给出相应满足最小二乘条件的自回归模型。通过该潜变量自回归模型不仅可以挖掘出原训练数据中的自相关特征，而且还可以消除潜变量自相关性的影响。因此，本发明方法与传统动态过程监测方法是存在显著不同的，而且模型的可解释性更强。可以说，本发明方法是一种更为优选的动态过程监测方法。

附图说明

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

图2为本发明方法中建立潜变量自相关模型的实施流程图。

图3为本发明方法与传统方法在监测TE过程故障上的监测详情对比图。

具体实施方式

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

如图1所示，本发明公开一种基于潜变量自回归模型的动态过程监测方法，具体实施方式如下所示。

表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	冷凝器冷却水流量

下面结合一个具体的工业过程的例子来说明本发明方法相对于现有方法的优越性与可靠性。该过程数据来自于美国田纳西-伊斯曼(TE)化工过程实验，原型是伊斯曼化工生产车间的一个实际工艺流程。目前，TE过程因其流程的复杂性，已作为一个标准实验平台被广泛用于故障检测研究。整个TE过程包括22个测量变量、12个操作变量、和19个成分测量变量。所采集的数据分为22组，其中包括1组正常工况下的数据集与21组故障数据。而在这些故障数据中，有16个是已知故障类型，如冷却水入口温度或进料成分的变化、阀门粘滞、反应动力学漂移等，还有5个故障类型是未知的。为了对该过程进行监测，选取如表1所示的33个过程变量，接下来结合该TE过程对本发明具体实施步骤进行详细的阐述。

首先，利用TE过程正常工况下的采样的960个数据建立故障检测模型，包括以下步骤：

步骤(1)：采集生产过程正常运行工况下的样本数据，组成训练数据矩阵X∈R^{960 ×33}，并计算矩阵X中各个列向量的均值μ₁，μ₂，…，μ₃₃以及标准差δ₁，δ₂，…，δ₃₃，对应组成均值向量μ＝[μ₁，μ₂，…，μ₃₃]^T与标准差对角矩阵Φ＝diag{δ₁，δ₂，…，δ₃₃}。

步骤(2)：根据公式对矩阵X实施标准化处理得到矩阵其中，Ξ∈R^960×33是由n个相同的均值向量μ组成的矩阵，即Ξ＝[μ，μ，…，μ]^T。

步骤(3)：记矩阵设置自相关阶数为D＝2后，再根据X_d＝[x_d，x_d+1，x_n+d+1]^T依次构造矩阵X₁，X₂，X₃，其中下标号d＝1，2，3，x_i为标准化后的第i个数据样本，i＝1，2，…，960。

步骤(4)：设置潜变量的个数为K＝6，根据如下所示步骤(4.1)至步骤(4.6)求取得到K个潜变量的投影向量w₁，w₂，…，w₆、K个自回归系数向量θ₁，θ₂，…，θ₆、和K个载荷向量p₁，p₂，…，p₆。图2展示了本发明方法中建立潜变量自相关模型的实施流程图，详细的实施过程具体为：

步骤(4.1)：初始化k＝1与初始化向量θ_k＝[1，1，…，1]^T，即向量θ_k中所有元素都等于1。

步骤(4.2)：根据如下所示公式(1)计算矩阵G后，求解广义特征值问题：中最小正实数特征值对应的特征向量w；

步骤(4.3)：根据公式计算得到投影向量w_k后，在根据公式(13)更新自回归系数向量θ_k。

步骤(4.4)：判断向量θ_k是否收敛？若否，则返回步骤(4.2)；若是，输出第k个自回归系数向量θ_k与投影向量w_k后执行步骤(4.5)。

步骤(4.5)：根据公式计算第k个潜变量的得分向量s_k后，再根据公式计算载荷向量p_k。

步骤(4.6)：判断是否满足条件：k＜K？若是，则置k＝k+1后根据公式更新矩阵再返回步骤(4.2)；若否，则得到K个潜变量的投影向量w₁，w₂，…，w₆、K个自回归系数向量θ₁，θ₂，…，θ₆、以及K个载荷向量p₁，p₂，…，p₆。

步骤(5)：根据公式Θ＝P(W^TP)^-1计算投影变换矩阵Θ，其中矩阵W＝[w₁，w₂，…，w₆]，载荷矩阵P＝[p₁，p₂，…，p₆]；

步骤(6)：以更新后的矩阵为训练数据，利用主元分析算法建立模型：具体的实施过程如步骤(6.1)至步骤(6.5)所示，其中U∈R^960×14、H∈R^33×14、和E∈R^960×33分别表示主元得分矩阵、主元载荷矩阵、和残差矩阵，η表示主元的个数。

步骤(6.1)：计算的协方差矩阵

步骤(6.2)：求解C所有特征值γ₁≥γ₂≥…≥γ_m所对应的特征向量h₁，h₂…，h_m。

步骤(6.3)：设置保留的主元个数η为满足如下所示条件的最小值，并将对应的η个特征向量组成主元载荷矩阵H＝[h₁，h₂…，h_η]。

步骤(6.4)：根据公式计算主元得分矩阵U∈R^960×14，那么PCA模型的残差矩阵为

步骤(7)：根据如下所示公式分别计算控制限ψ_lim、ξ_lim、和Q_lim

上式中，表示自由度为K、置信度为α的卡方分布所对应的值，F_{η，n-η，α}表示置信度为α、自由度分别为η与n-η的F分布所对应的值，表示自由度为h、置信度为α的卡方分布所对应的值，a和τ分别为Q统计量的估计均值和估计方差。F_{η，n-η，α}、和的具体数值可通过查询概率表获取。

当离线建模阶段完成后，采集TE过程中冷凝器冷却水进口温度阶跃故障工况下的960个采样数据，其中前160个数据是正常的，后800个数据才处于异常工况。在线故障监测的实施过程如下所示。

步骤(8)：采集最新采样时刻的样本数据x_t∈R^33×1，并找出其前D个采样时刻的样本数据x_t-1，x_t-2，…x_t-D。

步骤(9)：根据公式对x_t，x_t-1，x_t-2，…x_t-D分别实施标准化处理，对应得到其中j＝t，t-1，t-2，…，t-D。

步骤(10)：根据公式与公式分别计算潜变量的得分向量v_t，v_t-1，…，v_t-D与残差向量e。

步骤(11)：根据如下所示步骤(11.1)至步骤(11.3)消除变量的得分向量v_t中的自相关性。

步骤(11.1)初始化k＝1。

步骤(11.2)根据公式f_t，k＝v_t，k-V_kθ_k计算行向量f_t中的第k个元素f_t，k，其中v_t，k表示向量v_t中的第k个元素，V_k表示矩阵V＝[v_t-2 ^T，v_t-1 ^T]中的第k行向量。

步骤(11.3)判断是否满足条件：k＜K？若是，则置k＝k+1后返回步骤(11.2)；若否，则得到向量f_t。

步骤(12)根据如下所示公式分别计算统计量ψ、ξ、和Q的具体数值：

上式中，矩阵Λ＝U^TU/(n-1)。

步骤(13)判断是否满足条件ψ≤ψ_lim且ξ≤ξ_lim且Q≤Q_lim？若是，则当前样本采集自正常工况，返回步骤(8)继续监测下一时刻的样本数据；若否，则当前监测样本采集自故障工况。

如图3所示，传统动态PCA方法与本发明方法在监测TE过程冷凝器冷却水进口温度阶跃故障工况下的监测详情对比图。从图3中可以非常直接地看出本发明方法在故障发生后的漏报情况明显低于传统动态PCA方法。

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

Claims (1)

1.一种基于潜变量自回归模型的动态过程监测方法，其特征在于，包括以下步骤：

首先，本发明方法的离线建模阶段包括如下所示步骤(1)至步骤(7)；

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

步骤(2)：根据公式对矩阵X实施标准化处理得到矩阵其中，Ξ∈R^n×m是由n个相同的均值向量μ组成的矩阵，即Ξ＝[μ，μ，…μ]^T；

步骤(3)：记矩阵设置自相关阶数为D，根据X_d＝[x_d，x_d+1，…，x_n-D+d-1]^T依次构造矩阵X₁，X₂，…，X_D+1，其中下标号d＝1，2，…，D+1，x_t为标准化后的第i个数据样本，i＝1，2，…，n；

步骤(4)：设置潜变量的个数为K，根据如下所示步骤(4.1)至步骤(4.6)求取得到K个潜变量的投影向量w₁，w₂，…，w_K、K个自回归系数向量θ₁，θ₂，…，θ_K、和K个载荷向量p₁，p₂，…，p_K；

步骤(4.1)：初始化k＝1与初始化向量θ_k＝[1，1，…，1]^T，即向量θ_k中所有元素都等于1；

步骤(4.2)：根据如下所示公式(1)计算矩阵G后，再求解广义特征值问题：中最小正实数特征值对应的特征向量w；

上式(1)中，I_m表示m×m维的单位矩阵，符号表示Kronecker内积；

步骤(4.3)：根据公式计算得到投影向量w_k后，再根据如下所示公式更新自回归系数向量θ_k：

上式中，I_D表示D×D维的单位矩阵；

步骤(4.4)：判断向量θ_k是否收敛？若否，则返回步骤(4.2)；若是，则得到第k个自回归系数向量θ_k与投影向量w_k后执行步骤(4.5)；

步骤(4.5)：根据公式计算第k个潜变量的得分向量s_k后，再根据公式计算载荷向量p_k；

步骤(4.6)：判断是否满足条件：k＜K？若是，则置k＝k+1后根据公式更新矩阵再返回步骤(4.2)；若否，则得到K个潜变量的投影向量w₁，w₂，…，w_K、K个自回归系数向量θ₁，θ₂，…，θ_K、以及K个载荷向量p₁，p₂，…，p_K；

步骤(5)：根据公式Θ＝P(W^TP)^-1计算投影变换矩阵Θ，其中矩阵W＝[w₁，w₂，…，w_K]，载荷矩阵P＝[p₁，p₂，…，p_K]；

步骤(6)：利用主元分析算法为更新后的矩阵建立模型：具体的实施过程如步骤(6.1)至步骤(6.5)所示，其中U∈R^n×η、H∈R^m×η、和E∈R^n×m分别表示主元得分矩阵、主元载荷矩阵、和残差矩阵，η表示主元的个数；

步骤(6.1)：计算的协方差矩阵

步骤(6.2)：求解C所有特征值γ₁≥γ₂≥…≥γ_m所对应的特征向量h₁，h₂…，h_m；

步骤(6.3)：设置保留的主元个数η为满足如下所示条件的最小值，并将对应的η个特征向量组成主元载荷矩阵H＝[h₁，h₂…，h_η]；

步骤(6.4)：根据公式计算主元得分矩阵U∈R^n×η，那么残差矩阵为

步骤(7)：根据如下所示公式分别计算控制限ψ_lim、ξ_lim、和Q_lim

上式中，表示自由度为K、置信度为α的卡方分布所对应的值，F_{η，n-η，α}表示置信度为α、自由度分别为η与n-η的F分布所对应的值，表示自由度为h、置信度为α的卡方分布所对应的值，a和τ分别为Q统计量的估计均值和估计方差；

其次，完成上述离线建模阶段后，即可实施在线故障监测，具体包括如下所示步骤(8)至步骤(13)；

步骤(8)：采集最新采样时刻的样本数据x_t∈R^m×1，并找出其前D个采样时刻的样本数据x_t-1，x_t-2，…x_t-D，下标号t表示当前最新采样时刻；

步骤(9)：根据公式对x_t，x_t-1，x_t-2，…x_t-D分别实施标准化处理，对应得到其中j＝t，t-1，t-2，…，t-D；

步骤(10)：根据公式与公式分别计算潜变量的得分向量v_t，v_t-1，…，v_t-D与残差向量e；

步骤(11)：根据如下所示步骤(11.1)至步骤(11.3)消除潜变量的得分向量v_t中的自相关性；

步骤(11.1)初始化k＝1；

步骤(11.2)根据公式f_t，k＝v_t，k-V_kθ_k计算行向量f_t中的第k个元素f_t，k，其中v_t，k表示向量v_t中的第k个元素，V_k表示矩阵V＝[v_t-D ^T，v_t-D+1 ^T，…，v_t-1 ^T]中的第k行向量；

步骤(11.3)判断是否满足条件：k＜K？若是，则置k＝k+1后返回步骤(11.2)；若否，则得到向量f_t；

步骤(12)根据如下所示公式分别计算监测统计量ψ、ξ、和Q的具体数值：

上式中，矩阵Λ＝U^TU/(n-1)；

步骤(13)判断是否满足条件ψ≤ψ_lim且ξ≤ξ_lim且Q≤Q_lim？若是，则当前样本采集自正常工况，返回步骤(8)继续监测下一时刻的样本数据；若否，则当前监测样本采集自故障工况。