基于西储数据的轴承动力学建模与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)
外圈 BPFO	0.5 n (1 - d/D cosα) fr	107.4	103.4
内圈 BPFI	0.5 n (1 + d/D cosα) fr	162.2	156.2
滚动体 BSF	0.5 (D/d) (1 - (d/D cosα)²) fr	141.0	135.8
保持架 FTF	0.5 (1 - d/D cosα) fr	11.6	11.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/min，BPFO 相差约 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 之间是因为物理上滚动轴承系统的阻尼比很少超出这个范围；如果第一次迭代就顶到边界，说明共振频率还没对齐，回去调频率而不是继续迭代。