CN103336507A - 基于多模态协同时段自动划分的统计建模与在线监测方法 - Google Patents

Abstract

本发明公开了一种基于多模态协同时段自动划分的建模与在线监测方法，该方法兼顾批次轴与时间方向上的过程特性变化，同时考虑时段运行的时序性，对所有模态进行协同时段划分，在不同模态间得到统一的时段划分结果，并对每个模态内相似的过程特性建立统一的时段模型，简化了建模的复杂度。基于协同时段划分的结果，分析了各模态间的相对变化，针对不同波动类型建立多模态统计模型应用于在线监测，提高了在线监测性能。该方法易于实施，在注塑成型过程中得到成功应用，不仅有利于对具体过程特性的了解，而且增强了实际在线过程监测的可靠性和可信度，有助于判断工业过程运行状态，及时发现故障，从而保证实际生产的安全可靠运行和产品的高质量追求。

Description

基于多模态协同时段自动划分的统计建模与在线监测方法

技术领域

本发明属于间歇过程统计监测领域，特别是涉及一种基于多模态协同时段自动划分的统计建模与在线监测方法。

背景技术

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

传统的多向主元分析（MPCA）和多向最小偏二乘回归（MPLS）方法将一次间歇操作的所有数据当做一个整体来对待，不能反映变量的相关性在时间方向上的变化。此外，在线应用时必须要预估未知测量数据，因此很难实现在线实施，阻碍了其在实际生产中的广泛应用。间歇过程具有多时段特性。在同一批次中，歇操作中的过程变量相关关系并非随时间时刻变化，而是跟随过程操作进程或过程机理特性的变化发生规律性的改变，呈现分段性。在不同时段中，每个时段具有不同的过程变量轨迹、运行模式以及相关性特征，变量相关性有显著差异。在同一时段中，不同采样时刻过程变量的相关关系却近似一致。考虑到间歇过程的多时段特性，基于时段的建模方法将整个间歇过程划分为不同时段，从而可以建立基于时段的多个模型，并用于在线过程监测。

但是上述方法针对的是单一模态，即间歇过程处于相同的运行条件或单一设备中，而在控制技术及信息技术飞速发展、实现不断飞跃的今天，单一设备、单一模态的“单打独斗型”生产模式已无法满足企业和市场要求，而更倾向于多机、多生产品种、多运行模态的“集团军”协同作战。

以注塑成型过程为例，对于单一模态，注塑过程是个典型的多时段间歇生产过程。一个完整的注塑过程由闭模、注射座前进、注射、保压、塑化、冷却、开模、制件顶出等程序组成，而注射段，保压段和冷却段是决定制品品质的最终要阶段。在这三个时段中起关键作用的变量各不相同，其各自的过程特征也不相同，故需要分时段进行分析建模。而由于市场的需求，我们需要生产不同的注塑产品。对于不同的注塑产品，我们所关注的质量有所不同，如尺寸、外观、重量、耐温性、耐蚀性等等。因此，注塑过程生产的工艺参数、生产条件、生产设备也不尽相同，导致注塑成型过程具有多模态特性。如果将不同模态下的注塑过程考虑成一个整体，则不能反映模态间变量相关性的差异。考虑到在不同模态中，变量相关关系具有不同的特性，我们期望在进行时段划分时，同时对所有模态进行分析（即协同时段划分），这样可以得到一个统一的时段划分结果，从而可以基于时段对不同模态间关系进行有效分析，获取更多过程信息进行更高效的在线过程监测。

综合国内外的研究现状，现有关于批次过程监控的研究仍主要针对单机单工况进行。而大工厂环境下的集成化协同生产模式是批次过程的新发展趋势，多机多工况是集成化批次过程的显著特征。集成化批次生产条件下，针对不同设备、不同工况的海量数据，如何进行综合分析、处理，从中获取价值并加以合理有效利用；此外，不同工况间信息如何相互参考、进行有效共享，多工况条件下对过程多时段特性如何统筹考虑与分析，进而监测过程运行、指导批次生产。这些问题的思考和解决对于保证批次生产的稳定、可靠和高效运行有着极为重要的意义。遗憾的是，目前对于集成化批次过程多机多工况的研究还鲜有涉及。因此，开展面向集成化批次生产的多工况过程监测研究是十分必要的，而且恰逢其时。

本发明的内容以注塑成型这一典型间歇过程为实例，深入分析了间歇过程不同模态的过程时变特性，并能够兼顾各模态时段运行的时序性，基于此定量化模态间相对变化大小并确定相对变化的类型，从而针对不同波动类型建立不同的统计模型并应用于在线监测。到目前为止，尚未见与本发明相关报道。

发明内容

本发明的目的在于针对间歇过程多模态和多时段的特性下现有技术的不足，提供一种基于协同时段划分的多模态建模及在线监测方法。

本发明的目的是通过以下技术方案来实现的：

本发明的有益效果是，该方法以注塑成型过程为实例，分析了多模态间的相关性，对所有模态同时进行时段划分。并基于该时段划分结果，分析了多模态间的相对变化并建立不同的监测模型，提高了模型精度及在线监测故障诊断的效率，并加深了对多模态间歇过程运行状态的理解。最终可应用于多机多模态的集成化批次生产过程，提高生产效率，确保间歇生产的安全可靠运行以及产品的高质量追求。

附图说明

图1是多模态协同时段划分及统计建模流程图；

图2是基于多模态协同时段划分的在线监测流程图；

图3是多模态协同时段划分结果图；

图4是模态1独立时段划分结果图；

图5是模态2独立时段划分结果图；

图6是模态3独立时段划分结果图；

图7是相同模态监测图（建模对象：参考模态；监测对象：参考模态测试数据;监测指标：

）；

图8是相同模态监测图（建模对象：参考模态；监测对象：参考模态测试数据;监测指标：SPE_r）；

图9是不同模态监测图（建模对象：备选模态1；监测对象：参考模态测试数据;监测指标：

）；

图10是不同模态监测图（建模对象：备选模态1；监测对象：参考模态测试数据;监测指标：T_n ²）；

图11是不同模态监测图（建模对象：备选模态1；监测对象：参考模态测试数据;监测指标：T_e ²）；

图12是不同模态监测图（建模对象：备选模态1；监测对象：参考模态测试数据;监测指标：SPE_f）；

图13是不同模态监测图（建模对象：备选模态2；监测对象：参考模态测试数据;监测指标：）；

图14是不同模态监测图（建模对象：备选模态2；监测对象：参考模态测试数据;监测指标：T_n ²）；

图15是不同模态监测图（建模对象：备选模态2；监测对象：参考模态测试数据;监测指标：T_e ²）；

图16是不同模态监测图（建模对象：备选模态2；监测对象：参考模态测试数据;监测指标：SPE_f）；

图17是不同模态监测图（建模对象：参考模态；监测对象：备选模态1测试数据;监测指标：

）；

图18是不同模态监测图（建模对象：参考模态；监测对象：备选模态1测试数据;监测指标：SPE_r）；

图19是相同同模态监测图（建模对象：备选模态1；监测对象：备选模态1测试数据;监测指标：

）；

图20是相同模态监测图（建模对象：备选模态1；监测对象：备选模态1测试数据;监测指标：T_n ²）；

图21是相同同模态监测图（建模对象：备选模态1；监测对象：备选模态1测试数据;监测指标：T_e ²）；

图22是相同模态监测图（建模对象：备选模态1；监测对象：备选模态1测试数据;监测指标：SPE_f）；

图23是不同模态监测图（建模对象：备选模态2；监测对象：备选模态1测试数据;监测指标：

）；

图24是不同模态监测图（建模对象：备选模态2；监测对象：备选模态1测试数据;监测指标：T_n ²）；

图25是不同模态监测图（建模对象：备选模态2；监测对象：备选模态1测试数据;监测指标：T_e ²）；

图26是不同模态监测图（建模对象：备选模态2；监测对象：备选模态1测试数据;监测指标：SPE_f）；

图27是不同模态监测图（建模对象：参考模态；监测对象：备选模态2测试数据;监测指标：

）；

图28是不同模态监测图（建模对象：参考模态；监测对象：备选模态2测试数据;监测指标：SPE_r）；

图29是不同模态监测图（建模对象：备选模态1；监测对象：备选模态2测试数据;监测指标：

）；

图30是不同模态监测图（建模对象：备选模态1；监测对象：备选模态2测试数据;监测指标：T_n ²）；

图31是不同模态监测图（建模对象：备选模态1；监测对象：备选模态2测试数据;监测指标：T_e ²）；

图32是不同模态监测图（建模对象：备选模态1；监测对象：备选模态2测试数据;监测指标：SPE_f）；

图33是相同模态监测图（建模对象：备选模态2；监测对象：备选模态2测试数据;监测指标：

）；

图34是相同模态监测图（建模对象：备选模态2；监测对象：备选模态2测试数据;监测指标：T_n ²）；

图35是相同模态监测图（建模对象：备选模态2；监测对象：备选模态2测试数据;监测指标：T_e ²）；

图36是相同模态监测图（建模对象：备选模态2；监测对象：备选模态2测试数据;监测指标：SPE_f）；

具体实施方式

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

注塑成型过程是典型的多时段间歇生产过程，一般由注射、保压、冷却三个阶段构成。此外，塑化过程在冷却初始阶段完成。具体地说，在注射阶段，液压系统推动螺杆将塑料粘流体注入模腔中，直至模腔被流体充满。过程处于保压阶段时，仍有少量的粘流体被高压挤进模腔中，以补偿塑料粘流体在冷却和塑化时造成的体积收缩。保压阶段一直持续到模腔的浇口冻结，过程进入冷却段。当螺杆头部熔料逐渐增多，达到一定的注射量后，螺杆停止后退和转动，这段时间的过程状态也称塑化段。随着模腔中熔料继续冷却，塑件完全固化，模具打开，塑件被顶出，从而完成一个工作循环。此外，由于选用的生产设备不同、运行在不同生产操作条件下、采用不同的过程生产方案、生产不同的产品类型等等，注塑成型过程是一个多模态的的生产过程，不同模态下变量相关关系具有不同的特性。因此，不仅要考虑同一批次内的多时段性，还要考虑不同批次的多模态特性。

本发明基于多模态协同时段自动划分的建模与在线监测方法，包括以下步骤：

步骤1：获取过程分析数据：

设一个多模态的间歇操作过程具有J个测量变量和K个采样点，则每一个测量批次可得到一个K×J的矩阵X(K×J)。故对于每一个模态重复I批次的测量步骤后，得到的数据可以表述为一个三维矩阵

其中,m=1,2,…M,即一共M个模态。

本实例中，选取3组不同的机桶温度和保压压力，而保持其他工作条件不变，从而得到了3个工作模态，对于每个模态，采集了529个样本，11个过程变量（均通过传感器实时测量而得）以及28个批次用于建模，组成三维矩阵

其中，前20个批次用于建模，后8个批次用以测试。其中，11个过程变量分别为：压力阀开度、流量阀开度、螺杆行程、螺杆速度、注射压力、顶出行程、开合模行程、开合模速度、机桶头部温度、机桶中间温度和机桶尾部温度。

步骤2：数据预处理：

将每个三维矩阵

按照采集批次方向展开，即将一个操作批次内的各采样点上的变量按照时间顺序排开得到二维矩阵X_m(I×KJ)，由K个时间片矩阵X_m,k(I×J)组成，其中，下标k为时间指标；

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

$x_{k,i,j}=\frac{x_{k,i.j}-{\stackrel{\‾}{x}}_{k,j}}{S_{k,j}};---\left(1\right)$

其中：k是时间片指标。

是X_k矩阵任一列的均值，s_k,j是X_k矩相应列的标准差，

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

$s_{k,j}=\sqrt{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{(x_{k,i,j}-{\stackrel{\‾}{x}}_{k,j})}^2/(I-1)};---\left(2\right)$

步骤3：时间片PCA建模，该步骤由以下子步骤来实现：

（3.1）对步骤2标准化处理后的每一个时间片矩阵X_m,k(I×J)执行PCA分解，建立时间片PCA模型，其中PCA分解公式如下：

$X_{m,k}=T_{m,k}{P_{m,k}}^T=\underset{r=1}{\overset{J}{\mathrm{\Σ}}}{t}_{m,k,r}{p}_{m,k,r};---\left(3\right)$

其中：t_m,k,r为正交的主元向量，p_m,k,r为正交归一化的负载向量，r表示不同的PCA分解方向，上标T表示矩阵的转置；T_m,k(I×J)代表保留全部主元的得分矩阵，P_m,k(J×J)代表对应的负载矩阵。

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

$X_{m,k}=T_{m,k}{P_{m,k}}^T=\underset{r=1}{\overset{R_{m,k}}{\mathrm{\Σ}}}{t}_{m,k,r}{p}_{m,k,r}+E_{m,k}---\left(4\right)$

其中：r表示不同的PCA分解方向；T_k(I×R_m,k)与P_k(J×R_m,k)代表负载矩阵分别为保留R_m,k个主元后的得分矩阵和负载矩阵，E_m,k为残差矩阵。通过上述变换，多向主元分析法模型将原始数据空间分解为主元空间和残差空间，主元空间内代表主要的系统过程波动信息；这里所保留的主元个数R_m,k能够反映原过程中90％的过程波动信息。

（3.3）计算残差空间中各时间片k中对应各个批次的SPE指标：

SPE_m,k,i=e_m,k,i ^Te_m,k,i   （5）

其中，下标i表示时间片中的批次，e_m,k,i是对应每个模块k时刻第i批次的残差列向量。根据相同时刻上不同批次的SPE值服从带权重系数的χ²分布，从而确定出每个时间点上的控制限Ctr_m,k，它反应了时间片PCA模型的重构能力。

步骤4：确定基于变量展开模型的SPE指标控制限：从间歇过程初始点开始，依次将下一个时间片与之前的时间片组合在一起并按变量方式展开得到时间块

其中上标v代表变量展开方式。对新时间块矩阵进行PCA分析，提取出负载矩阵

计算其SPE值并根据相同时刻上不同批次的SPE值服从带权重系数的χ²分布，从而确定出每个时间点上的控制限Ctr_m,k。

步骤5：确定第一时段划分点k^*：比较在相同时间区域内的每个时间点上Ctr_m,k和

的大小，对于任意一个模态，如果发现连续三个样本呈现那么新加入的时间片对该时间块的PCA监测模型及相应的监测性能都有重大的影响。记加入新时间片前的时刻为k^*。其中，

是依附于Ctr_m,k的常数，称作缓和因子，它反映的是与时间片模型相比，时间块模型允许监测精度损失的程度。则本实例中对于所有3个模态，k^*时刻之前的时间片可认为是一个子时段。

步骤6：过程分析数据更新，确定所有划分时段：对于所有3个模态，根据步骤5中所获得的时刻k^*的指示，移除第一个子时段，把余下的间歇过程数据作为新的输入数据带入到第5步中并重复上述步骤5-6，划分不同时间段，直到没有数据余留。

本实例中，同时考虑了注塑过程3个不同模态的潜在过程特性并进行了协同时段划分。图3显示了3个模态的协同时段划分结果，取值为3。作为对比，图4至图6显示了对每个模态独立进行时段划分的结果。可以看出，对每个模态独立划分所得的结果则呈现显著的模态间的差异；而用本发明所述方法进行时段划分，联合考虑了3个模态的时变特性，故3个模态同时被划分为边界相同的子时段，不需要进行后处理，便于统计建模。

步骤7：基于时段划分结果的多模态统计建模：本实例中，选取第1个模态作为参考模态，其它2个作为备选模态。根据步骤6时段划分结果，对参考模态和备选模态每个时段内的时间片按变量方式展开组合成子时段代表性建模数据X^v _c,r(I_rK_c×J)和

其中，下标r和a表示参考模态和备选模态，c表示时段，上标v标是按变量方式展开，K_c表示时段持续的时间，I_r和I_a分别是模态的批次数。

（7.1）用步骤1所述方法分别对X^v _c,r(I_rK_c×J)和

进行数据标准化处理。

（7.2）对数据标准化处理后的X^v _c,r(I_rK_c×J)进行如下PCA分解:

$\begin{array}{c}{T}_{c,r}^v=X_{c,r}^v{P}_{c,r}\\ E_{c,r}^v=X_{c,r}^v{P}_{c,r}^e{P_{c,r}^e}^T\end{array},---\left(6\right)$

其中，和P_c,r(J×R_c,r)是主元空间中的主元和其相应的负载；和

是残差空间的残差和相应的负载，上标e表示残差空间。R_c,r是由累积的波动解释率所决定的提取的主元个数，

（7.3）将

分别投影到P_c,r和

方向，得到主元空间的主元得分

和残差空间的波动

$\begin{array}{c}{T}_{c,a}^v=X_{c,a}^v{P}_{c,r}\\ E_{c,a}^v=X_{c,a}^v{P}_{c,r}^e{P_{c,r}^e}^T=\underset{j=1}{\overset{R_{c,r}^e}{\mathrm{\Σ}}}{X}_{c,a}^v{p}_{c,a,j}^e{p_{c,a,j}^e}^T\end{array},---\left(7\right)$

在主元空间中定义一个R_c,r维向量Ratio_c,a，表示备选模态和参考模态波动的比值。其中，Ratio_c,a中的每一个元素

var()表示主元围绕参考模态中心的波动值，(:,i)表示矩阵的第i列。

同理，在残差空间中定义一个

维的向量Δ_c,a，表示备选模态和参考模态波动的差值。其中， ${\mathrm{\Δ}}_{c,a,i}=\left|\right|X_{c,a}^v{p}_{c,a,i}^e{p_{c,a,i}^e}^T|{|}^2-|\left|X_{c,r}^v{p}_{c,a,i}^e{p_{c,a,i}^e}^T\right|{|}^2,$ 表示欧式距离。

（7.4）分别选取Ratio_c,a指标和Δ_c,a指标。具体步骤如下：

选取出Ratio_c,a中值大于1的元素所对应的的p_c,r,i，组成

P_c,r中剩余的p_c,r,i组成

其中，

是所选取的主元的个数，

是剩余主元的个数， $R_{c,a}^n=R_{c,r}-R_{c,a}^{*}.$

同理，选取Δ_c,a中大于零的元素所对应的的

构成

中剩余的

构成

其中，

是所选取的主元的个数，

是剩余主元的个数， ${R_{c,a}^e}^n=R_{c,r}^e-{R_{c,a}^e}^{*}.$

（7.5）在主元空间和残差空间分别进行PCA重构，具体步骤如下：

在主元空间中，分别沿

和

方向对

进行重构：

其中，

和

分别表示对报警的T²监测量有贡献和无贡献的波动部分。在残差空间中，沿方向对

进行重构：

显而易见，SPE监测量报警是由导致的。

（7.6）再次PCA分解：

通过上述方法，在备选模态中，有增长的相关波动和无增长的相关波动被分离开来，并分别进行了PCA建模。

步骤8：计算实时监测统计指标：

（8.1）基于步骤7中多模态子空间分解建模的结果，时段c的备选模态与参考模态监测指标通过以下步骤计算出来：

参考模态：

$\begin{array}{c}{T}_{c,r}^v=X_{c,a}^v{P}_{c,r}\\ E_{c,r}^v=X_{c,a}^v{P}_{c,r}^e{P_{c,r}^e}^T\end{array},---\left(10\right)$

其中，

是表示系统波动的主元得分，

是残差。对于备选模态：

$\begin{array}{c}{T}_{c,a,f}=X_{c,a}^v{P}_{c,a,f}\\ T_{c,a,o}=X_{c,a}^v{P}_{c,a,o}\\ T_{c,a,f}^e=X_{c,a}^v{P}_{c,a,f}^e\\ E_{c,a}^f=X_{c,a}^v{P}_{c,r}^e{P_{c,r}^e}^T-X_{c,a}^v{P}_{c,a,f}^e{P_{c,a,f}^e}^T\end{array},---\left(11\right)$

（8.2）计算出各子空间时间片监测指标后，则可得到参考模态与备选模态每个时间点k的监测指标：

参考模态：

$\begin{array}{c}{T}_{r,k,i}^2={(t_{r,k,i}-{\stackrel{\‾}{t}}_{r,k})}^T{\mathrm{\Σ}}_{r,k}^{-1}(t_{r,k,i}-{\stackrel{\‾}{t}}_{r,k})\\ \mathrm{SP}{E}_{r,k,i}=e_{r,k,i}^T{e}_{r,k,i}\end{array}---\left(12\right)$

其中，下标r表示参考模态，下标i表示每个时间片的第i个批次，

表示从T_r,k计算得到的平均向量，∑_r,k是对角阵，每个元素表示时间片得分T_r,k的主元波动。

备选模态：

$\begin{array}{c}{T}_{a,f,k,i}^2={(t_{a,f,k,i}-{\stackrel{\‾}{t}}_{a,f,k})}^T{\mathrm{\Σ}}_{a,f,k}^{-1}(t_{a,f,k,i}-{\stackrel{\‾}{t}}_{a,f,k})\\ T_{a,o,k,i}^2={(t_{a,o,k,i}-{\stackrel{\‾}{t}}_{a,o,k})}^T{\mathrm{\Σ}}_{a,o,k}^{-1}(t_{a,o,k,i}-{\stackrel{\‾}{t}}_{a,o,k})\\ {T_{a,f,k,i}^e}^2={(t_{a,f,k,i}^e-{\stackrel{\‾}{t}}_{a,f,k}^e)}^T{{\mathrm{\Σ}}_{a,f,k}^e}^{-1}(t_{a,f,k,i}^e-{\stackrel{\‾}{t}}_{a,f,k}^e)\\ \mathrm{SP}{E}_{a,k}^f={e_{a,k,i}^f}^T{e}_{a,k,i}^f\end{array}---\left(13\right)$

其中，

分别表示T_a,f,k，T_a,o,k，

的平均向量，经过数据预处理后均为零向量。∑_a,f,k，∑_a,o,k，

为对角阵，其对角元素分别为T_a,f,k，T_a,o,k，

的每个主元的波动值。

（8.3）计算控制限：

T²统计指标用来测量各采样时刻过程变量偏离正常模态下平均轨迹的距离，并且服从带权重的F分布，从而计算得到控制限；同理，SPE近似服从加权χ²分布，从而可计算得控制限。

步骤9：基于多模态时段模型的在线过程监测：

基于步骤6划分的时段、步骤7建立的时段模型监测系统和步骤8所得的监测统计量控制限可以在线监测注塑成型等新运行间歇过程的状态。该步骤由以下子步骤来实现：

（9.1）采集新测量数据及新测量数据预处理

在线监测时，采集到k时刻的新过程测量数据x_new(J×1)（其中下标new代表新样本）后，首先需要进行数据预处理。根据参考模态由步骤2所获得的均值和标准差，根据过程时间的指示调用对应该时刻的均值和标准差对现有数据如公式（1）中所示进行标准化预处理。

本实例中，所采集的新测量数据为注塑过程中的压力阀开度、流量阀开度、螺杆行程、螺杆速度、注射压力、顶出行程、开合模行程、开合模速度、机桶头部温度、机桶中间温度和机桶尾部温度。

（9.2）计算新监测统计量

数据预处理后，根据公式（6）计算的PCA子时段模型，调用对应该新采样时刻所在时段的模型P_c,r(J×R_c,r)，按照如下方式计算得到主元得分，估计残差及其对应的Hotelling-T²与SPE两个监测统计指标：

$\begin{array}{c}{t_{r,\mathrm{new}}}^T={x_{\mathrm{new}}}^T{P}_{c,r}\\ {e_{r,\mathrm{new}}}^T={x_{\mathrm{new}}}^T{P}_{c,r}^e{P_{c,r}^e}^T\\ T_{r,\mathrm{new}}^2={(t_{r,\mathrm{new}}-{\stackrel{\‾}{t}}_{r,k})}^T{\mathrm{\Σ}}_{r,k}^{-1}(t_{r,\mathrm{new}}-{\stackrel{\‾}{t}}_{r,k})\\ \mathrm{SP}{E}_{r,\mathrm{new}}={e_{\mathrm{new}}}^T{e}_{\mathrm{new}}\end{array}---\left(14\right)$

其中,x_new是新的过程测量数据，

是前面根据训练数据获得的T_r,k的均值向量，∑_r,k是X_r,k的协方差矩阵。

（9.3）在线判断过程运行状态

比较两监测指标与各自的控制限，如果

和SPE_r,new监测指标均在控制限内，则过程可认为运行在参考模态正常状态。如果有超限报警发生，为了区分当前过程是处于备选模态正常运行状态还是故障状态，须按如下步骤进行区分：

（9.3.1）计算备选模态的新监测统计量

$\begin{array}{c}{t_{a,f,\mathrm{new}}}^T={x_{\mathrm{new}}}^T{P}_{c,a,f}\\ {t_{a,o,\mathrm{new}}}^T={x_{\mathrm{new}}}^T{P}_{c,a,o}\\ {t_{a,f,\mathrm{new}}^e}^T={x_{\mathrm{new}}}^T{P}_{c,a,f}^e{P_{c,a,f}^e}^T\\ {e_{a,\mathrm{new}}^f}^T={x_{\mathrm{new}}}^T{P}_{c,r}^e{P_{c,r}^e}^T-{x_{\mathrm{new}}}^T{P}_{c,a,f}^e{P_{c,a,f}^e}^T\\ {T_{a,f,\mathrm{new}}}^2={(t_{a,f,\mathrm{new}}-{\stackrel{\‾}{t}}_{a,f,k})}^T{\mathrm{\Σ}}_{a,f,k}^{-1}(t_{a,f,\mathrm{new}}-{\stackrel{\‾}{t}}_{a,f,k})\\ {T_{a,o,\mathrm{new}}}^2={(t_{a,o,\mathrm{new}}-{\stackrel{\‾}{t}}_{a,o,k})}^T{\mathrm{\Σ}}_{a,o,k}^{-1}(t_{a,o,\mathrm{new}}-{\stackrel{\‾}{t}}_{a,o,k})\\ {T_{a,f,\mathrm{new}}^e}^2={(t_{a,f,\mathrm{new}}^e-{\stackrel{\‾}{t}}_{a,f,k}^e)}^T{{\mathrm{\Σ}}_{a,f,k}^e}^{-1}(t_{a,f,\mathrm{new}}^e-{\stackrel{\‾}{t}}_{a,f,k}^e)\\ \mathrm{SP}{E}_{a,\mathrm{new}}^f={e_{a,\mathrm{new}}^f}^T{e}_{a,\mathrm{new}}^f\end{array}---\left(15\right)$

（9.3.2）判断过程运行状态

比较每个备选模态新的监测指标与控制限。如果所有新监测指标都处于某一个备选模态的控制限之内，则称新的样本与该模态匹配，即新过程运行在该备选模态正常状态下；相反，如果没有找到模态与新样本匹配，则新样本发生故障或者处于一种新的运行模态，该模态未包含在建模模型内。

本实例中，在获得了监测模型后，为了验证监测模型的精度，我们采取两种监测手段：

（a）相同模态监测：通过参考模态和备选模态分别建立监测模型，并用于监测与建模数据相同的模态；即建模所用的模态和应用对象所处的模态一致。

（b）不同模态监测：通过参考模态和备选模态分别建立监测模型，并用于监测其他不同模态；即建模所用的模态和应用对象处于不同模态。

对于每种监测手段，我们采用注塑过程正常批次数据来进行。在相同模态监测中，报警越少则模型的精度越高；而在不同模态监测时，报警越多则模型的精度越高。故在相同模态监测时，使用错报率（FAR，Fault Alarming Rate）来评价模型性能；在不同模态监测时，使用漏报率（MAR，Missing Alarming Rate）来评价模型性能。图7、8，图9至图12，图13至图16显示了相同模态监测和不同模态监测的结果。图7、8表示用参考模态进行相同模态监测的结果，显而易见，监测量均处于控制限之内；图9至图16分别表示用备选模态1和2建模对参考模态进行（即不同模态监测）监测的结果，可知大部分监测量都超出控制限。图17至图18、图19至图22和图23至图26是分别用参考模态、备选模态1和2建模，对备选模态1中一组正常测试数据进行监测的结果；图27、28，图29至32和图33至36是分别用参考模态、备选模态1和2建模，对备选模态2中一组正常测试数据进行监测的结果，两图结果都与期望结果相符合。综上，表明模型有良好的精度。

表1、表2和表3、表4分别是采用本发明所述的方法和基于独立模态划分的建模方法，对训练数据和测试数据进行监测的统计结果。可以看出，本发明所述方法与基于独立模态划分的建模方法具有相似的统计精度，表明了本发明所述方法的有效性。

表1采用本发明方法进行的相同模态监测和不同模态监测结果（Mean±MAD）表（监测对象为训练数据，Mean表示平均值，MAD表示平均波动值）

表1采用本发明方法进行的相同模态监测和不同模态监测结果（Mean±MAD）表（监测对象为正常测试数据，Mean表示平均值，MAD表示平均波动值）

其中，带阴影的值表示相同模态监测的结果，评价指标为FAR（%），不加阴影的为不同模态监测的结果，评价指标为MAR（%）。

表3采用单个模态独立建模方法进行的相同模态监测和不同模态监测结果（Mean±MAD）表（监测对象为训练数据，Mean表示平均值，MAD表示平均波动值）

表4采用单个模态独立建模方法进行的相同模态监测和不同模态监测结果（Mean±MAD）（监测对象为正常测试数据，Mean表示平均值，MAD表示平均波动值）

为了验证模型检测故障的能力，我们采取如下手段：

使用相应的模型进行监测，因而报警的原因只能是过程运行故障。这里，设置3种人为故障来阐明本发明所述方法的优越性，即分别对3个系统子空间（P_c,a,f,P_c,a,o and

）施加扰动。为了简化，我们对每个批次施加一个指数曲线的扰动（f(K×1)），指数曲线的形式为K(1-e^-t/τ)，其中时间常数τ=100，增益K=5。为了使每个故障只影响于一个子空间，则扰动的形式为：

x_i ^fT=x_i ^T+Δ_k ^T=x_i ^T+θ_c,1 ^Tf_k（i=1,2,3）

其中，x_i是每个时刻正常状态下的初始测量值构成的向量，Δ_k(J×1)表示施加的扰动，f_k表示指数曲线k时刻的值，θ_c,1表示三个系统子空间中的第一个方向。由于三个系统子空间相互正交，故三种故障互不影响。

表5对两种建模方法（两种方法均基于协同时段划分，一种采用本发明的多模态统计建模方法（BM，Between-mode modelling），另一种采用单个模态独立建模方法（MS，Mode-SPEcific separate modelling））的在线故障检测性能进行了比较，评价指标为相同模态监测中的第一次检测时间（FDT，First DetectionTime），FDT表示第一次连续三次检测到报警的时间。

表5多模态统计建模方法与单模态独立建模方法在相同模态监测中实时故障检测性能结果对比（Mean±MAD）（评价指标为FDT，检测对象为3个人为故障）

显而易见，本发明提出的多模态统计建模方法的FDT值小于单个模态独立建模方法，从而体现了多模态统计建模方法的优越性。

Claims (1)

1.一种基于多模态协同时段自动划分的建模与在线监测方法，其特征在于，该方法包括以下步骤：

步骤1：获取过程分析数据：设一个多模态的间歇操作过程具有J个测量变量和K个采样点，则每一个测量批次可得到一个K×J的矩阵X(K×J)。故对于每一个模态重复I批次的测量步骤后，得到的数据可以表述为一个三维矩阵

其中,m=1,2,…M,即一共M个模态。

步骤2：数据预处理：将每个三维矩阵

按照采集批次方向展开，即将一个操作批次内的各采样点上的变量按照时间顺序排开得到二维矩阵X_m(I×KJ)，由K个时间片矩阵X_m,k(I×J)组成，其中，下标k为时间指标；

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

$x_{k,i,j}=\frac{x_{k,i,j}-{\stackrel{\‾}{x}}_{k,j}}{S_{k,j}};---\left(1\right)$

其中：k是时间片指标。是X_k矩阵任一列的均值，s_k,j是X_k矩相应列的标准差，

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

$s_{k,j}=\sqrt{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{(x_{k,i,j}-{\stackrel{\‾}{x}}_{k,j})}^2/(I-1)};---\left(2\right)$

步骤3：时间片PCA建模，该步骤由以下子步骤来实现：

（3.1）对步骤2标准化处理后的每一个时间片矩阵X_m,k(I×J)执行PCA分解，建立时间片PCA模型，其中PCA分解公式如下：

$X_{m,k}=T_{m,k}{P_{m,k}}^T=\underset{r=1}{\overset{J}{\mathrm{\Σ}}}{t}_{m,k,r}{p}_{m,k,r};---\left(3\right)$

其中：t_m,k,r为正交的主元向量，p_m,k,r为正交归一化的负载向量，r表示不同的PCA分解方向，上标T表示矩阵的转置；T_m,k(I×J)代表保留全部主元的得分矩阵，P_m,k(J×J)代表对应的负载矩阵。

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

$X_{m,k}=T_{m,k}{P_{m,k}}^T=\underset{r=1}{\overset{R_{m,k}}{\mathrm{\Σ}}}{t}_{m,k,r}{p}_{m,k,r}+E_{m,k}---\left(4\right)$

其中：r表示不同的PCA分解方向；T_k(I×R_m,k)与P_k(J×R_m,k)代表负载矩阵分别为保留R_m,k个主元后的得分矩阵和负载矩阵，E_m,k为残差矩阵。通过上述变换，多向主元分析法模型将原始数据空间分解为主元空间和残差空间，主元空间内代表主要的系统过程波动信息；这里所保留的主元个数R_m,k能够反映原过程中90％的过程波动信息。

（3.3）计算残差空间中各时间片k中对应各个批次的SPE指标：

SPE_m,k,i=e_m,k,i ^Te_m,k,i   （5）

其中，下标i表示时间片中的批次，e_m,k,i是对应每个模块k时刻第i批次的残差列向量。根据相同时刻上不同批次的SPE值服从带权重系数的χ²分布，从而确定出每个时间点上的控制限Ctr_m,k，它反应了时间片PCA模型的重构能力。步骤4：确定基于变量展开模型的SPE指标控制限：从间歇过程初始点开始，依次将下一个时间片与之前的时间片组合在一起并按变量方式展开得到时间块

其中上标v代表变量展开方式。对新时间块矩阵进行PCA分析，提取出负载矩阵

计算其SPE值并根据相同时刻上不同批次的SPE值服从带权重系数的χ²分布，从而确定出每个时间点上的控制限Ctr_m,k。

步骤5：确定第一时段划分点k^*：比较在相同时间区域内的每个时间点上Ctr_m,k和的大小，对于任意一个模态，如果发现连续三个样本呈现

那么新加入的时间片对该时间块的PCA监测模型及相应的监测性能都有重大的影响。记加入新时间片前的时刻为k^*。其中，是依附于Ctr_m,k的常数，称作缓和因子，它反映的是与时间片模型相比，时间块模型允许监测精度损失的程度。则对于所有M个模态，k^*时刻之前的时间片可认为是一个子时段。

步骤6：过程分析数据更新，确定所有划分时段：对于所有M个模态，根据步骤5中所获得的时刻k^*的指示，移除第一个子时段，把余下的间歇过程数据作为新的输入数据带入到第5步中并重复上述步骤5-6，划分不同时间段，直到没有数据余留。

步骤7：基于时段划分结果的多模态统计建模：从M个模态中任意选取一个作为参考模态，其它作为备选模态。根据步骤6时段划分结果，对参考模态和备选模态每个时段内的时间片按变量方式展开组合成子时段代表性建模数据X^v _c,r(I_rK_c×J)和

其中，下标r和a表示参考模态和备选模态，c表示时段，上标v标是按变量方式展开，K_c表示时段持续的时间，I_r和I_a分别是模态的批次数。

（7.1）用步骤1所述方法分别对X^v _c,r(I_rK_c×J)和

进行数据标准化处理。

（7.2）对数据标准化处理后的X^v _c,r(I_rK_c×J)进行如下PCA分解:

$\begin{array}{c}{T}_{c,r}^v=X_{c,r}^v{P}_{c,r}\\ E_{c,r}^v=X_{c,r}^v{P}_{c,r}^e{P_{c,r}^e}^T\end{array},---\left(6\right)$

其中，

和P_c,r(J×R_c,r)是主元空间中的主元和其相应的负载；

和

是残差空间的残差和相应的负载，上标e表示残差空间。R_c,r是由累积的波动解释率所决定的提取的主元个数，

（7.3）将

分别投影到P_c,r和

方向，得到主元空间的主元得分

和残差空间的波动

$\begin{array}{c}{T}_{c,a}^v=X_{c,a}^v{P}_{c,r}\\ E_{c,a}^v=X_{c,a}^v{P}_{c,r}^e{P_{c,r}^e}^T=\underset{j=1}{\overset{R_{c,r}^e}{\mathrm{\Σ}}}{X}_{c,a}^v{p}_{c,a,j}^e{p_{c,a,j}^e}^T\end{array},---\left(7\right)$

在主元空间中定义一个R_c,r维向量Ratio_c,a，表示备选模态和参考模态波动的比值。其中，Ratio_c,a中的每一个元素var()表示主元围绕参考模态中心的波动值，(:,i)表示矩阵的第i列。

同理，在残差空间中定义一个

维的向量Δ_c,a，表示备选模态和参考模态波动的差值。其中， ${\mathrm{\Δ}}_{c,a,i}=\left|\right|X_{c,a}^v{p}_{c,a,i}^e{p_{c,a,i}^e}^T|{|}^2-|\left|X_{c,r}^v{p}_{c,a,i}^e{p_{c,a,i}^e}^T\right|{|}^2,$ 表示欧式距离。

（7.4）分别选取Ratio_c,a指标和Δ_c,a指标。具体步骤如下：

选取出Ratio_c,a中值大于1的元素所对应的的p_c,r,i，组成P_c,r中剩余的p_c,r,i组成其中，是所选取的主元的个数，

是剩余主元的个数， $R_{c,a}^n=R_{c,r}-R_{c,a}^{*}.$

同理，选取Δ_c,a中大于零的元素所对应的的

构成

中剩余的

构成

其中，

是所选取的主元的个数，

是剩余主元的个数， ${R_{c,a}^e}^n=R_{c,r}^e-{R_{c,a}^e}^{*}.$

（7.5）在主元空间和残差空间分别进行PCA重构，具体步骤如下：

在主元空间中，分别沿

和方向对进行重构：

其中，和

分别表示对报警的T²监测量有贡献和无贡献的波动部分。

在残差空间中，沿

方向对

进行重构：

显而易见，SPE监测量报警是由

导致的。

（7.6）再次PCA分解：

通过上述方式，在备选模态中，有增长的相关波动和无增长的相关波动被分离开来，并分别进行了PCA建模。

步骤8：计算实时监测统计指标：

（8.1）基于步骤7中多模态子空间分解建模的结果，时段c的备选模态与参考模态监测指标通过以下步骤计算出来：

参考模态：

$\begin{array}{c}{T}_{c,r}^v=X_{c,a}^v{P}_{c,r}\\ E_{c,r}^v=X_{c,a}^v{P}_{c,r}^e{P_{c,r}^e}^T\end{array},---\left(10\right)$

其中，是表示系统波动的主元得分，是残差。

对于备选模态：

$\begin{array}{c}{T}_{c,a,f}=X_{c,a}^v{P}_{c,a,f}\\ T_{c,a,o}=X_{c,a}^v{P}_{c,a,o}\\ T_{c,a,f}^e=X_{c,a}^v{P}_{c,a,f}^e\\ E_{c,a}^f=X_{c,a}^v{P}_{c,r}^e{P_{c,r}^e}^T-X_{c,a}^v{P}_{c,a,f}^e{P_{c,a,f}^e}^T\end{array},---\left(11\right)$

（8.2）计算出各子空间时间片监测指标后，则可得到参考模态与备选模态每个时间点k的监测指标：

参考模态：

$\begin{array}{c}{T}_{r,k,i}^2={(t_{r,k,i}-{\stackrel{\‾}{t}}_{r,k})}^T{\mathrm{\Σ}}_{r,k}^{-1}(t_{r,k,i}-{\stackrel{\‾}{t}}_{r,k})\\ \mathrm{SP}{E}_{r,k,i}=e_{r,k,i}^T{e}_{r,k,i}\end{array},---\left(12\right)$

其中，下标r表示参考模态，下标i表示每个时间片的第i个批次，表示从T_r,k计算得到的平均向量，∑_r,k是对角阵，每个元素表示时间片得分T_r,k的主元波动。

备选模态：

$\begin{array}{c}{T}_{a,f,k,i}^2={(t_{a,f,k,i}-{\stackrel{\‾}{t}}_{a,f,k})}^T{\mathrm{\Σ}}_{a,f,k}^{-1}(t_{a,f,k,i}-{\stackrel{\‾}{t}}_{a,f,k})\\ T_{a,o,k,i}^2={(t_{a,o,k,i}-{\stackrel{\‾}{t}}_{a,o,k})}^T{\mathrm{\Σ}}_{a,o,k}^{-1}(t_{a,o,k,i}-{\stackrel{\‾}{t}}_{a,o,k})\\ {T_{a,f,k,i}^e}^2={(t_{a,f,k,i}^e-{\stackrel{\‾}{t}}_{a,f,k}^e)}^T{{\mathrm{\Σ}}_{a,f,k}^e}^{-1}(t_{a,f,k,i}^e-{\stackrel{\‾}{t}}_{a,f,k}^e)\\ {T_{a,f,k,i}^e}^2={(t_{a,f,k,i}^e-{\stackrel{\‾}{t}}_{a,f,k}^e)}^T{{\mathrm{\Σ}}_{a,f,k}^e}^{-1}(t_{a,f,k,i}^e-{\stackrel{\‾}{t}}_{a,f,k}^e)\\ \mathrm{SP}{E}_{a,k}^f={e_{a,k,i}^f}^T{e}_{a,k,i}^f\end{array}---\left(13\right)$

其中，

分别表示T_a,f,k，T_a,o,k，

的平均向量，经过数据预处理后均未零向量。∑_a,f,k，∑_a,o,k，为对角阵，其对角元素分别为T_a,f,k，T_a,o,k，

的每个主元的波动值。

（8.3）计算控制限：

T²统计指标用来测量各采样时刻过程变量偏离正常模态下平均轨迹的距离，并且服从带权重的F分布，从而计算得到控制限；同理，SPE近似服从加权χ²分布，从而可计算得控制限。

步骤9：基于多模态时段模型的在线过程监测：

基于步骤6划分的时段、步骤7建立的时段模型监测系统和步骤8所得的监测统计量控制限可以在线监测注塑成型等新运行间歇过程的状态。该步骤由以下子步骤来实现：

（9.1）采集新测量数据及新测量数据预处理

在线监测时，采集到k时刻的新过程测量数据x_new(J×1)（其中下标new代表新样本）后，首先需要进行数据预处理。根据参考模态由步骤2所获得的均值和标准差，根据过程时间的指示调用对应该时刻的均值和标准差对现有数据如公式（1）中所示进行标准化预处理。

（9.2）计算新监测统计量

数据预处理后，根据公式（6）计算的PCA子时段模型，调用对应该新采样时刻所在时段的模型P_c,r(J×R_c,r)，按照如下方式计算得到主元得分，估计残差及其对应的Hotelling-T²与SPE两个监测统计指标：

$\begin{array}{c}{t_{r,\mathrm{new}}}^T={x_{\mathrm{new}}}^T{P}_{c,r}\\ {e_{r,\mathrm{new}}}^T={x_{\mathrm{new}}}^T{P}_{c,r}^e{P_{c,r}^e}^T\\ T_{r,\mathrm{new}}^2={(t_{r,\mathrm{new}}-{\stackrel{\‾}{t}}_{r,k})}^T{\mathrm{\Σ}}_{r,k}^{-1}(t_{r,\mathrm{new}}-{\stackrel{\‾}{t}}_{r,k})\\ \mathrm{SP}{E}_{r,\mathrm{new}}={e_{\mathrm{new}}}^T{e}_{\mathrm{new}}\end{array}---\left(14\right)$

其中,x_new是新的过程测量数据，

是前面根据训练数据获得的T_r,k的均值向量，∑_r,k是X_r,k的协方差矩阵。

（9.3）在线判断过程运行状态

比较两监测指标与各自的控制限，如果

和SPE_r,new监测指标均在控制限内，则过程可认为运行在参考模态正常状态。如果有超限报警发生，为了区分当前过程是处于备选模态正常运行状态还是故障状态，须按如下步骤进行区分：

（9.3.1）计算备选模态的新监测统计量

$\begin{array}{c}{t_{a,f,\mathrm{new}}}^T={x_{\mathrm{new}}}^T{P}_{c,a,f}\\ {t_{a,o,\mathrm{new}}}^T={x_{\mathrm{new}}}^T{P}_{c,a,o}\\ {t_{a,f,\mathrm{new}}^e}^T={x_{\mathrm{new}}}^T{P}_{c,a,f}^e{P_{c,a,f}^e}^T\\ {e_{a,\mathrm{new}}^f}^T={x_{\mathrm{new}}}^T{P}_{c,r}^e{P_{c,r}^e}^T-{x_{\mathrm{new}}}^T{P}_{c,a,f}^e{P_{c,a,f}^e}^T\\ {T_{a,f,\mathrm{new}}}^2={(t_{a,f,\mathrm{new}}-{\stackrel{\‾}{t}}_{a,f,k})}^T{\mathrm{\Σ}}_{a,f,k}^{-1}(t_{a,f,\mathrm{new}}-{\stackrel{\‾}{t}}_{a,f,k})\\ {T_{a,o,\mathrm{new}}}^2={(t_{a,o,\mathrm{new}}-{\stackrel{\‾}{t}}_{a,o,k})}^T{\mathrm{\Σ}}_{a,o,k}^{-1}(t_{a,o,\mathrm{new}}-{\stackrel{\‾}{t}}_{a,o,k})\\ {T_{a,f,\mathrm{new}}^e}^2={(t_{a,f,\mathrm{new}}^e-{\stackrel{\‾}{t}}_{a,f,k}^e)}^T{{\mathrm{\Σ}}_{a,f,k}^e}^{-1}(t_{a,f,\mathrm{new}}^e-{\stackrel{\‾}{t}}_{a,f,k}^e)\\ \mathrm{SP}{E}_{a,\mathrm{new}}^f={e_{a,\mathrm{new}}^f}^T{e}_{a,\mathrm{new}}^f\end{array}---\left(15\right)$

（9.3.2）判断过程运行状态

比较每个备选模态新的监测指标与控制限。如果所有新监测指标都处于某一个备选模态的控制限之内，则称新的样本与该模态匹配，即新过程运行在该备选模态正常状态下；相反，如果没有找到模态与新样本匹配，则新样本发生故障或者处于一种新的运行模态，该模态未包含在建模模型内。

Images