【实践任务】基于VMD的滚动轴承故障诊断

实践任务基于VMD的滚动轴承故障诊断

• 研究背景与意义

  1. 信号分解技术的发展历程

信号处理是工程科学中的核心技术领域。传统的傅里叶分析技术以正弦函数为基函数，适用于线性和稳态信号的分析。然而，傅里叶变换在分析时变非线性信号时存在无法表述信号时频局部特性的局限性。为处理非平稳信号，研究者相继发展了短时傅里叶变换（STFT）、Wigner-Ville分布、小波分析等方法。这些方法虽然在不同程度上改进了傅里叶分析的性能，但本质上仍属于全局分析方法，其分析性能取决于基函数的预先选取，存在一定的适应性问题。

1. 1. 自适应分解方法

1998年，美国工程院院士黄锷（Norden E. Huang）等人创造性地提出了经验模态分解方法，被认为是2000年以来以傅里叶变换为基础的线性和稳态频谱分析的一个重大突破。EMD方法的革命性在于：它是数据驱动的、自适应的，无需预先设定任何基函数，而是依据数据自身的时间尺度特征来进行信号分解。这一点与建立在先验性谐波基函数上的傅里叶分解和建立在先验性小波基函数上的小波分解方法具有本质性的差别。

2014年，Dragomiretskiy和Zosso提出了变分模态分解，这是一种全新的非递归信号分解方法。与EMD的递归“筛选”模式不同，VMD将信号分解转化为变分优化问题，具有坚实的数学理论基础，其实质是多个自适应维纳滤波组，表现出更好的噪声鲁棒性。

1. 1. 研究意义与应用价值

EMD和VMD作为两种代表性的自适应信号分解方法，在以下领域具有重要应用价值：

1. 机械故障诊断：通过对振动信号进行分解，提取故障特征频率，实现轴承、齿轮等关键零部件的状态监测与故障识别。
2. 生物医学信号处理：用于心电图（ECG）、脑电图（EEG）等生理信号的分析与疾病识别。研究表明，VMD在心电图分析中分类准确率可达98.30%。
3. 电力系统监测：应用于电力负荷预测、电力电缆故障定位等场景，具有较高的分析精度。
4. 地震资料处理：用于地震信号的去噪和烃类检测等地球物理勘探任务。

因此，深入理解EMD和VMD的基本原理，掌握其算法实现，并探索其在工程实际中的具体应用，具有重要的理论意义和实用价值。

• 技术原理与方法

2.1 经验模态分解（EMD）

2.1.1 本征模态函数（IMF）

EMD方法的核心是将复杂信号分解为有限个本征模态函数（Intrinsic Mode Function, IMF）和一个残余量。IMF必须满足以下两个条件：

1. 在整个数据序列中，极值点（局部极大值和局部极小值）的数量与过零点的数量相等或最多相差一个；
2. 在任何时间点，由局部极大值定义的上包络线和由局部极小值定义的下包络线的局部均值为零，即上下包络线关于时间轴局部对称。

2.1.2 EMD筛选过程

EMD通过称为筛选（Sifting）的迭代过程完成信号分解，具体步骤如下：

1. 找出信号x(t)的所有局部极大值点和局部极小值点；
2. 利用三次样条插值分别对极大值点和极小值点进行拟合，构造上包络线和下包络线；
3. 计算上下包络线的均值m(t)；
4. 从原信号中减去包络均值，得到；
5. 检查是否满足IMF的两个条件。若满足，则为第一个IMF分量；否则，将作为新的信号重复步骤1-4，直至满足IMF条件或达到预设的终止准则；
6. 从原信号中减去已提取的IMF分量，得到残余信号；
7. 将作为新的待分解信号，重复步骤1-6，依次分解出第二个、第三个等后续IMF分量。当残余信号为单调函数或变化足够小（极值点数少于预定阈值）时停止分解。

最终，原信号可表示为所有IMF分量与最终残余量之和。

这一分解过程表明EMD具有完备性，信号经分解后能够通过所有IMF分量及剩余分量被精确重构出来。

2.1.3 EMD的主要问题

EMD在实际应用中存在以下不足：

模态混叠：当信号中存在由异常事件（如间断信号、脉冲干扰和噪声等）引起的间歇现象时，同一个IMF中会包含不同尺度的信号成分，或者同一尺度的信号出现在不同的IMF中，导致IMF失去物理意义；

端点效应：在三次样条插值拟合包络线时，由于信号两端缺乏极值点约束，导致拟合误差在迭代中被逐步放大并向内部传播。

缺乏完整理论基础：算法缺少完整的数学框架支撑，在实际计算与应用中还存在筛分迭代停止标准等不确定性问题。

为克服上述问题，研究者提出了集合经验模态分解（EEMD）、完全集合经验模态分解（CEEMD）、自适应噪声完备集合经验模态分解（CEEMDAN）等多种改进方法。

2.2 变分模态分解（VMD）

2.2.1 基本思想

VMD是由Dragomiretskiy和Zosso于2014年提出的一种全新的信号分解方法。与EMD的递归筛选模式根本不同，VMD将信号分解问题转化为一个约束变分问题的求解过程。VMD建立在维纳滤波理论、希尔伯特变换理论和频率混合理论三大基础之上，假设每个“模态”是具有不同中心频率的有限带宽信号，其目标是使各模态的估计带宽之和最小。

2.2.2 变分模型的构造

VMD的分解过程可分为变分问题的构造和求解。模态定义：VMD将每个模态定义为一个调幅-调频信号。

变分问题的构建步骤如下：

1. 对每个模态进行希尔伯特变换，获得其解析信号的单边频谱；
2. 将每个模态的频谱通过乘以指数项解调到相应的基频带；
3. 通过解调信号的高斯平滑度来估计各模态的带宽。

2.2.3 变分模型的求解

为求解上述约束变分问题，引入二次惩罚因子和拉格朗日乘子，将约束问题转化为无约束问题。

采用交替方向乘子法（ADMM）求解该扩展拉格朗日表达式的鞍点，通过交替更新来逐步逼近最优解。

2.2.4 VMD的优点

VMD相比EMD具有以下显著优势：

1. 坚实的数学理论基础：VMD建立在变分优化和维纳滤波的严格数学框架之上，从根本上克服了EMD缺乏理论支撑的缺点；
2. 非递归分解模式：VMD采用变分优化而非递归筛选进行分解，有效避免了EMD的模态混叠问题；
3. 更好的噪声鲁棒性：其实质是多个自适应维纳滤波组，对噪声具有天然的抑制能力；
4. 更高的信号重构精度：研究表明，VMD在均方根误差（RMSE）、误差分析和信号重构等方面均优于EMD。

2.3 EMD与VMD的对比分析

在轴承故障诊断的实际对比中，VMD能够有效消除指数衰减直流偏移，而EMD在此方面存在困难。研究表明，VMD具有高精度和快速收敛的特性，在轴承故障诊断中表现更为优越，所以本次实践采用VMD方法进行仿真实验。

• 基于变分模态分解（VMD）的滚动轴承故障诊断

3.1 实验目的

（1）理解变分模态分解（VMD）在信号分解中的基本原理与优势；

（2）掌握利用VMD实现滚动轴承故障特征提取的方法；

（3）通过仿真故障信号，验证VMD结合包络解调在轴承故障诊断中的有效性；

（4）熟悉Python环境下信号处理与故障诊断的基本流程。

3.2 VMD原理简述

变分模态分解（Variational Mode Decomposition, VMD）是一种完全非递归的信号分解方法。它将信号分解问题转化为约束变分优化问题，假设每个模态都是具有不同中心频率的有限带宽信号，通过交替方向乘子法（ADMM）迭代求解各模态及其中心频率，从而实现信号的自适应分离。

对于输入信号，VMD的目标是分解出K个本征模态函数，使各模态的估计带宽之和最小，同时满足所有模态之和等于原信号。通过引入二次惩罚因子和拉格朗日乘子，将约束问题转化为无约束问题进行求解。

VMD具有良好的噪声鲁棒性和抗模态混叠能力，特别适合处理轴承振动信号这类包含冲击成分和强背景噪声的非平稳信号。

3.3 实验数据与代码

3.3.1 滚动轴承故障特征频率

滚动轴承常见故障包括内圈故障、外圈故障和滚动体故障。以某型号轴承为例，其理论特征频率公式如下：

内圈故障频率（BPFI）：

外圈故障频率（BPFO）：

本次实验仿真参数：滚珠数

，节圆直径

，滚珠直径

，接触角

，转频

计算得内圈故障特征频率：

3.3.2 仿真信号构造

模拟内圈故障的振动信号通常表现为周期性冲击衰减振荡，被转频调制，并叠加环境噪声。构造信号如下：

其中

为故障冲击周期；

为系统固有频率；=800为衰减系数；

为随机幅值，模拟冲击强度的微小波动；

为微小随机滑移，模拟实际转速波动；n(t)为高斯白噪声，信噪比约0dB；采样频率

，采样时长1s。

3.3.3 环境配置

依赖库：numpy、scipy、matplotlib，以及VMD的Python实现库vmdpy。

3.3.4 完整代码

import numpy as np import matplotlib.pyplot as plt from scipy.signal import hilbert from scipy.fft import fft, fftfreq from vmdpy import VMD plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False # ================== 1. 仿真内圈故障信号 ================== def simulate_bearing_fault(Fs=12000, duration=1.0, fr=30, fn=4000, beta=800): n_balls = 9 D = 39e-3 d = 7.9e-3 phi = 0 bpfi = n_balls / 2 * fr * (1 + d / D * np.cos(phi)) Tf = 1 / bpfi t = np.arange(0, duration, 1/Fs) N = len(t) num_impulses = int(duration / Tf) + 2 impulses = np.zeros(N) for i in range(num_impulses): t0 = i * Tf + np.random.normal(0, 1e-5) idx = (t >= t0) & (t < t0 + 5/beta) impulses[idx] += (np.sin(2 * np.pi * fn * (t[idx] - t0)) * np.exp(-beta * (t[idx] - t0)) * (0.8 + 0.2 * np.random.rand())) modulation = 1 + 0.5 * np.sin(2 * np.pi * fr * t) signal = impulses * modulation noise = np.random.randn(N) * 0.6 * np.std(signal) signal += noise return t, signal, bpfi, fr t, signal, bpfi, fr = simulate_bearing_fault() Fs = 12000 # ================== 2. VMD分解 ================== alpha = 2000 tau = 0.0 K = 4 DC = 0 init = 1 tol = 1e-7 u, u_hat, omega = VMD(signal, alpha, tau, K, DC, init, tol) if omega.ndim > 1: omega_final = omega[-1, :] # 最后一行 else: omega_final = omega if u.shape[0] == K and u.shape[1] != K: u = u.T # 转置为 (N, K) # 打印中心频率（转换为 Hz） center_freqs = omega_final * Fs / (2 * np.pi) print("各模态中心频率 (Hz):", np.round(center_freqs, 2)) # ================== 3. 选择最佳模态 ================== corr_coef = [np.corrcoef(signal.ravel(), u[:, k].ravel())[0, 1] for k in range(K)] best_imf_idx = np.argmax(np.abs(corr_coef)) print(f"与原始信号相关系数: {[round(c, 3) for c in corr_coef]}") print(f"选择 IMF{best_imf_idx+1} 进行包络分析") selected_imf = u[:, best_imf_idx] # ================== 4. 包络解调与包络谱 ================== analytic_signal = hilbert(selected_imf) envelope = np.abs(analytic_signal) envelope -= np.mean(envelope) N = len(envelope) freq = fftfreq(N, 1/Fs)[:N//2] env_spectrum = np.abs(fft(envelope))[:N//2] / N # ================== 5. 特征频率提取 ================== def find_peak(freq, spectrum, target_freq, bandwidth=5): idx = np.where((freq > target_freq - bandwidth) & (freq < target_freq + bandwidth))[0] if len(idx) == 0: return None, None peak_val = np.max(spectrum[idx]) peak_freq = freq[idx[np.argmax(spectrum[idx])]] return peak_freq, peak_val bpfi_peak, bpfi_amp = find_peak(freq, env_spectrum, bpfi) fr_peak, fr_amp = find_peak(freq, env_spectrum, fr) print(f"理论内圈故障频率: {bpfi:.2f} Hz") if bpfi_peak: print(f"包络谱检测频率: {bpfi_peak:.2f} Hz (相对误差: {abs(bpfi_peak-bpfi)/bpfi*100:.2f}%)") print(f"转频检测频率: {fr_peak:.2f} Hz" if fr_peak else "转频未检测到") # ================== 6. 可视化 ================== plt.figure(figsize=(14, 12)) plt.subplot(4, 1, 1) plt.plot(t, signal, linewidth=0.5) plt.title("原始仿真内圈故障信号") plt.xlabel("时间 (s)") plt.ylabel("幅值") plt.xlim(0, 0.2) plt.grid(alpha=0.3) for k in range(K): plt.subplot(4, 2, 3 + k) plt.plot(t, u[:, k], linewidth=0.5) plt.title(f"IMF {k+1} (中心频率 {center_freqs[k]:.1f} Hz)") plt.xlabel("时间 (s)") plt.ylabel("幅值") plt.xlim(0, 0.2) plt.grid(alpha=0.3) plt.tight_layout() plt.figure(figsize=(12, 5)) plt.subplot(1, 2, 1) plt.plot(t, selected_imf, 'b', linewidth=0.5, label='Selected IMF') plt.plot(t, envelope, 'r', linewidth=0.8, label='Envelope') plt.title("选取的IMF及其包络") plt.xlabel("时间 (s)") plt.ylabel("幅值") plt.xlim(0, 0.1) plt.legend() plt.grid(alpha=0.3) plt.subplot(1, 2, 2) plt.plot(freq, env_spectrum, 'k') plt.title("包络谱分析") plt.xlabel("频率 (Hz)") plt.ylabel("幅值") plt.xlim(0, 500) plt.grid(alpha=0.3) if bpfi_peak: plt.axvline(bpfi_peak, color='r', linestyle='--', label=f'检测 BPFI: {bpfi_peak:.1f} Hz') plt.axvline(bpfi, color='b', linestyle=':', label=f'理论 BPFI: {bpfi:.1f} Hz') if fr_peak: plt.axvline(fr_peak, color='g', linestyle='--', label=f'转频: {fr_peak:.1f} Hz') plt.legend() plt.tight_layout() plt.show()

3.4 实验结果与分析

3.4.1 VMD分解结果

采用VMD对仿真内圈故障信号进行分解，预设模态数K=4，惩罚因子=2000。各模态的最终收敛中心频率及与原始信号的相关系数如表1和图1所示。

表1 各IMF的中心频率与相关系数

图1 各IMF中心频率的可视化展示

从表1和图1可以看出：IMF1中心频率135.02 Hz，接近理论故障特征频率162.35 Hz，但相关系数较低（0.263），说明其包含了一定的故障信息，但可能受到噪声或其它成分的干扰，波形与原始信号的整体相似度不高。IMF2中心频率523.59 Hz，位于系统固有频率（4000 Hz）以下的中频段，相关系数中等。IMF3与IMF4中心频率分别为617.06 Hz和651.17 Hz，集中在较高频段，且相关系数依次增大（0.585与0.641）。这表明高频区模态承载了更多的信号能量和冲击特征，与仿真信号中系统固有频率4000 Hz经调制后的频带分布趋势相符。

根据相关系数最大化原则，选择IMF4作为承载故障冲击信息的最优模态，用于后续包络解调分析。

3.4.2 包络谱分析与故障特征提取

对IMF4进行希尔伯特变换提取包络，并计算包络谱，结果如图2所示。

图2包络谱分析图

在包络谱中搜索理论故障特征频率

及转频

，得到如下诊断结果：

1. 检测到的故障特征频率：162.00 Hz
2. 理论计算值：162.35 Hz
3. 相对误差：0.21%
4. 检测到的转频：30.00 Hz

包络谱中清晰呈现出162.00 Hz的谱峰，与理论BPFI值的误差仅为0.21%，远低于工程诊断中通常允许的±2%容差，有力地验证了故障特征频率的准确提取。同时，包络谱中明显检测到30.00 Hz的转频成分，这一现象完全符合内圈故障的物理特性：内圈随轴旋转，滚动体每经过缺陷位置产生的冲击强度会受到转频的正弦调制，因此在包络谱中除BPFI外还会出现转频和边频带。

3.4.3 VMD模态物理意义分析

结合各模态中心频率分布及诊断结果，可以进一步阐释VMD各分量的物理意义：

1. IMF1（135.02 Hz）在频率上紧邻故障特征频率，可能捕获了调制后的低频包络成分，但由于信噪比较低（仿真信号叠加了0 dB白噪声），其波形中噪声干扰相对严重，导致整体相关系数偏低；
2. IMF2–IMF4的中心频率均集中在500–700 Hz频带，远低于系统固有频率4000 Hz。这是因为仿真信号中高频振荡被转频调制和随机滑移展宽后，能量向中频段扩展，VMD自适应地将不同频段的冲击响应分离为多个模态。其中IMF4（651.17 Hz）能量最为集中，相关系数最高，说明其表征了冲击响应的主要频带成分，最适合用于包络解调。

这一分解结果表明，VMD成功将复杂的轴承故障信号解耦为具有明确物理意义的子带信号，有效避免了EMD中常见的模态混叠问题，同时展现出良好的噪声鲁棒性。

• 结论

基于VMD-包络谱的故障诊断流程，成功从强噪声背景下提取了内圈故障特征频率，检测值为162.00 Hz，与理论值误差仅0.21%，并检测到调制转频30.00 Hz，确诊为内圈故障。微小误差可能来源于：

1. 仿真时引入的随机滑移导致冲击间隔的细微抖动，使实际平均重复频率略微偏离理论值；
2. 频谱分辨率的限制，造成峰值频率估计存在量化误差；
3. 噪声干扰下峰值提取的微小波动。

尽管如此，0.21%的误差仍处于极低水平，充分验证了VMD在滚动轴承故障诊断中的有效性和精确性。该方法可进一步推广至实际工程振动信号的分析与自动诊断。