CN108572639B - 一种基于主成分自相关性剔除的动态过程监测方法 - Google Patents

Abstract

本发明公开一种基于主成分自相关性剔除的动态过程监测方法，旨在剔除传统动态主成分分析(DPCA)模型中主成分的自相关性，从而使相应过程监测模型具备更优越的动态过程监测效果。本发明方法首先为各测量变量引入在采样时间上的多个延时测量值；其次，建立DPCA模型；然后，利用延时测量值作为回归模型输入估计出对应的主成分信息；最后，利用剔除主成分自相关性的估计误差实施对动态过程的在线监测。相比于传统DPCA过程监测方法，本发明方法利用剔除主成分自相关性后的误差做为被监测对象，不仅不再受到自相关性的负面影响，而且误差的变化直接体现出主成分中自相关特征的异常变化。因此，本发明方法更适合于监测动态过程建模与监测。

Description

一种基于主成分自相关性剔除的动态过程监测方法

技术领域

本发明涉及一种数据驱动的过程监测方法，尤其涉及一种基于主成分自相关性剔除的动态过程监测方法。

背景技术

过程监测的目的在于及时而准确地发现故障，这对于保证安全生产与维持产品质量稳定是具有重要意义的。当前，过程监测的主流实施技术手段是数据驱动的方法，这主要得益于现代化工过程的大型化建设以及先进仪表与计算机技术的广泛应用，生产过程可以采集海量的数据。由于先进仪表技术的发展，采样时间间隔大为缩短，采样数据之间的时序自相关性是数据驱动的过程监测方法所必须考虑的一个问题。数据驱动的动态过程监测方法中最为典型的方法当属基于增广矩阵的动态主成分分析(Dynamic PrincipalComponent Analysis，DPCA)方法，其基本思想就是为各个训练样本数据引入延时测量值构成增广矩阵，然后对增广矩阵实施建模。使用增广矩阵可以同时将样本数据时序自相关性与变量之间的交叉相关性考虑进来，因此使用增广矩阵实施动态过程监测是最常见的技术手段。

然而，DPCA方法在建立过程监测模型时，是将自相关性与交叉相关性混淆在一起考虑。所提取的动态主成分信息虽然相互之间是不相关的，但是体现在采样时间上的自相关性却鲜有人提及。最近，有国外学者提出一种动态内部PCA(Dynamic-inner PCA，DiPCA)方法，可以剔除主成分的序列自相关性的干扰影响。但是，DiPCA方法针对的是静态PCA方法提取主成分的自相关性，未能在建模时将自相关性考虑进提取主成分信息的过程中。因此，DiPCA所能起到的改善动态过程监测性能是不明显的，这类改进思想还需进一步的深入研究。

一般而言，采样数据的自相关性主要表现在当前采样时刻的样本数据可以由前面几个时刻的样本数据预测出来。为了消除DPCA模型中动态主成分信息自相关性的负面影响，可以利用前面几个采样时刻的样本数据估计出当前采样时刻样本数据所对应的动态主成分信息，那么估计的动态主成分信息主要体现了真实动态主成分信息的自相关信息。两者之差不仅剔除了主成分信息的自相关性，而且还能保证在建模时同样如DPCA那般将自相关性与交叉相关性同时考虑进动态主成分信息的提取过程中，能起到改善动态过程监测效果的程度是可以预见的。

发明内容

本发明所要解决的主要技术问题是：如何剔除DPCA模型中主成分信息的自相关性，从而使相应过程监测模型具备更优越的动态过程监测效果。具体来讲，本发明方法首先为各测量变量引入在采样时间上的多个延时测量值；其次，按照传统DPCA建模的思想建立DPCA模型；然后，利用延时测量值作为回归模型输入估计出对应的主成分信息；最后，利用剔除主成分自相关性的估计误差实施对动态过程的在线监测。

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

(1)采集生产过程正常运行状态下的样本，组成训练数据集X∈R^N×m，并按照如下形式构造增广型矩阵X_a∈R^{(N-d)×m(d+1)}：

其中，N为训练样本数，m为过程测量变量数，R为实数集，R^N×m表示N×m维的实数矩阵，x_n∈R^1×m为第n个采样时刻的样本数据，下标号n＝1，2，…，N，d为引入的延时测量值的个数(一般取d＝1或2)。

(2)对矩阵X_a中每列实施标准化处理，得到均值为0，标准差为1的新数据矩阵

并将

表示成：

其中X_t∈R^(N-d)×m、X_t-1∈R^(N-d)×m、X_t-d∈R^(N-d)×m。

(3)利用主成分分析(Principal Component Analysis，PCA)算法为

建立相应的PCA模型：

其中T∈R^(N-d)×k、P∈R^m(d+1)×k、和E∈R^{(N-d)×m(d+1)}分别表示主成分矩阵、投影矩阵、和模型残差矩阵，上标号T表示矩阵或向量的转置，具体的实施过程如下所示：

①计算

的协方差矩阵

②求解C所有特征值γ₁≥γ₂≥…≥γ_m(d+1)所对应的特征向量p₁，p₂…，p_m(d+1)；

③设置保留的主成分个数k为满足如下所示条件的最小值，并将对应的k个特征向量组成载荷矩阵P＝[p₁，p₂…，p_k]；

④根据公式

与

分别计算主成分矩阵T∈R^(N-d)×k与模型残差矩阵E∈R^{(N-d)×m(d+1)}，那么相应的PCA模型为：

(4)将X_t-1，…，X_t-d组成一个新矩阵Y＝[X_t-1，…，X_t-d]∈R^(N-d)×md，并将主成分矩阵T按照公式T＝TΛ^-1/2进行归一化处理，其中Λ为γ₁，γ₂，…，γ_k组成的对角矩阵，

(5)利用偏最小二乘算法训练得到输入Y与输出T之间的回归模型，具体的实施过程如下所示：

①置i＝1与Z＝Y后，初始化向量u_i为T的第一列；

②根据公式w_i＝Zu_i/||Zu_i||、s_i＝Zw_i、q_i＝T^Ts_i/(s_i ^Ts_i)分别计算得到向量w_i、得分向量s_i、和系数向量q_i，其中||Zu_i||表示计算向量Zu_i的长度；

③根据公式u_new＝Tq_i/(q_i ^Tq_i)计算向量u_new；

④判断是否满足条件||u_i-u_new||＜10^-6；若否，则置u_i＝u_new后返回步骤②；若是，则执行⑤；

⑤根据公式p_i＝Z^Ts_i/(s_i ^Ts_i)计算第i个投影向量p_i，并保留向量p_i、向量w_i、和系数向量q_i；

⑥判断矩阵Y_i＝s_ip_i ^T中的最大元素是否大于0.01；若是，根据公式Z＝Z-s_ip_i ^T更新矩阵Z后执行步骤⑦；若否，则得到投影矩阵P＝[p₁，p₂，…，p_i]、矩阵W_i＝[w₁，w₂，…，w_i]、和系数矩阵Q＝[q₁，q₂，…，q_i]^T，并执行步骤⑧；

⑦判断i＜md；若是，则置i＝i+1后，返回步骤②；若否，则得到最终的投影矩阵P＝[p₁，p₂，…，p_i]、矩阵W_i＝[w₁，w₂，…，w_i]、和系数矩阵Q＝[q₁，q₂，…，q_i]^T；

⑧根据公式b_r＝W_r(P_r ^TW_r)^-1Q_r ^T计算在保留r个投影向量前提下的回归系数向量b_r，其中r＝1，2，…，i，P_r、W_r、和Q_r分别为矩阵P、W、和Q中前1至r列向量组成的矩阵；

⑨根据公式

计算输出T的预测值

后，计算保留不同个数投影向量所对应的模型预测误差

其中diag[]表示将矩阵对角线的元素单独组成向量的操作，sum{}表示求取向量中各元素之和；

⑩根据MSE₁，MSE₂，…，MSE_k的数值变化情况找出预测误差不再发生显著变化的收敛点位置，并将其下标号对应的具体数值作为偏最小二乘模型需保留的投影变量个数K，其中K≤i；

根据K确定最终的输入Y与输出T之间的回归模型，即：

上式中，回归矩阵Θ＝W_K(P_K ^TW_K)^-1Q_K ^T，

为回归模型误差。

(6)根据公式

计算

的协方差矩阵G，并根据如下所示公式计算监测统计量D的控制上限：

上式中，F(α，k，N-d-k)表示自由度为k与N-d-k的F分布在置信度α(一般取99％)下的取值。

上述步骤(1)～(6)为本发明方法的离线建模阶段，如下所示步骤(7)～(11)为在线动态过程监测的实施过程。

(7)收集新采样时刻的数据样本x_t∈R^1×m，引入其前d个采样时刻的样本得到增广向量x_a＝[x_t，x_t-1，…，x_t-d]，其中t表示当前采样时刻。

(8)对x_a实施与步骤(2)中同样的标准化处理得到

其中

分别对应于x_t，x_t-1，…，x_t-d标准化后的结果。

(9)根据公式

计算动态主成分向量z∈R^1×k，并根据公式

计算动态主成分向量的估计值

(10)根据公式z＝zΛ^-1/2归一化处理z后，计算剔除主成分自相关性后的误差

(11)根据公式D＝eG^-1e计算监测统计量D，并判断是否满足条件：D＞D_lim；若否，则当前样本为正常工况采样，返回步骤(7)继续实施对下一个样本数据的监测；若是，则当前采样数据来自故障工况。

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

首先，本发明方法是在传统DPCA方法基础上实施主成分自相关性剔除的，在提取主成分信息的过程中考虑到了采样数据的自相关性与交叉相关性；其次，本发明方法利用剔除主成分自相关性后的误差做为被监测对象，不仅使误差不再受到自相关性的负面影响，而且该误差的变化情况还能体现出主成分中自相关特征的异常变化。可以说，本发明方法更适合于动态过程建模与监测。

附图说明

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

图2为DPCA模型中主成分自相关性示意图

图3为误差中主成分自相关性剔除示意图。

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

具体实施方式

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

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

应用对象是来自于美国田纳西-伊斯曼(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	冷凝器冷却水流量

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

(1)收集生产过程正常运行状态下的数据样本，组成训练数据集X∈R^960×33，并构造增广型矩阵X_a∈R^958×99。

(2)对矩阵X_a中每列实施标准化处理，得到均值为0，标准差为1的新数据矩阵

并将

表示成：

(3)利用PCA算法为

建立相应的PCA模型：

(4)将X_t-1，X_t-2组成一个新矩阵Y＝[X_t-1，X_t-2]，并将主成分矩阵T按照公式T＝TΛ^-1/2进行归一化处理。

(5)利用偏最小二乘算法训练得到输入Y与输出T之间的回归模型：

为验证误差中主成分的自相关性已被剔除，首先将DPCA模型中主成分矩阵T的自相关性显示于图2中。从图2可以发现，主成分中存在很明显的自相关特性，尤其是第一个与第二个主成分，自相关性尤其明显。相比之下，从如图3所示的误差

的自相关性示意图中看看出，误差

中28个误差成分都不存在一阶及以上的自相关性。

(6)根据公式

计算

的协方差矩阵G，并根据如下所示公式计算监测统计量D的控制上限D_lim。

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

(7)收集新采样时刻的数据样本x_t∈R^1×33，引入其前d＝2个采样时刻的样本得到增广向量x_a＝[x_t，x_t-1，x_t-2]。

(8)对x_a实施与步骤(2)同样的标准化处理得到

(9)：根据公式

计算动态主成分向量z∈R^1×k，并根据公式

计算动态主成分向量的估计值

(10)根据公式z＝zΛ^-1/2归一化处理z后，计算剔除主成分自相关性后的误差

(11)根据公式D＝eC^-1e计算监测统计量D，并判断是否满足条件：D＞D_lim；若否，则当前样本为正常工况采样，返回步骤(7)继续实施对下一个样本数据的监测；若是，则当前采样数据自故障工况。

最后，将本发明方法与传统DPCA方法的过程监测详情对比于如图4中。值得指出的是，传统DPCA方法采用两个统计监测指标T²与Q实施过程监测，而本发明方法只使用一个监测统计指标D。从图4中可以发现，本发明方法对于该故障的监测效果要优越于传统DPCA方法，在故障发生后的故障漏报率显著低于传统DPCA方法的故障漏报率。

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

Claims (3)

1.一种基于主成分自相关性剔除的动态过程监测方法，其特征在于，包括以下步骤：

离线建模阶段的实施过程如下所示：

步骤(1)：采集生产过程正常运行状态下的样本，组成训练数据集X∈R^N×m，并按照如下形式构造增广型矩阵X_a∈R^{(N-d)×m(d+1)}：

上式中，N为训练样本数，m为过程测量变量数，R为实数集，R^N×m表示N×m维的实数矩阵，x_n∈R^1×m为第n个采样时刻的样本数据，下标号n＝1，2，…，N，d为引入的延时测量值的个数；步骤(2)：对矩阵X_a中每列实施标准化处理，得到均值为0，标准差为1的新数据矩阵

并将

表示成：

其中X_t∈R^(N-d)×m、X_t-1∈R^(N-d)×m、X_t-d∈R^(N-d)×m；

步骤(3)；利用主成分分析(Principal Component Analysis，PCA)算法为

建立相应的PCA模型：

其中T∈R^(N-d)×k、P∈R^m(d+1)×k和E∈R^{(N-d)×m(d+1)}分别表示主成分矩阵、投影矩阵和模型残差矩阵，上标号T表示矩阵或向量的转置；

步骤(4)：将X_t-1，…，X_t-d组成一个新矩阵Y＝[X_t-1，…，X_t-d]∈R^(N-d)×md，并将主成分矩阵T按照公式T＝TΛ^-1/2进行归一化处理，其中Λ＝T^TT/(N-d-1)；

步骤(5)：利用偏最小二乘算法训练得到输入Y与输出T之间的回归模型：

其中

为回归模型误差，Θ为回归矩阵；

步骤(6)：根据公式

计算

的协方差矩阵G，并根据如下所示公式计算监测统计量D的控制上限D_lim：

上式中，F(α，k，N-d-k)表示自由度为k与N-d-k的F分布在置信度α＝99％下的取值；在线过程监测的实施过程如下所示：

步骤(7)：收集新采样时刻的数据样本x_t∈R^1×m，引入其前d个采样时刻的样本得到增广向量x_a＝[x_t，x_t-1，…，x_t-d]，其中t表示当前采样时刻；

步骤(8)：对x_a实施与步骤(2)中同样的标准化处理得到

其中

分别对应于x_t，x_t-1，…，x_t-d标准化后的结果；

步骤(9)：根据公式

计算动态主成分向量z∈R^1×k，并根据公式

计算动态主成分向量的估计值

步骤(10)：根据公式z＝zΛ^-1/2归一化处理z后，计算剔除主成分自相关性后的误差

步骤(11)：根据公式D＝eG^-1e计算监测统计量D，并判断是否满足条件：D＞D_lim；若否，则当前样本为正常工况采样，返回步骤(7)继续实施对下一个样本数据的监测；若是，则当前采样数据来自故障工况。

2.根据权利要求1所述的一种基于主成分自相关性剔除的动态过程监测方法，其特征在于，所述步骤(3)中为

建立相应的PCA模型的具体实施过程如下所示：

①计算

的协方差矩阵

②求解C所有特征值γ₁≥γ₂≥…≥γ_m(d+1)所对应的特征向量p₁，p₂…，p_m(d+1)；

③设置保留的主成分个数k为满足如下所示条件的最小值，并将对应的k个特征向量组成载荷矩阵p＝[p₁，p₂…，p_k]；

④根据公式

与

分别计算主成分矩阵T∈R^(N-d)×k与模型残差矩阵E∈R^{(N-d)×m(d+1)}，那么相应的PCA模型为：

3.根据权利要求1所述的一种基于主成分自相关性剔除的动态过程监测方法，其特征在于，所述步骤(5)中利用偏最小二乘算法训练得到输入Y与输出T之间的回归模型的具体实施过程如下所示：

①置i＝1与Z＝Y后，初始化向量u_i为T的第一列；

②根据公式w_i＝Zu_i/||Zu_i||、s_i＝Zw_i、q_i＝T^Ts_i/(s_i ^Ts_i)分别计算得到向量w_i、得分向量s_i、和系数向量q_i，其中||Zu_i||表示计算向量Zu_i的长度；

③根据公式u_new＝Tq_i/(q_i ^Tq_i)计算向量u_new；

④判断是否满足条件||u_i-u_new||＜10^-6；若否，则置u_i＝u_new后返回步骤②；若是，则执行⑤；

⑤根据公式p_i＝Z^Ts_i/(s_i ^Ts_i)计算第i个投影向量p_i，并保留向量p_i、向量w_i、和系数向量q_i；

⑥判断矩阵Y_i＝s_ip_i ^T中的最大元素是否大于0.01？若是，根据公式Z＝Z-s_ip_i ^T更新矩阵Z后执行步骤⑦；若否，则得到投影矩阵p＝[p₁，p₂，…，p_i]、矩阵W_i＝[w₁，w₂，…，w_i]、和系数矩阵Q＝[q₁，q₂，…，q_i]^T，并执行步骤⑧；

⑦判断i＜md；若是，则置i＝i+1后，返回步骤②；若否，则得到最终的投影矩阵p＝[p₁，p₂…，p_i]、矩阵W_i＝[w_i，w₂，…，w_i]、和系数矩阵Q＝[q₁，q₂，…，q_i]^T；

⑧根据公式b_r＝W_r(P_r ^TW_r)^-1Q_r ^T计算在保留r个投影向量前提下的回归系数向量b_r，其中r＝1，2，…，i，P_r、W_r、和Qr分别为矩阵p、W、和Q中前1至r列向量组成的矩阵；

⑨根据公式

计算输出T的预测值

后，计算保留不同个数投影向量所对应的模型预测误差

其中diag[]表示将矩阵对角线的元素单独组成向量的操作，sum{}表示求取向量中各元素之和；

⑩根据MSE₁，MSE₂，…，MSE_k的数值变化情况找出预测误差不再发生显著变化的收敛点位置，并将其下标号对应的具体数值作为偏最小二乘模型需保留的投影变量个数K，其中K≤i；

根据K确定最终的输入Y与输出T之间的回归模型，即：

上式中，回归矩阵Θ＝W_K(P_K ^TW_K)^-1Q_K ^T，

为回归模型误差。

Images