CN104699075A - 一种基于不等长时段自动有序划分的过程监测方法 - Google Patents

Abstract

本发明公开了一种基于不等长时段自动有序划分的过程监测方法，该方法使用与监测性能相关的过程相似性衡量指标，有效捕捉了具有不等长操作周期的批次过程的不规则过程特性，从而将不等长批次过程自动有序地划分为多个不等长时段。其后，基于不等长子时段的划分结果建立了监测模型，并针对不同时间区域采取不同的监测策略，有效监测了过程的运行状态，区分了故障与时段转换。本发明简单易于实施，不仅有利于对具体过程特性的了解，而且增强了实际在线过程监测的可靠性和可信度，有助于工业工程师对过程运行状态做出准确判断，及时发现故障，从而保证实际生产的安全可靠运行和产品的高质量追求。

Description

一种基于不等长时段自动有序划分的过程监测方法

技术领域

本发明属于批次过程统计监测领域，特别是涉及一种基于不等长时段自动有序划分的过程监测方法。

背景技术

作为工业生产中一种重要的生产方式，批次过程与人们的生活息息相关，已被广泛应用于精细化工、生物制药、食品、聚合物反应、金属加工等领域。近年来，随着现代社会对多品种、多规格和高质量产品更迫切的市场需求，工业生产更加倚重于生产小批量、高附加值产品的间歇过程。批次生产的安全可靠运行以及产品的高质量追求已成为人们关注的焦点。基于数据的多元统计分析技术因其只需要正常工况下的过程数据来建立模型，同时它们在处理高维、高度耦合数据时具有独特的优势，越来越受到研究人员和现场工程师的青睐。批次过程的统计建模、在线监测、故障诊断及质量预测已成为广泛的研究课题。

传统的过程监测方法均基于一个理想化的假设：批次过程具有相同的操作周期且每个批次轨迹都是同步的。然而，在实际工业生产中，由于生产条件以及生产目标的不同，批次轨迹往往是不等长的。考虑到批次过程具有多时段特性，那么对于多时段批次过程来说，它的时段在批次方向也是不等长的。不等长的时段数据给过程统计分析和建模带来了困难。为了解决不等长问题，人们做出了一系列的研究，这些研究主要可以分为两大类：直接的信号同步方法和基于不等长时段划分的方法。

直接的信号同步方法指使用一些信号处理的方法将批次进行同步。如最简单的将长批次截断至与短批次长度相等的方法以及将短批次数据补充至与长批次相等的方法。前者只适用于不等长批次共同部分的轨迹是重合的情况，后者则需要大量的批次数据而不适用于实际工业生产过程。另一种方法是基于指示变量的方法，即用一个能指示过程开始与结束的变量代替时间指标，运用重采样与插值的手段对批次进行等长化处理。然而，重采样和插值往往改变了批次过程原有的过程特性，并且指示变量往往无法指示各个时段的开始与结束。此外，基于动态时间规整的方法也被用于批次轨迹的同步，该类方法虽然能将批次很好地等长化，但是同样也会改变原始的过程特性，导致过程监测能力下降。总的来说，直接的信号同步的方法没有考虑到批次过程的多时段特性，没有从本质上解决不等长批次的问题。

基于不等长时段划分的方法克服了直接信号同步方法的缺点，没有硬性将批次进行等长处理，而是先将不等长的时段划分出来，再针对每个时段建立监测模型，进行在线监测。然而，前人均是采用聚类的方法对不等长时段进行划分，该类方法有以下几个缺点：(1)过程相似性衡量指标与监测性能没有直接关系，导致监测模型精度欠缺；(2)没有考虑批次过程运行在时间上的有序性，从而导致时段划分结果呈现一种“混乱”状态，需要大量的后续处理；(3)采用时间片作为基本的分析单元，单个批次的过程特性容易被掩盖，使得时段划分结果不够精确，进而影响监测模型精度以及在线监测性能。

本发明的内容深入考虑了不等长批次过程的不规则的时变过程特性和实际过程运行的时序性以及时段划分结果对于之后监测性能的影响，提出了一种新的不等长子时段自动划分及过程监测方法。到目前为止，尚未见到与本发明相关的研究报道。

发明内容

本发明的目的在于针对现有针对不等长批次生产过程的时段划分及过程监测技术的不足，提供一种基于不等长时段自动有序划分的过程监测方法。该方法通过依次分析每个批次的潜在过程特性，有效捕捉了具有不等长操作周期的批次过程的不规则过程特性，从而将不等长批次过程自动有序地划分为多个不等长子时段；并基于不等长子时段建立了监测模型，针对不同时间区域采用不同的监测策略，有效监测了过程的运行状态，区分了故障与时段转换，最终可应用于实际工业生产现场，确保批次生产的安全可靠运行以及产品的高质量追求。

本发明的目的是通过以下技术方案来实现的：一种自动的步进式有序时段划分方法，该方法包括以下步骤：

(1)获取过程分析数据：设每一个不等长批次操作具有J个测量变量和K_i个采样点，则每一个测量批次可得到一个K_i×J的矩阵X_i(K_i×J)，重复I批次的测量步骤后，得到的数据表述为一个不规则的三维矩阵X，其中，i表示批次，测量变量为批次运行过程中可被测量的状态参数；

(2)数据预处理：将三维矩阵X按照采集变量方向展开，即保持三维矩阵X变量维度不变，依次将X_i(K_i×J)从上至下排列得到二维矩阵

设二维矩阵X内任意一点的变量为x_i,j，对该变量进行减均值、除以标准差的标准化处理，其中，下标i代表批次，j代表变量，标准化处理的计算公式如下：

$x_{i,j}=\frac{x_{i,j}-{\stackrel{\‾}{x}}_j}{s_j}---\left(1\right)$

其中：是二维矩阵X任一列的均值，s_j是二维矩阵X相应列的标准差，

$\begin{array}{c}{\stackrel{\‾}{x}}_j=\frac{1}{\mathrm{Sum}}\underset{i=1}{\overset{\mathrm{Sum}}{\mathrm{\Σ}}}{x}_{i,j}\\ s_j=\sqrt{\underset{i=1}{\overset{\mathrm{Sum}}{\mathrm{\Σ}}}{(x_{i,j}-{\stackrel{\‾}{x}}_j)}^2/(\mathrm{Sum}-1)}\\ \mathrm{Sum}=\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{K}_i\end{array}---\left(2\right)$

将已标准化的重新排列分离出时间片矩阵X_k(I_k×J)，即其中，k表示采样时间指标，I_k表示每个时间片所包含的批次操作次数；

(3)建立初始时段模型，该步骤包括以下子步骤：

(3.1)对步骤2标准化处理后得到的时间片矩阵X₁(I₁×J)表示为X_v,1(I₁×J)进行PCA分解，建立初始时段模型，其中PCA分解公式如下：

$X_{v,1}=T^{\′}{\left(P_{v,1}^{\′}\right)}^T=\underset{r=1}{\overset{J}{\mathrm{\Σ}}}{t}_r{p}_r;---\left(3\right)$

其中：下标v表示时段模型按变量方式展开，t_r为正交的主元向量，p_r为正交归一化的负载向量，r表示不同的PCA分解方向，上标T表示矩阵的转置；T′(I₁×J)代表保留全部主元的得分矩阵，P′_v,1(J×J)代表对应的负载矩阵；

(3.2)选取主元个数，将公式(3)重新表述成如下形式：

$X_{v,1}={{\mathrm{TP}}_{v,1}}^T+E=\underset{r=1}{\overset{R}{\mathrm{\Σ}}}{t}_r{p}_r+E---\left(4\right)$

其中：r表示不同的PCA分解方向，R为主元个数；T与P_v,1(J×R)分别保留R个主元后的得分矩阵和负载矩阵，E为残差矩阵；

(3.3)计算SPE指标，确定初始时段模型控制限：

SPE_i＝e_i ^Te_i   (5)

其中，下标i表示时间片中的批次，e_i是第i批次的残差列向量，即e_i为E的第i行向量的转置向量；确定出初始时段模型的控制限Ctr_v,1＝α*max(SPE_i)，其中，α是缓和因子；

(4)衡量变量相关性：调用初始时段模型依次计算下一个时间片每个批次的SPE指标：

t_i＝P_v,1x_i

$\begin{array}{c}{e}_i=x_i-P_{v,1}^T{t}_i\\ {\mathrm{SPE}}_i=e_i^T{e}_i\end{array}---\left(6\right)$

其中，下标i表示时间片中的批次，x_i为时间片中的每一个批次，e_i是第i批次的残差列向量，即e_i为E的第i行向量的转置向量，t_i表示第i批次的主元得分向量；如果发现SPE_i>Ctrl_v,1，则继续调用初始时段模型计算该批次下两个时刻的SPE指标，如果有连续三个样本呈现SPE_i>Ctrl_v,1，那么该批次变量相关性与初始时段模型所表征的特性差别显著，记录该批次第一次超限时间为k^*，则该批次k^*时刻之前的数据可认为是一个子时段；如果没有连续三个样本呈现SPE_i>Ctrl_v,1，则认为该批次仍处于当前时段没有进入下一时段；

(5)时段模型更新：将下一时间片中仍处于当前时段的批次加入初始时段模型的建模数据中并按变量展开，得到对重新按步骤3进行PCA建模，得到更新后的模型以及控制限Ctrl_v,2；

(6)确定所有批次的第一时段：重复步骤4～5，直至所有批次的第一时段均被划分出来；

(7)过程数据更新，确定所有时段：移除已划分出的不等长子时段数据，将剩余数据重新对齐作为下一次时段划分的数据；重复步骤3～6直到所有不等长时段均被划分出来，至此，每个不等长时段的最短批次长度k_s和最长批次长度k_l均被确定出来；

(8)基于不等长时段划分结果的统计建模，该步骤包括以下子步骤：

(8.1)将每个时段内k_s+1至k_l时刻的时间片X_k(I_k×J)按变量展开构造泛化时间片则k_s+1至k_l时刻的时间片X_k(I_k×J)均可由泛化时间片表示；

(8.2)根据步骤7的时段划分结果，每个时段内的时间片X_k(I_k×J)及泛化时间片按公式(1)进行标准化处理，获得时段内每个时刻的均值和标准差；

(8.3)将每个时段内1至k_s时刻的时间片X_k(I_k×J)和泛化时间片按变量展开得到对进行PCA建模：

T_c＝X_cP_c,s

$\begin{array}{c}{X}_c=T_c{P}_{c,s}^T+E_c=X_c{P}_{c,s}{P}_{c,s}^T+E_c\\ {\hat{X}}_c=T_c{P}_{c,s}^T=X_c{P}_{c,s}{P}_{c,s}^T\\ E_c=X_c-{\hat{X}}_c=X_c{P}_{c,e}{P}_{c,e}^T\end{array}---\left(7\right)$

其中，下标c表示时段指标，s表示主元空间，e表示残差空间；Kc表示时段所包含的采样时刻个数，P_c,s(J×R_c)是负载矩阵，是主元得分矩阵，R_c是所保留的主元个数，是由主元得分重构的时段数据矩阵，E_c是残差矩阵，P_c,e(J×R_r)(R_r＝J-R_c)反映了时段内次要的波动方向；

(9)计算实时监测统计指标并确定其控制限，该步骤包括以下子步骤：

(9.1)按公式(8)计算Hotelling-T²统计指标和SPE统计量；

$\begin{array}{c}{T_k}^2={(t_k-{\stackrel{\‾}{t}}_k)}^T{\mathrm{\Σ}}_c^{-1}(t_k-{\stackrel{\‾}{t}}_k)\\ {\mathrm{SPE}}_k=e_k^T{e}_k\end{array}---\left(8\right)$

其中，t_k(R_c×1)是每个批次每个时刻的主元得分向量，它可以从中分离得到；∑_c是由构成的协方差矩阵；表示所有t_k(R_c×1)的均值向量；e_k(J×1)是每个批次每个时刻的残差向量，它可以从中分离得到；

(9.2)确定控制限：根据过程数据服从多维高斯分布可知SPE统计量服从加权χ²分布，从而确定SPE统计量的控制限；将每个时刻的Hotelling-T²统计指标从小到大排列，选取第0.95I_k或者0.90I_k个Hotelling-T²统计指标作为其控制限；

(10)基于时段模型的在线过程监测：基于步骤7划分得到的不等长时段、步骤8建立的时段模型和步骤9所得的两个监测统计量可在线监测间歇过程的运行状态；该步骤由以下子步骤来实现：

(10.1)采集新测量数据：在线监测时，采集新的过程测量数据x_new(J×1)，其中，下标new代表新样本，J为测量变量，与步骤1中的测量变量相同；对于不同时刻的新数据，需要采取不同的监测策略；

(10.2)对于位于时段内[1,k_s]时刻的新数据，按如下子步骤进行在线监测：

(10.2.1)根据步骤8中获得的时段内每个时刻的均值和标准差，根据过程时间的指示调用当前时段相应时刻的均值和标准差，对采集到的过程测量数据如公式(1)中所示进行标准化预处理；

(10.2.2)计算新监测统计量：数据预处理后，调用对应该新采样时刻所在时段的模型P_c(J×R_c)，按照公式(9)估计对应的Hotelling-T²与SPE两个监测统计指标：

t_new ^T＝x_new ^TP_c

$\begin{array}{c}{\hat{x}}_{\mathrm{new}}=P_c{t}_{\mathrm{new}}\\ T_{\mathrm{new}}^2={(t_{\mathrm{new}}-{\stackrel{\‾}{t}}_k)}^T{\mathrm{\Σ}}_c^{-1}(t_{\mathrm{new}}-{\stackrel{\‾}{t}}_k)\\ {\mathrm{SPE}}_{\mathrm{new}}={(x_{\mathrm{new}}-{\hat{x}}_{\mathrm{new}})}^T(x_{\mathrm{new}}-{\hat{x}}_{\mathrm{new}})\end{array}---\left(9\right)$

其中，x_new是新的过程测量数据，t_new是新数据的主元得分向量，是重构的过程测量数据，是根据训练数据获得的主元得分的均值向量，∑_c是由构成的协方差矩阵；

(10.2.3)在线判断过程运行状态：实时比较两个监测指标与其各自的统计控制限，若两个监测指标都位于统计控制限之内，表明过程运行正常；若有一个以上监测指标超出正常控制限，表明过程有异常状况发生。

(10.3)对于位于时刻[k_s+1,k_l]的新数据，其时段隶属不能根据时刻来判断，则按如下子步骤进行在线监测：

(10.3.1)根据步骤8中获得的时段内每个时刻的均值和标准差，先根据过程时间的指示调用当前时段的相应时刻的均值和标准差，对采集到的过程测量数据如公式(1)中所示进行标准化预处理；

(10.3.2)计算新监测统计量：数据预处理后，调用对应该新采样时刻所在时段的模型P_c(J×R_c)，按照公式(9)计算得到主元得分，估计残差及其对应的Hotelling-T²与SPE两个监测统计指标；

(10.3.3)在线判断过程运行状态：实时比较两个监测指标与其各自的统计控制限，若两个监测指标都位于统计控制限之内，表明过程位于当前时段，且运行正常；若有一个以上监测指标超出正常控制限，表明过程有异常状况发生，此时需要判别过程是进入下一时段还是有故障发生；

(10.3.4)根据步骤8中所获得的时段内每个时刻的均值和标准差，调用下一时段相应时刻的均值和标准差对新数据进行标准化预处理；

(10.3.5)计算新监测统计量：数据预处理后，调用下一时段的模型P_c+1(J×R_c+1)，按照公式(9)计算得到主元得分，估计残差及其对应的Hotelling-T²与SPE两个监测统计指标；

(10.3.6)在线判断过程运行状态：实时比较两个监测指标与其各自的统计控制限，若两个监测指标都位于统计控制限之内，表明过程从当前时段进入下一时段，且运行正常；若有一个以上监测指标超出正常控制限，表明过程有异常状况发生。

进一步地，所述步骤1中，所述测量变量为以下9个：喷嘴温度、喷嘴压力、螺杆行程、螺杆速度、注射压力、塑化压力、压力阀开度、流量阀开度、模腔压力。

与现有技术相比，本发明的有益效果是：本发明为无过程先验知识条件下的不等长多阶段过程时段划分、建模和监测提供了新的研究思路。所提出的不等长子时段自动划分及过程监测方法可以应用于一类具有不等长多运行时段的批次生产过程，并能够提高建立的子时段模型的精度，改善后续过程在线监测的性能。所提出的方法在注塑工业过程中做了详细的实验研究，获得了成功应用，增强了对具体过程运行特性的了解，提高了过程监测过程的监测效率和故障检测结果的准确性，最终可应用于实际工业生产现场，确保间歇生产的安全可靠运行以及产品的高质量追求。

附图说明

图1是本发明所述自动的基于不等长时段自动有序划分的过程监测方法的流程图；

图2是本发明所述的三维数据展开模型示意图；

图3是本发明具体实施例中的本发明方法不等长时段划分结果与真实不等长时段结果对比图；

图4是本发明具体实施例中的本发明方法不等长时段划分结果与真实不等长时段结果对比图；

图5是α＝1.05的时段划分结果图；

图6是α＝1.15的时段划分结果图；

图7是α＝1.20的时段划分结果图；

图8是α＝2.00的时段划分结果图；

图9是基于聚类的方法第一次时段聚类结果图，(a)为θ＝0.1时聚类结果图；(b)为θ＝0.5时聚类结果图，(c)为θ＝1.0时聚类结果图，(d)为θ＝1.5时聚类结果图；

图10是时段转换监测图；

图11是保压段故障监测图。

具体实施方式

下面结合附图及具体实例，对本发明做进一步说明。

实施例1

注塑成型过程是典型的多时段间歇生产过程，一般由注射、保压、冷却三个阶段构成。此外，塑化过程在冷却初始阶段完成。具体地说，在注射阶段，液压系统推动螺杆将塑料粘流体注入模腔中，直至模腔被流体充满。过程处于保压阶段时，仍有少量的粘流体被高压挤进模腔中，以补偿塑料粘流体在冷却和塑化时造成的体积收缩。保压阶段一直持续到模腔的浇口冻结，过程进入冷却段。当螺杆头部熔料逐渐增多，达到一定的注射量后，螺杆停止后退和转动，这段时间的过程状态也称塑化段。随着模腔中熔料继续冷却，塑件完全固化，模具打开，塑件被顶出，从而完成一个工作循环。

如图1所示，本发明一种基于不等长时段自动有序划分的过程监测方法，包括以下步骤：

(1)获取过程分析数据：设每一个注塑成型过程中不等长批次操作具有J个测量变量和K_i个采样点，则每一个测量批次可得到一个K_i×J的矩阵X_i(K_i×J)，重复I批次的测量步骤后，得到的数据表述为一个不规则的三维矩阵X，其中，i表示批次，测量变量为注塑成型过程中可被测量的状态参数，包括喷嘴温度、喷嘴压力、螺杆行程、螺杆速度、注射压力、塑化压力、压力阀开度、流量阀开度、模腔压力等；

(2)数据预处理：将三维矩阵X按照采集变量方向展开，即保持三维矩阵X变量维度不变，依次将X_i(K_i×J)从上至下排列得到二维矩阵

设二维矩阵X内任意一点的变量为x_i,j，对该变量进行减均值、除以标准差的标准化处理，其中，下标i代表批次，j代表变量，标准化处理的计算公式如下：

$x_{i,j}=\frac{x_{i,j}-{\stackrel{\‾}{x}}_j}{s_j}---\left(1\right)$

其中：是二维矩阵X任一列的均值，s_j是二维矩阵X相应列的标准差，

$\begin{array}{c}{\stackrel{\‾}{x}}_j=\frac{1}{\mathrm{Sum}}\underset{i=1}{\overset{\mathrm{Sum}}{\mathrm{\Σ}}}{x}_{i,j}\\ s_j=\sqrt{\underset{i=1}{\overset{\mathrm{Sum}}{\mathrm{\Σ}}}{(x_{i,j}-{\stackrel{\‾}{x}}_j)}^2/(\mathrm{Sum}-1)}\\ \mathrm{Sum}=\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{K}_i\end{array}---\left(2\right)$

将已标准化的重新排列分离出时间片矩阵X_k(I_k×J)，即如图2所示，其中，k表示采样时间指标，I_k表示每个时间片所包含的批次操作次数；

(3)建立初始时段模型，该步骤包括以下子步骤：

(3.1)对步骤2标准化处理后得到的时间片矩阵X₁(I₁×J)表示为X_v,1(I₁×J)进行PCA分解，建立初始时段模型，其中PCA分解公式如下：

$X_{v,1}=T^{\′}{\left(P_{v,1}^{\′}\right)}^T=\underset{r=1}{\overset{J}{\mathrm{\Σ}}}{t}_r{p}_r;---\left(3\right)$

其中：下标v表示时段模型按变量方式展开，t_r为正交的主元向量，p_r为正交归一化的负载向量，r表示不同的PCA分解方向，上标T表示矩阵的转置；T′(I₁×J)代表保留全部主元的得分矩阵，P′_v,1(J×J)代表对应的负载矩阵；

(3.2)选取主元个数，将公式(3)重新表述成如下形式：

$X_{v,1}={{\mathrm{TP}}_{v,1}}^T+E=\underset{r=1}{\overset{R}{\mathrm{\Σ}}}{t}_r{p}_r+E---\left(4\right)$

其中：r表示不同的PCA分解方向，R为主元个数；T与P_v,1(J×R)分别保留R个主元后的得分矩阵和负载矩阵，E为残差矩阵；

(3.3)计算SPE指标，确定初始时段模型控制限：

SPE_i＝e_i ^Te_i   (5)

其中，下标i表示时间片中的批次，e_i是第i批次的残差列向量，即e_i为E的第i行向量的转置向量；确定出初始时段模型的控制限Ctr_v,1＝α*max(SPE_i)，其中，α是缓和因子；

(4)衡量变量相关性：调用初始时段模型依次计算下一个时间片每个批次的SPE指标：

t_i＝P_v,1x_i

$\begin{array}{c}{e}_i=x_i-P_{v,1}^T{t}_i\\ {\mathrm{SPE}}_i=e_i^T{e}_i\end{array}---\left(6\right)$

其中，下标i表示时间片中的批次，x_i为时间片中的每一个批次，e_i是第i批次的残差列向量，即e_i为E的第i行向量的转置向量，t_i表示第i批次的主元得分向量；如果发现SPE_i>Ctrl_v,1，则继续调用初始时段模型计算该批次下两个时刻的SPE指标，如果有连续三个样本呈现SPE_i>Ctrl_v,1，那么该批次变量相关性与初始时段模型所表征的特性差别显著，记录该批次第一次超限时间为k^*，则该批次k^*时刻之前的数据可认为是一个子时段；如果没有连续三个样本呈现SPE_i>Ctrl_v,1，则认为该批次仍处于当前时段没有进入下一时段；

(5)时段模型更新：将下一时间片中仍处于当前时段的批次加入初始时段模型的建模数据中并按变量展开，得到对重新按步骤3进行PCA建模，得到更新后的模型以及控制限Ctrl_v,2；

(6)确定所有批次的第一时段：重复步骤4～5，直至所有批次的第一时段均被划分出来；

(7)过程数据更新，确定所有时段：移除已划分出的不等长子时段数据，将剩余数据重新对齐作为下一次时段划分的数据；重复步骤3～6直到所有不等长时段均被划分出来，至此，每个不等长时段的最短批次长度k_s和最长批次长度k_l均被确定出来；

(8)基于不等长时段划分结果的统计建模，该步骤包括以下子步骤：

(8.1)将每个时段内k_s+1至k_l时刻的时间片X_k(I_k×J)按变量展开构造泛化时间片则k_s+1至k_l时刻的时间片X_k(I_k×J)均可由泛化时间片表示；(8.2)根据步骤7的时段划分结果，每个时段内的时间片X_k(I_k×J)及泛化时间片按公式(1)进行标准化处理，获得时段内每个时刻的均值和标准差；

(8.3)将每个时段内1至k_s时刻的时间片X_k(I_k×J)和泛化时间片按变量展开得到对进行PCA建模：

T_c＝X_cP_c，s

$\begin{array}{c}{X}_c=T_c{P}_{c,s}^T+E_c=X_c{P}_{c,s}{P}_{c,s}^T+E_c\\ {\hat{X}}_c=T_c{P}_{c,s}^T=X_c{P}_{c,s}{P}_{c,s}^T\\ E_c=X_c-{\hat{X}}_c=X_c{P}_{c,e}{P}_{c,e}^T\end{array}---\left(7\right)$

其中，下标c表示时段指标，s表示主元空间，e表示残差空间；K_c表示时段所包含的采样时刻个数，P_c,s(J×R_c)是负载矩阵，是主元得分矩阵，R_c是所保留的主元个数，是由主元得分重构的时段数据矩阵，E_c是残差矩阵，P_c,e(J×R_r)(R_r＝J-R_c)反映了时段内次要的波动方向；

(9)计算实时监测统计指标并确定其控制限，该步骤包括以下子步骤：

(9.1)按公式(8)计算Hotelling-T²统计指标和SPE统计量；

$\begin{array}{c}{T_k}^2={(t_k-{\stackrel{\‾}{t}}_k)}^T{\mathrm{\Σ}}_c^{-1}(t_k-{\stackrel{\‾}{t}}_k)\\ {\mathrm{SPE}}_k=e_k^T{e}_k\end{array}---\left(8\right)$

其中，t_k(R_c×1)是每个批次每个时刻的主元得分向量，它可以从中分离得到；∑_c是由构成的协方差矩阵；表示所有t_k(R_c×1)的均值向量；e_k(J×1)是每个批次每个时刻的残差向量，它可以从中分离得到；

(9.2)确定控制限：根据过程数据服从多维高斯分布可知SPE统计量服从加权χ²分布，从而确定SPE统计量的控制限；将每个时刻的Hotelling-T²统计指标从小到大排列，选取第0.95I_k或者0.90I_k个Hotelling-T²统计指标作为其控制限；

(10)基于时段模型的注塑成型过程在线过程监测：基于步骤7划分得到的不等长时段、步骤8建立的时段模型和步骤9所得的两个监测统计量可在线监测注塑成型过程的运行状态；该步骤由以下子步骤来实现：

(10.1)采集新测量数据：注塑成型过程在线监测时，采集新的过程测量数据x_new(J×1)，其中，下标new代表新样本，J为测量变量，与步骤1中的测量变量相同；对于不同时刻的新数据，需要采取不同的监测策略；

(10.2)对于位于时段内[1,k_s]时刻的新数据，按如下子步骤进行在线监测：

(10.2.1)根据步骤8中获得的时段内每个时刻的均值和标准差，根据过程时间的指示调用当前时段相应时刻的均值和标准差，对采集到的过程测量数据如公式(1)中所示进行标准化预处理；

(10.2.2)计算新监测统计量：数据预处理后，调用对应该新采样时刻所在时段的模型P_c(J×R_c)，按照公式(9)估计对应的Hotelling-T²与SPE两个监测统计指标：

t_new ^T＝x_new ^TP_c

$\begin{array}{c}{\hat{x}}_{\mathrm{new}}=P_c{t}_{\mathrm{new}}\\ T_{\mathrm{new}}^2={(t_{\mathrm{new}}-{\stackrel{\‾}{t}}_k)}^T{\mathrm{\Σ}}_c^{-1}(t_{\mathrm{new}}-{\stackrel{\‾}{t}}_k)\\ {\mathrm{SPE}}_{\mathrm{new}}={(x_{\mathrm{new}}-{\hat{x}}_{\mathrm{new}})}^T(x_{\mathrm{new}}-{\hat{x}}_{\mathrm{new}})\end{array}---\left(9\right)$

其中，x_new是新的过程测量数据，t_new是新数据的主元得分向量，是重构的过程测量数据，是根据训练数据获得的主元得分的均值向量，∑_c是由构成的协方差矩阵；

(10.2.3)在线判断注塑成型过程运行状态：实时比较两个监测指标与其各自的统计控制限，若两个监测指标都位于统计控制限之内，表明过程运行正常；若有一个以上监测指标超出正常控制限，表明过程有异常状况发生。

(10.3)对于位于时刻[k_s+1,k_l]的新数据，其时段隶属不能根据时刻来判断，则按如下子步骤进行在线监测：

(10.3.1)根据步骤8中获得的时段内每个时刻的均值和标准差，先根据过程时间的指示调用当前时段的相应时刻的均值和标准差，对采集到的过程测量数据如公式(1)中所示进行标准化预处理；

(10.3.2)计算新监测统计量：数据预处理后，调用对应该新采样时刻所在时段的模型P_c(J×R_c)，按照公式(9)计算得到主元得分，估计残差及其对应的Hotelling-T²与SPE两个监测统计指标；

(10.3.3)在线判断注塑成型过程运行状态：实时比较两个监测指标与其各自的统计控制限，若两个监测指标都位于统计控制限之内，表明过程位于当前时段，且运行正常；若有一个以上监测指标超出正常控制限，表明过程有异常状况发生，此时需要判别过程是进入下一时段还是有故障发生；

(10.3.4)根据步骤8中所获得的时段内每个时刻的均值和标准差，调用下一时段相应时刻的均值和标准差对新数据进行标准化预处理；

(10.3.5)计算新监测统计量：数据预处理后，调用下一时段的模型P_c+1(J×R_c+1)，按照公式(9)计算得到主元得分，估计残差及其对应的Hotelling-T²与SPE两个监测统计指标；

(10.3.6)在线判断注塑成型过程运行状态：实时比较两个监测指标与其各自的统计控制限，若两个监测指标都位于统计控制限之内，表明过程从当前时段进入下一时段，且运行正常；若有一个以上监测指标超出正常控制限，表明过程有异常状况发生。

实施例2

步骤与实施例1相同，在步骤1中，测量变量为9个；采集了23个注塑成型正常批次作为历史训练数据，12个正常批次以及7个注塑成型故障批次作为测试数据；每个批次的采样长度K_i从532个采样点到547个采样点不等。步骤3.2中，所保留主元个数R＝1。步骤8.3中，所保留主元个数R_c＝2。

根据基于不等长时段建立的监测模型以及相应的监测策略，工程师可以实时获得新过程的采样数据的在线监测结果，判断过程的运行状态，从而进行相应的调整措施，保证过程的安全生产和产品的高质量性。本实例中，为了更清楚阐述本发明的有效性，整个注塑过程只有注塑段是不等长的。通过调整注塑段的注塑速率，可以获得不同长度的注塑段数据。这里，设置3个不同速率22mm/s、24mm/s和26mm/s，那么，最长的注塑段长度为99个样本点，最短的注塑段为84个样本点。共采集35个正常批次，23个作为训练数据建立监测模型，12个批次作为测试数据。此外，采集了两类故障批次共7个作为测试数据。两种故障分别是注塑段故障和保压压力故障。图3～4显示了本发明所划分出的不等长注塑段与真实注塑段的对比结果(α＝1.2)，其中，横坐标代表批次，纵坐标表示时段的长度；图3的对象为训练批次；图4的对象为测试批次，点线均表示划分结果；圈线均表真实时段。可以看出，本发明所划分注塑段能很好与真实注塑段相吻合，说明了本发明时段划分的有效性。此外，为了更准确描述时段划分结果，定义一个衡量指标Δt＝|phase_real-phase_partition|，其中，phase_real表示真实时段，phase_partition表示本发明所划分出的时段。分别计算Δt的均值(Mean)和平均绝对偏差(MAD)，对于12个正常测试数据，Δt＝0.08±0.15；对于训练数据，Δt＝0.13±0.23。说明本发明有较好不等长时段划分的精度。为了说明本发明中α如何影响不等长时段划分的结果，图5～8展示了一个批次在不同α下的注塑段的划分结果。其中，横坐标表示样本点，纵坐标表示时段序列。可以看出，α越小，注塑段被划分成越多的子时段；随着α增大，所划分的子时段越来越少，当α增大到一定程度，注塑段和下一时段则不能被很好地区分开来。图9展示了传统的基于聚类的不等长时段划分方法的时段划分结果。其所采用的过程数据与本发明实例中完全一致。可以看出，所得时段划分结果杂乱无序，处于不同时间区域的样本点被划入同一时段，时段划分结果需要大量后续处理。此外，不等长的注塑段不能准确地被划分出来。为了验证本发明的监测性能的优越性，表1展示了本发明的在线监测结果与基于聚类的方法的在线监测结果。其中，FAR表示误报率；FAT表示第一次报警时间，FOT表示故障实际发生时间，ΔT＝FAT-FOT衡量了检测故障的能力。对于FAR和ΔT，计算其均值(Mean)和平均绝对偏差(MAD)。通过配对t检验可以看出，本发明的最优结果明显比基于聚类方法的最优结果有更好的精度，说明了本发明比传统的方法具有更好的监测性能。表2显示了本发明所划分的不等长注塑段和之后等长的保压段的正常测试数据的在线监测结果，通过配对t检验可以得知，注塑段与保压段的监测性能没有显著差别，从而说明了本发明有较好的不等长时段划分精度及在线监测性能。表3总结了本发明对于两类故障的在线监测结果，可以看出，对于两个监测统计量，ΔT的均值和平均绝对偏差都较小，说明本发明对故障检测的有效性。

表1.本发明方法与基于聚类方法的在线监测性能对比结果

表2.针对两个时段的正常测试数据的在线监测结果(α＝1.2，衡量指标FAR％)

表3.针对两类故障数据在线监测结果(α＝1.2，衡量指标ΔT)

图10显示了本发明区分时段转换和故障的能力，其中，粗虚线表示注塑段SPE控制限，细虚线表示保压段SPE控制限，叉线表示注塑段SPE统计量，点线表示保压段SPE统计量。可以看出正常数据的统计指标在注塑段很好地被控制限容纳，而进入保压段后，则超出了注塑段控制限而被保压段控制限容纳，说明时段转换被很好的检测出来。图11显示了本发明对故障的检测能力，其中，图中部虚线表示注塑段SPE控制限，图上部虚线表示保压段SPE控制限，黑色点线表示注塑段SPE统计量，灰色点线表示保压段SPE统计量。该故障发生在保压段，可以看出，当调用保压段监测模型后，新数据统计指标连续超出控制限，说明此时不是时段转换，而是发生了故障。

应该理解，本发明并不局限于上述具体实施例的注塑过程，凡是熟悉本领域的技术人员在不违背本发明精神的前提下还可做出等同变形或替换，这些等同的变型或替换均包含在本申请权利要求所限定的范围内。

Claims (2)

1.一种基于不等长时段自动有序划分的过程监测方法，其特征在于，该方法包括以下步骤：

(1)获取过程分析数据：设每一个不等长批次操作具有J个测量变量和K_i个采样点，则每一个测量批次可得到一个K_i×J的矩阵X_i(K_i×J)，重复I批次的测量步骤后，得到的数据表述为一个不规则的三维矩阵X，其中，i表示批次，测量变量为批次运行过程中可被测量的状态参数；

(2)数据预处理：将三维矩阵X按照采集变量方向展开，即保持三维矩阵X变量维度不变，依次将X_i(K_i×J)从上至下排列得到二维矩阵

设二维矩阵X内任意一点的变量为x_i,j，对该变量进行减均值、除以标准差的标准化处理，其中，下标i代表批次，j代表变量，标准化处理的计算公式如下：

$x_{i,j}=\frac{x_{i,j}-{\stackrel{\‾}{x}}_j}{s_j}---\left(1\right)$

其中：是二维矩阵X任一列的均值，s_j是二维矩阵X相应列的标准差，

${\stackrel{\‾}{x}}_j=\frac{1}{\mathrm{Sum}}\underset{i=1}{\overset{\mathrm{Sum}}{\mathrm{\Σ}}}{x}_{i,j}$

$s_j=\sqrt{\underset{i=1}{\overset{\mathrm{Sum}}{\mathrm{\Σ}}}{(x_{\mathrm{ij}}-{\stackrel{\‾}{x}}_j)}^2/(\mathrm{Sum}-1)}---\left(2\right)$

$\mathrm{Sum}=\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{K}_i$

将已标准化的重新排列分离出时间片矩阵X_k(I_k×J)，即其中，k表示采样时间指标，I_k表示每个时间片所包含的批次操作次数；

(3)建立初始时段模型，该步骤包括以下子步骤：

(3.1)对步骤2标准化处理后得到的时间片矩阵X₁(I₁×J)表示为X_v,1(I₁×J)进行PCA分解，建立初始时段模型，其中PCA分解公式如下：

$X_{v,1}=T^{\′}{\left(P_{v,1}^{\′}\right)}^T=\underset{r=1}{\overset{J}{\mathrm{\Σ}}}{t}_r{p}_r;---\left(3\right)$

其中：下标v表示时段模型按变量方式展开，t_r为正交的主元向量，p_r为正交归一化的负载向量，r表示不同的PCA分解方向，上标T表示矩阵的转置；T′(I₁×J)代表保留全部主元的得分矩阵，P′_v,1(J×J)代表对应的负载矩阵；

(3.2)选取主元个数，将公式(3)重新表述成如下形式：

$X_{v,1}={{\mathrm{TP}}_{v,1}}^T+E=\underset{r=1}{\overset{R}{\mathrm{\Σ}}}{t}_r{p}_r+E---\left(4\right)$

其中：r表示不同的PCA分解方向，R为主元个数；T与P_v,1(J×R)分别保留R个主元后的得分矩阵和负载矩阵，E为残差矩阵；

(3.3)计算SPE指标，确定初始时段模型控制限：

SPE_i＝e_i ^Te_i    (5)

其中，下标i表示时间片中的批次，e_i是第i批次的残差列向量，即e_i为E的第i行向量的转置向量；确定出初始时段模型的控制限Ctr_v,1＝α*max(SPE_i)，其中，α是缓和因子；

(4)衡量变量相关性：调用初始时段模型依次计算下一个时间片每个批次的SPE指标：

t_i＝P_v,1x_i

$e_i=x_i-P_{v,1}^T{t}_i---\left(6\right)$

${\mathrm{SPE}}_i=e_i^T{e}_i$

其中，下标i表示时间片中的批次，x_i为时间片中的每一个批次，e_i是第i批次的残差列向量，即e_i为E的第i行向量的转置向量，t_i表示第i批次的主元得分向量；如果发现SPE_i>Ctrl_v,1，则继续调用初始时段模型计算该批次下两个时刻的SPE指标，如果有连续三个样本呈现SPE_i>Ctrl_v,1，那么该批次变量相关性与初始时段模型所表征的特性差别显著，

记录该批次第一次超限时间为k^*，则该批次k^*时刻之前的数据可认为是一个子时段；如果没有连续三个样本呈现SPE_i>Ctrl_v,1，则认为该批次仍处于当前时段没有进入下一时段；

(5)时段模型更新：将下一时间片中仍处于当前时段的批次加入初始时段模型的建模数据中并按变量展开，得到对重新按步骤3进行PCA建模，得到更新后的模型以及控制限Ctrl_v,2；

(6)确定所有批次的第一时段：重复步骤4～5，直至所有批次的第一时段均被划分出来；

(7)过程数据更新，确定所有时段：移除已划分出的不等长子时段数据，将剩余数据重新对齐作为下一次时段划分的数据；重复步骤3～6直到所有不等长时段均被划分出来，至此，每个不等长时段的最短批次长度k_s和最长批次长度k_l均被确定出来；

(8)基于不等长时段划分结果的统计建模，该步骤包括以下子步骤：

(8.1)将每个时段内k_s+1至k_l时刻的时间片X_k(I_k×J)按变量展开构造泛化时间片则k_s+1至k_l时刻的时间片X_k(I_k×J)均可由泛化时间片表示；

(8.2)根据步骤7的时段划分结果，每个时段内的时间片X_k(I_k×J)及泛化时间片按公式(1)进行标准化处理，获得时段内每个时刻的均值和标准差；

(8.3)将每个时段内1至k_s时刻的时间片X_k(I_k×J)和泛化时间片按变量展开得到对进行PCA建模：

T_c＝X_cP_c，s

$\begin{array}{c}{X}_c=T_c{P}_{c,s}^T+E_c=X_c{P}_{c,s}{P}_{c,s}^T+E_c\\ {\hat{X}}_c=T_c{P}_{c,s}^T=X_c{P}_{c,s}{P}_{c,s}^T\end{array}---\left(7\right)$

$E_c=X_c-{\hat{X}}_c=X_c{P}_{c,e}{P}_{c,e}^T$

其中，下标c表示时段指标，s表示主元空间，e表示残差空间；Kc表示时段所包含的采样时刻个数，P_c,s(J×R_c)是负载矩阵，是主元得分矩阵，R_c是所保留的主元个数，是由主元得分重构的时段数据矩阵，E_c是残差矩阵，P_c,e(J×R_r)(R_r＝J-R_c)反映了时段内次要的波动方向；

(9)计算实时监测统计指标并确定其控制限，该步骤包括以下子步骤：

(9.1)按公式(8)计算Hotelling-T²统计指标和SPE统计量；

$\begin{array}{c}{T_k}^2={(t_k-{\stackrel{\‾}{t}}_k)}^T{\mathrm{\Σ}}_c^{-1}(t_k-{\stackrel{\‾}{t}}_k)\\ {\mathrm{SPE}}_k=e_k^T{e}_k\end{array}---\left(8\right)$

其中，t_k(R_c×1)是每个批次每个时刻的主元得分向量，它可以从中分离得到；∑_c是由构成的协方差矩阵；表示所有t_k(R_c×1)的均值向量；e_k(J×1)是每个批次每个时刻的残差向量，它可以从中分离得到；

(9.2)确定控制限：根据过程数据服从多维高斯分布可知SPE统计量服从加权χ²分布，从而确定SPE统计量的控制限；将每个时刻的Hotelling-T²统计指标从小到大排列，选取第0.95I_k或者0.90I_k个Hotelling-T²统计指标作为其控制限；

(10)基于时段模型的在线过程监测：基于步骤7划分得到的不等长时段、步骤8建立的时段模型和步骤9所得的两个监测统计量可在线监测间歇过程的运行状态；该步骤由以下子步骤来实现：

(10.1)采集新测量数据：在线监测时，采集新的过程测量数据x_new(J×1)，其中，下标new代表新样本，J为测量变量，与步骤1中的测量变量相同；对于不同时刻的新数据，需要采取不同的监测策略；

(10.2)对于位于时段内[1,k_s]时刻的新数据，按如下子步骤进行在线监测：

(10.2.1)根据步骤8中获得的时段内每个时刻的均值和标准差，根据过程时间的指示调用当前时段相应时刻的均值和标准差，对采集到的过程测量数据如公式(1)中所示进行标准化预处理；

(10.2.2)计算新监测统计量：数据预处理后，调用对应该新采样时刻所在时段的模型P_c(J×R_c)，按照公式(9)估计对应的Hotelling-T²与SPE两个监测统计指标：

t_new ^T＝x_new ^TP_c

$\begin{array}{c}{\hat{x}}_{\mathrm{new}}=P_c{t}_{\mathrm{new}}\\ T_{\mathrm{new}}^2={(t_{\mathrm{new}}-{\stackrel{\‾}{t}}_k)}^T{\mathrm{\Σ}}_c^{-1}(t_{\mathrm{new}}-{\stackrel{\‾}{t}}_k)\end{array}---\left(9\right)$

${\mathrm{SPE}}_{\mathrm{new}}={(x_{\mathrm{new}}-{\hat{x}}_{\mathrm{new}})}^T(x_{\mathrm{new}}-{\hat{x}}_{\mathrm{new}})$

其中，x_new是新的过程测量数据，t_new是新数据的主元得分向量，是重构的过程测量数据，是根据训练数据获得的主元得分的均值向量，∑_c是由构成的协方差矩阵；

(10.2.3)在线判断过程运行状态：实时比较两个监测指标与其各自的统计控制限，若两个监测指标都位于统计控制限之内，表明过程运行正常；若有一个以上监测指标超出正常控制限，表明过程有异常状况发生。

(10.3)对于位于时刻[k_s+1,k_l]的新数据，其时段隶属不能根据时刻来判断，则按如下子步骤进行在线监测：

(10.3.1)根据步骤8中获得的时段内每个时刻的均值和标准差，先根据过程时间的指示调用当前时段的相应时刻的均值和标准差，对采集到的过程测量数据如公式(1)中所示进行标准化预处理；

(10.3.2)计算新监测统计量：数据预处理后，调用对应该新采样时刻所在时段的模型P_c(J×R_c)，按照公式(9)计算得到主元得分，估计残差及其对应的Hotelling-T²与SPE两个监测统计指标；

(10.3.3)在线判断过程运行状态：实时比较两个监测指标与其各自的统计控制限，若两个监测指标都位于统计控制限之内，表明过程位于当前时段，且运行正常；若有一个以上监测指标超出正常控制限，表明过程有异常状况发生，此时需要判别过程是进入下一时段还是有故障发生；

(10.3.4)根据步骤8中所获得的时段内每个时刻的均值和标准差，调用下一时段相应时刻的均值和标准差对新数据进行标准化预处理；

(10.3.5)计算新监测统计量：数据预处理后，调用下一时段的模型P_c+1(J×R_c+1)，按照公式(9)计算得到主元得分，估计残差及其对应的Hotelling-T²与SPE两个监测统计指标；

(10.3.6)在线判断过程运行状态：实时比较两个监测指标与其各自的统计控制限，若两个监测指标都位于统计控制限之内，表明过程从当前时段进入下一时段，且运行正常；若有一个以上监测指标超出正常控制限，表明过程有异常状况发生。

2.根据权利要求1所述一种基于不等长时段自动有序划分的过程监测方法，其特征在于，所述步骤1中，所述测量变量为以下9个：喷嘴温度、喷嘴压力、螺杆行程、螺杆速度、注射压力、塑化压力、压力阀开度、流量阀开度、模腔压力。