基于多元退化數(shù)據(jù)的航空機電系統(tǒng)競爭故障預(yù)測?

孟 蕾 許愛強 董 超

（1.海軍航空工程學(xué)院飛行器工程系 煙臺 264001）（2.海軍航空工程學(xué)院科研部 煙臺 264001）（3.海軍航空工程學(xué)院基礎(chǔ)部 煙臺 264001）

1 引言

航空機電系統(tǒng)發(fā)生故障模式有兩種，退化故障和突發(fā)故障，其最終故障是兩種故障模式競爭的結(jié)果。目前國內(nèi)外文獻(xiàn)研究競爭故障問題，都是在突發(fā)故障與退化故障兩者不相連的情況下，著重于從單邊退化角度展開論述，是退化量、突發(fā)故障單一化的集中體現(xiàn)，退化量的多元化始終未能在學(xué)界占主導(dǎo)。對于這樣的問題，本文主要是通過航空機電系統(tǒng)提出了退化特性，構(gòu)建了多元退化量的航空機電系統(tǒng)的競爭和預(yù)測模型。

本文中，以最小二乘支持向量機預(yù)測法為基本計算工具，簡要的對退化參數(shù)做了基本測算，如非線性、小樣本等；退化量、故障間存在著非常特殊的聯(lián)系，為尋找隱含性質(zhì)，引入了“位置-尺度模型”，得出突發(fā)故障與退化量的相關(guān)參數(shù)，進(jìn)而根據(jù)競爭故障預(yù)測模型得到了航空機電系統(tǒng)未來一段時間內(nèi)的競爭故障概率。立足案例，在比較研究之后，認(rèn)為本文在方法論上具有優(yōu)勢，滿足精度要求。

2 構(gòu)建航空機電系統(tǒng)競爭故障預(yù)測模型

假設(shè)存在的性能特征參數(shù)數(shù)量是n，xi(t)、Ui和Li分別代表i(i=1，…，n)個性能特征參數(shù)的退化量與上、下閾值。測試起始時間為t，則有G(x，βi)（其中βi=(β1，β2，…，βk)，為該分布的參數(shù)向量），對應(yīng)的概率密度函數(shù)g(x，βi) 可表示為

工業(yè)部門對故障的概念進(jìn)行了界定，對性能特征參數(shù)的退化量的標(biāo)準(zhǔn)值和閾值范圍，有著嚴(yán)格的退化量限制，如果大于了規(guī)定的標(biāo)準(zhǔn)，則故障難免。假設(shè)第i個性能特征參數(shù)在發(fā)生了退化故障，則可依據(jù)下式計算其可能性：

因為航空機電系統(tǒng)本身就具有n個特定性能的參數(shù)，任一參數(shù)出現(xiàn)了退化故障，將會導(dǎo)致整個系統(tǒng)會產(chǎn)生故障。所以，可以把系統(tǒng)的故障視為由n個性能特征參數(shù)的退化故障通過競爭而產(chǎn)生的。假設(shè)所有的性能和特征參數(shù)之間保持相對獨立，則在t時刻航空機電系統(tǒng)的退化故障的概率可以通過下面的公式進(jìn)行估算。

突發(fā)故障這一指標(biāo)不但和時刻t等有關(guān)聯(lián)，而且與t時刻某性能特征參數(shù)的退化量xi(t) 等指標(biāo)有關(guān)聯(lián)。在航空機電系統(tǒng)中，若突發(fā)故障必然由退化量引起，故障出現(xiàn)時間是Th，且在產(chǎn)生關(guān)系的臨界時間Th會伴隨失效率（λi(t，x)）的同時出現(xiàn)，則Th的可靠性關(guān)系式為

航空機電系統(tǒng)在t時刻的突發(fā)故障概率可表示為

對于競爭故障出現(xiàn)的概率，其算式是

3 航空機電系統(tǒng)競爭故障的預(yù)測流程

航空機電系統(tǒng)進(jìn)行競爭故障預(yù)測的具體流程如下：首先確定性能退化參數(shù)與突發(fā)故障參數(shù)的分布類型；然后預(yù)測退化的未知參數(shù)，確定下一階段性能退化參數(shù)的分布函數(shù)；在考慮突發(fā)故障與退化量的相關(guān)度的情況下，求出突發(fā)故障預(yù)測模型的未知參數(shù)；最后對航空機電系統(tǒng)退化與突發(fā)競爭故障概率進(jìn)行預(yù)測。

3.1 分布函數(shù)的確定

性能退化參數(shù)分布類型的檢驗，本文采取Копмогоров檢驗方法。假設(shè)其分布函數(shù)G(x，β)的分布類型與時間t無關(guān)，只是參數(shù)向量β隨時間t的變化而改變，則G(x，β)即可看作β與時間t的函數(shù)，可表示為G(x，β(t))。在確定G(x，β(t))的分布種類前，不妨借鑒工程實踐經(jīng)驗，并且針對假設(shè)進(jìn)行分布擬合檢驗，最終確定出G(x，β(t))的分布類型。綜合考慮航空機電系統(tǒng)的雙向退化與單邊退化。通常情況下，都會將正態(tài)分布N(x;μ，σ2)作為理想化選擇，相應(yīng)的密度函數(shù)如式（6）所示：

式中，μ和σ2分別表示t時刻性能退化數(shù)據(jù)的均值和方差，則性能退化數(shù)據(jù)的分布參數(shù)向量β=(μ，σ2)，易知，β與時間t相關(guān)。

對航空機電系統(tǒng)性能退化數(shù)據(jù)分布類型的檢驗就變?yōu)橐粋€正態(tài)性檢驗過程，可利用樣本(X1，X2，…，Xn)對總體的分布是否服從N(μ，σ2)進(jìn)行檢驗。 樣本從小到大排列成：X(1)≤X(2)…≤X(n)。利用式（7）計算統(tǒng)計量W，式中ak(W)的值可查表獲得［1］。

然后在給定顯著性水平α與樣本容量n下，通過查表可得Wn，α。若W＜Wα，則拒絕H0，否則接受H0。

針對突發(fā)故障數(shù)據(jù)，本文根據(jù)工程經(jīng)驗假設(shè)其服從威布爾分布。檢驗方法仍采用Копмогоро對假設(shè)的航空機電系統(tǒng)突發(fā)故障數(shù)據(jù)分布類型進(jìn)行驗證，驗證過程與性能退化數(shù)據(jù)分布類型的檢驗類似。

3.2 參數(shù)的預(yù)測

3.2.1 退化故障參數(shù)預(yù)測

3.1節(jié)得出性能退化數(shù)據(jù)的分布類型，從而得出性能退化分布參數(shù)。根據(jù)歷史性能退化數(shù)據(jù)求解分布參數(shù)的歷史數(shù)據(jù)值，求解將來參數(shù)，確定分布函數(shù)，對退化故障、發(fā)生幾率做預(yù)測。

航空機電系統(tǒng)的性能退化數(shù)據(jù)的獲取主要基于狀態(tài)監(jiān)測系統(tǒng)或者是加速退化實驗，本文研究的性能退化數(shù)據(jù)主要是來源于狀態(tài)監(jiān)測系統(tǒng)獲取的數(shù)據(jù)。通過測試取得了相應(yīng)的性能數(shù)據(jù)是和時間相關(guān)的指標(biāo)，分布參數(shù)也將會隨著時間的變化而調(diào)整。為此，可從時間序列角度審視。根據(jù)上文相關(guān)論述，可以很清晰得知：從時間序列角度展開預(yù)測性能退化故障活動將顯得非常的簡便與合理。

在上述領(lǐng)域中，其主要的方法為神經(jīng)網(wǎng)絡(luò)、時間序列分析、支持向量機等［2～3］。由于航空機電系統(tǒng)性能大多數(shù)分布參數(shù)都集中體現(xiàn)出了非線性性質(zhì)，這決定了只有應(yīng)用統(tǒng)計學(xué)原理方可得出理想結(jié)果。支付向量機（SupportVector Machine，SVM）在處理小樣本、高維數(shù)的非線性問題方面具有很強的適用性，但是SVM算法隨著樣本量的增大，運算速度會降低［4］。最小二乘支持向量機（LV-SVM）算法既節(jié)約了計算時間又降低了計算難度。因此，本文選用LV-SVM預(yù)測模型對航空機電系統(tǒng)各性能退化數(shù)據(jù)的分布參數(shù)時間序列進(jìn)行預(yù)測［5］。

上式（8）中，αi為拉格朗日乘子，K為核函數(shù)，b?R為常數(shù)。由于，則預(yù)測模型可表示為

3.2.2 突發(fā)故障參數(shù)預(yù)測

本文將突發(fā)故障時間當(dāng)作響應(yīng)變量，在回歸函數(shù)中，退化量并不需要被排除。這說明退化量、故障出現(xiàn)時間的研究能夠被放置于回歸模型中討論。借助回歸分析，還需要尋找理想的模型——位置-尺度模型［6］，而且它在退化量x等參數(shù)的基礎(chǔ)上對Y=lnT的分布進(jìn)行了描述，模型可以使用下面的關(guān)系式進(jìn)行表述：

當(dāng)中，e的分布和x沒有聯(lián)系。H(e)-e的可靠度函數(shù)。μ(x)-位置參數(shù)，σ＞0-不變的尺度參數(shù)，x，Y的可靠度函數(shù)為

針對其中的突發(fā)故障時間Th，使用Yh=lnTh的分布，則Th=exрYh的可靠度函數(shù)可表示為

式中，α(x)=exр(μ(x))，δ=1/σ，S(w)=H(ln(w))。利用該模型，式（6）可進(jìn)一步表示為

由3.1節(jié)分析可知，航空機電突發(fā)故障時間Th服從威布爾分布［7］?，F(xiàn)對給定退化量x下Yh=lnTh的分布進(jìn)行分析，由于Th吻合于威布爾分布，則Yh=lnTh服可以吻合于極值分布［8］，其密度函數(shù)使用以下的公式進(jìn)行計算：

式中，位置參數(shù)μ(x)=lnη(x)，尺度參數(shù)σ=1/m，若采用標(biāo)準(zhǔn)極值分布［9］進(jìn)行表示，則Yh如式（15）所示：

式中，e的密度函數(shù)是exр(e-exр(e))。μ(x)以線性為主，令μ(x)=γ1+γ2·x，σ=γ3，則式（16）可表示為

在綜合研究了諸多的機電系統(tǒng)后，發(fā)現(xiàn)對象系統(tǒng)的故障并不是單一存在的，而是事件、時間的綜合表現(xiàn)，即tj(j=1，2，…，K)（退化量）必然在某個tj(j=1，2，…，K)（突發(fā)故障時間）中出現(xiàn)，此類關(guān)系被稱為似然關(guān)系［10-12］，由式（17）函數(shù)來闡述突發(fā)故障數(shù)據(jù)：

4 實例分析

某海軍航空飛行團生產(chǎn)的某型機電系統(tǒng)綜合ATS測試數(shù)據(jù)為例，隨機抽取10組機電系統(tǒng)某特征參數(shù)作為樣本進(jìn)行故障預(yù)測。其測試信息進(jìn)行測試的時間從2010年持續(xù)到2016年，主要是針對性能退化數(shù)據(jù)和突發(fā)故障數(shù)據(jù)進(jìn)行測試，因此可根據(jù)2010年～2015年的測試信息對機電系統(tǒng)2016年的故障概率進(jìn)行預(yù)測，而且把結(jié)果和機電系統(tǒng)的真實情況進(jìn)行對比，以以便可以檢驗本次系統(tǒng)設(shè)計的有效性和針對性。

結(jié)合在工程中的相關(guān)經(jīng)驗，可假定航空機電系統(tǒng)的一些性能特征的參數(shù)退化中的數(shù)據(jù)可以吻合于正態(tài)分布。以提取的十組機電系統(tǒng)中的特征性能參數(shù)2010年的相關(guān)的數(shù)據(jù)為實例，并最終得到以下退化數(shù)據(jù)（抽取而來）：9.93，10.61，11.00，10.43，12.10，10.22，10.53，10.55，10.67，11.12。按照分布擬合法對抽取值開展驗算，并最終得到分布狀態(tài)。

在集中整理了測試信息后，設(shè)該性能退化數(shù)據(jù)為X，則檢驗X服從正態(tài)分布與否就變?yōu)闄z驗假設(shè)H0:X～N(μ，σ2)是否成立。按升序排列性能退化數(shù)據(jù)，并對其進(jìn)行正態(tài)性檢驗，為方便計算，列表如下。

表1 性能退化數(shù)據(jù)的正態(tài)性檢驗

由表1可計算得到

將上述計算結(jié)果代入式（7），可得：W=0.8907。

若取顯著性水平α=0.05，則通過查閱W檢驗統(tǒng)計量W的α分位數(shù)表Wα可得W10，0.05=0.842。由于W＜W10，0.05，所以在顯著性水平α=0.05的情況下，可以接受原假設(shè)H0，即可確定2010年測試時某的一些性能參數(shù)的特征和正態(tài)分布相吻合。

按照提取的十組機電系統(tǒng)某的一些性能參數(shù)，針對2010～2015的性能參數(shù)進(jìn)行退化數(shù)據(jù)的正態(tài)性檢驗，結(jié)果表明所有的年份的數(shù)據(jù)都吻合正態(tài)分布。此外，正態(tài)性測試也適合剩余數(shù)據(jù)，相關(guān)操作方式如上。通過測算，發(fā)現(xiàn)如果結(jié)果出現(xiàn)了α=0.05的情況，其余性能特征參數(shù)的性能中的退化數(shù)據(jù)也吻合于正態(tài)分布。同理，結(jié)合工程中的經(jīng)驗，假定機電系統(tǒng)的突發(fā)故障數(shù)據(jù)仍然使用Копмогоров方法對其中的分布類型進(jìn)行證實，可得機電系統(tǒng)的突發(fā)故障參數(shù)服從威布爾分布。

由3.2節(jié)分析結(jié)果得出，航空機電系統(tǒng)性能退化參數(shù)可以吻合于正態(tài)分布，所以在針對性能參數(shù)進(jìn)行預(yù)測的過程中，也只是需要針對預(yù)測均值和方差等指標(biāo)進(jìn)行測試。

經(jīng)過對數(shù)據(jù)的分析處理，從2007～2015年中，提取了十組機電系統(tǒng)某性能特征參數(shù)，并對此類參數(shù)進(jìn)行特殊處理，得到方差、均值。其時間序列為：（10.7160，10.7530，10.7960，10.8000，10.8160，10.8240，10.9340，11.0020，11.0870）和（0.3534，0.3771，0.6871，0.5370，0.5003，0.4050，1.0234，0.5728，0.5980）。首先預(yù)測均值μ，選取m=1，相空間重構(gòu)均值序列的前5個數(shù)據(jù)，可得樣本對Xμ=[10.7160，10.7530，10.7960，10.8000，10.8160，10.8240]T、Yμ=[10.7530，10.7960，10.8000，10.8160，10.8240，10.9340]T。選取RBF函數(shù)作為核函數(shù)，設(shè)定參數(shù)γ和δ的搜索區(qū)域分別為[0.1，500]和[0.01，100]，采用交叉檢驗和網(wǎng)格搜索算法進(jìn)行求解，可得最優(yōu)參數(shù)δ2=3.2463，γ=116.4963。使用LS-SVM訓(xùn)練對樣本執(zhí)行輸入和輸出等操作，可以取得均值μ的預(yù)測的回歸直線，如下圖。得到訓(xùn)練完畢的LS-SVM回歸函數(shù)之后，即可開始預(yù)測輸入樣本10.7530， 10.7960， 10.8000， 10.8160， 10.8240，10.9340]T。預(yù)測結(jié)果顯示實測值與預(yù)測值基本吻合：同理，可以預(yù)測性能退化數(shù)據(jù)的方差時間序列，結(jié)果顯示絕對誤差小于-0.05，相對誤差小于0.009。

圖1 均值預(yù)測回歸曲線

分析2007年～2015年的故障數(shù)據(jù)可以發(fā)現(xiàn)，提取了十組機電系統(tǒng)的記錄，持續(xù)的時間為2008年～2015年，其故障的總數(shù)分別為2,2,2,2,1,1,1,1。因為對于機電系統(tǒng)進(jìn)行了實時的測試，所以需要以天為單位進(jìn)行測試，所以提取的十組數(shù)據(jù)的突發(fā)時間可以表示為具體的天數(shù)。以某退化特征參數(shù)的退化量為例，與突發(fā)故障時間的對應(yīng)關(guān)系可得，同時可以得出突發(fā)故障與退化量的相關(guān)參數(shù)的估計值。明確了數(shù)據(jù)的分類模型和相應(yīng)的參數(shù)，即可根據(jù)式（13）對機電系統(tǒng)在2015年和2016年的故障概率進(jìn)行預(yù)測。結(jié)果發(fā)現(xiàn)，實際進(jìn)行測評的結(jié)果和預(yù)測值是基本上一致的，所以建立的故障預(yù)測模型是合理而有效的。

5 結(jié)語

文章中，借助于航空機電系統(tǒng)中關(guān)系復(fù)雜的退化性質(zhì)，集中論述了退化發(fā)生頻率、故障出現(xiàn)時間的數(shù)學(xué)關(guān)聯(lián)，引入航空系統(tǒng)的故障和預(yù)測模型，而且針對退化數(shù)據(jù)中的分布參數(shù)的小樣本、非線性的特征，使用最小支持向量機預(yù)測的算法針對故障進(jìn)行了預(yù)測，基于退化量與突發(fā)故障的相關(guān)性，預(yù)測了突發(fā)故障預(yù)測模型相關(guān)參數(shù)。最后實例分析結(jié)果表明：本文設(shè)計的方法預(yù)測精度更高，更加合理有效。

［1］吳翊，李永樂，胡慶軍.應(yīng)用數(shù)理統(tǒng)計［M］.北京：國防科技大學(xué)出版社，2008.

［2］張廣明，袁宇浩，龔松建.基于改進(jìn)最小二乘支持向量機方法的短期風(fēng)速預(yù)測［J］.上海交通大學(xué)學(xué)報，2011，45（8）：1125-1129.

［3］叢林虎，徐廷學(xué)，楊繼坤，董琪.導(dǎo)彈退化故障預(yù)測方法研究［J］.電光與控制，2014，21（5）：78-82.

［4］洪杰，韓磊，苗學(xué)問，馬艷紅.基于支持向量機的滾動軸承狀態(tài)壽命評估［J］.北京航空航天大學(xué)學(xué)報，2010，36（8）：896-899.

［5］唐杰明，劉俊勇，楊可等.基于灰色模型和最小二乘支持向量機的電力短期負(fù)荷組合預(yù)測［J］.電網(wǎng)技術(shù)，2009，33（3）：63-65.

［6］Lawless J F.壽命數(shù)據(jù)中的統(tǒng)計模型與方法［M］.北京：中國統(tǒng)計出版社，1998.

［7］蘇春，張恒.基于性能退化數(shù)據(jù)和競爭失效分析的可靠性評估［J］.機械強度，2011，33（2）：196-200.

［8］孫紹輝.面向預(yù)知維修的航空發(fā)動機可靠性評估與預(yù)測研究［D］.南京：南京航空航天大學(xué)，2013.

［9］GOLDSTEIN M，WOOFF D.Bayes linear statistics，theo?ry and methods［M］.New York：Wiley Publisher，2007.

［10］喻天翔，孫玉秋，張祖明.多模式失效的機械零件可靠度計算新理論［J］.機械工程學(xué)報，2003（03）：34-35.

［11］蘇春，張恒.基于性能退化數(shù)據(jù)和競爭失效分析的可靠性評估［J］.機械強度，2011（02）：22-23.

［12］李方方，趙英凱，顏昕.基于Matlab的最小二乘支持向量機的工具箱及其應(yīng)用［J］.計算機應(yīng)用，2006（26）：358-360.