唐曦凌，梁 霖，高慧中，羅愛玲

(1.廣東電網(wǎng)公司電力科學(xué)研究院，廣州 510080;2.西安交通大學(xué) 機(jī)械工程學(xué)院，西安 710049)

在機(jī)械故障診斷中，能否準(zhǔn)確地提取故障特征是診斷的關(guān)鍵所在。由于機(jī)械振動信號常表現(xiàn)為復(fù)雜的非平穩(wěn)信號，且早期故障特征往往會被大量噪聲淹沒，這為故障特征提取帶來了很大的困難。

小波變換是目前分析非平穩(wěn)信號的一種有效工具，而尺度連續(xù)變化的連續(xù)小波變換能夠細(xì)致刻劃信號的局部形態(tài)，采用與特征波形類似的小波基對信號進(jìn)行連續(xù)小波變換，可通過時頻分布圖觀測到故障特征的周期性變化。但是由于機(jī)械設(shè)備結(jié)構(gòu)復(fù)雜，采集的振動信號復(fù)雜，且含有大量噪聲，僅通過時頻分布很難提取早期故障特征，且需要先驗(yàn)知識，對人員經(jīng)驗(yàn)依賴性強(qiáng)，容易誤判。

研究表明，故障特征信息主要存在于變換系數(shù)的部分尺度中，這使得變換的系數(shù)矩陣具有一定的稀疏性，同時，噪聲也包含在系數(shù)矩陣中。因此，有學(xué)者提出利用SVD對系數(shù)矩陣進(jìn)行特征約簡［1－2］，并提取故障特征。然而SVD只能提取系數(shù)矩陣中的線性特征，當(dāng)系數(shù)矩陣具有非線性分布時，則不能很好獲得其低維本質(zhì)結(jié)構(gòu)。

而非負(fù)矩陣分解(Nonnegative Matrix Factorization，NMF)為提取數(shù)據(jù)非線性特征提供了一個新思路。它是一種以非負(fù)性為約束的矩陣分解方法，能夠?qū)崿F(xiàn)非線性降維和基于局部的特征提取，且分解結(jié)果更稀疏。該方法已逐漸被應(yīng)用于多個領(lǐng)域，如Lee和Seung將NMF用于機(jī)器人對外界環(huán)境的感知，用來提取有意義的局部特征［3］;Kawamoto運(yùn)用NMF分解音樂信號的幅度譜，最終得到單個音調(diào)的特征［4］;蔡蕾等［5］結(jié)合SNMF和支持向量機(jī)，對機(jī)械運(yùn)行狀態(tài)的時頻圖進(jìn)行識別，進(jìn)而對機(jī)器運(yùn)行狀態(tài)進(jìn)行判斷。文獻(xiàn)［6］針對機(jī)械故障信號，利用正交性約束NMF在頻域?qū)崿F(xiàn)故障信號特征提取。

現(xiàn)有研究成果表明，由于非負(fù)性是較弱的約束性，這使得NMF分解結(jié)果常常并不唯一，因此需要根據(jù)實(shí)際應(yīng)用中加入其他約束，從而獲得滿意的唯一解。

本文針對早期微弱故障的特征提取問題，根據(jù)故障特征在小波分解系數(shù)矩陣中的稀疏性以及信噪比較低的特點(diǎn)，將基于稀疏性和光滑性約束的多約束NMF(multiple constraints NMF，mcNMF)方法應(yīng)用于連續(xù)小波變換(Continuous Wavelet Transform，CWT)的系數(shù)矩陣的分析中，提出了結(jié)合連續(xù)小波變換和多約束非負(fù)矩陣分解的特征提取方法CWT－mcNMF。通過仿真數(shù)據(jù)和實(shí)際案例表明該方法的有效性。

1 連續(xù)小波變換

設(shè)x(t)是具有有限能量的信號，則連續(xù)小波變換定義為x(t)與小波基函數(shù)的內(nèi)積［7］，即:

式中:Wx(τ，a)為小波變換系數(shù)矩陣;τ為時間定位參數(shù);a 為尺度參數(shù)，a ＞0;ψ(t)為小波函數(shù);ψτ，a(t)為小波基函數(shù)，是由小波母函數(shù)ψ(t)進(jìn)行伸縮和平移而形成的一組函數(shù)系列，即:

連續(xù)小波變換通過時間和尺度的連續(xù)變化將信號投影到時間—尺度平面上，系數(shù)矩陣Wx(τ，a)度量了信號與小波間的相似程度，反映了信號的特征信息。

2 非負(fù)矩陣分解方法原理

對一個T維的隨機(jī)向量v進(jìn)行了t次觀測，得到觀測矩陣 V=［v1，v2，…，vI］T，非負(fù)矩陣分解的問題可描述為:已知非負(fù)矩陣T，找出非負(fù)的W和非負(fù)的H，使式(3)成立［8］。

式中:V稱為觀測矩陣，是I×T的非負(fù)矩陣;W稱為基矩陣，是I×J的非負(fù)矩陣;H稱為系數(shù)矩陣，是J×T的非負(fù)矩陣。通常情況下I?J，這表示可對觀測矩陣進(jìn)行維數(shù)約簡，H則表示特征約簡后的矩陣。

為了解決該問題，目前的算法中，標(biāo)準(zhǔn)的ALS算法較其他算法，效率更高，更適合大數(shù)據(jù)的處理。但該方法對噪聲敏感，常獲得次優(yōu)解。而利用被處理信號的先驗(yàn)知識，對W和H施加更多的限制條件，如稀疏性、光滑性、正交性等，可獲得較滿意的解。加入相應(yīng)的約束性后，非負(fù)矩陣分解最優(yōu)化問題變?yōu)椋?］:

式中:·F—矩陣的 Frobenius范數(shù);αW為正則化系數(shù)，αW≥0;JW為對W施加的約束條件;αH為正則化系數(shù)，αH≥0;JH為對H施加的約束條件。

3 基于最優(yōu)小波系數(shù)的多約束NMF特征提取

3.1 特征提取思路

由于通過連續(xù)小波變換，使信號的故障特征信息的變換系數(shù)與其他分量不同，并具有稀疏性，另外，變換系數(shù)中仍包含噪聲分量，因此在利用非負(fù)矩陣分解對系數(shù)矩陣進(jìn)行分析時，為了提取故障特征分量，可對結(jié)果，即H施加稀疏性和光滑性約束:施加稀疏性約束，可獲得具有稀疏性的特征分量;施加光滑性約束，可提高特征分量的信噪比。最終實(shí)現(xiàn)系數(shù)矩陣的非線性降維;因缺乏W的先驗(yàn)知識，所以無法施加其他的約束性，本文采用以局部歐式距離為目標(biāo)函數(shù)的HALS算法進(jìn)行計(jì)算［10］。

小波變換中，為了最大程度的與故障特征匹配，需要根據(jù)不同的故障類型選擇合理的小波基，變換后的系數(shù)矩陣才能準(zhǔn)確的反映故障特征。在眾多故障中，沖擊類故障占有較高的比例，如缺陷軸承和齒輪產(chǎn)生的瞬時沖擊等，因此本文針對沖擊類故障，采用CWT－mcNMF方法提取微弱故障特征。

3.2 算法步驟

具體步驟如下:

步驟1:最優(yōu)系數(shù)矩陣的獲取。選擇與沖擊故障能夠較好匹配的Morlet小波［11］，對于機(jī)械振動信號，常采用其實(shí)部作為小波基函數(shù)，即:

式中:β為波形參數(shù)，控制小波波形的衰減率。為了達(dá)到最優(yōu)匹配，采用小波熵［12］確定波形參數(shù)β，并獲得最優(yōu)系數(shù)矩陣 Wβop(τ，a)。

步驟2:非負(fù)化。將系數(shù)矩陣每一行都加上一直流分量，使矩陣為非負(fù)。

步驟3:初始化。隨機(jī)初始化W和H，輸入約簡維數(shù)J。

步驟4:分別迭代計(jì)算W和H，直到收斂。

(1)對H的迭代規(guī)則

對H加入稀疏性約束和平滑性約束的最優(yōu)化問題如下:

式中:‖H‖2F=∑i‖hi‖22是對H施加光滑性約束，‖H‖l=∑jt|hjt|是對H施加稀疏性約束。

采用ALS獲得H的迭代規(guī)則為:

式中:［x］+=max(ε，x)。

對帶電粒子進(jìn)行如下分析，帶電粒子受到一個豎直向下的恒力電場力的作用，并且射入平行板電容器時帶有一個水平的初速度v0。帶電粒子在水平方向上做的是勻速直線運(yùn)動，在豎直方向上做的是初速度為零的勻變速直線運(yùn)動。這種運(yùn)動類似于物體在重力作用下只受到一個水平初速度的平拋運(yùn)動。所以可根據(jù)平拋運(yùn)動的相關(guān)知識來解有關(guān)題目，可以求得以下有關(guān)物理量或相關(guān)結(jié)論：

(2)對W的迭代規(guī)則

采用HALS對W進(jìn)行迭代，其獲得的迭代規(guī)則為:

式中:V(j)=E+wjhTj，另外，每次迭代需對W的每一列進(jìn)行歸一化處理:wij←wij/‖wj‖2。

步驟5:獲得特征約簡結(jié)果H。

在上述算法步驟中，參數(shù)τ按采樣周期離散，參數(shù)a的范圍和步長根據(jù)特征頻率和細(xì)化程度來選取;稀疏性正則化參數(shù) αHsp∈［0.01，0.5］;光滑性正則化參數(shù)αHsm有 三 種 取 值 方 式［9］:常 數(shù);αHcm(k)=20exp

4 仿真分析和實(shí)例分析

4.1 仿真分析

機(jī)械設(shè)備早期故障常常體現(xiàn)為微弱的沖擊信號，常常被淹沒在噪聲或其他較大的振動中，為了驗(yàn)證本文特征提取方法的有效性，現(xiàn)模擬一帶有噪聲的周期性沖擊衰減信號，信號波形如圖1所示。

圖1 含噪聲沖擊信號Fig.1 The impact noise simulation signal waveform

對該模擬信號采用CWT－mcNMF方法進(jìn)行處理。相應(yīng)參數(shù)為:尺度參數(shù) a∈［1，10］，步長為 0.1，β =0.3(如圖 2)，αHsp=0.09，αHsm，其中α0=0.001，J=2。考慮初始值對NMF算法的影響，重復(fù)運(yùn)行算法多次，獲得兩個特征基向量，通過它們間的加權(quán)組合，能夠構(gòu)成原始的最優(yōu)系數(shù)矩陣。其去均值后的波形如圖3。

圖2 沖擊信號小波熵隨波形參數(shù)的變化圖Fig.2 Wavelet entropy change chart with wavelet shape factor

從圖3可看出，第一個基向量含有包含沖擊成分和幅值較大的噪聲信息，而第二個基向量含有明顯的沖擊成分，說明被淹沒的沖擊特征被較好的提取出來。但是由于NMF屬于矩陣分解問題，其結(jié)果存在幅值不確定性，導(dǎo)致兩基向量的幅值與原始波形幅值相差較大。且由于NMF方法屬于非線性特征提取方法，而在分解過程中并未施加正交性等約束，因此獲得的基向量間具有冗余信息，在圖中反映為在沖擊脈沖較大的點(diǎn)，兩個基向量的波形幅值都較大。但以上因素不影響故障特征的提取。

圖3 CWT-mcNMF提取的微弱沖擊信號Fig.3 The impact signal waveform extracted by CWT-mcNMF

4.2 實(shí)例分析

(1)軸承外環(huán)故障

現(xiàn)選取信噪比較低的軸承故障信號進(jìn)行分析。實(shí)驗(yàn)對象為308型滾動軸承，模擬外環(huán)點(diǎn)蝕故障，采用加速度傳感器采集振動信號，軸承轉(zhuǎn)速為1 600 r/min，采樣頻率為40 k，采樣點(diǎn)數(shù)為4 096，故障特征頻率為82.05 Hz。信號時域波形圖4，由圖4可知，該故障信號信噪比較低，其輕微故障引起的周期性沖擊信號被較大的噪聲淹沒，無法直接從時域觀察到?jīng)_擊周期分量。從低頻區(qū)圖5中也無法識別故障特征頻率。

采用CWT－mcNMF提取波形中的周期性沖擊成分。相應(yīng)參數(shù)為:尺度參數(shù) a∈［1，5］，步長為0.1;β =0.1(如圖 6);αHsp=0.01;αHsm(k)=20exp(－ k/50);J=2。同樣重復(fù)運(yùn)行算法多次，獲得兩個特征基向量，去均值后波形如圖7。

分析圖7可知，第一個向量包含有沖擊成分在內(nèi)的軸承振動信息，第二個向量含有明顯的沖擊分量，且連續(xù)較大沖擊間的時間間隔約為12.6 ms，對應(yīng)特征頻率為79.37 Hz，在此頻率分辨率下，即為軸承故障特征頻率。這表明CWT-mcNMF方法可有效提取微弱的沖擊故障特征。同理，由于存在幅值不確定性，基向量與原始波形幅值相差較大;而第一個向量也存在微弱的沖擊信息，說明NMF分解結(jié)果間具有冗余性。

圖4 軸承外環(huán)故障原始波形Fig.4 The bearing outer ring fault signal waveform

圖5 軸承外環(huán)故障原始波形頻譜低頻區(qū)Fig.5 The low-frequency region of the bearing outer ring fault signal frequency spectrum

圖6 外環(huán)故障信號小波熵隨波形參數(shù)的變化圖Fig.6 Wavelet entropy change chart with wavelet shape factor

圖7 CWT－mcNMF提取的軸承外環(huán)故障沖擊信號Fig.7 The outer ring fault impact waveform extracted by CWT－mcNMF

(2)齒輪箱故障

為了進(jìn)一步驗(yàn)證基于稀疏性和光滑性約束的非負(fù)矩陣分解提取沖擊故障的能力，選擇某二級傳動的故障齒輪箱的輸出軸信號進(jìn)行分析。齒輪箱故障為中間軸的一級嚙合齒輪斷裂。電機(jī)轉(zhuǎn)速約為1 500 r/min，采樣頻率為10 240 Hz，采樣點(diǎn)數(shù)為4 096，從輸入軸到輸出軸間齒數(shù)分別為:26、64、40、85。經(jīng)計(jì)算，輸入軸轉(zhuǎn)頻z1=25 Hz，中間軸轉(zhuǎn)頻z2=10.15 Hz，輸出軸 z3=5 Hz。故障特征頻率為中間軸轉(zhuǎn)頻10.15 Hz。采集的輸出軸加速度振動信號時域波形如圖8。由于輸出軸測點(diǎn)距斷齒故障點(diǎn)較遠(yuǎn)，其沖擊故障完全被淹沒，而從其低頻譜區(qū)圖9中也無法直接找到故障特征頻率。

采用CWT－mcNMF方法分析輸出軸信號。相應(yīng)參數(shù)為:尺度參數(shù) a∈［1，10］，步長為 0.1，β =0.2(如圖10)，αHsp=0.09;αHsm(k)=20exp(－ k/50);J=2，獲得兩個特征基向量去均值后的波形，如圖11所示。

分析圖11可知，在兩個特征基向量中，第一個向量包含該輸出軸的嚙合等振動信息，而第二個向量具有明顯的沖擊分量，說明齒輪箱存在故障，但是仍無法直接判斷其故障類型，對其進(jìn)行FFT變換得到提取的故障波形頻譜圖，如圖12所示，頻譜圖的低頻區(qū)中非常清晰的顯示10 Hz的故障特征頻率，與齒輪箱中間軸一級嚙合的齒輪斷齒的故障吻合。

圖8 齒輪箱斷齒輸出軸原始波形Fig.8 The output shaft fault signal waveform of gearbox tooth breaking

圖9 齒輪箱斷齒輸出軸原始波形低頻譜區(qū)Fig.9 The low-frequency region of theoutput shaftsignal frequency spectrum

從以上分析可知:① 由于非負(fù)矩陣分解實(shí)質(zhì)上與獨(dú)立分量分析相同，屬于矩陣分解問題，因此得到的特征基向量具有幅值不確定性和順序不確定性，但并不影響故障特征的提取;② 由于NMF屬于非線性分解，且在分解過程中未施加正交性等約束，因此特征基向量間有冗余信息;③ 上述方法不僅適用于軸承和齒輪箱的故障信號，同樣適用于其他機(jī)械振動信號。

圖10 齒輪箱輸出軸小波熵隨波形參數(shù)的變化圖Fig.10 Wavelet entropy change chart with wavelet shape factor

圖11 CWT-mcNMF提取的齒輪箱輸出軸故障沖擊信號Fig.11 The tooth breaking fault impact waveform extracted by CWT-mcNMF

圖12 CWT-mcNMF提取的故障波形頻譜的低頻區(qū)Fig.12 The low-frequency region of the gearbox tooth breaking fault signal frequency spectrum extracted by CWT-mcNMF

5 結(jié)論

本文針對小波變換在早期故障特征提取中存在的不足，提出結(jié)合連續(xù)小波變換和多約束非負(fù)矩陣分解的故障特征提取方法。該方法根據(jù)故障特征在系數(shù)矩陣中的表現(xiàn)特點(diǎn)，采用基于光滑性和稀疏性約束的非負(fù)矩陣分解方法對連續(xù)小波變換的最優(yōu)系數(shù)矩陣進(jìn)行非線性特征約簡，獲得信噪比較高的故障特征。仿真信號與實(shí)際的軸承和齒輪箱故障樣本的應(yīng)用表明，該方法能夠獲得基于局部的故障特征，實(shí)現(xiàn)非線性降維，能夠在時域提取較明顯的故障信息，很好的滿足了故障診斷特征提取需求。

另外，非負(fù)矩陣分解方法是一種新興的信號處理方法，在機(jī)械故障診斷的領(lǐng)域還比較少，因此，需要進(jìn)一步研究適用于不同類型故障診斷的非負(fù)矩陣分解方法，和算法參數(shù)的控制，以及該方法與其他方法的結(jié)合等。另外，NMF的初始值對算法的收斂性也有較大影響，因此可結(jié)合初始化方法，以期得到最優(yōu)結(jié)果。

［1］梁 霖，徐光華，侯成剛.基于奇異值分解的連續(xù)小波消噪方法［J］.西安交通大學(xué)學(xué)報(bào)，2004，38(9):904 －908.LIANG Lin，XU Guang-hua，HOU Cheng-gang.Continuous wavelet transform denoising method based on singualr value decomposition［J］.Journal of Xi’an Jiaotong University，2004，38(9):904 －908.

［2］Brenner M J.Nonstationary dynamics data analysis with wavelet-SVD filtering［J］.Mechanical Systems and Signal Processing，2003，17(4):765 －786.

［3］Lee D D，Seung H S.Learning in intelligent embedded systems［EB/OL］.http://hebb.mit.edu/people/seung/papers/lee_es99.pdf.1999.

［4］Kawamoto T，Hotta K，Mishima T，et al.Estimation of single tones from chord sounds using non-negative matrix factorization［J］.Neural Network World，2000，3:429－436.

［5］蔡 蕾，朱永生.基于稀疏性非負(fù)矩陣分解和支持向量機(jī)的時頻圖像識別［J］.自動化學(xué)報(bào)，2009，35(10):1272－1277.CAI Lei， ZHU Yong-sheng. Time-frequency spectra recognition based on sparse non-negative matrix factorization and support vector machine［J］.Acta Automatica Sinica，2009，35(10):1272－1277.

［6］欒美潔.基于非負(fù)矩陣分解的機(jī)械故障信號特征分析的研究［D］.南京，東南大學(xué)，2009.

［7］Mallat S.A theory for multiresolution signal decomposition:the wavelet representation［J］. IEEE Trans Pattern Analysisand Machine Intelligence，1989，11(7):674－693.

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

［9］Cichocki A，Zdunek R，Phan A H，et al.Nonnogetive matrix and tensor factorizations［M］.United Kingdom:Wiley ＆Sons，Ltd，Publications，2009.

［10］Phan A H，Cichocki A.Multi-way nonnegative tensor factorization using fast hierarchical alternating leastsquares algorithm(HALS)［C］.Proc.of The 2008 International Symposium on NonlinearTheory and itsApplications，Budapest，Hungary，2008.

［11］栗茂林，梁 霖，王孫安，等.基于連續(xù)小波系數(shù)非線性流行學(xué)習(xí)的沖擊特征提取方法［J］.振動與沖擊，2012，31(1):106－111.LI Mao-lin，LIANG Lin，WANG Sun-an，et al.Mechanical impact feature extraction method based on nonlinearmanifold learning of continuous wavelet coefficients［J］.Journal of Vibration and Shock，2012，31(1):106 －111.

［12］林 京，屈梁生.基于連續(xù)小波變換的信號檢測技術(shù)與故障診斷［J］.機(jī)械工程學(xué)報(bào)，2000，36(2):95－100.LING Jin，QU Liang-sheng.Feature detection and fault diagnosis based on continuous wavelet transform［J］.Journal of Mechanical Engineering，2000，36(2):95－100.