一种基于报警信度融合的多向主元分析过程监测方法与流程

本发明涉及属于间歇过程多变量监测与故障监测领域，涉及一种基于报警信度融合的多向主元分析过程监测方法。

背景技术：

多向主元分析方法(mpca)，虽然被应用于间歇过程的状态监测以及故障诊断上，但是该方法在实际应用中做了几方面的假设：1)间歇过程的数据是线性的，2)间歇过程始终处于稳定运行状态，过程数据符合统一分布，3)整个间歇生产过程是单一的，等等。实际的工业生产并不能满足上述几个假设，且mpca使用单一不变的统计控制限很容易造成报警错误，例如，基于95％的概率统计控制限太严格，导致统计误报太多，而99.999％的概率统计控制限又过于宽松，导致过程异常监测不到，所以基于概率统计的报警控制限在工厂应用中受到很大的阻碍。

技术实现要素：

本发明的目的，就是克服现有技术的不足，提供一种基于信息融合的批次过程监测报警方法，该方法可以大大提高多向主元分析在间歇过程监测过程应用中的报警的准确度，为了达到上述目的，采用如下技术方案：

本发明包括以下步骤得到：

1)模型数据采集

设一个间歇操作具有j个测量变量和k个采样点，则每一个测量批次可得到一个j×k的矩阵，重复i批次的测量步骤后，得到的数据可以表述为一个三维矩阵x(i×j×k)，其中测量变量为温度、速度、压力、行程等批次运行过程中可被测量的状态参数。

2)三维数据展开

将三维矩阵x(i×j×k)按照采集批次方向展开，即将一个操作批次内的各采样点上的变量按照时间顺序排开得到二维矩阵显然矩阵为i行kj列。

3)二维矩阵标准化

设二维矩阵内任意一点的变量为对该变量进行减均值、除方差的标准化处理，标准化处理的计算公式如下：

其中：是矩阵任一列的均值，sjk是矩阵任一列的方差

4)多向主元分析法建模

对上一步经标准化后的二维矩阵执行pca分解，完成多向主元分析法的建模，其中pca分解过程如下：

s＝trace(t^tt/(i-1))；(3)

其中：ti为正交的主元向量，pi为正交归一化的负载向量，s是主元的协方差矩阵的迹，代表各个主元对于过程的解释度大小；

公式(2)中x分解得到得分矩阵t(i×jk)及负载矩阵p(jk×jk)。

5)选取主元个数r

将公式(2)重新表述成如下形式：

其中：tr(i×r)、pr(jk×r)分别为保留r个主元后的得分矩阵和负载矩阵，e为残差矩阵；

通过上述变换，多向主元分析法模型将原始数据空间分解为主元空间和残差空间，主元空间变量高度相关，一般来说足以描述数据的变异性。

主元个数r一般可根据用户的经验设定或者采用broken-stick准则，broken-stick的内容是当第r个主元的解释度s(r)占所有主元总贡献sum(s)的百分比大于g(r)的时候保留该主元，否则终止，其中g(r)的计算公式如下：

其中：s(r)是第r个主元的解释度，sum(s)是所有主元的贡献和。

6)计算模型数据的两个指标t2和spe以及对应的95％和99.999％控制限

a)t2的95％控制限：

fr,n-r,α表示自由度为r，n-r的f分布，α为5％；

t2的99.999％控制限：

fr,n-r,α表示自由度为r，n-r的f分布，α为0.001％；

b)spe的95％控制限：

cα是正态分布校准水平α下的临界值，α为5％，cα＝norminv(0.95,0,1)；

spe的99.999％控制限：

其中θi和h0的计算同上，此处α为0.001％，cα＝norminv(0.99999,0,1)。

7)计算样本数据两个空间的指标spe和t2

样本残差空间的指标spe计算公式(10)，其中i为单位矩阵，

spe＝ee^t＝x(i-prpr^t)x^t(10)

样本主元空间的指标t2计算公式(11)，

t2＝tr^tsr^-1tr(11)

8)基于模糊阈获取历史报警信度

分别以t2和spe的99.999％和95％统计控制限为模糊上下限，构造模糊隶属度函数来将统计量超过控制阈值的报警转化报警信度，代表spe或者t2的99.999％的概率统计控制限，代表spe或者t2的95％的概率统计控制限；

将步骤7)计算的spe、t2以及步骤6)中两个统计量的上下限分别带入函数得到两个报警信度序列bspe{bspe(1),bspe(2)……bspe(t)}，bt2{bt2(1),bt2(2)……bt2(t)},t为采样的时刻；

9)基于报警信度计算组合权重

首先计算报警信度序列bspe(t){bspe(1),bspe(2)……bspe(k)}的组合权重

(a)当t＝1时，因在t＝1之前没有相关的报警信息，所以t＝1时全局的报警信度是就是第一时刻的报警信度b0:1(na)＝bspe(1)；

(b)当t＝2时，全局的报警信度利用t＝1时刻的全局报警信度和t＝2时刻的报警信度进行加权融合，得到t＝2时刻的全局报警信度b0:2(na)，其计算如式所示：

b0:2(na)＝α2b0:1(na)+β2bspe(2)(13)

其中α2＝0.5表示对b0:1(na)的线性加权值，β2＝0.5表示当前时刻报警信度的加权值；

(c)当t≥3时,全局的报警信度利用前两时刻的全局报警信度和当前时刻的报警信度进行加权融合，

b0:t(na)＝αtb0:t-1(na)+βtbspe(t)(14)

其中权重系数αt表示上一时刻全局信度的权重因子，βt是当前时刻报警信度的加权值，二者满足关系式αt＝1-βt。其公式如下(15)-(16)

其中supt(bi)代表报警信度的支持度，其计算方法如下：

(i)当t≥3之后，可以得到t时刻的报警信度bspe和t-1以及t-2时刻的全局报警证据b0:t-1和b0:t-2，计算b0:t-1以及b0:t-2两两之间的距离采用欧式距离公式：

那么b0:t-1以及b0:t-2报警信度之间相似度，

sim(b0:t-1,b0:t-2)＝1-dst(b0:t-1,b0:t-2)(18)

同理可得，sim(b0:t-2,bspe)，sim(b0:t-1,bspe)

(ii)计算bspe和b0:t-1以及b0:t-2每个报警信度相对于其他两个报警信度的支持程度，所以每个证据的支持度(19)-(21)式所示：

supt(bt)＝sim(b0:t-1,bspe)+sim(b0:t-2,bspe)(19)

supt(b0:t-1)＝sim(b0:t-1,bspe)+sim(b0:t-2,b0:t-1)(20)

supt(b0:t-2)＝sim(b0:t-1,b0:t-2)+sim(b0:t-2,bspe)(21)

同理，可得报警信度序列bt2的组合权重

10)基于全局报警信度的报警决策

当b0:t(na)≥0.5则监测报警。

与现有技术相比，本发明的有益效果在于：传统的多向主元分析法模型用固定的统计控制限来判断过程异常与否，容易造成报警失误，本发明综合历史时刻与当前时刻的监测信息来进行报警决策，提高了过程监测故障诊断结果的准确性，为无过程先验知识条件下的过程监测提供了新的实用方法。

附图说明

图1是本发明所述模糊隶属度函数

图2是本发明所述基于多向主元分析法的批次过程的三维数据展开示意图。

图3是本发明具体实施例中注塑过程的某一批次的mpca监测结果。

具体实施方式

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

注塑成型是典型的批次过程，其一个批次主要包含注射、保压、塑化、冷却四个阶段，具体来说，在注射段，液压缸的高压推动螺杆向前将机桶内的熔融塑料推到模腔，当模腔被完全或者将近填充满的时候，过程切换到保压阶段，在该阶段，高压将少量材料继续填充到模腔中以补充由于冷却和固化带来的材料收缩；当胶口冷却，模腔中的材料不再被注射喷嘴影响的时候，保压段结束。螺杆旋转并后退，将足够量的熔融塑料推到螺杆前端。螺杆后退同时开始容积计算。头部熔料达到一定的注射量后，螺杆停止后退和转动，这段时间的过程状态称为塑化阶段。在保压段结束，塑化过程进行的时候，冷却阶段也同时进行着，直到模具内材料达到能被弹出的硬度，冷却阶段结束。

注塑过程通过不断重复这一过程来获得稳定一致的产品，以上述注塑成型过程为例，本发明所述的基于报警信度的多元统计监测方法，包括以下步骤得到：

(1)模型数据采集

设一个间歇操作具有j个测量变量和k个采样点，则每一个测量批次可得到一个j×k的矩阵，重复i批次的测量步骤后，得到的数据可以表述为一个三维矩阵x(i×j×k)。为了确保检测数据涵盖足够长时间的工作范围，一般工业上用来建模的数据批次i的取值大于100，测量变量为温度、速度、压力、行程等批次运行过程中可被测量的状态参数；基于过程时间长短、过程变化的快慢程度以及模型负担是否在合理的范围，采样点k个数一般小于1000。

本实施例中，测量变量实验室注塑机工作过程可获得的变量为8个：压力阀门开度，流量阀门开度，注射行程，注射速度，注射压力，机桶温度(3段)，操作批次i取100，每个批次保留的采样点k为488。

(2)三维数据展开

参见图2，将三维矩阵x按照采集批次方向展开，即将一个操作批次内的各采样点上的变量按照时间顺序排开得到二维矩阵矩阵为i行kj列；

(3)二维矩阵标准化

设二维矩阵内任意一点的变量为对该变量进行减均值、除方差的标准化处理，标准化处理的计算公式如下：

其中：是矩阵任一列的均值，sjk是矩阵任一列的方差，

本步骤的标准化处理相当于抽取了间歇过程中一次操作的平均运行轨迹，突出了间歇过程不同操作批次之间的一种正常随机波动。

(4)mpca建模

所谓的mpca建模就是先将三维矩阵展开成一个大的二维矩阵，再执行常规的pca分解的方法，本步骤对上一步经标准化处理后的二维矩阵x(i×jk)执行pca分解，其分解过程如下：

s＝trace(t^tt/(i-1))；(3)

其中：ti为正交的主元向量，pi为正交归一化的负载向量，s是主元的协方差矩阵的迹，代表各个主元对于过程的解释度大小。

公式(2)x分解得到得分矩阵t(i×jk)及负载矩阵p(jk×jk)。

(5)选取主元个数

一般来说，前几个主元一般包含着间隙过程的大部分变异信息，其他的主元可能主要包含噪声信息，因此公式(2)可以被重新表述成如下形式：

其中：tr(i×r)、pr(jk×r)分别为保留r个主元后的得分矩阵和负载矩阵，e为残差矩阵；

通过上述变换，mpca模型将原始数据空间分解为主元空间和残差空间，主元空间变量高度相关，一般来说足以描述数据的变异性。

主元个数r一般可根据用户的经验设定或者采用broken-stick准则，broken-stick的内容是当第r个主元的解释度s(r)占所有主元总贡献sum(s)的百分比大于g(r)的时候保留该主元，否则终止，其中g(r)的计算公式如下：

其中：s(r)是第r个主元的解释度，sum(s)是所有主元的贡献和,在本实施例中，主元r的个数选择3，对于过程的解释度为86.73。

(6)计算模型统计控制限

a)t2的95％控制限：

fr,n-r,α表示自由度为r，n-r的f分布，α为5％。

t2的99.999％控制限：

fr,n-r,α表示自由度为r，n-r的f分布，α为0.001％。

将r＝3，n＝100，带入上述公式以后，得到tsqcl95＝8.3447，tsqcl999＝30.4307；

b)spe的95％控制限：

cα是正态分布校准水平α下的临界值，α为5％，cα＝norminv(0.95,0,1)＝1.6449。

spe的99.999％控制限：

其中θi和h0的计算同上，此处α为0.001％，cα＝norminv(0.99999,0,1)＝4.2649。

(7)计算样本数据两个空间的统计量

spe＝ee^t＝x(i-prpr^t)x^t(10)

t2＝tr^tsr^-1tr(11)。

(8)基于模糊阈获取历史报警信度

利用公式(10)-(11)得到残差空间与主元空间两个报警信度序列bspe{bspe(1),bspe(2)……bspe(t)}，bt2{bt2(1),bt2(2)……bt2(t)},t为采样的时刻，其模糊隶属度函数可见图1。

(9)基于报警信度分别计算组合权重

首先计算报警信度序列bspe(t){bspe(1),bspe(2)……bspe(k)}的组合权重；

(a)当t＝1时，因在t＝1之前没有相关的报警信息，所以t＝1时全局的报警信度是就是第一时刻的报警信度b0:1(na)＝bspe(1)

(b)当t＝2时，全局的报警信度利用t＝1时刻的全局报警信度和t＝2时刻的报警信度进行加权融合，得到t＝2时刻的全局报警信度b0:2(na),其计算如式所示：

b0:2(na)＝α2b0:1(na)+β2bspe(2)(12)

其中α2＝0.5表示对b0:1(na)的线性加权值，β2＝0.5表示当前时刻报警信度的加权值；

(c)当t≥3时,全局的报警信度利用前两时刻的全局报警信度和当前时刻的报警信度进行加权融合，

b0:t(na)＝αtb0:t-1(na)+βtbspe(t)(13)

其中权重系数αt表示上一时刻全局信度的权重因子，βt是当前时刻报警信度的加权值，其计算方法如下：

(i)当t≥3之后，可以得到t时刻的报警信度bspe和t-1以及t-2时刻的全局报警证据b0:t-1和b0:t-2，计算b0:t-1以及b0:t-2两两之间的距离采用欧式距离公式：

那么b0:t-1以及b0:t-2报警信度之间相似度，

sim(b0:t-1,b0:t-2)＝1-dst(b0:t-1,b0:t-2)(15)

同理可得，sim(b0:t-2,bspe)，sim(b0:t-1,bspe)

(ii)计算bspe和b0:t-1以及b0:t-2每个报警信度相对于其他两个报警信度的支持程度，所以每个证据的支持度(16)-(18)式所示：

supt(bspe)＝sim(b0:t-1,bspe)+sim(b0:t-2,bspe)(16)

supt(b0:t-1)＝sim(b0:t-1,bspe)+sim(b0:t-2,b0:t-1)(17)

supt(b0:t-2)＝sim(b0:t-1,b0:t-2)+sim(b0:t-2,bspe)(18)

(iii)步骤的基础上得到每个报警信度相对于与其他两个报警信度的支持度，最后得到动态更新的权重因子βt，如下式(19)；

αt＝1-βt

同理，可以计算主元空间信度序列bt2{bt2(1),bt2(2)……bt2(t)}的报警信度组合权重。

(10)基于全局报警信度的报警决策

当b0:t(na)≥0.5则监测报警。

本实施例的结果如下，可见图3：

本发明综合历史时刻与当前时刻的报警信度来进行报警决策，提高了过程监测故障诊断结果的准确性有效性。

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