基于優化變分模態分解與混沌分形融合的滾動軸承故障識別

孫 康， 金江濤， 李 春,2， 許子非

(1.上海理工大學 能源與動力工程學院，上海 200093；2.上海市動力工程多相流動與傳熱重點實驗室，上海 200093)

風能作為應用最為廣泛的可再生能源之一，以其穩定、豐富及易于開發等優勢，對碳中和和碳達峰起到關鍵作用[1]。軸承作為應用最廣泛的機械傳動零部件之一，在高性能艦船推進電機、先進航空發動機及大型風力發電機組等精密旋轉機械中具有至關重要的作用，其健康狀況直接影響設備性能[2]。滾動軸承是機械傳動系統的重要部件之一，能否正常運轉關系到整臺設備的性能和使用壽命[3]。但惡劣的工作環境導致軸承隨服役時間增加出現各種故障。為保證設備正常運作，需要實時監測軸承工作時產生的振動信號，確保維護人員對設備的狀態作出判斷并對設備進行及時維修[3]。

軸承振動信號呈現出非線性與非平穩性特征[3]，傳統頻譜分析法存在故障判別方式與準確性不足等問題，無法滿足實際工程需求。

向玲等[4]基于小波分解提取故障軸承特征信息。隋文濤等[5]采用經驗模態分解(EMD)方法獲取軸承故障特征頻率。李洋等[6]利用輔助經驗模態分解(EEMD)方法對信號降噪，完成故障識別。但上述方法都存在缺乏自適應性、端點效應及模態混疊等問題，無法處理復雜非平穩信號，為此Dragomiretskiy等[7]提出一種變尺度信號處理方法——變分模態分解(VMD)方法，該方法避免了上述問題，具有良好的魯棒性，被國內外學者廣泛應用于故障診斷領域。

但VMD方法需人為設置懲罰因子c和模態分解數K。c取值過小，易造成分解所得本征模態函數(IMF)帶寬過大；反之過小。若K取值過大，原信號出現過分解現象；反之原信號無法完全分離，存在模態混疊問題[8]。唐貴基等[9]根據包絡譜特征因子選擇VMD參數，但未考慮轉頻對包絡信號峰值的影響。張俊等[10]利用粒子群算法優化VMD參數，雖避免了人為經驗影響，但粒子群算法易陷入局部最優，誤差較大。

當結構復雜的設備出現故障時，從統計學角度分析，其振動信號具有自相似性特征，故可用混沌分形理論來描述信號的復雜性和不規則程度，揭示隱藏于隨機現象中的真實規律[11]。

筆者以軸承振動信號為研究對象，將Lyapunov指數λ視為特征值。λ>0，表明信號具有混沌特性，且混沌特性與λ線性相關，以λ最大值為迭代尋優條件，選取VMD分解參數，提出優化變分模態分解(OLVMD)方法，以分解軸承原始振動信號。不同于傳統故障診斷方法，OLVMD方法完全脫離頻域分析，診斷時間更迅速，形式更簡潔，準確率更高，具有獨特優勢。原信號經OLVMD方法分解得到若干IMF分量，但故障特征僅存于少數IMF分量中，大多為噪聲信號。考慮到混沌具有分維性、敏感性及隨機性等特點,可量化分析振動信號[11]，將各故障Lyapunov指數最大的前3個分量作為含故障信息的特征信號，重構組成新的振動信號。為準確進行故障識別，筆者引入分形理論，計算軸承在不同狀態下的混沌關聯維數，并采用擬合偏差平方和方法對傳統的關聯維數計算方法進行改進。從分形角度分析軸承振動信號，以其關聯維數為特征量，定量描述軸承工作狀態的變化，并分析了損傷軸承的實驗數據。

1 變分模態分解

VMD將待分析信號視為若干模態分量線性疊加而成。各IMF分量被定義為調幅調頻信號，表達式如下：

uk(t)=Ak(t)cosφ(t)

(1)

式中：t為時間；uk(t)為第k個IMF；Ak(t)為瞬時幅值；φ(t)為信號相位。

通過希爾伯特變換計算各IMF分量單邊譜,再估計中心頻率，構造變分模型：

(2)

式中：k=1,2,…,K；σ(t)為單位脈沖函數；ωk為中心頻率；f(t)為輸入信號；?表示卷積運算；?t表示偏導運算；j為虛數。

引入懲罰因子及Lagrange乘法算子將約束變分問題變為非約束變分問題，得到擴展Lagrange表達式。

求解過程中各分量中心頻率與帶寬不斷更新，直至滿足迭代停止條件：

(3)

迭代終止時，信號頻域特性已完成自適應分離,然后通過傅里葉逆變換轉換到時域。

2 算法改進

2.1 優化變分模態分解方法

基于VMD方法分解振動信號時，懲罰因子c和模態分解數K直接影響分解精度[8]。Lyapunov指數λ可定量反映非線性振動信號的特征信息，λ越大，故障特征越顯著[11]。故本文以λ為特征值，選取最佳參數組合[c,K]，以分解所得IMF分量中λ最大的IMF分量為特征信號。基于最大Lyapunov指數提出OLVMD方法，流程如圖1所示。

圖1 OLVMD方法流程圖

2.2 有效性驗證

為驗證OLVMD方法可有效降噪并提取信號中的特征信息，以一段噪聲序列為研究對象，分別采用OLVMD方法和VMD方法對其進行預處理，以分解所得IMF分量中λ最大的IMF分量為特征信號，其時域圖如圖2所示。其中OLVMD方法最佳分解參數組合為c=1 200，K=10；借鑒劉洋等[12]所提方法，根據EEMD方法遞歸特性確定模態分解數K，c與采樣頻率相同，即c=1 000，K=3。

圖2 噪聲信號處理前后對比

由圖2可知，采用OLVMD方法處理后的信號降噪效果顯著，為定量表示，計算2種方法處理后信號的信噪比(SNR，SNR)，結果如表1所示。SNR越高，降噪效果越好[13]。由表1可知，原信號含有大量噪聲，信號特征被淹沒，信噪比最高。相比VMD方法，經OLVMD方法處理后的信號信噪比較低，降噪效果更明顯。OLVMD方法在處理非線性信號時，更易獲取純凈故障信息。

表1 有效性驗證

2.3 關聯維數計算

根據文獻[13]可得到關聯維數D(m,r)：

(4)

式中：m為嵌入維數；Cm(r)為關聯積分；r為超球半徑;r1

為此，采用擬合偏差平方和方法選取線性區域，以提高計算精確度，將lnCm(r)-lnr曲線中間部分視為直線，兩側為曲線。將整段曲線橫坐標記為[d1,d2]，其中直線段為[lnr1,lnr2]，滿足d1≤lnr1≤lnr2≤d2。具體擬合方法如下：

(1)d1≤lnr≤lnr1時，以多項式曲線y=f1(x)進行擬合，擬合偏差平方和記為ε1：

(5)

(2) lnr1

(6)

(3) lnr2≤lnr≤d2時，以直線y=f2(x)進行擬合，擬合偏差平方和記為ε3：

(7)

擬合偏差平方總和(即誤差)ε=ε1+ε2+ε3，當其值達到最小值時，[lnr1，lnr2]即為最佳線性區域。采用多項式曲線對lnCm(r)-lnr曲線進行擬合時，一般冪次取2較為合適。

2.3.1 仿真信號驗證

以經典混沌動力系統Lorenz吸引子為例，采用改進G-P算法求取關聯維數，其包含3個狀態變量的常微分方程：

(8)

式中：Pr為普朗特數；Ra為規范化的瑞利數；b的大小與吸引子軌跡形狀有關。

基于四階龍格庫塔法對上式求解，其中Pr=10，b=3，Ra=25，起始點(x,y,z)取(0,2,9)，迭代步長h為0.01，得到一個5 000×3的時間序列，選擇X變量進行分析，Lorenz吸引子軌跡圖及X變量時間序列如圖3所示。

(a) Lorenz吸引子軌跡圖

采用互信息函數法計算得到延遲時間為10，嵌入維數m從上至下依次為2,4,6，…，30所對應的雙對數曲線lnCm(r)-lnr如圖4所示。

圖4 ln Cm(r)-ln r曲線(無噪聲)

由圖4可知，尺度較小時，曲線分散且不規則；隨著尺度增大，曲線與lnr軸逐漸平行，并向同一方向收縮；處于中間尺度時，曲線接近于斜直線，可描述系統分形特征，分形維數在此區間具有意義。

基于擬合偏差平方和方法對圖4中線性較好的部分采用最小二乘法進行擬合，計算各嵌入維數對應的曲線斜率，以最佳嵌入維數對應的曲線斜率為該系統的關聯維數D(GP)，結果如圖5所示。

由圖5可知，嵌入維數較小時，嵌入維數與關聯維數呈正相關，當嵌入維數大于16時，關聯維數趨于穩定，不再增長。因此，最佳嵌入維數為16，其對應的曲線斜率3.231 0為Lorenz吸引子關聯維數。

圖5 無噪聲時關聯維數與嵌入維數的關系

2.3.2 抗噪性能分析

考慮到外部環境影響，實際軸承振動信號難免受到噪聲干擾，為驗證關聯維數具有強魯棒性。在Lorenz吸引子中加入不同大小的高斯白噪聲，計算系統關聯維數。為準確揭示關聯維數變化規律，規定若其波動范圍>0.01，則表明關聯維數受噪聲干擾，分析過程如下：

(1) 不斷減小信噪比，當SNR>6 dB時，Lorenz系統嵌入維數值穩定在16，關聯維數波動極小，可忽略不計，表明關聯維數不受噪聲干擾。

(2) 1 dB≤SNR≤6 dB時，系統嵌入維數和關聯維數發生細微改變。以SNR=4 dB為例，其lnCm(r)-lnr曲線與關聯維數如圖6所示。

(a) ln Cm(r)-ln r曲線

(b) D(GP)-m曲線

計算不同信噪比對應的最佳嵌入維數與關聯維數，結果如表2所示。

表2 不同信噪比時最佳嵌入維數與關聯維數的變化

由表2可知，隨著信噪比增大，關聯維數不斷減小，當信噪比增大到一定值時，關聯維數不再變化，表明關聯維數具有一定的抗噪性。

(3) 當SNR<1 dB時，lnCm(r)-lnr曲線斜率不再收斂于一個穩定值。以SNR=-4 dB為例，雙對數曲線與關聯維數如圖7所示。

(a) ln Cm(r)-ln r曲線

由圖7可知，曲線斜率不再收斂，無法判定最佳嵌入維數及其對應的關聯維數。

通過上述分析可知，關聯維數具有一定抗噪性，且其存在2個臨界值Q1和Q2，其中Q1≤Q2。當SNRQ2時，振動信號關聯維數不會產生波動，信號不受噪聲影響，關聯維數可準確表達信號的狀態信息。

3 結果與分析

3.1 實驗數據集

利用西安交通大學軸承數據中心提供的實驗數據[14]驗證所提分形故障診斷方法的有效性及準確性。實驗平臺如圖8所示，該裝置以交流電動機為動力源，通過轉速控制器改變軸承轉速，加速度傳感器位于軸承的水平及豎直方向。

圖8 軸承實驗平臺裝置圖

所研究的滾動軸承相關參數如表3[14]所示。

表3 軸承參數

設置采樣頻率為25.6 kHz，對轉速為2 100 r/min和2 250 r/min 2種狀態4種工況下的軸承振動數據進行采集。具體參數如表4所示。

表4 軸承加速壽命試驗工況

3.2 特征選取

原信號在4種工況下的時域及頻域圖見圖9，其中f為頻率。觀察圖9可知，各故障振動信號時域圖和頻域圖雖存在一定差異，但只能定性分辨幾種故障的不同，無法直接判斷故障類型。為提取有效故障特征，基于相空間重構法繪制4種故障三維吸引子軌跡，分析各振動信號的混沌特性，采用Cao方法與互信息函數法計算最佳延遲時間和嵌入維數，結果如表5所示，吸引子軌跡如圖10所示。

表5 最佳延遲時間和嵌入維數

由圖10可知，4種故障相圖均呈冗雜毛球狀，吸引子軌跡形態各異，呈非完全周期性也非完全隨機響應，表明各故障信號具有顯著混沌特性。但因早期振動信號微弱，且受噪聲干擾，無法通過吸引子軌跡判斷故障類型。

為對原信號降噪提取特征信息，采用OLVMD方法分解原信號，基于Lyapunov指數選取最佳分解參數[c,K]，其中K∈[3,10],c∈[1 000,2 000]，各信號Lyapunov指數如圖11所示。由圖11可知，以最大Lyapunov指數為尋優條件得到的最佳分解參數組合[c,K]如表6所示。

表6 最佳分解參數

分解各故障信號，各IMF分量時域圖見圖12。計算各IMF分量的Lyapunov指數，結果見圖13。

圖13 Lyapunov指數

由圖13可知，4種故障的Lyapunov指數各不相同，但均為正數，表明各故障信號具有不同程度的

混沌特性，以Lyapunov指數λ為特征值，λ越大，故障特征越豐富，取λ最大的前3個分量重構故障信號。

為定量表示OLVMD方法的降噪效果和特征提取能力，求取重組信號吸引子軌跡并采用流行學習方法進行降維，可視化故障信息，結果見圖14和圖15。

(c) 保持架磨損

(d) 外圈磨損

由圖14可知，相較原信號混沌圖，4種故障吸引子軌跡收縮，混沌特性增強，信號特征更顯著。計算各故障Lyapunov指數，可知由于保持架磨損對軸承運轉情況影響最小，其Lyapunov指數最小；多故障導致的混合損傷的Lyapunov指數最大。

圖15 降維可視化圖

由圖15可知，各故障類型間距增大，不存在模態混疊現象，4種故障可明顯區分，表明所提方法能很好地提取特征信息。

3.3 關聯維數

關聯維數描述了振動信號中各點之間的關聯性，對信號波動極為敏感，可量化分析故障的狀態特征，具有簡潔、高效、快速的特點，各故障關聯維數與嵌入維數的關系如圖16所示。

觀察圖16可得各故障最佳嵌入維數與關聯維數。為準確進行故障識別，需建立故障關聯維數區間，通過計算未知故障振動信號關聯維數，判斷其所處區間，完成故障診斷。

考慮到振動信號本身波動性影響，會出現偏差較大的關聯維數點，需剔除這些點以提高故障識別準確率。計算10次關聯維數并求其平均值以減小誤差，結果如表7所示。

表7 軸承各故障關聯維數

由表7可知，各故障關聯維數D的區間為：

為更直觀地表示各故障關聯維數區間，繪制狀態區間圖，如圖17所示。由圖17可知，各故障關聯維數對應的狀態區間存在顯著差別，互不相交,可準確完成故障識別。

圖17 軸承各故障關聯維數區間

綜上所述，對軸承進行故障識別的具體步驟如圖18所示,其中Dmin和Dmax分別為關聯維數的最小值和最大值。

圖18 故障識別流程圖

3.4 有效性驗證

為驗證所提方法的有效性，從4種故障中隨機選取數據進行驗證，為避免出現偏差，每種故障各選取3組數據，計算關聯維數，結果如表8所示。

表8 未知故障信號的關聯維數

對比表8中故障類型及其關聯維數區間，4種故障關聯維數平均值均在對應的關聯維數區間內，可準確完成故障識別。

3.5 對比驗證

為驗證所提方法的優越性，基于EEMD和VMD方法分解原信號并重構，采用相同方法求取各故障信號的關聯維數區間，結果如圖19所示。

圖19 不同方法各故障關聯維數區間

由圖19可知，采用EEMD與VMD方法所得關聯維數區間均出現混疊現象，無法準確判斷軸承故障類型，具有一定的局限性。

3.6 齒輪箱驗證

為驗證所采用的方法在實際工程機械軸承中的故障診斷效果，采用美國國家可再生能源實驗室(NREL)5 MW近海風力機齒輪箱數據進行驗證[15]，對齒輪斷齒面積為100%～25%(嚴重性分為A、B、C和D)4類不同故障情況進行分析，采用OLVMD方法進行降噪并選取敏感分量重構信號，計算各故障齒輪箱振動信號的關聯維數，取10次數據的平均值，結果如表9所示。

表9 齒輪箱各故障的關聯維數

由表9可知，各斷齒故障關聯維數區間為：100%斷齒D∈[5.200 2,5.318 4]；50%斷齒D∈[3.889 6，3.959 8]；25%斷齒D∈[4.375 7，4.459 3]；12.5%斷齒D∈[2.995 4，3.187 0]。為更直觀地表示各故障關聯維數區間，繪制狀態區間圖，如圖20所示。

圖20 齒輪箱各故障的關聯維數區間

由圖20可知，各斷齒故障關聯維數所對應的區間存在顯著差別，互不相交,表明所提方法對風力機齒輪箱的故障具有較好的診斷效果。

4 結 論

(1) 滾動軸承受噪聲干擾，特征信息被掩蓋，無法通過原信號判斷故障類型。

(2) 利用OLVMD方法可有效降噪并提取故障信號特征信息。

(3) 基于分形維數的故障診斷完全脫離頻域分析，利用振動信號自相似性可有效區分軸承運行狀態。

(4) 分形故障診斷操作簡單，診斷時間短，符合當前智能遠程故障診斷的發展趨勢。