ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

快速谱峭度(kurtogram)实现轴承故障诊断:从峭度到包络谱

快速谱峭度(kurtogram)实现轴承故障诊断:从峭度到包络谱 简介面向机械故障诊断与信号处理场景的快速谱峭度工具包以Kurtogram算法为核心提供完整Matlab实现与配套实验数据。内容覆盖峭度谱计算、快速峭度图绘制、谱峭度特征提取等关键环节适用于设备状态监测、滚动轴承故障识别及异常冲击特征提取等任务。资源共24个文件以17个.m源码脚本为主含Fast_Kurtogram.m、demo_Fast_Kurtogram.m、Find_stft_kurt.m等核心程序另有7个.mat轴承振动数据文件如OR0073_1、GGX224_DE_time等便于直接导入运行验证。压缩包整体仅4.64MB文件体量小适合需要快速复现快速峭度算法的工程师或初学者。已有1155人学习下载。通过源码梳理与数据实测可快速理解Kurtogram参数选择、峭度谱可视化与故障频带定位思路并迁移到自定义信号分析流程中提高故障诊断效率。1. 快速谱峭度与kurtogram轴承故障诊断为什么要先看峭度轴承早期故障的振动信号幅值未必有明显变化但波形会突然多出一串细小的冲击。峭度就是量化这种尖的程度正常振动接近高斯分布峭度约等于3出现剥落、压痕后每转一圈的冲击能把峭度推到10以上。麻烦的是冲击激发的结构共振往往在数千赫兹的高频段不在你平时习惯看的低频区。快速谱峭度kurtogram把峭度这个单一统计量扩展成中心频率—带宽—峭度的二维热图自动圈出峭度最大的共振频带。拿到这个频带再做包络谱轴承外圈、内圈、滚动体的特征频率才能清晰暴露出来。下面按原理、MATLAB实现、参数排错、选带验证四步展开代码直接改采样率和故障频率即可运行适合做故障诊断与设备状态监测的工程师参考。2. 峭度到快速峭度谱从四阶矩到kurtogram的分层结构2.1 峭度为什么对冲击敏感信号 x 的峭度定义为四阶中心矩与方差平方之比K E[(x − μ)⁴] / σ⁴对高斯白噪声K3对纯正弦K1.5对带冲击的振动K 明显大于 3。原因在于四阶矩对离群样本的放大权重远高于二阶矩一个样本偏离均值 d 倍对方差的贡献是 d²对峭度的贡献却是 d⁴。轴承故障产生的冲击在时域上正是少而尖的离群区所以峭度被公认为刻画冲击性最直接的统计量也常写作 excess kurtosisK−3正常信号应在 0 附近。但时域峭度是整个信号的一个标量它只能告诉你有没有冲击说不清冲击藏在哪个频带。实际轴承故障的冲击会激励起轴承座或传感器的固有共振共振频率高达几 kHz且随故障位置不同而变化。于是需要把峭度按频率展开这就引出了谱峭度。2.2 谱峭度每个频点上单独算峭度谱峭度的标准做法是对信号做短时傅里叶变换STFT得到时频矩阵 X(t, f) 后对每一根频率线沿时间轴求四阶矩与二阶矩之比SK(f) ⟨|X(t,f)|⁴⟩ / ⟨|X(t,f)|²⟩² − 2其中 ⟨·⟩ 表示对所有时间帧求平均。减 2 的作用是让平稳高斯噪声的谱峭度归零、纯正弦归到 −1这样谱峭度为正是这里存在脉冲性成分的明确标志。下面这个特性值得记住短时窗的窗长直接决定频率分辨率窗越长频率线越细但时间帧越少峭度估计的方差越大。这带来一个矛盾故障冲击的共振可能落在任意中心频率、任意带宽上。固定一种窗长只能扫描一种带宽扫不准就把冲击能量摊薄了峭度自然被压低。2.3 kurtogram 的分层结构解决带宽扫描问题kurtogram 的核心思路是对频带做树状划分第 k 层将 0fs/2 的整个频段切成 2^k 个等宽子带每个子带用与自身带宽匹配的窗长计算谱峭度得到一张以中心频率—带宽层数为坐标的峭度热图。层数越大带宽越窄能分辨的共振峰值越精确但每个子带内的统计方差也越大。实际工程里几乎不会用暴力遍历而是采用 Antoni 提出的快速算法用一组 1/3-二叉树结构的滤波器组逐层分解信号每个频带的峭度在分解过程中递推得到总计算量从 O(N²) 降到 O(N log N)这就是快速谱峭度名称的由来。MATLAB 里常见的 kurtogram 工具箱以及后续的快速峭度谱都基于这套结构输入一段时域信号输出以层数、频带为索引的峭度矩阵再换算成最优中心频率 fc 和带宽 bw。3. 用MATLAB跑通kurtogram快速谱峭度最小故障诊断代码3.1 先生成一串带共振的轴承故障信号写代码前先要有数据。模拟轴承外圈故障的经典模型是周期性冲击激励起高频共振再叠加噪声。% 模拟轴承外圈故障振动信号便于验证kurtogram选带 fs 25600; % 采样率 25.6 kHz t (0:fs-1)/fs; % 1 秒数据 fr 30; % 转频 30 Hz BPFO 5.4 * fr; % 外圈故障特征频率约162 Hz x zeros(size(t)); imp_idx 1 : round(fs/BPFO) : fs; % 每隔一个BPFO出现冲击 for k 1:length(imp_idx) n imp_idx(k) : min(imp_idx(k)127, fs); x(n) x(n) exp(-(n-imp_idx(k))/32) .* ... sin(2*pi*3200*(n-imp_idx(k))/fs); % 3200Hz共振衰减振荡 end x x 0.6*randn(size(t)); % 叠加强白噪声冲击间隔取 fs/BPFO每个冲击用指数衰减包络调制 3200 Hz 正弦模拟轴承座共振衰减常数 32 个采样点对应约 1.25 ms 衰减时间。共振频率请换成你设备实测的固有频率BPFO 换成按轴承参数算出的理论值。最后叠加的高斯白噪声把信噪比压到约 5 dB逼近早期故障的实际情况。3.2 调用 kurtogram 工具箱的最少命令用 Antoni 公开的 kurtogram 工具箱时计算和显示是分开的% 快速峭度谱计算与kurtogram可视化需先安装工具箱网上有公开版本 nlevel 4; % 分解层数一般取3~6 [K, fc, bw] kurtogram(x, fs, nlevel); % K:各层各带峭度值 tkurtogram(x, fs, nlevel) % 绘制kurtogram热图kurtogram 的输出 K 是 (nlevel1)×频带数的矩阵K(k,b) 表示第 k 层第 b 个频带的谱峭度fc 和 bw 按层存放各频带的中心频率与带宽不同版本返回结构略有差异用 whos 查看为准。最容易的做法是从热图上看到最亮区域所在层数与频带再到对应层的 fc、bw 数组里取值。上述故障代码在 fs25600、nlevel4 时通常能选中 30003400 Hz 附近的频带与你设定的共振频率吻合。3.3 不装工具箱也能跑一段简化的快速峭度谱实现有些环境装不了第三方工具箱或者你想弄清楚每一层到底做了什么。我自己留了一份教学版实现逐层切频带、逐带带通滤波后直接在时域算峭度。原理与 kurtogram 相同速度慢一个量级但选带结果一致。function kgram simple_kurtogram(x, fs, maxlevel) % 简化快速峭度谱逐层切带逐带滤波后计算峭度 % 输入: x 1D振动信号, fs 采样率, maxlevel 最大分解层数(建议3~6) % 输出: 表格列为峭度K、中心频率fc、带宽bw K []; fc []; bw []; for k 0:maxlevel nb 2^k; % 本层频带数 bw_k fs / (2^(k1)); % 每个子带带宽 for b 1:nb f_lo max((b-1)*bw_k, 1); % 频带下界避免0Hz边缘问题 f_hi min(b*bw_k, fs/2-1);% 频带上界避开Nyquist边缘 y bandpass(x, [f_lo, f_hi], fs, ... ImpulseResponse,fir,Steepness,0.85); K(end1) mean((y-mean(y)).^4) / (var(y)^2); % 峭度 fc(end1) (f_lof_hi)/2; bw(end1) bw_k; end end kgram table(K, fc, bw, VariableNames, {K,fc,bw}); end % 使用示例 res simple_kurtogram(x, fs, 4); [~, i] max(res.K); fc_opt res.fc(i); bw_opt res.bw(i);每层对 0fs/2 切成 2^k 段带宽按指数变窄bandpass 用 FIR 加 Steepness0.85 控制过渡带边缘频带上下限做了 1 Hz 与 fs/2−1 的截断避免 MATLAB 对边界频率报错。这个版本每层都重新滤波复杂度近似 O(层数×带宽数×N log N)1 秒数据、4 层约几秒可接受。信号明显非平稳启停机、变转速时这个简化版的估计方差会变大需要先按稳态段分窗处理。4. 快速峭度谱的工程参数与排错层数、带宽与噪声干扰4.1 影响选带结果的四个关键参数参数建议取值对结果的影响分解层数 nlevel / maxlevel3~6层数小带宽粗共振被抹平层数大峭度方差大易选到孤立窄噪声带信号长度≥20 个冲击周期帧数太少时峭度估计方差急剧增大选带不稳定带通滤波器 Steepness0.8~0.9过渡带过陡引入振铃Gibbs现象在峭度上制造假峰采样率与共振频率比fs ≥ 5×共振频率保证共振频带内有足够谱线选带结果可重复这四个参数里最容易出问题的是层数和信号长度。层数超过 6 后最窄子带可能只包含几十根谱线纯噪声段会因样本量不足而峭度虚高直接误导选带。信号长度不足时一整段数据只有十几次冲击任何一次幅值抖动都会显著改变峭度这就是同一组数据前后两次计算选带结果不一样的最常见原因。4.2 峭度峰值不等于最优频带的三个场景第一个场景是信号里混有单个野值比如电磁干扰脉冲或敲击。野值会瞬间拉高时域峭度kurtogram 会把该冲击所在的宽带整段判为最优。排查办法是把原始波形画出来真故障冲击是周期性的、宽度一致野值是单发的、与转频无关。第二个场景是强谐波成分会把谱峭度压低。平稳正弦的谱峭度趋近 −1齿轮啮合、轴频分量越强的频带峭度越低于是 kurtogram 的亮点容易跑到能量很低的噪声富集区。选带前先看该频带的原始功率谱总能量能量不足全频段峰值 1% 的频带即使峭度再高也建议放弃这是快速谱峭度最常见的误用点。第三个场景是转速波动。转速不稳时同一种故障的冲击间隔不断变化包络谱边带被抹平kurtogram 的峭度峰值也会漂移。处理方式是先用测速信号把数据分段每段内转速波动控制在 1% 以内再单独做快速峭度谱否则 MATLAB 里直接整段计算的结果会忽高忽低。4.3 与包络谱配合的正确姿势快速峭度谱只回答冲击在哪一段频率故障类型的判断要靠包络谱。选带、滤波、希尔伯特解调三步连在一起写% 对选中的最优带做带通滤波后取包络 y_bp bandpass(x, [fc_opt-bw_opt/2, fc_optbw_opt/2], fs, ... ImpulseResponse,fir,Steepness,0.85); env abs(hilbert(y_bp)); % 希尔伯特包络解调 N length(env); f_ax (0:N/2-1)/N*fs; A_env abs(fft(env)); A_env A_env(1:N/2); plot(f_ax, A_env); xlim([0 500]); grid on xlabel(频率 / Hz); ylabel(包络谱幅值)包络谱的横轴范围取到 5 倍转频就够重点核对 BPFO、BPFI、BSF 以及它们的 23 倍谐波。若包络谱能量集中在转频及其边带说明问题可能来自轴系不平衡或不对中而不是纯粹的轴承局部损伤。这一步是快速峭度谱最容易被跳过的环节也是把故障诊断误报成轴系故障的主要来源。5. 把快速峭度谱用进轴承故障诊断自动选带与包络谱验证最后这章给一个我常用的验证流程帮助轴承故障诊断入门阶段少走弯路。5.1 完整的诊断流程对一段实测振动数据按以下顺序处理截取转速稳定段建议 12 秒→ 去均值和趋势项 → 快速峭度谱选带 → 带通滤波 → 希尔伯特包络谱 → 对照轴承特征频率。前面各章的代码已经覆盖前五步最后一步只需按轴承几何参数算出理论频率外圈故障 BPFO n/2 · fr · (1 − d/D · cos α)内圈故障 BPFI n/2 · fr · (1 d/D · cos α)滚动体故障 BSF D/(2d) · fr · (1 − (d/D · cos α)²)其中 n 为滚动体数量d 为滚动体直径D 为节圆直径α 为接触角fr 为转频。5.2 两个落地的验证技巧第一个技巧是用滤波前后峭度对照确认选带有效。计算原始信号的峭度 K_raw 和带通后信号的峭度 K_bp若 K_bp 明显高于 K_raw比如从 4 提到 15 以上说明最优带确实集中了主要冲击能量如果几乎不变化就要回到 4.2 检查野值或噪声带误选。这个对照可以作为自动选带的置信度指标直接写进报警逻辑避免高位报警里混入大量假信号。第二个技巧是固定工况下对比健康与故障状态的 kurtogram。健康轴承的峭度谱峰值位置是随机的、层数不稳定故障发展后峰值会稳定收敛到某个共振频带。因此不要拿单次 kurtogram 的明亮区域直接当作故障证据而是连续采集多段数据看 fc_opt 和 bw_opt 的重复性——重复出现的位置才值得进入后续寿命预测模型的特征列表。建好基线后用fc_opt, bw_opt, 包络谱 BPFO 幅值三个量做趋势曲线比单纯盯时域峭度更能区分冲击变强和故障扩展两种情况这也是把 kurtogram 从离线分析升级成在线监控特征时最划算的一步。本文还有配套的精品资源点击获取
返回列表