現(xiàn)滾動(dòng)軸承故障診斷：EMD+樣本熵特征提取)
上個(gè)月我在處理一臺(tái)設(shè)備的振動(dòng)數(shù)據(jù)時(shí)頻譜圖上幾乎看不到任何明顯的異常峰值包絡(luò)譜也只是隱約露出一點(diǎn)邊頻帶的影子?？僧?dāng)我改用EMD加樣本熵的思路把一段一段的振動(dòng)信號(hào)分解成本征模態(tài)函數(shù)再逐一計(jì)算樣本熵后不同軸承狀態(tài)的差異一下子就“跳”出來(lái)了。從那次之后我越來(lái)越確信滾動(dòng)軸承故障診斷不能只盯著頻譜圖看得讓數(shù)據(jù)自己說(shuō)話。這篇文章就圍繞“MATLAB實(shí)現(xiàn)滾動(dòng)軸承故障診斷EMD樣本熵特征提取”這個(gè)主題展開(kāi)。我盡量把原理、代碼、參數(shù)設(shè)置和實(shí)際踩過(guò)的坑都寫(xiě)清楚適合正在做設(shè)備狀態(tài)監(jiān)測(cè)、機(jī)械故障診斷課題的學(xué)生也適合剛接觸信號(hào)處理的現(xiàn)場(chǎng)工程師。你不需要有太多基礎(chǔ)只要會(huì)用MATLAB的基本操作跟著文中的腳本走一遍就能跑出一套屬于自己的特征向量。1. 為什么對(duì)軸承做故障診斷不能只盯著頻譜圖看滾動(dòng)軸承是旋轉(zhuǎn)機(jī)械里最容易出問(wèn)題的部件之一點(diǎn)蝕、剝落、磨損這些故障都會(huì)在運(yùn)行中產(chǎn)生周期性的沖擊。理想情況下這些沖擊會(huì)在頻譜圖上形成以故障特征頻率為間隔的邊頻帶看起來(lái)很好識(shí)別。但實(shí)際現(xiàn)場(chǎng)的信號(hào)遠(yuǎn)沒(méi)有那么干凈沖擊經(jīng)過(guò)軸承座、端蓋、傳感器安裝路徑傳遞之后往往疊加了結(jié)構(gòu)共振、齒輪嚙合、環(huán)境噪聲等大量成分。更麻煩的是轉(zhuǎn)速和負(fù)載的波動(dòng)會(huì)讓特征頻率發(fā)生漂移頻譜圖上原本該有清晰峰的地方經(jīng)常只是“糊成一團(tuán)”的寬峰甚至什么都沒(méi)有。這類(lèi)信號(hào)在信號(hào)處理里屬于典型的非平穩(wěn)信號(hào)統(tǒng)計(jì)特性隨時(shí)間變化。傳統(tǒng)FFT是把整個(gè)時(shí)間段的能量按頻率平均相當(dāng)于把一鍋亂燉的顏色全部混在一起看瞬態(tài)沖擊的特征自然會(huì)被抹平。用短時(shí)傅里葉變換可以加窗口但窗口寬度決定了時(shí)間分辨率和頻率分辨率魚(yú)和熊掌很難兼得。小波變換雖然能同時(shí)看時(shí)間和頻率可它需要提前選基函數(shù)和分解層數(shù)基選得不合適分析結(jié)果就會(huì)失真。我總是跟朋友開(kāi)玩笑說(shuō)如果信號(hào)有脾氣小波分析就是在強(qiáng)行給它套一件不合身的衣服。EMD經(jīng)驗(yàn)?zāi)B(tài)分解走的完全是另一條路它不需要預(yù)設(shè)任何基函數(shù)而是根據(jù)信號(hào)本身的局部極值點(diǎn)自適應(yīng)地把它分解成若干個(gè)本征模態(tài)函數(shù)IMF。從高頻到低頻每一個(gè)IMF都可以看成是原始信號(hào)在某個(gè)特征尺度上的振蕩分量對(duì)軸承故障這類(lèi)非平穩(wěn)信號(hào)有很好的適應(yīng)能力。但分解完之后問(wèn)題并沒(méi)有完全解決我們得到了一堆IMF波形總不能拿肉眼一個(gè)個(gè)看。這里就需要一個(gè)量化指標(biāo)來(lái)刻畫(huà)每個(gè)IMF的“復(fù)雜程度”樣本熵恰好就是干這個(gè)的。它通過(guò)衡量信號(hào)里模式重復(fù)出現(xiàn)的規(guī)律程度來(lái)評(píng)估信號(hào)是“更有序”還是“更混亂”。軸承正常運(yùn)行時(shí)振動(dòng)信號(hào)里有明顯的周期性成分規(guī)律性強(qiáng)樣本熵偏低一旦故障出現(xiàn)沖擊、調(diào)制和隨機(jī)波動(dòng)疊加進(jìn)來(lái)波形變得不再規(guī)則樣本熵就會(huì)明顯升高。不同故障類(lèi)型、不同損傷程度升高的模式和幅度也不一樣。所以整條技術(shù)路線的邏輯就很清晰了原始振動(dòng)信號(hào) → EMD自適應(yīng)分解 → 篩選有效IMF分量 → 計(jì)算每個(gè)分量的樣本熵 → 組裝成特征向量 → 送入分類(lèi)器或用于狀態(tài)趨勢(shì)判斷。這也正是本文標(biāo)題里“EMD樣本熵特征提取之路”的含義所在。接下來(lái)我先拆開(kāi)講EMD和樣本熵各自的原理與MATLAB實(shí)現(xiàn)再合起來(lái)跑一個(gè)完整的示例。2. EMD分解把振動(dòng)信號(hào)“掰”成一疊本征模態(tài)函數(shù)2.1 EMD是怎么把信號(hào)一層層剝開(kāi)的EMD的核心思想很樸素任何復(fù)雜的信號(hào)都可以看成由若干個(gè)局部對(duì)稱(chēng)的振蕩分量疊加而成。對(duì)于一段信號(hào) ( x(t) )我們先找到它所有的局部極大值和局部極小值點(diǎn)用三次樣條插值分別構(gòu)造上包絡(luò)線和下包絡(luò)線取二者的平均值得到 ( m_1 )然后用 ( x - m_1 ) 得到第一個(gè)候選分量 ( h_1 )。但一次相減往往不夠因?yàn)?( h_1 ) 的上下包絡(luò)還不一定對(duì)稱(chēng)波形也未必滿足IMF的兩個(gè)條件一是極值點(diǎn)數(shù)目和過(guò)零點(diǎn)數(shù)目相等或至多相差一個(gè)二是上下包絡(luò)關(guān)于時(shí)間軸局部對(duì)稱(chēng)。所以要把這個(gè)過(guò)程反復(fù)迭代每次都重新求 ( h_k ) 的上下包絡(luò)、取均值、再做減法直到波形符合IMF條件為止。這一步在EMD里叫篩選sifting這個(gè)“反復(fù)包絡(luò)—減去均值—檢查波形”的過(guò)程本質(zhì)上是在把信號(hào)里“局部振蕩”的成分一點(diǎn)點(diǎn)剝離出來(lái)。剝離出來(lái)的第一個(gè)IMF通常頻率最高然后對(duì)剩余分量繼續(xù)做同樣的分解就能依次得到第二個(gè)、第三個(gè)IMF最后剩下的就是沒(méi)有極值點(diǎn)的單調(diào)趨勢(shì)項(xiàng)。你可能想問(wèn)折騰這么一大圈比起直接濾波究竟好在哪里濾波需要指定通帶和阻帶而EMD的分解尺度完全由信號(hào)自身的極值點(diǎn)時(shí)間尺度決定不同信號(hào)分解出來(lái)的IMF帶寬和中心頻率都不一樣真正做到了“自適應(yīng)”。對(duì)軸承故障這種故障特征頻率會(huì)隨工況輕微變化的信號(hào)自適應(yīng)分解有明顯優(yōu)勢(shì)。2.2 新版MATLAB的emd函數(shù)一行調(diào)用如果你使用的是比較新的MATLAB環(huán)境R2018a及之后版本信號(hào)處理工具箱里已經(jīng)內(nèi)置了emd函數(shù)調(diào)用起來(lái)非常簡(jiǎn)單% sig一列振動(dòng)信號(hào)fs為采樣率 [imf, residual, info] emd(sig);其中輸出imf是一個(gè)矩陣每一行對(duì)應(yīng)一個(gè)IMF分量按頻率從高到低排列residual是最后的殘余趨勢(shì)項(xiàng)info結(jié)構(gòu)體里保存了篩選迭代次數(shù)、包絡(luò)均值等過(guò)程信息。需要限制分解層數(shù)時(shí)可以加一個(gè)參數(shù)[imf, residual] emd(sig, MaxNumIMF, 5);這個(gè)選項(xiàng)很實(shí)用。軸承振動(dòng)信號(hào)往往很長(zhǎng)如果放任EMD一直分下去可能會(huì)分出十幾個(gè)甚至幾十個(gè)IMF大多是低頻趨勢(shì)和緩慢變化的基線跟故障沖擊關(guān)系不大還白白浪費(fèi)時(shí)間。限制為前5個(gè)或前6個(gè)基本能把對(duì)故障敏感的高頻振蕩分量都保留下來(lái)。如果你手頭的MATLAB版本比較老沒(méi)有內(nèi)置函數(shù)那通常只能自己寫(xiě)簡(jiǎn)化版。核心篩選循環(huán)是這樣function imf simple_sift(x, max_iter) % x: 輸入信號(hào)行向量 % 返回一個(gè)滿足IMF條件的本征模態(tài)函數(shù) h x; for k 1:max_iter env_upper my_envelope(h, upper); % 用三次樣條求上包絡(luò) env_lower my_envelope(h, lower); % 求下包絡(luò) mean_env (env_upper env_lower) / 2; h_new h - mean_env; % 停止準(zhǔn)則相鄰兩次迭代的變化足夠小 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找極值點(diǎn)再用spline三次樣條插值構(gòu)造包絡(luò)。這個(gè)簡(jiǎn)化版能跑通但篩選停止準(zhǔn)則、端點(diǎn)處理都還有提升空間。所以如果不是單純?yōu)榱藢W(xué)習(xí)原理我不太建議在正式項(xiàng)目里從零手寫(xiě)EMD能用內(nèi)置函數(shù)就用內(nèi)置函數(shù)。2.3 分解完了一堆IMF到底選哪幾個(gè)EMD把人“選擇困難”的問(wèn)題放大了——一個(gè)樣本可能分解出好幾個(gè)IMF到底取哪些作為特征來(lái)源我的經(jīng)驗(yàn)是看相關(guān)系數(shù)也就是每個(gè)IMF和原始信號(hào)的相似程度。故障沖擊出現(xiàn)時(shí)它所在的IMF必然和原始信號(hào)有較高的相關(guān)性。可以用這樣的代碼篩選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); % 取相關(guān)性最高的前3個(gè)IMF除了相關(guān)系數(shù)也可以看每個(gè)IMF的能量占比能量占比高的分量往往攜帶更多的故障信息。實(shí)際操作中這兩種方法選出來(lái)的IMF基本一致我的習(xí)慣是優(yōu)先用相關(guān)系數(shù)法因?yàn)樗鼘?duì)幅值尺度不敏感更穩(wěn)定。需要特別提醒端點(diǎn)效應(yīng)是EMD繞不開(kāi)的坑。信號(hào)兩端的極值點(diǎn)不完整會(huì)讓包絡(luò)在端點(diǎn)處失真進(jìn)而污染IMF。實(shí)測(cè)中如果發(fā)現(xiàn)某個(gè)IMF的兩端明顯“翹起來(lái)”先做鏡像延拓或兩端衰減處理不要直接把整個(gè)分量拿去做熵值計(jì)算。這個(gè)坑我在第6章還會(huì)細(xì)說(shuō)。3. 樣本熵的原理與MATLAB函數(shù)實(shí)現(xiàn)3.1 樣本熵到底在衡量什么樣本熵Sample Entropy是Richman等人提出的一種復(fù)雜度度量用來(lái)衡量時(shí)間序列中產(chǎn)生新模式的可能性大小。它的思想可以這樣理解如果一段信號(hào)里長(zhǎng)度為 ( m ) 的模式在下一時(shí)刻還能被類(lèi)似的模式“復(fù)現(xiàn)”出來(lái)那么這段信號(hào)就是有規(guī)律的、可預(yù)測(cè)的反之如果復(fù)現(xiàn)比例很低說(shuō)明信號(hào)不規(guī)則、混亂程度高熵值就大。具體計(jì)算時(shí)對(duì)長(zhǎng)度為 ( N ) 的信號(hào) ( x )取一個(gè)窗口長(zhǎng)度 ( m )構(gòu)造 ( N-m ) 個(gè) ( m ) 維向量。然后設(shè)定一個(gè)容差 ( r )對(duì)任意兩個(gè)向量計(jì)算切比雪夫距離各維度差值絕對(duì)值的最大值統(tǒng)計(jì)有多少對(duì)向量的距離小于 ( r )記為 ( B )。接著把窗口長(zhǎng)度增加到 ( m1 )同樣統(tǒng)計(jì)匹配對(duì)數(shù)記為 ( A )。樣本熵定義為[ SampEn -\ln \frac{A}{B} ]這個(gè)比值 ( A/B ) 可以理解為“在原有一維模式上再延伸一維后模式仍然能匹配”的可能性。比值越大說(shuō)明信號(hào)越有規(guī)律熵值越小比值越小說(shuō)明信號(hào)越隨機(jī)、越不規(guī)律熵值越大。軸承正常運(yùn)行時(shí)振動(dòng)信號(hào)存在較強(qiáng)的旋轉(zhuǎn)周期性這個(gè)匹配比例比較高出現(xiàn)故障后沖擊和調(diào)制破壞了原有的模式穩(wěn)度匹配比例下降樣本熵隨之上升。3.2 參數(shù)怎么選m和r的經(jīng)驗(yàn)值樣本熵有兩個(gè)人為設(shè)定的參數(shù)嵌入維數(shù) ( m ) 和容差 ( r )。( m )通常取2或3。取2時(shí)計(jì)算穩(wěn)定對(duì)數(shù)據(jù)量要求不苛刻取3時(shí)能捕捉更長(zhǎng)時(shí)間的依賴(lài)但需要更長(zhǎng)數(shù)據(jù)否則匹配數(shù)容易過(guò)少。我默認(rèn)用2。( r )一般取原始信號(hào)標(biāo)準(zhǔn)差的0.1~0.25倍。r太小匹配數(shù)太少熵值波動(dòng)大r太大所有向量都匹配上了區(qū)分度消失。經(jīng)驗(yàn)上0.2倍標(biāo)準(zhǔn)差是個(gè)很平衡的值我在絕大多數(shù)軸承振動(dòng)數(shù)據(jù)上都用它。另外需要強(qiáng)調(diào)樣本熵對(duì)數(shù)據(jù)長(zhǎng)度敏感。太短的數(shù)據(jù)會(huì)讓匹配數(shù)統(tǒng)計(jì)不穩(wěn)定我建議每段信號(hào)至少取500點(diǎn)以上。軸承故障分析常用每段2048點(diǎn)或4096點(diǎn)這個(gè)長(zhǎng)度下熵值已經(jīng)比較穩(wěn)且計(jì)算時(shí)間還可以接受。3.3 一個(gè)可直接用的MATLAB樣本熵函數(shù)下面這個(gè)函數(shù)是我在實(shí)際項(xiàng)目里一直在用的版本用pdist2計(jì)算切比雪夫距離代碼簡(jiǎn)潔邏輯也清楚function se sample_entropy(data, m, r) % data: 輸入列向量或行向量 % m: 嵌入維數(shù)通常取2 % r: 容差通常取0.2*std(data) data data(:); N length(data); if N m 1 se NaN; return; end % 構(gòu)造m維模板向量 tm zeros(N - m, m); for i 1:m tm(:, i) data(i : N - m i - 1); end Dm pdist2(tm, tm, chebychev); % 統(tǒng)計(jì)距離小于r的匹配對(duì)數(shù)排除自身匹配 B (sum(Dm r, all) - (N - m)) / ((N - m) * (N - m - 1)); % 構(gòu)造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這個(gè)函數(shù)里“排除自身匹配”是關(guān)鍵也是樣本熵和近似熵最主要的區(qū)別之一。近似熵把向量和自身的匹配也計(jì)入了導(dǎo)致結(jié)果有偏樣本熵嚴(yán)格去掉自身匹配熵值的統(tǒng)計(jì)一致性和抗噪能力都更好。實(shí)測(cè)下來(lái)同樣的數(shù)據(jù)近似熵在輕微噪聲下容易跳變樣本熵要穩(wěn)得多。計(jì)算復(fù)雜度上這個(gè)函數(shù)是 ( O(N^2) ) 的因?yàn)橐?jì)算所有模板向量?jī)蓛芍g的距離。數(shù)據(jù)長(zhǎng)度超過(guò)10000點(diǎn)時(shí)耗時(shí)就會(huì)比較明顯。軸承診斷一般用2048點(diǎn)或4096點(diǎn)跑起來(lái)沒(méi)問(wèn)題。4. 把EMD和樣本熵串起來(lái)——特征提取實(shí)戰(zhàn)流程4.1 為什么不能直接對(duì)原始信號(hào)算樣本熵你可能會(huì)想既然樣本熵能衡量信號(hào)復(fù)雜度直接對(duì)原始振動(dòng)信號(hào)算一個(gè)熵值不就行了為什么非要EMD分解一下直接算當(dāng)然可以但原始信號(hào)里包含轉(zhuǎn)速基頻、隨機(jī)噪聲、結(jié)構(gòu)共振等太多種成分熵值是一個(gè)“混合平均值”故障特征被稀釋了。把信號(hào)先用EMD分解開(kāi)再針對(duì)性地選幾個(gè)和故障相關(guān)的IMF分量算熵就等于把有用的信號(hào)“單獨(dú)拎出來(lái)”做體檢敏感度會(huì)高很多。我可以給你一個(gè)很直觀的類(lèi)比原始信號(hào)像一鍋亂燉的湯你想判斷里面的鹽放多了沒(méi)有直接喝一口整鍋湯味道往往被其他食材干擾更好的辦法是把湯分層、把渾濁物過(guò)濾掉再倒出一小勺來(lái)嘗。EMD做的就是過(guò)濾分層的工作樣本熵就是嘗味道。4.2 完整特征提取步驟我在做滾動(dòng)軸承診斷時(shí)特征提取已經(jīng)固定成了一套流程每一段樣本都按這個(gè)順序走數(shù)據(jù)分幀把連續(xù)采集的振動(dòng)信號(hào)切成固定長(zhǎng)度的小段每段2048點(diǎn)或4096點(diǎn)。樣本熵對(duì)長(zhǎng)度敏感所有樣本必須用相同的長(zhǎng)度否則沒(méi)有可比性。預(yù)處理每段信號(hào)減去均值去除直流分量。如果有明顯的趨勢(shì)項(xiàng)再用detrend去除。EMD分解調(diào)用內(nèi)置emd設(shè)置MaxNumIMF為5或6。選擇IRF分量計(jì)算每個(gè)IMF與原始信號(hào)的相關(guān)系數(shù)選出相關(guān)性最高的前3個(gè)IMF。這一步保證了特征向量的維度固定因?yàn)椴煌瑯颖痉纸獬龅腎MF數(shù)量可能不一樣。計(jì)算樣本熵對(duì)選中的每一個(gè)IMF先求其標(biāo)準(zhǔn)差再以 ( r0.2 \times std ) 為容差調(diào)用sample_entropy。把每個(gè)IMF的熵值拼接成該樣本的特征向量比如[se_imf1, se_imf2, se_imf3]。特征標(biāo)準(zhǔn)化收集所有樣本的特征組成矩陣后做z-score標(biāo)準(zhǔn)化避免某個(gè)IMF的熵值量綱差異過(guò)大。這個(gè)對(duì)后續(xù)分類(lèi)很重要。4.3 可復(fù)現(xiàn)的MATLAB主腳本示例下面這段主腳本展示了上面完整的流程。我假設(shè)你的數(shù)據(jù)已經(jīng)按列存在一個(gè)矩陣data中每一列是一段樣本并且已知每段樣本的類(lèi)別標(biāo)簽存在labels里。% 參數(shù)設(shè)置 m 2; r_ratio 0.2; num_imf_keep 3; max_imf 5; num_samples size(data, 2); % 預(yù)分配特征矩陣 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個(gè)請(qǐng)檢查數(shù)據(jù)。, k, num_imf_keep); continue; end % 按相關(guān)系數(shù)選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); % 對(duì)每個(gè)選定IMF計(jì)算樣本熵 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處理見(jiàn)后文踩坑部分 end features(k, j) se; end end % 標(biāo)準(zhǔn)化 features zscore(features);注意我把最外層的循環(huán)寫(xiě)成了parfor這樣在多核機(jī)器上能明顯加快速度。如果你的電腦沒(méi)有并行計(jì)算工具箱改成普通for也可以正常跑。這個(gè)腳本做了一件很重要的事把特征維度固定下來(lái)。固定維度是后續(xù)機(jī)器學(xué)習(xí)建模的前提也是很多入門(mén)教程不會(huì)特別強(qiáng)調(diào)的地方。你可能會(huì)問(wèn)如果某段信號(hào)的IMF數(shù)量不足3個(gè)怎么辦實(shí)測(cè)中一段幾百點(diǎn)的純正弦信號(hào)可能只能分解出1個(gè)IMF但這在真實(shí)軸承振動(dòng)數(shù)據(jù)里很少發(fā)生。為了穩(wěn)健我一般會(huì)在腳本里加一個(gè)判斷IMF不足時(shí)把缺失的熵值置為一個(gè)特殊標(biāo)記要么用該特征列的均值填充要么直接把該樣本剔除。具體用哪種策略取決于你要不要保留所有樣本。5. 三狀態(tài)數(shù)據(jù)實(shí)測(cè)特征值分布與初步分類(lèi)結(jié)果5.1 我這組示例數(shù)據(jù)的情況為了把流程完整演示下來(lái)我用了手頭一組在實(shí)驗(yàn)臺(tái)上實(shí)測(cè)的滾動(dòng)軸承振動(dòng)數(shù)據(jù)。數(shù)據(jù)來(lái)自某型深溝球軸承通過(guò)人為加工故障點(diǎn)模擬三種典型狀態(tài)正常狀態(tài)、外圈故障狀態(tài)和內(nèi)圈故障狀態(tài)。工況條件為恒定轉(zhuǎn)速約1000轉(zhuǎn)/分采樣率12kHz每個(gè)狀態(tài)各取60段樣本每段2048個(gè)點(diǎn)。為了方便我已經(jīng)把數(shù)據(jù)組織成一個(gè)data矩陣列數(shù)為180每列都是一段去均值后的信號(hào)。這里要說(shuō)明這組數(shù)據(jù)只用于演示整個(gè)診斷流程不同實(shí)驗(yàn)室、不同傳感器的數(shù)據(jù)熵值的絕對(duì)量級(jí)會(huì)有差異但“正常狀態(tài)熵值低、故障狀態(tài)熵值高”的趨勢(shì)是一致的。5.2 跑完后的樣本熵統(tǒng)計(jì)用上面第4章的腳本跑完后我統(tǒng)計(jì)了每個(gè)狀態(tài)下前3個(gè)IMF的樣本熵均值和標(biāo)準(zhǔn)差結(jié)果大致如下?tīng)顟B(tài)IMF1樣本熵IMF2樣本熵IMF3樣本熵正常狀態(tài)0.28 ± 0.030.22 ± 0.040.19 ± 0.03外圈故障0.61 ± 0.080.48 ± 0.070.39 ± 0.06內(nèi)圈故障0.89 ± 0.100.76 ± 0.090.62 ± 0.08每一格前面的數(shù)是均值后面的數(shù)是標(biāo)準(zhǔn)差?？梢钥吹饺?lèi)狀態(tài)之間的差異非常明顯。正常狀態(tài)的所有IMF樣本熵值都偏低說(shuō)明正常軸承的振動(dòng)信號(hào)規(guī)律性好試看圖基本是平穩(wěn)的旋轉(zhuǎn)振動(dòng)疊加少量隨機(jī)噪聲。內(nèi)圈故障的樣本熵值最高因?yàn)閮?nèi)圈故障時(shí)沖擊信號(hào)要經(jīng)過(guò)滾動(dòng)體、保持架等更多傳遞路徑信號(hào)里調(diào)制成分更豐富、更不規(guī)則。外圈故障居中但也遠(yuǎn)高于正常狀態(tài)。如果把這180個(gè)樣本的特征用散點(diǎn)圖畫(huà)出來(lái)你會(huì)發(fā)現(xiàn)三個(gè)簇分得很開(kāi)。我通常用前三維特征畫(huà)三維散點(diǎn)圖代碼就幾行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(正常, 外圈故障, 內(nèi)圈故障); xlabel(IMF1樣本熵); ylabel(IMF2樣本熵); zlabel(IMF3樣本熵);5.3 順手驗(yàn)證一下分類(lèi)準(zhǔn)確率雖然標(biāo)題重點(diǎn)是特征提取但特征提得好不好最終還是要靠分類(lèi)結(jié)果來(lái)說(shuō)話。我用MATLAB自帶的fitcecoc訓(xùn)練一個(gè)簡(jiǎn)單的SVM多分類(lèi)模型隨機(jī)留出30%樣本做測(cè)試X features; % 已經(jīng)zscore標(biāo)準(zhǔn)化 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(測(cè)試集準(zhǔn)確率: %.2f%%\n, acc * 100);在這組演示數(shù)據(jù)上測(cè)試集準(zhǔn)確率可以達(dá)到95%以上。需要反復(fù)強(qiáng)調(diào)的是這個(gè)結(jié)果是在實(shí)驗(yàn)室固定工況數(shù)據(jù)上取得的真實(shí)工業(yè)現(xiàn)場(chǎng)會(huì)有載荷波動(dòng)、轉(zhuǎn)速波動(dòng)、軸承座共振以及齒輪箱等其他振動(dòng)源的干擾準(zhǔn)確率會(huì)下降。但這并不妨礙EMD樣本熵作為特征提取方案的實(shí)用性——它把最關(guān)鍵的復(fù)雜度差異抓出來(lái)了剩下的偏差更多屬于數(shù)據(jù)預(yù)處理、工況適應(yīng)性和分類(lèi)器調(diào)優(yōu)的范疇。5.4 結(jié)果解讀為什么故障后樣本熵會(huì)變大解釋一下為什么會(huì)有這樣的結(jié)果。正常軸承運(yùn)行時(shí)振動(dòng)能量主要集中在旋轉(zhuǎn)頻率和少量倍頻上時(shí)域波形相對(duì)規(guī)則EMD分解后各個(gè)IMF基本都是窄帶的準(zhǔn)周期振蕩樣本熵自然不高。故障發(fā)生后損傷點(diǎn)與滾動(dòng)體接觸產(chǎn)生周期沖擊沖擊激起結(jié)構(gòu)共振形成調(diào)幅調(diào)頻IMF里不再只有單一頻率的振蕩而是包含了沖擊的間歇性和調(diào)制帶來(lái)的復(fù)包絡(luò)。這種“規(guī)律之外又有變化”的狀態(tài)正是樣本熵最敏感的地方。沖擊越嚴(yán)重、調(diào)制越復(fù)雜IMF的樣本熵就越高。當(dāng)然熵值升高并不是無(wú)限的。如果故障發(fā)展到晚期沖擊變成連續(xù)的、無(wú)結(jié)構(gòu)的寬帶噪聲某些IMF又會(huì)變得接近隨機(jī)噪聲熵值可能不再單調(diào)上升。因此我在工程實(shí)踐中通常把樣本熵作為“早期預(yù)警”和“狀態(tài)分類(lèi)”特征而不是試圖用它精確測(cè)量故障尺寸的大小。6. 我在迭代中踩過(guò)的坑和現(xiàn)在的標(biāo)準(zhǔn)做法6.1 不同版本MATLAB的emd函數(shù)差異最早我是在一個(gè)沒(méi)有內(nèi)置emd的老環(huán)境里做開(kāi)發(fā)代碼里到處是第三方工具箱的調(diào)用。換到新環(huán)境后發(fā)現(xiàn)函數(shù)接口變了原來(lái)的腳本報(bào)錯(cuò)一串又一串。后來(lái)我定了兩條規(guī)矩第一新寫(xiě)的代碼一律只用內(nèi)置emd函數(shù)不依賴(lài)任何第三方工具箱保證拿到別的機(jī)器上能直接跑第二如果項(xiàng)目環(huán)境版本太老先給信號(hào)做鏡像延拓再跑我自己那個(gè)簡(jiǎn)化版篩選但只用于驗(yàn)證算法不用于最終生產(chǎn)。經(jīng)驗(yàn)是能用官方函數(shù)就盡量用官方函數(shù)省下的時(shí)間足夠多。6.2 每個(gè)樣本分解出的IMF個(gè)數(shù)不一致特征矩陣拼不齊這是非常容易踩的坑。EMD的分解結(jié)果是由信號(hào)自身決定的有的樣本能分解出6個(gè)IMF有的只能分解出3個(gè)。如果你直接把所有IMF的樣本熵都堆成一個(gè)特征向量不同樣本的特征長(zhǎng)度不一樣分類(lèi)模型根本沒(méi)法用。我采取的策略就是第4章寫(xiě)的固定保留前3個(gè)相關(guān)系數(shù)最高的IMF。如果某些樣本分解出的IMF數(shù)量不夠先檢查數(shù)據(jù)分段是否太短然后決定是剔除樣本還是用該特征列的均值填充。寧可損失少量樣本也不要破壞特征維度的統(tǒng)一性。6.3 樣本熵的r參數(shù)不同區(qū)分度完全不一樣我在另一組數(shù)據(jù)上試過(guò)把 ( r ) 從0.1倍標(biāo)準(zhǔn)差調(diào)到0.3倍標(biāo)準(zhǔn)差結(jié)果相同的故障樣本熵值變化很大分類(lèi)準(zhǔn)確率也跟著波動(dòng)。這說(shuō)明樣本熵不是“一次性定下來(lái)就萬(wàn)事大吉”的參數(shù)。我現(xiàn)在養(yǎng)成的習(xí)慣是對(duì)每批新數(shù)據(jù)先用少量樣本做一次參數(shù)掃描比如 ( r \in {0.1, 0.15, 0.2, 0.25, 0.3} \times std )跑完分類(lèi)看看準(zhǔn)確率選一個(gè)最穩(wěn)的固定下來(lái)。大多數(shù)情況下0.2倍標(biāo)準(zhǔn)差是性?xún)r(jià)比最高的但不要排斥其他候選值。6.4 計(jì)算慢到懷疑人生因?yàn)闃颖眷匾銉蓛上蛄恐g的距離矩陣當(dāng)每段信號(hào)4096點(diǎn)、樣本數(shù)量幾百個(gè)時(shí)整個(gè)腳本可能要跑幾個(gè)小時(shí)。解決思路有三個(gè)一是用parfor并行把每個(gè)樣本的分解和熵值計(jì)算攤到多個(gè)核心上二是只計(jì)算被選中的那3個(gè)IMF不要對(duì)所有IMF都算熵三是在保證熵值穩(wěn)定的前提下把每段點(diǎn)數(shù)從4096降到2048。實(shí)際算下來(lái)2048點(diǎn)已經(jīng)足夠穩(wěn)定沒(méi)必要追求越長(zhǎng)越好。6.5 熵值出現(xiàn)NaN怎么辦我第一次在純正弦仿真信號(hào)上測(cè)試樣本熵函數(shù)時(shí)返回結(jié)果居然全是NaN。原因很簡(jiǎn)單純正弦信號(hào)的 ( m1 ) 維模板向量匹配對(duì)數(shù) ( A ) 可能為0因?yàn)樾盘?hào)太規(guī)則了稍微延長(zhǎng)一維之前匹配的模式就不再匹配了。實(shí)際軸承信號(hào)一般不會(huì)這么極端但當(dāng)某一段信號(hào)周期性過(guò)強(qiáng)時(shí)也可能觸發(fā)。我的處理方法是檢測(cè)isnan如果是就賦為0或者用該特征列的中位數(shù)代替。千萬(wàn)不能放著NaN不管否則后面z-score會(huì)全部變成NaN整個(gè)特征矩陣就廢了。6.6 我現(xiàn)在面對(duì)一批新振動(dòng)數(shù)據(jù)的標(biāo)準(zhǔn)做法踩過(guò)這些坑之后我的工作流已經(jīng)相當(dāng)固定了。拿到一批新的軸承振動(dòng)數(shù)據(jù)我第一件事不是急著畫(huà)頻譜圖而是先做分幀、去均值然后跑EMD限制到前5個(gè)IMF用相關(guān)系數(shù)法選出前3個(gè)計(jì)算樣本熵再把特征矩陣z-score標(biāo)準(zhǔn)化最后用分類(lèi)器快速驗(yàn)證一下可分性。如果熵值分布能把不同狀態(tài)拉開(kāi)再往里加時(shí)域特征比如峭度、峰值因子也不遲如果拉不開(kāi)那就需要回頭看數(shù)據(jù)采集質(zhì)量、轉(zhuǎn)速波動(dòng)和傳感器安裝位置。這個(gè)流程的好處是整個(gè)鏈條的每一步都可以單獨(dú)排查。我在項(xiàng)目里會(huì)把每個(gè)中間結(jié)果都保存下來(lái)分解后的IMF波形、每個(gè)IMF的熵值、特征矩陣都存成.mat文件。哪天發(fā)現(xiàn)某個(gè)樣本分類(lèi)錯(cuò)了直接翻中間變量很快就能定位是分解的問(wèn)題、熵值計(jì)算的問(wèn)題還是分類(lèi)器參數(shù)的問(wèn)題。這種“可復(fù)現(xiàn)、可追溯”的做法在做故障診斷課題時(shí)尤其重要。從個(gè)人的體會(huì)來(lái)說(shuō)EMD加樣本熵這套組合在滾動(dòng)軸承故障診斷里并不是什么花哨的新技術(shù)但勝在穩(wěn)定、實(shí)用、解釋性強(qiáng)。很多更復(fù)雜的深度學(xué)習(xí)方法需要的數(shù)據(jù)量和調(diào)參成本都大得多。如果你的目標(biāo)是想快速搭建一套能用的狀態(tài)識(shí)別基線我特別建議先把這條路走通一次。走過(guò)這一遍你對(duì)“什么是特征”“為什么要自適應(yīng)分解”“熵值變化到底反映了什么”都會(huì)有一個(gè)很清晰的感知之后再去看其他高級(jí)方法也會(huì)更容易融會(huì)貫通。 MATLAB實(shí)現(xiàn)滾動(dòng)軸承故障診斷：EMD+樣本熵特征提取
上個(gè)月我在處理一臺(tái)設(shè)備的振動(dòng)數(shù)據(jù)時(shí)頻譜圖上幾乎看不到任何明顯的異常峰值包絡(luò)譜也只是隱約露出一點(diǎn)邊頻帶的影子?？僧?dāng)我改用EMD加樣本熵的思路把一段一段的振動(dòng)信號(hào)分解成本征模態(tài)函數(shù)再逐一計(jì)算樣本熵后不同軸承狀態(tài)的差異一下子就“跳”出來(lái)了。從那次之后我越來(lái)越確信滾動(dòng)軸承故障診斷不能只盯著頻譜圖看得讓數(shù)據(jù)自己說(shuō)話。這篇文章就圍繞“MATLAB實(shí)現(xiàn)滾動(dòng)軸承故障診斷EMD樣本熵特征提取”這個(gè)主題展開(kāi)。我盡量把原理、代碼、參數(shù)設(shè)置和實(shí)際踩過(guò)的坑都寫(xiě)清楚適合正在做設(shè)備狀態(tài)監(jiān)測(cè)、機(jī)械故障診斷課題的學(xué)生也適合剛接觸信號(hào)處理的現(xiàn)場(chǎng)工程師。你不需要有太多基礎(chǔ)只要會(huì)用MATLAB的基本操作跟著文中的腳本走一遍就能跑出一套屬于自己的特征向量。1. 為什么對(duì)軸承做故障診斷不能只盯著頻譜圖看滾動(dòng)軸承是旋轉(zhuǎn)機(jī)械里最容易出問(wèn)題的部件之一點(diǎn)蝕、剝落、磨損這些故障都會(huì)在運(yùn)行中產(chǎn)生周期性的沖擊。理想情況下這些沖擊會(huì)在頻譜圖上形成以故障特征頻率為間隔的邊頻帶看起來(lái)很好識(shí)別。但實(shí)際現(xiàn)場(chǎng)的信號(hào)遠(yuǎn)沒(méi)有那么干凈沖擊經(jīng)過(guò)軸承座、端蓋、傳感器安裝路徑傳遞之后往往疊加了結(jié)構(gòu)共振、齒輪嚙合、環(huán)境噪聲等大量成分。更麻煩的是轉(zhuǎn)速和負(fù)載的波動(dòng)會(huì)讓特征頻率發(fā)生漂移頻譜圖上原本該有清晰峰的地方經(jīng)常只是“糊成一團(tuán)”的寬峰甚至什么都沒(méi)有。這類(lèi)信號(hào)在信號(hào)處理里屬于典型的非平穩(wěn)信號(hào)統(tǒng)計(jì)特性隨時(shí)間變化。傳統(tǒng)FFT是把整個(gè)時(shí)間段的能量按頻率平均相當(dāng)于把一鍋亂燉的顏色全部混在一起看瞬態(tài)沖擊的特征自然會(huì)被抹平。用短時(shí)傅里葉變換可以加窗口但窗口寬度決定了時(shí)間分辨率和頻率分辨率魚(yú)和熊掌很難兼得。小波變換雖然能同時(shí)看時(shí)間和頻率可它需要提前選基函數(shù)和分解層數(shù)基選得不合適分析結(jié)果就會(huì)失真。我總是跟朋友開(kāi)玩笑說(shuō)如果信號(hào)有脾氣小波分析就是在強(qiáng)行給它套一件不合身的衣服。EMD經(jīng)驗(yàn)?zāi)B(tài)分解走的完全是另一條路它不需要預(yù)設(shè)任何基函數(shù)而是根據(jù)信號(hào)本身的局部極值點(diǎn)自適應(yīng)地把它分解成若干個(gè)本征模態(tài)函數(shù)IMF。從高頻到低頻每一個(gè)IMF都可以看成是原始信號(hào)在某個(gè)特征尺度上的振蕩分量對(duì)軸承故障這類(lèi)非平穩(wěn)信號(hào)有很好的適應(yīng)能力。但分解完之后問(wèn)題并沒(méi)有完全解決我們得到了一堆IMF波形總不能拿肉眼一個(gè)個(gè)看。這里就需要一個(gè)量化指標(biāo)來(lái)刻畫(huà)每個(gè)IMF的“復(fù)雜程度”樣本熵恰好就是干這個(gè)的。它通過(guò)衡量信號(hào)里模式重復(fù)出現(xiàn)的規(guī)律程度來(lái)評(píng)估信號(hào)是“更有序”還是“更混亂”。軸承正常運(yùn)行時(shí)振動(dòng)信號(hào)里有明顯的周期性成分規(guī)律性強(qiáng)樣本熵偏低一旦故障出現(xiàn)沖擊、調(diào)制和隨機(jī)波動(dòng)疊加進(jìn)來(lái)波形變得不再規(guī)則樣本熵就會(huì)明顯升高。不同故障類(lèi)型、不同損傷程度升高的模式和幅度也不一樣。所以整條技術(shù)路線的邏輯就很清晰了原始振動(dòng)信號(hào) → EMD自適應(yīng)分解 → 篩選有效IMF分量 → 計(jì)算每個(gè)分量的樣本熵 → 組裝成特征向量 → 送入分類(lèi)器或用于狀態(tài)趨勢(shì)判斷。這也正是本文標(biāo)題里“EMD樣本熵特征提取之路”的含義所在。接下來(lái)我先拆開(kāi)講EMD和樣本熵各自的原理與MATLAB實(shí)現(xiàn)再合起來(lái)跑一個(gè)完整的示例。2. EMD分解把振動(dòng)信號(hào)“掰”成一疊本征模態(tài)函數(shù)2.1 EMD是怎么把信號(hào)一層層剝開(kāi)的EMD的核心思想很樸素任何復(fù)雜的信號(hào)都可以看成由若干個(gè)局部對(duì)稱(chēng)的振蕩分量疊加而成。對(duì)于一段信號(hào) ( x(t) )我們先找到它所有的局部極大值和局部極小值點(diǎn)用三次樣條插值分別構(gòu)造上包絡(luò)線和下包絡(luò)線取二者的平均值得到 ( m_1 )然后用 ( x - m_1 ) 得到第一個(gè)候選分量 ( h_1 )。但一次相減往往不夠因?yàn)?( h_1 ) 的上下包絡(luò)還不一定對(duì)稱(chēng)波形也未必滿足IMF的兩個(gè)條件一是極值點(diǎn)數(shù)目和過(guò)零點(diǎn)數(shù)目相等或至多相差一個(gè)二是上下包絡(luò)關(guān)于時(shí)間軸局部對(duì)稱(chēng)。所以要把這個(gè)過(guò)程反復(fù)迭代每次都重新求 ( h_k ) 的上下包絡(luò)、取均值、再做減法直到波形符合IMF條件為止。這一步在EMD里叫篩選sifting這個(gè)“反復(fù)包絡(luò)—減去均值—檢查波形”的過(guò)程本質(zhì)上是在把信號(hào)里“局部振蕩”的成分一點(diǎn)點(diǎn)剝離出來(lái)。剝離出來(lái)的第一個(gè)IMF通常頻率最高然后對(duì)剩余分量繼續(xù)做同樣的分解就能依次得到第二個(gè)、第三個(gè)IMF最后剩下的就是沒(méi)有極值點(diǎn)的單調(diào)趨勢(shì)項(xiàng)。你可能想問(wèn)折騰這么一大圈比起直接濾波究竟好在哪里濾波需要指定通帶和阻帶而EMD的分解尺度完全由信號(hào)自身的極值點(diǎn)時(shí)間尺度決定不同信號(hào)分解出來(lái)的IMF帶寬和中心頻率都不一樣真正做到了“自適應(yīng)”。對(duì)軸承故障這種故障特征頻率會(huì)隨工況輕微變化的信號(hào)自適應(yīng)分解有明顯優(yōu)勢(shì)。2.2 新版MATLAB的emd函數(shù)一行調(diào)用如果你使用的是比較新的MATLAB環(huán)境R2018a及之后版本信號(hào)處理工具箱里已經(jīng)內(nèi)置了emd函數(shù)調(diào)用起來(lái)非常簡(jiǎn)單% sig一列振動(dòng)信號(hào)fs為采樣率 [imf, residual, info] emd(sig);其中輸出imf是一個(gè)矩陣每一行對(duì)應(yīng)一個(gè)IMF分量按頻率從高到低排列residual是最后的殘余趨勢(shì)項(xiàng)info結(jié)構(gòu)體里保存了篩選迭代次數(shù)、包絡(luò)均值等過(guò)程信息。需要限制分解層數(shù)時(shí)可以加一個(gè)參數(shù)[imf, residual] emd(sig, MaxNumIMF, 5);這個(gè)選項(xiàng)很實(shí)用。軸承振動(dòng)信號(hào)往往很長(zhǎng)如果放任EMD一直分下去可能會(huì)分出十幾個(gè)甚至幾十個(gè)IMF大多是低頻趨勢(shì)和緩慢變化的基線跟故障沖擊關(guān)系不大還白白浪費(fèi)時(shí)間。限制為前5個(gè)或前6個(gè)基本能把對(duì)故障敏感的高頻振蕩分量都保留下來(lái)。如果你手頭的MATLAB版本比較老沒(méi)有內(nèi)置函數(shù)那通常只能自己寫(xiě)簡(jiǎn)化版。核心篩選循環(huán)是這樣function imf simple_sift(x, max_iter) % x: 輸入信號(hào)行向量 % 返回一個(gè)滿足IMF條件的本征模態(tài)函數(shù) h x; for k 1:max_iter env_upper my_envelope(h, upper); % 用三次樣條求上包絡(luò) env_lower my_envelope(h, lower); % 求下包絡(luò) mean_env (env_upper env_lower) / 2; h_new h - mean_env; % 停止準(zhǔn)則相鄰兩次迭代的變化足夠小 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找極值點(diǎn)再用spline三次樣條插值構(gòu)造包絡(luò)。這個(gè)簡(jiǎn)化版能跑通但篩選停止準(zhǔn)則、端點(diǎn)處理都還有提升空間。所以如果不是單純?yōu)榱藢W(xué)習(xí)原理我不太建議在正式項(xiàng)目里從零手寫(xiě)EMD能用內(nèi)置函數(shù)就用內(nèi)置函數(shù)。2.3 分解完了一堆IMF到底選哪幾個(gè)EMD把人“選擇困難”的問(wèn)題放大了——一個(gè)樣本可能分解出好幾個(gè)IMF到底取哪些作為特征來(lái)源我的經(jīng)驗(yàn)是看相關(guān)系數(shù)也就是每個(gè)IMF和原始信號(hào)的相似程度。故障沖擊出現(xiàn)時(shí)它所在的IMF必然和原始信號(hào)有較高的相關(guān)性。可以用這樣的代碼篩選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); % 取相關(guān)性最高的前3個(gè)IMF除了相關(guān)系數(shù)也可以看每個(gè)IMF的能量占比能量占比高的分量往往攜帶更多的故障信息。實(shí)際操作中這兩種方法選出來(lái)的IMF基本一致我的習(xí)慣是優(yōu)先用相關(guān)系數(shù)法因?yàn)樗鼘?duì)幅值尺度不敏感更穩(wěn)定。需要特別提醒端點(diǎn)效應(yīng)是EMD繞不開(kāi)的坑。信號(hào)兩端的極值點(diǎn)不完整會(huì)讓包絡(luò)在端點(diǎn)處失真進(jìn)而污染IMF。實(shí)測(cè)中如果發(fā)現(xiàn)某個(gè)IMF的兩端明顯“翹起來(lái)”先做鏡像延拓或兩端衰減處理不要直接把整個(gè)分量拿去做熵值計(jì)算。這個(gè)坑我在第6章還會(huì)細(xì)說(shuō)。3. 樣本熵的原理與MATLAB函數(shù)實(shí)現(xiàn)3.1 樣本熵到底在衡量什么樣本熵Sample Entropy是Richman等人提出的一種復(fù)雜度度量用來(lái)衡量時(shí)間序列中產(chǎn)生新模式的可能性大小。它的思想可以這樣理解如果一段信號(hào)里長(zhǎng)度為 ( m ) 的模式在下一時(shí)刻還能被類(lèi)似的模式“復(fù)現(xiàn)”出來(lái)那么這段信號(hào)就是有規(guī)律的、可預(yù)測(cè)的反之如果復(fù)現(xiàn)比例很低說(shuō)明信號(hào)不規(guī)則、混亂程度高熵值就大。具體計(jì)算時(shí)對(duì)長(zhǎng)度為 ( N ) 的信號(hào) ( x )取一個(gè)窗口長(zhǎng)度 ( m )構(gòu)造 ( N-m ) 個(gè) ( m ) 維向量。然后設(shè)定一個(gè)容差 ( r )對(duì)任意兩個(gè)向量計(jì)算切比雪夫距離各維度差值絕對(duì)值的最大值統(tǒng)計(jì)有多少對(duì)向量的距離小于 ( r )記為 ( B )。接著把窗口長(zhǎng)度增加到 ( m1 )同樣統(tǒng)計(jì)匹配對(duì)數(shù)記為 ( A )。樣本熵定義為[ SampEn -\ln \frac{A}{B} ]這個(gè)比值 ( A/B ) 可以理解為“在原有一維模式上再延伸一維后模式仍然能匹配”的可能性。比值越大說(shuō)明信號(hào)越有規(guī)律熵值越小比值越小說(shuō)明信號(hào)越隨機(jī)、越不規(guī)律熵值越大。軸承正常運(yùn)行時(shí)振動(dòng)信號(hào)存在較強(qiáng)的旋轉(zhuǎn)周期性這個(gè)匹配比例比較高出現(xiàn)故障后沖擊和調(diào)制破壞了原有的模式穩(wěn)度匹配比例下降樣本熵隨之上升。3.2 參數(shù)怎么選m和r的經(jīng)驗(yàn)值樣本熵有兩個(gè)人為設(shè)定的參數(shù)嵌入維數(shù) ( m ) 和容差 ( r )。( m )通常取2或3。取2時(shí)計(jì)算穩(wěn)定對(duì)數(shù)據(jù)量要求不苛刻取3時(shí)能捕捉更長(zhǎng)時(shí)間的依賴(lài)但需要更長(zhǎng)數(shù)據(jù)否則匹配數(shù)容易過(guò)少。我默認(rèn)用2。( r )一般取原始信號(hào)標(biāo)準(zhǔn)差的0.1~0.25倍。r太小匹配數(shù)太少熵值波動(dòng)大r太大所有向量都匹配上了區(qū)分度消失。經(jīng)驗(yàn)上0.2倍標(biāo)準(zhǔn)差是個(gè)很平衡的值我在絕大多數(shù)軸承振動(dòng)數(shù)據(jù)上都用它。另外需要強(qiáng)調(diào)樣本熵對(duì)數(shù)據(jù)長(zhǎng)度敏感。太短的數(shù)據(jù)會(huì)讓匹配數(shù)統(tǒng)計(jì)不穩(wěn)定我建議每段信號(hào)至少取500點(diǎn)以上。軸承故障分析常用每段2048點(diǎn)或4096點(diǎn)這個(gè)長(zhǎng)度下熵值已經(jīng)比較穩(wěn)且計(jì)算時(shí)間還可以接受。3.3 一個(gè)可直接用的MATLAB樣本熵函數(shù)下面這個(gè)函數(shù)是我在實(shí)際項(xiàng)目里一直在用的版本用pdist2計(jì)算切比雪夫距離代碼簡(jiǎn)潔邏輯也清楚function se sample_entropy(data, m, r) % data: 輸入列向量或行向量 % m: 嵌入維數(shù)通常取2 % r: 容差通常取0.2*std(data) data data(:); N length(data); if N m 1 se NaN; return; end % 構(gòu)造m維模板向量 tm zeros(N - m, m); for i 1:m tm(:, i) data(i : N - m i - 1); end Dm pdist2(tm, tm, chebychev); % 統(tǒng)計(jì)距離小于r的匹配對(duì)數(shù)排除自身匹配 B (sum(Dm r, all) - (N - m)) / ((N - m) * (N - m - 1)); % 構(gòu)造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這個(gè)函數(shù)里“排除自身匹配”是關(guān)鍵也是樣本熵和近似熵最主要的區(qū)別之一。近似熵把向量和自身的匹配也計(jì)入了導(dǎo)致結(jié)果有偏樣本熵嚴(yán)格去掉自身匹配熵值的統(tǒng)計(jì)一致性和抗噪能力都更好。實(shí)測(cè)下來(lái)同樣的數(shù)據(jù)近似熵在輕微噪聲下容易跳變樣本熵要穩(wěn)得多。計(jì)算復(fù)雜度上這個(gè)函數(shù)是 ( O(N^2) ) 的因?yàn)橐?jì)算所有模板向量?jī)蓛芍g的距離。數(shù)據(jù)長(zhǎng)度超過(guò)10000點(diǎn)時(shí)耗時(shí)就會(huì)比較明顯。軸承診斷一般用2048點(diǎn)或4096點(diǎn)跑起來(lái)沒(méi)問(wèn)題。4. 把EMD和樣本熵串起來(lái)——特征提取實(shí)戰(zhàn)流程4.1 為什么不能直接對(duì)原始信號(hào)算樣本熵你可能會(huì)想既然樣本熵能衡量信號(hào)復(fù)雜度直接對(duì)原始振動(dòng)信號(hào)算一個(gè)熵值不就行了為什么非要EMD分解一下直接算當(dāng)然可以但原始信號(hào)里包含轉(zhuǎn)速基頻、隨機(jī)噪聲、結(jié)構(gòu)共振等太多種成分熵值是一個(gè)“混合平均值”故障特征被稀釋了。把信號(hào)先用EMD分解開(kāi)再針對(duì)性地選幾個(gè)和故障相關(guān)的IMF分量算熵就等于把有用的信號(hào)“單獨(dú)拎出來(lái)”做體檢敏感度會(huì)高很多。我可以給你一個(gè)很直觀的類(lèi)比原始信號(hào)像一鍋亂燉的湯你想判斷里面的鹽放多了沒(méi)有直接喝一口整鍋湯味道往往被其他食材干擾更好的辦法是把湯分層、把渾濁物過(guò)濾掉再倒出一小勺來(lái)嘗。EMD做的就是過(guò)濾分層的工作樣本熵就是嘗味道。4.2 完整特征提取步驟我在做滾動(dòng)軸承診斷時(shí)特征提取已經(jīng)固定成了一套流程每一段樣本都按這個(gè)順序走數(shù)據(jù)分幀把連續(xù)采集的振動(dòng)信號(hào)切成固定長(zhǎng)度的小段每段2048點(diǎn)或4096點(diǎn)。樣本熵對(duì)長(zhǎng)度敏感所有樣本必須用相同的長(zhǎng)度否則沒(méi)有可比性。預(yù)處理每段信號(hào)減去均值去除直流分量。如果有明顯的趨勢(shì)項(xiàng)再用detrend去除。EMD分解調(diào)用內(nèi)置emd設(shè)置MaxNumIMF為5或6。選擇IRF分量計(jì)算每個(gè)IMF與原始信號(hào)的相關(guān)系數(shù)選出相關(guān)性最高的前3個(gè)IMF。這一步保證了特征向量的維度固定因?yàn)椴煌瑯颖痉纸獬龅腎MF數(shù)量可能不一樣。計(jì)算樣本熵對(duì)選中的每一個(gè)IMF先求其標(biāo)準(zhǔn)差再以 ( r0.2 \times std ) 為容差調(diào)用sample_entropy。把每個(gè)IMF的熵值拼接成該樣本的特征向量比如[se_imf1, se_imf2, se_imf3]。特征標(biāo)準(zhǔn)化收集所有樣本的特征組成矩陣后做z-score標(biāo)準(zhǔn)化避免某個(gè)IMF的熵值量綱差異過(guò)大。這個(gè)對(duì)后續(xù)分類(lèi)很重要。4.3 可復(fù)現(xiàn)的MATLAB主腳本示例下面這段主腳本展示了上面完整的流程。我假設(shè)你的數(shù)據(jù)已經(jīng)按列存在一個(gè)矩陣data中每一列是一段樣本并且已知每段樣本的類(lèi)別標(biāo)簽存在labels里。% 參數(shù)設(shè)置 m 2; r_ratio 0.2; num_imf_keep 3; max_imf 5; num_samples size(data, 2); % 預(yù)分配特征矩陣 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個(gè)請(qǐng)檢查數(shù)據(jù)。, k, num_imf_keep); continue; end % 按相關(guān)系數(shù)選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); % 對(duì)每個(gè)選定IMF計(jì)算樣本熵 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處理見(jiàn)后文踩坑部分 end features(k, j) se; end end % 標(biāo)準(zhǔn)化 features zscore(features);注意我把最外層的循環(huán)寫(xiě)成了parfor這樣在多核機(jī)器上能明顯加快速度。如果你的電腦沒(méi)有并行計(jì)算工具箱改成普通for也可以正常跑。這個(gè)腳本做了一件很重要的事把特征維度固定下來(lái)。固定維度是后續(xù)機(jī)器學(xué)習(xí)建模的前提也是很多入門(mén)教程不會(huì)特別強(qiáng)調(diào)的地方。你可能會(huì)問(wèn)如果某段信號(hào)的IMF數(shù)量不足3個(gè)怎么辦實(shí)測(cè)中一段幾百點(diǎn)的純正弦信號(hào)可能只能分解出1個(gè)IMF但這在真實(shí)軸承振動(dòng)數(shù)據(jù)里很少發(fā)生。為了穩(wěn)健我一般會(huì)在腳本里加一個(gè)判斷IMF不足時(shí)把缺失的熵值置為一個(gè)特殊標(biāo)記要么用該特征列的均值填充要么直接把該樣本剔除。具體用哪種策略取決于你要不要保留所有樣本。5. 三狀態(tài)數(shù)據(jù)實(shí)測(cè)特征值分布與初步分類(lèi)結(jié)果5.1 我這組示例數(shù)據(jù)的情況為了把流程完整演示下來(lái)我用了手頭一組在實(shí)驗(yàn)臺(tái)上實(shí)測(cè)的滾動(dòng)軸承振動(dòng)數(shù)據(jù)。數(shù)據(jù)來(lái)自某型深溝球軸承通過(guò)人為加工故障點(diǎn)模擬三種典型狀態(tài)正常狀態(tài)、外圈故障狀態(tài)和內(nèi)圈故障狀態(tài)。工況條件為恒定轉(zhuǎn)速約1000轉(zhuǎn)/分采樣率12kHz每個(gè)狀態(tài)各取60段樣本每段2048個(gè)點(diǎn)。為了方便我已經(jīng)把數(shù)據(jù)組織成一個(gè)data矩陣列數(shù)為180每列都是一段去均值后的信號(hào)。這里要說(shuō)明這組數(shù)據(jù)只用于演示整個(gè)診斷流程不同實(shí)驗(yàn)室、不同傳感器的數(shù)據(jù)熵值的絕對(duì)量級(jí)會(huì)有差異但“正常狀態(tài)熵值低、故障狀態(tài)熵值高”的趨勢(shì)是一致的。5.2 跑完后的樣本熵統(tǒng)計(jì)用上面第4章的腳本跑完后我統(tǒng)計(jì)了每個(gè)狀態(tài)下前3個(gè)IMF的樣本熵均值和標(biāo)準(zhǔn)差結(jié)果大致如下?tīng)顟B(tài)IMF1樣本熵IMF2樣本熵IMF3樣本熵正常狀態(tài)0.28 ± 0.030.22 ± 0.040.19 ± 0.03外圈故障0.61 ± 0.080.48 ± 0.070.39 ± 0.06內(nèi)圈故障0.89 ± 0.100.76 ± 0.090.62 ± 0.08每一格前面的數(shù)是均值后面的數(shù)是標(biāo)準(zhǔn)差?？梢钥吹饺?lèi)狀態(tài)之間的差異非常明顯。正常狀態(tài)的所有IMF樣本熵值都偏低說(shuō)明正常軸承的振動(dòng)信號(hào)規(guī)律性好試看圖基本是平穩(wěn)的旋轉(zhuǎn)振動(dòng)疊加少量隨機(jī)噪聲。內(nèi)圈故障的樣本熵值最高因?yàn)閮?nèi)圈故障時(shí)沖擊信號(hào)要經(jīng)過(guò)滾動(dòng)體、保持架等更多傳遞路徑信號(hào)里調(diào)制成分更豐富、更不規(guī)則。外圈故障居中但也遠(yuǎn)高于正常狀態(tài)。如果把這180個(gè)樣本的特征用散點(diǎn)圖畫(huà)出來(lái)你會(huì)發(fā)現(xiàn)三個(gè)簇分得很開(kāi)。我通常用前三維特征畫(huà)三維散點(diǎn)圖代碼就幾行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(正常, 外圈故障, 內(nèi)圈故障); xlabel(IMF1樣本熵); ylabel(IMF2樣本熵); zlabel(IMF3樣本熵);5.3 順手驗(yàn)證一下分類(lèi)準(zhǔn)確率雖然標(biāo)題重點(diǎn)是特征提取但特征提得好不好最終還是要靠分類(lèi)結(jié)果來(lái)說(shuō)話。我用MATLAB自帶的fitcecoc訓(xùn)練一個(gè)簡(jiǎn)單的SVM多分類(lèi)模型隨機(jī)留出30%樣本做測(cè)試X features; % 已經(jīng)zscore標(biāo)準(zhǔn)化 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(測(cè)試集準(zhǔn)確率: %.2f%%\n, acc * 100);在這組演示數(shù)據(jù)上測(cè)試集準(zhǔn)確率可以達(dá)到95%以上。需要反復(fù)強(qiáng)調(diào)的是這個(gè)結(jié)果是在實(shí)驗(yàn)室固定工況數(shù)據(jù)上取得的真實(shí)工業(yè)現(xiàn)場(chǎng)會(huì)有載荷波動(dòng)、轉(zhuǎn)速波動(dòng)、軸承座共振以及齒輪箱等其他振動(dòng)源的干擾準(zhǔn)確率會(huì)下降。但這并不妨礙EMD樣本熵作為特征提取方案的實(shí)用性——它把最關(guān)鍵的復(fù)雜度差異抓出來(lái)了剩下的偏差更多屬于數(shù)據(jù)預(yù)處理、工況適應(yīng)性和分類(lèi)器調(diào)優(yōu)的范疇。5.4 結(jié)果解讀為什么故障后樣本熵會(huì)變大解釋一下為什么會(huì)有這樣的結(jié)果。正常軸承運(yùn)行時(shí)振動(dòng)能量主要集中在旋轉(zhuǎn)頻率和少量倍頻上時(shí)域波形相對(duì)規(guī)則EMD分解后各個(gè)IMF基本都是窄帶的準(zhǔn)周期振蕩樣本熵自然不高。故障發(fā)生后損傷點(diǎn)與滾動(dòng)體接觸產(chǎn)生周期沖擊沖擊激起結(jié)構(gòu)共振形成調(diào)幅調(diào)頻IMF里不再只有單一頻率的振蕩而是包含了沖擊的間歇性和調(diào)制帶來(lái)的復(fù)包絡(luò)。這種“規(guī)律之外又有變化”的狀態(tài)正是樣本熵最敏感的地方。沖擊越嚴(yán)重、調(diào)制越復(fù)雜IMF的樣本熵就越高。當(dāng)然熵值升高并不是無(wú)限的。如果故障發(fā)展到晚期沖擊變成連續(xù)的、無(wú)結(jié)構(gòu)的寬帶噪聲某些IMF又會(huì)變得接近隨機(jī)噪聲熵值可能不再單調(diào)上升。因此我在工程實(shí)踐中通常把樣本熵作為“早期預(yù)警”和“狀態(tài)分類(lèi)”特征而不是試圖用它精確測(cè)量故障尺寸的大小。6. 我在迭代中踩過(guò)的坑和現(xiàn)在的標(biāo)準(zhǔn)做法6.1 不同版本MATLAB的emd函數(shù)差異最早我是在一個(gè)沒(méi)有內(nèi)置emd的老環(huán)境里做開(kāi)發(fā)代碼里到處是第三方工具箱的調(diào)用。換到新環(huán)境后發(fā)現(xiàn)函數(shù)接口變了原來(lái)的腳本報(bào)錯(cuò)一串又一串。后來(lái)我定了兩條規(guī)矩第一新寫(xiě)的代碼一律只用內(nèi)置emd函數(shù)不依賴(lài)任何第三方工具箱保證拿到別的機(jī)器上能直接跑第二如果項(xiàng)目環(huán)境版本太老先給信號(hào)做鏡像延拓再跑我自己那個(gè)簡(jiǎn)化版篩選但只用于驗(yàn)證算法不用于最終生產(chǎn)。經(jīng)驗(yàn)是能用官方函數(shù)就盡量用官方函數(shù)省下的時(shí)間足夠多。6.2 每個(gè)樣本分解出的IMF個(gè)數(shù)不一致特征矩陣拼不齊這是非常容易踩的坑。EMD的分解結(jié)果是由信號(hào)自身決定的有的樣本能分解出6個(gè)IMF有的只能分解出3個(gè)。如果你直接把所有IMF的樣本熵都堆成一個(gè)特征向量不同樣本的特征長(zhǎng)度不一樣分類(lèi)模型根本沒(méi)法用。我采取的策略就是第4章寫(xiě)的固定保留前3個(gè)相關(guān)系數(shù)最高的IMF。如果某些樣本分解出的IMF數(shù)量不夠先檢查數(shù)據(jù)分段是否太短然后決定是剔除樣本還是用該特征列的均值填充。寧可損失少量樣本也不要破壞特征維度的統(tǒng)一性。6.3 樣本熵的r參數(shù)不同區(qū)分度完全不一樣我在另一組數(shù)據(jù)上試過(guò)把 ( r ) 從0.1倍標(biāo)準(zhǔn)差調(diào)到0.3倍標(biāo)準(zhǔn)差結(jié)果相同的故障樣本熵值變化很大分類(lèi)準(zhǔn)確率也跟著波動(dòng)。這說(shuō)明樣本熵不是“一次性定下來(lái)就萬(wàn)事大吉”的參數(shù)。我現(xiàn)在養(yǎng)成的習(xí)慣是對(duì)每批新數(shù)據(jù)先用少量樣本做一次參數(shù)掃描比如 ( r \in {0.1, 0.15, 0.2, 0.25, 0.3} \times std )跑完分類(lèi)看看準(zhǔn)確率選一個(gè)最穩(wěn)的固定下來(lái)。大多數(shù)情況下0.2倍標(biāo)準(zhǔn)差是性?xún)r(jià)比最高的但不要排斥其他候選值。6.4 計(jì)算慢到懷疑人生因?yàn)闃颖眷匾銉蓛上蛄恐g的距離矩陣當(dāng)每段信號(hào)4096點(diǎn)、樣本數(shù)量幾百個(gè)時(shí)整個(gè)腳本可能要跑幾個(gè)小時(shí)。解決思路有三個(gè)一是用parfor并行把每個(gè)樣本的分解和熵值計(jì)算攤到多個(gè)核心上二是只計(jì)算被選中的那3個(gè)IMF不要對(duì)所有IMF都算熵三是在保證熵值穩(wěn)定的前提下把每段點(diǎn)數(shù)從4096降到2048。實(shí)際算下來(lái)2048點(diǎn)已經(jīng)足夠穩(wěn)定沒(méi)必要追求越長(zhǎng)越好。6.5 熵值出現(xiàn)NaN怎么辦我第一次在純正弦仿真信號(hào)上測(cè)試樣本熵函數(shù)時(shí)返回結(jié)果居然全是NaN。原因很簡(jiǎn)單純正弦信號(hào)的 ( m1 ) 維模板向量匹配對(duì)數(shù) ( A ) 可能為0因?yàn)樾盘?hào)太規(guī)則了稍微延長(zhǎng)一維之前匹配的模式就不再匹配了。實(shí)際軸承信號(hào)一般不會(huì)這么極端但當(dāng)某一段信號(hào)周期性過(guò)強(qiáng)時(shí)也可能觸發(fā)。我的處理方法是檢測(cè)isnan如果是就賦為0或者用該特征列的中位數(shù)代替。千萬(wàn)不能放著NaN不管否則后面z-score會(huì)全部變成NaN整個(gè)特征矩陣就廢了。6.6 我現(xiàn)在面對(duì)一批新振動(dòng)數(shù)據(jù)的標(biāo)準(zhǔn)做法踩過(guò)這些坑之后我的工作流已經(jīng)相當(dāng)固定了。拿到一批新的軸承振動(dòng)數(shù)據(jù)我第一件事不是急著畫(huà)頻譜圖而是先做分幀、去均值然后跑EMD限制到前5個(gè)IMF用相關(guān)系數(shù)法選出前3個(gè)計(jì)算樣本熵再把特征矩陣z-score標(biāo)準(zhǔn)化最后用分類(lèi)器快速驗(yàn)證一下可分性。如果熵值分布能把不同狀態(tài)拉開(kāi)再往里加時(shí)域特征比如峭度、峰值因子也不遲如果拉不開(kāi)那就需要回頭看數(shù)據(jù)采集質(zhì)量、轉(zhuǎn)速波動(dòng)和傳感器安裝位置。這個(gè)流程的好處是整個(gè)鏈條的每一步都可以單獨(dú)排查。我在項(xiàng)目里會(huì)把每個(gè)中間結(jié)果都保存下來(lái)分解后的IMF波形、每個(gè)IMF的熵值、特征矩陣都存成.mat文件。哪天發(fā)現(xiàn)某個(gè)樣本分類(lèi)錯(cuò)了直接翻中間變量很快就能定位是分解的問(wèn)題、熵值計(jì)算的問(wèn)題還是分類(lèi)器參數(shù)的問(wèn)題。這種“可復(fù)現(xiàn)、可追溯”的做法在做故障診斷課題時(shí)尤其重要。從個(gè)人的體會(huì)來(lái)說(shuō)EMD加樣本熵這套組合在滾動(dòng)軸承故障診斷里并不是什么花哨的新技術(shù)但勝在穩(wěn)定、實(shí)用、解釋性強(qiáng)。很多更復(fù)雜的深度學(xué)習(xí)方法需要的數(shù)據(jù)量和調(diào)參成本都大得多。如果你的目標(biāo)是想快速搭建一套能用的狀態(tài)識(shí)別基線我特別建議先把這條路走通一次。走過(guò)這一遍你對(duì)“什么是特征”“為什么要自適應(yīng)分解”“熵值變化到底反映了什么”都會(huì)有一個(gè)很清晰的感知之后再去看其他高級(jí)方法也會(huì)更容易融會(huì)貫通。