混沌奇異譜特性研究及在滾動軸承故障診斷中的應用

張淑清 賀 朋 左一格 陳榮飛 張 赟 劉 婉 姜萬錄

燕山大學電氣工程學院，秦皇島，066004

0 引言

機械故障診斷技術對產品質量和安全生產至關重要，其研究的關鍵在于信號特征提取和模式識別［1］。軸承故障診斷領域應用混沌理論的研究大致分為兩個方向。一是利用混沌振子檢測淹沒在強噪聲信號中的微弱信號，通過相變的檢測或是通過特征參數(shù)的檢測識別出信號中的微弱成分［2］。PATEL 等［3］利用基于包絡的 Hilbert變換輔以Duffing振蕩器成功識別出軸承的局部損傷。MING等［4］將小波分析思想和混沌振子檢測方法相結合組成聯(lián)合測量系統(tǒng)，準確判斷出軸承的故障。二是利用混沌動力學方法提取特征量，識別各種故障模式［2］，常見的特征統(tǒng)計量有Lyapunov指數(shù)、關聯(lián)維數(shù)和Kolmogovov熵。改進混沌關聯(lián)維數(shù)［5］的計算方法為軸承故障診斷提供了有力依據(jù)。

根據(jù)Takens定理，當重構相空間的維數(shù)m≥2d+1（d為吸引子的維數(shù)）時，動力系統(tǒng)的幾何結構可以完全分開［2］，所以混沌系統(tǒng)拓撲結構的特征量不應受嵌入維數(shù)的影響。事實上，關聯(lián)維數(shù)、Kolmogorov熵、Lyapunov指數(shù)等在無噪聲的情況下均不受嵌入維數(shù)影響，對噪聲是極度敏感的［2,6］。

本文提出一個新的特征量——混沌奇異譜，并提出一種基于混沌奇異譜特征提取的滾動軸承早期故障識別算法。

1 奇異譜理論及特性分析

假設混沌時間序列x(t)可以重構出m維相空間的吸引子，得到軌道矩陣：

其中，m為嵌入維數(shù)，τ為延時時間，N為相點個數(shù)，每一個相點坐標

設混沌動力系統(tǒng)的m個狀態(tài)變量相互獨立且相互正交，狀態(tài)空間分布在m維嵌入空間向量的主分量方向上，則存在變換矩陣A，使YT=ATXT的協(xié)方差矩陣為對角陣，即

式中，CX、CY分別是X和Y的延時-協(xié)變矩陣；sj為CY特征值，j=1,2,…,m。

將特征值由大到小排列，即s1≥s2≥???≥sm，記

稱 S1≥ S2≥ ???≥ Sm為系統(tǒng)的奇異譜［7?8］。將延時-協(xié)變矩陣CX定義為

其中，c(j)表示延時為j時x(n)的協(xié)方差，計算公式如下：

矩陣CX是一個正定對稱矩陣，特征值非負［8］，即e1≥ e2≥???≥ em＞ 0，這些特征向量Ek稱為經(jīng)驗正交函數(shù)（EOF），第k個主成分（PC）定義為混沌時間序列x(t)在第k個經(jīng)驗正交函數(shù)Ek上的正交投影系數(shù)［8］。

1.1 奇異譜的抗噪聲能力分析

設X(t)是零均值、連續(xù)、長度無限、各態(tài)遍歷的時間序列。定義積分算子如下：

其中，符號:=表示定義符號，ρ(t)∈L2(-τ,τ)，L2(-τ,τ)表示在區(qū)間(-τ,τ)上平方積分有限的函數(shù)的全體。CX(t)為時間序列X(t)在延時t的自相關函數(shù)。

若存在

即

則稱ρ(t)為CX(t)的本征函數(shù)，λ為本征函數(shù)對應的本征值。

假設CX(t)在整個時間軸上連續(xù)，并且屬于L2(-2τ,2τ)。由自相關函數(shù)計算公式可知，CX(t)是關于零點對稱函數(shù)，所以積分算子A是對稱正定算子，即(f,Ag)=(Af,g)。根據(jù)Hilbert?Schmidt理論，存在可數(shù)的 λ，且 λ1≥ λ2≥ ???≥ λk≥ ???≥ 0，使得式（9）成立［9］，此時與λ相關的{ρk}構成正交基：

一般稱{ρk}為經(jīng)驗正交函數(shù)。

將Mercer定理應用到Hilbert?Schmidt算子A，自相關函數(shù)可以展開為［8］

對于任意的t=s，若CX(t-s)有意義，則有

其中，CX(0)表示X(t)的總方差。在信號處理理論中，方差反映的是信號的總功率，由式（12）可知，{λk}反映信號功率的分布情況。

假設觀察序列X′(t)由一個內在部分X(t)和外在部分W(t)組成［8］，即

X(t)反映的是系統(tǒng)的內在屬性，既包括確定性成分，也包括隨機部分（混沌系統(tǒng)的偽隨機過程）；W(t)反映的是由誤差引起的純隨機過程，即白噪聲［8］。

觀察有噪聲的序列X′(t)，由于真實信號與噪聲不相關，由連續(xù)函數(shù)的自相關函數(shù)公式可知，序列X′(t)將會產生一個不連續(xù)的自相關函數(shù)。即在原點處，X′(t)的總方差為X(t)的總方差CX(0)和W(t)總方差之和，即

而在非零處，CX'(τ)=CX(τ)。

理論上白噪聲CW(0)→ ∞，CX′(τ)在零點處無意義，式（12）不成立。對特征值進行積分得

從特征值的積分公式（式（15））可以得出，有限間斷點并不影響積分結果。因此，特征值和特征向量在理想噪聲情況下并不受噪聲影響，且此時特征值之和只由內在部分X(t)的變量決定，顯示了當噪聲是理想噪聲時，信號的奇異譜并不存在噪聲平臺。

1.2 奇異譜穩(wěn)定性分析

考慮有限、均值為零、離散、采樣間隔為τs、序列長度為N的時間序列Xi=X(iτs)，自相關函數(shù)可以表示為

使用下標i代表離散時間t，j代表時延。奇異譜的離散形式可以表示為

其中，k=1,2,???,M。同理，自相關函數(shù)線性表示為

在離散情況下，當i=j時，式（18）可以表示為

即信號特征值之和等于信號的方差。由于信號與噪聲信號不相關，所以有自相關函數(shù)：

CX'(i-j)=E[(Yj-0)(Yi-0)]=E[(Xj+Wj)(Xi+Wi)]=E(XjXi)+E(XjWi)+E(WjXi)+E(WjWi)=

E(XjXi)+E(WjWi)=CX(i-j)+CW(i-j) (20)

又由于白噪聲是獨立同分布的隨機變量，其協(xié)方差矩陣為對角陣，即

由此可得

式中，M為協(xié)方差矩陣維數(shù)，也是嵌入維數(shù)。

分析式（22）可知，當τs固定，N趨于無窮，即嵌入維數(shù)M趨于無窮，此時，σ2w/M →0，λ(X′)k代表真實信號奇異譜估計；當τs固定，N有限，λ(X)k較小時，特征值有一個由σ2wM決定的噪聲平臺。

為了敏銳觀察偏差對奇異譜值的影響，引入對數(shù)函數(shù)，從而產生更加穩(wěn)定的噪聲平臺。奇異譜與特征值的關系式如下：

由式（23）可知，奇異譜值的變化由lg(λi+/M)決定，當λi較大時M的波動對奇異譜值的影響較小，可以忽略M的影響，奇異譜值基本保持穩(wěn)定；當λi接近M時，M的波動對奇異譜值的影響較大，奇異譜值會有所提高；當λi遠小于M時，此時奇異譜值幾乎完全由噪聲決定，因此會產生噪聲平臺。

2 奇異譜對混沌空間向量的幾何分析

考慮空間向量x=（x1,x2,…,xn）的一種線性組合：

式中，wk1為權重系數(shù)，組合成n維向量w1。

如果w1能夠使向量空間所有點的映射y1的方差最大，則稱y1為向量空間集合的第一個主成分。方差同時依賴于w1的權重系數(shù)和方向，并隨著w1范數(shù)的增大無限增大，因此通常約束其范數(shù)為1［10］。上述問題轉化為在滿足||w1||=1的范數(shù)條件下，極大化的數(shù)學期望，即

式中，Cx為x的協(xié)方差矩陣。

范數(shù)用歐幾里得范數(shù)表示：

式（25）的極大化問題的解由協(xié)方差矩陣Cx的特征向量e1,e2,…,en給出，特征向量的排序方式是由特征值λ1,λ2,…,λn按照由大到小的方式排列。極大化解為

同理可以推廣到m個主成分的情況，用ym=wTmx表示第m個主成分。wm是單位權重系數(shù)向量，約束關系是ym與前面所有主成分互不相關：

此時再極大化方差E(y2m)=E((wTmx)2)。換言之，在與Cx的前幾個特征向量正交的子空間中尋找y2m的最大數(shù)學期望。由文獻［11］可得

混沌吸引子具有無窮層次的自相似結構，奇異譜是這種空間結構的一種度量方式，反映空間點集整體的離散程度。主成分越大，表示在該方向的離散程度越大［12］。所以，未歸一化的、未取對數(shù)函數(shù)的奇異譜度量的是線性空間點集離散程度。

3 奇異譜抗噪聲及穩(wěn)定性數(shù)值驗證

針對典型的Lorenz系統(tǒng)，有

其中，σ=16，b=4，r=45.92，此時Lorenz系統(tǒng)為混沌系統(tǒng)。取嵌入維數(shù)m=15，時間延遲τ=6。將10對奇異譜對應的主成分分量放在同一個譜圖中，實驗結果見圖1。

由圖1a可以看出，對于一個動力系統(tǒng)，同一分量、不同時間段的時間序列具有極為相似的奇異譜，說明動力系統(tǒng)相空間重構后序列具有相對穩(wěn)定的奇異譜分布，且動力系統(tǒng)不同，其分布不同。針對上文的Lorenz系統(tǒng)，以x分量時間序列{xn}（n=1,2,???,N，N為時間序列長度）為研究對象。向時間序列{xn}添加10 dB噪聲。將10對奇異譜的主成分分量放在同一譜圖中，結果見圖1b。由圖1b可以看出，對于含噪聲的時間序列，其奇異譜分布基本在同一分布上，上下波動不大，說明奇異譜具有一定的抗噪聲干擾能力，魯棒性較好。

圖1 Lorenz系統(tǒng)x分量的10對奇異譜圖Fig.1 10 pairs of singular spectra ofxin Lorenz

對Lorenz系統(tǒng)中x方向的一條時間序列（點數(shù)為5 000）分別添加10～20 dB噪聲，計算其奇異譜，得到奇異譜圖（圖2）。

圖2 Lorenz系統(tǒng)x分量的不同噪聲含量的奇異譜圖Fig.2 Singular spectra with different noise ofxin Lorenz

由圖2可以看出，對于Lorenz系統(tǒng)的15個主成分，10 dB以上的噪聲不影響前7個分量能量大小。這表明噪聲只影響奇異譜能量較低的分量，而對奇異譜能量較高的分量影響不大。隨著噪聲的增強，后8個主成分能量不斷提高，且趨于水平線。

圖3為不同的嵌入維數(shù)下的白噪聲奇異譜圖，可以看出，噪聲的主分量大小相等，能量在主軸上分布均勻。隨著嵌入維數(shù)的增加，只是主分量降低，而均勻分布的狀況并沒有發(fā)生改變。將噪聲添加到時間序列中，噪聲會混亂主成分中能量較低的分量。較低主成分不具有能量下降的趨勢，而呈現(xiàn)出一種均勻分布。能量越大，擾亂的能量主成分個數(shù)越多。

圖3 10 dB噪聲的不同嵌入維數(shù)下的奇異譜圖Fig.3 Singular spectra with 10 dB noise of different embedded dimension

由圖1可看出，對于一個動力系統(tǒng)，同一分量、不同時間段的時間序列具有極為相似的奇異譜，說明動力系統(tǒng)相空間重構后序列具有相對穩(wěn)定的奇異譜分布，且動力系統(tǒng)不同，其分布不同。

4 故障診斷實例

以美國凱斯西儲大學的滾動軸承驅動端軸承的損傷部位（正常信號、內圈故障、滾珠體故障和外圈故障）為例來驗證本文方法的有效性。

4.1 軸承振動信號的混沌特性識別

每種故障選用10對樣本（樣本點數(shù)為2 000點），對于每條樣本序列，使用微分熵法［13］得到樣本序列的嵌入維數(shù)和延時時間，將它們作為小數(shù)據(jù)法的輸入，估計最大Ly?apunov指數(shù)，見圖4。

圖4 不同損傷部位的最大Lyapunov指數(shù)Fig.4 The largest Lyapunov exponent of different damage

由圖4可以看出，軸承不同損傷部位的最大Lyapunov指數(shù)均大于0，表明軸承振動信號具有混沌特性［14］，可以利用混沌相關理論對其進行分析。

4.2 軸承振動信號的奇異譜

由圖4可知，對于滾動軸承不同損傷部位，難以通過最大Lyapunov指數(shù)的閾值將故障分類，因此，采用本文提出的奇異譜方法，針對軸承不同損傷位置的振動信號計算奇異譜，繪制了4種振動信號的奇異譜圖（圖5）。由圖5看出，同一損傷程度的奇異譜分布較接近，上下波動不大，說明不同損傷程度的軸承振動信號都具有較為穩(wěn)定的奇異譜分布，因此，奇異譜可以作為表征振動信號內在屬性的特征量。

圖5 4種振動信號的奇異譜圖Fig.5 The singular spectra of four kinds of vibration signals

4.3 軸承故障診斷

實驗對各部位損傷數(shù)據(jù)分別取50組數(shù)據(jù)，每組時間序列長度為2 000。首先，對信號進行奇異譜分析，為不同損傷部位的數(shù)據(jù)貼上標簽，正常信號、內圈故障、滾珠體故障、外圈故障的標簽分別設定為1、2、3、4。隨機抽取其中的150組數(shù)據(jù)作為訓練數(shù)據(jù)訓練SVM［15］，SVM核函數(shù)選用徑向基函數(shù)，懲罰因子和徑向基函數(shù)參數(shù)通過網(wǎng)格法的交叉驗證獲得［15］。其余的50組作為驗證數(shù)據(jù)。

將嵌入維數(shù)和延時時間的取值范圍均設為3，4，…，22，圖6所示為不同嵌入維數(shù)和延時時間對不同損傷部位的判別準確度的影響。

圖6 嵌入維數(shù)和延時時間對準確度的影響Fig.6 The influence of embedded dimension and delay time on accuracy

由圖6可以看出，當嵌入維數(shù)大于5，延時時間小于13 s時，絕大部分準確度都大于等于96%。當嵌入維數(shù)小于等于7時，準確度較低，因為嵌入維數(shù)較小，混沌系統(tǒng)的相空間重構的拓撲結構沒有被完全打開；延時時間過長，準確度明顯降低。據(jù)此，嵌入維數(shù)和延時時間在合理范圍內，就可以滿足識別要求的準確度。這說明混沌奇異譜受嵌入維數(shù)和延時時間的影響較小。

圖7所示為嵌入維數(shù)m=15，時間延遲τ=10s時的檢測結果，50組待測數(shù)據(jù)的標簽檢測正確率為100%，說明基于奇異譜和支持向量機的判別方法對軸承不同損傷部位的故障診斷是有效的。

為了進一步說明該故障診斷方法的理論具有實際應用價值，筆者在上海寶鋼實際采集了軸承信號并進行了故障識別診斷。該軸承是由熱軋機中的軸承以及齒輪組成，采樣頻率為3 kHz，故障診斷結果見圖8。由圖8a可以看出，實際信號的奇異譜是穩(wěn)定的，圖8b顯示對實際的數(shù)據(jù)進行訓練分類的準確率為100%。實驗證明了基于混沌奇異譜的滾動軸承故障診斷方法的有效性和實用性。

圖7 不同損傷部位的測試樣本的識別結果Fig.7 The results of the test samples of different damage

圖8 寶鋼數(shù)據(jù)實驗結果Fig.8 Experimental results of the data of Baogang

5 結論

本文在研究混沌奇異譜特性的基礎上提出混沌奇異譜提取信號特征量的新方法，通過仿真實例證明了其穩(wěn)定性和抗噪聲能力，并結合SVM對實測軸承故障信號進行分析，驗證了所提方法的有效性。

［1］ 秦大力，于德介.基于本體的機械故障診斷貝葉斯網(wǎng)絡［J］.中國機械工程，2013，24（9）：1195?1209.QIN Dali，YU Dejie.Ontology ? based Diagnostic Bayesian Networks for Mechanical Fault Diagnosis［J］.China Mechanical Engineering，2013，24（9）：1195?1209.

［2］ 吳參，李興林，孫守遷，等.混沌理論在滾動軸承故障診斷中的應用［J］.軸承，2013（1）：60?64.WU Can，LI Xinglin，SUN Shouqian，et al.Applica?tion of Chaos Theory in Fault Diagnosis of Rolling Bearings［J］.Rolling，2013（1）：60?64.

［3］ PATEL V N，TANDON N，PANDEY R K.Defect Detection in Deep Groove Ball Bearing in Presence of External Vibration Using Envelope Analysis and Duffing Oscillator［J］.Measurement，2012，45（5）：960?970.

［4］ MING A B，ZHANG W，QIN Z Y，et al.Envelope Cal?culation of the Multi?component Signal and Its Appli?cation to the Deterministic Component Cancellation in Bearing Fault Diagnosis［J］.Mechanical Systems and Signal Processing，2015，50：70?100.

［5］ SENROY N，SURYANA RAYANAN S，RIBEIRO P F.An Improved Hilbert–Huang Method for Analy?sis of Time?varying Waveforms in Power Quality［J］.IEEE Transactions on Power Systems，2007，22（4）：1843?1850.

［6］ 黃國如，芮孝芳.流域降雨徑流時間序列的混沌識別及其預測研究進展［J］.水科學進展，2004（2）：255?260.HUANG Guoru，RUI Xiaofang.Study Advances in Diagnosis of Chaotic Behavior and Its Prediction for Rainfall and Stream Flow Time Series in Watershed［J］.Advances in Water Science，2004（2）：255?260.

［7］ SAFARI N，CHUNG C Y，PRICE G C D.A Novel Multi?stepShort?termWindPowerPredictionFrame?work Based on Chaotic Time Series Analysis and Singular Spectrum Analysis［J］.IEEE Transactions on Power Systems，2018，33（1）:590?601.

［8］ 劉元峰，趙玫.基于奇異譜分析的混沌序列降噪［J］.上海交通大學學報，2003，37（5）：778?780.LIU Yuanfeng，ZHAO Mei.De?noising of Chaotic Time Series Based on Singular Spectrum Analysis［J］.Journal of Shanghai Jiaotong University，2003，37（5）：778?780.

［9］ ELSNER J B，TSONIS A A.Singular Spectrum Analysis：a New Tool in Time Series Analysis［M］.Berlin：Springer Science&Business Media，2013：932?942.

［10］ LI Y，ZHANG C，WANG S.Low?complexity Non?uniform Penalized Affine Projection Algorithm for Sparse System Identification［J］.Circuits，Systems，and Signal Processing，2016，35（5）：1611?1624.

［11］ 楊竹青，李勇，胡德文.獨立成分分析方法綜述［J］.自動化學報，2002，28（5）：762?772.YANG Zhuqing，LI Yong，HU Dewen.Independent Component Analysis：a Survey［J］.Acta Automati?ca Sinica，2002，28（5）：762?772.

［12］ RAMOS G L，RODRIGUES G F，REGO C G.FDTD Near Field Data Reduction Using Principal Component Analysis［C］//Microwave and Opto?electronics Conference.Ipojuca，2015：1?4.

［13］ 張淑清，李新新，張立國，等.基于符號分析的極大聯(lián)合熵延遲時間求取方法［J］.物理學報,2013,62（11）：11056.ZHANG Shuqing，LI Xinxin，ZHANG Liguo，et al.Delay Time Obtain Method Using the Maximum Joint Entropy on the Basis of Symbolic Analysis［J］.Acta Physica Sinica，2013,62（11）：11056.

［14］ 劉劍鳴，楊霞，高躍龍，等.類Liu系統(tǒng)在水聲微弱信號檢測中的應用研究［J］.物理學報，2016，65（7）：70501.LIU Jianming，YANG Xia，GAO Yuelong，et al.Application of Similar Liu System in Underwater Weak Signal Detection［J］.Acta Physica Sinica，2016，65（7）：70501.

［15］ 杜紅樂，張燕.不均衡數(shù)據(jù)混合取樣分類算法［J］.燕山大學學報，2015（2）：158?164.DU Hongle，ZHANG Yan.A Classification Algo?rithm Based on Mixed Sampling for Imbalanced Dataset［J］.Journal of Yanshan University，2015（2）：158?164.*