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

(1.山西大學自動化與軟件學院，山西 太原 030013)(2.煤炭工業(yè)太原設計研究院，山西 太原 030001)

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

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

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

1 基本原理

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

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

(1)

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

(2)

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

(3)

(4)

(5)

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

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

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

(6)

1.2 樣本熵

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

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

(7)

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

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

(8)

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

2 故障特征提取流程

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

(9)

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

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

(10)

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

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

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

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

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

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

3 仿真信號驗證

為驗證本文所提方法的有效性，對齒輪箱行星輪故障仿真信號進行特征提取，不失一般性，本文不考慮諧波因素，只考慮載波頻率和調制頻率的基頻與振動傳遞路徑對振動信號的影響[9]，故齒輪箱行星輪故障可用如式(11)的調幅調頻信號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，為行星架的旋轉頻率；fs=75 Hz，為行星輪的故障特征頻率；fm=1 000 Hz，為齒輪嚙合頻率；A,B為信號的調制系數(shù)，均為1；初始相位φ=0；n(t)為高斯白噪聲。

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

圖2 仿真信號分析

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

表1 各模態(tài)分量指標計算

圖3 重構信號分析

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

4 實測數(shù)據(jù)驗證

4.1 試驗臺驗證

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

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

圖4 風電機組傳動系統(tǒng)故障試驗臺

表2 試驗臺參數(shù)

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

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

4.2 工程實例驗證

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

表3 風電場實測數(shù)據(jù)

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

表4 齒輪箱參數(shù)

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

5 結束語

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