許 雨 晨， 李 宏 坤， 馬 躍， 黃 剛 勁， 張 明 亮

( 大連理工大學(xué) 機(jī)械工程學(xué)院， 遼寧 大連 116024 )

0 引 言

工業(yè)領(lǐng)域中高端裝備的系統(tǒng)日趨復(fù)雜，綜合化、智能化程度不斷提升，運(yùn)行過程中維護(hù)和保障的成本也越來越高．由于設(shè)備組成環(huán)節(jié)的增加及各環(huán)節(jié)之間的高耦合性，發(fā)生故障和性能失效的可能性逐漸上升，因此重大裝備的預(yù)測和健康管理(prognostics and health management，PHM)技術(shù)逐漸成為研究熱點(diǎn)[1]．軸承作為工業(yè)領(lǐng)域的“心臟”，在各種機(jī)械設(shè)備中使用廣泛，是保障設(shè)備正常運(yùn)行的核心零件，一旦軸承發(fā)生故障，將對(duì)整個(gè)設(shè)備的運(yùn)行造成影響，導(dǎo)致安全事故和經(jīng)濟(jì)損失的發(fā)生，因此研究軸承的壽命預(yù)測方法具有重要價(jià)值．

文獻(xiàn)中指出大致可以將機(jī)械設(shè)備的壽命預(yù)測方法分為基于失效物理的方法、基于數(shù)據(jù)驅(qū)動(dòng)的方法和融合方法[2]．基于失效物理的方法根據(jù)機(jī)械的故障機(jī)制或損傷原理建立數(shù)學(xué)模型來描繪機(jī)械的退化狀態(tài)，具有較高的預(yù)測準(zhǔn)確度．但隨著機(jī)械設(shè)備日趨復(fù)雜化，很難根據(jù)故障機(jī)理建立精確的模型，因此這種方法的應(yīng)用領(lǐng)域具有較大的局限性．近年來隨著機(jī)器學(xué)習(xí)的發(fā)展，基于數(shù)據(jù)驅(qū)動(dòng)的方法成為壽命預(yù)測領(lǐng)域中熱門的研究方向．?dāng)?shù)據(jù)驅(qū)動(dòng)的方法不依賴專業(yè)的經(jīng)驗(yàn)知識(shí)，也不需要建立精確的數(shù)學(xué)模型，而是從歷史數(shù)據(jù)中自動(dòng)學(xué)習(xí)設(shè)備的退化信息，從而對(duì)未來趨勢進(jìn)行評(píng)估和預(yù)測．主流方法有神經(jīng)網(wǎng)絡(luò)、回歸模型、支持向量機(jī)、模糊理論、卡爾曼濾波和粒子濾波(particle filter，PF)等．融合方法是通過兩種及兩種以上方法的結(jié)合，將多種方法的優(yōu)勢進(jìn)行整合的預(yù)測方法．

雖然目前神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)等傳統(tǒng)數(shù)據(jù)驅(qū)動(dòng)的方法在壽命預(yù)測領(lǐng)域中已得到廣泛研究及應(yīng)用，但是在處理長期預(yù)測的不確定性問題上仍存在缺陷．針對(duì)這個(gè)問題，專家學(xué)者們提出使用貝葉斯濾波算法來解決[3]，其中粒子濾波是應(yīng)用最為廣泛的方法之一．基于貝葉斯濾波框架的壽命預(yù)測能夠給出剩余壽命的概率分布，改善長期預(yù)測的不確定性問題，并且相較于單一的預(yù)測值，更具參考價(jià)值[4]．Baraldi等[5]將粒子濾波與人工神經(jīng)網(wǎng)絡(luò)結(jié)合，對(duì)疲勞裂紋擴(kuò)展進(jìn)行了有效預(yù)測．Deng等[6]使用粒子濾波框架對(duì)鋰電池壽命進(jìn)行了預(yù)測，得到了較好的預(yù)測結(jié)果．Peng等[7]以飛機(jī)發(fā)動(dòng)機(jī)數(shù)據(jù)集為例，通過深度置信網(wǎng)絡(luò)構(gòu)建了健康指標(biāo)，采用粒子濾波算法對(duì)其壽命進(jìn)行了預(yù)測，驗(yàn)證結(jié)果表明該方法在壽命預(yù)測的不確定性表征方面具有較好的性能．雖然粒子濾波在處理長期預(yù)測的不確定性問題上有優(yōu)勢，但其算法本身仍有缺陷[8]，需要深入研究．

數(shù)據(jù)驅(qū)動(dòng)的軸承壽命預(yù)測方法中一個(gè)關(guān)鍵環(huán)節(jié)是確定開始預(yù)測的時(shí)刻[9](time to start prediction，TSP)，因?yàn)門SP對(duì)壽命預(yù)測的準(zhǔn)確性有重要影響．目前檢測TSP的常用方法有主觀確定、3σ方法、系統(tǒng)最長時(shí)間常數(shù)等，大多使用的是統(tǒng)計(jì)方法，不能針對(duì)某個(gè)系統(tǒng)自適應(yīng)檢測．

為解決上述問題，本文提出一種基于退化檢測和優(yōu)化粒子濾波[10]的軸承壽命預(yù)測方法．該方法首先利用軸承原始振動(dòng)信號(hào)構(gòu)建一個(gè)能夠全面描述軸承衰退過程的退化指標(biāo)，并通過所建立的退化指標(biāo)對(duì)軸承進(jìn)行自適應(yīng)退化檢測，確定開始預(yù)測的時(shí)刻，然后使用粒子群優(yōu)化算法對(duì)常規(guī)粒子濾波進(jìn)行改進(jìn)，提高算法的預(yù)測精度，應(yīng)用于軸承的壽命預(yù)測．

1 粒子濾波

1.1 粒子濾波算法

粒子濾波是一種基于蒙特卡羅采樣的貝葉斯估計(jì)算法．該算法的本質(zhì)是使用一系列具有對(duì)應(yīng)權(quán)值的隨機(jī)樣本點(diǎn)近似系統(tǒng)狀態(tài)的后驗(yàn)概率密度，從而得到系統(tǒng)狀態(tài)的最小方差估計(jì)[11]．該算法在非線性、非高斯系統(tǒng)的狀態(tài)參數(shù)估計(jì)中具有明顯優(yōu)勢，目前已經(jīng)廣泛應(yīng)用于機(jī)械設(shè)備及鋰電池等壽命預(yù)測研究中[12]．本文通過建立軸承的動(dòng)態(tài)模型，并基于粒子濾波算法對(duì)模型參數(shù)進(jìn)行估計(jì)，從而得到未來時(shí)刻軸承的系統(tǒng)狀態(tài)．

一般系統(tǒng)的動(dòng)態(tài)模型包括狀態(tài)轉(zhuǎn)移方程和量測方程，可以表示為

xk=fk(xk-1，vk-1)
zk=hk(xk，wk)

(1)

式中：xk表示系統(tǒng)在k時(shí)刻的狀態(tài)變量；fk：Rnx×Rnv→Rnx表示系統(tǒng)的狀態(tài)轉(zhuǎn)移方程；zk表示系統(tǒng)在k時(shí)刻的量測值；hk：Rnx×Rnw→Rnz表示系統(tǒng)的量測方程；vk和wk分別表示系統(tǒng)的過程噪聲和量測噪聲．

貝葉斯濾波算法是對(duì)系統(tǒng)狀態(tài)的后驗(yàn)分布進(jìn)行求解，首先利用動(dòng)態(tài)模型對(duì)k時(shí)刻系統(tǒng)狀態(tài)的先驗(yàn)分布進(jìn)行預(yù)測：

(2)

式中：p(xk|xk-1)由系統(tǒng)的狀態(tài)轉(zhuǎn)移方程定義．

當(dāng)新的量測值到來時(shí)，對(duì)先驗(yàn)分布進(jìn)行更新，得到k時(shí)刻的后驗(yàn)分布：

(3)

式中：p(zk|xk)由系統(tǒng)的量測方程定義；p(zk|z1：k-1)為歸一化常數(shù)，具有如下形式：

(4)

式(2)和(3)的計(jì)算過程涉及大量髙維積分運(yùn)算，難以求解，因此采用蒙特卡羅采樣來解決這個(gè)問題．通過蒙特卡羅方法可以將k時(shí)刻的后驗(yàn)分布離散加權(quán)為

(5)

(6)

(7)

由于常規(guī)算法中建議密度函數(shù)使用的是系統(tǒng)狀態(tài)的先驗(yàn)分布，舍棄了當(dāng)前時(shí)刻的量測值，因此出現(xiàn)了一些問題[10]．尤其當(dāng)似然函數(shù)很窄或者位于先驗(yàn)概率分布的尾部時(shí)，經(jīng)過幾次遞歸，大部分粒子權(quán)值將會(huì)變得很小，粒子權(quán)值的方差將會(huì)增大，即出現(xiàn)粒子退化現(xiàn)象．粒子退化現(xiàn)象是粒子濾波無法避免的問題，致使該方法在實(shí)際應(yīng)用中不能達(dá)到預(yù)期的濾波效果．為改善這個(gè)缺陷，方法之一是驅(qū)動(dòng)粒子向高似然區(qū)域移動(dòng)，因此本文引入粒子群優(yōu)化算法，對(duì)重要性采樣過程進(jìn)行優(yōu)化．

1.2 粒子群優(yōu)化粒子濾波算法

1.2.1 粒子群優(yōu)化算法 粒子群優(yōu)化算法是一種群智能尋優(yōu)算法，該算法模擬鳥類的捕食行為，通過鳥群中個(gè)體之間傳遞信息來搜索最優(yōu)解．該算法理論簡單、易于實(shí)現(xiàn)，因此在相關(guān)領(lǐng)域中的研究應(yīng)用極為廣泛．

粒子群優(yōu)化算法實(shí)現(xiàn)的基本思路為

(1)隨機(jī)初始化一群粒子，樣本容量為M，這群粒子僅具有速度v和位置x兩個(gè)屬性．

(2)定義適應(yīng)度函數(shù)，計(jì)算出每個(gè)粒子的適應(yīng)度值，找到每個(gè)粒子當(dāng)前的個(gè)體最優(yōu)值Pbest；通過信息共享找到整個(gè)群體中最優(yōu)的粒子，并將最優(yōu)粒子的個(gè)體最優(yōu)值Pbest作為全局最優(yōu)值Gbest．

(3)利用下式對(duì)粒子的速度和位置進(jìn)行更新：

vi=ω×vi+c1×r×(Pbesti-xi)+c2×r×

(Gbest-xi)

(8)

xi=xi+vi

(9)

式中：vi表示每個(gè)粒子的速度；ω稱為慣性權(quán)重；c1、c2統(tǒng)稱為學(xué)習(xí)因子；r是分布于(0，1)的隨機(jī)數(shù)；xi表示每個(gè)粒子的位置．

(4)返回步驟(2)，直至滿足結(jié)束條件．

為提高粒子群優(yōu)化算法的尋優(yōu)性能，本文使用參數(shù)動(dòng)態(tài)調(diào)整的速度更新公式．粒子群優(yōu)化算法傾向于在運(yùn)行初期進(jìn)行廣泛搜索，在運(yùn)行后期進(jìn)行局部求解．因此使用動(dòng)態(tài)調(diào)整的慣性權(quán)重，可以改善算法性能．本文使用線性慣性權(quán)重，其值根據(jù)迭代次數(shù)動(dòng)態(tài)變化，計(jì)算公式如下：

(10)

式中：ωmax、ωmin分別為慣性權(quán)重的最大值和最小值，t為當(dāng)前迭代次數(shù)，tmax為最大迭代次數(shù)．

學(xué)習(xí)因子分為個(gè)體學(xué)習(xí)因子c1和群體學(xué)習(xí)因子c2，c1表示粒子向自身經(jīng)歷的最優(yōu)值學(xué)習(xí)的能力，c2代表粒子向群體經(jīng)歷的最優(yōu)值學(xué)習(xí)的能力．對(duì)學(xué)習(xí)因子進(jìn)行動(dòng)態(tài)調(diào)整，也能改善算法的求解性能．因此本文算法中使用了動(dòng)態(tài)學(xué)習(xí)因子，迭代初期，要使粒子在全局范圍內(nèi)進(jìn)行廣泛搜索，增加粒子的多樣性，則需要設(shè)置較大的個(gè)體學(xué)習(xí)因子c1；而搜索后期則需要設(shè)置較大的群體學(xué)習(xí)因子c2，來加快算法的收斂速度．本文使用的動(dòng)態(tài)學(xué)習(xí)因子計(jì)算公式如下：

(11)

式中：cmax、cmin分別為學(xué)習(xí)因子的最大值和最小值．

1.2.2 粒子群優(yōu)化粒子濾波算法 通過1.1節(jié)內(nèi)容的介紹可知，常規(guī)粒子濾波算法中重要性采樣使用的建議密度函數(shù)是次優(yōu)的，導(dǎo)致算法濾波性能不夠精確，因此本文對(duì)重要性采樣過程進(jìn)行改進(jìn)，形成粒子群優(yōu)化粒子濾波算法．常規(guī)算法中采樣的粒子樣本分布與真實(shí)的后驗(yàn)概率分布差距較大，無法對(duì)系統(tǒng)狀態(tài)進(jìn)行有效估計(jì)．而粒子群優(yōu)化算法引入了當(dāng)前觀測值，能夠驅(qū)動(dòng)粒子向高似然區(qū)域移動(dòng)，所以優(yōu)化的結(jié)果是大部分粒子都能分布在真實(shí)狀態(tài)附近，從而提高了濾波精度[10]．

改進(jìn)算法的實(shí)現(xiàn)步驟如下：

(12)

利用1.2.1節(jié)中所描述的粒子群優(yōu)化算法，對(duì)從先驗(yàn)分布中采樣的粒子分布進(jìn)行優(yōu)化，然后根據(jù)下式計(jì)算重要性權(quán)值并歸一化：

(13)

(14)

歸一化權(quán)值為

(15)

(3)重采樣．重新采樣，并將所有粒子權(quán)值重置為1/M．

(4)輸出．使用下式對(duì)系統(tǒng)狀態(tài)進(jìn)行估計(jì)：

(16)

(5)k=k+1，新的量測值到來，返回第(2)步，直至退出算法．

2 軸承退化檢測

基于數(shù)據(jù)驅(qū)動(dòng)的軸承壽命預(yù)測過程中，有個(gè)關(guān)鍵問題需要解決，即確定開始預(yù)測的時(shí)間(TSP)．檢測軸承開始出現(xiàn)退化跡象的時(shí)刻，將其作為TSP．正確檢測TSP可以改善軸承壽命預(yù)測的性能，因此本文引入一種軸承退化檢測方法，分為以下兩步：構(gòu)建退化指標(biāo)和退化檢測．

2.1 構(gòu)建退化指標(biāo)

本文通過特征融合構(gòu)建退化指標(biāo)，具體步驟如下：

(1)特征提取

本文分別提取時(shí)域、時(shí)頻域和基于三角函數(shù)的特征組成候選特征集．提取的候選特征如下．時(shí)域特征：F1，均值；F2，標(biāo)準(zhǔn)差；F3，最大值；F4，最小值；F5，均方根；F6，峭度；F7，峰值；F8，峰峰值；F9，偏斜度；F10，熵值；F11，峰值指標(biāo)；F12，脈沖指標(biāo)；F13，峭度指標(biāo)；F14，裕度指標(biāo)；F15，波形指標(biāo)；F16，偏斜度指標(biāo)．時(shí)頻域特征：F17～F24，IMF能量值．三角函數(shù)的特征：F25，IHC標(biāo)準(zhǔn)差；F26，IHS標(biāo)準(zhǔn)差．時(shí)頻域特征提取方法是對(duì)原始振動(dòng)信號(hào)進(jìn)行集成經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition，EEMD)得到8個(gè)本征模態(tài)分量(intrinsic mode function，IMF)，計(jì)算每個(gè)IMF的能量作為特征值．兩個(gè)三角函數(shù)特征分別為反三角雙曲余弦(inverse hyperbolic cosine，IHC)標(biāo)準(zhǔn)差和反三角雙曲正弦(inverse hyperbolic sine，IHS)標(biāo)準(zhǔn)差[13]．此類三角函數(shù)可產(chǎn)生低尺度特征，并具有更好的趨勢性，計(jì)算公式如下：

(17)

(2)特征優(yōu)選

候選特征集中的退化特征對(duì)故障發(fā)展的敏感度有很大差別，為進(jìn)行有效預(yù)測，需選取敏感度高的特征用于構(gòu)建退化指標(biāo)．本文首先使用時(shí)間相關(guān)性、單調(diào)性和魯棒性對(duì)候選特征進(jìn)行評(píng)估．相關(guān)性表征了特征與時(shí)間序列的線性相關(guān)度，單調(diào)性用于描述特征持續(xù)增大或減小的趨勢，魯棒性反映特征對(duì)異常值或噪聲的容忍度．

采用移動(dòng)平均法對(duì)特征進(jìn)行平滑處理，得到退化特征的趨勢部分和殘差部分：

X(tk)=Xt(tk)+Xr(tk)

(18)

式中：X(tk)為退化特征，Xt(tk)為退化特征的趨勢值，Xr(tk)為退化特征的隨機(jī)值．

(19)

(20)

(21)

式中：K是原始數(shù)據(jù)長度，tk表示時(shí)間序列．

單一的評(píng)價(jià)指標(biāo)只能從特定角度衡量退化特征對(duì)故障發(fā)展的敏感度，并且對(duì)于同一個(gè)退化特征，評(píng)價(jià)指標(biāo)間也可能出現(xiàn)不一致的情況．因此需要將3個(gè)評(píng)價(jià)指標(biāo)綜合起來，對(duì)退化特征進(jìn)行優(yōu)選[14]．所以本文將3個(gè)指標(biāo)的線性組合作為綜合評(píng)價(jià)指標(biāo)，計(jì)算公式如下：

s.t.ωi>0

(22)

式中：J為綜合評(píng)價(jià)指標(biāo)；Ω為候選特征集；ωi為單一評(píng)價(jià)指標(biāo)的權(quán)重，本文中選取ω1=0.4，ω2=0.3，ω3=0.3．從上式中可以看出綜合評(píng)價(jià)指標(biāo)與單一評(píng)價(jià)指標(biāo)具有正相關(guān)關(guān)系，因此綜合指標(biāo)值越大，說明退化特征的敏感度越高，比較適用于壽命預(yù)測．

(3)特征融合

通過上文中特征選擇步驟，得到一個(gè)包含P個(gè)特征的優(yōu)選特征集，為構(gòu)建一個(gè)能夠較為全面反映軸承故障信息并具有較好趨勢性的指標(biāo)，本文使用加權(quán)融合方法將優(yōu)選特征集中的退化特征進(jìn)行融合，建立退化指標(biāo)．方法步驟如下：

步驟1計(jì)算每個(gè)優(yōu)選特征的權(quán)重．使用相關(guān)性評(píng)價(jià)指標(biāo)作為加權(quán)融合算法中的權(quán)重值，并將其歸一化：

(23)

式中：Ci是第i個(gè)特征的相關(guān)性評(píng)價(jià)指標(biāo)．

步驟2計(jì)算優(yōu)選特征集中每個(gè)特征的前s個(gè)點(diǎn)的均值Vi，作為軸承健康階段的特征值：

(24)

式中：Xi，k為第i個(gè)特征k時(shí)刻的值．

步驟3計(jì)算健康階段特征值和當(dāng)前時(shí)刻特征值之間的距離，然后使用步驟1中的歸一化權(quán)重進(jìn)行加權(quán)融合：

(25)

式中：Dt為t時(shí)刻退化指標(biāo)的值．

使用局部加權(quán)回歸算法平滑處理構(gòu)建的退化指標(biāo)，并將時(shí)間間隔設(shè)為1 min，進(jìn)行重采樣，得到最終的退化指標(biāo)．

2.2 退化檢測

上文中建立的退化指標(biāo)反映了軸承的運(yùn)行狀況，軸承運(yùn)行正常時(shí)，退化指標(biāo)變化平穩(wěn)，隨著運(yùn)行時(shí)間的增長，軸承逐漸出現(xiàn)損傷，退化指標(biāo)則會(huì)反映出軸承的故障演化過程．一旦軸承出現(xiàn)退化跡象，指標(biāo)值就會(huì)出現(xiàn)較大幅度的增長，說明軸承進(jìn)入了退化階段，需要從此時(shí)開始進(jìn)行壽命預(yù)測．退化檢測在壽命預(yù)測工作中是極為關(guān)鍵的環(huán)節(jié)，錯(cuò)誤識(shí)別的TSP會(huì)導(dǎo)致遺漏與早期故障相關(guān)的有用信息，或者將正常運(yùn)行的數(shù)據(jù)并入預(yù)測結(jié)果中[15]．這些情況都會(huì)影響預(yù)測結(jié)果的精度，因此確定準(zhǔn)確的TSP具有重要意義，然而目前檢測TSP的方法大多使用的是統(tǒng)計(jì)方法，常用方法有3σ方法、系統(tǒng)最長時(shí)間常數(shù)等，這些方法都無法針對(duì)特定系統(tǒng)自適應(yīng)調(diào)整，為此本文提出一種自適應(yīng)退化檢測法，通過衡量退化指標(biāo)的增長幅度來確定TSP．

首先確定一個(gè)窗口長度為l的滑動(dòng)窗，使用滑動(dòng)窗口中l(wèi)個(gè)樣本計(jì)算退化指標(biāo)的梯度值[16]，窗口中l(wèi)個(gè)樣本的梯度值可由下式得到：

(26)

式中：gi為第i個(gè)窗口的梯度值，Di為第i個(gè)窗口中l(wèi)個(gè)退化指標(biāo)的值．如果第i個(gè)窗口中計(jì)算的gi大于設(shè)置的閾值，就可以認(rèn)為軸承開始退化．圖1為自適應(yīng)檢測TSP的方法．

圖1 TSP檢測方法

3 基于PSO-PF的軸承剩余壽命預(yù)測方法

本文首先構(gòu)建了一個(gè)加權(quán)融合的退化指標(biāo)，并檢測開始預(yù)測的時(shí)刻(TSP)，然后利用PSO-PF算法從檢測到的TSP開始對(duì)軸承剩余壽命進(jìn)行預(yù)測．方法流程圖如圖2所示．

圖2 方法流程圖

本文使用Pairs-Erdogan裂紋擴(kuò)展模型[17]描述軸承的退化過程，故障尺寸增長率可由下式表示：

(27)

式中：x為故障尺寸；N為材料的疲勞壽命；C0和m為與材料相關(guān)的常數(shù)；ΔK為應(yīng)力強(qiáng)度因子，表達(dá)式如下：

(28)

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

將式(28)代入式(27)可得

(29)

假設(shè)量測值zk和狀態(tài)值xk存在線性關(guān)系，根據(jù)一般系統(tǒng)的動(dòng)態(tài)模型，將上式改寫為如下形式：

(30)

式中：A和n為需根據(jù)觀測值進(jìn)行估計(jì)的未知參數(shù)．

建立了軸承退化模型后，便可利用1.2節(jié)中提出的粒子群優(yōu)化粒子濾波算法對(duì)軸承剩余壽命進(jìn)行預(yù)測，預(yù)測的具體步驟如下：

(1)退化檢測．首先使用局部加權(quán)回歸算法平滑處理退化指標(biāo)，并進(jìn)行等距重采樣；然后利用自適應(yīng)退化檢測法確定軸承首次預(yù)測的時(shí)間點(diǎn)．

(2)將已知數(shù)據(jù)序列代入狀態(tài)轉(zhuǎn)移方程中，對(duì)模型參數(shù)A和n進(jìn)行非線性最小二乘擬合，得到模型參數(shù)的初始值．

(3)使用改進(jìn)的粒子濾波算法對(duì)模型參數(shù)不斷調(diào)整，待更新結(jié)束后，將最新的估計(jì)值作為模型參數(shù)，對(duì)軸承未來狀態(tài)進(jìn)行遞推預(yù)測，對(duì)剩余壽命分布進(jìn)行估計(jì)．

4 實(shí)驗(yàn)及結(jié)果分析

4.1 實(shí)驗(yàn)介紹

本文的實(shí)驗(yàn)數(shù)據(jù)來源于FEMTO-ST研究所的PRONOSTIA軸承實(shí)驗(yàn)臺(tái)進(jìn)行的軸承加速退化壽命實(shí)驗(yàn)，實(shí)驗(yàn)臺(tái)結(jié)構(gòu)如圖3所示．實(shí)驗(yàn)中為加快軸承退化速度，使用了加載裝置對(duì)其施加徑向力，將其全壽命周期縮短到數(shù)小時(shí)之內(nèi)．該實(shí)驗(yàn)中軸承轉(zhuǎn)速為1 800 r/min，采樣頻率為25.6 kHz，采樣時(shí)間為0.1 s，負(fù)載為4 000 N，每隔10 s進(jìn)行一次采樣．當(dāng)振動(dòng)幅值大于20g后，軸承視為失效．圖4為本文算法驗(yàn)證所用的軸承1和軸承2的振動(dòng)信號(hào)，從圖中可以看出軸承運(yùn)行前期振動(dòng)幅值很小且波動(dòng)平穩(wěn)，到運(yùn)行末期振動(dòng)幅值迅速增大，很短時(shí)間內(nèi)振動(dòng)幅值達(dá)到失效閾值．

圖3 PRONOSTIA軸承實(shí)驗(yàn)臺(tái)

(b) 軸承2振動(dòng)信號(hào)

4.2 軸承剩余壽命預(yù)測

從軸承1和軸承2的振動(dòng)信號(hào)中提取2.1節(jié)中所述的26個(gè)候選特征，使用綜合評(píng)價(jià)指標(biāo)J進(jìn)行評(píng)估，評(píng)估結(jié)果如表1所示．軸承1最優(yōu)特征和最差特征對(duì)比如圖5所示，軸承2最優(yōu)特征和最差特征對(duì)比如圖6所示．從圖中可以非常直觀地看出綜合評(píng)價(jià)較高的特征具有更好的時(shí)間相關(guān)性及單調(diào)性，而綜合評(píng)價(jià)較低的特征對(duì)故障的發(fā)展并不敏感，說明本文使用的特征選擇方法是有效的．對(duì)于軸承1和2分別選取前14個(gè)候選特征組成優(yōu)選特征集，對(duì)其進(jìn)行加權(quán)融合得到退化指標(biāo)；然后對(duì)其進(jìn)行平滑處理，并按時(shí)間間隔為1 min進(jìn)行重采樣．

表1 候選特征綜合評(píng)價(jià)

為驗(yàn)證融合后的特征更適用于壽命預(yù)測，對(duì)試驗(yàn)軸承的原始特征和融合特征的時(shí)間相關(guān)性進(jìn)行比較．從表2中可以看出融合特征的時(shí)間相關(guān)性較高，說明具有更好的趨勢性，同時(shí)融合特征又包含較全面的故障信息，因此相較于單一原始特征更適用于壽命預(yù)測．

(a) 最優(yōu)特征

(b) 最差特征

(a) 最優(yōu)特征

(b) 最差特征

表2 特征時(shí)間相關(guān)性對(duì)比

使用預(yù)處理后的退化指標(biāo)，利用前文所述的壽命預(yù)測方法，對(duì)軸承剩余壽命進(jìn)行預(yù)測．首先利用梯度值法檢測軸承開始退化的時(shí)間，確定首次預(yù)測時(shí)刻．根據(jù)梯度值計(jì)算結(jié)果，將梯度值出現(xiàn)較大幅度增長后的值設(shè)置為閾值，可得軸承1梯度閾值為0.001 5，確定的TSP為第321 min，軸承2梯度閾值為0.000 8，確定的TSP為第 271 min，然后從TSP開始對(duì)軸承進(jìn)行壽命預(yù)測．為證明本文提出的改進(jìn)粒子濾波預(yù)測算法的優(yōu)勢，使用常規(guī)粒子濾波算法進(jìn)行對(duì)比實(shí)驗(yàn)．

對(duì)于軸承1，首先使用軸承全壽命振動(dòng)數(shù)據(jù)構(gòu)建退化指標(biāo)，然后確定了TSP為第321 min，即將前320 min的退化指標(biāo)作為已知序列，將其代入狀態(tài)轉(zhuǎn)移方程，對(duì)模型中參數(shù)A和n進(jìn)行非線性最小二乘擬合，得到模型參數(shù)的初始值．然后使用本文提出的改進(jìn)粒子濾波算法對(duì)模型參數(shù)不斷調(diào)整，使用的訓(xùn)練數(shù)據(jù)為前320 min的退化指標(biāo)序列，待更新結(jié)束后，將最新的估計(jì)值作為模型參數(shù)，從第321 min起對(duì)軸承未來狀態(tài)進(jìn)行遞推預(yù)測，從而得到軸承的剩余壽命估計(jì)．軸承2壽命預(yù)測過程同理．

圖7和圖8分別為軸承1和軸承2的預(yù)測結(jié)果．圖7(a)和圖8(a)為使用常規(guī)粒子濾波進(jìn)行預(yù)測的結(jié)果，圖7(b)和圖8(b)為使用粒子群優(yōu)化粒子濾波算法預(yù)測的結(jié)果．通過對(duì)比兩圖，直觀上可以看出改進(jìn)算法的預(yù)測值更趨近于真實(shí)值，并且剩余壽命分布的置信區(qū)間更為集中，說明改進(jìn)算法的預(yù)測精度更高．

(a) PF預(yù)測結(jié)果

(b) PSO-PF預(yù)測結(jié)果

(a) PF預(yù)測結(jié)果

(b) PSO-PF預(yù)測結(jié)果

為量化說明改進(jìn)算法的優(yōu)勢，使用平均相對(duì)誤差(emap)、均方誤差(erms)和剩余壽命概率分布的區(qū)間寬度3個(gè)指標(biāo)對(duì)其進(jìn)行評(píng)價(jià)．分析構(gòu)建的退化指標(biāo)，將失效閾值設(shè)為0.5，則可以根據(jù)下式計(jì)算軸承的剩余使用壽命(Nr)：

Nr=tpre-tFPT

(31)

式中：tpre為預(yù)測的失效時(shí)刻，tFPT為開始預(yù)測的時(shí)刻．

emap和erms的計(jì)算公式如下：

(32)

(33)

剩余壽命概率分布的區(qū)間寬度是通過計(jì)算剩余壽命分布的95%置信區(qū)間得到的，區(qū)間寬度越小，說明預(yù)測的概率分布越集中，預(yù)測結(jié)果越準(zhǔn)確．

兩種預(yù)測方法的量化對(duì)比結(jié)果如表3所示．從預(yù)測誤差角度看，使用PSO-PF預(yù)測方法的平均相對(duì)誤差和均方誤差值都小于PF預(yù)測方法；并且從壽命分布來看，改進(jìn)算法預(yù)測的壽命概率分布的區(qū)間寬度更小，說明預(yù)測效果更好．因此本文使用的改進(jìn)粒子濾波算法能夠有效地提高軸承壽命預(yù)測的精度．

表3 預(yù)測結(jié)果對(duì)比

5 結(jié) 語

為準(zhǔn)確預(yù)測軸承剩余使用壽命，本文提出了一種粒子群優(yōu)化粒子濾波的軸承壽命預(yù)測方法．首先，為精確預(yù)測軸承剩余壽命，需建立能夠描繪軸承故障演化的退化指標(biāo)．本文通過特征提取、優(yōu)選、加權(quán)融合，獲得了全面反映軸承故障發(fā)展的退化指標(biāo)．并且通過自適應(yīng)退化檢測法得到首次預(yù)測時(shí)刻．最終使用粒子群優(yōu)化粒子濾波預(yù)測方法，從檢測到的首次預(yù)測時(shí)刻起對(duì)軸承進(jìn)行壽命預(yù)測．從與常規(guī)粒子濾波方法的對(duì)比結(jié)果來看，本文提出的改進(jìn)方法能夠有效提升粒子濾波算法的估計(jì)精度，更準(zhǔn)確地預(yù)測軸承剩余壽命．