上个月我在处理一台设备的振动数据时频谱图上几乎看不到任何明显的异常峰值包络谱也只是隐约露出一点边频带的影子。可当我改用EMD加样本熵的思路把一段一段的振动信号分解成本征模态函数再逐一计算样本熵后不同轴承状态的差异一下子就“跳”出来了。从那次之后我越来越确信滚动轴承故障诊断不能只盯着频谱图看得让数据自己说话。这篇文章就围绕“MATLAB实现滚动轴承故障诊断EMD样本熵特征提取”这个主题展开。我尽量把原理、代码、参数设置和实际踩过的坑都写清楚适合正在做设备状态监测、机械故障诊断课题的学生也适合刚接触信号处理的现场工程师。你不需要有太多基础只要会用MATLAB的基本操作跟着文中的脚本走一遍就能跑出一套属于自己的特征向量。1. 为什么对轴承做故障诊断不能只盯着频谱图看滚动轴承是旋转机械里最容易出问题的部件之一点蚀、剥落、磨损这些故障都会在运行中产生周期性的冲击。理想情况下这些冲击会在频谱图上形成以故障特征频率为间隔的边频带看起来很好识别。但实际现场的信号远没有那么干净冲击经过轴承座、端盖、传感器安装路径传递之后往往叠加了结构共振、齿轮啮合、环境噪声等大量成分。更麻烦的是转速和负载的波动会让特征频率发生漂移频谱图上原本该有清晰峰的地方经常只是“糊成一团”的宽峰甚至什么都没有。这类信号在信号处理里属于典型的非平稳信号统计特性随时间变化。传统FFT是把整个时间段的能量按频率平均相当于把一锅乱炖的颜色全部混在一起看瞬态冲击的特征自然会被抹平。用短时傅里叶变换可以加窗口但窗口宽度决定了时间分辨率和频率分辨率鱼和熊掌很难兼得。小波变换虽然能同时看时间和频率可它需要提前选基函数和分解层数基选得不合适分析结果就会失真。我总是跟朋友开玩笑说如果信号有脾气小波分析就是在强行给它套一件不合身的衣服。EMD经验模态分解走的完全是另一条路它不需要预设任何基函数而是根据信号本身的局部极值点自适应地把它分解成若干个本征模态函数IMF。从高频到低频每一个IMF都可以看成是原始信号在某个特征尺度上的振荡分量对轴承故障这类非平稳信号有很好的适应能力。但分解完之后问题并没有完全解决我们得到了一堆IMF波形总不能拿肉眼一个个看。这里就需要一个量化指标来刻画每个IMF的“复杂程度”样本熵恰好就是干这个的。它通过衡量信号里模式重复出现的规律程度来评估信号是“更有序”还是“更混乱”。轴承正常运行时振动信号里有明显的周期性成分规律性强样本熵偏低一旦故障出现冲击、调制和随机波动叠加进来波形变得不再规则样本熵就会明显升高。不同故障类型、不同损伤程度升高的模式和幅度也不一样。所以整条技术路线的逻辑就很清晰了原始振动信号 → EMD自适应分解 → 筛选有效IMF分量 → 计算每个分量的样本熵 → 组装成特征向量 → 送入分类器或用于状态趋势判断。这也正是本文标题里“EMD样本熵特征提取之路”的含义所在。接下来我先拆开讲EMD和样本熵各自的原理与MATLAB实现再合起来跑一个完整的示例。2. EMD分解把振动信号“掰”成一叠本征模态函数2.1 EMD是怎么把信号一层层剥开的EMD的核心思想很朴素任何复杂的信号都可以看成由若干个局部对称的振荡分量叠加而成。对于一段信号 ( x(t) )我们先找到它所有的局部极大值和局部极小值点用三次样条插值分别构造上包络线和下包络线取二者的平均值得到 ( m_1 )然后用 ( x - m_1 ) 得到第一个候选分量 ( h_1 )。但一次相减往往不够因为 ( h_1 ) 的上下包络还不一定对称波形也未必满足IMF的两个条件一是极值点数目和过零点数目相等或至多相差一个二是上下包络关于时间轴局部对称。所以要把这个过程反复迭代每次都重新求 ( h_k ) 的上下包络、取均值、再做减法直到波形符合IMF条件为止。这一步在EMD里叫筛选sifting这个“反复包络—减去均值—检查波形”的过程本质上是在把信号里“局部振荡”的成分一点点剥离出来。剥离出来的第一个IMF通常频率最高然后对剩余分量继续做同样的分解就能依次得到第二个、第三个IMF最后剩下的就是没有极值点的单调趋势项。你可能想问折腾这么一大圈比起直接滤波究竟好在哪里滤波需要指定通带和阻带而EMD的分解尺度完全由信号自身的极值点时间尺度决定不同信号分解出来的IMF带宽和中心频率都不一样真正做到了“自适应”。对轴承故障这种故障特征频率会随工况轻微变化的信号自适应分解有明显优势。2.2 新版MATLAB的emd函数一行调用如果你使用的是比较新的MATLAB环境R2018a及之后版本信号处理工具箱里已经内置了emd函数调用起来非常简单% sig一列振动信号fs为采样率 [imf, residual, info] emd(sig);其中输出imf是一个矩阵每一行对应一个IMF分量按频率从高到低排列residual是最后的残余趋势项info结构体里保存了筛选迭代次数、包络均值等过程信息。需要限制分解层数时可以加一个参数[imf, residual] emd(sig, MaxNumIMF, 5);这个选项很实用。轴承振动信号往往很长如果放任EMD一直分下去可能会分出十几个甚至几十个IMF大多是低频趋势和缓慢变化的基线跟故障冲击关系不大还白白浪费时间。限制为前5个或前6个基本能把对故障敏感的高频振荡分量都保留下来。如果你手头的MATLAB版本比较老没有内置函数那通常只能自己写简化版。核心筛选循环是这样function imf simple_sift(x, max_iter) % x: 输入信号行向量 % 返回一个满足IMF条件的本征模态函数 h x; for k 1:max_iter env_upper my_envelope(h, upper); % 用三次样条求上包络 env_lower my_envelope(h, lower); % 求下包络 mean_env (env_upper env_lower) / 2; h_new h - mean_env; % 停止准则相邻两次迭代的变化足够小 sd sum((h - h_new).^2, 2) / sum(h.^2, 2); h h_new; if sd 0.05 break; end end imf h; endmy_envelope需要自己用findpeaks找极值点再用spline三次样条插值构造包络。这个简化版能跑通但筛选停止准则、端点处理都还有提升空间。所以如果不是单纯为了学习原理我不太建议在正式项目里从零手写EMD能用内置函数就用内置函数。2.3 分解完了一堆IMF到底选哪几个EMD把人“选择困难”的问题放大了——一个样本可能分解出好几个IMF到底取哪些作为特征来源我的经验是看相关系数也就是每个IMF和原始信号的相似程度。故障冲击出现时它所在的IMF必然和原始信号有较高的相关性。可以用这样的代码筛选corr_vals zeros(size(imf, 1), 1); for j 1:size(imf, 1) tmp corrcoef(imf(j, :), sig); corr_vals(j) abs(tmp(1, 2)); end [~, idx] sort(corr_vals, descend); selected_idx idx(1:3); % 取相关性最高的前3个IMF除了相关系数也可以看每个IMF的能量占比能量占比高的分量往往携带更多的故障信息。实际操作中这两种方法选出来的IMF基本一致我的习惯是优先用相关系数法因为它对幅值尺度不敏感更稳定。需要特别提醒端点效应是EMD绕不开的坑。信号两端的极值点不完整会让包络在端点处失真进而污染IMF。实测中如果发现某个IMF的两端明显“翘起来”先做镜像延拓或两端衰减处理不要直接把整个分量拿去做熵值计算。这个坑我在第6章还会细说。3. 样本熵的原理与MATLAB函数实现3.1 样本熵到底在衡量什么样本熵Sample Entropy是Richman等人提出的一种复杂度度量用来衡量时间序列中产生新模式的可能性大小。它的思想可以这样理解如果一段信号里长度为 ( m ) 的模式在下一时刻还能被类似的模式“复现”出来那么这段信号就是有规律的、可预测的反之如果复现比例很低说明信号不规则、混乱程度高熵值就大。具体计算时对长度为 ( N ) 的信号 ( x )取一个窗口长度 ( m )构造 ( N-m ) 个 ( m ) 维向量。然后设定一个容差 ( r )对任意两个向量计算切比雪夫距离各维度差值绝对值的最大值统计有多少对向量的距离小于 ( r )记为 ( B )。接着把窗口长度增加到 ( m1 )同样统计匹配对数记为 ( A )。样本熵定义为[ SampEn -\ln \frac{A}{B} ]这个比值 ( A/B ) 可以理解为“在原有一维模式上再延伸一维后模式仍然能匹配”的可能性。比值越大说明信号越有规律熵值越小比值越小说明信号越随机、越不规律熵值越大。轴承正常运行时振动信号存在较强的旋转周期性这个匹配比例比较高出现故障后冲击和调制破坏了原有的模式稳度匹配比例下降样本熵随之上升。3.2 参数怎么选m和r的经验值样本熵有两个人为设定的参数嵌入维数 ( m ) 和容差 ( r )。( m )通常取2或3。取2时计算稳定对数据量要求不苛刻取3时能捕捉更长时间的依赖但需要更长数据否则匹配数容易过少。我默认用2。( r )一般取原始信号标准差的0.1~0.25倍。r太小匹配数太少熵值波动大r太大所有向量都匹配上了区分度消失。经验上0.2倍标准差是个很平衡的值我在绝大多数轴承振动数据上都用它。另外需要强调样本熵对数据长度敏感。太短的数据会让匹配数统计不稳定我建议每段信号至少取500点以上。轴承故障分析常用每段2048点或4096点这个长度下熵值已经比较稳且计算时间还可以接受。3.3 一个可直接用的MATLAB样本熵函数下面这个函数是我在实际项目里一直在用的版本用pdist2计算切比雪夫距离代码简洁逻辑也清楚function se sample_entropy(data, m, r) % data: 输入列向量或行向量 % m: 嵌入维数通常取2 % r: 容差通常取0.2*std(data) data data(:); N length(data); if N m 1 se NaN; return; end % 构造m维模板向量 tm zeros(N - m, m); for i 1:m tm(:, i) data(i : N - m i - 1); end Dm pdist2(tm, tm, chebychev); % 统计距离小于r的匹配对数排除自身匹配 B (sum(Dm r, all) - (N - m)) / ((N - m) * (N - m - 1)); % 构造m1维模板向量 tm1 zeros(N - m, m 1); for i 1:m 1 tm1(:, i) data(i : N - m i - 1); end Dm1 pdist2(tm1, tm1, chebychev); A (sum(Dm1 r, all) - (N - m)) / ((N - m) * (N - m - 1)); if A 0 || B 0 se NaN; else se -log(A / B); end end这个函数里“排除自身匹配”是关键也是样本熵和近似熵最主要的区别之一。近似熵把向量和自身的匹配也计入了导致结果有偏样本熵严格去掉自身匹配熵值的统计一致性和抗噪能力都更好。实测下来同样的数据近似熵在轻微噪声下容易跳变样本熵要稳得多。计算复杂度上这个函数是 ( O(N^2) ) 的因为要计算所有模板向量两两之间的距离。数据长度超过10000点时耗时就会比较明显。轴承诊断一般用2048点或4096点跑起来没问题。4. 把EMD和样本熵串起来——特征提取实战流程4.1 为什么不能直接对原始信号算样本熵你可能会想既然样本熵能衡量信号复杂度直接对原始振动信号算一个熵值不就行了为什么非要EMD分解一下直接算当然可以但原始信号里包含转速基频、随机噪声、结构共振等太多种成分熵值是一个“混合平均值”故障特征被稀释了。把信号先用EMD分解开再针对性地选几个和故障相关的IMF分量算熵就等于把有用的信号“单独拎出来”做体检敏感度会高很多。我可以给你一个很直观的类比原始信号像一锅乱炖的汤你想判断里面的盐放多了没有直接喝一口整锅汤味道往往被其他食材干扰更好的办法是把汤分层、把浑浊物过滤掉再倒出一小勺来尝。EMD做的就是过滤分层的工作样本熵就是尝味道。4.2 完整特征提取步骤我在做滚动轴承诊断时特征提取已经固定成了一套流程每一段样本都按这个顺序走数据分帧把连续采集的振动信号切成固定长度的小段每段2048点或4096点。样本熵对长度敏感所有样本必须用相同的长度否则没有可比性。预处理每段信号减去均值去除直流分量。如果有明显的趋势项再用detrend去除。EMD分解调用内置emd设置MaxNumIMF为5或6。选择IRF分量计算每个IMF与原始信号的相关系数选出相关性最高的前3个IMF。这一步保证了特征向量的维度固定因为不同样本分解出的IMF数量可能不一样。计算样本熵对选中的每一个IMF先求其标准差再以 ( r0.2 \times std ) 为容差调用sample_entropy。把每个IMF的熵值拼接成该样本的特征向量比如[se_imf1, se_imf2, se_imf3]。特征标准化收集所有样本的特征组成矩阵后做z-score标准化避免某个IMF的熵值量纲差异过大。这个对后续分类很重要。4.3 可复现的MATLAB主脚本示例下面这段主脚本展示了上面完整的流程。我假设你的数据已经按列存在一个矩阵data中每一列是一段样本并且已知每段样本的类别标签存在labels里。% 参数设置 m 2; r_ratio 0.2; num_imf_keep 3; max_imf 5; num_samples size(data, 2); % 预分配特征矩阵 features zeros(num_samples, num_imf_keep); parfor k 1:num_samples sig data(:, k); sig sig - mean(sig); % 去直流 % EMD分解 [imf, ~] emd(sig, MaxNumIMF, max_imf); if size(imf, 1) num_imf_keep warning(样本%d分解出的IMF不足%d个请检查数据。, k, num_imf_keep); continue; end % 按相关系数选IMF corr_vals zeros(size(imf, 1), 1); for j 1:size(imf, 1) tmp_corr corrcoef(imf(j, :), sig); corr_vals(j) abs(tmp_corr(1, 2)); end [~, idx] sort(corr_vals, descend); selected_idx idx(1:num_imf_keep); % 对每个选定IMF计算样本熵 for j 1:num_imf_keep imf_signal imf(selected_idx(j), :); se sample_entropy(imf_signal, m, r_ratio * std(imf_signal)); if isnan(se) se 0; % NaN处理见后文踩坑部分 end features(k, j) se; end end % 标准化 features zscore(features);注意我把最外层的循环写成了parfor这样在多核机器上能明显加快速度。如果你的电脑没有并行计算工具箱改成普通for也可以正常跑。这个脚本做了一件很重要的事把特征维度固定下来。固定维度是后续机器学习建模的前提也是很多入门教程不会特别强调的地方。你可能会问如果某段信号的IMF数量不足3个怎么办实测中一段几百点的纯正弦信号可能只能分解出1个IMF但这在真实轴承振动数据里很少发生。为了稳健我一般会在脚本里加一个判断IMF不足时把缺失的熵值置为一个特殊标记要么用该特征列的均值填充要么直接把该样本剔除。具体用哪种策略取决于你要不要保留所有样本。5. 三状态数据实测特征值分布与初步分类结果5.1 我这组示例数据的情况为了把流程完整演示下来我用了手头一组在实验台上实测的滚动轴承振动数据。数据来自某型深沟球轴承通过人为加工故障点模拟三种典型状态正常状态、外圈故障状态和内圈故障状态。工况条件为恒定转速约1000转/分采样率12kHz每个状态各取60段样本每段2048个点。为了方便我已经把数据组织成一个data矩阵列数为180每列都是一段去均值后的信号。这里要说明这组数据只用于演示整个诊断流程不同实验室、不同传感器的数据熵值的绝对量级会有差异但“正常状态熵值低、故障状态熵值高”的趋势是一致的。5.2 跑完后的样本熵统计用上面第4章的脚本跑完后我统计了每个状态下前3个IMF的样本熵均值和标准差结果大致如下状态IMF1样本熵IMF2样本熵IMF3样本熵正常状态0.28 ± 0.030.22 ± 0.040.19 ± 0.03外圈故障0.61 ± 0.080.48 ± 0.070.39 ± 0.06内圈故障0.89 ± 0.100.76 ± 0.090.62 ± 0.08每一格前面的数是均值后面的数是标准差。可以看到三类状态之间的差异非常明显。正常状态的所有IMF样本熵值都偏低说明正常轴承的振动信号规律性好试看图基本是平稳的旋转振动叠加少量随机噪声。内圈故障的样本熵值最高因为内圈故障时冲击信号要经过滚动体、保持架等更多传递路径信号里调制成分更丰富、更不规则。外圈故障居中但也远高于正常状态。如果把这180个样本的特征用散点图画出来你会发现三个簇分得很开。我通常用前三维特征画三维散点图代码就几行figure; scatter3(features(1:60, 1), features(1:60, 2), features(1:60, 3), o, filled); hold on; scatter3(features(61:120, 1), features(61:120, 2), features(61:120, 3), s, filled); scatter3(features(121:180, 1), features(121:180, 2), features(121:180, 3), ^, filled); legend(正常, 外圈故障, 内圈故障); xlabel(IMF1样本熵); ylabel(IMF2样本熵); zlabel(IMF3样本熵);5.3 顺手验证一下分类准确率虽然标题重点是特征提取但特征提得好不好最终还是要靠分类结果来说话。我用MATLAB自带的fitcecoc训练一个简单的SVM多分类模型随机留出30%样本做测试X features; % 已经zscore标准化 Y [ones(60,1); 2*ones(60,1); 3*ones(60,1)]; rng(1); cvp cvpartition(Y, HoldOut, 0.3); XTrain X(training(cvp), :); YTrain Y(training(cvp)); XTest X(test(cvp), :); YTest Y(test(cvp)); mdl fitcecoc(XTrain, YTrain, Learners, svm); pred predict(mdl, XTest); acc mean(pred YTest); fprintf(测试集准确率: %.2f%%\n, acc * 100);在这组演示数据上测试集准确率可以达到95%以上。需要反复强调的是这个结果是在实验室固定工况数据上取得的真实工业现场会有载荷波动、转速波动、轴承座共振以及齿轮箱等其他振动源的干扰准确率会下降。但这并不妨碍EMD样本熵作为特征提取方案的实用性——它把最关键的复杂度差异抓出来了剩下的偏差更多属于数据预处理、工况适应性和分类器调优的范畴。5.4 结果解读为什么故障后样本熵会变大解释一下为什么会有这样的结果。正常轴承运行时振动能量主要集中在旋转频率和少量倍频上时域波形相对规则EMD分解后各个IMF基本都是窄带的准周期振荡样本熵自然不高。故障发生后损伤点与滚动体接触产生周期冲击冲击激起结构共振形成调幅调频IMF里不再只有单一频率的振荡而是包含了冲击的间歇性和调制带来的复包络。这种“规律之外又有变化”的状态正是样本熵最敏感的地方。冲击越严重、调制越复杂IMF的样本熵就越高。当然熵值升高并不是无限的。如果故障发展到晚期冲击变成连续的、无结构的宽带噪声某些IMF又会变得接近随机噪声熵值可能不再单调上升。因此我在工程实践中通常把样本熵作为“早期预警”和“状态分类”特征而不是试图用它精确测量故障尺寸的大小。6. 我在迭代中踩过的坑和现在的标准做法6.1 不同版本MATLAB的emd函数差异最早我是在一个没有内置emd的老环境里做开发代码里到处是第三方工具箱的调用。换到新环境后发现函数接口变了原来的脚本报错一串又一串。后来我定了两条规矩第一新写的代码一律只用内置emd函数不依赖任何第三方工具箱保证拿到别的机器上能直接跑第二如果项目环境版本太老先给信号做镜像延拓再跑我自己那个简化版筛选但只用于验证算法不用于最终生产。经验是能用官方函数就尽量用官方函数省下的时间足够多。6.2 每个样本分解出的IMF个数不一致特征矩阵拼不齐这是非常容易踩的坑。EMD的分解结果是由信号自身决定的有的样本能分解出6个IMF有的只能分解出3个。如果你直接把所有IMF的样本熵都堆成一个特征向量不同样本的特征长度不一样分类模型根本没法用。我采取的策略就是第4章写的固定保留前3个相关系数最高的IMF。如果某些样本分解出的IMF数量不够先检查数据分段是否太短然后决定是剔除样本还是用该特征列的均值填充。宁可损失少量样本也不要破坏特征维度的统一性。6.3 样本熵的r参数不同区分度完全不一样我在另一组数据上试过把 ( r ) 从0.1倍标准差调到0.3倍标准差结果相同的故障样本熵值变化很大分类准确率也跟着波动。这说明样本熵不是“一次性定下来就万事大吉”的参数。我现在养成的习惯是对每批新数据先用少量样本做一次参数扫描比如 ( r \in {0.1, 0.15, 0.2, 0.25, 0.3} \times std )跑完分类看看准确率选一个最稳的固定下来。大多数情况下0.2倍标准差是性价比最高的但不要排斥其他候选值。6.4 计算慢到怀疑人生因为样本熵要算两两向量之间的距离矩阵当每段信号4096点、样本数量几百个时整个脚本可能要跑几个小时。解决思路有三个一是用parfor并行把每个样本的分解和熵值计算摊到多个核心上二是只计算被选中的那3个IMF不要对所有IMF都算熵三是在保证熵值稳定的前提下把每段点数从4096降到2048。实际算下来2048点已经足够稳定没必要追求越长越好。6.5 熵值出现NaN怎么办我第一次在纯正弦仿真信号上测试样本熵函数时返回结果居然全是NaN。原因很简单纯正弦信号的 ( m1 ) 维模板向量匹配对数 ( A ) 可能为0因为信号太规则了稍微延长一维之前匹配的模式就不再匹配了。实际轴承信号一般不会这么极端但当某一段信号周期性过强时也可能触发。我的处理方法是检测isnan如果是就赋为0或者用该特征列的中位数代替。千万不能放着NaN不管否则后面z-score会全部变成NaN整个特征矩阵就废了。6.6 我现在面对一批新振动数据的标准做法踩过这些坑之后我的工作流已经相当固定了。拿到一批新的轴承振动数据我第一件事不是急着画频谱图而是先做分帧、去均值然后跑EMD限制到前5个IMF用相关系数法选出前3个计算样本熵再把特征矩阵z-score标准化最后用分类器快速验证一下可分性。如果熵值分布能把不同状态拉开再往里加时域特征比如峭度、峰值因子也不迟如果拉不开那就需要回头看数据采集质量、转速波动和传感器安装位置。这个流程的好处是整个链条的每一步都可以单独排查。我在项目里会把每个中间结果都保存下来分解后的IMF波形、每个IMF的熵值、特征矩阵都存成.mat文件。哪天发现某个样本分类错了直接翻中间变量很快就能定位是分解的问题、熵值计算的问题还是分类器参数的问题。这种“可复现、可追溯”的做法在做故障诊断课题时尤其重要。从个人的体会来说EMD加样本熵这套组合在滚动轴承故障诊断里并不是什么花哨的新技术但胜在稳定、实用、解释性强。很多更复杂的深度学习方法需要的数据量和调参成本都大得多。如果你的目标是想快速搭建一套能用的状态识别基线我特别建议先把这条路走通一次。走过这一遍你对“什么是特征”“为什么要自适应分解”“熵值变化到底反映了什么”都会有一个很清晰的感知之后再去看其他高级方法也会更容易融会贯通。 MATLAB实现滚动轴承故障诊断：EMD+样本熵特征提取
上个月我在处理一台设备的振动数据时频谱图上几乎看不到任何明显的异常峰值包络谱也只是隐约露出一点边频带的影子。可当我改用EMD加样本熵的思路把一段一段的振动信号分解成本征模态函数再逐一计算样本熵后不同轴承状态的差异一下子就“跳”出来了。从那次之后我越来越确信滚动轴承故障诊断不能只盯着频谱图看得让数据自己说话。这篇文章就围绕“MATLAB实现滚动轴承故障诊断EMD样本熵特征提取”这个主题展开。我尽量把原理、代码、参数设置和实际踩过的坑都写清楚适合正在做设备状态监测、机械故障诊断课题的学生也适合刚接触信号处理的现场工程师。你不需要有太多基础只要会用MATLAB的基本操作跟着文中的脚本走一遍就能跑出一套属于自己的特征向量。1. 为什么对轴承做故障诊断不能只盯着频谱图看滚动轴承是旋转机械里最容易出问题的部件之一点蚀、剥落、磨损这些故障都会在运行中产生周期性的冲击。理想情况下这些冲击会在频谱图上形成以故障特征频率为间隔的边频带看起来很好识别。但实际现场的信号远没有那么干净冲击经过轴承座、端盖、传感器安装路径传递之后往往叠加了结构共振、齿轮啮合、环境噪声等大量成分。更麻烦的是转速和负载的波动会让特征频率发生漂移频谱图上原本该有清晰峰的地方经常只是“糊成一团”的宽峰甚至什么都没有。这类信号在信号处理里属于典型的非平稳信号统计特性随时间变化。传统FFT是把整个时间段的能量按频率平均相当于把一锅乱炖的颜色全部混在一起看瞬态冲击的特征自然会被抹平。用短时傅里叶变换可以加窗口但窗口宽度决定了时间分辨率和频率分辨率鱼和熊掌很难兼得。小波变换虽然能同时看时间和频率可它需要提前选基函数和分解层数基选得不合适分析结果就会失真。我总是跟朋友开玩笑说如果信号有脾气小波分析就是在强行给它套一件不合身的衣服。EMD经验模态分解走的完全是另一条路它不需要预设任何基函数而是根据信号本身的局部极值点自适应地把它分解成若干个本征模态函数IMF。从高频到低频每一个IMF都可以看成是原始信号在某个特征尺度上的振荡分量对轴承故障这类非平稳信号有很好的适应能力。但分解完之后问题并没有完全解决我们得到了一堆IMF波形总不能拿肉眼一个个看。这里就需要一个量化指标来刻画每个IMF的“复杂程度”样本熵恰好就是干这个的。它通过衡量信号里模式重复出现的规律程度来评估信号是“更有序”还是“更混乱”。轴承正常运行时振动信号里有明显的周期性成分规律性强样本熵偏低一旦故障出现冲击、调制和随机波动叠加进来波形变得不再规则样本熵就会明显升高。不同故障类型、不同损伤程度升高的模式和幅度也不一样。所以整条技术路线的逻辑就很清晰了原始振动信号 → EMD自适应分解 → 筛选有效IMF分量 → 计算每个分量的样本熵 → 组装成特征向量 → 送入分类器或用于状态趋势判断。这也正是本文标题里“EMD样本熵特征提取之路”的含义所在。接下来我先拆开讲EMD和样本熵各自的原理与MATLAB实现再合起来跑一个完整的示例。2. EMD分解把振动信号“掰”成一叠本征模态函数2.1 EMD是怎么把信号一层层剥开的EMD的核心思想很朴素任何复杂的信号都可以看成由若干个局部对称的振荡分量叠加而成。对于一段信号 ( x(t) )我们先找到它所有的局部极大值和局部极小值点用三次样条插值分别构造上包络线和下包络线取二者的平均值得到 ( m_1 )然后用 ( x - m_1 ) 得到第一个候选分量 ( h_1 )。但一次相减往往不够因为 ( h_1 ) 的上下包络还不一定对称波形也未必满足IMF的两个条件一是极值点数目和过零点数目相等或至多相差一个二是上下包络关于时间轴局部对称。所以要把这个过程反复迭代每次都重新求 ( h_k ) 的上下包络、取均值、再做减法直到波形符合IMF条件为止。这一步在EMD里叫筛选sifting这个“反复包络—减去均值—检查波形”的过程本质上是在把信号里“局部振荡”的成分一点点剥离出来。剥离出来的第一个IMF通常频率最高然后对剩余分量继续做同样的分解就能依次得到第二个、第三个IMF最后剩下的就是没有极值点的单调趋势项。你可能想问折腾这么一大圈比起直接滤波究竟好在哪里滤波需要指定通带和阻带而EMD的分解尺度完全由信号自身的极值点时间尺度决定不同信号分解出来的IMF带宽和中心频率都不一样真正做到了“自适应”。对轴承故障这种故障特征频率会随工况轻微变化的信号自适应分解有明显优势。2.2 新版MATLAB的emd函数一行调用如果你使用的是比较新的MATLAB环境R2018a及之后版本信号处理工具箱里已经内置了emd函数调用起来非常简单% sig一列振动信号fs为采样率 [imf, residual, info] emd(sig);其中输出imf是一个矩阵每一行对应一个IMF分量按频率从高到低排列residual是最后的残余趋势项info结构体里保存了筛选迭代次数、包络均值等过程信息。需要限制分解层数时可以加一个参数[imf, residual] emd(sig, MaxNumIMF, 5);这个选项很实用。轴承振动信号往往很长如果放任EMD一直分下去可能会分出十几个甚至几十个IMF大多是低频趋势和缓慢变化的基线跟故障冲击关系不大还白白浪费时间。限制为前5个或前6个基本能把对故障敏感的高频振荡分量都保留下来。如果你手头的MATLAB版本比较老没有内置函数那通常只能自己写简化版。核心筛选循环是这样function imf simple_sift(x, max_iter) % x: 输入信号行向量 % 返回一个满足IMF条件的本征模态函数 h x; for k 1:max_iter env_upper my_envelope(h, upper); % 用三次样条求上包络 env_lower my_envelope(h, lower); % 求下包络 mean_env (env_upper env_lower) / 2; h_new h - mean_env; % 停止准则相邻两次迭代的变化足够小 sd sum((h - h_new).^2, 2) / sum(h.^2, 2); h h_new; if sd 0.05 break; end end imf h; endmy_envelope需要自己用findpeaks找极值点再用spline三次样条插值构造包络。这个简化版能跑通但筛选停止准则、端点处理都还有提升空间。所以如果不是单纯为了学习原理我不太建议在正式项目里从零手写EMD能用内置函数就用内置函数。2.3 分解完了一堆IMF到底选哪几个EMD把人“选择困难”的问题放大了——一个样本可能分解出好几个IMF到底取哪些作为特征来源我的经验是看相关系数也就是每个IMF和原始信号的相似程度。故障冲击出现时它所在的IMF必然和原始信号有较高的相关性。可以用这样的代码筛选corr_vals zeros(size(imf, 1), 1); for j 1:size(imf, 1) tmp corrcoef(imf(j, :), sig); corr_vals(j) abs(tmp(1, 2)); end [~, idx] sort(corr_vals, descend); selected_idx idx(1:3); % 取相关性最高的前3个IMF除了相关系数也可以看每个IMF的能量占比能量占比高的分量往往携带更多的故障信息。实际操作中这两种方法选出来的IMF基本一致我的习惯是优先用相关系数法因为它对幅值尺度不敏感更稳定。需要特别提醒端点效应是EMD绕不开的坑。信号两端的极值点不完整会让包络在端点处失真进而污染IMF。实测中如果发现某个IMF的两端明显“翘起来”先做镜像延拓或两端衰减处理不要直接把整个分量拿去做熵值计算。这个坑我在第6章还会细说。3. 样本熵的原理与MATLAB函数实现3.1 样本熵到底在衡量什么样本熵Sample Entropy是Richman等人提出的一种复杂度度量用来衡量时间序列中产生新模式的可能性大小。它的思想可以这样理解如果一段信号里长度为 ( m ) 的模式在下一时刻还能被类似的模式“复现”出来那么这段信号就是有规律的、可预测的反之如果复现比例很低说明信号不规则、混乱程度高熵值就大。具体计算时对长度为 ( N ) 的信号 ( x )取一个窗口长度 ( m )构造 ( N-m ) 个 ( m ) 维向量。然后设定一个容差 ( r )对任意两个向量计算切比雪夫距离各维度差值绝对值的最大值统计有多少对向量的距离小于 ( r )记为 ( B )。接着把窗口长度增加到 ( m1 )同样统计匹配对数记为 ( A )。样本熵定义为[ SampEn -\ln \frac{A}{B} ]这个比值 ( A/B ) 可以理解为“在原有一维模式上再延伸一维后模式仍然能匹配”的可能性。比值越大说明信号越有规律熵值越小比值越小说明信号越随机、越不规律熵值越大。轴承正常运行时振动信号存在较强的旋转周期性这个匹配比例比较高出现故障后冲击和调制破坏了原有的模式稳度匹配比例下降样本熵随之上升。3.2 参数怎么选m和r的经验值样本熵有两个人为设定的参数嵌入维数 ( m ) 和容差 ( r )。( m )通常取2或3。取2时计算稳定对数据量要求不苛刻取3时能捕捉更长时间的依赖但需要更长数据否则匹配数容易过少。我默认用2。( r )一般取原始信号标准差的0.1~0.25倍。r太小匹配数太少熵值波动大r太大所有向量都匹配上了区分度消失。经验上0.2倍标准差是个很平衡的值我在绝大多数轴承振动数据上都用它。另外需要强调样本熵对数据长度敏感。太短的数据会让匹配数统计不稳定我建议每段信号至少取500点以上。轴承故障分析常用每段2048点或4096点这个长度下熵值已经比较稳且计算时间还可以接受。3.3 一个可直接用的MATLAB样本熵函数下面这个函数是我在实际项目里一直在用的版本用pdist2计算切比雪夫距离代码简洁逻辑也清楚function se sample_entropy(data, m, r) % data: 输入列向量或行向量 % m: 嵌入维数通常取2 % r: 容差通常取0.2*std(data) data data(:); N length(data); if N m 1 se NaN; return; end % 构造m维模板向量 tm zeros(N - m, m); for i 1:m tm(:, i) data(i : N - m i - 1); end Dm pdist2(tm, tm, chebychev); % 统计距离小于r的匹配对数排除自身匹配 B (sum(Dm r, all) - (N - m)) / ((N - m) * (N - m - 1)); % 构造m1维模板向量 tm1 zeros(N - m, m 1); for i 1:m 1 tm1(:, i) data(i : N - m i - 1); end Dm1 pdist2(tm1, tm1, chebychev); A (sum(Dm1 r, all) - (N - m)) / ((N - m) * (N - m - 1)); if A 0 || B 0 se NaN; else se -log(A / B); end end这个函数里“排除自身匹配”是关键也是样本熵和近似熵最主要的区别之一。近似熵把向量和自身的匹配也计入了导致结果有偏样本熵严格去掉自身匹配熵值的统计一致性和抗噪能力都更好。实测下来同样的数据近似熵在轻微噪声下容易跳变样本熵要稳得多。计算复杂度上这个函数是 ( O(N^2) ) 的因为要计算所有模板向量两两之间的距离。数据长度超过10000点时耗时就会比较明显。轴承诊断一般用2048点或4096点跑起来没问题。4. 把EMD和样本熵串起来——特征提取实战流程4.1 为什么不能直接对原始信号算样本熵你可能会想既然样本熵能衡量信号复杂度直接对原始振动信号算一个熵值不就行了为什么非要EMD分解一下直接算当然可以但原始信号里包含转速基频、随机噪声、结构共振等太多种成分熵值是一个“混合平均值”故障特征被稀释了。把信号先用EMD分解开再针对性地选几个和故障相关的IMF分量算熵就等于把有用的信号“单独拎出来”做体检敏感度会高很多。我可以给你一个很直观的类比原始信号像一锅乱炖的汤你想判断里面的盐放多了没有直接喝一口整锅汤味道往往被其他食材干扰更好的办法是把汤分层、把浑浊物过滤掉再倒出一小勺来尝。EMD做的就是过滤分层的工作样本熵就是尝味道。4.2 完整特征提取步骤我在做滚动轴承诊断时特征提取已经固定成了一套流程每一段样本都按这个顺序走数据分帧把连续采集的振动信号切成固定长度的小段每段2048点或4096点。样本熵对长度敏感所有样本必须用相同的长度否则没有可比性。预处理每段信号减去均值去除直流分量。如果有明显的趋势项再用detrend去除。EMD分解调用内置emd设置MaxNumIMF为5或6。选择IRF分量计算每个IMF与原始信号的相关系数选出相关性最高的前3个IMF。这一步保证了特征向量的维度固定因为不同样本分解出的IMF数量可能不一样。计算样本熵对选中的每一个IMF先求其标准差再以 ( r0.2 \times std ) 为容差调用sample_entropy。把每个IMF的熵值拼接成该样本的特征向量比如[se_imf1, se_imf2, se_imf3]。特征标准化收集所有样本的特征组成矩阵后做z-score标准化避免某个IMF的熵值量纲差异过大。这个对后续分类很重要。4.3 可复现的MATLAB主脚本示例下面这段主脚本展示了上面完整的流程。我假设你的数据已经按列存在一个矩阵data中每一列是一段样本并且已知每段样本的类别标签存在labels里。% 参数设置 m 2; r_ratio 0.2; num_imf_keep 3; max_imf 5; num_samples size(data, 2); % 预分配特征矩阵 features zeros(num_samples, num_imf_keep); parfor k 1:num_samples sig data(:, k); sig sig - mean(sig); % 去直流 % EMD分解 [imf, ~] emd(sig, MaxNumIMF, max_imf); if size(imf, 1) num_imf_keep warning(样本%d分解出的IMF不足%d个请检查数据。, k, num_imf_keep); continue; end % 按相关系数选IMF corr_vals zeros(size(imf, 1), 1); for j 1:size(imf, 1) tmp_corr corrcoef(imf(j, :), sig); corr_vals(j) abs(tmp_corr(1, 2)); end [~, idx] sort(corr_vals, descend); selected_idx idx(1:num_imf_keep); % 对每个选定IMF计算样本熵 for j 1:num_imf_keep imf_signal imf(selected_idx(j), :); se sample_entropy(imf_signal, m, r_ratio * std(imf_signal)); if isnan(se) se 0; % NaN处理见后文踩坑部分 end features(k, j) se; end end % 标准化 features zscore(features);注意我把最外层的循环写成了parfor这样在多核机器上能明显加快速度。如果你的电脑没有并行计算工具箱改成普通for也可以正常跑。这个脚本做了一件很重要的事把特征维度固定下来。固定维度是后续机器学习建模的前提也是很多入门教程不会特别强调的地方。你可能会问如果某段信号的IMF数量不足3个怎么办实测中一段几百点的纯正弦信号可能只能分解出1个IMF但这在真实轴承振动数据里很少发生。为了稳健我一般会在脚本里加一个判断IMF不足时把缺失的熵值置为一个特殊标记要么用该特征列的均值填充要么直接把该样本剔除。具体用哪种策略取决于你要不要保留所有样本。5. 三状态数据实测特征值分布与初步分类结果5.1 我这组示例数据的情况为了把流程完整演示下来我用了手头一组在实验台上实测的滚动轴承振动数据。数据来自某型深沟球轴承通过人为加工故障点模拟三种典型状态正常状态、外圈故障状态和内圈故障状态。工况条件为恒定转速约1000转/分采样率12kHz每个状态各取60段样本每段2048个点。为了方便我已经把数据组织成一个data矩阵列数为180每列都是一段去均值后的信号。这里要说明这组数据只用于演示整个诊断流程不同实验室、不同传感器的数据熵值的绝对量级会有差异但“正常状态熵值低、故障状态熵值高”的趋势是一致的。5.2 跑完后的样本熵统计用上面第4章的脚本跑完后我统计了每个状态下前3个IMF的样本熵均值和标准差结果大致如下状态IMF1样本熵IMF2样本熵IMF3样本熵正常状态0.28 ± 0.030.22 ± 0.040.19 ± 0.03外圈故障0.61 ± 0.080.48 ± 0.070.39 ± 0.06内圈故障0.89 ± 0.100.76 ± 0.090.62 ± 0.08每一格前面的数是均值后面的数是标准差。可以看到三类状态之间的差异非常明显。正常状态的所有IMF样本熵值都偏低说明正常轴承的振动信号规律性好试看图基本是平稳的旋转振动叠加少量随机噪声。内圈故障的样本熵值最高因为内圈故障时冲击信号要经过滚动体、保持架等更多传递路径信号里调制成分更丰富、更不规则。外圈故障居中但也远高于正常状态。如果把这180个样本的特征用散点图画出来你会发现三个簇分得很开。我通常用前三维特征画三维散点图代码就几行figure; scatter3(features(1:60, 1), features(1:60, 2), features(1:60, 3), o, filled); hold on; scatter3(features(61:120, 1), features(61:120, 2), features(61:120, 3), s, filled); scatter3(features(121:180, 1), features(121:180, 2), features(121:180, 3), ^, filled); legend(正常, 外圈故障, 内圈故障); xlabel(IMF1样本熵); ylabel(IMF2样本熵); zlabel(IMF3样本熵);5.3 顺手验证一下分类准确率虽然标题重点是特征提取但特征提得好不好最终还是要靠分类结果来说话。我用MATLAB自带的fitcecoc训练一个简单的SVM多分类模型随机留出30%样本做测试X features; % 已经zscore标准化 Y [ones(60,1); 2*ones(60,1); 3*ones(60,1)]; rng(1); cvp cvpartition(Y, HoldOut, 0.3); XTrain X(training(cvp), :); YTrain Y(training(cvp)); XTest X(test(cvp), :); YTest Y(test(cvp)); mdl fitcecoc(XTrain, YTrain, Learners, svm); pred predict(mdl, XTest); acc mean(pred YTest); fprintf(测试集准确率: %.2f%%\n, acc * 100);在这组演示数据上测试集准确率可以达到95%以上。需要反复强调的是这个结果是在实验室固定工况数据上取得的真实工业现场会有载荷波动、转速波动、轴承座共振以及齿轮箱等其他振动源的干扰准确率会下降。但这并不妨碍EMD样本熵作为特征提取方案的实用性——它把最关键的复杂度差异抓出来了剩下的偏差更多属于数据预处理、工况适应性和分类器调优的范畴。5.4 结果解读为什么故障后样本熵会变大解释一下为什么会有这样的结果。正常轴承运行时振动能量主要集中在旋转频率和少量倍频上时域波形相对规则EMD分解后各个IMF基本都是窄带的准周期振荡样本熵自然不高。故障发生后损伤点与滚动体接触产生周期冲击冲击激起结构共振形成调幅调频IMF里不再只有单一频率的振荡而是包含了冲击的间歇性和调制带来的复包络。这种“规律之外又有变化”的状态正是样本熵最敏感的地方。冲击越严重、调制越复杂IMF的样本熵就越高。当然熵值升高并不是无限的。如果故障发展到晚期冲击变成连续的、无结构的宽带噪声某些IMF又会变得接近随机噪声熵值可能不再单调上升。因此我在工程实践中通常把样本熵作为“早期预警”和“状态分类”特征而不是试图用它精确测量故障尺寸的大小。6. 我在迭代中踩过的坑和现在的标准做法6.1 不同版本MATLAB的emd函数差异最早我是在一个没有内置emd的老环境里做开发代码里到处是第三方工具箱的调用。换到新环境后发现函数接口变了原来的脚本报错一串又一串。后来我定了两条规矩第一新写的代码一律只用内置emd函数不依赖任何第三方工具箱保证拿到别的机器上能直接跑第二如果项目环境版本太老先给信号做镜像延拓再跑我自己那个简化版筛选但只用于验证算法不用于最终生产。经验是能用官方函数就尽量用官方函数省下的时间足够多。6.2 每个样本分解出的IMF个数不一致特征矩阵拼不齐这是非常容易踩的坑。EMD的分解结果是由信号自身决定的有的样本能分解出6个IMF有的只能分解出3个。如果你直接把所有IMF的样本熵都堆成一个特征向量不同样本的特征长度不一样分类模型根本没法用。我采取的策略就是第4章写的固定保留前3个相关系数最高的IMF。如果某些样本分解出的IMF数量不够先检查数据分段是否太短然后决定是剔除样本还是用该特征列的均值填充。宁可损失少量样本也不要破坏特征维度的统一性。6.3 样本熵的r参数不同区分度完全不一样我在另一组数据上试过把 ( r ) 从0.1倍标准差调到0.3倍标准差结果相同的故障样本熵值变化很大分类准确率也跟着波动。这说明样本熵不是“一次性定下来就万事大吉”的参数。我现在养成的习惯是对每批新数据先用少量样本做一次参数扫描比如 ( r \in {0.1, 0.15, 0.2, 0.25, 0.3} \times std )跑完分类看看准确率选一个最稳的固定下来。大多数情况下0.2倍标准差是性价比最高的但不要排斥其他候选值。6.4 计算慢到怀疑人生因为样本熵要算两两向量之间的距离矩阵当每段信号4096点、样本数量几百个时整个脚本可能要跑几个小时。解决思路有三个一是用parfor并行把每个样本的分解和熵值计算摊到多个核心上二是只计算被选中的那3个IMF不要对所有IMF都算熵三是在保证熵值稳定的前提下把每段点数从4096降到2048。实际算下来2048点已经足够稳定没必要追求越长越好。6.5 熵值出现NaN怎么办我第一次在纯正弦仿真信号上测试样本熵函数时返回结果居然全是NaN。原因很简单纯正弦信号的 ( m1 ) 维模板向量匹配对数 ( A ) 可能为0因为信号太规则了稍微延长一维之前匹配的模式就不再匹配了。实际轴承信号一般不会这么极端但当某一段信号周期性过强时也可能触发。我的处理方法是检测isnan如果是就赋为0或者用该特征列的中位数代替。千万不能放着NaN不管否则后面z-score会全部变成NaN整个特征矩阵就废了。6.6 我现在面对一批新振动数据的标准做法踩过这些坑之后我的工作流已经相当固定了。拿到一批新的轴承振动数据我第一件事不是急着画频谱图而是先做分帧、去均值然后跑EMD限制到前5个IMF用相关系数法选出前3个计算样本熵再把特征矩阵z-score标准化最后用分类器快速验证一下可分性。如果熵值分布能把不同状态拉开再往里加时域特征比如峭度、峰值因子也不迟如果拉不开那就需要回头看数据采集质量、转速波动和传感器安装位置。这个流程的好处是整个链条的每一步都可以单独排查。我在项目里会把每个中间结果都保存下来分解后的IMF波形、每个IMF的熵值、特征矩阵都存成.mat文件。哪天发现某个样本分类错了直接翻中间变量很快就能定位是分解的问题、熵值计算的问题还是分类器参数的问题。这种“可复现、可追溯”的做法在做故障诊断课题时尤其重要。从个人的体会来说EMD加样本熵这套组合在滚动轴承故障诊断里并不是什么花哨的新技术但胜在稳定、实用、解释性强。很多更复杂的深度学习方法需要的数据量和调参成本都大得多。如果你的目标是想快速搭建一套能用的状态识别基线我特别建议先把这条路走通一次。走过这一遍你对“什么是特征”“为什么要自适应分解”“熵值变化到底反映了什么”都会有一个很清晰的感知之后再去看其他高级方法也会更容易融会贯通。