从MATLAB代码到故障诊断：手把手教你分析风机CMS振动数据（附完整脚本）

工业风机振动数据诊断实战从MATLAB代码到故障特征解析在工业设备健康管理领域振动数据分析是识别潜在机械故障的黄金标准。特别是对于风力发电机组这类高价值资产提前发现轴承、齿轮箱等关键部件的异常状态能够避免数百万的意外停机损失。本文将带您完整走通CMS振动数据分析的全流程——从原始CSV数据导入到专业频谱图生成再到故障特征识别与诊断建议形成。1. 数据准备与环境配置1.1 数据文件结构规范工业现场采集的CMS振动数据通常以CSV格式存储规范的目录结构能大幅提升分析效率。建议按以下结构组织数据项目根目录/ ├── raw_data/ # 原始振动数据 │ ├── F020_0417_0424_6H/ # 测点编号日期范围 │ │ ├── 20240417_0900.csv │ │ └── 20240418_1430.csv ├── scripts/ # MATLAB分析脚本 └── results/ # 输出图表关键参数说明采样频率(fs)必须与采集设备设置一致示例中为25600Hz数据长度通常截取稳定运行段的1-2秒数据对应25600-51200个采样点测点命名如6H表示发电机驱动端轴承径向测点1.2 MATLAB环境初始化% 初始化脚本 CMS_Init.m clear; clc; close all; set(0,defaultfigurecolor,w); % 白色背景 set(0,DefaultAxesFontSize,12); % 统一字体大小 format compact; % 紧凑显示格式 % 添加工具包路径 addpath(genpath(./lib)); % 自定义函数目录 addpath(genpath(./data)); % 数据目录注意使用hua_fft和hua_baol等自定义频谱分析函数前需确保其已加入MATLAB路径。2. 数据加载与预处理2.1 批量读取CSV数据function [data_cell, file_info] loadCMSData(folder_path) % 获取目录下所有CSV文件 file_list dir(fullfile(folder_path, *.csv)); % 预分配存储空间 data_cell cell(length(file_list), 1); file_info struct(name,{}, date,{}, size,{}); % 并行读取提升效率 parfor i 1:length(file_list) file_path fullfile(folder_path, file_list(i).name); raw_data csvread(file_path); % 提取文件名中的时间信息假设格式为YYYYMMDD_HHMM [~, fname] fileparts(file_list(i).name); datetime_str extractBetween(fname, 1, 13); % 存储元数据 file_info(i).name file_list(i).name; file_info(i).date datetime(datetime_str,... InputFormat,yyyyMMdd_HHmm); file_info(i).size length(raw_data); % 存储振动数据去除前0.1秒瞬态过程 data_cell{i} raw_data(0.1*fs:end); end end2.2 数据质量检查执行以下检查确保数据可靠性采样完整性验证expected_points 10 * fs; % 假设每文件应含10秒数据 missing_files find([file_info.size] 0.9*expected_points);异常值检测function is_clean checkDataQuality(signal, fs) % 检查直流偏移 dc_offset mean(signal); % 检查峰值是否超量程 pk_value max(abs(signal)); % 检查噪声基底 noise_floor median(abs(signal))/0.6745; is_clean (dc_offset0.1) (pk_value9.5) (noise_floor0.001); end趋势项消除clean_signal detrend(raw_signal, linear);3. 时频分析技术实现3.1 时域波形可视化function plotTimeDomain(y, fs, plot_title) N length(y); time_axis (0:N-1)/fs; figure(Position, [100 100 800 400]); plot(time_axis, y, b, LineWidth, 1.2); % 专业图表格式设置 xlim([0 min(1, N/fs)]); % 默认显示1秒或全部 xlabel(Time (s), FontWeight,bold); ylabel(Amplitude (m/s²), FontWeight,bold); title([Time Domain: plot_title], Interpreter,none); grid on; % 添加统计信息标注 stats_str sprintf(Peak: %.3f\nRMS: %.3f\nKurtosis: %.2f,... max(abs(y)), rms(y), kurtosis(y)); annotation(textbox,[0.75 0.7 0.2 0.2],... String,stats_str,FitBoxToText,on); end3.2 快速傅里叶变换(FFT)实现function [freq, amp] calcFFT(y, fs, n_avg) % 分段平均计算提高信噪比 seg_length floor(length(y)/n_avg); window hann(seg_length); % 初始化存储 Pxx zeros(seg_length/21, n_avg); for k 1:n_avg segment y((k-1)*seg_length1 : k*seg_length); % 加窗处理 windowed segment .* window; % FFT计算 Y fft(windowed); Pxx(:,k) abs(Y(1:seg_length/21)); end % 平均处理 amp mean(Pxx, 2) * 2/sum(window); freq linspace(0, fs/2, seg_length/21); end关键参数说明表参数推荐值作用说明窗函数Hanning减少频谱泄漏平均次数8-16提高信噪比频率分辨率≥1Hz确保能分离邻近频率3.3 包络分析技术function [env, env_freq] calcEnvelope(y, fs, bp_freq) % 带通滤波 [b,a] butter(4, bp_freq/(fs/2), bandpass); filtered filtfilt(b,a,y); % 希尔伯特变换提取包络 analytic hilbert(filtered); env abs(analytic); % 对包络信号做FFT [env_freq, env_amp] calcFFT(env, fs, 8); % 绘制包络谱 figure; semilogy(env_freq, env_amp, r, LineWidth,1.5); xlabel(Frequency (Hz)); ylabel(Envelope Amplitude); title(Envelope Spectrum); grid on; end提示轴承故障特征频率通常出现在100-2000Hz范围内建议带通滤波设置在此区间。4. 故障特征识别与诊断4.1 轴承特征频率计算根据轴承几何参数计算理论故障频率function [bpfo, bpfi, bsf, ftf] calcBearingFreq(rpm, D, d, n) % rpm: 轴转速(rpm) % D: 轴承节径(mm) % d: 滚动体直径(mm) % n: 滚动体数量 RPS rpm / 60; % 转每秒 % 计算公式 bpfo n/2 * RPS * (1 - d/D * cos(0)); % 外圈故障频率 bpfi n/2 * RPS * (1 d/D * cos(0)); % 内圈故障频率 bsf D/d * RPS * (1 - (d/D)^2 * cos(0)^2); % 滚动体故障频率 ftf 1/2 * RPS * (1 - d/D * cos(0)); % 保持架故障频率 end4.2 典型故障频谱特征常见故障类型与频谱特征对照表故障类型时域特征频谱特征包络谱特征外圈损伤周期性冲击边频带间隔bpfo清晰bpfo谐波内圈损伤幅值调制边频带间隔bpfi转频与bpfi组合滚动体损伤随机冲击宽带噪声提升bsf及其谐波电腐蚀高频振荡高频离散峰电源频率相关4.3 案例电腐蚀故障诊断针对示例中361Hz间隔的频带特征匹配361Hz接近50Hz电源频率的7倍频7×50350Hz包络谱中出现转频12倍谐波12×29.68356Hz诊断结论diagnosis_result struct(... FaultType, Electrical Erosion,... Evidence, {361Hz spacing in spectrum,... 12X RPM in envelope},... Confidence, 0.85,... Maintenance, {Check grounding brushes,... Monitor bearing temperature}... );趋势监控建议建立361Hz成分的幅值趋势图设置三级报警阈值注意/警告/危险5. 自动化报告生成5.1 结果可视化整合function generateReport(file_info, diagnosis, fig_handles) % 创建PDF报告 import mlreportgen.dom.*; report Document(CMS_Report, pdf); open(report); % 添加标题页 title_page TitlePage; title_page.Title CMS Vibration Analysis Report; title_page.Author Diagnosis System; add(report, title_page); % 添加数据概览 section Section(Data Summary); table Table(4); table.Style {Border(solid), Width(100%)}; table_entry {Filename, file_info.name;... Recording Time, datestr(file_info.date);... Data Length, [num2str(file_info.size/fs) s];... Sampling Rate, [num2str(fs) Hz]}; for i 1:4 row TableRow; append(row, TableEntry(table_entry{i,1})); append(row, TableEntry(table_entry{i,2})); append(row, TableEntry(table_entry{i,3})); append(row, TableEntry(table_entry{i,4})); append(table, row); end append(section, table); add(report, section); % 插入图表 for i 1:length(fig_handles) fig Figure(fig_handles(i)); fig.Snapshot.Caption [Analysis Result num2str(i)]; add(report, fig); end % 添加诊断结论 section Section(Diagnosis Conclusion); ul UnorderedList; for i 1:length(diagnosis.Evidence) append(ul, ListItem(diagnosis.Evidence{i})); end append(section, ul); add(report, section); close(report); end5.2 历史数据对比分析function trendAnalysis(project_folder, target_freq) % 获取所有历史分析结果 result_files dir(fullfile(project_folder, results, *.mat)); % 提取目标频率成分历史值 trend_data zeros(length(result_files), 3); for i 1:length(result_files) data load(fullfile(result_files(i).folder, result_files(i).name)); [~, idx] min(abs(data.freq - target_freq)); trend_data(i,:) [datenum(data.file_info.date),... data.amp(idx),... data.env_amp(round(target_freq))]; end % 绘制趋势图 figure; subplot(2,1,1); plot(datetime(trend_data(:,1),ConvertFrom,datenum),... trend_data(:,2), o-); title([Spectral Component at num2str(target_freq) Hz]); subplot(2,1,2); plot(datetime(trend_data(:,1),ConvertFrom,datenum),... trend_data(:,3), s-r); title([Envelope Component at num2str(target_freq) Hz]); end