Python仿真轴承故障信号：从时域波形到频谱与包络谱分析实战-平芜编程栈

简介：本资源面向机械故障诊断、信号处理方向的本科生、研究生及工程技术人员，提供轴承典型故障（如内圈、外圈、滚动体缺陷）的仿真信号分析完整实践方案。资源包含时域波形图与频谱图的生成与可视化代码，辅以直观的操作录屏视频，帮助学习者理解故障特征频率提取、FFT变换原理及MATLAB信号分析流程。压缩包共3个文件（2个MATLAB脚本+1个AVI操作录像），总大小347KB，结构精炼：runme.m为主控入口，fft30.m封装核心频谱计算逻辑，录像视频详细演示路径设置、函数调用与图形解读全过程。已有2944人学习下载，特别适合零基础入门振动信号分析、需快速掌握故障特征识别方法的学习者，避免因路径错误或版本兼容问题导致运行失败。

1. 项目概述：从振动信号中“听”出轴承的健康状况

在工业设备运维领域，轴承就像设备的“关节”，它的健康状况直接决定了整台机器的运行寿命和稳定性。想象一下，一台大型风机或一台精密机床，如果其核心轴承出现了早期故障但未被发现，轻则导致停机停产，重则引发连锁反应，造成严重的安全事故和经济损失。因此，轴承故障的早期诊断与预警，是设备预测性维护的基石。

传统的故障诊断依赖老师傅“听音辨位”的经验，或者定期停机拆解检查，不仅效率低下，而且容易错过最佳维修窗口。随着传感器技术和信号处理算法的发展，基于振动信号的分析成为了主流的、非侵入式的诊断方法。本项目标题“轴承故障仿真信号时域波形图和频谱图+含代码操作演示视频”，其核心价值就在于提供一个从理论到实践、从数据到可视化的完整学习路径。它不仅仅是一张图或一段代码，而是一套完整的故障诊断“工具箱”入门指南。

通过生成模拟的轴承故障振动信号，我们可以不受真实实验条件和昂贵设备的限制，快速理解不同类型故障（如内圈、外圈、滚动体损伤）在信号上的特征表现。时域波形图让我们直观“看到”振动冲击的形态和周期，而频谱图（通常指功率谱密度图）则能将这些时域冲击“翻译”成频率成分，揭示故障的特征频率。配合代码和视频，学习者可以亲手复现这一分析过程，将抽象的数学公式和信号处理概念，转化为屏幕上跳动的波形和清晰的频谱峰，从而深刻掌握轴承故障诊断的核心技能。

这篇文章，我将以一个做过不少故障诊断项目的老工程师的身份，带你彻底拆解这个仿真分析流程。我会分享从信号生成原理、代码实现细节，到图形解读技巧和实际应用中的避坑经验，让你不仅能“跑通”代码，更能“读懂”信号，真正具备解决实际问题的能力。

2. 核心思路与方案设计：如何“无中生有”地模拟故障信号

要仿真轴承故障信号，我们不能凭空捏造数据，必须基于其物理机理和数学模型。一个健康的轴承在运行时，其振动信号主要包含转频及其倍频，以及一些随机噪声，整体比较“干净”。而当轴承的滚道或滚动体出现局部损伤（如点蚀、剥落）时，滚动体经过损伤点会产生周期性的冲击力，从而引发高频共振，这个物理过程是我们仿真的核心。

2.1 故障信号的数学模型解析

目前最经典且有效的模型是“周期性冲击响应模型”。我们可以把这个过程想象成用小锤子周期性地敲击一个钟。小锤子敲击的周期由故障特征频率决定，而钟被敲击后产生的“嗡鸣”声，就是系统的固有频率响应。数学模型上，它通常由以下几部分卷积或叠加而成：

周期性冲击序列：模拟故障点被周期性撞击的事件。这是一个间隔为故障周期（T = 1/故障特征频率）的狄拉克δ函数序列。
衰减正弦波（共振响应）：模拟轴承-传感器系统受到冲击后的瞬态振动响应。通常用一个指数衰减的正弦函数来表示：h(t) = exp(-β*t) * sin(2π*f_n*t)。其中，β是衰减系数，决定了共振衰减的快慢；f_n是系统的共振频率，通常在几千赫兹的范围。
背景噪声与谐波：为了更贴近真实信号，还需要加入：
- 高斯白噪声：模拟测量中的随机干扰。
- 转频谐波：模拟由转子不平衡、不对中等引起的低频振动成分。

因此，仿真信号的数学表达式可以概括为：x(t) = [周期性冲击序列] ⊗ [衰减正弦波] + 转频谐波 + 高斯白噪声（其中⊗表示卷积运算）。

为什么选择这个模型？因为它物理意义清晰，能很好地复现真实故障信号的两个核心特征：时域的周期性冲击和频域的共振频带。这对于初学者理解故障机理和后续的特征提取算法（如包络分析）至关重要。

2.2 关键参数的计算与选择依据

仿真不是乱真，参数的选择需要有理有据。这里有几个核心参数需要计算和设定：

轴承几何参数与故障特征频率：这是仿真的起点。你需要知道轴承的节径、滚子数、接触角等。故障特征频率（内圈、外圈、滚动体、保持架）有标准的计算公式。例如，外圈故障特征频率（BPFO）的计算公式为：BPFO = (n * f_r / 2) * (1 - (d/D) * cosα)其中，n是滚子数，f_r是轴旋转频率，d是滚子直径，D是节径，α是接触角。在仿真中，我们通常会直接设定一个特征频率值（如100Hz）来简化，但你必须清楚这个值的来源。
系统共振频率 (f_n)：这通常需要根据实际被测设备的经验值或实验数据来设定。对于中小型轴承，共振频率可能在1kHz到5kHz之间。在仿真中，我们常将其设定在2kHz-3kHz，以在频谱上产生明显且独立于故障频率的共振频带。
采样频率 (fs)：根据奈奎斯特采样定理，fs必须大于信号中最高频率成分的2倍。由于我们的信号包含共振频率f_n，因此fs至少需要大于2 * f_n。通常为了有好的波形显示效果，我们会设置fs为f_n的5-10倍以上，例如fs = 20kHz。
衰减系数 (β)：它控制了共振衰减的速度，直接影响时域冲击波的“宽度”和频域共振峰的“胖瘦”。β值越大，衰减越快，时域冲击越窄，频域共振峰越宽。一般通过调整它来匹配真实信号的衰减特性。

方案选型考量：为什么不直接用真实数据？因为对于初学者或算法开发初期，真实数据获取成本高、工况单一、故障标签难以精确获取。仿真数据成本低、可重复、参数可控，是验证算法流程、理解故障特征的绝佳工具。本方案平衡了模型的复杂度和真实性，足以支撑后续的时频分析和特征提取学习。

3. 核心细节解析与实操要点

理解了模型，我们就要动手用代码把它实现出来。这里我用Python（因其在科学计算和可视化方面的强大生态）来演示核心步骤，并穿插讲解每个环节的注意事项。

3.1 环境准备与库的选择

工欲善其事，必先利其器。我们主要依赖以下几个库：

import numpy as np import matplotlib.pyplot as plt from scipy import signal import warnings warnings.filterwarnings('ignore') # 可选，忽略一些不影响运行的警告

numpy：进行高效的数组运算和数学计算，生成时间序列和信号数据。
matplotlib：绘制时域波形图和频谱图的核心工具，其功能强大且定制灵活。
scipy.signal：提供信号处理相关的函数，如卷积、滤波、谱估计等。

实操心得：版本一致性。在进行项目复现或团队协作时，务必记录或使用requirements.txt文件锁定这些库的版本。不同版本间API的细微变化可能导致代码报错或图形显示异常。

3.2 信号生成函数的关键实现

我们将核心模型封装成一个函数。以下是关键代码段和逐行解析：

def generate_bearing_fault_signal(duration, fs, fault_freq, resonance_freq, decay_beta, snr_db, rpm): """ 生成轴承故障仿真信号。 参数: duration: 信号时长 (秒) fs: 采样频率 (Hz) fault_freq: 故障特征频率 (Hz)，如BPFO resonance_freq: 系统共振频率 (Hz) decay_beta: 共振衰减系数 snr_db: 信噪比 (dB) rpm: 轴转速 (转/分钟) 返回: t: 时间轴 x: 仿真信号 """ # 1. 生成时间轴 t = np.arange(0, duration, 1/fs) N = len(t) # 2. 计算轴旋转频率 shaft_freq = rpm / 60.0 # 单位转为 Hz # 3. 生成周期性冲击序列 (狄拉克梳状函数近似) # 思路：在冲击发生的时刻点，设置一个幅值，其他时间为0 impulse_train = np.zeros(N) # 计算冲击发生的时刻（每隔 1/fault_freq 秒一次） impulse_times = np.arange(0, duration, 1/fault_freq) # 将冲击时刻转换为最接近的采样点索引 impulse_indices = (impulse_times * fs).astype(int) # 确保索引不超出数组范围 impulse_indices = impulse_indices[impulse_indices < N] impulse_train[impulse_indices] = 1.0 # 设置冲击幅值为1 # 4. 生成衰减正弦波（单自由度系统冲击响应） # 时间轴从0开始，足够长以覆盖衰减过程 t_response = np.arange(0, 0.1, 1/fs) # 0.1秒的响应时间通常足够 h = np.exp(-decay_beta * t_response) * np.sin(2 * np.pi * resonance_freq * t_response) # 5. 卷积：冲击序列通过系统响应 fault_impact = np.convolve(impulse_train, h, mode='same') # 归一化，防止幅值过大 fault_impact = fault_impact / np.max(np.abs(fault_impact)) if np.max(np.abs(fault_impact)) > 0 else fault_impact # 6. 添加转频谐波（模拟低频振动） harmonic_signal = 0 for i in range(1, 4): # 添加前3阶谐波 harmonic_signal += 0.1 * (1/i) * np.sin(2 * np.pi * i * shaft_freq * t) # 7. 合成确定性信号 deterministic_signal = fault_impact + harmonic_signal # 8. 添加高斯白噪声 signal_power = np.mean(deterministic_signal ** 2) noise_power = signal_power / (10 ** (snr_db / 10)) noise = np.random.randn(N) * np.sqrt(noise_power) # 9. 最终信号 x = deterministic_signal + noise return t, x

关键点解析与避坑指南：

冲击序列的生成：我们使用单位脉冲来近似狄拉克函数。这里一个常见的坑是时间对齐。impulse_times * fs得到的是理论采样点位置，可能是小数，必须转换为整数索引。astype(int)是向下取整，可能引入微小的时间误差，但对于仿真演示可以接受。更精确的做法是用插值，但复杂度会增加。
卷积模式的选择：np.convolve的mode='same'保证了输出长度与输入impulse_train相同，方便后续处理。但卷积会使信号在时间上产生一定的延迟和展宽，这是物理过程的真实反映。
归一化处理：卷积后信号的幅值可能很大，进行归一化（/ np.max(np.abs(...))）是为了方便控制最终信号的幅值范围，并与后续添加的谐波、噪声在量级上匹配。注意判断分母为零的情况。
信噪比的控制：根据定义，信噪比SNR (dB) = 10 * log10(信号功率 / 噪声功率)。因此，噪声功率 = 信号功率 / (10^(SNR/10))。高斯噪声的标准差就是噪声功率的平方根。这是仿真中控制噪声水平的精确方法。
参数敏感度：decay_beta和resonance_freq是需要反复调试的参数。beta太小，冲击衰减慢，时域波形会混叠；beta太大，冲击不明显。resonance_freq决定了频谱图中共振频带的位置，需要与采样频率协调，避免混叠。

4. 实操过程与核心环节实现

有了信号生成函数，我们就可以生成数据并绘制核心的时域波形图和频谱图了。这部分是可视化的关键，直接决定了我们能否从图形中解读出故障信息。

4.1 参数设置与信号生成实例

我们设定一组“典型”参数来生成一个外圈故障信号：

# 仿真参数设置 duration = 1.0 # 信号时长 1秒 fs = 20000 # 采样频率 20kHz fault_freq = 100 # 故障特征频率 (BPFO) 100Hz resonance_freq = 3000 # 系统共振频率 3kHz decay_beta = 500 # 衰减系数 snr_db = 10 # 信噪比 10dB rpm = 1800 # 转速 1800转/分钟 # 生成信号 t, x = generate_bearing_fault_signal(duration, fs, fault_freq, resonance_freq, decay_beta, snr_db, rpm)

4.2 时域波形图的绘制与解读

时域图是我们观察信号的第一个窗口。

# 绘制时域波形图 plt.figure(figsize=(12, 4)) plt.plot(t, x, 'b-', linewidth=0.5, alpha=0.7) plt.xlabel('时间 [秒]') plt.ylabel('幅值') plt.title('轴承故障仿真信号 - 时域波形图') plt.grid(True, linestyle='--', alpha=0.5) plt.xlim([0, 0.2]) # 只看前0.2秒，细节更清晰 plt.tight_layout() plt.show()

图形解读与技巧：

观察周期性冲击：在时域波形中，你应该能看到一系列间隔大致相等的“冲击包络”。冲击之间的时间间隔就是故障周期T = 1/fault_freq = 0.01秒。在前0.2秒内，大约能看到20个冲击。
冲击的形态：每个冲击都不是一个单纯的尖脉冲，而是一个经过振荡衰减的“小波包”，这正是衰减正弦波响应卷积的结果。振荡的频率就是共振频率f_n。
调整显示范围：使用plt.xlim([0, 0.2])局部放大，是观察细节的关键。全时段1秒的图可能因为数据点太多而显得杂乱，无法看清周期结构。
图形美化：linewidth=0.5和alpha=0.7让曲线更细且半透明，在数据密集时避免糊成一团。网格线(grid)有助于更准确地读取时间点和幅值。

4.3 频谱图（功率谱密度）的绘制与解读

时域看周期，频域看成分。频谱图能将信号分解到频率域，揭示其内在结构。

# 计算并绘制频谱图 (使用Welch方法估计功率谱密度) plt.figure(figsize=(12, 4)) frequencies, psd = signal.welch(x, fs, nperseg=1024, noverlap=512, scaling='density') plt.semilogy(frequencies, psd, 'r-', linewidth=1) # 纵坐标用对数坐标 plt.xlabel('频率 [Hz]') plt.ylabel('功率谱密度 [V**2/Hz]') plt.title('轴承故障仿真信号 - 频谱图 (功率谱密度)') plt.grid(True, which='both', linestyle='--', alpha=0.5) # 主次网格 plt.xlim([0, fs/2]) # 显示奈奎斯特频率以下的部分 # 标记关键频率 plt.axvline(x=fault_freq, color='g', linestyle='--', alpha=0.7, label=f'故障频率 ({fault_freq}Hz)') plt.axvline(x=resonance_freq, color='orange', linestyle='--', alpha=0.7, label=f'共振频率 ({resonance_freq}Hz)') plt.legend() plt.tight_layout() plt.show()

图形解读与技巧：

谱估计方法的选择：直接对信号做FFT得到的周期图方差大、噪声高。scipy.signal.welch方法采用平均周期图法，将数据分段、加窗、求谱再平均，能获得更平滑、更稳定的频谱估计，是工程上的标准做法。
参数nperseg和noverlap：nperseg是每段的长度（点数），它决定了频率分辨率df = fs / nperseg。nperseg越大，分辨率越高，但方差也越大。noverlap是段之间的重叠点数，通常设为nperseg/2，可以在不减少段数的情况下增加平均次数，平滑频谱。这是一个需要权衡的参数。对于故障诊断，我们通常更关心频率位置而非绝对的谱线精细度，nperseg=1024是一个不错的起点。
对数坐标：plt.semilogy将纵坐标（功率）设置为对数坐标。这是因为信号成分（共振峰、故障频率边带）的功率可能相差好几个数量级，线性坐标下小成分会被淹没。对数坐标能同时清晰展示强成分和弱成分。
关键频率标记：使用plt.axvline在频谱图上用竖线标记出故障特征频率和共振频率，是解读频谱的“导航标”。在生成的频谱图中，你应该看到：
1. 在100Hz (fault_freq)处，可能有一个不太明显的小峰。这是因为周期性冲击的频谱在故障频率及其倍频处有谱线，但能量可能不高。
2. 在3000Hz (resonance_freq)附近，有一个非常突出的共振峰。这是故障诊断中更关键的标志！因为冲击激发了系统共振，大部分能量集中在这里。
3. 仔细观察共振峰周围，有时能看到以故障频率为间隔的边带（sidebands），这是调制现象的表现，是故障的又一强有力证据。在我们的仿真模型中，由于冲击是等幅的，边带可能不明显。更复杂的模型会引入幅值调制来模拟这一现象。

4.4 进阶分析：包络谱演示

为了更清晰地提取故障特征频率，在实际诊断中，我们常对共振频带进行解调分析，即包络分析。其步骤是：带通滤波（围绕共振频率）→ 希尔伯特变换求包络 → 对包络信号做频谱分析。这个频谱就叫包络谱，它能极大地增强故障频率成分。

# 简单的包络分析演示 (以共振频率为中心进行带通滤波) bp_low = resonance_freq - 500 # 带通下限 bp_high = resonance_freq + 500 # 带通上限 sos = signal.butter(N=4, Wn=[bp_low, bp_high], btype='bandpass', fs=fs, output='sos') x_filtered = signal.sosfilt(sos, x) # 滤波 # 希尔伯特变换求包络 analytic_signal = signal.hilbert(x_filtered) amplitude_envelope = np.abs(analytic_signal) # 包络线 # 计算包络信号的频谱 freq_env, psd_env = signal.welch(amplitude_envelope, fs, nperseg=1024, noverlap=512, scaling='density') # 绘制包络谱 plt.figure(figsize=(12, 4)) plt.plot(freq_env, psd_env, 'b-', linewidth=1) plt.xlabel('频率 [Hz]') plt.ylabel('功率谱密度') plt.title('包络谱 (重点关注低频段)') plt.grid(True, linestyle='--', alpha=0.5) plt.xlim([0, 500]) # 包络谱主要看低频部分，聚焦故障频率及其倍频 # 标记故障频率及其倍频 for i in range(1, 6): plt.axvline(x=i*fault_freq, color='r', linestyle=':', alpha=0.5, linewidth=0.8) plt.text(i*fault_freq, plt.ylim()[1]*0.9, f'{i}×', fontsize=8, color='r', ha='center') plt.tight_layout() plt.show()

在绘制的包络谱中，将X轴范围限制在0-500Hz。你会清晰地看到在100Hz、200Hz、300Hz等处出现明显的谱峰，它们分别对应故障特征频率的1倍频、2倍频、3倍频。这是诊断轴承局部故障最直接的频域证据。相比于直接观察原始频谱，包络谱能有效抑制无关的低频转频成分和高频噪声，将故障特征频率成分凸显出来。

5. 常见问题与排查技巧实录

在实际操作和代码调试过程中，你一定会遇到各种问题。下面是我总结的一些典型问题及其解决方法。

5.1 时域图看不到明显的周期性冲击

可能原因1：信噪比(snr_db)设置过低。
- 排查：检查生成的信号x，计算其标准差。尝试将snr_db调高（如20或30），重新生成信号并绘图。
- 解决：在仿真学习阶段，可以先用高信噪比观察理想波形，再逐步降低信噪比以模拟更真实的恶劣环境。
可能原因2：衰减系数(decay_beta)过大。
- 排查：decay_beta过大导致单个冲击衰减过快，能量集中在极短时间内，在图形上可能表现为难以察觉的细刺。
- 解决：减小decay_beta（如从500降到100），让冲击振荡衰减过程更长，在时域上更“宽”，更容易被观察到。
可能原因3：显示时间范围太宽。
- 排查：你是否在绘制全时长（如1秒）的图形？冲击间隔是0.01秒，在1000个数据点中只占10个点，在全局图中就像一粒沙子。
- 解决：使用plt.xlim([0, 0.1])或更小的范围进行局部放大，这是必须掌握的看图技巧。

5.2 频谱图中故障频率处没有峰值

可能原因1：故障频率能量太低，被噪声和谐波淹没。
- 现象：这是正常且符合物理实际的。在早期故障或高噪声环境下，故障特征频率的谱线往往很微弱。
- 解决：不要只盯着原始频谱！这正是引入包络分析的原因。按照4.4节的步骤计算包络谱，故障频率及其倍频通常会变得非常明显。
可能原因2：频率分辨率不足。
- 排查：检查signal.welch中的nperseg参数。频率分辨率df = fs / nperseg。如果df大于故障频率的间隔，谱线可能无法落在精确的频率点上，导致能量泄露，峰值不明显。
- 解决：适当增加nperseg（如2048或4096），但注意这会降低谱估计的平滑度。更好的方法是计算包络谱，它对分辨率的要求相对较低。
可能原因3：冲击序列的生成存在误差。
- 排查：检查impulse_indices的生成。打印前几个冲击索引，看它们是否均匀间隔fs / fault_freq个点。
- 解决：确保fault_freq能被fs整除，或者使用更精确的插值方法生成冲击序列，但这对演示影响通常不大。

5.3 共振峰在频谱上不明显或位置不对

可能原因1：共振频率(resonance_freq)设置过高，接近或超过奈奎斯特频率(fs/2)。
- 排查：确认resonance_freq < fs / 2。例如，fs=10000Hz，resonance_freq就不能超过5000Hz，否则会发生混叠。
- 解决：降低resonance_freq或提高fs。
可能原因2：衰减系数(decay_beta)过大或过小。
- 排查：decay_beta影响共振峰的宽度和高度。beta过大，共振响应持续时间短，频域上共振峰宽而矮；beta过小，共振峰窄而高，但在时域可能与其他冲击混叠。
- 解决：这是一个需要根据“感觉”调整的参数。通常beta在resonance_freq * 0.1到resonance_freq * 0.3之间尝试。例如resonance_freq=3000，beta可以尝试300到900。

5.4 代码运行报错或图形异常

ValueError: operands could not be broadcast together...
- 排查：检查时间轴t和信号x的长度是否一致。在卷积或滤波操作后，数组长度可能发生变化。
- 解决：确保在绘图时plt.plot(t, x)中的t和x维度相同。对于mode='same'的卷积，长度不变。对于滤波，输入输出长度通常一致。
图形显示空白或只有坐标轴
- 排查：检查绘图代码中plt.plot()的参数是否正确，数据x是否全为0或NaN。
- 解决：在生成信号后，打印print(x.min(), x.max(), x.mean())查看数据范围。确保plt.show()被正确调用（在脚本中需要，在Jupyter Notebook中可能自动显示）。
频谱图纵坐标值极小（如1e-30）
- 排查：检查信号x的幅值是否过小（例如未经归一化）。
- 解决：在生成信号时进行适当的归一化，或检查signal.welch的scaling参数是否正确设置为'density'或'spectrum'。

5.5 从仿真到实战的思维转换

最后，也是最关键的一点：仿真和现实的差距。仿真信号是纯净、可控的，而真实信号要复杂得多。

真实信号的复杂性：真实环境中存在大量的背景噪声、其他部件的振动干扰、转速波动、传递路径变化等，故障特征可能非常微弱且被淹没。
仿真参数如何对应现实：仿真中的resonance_freq和decay_beta需要通过对健康设备或已知故障设备的信号进行“系统辨识”来大致估计。这不是一个精确值，而是一个范围。
下一步该学什么：当你通过本项目掌握了故障信号的基本形态和时频分析方法后，下一步应该学习：
1. 处理真实数据：学习如何从振动传感器（如加速度计）采集数据，了解数据采集卡、采样定理在实际中的应用。
2. 信号预处理：去趋势、降噪（如小波降噪）、滤波等技术的实际应用。
3. 特征提取的进阶方法：除了频谱和包络谱，还有小波变换、经验模态分解(EMD)、峭度、峰值因子等时域、频域、时频域特征。
4. 智能诊断算法：如何利用机器学习（如SVM、随机森林）或深度学习（CNN、LSTM）对提取的特征进行自动分类和故障识别。