简介这是一套面向机械故障诊断与信号处理研究者的快速谱峭度Kurtogram算法工具包适合开展轴承、齿轮等旋转机械故障特征提取与早期预警分析。资源以Matlab脚本.m和实验数据.mat为核心共24个文件压缩包大小仅4.64MB便于快速下载与运行。其中包含Fast_Kurtogram、Find_wav_kurt、Find_stft_kurt等核心算法实现以及Kurtogram可视化与多种谱峭度计算函数并配有不同工况下的轴承振动数据如OR007、NQX211、GGX224等可用于复现典型故障案例、验证算法效果。已有1155人学习下载适合具备一定信号处理基础、希望将峭度谱方法落地到实际数据集上的工程师和研究人员。通过研读这些相互调用的脚本可系统理解快速谱峭度的计算流程、频带选择逻辑及故障特征定位思路为后续自定义诊断算法提供可直接修改的代码框架。1. 快速谱峭度故障诊断里先“定位”再“解调”的那个前置步骤给滚动轴承做故障诊断时最常见的困境不是“有没有故障”而是“故障冲击藏在哪个频段”。转频、齿轮啮合频率、结构共振和随机噪声叠在一起直接看频谱或做包络分析特征频率往往被淹没。快速谱峭度Fast Kurtogram解决的就是这个问题它把信号按频带切分计算每个频带的峭度找到瞬态冲击最集中的共振带再把带通滤波后的信号交给包络谱做解调。Kurtogram的横轴是频率纵轴是分解层数颜色深浅代表谱峭度值最大点对应的中心频率和带宽就是后续诊断的入口。这篇文章面向做机械故障诊断、状态监测和PHM的工程师也适合刚接触轴承故障诊断入门的研究生。我会从谱峭度的定义、快速算法的分解策略讲起给出可在 MATLAB 里直接跑通的最小实现再结合仿真数据演示“Kurtogram 定带 → 带通滤波 → 包络谱解调”的完整链路。最后落在一个常被忽略的验证技巧上如何确认你选中的共振带不是偶然的峰值。2. 从峭度到 Kurtogram谱峭度为什么能定位故障共振带2.1 时域峭度的局限与谱峭度的定义2.1.1 时域峭度只能告诉你“有没有”不能告诉你“在哪”峭度是四阶标准化矩度量信号分布相对于高斯分布的尖锐程度。对滚动轴承而言外圈或内圈出现局部剥落时滚珠经过缺陷位置会产生周期性冲击时域波形上出现稀疏的瞬态尖峰峭度值明显升高。这也是很多诊断系统把时域峭度作为报警指标的原因。但时域峭度有个致命缺陷它把整个频段的信息揉成了一个标量。背景噪声强一点、齿轮振动大一点冲击就会被稀释峭度值可能只在故障早期短暂抬升随后下降甚至低于正常值。更关键的是即使峭度升高了你依然不知道下一步该对哪个频段做带通滤波。峭度是“有没有异常”的粗筛不是“异常在哪”的定位工具。2.1.2 谱峭度把峭度摊到每个频率上谱峭度的思路很简单先用短时傅里叶变换把时域信号切成时间-频率平面然后在每个频率点 f 上沿着时间方向计算规范化四阶矩。定义如下K_x(f) E{|X(t,f)|^4} / (E{|X(t,f)|^2})^2 - 2式中 X(t,f) 是信号 x(t) 的短时傅里叶变换结果E 表示对时间帧做平均。减 2 是让平稳高斯噪声在该频率上的谱峭度值为 0这样谱峭度图中所有显著大于 0 的频带都对应“偏离高斯”的成分——也就是瞬态冲击所在的位置。这里面的物理含义值得多说一句。轴承故障冲击不是平稳信号它的能量在频域上不连续而是激发结构共振形成一段局部放大的频带。这段频带的幅值随时间呈衰减震荡波形统计分布明显非高斯因此谱峭度值高。而连续旋转的转频成分、齿轮啮合成分随时间变化平稳谱峭度接近 0。所以谱峭度天然是“瞬态冲击探测器”和故障诊断的需求高度契合。2.2 Fast Kurtogram用 1/3-二叉树把计算量降下来2.2.1 全 STFT 谱峭度为什么算不动直接按定义计算谱峭度有个工程问题要得到足够的频率分辨率短时傅里叶变换的窗长必须足够长而频率点数一多每个频率点都要对全部时间帧求四阶矩计算量迅速膨胀。工业振动数据动辄几十万点设备在线监测又要求快速出结果全谱峭度的实时性不达标。Antoni 在 2007 年提出了快速谱峭度算法核心思想是用滤波器组对信号做逐级二分和三分分解用不同带宽的子带信号替代短时傅里叶变换的频率切片。每个子带的谱峭度计算都是直接在时域信号上做的先带通滤波再分帧估计四阶矩。这样频率分辨率和计算复杂度之间不再互相绑定整体效率提升几个数量级。2.2.2 1/3-二叉树与频带划分Fast Kurtogram 的分解结构是把频带按“先粗后细”的方式逐步切分。第一层把整个奈奎斯特频带二等分第二层有两个选择既可以把子带再二等分也可以按三分之一处切分然后把这两种切分的结果合并成一个统一的层级-频率网格。说白了就是在频率分辨率和估计精度之间做折中带宽窄的子带频率定位准但每个子带内包含的瞬态事件数少峭度估计方差大带宽宽的频带峭度估计稳定但中心频率定位粗糙。Kurtogram 图就是这个网格的可视化纵轴是分解层数每往下一层代表带宽减半或减为三分之一横轴是中心频率颜色是谱峭度值。算法最终返回网格中的最大值点包含三个关键参数中心频率 fc、带宽 Bw、谱峭度值 Kv。这三个参数就是故障定位的核心输出直接喂给后续的带通滤波器。2.3 谱峭度、包络谱、功率谱在故障诊断里的分工工具看的对象输出局限功率谱平稳周期成分转频、啮合频率及谐波瞬态冲击能量被平均早期故障看不清包络谱调幅信号解调后的故障特征频率必须预先知道带通频率否则结果无意义谱峭度 / Kurtogram非平稳瞬态共振带中心频率、带宽、峭度值只定位不含“诊断结论”需与包络谱联动在故障诊断流程里Kurtogram 通常不出最终结论它负责定频带包络谱负责定频率。两者组合起来才能在噪声背景下识别出故障特征频率。这也是“Kurtogram 包络谱”被写进大量轴承故障诊断论文里的原因。3. 用 MATLAB 跑通快速谱峭度的最小实现3.1 先造一个带冲击的仿真信号没有合适实验数据时最好先构造已知故障特征的仿真信号这样能精确验证算法输出是否和设定一致。下面这段代码生成包含重复瞬态冲击、转频正弦和噪声的仿真振动信号% 仿真轴承外圈故障振动信号 fs 20000; % 采样率 20 kHz t (0:0.5*fs-1)/fs; % 0.5 秒信号 fr 25; % 转频 25 Hz bpf 130; % 外圈故障特征频率 130 Hz % 转频分量转子不平衡产生的平稳正弦 x 0.8 * sin(2*pi*fr*t); % 周期性冲击序列衰减正弦模拟故障激励共振 for k 1:60 ti (k-1) / bpf; % 冲击发生时刻 idx find(t ti, 1); % 对应采样点索引 idx_range idx:min(idx250, length(t)); x(idx_range) x(idx_range) 3 * exp(-2000*t(1:length(idx_range))) .* sin(2*pi*3500*t(1:length(idx_range))); end % 叠加高斯白噪声 x x 0.6 * randn(size(t));这段代码的关键参数有三个冲击周期1/bpf决定了故障特征频率3500 Hz 是结构共振频率指数衰减系数exp(-2000*t)模拟了冲击的能量衰减速度。实际运行后Kurtogram 应当把共振频带定位在 3500 Hz 附近而不能是噪声所在的低频段或高频段这是后续验证算法正确性的基准。3.2 自写 Fast Kurtogram 核心循环直接调用现成工具箱当然可以但为了理解参数含义有必要从零写一个精简版实现。核心逻辑是递归二分滤波频带并计算各层峭度function [fc, bw, kv] fk_demo(x, nlevel, overlap) % 简化版快速谱峭度二分滤波器组 滑窗峭度 % 输入x 信号列向量nlevel 分解层数overlap 帧信号重叠率 % 输出最大峭度对应的中心频率 fc、带宽 bw、峭度值 kv fs 20000; % 采样率实际使用中作为参数传入更合理 N length(x); nband 2^nlevel; % 最终子带数 bw fs / 2 / nband; % 每层滤波器带宽 kurt zeros(nband, 1); % 每个频带的峭度值 for m 0:nband-1 % 设计带通滤波器用 6 阶 Butterworth过渡带由频率分辨率决定 fl m * bw / (fs/2); % 归一化低频截止 fh (m1) * bw / (fs/2); % 归一化高频截止 [b, a] butter(6, [fl fh], bandpass); y filter(b, a, x); % 带通滤波 % 分帧帧长取基频周期的 4 倍重叠率默认 0.5 flen round(fs / bw * 0.05); % 原则每帧至少包含 4-6 个冲击周期 step round(flen * (1 - overlap)); frames buffer(y, flen, flen - step, nodelay); % 计算每帧能量和四阶矩然后按时间平均 seg_energy sum(frames.^2, 1); seg_power2 sum(frames.^4, 1); kurt(m1) mean(seg_power2) / mean(seg_energy).^2 - 2; end % 取峭度最大值对应的频带 [kv, idx] max(kurt); fc (idx - 0.5) * bw; end代码逻辑分三步先做带通滤波再分帧估计能量和高阶矩最后按时间平均得到该子带的谱峭度。几个参数需要说明butter(6, …)用 6 阶巴特沃斯滤波器阶数太低则子频带间泄漏严重峭度值被污染帧长flen的选取原则是每帧至少包含 4-6 个冲击周期太长则冲击被平摊失去意义太短则四阶矩估计方差过大- 2的修正和谱峭度定义保持一致让平稳噪声段接近 0。实际生产环境里这个二分法只是 Fast Kurtogram 的简化版真正的算法按 1/3-二叉树同时做二等分和三等分在每一层比较两种切分方式哪个峭度更高。但对于理解“滤波器组 分帧 四阶矩”这套机制二分版本足够直观调参逻辑也完全兼容。3.3 Kurtogram 关键参数表参数默认值调大调小分解层数 nlevel4频率定位更细但峭度估计方差大更快适合低采样率或在线实时帧重叠率 overlap0.5峭度估计平滑计算量增大速度快但方差上升滤波器阶数6带外抑制好计算稍慢过渡带宽相邻频带串扰大帧长 flen20 ms对低频冲击友好时间分辨率变差高频瞬态定位准但统计不稳定参数之间是互相牵连的。nlevel 不是越大越好层数过深到带宽小于故障特征频率的 2 倍时一个冲击就会被切到两个相邻频带里导致两个子带的峭度同时下降图案出现“断层”或“棋盘格”。发现这种情况时优先减小 nlevel不要盲目追求频率分辨率。4. 快速峭度谱实战定带、滤波、包络谱解调全流程4.1 从 Kurtogram 输出到包络谱识别故障频率拿到 Kurtogram 的输出[fc, bw]后完整的诊断链路分为四步带通滤波、Hilbert 变换求包络、包络谱分析、对照故障特征频率表。% 输入x 原始振动信号fc 和 bw 来自 fk_demo 输出 fl max(0, (fc - bw/2)) / (fs/2); fh min(1, (fc bw/2)) / (fs/2); [b, a] butter(4, [fl fh], bandpass); y_band filter(b, a, x); % 带通后的共振信号 % Hilbert 包络解调 env abs(hilbert(y_band)); % 解析信号取模得到包络 % 包络谱 N length(env); w hann(N, periodic); env_fft fft(env .* w); f_axis (0:N/2-1) * fs / N; env_spec 2 * abs(env_fft(1:N/2)) / sum(w); % 找包络谱前 5 个峰值并输出 [pks, locs] findpeaks(env_spec, MinPeakHeight, max(env_spec)*0.3, MinPeakDistance, 5); [~, idx] sort(pks, descend); top5 sort(f_axis(locs(idx(1:min(5, length(idx))))));带通滤波这里用的 4 阶巴特沃斯比峰值定位阶段低原因是包络谱对带外残留噪声敏感带通越陡峭包络谱越干净。但阶数过高会引入相位畸变工程上 4-6 阶是折中区间。包络谱的峰值位置要对照轴承参数计算的理论故障特征频率来判断。仿真信号中理论 BPFO 是 130 Hz那么包络谱在 130 Hz、260 Hz、390 Hz 处应出现衰减的峰族。如果峰值和理论值偏差在 1% 以内基本可以确认该轴承外圈存在局部缺陷。4.2 共振带不明显时的调参策略实际工业数据往往比仿真信号复杂。齿轮啮合频率及其边带也会产生调幅效应形成伪瞬态和轴承冲击竞争高峭度频带。这种情况俗称为“齿轮淹没”或“调制干扰”。常规做法是看 Kurtogram 的前几个局部极大值而只锁定全局最大点。具体实现上把fk_demo输出的所有频带峭度值排序取前三个局部最大峰对应的[fc, bw]分别做带通和包络谱解调观察哪个频带解调出的特征频率与轴承故障频率吻合。齿轮啮合带来的调制边带通常对应分辨率低的宽频带而轴承冲击对应窄的共振峰两者在带宽上会露出马脚。经验法则是优先选带宽较窄且峭度值超过背景均值 3 倍以上的频带。另一个常见问题是冲击重复频率恰好和某个转子谐波成分接近导致包络谱里真假峰值混叠。此时需要提高包络谱的频率分辨率把原始信号截断为更长的帧或者对包络做 rms 平滑后再做 FFT。注意不要仅凭一次分析下结论最好对比两个不同转速工况下的包络谱故障特征频率与转频之比在变转速下保持不变而干扰频率会随转速等比漂移。4.3 Kurtogram 排错清单现象可能原因应对颜色图平坦没有突出峰信号噪声过大或冲击太弱增加抗混叠预处理降低 nlevel 并用更长帧最大峭度频带在采样率上限附近高频电噪声或虚假瞬态先做低通抗混叠滤波降到 10 kHz 再分析相邻层峭度值剧烈跳变帧长太短四阶矩估计不稳定增大 overlap 到 0.75 或帧长翻倍频带窄到包络谱无周期峰值nlevel 过深撕裂冲击序列减小 nlevel确保带宽大于 3 倍故障特征频率两个工况下最大频带漂移很大转速波动频率不再是准平稳切换到阶次跟踪或同步平均后再做 Kurtogram最容易踩的坑是在原始信号上直接调用 Kurtogram 而不做任何预处理。实际工业测点往往叠加了 50 Hz 工频干扰和结构固有频率造成的非故障冲击这些成分同样会在 Kurtogram 上形成高峭度区域。建议把变速器壳体振动、油膜涡动等干扰源单独建模排除后再对残差信号做谱峭度分析。5. 最后留一手用分块重计算验证 Kurtogram 选带的稳定性前面提到Kurtogram 最大值本身只是一个点估计它受帧长、重叠率、分解层数的共同影响。同一个信号稍微改动重叠率最大峭度对应的中心频率就可能跳几百赫兹。这在工程上是不可接受的你按这个频带设计了带通滤波器在线跑几个月结果发现初始选带有偏差。我常用的验证方法叫“分块稳定性检验”。把原始信号切成互不重叠的四段每段分别调用fk_demo得到各自的最大峭度中心频率。如果四个结果的中心频率偏移量小于一个带宽说明选带稳定可以放心用于在线监测如果偏移量超过两倍带宽说明该信号存在明显的非平稳性要么时变转速要么干扰源随机出现这时单次 Kurtogram 结果不可靠需要改用角度域重采样后的阶比 Kurtogram。% 分块稳定性验证四段独立计算中心频率 seg_len floor(N/4); fc_list zeros(4, 1); for k 1:4 seg x((k-1)*seg_len1 : k*seg_len); [fc_list(k), ~, ~] fk_demo(seg, 4, 0.5); end fc_span max(fc_list) - min(fc_list); bw_est fs / 2 / 2^4; % nlevel4 时的子带带宽 stable fc_span bw_est; % true 表示选带稳定分块之外还有更严格的双重验证拿走 10% 的随机样本重算一次统计中心频率的置信区间。实现上就是二次重采样重复 20 次看中心频率落在哪个带宽范围内的比例。占比高于 80%该频带可以写入设备档案低于 50%就换用多个子带同时监测把五个最高峭度频带全部接入在线诊断逻辑稍后让特征频率识别来仲裁它们各自的输出。这层验证存在的意义是快速谱峭度给出的不是一个“最优带”而是一个“候选带”。在设备和工况都相对稳定的前提下分块重算的结果应当保持一致一旦出现不一致优先怀疑不是算法问题而是信号本身的时变性。带着这个筛选结果再去设计带通滤波和包络谱阈值设备诊断系统的误报率会低很多。本文还有配套的精品资源点击获取
