失效物理與數據調制融合的航空液壓泵壽命估計

王少萍，耿藝璇，石存

1. 北京航空航天大學 自動化科學與電氣工程學院，北京 100191 2. 北京航空航天大學 寧波創新研究院，寧波 315800

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

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

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

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

1 混合潤滑下多場耦合失效機理建模

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

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

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

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

(1)

式中：h(r,θ)為任意區間頂點的油膜厚度，其表達式為

(2)

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

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

(3)

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

(4)

假設航空液壓泵內共有N個滑靴組件，而缸體每旋轉一周各滑靴的行程為

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

(5)

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

(6)

航空液壓泵滑靴在單位旋轉時間(s)內的磨損體積為

(7)

基于式(7)并對動態潤滑偏磨磨損模型進行適當簡化，可推導得到單位時間內的滑靴磨損高度和單位時間內的滑靴副泄漏間隙增量：

(8)

(9)

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

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

(10)

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

(11)

對動態方程(11)進行適當簡化，得到離散系統動態方程：

(12)

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

(13)

2 基于粒子濾波的數據更新算法

貝葉斯濾波估計為產品性能退化狀態估計和剩余壽命預測提供了一套完整的理論體系框架，目前已被成功應用于液壓領域[27]。在貝葉斯濾波方法中，系統當前時刻的退化狀態可由系統上一時刻的退化狀態遞推得到，然后利用當前時刻的觀測數據完成貝葉斯動態更新，得到狀態的最優后驗估計?；趧討B系統退化狀態的不斷遞推估計，就可以實現長周期內的系統剩余壽命預測。

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

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

2.1 通用粒子濾波算法

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

對任意一個動態系統，可定義如下狀態方程：

(14)

式中：fk(·)為該動態系統的狀態轉移方程；hk(·)為系統的測量方程；fk(·)和hk(·)都可能是非線性函數；xk和xk-1分別表示該系統在k時刻和k-1時刻的狀態量；uk-1和uk分別為系統在k-1和k時刻的輸入；ωk-1為系統k-1時刻的隨機過程噪聲;zk為k時刻的系統狀態觀測量；υk為系統在k時刻的隨機測量噪聲。在粒子濾波器中，過程噪聲ωk-1和測量噪聲υk都可以是非Gaussian噪聲。

基于上述狀態方程(14)，粒子濾波器的退化狀態估計過程可以分為預測和更新兩個階段。假設當前時刻為k，在預測階段，粒子濾波器利用系統1:k-1時刻的所有觀測值z1:k-1，基于系統的退化物理模型f(·)計算當前時刻退化狀態值xk的先驗概率p(xk|z1:k-1)；在更新階段，粒子濾波器利用k時刻的最新觀測值zk對先驗概率進行更新，得到當前時刻的退化狀態后驗概率p(xk|z1:k)。

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

(15)

(16)

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

(17)

式中：Wk(xk)為粒子權重，計算公式為

(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)

則根據貝葉斯定理可以得到粒子k時刻重要性權值和k-1時刻重要性權值的遞歸形式為

(23)

進一步假設動態系統k時刻狀態xk的重要性概率僅由k-1時刻的狀態xk-1和當前時刻的觀測值zk決定，則公式可以進一步簡化為

(24)

然而在實際應用中，上述通用粒子濾波框架常面臨重要性概率密度函數的選擇和樣本粒子枯竭問題，因此需要對通用粒子濾波算法進行改進。

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

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

重要性概率密度函數選擇的核心原則是要使粒子的權值方差盡可能小，基于該原則可以得到最優重要性概率密度函數：

(25)

而根據貝葉斯公式：

(26)

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

(27)

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

(28)

使用式(27)的最優重要性概率密度函數進行重采樣，則新粒子的產生既考慮了已有的先驗粒子信息，又考慮了當前時刻的最新觀測數據。相當于將先驗分布中權值較大的粒子移動到似然函數值較大區域，使得產生的新粒子群中有效粒子數變多、所需粒子總數減小，從而降低計算成本。

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

為解決上述粒子枯竭問題，根據離散的后驗概率密度構建其連續分布形式，從重構的連續分布中進行采樣，其采樣過程可近似為

(29)

(30)

Kh(x)表示對核密度K(·)進行重新標定后的核密度函數，h(h>0)表示核帶寬；nx為狀態向量維數。重新標定的原則是通過選取核帶寬h和核函數K(·)，使以下正則經驗密度和真實后驗密度間的均方誤差期望值最?。?/p>

(31)

選擇Epanechnikov核作為最佳核密度函數，其表達式為

(32)

式中：cnx表示半徑為Rnx的單位球的體積。假設粒子濾波所估計狀態的先驗密度是一個協方差矩陣為單位陣的高斯分布，則對應的最佳核帶寬為

(33)

(34)

綜上所述，采用最優重要性概率密度函數和正則化變換重采樣改進粒子濾波算法步驟如下：

步驟1粒子初始化：初始0時刻，根據系統先驗分布生成初始粒子群，粒子數量為Ns。

(35)

步驟2粒子值更新：當前k(k≥1)時刻，根據式(13)建立的失效物理退化模型更新粒子值及權重：

i=1,2,…,Ns

(36)

步驟3粒子權重更新：當前k時刻，更新粒子權重并進行歸一化處理：

(37)

可得k時刻系統狀態xk的最優估計：

(38)

(39)

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

3 粒子濾波數據調制液壓泵壽命估計

上述粒子濾波數據更新算法考慮了重要性概率密度函數的選取，同時也解決了重采樣帶來的樣本枯竭問題，可以提高動態系統的狀態估計精度。本部分將上述算法應用于航空液壓泵，以求實現精確的航空液壓泵狀態估計和壽命預測。

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

RULk=tk+r-tk

(40)

式中：tk+r為粒子濾波器預測的退化狀態超過失效閾值的時刻點;tk為系統當前時刻點，r為從當前時刻到預測的退化狀態首次超過失效閾值的預測步數。

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

步驟1根據航空液壓泵工作要求設定動態系統的失效閾值xth，若累積退化量超過該閾值，認為航空液壓泵發生失效。

步驟2基于所建立的航空液壓泵物理退化模型和當前階段觀測數據，確定當前階段k時刻的退化狀態xk，作為該時刻剩余壽命預測起始點。

步驟4在k=k+1時刻，引入最新外場觀測數據，對航空液壓泵的失效物理退化模型進行調制，計算該時刻系統的退化狀態。

步驟5重復步驟2～步驟4，直到當前時刻基于外場數據調制更新得到的液壓泵性能退化狀態超過設定的失效閾值(即xk>xth)，停止壽命預測。

步驟6根據各時刻k對應的退化狀態估計結果和剩余使用壽命預測值RULk，繪制退化狀態估計曲線和剩余使用壽命曲線。

4 實例驗證

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

從圖5中可以看到，航空液壓泵磨損退化過程具有較大的隨機性，這主要是因為液壓泵在工作過程中關鍵摩擦副的磨損率隨載荷工況實時變化，因此即便同一型號的液壓泵其退化過程也不完全相同。對比試驗中的實際觀測數據與本文方法可以看到，本文所提出的數據調制方法可以實現長周期內航空液壓泵性能退化狀態的準確估計。

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

以當前液壓泵運行時間(750 h)作為起始點，隨運行時間繼續增加可以迭代計算出未來任意時刻的剩余使用壽命，其結果如圖7所示。

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

最后，考慮到航空液壓泵在全壽命周期中退化過程存在多階段特性，不同階段適合采用的失效模型可能存在差別，后續將開展多階段失效物理建模，探索有效的基于數據的退化階段突變點檢測算法，將本文方法推廣到多階段退化下的壽命預測；同時，進一步收集其他型號航空液壓泵的在線退化數據，對該方法的在線壽命預測效果開展進一步驗證。

5 結 論

針對航空液壓泵失效機理分析不完備、運行中載荷工況動態變化且外場觀測數據具有不確定性，本文提出了一種全息數據調制更新的航空液壓泵壽命估計方法。通過實際航空液壓泵試驗驗證，結論如下：

1) 基于失效物理與動態數據調制更新的壽命估計理論在挖掘航空液壓泵退化失效機理的同時充分利用觀測信息，通過粒子濾波將外場信息有機調制到性能退化模型中，實現了航空液壓泵全生命周期準確的壽命估計。

2) 針對粒子濾波采樣樣本枯竭問題，基于最優重要性概率密度函數遞推和正則化改進粒子濾波算法，提高了粒子群中的有效粒子比例，增強了算法對觀測不確定性的處理能力。

3) 在狀態估計和RUL預測階段，利用觀測數據對所建立的失效物理模型進行調制，能較好地刻畫航空液壓泵的退化隨機性，RUL預測結果與實際觀測結果高度吻合，可為航空液壓泵維修決策等健康管理活動提供支持。