CN109669414B - 一种基于自相关特征分解的动态过程监测方法 - Google Patents

Abstract

本发明公开一种基于自相关特征分解的动态过程监测方法，旨在从数据中分解出自相关特征成分，并基于此实施动态过程监测。本发明方法首先利用一种全新的自相关特征分解算法提取出潜藏的自相关特征；然后，将自相关特征成分中的自相关性过滤掉；最后，利用自相关特征成分误差以及模型误差实施动态过程监测。本发明方法的优势在于，首先本发明方法涉及的自相关分解算法是一种全新的特征提取算法，旨在挖掘出潜在特征成分的时间序列回归关系。其次，本发明方法在分解出训练数据中的自相关特征后，对各个特征成分逐一实施自相关性消除，并利用不包含显著自相关特征的残差实施动态过程监测。因此，本发明方法更适合于动态过程监测。

Description

一种基于自相关特征分解的动态过程监测方法

技术领域

本发明涉及一种数据驱动的过程监测方法，尤其涉及一种基于自相关特征分解的动态过程监测方法。

背景技术

随着传感器技术与计算技术在工业领域的广泛应用，对工业“大数据”的研究体现了现代工业过程数字化管理的水准。在这个背景下，利用采样数据实施过程运行状态的监测于是乎得到了工业界与学术界的高度重视。近年来，学术界与工业界都投入了大量的人力与物力研究以故障检测与诊断为核心任务的数据驱动的过程监测方法。在数据驱动的过程监测研究领域，诸多的数据挖掘算法都找到了用武之地。这其中当以统计过程监测为主，它是被研究得最多的方法。而作为一种经典的多变量统计分析算法，主成分分析(Principal Component Analysis，PCA)是最主流的实施技术手段，已发展出各式各样的过程监测方法技术体系。

由于先进仪表技术的发展，采样时间间隔大为缩短，采样数据之间的时序自相关性是数据驱动的过程监测方法所必须考虑的一个问题。数据驱动的动态过程监测方法中最为典型的方法当属基于增广矩阵的动态主成分分析(Dynamic Principal ComponentAnalysis，DPCA)方法，其基本思想就是为各个训练样本数据引入延时测量值构成增广矩阵，然后对增广矩阵实施建模。使用增广矩阵可以同时将样本数据时序自相关性与变量之间的交叉相关性考虑进来，因此使用增广矩阵实施动态过程监测是最常见的技术手段。然而，DPCA方法在建立过程监测模型时，是将自相关性与交叉相关性混淆在一起考虑。若是能将自相关性与交叉相关特征信息分开提取，那么相应模型就具备更强的可解释性。

一般而言，采样数据的自相关性主要体现在时间序列上，也就说当前时刻的样本数据与其前面几个时刻的采样数据是存在相互关联的。在某些情况下，故障工况下的采样数据改变的恰好是这个序列样本的自相关性。若不对自相关性特征进行单独的挖掘，这类故障可能就无法被识别出来。由此可见，对数据进行自相关特征挖掘在动态过程监测领域有着重要的意义与作用。一般而言，自相关特征的外在表现形式可简单描述成：当前时刻的样本数据可由前面几个采样时刻的样本数据预测出来。因此，如何挖掘这种时间序列关系是实施动态过程监测的关键。

发明内容

本发明所要解决的主要技术问题是：如何从数据中分解出自相关特征成分，并基于此实施动态过程监测。具体来讲，本发明方法首先利用一种新的特征提取算法通过投影变换提取出潜藏的自相关特征；然后，将自相关特征成分中的自相关性过滤掉；最后，利用自相关特征成分误差以及模型误差实施动态过程监测。

本发明解决上述技术问题所采用的技术方案为：一种基于自相关特征分解的动态过程监测方法，包括以下步骤：

(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)以

为新训练数据矩阵，设置自相关样本的个数A，按照如下所示公式构建A+1个数据矩阵X₁，X₂，…，X_A+1：

X_a＝[x_a，x_a+1，…，x_n-A+a-1]^T (2)

其中，下标号a＝1，2，…，A+1，x_i∈R^m×1，i＝1，2，…，n。

(4)依照本发明方法所涉及的自相关特征分解算法的实施过程，求解得到m个特征向量w₁，w₂，…，w_m，且保证各特征向量w_j都满足条件w_j ^TX_A+1 ^TX_A+1w_j＝1，其中j＝1，2，…，m。

本发明所涉及的自相关特征分解算法旨在搜寻投影变换向量w∈R^m×1从而使如下所示的目标函数达到最优，即：

约束条件：w^TX_A+1 ^TX_A+1w＝1

上式(3)中，矩阵Z＝[X_A，X_A-1，…，X₁]，列向量β＝[β₁，β₂，…，β_A]^T中各元素β₁，β₂，…，β_A为权重系数，

符号|| ||表示计算向量的长度

求解公式(3)中的带约束条件的最小化问题可通过拉格朗日乘子法完成，即首先定义如下所示的拉格朗日目标函数L：

然后，根据如下所示公式分别计算函数L相对于w与β的偏导数：

上两式中，λ为拉格朗日乘子法引入的乘子因子。根据拉格朗日乘子法的原则，需令偏导数

与

由此可得到如下所示的表达式：

很显然，上式(7)定义了一个广义特征值问题，上式(8)定义了向量w与向量β之间的关系。由于w与β都是待求解的向量，无法直接分开求解公式(7)与公式(8)中定义的等式关系。考虑到公式(7)与公式(8)这种相互耦合，且若是知道了向量w，则可通过公式(8)求解向量β，可采取如下所示的迭代过程计算向量w与β：

①初始化特征向量w＝[1，1，…，1]^T，并根据公式

更新特征向量w。

②根据公式β＝(Y^TY)^-1Y^TX_A+1w计算向量β，其中Y＝[X_Aw，X_A-1w，…，X₁w]。

③求解广义特征值问题：Φ₁w＝λΦ₂w，计算最小特征值所对应的特征向量w，其中Φ₁＝X_A+1 ^TX_A+1-X_A+1 ^TZ-Z^TX_A+1+Z^TZ^T，Φ₂＝X_A+1 ^TX_A+1，Z＝β₁X_A+β₂X_A-1+…+β_AX₁。

④根据公式

计算向量w_new，并判断是否满足收敛条件：||w-w_new||＜10^-8；若否，则置w＝w_new后，返回步骤②；若是，则输出收敛后的向量β与向量w。

⑤在广义特征值问题Φ₁w＝λΦ₂w中，最多可求得m个特征值，按大小进行升序排列为λ₁，λ₂，…，λ_m，对应的特征向量为w₁，w₂，…，w_m。值得注意的是，这里的特征向量w₁，w₂，…，w_m都是根据公式

处理后的向量。最小特征值λ₁及其对应的特征向量w₁为公式(7)的最优解，其余特征值及其对应的特征向量为次优解。

(5)设置自相关特征成分个数K后，可利用投影变换矩阵W＝[w₁，w₂，…，w_K]从矩阵

中分解出相应的自相关特征成分

并将其余的特征向量组成矩阵

由于训练数据中的自相关性特征可通过投影变换矩阵W分解出来，自相关特征成分矩阵S中各成分包含有显著的自相性，因此需要进一步描述各成分的自相关性。

(6)初始化k＝1，并根据如下所示步骤①至步骤③描述S中各成分的自相关性特征，从而得到K个最小二乘系数向量B₁，B₂，…，B_K。

①设置列向量s为矩阵S中的第k列，并将向量s中第a个元素至第n-A+a-1个元素组成列向量s_a，由于a可取值1，2，…，A+1，依次可得到A+1个列向量s₁，s₂，…，s_A+1。

②根据公式计算B_k＝(S₀ ^TS₀)^-1S₀ ^Ts_A+1计算最小二乘系数向量B_k∈R^A×1，其中矩阵S₀＝[s_A，s_A-1，…，s₁]。

③判断是否满足条件k＜K？若是，则置k＝k+1后，返回步骤①；若否，则保留得到的K个最小二乘系数向量Ｂ₁，B₂，…，B_K。

可以发现，上述步骤①至步骤③对K个自相关特征成分的时间序列关系实施了逐个逐个的描述，且时间序列上的关系是通过最小二乘法来描述的。

(7)根据公式

与公式

分别计算监测统计量D与Q的控制上限D_lim与Q_lim，其中

表示自由度为K的卡方分布在置信度α(一般取99％)下的取值，

表示自由度为m-K的卡方分布在置信度α下的取值，可通过查询概率统计表获取。

上述步骤(1)～(7)为本发明方法的离线建模阶段，在完成该阶段后，需保留的模型参数有：步骤(1)中均值向量μ，步骤(2)中的对角矩阵

步骤(5)中的投影变换矩阵W与矩阵

步骤(6)中的K个最小二乘系数向量B₁，B₂，…，B_K，和步骤(7)中的控制上限D_lim与Q_lim。

(8)采集当前时刻的样本数据x_t∈R^1×m，引入其前A个采样时刻的样本数据x_t-1，x_t-2，…，x_t-A，其中t表示当前采样时刻。

(9)根据公式

对样本数据x_t，x_t-1，…，x_t-A实施标准差处理对应得到样本数据

其中，下标号d指代t，t-1，…，t-A。

(10)根据公式

分解出自相关特征向量y_t，y_t-1，…，y_t-A。

(11)初始化k＝1后，根据如下所示步骤①至步骤③消除自相关特征向量y_t中的自相关特性。

①将向量y_t，y_t-1，…，y_t-A中第k个元素分别记录为

并将

组成一个行向量

②根据公式

计算向量y_t中第k个元素消除自相关性后的残差e_k。

③判断是否满足条件k＝K；若是，则置k＝k+1后返回步骤①；若否，则将得到的K个误差e₁，e₂，…，e_K组成一个行向量e＝[e₁，e₂，…，e_K]。

(12)根据公式

计算模型残差向量f后，再根据公式D＝ee^T与公式Q＝ff^T分别计算监测统计量D与Q。

(13)判断是否满足条件：D≤D_lim且Q≤Q_lim；若是，则当前样本采集自正常工况，返回步骤(8)继续实施对下一个样本数据的监测；若否，则当前采样数据为故障样本数据。

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

首先，本发明方法所涉及的自相关特征分解算法是一种全新的算法，旨在通过投影变换挖掘出潜在特征成分的时间序列回归关系。从这点上看，本发明方法提取的特征成分是充分考虑了时间序列上的自相关特征。其次，本发明方法在分解出训练数据中的自相关特征后，对各个特征成分逐一实施自相关性消除，并利用不包含显著自相关特征的残差实施动态过程监测。可以说，本发明方法更适合于动态过程建模与监测。

附图说明

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

图2为第1个自相关特征成分至第10个自相关特征成分的自相关性示意图。

图3为第11个自相关特征成分至第16个自相关特征成分的自相关性示意图。

图4为TE过程物料进口温度故障的监测详情对比图。

具体实施方式

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

本发明公开一种基于自相关特征分解的动态过程监测方法，其离线建模的实施流程如图1所示。下面结合一个具体的工业过程的例子来说明本发明方法的具体实施过程，以及相对于现有方法的优越性。

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

表1：TE过程监测变量。

首先，利用TE过程正常工况下采样的960个样本数据离线建立动态过程监测模型，包括以下步骤：

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

步骤(2)：根据公式

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

步骤(3)：以

为新训练数据矩阵，设置自相关样本的个数A＝2，根据公式X_a＝[x_a，x_a+1，…，x_n-A+a-1]^T构建A+1个数据矩阵X₁，X₂，…，X_A+1，其中a＝1，2，3。

步骤(4)：根据自相关特征分解算法的实施过程，求解得到33个特征向量w₁，w₂，…，w₃₃，且保证各特征向量w_j都满足条件w_j ^TX_A+1 ^TX_A+1w_j＝1，其中j＝1，2，…，33。

步骤(5)：设置自相关特征成分个数K＝10后，可利用投影变换矩阵W＝[w₁，w₂，…，w₁₀]从训练数据

中分解出相应的自相关特征成分

并将其余的特征向量组成矩阵

为验证分解出的自相关特征成分S中存在显著的自相关性，现将第1个自相关特征成分至第10个自相关特征成分的自相关性显示于图2中。从图2中可以发现，矩阵S中各列向量存在较强的自相关性。

此外，额外将第11个自相关特征成分至第16个自相关特征成分的自相关性显示于图3中。从图3中可以看出，所列举的6个自相关特征成分不存在显著的自相关特性。因此，本实施案例中设置K＝10个自相关特征成分是合理的。

步骤(6)：初始化k＝1，并根据如下所示步骤(6.1)至步骤(6.3)计算得到K个最小二乘系数向量B₁，B₂，…，B₁₀。

步骤(6.1)：设置列向量s为矩阵S中的第k列，并将向量s中第a个元素至第n-A+a-1个元素组成列向量s_a，由于a可取值1，2，3，依次可得到3个列向量s₁，s₂，s₃。

步骤(6.2)：根据公式计算B_k＝(S₀ ^TS₀)^-1S₀ ^Ts₃计算最小二乘系数向量B_k∈R^2×1，其中矩阵S₀＝[s₂，s₁]。

步骤(6.3)：判断是否满足条件k＜10；若是，则置k＝k+1后，返回步骤(6.1)；若否，则保留得到的K个最小二乘系数向量B₁，B₂，…，B₁₀。

步骤(7)：根据公式

与公式

分别计算监测统计量D与Q的控制上限D_lim与Q_lim。

其次，采集TE过程物料进口温度故障条件下的测试数据集，实施在线过程监测。值得指出的是，该测试数据集前160个样本数据采集自正常工况，故障工况从161个时刻起引入。

步骤(8)：采集当前时刻的样本数据x_t∈R^1×m，引入其前A＝2个采样时刻的样本数据x_t-1，x_t-2，其中t表示当前采样时刻。

步骤(9)：根据公式

对样本数据x_t，x_t-1，x_t-2实施标准差处理对应得到样本数据

其中，下标号d指代t，t-1，t-2。

步骤(10)：根据公式

分解出自相关特征向量y_t，y_t-1，y_t-2。

步骤(11)：根据最小二乘系数向量B₁，B₂，…，B_K逐一消除自相关特性以得到残差向量e。

步骤(12)：根据公式

计算模型残差向量f后，再根据公式D＝ee^T与公式Q＝ff^T分别计算监测统计指标D与Q。

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

最后，将本发明方法与传统DPCA方法的过程监测详情对比于如图4中。从图4中可以发现，本发明方法对于该故障的监测效果要优越于传统DPCA方法，在故障发生后的故障漏报率显著低于传统DPCA方法的故障漏报率。

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

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)：以

为新训练数据矩阵，设置自相关样本的个数为A，按照如下所示公式构建A+1个数据矩阵X₁，X₂，…，X_A+1：

X_a＝[x_a，x_a+1，…，x_n-A+a-1]^T (2)

其中，下标号a＝1，2，…，A+1，x_i∈R^m×1，i＝1，2，…，n；

步骤(4)：根据如下所示的步骤(4.1)至步骤(4.5)求解得到m个特征向量w₁，w₂，…，w_m，且保证各特征向量w_j都满足条件w_j ^TX_A+1 ^TX_A+1w_j＝1，其中j＝1，2，…，m；

步骤(4.1)：初始化特征向量w＝[1，1，…，1]^T，并根据公式

更新特征向量w；

步骤(4.2)：根据公式β＝(Y^TY)^-1Y^TX_A+1w计算向量β，其中矩阵Y＝[X_Aw，X_A-1w，…，X₁w]；

步骤(4.3)：求解广义特征值问题：Φ₁w＝λΦ₂w，计算最小特征值所对应的特征向量w，其中，矩阵Φ₁＝X_A+1 ^TX_A+1-X_A+1 ^TZ-Z^TX_A+1+Z^TZ^T，矩阵Z＝β₁X_A+β₂X_A-1+…+β_AX₁，矩阵Φ₂＝X_{A+ 1} ^TX_A+1；

步骤(4.4)：根据公式

计算向量w_new，并判断是否满足收敛条件：||w-w_new||＜10^-8；若否，则置w＝w_new后，返回步骤(4.2)；若是，则输出收敛后的向量β；

步骤(4.5)：求解广义特征值问题：Φ₁w＝λΦ₂w，可得到m个特征值，按大小进行升序排列为λ₁≤λ₂≤…≤λ_m，对应的特征向量为w₁，w₂，…，w_m，特征向量w₁，w₂，…，w_m都需根据公式

进行处理；

步骤(5)：设置自相关特征成分个数K后，可利用投影变换矩阵W＝[w₁，w₂，…，w_K]从矩阵

中分解出相应的自相关特征成分

并将其余的特征向量组成矩阵

步骤(6)：初始化k＝1，并根据如下所示步骤(6.1)至步骤(6.3)计算得到K个最小二乘系数向量B₁，B₂，…，B_K；

步骤(6.1)：设置列向量s为矩阵S中的第k列，并将向量s中第a个元素至第n-A+a-1个元素组成列向量s_a，由于a可取值1，2，…，A+1，依次可得到A+1个列向量s₁，s₂，…，s_A+1；

步骤(6.2)：根据公式计算B_k＝(S₀ ^TS₀)^-1S₀ ^Ts_A+1计算最小二乘系数向量B_k∈R^A×1，其中矩阵S₀＝[s_A，s_A-1，…，s₁]；

步骤(6.3)：判断是否满足条件k＜K；若是，则置k＝k+1后，返回步骤(6.1)；若否，则保留得到的K个最小二乘系数向量B₁，B₂，…，B_K；

步骤(7)：根据公式

与公式

分别计算监测统计量D与Q的控制上限D_lim与Q_lim，其中

表示自由度为K的卡方分布在置信度α下的取值，

表示自由度为m-K的卡方分布在置信度α下的取值，可通过查询概率统计表获取；

离线建模阶段包括上述步骤(1)至步骤(7)，需保留如下所示的模型参数以便实施在线监测时调用：步骤(1)中均值向量μ，步骤(2)中的对角矩阵

步骤(5)中的投影变换矩阵W与矩阵

步骤(6)中的K个最小二乘系数向量B₁，B₂，…，B_K，和步骤(7)中的控制上限D_lim与Q_lim；

步骤(8)：采集当前时刻的样本数据x_t∈R^1×m，引入其前A个采样时刻的样本数据x_t-1，x_t-2，…，x_t-A，其中t表示当前采样时刻；

步骤(9)：根据公式

对样本数据x_t，x_t-1，…，x_t-A实施标准差处理对应得到样本数据

其中，下标号d指代t，t-1，…，t-A；

步骤(10)：根据公式

分解出自相关特征向量y_t，y_t-1，…，y_t-A；

步骤(11)：初始化k＝1后，根据如下所示步骤(11.1)至步骤(11.3)消除自相关特性以得到残差向量e；

步骤(11.1)：将向量y_t，y_t-1，…，y_t-A中第k个元素分别记录为

并将其组成一个行向量

步骤(11.2)：根据公式

计算向量y_t中第k个元素消除自相关性后的残差e_k；

步骤(11.3)：判断是否满足条件k＝K；若是，则置k＝k+1后返回步骤(11.1)；若否，则将得到的K个误差e₁，e₂，…，e_K组成一个行向量e＝[e₁，e₂，…，e_K]；

步骤(12)：根据公式

计算模型残差向量f后，再根据公式D＝ee^T与公式Q＝ff^T分别计算监测统计量D与Q；

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

Images