MATLAB实现滚动轴承故障诊断（外圈故障）-阿里云开发者社区

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

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

1. 外圈故障特征频率

外圈故障特征频率公式：

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

其中：

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

2. 诊断流程

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

二、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. 工程应用建议

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

3. 扩展方向

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

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

五、总结

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

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

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

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

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

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

1. 外圈故障特征频率

2. 诊断流程

二、MATLAB实现步骤

1. 参数设置与信号模拟

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

3. 特征提取

(1) 时域特征

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

(3) 包络分析（解调）

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

(1) 特征矩阵构建

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

(3) 诊断结果可视化

三、实际应用与扩展

1. 真实数据采集建议

2. 高级特征提取方法

3. 深度学习诊断方法

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

1. 典型诊断结果

2. 工程应用建议

3. 扩展方向

五、总结

热门文章

最新文章

相关电子书

探索云世界

热门

云计算

大数据

云原生

人工智能

数据库

开发与运维

活动广场

任务中心

训练营

直播

乘风者计划

下载

镜像站

技术资料

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

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

1. 外圈故障特征频率

2. 诊断流程

二、MATLAB实现步骤

1. 参数设置与信号模拟

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

3. 特征提取

(1) 时域特征

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

(3) 包络分析（解调）

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

(1) 特征矩阵构建

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

(3) 诊断结果可视化

三、实际应用与扩展

1. 真实数据采集建议

2. 高级特征提取方法

3. 深度学习诊断方法

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

1. 典型诊断结果

2. 工程应用建议

3. 扩展方向

五、总结

热门文章

最新文章

相关电子书