尧图网站建设 
简介这套电机轴承故障诊断程序基于MATLAB实现采用经验模态分解EMD方法处理非线性、非平稳振动信号面向机械故障诊断研究人员、电机维护工程师及信号处理相关专业学生。包内包含8个文件以7个.m脚本和1个.mat数据文件为主脚本涵盖EMD分解、希尔伯特包络谱、能量谱计算以及内圈/外圈故障分析等完整流程数据文件可直接用于算法验证。整套程序压缩后仅1.14MB结构简洁便于理解与应用。已有247人学习下载实用性强。通过学习这份代码使用者能掌握从振动信号采集到EMD分解、IMF分量提取、包络谱与能量谱构建、故障特征频率识别的全链路方法并根据频谱峰值判断轴承健康状态为滚动体缺陷、内外圈故障等典型失效模式提供诊断思路。同时该案例也是非线性信号处理在工程中落地的良好范本适合作为课题研究或课程设计的参考实现。1. 电机轴承故障诊断为什么要先用 EMD 分解电机轴承出现早期缺陷时传感器采到的振动信号并不是干净的周期冲击。转频及其谐波、电机电气干扰、现场噪声全混在一起故障冲击又会激起轴承座和传感器安装位置的机械共振能量散落在一片高频带上。直接对原始信号做 FFT故障特征频率的谱线往往被共振带和转频成分淹没肉眼很难锁定。EMD经验模态分解的好处在于不需要预设基函数它按信号自身的时间尺度从高频到低频逐层剥离出本征模态函数IMF把「故障冲击调制的共振成分」和「转频、噪声」分开再对含故障信息的 IMF 做包络谱和能量谱故障频率就能以明显峰值的形式暴露出来。下面就把这条链路上的 EMD 分解程序、包络谱计算、能量谱统计和特征频率定位一次讲清楚。适合正在做设备故障诊断、状态监测和 PHM 的工程师也适合用 MATLAB 做故障诊断方向实验和课题的同学。2. EMD 分解的原理与 MATLAB 里能直接跑的最小程序2.1 EMD 分解在做什么把非平稳信号剥成 IMFEMD 的核心假设是任何复杂信号都能拆成有限个本征模态函数加一个残余项。每个 IMF 必须满足两个条件第一整个数据段内极值点个数与过零点个数相等或最多相差一个第二任意时刻由上下包络线确定的均值接近零。分解过程就是反复的「筛分sifting」先找信号的局部极大值和极小值用三次样条拟合出上包络、下包络取均值后从原信号中减掉再对剩余信号重复上述操作直到满足 IMF 条件。这样剥出的第一个 IMF 是原始信号里最高频的成分然后继续对残余信号做同样处理下一层频率更低直到残余项变成单调趋势。这个自适应分解是 EMD 区别于小波和滤波器组的根本点。小波要选基函数和分解层数选得不合适特征就被改变了EMD 不需要完全由数据驱动。对电机轴承这种「转频 共振冲击 噪声」的混合信号冲击成分往往落在第一个或第二个 IMF 里而转频、谐波会留在低频 IMF 或残余项中这正好为后续包络谱定位特征频率提供了干净的输入。2.2 MATLAB 里跑通 EMD 的最小命令从 R2018a 开始Signal Processing Toolbox 就自带emd()函数不需要找第三方代码包。更早版本中常用的第三方 EMD 工具箱写法稍有区别但核心筛分逻辑一样。如果你的环境里help emd能正常出文档就直接按下面的方式来。先用一段合成的轴承外圈故障信号做演示方便对照结果是否正常fs 12000; % 采样率 12 kHz t 0:1/fs:1.5; % 1.5 秒振动数据 fr 30; % 轴频 30 Hz约 1800 r/min x_rot sin(2*pi*fr*t) 0.2*sin(2*pi*2*fr*t); % 转频及其二倍频 f_bpfo 107.5; % 6205 轴承外圈故障特征频率30 Hz 下计算值 tau 1/f_bpfo; % 冲击间隔 impact 0.6*exp(-700*mod(t, tau)) .* sin(2*pi*1500*mod(t, tau)); x x_rot impact 0.15*randn(size(t)); [imf, residual, info] emd(x);合成信号的思路是mod(t, tau)生成一个每隔tau秒归零的锯齿波exp(-700*mod(...))让每个周期内的冲击按指数衰减再乘上 1500 Hz 的正弦去模拟冲击激励起结构共振后的衰减振荡。这是轴承故障诊断实验里最常用的冲击模型真实传感器数据里能看到完全一样的波形形态只是噪声更大、混叠更多。emd(x)返回三个量。imf是分解出的本征模态函数矩阵每一行是一个 IMF顺序按频率从高到低residual是最后剩下的残余项一般对应信号的趋势或直流成分info里记录了筛分迭代次数、极值点数、实际分解出的 IMF 数量等过程信息。把这三个返回值搞清楚后面所有 IMF 筛选和特征提取都是在imf上进行的。2.3 分解结果的查看与常用展示写法EMD 跑完后第一件事是看图。常见做法是把前几个 IMF 和残余项竖着叠在一起用同一时间轴对比一行一个分量n_plot min(6, size(imf, 1)); figure(Color, w); for k 1:n_plot subplot(n_plot 1, 1, k); plot(t, imf(k, :), LineWidth, 0.6); ylabel([IMF, num2str(k)]); xlim([0.2, 0.6]); % 只看局部一段否则冲击细节看不清 end subplot(n_plot 1, 1, n_plot 1); plot(t, residual, LineWidth, 0.6); ylabel(residual); xlabel(时间 (s));这里的xlim([0.2 0.6])是很多人忽略的一步。1.5 秒信号里包含上百个冲击全部画出来时每个冲击只占几个采样点根本看不出衰减波形放大到 0.4 秒窗口后IMF1 里的冲击序列、IMF2 和 IMF3 里的转频正弦都能直观区分。如果发现某个 IMF 两端出现大幅飞翼整条线抖得离谱那是端点效应处理方法在第 4 章细说。对于长数据emd()的耗时会随信号长度明显上升。如果只是为了诊断轴承故障不需要把几十万点的数据整段丢进去截取 25 秒、包住至少几十个冲击周期就足够后面算包络谱时频率分辨率也够用。提示emd()依赖 Signal Processing Toolbox。运行前用license(test, Signal_Toolbox)确认工具箱可用返回 1 再执行分解。3. 包络谱与能量谱EMD 分解输出怎么变成故障判据3.1 包络谱为什么能定位故障频率轴承故障冲击以特征频率重复出现但它不是直接体现在振动幅值上而是对结构共振频率进行幅值调制。直接对原始信号做 FFT看到的是共振频率附近的一大片边带特征频率本身并不单独成线。包络分析的思路恰恰是解调先把信号的幅值包络提取出来去掉高频载波只剩下「冲击重复频率」这个低频调制信息再对这个包络做 FFT特征频率处就会形成明显的单根谱峰。在 MATLAB 里提取包络最常用的工具是希尔伯特变换。函数hilbert(sig)返回解析信号实部是原信号虚部是希尔伯特变换结果取绝对值就得到瞬时幅值包络。对轴承故障诊断来说包络谱的物理含义就是「每秒多少个冲击」的频谱而冲击重复频率正好对应外圈、内圈、滚动体或保持架的特征频率因此包络谱峰值定位是整个诊断流程中最关键的一步。3.2 包络谱的 MATLAB 计算步骤把包络谱计算写成独立函数后续对任意 IMF 都能直接调用function [f_axis, env_spec] envelopeSpectrum(sig, fs) n length(sig); analytic hilbert(sig); % 解析信号 env abs(analytic); % 幅值包络 env env - mean(env); % 去直流避免 0 Hz 处能量堆积 spec fft(env); spec abs(spec / n); % 归一化幅值 one_sided spec(1:fix(n/2) 1); one_sided(2:end-1) 2 * one_sided(2:end-1); % 单边谱幅值还原 f_axis fs * (0:(n/2)) / n; env_spec one_sided; end逐行解释一下。hilbert(sig)得到的是复信号abs()取模后就是包络这一步同时完成了正交解调。减去均值是因为包络恒为正含有很大的直流分量不去掉的话 0 Hz 处会出现一个压倒性的谱峰把故障频率附近的细节全部遮住。后半段是做单边频谱的常规处理FFT 结果关于 fs/2 对称只取前半段中间频率以外的谱线幅值乘 2 还原真实幅值。对分解出的 IMF1 调用这个函数再在 0500 Hz 范围内找峰值[f_axis, env_spec] envelopeSpectrum(imf(1, :), fs); [pks, locs] findpeaks(env_spec, f_axis, ... MinPeakHeight, 0.1*max(env_spec), ... MinPeakDistance, 2); figure(Color, w); plot(f_axis, env_spec, LineWidth, 0.8); xlim([0 400]); xlabel(频率 (Hz)); ylabel(包络谱幅值); hold on; plot(locs, pks, rv, MarkerSize, 6);findpeaks的两个参数值得说明MinPeakHeight设成最大幅值的 10%用来滤掉噪声底上的小毛刺MinPeakDistance设成 2 Hz保证同一频率附近的边带峰不会被重复计入。运行后会在 107.5 Hz 附近看到一个明显主峰这就是外圈故障特征频率。如果对原始信号直接做这个操作同样的峰也存在但幅值和信噪比会明显差一截这正是先做 EMD 分离的价值。3.3 能量谱统计每个 IMF 的能量占比包络谱解决的是「故障频率落在哪里」的问题能量谱解决的是「故障能量分散在哪些尺度」的问题。每个 IMF 是相互正交近似独立的窄带分量对其求平方和再归一化就能得到信号能量在不同频率尺度的分布energy sum(imf.^2, 2); % 每一行 IMF 的能量 energy_ratio energy / sum(energy) * 100; figure(Color, w); bar(1:size(imf, 1), energy_ratio, FaceColor, [0.2 0.4 0.8]); xlabel(IMF 序号); ylabel(能量占比 (%)); grid on;健康状态下振动能量主要分布在转频及其低次谐波对应的低频 IMF 里高频 IMF 能量占比很小。轴承出现剥落、点蚀后冲击能量被激励到共振频段对应的高频 IMF 能量占比会明显抬升。实际做监测时通常是采集同型号电机健康与故障两组数据对比同一 IMF 序号的能量占比变化这个比值比单看幅值稳定得多也适合做趋势预警。更细致一点可以对每个 IMF 的频谱平方做积分得到频带能量但对轴承诊断而言时域平方和加包络谱已经能覆盖绝大多数场景。输出量计算位置在诊断里的用途IMFemd()第一返回值分离冲击调制成分供后续分析包络谱对选定 IMF 做hilbertfft在特征频率处找峰值判定故障类型能量谱对每个 IMF 求平方和再归一化判断能量集中在哪个尺度对比健康/故障样本4. 特征频率定位、IMF 筛选与 EMD 分解的必调参数4.1 轴承故障特征频率怎么算包络谱上的峰值要能和理论特征频率对上诊断结论才有意义。特征频率由轴承几何参数和转频决定计算前需要知道滚动体个数 Z、滚动体直径 d、节圆直径 D、接触角 α 和当前转频 fr。电机轴承的转频直接用fr rpm / 60换算四项特征频率按下面公式计算故障部位特征频率计算公式外圈 BPFOZ/2 * fr * (1 - d/D * cos(alpha))内圈 BPFIZ/2 * fr * (1 d/D * cos(alpha))滚动体 BSFD/(2*d) * fr * (1 - (d/D)^2 * cos(alpha)^2)保持架 FTFfr/2 * (1 - d/D * cos(alpha))以电机测试里常用的 6205 深沟球轴承为例Z9d7.94 mmD39.04 mm接触角 α0°转频 fr30 Hz。代入外圈公式得到 BPFO 9/2 × 30 × (1 − 7.94/39.04) ≈ 107.5 Hz前面的合成信号就是按这个值生成的。实际设备中转速会有波动包络谱峰可能偏移 12 Hz所以找峰值做匹配时要留容差。把轴承型号参数和转频封装成函数批量诊断时每个样本都能自动算出理论特征频率不用每次手算。4.2 用峭度和相关系数筛选有效 IMFEMD 分解出的 IMF 不是每一层都有诊断价值。第一层往往噪声占主导最后几层基本是低频趋势直接拿 IMF1 做包络谱可能漏掉故障信息。常用的筛选指标是峭度加相关系数两者做乘积排序kurt_val kurtosis(imf, 0, 2); % 每个 IMF 的峭度 corr_val zeros(size(imf, 1), 1); for k 1:size(imf, 1) corr_val(k) corr(imf(k, :), x); % 与原始信号的相关系数 end score kurt_val .* corr_val; [~, best] max(score); [fb, sb] envelopeSpectrum(imf(best, :), fs);峭度度量的是信号冲击性正常振动近似高斯分布峭度接近 3轴承出现早期缺陷时冲击成分会使峭度升到 410且故障越早期越明显。但峭度高的 IMF 也可能是纯噪声尖峰所以乘上与原始信号的相关系数保留那些既像原信号又带强冲击性的分量。实际项目里用这个乘积排序命中率比单看峭度高不少。筛选之后还有个常用动作把选中的 IMF 进行信号重构即把筛选出的几个 IMF 相加替代原始信号做后续分析。这样既保留了故障冲击成分又丢掉了转频干扰和部分噪声包络谱的底噪会更平特征频率峰更突出。4.3 emd() 的参数怎么调与三个常见坑emd()提供几个关键参数默认值能跑通大多数数据但遇到长信号、碎片化信号时值得显式设置参数名默认值作用与调节建议MaxNumIMF不限限制最大分解层数防止算出过多碎片 IMF诊断场景设 58SiftRelativeTolerance0.1/0.2筛分停止容差调小分解更精细但更慢噪声大时反而该调大MaxNumSiftIterations100单次筛分迭代上限信号碎片多时调大到 200Display0置 1 可观察分解进度和每层迭代次数排错时打开调用方式是emd(x, MaxNumIMF, 6, SiftRelativeTolerance, 0.2, Display, 1)。诊断场景中限制MaxNumIMF是最实用的一招默认情况下 EMD 可能分解出十几个 IMF后几层只是把残余趋势拆得更碎既拖慢速度又干扰筛选排序限制层数后计算量明显下降诊断精度基本不受影响。三个常见坑需要在实际数据上特别注意。第一是端点效应信号两端缺少极值点约束包络拟合发散导致 IMF 两端出现大幅飞翼。缓解办法是截取分析段时前后各多留 0.20.5 秒数据分解完丢掉两端限制MaxNumIMF也能减小低层 IMF 的端点畸变。第二是模式混叠故障冲击如果间歇出现同一个特征频率可能被拆到两个 IMF 里包络谱上表现为特征频率两侧散开的多个小峰主峰幅值明显低于预期。如果确认波形里有冲击但包络谱峰值弱可以按集合经验模态分解EEMD的思路处理对叠加不同白噪声的多个副本分别做 EMD 再平均代价是计算量成倍上升一般诊断场景先用普通 EMD。第三是采样率过低特征频率本身多在 1 kHz 以内但冲击激励的共振常到几 kHz 甚至更高采样率低于 10 kHz 时包络波形会失真建议电机轴承诊断的数据至少用 1020 kHz 采样。5. 把 EMD 轴承故障诊断流程封装成可复用的 MATLAB 函数5.1 封装函数一次跑完分解、筛选、包络谱与匹配前面的步骤单独跑没问题但换成批量数据就会手忙脚乱。常见做法是把整条链路封装成一个函数输入振动信号、采样率和待核查特征频率输出诊断中间量function out bearingEmdDiagnosis(x, fs, f_targets, varargin) % x 振动信号列向量 % fs 采样率 % f_targets 待核查的故障特征频率如 [107.5 162.2] p inputParser; addParameter(p, MaxNumIMF, 6, (v) isnumeric(v) v 0); addParameter(p, SiftRelativeTolerance, 0.2, (v) isnumeric(v) v 0); parse(p, varargin{:}); [imf, ~, ~] emd(x, MaxNumIMF, p.Results.MaxNumIMF, ... SiftRelativeTolerance, p.Results.SiftRelativeTolerance, ... Display, 0); kurt_val kurtosis(imf, 0, 2); corr_val zeros(size(imf, 1), 1); for k 1:size(imf, 1) corr_val(k) corr(imf(k, :), x); end score kurt_val .* corr_val; [~, best] max(score); [f_axis, spec] envelopeSpectrum(imf(best, :), fs); [pks, locs] findpeaks(spec, f_axis, ... MinPeakHeight, 0.3*max(spec), MinPeakDistance, 2); match zeros(size(f_targets)); for k 1:numel(f_targets) delta abs(locs - f_targets(k)); if any(delta 2) match(k) 1; end end out.bestIMF best; out.kurtosis kurt_val; out.featureTable table(f_targets, match, ... VariableNames, {FaultHz, Matched}); out.envelopeFreq f_axis; out.envelopeSpec spec; out.peaks table(locs, pks); end函数里的inputParser用来处理可选参数调用方可以写bearingEmdDiagnosis(x, fs, [107.5], MaxNumIMF, 8)覆盖默认值。findpeaks的MinPeakHeight用最大幅值的比例避免不同样本幅值差异导致阈值失效。匹配逻辑是计算包络谱峰位置与理论特征频率的差值小于 2 Hz 即判为命中对应前面说的转速波动容差。5.2 用已知故障频率验证整套流程封装之后必须用已知答案验证一遍。拿第 2 章的合成信号理论外圈特征频率就是 107.5 Hz运行out bearingEmdDiagnosis(x, fs, [107.5]); disp(out.featureTable);正常情况下Matched列返回 1。再拿一段不含冲击的纯转频加噪声信号跑同一个函数Matched应为 0。两个结果都符合预期说明分解、筛选、包络谱、匹配整条链路是通的之后换成实测数据才有可信度。实测数据验证时建议先取一个已知健康状态的样本做基线记录包络谱噪声中位数把故障判定的阈值从「固定比例」改成「超过基线中位数的若干倍」这样对不同机组、不同测点位置的适应性强得多。对批量数据做筛查时把每次调用的out.featureTable逐行追加成总表几万组样本的筛选结果就能直接和故障频率统计或机器学习分类器的输入特征对接。本文还有配套的精品资源点击获取 EMD分解与包络谱在轴承故障诊断中的应用
简介这套电机轴承故障诊断程序基于MATLAB实现采用经验模态分解EMD方法处理非线性、非平稳振动信号面向机械故障诊断研究人员、电机维护工程师及信号处理相关专业学生。包内包含8个文件以7个.m脚本和1个.mat数据文件为主脚本涵盖EMD分解、希尔伯特包络谱、能量谱计算以及内圈/外圈故障分析等完整流程数据文件可直接用于算法验证。整套程序压缩后仅1.14MB结构简洁便于理解与应用。已有247人学习下载实用性强。通过学习这份代码使用者能掌握从振动信号采集到EMD分解、IMF分量提取、包络谱与能量谱构建、故障特征频率识别的全链路方法并根据频谱峰值判断轴承健康状态为滚动体缺陷、内外圈故障等典型失效模式提供诊断思路。同时该案例也是非线性信号处理在工程中落地的良好范本适合作为课题研究或课程设计的参考实现。1. 电机轴承故障诊断为什么要先用 EMD 分解电机轴承出现早期缺陷时传感器采到的振动信号并不是干净的周期冲击。转频及其谐波、电机电气干扰、现场噪声全混在一起故障冲击又会激起轴承座和传感器安装位置的机械共振能量散落在一片高频带上。直接对原始信号做 FFT故障特征频率的谱线往往被共振带和转频成分淹没肉眼很难锁定。EMD经验模态分解的好处在于不需要预设基函数它按信号自身的时间尺度从高频到低频逐层剥离出本征模态函数IMF把「故障冲击调制的共振成分」和「转频、噪声」分开再对含故障信息的 IMF 做包络谱和能量谱故障频率就能以明显峰值的形式暴露出来。下面就把这条链路上的 EMD 分解程序、包络谱计算、能量谱统计和特征频率定位一次讲清楚。适合正在做设备故障诊断、状态监测和 PHM 的工程师也适合用 MATLAB 做故障诊断方向实验和课题的同学。2. EMD 分解的原理与 MATLAB 里能直接跑的最小程序2.1 EMD 分解在做什么把非平稳信号剥成 IMFEMD 的核心假设是任何复杂信号都能拆成有限个本征模态函数加一个残余项。每个 IMF 必须满足两个条件第一整个数据段内极值点个数与过零点个数相等或最多相差一个第二任意时刻由上下包络线确定的均值接近零。分解过程就是反复的「筛分sifting」先找信号的局部极大值和极小值用三次样条拟合出上包络、下包络取均值后从原信号中减掉再对剩余信号重复上述操作直到满足 IMF 条件。这样剥出的第一个 IMF 是原始信号里最高频的成分然后继续对残余信号做同样处理下一层频率更低直到残余项变成单调趋势。这个自适应分解是 EMD 区别于小波和滤波器组的根本点。小波要选基函数和分解层数选得不合适特征就被改变了EMD 不需要完全由数据驱动。对电机轴承这种「转频 共振冲击 噪声」的混合信号冲击成分往往落在第一个或第二个 IMF 里而转频、谐波会留在低频 IMF 或残余项中这正好为后续包络谱定位特征频率提供了干净的输入。2.2 MATLAB 里跑通 EMD 的最小命令从 R2018a 开始Signal Processing Toolbox 就自带emd()函数不需要找第三方代码包。更早版本中常用的第三方 EMD 工具箱写法稍有区别但核心筛分逻辑一样。如果你的环境里help emd能正常出文档就直接按下面的方式来。先用一段合成的轴承外圈故障信号做演示方便对照结果是否正常fs 12000; % 采样率 12 kHz t 0:1/fs:1.5; % 1.5 秒振动数据 fr 30; % 轴频 30 Hz约 1800 r/min x_rot sin(2*pi*fr*t) 0.2*sin(2*pi*2*fr*t); % 转频及其二倍频 f_bpfo 107.5; % 6205 轴承外圈故障特征频率30 Hz 下计算值 tau 1/f_bpfo; % 冲击间隔 impact 0.6*exp(-700*mod(t, tau)) .* sin(2*pi*1500*mod(t, tau)); x x_rot impact 0.15*randn(size(t)); [imf, residual, info] emd(x);合成信号的思路是mod(t, tau)生成一个每隔tau秒归零的锯齿波exp(-700*mod(...))让每个周期内的冲击按指数衰减再乘上 1500 Hz 的正弦去模拟冲击激励起结构共振后的衰减振荡。这是轴承故障诊断实验里最常用的冲击模型真实传感器数据里能看到完全一样的波形形态只是噪声更大、混叠更多。emd(x)返回三个量。imf是分解出的本征模态函数矩阵每一行是一个 IMF顺序按频率从高到低residual是最后剩下的残余项一般对应信号的趋势或直流成分info里记录了筛分迭代次数、极值点数、实际分解出的 IMF 数量等过程信息。把这三个返回值搞清楚后面所有 IMF 筛选和特征提取都是在imf上进行的。2.3 分解结果的查看与常用展示写法EMD 跑完后第一件事是看图。常见做法是把前几个 IMF 和残余项竖着叠在一起用同一时间轴对比一行一个分量n_plot min(6, size(imf, 1)); figure(Color, w); for k 1:n_plot subplot(n_plot 1, 1, k); plot(t, imf(k, :), LineWidth, 0.6); ylabel([IMF, num2str(k)]); xlim([0.2, 0.6]); % 只看局部一段否则冲击细节看不清 end subplot(n_plot 1, 1, n_plot 1); plot(t, residual, LineWidth, 0.6); ylabel(residual); xlabel(时间 (s));这里的xlim([0.2 0.6])是很多人忽略的一步。1.5 秒信号里包含上百个冲击全部画出来时每个冲击只占几个采样点根本看不出衰减波形放大到 0.4 秒窗口后IMF1 里的冲击序列、IMF2 和 IMF3 里的转频正弦都能直观区分。如果发现某个 IMF 两端出现大幅飞翼整条线抖得离谱那是端点效应处理方法在第 4 章细说。对于长数据emd()的耗时会随信号长度明显上升。如果只是为了诊断轴承故障不需要把几十万点的数据整段丢进去截取 25 秒、包住至少几十个冲击周期就足够后面算包络谱时频率分辨率也够用。提示emd()依赖 Signal Processing Toolbox。运行前用license(test, Signal_Toolbox)确认工具箱可用返回 1 再执行分解。3. 包络谱与能量谱EMD 分解输出怎么变成故障判据3.1 包络谱为什么能定位故障频率轴承故障冲击以特征频率重复出现但它不是直接体现在振动幅值上而是对结构共振频率进行幅值调制。直接对原始信号做 FFT看到的是共振频率附近的一大片边带特征频率本身并不单独成线。包络分析的思路恰恰是解调先把信号的幅值包络提取出来去掉高频载波只剩下「冲击重复频率」这个低频调制信息再对这个包络做 FFT特征频率处就会形成明显的单根谱峰。在 MATLAB 里提取包络最常用的工具是希尔伯特变换。函数hilbert(sig)返回解析信号实部是原信号虚部是希尔伯特变换结果取绝对值就得到瞬时幅值包络。对轴承故障诊断来说包络谱的物理含义就是「每秒多少个冲击」的频谱而冲击重复频率正好对应外圈、内圈、滚动体或保持架的特征频率因此包络谱峰值定位是整个诊断流程中最关键的一步。3.2 包络谱的 MATLAB 计算步骤把包络谱计算写成独立函数后续对任意 IMF 都能直接调用function [f_axis, env_spec] envelopeSpectrum(sig, fs) n length(sig); analytic hilbert(sig); % 解析信号 env abs(analytic); % 幅值包络 env env - mean(env); % 去直流避免 0 Hz 处能量堆积 spec fft(env); spec abs(spec / n); % 归一化幅值 one_sided spec(1:fix(n/2) 1); one_sided(2:end-1) 2 * one_sided(2:end-1); % 单边谱幅值还原 f_axis fs * (0:(n/2)) / n; env_spec one_sided; end逐行解释一下。hilbert(sig)得到的是复信号abs()取模后就是包络这一步同时完成了正交解调。减去均值是因为包络恒为正含有很大的直流分量不去掉的话 0 Hz 处会出现一个压倒性的谱峰把故障频率附近的细节全部遮住。后半段是做单边频谱的常规处理FFT 结果关于 fs/2 对称只取前半段中间频率以外的谱线幅值乘 2 还原真实幅值。对分解出的 IMF1 调用这个函数再在 0500 Hz 范围内找峰值[f_axis, env_spec] envelopeSpectrum(imf(1, :), fs); [pks, locs] findpeaks(env_spec, f_axis, ... MinPeakHeight, 0.1*max(env_spec), ... MinPeakDistance, 2); figure(Color, w); plot(f_axis, env_spec, LineWidth, 0.8); xlim([0 400]); xlabel(频率 (Hz)); ylabel(包络谱幅值); hold on; plot(locs, pks, rv, MarkerSize, 6);findpeaks的两个参数值得说明MinPeakHeight设成最大幅值的 10%用来滤掉噪声底上的小毛刺MinPeakDistance设成 2 Hz保证同一频率附近的边带峰不会被重复计入。运行后会在 107.5 Hz 附近看到一个明显主峰这就是外圈故障特征频率。如果对原始信号直接做这个操作同样的峰也存在但幅值和信噪比会明显差一截这正是先做 EMD 分离的价值。3.3 能量谱统计每个 IMF 的能量占比包络谱解决的是「故障频率落在哪里」的问题能量谱解决的是「故障能量分散在哪些尺度」的问题。每个 IMF 是相互正交近似独立的窄带分量对其求平方和再归一化就能得到信号能量在不同频率尺度的分布energy sum(imf.^2, 2); % 每一行 IMF 的能量 energy_ratio energy / sum(energy) * 100; figure(Color, w); bar(1:size(imf, 1), energy_ratio, FaceColor, [0.2 0.4 0.8]); xlabel(IMF 序号); ylabel(能量占比 (%)); grid on;健康状态下振动能量主要分布在转频及其低次谐波对应的低频 IMF 里高频 IMF 能量占比很小。轴承出现剥落、点蚀后冲击能量被激励到共振频段对应的高频 IMF 能量占比会明显抬升。实际做监测时通常是采集同型号电机健康与故障两组数据对比同一 IMF 序号的能量占比变化这个比值比单看幅值稳定得多也适合做趋势预警。更细致一点可以对每个 IMF 的频谱平方做积分得到频带能量但对轴承诊断而言时域平方和加包络谱已经能覆盖绝大多数场景。输出量计算位置在诊断里的用途IMFemd()第一返回值分离冲击调制成分供后续分析包络谱对选定 IMF 做hilbertfft在特征频率处找峰值判定故障类型能量谱对每个 IMF 求平方和再归一化判断能量集中在哪个尺度对比健康/故障样本4. 特征频率定位、IMF 筛选与 EMD 分解的必调参数4.1 轴承故障特征频率怎么算包络谱上的峰值要能和理论特征频率对上诊断结论才有意义。特征频率由轴承几何参数和转频决定计算前需要知道滚动体个数 Z、滚动体直径 d、节圆直径 D、接触角 α 和当前转频 fr。电机轴承的转频直接用fr rpm / 60换算四项特征频率按下面公式计算故障部位特征频率计算公式外圈 BPFOZ/2 * fr * (1 - d/D * cos(alpha))内圈 BPFIZ/2 * fr * (1 d/D * cos(alpha))滚动体 BSFD/(2*d) * fr * (1 - (d/D)^2 * cos(alpha)^2)保持架 FTFfr/2 * (1 - d/D * cos(alpha))以电机测试里常用的 6205 深沟球轴承为例Z9d7.94 mmD39.04 mm接触角 α0°转频 fr30 Hz。代入外圈公式得到 BPFO 9/2 × 30 × (1 − 7.94/39.04) ≈ 107.5 Hz前面的合成信号就是按这个值生成的。实际设备中转速会有波动包络谱峰可能偏移 12 Hz所以找峰值做匹配时要留容差。把轴承型号参数和转频封装成函数批量诊断时每个样本都能自动算出理论特征频率不用每次手算。4.2 用峭度和相关系数筛选有效 IMFEMD 分解出的 IMF 不是每一层都有诊断价值。第一层往往噪声占主导最后几层基本是低频趋势直接拿 IMF1 做包络谱可能漏掉故障信息。常用的筛选指标是峭度加相关系数两者做乘积排序kurt_val kurtosis(imf, 0, 2); % 每个 IMF 的峭度 corr_val zeros(size(imf, 1), 1); for k 1:size(imf, 1) corr_val(k) corr(imf(k, :), x); % 与原始信号的相关系数 end score kurt_val .* corr_val; [~, best] max(score); [fb, sb] envelopeSpectrum(imf(best, :), fs);峭度度量的是信号冲击性正常振动近似高斯分布峭度接近 3轴承出现早期缺陷时冲击成分会使峭度升到 410且故障越早期越明显。但峭度高的 IMF 也可能是纯噪声尖峰所以乘上与原始信号的相关系数保留那些既像原信号又带强冲击性的分量。实际项目里用这个乘积排序命中率比单看峭度高不少。筛选之后还有个常用动作把选中的 IMF 进行信号重构即把筛选出的几个 IMF 相加替代原始信号做后续分析。这样既保留了故障冲击成分又丢掉了转频干扰和部分噪声包络谱的底噪会更平特征频率峰更突出。4.3 emd() 的参数怎么调与三个常见坑emd()提供几个关键参数默认值能跑通大多数数据但遇到长信号、碎片化信号时值得显式设置参数名默认值作用与调节建议MaxNumIMF不限限制最大分解层数防止算出过多碎片 IMF诊断场景设 58SiftRelativeTolerance0.1/0.2筛分停止容差调小分解更精细但更慢噪声大时反而该调大MaxNumSiftIterations100单次筛分迭代上限信号碎片多时调大到 200Display0置 1 可观察分解进度和每层迭代次数排错时打开调用方式是emd(x, MaxNumIMF, 6, SiftRelativeTolerance, 0.2, Display, 1)。诊断场景中限制MaxNumIMF是最实用的一招默认情况下 EMD 可能分解出十几个 IMF后几层只是把残余趋势拆得更碎既拖慢速度又干扰筛选排序限制层数后计算量明显下降诊断精度基本不受影响。三个常见坑需要在实际数据上特别注意。第一是端点效应信号两端缺少极值点约束包络拟合发散导致 IMF 两端出现大幅飞翼。缓解办法是截取分析段时前后各多留 0.20.5 秒数据分解完丢掉两端限制MaxNumIMF也能减小低层 IMF 的端点畸变。第二是模式混叠故障冲击如果间歇出现同一个特征频率可能被拆到两个 IMF 里包络谱上表现为特征频率两侧散开的多个小峰主峰幅值明显低于预期。如果确认波形里有冲击但包络谱峰值弱可以按集合经验模态分解EEMD的思路处理对叠加不同白噪声的多个副本分别做 EMD 再平均代价是计算量成倍上升一般诊断场景先用普通 EMD。第三是采样率过低特征频率本身多在 1 kHz 以内但冲击激励的共振常到几 kHz 甚至更高采样率低于 10 kHz 时包络波形会失真建议电机轴承诊断的数据至少用 1020 kHz 采样。5. 把 EMD 轴承故障诊断流程封装成可复用的 MATLAB 函数5.1 封装函数一次跑完分解、筛选、包络谱与匹配前面的步骤单独跑没问题但换成批量数据就会手忙脚乱。常见做法是把整条链路封装成一个函数输入振动信号、采样率和待核查特征频率输出诊断中间量function out bearingEmdDiagnosis(x, fs, f_targets, varargin) % x 振动信号列向量 % fs 采样率 % f_targets 待核查的故障特征频率如 [107.5 162.2] p inputParser; addParameter(p, MaxNumIMF, 6, (v) isnumeric(v) v 0); addParameter(p, SiftRelativeTolerance, 0.2, (v) isnumeric(v) v 0); parse(p, varargin{:}); [imf, ~, ~] emd(x, MaxNumIMF, p.Results.MaxNumIMF, ... SiftRelativeTolerance, p.Results.SiftRelativeTolerance, ... Display, 0); kurt_val kurtosis(imf, 0, 2); corr_val zeros(size(imf, 1), 1); for k 1:size(imf, 1) corr_val(k) corr(imf(k, :), x); end score kurt_val .* corr_val; [~, best] max(score); [f_axis, spec] envelopeSpectrum(imf(best, :), fs); [pks, locs] findpeaks(spec, f_axis, ... MinPeakHeight, 0.3*max(spec), MinPeakDistance, 2); match zeros(size(f_targets)); for k 1:numel(f_targets) delta abs(locs - f_targets(k)); if any(delta 2) match(k) 1; end end out.bestIMF best; out.kurtosis kurt_val; out.featureTable table(f_targets, match, ... VariableNames, {FaultHz, Matched}); out.envelopeFreq f_axis; out.envelopeSpec spec; out.peaks table(locs, pks); end函数里的inputParser用来处理可选参数调用方可以写bearingEmdDiagnosis(x, fs, [107.5], MaxNumIMF, 8)覆盖默认值。findpeaks的MinPeakHeight用最大幅值的比例避免不同样本幅值差异导致阈值失效。匹配逻辑是计算包络谱峰位置与理论特征频率的差值小于 2 Hz 即判为命中对应前面说的转速波动容差。5.2 用已知故障频率验证整套流程封装之后必须用已知答案验证一遍。拿第 2 章的合成信号理论外圈特征频率就是 107.5 Hz运行out bearingEmdDiagnosis(x, fs, [107.5]); disp(out.featureTable);正常情况下Matched列返回 1。再拿一段不含冲击的纯转频加噪声信号跑同一个函数Matched应为 0。两个结果都符合预期说明分解、筛选、包络谱、匹配整条链路是通的之后换成实测数据才有可信度。实测数据验证时建议先取一个已知健康状态的样本做基线记录包络谱噪声中位数把故障判定的阈值从「固定比例」改成「超过基线中位数的若干倍」这样对不同机组、不同测点位置的适应性强得多。对批量数据做筛查时把每次调用的out.featureTable逐行追加成总表几万组样本的筛选结果就能直接和故障频率统计或机器学习分类器的输入特征对接。本文还有配套的精品资源点击获取