代码是一个用于生成工业轴承故障振动信号的综合性仿真系统,主要目的是为轴承故障诊断算法开发和测试提供高质量、真实感强的训练和验证数据。
系统通过物理建模和随机过程相结合的方式,模拟了包括健康状态、8种单一故障和3种混合故障在内的多种轴承运行状态。整个算法流程首先进行系统初始化,包括参数配置、随机种子设置和输出目录创建;然后根据用户配置构建故障类型列表,包括健康状态、单一故障和混合故障;接着进入主生成循环,对每种故障类型生成指定数量的信号样本。
对于每个信号样本,系统首先生成基础振动信号,然后根据故障类型注入相应的故障特征,这些故障特征基于旋转机械动力学原理,如不对中故障表现为2倍频和3倍频谐波,不平衡故障表现为1倍频主导且与速度平方相关,轴承游隙故障产生次同步振动等。系统还考虑了多种实际运行条件的影响,包括转速变化、负载波动和温度变化,并基于这些条件物理正确地计算Sommerfeld数来反映轴承润滑状态。
为了增强信号的真实性,系统引入了8种不同类型的噪声源,包括测量噪声、电磁干扰、粉红噪声、环境漂移、量化噪声、传感器漂移、混叠伪影和脉冲噪声。此外,系统还支持多严重度级别的故障模拟,包括初期、轻度、中度和严重四个等级,并且部分信号会模拟故障的时序演进过程。
系统还包含了瞬态行为模拟,如速度斜坡、负载阶跃和热膨胀等非平稳过程。最后,系统通过数据增强技术(包括时间偏移、幅度缩放和噪声注入)来扩展数据集多样性,并为每个信号保存详细的元数据信息,包括故障类型、严重度、运行条件、物理参数和噪声配置等,确保生成的数据集既具有物理真实性又具有统计多样性,可直接用于机器学习模型的训练和评估。
算法步骤如下:
1.系统初始化阶段
设置信号生成的基本参数,包括采样频率、信号持续时间、基础旋转速度
配置数据集生成参数,确定每种故障类型要生成的信号数量
建立输出目录结构用于存储生成的数据文件
设置随机数种子以确保结果的可重现性
2. 故障类型配置阶段
选择要包含的健康状态、单一故障类型和混合故障组合
构建完整的故障类型列表,包括8种单一故障和3种混合故障
配置故障严重度分级系统,分为初期、轻度、中度和严重四个等级
设置故障严重度的时间演进特性,模拟故障随时间的渐进发展
3. 运行条件模拟阶段
生成变化的运行速度,在标称速度基础上添加随机波动
模拟负载变化,在额定负载的30%到100%范围内随机选择
设置工作温度变化,在40°C到80°C范围内随机分布
计算综合运行因子,将速度、负载和温度的影响进行组合
4. 物理参数计算阶段
基于运行条件物理正确地计算轴承的Sommerfeld数
考虑粘度随温度的变化关系、速度因子和负载因子的影响
生成Reynolds数和轴承间隙比等关键物理参数
计算物理因子来反映物理参数对振动信号的整体影响
5. 瞬态行为建模阶段
模拟速度斜坡变化,在信号中间段实现速度的线性增加或减少
生成负载阶跃变化,在特定时间点突然改变负载水平
模拟热膨胀效应,展现温度变化引起的缓慢振动特性变化
应用瞬态调制函数,将瞬态行为与基础信号进行融合
6. 基础信号生成阶段
创建基础振动信号,包含低水平的高斯白噪声
根据运行条件和物理参数调整基础信号的幅度水平
确保基础信号具有真实工业振动信号的统计特性
7. 噪声模型应用阶段
添加传感器电子噪声,模拟测量系统的不完美性
引入电磁干扰,模拟电源线和其他电气设备的干扰
生成粉红噪声,反映1/f频率特性在真实信号中的存在
添加环境漂移,模拟温度和环境变化引起的缓慢偏移
应用量化噪声,模拟模数转换过程中的精度限制
加入传感器漂移,反映传感器性能随时间的变化
模拟混叠伪影,展现采样不足时的高频信号折叠
生成脉冲噪声,模拟偶发的冲击和干扰事件
8. 故障特征注入阶段
对于不对中故障,注入强烈的2倍频和3倍频谐波分量
对于不平衡故障,生成以1倍频为主的振动且与速度平方相关
对于轴承游隙故障,产生次同步振动并伴随谐波成分
对于润滑问题,模拟粘滑运动现象并添加金属接触冲击
对于气蚀故障,生成随机出现的高频能量突发信号
对于磨损故障,产生宽带噪声和幅度调制效应
对于油涡动故障,形成稳定的次同步振动成分
对于混合故障,将相关单一故障的特征进行叠加组合
9. 数据增强处理阶段
应用时间偏移增强,通过循环移位改变信号的起始位置
进行幅度缩放增强,随机调整信号的总体幅度水平
实施噪声注入增强,添加额外噪声以提高模型鲁棒性
确保增强后的信号仍保持原有的物理特性和统计分布
10. 元数据组装阶段
记录故障类型、严重度级别和演进特性信息
保存运行条件参数,包括速度、负载、温度等具体数值
存储物理参数计算结果,如Sommerfeld数、Reynolds数等
记录瞬态行为类型和相关参数配置
保存信号基本属性,包括RMS值、峰值、峰值因子等统计特征
记录应用的噪声源类型和增强处理方法
添加生成时间戳、版本信息和随机种子等系统元数据
11. 数据保存阶段
按照故障类型和信号序号组织文件名
区分基础信号和增强信号的文件命名
将信号数据、采样频率和完整元数据保存为MAT文件格式
确保数据文件的结构化和标准化,便于后续处理分析
12. 质量评估阶段
统计总生成信号数量和各类故障分布
计算生成效率,评估系统性能表现
验证信号质量,确保符合预期的物理真实性和统计特性
生成配置总结报告,确认所有功能模块的正确运行
13. 结果输出阶段
显示详细的生成统计信息
提供数据集质量评估和预期应用效果
输出完整的配置总结和使用建议
确保生成的数据集可直接用于故障诊断算法的训练和测试
流程图可以适当参考
开始 ├─ 系统初始化 │ ├─ 加载配置参数 │ ├─ 设置随机种子 │ └─ 创建输出目录 │ ├─ 故障类型配置 │ ├─ 健康状态 │ ├─ 单一故障(8种) │ └─ 混合故障(3种) │ ├─ 主生成循环(按故障类型) │ │ │ ├─ 信号生成循环(按信号数量) │ │ │ │ │ ├─ 严重度配置 │ │ │ ├─ 选择严重度级别 │ │ │ └─ 时序演进设置 │ │ │ │ │ ├─ 运行条件设置 │ │ │ ├─ 速度变化 │ │ │ ├─ 负载波动 │ │ │ └─ 温度影响 │ │ │ │ │ ├─ 物理参数计算 │ │ │ ├─ Sommerfeld数 │ │ │ ├─ Reynolds数 │ │ │ └─ 间隙比 │ │ │ │ │ ├─ 瞬态行为模拟 │ │ │ ├─ 速度斜坡 │ │ │ ├─ 负载阶跃 │ │ │ └─ 热膨胀 │ │ │ │ │ ├─ 基础信号生成 │ │ │ └─ 低水平噪声 │ │ │ │ │ ├─ 噪声模型应用 │ │ │ ├─ 测量噪声 │ │ │ ├─ 电磁干扰 │ │ │ ├─ 粉红噪声 │ │ │ ├─ 环境漂移 │ │ │ ├─ 量化噪声 │ │ │ ├─ 传感器漂移 │ │ │ ├─ 混叠伪影 │ │ │ └─ 脉冲噪声 │ │ │ │ │ ├─ 故障特征注入 │ │ │ ├─ 基于故障类型 │ │ │ ├─ 物理正确建模 │ │ │ └─ 严重度调整 │ │ │ │ │ ├─ 数据增强处理 │ │ │ ├─ 时间偏移 │ │ │ ├─ 幅度缩放 │ │ │ └─ 噪声注入 │ │ │ │ │ ├─ 元数据组装 │ │ │ └─ 完整信号信息 │ │ │ │ │ └─ 信号数据保存 │ │ └─ 文件输出 │ │ │ └─ 进度显示 │ ├─ 数据增强处理 │ └─ 额外样本生成 │ └─ 生成总结报告 ├─ 统计信息 ├─ 性能指标 └─ 质量评估部分代码如下
function plot_single_signal_analysis(signal, fs, fault_type, metadata, config) % 单信号综合分析绘图 t = (0:length(signal)-1)' / fs; % 创建图形 figure('Position', [100, 100, 1600, 1200], 'Color', 'white'); % 1. 时域波形 subplot(3, 3, 1); plot(t, signal, 'b-', 'LineWidth', config.style.line_width); xlabel('时间 (s)'); ylabel('振幅'); title(sprintf('%s - 时域波形', fault_type), 'FontSize', config.style.font_size+2, 'FontWeight', 'bold'); grid on; xlim([0, min(1.0, max(t))]); % 添加统计信息 stats_text = sprintf('RMS: %.4f\n峰值: %.4f\n峰值因子: %.2f', ... rms(signal), max(abs(signal)), max(abs(signal))/(rms(signal)+eps)); text(0.02, 0.98, stats_text, 'Units', 'normalized', ... 'VerticalAlignment', 'top', 'FontSize', config.style.font_size-2, ... 'BackgroundColor', 'white', 'EdgeColor', 'black'); % 2. 概率密度函数 subplot(3, 3, 2); [pdf_values, pdf_edges] = histcounts(signal, 100, 'Normalization', 'pdf'); pdf_centers = (pdf_edges(1:end-1) + pdf_edges(2:end)) / 2; plot(pdf_centers, pdf_values, 'r-', 'LineWidth', config.style.line_width); xlabel('振幅'); ylabel('概率密度'); title('概率密度函数 (PDF)'); grid on; % 3. 自相关函数 subplot(3, 3, 3); [acf, lags] = xcorr(signal, 200, 'coeff'); lag_time = lags / fs; plot(lag_time, acf, 'g-', 'LineWidth', config.style.line_width); xlabel('时延 (s)'); ylabel('自相关系数'); title('自相关函数'); grid on; xlim([-0.02, 0.02]); % 4. 功率谱密度 subplot(3, 3, 4); [psd, freq_psd] = pwelch(signal, hanning(1024), 512, 1024, fs); semilogy(freq_psd, psd, 'm-', 'LineWidth', config.style.line_width); xlabel('频率 (Hz)'); ylabel('功率谱密度'); title('功率谱密度 (PSD)'); grid on; xlim(config.freq_range); % 5. 频谱图 subplot(3, 3, 5); [s, f, t_spec] = spectrogram(signal, 256, 250, 512, fs, 'yaxis'); imagesc(t_spec, f, 10*log10(abs(s))); axis xy; xlabel('时间 (s)'); ylabel('频率 (Hz)'); title('频谱图'); colorbar; ylim(config.spectrogram_freq_range); % 6. 包络谱 subplot(3, 3, 6); analytic_signal = hilbert(signal); envelope_signal = abs(analytic_signal); [env_psd, env_freq] = pwelch(envelope_signal, hanning(1024), 512, 1024, fs); plot(env_freq, env_psd, 'c-', 'LineWidth', config.style.line_width); xlabel('频率 (Hz)'); ylabel('包络谱幅度'); title('包络谱'); grid on; xlim([0, 200]); % 7. 时域统计特征 subplot(3, 3, 7); features = calculate_time_domain_features(signal); feature_names = {'RMS', '峰值', '峰度', '偏度', '波形因子', '脉冲因子'}; bar(features); set(gca, 'XTickLabel', feature_names); ylabel('数值'); title('时域统计特征'); grid on; rotateXLabels(gca, 45); % 8. 运行条件信息 subplot(3, 3, 8); if isfield(metadata, 'operating_factor') op_params = [metadata.speed_rpm, metadata.load_percent, ... metadata.temperature_C, metadata.operating_factor]; param_names = {'转速 (RPM)', '负载 (%)', '温度 (°C)', '运行因子'}; bar(op_params); set(gca, 'XTickLabel', param_names); ylabel('数值'); title('运行条件'); grid on; rotateXLabels(gca, 45); else text(0.5, 0.5, '无运行条件数据', 'HorizontalAlignment', 'center', ... 'FontSize', config.style.font_size); title('运行条件'); end % 9. 物理参数信息 subplot(3, 3, 9); if isfield(metadata, 'sommerfeld_number') physics_params = [metadata.sommerfeld_number, metadata.reynolds_number, ... metadata.clearance_ratio, metadata.physics_factor]; physics_names = {'Sommerfeld数', 'Reynolds数', '间隙比', '物理因子'}; bar(physics_params); set(gca, 'XTickLabel', physics_names); ylabel('数值'); title('物理参数'); grid on; rotateXLabels(gca, 45); else text(0.5, 0.5, '无物理参数数据', 'HorizontalAlignment', 'center', ... 'FontSize', config.style.font_size); title('物理参数'); end % 添加总体标题 sgtitle(sprintf('轴承故障诊断分析 - %s (严重度: %s)', ... fault_type, metadata.severity), 'FontSize', 16, 'FontWeight', 'bold'); % 保存图形 if config.save_figures filename = fullfile(config.output_dir, ... sprintf('detailed_analysis_%s.png', fault_type)); saveas(gcf, filename, config.figure_format); fprintf('保存详细分析图: %s\n', filename); end end function plot_comparative_analysis(signals, fault_types, fs, config) % 多故障类型对比分析 figure('Position', [100, 100, 1800, 1000], 'Color', 'white'); % 1. 时域波形对比 subplot(2, 3, 1); hold on; colors = config.style.colors; for i = 1:length(signals) t = (0:length(signals{i})-1)' / fs; plot(t(1:min(2000, length(t))), signals{i}(1:min(2000, length(signals{i}))), ... 'Color', colors(i,:), 'LineWidth', config.style.line_width, ... 'DisplayName', fault_types{i}); end xlabel('时间 (s)'); ylabel('振幅'); title('时域波形对比'); legend('show', 'Location', 'best'); grid on; xlim([0, 0.1]); % 2. 功率谱密度对比 subplot(2, 3, 2); hold on; for i = 1:length(signals) [psd, freq] = pwelch(signals{i}, hanning(1024), 512, 1024, fs); plot(freq, 10*log10(psd), 'Color', colors(i,:), ... 'LineWidth', config.style.line_width, 'DisplayName', fault_types{i}); end xlabel('频率 (Hz)'); ylabel('功率谱密度 (dB)'); title('功率谱密度对比'); legend('show', 'Location', 'best'); grid on; xlim(config.freq_range); % 3. 特征雷达图 subplot(2, 3, 3); features_matrix = []; feature_names = {'RMS', '峰值', '峰度', '偏度', '波形因子', '脉冲因子'}; for i = 1:length(signals) features = calculate_time_domain_features(signals{i}); features_matrix(i, :) = features; end % 归一化特征 features_norm = zeros(size(features_matrix)); for j = 1:size(features_matrix, 2) features_norm(:, j) = features_matrix(:, j) / max(features_matrix(:, j)); end % 绘制雷达图 polarplot_radar(features_norm', feature_names, fault_types, colors); title('特征雷达图对比'); % 4. 包络谱对比 subplot(2, 3, 4); hold on; for i = 1:length(signals) analytic_signal = hilbert(signals{i}); envelope_signal = abs(analytic_signal); [env_psd, env_freq] = pwelch(envelope_signal, hanning(1024), 512, 1024, fs); plot(env_freq, 10*log10(env_psd), 'Color', colors(i,:), ... 'LineWidth', config.style.line_width, 'DisplayName', fault_types{i}); end xlabel('频率 (Hz)'); ylabel('包络谱 (dB)'); title('包络谱对比'); legend('show', 'Location', 'best'); grid on; xlim([0, 200]); % 5. 概率密度函数对比 subplot(2, 3, 5); hold on; for i = 1:length(signals) [pdf_values, pdf_edges] = histcounts(signals{i}, 50, 'Normalization', 'pdf'); pdf_centers = (pdf_edges(1:end-1) + pdf_edges(2:end)) / 2; plot(pdf_centers, pdf_values, 'Color', colors(i,:), ... 'LineWidth', config.style.line_width, 'DisplayName', fault_types{i}); end xlabel('振幅'); ylabel('概率密度'); title('概率密度函数对比'); legend('show', 'Location', 'best'); grid on; % 6. 特征散点图 subplot(2, 3, 6); rms_vals = cellfun(@(x) rms(x), signals); kurtosis_vals = cellfun(@(x) kurtosis(x), signals); scatter(rms_vals, kurtosis_vals, 100, colors(1:length(signals), :), 'filled'); % 添加标签 for i = 1:length(signals) text(rms_vals(i), kurtosis_vals(i), fault_types{i}, ... 'FontSize', config.style.font_size-2, 'HorizontalAlignment', 'center'); end xlabel('RMS值'); ylabel('峰度'); title('特征散点图 (RMS vs 峰度)'); grid on; sgtitle('轴承故障多类型对比分析', 'FontSize', 16, 'FontWeight', 'bold'); % 保存图形 if config.save_figures filename = fullfile(config.output_dir, 'comparative_analysis.png'); saveas(gcf, filename, config.figure_format); fprintf('保存对比分析图: %s\n', filename); end end参考文章:
基于多物理场耦合与动态演进建模的工业轴承故障振动信号生产级仿真生成(MATLAB) - 哥廷根数学学派的文章
https://zhuanlan.zhihu.com/p/1978074463636562216
工学博士,担任《Mechanical System and Signal Processing》审稿专家,担任
《中国电机工程学报》优秀审稿专家,《控制与决策》,《系统工程与电子技术》,《电力系统保护与控制》,《宇航学报》等EI期刊审稿专家。
擅长领域:现代信号处理,机器学习,深度学习,数字孪生,时间序列分析,设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。