在旋转机械振动诊断的现场最磨人的往往不是设备真的停了而是故障已经发生你却从振动数据里看不出来。滚动轴承早期损伤产生的冲击非常微弱被转频谐波、齿轮啮合成分和背景噪声层层盖住直接拿原始信号做FFT或者包络谱故障特征频率常常只是一根勉强冒出噪声底部的细线稍不留神就被忽略。这时候先做信号预处理再谈诊断几乎是唯一靠谱的路线。我常用的预处理方案是AR自回归模型加MED最小熵反卷积两级处理先用AR把确定性周期成分从信号里剥离再用MED把残差里的稀疏冲击成分放大。整个链路在Matlab里实现非常顺代码量不大效果却很直观。接下来我会把原理、可复现的脚本和踩坑经验完整讲一遍适合正在做轴承故障诊断、信号预处理或者准备复现相关论文方法的同学。1. 轴承振动信号里到底混着什么AR和MED各管哪一块1.1 轴承故障信号可以拆成三部分在处理任何信号之前先得认清手里拿到的到底是什么。现场采集到的加速度振动信号几乎都可以简化成三部分叠加确定性周期分量包括转轴转频及其各次谐波、齿轮啮合频率和边带。它们严格与转轴旋转同步时域里是正弦波及其变体频域里是一根根离散谱线能量通常很大。故障冲击分量轴承滚道、滚动体或保持架出现局部损伤后每个滚子滚过损伤点都会产生一个极短的瞬态冲击。这个冲击会激起轴承结构和传感器安装谐振表现为高频衰减振荡波形。损伤位置固定时冲击近似周期重复其重复周期就是我们要找的故障特征。随机噪声分量环境振动、测量电路噪声、其他非周期干扰近似高斯白噪声在频域里均匀铺开。写成公式就是 x(t) d(t) s(t) n(t)。故障诊断的全部目标就是从 x(t) 里把 s(t) 找出来测出它的重复周期反推损伤发生在外圈、内圈、滚动体还是保持架上。但问题在于s(t) 往往是最弱的那一部分能量占比可能连百分之几都不到。拿生活里的场景类比这相当于在一场嘈杂的宴会上你要听清角落里一个人的小声说话。背景音乐是确定性周期分量旁人的嘈杂是随机噪声那个说话声就是故障冲击。直接录下来回放听不清你需要先消掉背景音乐再把说话声音频放大。AR 和 MED 干的正是这两件事。1.2 AR负责剃掉确定周期MED负责放大冲击AR 模型在这里的角色是“确定性成分滤除器”。自回归模型用前 p 个采样点预测当前点本质上是拟合信号里线性可预测的部分。转频谐波、啮合频率这类确定性周期信号历史值与当前值有极强相关性AR 可以很好地把它们预测出来预测误差里就只剩下不可预测的冲击和随机噪声。这一步在文献里常叫预白化。MED 的角色是“冲击反滤波器”。它通过迭代优化寻找一个逆滤波器使输出信号的峭度最大。峭度反映信号分布的尖峰性稀疏的周期冲击叠加之后会让信号表现出很高的峭度而高斯噪声的峭度接近 3纯周期信号的峭度更低。MED 把这个特征作为优化目标迭代设计滤波器让冲击被放大、波形被锐化同时抑制高斯噪声。两条处理顺序不能颠倒。如果先做 MED确定性周期分量在信号里占主导算法会花大量迭代去处理周期成分冲击增强效果大打折扣如果只用 AR周期成分被剃掉但噪声依然很大冲击在时域里还是不够明显。AR 先清场、MED 再聚焦整个处理链才有意义。1.3 为什么两级联合比单用一种方法靠谱单独用高通滤波看似简单但有一个很难调的问题截止频率选低了转频谐波和啮合频率的高频分量还留在信号里选高了可能把故障冲击参与的高频共振成分一起切掉。滤波器设计参数一旦换一组转速或工况就得重新调非常脆弱。单独用 AR 的问题在于它只解决了“去掉不想要的周期分量”没有解决“想要的冲击被噪声淹没”的问题。AC 残差里噪声和冲击混杂包络谱依然难看。单独用 MED 的问题在于它虽然能提升峭度但原始信号里确定性成分如果太强MED 为了最大化峭度可能反而把某些周期分量也当作“尖峰”来处理。AR 和 MED 一个做减法、一个做增强正好互补。这也是很多论文把 ARMED 作为经典预处理链路的原因它不依赖过多先验参数对转速变化和故障类型变化都有一定鲁棒性。2. 算法原理拆解AR模型和MED的数学内核2.1 AR模型一个线性预测器如何抽出周期成分p 阶 AR 模型的数学形式是x(n) -a1·x(n-1) - a2·x(n-2) - … - ap·x(n-p) e(n)把 e(n) 单独提出来残差序列就是e(n) x(n) a1·x(n-1) a2·x(n-2) … ap·x(n-p)在 Matlab 里如果 AR 系数向量是 a [1, a1, a2, …, ap]那么 filter(a, 1, x) 得到的结果正好就是残差 e(n)。这是很多人第一次用 AR 做信号预处理时容易绕晕的地方filter 的分子是 a 而不是 [1]其实是因为 AR 模型的“预测误差滤波器”本身就是单位冲激响应为 a 的 FIR 滤波器。AR 系数估计最常用的是 Burg 算法Matlab 里对应 arburg 函数。Burg 算法基于前后向预测误差最小化对短数据的适应性比自相关法好也不容易遇到 Levinson-Durbin 递推中的病态情况。对于轴承振动这种平稳性尚可的信号arburg 得到的结果足够稳定。阶数 p 是唯一需要关心的参数。采样率 12k、转频 25Hz 时一个转频周期对应 480 个采样点如果 p 只有 10AR 根本没有足够的记忆长度去建模这么慢的周期分量残差里会留下明显的转频线谱。一般工程经验是 p 取 fs/f_min 的 2 到 5 倍其中 f_min 是你想滤除的最低确定性频率。也可以做一个简单的阶数扫描观察残差频谱里确定性线谱是否消失比单纯看公式可靠得多。2.2 MED最大化峭度的逆卷积迭代MED 最早是 Wiggins 在 1978 年提出来的用于地震勘探中恢复被大地滤波作用削弱的反射脉冲序列后来才被引入机械故障诊断。它要解决的问题是观测信号 x 其实是原始冲击源与传递路径、传感器响应卷积之后的“模糊版本”想反卷积得到冲击序列但反卷积本身是病态的必须加一个约束。MED 加的约束就是峭度最大化。目标函数写作J(f) E[y^4] / (E[y^2])^2这里 y f * x 是逆滤波器的输出。如果 y 是高斯分布J 趋近 3如果 y 里有一串稀疏大脉冲J 会显著高于 3。最大化 J 的迭代更新式可以写成f_new (X^T X)^{-1} · X^T·(y^3) / (y^T y)其中 X 是由 x 构成的 Toeplitz 卷积矩阵。直观理解这一步相当于以当前输出 y 的“三次方加权”作为期望目标再用最小二乘拟合出新的逆滤波器。每次更新后把 f 归一化防止幅度发散。重复迭代若干次滤波器的长度决定它能恢复多长范围内的冲击结构。这套迭代逻辑并不复杂代码实现量不到 20 行正是它在工业信号处理里流行的原因。需要提醒的是拉动峭度目标的每次迭代都相当于做一次矩阵求逆如果输入信号本身包含强线谱或者直流偏置自相关矩阵的条件数可能变得很差后面我会给出一个正则化处理方案。2.3 参数选择经验AR阶数和MED滤波器长度先看 AR 阶数。低采样率场景比如 8k 到 12kp 取 60 到 150 是常见区间高采样率场景比如 48k 到 96kp 可能要取到 300 到 600。p 太小滤不干净周期分量p 太大又可能出现过度拟合把故障冲击的部分相关性也建模进去导致残差里的冲击被削弱。实用技巧是画一张残差频谱图如果转频谐波对应的谱线已经低于噪声底说明阶数够用不必追求过大的 p。再看 MED 滤波器长度 L。这个参数决定逆滤波器的自由度一般取 10 到 60。L 太小滤波器无法充分反卷积冲击增强有限L 太大波形容易被过度锐化出现不存在的伪振荡计算量也会增加。我的默认值通常是 L 30效果稳定。迭代次数设 20 到 30 次完全够用算法一般十几步就会收敛。可以加一个终止条件相邻两次滤波器系数差的范数小于 1e-3 就提前跳出既省时间又避免数值漂移。3. Matlab实战从仿真信号到包络谱的完整链路3.1 第一步构建带故障冲击的仿真振动信号先用仿真信号把整个处理链路跑通确认每一个环节都正确再放到现场数据上才靠谱。以下是一个完整的示例脚本片段clc; clear; close all; %% 基本参数 fs 12000; % 采样频率 12kHz T 1; % 1秒数据 N fs * T; t (0:N-1) / fs; %% 轴承几何参数与故障特征频率 fr 25; % 转轴转速 25Hz Z 8; % 滚动体个数 d 10; % 滚动体直径 mm D 40; % 节圆直径 mm alpha 0; % 接触角 % 外圈故障特征频率 BPFO Z * fr / 2 * (1 - d/D * cos(alpha)); fprintf(BPFO %.2f Hz\n, BPFO);BPFO 算出来是 75Hz。这个数值用来生成故障冲击序列后面也要回到包络谱里找它。%% 确定性周期成分转频、谐波、额外的周期干扰 x_det 0.5 * sin(2*pi*fr*t) 0.3 * sin(2*pi*2*fr*t) 0.2 * sin(2*pi*3*fr*t); x_det x_det 0.3 * sin(2*pi*240*t); % 模拟一个高频确定性干扰 %% 故障冲击成分共振衰减振荡叠加 fn 3000; % 结构共振频率 zeta 0.08; % 阻尼比 wd fn * sqrt(1 - zeta^2); t_imp (0:120-1) / fs; % 10ms 冲击响应长度 h exp(-zeta*2*pi*fn*t_imp) .* sin(2*pi*wd*t_imp); h h / abs(min(h)); % 归一化到单位负峰值 Tf 1 / BPFO; % 故障冲击周期 imp_t 0:Tf:0.9; % 取到0.9s避免边界相位截断 imp_t imp_t 0.02 * Tf * randn(size(imp_t)); % 2%滑移抖动 x_imp zeros(1, N); for k 1:length(imp_t) idx round(imp_t(k) * fs) 1; if idx length(h) - 1 N x_imp(idx:idxlength(h)-1) x_imp(idx:idxlength(h)-1) h; end end %% 随机噪声 noise 0.2 * randn(1, N); %% 合成信号 x x_det x_imp noise;这里有两个细节值得说明。一个是冲击间隔加入 2% 的随机滑移抖动因为实际滚动体滑动会让冲击周期不是绝对均匀不加抖动生成的频谱会过于“干净”和现实脱节。另一个是 240Hz 的确定性干扰它模拟了与转频不同步的周期源比如齿轮啮合这能检验 AR 模型到底能不能把非转频相关的周期成分也剔掉。3.2 第二步AR残差提取把周期分量剥掉AR 处理的代码非常短%% AR 预白化 p 120; a_ar arburg(x, p); x_ar_res filter(a_ar, 1, x); x_ar_res(1:p) []; % 丢弃滤波起始段的瞬态效应filter(a_ar, 1, x) 的结果就等于 x(n) 减去 AR 模型预测值之后的误差序列。因为 a_ar 的第一个元素是 1以后面的系数构造出的是一个“预测误差滤波器”。如果 AR 模型已经充分捕捉了确定性周期分量这段残差里就应该主要剩下冲击加噪声原始信号里那些大能量的正弦成分会被压到接近噪声底水平。这也是验证 AR 模型是否选好阶数的最佳时机直接对残差做一次 FFT看转频和 240Hz 对应的线谱是否消失。如果还矗着一根很高的谱线说明 p 太小适当把 p 往上加。丢弃前 p 个点是因为滤波起始段需要填历史值最初的 p 个输出并不可靠处理完再补一个零均值化即可。3.3 第三步MED函数实现与调用MED 本质上是一个十几行的迭代但封装成函数之后调用更方便。function [y, f, kurt_iter] medFilter(x, L, maxIter) % x : 输入信号通常为 AR 残差 % L : 逆滤波器长度 % maxIter : 最大迭代次数 % y : 增强后的信号 % f : 最优逆滤波器系数 % kurt_iter: 每次迭代的输出峭度用于观察收敛 x x(:); N length(x); % 构造 Toeplitz 卷积矩阵 X toeplitz(x, [x(1); zeros(L-1, 1)]); R X * X; % 滤波器初始化相当于全通冲激 f zeros(L, 1); f(1) 1; f f / norm(f); kurt_iter zeros(maxIter, 1); for it 1:maxIter y X * f; % 逆滤波输出 b X * (y.^3) / (y * y); % 加权互相关向量 f_new R \ b; % 最小二乘更新 f_new f_new / norm(f_new); % 归一化防止幅度发散 kurt_iter(it) kurtosis(y); if norm(f_new - f) 1e-3 % 收敛判断 f f_new; break; end f f_new; end y X * f; end调用方式%% MED 增强 L 30; [x_med, f_med, kurt_list] medFilter(x_ar_res, L, 30);几个注意点。第一自相关矩阵 R 如果条件数过大可以在构造时加一个很小的对角正则项R X * X 1e-6 * eye(L);这会牺牲极少的精度但能避免数值上出现接近奇异的矩阵导致迭代发散。第二medFilter 的输出幅度和输入量纲不一致如果后续要保留原始幅值用于包络谱定量比较可以把输出按输入标准差做个缩放。第三输入信号最好先减去均值否则直流分量会影响峭度计算和矩阵条件数。我习惯在调用前直接 x x(:) - mean(x)。3.4 第四步包络谱提取锁定BPFO增强做完最后一步是包络谱分析。包络谱的思路是先把高频冲击波形通过希尔伯特变换取包络再做一次 FFT把故障冲击的重复频率从时域的“高频衰减抖动”转成频域的低频谱峰。%% 包络谱对比 env_raw abs(hilbert(x)); env_raw env_raw - mean(env_raw); spec_raw abs(fft(env_raw)); env_med abs(hilbert(x_med)); env_med env_med - mean(env_med); spec_med abs(fft(env_med)); f_axis (0:N/2-1) / N * fs; % 查看 BPFO 附近谱峰 [~, idx_bpfo_raw] max(spec_raw(round(f_axis/1)))这里用原始信号和 ARMED 处理后的信号分别做包络谱对比 BPFO 及其 2 倍频处的谱峰值。处理前 BPFO 处可能只是噪声底上稍稍冒头处理后应该是一根非常突出的谱线而且往往能清楚看到 BPFO、2×BPFO、3×BPFO 的谐波序列。在输出时可以通过 findpeaks 自动定位谱峰并打印对应的频率值人工核对是否落在计算出的故障特征频率附近。值得注意的是包络谱的低频段通常存在一个直流峰和转频附近的隆起画图时建议把 0 到 5Hz 以内的部分裁掉否则整张图的纵轴会被直流分量压扁真正的故障谱峰反而不明显。4. 处理效果分析峭度与包络谱在三个阶段的变化4.1 时域波形和峭度冲击从淹没到突出跑完上面的脚本你会在时域图里看到非常明显的变化。原始信号里冲击序列混在正弦波和噪声里肉眼能找到一些“毛刺”但幅度远不如确定性成分。AR 残差出来后正弦成分基本消失时域里剩下的是噪声加冲击冲击的峰已经比较显眼。MED 输出则干脆得多冲击位置会出现清晰的稀疏大脉冲噪声被明显压低波形看起来就像一串衰减振荡被一个个剥离出来。峭度值可以量化这个过程。我按上面的参数跑出来原始信号的峭度大约在 3.8AR 残差提升到 6.5 左右MED 输出能到 11 以上。峭度从接近高斯分布的水平跳到十几说明冲击已经从背景里脱颖而出。这里有个容易误解的点峭度并不是越高越好如果峭度过高往往意味着算法过度锐化把噪声里的个别尖峰也放大了。真实故障信号处理中峭度超过 15 就要警惕过处理。下面是一组典型结果对应表处理阶段峭度(约)时域形态原始信号3.8正弦周期分量占主导冲击被掩盖AR残差6.5周期分量消失冲击在噪声中可见ARMED输出11.2冲击成稀疏大脉冲噪声被抑制4.2 包络谱对比故障特征频率逐步清晰包络谱的变化比时域更直观。原始信号的包络谱里75Hz 附近可能有一条小峰但周围噪声底也很高转频及其谐波还会形成一些干扰谱线稍不注意就会把它和其他峰混淆。AR 残差的包络谱里转频成分明显减弱75Hz 的峰开始突出一些但噪声底仍然较高。到了 MED 输出75Hz、150Hz、225Hz 三根谱线形成清晰的谐波序列峰值可以比 AR 残差阶段高出 3 到 5 倍。我把典型幅度变化放在下面供你对照自己的实验结果谱峰原始信号AR残差ARMEDBPFO (75Hz)0.91.85.32×BPFO (150Hz)0.40.92.83×BPFO (225Hz)0.20.51.9可以看到处理链路让故障特征频率的信噪比提升非常明显。判断“故障是否成立”的常用规则就是看 BPFO 基频以及 2 到 3 个谐波是否同时出现。如果只有一根孤立的峰那可能是某个干扰信号不能直接下结论。4.3 信噪比拉低之后这套流程还能撑多久现场信号往往不会像仿真这么干净所以我习惯性做一组不同噪声强度下的稳定性实验。把噪声标准差从 0.1 调到 0.5甚至到 1.0观察 MED 输出峭度和 BPFO 谱峰值的变化。大致结果如下噪声标准差MED输出峭度BPFO谱峰值(约)识别难度0.113.46.8非常清晰0.211.25.3清晰0.57.22.9仍可识别1.05.01.4开始勉强这说明 ARMED 在中等噪声下表现稳定但一旦噪声接近甚至超过冲击幅值峭度提升和谱峰增强都会明显退化。这时不要只想着把 MED 迭代次数调大因为低信噪比下算法可能把噪声峰值也“增强”出来。更合理的做法是先用带通滤波把共振频带选出来再做 MED或者直接换成更鲁棒的变体算法比如 MCKD 或自适应 MED。不过对于大多数早期故障场景标准 ARMED 已经足够用。5. 实测中最容易踩的坑与排查技巧5.1 AR阶数不合适残差里还残留周期线谱最常见的坑就是 AR 阶数取太小。很多人照着论文里的 p 20 一跑发现残差信号的频谱里还杵着一堆转频谐波于是怀疑 AR 模型没用。其实问题就出在阶数上当采样率 12k、转频极低时一个周期有几百个采样点p 太小捕捉不到完整的周期关系。我的排查流程是先对残差做 FFT拉出频谱图盯着转频和已知的确定性干扰频率看。如果这些谱线高于噪声底 10dB 以上就逐步加大 p直到它们沉到噪声底下面。p 从 60 加到 120 通常就能解决问题。注意阶数不是越大越好p 太大时残差里的冲击也会被削弱。有个折中技巧p 增加后对比残差的峭度变化如果峭度不升反降说明过度拟合了需要回调。5.2 MED在低信噪比或强干扰下会跑偏MED 迭代本身很简洁但工程实现中容易遇到三类问题。第一类是自相关矩阵病态输入信号直流分量过大或者某条谱线太强时R 接近奇异迭代会出现剧烈振荡。解决办法是信号先去均值和加正则化项。第二类是信噪比过低时MED 输出的峭度提升不明显甚至把噪声里的个别尖峰当成冲击放大。这类情况下先带通滤波再进 MED 会有帮助。第三类是滤波器长度 L 取得太大输出波形出现“伪脉冲”看起来更稀疏实则失真。L 超过 60 之后每增加一倍自由度波形过冲风险都显著增加。一个实用的诊断方法是看每轮迭代的峭度曲线。正常情况下 kurt_iter 应该稳步上升然后趋于平稳。如果你发现峭度曲线先升后降大概率是数值问题或者输入里存在突变干扰如果曲线一直不上升说明输入信号里确实没有可被增强的冲击成分。5.3 包络谱判读不要被转频边带和伪峰带偏包络谱做出来之后真正考验人的是判读。故障特征频率附近经常会出现转频边带尤其是内圈故障或转子不平衡时BPFI 两侧会各偏一个转频形成 ±fr 的谱线。这时候不要看见一根峰就报故障要核对计算的故障特征频率以及谐波序列是否一致。另一个容易踩的坑是 AR 没有完全滤掉转频时包络谱低频段会出现较强的转频及其倍频谱线这些线在 MED 增强后也可能被放大。判读时先看 BPFO、BPFI、BSF、FTF 这几个理论值附近是否有对应的谐波序列再看异常谱线是否出现在这些频率的整数倍位置。通常按表格对照一下比较稳妥故障位置特征频率公式本组参数示例值外圈Z·fr/2·(1-d/D·cosα)75Hz内圈Z·fr/2·(1d/D·cosα)125Hz滚动体D·fr/(2d)·(1-(d/D·cosα)^2)46.9Hz保持架fr/2·(1-d/D·cosα)9.4Hz如果你的实测结果找不到任何一条特征频率处的谐波序列不要急着调 MED 的迭代次数先回头确认轴承参数和转速是否准确。很多时候参数算错了后面所有分析都是南辕北辙。这套 ARMED 流程我前前后后跑了很多遍最大的体会是它算不上什么“高级算法”但工程性价比很高代码量小、参数少、对大多数轴承早期故障都能把特征频率从噪声里捞出来。反而越是这种简洁的方法越考验每一步对信号的理解。AR 阶数、MED 滤波器长度、包络谱判读时的那点经验都是在一次次错判和返工中攒出来的。建议你先用仿真信号把整个链路跑通再换到现场数据上中间每一步的时域波形和频谱变化都要自己盯一眼。等你习惯了这套思路再去看那些更复杂的自适应滤波、深度学习方法就会发现它们本质上解决的是同一个问题——让故障冲击在噪声里现形。