張旭幫，王靈梅，孟恩隆，賈成真，原升耀，王 平

(1.山西大學(xué)自動(dòng)化與軟件學(xué)院，山西 太原 030013)(2.煤炭工業(yè)太原設(shè)計(jì)研究院，山西 太原 030001)

風(fēng)電機(jī)組通常在復(fù)雜多變的惡劣環(huán)境中運(yùn)行，傳動(dòng)鏈中的齒輪、軸承等部件容易出現(xiàn)故障，影響正常的能量傳遞，造成發(fā)電損失，嚴(yán)重時(shí)會(huì)損壞齒輪箱、發(fā)電機(jī)等關(guān)鍵部件，形成巨大的運(yùn)維成本。因此，行之有效的故障診斷方法對于確保風(fēng)電機(jī)組的可靠和安全運(yùn)行至關(guān)重要。由于機(jī)組運(yùn)行受環(huán)境噪聲影響，信號表現(xiàn)為非平穩(wěn)、非線性，而目前的時(shí)域、頻域分析方法效果不佳，因此時(shí)頻分析法被提出。Huang等[1]提出了經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition，EMD)方法，在非平穩(wěn)信號處理方面得到了廣泛應(yīng)用，但它依賴基函數(shù)的選取，在遞歸過程中誤差會(huì)不斷疊加，存在端點(diǎn)效應(yīng)、虛假分量、模態(tài)混疊、缺乏理論支撐等問題；Dragomiretskiy等[2]提出了變分模態(tài)分解(variational mode decomposition，VMD)方法，其本質(zhì)上是通過變分問題使分解結(jié)果穩(wěn)定，雖然克服了模態(tài)混疊、端點(diǎn)效應(yīng)等缺點(diǎn)，但VMD方法受預(yù)設(shè)參數(shù)模態(tài)分量數(shù)K與懲罰因子α的制約。

近年來，王朝閣等[3]提出利用粒子群算法對最優(yōu)最小反褶積算法參數(shù)進(jìn)行優(yōu)化，確定濾波器長度與故障周期搜索范圍，對軸承故障進(jìn)行識(shí)別與診斷；夏俊等[4]利用松鼠優(yōu)化算法優(yōu)化共振稀疏分解的品質(zhì)因子與權(quán)重系數(shù)，得到高低共振分量，經(jīng)Hilbert包絡(luò)分析后用于提取軸承微弱故障信息有一定的效果；王奉濤等[5]提出在VMD算法的不同K值下進(jìn)行信號分解，通過計(jì)算各分量與信號總能量，選取最佳本征模態(tài)分量(intrinsic mode function，IMF)，避免了信號分解不足或過分解問題；李華等[6]為使VMD算法取得較好的分解效果，對比不同K值下分解所得到的IMF峭度值，確定最優(yōu)K值，最后進(jìn)行濾波、包絡(luò)解調(diào)分析判斷軸承故障。在上述研究中，僅考慮了模態(tài)分量數(shù)K值對VMD算法的影響，而忽略了懲罰因子的影響。李華等[7]提出了基于信息熵優(yōu)化VMD參數(shù)的方法，但優(yōu)化過程是順序優(yōu)化，計(jì)算不同K值下所得分量的信息熵，確定最小信息熵對應(yīng)K值，通過既定K值去確定懲罰因子，因此所得參數(shù)組合不一定是全局最優(yōu)解。

針對VMD優(yōu)化參數(shù)難以確定的問題，本文提出了基于改進(jìn)參數(shù)優(yōu)化VMD的故障特征提取方法?？紤]到旋轉(zhuǎn)機(jī)械發(fā)生故障時(shí)所具有的強(qiáng)噪聲、周期性特征，用對周期性沖擊較敏感的峭度指標(biāo)和衡量信號時(shí)間序列復(fù)雜度的樣本熵構(gòu)造綜合性能指標(biāo)，作為優(yōu)化算法的適應(yīng)度函數(shù)。將振動(dòng)信號經(jīng)參數(shù)優(yōu)化的變分模態(tài)分解后，篩選綜合性能指標(biāo)較小的模態(tài)分量進(jìn)行重構(gòu)，此時(shí)所得到的重構(gòu)信號既保留了原始信號的特征，也突出了故障周期沖擊特征，最后進(jìn)行Hilbert包絡(luò)分析，判斷故障部位。通過對山西某風(fēng)電場的齒輪箱故障數(shù)據(jù)與風(fēng)電機(jī)組傳動(dòng)系統(tǒng)故障試驗(yàn)臺(tái)采集的實(shí)測故障數(shù)據(jù)進(jìn)行對比驗(yàn)證，表明該方法能夠在強(qiáng)噪聲背景下準(zhǔn)確提取故障特征信息。

1 基本原理

1.1 變分模態(tài)分解

VMD算法假設(shè)實(shí)際信號分解所得IMF是具有中心頻率的有限帶寬，通過循環(huán)迭代確定最優(yōu)解，尋找K個(gè)模態(tài)分量，經(jīng)過計(jì)算信號梯度平方的L-2范數(shù)來估計(jì)各分量的帶寬，使每個(gè)模態(tài)的估計(jì)帶寬之和最小，實(shí)現(xiàn)對信號的自適應(yīng)分解。其對應(yīng)的約束變分模型表達(dá)式為：

(1)

引入二次懲罰因子α和Lagrange乘子λ，利用增廣拉格朗日乘子法將式(1)轉(zhuǎn)化為無約束優(yōu)化模型，表達(dá)式如下：

(2)

式中：L(·)表示利用增廣拉格朗日乘數(shù)法構(gòu)建的無約束方程；y(t)為式(1)中的所有模態(tài)之和；λ,λ(t)為拉格朗日乘子。

(3)

(4)

(5)

具體實(shí)現(xiàn)過程如下：

3)根據(jù)式(5)更新λ。

4)根據(jù)式(6)判斷是否收斂，其中ε為任意大于0的正數(shù)，代表分解精度。若滿足收斂條件則信號分解過程結(jié)束；否則迭代次數(shù)n=n+1,返回步驟2)繼續(xù)分解。

(6)

1.2 樣本熵

為了減小環(huán)境噪聲對提取故障特征的影響，本文引入樣本熵對所采集的實(shí)際信號進(jìn)行評價(jià)。樣本熵由非線性動(dòng)力學(xué)研究人員在20世紀(jì)末提出[8]，可以衡量實(shí)際信號的隨機(jī)性和復(fù)雜性。

假設(shè)實(shí)際信號是由N個(gè)數(shù)據(jù)組成的時(shí)間序列{x(n)}=x(1),x(2),…,x(N)，則序列樣本熵SampEn(m,r)為：

(7)

式中：r為相似容限閾值；m為初始分段的維數(shù)；Bm(r)和Am+1(r)分別是由時(shí)間序列組成的m維與(m+1)維向量序列，經(jīng)過相空間重構(gòu)、重構(gòu)向量距離計(jì)算，滿足相似容限r(nóng)下匹配m與(m+1)個(gè)數(shù)據(jù)點(diǎn)的概率。

當(dāng)N為有限值時(shí)，實(shí)際信號經(jīng)VMD算法分解后模態(tài)分量IMFi的樣本熵SampEn(m,r,IMFi)可用式(8)估計(jì)：

(8)

式中：i為信號分解得到的模態(tài)分量序號；r一般取輸入時(shí)間序列IMFi標(biāo)準(zhǔn)差的0.2倍。若模態(tài)分量的自相似性越高，信號復(fù)雜度越低，包含的周期性故障信息越多，則樣本熵值越小；反之，模態(tài)分量的自相似性較低，信號復(fù)雜度越高，包含周期性特征越少，則樣本熵值越大。

2 故障特征提取流程

由于故障信息往往包含在由故障沖擊成分引起的幅值調(diào)制信號中，為使所分解的模態(tài)分量包含較多的故障特征信息，本文所選目標(biāo)函數(shù)在樣本熵的基礎(chǔ)上引入了對早期沖擊故障較為敏感的峭度指標(biāo)，各模態(tài)分量的峭度qdi定義為：

(9)

式中:E(·)表示括號中變量值的數(shù)學(xué)期望；μi，σi分別為IMFi的均值與方差；IMFi(j)為第i個(gè)分量中的第j個(gè)元素。故綜合性能指標(biāo)C作為遺傳算法的目標(biāo)函數(shù)可表示為：

C=SampEn(m,r,IMFi)+1/abs(qdi-3)

(10)

具體優(yōu)化步驟如下：

步驟1，種群初始化。種群中的個(gè)體由模態(tài)分解個(gè)數(shù)K和懲罰因子α組成，不同個(gè)體的[K,α]組合不同。

步驟2，對初始種群進(jìn)行選擇、交叉和變異。按照目標(biāo)函數(shù)式(10)對種群中每個(gè)個(gè)體的目標(biāo)函數(shù)值進(jìn)行求解，得到目標(biāo)函數(shù)最小值Cmin所對應(yīng)的個(gè)體，并記錄每一代個(gè)體進(jìn)化的最佳適應(yīng)度與平均適應(yīng)度。

步驟3，根據(jù)[K,α]組合的最優(yōu)值，由VMD算法對振動(dòng)信號進(jìn)行分解，選取各模態(tài)分量與原始信號相關(guān)性較大且綜合性能指標(biāo)較小的分量進(jìn)行Hilbert包絡(luò)譜分析，得出診斷結(jié)果。

基于以上理論，本文所提出的基于改進(jìn)參數(shù)優(yōu)化VMD的故障特征提取流程如圖1所示。

圖1 改進(jìn)參數(shù)優(yōu)化VMD故障特征提取流程

3 仿真信號驗(yàn)證

為驗(yàn)證本文所提方法的有效性，對齒輪箱行星輪故障仿真信號進(jìn)行特征提取，不失一般性，本文不考慮諧波因素，只考慮載波頻率和調(diào)制頻率的基頻與振動(dòng)傳遞路徑對振動(dòng)信號的影響[9]，故齒輪箱行星輪故障可用如式(11)的調(diào)幅調(diào)頻信號x(t)表示[10]。

x(t)=[1-cos(2πfrt)][1+Acos(2πfst)]·cos[2πfmt+Bsin(2πfst+φ)]+n(t)

(11)

式中：fr=30 Hz，為行星架的旋轉(zhuǎn)頻率；fs=75 Hz，為行星輪的故障特征頻率；fm=1 000 Hz，為齒輪嚙合頻率；A,B為信號的調(diào)制系數(shù)，均為1；初始相位φ=0；n(t)為高斯白噪聲。

為使仿真信號更接近實(shí)際現(xiàn)場工況，仿真信號中加入信噪比為-12 dB的高斯白噪聲n(t)，采樣頻率為12 kHz，分析數(shù)據(jù)點(diǎn)數(shù)為8 196個(gè)。齒輪箱行星輪故障仿真信號分析結(jié)果如圖2所示。

圖2 仿真信號分析

圖2(a)為Fourier頻譜圖，由圖可以看出，嚙合頻率fm、嚙合頻率與行星輪故障頻率的組合(fm+fs)以及在其兩側(cè)出現(xiàn)的以行星架轉(zhuǎn)頻的邊頻帶占據(jù)主導(dǎo)地位，符合行星輪發(fā)生故障時(shí)具有的特征；圖2(b)為包絡(luò)譜圖，圖中并未出現(xiàn)任何特征頻率成分，特征信息已被噪聲淹沒；圖2(c)為對VMD算法的參數(shù)尋優(yōu)的過程，設(shè)置尋優(yōu)范圍K∈[3,10]，α∈[500,2 000]，得到模態(tài)分量與懲罰因子的最優(yōu)組合為[8,1 397]。對VMD所得到的各模態(tài)分量，分別計(jì)算加入噪聲信號的相關(guān)性系數(shù)與綜合性能指標(biāo)C，結(jié)果見表1。從表中可以看出，IMF2分量與原始信號的相關(guān)性系數(shù)最大，并且綜合性能指標(biāo)最小，故將IMF2分量作為重構(gòu)信號進(jìn)行分析，結(jié)果如圖3所示。

表1 各模態(tài)分量指標(biāo)計(jì)算

圖3 重構(gòu)信號分析

圖3(a)與圖2(a)對比可以看出，除了有圖2(a)中的主導(dǎo)頻率以外，出現(xiàn)了額外的特征組合頻率(fm+2fs+fr)和(fm+2fs-fr)，并對信號中的噪聲有一定抑制作用，提高了信號信噪比。圖3(b)與圖2(b)對比可以看出，由于噪聲影響，圖2(b)中并未出現(xiàn)任何故障特征信息，而圖3(b)中可以清晰地看出行星架的旋轉(zhuǎn)頻率fr、行星輪的故障特征頻率fs以及(fs+fr)和(fs-fr)的組合頻率，驗(yàn)證了本文所提方法對故障特征提取的有效性。

4 實(shí)測數(shù)據(jù)驗(yàn)證

4.1 試驗(yàn)臺(tái)驗(yàn)證

如圖4所示，采用實(shí)驗(yàn)室風(fēng)電機(jī)組傳動(dòng)系統(tǒng)故障試驗(yàn)臺(tái)驗(yàn)證本文所提方法的有效性。在齒輪箱箱體及所在軸安裝美國PCB公司生產(chǎn)的加速度傳感器采集振動(dòng)信號，在高、低速軸分別配備光電編碼器采集轉(zhuǎn)速信號，并由變頻器控制輸出轉(zhuǎn)速，利用自主研制的數(shù)據(jù)采集儀采集高速軸振動(dòng)信號，其中試驗(yàn)臺(tái)參數(shù)見表2。

本文在轉(zhuǎn)速為1 200 rad/min、采樣頻率為12.8 kHz下采集1 s高速軸小齒輪斷齒的振動(dòng)加速度，由試驗(yàn)臺(tái)參數(shù)可計(jì)算得到高速軸小齒輪轉(zhuǎn)頻fr=1 200/60=20 Hz，齒輪嚙合頻率fm=20×18=360 Hz。采用遺傳算法尋優(yōu)得到模態(tài)分量與懲罰因子的最優(yōu)組合為[8,1 668]，對VMD分解算法得到的模態(tài)分量進(jìn)行選取、重構(gòu)，最后進(jìn)行包絡(luò)譜分析，結(jié)果如圖5所示。

圖4 風(fēng)電機(jī)組傳動(dòng)系統(tǒng)故障試驗(yàn)臺(tái)

表2 試驗(yàn)臺(tái)參數(shù)

圖5 齒輪斷齒數(shù)據(jù)分析

圖5(a)是齒輪斷齒狀態(tài)下未經(jīng)本文方法優(yōu)化重構(gòu)的信號包絡(luò)譜，從圖中可以看出，譜圖中并未出現(xiàn)斷齒齒輪所在軸的轉(zhuǎn)頻與嚙合頻率。經(jīng)本文方法處理的包絡(luò)譜如圖5(b)所示，從圖中可以清晰地看到斷齒齒輪的2倍轉(zhuǎn)頻、6倍轉(zhuǎn)頻占主導(dǎo)地位，除此之外還有很多的倍轉(zhuǎn)頻，因此可以據(jù)此判斷高速軸齒輪發(fā)生了故障。

4.2 工程實(shí)例驗(yàn)證

為驗(yàn)證本文所提方法在實(shí)際工程中的有效性，利用筆者所在團(tuán)隊(duì)自主研發(fā)并已實(shí)際應(yīng)用于山西某風(fēng)電場的狀態(tài)監(jiān)測與故障診斷系統(tǒng)，在齒輪箱表面安裝加速度傳感器，采集風(fēng)電機(jī)組傳動(dòng)鏈振動(dòng)數(shù)據(jù)。已知該風(fēng)電場2#風(fēng)機(jī)在2017年3月初因齒輪箱故障停機(jī)，本文選取的風(fēng)電場監(jiān)視控制與數(shù)據(jù)采集(supervisory control and data acquisition，SCADA)系統(tǒng)實(shí)測數(shù)據(jù)為故障前4個(gè)月中每月某時(shí)刻的數(shù)據(jù)，見表3，對其進(jìn)行分析，根據(jù)齒輪箱結(jié)構(gòu)參數(shù)和高速軸轉(zhuǎn)速，計(jì)算得到齒輪箱二級平行級齒輪的嚙合頻率為fm=Z2fz/60，其中Z2和fz分別為二級平行級高速軸小齒輪齒數(shù)和發(fā)電機(jī)轉(zhuǎn)速。機(jī)組額定轉(zhuǎn)速為1 800 rad/min，齒輪箱由行星級、一級平行級、二級平行級組成，其參數(shù)見表4。

表3 風(fēng)電場實(shí)測數(shù)據(jù)

對該風(fēng)電場4個(gè)時(shí)間點(diǎn)的齒輪箱高速軸齒輪振動(dòng)加速度信號進(jìn)行分析，結(jié)果如圖6所示。經(jīng)過分析發(fā)現(xiàn)，在圖6(a)中振動(dòng)信號包絡(luò)譜較為整齊，且幅值較小，并未出現(xiàn)額定轉(zhuǎn)速下的高速軸嚙合頻率級轉(zhuǎn)頻；而圖6(b)、圖6(c)、圖6(d)均出現(xiàn)對應(yīng)時(shí)間點(diǎn)數(shù)據(jù)下齒輪嚙合頻率的1、2、3倍頻，且占據(jù)主導(dǎo)地位，包絡(luò)譜幅值較2016年11月22日劇增，說明在2016年11月22日至12月1日之間二級平行級齒輪處出現(xiàn)故障，直至停機(jī)，與現(xiàn)場的實(shí)際情況相符。

表4 齒輪箱參數(shù)

圖6 風(fēng)電場實(shí)測故障數(shù)據(jù)分析

5 結(jié)束語

針對風(fēng)電機(jī)組傳動(dòng)鏈常處于強(qiáng)噪聲環(huán)境中，而現(xiàn)有VMD算法仍依靠人為經(jīng)驗(yàn)設(shè)定，導(dǎo)致過分解、欠分解、特征提取效果不佳等問題，本文引入樣本熵與峭度對VMD算法預(yù)設(shè)參數(shù)進(jìn)行自適應(yīng)優(yōu)化，有效提高輸入信號的信噪比，實(shí)現(xiàn)故障處的特征頻率幅值增強(qiáng)。實(shí)驗(yàn)結(jié)果表明，本文所提方法對強(qiáng)噪聲環(huán)境下提取故障特征具有較好的效果，在實(shí)際風(fēng)電機(jī)組傳動(dòng)鏈故障診斷領(lǐng)域有較高的實(shí)用價(jià)值。