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

空軍工程大學(xué) 航空航天工程學(xué)院，陜西 西安 710038

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

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

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

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

1 機載燃油傳輸系統(tǒng)試驗平臺

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

表1 產(chǎn)品主要工作參數(shù)Table 1 Product main working parameters

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

圖1 機載燃油轉(zhuǎn)輸系統(tǒng)試驗平臺結(jié)構(gòu)示意圖Fig.1Structure diagram of onboard fuel transfer system experimental platform

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

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

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

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

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

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

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

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

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

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

因此，在隨機效應(yīng)下的Wiener過程模型未知參數(shù)的似然函數(shù)可表示為：

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

3 Bayes剩余壽命在線預(yù)測

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

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

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

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

在獲得不同時間點的性能退化數(shù)據(jù)后，可以得到實時監(jiān)測的性能退化數(shù)據(jù)x的似然函數(shù) 為：

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

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

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

即可得出：

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

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

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

4 基于EM算法的超參數(shù)先驗估計

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

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

采用傳統(tǒng)的極大似然估計法可解得：

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

（1）E步

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

（2）M步

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

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

5 實例分析

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

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

圖3 原始信號監(jiān)測圖Fig.3 Original signal monitoring diagram

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

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

圖5 轉(zhuǎn)換后性能退化數(shù)據(jù)曲線Fig.5 Performance degradation data curve after conversion

首先針對上述產(chǎn)品退化過程數(shù)據(jù)，進行顯著水平為0.1的正態(tài)分布假設(shè)檢驗，結(jié)果接受假設(shè)，表明產(chǎn)品性能退化過程服從Wiener過程。

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

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

表2 參數(shù)更新值Table 2 Parameter update value

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

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

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

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

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

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

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

6 結(jié)論

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

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

[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] 柴天佑. 生產(chǎn)制造全流程優(yōu)化控制對控制與優(yōu)化理論方法的挑戰(zhàn) [J].自動化學(xué)報，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].長沙：國防科技大學(xué)，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] 彭喜元，彭宇，劉大同. 數(shù)據(jù)驅(qū)動的故障預(yù)測[M].哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2015.PENG Xiyuan， PENG Yu， LIU Datong. Data driven fault prediction[M]. Harbin： Harbin Institute of Technology Press，2015. （in Chinese）