郭潤夏，張國良，王 雨

（中國民航大學(xué)電子信息與自動化學(xué)院，天津 300300）

軸承的剩余壽命（RUL，remaining useful life）主要是通過提取振動信號的時域和頻域特征來進行分析和判斷的[1]，主要有兩種方法：物理模型方法和數(shù)據(jù)驅(qū)動方法[2]。物理模型方法是指根據(jù)物理原則建立軸承狀態(tài)數(shù)學(xué)模型[3]，再根據(jù)以往的觀測值預(yù)測軸承未來的狀態(tài)，如用自回歸的方法來預(yù)測回轉(zhuǎn)軸承的RUL[4]。數(shù)據(jù)驅(qū)動方法是通過大量的趨勢數(shù)據(jù)學(xué)習(xí)軸承的狀態(tài)，并預(yù)測其未來狀態(tài)，如引入深度神經(jīng)網(wǎng)絡(luò)跟蹤軸承狀態(tài)并預(yù)測RUL[5-6]。然而，物理模型方法和數(shù)據(jù)驅(qū)動方法都有其自身的局限性，一般而言，物理模型難以建立，而數(shù)據(jù)驅(qū)動方法需要大量的數(shù)據(jù)，浪費計算資源。

粒子濾波算法適用于基于物理模型的非線性、非高斯系統(tǒng)的狀態(tài)估計，其通過一組加權(quán)粒子來近似后驗概率密度，主要包括3 個步驟[7]：重要性采樣、權(quán)重計算、重采樣。到目前為止，粒子濾波算法已成功應(yīng)用于許多領(lǐng)域的RUL 預(yù)測[8-9]。但粒子濾波算法存在一個嚴(yán)重的問題，即粒子退化問題[10]，為了將粒子濾波算法有效地應(yīng)用于實際工程中，首先要解決粒子退化問題。

為解決上述問題，提出了一種改進的粒子濾波算法——基于改進煙花算法的粒子濾波（FPF，fireworks algorithm-based particle filter）算法，并通過實驗驗證了該算法的有效性。FPF 采用了新的重采樣過程，即對煙花算法進行改進，然后利用改進的煙花算法優(yōu)化粒子權(quán)值，最后通過高斯近似方法產(chǎn)生新粒子。改進后的粒子濾波算法能有效抑制粒子退化。在此基礎(chǔ)上，根據(jù)Paris-Erdogan 模型預(yù)測了軸承的RUL。

1 煙花算法

1.1 標(biāo)準(zhǔn)煙花算法

煙花算法[11]是一種新型的群智能尋優(yōu)算法，每個初始煙花將被視為可行解空間的一個解, 通過爆炸操作使每個煙花爆炸產(chǎn)生爆炸火花，然后選擇一部分煙花隨機變異，產(chǎn)生變異火花，最后,從煙花、爆炸火花和變異火花的集合中選擇下一代煙花。

煙花算法主要由4 個基本部分組成：爆炸操作、變異操作、映射操作和選擇策略。不失一般性，假設(shè)待解決的優(yōu)化問題目標(biāo)函數(shù)為

式中：X 是函數(shù)的獨立變量；Ω 是解空間的可行域。

1）爆炸操作

根據(jù)適應(yīng)度函數(shù)計算每個煙花Xi的適應(yīng)度值f（Xi）。適應(yīng)度越好的煙花產(chǎn)生的火花越多，適應(yīng)度越差的煙花產(chǎn)生的火花越少，從而達到局部搜索和全局搜索的平衡。煙花Xi的爆炸半徑Ai和火花數(shù)Si的計算式分別為

式中：A 和M 分別是用來調(diào)整爆炸幅度、爆炸數(shù)目的兩個常數(shù)；уmin和уmax是煙花的最小適應(yīng)度值和最大適應(yīng)度值；ε 是一個機器數(shù)。

同時，為了避免產(chǎn)生過多或過少的火花，采用以下規(guī)則來限制火花的數(shù)量，即

式中：a 和b 是兩個常數(shù)，a ＜b ＜1；round()是取整函數(shù)。

根據(jù)式（2）～式（4）計算每個煙花Xi在第k 維度爆炸產(chǎn)生的火花為

2）變異操作

隨機挑選部分煙花變異以增強種群多樣性，煙花Xi在第k 維變異產(chǎn)生的火花計算如下

式中g(shù)auss()為高斯函數(shù)。

3）映射操作

當(dāng)煙花Xi在第k 維進行爆炸、變異操作時，如果超出了煙花種群在該維度的界限，則進行如下操作

4）選擇策略

從原始煙花、火花中選擇出下一代煙花，在候選集合K 中，適應(yīng)度值最小的原始煙花或火花必被選擇成為下一代煙花，其余N-1 個煙花以輪盤賭的方式選擇，每個煙花或火花被選擇的概率為

1.2 改進的煙花算法

重采樣過程是粒子優(yōu)化過程，而可以用煙花算法來解決粒子優(yōu)化問題。為了使煙花算法更好地設(shè)計重采樣過程，將有針對性地對煙花算法進行改進，具體改進操作如下。

1）選擇性的爆炸操作

除了適應(yīng)度值最小的煙花之外，根據(jù)式（1）～式（5）對剩余的煙花執(zhí)行爆炸操作。

2）選擇初始下一代煙花

執(zhí)行完上一步操作后，根據(jù)適應(yīng)度值從初始煙花、火花中選出N 個初始下一代煙花，被選出的煙花需要滿足以下條件

3）選擇需變異的煙花

用煙花集合C 表示選出來的初始下一代煙花，將集合中的元素按適應(yīng)度值大小進行排序，并按照黃金分割比將N 個煙花劃分為兩個部分，選擇小比例部分的煙花，分割示意圖如圖1 所示。

圖1 分割示意圖Fig.1 Splitting diagram

假設(shè)煙花集合中有8 個煙花，先將8 個煙花按適應(yīng)度值從小到大排序，根據(jù)黃金分割比將其劃分為紅、綠兩部分，3 個紅色的煙花即是需要變異的煙花。

4）變異操作

對步驟3）中選擇出的煙花執(zhí)行變異操作，變異后的煙花與未變異的煙花最終形成下一代煙花，變異操作如下

當(dāng)爆炸火花、變異火花超出界限時，根據(jù)式（7）執(zhí)行映射操作。

2 粒子濾波

本節(jié)主要介紹基礎(chǔ)粒子濾波（Basic PF）[2，12-13]及引入改進的煙花算法來設(shè)計粒子濾波算法的重采樣過程。

2.1 基礎(chǔ)粒子濾波算法

一般而言，系統(tǒng)的空間模型如下

式中：xk，zk分別表示k 時刻的系統(tǒng)狀態(tài)和測量數(shù)據(jù)；vk-1，wk分別表示已知分布的k-1 時刻的過程噪聲和k 時刻的測量噪聲, 并假設(shè)過程噪聲與測量噪聲獨立不相關(guān)。

已知初始狀態(tài)分布為p（x0），根據(jù)式（11）和式（12）分別計算系統(tǒng)的狀態(tài)轉(zhuǎn)移概率密度函數(shù)、似然概率密度函數(shù)為

式中：p（xk|xk-1）表示在狀態(tài)xk-1下xk的狀態(tài)轉(zhuǎn)移概率密度函數(shù)；p（zk| xk）是在狀態(tài)xk下zk的似然概率密度函數(shù)；pv()為過程噪聲的概率密度函數(shù)；pw()為測量噪聲的概率密度函數(shù)。

粒子濾波算法[12-13]基于蒙特卡羅模擬思想，從重要性概率密度函數(shù)q（xk| xk-1，zk）中采集N 個樣本粒子，并計算這些粒子的權(quán)值集，然后用該組粒子及其對應(yīng)的權(quán)值來近似系統(tǒng)狀態(tài)的后驗概率密度，每個粒子的對應(yīng)權(quán)值計算式如下

根據(jù)式（15）和式（16），可以得到權(quán)值的迭代公式如下

計算得到粒子的權(quán)值后，歸一化權(quán)值

再計算系統(tǒng)的后驗概率密度

最后，可以得到系統(tǒng)的狀態(tài)估計為

由式（17）可看出，隨著迭代次數(shù)的增加，大部分粒子的權(quán)值會越來越小，這就是粒子退化現(xiàn)象，因此，有必要重新設(shè)計重采樣過程解決粒子的退化問題。

2.2 基于改進煙花算法的粒子濾波算法

改進煙花算法設(shè)計重采樣過程并給出了FPF 的具體操作步驟。

2.2.1 重采樣

重采樣是為了抑制粒子的退化，通常用有效粒子數(shù)Neff來衡量粒子的退化程度，有效粒子數(shù)越少，退化越嚴(yán)重，有效粒子數(shù)計算如下

重采樣的基本過程是從離散近似的后驗概率密度中重采樣產(chǎn)生一組新的粒子。重采樣的原則是在權(quán)重大的粒子附近產(chǎn)生新粒子，使新粒子也具有較大的權(quán)值。當(dāng)粒子出現(xiàn)退化，即有效粒子數(shù)小于設(shè)置的重采樣閥值Nthr時，經(jīng)過改進的煙花算法進行操作，在權(quán)重大的粒子附近產(chǎn)生較大權(quán)重的新粒子，從而增加了有效粒子數(shù)，滿足重采樣的原則要求，解決了基礎(chǔ)粒子濾波算法的粒子退化問題，提高了算法的估計精度。

粒子的權(quán)值稱為煙花權(quán)值，當(dāng)新的觀測數(shù)據(jù)到達時，得到帶有煙花權(quán)值的粒子集合當(dāng)出現(xiàn)粒子退化現(xiàn)象時，根據(jù)改進的煙花算法對煙花權(quán)值進行優(yōu)化，獲得下一代煙花權(quán)值，并通過高斯近似產(chǎn)生新的粒子。

在新的重采樣過程中，適應(yīng)度函數(shù)為

重采樣流程圖如下。

圖2 重采樣步驟Fig.2 Resampling steps

2.2.2 操作步驟

圖3 FPF 的執(zhí)行過程Fig.3 The execution process of FPF

步驟1初始化：

（1）從初始概率密度p（x0）采樣；

步驟2序慣重要性采樣：

步驟3當(dāng)Neff＜Nthr，重采樣及狀態(tài)估計：

2.3 FPF 仿真驗證

在本節(jié)中，采用經(jīng)典的非線性模型

驗證FPF 的有效性和可靠性[13]。式中wk服從gauss（0，1）分布。

初始粒子服從gauss（0.5，1）分布。粒子數(shù)N=200，重采樣閥值Nthr=N/2。同時，以均方根誤差（RMSE）來評價算法的性能，即

假設(shè)系統(tǒng)模型是精確的，根據(jù)狀態(tài)方程及初始真值x0=0.5 計算下一時刻的狀態(tài)值。為了展示FPF 能有效抑制粒子退化，提高狀態(tài)估計精度，將使用Basic PF 與之進行比較。圖4 為狀態(tài)估計仿真圖，圖5 和圖6 分別為有效粒子數(shù)和均方根誤差仿真圖。

圖4 狀態(tài)估計Fig.4 State estimation

圖5 有效粒子數(shù)Fig.5 Number of effective particles

由圖4 可看出，F(xiàn)PF 的狀態(tài)估計更接近真實值。圖5 中，Basic PF 的有效粒子數(shù)幾乎變?yōu)?，但FPF 的有效粒子數(shù)仍在100～200 之間。圖6 中，與Basic PF相比，F(xiàn)PF 的均方根誤差較小。綜上，F(xiàn)PF 可以在抑制粒子退化的同時增加有效粒子的數(shù)量，狀態(tài)估計更準(zhǔn)確。

圖6 均方根誤差Fig.6 Root mean square error

3 電液舵機軸承的RUL 預(yù)測

RUL 定義為從預(yù)測起點開始到預(yù)測損傷指標(biāo)達到某一閥值的時間。當(dāng)損傷指標(biāo)超過閥值時，認(rèn)為該部件不能再繼續(xù)長期使用。一般來說，在不同的領(lǐng)域，RUL的計算方法不同，因此，RUL 的計算定義如下

式中：d()是距離函數(shù)；xi是預(yù)測起點的狀態(tài)；xj是預(yù)測終點的狀態(tài)；Δt 是采樣間隔；xthreshold是預(yù)定義的閥值；ξ是一個可調(diào)的較小的常數(shù)。利用FPF 算法對軸承狀態(tài)進行估計后，可根據(jù)式（25）得到預(yù)測的RUL。

為了有效地預(yù)測軸承的RUL，提取振動加速度信號的時域特征均方根（RMS）作為損傷指標(biāo)[2，14]（該指標(biāo)也是軸承在k 時刻的健康狀態(tài)xk）；定義為

根據(jù)提取出的軸承特征狀態(tài)，建立軸承的退化模型[14]如下

式中β 和m 是模型參數(shù)。

當(dāng)加速度信號均方根值達到1.5 g 時，軸承振動信號波動劇烈，說明軸承磨損嚴(yán)重，因此，設(shè)置xthreshold=1.5 g。

實驗平臺由加速疲勞裝置、電控可編程邏輯控制器（PLC，programmable logic controller）裝置、電機、齒輪箱和傳感器組成，如圖7 所示。傳感器主要包括加速度傳感器、速度傳感器和扭矩傳感器。信號（如振動）由傳感器測量，并通過采集卡傳遞到計算機。

圖7 實驗平臺Fig.7 Experimental platform

實驗平臺給出不同時間段的FPF 的預(yù)測結(jié)果，同時給出Basic PF、相對熵粒子濾波（MREIS-PF）的預(yù)測結(jié)果對比。圖8 是FPF 對3 個時間點損傷指標(biāo)的預(yù)測；圖9 為FPF、Basic PF 和MREIS-PF 在同一時間起點的狀態(tài)預(yù)測；圖10 為FPF、Basic PF 和MREIS-PF在整個過程中多個時間起點預(yù)測的RUL，并與真實RUL 進行對比。

圖8 FPF 在3 個時間起點的狀態(tài)預(yù)測Fig.8 State prediction by FPF at three time starting points

圖9 不同算法在同一時間起點的狀態(tài)預(yù)測Fig.9 State prediction by different algorithms at the same time starting point

圖10 剩余壽命預(yù)測結(jié)果Fig.10 Prediction of RUL

表1 給出了基于不同算法的多組軸承的RUL 預(yù)測和誤差率指標(biāo)，誤差率指標(biāo)定義如下

表1 不同算法的多組軸承的預(yù)測結(jié)果Tab.1 Predicted results of multiple sets of bearings obtained from different algorithms

由圖8～圖10 和表1 可看出：由于粒子的退化，Basic PF 在狀態(tài)預(yù)測期間波動較大且不夠準(zhǔn)確，導(dǎo)致預(yù)測的軸承RUL 與真實值相差較大；由于引入了設(shè)計的重采樣過程，F(xiàn)PF 在狀態(tài)跟蹤過程中更加穩(wěn)定，比Basic PF 和MREIS-PF 有更好的預(yù)測結(jié)果；在預(yù)測的早期階段，由于軸承退化過程緩慢，F(xiàn)PF 預(yù)測的RUL時間要比實際RUL 時間長，在預(yù)測的中間階段，軸承的磨損速率比早期階段快，F(xiàn)PF 預(yù)測的結(jié)果更接近真實值，當(dāng)軸承進入預(yù)測的后期階段，軸承退化加劇,由于模型包含了緩慢退化過程的信息，不能完全適應(yīng)軸承的快速退化，導(dǎo)致預(yù)測值大于真實值。綜合分析各階段預(yù)測結(jié)果，F(xiàn)PF 能較好地跟蹤軸承的磨損狀態(tài)，預(yù)測軸承的RUL 與真實值比較接近，在預(yù)測軸承RUL 方面顯示出較強的優(yōu)勢。

4 結(jié)語

針對粒子濾波算法中的粒子退化問題，采用改進的煙花算法設(shè)計重采樣過程，提出了一種改進的粒子濾波算法，可以減小粒子權(quán)值方差，抑制粒子退化。根據(jù)Paris-Erdogan 模型建立的物理模型能較好地適應(yīng)軸承的退化過程。同時，基于振動加速度信號提取的軸承狀態(tài)的時域特征均方根值能較好地反映軸承的磨損程度。在此基礎(chǔ)上，利用所提出的粒子濾波算法對電液舵機軸承的RUL 進行預(yù)測，取得了良好的效果。所提出的基于物理模型的FPF 預(yù)測電液舵機軸承RUL 的方法具有較好的實時性。然而，物理模型的參數(shù)會影響預(yù)測RUL 的精度，如何更好地辨識模型參數(shù)，從而獲得更精確的軸承RUL 值，需要進一步的研究。