李宏坤，陳禹臻，張志新，周 帥

(1.大連理工大學(xué) 機(jī)械工程學(xué)院，大連 116023;2.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，大連 116023)

時(shí)頻分析方法是對(duì)非平穩(wěn)、非線性信號(hào)進(jìn)行分析的有效方法。機(jī)械振動(dòng)所構(gòu)成的時(shí)頻圖像含有描述機(jī)械工作狀態(tài)的特征信息［1］。為高效率地處理存放于像素矩陣中的數(shù)據(jù)信息，一般采取將矩陣進(jìn)行分解的方法。分解后可以對(duì)原始矩陣中存放的大量數(shù)據(jù)進(jìn)行壓縮和概括，進(jìn)而對(duì)壓縮矩陣進(jìn)行特征提取，得到機(jī)械設(shè)備的狀態(tài)信息。非負(fù)矩陣分解(Non-negative Matrix Factorization，NMF)［2］是目前國際上新的矩陣分解方法，并在理論和應(yīng)用方面已相對(duì)成熟。較完整的NMF理論由Lee等［2］提出，并應(yīng)用于人臉識(shí)別，這引起了科學(xué)界的廣泛關(guān)注。常用的矩陣分解方法有:主成分分析(PCA)、矢量量化(VQ)等。在所有這些方法中，原始的大矩陣V通過V=WH被近似分解為低秩的形式。這些方法的共同特點(diǎn)是:計(jì)算機(jī)的“認(rèn)知”過程是基于數(shù)據(jù)整體的過程，即因子W和H在線性組合時(shí)可進(jìn)行加運(yùn)算，亦可進(jìn)行減運(yùn)算。這樣即使輸入的初始矩陣元素是全正的，分解后的結(jié)果也不能保證非負(fù)性。然而負(fù)值元素在實(shí)際問題中往往是沒有意義的［3］。例如圖像數(shù)據(jù)中不可能有負(fù)值的像素點(diǎn)，在文檔統(tǒng)計(jì)中負(fù)值也是無法解釋的。因此，探索矩陣的非負(fù)分解方法一直是很有意義的研究問題。正是如此，NMF方法才得到人們的廣泛關(guān)注。非負(fù)矩陣分解以計(jì)算機(jī)的“認(rèn)知”過程是基于部分的過程為前提，通過添加“線性組合時(shí)只允許加運(yùn)算不允許減運(yùn)算”的限制條件，保證了分解結(jié)果的可解釋性。與PCA和VQ相比，NMF數(shù)學(xué)實(shí)現(xiàn)簡便、分解形式和分解結(jié)果上的可解釋、以及占用存儲(chǔ)空間少［4］，更重要的是非負(fù)限制與聯(lián)系部分形成整體的直覺概念是協(xié)調(diào)的，這與智慧生命對(duì)事物感知的生理過程是相符的［2］。目前它已廣泛應(yīng)用于諸多領(lǐng)域。

蔡蕾等［5］提出了基于圖像識(shí)別的設(shè)備運(yùn)行狀態(tài)識(shí)別方法，并初步研究了稀疏性非負(fù)矩陣分解中各參數(shù)對(duì)識(shí)別率的影響。本文根據(jù)設(shè)備故障診斷中實(shí)際存在的問題，研究采用Hilbert時(shí)頻譜構(gòu)建圖像，應(yīng)用非負(fù)矩陣分解和主元分析對(duì)其進(jìn)行特征向量提取，完成狀態(tài)的識(shí)別。并以滾動(dòng)軸承在不同狀態(tài)下的識(shí)別為例驗(yàn)證方法的有效性。

1 非負(fù)矩陣分解及PCA降維

1.1 NMF定義及算法實(shí)現(xiàn)

非負(fù)矩陣分解使計(jì)算機(jī)對(duì)整體的認(rèn)知是基于其部分的認(rèn)知。其定義為［2，6］:給定一個(gè)n×m圖像矩陣V，其每一列都是一組圖像數(shù)據(jù)，并且包含n個(gè)非負(fù)像素值。為了估計(jì)原始圖像矩陣，需要尋找兩個(gè)非負(fù)矩陣因子 W∈Rn×r和 H∈Rr×m，使得:

式中:W中的r列數(shù)據(jù)構(gòu)成基本圖像序列組，H中的每一列被認(rèn)為是一個(gè)由構(gòu)成V中每一個(gè)圖像的基本圖像的線性組合系數(shù)組成的編碼且與在V中的每一個(gè)圖像序列一一對(duì)應(yīng)。因數(shù)分解中r值的選擇通常滿足:

WH的乘積即可被認(rèn)為是V中數(shù)據(jù)的壓縮形式。線性組合只允許加法運(yùn)算而不允許減法運(yùn)算的約束將導(dǎo)致分解具有非負(fù)性。

NMF的求解問題實(shí)際上是一個(gè)最優(yōu)化問題，并利用乘性迭代的方法求解W和H。以W和H的非負(fù)性為初始條件，通過使目標(biāo)函數(shù)F收斂至局部最大，得到非負(fù)矩陣V的近似分解［7］。本文使用的目標(biāo)函數(shù)定義為 KL 散度(Kullback-Leibler divergence)［8］:

其中:分解矩陣W和H的更新法則如下:

當(dāng)該目標(biāo)函數(shù)達(dá)到局部最大值時(shí)，得到非負(fù)矩陣V的最優(yōu)近似分解。分解的精確性通過控制Viu與(WH)iu的商值來參與更新過程。更新法則不但維持著W和H的非負(fù)性，而且能有效杜絕W列上數(shù)據(jù)的“退化”現(xiàn)象。

1.2 PCA降維簡介

PCA方法是一種基于二階統(tǒng)計(jì)的數(shù)據(jù)分析方法。該方法通過對(duì)各個(gè)變量之間相關(guān)關(guān)系進(jìn)行分析，用一組較少的、互不相關(guān)的、盡可能多地保留原來復(fù)雜變量所反映信息的主元變量來代替原來較多的變量，從而實(shí)現(xiàn)對(duì)數(shù)據(jù)的降維。這種方法可以有效地找出數(shù)據(jù)中最“主要”的元素和結(jié)構(gòu)，去除噪音和冗余，將原有的復(fù)雜數(shù)據(jù)降維，揭示隱藏在復(fù)雜數(shù)據(jù)背后的簡單結(jié)構(gòu)。它的優(yōu)點(diǎn)是簡單，而且無參數(shù)限制，可以方便的應(yīng)用于各個(gè)場(chǎng)合［9－10］。

2 基于NMF和PCA的狀態(tài)識(shí)別

傅里葉變換對(duì)非平穩(wěn)與非線性信號(hào)進(jìn)行頻譜分析往往得不到理想的效果。因此，為精確的描述這些非穩(wěn)定信號(hào)的局部信息，時(shí)頻分析方法得到廣泛的重視和研究。時(shí)頻方法得到的結(jié)果往往可以用圖像的形式(如輪廓圖像、灰度圖像、浮雕圖像)表達(dá)出來。因此故障診斷實(shí)質(zhì)上是對(duì)圖像的分類:如果圖像分類正確，則可正確的識(shí)別故障。關(guān)于對(duì)圖像的分類，傳統(tǒng)的方法是人工判斷，這將在故障診斷中引入主觀因素。而用計(jì)算機(jī)自動(dòng)處理(如智能識(shí)別系統(tǒng))時(shí)，往往由于圖像數(shù)據(jù)的尺寸和維數(shù)太大而無法進(jìn)行處理。因此，對(duì)圖像在其信息盡可能完整的條件下進(jìn)行降維是非常必要的。

由于圖像像素的非負(fù)性，傳統(tǒng)的降維方法如主成分分析(PCA)、矢量量化(VQ)等方法的應(yīng)用效果往往不是很理想。而加入非負(fù)限定條件的NMF能很好的解決上述問題。根據(jù)設(shè)備故障診斷的實(shí)際特點(diǎn)，本文采用Hilbert時(shí)頻譜構(gòu)建時(shí)頻圖像，應(yīng)用非負(fù)矩陣分解與主元分析方法進(jìn)行特征提取和狀態(tài)識(shí)別。本方法的流程圖如圖1所示。

圖1 基于NMF和PCA的設(shè)備狀態(tài)識(shí)別方法流程圖Fig.1 Flow chart for classification based on NMF and PCA

此方法能夠有效對(duì)包含設(shè)備故障信息的振動(dòng)信號(hào)進(jìn)行表征。采用NMF方法能夠獲取其主要特征向量從而提高識(shí)別的可靠性，PCA能夠進(jìn)行有效降維，便于實(shí)際設(shè)備故障診斷的應(yīng)用。

3 實(shí)驗(yàn)測(cè)試

3.1 實(shí)驗(yàn)數(shù)據(jù)

本文采用的數(shù)據(jù)來自于大連理工大學(xué)振動(dòng)工程研究所的多功能故障診斷試驗(yàn)臺(tái)進(jìn)行的滾動(dòng)軸承故障模擬實(shí)驗(yàn)。模擬的工作狀態(tài)主要有正常、外環(huán)磨損、內(nèi)環(huán)磨損和滾動(dòng)體磨損四種。每個(gè)狀態(tài)共采集了125組數(shù)據(jù)，一組數(shù)據(jù)包含4個(gè)周期，即每個(gè)狀態(tài)包含了500個(gè)周期的數(shù)據(jù)。滾動(dòng)軸承的型號(hào)為N205，轉(zhuǎn)速為1 200 r/min。關(guān)于本實(shí)驗(yàn)對(duì)振動(dòng)信號(hào)采集的要求、獲得單周期信號(hào)Yi(t)的方法以及循環(huán)平均以得到去一階信號(hào)yi(t)的過程見文獻(xiàn)［11］。

3.2 Hilbert時(shí)頻譜分析

Hilbert時(shí)頻譜是利用瞬時(shí)頻率進(jìn)行信號(hào)時(shí)頻分析的方法［12］。將濾去一階的信號(hào)yi(t)進(jìn)行Hilbert時(shí)頻譜分析。為了削弱噪聲的干擾、突出故障信息，可以對(duì)時(shí)頻譜進(jìn)行平均，本文采用5幅Hilbert圖像平均的方法。這樣每種情況將得到100組Hilbert時(shí)頻譜圖像。圖2為平均后不同狀態(tài)下的Hilbert時(shí)頻譜。

圖2 不同狀態(tài)的Hilbert時(shí)頻圖像Fig.2 Different Hilbert Time-frequency image

由圖2可知，各工作狀態(tài)下的Hilbert時(shí)頻圖像并無明顯可區(qū)分特點(diǎn)，因此單單依靠Hilbert時(shí)頻譜是無法進(jìn)行狀態(tài)識(shí)別的。鑒于Hilbert時(shí)頻圖像含有豐富的狀態(tài)信息，為提取有用信息，需對(duì)其進(jìn)行進(jìn)一步的圖像處理。降低;r值選取過大時(shí)，則可能帶入多余的冗余信息和噪聲。當(dāng)r值為45時(shí)，平均識(shí)別成功率達(dá)到最大［13］。因此本文r選擇45。則每種情況將得到100組維數(shù)為45的特征向量ei。

4 基于NMF和PCA的狀態(tài)識(shí)別

4.1 特征向量的提取

在圖像處理之前，要對(duì)Hilbert時(shí)頻圖像(為201×479的矩陣)進(jìn)行真彩化和歸一化處理(這樣可以減小特征向量不同量級(jí)之間差異的影響，提高訓(xùn)練和預(yù)測(cè)的準(zhǔn)確性)。

將真彩化和歸一化處理后的矩陣(201×479)表示成一個(gè)96279維的列向量。每種狀態(tài)(正常、外環(huán)磨損、內(nèi)環(huán)磨損和滾動(dòng)體磨損)都有100個(gè)圖像數(shù)據(jù)。由100組96 279維的列向量構(gòu)成96279×100矩陣V作為非負(fù)矩陣。經(jīng)非負(fù)矩陣分解后，得到矩陣因子W和H，滿足V≈WH。W的列向量可認(rèn)為是組成原圖像矩陣V的基本圖像。將代表信號(hào)圖像的96 279維列向量vi與W進(jìn)行線性組合，即ei=WTvi。這樣，r維向量ei即可認(rèn)為是該幅信號(hào)圖像vi的特征向量。其中r的值滿足式(2)，即 r＜99.896。

非負(fù)矩陣分解r值的選取至關(guān)重要。r為描述原始圖像矩陣基本圖像的個(gè)數(shù)以及得到的特征向量的維數(shù)。r值選取過小時(shí)，得到的基本圖像不足以描述原始圖像矩陣，因而向基投影后部分信息丟失，導(dǎo)致識(shí)別率

4.2 利用歐式距離分離各種狀態(tài)

歐式距離是圖像識(shí)別中經(jīng)常采用的方法，因此本文首先采用歐式距離對(duì)不同狀態(tài)下的故障進(jìn)行判別。從上述四種情況的100個(gè)特征向量中隨機(jī)各選一個(gè)向量作為參考向量，分別命名為 ball，inner，normal，outer。將特征向量組矩陣NMF(NMF矩陣由45×400構(gòu)成，其中1～100列、101～200列、201～300列、301～400列分別是滾動(dòng)體磨損、內(nèi)環(huán)磨損、正常、外環(huán)磨損的特征向量)的每一列與上述四個(gè)參考向量分別作歐式距離。當(dāng)對(duì)應(yīng)情況的歐式距離最小時(shí)則識(shí)別成功，圖3為識(shí)別結(jié)果。

圖3 歐式距離識(shí)別結(jié)果Fig.3 Classification result by using euclidean distance

從圖3可以看出，除對(duì)正常情況識(shí)別成功以外，其余都不理想。這說明雖然特征向量的維數(shù)較大時(shí)，可以較全面的反應(yīng)原始圖像的特征信息，但維數(shù)的增加并不利于利用特征向量進(jìn)行狀態(tài)識(shí)別。所以在狀態(tài)識(shí)別之前對(duì)特征向量進(jìn)行降維變得非常必要。

4.3 基于PCA的降維處理

由于PCA的諸多優(yōu)點(diǎn)(如前所述)，本文采用PCA方法對(duì)特征向量進(jìn)行降維處理。經(jīng)過PCA處理，將NMF處理得到的400組維數(shù)為45的特征向量ei降維成維數(shù)為3的特征向量ei'。通過上述的數(shù)據(jù)處理，得到正常、外環(huán)磨損、內(nèi)環(huán)磨損和滾動(dòng)體磨損四種情況分別對(duì)應(yīng)的特征向量組。每種情況對(duì)應(yīng)100組維數(shù)為3的特征向量。將降維后的特征向量組矩陣命名為PCA(PCA矩陣由3×400構(gòu)成，其中1～100列、101～200列、201～300列、301～400列分別是滾動(dòng)體磨損、內(nèi)環(huán)磨損、正常、外環(huán)磨損降維后的特征向量)。本文將一個(gè)三維特征向量當(dāng)作三維空間中的一個(gè)點(diǎn)并利用三維坐標(biāo)系來分離不同的狀態(tài)，圖4、圖5為特征向量在三維坐標(biāo)系上的反映?？梢钥闯觯?、外環(huán)磨損、內(nèi)環(huán)磨損和滾動(dòng)體磨損四種狀態(tài)被識(shí)別的非常清楚。正常、外環(huán)磨損和滾動(dòng)體磨損的分布區(qū)域依次鏈狀分布，而內(nèi)環(huán)磨損的點(diǎn)則成環(huán)狀圍繞在鏈的周圍。因而解決了只采用Hilbert時(shí)頻譜識(shí)別故障狀態(tài)不理想的問題，也解釋了歐式距離無法識(shí)別狀態(tài)的原因。

因此經(jīng)過主元分析之后可以有效獲取識(shí)別向量作為參數(shù)進(jìn)行設(shè)備的狀態(tài)識(shí)別，此方法可以為以后的狀態(tài)辨識(shí)研究打下堅(jiān)實(shí)的基礎(chǔ)。

5 結(jié)論

本文研究了基于NMF的時(shí)頻圖像特征提取與PCA降維進(jìn)行振動(dòng)信號(hào)特征提取與狀態(tài)識(shí)別的方法，并采用滾動(dòng)軸承信號(hào)的故障診斷進(jìn)行方法的有效性驗(yàn)證。通過對(duì)比，此方法在識(shí)別滾動(dòng)體磨損等沖擊現(xiàn)象不明顯的故障狀態(tài)時(shí)比Hilbert時(shí)頻譜更為有效。此方法的研究將有助于設(shè)備故障診斷和預(yù)知維修的發(fā)展。

［1］ Newton R E.Fragility assessment theory and practice［R］.Monterey Research Laboratory， Inc. Monterey，Califonia，1968.

［2］ Lee D D，Seung H S.Learning the parts of objects by nonnegative matrix factorization ［J］.Nature，1999，(401):788－791.

［3］ Paatero P，Tapper U.Positive matrix factorization:a nonnegative factor model with optimal utilization of error estimates of data values［J］.Environmetrics，1994，5(2):111－126.

［4］ 汪 鵬.非負(fù)矩陣分解:數(shù)學(xué)的奇妙力量［J］.計(jì)算機(jī)教育，2004(10):38－40.

［5］ 蔡 蕾，朱永生.基于稀疏性非負(fù)矩陣分解和支持向量機(jī)的時(shí)頻圖像識(shí)別［J］.自動(dòng)化學(xué)報(bào)，2009，35(10):1272－1277.

［6］ Lee D， Seung H. Algorithmsfornon-negative matrix factorization［M］.USA:NeuralInformationProcessing Systems(NIPS).Denver，CO，2000，7.

［7］ Lee D D，Seung H S.Algorithms for non-negative matrix factorization［C］.Advances in Neural Information Processing Systems.Cambridge:MTT Press，2000:556－562.

［8］ 李 樂，章毓晉.非負(fù)矩陣分解算法綜述［J］.電子學(xué)報(bào)，2008，36(4):737－743.

［9］ Jolliffe I T.Principal component analysis(2nded)［M］.New York:Springer，2002.

［10］ He Q B，Yan R Q，Kong F R，et al.Machine condition monitoring using principal component representations［J］.Mechanical Systems and Signal Processing，2009，23(2):446－466.

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

［12］ Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.Proceedings of The Royal Society of London Series A-Mathematical Physical and Engineering Sciences，1998，454:903－ 995.

［13］ 宦若虹，楊汝良.基于小波域NMF特征提取的SAR圖像目標(biāo)識(shí)別方法［J］.電子與信息學(xué)報(bào)，2009，31(3):588－591.