基于Wiener過(guò)程的航空燃油泵壽命預(yù)測(cè)

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

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

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

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

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

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

1 機(jī)載燃油傳輸系統(tǒng)試驗(yàn)平臺(tái)

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

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

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

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

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

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

2 基于Wiener過(guò)程的性能退化模型

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

式中：μ為漂移系數(shù)，σ為擴(kuò)散系數(shù)，B（t）為標(biāo)準(zhǔn)布朗運(yùn)動(dòng)。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

即可得出：

μ，ω的后驗(yàn)邊緣密度分別為：

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

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

4 基于EM算法的超參數(shù)先驗(yàn)估計(jì)

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

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

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

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

（1）E步

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

（2）M步

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

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

5 實(shí)例分析

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

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

圖3 原始信號(hào)監(jiān)測(cè)圖Fig.3 Original signal monitoring diagram

由圖4可看出，所提取的特征信號(hào)值在BüC段具有明顯的線性下降趨勢(shì)，因此，提取BüC段特征信號(hào)數(shù)據(jù)，選取性能退化閾值為61Pa。進(jìn)行剩余壽命預(yù)測(cè)分析建模。以同批9個(gè)產(chǎn)品常規(guī)退化試驗(yàn)數(shù)據(jù)作為個(gè)體剩余壽命預(yù)測(cè)先驗(yàn)信息，取同批單獨(dú)一個(gè)產(chǎn)品退化數(shù)據(jù)作為其實(shí)時(shí)更新信息。為方便模型計(jì)算，對(duì)原始數(shù)據(jù)進(jìn)行初始值歸零轉(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

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

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

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

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

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

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

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

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

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

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

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

6 結(jié)論

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

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

[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)化控制對(duì)控制與優(yōu)化理論方法的挑戰(zhàn) [J].自動(dòng)化學(xué)報(bào)，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過(guò)程的可靠性建模方法研究[D].長(zhǎng)沙：國(guó)防科技大學(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ū)動(dòng)的故障預(yù)測(cè)[M].哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2015.PENG Xiyuan， PENG Yu， LIU Datong. Data driven fault prediction[M]. Harbin： Harbin Institute of Technology Press，2015. （in Chinese）