孟 蕾 許愛強 董 超

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

1 引言

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

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

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

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

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

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

突發(fā)故障這一指標不但和時刻t等有關(guān)聯(lián)，而且與t時刻某性能特征參數(shù)的退化量xi(t) 等指標有關(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)競爭故障的預測流程

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

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

性能退化參數(shù)分布類型的檢驗，本文采取Копмогоров檢驗方法。假設其分布函數(shù)G(x，β)的分布類型與時間t無關(guān)，只是參數(shù)向量β隨時間t的變化而改變，則G(x，β)即可看作β與時間t的函數(shù)，可表示為G(x，β(t))。在確定G(x，β(t))的分布種類前，不妨借鑒工程實踐經(jīng)驗，并且針對假設進行分布擬合檢驗，最終確定出G(x，β(t))的分布類型。綜合考慮航空機電系統(tǒng)的雙向退化與單邊退化。通常情況下，都會將正態(tài)分布N(x;μ，σ2)作為理想化選擇，相應的密度函數(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)進行檢驗。 樣本從小到大排列成：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)驗假設其服從威布爾分布。檢驗方法仍采用Копмогоро對假設的航空機電系統(tǒng)突發(fā)故障數(shù)據(jù)分布類型進行驗證，驗證過程與性能退化數(shù)據(jù)分布類型的檢驗類似。

3.2 參數(shù)的預測

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

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

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

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

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

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

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

當中，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）可進一步表示為

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

式中，位置參數(shù)μ(x)=lnη(x)，尺度參數(shù)σ=1/m，若采用標準極值分布［9］進行表示，則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ù)作為樣本進行故障預測。其測試信息進行測試的時間從2010年持續(xù)到2016年，主要是針對性能退化數(shù)據(jù)和突發(fā)故障數(shù)據(jù)進行測試，因此可根據(jù)2010年～2015年的測試信息對機電系統(tǒng)2016年的故障概率進行預測，而且把結(jié)果和機電系統(tǒng)的真實情況進行對比，以以便可以檢驗本次系統(tǒng)設計的有效性和針對性。

結(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ù)據(jù)為X，則檢驗X服從正態(tài)分布與否就變?yōu)闄z驗假設H0:X～N(μ，σ2)是否成立。按升序排列性能退化數(shù)據(jù)，并對其進行正態(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的情況下，可以接受原假設H0，即可確定2010年測試時某的一些性能參數(shù)的特征和正態(tài)分布相吻合。

按照提取的十組機電系統(tǒng)某的一些性能參數(shù)，針對2010～2015的性能參數(shù)進行退化數(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ù)仍然使用Копмогоров方法對其中的分布類型進行證實，可得機電系統(tǒng)的突發(fā)故障參數(shù)服從威布爾分布。

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

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

圖1 均值預測回歸曲線

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

5 結(jié)語

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

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

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

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

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

［5］唐杰明，劉俊勇，楊可等.基于灰色模型和最小二乘支持向量機的電力短期負荷組合預測［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］孫紹輝.面向預知維修的航空發(fā)動機可靠性評估與預測研究［D］.南京：南京航空航天大學，2013.

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

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

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

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