💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
💥1 概述
数据包含取自模拟流程示例的二维数据集。此数据用于训练和测试内核 PCA 以进行故障检测。训练后，针对输出数据空间中的每个位置计算广泛用于故障检测的T2和Q统计指标，从而生成等高线图。然后将 2% 显著性水平检测限叠加在地图上，作为数据空间的正常（绿色）和错误（洋红色）区域之间的边界。
使用等高线图，可以将各种核类型和参数选择对正常和错误过程状态之间的决策边界的影响可视化。
一、引言
在现代工业生产中，故障检测与诊断是确保设备安全、稳定运行的关键环节。随着工业过程日趋复杂，传统基于模型的故障检测方法往往难以满足需求。核主成分分析（Kernel Principal Component Analysis, KPCA）作为主成分分析（PCA）的非线性扩展，通过核函数将原始数据映射到高维特征空间，从而提取数据中的非线性特征，提高故障检测的准确性和效率。本文旨在探讨基于KPCA的故障检测方法，重点阐述T²和Q统计指数的可视化分析，并探讨其在故障诊断中的应用。
二、KPCA基本原理
KPCA的核心思想是利用核函数将原始数据映射到高维特征空间，并在该特征空间中进行PCA分析。具体步骤如下：
- 数据预处理：对输入数据进行标准化处理，以消除量纲和数值范围的影响。
- 核函数选择：常用的核函数包括线性核、多项式核、高斯核（径向基函数，RBF）和Sigmoid核等。选择合适的核函数是KPCA成功的关键，它需要根据数据的特性进行选择。
- 计算核矩阵：通过核函数计算数据之间的相似性，得到核矩阵K。
- 中心化核矩阵：去除核矩阵中的均值，以简化后续计算。
- 特征值分解：对中心化后的核矩阵进行特征值分解，得到特征值和特征向量。
- 主成分选择：选择前k个主成分构建低维子空间，用于后续的故障检测。
三、T²和Q统计指数
在基于KPCA的故障检测中，常用的统计指数包括T²统计指数（Hotelling's T-squared statistic）和Q统计指数（Squared Prediction Error, SPE），也称为Q统计量或重构误差。
-
T²统计指数：衡量样本在新主成分空间中的距离，反映了样本相对于正常样本的偏离程度。T²值越大，表明样本偏离正常范围越远，发生故障的可能性越高。其计算公式为：
T2=i=1∑k(λizi2)
其中，zi是样本在高维特征空间中第i个主成分上的投影，λi是第i个主成分对应的特征值，k是选取的用于重构的主成分个数。
-
Q统计指数：衡量样本在残差空间的距离，反映了模型重构样本的能力。Q值越大，表明模型重构样本的能力越差，样本偏离模型的程度越大，发生故障的可能性越高。其计算公式为：
Q=∥Φ(x)−Φ^(x)∥2
其中，Φ(x)是样本在高维特征空间中的投影，Φ^(x)是样本在高维特征空间中的重构。Q统计量可以由核矩阵K计算得到。
四、T²和Q统计指数的可视化分析
可视化分析在故障检测中具有重要意义，它可以帮助操作人员直观地了解系统的运行状态，及时发现异常情况。以下是T²和Q统计指数的可视化分析方法：
-
时序图：将T²和Q统计指数随时间变化绘制成时序图，可以直观地观察统计指数的变化趋势。通过设定控制限，可以检测出超过控制限的异常点，从而判断系统是否发生故障。
-
散点图：以T²和Q统计指数为坐标轴，绘制散点图，可以将数据点分布在二维空间中。正常数据点通常聚集在某个区域，而故障数据点则偏离该区域。通过分析散点图的分布，可以初步判断系统是否存在故障。
-
控制图：基于T²和Q统计指数，可以构建控制图，如Shewhart控制图、CUSUM控制图和EWMA控制图等。控制图可以用来监控系统的稳定性，及时发现系统的异常波动。
五、基于KPCA的故障检测方法实现
以下是基于KPCA的故障检测方法的Matlab实现步骤：
-
数据加载与预处理：加载训练数据和测试数据，并进行标准化处理。
-
KPCA分析：选择合适的核函数和参数，计算核矩阵并进行中心化处理。然后，对中心化后的核矩阵进行特征值分解，选择前k个主成分构建低维子空间。
-
计算T²和Q统计指数：根据KPCA分析的结果，计算训练数据和测试数据的T²和Q统计指数。
-
设定控制限：根据训练数据的T²和Q统计指数，设定控制限。常用的方法包括使用卡方分布或百分位数法。
-
故障检测：将测试数据的T²和Q统计指数与控制限进行比较，判断是否存在故障。
六、案例分析
以某化工过程为例，使用基于KPCA的故障检测方法进行故障检测。首先，收集正常运行状态下的数据样本作为训练集，利用KPCA将训练集数据映射到高维特征空间中，并计算主成分和故障检测模型。然后，将实时采集的数据样本映射到特征空间中，并利用故障检测模型进行故障检测。通过可视化分析T²和Q统计指数的变化，可以直观地了解系统的运行状态，并及时发现潜在的故障点。
七、优势与挑战
优势：
- 处理非线性关系：KPCA能够将数据映射到高维特征空间中，从而能够处理非线性关系，提高故障检测的准确性。
- 提取更有区分性的特征：KPCA通过映射到高维特征空间，能够提取更有区分性的特征，从而提高故障检测的灵敏度。
- 处理高维数据：KPCA能够处理高维数据，避免了维度灾难问题。
挑战：
- 计算复杂度高：KPCA需要计算核矩阵，计算复杂度较高，尤其是在处理大规模数据时。
- 参数选择困难：KPCA中的核函数选择和参数调节对故障检测的效果有很大影响，但如何选择合适的核函数和参数是一个挑战。
- 数据样本不平衡：在实际应用中，正常样本往往比故障样本多，导致故障检测模型容易受到正常样本的影响，降低了故障检测的准确性。
📚2 运行结果
部分代码：
%% Get 2D data
close all; clc; tic;
if nargin == 0
load dataset.mat p;
train = p{1}; test = p{2};
% Kernel types and parameters:
ktype = 'rbf'; kpar = 1; % RBF kernel
%ktype = 'rbf'; kpar = 10; % RBF kernel
%ktype = 'rbf'; kpar = 0.9; % RBF kernel
%ktype = 'rbfpoly'; kpar = [1 1 0.65]; % mixed kernel
%ktype = 'poly'; kpar = 2; % polynomial kernel
%ktype = 'imquad'; kpar = 10; % inverse multiquadric kernel
%ktype = 'cauchy'; kpar = 5; % Cauchy kernel
end
%lax = [-15 15 -15 15];
lax = [-4 10 -3 6]; % Axes limits
N = length(train); M = length(test);
z0T = train; z1T = test; % Training and Test data
[xx,yy] = meshgrid(lax(1):0.05:lax(2),... % Meshgrid for contours
lax(3):0.05:lax(4));
z2T = [xx(:) yy(:)]; L = length(z2T); % Vectorize meshgrid points
K.type = ktype; K.p = kpar; % Kernel type and parameters
set(0,'defaultfigurecolor',[1 1 1]); % Set fig color to w
conf = 0.99; % Significance level (*100%)
% Normalize 2D Data
zm = mean(z0T); zs = std(z0T);
z0 = (z0T - zm(ones(N,1),:))./zs(ones(N,1),:); % Normalize training z
z1 = (z1T - zm(ones(N,1),:))./zs(ones(M,1),:); % Normalize test z
z2 = (z2T - zm(ones(L,1),:))./zs(ones(L,1),:); % Normalize surf z
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% KERNEL PRINCIPAL COMPONENTS ANALYSIS %
[K0c,K0,U0] = kerneltrain(z0,K); % Populate kernel matrix
K1c = kerneltest(z1,z0,K0,U0,K); % Project test data to RKHS
K2c = kerneltest(z2,z0,K0,U0,K); % Project surf data to RKHS
[V,D] = eig(K0c/N); % Eigenvalue decomposition
[S,sj] = sort(diag(D),'descend'); % Sort eigenvalues
V = V(:,sj); S = S'; % Re-arrange eigenvectors
S(S < 1e-7) = []; % Remove eigenvalues <= 0
P = V(:,1:length(S))*diag(S.^-0.5); % Projection matrix
if ~isreal(S)
disp('Complex eigenvalues detected.'); % Warn about complex eigs
end
%% Perform KPCA Monitoring
CS = cumsum(S)/sum(S)*100;
RP = find(CS >= 99.9,1); % Get eigenvalues by %CPV
disp([num2str(RP) ' principal'...
' components chosen.']);
t0 = K0c*P(:,1:RP); % Kernel subspace (train)
t1 = K1c*P(:,1:RP); % Kernel subspace (test)
t2 = K2c*P(:,1:RP); % Kernel subspace (surf)
T2 = sum((t0.^2)./S(ones(N,1),1:RP),2); % T2 statistics (train)
t0n = K0c*P; % Full kernel space
Q = abs(sum(t0n.^2,2) - sum(t0.^2,2)); % Q statistics (train)
if strcmp(ktype,'rbf') == 1
fprintf('\n At infinite fault magnitude:\n');
U1 = ones(1,N)/N;
tt = U1*K0*(U0 - eye(N))*P(:,1:RP);
fprintf(' T2 limit: %.2f\n',...
sum((tt.^2)./S(1:RP),2)); % Limit of T2 for RBF
tu = U1*K0*(U0 - eye(N))*P;
fprintf(' Q limit: %.2f\n\n',...
abs(sum(tu.^2,2) - sum(tt.^2,2))); % Limit of Q for RBF
end
T2t = sum((t1.^2)./S(ones(M,1),1:RP),2); % T2 statistics (test)
T2u = sum((t2.^2)./S(ones(L,1),1:RP),2); % T2 statistics (surf)
t1n = K1c*P; t2n = K2c*P;
Qt = abs(sum(t1n.^2,2) - sum(t1.^2,2)); % Q statistics (test)
Qu = abs(sum(t2n.^2,2) - sum(t2.^2,2)); % Q statistics (surf)
%% Plot monitoring charts
figure(3); subplot(211);
semilogy(1:N,T2,'b',1:M,T2t,'m','linewidth',1.2); % T2 monitoring chart
xlabel('Time'); ylabel('T^2'); subplot(212);
semilogy(1:N,Q,'b',1:M,Qt,'m','linewidth',1.2); % Q monitoring chart
xlabel('Time'); ylabel('Q');
🎉3 参考文献
部分理论来源于网络，如有侵权请联系删除。
[1] K.E.S. Pilario, Y. Cao, and M. Shafiee. Mixed Kernel Canonical Variate Dissimilarity Analysis for Incipient Fault Monitoring in Nonlinear Dynamic Processes. Comput. and Chem. Eng., 123, 143-154. 2019. doi: 10.1016/j.compchemeng.2018.12.027
扫一扫
7839
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
