基于連續(xù)小波系數非線性流形學習的沖擊特征提取方法

栗茂林，梁 霖，王孫安，莊 健

(西安交通大學 機械工程學院，西安 710049)

基于連續(xù)小波系數非線性流形學習的沖擊特征提取方法

栗茂林，梁 霖，王孫安，莊 健

(西安交通大學 機械工程學院，西安 710049)

為了提取機械設備故障引發(fā)的沖擊成分，提出了一種基于連續(xù)小波系數非線性流形學習的沖擊故障特征提取方法。首先，基于小波熵方法優(yōu)化出最優(yōu)的Morlet小波波形參數，實現(xiàn)與沖擊特征成分的最佳匹配，獲取包含沖擊特征信息的最優(yōu)小波系數矩陣。其次，采用局部切空間排列算法對最優(yōu)小波系數矩陣進行非線性約簡，并基于峭度指標最大化原則，確定出特征空間中的有效低維嵌入，從而提取出最優(yōu)的沖擊故障特征。最后，通過仿真數據和工程實際的應用對比分析，表明該方法采用了局部線性化和全局排列的思想，與線性奇異值分解方法相比，不僅在時域上提取出峭度更大的微弱沖擊特征成分，而且在頻譜中還提取出了相應的低頻故障特征。

特征提取;連續(xù)小波變換;非線性流形學習;沖擊故障

在機械設備的故障中，沖擊類故障占有較高的比例，如動靜件間的周期性碰撞、缺陷軸承或齒輪運行時的瞬時沖擊等。沖擊響應中包含了周期性的反映固有頻率的有阻尼簡諧振動，其中，沖擊力大小、固有特征及沖擊重復周期等都是診斷的重要信息。然而由于設備的早期故障，往往使得這些周期沖擊簡諧振動特征被噪聲淹沒而不易識別。小波變換作為一種時頻分析方法是研究非平穩(wěn)沖擊故障特征的熱門工具。

在小波變換中，連續(xù)小波變換的尺度連續(xù)變化可以充分細致地刻劃信號的局部形態(tài)。采用與沖擊特征波形近似的基小波進行連續(xù)小波變換時，沖擊特征分量的變換系數與噪聲等其他分量不同，因此通過時頻分布圖可以觀察到變換系數的周期性變化情況，如文獻［1］通過最優(yōu)Morlet小波變換在時間-尺度分布中反映出軸承和齒輪振動信號中的周期分量;文獻［2］通過Haar小波變換的灰度圖檢測出變速箱周期性的沖擊特征。但另一方面，時頻分布圖雖然可以直觀形象地反映出各尺度下的小波系數變化情況，但這種表示方法不利于后續(xù)的診斷識別處理。因此，文獻［3］對變換系數進行處理，通過統(tǒng)計尺度-小波能量譜的分布特性曲線提取出更好的信號特征;文獻［4］認為瞬態(tài)成分和噪聲的小波變換系數具有不同的統(tǒng)計分布特征，因此通過統(tǒng)計檢驗方法提取出小波系數中的瞬態(tài)成分。這些研究成果直接對部分尺度的小波系數分析處理，對人員依賴性強，需要先驗知識，容易導致診斷信息的丟失。

研究表明，沖擊特征信息主要被分解在小波變換系數的部分尺度中，使小波變換系數矩陣具有一定的稀疏性。因此，對矩陣進行稀疏性分解的奇異值分解(singular value decomposition，SVD)方法被應用于連續(xù)小波系數的特征提取中［5-6］，其思路是對小波系數矩陣用線性方法來逼近系統(tǒng)的空間流形，通過分解的主要特征向量得到系數矩陣中的固有成分，即沖擊特征。這樣既體現(xiàn)了連續(xù)小波變換的高分辨率分析的優(yōu)點，又克服了連續(xù)小波變換各尺度系數間的冗余性。但是，研究表明奇異值分解本身屬于線性學習方法，忽略了系數矩陣數據集存在的凸凹性，僅能提取出其中的線性特征，不能找到真正的數據分布結構。

而非線性流形學習方法為基于數據分布的內在結構分析提供了一種新的研究途徑，其基本思想是:高維觀測空間中的點由少數獨立變量的共同作用在觀測空間張成一個流形，在盡可能地保證數據間的幾何關系和距離測度不變的前提下，有效地展開觀測空間卷曲的流形來發(fā)現(xiàn)內在的主要變量，實現(xiàn)對數據集的約簡。這就意味著流形學習比傳統(tǒng)線性約簡方法更能體現(xiàn)事物的本質，更有利于對數據的理解和進一步處理。如黎敏等人對高維相空間中重構的流形拓撲結構進行綜合信息提取［7］;蔣全勝等人運用拉普拉斯特征映射的非線性降維方法保留了振動信號中內含的整體幾何結構信息［8］，這些研究為非線性流形學習在故障診斷中的應用開辟了思路。

本文針對早期故障診斷中微弱沖擊特征提取的問題，將適用于故障特征流形空間的局部切空間排列算法(local tangent space alignment，LTSA)應用于連續(xù)小波變換(continuous wavelet transform，CWT)最優(yōu)系數矩陣的稀疏性分析中，提出基于連續(xù)小波系數非線性流形學習的沖擊特征提取方法CWT-LTSA。其中，基于峭度指標最大化原則確定有效的低維嵌入向量，提取最優(yōu)的沖擊故障特征。通過仿真數據和實際案例的應用對比分析，表明方法的有效性。

1 連續(xù)小波變換CWT

對于能量有限的信號x(t)，其連續(xù)小波變換可以表示為與小波函數的內積［9］，即:

式中:τ、a∈L2;τ為時間位置參數;a為尺度參數;ψ(t)為小波函數;ψτ，a(t)為基小波函數。

對于存在傅里葉變換對的基小波函數，由傅里葉變換的卷積性質，式(1)又可表示為:

式中:F-1表示傅里葉逆變換，X(f)和ψ*(f)是x(t)和ψ(t)的傅里葉變換。

式(2)變換后的小波系數矩陣W(τ，a)通過時間和尺度的變化度量了信號x(t)與小波ψ(t)之間的相似程度。信號中的特征信息(如頻率和周期重復性)反映在系數矩陣的變化中。

2 局部切空間排列算法LTSA

局部切空間排列算法LTSA通過逼近每個樣本的切空間來構建低維流形的局部幾何，并利用局部切空間排列求出整體低維嵌入坐標，恢復出流形等距的低維空間子集［10］。算法包括了兩個過程。

2.1 局部線性擬合

2.2 局部坐標全局排列

為了唯一地確定T，引入約束TTT=Id。

由于全1向量e是矩陣B的零特征值對應的特征向量，所以取矩陣B的第2到第(d+1)個最小的特征值對應的特征向量所組成的矩陣就是所求的T，即為描述了樣本空間矩陣X中非線性主流形的正交低維全局坐標映射矩陣。

3 基于流形的最優(yōu)小波系數特征提取

基于全局思想的流形學習方法，如ISOMAP，要求流形所對應低維空間的子集是凸的，才能保證構造的流形上距離遠的點之間的聯(lián)系是準確的。然而，基于局部思想的LTSA方法只考慮流形鄰近點之間的關系，不要求流形所對應低維空間的子集是凸的，有著更廣泛的適用對象。因此，本文結合沖擊特征的特點，提出了一種連續(xù)小波系數局部切空間流形學習的沖擊特征提取方法CWT-LTSA。通過對最優(yōu)連續(xù)小波系數矩陣流形的學習，結合峭度指標最大化原則，提取出時頻分布中的低維嵌入，實現(xiàn)沖擊特征的提取。具體的步驟如下:

步驟1基于Morlet小波的樣本空間獲取

由于Morlet小波與機械振動中的瞬態(tài)沖擊特征比較近似，可以保證與信號較好的匹配性，因此是提取沖擊特征的有效小波函數，其數學表達如下:

對于機械振動信號來說，通常采用其實部作為基小波，即:

式中:β為波形參數，控制著小波波形的衰減率。為了與沖擊信號取得最優(yōu)匹配，采用小波熵［1］確定最優(yōu)波形參數 βop，獲取最優(yōu)系數矩陣 Wβop(τ，a)。

式中:j為波形參數β的變化，通常在［0.1，2］范圍內取值，S(Wj)為小波系數矩陣的信息熵。

步驟2樣本空間的局部鄰域選取

研究表明，不同的鄰域點k將獲得不同的嵌入效果:鄰域太小使得切空間易非滿秩，而過大又容易損失局部幾何信息［10］，本文依據經驗取值。

步驟3局部線性投影

對每個樣本向量點的鄰域，計算中心化矩陣XixieT的d個最大奇異值對應的左奇異向量，并將這d個左奇異向量組成矩陣Qi。

步驟4低維特征的提取

為了提取最優(yōu)的沖擊特征，以反映沖擊特征較好的峭度指標作為衡量指標，即最優(yōu)的低維嵌入波形uop將具有最大的峭度值。

式中:Kurtosis為峭度指標。

在這里，需要指出的是，由于尺度的連續(xù)變化，沖擊特征分量的小波變換系數在尺度范圍內將具有一定的相似性，因此，提取的低維全局坐標矩陣T=［u2，u3，…，ud+1］T中的某些維度上會出現(xiàn)類似文獻［11］中提到的坍塌現(xiàn)象，即由于相似度較高，數據在約簡后的特征空間中表現(xiàn)為前幾維的坐標值具有近似性，但由于沖擊分量幅值、初相有差異，并不影響其正交性。另外，由式(10)引入的約束TTT=Id可知，低維全局坐標矩陣T是歸一化的特征向量，提取的波形顯示歸一化幅值，采用與幅值大小無關的峭度指標有效地反映歸一化波形的沖擊衰減程度。

4 仿真試驗與工程應用分析

4.1 仿真試驗分析

為驗證方法的有效性，模擬一周期的沖擊衰減信號，并包含一定的噪聲成分，其波形如圖1所示。對模擬數據采用CWT-LTSA方法進行計算，其中最優(yōu)波形參數 βop為0.3，尺度參數 a∈［1，20］，步長為0.1，鄰域點k依經驗取為10，嵌入維數d取為5。

圖1 加噪的沖擊模擬信號波形Fig.1 The impact noise simulation signal waveform

經計算，提取低維全局坐標下的波形峭度變化情況如圖2所示，最大的峭度值為15.92，對應最優(yōu)嵌入維數dop=1，提取的波形如圖3所示。為了對比學習效果，采用CWT-SVD方法提取的沖擊信號波形如圖4所示，波形的峭度指標值僅為7.56。對比可知，在時域波形中，CWT-SVD方法也提取出了信號的沖擊成分，但CWT-LTSA的效果更好，具有更小的噪聲干擾，且沖擊特征的周期性更好。

圖2 低維嵌入空間下的峭度值Fig.2 The kurtosis index in low-dimensional embedding space

圖3 CWT-LTSA提取的沖擊信號波形Fig.3 The impact signal waveform extracted by CWT-LTSA

圖4 CWT-SVD提取的沖擊信號波形Fig.4 The impact signal waveform extracted by CWT-SVD

分析可知，SVD分解的是兩兩正交的線性子空間，而LTSA提取的全局坐標是非線性子空間。由于SVD將小波系數矩陣作為整體線性化考慮，使得噪聲信息在各個子空間中均有一定的投影，沖擊分量的投影也因丟掉一些細節(jié)信息受到一定程度的影響。相比而言，LTSA算法通過局部鄰域的線性化處理以及全局排列的策略獲取非線性子空間，非線性坐標更加符合由沖擊分量形成的連續(xù)變化的局部幾何結構，而噪聲分量由于隨機波動形成變化不連續(xù)的局部結構，通過局部線性化處理平滑掉大部分噪聲分量。因此，沖擊分量在非線性子空間中的投影幅值大于線性子空間的幅值，且相當一部分的噪聲信息被分解在其他子空間中，從而使得沖擊分量所在的子空間中具有較小的噪聲干擾。

4.2 滾動軸承故障特征的提取分析

現(xiàn)選取信噪比較低的軸承故障測試信號進行應用分析。由加速度傳感器在軸承座上拾取308軸承的輕微外環(huán)故障得到的信號如圖5所示，軸承轉速為1 600 r/min，采樣頻率為 40 k，故障特征頻率為 82.05 Hz。波形中由故障引起的周期沖擊衰減成分受到較大的噪聲干擾而被淹沒，無法觀察到周期分量。圖6所示的對應頻譜中也無法體現(xiàn)故障特征頻率。

首先，基于Morlet小波變換對軸承故障振動信號進行分析處理，尺度參數 a∈［1，20］，步長為 0.1，最優(yōu)波形參數βop為0.2。在最優(yōu)包絡的小波變化系數中，根據系數能量最大的尺度5.5對小波系數矩陣進行切片提取的波形如圖7所示。從圖中可以看出，小波系數中包含的特征成分及周期性變化不明顯，需要進一步的分析處理。

其次，采取CWT-LTSA來提取小波系數矩陣中的沖擊特征，獲取的低維全局坐標下的波形峭度變化情況如圖8所示，選取峭度指標值19.11對應的最優(yōu)低維嵌入dop=3作為提取的波形，如圖9所示。圖10所示為其他3個低維嵌入(d=1，d=2和d=4)的波形。與圖9所示的最優(yōu)波形具有一定的相似性，即提取出一系列的沖擊衰減分量，但提取質量較差。另外，采用CWT-SVD方法進行對比分析，其提取的沖擊波形如圖11所示，峭度指標值為12.20。由此可見，僅從提取波形的沖擊成分來看，顯然非線性流形學習的特征提取能力更強，提取的沖擊成分更明顯，噪聲干擾更小。

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

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

圖7 CWT系數的切片波形Fig.7 The slice waveform of CWT coefficients

圖8 不同嵌入維數下的峭度值Fig.8 The kurtosis index in embedding space

圖9 CWT-LTSA提取的外環(huán)故障沖擊波形Fig.9 The outer ring fault impact waveform extracted by CWT-LTSA

圖10 其他嵌入維數中的外環(huán)故障沖擊波形Fig.10 The outer ring fault impact waveform in the other embedding space

圖11 CWT-SVD提取的外環(huán)故障沖擊波形Fig.11 The outer ring fault impact waveform extracted by CWT-SVD

另外，從頻域角度來看，圖12所示為非線性流形學習和奇異值分解提取波形的頻譜同原始波形頻譜的對比圖。由頻率組成可見，將小波系數矩陣作為整體線性化考慮的CWT-SVD方法，將信息分解至不同的子空間中。由于各子空間包含的信息成分不同，重構提取過程相當于帶通濾波，頻率成分集中在2 kHz-5 kHz范圍，沒有提取出82.05 Hz低頻故障特征頻率。而CWT-LTSA方法通過局部鄰域的線性化處理以及全局排列后的特征提取策略，提取的低維嵌入波形的頻譜中不僅包含沖擊特征對應的主要頻帶，還提取出系數矩陣中局部鄰域間的變化信息，即小波系數值的周期性波動，使得如圖13所示的頻譜低頻區(qū)存在82.05 Hz的外環(huán)故障特征頻率，更好地滿足特征提取的需要。

圖12 兩種方法提取波形的頻譜與原始波形頻譜的對比Fig.12 The frequency spectrum comparison with the extracted signals by two methods and the original waveform

圖13 CWT-LTSA提取外環(huán)故障波形頻譜的低頻區(qū)Fig.13 The low-frequency region of the bearing outer ring fault signal frequency spectrum extracted by CWT-LTSA

為了進一步驗證非線性流形學習提取沖擊特征的能力，選擇308軸承的內環(huán)故障數據進行分析。轉速為1 600 r/min，采樣頻率為40 k，原始波形和頻譜分別如圖14和15所示，其中，原始波形的峭度值為9.38，故障的特征頻率為131.3 Hz。在原始波形的頻譜中，故障特征成分被107.4 Hz和156.2 Hz的兩個分量干擾無法識別。

圖14 內環(huán)故障原始波形Fig.14 The bearing inner ring fault signal waveform

圖15 內環(huán)故障原始波形頻譜的低頻區(qū)Fig.15 The low-frequency region of the bearing inner ring fault signal frequency spectrum

采用CWT-LTSA方法提取的沖擊波形如圖16所示，峭度值為14.58(對應第2維低維嵌入)，噪聲成分明顯減少，頻譜如圖17所示，131.3 Hz的故障特征頻率被有效地提取出來。CWT-SVD提取的效果如圖18所示，故障特征頻率沒有被有效分離出來，由此可見本文所提方法的有效性。

圖16 CWT-LTSA提取的內環(huán)故障沖擊波形Fig.16 The inner ring fault impact waveform extracted by CWT-LTSA

圖17 CWT-LTSA提取內環(huán)故障波形頻譜的低頻區(qū)Fig.17 The low-frequency region of the bearing inner ring fault signal frequency spectrum extracted by CWT-LTSA

圖18 CWT-SVD提取的內環(huán)故障波形頻譜的低頻區(qū)Fig.18 The low-frequency region of the bearing inner ring fault signal frequency spectrum extracted by CWT-SVD

5 結論

本文針對連續(xù)小波變換在沖擊特征提取應用中存在的不足，采樣局部切空間排列算法，提出了一種基于連續(xù)小波系數非線性流形學習的沖擊故障特征提取方法。該方法利用局部切空間排列算法，對連續(xù)小波變換的最優(yōu)高維系數矩陣進行非線性變換，基于峭度指標最大化原則，在特征空間中提取出反映沖擊故障的低維嵌入。仿真數據和滾動軸承內外環(huán)故障樣本的應用表明，與線性奇異值分解方法相比，由于非線性切空間排列流形學習方法考慮了樣本空間中數據集的局部變化信息，因此，不僅在時域上提取出更完整的微弱沖擊特征成分，而且在頻譜中提取出了低頻故障特征分量，更好地滿足了診斷特征的提取需求。

另外，局部切空間排列算法的提取效果較奇異值分解算法好，但在一定程度上受到鄰域點等參數選取和樣本均勻性的影響，因此，后續(xù)將重點研究鄰域點的自適應選取及樣本對嵌入結果的影響。

［1］林 京，屈梁生.基于連續(xù)小波變換的信號檢測技術與故障診斷［J］.機械工程學報，2000，36(2):95-100.

［2］張梅軍，唐 建，陳江海.基于連續(xù)小波灰度圖的變速箱故障診斷［J］.振動、測試與診斷，2007，27(1):65-66.

［3］劉 剛，屈梁生.應用連續(xù)小波變換提取機械故障的特征［J］.西安交通大學學報，2000，34(11):74-79.

［4］朱忠奎，顧 軍，芮延年，等.基于連續(xù)小波和統(tǒng)計檢驗的瞬態(tài)成分檢測與應用［J］.振動工程學報，2006，19(4):559-565.

［5］梁 霖，徐光華，侯成剛.基于奇異值分解的連續(xù)小波消噪方法［J］.西安交通大學學報，2004，38(9):904-908.

［6］Brenner M J.Nonstationary dynamics data analysis with wavelet-SVD filtering［J］.Mechanical systems and signal processing，2003，17(4):765-786.

［7］黎 敏，徐金梧，陽建宏，等.一種基于流形拓撲結構的軸承故障分類方法［J］.控制工程，2009，16(3):358-362.

［8］蔣全勝，賈民平，胡建中，等.基于拉普拉斯特征映射的故障模式識別方法［J］.系統(tǒng)仿真學報，2008，20(20):5710-5713.

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

［10］ Zhang Z Y，Zha H Y.Principal manifolds and nonlinear dimensionality reduction via tangent space alignment［J］.SIAM JournalofScientific Computing，2004，26(1):313-338.

［11］詹德川，周志華.基于流形學習的多示例回歸算法［J］.計算機學報，2006，29(11):1948-1955.

Mechanical impact feature extraction method based on nonlinear manifold learning of continuous wavelet coefficients

LI Mao-lin，LIANG Lin，WANG Sun-an，ZHUANG Jian

(School of Mechanical Engineering，Xi'an Jiaotong University，Xi'an 710049，China)

To acquire an impact component aroused by mechanical fault，a novel feature extraction method based on nonlinear manifold learning of continuous wavelet coefficients was put forward.Firstly，the wavelet entropy method was adopted to optimize the Morlet wavelet shape factor in order to match with the impact components to obtain the optimal continuous wavelet coefficients.Secondly，the nonlinear manifold learning algorithm named local tangent space alignment was used to reduce the optimal wavelet coefficients matrix，and according to the principle of the maximum kurtosis index，the low-dimensional embedded vectors introduced to reflect impact failures were extracted from the global coordinate feature matrix. Finally， simulations and industrial applications showed that compared with the singular value decomposition，this approach is effective to extract not only the weak impacts with the greater kurtosis in time waveform，but also the fault feature frequencies in frequency spectrum.

feature extraction;continuous wavelet transformation;nonlinear manifold learning;impact fault

TH165.3

A

國家自然科學基金項目(51075323，50705073);中央高?；究蒲袠I(yè)務費專項資金資助(xjj20100066);北京交通大學軌道車輛結構可靠性與運用檢測技術教育部工程研究中心開放課題(SROM RGV(BJTU)2010-002)

2010-06-17 修改稿收到日期:2010-12-13

栗茂林 女，博士生，講師，1978年12月生