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 ❤️）
• 频域特征：故障信号在外圈特征频率处有明显谱峰
• 包络分析：故障信号包络谱在特征频率处出现明显峰值
• 分类准确率：SVM分类器在模拟数据上可达95%以上准确率

2. 工程应用建议

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

3. 扩展方向

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

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

五、总结

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

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

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

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