为了提高现代工业过程的生产效率，近年来故障检测越来越重要.故障检测有三种主要方法：基于知识的方法、基于模型的方法、基于数据驱动的方法^[1].基于知识的方法的缺点是需要大量不容易获取的过程知识.尽管数据驱动方法已被广泛应用，但在过去的二十年中基于子空间辨识模型的方法吸引了过程建模和监控领域的学者和专家的关注.常规的子空间辨识方法包括规范变量分析(CVA)，子空间状态空间系统辨识的数值算法(N4SID)，多变量输出误差状态空间(MOESP)^[2]. CVA对线性过程静态操作条件表现出了非常好的性能，但是，许多工业过程需要进行预定的正常改变，比如设定值改变、化学反应器进料比的改变等.因此，针对时变过程开发在线故障检测方法是非常必要的.

递推辨识方法存在的关键问题是较高的计算负荷.根据获取状态变量和可观测矩阵的不同，递推子空间辨识方法可以分为两类：一类是避免使用奇异值分解.Gustafuson^[3]首次提出了基于投影近似子空间跟踪(PAST)的递推子空间辨识方法.Oku等^[4]提出了使用梯度类型子空间跟踪的递推4SID算法.Houtzager等^[5]提出了基于predictor的递推辨识方法.另一类是直接使用奇异值分解或降低奇异值的计算负荷.Choi等^[6]提出了基于CVA更新状态空间模型的自适应故障检测方法.Ding等^[7]提出了基于等价空间和观测器的故障检测和隔离系统.

为了降低Hankel矩阵奇异值分解的计算负荷，引入一阶干扰理论实现Hankel矩阵的递推奇异值分解更新.与常规奇异值分解相比，基于一阶干扰理论的递推奇异值分解显著降低了计算负荷.本文提出了基于一阶干扰理论的递推规范变量分析(RCVA-FOP)时变过程故障检测方法，并将该方法用于田纳西伊斯曼化工过程的在线故障检测.

1 规范变量分析

多变量状态空间模型为

其中：x_(k)∈Rⁿ和y_(k)∈R^m_y是采样时刻k的状态和观测矢量；A∈R^n×n，B∈R^n×n，C∈R^m_y×n，D∈R^n×n是系统矩阵；w_(k)和v_(k)是随机干扰，假定是零均值正态分布白噪声.

在CVA方法中，测量矢量由过去和未来的测量值扩展而成，分别形成过去、未来矢量z_{p (r)}，y_{f (r)}，其中r表示一类序数.

其中，p代表过去和未来观测值的长度.一般情况，n/m_y≤p.

设置r=p+1, p+2, …, p+N，过去和未来Hankel矩阵Y_p∈R^mp×N和Y_f∈R^mp×N定义如下：

一组由l个观测值组成的变量，y_{p (p+1)}最后的两个元素是y₍₁₎，u₍₁₎，然而y_{f (p+N)}最后的元素是y_(l).因此，这些Hankel矩阵的最大列数是N=l-p-f+1.

过去和未来观测值的协方差和互方差矩阵可以使用如下公式计算：

CVA试图找到a^Ty_{f (r)}和b^Ty_{p (r)}之间的最优线性组合，以便于过去和未来观测值之间的相关关系最大化.该最优化问题的解决方案可由Hankel矩阵H的奇异值分解求得.

对上式进行SVD，规范变量z_(r)可基于过去测量值求得：

其中J=V^TΣ_pp^-1/2是转换过去测量值到规范变量空间的转换矩阵.规范变量空间可以被分成两个部分：主空间和残差空间，它们的维数由奇异值的维数确定.前n个主要奇异值确定主空间x_(r)∈Rⁿ，并且剩余的奇异值确定残差空间r_(r)∈R^mp-n，其中n < mp，奇异值个数n应事先给出.因此，规范变量可以写成：

其中，J_x=V_x^TΣ_pp^-1/2，V_x包含V的前n列.状态变量不仅是估计的规范变量的子集，也是过去观测矢量的线性组合.因此，像规范变量一样状态变量定义为过去观测矢量的线性组合.在系统状态被确定之后，矩阵A, B, C, D可以通过线性最小二乘递推来估计.最后，可以重新计算系统的输出.

2 基于递推规范变量分析的故障检测

对多模过程，每个模式有不同的特性，例如变量之间均值、方差和相关关系结构的改变导致主要状态变量个数的改变.每次采集到新的测量变量时，均值和协方差随着模型结构改变程度进行更新.本文中，采用指数权重滑动平均而不需要回调过去的训练数据.

过去观测矢量的协方差矩阵可以通过指数权重滑动平均来估计.

其中β是遗忘因子.

2.1 基于一阶干扰理论更新状态矢量

基于干扰理论更新奇异值和奇异值矢量的根源可以参考文献[8]中提到的结果.在例子中，奇异值和奇异值矢量的初始集来自离线辨识.在线更新状态矢量开始于当前测量值.z_k表示为z_k=z_{f (r)}z_{f (r)}^Tz_{f (r)}，其中，z_{f (r)}是列矢量包含当前和未来测量值.相似地，增补变量矢量可以利用以前的测量值z_τ=z_{p (r)}z_{p (r)}^Tz_{p (r)}.

定义如下过渡变量，和p_{u, 0}=，q_{u, 0}=[0 … 0].

首先，更新左奇异值矢量执行如下操作：

下标1≤i≤q代表矩阵列或矢量的第i元素.

然后，更新右奇异值矢量，定义两个过渡变量:p_{v, 0}=，q_{v, 0}=[0 … 0].

右奇异值矢量可计算如下：

奇异值矢量更新如下：

最后标准化后如下：

给定更新的奇异值矢量，新的规范变量可以获得.应该事先给出递推模型的阶数，计算状态和残差矢量的方程如下：

其中：x_(r)是在时间点r的状态矢量; V_x包含V_r在时间点r的前n列; z_{p (r)}是递推建模过程中新的过去观测矢量.

仅对Φ_z∈R^q×q的完整更新进行复杂度和计算成本计算.如果表示想要估计的子空间的秩，通常l < < q.许多可供选择的要求较少操作的方案被开发来避免极度高的计算复杂度.几种常见算法SVD、逆迭代、秩一更新、投影近似子空间跟踪(PAST)和RSVD-FOP的复杂度分别为O(q³)，O(q²l)，O(q³)，O(q²)，O(q²).与SVD相比，基于一阶干扰理论的递推SVD的计算成本大大减小.

2.2 确定主空间的维数

确定CVA模型阶数的方法有很多.累积方差百分比(CPV)、Akaike信息评价准则(AIC)、最大似然比检验和交叉验证等.然而，不是所有方法都适合确定递推CVA的阶数，比如交叉验证法.本文采用了累积方差百分比来确定主要的奇异值个数：

2.3 故障检测统计量及其控制限

为了监控过程故障，本文采用两个分别来自状态空间和噪声空间的统计量.基于状态空间定义T²统计量为T²=x_(r)^Tx_(r).

其置信水平α的控制限计算如下:

其中，F_{n, N-n(α)}在置信水平为α时自由度为n和N-n的F分布的临界值.

总的过程噪声统计量计算公式如下：

基于正态分布假设，Q统计量的控制限计算如下：

其中：θ_i=h_r-n，h_r是Hankel矩阵的秩；h₀=1-(2θ₁θ₃/3θ₂²)；c_α是对应1-α百分位的正态偏差.

3 应用于田纳西伊斯曼化工过程

采用田纳西伊斯曼(TE)化工过程仿真平台产生的数据进行仿真研究和分析.测试数据集使用TE过程的simulink仿真代码产生.在第300采样时刻将产品设定值从22.89改变为24.89，并持续到第800采样时刻.为了验证提出的方法，假定设定值改变属于正常操作.两种类型的过程变量故障F4，F10在第800采样时刻发生并持续到仿真时间停止.首先使用测试数据中前200个采样时刻的数据来建立CVA模型，并取得k时刻的过去观测矢量的协方差矩阵.然后，再使用测试数据进行模型更新.因为计算时间和p²成正比，在本次仿真研究中将p的值设置为3.

首先，对故障4进行仿真对比研究.图 1和图 2分别给出了基于CVA和基于RCVA-FOP方法的在线故障检测结果.

因为基于常规CVA的检测方法没有及时更新模型，所以产生了误报.如图 1所示，其监控统计量从第300采样时刻以后开始逐渐超过控制限，产生了较高的误报率.与基于CVA的检测方法对比，如图 2所示基于RCVA-FOP的在线检测方法不仅能够通过模型更新适应过程的正常改变，而且还能检测到故障4 (反应器冷却水入口温度阶跃变化)，降低了误报率，提高了过程的在线监控性能.

然后，对故障10进行仿真对比研究.图 3和图 4分别给出了基于CVA和基于RCVA-FOP方法的在线故障检测结果.

因为基于CVA的检测方法没有根据过程的正常变化及时更新模型，如图 3所示其监控统计量从第300采样时刻后开始逐渐超过控制限，增大了误报率.与基于CVA的检测方法对比，如图 4所示提出的RCVA-FOP能够通过模型更新适应过程的正常改变，同时还能检测到故障10 (流2中C的进料温度随机变化)，大大降低了误报率，提高了过程的在线监控性能.由于故障产生初期变化比较缓慢，所以存在一定的检测延时.

4 结语

将提出的递推规范变量分析应用于田纳西伊斯曼化工过程中进行仿真研究.基于一阶干扰理论的递推奇异值分解降低了算法的计算负荷.仿真结果表明，提出的方法不仅可以有效适应过程的正常改变，而且可以检测到两种类型的故障.