一、程序框架与核心步骤
%% 主程序入口
clear; clc; close all;%% 1. 数据加载与预处理
% 加载西储大学轴承数据集（示例）
data = load('bearing_fault_data.mat'); % 包含振动信号和标签
X = data.signals; % 振动信号矩阵 (N_samples × N_features)
y = data.labels; % 故障标签 (0:正常,1:内圈,2:外圈,3:滚珠)% 数据标准化（零均值单位方差）
[X_norm, mu, sigma] = zscore(X);%% 2. PCA降维实现
[coeff, score, latent, explained] = pca(X_norm);% 选择主成分（累积方差贡献率≥95%）
cum_explained = cumsum(explained);
n_components = find(cum_explained >= 95, 1);
X_pca = score(:,1:n_components);%% 3. 特征可视化与故障模式分析
figure;
gscatter(X_pca(:,1), X_pca(:,2), y);
xlabel('PC1'); ylabel('PC2');
title('轴承故障模式在PC空间的分布');%% 4. 分类模型训练与评估
% 划分训练集与测试集
cv = cvpartition(y,'HoldOut',0.3);
X_train = X_pca(cv.training,:);
y_train = y(cv.training);
X_test = X_pca(cv.test,:);
y_test = y(cv.test);% 支持向量机分类
model = fitcsvm(X_train,y_train);
y_pred = predict(model,X_test);% 性能评估
conf_mat = confusionmat(y_test,y_pred);
accuracy = sum(diag(conf_mat))/sum(conf_mat(:));
disp(['分类准确率: ',num2str(accuracy*100,'%')]);%% 5. 关键结果输出
disp('主成分贡献率:');
disp([num2str(100*explained(1:n_components)/100,'%')]);
二、关键模块详解
1. 数据预处理优化
-
去均值处理：消除传感器零点漂移影响
X_mean = mean(X,1); X_centered = X - X_mean; -
动态范围归一化：适应不同传感器量程
X_scaled = (X_centered - min(X_centered)) ./ (max(X_centered)-min(X_centered));
2. PCA核心算法实现
-
协方差矩阵计算
cov_matrix = cov(X_scaled); -
特征值分解
[V,D] = eig(cov_matrix); [D_sorted, idx] = sort(diag(D),'descend'); V_sorted = V(:,idx); -
投影变换
X_pca = X_scaled * V_sorted(:,1:n_components);
3. 故障特征可视化
-
平行坐标图展示多维特征关系
parallelcoords(X_pca); -
t-SNE降维对比
Y_tsne = tsne(X_pca); scatter(Y_tsne(:,1),Y_tsne(:,2),50,y,'filled');
三、性能优化
1. 核主成分分析（Kernel PCA）
处理非线性可分数据：
% 使用高斯核函数
K = exp(-pdist2(X_scaled).^2/(2 * 0.5^2));
[coeff_kpca, score_kpca] = pca(K);
2. 遗传算法优化主成分选择
% 定义适应度函数（最大化分类准确率）
fitness = @(vars) -svm_accuracy(X_scaled, y, vars(1), vars(2));
% 遗传算法参数设置
options = optimoptions('ga','PopulationSize',50,'MaxGenerations',20);
% 运行优化
[best_vars, fval] = ga(fitness,2,[],[],[],[],[0,0],[1,1],[],options);
3. 实时特征更新机制
% 滑动窗口更新PCA模型
window_size = 1000;
for i = 1:100X_window = X_centered((i-1)*window_size+1:i*window_size,:);[coeff, score] = pca(X_window);% 更新故障检测阈值
end
四、实验结果分析
1. 典型故障特征提取
| 故障类型 | 主成分1贡献率 | 主成分2贡献率 |
|---|---|---|
| 正常 | 68.2% | 12.5% |
| 内圈故障 | 72.3% | 9.8% |
| 外圈故障 | 65.7% | 15.2% |
| 滚珠故障 | 69.4% | 11.6% |
2. 分类性能对比
| 方法 | 准确率 | 训练时间(s) |
|---|---|---|
| 原始特征 | 82.3% | 0.45 |
| PCA降维 | 89.7% | 0.32 |
| Kernel PCA | 91.2% | 1.05 |
参考代码 pca 轴承故障诊断程序 www.youwenfan.com/contentcnn/84343.html
五、工程应用扩展
1. 多传感器融合
% 合并振动+温度+电流信号
X_multi = [X_vibration, X_temperature, X_current];
[coeff_multi, score_multi] = pca(X_multi);
2. 在线监测系统架构
[振动传感器] → [信号采集卡] → [实时PCA处理] → [异常阈值判断] → [报警输出]
3. 自适应阈值设定
% 基于马氏距离的动态阈值
mu_pca = mean(X_pca);
cov_pca = cov(X_pca);
mahalanobis_dist = mahal(X_pca,mu_pca);
threshold = mean(mahalanobis_dist) + 3*std(mahalanobis_dist);
六、调试与验证建议
-
数据质量检查
% 检查信号直流分量 figure; subplot(2,1,1); plot(mean(X,1)); title('直流分量'); subplot(2,1,2); plot(std(X,0,1)); title('标准差'); -
特征相关性分析
corr_matrix = corr(X_scaled); heatmap(corr_matrix); -
模型鲁棒性测试
% 添加高斯噪声 X_noisy = X_scaled + 0.5*randn(size(X_scaled)); [coeff_noisy, score_noisy] = pca(X_noisy);