基于Wiener過程的航空燃油泵壽命預測

張先航，李曙林，常飛，杜旭，尹俊杰，肖堯

空軍工程大學 航空航天工程學院，陜西 西安 710038

目前，故障預測與健康管理（PHM）技術已成為我軍航空裝備重點發展的關鍵技術，其核心功能——長壽命、高可靠性的航空裝備關鍵部位剩余壽命預測能力決定了PHM技術在維修保障中的價值。對于剩余壽命預測的研究，國內外學者在過去十年中給予了充分的關注，并取得了長足的發展[1]。Percht在關于PHM的專著[2]中，將現有的剩余壽命預測方法分為三類：基于機理模型的方法、數據驅動的方法和前兩者融合的方法。參考文獻[3]指出“復雜的工業過程往往具有多變量、強耦合、強非線性、大延時、生產邊界條件變化頻繁、動態特性隨工況變化、難以用數學模型描述等綜合復雜特性”，進而認為基于數據驅動的方法為這類設備的控制、決策、優化提供了可行的途徑。

數據驅動的方法主要包括人工智能和統計數據驅動的方法（基于隨機模型建模數據的方法）[1]。人工智能方法一般利用得到的數據，通過機器學習擬合設備性能變量的演化規律，進而通過外推得到失效閾值以實現剩余壽命預測[4]，但是其難以得到體現剩余壽命預測隨機不確定性的概率分布函數。從另一個角度來看，剩余壽命預測是基于當前獲得的監測信息對未來的失效事件進行預測，因此，預測的結果不可避免地具有不確定性。相比之下，基于統計數據驅動的方法以概率論統計理論為基礎，通過統計或隨機模型得到剩余壽命的概率分布，便于量化剩余壽命預測結果的不確定性。

由于維納（Wiener）過程具有良好的數學特性，自20世紀90年代以來，該隨機過程已被廣泛應用于設備可靠性預測與壽命分析[5，6]。Doksum和Hoyland[6]首次將Wiener過程用于加速模型建模，利用Wiener過程建模加速退化數據，推斷正常應力水平下的設備壽命；Tseng[7]等基于Wiener過程對LED燈的剩余壽命預測問題進行了研究，但該文獻只考慮了設備當前時刻的狀態監測信息；Peng[8]考慮了Wiener過程的漂移系數為隨機變量時的退化建模的問題，推導出壽命分布的顯示解，但該方法未利用設備運行過程中的監測數據。

航空燃油泵是飛機燃油系統的核心部件，其性能好壞直接影響著飛行安全。但作為高可靠、長壽命產品，航空燃油泵在短期內難以觀察到失效時間數據。在航空燃油泵的故障診斷研究中， 國內外廣泛采用振動信號或者壓力信號[9]作為狀態監測和診斷的信息源，2013年，印度學者Muralidharan 和Sugumaran[10]使用振動信號，結合小波分析進行特征提取，運用模糊邏輯分類方法，實現對故障的分類識別。因此，選擇與產品壽命相關的出口壓力作為性能退化特征量，采用Wiener過程建立性能退化模型，運用貝葉斯（Bayes）統計推斷的方法實現模型的實時在線更新，得到燃油泵運行過程中的實時剩余壽命預測分布情況，對燃油系統和飛機[1，2]PHM具有重要意義。

1 機載燃油傳輸系統試驗平臺

采用南京機電液壓工程研究中心提供的某型離心式交流電動泵為研究對象，該型號燃油泵是機載燃油系統的核心附件之一，主要為散熱分系統和供油箱供油。該燃油泵工作介質：航空噴氣燃料RP-3；介質溫度：-60~ +85℃；工作參數見表1。

表1 產品主要工作參數Table 1 Product main working parameters

搭建的機載燃油轉輸系統試驗平臺如圖1所示，主要包括儲油箱、供油箱、離心式燃油泵、三個振動加速度計傳感器（CA-YD-182-10）、一個壓力傳感器（CY-YZ-001）、數據采集設備（億恒MI-7008）等，其中，振動傳感器采用磁吸座的方式安裝于電機殼上相互垂直的三個位置。壓力傳感器采用轉接管安裝在燃油泵出口60~80cm處，如圖2所示。

圖1 機載燃油轉輸系統試驗平臺結構示意圖Fig.1Structure diagram of onboard fuel transfer system experimental platform

圖2 壓力傳感器安裝實物圖Fig.2 Physical diagram of pressure sensor installation

測試時，首先打開閥門，將儲油箱注滿燃油，接通泵電源，泵啟動后連續工作，待泵運行穩定后，采集泵振動信號和出口的壓力信號，信號采集結束后，關閉電源和閥門。

2 基于Wiener過程的性能退化模型

采用隨機效應下的Wiener過程對燃油泵性能退化數據進行建模假設產品在t時刻的性能退化量為X(t)，X(t)為一維連續隨機過程且服從Wiener過程，則其數學模型可以表示為：

式中：μ為漂移系數，σ為擴散系數，B（t）為標準布朗運動。

考慮同批產品不同個體之間的隨機性，認為在此條件下μ和σ都是隨機變量，聯合共軛先驗分布為正態—逆Gamma分布[11]，記 ω=1/σ2， 則已知 ω 服從Gamma分布，設ω～Gamma(a,b)，則在給定 ω 的條件下，u|ω～Normal(c,d/ω)。由Wiener過程的定義[11]及式（1），在給定 μ，σ2的條件下，可知：

則對任意給定的時刻t>s>0，在t-s時間段內的產品性能退化量的增量X(t)-X(s)可表示為：

假設有N個產品進行退化試驗，在M個不同時刻退化量進行觀測。記X(tij)表示第i個產品在時刻j處的測量值（i=1，2，…，N；j=1，2，…，M），通過觀測可得到如下性能退化量的測量數據：

為方便計算， ，其中Δtij表示第i個樣品在時刻tij-1到tij之間的時間增量，ΔXij表示第i個樣品在時刻tij-1到tij之間的退化量的增量，則可得到性能退化增量數據為：

由Wiener過程的性質可知，在給定μ和σ2的條件下，有，因此，性能退化量增量ΔXij的概率密度函數為：

因此，在隨機效應下的Wiener過程模型未知參數的似然函數可表示為：

設產品的退化失效閾值為l，產品首次達到l的壽命T服從逆高斯分布[13]，當產品退化量X(t)小于l時，令

3 Bayes剩余壽命在線預測

為提高預測的準確性和合理性，采用Bayes方法對退化模型的參數進行實時更新。

由Bayes定義可知，參數空間θ上任一概率分布都稱作先驗分布。用π(θ)表示θ的先驗分布，這里π(θ)是隨機變量θ的概率函數。θ為連續型隨機變量，π(θ)為θ的密度函數。

在獲得樣本X后，θ的后驗分布就是在給定X=x條件下θ的條件分布，記為π(θ|x)。對于連續性隨機變量，其概率密度函數為：

式中： 為X，θ的 聯 合 分 布 ，而為x的邊緣分布，已知參數μ，ω聯合分布情況，可得退化模型的參數先驗分布為：

在獲得不同時間點的性能退化數據后，可以得到實時監測的性能退化數據x的似然函數 為：

根據 Bayes公式可以將 π(μ，ω)更新為 π(μ，ω|X)，由式（10）可知，π(θ，x)正比例于似然和先驗的乘積，即：

將式（11）和式（12）代入式（13），可得其參數μ，ω后驗分布：

由于分布參數的后驗分布與其共軛先驗分布的函數形式相同，通過對比式（11）與式（14）的函數形式，得到其超參數后驗估計值 a*，b*，c*，d*分別為：

即可得出：

μ，ω的后驗邊緣密度分別為：

在平方損失函數下，μ，ω的Bayes估計為 ：

將 代入式（8）即可得到產品的剩余壽命分布，剩余壽命期望為 。

4 基于EM算法的超參數先驗估計

由式（13）可知，要想求得μ，ω后驗均值，首先要求出參數a，b，c，d的先驗估計值。這里采用EM方法求解其超參數估計值。

以收集到的N個燃油泵在M個監測時刻所得到的性能退化數據值為先驗信息，設第i個產品在第j次監測得到的性能退化數據為Xij，監測時間點為tij，引入兩組潛在數據(μ1，μ2，…，μM)和 (ω1，ω2，…，ωM)。建立包含參數 μi，ωi和超參數a，b，c，d的完全似然函數。

采用傳統的極大似然估計法可解得：

由式（12）~式（15）中可知， 含有隱含參數μ，ω的項，即 項。本文采用EM 算法求解帶隱含變量的上述難題。

（1）E步

根據EM方法，E步驟涉及計算，根據 ωi服從 Gamma分布，ui|ωi服從正態分布，可推出ui服從非中心t分布，各項的期望值可根據以上分布函數的統計特性[12]推導得出。設{at，bt，ct，dt}為迭代t次后超參數估計值，解得在迭代t+1步后各項包含潛在參數 μi，ωi的期望為：

（2）M步

將在E步中計算得出的各項包含潛在數據μi，ωi的參數項代入式（24）~式（27），可得t+1步迭代后各參數{a，b，c，d}估計值：

當相鄰兩步的誤差達到預設精度時，即可停止迭代。輸出最后一次的參數估計值，即為模型的超參數估計值。

5 實例分析

在實驗室對某型航空燃油泵進行模擬工作電壓下的試驗，采集試驗過程中記錄燃油泵出口處壓力傳感器的輸出值，從其性能退化監測數據時間序列圖來看，其性能退化呈現出分階段、震蕩退化的形式，且原始監測數據中噪聲較多，非常不易對其性能退化過程進行預測分析。

因此，采用小波分析的方法，提取其性能下降趨勢特征信號，以db5小波對原始信號進行5層分解，提取其具有性能退化趨勢的特征信號，如圖3所示。

圖3 原始信號監測圖Fig.3 Original signal monitoring diagram

由圖4可看出，所提取的特征信號值在BüC段具有明顯的線性下降趨勢，因此，提取BüC段特征信號數據，選取性能退化閾值為61Pa。進行剩余壽命預測分析建模。以同批9個產品常規退化試驗數據作為個體剩余壽命預測先驗信息，取同批單獨一個產品退化數據作為其實時更新信息。為方便模型計算，對原始數據進行初始值歸零轉換的預處理，得到轉化后的性能退化曲線如圖5所示。

圖4 小波分析性能特征提取圖Fig.4 Wavelet analysis performance feature extraction drawing

圖5 轉換后性能退化數據曲線Fig.5 Performance degradation data curve after conversion

首先針對上述產品退化過程數據，進行顯著水平為0.1的正態分布假設檢驗，結果接受假設，表明產品性能退化過程服從Wiener過程。

利用采集的9組燃油泵性能退化數據，運用EM算法進行超參數估計，得到a，b，c，d超參數先驗值。設超參數 a，b，c，d初始值為 3.0，0.01，0.01，0.01，當精度 e小于10-8時則停止EM迭代。得到超參數先驗分布為：a=3.127；b=0.000371；c=0.00337；d=0.00132。

通過實時監測的性能退化數據運用Bayes方法對參數進行實時更新，在得到的a，b，c，d后驗均值的基礎上，通過式（19）得到μ，ω的后驗均值，從而實現對分布參數μσ的在線更新，見表 2。

表2 參數更新值Table 2 Parameter update value

將μ，σ2代入式（8）可得到其剩余壽命概率密度函數，如圖6、圖7所示。

圖6 剩余壽命預測效果三維圖Fig.6 Residual life prediction effect 3D diagram

圖7 剩余壽命預測效果二維圖Fig.7 Residual life prediction effect 2D diagram

由圖8、圖9可知，基于Wiener過程的航空燃油泵壽命預測均值路徑與實際退化情況基本相吻合，隨著監測時間的推進，其預測結果的概率密度置信區間越來越小，即其預測不確定性也越來越小，精度越來越高。

通過對剩余壽命預測不確定性進行量化分析，得出其剩余壽命預測水平在不同置信度[13]下的置信區間，見表3。

表3 不同壽命預測置信水平下置信區間Table 3 Confidence interval under different life expectancy confidence levels

由表3可知，剩余壽命預測的置信區間99%置信度下的壽命預測區間為（20.57，40.62）。在監測后期可以保持在20h左右，預測區間與預測剩余壽命均值之比在10%左右?？梢愿鶕O測航空燃油泵實際維修保障情況，為其選擇設計單位時間消耗費用最低且使用度最高的剩余壽命預測精度，從而達到維修保障的效能最佳。

6 結論

本文采用小波分析的方法提取了燃油泵具有性能退化趨勢的信號，運用基于隨機效應的Wiener過程對其性能退化數據進行建模分析，通過Bayes方法實現預測模型參數的在線更新，得到其剩余壽命預測結果。

結果表明，基于隨機效應的Wiener過程的壽命預測的結果與實際退化過程基本吻合，且預測結果不確定性隨運行時間的推移越來越低，在接近壽命閾值的預測階段其不確定率可以保持在10%左右。本文所采用的分析方法可以為開展PHM剩余壽命預測效能分析提供一定的思路與借鑒，得出的預測結果對于視情維修保障活動的開展具有一定的指導意義與價值。

[1] Si X， Wang W， Hu C， et al. Remaining useful life estimation–A review on the statiscal data driven approaches[J].European Journal of Operational Research， 2011，213（1）：1-14.

[2] Pecht M. Prognostics and health management of electronics[M].New Jersey：Wiley Online Library，2008.

[3] 柴天佑. 生產制造全流程優化控制對控制與優化理論方法的挑戰 [J].自動化學報，2009，35（6）：641-649.CHAI Tianyou. Challenges of optimal control for plantwide prdution processes in terms of control and optimization theories[J].Journal of Automation， 2009， 35（6）： 641-649. （in Chinese）

[4] Jardine A， Lin D， Banjevic D. A review on marchinery diagnostic and prognostics implementing condition-based maintenance[J]. Mechanical Systems and Signal Processing，2006， 20（07）： 1483-1510.

[5] Heng A， Zhang S， Tan A， et al. Rotaing machinery prognostics，state of the art， challenges and opportunities[J]. Mechanical Systems and Signal Processing， 2009， 23（3）：724-739.

[6] Doksum K， Hoyland A. Models for variable-stress accelerated life testing experiments based on wiener processes and the inverse Gaussian distribution [J]. Technometrics， 1992， 34（1）：74-82.

[7] Tseng S， Tang J，Ku I. Determination of burn-in parameters and residual life for highly reliable products[J]. Navel Research Logistics， 2003， 50（1）： 1-4.

[8] Peng C，Tseng S. Mis-specification analysis of linear degration models[J]. IEEE Transactions on Reliability， 2009， 58（3）：1161-1170.

[9] Hancoke K M， Zhang Q. A hybrid approach to hydraulic vane pump condition monitoring and fault detection[J]. Transactions of a the Asabe， 2006， 49（4）：1203-1211.

[10] Muralidharan V， Sugumaran V. Rough set based rule learning and fuzzy classification of wavelet features for fault diagnosis of monoblock centrifugal pump[J]. Measurement， 2013， 46（9）：3057-3063.

[11] 彭寶華. 基于Wiener過程的可靠性建模方法研究[D].長沙：國防科技大學，2010.PENG Baohua. Rearch on reliability modeling method based on Wiener process[D]. Changsha： National Defense University，2010. （in Chinese）

[12] Wang X L， Cheng Z J， Guo B. Residual life forecasting of degradation data[J]. Journal of Multivariate Analysis， 2010， 101（2）：340-351.

[13] 彭喜元，彭宇，劉大同. 數據驅動的故障預測[M].哈爾濱：哈爾濱工業大學出版社，2015.PENG Xiyuan， PENG Yu， LIU Datong. Data driven fault prediction[M]. Harbin： Harbin Institute of Technology Press，2015. （in Chinese）