基于EEMD和IHHO-LSSVM滾動(dòng)軸承剩余壽命預(yù)測(cè)方法的研究

胡豁然，李亞莎

(三峽大學(xué)電氣與新能源學(xué)院，湖北 宜昌 443002)

滾動(dòng)軸承是電力設(shè)備中應(yīng)用最廣的一種典型重要零部件，滾動(dòng)軸承的運(yùn)行狀況直接影響電氣設(shè)備的安全、可靠運(yùn)行，壽命是衡量軸承的重要指標(biāo)[1-2]。隨著機(jī)器學(xué)習(xí)的發(fā)展，基于數(shù)據(jù)驅(qū)動(dòng)的軸承壽命預(yù)測(cè)逐漸成為國(guó)內(nèi)外重點(diǎn)研究方向。但軸承原始信號(hào)中不可避免會(huì)有大量噪聲，影響了軸承壽命預(yù)測(cè)的準(zhǔn)確性，因此需要對(duì)原始信號(hào)進(jìn)行處理，文獻(xiàn)[3]針對(duì)原始信號(hào)中的大量噪聲，采用EEMD對(duì)原始數(shù)據(jù)進(jìn)行分解、重構(gòu)。為了進(jìn)行退化特征的提取和壽命預(yù)測(cè)，文獻(xiàn)[4]提出了一種基于多頻率尺度模糊熵和ELM的滾動(dòng)軸承剩余壽命預(yù)測(cè)；文獻(xiàn)[5-6]提出了一種基于改進(jìn)HMM和Pearson相似度分析的滾動(dòng)軸承自適應(yīng)壽命預(yù)測(cè)方法；文獻(xiàn)[7-8]提出了一種基于小波包能量譜和改進(jìn)FOA-GRNN的軸承壽命預(yù)測(cè)；文獻(xiàn)[9-10]提出了基于t-SNE和LSTM對(duì)滾動(dòng)軸承剩余壽命的研究。上述方法為軸承退化特征的提取和壽命預(yù)測(cè)提供了有效借鑒，但仍有考慮不周的地方，例如對(duì)原始信號(hào)進(jìn)行降噪的方法仍存在不足，壽命預(yù)測(cè)只能預(yù)測(cè)短期的退化趨勢(shì)且預(yù)測(cè)結(jié)果較為保守。

針對(duì)以上問(wèn)題，提出一套新的滾動(dòng)軸承壽命預(yù)測(cè)方法，先用EEMD對(duì)原始信號(hào)進(jìn)行分解，然后根據(jù)相關(guān)性系數(shù)指標(biāo)對(duì)IMF進(jìn)行篩選并重構(gòu)信號(hào)，并對(duì)特征指標(biāo)集進(jìn)行核主成分分析，選取累計(jì)貢獻(xiàn)率符合要求的核主成分進(jìn)行分析，建立IHHO-LSSVM模型，采用改進(jìn)的哈里斯鷹算法優(yōu)化LSSVM參數(shù)，提高滾動(dòng)軸承剩余壽命預(yù)測(cè)模型的性能和精度。

1 軸承性能退化評(píng)估模型構(gòu)建

1.1 EEMD原理

WU[11]等提出的EEMD是EMD的一種改進(jìn)算法，EMD在使用過(guò)程中存在因噪聲或脈沖產(chǎn)生的模態(tài)混疊現(xiàn)象，為了解決該問(wèn)題，EEMD利用白噪聲的頻譜均勻性，通過(guò)每次給信號(hào)加入不同幅值的高斯白噪聲改變信號(hào)極值點(diǎn)特性，形成基本無(wú)擬合誤差的上、下包絡(luò)線，從而消除模態(tài)混疊現(xiàn)象。

EEMD具體分解過(guò)程如下。

a.在信號(hào)中加入一定幅值的白噪聲序列得到新的信號(hào)x(t)。

b.對(duì)新信號(hào)x(t)進(jìn)行EEMD分解，獲得k個(gè)IMF分量cj(t)和一個(gè)余量r(t)，(j=1,2，…，k)。

(1)

c.每次加入相同振幅的不同白噪聲序列，依次重復(fù)步驟a、bN次。

d.對(duì)分解得到的各個(gè)IMF分量均值作為最終結(jié)果，最終得到的EEMD分解后的IMF為

(2)

e.選取合適的IMF分量進(jìn)行重構(gòu)數(shù)據(jù)。

1.2 核主成分分析

核主成分分析屬于非線性主元分析方法，采用非線性的方法提取主成分。

首先對(duì)時(shí)域、頻域、小波包能量譜所構(gòu)成的樣本進(jìn)行非線性變化，將其映射到一個(gè)高維特征空間F，協(xié)方差矩陣為

(3)

式中：xj為輸入空間的樣本(j= 1,2,…,m)；Φ(x)為在特征空間中的樣本點(diǎn)；C為協(xié)方差矩陣；m為輸入空間的樣本個(gè)數(shù)。

以上矩陣的特征向量和特征值滿足以下條件：

λV-CV=0

(4)

式中：V為協(xié)方差矩陣特征向量；λ為協(xié)方差矩陣特征值。

引入非線性函數(shù)Φ(xj)，可得：

λΦ(xj)V-Φ(xj)CV=0

(5)

式 (5)中得特征向量V可由Φ(x)線性表示，即：

(6)

將式(4)和式(6)代入式(5)，引入核函數(shù)Kij=K(xi,xj)=Φ(xi)Φ(xj),簡(jiǎn)化后有：

mλα-Kα=0

(7)

式中：K為核矩陣；α為核矩陣K的特征向量。

對(duì)于任意樣本，在特征空間F中主元Φ(x)上的投影為

(8)

1.3 相關(guān)系數(shù)指標(biāo)

分解結(jié)果中IMF分量容易產(chǎn)生虛假分量，因此引入相關(guān)系數(shù)辨別剔除虛假分量，設(shè)y1=(y11+y12+…+y1n)T，y2=(y21+y22+…+y2n)T為2個(gè)數(shù)據(jù)向量。則y1(t)、y2(t)的相關(guān)系數(shù)定義為

(9)

通過(guò)式(9)可知，計(jì)算各IMF與原始數(shù)據(jù)之間的相關(guān)系數(shù)，絕對(duì)值越大表示在該IMF分量中所包含的原始信息量越大。峭度指標(biāo)是一個(gè)無(wú)量綱參數(shù)，它對(duì)于振動(dòng)沖擊信號(hào)非常敏感，通過(guò)峭度指標(biāo)可以反映軸承振動(dòng)信號(hào)中的沖擊成分，因此適用于對(duì)機(jī)械系統(tǒng)進(jìn)行故障診斷。當(dāng)滾動(dòng)軸承零部件在正常運(yùn)行的過(guò)程中，對(duì)采集后的振動(dòng)信號(hào)進(jìn)行峭度值的計(jì)算。

經(jīng)過(guò)峭度指標(biāo)和有效系數(shù)的辨別篩選，選取合適的IMF分量重構(gòu)信號(hào)。在重構(gòu)信號(hào)中選取合適的特征值能夠反映軸承的退化性能評(píng)估指標(biāo)。

2 軸承壽命預(yù)測(cè)模型的建立

2.1 改進(jìn)SVM預(yù)測(cè)模型

在支持向量機(jī)的基礎(chǔ)上引入最小二乘線性理論就是最小二乘支持向量機(jī)，它是SVM的改進(jìn)算法，具有和標(biāo)準(zhǔn)支持向量機(jī)互不相同的約束條件，它是用等式約束代替SVM中的不等式約束，將不等式約束的二次規(guī)劃問(wèn)題轉(zhuǎn)化成線性矩陣求解問(wèn)題，與標(biāo)準(zhǔn)的支持向量機(jī)相比，改進(jìn)的LSSVM求解速度得到改善，精度得到提高。LSSVM的函數(shù)估計(jì)問(wèn)題如下。

(10)

此時(shí)約束條件：

s.t.yi-ωΦ(xj)-b=0

(11)

式中：ω為權(quán)向量；γ為懲罰參數(shù)；b為偏差向量；ei為松弛變量，表示第i個(gè)數(shù)據(jù)的預(yù)測(cè)輸出和實(shí)際輸出的誤差值。

式(11)和式(12)對(duì)應(yīng)的拉格朗日函數(shù)為

b+ei-yi]

(12)

式中：αi為L(zhǎng)agrange乘子，依據(jù)KKT條件進(jìn)行求解，即對(duì)式(12)中的ω、b、e、αi求偏導(dǎo)數(shù)，并令其偏導(dǎo)數(shù)為零。

消掉變量ω、e，得到方程組如下。

(13)

yn]T。

σ為核函數(shù)寬密度。

最終得到LSSVM預(yù)測(cè)回歸模型為

(14)

建立LSSVM輸入x=(xn-m,xn-m+1,…,xn-1)與輸出y=(xn)一種非線性的映射關(guān)系：Rm→R,m為輸入函數(shù)，然后得到LSSVM的訓(xùn)練樣本和測(cè)試樣本對(duì)為

式中：l為訓(xùn)練樣本的對(duì)數(shù)；s為訓(xùn)練數(shù)據(jù)的起點(diǎn)；n為預(yù)測(cè)點(diǎn)數(shù)。輸入為Xtrain，輸出為Ytrain，從而訓(xùn)練出預(yù)測(cè)模型，求出LSSVM中預(yù)測(cè)函數(shù)(14)的具體數(shù)學(xué)表達(dá)式，將測(cè)試樣本Xtest帶入到式(14)中可以得到對(duì)應(yīng)的n個(gè)預(yù)測(cè)值。

2.2 改進(jìn)哈里斯鷹優(yōu)化算法(IHHO)

2.2.1 HHO優(yōu)化算法

HHO算法[12]是Heidari在2019年受哈里斯鷹捕食獵物啟發(fā)所提出的一種新型智能優(yōu)化算法，主要由搜索階段、搜索與開發(fā)轉(zhuǎn)化階段和開發(fā)階段組成。

哈里斯鷹棲息在搜索空間[lb,ub]觀察某位置圍捕獵物；在迭代時(shí)根據(jù)其他個(gè)體與獵物的位置以概率p選擇更新。

(15)

式中：X(t+1)和X(t)表示第t+1次和t次迭代時(shí)個(gè)體的位置；Xt為獵物時(shí)的位置；Xd為個(gè)體的位置；ri(i=1,…,4)和p均為[0,1]內(nèi)的一個(gè)隨機(jī)數(shù)，Xm(t)為t次迭代時(shí)的N只個(gè)體表達(dá)平均位置。

HHO算法根據(jù)獵物的逃逸能量在搜索和不同的開發(fā)行為之間轉(zhuǎn)換，用能量E線性動(dòng)態(tài)遞減調(diào)控算法的全局搜索和局部開發(fā)性能，其定義式為

(16)

式中：E0為獵物的初始能量，為[-1,1]之間的隨機(jī)數(shù)；T為最大迭代數(shù)。當(dāng)|E|≥1，HHO進(jìn)入全局搜索階段；當(dāng)|E|<1時(shí)，HHO算法則進(jìn)入局部開發(fā)過(guò)程。

在HHO算法中，通過(guò)因子β∈[0,1]來(lái)描述獵物是否逃脫成功：當(dāng)β<0.5時(shí)表示獵物逃脫成功，反之則失效：同時(shí)哈里斯鷹根據(jù)獵物能量E與0.5的相對(duì)大小來(lái)采取不同的圍捕策略，包含軟圍攻、硬圍攻、漸進(jìn)式快速俯沖圍捕和漸進(jìn)式快速俯沖硬捕4種情形。

Case1：軟圍攻，當(dāng)β≥0.5且|E|≥0.5，表示獵物的能量E充足，并且以隨機(jī)跳躍的方式逃脫圍捕，但最終失敗并被捕獲，哈里斯鷹的位置更新為

X(t+1)=ΔXt(t)-E|JXt(t)-Xt(t)|

(17)

式中：ΔXt為第t次迭代時(shí)鷹與獵物的位置偏差；J=2(1-r5) 為獵物逃跑時(shí)隨機(jī)跳躍能力，r5∈(0,1)的一個(gè)隨機(jī)數(shù)。

Case2：硬圍攻，當(dāng)β≥0.5且 |E|<0.5，表示獵物的能量E較低而被鷹直接捕獲，位置更新為

X(t+1)=Xt(t)-|EΔX(t)|

(18)

Case3：漸進(jìn)式快速俯沖圍捕，當(dāng)β<0.5且|E|≥0.5，表示獵物能量E充足可保證逃脫成功，但鷹以最優(yōu)方向俯沖并軟圍捕獵物，其位置更新為

(19)

Y=Xt(t)-E|JXt(t)-Xt(t)|

(20)

Z=Y+S+LF(D)

(21)

式中:D為問(wèn)題維度；S為D維隨機(jī)向量；LF為L(zhǎng)evy函數(shù)。

Case4：漸進(jìn)式快速俯沖硬捕，當(dāng)β<0.5且|E|<0.5，表示獵物能量E較低，而鷹俯沖硬圍捕獵物，其位置更新式為

(22)

2.2.2 改進(jìn)搜索階段

為了促進(jìn)搜索階段的全局搜索，式(8)哈里斯鷹移動(dòng)取決于獵物和其他個(gè)體的位置，而對(duì)于IHHO，第1種策略中(p<0.5)，其他個(gè)體鷹捕食獵物的移動(dòng)被認(rèn)為是模仿哈里斯鷹曲折運(yùn)動(dòng)飛行，采用Levy飛行模擬。而在第2種策略中(p>0.5)，哈里斯鷹個(gè)體移動(dòng)在內(nèi)部是隨機(jī)的。因此為了改善這一階段，將哈里斯鷹與其他3只隨機(jī)鷹初始值的平均位置來(lái)升級(jí)更新位置。

(23)

式中：Xj、Xk、Xl為哈里斯鷹隨機(jī)初始選擇的位置；Xm(t)為獵物的位置；LF為L(zhǎng)evy函數(shù)。當(dāng)概率因子p<0.5，哈里斯鷹的位置基于Levy飛行；當(dāng)概率因子p≥0.5，哈里斯鷹遵循式(23)中第2行公式，隨機(jī)更改內(nèi)部的位置。

2.2.3 改進(jìn)能量調(diào)控機(jī)制

在HHO算法中，獵物能量E的大小反映哈里斯鷹對(duì)模型最優(yōu)解的搜尋捕獲能力，并不能有效描述哈里斯鷹經(jīng)多輪協(xié)同圍捕獵物的情況[13]，本文提出一種能量周期性調(diào)控機(jī)制，實(shí)現(xiàn)了“全局+局部”的尋優(yōu)搜索。以余弦函數(shù)表示獵物能量E的周期遞變性，定義式為

(24)

式中：k=0,1,2,…，為獵物能量E的周期遞減數(shù)。

3 IHHO-LSSVM模型流程圖

結(jié)合IHHO 算法和改進(jìn)LSSVM預(yù)測(cè)模型，其IHHO-LSSVM的流程如圖1所示。

圖1 基于IHHO-LSSVM的流程

選取KPCA降維后的軸承退化特征作為輸入向量，軸承剩余壽命值作為輸出向量，以預(yù)測(cè)壽命均方誤差(MSE)最小化為IHHO尋優(yōu)目標(biāo)，以核參數(shù)δ和正規(guī)化參數(shù)C為優(yōu)化參數(shù)。

(25)

4 基于軸承全壽命預(yù)測(cè)的驗(yàn)證

4.1 試驗(yàn)裝置

采用來(lái)自美國(guó)凱西西儲(chǔ)大學(xué)電氣工程實(shí)驗(yàn)室的滾動(dòng)軸承全壽命周期加速軸承性能退化數(shù)據(jù)。試驗(yàn)平臺(tái)如圖2所示。

圖2 滾動(dòng)軸承全壽命試驗(yàn)平臺(tái)

軸承試驗(yàn)臺(tái)上的主軸上安裝了4個(gè)Rexnord-ZA-2115型號(hào)的雙列滾動(dòng)軸承。試驗(yàn)過(guò)程中，每個(gè)軸承將承受26 670 N的徑向載荷，每個(gè)軸承的X和Y方向各安裝了1個(gè)加速度傳感器，采樣頻率為20 kHz，每隔10 min采集1次滾動(dòng)軸承振動(dòng)信號(hào)，一共采集984個(gè)樣本，每個(gè)樣本長(zhǎng)度20 480個(gè)點(diǎn)。共20 152 320個(gè)數(shù)據(jù)，軸承運(yùn)行了1周，直到軸承退化失效，如圖3所示。

圖3 軸承原始信號(hào)

4.2 基于EEMD的軸承信號(hào)降噪

采集的信號(hào)中一般含有噪聲，不處理將會(huì)對(duì)軸承的性能退化指標(biāo)產(chǎn)生較大影響。文中采用EEMD方法分解、篩選、重構(gòu)，得到重構(gòu)信號(hào)(每隔10個(gè)數(shù)據(jù)抽取1個(gè)，每個(gè)樣本長(zhǎng)度2048個(gè)點(diǎn))。具體操作：對(duì)原始振動(dòng)信號(hào)進(jìn)行EEMD計(jì)算，其中EEMD算法設(shè)置為Rstd=0.2，添加高斯白噪聲次數(shù)為1000，分解得到如圖4所示的21個(gè)IMF向量，計(jì)算21個(gè)IMF向量的峭度值和相關(guān)系數(shù)。

從圖4可知，IMF1-IMF12振動(dòng)沖擊信號(hào)較多，而后續(xù)IMF缺少信息，利用相關(guān)系數(shù)對(duì)12個(gè)IMF分量進(jìn)行篩選，前12個(gè)相關(guān)系數(shù)和峭度值的計(jì)算值，如表1所示。

圖4 各個(gè)IMF分量

表1 前12個(gè)相關(guān)系數(shù)和峭度值的計(jì)算

利用相關(guān)系數(shù)和峭度值的計(jì)算，峭度值偏離越大，說(shuō)明軸承的故障信息越多。因此選取了1、3、4、5、6、7、8、9總計(jì)8個(gè)作為重構(gòu)信號(hào)分量。重構(gòu)信號(hào)中沖擊的成分更明顯，同時(shí)高頻和部分低頻中噪聲被去除。

4.3 驗(yàn)證結(jié)果與分析

針對(duì)軸承時(shí)域，頻域中每個(gè)特征所呈現(xiàn)出的退化能力各不相同，選取能夠較好地反映出軸承在全壽命周期過(guò)程中退化趨勢(shì)的特征，剔除在全壽命周期過(guò)程中幾乎沒(méi)有任何變化的特征。因此，在時(shí)域退化特征中，篩選了均方根值、標(biāo)準(zhǔn)差、峰峰值為時(shí)域特征，均值頻率、標(biāo)準(zhǔn)差為頻域特征?？偣策x取了22個(gè)退化特征，利用核主成分分析方法(KPCA方法)將特征退化指標(biāo)加權(quán)融合，提取KPCA的第一主成分作為性能退化評(píng)估指標(biāo)，既能夠得到軸承性能退化趨勢(shì)，又使維數(shù)得到簡(jiǎn)化。將KPCA第一主成分進(jìn)行標(biāo)準(zhǔn)化處理，結(jié)果如圖5所示。

圖5 KPCA第一主成分

對(duì)比KPCA第一主成分退化趨勢(shì)和時(shí)域、頻域、小波包能量譜可以發(fā)現(xiàn)，KPCA第一主成分在500點(diǎn)(早期故障)后就已經(jīng)呈現(xiàn)優(yōu)良的上升趨勢(shì)，而單一的時(shí)域、頻域、小波包能量譜則不具有這種優(yōu)勢(shì)。因此，利用KPCA加權(quán)融合提取主成分呈現(xiàn)軸承退化性能趨勢(shì)比單一的特征要全面，既有良好的上升趨勢(shì)，又對(duì)軸承故障早期具有敏感性。將KPCA各個(gè)成分的貢獻(xiàn)率從大到小依次排列，選取累計(jì)貢獻(xiàn)率達(dá)到85%的前n個(gè)主成分，如表2所示。

表2 部分KPCA的貢獻(xiàn)率

表2中，前3個(gè)主成分已經(jīng)達(dá)到85%，故選取軸承退化期的第1、2、3個(gè)主成分和對(duì)應(yīng)的剩余壽命組成訓(xùn)練樣本，建立IHHO-LSSVM壽命預(yù)測(cè)模型，利用500點(diǎn)后的數(shù)據(jù)訓(xùn)練，等間隔取7個(gè)作為預(yù)測(cè)樣本模型輸入，將LSSVM的核參數(shù)δ和正規(guī)化參數(shù)C通過(guò)改進(jìn)的哈里斯鷹算法來(lái)進(jìn)行優(yōu)化，設(shè)置種群空間維度dim=2，最大迭代次數(shù)maxgen=100,設(shè)置哈里斯鷹種群數(shù)目sizepop=30，將未改進(jìn)的HHO算法與IHHO算法的適應(yīng)度值(適應(yīng)度值選擇均方誤差)曲線迭代進(jìn)行對(duì)比，如圖6所示。

由圖6可知，IHHO算法在收斂精度上提高了許多，為了對(duì)比和評(píng)估預(yù)測(cè)軸承壽命預(yù)測(cè)效果，將未采用EEMD-KPCA處理的IHHO-LSSVM模型和采用EEMD-KPCA處理的LSSVM模型進(jìn)行壽命預(yù)測(cè)對(duì)比。如圖7所示，各個(gè)模型的預(yù)測(cè)誤差對(duì)比如表3所示。

圖6 2種方法尋優(yōu)LSSVM參數(shù)的適應(yīng)度值收斂趨勢(shì)

(a)EEMD-LSSVM模型

(b)IHHO-LSSVM模型

(c)EEMD-KPCA和IHHO-LSSVM模型圖7 各個(gè)模型壽命預(yù)測(cè)結(jié)果

表3 各個(gè)模型的預(yù)測(cè)誤差對(duì)比

由圖7可知，縱坐標(biāo)表示每個(gè)樣本對(duì)應(yīng)的剩余壽命時(shí)間序列，橫坐標(biāo)表示軸承的使用時(shí)間，3種模型預(yù)測(cè)的軸承剩余壽命與真實(shí)壽命變化趨勢(shì)一致。對(duì)照?qǐng)D7、表3可知，利用EEMD降噪去除噪聲、分解、重構(gòu)信號(hào)，用KPCA融合前3個(gè)主成分作為軸承性能退化指標(biāo)，最后采用IHHO-LSSVM作為壽命預(yù)測(cè)模型的3類誤差均小于前2種，說(shuō)明EEMD分解重構(gòu)具有消除噪聲的影響；KPCA消除了時(shí)域、頻域等各特征指標(biāo)間的相關(guān)性導(dǎo)致的冗余數(shù)據(jù)；IHHO-LSSVM模型預(yù)測(cè)準(zhǔn)確率對(duì)于其他支持向量機(jī)模型有所提高。

5 結(jié)論

a.采用EEMD-KPCA方法提取軸承性能退化指標(biāo)，用 EEMD降噪重構(gòu)信號(hào)，針對(duì)時(shí)域、頻域及小波包能量譜各特征指標(biāo)間相關(guān)性而導(dǎo)致的數(shù)據(jù)冗余現(xiàn)象，提出采用KPCA提取主成分。該指標(biāo)能夠兼顧軸承早期故障敏感性和全面表征軸承退化情況，為建立預(yù)測(cè)模型提高了良好的基礎(chǔ)。

b.提出了一種改變搜索階段和能量調(diào)控機(jī)制的新型哈里斯鷹群智能算法，與原算法相比，更容易找到LSSVM預(yù)測(cè)模型的核參數(shù)δ和正規(guī)化參數(shù)C的最優(yōu)解。

c.建立一種新的IHHO-LSSVM滾動(dòng)軸承壽命預(yù)測(cè)模型，與EEMD-LSSVM和HHO-LSSVM預(yù)測(cè)模型相比，滾動(dòng)軸承的壽命預(yù)測(cè)精度得到了提高，為滾動(dòng)軸承壽命預(yù)測(cè)提供了一種新的方法，確保電力設(shè)備全壽命周期健康管理，保證生產(chǎn)正常連續(xù)運(yùn)行。