程軍圣, 李夢君, 歐龍輝, 楊 宇

(湖南大學(xué) 汽車車身先進設(shè)計制造國家重點實驗室，長沙 410082)

齒輪作為旋轉(zhuǎn)機械中不可或缺的一部分，其故障信號往往表現(xiàn)出非線性非平穩(wěn)特性，針對這種非線性非平穩(wěn)特性，常采用時頻分析方法對這類信號進行預(yù)處理。常用的時頻分析方法有經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)等方法[1-3]。但是EMD方法是基于“經(jīng)驗性”算法，缺少堅實的數(shù)學(xué)基礎(chǔ)[4]，并且存在模態(tài)混淆等缺陷[5]。Dragomiretskiy等[6]最近提出了一種新的分解方法—變分模態(tài)分解(Variational Mode Decomposition，VMD)。它是一種具有堅實的數(shù)學(xué)基礎(chǔ)的分解方法，通過建立對應(yīng)的約束變分模型表達式，引入二次懲罰因子和拉格朗日乘法算子，將約束性變分問題變?yōu)榉羌s束性變分問題，最終獲得具有數(shù)學(xué)意義并且能有效抑制模態(tài)混疊現(xiàn)象的帶限內(nèi)稟模態(tài)函數(shù)。

VMD相比EMD等方法，在很多方面表現(xiàn)出更加優(yōu)異的性能[5-9]。Zheng等將VMD方法應(yīng)用在碰撞信號的去噪上，效果明顯優(yōu)于EMD方法，文獻[7]論證了VMD方法有很好的濾波特性，文獻[8]等將VMD方法應(yīng)用到旋轉(zhuǎn)機械的碰摩故障診斷，唐貴基等將VMD方法應(yīng)用到軸承的早期故障診斷。但是VMD方法自身也存在著一些缺陷。首先，在用VMD方法分解信號之前，必須預(yù)先設(shè)定模態(tài)數(shù)和懲罰參數(shù)，文獻[7]將去趨勢波動分析(Detrended Fluctuation Analysis，DFA)方法與VMD方法相結(jié)合，建立光纖陀螺信號的長相關(guān)系數(shù)α0與模態(tài)數(shù)k的關(guān)系式，以計算α0獲得最佳的模態(tài)數(shù)，但是忽略了懲罰參數(shù)α對對于VMD分解效果的影響，且僅應(yīng)用于光纖陀螺信號。文獻[8]雖然考慮了懲罰參數(shù)α對于VMD分解效果的影響，但是只能將VMD算法運用到簡單的諧波信號和噪聲的合成信號。其次，在用VMD方法對信號分解后，進行濾波或提取特征值等不同處理時，需要選擇合適的分量和排序，但是VMD方法存在分解出來的分量排列不規(guī)律的缺陷，從而給分量的選擇和排序帶來困難。

針對上述缺陷，本文首先把主模態(tài)分析(Principle Mode Analysis,PMA)與VMD方法結(jié)合，使帶限內(nèi)稟模態(tài)函數(shù)(Band-Limited Intrinsic Mode Function, BIMF)分量按照一定的規(guī)則排序。然后用螢火蟲算法對變分模態(tài)分解方法的最佳影響參數(shù)[k,α]組合進行搜索，并提出一個新的特征值—正交低峰值作為螢火蟲算法優(yōu)化[k,α]的目標。搜索得到的最佳結(jié)果設(shè)定變分模態(tài)分解算法的懲罰參數(shù)α和分量個數(shù)k，并根據(jù)PMA處理預(yù)先設(shè)定的故障特征參數(shù)綜合特征將BIMF分量進行排序。最后，用復(fù)雜調(diào)幅調(diào)頻仿真信號對FA-PMA-VMD方法進行仿真分析，并將FA-PMA-VMD方法運用到實際齒根裂紋故障診斷中，驗證改進VMD方法的實際有效性。

1 FA-PMA-VMD方法

1.1 PMA-VMD方法

式(1)信號其時域波形如圖1所示，由中心頻率分別為200 Hz,400 Hz,20 Hz三個分量合成。

(1)

圖1 輸入信號的時域波形圖

直接用VMD分解后各個BIMF分量及其頻譜圖，如圖2和圖3所示，可觀察到，這些BIMF分量并沒有按照中心頻率高低的順序這類的規(guī)律排列。

當(dāng)實際信號的BIMF分量較多，對于多個排列無規(guī)律的BIMF分量，需要重點研究信號的一個或多個特征時，無法直接判斷哪一個BIMF分量的綜合特征最明顯，也無法直接判斷中心頻率大小相鄰的幾個BIMF分量之間的關(guān)系。因此需要一種針對特征參數(shù)分析的方法對VMD方法分解的BIMF進行分析處理。本文引入PMA方法對VMD方法分解的分量進行排序。PMA是受主分量分析啟發(fā)，根據(jù)故障特征參數(shù)，把主模態(tài)用于振動信號BIMF分量分析的方法。

本文把某種故障特征最明顯的BIMF分量稱為主模態(tài)分量，主模態(tài)分量對應(yīng)的模態(tài)稱為主模態(tài)。PMA方法為：

步驟1故障信號X(t)經(jīng)過信號分解后為個n個BIMF與余量res之和即

(2)

步驟2提取第k個BIMF的n種故障特征參數(shù)MK1,MK2,MK3,…,MKn,將這些故障參數(shù)組成一個故障特征參數(shù)向量FK=[MK1,MK2,MK3,…,MKn](k=1,…,n)。FK可以組成一個故障特征參數(shù)矩陣F

F=(F1,F2,F3,…,FN)T

(3)

m=(1,2,…,n)

(4)

m=(1,2,…,n)

(5)

步驟5式(5)可以看到協(xié)方差矩陣L為方陣，求解L的特征值λ與特征向量v

Lv=λv

(6)

步驟6將特征值按大小排列λ1≥λ2≥…≥λn,每個特征值λk的單位化的特征向量組成的空間V

V=(v1,v2,…,vn)T

(7)

步驟7F(E)到空間V的投影表示為S(S=(s1,s2,…,sn)T)。

(8)

步驟8計算各sk對應(yīng)的故障特征參數(shù)列的投影和Bk

(9)

步驟9按照投影絕對和Bk的大小順序，把Bk對應(yīng)的BIMF分量進行排序，得到對一個信號的BIMF分量進行主分量分析的結(jié)果。

含(E)表示單位化后的變量，故障特征參數(shù)矩陣通過步驟3～步驟6的變換，消除了各個特征參數(shù)向量的相關(guān)性，并得到了一組故障特征參數(shù)矩陣的單位正交基，保證在步驟7中的單位化故障特征參數(shù)投影信息獨立性，步驟8以投影和的方式評判綜合信息量的大小。第k個BIMF分量對應(yīng)的特征參數(shù)向量在空間V中的投影和sk越大，說明該BIMF所包含的綜合故障信息量越大，其對應(yīng)的模態(tài)對故障越敏感。

PMA-VMD方法是通過VMD將輸入信號分解后，把PMA用于對BIMF的分析的方法。具體步驟如下：首先通過VMD方法對信號進行分解，得到BIMF分量；其次提取各BIMF分量的故障特征參數(shù)，得到故障特征參數(shù)向量；然后用PMA對故障特征參數(shù)進行分析，得到綜合故障特征的大小順序；最后按照故障特征的大小順序?qū)IMF重新排序。

通過上面分析可以發(fā)現(xiàn)，PMA-VMD方法對故障特征參數(shù)綜合排序，以達到分解出的BIMF分量按照所選的故障特征參數(shù)綜合特征從強到弱排列。

本文此處以頻率重心為特征參數(shù)，模態(tài)的個數(shù)k依據(jù)仿真信號的情況設(shè)置為3，懲罰參數(shù)α設(shè)置為1 200。

利用PMA-VMD對于式(1)對應(yīng)的輸入信號分解后的結(jié)果如圖4所示。對比圖4與圖2，圖3，可以發(fā)現(xiàn)用PMA-VMD對輸入信號進行分解后BIMF分量按照中心頻率的大小，從小到大排列。

圖2 各個分量的時域波形圖

Fig.2 The time domain waveform of each BIMF signal

圖3 各個分量的頻譜

Fig.3 The frequency spectrum of each BIMF signal

圖4 各個分量的時域波形圖

Fig.4 The time domain waveform of each BIMF signal

1.2 FA-PMA-VMD方法

VMD在對信號分解前，要給出模態(tài)的個數(shù)，和懲罰參數(shù)，即人為設(shè)定k的值和α值。根據(jù)上文可以知道，模態(tài)個數(shù)選取的多少會直接影響最終的分離結(jié)果。

本文以式(1)為輸入信號，來試驗不同[k,α]參數(shù)的選取對于VMD算法的分解效果的影響。

當(dāng)k值取4，同時取一個極大的懲罰參數(shù)值α=3×106。此時中心頻率為200 Hz的調(diào)幅調(diào)頻信號完全被舍棄,如圖5和圖6所示。

當(dāng)k值取4，同時取一個較小的懲罰參數(shù)值α=300，此時三個中心頻率為20 Hz,200 Hz和400 Hz都存在三個BIMF分量中，但是中心頻率為200 Hz的調(diào)幅調(diào)頻信號被分到了BIMF2,BIMF3兩分量中，如圖7和圖8所示。

圖5k=4,α=3×106時域波形圖

Fig.5 The time domain waveform of each BIMF signalk=4,α=3×106

圖6 k=4,α=3×106頻譜

Fig.6 The frequency spectra of each BIMF signalk=4,α=3×106

圖7k=4,α=300時域波形圖

Fig.7 The time domain waveform of each BIMF signalk=4,α=300

圖8 k=4，α=300頻譜

Fig.8 The frequency spectra of each BIMF signalk=4，α=300

由于k值和α的選取對于故障的提取至關(guān)重要。而對一個參數(shù)單獨優(yōu)化往往會忽略這兩個參數(shù)之間的相互影響，只能獲得相對最優(yōu)的結(jié)果。因此，本文引入螢火蟲算法對兩個影響參數(shù)同時進行優(yōu)化，以避免人為主觀因素的干預(yù)，自動篩選出最佳的影響參數(shù)組合。

利用螢火蟲算法搜尋VMD算法的模態(tài)數(shù)和懲罰參數(shù)α?xí)r，需要一個評價參數(shù)來評價螢火蟲每次更新位置的適應(yīng)度。本文引入正交低峰值作為優(yōu)化目標，正交低峰值是一個新的評價參數(shù)，其計算如式(10)所示，式(10)中k為分量的個數(shù)，m為螢火蟲的個數(shù)。

正交低峰值的計算條件是BIMFi,BIMFi+1為中心頻率相鄰的兩個BIMF分量。

(10)

I=Zm=F(k,α)

(11)

β(r)=β0e-γrm

(12)

式中:k為分量的個數(shù)，α為懲罰參數(shù)，這兩個參數(shù)作為函數(shù)的輸入變量，分別為螢火蟲算法中的x軸和y軸上的變量;I為函數(shù)的輸出變量，為螢火蟲算法中z軸的變量。其計算原理如圖9所示。

圖9 目標函數(shù)計算原理圖

式(10)和式(11)對目標函數(shù)的定義可知，在用螢火蟲算法搜尋VMD算法的模態(tài)數(shù)k和懲罰參數(shù)α?xí)r，首先需要以頻率重心為特征參數(shù)，通過PMA對BIMF分量按照中心頻率的大小，從小到大排列，以便求出對應(yīng)分量的個數(shù)和懲罰參數(shù)的目標函數(shù)I。

FA-PMA-VMD算法的計算步驟為：

步驟1設(shè)置故障特征參數(shù)(頻率重心為默認設(shè)置，用于計算螢火蟲目標函數(shù)前通過PMA對獲得的BIMF分量進行排序。其他故障特征參數(shù)根據(jù)實際的故障來設(shè)置)。

步驟2初始化螢火蟲算法參數(shù)，最大吸引力β0、光吸收系數(shù)γ、步長因子α，給定群體規(guī)模m，指定迭代次數(shù)T。

步驟3隨機初始化螢火蟲的位置，計算螢火蟲的目標函數(shù)I=Zm=F(k,α)作為各自最大熒光亮度In。

步驟4計算螢火蟲的相對亮度和吸引度β[10]，吸引度公式如式(12)所示，根據(jù)相對亮度決定螢火蟲的移動。

步驟5更新螢火蟲的空間位置，對處于最佳位置的螢火蟲進行隨機擾動。

步驟6根據(jù)螢火蟲的位置從新計算螢火蟲的亮度。

步驟7判斷是否滿足終止條件，如果滿足，跳出循環(huán)獲得最優(yōu)的全局極值點Zm和個體[k,α]。如不滿足，則重復(fù)步驟3～步驟5直到滿足循環(huán)終止條件。

步驟8根據(jù)前面確定的[k,α]，確定BIMF最終分量的個數(shù)，最后通過PMA獲得按照所選的故障特征參數(shù)綜合特征從強到弱排列BIMF分量。

2 仿真信號分析

在式(1)的信號中加入一個正弦信號n(t)=sin(2π120t)，如式(13)所示，作為輸入信號,如圖10所示。

S(t)=X(t)+n(t)

(13)

對輸入信號S(t)用EEMD進行分解，對于幅值系數(shù)K和總體平均次數(shù)M的選取，Wu等[11]建議K由原信號的標準差乘以一個分數(shù)來定義， 一般取原信號標準差的0.2倍。這樣當(dāng)M為一、兩百次的時候，殘留噪聲引起的誤差一般會處在一個比較低的水平(不足1%)。同時，他們還認為：在噪聲水平合適的條件下，一味地增加執(zhí)行次數(shù)對結(jié)果的改善并不顯著。而反過來在一定范圍內(nèi)提高噪聲水平對結(jié)果的影響也甚微，當(dāng)然噪聲太大也可能完全淹沒信號的本質(zhì)特征，使得分解毫無意義。本文中K取原信號標準差的0.2倍，M取值為100。

圖10 加入諧波信號后的時域波形圖

其結(jié)果如圖11所示?？梢园l(fā)現(xiàn)分量之間出現(xiàn)了較為嚴重的頻率混疊現(xiàn)象。

用FA-PMA-VMD方法對輸入信號進行處理，首先初始化螢火蟲算法的相關(guān)參數(shù)如表1所示。表1中：β0為最大吸引力;γ為光吸收系數(shù);α為步長因子;n為群體的規(guī)模;T為迭代的次數(shù)。

表1 螢火蟲算法的相關(guān)參數(shù)

圖11 EEMD分解時域波形圖

齒輪齒根裂紋故障會對齒輪的信號帶來沖擊，文獻[12]表明峭度(Kurtosis，Kurt)對于沖擊信號特別敏感，而均方根值(Root Mean Square，RMS) 能表明信號的強弱，多尺度模糊熵(Multiscale Fuzzy Entropy,MFE)為信號復(fù)雜度的量化統(tǒng)計指標[13]，因此都對齒根裂紋故障比較敏感。齒根裂紋故障信號為調(diào)幅調(diào)頻信號，故本文將Kurt,RMS和MFE為故障特征參數(shù)。

按照FA-PMA-VMD方法的計算步驟，經(jīng)過計算可以得到全局極點min(Imax)=0.009 2和最優(yōu)個體值[k,α]=[4,2 325],計算結(jié)果如圖12所示。多次計算后的結(jié)果如表2所示。

圖12 螢火蟲算法對k和α優(yōu)化后結(jié)果

次數(shù)12345α2 3252 2972 3142 3312 330K44444

取五次計算后的平均值[4,2 319.4]作為優(yōu)化后的最終結(jié)果，F(xiàn)A-PMA-VMD分解的最終分量的故障特征參數(shù)列的投影和Bk如表3所示。其分解結(jié)果如圖13和圖14所示。

對比圖11和圖13、圖14，可以觀察到FA-PMA-VMD對仿真信號的分解后，在效果優(yōu)于EEMD，同時也可以觀察到EEMD分解后的分量只是由能量的大小按照高頻到低頻排列，而FA-PMA-VMD分解后，將含有豐富故障特征的調(diào)幅調(diào)頻信息的BIMF分量排在前面。

表3最終分量的故障特征參數(shù)列的投影和Bk值

Tab.3Bkvaluesofthefaultcharacteristicparametersofthefinalcomponent

BIMFii=1i=2i=3i=4Bi1.122 30.803 30.778 20.672 1

圖13 FA-PMA-VMD分解分量時域波形圖

圖14 FA-PMA-VMD分解分量頻譜

Fig.14 Frequency spectra of the results decomposed by FA-PMA-VMD

3 實驗信號分析

為驗證本文方法對于振動信號分析的有效性，利用實際含有齒根裂紋的齒輪振動信號進行驗證。

實驗中齒輪為模數(shù)2.5 mm、齒數(shù)37的標準齒輪，齒根切割的裂紋寬度為0.15 mm，深度為20%。實驗伺服電機主軸轉(zhuǎn)速為 600 r/min，即轉(zhuǎn)頻為fr=10 Hz,嚙合頻率fz=370 Hz。實驗采樣頻率為1 024 Hz，每次采樣時間持續(xù)1 s。采集到的實驗信號如圖15所示。

圖15 齒輪裂紋故障下的實驗信號

用FA-PMA-VMD方法對實驗信號進行分解，根據(jù)第2節(jié)仿真信號分析，選取RMS,Kurt,MFE為故障特征參數(shù)。最終分量的個數(shù)以及對應(yīng)分量的故障特征參數(shù)列的投影和Bk如表4所示。分解結(jié)果如圖16和圖17所示。繼續(xù)對BIMF2分量進行包絡(luò)解調(diào)運算，并作出其包絡(luò)譜如圖18所示，其中特征頻率及倍頻處譜線幅值十分突出，說明特征頻率信息被提取了出來并包含在BIMF2分量中。

表4最終分量的故障特征參數(shù)列的投影和Bk值

Tab.4Bkvaluesofthefaultcharacteristicparametersofthefinalcomponent

BIMFii=1i=2i=3i=4i=5i=6Bi0.671 10.328 90.137 20.134 40.100 00.029 3

圖16 BIMF分量的時域波形圖

圖17 BIMF分量的頻譜

圖18 BIMF2分量的包絡(luò)譜

圖19 20組實驗信號中故障信息所在的BIMF分量及其個數(shù)

Fig.19 The BIMF components and the numbers of the fault information in the 20 experimental signals

為了進一步驗證FA-PMA-VMD方法可以將含豐富故障信息的BIMF分量靠前排列，選取與實驗信號同一測點但不同時間段獲取的20組信號。用FA-PMA-VMD方法對這些信號進行分解。然后以故障信息所在的分量的排列順序號為x軸，x對應(yīng)的信號的個數(shù)(含有故障信息的同一順序號的BIMF分量個數(shù))為y軸建立直方圖如圖19所示?？梢杂^察到這20組試驗信號經(jīng)過FA-PMA-VMD方法分解后獲得的20組BIMF分量中，故障信息所在的分量主要集中在每組分量的前三項，表明FA-PMA-VMD方法可以將含豐富故障信息的BIMF分量靠前排列。

4 結(jié) 論

(1) 本文針對VMD方法分解所得BIMF分量排序不規(guī)律的問題，提出了PMA-VMD的方法，把主模態(tài)分析與VMD方法結(jié)合起來，根據(jù)故障特征來對BIMF分量進行排序，可將含豐富故障信息的BIMF分量靠前排列。對于大量BIMF分量組中含故障信息的BIMF分量的甄別具有指導(dǎo)意義。

(2) 針對PMA-VMD中[k,α]參數(shù)的選取對分解結(jié)果影響較大的問題，引入螢火蟲算法，提出FA-PMA-VMD方法對參數(shù)[k,α]進行組合優(yōu)化。并在仿真信號中證明了此方法的可行性。

(3) 用FA-PMA-VMD方法對與帶齒根裂紋的齒輪振動實際信號進行分析，提取了有效的故障特征，得到了準確的診斷結(jié)果。具有一定的實用性。