基于比例風(fēng)險(xiǎn)模型的飛機(jī)部件剩余壽命預(yù)測

何 浪，魯 亮+，陸錦鶴，范圣嫻，趙立新，張宇翔

(1.中國民航大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，天津 300300；(2.北京飛機(jī)維修工程有限公司 附件/起落架大修產(chǎn)品事業(yè)部，北京 100621)

0 引 言

全球的空中交通運(yùn)輸量正在快速增加，預(yù)計(jì)到2037年民航客機(jī)總數(shù)將達(dá)到45 265架[1]。在商業(yè)運(yùn)營中，民用飛機(jī)的維修成本通常占航空公司直接運(yùn)營成本的10%～20%[2]，準(zhǔn)確預(yù)測飛機(jī)部件從當(dāng)前狀態(tài)到失效狀態(tài)的剩余壽命可以提早預(yù)警飛機(jī)部件故障，并合理制定維修計(jì)劃，從而有效降低維修成本，這在航空公司的安全運(yùn)行和成本控制方面至關(guān)重要。

在現(xiàn)有工業(yè)設(shè)備的剩余壽命預(yù)測方法中[3，4]，基于數(shù)據(jù)驅(qū)動(dòng)的方法不依賴設(shè)備的失效原理，而是通過對設(shè)備運(yùn)行數(shù)據(jù)進(jìn)行建模分析來預(yù)測設(shè)備剩余壽命，因而被廣泛使用。該方法可大致分為基于退化數(shù)據(jù)的方法和基于失效數(shù)據(jù)的方法兩類[5]。前者雖然可以通過多傳感器信息提高剩余壽命預(yù)測的精度[6-8]，但是增加了系統(tǒng)的復(fù)雜性及成本[9]；而后者則需選擇合適的壽命分布，且需要足夠的失效數(shù)據(jù)進(jìn)行分布擬合。文獻(xiàn)[10，11]分別解決了壽命分布的選擇問題和分布擬合的數(shù)據(jù)量問題，但沒有考慮飛機(jī)各部件壽命的個(gè)體差異。事實(shí)上，一些飛機(jī)部件在進(jìn)行拆換時(shí)會(huì)從備件倉庫中選出調(diào)試好的部件進(jìn)行更換，這些部件可能經(jīng)歷了不同的飛行時(shí)長、維修模式以及維修次數(shù)等。此外，文獻(xiàn)[12]中指出，飛機(jī)年齡、飛機(jī)發(fā)動(dòng)機(jī)型號等飛機(jī)參數(shù)也可能成為影響飛機(jī)故障概率的潛在因素。若將這些因素與部件的失效數(shù)據(jù)結(jié)合，可進(jìn)一步提高壽命預(yù)測結(jié)果的精確度，并使其具有良好的可解釋性。

比例風(fēng)險(xiǎn)模型(proportional hazard model，PHM)[13]作為一種可以分析多個(gè)變量與生存時(shí)間之間關(guān)系的統(tǒng)計(jì)回歸模型，能處理可靠性數(shù)據(jù)且不需要提前對生存分布進(jìn)行假設(shè)，正逐漸被應(yīng)用于可靠性領(lǐng)域[14，15]與航空發(fā)動(dòng)機(jī)的剩余壽命預(yù)測研究中。文獻(xiàn)[16-18]基于該模型分析了多個(gè)參數(shù)對飛機(jī)發(fā)動(dòng)機(jī)、輔助動(dòng)力裝置(auxiliary power unit，APU)等部件壽命的影響，但仍利用部件的狀態(tài)監(jiān)測數(shù)據(jù)進(jìn)行壽命預(yù)測，忽略了部件狀態(tài)以及部件所處環(huán)境對其壽命的影響。此外，飛機(jī)由大量的子系統(tǒng)和部件組成，很難直接獲取到部分部件可用的狀態(tài)監(jiān)測數(shù)據(jù)。

本文彌補(bǔ)了基于失效數(shù)據(jù)方法在考慮部件個(gè)體差異時(shí)存在的不足，同時(shí)規(guī)避了基于狀態(tài)監(jiān)測數(shù)據(jù)的方法存在的問題，綜合運(yùn)用飛機(jī)部件的歷史維修數(shù)據(jù)和歷史失效數(shù)據(jù)，提出了一種基于比例風(fēng)險(xiǎn)模型的飛機(jī)部件剩余壽命預(yù)測方法，并通過對比實(shí)驗(yàn)驗(yàn)證其有效性。

1 飛機(jī)部件剩余壽命預(yù)測

1.1 可靠性分析

可靠性分析是剩余壽命預(yù)測的基礎(chǔ)[19]?？煽啃允侵覆考谝?guī)定條件和規(guī)定時(shí)間間隔內(nèi)執(zhí)行指定功能的能力，涉及失效概率F、可靠度R和失效率h等統(tǒng)計(jì)指標(biāo)。

設(shè)壽命T為部件從完好狀態(tài)到失效狀態(tài)所經(jīng)歷的時(shí)間，失效(故障)概率F表示部件工作時(shí)間T小于t時(shí)失效的概率，可表示為

(1)

其中，f(x) 表示T的概率密度函數(shù)。在工程實(shí)踐中，該變量通常通過實(shí)驗(yàn)評估確定。

可靠度R指設(shè)備工作時(shí)間T大于t的概率，可用式(2)表示

(2)

失效率h描述為設(shè)備工作時(shí)間到t，并在時(shí)刻t后單位時(shí)間內(nèi)失效的概率

(3)

由式(2)可知f(x)=-dR(t)/dt， 代入式(3)中，則有

(4)

對等式(4)左右兩邊同時(shí)積分可推出可靠度R與失效率h之間的關(guān)系

(5)

1.2 比例風(fēng)險(xiǎn)模型

比例風(fēng)險(xiǎn)模型的基本形式定義為

h(t|X)=h0(t)exp(βX)

(6)

其中，h0(t) 表示各特征變量均為0時(shí)的基本失效率函數(shù)，h(t|X) 表示各特征變量的值均固定時(shí)的失效率函數(shù)，X=(X1，X2，…，Xp)T為影響部件剩余壽命的特征向量，β=(β1，β2，…，βp) 為相關(guān)特征對應(yīng)的系數(shù)向量，表示特征變量對部件壽命的影響程度，p表示特征的個(gè)數(shù)。

1.2.1 特征選擇

為了保留每個(gè)特征的實(shí)際意義并避免由自變量之間存在的相關(guān)關(guān)系，即多重共線性造成的特征冗余，采用逐步回歸法對主要特征進(jìn)行篩選。逐步回歸法根據(jù)檢驗(yàn)的顯著性對特征變量逐個(gè)分析，通過與設(shè)定的顯著性水平對比，以此來決定對該特征變量進(jìn)行引入或者剔除。其步驟如下：

(1)對于未引入的特征變量，從中選擇與因變量相關(guān)程度最高的特征變量并建立回歸方程，進(jìn)行顯著性檢驗(yàn)。若通過檢驗(yàn)，則重復(fù)該步驟；否則轉(zhuǎn)步驟(2)。

(2)對于已引入的特征變量，進(jìn)行顯著性檢驗(yàn)，若檢驗(yàn)不通過，則剔除該特征變量，確保每次引入新的特征變量之前回歸方程只包含顯著的變量，轉(zhuǎn)步驟(1)；否則結(jié)束。

此外，在引入特征變量的過程中，其先后順序并不代表其在最終回歸方程中的重要程度，后引入的相對重要的特征變量則會(huì)因?yàn)槎嘀毓簿€性而被剔除。因此，需要在進(jìn)行逐步回歸前進(jìn)行多重共線性的診斷。

多重共線性的診斷方法主要分為兩類，一類是經(jīng)驗(yàn)式的診斷方法，包括觀察特征變量的相關(guān)系數(shù)矩陣、通過專業(yè)知識或一般經(jīng)驗(yàn)分析回歸系數(shù)的代數(shù)符號以及某一個(gè)特征變量是另一部特征變量的近似線性組合等。另一類為統(tǒng)計(jì)診斷方法，其最主要的方法為使用方差膨脹因子。特征變量Xp的方差膨脹因子記為VIFp，公式如下

(7)

1.2.2 比例風(fēng)險(xiǎn)模型假定的檢驗(yàn)

比例風(fēng)險(xiǎn)模型要求各特征變量的作用效應(yīng)不隨時(shí)間的變化而變化，故在建立比例風(fēng)險(xiǎn)模型前需要進(jìn)行比例風(fēng)險(xiǎn)假設(shè)驗(yàn)證?，F(xiàn)有比例風(fēng)險(xiǎn)假設(shè)方法可分為圖示法和假設(shè)檢驗(yàn)法，其中圖示法具有清晰、實(shí)用等特點(diǎn)，主要包括Cox &K-M比較法[13]、累積風(fēng)險(xiǎn)函數(shù)圖示法[20]以及Schoenfeld殘差法[21]。

Cox &K-M比較法最早由Cox提出，通過將比例風(fēng)險(xiǎn)模型和Kaplan-Meier估計(jì)的各組生存曲線進(jìn)行比較，若曲線的趨勢一致且無交叉，則滿足比例風(fēng)險(xiǎn)假定。

累積風(fēng)險(xiǎn)函數(shù)圖示法可通過觀察對數(shù)累積風(fēng)險(xiǎn)函數(shù)對時(shí)間的趨勢圖，若各組曲線互相平行，則滿足比例風(fēng)險(xiǎn)假定。

根據(jù)式(6)以及累積風(fēng)險(xiǎn)函數(shù)可知

(8)

對等式(8)兩邊取對數(shù)

lnH(t|X)=lnH0(t)+βX

(9)

由上式可知，若不同組間的對數(shù)累積風(fēng)險(xiǎn)函數(shù)差值lnH(t|X1)-lnH(t|X2)=β(X1-X2) 為常數(shù)，則對應(yīng)圖像中各組曲線應(yīng)保持平行或重疊。

Schoenfeld殘差法最早由Schoenfeld提出，在ti時(shí)刻部件的偏殘差ri定義如下

ri=(ri1，ri2，…，rik)

(10)

(11)

1.2.3 參數(shù)估計(jì)

由式(6)可知，比例風(fēng)險(xiǎn)模型可以看作基本失效率函數(shù)的非參數(shù)部分和一個(gè)與失效因素線性相關(guān)的參數(shù)部分的乘積。

對于非參數(shù)部分，不同部件的生存時(shí)間對應(yīng)的分布不同。常見的壽命分布函數(shù)有正態(tài)分布、指數(shù)分布、威布爾分布、對數(shù)正態(tài)分布、伽馬分布等，需選擇最適合部件的分布函數(shù)對其生存時(shí)間分布進(jìn)行擬合。首先通過極大似然估計(jì)得到相應(yīng)分布的參數(shù)，然后根據(jù)AD檢驗(yàn)統(tǒng)計(jì)量的大小選擇合適的分布類型。

對于參數(shù)部分，協(xié)變量的參數(shù)估計(jì)通常可以用以下偏似然方程進(jìn)行求解

(12)

其中，n表示觀測樣本的數(shù)量，δi為狀態(tài)變量。當(dāng)δi=0時(shí)表示部件在故障前被拆除，δi=1時(shí)則表示部件故障時(shí)被拆除。

式(12)的對數(shù)方程為

(13)

lnL(β) 的一階導(dǎo)數(shù)為

(14)

根據(jù)式(13)和式(14)，可采用梯度下降法對參數(shù)進(jìn)行估計(jì)。另一種辦法則是繼續(xù)求出式(13)的二階偏導(dǎo)，然后利用Newton-Raphson方法進(jìn)行參數(shù)估計(jì)。該方法雖然收斂快，但是需要計(jì)算二階偏導(dǎo)以及矩陣的逆，計(jì)算復(fù)雜度偏大。

1.3 剩余壽命預(yù)測

在獲得參數(shù)β的估計(jì)值后，將式(6)代入式(5)中可以推出在X條件下的可靠度函數(shù)

(15)

對等式(15)兩邊取對數(shù)

(16)

剩余壽命的平均值為

(17)

根據(jù)式(16)、式(17)可得到t時(shí)刻在X條件下的平均剩余壽命

(18)

本文采取固定時(shí)間點(diǎn)預(yù)測的方式預(yù)測部件的剩余壽命。如圖1所示，根據(jù)部件的平均更換間隔時(shí)間設(shè)置預(yù)測點(diǎn)t1，t2，…，tn， 在各預(yù)測點(diǎn)對飛機(jī)輪胎進(jìn)行剩余壽命的預(yù)測。

圖1 預(yù)測間隔

2 數(shù)據(jù)集描述

本文使用的數(shù)據(jù)集收集自北京飛機(jī)維修工程有限公司。數(shù)據(jù)集包含858個(gè)部件，共85 259條維修記錄，涉及163架B737 NG型號的飛機(jī)，其時(shí)間跨度為2016-2019年。數(shù)據(jù)屬性包括收貨日期、部件編號、維修模式代碼、維修費(fèi)用、裝機(jī)日期、拆卸日期、部件自出廠到現(xiàn)在的累計(jì)飛行小時(shí)數(shù)(time since new，TSN)、部件自上次大修以來的累計(jì)飛行小時(shí)數(shù)(time since overhaul，TSO)、部件本次裝機(jī)到拆卸之間的累計(jì)飛行小時(shí)數(shù)(time since install，TSI)，部件本次裝機(jī)到拆卸之間的累計(jì)飛行小時(shí)數(shù)、飛機(jī)編號、發(fā)動(dòng)機(jī)型號、飛機(jī)制造日期等。其中，TSI為因變量，表示部件的壽命。根據(jù)數(shù)據(jù)的可用性以及部件的經(jīng)濟(jì)效益，在本文實(shí)驗(yàn)中選擇主輪的裝拆記錄進(jìn)行數(shù)據(jù)分析，并從部件自身以及部件所處環(huán)境兩個(gè)方面進(jìn)行特征提取。對于主輪自身，TSN、TSO、部件累計(jì)維修次數(shù)(TIMES)以及維修的模式(MODE)都可能會(huì)影響部件的壽命。飛機(jī)在減速剎車時(shí)，主輪要面臨急劇的溫度變化，可能對輪胎的壽命造成影響，因此可以考慮環(huán)境溫度變化的影響，故將季度(Quarter)作為主輪的備選特征。此外，考慮飛機(jī)的結(jié)構(gòu)的老化以及發(fā)動(dòng)機(jī)推力的不同，也有可能影響主輪的壽命，故將主輪所在飛機(jī)的年齡(AGE)和發(fā)動(dòng)機(jī)型號(Engine)作為備選特征。表1給出了主輪備選特征描述。

表1 特征描述

表2中為各連續(xù)型變量的描述性統(tǒng)計(jì)。為了更好地評估這些參數(shù)對TSI的影響，需要對這些連續(xù)變量進(jìn)行歸一化。歸一化方法如下

表2 連續(xù)變量的描述性統(tǒng)計(jì)

(19)

其中，Z1max、Z1min、Z2max、Z2min、Z3max、Z3min、Z4max和Z4min分別表示TSN、TSO、AGE以及TIMES的最大最小值。

對于分類變量，采用箱形圖展示各分類變量關(guān)于TSI的數(shù)據(jù)分布，如圖2所示。

圖2 分類變量的數(shù)據(jù)分布

在圖2(a)中，維修模式包括新件(NEW)、修理(repair，RP)和大修(overhaul，OH)。

3 模型驗(yàn)證與分析

3.1 主輪與TSI分布分析

通過常見壽命分布函數(shù)對主輪的TSI進(jìn)行擬合，然后根據(jù)AD檢驗(yàn)統(tǒng)計(jì)量選用AD值較小的分布函數(shù)。表3中給出了擬合效果較好的分布的參數(shù)估計(jì)結(jié)果。

表3 各分布參數(shù)估計(jì)結(jié)果

圖3為各分布的擬合結(jié)果，圖中的直方圖為主輪TSI的數(shù)據(jù)分布。從圖中可以看出這幾種分布曲線與主輪的數(shù)據(jù)分布都較為接近，擬合的效果較好。

圖3 分布曲線擬合

3.2 檢驗(yàn)與回歸結(jié)果

在逐步回歸篩選的過程中，TIMES與MODE被剔除。其中，TIMES與TSN存在多重共線性，而MODE的顯著性檢驗(yàn)未通過。通過逐步回歸篩選后，TSN、TSO、AGE、Engine和Quarter是篩選出的對TSI影響最顯著的特征變量。

對主輪進(jìn)行比例風(fēng)險(xiǎn)回歸建模之前，需要對所選特征變量進(jìn)行比例風(fēng)險(xiǎn)假設(shè)驗(yàn)證。選擇基于累積風(fēng)險(xiǎn)函數(shù)圖示法和Schoenfeld殘差法分別對分類變量和連續(xù)變量進(jìn)行比例風(fēng)險(xiǎn)假設(shè)驗(yàn)證。由圖4可以看出發(fā)動(dòng)機(jī)型號、季度等分類變量的對數(shù)累積風(fēng)險(xiǎn)函數(shù)關(guān)于TSI的曲線都基本平行或等距，說明這幾個(gè)參數(shù)滿足比例風(fēng)險(xiǎn)假定。

圖4 分類變量的對數(shù)累積風(fēng)險(xiǎn)函數(shù)

圖5為各連續(xù)變量的Schoenfeld殘差圖，可以看出圖中的擬合曲線基本以0為中心波動(dòng)，滿足比例風(fēng)險(xiǎn)假定。

圖5 連續(xù)變量的殘差

表4中，維修模式、發(fā)動(dòng)機(jī)和季度類型用虛擬變量表示。虛擬變量是指取值為0或1的變量，用于表示分類變量中某些屬性狀態(tài)是否存在。根據(jù)虛擬變量的設(shè)置規(guī)則，需要選取一種屬性作為分析的基礎(chǔ)，其它屬性則作為虛擬變量進(jìn)行分析。以“發(fā)動(dòng)機(jī)型號”為例，CFM56-7B、CFM56-7B22、CFM56-7B24和CFM56-7B26這4種型號分別用Engine_1、Engine_2、Engine_3、Engine_4表示，將Engine_1作為參照變量，則Engine_2、Engine_3、Engine_4作為虛擬變量。各變量的系數(shù)為β中對應(yīng)的值，P值表示該變量在設(shè)定的顯著性水平(α=0.05)下是否顯著，即該變量對基本失效率h0(t) 的影響是否顯著。HR表示風(fēng)險(xiǎn)比，表示在其它變量保持一定時(shí)，該變量對基本失效率的影響程度。

表4 比例風(fēng)險(xiǎn)回歸結(jié)果

從表4中可以看出飛機(jī)機(jī)齡會(huì)使輪胎的故障風(fēng)險(xiǎn)增加；裝有CFM56-7B22和CFM56-7B24的飛機(jī)上的輪胎更容易故障；第二季度和第三季度輪胎故障的風(fēng)險(xiǎn)更大。飛機(jī)輪胎值得注意的是TSN和TSO的HR值意味著當(dāng)它們增加時(shí)輪胎的故障風(fēng)險(xiǎn)會(huì)降低，其主要的原因可能是由數(shù)據(jù)篩選的結(jié)果造成的。為了完整跟蹤飛機(jī)輪胎從新件到經(jīng)過數(shù)次維修后的連續(xù)狀態(tài)，本文在數(shù)據(jù)篩選時(shí)選擇了2016-2019期間使用的新輪胎的維修記錄，因此數(shù)據(jù)中輪胎的累計(jì)飛行小時(shí)數(shù)不足。其次，航空公司現(xiàn)有的維修策略會(huì)要求一些部件在經(jīng)過幾次修理后就會(huì)經(jīng)歷一次大修，使得部件維修記錄中的TSO值也不會(huì)太大。

3.3 剩余壽命預(yù)測結(jié)果分析

3.3.1 比例風(fēng)險(xiǎn)模型預(yù)測結(jié)果

表5 比例風(fēng)險(xiǎn)模型預(yù)測結(jié)果

3.3.2 對比方法預(yù)測結(jié)果

為了更好地評估比例風(fēng)險(xiǎn)模型在預(yù)測部件剩余壽命時(shí)的效果，采用了其它回歸模型進(jìn)行對比，即將比例風(fēng)險(xiǎn)模型中的特征作為其它回歸模型的參數(shù)。由于比例風(fēng)險(xiǎn)模型假設(shè)特征變量以乘積的形式影響基本失效率，包含了部件的壽命分布，因此本文在與其它傳統(tǒng)回歸方法對比時(shí)，考慮將部件的壽命信息加入到回歸模型中。在文獻(xiàn)[22]中，Kadhem等用威布爾分布對歷史風(fēng)速數(shù)據(jù)進(jìn)行擬合，并用威布爾分布進(jìn)行隨機(jī)采樣生成一個(gè)風(fēng)速數(shù)據(jù)，通過將該數(shù)據(jù)與其它影響風(fēng)速的特征作為人工神經(jīng)網(wǎng)絡(luò)的輸入特征對風(fēng)速進(jìn)行預(yù)測。此外，在文獻(xiàn)[23]中，He等提出了一種加性平均剩余壽命模型，分析了潛在風(fēng)險(xiǎn)因素對平均剩余壽命的影響。結(jié)合上述文獻(xiàn)的數(shù)據(jù)處理思想，本文對分類回歸樹、多層感知機(jī)回歸、多元線性回歸、支持向量回歸以及隨機(jī)森林回歸等其它回歸模型分別嘗試了如下兩種不同的改進(jìn)方法。

根據(jù)式(17)，在t1，t2，…，tn時(shí)刻計(jì)算出對應(yīng)的平均剩余壽命m(ti)， 其中i=1，2，…，n。

(1)將m(ti) 加入到特征向量X中，作為第p+1個(gè)特征。然后將新的特征向量作為其它回歸模型的輸入。

(2)計(jì)算m(ti) 與實(shí)際剩余壽命的差值并作為因變量，然后將特征向量X作為其它回歸模型的輸入。

按照上述改進(jìn)方法，使用訓(xùn)練集對各模型進(jìn)行訓(xùn)練，測試集來驗(yàn)證各模型的預(yù)測能力。本文將均方根誤差(RMSE)、平均絕對誤差(MAE)、平均絕對百分誤差(MAPE)作為評估本文方法以及其它方法的指標(biāo)。這些指標(biāo)是預(yù)測任務(wù)中常用的評價(jià)指標(biāo)，數(shù)學(xué)表達(dá)式如下

(20)

(21)

(22)

表6給出了比例風(fēng)險(xiǎn)模型與其它回歸模型壽命預(yù)測的結(jié)果。其中L表示未改進(jìn)的回歸模型，L1表示第一種改進(jìn)，即將m(ti) 加入到特征向量X中，將新的特征向量作為其它回歸模型的輸入。L2表示第二種改進(jìn)，即計(jì)算m(ti) 與實(shí)際剩余壽命的差值作為因變量，將特征向量X作為輸入。

表6 各方法對比結(jié)果

可以看出改進(jìn)后的回歸模型預(yù)測精度有了明顯提升。L1和L2相較L的RMSE分別降低了40.63以及46.83，MAE則分別降低了34.23以及38.24。從MAPE來看，L1和L2方法相較L在預(yù)測精度上分別平均提升了8.44%以及9.4%。而比例風(fēng)險(xiǎn)模型在3個(gè)指標(biāo)上都要優(yōu)于其它回歸模型，在整體預(yù)測效果上，比例風(fēng)險(xiǎn)模型>L2>L1>L。

4 結(jié)束語

本文通過對飛機(jī)部件的維修數(shù)據(jù)進(jìn)行分析，提出了基于比例風(fēng)險(xiǎn)模型的飛機(jī)部件剩余壽命預(yù)測方法。通過對飛機(jī)主輪進(jìn)行實(shí)例驗(yàn)證，分析了影響飛機(jī)主輪壽命的主要因素以及影響程度，可為部件維修策略的調(diào)整提供參考。此外，通過對其它回歸模型進(jìn)行改進(jìn)，提供了比例風(fēng)險(xiǎn)模型的對比方法。實(shí)驗(yàn)結(jié)果表明，改進(jìn)后的回歸模型的預(yù)測精度平均提升了8.92%，同時(shí)比例風(fēng)險(xiǎn)模型的整體預(yù)測效果要優(yōu)于改進(jìn)后的回歸模型。