Python信号处理实战：用SciPy和NumPy给振动信号做个‘高阶体检’（双谱图入门）
Python信号处理实战：用SciPy和NumPy给振动信号做个‘高阶体检’（双谱图入门）
轴承振动信号里藏着什么秘密？当传统频谱分析只能告诉你"有什么频率"时，双谱图却能揭示"这些频率如何相互作用"。想象一下，你手里有一段来自工业设备的振动数据，常规分析显示频谱正常，但设备依然异常发热——这时候，就需要请出我们今天的主角：双谱图这个能看穿非线性关系的"X光机"。
1. 从振动信号到双谱图：为什么要做高阶体检？
去年检修一台离心泵时，我遇到个典型案例：FFT频谱显示各轴频幅值都在安全范围内，但设备运行时总伴随不规则异响。直到用双谱图分析，才发现20Hz和35Hz两个看似正常的振动分量，其实存在强烈的二次相位耦合——这种隐藏在常规频谱背后的"暗流"，正是轴承座松动的特征指纹。
双谱图相比传统功率谱有三大不可替代的优势：
- 相位关系可视化：能捕捉频率分量间的相位耦合（比如f1+f2=f3时是否存在固定相位关系）
- 非线性检测：对高斯噪声天然免疫，特别适合检测轴承裂纹、齿轮啮合异常等非线性特征
- 调制分析：清晰显示载荷波动导致的幅值调制现象（常见于不对中故障）
# 模拟轴承外圈故障的典型调制信号 import numpy as np fs = 4096 # 采样频率 t = np.linspace(0, 1, fs) carrier = 120 # 轴承特征频率 modulation = 30 # 转频 x = (1 + 0.5*np.sin(2*np.pi*modulation*t)) * np.sin(2*np.pi*carrier*t)注意：实际工业信号往往还包含随机冲击成分，上述代码仅展示最简单的幅值调制现象
2. 双谱图实验室：亲手搭建分析管道
理解双谱图最好的方式就是动手实现它。我们将使用SciPy的signal模块构建完整分析流程，关键步骤包括：
2.1 信号预处理：给振动数据"降噪"
真实传感器数据总免不了噪声污染。先加载示例数据并做初步处理：
from scipy import signal import matplotlib.pyplot as plt # 生成带噪声的仿真信号（含非线性相互作用） t = np.linspace(0, 1, 8000) f1, f2 = 50, 85 signal_clean = np.sin(2*np.pi*f1*t) * np.sin(2*np.pi*f2*t) noise = 0.2 * np.random.randn(len(t)) x = signal_clean + noise # 带通滤波设计 b, a = signal.butter(4, [40, 90], 'bandpass', fs=8000) x_filtered = signal.filtfilt(b, a, x)参数选择技巧：
- 滤波器阶数不宜过高（避免相位失真）
filtfilt实现零相位滤波- 截止频率应根据设备特征频率设置
2.2 双谱计算核心：从STFT到双谱矩阵
双谱本质上是三阶累积量的二维傅里叶变换。实际操作中常用间接法计算：
# 计算短时傅里叶变换(STFT) f, t, Zxx = signal.stft(x_filtered, fs=8000, nperseg=256, noverlap=128) # 构建双谱矩阵 B = np.zeros((len(f), len(f)), dtype=complex) for i in range(len(f)): for j in range(len(f)): k = np.argmin(np.abs(f - (f[i]+f[j]))) # 寻找f1+f2对应的频点 B[i,j] = Zxx[i] * Zxx[j] * np.conj(Zxx[k])关键参数解析：
| 参数 | 作用 | 典型取值 |
|---|---|---|
nperseg | 分段长度 | 256-1024 |
noverlap | 重叠样本数 | 50%-75%分段长度 |
window | 窗函数类型 | 'hann'最常用 |
2.3 可视化：把矩阵变成诊断热图
双谱图的解读关键在于观察对称性和亮点分布：
plt.figure(figsize=(10,8)) plt.imshow(np.log10(np.abs(B)), extent=[f[0], f[-1], f[0], f[-1]], origin='lower', aspect='auto', cmap='jet') plt.colorbar(label='Magnitude (dB)') plt.plot([f[0], f[-1]], [f[0], f[-1]], 'w--', linewidth=0.5) # 对角线参考线 plt.title('Bispectrum of Bearing Vibration') plt.xlabel('Frequency f1 (Hz)') plt.ylabel('Frequency f2 (Hz)')健康信号的双谱图应呈现：
- 主对角线附近能量集中
- 非对角线区域均匀分布
- 无明显孤立亮点
3. 实战诊断：从双谱图中读出故障密码
某风机轴承振动信号分析案例：
3.1 正常状态基准
健康轴承的双谱图呈现典型特征：
- 能量集中在对角线附近（线性主导）
- 整体分布均匀对称
- 积分双谱值低于阈值（示例中为0.15）
# 计算归一化双谱积分 B_norm = np.abs(B) / np.max(np.abs(B)) bispec_integral = np.sum(B_norm) / (len(f)**2)3.2 故障信号特征提取
当出现外圈故障时，双谱图会显示：
- 非对角线出现明显亮点（如f1=85Hz, f2=120Hz位置）
- 存在对称的相位耦合点
- 积分值显著上升（>0.3需预警）
# 自动检测显著耦合点 from skimage.feature import peak_local_max peaks = peak_local_max(np.log10(np.abs(B)), min_distance=5, threshold_rel=0.2) # 标记可疑点 for peak in peaks: if not (peak[0] == peak[1]): # 排除对角线上的点 plt.scatter(f[peak[1]], f[peak[0]], s=100, edgecolors='w', facecolors='none') print(f'发现非对角耦合点：({f[peak[1]]:.1f}Hz, {f[peak[0]]:.1f}Hz)')4. 进阶技巧：提升双谱分析效果的五个关键
分段策略优化：
- 对于瞬态冲击信号，采用更短的分段（如128点）
- 平稳信号可用较长分段提升频率分辨率
窗函数选择：
# 不同窗函数对比 windows = ['hann', 'blackman', 'flattop'] fig, axes = plt.subplots(1, 3, figsize=(15,4)) for ax, win in zip(axes, windows): f, t, Zxx = signal.stft(x, window=win, nperseg=256) # ...计算并绘制双谱... ax.set_title(f'{win} window')多分辨率分析：
- 先用大nperseg定位可疑频段
- 再对小段数据精细分析
时变信号处理：
# 滑动窗口分析 for i in range(0, len(x)-800, 200): segment = x[i:i+800] # 计算该段双谱 # 监控特征点能量变化与其他特征参数联动：
- 结合包络谱分析
- 对比小波变换结果
- 参考振动烈度趋势