在旋转机械振动诊断的现场，最磨人的往往不是设备真的停了，而是故障已经发生，你却从振动数据里看不出来。滚动轴承早期损伤产生的冲击非常微弱，被转频谐波、齿轮啮合成分和背景噪声层层盖住，直接拿原始信号做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 一个做减法、一个做增强，正好互补。这也是很多论文把 AR+MED 作为经典预处理链路的原因：它不依赖过多先验参数，对转速变化和故障类型变化都有一定鲁棒性。
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 只有 10，AR 根本没有足够的记忆长度去建模这么慢的周期分量，残差里会留下明显的转频线谱。一般工程经验是 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 到 12k，p 取 60 到 150 是常见区间；高采样率场景，比如 48k 到 96k，p 可能要取到 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:idx+length(h)-1) = x_imp(idx:idx+length(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)))这里用原始信号和 AR+MED 处理后的信号分别做包络谱，对比 BPFO 及其 2 倍频处的谱峰值。处理前 BPFO 处可能只是噪声底上稍稍冒头，处理后应该是一根非常突出的谱线，而且往往能清楚看到 BPFO、2×BPFO、3×BPFO 的谐波序列。在输出时，可以通过 findpeaks 自动定位谱峰并打印对应的频率值，人工核对是否落在计算出的故障特征频率附近。
值得注意的是，包络谱的低频段通常存在一个直流峰和转频附近的隆起，画图时建议把 0 到 5Hz 以内的部分裁掉，否则整张图的纵轴会被直流分量压扁，真正的故障谱峰反而不明显。
4. 处理效果分析：峭度与包络谱在三个阶段的变化
4.1 时域波形和峭度：冲击从淹没到突出
跑完上面的脚本，你会在时域图里看到非常明显的变化。原始信号里，冲击序列混在正弦波和噪声里，肉眼能找到一些“毛刺”，但幅度远不如确定性成分。AR 残差出来后，正弦成分基本消失，时域里剩下的是噪声加冲击，冲击的峰已经比较显眼。MED 输出则干脆得多，冲击位置会出现清晰的稀疏大脉冲，噪声被明显压低，波形看起来就像一串衰减振荡被一个个剥离出来。
峭度值可以量化这个过程。我按上面的参数跑出来，原始信号的峭度大约在 3.8，AR 残差提升到 6.5 左右，MED 输出能到 11 以上。峭度从接近高斯分布的水平跳到十几，说明冲击已经从背景里脱颖而出。这里有个容易误解的点：峭度并不是越高越好，如果峭度过高，往往意味着算法过度锐化，把噪声里的个别尖峰也放大了。真实故障信号处理中，峭度超过 15 就要警惕过处理。
下面是一组典型结果对应表：
| 处理阶段 | 峭度(约) | 时域形态 |
|---|---|---|
| 原始信号 | 3.8 | 正弦周期分量占主导，冲击被掩盖 |
| AR残差 | 6.5 | 周期分量消失，冲击在噪声中可见 |
| AR+MED输出 | 11.2 | 冲击成稀疏大脉冲，噪声被抑制 |
4.2 包络谱对比：故障特征频率逐步清晰
包络谱的变化比时域更直观。原始信号的包络谱里，75Hz 附近可能有一条小峰，但周围噪声底也很高，转频及其谐波还会形成一些干扰谱线，稍不注意就会把它和其他峰混淆。AR 残差的包络谱里，转频成分明显减弱，75Hz 的峰开始突出一些，但噪声底仍然较高。到了 MED 输出，75Hz、150Hz、225Hz 三根谱线形成清晰的谐波序列，峰值可以比 AR 残差阶段高出 3 到 5 倍。
我把典型幅度变化放在下面，供你对照自己的实验结果：
| 谱峰 | 原始信号 | AR残差 | AR+MED |
|---|---|---|---|
| BPFO (75Hz) | 0.9 | 1.8 | 5.3 |
| 2×BPFO (150Hz) | 0.4 | 0.9 | 2.8 |
| 3×BPFO (225Hz) | 0.2 | 0.5 | 1.9 |
可以看到，处理链路让故障特征频率的信噪比提升非常明显。判断“故障是否成立”的常用规则，就是看 BPFO 基频以及 2 到 3 个谐波是否同时出现。如果只有一根孤立的峰，那可能是某个干扰信号，不能直接下结论。
4.3 信噪比拉低之后，这套流程还能撑多久
现场信号往往不会像仿真这么干净，所以我习惯性做一组不同噪声强度下的稳定性实验。把噪声标准差从 0.1 调到 0.5，甚至到 1.0，观察 MED 输出峭度和 BPFO 谱峰值的变化。大致结果如下：
| 噪声标准差 | MED输出峭度 | BPFO谱峰值(约) | 识别难度 |
|---|---|---|---|
| 0.1 | 13.4 | 6.8 | 非常清晰 |
| 0.2 | 11.2 | 5.3 | 清晰 |
| 0.5 | 7.2 | 2.9 | 仍可识别 |
| 1.0 | 5.0 | 1.4 | 开始勉强 |
这说明 AR+MED 在中等噪声下表现稳定，但一旦噪声接近甚至超过冲击幅值，峭度提升和谱峰增强都会明显退化。这时不要只想着把 MED 迭代次数调大，因为低信噪比下算法可能把噪声峰值也“增强”出来。更合理的做法是先用带通滤波把共振频带选出来，再做 MED，或者直接换成更鲁棒的变体算法，比如 MCKD 或自适应 MED。不过对于大多数早期故障场景，标准 AR+MED 已经足够用。
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·(1+d/D·cosα) | 125Hz |
| 滚动体 | D·fr/(2d)·(1-(d/D·cosα)^2) | 46.9Hz |
| 保持架 | fr/2·(1-d/D·cosα) | 9.4Hz |
如果你的实测结果找不到任何一条特征频率处的谐波序列，不要急着调 MED 的迭代次数，先回头确认轴承参数和转速是否准确。很多时候参数算错了，后面所有分析都是南辕北辙。
这套 AR+MED 流程我前前后后跑了很多遍，最大的体会是它算不上什么“高级算法”，但工程性价比很高：代码量小、参数少、对大多数轴承早期故障都能把特征频率从噪声里捞出来。反而越是这种简洁的方法，越考验每一步对信号的理解。AR 阶数、MED 滤波器长度、包络谱判读时的那点经验，都是在一次次错判和返工中攒出来的。建议你先用仿真信号把整个链路跑通，再换到现场数据上，中间每一步的时域波形和频谱变化都要自己盯一眼。等你习惯了这套思路，再去看那些更复杂的自适应滤波、深度学习方法，就会发现它们本质上解决的是同一个问题——让故障冲击在噪声里现形。
