基于無跡卡爾曼濾波的離心泵剩余壽命研究*

王炳波 劉賡傳

(中國石油工程建設(shè)有限公司北京設(shè)計分公司)

0 引 言

泵類設(shè)備廣泛應(yīng)用于石油石化行業(yè)，為液體或氣體物料的輸送提供動力。為保證系統(tǒng)的安全穩(wěn)定運行，應(yīng)對泵類設(shè)備性能退化和壽命預(yù)測進行研究[1]?，F(xiàn)有的泵類設(shè)備壽命預(yù)測方法大致可以分為基于物理失效模型的方法和數(shù)據(jù)驅(qū)動的方法，基于數(shù)據(jù)驅(qū)動的方法又可分為機器學(xué)習(xí)方法、多元統(tǒng)計分析方法、特征量提取方法和信息融合方法等[2-3]。郭文琪等[4]綜述了礦井主排水泵壽命預(yù)測方法，指出泵類壽命衰退指標(biāo)應(yīng)該多樣化，只有多種預(yù)測模型有機結(jié)合才能保證預(yù)測模型的準(zhǔn)確性；劉銳等[5]從離心泵流量公式出發(fā)，結(jié)合流量與間隙的物理關(guān)系模型，建立了離心泵壽命評估物理模型；韓可[6]基于變分模態(tài)分解(VMD)和支持向量數(shù)據(jù)描述(SVDD)方法，建立了液壓泵性能退化模型；何慶飛等[7]建立了基于灰色支持向量機的液壓泵壽命預(yù)測方法；馬濟喬等[8]搭建了液壓泵加速退化試驗平臺，建立了液壓泵壽命預(yù)測方法，并進行了預(yù)測方法的可靠性分析?；谖锢砟Ｐ偷碾x心泵壽命預(yù)測方法，由于模型簡化和失效機理的多樣性，所以沒有普遍適用性；基于數(shù)據(jù)驅(qū)動的泵類設(shè)備預(yù)測方法不能根據(jù)實時監(jiān)測的數(shù)據(jù)狀態(tài)對退化模型進行實時更新。因此，本文基于無跡卡爾曼濾波方法，提出了一種結(jié)合已有同類設(shè)備的退化數(shù)據(jù)，并有機結(jié)合特定離心泵實時退化數(shù)據(jù)的離心泵剩余壽命預(yù)測方法。

1 理論基礎(chǔ)

1.1 離心泵性能退化模型的構(gòu)建

離心泵的壽命與多種因素有關(guān)，如葉輪和軸承磨損、密封失效等，在性能退化期間其揚程、功率和效率等參數(shù)會發(fā)生改變，因此可以用揚程、功率和效率等參數(shù)來表征離心泵的健康狀況。離心泵揚程H與流量Q的關(guān)系可以近似表示為如下的二次函數(shù)：

H=K1+K2Q-K3Q2

(1)

式中：K1、K2、K3為常數(shù)。

葉輪以及軸承磨損是離心泵常見的故障，磨損速率直接影響了離心泵的使用壽命[9]。葉輪磨損主要影響公式(1)中的常數(shù)項K1，軸承磨損主要影響泵轉(zhuǎn)速n，上述磨損規(guī)律可用如下公式表示：

(2)

同時考慮兩種磨損，即把公式(2)帶入公式(1)，經(jīng)過變換后可得[10]：

(3)

綜上，離心泵揚程的壽命衰退關(guān)系可表示為：

H=K1+K2Q-K3Q2-At2-Bt

(4)

式中：A、B為常數(shù)。

由于離心泵的功率和效率與離心泵揚程成正比，所以其功率和效率的性能衰退關(guān)系也可表示為二次函數(shù)的形式。為了綜合考慮多個性能衰退指標(biāo)，引入離心泵壽命健康指標(biāo)Hi。該指標(biāo)根據(jù)熵值法綜合考慮離心泵的揚程、功率和效率，因此離心泵運行時間t與健康指標(biāo)的關(guān)系可用如下二次函數(shù)來表示：

Hi=a′t2+b′t+c′

(5)

式中：a′、b′、c′為常數(shù)，a′<0。

1.2 支持向量機基本原理

(6)

根據(jù)KKT條件[14]，上式可轉(zhuǎn)化為凸二次優(yōu)化問題：

(7)

當(dāng)x與y為非線性關(guān)系時，可通過滿足Mercer條件的核函數(shù)將該非線性問題轉(zhuǎn)化為更高維度的線性問題，此時只需將式(7)中的xi·xj替換為核函數(shù)K(xi·xj)。常用的核函數(shù)為RBF核函數(shù)和Ploy核函數(shù)，本文采用RBF核函數(shù)(p為核函數(shù)的指數(shù))，其表達式為：

(8)

1.3 卡爾曼濾波基本原理

卡爾曼濾波(Kalman filtering，KF)是一種利用線性系統(tǒng)狀態(tài)方程，通過系統(tǒng)輸入輸出觀測數(shù)據(jù)，對系統(tǒng)狀態(tài)進行最優(yōu)估計的算法，經(jīng)過幾十年的發(fā)展已拓展至非線性系統(tǒng)的狀態(tài)評估，可應(yīng)用在軸承壽命預(yù)測、交通流預(yù)測、電池剩余壽命預(yù)測和機械零件壽命預(yù)測等方面[15-18]。無跡卡爾曼濾波(Unscented Kalman Filter，UKF)是一種利用采樣策略逼近非線性分布的濾波方法，其核心是無跡變換(Unscented Transform，UT)，即一種計算非線性變換中變量的統(tǒng)計特征的新方法，UKF既突破了卡爾曼濾波不能解決非線性問題的局限，又克服了拓展卡爾曼濾波(EKF)計算量較大的缺點。

卡爾曼濾波是一種特殊的貝葉斯濾波，是基于所有可用信息來構(gòu)造系統(tǒng)狀態(tài)的概率密度函數(shù)，即根據(jù)給定隨機變量Zk=[z1……zk]T、Uk=[u0……uk]和初始狀態(tài)x0的概率密度ρ0，估計出后驗概率密度ρ(xk|Zk,Uk)，然后根據(jù)觀測數(shù)據(jù)來計算狀態(tài)量取不同值時的置信度，進而得到狀態(tài)的最佳估計。其預(yù)測和更新原理為：假設(shè)k時刻狀態(tài)向量的后驗概率密度函數(shù)已知，則通過系統(tǒng)狀態(tài)方程可得到k+1時刻系統(tǒng)狀態(tài)向量的先驗概率密度函數(shù)[19]：

(9)

當(dāng)k+1時刻觀測向量已知時，可根據(jù)貝葉斯準(zhǔn)則用zk+1更新狀態(tài)向量的先驗概率密度函數(shù)：

(10)

式中：xk為系統(tǒng)的狀態(tài)向量，uk為系統(tǒng)的輸入向量。

2 預(yù)測模型的建立與求解

2.1 預(yù)測模型算法流程

基于UKF的離心泵壽命算法流程如圖1所示。根據(jù)同類離心泵的退化歷史數(shù)據(jù)，采用熵值法確定復(fù)合健康指標(biāo)，并確定離心泵失效時的健康指標(biāo)閾值，通過數(shù)據(jù)擬合得到離心泵卡爾曼濾波狀態(tài)方程，并得到k時刻與k-1時刻狀態(tài)方程的傳遞關(guān)系。對于需要預(yù)測壽命的離心泵，根據(jù)其歷史運行數(shù)據(jù)和出廠設(shè)定的壽命，通過支持向量機回歸得到剩余壽命與健康指標(biāo)之間的非線性關(guān)系，將其作為卡爾曼濾波觀測方程對原來的狀態(tài)方程進行更新，最后把根據(jù)狀態(tài)方程得到的先驗估計剩余壽命和根據(jù)觀測得到的剩余壽命加權(quán)(卡爾曼增益)，得到運行狀態(tài)下離心泵的剩余壽命。

圖1 基于UKF的離心泵壽命算法流程圖Fig.1 Centrifugal pump life algorithm flowchart based on UKF

2.2 熵值法構(gòu)建復(fù)合健康指標(biāo)

根據(jù)壽命退化模型可知揚程H、功率P以及效率η均與離心泵的壽命有關(guān)，且這3個參數(shù)比較容易獲取。采用多個指標(biāo)更容易獲得離心泵性能退化的真實情況，因此通過熵值法把這3個因素綜合起來構(gòu)建離心泵健康指標(biāo)Hi。熵值法的原理是根據(jù)各個因素的離散程度來確定各因素的重要程度，熵值越大，該因素在所有因素中的權(quán)重就越大。構(gòu)建的復(fù)合健康指標(biāo)如下：

Hi=wHH+wPP+wηη

(11)

式中：w表示權(quán)重，其下標(biāo)表示各個指標(biāo)。

權(quán)重的計算方法為：

(12)

式中：ej為第j個指標(biāo)的信息熵值，具體公式見文獻[20]。

根據(jù)已有同類離心泵相關(guān)歷史數(shù)據(jù)，可以較容易獲得離心泵失效時的性能指標(biāo)，該值即為健康指標(biāo)的閾值。

根據(jù)該閾值和當(dāng)前指標(biāo)可以計算離心泵的剩余壽命。但該復(fù)合健康指標(biāo)只有在符合高斯分布的情況下才能更好地應(yīng)用卡爾曼濾波算法，因此如果構(gòu)造的初始復(fù)合健康指標(biāo)不能近似為高斯分布，則可通過Box-Cox變換將其轉(zhuǎn)化為高斯分布[21]。

2.3 UKF濾波剩余壽命模型方程構(gòu)建與求解

健康指標(biāo)與離心泵運行時間之間的關(guān)系可表示為公式(5)所示的二次函數(shù)。當(dāng)健康指標(biāo)閾值確定后，對應(yīng)的離心泵壽命設(shè)置為0，若已知離心泵的總壽命，則公式(5)可轉(zhuǎn)化為健康指標(biāo)和剩余壽命的函數(shù)關(guān)系。該二次函數(shù)的系數(shù)可通過同類離心泵的全壽命試驗數(shù)據(jù)擬合得到。由于離心泵退化曲線函數(shù)光滑且單調(diào)(不考慮突發(fā)故障)，所以剩余壽命亦可表示為健康指標(biāo)的二次函數(shù)，進而可以得到k時刻與k-1時刻離心泵剩余壽命與健康指標(biāo)Hi之間的關(guān)系式：

Lk=Lk-1+A′·Hi+wk

(13)

公式(13)即為卡爾曼濾波的狀態(tài)方程。式中的A′為常數(shù)，wk為符合正態(tài)分布的噪聲，受不同廠家不同批次產(chǎn)品的影響。

如果要對某一特定的離心泵進行壽命預(yù)測，則除了得到上述由同類離心泵設(shè)備擬合得到的狀態(tài)方程外，還要考慮該設(shè)備實時運行的參數(shù)。因此，還需建立觀測方程來對狀態(tài)方程進行實時更新，以獲得運行中離心泵的剩余壽命。觀測方程為一個非線性函數(shù)，該非線性關(guān)系通過支持向量機對離心泵的歷史運行數(shù)據(jù)學(xué)習(xí)得到，具體過程如下：

對于樣本序列S={(Hi(k),L(k))，Hi(k)∈Rl,L(k)∈R,k=1,2,……,N}，前n個樣本數(shù)據(jù)用于支持向量機訓(xùn)練，后N-n個數(shù)據(jù)樣本用于預(yù)測。構(gòu)造如式(14)和式(15)所示的訓(xùn)練樣本對Xtrain、Ytrain和預(yù)測樣本輸入Xtest。具體到本文，將收集的離心泵剩余壽命和健康指標(biāo)數(shù)據(jù)分別作為輸入X和輸出Y，歸一化和隨機順序后，挑選前n個數(shù)據(jù)代入式(7)，則可求解拉格朗日乘子，得到非線性關(guān)系h(·)，然后將拉格朗日乘子與剩余壽命預(yù)測樣本輸入回歸函數(shù)，再計算得到預(yù)測樣本，輸出健康指標(biāo)Ytest[22-23]。

Xtrain=[L(1),L(1),……,L(n)]

(14)

Ytrain=[Hi(1),Hi(2),……,Hi(n)]

(15)

Xtest=[L(n+1),L(n+2),……,L(N-n)]

(16)

上述訓(xùn)練過程得到的非線性函數(shù)即為卡爾曼濾波算法中的觀測方程，可表示為：

Hi,k-1=h(Lk-1)+vk

(17)

式中：vk為SVM模型殘差，假設(shè)為高斯白噪聲，來源于測量儀器儀表的噪聲。

2.3.1 初始向量(剩余壽命)

剩余壽命計算如下：

(18)

(19)

2.3.2 計算Sigma點

對2n+1個Sigma點進行如下計算：

(20)

2.3.3 狀態(tài)更新

狀態(tài)方程(公式(13))基于Sigma點下為：

χk|k-1=f(χk-1)

(21)

將Sigma點進行加權(quán)求和：

(22)

求Sigma點的加權(quán)預(yù)測方差：

(23)

式中：Qk為系統(tǒng)噪聲向量wk的協(xié)方差矩陣。

進行Sigma點的觀測方程(公式17)預(yù)測值求解。

yk|k-1=h(χk-1)

(24)

(25)

2.3.4 進行觀測更新

協(xié)方差估計進行第k次更新：

(26)

觀測中協(xié)方差進行第k次更新：

(27)

式中：Rk為測量噪聲向量vk的協(xié)方差矩陣。

卡爾曼增益進行第k次預(yù)測更新：

(28)

最優(yōu)狀態(tài)估計進行第k次更新：

(29)

最優(yōu)預(yù)測誤差的估計進行第k次更新：

(30)

3 模型仿真與驗證分析

3.1 同類離心泵退化模型數(shù)據(jù)擬合

用MATLAB實現(xiàn)上述模型算法，采用郭文琪等的離心泵加速退化試驗數(shù)據(jù)對模型進行驗證和分析。離心泵加速退化試驗所用離心泵型號為IS100-80-125，揚程為20 m。選取石英砂作為磨料加入水中，石英砂顆粒尺寸為2～4目，避免顆粒較大引起管路堵塞，從而對離心泵的運行數(shù)據(jù)造成干擾，離心泵每天運行8 h，每2 h采集一次試驗數(shù)據(jù)，試驗時長25 d，共有200 h、100組壽命預(yù)測指標(biāo)時間序列數(shù)據(jù)[7]。采用二次曲線的形式對原始數(shù)據(jù)擬合，擬合曲線如圖2所示，方程如式(31)所示。

y=-0.000 203 8x2+0.001 095x+8.835

(31)

該擬合模型均方根誤差為0.098 5，R-square系數(shù)為0.997 7，擬合效果較好?？梢娪枚魏瘮?shù)建立離心泵剩余壽命和健康指標(biāo)的方法較為可行，同時也驗證了模型中卡爾曼狀態(tài)方程的可行性，擬合噪聲的分布可由壽命指標(biāo)的真實值和擬合曲線之間的誤差得到，并變換為高斯白噪聲。

圖2 壽命健康指標(biāo)與時間的關(guān)系Fig.2 Relationship between life health index and time

3.2 基于SVM方法的觀測方程可行性驗證

上述擬合曲線只能反映類似機理同類型離心泵性能退化趨勢，針對某一特定離心泵，除了考慮共性因素外，還需要考慮其個性因素，即對性能退化模型中的系數(shù)及時更新。對于某一處在退化期的同類但應(yīng)用于不同工況的離心泵，假設(shè)已知其前期的80組退化數(shù)據(jù)，計算健康指標(biāo)后，歸一化處理并打亂數(shù)據(jù)順序，選取前60組數(shù)據(jù)對SVM進行訓(xùn)練，采用后20組數(shù)據(jù)進行預(yù)測。模型的初始參數(shù)如下：誤差初值為 0.000 1，懲罰因子初值為10 000，采用10-fold交叉驗證度量SVM精度及推廣能力。最終得到SVM方法預(yù)測的健康指標(biāo)如圖3所示。健康指標(biāo)的平均絕對誤差為-0.004 7，均方根誤差為0.054 1?？梢钥闯鯯VM方法對同類離心泵的歷史數(shù)據(jù)回歸和預(yù)測比較準(zhǔn)確，可以用于模型中卡爾曼濾波觀測方程的建立。

圖3 壽命健康指標(biāo)實際值與SVM預(yù)測值Fig.3 Actual value vs SVM predicted value of life health index

3.3 無跡卡爾曼濾波仿真結(jié)果

離心泵在衰退期內(nèi)的性能退化過程可以利用上述建立的模型進行跟蹤，跟蹤結(jié)果如圖4所示。退化初期由于可利用的測量數(shù)據(jù)較少，所以跟蹤結(jié)果不理想，與真實值之間的誤差較大；隨著離心泵可利用數(shù)據(jù)的增加，UKF跟蹤結(jié)果與真實值之間的誤差慢慢變小。壽命健康指標(biāo)跟蹤誤差如圖5所示。由圖5可以看出，當(dāng)跟蹤結(jié)果在58 h之后，所建立的方法能夠較為準(zhǔn)確地反映離心泵的性能退化情況。假設(shè)壽命健康指標(biāo)<3時，離心泵失效，則進入衰退期后，離心泵可運行約178 h。上述卡爾曼濾波模型和試驗方法得到的剩余壽命和衰退期內(nèi)運行時間的關(guān)系如圖6所示。由圖6可以看出，當(dāng)模型初值與真實值誤差較大時，隨著運行時間的進行，建立的卡爾曼濾波方法能夠較快收斂到真實剩余壽命附近。但由于目前試驗數(shù)據(jù)不足，所以在未來的工作中需要獲取更多的試驗數(shù)據(jù)來驗證和完善建立模型的穩(wěn)定性和預(yù)測精度。

圖4 UKF方法離心泵衰退期跟蹤結(jié)果Fig.4 Centrifugal pump degradation track by UKF method

圖5 UKF方法壽命健康指標(biāo)跟蹤誤差Fig.5 Life health index tracking error of UKF method

圖6 剩余壽命與衰退期內(nèi)運行時間的關(guān)系Fig.6 Relationship between residual life and running time in degradation period

4 結(jié) 論

根據(jù)理論分析，考慮離心泵的葉輪和軸承磨損，推導(dǎo)了離心泵揚程的衰退模型，并利用支持向量機和無跡卡爾曼濾波方法建立了離心泵剩余壽命預(yù)測模型，得出如下結(jié)論。

(1)考慮葉輪和軸承磨損失效時，推導(dǎo)了離心泵壽命與其揚程、功率及效率的關(guān)系，該關(guān)系可用二次函數(shù)的形式來表示。根據(jù)離心泵加速退化試驗數(shù)據(jù)，并對試驗數(shù)據(jù)進行擬合，擬合結(jié)果比較符合二次函數(shù)規(guī)律。

(2)根據(jù)在線監(jiān)測的部分離心泵性能退化健康指標(biāo)數(shù)據(jù)，支持向量機能夠較為準(zhǔn)確地預(yù)測離心泵性能退化過程，可以為卡爾曼濾波提供較為準(zhǔn)確的觀測值。

(3)MATLAB仿真結(jié)果表明，建立的模型在生命健康指標(biāo)跟蹤初期誤差較大，但當(dāng)跟蹤一段時間后，隨著可用數(shù)據(jù)的增加，跟蹤誤差逐漸減小，可以較為準(zhǔn)確地描述離心泵性能退化過程，并對剩余壽命進行預(yù)測。