《MECHANICAL SYSTEMS AND SIGNAL PROCESSING》：Fault detection based on a Quality-associated virtual variable and its detectability analysis
•当质量变量不可测量时的质量相关故障检测方法。
•构建质量关联虚拟变量作为质量指标。
•扩展质量相关故障检测方法的应用范围。
•理论上，推导出故障可检测性的充分条件。
符号表
符号
B 由FastICA估计的正交矩阵
fj 第j个变量的故障幅度
gj k?(x)关于xj的梯度
I 单位矩阵
I2 (K)ICA的系统统计量
I(?) 指示函数
k(xi,xj) 径向基核函数
k(x) 样本x的核向量
k?(x) 样本x的中心化和缩放后的核向量
k?(xj) 变量xj的核向量
K 样本格拉姆矩阵
K? K的中心化矩阵
Kscl K的中心化和缩放矩阵
K? 变量格拉姆矩阵
L 块重抽样的滑动块大小
Lj 第j个变量的滞后阶数
M 自助法或置换重复次数
m 过程变量数
n X中的样本数
Num{?} 满足给定条件的样本数
pj 第j个变量的经验P值
pBH,j 经错误发现率调整的P值
P 白化矩阵
Q (K)ICA的非系统统计量
snew 新样本的独立成分
S 独立成分矩阵
tr(?) 矩阵的迹
U 格拉姆矩阵K的特征向量矩阵
V 样本协方差矩阵的特征向量矩阵
xi X中的第i个样本向量
xj X中的第j个变量向量
xnew 新的在线测试样本
xnew,r xnew中由质量相关变量组成的子向量
xnew,u xnew中由质量无关变量组成的子向量
x? 故障样本的无故障部分
X 过程变量矩阵
X? 残差矩阵
Xr X中由质量相关变量组成的子矩阵
Xu X中由质量无关变量组成的子矩阵
y 质量变量
希腊字母
α 构建的质量关联虚拟变量 (Qav)
β Lasso的回归系数向量
β?j 第j个变量的Lasso估计
β?j,b? 第b个块自助样本中第j个变量的Lasso系数
β?j,b? 中心化自助系数，β?j,b?=β?j,b??β?j
δI,δQ I2和Q统计量的控制限
δIr,δQr I2和Q统计量中质量相关部分的控制限
δIu,δQu I2和Q统计量中质量无关部分的控制限
δψ 组合统计量ψ的控制限
Λ 格拉姆矩阵Kscl的特征值矩阵
ξj 第j个变量的故障方向向量
ρj(k) 第j个变量在滞后k处的自相关系数
σ 径向基函数的核参数
ΣI I2统计量的系数矩阵
ΣQ Q统计量的系数矩阵
Σψ ψ统计量的系数矩阵
?(?) 样本级非线性映射
??(?) 变量级非线性映射
Φ 核特征空间中的数据矩阵
ψ 组合统计量，ψ=I2/δI+Q/δQ
ψr,ψu QR和QU部分的组合统计量
Ω 样本协方差矩阵的特征值矩阵
缩写
ACF 自相关函数
BH Benjamini–Hochberg
FAR 误报率
FDR 故障检测率
EHSM 电液伺服机构
ICA 独立成分分析
ICs 独立成分
IDV TEP基准的故障编号
KDE 核密度估计
KICA 核独立成分分析
KPI 关键性能指标
Lasso 最小绝对收缩和选择算子
MI 互信息
MSPM 多变量统计过程监控
Qav 质量关联虚拟变量
QavICA 基于Qav的ICA
QavKICA 基于Qav的KICA
QR 质量相关
QU 质量无关
QV 质量变量
QVs 质量变量
TEP Tennessee Eastman Process
XMEAS TEP的测量变量
XMV TEP的操作变量
1. 引言
故障检测与诊断是确保大型设备运行、工业过程安全以及维持高产品质量的关键技术[1]。与基于模型或基于知识的方法相比，数据驱动的多变量统计过程监控(MSPM)能有效从物联网数据中提取信息，减轻了对精确模型和专家知识的依赖[2]。此外，使用成熟的多变量统计分析方法使MSPM特别适合高度相关的多变量工业过程。这些优势使得MSPM在工程和学术领域都越来越受欢迎。
在大型工业过程或复杂设备中，存在各种测量变量。这些变量大致可分为两类：过程变量和质量变量(QVs)。过程变量涉及直接可观测的物理和化学量，如电压、电流、流量等，这些变量反映设备或生产的实时状态。另一方面，质量变量是指与设备性能或质量指标相关的变量，如火箭伺服喷管摆角、轧钢厚度等，这些质量变量是评估系统稳定性或产品质量的重要依据。
工业复杂系统包含各种物理、化学和生化反应，伴随物质和能量转换，因此呈现出高度复杂的特性。在系统运行过程中，收集到多种多样的过程变量。然而，并非所有过程变量都与系统的关键性能或质量指标密切相关。许多过程变量仅作为辅助操作参数，其波动对系统关键性能没有直接的因果影响；这些被称为质量无关(QU)变量。相反，直接驱动或结构性影响系统性能的变量——同时排除由共同潜在因素引起的虚假相关性——被定义为质量相关(QR)变量。对所有过程变量直接建模容易因质量无关变量的微小波动而引发误报，导致不必要的维护并降低生产效率。因此，现代工业的发展和日益激烈的竞争对故障检测提出了更严格的要求——即能够区分QR和QU故障。
QR故障检测的研究已广泛开展，大致可分为三类：(i) 预分类方法根据过程变量与QVs的相关强度将过程变量分为QR和QU组，然后分别对每组建模。代表性方法包括互信息(MI)[3], [4], [5], [6]、最小冗余最大相关性[7]、转移熵[8], [9]和投影变量重要性[10], [11]等。(ii) 后分解方法，也称为基于质量相关的回归方法，首先对所有过程变量建模形成投影空间，建立过程变量与QVs之间的回归关系，然后将投影空间分解为QR和QU子空间。代表性研究包括全偏最小二乘(TPLS)[12]、改进偏最小二乘(IPLS)[13]和全主成分回归(TPCR)[14]。最近，这种后分解方法取得了很大进展。Si等人[15]通过奇异值分解提取关键性能指标(KPI)，提出了KPI-KPLS方法，而Wang等人[16]基于过程变量与核样本之间的比例关系引入了KPCA-CCA方法。(iii) 基于深度学习的方法利用深度神经网络将过程变量和质量变量投影到高维特征空间，通过在潜在空间中的相关分析或回归提取质量相关特征。例如，Chen等人[17]提出人工神经网络相关分析，以最大化过程数据和质量指标之间的相关性，从而捕捉非线性质量关联。为解决标签稀缺问题，Ma等人[18]开发了一种半监督变分自编码器，利用未标记数据将潜在空间解耦为QR和QU子空间。Chen等人[19]进一步进行了比较研究，系统评估了各种深度神经网络辅助的典型相关分析。此外，针对特定应用场景出现了一些变体，如基于块的方法[20]、分层方法[21]和分布式方法[22]等[23], [24]。
尽管各有优势，但这三类方法的应用都依赖于一个关键假设：QVs是可测量且可靠的。具体来说，QVs对于建立数学相关性、指导空间分解或作为网络训练中的监督标签是必不可少的。然而，在实际工业场景中，恶劣的生产环境或测量技术的限制往往使某些QVs难以甚至无法获取[25]。例如，脱丁烷塔塔底丁烷浓度[26]和竖炉磁管回收率[27]。此外，虽然基于深度学习的方法可以处理强非线性，但它们通常作为“黑箱”运行，严重缺乏物理和统计可解释性。由于现有的质量相关故障检测方法依赖质量变量作为监督信号来建立过程变量与质量变量之间的关联，实时质量测量的缺失导致监控模型失去输入参考，使这些方法在实际部署中不可行。这引出了一个关键研究问题：在完全缺乏可测量QVs的情况下，是否有可能建立一个可解释的替代框架来进行变量分类和质量相关故障检测？据我们所知，针对这种未测量QV场景同时保持物理透明性的定向研究仍然相当有限，这使得它成为一个值得深入研究的问题，也是本工作的主要焦点。
此外，深入研究故障可检测性理论分析的研究也较少，特别是在非线性系统中。这种缺乏是由于引入核函数，将故障检测转移到核空间。然而，故障可检测性分析应在原始数据空间中进行，在该空间中系统的物理属性和行为更加直观和可解释。由于预映像问题[28]，核映射无法从核空间特征精确重建原始数据[29]。这种重建对于可检测性分析至关重要，因为它允许检测统计量与系统中的物理传感器直接关联。因此，预映像问题对非线性系统中的故障可检测性分析构成了重大挑战。这种理论分析的缺乏使得难以严格控制检测率和误报率之间的权衡，从而削弱了在安全关键工业中故障检测系统的可靠性。
为解决这些问题，本文提出了一种基于质量关联虚拟变量(Qav)的故障检测方法。该方法开发了一个全面的故障检测框架，包括基于Qav的过程变量分类、QR故障检测策略和故障可检测性分析。
本文的主要贡献总结如下：
(1) 提出了一个质量相关故障检测框架，消除了对在线质量变量的依赖。受工业过程物理机制的启发，我们假设存在一个主导的共享质量相关潜在因子，并讨论两种违规场景。在此假设下，构建Qav以捕获主导协同趋势。然后执行结合惩罚线性回归的筛选方法以识别QR变量组，在此基础上构建基于ICA的故障检测方案，从而拓宽了质量相关故障检测的适用性。
(2) 此外，通过引入核技巧和MI，将Qav从线性系统扩展到非线性系统。这种扩展允许开发更全面的非线性系统QR故障检测，称为基于Qav的核ICA(QavKICA)。
(3) 遵循故障重构的范式，分别在理论和实践中为线性和非线性系统推导出故障可检测性的充分条件。
本文的其余部分结构如下。第2节概述了ICA和KICA的基础知识。第3节详细说明了所提出方法的动机、Qav构建和故障检测，以及故障可检测性分析。第4节展示了在数值示例和Tennessee Eastman Process (TEP)上进行的仿真和应用。最后，第5节总结并展望了未来研究方向。
2. 预备知识
对于一个具有n个样本和m个变量的系统，假设正常操作下的测量数据为：
(1) X=[x1,x2,…,xn]∈Rm×n,
其中xi表示第i个样本，X已被标准化为零均值和单位方差。
2.1. 基于ICA的故障检测
ICA的核心思想是测量数据X可以近似为一系列独立成分(ICs)的混合，如下所示：
(2) X=AS+X?,
其中A是混合矩阵，S和X?分别是ICs矩阵和残差矩阵。
为去相关，ICA的第一步是白化，白化矩阵定义如下：
(3) P=Ω?1/2VT,
其中Ω和V分别是协方差矩阵XXT/(n?1)的特征值矩阵和特征向量矩阵。因此，测量数据X可以白化为：
(4) Z=PX=PAS=BS,
其中B是待确定的正交矩阵。由(4)式，ICs可以表示为：
(5) S=BTZ。
通常，使用FastICA来计算B，详细过程可在[30]中找到。
2.2. 基于KICA的故障检测
KICA的核心在于通过非线性变换?(.)将非线性测量xi映射到高维线性空间。然后，测量数据X被映射为：
(6) Φ=?(x1),?(x2),…,?(xn)。
由于?(.)难以显式表达，使用核技巧和径向基核函数来计算格拉姆矩阵K，其(i,j)th元素可计算如下：
(7) [K]i,j=k(xi,xj)=exp?‖xi?xj‖2/σ,
其中σ是核参数，我们有K=ΦTΦ。执行中心化和缩放如下：
(8) K?=K?1nK?K1n+1nK1n,
(9) Kscl=K?trace(K?)/n,
其中1n=1n1?1???1?1∈Rn×n。在均值中心化和缩放后，对Kscl进行特征值分解，得到特征值矩阵Λ和特征向量矩阵V，然后白化矩阵可定义为：
(10) P=nΦUΛ?1,
并且Φ可以被白化为：
(11) Z=PTΦ=nΛ?1UTKscl。
白化后，使用FastICA计算B，以最大化线性组合sj=bjTZ的非高斯性，样本xnew的ICs可表示为：
(12) snew=D1/2BTPT?(xnew),
其中D是Kscl特征值矩阵中的d阶顺序主子矩阵。KICA的详细描述可参考文献[31]。
3. 提出的方法
3.1. 问题描述与动机
在复杂设备或工业过程中，质量相关(QR)变量并非相互独立。通过系统内物质、能量和信息的不断传递与反馈，这些变量在正常运行过程中发生动态耦合，并趋向于协调演化。相比之下，质量无关（QU）变量通常不参与这一共同模式。这种结构耦合具有重要的实际意义：产品质量的影响通常表现为在大量被测过程变量中共享的主导性共同运动。
当质量相关变量（QV）难以直接测量时，一个自然的问题随之而来：能否构造一个单一的虚拟变量，作为系统质量特性的综合度量？我们将此虚拟变量称为Qav，它与所有QR变量显著相关，同时与所有QU变量不相关。如果能够构造出这样的Qav，那么质量相关故障检测就不再依赖于QV本身：我们可以通过评估每个过程变量与Qav之间相关性的强度来准确分类QR变量，然后对QR组进行质量相关故障检测。这一想法的有效性基于以下物理假设。
假设
在正常运行条件下，过程变量共享一个与产品质量相关的主导潜在因子，而QU变量和测量噪声则在多个分散方向上波动，因此它们不参与这一主导模式。
该假设描述了过程固有的结构特性。它并非普遍成立，至少在两种情况下可能被违反：
情况1：一组相互强相关的QU变量表现出比质量相关部分更大的组合变异性。在这种情况下，主导共同运动被QU块捕获，Qav反映的是工况变化而非产品质量。
情况2：质量相关的变异性相对于测量噪声或缓慢漂移而言较弱。此时主导模式变得不稳定，可能被噪声主导，导致所得变量划分不可靠。在实践中，在稳态工况下收集训练数据并施加轻度平滑可以缓解此问题，假设检验的显著性水平也可相应调整。
基于这一思想，所提出的基于Qav的故障检测方法的总体框架如图1所示。
下载：下载高分辨率图像（567KB）
下载：下载全尺寸图像
图1. 基于Qav的质量相关故障检测的总体框架。
3.2. 基于QavICA的线性系统QR故障检测
3.2.1. Qav的构造
对于具有n个样本和m个变量的系统，假设正常运行下的测量数据为：
(13) X=[x1,x2,…,xn]=[x1;x2;...;xm]∈Rm×n.
其中xi表示第i个样本，xj表示第j个变量。为了构造Qav α∈Rn，目标函数是最大化α与所有变量的内积平方和，即：
(14) α=argmax∑j=1m(xjαn?1)2, s.t.αTα=n?1.
使用拉格朗日乘数法，式(14)可转化为：
(15) L(α)=∑j=1m(xjαn?1)2+λ(n?1?αTα),
其中λ为拉格朗日乘数。对L(α)关于α求偏导：
(16) ?L(α)?α=2(n?1)2∑j=1m(xjα)(xj)T?2λα=2(n?1)2[(x1α)(x1)T+?+(xmα)(xm)T]?2λα.
令?L(α)/?α=0，式(14)中的优化问题可转化为以下特征值分解问题：
(17) XTXα(n?1)2=λα, (s.t.αTα=n?1).
选取最大特征值对应的特征向量作为Qav。尽管Qav在数学上与归一化第一主成分成比例，但其物理和操作角色有所不同：它不是表示用于降维的传统样本主成分得分，而是以变量方式从过程数据矩阵中导出。具体而言，它充当捕获过程变量间主导共同运动模式的变量空间载荷方向，而非减少样本维度。为了直观理解这种等价性，我们提供了过程数据结构的几何解释，如图2所示。
考虑高维向量空间中过程变量的波动方向。质量相关波动（红色箭头）由潜在的耦合机制驱动，趋于协调变化。从几何上看，它们的变异向量高度对齐，形成建设性叠加，产生一个具有大幅度的主导合向量（向量r）。相比之下，质量无关波动和噪声（黑色箭头）通常呈现随机或弱相关；它们的变异向量分散在各个方向上，求和时大部分相互抵消，产生一个幅度可忽略的合向量（向量u）。因此，整体波动方向（向量v）——即第一主成分——绝大多数由协调的质量相关簇决定。这种几何主导性确保了第一主成分有效捕获主要的质量相关潜在因子，使其与背景噪声区分开来。
值得注意的是，如果QU变量彼此强相关（例如由负荷或工况驱动的变量），它们的合向量将不再微小，并可能接管主导方向——这是第3.1节讨论的失效情况。
下载：下载高分辨率图像（345KB）
下载：下载全尺寸图像
图2. Qav作为主导方向的几何示意图。
3.2.2. 基于Qav的变量分类和QR故障检测
所提出的QR故障检测的第一步是将过程变量分类为QR和QU组。为了消除由共享潜在变量引起的伪相关，同时严格控制时间自相关和多重检验偏差，我们提出了一种集成数据驱动自相关分析、多元惩罚回归和非参数重抽样的变量筛选框架。
工业时间序列数据通常违反独立同分布假设。为了在后续非参数检验中准确保留局部时间依赖结构，首先基于自相关函数（ACF）的衰减特性，以无监督、数据驱动的方式确定最优滑动块大小L。对于中心化过程变量矩阵X∈Rm×n中的每个特征xj，计算其在滞后阶l的样本自相关系数：
(18) ρj(l)=∑i=1n?l(xii?x?j)(xi+lj?x?j)∑i=1n(xii?x?j)2.
对于无自相关的白噪声序列，其ACF值应以95%的概率落在±1.96/n的渐近置信区间内。定位每个特征首次落入该置信区间的滞后阶Lj：
(19) Lj=min{l | |ρj(k)|<1.96/n}.
提取所有特征截断阶数的中位数作为全局最优滑动块大小L=ceil(median(L1,L2,…,Lm))，确保该块长度充分覆盖系统的主导时间动态记忆周期。
为了避免单变量边际相关检验导致的严重伪相关和循环推理陷阱，采用Lasso回归对所有变量进行多元联合评估。Lasso通过L1正则化评估每个特征对Qav的条件结构贡献，并将冗余特征的效应量压缩至零。目标函数定义为：
(20) β?=argminβ(12n∑i=1n(αi?xiTβ)2+λ∑j=1m|βj|),
其中最优惩罚参数λ通过K折交叉验证获得。所得原始稀疏回归系数向量记为β?。
基于块大小L，对样本序列执行移动块Bootstrap方法，构造保留原始时间特征的伪样本集。从有效起始索引中有效回地随机抽取[n/L]个索引，提取相应的连续数据块，按顺序拼接并截断至原始样本长度n，生成具有相同自相关强度的伪样本对(Xb*,yb*)。该重采样过程重复M次。在每次重采样迭代中，使用固定的优化λ对伪样本拟合Lasso模型，获得重采样系数向量β?b*。
为进行统计假设检验（零假设H0:βj=0），需要对Bootstrap经验分布进行中心化，使其均值为零，模拟零假设成立的状态：
(21) β?j,b*=β?j,b*?β?j.
对每个变量xj，计算原始观测系数β?j在中心化零分布中的双尾概率作为经验P值：
(22) pj=1M∑b=1MI(|β?j,b*|≥|β?j|),
其中I(?)为示性函数。
为控制多重假设检验引起的假阳性率膨胀，采用Benjamini-Hochberg（BH）程序对m个变量的P值进行错误发现率校正。将P值按升序排列p(1)≤p(2)≤?≤p(m)，并从后向前迭代计算校正P值：
(23) pBH,(j)=min{k≥i}min(1, m/k p(k)).
若pBH,(j)<γ（例如显著性水平γ=0.05），则拒绝零假设，最终将该变量识别为qr变量；其余变量标记为qu变量。 为此，分别对qr部分xr和qu部分xu构建第2.1节所述的ica模型。在线检测时，对于给定样本xnew，可将其分类为qr部分xnew,r和qu部分xnew,u。通常定义两个统计量i2和q来监测系统性和非系统性变化。以qr部分为例，xnew,r的i2和q统计量计算如下： (24) i2(xnew,r)=xnew,rTPTBBTPxnew,r, (25) q(xnew,r)=xnew,rTPT(I?BBT)Pxnew,r, 其中i表示单位矩阵。考虑到分布的不确定性，采用核密度估计（kde）确定控制限δir和δqr。 显然，i2和q以互补方式运行。为简化检测逻辑，引入组合统计量ψ： (26) ψr(xnew)=I2(xnew,r)δIr+Q(xnew,r)δQr. 类似地，可以获得qu统计量ψu(xnew)。其控制限δψr和δψu可通过kde获得。现在我们将qavica的详细步骤总结为算法1。 下载：下载高分辨率图像（321kb） 下载：下载全尺寸图像 3.2.3. 故障可检测性分析 故障可检测性分析基于几何描述[32]，其中故障样本由无故障部分x*叠加故障部分组成。发生在第j个传感器中的故障公式化为： (27) x=x*+ξjfj, 其中ξj=[0,…,1,…,0]T∈Rm。 定理1 对于线性系统，式(27)描述的样本x被ψ统计量检测到的充分条件如下： (28) fj>2δψ/Σψ,jj 或 fj其中
(29) Σψ=PTBBTPδI+PT(I?BBT)PδQ,
ψ∈{ψr,ψu}，δI∈{δIr,δIu}，δQ∈{δQr,δQu}，Σψ,jj指矩阵Σψ的第j行第j列元素。
证明
对于样本x，其ψ统计量可简化为：
(30) ψ(x)=I2(x)δI+Q(x)δQ=xT(PTBBTPδI+PT(I?BBT)PδQ)x=‖Σψ1/2x‖2.
x被检测到的充分条件是minψ(x)>δψ。将(27)代入(30)得：
(31) ‖Σψ1/2x‖2=‖Σψ1/2(x*+ξjfj)‖2.
利用向量范数的逆三角不等式，可得：
(32) ‖Σψ1/2(x*+ξjfj)‖≥|‖Σψ1/2ξjfj‖?‖Σψ1/2x*‖|.
考虑故障条件，故障项的主导效应远大于正常状态。当故障幅值fj足够大（满足‖Σψ1/2ξjfj‖≥‖Σψ1/2x*‖）时，上述表达式可放宽为：
(33) ‖Σψ1/2(x*+ξjfj)‖≥‖Σψ1/2ξjfj‖?‖Σψ1/2x*‖.
由于上述不等式两边均为非负，两边平方得：
(34) ‖Σψ1/2(x*+ξjfj)‖2≥(|fj|Σψjj?‖Σψ1/2x*‖)2=(|fj|Σψjj?‖Σψ1/2x*‖)2.
由于x*是x的无故障部分，满足‖Σψ1/2x*‖≤δψ。因此：
(35) ‖Σψ1/2(x*+ξjfj)‖2≥(|fj|Σψjj?δψ)2.
因此，样本x=x*+ξjfj被可靠检测的充分条件是ψ(x)>δψ，即：
(36) (|fj|Σψjj?δψ)2>δψ.
考虑到|fj|Σψjj>δψ且δψ>0，两边开平方并整理得：|fj|Σψ,jj>2δψ。最终，故障幅值必须满足的充分条件为：
fj>2δψ/Σψ,jj 或 fj因此，定理1得证。 ■
3.3. 基于QavKICA的非线性系统QR故障检测
3.3.1. Qav的构造
式(14)–(17)默认变量之间存在线性关系，这显然不适用于非线性系统。对于非线性系统，最直接的思路是通过非线性映射?*(·)将变量xj投影到线性空间，即xj→?*(xj)。这样，?*(xj)呈现线性关系，式(14)–(17)转化为：
(37) α=argmax∑j=1m(?*(xj)αm?1)2, s.t.αTα=m?1,
(38) L(α)=∑j=1m(?*(xj)αm?1)2+λ(m?1?αTα),
(39) ?L(α)?α=2(m?1)2∑j=1m(?*(xj)α)?*(xj)T?2λα,
(40) Φ*TΦ*α(m?1)2=λα, (s.t.αTα=m?1),
其中Φ*=[?*(x1)T,…,?*(xm)T]。由于难以显式表达?(·)，通常采用核技巧。引入Gram矩阵K*=Φ*TΦ*，其元素定义为：
(41) kij*=??*(xi),?*(xj)?=exp(?‖xi?xj‖2σ).
注意，此处K*∈Rm×m是变量方向的Gram矩阵。为避免核方法中固有的原像问题，Qav的构造和后续变量分类均在Hilbert空间中进行。因此，在非线性系统中，Qav α是m维向量。式(37)转化为：
(42) K*α(m?1)2=λα, (s.t.αTα=m?1).
与线性系统一致，选取最大特征值对应的特征向量作为Qav α。
3.3.2. 基于Qav的变量分类和QR故障检测
在处理工业过程数据中的非线性特征选择时，为消除时间自相关同时克服复杂非线性约束，提出了一种基于互信息（MI）和块置换检验的非线性筛选方法。MI可表示为两个变量的联合分布p(xi,xj)与其边缘分布p(xi)p(xj)之积的相对熵：
(43) I(xi,xj)=?p(xi,xj)logp(xi,xj)p(xi)p(xj)dxidxj.
引入信息熵后，式(43)可限制为：
(44) I(xi,xj)=H(xi)+H(xj)?H(xi,xj),
(45) H(xi)=?∫p(xi)logp(xi)dxi,
(46) H(xi,xj)=??p(xi,xj)logp(xi,xj)dxidxj.
MI越大，变量间相关性越强；反之，MI越小，变量间相关性越弱。为稳健捕获非线性依赖而不做限制性分布假设，采用Moon等人[33]提出的KDE方法估计MI。
与线性系统类似，采用ACF衰减准则确定最优块大小L。由于通过核方法映射回原始数据空间涉及原像问题，我们的方法直接利用核向量 k?(xj)=[k1j?,k2j?,…,kmj?]T∈Rm 来评估核域中与 α∈Rm（Qav）之间的非线性依赖（通过互信息）。由于工业过程变量在连续时间步之间通常表现出强自相关，因此采用块置换来构建统计上有效的零分布。具体而言，执行 M 次块置换以生成置换后的互信息序列 I(i)k?(xj);αi=1M，由此计算经验 p 值：
(47)pj=1M∑Mi=1II(i)k?(xj);α≥Ik?(xj);α。
与线性系统类似，采用 BH 准则根据 (23) 校正 p 值，得到调整后的值 pBH,j。如果 pBH,j 小于显著性水平 γ（例如 γ=0.05），则拒绝零假设，变量 xj 被分类为 QR 变量；否则，被分类为 QU 变量。
此时，测量数据 X 可以明确分为两组：QR 部分 Xr 和 QU 部分 Xu。随后，分别对 Xr 和 Xu 构建第 2.2 节所述的 KICA 模型。给定一个标准测试样本 xnew，按照上述方法将其分解为 xnew,r 和 xnew,u。以 xnew,r 为例，我们计算核向量 knew,r=[k(x1,r,xnew,r),…,k(xn,r,xnew,r)]T（其中 xi,r 是 Xr 中的第 i 个样本），并执行中心化和缩放如下：
(48)k?new,r=knew,r?K1new?1nknew,r+1nK1newtrace(K?)/n，
其中 1new=1n1,…,1∈Rn×1。定义统计量 I2 和 Q 以监测系统性和非系统性变化：
(49)I2(xnew,r)=k?new,rTnUΛ?1BBTΛ?1UTk?new,r，
(50)Q(xnew,r)=k?new,rTnUΛ?1I?BBTΛ?1UTk?new,r，
其中 I 是单位矩阵。采用核密度估计确定控制限 δIr 和 δQr。组合统计量 ψ 可按 (26) 计算，其控制限 δψr 可通过 KDE 获得。类似地，可得到 QU 统计量 ψu(xnew) 及其控制限 δψu。
现在我们将 QavKICA 的详细步骤总结为算法 2。
下载：下载高分辨率图像 (302KB)
下载：下载全尺寸图像
3.3.3. 故障可检测性分析
定理 2
对于非线性系统，式 (22) 中描述的样本 X 能被 ψ 统计量检测到的充分条件如下：
(51)fj?b+b2?4ac2a，
其中 a=gjTΣψgjgj?null(Σψ)，b=2k?(x?)TΣψgj，c=k?(x?)TΣψk?(x?)?δψ。这里，Σψ=nUΛ?1BBTΛ?1UTδIr+nUΛ?1(I?BBT)Λ?1UTδQr 是一个半正定系数矩阵；gj 表示 k?(x) 关于第 j 个变量 xj 的梯度；k?(x) 和 k?(x?) 分别是 x 和 x? 的居中且缩放后的核向量；gj 的第 t 个元素可计算为：
(52)(gj)t=trace(K?)n2(xt)j?(x)jσk(x,xt)?∑i=1n1new2(xt)j?(x)jσk(x,xi)，
其中 xt 是训练数据集 X 中的第 t 个样本，而 (?)j 表示向量的第 j 个元素。
证明
对于样本 x，其 ψ 统计量可简化为：
(53)ψ(x)=I2(x)δI+Q(x)δQ=k?(x)TΣψk?(x)。
x 能被检测到的充分条件是 ψ(x)>δψ。
为了评估第 j 个变量中故障幅值 fj 的影响（x=x?+ξjfj），我们使用一阶泰勒级数在 x? 的邻域内展开 k?(x)：
(54)k?(x?+ξjfj)≈k?(x?)+gjfj。
将此近似代入 (53) 得到：
(55)ψ(x)≈k?(x?)+gjfjTΣψk?(x?)+gjfj=k?(x?)TΣψk?(x?)+2k?(x?)TΣψgjfj+gjTΣψgjfj2。
代入系数 a、b 和 c 的定义，检测条件 ψ(x)>δψ 转化为标准一元二次不等式：
(56)afj2+bfj+c>0。
由于 x? 是无故障基线部分，其统计量满足 k?(x?)TΣψk?(x?)<δψ，这保证了常数项 c<0。此外，在假设 gj?null(σψ) 下，确保 a>0。因此，b2?4ac>0 成立。
因此，要使不等式 ψ(x)>δψ 成立，故障幅值 fj 必须满足：
fj?b+b2?4ac2a。
因此，定理 2 得证。 ■
注记
在上述理论推导中，x? 作为故障样本 x 的无故障基线。在实际工业部署中，x? 可以通过计算离线正常训练数据的均值向量（适用于静态过程）或采用最近邻基线匹配方案（适用于多模式过程）来获得。
4. 仿真与应用研究
在本节中，通过数值仿真和 TEP 评估所提方法的有效性。在 QR 故障检测实验中，采用误报率（FAR）和故障检测率（FDR）来评估不同方法的性能，其计算公式如下：
(57)FAR=Numψr(x)>δψrx∈QU faultNumx∈QU fault，
(58)FDR=Numψr(x)>δψrx∈faultNumx∈fault， 对于 QR 故障
Numψu(x)>δψux∈faultNumx∈fault， 对于 QU 故障，
其中 Num{.} 表示满足条件的样本数量。
4.1. 数值示例仿真
考虑以下线性数值示例 [12]：
(59)x=Az+e,y=Cx+v，
其中 z∈R3，zj～U[0,1]，e～N(0,0.052)，v～N(0,0.012)，e 和 v 是随机噪声，x∈R7 和 y∈R1 分别是过程变量和质量变量。A=[1 3 4 4 0 0 0; 3 1 1 4 1 1 1; 1 1 3 0 0 2 3]T，C=[2 3 1 1 0 0 0]。因此，x1?x4 是 QR 变量，而 x5?x7 是 QU 变量。
4.1.1. 相关性分析和变量分类
根据 (59)，生成 500 个正常样本。为了验证 (17) 中构建的 Qav α 是否能有效反映 (59) 中质量变量 y 的变化，我们对两者进行了标准化并绘制了它们的趋势。图 3 显示，虽然 α 的曲线与 y 的曲线不完全重叠，但总体趋势显示出高度一致性。该结果表明 Qav 能有效跟踪质量变量的变化。此外，为了量化每个变量的时间依赖性并确定最优滑动块大小，我们对 (59) 中的过程变量 x1 到 x7 进行了自相关函数（ACF）衰减分析。结果如图 4 所示。红色虚线表示白噪声序列的 95% 渐近置信界；一旦变量的 ACF 曲线首次落入这些界内，其当前观测值在某种意义上不再显著受过去值影响。如图 4 所示，所有七个过程变量在低滞后阶数（Lag < 5）时都表现出强自相关，ACF 值远高于 0.2。到 Lag=7 时，大多数变量 x1,x2,x3,x4,x6,x7 的自相关系数首次进入并稳定在 95% 置信区间内。因此，后续移动块自助法的块大小设置为 L=7，如图中黑色实线所示。
为了验证所提方法中变量划分策略的有效性，我们进行了一组对比实验；结果如图 5 所示。在这些实验中，零假设下的零分布通过自助重采样构建，每个过程变量的显著性通过结合假发现率控制的假设检验来评估。作为基于真实质量变量的基线，图 5(a) 显示了真实质量变量与过程变量之间相关性检验的基本结果。如图所示，x1 到 x4 的观测系数（红色实线）明显位于双侧置信界（红色虚线）之外，统计量显示 P=0.000。因此，这四个变量被正确划分为 QR 变量。相反，x5 到 x7 的系数被方法直接收缩为零（黑色线），给出 P=1.000，它们被分类为 QU 变量。
下载：下载高分辨率图像 (373KB)
下载：下载全尺寸图像
图 3. 质量变量和 Qav 的趋势。
下载：下载高分辨率图像 (438KB)
下载：下载全尺寸图像
图 4. 过程变量的 ACF 衰减。
图 5(b) 显示了验证结果：使用 Qav 运行相同的假设检验过程。变量 x1 到 x4 仍表现出强显著性（P=0.000）并被分类为 QR 变量，而 x5 到 x7 均未通过显著性检验——例如，x5 的观测系数落在零假设的高概率区域内（P=0.465），而 x6 和 x7 被收缩为零。因此，三者都被正确分配到 QU 组。比较图 5(a) 和 (b) 的统计结果，由 Qav 得出的变量分类（x1?x4 为 QR，x5?x7 为 QU）与使用质量变量获得的基线完全匹配。这种一致性表明 Qav 可以作为质量变量的替代，从而支持后续的质量相关故障检测。
下载：下载高分辨率图像 (694KB)
下载：下载全尺寸图像
图 5. 数值示例中的划分结果。(a) 直接使用质量变量 y。(b) 使用所提出的 Qav。
4.1.2. 变量分类的稳定性实验
所提方法的一个关键前提是过程变量表现出单一主导的共同变化方向，该方向可以可靠地识别。为了验证这一假设，我们在图 6 中展示了由正常条件训练数据计算的 XTX 的特征值谱。最大特征值 λ1=8.82×10?3 明显占主导：它大约是第二特征值 λ2=1.78×10?3 的 4.95 倍，对应的相对特征间隙 (λ1?λ2)/λ1 约为 80%。第三特征值之后迅速衰减至零。这种谱模式表明过程变化集中在少数几个主导方向上，第一个方向捕捉主导模式。该主导模式对应于质量相关的共享潜在因子；λ1 和 λ2 之间的清晰分离证实了主导方向在正常操作条件下是可识别且稳定的，这证明了使用 Qav 作为不可测量质量变量的替代是合理的，并支持了变量划分的可靠性。
此外，宽特征间隙表明主导方向对正常操作条件下训练数据的小扰动不敏感，这在下文敏感性分析中进一步量化。
为了评估所提方法对样本变化的鲁棒性，对 Qav 进行了自助稳定性分析。我们执行了有放回的自助重采样，并计算了原始向量与自助向量之间的余弦相似度。
下载：下载高分辨率图像 (198KB)
下载：下载全尺寸图像
图 6. 正常操作条件下 XTX 的特征值谱。
如图 7 所示，余弦相似度的分布严重偏向于 1，平均值为 0.9992。几何上，这对应于约 2.15° 的平均角度偏差。相似度值的高度集中（范围从 0.995 到 1.0）表明估计的主成分方向非常稳定，对采样扰动不敏感。这些结果证实了提取的特征结构对于实际应用是稳健且可靠的。
下载：下载高分辨率图像 (190KB)
下载：下载全尺寸图像
图 7. 余弦相似度的分布。
4.1.3. 基于 Qav 的质量相关故障检测
为了测试所提 QavICA 的故障检测性能，我们选择了 TPCR [14] 和 IPLS [13] 作为比较方法。TPCR 和 IPLS 都是该领域广泛认可的经典 QR 故障检测方法。TPCR 的主成分数量基于累积方差贡献率确定，而 IPLS 则通过交叉验证确定。为了保持一致性和公平性，所有参与比较的方法的置信水平统一设置为 0.95。有两种故障场景描述如下：
故障 1：x3 上的阶跃偏置：x=x?+00120000。
故障 2：x5 上的阶跃偏置：x=x?+0000200。
对于每个故障场景，我们生成了 1000 个故障样本，并使用三种不同的方法进行检测。为了全面评估这些方法的性能，我们进行了 200 次重复实验，并对每种方法的 FAR 和 FDR 进行了直方图分析，结果如图 8 所示。对于 QR 故障（故障 1），图 8(a) 显示在绝大多数重复实验中，所提出的 QavICA 方法实现了接近 100% 的 FDR，显著优于其他两种方法，三种方法的平均 FDR 分别为 94.36%、97.16% 和 99.50%。相反，对于 QU 故障（故障 2），图 8(b) 表明 QavICA 方法的 FAR 极低，接近 0%，三种方法的平均 FAR 分别为 13.06%、1.70% 和 0.85%。直方图分析结果清楚地表明，所提出的 QavICA 方法在 QU 故障中具有比其他方法更低的 FAR，而在 QR 故障中其 FDR 显著高于其他方法，从而验证了 QavICA 方法的高鲁棒性。
下载：下载高分辨率图像 (281KB)
下载：下载全尺寸图像
图 8. 重复实验中的分布直方图。(a) 故障 1 的 FDR。(b) 故障 2 的 FAR。
4.1.4. 可检测性分析的验证
根据第 3.2.2 节的结果，将包含 500 个正常样本的矩阵 X 划分为 QR 部分 Xr 和 QU 部分 Xu，并按照第 2.1 节所述分别建立 ICA 模型。然后，根据 (29)，系数矩阵 Σψr 和 Σψu 可计算为：
Σψr=4.132511.1746?7.7221?7.672711.174633.9469?23.5310?22.0967?7.7221?23.531016.553915.1349?7.6727?22.096715.134914.9642，
Σψu=0.4954?2.28172.1105?2.281715.5035?14.64082.1105?14.640813.9857。
控制限为 δψr=2.0403 和 δψu=1.8007。根据定理 1，对于变量 x1，式 (27) 中描述的样本 x 能被 ψr 统计量检测到的充分条件是 f1>1.4053 或 f1同样地，对于变量x2，充分条件为f2>0.4903或f2
