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