1. 为什么必须先把轴承“算明白”故障建模的工程价值先说个让我印象挺深的场景。之前帮一个做旋转机械状态监测的团队处理过一个问题他们采集了大量真实轴承振动数据喂给一个卷积神经网络做故障分类实验室验证精度接近98%但一上现场就掉到80%出头。折腾了两个月最后发现瓶颈根本不在算法而在训练数据——现场工况的转速波动、负载变化、噪声水平远比实验室“干净”数据复杂模型见过的变化太少泛化能力自然上不去。那会儿我们做的一件事就是用MATLAB把轴承动力学模型搭起来先仿真生成各种工况下的故障振动信号跟实测数据混合后重新训练模型最后现场精度提到了92%。这次经历让我彻底意识到一件事轴承故障仿真不是学术圈自娱自乐它是工程诊断算法落地过程中非常关键的一环。这篇东西我想围绕“轴承动力学”这条线完整走一遍从故障机理分析、特征频率计算、动力学建模到用MATLAB生成仿真振动数据的全过程。适合这几类读者做设备状态监测、故障诊断算法开发的工程师需要大量带标签数据但实测采集成本太高研究旋转机械动力学的学生想把理论公式变成可跑的仿真程序刚接触MATLAB/Simulink想用一个完整案例串联信号处理、数值积分、数据分析的初学者。你不需要先成为动力学专家只要懂一点振动的基本概念跟着把代码跑通就能生成很接近真实传感器信号的轴承故障数据。2. 故障特征频率一切仿真的起点2.1 轴承几何参数与四个特征频率轴承故障为什么能在振动信号里被识别出来核心机理是滚动体在滚道上运动时一旦经过局部损伤比如外圈剥落、内圈裂纹、滚动体点蚀就会产生周期性的冲击激励激起轴承座和传感器安装结构的固有频率振动。这个“周期性”不是随便来的它由轴承的几何尺寸和转速严格决定。我们常说的特征频率有四类故障类型特征频率公式物理含义外圈故障BPFOBPFO (n/60) × (N/2) × (1 - (d/D) × cos α)滚动体经过外圈某一损伤点的频率内圈故障BPFIBPFI (n/60) × (N/2) × (1 (d/D) × cos α)滚动体经过内圈损伤点的频率滚动体故障BSFBSF (n/60) × (D/(2d)) × (1 - (d/D)² × cos² α)滚动体自转频率局部缺陷每次接触滚道产生两次冲击保持架故障FTFFTF (n/60) × 0.5 × (1 - (d/D) × cos α)保持架旋转频率通常较低公式里的参数含义n转轴转速rpmN滚动体个数d滚动体直径mmD节圆直径mmα接触角度。这里有个特别容易搞混的点**内圈故障的特征频率不是简单的“转频乘以某倍数”它还会被转频调制。**因为内圈随轴旋转损伤点相对于载荷区周期性进出导致冲击幅值时大时小。所以仿真内圈故障时光有周期冲击还不够必须叠加上转频调制否则生成信号的包络谱会缺关键边频带跟实测对不上。2.2 用MATLAB算一组真实参数我常用SKF 6205-2RS深沟球轴承做演示因为公开论文里大量使用这个型号比如著名的CWRU轴承数据集参数齐全方便对照验证。滚动体个数 N 9滚动体直径 d 7.94 mm节圆直径 D 39.04 mm接触角 α 0°深沟球轴承近似按0算在MATLAB里定义一个脚本% 轴承参数 N 9; % 滚动体个数 d 7.94; % 滚动体直径 mm D 39.04; % 节圆直径 mm alpha 0; % 接触角度 n_rpm 1750; % 转速 rpm fr n_rpm / 60; % 转频 Hz % 特征频率计算 cos_alpha cosd(alpha); BPFO fr * (N / 2) * (1 - (d / D) * cos_alpha); BPFI fr * (N / 2) * (1 (d / D) * cos_alpha); BSF fr * (D / (2 * d)) * (1 - (d / D)^2 * cos_alpha^2); FTF fr * 0.5 * (1 - (d / D) * cos_alpha); fprintf(转频 fr %.2f Hz\n, fr); fprintf(外圈 BPFO %.2f Hz\n, BPFO); fprintf(内圈 BPFI %.2f Hz\n, BPFI); fprintf(滚动体 BSF %.2f Hz\n, BSF); fprintf(保持架 FTF %.2f Hz\n, FTF);转速1750 rpm下计算结果大致是转频 fr ≈ 29.17 HzBPFO ≈ 104.56 HzBPFI ≈ 157.94 HzBSF ≈ 68.72 HzFTF ≈ 11.65 Hz外圈故障频率约等于9倍转频的区域内圈故障接近5.4倍转频滚动体故障居中。拿到这些数值后面才有基准去验证仿真信号的包络谱峰值是否落在正确位置。2.3 参数敏感性为什么小数第一位都影响判断很多人觉得特征频率算个大概就行实际上工程里差之毫厘谬以千里。原因在于包络谱的频率分辨率 采样频率 / FFT点数如果你的采样时长只有1秒频率分辨率只有1 HzBPFO的预测值差0.5 Hz会直接导致谱峰偏移造成误判。另外轴承实际工作时会有温升导致内外圈和滚动体热膨胀节圆直径和接触角都会微变。重载低速工况下接触角可能偏移实测特征频率和理论值之间常有1%2%的偏差。所以行业里的做法是理论值算出来只是参考诊断时在理论值附近搜索局部谱峰。仿真时也要考虑这个偏差。我通常在理论上加一个微小的随机频率扰动模拟工况波动这样的数据训练出来的模型抗干扰能力更强这个后面章节细说。3. 从单自由度到分段Hertz接触动力学模型搭建3.1 单自由度模型最基础的“质量-弹簧-阻尼”特征频率解决的是“冲击多久来一次”动力学模型解决的是“冲击来了系统怎么响应”。最简单也最常用的是单自由度振动模型m·x c·x k·x F(t)m等效质量轴承座和传感器结构贡献c阻尼系数k接触刚度x振动位移F(t)外部激励力这个模型物理意义很直观传感器测到的振动本质上是结构受到激励后的受迫响应。故障冲击是激励源轴承座和传感器安装座的固有振动特性决定信号长什么样。MATLAB里用ode45解这个二阶微分方程。先把方程改写成一阶状态空间x₁ x x₂ x则x₁ x₂ x₂ (F(t) - c·x₂ - k·x₁) / m% 单自由度轴承振动模型 clear; clc; % 系统参数 m 0.6; % 等效质量 kg k 2.5e7; % 等效刚度 N/m zeta 0.02; % 阻尼比 c 2 * zeta * sqrt(m * k); % 阻尼系数 fs 51200; % 采样率 51.2kHz工业采集卡常用 t_end 2; % 仿真时长 2s t (0:1/fs:t_end); fr 29.17; % 转频 Hz BPFO 104.56; % 外圈故障特征频率 Hz % 构造冲击力序列简化为狄拉克脉冲序列 impulse_interval 1 / BPFO; impulse_times 0:impulse_interval:t_end; F zeros(size(t)); for ii 1:length(impulse_times) [~, idx] min(abs(t - impulse_times(ii))); F(idx) 100; % 每个冲击力幅值 100N end % 转换为状态空间求解 [t_ode, y] ode45((t, y) single_dof_ode(t, y, m, c, k, F, t), ... [0 t_end], [0 0]); x y(:, 1); % 位移 v y(:, 2); % 速度 % 重采样到等间隔时间向量 x_resamp interp1(t_ode, x, t); v_resamp interp1(t_ode, v, t); figure; subplot(2,1,1); plot(t, x_resamp); xlim([0 0.2]); title(位移响应外圈故障单自由度模型); xlabel(时间 (s)); ylabel(位移 (m)); subplot(2,1,2); % 包络谱验证 env abs(hilbert(v_resamp)); Nfft 8192; win hann(Nfft); [pxx, f] pwelch(env, win, Nfft/2, Nfft, fs); plot(f, 10*log10(pxx)); xlim([0 500]); title(速度包络谱); xlabel(频率 (Hz)); ylabel(功率谱 (dB)); % 局部函数 function dydt single_dof_ode(t, y, m, c, k, F, t_grid) F_interp interp1(t_grid, F, t, linear, 0); dydt zeros(2,1); dydt(1) y(2); dydt(2) (F_interp - c*y(2) - k*y(1)) / m; end运行后位移响应是典型的衰减振荡波形每次冲击后结构以固有频率 ≈ sqrt(k/m) / (2π) 振动并逐渐衰减。包络谱上能看到BPFO处的峰值。这个模型最大的问题**冲击力被简单当成等幅脉冲序列没有考虑滚动体与损伤区域的真实接触过程。**实际故障冲击的力幅和波形跟Hertz接触理论直接相关不是简单给个常数。3.2 分段Hertz接触刚度让模型贴近物理Hertz接触理论描述两个弹性体接触时的应力与变形关系。对轴承来说滚动体和滚道接触时法向接触力与变形量的关系是F K_h · δ^(3/2)K_hHertz接触刚度系数δ接触变形量弹性趋近量这个关系是非线性的意味着等效刚度并不是常数而是随变形量变化。仿真时需要在每个时间步根据当前位置计算接触力不能再用固定刚度k。更关键的是当滚动体进入损伤区域时接触状态会突变。损伤区域内的接触刚度显著下降表现为一个瞬时的“卸力—再冲击”过程。这也是为什么真实故障振动信号里冲击波形不是标准正弦或指数衰减而是带有复杂高频成分的瞬态波形。用MATLAB实现分段Hertz接触模型思路是计算每个滚动体相对于损伤点的角位置判断是否落入损伤区域损伤角宽度落入区域前正常接触落入后按退化刚度计算接触力对所有滚动体的接触力求和作为总激励。% 外圈故障的分段Hertz接触仿真 clear; clc; % 轴承与工况参数 N 9; d_ball 7.94e-3; % 滚动体直径 m D_pitch 39.04e-3; % 节圆直径 m m_eff 0.6; % 等效质量 kg zeta 0.02; fs 51200; t_end 1; t (0:1/fs:t_end); n_rpm 1750; fr n_rpm / 60; % 外圈特征频率与损伤参数 BPFO_Hz fr * (N/2) * (1 - d_ball/D_pitch); damage_angle 2; % 损伤角宽度度 damage_center 0; % 损伤中心在角度0处 % Hertz接触刚度按钢制轴承近似 K_h 1.2e10; % N/m^(3/2)经验值 % 模拟每个滚动体的角位置 dphi 2*pi / N; phi0 0; % 主循环计算等效接触力 F_contact zeros(size(t)); delta0 1e-5; % 预紧变形量 m for ii 1:length(t) shaft_angle mod(2*pi*fr*t(ii), 2*pi); F_total 0; for jj 1:N phi_j phi0 shaft_angle (jj-1)*dphi; phi_j mod(phi_j, 2*pi); % 判断是否落入损伤区域外圈固定损伤中心角度已知 delta_phi abs(phi_j - deg2rad(damage_center)); delta_phi min(delta_phi, 2*pi - delta_phi); if delta_phi deg2rad(damage_angle/2) % 落入损伤区域接触刚度退化 K_local 0.15 * K_h; else K_local K_h; end % 变形量简化假设每个滚动体承载相同但实际有载荷分布 delta delta0; F_total F_total K_local * delta^1.5; end F_contact(ii) F_total; end % 加速度响应等效为F/m accel F_contact / m_eff; figure; subplot(2,1,1); plot(t, F_contact); xlim([0 0.1]); title(接触力时域波形外圈故障); xlabel(时间 (s)); ylabel(接触力 (N)); subplot(2,1,2); env abs(hilbert(accel)); [pxx, f] pwelch(env, hann(8192), 4096, 8192, fs); plot(f, 10*log10(pxx)); xlim([0 300]); title(加速度包络谱); xlabel(频率 (Hz)); ylabel(功率谱 (dB));这个模型比单自由度朴素版真实得多因为当滚动体经过损伤区时接触力会明显跌落再恢复形成特有的“卸荷—冲击”波形包络谱中的BPFO谐波更加丰富。3.3 多自由度扩展轴承-转子-基座耦合单自由度模型在信号形态上够用但如果你需要研究“哪个传感器位置能更好捕捉故障特征”“不同测点之间信号差异有多大”就得把模型扩展到多自由度。典型的做法是建立轴承-转子-基座耦合模型大致结构是转子轴用一个转动惯量和两个支撑轴承表示每个轴承用两个方向的刚度-阻尼单元表示水平、垂直基座用低刚度低阻尼的附加质量块表示模拟传感器安装结构故障冲击通过Hertz接触模型施加在对应轴承的内部节点上。这个系统用状态空间描述自由度数量从1增加到610求解时间明显变长但能模拟出水平/垂直方向振动差异、测点距离导致的衰减和相位差这些信息对传感器布局选型很有参考价值。MATLAB里用Simulink搭建这种模型比手写ode45更直观。每个轴承用一个“Translational Spring-Damper”模块转子用“Inertia”模块故障冲击用“MATLAB Function”模块自定义Hertz接触力。跑一遍后可以用“Scope”看波形也可以把数据导回工作区做包络谱分析。Simulink的优势是改参数方便改个刚度、加个测点改连线就行不用每次都改代码。4. 故障冲击注入与振动信号模拟从“像”到“是”4.1 为什么要单独做信号级模拟动力学模型生成的是物理层面的振动响应但传感器实际采集到的信号还叠加了很多环节的影响传感器的频响特性、传输路径衰减、电磁噪声、工频干扰、转速波动带来的频率漂移等。动力学模型直接输出拿来当训练数据跟实测数据之间的域差异还是偏大。所以工程上常用“混合建模”思路用动力学模型生成冲击响应波形作为信号模板再用信号处理手段把它调制到完整的振动信号中这样可以灵活控制噪声、调幅、畸变等参数批量生成多样本。4.2 故障振动信号的完整生成流程一个较完整的外圈故障仿真信号流程如下生成单位冲击响应 → 按特征频率排布冲击序列 → 乘上周期性幅值调制 → 叠加转频谐波成分 → 叠加高斯白噪声 → 加传感器频响滤波 → 输出第3章已经讲过怎么用动力学模型生成冲击响应。如果有现成的实测故障冲击片段也可以直接从实测数据里截取一段作为模板这样的仿真信号保真度更高。单位冲击响应用MATLAB生成指数衰减正弦波作为单位冲击响应模板% 单位冲击响应模板 fn 2400; % 结构固有频率 Hz典型轴承座频率 zeta_c 0.03; % 阻尼比 tau 1/(2*pi*fn*zeta_c); % 衰减时间常数 T_impulse 0.005; % 冲击响应长度 5ms t_imp (0:1/fs:T_impulse); h_imp exp(-t_imp/tau) .* sin(2*pi*fn*t_imp); % 归一化 h_imp h_imp / max(abs(h_imp));按特征频率排布冲击序列% 生成冲击序列 t_total 2; t (0:1/fs:t_total); x_fault zeros(size(t)); % 外圈故障等间隔冲击 T_interval 1 / BPFO_Hz; imp_start 0.05; imp_times imp_start:T_interval:t_total; for ii 1:length(imp_times) idx_start round(imp_times(ii) * fs) 1; if idx_start length(h_imp) - 1 length(t) x_fault(idx_start:idx_startlength(h_imp)-1) ... x_fault(idx_start:idx_startlength(h_imp)-1) h_imp; end end内圈故障转频调制内圈故障时损伤点周期性地进出载荷区冲击幅值受转频调制。叠加方式% 内圈故障调制函数 mod_depth 0.8; % 调制深度 0~1 mod_signal 1 - mod_depth * (0.5 0.5 * cos(2*pi*fr*t)); x_inner x_fault .* mod_signal;叠加背景振动与噪声实测信号里除了故障冲击还有转子不平衡、不对中引起的转频及其谐波、随机噪声。这些都要加进去否则信号太干净跟实际差距太大。% 转频谐波不平衡/不对中成分 x_rotor 0.2 * sin(2*pi*fr*t) 0.08 * sin(4*pi*fr*t); % 白噪声通过高通滤波模拟传感器噪声特性 noise_power 0.05; x_noise noise_power * randn(size(t)); % 传感器频响特性一阶高通低通模拟加速度传感器 [b_hp, a_hp] butter(2, 10/(fs/2), high); [b_lp, a_lp] butter(4, 10000/(fs/2), low); x_noise filter(b_lp, a_lp, filter(b_hp, a_hp, x_noise)); % 合成 x_total x_fault x_rotor x_noise;随机工况参数扰动为了让生成的数据不像“复制粘贴”可以把每次仿真的特征间隔添加微小随机扰动。模拟实际转速波动% 转速波动随机游走 speed_fluctuation 0.01; fr_inst fr * (1 speed_fluctuation * randn(size(t))); phase 2*pi * cumsum(fr_inst) / fs;这样生成的数据即使来自同一批参数信号形态也有自然差异用于训练模型能有效提升泛化能力。4.3 故障仿真信号质量评价清单生成数据之后别急着拿去训练模型先按这个清单检查信号是否合理检查项方法达标标准特征频率是否准确包络谱峰值比对峰值频率与理论BPFO/BPFI偏差 1%冲击间隔是否均匀时域波形峰值间隔直方图间隔均值约等于1/特征频率方差不过大调制是否存在包络谱边频带内圈故障应有转频间隔的边带信噪比是否符合需求SNR计算根据训练需求调整噪声功率转速波动是否合理瞬时频率估计波动范围在设定范围这一步花时间磨信号质量后面所有下游工作都会受益。5. 仿真数据怎么用才有说服力验证与应用5.1 包络谱验证先证明信号对了仿真数据的第一道关是包络谱验证。把第4章生成的信号做Hilbert变换求包络再对包络做FFT看谱峰是否落在理论特征频率处。% 包络谱分析验证 env abs(hilbert(x_total)); Nfft 16384; [pxx, f] pwelch(env, hann(Nfft), Nfft/2, Nfft, fs); % 找谱峰 [pks, locs] findpeaks(10*log10(pxx), MinPeakHeight, -30); [~, sort_idx] sort(pks, descend); top5_freqs f(locs(sort_idx(1:min(5,length(sort_idx))))); fprintf(理论BPFO: %.2f Hz\n, BPFO_Hz); fprintf(包络谱前5个峰值频率: ); fprintf(%.2f , top5_freqs); fprintf(\n);如果峰值没有出现在理论值附近大概率是冲击排布有问题或调制参数设错回头检查。5.2 数据增强少量故障样本的扩充故障诊断领域最常见的问题**正常样本好几百组故障样本只有几十组。**仿真数据最大的应用价值就是扩充故障类样本。我常用的扩充策略参数微调每次生成时在特征频率±1%、幅值±20%、噪声±3dB范围内随机扰动生成上千组多测点多通道同一个故障在不同传感器位置的信号每个位置单独建模扩大通道维度混合工况改变转速区间比如500-3000rpm分段、不同负载等级生成工况覆盖更广的训练集仿真实测混合训练集里仿真数据占70%实测占30%实测占比太低模型容易偏移太高样本量不够。亲身测试过用纯仿真数据训练的CNN模型在同一个工况的实测数据上测试准确率能到70%80%如果再混入20%30%实测数据准确率能冲到90%以上。完全不用实测数据只靠仿真还没有办法完全消除域差异这是个事实。5.3 模型验证的双重保险如果你的目标是发论文或者做算法对比建议做两个层面的验证仿真到仿真同一模型用A组仿真参数训练B组仿真参数测试。这个验证的是算法在数据变化下的稳定性仿真到实测模型用了多少仿真数据训练最后一定要有实测数据做盲测。实测数据至少要有同型号轴承的正常和故障样本否则论文结论没有说服力。我之前见过一个研究用纯仿真数据训练联邦诊断模型仿真指标非常漂亮但是因为没有实测盲测审稿人直接质疑可行性。后来补了一组实测验证才通过。仿真数据是放大器不是替代品——它放大真实数据的有效信息不能凭空创造物理规律之外的信息。6. 实操中的坑与个人经验6.1 最容易忽略的几个问题采样率设置不合理。轴承固有频率通常在200010000 Hz如果采样率只有8 kHz很多冲击响应细节根本采不到包络谱也会出现混叠。工业采集卡常用12.8kHz、25.6kHz、51.2kHz。仿真时建议不低于25.6kHz否则信号质量与实测严重脱节。固有频率拍脑袋取。结构固有频率决定了冲击响应的主频取错会导致仿真信号频谱形态和实测差异很大。如果你手头有实测信号直接对故障冲击段做频谱分析取主导峰值作为fn没有实测就用经验值20005000 Hz并做多组对比别只取一个固定值。阻尼比设太大或太小。阻尼比决定冲击响应衰减快慢。zeta0.01时衰减太慢冲击会重叠包络谱模糊zeta0.1时衰减太快冲击变成孤立的窄脉冲高频成分过重。轴承钢结构的典型阻尼比在0.010.05之间建议从这个范围起步。没有转频调制硬说内圈故障。内圈故障信号必须包含转频调制成分否则包络谱上边频带消失跟外圈故障特征混淆。这相当于造了一个“四不像”数据模型学了反而有害。6.2 我踩过的一个具体坑我初期做仿真时直接把冲击响应用了同一个相位每个冲击都从正弦波的初始相位开始。结果生成信号的包络谱里出现大量与特征频率无关的精细谐波跟实测对不上。后来意识到**真实故障冲击的初始相位是由滚动体进入损伤区瞬间的接触状态决定的每次冲击相位并不相同。**每一次冲击都应该随机化正弦衰减波形的初始相位。% 修正每次冲击随机初始相位 phi_random 2*pi * rand(size(imp_times)); for ii 1:length(imp_times) idx_start round(imp_times(ii) * fs) 1; if idx_start length(t_imp) - 1 length(t) h_imp_i exp(-t_imp/tau) .* sin(2*pi*fn*t_imp phi_random(ii)); h_imp_i h_imp_i / max(abs(h_imp_i)); x_fault(idx_start:idx_startlength(t_imp)-1) ... x_fault(idx_start:idx_startlength(t_imp)-1) h_imp_i; end end这个改动之后包络谱的形态立刻正常了很多边频带分布跟实测数据的一致性也明显提高了。6.3 数据量多少合适怎么判断仿真数据并非越多越好。我的经验是先按5.1的清单验证200组信号确认质量合格后再批量生成如果验证不过生成5万组也是错的。批量生成后用PCA或t-SNE画特征分布看故障类之间是否重叠。如果不同故障类的特征分布完全重叠说明物理参数或信号参数设置有问题先排查参数而不是继续叠加样本。一个值得记住的对比CWRU公开数据集里每个工况每个故障类型的样本并不多但学术界用了十几年原因是每个样本都经过严格物理校验。数据质量永远是第一位的数量第二。扩展到更丰富的场景时可以考虑把振动仿真和电流信号仿真结合起来做电机电流特征分析MCSA能识别轴承故障通过负载转矩波动在电流信号中留下的痕迹。这条线需要把轴承动力学模型和电机电磁模型耦合MATLAB里可以用Simscape Electrical配合自定义机械负载实现是个有深度也很有意思的方向。