MATLAB振动故障诊断：峭度、小波与希尔伯特包络分析的工程实现

简介振动故障诊断是机械设备健康监测与故障预判的核心方法压缩包内含MATLAB信号分析与诊断建模代码适合机械、自动化方向学生、设备维护工程师及相关科研人员学习参考既能用于理解理论原理也能辅助动手实践。压缩包共26个文件容量49KB包含25个.m脚本和1个.txt说明覆盖振动信号时间序列构造、傅里叶级数与FFT频谱分析、小波分解与降噪、希尔伯特变换与包络检波以及概率密度、自协方差、峭度、峰值因子等特征量计算基本构成一套完整的振动特征提取与趋势分析链路。目前已有337人学习。通过运行和修改这些脚本、对照txt说明可直观理解从原始时域波形到频域特征、从统计指标到故障判别的完整过程掌握数据预处理、特征参数计算和诊断结果的工程实现方法同时也能梳理代码运行逻辑适合作为课程实验、毕业设计或科研预研的参考代码库。1. 振动故障诊断的 MATLAB 实现思路倾斜故障为什么最难缠正常运转的设备振动信号是平稳随机过程而轴承座偏斜、转子不对中、轴颈倾斜这类“倾斜故障”在时域波形上跟正常状态几乎看不出差异但峭度、包络谱、小波细节系数会提前数百小时开始漂移。振动故障诊断技术的核心不是看图说话而是把隐藏在噪声里的周期冲击提取出来MATLAB 的 FFT、小波、希尔伯特与统计特征工具箱正好把这条链路变成可运行脚本。这篇博文拆解的 MATLAB 工程包围绕倾斜场景组织qingxieshipin.m 从视频帧提取位移qingxieshijianxulie.m 构造时间序列qingxiexiaobofenxi.m、qingxiehilbert.m、qingxiepianduqiaodu.m 完成小波、希尔伯特与统计特征计算最后用模式识别输出故障类型。适合设备健康监测工程师、写毕业设计的机械或自动化学生以及想把信号处理理论变成代码的研究者还没入门 MATLAB 也不影响阅读每段代码都给了参数注释。2. 时域统计特征计算峭度、偏度、方差与自协方差在 MATLAB 中的实现2.1 为什么先算时域特征倾斜故障的物理表现倾斜故障在物理上表现为转子轴线与轴承座轴线产生夹角运行时油膜厚度不均、轴颈与轴瓦发生局部接触产生周期性冲击。这种冲击在原始波形里只占几个采样点幅值不大却会把概率密度函数的尾部明显拉长。时域统计特征的作用就是把“尾部拉长”量化成峭度、偏度、方差、自协方差等几个数字作为后续诊断模型的输入。对振动工程师来说最直观的经验值是健康滚动轴承的峭度稳定在 3 附近出现剥落、倾斜碰摩时峭度迅速升到 48概率密度曲线出现双峰或明显不对称。偏度描述波形上下不对称倾斜故障常伴随单向摩擦偏度偏移尤其明显。工程包里 qingxiepianduqiaodu.m 对应的就是下面这套计算逻辑。2.2 峭度与偏度pianduqiaodu.m 的计算逻辑与阈值判断function [kurt, skew] pianduqiaodu(x) % 输入 x振动加速度时间序列列向量单位 m/s^2 % 输出 kurt峭度无量纲skew偏度无量纲 x x(:) - mean(x); % 去均值消除传感器零漂 n length(x); sigma std(x, 1); % 总体标准差分母为 n if sigma 0 kurt 0; skew 0; % 常数信号直接返回避免除零 return; end kurt sum(x.^4) / n / sigma^4; % 四阶中心矩 / 标准差^4 skew sum(x.^3) / n / sigma^3; % 三阶中心矩 / 标准差^3 end逻辑说明先做去均值和总体标准差归一化。峭度用四次方放大离群冲击点所以对脉冲类故障最敏感偏度用三次方保留正负号专门捕捉波形上下不对称。std(x,1) 除以 n 与 std(x) 除以 n-1 的区别在长序列上很小但截取短窗计算峭度时推荐统一用 std(x,1)否则峭度基线会有偏移。参数说明窗长取 2048 或 4096 个点约等于 20 倍以上转频周期窗太短峭度抖动剧烈窗太长故障冲击被平均掉。报警阈值不建议用固定值 3 一刀切正确做法是在设备正常阶段统计 100 组峭度取均值加 3 倍标准差作为报警线。特征参数对应函数文件正常参考范围倾斜或碰摩趋势峭度pianduqiaodu.m / qingxiepianduqiaodu.m2.53.5升至 48偏度pianduqiaodu.m-0.30.3明显偏向某一侧方差junzhifangcha.m / qingxiejunzhifangcha.m基线稳定随损伤缓慢上升自协方差峰值zixiefangcha.m / qingxiezixiefangcha.m小于 0.2周期性峰值超过 0.3概率密度熵gailvmidu.m / qingxiegailvmidu.m56 bit冲击成分越多熵越低2.3 方差、自协方差与概率密度估计的实现方差体现振动能量整体水平junzhifangcha.m 就是 var(x) 的手写版不再展开自协方差的作用是找出冲击的重复周期因为故障冲击的间隔会在自协方差序列上形成稳定的峰。工程包里的 zixiefangcha.m 实现如下function [rxx, tau] zixiefangcha(x, maxLag) % 自协方差序列用于检测周期性冲击间隔 x x(:) - mean(x); n length(x); maxLag min(maxLag, n - 1); rxx zeros(maxLag 1, 1); for tau 0:maxLag rxx(tau 1) x(1:n - tau) * x(tau 1:n) / (n - tau); end rxx rxx / rxx(1); % 归一化0 时刻为 1 tau (0:maxLag); end逻辑说明零时刻的自协方差就是方差归一化后得到自相关系数。周期性故障的冲击间隔 T 会在 lagT 处出现一个峰峰值大于 0.2 就可以作为故障周期证据如果峰值均匀衰减则说明信号近似白噪声周期性较弱。概率密度估计对应 gailvmidu.m实际计算可以直接调 MATLAB 的 ksdensity[f, xi] ksdensity(x, Bandwidth, 0.05);带宽参数直接决定曲线平滑程度带宽取 0.05 时能保留双峰细节取 0.5 会把两个峰合并成一个故障特征就被抹平了。现场信号如果转速波动大带宽要适当放大否则画出来的概率密度曲线全是毛刺。2.4 把时域特征组装成模型输入单看一个特征很难下结论工程上的做法是把峭度、偏度、方差、自协方差峰值、概率密度熵拼成特征矩阵再送进模式识别模型。信息熵的计算可以自己写也可以借 MATLAB Statistics 工具箱里的 histcounts 先统计分布再按信息熵公式计算故障时概率密度变宽或出现多峰熵值下降。特征之间量纲差异很大记得做 zscore 归一化否则距离类算法会被方差主导。样本多、特征维度高时可以先用 MATLAB 优化工具箱跑一遍逐步回归或遗传算法特征筛选把不敏感特征丢掉模型训练速度和准确率都会有提升。3. 频域与复频域傅里叶级数系数、希尔伯特变换与包络解调的实现3.1 FFT 的频谱泄漏与窗函数选择很多初学者拿振动波形直接 fft得到的是整段数据的平均频谱对倾斜故障这类非平稳信号会漏掉瞬态冲击。FFT 隐含周期性延拓当窗长不是信号周期的整数倍时频谱发生泄漏故障特征频率被旁瓣淹没。常见做法是在 fft 前加汉宁窗 hann 或平顶窗 flattopwin汉宁窗主瓣窄适合转频附近的边带分析平顶窗幅值精度高适合振动总量校准。加窗后幅值会衰减一半左右回乘系数 2 即可。即使现在能用自然语言让 AI 编程助手直接生成 MATLAB 脚本窗函数和采样率这些概念不搞清楚生成的代码跑出来的谱照样是坏的。采样率方面工程惯例是取最大分析频率的 2.56 倍以上不是教科书的 2 倍因为硬件抗混叠滤波器存在过渡带。比如要看到 500Hz 的故障特征频率采样率建议设 1280Hz 或更高。3.2 fuliyejishuxishu.m傅里叶级数系数与故障特征频率对照傅里叶级数系数把周期信号分解成基频整数倍分量机械上对应 1X、2X、3X 转频成分。倾斜不对中的典型标志是 2X 分量异常升高油膜涡动则表现为 0.38X0.48X 低频分量。工程包里的 fuliyejishuxishu.m 按转频窄带做能量求和function [amp, fr] fuliyejishuxishu(x, Fs, N) % 输入 x时间序列Fs采样率N需要考察的谐波阶数 % 输出 ampN 维谐波幅值向量fr基频估计值 L length(x); f Fs * (0:L-1) / L; % 频率轴 X abs(fft(x - mean(x))) * 2 / L; % 单边幅值谱去直流 X X(1:floor(L/2)1); [~, i0] max(X(2:end)); % 找全局最大峰作为转频 fr f(i0 1); amp zeros(N, 1); for k 1:N band (f k*fr*0.98) (f k*fr*1.02); % 谐波窄带 amp(k) sum(X(band)); % 窄带能量求和 end end逻辑说明直接取单根谱线容易被泄漏和噪声干扰对 1X、2X 各取转频附近 ±2% 的窄带能量求和相当于对频谱做平滑积分工程上更稳定。参数说明转速波动大时把 0.98 和 1.02 放宽到 0.95 和 1.05N 取 3 或 4 就够高阶谐波幅值太低多算只会引入噪声。故障类型特征频率主要表现转子不对中2X部分伴随 1X2X 谐波幅值骤升轴倾斜或局部碰摩0.5X、1X 边带分数倍频成分出现滚动轴承外圈故障BPFO 约 0.4nX共振频带出现包络尖峰油膜涡动0.38X0.48X低频缓慢波动3.3 hilbertmotai.m 与希尔伯特包络解调让故障频率浮现出来希尔伯特变换的作用是构造解析信号实部是原信号虚部是原信号的 90 度相移包络就是解析信号的模。故障冲击激起的高频共振是载波故障特征频率是调制信号直接做 FFT 只能看到共振频带上一大片能量先带通滤波再做希尔伯特包络最后对包络做 FFT故障特征频率才会以清晰谱峰的方式出现。hilbertmotai.m 的完整流程如下function [env_spec, f_env] hilbertmotai(x, Fs, fc, bw) % fc带通中心频率bw带通宽度 d designfilt(bandpassiir, FilterOrder, 4, ... HalfPowerFrequency1, fc-bw/2, HalfPowerFrequency2, fcbw/2, ... SampleRate, Fs); x_filtered filtfilt(d, x); % 零相位滤波避免相位畸变 env abs(hilbert(x_filtered)); % 希尔伯特变换求包络 L length(env); f_env Fs * (0:L/2) / L; env_spec abs(fft(env - mean(env))) * 2 / L; env_spec env_spec(1:floor(L/2)1); end逻辑说明filtfilt 做零相位滤波正向和反向各过一遍包络不会出现相位偏移。中心频率 fc 的选择要基于原始信号功率谱密度的共振峰把 fc 设在峰中心bw 按共振峰宽度取 5002000Hz 都常见太窄会切掉冲击能量太宽会放进无关噪声。3.4 xiangjiapingjun.m 与相角平均同步抑制非周期干扰相角平均也叫时间同步平均利用键相脉冲确定每一转的起点把多转波形按角度对齐后平均。随机噪声和与转频不同步的振动分量在平均中被抵消与转频严格同步的故障分量被保留增强。包里的 qingxiexiangjiapingjun.m 和 xiangjiapingjun.m 就是干这个的function x_avg xiangjiapingjun(x, key_phase, samples_per_rev) % key_phase键相脉冲对应的采样点索引 n length(key_phase) - 1; x_avg zeros(samples_per_rev, 1); for i 1:n seg x(key_phase(i):key_phase(i)samples_per_rev-1); x_avg x_avg seg; end x_avg x_avg / n; % 平均后非同步成分衰减 end参数说明samples_per_rev Fs / 转频取整后作为每转采样点数平均圈数 n 越多噪声按 1/sqrt(n) 衰减现场一般取 32 转以上效果才明显。注意键相脉冲要有足够幅值否则索引抖动会把相位平均的结果再次破坏。4. 小波分析与小波降噪非平稳振动信号的时频局部化处理4.1 为什么倾斜故障的瞬态冲击必须用小波FFT 是全局变换得到的是整段信号的平均频谱没法回答“这个冲击发生在哪一转、哪一刻”。倾斜碰摩、轴承滚道剥落产生的瞬态冲击只持续几个毫秒能量在 FFT 里被平均到整个时间窗峰值明显变矮。小波变换用可变窗口处理这个问题低频处窗口长、频率分辨率高高频处窗口短、时间定位准。cwt 适合可视化时频谱dwt 和 wavedec 适合分解重构与小波降噪。4.2 xiaobofenxi.m连续小波变换的 MATLAB 实现工程包里 qingxiexiaobofenxi.m 和 xiaobofenxi.m 的核心就是一行 cwtfunction [cfs, frq] xiaobofenxi(x, Fs) % 连续小波变换输出时频矩阵 [cfs, frq] cwt(x, amor, Fs); % 复 Morlet 小波兼顾幅值与相位 end逻辑说明cfs 是复小波系数矩阵行对应频率列对应时间画图时用 imagesc(t, frq, abs(cfs)) 显示幅值时频谱相位信息则保留在 angle(cfs) 里可做进一步包络相位分析。小波基的选择原则amor 即解析 Morlet 小波频率分辨率好适合旋转机械的转频和边带分析morse 是通用性强的新默认小波bump 小波频带隔离度极高适合把两个相邻的故障特征频率剥离开但对瞬态冲击的时间定位会模糊。现场先跑一遍 cwt 看全貌再根据频带的分离难度换小波基。4.3 xiaobojiangzao.m小波阈值降噪的参数组合小波降噪的目的是保留冲击成分、抑制平稳噪声不能用默认参数一把梭。需要确定的参数有三组小波基、分解层数、阈值规则。工程包里的 xiaobojiangzao.m 封装成function x_den xiaobojiangzao(x, Fs, wname, level) % wname小波基如 sym5level分解层数经验取 4~6 x_den wdenoise(x, level, Wavelet, wname, ... DenoisingMethod, Bayes, ... ThresholdRule, Median, ... NoiseEstimate, LevelDependent); end逻辑说明Bayes 方法适合现场噪声强度未知的场景噪声接近高斯白噪声时rigrsure 无偏风险估计更准。ThresholdRule 选 Median 比 Mean 更抗离群值干扰。NoiseEstimate 选 LevelDependent 会对每一层单独估计噪声水平比全局估计更贴合真实信号。新版本 MATLAB 的 wdenoise 接口在不同发行版之间有细节变化显式指定 Wavelet 和 ThresholdRule 就可避免默认值不一致的坑。验证降噪效果有个很实用的办法对比降噪前后的峭度。真实冲击成分应保留所以降噪后峭度不应明显下降如果峭度从 6 掉到 2.8说明把小波细节系数里的冲击当成噪声削掉了此时应减小分解层数或改用更保守的阈值规则。小波基特性适用场景sym5对称性好、消失矩 5通用振动降噪db4紧支撑、计算快实时在线监测amor / morse复数解析小波cwt 时频谱、包络形态分析bump频带隔离强分离相邻故障特征频率4.4 xinhaojifen.m 与 njdaoshu.m积分和微分在预处理中的作用加速度信号积分成速度或位移时低频漂移会被积分放大标准流程是小波降噪后再积分同时用高通滤波去掉 0.5Hz 以下的分量。工程包里的 xinhaojifen.m 核心是 cumtrapz 数值积分v cumtrapz(t, x_den); % 梯形法数值积分 v v - polyfit(t, v, 1) * [t; ones(size(t))]; % 去掉线性趋势漂移说明cumtrapz 是梯形法积分前必须先降噪polyfit 拟合一次项并减去消除积分过程累积的直流漂移。njdaoshu.m 对应 n 阶导数计算求导会放大高频噪声先小波降噪再做差分才稳定实际工程中也可以直接把某一层小波细节系数当作导数的估计量效果比直接 diff 平滑得多。5. 从倾斜视频到时间序列完整跑通振动故障诊断链路5.1 qingxieshipin.m从视频帧提取振动位移工程包里的 qingxieshipin.m 读取倾斜圆柱体或转轴视频对每一帧做二值化和质心跟踪输出目标中心像素位移再乘以标定系数换算成真实位移。这是非接触式视觉测振最简实现。普通 30fps 视频只能分析 15Hz 以内的低频晃动想看到几百赫兹的故障冲击需要高速相机。代码思路是逐帧取质心for k 1:numFrames img rgb2gray(read(v, k)); bw imbinarize(img, graythresh(img)); % Otsu 阈值分割 s regionprops(bw, Centroid); pos(k) s(1).Centroid(1); end pos pos - movmean(pos, 15); % 去除手持拍摄的慢漂移灰度阈值用 graythresh 自动估计光照变化大时改用 imbinarize 的 adaptive 模式否则质心会跳动。5.2 qingxieshijianxulie.m重采样与 CSV 导入 MATLAB 做 FFT 验证视频帧率不稳定会导致时间序列非等间隔qingxieshijianxulie.m 用 interp1 做样条重采样得到固定采样率数据后存成 CSV。把 CSV 导入到 MATLAB 中进行 FFT 仿真验证是诊断流程里反复出现的高频操作data readmatrix(vibration_signal.csv); t data(:, 1); x data(:, 2) - mean(data(:, 2)); Fs 1 / median(diff(t)); % 从时间列反推采样率 [pxx, f] pwelch(x, hann(1024), 512, 4096, Fs); [~, i] max(pxx(2:end)); fprintf(主峰值频率: %.2f Hz\n, f(i));pwelch 是 Welch 平均周期图法比直接 FFT 稳定适合现场含噪信号。hann(1024) 是窗长512 是重叠点数4096 是 FFT 点数采样率 1000Hz 时1024 点窗能区分约 1Hz 的相邻频率足够看清 1X、2X 和 0.4X 的差别。算出的主峰频率要和理论转频对比先由转速算 fr_rot再看峰值是否落在 1X、2X 或 0.4X 附近以此判断振动故障诊断技术链路是否跑通。5.3 taylor.m 与 youren.txt非线性特征与启动顺序工程包里的 taylor.m 对非线性刚度段的振动信号做泰勒展开用多项式系数描述非线性程度。倾斜故障导致配合间隙处的恢复力出现明显非线性展开后的三阶项系数变化比一阶项灵敏得多适合作为补充特征。实际使用时先读根目录下 youren.txt 里的运行说明再跑 mainprogram 主脚本生成基准结果然后用自己采集的振动 CSV 替换输入文件逐段执行 qingxie 系列函数。哪个环节的结果异常就回到对应章节检查参数峭度异常看窗长包络谱没峰看中心频率降噪后波形发虚就调分解层数。这套代码不用改架构就能从演示项目变成你自己的设备健康监测脚本。本文还有配套的精品资源点击获取