基于改進(jìn)Fisher準(zhǔn)則、VMD、距離相關(guān)系數(shù)和 核極限學(xué)習(xí)機(jī)的軸承故障診斷

宋坤駿， 丁建明， 林建輝

(西南交通大學(xué)　牽引動(dòng)力國(guó)家重點(diǎn)實(shí)驗(yàn)室， 成都 610031)

滾動(dòng)軸承作為機(jī)械設(shè)備中常見且高故障率的部件，其狀態(tài)的及時(shí)監(jiān)測(cè)和故障的準(zhǔn)確診斷具有十分重要的意義。當(dāng)前振動(dòng)信號(hào)故障特征提取手法有很多，包括快速傅里葉變換，小波變換，Wigner-Ville分布，由Huang等人[1]提出的經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)和Dragomiretskiy等人[2]提出的變分模態(tài)分解(VMD)等。就EMD和VMD兩者的特點(diǎn)而言，EMD算法因?yàn)槿狈?jiān)實(shí)的數(shù)學(xué)基礎(chǔ)面臨著模態(tài)混疊，端點(diǎn)效應(yīng)等問題。VMD算法源自嚴(yán)格的數(shù)學(xué)推導(dǎo)，基本消除了模態(tài)混疊現(xiàn)象，同時(shí)由于VMD可看成一組自適應(yīng)的維納濾波器，因而對(duì)噪聲具有較高的魯棒性。VMD在實(shí)際工程應(yīng)用中面臨的主要問題在于算法參數(shù)的選取。針對(duì)這一問題，文獻(xiàn)[3]提出了以去噪后的信噪比和均方差為優(yōu)化目標(biāo)采用遺傳算法優(yōu)化VMD參數(shù)。文獻(xiàn)[4]提出了以模態(tài)分量的包絡(luò)熵值最小化為優(yōu)化目標(biāo)采用遺傳算法優(yōu)化VMD參數(shù)。上述兩種優(yōu)化方法均從信號(hào)預(yù)處理的角度出發(fā)，沒有結(jié)合分類算法，而參數(shù)選擇對(duì)于分類算法的效果實(shí)際上也有很大影響，典型例子當(dāng)屬支持向量機(jī)(SVM)。文獻(xiàn)[5]采用差分進(jìn)化算法通過最小化SVM在驗(yàn)證集上的泛化誤差優(yōu)化SVM的參數(shù)，該方法的特點(diǎn)是原理較為簡(jiǎn)單，優(yōu)化結(jié)果準(zhǔn)確，缺點(diǎn)是計(jì)算量偏大。文獻(xiàn)[6]為減小計(jì)算量，提出了基于Fisher準(zhǔn)則和最大熵原理的SVM核參數(shù)選取方法。文中借鑒其思想并加以提升，提出了以改進(jìn)的核空間多類別Fisher準(zhǔn)則為優(yōu)化目標(biāo)，優(yōu)化小波核極限學(xué)習(xí)機(jī)的核參數(shù)ω和a。

故障診斷流程如下：首先通過模態(tài)中心頻率錯(cuò)開原則以及信噪比最大原則確定VMD預(yù)設(shè)模態(tài)數(shù)K和二次懲罰參數(shù)α。然后將信號(hào)通過VMD分解為若干模態(tài)分量，提取各分量的峭度，中心頻率和同原信號(hào)的距離相關(guān)系數(shù)作為輸入核極限學(xué)習(xí)機(jī)的特征向量，并使用模擬退火粒子群算法基于改進(jìn)的核空間Fisher準(zhǔn)則優(yōu)化小波核極限學(xué)習(xí)機(jī)的核參數(shù)ω和a。最后使用參數(shù)優(yōu)化的小波核極限學(xué)習(xí)機(jī)對(duì)軸承的4種故障進(jìn)行分類辨識(shí)。

1　變分模態(tài)分解(VMD)

VMD算法能夠非迭代地求出非平穩(wěn)信號(hào)的各個(gè)本征模態(tài)函數(shù)，不僅避免了EMD等迭代型算法的篩分停止準(zhǔn)則選取，模態(tài)混疊等問題，也有效提高了分解效率。同EMD不同的是，VMD需要預(yù)先設(shè)定希望獲得的本征模態(tài)函數(shù)(IMF)分量數(shù)K，它通過最小化K個(gè)模態(tài)函數(shù)的帶寬和構(gòu)造變分問題。求取模態(tài)函數(shù)帶寬和的手段是：先進(jìn)行Hilbert變換求出邊際譜，然后乘以指數(shù)因子將模態(tài)函數(shù)的頻譜移動(dòng)到待求的中心頻率，最后計(jì)算上述信號(hào)的梯度平方L2范數(shù)作為帶寬的估計(jì)。整個(gè)需要求解的約束變分模型如下：

(1)

式中uk(t)為第k個(gè)IMF分量(1≤k≤K),ωk為其待求的中心頻率。Dragomiretskiy和Zosso通過引入二次懲罰因子α和Lagrange乘子λ消除約束，將上述極小值問題轉(zhuǎn)化為求解下述增廣Lagrange函數(shù)的鞍點(diǎn)

(2)

這里二次懲罰項(xiàng)可以提升收斂性能和對(duì)噪聲的魯棒性，拉氏乘子使得重構(gòu)約束條件得到嚴(yán)格滿足?？赏ㄟ^交替方向乘子法(ADMM)求解上述問題，每次求解一個(gè)變量，其他變量視為已知，可以在頻域求得第n+1次迭代的結(jié)果為

(3)

(4)

(2)n=n+1迭代開始

2　小波核極限學(xué)習(xí)機(jī)[7]

極限學(xué)習(xí)機(jī)[8]是一種單隱層前饋神經(jīng)網(wǎng)絡(luò)的特殊訓(xùn)練算法，該算法最大的特點(diǎn)在于大大減少了需要迭代計(jì)算求解的權(quán)值和偏置個(gè)數(shù)，將原兩層權(quán)值的迭代求解問題轉(zhuǎn)換為通過最小二乘法求解最后一層的輸出權(quán)值。因而不僅大大加快了訓(xùn)練速度，而且非但沒有喪失泛化能力反而有所提升。單隱層前饋神經(jīng)網(wǎng)絡(luò)的輸出層公式為

=h(x)β

(5)

這里h(x)=[h1(x),h2(x),…,hL(x)]T和β=[β1,β2,…,βL]T分別是隱藏層激活函數(shù)向量和隱藏層輸出權(quán)值，ai和bi分別是輸入層權(quán)值和偏置，β由正則化的Moore-Penorse廣義逆給出為

βt=HT(I/C+HHT)-1T

(6)

上式中關(guān)于N個(gè)輸入樣本[x1,x2,…,xN]的隱層輸出矩陣H定義為

(7)

其中L是隱層神經(jīng)元數(shù)，而關(guān)于N個(gè)輸入樣本的訓(xùn)練數(shù)據(jù)目標(biāo)矩陣則由下式定義

(8)

其中m是輸出維數(shù)。由以上公式可以得出ELM算法的訓(xùn)練流程圖如下：

圖1　極限學(xué)習(xí)機(jī)訓(xùn)練流程

由圖1知，原始極限學(xué)習(xí)機(jī)算法需要隨機(jī)選取輸入層權(quán)值ai和輸入層偏置bi,手動(dòng)選取隱層神經(jīng)元數(shù)L，給算法帶來很大不確定性。仔細(xì)分析可以發(fā)現(xiàn)，這些參數(shù)的作用僅在于設(shè)定隱層特征映射h(x)的具體形式，而實(shí)際上由式(5)和式(6)可知，在輸出函數(shù)中無需知道h(x)的具體形式,只需知道內(nèi)積h(xi)·h(xj)的形式即可，如下式所示：

fL(x)=h(x)β=[K(x,x1),…,K(x,xN)](1/C+

G)-1T

(9)

(10)

小波核函數(shù)中的參數(shù)ω和a將通過文中的改進(jìn)Fisher準(zhǔn)則采用模擬退火粒子群算法予以聯(lián)合優(yōu)化。

3　改進(jìn)的Fisher準(zhǔn)則

(11)

文中Fisher準(zhǔn)則用于指導(dǎo)核函數(shù)參數(shù)的選擇，也就是用來挑選隱層非線性特征映射φ，為此需要將核空間的類間離散度和類內(nèi)離散度用核函數(shù)表達(dá)出來，這是可以辦到的，結(jié)果為

(12)

(13)

(14)

以上S2/S1的表達(dá)式中，xir代表第i類的第r個(gè)樣本點(diǎn)。

實(shí)際運(yùn)用時(shí)為了計(jì)算方便和精度考慮，可以等價(jià)地最小化S1-S2。

4　模擬退火粒子群算法

粒子群算法(PSO)[10]是一種在復(fù)雜解空間中通過粒子群的運(yùn)動(dòng)近似搜索最優(yōu)解的啟發(fā)式算法，具有易于實(shí)現(xiàn)和收斂速度快的優(yōu)點(diǎn)，但在進(jìn)化后期容易出現(xiàn)粒子趨同，收斂停滯，優(yōu)化精度不理想等現(xiàn)象?？紤]到模擬退火算法具有較強(qiáng)的概率突變跳出局部最優(yōu)解的能力，提出一種結(jié)合模擬退火和粒子群算法的混合優(yōu)化算法。基于模擬退火的粒子群算法的步驟如下

① 隨機(jī)初始化種群中各個(gè)微粒的位置xi和速度vi

② 初始化pi為粒子i當(dāng)前最優(yōu)適應(yīng)度值對(duì)應(yīng)的位置，初始化pg為當(dāng)前所有粒子中最優(yōu)個(gè)體的位置

③ 確定初始溫度，一般采用如下初始溫度和退溫方式：

T(t+1)=λT(t)T(0)=f(pg)/ln5

(15)

④ 按照下列公式更新各個(gè)微粒的速度和位置：

vij(t+1)=φ{(diào)vij(t)+c1r1[pij-xij(t)]+

c2r2[pgj-xij(t)]}

(16)

xij(t+1)=xij(t)+vij(t+1)

(17)

⑥ 按照各微粒的新位置更新個(gè)體最優(yōu)位置pi和群體最優(yōu)位置pg

⑦ 在同一溫度下，步驟④～⑥重復(fù)L次

⑧ 更新溫度

⑨ 若達(dá)到最大迭代次數(shù)，則停止搜索，輸出群體最優(yōu)值，否則轉(zhuǎn)④

5　故障特征提取

距離相關(guān)系數(shù)[11]是傳統(tǒng)Pearson線性相關(guān)系數(shù)的一種改進(jìn)，它基于兩個(gè)隨機(jī)向量X和Y的聯(lián)合分布fXY(x,y)和各自分布的乘積fX(x)fY(y)之間的函數(shù)距離，可以度量?jī)蓚€(gè)隨機(jī)向量X和Y之間的非線性相關(guān)關(guān)系。隨機(jī)向量X和Y的樣本Xk和Xk，(k=1，2，…，n)之間的經(jīng)驗(yàn)距離相關(guān)系數(shù)定義為

(18)

上式中若分母為0，則dCor(X,Y)=0，其中經(jīng)驗(yàn)距離協(xié)方差dCov(X,Y)定義為

(19)

和式中的有關(guān)項(xiàng)定義為

6　試驗(yàn)數(shù)據(jù)分析

文中所用故障軸承振動(dòng)加速度數(shù)據(jù)來自美國(guó)Case Western Reserve University，數(shù)據(jù)的采樣頻率為12 kHz,軸承轉(zhuǎn)速為1 797 r/min，對(duì)應(yīng)旋轉(zhuǎn)頻率30 Hz。軸承損傷點(diǎn)是用電火花機(jī)人工加工形成的，其直徑為0.017 8 cm。測(cè)試軸承是SKF6205型深球溝軸承的驅(qū)動(dòng)端軸承。采集了外圈，內(nèi)圈，滾動(dòng)體故障和正常情況共4種信號(hào)。每種信號(hào)取長(zhǎng)度為2 048的一段作為樣本進(jìn)行特征提取。

以內(nèi)圈信號(hào)為例，首先基于任意兩個(gè)模態(tài)分量的中心頻率不能太過接近，即最小中心頻率間隔不能小于最大中心頻率間隔的20%定出VMD算法的模態(tài)數(shù)K=4。K=2～5的情況下根據(jù)采樣頻率可計(jì)算出內(nèi)圈故障信號(hào)各模態(tài)中心頻率如表1所示，由表1可見當(dāng)K=5時(shí)最后兩個(gè)模態(tài)的中心頻率較為接近，可能發(fā)生模態(tài)混淆，故取K=4。

VMD的二次懲罰參數(shù)α的最優(yōu)值可以通過信噪比和模態(tài)帶寬方面的考慮求得。首先來考慮信噪比，這里采用的降噪方案基于上節(jié)的距離相關(guān)系數(shù)，具體方法是求出各個(gè)模態(tài)分量同原信號(hào)的距離相關(guān)系數(shù)，舍去小于最大相關(guān)系數(shù)25%的較小相關(guān)系數(shù)對(duì)應(yīng)的IMF分量，剩余的IMF分量求和重構(gòu)即完成降噪過程。例如表2給出了內(nèi)圈故障信號(hào)在α=2 000，K=4，τ=0的情況下各模態(tài)同原信號(hào)的距離相關(guān)系數(shù)dCor值。

表1　內(nèi)圈故障信號(hào)各IMF中心頻率

表2　內(nèi)圈故障信號(hào)各IMF同原 信號(hào)距離相關(guān)系數(shù)

由于表中最小的相關(guān)系數(shù)大于最大相關(guān)系數(shù)的25%，所以沒有模態(tài)被舍去，即便如此，這4個(gè)模態(tài)重構(gòu)求和后并不精確等于原信號(hào)。這是因?yàn)橐匀ピ攵侵貥?gòu)原信號(hào)為目的，為了將噪聲排除在重構(gòu)和外，已將施加重構(gòu)約束的拉氏乘子λ(t)的更新系數(shù)τ的值設(shè)為0(見VMD算法流程第5步)。根據(jù)重構(gòu)前后信號(hào)u(t)和f(t)，可以算出信噪比

(20)

計(jì)算發(fā)現(xiàn)信噪比隨著懲罰參數(shù)增大而減小且遞減速率漸趨緩慢(圖2)。由于VMD可以看成自適應(yīng)的維納濾波器組，因而信噪比跟模態(tài)分量帶寬一樣均隨懲罰參數(shù)變小而變大，兩者變化趨勢(shì)是一致的。

圖2　不同懲罰參數(shù)α下VMD濾波后信噪比

由于參數(shù)α越小信噪比越大且各模態(tài)的帶寬越大，因此應(yīng)當(dāng)選取各模態(tài)無頻帶交疊時(shí)的最小α值。各模態(tài)頻帶不發(fā)生交疊的條件是任兩模態(tài)中心頻率間隔大于兩模態(tài)帶寬之和的一半。由此可以計(jì)算出滿足上述條件的最小α取值為300。

獲得VMD算法關(guān)鍵參數(shù)K和α之后，將信號(hào)做VMD分解,分解所得IMF分量的時(shí)域波形如圖3所示。圖3中模態(tài)的沖擊成分以200個(gè)采樣點(diǎn)為周期，對(duì)應(yīng)頻率為60 Hz,這不是故障的特征頻率而是軸承旋轉(zhuǎn)頻率的倍頻，所以難以從時(shí)域信號(hào)判斷故障類型。隨后提取各模態(tài)的峭度，模態(tài)中心頻率和同原信號(hào)的距離相關(guān)系數(shù)作為輸入小波核極限學(xué)習(xí)機(jī)的特征向量。核極限學(xué)習(xí)機(jī)的核參數(shù)ω和a使用核空間的改進(jìn)Fisher準(zhǔn)則予以優(yōu)化。求解此優(yōu)化問題用到模擬退火粒子群算法，算法參數(shù)為粒子數(shù)N=40，學(xué)習(xí)因子c1=2.05，c2=2.05，收縮因子φ=0.729，最大迭代次數(shù)M=100，退火常數(shù)λ=0.5。優(yōu)化計(jì)算得到的小波核參數(shù)的最優(yōu)值為ω=1.232，a=34.69。小波核函數(shù)指數(shù)因子中的參數(shù)a的增加會(huì)增加該指數(shù)因子，從而增大了核空間的內(nèi)積，直觀上提高了核空間的線性區(qū)分度。

為說明文中提出的算法的效果，運(yùn)行3組對(duì)比試驗(yàn)：(1)同樣的粒子群算法參數(shù)但使用原始Fisher準(zhǔn)則，優(yōu)化得到的小波核參數(shù)的最優(yōu)值為ω=0.002 55，a=4.275。(2)將同一信號(hào)進(jìn)行EMD分解，提取分解后同原信號(hào)距離相關(guān)系數(shù)最大的前4個(gè)IMF分量的峭度，平均瞬時(shí)頻率以及距離相關(guān)系數(shù)后輸入?yún)?shù)同為ω=1.232，a=34.69的核極限學(xué)習(xí)機(jī)中。(3)同一組信號(hào)VMD分解，但是用普通的相關(guān)系數(shù)而非距離相關(guān)系數(shù)作為特征。

全部4組試驗(yàn)條件下測(cè)試集準(zhǔn)確率隨正則化參數(shù)變化的曲線如圖4所示，由圖4看出，4組試驗(yàn)條件下測(cè)試集準(zhǔn)確率均隨著正則化參數(shù)C增大而增大，且增幅漸趨平緩。使用VMD分解且用改進(jìn)的Fisher準(zhǔn)則優(yōu)化的小波核極限學(xué)習(xí)機(jī)的分類準(zhǔn)確率最高，說明Fisher準(zhǔn)則改進(jìn)后確實(shí)起到了進(jìn)一步改良非線性特征映射，使得映射后不同類別的樣本點(diǎn)在核空間中更為有效分離的作用。造成這種現(xiàn)象的原因是改進(jìn)的Fisher準(zhǔn)則通過類間最小差異度和類內(nèi)最大差異度定義的類內(nèi)和類間離散度比原始Fisher準(zhǔn)則通過樣本均值定義的類內(nèi)和類間離散度更為合理，考慮了類中樣本點(diǎn)分布的極差，顧及了分布的全局性質(zhì)。而以樣本均值為代表計(jì)算的離散度體現(xiàn)的是個(gè)別樣本點(diǎn)的局部性質(zhì)。因而改進(jìn)的Fisher準(zhǔn)則能夠在核空間獲得更為優(yōu)良的線性可分性。使用普通相關(guān)系數(shù)而非距離相關(guān)系數(shù)作為特征的試驗(yàn)組準(zhǔn)確率排在最后，這是因?yàn)榫嚯x相關(guān)系數(shù)比起通常的只能度量線性關(guān)系的線性相關(guān)系數(shù)更好地反映了隨機(jī)變量間的非線性依賴關(guān)系，是一種更為有效的刻畫隨機(jī)變量間獨(dú)立性的統(tǒng)計(jì)量。

圖3　α=300 K=4時(shí)VMD 4個(gè)模態(tài)的時(shí)域圖

圖4　測(cè)試集準(zhǔn)確率隨正則化 參數(shù)C變化曲線

為了研究文中所提算法對(duì)噪聲的魯棒性，對(duì)4種工況振動(dòng)信號(hào)加入信噪比為5 dB的高斯白噪聲后再進(jìn)行VMD分解，特征抽取和狀態(tài)識(shí)別。表3和表4分別為加入噪聲前后的混淆矩陣數(shù)據(jù)。

表3　算法在添加高斯白噪聲前混淆矩陣

表4　算法在添加高斯白噪聲后混淆矩陣

由表4可知，加入高斯白噪聲后，文中算法對(duì)正常工況，外圈故障，內(nèi)圈故障和滾動(dòng)體故障的識(shí)別率分別為97.5%，95%，100%和95%,平均識(shí)別率仍然達(dá)到了96.875%，表明VMD分解方法具有良好的噪聲魯棒性。作為對(duì)比，表5和表6是采用EMD分解方法的相應(yīng)混淆矩陣，可以看到噪聲的加入嚴(yán)重影響了EMD的分類準(zhǔn)確率。EMD分解方法不能很好地應(yīng)用于實(shí)際中富含噪聲的真實(shí)信號(hào)。

表5　采用EMD分解的算法在添加 高斯白噪聲前混淆矩陣

表6　采用EMD分解的算法在添加 高斯白噪聲后混淆矩陣

為了說明核極限學(xué)習(xí)機(jī)相對(duì)同樣基于核方法的支持向量機(jī)學(xué)習(xí)算法的優(yōu)越性，表7給出了基于VMD分解和SVM分類算法的混淆矩陣數(shù)據(jù)。這里在4類數(shù)據(jù)的每2類之間訓(xùn)練一個(gè)Morlet小波核SVM分類器，核參數(shù)采用改進(jìn)的Fisher準(zhǔn)則予以優(yōu)化，樣本最后的分類結(jié)果是在2類分類中“勝出”數(shù)最大的類別。

從表7可知SVM平均識(shí)別率為95.625%，同表3對(duì)比可知在準(zhǔn)確率方面略低于核極限學(xué)習(xí)機(jī)。表8給出了不同的訓(xùn)練樣本數(shù)占總樣本數(shù)百分比情況下基于交叉驗(yàn)證和基于改進(jìn)的核空間Fisher準(zhǔn)則的SVM 核參數(shù)尋優(yōu)耗時(shí)和準(zhǔn)確率數(shù)據(jù)。

表7　采用VMD分解的SVM算法的混淆矩陣

表8　基于交叉驗(yàn)證和基于Fisher準(zhǔn)則的 SVM核參數(shù)優(yōu)化的準(zhǔn)確率和耗時(shí)數(shù)據(jù)

從表8看出改進(jìn)的核空間的Fisher準(zhǔn)則相比交叉驗(yàn)證的參數(shù)優(yōu)化方法而言在準(zhǔn)確率上相差不多但是參數(shù)優(yōu)化的耗時(shí)有明顯下降。所以就SVM算法來說，改進(jìn)的Fisher準(zhǔn)則不失為一種快速的優(yōu)選參數(shù)的手段。對(duì)于核極限學(xué)習(xí)算法而言，由于核極限學(xué)習(xí)機(jī)的訓(xùn)練速度遠(yuǎn)快于SVM，因而改進(jìn)的Fisher準(zhǔn)則在耗時(shí)上的優(yōu)勢(shì)無法發(fā)揮，但是對(duì)于小樣本訓(xùn)練集的情形，交叉驗(yàn)證方法不再適用，而改進(jìn)的Fisher準(zhǔn)則仍然可以在小樣本情形下指導(dǎo)SVM和核極限學(xué)習(xí)機(jī)的核參數(shù)優(yōu)化。

7　結(jié)　論

針對(duì)核極限學(xué)習(xí)機(jī)在軸承故障診斷中遇到的參數(shù)選取問題提出了一種采用改進(jìn)后核空間中的Fisher準(zhǔn)則進(jìn)行參數(shù)優(yōu)選的方案，針對(duì)振動(dòng)信號(hào)信噪比低的問題提出了一種基于變分模態(tài)分解和距離相關(guān)系數(shù)的濾波和特征提取方法。通過實(shí)驗(yàn)數(shù)據(jù)的分析，得出以下結(jié)論：

(1) 距離相關(guān)系數(shù)能夠很好地度量隨機(jī)變量間的非線性相關(guān)關(guān)系，在濾波和特征提取方面的效果要優(yōu)于線性相關(guān)系數(shù)。

(2) VMD分解方法有效避免了EMD分解的模態(tài)混疊現(xiàn)象，而且對(duì)噪聲有較高魯棒性，在軸承故障信號(hào)分析領(lǐng)域有廣闊的應(yīng)用前景。

(3) 提出的改進(jìn)的Fisher準(zhǔn)則相比原始的Fisher準(zhǔn)則提升了分類算法的泛化性能，相比交叉驗(yàn)證法降低了參數(shù)優(yōu)化耗時(shí)，適用于SVM算法的參數(shù)優(yōu)化以及小樣本情形下核極限學(xué)習(xí)算法的參數(shù)優(yōu)化。