基于煙花粒子濾波的電液舵機(jī)軸承剩余壽命預(yù)測(cè)

郭潤(rùn)夏，張國(guó)良，王 雨

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

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

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

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

1 煙花算法

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

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

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

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

1）爆炸操作

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

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

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

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

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

2）變異操作

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

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

3）映射操作

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

4）選擇策略

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

1.2 改進(jìn)的煙花算法

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

1）選擇性的爆炸操作

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

2）選擇初始下一代煙花

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

3）選擇需變異的煙花

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

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

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

4）變異操作

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

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

2 粒子濾波

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

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

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

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

已知初始狀態(tài)分布為p（x0），根據(jù)式（11）和式（12）分別計(jì)算系統(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()為測(cè)量噪聲的概率密度函數(shù)。

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

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

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

再計(jì)算系統(tǒng)的后驗(yàn)概率密度

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

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

2.2 基于改進(jìn)煙花算法的粒子濾波算法

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

2.2.1 重采樣

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

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

粒子的權(quán)值稱為煙花權(quán)值，當(dāng)新的觀測(cè)數(shù)據(jù)到達(dá)時(shí)，得到帶有煙花權(quán)值的粒子集合當(dāng)出現(xiàn)粒子退化現(xiàn)象時(shí)，根據(jù)改進(jìn)的煙花算法對(duì)煙花權(quán)值進(jì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)估計(jì)：

2.3 FPF 仿真驗(yàn)證

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

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

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

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

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

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

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

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

3 電液舵機(jī)軸承的RUL 預(yù)測(cè)

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

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

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

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

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

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

實(shí)驗(yàn)平臺(tái)由加速疲勞裝置、電控可編程邏輯控制器（PLC，programmable logic controller）裝置、電機(jī)、齒輪箱和傳感器組成，如圖7 所示。傳感器主要包括加速度傳感器、速度傳感器和扭矩傳感器。信號(hào)（如振動(dòng)）由傳感器測(cè)量，并通過采集卡傳遞到計(jì)算機(jī)。

圖7 實(shí)驗(yàn)平臺(tái)Fig.7 Experimental platform

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

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

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

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

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

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

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

4 結(jié)語

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