做信号分解的人大概率都绕不开这么个场景手里拿到一段轴承振动数据想看故障频率但原始信号里既有转频、又有轴承固有谐振的衰减振荡还叠着随机噪声。直接做FFT会发现频谱乱成一锅粥能量都铺在宽频带上。传统做法是先做带通滤波再包络谱但滤波频带怎么选本身就是一门玄学。LMDLocal Mean Decomposition局部均值分解提供了一个很直观的思路——把信号分解成若干个具有物理意义的调幅调频分量PF分量逐个分析瞬时频率和瞬时幅值。这篇博文我会把LMD的数学原理用通俗的方式拆开讲清楚给出完整可运行的MATLAB实现再用滚动轴承故障仿真信号做一次全流程演示最后把端点效应、模态混叠、窗口选择这些实际工程里的坑都过一遍。适合正在做机械故障诊断、非平稳信号分析的研究生和工程师参考。1. LMD和EMD有什么本质区别先搞懂它解决什么问题1.1 非平稳信号为什么不能直接分析想象你录了一段人声想从中提取说话者音调随时间的变化。直接把整段语音做FFT只能看到300Hz到3400Hz的一个宽频包络根本看不出每个字音调是高是低。因为FFT假设信号是平稳的而现实中的语音、振动、心电信号到处都是突变和频率变化。更麻烦的是Hilbert变换。它虽然能算瞬时频率但有个前提信号在任意时刻只能有一个主导频率也就是“单分量信号”。真实信号哪有那么听话轴承振动里同时有转频、啮合频率、故障冲击引发的共振衰减波这几个成分叠在一起直接做Hilbert变换得到的瞬时频率会在多个频率之间来回跳完全没有物理意义。所以思路很自然先做“分解”把复杂信号拆成多个单分量再逐一对每个分量做瞬时幅值和瞬时频率分析。这也是Hilbert-Huang变换那套思路的核心。LMD就是这条路上一个非常有特点的算法。1.2 LMD的产品化表达PF分量是自带瞬时幅值的AM-FM信号LMD由Smith在2005年提出它的输出不是EMD那种IMF而是一组乘积函数Product Function简称PF。每个PF分量理论上都可以写成PF(t) a(t) · s(t)其中a(t)是慢变的瞬时幅值包络s(t)是瞬时频率随时间变化的纯调频信号幅值恒为1。这种“包络乘载波”的表达方式对工程信号特别友好。举个例子齿轮箱振动里常见的幅值调制一个齿轮齿面出现局部磨损每旋转一圈会产生一次冲击这个冲击会激发高频结构共振同时冲击强度又随载荷缓慢波动。在波形上看到的就是“高频振荡的幅度被低频信号包住”。LMD能把这种信号直接拆成“包络×调频项”包络对应故障冲击强度调频项对应振动固有频率物理意义非常清晰。1.3 和EMD的关键对比EMD在实际中用得最多所以很多人的第一反应是拿LMD和EMD比较。两者的根本差异在于筛分构建方式EMD是找上下包络然后取均值LMD是计算相邻极值点的局部均值和包络估计再做迭代解调。对比维度EMDLMD分量类型IMF本征模态函数PF乘积函数分量结构调幅调频但不严格要求显式表达为包络×纯调频瞬时频率获取对IMF做Hilbert变换对纯调频项的相位求导核心构建方式上下包络均值局部均值包络估计迭代常见痛点模态混叠、端点效应局部均值构造方式敏感结果可解释性依赖后处理物理意义更直观我自己做轴承故障诊断时最大的感受是EMD分解出来的IMF在端点处经常出现上下包络交叉、瞬时频率为负值这类怪现象。LMD因为做了“除以包络”这一步把信号强制压到单位幅值附近再求瞬时频率时稳定性会好不少。当然LMD也不是没有代价它对“局部均值函数怎么构造”这件事非常敏感这也是后面要重点讲的部分。2. LMD算法原理拆解局部均值、包络估计与迭代解调2.1 从极值点出发局部均值和包络估计为什么这么算LMD的第一步是找信号s(t)的所有局部极值点n_i包括极大值和极小值。然后对相邻两个极值点n_i和n_{i1}做两个简单运算局部均值m_i (n_i n_{i1}) / 2包络估计a_i |n_i - n_{i1}| / 2这个式子的直觉其实特别好理解。相邻一个极大值一个极小值它们的中点大体就是这段局部波形围绕的“中心位置”所以叫局部均值而极大值和极小值幅度差的一半就是这半个周期里信号振幅的大致估计也就是包络。举个数字例子如果一段信号的极值点依次是1.2、-0.8、1.0、-0.7那么第一个局部均值是0.2第一个包络估计是1.0第二个局部均值是0.1第二个包络估计是0.85。把这些离散的m_i和a_i值连成随时间的连续曲线就得到了局部均值函数m_11(t)和包络估计函数a_11(t)。连成曲线这一步就是LMD实现里最核心也是最容易出问题的地方。Smith原始论文用的是滑动平均但后来的实践表明用三次样条插值通常更稳定。后面第4章会详细讲两者的差异。2.2 迭代解调把幅值变化“除”掉逼近纯调频信号得到局部均值函数和包络估计函数之后LMD进入内层迭代从原始信号中减去局部均值函数h_11(t) s(t) - m_11(t)用包络估计函数归一化s_11(t) h_11(t) / a_11(t)第二步是LMD的精髓。除以包络的目的是把信号的幅值变化“压平”让信号变成幅度恒为1的纯调频信号。如果做完一次之后信号还是带有明显的幅值波动说明包络没剥离干净那就用s_11(t)作为新的输入重复找极值点、算局部均值、算包络估计、再减均值、再除以包络一直迭代下去。这就像把一段录音先做自动增益控制不管音量是忽大忽小先把响度拉平只留下音调和节奏信息。工程信号往往包含多级调制所以一次“拉平”不够必须反复迭代。内层迭代什么时候停看包络估计函数a_1(n1)(t)是否在整个时间范围内都接近1。通常的判断条件是max|a_1(n1)(t) - 1| ΔΔ一般取0.001到0.01。如果取得太大比如0.05包络没剥离干净PF分量的幅值包络会残留毛刺如果取得太小迭代次数会暴增甚至因为数值精度问题永远不收敛。2.3 累积包络与PF输出残差剥离内层迭代结束时把每次迭代得到的所有包络估计函数乘起来就得到这个PF分量的瞬时幅值a_1(t) a_11(t) · a_12(t) · ... · a_1n(t)而最后一次迭代得到的s_1n(t)就是纯调频信号。两者相乘得到第一个PF分量PF_1(t) a_1(t) · s_1n(t)然后从原始信号中减去PF_1得到残差u_1(t) x(t) - PF_1(t)对残差重复整个分解过程得到PF_2、PF_3……直到残差信号没有足够多的极值点或者能量足够小。最终原始信号可以重构为x(t) Σ PF_i(t) u_k(t)整个分解过程是自适应的不需要提前指定要分解出多少个分量算法会按照信号本身的复杂度逐层剥离。2.4 内外两层终止条件的设计意图LMD有两层循环每层循环都需要注意终止条件。内层循环的终止条件是包络估计函数趋近1表示信号已经变成纯调频信号。这个条件的物理含义是“幅值调制已经被完全剥离”。外层循环的终止条件是残差信号极值点数量不足或者信号能量低于预设阈值表示“剩下的成分已经无法再分解出有意义的调幅调频分量”。我在实现代码时给内层循环设了最大迭代次数上限通常取50到200。这样做是为了防止算法在某些极端信号下发散或进入死循环比如信号幅值接近0、包络估计函数出现极小值等情况。代码里一旦触发上限就强制退出取当前结果作为近似PF后续可以通过可视化判断这个分量是否可信。3. MATLAB代码实现从零手写一个LMD分解器3.1 工具箱依赖与函数选型在MATLAB中实现LMD核心依赖并不复杂findpeaks函数来自Signal Processing Toolbox用来提取局部极值点spline函数是MATLAB基础函数用来做三次样条插值。如果没有Signal Processing Toolbox可以自己写一个基于diff符号变化的极值点检测函数逻辑不复杂但要注意处理平台段和端点情况。主函数我命名为lmd_decompose输入待分解信号、最大PF数量、纯调频判断阈值输出PF分量矩阵和残差信号。3.2 完整主函数代码function [PFs, residue] lmd_decompose(x, max_pf, threshold) % LMD: Local Mean Decomposition (局部均值分解) % 输入: % x - 待分解的一维信号(双精度向量) % max_pf - 最大PF分量数量, 默认8 % threshold - 纯调频判断阈值, 默认0.001 % 输出: % PFs - PF分量矩阵, 每行一个分量 % residue - 分解后的残差信号 if nargin 2 || isempty(max_pf), max_pf 8; end if nargin 3 || isempty(threshold), threshold 0.001; end x x(:); N length(x); PFs zeros(0, N); residue x; for p 1:max_pf s residue; a_product ones(size(s)); % 累积包络: 内层各次包络估计函数的乘积 s_new s; % 初始化, 防止第一轮就终止时未定义 for iter 1:200 % ---------- 1. 提取局部极值点 ---------- % findpeaks找极大值; 对-s找findpeaks即为极小值 [max_locs, max_vals] findpeaks(s); [min_locs, min_vals] findpeaks(-s); min_vals -min_vals; % 合并并按位置排序 all_locs [max_locs, min_locs]; all_vals [max_vals, min_vals]; [all_locs, idx] sort(all_locs); all_vals all_vals(idx); if length(all_locs) 4 break; % 极值点数量不足, 无法继续 end % ---------- 2. 相邻极值点的局部均值与包络估计 ---------- loc_m (all_vals(1:end-1) all_vals(2:end)) / 2; loc_a abs(all_vals(1:end-1) - all_vals(2:end)) / 2; if length(loc_m) 3 break; % 点数太少, 样条插值不稳定 end % 局部均值/包络估计的位置取相邻极值点中点 t_m (all_locs(1:end-1) all_locs(2:end)) / 2; t_a t_m; % ---------- 3. 三次样条插值构造连续函数 ---------- % 端点镜像拓延, 缓解边界效应 t_m_ext [2*t_m(1)-t_m(2), t_m, 2*t_m(end)-t_m(end-1)]; loc_m_ext [loc_m(1), loc_m, loc_m(end)]; t_a_ext [2*t_a(1)-t_a(2), t_a, 2*t_a(end)-t_a(end-1)]; loc_a_ext [loc_a(1), loc_a, loc_a(end)]; m_interp spline(t_m_ext, loc_m_ext, 1:N); a_interp spline(t_a_ext, loc_a_ext, 1:N); a_interp max(a_interp, eps); % 防止除零/负幅值 % ---------- 4. 减均值并除以包络 ---------- h s - m_interp; s_new h ./ a_interp; % ---------- 5. 累积包络 ---------- a_product a_product .* a_interp; % ---------- 6. 纯调频判断 ---------- if max(abs(a_interp - 1)) threshold break; end s s_new; end % ---------- 生成当前PF分量 ---------- pf a_product .* s_new; PFs(p, :) pf; % ---------- 更新残差 ---------- residue residue - pf; % ---------- 残差极值点过少时退出 ---------- [max_locs, ~] findpeaks(residue); [min_locs, ~] findpeaks(-residue); if length(max_locs) length(min_locs) 4 break; end end end代码里有几个地方值得说明。一是findpeaks(-s)这个技巧MATLAB的findpeaks只能找局部极大值想找极小值就取负号再找极大值得到结果再取负还原。二是端点镜像拓延这个操作会在第4章详细讲目的是不让spline在边界处出现大幅摆动。3.3 仿真验证调幅-调频叠加信号写一段仿真信号来验证代码是否正常工作。构造三个叠加成分一个80Hz载波、受8Hz幅值调制和20Hz相位调制的信号一个150Hz正弦一个衰减的280Hz振荡再加少量噪声。fs 1000; t (0:999)/fs; % 三个叠加成分 x1 (1 0.5*cos(2*pi*8*t)) .* cos(2*pi*80*t 2*sin(2*pi*20*t)); x2 0.2 * sin(2*pi*150*t); x3 0.8 * exp(-8*t) .* cos(2*pi*280*t); x x1 x2 x3 0.01*randn(size(t)); [PFs, residue] lmd_decompose(x, 5, 0.001); figure; for k 1:size(PFs,1) subplot(size(PFs,1)1, 1, k); plot(t, PFs(k,:)); ylabel(sprintf(PF%d, k)); end subplot(size(PFs,1)1, 1, size(PFs,1)1); plot(t, residue); ylabel(residue); xlabel(Time (s));运行之后可以看到前几个PF分量分别对应三个主要成分最后一个PF加残差对应噪声。这里有个注意点LMD的分解顺序不是按输入信号的频率高低严格排列的哪个分量先被剥离取决于每层残差中哪个振荡成分占主导。分析时不要想当然地认为PF1就一定是最高频分量要结合波形和频谱来看。如果想验证PF1的瞬时幅值是否解调正确可以这样analytic hilbert(PFs(1,:)); env abs(analytic); inst_phase unwrap(angle(analytic)); inst_freq diff(inst_phase) / (2*pi) * fs; figure; subplot(2,1,1); plot(t, env); ylabel(瞬时幅值); subplot(2,1,2); plot(t(2:end), inst_freq); ylabel(瞬时频率(Hz));对于第一个PF分量瞬时幅值应该近似10.5cos(2π·8t)瞬时频率应该近似8040cos(2π·20t) Hz。如果这两个曲线都符合预期说明LMD的核心逻辑没问题。4. 工程中的坑端点效应、局部均值构造与模态混叠4.1 端点效应为什么两端总是先“飞”LMD和EMD一样最让人头疼的就是端点效应。信号两端没有完整的极值点信息三次样条在插值边界时就会出现较大误差这个误差会向内传播导致分解结果在时间轴两端明显失真。我第一次跑LMD的时候就发现了这个现象PF分量中间段还挺正常头部和尾部却出现大幅摆动幅值甚至比原始信号还大几倍。原因就是极值点在端点处缺失样条插值在边界处外延得不到控制。缓解端点效应的常用手段是镜像拓延。思路是把信号两端的波形做镜像对称向外延长一段让边界处有足够的极值点参与插值。具体的实现可以有两种方式在极值点层面拓延对局部均值/包络估计点做镜像延拓我代码里采用的方式在信号层面拓延先对原始信号做镜像延拓再找极值点分解完再去掉延拓部分方式二效果通常更好但计算量稍大。实践中的做法是每端多延拓2到3个极值点间距的长度。延拓太短端点效应压不住延拓太长计算浪费且可能引入不相关的失真。4.2 局部均值函数构造滑动平均还是三次样条这是LMD复现中最大的一个坑。Smith原始论文里用滑动平均来构造局部均值函数和包络估计函数具体做法是对离散的m_i和a_i序列做多窗口移动平均然后再插值到全时间轴。问题在于滑动平均的窗口长度怎么选论文里没有给严格的自适应原则不同窗口对分解结果影响巨大。窗口太小均值函数和包络估计函数带有大量毛刺迭代解调容易发散窗口太大信号被过度平滑细节丢失分解出来的PF分量模糊。更麻烦的是极值点分布通常不均匀固定窗口长度很难同时适应密集段和稀疏段。实际使用中三次样条插值已经成为替代滑动平均的主流做法。它通过所有离散均值点生成光滑的连续曲线不需要人为设置窗口长度实现也更简洁。我在第3章代码里采用的就是三次样条。如果你看到某篇论文的LMD复现结果奇奇怪怪先检查它是不是用了固定窗口滑动平均——很多复现失败的根源都在这里。4.3 模态混叠与过分解判断与处理模态混叠指本应属于同一物理成分的信号被拆到多个PF分量中或者不同频率成分混进同一个PF分量。LMD的迭代解调机制比EMD稍微抗混叠但遇到频率成分较近、或者噪声能量较强时同样会翻车。过分解则是另一个极端算法把噪声也分解成看似规律的PF分量。这在中高频段尤其常见。判断方法很简单把分解结果和原始信号放在一起看如果某个PF分量幅值很小、波形杂乱、没有清晰的频谱主峰基本可以判断是过分解出来的噪声分量。处理手段有几个方向分解前做轻度的带通预滤波去掉明显无关的频带把纯调频阈值Δ调大一点让内层迭代早点收敛减少无效分解限制外层最大PF数量防止算法无限拆下去如果模态混叠严重可以考虑加白噪声的集合平均策略类EEMD思路但计算代价会高很多4.4 参数调节建议我整理了一份日常调参时使用的参考表不同信号类型可以根据实际效果上下浮动参数参考取值作用与调节思路纯调频阈值Δ0.001~0.01控制内层迭代精度噪声大时调大内层最大迭代次数50~200防止死循环信号平稳时可用较小值外层最大PF数5~10防止过分解按物理机理估计分量数端点拓延宽度2~3个极值点间距缓解端点效应信号越长可以越宽松包络下限eps或1e-6防止除以零数值保护如果你发现分解结果对参数非常敏感建议先固定阈值和迭代次数只调节最大PF数减少调参维度。等流程跑通后再回来微调其他参数。5. 实战演练滚动轴承故障信号的特征提取5.1 故障特征频率与仿真信号滚动轴承外圈故障的特征频率BPFO近似为BPFO (n_r / 2) · f_r · (1 - d/D · cosα)其中n_r是滚动体数量f_r是转频d是滚动体直径D是节径α是接触角。为了演示方便我直接设定转频30Hz外圈故障特征频率118Hz结构共振频率2000Hz。仿真信号构造思路每个故障周期产生一次冲击冲击在结构共振频率处引发衰减振荡同时冲击幅度受转频调制。最终再加上一点白噪声fs 10000; t (0:0.5*fs-1)/fs; fr 30; f_bpfo 118; fn 2000; zeta 0.05; impulse zeros(size(t)); T 1/f_bpfo; for k 0:floor(t(end)/T) tk k * T; idx find(t tk, 1); % 冲击起始点 time_local t(idx:end) - tk; % 冲击后局部时间 damped exp(-2*pi*fn*zeta*time_local) .* cos(2*pi*fn*time_local); len min(length(damped), length(t)-idx1); impulse(idx:idxlen-1) impulse(idx:idxlen-1) damped(1:len); end x (1 0.3*cos(2*pi*fr*t)) .* impulse 0.02*randn(size(t));这个信号的时域波形上可以看到周期性的冲击衰减频谱上则是以2000Hz为中心的一大片高频能量。直接用FFT很难看出118Hz故障特征因为故障特征体现在冲击的重复频率上而不是载波频率上。5.2 LMD分解与包络谱诊断用LMD把仿真信号分解成PF分量然后对第一个占主导的PF分量做包络谱分析[PFs, ~] lmd_decompose(x, 5, 0.005); env abs(hilbert(PFs(1,:))); Nfft 2^nextpow2(length(env)); S abs(fft(env, Nfft)); f_axis (0:Nfft/2-1) * fs / Nfft; figure; plot(f_axis(1:400), S(1:400)); xlim([0 400]); xlabel(频率(Hz)); ylabel(幅值);在包络谱中118Hz及其倍频236Hz、354Hz处会出现明显的谱峰这就是外圈故障的典型特征。如果进一步对多个PF分量分别做包络谱观察哪个分量在故障特征频率处能量最突出就能判断故障冲击主要调制在哪个频带上。这个方法比直接对原始信号做包络谱更稳健。因为LMD已经把最相关的调幅调频分量从其他干扰中分离出来包络谱的谱峰更干净信噪比更高。5.3 实测中要注意的细节用LMD做轴承故障诊断有几个细节非常重要采样率要足够高。LMD要提取冲击激发的高频共振如果采样率只有1kHz而共振频率在2kHz以上信号已经被严重混叠分解结果毫无意义。建议采样率至少是最高关注频率的5倍以上。信号长度要覆盖足够多的故障周期。如果只截取两三个冲击周期LMD分解和包络谱都会因为样本太少而失真。一般至少保证50个故障周期的长度。噪声不能完全无视。LMD对强噪声比较敏感实测信号通常比仿真信号脏得多。我的经验是先做轻度带通滤波把明显无关的频带去掉再做LMD效果会好很多。但如果滤波带宽本身选错了又回到了传统方法的玄学问题。所以滤波带宽宜宽不宜窄只去掉极端高频噪底和极低频趋势即可。6. LMD的适用范围、改进方向与我的经验6.1 什么时候用LMD什么时候换VMD/EMDLMD不是万能的选型时要看信号特点。信号场景推荐算法原因强调幅调频、故障冲击明显LMDPF分量物理意义清晰包络解调方便多通道同步信号多元LMD或VMD保证各通道分量一致性无调制的平稳叠加信号VMDLMD可能过度分解宽频强噪声VMD或EEMDLMD对噪声敏感需要严格数学最优VMD变分框架更规范我自己的习惯是EMD、LMD、VMD各跑一遍对比分解结果后选择物理意义最清晰的那个。听起来麻烦但实际多写几行脚本的事却能让结果分析可靠很多。6.2 值得尝试的改进方向如果要把LMD用到科研项目里下面几个方向值得探索第一用三次样条插值替代滑动平均这是目前最实用的改进也是我在代码里采用的方式。第二针对噪声较强的情况引入集合平均思想类似EEMD那样加入有限次白噪声再取平均可以有效缓解模态混叠。第三开发自适应阈值选择机制根据信号噪声能量估计自动确定内层迭代停止阈值。第四推广到多元LMD同时处理多通道振动信号避免各通道独立分解导致的分量错位问题。6.3 几个容易被忽视的实操细节最后分享几个我自己踩过的坑。第一次实现LMD时我给内层循环设了最大迭代次数但设得很大结果遇到一段幅值接近零的低能量信号迭代一两百次都停不下来整个脚本卡死。后来把上限压到200并结合一个残差能量判断问题才解决。现在我的代码里既能按阈值收敛也能强制退出不会因为个别信号把整个分析流程卡掉。阈值方面我也做过对比实验Δ取0.01时PF分量的包络曲线会显得比较毛糙瞬时幅值上能看到细微的锯齿调到0.001后包络光滑很多但迭代时间几乎翻倍。对于一般工程诊断0.005是一个不错的折中选择如果追求精细分析再用0.001。极值点提取这一步看似简单实际也容易出问题。比如信号带有直流偏置或趋势项时极值点的均匀性会变差分解出的第一个PF可能把趋势当成分量。建议在分解前先去趋势或做一次高通滤波把直流和超低频趋势去掉LMD的分解质量会明显提升。说实话LMD在MATLAB里手写并不难真正难的是参数选择和端点处理很多论文复现不了往往就是这几个细节没处理好。如果你也遇到分解结果乱跑的怪现象先检查端点拓延和阈值再检查极值点提取是否太粗糙。调通之后拿来做故障特征提取还是相当顺手的——至少在我处理轴承振动信号的经验里它比EMD稳定不少结果也更接近物理直觉。