王少萍，耿藝璇，石存

1. 北京航空航天大學(xué) 自動(dòng)化科學(xué)與電氣工程學(xué)院，北京 100191 2. 北京航空航天大學(xué) 寧波創(chuàng)新研究院，寧波 315800

柱塞式航空液壓泵是航空液壓系統(tǒng)的關(guān)鍵部件，具有高速、高壓、長(zhǎng)壽命的特征，其可靠性和壽命對(duì)飛機(jī)任務(wù)完成和安全服役至關(guān)重要。航空液壓泵結(jié)構(gòu)復(fù)雜、使用工況惡劣，在交變載荷作用下其關(guān)鍵摩擦副(滑靴副、配流副和柱塞副)會(huì)發(fā)生不同程度的摩擦磨損，進(jìn)而導(dǎo)致泄漏增加，當(dāng)性能難以滿足液壓系統(tǒng)的流量需求時(shí)會(huì)發(fā)生失效[1]。鑒于航空液壓泵要求高可靠長(zhǎng)壽命，出廠前很難通過(guò)試驗(yàn)確定準(zhǔn)確的壽命指標(biāo)；而航空液壓泵外場(chǎng)飛行數(shù)據(jù)雖然非常珍貴，但存在大量不確定性因素，也尚未應(yīng)用于其壽命估計(jì)。因此，目前亟需摸清航空液壓泵的失效物理以建立宏微觀性能退化規(guī)律，并將源源不斷的外場(chǎng)飛行數(shù)據(jù)有機(jī)融合到退化關(guān)系模型中，準(zhǔn)確描述其全生命周期的性能退化關(guān)系，進(jìn)而實(shí)現(xiàn)液壓泵的準(zhǔn)確壽命估計(jì)。

目前壽命估計(jì)方法可分為基于模型驅(qū)動(dòng)的方法、基于數(shù)據(jù)驅(qū)動(dòng)的方法和模型-數(shù)據(jù)混合方法這3類[2]?；谖锢砟Ｐ偷姆椒ㄍㄟ^(guò)選擇與產(chǎn)品壽命高度相關(guān)的物理變量，在分析產(chǎn)品退化失效機(jī)理的前提下建立數(shù)學(xué)模型刻畫(huà)產(chǎn)品的性能演變過(guò)程，通過(guò)對(duì)性能特征量的監(jiān)測(cè)和估計(jì)預(yù)測(cè)產(chǎn)品的剩余壽命。該類方法已成功應(yīng)用于機(jī)械零部件(如潤(rùn)滑系統(tǒng)[3]、軸承[4-5])、半導(dǎo)體器件(如功率MOS器件[6])、機(jī)電部件(如感應(yīng)電動(dòng)機(jī)[7])、光電部件(如激光器[8]、LED[9])和電子器件(如電容[10]、蓄電池[11])等產(chǎn)品的失效建模與壽命預(yù)測(cè)。然而，在失效機(jī)理分析不完備或部分模型參數(shù)很難獲取的情況下，完全基于物理模型的研究方法也會(huì)在模型外推時(shí)產(chǎn)生較大誤差，導(dǎo)致壽命預(yù)測(cè)精度下降?；跀?shù)據(jù)驅(qū)動(dòng)的壽命估計(jì)方法通過(guò)搜集大量樣本對(duì)數(shù)據(jù)模型進(jìn)行訓(xùn)練，進(jìn)而得到可以描述產(chǎn)品性能退化規(guī)律的模型，如神經(jīng)網(wǎng)絡(luò)模型[12-13]或隨機(jī)過(guò)程模型等[14-20]，但該類方法無(wú)法進(jìn)行物理解釋，且在樣本量較少時(shí)很難使用。

為了對(duì)樣本和監(jiān)測(cè)數(shù)據(jù)有限的產(chǎn)品進(jìn)行準(zhǔn)確的壽命估計(jì)，研究人員發(fā)現(xiàn)必須在深入挖掘其失效物理機(jī)理的同時(shí)，盡可能充分利用設(shè)備的運(yùn)維數(shù)據(jù)，克服單一壽命預(yù)測(cè)方法對(duì)模型或數(shù)據(jù)的強(qiáng)依賴性，于是提出了一系列模型-數(shù)據(jù)混合方法，并對(duì)其可行性進(jìn)行了驗(yàn)證。孫宏達(dá)等[21]提出了一種失效物理-數(shù)據(jù)融合方法，建立了燃油泵出口壓力退化模型，使用無(wú)跡卡爾曼濾波算法實(shí)現(xiàn)模型參數(shù)更新和壽命估計(jì)，并基于機(jī)載燃油泵性能退化數(shù)據(jù)集驗(yàn)證了該方法的有效性。范立明等[22]建立了鋰電池退化經(jīng)驗(yàn)?zāi)Ｐ?，通過(guò)對(duì)容量等測(cè)量數(shù)據(jù)進(jìn)行濾波動(dòng)態(tài)更新模型參數(shù)，實(shí)現(xiàn)了對(duì)鋰電池剩余壽命的準(zhǔn)確預(yù)測(cè)。馬里蘭大學(xué)的Chen和Pecht[23]考慮了壽命預(yù)測(cè)過(guò)程中的不確定性，采用粒子濾波算法對(duì)鋰電池退化模型參數(shù)進(jìn)行更新，并基于試驗(yàn)數(shù)據(jù)驗(yàn)證了該方法的有效性。Zhao等[24]將有限元模型與狀態(tài)監(jiān)測(cè)數(shù)據(jù)融合，利用支持向量機(jī)(SVM)預(yù)測(cè)模型更新由有限元分析得到的退化模型參數(shù)，實(shí)現(xiàn)了對(duì)齒輪裂紋擴(kuò)展長(zhǎng)度和剩余壽命的估計(jì)。Zang等[25]將基于Paris裂紋擴(kuò)展模型的粒子濾波算法與前饋神經(jīng)網(wǎng)絡(luò)融合，提出了一種適用于軌道電纜的剩余使用壽命預(yù)測(cè)方法。總體來(lái)說(shuō)，現(xiàn)階段模型-數(shù)據(jù)融合壽命估計(jì)方法尚未真正將失效物理與數(shù)據(jù)驅(qū)動(dòng)有機(jī)結(jié)合，相關(guān)理論與方法仍有待進(jìn)一步探索和驗(yàn)證。

基于上述原因，本文提出一種基于失效物理模型和數(shù)據(jù)(內(nèi)場(chǎng)研制、外場(chǎng)服役、維修保障等)調(diào)制融合的航空液壓泵壽命估計(jì)方法，構(gòu)建了混合潤(rùn)滑條件下多場(chǎng)耦合的液壓泵失效機(jī)理模型，并將其性能退化用隨機(jī)過(guò)程表征，然后將外場(chǎng)和維修數(shù)據(jù)作為實(shí)測(cè)數(shù)據(jù)通過(guò)粒子濾波不斷更新調(diào)制到性能退化過(guò)程中，期間考慮不同信息的融合，實(shí)現(xiàn)航空液壓泵的準(zhǔn)確退化狀態(tài)估計(jì)和全生命周期剩余壽命預(yù)測(cè)。

1 混合潤(rùn)滑下多場(chǎng)耦合失效機(jī)理建模

航空液壓泵存在3對(duì)關(guān)鍵摩擦副，即滑靴副、配流副和柱塞副。這里以滑靴副為例介紹其混合潤(rùn)滑條件下的多場(chǎng)耦合失效機(jī)理建模思路，其他兩對(duì)摩擦副的性能退化建模與之類似。

圖1所示為液壓泵滑靴副剖面圖，當(dāng)航空液壓泵高速運(yùn)轉(zhuǎn)時(shí)，滑靴在慣性力、離心力矩、摩擦力作用下可能發(fā)生傾覆，局部潤(rùn)滑油膜厚度變小，滑靴和斜盤不再處于平行潤(rùn)滑狀態(tài)，使得滑靴/斜盤間形成楔形油膜，引發(fā)滑靴偏磨。

考慮到滑靴表面相對(duì)粗糙、材料硬度低，斜盤表面相對(duì)光滑、材料硬度高，建模中假設(shè)斜盤表面光滑，滑靴表面具有一定的粗糙度，因此滑靴副具有不規(guī)則的油膜分布。為了推導(dǎo)滑靴副的偏磨磨損模型，對(duì)滑靴/斜盤間的楔形油膜進(jìn)行網(wǎng)格劃分(如圖2所示)。在每個(gè)網(wǎng)格區(qū)間里，可近似認(rèn)為油膜均勻分布，平均油膜的厚度為

h(r,θ+dθ)+h(r+dr,θ+dθ))

(1)

式中：h(r,θ)為任意區(qū)間頂點(diǎn)的油膜厚度，其表達(dá)式為

(2)

其中:routS和rinS分別是滑靴底面密封帶的外徑和內(nèi)徑;ΔhT和ΔhP、ΔhC分別是考慮金屬表面熱變形、底面流體壓力引起的彈性變形、滑靴磨損輪廓對(duì)油膜厚度影響的修正量。

基于式(1)、式(2)和文獻(xiàn)[26]提出的磨粒磨損模型，可得到滑靴底面(r,θ)區(qū)域內(nèi)單位行程的磨損量：

(3)

式中:Rw(r,θ)為(r,θ)區(qū)域內(nèi)滑靴相對(duì)斜盤表面單位行程的磨損量；v為滑靴底面相對(duì)斜盤表面的相對(duì)運(yùn)動(dòng)速度；δ0為接觸表面粗糙峰的高度；Δ為接觸面粗糙峰的長(zhǎng)度；τ=ηm/Gc為延遲時(shí)間，其中ηm為材料的動(dòng)力黏度，Gc為材料的剪切模量。將式(3)在整個(gè)求解域內(nèi)進(jìn)行積分，可以得到整個(gè)滑靴底面的磨損量：

(4)

假設(shè)航空液壓泵內(nèi)共有N個(gè)滑靴組件，而缸體每旋轉(zhuǎn)一周各滑靴的行程為

L=2πR+4R(secβ-1)

(5)

式中：R為柱塞孔分度圓半徑；β為斜盤傾角。將式(4)沿滑靴運(yùn)動(dòng)軌跡進(jìn)行曲線積分，可得到缸體每旋轉(zhuǎn)一周泵內(nèi)所有滑靴的總磨損量：

(6)

航空液壓泵滑靴在單位旋轉(zhuǎn)時(shí)間(s)內(nèi)的磨損體積為

(7)

基于式(7)并對(duì)動(dòng)態(tài)潤(rùn)滑偏磨磨損模型進(jìn)行適當(dāng)簡(jiǎn)化，可推導(dǎo)得到單位時(shí)間內(nèi)的滑靴磨損高度和單位時(shí)間內(nèi)的滑靴副泄漏間隙增量：

(8)

(9)

對(duì)于航空液壓泵，以殼體回油流量作為系統(tǒng)的觀測(cè)量，基于得到的滑靴副磨損間隙寬度可得到滑靴副流向殼體的泄漏模型：

(pr-pL)=10 000·θ10·kq·

(10)

綜上所述，將滑靴副間隙寬度定義為航空液壓泵的退化狀態(tài)，即x(t)=δ(t)，觀測(cè)量選為殼體回油流量y(t)=Qcase(t)，系統(tǒng)輸入u(t)=[u1,u2]=[n,ps]，則系統(tǒng)性能退化狀態(tài)方程為

(11)

對(duì)動(dòng)態(tài)方程(11)進(jìn)行適當(dāng)簡(jiǎn)化，得到離散系統(tǒng)動(dòng)態(tài)方程：

(12)

方程(12)從失效物理角度考慮了航空液壓泵滑靴副退化狀態(tài)的動(dòng)態(tài)演變規(guī)律，但考慮到外場(chǎng)使用中，系統(tǒng)工況和外部載荷通常難以準(zhǔn)確獲取，系統(tǒng)退化過(guò)程存在固有不確定性，且退化狀態(tài)和觀測(cè)量會(huì)受到過(guò)程噪聲和觀測(cè)噪聲的影響，還需對(duì)航空液壓泵性能退化模型的狀態(tài)轉(zhuǎn)移方程和觀測(cè)方程進(jìn)一步添加噪聲，得到：

(13)

2 基于粒子濾波的數(shù)據(jù)更新算法

貝葉斯濾波估計(jì)為產(chǎn)品性能退化狀態(tài)估計(jì)和剩余壽命預(yù)測(cè)提供了一套完整的理論體系框架，目前已被成功應(yīng)用于液壓領(lǐng)域[27]。在貝葉斯濾波方法中，系統(tǒng)當(dāng)前時(shí)刻的退化狀態(tài)可由系統(tǒng)上一時(shí)刻的退化狀態(tài)遞推得到，然后利用當(dāng)前時(shí)刻的觀測(cè)數(shù)據(jù)完成貝葉斯動(dòng)態(tài)更新，得到狀態(tài)的最優(yōu)后驗(yàn)估計(jì)。基于動(dòng)態(tài)系統(tǒng)退化狀態(tài)的不斷遞推估計(jì)，就可以實(shí)現(xiàn)長(zhǎng)周期內(nèi)的系統(tǒng)剩余壽命預(yù)測(cè)。

目前常用的非線性系統(tǒng)貝葉斯濾波方法主要包括擴(kuò)展卡爾曼濾波算法、無(wú)跡卡爾曼濾波算法和粒子濾波算法3種[28-30]。其中，擴(kuò)展卡爾曼濾波的基本思想是利用一階Taylor展開(kāi)，首先將非線性系統(tǒng)在狀態(tài)估計(jì)值處進(jìn)行線性化，在此基礎(chǔ)上再使用卡爾曼濾波公式對(duì)狀態(tài)變量進(jìn)行估計(jì)。無(wú)跡卡爾曼濾波采用統(tǒng)計(jì)線性化的方法近似得到非線性模型中的狀態(tài)變量概率分布，與擴(kuò)展卡爾曼濾波方法相比不需要對(duì)雅可比矩陣進(jìn)行求導(dǎo)和忽略高階項(xiàng)，因此可以達(dá)到三階Taylor展開(kāi)的近似精度，但同樣只能使用高斯分布對(duì)真實(shí)的后驗(yàn)分布進(jìn)行近似。

與以上兩種非線性濾波算法相比，粒子濾波(Particle Filter, PF)在貝葉斯濾波方法的基礎(chǔ)上結(jié)合了蒙特卡洛模擬思想，其基本原理是采用一系列隨機(jī)樣本粒子，基于退化狀態(tài)空間模型和系統(tǒng)不斷更新的觀測(cè)數(shù)據(jù)，對(duì)樣本粒子值和粒子權(quán)重值進(jìn)行遞推更新，從而實(shí)現(xiàn)系統(tǒng)直到當(dāng)前時(shí)刻的退化狀態(tài)最優(yōu)后驗(yàn)估計(jì)。該方法不需要對(duì)非線性系統(tǒng)的動(dòng)態(tài)方程進(jìn)行簡(jiǎn)化，對(duì)后驗(yàn)概率分布的特定形式?jīng)]有要求，可以解決非高斯噪聲下的狀態(tài)估計(jì)問(wèn)題，具備較好的處理退化過(guò)程不確定性的能力。經(jīng)調(diào)研發(fā)現(xiàn)，工程中的實(shí)際航空液壓泵性能退化系統(tǒng)是非線性且非高斯的，因此，粒子濾波方法在解決航空液壓泵的動(dòng)態(tài)退化狀態(tài)估計(jì)與壽命預(yù)測(cè)方面具有一定優(yōu)勢(shì)。

2.1 通用粒子濾波算法

粒子濾波器又稱為序貫蒙特卡羅(Sequential Monte Carlo, SMC)方法，可以基于蒙特卡羅仿真近似逼近貝葉斯濾波中非線性函數(shù)的積分求解，從而得到退化狀態(tài)估計(jì)的近似最優(yōu)解。

對(duì)任意一個(gè)動(dòng)態(tài)系統(tǒng)，可定義如下?tīng)顟B(tài)方程：

(14)

式中：fk(·)為該動(dòng)態(tài)系統(tǒng)的狀態(tài)轉(zhuǎn)移方程；hk(·)為系統(tǒng)的測(cè)量方程；fk(·)和hk(·)都可能是非線性函數(shù)；xk和xk-1分別表示該系統(tǒng)在k時(shí)刻和k-1時(shí)刻的狀態(tài)量；uk-1和uk分別為系統(tǒng)在k-1和k時(shí)刻的輸入；ωk-1為系統(tǒng)k-1時(shí)刻的隨機(jī)過(guò)程噪聲;zk為k時(shí)刻的系統(tǒng)狀態(tài)觀測(cè)量；υk為系統(tǒng)在k時(shí)刻的隨機(jī)測(cè)量噪聲。在粒子濾波器中，過(guò)程噪聲ωk-1和測(cè)量噪聲υk都可以是非Gaussian噪聲。

基于上述狀態(tài)方程(14)，粒子濾波器的退化狀態(tài)估計(jì)過(guò)程可以分為預(yù)測(cè)和更新兩個(gè)階段。假設(shè)當(dāng)前時(shí)刻為k，在預(yù)測(cè)階段，粒子濾波器利用系統(tǒng)1:k-1時(shí)刻的所有觀測(cè)值z(mì)1:k-1，基于系統(tǒng)的退化物理模型f(·)計(jì)算當(dāng)前時(shí)刻退化狀態(tài)值xk的先驗(yàn)概率p(xk|z1:k-1)；在更新階段，粒子濾波器利用k時(shí)刻的最新觀測(cè)值z(mì)k對(duì)先驗(yàn)概率進(jìn)行更新，得到當(dāng)前時(shí)刻的退化狀態(tài)后驗(yàn)概率p(xk|z1:k)。

蒙特卡羅模擬的基本思想是用樣本均值代替期望值，假設(shè)可以從后驗(yàn)概率中采樣到N個(gè)樣本，則基于蒙特卡羅近似得到的后驗(yàn)概率為

(15)

(16)

然而，由于狀態(tài)量的后驗(yàn)概率未知，通常不能從后驗(yàn)概率分布中直接進(jìn)行采樣，因此需要選擇一個(gè)建議分布q(x|z)以替代上述p(x|z)。假設(shè)q(x|z)的支撐集涵蓋了p(x|z)的支撐集，將其代入式(16)可得：

(17)

式中：Wk(xk)為粒子權(quán)重，計(jì)算公式為

(18)

將式(17)中的分母改寫為

(19)

整理可得：

E[f(xk)]=

(20)

(21)

q(x0:k|z1:k)=q(xk|x0:k-1,z1:k)q(x0:k-1|z1:k-1)

(22)

則根據(jù)貝葉斯定理可以得到粒子k時(shí)刻重要性權(quán)值和k-1時(shí)刻重要性權(quán)值的遞歸形式為

(23)

進(jìn)一步假設(shè)動(dòng)態(tài)系統(tǒng)k時(shí)刻狀態(tài)xk的重要性概率僅由k-1時(shí)刻的狀態(tài)xk-1和當(dāng)前時(shí)刻的觀測(cè)值z(mì)k決定，則公式可以進(jìn)一步簡(jiǎn)化為

(24)

然而在實(shí)際應(yīng)用中，上述通用粒子濾波框架常面臨重要性概率密度函數(shù)的選擇和樣本粒子枯竭問(wèn)題，因此需要對(duì)通用粒子濾波算法進(jìn)行改進(jìn)。

2.2 基于重要性重采樣的粒子濾波算法

通用粒子濾波的重要性密度函數(shù)通常為一般采樣退化狀態(tài)的轉(zhuǎn)移概率密度函數(shù)p(xk|xk-1)，其優(yōu)點(diǎn)是計(jì)算簡(jiǎn)單且容易實(shí)現(xiàn)，但由于在粒子采樣中沒(méi)有考慮系統(tǒng)最新的測(cè)量數(shù)據(jù)，得到的粒子群往往與真實(shí)后驗(yàn)分布存在一定程度的偏差。特別是在粒子的似然函數(shù)分布形狀比較狹窄或者似然函數(shù)位于粒子先驗(yàn)分布的尾部時(shí)，該方法的計(jì)算誤差將會(huì)明顯增加。

重要性概率密度函數(shù)選擇的核心原則是要使粒子的權(quán)值方差盡可能小，基于該原則可以得到最優(yōu)重要性概率密度函數(shù)：

(25)

而根據(jù)貝葉斯公式：

(26)

將式(26)代入式(25)，可得最優(yōu)重要性概率密度函數(shù)：

(27)

將式(27)代入式(24)，可確定基于最優(yōu)重要性概率密度函數(shù)更新后的粒子權(quán)重為

(28)

使用式(27)的最優(yōu)重要性概率密度函數(shù)進(jìn)行重采樣，則新粒子的產(chǎn)生既考慮了已有的先驗(yàn)粒子信息，又考慮了當(dāng)前時(shí)刻的最新觀測(cè)數(shù)據(jù)。相當(dāng)于將先驗(yàn)分布中權(quán)值較大的粒子移動(dòng)到似然函數(shù)值較大區(qū)域，使得產(chǎn)生的新粒子群中有效粒子數(shù)變多、所需粒子總數(shù)減小，從而降低計(jì)算成本。

另一方面，通用粒子濾波算法雖然可以通過(guò)重采樣在一定程度上克服粒子權(quán)值退化，但同時(shí)也帶來(lái)了粒子多樣性枯竭的問(wèn)題，其原理及改進(jìn)方法如圖3所示?？梢钥吹?，隨著迭代計(jì)算次數(shù)增加，只有少數(shù)粒子權(quán)重較大而其余粒子權(quán)重很小，重采樣后獲得的新粒子集基本復(fù)制了少數(shù)的大權(quán)重粒子，導(dǎo)致大部分粒子所包含的統(tǒng)計(jì)信息丟失，樣本多樣性降低，狀態(tài)估計(jì)精度下降。

為解決上述粒子枯竭問(wèn)題，根據(jù)離散的后驗(yàn)概率密度構(gòu)建其連續(xù)分布形式，從重構(gòu)的連續(xù)分布中進(jìn)行采樣，其采樣過(guò)程可近似為

(29)

(30)

Kh(x)表示對(duì)核密度K(·)進(jìn)行重新標(biāo)定后的核密度函數(shù)，h(h>0)表示核帶寬；nx為狀態(tài)向量維數(shù)。重新標(biāo)定的原則是通過(guò)選取核帶寬h和核函數(shù)K(·)，使以下正則經(jīng)驗(yàn)密度和真實(shí)后驗(yàn)密度間的均方誤差期望值最小：

(31)

選擇Epanechnikov核作為最佳核密度函數(shù)，其表達(dá)式為

(32)

式中：cnx表示半徑為Rnx的單位球的體積。假設(shè)粒子濾波所估計(jì)狀態(tài)的先驗(yàn)密度是一個(gè)協(xié)方差矩陣為單位陣的高斯分布，則對(duì)應(yīng)的最佳核帶寬為

(33)

(34)

綜上所述，采用最優(yōu)重要性概率密度函數(shù)和正則化變換重采樣改進(jìn)粒子濾波算法步驟如下：

步驟1粒子初始化：初始0時(shí)刻，根據(jù)系統(tǒng)先驗(yàn)分布生成初始粒子群，粒子數(shù)量為Ns。

(35)

步驟2粒子值更新：當(dāng)前k(k≥1)時(shí)刻，根據(jù)式(13)建立的失效物理退化模型更新粒子值及權(quán)重：

i=1,2,…,Ns

(36)

步驟3粒子權(quán)重更新：當(dāng)前k時(shí)刻，更新粒子權(quán)重并進(jìn)行歸一化處理：

(37)

可得k時(shí)刻系統(tǒng)狀態(tài)xk的最優(yōu)估計(jì)：

(38)

(39)

上述粒子濾波算法引入最新觀測(cè)數(shù)據(jù)對(duì)粒子權(quán)值進(jìn)行更新調(diào)整，從基于離散后驗(yàn)概率密度重構(gòu)的近似連續(xù)分布中對(duì)粒子進(jìn)行重采樣，得到的粒子后驗(yàn)分布雖然較為分散，但新粒子群中有效粒子數(shù)變多、所需粒子總數(shù)減小，降低了計(jì)算成本，而且能在一定程度上解決粒子多樣性的快速枯竭問(wèn)題，避免算法對(duì)少數(shù)奇異點(diǎn)過(guò)于敏感，從而保證狀態(tài)估計(jì)和壽命預(yù)測(cè)的效果。

3 粒子濾波數(shù)據(jù)調(diào)制液壓泵壽命估計(jì)

上述粒子濾波數(shù)據(jù)更新算法考慮了重要性概率密度函數(shù)的選取，同時(shí)也解決了重采樣帶來(lái)的樣本枯竭問(wèn)題，可以提高動(dòng)態(tài)系統(tǒng)的狀態(tài)估計(jì)精度。本部分將上述算法應(yīng)用于航空液壓泵，以求實(shí)現(xiàn)精確的航空液壓泵狀態(tài)估計(jì)和壽命預(yù)測(cè)。

航空液壓泵的整個(gè)壽命周期可分為3個(gè)階段：健康狀態(tài)、故障退化狀態(tài)和嚴(yán)重故障狀態(tài)，其剩余壽命的預(yù)測(cè)一般從故障退化狀態(tài)下開(kāi)始進(jìn)行，通過(guò)壽命預(yù)測(cè)模型得到系統(tǒng)從當(dāng)前時(shí)刻發(fā)展到嚴(yán)重故障狀態(tài)所需要的時(shí)間。根據(jù)動(dòng)態(tài)系統(tǒng)性能退化模型和不斷更新的觀測(cè)數(shù)據(jù)，利用改進(jìn)的粒子濾波方法不斷向后預(yù)測(cè)系統(tǒng)退化狀態(tài)，直到某一時(shí)刻的預(yù)測(cè)狀態(tài)超過(guò)了設(shè)定的失效閾值xth，認(rèn)為液壓泵出現(xiàn)嚴(yán)重故障狀態(tài)，不能達(dá)到正常工作所需標(biāo)準(zhǔn)，航空液壓泵發(fā)生失效。因此，航空液壓泵從當(dāng)前時(shí)刻k到嚴(yán)重故障狀態(tài)的剩余使用壽命(Remaining Useful Life, RUL)計(jì)算公式為

RULk=tk+r-tk

(40)

式中：tk+r為粒子濾波器預(yù)測(cè)的退化狀態(tài)超過(guò)失效閾值的時(shí)刻點(diǎn);tk為系統(tǒng)當(dāng)前時(shí)刻點(diǎn)，r為從當(dāng)前時(shí)刻到預(yù)測(cè)的退化狀態(tài)首次超過(guò)失效閾值的預(yù)測(cè)步數(shù)。

剩余壽命預(yù)測(cè)算法的流程如圖4所示。其具體步驟如下：

步驟1根據(jù)航空液壓泵工作要求設(shè)定動(dòng)態(tài)系統(tǒng)的失效閾值xth，若累積退化量超過(guò)該閾值，認(rèn)為航空液壓泵發(fā)生失效。

步驟2基于所建立的航空液壓泵物理退化模型和當(dāng)前階段觀測(cè)數(shù)據(jù)，確定當(dāng)前階段k時(shí)刻的退化狀態(tài)xk，作為該時(shí)刻剩余壽命預(yù)測(cè)起始點(diǎn)。

步驟4在k=k+1時(shí)刻，引入最新外場(chǎng)觀測(cè)數(shù)據(jù)，對(duì)航空液壓泵的失效物理退化模型進(jìn)行調(diào)制，計(jì)算該時(shí)刻系統(tǒng)的退化狀態(tài)。

步驟5重復(fù)步驟2～步驟4，直到當(dāng)前時(shí)刻基于外場(chǎng)數(shù)據(jù)調(diào)制更新得到的液壓泵性能退化狀態(tài)超過(guò)設(shè)定的失效閾值(即xk>xth)，停止壽命預(yù)測(cè)。

步驟6根據(jù)各時(shí)刻k對(duì)應(yīng)的退化狀態(tài)估計(jì)結(jié)果和剩余使用壽命預(yù)測(cè)值RULk，繪制退化狀態(tài)估計(jì)曲線和剩余使用壽命曲線。

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

為了驗(yàn)證本文提出的基于失效物理退化模型和粒子濾波全息調(diào)制的剩余壽命估計(jì)算法的有效性，選取4臺(tái)某型號(hào)航空液壓泵進(jìn)行了壽命試驗(yàn)。為了加快磨損退化速度采用全流量工況，監(jiān)測(cè)液壓泵回油流量，每隔5 h觀測(cè)一次，得到的退化狀態(tài)估計(jì)結(jié)果如圖5所示。根據(jù)該型號(hào)泵的設(shè)計(jì)參數(shù)，在全流量工況下，回油流量超過(guò)2.8 L/min時(shí)認(rèn)為該泵完全失效，即xth=2.8 L/min。

從圖5中可以看到，航空液壓泵磨損退化過(guò)程具有較大的隨機(jī)性，這主要是因?yàn)橐簤罕迷诠ぷ鬟^(guò)程中關(guān)鍵摩擦副的磨損率隨載荷工況實(shí)時(shí)變化，因此即便同一型號(hào)的液壓泵其退化過(guò)程也不完全相同。對(duì)比試驗(yàn)中的實(shí)際觀測(cè)數(shù)據(jù)與本文方法可以看到，本文所提出的數(shù)據(jù)調(diào)制方法可以實(shí)現(xiàn)長(zhǎng)周期內(nèi)航空液壓泵性能退化狀態(tài)的準(zhǔn)確估計(jì)。

試驗(yàn)中，某一航空液壓泵在連續(xù)監(jiān)測(cè)750 h后的剩余壽命預(yù)測(cè)結(jié)果如圖6所示，其中實(shí)際回油流量觀測(cè)數(shù)據(jù)截止到750 h。圖中黑色曲線表示該液壓泵從開(kāi)始服役到當(dāng)前時(shí)刻的歷史觀測(cè)數(shù)據(jù)，紅色曲線表示基于本文所提出數(shù)據(jù)調(diào)制方法得到的液壓泵性能退化狀態(tài)估計(jì)值及對(duì)未來(lái)退化狀態(tài)的預(yù)測(cè)值，虛線為該液壓泵的失效閾值。從圖6中可以看到，根據(jù)退化狀態(tài)預(yù)測(cè)結(jié)果和設(shè)定的失效閾值，當(dāng)前時(shí)刻(750 h)該液壓泵的剩余使用壽命預(yù)測(cè)值為279 h。

以當(dāng)前液壓泵運(yùn)行時(shí)間(750 h)作為起始點(diǎn)，隨運(yùn)行時(shí)間繼續(xù)增加可以迭代計(jì)算出未來(lái)任意時(shí)刻的剩余使用壽命，其結(jié)果如圖7所示。

圖7中，黑色實(shí)線代表該航空液壓泵的真實(shí)剩余壽命，藍(lán)色折線表示基于本文所提出的數(shù)據(jù)調(diào)制方法得到的剩余壽命預(yù)測(cè)值，可以看到二者高度吻合，進(jìn)一步驗(yàn)證了本文所提出壽命預(yù)測(cè)方法的有效性和準(zhǔn)確性。此外，使用蒙特卡洛方法計(jì)算后驗(yàn)概率通常面臨計(jì)算量較大、計(jì)算時(shí)間相對(duì)較長(zhǎng)的問(wèn)題，為了加快所用航空液壓泵的磨損退化程度，本試驗(yàn)采用了全流量工況，每隔5 h記錄一次泵的回油流量數(shù)據(jù)。經(jīng)過(guò)驗(yàn)證，試驗(yàn)采樣周期5 h遠(yuǎn)遠(yuǎn)大于使用粒子濾波算法進(jìn)行單步壽命預(yù)測(cè)的平均時(shí)間，因此使用所提出方法對(duì)該型號(hào)的航空液壓泵進(jìn)行在線壽命預(yù)測(cè)同樣可行。

最后，考慮到航空液壓泵在全壽命周期中退化過(guò)程存在多階段特性，不同階段適合采用的失效模型可能存在差別，后續(xù)將開(kāi)展多階段失效物理建模，探索有效的基于數(shù)據(jù)的退化階段突變點(diǎn)檢測(cè)算法，將本文方法推廣到多階段退化下的壽命預(yù)測(cè)；同時(shí)，進(jìn)一步收集其他型號(hào)航空液壓泵的在線退化數(shù)據(jù)，對(duì)該方法的在線壽命預(yù)測(cè)效果開(kāi)展進(jìn)一步驗(yàn)證。

5 結(jié) 論

針對(duì)航空液壓泵失效機(jī)理分析不完備、運(yùn)行中載荷工況動(dòng)態(tài)變化且外場(chǎng)觀測(cè)數(shù)據(jù)具有不確定性，本文提出了一種全息數(shù)據(jù)調(diào)制更新的航空液壓泵壽命估計(jì)方法。通過(guò)實(shí)際航空液壓泵試驗(yàn)驗(yàn)證，結(jié)論如下：

1) 基于失效物理與動(dòng)態(tài)數(shù)據(jù)調(diào)制更新的壽命估計(jì)理論在挖掘航空液壓泵退化失效機(jī)理的同時(shí)充分利用觀測(cè)信息，通過(guò)粒子濾波將外場(chǎng)信息有機(jī)調(diào)制到性能退化模型中，實(shí)現(xiàn)了航空液壓泵全生命周期準(zhǔn)確的壽命估計(jì)。

2) 針對(duì)粒子濾波采樣樣本枯竭問(wèn)題，基于最優(yōu)重要性概率密度函數(shù)遞推和正則化改進(jìn)粒子濾波算法，提高了粒子群中的有效粒子比例，增強(qiáng)了算法對(duì)觀測(cè)不確定性的處理能力。

3) 在狀態(tài)估計(jì)和RUL預(yù)測(cè)階段，利用觀測(cè)數(shù)據(jù)對(duì)所建立的失效物理模型進(jìn)行調(diào)制，能較好地刻畫(huà)航空液壓泵的退化隨機(jī)性，RUL預(yù)測(cè)結(jié)果與實(shí)際觀測(cè)結(jié)果高度吻合，可為航空液壓泵維修決策等健康管理活動(dòng)提供支持。