基于包络分析的振动故障诊断：从原理到MATLAB实现

简介：这份振动故障诊断MATLAB源码包面向机械健康监测和故障预测方向的工程师、研究人员及学生，以实际可运行的m脚本和说明文档，演示从振动信号预处理、时域/频域/复频域分析到特征提取与模型训练识别的完整流程。压缩包共26个文件，包含25个m文件与1个txt说明，大小仅49KB，其中m文件覆盖快速傅里叶变换、小波分析、希尔伯特变换、峭度/偏度/方差等统计量计算及支持向量机等分类算法的建模调用，txt用于参数与操作指引，适合对照理论学习并动手复现。目前已有337人学习下载。通过运行其中的主程序与各功能模块，读者可直观掌握故障特征频率识别、包络分析、概率密度估计等关键技巧，并能借助代码快速修改参数，迁移到轴承、齿轮等旋转机械的振动诊断场景，提升实际工程中的故障预判与定位能力。

1. 振动故障诊断为什么不能只看频谱图

机械设备出故障时，最先响应的往往不是温度、不是电流，而是振动。滚动轴承内圈出现剥落、齿轮齿面发生点蚀、转子产生不对中时，每转一圈都会激发一次冲击，这些冲击以极短的脉冲形式叠加在机器正常运行的高频振动里。问题在于：直接对原始振动信号做FFT，看到的频谱图往往只有一条宽泛的共振峰，故障频率被淹没在背景噪声和结构共振中，根本无法直接读出“内圈故障BPFI=98.7Hz”这类结论。这恰恰是振动故障诊断最反直觉的地方——故障信息不在振动幅度最大的地方，而是藏在幅值调制形成的边带里。

本文要解决的就是这件事：从振动故障诊断的物理原理出发，讲清楚为什么冲击性故障会产生调制、为什么包络解调能把故障频率从高频载波中剥离出来，然后给出用MATLAB从数据导入、滤波、特征提取到包络谱识别的完整落地路径。适合正在做设备状态监测的工程师、研究生以及需要从零搭建诊断算法的开发者。全文不依赖任何商业诊断软件，只用一个MATLAB脚本就能在真实轴承数据上跑通。

2. 振动故障诊断的底层原理：从冲击调制到包络解调

2.1 故障为什么“藏”在调制信号里

滚动轴承的滚动体滚过局部缺陷时，会产生一个极短时宽的冲击脉冲。这个脉冲的频谱非常宽，能覆盖到数千甚至上万赫兹，会激励起轴承座、端盖、传感器安装结构的高频固有共振。传感器拾取的信号因此由两部分叠加而成：一是跟转速同步的周期性冲击序列，频率为故障特征频率（几十到几百赫兹）；二是被冲击激发的高频共振衰减振荡（几千赫兹以上）。

机械系统可以看作一个线性系统，冲击是激励，共振是响应。响应幅度受到冲击强度的“调制”——故障越严重、间隙越大、载荷越高，共振波的幅值变化越明显。从数学上看，传感器信号近似为：

x(t) = A(t) · cos(2πf_res t)

其中 f_res 是结构共振频率，A(t) 是随时间变化的幅值包络，它周期性重复，重复频率等于故障特征频率。这个模型就是幅值调制：低频的故障特征“骑”在高频载波上。直接看原始波形，可以看到高频振荡的幅度一节一节变化，但直接用FFT却很难看清低频的调制频率——因为FFT输出的能量几乎全部集中在高频共振峰附近，低频分量相对太弱。

2.2 傅里叶变换能看到什么，看不到什么

FFT适合处理平稳、周期信号。当机器稳定运转时，振动信号近似平稳，FFT能把能量按频率分布展示出来，这是它的优势。但对冲击性故障信号，FFT有两个先天不足。

其一，冲击信号不是正弦波，它包含丰富的谐波成分。频谱图上，冲击的周期性会转化为故障特征频率及其整数倍谐波，但这些谱线能量很低，容易淹没在噪声里。其二，共振峰本身很宽，故障冲击激发的是衰减振荡，相当于给共振峰“镀”了一层额外带宽，谱线变糊，附近出现密集的边带。实际工程中，如果转速波动、载荷变化，频谱这种“全局平均”的手段会把瞬态冲击抹平——故障特征彻底消失。

所以在振动故障诊断实践中，直接对原始信号做FFT只能作为粗筛手段，用来判断“有没有异常”，比如总能量是否变大、高频段是否有明显峰起。要想定位故障部位（内圈、外圈、滚动体还是齿轮齿面），必须把调制信号解调出来，这就是包络分析的核心价值。

2.3 包络谱为什么能把故障频率“搬”到低频区

包络分析的本质是去掉高频载波，只保留幅值调制信号。最常用的数学工具是希尔伯特变换：对原始信号 x(t)，构造解析信号

z(t) = x(t) + j·H{x(t)}

其中 H{x(t)} 是 x(t) 的希尔伯特变换。解析信号的幅值

|z(t)| = sqrt(x²(t) + H{x(t)}²)

就是信号在每一时刻的瞬时幅值，也就是包络。实现时用MATLAB内置的hilbert函数，一行代码即可得到复数解析信号，再取abs就得到包络。

对包络信号再做一次FFT，得到的就是包络谱。这个包络谱的低频段会出现明显的谱峰，位于故障特征频率 f_BPFO、f_BPFI、f_BSF 及其谐波处。之所以说“搬”到低频，是因为包络信号已经不含载波成分，频率轴上关注的区间从几kHz降到了几百Hz以内。这个频段的谱线识别难度远低于直接从原始频谱里找边带，而且对转速波动不敏感，工程上非常实用。

这里要特别强调一个细节：包络分析要求先做带通滤波，把共振频带单独切出来，再做希尔伯特变换。如果对全频段信号直接求包络，会把各种频率的调制混叠在一起，包络谱上出现大量无关峰值。所以实际流程是：带通滤波 → 包络提取 → FFT，三步缺一不可。

3. 用MATLAB搭建振动故障诊断的预处理与特征提取流水线

3.1 把采集器数据导入工作区并做有效性检查

振动故障诊断的第一步永远是数据质量检查，而不是直接滤波。我从现场拿到的数据常来自NI采集卡、LMS或者工业级状态监测系统，格式各异，但导入MATLAB后第一件事都是确认三件事：采样率、信号长度、单位。

% 清空工作区并导入数据 data = csvread('bearing_vibration.csv', 1, 0); % 跳过首行表头，读取振动数据 fs = 25600; % 采样率，单位Hz，与采集器设置一致 t = (0:length(data)-1) / fs; % 构造时间轴 % 有效性检查：是否有NaN、是否饱和、是否有异常尖峰 if any(isnan(data)) error('数据包含NaN值，请检查采集通道'); end if max(abs(data)) > 10 warning('信号幅值偏大，请确认单位是否为g或m/s^2'); end % 去掉直流分量：加速度传感器零漂会引起均值偏移 data = data - mean(data); % 快速可视化确认信号形态 figure; plot(t(1:2000), data(1:2000)); xlabel('时间 (s)'); ylabel('加速度 (g)'); title('原始振动信号前0.08秒');

这段脚本的关键在于采样率必须与硬件一致。振动诊断中采样率决定可分析的最高频率，如果要分析轴承外圈故障特征频率（通常在几十到几百Hz），采样率不需要太高，但要能覆盖共振频带（一般需要至少5kHz以上），而共振频带的截止频率直接决定带通滤波的参数选择。

3.2 去趋势、去均值与抗混叠滤波的先后顺序

预处理顺序有讲究。先去均值、再去趋势、最后滤波，这是标准做法。

去趋势用detrend函数，它消除的是信号中的线性或多项式趋势项。在工业现场，由于温度漂移、传感器应变、基座松动，采集到的信号往往叠加上一个缓变的基线漂移。如果不处理，FFT的低频段会出现超大的虚假分量，包络谱低频段也会被“抬高”，导致故障特征频率附近的弱峰值被掩盖。

% 去趋势：去除缓变基线漂移 data_detrended = detrend(data, 'linear'); % 带通滤波：先定位共振频带，再滤波 low_freq = 2000; % 带通下限，单位Hz，根据共振频带调整 high_freq = 6000; % 带通上限，单位Hz [b, a] = butter(4, [low_freq high_freq]/(fs/2), 'bandpass'); data_filtered = filtfilt(b, a, data_detrended);

这里过滤波器截止频率选择的物理依据是：必须先确认结构共振发生在哪个频段。实际做法是先对原始信号做一次全频段频谱分析，找到能量显著抬升的频带，把这个频带作为带通范围。filtfilt做零相位滤波，避免滤波本身引入相位偏移导致包络的时域定位不准。但这个滤波器只是整个诊断流水线中的第3步，做完后信号已经变成窄带的共振响应，可以放心做包络分析。

3.3 时域特征公式的MATLAB实现：峰值因子、峭度与均方根

包络谱是诊断的利器，但工程现场不能每个测点都立刻做谱分析。时域特征参数作为快速筛查手段，0.1秒就能算完一个测点的健康程度，适合在线状态监测系统做报警判定。工程上最常用的是三个互补指标：均方根值、峰值因子、峭度。

function features = time_features(signal) % 输入signal为预处理后的振动信号 % 输出features = [rms, peak_factor, kurtosis] sig_rms = rms(signal); sig_peak = max(abs(signal)); peak_factor = sig_peak / sig_rms; sig_kurtosis = kurtosis(signal); features = [sig_rms, peak_factor, sig_kurtosis]; end

这三个参数的物理含义完全不同：均方根值衡量振动能量总体水平，对磨损类故障敏感但早期故障几乎无反应；峰值因子是峰值与均方根的比值，对单次冲击敏感，轴承早期剥落往往表现为峰值因子突然增大；峭度是四阶统计矩归一化结果，反映信号分布的“重尾”程度。正常振动近似高斯分布，峭度约等于3；轴承早期故障冲击是少数极端事件，峭度能飙升到5以上。

实际使用时要组合判断：如果均方根上升明显且峭度同步升高，大概率是磨损加剧；如果均方根变化不大但峰值因子和峭度升高，更可能是早期冲击性故障。这个区分比单独看某个指标可靠得多，而且计算成本极低，可以在嵌入式系统上跑。

4. 滚动轴承故障诊断实战：包络谱与特征频率对照

4.1 计算轴承故障特征频率的4个公式与参数表

滚动轴承的局部缺陷在不同元件上通过时会产生不同的冲击重复频率。工程上最常用的是四个特征频率公式，全部由几何参数和转速推导：外圈固定时滚动体通过外圈缺陷的频率为BPFO，通过内圈缺陷的频率为BPFI，滚动体自身缺陷频率为BSF，保持架转频为FTF。

其中 fr 为转频（Hz），Z 为滚动体数量，d 为滚动体直径，D 为节圆直径，α 为接触角。这些参数在轴承型号手册里都能查到，不需要高精度测量。以某型号深沟球轴承为例，Z=9，d=7.94mm，D=39.04mm，α=0°，转频25Hz时，BPFO≈92.1Hz，BPFI≈132.9Hz。

% 轴承特征频率计算函数 function freq = bearing_frequencies(fr, Z, d, D, alpha) % 参数说明： % fr : 转轴频率 (Hz) % Z : 滚动体数量 % d : 滚动体直径 (mm) % D : 节圆直径 (mm) % alpha : 接触角 (rad)，无接触角时填0 cos_angle = cos(alpha); FTF = (fr/2) * (1 - (d/D)*cos_angle); BPFO = (Z*fr/2) * (1 - (d/D)*cos_angle); BPFI = (Z*fr/2) * (1 + (d/D)*cos_angle); BSF = (D*fr/(2*d)) * (1 - (d/D)^2 * cos_angle^2); freq = [FTF, BPFO, BPFI, BSF]; end

这里有个必须靠经验才会注意的陷阱：公式里的 d/D 是个无量纲比值，只要直径用同一单位就成立，但不同轴承手册给的接触角不完全一致，计算结果的误差在±1Hz内都算正常。包络谱找峰时给5Hz左右的容差区间比死盯单个频率点要稳得多。

4.2 用envelope与FFT提取包络谱并识别故障特征

实际环境采集的振动信号噪声大、干扰成分复杂，包络谱分析前必须先滤波。带着共振频带的尖峰，包络谱上会看到鲜明的调制谱。我一般用下列脚本处理一组信号：

function [env_spec, f_axis] = envelope_spectrum(signal, fs, fl, fh) % 输入： % signal : 预处理后的振动信号 % fs : 采样率 % fl : 带通滤波下限 % fh : 带通滤波上限 % 输出： % env_spec : 归一化包络谱幅值 % f_axis : 频率轴 % 带通滤波提取共振频带 [b, a] = butter(4, [fl fh]/(fs/2), 'bandpass'); sig_band = filtfilt(b, a, signal); % 提取包络 env = abs(hilbert(sig_band)); % 对包络做FFT，取单边谱 L = length(env); NFFT = 2^nextpow2(L); Y = fft(env, NFFT) / L; f_axis = fs/2 * linspace(0, 1, NFFT/2+1); env_spec = 2 * abs(Y(1:NFFT/2+1)); % 平滑处理，方便找峰 env_spec = movmean(env_spec, 5); end

这段代码的工程要点有三个。第一，hilbert取abs后得到的包络包含直流分量和故障特征频率分量，FFT后直流分量会非常高，所以要看0Hz以上的谱峰时必须忽略第一个点。第二，包络谱的谱线密度取决于数据长度，数据越短，频率分辨率越粗，短数据下BPFO和BPFI可能完全分不开，所以要尽量截取足够长的平稳段，建议至少采集40个转频周期以上。第三，movmean做5点平滑可以去掉窄带毛刺，但窗口太大会把相邻谱峰合并，轴承故障谱峰本来就窄，不宜过度平滑。

结束后在包络谱上搜索峰值，逐个与理论BPFO、BPFI、BSF对比。外圈故障在BPFO处出现主峰，内圈故障在BPFI处出现主峰外还会在BPFI两侧各间隔一个转频处出现边带，滚动体故障在BSF处同样有边带且常伴随2倍BSF谐波。如果所有理论频率处都看不到明显峰值、只有宽泛的谱线抬高，那更可能是润滑不良或松动而非局部缺陷。

4.3 边频带识别：为什么只看主峰值会漏掉内圈故障

内圈故障和转轴的旋转有关系，缺陷位置每转一圈会因为承载区变化而改变冲击强度，于是内圈故障的包络谱在BPFI主峰两侧会出现以转频fr为间距的边带。只看主峰而忽视边带，很容易将内圈故障误判为滚动体故障——因为滚动体故障也会有类似但更窄的边带结构。

% 在包络谱中查找边带 % f_bpfi 为理论内圈特征频率 fr = 25; % 转频 search_range = f_bpfi + [-2, 2]; % 主峰搜索范围 [~, idx_main] = findpeaks(env_spec(f_axis >= search_range(1) & f_axis <= search_range(2)), ... f_axis(f_axis >= search_range(1) & f_axis <= search_range(2)), ... 'NPeaks', 1, 'SortStr', 'descend'); % 边带搜索：主峰左右±fr处 left_band = f_bpfi - fr; right_band = f_bpfi + fr; band_tolerance = 2; if any(abs(f_axis - left_band) < band_tolerance) && any(abs(f_axis - right_band) < band_tolerance) disp('内圈故障特征明显：发现BPFI主峰及两侧转频边带'); else disp('未检测到完整边带结构，需结合时域冲击间隔再做判断'); end

判断边带存在的代码逻辑谈不上复杂，难点在于噪声导致某一边的边带峰值幅度很低。经验做法是降低阈值只搜索某一侧，比如先找左侧边带，存在左侧再找右侧，同时边带峰值只要高于周围局部噪声底的三倍即可认为有效。单纯机械地找左右两侧峰值会漏掉诊断价值极大的弱边带。

5. 让诊断结果更可靠：谱峰稳定性验证与误报排查

5.1 多段数据循环验证，剔除随机冲击干扰

一次FFT得到的主峰值不一定是故障特征，也可能是某个无关的随机冲击。工业现场的锤击、电磁干扰、负载突变都会在包络谱上留下单次高峰。最稳妥的验证手段是：采集多段独立数据，每段做一次包络谱，比较各段的峰值频率是否稳定出现在同一位置。

% 将长信号切分为多段，每段独立做包络谱 seg_len = fs * 2; % 每段2秒 n_seg = floor(length(data) / seg_len); peak_matrix = zeros(n_seg, 5); % 记录每段的前5个峰 for i = 1:n_seg seg = data((i-1)*seg_len+1 : i*seg_len); [spec, f_axis] = envelope_spectrum(seg, fs, 2000, 6000); [pks, locs] = findpeaks(spec, f_axis, 'SortStr', 'descend', 'NPeaks', 5); peak_matrix(i, 1:length(locs)) = locs; end % 统计每段峰值的标准差 peak_std = std(peak_matrix, 0, 1);

正常运行状态下，各段的峰值频率应当保持一致。如果第1段出现一个高峰但第2、3段相同频率处没有对应能量，那基本可以确定是偶发干扰而非固定故障。另一种角度是看转速变化：如果转速从900rpm升到1200rpm，特征频率理论值也应按比例升高，如果实测峰值不随转速等比移动，说明那个峰不是故障产生的调制分量。

5.2 基线漂移与阈值自动校准问题

固定阈值报警是状态监测系统最常见的误报原因。设备在不同工况下载荷不同、转速不同，振动能量差异巨大。今天定死一个1g的报警阈值，明天换一个工况就天天误报。工程上更合理的是“相对基线法”：设备状态良好时采集一段数据建立基线特征库，包括各倍频幅值、包络谱各特征频率处的幅值分布，之后实时数据的诊断结果与基线对比，超过基线均值加3倍标准差的频点才判异常。

基线自适应的关键在持续更新。我一般会将最近100组正常数据做加权平均作为滑动基线，而不是永远用投产时的那组数据。设备在磨合期、磨损期的正常振动水平不同，固定基线无法适应这种缓变过程。滑动基线加上突变检测，才能在早期故障的“突变”与实际工况造成的“缓变”之间做出区分。

5.3 用spectrogram定位共振频带再定带通滤波参数

前文提到的带通滤波频率范围是包络分析成败的关键，很多时候初判不准带通范围，导致包络谱一片平坦。快速定位共振频带的工具是短时傅里叶变换：

% 振动信号共振频带可视化 [s, f, t] = spectrogram(data, hann(1024), 512, 4096, fs); surf(t, f, 10*log10(abs(s)), 'EdgeColor', 'none'); axis xy; view(0, 90); xlabel('时间 (s)'); ylabel('频率 (Hz)'); colorbar;

从时频图中可以直观看到能量集中带的范围。对滚动轴承，共振频带一般出现在2kHz到8kHz之间，具体取决于结构刚度和质量。能看出明显的水平亮带就是共振频带，将带通滤波器的通带设置在这个范围内。一个被反复验证的技巧是：带通范围不要覆盖整个共振峰，窄带（比共振带宽小20%~30%）往往能获得更干净的包络谱——过宽的带通会把邻近的非调制噪声一起放进来。

时频图看共振带位置的操作每次诊断都值得做，但不需要每次完整跑。建立几种常见结构的共振频率库后，同一型号设备可以直接复用原有带通参数，出现新机型再专门标定一次，现场效率会快很多。

本文还有配套的精品资源，点击获取