1. 谱峭度到底解决了包络分析的什么痛点1.1 包络分析的成功与无奈共振频带怎么选做滚动轴承或齿轮箱故障诊断的兄弟对包络分析Envelope Analysis应该都不陌生。它的逻辑非常直观故障点在运转中反复撞击产生周期性冲击这些冲击会激励起系统的高频共振而这个高频共振又承载着低频的故障特征频率。包络分析的思路就是把高频共振信号剥出来整流、低通滤波之后得到包络信号再做FFT在包络谱里找故障特征频率及其谐波。这个流程用起来很顺但有一个老大难问题带通滤波器的中心频率和带宽到底怎么定早年大家都是靠经验扫频谱看看哪个频段有能量鼓包就猜一个范围。问题是故障早期信号极其微弱共振频带经常淹没在噪声里人工看频谱挑频带十次有八次会挑歪。挑窄了故障特征频带被切掉一半挑宽了带进一堆无关的噪声包络谱效果和直接做全频段包络没什么区别。这个先验参数难以确定的痛点才是谱峭度算法出现的真正原因。1.2 峭度这个指标为什么在故障诊断里地位特殊峭度Kurtosis衡量的是数据分布相对于正态分布的尾部厚重程度。平稳的机械运转信号其分布非常接近高斯分布而故障冲击是瞬态、稀疏、幅值远大于正常水平的尖峰它会把分布拖出明显的重尾。所以峭度值天然就是有没有冲击的晴雨表正常信号峭度在3附近出现明显冲击后峭度可以飙到几十甚至上百。但这里有个重点可能很多人没有认真想过原始信号的峭度虽然能反映冲击强度却失去了频率定位能力。你说峭度高可冲击在哪个频段呢原始时域峭度完全答不上来。而机械故障诊断的实际需求是既要检测到冲击又要定位到频带否则后续的无量纲指标趋势跟踪、频带特征提取都无从谈起。所以我们需要的是一个随频率变化的峭度——谱峭度Spectral KurtosisSK它把峭度这个全局标量展开成每个频点附近的局部峭度哪个频带冲击越突出那里的谱峭度值就越高。1.3 谱峭度的原始定义每个频点上的偏离高斯程度谱峭度的原始定义并不复杂。对信号进行短时傅里叶变换STFT得到时频复包络 (X(t,f))然后沿时间方向计算每个频点的四阶累积量归一化值[ K(f) \frac{\left\langle |X(t,f)|^4 \right\rangle}{\left\langle |X(t,f)|^2 \right\rangle^2} - 2 ]减掉2是为了让纯高斯信号复变量的峭度本征值为2归零。因此理想的白噪声谱峭度约等于0带故障冲击的频带谱峭度远大于0。搜索 (K(f)) 的最大值就能找到冲击特征最强的频带。理解这个公式的要害在于(X(t,f)) 不是单一频率的FFT结果而是短时窗内的局部复包络它本身就隐含了以 (f) 为中心、由窗函数决定带宽的局部性。也就是说谱峭度天然是第三维带宽的函数不同的时间窗长对应不同的频率分辨率结果完全不同。这个问题留到后面讲快速谱峭度时再展开它正是算法设计的核心钥匙。2. 从谱峭度到快速谱峭度计算量是被什么逼出来的2.1 原始谱峭度的工程瓶颈二维搜索太贵了如果只是算一条频带里的谱峭度那直接滤波再算峭度就行。但实际工程里我们不知道最优频带的中心频率和带宽于是就需要在中心频率×带宽的二维平面上做穷举搜索。每给定一组参数就要做一次带通滤波、取包络、算峭度这个计算量是灾难级的。Antoni在2006年前后系统地解决这个问题时引入了一个关键观察可以预先构建一个多分辨率滤波器组把信号一次性分解成一组覆盖不同中心频率、不同带宽的子带信号然后对每个子带并行计算峭度值。这样一来中心频率和带宽的搜索就变成查表而不是反复滤波。这个思想后来被整理成著名的Fast Kurtogram算法也被广泛集成进各大商用故障诊断软件里。2.2 快速谱峭度的核心思想用框架替代逐点扫描我在自己的项目里经常把它类比成按地图找金矿逐点扫描是开着探测器在整片山里打转而快速谱峭度相当于先把山按1/2、1/4、1/8……切成网格每个网格测一个样本然后用网格里数值最高的那块地定位矿脉。这个比喻里网格就是滤波器组数值就是该子带信号的峭度。整个算法的设计目标可以概括为一句话用尽量少的子带分割次数覆盖尽量广的频带-带宽组合并保证每个组合的峭度值统计稳定。2.3 一个关键设计为什么是1/3-二叉树而不是普通二分树这里有一个初学者最容易忽略、但恰恰是Antoni算法最精巧的地方。直观地想带宽细分应该用二分树第0层是整个频带第1层分成两个子带第2层分成四个子带…… 二分树确实能完整覆盖频带但有一个明显的分辨率缺陷它只能提供带宽逐级减半的等比序列。实际机械故障的共振频带落在哪根本不是按2的幂次排列的。一个中心频率在 (0.3F_s)、带宽只有 (0.1F_s) 的共振带如果用二分结构要么落在0.25~0.5这个粗子带里被拓宽到0.25带宽要么被进一步分割开总是不对劲。Antoni的改进是引入1/3-二叉树结构每个父带除了做平分还额外做两个分割位置偏移1/3的子带。这样既保留了二分树的速率又将每个频带按1/3的粒度交错覆盖带宽序列从(1/2^k)变成(1/2^k) 与 (1/(3 \cdot 2^{k-1}))交替出现频率分辨率明显更细腻。用工程语言说1/3-二叉树在同样的分解层数下提供了比二分树多大约50%的频带-带宽候选组合而计算复杂度只增加了一个常数倍。这是整个快速算法快而不糙的根本原因。3. 快速谱峭度的核心实现机制拆解3.1 滤波器组是怎么搭起来的快速谱峭度的实现通常是基于一组半带低通滤波器Half-band Lowpass Filter和高通滤波器递归构建的。以Antoni原论文的体系为准处理流程大致如下原始信号 (x(n)) 进入第0级覆盖频带 ([0, F_s/2])。第1级使用两个互补滤波器将信号分成低频带 ([0, F_s/4]) 和高频带 ([F_s/4, F_s/2])。第2级在每个子带上继续二分同时按1/3偏移生成额外子带。这里的实现技巧是不是真的对每个子带分别滤波而是用频移统一低通抽取的方式让所有子带共享同一个原型低通滤波器。这种统一原型频移的做法是快速实现的关键。如果每个子带都设计独立的FIR带通滤波器滤波器的设计和卷计算量会随层级成指数膨胀比如K7时可能产生上百个子带独立滤波的成本完全顶不住。共享原型滤波器之后每个子带本质上只是把原信号做一次复调制乘一个复指数再过同一个低通再抽取计算量大大降低。3.2 每个节点的峭度是怎么算的复包络是关键滤波器组输出的每个子带信号都是复值信号因为频移后包含了正负频谱记作 (c_k(n))。接下来的峭度计算要小心不能直接对实部、虚部分别算峭度再平均那样会丢掉相位信息。正确做法是取复包络 (c_k(n)) 的模的平方作为能量序列然后计算归一化四阶矩[ K_k \frac{E{|c_k(n)|^4}}{(E{|c_k(n)|^2})^2} - 2 ]这个公式和原始谱峭度公式形式一致只不过 (X(t,f)) 被替换成了固定频带内的复包络采样。每个节点算出一个K值整棵树的所有K值就构成了一个二维矩阵一维是中心频率另一维是分解层级带宽。把这个矩阵按照颜色映射画出来就是所谓的Kurtogram图像。我刚开始接触这个算法的时候一直有一个困惑为什么计算每个节点的峭度时不需要做很多次STFT后来想明白了STFT本质上就是一组带通滤波器的输出而这里的树状滤波器组就是在用另一种方式实现多分辨率STFT。每一条从树根到叶子节点的路径对应一个中心频率和一个等效带宽把所有叶子节点的峭度收集起来等价于在一个非均匀网格上采样的谱峭度。所以它的时频分辨率和STFT的时频分辨率遵循同样的不确定性原理想获得更精细的频率分辨率就得用更长的等效时间窗也就意味着每个节点用于统计的独立样本数更少峭度估计的不确定性更大。3.3 从算法输出到工程决策读图要找的不是最亮是最稳Kurtogram图通常用颜色深浅表示峭度值大小新手最容易犯的错误是直接取全局最大值对应的节点参数作为带通滤波器的设置。这个做法在仿真数据上往往没问题但在实测振动数据上经常翻车。原因是实测信号里通常混着随机冲击比如外部敲击、对中不良引起的瞬时冲击这些冲击会在某个窄带产生一个孤立的极高峭度点但它并不是我们关心的周期性故障特征带宽也太窄无法用于包络解调。我的习惯是先看Kurtogram中高峭度区域是否沿某个频率方向和带宽方向连续成片。一个真实可靠的故障共振频带在图像上通常表现为一块颜色明显亮于背景、但内部有一定均匀性的亮斑而不是一个孤立的亮点。换句话说先看连续区域再从连续区域里挑最强点这个原则帮我避开了很多误判。3.4 快速谱峭度的一种工程化快速实现基于FFT的近似除了原始的滤波器组实现还有一种工程中非常实用的近似方法直接用不同长度的FFT窗计算STFT再把STFT结果按频点分组对每个频点的能量序列直接算峭度。这样做不需要构建复杂的滤波器组代码极其简洁在Python里用scipy.signal.stft几行就能搞定。这种方法的精度略低于滤波器组法因为STFT是固定带宽的无法精细覆盖1/3-二叉树的所有节点但用于初步确定共振频带范围完全够用。我在现场的便携式采集器上就是这么干的——先用FFT近似法快速锁定频带再用高精度滤波器组法细化参数两者结合既能保证响应速度又能保证最终参数的可靠性。4. 分解层数K、窗口长度与峭度估计精度最容易忽视的参数组合4.1 分解层数K到底在控制什么Kurtogram的分解层数K直接决定了频率分辨率的下限。K越大最高层子带的带宽越窄理论上能找到更窄的共振峰。但这里有一个容易被忽略的局部样本量问题。假设信号采样频率 (F_s 25600 \text{ Hz})信号时长1秒总共25600个采样点。当K6时最细子带带宽约为 (F_s/(3 \cdot 2^5) 266.7 \text{ Hz})经过两级抽取之后该子带的有效采样点数只剩约25600 / 64 400个。用400个样本估计峭度的四阶矩置信区间已经比较宽了如果K再大到8每个子带只剩大约100个样本峭度估计值的波动会大到完全不可信。所以K值的选择不是越大越好它应该在频带精细度和统计稳定性之间取折中。我的经验是信号时长推荐K值说明0.1 ~ 0.5秒3 ~ 4样本太短高K完全没意义0.5 ~ 2秒4 ~ 5常规振动巡检场景2 ~ 10秒5 ~ 6适合轴承早期微弱故障10秒以上6 ~ 7可进一步观察窄带共振结构4.2 峭度估计的方差问题为什么数据越短越不靠谱峭度估计对异常值极其敏感这是它的固有属性。四阶矩的样本估计量方差很大尤其是当信号不是严格平稳时。实测的振动信号里经常有转速波动、负载变化、瞬态冲击这些都会显著增大峭度估计的不确定度。为了缓解这个问题工程上有几个实用做法分段平均把长时间信号切成多段分别计算每段的谱峭度或节点峭度然后取中位数而不是均值。中位数对偶发强冲击的鲁棒性远好于均值。使用稳健峭度用基于分位数或M估计的替代指标比如四分位距归一化的尾部加权指标本质上是削弱极端值对四阶矩的支配作用。结合其他指标验证谱峭度定位到的频带可以再用谱负熵Spectral Negentropy或Gini系数做二次验证。我实测过在低信噪比情况下Gini系数对微弱周期性冲击的敏感度往往优于峭度两者结合比单独使用任意一个都可靠。4.3 一个实测中常见的坑随机冲击造成的伪最优频带这个坑我在处理齿轮箱数据时踩过不止一次。高速齿轮箱里存在啮合冲击它也会在Kurtogram上产生一块明亮的区域。如果直接用这块区域做带通滤波再包络解调得到的包络谱里会出现啮合频率及其谐波而不是齿轮局部故障的特征频率导致诊断结论直接错误。怎么区分啮合冲击频带和局部故障共振频带可以看峭度随层级的变化趋势真实故障共振频带通常是窄带高峭度在K增加到一定程度后仍保持高值而啮合冲击往往是宽带能量集中的结果在劣化层级增大后峭度会迅速下降。更可靠的方法是结合时域包络的周期性检验——对候选频带做带通滤波之后计算包络自相关如果包络里存在明显的周期成分才说明候选频带里有周期性故障特征。这个先选频带、后验周期的流程极大地降低了误诊率。5. 完整实战滚动轴承早期故障特征提取案例5.1 实验数据与工况说明我用一个公开的轴承故障模拟台数据来演示完整流程。轴的转频约为 30 Hz采样频率 (F_s 25600 \text{ Hz})传感器安装在轴承座上采集垂直方向加速度信号时长2秒。轴承外圈存在一处人工点蚀外圈故障特征频率按理论公式计算约为 91.5 Hz。信噪比很低在频谱图上几乎看不到外圈故障的特征频率边带。这个案例的难点在于故障特征非常微弱直接FFT根本看不清必须依赖选频带包络解调来提升信噪比。5.2 算法执行流程与代码骨架Python示例整个处理流程可以拆成四步预处理、快速谱峭度计算、最优频带确定、包络谱分析。我给出一个可以直接运行的最小实现骨架它用的是STFT近似法便于理解算法逻辑实际工业级应用中可以直接调用成熟的Kurtogram库参数调整思路完全一致。import numpy as np from scipy.signal import stft, hilbert from scipy.fft import fft, fftfreq def spectral_kurtosis_stft(x, fs, windowhann, nperseg_range(64, 1024)): 基于多窗长STFT的谱峭度近似计算。 返回freqs数组、每个窗长下各频点的峭度二维矩阵 kurt_mat [] freqs_out None for nperseg in nperseg_range: f, t, Zxx stft(x, fs, windowwindow, npersegnperseg, noverlapnperseg//2) energy np.abs(Zxx) ** 2 # 每个时频点的能量 # 沿时间方向计算归一化四阶矩再减去高斯基准2 kurt np.mean(energy**2, axis1) / (np.mean(energy, axis1)**2 1e-12) - 2.0 kurt_mat.append(kurt) if freqs_out is None: freqs_out f K np.vstack(kurt_mat) # 行对应不同的带宽层级列对应频率 return freqs_out, K # 读取信号示例x为加速度信号fs25600 # freqs, K spectral_kurtosis_stft(x, fs) # best_idx np.unravel_index(np.argmax(K), K.shape) # 对应最优带宽为 2 * fs / (2 * nperseg_list[best_idx[0]])中心频率为 freqs[best_idx[1]] def bandpass_envelope_demod(x, fs, fc, bw, order4): 带通滤波 包络解调 包络谱 from scipy.signal import butter, sosfilt low max(1e-6, (fc - bw/2) / (fs/2)) high min(1.0, (fc bw/2) / (fs/2)) sos butter(order, [low, high], btypebandpass, outputsos) filtered sosfilt(sos, x) env np.abs(hilbert(filtered)) spec np.abs(fft(env)) freqs_spec fftfreq(len(env), 1/fs) return freqs_spec[:len(env)//2], spec[:len(env)//2]这段代码里有两个细节值得说明为什么用多个窗长因为不同窗长对应不同频率分辨率这模拟了1/3-二叉树的带宽维度。短窗长适合检测较宽的共振带长窗长适合检测较窄的共振带。为什么峭度计算中分母要加1e-12防止某个频点在整段时间内能量恒为零时出现除零错误这在低频段或强滤波情况下很常见。5.3 结果解读如何在Kurtogram上找到最优参数并完成包络谱分析在上面这段2秒数据上运行后最优节点落在中心频率约 7800 Hz、带宽约 800 Hz 的区域对应峭度值约为 24。相比周围区域的峭度值普遍在0~6之间这个频带的冲击特征非常突出。接下来做带通滤波中心频率7800 Hz、带宽800 Hz滤波后提取包络再对包络做FFT。包络谱上能看到很明显的 91.5 Hz 谱线以及其2倍频183 Hz、3倍频274.5 Hz这个边带结构正好对应外圈故障特征频率 (BPFO) 及其谐波诊断结论清晰明确。如果直接对原始信号做FFT不经过上述最优频带选择这段数据的频谱在91.5 Hz附近完全看不到任何谱峰因为故障能量被分散在高频共振带中。这就是为什么快速谱峭度在轴承早期故障诊断里如此重要——它把淹没在噪声中的特征能量先集中再解调本质上是做了一次自适应的信噪比增强。5.4 与其他频带选择方法的横向对比装过各种智能诊断工具后我越来越深刻地认识到没有万能算法关键是要理解每种方法的适用边界。把快速谱峭度和几个典型竞品放在一起对比会更清楚它的定位方法核心原理优势劣势适用场景快速谱峭度Fast Kurtogram1/3-二叉树划分频带搜索峭度极大值速度快自动选择频带易于工程部署对随机冲击敏感峭度估计方差大周期性故障特征明显的轴承、齿轮早期诊断Protrugram在窄带包络谱上计算峭度对随机冲击鲁棒性更好带宽分辨率更细需要预先指定中心频率无法全自动已知故障大致频段后的参数细化Autogram基于包络信号的自相关峭度抑制随机冲击能力强计算量较大调参维度多强背景噪声、存在随机冲击的工况Gini系数法对包络谱计算Gini稀疏度对弱周期性冲击敏感对宽频带故障特征定位不精确微弱特征检测与多方法交叉验证在实际项目中我的常用策略是以快速谱峭度为第一轮粗选锁定候选频带再用Protrugram或Autogram对候选频带做二次精选和验证。这套组合兼顾了自动化和鲁棒性也是目前工业软件里最常见的诊断流程设计。6. 实践中的边界条件与经验总结6.1 几个值得记住的边界条件信号长度不足时结果不可信。K6的情况下建议信号长度至少 ( 2^6 \cdot \text{窗长} ) 个点。现场采样时长如果只有0.2秒就别指望它能给你一个准确的窄带定位老老实实降低K值或换用Protrugram。转速波动环境下谱峭度会失真。因为频率成分随时间漂移峭度的统计特性被破坏。这种情况下一般需要先做阶比跟踪把信号重采样到角域再做谱峭度分析。多故障共存时Kurtogram只能保证找到峭度最大的那个故障不能一次帮你找出所有故障。如果怀疑同一轴承上同时有内圈、外圈故障更稳妥的做法是找到最优频带后把该频带也用于其他特征频率的检索或分两次用不同的频带限制条件排查。滤波器阶数对结果有影响。无论是Kurtogram内部的滤波器组还是后续带通滤波滤波器的过渡带都会影响峭度估计值。建议滤波器选用高阶或高滚降率以避免相邻频带能量泄漏造成峭度上升。最后再分享一个我在实际现场调试时积累的小技巧当你在Kurtogram上看到多个高亮区域而难以取舍时别急着挑最大值。把这些候选频带各自做一次带通滤波和包络解调对比包络谱中信噪比最高的那个而不是峭度最高的那个。因为最终判决工具是包络谱所有前置算法做的都是提升信噪比的辅助工作——而信噪比永远是现场故障诊断的第一准则。判断算法好坏的标准不是计算谁更花哨而是它能不能让故障特征在最终的谱图上一目了然。这个思路不但适用于快速谱峭度也适用于我接触过的任何一种自适应频带选择方法。