做信號分解的人大概率都繞不開這么個場景手里拿到一段軸承振動數據想看故障頻率但原始信號里既有轉頻、又有軸承固有諧振的衰減振蕩還疊著隨機噪聲。直接做FFT會發(fā)現頻譜亂成一鍋粥能量都鋪在寬頻帶上。傳統(tǒng)做法是先做帶通濾波再包絡譜但濾波頻帶怎么選本身就是一門玄學。LMDLocal Mean Decomposition局部均值分解提供了一個很直觀的思路——把信號分解成若干個具有物理意義的調幅調頻分量PF分量逐個分析瞬時頻率和瞬時幅值。這篇博文我會把LMD的數學原理用通俗的方式拆開講清楚給出完整可運行的MATLAB實現再用滾動軸承故障仿真信號做一次全流程演示最后把端點效應、模態(tài)混疊、窗口選擇這些實際工程里的坑都過一遍。適合正在做機械故障診斷、非平穩(wěn)信號分析的研究生和工程師參考。1. LMD和EMD有什么本質區(qū)別先搞懂它解決什么問題1.1 非平穩(wěn)信號為什么不能直接分析想象你錄了一段人聲想從中提取說話者音調隨時間的變化。直接把整段語音做FFT只能看到300Hz到3400Hz的一個寬頻包絡根本看不出每個字音調是高是低。因為FFT假設信號是平穩(wěn)的而現實中的語音、振動、心電信號到處都是突變和頻率變化。更麻煩的是Hilbert變換。它雖然能算瞬時頻率但有個前提信號在任意時刻只能有一個主導頻率也就是“單分量信號”。真實信號哪有那么聽話軸承振動里同時有轉頻、嚙合頻率、故障沖擊引發(fā)的共振衰減波這幾個成分疊在一起直接做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。這種“包絡乘載波”的表達方式對工程信號特別友好。舉個例子齒輪箱振動里常見的幅值調制一個齒輪齒面出現局部磨損每旋轉一圈會產生一次沖擊這個沖擊會激發(fā)高頻結構共振同時沖擊強度又隨載荷緩慢波動。在波形上看到的就是“高頻振蕩的幅度被低頻信號包住”。LMD能把這種信號直接拆成“包絡×調頻項”包絡對應故障沖擊強度調頻項對應振動固有頻率物理意義非常清晰。1.3 和EMD的關鍵對比EMD在實際中用得最多所以很多人的第一反應是拿LMD和EMD比較。兩者的根本差異在于篩分構建方式EMD是找上下包絡然后取均值LMD是計算相鄰極值點的局部均值和包絡估計再做迭代解調。對比維度EMDLMD分量類型IMF本征模態(tài)函數PF乘積函數分量結構調幅調頻但不嚴格要求顯式表達為包絡×純調頻瞬時頻率獲取對IMF做Hilbert變換對純調頻項的相位求導核心構建方式上下包絡均值局部均值包絡估計迭代常見痛點模態(tài)混疊、端點效應局部均值構造方式敏感結果可解釋性依賴后處理物理意義更直觀我自己做軸承故障診斷時最大的感受是EMD分解出來的IMF在端點處經常出現上下包絡交叉、瞬時頻率為負值這類怪現象。LMD因為做了“除以包絡”這一步把信號強制壓到單位幅值附近再求瞬時頻率時穩(wěn)定性會好不少。當然LMD也不是沒有代價它對“局部均值函數怎么構造”這件事非常敏感這也是后面要重點講的部分。2. LMD算法原理拆解局部均值、包絡估計與迭代解調2.1 從極值點出發(fā)局部均值和包絡估計為什么這么算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值連成隨時間的連續(xù)曲線就得到了局部均值函數m_11(t)和包絡估計函數a_11(t)。連成曲線這一步就是LMD實現里最核心也是最容易出問題的地方。Smith原始論文用的是滑動平均但后來的實踐表明用三次樣條插值通常更穩(wěn)定。后面第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)作為新的輸入重復找極值點、算局部均值、算包絡估計、再減均值、再除以包絡一直迭代下去。這就像把一段錄音先做自動增益控制不管音量是忽大忽小先把響度拉平只留下音調和節(jié)奏信息。工程信號往往包含多級調制所以一次“拉平”不夠必須反復迭代。內層迭代什么時候?？窗j估計函數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有兩層循環(huán)每層循環(huán)都需要注意終止條件。內層循環(huán)的終止條件是包絡估計函數趨近1表示信號已經變成純調頻信號。這個條件的物理含義是“幅值調制已經被完全剝離”。外層循環(huán)的終止條件是殘差信號極值點數量不足或者信號能量低于預設閾值表示“剩下的成分已經無法再分解出有意義的調幅調頻分量”。我在實現代碼時給內層循環(huán)設了最大迭代次數上限通常取50到200。這樣做是為了防止算法在某些極端信號下發(fā)散或進入死循環(huán)比如信號幅值接近0、包絡估計函數出現極小值等情況。代碼里一旦觸發(fā)上限就強制退出取當前結果作為近似PF后續(xù)可以通過可視化判斷這個分量是否可信。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; % 極值點數量不足, 無法繼續(xù) 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; % 點數太少, 樣條插值不穩(wěn)定 end % 局部均值/包絡估計的位置取相鄰極值點中點 t_m (all_locs(1:end-1) all_locs(2:end)) / 2; t_a t_m; % ---------- 3. 三次樣條插值構造連續(xù)函數 ---------- % 端點鏡像拓延, 緩解邊界效應 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. 工程中的坑端點效應、局部均值構造與模態(tài)混疊4.1 端點效應為什么兩端總是先“飛”LMD和EMD一樣最讓人頭疼的就是端點效應。信號兩端沒有完整的極值點信息三次樣條在插值邊界時就會出現較大誤差這個誤差會向內傳播導致分解結果在時間軸兩端明顯失真。我第一次跑LMD的時候就發(fā)現了這個現象PF分量中間段還挺正常頭部和尾部卻出現大幅擺動幅值甚至比原始信號還大幾倍。原因就是極值點在端點處缺失樣條插值在邊界處外延得不到控制。緩解端點效應的常用手段是鏡像拓延。思路是把信號兩端的波形做鏡像對稱向外延長一段讓邊界處有足夠的極值點參與插值。具體的實現可以有兩種方式在極值點層面拓延對局部均值/包絡估計點做鏡像延拓我代碼里采用的方式在信號層面拓延先對原始信號做鏡像延拓再找極值點分解完再去掉延拓部分方式二效果通常更好但計算量稍大。實踐中的做法是每端多延拓2到3個極值點間距的長度。延拓太短端點效應壓不住延拓太長計算浪費且可能引入不相關的失真。4.2 局部均值函數構造滑動平均還是三次樣條這是LMD復現中最大的一個坑。Smith原始論文里用滑動平均來構造局部均值函數和包絡估計函數具體做法是對離散的m_i和a_i序列做多窗口移動平均然后再插值到全時間軸。問題在于滑動平均的窗口長度怎么選論文里沒有給嚴格的自適應原則不同窗口對分解結果影響巨大。窗口太小均值函數和包絡估計函數帶有大量毛刺迭代解調容易發(fā)散窗口太大信號被過度平滑細節(jié)丟失分解出來的PF分量模糊。更麻煩的是極值點分布通常不均勻固定窗口長度很難同時適應密集段和稀疏段。實際使用中三次樣條插值已經成為替代滑動平均的主流做法。它通過所有離散均值點生成光滑的連續(xù)曲線不需要人為設置窗口長度實現也更簡潔。我在第3章代碼里采用的就是三次樣條。如果你看到某篇論文的LMD復現結果奇奇怪怪先檢查它是不是用了固定窗口滑動平均——很多復現失敗的根源都在這里。4.3 模態(tài)混疊與過分解判斷與處理模態(tài)混疊指本應屬于同一物理成分的信號被拆到多個PF分量中或者不同頻率成分混進同一個PF分量。LMD的迭代解調機制比EMD稍微抗混疊但遇到頻率成分較近、或者噪聲能量較強時同樣會翻車。過分解則是另一個極端算法把噪聲也分解成看似規(guī)律的PF分量。這在中高頻段尤其常見。判斷方法很簡單把分解結果和原始信號放在一起看如果某個PF分量幅值很小、波形雜亂、沒有清晰的頻譜主峰基本可以判斷是過分解出來的噪聲分量。處理手段有幾個方向分解前做輕度的帶通預濾波去掉明顯無關的頻帶把純調頻閾值Δ調大一點讓內層迭代早點收斂減少無效分解限制外層最大PF數量防止算法無限拆下去如果模態(tài)混疊嚴重可以考慮加白噪聲的集合平均策略類EEMD思路但計算代價會高很多4.4 參數調節(jié)建議我整理了一份日常調參時使用的參考表不同信號類型可以根據實際效果上下浮動參數參考取值作用與調節(jié)思路純調頻閾值Δ0.001~0.01控制內層迭代精度噪聲大時調大內層最大迭代次數50~200防止死循環(huán)信號平穩(wěn)時可用較小值外層最大PF數5~10防止過分解按物理機理估計分量數端點拓延寬度2~3個極值點間距緩解端點效應信號越長可以越寬松包絡下限eps或1e-6防止除以零數值保護如果你發(fā)現分解結果對參數非常敏感建議先固定閾值和迭代次數只調節(jié)最大PF數減少調參維度。等流程跑通后再回來微調其他參數。5. 實戰(zhàn)演練滾動軸承故障信號的特征提取5.1 故障特征頻率與仿真信號滾動軸承外圈故障的特征頻率BPFO近似為BPFO (n_r / 2) · f_r · (1 - d/D · cosα)其中n_r是滾動體數量f_r是轉頻d是滾動體直徑D是節(jié)徑α是接觸角。為了演示方便我直接設定轉頻30Hz外圈故障特征頻率118Hz結構共振頻率2000Hz。仿真信號構造思路每個故障周期產生一次沖擊沖擊在結構共振頻率處引發(fā)衰減振蕩同時沖擊幅度受轉頻調制。最終再加上一點白噪聲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分量分別做包絡譜觀察哪個分量在故障特征頻率處能量最突出就能判斷故障沖擊主要調制在哪個頻帶上。這個方法比直接對原始信號做包絡譜更穩(wěn)健。因為LMD已經把最相關的調幅調頻分量從其他干擾中分離出來包絡譜的譜峰更干凈信噪比更高。5.3 實測中要注意的細節(jié)用LMD做軸承故障診斷有幾個細節(jié)非常重要采樣率要足夠高。LMD要提取沖擊激發(fā)的高頻共振如果采樣率只有1kHz而共振頻率在2kHz以上信號已經被嚴重混疊分解結果毫無意義。建議采樣率至少是最高關注頻率的5倍以上。信號長度要覆蓋足夠多的故障周期。如果只截取兩三個沖擊周期LMD分解和包絡譜都會因為樣本太少而失真。一般至少保證50個故障周期的長度。噪聲不能完全無視。LMD對強噪聲比較敏感實測信號通常比仿真信號臟得多。我的經驗是先做輕度帶通濾波把明顯無關的頻帶去掉再做LMD效果會好很多。但如果濾波帶寬本身選錯了又回到了傳統(tǒng)方法的玄學問題。所以濾波帶寬宜寬不宜窄只去掉極端高頻噪底和極低頻趨勢即可。6. LMD的適用范圍、改進方向與我的經驗6.1 什么時候用LMD什么時候換VMD/EMDLMD不是萬能的選型時要看信號特點。信號場景推薦算法原因強調幅調頻、故障沖擊明顯LMDPF分量物理意義清晰包絡解調方便多通道同步信號多元LMD或VMD保證各通道分量一致性無調制的平穩(wěn)疊加信號VMDLMD可能過度分解寬頻強噪聲VMD或EEMDLMD對噪聲敏感需要嚴格數學最優(yōu)VMD變分框架更規(guī)范我自己的習慣是EMD、LMD、VMD各跑一遍對比分解結果后選擇物理意義最清晰的那個。聽起來麻煩但實際多寫幾行腳本的事卻能讓結果分析可靠很多。6.2 值得嘗試的改進方向如果要把LMD用到科研項目里下面幾個方向值得探索第一用三次樣條插值替代滑動平均這是目前最實用的改進也是我在代碼里采用的方式。第二針對噪聲較強的情況引入集合平均思想類似EEMD那樣加入有限次白噪聲再取平均可以有效緩解模態(tài)混疊。第三開發(fā)自適應閾值選擇機制根據信號噪聲能量估計自動確定內層迭代停止閾值。第四推廣到多元LMD同時處理多通道振動信號避免各通道獨立分解導致的分量錯位問題。6.3 幾個容易被忽視的實操細節(jié)最后分享幾個我自己踩過的坑。第一次實現LMD時我給內層循環(huán)設了最大迭代次數但設得很大結果遇到一段幅值接近零的低能量信號迭代一兩百次都停不下來整個腳本卡死。后來把上限壓到200并結合一個殘差能量判斷問題才解決。現在我的代碼里既能按閾值收斂也能強制退出不會因為個別信號把整個分析流程卡掉。閾值方面我也做過對比實驗Δ取0.01時PF分量的包絡曲線會顯得比較毛糙瞬時幅值上能看到細微的鋸齒調到0.001后包絡光滑很多但迭代時間幾乎翻倍。對于一般工程診斷0.005是一個不錯的折中選擇如果追求精細分析再用0.001。極值點提取這一步看似簡單實際也容易出問題。比如信號帶有直流偏置或趨勢項時極值點的均勻性會變差分解出的第一個PF可能把趨勢當成分量。建議在分解前先去趨勢或做一次高通濾波把直流和超低頻趨勢去掉LMD的分解質量會明顯提升。說實話LMD在MATLAB里手寫并不難真正難的是參數選擇和端點處理很多論文復現不了往往就是這幾個細節(jié)沒處理好。如果你也遇到分解結果亂跑的怪現象先檢查端點拓延和閾值再檢查極值點提取是否太粗糙。調通之后拿來做故障特征提取還是相當順手的——至少在我處理軸承振動信號的經驗里它比EMD穩(wěn)定不少結果也更接近物理直覺。