王雙海,米大斌,蘆浩,姜文,龔思遠(yuǎn),梁濤
(1.河北建投能源投資股份有限公司,河北石家莊 050011;2.河北建設(shè)投資集團(tuán)有限責(zé)任公司,河北石家莊 050011;3.河北工業(yè)大學(xué)人工智能與數(shù)據(jù)科學(xué)學(xué)院,天津 300131)
汽輪機(jī)是火力發(fā)電的動力機(jī)械,而軸承又是汽輪機(jī)的重要組成部件。軸承一旦發(fā)生故障,必將影響整個機(jī)械設(shè)備的運(yùn)行,甚至造成經(jīng)濟(jì)損失以及人員傷亡。所以,在實際生產(chǎn)中及時診斷出軸承的運(yùn)行狀態(tài)至關(guān)重要。當(dāng)軸承發(fā)生故障時,振動信號特性也會隨之改變,若能提取出各狀態(tài)的信號特征,則可進(jìn)行軸承故障診斷。
文獻(xiàn)[1]中采用經(jīng)驗?zāi)B(tài)分解(EMD)把原始信號分解為多個平穩(wěn)的固有模態(tài)函數(shù),并將其峭度和裕度因子作為神經(jīng)網(wǎng)絡(luò)的輸入?yún)?shù)來識別滾動軸承故障模式。文獻(xiàn)[2]中將EMD和奇異值分解相結(jié)合提取軸承故障特征向量,并用模糊C均值聚類來識別滾動軸承故障類型。文獻(xiàn)[3]中針對EMD的模態(tài)混疊,采用集成經(jīng)驗?zāi)B(tài)分解(EEMD)來提取發(fā)動機(jī)曲軸故障特征,最后識別出不同種類的故障。文獻(xiàn)[4]中針對EMD方法存在端點(diǎn)效應(yīng)、模態(tài)混疊和虛假分量的問題,采用局部均值分解(LMD)和樣本熵相結(jié)合來提取軸承故障特征向量,并用貝葉斯網(wǎng)絡(luò)來識別滾動軸承的故障類型。雖然以上兩種方法相對EMD效果要好,但依然無法從根本上解決模態(tài)混淆和端點(diǎn)效應(yīng)等問題。2014年, DRAGOMIRETSKIY提出了一種新型的自適應(yīng)信號處理方法——變分模態(tài)分解(Variational Mode Decomposition,VMD)方法[5],作為一種新型的自適應(yīng)信號處理方法,VMD方法引入變分模型,將信號的分解轉(zhuǎn)換為約束模型最優(yōu)解的尋優(yōu)問題,可以避免端點(diǎn)效應(yīng)、抑制模態(tài)混淆,并且具有很高的分解效率[6]。
排列熵作為一種非線性理論的信號分析方法,能夠較好地度量時間序列的復(fù)雜性,可以將排列熵作為不同復(fù)雜度信號的特征參數(shù)。高維特征向量包含了大量故障特征信息,但是維數(shù)過高不僅會導(dǎo)致計算的復(fù)雜度過高的問題,還會對分類器性能造成不良影響[8]。因此,需要利用降維方法進(jìn)行二次特征提取,去除冗余特征。傳統(tǒng)的降維方法如主成分分析法(PCA)、線性判別分析法(LDA),雖然對線性結(jié)構(gòu)的高維數(shù)據(jù)有很好的降維效果,但對實際獲取的非線性結(jié)構(gòu)數(shù)據(jù),降維效果卻是差強(qiáng)人意。近年來新提出的降維方法t分布隨機(jī)鄰近嵌入(t-SNE)是一種基于流行學(xué)習(xí)的非線性數(shù)據(jù)降維方法,具有優(yōu)異的降維能力[8]。
本文作者以軸承振動信號為研究對象,提出了基于多島遺傳算法優(yōu)化VMD參數(shù)、排列熵和t-SNE的特征提取方法。仿真結(jié)果表明,該方法可有效提高軸承故障診斷的正確率。
變分模態(tài)分解過程實質(zhì)上是構(gòu)造變分問題與求解變分問題的過程,從構(gòu)造和分解兩個方面介紹VMD方法。
VMD算法將分解后的IMF(Intrinsic Modal Functions)分量重新定義為一個調(diào)幅-調(diào)頻信號,假設(shè)原始信號被分解為K個IMF分量,則第k個IMF分量的表達(dá)式為
uk(t)=Ak(t)cos[φk(t)]
(1)
式中:k∈{1,…,K};φk(t)是一個非遞減的相位函數(shù);Ak(t)表示包絡(luò)函數(shù)。
首先,對每一個IMF分量,利用Hilbert變換求其解析信號以獲得單邊頻譜,其頻譜表達(dá)式為
(2)
其次,根據(jù)混合預(yù)估中心頻率,調(diào)制頻譜至對應(yīng)基頻帶,記為
(3)
最后,計算解調(diào)信號的時間梯度L2范數(shù)的平方,估計出模態(tài)分量的帶寬,并引入約束條件,構(gòu)造出約束的變分模型為
(4)
式中:{uk}代表VMD分解得到的k個IMF分量;{ωk}代表IMF分量對應(yīng)的中心頻率;δ(t)為狄利克雷函數(shù);?為卷積運(yùn)算;x代表原始信號。
引入二次懲罰因子α和拉格朗日乘法系數(shù)λ(t),將約束變分問題變?yōu)榉羌s束變分問題。其中二次懲罰因子可在高斯噪聲存在的情況下保證信號的重構(gòu)精度,拉格朗日系數(shù)使得約束條件保持嚴(yán)格性,增廣的拉格朗日函數(shù)如下:
L({uk},{ωk},λ)=〈λ(t),x(t)-
(5)
然后用交替方向乘子方法(ADMM)對上式進(jìn)行極值求解得到IMF分量的頻域表達(dá)式如下所示:
(6)
(7)
(8)
通過以上分析,VMD的算法流程如下:
(2)根據(jù)式(6)和式(7)更新uk和ωk;
(3)根據(jù)式(8)更新λ;
在VMD分解過程中,分解的IMF分量個數(shù)K和懲罰參數(shù)α對分解結(jié)果有較大影響。若K值過大,分解結(jié)果中會產(chǎn)生虛假分量,誤導(dǎo)對實驗結(jié)果的判斷;若K值較小,則會產(chǎn)生模態(tài)混疊現(xiàn)象。若懲罰參數(shù)α越小,分解后各個IMF分量的帶寬就越小;相反,若懲罰參數(shù)α越大,分解后各個IMF分量的帶寬就越大。除K、α以外,其他參數(shù)對分解效果影響較小,設(shè)置為經(jīng)驗值,即tau=0,init=1,DC=0,ε=1×10-7。因此,在對信號進(jìn)行VMD分解前,需要選取適當(dāng)?shù)腒和α。
在VMD參數(shù)的優(yōu)化過程中,采用粒子群優(yōu)化算法對VMD的參數(shù)進(jìn)行優(yōu)化[9]。但是粒子群算法容易陷入局部最優(yōu),搜索路徑也比較復(fù)雜?!斑z傳算法”的概念是由密西根大學(xué)的HOLLAND于1975年首次提出的,模仿自然界中“優(yōu)勝劣汰,適者生存”的法則進(jìn)行優(yōu)化過程[7]。在參數(shù)優(yōu)化過程中,采用遺傳算法對VMD的參數(shù)進(jìn)行優(yōu)化。它通過染色體選擇、交叉、變異等迭代操作得到參數(shù)的優(yōu)化解。該算法具有隱式并行搜索的特點(diǎn),在求解過程中具有很強(qiáng)的全局搜索能力,但是遺傳算法容易過早收斂和陷入局部最優(yōu)。針對遺傳算法經(jīng)常出現(xiàn)的過早收斂、產(chǎn)生局部最優(yōu)解的問題,有學(xué)者提出了多島遺傳算法(Multi-Island Genetic Algorithm, MIGA)。多島遺傳算法通過反復(fù)恰當(dāng)?shù)厥褂眠z傳算法的算子和選擇原則,從親代到子代,從子代到孫代,從孫代到重孫代,不停地繁衍,使得種群對環(huán)境的適應(yīng)性不斷升高[10]。本文作者采用多島遺傳算法進(jìn)行參數(shù)優(yōu)化,獲取最優(yōu)的參數(shù)組合。多島遺傳算法尋優(yōu)的具體步驟如下:
步驟1,初始化種群P0。
步驟2,將初始種群P0劃分為多個“島”,即子種群。
步驟3,計算各個“島”上每個個體的適應(yīng)度值。
步驟4,在每個島上進(jìn)行標(biāo)準(zhǔn)遺傳算法中的選擇、交叉和變異等操作。
步驟5,若滿足遷移條件,則從一個島遷移到另一個島。
步驟6,若當(dāng)前迭代次數(shù)n小于最大迭代次數(shù)MaxGen,則繼續(xù)進(jìn)行步驟3,否則進(jìn)行步驟7。
步驟7,輸出種群中適應(yīng)度值最優(yōu)的個體作為最優(yōu)解。
多島遺傳算法的參數(shù)設(shè)計如下:初始種群劃分為10個島,每個島的種群數(shù)設(shè)為10,優(yōu)化10代,則計算個體數(shù)為1 000個。變異概率設(shè)為0.01,交叉概率設(shè)為1,遷移概率設(shè)為0.01,精英個體數(shù)量設(shè)為1。
為了驗證MIGA的收斂性和優(yōu)化性能,測試了不同的算法,包括PSO算法、GA算法、MIGA算法。通過以下函數(shù)作為適應(yīng)度函數(shù)比較它們的優(yōu)化性能
(9)
對于PSO算法,最大迭代次數(shù)設(shè)為50,種群規(guī)模為100。對于遺傳算法,最大迭代次數(shù)設(shè)為50,種群規(guī)模為100,交叉概率為1,變異概率為0.01。對于多島遺傳算法,初始種群劃分為10個島,每個島的種群數(shù)設(shè)為10,最大迭代次數(shù)為50,交叉概率為1,變異概率為0.01,遷移概率為0.01,精英個體數(shù)量為1。優(yōu)化過程的迭代曲線如圖1所示??梢钥闯觯簩τ诘仁?9)中的適應(yīng)度函數(shù),PSO算法的收斂速度優(yōu)于GA算法,但是PSO算法和GA算法都陷入了局部最優(yōu)解;文中所用的MIGA算法不僅收斂速度優(yōu)于PSO算法和GA算法,并且獲得的最優(yōu)值更接近理論最優(yōu)值。結(jié)果表明:MIGA算法具有較快的收斂速度、較強(qiáng)的全局優(yōu)化能力。
圖1 優(yōu)化過程的迭代曲線
用多島遺傳算法進(jìn)行參數(shù)優(yōu)化時,需選定一個適應(yīng)度函數(shù)。文獻(xiàn)[11]中提出了包絡(luò)熵Ep的概念,故障信號經(jīng)VMD分解后,IMF分量的包絡(luò)熵值Ep可反映出該分量的稀疏特性。如果分解后IMF分量中包含的噪聲較多,從而掩蓋故障沖擊特征,則該IMF分量的稀疏性較弱,包絡(luò)熵較大;相反,若IMF分量中包含較多的規(guī)律性故障沖擊,則該IMF分量具有較強(qiáng)的稀疏性,包絡(luò)熵較小。在K與α的影響下,K個分量就會有K個包絡(luò)熵值,選擇K個包絡(luò)熵值中最小的一個作為局部最小包絡(luò)熵minEp,該最小熵值對應(yīng)的分量有著豐富的故障特征信息[6]。將局部最小熵作為多島遺傳算法中的適應(yīng)度函數(shù),整個搜索過程就是要找到全局最小包絡(luò)熵以及對應(yīng)的最佳分量所在的最優(yōu)參數(shù)組合[K,α]。所以,構(gòu)建如下的適應(yīng)度函數(shù):
minF=minEp
(10)
式中包絡(luò)熵Ep的計算公式為
(11)
式中:j=1,2,3,…,N;a(j)是x(j)經(jīng)Hilbert解調(diào)得到的包絡(luò)信號;ej是對a(j)進(jìn)行歸一化得到的;Ep是根據(jù)信息熵的計算規(guī)則得到的。
為了驗證局部包絡(luò)熵作為適應(yīng)度函數(shù)的有效性,構(gòu)造了一個滾動軸承的模擬信號。其表達(dá)式如下:
(12)
式中:c(t)表示加噪聲的軸承模擬信號;y0為位移常數(shù);fn為軸承的固有頻率;ξ為阻尼系數(shù);n(t)為趨于真實噪聲的高斯白噪聲;t為采樣時間;采樣頻率fs=20 kHz,采樣點(diǎn)數(shù)量N=4 096。這里,設(shè)置y0=3,fn=3 000 Hz,ξ=0.09。
圖2展示了在無噪聲、噪聲強(qiáng)度為1、噪聲強(qiáng)度為2和噪聲強(qiáng)度為3的情況下模擬信號的波形。通過4幅圖的比較,可以發(fā)現(xiàn)信號中的周期性脈沖完全被強(qiáng)背景噪聲所淹沒;噪聲強(qiáng)度與包絡(luò)熵呈正相關(guān):隨著噪聲的增強(qiáng)(減弱),周期脈沖變得模糊(顯著),包絡(luò)熵變大(變小)。這也證明了信號稀疏度越大、包絡(luò)熵越小的規(guī)律。因此,包絡(luò)熵最小的模態(tài)分量可以確定為最佳分量。
圖2 不同噪聲強(qiáng)度下模擬信號的時域波形
相關(guān)系數(shù)是反映變量之間相關(guān)關(guān)系密切程度的統(tǒng)計指標(biāo)[16]。文中采用皮爾遜相關(guān)系數(shù)分析各IMF分量與原始信號的相關(guān)性程度。設(shè)樣本X和樣本Y,則兩者的相關(guān)系數(shù)為
(13)
式中:r為相關(guān)系數(shù);Cov(X,Y)為樣本X和樣本Y的協(xié)方差;D(X)為樣本X的方差;D(Y)為樣本Y的方差。
相關(guān)系數(shù)r的值越大,代表兩個樣本的相關(guān)性越高。通過計算IMF分量與原始信號的相關(guān)系數(shù)來選取最能代表軸承振動信號的部分分量。
排列熵是描述時間序列復(fù)雜性和不規(guī)律程度的一種參數(shù),具有計算簡單、抗噪能力強(qiáng)等優(yōu)點(diǎn)[14]。排列熵的算法簡介如下。
對于長度為n的時間序列{[x(i),i=1,2,…,n]},參考相空間延遲坐標(biāo)法對其進(jìn)行重構(gòu),每個采樣點(diǎn)取其連續(xù)的m個樣本點(diǎn),得到點(diǎn)x(i)的m維重構(gòu)向量[12]:
X(i)={x(i),x(i+τ),…,x[i+(m-1)τ]}
(14)
其中:m是嵌入維數(shù);τ是延遲時間;i=1,2,…,n-(m-1)τ,即X(i)向量的個數(shù)為n-(m-1)τ。將X(i)中的元素按從大到小的順序排列,即
x[i+(j1-1)τ]≤x[i+(j2-1)τ]≤…≤
x[i+(jm-1)τ]
(15)
式中:j1、j2,…,jm為X(i)中各元素的位置。
如果X(i)中出現(xiàn)大小相等的兩個元素,如x[i+(j1-1)τ]=x[i+(j2-1)τ],則按照j的大小來排序,即x[i+(j1-1)τ]≤x[i+(j2-1)τ]。
所以任一向量X(i)一定有一種符號序列:
S(l)=(j1,j2,…,jm)
(16)
其中:l=1,2,…,k且k≤m!。根據(jù)排列組合定理,m個元素最多有m!種排列方式,S(l)只是其中的一種。若S(l)這種排列方式出現(xiàn)的概率為Pl,則序列{x(i),i=1,2,…,n}的排列熵可定義為
(17)
若H(m,τ)越小說明時間序列越規(guī)則,若H(m,τ)越大說明時間序列無序程度越高,信號復(fù)雜度越大。
排列熵的計算值大小主要受3個參數(shù)的影響,分別是信號序列長度n、嵌入維度m、延遲時間τ。作者依據(jù)文獻(xiàn)[12]分別設(shè)置嵌入維度m=6,延遲時間τ=1。
為了證明排列熵具有表征信號特征的能力,下面計算6組具有代表性信號的排列熵值。
x1為長度為3 000、服從均勻分布白噪聲信號。
x2(t)=0.5t2+0.9t-0.3
x3(t)=sin(2π·900t)
x4(t)=sin(2π·30t)
x5(t)=sin(2π·30t)·sin(2π·90t)
x6(t)=0.5πt+5
經(jīng)計算,x1白噪聲的排列熵值最大為6.456 2;其次x5高頻正弦信號排列熵值為3.021 4;x3、x4非高頻正弦信號、調(diào)幅信號為1.099 2,1.671 8;x2、x6時域的曲線和直線信號排列熵值最小都為0。由此分析,隨著信號復(fù)雜度越大,排列熵值越大,排列熵可以作為信號分類的特征。
t分布隨機(jī)鄰近嵌入(t-SNE)算法是一種深度學(xué)習(xí)的非線性流行學(xué)習(xí)算法,對高維非線性數(shù)據(jù)集具有優(yōu)異的降維效果[13]。t-SNE算法具體步驟[17]如下:
(1)設(shè)高維樣本X={x1,x2,…,xn},任意兩個樣本數(shù)據(jù)點(diǎn)xi和xj相似的條件概率為
(18)
式中:δi為以xi為中心點(diǎn)的高斯分布方差,可根據(jù)用戶指定的困惑度Perp經(jīng)二分搜索確定,困惑度定義為
(19)
(2)計算高維樣本數(shù)據(jù)的聯(lián)合概率密度:
(20)
(3)初始化低維數(shù)據(jù)解y0={y1,y2,…,yn}。
(4)計算低維數(shù)據(jù)相似度:
(21)
(5)使用KL散度度量高維數(shù)據(jù)分布P和低維數(shù)據(jù)分布Q的相似度:
(22)
(6)計算梯度
(23)
(7)獲取低維數(shù)據(jù):
(24)
式中:η為學(xué)習(xí)率;α(t)為優(yōu)化參數(shù);t代表第t次迭代。
(8)重復(fù)步驟(4)—(7),直到滿足迭代次數(shù)T。經(jīng)T次迭代得到低維數(shù)據(jù)y(T)={y1,y2,…,yn}。
基于上述理論,基于MIGA-VMD分解、排列熵和t-SNE的軸承故障診斷模型如圖3所示。具體步驟如下:
圖3 故障診斷流程
(1)獲取各狀態(tài)的振動信號。
(2)用MIGA算法搜索VMD方法在各狀態(tài)下的最佳參數(shù)組合[K0,α0],采用經(jīng)參數(shù)優(yōu)化的VMD方法分解各狀態(tài)下的振動信號得到K0個IMF分量。
(3)計算各IMF分量與原始信號的相關(guān)性系數(shù)。
(4)計算與原始信號相關(guān)性較高的6個IMF分量的排列熵,組成特征向量。
(5)利用t-SNE方法降維得到三維特征向量。
(6)將訓(xùn)練集的特征向量輸入到SVM進(jìn)行訓(xùn)練,得到SVM分類模型。
(7)將測試集樣本的特征向量輸入到訓(xùn)練好的SVM模型,實現(xiàn)故障診斷。
采用凱斯西儲大學(xué)電氣工程實驗室的數(shù)據(jù)進(jìn)行仿真。選取軸承的負(fù)載為0,轉(zhuǎn)速為1 797 r/min,損傷直徑為0.177 8 mm,采樣頻率為12 kHz的正常狀態(tài)、內(nèi)圈故障、滾珠故障、外圈故障振動信號。每種狀態(tài)各取前10 s的數(shù)據(jù),樣本長度設(shè)為2 000,即每種狀態(tài)各有60組樣本,共240組樣本,其中160組樣本作為訓(xùn)練樣本,80組樣本作為測試樣本。
首先,利用MIGA算法搜索VMD算法的最優(yōu)參數(shù)組合。以內(nèi)圈故障信號為例,VMD參數(shù)尋優(yōu)過程的迭代曲線如圖4所示。經(jīng)過8次迭代搜索到了全局的最優(yōu)解,此時的包絡(luò)熵值最小,為4.238。該全局最優(yōu)解對應(yīng)的參數(shù)α=2 999、K=6,將其引入到VMD算法的參數(shù)設(shè)置中。圖5、圖6為對內(nèi)圈故障信號進(jìn)行VMD處理得到的6個模態(tài)分量對應(yīng)的信號波形和頻譜。
圖4 尋優(yōu)過程的迭代曲線
圖5 內(nèi)圈故障信號波形
對4種狀態(tài)信號分別進(jìn)行VMD參數(shù)優(yōu)化,對應(yīng)的最佳參數(shù)組合如表1所示。
根據(jù)表1中的參數(shù)組合[K0,α0]設(shè)置VMD方法的K、α,利用參數(shù)優(yōu)化的VMD方法分解樣本得到多個IMF分量。
圖7為4種狀態(tài)下各分量與原始信號的相關(guān)程度。
圖7 4種狀態(tài)下各IMF分量與原始信號的相關(guān)程度
選取與原始信號相關(guān)性較高的6個IMF分量并計算其排列熵做為特征向量。
利用t-SNE方法對六維的特征向量R240×6進(jìn)行降維得到低維特征向量R240×3,為了驗證t-SNE的降維效果,同時采用PCA、LDA、t-SNE 3種方法進(jìn)行降維處理,降維結(jié)果如圖8所示。
圖8 3種方法降維結(jié)果
從圖8可以看出:在3種降維方法中,經(jīng)t-SNE降維后4種類型的數(shù)據(jù)特征聚集性最好,特征區(qū)分最為明顯。因此,本文作者采用t-SNE方法進(jìn)行降維,得到三維的特征向量。
將160組訓(xùn)練樣本的特征向量輸入到SVM分類器進(jìn)行訓(xùn)練,得到用于故障診斷的SVM模型。通過訓(xùn)練好的SVM模型對4種狀態(tài)下的80組測試樣本分類,測試樣本的分類結(jié)果如圖9所示。在SVM的訓(xùn)練和測試中,采用數(shù)字標(biāo)簽代表軸承的4種運(yùn)行狀態(tài):正常狀態(tài)(標(biāo)簽1),內(nèi)圈故障狀態(tài)(標(biāo)簽2),滾珠故障狀態(tài)(標(biāo)簽3),外圈故障狀態(tài)(標(biāo)簽4)。
圖9 測試樣本的分類結(jié)果
由圖9可以看出:采用基于MIGA-VMD和排列熵、t-SNE的特征提取方法,軸承狀態(tài)的診斷精度可達(dá)到100%。
傳統(tǒng)的VMD參數(shù)設(shè)置方法參考文獻(xiàn)[15]的中心頻率選擇法,通過觀察不同模態(tài)數(shù)K下對應(yīng)的中心頻率變化,來確定K的取值。若中心頻率變化較大,說明出現(xiàn)了欠分解,則增大K值;若中心頻率變化較小,說明出現(xiàn)了過分解,則減小K值。
為證明所提方法的有效性,選用同樣的樣本數(shù)據(jù),分別采用EMD、EEMD、LMD、傳統(tǒng)VMD(K=6,α=2 000)、MIGA-VMD 5種方法進(jìn)行信號分解,同樣采用排列熵構(gòu)建特征向量,并用t-SNE降維處理得到低維特征向量來訓(xùn)練SVM模型,用訓(xùn)練好的SVM模型對4種狀態(tài)的測試樣本分類。5種方法得到的低維特征分布與SVM分類準(zhǔn)確率如圖10、表2所示。
由圖10、表2可以看出:基于MIGA-VMD和排列熵、t-SNE的特征提取方法,診斷正確率優(yōu)于其他4種方法,可以從干擾信號中提取出微小故障特征,實現(xiàn)對不同故障正確診斷識別。
圖10 低維特征分布
表2 5種方法的SVM分類準(zhǔn)確率
提出基于MIGA-VMD和排列熵、t-SNE的特征提取方法,借助參數(shù)優(yōu)化的VMD方法分解振動信號,計算分解后與原始信號相關(guān)性較高的部分IMF分量的排列熵作為故障特征,利用t-SNE方法進(jìn)行降維得到低維特征向量并將其輸入SVM實現(xiàn)軸承的故障診斷。經(jīng)MIGA算法優(yōu)化參數(shù)的VMD方法與EMD、EEMD、LMD、傳統(tǒng)VMD等方法相比,提取的故障特征更為豐富,有利于對軸承狀態(tài)的診斷識別。t-SNE非線性降維方法與傳統(tǒng)線性降維方法PCA、LDA相比,在軸承診斷中具有更佳的降維效果。仿真實驗表明,采用此MIGA-VMD和排列熵、t-SNE的特征提取方法,軸承的故障診斷精度有明顯的提高。