当前位置：首页 > news >正文

经典谱估计实战：从BT法到Welch法的演进与权衡

news 2026/5/13 13:16:21

1. 经典谱估计的工程困局：为什么我们需要改进周期图法？

记得第一次用周期图法分析电机振动信号时，我盯着屏幕上锯齿状的频谱曲线愣住了——理论上50Hz的工频分量，在实际谱图上却像心电图一样剧烈波动。这种"毛刺现象"正是经典谱估计最典型的痛点：方差大得让工程师怀疑人生。

传统周期图法直接对信号做傅里叶变换求功率谱，数学表达很简单：

import numpy as np def periodogram(x): N = len(x) fft_result = np.fft.fft(x) return np.abs(fft_result)**2 / N

但实际测试时会发现，同一设备在相同工况下，每次采集信号算出的谱图都像不同的指纹。这是因为周期图法本质是功率谱的极大似然估计，当数据长度N有限时，估计方差根本不会随着N增大而减小。我在某风机故障诊断项目中做过统计：10次重复测量中，100Hz分量幅值的波动范围竟达到±30%，这显然无法用于精确诊断。

更麻烦的是频率分辨率与方差不可兼得的悖论。理论上分辨率取决于采样时长T，T越长分辨率越高。但实际增加T会导致：

计算量指数级增长（还记得被1024点FFT支配的恐惧吗？）
长时信号更难满足平稳性假设
方差依然顽固存在

这种困局催生了BT法的出现——先计算自相关函数再做傅里叶变换。虽然数学上等价于周期图法，但计算过程更稳定。不过实测发现，当信噪比较低时，自相关函数的估计误差会被傅里叶变换放大，最终谱图会出现"伪峰"。我曾用BT法分析轴承故障信号时，就误把噪声引起的伪峰当成了故障特征频率，差点引发误停机事故。

2. 从Bartlett到Welch：工程师的渐进式改良之路

2.1 Bartlett平均法：用分辨率换稳定性的第一代方案

第一次成功降服频谱方差，是在某汽车NVH测试项目中。当时尝试了Bartlett的分段平均策略：把10秒振动信号切成20段0.5秒数据，分别计算周期图再平均。效果立竿见影——频谱曲线终于不再"跳舞"了。

Bartlett法的核心公式：

def bartlett_psd(x, L): M = len(x) // L # 每段长度 psd_sum = np.zeros(M) for i in range(L): segment = x[i*M : (i+1)*M] psd_sum += periodogram(segment) return psd_sum / L

但代价也很明显：

频率分辨率从原来的0.1Hz降到了2Hz（Δf=1/T=1/0.5s）
丢失了200Hz以上的高频成分（M=500对应Nyquist频率=250Hz）

这个案例让我深刻理解了估计方差与分辨率的Trade-off。当处理转子启停机信号时，这种缺陷尤为致命——转速变化导致特征频率移动，需要高分辨率追踪，但Bartlett法会模糊这些关键细节。

2.2 Welch法：给工程师的"瑞士军刀"

真正让我工作效率飞跃的，是掌握了Welch法的三项关键改良：

50%重叠分段：把信号像瓦片一样交错分割，使数据利用率提升近一倍。实测显示，在相同方差水平下，Welch法比Bartlett法分辨率提高约30%
汉宁窗应用：用这个余弦窗函数压制频谱泄漏。还记得第一次看到加窗前后的对比：原本被泄漏淹没的齿轮啮合边频带，终于清晰可辨
分段长度动态调整：通过M值灵活平衡实时性与精度。在在线监测系统中，我通常设置M=2048（4.096秒@500Hz采样），L=8~12段

Python实现Welch法的核心参数：

from scipy import signal f, Pxx = signal.welch( x, # 输入信号 fs=500, # 采样率 window='hann', # 汉宁窗 nperseg=1024, # 段长度 noverlap=512, # 重叠点数 nfft=2048, # FFT点数 scaling='density' # 功率谱密度 )

在某风电齿轮箱诊断中，对比三种方法的效果：