怎么用Python实现振动信号的包络谱故障诊断

一、为什么直接 FFT 难以发现轴承早期损伤？

设想这样一种典型工况：一台转速为 1800 rpm（转频 30 Hz）的工业离心风机，其滚动轴承外圈出现了一个微小的点蚀斑点。每当滚动体碾过该缺陷点时，都会产生一次持续时间仅微秒级的机械瞬态撞击。

这种冲击在频域上具有宽带激励特性，会激发出传感器与轴承支撑结构的高频固有共振（通常在几 kHz 甚至十几 kHz 范围内）。因此，传感器拾取到的信号本质上是一个受到周期性冲击调制的衰减高频振荡波。如果直接对这段原始时域波形做快速傅里叶变换（FFT），常常会遇到两个致命困难：

• 低频段工频能量遮蔽：电机的动不平衡、轴对中不良、偏心等宏观几何运动会产生极大的低频工频振动（1X、2X 转频），其能量常常比早期点蚀微弱冲击高出数十甚至数百倍，直接把低频处的轴承故障特征谱线彻底淹没。
• 高频共振带能量分散：虽然冲击激发的共振能量集中在高频带，但它表现为以固有共振频率为载波中心、以故障特征频率为间隔的密集边频带。在强烈背景白噪声干扰下，各条边频微弱且信噪比低，根本无法肉眼辨识或阈值告警。

包络谱分析的破局思路在于“解调”：既然故障信息调制在高频共振载波的外轮廓（即幅值包络）中，我们首先使用带通滤波器扣出共振频段，再通过数学变换剥离掉高频振荡载波，仅保留随时间周期性涨落的“外包络线”。对这条纯净的低频包络做 FFT，原本隐藏在载波内部的故障特征频率及其倍频就会以高幅值谱峰的形式干净利落地显现出来。

二、轴承四类理论故障特征频率的计算基准

在执行 Python 解调之前，必须先根据轴承几何参数与工作转速，计算出构件出现点蚀损伤时的理论撞击频率基准。滚动轴承的核心几何参数包括：

• fr：转轴旋转频率（Hz），若转速为 N rpm，则 fr = N / 60；
• Z：单列滚动体（滚珠或滚子）数量；
• d：滚动体平均直径（mm）；
• D：轴承滚道节径（Pitch Diameter，mm）；
• alpha：接触角（度数，向心推力轴承取非零角度，深沟球轴承纯径向取 0°）。

在 Python 中，可以通过一个轻量函数快速封装上述计算公式：

import numpy as np

def calculate_bearing_frequencies(fr, Z, d, D, alpha_deg=0.0):
    alpha = np.radians(alpha_deg)
    gamma = (d / D) * np.cos(alpha)
    
    bpfo = (Z * fr / 2.0) * (1.0 - gamma)
    bpfi = (Z * fr / 2.0) * (1.0 + gamma)
    bsf  = (D * fr / (2.0 * d)) * (1.0 - gamma**2)
    ftf  = (fr / 2.0) * (1.0 - gamma)
    
    return {"BPFO": bpfo, "BPFI": bpfi, "BSF": bsf, "FTF": ftf}

# 示例：转速 1800 rpm (fr = 30 Hz), 滚珠数 Z = 9, 滚珠直径 d = 7.94 mm, 节径 D = 39.04 mm
freqs = calculate_bearing_frequencies(fr=30.0, Z=9, d=7.94, D=39.04)
# freqs['BPFO'] 约等于 107.5 Hz

三、包络解调的核心数学实现：Hilbert 变换

包络提取的核心技术是利用希尔伯特变换（Hilbert Transform）构建解析信号（Analytic Signal）。

对于一个经过带通滤波提取出的带通冲击信号 x(t)，其希尔伯特变换定义为 x(t) 与 1/(πt) 的卷积积分：

希尔伯特变换在频域上的物理意义非常直观：它将输入信号的正频分量全部顺时针旋转 90°（乘以 -j），将负频分量全部逆时针旋转 90°（乘以 +j），而保持幅频特性不变。将原实信号作为实部，希尔伯特变换后的正交信号作为虚部构造成复数信号 z(t)，就构成了单边谱解析信号。此时对 z(t) 求模值 |z(t)|，便完美抽离出了调制载波上的瞬时包络曲线 A(t)。

四、完整的 Python 实战代码与仿真验证

为了让每一位读者无需外接工业传感器就能直接调试并验证算法，下面的 Python 脚本构造了一个贴合真实工况的轴承外圈点蚀仿真信号：包含旋转工频低频杂散、周期性外圈冲击（BPFO = 105.0 Hz）、结构共振高频响应（中心频率 3000 Hz、阻尼指数衰减）以及高斯背景白噪声。

环境依赖仅需 numpy、scipy 和 matplotlib，在终端使用 pip install numpy scipy matplotlib 即可安装。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# ==========================================
# 1. 模拟生成具有微小外圈点蚀的振动信号
# ==========================================
def generate_synthetic_vibration(fs=20000, duration=1.0, fr=30.0, bpfo=105.0, f_res=3000.0):
    t = np.linspace(0, duration, int(fs * duration), endpoint=False)
    N = len(t)
    
    # (A) 转轴旋转工频干扰 (1X, 2X 转频，能量较大)
    rotational_noise = 1.5 * np.sin(2 * np.pi * fr * t) + 0.8 * np.sin(2 * np.pi * (2 * fr) * t)
    
    # (B) 构造外圈周期性冲击脉冲序列 (周期 T = 1 / BPFO)
    impact_train = np.zeros(N)
    impact_interval = int(fs / bpfo)
    for idx in range(0, N, impact_interval):
        impact_train[idx] = 1.0
        
    # (C) 结构高频共振冲激响应 (3000 Hz 衰减正弦波)
    decay_len = int(fs * 0.008) # 衰减时间常数约 8ms
    t_decay = np.linspace(0, 0.008, decay_len, endpoint=False)
    impulse_response = np.exp(-600 * t_decay) * np.sin(2 * np.pi * f_res * t_decay)
    
    # 冲击卷积共振响应
    fault_signal = signal.convolve(impact_train, impulse_response, mode='same') * 0.6
    
    # (D) 添加随机高斯背景白噪声
    white_noise = np.random.normal(0, 0.5, N)
    
    # 合成最终传感器拾取信号
    raw_signal = rotational_noise + fault_signal + white_noise
    return t, raw_signal

# ==========================================
# 2. 包络谱分析核心诊断流水线
# ==========================================
def compute_envelope_spectrum(x, fs, band_low=2000.0, band_high=4000.0, filter_order=3):
    # 步骤 A：设计 Butterworth 二阶/四阶带通滤波器 (使用 SOS 格式提高数值稳定性)
    sos = signal.butter(filter_order, [band_low, band_high], btype='bandpass', fs=fs, output='sos')
    # 双向零相位失真滤波
    x_filtered = signal.sosfiltfilt(sos, x)
    
    # 步骤 B：执行希尔伯特变换求取解析信号
    analytic_signal = signal.hilbert(x_filtered)
    
    # 步骤 C：提取幅值包络并去除直流偏置
    envelope = np.abs(analytic_signal)
    envelope_ac = envelope - np.mean(envelope)
    
    # 步骤 D：计算包络信号的单边快速傅里叶变换 (FFT)
    N = len(envelope_ac)
    fft_vals = np.fft.rfft(envelope_ac)
    freqs = np.fft.rfftfreq(N, d=1.0 / fs)
    amplitude_spectrum = np.abs(fft_vals) * (2.0 / N)
    
    return x_filtered, envelope_ac, freqs, amplitude_spectrum

# ==========================================
# 3. 运行执行与结果验证
# ==========================================
if __name__ == "__main__":
    fs = 20000         # 采样率 20 kHz
    duration = 1.5     # 采样时长 1.5 秒 (提供 0.67 Hz 的频率分辨率)
    target_bpfo = 105.0# 设定理论外圈故障特征频率 105 Hz
    
    # 生成测试振动信号
    t, raw_vib = generate_synthetic_vibration(fs=fs, duration=duration, fr=30.0, bpfo=target_bpfo, f_res=3000.0)
    
    # 执行包络谱算法 (选取 2000~4000 Hz 覆盖 3000 Hz 的结构共振带)
    x_band, env, f_axis, env_spec = compute_envelope_spectrum(
        raw_vib, fs=fs, band_low=2000.0, band_high=4000.0
    )
    
    # 在低频诊断区 [0 ~ 500 Hz] 查找最高谱峰
    diag_mask = (f_axis >= 10.0) & (f_axis <= 400.0)
    peak_idx = np.argmax(env_spec[diag_mask])
    detected_freq = f_axis[diag_mask][peak_idx]
    
    print("=== 诊断比对结果 ===")
    print(f"理论故障特征频率 BPFO: {target_bpfo:.2f} Hz")
    print(f"包络谱检出最高特征峰: {detected_freq:.2f} Hz")
    print(f"诊断判定: {'外圈局部剥落故障明确' if abs(detected_freq - target_bpfo) < 2.0 else '未见明显特征'}")

示例执行输出说明

运行上述代码后，控制台将输出如下诊断结论：

=== 诊断比对结果 ===
理论故障特征频率 BPFO: 105.00 Hz
包络谱检出最高特征峰: 105.33 Hz
诊断判定: 外圈局部剥落故障明确

检出频率（105.33 Hz）与理论值（105.00 Hz）的细微微差来自于 1.5 秒采样时长下的离散 FFT 频率栅格分辨率（Δf = 1 / 1.5 ≈ 0.67 Hz），诊断精度完全在工程允许的误差范围之内。

五、工程落地避坑指南与参数调优排查

在实验室仿真中，因为参数已知，诊断往往一击即中；但在工程现场实测时，许多工程师却常遇到“包络谱做出来全是一片平坦白噪声”或“谱峰对不上”的问题。结合现场实施经验，梳理如下四维排查决策指南：

六、总结与实践建议

利用 Python 实现振动信号的包络谱故障诊断，本质上是一个清晰的“物理滤波 + 数学解调 + 谱分析”闭环：

1. 定位共振带：借由结构共振的高频载波放大微弱冲击能量，通过带通滤波剥离低频工频机械运动和大能量背景干扰；
2. 提取外包络：使用 scipy.signal.hilbert 构造正交解析信号并取模，捕获由冲击引发的瞬态衰减包络曲线；
3. 必须去直流：减除均值以消除幅值模值的恒正偏置，避免 0 Hz 直流巨峰遮盖；
4. 对齐特征峰：执行 FFT 并与几何推导的 BPFO、BPFI、BSF 及其谐波进行对齐校核，结合滑差修正与边频带形态做出明确判断。

掌握这套简洁严谨的 Python 流水线，不仅能让你在离线分析科研中快速验证各种故障轴承数据集（如西储大学 CWRU 轴承数据、江南大学轴承数据等），也能便捷地集成到工业物联网（IIoT）边缘网关与自动化设备健康监测平台中。