做轴承故障诊断这些年，我最大的体会是：算法模型排着队等着验证，可数据源却常常让人头疼。实测数据要去现场采集，成本高、工况杂，故障类型还不好控制。所以我很早就养成了一个习惯——先用基于故障机理的动力学仿真代码生成信号，像搭积木一样把正常、外圈故障、内圈故障、滚动体故障这四类典型工况全部模拟出来，再拿这些信号去训练算法、验证理论。
这篇文章要聊的，就是一套用MATLAB实现的滚动轴承动力学建模代码。它根据轴承的几何参数、转速和故障机理，直接仿真出振动信号，输出时域波形、频谱和包络谱，既能用于故障诊断算法的前期验证，也能支撑教学演示和论文仿真。不管你是刚接触故障诊断的研究生，还是被数据问题折腾的工程师，这套代码都能给你一个可控、可复现的实验环境，让你把精力放在算法本身。
1. 为什么需要轴承动力学仿真代码
1.1 故障诊断研究中的信号来源困局
做故障诊断的人都懂，实测数据是金标准，但获取成本极高。某高校实验室搭一个轴承试验台，电机、联轴器、轴承座、加载机构一套下来少说几万，还要人工加工故障——电火花、激光切割、线切割，工序麻烦不说，故障尺寸还很难精确控制。更头疼的是，同一型号轴承在不同转速、不同载荷下的故障响应差异很大，今天采集的数据和明天采集的波形完全对不上，实验重复性很差。
如果用纯信号合成来凑数，比如直接生成几个正弦波加白噪声，那又缺少物理意义。真实故障信号不是简单的正弦叠加，而是周期性冲击与系统共振调制后的复杂波形，冲击会激起轴承–转子系统的固有频率，产生衰减振荡，再叠加转频、谐波、边频带。纯合成信号做不到这种层次，拿去验证算法，很多时候是“看起来能用，一上真实数据就翻车”。
1.2 动力学建模相比纯信号合成的优势
这里说的动力学模型，是把轴承–转子系统简化成质量–弹簧–阻尼结构，用微分方程描述振动响应。故障在物理上表现为对系统的周期性冲击激励，不同的故障位置对应不同的冲击频率和幅值调制方式。用这套模型生成仿真信号，有几个明显好处。
第一是参数可调、规律可控。改变转速、载荷、故障尺寸，信号的频带宽度、调制特征会按物理规律变化，不是人为拼凑出来的。第二是特征丰富。仿真信号天然包含冲击衰减、共振调制、谐波和边频，拿去验证包络谱、小波分解、经验模态分解甚至深度学习模型，都很有说服力。第三是完全可复现。代码跑出来的结果一致，论文里写“参数如上、结果如下”经得起推敲，不用担心实验环境漂移。
2. 滚动轴承故障机理与特征频率计算
2.1 内外圈与滚动体故障的物理模型
滚动轴承主要由内圈、外圈、滚动体和保持架组成。正常状态下，滚动体在滚道里滚动，振动信号主要是低频的转频振动和微小的随机噪声。一旦某个部件表面出现剥落、裂纹或压痕，滚动体滚过缺陷处就会产生一个瞬时冲击，这个冲击的频率就是所谓的故障特征频率。
外圈故障时，缺陷位置相对传感器固定，冲击每隔一个外圈通过频率出现一次，间隔均匀。内圈故障时，缺陷随内圈旋转，冲击频率是内圈通过频率，而且冲击幅值会随载荷区位置的变化被转频调制。滚动体故障时，滚动体自转一周会与内外滚道各接触一次，产生两次冲击，冲击间隔等于滚动体故障特征频率，同时冲击幅值还会受到保持架旋转频率的调制。理解这些物理过程，是建模的基础。
2.2 特征频率公式与参数设定
以最常见的深沟球轴承为例，设节圆直径 (D)、滚动体直径 (d)、接触角 (\alpha)、滚动体个数 (z)、转频 (f_r)（单位Hz）。故障特征频率的计算公式如下：
- 外圈故障频率：(BPFO = \frac{z f_r}{2} \left(1 - \frac{d}{D} \cos\alpha\right))
- 内圈故障频率：(BPFI = \frac{z f_r}{2} \left(1 + \frac{d}{D} \cos\alpha\right))
- 滚动体故障频率：(BSF = \frac{D f_r}{2d} \left(1 - \left(\frac{d}{D} \cos\alpha\right)^2\right))
- 保持架故障频率：(FTF = \frac{f_r}{2} \left(1 - \frac{d}{D} \cos\alpha\right))
注意，转速给定通常用RPM，要换算成Hz需要除以60。比如转速1500RPM，转频就是(f_r = 1500 / 60 = 25)Hz。下面用一组典型深沟球轴承参数做示例：滚动体个数(z=9)，节圆直径(D=39)mm，滚动体直径(d=7.5)mm，接触角(\alpha=0^\circ)。
把参数代入公式算一遍：
- 外圈故障频率：(BPFO = \frac{9 \times 25}{2} \times \left(1 - \frac{7.5}{39}\right) = 112.5 \times 0.8077 \approx 90.9) Hz
- 内圈故障频率：(BPFI = 112.5 \times \left(1 + \frac{7.5}{39}\right) = 112.5 \times 1.1923 \approx 134.1) Hz
- 滚动体故障频率：(BSF = \frac{39 \times 25}{2 \times 7.5} \times \left(1 - \left(\frac{7.5}{39}\right)^2\right) = 65 \times 0.963 \approx 62.6) Hz
- 保持架故障频率：(FTF = 12.5 \times 0.8077 \approx 10.1) Hz
这个计算结果后面要直接写进代码里生成冲击序列。
3. MATLAB代码整体架构与核心模块
3.1 参数初始化与数据生成流程
整套代码采用主脚本加函数的模块化结构，主要分四步：定义轴承参数和工况参数、生成故障冲击序列、求解动力学方程得到振动响应、绘图输出。主脚本清晰列出所有可调参数，方便批量实验。
参数初始化要包括：采样率(f_s)、采样时长(T)、转速、载荷、系统固有频率、阻尼比、轴承几何参数。这些参数直接决定仿真信号的质量，必须放在一起方便统一调整。
数据生成流程的核心是冲击序列的构造。冲击序列是一串等间隔的脉冲，间隔就是对应故障特征频率的倒数。比如外圈故障，冲击周期就是(1/BPFO)，也就是每经过0.011秒产生一次冲击。内圈故障除了冲击间隔是(1/BPFI)，还要让冲击幅值随转频变化。滚动体故障要额外叠加保持架调制。
3.2 动力学方程离散化与求解
滚动轴承的振动系统可以简化成单自由度质量–弹簧–阻尼系统，动力学方程为：
[ m\ddot{x} + c\dot{x} + kx = F(t) ]
其中(m)是等效质量，(c)是阻尼，(k)是刚度，(F(t))是故障冲击力序列。固有频率(\omega_n = \sqrt{k/m})，阻尼比(\zeta = c/(2\sqrt{km}))。实际冲击会激起系统共振，所以响应表现为衰减振荡，这正是实测信号中常见的“冲击加衰减”波形。
我用四阶龙格库塔法对微分方程离散求解，采样频率设为(f_s)，时间步长(\Delta t = 1/f_s)。核心代码如下：
% 动力学响应求解主函数 function x = bearing_response(impact_train, fs, fn, zeta) dt = 1/fs; N = length(impact_train); x = zeros(1, N); v = 0; % 初始速度 x0 = 0; % 初始位移 wn = 2*pi*fn; % 数值积分 for i = 2:N F = impact_train(i); a = (F - 2*zeta*wn*v - wn^2*x0) / 1; % 质量归一化为1 v = v + a*dt; x0 = x0 + v*dt; x(i) = x0; end end这里把质量归一化为1，共振响应频率(\omega_n)直接表示系统固有频率。阻尼比(\zeta)通常取0.02到0.05，太大会把冲击衰减抹平，太小又会造成无限振荡，和真实信号不符。
3.3 四种工况的差异化建模
正常状态最简单：只包含一个低频转频正弦信号、一些随机白噪声和微小冲击。实际轴承正常工作时也有一点润滑搅动和结构噪声，但整体是平稳随机过程。
外圈故障：冲击序列是周期脉冲，周期为(1/BPFO)，幅值恒定。因为外圈固定，冲击传递路径不变，所以信号包络基本均匀。
内圈故障：冲击序列周期为(1/BPFI)，但冲击幅值必须乘上一个与转频同步的调制系数。调制系数用(\cos(2\pi f_r t))实现，模拟载荷区集中在内圈单点时的周期性强弱变化。
滚动体故障：冲击周期为(1/BSF)，理论上每个周期内产生两次冲击（分别接触内外圈），实现时可以用两个脉冲间隔半个周期的组合。同时由于滚动体故障位置随保持架旋转，冲击幅值还受(FTF)频率调制。
四种工况的差异都建立在同一个动力学模型上，只是冲击序列的生成规则不同。切换工况只要换一段生成冲击序列的函数，非常方便。
4. 实操过程：从运行到结果可视化
4.1 运行环境与文件结构
开发环境建议MATLAB R2016及以上版本，工具箱只用基础版和信号处理工具箱（Signal Processing Toolbox），不需要额外安装。文件结构建议如下：
main_bearing_sim.m：主脚本，定义参数、调用函数、绘图generate_impact.m：生成不同工况的冲击序列bearing_response.m：动力学响应求解函数freq_calculation.m：故障特征频率计算函数
这种拆分方式便于调试。主脚本控制全局，子函数各管一摊，换参数不用翻大文件。
4.2 参数调节经验与信号采样注意
采样率是第一个要关注的参数。故障特征频率本身不高，也就几十到几百赫兹，但冲击激起的共振可能在几千赫兹。按耐奎斯特采样定理，采样率至少要是最高分析频率的2倍。我通常设(f_s=20\text{kHz})或(50\text{kHz})，保证共振衰减段有足够采样点，后续做包络谱分析也游刃有余。
采样时长也要注意。至少包含20到30个冲击周期，才能看到稳定的频谱特征。以本文算例(BPFO=90.9)Hz，周期约0.011秒，采样0.2秒就够。但如果是滚动体慢转的情况，保持架频率才10Hz，就得多采几秒数据，不然看不到完整的调制周期。
阻尼比的经验值：轻载时取0.01到0.02，重载时取0.03到0.05。这个参数决定了冲击衰减速度，直接影响时域波形看上去“清不清晰”。调模型的时候先固定其他参数，单看时域波形，慢慢把阻尼比调到既有明显冲击、又不过度振荡的状态。
4.3 时域图、包络谱与特征提取
运行完主脚本后，直接输出四张图：时域波形、幅值谱、包络谱和时频图。时域波形上，外圈故障能看到均匀的周期冲击；内圈故障能看到幅值忽大忽小的冲击；滚动体故障能看出成对出现的冲击特征。
包络谱是验证故障频率的关键。对时域信号做希尔伯特变换求包络，再对包络信号做FFT，能在包络谱中清晰看到故障特征频率及其谐波。我用这段代码实现包络谱：
% 包络谱计算 y_env = abs(hilbert(x)); % 希尔伯特包络 y_env_detrend = y_env - mean(y_env); Y_env = fft(y_env_detrend); f_axis = (0:length(Y_env)-1) * fs / length(Y_env); plot(f_axis(1:fs/2), abs(Y_env(1:fs/2)));把频谱里的幅值峰值和计算得到的故障特征频率对照，能直接验证代码是否对了。以内圈故障为例，包络谱峰值应出现在(BPFI=134.1)Hz，同时能看到(2BPFI)、(3BPFI)等谐波，以及被转频(f_r)调制的边频带。
5. 常见问题与排查技巧
5.1 故障特征频率算不对怎么办
这是最常踩的坑。第一，单位一定要换算。转速给的RPM，转频要除以60。第二，接触角不能凭感觉忽略。圆柱滚子轴承接触角为0没问题，角接触球轴承一定要把实际接触角代进公式。第三，滚动体故障频率的公式里有一个平方项，手算时很容易漏掉。第四，内圈故障的调制信号频率是转频(f_r)，不是保持架频率，别搞混。
排查方法是先单独把特征频率计算模块跑一遍，给定一组已知参数，用公式手算对比，确认无误后再去生成信号。我习惯把这几个频率值打印在命令行窗口里，一边看信号一边对照，一眼就能看出毛病。
5.2 仿真信号与实测信号差距大怎么解决
很多新手生成时域信号后，发现冲击就是一个个孤立脉冲，整体看起来像“棒棒糖”，和实测那种连续振荡的波形完全不同。原因在于冲击激励不能只是单点脉冲，而应该激起系统的衰减振荡。
解决方法是给冲击序列加上一个共振衰减包络。每个冲击不是理想脉冲，而是一个高频振荡乘以指数衰减项。实现时可以用一个冲击响应函数模板，和冲击序列做卷积：
% 冲击响应模板 impulse_t = 0:1/fs:0.005; % 观察5ms h = exp(-zeta*wn*impulse_t) .* sin(wn*impulse_t); % 用卷积模拟共振响应 x_response = conv(impact_train, h, 'same');这样生成出来的信号，冲击后面拖着一串衰减尾巴，时域波形立刻就有实感了，包络谱的信噪比也大幅提升。
5.3 代码运行效率优化建议
仿真时长较长时，循环求解速度可能很慢。第一优先考虑用卷积替代循环，把动力学响应看作线性系统对冲击序列的响应。第二用滤波法，将冲击序列通过一个二阶带通滤波器，也能模拟共振环节。第三，如果必须用循环，也要在MATLAB里尽量避免逐点赋值，改用矩阵运算或并矢操作。
还有一个小技巧：可以把高频共振响应做成分段常系数状态空间模型，再利用MATLAB的lsim函数求解，速度比手写积分快得多。不过对大部分故障诊断实验来说，原始四阶龙格库塔已经够用，不必过度优化。
6. 个人实操心得与扩展方向
我把这套代码用在某设备健康监测项目的数据增强环节，效果很明显。原始实测故障样本不够，我就用仿真信号扩充训练集，特别是在卷积神经网络里做迁移学习，先用仿真数据预训练，再用少量实测数据微调，识别准确率提升了一大截。可见这套模型的价值不止在“画波形好看”，它确实能作为真实系统的近似替身。
不过要说清楚，仿真永远替代不了实测。动力学模型把系统简化成了线性时不变结构，但真实轴承存在间隙、非线性接触刚度和冲击耦合，这些因素在极端工况下会占据主导。我的经验是：用仿真信号做算法预研、做理论验证是靠谱的，但最后一定要用实测数据做终检，两者结合才是稳妥的工程路径。
最后再分享一个小技巧：代码里可以加一个“故障严重度”参数，通过控制冲击幅值和冲击宽度的比例来模拟早期和晚期故障。早期故障冲击幅值小、脉冲窄，晚期故障冲击幅值大、脉冲宽。调整这个参数再配合不同信噪比，就能生成一个故障演化梯度。想研究寿命预测和退化建模的朋友，这比到处找实测退化数据方便得多。后面有时间，我准备把变转速和变荷载的工况加进去，让这套仿真环境真正贴近工业现场的复杂状态。
