简介面向机械工程领域研究人员与技术人员的西储数据轴承动力学建模与仿真资源包围绕MATLAB/Simulink实现数据驱动建模适用于轴承性能评估、故障诊断与优化设计等场景。内容从轴承力学特性与运动规律的基础理论切入完整覆盖西储数据读取、模型参数定义、动力学方程构建以及Simulink仿真输出结果图表的全过程并附有可运行的简化MATLAB代码示例帮助入门者快速上手并根据实际工况调整参数。资源包共4个文件包含PDF与Word版完整技术文档、Markdown阅读版以及HTML网页版便于按需选择阅读方式整体仅287KB。已有610人学习这套轻量资料能在低成本条件下助力机械工程师理解数据驱动建模流程、减少实验消耗是一份兼顾理论解析与工程实践的实用参考资料。1. 先想明白西储数据驱动的轴承动力学建模与仿真到底驱动了什么西储数据驱动的轴承动力学建模与仿真并不是拿真实数据去拟合一个黑匣子而是把公开振动信号当成“参数来源”从数据里反推出轴承系统的共振频率、阻尼比、转速和故障特征频率再把这些参数写进一个明确的动力学模型用 Matlab 程序生成带可控标签的仿真信号。这个套路最典型的用途有两个一是你要验证诊断算法但公开数据样本不够、故障类型不全需要按需生成“长得很像西储数据”的仿真样本二是不满足于“随机脉冲叠加正弦波”那种玩具信号希望仿真信号在包络谱形态、共振衰减行为上和真实采集数据同口径。适合手里有 Matlab、想快速搭起一套可重复参数扫描仿真平台的工程师也适合刚接触故障诊断、想搞懂轴承故障信号是怎么形成的初学者。本章先把结论说清楚数据驱动的是参数动力学模型负责把参数变成信号。2. 从西储数据里抠出建模必需的硬参数2.1 读入 .mat 文件文件结构、采样率与转速确认西储数据的 .mat 文件里通常只有一个结构体字段字段内容是一段很长的一维振动信号。我一般不用硬编码的字段名去取数一是因为不同版本的文件字段名不完全一致二是用fieldnames动态取最省事也防止脚本换一台电脑就报错。function [sig, fs, rpm] load_bearing_dataset(matpath) data load(matpath); fld fieldnames(data); if numel(fld) ~ 1 error(该数据文件包含多个字段请手动指定); end raw data.(fld{1}); raw raw(:); % 统一成列向量 % 12k 与 48k 的区分优先按路径里的关键字判断 if contains(matpath, 48k) fs 48000; else fs 12000; end % 去掉前 0.1 秒的启动段取稳态段 sig raw(round(0.1 * fs) 1 : end); % 转速不要从文件名猜先给占位值后面按负载确认 rpm 1797; end这段代码的关键点是采样率判断。西储数据分 12 kHz 和 48 kHz 两组路径或文件名里通常带着关键字按长度粗判只在数据量差异明显时可靠万一你拿到的是一个截断样本按长度判断就会翻车。contains在旧版 Matlab 里不可用老环境换成strfind(matpath, 48k)即可。转速这里要单独强调不同负载下电机转速并不固定常见的是 1797/1772/1750/1730 r/min 这四档分别对应 0 到 3 hp 负载。不要把文件末尾的三位数字默认成转速那个数字可能是样本编号。我处理西储数据的习惯是先把负载对应的转速做成一个 switch 映射表人工确定好再用否则后面算特征频率全部偏掉。2.2 用 Wel ch 功率谱与半功率带宽估计共振频率和阻尼比故障冲击激起的是轴承系统的结构共振而不是故障特征频率本身。所以动力学模型里最重要的两个物理量是共振频率 fr 和阻尼比 ζ前者决定冲击响应振荡多快后者决定衰减多快。建模之前必须先估计这两个量否则仿真信号的时域形态和真实数据完全不像。共振频率的估计并不复杂对信号做 Wel ch 功率谱在 1 kHz 到 6 kHz 范围里找最大峰。这个峰通常对应传感器通道的共振频带。阻尼比再用半功率带宽法求找到峰值左右两侧功率降到一半处的频率点带宽除以两倍共振频率就是阻尼比。function [fr, zeta] estimate_resonance(sig, fs) [pxx, f] pwelch(sig, hann(4096), 2048, 8192, fs); pxx movmedian(pxx, 5); % 平滑防止窄带干扰冒充共振峰 band f 800 f 6000; pxx_band pxx(band); f_band f(band); [~, ipk] max(pxx_band); fr f_band(ipk); % 半功率带宽法求阻尼比 half pxx_band(ipk) / 2; left find(pxx_band(1:ipk) half, 1, last); right ipk find(pxx_band(ipk 1:end) half, 1, first); if isempty(left) || isempty(right) zeta 0.02; % 边界被截断先用经验初值 else zeta (f_band(right) - f_band(left)) / (2 * fr); end end这里有个常见的玄学点移动中位数窗口取 5 或 7 差别不大但如果窗口太大会把真实的共振峰也磨平。另一个坑是 800 Hz 以下往往有转频和低次谐波不能纳入共振峰搜索范围6 kHz 以上又可能进入高频结构模态和轴承故障诊断关心的频段关系不大。若左右半功率点找不到说明这个峰被邻峰干扰或截断了此时硬算阻尼比没有意义先给 0.02 的经验值后面用标定流程修正。2.3 故障特征频率计算公式、典型参数表与转速修正有了转速和轴承几何参数就能把各类故障的特征频率算出来。西储数据里大量使用的是某型 6205 深沟球轴承滚动体数 9 个节圆直径约 39.04 mm滚动体直径约 7.94 mm接触角为 0°。下面这份公式和频率表建议贴在你 Matlab 脚本的注释里每次对结果前先扫一眼。故障位置特征频率公式1797 r/min 下频率(Hz)1730 r/min 下频率(Hz)外圈 BPFO0.5 n (1 - d/D cosα) fr107.4103.4内圈 BPFI0.5 n (1 d/D cosα) fr162.2156.2滚动体 BSF0.5 (D/d) (1 - (d/D cosα)²) fr141.0135.8保持架 FTF0.5 (1 - d/D cosα) fr11.611.2这里的 n 是滚动体个数d 是滚动体直径D 是节圆直径α 是接触角。换成其它轴承型号时一定要替换几何参数不要拿着 6205 的参数硬套。Matlab 里我习惯写成独立函数function fc bearing_fault_frequencies(n, d, D, alpha, rpm) fr_hz rpm / 60; c cos(alpha); ratio d / D; fc [ fr_hz, 0.5 * n * (1 - ratio * c) * fr_hz, % BPFO 0.5 * n * (1 ratio * c) * fr_hz, % BPFI 0.5 * (D / d) * (1 - (ratio * c)^2) * fr_hz, % BSF 0.5 * (1 - ratio * c) * fr_hz % FTF ]; end逻辑上这组公式不需要仿真模型就能算但它是后续冲击序列的周期依据。实测中包络谱上看到的特征峰位通常和理论值差 0.5 Hz 以内如果差出 2 Hz 以上先查转速不要怀疑公式。还要知道一个事实滚动体在跑道上存在微小滑移所以特征峰附近一定有一两个小边带这是物理现象不是你程序写错了。3. 把故障冲击注入动力学方程Matlab 仿真核心实现3.1 生成故障冲击序列抖动、幅值调制与最小实现轴承故障信号的激励源是一串周期性的冲击力但冲击间隔并不严格等周期滚动体经过缺陷位置时存在随机滑移。如果完全等周期仿真包络谱会过于干净和真实数据差异极大。所以生成冲击序列时一定要加时间抖动。同时内圈故障的冲击点随内圈旋转故障点经过负载区时振动传递路径在变需要叠加转频幅值调制。function F build_impulse_sequence(t, fs, fc_fault, fr_rot, mode, A, jitter_ratio) % t : 时间向量(s) % fs : 采样率(Hz) % fc_fault : 故障特征频率(Hz) % fr_rot : 转轴频率(Hz) % mode : outer / inner / ball % A : 平均脉冲幅值 % jitter_ratio : 抖动标准差占故障周期的比例建议 0.005~0.02 T_fault 1 / fc_fault; n_pulse floor((t(end) - t(1)) / T_fault); % 每个冲击时刻 理论周期 随机抖动 t_imp t(1) 0.1 * T_fault (0:n_pulse - 1) * T_fault ... jitter_ratio * T_fault .* randn(1, n_pulse); t_imp t_imp(t_imp t(1) t_imp t(end)); A_pulse A * ones(size(t_imp)); switch mode case inner % 内圈故障故障点随机圈旋转产生转频幅值调制 A_pulse A_pulse .* (1 0.8 * cos(2 * pi * fr_rot .* t_imp)); case ball % 滚动体故障故障点随保持架旋转调制频率约 0.5 倍转频 A_pulse A_pulse .* (1 0.5 * cos(2 * pi * 0.5 * fr_rot .* t_imp)); case outer % 外圈故障固定于负载区幅值调制较弱保持常数 otherwise error(mode 只能是 outer/inner/ball); end % 将每个冲击放到最近的采样点上 F zeros(size(t)); idx round((t_imp - t(1)) * fs) 1; for i 1:numel(idx) F(idx(i)) F(idx(i)) A_pulse(i); end end抖动比例 jitter_ratio 是这里最关键的一个参数。设成 0包络谱会尖锐得像一根根针设到 0.05 以上故障特征峰可能被抹平严重时包络谱根本看不出周期。我一般从 0.01 起步再根据真实数据包络谱的谱峰宽度微调。内圈调制深度取 0.8 是经验值不要做成满调制 1.0否则信号会在某些时刻完全消失和实测不符。滚动体故障的调制频率用 0.5 倍转频近似实际中保持架转速还需要乘一个滑动因子但西储数据规模下 0.5 倍已经够用。3.2 单自由度动力学模型状态空间与 lsim 求解冲击序列是激励力把它作用到一个质量-弹簧-阻尼系统上求解二阶微分方程就得到传感器位置的振动响应。这里用 mx cx k*x F(t) 描述其中 m 等效质量k 是刚度c 是阻尼系数。等效质量取 1 kg 即可因为脉冲幅值 A 本身可以任意缩放最后再统一匹配真实信号的能量。function [acc, yt] sim_bearing_response(t, fs, F, fr_res, zeta) m 1; k (2 * pi * fr_res)^2 * m; % 由目标共振频率反推刚度 c 2 * zeta * sqrt(k * m); % 由阻尼比反推阻尼系数 Amat [0, 1; -k/m, -c/m]; Bmat [0; 1/m]; Cmat [1, 0]; sys ss(Amat, Bmat, Cmat, 0); % lsim 求解状态空间返回位移和速度 yt lsim(sys, F, t, [0; 0]); % 加速度计测的是加速度由状态方程直接求出 acc - (k * yt(:,1) c * yt(:,2)) / m F / m; % 低通消除脉冲注入带来的数值毛刺 cut min(0.95, 1.2 * fr_res / (fs / 2)); [b, a] butter(2, cut, low); acc filtfilt(b, a, acc); end这里用状态空间而不是直接叠加冲击响应好处是今后如果想在模型里加入间隙、游隙等非线性因素只要改 Amat 和 Bmat 就可以了冲击序列生成代码完全不用动。lsim返回的是系统状态第一列是位移第二列是速度加速度必须由方程x F/m - k*x/m - c*x/m计算而不是对位移做两次差分直接差分会放大脉冲时刻的数值突变。低通截止频率取共振频率的 1.2 倍既能保留共振振荡又能削掉激励导致的毛刺。3.3 加噪声与传感器通带让仿真信号和西储数据同口径裸的动力学响应太过理想工程诊断算法在无噪声仿真上表现好一上真实数据就崩这是常见翻车现场。所以后处理必须给信号加噪声并做与真实数据一致的带通滤波。我一般把信噪比控制在 5 到 15 dB具体看你要验证的算法抗噪水平。function y finalize_sim_signal(acc, fs, SNR_dB) rms_sig rms(acc); noise randn(size(acc)); rms_noise rms(noise); noise_scaled noise * rms_sig / rms_noise / 10^(SNR_dB / 20); y acc noise_scaled; % 800~6000 Hz 带通和真实信号同频段 f_lo max(0.001, 800 / (fs / 2)); f_hi min(0.99, 6000 / (fs / 2)); [b, a] butter(4, [f_lo, f_hi], bandpass); y filtfilt(b, a, y); end噪声按 RMS 比例缩放而不是按峰值缩放这样不同负载下信号强弱不同时信噪比定义保持稳定。带通上下限要和你在 2.2 节估计共振频带时用的范围一致否则后续标定会产生系统性偏差。完整调用顺序是这样的生成时间轴构造冲击序列求解动力学模型加噪滤波。仿真信号和西储真实数据后续走同一个包络谱函数才谈得上对比。fs 12000; t (0:round(4 * fs) - 1) / fs; F build_impulse_sequence(t, fs, 162.2, 29.95, inner, 1.0, 0.01); [acc, ~] sim_bearing_response(t, fs, F, 2500, 0.03); y_sim finalize_sim_signal(acc, fs, 10);这里 2500 Hz 是共振频率初值0.03 是阻尼比初值它们来自 2.2 节的估计。注意内圈故障特征频率 162.2 Hz 对应的是 1797 r/min 工况如果你要仿别的负载得把所有频率重新算一遍。4. 用西储真实数据反标定模型参数4.1 标定目标与误差函数定义仿真参数不是拍脑袋定的需要拿西储真实数据做反标定。标定目标不是让两个波形逐点吻合那不现实而是让包络谱在故障特征频率的基频和几个主要谐波处形态接近。误差函数定义如下取故障特征频率 f_c 的 1 到 4 次谐波处包络谱幅值组成向量分别归一化后计算归一化均方根误差。峰位偏差超过 1 Hz 时直接加 100 的大惩罚强制搜索往正确峰位走。function [env, f] envelope_spectrum(x, fs) env_time abs(hilbert(x)); [env, f] pwelch(env_time, hann(4096), 2048, 8192, fs); end function err envelope_error(y_sim, y_exp, fs, fc) [env_exp, f] envelope_spectrum(y_exp, fs); [env_sim, ~] envelope_spectrum(y_sim, fs); harmonics 1:4; P_exp zeros(1, 4); P_sim zeros(1, 4); for j 1:4 [~, i_exp] min(abs(f - fc * j)); [~, i_sim] min(abs(f - fc * j)); P_exp(j) env_exp(i_exp); P_sim(j) env_sim(i_sim); end % 用一阶峰做归一化消除整体幅值影响 P_exp P_exp / P_exp(1); P_sim P_sim / P_sim(1); err norm(P_sim - P_exp) / norm(P_exp); % 峰位偏差惩罚 band f (fc - 4) f (fc 4); [~, i_exp] max(env_exp(band)); [~, i_sim] max(env_sim(band)); if abs(f(band(i_exp)) - f(band(i_sim))) 1.0 err err 100; end end这个误差函数只看前 4 次谐波因为更高阶谐波幅值低受噪声影响大放进目标函数只会把搜索带偏。归一化处理后整体幅值倍率被消除标定参数就集中在共振频率、阻尼比、抖动比例和调制深度上这些才是影响包络谱形态的量。4.2 网格搜索参数范围、步长与运行量误差函数有了接下来最稳的做法是网格搜索而不是直接上 fminsearch。误差面对参数不一定是凸的直接给优化器一个初值很容易掉到局部最小值。我一般先固定抖动比例和调制深度只对共振频率 f_res 和阻尼比 ζ 做两层网格。fr0 fr_res_from_data; % 2.2 节估计值 f_res_list fr0 - 300 : 50 : fr0 300; zeta_list 0.005 : 0.005 : 0.08; err_map zeros(numel(f_res_list), numel(zeta_list)); for i 1:numel(f_res_list) for j 1:numel(zeta_list) [acc, ~] sim_bearing_response(t, fs, F, f_res_list(i), zeta_list(j)); y_sim finalize_sim_signal(acc, fs, SNR_dB); err_map(i, j) envelope_error(y_sim, y_exp, fs, fc); end end [min_err, idx_lin] min(err_map(:)); [i_best, j_best] ind2sub(size(err_map), idx_lin); best_fr f_res_list(i_best); best_zeta zeta_list(j_best);网格规模约 13×15也就是 195 次仿真每次 4 秒信号在普通台式机上大概零点几秒整个搜索几分钟内跑完。这个运行量完全可接受没必要一开始就优化算法。搜索完成后先看一眼min_err是不是远大于其它网格点如果误差面一片平坦说明你对齐的是噪声峰而不是真实特征峰得回头检查转速。4.3 验证包络谱叠画与判定标准参数标定完不是直接收工要叠画验证。把仿真信号和西储真实信号放进同一个包络谱函数同尺度画出来用眼睛确认谱峰位置和形态。[env_sim, fenv] envelope_spectrum(y_sim, fs); [env_exp, ~] envelope_spectrum(y_exp, fs); figure; semilogy(fenv, env_sim, b); hold on; semilogy(fenv, env_exp, r, LineWidth, 1.5); xlabel(频率 (Hz)); ylabel(包络谱幅值); legend(仿真信号, 西储真实数据); xlim([0, 800]);判定指标可接受范围检查位置特征峰频率偏差小于 1 Hz包络谱基频峰位一阶谐波幅值比0.8 ~ 1.2仿真/真实幅值比前 3 阶谐波衰减趋势单调衰减4.1 的 P 向量时域包络间隔变异系数5% ~ 20%抖动参数相关频率偏差超过 1 Hz 时不要继续调幅值先回去查转速谐波衰减趋势不对优先调阻尼比和调制深度。同型号轴承、不同负载下标定过的共振频率和阻尼比基本可以复用只需要重新计算转速相关的特征频率这也是这套数据驱动建模最大的工程价值。5. 建模对标过程中最容易翻车的五个坑5.1 采样率误判导致带通滤波通带错位现象仿真信号做带通滤波后能量明显比真实数据弱包络谱共振峰消失。原因load 函数按数据长度粗判采样率把 48 kHz 的数据判成了 12 kHz归一化频率全部错位。解决不要用长度判断采样率改用路径或文件名中的关键字如果文件命名不规范至少用数据说明文档确认后再写死。这个坑隐蔽的地方在于峰位可能没错但能量分布完全不对排查时容易怀疑滤波参数。5.2 转速用名义值BPFO 差出 2 到 4 Hz现象仿真和真实数据的包络谱都能看到特征峰但峰位对不齐误差函数的惩罚项一直触发。原因仿真里固定用 1797 r/min 计算特征频率而真实数据是 3 hp 负载实际转速只有 1730 r/minBPFO 相差约 3.7 Hz。解决仿真输入和真实数据必须用同一个转速建议把转速表做成语料文件或 switch 映射按负载读取。转速是整套仿真里最不该猜的参数猜错之后所有标定都会变成玄学。5.3 内圈故障忘记转频调制包络谱形态完全不对现象内圈故障的仿真包络谱在基频和谐波处非常整齐但真实数据的特征峰两侧有明显的转频边带。原因内圈故障点随内圈旋转传递路径周期性变化天然产生转频幅值调制仿真里直接用等幅脉冲序列漏掉了这个调制。解决在build_impulse_sequence里对inner模式加1 0.8 * cos(2*pi*fr_rot*t_imp)调制。这个坑我自己踩过当时白折腾了一晚上最后发现就是少了一个调制项。5.4 抖动比例设成 0仿真信号“过干净”现象仿真信号在包络谱上峰又高又窄故障分类算法跑仿真精度到 99%换到真实数据立刻掉到 70%。原因滚动体在跑道上的运动存在随机滑移冲击间隔不是严格等周期抖动设 0 等于把真实物理中的随机性完全抹掉了。解决jitter_ratio 取 0.005 到 0.02取值依据是真实数据包络谱峰宽峰太窄就加大峰被抹平就减小。这一步是决定仿真信号有没有工程价值的分水岭。5.5 把高次谐波当成共振频率现象用 Wel ch 谱找共振峰时选了一个又尖又高的峰但把它当共振频率建模后时域包络完全看不出指数衰减振荡。原因那个峰可能是故障特征频率的高次谐波或者是变频器干扰不是结构固有频率。解决对候选峰做窄带滤波后看时域包络真正的共振响应会在每个冲击后呈现明显的指数衰减正弦形态或者把多个候选峰分别代入模型选标定误差最小的那个。只看功率谱幅值挑峰是最容易踩的视觉陷阱。6. 收尾技巧用幅值比快速锁定阻尼比剩下交给系统谐振峰6.1 阻尼比的反比例修正公式网格搜索能锁定参数但每次都要跑上百次仿真参数一变就得重跑效率不高。我常用的技巧是利用单自由度系统在共振频率处的增益特性幅值增益近似等于 1/(2ζ)。也就是说包络谱故障峰的幅值和阻尼比近似成反比。若当前仿真峰幅值是 A_cur目标是 A_target那么阻尼比修正为ζ_new ζ_cur × (A_cur / A_target)用这个公式迭代三四次阻尼比就收敛到目标值附近。前提是共振频率已经对齐否则峰幅值会受到频偏影响修正不再可靠。function p envelope_peak_at(y, fs, fc) [env, f] envelope_spectrum(y, fs); band f fc - 4 f fc 4; p max(env(band)); end % 标定阻尼比的快速迭代 zeta_cur 0.03; A_target envelope_peak_at(y_exp, fs, fc); % 真实数据目标峰幅值 for iter 1:4 [acc, ~] sim_bearing_response(t, fs, F, best_fr, zeta_cur); y_sim finalize_sim_signal(acc, fs, SNR_dB); A_cur envelope_peak_at(y_sim, fs, fc); zeta_cur zeta_cur * (A_cur / A_target); zeta_cur max(0.005, min(zeta_cur, 0.1)); % 限幅防震荡 if abs(zeta_cur / zeta_old - 1) 0.01 break; end zeta_old zeta_cur; end注意这里有三个细节A_target 要从真实包络谱里取但取之前要先把噪声底去掉否则噪声会把幅值抬高导致阻尼比越迭代越偏。噪声底可以用特征峰两侧 5 Hz 范围内包络谱的中位数近似。限幅到 0.005 到 0.1 之间是因为物理上滚动轴承系统的阻尼比很少超出这个范围如果第一次迭代就顶到边界说明共振频率还没对齐回去调频率而不是继续迭代。6.2 锁峰位在前调阻尼在后调参顺序永远固定先锁转速再锁共振频率最后调阻尼比抖动比例放最后微调。共振频率的直接依据是包络谱主峰位置把仿真结果移动到你目标峰位的 1 Hz 带宽内阻尼比的反比例修正才有效。抖动比例对基峰幅值影响不大但对谱峰宽度影响明显等幅值收敛了再动它。我自己的习惯是拿到一个新样本第一件事永远不是看波形而是先确认转速、画出共振频带。那次内圈故障忘掉调制白折腾一晚上之后我把这条经验刻进了工作流。希望这套方法也能帮你少走几步弯路。本文还有配套的精品资源点击获取 
长河君 · 2026/10/16 20:37:54 · 阅读 1.2k 基于西储数据的轴承动力学建模与Matlab仿真实现
简介面向机械工程领域研究人员与技术人员的西储数据轴承动力学建模与仿真资源包围绕MATLAB/Simulink实现数据驱动建模适用于轴承性能评估、故障诊断与优化设计等场景。内容从轴承力学特性与运动规律的基础理论切入完整覆盖西储数据读取、模型参数定义、动力学方程构建以及Simulink仿真输出结果图表的全过程并附有可运行的简化MATLAB代码示例帮助入门者快速上手并根据实际工况调整参数。资源包共4个文件包含PDF与Word版完整技术文档、Markdown阅读版以及HTML网页版便于按需选择阅读方式整体仅287KB。已有610人学习这套轻量资料能在低成本条件下助力机械工程师理解数据驱动建模流程、减少实验消耗是一份兼顾理论解析与工程实践的实用参考资料。1. 先想明白西储数据驱动的轴承动力学建模与仿真到底驱动了什么西储数据驱动的轴承动力学建模与仿真并不是拿真实数据去拟合一个黑匣子而是把公开振动信号当成“参数来源”从数据里反推出轴承系统的共振频率、阻尼比、转速和故障特征频率再把这些参数写进一个明确的动力学模型用 Matlab 程序生成带可控标签的仿真信号。这个套路最典型的用途有两个一是你要验证诊断算法但公开数据样本不够、故障类型不全需要按需生成“长得很像西储数据”的仿真样本二是不满足于“随机脉冲叠加正弦波”那种玩具信号希望仿真信号在包络谱形态、共振衰减行为上和真实采集数据同口径。适合手里有 Matlab、想快速搭起一套可重复参数扫描仿真平台的工程师也适合刚接触故障诊断、想搞懂轴承故障信号是怎么形成的初学者。本章先把结论说清楚数据驱动的是参数动力学模型负责把参数变成信号。2. 从西储数据里抠出建模必需的硬参数2.1 读入 .mat 文件文件结构、采样率与转速确认西储数据的 .mat 文件里通常只有一个结构体字段字段内容是一段很长的一维振动信号。我一般不用硬编码的字段名去取数一是因为不同版本的文件字段名不完全一致二是用fieldnames动态取最省事也防止脚本换一台电脑就报错。function [sig, fs, rpm] load_bearing_dataset(matpath) data load(matpath); fld fieldnames(data); if numel(fld) ~ 1 error(该数据文件包含多个字段请手动指定); end raw data.(fld{1}); raw raw(:); % 统一成列向量 % 12k 与 48k 的区分优先按路径里的关键字判断 if contains(matpath, 48k) fs 48000; else fs 12000; end % 去掉前 0.1 秒的启动段取稳态段 sig raw(round(0.1 * fs) 1 : end); % 转速不要从文件名猜先给占位值后面按负载确认 rpm 1797; end这段代码的关键点是采样率判断。西储数据分 12 kHz 和 48 kHz 两组路径或文件名里通常带着关键字按长度粗判只在数据量差异明显时可靠万一你拿到的是一个截断样本按长度判断就会翻车。contains在旧版 Matlab 里不可用老环境换成strfind(matpath, 48k)即可。转速这里要单独强调不同负载下电机转速并不固定常见的是 1797/1772/1750/1730 r/min 这四档分别对应 0 到 3 hp 负载。不要把文件末尾的三位数字默认成转速那个数字可能是样本编号。我处理西储数据的习惯是先把负载对应的转速做成一个 switch 映射表人工确定好再用否则后面算特征频率全部偏掉。2.2 用 Wel ch 功率谱与半功率带宽估计共振频率和阻尼比故障冲击激起的是轴承系统的结构共振而不是故障特征频率本身。所以动力学模型里最重要的两个物理量是共振频率 fr 和阻尼比 ζ前者决定冲击响应振荡多快后者决定衰减多快。建模之前必须先估计这两个量否则仿真信号的时域形态和真实数据完全不像。共振频率的估计并不复杂对信号做 Wel ch 功率谱在 1 kHz 到 6 kHz 范围里找最大峰。这个峰通常对应传感器通道的共振频带。阻尼比再用半功率带宽法求找到峰值左右两侧功率降到一半处的频率点带宽除以两倍共振频率就是阻尼比。function [fr, zeta] estimate_resonance(sig, fs) [pxx, f] pwelch(sig, hann(4096), 2048, 8192, fs); pxx movmedian(pxx, 5); % 平滑防止窄带干扰冒充共振峰 band f 800 f 6000; pxx_band pxx(band); f_band f(band); [~, ipk] max(pxx_band); fr f_band(ipk); % 半功率带宽法求阻尼比 half pxx_band(ipk) / 2; left find(pxx_band(1:ipk) half, 1, last); right ipk find(pxx_band(ipk 1:end) half, 1, first); if isempty(left) || isempty(right) zeta 0.02; % 边界被截断先用经验初值 else zeta (f_band(right) - f_band(left)) / (2 * fr); end end这里有个常见的玄学点移动中位数窗口取 5 或 7 差别不大但如果窗口太大会把真实的共振峰也磨平。另一个坑是 800 Hz 以下往往有转频和低次谐波不能纳入共振峰搜索范围6 kHz 以上又可能进入高频结构模态和轴承故障诊断关心的频段关系不大。若左右半功率点找不到说明这个峰被邻峰干扰或截断了此时硬算阻尼比没有意义先给 0.02 的经验值后面用标定流程修正。2.3 故障特征频率计算公式、典型参数表与转速修正有了转速和轴承几何参数就能把各类故障的特征频率算出来。西储数据里大量使用的是某型 6205 深沟球轴承滚动体数 9 个节圆直径约 39.04 mm滚动体直径约 7.94 mm接触角为 0°。下面这份公式和频率表建议贴在你 Matlab 脚本的注释里每次对结果前先扫一眼。故障位置特征频率公式1797 r/min 下频率(Hz)1730 r/min 下频率(Hz)外圈 BPFO0.5 n (1 - d/D cosα) fr107.4103.4内圈 BPFI0.5 n (1 d/D cosα) fr162.2156.2滚动体 BSF0.5 (D/d) (1 - (d/D cosα)²) fr141.0135.8保持架 FTF0.5 (1 - d/D cosα) fr11.611.2这里的 n 是滚动体个数d 是滚动体直径D 是节圆直径α 是接触角。换成其它轴承型号时一定要替换几何参数不要拿着 6205 的参数硬套。Matlab 里我习惯写成独立函数function fc bearing_fault_frequencies(n, d, D, alpha, rpm) fr_hz rpm / 60; c cos(alpha); ratio d / D; fc [ fr_hz, 0.5 * n * (1 - ratio * c) * fr_hz, % BPFO 0.5 * n * (1 ratio * c) * fr_hz, % BPFI 0.5 * (D / d) * (1 - (ratio * c)^2) * fr_hz, % BSF 0.5 * (1 - ratio * c) * fr_hz % FTF ]; end逻辑上这组公式不需要仿真模型就能算但它是后续冲击序列的周期依据。实测中包络谱上看到的特征峰位通常和理论值差 0.5 Hz 以内如果差出 2 Hz 以上先查转速不要怀疑公式。还要知道一个事实滚动体在跑道上存在微小滑移所以特征峰附近一定有一两个小边带这是物理现象不是你程序写错了。3. 把故障冲击注入动力学方程Matlab 仿真核心实现3.1 生成故障冲击序列抖动、幅值调制与最小实现轴承故障信号的激励源是一串周期性的冲击力但冲击间隔并不严格等周期滚动体经过缺陷位置时存在随机滑移。如果完全等周期仿真包络谱会过于干净和真实数据差异极大。所以生成冲击序列时一定要加时间抖动。同时内圈故障的冲击点随内圈旋转故障点经过负载区时振动传递路径在变需要叠加转频幅值调制。function F build_impulse_sequence(t, fs, fc_fault, fr_rot, mode, A, jitter_ratio) % t : 时间向量(s) % fs : 采样率(Hz) % fc_fault : 故障特征频率(Hz) % fr_rot : 转轴频率(Hz) % mode : outer / inner / ball % A : 平均脉冲幅值 % jitter_ratio : 抖动标准差占故障周期的比例建议 0.005~0.02 T_fault 1 / fc_fault; n_pulse floor((t(end) - t(1)) / T_fault); % 每个冲击时刻 理论周期 随机抖动 t_imp t(1) 0.1 * T_fault (0:n_pulse - 1) * T_fault ... jitter_ratio * T_fault .* randn(1, n_pulse); t_imp t_imp(t_imp t(1) t_imp t(end)); A_pulse A * ones(size(t_imp)); switch mode case inner % 内圈故障故障点随机圈旋转产生转频幅值调制 A_pulse A_pulse .* (1 0.8 * cos(2 * pi * fr_rot .* t_imp)); case ball % 滚动体故障故障点随保持架旋转调制频率约 0.5 倍转频 A_pulse A_pulse .* (1 0.5 * cos(2 * pi * 0.5 * fr_rot .* t_imp)); case outer % 外圈故障固定于负载区幅值调制较弱保持常数 otherwise error(mode 只能是 outer/inner/ball); end % 将每个冲击放到最近的采样点上 F zeros(size(t)); idx round((t_imp - t(1)) * fs) 1; for i 1:numel(idx) F(idx(i)) F(idx(i)) A_pulse(i); end end抖动比例 jitter_ratio 是这里最关键的一个参数。设成 0包络谱会尖锐得像一根根针设到 0.05 以上故障特征峰可能被抹平严重时包络谱根本看不出周期。我一般从 0.01 起步再根据真实数据包络谱的谱峰宽度微调。内圈调制深度取 0.8 是经验值不要做成满调制 1.0否则信号会在某些时刻完全消失和实测不符。滚动体故障的调制频率用 0.5 倍转频近似实际中保持架转速还需要乘一个滑动因子但西储数据规模下 0.5 倍已经够用。3.2 单自由度动力学模型状态空间与 lsim 求解冲击序列是激励力把它作用到一个质量-弹簧-阻尼系统上求解二阶微分方程就得到传感器位置的振动响应。这里用 mx cx k*x F(t) 描述其中 m 等效质量k 是刚度c 是阻尼系数。等效质量取 1 kg 即可因为脉冲幅值 A 本身可以任意缩放最后再统一匹配真实信号的能量。function [acc, yt] sim_bearing_response(t, fs, F, fr_res, zeta) m 1; k (2 * pi * fr_res)^2 * m; % 由目标共振频率反推刚度 c 2 * zeta * sqrt(k * m); % 由阻尼比反推阻尼系数 Amat [0, 1; -k/m, -c/m]; Bmat [0; 1/m]; Cmat [1, 0]; sys ss(Amat, Bmat, Cmat, 0); % lsim 求解状态空间返回位移和速度 yt lsim(sys, F, t, [0; 0]); % 加速度计测的是加速度由状态方程直接求出 acc - (k * yt(:,1) c * yt(:,2)) / m F / m; % 低通消除脉冲注入带来的数值毛刺 cut min(0.95, 1.2 * fr_res / (fs / 2)); [b, a] butter(2, cut, low); acc filtfilt(b, a, acc); end这里用状态空间而不是直接叠加冲击响应好处是今后如果想在模型里加入间隙、游隙等非线性因素只要改 Amat 和 Bmat 就可以了冲击序列生成代码完全不用动。lsim返回的是系统状态第一列是位移第二列是速度加速度必须由方程x F/m - k*x/m - c*x/m计算而不是对位移做两次差分直接差分会放大脉冲时刻的数值突变。低通截止频率取共振频率的 1.2 倍既能保留共振振荡又能削掉激励导致的毛刺。3.3 加噪声与传感器通带让仿真信号和西储数据同口径裸的动力学响应太过理想工程诊断算法在无噪声仿真上表现好一上真实数据就崩这是常见翻车现场。所以后处理必须给信号加噪声并做与真实数据一致的带通滤波。我一般把信噪比控制在 5 到 15 dB具体看你要验证的算法抗噪水平。function y finalize_sim_signal(acc, fs, SNR_dB) rms_sig rms(acc); noise randn(size(acc)); rms_noise rms(noise); noise_scaled noise * rms_sig / rms_noise / 10^(SNR_dB / 20); y acc noise_scaled; % 800~6000 Hz 带通和真实信号同频段 f_lo max(0.001, 800 / (fs / 2)); f_hi min(0.99, 6000 / (fs / 2)); [b, a] butter(4, [f_lo, f_hi], bandpass); y filtfilt(b, a, y); end噪声按 RMS 比例缩放而不是按峰值缩放这样不同负载下信号强弱不同时信噪比定义保持稳定。带通上下限要和你在 2.2 节估计共振频带时用的范围一致否则后续标定会产生系统性偏差。完整调用顺序是这样的生成时间轴构造冲击序列求解动力学模型加噪滤波。仿真信号和西储真实数据后续走同一个包络谱函数才谈得上对比。fs 12000; t (0:round(4 * fs) - 1) / fs; F build_impulse_sequence(t, fs, 162.2, 29.95, inner, 1.0, 0.01); [acc, ~] sim_bearing_response(t, fs, F, 2500, 0.03); y_sim finalize_sim_signal(acc, fs, 10);这里 2500 Hz 是共振频率初值0.03 是阻尼比初值它们来自 2.2 节的估计。注意内圈故障特征频率 162.2 Hz 对应的是 1797 r/min 工况如果你要仿别的负载得把所有频率重新算一遍。4. 用西储真实数据反标定模型参数4.1 标定目标与误差函数定义仿真参数不是拍脑袋定的需要拿西储真实数据做反标定。标定目标不是让两个波形逐点吻合那不现实而是让包络谱在故障特征频率的基频和几个主要谐波处形态接近。误差函数定义如下取故障特征频率 f_c 的 1 到 4 次谐波处包络谱幅值组成向量分别归一化后计算归一化均方根误差。峰位偏差超过 1 Hz 时直接加 100 的大惩罚强制搜索往正确峰位走。function [env, f] envelope_spectrum(x, fs) env_time abs(hilbert(x)); [env, f] pwelch(env_time, hann(4096), 2048, 8192, fs); end function err envelope_error(y_sim, y_exp, fs, fc) [env_exp, f] envelope_spectrum(y_exp, fs); [env_sim, ~] envelope_spectrum(y_sim, fs); harmonics 1:4; P_exp zeros(1, 4); P_sim zeros(1, 4); for j 1:4 [~, i_exp] min(abs(f - fc * j)); [~, i_sim] min(abs(f - fc * j)); P_exp(j) env_exp(i_exp); P_sim(j) env_sim(i_sim); end % 用一阶峰做归一化消除整体幅值影响 P_exp P_exp / P_exp(1); P_sim P_sim / P_sim(1); err norm(P_sim - P_exp) / norm(P_exp); % 峰位偏差惩罚 band f (fc - 4) f (fc 4); [~, i_exp] max(env_exp(band)); [~, i_sim] max(env_sim(band)); if abs(f(band(i_exp)) - f(band(i_sim))) 1.0 err err 100; end end这个误差函数只看前 4 次谐波因为更高阶谐波幅值低受噪声影响大放进目标函数只会把搜索带偏。归一化处理后整体幅值倍率被消除标定参数就集中在共振频率、阻尼比、抖动比例和调制深度上这些才是影响包络谱形态的量。4.2 网格搜索参数范围、步长与运行量误差函数有了接下来最稳的做法是网格搜索而不是直接上 fminsearch。误差面对参数不一定是凸的直接给优化器一个初值很容易掉到局部最小值。我一般先固定抖动比例和调制深度只对共振频率 f_res 和阻尼比 ζ 做两层网格。fr0 fr_res_from_data; % 2.2 节估计值 f_res_list fr0 - 300 : 50 : fr0 300; zeta_list 0.005 : 0.005 : 0.08; err_map zeros(numel(f_res_list), numel(zeta_list)); for i 1:numel(f_res_list) for j 1:numel(zeta_list) [acc, ~] sim_bearing_response(t, fs, F, f_res_list(i), zeta_list(j)); y_sim finalize_sim_signal(acc, fs, SNR_dB); err_map(i, j) envelope_error(y_sim, y_exp, fs, fc); end end [min_err, idx_lin] min(err_map(:)); [i_best, j_best] ind2sub(size(err_map), idx_lin); best_fr f_res_list(i_best); best_zeta zeta_list(j_best);网格规模约 13×15也就是 195 次仿真每次 4 秒信号在普通台式机上大概零点几秒整个搜索几分钟内跑完。这个运行量完全可接受没必要一开始就优化算法。搜索完成后先看一眼min_err是不是远大于其它网格点如果误差面一片平坦说明你对齐的是噪声峰而不是真实特征峰得回头检查转速。4.3 验证包络谱叠画与判定标准参数标定完不是直接收工要叠画验证。把仿真信号和西储真实信号放进同一个包络谱函数同尺度画出来用眼睛确认谱峰位置和形态。[env_sim, fenv] envelope_spectrum(y_sim, fs); [env_exp, ~] envelope_spectrum(y_exp, fs); figure; semilogy(fenv, env_sim, b); hold on; semilogy(fenv, env_exp, r, LineWidth, 1.5); xlabel(频率 (Hz)); ylabel(包络谱幅值); legend(仿真信号, 西储真实数据); xlim([0, 800]);判定指标可接受范围检查位置特征峰频率偏差小于 1 Hz包络谱基频峰位一阶谐波幅值比0.8 ~ 1.2仿真/真实幅值比前 3 阶谐波衰减趋势单调衰减4.1 的 P 向量时域包络间隔变异系数5% ~ 20%抖动参数相关频率偏差超过 1 Hz 时不要继续调幅值先回去查转速谐波衰减趋势不对优先调阻尼比和调制深度。同型号轴承、不同负载下标定过的共振频率和阻尼比基本可以复用只需要重新计算转速相关的特征频率这也是这套数据驱动建模最大的工程价值。5. 建模对标过程中最容易翻车的五个坑5.1 采样率误判导致带通滤波通带错位现象仿真信号做带通滤波后能量明显比真实数据弱包络谱共振峰消失。原因load 函数按数据长度粗判采样率把 48 kHz 的数据判成了 12 kHz归一化频率全部错位。解决不要用长度判断采样率改用路径或文件名中的关键字如果文件命名不规范至少用数据说明文档确认后再写死。这个坑隐蔽的地方在于峰位可能没错但能量分布完全不对排查时容易怀疑滤波参数。5.2 转速用名义值BPFO 差出 2 到 4 Hz现象仿真和真实数据的包络谱都能看到特征峰但峰位对不齐误差函数的惩罚项一直触发。原因仿真里固定用 1797 r/min 计算特征频率而真实数据是 3 hp 负载实际转速只有 1730 r/minBPFO 相差约 3.7 Hz。解决仿真输入和真实数据必须用同一个转速建议把转速表做成语料文件或 switch 映射按负载读取。转速是整套仿真里最不该猜的参数猜错之后所有标定都会变成玄学。5.3 内圈故障忘记转频调制包络谱形态完全不对现象内圈故障的仿真包络谱在基频和谐波处非常整齐但真实数据的特征峰两侧有明显的转频边带。原因内圈故障点随内圈旋转传递路径周期性变化天然产生转频幅值调制仿真里直接用等幅脉冲序列漏掉了这个调制。解决在build_impulse_sequence里对inner模式加1 0.8 * cos(2*pi*fr_rot*t_imp)调制。这个坑我自己踩过当时白折腾了一晚上最后发现就是少了一个调制项。5.4 抖动比例设成 0仿真信号“过干净”现象仿真信号在包络谱上峰又高又窄故障分类算法跑仿真精度到 99%换到真实数据立刻掉到 70%。原因滚动体在跑道上的运动存在随机滑移冲击间隔不是严格等周期抖动设 0 等于把真实物理中的随机性完全抹掉了。解决jitter_ratio 取 0.005 到 0.02取值依据是真实数据包络谱峰宽峰太窄就加大峰被抹平就减小。这一步是决定仿真信号有没有工程价值的分水岭。5.5 把高次谐波当成共振频率现象用 Wel ch 谱找共振峰时选了一个又尖又高的峰但把它当共振频率建模后时域包络完全看不出指数衰减振荡。原因那个峰可能是故障特征频率的高次谐波或者是变频器干扰不是结构固有频率。解决对候选峰做窄带滤波后看时域包络真正的共振响应会在每个冲击后呈现明显的指数衰减正弦形态或者把多个候选峰分别代入模型选标定误差最小的那个。只看功率谱幅值挑峰是最容易踩的视觉陷阱。6. 收尾技巧用幅值比快速锁定阻尼比剩下交给系统谐振峰6.1 阻尼比的反比例修正公式网格搜索能锁定参数但每次都要跑上百次仿真参数一变就得重跑效率不高。我常用的技巧是利用单自由度系统在共振频率处的增益特性幅值增益近似等于 1/(2ζ)。也就是说包络谱故障峰的幅值和阻尼比近似成反比。若当前仿真峰幅值是 A_cur目标是 A_target那么阻尼比修正为ζ_new ζ_cur × (A_cur / A_target)用这个公式迭代三四次阻尼比就收敛到目标值附近。前提是共振频率已经对齐否则峰幅值会受到频偏影响修正不再可靠。function p envelope_peak_at(y, fs, fc) [env, f] envelope_spectrum(y, fs); band f fc - 4 f fc 4; p max(env(band)); end % 标定阻尼比的快速迭代 zeta_cur 0.03; A_target envelope_peak_at(y_exp, fs, fc); % 真实数据目标峰幅值 for iter 1:4 [acc, ~] sim_bearing_response(t, fs, F, best_fr, zeta_cur); y_sim finalize_sim_signal(acc, fs, SNR_dB); A_cur envelope_peak_at(y_sim, fs, fc); zeta_cur zeta_cur * (A_cur / A_target); zeta_cur max(0.005, min(zeta_cur, 0.1)); % 限幅防震荡 if abs(zeta_cur / zeta_old - 1) 0.01 break; end zeta_old zeta_cur; end注意这里有三个细节A_target 要从真实包络谱里取但取之前要先把噪声底去掉否则噪声会把幅值抬高导致阻尼比越迭代越偏。噪声底可以用特征峰两侧 5 Hz 范围内包络谱的中位数近似。限幅到 0.005 到 0.1 之间是因为物理上滚动轴承系统的阻尼比很少超出这个范围如果第一次迭代就顶到边界说明共振频率还没对齐回去调频率而不是继续迭代。6.2 锁峰位在前调阻尼在后调参顺序永远固定先锁转速再锁共振频率最后调阻尼比抖动比例放最后微调。共振频率的直接依据是包络谱主峰位置把仿真结果移动到你目标峰位的 1 Hz 带宽内阻尼比的反比例修正才有效。抖动比例对基峰幅值影响不大但对谱峰宽度影响明显等幅值收敛了再动它。我自己的习惯是拿到一个新样本第一件事永远不是看波形而是先确认转速、画出共振频带。那次内圈故障忘掉调制白折腾一晚上之后我把这条经验刻进了工作流。希望这套方法也能帮你少走几步弯路。本文还有配套的精品资源点击获取