VMD樣本熵特征提取方法及其在行星變速箱故障診斷中的應用

楊大為， 馮輔周， 趙永東， 江鵬程， 丁 闖

(陸軍裝甲兵學院 機械工程系，北京 100072)

行星變速箱廣泛應用于大型工程機械中，具有傳動比大，機構緊湊，傳遞效率高等優點，其工作環境惡劣，齒輪故障時有發生，不及時進行處理，將嚴重影響變速箱功能。行星變速箱傳動機理復雜，振動傳感器測得信號干擾和衰減嚴重，為典型的非線性非平穩信號，行星輪故障和太陽輪故障難于區分。而現有的很多研究還停留在套用定軸變速箱故障診斷方法，不能針對行星變速箱特點進行有效的分析[1]。某型行星變速箱為多檔位行星變速箱，同時存在多個定軸輪系和行星輪系，行星輪系還存在復合框架和兩級復合行星排，多數檔位需要由2～3個行星排同時工作，而且每個行星排至少有3個行星齒輪，并且結構非常緊湊，箱體為圓柱體并且表面存在大量散熱筋，內部有三個離合器和兩個制動器共五個操縱件，傳感器測點選擇十分困難。這些自身原因給該型行星變速箱的故障診斷工作帶來了更大的困難。

Dragomiretskiy等[2]提出一種全新的自適應非遞歸信號處理算法變分模態分解，迭代尋找變分問題最優解以確定每個IMF分量的中心頻率和帶寬，實現信號頻譜的劃分及各分量的有效自適應分解。不同于EMD、EEMD和LMD傳統遞歸算法，VMD算法擺脫了傳統信號分解的遞歸篩分剝離模式的約束，有著堅實的數學基礎，可以緩解模態混疊和邊界效應，還有運算效率高和魯棒性強的優勢[3]。樣本熵是一種對近似熵進行改進的計算時間序列復雜度的算法，非常適合處理非線性非平穩信號[4]。與近似熵對比，樣本熵提高了統計精度，降低了對時間序列長度的依賴，具有更好的一致性。系統不同的運行狀態對應著不同的樣本熵值，所以樣本熵可用于表征系統的運行狀態[5]。

本文針對傳統樣本熵的不能有效提取行星變速箱故障特征的問題，提出結合VMD與樣本熵的特征提取方法，并將其應用于某型行星變速箱。首先依據中心頻率觀察法選取VMD的分解尺度和依據樣本熵最小原則選取二次懲罰因子，再利用VMD將振動信號分解為若干個IMF分量，然后依據敏感度最大原則確定VMD分解各IMF與原信號相關系數閾值，選取大于相關系數閾值的IMF分量重構信號，最后對重構后的信號計算樣本熵，依據熵值判斷行星變速箱的運行狀態。

1 VMD樣本熵特征提取方法

1.1 VMD算法原理

VMD是一種區別于EMD、EEMD和LMD的非遞歸新算法，能夠自適應地分解信號，得到k個中心頻率為ωk的模態函數uk，其中k為預設尺度數。VMD分解過程的實質就是構造和求解變分問題的過程[6]。

為估計各模態函數uk的帶寬，如下步驟構造變分問題：

(1)通過對每個uk進行Hilbert變換，獲得其相應的單邊頻譜；

(2)通過指數混合調制方法，將每個uk的頻譜移動到各自估算的中心頻率；

(3)依據高斯平滑度和梯度平方準則對信號進行解調估計各uk的帶寬。

約束變分問題則可表示為：

(1)

(2)

式中：uk：={u1,u2,…,uk}為各模態函數；ωk：={ω1,ω2,…,ωk}為各模態函數的中心頻率。

為求解該約束變分問題，引入拉格朗日乘子λ(t)和二次懲罰因子α，將其變為非約束變分問題，得到增廣拉格朗日表達式為：

L({uk},{ωk},λ)：

(3)

迭代步驟如下：

(2)迭代次數n=n+1；

(3)fork=1∶K

(4)

(5)

(4) 根據式(6)，對于所有ω≥0，進行雙重提升，更新λ；

(6)

式中：τ為噪聲容限，可設為0以達良好去噪效果。

1.2 樣本熵原理

樣本熵適用于了量化非線性非平穩信號的復雜度，具有無需自我匹配度、計算快、精度高的優點[7]。樣本熵可以衡量行星變速箱振動信號復雜程度，如果振動信號成分單一，周期性越明顯，信號噪聲干擾越少，信號復雜程度越低，樣本熵值越小，反之信號噪聲干擾越多，信號復雜程度越高，樣本熵值越大。

給定序列X{x(n),n=1,2,…,N}，計算樣本熵步驟如下：

(1) 對序列X進行相空間重構，獲得矩陣Y,m為模式維數。

(7)

1≤i,j,K≤N-m+1

(2) 計算向量Y(i)與Y(j)中對應元素最大差值，將其絕對值定義為兩者的距離d(i,j)

d(i,j)=maxx(i+k)-x(j+k)
0≤k≤m-1
1≤i,j≤N-m+1,j≠i

(8)

(9)

(10)

(4) 維數增加到m+1，獲得一組m+1維向量，重復(1)～(3),得到Bm+1(r)，

(5) 理論上，樣本熵定義為

(11)

當N為有限值時，樣本熵的估計值為

(12)

由式(12)可知，模式維數m和相似容限r的取值都會影響最后計算所得樣本熵值。根據文獻[8]的研究結果和試驗對比，應取m=2，r=0.15×Std(Std為信號的標準差)。

1.3 特征提取性能評估

為量化特征參數對行星變速箱狀態的分類能力，引入雙樣本Z值檢驗法。雙樣本Z值檢驗法能有效地對兩類樣本在統計上的差異進行評估，計算結果Z值可以作為特征參數的分類能力進行的評價指標，在相同情況下，Z值越大，兩類樣本的分類距離越大，說明特征參數的分類能力越強。

X1{x11,x12,…,x1j}和X2{x21,x22,…,x2j}分別是兩類樣本，定義特征參量的Z值(敏感度)為

(13)

1.4 結合VMD與樣本熵的特征提取

結合VMD與樣本熵的特征提取具體方法流程如下：

(1) 依據基于信號頻譜特性和VMD分解能力的中心頻率觀察法確定VMD分解尺度，依據重構信號樣本熵值最小原則確定VMD二次懲罰因子，使用VMD分解振動信號。

(2) 計算分解所得的各IMF分量與原信號的相關系數，依據敏感度最大原則確定相關系數閾值，選取大于等于閾值的IMF分量進行重構。

(3) 計算重構信號的樣本熵值，作為故障特征參數，使用敏感度指標對特征參數分類能力進行評估，驗證特征參數的有效性。

2 實例分析

2.1 行星變速箱故障模擬試驗介紹

實驗對象為某型三軸式離合器換擋三自由度行星變速箱，其主動軸齒輪和中間軸齒輪為定軸齒輪，主軸上有三個行星排，K1排為外嚙合雙行星排，其余兩排為簡單行星排。故障模擬試驗臺的組成如圖1(a)所示，試驗工況設定輸入轉速為1 500 r/min，負載為900 N·m，均為試驗臺所能達到最大轉速和最大加載，此時齒輪振動最為劇烈，便于提取故障特征。變速箱擋位設置在Ⅳ檔，此時K3太陽輪固定，由內齒圈傳入動力，內齒圈帶動行星輪轉動，從而使行星架轉動將動力輸出。行星變速箱此時只有K3行星排傳動方式為行星傳動， K1和K2行星排整體旋轉，可以排除其他行星排的干擾。傳感器測點安置在右側輸出軸端的箱體內部靠近K3內齒圈的位置，如圖1(b)，該測點距振源近，受傳遞路徑干擾小，能更好地采集齒輪故障的振動信號。故障設置在Ⅳ檔K3行星排的太陽輪Z=30的某輪齒齒面上，利用線切割將其某個齒切除來模擬斷齒故障，實物圖為圖1(c)。

圖1 行星變速箱故障模擬試驗有關情況Fig.1 Introduction of planetary gearbox fault simulation experiment

行星輪系的運動方式不同于定軸輪系為復合運動，行星輪既自轉又公轉，采集的數據應該包括相同個行星齒輪回到初始位置的完整周期，這樣得到的相同工況的樣本信號對應著行星輪系運行的相同物理過程，此時計算的樣本熵才具有意義，否則采樣所取完整周期數不同，測得信號對應的物理過程不同，故障產生的沖擊數也不同，這樣將多個振動信號的樣本熵做橫向比較將失去意義。所以針對行星輪系的特殊情況，需要考慮行星齒輪回到初始位置周期的問題。

表1 K3齒輪參數

依據齒輪相關參數和理論計算可得K3行星排參數如表1，由表知，行星輪自轉頻率為行星架轉頻的2倍，即行星架每轉一圈，行星輪就要自轉兩圈，行星輪與太陽輪嚙合兩個周期30個齒，剛好是太陽輪齒數的一倍，太陽輪嚙合一個周期，行星輪剛好旋轉到初始位置。即行星架轉一圈，行星輪剛好旋轉到初始位置，一個完整周期歷時t=1/17.65≈0.056 6 s。為保證采集信號的可用性，本次試驗采樣頻率20 kHz，每次采樣取30個周期，間隔30個周期，即每次采樣1.7 s，間隔1.7 s，連續采集52組數據用于處理。

實測正常和太陽輪斷齒信號時域波形如圖2所示，從時域看，故障信號沖擊脈沖不凸顯，并且沒有明顯的周期；從頻域看，實測信號干擾較多，頻帶成分復雜，低頻段定軸齒輪嚙合頻率及其邊頻帶和噪聲干擾嚴重，找不到太陽輪斷齒的故障特征頻率，僅從時域和頻域都不能區分正常和故障信號。

圖2 實測信號時域及頻譜圖Fig.2 Measured signals in time domain and corresponding spectrum

2.2 VMD分解重構參數優化選取

使用VMD時需要選定預設尺度數K和二次懲罰因子α，預設尺度數K和二次懲罰因子α都會對分解結果產生影響，由于實測信號復雜多變，使用VMD算法的難點和關鍵在于如何選定合適的Κ和α值。

目前，使用VMD時一般選擇默認的Κ和α值，不能保證VMD參數選取最優，極大限制了VMD的性能。文獻[11]采取粒子群算法對VMD參數進行優化選取，取得了比較好的效果，但該優化過程繁瑣復雜，粒子群尋優需要很長的時間，不利于在線監測以及實時診斷，且所處理信號分解后不涉及信號重構，沒有考慮相關系數閾值對信號重構的影響。因此，本文提出一種簡單快捷且考慮相關系數閾值的方法對VMD算法參數進行優化選取。

2.2.1 基于中心頻率觀察法的分解尺度選取

由于VMD分解得到的IMF分量的中心頻率是由低頻至高頻分布的，如預設尺度數K值從小到大取值，則最后一層IMF分量的中心頻率首次達到最大值時，將不會出現分解不足的問題，若隨著K值增大，最大中心頻率仍然保持相對穩定，則可認為此時K的取值為最佳取值[12]。

當K取不同值時，取一組太陽輪斷齒故障信號進行VMD分解，得到IMF分量的中心頻率如表2所示，由表2得，當K=5時，IMF1取到中心頻率最小值950 Hz并當K≥5時，IMF分量的中心頻率最小值趨于穩定。當K=8時，IMF8取到中心頻率最大值8 028 Hz并當K≥8時中心頻率最大值趨于穩定。當時，IMF1中心頻率對應行星變速箱主動軸定軸齒輪的嚙合頻率，是頻譜低頻段的主要特征頻率， IMF8中心頻率對應著信號頻譜最后一個明顯峰值，符合信號實際情況。當K≥8時中心頻率最大值和最小值都保持相對穩定，不再出現新的中心頻率最大值或者最小值，從而保證VMD分解不會遺漏更高或者更低的中心頻率，則認為此時VMD分解能力最佳，故本次取預設尺度K=8。

2.2.2 基于樣本熵值最小原則的二次懲罰因子選取

VMD算法中二次懲罰因子α越小，分解得到的各IMF分量的帶寬越大，α越大則帶寬越小[12]。使用VMD算法處理信號后，如果信號包含的故障信息較少，與故障相關的周期性沖擊脈沖不凸顯，信號干擾和噪聲較強，樣本熵值較大。如果包含的故障信息較多，出現與故障相關的周期性沖擊脈沖，信號的干擾和噪聲較弱，呈現出較強的規律性，樣本熵值較小，故選取最優的二次懲罰因子α應使重構信號的樣本熵值最小。在使用VMD算法時，采用不同的懲罰因子對信號進行VMD重構后求樣本熵將得到不同的結果。本文在預設尺度K=8的條件下采用不同的懲罰因子處理試驗52組

表2 取不同K值時VMD分解后各IMF分量中心頻率

實驗數據，所得樣本熵均值如圖3，可知當α=2 500時，重構信號的稀疏性最強，樣本熵值到達最小值，故本文取二次懲罰因子α=2 500。

圖3 二次懲罰因子和VMD重構信號樣本熵值關系圖Fig.3 Relationship of secondary penalty factor and sampleentropy of VMD reconstructed signal

2.2.3 基于敏感度最大原則的相關系數閾值選取

為能更好地提取特征，需要選取分解結果中有效的IMF進行重構。相關系數可以衡量各IMF分量和原信號的相關性，各IMF分量與原信號的相關系數越小者，表明它與原信號的相關性越低，包含原信號的可用信息越少，可能屬于噪聲干擾成分，反之，表明該IMF中包含原信號的可用信息越多[13]。相關系數閾值過小不能有效去除噪聲干擾或對故障不敏感的IMF分量，閾值過大可能會遺漏對故障較為敏感但比較微弱的IMF分量，因此相關系數閾值選取對信號重構結果有很大影響。根據與原信號的相關系數大于等于給定閾值的原則，篩選出有用的IMF分量，從而進行信號重構，采用區分故障的敏感度指標最大的原則確定相關系數閾值，使重構信號的樣本熵對故障狀態最為敏感。

圖4為斷齒信號VMD分解所得8個IMF分量，表3為各IMF分量與原信號的相關系數。本文在分解層數K=8和懲罰因子α=2 500的條件下，采取不同的相關系數閾值對VMD處理后的52組信號進行重構，然后對重構信號求樣本熵，并用敏感度指標進行評估。由表4可知，在相關系數閾值取r=0.3時，敏感度指標最大，說明此時VMD樣本熵區分斷齒故障狀態與正常狀態的能力最強。因此，取相關系數閾值r=0.3。

圖4 VMD分解結果Fig.4 The result of VMD

分量相關系數IMF10.332 5IMF20.681 6IMF30.520 5IMF40.315 4IMF50.284 0IMF60.184 2IMF70.136 7IMF80.130 1

表4 不同相關系數閾值重構樣本熵

依據確定的VMD參數對太陽輪斷齒故障信號進行重構，得重構前后信號局部頻譜如圖5。由圖5可知，重構后的信號頻譜更加清晰，噪聲干擾得到有效抑制，行星輪系嚙合頻率(529 Hz)及其二倍頻(1 058 Hz)、三倍頻(1 587 Hz)更加凸顯，證明了VMD重構方法的有效性。但其頻譜主要頻率成分仍為變速箱定軸部分特征頻率，行星輪系故障特征頻率仍然較為微弱，需要結合樣本熵進一步分析處理。

圖5 信號頻譜圖Fig.5 Signals frequency spectrum

2.3 結果分析

使用VMD算法和EEMD算法[14-17]對52組實測信號進行重構，然后分別計算齒輪正常狀態及齒輪斷齒故障狀態的樣本熵值(直接樣本熵)、EEMD重構后的樣本熵值(EEMD樣本熵)、VMD重構后的樣本熵值(VMD樣本熵)。圖6和表5給出了輸入轉速為1 500 r/min和負載為900 N·m時齒輪正常狀態和太陽齒輪斷齒狀態的樣本熵計算結果。

表5 實測信號樣本熵對比

圖6 不同方法樣本熵對比圖Fig.6 Comparison of sample entropy with different methods

由圖6和表5知，斷齒故障狀態的樣本熵值小于正常狀態的樣本熵值。因為在行星變速箱正常工作狀態時，齒輪振動幅度較小，振動信號相對復雜，并無明顯規律性，故其樣本熵值較大；當發生齒輪斷齒故障時，振動信號中會產生與齒輪故障相關的有規律的沖擊脈沖，導致信號復雜程度降低，樣本熵減小，證明樣本熵能衡量振動信號的復雜程度，描述行星變速箱齒輪運行狀態。

行星變速箱振動信號成分復雜，信號的直接樣本熵通常無法準確區分故障，故樣本熵不能直接作為判斷故障的依據。而VMD樣本熵和EEMD樣本熵值較直接樣本熵值小，表明VMD樣本熵和EEMD樣本熵在保留不同狀態信號特征差異的基礎上有效降低了信號中的隨機成分，同時VMD樣本熵的均值和標準差更小，說明了VMD樣本熵較EEMD樣本熵具有更好的穩定性，在更好地去除背景噪聲的同時凸顯了沖擊成分。

通過敏感度指標評估三種方法的分類能力，結果表明，相同試驗工況下VMD樣本熵敏感度最大，說明VMD樣本熵對斷齒故障狀態和正常狀態的分類能力最強，對齒輪故障狀態最敏感，可以作為評判齒輪狀態的依據。

為檢驗VMD算法在低采樣頻率的優勢，對實測數據進行重采樣，分別使用VMD算法和EEMD算法處理采樣頻率為625、1 250、2 500、5 000、10 000、20 000 Hz的信號，結果如圖7所示。從對相同信號的處理結果可以看出，在各個采樣頻率下，VMD算法得到的輸出信號信噪比和敏感度都明顯高于EEMD算法，且VMD算法處理信號的時間明顯少于EEMD算法，這就意味著VMD算法可在采樣頻率較低的情況下仍舊能夠獲得高信噪比的輸出信號，結合樣本熵算法能更好地區分正常和斷齒故障，降低對硬件的要求，減少運算量。

圖7 不同采樣頻率下兩種重構方法對比圖Fig.7 Comparison of results from two reconstruction methodsusing different sampling frequencies

為檢驗VMD樣本熵方法的適用性，使用同樣方法設置K3行星排Z15行星齒輪斷齒故障，選取處理行星齒輪故障信號的VMD參數，分別計算信號的VMD樣本熵和EEMD樣本熵，如圖8所示。由圖8可知，VMD樣本熵對行星齒輪斷齒故障有較好的區分能力，而EEMD樣本熵方法則不能有效區分，存在一定的局限性，證明了VMD算法處理非線性非平穩復雜信號相對于EEMD算法的優勢。由于K3行星排有6個行星輪，每個行星輪斷齒產生的沖擊相對于太陽輪斷齒產生的沖擊弱，加上行星輪復合運動的影響，使信號復雜程度增加，樣本熵值增大，且樣本熵的穩定性相對降低。

表6 不同工況VMD樣本熵

圖8 不同狀態信號的VMD和EEMD樣本熵圖Fig.8 VMD and EEMD sample entropy of different condition signals

使用本文方法處理不同工況下的太陽齒輪斷齒信號，結果如表6。由表6可知，在行星變速箱不同轉速、不同負載的工況下，太陽齒輪斷齒信號的VMD樣本熵均小于正常信號的VMD樣本熵，具有較高的敏感度。并且轉速越高，負載越大，齒輪斷齒產生的沖擊越明顯，采集到的信號越規律，樣本熵值越小，與正常信號的區別越大，VMD樣本熵區分齒輪狀態的敏感度越高，進一步證明了本文方法可以有效應用于行星變速箱故障特征提取。

3 結 論

本文采用VMD樣本熵方法提取了行星變速箱齒輪故障特征，通過行星變速箱齒輪故障模擬實驗和數據處理結果表明：

(1)本文采用的VMD參數選取方法快捷有效，可以應用于VMD分解參數選取，結合VMD和樣本熵的特征提取方法相對于EEMD樣本熵具有敏感度高，計算耗時少，信噪比高，不依賴采樣頻率的優勢，值得深入研究。

(2)本文考慮了行星齒輪的周期問題，保證了數據的可用性，并在不同試驗工況下均可獲得良好的效果，證明了方法的有效性，并初步區分了太陽輪斷齒和行星輪斷齒故障，可以應用于行星變速箱運行狀態的特征提取，為下一步的故障識別打下基礎。