前阵子帮一个产线朋友处理减速机振动数据他第一时间把FFT频谱发过来问“为什么频谱上找不到外圈故障的边带”其实这是很多刚接触滚动轴承故障诊断的人都会卡住的地方——不是FFT算错了而是选错了分析频段。滚动轴承早期故障产生的冲击能量非常小它的特征被调制在高频共振带上直接对原始信号做FFT那点冲击成分早被淹没在宽频噪声里了。要解开这个结我们就得用谱峭度先把整个频带“扫描”一遍锁定故障冲击最集中的共振频带再对该频带做包络谱分析。这篇就是基于MATLAB快速谱峭度包络谱整套流程的完整记录里面包含我实际调试时踩过的坑和总结出的参数选择经验适合正在做故障诊断课程设计、毕业论文以及刚入门设备状态监测的工程师参考。1. 为什么要用谱峭度和包络谱轴承故障信号的“隐身术”与“现形记”1.1 直接做FFT为什么看不到故障频率很多教材里有个典型说法轴承外圈故障振动信号中会出现BPFO及其谐波成分。于是大家拿到数据后第一时间做FFT期望在频谱里看到一根根清晰的谱线结果往往是整个低频段乱糟糟一片连转频都找不仔细更别提BPFO了。这不是FFT代码写错了而是实际轴承振动信号的结构远比理想模型复杂。轴承故障的物理过程是滚动体滚过缺陷时产生一个极短的冲击脉冲这个脉冲的频带非常宽会激励起轴承座、端盖、传感器安装部位的固有共振。我们实际采集到的信号可以近似看作一串窄带高频衰减振荡振荡频率即系统共振频率重复频率即故障特征频率。用调幅信号的语言来说就是“高频载波被低频脉冲调制”。直接对这样的信号做FFT我们能看到的频谱主体集中在共振频率附近而故障特征频率由于调制效应会以边频带形式出现但这些边频带的幅值相比共振峰小得多很容易被噪声掩盖在频谱上看就是“什么都没有”。1.2 包络分析的本质把高频共振信号搬回低频包络谱分析Envelope Spectrum Analysis的思路很直接既然故障信息藏在调制包络里我们就把包络提取出来再对包络做FFT。这样调制边带被解调回故障特征频率及其谐波在低频段形成清晰的谱线。实现上最常用的工具是Hilbert变换得到信号的解析信号后取模得到的包络信号再FFT就是包络谱。对于包络分析关键前置步骤是带通滤波——只有把分析频带对准包含故障冲击能量的共振频带解调出来的信号才足够干净。1.3 谱峭度扮演的角色自动寻找最佳共振频带既然包络分析需要一个“合适的频带”那问题就变成了这个频带在哪传统做法是用经验反复尝试比如在2kHz到10kHz之间几个频带挨个试看哪个包络谱效果最好。试过的人都知道这非常依赖经验而且不同转速、不同工况下共振频带会改变换一台设备就得重新试。谱峭度Spectral Kurtosis就是来解决这个问题的。峭度是四阶统计量用来衡量信号偏离高斯分布的程度。正常轴承运行时的振动信号接近高斯分布峭度接近0一旦出现冲击性故障信号中会出现大量“尖峰”峭度显著升高。谱峭度则是把峭度概念推广到每个频率点计算每个频带上的峭度值然后画成一张“频率-峭度”图。峭度最大的频带就是包含最多冲击能量、最适合做包络分析的频带。快速谱峭度Fast Kurtogram由Antoni在2007年提出通过1/3-二叉树结构高效划分频带计算每个子带的谱峭度值形成一张二维色图。这张图可以直观地告诉我们故障冲击能量集中在哪个频率范围、哪个频带宽度最合适。有了它包络谱分析就从一个“试错游戏”变成了一个有明确依据的工程流程。2. 快速谱峭度的核心逻辑从峭度的数学定义到1/3二叉树频带划分2.1 峭度的直觉理解与数学形式先聊透峭度。对一段时域信号x其峭度定义为K E[(x - μ)⁴] / σ⁴其中μ是均值σ是标准差。正态分布的峭度刚好是3所以工程上常用超额峭度K-3把高斯信号归零。峭度大说明信号中存在大量远离均值的“尖峰”。你可以把它想象成一个“尖峰探测器”一个班的学生体重如果大家都差不多分布接近正态如果突然有极少数人胖得离谱分布就会出现明显的“重尾”峭度就会明显上升。轴承故障信号里滚动体周期性滚过缺陷产生的冲击脉冲就是这些“极少数胖得离谱的学生”。峭度越大说明冲击成分在信号中占有的能量比例越高。这也是为什么很多状态监测系统直接把时域峭度作为早期报警指标——简单、快、无需任何先验频率信息。2.2 谱峭度在每个频带上分别计算峭度时域峭度只能反映整个信号的整体冲击性无法告诉我们“冲击能量集中在哪个频带”。谱峭度的改进正是针对这一点先把信号做短时傅里叶变换STFT在每个频率点f上得到一段复值时间序列X(t, f)然后计算这个序列的四阶累积量K(f) E{|X(t,f)|⁴} / (E{|X(t,f)|²})² - 2这里的核心思想是如果某个频带内信号只是平稳高斯噪声那么X(t,f)的幅度分布近似瑞利分布K(f)约等于0如果该频带内叠加了周期性冲击引起的共振衰减振荡X(t,f)的幅度会出现明显的闪变K(f)就会显著增大。所以谱峭度本质上是在逐个频带上做“冲击性检验”把所有可能包含故障信息的频带都筛出来。2.3 快速算法的关键1/3二叉树结构与工程意义直接在每个频率点都做一次四阶矩估计计算量非常大不适合工程应用。Antoni的快速谱峭度算法巧妙地用带通滤波器组代替逐点计算采用1/3二叉树结构把频带逐层细分第0层是整个频带第1层分为两个子带第2层每个子带再分两个依此类推同时在每层之间插入1/3倍频程的频带划分保证频率分辨率与带宽的合理平衡。每个子带的中心频率和带宽确定后计算滤波输出信号的复包络峭度最终生成Kurtogram色图。这个结构有两个工程上的直接好处。第一是计算速度极快实测一段10秒、采样率25.6kHz的信号在普通笔记本上完成6层分解加绘图的耗时通常不超过两秒完全满足离线分析和在线快速诊断的需要。第二是频带信息直观色图的横轴是频率纵轴是分解层数对应不同带宽颜色代表峭度值。算法自动推荐峭度最大的子带但工程师也能通过色图看到次优频带、多频带为后续人工判断提供依据。2.4 快速谱峭度工具箱的基本用法在MATLAB中做快速谱峭度可以直接使用John Antoni公开的Kurtogram工具箱。我自己的用法是把工具箱文件全部加入当前工程目录然后用下面这段代码调用核心函数% x: 振动加速度信号 % fs: 采样频率 % level: 最大分解层数经验取4~6 [K, level_f, freq_f] fast_kurtogram(x, fs, level); % 找到峭度最大的子带 [Kmax, idx] max(K(:));其中level_f和freq_f分别记录每个子带的分解层数对应带宽和中心频率。拿到Kmax对应的索引后就可以换算最优带通滤波器的上下限频率。需要注意的是Kurtogram工具箱输出的是特征频率向量与峭度矩阵使用时务必检查矩阵维度与频带边界对应关系别在索引换算时弄错。3. MATLAB实现全流程从振动数据到最优共振频带的完整代码3.1 工程问题描述与参数准备我们以一套典型的滚动轴承故障诊断为场景电机转速约为1800 rpm转频fr30 Hz加速度传感器安装在轴承座上采样频率fs25600 Hz采样时长10秒。被测轴承为深沟球轴承6205节圆直径D39.04 mm滚动体直径d7.94 mm滚动体个数Z9接触角α0°。先用轴承特征频率公式算出各故障频率作为后续判别的基准%% 轴承参数与特征频率计算 fs 25600; % 采样率 Hz fr 30; % 转频 Hz D 39.04e-3; % 节圆直径 m d 7.94e-3; % 滚动体直径 m Z 9; % 滚动体数量 alpha 0; % 接触角 rad BPFO Z * fr / 2 * (1 - d/D * cos(alpha)); % 外圈故障特征频率 BPFI Z * fr / 2 * (1 d/D * cos(alpha)); % 内圈故障特征频率 BSF D / (2*d) * fr * (1 - (d/D * cos(alpha))^2); % 滚动体故障特征频率 FTF fr / 2 * (1 - d/D * cos(alpha)); % 保持架故障特征频率代入数值后得到BPFO约107.5 HzBPFI约162.5 HzBSF约70.2 HzFTF约11.9 Hz。这些频率在后面对照包络谱谱线时非常关键尤其是BPFO和BPFI它们之间相差大约55 Hz与转频的调制关系也能帮助我们区分内圈还是外圈故障。3.2 加载信号并通过快速谱峭度定位最优频带读入振动数据后我习惯先绘制时域波形的粗览图确认数据中没有明显的仪器饱和、断点或异常冲击再进入谱峭度计算。这样能避免把传感器敲击、噪声脉冲误当成轴承故障信号。%% 加载信号示例代码按实际文件格式调整 load(bearing_fault.mat); % 内含 x 和 fs t (0:length(x)-1) / fs; %% 快速谱峭度计算 level 6; [K, level_f, freq_f] fast_kurtogram(x, fs, level); %% 找到最优频带 [Kmax, idx_max] max(K(:)); [lv_opt, f_opt] ind2sub(size(K), idx_max); % 根据level_f和freq_f矩阵换算最优频带的中心频率和带宽 bw_opt fs / 2^(level_f(lv_opt, f_opt) 1); % 简化估算具体需按工具箱输出换算 fc_opt freq_f(lv_opt, f_opt); % 最优频带中心频率 f_low fc_opt - bw_opt/2; f_high fc_opt bw_opt/2; fprintf(最优频带: %.1f Hz ~ %.1f Hz, 峭度: %.2f\n, f_low, f_high, Kmax);这里有一个容易犯的坑level_f和freq_f在不同版本工具箱中的结构略有差异有的记录的是子带边界频率有的记录中心频率直接套用索引可能得到错误频带。我自己的做法是先打印几个频带的freq_f数值与色图横轴对比确认再决定如何换算带宽。如果嫌麻烦也可以用kurtogram函数直接输出Kurtogram图图上标注的最优频带参数是经过换算的可以直接读取。3.3 用带通滤波器提取共振信号并用Hilbert求包络谱得到最优频带后下一步就是对原始信号做带通滤波、提取包络、计算包络谱。滤波器我推荐用Butterworth加filtfilt零相位滤波避免普通滤波引起的相位失真干扰冲击位置判断。滤波阶数4到6阶足够太高会把共振频带削出振铃太低则阻带衰减不足。%% 带通滤波截取最优共振频带 ord 4; [b, a] butter(ord, [f_low, f_high]/(fs/2), bandpass); x_f filtfilt(b, a, x); %% 希尔伯特变换提取包络 env abs(hilbert(x_f)); %% 包络谱计算 N length(env); env env - mean(env); % 去除直流分量 EnvSpec abs(fft(env)); f_axis (0:N/2-1) * fs / N; EnvSpec EnvSpec(1:N/2); %% 绘制包络谱与特征频率标注 figure; plot(f_axis, EnvSpec); xlim([0, 500]); hold on; for k 1:4 xline(k*BPFO, --r, sprintf(%dBPFO, k)); xline(k*BPFI, --g, sprintf(%dBPFI, k)); end xlabel(频率 (Hz)); ylabel(幅值); title(最优频带包络谱); legend(包络谱, BPFO谐波, BPFI谐波); grid on;运行后如果故障特征明显包络谱中会在BPFO及其2倍、3倍频处出现清晰谱峰对应的就是我们最关心的故障指示。这里强调一个操作细节包络谱的纵轴是线性幅值还是对数幅值会影响对低频小幅值谱线的判断。我一般先用线性坐标看主峰再用对数坐标观察微弱的高次谐波两者配合能更全面评估故障严重程度。3.4 没有工具箱时的替代方案分频带峭度扫描有的情况下不方便下载第三方工具箱此时可以自写一个简化的分频带峭度扫描程序作为快速谱峭度的轻量替代。思路不复杂把整个频率范围均匀划分成若干候选频带用带通滤波器逐个提取后计算时域峭度峭度最高的频带就是近似最优频带。这个方法的频率分辨率不如1/3二叉树精细但对于早期故障已经足够。%% 简化版分频带峭度扫描 nband 48; % 候选频带数量可调 bw fs / nband; % 每个频带的宽度 kurt_band zeros(1, nband); for k 1:nband fL (k-1) * bw; fH k * bw; [bb, aa] butter(6, [fL, fH]/(fs/2), bandpass); y filtfilt(bb, aa, x); kurt_band(k) kurtosis(y); end [~, best] max(kurt_band); fL_best (best-1) * bw; fH_best best * bw; fprintf(简化扫描最优频带: %.0f ~ %.0f Hz\n, fL_best, fH_best);但这个简化方法有一个明显局限频带宽度是固定的如果真实共振频带比划分带宽窄或宽峭度值会被平均稀释。所以我的建议是先用粗划分比如32个频带看大致位置再用更细的划分在候选位置附近加密扫描结合包络谱效果综合判断。4. 谱峭度包络谱联合诊断从谱线到故障类型的判别方法4.1 特征频率与转频的调制关系拿到包络谱后如何判断是外圈、内圈还是滚动体故障最直接的依据是看谱峰是否出现在对应的理论特征频率及其谐波处。但在实际谱图中还需要关注以下几点外圈故障的特征是BPFO及其谐波处有谱峰且谱峰处通常伴随转频fr的边带。由于外圈固定安装于轴承座载荷分布相对稳定所以BPFO频段的谱线通常尖锐、能量集中边带相对较弱。内圈故障则不同。内圈随转子一起旋转故障点周期性进入和退出承载区导致冲击幅度周期性变化因此BPFI处不仅有谱峰其两侧还会出现明显的转频fr及其倍频调制边带。换句话说“BPFI谱峰带fr边带”是内圈故障的典型指纹。滚动体故障的包络谱特征相对复杂。由于滚动体自转频率较高且故障点有时与内圈接触、有时与外圈接触冲击幅度不稳定BSF处谱峰往往较矮且伴随保持架特征频率FTF或转频的边带。实际诊断中BSF附近噪声偏大需要结合更高阶谐波综合判断。保持架故障的特征频率FTF通常低于10 Hz且幅值很小容易被频谱低频段的高能量漏过掩盖因此包络谱分析对保持架故障的灵敏度有限这里不做重点。4.2 内圈故障和外圈故障的判别实例我拿一组典型的内圈故障数据举例。谱峭度选出的最优频带在4.5kHz到7kHz附近该频带包络谱在162.5 Hz处出现绝对主峰同时在107.5 Hz处几乎看不到谱线。乍一看像是单纯的内圈故障BPFI但仔细观察162.5 Hz主峰两侧能看到间隔约30 Hz的一对边带也就是132.5 Hz和192.5 Hz。转频是30 Hz所以这个模式与“内圈故障、周期性进入承载区调制”完全吻合判断为内圈故障无疑。换一组外圈故障数据包络谱主峰出现在107.5 Hz处而且2倍频215 Hz、3倍频322.5 Hz处都有可辨识的谱峰主峰两侧的边带减弱到几乎不可见。即使不看时域波形和包络波形单凭这个谱图也能做出外圈故障的判断。我的体会是特征频率绝对值重要但谱峰之间的“相对关系”——谐波序列是否完整、边带间距是否等于转频、主峰与边带的幅值比——才是区分故障类型的关键。只看单一谱线位置很容易误判。4.3 多频带对比为什么只看最优频带不够工程中我常遇到一种情况Kurtogram色图上不只一个高峭度区域最优频带包络谱很干净但次优频带包络谱能揭示更多信息。举一个具体例子某次分析中算法选出的最优频带在6kHz附近包络谱显示出明显的BPFO峰但另一个峭度次高的频带在2.5kHz附近包络谱中除了BPFO还出现了明显的啮合频率边带。后查得知该设备齿轮箱与轴承座之间存在结构耦合齿轮故障的冲击也调制到了这一频带。这意味着谱峭度选出的最优频带反映的是“冲击能量最集中的频带”不一定是“唯一包含诊断信息的频带”。因此我的做法是在自动推荐的基础上对Kurtogram色图中峭度前3到5名的频带都做一次包络谱把多频带包络谱并排显示对比。如果所有高峭度频带都指向同一故障频率结论可靠性会大幅提高如果不同频带指向不同故障则应警惕多重故障或来自其他部件的干扰。5. 实测中的参数陷阱与调参经验那些会让结果“翻车”的细节5.1 数据长度不足导致频率分辨率不够包络谱的频率分辨率取决于包络信号时长Δf fs/(N_band)其中N_band是包络信号长度。假设信号时长只有0.25秒那么包络谱频率分辨率约4 Hz。对转频30Hz、BPFO约107.5Hz的轴承来说4 Hz的分辨率勉强够用但如果故障特征频率之间相差很小比如BSF与某个边带仅差2 Hz就完全无法区分。我给初学者的建议是至少采集20到50个故障周期。外圈故障周期为1/BPFO约0.0093秒50个周期约0.47秒所以1秒以上的数据对BPFO级别的故障足够但如果要观察转频边带尤其是保持架故障FTF约12Hz则需要至少5到10秒数据才能清晰分辨。压箱底的经验是故障诊断宁可多采几秒数据也不要为了节省存储空间缩短采样时长后期发现分辨率不够再补采成本更高。5.2 转速波动对谱线模糊化的影响如果设备在分析期间存在明显转速波动故障冲击的重复周期不再严格恒定包络谱中的谱峰会被“抹开”幅值降低、带宽变宽甚至难以识别。这对谱峭度算法本身是致命打击——它假设故障频率稳定一旦转速漂移同一频带上计算的峭度值会被多段不同频率的冲击平均掉导致最优频带选择失准。遇到转速波动工况我会优先使用阶比跟踪或重采样技术先通过键相脉冲把信号从时间域映射到角度域再做谱峭度和包络分析。如果没有键相信号至少应检查一段数据的瞬时转速比如用转速计或基于转频的STFT估计如果波动超过1%就要谨慎解释包络谱结果。5.3 最优频带接近Nyquist频率时的“假阳性”快速谱峭度有时会把最优频带选在接近fs/2的高频区。这时候要警惕因为它可能不是轴承共振而是信号采集链路的伪迹或结构噪声。例如传感器的固有谐振峰值较高或电缆屏蔽不佳引入的电磁干扰都会产生高频冲击性成分拉高高频段峭度。排查方法很简单先看时域波形如果“冲击”出现的时间间隔完全随机、与轴转周期毫无关联大概率不是轴承故障再看原始频谱中该频带附近是否有工频及其谐波如果有应考虑电源干扰。我处理过一批现场数据谱峭度把最优频带定位到12kHz包络谱却没有任何轴承特征频率反而出现了100Hz电源纹波的谐波最终确认是传感器线缆屏蔽层破损。这类假阳性在实验室里不常见但在现场非常常见。5.4 带通滤波参数不当导致的信号失真很多人在选择带通滤波器时随手用bandpass函数默认阶数不够高结果带外噪声泄漏严重也有人把阶数拉到20以上结果滤波后的信号出现明显的振铃现象包络谱中产生大量虚假伪峰。我推荐的原则是Butterworth滤波器阶数4~8阶使用filtfilt做零相位滤波。中心频带较窄时适当提高阶数到8较宽时4到6阶足够。滤波完成后必须对比滤波前后的时域波形确保冲击序列依旧保持原有的周期性。如果发现滤波后冲击之间出现了额外的“尾巴”或振荡多半是阶数过高或频带边界选择不当需要调整。5.5 峭度值下降不等于故障消除晚期故障的识别问题一个鲜为人知的坑是随着故障从早期到中晚期发展时域峭度和谱峭度反而可能下降。原因是故障严重时冲击变成连续的大面积接触不再是单一“尖峰”同时多个冲击叠加后的分布趋于平滑峭度回落到接近正常水平。也就是说如果你只在峭度值最高的时段采集数据可能恰好错过最紧急的报警窗口。这就是为什么我不建议单独依赖谱峭度做趋势监测而应该同时跟踪包络谱特征频率处的幅值、时域峭度、有效值RMS等多维指标。实际项目中我习惯记录每个时刻的最优频带峭度值和对应特征频率幅值把它们存入历史数据库观察趋势变化。如果峭度值下降但BPFO幅值仍在上升说明故障在恶化但冲击性减弱此时更应保持警惕。6. 程序扩展思路从单段信号诊断到批量状态监测6.1 批处理脚本一次处理几十个数据文件实际项目里很少只分析一段信号最常见的需求是把一批文件全部跑一遍输出一个汇总表。我的做法是写一个批处理函数外层用dir遍历文件夹内层调用前面实现的单段分析函数最后把故障频率幅值、最优频带、峭度值写入Excel。files dir(data/*.mat); results table(); for i 1:length(files) load(fullfile(files(i).folder, files(i).name)); % 调用谱峭度包络谱分析函数 [bpfo_amp, bfi_amp, kurt_max, f_low, f_high] bearing_diag(x, fs, fr, D, d, Z); results [results; table({files(i).name}, bpfo_amp, bfi_amp, kurt_max, f_low, f_high)]; end writetable(results, diagnosis_results.xlsx);批处理的关键在于容错。个别文件可能因传感器松动、数据损坏等原因产生异常波形如果不加保护整个循环会中断前面的结果全部丢失。所以我的循环里必加try-catch异常文件单独记录继续处理后续文件。6.2 与包络谱结合构建健康趋势指标单次诊断只能回答“现在是否故障”回答不了“故障如何发展”。要构建趋势需要把每次分析得到的特征频率幅值、峭度、RMS等指标串联成时间序列。以BPFO幅值为例健康状态时几乎为零早期故障出现后开始升高中期可能保持稳定或波动严重时达到峰值甚至下降。以RMS为例故障初期变化不大中后期快速上升。实际监测项目里我会把多个指标放进同一个趋势图中用不同颜色区分并叠加报警阈值。阈值设定方法有两种一是基于健康样本统计的3σ准则二是基于行业标准经验值。前者更客观但前提是故障前有多组健康数据后者更省事但需要根据具体设备调整。6.3 与其他诊断方法的融合谱峭度包络谱是非常有效的方法但并非万能。现场诊断时我还会配合时域统计指标RMS、峭度、峰值因子、频段能量比、倒频谱、深度学习的频谱图像识别等方法交叉验证。谱峭度的最大价值在于“降维”——把宽频段的搜索问题变成少数几个候选频带大幅节省人工试错时间。而包络谱的价值在于“精判”——在选定的频带内给出物理意义明确的特征频率谱线。两者结合构成了我目前最顺手的一套分析流程。已跑通的这套MATLAB程序我已经在实验室故障注入数据和现场采集数据上反复验证过。我自己优化后的流程是先用Kurtogram快速定位候选频带再用多频带包络谱交叉验证最后结合时域特征核对。整个过程从加载数据到输出诊断结论通常不到5分钟。最后分享一个小技巧把fast_kurtogram的色图和包络谱图保存在同一个PDF中诊断报告就能非常直观地向非专业人员展示“为什么在这个频带分析”“故障频率在哪里”——这比单纯贴一张频谱图有说服力得多。
