張 鋼，譚 波，梁偉閣，田福慶

(海軍工程大學(xué)兵器工程學(xué)院,湖北 武漢 430033)

0 引言

滾動(dòng)軸承是各類旋轉(zhuǎn)機(jī)械中最常用的通用零部件之一，也是旋轉(zhuǎn)機(jī)械易損件之一。滾動(dòng)軸承一旦發(fā)生故障，將導(dǎo)致災(zāi)難性事故。因此，預(yù)測(cè)滾動(dòng)軸承剩余壽命具有重要意義[1]。

智能預(yù)測(cè)技術(shù)是滾動(dòng)軸承剩余壽命預(yù)測(cè)的主要手段之一，其中相關(guān)向量機(jī)(relevance vector machine, RVM)具有較好的稀疏性、泛化能力，逐漸成為剩余壽命預(yù)測(cè)領(lǐng)域的研究熱點(diǎn)。文獻(xiàn)[2]針對(duì)風(fēng)電功率序列非線性、非平穩(wěn)性等特點(diǎn)，提出一種基于混沌布谷鳥搜索算法優(yōu)化相關(guān)向量機(jī)的短期風(fēng)電功率預(yù)測(cè)方法。文獻(xiàn)[3]提出一種基于相關(guān)向量機(jī)的電力負(fù)荷中期預(yù)測(cè)方法。RVM雖然具有較高的中短期預(yù)測(cè)精度，但是由于機(jī)械部件性能退化復(fù)雜性，導(dǎo)致RVM預(yù)測(cè)算法出現(xiàn)長(zhǎng)期趨勢(shì)預(yù)測(cè)精度不高的問題[4]。同時(shí)，相關(guān)向量機(jī)模型在預(yù)測(cè)過程中包括數(shù)據(jù)重構(gòu)和迭代預(yù)測(cè)兩個(gè)過程，計(jì)算效率較低[5-7]。本文針對(duì)此問題，提出了基于改進(jìn)相關(guān)向量機(jī)的滾動(dòng)軸承剩余壽命預(yù)測(cè)方法。

1 改進(jìn)相關(guān)向量機(jī)原理

1.1 相關(guān)向量機(jī)原理

ti=y(xi)+ε

(1)

式(1)中，N是樣本數(shù)，y(xi)表示非線性函數(shù)，ε是獨(dú)立同分布的高斯噪聲，且ε～N(0,σ2)。

非線性函數(shù)y(x)可以表示為核函數(shù)乘以相應(yīng)權(quán)重系數(shù)后的求和：

(2)

式(2)中，w=[w1,w2,…,wN]T是權(quán)重向量，wi表示第i個(gè)訓(xùn)練數(shù)據(jù)對(duì)應(yīng)的權(quán)值，w0表示偏置參數(shù)，K(x,xi)是核函數(shù)。因此，目標(biāo)輸出值ti∈R的概率分布是均值為y(x,w)，方差為σ2的高斯分布。在整個(gè)訓(xùn)練數(shù)據(jù)集上的似然函數(shù)可以表示為：

(3)

式(3)中，t=[t1,t2,…,tN]T；ω為全部權(quán)重系數(shù)和偏置組成的向量，ω=[w0,w1,…,wN]T；Φ=[φ1,φ2,…,φN]T表示輸入向量帶入核函數(shù)后得到的N×(N+1)維矩陣，φi=[1,K(xi,x1),…,K(xi,xN)]，i=1,2,…,N。

假設(shè)ωi概率分布為0附近的高斯分布，即滿足：

(4)

式(4)中，αi為權(quán)重ωi對(duì)應(yīng)的超參數(shù)，決定權(quán)重的分布范圍。各權(quán)重相互獨(dú)立，因此，對(duì)所有權(quán)重ω有：

(5)

式(5)中，α為由N+1個(gè)超參數(shù)αi組成的向量，α=[α0,α1,…,αN]T。

根據(jù)參數(shù)的先驗(yàn)分布，結(jié)合貝葉斯推理可知，所有未知參數(shù)的后驗(yàn)分布可以表示為：

(6)

因此，根據(jù)新的輸入樣本x*，預(yù)測(cè)輸出值t的概率分布為：

(7)

式(7)的近似解為：

(8)

(9)

式(9)中，C=σ2I+ΦA(chǔ)-1ΦT，μ和Σ是權(quán)重后驗(yàn)分布的均值和方差：

μ=σ2ΣΦTt

(10)

Σ=(A+σ-2ΦTΦ)-1

(11)

式(10)中，A=diag(α0,α1,…,αN)。

(12)

(13)

式(12)、式(13)中，γi=1-αiΣii，γi∈[0,1]，Σii表示權(quán)重方差矩陣的對(duì)角線元素。

當(dāng)給定新的輸入樣本x*，輸出值t*的概率分布為：

(14)

由式(8)和式(9)可知，式(14)積分項(xiàng)的兩部分都是高斯函數(shù)乘積形式，因此積分后的結(jié)果為：

(15)

其中，

y*=μTΦ(x*)

(16)

(17)

由式(16)、式(17)可知，相關(guān)向量機(jī)不僅能夠給出預(yù)測(cè)值的期望，同時(shí)能夠給定預(yù)測(cè)值的置信區(qū)間。

1.2 改進(jìn)相關(guān)向量機(jī)模型

相關(guān)向量機(jī)雖然能夠提供預(yù)測(cè)結(jié)果的置信區(qū)間，但是面對(duì)大量訓(xùn)練數(shù)據(jù)，存在計(jì)算效率低、長(zhǎng)期預(yù)測(cè)能力較差、預(yù)測(cè)結(jié)果受核函數(shù)的影響較大等問題。本節(jié)提出一種結(jié)合多項(xiàng)式函數(shù)的改進(jìn)相關(guān)向量機(jī)方法，并將其應(yīng)用于滾動(dòng)軸承剩余壽命預(yù)測(cè)中。

多項(xiàng)式回歸模型具有較強(qiáng)的曲線擬合能力，長(zhǎng)期預(yù)測(cè)精度較高，但是屬于點(diǎn)估計(jì)方法，無(wú)法根據(jù)歷史數(shù)據(jù)獲得預(yù)測(cè)結(jié)果的置信區(qū)間。本節(jié)將多項(xiàng)式回歸模型較強(qiáng)的長(zhǎng)期預(yù)測(cè)能力與相關(guān)向量機(jī)相結(jié)合，提出一種改進(jìn)相關(guān)向量機(jī)的剩余壽命預(yù)測(cè)模型。

多項(xiàng)式回歸模型的數(shù)學(xué)表達(dá)式為：

y=a0+a1x+a2x2+…+anxn

(18)

式(18)中，y是目標(biāo)函數(shù)值，x是輸入值，a0,a1,a2,…,an是模型系數(shù)。

2 基于改進(jìn)相關(guān)向量機(jī)的剩余壽命預(yù)測(cè)方法

基于改進(jìn)相關(guān)向量機(jī)的剩余壽命預(yù)測(cè)算法流程如下：

輸出：預(yù)測(cè)值y*和預(yù)測(cè)值的置信區(qū)間；

步驟2 初始化參數(shù){αi}和σ2；

步驟4 根據(jù)式(12)和式(13)更新αi和σ2；

步驟6 利用相關(guān)向量xRV確定多項(xiàng)式回歸模型系數(shù)，計(jì)算各訓(xùn)練樣本輸入值xi到多項(xiàng)式曲線的距離；

步驟8 給定預(yù)測(cè)數(shù)據(jù)x*，代入式(18)，計(jì)算得到預(yù)測(cè)值y*，利用式(17)計(jì)算得到預(yù)測(cè)值的置信區(qū)間。

3 試驗(yàn)與分析

3.1 試驗(yàn)介紹

XJTU-SY滾動(dòng)軸承試驗(yàn)數(shù)據(jù)是由西安交通大學(xué)設(shè)計(jì)科學(xué)與基礎(chǔ)構(gòu)件研究所和長(zhǎng)興昇陽(yáng)科技有限公司合作試驗(yàn)采集得到[8]。

XJTU-SY滾動(dòng)軸承試驗(yàn)平臺(tái)由一臺(tái)交流感應(yīng)電機(jī)、一臺(tái)電機(jī)轉(zhuǎn)速控制器、一個(gè)支撐軸、兩個(gè)支撐軸承(重型滾動(dòng)軸承)和一個(gè)液壓加載系統(tǒng)組成，如圖1所示。該試驗(yàn)平臺(tái)旨在開展?jié)L動(dòng)軸承在不同工況下(如不同徑向力和轉(zhuǎn)速)的加速退化試驗(yàn)。液壓加載系統(tǒng)產(chǎn)生徑向力并施加于被測(cè)軸承殼體上，交流感應(yīng)電機(jī)的轉(zhuǎn)速控制器可以控制滾動(dòng)軸承轉(zhuǎn)速。為了采集被測(cè)軸承的振動(dòng)信號(hào)，在被測(cè)軸承殼體上放置兩個(gè)PCB 352C33型加速度計(jì)，一個(gè)安裝在水平軸承，另一個(gè)安裝在垂直軸上。加速度計(jì)的采樣頻率設(shè)置為25.6 kHz，采樣時(shí)間為1.28 s，采樣周期為1 min。

圖1 XJTU-SY滾動(dòng)軸承試驗(yàn)平臺(tái)Fig.1 XJTU-SY bearing experimental platform

試驗(yàn)所用軸承型號(hào)為L(zhǎng)DK UER204軸承，共有三種運(yùn)行工況，每類工況下做五組試驗(yàn)。各工況條件如表1所示。

表1 試驗(yàn)工況

由于負(fù)荷加載在水平方向，水平方向的傳感器能采集到更多的運(yùn)行狀態(tài)信息。因此，采用水平方向的振動(dòng)加速度信號(hào)開展剩余壽命預(yù)測(cè)研究。水平方向的原始振動(dòng)加速度信號(hào)如圖2所示。

圖2 水平方向的原始振動(dòng)加速度信號(hào)Fig.2 The original vibration acceleration signal of horizontal direction

由圖2可知，在正常運(yùn)行工況下，振動(dòng)加速度信號(hào)僅有小范圍的波動(dòng)，因此很難利用這類歷史數(shù)據(jù)預(yù)測(cè)剩余壽命。當(dāng)滾動(dòng)軸承進(jìn)入性能退化階段時(shí)，振動(dòng)加速度值隨著運(yùn)行時(shí)間的增加而不斷增加，其包含豐富的性能退化信息，因此選用這一階段的振動(dòng)加速度信號(hào)作為歷史監(jiān)測(cè)數(shù)據(jù)預(yù)測(cè)剩余壽命。預(yù)測(cè)起始點(diǎn)的選擇采用文獻(xiàn)[9]介紹的自適應(yīng)性能退化檢測(cè)方法。當(dāng)滾動(dòng)軸承的振動(dòng)加速度幅值超過20g時(shí)，將存在嚴(yán)重的安全問題，本文將振動(dòng)加速度為20g時(shí)的滾動(dòng)軸承定義為失效軸承。

3.2 剩余壽命預(yù)測(cè)

確定滾動(dòng)軸承起始預(yù)測(cè)時(shí)間點(diǎn)Ts后，利用訓(xùn)練好的改進(jìn)向量機(jī)模型性能退化曲線，當(dāng)性能退化曲線達(dá)到失效閾值(20g)時(shí)，對(duì)應(yīng)的時(shí)間即為失效時(shí)刻Tfail_pred。當(dāng)前時(shí)刻到預(yù)測(cè)失效時(shí)刻的差值即為滾動(dòng)軸承剩余壽命，根據(jù)文獻(xiàn)[10]，滾動(dòng)軸承剩余壽命定義為：

RUL(Ts)=inf{t:f(t+Ts)≥η|f}

(19)

式(9)中，RUL(Ts)是自預(yù)測(cè)起始時(shí)刻Ts的剩余壽命，f(t+Ts)是t+Ts時(shí)刻預(yù)測(cè)的性能退化狀態(tài)，f所有的健康因子，η是失效閾值。

本文提取最大幅值(maximum amplitude, MA)作為滾動(dòng)軸承健康因子，得到不同工況下比較典型的全壽命周期健康因子曲線如圖3所示。

圖3 不同工況下典型的全壽命周期健康因子曲線Fig.3 Typical life cycle health factor curves under different working conditions

以工況3下bearing3_1為例，利用本文提出的基于改進(jìn)相關(guān)向量機(jī)的剩余壽命預(yù)測(cè)方法得到性能退化曲線如圖4所示。

圖4 Bearing3_1的性能退化預(yù)測(cè)曲線Fig.4 The predicted degradation curve of Bearing3_1

圖4中，折線表示健康因子曲線，圓點(diǎn)表示相關(guān)向量，通過折線和圓點(diǎn)的光滑曲線表示多項(xiàng)式回歸模型擬合曲線，最左側(cè)的豎直點(diǎn)劃線后的光滑曲線表示性能退化預(yù)測(cè)曲線，兩條點(diǎn)線之間表示預(yù)測(cè)值的置信區(qū)間。豎直方向，由左至右，依次表示預(yù)測(cè)起始時(shí)刻，預(yù)測(cè)失效時(shí)刻，預(yù)測(cè)失效時(shí)刻的上置信區(qū)間點(diǎn)，預(yù)測(cè)失效時(shí)刻下置信區(qū)間點(diǎn)。

預(yù)測(cè)起始時(shí)刻為Ts=2 430 min，預(yù)測(cè)軸承失效時(shí)刻為Tfail_pred=2 535 min，實(shí)際的失效時(shí)刻為Tfail_real=2 537 min。因此，帶入剩余壽命計(jì)算公式(19)，預(yù)測(cè)的剩余壽命為RULpred=105 min，實(shí)際剩余壽命為RULreal=107 min，剩余壽命預(yù)測(cè)結(jié)果的置信區(qū)間為[103,108]min。由圖4和以上數(shù)據(jù)可以看出，本文提出的方法能夠有效預(yù)測(cè)滾動(dòng)軸承剩余壽命，實(shí)際剩余壽命包含在置信區(qū)間中，預(yù)測(cè)的剩余壽命小于實(shí)際剩余壽命，在工業(yè)應(yīng)用中，有利于提前預(yù)警。

3.3 對(duì)比分析

為了進(jìn)一步說明本文所提方法的有效性，分別利用原始相關(guān)向量機(jī)、深度置信網(wǎng)絡(luò)和粒子濾波方法預(yù)測(cè)滾動(dòng)軸承剩余壽命。以bearing3_1為例，起始預(yù)測(cè)點(diǎn)為Ts=2 430 min，各類方法得到的剩余壽命預(yù)測(cè)結(jié)果如圖5所示。

圖5 不同方法得到的剩余壽命預(yù)測(cè)結(jié)果Fig.5 The remaining useful life prediction results of different methods

圖5中，曲線分別表示實(shí)際剩余壽命，基于改進(jìn)RVM預(yù)測(cè)的剩余壽命，基于RVM預(yù)測(cè)的剩余壽命[11]，基于深度置信網(wǎng)絡(luò)(deep belief network, DBN)預(yù)測(cè)的剩余壽命[12]，基于粒子濾波(particle filtering, PF)算法得到的剩余壽命[13]。其中深度置信網(wǎng)絡(luò)是一類目前較火的基于深度學(xué)習(xí)網(wǎng)絡(luò)的剩余壽命預(yù)測(cè)模型，屬于數(shù)據(jù)驅(qū)動(dòng)的剩余壽命預(yù)測(cè)方法；粒子濾波屬于典型的機(jī)理驅(qū)動(dòng)的剩余壽命預(yù)測(cè)方法。

由圖5可知，基于改進(jìn)相關(guān)向量機(jī)的剩余壽命預(yù)測(cè)方法能夠較快的收斂，且預(yù)測(cè)精度較其他三種預(yù)測(cè)方法有提高。相較于RVM，長(zhǎng)期預(yù)測(cè)精度得到有效提高。相對(duì)于RVM和PF，基于DBN的剩余壽命預(yù)測(cè)方法預(yù)測(cè)精度提高，但是其預(yù)測(cè)結(jié)果屬于點(diǎn)估計(jì)，無(wú)法提供預(yù)測(cè)置信區(qū)間。改進(jìn)RVM在預(yù)測(cè)點(diǎn)Ts=2 460時(shí)，預(yù)測(cè)結(jié)果偏差較大，主要原因是當(dāng)起始預(yù)測(cè)點(diǎn)為Ts=2 460時(shí)，由圖4可知滾動(dòng)軸承健康因子曲線突然變陡峭，從而預(yù)測(cè)得到的剩余壽命出現(xiàn)較大誤差，但隨著數(shù)據(jù)量的增加，采用改進(jìn)相關(guān)向量機(jī)的剩余壽命預(yù)測(cè)結(jié)果迅速收斂于實(shí)際剩余壽命，說明該方法具有較好的魯棒性。

在運(yùn)算量方面，本文所提的改進(jìn)方法主要包括兩部分的計(jì)算內(nèi)容：一是訓(xùn)練階段；二是預(yù)測(cè)階段。在訓(xùn)練階段，利用相關(guān)向量機(jī)模型對(duì)訓(xùn)練數(shù)據(jù)稀疏化處理，得到相關(guān)向量，利用相關(guān)向量確定多項(xiàng)式回歸模型參數(shù)，該階段計(jì)算時(shí)間要長(zhǎng)于單獨(dú)訓(xùn)練相關(guān)向量機(jī)模型的時(shí)間。預(yù)測(cè)階段，利用多項(xiàng)式回歸模型直接預(yù)測(cè)性能退化趨勢(shì)，省略數(shù)據(jù)重構(gòu)及迭代預(yù)測(cè)過程，因此該階段的計(jì)算時(shí)間比利用相關(guān)向量機(jī)模型進(jìn)行預(yù)測(cè)的時(shí)間短。以文中bearing3_1訓(xùn)練、測(cè)試過程為例，總得計(jì)算時(shí)間為21.4 s，而單獨(dú)使用相關(guān)向量機(jī)模型的計(jì)算時(shí)間為23 s。綜上，本文所提方法的運(yùn)算量比相關(guān)向量機(jī)模型運(yùn)算量要小，提高了計(jì)算效率。

4 結(jié)論

本文提出了基于改進(jìn)相關(guān)向量機(jī)的滾動(dòng)軸承剩余壽命預(yù)測(cè)方法。該方法綜合多項(xiàng)式回歸模型長(zhǎng)期趨勢(shì)預(yù)測(cè)精度高，相關(guān)向量機(jī)能夠提供稀疏的相關(guān)向量和提供預(yù)測(cè)結(jié)果置信區(qū)間的優(yōu)勢(shì)，有效提高了滾動(dòng)軸承壽命預(yù)測(cè)精度，且所提方法在預(yù)測(cè)階段不需要數(shù)據(jù)重構(gòu)和迭代預(yù)測(cè)過程，提高了計(jì)算效率。試驗(yàn)結(jié)果表明，本文所提方法能夠有效提高剩余壽命精度和計(jì)算效率，并提供預(yù)測(cè)置信區(qū)間，有利于制定維修計(jì)劃。