搞滚动轴承故障诊断这行最烦的事就是现场采集回来的振动信号乱七八糟。转频、齿轮啮合、轴承冲击、环境噪声全混在一起直接做FFT频谱特征频率经常被噪声淹没或者跟别的频率成分叠在一起根本看不出来。后来我把VMD变分模态分解用到项目里才算是找到一套比较顺手的处理链路。这篇就把我把VMD用在滚动轴承故障诊断上的完整经验写出来包括原理、参数怎么调、代码怎么写、还有我踩过的坑给同样做信号处理的朋友做个参考。1. 为什么是VMD轴承故障信号的复杂程度超预期1.1 轴承振动信号里都有什么一个正常的滚动轴承振动信号里主要成分就是转频及低倍频。一旦轴承产生局部故障点蚀、剥落、裂纹滚动体经过缺陷位置时会产生周期性冲击。这种冲击有两个特点一是能量在极短时间集中释放对应到频谱上是宽频带激励二是冲击的幅值会受到转频或保持架频率的调制也就是幅值调制现象。我经常给刚入行的同事打比方轴承故障信号就像一个人站在嘈杂的集市里喊话喊的内容故障特征频率被周围的叫卖声齿轮啮合、转频谐波、随机噪声盖住了而且他喊的时候还忽大忽小调制。我们要做的就是把这个人的声音从集市里拎出来还得保住他喊的内容不跑调。1.2 传统方法的痛点早几年大家用EMD经验模态分解比较多。EMD是纯数据驱动的自适应分解方法不用预设基函数按频率从高到低层层剥离。听起来很完美用起来很痛苦模态混叠严重一个固有模态函数IMF里混着不同频率的成分有时同一个频率片段被拆到两个IMF里端点效应信号两端容易产生飞翼幅度能飞到吓人的程度缺乏数学理论支撑EMD的分解依据是极值点包络矩阵形式都给不出来论文里不太好写有人会针对模态混叠提出EEMD加入白噪声辅助但代价是计算量上去了残余噪声也上去了分解结果每次跑还不完全一样复现性是搞工程的痛点。VMD恰好在这几个问题上打了个翻身仗。1.3 VMD的基本思路VMD是2014年Dragomiretskiy和Zosso发表在IEEE Transactions on Signal Processing上的工作。核心思路是把信号分解变成一个约束变分问题的求解假设原始信号由K个有限带宽的模态分量组成每个模态围绕一个中心频率在频域上紧凑分布。求解目标就是让每个模态的带宽总和最小同时所有模态加起来能精确重构原始信号。这个思路相当于给每个频率成分划定势力范围谁的带宽窄谁就在自己的频带里待着不准越界这就从根源上抑制了模态混叠。而且整个求解过程有严谨的数学框架理论站得住脚。2. VMD的数学逻辑与核心参数2.1 变分模型到底在优化什么VMD把信号分解问题写成一个带约束的优化问题min Σₖ ‖ ∂ₜ [(δ(t) j/πt) * uₖ(t)] e^(-jωₖt) ‖²约束条件是所有模态加起来等于原始信号f(t)。每个uₖ(t)是分解出的一个模态分量对应一个IMFωₖ是这个模态的中心频率。∂ₜ是时间导数δ(t)j/πt对应的希尔伯特变换把信号变换成解析信号约束到单边频谱乘上e^(-jωₖt)把频谱搬移到基带然后对时间求导测带宽。整个过程的目标是让每个模态的频带尽量窄同时所有模态的叠加又能还原原始信号。2.2 求解过程简述这个约束问题用拉格朗日乘子法转化成无约束问题加入惩罚项α控制带宽约束的严格程度再引入对偶变量λ保证信号重构精度L({uₖ}, {ωₖ}, λ) α Σₖ ‖∂ₜ [(δ(t) j/πt) * uₖ(t)] e^(-jωₖt)‖² ‖f(t) - Σₖ uₖ(t)‖² ⟨λ(t), f(t) - Σₖ uₖ(t)⟩然后通过ADMM交替方向乘子法迭代求解。整个过程分为三步每一步固定其他两个变量更新当前变量更新模态uₖ在频域里是一个维纳滤波结构相当于把信号在中心频率ωₖ附近的成分滤出来更新中心频率ωₖ取当前模态频谱的质心位置作为新的中心频率更新对偶变量λ让最终重构信号和原始信号的差逼向零这样迭代几十次之后K个模态和对应的中心频率就稳定下来了。ADMM迭代对初值有一定敏感性但整体是收敛的。提示VMD的频域更新公式分母里有1 2α(ω - ωₖ)²本质上就是一个中心频率为ωₖ的自适应维纳滤波器。当α越大通带越窄频谱越锐利。2.3 参数含义全解VMD需要设置的参数一共有六个但真正影响结果的主要是三个参数常用值含义与影响K3~10模态分解个数。取小了欠分解取大了过分解。核心参数alpha默认2000带宽惩罚因子。控制每个模态的频带宽度tau0噪声容限。设为0表示信号分解过程中不额外考虑噪声项结果更干净DC0是否提取直流分量。轴承信号一般设0init1中心频率初始化方式。1为均匀初始化0为全零初始化tol1e-6收敛精度。默认即可我常用的写法是调用vmdpy库一行搞定from vmdpy import VMD u, u_hat, omega VMD(f, alpha2000, tau0, K5, DC0, init1, tol1e-7)返回的三个结果中u是分解出的K个模态每行一个模态的时域信号u_hat是模态的频域表达omega是每个模态的中心频率迭代收敛过程。主要用的是u和omega。3. 参数怎么选中心频率观测法与峭度准则3.1 K值选多少中心频率是照妖镜很多朋友上来就问我K怎么选说实话没有通用答案但有通用方法。VMD跑完后返回的omega是最有价值的参考信息。把K从小到大取跑完看每个模态的中心频率K取小了最后两个模态的中心频率相差很远中间明显有一段频带没有模态覆盖信号被欠分解K取合适各模态中心频率在频带上均匀分布没有中心频率特别接近的两个模态K取大了出现两个或以上模态的中心频率非常接近差值小于3倍频率分辨率这就是过分解出现了虚假模态实际项目里我的做法是先粗暴地把K从3取到9每个K值跑一次把中心频率列个表。如果加了K之后新增模态的中心频率跟现有模态挤在一起那就是K过头了。比如滚动轴承外圈故障数据K4时中心频率依次是98Hz、327Hz、982Hz、2658Hz间隔较大说明4个模态区分度好。3.2 alpha怎么调根据数据长度和频谱复杂度来alpha默认2000在大多数时候够用但如果数据点很少比如1秒采样率8kHz只有8000个点带宽就会分得比较粗糙这时可以适当提高alpha到3000~5000让模态更聚焦。如果信号本身频带很宽比如包含宽频冲击alpha可以降到1000左右。我的经验是alpha对结果的影响没有K那么敏感差个1000~2000结果基本稳定不用过度纠结。K选错了才是灾难性的。3.3 选哪个模态卷积峭度最大化准则VMD分解完不是把K个模态都拿去分析要挑出含有故障冲击成分的模态。冲击成分有一个特点峭度高。峭度是四阶统计量对冲击特别敏感正常振动信号的峭度接近3出现轴承故障时峭度能到20以上。最省事的办法对每个模态计算时域峭度取峭度最大的那个模态做包络谱分析。更严谨的做法是计算包络谱峭度因为即便模态时域峭度高包络谱里也不一定能看到清晰的故障特征频率。我一般两个都看先筛峭度大的再做包络谱验证。3.4 一个完整的参数选择流程我给自己定的流程是这样的先定带宽惩罚因子alpha2000中心频率初始化为均匀分布K从4开始跑一次打印中心频率看中心频率是否均匀覆盖主要频段有没有明显欠分解的迹象逐步增加K每次看新增模态的中心频率是否跟已有模态靠得太近确定K后若某个模态带宽过宽或过窄微调alpha重复跑一次计算各模态峭度锁定故障模态对故障模态做包络谱验证特征频率这套流程在工程上够用不用去碰那些复杂的优化算法。4. 滚动轴承故障诊断的完整代码实现4.1 先算故障特征频率轴承故障特征频率是诊断判别的基准。数学公式长这样外圈故障特征频率BPFO fr × N / 2 × (1 - d/D × cosα)内圈故障特征频率BPFI fr × N / 2 × (1 d/D × cosα)滚动体故障特征频率BSF fr × D / (2d) × (1 - (d/D × cosα)²)保持架故障特征频率FTF fr / 2 × (1 - d/D × cosα)其中fr是轴的转频N是滚动体数量d是滚动体直径D是轴承节圆直径α是接触角。这几个参数都可以从轴承手册或者网上公开的数据表里查。举个例子6205深沟球轴承滚动体数量9颗滚动体直径7.94mm节圆直径39.04mm接触角0度。轴转速1750rpm时转频约29.17Hz算出来外圈故障特征频率约104.25Hz内圈故障特征频率约158.38Hz滚动体故障特征频率约68.71Hz保持架故障特征频率约11.58Hz4.2 信号分解与包络谱分析拿到原始振动信号后完整流程是这样的对原始信号做预处理去除直流分量、均值消除趋势用VMD把信号分解成K个模态计算各模态时域峭度锁定故障模态对故障模态做Hilbert变换求包络信号对包络信号做FFT得到包络谱在包络谱里找故障特征频率及其2倍、3倍频为什么做包络谱而不是直接看频谱因为滚动轴承故障信号是冲击诱发的幅值调制信号调制源是故障特征频率。包络解调就是把调制信号解出来在包络谱里能看到清晰的调制频率峰这正是故障证据。直接在原信号频谱里看冲击能量散布在很宽的频带上特征频率周围全是旁瓣根本看不清。4.3 Python实操代码我用的是西储大学公开的轴承数据做演示采样率12kHz故障深度0.014英寸外圈故障转速1750rpm。完整代码如下import numpy as np import pandas as pd from scipy.signal import hilbert from scipy.fft import fft from vmdpy import VMD import matplotlib.pyplot as plt # 1. 加载数据 data pd.read_csv(outer_race.csv, headerNone) f data.iloc[:, 0].values # 第一列是振动加速度信号 fs 12000 # 采样率 # 2. 简单预处理去均值 f f - np.mean(f) # 3. VMD参数 alpha 2000 tau 0 K 6 DC 0 init 1 tol 1e-7 # 4. 执行VMD分解 u, u_hat, omega VMD(f, alpha, tau, K, DC, init, tol) # 5. 计算每个模态的峭度 def kurtosis(x): x x - np.mean(x) return np.mean(x**4) / (np.mean(x**2)**2) kurt_values [kurtosis(u[i, :]) for i in range(K)] fault_mode_idx int(np.argmax(kurt_values)) # 峭度最大的模态 # 6. 对故障模态做Hilbert包络解调 analytic_signal hilbert(u[fault_mode_idx, :]) envelope np.abs(analytic_signal) # 7. 包络谱分析 N len(envelope) freqs np.fft.rfftfreq(N, 1/fs) envelope_fft np.abs(fft(envelope))[:len(freqs)] # 8. 打印结果 print(f各模态峭度{np.round(kurt_values, 2)}) print(f故障模态索引{fault_mode_idx}) print(f故障模态中心频率{omega[-1, fault_mode_idx]:.2f} Hz)跑完这个代码我截图记录下来的典型结果是峭度最大的模态中心频率在2800Hz附近包络谱里清晰地看到104Hz、208Hz、312Hz三根谱线。对照BPFO104.25Hz三根线的间隔刚好一个故障特征频率这就是典型的外圈故障包络谱特征。4.4 结果解读的实操经验看到包络谱里的特征频率还要排除一个干扰如果只在1倍频处有峰有可能是转频谐波和故障特征频率碰巧接近但如果在1倍频、2倍频、3倍频处都有等间隔的峰那基本就是实锤故障。再退一步看原始信号的时域波形如果能够看到周期性的冲击冲击间隔约等于特征频率的倒数就可以交叉验证。如果包络谱里特征频率处没有峰或者峰值远低于底噪那这个轴承还处在正常状态或者故障非常早期不要硬从数据里找故障证据这是很多人容易误判的地方。5. 现场踩过的VMD深坑与排查链路5.1 过分解的假象K取大了的灾难体验我第一次上手VMD时想着模态多一点更精细把K直接设成8。结果分解出的好几个模态中心频率挤在200Hz附近其中两个模态波形长得几乎一样峭度还都特别高。我拿那个峭度最高的模态做包络谱谱线又杂又乱完全找不到故障特征频率。排查的过程是这样的先把模态两两做互相关发现有两个模态的互相关系数高达0.93这说明它们本质上是同一个成分被劈成了两半再看看它们的中心频率一个是96Hz一个是103Hz在频谱上相隔不到一个频率分辨率把K从8降到5重跑模态之间的相关系数都降到0.3以下边界清晰了包络谱里104Hz的故障线立刻跳出来了这就是教科书里常说的过分解导致虚假模态实际代码跑出来的表现就是中心频率撞车、模态波形高度相似、包络谱多出很多杂乱谱线。5.2 端点效应数据两端的伪冲击VMD虽然数学基础扎实但分解过程中同样存在端点效应。信号两端因为窗函数截断、滤波器边缘不连续会产生虚假的高低频成分。表现是分解出的第一个或最后一个模态在信号两端有明显的幅度振荡看起来像冲击。我的排查方法是把故障模态时域波形画出来看冲击是否均匀分布在全时段如果冲击只出现在两端而不是整个时间段都有那就可疑对信号做镜像延拓或多项式拟合延拓后再做VMD对比两端结果是否改善这个方法花不了多少时间但能避免很多误判。现场数据里那些两端飞起的大冲击很多时候都是处理出来的假象。5.3 采样率与数据长度的限制VMD有个容易被忽略的限制条件模态中心频率不能超过采样率的四分之一否则高频模态会失效。这是因为中心频率估计依赖频谱质心频谱上限是Nyquist频率超过之后信号频谱被折叠质心会往回偏移。另外数据太短也麻烦。VMD的频域滤波对频率分辨率有要求如果只有几百个采样点频域带宽分不开模态就会在低频谱段混成一团。我的经验是轴承故障诊断的数据长度至少要有40到60个完整的轴承转频周期19200Hz采样率下取个1秒以上的数据比较稳。5.4 跟EMD系列对比的失败与逆袭之前用EEMD处理同一组外圈故障数据分解出的IMF有9个其中三个都包含100Hz附近的成分模态混叠虽然减轻了但没完全消失包络谱里104Hz附近有几个靠得很近的谱峰分辨起来费力。VMD只用了6个模态每个模态都带明显的中心频率标签故障特征一目了然。这里不是要踩EEMD它有它的适应场景。但在滚动轴承这类频带分割需求明确的场景下VMD的边界清晰度和结果一致性确实好用。6. 功率分解视角VMD如何拆解频带能量6.1 功率分解的含义标题里有个说法叫VMD功率分解我理解是这个意思VMD可以把一个复杂信号按频带能量拆成多个模态每个模态近似代表一个频段内的信号分量。对每个模态做功率分析就能知道原始信号的能量主要分布在哪几个频带每个模态对总功率的贡献占比是多少。这在工程上很有用。正常轴承运行状态下信号能量主要集中在转频和低倍频段外圈故障发生后故障冲击激发的共振频带能量占比会明显上升且包含包络解调后的特征频率线。通过对比各模态的功率占比变化可以定位故障发生的位置和严重程度。6.2 功率占比的计算方法代码实现上很简单# 计算每个模态的功率占比 mse np.array([np.mean(u[i, :]**2) for i in range(K)]) power_total np.sum(mse) power_ratio mse / power_total print(f各模态功率占比{np.round(power_ratio, 4)})各模态功率占比干脆利落地告诉你哪一部分频带占了主导。比如外圈故障数据里第1个模态中心频率98Hz占12%第3个模态中心频率982Hz附近恰好接近轴承结构共振区占58%故障冲击能量大多被这个模态收走了。做故障分类时直接把各模态功率占比、峭度、中心频率做成一个特征向量喂给分类器效果比直接给原始波形好得多。6.3 一个真实的诊断案例数据来自某离心风机输入端滚动轴承转速980rpm采样率25600Hz。我用VMD分解K5alpha3000第1模态中心频率45Hz功率占比8%对应转频基频第2模态中心频率89Hz功率占比6%对应2倍转频第3模态中心频率221Hz功率占比15%是保持架和滚动体特征频率附近的频带第4模态中心频率1084Hz功率占比48%是轴承结构共振响应的主频带第5模态中心频率4212Hz功率占比23%是齿轮啮合频带对第4模态做包络谱在1084Hz谱峰两侧找到边频带边频间隔为7.9Hz。对照保持架故障特征频率7.86Hz判断为保持架磨损。后来拆机检查证实了判断。这个案例说明功率分解帮我们快速圈定重点频带包络解调则负责从中提取具体的故障证据两者配合效率很高。6.4 功率分解的局限性功率占比分析有一个注意点如果两个故障同时存在它们的冲击能量可能落在同一个共振频带里功率占比只能告诉你这个频带能量大但是分不清是外圈故障还是滚动体故障。这时候就必须看包络谱里的特征频率线。另外功率占比是相对量信号整体幅值变大时比如载荷增加即使轴承状态不变绝对功率也会上升。因此跨工况对比时要么保证工况一致要么把特征频率对应的功率单独归一化处理。我在做状态趋势监测时都是按同一个工况采集数据再比较功率占比的变化趋势。7. 不同工况下的参数调整与经验总结7.1 变转速工况轴承故障特征频率与转频成正比。转速变化时故障特征频率会来回漂移VMD的中心频率也会跟着移动。处理变转速信号时K要适当取大一些7~10否则转速升高时模态带宽不够部分频带会漏分析。alpha可以适当调大到3000让每个模态更聚焦但前提是模态中心频率不能跟转速变化范围撞上。7.2 强噪声环境现场环境噪声大尤其在化工厂、电厂等设备旁边。VMD在强噪声下还能工作但K的选择要保守一些。噪声频带没有明确的谱峰VMD会把噪声底也分割到某个模态里。我给这类信号的处理建议是先对原始信号做带通滤波抑制明显的外来干扰再用VMD做分解对分解结果做包络谱时只看故障特征频率及其倍频附近的窄带不要全谱观察7.3 模态数目的最终建议我这里给出一个我日常使用频率最高的参数组合K6alpha2000tau0init1tol1e-7。以这个为起点跑一次看中心频率中心频率均匀分布就直接用如果欠分解或者过分解再对K加减1重跑。这是效率和准确率的平衡点。8. 最后分享一点自己的体会做信号分解这几年最大感受是VMD不是万灵药但用对了绝对是好工具。很多人一上来就追求复杂算法什么优化K值、自适应alpha研究半天论文结果在实战里还不如老老实实把K从4试到8看中心频率分布来得直接。算法是为人服务的花三分钟跑个实验比花三小时去推导最优参数理论工程价值大得多。还有一个实用的技巧VMD一旦跑出结果先去看omega这个返回值模态太多、太少、撞频、偏置全部在中心频率上写着。你盯着中心频率看几眼比去琢磨各种无效指标都好使。这套流程我目前仍在用从实验室验证搬到现场数据上可靠性一直在线。轴承故障诊断本质上是证据链的分析过程VMD负责把原始信号整理成有清晰语义的模态分量包络谱负责从模态中提取特征频率。链条两侧都完整了判断就会准确。希望这篇能少让你走点弯路。