栗子旋，高丙朋

(新疆大學(xué)電氣工程學(xué)院，新疆烏魯木齊 830017)

0 前言

滾動(dòng)軸承作為旋轉(zhuǎn)機(jī)械的關(guān)鍵部件，決定著機(jī)械設(shè)備能否安全平穩(wěn)運(yùn)行。在復(fù)雜的運(yùn)行工況下，隨著服役時(shí)間的增加，軸承性能逐步退化，發(fā)生故障甚至導(dǎo)致失效。對滾動(dòng)軸承的剩余使用壽命(Remaining Useful Life，RUL)進(jìn)行預(yù)測，可以確定最佳維修點(diǎn)，實(shí)現(xiàn)預(yù)測性維護(hù)，減少事故的發(fā)生和經(jīng)濟(jì)損失，提高設(shè)備的可靠性[1]。

基于數(shù)據(jù)驅(qū)動(dòng)的滾動(dòng)軸承剩余壽命預(yù)測技術(shù)的研究包含2個(gè)重要方面：運(yùn)行狀態(tài)的特征提取和剩余壽命預(yù)測模型的建立。特征提取是關(guān)鍵，壽命預(yù)測是目標(biāo)。良好的特征能夠反映軸承的運(yùn)行狀態(tài)和退化趨勢，而壽命預(yù)測模型則能夠根據(jù)退化特征挖掘軸承的退化規(guī)律，實(shí)現(xiàn)準(zhǔn)確地預(yù)測滾動(dòng)軸承的剩余使用壽命。文獻(xiàn)[2]通過雙樹復(fù)小波重構(gòu)原始信號，得到去噪后的信號，提取相應(yīng)的特征，采用BP網(wǎng)絡(luò)對破碎機(jī)轉(zhuǎn)軸進(jìn)行壽命預(yù)測。文獻(xiàn)[3]提取時(shí)域特征描述軸承的退化過程，并采用濾波器進(jìn)行降噪處理，然后采用優(yōu)化算法得到支持向量回歸的參數(shù)，訓(xùn)練多個(gè)模型進(jìn)行軸承壽命預(yù)測。為了能夠更好地描述軸承的退化特征，提取時(shí)域、頻域以及時(shí)頻域特征共同構(gòu)成特征矩陣。然而多維特征容易導(dǎo)致維數(shù)災(zāi)難。為了兼顧全面反映軸承的運(yùn)行狀態(tài)和計(jì)算效率，對特征矩陣進(jìn)行降維處理。JENSSEN[4]提出了核熵成分分析(Kernel Entropy Component Analysis，KECA)，該算法將原始非線性數(shù)據(jù)經(jīng)過核函數(shù)映射到高維特征空間中，在最大保留原始信息的情況下進(jìn)行降維融合。

壽命預(yù)測階段，有學(xué)者從軸承初期開始進(jìn)行壽命預(yù)測[5-6]，但是軸承初始運(yùn)行階段狀態(tài)平穩(wěn)，對它進(jìn)行壽命預(yù)測沒有實(shí)際工程意義，并且非退化數(shù)據(jù)的引入將造成壽命預(yù)測的精度下降。文獻(xiàn)[7]采用退化階段前期的數(shù)據(jù)訓(xùn)練模型，然后預(yù)測軸承的剩余使用壽命。文獻(xiàn)[8]采用時(shí)變卡爾曼濾波進(jìn)行壽命預(yù)測，構(gòu)造基于線性和二次函數(shù)的卡爾曼濾波模型，自動(dòng)適應(yīng)切換濾波模型，以處理軸承的正常運(yùn)行階段和加速退化階段的監(jiān)測數(shù)據(jù)，從而達(dá)到有效預(yù)測剩余壽命的目的。

相關(guān)向量機(jī)(Relevance Vector Machine，RVM)是在貝葉斯框架下發(fā)展起來的支持向量機(jī)的另一種改進(jìn)形式[9]。該模型結(jié)合了貝葉斯原理、自動(dòng)相關(guān)決策機(jī)制和最大似然等理論[10]，能夠處理高維，非線性、小樣本等問題；具有良好的稀疏性、泛化能力和較高的預(yù)測精度，RVM的核函數(shù)不受限于Mercer條件，在核函數(shù)的選擇和構(gòu)建中更為靈活[11]。核函數(shù)在很大程度上決定著RVM模型的預(yù)測性能[12]，因此采用獵食者-獵物優(yōu)化(Hunter-Prey Optimization，HPO)算法去尋優(yōu)得到RVM的核參數(shù)，該方法主要思想結(jié)合了獵食者向著獵物的方向調(diào)整自己的位置，而獵物則是向著安全的地方調(diào)整自己的位置[13]。

綜上，本文作者提出一種基于特征融合和HPO優(yōu)化RVM的滾動(dòng)軸承剩余壽命預(yù)測方法。提取多域特征，利用特征篩選和核熵成分分析相結(jié)合，提取軸承的退化狀態(tài)信息。利用HPO優(yōu)化算法自適應(yīng)地確定相關(guān)向量機(jī)(Relevance Vector Machine，RVM)的最佳核寬度，通過訓(xùn)練軸承得到壽命預(yù)測模型，將測試軸承的退化特征輸入到模型中，實(shí)現(xiàn)軸承的壽命預(yù)測。

1 特征提取和融合

1.1 多域特征提取

1.1.1 時(shí)域特征

時(shí)域特征包括有量綱和量綱一化2種，特征值的大小可以反映設(shè)備的運(yùn)行狀態(tài)，提取均值、方差、均方根、標(biāo)準(zhǔn)差、最大值、最小值、峰峰值、波形指標(biāo)、峭度指標(biāo)等16個(gè)時(shí)域特征，記為F1-F16。

1.1.2 頻域特征

對原始振動(dòng)信號進(jìn)行快速傅里葉變換，然后計(jì)算相關(guān)的頻域特征，包括頻率均值、頻率標(biāo)準(zhǔn)差、頻率重心等能夠反映頻域振動(dòng)能量大小、頻譜分散或集中程度以及頻帶位置變化的13個(gè)頻域特征。記為F17-F29。

1.1.3 時(shí)頻域特征

小波包分解是一種時(shí)頻域分析方法，可以實(shí)現(xiàn)信號的高頻部分和低頻部分同時(shí)分解，采用這種分解方式既不會(huì)造成信息冗余，也不會(huì)造成信息疏漏。計(jì)算3層小波包分解的能量熵、奇異譜熵和8個(gè)頻帶的能量，共10個(gè)特征。記為F30-F39。

1.2 特征篩選

良好的退化特征應(yīng)具有較強(qiáng)的魯棒性和時(shí)間相關(guān)性。提取的多個(gè)特征并不全都能滿足這2個(gè)特性，因此利用魯棒性和相關(guān)性評價(jià)指標(biāo)對特征進(jìn)行篩選。

魯棒性：反映特征對噪聲或異常值的抗干擾能力，表達(dá)式為

(1)

式中：K為樣本個(gè)數(shù)；XT(tk)是經(jīng)過平滑處理后的趨勢項(xiàng)；XR(tk)是波動(dòng)項(xiàng)；X(tk)=XT(tk)+XR(tk)。

相關(guān)性：表示特征與時(shí)間序列的相關(guān)度，采用Spearman相關(guān)系數(shù)，計(jì)算公式為

(2)

di=rankXi-rankTi(i=1,…,K)

(3)

式中：X、T分別表示2個(gè)時(shí)間序列；di代表2個(gè)序列變量在排序后的秩差。

由于單個(gè)指標(biāo)只能夠評價(jià)特征在某個(gè)方面的特性，因此使用自適應(yīng)加權(quán)方法構(gòu)建綜合評價(jià)指標(biāo)：

J=ω1Rob(X)+ω2Cor(X,T)

(4)

式中：ωi為評價(jià)指標(biāo)的權(quán)重。

采用Critic方法自適應(yīng)確定權(quán)重。該方法綜合考慮了指標(biāo)間的差異性和相關(guān)性，如式(5)所示：

(5)

式中：σj為指標(biāo)j的方差；ζij為相關(guān)系數(shù)。然后計(jì)算權(quán)重，如式(6)所示：

(6)

式中：ωj表示指標(biāo)j的權(quán)重。

1.3 核熵成分分析

KECA是基于Renyi二次熵和Parzen窗函數(shù)提出來的。Renyi二次熵?cái)?shù)學(xué)形式如式(7)所示：

(7)

式中：p(x)為樣本集D={x1，x2，…，xN}的概率密度函數(shù)。將式(7)中的積分部分定義為

(8)

式(8)可以改寫為以下等式：

(9)

(10)

式中：K為N×N的核矩陣；I是N×1的向量。對K進(jìn)行特征分解：

K=EDET

(11)

式中：D=diag(λ1,λ2,…,λN)；E是特征向量矩陣，E={e1,e2,…,eN}。將式(11)代入式(9)可以得到：

(12)

式中：λi和ei是按照二次熵的估值降序排列得到的。選擇滿足要求的前k個(gè)特征值與特征向量進(jìn)行投影變換可以得到Φeca：

(13)

2 HPO算法優(yōu)化相關(guān)向量機(jī)

2.1 相關(guān)向量機(jī)

tn=y(xn;w)+εn

(14)

式中：t=(t1,…,tN)；w=(w0,…,wN)；y(xn;w)為非線性函數(shù)；εn為高斯噪聲。假設(shè)tn服從高斯分布：

p(t|x)=N(t|y(x),σ2)

(15)

由貝葉斯定理，數(shù)據(jù)集的似然估計(jì)可以寫成：

(16)

式中：Φ是核矩陣；Φnm=K(xn,xm-1)，Φn1=1。

根據(jù)貝葉斯規(guī)則獲得權(quán)重的后驗(yàn)概率：

p(w|t,α,σ2)=(2π)-(N+1)/2|Σ|-1/2·

(17)

Σ=(ΦTBΦ+A)-1

(18)

μ=ΣTBt

(19)

式中：A=diag(α0，α1，…，αN)；B=σ-2IN；Σ為后驗(yàn)方差；μ為均值。

對權(quán)重積分，獲得超參數(shù)邊緣概率分布：

p(t|α,σ2)=(2π)-N/2|B-1+ΦA(chǔ)-1ΦT|-1/2·

(20)

通過迭代計(jì)算得到超參數(shù)α和σ2的最優(yōu)解：

(21)

(22)

式中：Σii表示協(xié)方差矩陣Σ中第i個(gè)對角元素。

相關(guān)向量機(jī)的核函數(shù)主要分為局部核函數(shù)和全局核函數(shù)兩類：局部核函數(shù)學(xué)習(xí)能力強(qiáng)，泛化能力弱；全局核函數(shù)泛化能力強(qiáng)，學(xué)習(xí)能力弱。采取加權(quán)組合核函數(shù)，選取一個(gè)全局核函數(shù)多項(xiàng)式核函數(shù)和局部核函數(shù)高斯核函數(shù)進(jìn)行加權(quán)組合，來提高模型的預(yù)測精度和魯棒性，形式如下：

(23)

2.2 HPO算法原理

獵食者-獵物優(yōu)化算法是于2022年提出的一種新的優(yōu)化算法，具有收斂速度快、尋優(yōu)能力強(qiáng)的特點(diǎn)，并且人工設(shè)置參數(shù)少。

隨機(jī)生成種群的初始位置。

xi=rand(1,d)·(ub-lb)+lb

(24)

式中：d為目標(biāo)個(gè)數(shù)；ub、lb分別為搜索空間的上邊界和下邊界。

獵食者根據(jù)下式進(jìn)行搜索：

xi,j(t+1)=xi,j(t)+0.5{[2CZPpos(j)-xi,j(t)]+

[2(1-C)Zμ(j)-xi,j(t)]}

(25)

其中：x(t)為當(dāng)前獵食者的位置；x(t+1)為獵食者下一個(gè)位置；Ppos為獵物的位置；μ為所有位置的平均值；Z是由下式計(jì)算的自適應(yīng)參數(shù)。

P=R1

Z=R2?δIDX+R3?(～δIDX)

(26)

其中：R1和R3是[0，1]內(nèi)的隨機(jī)向量；δIDX是滿足條件P=0向量R1的索引值；C是探索和開發(fā)之間的平衡參數(shù)，計(jì)算如下：

(27)

其中：εit是當(dāng)前迭代值；εitmax是最大迭代次數(shù)。

計(jì)算所有位置的平均值：

(28)

計(jì)算歐幾里得距離：

(29)

距離位置平均值最大的搜索代理視為獵物Ppos：

Ppos=xi|iis index of Max(end)sort(Deuc)

(30)

為了解決則算法延遲收斂性的問題，考慮一種遞減機(jī)制，如下式所示：

kbest=round(C×N)

(31)

其中：N是搜索代理的數(shù)量。因此，獵物的位置計(jì)算公式可寫為

Ppos=xi|iis soredDeuc(kbest)

(32)

為了選擇獵人和獵物，結(jié)合式(25)和式(30)提出式(33)：

xi(t+1)=

(33)

其中：R5是[0，1]內(nèi)的隨機(jī)數(shù)；β是調(diào)節(jié)參數(shù)。

2.3 RVM參數(shù)優(yōu)化

采用HPO算法對混合核參數(shù)進(jìn)行尋優(yōu)。優(yōu)化的具體流程如圖1所示，優(yōu)化步驟如下：

圖1 參數(shù)優(yōu)化流程

步驟1，設(shè)置種群數(shù)、迭代次數(shù)以及優(yōu)化參數(shù)的上下界、目標(biāo)函數(shù)。

步驟2，計(jì)算適應(yīng)度值，記錄最優(yōu)位置。

步驟3，利用式(26)和(27)更新參數(shù)Z和C。

步驟4，根據(jù)式(33)更新獵食者或獵物的位置，并計(jì)算適應(yīng)度值和最優(yōu)位置。

步驟5，判斷是否滿足停止條件，如果滿足，則輸出最優(yōu)解，跳轉(zhuǎn)至步驟6，否則重復(fù)步驟3-4。

步驟6，利用最優(yōu)解構(gòu)建RVM預(yù)測模型。

3 壽命預(yù)測流程

剩余使用壽命定義為

RUL=Tend-Ti

(34)

(35)

式中：Tend為軸承失效時(shí)刻；Ti為當(dāng)前監(jiān)測點(diǎn)。為減少不同軸承間的差異，采用式(35)將軸承的剩余壽命進(jìn)行歸一化處理。

基于特征融合與HPO-RVM的軸承剩余壽命預(yù)測步驟如下：

步驟1，提取時(shí)域、頻域、時(shí)頻域特征。

步驟2，計(jì)算每個(gè)特征的魯棒性和相關(guān)性，采用Critic方法確定評價(jià)指標(biāo)的權(quán)值，然后計(jì)算每個(gè)特征的綜合得分。

步驟3，選取得分最高的前15個(gè)特征采用KECA降維融合，并以貢獻(xiàn)率大于90%的主成分作為退化特征。

步驟4，采用3σ方法確定軸承的初始故障時(shí)刻，當(dāng)故障發(fā)生以后啟動(dòng)剩余壽命預(yù)測機(jī)制。

步驟5，根據(jù)訓(xùn)練數(shù)據(jù)，采用HPO尋優(yōu)得到RVM核參數(shù)，代入RVM中得到HPO-RVM壽命預(yù)測模型。

步驟6，將測試集代入訓(xùn)練好的模型中進(jìn)行剩余壽命預(yù)測，如圖2所示。

圖2 壽命預(yù)測流程

步驟7，模型性能評價(jià)。

評價(jià)指標(biāo)：均方根誤差δRMSE、平均絕對誤差δMAE、決定系數(shù)R2。前2個(gè)評價(jià)指標(biāo)越小說明模型的預(yù)測能力越好，第3個(gè)評價(jià)指標(biāo)越大說明模型的預(yù)測能力越好[14]。

(36)

(37)

(38)

4 實(shí)驗(yàn)仿真驗(yàn)證

4.1 數(shù)據(jù)集介紹

文中實(shí)驗(yàn)驗(yàn)證采用的是XJTU-SY全壽命滾動(dòng)軸承數(shù)據(jù)集，圖3所示為搭建的實(shí)驗(yàn)平臺。測試所采用滾動(dòng)軸承的型號為LDK UER204。實(shí)驗(yàn)共測試了3種工況共計(jì)15個(gè)軸承，具體工況信息如表1所示。傳感器為PCB 352C33加速度傳感器，采樣頻率25.6 kHz，每間隔1 min，采樣1.28 s。詳細(xì)實(shí)驗(yàn)過程和信息參考文獻(xiàn)[15]。

以工況1中的5個(gè)軸承為研究對象進(jìn)行滾動(dòng)軸承的剩余壽命預(yù)測研究，分別命名為軸承1-1、軸承1-2、軸承1-3、軸承1-4、軸承1-5。

4.2 退化特征的構(gòu)建

以工況1為例，提取時(shí)域、頻域、時(shí)頻域39個(gè)特征，然后歸一化。軸承1-1的特征如圖4—6所示。然后根據(jù)評價(jià)指標(biāo)魯棒性和相關(guān)性計(jì)算每個(gè)特征的得分，采用Critic方法確定2個(gè)評價(jià)指標(biāo)的權(quán)重，然后得到綜合評價(jià)指標(biāo)。

圖4 軸承1-1時(shí)域特征

圖5 軸承1-1頻域特征

圖6 軸承1-1時(shí)頻域特征

根據(jù)式(4)計(jì)算工況1下5個(gè)軸承各個(gè)特征綜合得分的平均值。如圖7所示，選出前15個(gè)得分最高的特征，依次為37、36、29、22、2、10、1、3、17、34、6、33、32、26、35?？梢钥闯觯汉Y選出的特征包括時(shí)域特征、頻域特征和時(shí)頻域特征。

圖7 工況1軸承特征篩選

利用KECA方法將篩選出來的特征指標(biāo)進(jìn)行降維融合，從圖8中可以看出：工況1的各個(gè)軸承第一主成分的貢獻(xiàn)率都大于90%，因此選取第一主成分作為軸承的健康指標(biāo)。為降低不同軸承間的差異性，對健康指標(biāo)進(jìn)行歸一化處理。為了降低噪聲或異常值引起的指標(biāo)波動(dòng)，采用局部加權(quán)回歸方法進(jìn)行平滑處理，得到工況1的健康指標(biāo)如圖9所示。

圖8 工況1各軸承主成分貢獻(xiàn)率

圖9 工況1軸承健康指標(biāo)

軸承運(yùn)行初期，退化特征不明顯，特征指標(biāo)平穩(wěn)，對它進(jìn)行壽命預(yù)測工程意義不大。當(dāng)發(fā)生早期故障以后，軸承開始退化，啟動(dòng)剩余壽命預(yù)測。在進(jìn)行壽命預(yù)測時(shí)只是用退化階段的數(shù)據(jù)，避免了非退化數(shù)據(jù)對壽命預(yù)測的精度造成干擾。文中采用3σ方法檢測軸承的早期故障時(shí)刻，從而確定壽命預(yù)測的初始時(shí)間，具體方法參考文獻(xiàn)[16]。得到軸承1-1到軸承1-5發(fā)生故障的時(shí)刻為70、36、59、121、34 min。因此將發(fā)生早期故障以后的退化指標(biāo)歸一化作為剩余壽命預(yù)測的數(shù)據(jù)，從而提高預(yù)測的精度。

4.3 滾動(dòng)軸承壽命預(yù)測

壽命預(yù)測階段，以前4個(gè)軸承作為訓(xùn)練集訓(xùn)練壽命預(yù)測模型，最后1個(gè)軸承作為測試集，驗(yàn)證模型的預(yù)測性能。采用HPO算法尋優(yōu)得到RVM的組合核參數(shù)，HPO算法的種群數(shù)設(shè)置為100，最大迭代次數(shù)為100，得到的最優(yōu)參數(shù)用于RVM預(yù)測模型的訓(xùn)練。將待測數(shù)據(jù)輸入到模型中進(jìn)行預(yù)測，得到的壽命預(yù)測結(jié)果如圖10所示?？梢钥闯觯侯A(yù)測的總體趨勢與實(shí)際壽命百分比接近，大部分實(shí)際值都在預(yù)測結(jié)果95%的置信區(qū)間內(nèi)。

圖10 所提出的預(yù)測模型預(yù)測結(jié)果

將文中提出的預(yù)測模型與BP神經(jīng)網(wǎng)絡(luò)、極限學(xué)習(xí)機(jī)ELM、支持向量機(jī)SVM、最小二乘支持向量機(jī)LSSVM、相關(guān)向量機(jī)RVM、優(yōu)化支持向量機(jī)、優(yōu)化最小二乘支持向量機(jī)、優(yōu)化高斯核函數(shù)的相關(guān)向量機(jī)等預(yù)測模型進(jìn)行對比。預(yù)測結(jié)果如圖11所示。

由圖11可以看出：文中所提的模型預(yù)測結(jié)果更加接近真實(shí)值，文中模型具有一定的優(yōu)勢。然后計(jì)算各個(gè)模型預(yù)測性能評價(jià)指標(biāo)。δRMSE和δMAE值越小越好，R2越大越好。由圖12可以看出：文中所提的模型δRMSE和R2均為最優(yōu)結(jié)果，δMAE指標(biāo)僅次于基于高斯核函數(shù)的相關(guān)向量機(jī)模型，但結(jié)果相近，總體上文中所提的模型更具優(yōu)勢。

圖12 不同預(yù)測模型預(yù)測結(jié)果評價(jià)

將文中采用的HPO優(yōu)化算法與常用的經(jīng)典優(yōu)化算法粒子群優(yōu)化算法(PSO)、遺傳優(yōu)化算法(GA)，以及較新的優(yōu)化算法哈里斯鷹優(yōu)化算法(HHO)、晶體結(jié)構(gòu)優(yōu)化算法(CRYSTAL)等進(jìn)行對比試驗(yàn)，得到預(yù)測結(jié)果的評價(jià)指標(biāo)如圖13所示?？梢钥闯觯篐PO優(yōu)化算法可以取得較好的預(yù)測結(jié)果。

圖13 不同算法優(yōu)化模型預(yù)測結(jié)果評價(jià)

文中是以最小均方根誤差作為HPO優(yōu)化RVM核參數(shù)的目標(biāo)函數(shù)，實(shí)際應(yīng)用中可以根據(jù)需要設(shè)置不同的目標(biāo)函數(shù)，如δRMSE、δMAE、1-R2、δRMSE+δMAE+(1-R2)等最小值作為目標(biāo)函數(shù)。當(dāng)采用以上目標(biāo)函數(shù)訓(xùn)練模型，測試模型的預(yù)測性能如圖14所示?？梢钥闯觯涸谟?xùn)練模型時(shí)采用不同的目標(biāo)函數(shù)可以使相應(yīng)的指標(biāo)達(dá)到最小值，以滿足相應(yīng)的需求。

圖14 不同目標(biāo)函數(shù)優(yōu)化模型預(yù)測結(jié)果評價(jià)

5 結(jié)論

文中提出基于特征融合與HPO-RVM的滾動(dòng)軸承剩余壽命預(yù)測方法，利用軸承全壽命數(shù)據(jù)集進(jìn)行方法驗(yàn)證，得出以下結(jié)論：

(1)采用綜合評價(jià)指標(biāo)篩選出來的特征包含時(shí)域、頻域、時(shí)頻域特征，將篩選出來的特征進(jìn)一步用KECA降維融合，選取能夠反映90%以上特征信息的主成分作為退化特征，不僅降低了計(jì)算量，而且特征信息更加豐富，避免了單一特征不能較好地反映軸承的退化特征。

(2)采用檢測到軸承發(fā)生故障以后再進(jìn)行壽命預(yù)測的策略，避免了非退化特征對剩余壽命預(yù)測精度的影響。

(3)HPO優(yōu)化混合核RVM的預(yù)測模型，其預(yù)測精度優(yōu)于BP、ELM、SVM等模型；所采用的優(yōu)化算法在此實(shí)驗(yàn)中優(yōu)于PSO、GA等算法；并討論了不同的目標(biāo)函數(shù)對預(yù)測模型的影響。綜合說明文中所提的壽命預(yù)測模型具有一定的優(yōu)勢。