岳應(yīng)娟，王 旭，蔡艷平，牟偉杰

（1.火箭軍工程大學理學院，西安 710025； 2.火箭軍工程大學五系，西安 710025）

前言

柴油機是非常復(fù)雜的機械系統(tǒng)，所測的振動信號屬于典型的非線性、非平穩(wěn)信號?；谡駝有盘柕牟裼蜋C故障診斷的本質(zhì)是模式識別問題，其中，如何從振動信號中準確提取出故障特征是問題關(guān)鍵。

時頻分析的方法能很好地反映信號時間和頻率的局部化信息，處理非線性、非平穩(wěn)的振動信號時效果較好[1]。采用時頻分析是將獲得的時域信號轉(zhuǎn)化為時頻圖像，時頻圖像提供了足夠的信號故障分類信息，可通過分析圖像實現(xiàn)故障的特征提取和智能診斷[2]。如文獻[3]中通過灰度共生矩陣提取出時頻圖像的紋理特征；文獻[4]中通過提取時頻圖像的三維質(zhì)心與信息熵來進行設(shè)備狀態(tài)識別；文獻[5]中提出從時頻圖像的形狀特征、灰度統(tǒng)計特征和紋理特征來綜合提取特征參數(shù)的思想；文獻[6]中提出對時頻圖像進行分割并選取特征體質(zhì)心位置、面積、數(shù)目和熵作為特征參數(shù)。

然而，在利用圖像分析方法對時頻圖像進行特征提取時，大多是提取圖像的紋理特征、灰色統(tǒng)計特征或形狀特征。這些方法都忽略了時頻圖像中的局部特征關(guān)系，難以區(qū)分出在不同時刻具有相同頻率分量的時頻圖像。若將時頻圖像矩陣元素直接構(gòu)成特征向量進行模式識別，又會由于矩陣數(shù)據(jù)維數(shù)過大帶來“維數(shù)災(zāi)難”，故有必要進一步獲取蘊含在時頻圖像內(nèi)部的低維特征。非負矩陣分解算法(nonnegative matrix factorization，NMF)[7]是一種局部特征提取方法，能將存放數(shù)據(jù)的矩陣維數(shù)進行約簡并保留關(guān)鍵信息。文獻[8]和文獻[9]中將NMF應(yīng)用于齒輪和軸承時頻圖像矩陣的低維局部特征提取，取得了很好的效果。但NMF算法存在的問題是計算效率過低，一定程度上限制了其在機械故障診斷領(lǐng)域的應(yīng)用。文獻[10]中提出了一種分塊非負矩陣分解算法(block NMF，BNMF)，有效提高了NMF算法的計算效率。但分塊運算將時頻矩陣分解為若干塊小矩陣，分別運算后再進行合成，特征矩陣的稀疏性又難以保證，不利于故障類型的識別。

針對以上問題，本文中提出一種新的時頻圖像特征約簡方法(記為C-BSNMF方法)。利用三次卷積插值法將振動譜圖像進行維度壓縮，通過矩陣分塊計算方法結(jié)合稀疏非負矩陣分解，將壓縮后的圖像矩陣分解為一系列代表了不同時間和頻率范圍內(nèi)的時頻分量的基圖像的線性組合，直接提取出時頻圖像中的故障特征。實驗驗證表明，該方法能夠準確地提取出柴油機氣門機構(gòu)的故障特征，并且大幅提高了特征提取的計算效率。

1 分塊稀疏非負矩陣算法

矩陣的稀疏計算由向量的稀疏度計算推演而來。給定一個向量x∈Rn，文獻[11]中提出使用L1范數(shù)與L2范數(shù)之比作為該向量的稀疏度測度:

一個向量的稀疏度介于0～1之間。圖1給出了對稀疏度的直觀理解，稀疏度越大，特征量就越少，主要特征就越明顯。

圖1 不同稀疏度下的稀疏向量

1.1 稀疏非負矩陣分解

SNMF是在稀疏編碼和NMF相結(jié)合的基礎(chǔ)上提出的[12]。SNMF通過對NMF中的基矩陣或系數(shù)矩陣施加稀疏約束來進一步減少約簡后信息的冗余。SNMF的本質(zhì)是一種矩陣分解和投影技術(shù)，其基本原理如下。

對于非負觀測矩陣V，找到兩個合適的非負矩陣因子W和H，滿足:

式中:n為數(shù)據(jù)樣本的維數(shù)；m為集合中數(shù)據(jù)樣本的個數(shù)；r為特征維數(shù)，一般情況下，r遠小于n且滿足f(i＋u，j＋v)。

為了描述V≈W·H的近似效果，利用矩陣V與W·H間的K－L散度作為近似誤差，在NMF的基礎(chǔ)上加入了稀疏懲罰項α對應(yīng)的目標函數(shù):

其對應(yīng)的優(yōu)化問題為

根據(jù)此目標函數(shù)，文獻[11]中推導(dǎo)出了W和H的更新規(guī)則，推導(dǎo)過程參照經(jīng)典的NMF基本算法，也即當α＝0時SNMF退化為經(jīng)典NMF算法。W和H的迭代公式為

1.2 分塊稀疏非負矩陣分解

分塊計算的思想最早在文獻[13]中提出。雖然SNMF算法符合“局部構(gòu)成整體”的認知規(guī)律，但它開始是將矩陣看成整體而非局部之和。如果從開始就將矩陣視為“整體為局部之和”，對矩陣分塊各自分解后再將分解后的各子塊按原序拼接，而分解過程中矩陣的各塊之間并不運算，則可以有效提高算法的收斂性能。本文中將這種分塊運算方法應(yīng)用于SNMF，提出一種分塊稀疏非負矩陣分解算法(block SNMF，BSNMF)，其基本原理如下。

將非負觀測矩陣按列平均分為一系列塊矩陣的形式 Vm×n＝[V1V2… Vb]，其中 Vi∈Rm×n0，b 為分塊數(shù)，n0＝n/b。對每一個塊矩陣按照式(5)～式(7)進行迭代運算，得到相應(yīng)的基矩陣和系數(shù)矩陣。

式中r0＝r/b為塊矩陣的特征維數(shù)。按照式(9)和式(10)將分解得到的基矩陣和系數(shù)矩陣按原順序進行合成，即可得到原非負觀測矩陣的基矩陣和系數(shù)矩陣。

1.3 算法計算復(fù)雜度分析

以算法中的乘法次數(shù)來表征計算復(fù)雜度。由式(5)，在矩陣H的迭代計算中，SNMF算法的計算復(fù)雜度為

由式(6)和式(7)，在矩陣W的迭代計算中，SNMF算法的計算復(fù)雜度為

一次迭代過程SNMF算法的計算復(fù)雜度:

對于本文中所提出的BSNMF算法，在Hi與Wi的迭代計算中，計算復(fù)雜度分別為

則一次迭代中BSNMF算法的計算復(fù)雜度為

比較式(13)與式(16)可知，BSNMF的計算量要比SNMF的計算量要小得多。

2 三次卷積插值法

雖然BSNMF算法可有效減少SNMF算法的計算量，但由于分塊計算時要滿足r0＝r/b，使r0為正整數(shù)，故分塊數(shù)b也要滿足(m＋n)b＜mn，對于較高維度的振動譜圖像來說，算法耗時仍舊很長。所以，考慮利用三次卷積插值法對每個振動譜圖像進行維度壓縮。

三次卷積插值[14]又稱為立方卷積插值，其本質(zhì)是利用采樣點的周圍16個點的灰度值作三次內(nèi)插運算，如圖2所示。

圖2 待插值點與16個鄰近點排布

由連續(xù)信號采樣定理可知，若對采樣值用插值函數(shù)sinc(w)＝sin(w)/(w)進行插值，則可準確地恢復(fù)原函數(shù)，得到采樣點間任意點的值。三次卷積插值法實質(zhì)上就是利用一個三次多項式來近似理論上最佳插值函數(shù)S(w):

目標值 f(i＋u，j＋v)可由如下插值公式得到:

其中:

三次卷積插值法帶有邊緣增強的效果，能較好地保持圖像的細微結(jié)構(gòu)。因此利用三次卷積插值算法與BSNMF相結(jié)合，能夠?qū)Ω呔S的振動譜圖像矩陣進行有效的特征約簡，同時具有較高的執(zhí)行效率，使算法的泛化性能得到進一步增強。

3 柴油機氣門故障特征提取實例

3.1 實驗流程

基于時頻圖像特征約簡方法進行柴油機故障特征提取和效果驗證的流程如圖3所示。

圖3 故障特征提取及效果驗證流程

實驗對比方法設(shè)置如下。

(1)時頻圖像生成，分別用線性時頻分析中具有代表性的短時傅里葉變換(STFT)和非線性時頻分析中的魏格納分布(WVD)、偽魏格納分布(PWVD)、平滑偽魏格納分布(SPWVD)、喬威廉姆斯分布(CWD)以及參數(shù)化時頻分析方法中利用匹配追蹤算法(MP)與WVD相結(jié)合生成的時頻分布，共6種時頻分析方法產(chǎn)生振動譜圖像來驗證本文方法。以工況1的時頻分布為例，用6種分析方法生成時頻圖像，為便于進一步的圖像分類與識別，將圖像轉(zhuǎn)化為維度為420×560的灰度圖的形式，如圖4所示。

圖4 工況1的6類時頻表示振動譜圖

STFT描述的是信號的頻譜隨時間變化的情況，分辨率自適應(yīng)性差；其余5類分布都是描述的信號的能量密度隨時間變化的情況，WVD具有最佳的時頻分辨率，但是交叉干擾項嚴重，PWVD，SPWVD和CWD是通過加窗或引入核函數(shù)的方法來抑制WVD的交叉干擾項，但都不可避免地降低了時頻分辨率；利用MP算法將信號分解為一系列單分量原子，再疊加每個原子分量的WVD，得到的時頻分布保持了WVD良好的時頻聚集性，并且消除了WVD中交叉項的干擾。

(2)特征提取，本文方法運用三次卷積插值算法對時頻圖像矩陣進行維度壓縮，用分塊稀疏非負矩陣分解算法對時頻表示特征約簡(C-BSNMF特征提取)。將這種方法與現(xiàn)有的NMF，BNMF和SNMF 3種方法進行對比。以特征提取時間來表征計算時間復(fù)雜度，對比方法的計算效率；以分類器的識別率來表征特征提取效果，對比方法的特征提取性能。分類器分別選擇模式識別性能較為穩(wěn)定的K-最近鄰分類器(KNNC)和在處理小樣本、非線性數(shù)據(jù)時效果較好的支持向量機(SVM)。

3.2 實驗設(shè)置

以6135型柴油機為研究對象，實驗平臺由柴油機、傳動軸、電機和控制臺4部分組成，如圖5所示。取柴油機第2缸蓋表面振動信號對柴油機進行故障診斷，采樣頻率25kHz，轉(zhuǎn)速為1 500r/min，測試過程中，柴油機空載運行。共設(shè)置了4種氣門間隙狀況，具體情況如表1所示。其中0.06，0.3和0.5mm分別對應(yīng)氣門間隙過小、正常和過大，開口表示在氣門上開4mm×1mm孔來模擬嚴重漏氣故障。共采集柴油機氣門4種故障狀態(tài)下各60組振動信號樣本，總計240個。

圖5 實驗平臺

3.3 實驗工況分析

根據(jù)6135型柴油機的工作原理可知，引起缸蓋振動的原因主要是氣體燃燒時產(chǎn)生的爆壓，氣門落座撞擊和排氣門開啟所引起的氣流沖擊等，鄰缸的各種振動激勵源也會對缸蓋的振動產(chǎn)生較大影響[15]。圖6示出進排氣門開閉和點火時刻與曲軸轉(zhuǎn)角的關(guān)系。進氣門開啟的角度在排氣上止點前20°CA附近，關(guān)閉的角度在進氣下止點后48°CA附近；排氣門開啟的角度在做功下止點前48°CA附近，關(guān)閉的角度在進氣下止點后20°CA附近；柴油機在0°CA點火。圖7和圖8給出了4種工況下振動信號的時域和功率譜圖像。

表1 4種實驗工況設(shè)置 mm

圖6 柴油機燃燒和氣門開閉轉(zhuǎn)角圖

圖7 4類工況時域波形

圖8 4類工況功率譜圖

從圖7的時域波形圖可以看出，氣門間隙狀態(tài)的變化對振動信號有較大的影響，主要體現(xiàn)在幾個劇烈振動分量產(chǎn)生與消亡的時刻發(fā)生了變化，各劇烈振動分量的幅值變化也特別明顯。比較圖8中4種狀態(tài)下的功率譜圖可以看出，在正常狀態(tài)下，缸蓋振動信號的能量主要集中在6.5～8.5kHz之間的高頻區(qū)域。在非正常狀態(tài)，也即氣門間隙出現(xiàn)故障時，主要能量向更高的頻段移動，一般集中在8.5～12.0kHz之間。

3.4 實驗結(jié)果與分析

將每種工況的60組振動信號中隨機選出30組信號組成訓練樣本，其余信號組成測試樣本。第3.1節(jié)中所述6類時頻表征方法生成圖像的維度為420×560，利用三次卷積插值算法對時頻矩陣進行0.1倍壓縮，得到維度為42×5642×56的6類時頻矩陣，如圖9所示。由于實際應(yīng)用中估計低維空間維數(shù)R十分困難，且依據(jù)應(yīng)用目的不同維數(shù)設(shè)置準則也存在差別[16]。同時，受分塊計算算法的限制，根據(jù)實際的矩陣維度情況，以16為步長，R從16到128取值，在不同的特征維數(shù)下分別迭代計算，并限定總的迭代次數(shù)為160(分塊運算時每一個塊矩陣迭代10次)。

表2給出了4種算法對WVD時頻圖像進行特征提取的計算效率，均不包含圖像載入時間。實驗環(huán)境為Matlab R2012b，AMD A8處理器，CPU主頻1.90GHz，內(nèi)存為4GB，操作系統(tǒng)為Windows 7。

圖9 利用三次卷積插值算法壓縮后的6類時頻圖像

表2 4種方法對WVD時頻特征提取的耗時情況s

由表可見:從計算耗時的角度分析，縱向來看，用同一種矩陣分解方法時，計算耗時隨選取的特征維數(shù)的增大而增加；橫向來看，分塊計算的方法能夠提高運算效率，分塊數(shù)為16時，BNMF的運算效率較NMF提高了13倍左右，印證了第1.3節(jié)計算復(fù)雜度的推導(dǎo)。通過運用三次卷積插值法對時頻圖像進行維度壓縮，本文中所提的 C-BSNMF方法在SNMF的基礎(chǔ)上，運算速度提高了1 500～2 000倍。

表3給出了不同特征維數(shù)下，運用KNNC與SVM分別作為分類器，4種方法對WVD時頻圖像特征提取的識別效果。將生成的6種振動時頻圖像分布用4種約簡方法進行處理，并用KNNC和SVM兩種分類器進行模式識別，得到不同時頻圖像應(yīng)用不同約簡方法的識別結(jié)果，如圖10所示。

圖10 6類時頻圖像識別正確率比較

由圖10可以看到，對于同樣的一批振動數(shù)據(jù)，在同樣的分類指標下，采用的時頻分析方法不同，其識別正確率也不一樣。其中，以MP-WVD和WVD得到的時頻分布的識別正確率相對較高，最高識別率都達到了100%，而STFT和CWD分布的識別正確率相對較低。這表明，時頻分析方法對識別結(jié)果有一定的影響。WVD分布具有最好的時頻聚集特性，MP-WVD算法得到的時頻分布是將單分量原子信號WVD線性疊加，也保持了最好的時頻聚集特性。能量本身是一種二次型的表示，這種描述具有更好的邊緣特性和有限支撐性質(zhì)。這種良好的時頻聚集特性使各狀態(tài)間的差異性相對更為明顯。不論是用SVM還是KNNC作為分類器，都能保證很高的識別正確率。相比下，PWVD，SPWVD和CWD 3種分布雖然抑制了WVD分布的交叉項的干擾，但也降低了時頻分布的聚集性能，利用本文方法得到的特征提取效果反而不如WVD。

表3 4種方法對WVD時頻圖像特征提取效果對比

對比分析NMF，BNMF，SNMF和 C-BSNMF 4種方法的特征提取效果。以SVM作為分類器時，可以看到，本文所提方法C-BSNMF的效果最優(yōu)(采用CWD時頻分布時除外)；以KNNC作為分類器時，本文所提方法的效果也最優(yōu)(采用PWVD時頻分布時除外)，充分說明了本文方法更易獲取蘊含在高維時頻矩陣內(nèi)部的低維特征。

4 結(jié)論

(1)將基于時頻圖像的特征約簡提取方法應(yīng)用于6135型柴油機氣門機構(gòu)的故障自動診斷中，取得了較高的識別正確率。說明該方法對于柴油機的復(fù)雜非平穩(wěn)振動信號的時頻特征提取是有效的。

(2)結(jié)合三次卷積插值法與分塊稀疏非負矩陣分解算法的特征約簡方法在低維特征提取效果和計算效率上較現(xiàn)有的幾種特征約簡方法有很大程度的提高，為非線性、非平穩(wěn)信號的特征提取提供了一種有效的方法。

(3)利用時頻圖像特征約簡的方法提取柴油機氣門機構(gòu)故障特征時，基于匹配追蹤算法的時頻分布和魏格納時頻分布的特征提取效果要優(yōu)于其他幾種時頻表示方法，說明利用本文時頻圖像良好的時頻聚集特性能使特征提取效果更優(yōu)。此結(jié)果對時頻分析方法的研究方向有一定的借鑒意義。

[1] 張俊紅，王健，畢鳳榮，等.基于EMD和時頻分析的低振動機體結(jié)構(gòu)優(yōu)化研究[J].振動與沖擊，2014，33(3):117－121.

[2] FENG Zhipeng， LIANG Min， CHU Fulei.Recent advances in time-frequency analysis methods for machinery fault diagnosis:a review with application examples[J].Mechanical Systems and Signal Processing，2013，38(1):165－205.

[3] 劉建敏，劉遠宏，江鵬程，等.基于包絡(luò)S變換時頻圖像提取齒輪故障特征[J].振動與沖擊，2014，33(1):165－169.

[4] 李宏坤，周帥，黃文宗.基于時頻圖像特征提取的狀態(tài)識別方法研究與應(yīng)用[J].振動與沖擊，2010，29(7):184－188.

[5] 蔡艷平，李艾華，何艷萍，等.基于振動譜時頻圖像特征及SVM參數(shù)同步優(yōu)化識別的內(nèi)燃機故障診斷[J].內(nèi)燃機學報，2012，30(4):377－383.

[6] 蔡艷平，李艾華，王濤，等.基于時頻譜圖與圖像分割的柴油機故障診斷[J].內(nèi)燃機學報，2011，29(2):181－186.

[7] LEE D D，SEUNG H S.Learning the parts of objects by non-negative matrix factorization[J].Nature，1999，401(21):788－791.

[8] 李宏坤，陳禹臻，張志新，等.基于非負矩陣分解與主元分析的時頻圖像識別研究方法[J].振動與沖擊，2012，31(18):169－172.

[9] LI Bing， ZHANG Peilin， TIAN Tao， et al.A new feature extraction and selection scheme for hybird fault diagnosis of gearbox[J].Expert Systems with Applications，2011，38(8):10000－10009.

[10] 潘彬彬，陳文勝，徐晨.基于分塊非負矩陣分解人臉識別增量學習[J].計算機應(yīng)用研究，2009，26(1):117－120.

[11] HOYER P O.Non-negative matrix factorization with sparseness constraints[J].Journal of Machine Learning Research，2004(5):1457－1469.

[12] LIU Weixiang， ZHENG Nanning， LU Xiaofeng.Non-negative matrix factorization for visual coding[C].Washington:Proc of IEEE Int Conf on Acoustics， Speech and Signal Processing，2003:293－296.

[13] SRA S.Block-iterative algorithms for non-negative matrix approximation[C].Eighth IEEE International Conference on Data Mining，2008:1037－1042.

[14] ZHOU D，SHEN X，DONG W.Image zooming using directional cubic convolution interpolation[J].IET Image Processing，2012，6(6):627－634.

[15] WANG Qinghua， ZHANG Youyun， CAI Lei， et al.Fault diagnosis for diesel valve trains based on non-negative matrix factorization and neural network ensemble[J].Mechanical Systems and Signal Processing，2009，23(5):1683－1695.

[16] GAO Huizhong， LIANG Lin， CHEN Xiaoguang， et al.Feature extraction and recognition for rolling element bearing fault utilizing short-time Fourier transform and non-negative matrix factorization[J].Chinese Journal of Mechanical Engineering，2015，28(1):96－105.