一种基于子时段MPCA-SVM的间歇过程故障诊断方法

摘要

一种基于子时段MPCA-SVM的间歇过程故障诊断方法涉及基于模式识别的故障诊断领域。本发明首先对发酵过程的三维数据进行展开，并沿时间方向切片；其次利用MPCA对间歇过程进行时段粗划分和细化分；最后在每个子时段内建立MPCA监测模型和SVM诊断模型。在线故障诊断包括：对采集到的数据按照模型进行处理，计算其统计量并与控制限进行比较。若未超限，则生产正常运行；若超限，则将数据代入到相应时段的SVM诊断模型中进行故障诊断。本发明只填充发生故障时段的数据，减小了由于人为的填充过多的未知数据对SVM故障诊断的准确率带来的影响。同时，也减少建模的数量，从而解决了因频繁地更新模型而造成的诊断过程复杂的问题。

说明书

技术领域

本发明涉及基于模式识别的故障诊断技术领域，特别涉及一种针对间歇 过程在线故障诊断技术。本发明的基于模式识别的方法即是在典型的间歇过 程—青霉素发酵过程故障监测方面的具体应用。

背景技术

间歇过程因其产品具有特定功能、高附加值、小批量、多品种等特点， 使得其在生产中所占的比重越来越大，因此对间歇过程的故障诊断也越来越 重要。但间歇过程具有动态性、强非线性和时段特性等特点，且操作复杂度 高和产品质量易受到环境等因素的影响，使得对其故障诊断的研究面临更大 的挑战。

对于间歇过程的故障诊断，目前比较常用的方法有贡献图和模式识别的 方法，有些学者利用多向主成分分析(Multi-wayprincipalcomponent analysis,MPCA)对间歇过程进行分时段在线监控，并利用贡献图的方法追溯 故障变量，但该方法采用正常的数据进行故障诊断，不能够真正反应故障的 信息，且忽略了变量之间的相关性，只能诊断单变量故障，而模式识别的方 法可以克服MPCA贡献图的不足。支持向量机(SupportVectorMachine,SVM) 因其在小样本下具有较强的非线性系统学习能力，使得SVM成为一种被广泛 应用的模式识别方法。有些学者利用MPCA与SVM对整个间歇过程建立离线故 障诊断模型，不可避免的在诊断阶段需人为的填充过程的未知数据，但填充 的数据往往与真实的过程数据存在差异，在一定程度上会影响到故障诊断的 准确率；为了解决这个问题，另外一些学者将滑动时间窗技术与SVM结合应 用于间歇过程的故障诊断，虽可解决因大量填充数据带来诊断准确率的问题， 但存在模型更新频繁的缺点，使得诊断过程繁琐。

发明内容

为了克服现有技术的不足，本发明提供了一种基于子时段MPCA-SVM的间 歇过程在线故障诊断方法。利用MPCA划分间歇过程，将故障诊断的搜索空间 局限到一个个特定的时段内，并对每个时段建立故障监测和诊断模型。该方 法无需填充整个生产过程的未知数据，而只填充发生故障时段的数据，因此 可减小由于人为的填充过多的未知数据对SVM故障诊断的准确率带来的影响。 同时，也可以减少建模的数量，从而解决了因频繁地更新模型而造成的诊断 过程复杂的问题。

本发明采用了如下的技术方案及实现步骤：

A.离线建模阶段

1)采集发酵过程正常工况下的历史数据，所述的历史数据X由离线测试 得到的同一发酵过程相同工艺下的I批次数据构成，X＝(X₁,X₂,...,X_I)^T，其中 X_i(i＝1,2,…,I)表示第i批次数据。每个批次包含K个采样时刻，即 X_i＝(X_i,1,X_i,2,...,X_i,K)，其中X_i,k表示第i批次第k采样时刻采集的数据。每个采 样时刻采集J个过程变量，即X_i,k＝(x_i,k,1,x_i,k,2,...,x_i,k,J)，其中x_i,k,j表示第i批次中 第k采样时刻的第j个过程变量的测量值；

2)对历史数据X进行标准化处理，处理方式如下：

首先计算历史数据X的所有时刻上所有过程变量的均值和标准方差，其 中第k采样时刻的第j个过程变量的均值的计算公式为其中 x_i,k,j表示第i批次中第k采样时刻的第j个过程变量的测量值，k＝1,...,K， j＝1,...,J；第k采样时刻的第j个过程变量的标准方差s_k,j的计算公式为， $s_{k,j}=\sqrt{\frac{1}{I-1}\underset{i=1}{\overset{I}{\Sigma }}{(x_{i,k,j}-{\bar{x}}_{k,j})}^2},$k＝1,...,K，j＝1,...,J；

然后对历史数据X进行标准化，其中第i批次中第k采样时刻的第j个过程 变量的标准化计算公式如下：

${\tilde{x}}_{i,k,j}=\frac{x_{i,k,j}-{\bar{x}}_{k,j}}{s_{k,j}}---\left(1\right)$

其中，i＝1,...,I，j＝1,...,J，k＝1,...,K；

3)X通过步骤2)的标准化得到新的二维矩阵X'，该矩阵共有(K×J)个列 向量，即X'＝(X₁',X'₂,...,X'_K×J)，其中第j个列向量X'_j＝(X'_j,1,...,X'_j,K)^T， X'_j,k＝(X'_j,k,1,...,X'_j,k,I)^T，其中X'_j,k,i表示经过步骤2)标准化处理后的第j个过程变 量第k个采样时刻在第i个批次中对应的值，其中i＝1,...,I，j＝1,...,J，k＝1,...,K；

4)利用多向主元分析MPCA方法提取X'中每个时刻的主成分。若提取第k 个采样时刻数据X'_k主成分，具体的步骤如下：

4.1)求出二维矩阵X'_k的协方差矩阵COV；

$COV=\frac{1}{I-1}{X}_k^{\prime T}{X}_k^{\prime }---\left(2\right)$

4.2)对矩阵COV的特征值分解；

COV＝VΛV^T(3)

式(3)揭示了协方差矩阵的关联关系，其中Λ＝(λ₁,λ₂,…,λ_v)为对角阵，v为二维矩 阵X'_k的特征值数，Λ包含m个幅值递减的非负实特征值(λ₁≥λ₂≥…≥λ_m≥0)。V 是正交阵(V^ΤV＝E,E为单位阵)。

定义上，通过MPCA方法，矩阵X'_k等价于X'_k＝QP^Τ+R，并将其代入到式 (3)得到：

$S=\frac{1}{I-1}{\mathrm{PQ}}^{T}QPQ---\left(4\right)$

其中，Q为得分矩阵，P为负载矩阵，R为残差矩阵。

对应式(3)和(4)各项，可得：

$\begin{array}{c}P=V\\ \Lambda =\frac{1}{I-1}{Q}^{T}Q\end{array}---\left(5\right)$

若前A个主成分的累加方程超过一个阈值0.85，那么就提取前A个主元 来作为综合指标，则原始的J维空间就变为A维，且A≤J。

4.3)求出得分矩阵Q；

为了最优地获取数据的变化量，同时最小化随机噪声对PCA产生的影响， 保留与A个最大特征值相对应的负载向量，则X'_k在低维空间的投影信息就包 含在得分矩阵：

Q＝X'_kP(6)

5)利用MPCA对间歇过程进行时段划分：

5.1)时段粗划分：

当相邻采样点具有相同的主成分个数(Q矩阵的列数)时，就把这些采样点 划分到同一个时段，若具有相同主成分个数的时段长度L小于整个间歇过程的 1/10时，将此时段数据归到左右相邻的主成分个数相差相对小的一段内。若 此时段长度小于整个发酵过程的1/10且该时段的主成分个数与相邻的两个时 段的主成分个数差值相等时，则分别取前后时段主成分个数为该段的主成分 个数并计算其贡献率，把该段归到贡献率变化相对小的一个时段内。

最终粗划分可以得到F个时段，分别表示为S₁,S₂,…,S_F。

5.2)时段细划分：在任意粗划分的时段S_f(f＝1,2,…,F)内，利用负载矩 阵之间角度信息和距离信息来定义相似性度量公式(7)。

$d_{b,c}^{S_f}=\underset{l=1}{\overset{a_f}{\Sigma }}{\gamma }_l\frac{\left|\right|P_{bl}^T{P}_{cl}\left|\right|}{\left|\right|P_{bl}\left|\right|·\left|\right|P_{cl}\left|\right|}+(1-{\gamma }_l)e^{-\left|\right|P_{bl}-P_{cl}\left|\right|}---\left(7\right)$

${\gamma }_l=\frac{1}{l}/\underset{h=1}{\overset{a_f}{\Sigma }}\frac{1}{h},l=h=1,2,...,a_f---\left(8\right)$

其中b和c表示在时段S_f内任意相邻的两个采样时间，a_f为时段S_f内的 采样点个数，l和h表示采样点1到a_f中的任意一个值，||P_bl||和||P_cl||分别是b 和c采样时刻的负载矩阵的模值，表示b和c采样时刻负载矩阵的距离， r_l为加权系数，用以强调不同投影方向的不同重要性，表示两个时间片负 载矩阵中a_f个投影方向的夹角余弦值的加权和，表示b和c采样时刻负载 矩阵的相似度。

最终细化分得到G个时段，分别表示为S₁',S'₂,...,S'_G，并计算任意一个子时 段S'_g(g＝1,2,…,G)得均值负载矩阵

6)在S₁',S'₂,...,S'_G中建立各子时段MPCA监控模型，如式(9)。

对于每个子时段均具有相近的负载矩阵，因此各个子时段的MPCA模型均 可采用均值负载矩阵来描述。对于任意子时段S'_g中任意第t时刻的MPCA模 型可以表示为：

$\left(\begin{array}{c}{q}_t=S_g^{\prime }{\bar{P}}_{S_g^{\prime }}\\ \widehat{S_g^{\prime }}=q_t{\left({\bar{P}}_{S_g^{\prime }}\right)}^T\\ r_t=S_g^{\prime }-\widehat{S_g^{\prime }}\end{array}\right)---\left(9\right)$

其中q_t是子时段S'_g中第t时刻的负载向量，子时段S'_g的特征向量，r_t是 残差向量。

7)利用MPCA建立监控模型时，必须先要确定两个控制限，控制限用来 判断当前数据是否处于正常运行状态，这两个控制限分别称为T²统计量控制 限和SPE统计量控制限，计算T²统计量的控制限满足F分布式，如式(10)， 其中β为主成分个数，I为建模的批次数，α为显著性水平。

$T_{t,\alpha }^2:\frac{\beta (I-1)}{I-\beta }{F}_{\beta ,I-\beta ,a}---\left(10\right)$

SPE统计量的控制限由式(11)计算，其中mean_t，var_t分别表示建模数据 中每个子时段批次测量数据在t时刻SPE的均值和方差。

$\left(\begin{array}{c}{\mathrm{SPE}}_{t\alpha }~\frac{{\mathrm{var}}_t}{2{\mathrm{mean}}_t}{\chi }_{t,h_t,\alpha }^2\\ h_t=\frac{2{\mathrm{mean}}_t^2}{{\mathrm{var}}_t}\end{array}\right)---\left(11\right)$

8)在每个子时段内，采用1)和2)中的方式，分别取不同类型故障运行 状态下的数据各I_fault组，利用3)方法标准化处理，并利用4)中方法提取出各 个类型的故障数据的主成分Q₁,Q₂,…,Q_fault；

9)将提取出的主成分Q₁,Q₂,…,Q_fault作为支持向量机模型的输入，在每个子 时段内建立一个故障诊断模型，最后可得到G个模型，分别是SVM₁，SVM₂，…， SVM_G，其中每个SVM模型都采用一对一的多分类器构造方式。

9.1)任意一个SVM_g(g＝1,2,…,G)模型都是一个决策超平面。

9.2)一对一的多分类器构造方法，该方法针对u类训练样本，将不同类 样本两两组合用来训练个二分类器。然后将测试样本代入所有分类器， 采用投票法进行决策，如果分类器SVM_1,2认为样本属于第1类，则在第1类的 投票上加1，否则在第2类的投票上加1，最后统计各类的票数，得票最多的类 即为测试样本所属类别。

B.在线故障诊断阶段

10)采集当前时刻数据X_t(1×J)，并判断当前时刻所处的操作子时 S'_g(g＝1,2,…,G)段；

11)按照式(9)计算X_t(1×J)所对应的得分向量Q_t(1×β)以及预测误差向 量R_t(1×J)，式中E为(J×J)的单位矩阵；

$\left(\begin{array}{c}{Q}_t=X_t{\bar{P}}_{S_g^{\prime }}\\ R_t=X(E-{\bar{P}}_{S_g^{\prime }}{\left({\bar{P}}_{S_g^{\prime }}\right)}^T)\end{array}\right)---\left(12\right)$

12)计算X_t(1×J)对应的和SPE_t统计量，计算公式如式(10)；

$\left(\begin{array}{c}{T}^2=Q_t{\left(S_g^{\prime }\right)}^{-1}{Q_t}^T\\ SPE=R_t{R_t}^T\end{array}\right)---\left(13\right)$

其中，R_t是矩阵X_t(1×J)的预测误差向量。

13)判断此时的或SPE_t统计量否超出控制限T²或SPE，若两者均未超 出，则当前过程测量数据正常，返回1)；否则当前操作过程有故障，进行14)；

14)在子时段S'_g内，对新批次未知的数据进行数据填充，填充数据为当 前时刻X_t(1×J)，最后可以得到该时段的数据矩阵为

15)对数据矩阵标准化处理并提取出该数据矩阵的主成分，之后代 入到该子时段的SVM_g(g＝1,2,…,G)模型中，进行故障诊断。

有益效果

与现有技术相比，本发明利用MPCA划分间歇过程，在每个子时段建立故 障监测和诊断模型，缩小了故障诊断的搜索空间。运用该方法在线故障诊断 时，无需填充整个过程的未知数据，而只是填充发生故障时段的未知数据， 这样可减小过多的填充未知数据带来的影响，从而可提高故障诊断的准确率。 同时无需频繁的更新故障诊断模型，而只需选出发生故障时段的诊断模型， 降低了故障诊断的复杂性。

附图说明

图1为本发明的离线建模流程图；

图2为在线故障诊断流程图；

图3为正常数据X(30×10×400)组成形式示意图；

图4子时段划分结果；

图5子时段5中MPCA对故障批次的T²和SPE监控图；

图6不分时段SVM对测试数据分类的效果图；

图7本发明对测试数据分类的效果图；

具体实施方式

美国Illinois州立理工学院过程监控与技术小组开发的青霉素仿真平台 PenSim2.0，为间歇过程的监控和故障诊断提供了一个标准平台，该平台已 成为国际上较有影响力的青霉素仿真平台。

本发明以该平台为仿真研究对象，设置青霉素发酵每个批次的反应时间 为400h，采样间隔为1小时，选取10个过程变量进行仿真研究，如表1所示。 同时，该平台可以设定三种故障：1.空气流量，2.搅拌功率，3.底物流加速 率。故障扰动的类型有阶跃扰动和斜坡扰动两种，可进一步设定两种扰动的 幅度、扰动的引入时间和终止时间。

表1建立模型所用变量

本实验仿真30批正常数据X(30×10×400)作为划分时段和建立MPCA故 障监测模型的数据库，其中每批次的初始条件略有变化，同时对所有测量变 量均加入了随机测量噪声。

另外，也仿真了故障1(变量1阶跃故障数据)、故障2(变量2阶跃故障 数据)和故障3(变量3阶跃故障数据)各30批用于训练SVM故障诊断模型，其 中每个故障选用15批用来建立故障诊断模型，剩下的15批用于测试该诊断 模型。为了使每一种故障数据包含更全面的故障信息，同一种故障的15批训 练数据分别是在不同的反应时间引入扰动且均延续到反应结束。

将本发明方法应用到上述青霉素发酵过程仿真对象包括离线建模和在线 诊断两大步骤，具体陈述如下：

A.离线建模阶段：

步骤1：将上述30批正常数据X(30×10×400)沿批次方向展开，具体形 式见图2。图中为一个30行4000列的矩阵。

步骤2：对X(30×10×400)进行标准化处理。

首先按公式计算第k采样时刻的第j个过程变量在所有批 次上的均值，其中x_i,k,j为X_30×10×400第i批次中第k采样时刻的第j个过程变量的测 量值，k＝1,...,400，j＝1,...,10；第k采样时刻的第j个过程变量的标准方差s_k,j的 计算公式为，$s_{k,j}=\sqrt{\frac{1}{30-1}\underset{i=1}{\overset{30}{\Sigma }}{(x_{i,k,j}-{\bar{x}}_{k,j})}^2},$k＝1,...,400，j＝1,...,10。

之后对X(30×10×400)进行标准化，其中第i批次中第k采样时刻的第j 个过程变量的标准化计算公式如下：

${\tilde{x}}_{i,k,j}=\frac{x_{i,k,j}-{\bar{x}}_{k,j}}{s_{k,j}}$

其中，i＝1,...,30，j＝1,...,10，k＝1,...,400；

步骤3：X(30×10×400)通过步骤2)标准化后得到二维矩阵为 X'(30×(10×400))；

步骤4：将X'(30×(10×400))沿时刻方向划分为400个时间片矩阵 X_k'(30×10)，k＝1,...,400。

步骤5：利用MATLAB中主成分分析函数(princomp)提取出二维矩阵 X_k'(30×10)的主成分Q_k。

步骤6：当Q_δ＝Q_θ(δ＝θ+1,δ＜θ且θ＝2,…,400)时,δ和θ即可划分到同一 个时段，若时段长度L达不到采样点时，可以将此段数据归到左右 相邻的主成分个数相差相对小的一段内。若时段长度L少于采样点 且与相邻的主成分个数差值相等时，则分别取该时段的前后两个时段的主成 分个数为该时段的主成分个数并计算其贡献率，把该时段归到贡献率变化相 对小的一段内。最终划分出4个粗时段S_f(f＝1,…,4)。

步骤7：在时段S_f内，利用公式(7)，当计算出来大于0.8时，就把b 和c采样时间放到同一个细化分的时段内。最终将间歇过程划分为7个子时 段S_g(g＝1,…,7)，分别为时段1(1-37h)，时段2(38-85h)，时段3(86-174h)， 时段4(175-215h)，时段5(216-296h)，时段6(297-377h)，时段7(378-400h)。

步骤8：分阶段离线建模

该部分以时段5为例，对子时段的MPCA监控模型和SVM故障诊断模型建 立进行描述。

步骤8.1：取出X'(30×(10×400))中第5时段相对应的数据建立时段5的 MPCA模型，并且计算出时段5的均值负载矩阵

步骤8.2：利用式10)和11)确定时段5的T²和SPE统计量控制限，控制 限图如图5。

步骤8.3：利用步骤1，步骤2，步骤3和步骤4的方法，对时段5中故 障1,2和3数据进行预处理，标准化后得到矩阵X_5f1,X_5f2和X_5f3。

步骤8.4：利用Matlab中的主成分分析函数(princomp)对X_5f1,X_5f2和 X_5f3进行分解，可得到主成分Q_5f1,Q_5f2和Q_5f3。

步骤8.5：利用Matlab中libsvm工具箱的svmtrain函数对Q_5f1,Q_5f2和 Q_5f3进行训练，最终可得到决策超平面F₅(□)，该决策超平面就是故障诊断模型 SVM₅。

同理，按照步骤8.2到8.5的方法，建立其他时段的故障诊断模型，最 终得到模型SVM₁，SVM₂，SVM₃，SVM₄，SVM₆，SVM₇。

B.在线诊断阶段：

该部分以变量1在220个时刻引入5％的阶跃故障为测试数据。

步骤9：采集当前发酵过程第t采样时刻的10个过程变量的数据x_t，并根 据步骤2中得到的t时刻的均值和标准方差对其进行标准化得到 其中第t采样时刻的第j个过程变量的标准化公式如下：

${\tilde{x}}_{t,j}=\frac{x_{t,j}-{\bar{x}}_{t,j}}{s_{t,j}}$

其中，x_t,j为当前第t采样时刻所采集发酵数据中的第j个过程变量，为第t采样时刻的第j个过程变量的平均值，s_t,j为第t采样时刻的第j个过程变 量的标准方差，j＝1,...,10；

步骤10：依据公式15)，计算标准化后的t时刻采集数据的对应的主成 分向量Q_5t和预测误差向量R_5t。

步骤11：根据步骤8.3中的方法计算t时刻和SPE_t统计量。

步骤12：将上述计算得到的监控统计量和SPE_t与建模步骤8.3确定的 控制限进行比较，如果无超限，即为正常，则采集下一个时刻的数据并返回 步骤9,10,11和12，直至发酵过程完毕；否则认为发生故障。

步骤13：若有故障，则首先判断故障所发生的阶段为第5时段，并调取 出该时段的故障诊断模型SVM₅；

步骤14：在第5时段内，利用当前时刻采样数据x_t将未知的数据填充满， 则得到该时段的数据矩阵为X₅；

步骤15：根据步骤2中的数据标准化方法和步骤5中PCA提出主成分方 法对X₅进行数据预处理，最终可得到主成分Q₅。

步骤16：在Matlab中，svmpredict函数利用步骤8.5中svmtrain训练 的模型model参数，计算Q₅故障数据所被判别的标签为1和判别的准确率 accurate为86.67％。

上述步骤即为本发明方法在青霉素发酵仿真平台故障监测和诊断领域的 具体应用。为了验证本发明的有效性，对三种故障各15批数据进行实验。图 5是故障发生在第5时段的故障监控图，其他时段的故障监控图与图5相似这 里不再具体说明。图中与横坐标平行的线表示通过核密度估计方法确定的控 制限，曲线为实时监控值。若曲线高于控制限说明有故障发生，否则为正常 运行状态。

图6和图7分别表示故障发生在子时段5时，在不分阶段方法和本发明 方法下的故障诊断效果图，图中圆圈代表实际测试集的分类，星号代表预测 测试集分类。在同一个测试集样本下，当两种标志相重合，那么说明测试结 果与实际结果相同，否则测试结果错误。

为了更加有效的证明本发明方法应用于发酵过程故障诊断上的 优越性，本发明又分别对其他时段进行了故障诊断研究，诊断结果列表对比 如下：

表2各个时段故障诊断结果

从上表不难发现，当故障发生在前几个时段时使用本发明方法可以大大 的提高故障的准确率，但随着生产的继续，此方法虽比传统的方法诊断效果 好，但是准确率提高的幅度确实有限的。从整体来看，本发明方法可以明显 的提高故障准确率。