PCA故障诊断步骤剖析

基于PCA算法的故障诊断步骤

离线PCA监测模型的计算步骤：

（1） 选择监控变量，收集正常工况下的各变量的样本，记为训

练样本数据X_train和检验数据X_test；

X_train为n×m矩阵，即n个样本，m个观测变量（即以列向量来看

的话，为一个观测变量各个采样点的值）

对样本数据X_train和检验数据X_test进行中心化和标准化处理

1 / 12

PCA故障诊断步骤剖析

得到和；

中心化处理：按列对X_train减去观测变量的均值

观测变量某一采样点的值减去这一观测变量所有采样点的平均值

求取一列（即某一观测变量）的平均值

标准化处理：对X_train除以观测变量的标准差（按列（观测变量）

进行）

标准差

求出标准化矩阵的协方差矩阵∑；

的协方差矩阵对∑为：

（2） 对∑进行特征分解,求得特征值

（

荷向量）；

（3） 确定主元个数, 确定了主元个数k,就得到了k个特征值

2 / 12

）及其对应的特征向量（负

PCA故障诊断步骤剖析

,及其对应的特征向量

A:累计贡献率法：

；

前k个主元的累积方差贡献率为：

当前k个主元的累积方差贡献率达到85%，则主元个数取k值 B:交叉检验估计法：

将采集到的数据分成k个部分，1部分数据用来建立主元模型,剩下的k-1部分用来作为检验数据去检验所建的模型。如此,建立若干个不同主元个数的模型,并测试所建立的模型,从中选取一个通过检验后误差最小的模型的主元个数作为系统主元个数。 （4） 建立PCA主元模型，并进行交叉验证以确定误差最小 按照

,求出第i个主元，并依据

求出其主元模型 用

带入得到另一主元模型

，依据

,求出模型误

差，确定模型误差最小的那个模型即为主元模型。 （5） 计算T统计量控制限和SPE统计量控制限；

对于样本个数为n,主元个数为k的过程变量X_train, T统计量服从自由度为k和n一k的F分布,则置信度为а的T统计量控

3 / 12

2

2

2

PCA故障诊断步骤剖析

制上限为:

或

检验水平为а的SPE统计量控制上限为：

,,,

是与（1-）分位点对应的标准差 在线过程监测与故障诊断步骤： （1） 采集第i时刻的在线实时数据

进行中心化和标准化处理得到

（2） 按照

得分向量，依据计量

（3） 计算

,这里

2

（为1×m矩阵），并；

,求出

的

,求出PCA模型估；

2

的T统计量和SPE统计量，并画出T统计量和SPE

统计量的控制图；

（4） 将上述计算结果与T统计量控制限和SPE统计量控制限比

较,以检测过程运行有无异常，当有异常状态发生时,绘

4 / 12

2

PCA故障诊断步骤剖析

制贡献图，找出与故障相关的系统变量：

1) 检查每个观测值x的标准化得分, 并确定造成失控

状态的r(r2) 计算每个变量相对于失控得分的贡献率是:

3) 当

是负时,设它为零；

4) 计算第j个过程变量的总贡献率：

5) 把所有m个过程变量的

画在一个曲线图上。

PCA_TE仿真程序：

%%TE过程的传统主元分析在Matlab中的仿真程序 %建立模型：

%载入模型数据，以故障11为例

Xtrain = load('G:d11.dat'); Xtrain = double(Xtrain);

%载入测试数据

Xtest = load('G:d11_te.dat');

5 / 12

PCA故障诊断步骤剖析

Xtest = double(Xtest);

%标准化处理：

X_mean = mean(Xtrain); %按列求Xtrain平均值 X_std = std(Xtrain); %求标准差

[X_row,X_col] = size(Xtrain); %求Xtrain行、列数 % for i = 1:X_col

%Xtrain(:,i)=(Xtrain(:,i) - X_mean(i)./X_std(i)); %Xtest(:,i) = (Xtest(:,i) - X_mean(i)./X_std(i)); % end Xtrain=(Xtrain-repmat(X_mean,X_row,1))./repmat(X_std,X_row,1);

%求协方差矩阵

sigmaXtrain = cov(Xtrain);

%对协方差矩阵进行特征分解，lamda为特征值构成的对角阵，T的列为单位特征向量，且与lamda中的特征值一一对应：

[T,lamda] = eig(sigmaXtrain); % disp('特征根（由小到大）'); % disp(lamda); % disp('特征向量：');

6 / 12

PCA故障诊断步骤剖析

% disp(T);

%取对角元素(结果为一列向量)，即lamda值，并上下反转使其从大到小排列，主元个数初值为1，若累计贡献率小于90%则增加主元个数

D = flipud(diag(lamda)); num_pc = 1; while sum(D(1:num_pc))/sum(D) < 0.9

num_pc = num_pc +1;

end

%取与lamda相对应的特征向量

P = T(:,X_col-num_pc+1:X_col);

%求置信度为99%、95%时的T2统计控制限 T2UCL1=num_pc*(X_row-1)*(X_row+1)*finv(0.99,num_pc,X_row - num_pc)/(X_row*(X_row - num_pc)); T2UCL2=num_pc*(X_row-1)*(X_row+1)*finv(0.95,num_pc,X_row - num_pc)/(X_row*(X_row - num_pc));

%置信度为99%的Q统计控制限 for i = 1:3

7 / 12

PCA故障诊断步骤剖析

theta(i) = sum((D(num_pc+1:X_col)).^i); end

h0 = 1 - 2*theta(1)*theta(3)/(3*theta(2)^2); ca = norminv(0.99,0,1);

QUCL = theta(1)*(h0*ca*sqrt(2*theta(2))/theta(1)

+ 1 + theta(2)*h0*(h0 - 1)/theta(1)^2)^(1/h0);

%在线监测： %标准化处理 n = size(Xtest,1);

Xtest=(Xtest-repmat(X_mean,n,1))./repmat(X_std,n,1);

%求T2统计量，Q统计量 [r,y] = size(P*P'); I = eye(r,y);

T2 = zeros(n,1); Q = zeros(n,1); for i = 1:n

T2(i)=Xtest(i,:)*P*inv(lamda(52-num_pc+1:52,

52-num_pc+1:52))*P'*Xtest(i,:)';

8 / 12

PCA故障诊断步骤剖析

Q(i) = Xtest(i,:)*(I - P*P')*Xtest(i,:)'; end %绘图 figure

subplot(2,1,1);

plot(1:n,T2,'k'); title('主元分析统计量变化图'); xlabel('采样数'); ylabel('T^2'); hold on;

line([0,n],[T2UCL1,T2UCL1],'LineStyle','--','Color','r');

line([0,n],[T2UCL2,T2UCL2],'LineStyle','--','Color','g');

subplot(2,1,2); plot(1:n,Q,'k'); xlabel('采样数'); ylabel('Q');

9 / 12

PCA故障诊断步骤剖析

hold on;

line([0,n],[QUCL,QUCL],'LineStyle','--','Color','r'); %贡献图

%1.确定造成失控状态的得分 S = Xtest(400,:)*P(:,1:num_pc); r = [ ]; for i = 1:num_pc

if S(i)^2/lamda(i) > T2UCL1/num_pc r = cat(2,r,i); end end

%2.计算每个变量相对于上述失控得分的贡献cont = zeros(length(r),52); for i = length(r) for j = 1:52 cont(i,j)

abs(S(i)/D(i)*P(j,i)*Xtest(400,j)); end

10 / 12

=

PCA故障诊断步骤剖析

end

%3.计算每个变量的总贡献 CONTJ = zeros(52,1); for j = 1:52

CONTJ(j) = sum(cont(:,j)); end

%4.计算每个变量对Q的贡献 e = Xtest(400,:)*(I - P*P'); contq = e.^2;

%5. 绘制贡献图 figure;

subplot(2,1,1); bar(CONTJ,'k'); xlabel('变量号'); ylabel('T^2贡献率 %');

subplot(2,1,2); bar(contq,'k'); xlabel('变量号');

11 / 12

PCA故障诊断步骤剖析

ylabel('Q贡献率 %');

12 / 12