滚动轴承故障诊断实战:MATLAB包络谱分析与特征频率计算
简介一份面向滚动轴承故障诊断与MATLAB信号处理学习者的PDF讲稿以真实轴承实测数据为案例完整演示从参数确定、特征频率计算到时域分析、EMD分解及Hilbert包络谱识别的技术路线。内容先依据轴承转速、滚珠个数、滚动体直径等参数推导外圈、内圈、滚动体、保持架四类故障特征频率再提取有效值、峰值因子、峭度等时域指标做小波去噪并通过EMD分解筛选相关IMF由包络峰值164.1Hz与内圈理论频率162.21Hz接近判定内圈故障。资源共1个PDF文件大小267KB已有272人学习浏览附两个MATLAB程序分别用于小波去噪时域分析、EMD分解与包络谱绘制配有结果图与代码注释便于复现。通篇串联理论计算、信号处理与程序实现适合机械专业学生、故障诊断入门者及设备维护工程师参考。1. 滚动轴承故障诊断到底在做什么从一段振动信号到一句“该换轴承了”拿到“滚动轴承故障诊断(附MATLAB程序)讲课讲稿.pdf”这个标题时我第一反应是这应该是很多设备维护工程师和研究生都想要的那类资料——不是讲理论讲到天上去而是能拿着振动信号、打开MATLAB、跑出一个故障结论的完整流程。滚动轴承故障诊断这件事本质上就是回答三个问题轴承坏没坏、坏在哪个部件、坏到什么程度。而这三个问题的答案全部藏在振动信号的频域结构里外圈故障有外圈的特征频率内圈故障有内圈的边带滚动体故障有滚动体的冲击周期。读懂这些频率成分就相当于给轴承做了一次心电图。很多刚接触这个方向的人容易误以为故障诊断等于机器学习——采集数据、打标签、训练分类器。但实际工程中绝大多数可靠方案的第一步还是信号处理和特征提取也就是这篇文章标题里说的“附MATLAB程序”要解决的事用程序把原始振动信号转成能直接判读的谱图和指标。适合读这篇文章的人是手里有振动数据、想用MATLAB把诊断流程跑通、又不想被各类论文里花哨算法带偏的工程师和研究生。这篇文章就按我自己做项目的路径把这个流程拆开讲清楚。2. 包络谱分析为什么是轴承诊断的主力希尔伯特解调与特征频率计算2.1 轴承故障为什么会在频域里“藏起来”滚动轴承运行时的振动信号是典型的调制信号。轴承某个部件出现局部损伤后每次滚过损伤点都会产生一个冲击这个冲击的幅度很小而且被高频的共振响应“包裹”住。直接对原始信号做FFT看到的频谱往往是一片宽频噪声故障特征频率被淹没在背景里很难直接读出峰值。解决思路是把“包裹”信号的高频载波去掉只留下冲击的包络。这个操作在信号处理里叫包络解调。常见做法是先对原始信号做带通滤波把能量集中的共振频带取出来然后做希尔伯特变换得到解析信号再取解析信号的模就得到了包络信号。最后对包络信号做FFT得到的包络谱里故障特征频率及其谐波会以明显的谱峰形式出现。这套流程在MATLAB里用四五行代码就能实现但每一步的参数选择都会直接影响最终谱图的质量。包络谱分析之所以是轴承诊断的绝对主力是因为它对局部损伤类故障极其敏感。磨损类故障的振动特征不明显包络谱未必好用但点蚀、剥落、裂纹这类冲击型故障包络谱几乎一抓一个准。这也是我平时做项目时默认首选它的原因。2.2 特征频率怎么算四个部件各自对应一组公式要判断包络谱里哪个峰是故障峰必须先算出理论上故障特征频率是多少。轴承的几何参数决定特征频率常用的四个特征频率如下特征频率计算公式物理含义转频FTFT n / 60轴的旋转频率基频外圈故障特征频率BPFOBPFO (Z/2) × FT × (1 - d × cos(α) / D)滚子滚过外圈损伤点的频率内圈故障特征频率BPFIBPFI (Z/2) × FT × (1 d × cos(α) / D)滚子滚过内圈损伤点的频率滚动体故障特征频率BSFBSF (D/(2d)) × FT × (1 - (d × cos(α) / D)^2)滚动体自转一周经过损伤点的频率公式中Z是滚动体个数d是滚动体直径D是节圆直径α是接触角。这些几何参数一般在轴承型号手册里能查到但实际做诊断时最靠谱的办法是拿一个同型号的新轴承用千分尺量出滚动体直径和节圆直径再数一下滚动体个数。我见过不少翻车案例都是工程师直接抄手册里近似值导致特征频率算出来差了0.5Hz左右谱图上峰值对不上误判成没有故障。2.3 仿真信号先行先证明程序逻辑对再碰真实数据我刚接触轴承诊断时吃过一个亏直接拿现场采集的真实信号调程序谱图一团乱麻根本分不清哪个峰是特征频率、哪个是噪声随机峰值。后来养成一个习惯——先用仿真信号验证程序逻辑正确再上真实数据。仿真信号的数学模型不复杂周期性冲击信号叠加高频共振衰减和噪声。%% 参数定义 Fs 25600; % 采样率单位Hz T 10; % 信号时长单位秒 t (0:1/Fs:T-1/Fs); % 时间序列 fr 30; % 转频单位Hz对应1800rpm BPFO 123.5; % 外圈故障特征频率单位Hz fc 3500; % 共振频率单位Hz tau 0.002; % 冲击衰减时间常数单位秒 %% 生成周期性冲击序列 impulse_positions round((0:1/BPFO:T-1/BPFO) * Fs) 1; % 每隔1/BPFO秒产生一个冲击换算成采样点位置 impulse_train zeros(length(t), 1); impulse_train(impulse_positions) 1; %% 生成冲击响应指数衰减正弦波 impulse_response exp(-t/tau) .* sin(2*pi*fc*t); % 每个冲击在共振频率fc上衰减振荡 %% 卷积得到冲击响应信号 signal conv(impulse_train, impulse_response, same); %% 添加噪声 signal_noise signal 0.5 * randn(size(signal)); % 噪声幅度取冲击信号RMS的一半模拟真实工况的底噪这个代码的核心逻辑是用单位冲击序列与单冲击响应做卷积生成一段“理论外圈故障信号”。impulse_positions里每个冲击的间隔是 1/BPFO 秒这是整个仿真信号里最重要的参数——它决定了包络谱里故障峰的位置。fc和tau共同决定冲击的波形形态实际数值不需要精确只要共振频率落在带通滤波器通带内即可。噪声幅度的调整可以让信号更接近真实采集情况方便检验后面诊断程序在不同信噪比下的表现。仿真信号跑通后用来检验包络谱分析程序是否正确如果程序把 BPFO 峰值解调出来了而且谐波次数一清二楚程序逻辑就对。这个过程建议每个做轴承诊断的人都走一遍省下的排查时间远超写这段代码的时间。3. 用MATLAB跑通最小诊断流程从仿真信号到真实数据的完整代码3.1 完整的包络谱分析程序滤波、希尔伯特变换、FFT一气呵成仿真信号生成之后下一步是写包络谱分析程序。这个程序就是整个诊断流程的核心也是我平时项目里复用率最高的代码段。带通滤波的边界选择是这里最需要经验的地方直接决定包络谱的信噪比。%% 包络谱分析程序 function [env_spectrum, freq_axis] envelope_spectrum(x, Fs, bp_low, bp_high) % x : 原始振动信号列向量 % Fs : 采样率单位Hz % bp_low : 带通滤波下截止频率单位Hz % bp_high : 带通滤波上截止频率单位Hz %% 1. 带通滤波 [b, a] butter(4, [bp_low, bp_high]/(Fs/2), bandpass); % 4阶巴特沃斯带通滤波器需要归一化频率 x_filtered filtfilt(b, a, x); % 使用零相位滤波避免相位偏移影响冲击位置 %% 2. 希尔伯特变换求包络 analytical hilbert(x_filtered); % 构造解析信号实部为原信号虚部为希尔伯特变换结果 envelope abs(analytical); % 取复数的模得到包络信号 %% 3. 去除包络信号中的直流分量 envelope envelope - mean(envelope); %% 4. 对包络信号做FFT N length(envelope); freq_axis (0:N/2-1) * Fs / N; env_spectrum abs(fft(envelope)); env_spectrum env_spectrum(1:N/2); env_spectrum(2:end-1) 2 * env_spectrum(2:end-1); % 单边谱除直流外幅值加倍 %% 5. 归一化按最大峰值缩放 env_spectrum env_spectrum / max(env_spectrum); end这段程序里最需要注意的参数是bp_low和bp_high。带通范围选得不对包络谱直接失真选窄了会把边带截掉选宽了会混入无关的共振峰。我的经验是先用原始信号的频谱找能量集中的共振频带以谱峰为中心左右各取500-1000Hz作为带宽。具体到轴承诊断一般取中心频率的0.7到1.3倍作为通带边界效果比较稳定。filtfilt代替filter也是一处关键选择。filter会引入相位延迟导致冲击位置在时间域里偏移虽然包络谱的峰值位置不受影响但后续如果有做时域同步平均的需求相位失真就会累积出问题。filtfilt零相位处理虽然计算量翻倍但信号质量更好值得。3.2 真实信号处理流程数据分段、转速同步和工况记录真实采集的轴承信号和仿真信号有个本质区别真实信号里转频和故障特征频率不是恒定的。电机转速会有微小波动负载变化也会改变轴承的受力状态。所以真实数据处理的完整流程要比仿真多两步——数据分段和工况同步。数据分段的逻辑是按固定时长把长信号切成若干段每段独立做包络谱分析然后看各段特征频率处幅值的稳定性。如果某段出现明显高于其他段的峰值说明这段时间内轴承状态发生了突变。分段长度一般取1到2秒对应转频的整数倍周期这样FFT结果里转频及其谐波不会因为分段不整周期而出现频谱泄漏。转速同步是另一个容易忽视的点。真实场景里如果只有加速度信号没有同时采集转速信号那么特征频率的理论值只是个近似。我的习惯是诊断报告里标注清楚“按名义转速计算的特征频率”同时把包络谱里转频及其谐波的实测位置标出来用实际转频修正特征频率的理论值。这个做法虽然多花两分钟但能避免因为转速偏差产生的误判。3.3 程序跑通后必须看的三个输出包络谱、时域波形、特征频率标记图程序跑完不能只看一个包络谱图就下结论。我一般同时看三样东西原始时域波形、包络谱全貌、特征频率标记图。时域波形里能直接看到冲击的周期性和幅值这是判断“故障是否处于早期”最直观的证据包络谱全貌能确认特征频率及其二倍频、三倍频是否依次衰减——这是典型的局部损伤特征特征频率标记图则是在包络谱上竖线标出BPFO、BPFI、BSF的理论位置一眼就能看出谱峰和理论频率是否重合。%% 绘制诊断结果图 figure(Position, [100, 100, 1200, 800]); subplot(3, 1, 1); plot(t, signal_noise); title(原始振动信号时域波形); xlabel(时间 (s)); ylabel(加速度 (m/s^2)); subplot(3, 1, 2); [env_spectrum, freq_axis] envelope_spectrum(signal_noise, Fs, 2500, 4500); plot(freq_axis, env_spectrum); xlim([0, 500]); % 只显示低频段故障特征频率通常集中在这里 title(包络谱); xlabel(频率 (Hz)); ylabel(归一化幅值); subplot(3, 1, 3); plot(freq_axis, env_spectrum); hold on; xline(BPFO, r--, BPFO); xline(2*BPFO, g--, 2×BPFO); xline(3*BPFO, b--, 3×BPFO); xlim([0, 500]); title(包络谱与特征频率标记); xlabel(频率 (Hz)); ylabel(归一化幅值); legend(包络谱, BPFO, 2×BPFO, 3×BPFO);这里xline的作用是在谱图上画理论特征频率的标记线实际特征峰与标记线重合度越高诊断结论越可靠。我见过有人只看包络谱里最高峰就下结论结果最高峰其实是转频的二倍频跟轴承故障一点关系没有。有了标记图这种低级误判能直接避免。4. 从“看谱图”到“自动判病”时域特征指标与阈值诊断的MATLAB实现4.1 为什么需要时域指标谱图判断有主观性指标判断可复现包络谱分析需要人来看图判定这在实际工程里有三个问题一是看图判定需要经验新手容易把边带当成故障峰二是人工判定无法批量处理连续监测数据三是不同人看图结论可能不一致。所以真正落地到设备状态监测系统里一定需要自动化的判定逻辑。时域特征指标就是这类逻辑的基础用几个数值描述信号特征再用阈值判断状态。常用的时域指标包括均方根值RMS、峰值因子Crest Factor、峭度Kurtosis。RMS反映振动能量总体水平磨损类故障和润滑不良时RMS会持续爬升峰值因子反映冲击的尖锐程度早期点蚀时峰值因子升高明显峭度对冲击型故障极其敏感正常轴承峭度接近3出现局部损伤后峭度可以涨到5甚至更高。三个指标各有盲区组合起来使用才能覆盖不同故障阶段。4.2 指标计算代码一个函数算出全套特征%% 时域特征指标计算函数 function features bearing_time_features(x) % x: 原始振动信号列向量 features.RMS sqrt(mean(x.^2)); % 均方根值衡量振动能量 features.peak max(abs(x)); % 绝对峰值 features.crest_factor features.peak / features.RMS; % 峰值因子反映冲击尖锐程度 features.kurtosis mean((x - mean(x)).^4) / (std(x)^4); % 峭度高斯分布时为3冲击型故障时显著升高 features.skewness mean((x - mean(x)).^3) / (std(x)^3); % 偏度反映波形不对称性 features.impulse_factor features.peak / mean(abs(x)); % 脉冲因子对早期冲击敏感 features.clearance_factor features.peak / (mean(sqrt(abs(x))))^2; % 裕度因子对磨损类故障敏感 end这段程序里的kurtosis计算用的是定义式不是MATLAB自带的kurtosis函数两者对“是否减3”的处理不同。定义式算出的正态信号峭度是3自带函数算出的正态信号峭度是0。如果你跟同行对比数据必须先确认双方用的哪种定义否则数值对不上会闹笑话。crest_factor对早期故障最敏感但也最容易受随机噪声干扰单次计算不稳定建议做多次平均。这些指标在MATLAB里都有现成函数但自己写一遍的好处是清楚每个值的量纲和物理含义后续改进算法时改起来更快。4.3 阈值怎么定从“固定阈值”到“相对阈值”的工程经验指标算出来后判定逻辑的核心是阈值。最粗糙的做法是拍脑袋定一个固定值——比如峭度大于4就报警。这在固定工况、单一型号轴承上可能够用但换了转速、换了负载、换了轴承型号阈值立刻失效。实际工程里更可靠的思路是“相对阈值”法先采集设备正常运行状态下的振动信号建立基线指标值然后用实时指标与基线值做比值。比值超过1.5倍提示关注超过2倍提示报警超过3倍提示停机检查。这个方法的优点是不需要知道轴承的绝对振动水平每台设备用自己“健康状态”做参照天然的工况适应能力。相对阈值的具体实现不复杂把设备首次运行时采集的数据存一份基线之后每次采集数据计算指标后除以基线值。要注意的是基线数据必须覆盖设备不同负载工况我见过有工程师只采集了空载状态的基线带载一开就误报不断。采集基线时至少记录三个负载点下的指标平均值作为查询表保存。5. 故障诊断实战避坑采样参数、频谱泄露与误判的5条排查记录5.1 采样率和采样时长不匹配特征频率根本“看不见”现象包络谱做完后BPFO附近是一片平缓的隆起没有尖锐谱峰。原因采样率或者FFT点数不足频率分辨率太差特征频率和相邻干扰频率无法分离。解决先算频率分辨率 Δf Fs / N保证 Δf 小于特征频率间距的五分之一。比如 BPFO 是123.5Hz转频是30Hz边带间距就是30Hz这时 Δf 需要小于6Hz也就是说 N/Fs 要大于0.167秒采样时长必须足够长。这个坑我踩过最狠的一次用的现有数据文件只有0.1秒时长做出来的包络谱里 BPFO 和它的边带糊成一团当时还以为是滤波参数不对调了一整天滤波器最后才发现是数据长度不够白白浪费时间。5.2 带通滤波范围选错包络谱里全是假峰值现象包络谱里出现很多不明来源的尖锐谱峰且峰值位置随滤波边界变化而变化。原因带通范围包含了多个共振频带包络里混入了不同频率的调制成分。解决先看原始信号频谱图确定能量集中的主共振带只保留这一个频带做包络解调。如果频谱上有两个明显的共振峰优先选择能量更大的那个如果两个能量接近分别做两次包络谱分析对比哪个结果里故障特征更清晰。5.3 频谱泄露让特征频率峰值“长胖”窄带搜索自动判错现象用峰值搜索算法在指定频率附近找最大值时找到的频率和理论BPFO偏差很大。原因信号做了加窗后发现峰值位置偏移或者没做整周期截断导致频谱能量分散。解决做FFT前先用detrend去除趋势项再选择汉宁窗或布莱克曼窗在峰值搜索时搜索范围放宽到理论值±3%以内取窗口内能量重心作为特征频率实测值而不是直接取最大幅值对应的频率。能量重心法在MATLAB里实现就三行循环逻辑但非常实用先定位最大峰然后在峰值左右各取几个谱线按幅值加权求平均频率。这个方法抗频谱泄漏能力比直接取峰值强很多。5.4 转速不稳定导致边带消失误判为“轴承没问题”现象设备带载运行时包络谱里只有BPFO基频二倍频和三倍频几乎没有幅值也不稳定。原因转速波动导致冲击周期不恒定包络谱里能量被分散到多个相邻谱线。解决采集数据的同时记录转速信号如果转速波动超过±1%改用角度域分析方法也就是把等时间间隔采样重采样为等角度间隔再做包络谱分析。MATLAB里可以用resample配合转速脉冲信号实现阶次跟踪重采样后的阶次谱不受转速波动影响。这个坑在做变频设备时几乎必踩。现场电机转速从30Hz变到50Hz如果还按固定时间采样分析结果完全对不上。加一路转速通道做阶次跟踪能省掉后续大量返工。5.5 阈值诊断只看单一指标早期故障被噪声数据“骗过”现象峭度指标忽高忽低有时超过报警阈值触发停机有时又回落到正常区间。原因单次时域计算的峭度对方波干扰和随机脉冲极为敏感一个外界的电磁干扰尖峰就能让峭度翻倍。解决判定逻辑改为“连续三段数据中至少两段指标超阈值才报警”同时把RMS指标纳入联合判定——RMS未超限而峭度超限时提示“关注”两个指标同时超限才触发报警。此外对原始数据先做中值滤波或剔除离群点能有效减少单点干扰的影响。我现在做诊断程序时报警逻辑一定写进“参数配置模块”方便现场调整报警阈值和确认窗口长度。这个模块看起来不起眼实际用起来比算法本身还重要——现场工程师每天面对的是不同工况没有可调的阈值界面再好的算法也会因为误报太多被关掉。6. 进阶用法把包络谱、特征指标和趋势曲线串成一册“轴承体检报告”诊断流程跑通后真正让方案产生长期价值的是趋势跟踪。单次诊断只能回答“此刻有没有故障”而趋势跟踪能回答“故障在加速还是稳定预计还能跑多久”。我的做法是将每日采集的数据做包络谱分析和特征指标计算把RMS、峭度、BPFO处幅值三个量分别存储绘制成趋势曲线。RMS曲线反映整体劣化过程峭度曲线反映冲击强度的短期波动BPFO幅值反映故障源的能量变化三张曲线图叠加工况记录设备状态一目了然。趋势曲线里最值得看的是变化率而不是绝对值。某些轴承的RMS可能在三个月里从2缓慢升到3而另一些轴承两周内从2暴涨到5。前者是均匀磨损后者是快速扩展的疲劳剥落。变化率的阈值设定同样用相对法——当前值除以前一周移动平均值比值超过1.8就该安排停机检修窗口了。最后一个习惯是每次诊断完成后把谱图、时域波形、指标数据和工况信息一起归档到以设备编号和日期命名的文件夹里方便后续追溯。我见过太多项目因为现场人员流动导致之前积累的诊断经验全部丢失新来的人要从零摸索阈值和滤波参数。加上这套归档习惯后哪怕换了操作者也能按历史报告快速对齐参数。做轴承诊断这几年我最大的教训是别迷信复杂的算法先把手上的信号处理基本功做扎实。一套简单的包络谱加三个时域指标在绝大多数工业场景里已经能解决八成的诊断需求。把这段MATLAB程序和归档习惯固化下来比反复尝试各种新模型实在得多。希望帮到你。本文还有配套的精品资源点击获取