張?zhí)镉?王慶鋒* 舒 悅 肖 旺

(1.北京化工大學(xué)機(jī)電工程學(xué)院,北京 100029;2.北京化工大學(xué)高端機(jī)械裝備健康監(jiān)控及自愈化北京市重點(diǎn)實(shí)驗(yàn)室,北京 100029;3.合肥通用機(jī)械研究院有限公司, 合肥 230031;4.國家管網(wǎng)集團(tuán)聯(lián)合管道有限責(zé)任公司西部分公司,烏魯木齊 830013)

引 言

滾動軸承廣泛應(yīng)用于旋轉(zhuǎn)類機(jī)械設(shè)備中,是關(guān)鍵易損件之一。 對智能傳感器所采集的滾動軸承實(shí)時振動數(shù)據(jù)的分析,可為預(yù)防性維修決策提供依據(jù),延長滾動軸承的使用壽命,大大提高生產(chǎn)效率,因此關(guān)于滾動軸承剩余使用壽命(RUL)預(yù)測的研究近年來越來越受到人們的重視[1-2]。

目前關(guān)于剩余使用壽命預(yù)測的方法可歸納為兩種[3]:第一種是基于數(shù)據(jù)驅(qū)動的方法,第二種是基于模型驅(qū)動的方法。 基于數(shù)據(jù)驅(qū)動的方法是在獲取一定數(shù)據(jù)量的基礎(chǔ)上,運(yùn)用統(tǒng)計分析或機(jī)器學(xué)習(xí)方法對剩余使用壽命進(jìn)行預(yù)測。 該種方法對數(shù)據(jù)量的大小有要求,所以往往很難符合實(shí)際工程化應(yīng)用需求,存在一定的局限性。 基于模型驅(qū)動的方法則是為描述系統(tǒng)的衰退趨勢建立數(shù)學(xué)模型,并根據(jù)觀測數(shù)據(jù)對數(shù)學(xué)模型的參數(shù)進(jìn)行實(shí)時更新,進(jìn)而對剩余使用壽命進(jìn)行預(yù)測。 相較于基于數(shù)據(jù)驅(qū)動的方法,基于模型驅(qū)動的方法可以充分利用已知數(shù)據(jù)包含的故障信息,并得出相對穩(wěn)定的剩余使用壽命預(yù)測結(jié)果。

常見的基于模型驅(qū)動的方法包括隱馬爾可夫模型、卡爾曼濾波、粒子濾波等。 其中粒子濾波預(yù)測方法基于貝葉斯濾波算法提出,可改善其他方法在處理長期預(yù)測時的不確定性問題,現(xiàn)已廣泛應(yīng)用于剩余使用壽命預(yù)測領(lǐng)域。 Li 等[4]提出了基于粒子濾波的預(yù)測方法,使用相關(guān)矩陣聚類和加權(quán)算法計算健康指標(biāo),并運(yùn)用粒子濾波算法預(yù)測健康指標(biāo)的變化趨勢,最后得到軸承的剩余使用壽命。 祝志遠(yuǎn)等[5]提出將粒子濾波算法應(yīng)用于疲勞裂紋擴(kuò)展的參數(shù)估計和剩余壽命預(yù)測上,預(yù)測得到了剩余壽命中值和置信區(qū)間。 Lei 等[6]提出一種基于模型的滾動軸承剩余使用壽命預(yù)測方法,該方法包括健康指標(biāo)構(gòu)建和RUL 預(yù)測兩個模塊,在第一個部分構(gòu)建了一個名為加權(quán)最小量化誤差的新健康指標(biāo),在第二個部分使用粒子濾波預(yù)測RUL,減小了預(yù)測中的累計誤差。 這些采用粒子濾波算法預(yù)測剩余壽命的研究可獲得更加準(zhǔn)確的預(yù)測結(jié)果,充分體現(xiàn)了粒子濾波算法在解決長期不確定性預(yù)測方面的優(yōu)勢,但其本身仍存在多次迭代后粒子權(quán)重退化問題。

為解決上述粒子濾波算法存在的缺陷,余臻等[7]提出使用無跡粒子濾波方法預(yù)測航天發(fā)動機(jī)的排氣溫度,獲得了比粒子濾波算法更優(yōu)良的效果。許雨晨等[8]提出粒子群優(yōu)化算法改進(jìn)的粒子濾波預(yù)測算法,可以較為準(zhǔn)確地預(yù)測出滾動軸承的剩余使用壽命。 賀寧等[9]提出使用天牛須搜索算法優(yōu)化粒子濾波的重采樣過程,解決了粒子多樣性喪失問題,并將該方法應(yīng)用于電池的剩余使用壽命預(yù)測,獲得了接近真實(shí)電池壽命的結(jié)果,預(yù)測中的累計誤差明顯減小。 徐仁義等[10]提出一種基于均方諧噪比的指標(biāo),之后運(yùn)用改進(jìn)正則化的粒子濾波算法優(yōu)化重采樣過程,該算法利用基于歐式距離的核函數(shù)改進(jìn)了粒子濾波算法中的重采樣過程,在預(yù)測滾動軸承的剩余使用壽命中獲得了更接近真實(shí)剩余使用壽命的結(jié)果。

此外,在預(yù)測過程中,還會存在特征指標(biāo)波動性大、無法確定合適的置信區(qū)間進(jìn)而導(dǎo)致預(yù)測精度下降的問題,現(xiàn)在通常采用設(shè)置95%或99%置信區(qū)間的手段,但該方法主觀選擇的程度較大,無法適用于全部的特征指標(biāo),泛化性較差。 為解決上述問題,本文提出一種融合Hodrick-Prescott(HP)趨勢濾波-邊界線(HPTF-BL)、獵人獵物優(yōu)化算法(HPO)改進(jìn)粒子濾波(PF)的滾動軸承剩余使用壽命預(yù)測方法。以振動信號為例,首先提取表達(dá)設(shè)備退化趨勢的特征,然后通過HPTF-BL 對特征進(jìn)行處理,得到特征的上下退化邊界線與主要退化趨勢,解決了較難選擇合適的置信區(qū)間且選擇主觀性大的問題,并減少了預(yù)測誤差;最后使用HPO-PF 算法對滾動軸承的剩余使用壽命進(jìn)行預(yù)測,解決了粒子濾波算法中粒子權(quán)重退化問題,同時提高了預(yù)測精度。

1 滾動軸承剩余使用壽命預(yù)測模型

1.1 基于小波信息熵的特征提取方法

何正嘉等[11]提出用歸一化的特征指標(biāo)來評判機(jī)械設(shè)備的健康狀態(tài),例如相關(guān)系數(shù)法、凝聚函數(shù)法、小波信息熵法等,這些方法計算出的特征指標(biāo)均在0 ～1 之間。 由于小波信息熵法使用了熵值這一復(fù)雜性衡量指標(biāo),該指標(biāo)可反映出機(jī)械設(shè)備更多的故障信息,因此對于一臺包含滾動軸承的機(jī)械設(shè)備,適合使用小波信息熵法計算特征指標(biāo),來評判其運(yùn)行狀態(tài)。

設(shè)備運(yùn)行中產(chǎn)生振動信號,通過小波包分解l次并對每個頻帶重構(gòu)得到2l個分解信號xli(k),第i個分解頻帶信號的能量Eli和相對能量分別為

式中,i=1,2,…,2l;k=1,2,…,n,n∈Z。 相對能量總和

小波信息熵Ent定義為

式中,對數(shù)的底取2l,則Ent∈[0,1]。 設(shè)備運(yùn)行狀態(tài)數(shù)增加,不確定性增強(qiáng),小波信息熵值必然增大,則該時刻旋轉(zhuǎn)機(jī)械的運(yùn)行狀態(tài)越差;小波信息熵值越小,說明該時刻旋轉(zhuǎn)機(jī)械的運(yùn)行狀態(tài)越好[12]。

由于滾動軸承的性能退化具有單調(diào)不可逆性,可將單調(diào)性作為特征指標(biāo)與滾動軸承性能退化一致性的評判標(biāo)準(zhǔn)[13]。 當(dāng)特征指標(biāo)隨時間單調(diào)遞增或者遞減時,單調(diào)性為1;當(dāng)特征指標(biāo)隨機(jī)波動時,單調(diào)性為0。

單調(diào)性計算公式如下。

式中,X是時間序列,ε(x)為單位階躍函數(shù),x≥0時,ε(x)為1;x＜0 時,n為總樣本數(shù)。

采用在美國辛辛那提大學(xué)智能維護(hù)系統(tǒng)(IMS)中心實(shí)驗(yàn)平臺上采集的數(shù)據(jù)集2(IMS2)的1 號軸承數(shù)據(jù)[14]作為驗(yàn)證數(shù)據(jù),計算當(dāng)下常用的4 種用于滾動軸承剩余使用壽命預(yù)測的特征指標(biāo)的單調(diào)性,4種指標(biāo)分別為有效值(標(biāo)簽為1)、整流平均值(標(biāo)簽為2)、峭度值(標(biāo)簽為3)和譜距離指標(biāo)(標(biāo)簽為4),并計算基于小波信息熵(標(biāo)簽為5)的特征指標(biāo)的單調(diào)性,對比結(jié)果如圖1 所示。

圖1 幾種特征指標(biāo)的單調(diào)性對比Fig.1 Monotonicity comparison of several indexes

由圖1 的對比可知,小波信息熵指標(biāo)比其他4種特征指標(biāo)具有更好的單調(diào)性,適用于作為后續(xù)滾動軸承剩余使用壽命預(yù)測的指標(biāo)。

1.2 基于HPTF-BL 的特征處理方法

1.2.1 HP 濾波

HP 趨勢濾波由Hodrick 等[15]提出,是經(jīng)濟(jì)學(xué)中常用的數(shù)據(jù)分析方法,可將數(shù)據(jù)分解為長期趨勢項(xiàng)和短期波動項(xiàng)。 由于HP 趨勢濾波具有降低數(shù)據(jù)噪聲影響的作用,現(xiàn)常被用于提取各種時間序列的趨勢,在數(shù)據(jù)預(yù)測及產(chǎn)品性能退化的可靠性分析等領(lǐng)域應(yīng)用廣泛。

HP 趨勢濾波是一種高通濾波器,利用了最小二乘損失函數(shù),采用l2 范數(shù)對二次差分矩陣進(jìn)行計算,可將時間序列xt分解為長期趨勢項(xiàng)xk和具有隨機(jī)波動特性的波動項(xiàng)xc。 其中,趨勢項(xiàng)xk被定義為式(5)的最小化函數(shù)解。

式中,i為時間序列數(shù)據(jù)的序號,n為時間序列的樣本個數(shù),右側(cè)第一項(xiàng)表示趨勢項(xiàng)對原序列的跟蹤程度,第二項(xiàng)表示趨勢項(xiàng)的光滑程度,λ為平滑參數(shù),用于控制趨勢項(xiàng)的平滑程度。 對式(5)求xki序列的一階偏導(dǎo),即可獲得趨勢序列xk。

式中G為系數(shù)矩陣,I為單位矩陣。 當(dāng)λ→0 時,趨勢項(xiàng)對序列的跟蹤程度達(dá)到最大;當(dāng)λ→+∞時,趨勢項(xiàng)序列光滑程度達(dá)到最大;當(dāng)λ=0 時,HP 濾波方法即退化為最小二乘法。 通過上述HP 濾波方法即可將時間序列分解為周期項(xiàng)和波動項(xiàng)。

1.2.2 HPTF-BL 處理方法

本文提出一種針對波動性大的健康指標(biāo)的處理方法—HP 趨勢濾波-邊界線(HPTF-BL)方法,可以很好地解決健康指標(biāo)波動性大、置信區(qū)間難以確定的問題,減少最后預(yù)測結(jié)果的誤差。

HPTF-BL 方法是在HP 濾波的基礎(chǔ)上對其公式進(jìn)行形式改寫的一種方法,其主要步驟如下:

1)窗口劃分,即對獲得的已知時間序列xt進(jìn)行分割,獲得每組包含m個數(shù)據(jù)的多個窗口;

2)對每個窗口內(nèi)的數(shù)據(jù)大小進(jìn)行排序,獲得該窗口內(nèi)的最大值與最小值;

3)用最大值代替窗口內(nèi)所有的數(shù)據(jù),得到上邊界xup,用最小值代替窗口內(nèi)所有數(shù)據(jù),得到下邊界xlow;

4)用HP 濾波對時間序列xt、步驟3)中得到的上邊界xup和下邊界xlow進(jìn)行處理,得到主要趨勢項(xiàng)xk,xk的計算公式如式(5)所示。

處理后的健康指標(biāo)HI如式(7)所示。

1.3 改進(jìn)粒子濾波算法

1.3.1 粒子濾波算法

粒子濾波是在傳統(tǒng)非線性濾波方法如卡爾曼濾波、擴(kuò)展卡爾曼濾波的基礎(chǔ)上發(fā)展而來的,在非線性、非高斯系統(tǒng)的模型參數(shù)估計中表現(xiàn)出明顯的優(yōu)勢,并且已被應(yīng)用于壽命預(yù)測領(lǐng)域[16]。

假設(shè)系統(tǒng)在離散時間序列tk=kΔt的狀態(tài)可以用式(8)的狀態(tài)傳遞函數(shù)描述。

式中,xk是系統(tǒng)在tk時刻的狀態(tài),fk是系統(tǒng)狀態(tài)傳遞函數(shù),θk是模型參數(shù)向量,ωk是系統(tǒng)噪聲。 系統(tǒng)狀態(tài)值與觀測值之間的關(guān)系為

式中,zk是系統(tǒng)在tk時刻的觀測值,hk是系統(tǒng)的觀測函數(shù),vk是觀測噪聲。 首先根據(jù)k-1 時刻的狀態(tài)對k時刻的先驗(yàn)概率密度函數(shù)進(jìn)行預(yù)測

當(dāng)測得新的觀測值后,對狀態(tài)概率密度函數(shù)進(jìn)行更新,得到k時刻的后驗(yàn)概率密度函數(shù)

用重要密度函數(shù)的離散采樣點(diǎn)和對應(yīng)權(quán)值來描述p(xk/z1:k),則式(11)轉(zhuǎn)化為

式中,δ(·)表示離散采樣點(diǎn),重要性權(quán)值wk采用式(13)進(jìn)行更新。

將式(13)條件概率增加參數(shù)項(xiàng)θk

利用新觀測值zk實(shí)現(xiàn)對系統(tǒng)狀態(tài)xk和模型參數(shù)θk所對應(yīng)權(quán)值的不斷更新。

1.3.2 獵人獵物優(yōu)化算法

獵人獵物優(yōu)化(hunter-prey optimizer, HPO)算法[17]通過模擬動物獵食過程,在搜索空間中按照一定的規(guī)則和策略對種群進(jìn)行引導(dǎo)和控制,不斷更新群體中每個獵人或獵物的位置,并用目標(biāo)函數(shù)評估出新的位置,對一個問題進(jìn)行尋優(yōu)。 該算法具有收斂快、尋優(yōu)能力強(qiáng)的特點(diǎn)。

首先獵人或獵物在搜索范圍內(nèi)按照式(15)隨機(jī)初始化位置。

式中,xi是獵人或獵物的位置,lb是問題變量的最小值,ub是問題變量的最大值,d是問題變量的維度。式(16)定義了搜索空間的下界和上界。

生成初始種群后,使用目標(biāo)函數(shù)計算每個解的適應(yīng)度值。 對于獵人的搜索機(jī)制,式(17)給出了其數(shù)學(xué)模型

式中,x(t)是當(dāng)前獵人位置,x(t+1)是獵人的下一次迭代位置,Ppos是獵物的位置,μ是所有位置的平均值,Z是由式(18)計算出的自適應(yīng)參數(shù)

其中,R1和R3是[0,1]內(nèi)的隨機(jī)向量,P是R1＜C的索引值,R2是[0,1]內(nèi)的隨機(jī)數(shù),IDX是滿足條件(P= =0)的向量R1的索引值,C是探索和開發(fā)之間的平衡參數(shù),其值在迭代過程中從1 減小到0.02,具體計算如下。

其中,it是當(dāng)前迭代次數(shù),Max是最大迭代次數(shù)。 計算獵物的位置Ppos,再根據(jù)式(20)計算所有位置的平均值μ,然后計算與該平均位置的距離。

距離位置平均值最大的位置被視為獵物位置Ppos。 假設(shè)最佳安全位置是最佳全局位置,這將使獵物有更好的生存機(jī)會,獵人可能會選擇另一個獵物,由式(21)更新獵物位置

式中,x(t)是獵物當(dāng)前的位置,x(t+1)是獵物的下一次迭代位置,Tpos是全局最優(yōu)位置,Z是由式(18)計算出的自適應(yīng)參數(shù),R4是[-1,1]內(nèi)的隨機(jī)數(shù);cos 函數(shù)及其輸入?yún)?shù)允許下一個獵物位置在不同半徑和角度的全局最優(yōu)位置,并提高開發(fā)階段的性能。

為了選擇獵人和獵物,結(jié)合式(17)和式(21),R5是[0,1]內(nèi)的隨機(jī)數(shù),β為調(diào)節(jié)參數(shù),如果R5值小于β,搜索到的位置被視為獵人,下一個位置將用式(17)更新;如果R5值大于β,搜索到的位置被視為獵物,下一個位置將用式(21)更新。

1.3.3 獵人獵物優(yōu)化算法改進(jìn)粒子濾波

利用獵人獵物優(yōu)化算法改進(jìn)粒子濾波的重采樣過程,克服了粒子經(jīng)過多次迭代后權(quán)重退化的問題。將粒子濾波中的粒子先驗(yàn)狀態(tài)作為獵人獵物初始種群個體位置,利用迭代尋優(yōu)過程改善粒子的分布情況,將退化的粒子集優(yōu)化至高似然值,使大部分粒子都能集中在真實(shí)狀態(tài)附近,解決了常規(guī)粒子濾波算法中重采樣使用的函數(shù)是次優(yōu)的問題。

HPO-PF 算法的實(shí)現(xiàn)步驟如下。

1)設(shè)置HPO-PF 算法的參數(shù),包括粒子個數(shù)、種群規(guī)模、最大迭代次數(shù)。

2)初始化粒子集。 從重要性概率密度函數(shù)中隨機(jī)抽取N個粒子。

3)更新粒子位置計算各個粒子的適應(yīng)度值,將粒子集作為獵人或獵物,更新當(dāng)前最優(yōu)位置。 定義適應(yīng)度函數(shù)為

式中,zR當(dāng)前狀態(tài)的實(shí)際量測值,zP為預(yù)測值。

4)種群移動,更新獵人或獵物的位置。 根據(jù)式(17)或式(21)計算獵人或獵物當(dāng)前的位置。

5)判斷循環(huán)是否停止。 當(dāng)滿足最大迭代次數(shù)或達(dá)到最優(yōu)適應(yīng)度值時,循環(huán)結(jié)束,否則重復(fù)步驟4)。

6)根據(jù)式(14)計算各個粒子權(quán)重,最后得出估計狀態(tài)。

1.4 預(yù)測模型構(gòu)建與參數(shù)選擇

1.4.1 模型構(gòu)建

本文提出的滾動軸承剩余使用壽命預(yù)測方法流程圖見圖2,方法分為兩個階段:構(gòu)建退化預(yù)測指標(biāo)以及預(yù)測剩余壽命。

圖2 剩余使用壽命預(yù)測模型Fig.2 Predictive model of RUL

本文使用Paris-Erdogan 裂紋擴(kuò)展模型[18-19]作為狀態(tài)空間模型,該狀態(tài)空間模型能很好地描述軸承的退化過程,并且易于構(gòu)建,可滿足工程上關(guān)于預(yù)測的可操作性的要求。 其中故障尺寸增長率可表示為

式中,x為故障尺寸;N為材料的疲勞壽命;C0和m為與材料相關(guān)的常數(shù);ΔK為應(yīng)力強(qiáng)度因子,計算式為

式中,β為與材料相關(guān)的參數(shù)。

將式(24)代入式(23)可得到

假設(shè)量測值zk和狀態(tài)值xk存在線性關(guān)系,改寫式(25)成如下形式。

式中,A和n為需要根據(jù)觀測值進(jìn)行估計的未知參數(shù),vk為觀測噪聲,ωk為系統(tǒng)噪聲,式(26) 與式(27)即為狀態(tài)空間模型表達(dá)式。

預(yù)測滾動軸承剩余使用壽命的具體步驟如下。

1)構(gòu)建退化預(yù)測指標(biāo)

①獲取滾動軸承的原始振動信號,計算輸入的原始數(shù)據(jù)的信息熵指標(biāo)。

②使用1.2 節(jié)的HPTF-BL 特征處理方法對計算出的信息熵指標(biāo)進(jìn)行處理,得到退化預(yù)測指標(biāo)。

2)預(yù)測剩余使用壽命

①使用最小二乘擬合初始化狀態(tài)空間模型的參數(shù)A和n,運(yùn)用HPO-PF 算法對參數(shù)進(jìn)行更新迭代,軸承失效退化曲線由算法遞推預(yù)測得到。

②將預(yù)測出的失效退化曲線結(jié)合預(yù)先設(shè)定好的失效閾值曲線,得到滾動軸承在該時刻預(yù)測的最長壽命、中間壽命以及最短壽命,分別對應(yīng)預(yù)測值上限、參考預(yù)測值與預(yù)測值下限[20]。

1.4.2 預(yù)測結(jié)果評估

使用均方根誤差(RMSE)以及擬合優(yōu)度R2作為預(yù)測性能的評估指標(biāo),來驗(yàn)證預(yù)測模型的有效性。均方根誤差RMSE 定義如下。

用相關(guān)系數(shù)(即擬合優(yōu)度)R2來評估預(yù)測結(jié)果與滾動軸承的真實(shí)退化曲線的擬合程度。R2的值在0 ～1 之間,其值越接近于1,表明預(yù)測出的退化曲線的擬合效果越好,預(yù)測模型的準(zhǔn)確度越高。R2定義如下。

1.4.3 參數(shù)選擇

在通過HPTF-BL 方法處理特征指標(biāo)的過程中,需要對λ和m這兩個參數(shù)進(jìn)行選擇。λ是一個控制時間序列平滑性的指標(biāo),在進(jìn)行特征指標(biāo)處理時需要獲取特征指標(biāo)的主要退化趨勢,根據(jù)文獻(xiàn)[21],λ選擇為16。m表示窗口內(nèi)數(shù)據(jù)的多少,至少可選擇為3,至多可選擇為與已知序列內(nèi)數(shù)據(jù)量相同的值,m的值越大,得到的結(jié)果越具有代表性,但往往會忽略一些數(shù)據(jù)信息;m的值越小,得到的結(jié)果則更能反映出特征指標(biāo)的邊界信息,具有更好的準(zhǔn)確性,本文中m選擇為5。

在改進(jìn)的粒子濾波算法中,本文設(shè)置HPO 算法的種群大小為30,最大迭代次數(shù)為500。 粒子濾波算法中粒子數(shù)的取值對預(yù)測結(jié)果有較大影響,使用IMS2 的1 號軸承數(shù)據(jù)[14]作為驗(yàn)證數(shù)據(jù),討論不同粒子數(shù)對預(yù)測結(jié)果的影響。 選用已知數(shù)據(jù)量為367的數(shù)據(jù),使用上述兩種評估指標(biāo)對預(yù)測結(jié)果進(jìn)行評估,結(jié)果如表1 所示。

表1 不同粒子數(shù)對預(yù)測結(jié)果的影響Table 1 Effect of different particle numbers on the prediction results

通過表1 可知,當(dāng)粒子數(shù)為500 時,粒子濾波可獲得最優(yōu)的預(yù)測效果,故粒子數(shù)設(shè)置為500。

2 實(shí)例驗(yàn)證

2.1 實(shí)驗(yàn)數(shù)據(jù)

選用美國辛辛那提大學(xué)IMS 中心實(shí)驗(yàn)平臺上采集的兩組數(shù)據(jù)[14]對本文提出的方法進(jìn)行驗(yàn)證。兩組數(shù)據(jù)分別如下:數(shù)據(jù)集1(IMS1)中的3 號軸承運(yùn)行至失效的數(shù)據(jù),該組數(shù)據(jù)總共包含2 156 個數(shù)據(jù)文件,3 號軸承在失效實(shí)驗(yàn)結(jié)束時出現(xiàn)內(nèi)圈故障;IMS2 中的1 號軸承運(yùn)行至失效的數(shù)據(jù),該組數(shù)據(jù)共包含984 個數(shù)據(jù)文件,1 號軸承在失效實(shí)驗(yàn)結(jié)束時出現(xiàn)外圈故障。 每一個數(shù)據(jù)文件均對應(yīng)一個文件序號。

2.2 數(shù)據(jù)驗(yàn)證及分析

首先,對IMS1 中3 號軸承和IMS2 中1 號軸承x軸方向的原始振動信號進(jìn)行特征提取,計算出信息熵指標(biāo),如圖3 所示。

圖3 信息熵指標(biāo)Fig.3 Information entropy index

由文獻(xiàn)[22]可知在IMS1 中的數(shù)據(jù)文件序號為1 608、IMS2 中的數(shù)據(jù)文件序號為533 時,軸承發(fā)生性能退化。 工程應(yīng)用場景下,一旦滾動軸承發(fā)生早期故障,之后得到的數(shù)據(jù)便可作為剩余使用壽命預(yù)測的已知訓(xùn)練數(shù)據(jù),故本文分別選取IMS1 在數(shù)據(jù)文件序號為1 608、IMS2 在數(shù)據(jù)文件序號為533 之后的軸承數(shù)據(jù)的信息熵指標(biāo),采用HPTF-BL 方法進(jìn)行處理,處理后的結(jié)果如圖4 所示。 可以看出一維的健康指標(biāo)被分為上下邊界以及信息熵主要趨勢,數(shù)據(jù)長度與原健康指標(biāo)的數(shù)據(jù)長度保持一致,并將原本的健康指標(biāo)包含其中,消除了原本健康指標(biāo)存在的局部波動,解決了隨機(jī)噪聲對預(yù)測結(jié)果的影響,具有單調(diào)性好、波動小的特點(diǎn)。 得到的3 條曲線在之后的步驟中同時進(jìn)行預(yù)測,可得到基于信號自身的預(yù)測置信區(qū)間,增強(qiáng)了在預(yù)測趨勢步驟中置信區(qū)間設(shè)置的可解釋性,進(jìn)一步提高了RUL 預(yù)測結(jié)論的準(zhǔn)確性。

圖4 處理后的信息熵指標(biāo)Fig.4 Information entropy index after treatment

當(dāng)IMS1 的已知數(shù)據(jù)分別為數(shù)據(jù)文件序號1 608至數(shù)據(jù)文件序號1 830、1 930、2 030、2 130 時,使用HPO-PF 算法對已知數(shù)據(jù)進(jìn)行訓(xùn)練,迭代預(yù)測得到失效退化曲線。 根據(jù)文獻(xiàn)[11]設(shè)定從滾動軸承發(fā)生性能退化開始,小波信息熵指標(biāo)變化值超過0.2時,即為達(dá)到運(yùn)行狀態(tài)不滿意階段,IMS1 軸承在信息熵指標(biāo)為0.75 時開始退化,故IMS1 軸承的失效閾值設(shè)定為0.95。 IMS1 的3 號軸承在已知數(shù)據(jù)量不同情況下的預(yù)測結(jié)果如圖5 所示。

圖5 IMS1 的3 號軸承預(yù)測結(jié)果Fig.5 Predicted results for bearing No.3 in data set 1

當(dāng)IMS2 的已知數(shù)據(jù)分別為數(shù)據(jù)文件序號533至數(shù)據(jù)文件序號750、800、850、900 時,使用HPOPF 算法對已知數(shù)據(jù)進(jìn)行訓(xùn)練,迭代預(yù)測得到失效退化曲線。 IMS2 軸承在信息熵指標(biāo)達(dá)到0.08 時開始退化,故相應(yīng)的失效閾值設(shè)定為0.3。 IMS2 的1 號軸承在已知數(shù)據(jù)量不同情況下的預(yù)測結(jié)果如圖6所示。

由圖5、6 預(yù)測結(jié)果可知,隨著已知數(shù)據(jù)文件的增加,HPO-PF 算法預(yù)測出的趨勢越來越接近真實(shí)的信息熵指標(biāo),最終該滾動軸承的剩余使用壽命L可用式(30)計算得出。

式中,Np為預(yù)測出的軸承失效時刻對應(yīng)的數(shù)據(jù)文件序號;NF為開始預(yù)測時刻對應(yīng)的數(shù)據(jù)文件序號;tper為采樣時間,在本文案例中為10 min。 最后可根據(jù)式(30)計算出在已知數(shù)據(jù)量大小不同時軸承的剩余使用壽命,具體結(jié)果如表2、3 所示。

由表2 和圖5 可知,當(dāng)IMS1 的已知數(shù)據(jù)為數(shù)據(jù)文件序號1 608 ～1 830 時,代入模型訓(xùn)練的數(shù)據(jù)有向上的趨勢,所以預(yù)測得出的結(jié)果超前真實(shí)壽命值800 min 左右;但當(dāng)已知數(shù)據(jù)為數(shù)據(jù)文件序號1 608 ～2 030 時,HPO-PF 算法的預(yù)測值上限為2 780 min、參考預(yù)測值為2 700 min、預(yù)測值下限為2 620 min;當(dāng)已知數(shù)據(jù)為數(shù)據(jù)文件序號1 608 ～2 130時,預(yù)測結(jié)果與真實(shí)剩余壽命僅相差400 min 左右,越來越接近于真實(shí)壽命。 所得出的預(yù)測值上限、參考預(yù)測值和預(yù)測值下限可共同為企業(yè)預(yù)測維修時間點(diǎn)提供建議。

由表3 和圖6 可知,當(dāng)IMS2 的已知數(shù)據(jù)為數(shù)據(jù)文件序號533 ～750 和533 ～800 時,代入模型訓(xùn)練的數(shù)據(jù)略有向上的趨勢,但總體上看還是趨于平緩,所以預(yù)測得出的結(jié)果延遲真實(shí)壽命值都在4 000 min 以上,不超過5 000 min;當(dāng)已知數(shù)據(jù)為數(shù)據(jù)文件序號533 ～850 時,HPO-PF 算法的預(yù)測結(jié)果開始趨近于真實(shí)剩余壽命,只相差400 min 左右;當(dāng)已知數(shù)據(jù)為數(shù)據(jù)文件序號533 ～900 時,預(yù)測結(jié)果與真實(shí)剩余壽命僅相差100 min 左右,越來越接近真實(shí)壽命。 由以上結(jié)果可知,本文預(yù)測模型得出的結(jié)論可根據(jù)已知數(shù)據(jù)量的變化實(shí)時動態(tài)改變。

2.3 性能對比

為了進(jìn)一步驗(yàn)證本文提出的剩余壽命預(yù)測方法的優(yōu)越性,選取5 種方法在IMS2 上進(jìn)行比較和分析,即使用指數(shù)退化模型和Paris-Erdogan 模型作為狀態(tài)空間函數(shù)的粒子濾波算法(PF-E 算法和PF-P算法)、使用指數(shù)退化模型和Paris-Erdogan 模型作為狀態(tài)空間函數(shù)的粒子群算法優(yōu)化粒子濾波算法(PSO-PF-E 算法和PSO-PF-P 算法)以及使用指數(shù)退化模型作為狀態(tài)空間函數(shù)的獵人獵物算法優(yōu)化粒子濾波算法(HPO-PF-E 算法)。 以上5 種方法的基本參數(shù)都與本文所提方法保持一致,采用已知數(shù)據(jù)量為367 的組別進(jìn)行驗(yàn)證,只使用信息熵指標(biāo)的主要趨勢進(jìn)行預(yù)測,最后通過RMSE 和R2評價預(yù)測效果,得到的對比結(jié)果如表4 所示。

表4 本文所提方法與其他5 種方法的對比結(jié)果Table 4 Comparison between proposed method and five other methods

對比結(jié)果表明,使用Paris-Erdogan 模型作為狀態(tài)空間函數(shù)的方法整體上的預(yù)測效果都優(yōu)于使用指數(shù)退化模型作為狀態(tài)空間函數(shù)的方法;而在使用Paris-Erdogan 模型的方法中,本文所提方法不論從均方根誤差上還是擬合優(yōu)度(相關(guān)系數(shù))上的表現(xiàn)都優(yōu)于其他方法的預(yù)測效果。

3 結(jié)論

考慮到實(shí)際工廠的剩余使用壽命預(yù)測需求,本文提出了一種融合HP 趨勢濾波-邊界線(HPTFBL)、獵人獵物優(yōu)化算法改進(jìn)粒子濾波(HPO-PF)的滾動軸承剩余使用壽命預(yù)測方法,該方法從健康指標(biāo)上下邊界刻畫、剩余使用壽命預(yù)測方法兩個方面進(jìn)行了優(yōu)化。 實(shí)驗(yàn)與對比分析的結(jié)果表明,該方法針對不同時刻以前的歷史時間序列數(shù)據(jù),可預(yù)測獲得不同時刻下滾動軸承剩余使用壽命的預(yù)測值上限、參考預(yù)測值和預(yù)測值下限,具有預(yù)測累計誤差小、預(yù)測精度高的特點(diǎn),可為企業(yè)工廠的預(yù)測性維修決策提供依據(jù)。