基于一种改进稀疏动态慢特征分析的高速公路服务区空调空气处理单元故障检测研究

Fault Detection of Air Conditioning Unit in Expressway Service Area Based on an Improved Sparse Dynamic Slow Feature Analysis

摘要: 高速公路服务区暖通空调空气处理系统表现出很强的双向动态特性，为了通过处理空气处理系统的双向动态特性和对提取的潜在变量施加稀疏性，本文提出一种改进的稀疏动态慢特征分析策略来检测空气处理系统的故障。在提出的稀疏动态慢特征分析中，采用自回归移动平均模型来揭示变量之间的自相关关系。然后应用多路数据分析，通过将扩充的三维数据集转换为展开的矩阵，计算出在多个批处理运行中的分批动态特性。进一步建立动态慢特征分析模型充分处理批运行中的时间动态特性。最后，融入特征稀疏表示技术，通过对负载向量进行稀疏约束，消除了无意义变量之间的耦合。在ASHRAE研究项目RP-1312实验数据集上进行的案例研究验证了所提出的故障检测方案的有效性。

Abstract: HVAC air handling systems in expressway service areas exhibit strong bidirectional dynamics. In order to detect the faults of the air handling system by dealing with the bidirectional dynamics of the air handling system and imposing the sparsity on the extracted latent variables, this paper proposes an improved sparse dynamic slow feature analysis strategy. In the proposed sparse dynamic slow feature analysis, an autoregressive moving average model is employed to reveal the autocorrelation among the variables. Multi-way data analysis is then applied to handle the batch dynamics over multiple batch runs by converting the augmented 3D dataset into an unfolded matrix. A dynamic slow feature analysis model is further established to fully deal with the time dynamic characteristics in batch operation. Finally, feature sparse representation technology is incorporated, and the coupling between meaningless variables is eliminated by sparsely constraining the load vector. A case study conducted on the ASHRAE research project RP-1312 experimental dataset verifies the effectiveness of the proposed fault detection scheme.

1. 介绍

作为高速公路服务区暖通空调的重要模块，暖通空调系统的空气处理系统发生故障会降低室内舒适度，缩短设备寿命，浪费建筑能源 [1] 。因此，为了监测和排除空气处理系统的故障，研究一种有效的方法对空气处理系统进行故障检测是非常必要的。随着空气处理机组系统记录的运行数据越来越丰富，基于数据驱动的故障检测方法已经成为最近的研究热点 [2] 。

基于数据驱动的空调处理系统故障检测方法主要有两类：人工智能方法和多元统计分析方法 [1] [2] 。人工智能方法已经引起了越来越多的关注，因为它们拥有强大的能力能将原始数据样本转换到一个新的基于内部非线性映射的特征空间。例如，Lee等人 [3] 提出了一种基于深度神经网络的故障检测技术来识别空气处理系统故障。考虑到训练样本不足的问题，Yan等人 [4] 提出了一种改进的支持向量机来处理不平衡数据。然而，在模型训练阶段，大多数人工智能方法使用整个原始变量作为模型输入。这种缺陷会导致故障检测效果不够准确，会包括一些冗余变量和过程干扰。相反，多元统计分析方法从原始过程数据中提取潜在变量信息，捕获过程变化的本质原因，具有使用降维技术消除噪声和冗余变量的优势。因此，近年来，多元统计分析方法已被用于实现空气处理系统故障检测，如基于线性因果模型的主成分分析 [5] ，基于贝叶斯网络的非线性主成分分析 [6] 等。然而，上述大多数故障诊断方法的模型都是建立在空气处理系统稳态运行的基础上 [1] [2] [4] 。这些故障检测方法具有捕获动态空气处理系统的时变特性的缺点。

由于室外和室内环境条件的变化以及空气处理系统在相邻的几天内反复运行特性，空气处理系统具有双向动态特性 [7] [8] ：时间型和批次型动态。因此，迫切需要提出一种动态故障诊断方法来处理空气处理系统的双向时变动态特性。为了解决时间上的动态性质，出现了一些动态故障检测技术。丁等 [9] 提出了一种改进的加权移动平均技术来解决多故障发生的情况。Karami等人 [10] 讨论了一种增强的高斯混合模型来处理时变动态特性。然而，现有的故障检测技术不足以解决空气处理系统的时间动态特性，因为导出的潜在变量在不同的时间间隔内是独立的 [8] 。此外，由于训练数据总是在短时间的操作中被记录，而不是在长时间的运行中被收集，所以如何处理空气处理系统的批次间动态特性目前还处于研究阶段。

近年来，基于慢特征分析的方法已经被用来处理批次过程的时间动态特性，因为它能够挖掘慢变化的潜在变量，这些变量可以捕捉潜在的驱动力。例如，基于增强核慢特征分析的方法 [8] ，改进的统计慢特征分析算法 [11] 等。这些方法更适合挖掘时间动态特性，因为它们可以表征过程的潜在驱动力。因此通过采用慢特征分析来处理空气处理系统的时间动态，故障检测的性能能够得到显著增强。为了进一步揭示变量的自相关关系，自回归移动平均模型 [12] [13] 被集成到基于慢特征分析的方法中，以使用其先前的观察值来增加每个输入样本。另一方面，为了有效地计算出分批动态变化和偏差，多路数据分析 [14] 被融合到建立的检测模型中。基于展开的二维训练矩阵，构造了动态慢特征分析模型，以更充分地处理空气处理系统的时间动态。因此，通过结合动态慢特征分析模型和多路数据分析，提出基于三维数据的动态慢特征分析方法，可以有效地处理空气处理系统的两维时变动态特性。

通过进一步分析，潜在变量的稀疏性对提高监控模型的故障检测有效性具有关键意义 [15] 。不具有稀疏性的潜在变量应用于空气处理系统故障检测时，它们不能揭示原始变量之间的显著相关性。这一缺点不仅阻碍了故障检测结果的物理解释，而且无法消除原始变量之间的冗余耦合关系。作为一种新的表示学习方法，稀疏表示技术 [16] 最近在故障检测领域获得了很多关注。例如，Luo等人 [16] 提出了稀疏稳健主成分分析，以产生可解释的主成分分析，并增强模型的稳健性。为了选择关键故障变量，Yu等人 [17] 提出了一种基于稀疏指数判别分析的方法。Zhang等人 [18] 在卷积自动编码器网络中嵌入了一个稀疏正则化层，以消除提取的特征冗余信息。受上述基于稀疏表示的故障检测方法的优越性和成功应用的启发，论文对提出的三维数据动态慢特征分析模型施加了稀疏限制，以有效地保持故障区分特征并改善空气处理系统故障检测性能。

为了有效处理空气处理系统的两维时变动态特性，本文提出了一种新的稀疏动态慢特征分析方法，检测高速公路服务区空气处理系统故障，进一步提高故障检测的有效性。在提出的方法中，三维为了揭示变量的自相关关系，首先集成自回归移动平均外生模型，在每个批次数据集中扩充每个样本，然后应用多路数据分析来解决多批次运行中的批次动态性质三维。最后，在展开的二维训练数据集上，建立动态慢特征分析模型。为了提高模型的可解释性和故障检测能力，进一步构建稀疏动态慢特征分析算法来提取稀疏慢特征，利用索套惩罚将稀疏约束施加到动态慢特征分析模型中，对负荷向量的非零元素施加惩罚，以将无关紧要的变量的系数收缩为零。然后利用稀疏广义特征值技术求解动态慢特征分析的优化问题。在项目RP-1312的实验数据集上对所提出的空气处理系统故障检测方案进行了详细的实验，并与相关方法进行了比较。

2. 预备知识

2.1. 空气处理系统

作为供暖、通风和空调系统的一个重要设备，变风量空气处理系统 [19] 目前被广泛使用。在变风量空气处理系统中，使用加热和冷却盘管对送风进行加热、冷却和除湿。然后，供应空气通过管道系统输入到多个区域的空气处理终端。通过回风管，将多个房间的空气送回。送风和回风风扇用于促进管道系统中的空气循环。送风静压控制器、风门控制器、送风温度控制器和回风流量控制器用于保持合适的送风压力和温度以及充足的室外空气通风。通过一些合理的控制策略，送风静压、送风温度、回风流量差控制回路的设定值可以根据实际运行情况适时地重新设定。有缺陷的系统设计和不正确的工艺操作会导致控制器和传感器出现故障。根据部位和原因不同，典型故障有风门卡滞、传感器测量偏差和阀门泄漏等。然而，目前很难进行有效的空气处理系统故障检测，因为空气处理系统的批处理运行中存在时间动态行为 [8] 。此外，在连续几天内，空气处理系统的特点是不同批次运行之间的分批动态。因此，需要解决空气处理系统故障检测的两维时变动态问题。此外，在传统的故障检测模型中，通常应用不具有稀疏性的潜在变量来检测和检测空气处理系统故障，削弱了监控模型的故障敏感性和检测能力 [20] 。

2.2. 慢特征分析算法

慢特征分析的原理是产生相对于模型输入具有缓慢变化的输出。对于 $x\left(t\right)={\left[x_1\left(t\right),\cdots ,x_m\left(t\right)\right]}^T$ ，慢特征分析利用变换函数 $g\left(x\right)={\left[g_1\left(x\right),\cdots ,g_J\left(x\right)\right]}^{\text{T}}$ 以确保输出的 $y\left(t\right)={\left[y_1\left(t\right),\cdots ,y_J\left(t\right)\right]}^{\text{T}}$ 具有最慢的变化。其中， $y_j\left(t\right)=g_j\left(x\left(t\right)\right)$ 。慢特征分析的目标函数表达式为：

$\text{min}\langle {\left({\stackrel{\dot{}}{y}}_j\left(t\right)\right)}^2\rangle$ (1)

s.t.

$\langle y_j\left(t\right)\rangle =0$ (2)

$\langle y_j^2\left(t\right)\rangle =1$ (3)

$\forall i\ne j:\langle y_i\left(t\right)y_j\left(t\right)\rangle =0$ (4)

其中， $y_j\left(t\right)$ 表示第j个慢速特征， ${\stackrel{\dot{}}{y}}_j\left(t\right)$ 表示 $y_j\left(t\right)$ 的导数，运算符 $\langle \cdot \rangle$ 计算如下：

$\langle f\left(t\right)\rangle =\frac{1}{t_1-t_0}{\displaystyle {\int }_{t_0}^{t_1}f\left(t\right)dt}$ (5)

等式(1)保证输出变量具有最慢的变化，等式(2)被设计成确保输出变量具有零均值。等式(3)用于去除不变解，而等式(4)确保输出变量相互正交。 $y_1\left(t\right)$ 拥有最慢的变化， $y_2\left(t\right)$ 拥有第二慢的变化，依此类推。

在 $g_j(\cdot )$ 是线性映射的情况下， $y_j\left(t\right)$ 被公式化为 $y_j\left(t\right)=w_j^{\text{T}}x\left(t\right)$ ，其中 $w_j$ 代表载荷向量。慢特征分析优化后被重写为：

$\begin{array}{l}\min \langle {\stackrel{\dot{}}{y}}_j^2\left(t\right)\rangle =\min w_j^{\text{T}}\langle \stackrel{\dot{}}{x}\left(t\right){\stackrel{\dot{}}{x}}^{\text{T}}\left(t\right)\rangle w_j=\min w_j^{\text{T}}A{w}_j\\ s.t.\langle y_j^2\left(t\right)\rangle =w_j^{\text{T}}\langle x\left(t\right)x^{\text{T}}\left(t\right)\rangle w_j=w_j^{\text{T}}B{w}_j=1\end{array}$ (6)

其中，A和B分别表示 $\stackrel{\dot{}}{x}\left(t\right)$ 和 $x\left(t\right)$ 的协方差。 $\stackrel{\dot{}}{x}\left(t\right)$ 近似计算为 [9] ：

$\stackrel{\dot{}}{x}\left(t\right)=\frac{dx\left(t\right)}{dt}\approx x\left(t\right)-x\left(t-1\right),t=2,3,\cdots ,T$ (7)

所定义的慢特征分析模型的目标函数被进一步转换成特征值分解问题，表达示为：

$A{w}_j={\lambda }_{j}B{w}_j$ (8)

3. 改进的稀疏动态慢特征分析故障检测模型

通过对空气处理系统运行的分析，我们发现空气处理系统的动态特征表现为时间和批次方向上的两维动态特性 [8] 。然而，如前所述，针对传统的空气处理系统故障检测方法无法解决两维时变动态特性，以及提取的潜在变量是所有原始变量的线性组合，不具有稀疏性的问题，我们提出一种新的稀疏动态慢特征分析方法。该方法首先采用自回归移动平均外生模型来扩充每个批处理运行数据集，以便捕获变量的自相关关系；然后对多个扩充的正常分批运行数据集执行多路数据分析以处理批次间动态特性，并使用展开的训练矩阵建立动态慢特征分析模型以解决时间动态特性；最后将特征稀疏表示技术集成到动态慢特征分析方法中，进一步构建基于稀疏动态慢特征分析模型，以获得用于故障检测的稀疏慢特征。

三维自回归移动平均外生模型的目标是用以前的观察值增加一个批处理运行数据集中的每个样本。具体而言，第i批运行数据集X_i的增广矩阵是通过用其先前的d个样本扩展当前样本来构建的，给出如下 $X_i^d$

$X_i^d=\left[\begin{array}{cccc}{x}^{\text{T}}\left(k\right)& x^{\text{T}}\left(k-1\right)& \cdots & x^{\text{T}}\left(k-d\right)\\ x^{\text{T}}\left(k-1\right)& x^{\text{T}}\left(k-2\right)& \cdots & x^{\text{T}}\left(k-d-1\right)\\ ⋮& ⋮& \ddots & ⋮\\ x^{\text{T}}\left(k+d-K\right)& x^{\text{T}}\left(k+d-K-1\right)& \cdots & x^{\text{T}}\left(k-K\right)\end{array}\right]$ (9)

其中， $x\left(k\right)$ 表示第k个样本，K是样本数，d表示时滞。在文献 [21] 和 [22] 中，时滞d = 2被证明适合于揭示变量的自相关关系。因此，在本文中，时滞d也根据经验确定为2。

3.1. 基于三维数据的动态慢特征分析

3.1.1. 三维训练数据分析

假设空气处理系统存储了相邻几天的正常运行数据，获取的多批数据集表示为X_i (J × K)，其中i = 1，2，……，I，K和J分别表示每批运行的样本数和变量数，I表示批运行的次数。I个正常批量数据集X_i (J × K)被自回归移动平均外生模型扩充，以获得扩展矩阵X_i^d (J_d × K)，i = 1，2，……，I，其中J_d表示具有时滞d的扩充变量的数量。

Figure 1. Data dynamic expansion and 3D data analysis diagram

图1. 数据动态扩充和三维数据分析图

如图1所示，首先利用所有这些扩展矩阵X_i^d (J_d × K)来构造三维矩阵X^d (I × J_d × K)。然后通过多路数据分析将三维矩阵转换成两路矩阵，构建改进的基于稀疏三维数据的动态慢特征分析故障检测模型。

首先，将矩阵X^d (I × J_d × K)划分成k = 1，2，……，K的矩阵X_k (I × J_d)，然后，通过按照分批方向放置K个矩阵X_k (I × J_d)，将三维训练数据集X^d (I × J_d × K)展开成两路数据集 ${\tilde{X}}^d\left(I\times K{J}_d\right)$ 。其次，对训练数据集 ${\tilde{X}}^d\left(I\times K{J}_d\right)$ 进行归一化处理，给定 ${\tilde{X}}^d\left(I\times K{J}_d\right)=\left[{\tilde{x}}_1^{\text{T}},{\tilde{x}}_2^{\text{T}},\cdots ,{\tilde{x}}_I^{\text{T}}\right]$ ，其中 ${\tilde{x}}_i^{}\in {ℝ}^{K{J}_d\times 1}$ 表示数据集 ${\tilde{X}}^d(I\times K{J}_d)$ 的第i个样本，该样本 $\tilde{x}{}_i$ 被归一化为：

$x_i\text{=}\frac{{\tilde{x}}_i-mean\left({\tilde{X}}^d\left(I\times K{J}_d\right)\right)}{std\left({\tilde{X}}^d\left(I\times K{J}_d\right)\right)},i=1,2,\cdots ,I$ (10)

其中，mean ()表示平均值运算符，std ()表示标准差运算符。

通过这样的归一化向量，归一化矩阵被建立为 $X^d\left(I\times K{J}_d\right)=\left[x_1^{\text{T}},x_2^{\text{T}},\cdots ,x_I^{\text{T}}\right]$ 。最后一步，利用变量展开操作将归一化数据集X^d (I × KJ_d)转换为新的数据集X^d (IK × J_d)。

3.1.2. 建立动态慢特征分析模型

由于慢特征分析模型在揭示潜在驱动力方面的出色能力，使用展开的双向数据集X^d (IK × J_d)建立了动态慢特征分析模型，以充分解决空气处理系统的时间动态特性。

设 $x_k\left(i\right)\in {ℝ}^{J_d}$ 代表矩阵X_k (I × J_d)中的第i个样本。通过使用K-最近邻(K-NN)准则从矩阵X_k (I × J_d)中选择其l个最近邻样本，记为x_k^j (i)，其中j = 1，2，……，l。与样本x_k (i)相关的伪时间序列被建立为：

$t_k\left(i\right)={\left[x_k\left(i\right),x_k^1\left(i\right),x_k\left(i\right),x_k^2\left(i\right),\cdots ,x_k\left(i\right),x_k^l\left(i\right)\right]}^{\text{T}}\in {ℝ}^{2l\times J_d}$ (11)

矩阵X_k (I × J_d)的伪时间序列可以构建为 $T_k\left(2Il\times J_d\right)={\left[t_k\left(1\right),t_k\left(2\right),\cdots ,t_k\left(I\right)\right]}^{\text{T}}$ 。同样，可以重建如下数据集X^d (IK × J_d)的伪时间序列T (2IKl × J_d)：

$T(2IKl\times J_d)={[T_1^{\text{T}},T_2^{\text{T}},\cdots ,T_K^{\text{T}}]}^{\text{T}}$ (12)

矩阵T (2IKl ×J_d)中的第(2i − 1)个和第2i个样本分别表示为 $\tau \left(2i-1\right)\in {ℝ}^{J_d}$ 和 $\tau \left(2i\right)\in {ℝ}^{J_d}$ 。因此，时间变化信息被认为出现在 $\tau \left(2i-1\right)$ 和 $\tau \left(2i\right)$ 中。在下面的小节中，矩阵T (2IKl × J_d)和X^d (IK × J_d)将用于建立三维数据的动态慢特征分析模型故障检测模型。

利用构造的矩阵T (2IKl × J_d)和展开的训练矩阵X^d (IK × J_d)，首先通过求解下面的优化问题建立基于三维数据的动态慢特征分析模型。

$\begin{array}{l}{J}_{TBDSFA}=\underset{w}{\text{min}}\frac{1}{IKl-1}{w}_{}^{\text{T}}{\displaystyle \underset{i=1}{\overset{IKl}{\sum }}\left(\tau \left(2i\right)-\tau \left(2i-1\right)\right){\left(\tau \left(2i\right)-\tau \left(2i-1\right)\right)}^{\text{T}}}w\\ s.t.\frac{1}{IK-1}{w}_{}^{\mathrm{T}}{\displaystyle \underset{i=1}{\overset{IK}{\sum }}x\left(i\right)x_{}^{\mathrm{T}}\left(i\right)}w=1\end{array}$ (13)

其中，向量w表示负载向量，x_i表示第i个样本。

矩阵T (2IKl × J_d)的时间变化矩阵 $\Delta T$ 计算如下：

$\Delta T=\left[\Delta \tau \left(1\right),\Delta \tau \left(2\right),\cdots ,\Delta \tau \left(IKl\right)\right]$ (14)

其中， $\Delta \tau \left(i\right)=\tau \left(2i\right)-\tau \left(2i-1\right)$ $i=1,2,\mathrm{...},IKl$。

进一步计算时间变化协方差矩阵 $L=\frac{1}{IKl-1}\Delta T\Delta T^{\text{T}}$ 和协方差矩阵 $C=\frac{1}{IK-1}{X}^d{}^{\text{T}}{X}^d$ 。等式(13)中的动

态慢特征分析模型的目标函数改写为：

$\begin{array}{l}{J}_{TBDSFA}=\underset{w}{\text{min}}\frac{1}{IKl-1}{w}_{}^{\text{T}}{\displaystyle \underset{i=1}{\overset{IKl}{\sum }}\Delta \tau (i)\Delta \tau {(i)}^{\text{T}}}w=\underset{w_{}}{\text{min}}{w}_{}^{\text{T}}\frac{1}{IKl-1}\Delta T\Delta T^{\text{T}}w=\underset{w_{}}{\text{min}}{w}_{}^{\text{T}}Lw\\ s.t.\frac{1}{IK-1}{w}_{}^{\mathrm{T}}{\displaystyle \underset{i=1}{\overset{IK}{\sum }}x(i)x_{}^{\mathrm{T}}(i)}w=w_{}^{\mathrm{T}}\frac{1}{IK-1}{X}^d{}^{\text{T}}{X}^{d}w=w_{}^{\mathrm{T}}Cw=1\end{array}$ (15)

求解下面的特征值分解问题解决上述优化目标：

$Lw=\lambda Cw$ (16)

收集p个最小特征值的特征向量建立负荷矩阵。

3.2. 建立稀疏动态慢特征分析的故障检测模型

为了消除变量之间无意义的耦合关系，进一步提高故障检测能力，我们将特征稀疏表示技术引入到动态慢特征分析模型中，提取稀疏慢特征信息。在目标函数中对负荷向量的元素施加惩罚，把无关紧要变量的系数收缩为零，把系数非零的变量确定为关键变量。因此，构造的稀疏动态慢特征分析算法可以生成一系列包含较少非零元素的稀疏负荷向量，以探索空气处理系统的潜在驱动力。

在构建的基于三维数据的动态慢特征分析的优化目标上施加额外的约束 $\rho {\Vert w\Vert }_0$ ，计算稀疏负荷向量。基于稀疏动态慢特征分析的目标函数构造如下：

$\begin{array}{l}{J}_{STBDSFA}=\underset{w}{\max }{w}_{}^{\text{T}}\hat{L}w-\rho {\Vert w\Vert }_0\\ s.t.w_{}^{\mathrm{T}}Cw=1\end{array}$ (17)

其中， $\hat{L}=-L=-\frac{1}{IKl-1}\Delta T\Delta T^{\text{T}}$ ， $\rho >0$ 是用于调整负荷向量w_j稀疏度的松弛因子。

为了解决方程(17)中的优化是不连续和非凸问题，利用文献 [16] 中讨论的松弛技术将方程(17)转换成如下的连续凸优化问题：

$\begin{array}{l}\underset{w}{\max }{w}_{}^{\text{T}}\hat{L}{w}_{}-{\rho }_{\epsilon }{\displaystyle \underset{i=1}{\overset{J_d}{\sum }}\log \left(\epsilon +\left|w_i\right|\right)}\\ s.t.w_{}^{\mathrm{T}}Cw\le 1\end{array}$ (18)

其中， ${\rho }_{\epsilon }=\rho /\log \left(1+{\epsilon }^{-1}\right)$ ，w_i表示负载向量w的第i个元素。

选择一个参数 $\beta \ge \max \left(0,-{\lambda }_{\min }\left(\hat{L}\right)\right)$ 满足 $\hat{L}+\beta I_{J_d}\in S_{+}^{J_d}$ ，等式(18)可以重新表述为：

$\begin{array}{l}\underset{w}{\min }\beta {\Vert w\Vert }_2^2-w_{}^{\text{T}}\left(\hat{L}+\beta I_{J_d}\right)w+{\rho }_{\epsilon }{\displaystyle \underset{i=1}{\overset{J_d}{\sum }}\log \left(\epsilon +\left|w_i\right|\right)}\\ s.t.w_{}^{\mathrm{T}}Cw\le 1\end{array}$ (19)

根据文献 [16] ，方程(19)的解可以通过计算下面的二次规划问题来计算。

$\begin{array}{l}{w}^{\left(j+1\right)}=\arg \underset{w}{\min }\beta {\Vert w\Vert }_2^2-2w_{}^{\text{T}}\left(\hat{L}+\beta I_{J_d}\right)w^{\left(j\right)}+{\rho }_{\epsilon }{\displaystyle \underset{i=1}{\overset{J_d}{\sum }}\frac{\left|w_i\right|}{\left|w_i^{\left(j\right)}+\epsilon \right|}}\\ s.t.w_j^{\mathrm{T}}C{w}_j\le 1\end{array}$ (20)

其中，j代表迭代次数。

稀疏动态慢特征分析的第一个稀疏负荷向量可以由稀疏广义特征值算法直接计算。为了求解其它稀疏负荷向量，采用了正交投影收缩技术，基于先前的稀疏负荷向量对矩阵 $\hat{L}$ 进行收缩 [17] 。

$\begin{array}{l}{\hat{L}}_{j+1}=\left(I-v_j{v}_j^{\text{T}}\right){\hat{L}}_j\left(I-v_j{v}_j^{\text{T}}\right)\\ v_{j+1}=\frac{\left(I-V_j{V}_j^{\text{T}}\right)w_j}{\Vert \left(I-V_j{V}_j^{\text{T}}\right)w_j\Vert }\end{array}$ (21)

其中， $w_1,w_2,\mathrm{...},w_j$ 表示通过对收缩矩阵执行稀疏广义特征值方法计算稀疏负荷向量。基于三维数据的动态慢特征分析模型的特征向量被确定为实施过程和稀疏广义特征值算法的初始解。

改进的稀疏动态慢特征分析模型生成了一组对应广义特征值 ${\lambda }_1^{}<{\lambda }_2^{}<\cdots <{\lambda }_{J_d}^{}$ 的稀疏广义特征向量 $A_{IK}=\left[w_1,w_2,\cdots ,w_{J_d}\right]$ 。这些广义特征向量可以表征稀疏慢特征的变化趋势。保留太多稀疏的慢特征会将不必要的注意力集中在短期系统扰动上，导致基于稀疏动态慢特征分析的监控效果恶化。因此，为了保持空气处理系统的关键特征信息，我们只选择前p个重要的广义特征向量 $A_p=\left[w_1,w_2,\cdots ,w_p\right]$ 来导出变化最慢的稀疏慢特征。通过计算稀疏慢特征的累积慢度贡献率 $R_p$ 来确定参数p [8] 。

$R_p={\displaystyle \underset{j=1}{\overset{p}{\sum }}{\mu }_j}/{\displaystyle \underset{j=1}{\overset{J_d}{\sum }}{\mu }_j}$ (22)

其中， ${\mu }_j$ 定义为 ${\mu }_j=1/{\lambda }_j^{}$ 。在改进的稀疏动态慢特征分析模型中，p值是根据95%的累积满度贡献率确定的。

给定新批次运行中的标准化测试样本 $x_{te}$ ，首先使用自回归移动平均模型构建其时间滞后向量 $x_{te}^d$ 。然后计算稀疏慢特征 $y_{te}=A_p^{\text{T}}{x}_{te}^d$ ，依据扩展的测试样本 $x_{te}^d$ ，计算T²统计量衡量空气处理系统的主要变化趋势；计算SPE统计量捕获短期系统扰动。

$T^2=y_{te}^{\text{T}}{y}_{te}=x_{te}^{d\text{T}}{A}_p{A}_p^{\text{T}}{x}_{te}^d$ (23)

$SPE=y_{total}^{\text{T}}{y}_{total}-y_{te}^{\text{T}}{y}_{te}=x_{te}^{d\text{T}}{A}_{IK}{A}_{IK}^{\text{T}}{x}_{te}^d-x_{te}^{d\text{T}}{A}_p{A}_p^{\text{T}}{x}_{te}^d$ (24)

其中， $y_{total}=A_{IK}^{\text{T}}{x}_{te}^d$ 。

4. 实验和比较

4.1. 实验数据和参数设置介绍

基于改进的稀疏动态慢特征分析模型的故障检测方案在ASHRAE研究项目RP-1312提供的空气处理系统数据集上进行了检验。该项目进行了一系列现场实验，以模拟服务于四个区域的空气处理系统的时变动态特性。因此，ASHRAE项目RP-1312可以提供多种故障模式和正常运行的实验数据集。关于项目RP-1312的更多细节可在文献 [21] 中找到。

Table 1. Introduction of normal data sets

表1. 正常数据集的介绍

Table 2. Basic information about the fault mode of the air handling system

表2. 空气处理系统故障模式的基本信息

包含空气处理系统变量、设定值、控制信号和致动器状态的16个关键变量被选为监控变量。为了处理空气处理系统在长期运行期间的分批动态，利用表1中给出的13个正常数据集N₁~N₁₃建立三维建模矩阵，这些正常运行批次是从十三个不同的运行日收集的。如表2所示，在本实验中采用了多种故障模式来评估所开发的故障检测方案的有效性。这些不同的故障模式包括冷却盘管阀卡滞、加热盘管阀泄漏和室外空气风门卡滞或泄漏。

原始正常和故障实验数据的采样间隔为1分钟，我们以3分钟的间隔对正常和故障实验数据进行重采样，以降低故障检测方法的计算复杂度。在数据重采样阶段，16个监测变量被添加了高斯测量噪声。为了评估构建的稀疏动态慢特征分析方法的故障检测效果，将相关故障批次数据集中的前120个正常样本与前360个故障样本组合成测试数据集。因此，表2所示的故障模式出现在测试数据集中的第121个样本。

将所提出改进的稀疏动态慢特征分析方法的有效性与稀疏主成分分析和稀疏慢特征分析进行对比。对于稀疏主成分分析和稀疏慢特征分析模型，仅使用一个在运行日存储的正常数据集来建立训练数据。对于改进的稀疏动态慢特征分析方法，自回归移动平均外生模型中使用的时滞d选择为2。在改进的稀疏动态慢特征分析模型、稀疏慢特征分析模型和稀疏主成分分析模型中，保留了占95%贡献的主成分，T²和SPE的控制限由99%的置信水平决定。

为了评估方法的故障检测效果，使用故障检测率和故障检测时间作为评价指标。故障检测率指数计算为检测到的故障数据点与实际故障数据点总数的比率。连续3个样本超出监测统计的控制极限后，发现故障。故障检测时间指数被定义为最早检测到的3个连续故障样本的第一个数量。故障检测时间越早，故障检测率越高，检测效果越好。

4.2. 故障检测结果的比较

(1) 故障F5的故障检测结果

Figure 2. Fault detection diagram of three methods for fault F5

图2. 对于故障F5的三种方法的故障检测图

故障F5的改进稀疏动态慢特征分析、稀疏慢特征分析和稀疏主成分分析的监控图如图2所示。图2(a)中的稀疏主成分分析模型的故障检测结果揭示了T²和SPE统计分别在第130和第127样本检测到故障F5。然而，稀疏主成分分析模型的T²和SPE统计都具有最低的故障检测率，因为在第130次和第127次采样后，大量故障数据点低于其相应的控制极限。根据图2(b)，稀疏慢特征分析模型获得改进的监控结果，其中其统计T²在第241个样本警告故障F5，而其统计SPE在第124个样本检测故障F5。然而，稀疏慢特征分析模型的统计SPE仍然拥有相对较大的漏检测故障样本，因为许多故障数据点在第124个数据点之后下降到控制极限以下。与上述两种方法的故障检测图相比，所提出的改进稀疏动态慢特征分析模型显示了更好的监控结果。从图2(c)可以看出，其T²和SPE统计在第121个数据点同时发现故障F5，漏检故障数据点最少，表明改进的稀疏动态慢特征分析模型实现了最早和最好的故障检测效果。因此，在三种监测方法中，改进的稀疏动态慢特征分析模型是检测故障F5最有效的方法。

(2) 故障F10的故障检测结果

Figure 3. Fault detection diagram of three methods for fault F10

图3. 对于故障F10的三种方法的故障检测图

图3显示了故障F10的检测结果。如图3(a)所示，稀疏主成分分析模型获得了最差的故障检测图表，其中其T²统计在第237个样本处警告故障F10，而SPE统计在第121个样本处检测到故障F10。此外，稀疏主成分分析模型的T²和SPE统计都带来很低的故障检测率，因为在检测到故障后，许多故障数据点低于控制极限。相反，从图3(b)中，稀疏慢特征分析模型的T²和SPE统计分别在第121和第171个样本处警告故障F10。然而，由于大量真实故障数据点被错误地认为是正常数据点，它们仍然导致不好的故障检测结果。图3(c)中所示的改进稀疏动态慢特征分析模型的T²和SPE统计数据在第121个数据点同时发出F10故障警告，而没有遗漏检测到的故障数据点，这在三种检测方法中最早且最精确地对F10故障发出报警。仿真结果证明了改进的稀疏动态慢特征分析模型相对于其他两种模型具有良好的故障检测效果。

4.3. 故障F1~F11的故障检测结果

表3和表4评估了11个故障的改进稀疏动态慢特征分析模型、稀疏慢特征分析模型和稀疏主成分分析模型的检测结果，所有测试故障都是阶跃变化故障。因此，如表3所示，这三种方法具有在第121个样本引入后立即警告大多数测试故障的能力。更具体地说，稀疏主成分分析模型导致故障F1、F2、F6和F8~F11的一些故障检测延迟，而它没有故障检测延迟来警告其余故障。稀疏慢特征分析模型的T²统计量仅具有故障检测延迟以警告故障F3、F6和F8，而其SPE统计量具有时间延迟所以给出故障F2~F8、F10和F11的警报。相反，改进的稀疏动态慢特征分析模型实现了几乎没有故障检测延迟来警告这11个测试故障，这与稀疏慢特征分析模型和稀疏主成分分析模型相比更早的显示了故障检测效果。

然而，表4中列出这三种方法的故障检测率指数的不同值。由于无法处理空气处理系统的两维动态特性，稀疏慢特征分析模型和稀疏主成分分析模型在检测到故障后都错误地将真实故障数据点视为正常数据点。这导致T²和SPE统计数据的故障检测率较低。改进的稀疏动态慢特征分析模型的T²和SPE统计实现了几乎100%的故障检测率来检测11个测试故障，这是三种检测方法中最高的。这是因为改进的稀疏动态慢特征分析模型可以解决两维动态、执行空气处理系统的特征稀疏表示。表3和表4的结果也再次证明了改进的稀疏动态慢特征分析模型相对于其他两种方法的显著检测性能。

Table 3. Comparison of fault detection time of the three algorithms

表3. 三种算法的故障检测时间比较

Table 4. Comparison of fault detection rates of the three algorithms

表4. 三种算法的故障检测率比较

5. 结论

本文提出了一种改进的稀疏动态慢特征分析模型的高速公路服务区空气处理系统故障检测方法。这是第一次将特征稀疏表示技术集成到基于慢特征分析算法的方法中来监测空气处理系统。在提出的方法中，首先针对空气处理系统的两维动态特性，提出了一种基于三维数据的动态慢特征分析模型。具体而言，多路数据分析用于处理分批动态特性，动态慢特征分析融合自回归移动平均模型用于计算时间动态特性以及自相关关系。进一步提出了基于稀疏三维数据的动态慢特征分析算法，增强了模型的可解释性，提高了故障检测能力。通过将特征稀疏表示技术集成到三维数据的动态慢特征分析模型中，对负荷向量的非零元素施加惩罚，将无关变量的系数收缩为零。实验和比较结果证明了所提出的改进稀疏动态慢特征分析算法具有良好的监控效果。

NOTES

^*第一作者。

^#通讯作者。