謝棕，張麗萍

(福州大學(xué) 機(jī)械工程及自動(dòng)化學(xué)院，福建 福州 350108)

0 引言

齒輪及齒輪箱作為機(jī)械設(shè)備常用的調(diào)節(jié)轉(zhuǎn)速和傳遞轉(zhuǎn)矩的旋轉(zhuǎn)機(jī)械設(shè)備，不僅能夠傳遞較大的功率和載荷，而且具有較好的可靠性[1]。但是在高精度的切削加工中，當(dāng)齒輪在變轉(zhuǎn)速、變載荷等復(fù)雜工況下工作，極易受到損傷、產(chǎn)生磨損、斷齒等情況，使得加工精度大打折扣。因此，齒輪的狀態(tài)監(jiān)測和故障診斷變得尤為重要。

齒輪的故障信號(hào)多為非線性、非平穩(wěn)信號(hào)。經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)和集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)在分析非線性、非平穩(wěn)信號(hào)時(shí)具有良好的優(yōu)越性，但是存在嚴(yán)重的端點(diǎn)效應(yīng)和模態(tài)混疊問題[2]。基于此，DRAGOMIRETSKIY K[3]在2014年提出一種自適應(yīng)、非遞歸模式的信號(hào)分解方法——變分模態(tài)分解(variational mode decomposition,VMD)，它很好地解決了EMD算法存在的問題，并且具有更堅(jiān)實(shí)的數(shù)學(xué)理論基礎(chǔ)。但VMD算法在分解過程中會(huì)同時(shí)受分解參數(shù)K和懲罰因子α這兩個(gè)參數(shù)的影響。因此，如何快速、準(zhǔn)確地確定這兩個(gè)參數(shù)是研究的重點(diǎn)。劉長良等[4]提出一種觀察中心頻率的方法確定分解參數(shù)K值，當(dāng)開始出現(xiàn)中心頻率相近時(shí)，認(rèn)為出現(xiàn)了過分解，以此來確定K值，但卻忽略了懲罰因子α對(duì)分解結(jié)果的影響。瞿紅春等[5]利用手動(dòng)尋優(yōu)對(duì)VMD參數(shù)進(jìn)行優(yōu)化，先假定α不變，不斷尋優(yōu)K值使得包絡(luò)熵最小，再將此K值設(shè)為定值，不斷尋優(yōu)α值使得包絡(luò)熵最小，得出分解參數(shù)K和懲罰因子α。但這樣得出的K和α并不是最優(yōu)值，且工作量巨大?；谝陨戏治觯疚奶岢鲆訧MF分量局部包絡(luò)熵最小為優(yōu)化目標(biāo)，利用改進(jìn)的蝙蝠算法優(yōu)化VMD分解參數(shù)K和懲罰因子α，進(jìn)而得到最優(yōu)的分解參數(shù)，通過計(jì)算K個(gè)IMF分量的多尺度排列熵(multiscale permutation entropy,MPE)構(gòu)造故障特征向量，輸入到極限學(xué)習(xí)機(jī)(extreme learning machine,ELM)進(jìn)行訓(xùn)練和識(shí)別，進(jìn)而實(shí)現(xiàn)齒輪的故障診斷。

1 參數(shù)優(yōu)化VMD

1.1 VMD原理

VMD分解是一種非遞歸模式的信號(hào)分解方法，通過迭代搜尋約束變分模型最優(yōu)解，將輸入信號(hào)f分解成k個(gè)離散的模態(tài)uk(t)，使得每個(gè)模態(tài)的估計(jì)帶寬之和最小，約束條件為各模態(tài)之和等于輸入信號(hào)f。

構(gòu)造相應(yīng)的約束變分模型如下：

其中：{uk}={u1,…,uk}，表示分解得到的k個(gè)模態(tài)函數(shù)；{wk}={w1,…,wk}，表示各模態(tài)分量的中心頻率[6]。

為求解式(1)的最優(yōu)解，引入拉格朗日乘法算子λ(t)以及二次懲罰因子α，將約束性變分問題變?yōu)榉羌s束性變分問題。λ(t)讓約束條件保持嚴(yán)格性，α能充分減少高斯噪聲對(duì)信號(hào)的干擾：

(2)

利用交替方向乘子算法不斷迭代更新各個(gè)模態(tài)及其中心頻率[7]，尋求式(2)的靶點(diǎn)。具體實(shí)現(xiàn)算法如下(n為迭代次數(shù))：

2)迭代k=1∶K并更新各模態(tài)信號(hào)uk和中心頻率ωk，i∈[1,K]。

更新uk：

(3)

更新ωk：

(4)

3)根據(jù)下式更新λ，Γ為更新因子

(5)

4)若滿足下方收斂條件則停止迭代，否則重復(fù)1)、2)，ε為>0的正數(shù)。

(6)

結(jié)束迭代，得到K個(gè)IMF分量。

1.2 改進(jìn)蝙蝠算法

蝙蝠算法(bat algorith,BA)是由英國學(xué)者YANG X S[8]在2010年提出的一種群智能優(yōu)化算法。BA算法把求解問題目標(biāo)函數(shù)的適應(yīng)度值度量成蝙蝠個(gè)體所處位置的優(yōu)劣。最優(yōu)解類比蝙蝠的獵物，脈沖響度和頻率的變化在一定程度上表明與最優(yōu)解的靠近程度。

在BA算法中，搜索頻率、位置和速度的更新公式如下：

fi=fmin+(fmax-fmin)β

(7)

(8)

(9)

其中：fi、fmax、fmin分別表示第i只蝙蝠當(dāng)前時(shí)刻發(fā)出的聲波頻率以及聲波頻率的最大值與最小值；β∈[0,1]的隨機(jī)數(shù)；x*表示當(dāng)前最優(yōu)位置。

在局部搜索中，如果一個(gè)蝙蝠選擇了一個(gè)最優(yōu)解，那么將在此最優(yōu)解附近隨機(jī)產(chǎn)生一個(gè)新解：

Xnew=Xold+εAt

(10)

其中：Xold為從當(dāng)前最優(yōu)解中選擇的一個(gè)解；At為t時(shí)刻所有蝙蝠響度的平均值；ε是[-1,1]的一個(gè)隨機(jī)數(shù)。

蝙蝠脈沖的響度Ai和速率ri通過下式計(jì)算：

(11)

(12)

其中：α和γ是常數(shù)，并且0<α<1，γ>0。響度Ai和速率ri的持續(xù)迭代更新表明與最優(yōu)解的靠近程度。

蝙蝠算法能夠在算法運(yùn)行的前期通過將全局優(yōu)化轉(zhuǎn)換到局部優(yōu)化來實(shí)現(xiàn)算法的快速收斂，這同時(shí)也會(huì)導(dǎo)致該算法過早地處于停滯階段。為改善這樣的問題，本文在速度更新公式(8)中引入了線性遞減慣性權(quán)重ω(k)的概念：

(13)

ω(k)=wmax(wmax-wmin)(Tmax-k)/Tmax

(14)

其中：ωmax為初始慣性權(quán)重；ωmin為迭代至最大次數(shù)時(shí)的慣性權(quán)重；k為當(dāng)前迭代次數(shù)；Tmax為最大迭代次數(shù)。查閱相關(guān)文獻(xiàn)知，ωmax=0.9，ωmin=0.4時(shí)算法性能最好[9]。迭代初期較大的初始慣性權(quán)重使算法保持了較強(qiáng)的全局搜索能力，而迭代后期較小的慣性權(quán)重有利于算法進(jìn)行更精確的局部搜索。

1.3 VMD參數(shù)優(yōu)化

VMD分解前，需要預(yù)先設(shè)定分量個(gè)數(shù)K和懲罰因子α的值。K和α的取值不同，會(huì)有不同的分解結(jié)果。K值過小，會(huì)造成模態(tài)混疊；K值過大，會(huì)產(chǎn)生虛假分量。α值過小，各IMF分量的帶寬越大；α值過大，各IMF分量的帶寬越小。因此，如何確定合適的K和α的值是研究的關(guān)鍵所在。

利用蝙蝠算法對(duì)VMD進(jìn)行參數(shù)優(yōu)化時(shí)，需要選定一個(gè)適應(yīng)度函數(shù)。本文采用VMD分解后的局部極小包絡(luò)熵作為蝙蝠算法的適應(yīng)度函數(shù)，搜尋最優(yōu)的VMD參數(shù)組合。信號(hào)pj的包絡(luò)熵Ep可以表示成：

(15)

(16)

(17)

式中：pj是a(j)的概率分布序列；a(j)為信號(hào)x(j)經(jīng)過Hilbert解調(diào)后得到的包絡(luò)信號(hào)。

改進(jìn)BA算法優(yōu)化VMD參數(shù)步驟如下：

步驟1 BA算法基本參數(shù)初始化。主要有種群大小sizepop、最大迭代次數(shù)maxgen、參數(shù)α和K的搜索范圍、最大脈沖響度A0、最大脈沖速率r0、脈沖頻率范圍、音量的衰減系數(shù)α以及搜索頻率的增強(qiáng)系數(shù)γ。

步驟2 初始化種群，以局部極小包絡(luò)熵作為適應(yīng)度函數(shù)求出一組最優(yōu)的α和K。

步驟3 根據(jù)式(7)、式(13)、式(9)不斷迭代更新種群的搜索脈沖頻率、速度和位置。

步驟4 生成均勻分布隨機(jī)數(shù)rand，如果rand>ri，則對(duì)當(dāng)前最優(yōu)解進(jìn)行隨機(jī)擾動(dòng)，產(chǎn)生一個(gè)新的解。

步驟5 如果計(jì)算出來的包絡(luò)熵小于所得出的極小包絡(luò)熵，對(duì)最優(yōu)的α和K以及最小的包絡(luò)熵進(jìn)行更新，然后按式(11)、式(12)對(duì)響度和頻率進(jìn)行更新。

步驟6 重復(fù)步驟3-步驟5，直至達(dá)到最大迭代次數(shù)以及確定最小的包絡(luò)熵值，得出最佳參數(shù)α和K。

2 多尺度排列熵

排列熵是由BANDT C等[10]提出的一種衡量一維時(shí)間序列復(fù)雜度的平均熵參數(shù)，常用來提取機(jī)械故障的特征。但齒輪及齒輪箱運(yùn)轉(zhuǎn)過程中的故障特征信息分布在多尺度中，僅用單一尺度的排列熵進(jìn)行分析，會(huì)遺漏其余尺度上的故障特征信息。

針對(duì)排列熵的不足，COSTA M等[11]提出了多尺度排列熵的概念。多尺度排列熵是在排列熵的基礎(chǔ)上將時(shí)間序列進(jìn)行多尺度粗?；缓笥?jì)算不同尺度下粗?；蛄械呐帕徐兀唧w計(jì)算過程如下：

(18)

式中：s為尺度因子；[N/s]表示取整。

(19)

式中：l表示第l個(gè)重構(gòu)的分量；τ表示延遲時(shí)間；m表示嵌入維數(shù)。

3)將式(19)得到的重構(gòu)時(shí)間序列按升序排列可得

S(g)=(j1,j2,…,jm)

(20)

(21)

4)當(dāng)Pg=1/m!時(shí)，每種符號(hào)序列的概率都相等，此時(shí)時(shí)間序列的復(fù)雜程度最高，排列熵最大，為lnm!。為方便表示，通常會(huì)將Hp(m)進(jìn)行歸一化處理，即

Hp=Hp(m)/ln(m!)

(22)

Hp的取值范圍為[0,1]，Hp的大小反映了時(shí)間序列的復(fù)雜和隨機(jī)程度。Hp越大，說明時(shí)間序列越隨機(jī)，反之，說明時(shí)間序列越規(guī)則。

3 算法流程

本文提出基于參數(shù)優(yōu)化VMD和多尺度排列熵的齒輪故障診斷方法，具體步驟如下：

1)利用改進(jìn)的蝙蝠算法優(yōu)化VMD分解參數(shù)K和懲罰因子α，進(jìn)而得到最優(yōu)的分解參數(shù)，得到K個(gè)模態(tài)分量IMF1、IMF2、…、IMFK。

2)計(jì)算被選模態(tài)分量的多尺度排列熵，選取合適尺度的排列熵值作為故障特征向量。

3)將故障特征向量輸入到已經(jīng)訓(xùn)練好的極限學(xué)習(xí)機(jī)中進(jìn)行訓(xùn)練和識(shí)別。

故障診斷流程如圖1所示。

圖1 齒輪故障診斷流程

4 實(shí)驗(yàn)分析

4.1 實(shí)驗(yàn)說明

為了驗(yàn)證本文所提方法的有效性，在實(shí)驗(yàn)平臺(tái)下(圖2)，采集試驗(yàn)臺(tái)減速箱齒輪不同狀態(tài)下的實(shí)時(shí)振動(dòng)信號(hào)。齒輪減速箱型號(hào)為JZQ200，采樣頻率為8 200 Hz。選擇齒輪正常、不平衡、松動(dòng)、錯(cuò)位、斷齒、裂痕這6種不同狀態(tài)的數(shù)據(jù)進(jìn)行實(shí)驗(yàn)驗(yàn)證。

圖2 齒輪實(shí)驗(yàn)平臺(tái)

4.2 故障特征提取

以減速箱齒輪出現(xiàn)裂痕為分析對(duì)象，本文在以下的實(shí)驗(yàn)中選取的蝙蝠算法參數(shù)均為：sizepop=36，參數(shù)α和K的搜索范圍low=[1 500 3]，high=[2 500 8]，速度v的范圍為[-100 100]，A0=0.5，r0=0.25，脈沖頻率范圍[0 1]，α=0.8，γ=0.5。蝙蝠算法不同參數(shù)組合的優(yōu)化如表1所示。

表1 蝙蝠算法不同參數(shù)組合的優(yōu)化

圖3為局部極小包絡(luò)熵值隨尋優(yōu)迭代次數(shù)變化曲線，優(yōu)化結(jié)果如表2所示。

圖3 尋優(yōu)迭代圖

表2 優(yōu)化結(jié)果

作為對(duì)比，同時(shí)采用手動(dòng)尋優(yōu)的方法來求解最優(yōu)的(K,α)的值。VMD參數(shù)設(shè)置如下：初始化alpha=2 000，tau=0，DC=0，init=1，tol=le-7。

表3 不同K值下的包絡(luò)熵值

由于篇幅原因，上述表格每個(gè)K值僅陳列5個(gè)極小包絡(luò)熵值，觀察到最小的包絡(luò)熵值僅為5.554 2。對(duì)比蝙蝠算法的迭代尋優(yōu)結(jié)果，通過手動(dòng)尋優(yōu)求得的最小包絡(luò)熵值并不理想。

綜合分析，改進(jìn)蝙蝠算法迭代尋優(yōu)的最優(yōu)參數(shù)K=6，α=247 4。選取裂痕齒輪一組數(shù)據(jù)中的4 500個(gè)樣本點(diǎn)進(jìn)行VMD分解，得到6個(gè)IMF分量，分解結(jié)果如圖4所示。

圖4 VMD分解結(jié)果

MPE的計(jì)算需要設(shè)置嵌入維數(shù)m，延遲時(shí)間τ和尺度因子s。其中維數(shù)m取值范圍通常是3～7，m取值太小，算法的檢測性能降低，若m取值太大，將無法反映時(shí)間序列的細(xì)微變化。延遲時(shí)間τ對(duì)時(shí)間序列的計(jì)算影響較小，尺度因子s最大值一般取>10即可，本文選取m=3，τ=1，s=20。

于是就能求出前面VMD分解后6個(gè)IMF分量的排列熵，如圖5所示。

圖5 裂痕齒輪各IMF分量的MPE值

4.3 齒輪的故障診斷

對(duì)測得的正常、不平衡、松動(dòng)、錯(cuò)位、斷齒、裂痕這6種不同狀態(tài)的數(shù)據(jù)，每種數(shù)據(jù)取50組，每組4 500個(gè)點(diǎn)。對(duì)這50組數(shù)據(jù)進(jìn)行VMD分解后并計(jì)算模態(tài)分量的MPE值，6個(gè)MPE分量，每個(gè)分量提取20尺度的因子，所以每個(gè)樣本就得到維度為120的特征。將故障特征向量輸入到已經(jīng)訓(xùn)練好的ELM中，取前40組數(shù)據(jù)作為已知故障樣本，后10組數(shù)據(jù)作為待識(shí)別故障樣本，診斷情況如圖6所示。

圖6 ELM故障識(shí)別

圖6中，橫坐標(biāo)分別代表訓(xùn)練樣本和測試樣本，縱坐標(biāo)代表6種故障類別的標(biāo)簽。測試樣本中的空心點(diǎn)代表預(yù)測類別。實(shí)心點(diǎn)代表實(shí)際類別。當(dāng)測試的數(shù)據(jù)符合同類故障的標(biāo)準(zhǔn)時(shí)，空心點(diǎn)和實(shí)心點(diǎn)會(huì)重合，代表識(shí)別正確。本文測試集總的準(zhǔn)確率達(dá)到了96.7%，對(duì)于絕大多數(shù)的故障情況，都能達(dá)到很好的識(shí)別效果。

5 結(jié)語

針對(duì)齒輪故障特征在單一尺度難以全面提取的問題，提出一種基于參數(shù)優(yōu)化的VMD和多尺度排列熵的齒輪故障診斷方法。相比較手動(dòng)尋優(yōu)，該方法可以搜尋更優(yōu)的VMD參數(shù)組合，更有效地提取出不同故障狀態(tài)下的信號(hào)特征參數(shù)。分析結(jié)果表明，通過參數(shù)優(yōu)化VMD和多尺度排列熵的齒輪故障診斷方法，極大地提高了齒輪故障診斷的準(zhǔn)確性。