【实践任务】基于VMD的滚动轴承故障诊断
实践任务基于VMD的滚动轴承故障诊断研究背景与意义信号分解技术的发展历程信号处理是工程科学中的核心技术领域。传统的傅里叶分析技术以正弦函数为基函数适用于线性和稳态信号的分析。然而傅里叶变换在分析时变非线性信号时存在无法表述信号时频局部特性的局限性。为处理非平稳信号研究者相继发展了短时傅里叶变换STFT、Wigner-Ville分布、小波分析等方法。这些方法虽然在不同程度上改进了傅里叶分析的性能但本质上仍属于全局分析方法其分析性能取决于基函数的预先选取存在一定的适应性问题。自适应分解方法1998年美国工程院院士黄锷Norden E. Huang等人创造性地提出了经验模态分解方法被认为是2000年以来以傅里叶变换为基础的线性和稳态频谱分析的一个重大突破。EMD方法的革命性在于它是数据驱动的、自适应的无需预先设定任何基函数而是依据数据自身的时间尺度特征来进行信号分解。这一点与建立在先验性谐波基函数上的傅里叶分解和建立在先验性小波基函数上的小波分解方法具有本质性的差别。2014年Dragomiretskiy和Zosso提出了变分模态分解这是一种全新的非递归信号分解方法。与EMD的递归“筛选”模式不同VMD将信号分解转化为变分优化问题具有坚实的数学理论基础其实质是多个自适应维纳滤波组表现出更好的噪声鲁棒性。研究意义与应用价值EMD和VMD作为两种代表性的自适应信号分解方法在以下领域具有重要应用价值机械故障诊断通过对振动信号进行分解提取故障特征频率实现轴承、齿轮等关键零部件的状态监测与故障识别。生物医学信号处理用于心电图ECG、脑电图EEG等生理信号的分析与疾病识别。研究表明VMD在心电图分析中分类准确率可达98.30%。电力系统监测应用于电力负荷预测、电力电缆故障定位等场景具有较高的分析精度。地震资料处理用于地震信号的去噪和烃类检测等地球物理勘探任务。因此深入理解EMD和VMD的基本原理掌握其算法实现并探索其在工程实际中的具体应用具有重要的理论意义和实用价值。技术原理与方法2.1 经验模态分解EMD2.1.1 本征模态函数IMFEMD方法的核心是将复杂信号分解为有限个本征模态函数Intrinsic Mode Function, IMF和一个残余量。IMF必须满足以下两个条件在整个数据序列中极值点局部极大值和局部极小值的数量与过零点的数量相等或最多相差一个在任何时间点由局部极大值定义的上包络线和由局部极小值定义的下包络线的局部均值为零即上下包络线关于时间轴局部对称。2.1.2 EMD筛选过程EMD通过称为筛选Sifting的迭代过程完成信号分解具体步骤如下找出信号x(t)的所有局部极大值点和局部极小值点利用三次样条插值分别对极大值点和极小值点进行拟合构造上包络线和下包络线计算上下包络线的均值m(t)从原信号中减去包络均值得到检查是否满足IMF的两个条件。若满足则为第一个IMF分量否则将作为新的信号重复步骤1-4直至满足IMF条件或达到预设的终止准则从原信号中减去已提取的IMF分量得到残余信号将作为新的待分解信号重复步骤1-6依次分解出第二个、第三个等后续IMF分量。当残余信号为单调函数或变化足够小极值点数少于预定阈值时停止分解。最终原信号可表示为所有IMF分量与最终残余量之和。这一分解过程表明EMD具有完备性信号经分解后能够通过所有IMF分量及剩余分量被精确重构出来。2.1.3 EMD的主要问题EMD在实际应用中存在以下不足模态混叠当信号中存在由异常事件如间断信号、脉冲干扰和噪声等引起的间歇现象时同一个IMF中会包含不同尺度的信号成分或者同一尺度的信号出现在不同的IMF中导致IMF失去物理意义端点效应在三次样条插值拟合包络线时由于信号两端缺乏极值点约束导致拟合误差在迭代中被逐步放大并向内部传播。缺乏完整理论基础算法缺少完整的数学框架支撑在实际计算与应用中还存在筛分迭代停止标准等不确定性问题。为克服上述问题研究者提出了集合经验模态分解EEMD、完全集合经验模态分解CEEMD、自适应噪声完备集合经验模态分解CEEMDAN等多种改进方法。2.2 变分模态分解VMD2.2.1 基本思想VMD是由Dragomiretskiy和Zosso于2014年提出的一种全新的信号分解方法。与EMD的递归筛选模式根本不同VMD将信号分解问题转化为一个约束变分问题的求解过程。VMD建立在维纳滤波理论、希尔伯特变换理论和频率混合理论三大基础之上假设每个“模态”是具有不同中心频率的有限带宽信号其目标是使各模态的估计带宽之和最小。2.2.2 变分模型的构造VMD的分解过程可分为变分问题的构造和求解。模态定义VMD将每个模态定义为一个调幅-调频信号。变分问题的构建步骤如下对每个模态进行希尔伯特变换获得其解析信号的单边频谱将每个模态的频谱通过乘以指数项解调到相应的基频带通过解调信号的高斯平滑度来估计各模态的带宽。2.2.3 变分模型的求解为求解上述约束变分问题引入二次惩罚因子和拉格朗日乘子将约束问题转化为无约束问题。采用交替方向乘子法ADMM求解该扩展拉格朗日表达式的鞍点通过交替更新来逐步逼近最优解。2.2.4 VMD的优点VMD相比EMD具有以下显著优势坚实的数学理论基础VMD建立在变分优化和维纳滤波的严格数学框架之上从根本上克服了EMD缺乏理论支撑的缺点非递归分解模式VMD采用变分优化而非递归筛选进行分解有效避免了EMD的模态混叠问题更好的噪声鲁棒性其实质是多个自适应维纳滤波组对噪声具有天然的抑制能力更高的信号重构精度研究表明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(Fs12000, duration1.0, fr30, fn4000, beta800): 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_idx1} 进行包络分析) 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, bandwidth5): 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, linewidth0.5) plt.title(原始仿真内圈故障信号) plt.xlabel(时间 (s)) plt.ylabel(幅值) plt.xlim(0, 0.2) plt.grid(alpha0.3) for k in range(K): plt.subplot(4, 2, 3 k) plt.plot(t, u[:, k], linewidth0.5) plt.title(fIMF {k1} (中心频率 {center_freqs[k]:.1f} Hz)) plt.xlabel(时间 (s)) plt.ylabel(幅值) plt.xlim(0, 0.2) plt.grid(alpha0.3) plt.tight_layout() plt.figure(figsize(12, 5)) plt.subplot(1, 2, 1) plt.plot(t, selected_imf, b, linewidth0.5, labelSelected IMF) plt.plot(t, envelope, r, linewidth0.8, labelEnvelope) plt.title(选取的IMF及其包络) plt.xlabel(时间 (s)) plt.ylabel(幅值) plt.xlim(0, 0.1) plt.legend() plt.grid(alpha0.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(alpha0.3) if bpfi_peak: plt.axvline(bpfi_peak, colorr, linestyle--, labelf检测 BPFI: {bpfi_peak:.1f} Hz) plt.axvline(bpfi, colorb, linestyle:, labelf理论 BPFI: {bpfi:.1f} Hz) if fr_peak: plt.axvline(fr_peak, colorg, linestyle--, labelf转频: {fr_peak:.1f} Hz) plt.legend() plt.tight_layout() plt.show()3.4 实验结果与分析3.4.1 VMD分解结果采用VMD对仿真内圈故障信号进行分解预设模态数K4惩罚因子2000。各模态的最终收敛中心频率及与原始信号的相关系数如表1和图1所示。表1 各IMF的中心频率与相关系数模态中心频率 (Hz)相关系数IMF1135.020.263IMF2523.590.364IMF3617.060.585IMF4651.170.641图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包络谱分析图在包络谱中搜索理论故障特征频率及转频得到如下诊断结果检测到的故障特征频率162.00 Hz理论计算值162.35 Hz相对误差0.21%检测到的转频30.00 Hz包络谱中清晰呈现出162.00 Hz的谱峰与理论BPFI值的误差仅为0.21%远低于工程诊断中通常允许的±2%容差有力地验证了故障特征频率的准确提取。同时包络谱中明显检测到30.00 Hz的转频成分这一现象完全符合内圈故障的物理特性内圈随轴旋转滚动体每经过缺陷位置产生的冲击强度会受到转频的正弦调制因此在包络谱中除BPFI外还会出现转频和边频带。3.4.3 VMD模态物理意义分析结合各模态中心频率分布及诊断结果可以进一步阐释VMD各分量的物理意义IMF1135.02 Hz在频率上紧邻故障特征频率可能捕获了调制后的低频包络成分但由于信噪比较低仿真信号叠加了0 dB白噪声其波形中噪声干扰相对严重导致整体相关系数偏低IMF2–IMF4的中心频率均集中在500–700 Hz频带远低于系统固有频率4000 Hz。这是因为仿真信号中高频振荡被转频调制和随机滑移展宽后能量向中频段扩展VMD自适应地将不同频段的冲击响应分离为多个模态。其中IMF4651.17 Hz能量最为集中相关系数最高说明其表征了冲击响应的主要频带成分最适合用于包络解调。这一分解结果表明VMD成功将复杂的轴承故障信号解耦为具有明确物理意义的子带信号有效避免了EMD中常见的模态混叠问题同时展现出良好的噪声鲁棒性。结论基于VMD-包络谱的故障诊断流程成功从强噪声背景下提取了内圈故障特征频率检测值为162.00 Hz与理论值误差仅0.21%并检测到调制转频30.00 Hz确诊为内圈故障。微小误差可能来源于仿真时引入的随机滑移导致冲击间隔的细微抖动使实际平均重复频率略微偏离理论值频谱分辨率的限制造成峰值频率估计存在量化误差噪声干扰下峰值提取的微小波动。尽管如此0.21%的误差仍处于极低水平充分验证了VMD在滚动轴承故障诊断中的有效性和精确性。该方法可进一步推广至实际工程振动信号的分析与自动诊断。