柴油機振動信號快速稀疏分解與二維特征編碼*

王 旭， 岳應娟， 蔡艷平

(火箭軍工程大學理學院 西安，710025)

引 言

柴油機運行中既有旋轉運動，又有往復運動，其振動響應信號復雜，耦合嚴重，具有較強的非線性、非平穩時變特征[1]。采用時頻分析是通過對柴油機振動時頻圖像的識別來實現對故障的診斷，其技術關鍵在于準確的時頻表征和有效的特征提取，眾多學者對此進行了大量的探索。文獻[2]提出利用灰度共生矩陣提取柴油機振動信號三階累積量圖像的故障診斷方法。文獻[3]將Wigner分布與分形維數相結合，對柴油機進行故障診斷。文獻[4]將圖像分割理論引入柴油機故障診斷之中。文獻[5]提出基于經驗模態分解和維格納分布的柴油機振動信號時頻表征方法。以上方法將柴油機的故障診斷問題轉化為柴油機振動時頻圖像的模式識別問題，取得了較好的效果。

時頻表征方面，利用小波分析(wavelet analysis，簡稱WT)、維格納分布(Wigner-Ville distribution ，簡稱WVD)及短時傅里葉變換(short time Fourier transformation，簡稱STFT)等方法對柴油機缸蓋振動信號進行時頻表征是應用較廣泛的信號處理手段[6]。匹配追蹤算法(matching pursuit，簡稱MP)是一種基于過完備冗余時頻字典對信號進行稀疏分解的方法，相比于以上時頻分析方法，MP算法自適應性更好，對時頻分布中各分量的刻畫能力更強。但是MP算法計算量和存儲量相當大，限制了其在柴油機振動信號處理中的應用[7]。

特征提取方面，為了不依靠先驗知識實現對不同工況時頻圖像的自動分類，需要對時頻圖像進行特征提取，其本質是對時頻圖像數據降維的同時最大限度保留不同樣本間的差異化信息。非負矩陣分解算法(non-negative matrix factorization, 簡稱NMF)是一種基于局部特征的矩陣分解方法，由于添加了非負的限制條件，能夠保證分解結果的可解釋性。文獻[8-9]將NMF應用于柴油機和軸承時頻圖像矩陣的特征提取，取得了較好的效果。但是NMF算法在降維前需將圖像矩陣向量化處理，破壞了圖像矩陣間的空間位置信息，并且向量化處理后的圖像矩陣往往維數較高，使得NMF的計算效率較低。對此，文獻[10-12]從不同角度提出了二維非負矩陣分解(2-dimensional non-negative matrix factorization, 簡稱2DNMF)方法，分解前不需將圖像矩陣向量化，提高了特征提取效果，被成功應用于故障診斷領域[13-14]。但是2DNMF在分解前都是將所有訓練圖像矩陣按行或列拼合，拼合后的初始分解矩陣維度依然較大，并且2DNMF算法沒有考慮到不同類圖像間的差異化信息，將所有圖像樣本一致對待，統一求解投影矩陣，這顯然是不利于模式識別的。

基于以上分析，筆者著眼于柴油機故障診斷中的時頻表征與特征提取過程，針對MP算法分解效率問題和2DNMF特征提取性能問題，提出了基于AMP算法和雙向二維非負矩陣分解的柴油機振動信號時頻分析方法。將該方法應用于氣門間隙正常、氣門間隙過小、氣門間隙過大和氣門漏氣4種不同氣門狀態信號的診斷試驗中，故障識別正確率最高可達100%，充分證明該方法用于柴油機自動故障診斷的有效性。

1 信號的稀疏分解及時頻表征

1.1 改進的自適應MP算法實現

匹配追蹤算法在計算和存儲上的瓶頸主要在于過完備字典的制備，在分解時需要遍歷一個龐大的原子字典中的所有原子后再找到最匹配的原子組合。柴油機振動信號通常會聚集于有限寬度的頻帶中，而Gabor原子具有最好的時頻聚集性，所以利用Gabor原子字典分解信號時對其他頻率范圍內原子的搜索過程其實是無意義的。如果令字典隨殘余信號的功率譜分布自適應更新，縮小原子搜索范圍，在保證算法稀疏性的基礎上可以大幅提高算法的計算效率。基于此，筆者提出了自適應匹配追蹤算法。

設H為Hilbert空間，D={gγ(t)}γ∈Γ為H中過完備字典。原子gγ(t)由參數γ描述。Gabor原子[15]可表示為

(1)

其中：g(t)=e-πt2為高斯窗函數；原子參數γ=(s,u,ξ)；s為尺度因子；u為位移因子；ξ為頻率因子。

原子經歸一化處理后‖gγ(t)‖=1。設f為待分解信號，且f∈H，有

(2)

(3)

(4)

(5)

Gabor原子在頻域的能量主要集中在以調制頻率為中心的頻率區域內，自適應字典Dn中所有的原子主要能量均聚集于中心頻率ξ0附近，使得該字典中的原子和殘留信號能有較好的匹配性。新字典為Dn的參數集合Γn∈R+×R中包含γ=aj,pajΔu,kajΔξ，a=2，Δu=1/2，Δξ=π，0

AMP算法流程如下。

1) 定義s=2j，其中：j∈(0,log2N)；u=s。計算原子gγ的包絡g((t-u)/s)，并保存為包絡庫。由于Gabor原子的能量主要集中在區域0,2s內，為了減少計算量，每個包絡可以只計算區域0,2s內的數值點。

4) 找到上步運算中的互相關系數最大值，確定相應的匹配到的原子參數。

1.2 仿真分析

為分析AMP算法的性能，建立一個多分量仿真信號x(t)，信號長度為256，由5個具有高斯包絡的原子信號分量疊加而成，采樣頻率歸一化為1 Hz，5個信號分量的時、頻域分布中心分別位于(t1,f1)=(60,0.1)，(t2,f2)=(60,0.4)，(t3,f3)=(130,0.25)，(t4,f4)=(200,0.1)和(t5,f5)=(200,0.4)，仿真信號表達式為式(6)。信號的時域波形與時頻域分布對應如圖1所示。

(6)

圖1 仿真信號時域波形及時頻域分布對應圖Fig.1 The time and frequency domain waveform

表16次迭代殘余信號能量及耗時

Tab.1Residualsignalenergyandtimecostin6iterations

迭代次數MP耗時/msAMP耗時/msMP殘余能量/%AMP殘余能量/%1504.0434.64982.6479.072394.4732.43664.6258.173382.0442.13244.0437.274397.4592.00223.3119.175418.7642.5012.951.076481.9362.1092.150.82

表1為6次迭代過程中MP與AMP兩種算法對應的原子匹配時間和殘余信號能量百分比。理論上，仿真信號由5個原子分量復合而成，所以無論MP算法還是AMP算法都應在5次迭代后匹配出所有原子分量。MP與AMP算法在第5次迭代后，殘余信號的能量分別達到了2.15%和0.82%，均較好地匹配出了信號的原子分量，實現了信號的稀疏分解。由于本研究中AMP算法根據殘余信號功率譜分布將搜索字典限定在固定頻帶中，其搜索到的原子對于分解信號的匹配性更好，所以每次迭代中AMP算法的殘余能量都低于MP算法。同時，由于AMP算法中自適應匹配的原子字典較MP算法中的過完備冗余字典維度有較大幅的壓縮，使得AMP算法較MP算法的計算效率提高了150～200倍，從而能夠利用AMP算法對高維度、強耦合的信號進行分析。

為了獲得信號的時頻分布，利用AMP算法將仿真信號分解為5個原子分量，通過二次疊加每個原子的Wigner-Ville分布得到原信號的時頻分布。仿真信號的WVD分布與AMP稀疏分解時頻分布情況如圖2所示。

圖2 仿真信號時頻分布Fig.2 Time-frequency distribution of the simulation signal

由圖2看出，直接對仿真信號進行WVD分析，每兩個原子之間都會產生一個交干擾叉項，嚴重影響對原信號分量的分析。利用AMP算法得到的信號稀疏分解時頻分布中交叉項的問題得到了很好的解決，所有原子分量能量分布均衡，且保持了優良的時頻聚集性，所以利用AMP算法對柴油機振動信號進行稀疏分解，能夠使各時頻分量物理意義更加明確，有利于對柴油機故障的判別。

2 基于雙向2DNMF的時頻特征編碼方法

為了獲得更好的特征提取效果，筆者提出一種雙向二維非負矩陣分解算法，其與NMF，2DNMF有著相同的目標函數及迭代規則，但是架構不同。TD2DNMF算法分別對行基、列基矩陣進行求解，進一步得到TD2DNMF的二維基，用于對時頻圖像進行特征編碼。

首先，對行基進行求解，求解時考慮到不同類別樣本的差異性，對各類別樣本矩陣并行運算得到相應的基矩陣和系數矩陣；然后，組合成整個訓練樣本的行基投影矩陣和系數矩陣，由此得出的系數矩陣具有更好的稀疏度。假設有k類模式ω1,ω2,…,ωk，每類模式有m個訓練樣本圖像Aa,b,a=1,2,…,k;b=1,2,…,m，每個圖像對應大小為p×q的矩陣。將所有圖像用矩陣X表示為

Xp×qmk=[X1,X2,…,Xk]

(7)

其中：Xa=[Aa,1,Aa,2,…，Aa,m]；a=1,2,…,k。

對于非負矩陣Xa，求解非負矩陣La和Ha，滿足

(8)

其中：Xa∈Rp×qm,La∈Rp×r,Ha∈Rr×qm；r為特征維數，滿足(p+qm)r

為了描述X≈L·H的近似效果，利用矩陣X與L·H間的K-L 散度作為近似誤差，對應的目標函數[16]為

(9)

其對應的優化問題為

minD(V‖WH)s.t.W,H≥0

(10)

相應的迭代規則為

Hrj←Hrj∑iLirXij/(LH)ij

(11)

Lir←Lir∑jHrjXij/(LH)ij

(12)

(13)

由此得到每一類別圖像數據矩陣對應的分解矩陣因子，對于矩陣Xp×qmk對應的TD2DNMF分解因子可由每一類別的分解因子組合而成

(17)

(18)

與行基的求解過程類似，根據K-L散度構造目標函數進行迭代運算，確定列基分解的特征維數d，得到矩陣HT對應的TD2DNMF分解因子

(19)

任意圖像矩陣Ap×q在行基和列基上的投影系數矩陣為

(20)

TD2DNMF的二維基定義為

B=LRT

(21)

通過將柴油機振動時頻圖像向TD2DNMF的二維基進行投影，可得到相對應的特征編碼。利用分類器對特征編碼進行識別，可實現時頻圖像的自動分類，進一步可實現對故障的判別。

3 柴油機故障診斷實例分析

基于AMP振動信號稀疏分解與TD2DNMF二維時頻特征編碼的柴油機故障診斷方法流程如圖3所示。

圖3 故障診斷流程框圖Fig.3 The flowchart of fault diagnosis

3.1 試驗設備與工況設置

以6135G型柴油機為研究對象，用加速度傳感器測量缸蓋振動信號，試驗平臺如圖4所示。分別將氣門間隙設置成不同狀態，用以模擬氣門的不同間隙異常故障。在轉速穩定于1 500 r/min時測量第2缸的缸蓋振動信號，單通道采樣頻率為25 kHz，試驗具體工況設置如表2所示。試驗中運行狀態為空載，氣門正常間隙值為0.30 mm，使用0.06 mm模擬氣閥間隙過小，使用0.50 mm模擬氣閥間隙過大，在氣閥上開4 mm×1 mm孔來模擬嚴重漏氣故障。共采集柴油機氣門4種狀態下各60組振動信號，總計240個樣本。

表2 4種試驗工況設置

圖4 試驗平臺及傳感裝置設置Fig.4 Experimental Platform and sensing device

3.2 信號稀疏分解生成時頻圖像

圖5為4種工況下振動信號的時域波形。進氣門開啟的角度在310°附近，關閉的角度在-120°附近；排氣門開啟的角度在120°附近，關閉的角度在-310°附近，柴油機在0°產生燃燒激振。

圖5 信號時域波形Fig.5 Time domain waveform of signal

圖6 4類工況信號分解重構信號及殘余信號Fig.6 Signal construction and residues of four states

用AMP算法對振動響應信號進行分解和重構，迭代200次后殘余信號能量達到原信號能量的5%以內。4種工況分解的重構信號和殘余信號如圖6所示。可以看到,原信號中的所有沖擊分量均得到了很好的匹配。在時頻域上對振動信號進行分析，將AMP分解原子的Wigner-Ville分布進行疊加得到各工況的時頻分布圖像，如圖7所示。

對柴油機振動信號進行AMP時頻表征，4類工況信號分別耗時9.022，8.678，8.869和8.955s，平均耗時8.881s，而相同條件下對振動信號直接進行WVD時頻表征，4類工況信號平均耗時814.65s，說明AMP時頻分析方法更能夠滿足柴油機故障診斷的時效性。從圖7可以看出，正常工況對應的頻率能量較大，曲軸轉角0°處對應能量更高，說明混合氣體燃燒效率更高。狀態2，3，4對應的頻率能量小，說明氣體燃燒做功不充分。對于狀態2，4，由于氣門漏氣會導致燃燒不充分，狀態4的漏氣最嚴重，所以狀態4的燃燒功率最低。氣門間隙故障不但影響混合氣體的燃燒效率，而且影響氣門落座對缸蓋的沖擊。對比狀態2，3和4曲軸轉角-310°處頻率能量，狀態3的頻率能量明顯大于狀態2，可以知道狀態3的排氣門間隙大，而狀態2的排氣門間隙小。

圖7 4類工況的AMP時頻分布情況Fig.7 AMP time-frequency distribution of four states

通過以上分析，時頻圖像清楚地表明不同間隙的氣門落座沖擊以及混合氣體的燃燒效率所占的頻率分量是不同的。因此，柴油機振動信號的稀疏分解時頻圖像能夠較好地反映出故障特征。但是，以上分析需要一定的先驗知識做前提，要依靠計算機實現故障的自動診斷，需進一步提取時頻圖像中的特征參量。

3.3 TD2DNMF時頻特征編碼提取與故障識別

取采集到的240個信號作為研究對象并分別繪制時頻圖像，由于彩色圖像對應的是三維矩陣，為了便于計算機處理，對圖像進行預處理，采用閾值平均法將其轉化為灰度圖像，相應得到240個420×560像素點的矩陣樣本。

圖8為特征維數R=d=8時，柴油機4類工況振動信號時頻分布圖像測試集對應的特征編碼。圖中每個像素的幅值嚴格與樣本系數值對應，由于文章篇幅有限，每種工況下選取5個樣本顯示，每個樣本編碼圖的橫、縱坐標僅代表矩陣維度。圖中每一行代表一種柴油機工況，從上到下依次為氣門間隙正常、過小、過大和嚴重漏氣工況。可以看出，TD2DNMF對數據進行了有效降維，將420×560維數據壓縮到8×8維，大大降低了模式識別的計算復雜度。同種工況的編碼矩陣類內差異小，不同工況間編碼矩陣的類間差異較大。對圖像進行TD2DNMF特征提取，在對數據進行有效降維的同時，不同工況對應時頻圖像間的差異化信息得到了較大程度的保留，有利于提高故障模式識別的正確率。

圖8 4類工況時頻圖像特征編碼Fig.8 Feature code of time-frequency distribution of 4 states

為比較不同算法的計算效率，分別利用傳統非負矩陣分解算法、文獻[10]提出的二維非負矩陣分解算法(2-dimensional NMF-Zhang，簡稱2DNMF-Z)、文獻[11]提出的二維非負矩陣分解算法(2-dimensional NMF-Gu，簡稱2DNMF-G)、文獻[12] 提出的并行二維非負矩陣分解算法((2D)2NMF)以及筆者提出的TD2DNMF算法對4種工況所有樣本進行特征提取。NMF特征維數設定為二維方法編碼矩陣行、列維數的乘積。4種算法統一最大迭代步長為100，目標函數容忍誤差為10-5。表3為5種算法特征提取的計算時間，均不包含圖像載入的時間。試驗環境為 Matlab R2012b，AMD A8處理器，1.90 GHz主頻CPU，4GB內存，Win7操作系統。

由表3可見，隨著特征維數的增加，5種算法特征提取時間整體呈增長的趨勢。二維NMF的計算效率明顯要高于一維NMF，一維NMF對圖像矩陣進行向量化后數據矩陣維度過大。例如，筆者所用圖像維度為420×560，120幅訓練圖像組成的數據矩陣維度為120×235 200，計算任務量十分繁重。二維NMF算法中，2DNMF-G的計算效率略低于其他幾種，這是由于2DNMF-G在計算時將原始圖片分別按行、列拼合，初始分解矩陣維數分別為420×67 200和560×50 400，數據維度依然很大。(2D)2NMF算法在初始分解矩陣的組成上與2DNMF-G相同，但是采用并行運算的方法，計算效率有所提高。2DNMF-Z算法的初始分解矩陣僅將圖像按行拼接，再將一次NMF分解后得到的系數矩陣分塊轉置再拼接后作為新的初始分解矩陣進行求解，兩個初始分解矩陣的維度分別為420×67 200和560×960，維度小于2DNMF-G與(2D)2NMF，計算效率與后者相當。本研究方法TD2DNMF將圖像矩陣按行拼接后按照各自的類別進行分解，將得到的4個系數矩陣各自分塊轉置再拼接，組成新的初始分解矩陣，兩次分解的數據矩陣維度分別為420×67 200和560×120。與其他幾類算法相比，本研究方法初始分解矩陣維度更低，有效提高了算法的計算效率。

為比較5種方法特征提取效果，采用最近鄰分類器(nearest neighbor classifier, 簡稱NNC)對提取到的特征參數進行模式識別。診斷時分別采用柴油機振動信號的Wigner-Ville時頻分布圖和AMP稀疏分解時頻分布圖的240個樣本，每種工況60個樣本，從中隨機抽取30個作為訓練樣本，其余作為測試樣本。以NNC分類正確率評價5種方法的特征提取效果，為減小隨機誤差，試驗重復10次取平均準確率作為最終結果，如圖9所示。

表3不同方法特征提取效率對比

Tab.3Timecostindifferentfeatureextractionmethodss

特征維數特征提取方法NMF2DNMF-Z2DNMF-G(2D)2NMFTD2DNMF8(4×2)304.884.8105.980.0646.216(4×4)367.286.9106.381.4347.324(8×3)408.694.1133.787.2657.932(8×4)463.898.1135.293.559.840(8×5)622.298.4134.795.462.448(8×6)882.898.7136.996.363.956(8×7)921.499.4137.298.662.564(8×8)960.999.6136.599.764.872(8×9)1002.599.3140.9100.366.580(16×5)1021.8118.9164.5114.872.6

圖9 不同特征提取方法分類性能對比Fig.9 Performance of classifier based on different feature extraction methods

圖9(a)可以看出，使用3種分類方法對WVD時頻分布圖進行特征提取，其中NMF的識別率相對較低，2DNMF-G和(2D)2NMF的識別準確率相差不大，各個特征維度的識別準確率都在90%～98%之間，特征維數大于32時2DNMF-Z的識別率稍高，整體來看TD2DNMF特征提取效果最好。但是由于WVD分布中交叉項的干擾，會對故障類別帶來較大影響，所以其識別正確率最高為99.17%，穩定于98%上下。圖9(b)中，3種分類方法對AMP稀疏分解時頻分布圖進行特征提取，識別率與圖9(a)中相比，5種方法均有所提高。這是由于采用AMP對柴油機振動信號進行時頻分析時，生成的時頻分布時頻聚集性更好，各個工況間的差異更明顯，也更利于分類器的分類。TD-2DPCA在特征矩陣維度為8×4和更高維度時，識別準確率穩定在100%，直到特征維數為80時識別率才降為99.83%。對比圖9(a)和(b)可以發現，采用基于AMP稀疏分解與TD2DNMF特征編碼的故障診斷方法更適用于柴油機氣門間隙故障的診斷，并具有較高的診斷精度，且對于一個測試信號樣本而言，整個診斷過程運行時間平均為9.483s，運行效率較高。

為了與圖像方法形成對比，筆者在振動信號的頻域進行故障特征的提取。利用經驗模態分解方法(empirical mode decomposition，簡稱EMD)將信號分解為有限階的內稟模態函數(intrinsic mode function，簡稱IMF)之和，剔除偽分量后保留了前三階IMF。利用AR參數模型求解前三階IMF的功率譜，采用最終預測誤差準則判斷最優階次，使用Yule-Walker方法估計AR模型參數，將每階IMF的前8個AR參量作為特征向量，仍使用NNC作為分類器進行模式識別，訓練集、測試集樣本個數與前文保持一致，平均識別正確率為86.67%。相比而言，振動信號的時頻圖像中既包含了信號時域中各時頻分量的產生及消亡時刻的信息，又包含了頻域中信號的能量分布信息，對柴油機振動信號非平穩時變特征具有較強的刻畫能力。

4 結 論

1) AMP算法利用隨分解殘差信號自適應更新的Gabor字典，通過原子與殘留信號的互相關運算選擇最優原子，與原始MP算法相比，運算速度得到了大幅提高，同時也保持了良好的稀疏性能。

2) TD2DNMF算法將數據矩陣行、列維信息融合到一個判別分析框架中，將不同類別的數據信息并行運算，對柴油機振動時頻圖像樣本進行特征編碼。與現有幾種二維NMF分解算法相比，TD2DNMF算法能更有效地提取柴油機振動時頻圖像的差異特征，計算效率也得到了提高。

3) 利用AMP算法分解振動信號，疊加匹配原子的Wigner-Ville分布所得到的時頻分布具有很高的時、頻分辨率且無交叉項的干擾。用該方法對柴油機不同氣門間隙工況的振動信號進行分析，各工況的時頻分布特征明顯，時頻分量物理意義明確。AMP稀疏分解與TD2DNMF特征編碼相結合的故障診斷方法用于柴油機故障診斷中可獲得較高的診斷精度，是一種有效的機械設備故障診斷方法。