王雙海，米大斌，蘆浩，姜文，龔思遠，梁濤

(1.河北建投能源投資股份有限公司，河北石家莊 050011;2.河北建設(shè)投資集團有限責任公司，河北石家莊 050011;3.河北工業(yè)大學(xué)人工智能與數(shù)據(jù)科學(xué)學(xué)院，天津 300131)

0 前言

汽輪機是火力發(fā)電的動力機械，而軸承又是汽輪機的重要組成部件。軸承一旦發(fā)生故障，必將影響整個機械設(shè)備的運行，甚至造成經(jīng)濟損失以及人員傷亡。所以，在實際生產(chǎn)中及時診斷出軸承的運行狀態(tài)至關(guān)重要。當軸承發(fā)生故障時，振動信號特性也會隨之改變，若能提取出各狀態(tài)的信號特征，則可進行軸承故障診斷。

文獻[1]中采用經(jīng)驗?zāi)B(tài)分解(EMD)把原始信號分解為多個平穩(wěn)的固有模態(tài)函數(shù)，并將其峭度和裕度因子作為神經(jīng)網(wǎng)絡(luò)的輸入?yún)?shù)來識別滾動軸承故障模式。文獻[2]中將EMD和奇異值分解相結(jié)合提取軸承故障特征向量，并用模糊C均值聚類來識別滾動軸承故障類型。文獻[3]中針對EMD的模態(tài)混疊，采用集成經(jīng)驗?zāi)B(tài)分解(EEMD)來提取發(fā)動機曲軸故障特征，最后識別出不同種類的故障。文獻[4]中針對EMD方法存在端點效應(yīng)、模態(tài)混疊和虛假分量的問題,采用局部均值分解(LMD)和樣本熵相結(jié)合來提取軸承故障特征向量，并用貝葉斯網(wǎng)絡(luò)來識別滾動軸承的故障類型。雖然以上兩種方法相對EMD效果要好，但依然無法從根本上解決模態(tài)混淆和端點效應(yīng)等問題。2014年, DRAGOMIRETSKIY提出了一種新型的自適應(yīng)信號處理方法——變分模態(tài)分解(Variational Mode Decomposition,VMD)方法[5]，作為一種新型的自適應(yīng)信號處理方法，VMD方法引入變分模型,將信號的分解轉(zhuǎn)換為約束模型最優(yōu)解的尋優(yōu)問題,可以避免端點效應(yīng)、抑制模態(tài)混淆,并且具有很高的分解效率[6]。

排列熵作為一種非線性理論的信號分析方法，能夠較好地度量時間序列的復(fù)雜性，可以將排列熵作為不同復(fù)雜度信號的特征參數(shù)。高維特征向量包含了大量故障特征信息，但是維數(shù)過高不僅會導(dǎo)致計算的復(fù)雜度過高的問題，還會對分類器性能造成不良影響[8]。因此，需要利用降維方法進行二次特征提取,去除冗余特征。傳統(tǒng)的降維方法如主成分分析法(PCA)、線性判別分析法(LDA),雖然對線性結(jié)構(gòu)的高維數(shù)據(jù)有很好的降維效果，但對實際獲取的非線性結(jié)構(gòu)數(shù)據(jù)，降維效果卻是差強人意。近年來新提出的降維方法t分布隨機鄰近嵌入(t-SNE)是一種基于流行學(xué)習(xí)的非線性數(shù)據(jù)降維方法，具有優(yōu)異的降維能力[8]。

本文作者以軸承振動信號為研究對象，提出了基于多島遺傳算法優(yōu)化VMD參數(shù)、排列熵和t-SNE的特征提取方法。仿真結(jié)果表明，該方法可有效提高軸承故障診斷的正確率。

1 VMD簡介

變分模態(tài)分解過程實質(zhì)上是構(gòu)造變分問題與求解變分問題的過程，從構(gòu)造和分解兩個方面介紹VMD方法。

1.1 變分問題的構(gòu)造

VMD算法將分解后的IMF(Intrinsic Modal Functions)分量重新定義為一個調(diào)幅-調(diào)頻信號，假設(shè)原始信號被分解為K個IMF分量，則第k個IMF分量的表達式為

uk(t)=Ak(t)cos[φk(t)]

(1)

式中:k∈{1,…,K};φk(t)是一個非遞減的相位函數(shù)；Ak(t)表示包絡(luò)函數(shù)。

首先，對每一個IMF分量，利用Hilbert變換求其解析信號以獲得單邊頻譜，其頻譜表達式為

(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ù)；?為卷積運算；x代表原始信號。

1.2 變分問題的求解

引入二次懲罰因子α和拉格朗日乘法系數(shù)λ(t)，將約束變分問題變?yōu)榉羌s束變分問題。其中二次懲罰因子可在高斯噪聲存在的情況下保證信號的重構(gòu)精度，拉格朗日系數(shù)使得約束條件保持嚴格性，增廣的拉格朗日函數(shù)如下：

L({uk},{ωk},λ)=〈λ(t),x(t)-

(5)

然后用交替方向乘子方法(ADMM)對上式進行極值求解得到IMF分量的頻域表達式如下所示：

(6)

(7)

(8)

通過以上分析,VMD的算法流程如下：

(2)根據(jù)式(6)和式(7)更新uk和ωk；

(3)根據(jù)式(8)更新λ；

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

在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。因此，在對信號進行VMD分解前，需要選取適當?shù)腒和α。

在VMD參數(shù)的優(yōu)化過程中，采用粒子群優(yōu)化算法對VMD的參數(shù)進行優(yōu)化[9]。但是粒子群算法容易陷入局部最優(yōu)，搜索路徑也比較復(fù)雜?！斑z傳算法”的概念是由密西根大學(xué)的HOLLAND于1975年首次提出的，模仿自然界中“優(yōu)勝劣汰，適者生存”的法則進行優(yōu)化過程[7]。在參數(shù)優(yōu)化過程中，采用遺傳算法對VMD的參數(shù)進行優(yōu)化。它通過染色體選擇、交叉、變異等迭代操作得到參數(shù)的優(yōu)化解。該算法具有隱式并行搜索的特點，在求解過程中具有很強的全局搜索能力，但是遺傳算法容易過早收斂和陷入局部最優(yōu)。針對遺傳算法經(jīng)常出現(xiàn)的過早收斂、產(chǎn)生局部最優(yōu)解的問題，有學(xué)者提出了多島遺傳算法(Multi-Island Genetic Algorithm, MIGA)。多島遺傳算法通過反復(fù)恰當?shù)厥褂眠z傳算法的算子和選擇原則，從親代到子代，從子代到孫代，從孫代到重孫代，不停地繁衍，使得種群對環(huán)境的適應(yīng)性不斷升高[10]。本文作者采用多島遺傳算法進行參數(shù)優(yōu)化，獲取最優(yōu)的參數(shù)組合。多島遺傳算法尋優(yōu)的具體步驟如下：

步驟1，初始化種群P0。

步驟2,將初始種群P0劃分為多個“島”，即子種群。

步驟3,計算各個“島”上每個個體的適應(yīng)度值。

步驟4,在每個島上進行標準遺傳算法中的選擇、交叉和變異等操作。

步驟5,若滿足遷移條件，則從一個島遷移到另一個島。

步驟6,若當前迭代次數(shù)n小于最大迭代次數(shù)MaxGen，則繼續(xù)進行步驟3,否則進行步驟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算法具有較快的收斂速度、較強的全局優(yōu)化能力。

圖1 優(yōu)化過程的迭代曲線

用多島遺傳算法進行參數(shù)優(yōu)化時，需選定一個適應(yīng)度函數(shù)。文獻[11]中提出了包絡(luò)熵Ep的概念，故障信號經(jīng)VMD分解后，IMF分量的包絡(luò)熵值Ep可反映出該分量的稀疏特性。如果分解后IMF分量中包含的噪聲較多，從而掩蓋故障沖擊特征，則該IMF分量的稀疏性較弱，包絡(luò)熵較大；相反，若IMF分量中包含較多的規(guī)律性故障沖擊，則該IMF分量具有較強的稀疏性，包絡(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)進行歸一化得到的；Ep是根據(jù)信息熵的計算規(guī)則得到的。

為了驗證局部包絡(luò)熵作為適應(yīng)度函數(shù)的有效性，構(gòu)造了一個滾動軸承的模擬信號。其表達式如下：

(12)

式中:c(t)表示加噪聲的軸承模擬信號;y0為位移常數(shù);fn為軸承的固有頻率;ξ為阻尼系數(shù);n(t)為趨于真實噪聲的高斯白噪聲;t為采樣時間;采樣頻率fs=20 kHz，采樣點數(shù)量N=4 096。這里，設(shè)置y0=3，fn=3 000 Hz，ξ=0.09。

圖2展示了在無噪聲、噪聲強度為1、噪聲強度為2和噪聲強度為3的情況下模擬信號的波形。通過4幅圖的比較，可以發(fā)現(xiàn)信號中的周期性脈沖完全被強背景噪聲所淹沒；噪聲強度與包絡(luò)熵呈正相關(guān)：隨著噪聲的增強(減弱)，周期脈沖變得模糊(顯著)，包絡(luò)熵變大(變小)。這也證明了信號稀疏度越大、包絡(luò)熵越小的規(guī)律。因此，包絡(luò)熵最小的模態(tài)分量可以確定為最佳分量。

圖2 不同噪聲強度下模擬信號的時域波形

3 相關(guān)系數(shù)準則

相關(guān)系數(shù)是反映變量之間相關(guān)關(guān)系密切程度的統(tǒng)計指標[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ù)來選取最能代表軸承振動信號的部分分量。

4 排列熵

4.1 排列熵原理

排列熵是描述時間序列復(fù)雜性和不規(guī)律程度的一種參數(shù)，具有計算簡單、抗噪能力強等優(yōu)點[14]。排列熵的算法簡介如下。

對于長度為n的時間序列{[x(i),i=1,2,…,n]}，參考相空間延遲坐標法對其進行重構(gòu)，每個采樣點取其連續(xù)的m個樣本點，得到點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ù)雜度越大。

4.2 參數(shù)選擇

排列熵的計算值大小主要受3個參數(shù)的影響，分別是信號序列長度n、嵌入維度m、延遲時間τ。作者依據(jù)文獻[12]分別設(shè)置嵌入維度m=6，延遲時間τ=1。

4.3 幾組典型信號的排列熵值

為了證明排列熵具有表征信號特征的能力，下面計算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ù)雜度越大，排列熵值越大，排列熵可以作為信號分類的特征。

5 t-SNE流行學(xué)習(xí)算法

t分布隨機鄰近嵌入(t-SNE)算法是一種深度學(xué)習(xí)的非線性流行學(xué)習(xí)算法，對高維非線性數(shù)據(jù)集具有優(yōu)異的降維效果[13]。t-SNE算法具體步驟[17]如下:

(1)設(shè)高維樣本X={x1,x2,…,xn}，任意兩個樣本數(shù)據(jù)點xi和xj相似的條件概率為

(18)

式中：δi為以xi為中心點的高斯分布方差，可根據(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}。

6 基于MIGA-VMD和t-SNE的滾動軸承故障診斷流程

基于上述理論，基于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進行訓(xùn)練，得到SVM分類模型。

(7)將測試集樣本的特征向量輸入到訓(xùn)練好的SVM模型，實現(xiàn)故障診斷。

7 仿真實驗

7.1 實驗數(shù)據(jù)描述

采用凱斯西儲大學(xué)電氣工程實驗室的數(shù)據(jù)進行仿真。選取軸承的負載為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組樣本作為測試樣本。

7.2 故障特征提取

首先，利用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)圈故障信號進行VMD處理得到的6個模態(tài)分量對應(yīng)的信號波形和頻譜。

圖4 尋優(yōu)過程的迭代曲線

圖5 內(nèi)圈故障信號波形

對4種狀態(tài)信號分別進行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分量并計算其排列熵做為特征向量。

7.3 故障特征簡約

利用t-SNE方法對六維的特征向量R240×6進行降維得到低維特征向量R240×3，為了驗證t-SNE的降維效果，同時采用PCA、LDA、t-SNE 3種方法進行降維處理，降維結(jié)果如圖8所示。

圖8 3種方法降維結(jié)果

從圖8可以看出：在3種降維方法中，經(jīng)t-SNE降維后4種類型的數(shù)據(jù)特征聚集性最好，特征區(qū)分最為明顯。因此，本文作者采用t-SNE方法進行降維，得到三維的特征向量。

7.4 故障分類

將160組訓(xùn)練樣本的特征向量輸入到SVM分類器進行訓(xùn)練，得到用于故障診斷的SVM模型。通過訓(xùn)練好的SVM模型對4種狀態(tài)下的80組測試樣本分類，測試樣本的分類結(jié)果如圖9所示。在SVM的訓(xùn)練和測試中，采用數(shù)字標簽代表軸承的4種運行狀態(tài)：正常狀態(tài)(標簽1)，內(nèi)圈故障狀態(tài)(標簽2)，滾珠故障狀態(tài)(標簽3)，外圈故障狀態(tài)(標簽4)。

圖9 測試樣本的分類結(jié)果

由圖9可以看出：采用基于MIGA-VMD和排列熵、t-SNE的特征提取方法，軸承狀態(tài)的診斷精度可達到100%。

7.5 方法對比

傳統(tǒng)的VMD參數(shù)設(shè)置方法參考文獻[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種方法進行信號分解，同樣采用排列熵構(gòu)建特征向量，并用t-SNE降維處理得到低維特征向量來訓(xùn)練SVM模型，用訓(xùn)練好的SVM模型對4種狀態(tài)的測試樣本分類。5種方法得到的低維特征分布與SVM分類準確率如圖10、表2所示。

由圖10、表2可以看出:基于MIGA-VMD和排列熵、t-SNE的特征提取方法，診斷正確率優(yōu)于其他4種方法，可以從干擾信號中提取出微小故障特征，實現(xiàn)對不同故障正確診斷識別。

圖10 低維特征分布

表2 5種方法的SVM分類準確率

8 結(jié)論

提出基于MIGA-VMD和排列熵、t-SNE的特征提取方法,借助參數(shù)優(yōu)化的VMD方法分解振動信號,計算分解后與原始信號相關(guān)性較高的部分IMF分量的排列熵作為故障特征，利用t-SNE方法進行降維得到低維特征向量并將其輸入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的特征提取方法,軸承的故障診斷精度有明顯的提高。