基于HLS-SVDR和SPPCS的CEEMD的滾動軸承微故障特征提取*

徐 波， 周鳳星， 馬婭婕， 嚴?？?， 黎會鵬

(1.武漢科技大學信息科學與工程學院 武漢，430081) (2.黃岡師范學院電子信息學院 黃岡，438000)

引 言

Huang等[1]提出了Hilbert-Huang變換算法(Hilbert-Huang transform, 簡稱HHT)，適合處理非平穩(wěn)隨機信號，在機械故障診斷中得到廣泛應用。經驗模態(tài)分解(empirical mode decomposition, 簡稱EMD)作為HHT的核心部分，算法直觀簡單，具有正交性、完備性和自適應性等特點，然而存在包絡擬合過沖/欠沖、端點效應和模態(tài)混疊等問題，嚴重影響其推廣和應用。針對模態(tài)混疊問題，胡愛軍等[2]提出了加入高頻諧波再進行EMD分解的模態(tài)混疊的消除方法。曹瑩等[3]提出了基于形態(tài)濾波預處理與特征尺度匹配的模態(tài)混疊抑制方法。但上述方法經驗性較強，缺乏自適應性。湯寶平等[4]提出了基于獨立分量分析的模態(tài)混疊消除方法，但該方法沒有考慮端點效應嚴重影響EDM分解結果的有效性。Yeh等[5]提出了互補集總經驗模態(tài)分解方法，通過向原始信號中添加兩對相反的白噪聲信號分別進行EMD分解，克服了EEMD的白噪聲殘留的問題，有效消除了模態(tài)混疊問題，對端點效應也有一定的抑制作用，是目前消除模態(tài)混疊簡單且有效的方法。

針對包絡擬合過沖/欠沖問題，相關學者提出了如三次B樣條插值[6]、分段三次Hermite樣條插值[7]、分段拋物線插值[8]、單調分段三次樣條插值法[9]和分段約束三次樣條插值[10]等改進方法。上述方法中，三次B樣條插值的擬合特性沒有明顯提高；分段三次Hermite樣條插值的平滑性不好；單調分段三次樣條插值法的計算量太大；分段約束三次樣條插值的包裹性差。

針對端點效應問題，相關學者提出了如鏡像延拓[11]、最大相關波形延拓[12]、極值波延拓[13]、最大斜率再優(yōu)化極值延拓[14]、PDE信號修補[15]、均生函數(shù)周期疊加外推[16]、最小平方距離相關[17]及加窗函數(shù)[18]等基于信號特征或極值特征的延拓方法。但這些方法具有很強的經驗性，泛化能力差。另外一些學者提出了如支持向量回歸預測[19]、支持向量機[20]、神經網絡延拓[21]及自回歸綜合滑動平均模型延拓[22]等基于數(shù)據預測延拓的方法，能夠較為準確地進行端點延拓，然而算法復雜、效率低、參數(shù)難以選擇，工程實現(xiàn)度低。

筆者在總結上述改進方法的基礎上，提出同倫-最小二乘支持向量雙回歸的保形分段三次樣條CEEMD方法，獲得一個完備的CEEMD算法。通過對實際振動信號的分析及滾動軸承微故障頻率特征的提取來驗證該方法能夠從整體上解決EMD存在的擬合過沖/欠沖、端點效應和模態(tài)混疊問題。

1 基本原理

1.1 CEEMD簡介

在對信號x(t)做EEMD分解時，令集總平均的次數(shù)為N，對x(t)第j次添加白噪聲uj(t)后得到

(1)

令x(t)減去uj(t)，得到

(2)

(3)

對式(3)求集總平均，獲得

(4)

由于對x(t)分別添加同一白噪聲，但符號相反，使得相加后殘留在IMFj+和IMFj-中的白噪聲基本互相抵消，從而使最后IMF分量中的白噪聲基本互相抵消，上述添加符號相反的白噪聲正是“互補”的本質。CEEMD不僅消除了模態(tài)混疊問題，而且基本上消除了白噪聲殘留問題。

1.2 包絡擬合過沖/欠沖分析

Huang在文章中沒有給出構造上、下包絡線時產生包絡擬合過沖/欠沖問題的嚴格定義。龍思勝等[23]認為“擬合過沖/欠沖”現(xiàn)象是在待處理數(shù)據序列的每個子段內，所得擬合函數(shù)的部分數(shù)值大于這兩個數(shù)據的大者或小于這兩個數(shù)據的小者，并提出了一階導數(shù)零點法，一定程度上解決了擬合過沖/欠沖問題。由于沒有嚴格的數(shù)學依據及證明，而是根據處理數(shù)據時的經驗，對于任意信號的有效性有待驗證。使用最常用的三次樣條插值對取自滾動軸承實驗臺的故障信號x(t)(如圖1所示，其中x(t)為時間數(shù)據序列為無量綱單位)進行一次上、下包絡線插值擬合，得到包絡線如圖2所示。

圖1 滾動軸承振動信號x(t)Fig.1 Rolling bearing′s vibration signal x(t)

圖2 三次樣條插值的包絡線Fig.2 The envelope curve of spline

為了便于觀察，將圖2中包絡線進行局部放大，其放大圖如圖3，4所示，可以觀測到包絡線出現(xiàn)了明顯的擬合過沖/欠沖問題。

1.3 端點效應分析

由于實際采集的信號的長度有限，信號的兩端在極少數(shù)情況下能夠獲得有效的極值點，因此對極值數(shù)據序列進行插值擬合時，在上、下包絡線在的兩端出現(xiàn)發(fā)散現(xiàn)象，稱為端點效應，且這種“發(fā)散”會隨著分解過程地繼續(xù)進行而不斷向信號內部傳播，最終“污染”整個信號而使分解結果嚴重失真而失去物理意義。采用圖1中滾動軸承的振動信號x(t)進行端點效應的分析，并使用三次樣條插值法對x(t)進行包絡線擬合，結果如圖5所示。由于信號x(t)兩端不存在極值點，獲得的上、下包絡線在曲線兩端出現(xiàn)嚴重的端點效應。

圖3 三次樣條插值的過沖問題的局部放大圖Fig.3 The enlarged partial view of the overshoot problem of spline

圖4 三次樣條插值欠沖問題的局部放大圖Fig.4 The enlarged partial view of the undershoot problem of spline

圖5 三次樣條插值端點效應的局部放大圖Fig.5 The enlarged partial view of the left end effect of spline

2 保形分段三次樣條插值算法

2.1 保形分段插值的構造

給定平面點集Ti(xi,yi)(i=1,2,…,n)，將相鄰兩點使用直線依次連接成一個多邊形{T0,T1,…,Tn}，并令多邊形的邊向量為

Si=ti-ti-1(i=1,2,…,n)

(5)

如果ti×ti-1與ti×ti+1平行且方向相反，則Ti為{T0,T1,…,Tn}的拐點。

2.2 插值曲線在型值點處的切矢

由多邊形的邊矢定義可知，插值曲線在Ti處的切矢Ki定義為：如果Ti是多邊形{T0,T1,…,Tn}的拐點，則令

Ki=kiSi+(1+ki)(-Si+1)

(6)

否則令

Ki=kiSi+(1-ki)(Si+1)

(7)

其中：Ki(0

如果Ti-1，Ti，Ti+1三點共線，則定義Ti-1=Ti=Ti+1=Si，即選取ti-1=0，ti+1=1。通過切矢的定義可知，KiKi+1總是平行KiSi且有相同方向，故分別過Ti和Ti+1且平行于Ki和Ki+1直線li和li+1的交點Vi為

(i=0,1,2,…,n-1)

(8)

2.3 三次樣條插值的保形性

如果Vi-1，Vi，Vi+1，Vi+2構成的三邊形是凸的，則令控制多邊形{v-1v0…vnvn+1}的邊矢量為

ai=Vi-Vi-1(i=0,1,…,n+1)

(9)

給定m+1個控制頂點ri(i=1,2,…,m)，相應的三次樣條曲線的分段函數(shù)表示為

(0≤u≤1;k=0,1,2,…,m-3)

(10)

其中:基函數(shù)為

(11)

從基函數(shù)可以得到端點的性質

(12)

則三次樣條曲線的deBoor點定義為

(13)

由式(9)及式(13)可得

a0=(a1-a2)/4

(14)

a0a1=a1a2/4

(15)

Vi-1，Vi，Vi+1，Vi+2構成的三邊形是凸的，可保證插值曲線在兩端不會出現(xiàn)多余的拐點，并設定λi的取值范圍，使得三次樣條插值曲線保形。如果aiai+1，ai+1ai+2方向一致，則Vi-1，Vi，Vi+1，Vi+2構成的三邊形是凸的，所得交點為

(16)

如果aiai+1，ai+1ai+2方向相反時，則Vi-1，Vi，Vi+1，Vi+2構成另一個拐向，記為

(17)

令

(i=0,1,…,n)

(18)

根據圖6所示的控制多邊形，如果選取

0<λi≤ei(i=0,1,…,n)

(19)

則通過調節(jié)參數(shù)λi滿足式(19)，構造出保形分段三次樣條插值法。

圖6 控制多邊形Fig.6 Controlled polygon

2.4 實驗驗證

為了驗證筆者提出的SPPS插值方法能夠有效地消除包絡擬合過沖/欠沖問題，使用圖1中的振動信號x(t)進行實驗分析，結果如圖7～9所示。使用改進的三次B-樣條插值法進行包絡線插值擬合，擬合結果如圖7所示。對比圖7與圖2的擬合結果可知，改進的三次B-樣條插值得到的上、下包絡曲線依然出現(xiàn)嚴重的擬合過沖/欠沖，并沒有消除擬合過沖/欠沖問題。

為了進一步對比其他插值方法，采用分段約束三次樣條插值方法進行包絡線擬合，擬合結果如圖8所示。對比圖8與圖7的擬合結果可以看到，采用分段約束三次樣條插值方法擬合的上、下包絡曲線沒有出現(xiàn)擬合過沖/欠沖現(xiàn)象，但從圖中A，B兩處擬合的曲線段可以觀察到曲線擬合過陡，曲線不夠平滑，使得擬合曲線的有效性不能得到保證，限于篇幅，此處不對上述其他插值方法進行對比分析。

針對圖8中的擬合曲線出現(xiàn)不夠平滑的問題，筆者提出了保形分段三次樣條插值法擬合包絡曲線，擬合結果如圖9所示。擬合結果表明，該方法不但能夠有效抑制擬合過沖/欠沖問題，而且擬合的包絡曲線與圖8相比，更平滑、有效，為后續(xù)的CEEMD分解過程奠定了基礎。

圖7 三次B-樣條插值的包絡線Fig.7 The envelope of B-spline

圖8 分段約束三次樣條插值的包絡線Fig.8 The envelope of constrained piecewise spline

圖9 保形分段三次樣條插值的包絡線Fig.9 The envelope of shape-preserving piecewise spline

3 同倫-最小二乘支持向量雙回歸算法

支持向量回歸(support vector machine regression，簡稱SVR)算法具有十分出色的預測能力。Suykens等[24]提出了最小二乘支持向量回歸(least squares support vector regression，簡稱LS-SVR)算法，提高了SVR的求解效率，降低了SVR的學習難度。另外，在實際的工程應用中，許多的研究對象是具有多輸入多輸出的復雜系統(tǒng)。包絡均值曲線左、右預測覆蓋方法是一個兩輸入兩輸出的非線性系統(tǒng)。因此，針對多輸入多輸出，提出了同倫-最小二乘支持向量雙回歸(HLS-SVDR)算法，算法如下。

設已知訓練集X={(x1,y1),…,(xk,yk)}，其中xk∈Rd，yk∈Rm，即有d維輸入，m維輸出。假設有l(wèi)個訓練樣本，即k=1,2,…,l。令yk,m表示yk中第m維輸出的值(m=1,2,…,N)。模型的目標函數(shù)和約束條件為

(m=1,2,…,n;k=1,2,…,l)

(20)

其中：ηk為第K個輸入變量所對應的所有維數(shù)輸出的擬合誤差和；ξk,m為第K個輸入變量所對應的第m維輸出的擬合誤差；C0為所有維輸出擬合誤差的懲罰系數(shù)；Cm為第m維輸出擬合誤差的懲罰系數(shù)。

其Lagrange函數(shù)形式為

(21)

引入同倫算法，即對于兩個函數(shù)f(x)和g(x)引入一個嵌入變量λ，構造一個新的函數(shù)φ(x,λ)=(1-λ)f(x)+λg(x)，其中，參數(shù)λ∈(0,1)。由于具有良好的收斂性，能夠有效地求解非線性代數(shù)方程組，在科學和工程中得到廣泛應用[25]。

在多輸出最小二乘支持向量回歸算法中，引入同倫嵌入變量λ可以將懲罰參數(shù)區(qū)間(0,∞)轉化為有限區(qū)間(0,1)，大大縮短懲罰參數(shù)的優(yōu)化時間。以此構造HLS-SVDR算法，目標函數(shù)和約束條件為

(λm∈(0,1);m=1,2;k=1,2,…,l)

(22)

其Lagrange函數(shù)形式為

(23)

分別對ωm，bm，ξk,m，ηk，βk和αk,m求偏導并令其為零，可得

(24)

根據式(24)，消除ωm,bm,ξk,m和ηk，可得Lagrange乘子式αm=(α1,m,α2,m,…,αl,m)T，β=(β1,β2,…,βl)T和偏置bm的求取方程，即

其中：K(·,·)為第m維輸出的核函數(shù)。

令Km=(Km(xk,xk′))l×l，E=(1,1,…,1)T∈Rl，ym=(y1,m,y2,m,…,yl,m)T，I為l×l的單位矩陣，則有

(26)

求解上式，得到同倫-多輸出最小二乘支持向量回歸第m維輸出的函數(shù)形式為

(27)

圖10 向量轉換為張量Fig.10 Converting the vector to tensor

依據一維向量轉為二階張量的方法，將訓練集X重構為兩輸入兩輸出的訓練集X′

(xi∈Rd-1;l=i×d)

(28)

將重構的訓練集X′和X″分別作為同倫-最小二乘支持向量雙回歸算法的訓練樣本集和測試樣本集。此外，核函數(shù)的選擇也決定了訓練模型的精確度，這里將使用最常用Gauss徑向基核

Km(xk,xk′)=exp(-‖xk-xk′‖2/2σ2)

(29)

由于σ∈(0,∞)，難以優(yōu)化。在此，同樣采用同倫算法，對核函數(shù)進行改進

K1(xk,xk′)=exp(-‖xk-xk′‖2/2h2)

K2(xk,xk′)=exp(-‖xk-xk′‖2/2(1-h)2)

(30)

其中：h∈(0,1)，使核參數(shù)h的優(yōu)化時間大大縮短，同時也提高了最優(yōu)核參數(shù)值的選擇精度。

將同倫-最小二乘支持向量雙回歸與最小二乘支持向量雙回歸對振動信號CEEMD分解的包絡均值m1(t)做實驗對比。使用張量方法，將m1(t)數(shù)據序列重構為X′和X″。使用量子遺傳(QGA)算法對最小支持向量雙回歸的懲罰參數(shù)C0，Cm及核參數(shù)σ及同倫-最小支持向量雙回歸的懲罰參數(shù)C0，λm及核參數(shù)h進行最優(yōu)選擇。實驗結果如圖11，12和表1所示。表1是對圖11中的HLSSVDR算法和圖12中LSSVDR算法在訓練和測試精度、運行時間進行的對比分析。根據實驗結果可知，基于HLSSVDR算法在預測精度及運算時間上都優(yōu)于LSSVDR算法。

表1LSSVDR和HLSSVDR運行時間及預測誤差

Tab.1TheelapsedtimeandpredictionerrorofLSSVDRandHLSSVDR

預測方法位置樣本均方誤差運行時間/sLSSVDR左端點右端點訓練集8.814 9×10-70.012 0測試集1.220 3×10-60.011 7訓練集1.240 9×10-60.011 4測試集1.279 7×10-60.012 2HLSSVDR左端點右端點訓練集1.241 3×10-80.004 21測試集4.731 1×10-70.003 37訓練集8.721 4×10-70.003 07測試集2.641 2×10-70.003 51

圖11 LSSVDR左、右雙向預測結果及誤差結果Fig.11 The left and right prediction and error results by LSSVDR

為了進一步驗證HLSSVDR在抑制端點效應方面的有效性，使用該方法對圖1中x(t)兩端的包絡均值進行預測覆蓋，實驗結果如圖13所示。圖13(a)使用HLSSVDR方法對信號兩端的包絡均值進行預測，獲得預測值，為了便于觀察兩端點的包絡均值預測值，將圖13(a)中的左、右端點的包絡均值預測值進行局部放大，如圖13(b)和圖13(c)所示。其中，圖13(b)中的A-B和圖13(c)中的D-E是存在的極值點插值獲得包絡均值曲線段和包絡均值預測曲線段?？梢钥吹?，左、右兩端的包絡均值預測值與真實值誤差很小，預測精度高。另外，圖13(b)中的B-C和圖13(c)中的E-F分別對應的是不存在極值點的包絡均值曲線段和包絡均值預測曲線段，將其與沒有進行端點效應處理的包絡均值曲線進行對比，從圖13(b)和圖13(c)中可看到存在很大的誤差。鑒于文中已經驗證HLSSVDR方法有較高的預測精度，能夠有效保證預測信號的準確度，因此EMD過程中保證有效的端點值才能保證信號分解的有效性。

圖12 HLSSVDR左、右雙向預測結果及誤差結果Fig.12 The left and right prediction and error results by HLSSVDR

圖13 包絡線均值的預測Fig.13 The mean prediction values of the envelope

4 滾動軸承外圈微故障實例分析

滾動軸承作為旋轉機械設備的關鍵組件，易發(fā)生軸承內/外圈裂紋或斷裂、滾動體磨損或破損等故障，對其進行早期檢測與診斷，能夠有效阻止故障發(fā)生，然而早期的故障特征頻率極其微弱和難以提取。筆者提出了完備CEEMD方法用于滾動軸承微故障特征提取。

圖14 滾動軸承故障實驗與采集裝置Fig.14 The fault testing and collecting device of rolling bearing

4.1 滾動軸承實驗臺簡介

如圖14所示，筆者選用美國凱斯西儲大學滾動軸承故障實驗臺、驅動端軸承型號為SKF6205的軸承外圈故障的振動信號作為分析對象。軸承外圈故障采用電火花加工單點損傷，直徑為0.173 4 mm(模擬軸承外圈微故障)，采樣頻率為12 kHz，采樣點數(shù)為4 096，電機轉速為1 797 r/min，軸承外圈故障頻率理論計算值為3.594 8×1 797/60=107.364 Hz。按照以上參數(shù)采集相關振動信號，如圖15所示。

圖15 滾動軸承外圈微故障振動信號Fig.15 The outer fault vibration signal of rolling bearing

4.2 微故障特征提取

為了進一步驗證完備CEEMD方法的有效性和優(yōu)越性，將其與傳統(tǒng)的EEMD，CEEMD方法進行對比分析。首先，使用多小波包對原始振動信號進行降噪預處理，提高信號的信噪比。如圖16所示。

圖16 基于DDWP降噪后的振動信號Fig.16 Denoising vibration signal by DDWP

然后，分別利用EEMD，CEEMD及完備的CEEMD方法對降噪后的滾動軸承外圈微故障振動信號進行分解，獲得IMF分量和殘余量。實驗結果如圖17～19所示。經過CEEMD、完備CEEMD方法分解得到的IMF分量的個數(shù)少于EEMD方法獲得的IMF分量個數(shù)。另外，CEEMD與完備CEEMD的IMF分量的個數(shù)相同，對比圖17中的IMF分量C7，C9，C10，C11與圖18，圖19中的IMF分量C7，C8，C9，C10，可以觀察到完備的CEEMD相應的IMF分量兩端的端點效應，相對于傳統(tǒng)的EEMD和CEEMD方法得到了抑制。

圖17 EEMD分解的IMF分量Fig.17 The IMF component by EEMD

圖18 CEEMD分解的IMF分量Fig.18 The IMF component by CEEMD

圖19 完備CEEMD分解的IMF分量Fig.19 The IMF component by Complete CEEMD

為了準確說明完備CEEMD算法獲得的IMF分量的特征信息，將分別對3種不同方法獲得的IMF分量進行Hilbert變換，獲得對應的滾動軸承外圈微故障振動信號的Hilbert邊際譜如圖20～22所示，分別為EEMD，CEEMD、完備CEEMD方法對應的滾動軸承外圈微故障振動信號的Hilbert邊際譜。圖20表明滾動軸承外圈微故障的特征頻率107.3Hz及相應的倍頻分量基本上被其他無關的特征頻率所淹沒，無法區(qū)分故障特征頻率與無關特征頻率。從圖21可觀察到表征滾動軸承外圈微故障的特征頻率為107.3 Hz，但相應的二倍頻特征頻率不能很好地辨析。從圖22中也可清晰地觀察到表征滾動軸承外圈微故障的特征頻率(107.3 Hz)及相應的二倍頻、三倍頻特征頻率。

圖20 EEMD分解的Hilbert邊際譜Fig.20 The Hilbert Marginal Spectrum by EEMD

圖21 CEEMD分解的Hilbert邊際譜Fig.21 The Hilbert Marginal Spectrum by CEEMD

圖22 完備CEEMD分解的Hilbert邊際譜Fig.22 The Hilbert Marginal Spectrum by Complete CEEMD

上述實驗分析表明，筆者提出的同倫-最小二乘支持向量雙回歸的保形分段三次樣條CEEMD的改進方法能夠更有效地提取滾動軸承外圈微故障的故障特征，能夠整體上有效解決EMD的包絡擬合過沖/欠沖、端點效應及模態(tài)混疊問題。

5 結 論

1) 通過對EEMD的包絡擬合過沖/欠沖、端點效應、模態(tài)混疊問題進行分析。在CEEMD方法的基礎上提出了同倫-最小二乘支持向量雙回歸的保形分段三次樣條的CEEMD方法。實驗結果表明，筆者所提方法能夠更有效解決EMD存在的3個主要問題。

2) 針對包絡擬合過沖/欠沖問題，提出了保形分段三次樣條插值的EMD改進方法，該方法利用保形分段法來構造具有二階逼近精度、分段少、運算量小的三次樣條插值方法來抑制包絡擬合過沖/欠沖問題，避免傳統(tǒng)插值方法導致的插值誤差隨著分解過程的持續(xù)進行而出現(xiàn)誤差不斷累積，造成嚴重誤差。實驗結果表明，與其他插值法相比，筆者提出的插值方法能更有效地解決該問題。

3) 針對端點效應問題，提出了基于數(shù)據驅動的同倫-最小二乘支持向量雙回歸的包絡均值預測覆蓋的方法。同其他延拓方法相比，該方法不依賴信號本身特性，而是采用機器學習的方法對信號兩端的包絡均值進行預測，將預測值覆蓋原有的隨機插值。實驗結果表明，該方法有較高的預測精度，收斂速度快、運算時間短，有效克服了基于數(shù)據預測抑制端點效應方法難以工程實現(xiàn)的問題。

4) CEEMD方法能夠很好地解決模態(tài)混疊問題，但也存在對原始信號添加的正、負白噪聲分解出的IMF+和IMF-分量個數(shù)不一定相等、導致少量模態(tài)混疊的問題。