CN102662321B - 一种主元分析监控模型的在线更新方法 - Google Patents

Abstract

本发明涉及一种主元分析监控模型的在线更新方法，包括以下步骤：1)在工业现场设置一包括数据采集设备和监控计算机的模型在线更新***；2)传统PCA建模模块利用历史数据建立PCA初始监控模型；3)监控开始后，均值方差更新模块根据实时过程数据以及现有PCA模型计算出新模型中的均值

和标准差σ′；4)投影点计算模块计算出新样本的残差向量，传送给残差判定模块；5)残差判定模块根据残差向量模的大小决定投影方向的更新方法：如果残差较大，则调用主元空间调整模块；如果残差值较小，则调用主元方向微调模块；最后得到新模型的负荷向量P′_nk和特征值矩阵Λ′_kk；6)控制限更新模块对模型的统计量控制限

和进行更新；***最终输出新模型Ω′，用于工业过程中的在线监控及故障诊断。

Description

一种主元分析监控模型的在线更新方法

技术领域

本发明涉及一种用于多变量统计过程监控模型的在线更新方法，特别是关于一种基于增量主元分析方法的多变量统计过程监控模型在线更新方法。

背景技术

PCA(Principal Components Analysis，主元分析)是一种多变量统计过程监控建模方法，通过挖掘过程变量之间的线性相关性，建立反映***内在规律的监控模型，提供给多变量统计过程监控***，从而实现对生产过程的有效监控。传统的基于PCA的监控模型建模过程首先需要收集大量能代表生产过程特征的生产数据，然后基于这些历史数据建立监控模型。然而，在实际生产过程中，催化剂退化、设备老化积灰等外界环境因素容易造成过程变量的缓慢漂移现象。这将改变过程变量的均值、标准差等数据分布属性，使得基于历史数据建立的PCA监控模型与生产过程的实际情况逐渐产生偏差，导致监控模型失效，致使监控***产生大量的虚警。

现有的能够进行模型在线更新以适应过程变量缓慢漂移的PCA监控方法大致有两类。第一类方法将指数加权滑动平均(Exponentially Weighted MovingAverage，EWMA)滤波器与PCA模型相结合，仅对过程变量的均值和标准差进行更新。该方法的优点是计算速度较快，存储需求较低；缺点是无法更新主元投影方向等PCA模型结构信息，使得监控模型与实际生产过程之间仍然存在偏差。第二类方法能够在线更新整个PCA模型，包括变量的均值、标准差、主元投影方向和主元个数等。其优点可以随时保证监控模型与实际生产过程的一致性，缺点是对协方差矩阵的递推更新和“特征值分解”计算将导致较高的计算复杂度，以及较高的存储需求。综上，现有的自适应模型监控方法或者是监控效果不够精确，或者是计算复杂性较高，都存在着一定的不足之处。

近年来，有一种应用于模式识别领域的在线数据压缩算法——增量PCA(Incremental PCA)方法引起了一些学者的关注。该方法具有计算效率高，且无需额外存储空间的优点，并且能够在线更新模型的均值、标准差以及主元投影方向。但将该方法应用到多变量统计过程监控时，数据的时变性将导致监控模型的主元个数的不断增加，甚至趋近于与变量个数相同，使得SPE(Squared PredictionError，平方预测误差)统计量近乎失效，因而无法用于实际监控。

发明内容

针对上述问题，本发明的目的是提供一种应用增量主元分析算法的统计过程监控模型在线更新方法，该方法能够实时更新模型的均值、标准差以及主元投影方向，对工业过程中的过程变量缓慢漂移现象具有鲁棒性，还具有可靠性高，计算速度快，对***的存储需求低等优点，适用于实际生产过程多变量统计监控***的模型在线更新。

为实现上述目的，本发明采取以下技术方案：一种主元分析监控模型的在线更新方法，其包括以下步骤：1)在工业现场设置一包括数据采集设备和监控计算机的模型在线更新***；所述监控计算机内预置有一传统PCA建模模块，一均值方差更新模块，一投影点计算模块，一残差判定模块，一主元空间调整模块，一主元方向微调模块和一控制限更新模块；2)在模型在线更新***启动前，数据采集设备收集能代表生产过程特征的历史生产数据作为训练样本x输入给监控计算机中的传统PCA建模模块，建立PCA初始监控模型Ω，并将PCA初始监控模型传送给均值方差更新模块；PCA初始监控模型Ω为：

$\mathrm{\Ω}=(\stackrel{\‾}{x},\mathrm{\σ},P_{\mathrm{nk}},{\mathrm{\Λ}}_{\mathrm{kk}},N,{\mathrm{\δ}}_{\mathrm{\α}}^2,T_{\mathrm{\α}}^2)$

式中， $\stackrel{\‾}{x}={\mathrm{\Σ}}_{i=1}^N{x}_i/N$ 表示模型均值； $\mathrm{\σ}=\sqrt{{\mathrm{\Σ}}_{i=1}^N{(x_i-\stackrel{\‾}{x})}^2/(N-1)}$ 表示标准差；P_nk表示特征值分解后最大的k个主元得分所对应的负荷向量，即主元投影向量；Λ_kk表示特征值矩阵，其对角元素由最大的k个主元得分构成；N表示选取的训练样本个数；

表示SPE统计量在显著性水平α下的控制限；表示T²统计量在显著性水平α下的控制限；x_i表示第i个训练样本(i＝1，…，N)；k表示保留的主元个数；n表示过程变量的维数；3)数据采集设备将工业生产过程中的实时过程变量转换为实时过程数据x_N+m，传送给均值方差更新模块；同时，模型在线更新***将所需的更新参数也传送给均值方差更新模块；而均值方差更新模块中的现有PCA监控模型为

均值方差更新模块根据实时过程数据x_N+m和更新参数，计算出新的PCA监控模型

中的模型均值

和标准差σ′，将计算结果进行存储，并将现有PCA监控模型以及相关参数传送给投影点计算模块；4)投影点计算模块计算出实时过程数据x_N+m在主元空间上的投影点g的坐标，以及实时过程数据x_N+m到投影点g的距离，即残差向量h，并将计算结果、现有PCA监控模型以及相关参数传送给残差判定模块，用于在后续步骤中对主元投影方向进行调整；具体计算公式如下：

$g=P_{\mathrm{nk}}^T(x_{N+m}-{\stackrel{\‾}{x}}^{\′}){\mathrm{\Σ}}_{N+m}^{-1}$

$h=(x_{N+m}-{\stackrel{\‾}{x}}^{\′}){\mathrm{\Σ}}_{N+m}^{-1}-P_{\mathrm{nk}}g$

其中，∑_N+m＝diag(σ′)，是以更新后模型的标准差向量σ′为对角元素的对角阵；5)残差判定模块根据投影点的残差向量的模的大小决定调用哪一种模型投影方向更新方法：如果残差值较大，则调用主元空间调整模块，进入步骤6)；如果残差值较小，则调用主元方向微调模块，进入步骤7)；主元空间调整模块和主元方向微调模块的输出形式相同，均为更新后模型的主元投影向量以及特征值矩阵；残差判定模块通过比较残差向量h的模与增量阈值η的大小来选择当前样本的更新方法：①当残差h的模较大时，即||h||＞η，将优化的残差向量现有PCA模型和相关参数传送给主元空间调整模块，进入步骤6)，对整个主元空间进行调整计算；②当残差h的模较小时，即||h||＜η，将现有PCA模型和相关参数传送给主元方向微调模块，进入步骤7)，进行主元投影方向的微调计算；6)主元空间调整模块基于过程变量相关性本质不变的思想，利用主元分析方法对扩展的主元投影向量进行降维，从而在保持主元个数不变的同时，对现有主元投影方向进行更新，得到新模型的负荷向量P′_nk和特征值矩阵Λ′_kk；特征值矩阵Λ′_kk的求解公式为：

$(\mathrm{\μ}\left[\begin{array}{cc}{\mathrm{\Λ}}_{\mathrm{kk}}& \underset{\‾}{0}\\ {\underset{\‾}{0}}^T& 0\end{array}\right]+(1-\mathrm{\μ})\left[\begin{array}{cc}{\mathrm{qq}}^T& \mathrm{\ρq}\\ {\mathrm{\ρq}}^T& {\mathrm{\ρ}}^2\end{array}\right])\≈{\mathrm{R\Λ}}_{\mathrm{kk}}^{\′}{R}^T$

式中， $\mathrm{\ρ}={\hat{h}}^T\left[{\mathrm{\Σ}}_{N+m}^{-1}(x_{N+m}-\stackrel{\‾}{x})\right],$ $q=P_{\mathrm{nk}}^T\left[{\mathrm{\Σ}}_{N+m}^{-1}(x_{N+m}-\stackrel{\‾}{x})\right],$ R表示旋转矩阵，对左式中的矩阵进行特征值分解，最大的k个特征值构成了更新后的特征值矩阵Λ′_kk，其对应的特征值向量就是旋转矩阵R；由于现有的负荷向量P_nk中的每个列向量与残差向量h彼此正交，构成扩展的主元投影向量，则新模型的负荷向量P′_nk是扩展的主元投影向量的线性组合；将上式中的旋转矩阵R代入，得到新模型的负荷向量P′_nk的计算公式为：

$P_{\mathrm{nk}}^{\′}=[P_{\mathrm{nk}},\hat{h}]R$

该模块将计算得到的新模型的主元投影向量P′_nk进行存储，并将特征值矩阵Λ′_kk传送到控制限更新模块，进入步骤8)，继续完成对主元得分的调整以及对控制限的更新工作；7)主元方向微调模块同样利用现有PCA模型和之前步骤得到的结果对新模型的负荷向量P′_nk和特征值矩阵Λ′_kk进行计算；由于只是在主元空间内部对主元投影方向进行微调，仅需要在主元空间内部重新进行主元分析，即得到更新后的主元投影方向；特征值矩阵Λ′_kk的求解公式为：

(λΛ_kk+(1-λ)qq^T)≈RΛ′_kkR^T

对左式中的矩阵进行特征值分解，所有k个特征值构成了更新后的特征值矩阵Λ′_kk，其对应的特征值向量就是旋转矩阵R；新模型的主元投影向量P′_nk是由现有的主元投影向量P_nk中的列向量所构成的一个新的线性组合，由以下公式计算得出：

P′_nk＝P_nkR

该模块将计算得到的新模型的主元投影向量P′_nk进行存储，并将特征值矩阵Λ′_kk传送到控制限更新模块，进入步骤8)，继续完成对主元得分的调整以及对控制限的更新工作；8)控制限更新模块基于前面的计算结果对更新后的特征值矩阵Λ′_kk进行规范化，从而完成对模型监控统计量指标的控制限

和

的更新，最终得到完全更新后的新模型；控制限更新模块将上述计算结果进行存储，将最终得到的完全更新后的新模型

进行输出，用于工业过程中的在线监控以及故障诊断；对更新后的特征值矩阵Λ′_kk进行规范化的方法为：定义现有模型中前k个主元的特征值累积贡献率为Cum(Λ_kk)，即前k个保留的特征值的累加和占所有特征值累加和的比例为Cum(Λ_kk)；由于新模型仅对前k个主元的特征值进行了调整，因而记新模型的前k个主元的特征值累积贡献率为Cum(Λ′_kk)，则新的特征值矩阵Λ′_kk的每个对角线上元素应调整为：

Λ′_i，i＝Λ′_i，i×Cum(Λ)/Cum(Λ′)

式中，Λ′_i，i表示新的特征值矩阵Λ′_kk的第(i，i)个元素；经过以上调整，在物理关联性不变的前提下，新模型的SPE控制限

无需进行更新调整也能满足绝大多数情况下的监控要求；因此，新模型的SPE控制限

所述步骤3)中，新的PCA监控模型的模型均值和标准差σ′的更新公式如下：

${\stackrel{\‾}{x}}^{\′}=\mathrm{\λ}\stackrel{\‾}{x}+(1-\mathrm{\λ})x_{N+m}$

${\mathrm{\σ}}^{\′}=\sqrt{{\mathrm{\λ\σ}}^2+(1-\mathrm{\λ}){(x_{N+m}-\stackrel{\‾}{x})}^2}$

式中，x_N+m表示在实际监控过程中采集到的第m个实时过程变量的实时过程数据；m＝1，2，3…；定义监控过程中，均值方差更新模块中的现有PCA监控模型为当监测第一个实时过程变量时，即m＝1时，均值方差更新模块中的现有PCA监控模型就是传统PCA建模模块给出的初始监控模型

λ为遗忘因子，属于更新参数，其取值范围为(0，1)；增量阈值η也属于更新参数，取值范围为(0，δ_α)。

所述步骤8)中，新的PCA监控模型的T²控制限

采用F分布进行计算，公式为：

$T_a^{2\′}=\frac{k({(N+m)}^2-1)}{(N+m)(N+m-k)}{F}_{k,N+m-k;\mathrm{\α}}$

其中，k为保留的主元个数，N为初始建模模块收集的训练样本总数，m为当前在监控阶段采集到的样本总数，F_{k，N+m-k；α}对应于检验水平为α，自由度为k，N+m-k条件下的F分布临界值，通过查表得到，因此，仅需要及时更新训练样本数N的数值，即得到新的PCA监控模型的T²控制限

本发明由于采取以上技术方案，其具有以下优点：1、本发明通过对模型的投影方向进行更新，从而保证了模型结构与实际***的一致性，与仅更新模型均值、标准差的PCA模型更新方法相比，具有更低的虚警率和更加准确的监控效果。2、本发明基于过程变量相关性本质不变的思想，控制主元个数保持不变，因而无需使用整个数据协方差矩阵进行更新，仅需要对反映***性质的少数主元投影方向进行微调；因此在实际应用中，本发明与其他需要更新整个模型结构的PCA模型更新方法相比，具有更快的计算速度和更低的存储需求，并且能够得到相似甚至更好的监控效果。3、本发明对工业过程中经常发生的数据缓慢漂移现象具有鲁棒性，能够对真实故障所引起变量突变及变量相关性消失等现象及时进行报警。本发明构思巧妙，精确实用，可广泛用于实际的工业过程监控过程中。

附图说明

图1是本发明结构示意图

图2是本发明模块结构示意图

图3是使用传统PCA监控模型的实施例监控结果示意图

图4是本发明实施例结果示意图

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明基于以下思想：在正常运行的工业***中，过程变量间的物理关联性不会随着时间发生变化，因而，变量间的相关性本质不会发生改变，主元个数应当维持不变。本发明通过实时更新模型的均值和标准差，并对模型的主元投影方向进行调整，仅仅利用PCA模型的基本主元信息，就能够对监控***进行快速而有效的更新，从而在实施工业***在线监控的同时，实现了根据实时过程数据对PCA模型进行在线更新的功能。

本发明方法可以概括为：如图1所示，根据工业现场实际情况，设置一包括数据采集设备1、监控计算机2的模型在线更新***。数据采集设备1将工业生产过程中的实时过程变量和历史过程变量转换为数字信号，并传送给监控计算机2；监控计算机2通过处理历史过程数据建立初始PCA模型，并根据实时过程数据以及更新参数进行监控模型的在线更新，实时输出更新后的PCA监控模型。

本发明方法具体包括以下步骤：

1)如图1所示，根据工业现场实际情况，设置一包括数据采集设备1、监控计算机2的模型在线更新***。如图2所示，监控计算机2内设置有一传统PCA建模模块21，一均值方差更新模块22，一投影点计算模块23，一残差判定模块24，一主元空间调整模块25，一主元方向微调模块26和一控制限更新模块27。

2)在模型在线更新***启动前，首先将数据采集设备1收集的历史数据进行人工筛选，根据一定的工业现场经验挑选出若干充分且具有代表性的历史数据作为训练样本传送给传统PCA建模模块21。传统PCA建模模块21仅在模型在线更新***启动前被调用一次，根据训练样本建立PCA初始监控模型，并将PCA初始监控模型传送给均值方差更新模块22。当在线更新开始后，传统PCA建模模块21将处于休眠状态，不再参与在线更新工作。第一次更新计算完成后，新的监控模型将由***自身进行存储，不再需要外界进行输入。

传统PCA建模模块21采用传统的主元分析建模方法建立PCA模型：对标准化后的协方差矩阵进行特征值分解，以得到模型的负荷向量和主元得分。通过累积贡献率阈值法，或交叉验证法确定保留的主元个数，并根据统计学公式计算出相应的统计量控制限，最终得到完整的PCA初始监控模型Ω：

$\mathrm{\Ω}=(\stackrel{\‾}{x},\mathrm{\σ},P_{\mathrm{nk}},{\mathrm{\Λ}}_{\mathrm{kk}},N,{\mathrm{\δ}}_{\mathrm{\α}}^2,T_{\mathrm{\α}}^2)$

式中， $\stackrel{\‾}{x}={\mathrm{\Σ}}_{i=1}^N{x}_i/N$ 表示模型均值； $\mathrm{\σ}=\sqrt{{\mathrm{\Σ}}_{i=1}^N{(x_i-\stackrel{\‾}{x})}^2/(N-1)}$ 表示标准差；P_nk表示特征值分解后最大的k个主元得分所对应的负荷向量，即主元投影向量；Λ_kk表示特征值矩阵，其对角元素由最大的k个主元得分构成；N表示选取的训练样本个数；

表示SPE(Squared Prediction Error，平方预测误差)统计量在显著性水平α下的控制限；

表示T²统计量(描述主元子空间上的样本变化情况)在显著性水平α下的控制限；x_i表示第i个训练样本(i＝1，…，N)；k表示保留的主元个数；n表示过程变量的维数。

3)数据采集设备1将工业生产过程中的实时过程变量转换为实时过程数据，传送给监控计算机2中的均值方差更新模块22。同时，模型在线更新***将所需的更新参数也传送给均值方差更新模块22。更新参数包括：遗忘因子λ，取值范围为(0，1)；增量阈值η，取值范围(0，δ_α)，其中，δ_α为PCA初始监控模型Ω中SPE控制限的根号值。上述更新参数的初值一般由人为设定，并且可以在实际监控过程中根据具体情况进行人为或机器自动的调整。

假设在实际监控过程中采集到的第m个样本的实时过程数据为x_N+m，m＝1，2，3…，并以该样本的监控过程为例说明PCA监控模型的更新流程。定义监控过程中的现有PCA监控模型为

那么，当监测第一个样本的实时过程数据x_N+1时，现有PCA监控模型就是传统PCA建模模块21给出的初始监控模型 $\mathrm{\Ω}=(\stackrel{\‾}{x},\mathrm{\σ},P_{\mathrm{nk}},{\mathrm{\Λ}}_{\mathrm{kk}},N,{\mathrm{\δ}}_{\mathrm{\α}}^2,T_{\mathrm{\α}}^2).$

均值方差更新模块22根据所采集到的新样本的实时过程数据x_N+m，以及所需的更新参数，计算出新的PCA监控模型

中的模型均值

和标准差σ′，对计算结果进行存储，并将现有的PCA监控模型Ω以及相关参数传送给投影点计算模块23。其中，新的PCA监控模型的模型均值

和标准差σ′的更新公式如下：

${\stackrel{\‾}{x}}^{\′}=\mathrm{\λ}\stackrel{\‾}{x}+(1-\mathrm{\λ})x_{N+m}$

${\mathrm{\σ}}^{\′}=\sqrt{{\mathrm{\λ\σ}}^2+(1-\mathrm{\λ}){(x_{N+m}-\stackrel{\‾}{x})}^2}$

对模型均值和标准差进行实时的归一化可以确保PCA算法的持续有效性，并减小量纲对模型准确度的影响。

4)投影点计算模块23主要用于计算新样本x_N+m在主元空间上的投影点g的坐标，以及残差向量h(用来描述x_N+m到投影点g的距离)，并将计算结果、现有PCA监控模型Ω以及相关参数传送给残差判定模块24，用于在后续步骤中对主元投影方向进行调整。具体计算公式如下：

$g=P_{\mathrm{nk}}^T(x_{N+m}-{\stackrel{\‾}{x}}^{\′}){\mathrm{\Σ}}_{N+m}^{-1}$

$h=(x_{N+m}-{\stackrel{\‾}{x}}^{\′}){\mathrm{\Σ}}_{N+m}^{-1}-P_{\mathrm{nk}}g$

其中，∑_N+m＝diag(σ′)，是以更新后模型的标准差向量σ′为对角元素的对角阵。

5)残差判定模块24根据残差向量的模||h||的大小决定调用哪一种模型投影方向更新方法：如果残差值较大，则调用主元空间调整模块25，进入步骤6)；如果残差值较小，则调用主元方向微调模块26，进入步骤7)；其中，主元空间调整模块25和主元方向微调模块26的输出形式相同，均为更新后模型的主元投影向量以及特征值矩阵。根据残差向量决定调节方法可以避免由噪声所引起的过度调节，增强***对噪声的鲁棒性。

残差判定模块24通过比较残差向量h的模与增量阈值η的大小来选择当前样本的更新方法：

①当残差h的模较大时，即||h||＞η，意味着新样本点到原有的主元投影空间的距离较远，需要对整个主元空间的投影方向进行调整，因此将优化的残差向量

现有PCA模型和相关参数传送给主元空间调整模块25，进入步骤6)，对整个主元空间的投影方向进行调整计算；

②当残差h的模较小时，即||h||＜η，新的样本基本符合现有模型的变量相关性假设，只需要在原主元空间内部对现有的主元方向进行微调即可，因此将原数据、现有PCA模型和相关参数传送给主元方向微调模块26，进入步骤7)，进行主元投影方向的微调计算。

6)主元空间调整模块25基于过程变量相关性本质不变的思想，利用主元分析方法对扩展的主元投影向量进行降维，从而在保持主元个数保持不变的同时，对现有主元投影方向进行更新，得到新模型的负荷向量P′_nk和特征值矩阵Λ′_kk。

特征值矩阵Λ′_kk的求解公式为：

$(\mathrm{\μ}\left[\begin{array}{cc}{\mathrm{\Λ}}_{\mathrm{kk}}& \underset{\‾}{0}\\ {\underset{\‾}{0}}^T& 0\end{array}\right]+(1-\mathrm{\μ})\left[\begin{array}{cc}{\mathrm{qq}}^T& \mathrm{\ρq}\\ {\mathrm{\ρq}}^T& {\mathrm{\ρ}}^2\end{array}\right])\≈{\mathrm{R\Λ}}_{\mathrm{kk}}^{\′}{R}^T---\left(1\right)$

式中， $\mathrm{\ρ}={\hat{h}}^T\left[{\mathrm{\Σ}}_{N+m}^{-1}(x_{N+m}-\stackrel{\‾}{x})\right],$ $q=P_{\mathrm{nk}}^T\left[{\mathrm{\Σ}}_{N+m}^{-1}(x_{N+m}-\stackrel{\‾}{x})\right],$ (k+1)×k维矩阵R表示待求的旋转矩阵。对式(1)左式中的矩阵进行特征值分解，挑选最大的k个特征值构成更新后的特征值矩阵Λ′_kk，其对应的特征向量就是旋转矩阵R。

由于现有的负荷向量P_nk中的每个列向量与残差向量h彼此正交，可以构成扩展的主元投影向量，则新模型的负荷向量P′_nk是扩展的主元投影向量的线性组合。将上式中的旋转矩阵R代入，可以得到新模型的负荷向量P′_nk的计算公式为：

$P_{\mathrm{nk}}^{\′}=[P_{\mathrm{nk}},\hat{h}]R$

该模块将计算得到的新模型的负荷向量P′_nk进行存储，并将特征值矩阵Λ′_kk、现有PCA模型和相关参数传送到控制限更新模块27，进入步骤8)，继续完成对主元得分的调整以及对控制限的更新工作。

7)主元方向微调模块26同样利用现有PCA模型和之前步骤得到的结果对新模型的负荷向量P′_nk和特征值矩阵Λ′_kk进行计算。由于只是在主元空间内部对主元投影方向进行微调，仅需要在主元空间内部重新进行主元分析，就可以得到更新后的主元投影方向。

特征值矩阵Λ′_kk的求解公式为：

(λΛ_kk+(1-λ)qq^T)≈RΛ′_kkR^T    (2)

对式(2)左式中的矩阵进行特征值分解，所有k个特征值构成了更新后的特征值矩阵Λ′_kk，其对应的特征向量就是旋转矩阵R(R是k×k维矩阵)。新模型的负荷向量P′_nk应当是由现有的负荷向量P_nk中的列向量所构成的一个新的线性组合，可由以下公式计算得出：

P′_nk＝P_nkR

该模块将计算得到的新模型的主元投影向量P′_nk进行存储，并将特征值矩阵Λ′_kk、现有PCA模型和相关参数传送到控制限更新模块27，进入步骤8)，继续完成对主元得分的调整以及对控制限的更新工作。

8)控制限更新模块27基于前面的计算结果对更新后的特征值矩阵Λ′_kk进行规范化，从而完成对模型监控统计量指标的控制限和的更新，最终得到完全更新后的新模型。

对更新后的特征值矩阵Λ′_kk进行规范化的方法为：定义现有模型中前k个主元的特征值累积贡献率为Cum(Λ_kk)，即前k个保留的特征值的累加和占所有特征值累加和的比例为Cum(Λ_kk)。由于新模型仅对前k个主元的特征值进行了调整，因而记新模型的前k个主元的特征值累积贡献率为Cum(Λ′_kk)，则新的特征值矩阵Λ′_kk的每个对角线上元素应调整为：

Λ′_i，i＝Λ′_i，i×Cum(Λ)/Cum(Λ′)

式中，Λ′_i，i表示新的特征值矩阵Λ′_kk的第(i，i)个元素。

经过以上调整，主元空间和残差空间的比例将基本保持不变，在物理关联性不变的前提下，新模型的SPE控制限

无需进行更新调整也能满足绝大多数情况下的监控要求。因此，新模型的SPE控制限

另一方面，新的PCA监控模型的T²控制限

可以采用传统算法计算得到，例如可以利用F分布进行计算，公式为：

$T_a^{2\′}=\frac{k({(N+m)}^2-1)}{(N+m)(N+m-k)}{F}_{k,N+m-k;\mathrm{\α}}$

其中，k为保留的主元个数，N为初始建模模块收集的训练样本总数，m为当前在监控阶段采集到的样本总数，F_{k，N+m-k；α}是对应于显著性水平为α，自由度为k，N+m-k条件下的F分布临界值，可查表得到。显然，本发明仅需要及时更新监控阶段采集到的样本数以得到新的PCA监控模型的T²控制限

该模块将计算得到的新特征值矩阵Λ′_kk以及T²控制限

进行存储，完成了新模型的全部更新工作，将完整的新模型

进行输出，用于工业过程中的在线监控以及故障诊断。

下面列举一具体实施例对本发明的应用进行详细描述。

本实例采用北方微电子公司等离子体刻蚀***的一组实际生产数据作为实施样例。该实验数据包含339个样本，样本维数为16维，涵盖了反映***正常变化特性的16个非设定值的工程变量，所有样本均为正常样本，不存在故障样本。由于设备积灰等原因，实验数据具有较为明显的数据漂移现象。采用传统PCA监控模型在该数据集上的监控结果如图3所示，建模样本为前100组样本。图中横坐标为样本序号(1～339)，严格按照样本(wafer)的加工顺序排列。纵坐标代表各统计量的监控数值，实线(图中蓝色线)是各统计量的控制限，圆圈(图中红色圆圈)代表每个样本的实际统计量监控值。在一般情况下，如果某个样本的圆圈在任意一个控制图中位于控制限以上(SPE或T²统计量超出控制限)，则***应发出报警。观察可知，在整个监控过程中，监控误差不断增大，SPE统计量逐渐超出控制限，虚警率高达64.02％，监控模型渐趋失效。

使用上述实际数据对本发明中的方法进行仿真，具体的仿真步骤和实验结果如下：

(1)设定SPE和T²统计量控制限的显著性水平为95％，指定更新参数中的遗忘

因子为λ＝0.95，增量阈值为η＝0.5。

(2)调用传统PCA建模模块21对前100组样本进行初始建模，依次读取所有剩下的正常样本作为实时过程数据进行模型在线更新。将模型在线更新***输出的监控模型用于进行在线监控。

(3)给出本发明中方法对所有样本的监控结果，包括SPE控制图和T²控制图，

如图4所示(图中标记含义与图3中完全相同)。在整个监控阶段，本发明中方法的虚警率仅为27.62％，监控效果较好。

由实施例可知，采用本发明进行模型在线更新得到的监控模型，具有较低的虚警率，对大多数工业过程中容易出现的过程变量缓慢漂移现象具有鲁棒性。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (2)

1.一种主元分析监控模型的在线更新方法，其包括以下步骤：

1）在工业现场设置一包括数据采集设备和监控计算机的模型在线更新***；所述监控计算机内预置有一传统主元分析建模模块，一均值方差更新模块，一投影点计算模块，一残差判定模块，一主元空间调整模块，一主元方向微调模块和一控制限更新模块；

2）在模型在线更新***启动前，数据采集设备收集能代表生产过程特征的历史生产数据作为训练样本x输入给监控计算机中的传统主元分析建模模块，建立主元分析初始监控模型Ω，并将主元分析初始监控模型传送给均值方差更新模块；

主元分析初始监控模型Ω为：

$\mathrm{\Ω}=(\stackrel{\‾}{x},\mathrm{\σ},P_{\mathrm{nk}},{\mathrm{\Λ}}_{\mathrm{kk}},N,{\mathrm{\δ}}_{\mathrm{\α}}^2,T_{\mathrm{\α}}^2)$

式中，表示模型均值；

表示标准差；P_nk表示特征值分解后最大的k个主元得分所对应的负荷向量，即主元投影向量；Λ_kk表示特征值矩阵，其对角元素由最大的k个主元得分构成；N表示选取的训练样本个数；

表示平方预测误差统计量在显著性水平α下的控制限；

表示T²统计量在显著性水平α下的控制限；x_i表示第i个训练样本，i＝1,…,N；k表示保留的主元个数；n表示过程变量的维数；

3）数据采集设备将工业生产过程中的实时过程变量转换为实时过程数据x_N+m，传送给均值方差更新模块；同时，模型在线更新***将所需的更新参数也传送给均值方差更新模块；而均值方差更新模块中的现有主元分析监控模型为

均值方差更新模块根据实时过程数据x_N+m和更新参数，计算出新的主元分析监控模型

中的模型均值

和标准差σ′，将计算结果进行存储，并将现有主元分析监控模型以及相关参数传送给投影点计算模块；其中，新的主元分析监控模型的模型均值

和标准差σ′的更新公式如下：

${\stackrel{\‾}{x}}^{\′}=\mathrm{\λ}\stackrel{\‾}{x}+(1-\mathrm{\λ})x_{N+m}$

${\mathrm{\σ}}^{\′}=\sqrt{{\mathrm{\λ\σ}}^2+(1-\mathrm{\λ}){(x_{N+m}-\stackrel{\‾}{x})}^2}$

式中，x_N+m表示在实际监控过程中采集到的第m个实时过程变量的实时过程数据；m＝1,2,3…；定义监控过程中，均值方差更新模块中的现有主元分析监控模型为

当监测第一个实时过程变量时，即m＝1时，均值方差更新模块中的现有主元分析监控模型就是传统主元分析建模模块给出的初始监控模型

λ为遗忘因子，属于更新参数，其取值范围为（0,1）；增量阈值η也属于更新参数，取值范围为（0,δ_α）；

4）投影点计算模块计算出实时过程数据x_N+m在主元空间上的投影点g的坐标，以及实时过程数据x_N+m到投影点g的距离，即残差向量h，并将计算结果、现有主元分析监控模型以及相关参数传送给残差判定模块，用于在后续步骤中对主元投影方向进行调整；

具体计算公式如下：

$g=P_{\mathrm{nk}}^T(x_{N+m}-{\stackrel{\‾}{x}}^{\′}){\mathrm{\Σ}}_{N+m}^{-1}$

$h=(x_{N+m}-{\stackrel{\‾}{x}}^{\′}){\mathrm{\Σ}}_{N+m}^{-1}-P_{\mathrm{nk}}g$

其中，Σ_N+m＝diag(σ′)，是以更新后模型的标准差向量σ′为对角元素的对角阵；

5）残差判定模块根据投影点的残差向量的模的大小决定调用哪一种模型投影方向更新方法：如果残差值较大，则调用主元空间调整模块，进入步骤6）；如果残差值较小，则调用主元方向微调模块，进入步骤7）；主元空间调整模块和主元方向微调模块的输出形式相同，均为更新后模型的主元投影向量以及特征值矩阵；

残差判定模块通过比较残差向量h的模与增量阈值η的大小来选择当前样本的更新方法：

①当残差h的模较大时，即||h||＞η，将优化的残差向量现有主元分析模型和相关参数传送给主元空间调整模块，进入步骤6），对整个主元空间进行调整计算；

②当残差h的模较小时，即||h||＜η，将原数据、现有主元分析模型和相关参数传送给主元方向微调模块，进入步骤7），进行主元投影方向的微调计算；

6）主元空间调整模块基于过程变量相关性本质不变的思想，利用主元分析方法对扩展的主元投影向量进行降维，从而在保持主元个数不变的同时，对现有主元投影方向进行更新，得到新模型的负荷向量

和特征值矩阵

特征值矩阵

的求解公式为：

$(\mathrm{\λ}\left[\begin{array}{cc}{\mathrm{\Λ}}_{\mathrm{kk}}& \underset{\‾}{0}\\ {\underset{\‾}{0}}^T& 0\end{array}\right]+(1-\mathrm{\λ})\left[\begin{array}{cc}{\mathrm{qq}}^T& \mathrm{\ρq}\\ {\mathrm{\ρq}}^T& {\mathrm{\ρ}}^2\end{array}\right])\≈{\mathrm{R\Λ}}_{\mathrm{kk}}^{\′}{R}^T$

式中， $\mathrm{\ρ}={\hat{h}}^T\left[{\mathrm{\Σ}}_{N+m}^{-1}(x_{N+m}-\stackrel{\‾}{x})\right],q=P_{\mathrm{nk}}^T\left[{\mathrm{\Σ}}_{N+m}^{-1}(x_{N+m}-\stackrel{\‾}{x})\right],$ R表示旋转矩阵，对左式中的矩阵进行特征值分解，最大的k个特征值构成了更新后的特征值矩阵其对应的特征值向量就是旋转矩阵R；

由于现有的负荷向量P_nk中的每个列向量与残差向量h彼此正交，构成扩展的主元投影向量，则新模型的负荷向量

是扩展的主元投影向量的线性组合；将上式中的旋转矩阵R代入，得到新模型的负荷向量

的计算公式为：

$P_{\mathrm{nk}}^{\′}=[P_{\mathrm{nk}},\hat{h}]R$

该模块将计算得到的新模型的主元投影向量

进行存储，并将特征值矩阵

以及相关参数传送到控制限更新模块，进入步骤8），继续完成对主元得分的调整以及对控制限的更新工作；

7）主元方向微调模块同样利用现有主元分析模型和之前步骤得到的结果对新模型的负荷向量

和特征值矩阵

进行计算；由于只是在主元空间内部对主元投影方向进行微调，仅需要在主元空间内部重新进行主元分析，即得到更新后的主元投影方向；

特征值矩阵

的求解公式为：

$({\mathrm{\λ\Λ}}_{\mathrm{kk}}+{(1-\mathrm{\λ})\mathrm{qq}}^T)\≈{\mathrm{R\Λ}}_{\mathrm{kk}}^{\′}{R}^T$

对左式中的矩阵进行特征值分解，所有k个特征值构成了更新后的特征值矩阵

，其对应的特征值向量就是旋转矩阵R；新模型的主元投影向量是由现有的主元投影向量P_nk中的列向量所构成的一个新的线性组合，由以下公式计算得出：

$P_{\mathrm{nk}}^{\′}=P_{\mathrm{nk}}R$

该模块将计算得到的新模型的主元投影向量

进行存储，并将特征值矩阵

以及相关参数传送到控制限更新模块，进入步骤8），继续完成对主元得分的调整以及对控制限的更新工作；

8）控制限更新模块基于前面的计算结果对更新后的特征值矩阵

进行规范化，从而完成对模型监控统计量指标的控制限

和

的更新，最终得到完全更新后的新模型；控制限更新模块将上述计算结果进行存储，将最终得到的完全更新后的新模型

进行输出，用于工业过程中的在线监控以及故障诊断；

对更新后的特征值矩阵

进行规范化的方法为：定义现有模型中前k个主元的特征值累积贡献率为Cum(Λ_kk)，即前k个保留的特征值的累加和占所有特征值累加和的比例为Cum(Λ_kk)；由于新模型仅对前k个主元的特征值进行了调整，因而记新模型的前k个主元的特征值累积贡献率为

则新的特征值矩阵的每个对角线上元素应调整为：

${\mathrm{\Λ}}_{i,i}^{\′}={\mathrm{\Λ}}_{i,i}^{\′}\×\mathrm{Cum}\left(\mathrm{\Λ}\right)/\mathrm{Cum}\left({\mathrm{\Λ}}^{\′}\right)$

式中，表示新的特征值矩阵

的第(i,i)个元素；

经过以上调整，在物理关联性不变的前提下，新模型的平方预测误差控制限

无需进行更新调整也能满足绝大多数情况下的监控要求；因此，新模型的平方预测误差控制限

2.如权利要求1所述的一种主元分析监控模型的在线更新方法，其特征在于：所述步骤8）中，新的主元分析监控模型的T²控制限采用F分布进行计算，公式为：

$T_a^{2\′}=\frac{k({(N+m)}^2-1)}{(N+m)(N+m-k)}{F}_{k,N+m-k;\mathrm{\α}}$

其中，k为保留的主元个数，N为初始建模模块收集的训练样本总数，m为当前在监控阶段采集到的样本总数，F_k,N+m-k;α对应于检验水平为α，自由度为k，N+m-k条件下的F分布临界值，通过查表得到，因此，仅需要及时更新训练样本数N的数值，即得到新的主元分析监控模型的T²控制限

Images