MATLAB实现滚动轴承故障诊断（外圈故障）

一、滚动轴承故障诊断原理

滚动轴承故障诊断是通过分析振动信号来识别轴承运行状态的技术。当轴承出现故障时，会产生周期性冲击信号，其特征频率与故障类型和轴承几何参数相关。

1. 外圈故障特征频率

外圈故障特征频率公式：

$f_{out}=\frac{n}{2}(1−\frac{d}{D}cosα)f_r$

其中：

• n：滚动体数量
• d：滚动体直径
• D：节圆直径
• α：接触角
• $f_r$：轴旋转频率

2. 诊断流程

1. 信号采集：获取轴承振动信号
2. 预处理：去噪、滤波
3. 特征提取：时域、频域、时频域特征
4. 故障识别：模式分类（正常/外圈故障）

二、MATLAB实现步骤

1. 参数设置与信号模拟

clear; clc; close all; % ==================== 轴承参数设置 ==================== n = 8; % 滚动体数量 d = 7.94e-3; % 滚动体直径(m) D = 39.04e-3; % 节圆直径(m) alpha = 0; % 接触角(度) fr = 25; % 轴旋转频率(Hz) - 1500 RPM fs = 10000; % 采样频率(Hz) T = 1; % 信号时长(s) % 计算外圈故障特征频率 f_out = (n/2)*(1 - d/D*cosd(alpha))*fr; % ==================== 模拟振动信号 ==================== t = 0:1/fs:T-1/fs; % 时间向量 N = length(t); % 正常轴承信号（基频振动 + 噪声） x_normal = 0.5*sin(2*pi*30*t) + 0.2*randn(1, N); % 外圈故障信号（基频振动 + 周期性冲击 + 噪声） impulse_period = 1/f_out; % 冲击周期 impulse_times = 0:impulse_period:T-impulse_period; x_fault = 0.5*sin(2*pi*30*t); % 基频成分 % 添加周期性冲击 for i = 1:length(impulse_times) idx = round(impulse_times(i)*fs); if idx > 0 && idx <= N % 冲击波形（高斯脉冲） pulse_width = 0.005; % 脉冲宽度(s) pulse_samples = round(pulse_width*fs); t_pulse = linspace(-pulse_width/2, pulse_width/2, pulse_samples); pulse = exp(-50*(t_pulse).^2); % 高斯脉冲 % 添加冲击 start_idx = max(1, idx-floor(pulse_samples/2)); end_idx = min(N, idx+floor(pulse_samples/2)-1); seg_len = end_idx - start_idx + 1; x_fault(start_idx:end_idx) = x_fault(start_idx:end_idx) + 0.8*pulse(1:seg_len); end end x_fault = x_fault + 0.2*randn(1, N); % 添加噪声 % 可视化信号 figure; subplot(2,1,1); plot(t, x_normal); title('正常轴承振动信号'); xlabel('时间(s)'); ylabel('振幅'); xlim([0, 0.2]); subplot(2,1,2); plot(t, x_fault); title('外圈故障轴承振动信号'); xlabel('时间(s)'); ylabel('振幅'); xlim([0, 0.2]); 

2. 信号预处理（小波去噪）

% ==================== 小波去噪 ==================== % 正常信号去噪 [cA_normal, cD_normal] = dwt(x_normal, 'db4'); thresh_normal = thselect(cD_normal, 'sqtwolog'); % 阈值选择 cD_normal_thresh = wthresh(cD_normal, 's', thresh_normal); x_normal_denoised = idwt(cA_normal, cD_normal_thresh, 'db4'); % 故障信号去噪 [cA_fault, cD_fault] = dwt(x_fault, 'db4'); thresh_fault = thselect(cD_fault, 'sqtwolog'); cD_fault_thresh = wthresh(cD_fault, 's', thresh_fault); x_fault_denoised = idwt(cA_fault, cD_fault_thresh, 'db4'); % 可视化去噪效果 figure; subplot(2,1,1); plot(t, x_normal, 'b', t, x_normal_denoised, 'r'); legend('原始信号', '去噪信号'); title('正常信号去噪效果'); xlim([0, 0.1]); subplot(2,1,2); plot(t, x_fault, 'b', t, x_fault_denoised, 'r'); legend('原始信号', '去噪信号'); title('故障信号去噪效果'); xlim([0, 0.1]); 

3. 特征提取

(1) 时域特征

% ==================== 时域特征提取 ==================== features = @(x) [rms(x), % 均方根 peak2peak(x), % 峰峰值 kurtosis(x), % 峭度 skewness(x), % 偏度 crestfactor(x)]; % 峰值因子 % 计算正常信号特征 feat_normal = features(x_normal_denoised); % 计算故障信号特征 feat_fault = features(x_fault_denoised); % 显示特征 fprintf('正常信号时域特征:\n'); fprintf('RMS: %.4f, Peak2Peak: %.4f, Kurtosis: %.4f, Skewness: %.4f, CrestFactor: %.4f\n', feat_normal); fprintf('\n故障信号时域特征:\n'); fprintf('RMS: %.4f, Peak2Peak: %.4f, Kurtosis: %.4f, Skewness: %.4f, CrestFactor: %.4f\n', feat_fault); 

(2) 频域特征（FFT分析）

% ==================== 频域特征提取 ==================== NFFT = 2^nextpow2(N); % FFT点数 f = fs/2*linspace(0,1,NFFT/2+1); % 频率向量 % 正常信号频谱 Y_normal = fft(x_normal_denoised, NFFT)/N; P_normal = 2*abs(Y_normal(1:NFFT/2+1)); % 故障信号频谱 Y_fault = fft(x_fault_denoised, NFFT)/N; P_fault = 2*abs(Y_fault(1:NFFT/2+1)); % 计算特征频率处的幅值 [~, idx_out] = min(abs(f - f_out)); % 外圈故障特征频率索引 feature_freq_amplitude = [P_normal(idx_out), P_fault(idx_out)]; % 可视化频谱 figure; subplot(2,1,1); plot(f, P_normal); title('正常信号频谱'); xlabel('频率(Hz)'); ylabel('幅值'); xlim([0, 500]); subplot(2,1,2); plot(f, P_fault); hold on; plot([f_out, f_out], [0, max(P_fault)], 'r--'); % 标记故障特征频率 title('故障信号频谱'); xlabel('频率(Hz)'); ylabel('幅值'); legend('频谱', '外圈故障特征频率'); xlim([0, 500]); 

(3) 包络分析（解调）

% ==================== 包络分析 ==================== % 希尔伯特变换提取包络 env_normal = abs(hilbert(x_normal_denoised)); env_fault = abs(hilbert(x_fault_denoised)); % 包络谱分析 Y_env_normal = fft(env_normal, NFFT)/N; P_env_normal = 2*abs(Y_env_normal(1:NFFT/2+1)); Y_env_fault = fft(env_fault, NFFT)/N; P_env_fault = 2*abs(Y_env_fault(1:NFFT/2+1)); % 可视化包络谱 figure; subplot(2,1,1); plot(f, P_env_normal); title('正常信号包络谱'); xlabel('频率(Hz)'); ylabel('幅值'); xlim([0, 500]); subplot(2,1,2); plot(f, P_env_fault); hold on; plot([f_out, f_out], [0, max(P_env_fault)], 'r--'); % 标记故障特征频率 title('故障信号包络谱'); xlabel('频率(Hz)'); ylabel('幅值'); legend('包络谱', '外圈故障特征频率'); xlim([0, 500]); 

4. 故障诊断（模式识别）

(1) 特征矩阵构建

% ==================== 构建特征矩阵 ==================== % 生成更多样本（实际应用中应从实验获取） num_samples = 50; % 每种状态样本数 % 正常样本特征 features_normal = zeros(num_samples, 6); for i = 1:num_samples % 添加随机噪声模拟不同工况 noise_level = 0.1 + 0.1*rand(); x = x_normal_denoised + noise_level*randn(1, N); features_normal(i, 1:5) = features(x); features_normal(i, 6) = feature_freq_amplitude(1) * (0.8 + 0.4*rand()); end % 故障样本特征 features_fault = zeros(num_samples, 6); for i = 1:num_samples noise_level = 0.1 + 0.1*rand(); x = x_fault_denoised + noise_level*randn(1, N); features_fault(i, 1:5) = features(x); features_fault(i, 6) = feature_freq_amplitude(2) * (0.8 + 0.4*rand()); end % 合并特征矩阵 features_all = [features_normal; features_fault]; labels = [zeros(num_samples, 1); ones(num_samples, 1)]; % 0=正常, 1=故障 

(2) SVM分类器训练与测试

% ==================== SVM分类器 ==================== % 划分训练集和测试集 rng(42); % 设置随机种子 cv = cvpartition(size(features_all, 1), 'HoldOut', 0.3); trainIdx = training(cv); testIdx = test(cv); X_train = features_all(trainIdx, :); y_train = labels(trainIdx); X_test = features_all(testIdx, :); y_test = labels(testIdx); % 训练SVM分类器 svm_model = fitcsvm(X_train, y_train, 'KernelFunction', 'rbf', 'BoxConstraint', 1); % 测试分类器 y_pred = predict(svm_model, X_test); % 计算准确率 accuracy = sum(y_pred == y_test) / numel(y_test); conf_mat = confusionmat(y_test, y_pred); fprintf('分类准确率: %.2f%%\n', accuracy*100); disp('混淆矩阵:'); disp(conf_mat); 

(3) 诊断结果可视化

% ==================== 可视化结果 ==================== % 特征空间可视化（前两维） figure; gscatter(features_all(:,1), features_all(:,2), labels); hold on; plot(svm_model.SupportVectors(:,1), svm_model.SupportVectors(:,2), 'ko', 'MarkerSize', 10, 'LineWidth', 2); title('特征空间分布 (RMS vs 峰峰值)'); xlabel('均方根(RMS)'); ylabel('峰峰值'); legend('正常', '故障', '支持向量'); grid on; % ROC曲线 [y_score, scores] = predict(svm_model, X_test); [X_roc, Y_roc, T_roc, AUC] = perfcurve(y_test, scores(:,2), 1); figure; plot(X_roc, Y_roc, 'LineWidth', 2); title(sprintf('ROC曲线 (AUC = %.2f)', AUC)); xlabel('假阳性率'); ylabel('真阳性率'); grid on; 

三、实际应用与扩展

1. 真实数据采集建议

% ==================== 真实数据采集 ==================== % 使用加速度传感器采集振动信号 % 示例代码（需硬件支持）： % % % 创建数据采集对象 % daqreset; % s = daq.createSession('ni'); % addAnalogInputChannel(s, 'Dev1', 0, 'Voltage'); % s.Rate = fs; % s.DurationInSeconds = T; % % % 采集正常轴承信号 % data_normal = s.startForeground(); % % % 采集外圈故障轴承信号（人工模拟） % data_fault = s.startForeground(); 

2. 高级特征提取方法

% ==================== 小波包能量特征 ==================== % 小波包分解 wpt = wpdec(x_fault_denoised, 4, 'db4'); % 计算各节点能量 nodes = [5, 6, 9, 10, 13, 14, 17, 18]; % 选择节点 energy = zeros(1, length(nodes)); for i = 1:length(nodes) node_data = wprcoef(wpt, nodes(i)); energy(i) = sum(node_data.^2); end % 能量占比作为特征 energy_ratio = energy / sum(energy); 

3. 深度学习诊断方法

% ==================== CNN故障诊断 ==================== % 构建CNN模型 layers = [ imageInputLayer([1 1024 1]) convolution2dLayer([1 64], 16, 'Padding', 'same') batchNormalizationLayer reluLayer maxPooling2dLayer([1 4], 'Stride', [1 4]) convolution2dLayer([1 32], 32, 'Padding', 'same') batchNormalizationLayer reluLayer maxPooling2dLayer([1 4], 'Stride', [1 4]) fullyConnectedLayer(64) reluLayer fullyConnectedLayer(2) softmaxLayer classificationLayer ]; % 训练选项 options = trainingOptions('adam', ... 'MaxEpochs', 30, ... 'MiniBatchSize', 32, ... 'ValidationData', { X_val, y_val}, ... 'Plots', 'training-progress'); % 训练网络 net = trainNetwork(X_train, y_train, layers, options); 

四、诊断结果分析与工程应用

1. 典型诊断结果

• 时域特征：故障信号峭度值显著高于正常信号（>3 vs <3）
• 频域特征：故障信号在外圈特征频率处有明显谱峰
• 包络分析：故障信号包络谱在特征频率处出现明显峰值
• 分类准确率：SVM分类器在模拟数据上可达95%以上准确率

2. 工程应用建议

1. 传感器布置：在轴承座径向安装加速度传感器
2. 采样参数：采样频率≥10倍最高分析频率（通常5-20kHz）
3. 诊断流程： 实时采集振动信号 预处理（去噪、滤波） 特征提取（时域+频域+时频域） 模式识别（阈值判断/SVM/神经网络）
4. 预警阈值： 峭度>5 → 潜在故障 特征频率幅值>3倍基线 → 明确故障

3. 扩展方向

• 复合故障诊断：同时识别内圈、外圈、滚动体故障
• 变工况诊断：考虑转速变化对特征频率的影响
• 迁移学习：利用预训练模型解决小样本问题
• 边缘计算：部署轻量级模型到嵌入式设备

参考代码 实现滚动轴承故障诊断 www.youwenfan.com/contentalg/83047.html

五、总结

本MATLAB实现展示了滚动轴承外圈故障诊断的完整流程：

1. 信号模拟：生成正常和外圈故障振动信号
2. 预处理：小波去噪提高信噪比
3. 特征提取：时域（峭度、峰值因子）、频域（特征频率幅值）、时频域（包络谱）
4. 故障识别：SVM分类器实现高精度诊断

通过调整轴承参数和故障严重程度，可模拟不同工况下的故障特征。实际应用中，建议使用凯斯西储大学轴承数据中心（CWRU Bearing Data Center）的真实数据验证算法性能。

注意：人工模拟信号仅用于算法验证，实际诊断需使用真实采集的振动信号。工业应用中建议结合温度、油液分析等多传感器信息进行综合诊断。