基于核主成分分析KPCA故障诊断，SPE统计量 T2统计量Matlab代码.rar

clc clear all close all %% KPCA %% 读取数据（数据的行数为采样点个数，列数为传感器个数） load trainData load testData %% 参数设置 beta =0.99; % 置信系数 PCContributionRate = 0.85; % 主元贡献率 KernelWidth = 20; % 高斯核函数参数：核宽度, ker = struct('type','gauss','width',KernelWidth); %% 预处理 trainData_mean = mean(trainData,1); % 训练数据的均值 trainData_std = std(trainData); % 训练数据的标准差 trainData_zscore = zscore(trainData);% 标准化 trainDataLength = length(trainData); %% KPCA K = computeKernelMatrix(ker,trainData_zscore,trainData_zscore); % 核矩阵 UnitTrain = ones(trainDataLength,trainDataLength)/trainDataLength; K_Centering = K-UnitTrain*K-K*UnitTrain+UnitTrain*K*UnitTrain; % 核矩阵中心化 K_Scaling = K_Centering*(trainDataLength -1)/(trace(K_Centering)); % 核矩阵方差归一化 [V,D] = eigs(K_Centering);% 特征值分解并排序 lamda = diag(D); lamdaLength = length(lamda); % 特征值个数 Vseq_Zscore = cell2mat(arrayfun(@(n) V(:,n)/(norm(V(:,n))*sqrt(lamda(n,1))), ... 1:lamdaLength , 'UniformOutput', 0)); % 特征向量标准化 % 求取主元个数 PCNumber = find(cumsum(lamda/sum(lamda)) >= PCContributionRate,1, 'first'); % PCNumber = 1; LoadingMatrix = Vseq_Zscore(:,1:PCNumber); % 求取主成分 ScoreMatrix = K_Scaling*LoadingMatrix ; % 训练样本的SPE统计量和T2统计量 % T2 统计量 T2Train = diag(ScoreMatrix*inv(diag(lamda(1:PCNumber)))*ScoreMatrix'); % SPE 统计量 SPETrain = sum((K_Scaling*Vseq_Zscore(:,PCNumber+1:lamdaLength)).^2,2); % SPE统计量和T2统计量的控制限计算（核密度估计KDE方法） [SPEf,SPExi]=ksdensity(SPETrain,'width',1,'function','cdf'); SPEIndex = find(SPEf >= PCContributionRate,1, 'first'); SPEControlLimit = SPExi(SPEIndex)*ones(trainDataLength,1); [T2f,T2xi]=ksdensity(T2Train,'width',1,'function','cdf'); T2Index = find(T2f >= PCContributionRate,1, 'first'); T2ControlLimit = T2xi(T2Index)*ones(trainDataLength,1); %% 测试数据 % 利用离线模型的均值与标准差对测试批次进行标准化处理 testData_Zscore = bsxfun(@rdivide,bsxfun(@minus,testData,trainData_mean),trainData_std); testDataLength = length(testData_Zscore); KTest = computeKernelMatrix(ker,testData_Zscore,trainData_zscore); % 测试数据与训练数据的核矩阵 UnitTest = ones(testDataLength,trainDataLength )/trainDataLength; % 测试数据的单位矩阵 KTest_Centering = KTest-UnitTest*K-KTest*UnitTrain+UnitTest*K*UnitTrain; % 核矩阵中心化 KTest_Scaling = KTest_Centering*(trainDataLength-1)/(trace(K_Centering )); % 核矩阵方差归一化 % 求取测试数据的主成分 testScoreMatrix = KTest_Scaling*LoadingMatrix ; % 测试样本的SPE统计量和T2统计量 % T2 统计量 T2 = diag(testScoreMatrix*inv(diag(lamda(1:PCNumber)))*testScoreMatrix'); % SPE 统计量 SPE = sum((KTest_Scaling*Vseq_Zscore).^2,2)-sum(testScoreMatrix.^2,2); %% 绘图 % T2 figure plot(T2ControlLimit,'k-','Linewidth',1,'MarkerSize',10,'MarkerEdgeColor','k', 'MarkerFaceColor','k') hold on plot(T2,'r-o','Linewidth',1,'MarkerSize',2,'MarkerEdgeColor','r', 'MarkerFaceColor','r'),title('T2 统计量') % SPE figure plot(SPEControlLimit,'k-','Linewidth',1,'MarkerSize',10,'MarkerEdgeColor','k', 'MarkerFaceColor','k') hold on plot(SPE,'r-o','Linewidth',1,'MarkerSize',2,'MarkerEdgeColor','r', 'MarkerFaceColor','r'),title('SPE统计量')