孔祥芬,劉敬赟,王杰,唐淑珍

(中國民航大學 航空工程學院,天津 300300)

2020年初新冠肺炎疫情的突然爆發(fā),對全球航空運輸業(yè)造成了極大的負面影響,與航空運輸業(yè)密切相關的航空維修業(yè)也受到了前所未有的沖擊和挑戰(zhàn)[1]。為落實我國對當前民航 “保安全運行、保應急運輸、保風險可控、保精細施策”的要求以及為后續(xù)航空市場的恢復做準備,不少維修單位開始調(diào)整對整機和航空附部件的維修計劃。其中準確獲取航空附部件的可靠性參數(shù)對調(diào)整維修計劃起關鍵性作用。整體驅(qū)動發(fā)電機(IDG)作為飛機電源系統(tǒng)的核心部件[2-3],為機載用電設備提供恒頻交流電,其可靠運行是飛機安全運行的重要保障。因此,結合飛機IDG的故障數(shù)據(jù)特點,尋找可更加客觀準確估計飛機IDG可靠性參數(shù)的方法,對提高飛機IDG的使用可靠性和制定經(jīng)濟合理的維修計劃具有重要意義。

針對飛機IDG這類關鍵航空附部件產(chǎn)品,可靠性要求高而能獲取的歷史故障數(shù)據(jù)較少,屬于小樣本范疇(n≤30)[4]。目前工程實際中普遍采用最小二乘法(Least square regression,LSR) 對具有小樣本故障數(shù)據(jù)特征的部件進行可靠性分析。文獻[5]對競爭失效可靠度模型中的參數(shù)采用最小二乘法進行估計。文獻[6]通過對比分析中位秩公式和平均秩次法的擬合效果,運用最小二乘法對航空裝備進行參數(shù)估計分析;隨著現(xiàn)代應用統(tǒng)計學和計算機科學的發(fā)展,出現(xiàn)了不少針對小樣本數(shù)據(jù)進行可靠分析的方法,文獻[7]針對某機械產(chǎn)品的故障數(shù)據(jù)小樣本刪失的特點,建立了分布識別和參數(shù)估計模型利用支持向量機訓練樣本進而對數(shù)據(jù)進行預測擴充,最后利用擴充后的樣本進行可靠性參數(shù)估計;文獻[8]使用支持向量回歸模型對小樣本故障數(shù)據(jù)進行威布爾分布MTBF估計;文獻[9]將支持向量回歸應用到風電機組的可靠性參數(shù)估計中,并提出該方法更適用于小樣本數(shù)據(jù)的情況;文獻[10]提出一種支持向量機的改進算法——最小二乘支持向量機對服從威布爾分布的發(fā)動機系統(tǒng)進行參數(shù)估計;文獻[11]通過模擬仿真得到滾動軸承的數(shù)據(jù),采用最小二乘支持向量機進行回歸預測其剩余使用壽命。

目前在航空維修可靠性分析方面,較少考慮到部件的故障數(shù)據(jù)為小樣本時對可靠性參數(shù)估計的影響。因此,以某航空公司B737-800飛機IDG部件的故障數(shù)據(jù)(可靠性指標故障間隔時間)為研究對象,以威布爾分布為例,對比分析最小二乘支持向量回歸機、支持向量回歸機和最小二乘法3種方法的可靠性參數(shù)估計結果,找到更加客觀準確的飛機IDG可靠性參數(shù)估計方法。

1 飛機IDG的可靠性參數(shù)估計

1.1 飛機IDG可靠性分布模型

威布爾分布可以描述產(chǎn)品部件處于浴盆曲線各個階段的分布規(guī)律的特征,比其他的分布模型的適用范圍更廣,尤其適用于機電類產(chǎn)品失效的分布規(guī)律描述。研究和實踐結果表明,飛機機電、液壓等系統(tǒng)及其子系統(tǒng)、零部件大都適用威布爾分布來擬合[3,12-14]。

兩參數(shù)威布爾分布的累計失效概率函數(shù)為

(1)

式中:t為部件產(chǎn)品的故障時間;β為形狀參數(shù)；η為尺寸參數(shù)。

一般地,當樣本量n≤20時通常使用中位秩公式來求解累計失效概率密度

(2)

式中:n為部件故障數(shù)據(jù)的樣本量；i為部件故障時間從小到大排序序號。

將威布爾分布模型線性化為

ln(ln[1/(1-F(t))])=βln(t)-βln(η) (3)

設:

(4)

可得線性關系式為

y=ωx+b

(5)

對于一組故障數(shù)據(jù)樣本{ti},i=1,…，n,將數(shù)據(jù)線性化后得到數(shù)據(jù)集T={(x1,y1),(x2,y2),…(xn,yn)}。數(shù)據(jù)集T={(x1,y1),(x2,y2),…，(xn,yn)}線性相關程度通常用相關系數(shù)ρxy來度量，即

(6)

1.2 最小二乘法

最小二乘法LSR通過最小化觀測值和擬合值之間誤差的平方和,來尋找數(shù)據(jù)的最佳函數(shù)匹配,該方法對觀測樣本量無特定要求[15]。當研究對象服從威布爾分布時,參數(shù)估計思路如下:

對于平面直線l:y=ωx+b與數(shù)據(jù)集T={(x1,y1),(x2,y2),…(xn,yn)},量|yi-(b+ωxi)|描述了點(xi,yi)沿y軸方向到直線l的距離。令

(7)

描述數(shù)據(jù)集T={(x1,y1),(x2,y2),…，(xn,yn)}與直線l的偏離程度。顯然,偏離越小越好。

(8)

(9)

解該方程組得:

(10)

由式(4)得β和η的估計值為:

(11)

1.3 支持向量回歸機

支持向量回歸機(Support vector regression,SVR)是通過一個非線性函數(shù)變換,將低維空間映射到高維空間,在高維空間找到線性可分的超平面,使得擬合點到超平面中最遠的樣本點的“距離”最小。該算法在解決小樣本問題中表現(xiàn)出了獨特的優(yōu)勢[16]。

對于給定的數(shù)據(jù)集

T={(xi,yi)，…，(xl,yl)}∈(Rn×y)l

(12)

式中:xi∈Rn；yi∈y=R；i=1,…,l。

在高維空間中,輸入樣本點和輸出樣本點回歸函數(shù)關系表示為

f(x)=ωTφ(x)+bφ:Rn→F

(13)

式中:ω為權向量；b為偏差量;φ(x)為高維空間F中非線性變換函數(shù)。

選擇適當?shù)臄?shù)ε>0和懲罰參數(shù)C>0,支持向量回歸機模型可表示為求解以下優(yōu)化問題:

(14)

(15)

一般引入Lagrange函數(shù)求解此類問題,

(16)

進而得到優(yōu)化問題的對偶形式

(17)

解得

(18)

解得線性方程

(19)

式中k(xi,xj)為滿足Mercer條件的核函數(shù)。

1.4 最小二乘支持向量回歸機

最小二乘支持向量回歸機(Least squares support vector regression,LSSVR)是在SVR的基礎上發(fā)展而來。與SVR不同的是,LSSVR在優(yōu)化問題用平方誤差項代替了SVR中的ε-不敏感損失函數(shù),降低了模型優(yōu)化過程中的復雜性,克服了模型訓練時間較長的問題[17]。

對于給定的數(shù)據(jù)集

LSSVR可表示為求解以下優(yōu)化問題:

(20)

s.t.y-((ω·φ(x)+b)=ξii=1,…,l

式中:γ為正則化常數(shù)；ξi為誤差項。

一般用拉格朗日法求解這類優(yōu)化問題:

(21)

式中αi是拉格朗日乘子。

(22)

引入核函數(shù)K(xi,xj)=φ(xi)·φ(xj)。根據(jù)式(22),優(yōu)化問題轉(zhuǎn)化為求解線性方程:

(23)

最后可得回歸方程為

(24)

1.5 SVR和LSSVR方法下的Weibull分布參數(shù)估計

步驟1 核函數(shù)及其參數(shù)的選擇。SVR和LSSVR模型性能的優(yōu)劣主要取決于所選擇的核函數(shù)及其參數(shù)。其中主要的核函數(shù)有線性核函數(shù),多項式核函數(shù)、RBF核函數(shù)、sigmoid核函數(shù)。本文采用常用的RBF核函數(shù)。

K(xi,xj)=exp(-‖xi-xj‖/σ2),γ>0

(25)

(26)

步驟2 Weibull分布參數(shù)估計。分別使用MATLAB的LIBSVM、LS-SVM工具箱,對SVR和LSSVR模型進行參數(shù)尋優(yōu),得到最優(yōu)回歸直線。然后,結合圖估計法對飛機IDG所服從的威布爾分布進參數(shù)估計:根據(jù)式(4)和式(5)可知,形狀參數(shù)β為所得最優(yōu)回歸直線的斜率;尺寸參數(shù)η又稱特征壽命,即當累計故障概率分布函數(shù)F(t)=63.2%時所對應的壽命。根據(jù)式(4),當F(t)=63.2%時相應的yi=0,故直線yi=0與所得回歸直線交點的橫坐標即為lnη的值。

1.6 參數(shù)估計結果評價

通常,可采用統(tǒng)計學中的統(tǒng)一化的均方根誤差(NRMSE)對不同參數(shù)估計方法所得結果進行定量評價。

(27)

2 案例分析

2.1 飛機IDG部件介紹

飛機整體驅(qū)動發(fā)電機(IDG)又稱組合驅(qū)動發(fā)動機,是目前國內(nèi)外一些大中型客機供電系統(tǒng)中普遍采用的主要電源,如B737系列,A320等。以B737-800機型作為研究對象,其中每架B737-800飛機配備兩臺飛機IDG,分別位于左、右發(fā)動機附件輪箱上。飛機IDG主要由恒速傳動裝置(Constant speed drive,CSD)和交流發(fā)電機兩大部分組成,主要工作原理為:恒速傳動裝置將變化的發(fā)動機轉(zhuǎn)速轉(zhuǎn)變?yōu)楹愣ǖ霓D(zhuǎn)速,然后通過交流發(fā)電機產(chǎn)生115 V/400 Hz的恒頻交流電,以保障機載用電設備的正常運行。飛機IDG的工作原理如圖1所示。

圖1 飛機IDG的工作原理

2.2 故障數(shù)據(jù)預處理

本文統(tǒng)計的故障數(shù)據(jù)來源于國內(nèi)某航空公司45架B737-800飛機電源系統(tǒng)的故障數(shù)據(jù)電子記錄,記錄時間為2015年2月到2018年7月。通過對該組故障數(shù)據(jù)記錄歸類、篩選預處理,從中提取出飛機IDG部件的故障發(fā)生時間、故障描述、維修措施和飛行時間等故障數(shù)據(jù)信息。最終得到飛機IDG部件的故障數(shù)據(jù)即同一機號同一位置的IDG部件相鄰兩次故障間的飛行時間(可靠性指標故障間隔時間),共計15條,部分飛機IDG故障數(shù)據(jù)記錄如表1所示。

表1 飛機IDG的部分故障數(shù)據(jù)記錄

然后,將飛機IDG故障數(shù)據(jù)進行從小到大排序,如表2第1列所示。

表2 故障數(shù)據(jù)、訓練樣本及參數(shù)估計結果

采用式(3)將得到的飛機IDG故障數(shù)據(jù)線性化處理,得到數(shù)據(jù)集T={(x1,y1),(x2,y2),…(xn,yn)}。計算x與y之間的相關系數(shù)得

(28)

根據(jù)結果可知,x與y近似線性關系,說明該B737-800機隊飛機IDG的壽命分布較好地服從二參數(shù)威布爾分布。

2.3 飛機IDG可靠性參數(shù)估計

在MATLAB平臺上分別對SVR和LSSVR進行編程。SVR和LSSVR采用網(wǎng)格搜索法和交叉驗證法對參數(shù)進行優(yōu)化選擇。同時得到最優(yōu)回歸擬合線,分別如圖2和圖3所示,擬合結果如表2第5和第6列所示。

圖2 SVR方法下的最優(yōu)擬合直線

圖3 LSSVR方法下的最優(yōu)擬合直線

根據(jù)圖估計法,SVR和LSSVR所得擬合直線的斜率即為形狀參數(shù)β,通過圖2和圖3中的虛線所示,直線y=0與直線的交點A的橫坐標x=lnη,從而得到尺寸參數(shù)η的值。同時,運用LSR得到參數(shù)估計結果,3種方法的結果如表3所示。

表3 飛機IDG的可靠性參數(shù)估計結果

2.4 不同參數(shù)估計方法的結果對比研究

為了評價3種方法的可靠性參數(shù)估計效果,采用式(22)指標進行可靠性參數(shù)估計結果評價,所得結果如表4第2列所示。

表4 飛機IDG(實際數(shù)據(jù))的可靠性參數(shù)估計結果評價

從以上分析結果可知,LSSVR的誤差最小,參數(shù)估計精度最高;SVR其次;LSR的誤差最大,參數(shù)估計精度最低。

其次,3種可靠性參數(shù)估計方法在MATLAB平臺的運行時間如表4第3列所示。由于LSR無需進行模型內(nèi)的參數(shù)尋優(yōu),速度最快。LSSVR相比SVR運行速度較快。

2.5 不同樣本量下的3種參數(shù)估計方法對比研究

在飛機IDG故障數(shù)據(jù)為小樣本情況下,為觀察樣本量變化時3種可靠性參數(shù)估計方法的穩(wěn)定性[18]。將采集到的2015年～2018年飛機IDG故障數(shù)據(jù),分成一年內(nèi)的故障數(shù)據(jù)(即2015年～2016年)、兩年內(nèi)的故障數(shù)據(jù)(即2015年～2017年)和三年內(nèi)(即2015年～2018)的故障數(shù)據(jù),對各組的故障數(shù)據(jù)分別使用LSR,SVR,LSSVR這3種方法估計飛機IDG可靠性參數(shù),并以NRMSE為指標進行參數(shù)估計結果評價,評價結果如圖4所示。

圖4 不同樣本量下3種方法的參數(shù)估計誤差結果(實際數(shù)據(jù))

由圖4可知,隨著樣本量的減小,3種可靠性參數(shù)估計方法的誤差均呈現(xiàn)出增大的趨勢。但LSSVR的參數(shù)估計誤差增長速度最慢,受樣本量變化的影響最小,具有較高的穩(wěn)定性。當樣本量相同時,LSSVR的參數(shù)估計誤差普遍小于SVR和LSR。即在飛機IDG故障數(shù)據(jù)為小樣本情況下,相比LSR和SVR,LSSVR所得參數(shù)估計結果普遍具有較高的精度。

3 蒙特卡羅仿真驗證

為進一步驗證3種可靠性參數(shù)估計方法對比分析結果的準確性,對服從兩參數(shù)威布爾分布的飛機IDG部件進行蒙特卡羅仿真[19],其中形狀參數(shù)和尺度參數(shù)分別設定為β=1.5,η=10 000。

首先,對3種方法的可靠性參數(shù)估計結果及其評價進行驗證。將飛機IDG部件的運行時間設定為20 000FH進行多次仿真,得到小樣本故障數(shù)據(jù){ti}。將該組故障數(shù)據(jù){ti}線性化后得到小樣本數(shù)據(jù)集T={(x1,y1),(x2,y2),…,(xn,yn)}后,采用3種方法進行可靠性參數(shù)估計。3種方法下的可靠性參數(shù)估計結果及其評價如表5所示。

表5 飛機IDG(仿真數(shù)據(jù))的可靠性參數(shù)估計結果及其評價

然后,對樣本量變化時3種方法的可靠性參數(shù)估計結果的穩(wěn)定性進行驗證。將故障樣本量分別設定為5、10、15、20進行蒙特卡羅仿真。根據(jù)多次仿真抽樣得到故障數(shù)據(jù)樣本,采用3種參數(shù)估計方法計算得到不同樣本量下3種方法的可靠性參數(shù)估計評價結果,如圖5所示。

圖5 不同樣本量下3種方法的參數(shù)估計誤差結果(仿真數(shù)據(jù))

由兩組仿真結果可知,當飛機IDG故障數(shù)據(jù)為小樣本時,LSSVR算法下的可靠性參數(shù)估計誤差最小,SVM次之,LSR較大,與案例分析中所得結果一致;3種可靠性參數(shù)估計方法在MTALAB平臺的運行時間、當樣本量變化時的穩(wěn)定性所得結果也與案例分析中所得結果基本一致。

4 結論

本文針對飛機IDG可靠性高而故障數(shù)據(jù)為小樣本的特點,分別采用LSSVR、SVR和LSR對飛機IDG進行可靠性參數(shù)估計。通過實際案例分析以及蒙特卡羅仿真驗證,對比分析了3種方法的可靠性參數(shù)估計結果?？傻靡韵陆Y論:

在飛機IDG 可靠性參數(shù)估計中,當獲取的故障數(shù)據(jù)屬于小樣本(n≤30)、壽命分布服從兩參數(shù)威布爾分布時,LSSVR的參數(shù)估計精度最高,SVR次之,LSR最低;在運行時間方面,LSR由于無需進行模型內(nèi)的參數(shù)尋優(yōu)速度較快,LSSVR次之,SVR 較差。當樣本量變化時,隨著樣本量的減小,3種參數(shù)估計方法的精度均呈現(xiàn)減小的趨勢。但LSSVR的參數(shù)估計精度受樣本量變化的影響最小,具有較高的穩(wěn)定性。可知,當飛機IDG的故障數(shù)據(jù)為小樣本時,采用LSSVR可獲得較為客觀準確的可靠性參數(shù)估計結果。

根據(jù)本文研究結果可為航空維修單位制定飛機IDG相關維修計劃提供理論依據(jù),以節(jié)約維修成本提高飛機的使用效率;此外,該小樣本故障數(shù)據(jù)情況下的可靠性參數(shù)估計方法還可以為飛機其他系統(tǒng)、附部件以及其他領域的高可靠性裝備提供參考依據(jù)。