拓冰网站定制 
轴承故障诊断这个方向我陆陆续续做了快八年。从最早拿着测振笔在现场一遍遍测数据到后来用MATLAB做动力学模拟反推故障机理这中间最大的感悟是现场实测信号太“脏”了你很难分清哪些特征来自轴承本身哪些是齿轮啮合、工频干扰、结构共振混进来的。这也就是为什么一定要做故障机理建模和数值计算——你先把“干净”的故障信号在仿真里搞明白再去认现场信号就有底气了。这篇内容围绕滚动轴承动力学模拟展开覆盖故障机理建模、数值计算、多故障类型模拟以及最终数据分析报告怎么出。适合正在做故障诊断相关课题的学生、刚接触设备状态监测的工程师或者想用MATLAB把轴承振动信号摸透的入门研究者。我会把自己踩过的坑、验证过的参数、调试思路都写在里面照着做你至少能少走三个月弯路。1. 内容整体设计与思路拆解1.1 为什么选MATLAB做轴承动力学模拟做轴承故障模拟市面上其实有几条路。一是纯数学推导用微分方程描述滚动体与滚道之间的接触力然后自己写数值积分求解这活儿用C、Python都能干但开发周期长二是用商业有限元软件做瞬态动力学分析精度高但建模复杂一个模型折腾一两周很正常三是用MATLAB/Simulink既有成熟的ODE求解器又有强大的信号分析工具箱一套流程下来不用来回切换工具。我自己用过Python试过类似的工作NumPySciPy组合确实灵活但到了信号处理环节MATLAB的Signal Processing Toolbox里现成的函数太方便了尤其是后面做包络谱分析、时频分析的时候一行命令能顶Python十几行代码。还有一个现实因素是很多高校和企业都在用MATLAB你写出来的东西别人能直接跑沟通成本低。就我个人的实践来说MATLAB做轴承动力学模拟的核心优势在于“算法验证”环节。模拟信号生成只是第一步后面要验证包络解调、谱峭度、小波变换这些诊断算法靠不靠谱MATLAB里调试这些特别顺手。后面我会专门讲怎么用这套模拟信号去给算法“打分”。1.2 整体建模方案与数值计算路线完整的轴承动力学模拟流程大致分四步几何与运动学参数定义、故障激励建模、系统振动响应求解、信号分析与故障特征提取。第一步是把自己要模拟的轴承型号参数吃透节圆直径、滚动体直径、滚动体数量、接触角这些直接决定故障特征频率。第二步是建模故障激励滚动轴承的故障本质上是滚动体经过损伤部位时产生周期性冲击这个冲击的幅值、频率、衰减特性是核心。第三步是把冲击激励加到系统模型里求解振动响应。第四步是信号处理把仿真得到的时域信号做FFT、包络谱验证故障特征频率是否和理论值吻合。数值计算路线上我对新手最推荐的方式是“先离散后求解”。就是把振动微分方程用数值积分方法逐步推进MATLAB的ode45是自适应步长的龙格库塔法在处理这类含间歇性冲击的系统时稳定性不错。但要注意冲击信号作用时间极短若步长控制不当很容易漏掉冲击峰值后面我会详细讲怎么处理这个坑。1.3 为什么先建立“简化模型”而不是全三维有限元模型很多初学者上来就想做全三维的有限元模型轴承所有零件全部建出来滚动体、保持架、内外圈一个不落然后跑瞬态动力学分析。坦白说这条路对纯学术研究有价值但对于绝大多数工程诊断需求是过度设计。原因很简单你建模的目的是搞清楚“故障会在振动信号里产生什么特征”而不是“轴承内部应力分布到底有多复杂”。在故障诊断这个场景里我们关心的是振动信号中周期性冲击的出现规律、调制特征、频率分布而这些主要由轴承的运动学关系和接触冲击的周期性决定。所以用集中参数模型或者等效激励模型就足够了。我用过Matlab里自己写简化的“外圈固定、内圈旋转”模型自由度不需要太多本质上是把轴承系统简化为弹簧-质量-阻尼系统滚动体通过故障区域时产生冲击力这个冲击力作用到系统上激起固有频率的瞬态振动。这种模型有明确的物理意义又不需要巨量的计算资源一套参数跑几千转的模拟数据也就是几分钟的事。先把简化模型吃透再去拓展更复杂的模型这个顺序不会错。2. 滚动轴承故障机理与特征频率计算2.1 滚动轴承的几何结构与故障类型划分滚动轴承的基本结构说白了就是内外两个圈中间一堆滚动体滚动体之间还有个保持架防止它们挤在一起。运转的时候内圈随轴旋转滚动体在内外滚道之间做滚动兼滑动运动把轴承载荷传递给外圈和轴承座。在故障诊断里我们主要关注四类故障外圈故障、内圈故障、滚动体故障、保持架故障。这四类故障在振动信号里产生的特征完全不同。外圈故障的信号最好认因为外圈通常固定在轴承座上故障点位置相对传感器不变所以脉冲序列是均匀的内圈故障就复杂了因为故障点跟随转轴旋转冲击力传到传感器时要经过滚动体和外圈路径一直在变信号会出现明显的幅值调制滚动体故障则更隐蔽滚动体自转的同时还在公转故障点会周期性接触内外滚道接触的位置和载荷不断变化信号特征最不容易提取。2.2 故障特征频率的理论计算从运动学关系说起要理解故障特征频率就得先知道滚动体怎么运动。以深沟球轴承为例设节圆直径为 (D)滚动体直径为 (d)滚动体数量为 (Z)接触角为 (\alpha)转频为 (f_r)单位Hz。把轴承当作纯滚动接触来分析内外圈与滚动体之间没有相对滑动就能推出四类故障的特征频率公式。外圈故障特征频率 (BPFO) 的计算公式是[ BPFO \frac{Z}{2} \cdot f_r \cdot \left(1 - \frac{d}{D}\cos\alpha\right) ]内圈故障特征频率 (BPFI) 是[ BPFI \frac{Z}{2} \cdot f_r \cdot \left(1 \frac{d}{D}\cos\alpha\right) ]滚动体故障特征频率 (BSF) 是[ BSF \frac{D}{2d} \cdot f_r \cdot \left(1 - \left(\frac{d}{D}\cos\alpha\right)^2\right) ]保持架故障特征频率 (FTF) 是[ FTF \frac{1}{2} \cdot f_r \cdot \left(1 - \frac{d}{D}\cos\alpha\right) ]注意(BSF) 公式里有时会乘个系数2这是因为滚动体每转一圈会和内外滚道各接触一次但工程上通常取上面这个值作为故障冲击的主频率。实际算的时候这些公式算出来的频率只是理论值真实轴承会因为载荷、游隙、润滑状态等因素产生微小的滑差频率有1%到2%的波动是正常现象。2.3 参数敏感性分析哪些参数对特征频率影响最大我刚刚学这个的时候对着公式看了半天只觉得每个参数都有影响但不知道哪个影响大。后来做模拟做多了慢慢摸出了一些门道。滚动体数量 (Z) 的影响是线性的也是最直接的。其他参数不变的情况下(Z) 越多特征频率越高频谱上的谱线间距也越大故障识别越容易。节圆直径 (D) 和滚动体直径 (d) 的比值也很关键它决定了内外圈特征频率的差异程度。接触角 (\alpha) 在小角度范围内对结果影响不大但到了角接触轴承那种大接触角的情况就要重新考虑了。这里我一般会在MATLAB里写一个参数敏感性分析脚本逐个变量小幅扰动看看特征频率变化率。做这个不是为了写论文凑图而是为了搞清楚如果现场给的轴承参数不太准算出来的特征频率误差到底有多大诊断结论靠不靠谱。这个工作做完你对公式的理解会深很多。2.4 用MATLAB实现特征频率的自动计算既然参数和公式都明确了写一个特征频率计算函数就水到渠成。我自己常用的是这个版本代码不多但足够实用function [BPFO, BPFI, BSF, FTF] bearing_fault_frequencies(D, d, Z, alpha_deg, fr) % 计算滚动轴承四类故障特征频率 % 输入 % D - 节圆直径 (mm) % d - 滚动体直径 (mm) % Z - 滚动体数量 % alpha_deg - 接触角 (度) % fr - 转频 (Hz) % 输出四类故障特征频率 alpha deg2rad(alpha_deg); cos_a cos(alpha); term1 1 - (d / D) * cos_a; term2 1 (d / D) * cos_a; BPFO (Z / 2) * fr * term1; BPFI (Z / 2) * fr * term2; BSF (D / (2 * d)) * fr * (1 - ((d / D) * cos_a)^2); FTF (1 / 2) * fr * term1; end调用的时候比如模拟6205深沟球轴承节圆直径约39.04mm滚动体直径约7.94mm滚动体数量9个接触角取0度转频取29.5Hz约1770rpm算出来BPFO约为106HzBPFI约为159.5Hz。后面模拟信号的时候这些频率就是我要在频谱和包络谱里寻找的目标。注意不同轴承型号的几何参数差异很大同一型号不同厂家的加工误差也不一样。有条件的话优先查轴承手册拿不到精确参数就用近似值但一定得在分析报告里注明。3. 多故障类型的动力学建模与数值计算3.1 冲击激励的数学描述周期脉冲与衰减振荡轴承局部故障的振动响应在机理上可以分解为两个过程一个是故障点经过接触区时产生的冲击力另一个是这个冲击力激起的系统固有振动。先说冲击力。一个理想的局部缺陷比如外滚道上的一个剥落坑滚动体每次滚过这个坑的时候都会产生一个极短时长的冲击脉冲。这个脉冲在频域里能量分布很宽能够激起轴承系统从几百赫兹到几千赫兹的固有频率。工程上常用指数衰减正弦波来描述系统的冲击响应也就是[ x_{i}(t) A \cdot e^{-\zeta \omega_n (t - t_i)} \cdot \sin\left(\omega_d (t - t_i)\right), \quad t \ge t_i ]其中 (\omega_n 2\pi f_n) 是系统固有角频率(\omega_d \omega_n\sqrt{1-\zeta^2}) 是有阻尼固有角频率(\zeta) 是阻尼比(f_n) 是系统固有频率(A) 是冲击幅值(t_i) 是第 (i) 次冲击发生的时刻。这段描述很重要它把轴承故障信号的长相给定下来了一串衰减振荡周期性出现周期等于特征频率的倒数。比如外圈故障BPFO是106Hz那么每1/106秒大概0.00943秒就出现一次冲击。3.2 外圈故障模拟实现最“规矩”的故障信号外圈故障的模拟在所有故障里最简单、也最适合入门。原因正如前面说的外圈固定在轴承座上故障点位置相对传感器恒定每次冲击的传递路径基本不变所以冲击幅值均匀、间隔均匀。在MATLAB里我一般这么构造信号fs 25600; % 采样率 25.6kHz T 2; % 信号时长 2秒 t 0:1/fs:T-1/fs; t_total length(t); x zeros(1, t_total); fr 29.5; % 转频 BPFO 106; % 外圈特征频率 T_fault 1/BPFO; % 故障冲击周期 fn 2000; % 系统固有频率 2kHz zeta 0.05; % 阻尼比 A 1.0; % 冲击幅值 % 生成冲击序列 for k 1:floor(T/T_fault) t_k (k-1) * T_fault; idx_k round(t_k * fs) 1; if idx_k t_total break; end n_pts min(round(0.01 * fs), t_total - idx_k 1); % 单次冲击响应长度约10ms tt (0:n_pts-1)/fs; wd 2*pi*fn*sqrt(1-zeta^2); impact A * exp(-zeta*2*pi*fn*tt) .* sin(wd*tt); idx_end min(idx_k n_pts - 1, t_total); x(idx_k:idx_end) x(idx_k:idx_end) impact(1:(idx_end-idx_k1)); end这样得到的就是一个干净的、没有噪声的外圈故障模拟信号。你会发现时域波形里每隔约9.4ms出现一个衰减振荡幅值基本一致连续性好。把这个信号加一点白噪声做FFT频谱上能看到2000Hz附近有明显的共振峰而做包络谱之后频谱峰就落到了106Hz和它的倍频上这就是外圈故障的“身份证”。3.3 内圈故障模拟传递路径调制与幅值变化内圈故障比外圈麻烦不少。麻烦的核心是传递路径变化。内圈故障点随轴一起旋转当故障点处于承载区时冲击通过滚动体、外圈传到传感器路径较短幅值大当故障点转到承载区对面、靠近传感器远端时信号传过来要绕一大圈幅值就小了。这个变化在时域里看起来就是冲击幅值被“调幅”了调制的频率正好是转频 (f_r)。所以内圈故障信号的构造要在这个地方做文章。我还是沿用前面冲击响应的框架但给每个冲击乘以一个幅值调制函数BPFI 159.5; T_fault 1/BPFI; x zeros(1, t_total); A_base 1.0; mod_depth 0.7; % 调制深度 0~1之间 for k 1:floor(T/T_fault) t_k (k-1) * T_fault; % 计算幅值调制系数随转频做余弦变化 mod_angle 2 * pi * fr * t_k; amp_factor 1 - mod_depth * 0.5 * (1 cos(mod_angle)); % 也可直接用 amp_factor 1 - mod_depth * 0.5 * (1 cos(mod_angle)); A_k A_base * amp_factor; idx_k round(t_k * fs) 1; if idx_k t_total break; end n_pts min(round(0.01 * fs), t_total - idx_k 1); tt (0:n_pts-1)/fs; impact A_k * exp(-zeta*2*pi*fn*tt) .* sin(wd*tt); idx_end min(idx_k n_pts - 1, t_total); x(idx_k:idx_end) x(idx_k:idx_end) impact(1:(idx_end-idx_k1)); end这个调制函数让冲击在转频的每个周期内“忽大忽小”。做包络谱分析的时候你会看到除了BPFI本身及其谐波之外BPFI旁边还出现了以转频为间隔的边带也就是BPFI±fr、BPFI±2fr这些。这个边带结构是内圈故障的重要识别特征也是和外圈故障区分开来的关键。注意调制深度是一个可以调节的参数。实际工况中调制深度取决于载荷分布、径向游隙和传感器位置。模拟的时候建议用0.5到0.8之间的值太浅了边带不明显太深了主冲击几乎消失都不利于后续算法验证。3.4 滚动体故障模拟双重调制与非线性特征滚动体故障是四种故障里最难模拟、也最难识别的。难点在于滚动体一边自转一边公转故障点会交替接触内滚道和外滚道。接触外滚道时轴承座径向方向上的冲击响应比较直接接触内滚道时冲击要经过内圈、滚动体、外圈层层传递信号特征完全不同。所以滚动体故障的冲击不是“每转一圈打一下”而是“每转一圈打两下”——一下打在外滚道一下打在内滚道而且这两下的幅值和传递路径都不一样。这就构成了双重调制一是转频调制二是保持架频率调制。前者来自内圈旋转引起的路径变化后者来自滚动体公转时交替接触内外滚道。工程上实现的时候我习惯把冲击序列分成两组BSF 69.5; % 滚动体特征频率 T_fault 1/BSF; x zeros(1, t_total); fn 3000; % 滚动体故障激起的固有频率通常更高 zeta 0.04; for k 1:floor(T/T_fault) t_k (k-1) * T_fault; % 外滚道接触冲击 A_k1 0.8 * (1 0.5 * cos(2*pi*fr*t_k)); % 内滚道接触冲击幅值和相位都略有不同 A_k2 0.6 * (1 0.5 * cos(2*pi*fr*t_k pi/3)); % 分别叠加两种冲击 % ... 类似前面代码分别用 A_k1 和 A_k2 生成两个冲击并叠加 end两组冲击间隔上是交替的时间间隔是BSF周期的一半。你把这个信号画出来能直观看到冲击幅值呈现“一高一低”的交替模式。包络谱上除了BSF及其谐波还会出现非常丰富的边带频率间隔是转频和保持架频率的叠加。如果不做细致的谱分析滚动体故障很容易被误判为随机冲击干扰。3.5 保持架故障模拟低频调制及其他注意事项保持架故障在工程实际中相对少见但一旦发生往往意味着润滑不良、磨损或者装配问题早期识别很有价值。保持架故障的特征频率FTF很低通常只有转频的0.4倍左右比如30Hz转频下FTF大约12Hz。这么低的频率在频谱上很容易被转频的谐波淹没所以提取难度很大。模拟保持架故障时建模重点是频率调制和幅值不稳定。保持架损伤会导致滚动体运动不规律冲击间隔不再是严格等间距表现为随机性抖动的低频冲击。我一般是在严格周期冲击的基础上引入随机时移FTF 12.0; T_fault 1/FTF; jitter_ratio 0.05; % 随机时移比例 for k 1:floor(T/T_fault) t_k (k-1) * T_fault (rand - 0.5) * jitter_ratio * T_fault; % 其余操作与前面类似幅值可适度随机变化 end这个随机时移让信号变得更加“混沌”也更贴近实际情况。但代价是包络谱上FTF峰会变矮变宽甚至分裂成好几个小峰识别困难。所以做保持架故障分析时不能只看包络谱建议结合时域信号的周期性观察统计脉冲间隔的均值与方差那更有说服力。4. 振动数据分析与故障特征提取实操4.1 模拟信号生成之后的“体检”时域指标怎么看信号生成之后千万别急着做频谱。得先做时域“体检”也就是算一些基础统计量心里有数再继续。我常用的时域指标有峰值、均方根值RMS、峭度Kurtosis、波形因子和峰值因子。早期诊断中峰值因子和峭度对脉冲类故障特别敏感——健康轴承的振动信号接近高斯分布峭度约为3一旦出现局部故障周期性冲击让信号的非高斯性增强峭度会明显升高。外圈故障模拟信号的峭度通常能到15以上内圈受阻于幅值调制峭度略低但也在10附近。如果加噪之后峭度降到5以下就要小心了说明冲击被噪声掩盖常规方法可能看不出来了。在MATLAB里计算峭度的函数现成得很k kurtosis(x); rms_val rms(x); peak_val max(abs(x)); crest_factor peak_val / rms_val;这几个值先记下来为后面写数据分析报告提供数字素材。报告里不放空话用这组数据说明“该模拟信号具有显著脉冲特性”比写一万字都管用。4.2 频域分析FFT、频谱泄漏与窗函数选择时域指标只能告诉你“可能有故障”具体是什么故障还得靠频谱或包络谱。最基础的频域分析就是把时域信号做FFT转换到频域。但这里有个新手必踩的坑频谱泄漏。如果你的采样时长不是信号周期的整数倍或者特征频率与FFT频率分辨率不对齐频谱上就会出现谱峰变宽、幅值变矮的现象甚至出现“假峰”。用一个合适窗函数能有效缓解。我自己的习惯是看整体频谱时用汉宁窗hann兼顾主瓣宽度和旁瓣衰减做冲击类信号的特征提取时用矩形窗反而更直接因为窗函数本身会降低冲击特性做包络谱时先带通滤波再包络解调这时候窗函数的选择影响不大但采样长度要够。N length(x); win hann(N); xf fft(x .* win, N); xf_single 2 * abs(xf(1:N/2)) / N; f_axis (0:N/2-1) * fs / N; plot(f_axis, xf_single);注意加窗之后幅值会有衰减所以谱图上的幅值已经不是真实物理幅值。要做幅值校正就得乘以窗函数的能量修正系数。但做故障诊断一般主要看谱峰位置和相对高低幅值绝对大小不是最重要的事。4.3 包络谱分析希尔伯特变换与故障特征提取包络谱分析是滚动轴承故障诊断的核心手段。思路非常简单轴承故障信号里的冲击成分激起的是高频固有振动故障特征频率本身频率很低几十到几百Hz直接看原始频谱非常吃力因为能量都集中在高频共振区。我们先把高频载波去掉提取它的包络再对包络做FFT低频的故障特征频率就清晰可见了。最常用的包络提取方法是希尔伯特变换。对实信号 (x(t))其解析信号为 (z(t) x(t) j \cdot H[x(t)])其中 (H[x]) 是希尔伯特变换。包络就是解析信号的模。MATLAB里一行代码就能做x_env abs(hilbert(x_filt)); % x_filt 是带通滤波后的信号但要注意直接对整个信号做希尔伯特变换提取包络不一定是最好选择。更规范的做法是先做一个针对固有频率所在频段的带通滤波把共振峰周围的那部分能量保留下来滤掉转频、齿轮啮合等低频干扰和噪声再做希尔伯特变换提取包络最后对包络做FFT得到包络谱。这个方法对于外圈、内圈、滚动体故障都有非常好的效果。我做过一个对比实验同样一组内圈故障模拟信号原始频谱上只能看到2000Hz附近的共振峰和一些杂乱的边带而带通滤波加包络解调之后159.5Hz的BPFI清晰可见旁边还有间隔约30Hz的边带族整个“调制结构”一目了然。4.4 不同故障类型的频谱特征对比一张表说清做完四种故障的模拟和谱分析之后把结果汇总成一张对比表写作报告的必备输出。我自己常用的表格逻辑如下故障类型特征频率调制特征包络谱主要峰位置识别难易程度外圈故障BPFO无明显调制BPFO、2×BPFO、3×BPFO容易内圈故障BPFI转频调制BPFI±n×fr边带中等滚动体故障BSF双重调制交替冲击BSF、BSF±n×fr可能伴随保持架频率成分较难保持架故障FTF频率随机抖动FTF常被淹没观察包络时域周期性困难做区别的时候有几个经验外圈故障的包络谱干干净净峰都落在BPFO的整数倍上内圈故障“峰旁边有卫星峰”边带间隔等于转频滚动体故障的包络谱能量分散有时你会发现BSF的1/2处也会出现峰因为滚动体交替接触内外滚道时每个冲击时间间隔是BSF周期的一半要从时域上验证保持架故障看起来更像随机冲击频率特征不稳定。4.5 仿真“变脏”噪声与干扰的加入做过模拟的人都知道干净的仿真信号用来讲课没问题但如果你要做诊断算法的性能评测或者要验证自己的方法在实际数据上能不能用纯信号反而帮不上忙。所以我在确认模型正确之后一定会做一步加噪声。加噪分几种层次。第一种是白噪声模拟传感器本底噪声和环境随机干扰用MATLAB的awgn函数就能实现信噪比从5dB到20dB逐个试。第二种是周期性干扰比如转频的谐波干扰相当于轴系不平衡或对中不良的干扰成分在时域上表现为低频正弦波。第三种是随机脉冲干扰模拟装配松动或外界敲击产生的偶发冲击这是最难处理的因为它在时域上长得和轴承故障冲击很像但间隔不规律。加噪之后你再拿上面的流程去分析就能体会到从“理想环境”到“实际工程”的差距。我实测下来的数据是信噪比15dB时外圈故障包络谱依然清晰信噪比降到5dB外圈故障还勉强认得出内圈故障的边带基本糊成一片滚动体故障已经几乎没法肉眼识别了。这个“可用性边界”对参数优化和算法选型非常关键。5. 参数敏感性分析与数据分析报告输出要点5.1 转速变化对特征频率的影响规律做完单点工况的模拟之后建议大家一定要做变转速分析。转速这个参数在整个故障诊断里相当于万物的“尺子”。转速一旦变了所有特征频率成比例变化但系统固有频率不变调制边的间隔也随之变化。我在模拟中做过一组转速扫描实验从600rpm到3000rpm每隔200rpm算一次特征频率并同步生成模拟信号。结果非常直观故障特征频率与转速严格线性关系你拿拟合直线的斜率除以理论系数能反过来验证你算的特征频率公式有没有写错。这个交叉验证的方法推荐给所有刚建好模型的人花半个小时能省下后面无数小时的排查时间。5.2 故障尺寸与载荷对幅值的影响故障尺寸主要影响冲击幅值和冲击的持续时间。剥落坑越大滚动体滚过的时候产生的冲击能量越大但冲击持续时间也可能变长激起的高频成分会有所变化。载荷的影响更直接载荷越大滚动体与滚道之间的接触力越大冲击幅值也越大。模拟的时候我一般把故障尺寸映射到冲击幅值 (A) 和阻尼比 (\zeta) 上。小的点蚀对应较高的固有频率激励因为接触刚度变化剧烈大的剥落对应较低频的瞬态响应因为接触过程缓和一些。这些参数之间的关系很难用精确公式描述工程上更多是靠标定数据修正所以我在报告里会明确标注“幅值模拟参数基于经验设定用于趋势研究而非绝对精度预测”。5.3 采样率与数据长度选择频率分辨率是怎么算出来的这个问题看着简单实际栽跟头的人特别多。采样率决定你能看到的最高频率也就是奈奎斯特频率等于采样率的一半。数据长度决定频率分辨率等于采样率除以数据点数。所以采样率越高、数据越长频率分辨率越高能把靠得很近的谱峰区分开。但盲目追求高采样率也是有代价的数据量变大计算变慢而且高频部分往往是噪声主导不利于低频特征提取。轴承冲击信号的固有频率通常在1kHz到10kHz之间所以我一般选25.6kHz采样率对应有效分析带宽12.8kHz足够覆盖轴承的固有频率范围又不至于为了存储和计算发愁。频率分辨率的经验公式是 ( \Delta f fs / N )。如果你需要分辨两根相距5Hz的谱线那数据长度至少要 (fs / 5) 个点。比如25.6kHz采样要分辨5Hz就得至少5120个点加上窗函数的影响适当加长到8000以上。这个逻辑每次建模前先算清楚别等分析完才发现分辨率不够。5.4 数据分析报告怎么写图表组织与结论逻辑最后收口到报告输出。一篇合格的轴承动力学模拟数据分析报告不是把图表堆上去就完了要有清晰的分析链条。我的报告结构一般是这样的问题定义与对象描述轴承型号、工况参数、故障类型说明模型建立过程几何参数、特征频率理论计算、建模假设时域分析结果波形特征、统计指标、脉冲周期观察频域与包络谱分析频谱图、包络谱图、边带结构识别不同故障类型对比按前面那张表组织结论与不足识别出哪些特征、哪些地方受限于模型精度、后续怎么改进。图表方面我强烈建议“每张图只表达一个主要结论”。比如时域图就只展现冲击周期性包络谱图只展现特征频率峰不要一张图上塞好几个信号对比。图注写清楚仿真参数、采样率、滤波范围别让人猜。报告里面所有结论都要能对应到图谱上具体的峰或者结构读者拿着你的报告去对照自己的数据能快速复现这才是有价值的报告。注意写报告时别忘了标注“这是基于理想化假设的动力学模拟结果”尤其是幅值信息不能直接当作现场定量诊断的依据。这个声明不是甩锅而是基本的工程素养。6. 常见问题与排查技巧实录6.1 模拟信号频谱上找不到特征频率最可能的原因这个问题我在这里统一回复因为太常遇到了。排查顺序从简单到复杂第一检查特征频率公式里的参数是不是代错了特别是节圆直径和滚动体直径的单位。这个听起来像废话但我真的见过有人把毫米当成米代进去结果特征频率差了3000倍。第二检查采样时长够不够长。频率分辨率不够特征频率小于一个频率分辨率时谱峰直接消失在背景里尤其保持架故障最容易出现这种情况。第三检查DFT/FFT的频率轴映射。MATLAB的fft输出默认从0到fs之间不处理的话你看到的结果横坐标整体偏移峰值定位不准。第四检查是否用了合适的解调方法。滚动体和内圈故障直接看原始频谱基本看不出名堂必须做包络谱。如果包络谱里也没有多半是特征频率太小、谱线被低频段的巨大能量覆盖了这时需要做高通滤波或者带通滤波。6.2 冲击波形“拖尾”时间太长可能是阻尼比设置出了偏差我在调参的时候经常出现这种情况模拟出的冲击响应衰减太慢前一个冲击还没衰减完后一个冲击就来了时域波形糊成一团。阻尼比 (\zeta) 是这个问题的直接控制参数。取值太小比如0.005衰减极慢冲击叠加严重太大比如0.2衰减太快冲击短得几乎看不到。对不同固有频率合适的阻尼比范围不太一样低频几百赫兹适合稍小一点的阻尼比高频3000赫兹以上可以适度加大。我一般的起点是0.03~0.06然后看着波形调把“衰减到原始幅值的10%所需要的时间”控制在冲击周期的0.5倍左右这样时域波形既清楚又不重叠。6.3 内圈故障的边带辨识有争议怎么用谱图让人信服内圈故障包络谱上的边带是我见过最容易引战的内容因为不同仿真参数下边带数量、间隔、幅值差异很大有的人觉得左边带比右边带高有的人得出相反的结论。实际上这不是代码错了而是调制方式不同。幅值调制本身是非对称边带的边带幅值分布还受传递函数、滤波范围影响。我的建议是报告里不只是贴一张包络谱图而是同时放一段解调后的时域包络波形让大家直观看到“幅值在周期性地忽大忽小”再用包络谱说明边带间隔等于转频。图谱结合结论就扎实了。6.4 随机滑差该怎么引入信号“太完美”也不真实纯周期信号做出来太“教科书”了和实际测量数据差距很大。实际轴承运转时滚动体会产生微小的随机滑移冲击间隔不会严格等间距。为此我在所有模拟信号里都预留了一个“随机时移”选项让每次冲击的时刻在小范围内随机变动jitter jitter_ratio * T_fault * randn; % 随机时移 t_k (k-1) * T_fault jitter;这里的jitter_ratio取0.005到0.02之间也就是说冲击间隔最多变化约2%。这个随机时移对时域波形影响不大但对包络谱有一个显著的“副作用”峰会变宽变矮这和实测信号的频谱特性更像。如果你在做算法测试一定要引入这个环节否则你的算法在仿真数据上表现优秀一到实测数据就崩。6.5 MATLAB版本兼容性中文注释乱码与脚本迁移问题最后分享一个非常实际的问题。不少读者跟我说跑我的脚本遇到中文乱码尤其MATLAB 2023之后的中文注释在部分系统默认编码下会显示成乱码。这个问题本质是文件编码不一致MATLAB在Windows下默认用GBK读文件而源代码是UTF-8就会出乱码。解决办法很简单要么把MATLAB的编码设置改为UTF-8预设项-常规-文件编码要么把所有中文注释统一改成英文要么用专门的代码编辑器把脚本转存成带BOM的UTF-8格式。我个人建议代码里尽量用英文注释因为跑数据的时候中英文输入法切换极其恶心而且分享代码时省去一堆编码纠纷。数据报告里写中文没问题代码就让它格干净一些吧。结尾一点个人的经验体会做轴承动力学模拟这么久我最深的一点体会是仿真模型没有“完不完成”这个说法只有“够不够用”的区别。今天你建一个简单的外圈故障模型只需要半天明天你要分析变转速工况下的轴承故障今天的模型就得加上转速关于时间的函数后天你要做早期微弱故障诊断模型又得考虑滑差和滚动体打滑的细枝末节。这就像搭积木基础打牢了后面往上面加东西就不会乱。如果你打算照着这篇内容做一次完整的轴承动力学模拟我个人建议的起步动作是先拿一组已知轴承型号的数据比如6205深沟球轴承手工算一遍四个特征频率再用MATLAB生成一组无噪的外圈故障信号做一遍完整的时域、频域、包络谱分析把流程走通之后再去碰内圈和滚动体故障。一步一步来每个环节都搞明白“为什么”远远好过直接把代码跑完就觉得自己会了。最后分享一个小技巧做完仿真后去公开数据集比如凯斯西储大学轴承数据中心找同型号轴承的实测数据把实测信号放进你验证过的分析流程里跑一遍。仿真的价值不在仿真本身而在于你拿着仿真训练出来的“眼睛”去看真实世界的时候真的能认出故障来。那一瞬间你会发现前面花的大量时间全都值了。 基于MATLAB的滚动轴承动力学模拟与故障诊断全流程解析
轴承故障诊断这个方向我陆陆续续做了快八年。从最早拿着测振笔在现场一遍遍测数据到后来用MATLAB做动力学模拟反推故障机理这中间最大的感悟是现场实测信号太“脏”了你很难分清哪些特征来自轴承本身哪些是齿轮啮合、工频干扰、结构共振混进来的。这也就是为什么一定要做故障机理建模和数值计算——你先把“干净”的故障信号在仿真里搞明白再去认现场信号就有底气了。这篇内容围绕滚动轴承动力学模拟展开覆盖故障机理建模、数值计算、多故障类型模拟以及最终数据分析报告怎么出。适合正在做故障诊断相关课题的学生、刚接触设备状态监测的工程师或者想用MATLAB把轴承振动信号摸透的入门研究者。我会把自己踩过的坑、验证过的参数、调试思路都写在里面照着做你至少能少走三个月弯路。1. 内容整体设计与思路拆解1.1 为什么选MATLAB做轴承动力学模拟做轴承故障模拟市面上其实有几条路。一是纯数学推导用微分方程描述滚动体与滚道之间的接触力然后自己写数值积分求解这活儿用C、Python都能干但开发周期长二是用商业有限元软件做瞬态动力学分析精度高但建模复杂一个模型折腾一两周很正常三是用MATLAB/Simulink既有成熟的ODE求解器又有强大的信号分析工具箱一套流程下来不用来回切换工具。我自己用过Python试过类似的工作NumPySciPy组合确实灵活但到了信号处理环节MATLAB的Signal Processing Toolbox里现成的函数太方便了尤其是后面做包络谱分析、时频分析的时候一行命令能顶Python十几行代码。还有一个现实因素是很多高校和企业都在用MATLAB你写出来的东西别人能直接跑沟通成本低。就我个人的实践来说MATLAB做轴承动力学模拟的核心优势在于“算法验证”环节。模拟信号生成只是第一步后面要验证包络解调、谱峭度、小波变换这些诊断算法靠不靠谱MATLAB里调试这些特别顺手。后面我会专门讲怎么用这套模拟信号去给算法“打分”。1.2 整体建模方案与数值计算路线完整的轴承动力学模拟流程大致分四步几何与运动学参数定义、故障激励建模、系统振动响应求解、信号分析与故障特征提取。第一步是把自己要模拟的轴承型号参数吃透节圆直径、滚动体直径、滚动体数量、接触角这些直接决定故障特征频率。第二步是建模故障激励滚动轴承的故障本质上是滚动体经过损伤部位时产生周期性冲击这个冲击的幅值、频率、衰减特性是核心。第三步是把冲击激励加到系统模型里求解振动响应。第四步是信号处理把仿真得到的时域信号做FFT、包络谱验证故障特征频率是否和理论值吻合。数值计算路线上我对新手最推荐的方式是“先离散后求解”。就是把振动微分方程用数值积分方法逐步推进MATLAB的ode45是自适应步长的龙格库塔法在处理这类含间歇性冲击的系统时稳定性不错。但要注意冲击信号作用时间极短若步长控制不当很容易漏掉冲击峰值后面我会详细讲怎么处理这个坑。1.3 为什么先建立“简化模型”而不是全三维有限元模型很多初学者上来就想做全三维的有限元模型轴承所有零件全部建出来滚动体、保持架、内外圈一个不落然后跑瞬态动力学分析。坦白说这条路对纯学术研究有价值但对于绝大多数工程诊断需求是过度设计。原因很简单你建模的目的是搞清楚“故障会在振动信号里产生什么特征”而不是“轴承内部应力分布到底有多复杂”。在故障诊断这个场景里我们关心的是振动信号中周期性冲击的出现规律、调制特征、频率分布而这些主要由轴承的运动学关系和接触冲击的周期性决定。所以用集中参数模型或者等效激励模型就足够了。我用过Matlab里自己写简化的“外圈固定、内圈旋转”模型自由度不需要太多本质上是把轴承系统简化为弹簧-质量-阻尼系统滚动体通过故障区域时产生冲击力这个冲击力作用到系统上激起固有频率的瞬态振动。这种模型有明确的物理意义又不需要巨量的计算资源一套参数跑几千转的模拟数据也就是几分钟的事。先把简化模型吃透再去拓展更复杂的模型这个顺序不会错。2. 滚动轴承故障机理与特征频率计算2.1 滚动轴承的几何结构与故障类型划分滚动轴承的基本结构说白了就是内外两个圈中间一堆滚动体滚动体之间还有个保持架防止它们挤在一起。运转的时候内圈随轴旋转滚动体在内外滚道之间做滚动兼滑动运动把轴承载荷传递给外圈和轴承座。在故障诊断里我们主要关注四类故障外圈故障、内圈故障、滚动体故障、保持架故障。这四类故障在振动信号里产生的特征完全不同。外圈故障的信号最好认因为外圈通常固定在轴承座上故障点位置相对传感器不变所以脉冲序列是均匀的内圈故障就复杂了因为故障点跟随转轴旋转冲击力传到传感器时要经过滚动体和外圈路径一直在变信号会出现明显的幅值调制滚动体故障则更隐蔽滚动体自转的同时还在公转故障点会周期性接触内外滚道接触的位置和载荷不断变化信号特征最不容易提取。2.2 故障特征频率的理论计算从运动学关系说起要理解故障特征频率就得先知道滚动体怎么运动。以深沟球轴承为例设节圆直径为 (D)滚动体直径为 (d)滚动体数量为 (Z)接触角为 (\alpha)转频为 (f_r)单位Hz。把轴承当作纯滚动接触来分析内外圈与滚动体之间没有相对滑动就能推出四类故障的特征频率公式。外圈故障特征频率 (BPFO) 的计算公式是[ BPFO \frac{Z}{2} \cdot f_r \cdot \left(1 - \frac{d}{D}\cos\alpha\right) ]内圈故障特征频率 (BPFI) 是[ BPFI \frac{Z}{2} \cdot f_r \cdot \left(1 \frac{d}{D}\cos\alpha\right) ]滚动体故障特征频率 (BSF) 是[ BSF \frac{D}{2d} \cdot f_r \cdot \left(1 - \left(\frac{d}{D}\cos\alpha\right)^2\right) ]保持架故障特征频率 (FTF) 是[ FTF \frac{1}{2} \cdot f_r \cdot \left(1 - \frac{d}{D}\cos\alpha\right) ]注意(BSF) 公式里有时会乘个系数2这是因为滚动体每转一圈会和内外滚道各接触一次但工程上通常取上面这个值作为故障冲击的主频率。实际算的时候这些公式算出来的频率只是理论值真实轴承会因为载荷、游隙、润滑状态等因素产生微小的滑差频率有1%到2%的波动是正常现象。2.3 参数敏感性分析哪些参数对特征频率影响最大我刚刚学这个的时候对着公式看了半天只觉得每个参数都有影响但不知道哪个影响大。后来做模拟做多了慢慢摸出了一些门道。滚动体数量 (Z) 的影响是线性的也是最直接的。其他参数不变的情况下(Z) 越多特征频率越高频谱上的谱线间距也越大故障识别越容易。节圆直径 (D) 和滚动体直径 (d) 的比值也很关键它决定了内外圈特征频率的差异程度。接触角 (\alpha) 在小角度范围内对结果影响不大但到了角接触轴承那种大接触角的情况就要重新考虑了。这里我一般会在MATLAB里写一个参数敏感性分析脚本逐个变量小幅扰动看看特征频率变化率。做这个不是为了写论文凑图而是为了搞清楚如果现场给的轴承参数不太准算出来的特征频率误差到底有多大诊断结论靠不靠谱。这个工作做完你对公式的理解会深很多。2.4 用MATLAB实现特征频率的自动计算既然参数和公式都明确了写一个特征频率计算函数就水到渠成。我自己常用的是这个版本代码不多但足够实用function [BPFO, BPFI, BSF, FTF] bearing_fault_frequencies(D, d, Z, alpha_deg, fr) % 计算滚动轴承四类故障特征频率 % 输入 % D - 节圆直径 (mm) % d - 滚动体直径 (mm) % Z - 滚动体数量 % alpha_deg - 接触角 (度) % fr - 转频 (Hz) % 输出四类故障特征频率 alpha deg2rad(alpha_deg); cos_a cos(alpha); term1 1 - (d / D) * cos_a; term2 1 (d / D) * cos_a; BPFO (Z / 2) * fr * term1; BPFI (Z / 2) * fr * term2; BSF (D / (2 * d)) * fr * (1 - ((d / D) * cos_a)^2); FTF (1 / 2) * fr * term1; end调用的时候比如模拟6205深沟球轴承节圆直径约39.04mm滚动体直径约7.94mm滚动体数量9个接触角取0度转频取29.5Hz约1770rpm算出来BPFO约为106HzBPFI约为159.5Hz。后面模拟信号的时候这些频率就是我要在频谱和包络谱里寻找的目标。注意不同轴承型号的几何参数差异很大同一型号不同厂家的加工误差也不一样。有条件的话优先查轴承手册拿不到精确参数就用近似值但一定得在分析报告里注明。3. 多故障类型的动力学建模与数值计算3.1 冲击激励的数学描述周期脉冲与衰减振荡轴承局部故障的振动响应在机理上可以分解为两个过程一个是故障点经过接触区时产生的冲击力另一个是这个冲击力激起的系统固有振动。先说冲击力。一个理想的局部缺陷比如外滚道上的一个剥落坑滚动体每次滚过这个坑的时候都会产生一个极短时长的冲击脉冲。这个脉冲在频域里能量分布很宽能够激起轴承系统从几百赫兹到几千赫兹的固有频率。工程上常用指数衰减正弦波来描述系统的冲击响应也就是[ x_{i}(t) A \cdot e^{-\zeta \omega_n (t - t_i)} \cdot \sin\left(\omega_d (t - t_i)\right), \quad t \ge t_i ]其中 (\omega_n 2\pi f_n) 是系统固有角频率(\omega_d \omega_n\sqrt{1-\zeta^2}) 是有阻尼固有角频率(\zeta) 是阻尼比(f_n) 是系统固有频率(A) 是冲击幅值(t_i) 是第 (i) 次冲击发生的时刻。这段描述很重要它把轴承故障信号的长相给定下来了一串衰减振荡周期性出现周期等于特征频率的倒数。比如外圈故障BPFO是106Hz那么每1/106秒大概0.00943秒就出现一次冲击。3.2 外圈故障模拟实现最“规矩”的故障信号外圈故障的模拟在所有故障里最简单、也最适合入门。原因正如前面说的外圈固定在轴承座上故障点位置相对传感器恒定每次冲击的传递路径基本不变所以冲击幅值均匀、间隔均匀。在MATLAB里我一般这么构造信号fs 25600; % 采样率 25.6kHz T 2; % 信号时长 2秒 t 0:1/fs:T-1/fs; t_total length(t); x zeros(1, t_total); fr 29.5; % 转频 BPFO 106; % 外圈特征频率 T_fault 1/BPFO; % 故障冲击周期 fn 2000; % 系统固有频率 2kHz zeta 0.05; % 阻尼比 A 1.0; % 冲击幅值 % 生成冲击序列 for k 1:floor(T/T_fault) t_k (k-1) * T_fault; idx_k round(t_k * fs) 1; if idx_k t_total break; end n_pts min(round(0.01 * fs), t_total - idx_k 1); % 单次冲击响应长度约10ms tt (0:n_pts-1)/fs; wd 2*pi*fn*sqrt(1-zeta^2); impact A * exp(-zeta*2*pi*fn*tt) .* sin(wd*tt); idx_end min(idx_k n_pts - 1, t_total); x(idx_k:idx_end) x(idx_k:idx_end) impact(1:(idx_end-idx_k1)); end这样得到的就是一个干净的、没有噪声的外圈故障模拟信号。你会发现时域波形里每隔约9.4ms出现一个衰减振荡幅值基本一致连续性好。把这个信号加一点白噪声做FFT频谱上能看到2000Hz附近有明显的共振峰而做包络谱之后频谱峰就落到了106Hz和它的倍频上这就是外圈故障的“身份证”。3.3 内圈故障模拟传递路径调制与幅值变化内圈故障比外圈麻烦不少。麻烦的核心是传递路径变化。内圈故障点随轴一起旋转当故障点处于承载区时冲击通过滚动体、外圈传到传感器路径较短幅值大当故障点转到承载区对面、靠近传感器远端时信号传过来要绕一大圈幅值就小了。这个变化在时域里看起来就是冲击幅值被“调幅”了调制的频率正好是转频 (f_r)。所以内圈故障信号的构造要在这个地方做文章。我还是沿用前面冲击响应的框架但给每个冲击乘以一个幅值调制函数BPFI 159.5; T_fault 1/BPFI; x zeros(1, t_total); A_base 1.0; mod_depth 0.7; % 调制深度 0~1之间 for k 1:floor(T/T_fault) t_k (k-1) * T_fault; % 计算幅值调制系数随转频做余弦变化 mod_angle 2 * pi * fr * t_k; amp_factor 1 - mod_depth * 0.5 * (1 cos(mod_angle)); % 也可直接用 amp_factor 1 - mod_depth * 0.5 * (1 cos(mod_angle)); A_k A_base * amp_factor; idx_k round(t_k * fs) 1; if idx_k t_total break; end n_pts min(round(0.01 * fs), t_total - idx_k 1); tt (0:n_pts-1)/fs; impact A_k * exp(-zeta*2*pi*fn*tt) .* sin(wd*tt); idx_end min(idx_k n_pts - 1, t_total); x(idx_k:idx_end) x(idx_k:idx_end) impact(1:(idx_end-idx_k1)); end这个调制函数让冲击在转频的每个周期内“忽大忽小”。做包络谱分析的时候你会看到除了BPFI本身及其谐波之外BPFI旁边还出现了以转频为间隔的边带也就是BPFI±fr、BPFI±2fr这些。这个边带结构是内圈故障的重要识别特征也是和外圈故障区分开来的关键。注意调制深度是一个可以调节的参数。实际工况中调制深度取决于载荷分布、径向游隙和传感器位置。模拟的时候建议用0.5到0.8之间的值太浅了边带不明显太深了主冲击几乎消失都不利于后续算法验证。3.4 滚动体故障模拟双重调制与非线性特征滚动体故障是四种故障里最难模拟、也最难识别的。难点在于滚动体一边自转一边公转故障点会交替接触内滚道和外滚道。接触外滚道时轴承座径向方向上的冲击响应比较直接接触内滚道时冲击要经过内圈、滚动体、外圈层层传递信号特征完全不同。所以滚动体故障的冲击不是“每转一圈打一下”而是“每转一圈打两下”——一下打在外滚道一下打在内滚道而且这两下的幅值和传递路径都不一样。这就构成了双重调制一是转频调制二是保持架频率调制。前者来自内圈旋转引起的路径变化后者来自滚动体公转时交替接触内外滚道。工程上实现的时候我习惯把冲击序列分成两组BSF 69.5; % 滚动体特征频率 T_fault 1/BSF; x zeros(1, t_total); fn 3000; % 滚动体故障激起的固有频率通常更高 zeta 0.04; for k 1:floor(T/T_fault) t_k (k-1) * T_fault; % 外滚道接触冲击 A_k1 0.8 * (1 0.5 * cos(2*pi*fr*t_k)); % 内滚道接触冲击幅值和相位都略有不同 A_k2 0.6 * (1 0.5 * cos(2*pi*fr*t_k pi/3)); % 分别叠加两种冲击 % ... 类似前面代码分别用 A_k1 和 A_k2 生成两个冲击并叠加 end两组冲击间隔上是交替的时间间隔是BSF周期的一半。你把这个信号画出来能直观看到冲击幅值呈现“一高一低”的交替模式。包络谱上除了BSF及其谐波还会出现非常丰富的边带频率间隔是转频和保持架频率的叠加。如果不做细致的谱分析滚动体故障很容易被误判为随机冲击干扰。3.5 保持架故障模拟低频调制及其他注意事项保持架故障在工程实际中相对少见但一旦发生往往意味着润滑不良、磨损或者装配问题早期识别很有价值。保持架故障的特征频率FTF很低通常只有转频的0.4倍左右比如30Hz转频下FTF大约12Hz。这么低的频率在频谱上很容易被转频的谐波淹没所以提取难度很大。模拟保持架故障时建模重点是频率调制和幅值不稳定。保持架损伤会导致滚动体运动不规律冲击间隔不再是严格等间距表现为随机性抖动的低频冲击。我一般是在严格周期冲击的基础上引入随机时移FTF 12.0; T_fault 1/FTF; jitter_ratio 0.05; % 随机时移比例 for k 1:floor(T/T_fault) t_k (k-1) * T_fault (rand - 0.5) * jitter_ratio * T_fault; % 其余操作与前面类似幅值可适度随机变化 end这个随机时移让信号变得更加“混沌”也更贴近实际情况。但代价是包络谱上FTF峰会变矮变宽甚至分裂成好几个小峰识别困难。所以做保持架故障分析时不能只看包络谱建议结合时域信号的周期性观察统计脉冲间隔的均值与方差那更有说服力。4. 振动数据分析与故障特征提取实操4.1 模拟信号生成之后的“体检”时域指标怎么看信号生成之后千万别急着做频谱。得先做时域“体检”也就是算一些基础统计量心里有数再继续。我常用的时域指标有峰值、均方根值RMS、峭度Kurtosis、波形因子和峰值因子。早期诊断中峰值因子和峭度对脉冲类故障特别敏感——健康轴承的振动信号接近高斯分布峭度约为3一旦出现局部故障周期性冲击让信号的非高斯性增强峭度会明显升高。外圈故障模拟信号的峭度通常能到15以上内圈受阻于幅值调制峭度略低但也在10附近。如果加噪之后峭度降到5以下就要小心了说明冲击被噪声掩盖常规方法可能看不出来了。在MATLAB里计算峭度的函数现成得很k kurtosis(x); rms_val rms(x); peak_val max(abs(x)); crest_factor peak_val / rms_val;这几个值先记下来为后面写数据分析报告提供数字素材。报告里不放空话用这组数据说明“该模拟信号具有显著脉冲特性”比写一万字都管用。4.2 频域分析FFT、频谱泄漏与窗函数选择时域指标只能告诉你“可能有故障”具体是什么故障还得靠频谱或包络谱。最基础的频域分析就是把时域信号做FFT转换到频域。但这里有个新手必踩的坑频谱泄漏。如果你的采样时长不是信号周期的整数倍或者特征频率与FFT频率分辨率不对齐频谱上就会出现谱峰变宽、幅值变矮的现象甚至出现“假峰”。用一个合适窗函数能有效缓解。我自己的习惯是看整体频谱时用汉宁窗hann兼顾主瓣宽度和旁瓣衰减做冲击类信号的特征提取时用矩形窗反而更直接因为窗函数本身会降低冲击特性做包络谱时先带通滤波再包络解调这时候窗函数的选择影响不大但采样长度要够。N length(x); win hann(N); xf fft(x .* win, N); xf_single 2 * abs(xf(1:N/2)) / N; f_axis (0:N/2-1) * fs / N; plot(f_axis, xf_single);注意加窗之后幅值会有衰减所以谱图上的幅值已经不是真实物理幅值。要做幅值校正就得乘以窗函数的能量修正系数。但做故障诊断一般主要看谱峰位置和相对高低幅值绝对大小不是最重要的事。4.3 包络谱分析希尔伯特变换与故障特征提取包络谱分析是滚动轴承故障诊断的核心手段。思路非常简单轴承故障信号里的冲击成分激起的是高频固有振动故障特征频率本身频率很低几十到几百Hz直接看原始频谱非常吃力因为能量都集中在高频共振区。我们先把高频载波去掉提取它的包络再对包络做FFT低频的故障特征频率就清晰可见了。最常用的包络提取方法是希尔伯特变换。对实信号 (x(t))其解析信号为 (z(t) x(t) j \cdot H[x(t)])其中 (H[x]) 是希尔伯特变换。包络就是解析信号的模。MATLAB里一行代码就能做x_env abs(hilbert(x_filt)); % x_filt 是带通滤波后的信号但要注意直接对整个信号做希尔伯特变换提取包络不一定是最好选择。更规范的做法是先做一个针对固有频率所在频段的带通滤波把共振峰周围的那部分能量保留下来滤掉转频、齿轮啮合等低频干扰和噪声再做希尔伯特变换提取包络最后对包络做FFT得到包络谱。这个方法对于外圈、内圈、滚动体故障都有非常好的效果。我做过一个对比实验同样一组内圈故障模拟信号原始频谱上只能看到2000Hz附近的共振峰和一些杂乱的边带而带通滤波加包络解调之后159.5Hz的BPFI清晰可见旁边还有间隔约30Hz的边带族整个“调制结构”一目了然。4.4 不同故障类型的频谱特征对比一张表说清做完四种故障的模拟和谱分析之后把结果汇总成一张对比表写作报告的必备输出。我自己常用的表格逻辑如下故障类型特征频率调制特征包络谱主要峰位置识别难易程度外圈故障BPFO无明显调制BPFO、2×BPFO、3×BPFO容易内圈故障BPFI转频调制BPFI±n×fr边带中等滚动体故障BSF双重调制交替冲击BSF、BSF±n×fr可能伴随保持架频率成分较难保持架故障FTF频率随机抖动FTF常被淹没观察包络时域周期性困难做区别的时候有几个经验外圈故障的包络谱干干净净峰都落在BPFO的整数倍上内圈故障“峰旁边有卫星峰”边带间隔等于转频滚动体故障的包络谱能量分散有时你会发现BSF的1/2处也会出现峰因为滚动体交替接触内外滚道时每个冲击时间间隔是BSF周期的一半要从时域上验证保持架故障看起来更像随机冲击频率特征不稳定。4.5 仿真“变脏”噪声与干扰的加入做过模拟的人都知道干净的仿真信号用来讲课没问题但如果你要做诊断算法的性能评测或者要验证自己的方法在实际数据上能不能用纯信号反而帮不上忙。所以我在确认模型正确之后一定会做一步加噪声。加噪分几种层次。第一种是白噪声模拟传感器本底噪声和环境随机干扰用MATLAB的awgn函数就能实现信噪比从5dB到20dB逐个试。第二种是周期性干扰比如转频的谐波干扰相当于轴系不平衡或对中不良的干扰成分在时域上表现为低频正弦波。第三种是随机脉冲干扰模拟装配松动或外界敲击产生的偶发冲击这是最难处理的因为它在时域上长得和轴承故障冲击很像但间隔不规律。加噪之后你再拿上面的流程去分析就能体会到从“理想环境”到“实际工程”的差距。我实测下来的数据是信噪比15dB时外圈故障包络谱依然清晰信噪比降到5dB外圈故障还勉强认得出内圈故障的边带基本糊成一片滚动体故障已经几乎没法肉眼识别了。这个“可用性边界”对参数优化和算法选型非常关键。5. 参数敏感性分析与数据分析报告输出要点5.1 转速变化对特征频率的影响规律做完单点工况的模拟之后建议大家一定要做变转速分析。转速这个参数在整个故障诊断里相当于万物的“尺子”。转速一旦变了所有特征频率成比例变化但系统固有频率不变调制边的间隔也随之变化。我在模拟中做过一组转速扫描实验从600rpm到3000rpm每隔200rpm算一次特征频率并同步生成模拟信号。结果非常直观故障特征频率与转速严格线性关系你拿拟合直线的斜率除以理论系数能反过来验证你算的特征频率公式有没有写错。这个交叉验证的方法推荐给所有刚建好模型的人花半个小时能省下后面无数小时的排查时间。5.2 故障尺寸与载荷对幅值的影响故障尺寸主要影响冲击幅值和冲击的持续时间。剥落坑越大滚动体滚过的时候产生的冲击能量越大但冲击持续时间也可能变长激起的高频成分会有所变化。载荷的影响更直接载荷越大滚动体与滚道之间的接触力越大冲击幅值也越大。模拟的时候我一般把故障尺寸映射到冲击幅值 (A) 和阻尼比 (\zeta) 上。小的点蚀对应较高的固有频率激励因为接触刚度变化剧烈大的剥落对应较低频的瞬态响应因为接触过程缓和一些。这些参数之间的关系很难用精确公式描述工程上更多是靠标定数据修正所以我在报告里会明确标注“幅值模拟参数基于经验设定用于趋势研究而非绝对精度预测”。5.3 采样率与数据长度选择频率分辨率是怎么算出来的这个问题看着简单实际栽跟头的人特别多。采样率决定你能看到的最高频率也就是奈奎斯特频率等于采样率的一半。数据长度决定频率分辨率等于采样率除以数据点数。所以采样率越高、数据越长频率分辨率越高能把靠得很近的谱峰区分开。但盲目追求高采样率也是有代价的数据量变大计算变慢而且高频部分往往是噪声主导不利于低频特征提取。轴承冲击信号的固有频率通常在1kHz到10kHz之间所以我一般选25.6kHz采样率对应有效分析带宽12.8kHz足够覆盖轴承的固有频率范围又不至于为了存储和计算发愁。频率分辨率的经验公式是 ( \Delta f fs / N )。如果你需要分辨两根相距5Hz的谱线那数据长度至少要 (fs / 5) 个点。比如25.6kHz采样要分辨5Hz就得至少5120个点加上窗函数的影响适当加长到8000以上。这个逻辑每次建模前先算清楚别等分析完才发现分辨率不够。5.4 数据分析报告怎么写图表组织与结论逻辑最后收口到报告输出。一篇合格的轴承动力学模拟数据分析报告不是把图表堆上去就完了要有清晰的分析链条。我的报告结构一般是这样的问题定义与对象描述轴承型号、工况参数、故障类型说明模型建立过程几何参数、特征频率理论计算、建模假设时域分析结果波形特征、统计指标、脉冲周期观察频域与包络谱分析频谱图、包络谱图、边带结构识别不同故障类型对比按前面那张表组织结论与不足识别出哪些特征、哪些地方受限于模型精度、后续怎么改进。图表方面我强烈建议“每张图只表达一个主要结论”。比如时域图就只展现冲击周期性包络谱图只展现特征频率峰不要一张图上塞好几个信号对比。图注写清楚仿真参数、采样率、滤波范围别让人猜。报告里面所有结论都要能对应到图谱上具体的峰或者结构读者拿着你的报告去对照自己的数据能快速复现这才是有价值的报告。注意写报告时别忘了标注“这是基于理想化假设的动力学模拟结果”尤其是幅值信息不能直接当作现场定量诊断的依据。这个声明不是甩锅而是基本的工程素养。6. 常见问题与排查技巧实录6.1 模拟信号频谱上找不到特征频率最可能的原因这个问题我在这里统一回复因为太常遇到了。排查顺序从简单到复杂第一检查特征频率公式里的参数是不是代错了特别是节圆直径和滚动体直径的单位。这个听起来像废话但我真的见过有人把毫米当成米代进去结果特征频率差了3000倍。第二检查采样时长够不够长。频率分辨率不够特征频率小于一个频率分辨率时谱峰直接消失在背景里尤其保持架故障最容易出现这种情况。第三检查DFT/FFT的频率轴映射。MATLAB的fft输出默认从0到fs之间不处理的话你看到的结果横坐标整体偏移峰值定位不准。第四检查是否用了合适的解调方法。滚动体和内圈故障直接看原始频谱基本看不出名堂必须做包络谱。如果包络谱里也没有多半是特征频率太小、谱线被低频段的巨大能量覆盖了这时需要做高通滤波或者带通滤波。6.2 冲击波形“拖尾”时间太长可能是阻尼比设置出了偏差我在调参的时候经常出现这种情况模拟出的冲击响应衰减太慢前一个冲击还没衰减完后一个冲击就来了时域波形糊成一团。阻尼比 (\zeta) 是这个问题的直接控制参数。取值太小比如0.005衰减极慢冲击叠加严重太大比如0.2衰减太快冲击短得几乎看不到。对不同固有频率合适的阻尼比范围不太一样低频几百赫兹适合稍小一点的阻尼比高频3000赫兹以上可以适度加大。我一般的起点是0.03~0.06然后看着波形调把“衰减到原始幅值的10%所需要的时间”控制在冲击周期的0.5倍左右这样时域波形既清楚又不重叠。6.3 内圈故障的边带辨识有争议怎么用谱图让人信服内圈故障包络谱上的边带是我见过最容易引战的内容因为不同仿真参数下边带数量、间隔、幅值差异很大有的人觉得左边带比右边带高有的人得出相反的结论。实际上这不是代码错了而是调制方式不同。幅值调制本身是非对称边带的边带幅值分布还受传递函数、滤波范围影响。我的建议是报告里不只是贴一张包络谱图而是同时放一段解调后的时域包络波形让大家直观看到“幅值在周期性地忽大忽小”再用包络谱说明边带间隔等于转频。图谱结合结论就扎实了。6.4 随机滑差该怎么引入信号“太完美”也不真实纯周期信号做出来太“教科书”了和实际测量数据差距很大。实际轴承运转时滚动体会产生微小的随机滑移冲击间隔不会严格等间距。为此我在所有模拟信号里都预留了一个“随机时移”选项让每次冲击的时刻在小范围内随机变动jitter jitter_ratio * T_fault * randn; % 随机时移 t_k (k-1) * T_fault jitter;这里的jitter_ratio取0.005到0.02之间也就是说冲击间隔最多变化约2%。这个随机时移对时域波形影响不大但对包络谱有一个显著的“副作用”峰会变宽变矮这和实测信号的频谱特性更像。如果你在做算法测试一定要引入这个环节否则你的算法在仿真数据上表现优秀一到实测数据就崩。6.5 MATLAB版本兼容性中文注释乱码与脚本迁移问题最后分享一个非常实际的问题。不少读者跟我说跑我的脚本遇到中文乱码尤其MATLAB 2023之后的中文注释在部分系统默认编码下会显示成乱码。这个问题本质是文件编码不一致MATLAB在Windows下默认用GBK读文件而源代码是UTF-8就会出乱码。解决办法很简单要么把MATLAB的编码设置改为UTF-8预设项-常规-文件编码要么把所有中文注释统一改成英文要么用专门的代码编辑器把脚本转存成带BOM的UTF-8格式。我个人建议代码里尽量用英文注释因为跑数据的时候中英文输入法切换极其恶心而且分享代码时省去一堆编码纠纷。数据报告里写中文没问题代码就让它格干净一些吧。结尾一点个人的经验体会做轴承动力学模拟这么久我最深的一点体会是仿真模型没有“完不完成”这个说法只有“够不够用”的区别。今天你建一个简单的外圈故障模型只需要半天明天你要分析变转速工况下的轴承故障今天的模型就得加上转速关于时间的函数后天你要做早期微弱故障诊断模型又得考虑滑差和滚动体打滑的细枝末节。这就像搭积木基础打牢了后面往上面加东西就不会乱。如果你打算照着这篇内容做一次完整的轴承动力学模拟我个人建议的起步动作是先拿一组已知轴承型号的数据比如6205深沟球轴承手工算一遍四个特征频率再用MATLAB生成一组无噪的外圈故障信号做一遍完整的时域、频域、包络谱分析把流程走通之后再去碰内圈和滚动体故障。一步一步来每个环节都搞明白“为什么”远远好过直接把代码跑完就觉得自己会了。最后分享一个小技巧做完仿真后去公开数据集比如凯斯西储大学轴承数据中心找同型号轴承的实测数据把实测信号放进你验证过的分析流程里跑一遍。仿真的价值不在仿真本身而在于你拿着仿真训练出来的“眼睛”去看真实世界的时候真的能认出故障来。那一瞬间你会发现前面花的大量时间全都值了。