2026/9/11 7:05:33
FEATURED · 精选文章 
1. 项目概述轴承故障信号模拟与分析的工程价值作为一名在设备状态监测领域工作多年的工程师我深知滚动轴承故障诊断的重要性。轴承作为旋转机械的核心部件其健康状况直接影响设备运行安全。传统故障诊断依赖物理实验成本高、周期长而MATLAB提供的信号模拟与分析方法让我们能在计算机上高效完成故障特征研究。这个项目使用MATLAB R2018a实现了两大核心功能滚动轴承典型故障内圈、外圈、滚动体损伤的振动信号模拟生成基于时频联合分析方法的故障特征提取特别说明虽然示例使用R2018a版本但文中所有代码和方法均兼容MATLAB 2016b及以上版本部分新版本特性会额外标注。2. 故障信号建模原理与实现2.1 轴承故障的动力学模型滚动轴承故障振动信号具有明显的调制特性主要包含以下成分承载区周期性冲击故障特征频率高频共振轴承结构固有频率幅值调制旋转频率对冲击的调制数学表达式可表示为x(t) ∑A_i·s(t - iT - τ_i)·cos(2πf_n(t - iT) φ_i) n(t)其中A_i第i次冲击的幅值T故障特征周期1/故障频率f_n共振频率通常2-10kHzn(t)背景噪声2.2 MATLAB实现代码解析function [t, vibration] generateBearingFault(... faultType, fs, duration, rpm, bearingParams) % 参数说明 % faultType: inner/outer/ball 对应不同故障类型 % fs: 采样频率(Hz) % duration: 信号时长(s) % rpm: 轴转速(rpm) % bearingParams: 结构体包含轴承几何参数 t 0:1/fs:duration; N length(t); % 计算特征频率 fc getFaultFrequency(faultType, rpm, bearingParams); % 生成冲击序列 T 1/fc; impulseTimes 0:T:duration; impulseTrain zeros(size(t)); for i 1:length(impulseTimes) idx round(impulseTimes(i)*fs) 1; if idx N % 指数衰减的冲击模型 decay exp(-2000*(t - impulseTimes(i)).^2); impulseTrain impulseTrain decay.*(t impulseTimes(i)); end end % 添加共振调制 fn 4500; % 典型轴承共振频率 carrier sin(2*pi*fn*t); % 添加幅值调制反映旋转频率 fm rpm/60; amplitudeMod 1 0.3*sin(2*pi*fm*t); % 合成信号并添加噪声 vibration amplitudeMod .* impulseTrain .* carrier; vibration vibration 0.1*randn(size(vibration)); end关键技巧实际工程中共振频率fn需要通过轴承型号查询或实验测定。常用轴承的共振频率范围通常在2-10kHz之间。2.3 不同故障类型的特征频率计算故障特征频率与轴承几何参数强相关function fc getFaultFrequency(faultType, rpm, bp) % bp应包含以下字段 % bp.d: 滚动体直径 % bp.D: 节圆直径 % bp.n: 滚动体数量 % bp.contactAngle: 接触角(度) omega rpm/60; % 转/秒 alpha bp.contactAngle * pi/180; switch faultType case inner fc bp.n/2 * omega * (1 bp.d/bp.D * cos(alpha)); case outer fc bp.n/2 * omega * (1 - bp.d/bp.D * cos(alpha)); case ball fc bp.D/(2*bp.d) * omega * (1 - (bp.d/bp.D * cos(alpha))^2); end end3. 时频分析方法深度解析3.1 短时傅里叶变换(STFT)实现STFT是分析非平稳信号的经典方法通过加窗分段计算频谱function [s, f, t] mySTFT(x, fs, window, noverlap, nfft) % 参数说明 % window: 窗函数向量 % noverlap: 重叠样本数 % nfft: FFT点数 wlen length(window); step wlen - noverlap; frames fix((length(x) - noverlap)/step); % 预分配矩阵 s zeros(nfft/21, frames); % 逐帧计算 for i 1:frames idx (i-1)*step 1; xw x(idx:idxwlen-1) .* window; X abs(fft(xw, nfft)).^2; s(:,i) X(1:nfft/21); end f (0:nfft/2)*fs/nfft; t ((0:frames-1)*step wlen/2)/fs; end工程经验对于轴承故障分析推荐使用汉宁窗窗长选择2-3个故障周期重叠率75%可获得最佳时频分辨率平衡。3.2 时频图优化显示技巧% 生成时频图 [s, f, t] mySTFT(vibration, fs, hann(1024), 768, 1024); % 对数缩放增强细节 s_log 10*log10(s eps); % 智能频率范围截取 max_freq 10000; % 根据实际轴承特性调整 f_idx f max_freq; % 绘制时频图 figure imagesc(t, f(f_idx), s_log(f_idx,:)) axis xy colormap jet colorbar xlabel(Time (s)) ylabel(Frequency (Hz)) title(STFT Time-Frequency Analysis) % 叠加特征频率参考线 hold on plot(t, repmat(fc, size(t)), w--, LineWidth, 1.5) plot(t, repmat(fcfn, size(t)), g--, LineWidth, 1.5)3.3 其他时频分析方法对比方法优点缺点适用场景STFT计算简单实现直观时频分辨率固定快速分析故障初步筛查小波变换多分辨率分析适合瞬态特征基函数选择依赖经验冲击特征明显的早期故障HHT自适应分解适合非线性信号端点效应严重计算量大复杂调制信号分析Wigner-Ville时频分辨率高存在交叉项干扰高精度分析信噪比较高时4. 工程应用中的关键问题与解决方案4.1 噪声环境下的特征增强实际工业现场信噪比(SNR)往往低于实验室条件。我们采用包络分析增强特征function envelope hilbertEnvelope(x) % 希尔伯特变换提取包络 analytic hilbert(x); envelope abs(analytic); % 低通滤波平滑包络 [b,a] butter(4, 0.1, low); envelope filtfilt(b, a, envelope); end应用示例env hilbertEnvelope(vibration); env env - mean(env); [s_env, ~, ~] mySTFT(env, fs, hann(512), 384, 512);4.2 多故障耦合情况分析当多个故障同时存在时可采用以下策略带通滤波分离不同共振频带对每个频带单独进行包络分析比较各频带包络谱中的特征频率成分% 设计带通滤波器组 bpFreqs [2000 4000; 3500 5500; 5000 7500]; % 示例频带 figure for i 1:size(bpFreqs,1) [b,a] butter(4, bpFreqs(i,:)/(fs/2), bandpass); x_band filtfilt(b, a, vibration); env_band hilbertEnvelope(x_band); % 计算包络谱 [Penv, fenv] pwelch(env_band, hann(1024), 512, 1024, fs); subplot(size(bpFreqs,1),1,i) plot(fenv, 10*log10(Penv)) xlim([0 1000]) title(sprintf(Envelope Spectrum %.0f-%.0f Hz, bpFreqs(i,:))) end4.3 实际工程中的参数选择指南采样频率最低要求≥2.56 × 最高感兴趣频率推荐值通常12.8-25.6kHz覆盖常见轴承共振频率分析时长至少包含10个完整的故障周期示例计算对于100Hz故障频率至少0.1s数据STFT参数% 自适应参数设置示例 fc 85; % 预估故障频率(Hz) T 1/fc; window_len round(3 * T * fs); % 覆盖3个故障周期 window hann(window_len); noverlap round(0.75 * window_len); nfft max(1024, 2^nextpow2(window_len));5. 完整工作流示例5.1 从参数设置到结果可视化的完整流程%% 1. 轴承参数设置 bearingParams.d 7.94e-3; % 滚动体直径(m) bearingParams.D 39.04e-3; % 节圆直径(m) bearingParams.n 9; % 滚动体数量 bearingParams.contactAngle 15; % 接触角(度) %% 2. 生成故障信号 fs 25600; % 采样频率(Hz) duration 0.5; % 信号时长(s) rpm 1800; % 轴转速(rpm) faultType outer; % 故障类型 [t, vibration] generateBearingFault(... faultType, fs, duration, rpm, bearingParams); %% 3. 时频分析 window hann(2048); noverlap 1536; nfft 4096; [s, f, t_stft] mySTFT(vibration, fs, window, noverlap, nfft); s_log 10*log10(s eps); %% 4. 包络分析 env hilbertEnvelope(vibration); [Penv, fenv] pwelch(env, hann(1024), 512, 1024, fs); %% 5. 可视化 figure(Position, [100 100 900 600]) % 时域波形 subplot(3,1,1) plot(t, vibration) xlabel(Time (s)) ylabel(Amplitude) title(Time Domain Waveform) xlim([0 0.1]) % 时频分析 subplot(3,1,2) imagesc(t_stft, f(f10000), s_log(f10000,:)) axis xy colormap jet colorbar xlabel(Time (s)) ylabel(Frequency (Hz)) title(STFT Time-Frequency Analysis) % 包络谱 subplot(3,1,3) plot(fenv, 10*log10(Penv)) xlabel(Frequency (Hz)) ylabel(Power/frequency (dB/Hz)) title(Envelope Spectrum) xlim([0 1000]) grid on % 标记理论故障频率 fc getFaultFrequency(faultType, rpm, bearingParams); hold on plot([fc fc], ylim, r--) text(fc, max(ylim)-5, sprintf(%.1f Hz, fc), Color, r)5.2 结果解读要点时域波形观察周期性冲击的存在性评估信噪比水平时频图查找与理论故障频率一致的垂直条纹注意高频区域共振频带的能量变化包络谱确认故障频率及其谐波成分比较各频率成分的相对幅值诊断决策标准当包络谱中故障频率成分比相邻频率成分高6dB以上时可判定存在对应故障。6. 工程实践中的经验总结信号采集注意事项加速度传感器应尽量靠近轴承安装确保采样频率设置正确避免混叠记录完整的工况参数转速、负载等分析优化技巧% 使用重采样技术分析变速工况 original_fs 25600; resample_ratio rpm/nominal_rpm; vibration_resampled resample(vibration, original_fs, round(original_fs*resample_ratio));常见问题排查问题时频图中看不到明显故障特征检查确认传感器安装是否松动调整尝试不同共振频带的包络分析问题包络谱中出现未知频率成分检查核对设备其他旋转部件的特征频率验证改变转速观察频率是否线性变化性能优化建议对于长时间信号分析可采用分段处理block_size 60 * fs; % 60秒每块 for i 0:ceil(length(vibration)/block_size)-1 idx i*block_size (1:min(block_size, length(vibration)-i*block_size)); processBlock(vibration(idx)); end扩展应用方向结合机器学习进行自动故障分类开发实时监测系统需转为C/C代码集成到SCADA系统中进行在线监测
MATLAB轴承故障信号模拟与时频分析技术详解
1. 项目概述轴承故障信号模拟与分析的工程价值作为一名在设备状态监测领域工作多年的工程师我深知滚动轴承故障诊断的重要性。轴承作为旋转机械的核心部件其健康状况直接影响设备运行安全。传统故障诊断依赖物理实验成本高、周期长而MATLAB提供的信号模拟与分析方法让我们能在计算机上高效完成故障特征研究。这个项目使用MATLAB R2018a实现了两大核心功能滚动轴承典型故障内圈、外圈、滚动体损伤的振动信号模拟生成基于时频联合分析方法的故障特征提取特别说明虽然示例使用R2018a版本但文中所有代码和方法均兼容MATLAB 2016b及以上版本部分新版本特性会额外标注。2. 故障信号建模原理与实现2.1 轴承故障的动力学模型滚动轴承故障振动信号具有明显的调制特性主要包含以下成分承载区周期性冲击故障特征频率高频共振轴承结构固有频率幅值调制旋转频率对冲击的调制数学表达式可表示为x(t) ∑A_i·s(t - iT - τ_i)·cos(2πf_n(t - iT) φ_i) n(t)其中A_i第i次冲击的幅值T故障特征周期1/故障频率f_n共振频率通常2-10kHzn(t)背景噪声2.2 MATLAB实现代码解析function [t, vibration] generateBearingFault(... faultType, fs, duration, rpm, bearingParams) % 参数说明 % faultType: inner/outer/ball 对应不同故障类型 % fs: 采样频率(Hz) % duration: 信号时长(s) % rpm: 轴转速(rpm) % bearingParams: 结构体包含轴承几何参数 t 0:1/fs:duration; N length(t); % 计算特征频率 fc getFaultFrequency(faultType, rpm, bearingParams); % 生成冲击序列 T 1/fc; impulseTimes 0:T:duration; impulseTrain zeros(size(t)); for i 1:length(impulseTimes) idx round(impulseTimes(i)*fs) 1; if idx N % 指数衰减的冲击模型 decay exp(-2000*(t - impulseTimes(i)).^2); impulseTrain impulseTrain decay.*(t impulseTimes(i)); end end % 添加共振调制 fn 4500; % 典型轴承共振频率 carrier sin(2*pi*fn*t); % 添加幅值调制反映旋转频率 fm rpm/60; amplitudeMod 1 0.3*sin(2*pi*fm*t); % 合成信号并添加噪声 vibration amplitudeMod .* impulseTrain .* carrier; vibration vibration 0.1*randn(size(vibration)); end关键技巧实际工程中共振频率fn需要通过轴承型号查询或实验测定。常用轴承的共振频率范围通常在2-10kHz之间。2.3 不同故障类型的特征频率计算故障特征频率与轴承几何参数强相关function fc getFaultFrequency(faultType, rpm, bp) % bp应包含以下字段 % bp.d: 滚动体直径 % bp.D: 节圆直径 % bp.n: 滚动体数量 % bp.contactAngle: 接触角(度) omega rpm/60; % 转/秒 alpha bp.contactAngle * pi/180; switch faultType case inner fc bp.n/2 * omega * (1 bp.d/bp.D * cos(alpha)); case outer fc bp.n/2 * omega * (1 - bp.d/bp.D * cos(alpha)); case ball fc bp.D/(2*bp.d) * omega * (1 - (bp.d/bp.D * cos(alpha))^2); end end3. 时频分析方法深度解析3.1 短时傅里叶变换(STFT)实现STFT是分析非平稳信号的经典方法通过加窗分段计算频谱function [s, f, t] mySTFT(x, fs, window, noverlap, nfft) % 参数说明 % window: 窗函数向量 % noverlap: 重叠样本数 % nfft: FFT点数 wlen length(window); step wlen - noverlap; frames fix((length(x) - noverlap)/step); % 预分配矩阵 s zeros(nfft/21, frames); % 逐帧计算 for i 1:frames idx (i-1)*step 1; xw x(idx:idxwlen-1) .* window; X abs(fft(xw, nfft)).^2; s(:,i) X(1:nfft/21); end f (0:nfft/2)*fs/nfft; t ((0:frames-1)*step wlen/2)/fs; end工程经验对于轴承故障分析推荐使用汉宁窗窗长选择2-3个故障周期重叠率75%可获得最佳时频分辨率平衡。3.2 时频图优化显示技巧% 生成时频图 [s, f, t] mySTFT(vibration, fs, hann(1024), 768, 1024); % 对数缩放增强细节 s_log 10*log10(s eps); % 智能频率范围截取 max_freq 10000; % 根据实际轴承特性调整 f_idx f max_freq; % 绘制时频图 figure imagesc(t, f(f_idx), s_log(f_idx,:)) axis xy colormap jet colorbar xlabel(Time (s)) ylabel(Frequency (Hz)) title(STFT Time-Frequency Analysis) % 叠加特征频率参考线 hold on plot(t, repmat(fc, size(t)), w--, LineWidth, 1.5) plot(t, repmat(fcfn, size(t)), g--, LineWidth, 1.5)3.3 其他时频分析方法对比方法优点缺点适用场景STFT计算简单实现直观时频分辨率固定快速分析故障初步筛查小波变换多分辨率分析适合瞬态特征基函数选择依赖经验冲击特征明显的早期故障HHT自适应分解适合非线性信号端点效应严重计算量大复杂调制信号分析Wigner-Ville时频分辨率高存在交叉项干扰高精度分析信噪比较高时4. 工程应用中的关键问题与解决方案4.1 噪声环境下的特征增强实际工业现场信噪比(SNR)往往低于实验室条件。我们采用包络分析增强特征function envelope hilbertEnvelope(x) % 希尔伯特变换提取包络 analytic hilbert(x); envelope abs(analytic); % 低通滤波平滑包络 [b,a] butter(4, 0.1, low); envelope filtfilt(b, a, envelope); end应用示例env hilbertEnvelope(vibration); env env - mean(env); [s_env, ~, ~] mySTFT(env, fs, hann(512), 384, 512);4.2 多故障耦合情况分析当多个故障同时存在时可采用以下策略带通滤波分离不同共振频带对每个频带单独进行包络分析比较各频带包络谱中的特征频率成分% 设计带通滤波器组 bpFreqs [2000 4000; 3500 5500; 5000 7500]; % 示例频带 figure for i 1:size(bpFreqs,1) [b,a] butter(4, bpFreqs(i,:)/(fs/2), bandpass); x_band filtfilt(b, a, vibration); env_band hilbertEnvelope(x_band); % 计算包络谱 [Penv, fenv] pwelch(env_band, hann(1024), 512, 1024, fs); subplot(size(bpFreqs,1),1,i) plot(fenv, 10*log10(Penv)) xlim([0 1000]) title(sprintf(Envelope Spectrum %.0f-%.0f Hz, bpFreqs(i,:))) end4.3 实际工程中的参数选择指南采样频率最低要求≥2.56 × 最高感兴趣频率推荐值通常12.8-25.6kHz覆盖常见轴承共振频率分析时长至少包含10个完整的故障周期示例计算对于100Hz故障频率至少0.1s数据STFT参数% 自适应参数设置示例 fc 85; % 预估故障频率(Hz) T 1/fc; window_len round(3 * T * fs); % 覆盖3个故障周期 window hann(window_len); noverlap round(0.75 * window_len); nfft max(1024, 2^nextpow2(window_len));5. 完整工作流示例5.1 从参数设置到结果可视化的完整流程%% 1. 轴承参数设置 bearingParams.d 7.94e-3; % 滚动体直径(m) bearingParams.D 39.04e-3; % 节圆直径(m) bearingParams.n 9; % 滚动体数量 bearingParams.contactAngle 15; % 接触角(度) %% 2. 生成故障信号 fs 25600; % 采样频率(Hz) duration 0.5; % 信号时长(s) rpm 1800; % 轴转速(rpm) faultType outer; % 故障类型 [t, vibration] generateBearingFault(... faultType, fs, duration, rpm, bearingParams); %% 3. 时频分析 window hann(2048); noverlap 1536; nfft 4096; [s, f, t_stft] mySTFT(vibration, fs, window, noverlap, nfft); s_log 10*log10(s eps); %% 4. 包络分析 env hilbertEnvelope(vibration); [Penv, fenv] pwelch(env, hann(1024), 512, 1024, fs); %% 5. 可视化 figure(Position, [100 100 900 600]) % 时域波形 subplot(3,1,1) plot(t, vibration) xlabel(Time (s)) ylabel(Amplitude) title(Time Domain Waveform) xlim([0 0.1]) % 时频分析 subplot(3,1,2) imagesc(t_stft, f(f10000), s_log(f10000,:)) axis xy colormap jet colorbar xlabel(Time (s)) ylabel(Frequency (Hz)) title(STFT Time-Frequency Analysis) % 包络谱 subplot(3,1,3) plot(fenv, 10*log10(Penv)) xlabel(Frequency (Hz)) ylabel(Power/frequency (dB/Hz)) title(Envelope Spectrum) xlim([0 1000]) grid on % 标记理论故障频率 fc getFaultFrequency(faultType, rpm, bearingParams); hold on plot([fc fc], ylim, r--) text(fc, max(ylim)-5, sprintf(%.1f Hz, fc), Color, r)5.2 结果解读要点时域波形观察周期性冲击的存在性评估信噪比水平时频图查找与理论故障频率一致的垂直条纹注意高频区域共振频带的能量变化包络谱确认故障频率及其谐波成分比较各频率成分的相对幅值诊断决策标准当包络谱中故障频率成分比相邻频率成分高6dB以上时可判定存在对应故障。6. 工程实践中的经验总结信号采集注意事项加速度传感器应尽量靠近轴承安装确保采样频率设置正确避免混叠记录完整的工况参数转速、负载等分析优化技巧% 使用重采样技术分析变速工况 original_fs 25600; resample_ratio rpm/nominal_rpm; vibration_resampled resample(vibration, original_fs, round(original_fs*resample_ratio));常见问题排查问题时频图中看不到明显故障特征检查确认传感器安装是否松动调整尝试不同共振频带的包络分析问题包络谱中出现未知频率成分检查核对设备其他旋转部件的特征频率验证改变转速观察频率是否线性变化性能优化建议对于长时间信号分析可采用分段处理block_size 60 * fs; % 60秒每块 for i 0:ceil(length(vibration)/block_size)-1 idx i*block_size (1:min(block_size, length(vibration)-i*block_size)); processBlock(vibration(idx)); end扩展应用方向结合机器学习进行自动故障分类开发实时监测系统需转为C/C代码集成到SCADA系统中进行在线监测 
