从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 end

2.2 数据质量检查

执行以下检查确保数据可靠性：

1. 采样完整性验证：

   expected_points = 10 * fs; % 假设每文件应含10秒数据 missing_files = find([file_info.size] < 0.9*expected_points);

2. 异常值检测：

   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_offset<0.1) & (pk_value<9.5) & (noise_floor>0.001); end

3. 趋势项消除：

   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'); end

3.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/2+1, n_avg); for k = 1:n_avg segment = y((k-1)*seg_length+1 : k*seg_length); % 加窗处理 windowed = segment .* window; % FFT计算 Y = fft(windowed); Pxx(:,k) = abs(Y(1:seg_length/2+1)); end % 平均处理 amp = mean(Pxx, 2) * 2/sum(window); freq = linspace(0, fs/2, seg_length/2+1)'; end

关键参数说明表：

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)); % 保持架故障频率 end

4.2 典型故障频谱特征

常见故障类型与频谱特征对照表：

4.3 案例：电腐蚀故障诊断

针对示例中361Hz间隔的频带：

1. 特征匹配：

   • 361Hz接近50Hz电源频率的7倍频（7×50=350Hz）
   • 包络谱中出现转频12倍谐波（12×29.68=356Hz）

2. 诊断结论：

   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'}... );

3. 趋势监控建议：

   • 建立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); end

5.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