王巖韜，李景良，谷潤平

中國民航大學(xué)空管學(xué)院，天津 300300

近年不安全事件表明，民航工作高質(zhì)量發(fā)展依賴于航班運(yùn)行風(fēng)險(xiǎn)防控能力的提升[1].針對(duì)航班運(yùn)行風(fēng)險(xiǎn)影響因素繁雜、隨時(shí)間變化、等級(jí)界定困難等問題，歐美國家最早從1991年開始航班風(fēng)險(xiǎn)評(píng)價(jià)研究[2].從近年研究成果看，研究理論和手段仍使用事故樹等傳統(tǒng)方法[3-4].而在實(shí)踐方面，自2007—2012年間美國、加拿大民航局研發(fā)了Flight risk assessment tool (FRAT)[5]、RA check[6]等飛行風(fēng)險(xiǎn)評(píng)估工具后，至今未有可信的使用報(bào)告和精度分析；除此之外，近年也未見其他突破.在國內(nèi)，航班風(fēng)險(xiǎn)評(píng)估研究最早從1999年起步[7]，而后王巖韜等[8-9]建立了具有指導(dǎo)性、實(shí)用性的航班運(yùn)行風(fēng)險(xiǎn)評(píng)估體系，并運(yùn)用多算法協(xié)作的方式將航班風(fēng)險(xiǎn)分類正確率提高至95%.對(duì)于民航非動(dòng)態(tài)風(fēng)險(xiǎn)評(píng)估，上述研究已取得較好成果.但將前續(xù)成果應(yīng)用于民航一線運(yùn)行后發(fā)現(xiàn)，航班評(píng)估技術(shù)只能做到隨航班要素快速變化而變化，評(píng)估結(jié)果代表已有狀態(tài)，而風(fēng)險(xiǎn)管控工作需要時(shí)間提前量，因此風(fēng)險(xiǎn)預(yù)測技術(shù)才是風(fēng)險(xiǎn)防控實(shí)質(zhì)所亟需.

2012年，歐盟開展航空運(yùn)行主動(dòng)性安全管理項(xiàng)目（Proactive safety performance for operations,PROSPERO），旨在通過對(duì)航空運(yùn)輸系統(tǒng)的風(fēng)險(xiǎn)主動(dòng)預(yù)測，降低人為差錯(cuò)和空中交通事故率，但詳細(xì)方案和結(jié)果報(bào)告一直未見報(bào)道，并且最新消息表明該項(xiàng)目已于2015年戛然而止[10].2018年Lalis等[11]使用線性回歸模型和ARMA模型對(duì)航空安全績效指標(biāo)進(jìn)行預(yù)測，但以月為單位的、本身就波動(dòng)很小的時(shí)間序列嚴(yán)重削弱了其實(shí)用價(jià)值.2019年，Zhang與Mahadevans[12]融合支持向量機(jī)和深度神經(jīng)網(wǎng)絡(luò)兩種算法應(yīng)用于航空風(fēng)險(xiǎn)預(yù)測中，對(duì)單一事件預(yù)測精度達(dá)到了81%，但其過程實(shí)質(zhì)上與風(fēng)險(xiǎn)評(píng)估類似，預(yù)測周期過短，預(yù)測能力不足.在國內(nèi)，王巖韜團(tuán)隊(duì)于2019年提出了基于動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)的航班運(yùn)行風(fēng)險(xiǎn)預(yù)測模型[13]，但該模型主要針對(duì)單個(gè)航班的進(jìn)近、著陸階段的風(fēng)險(xiǎn)預(yù)測，在預(yù)測周期很短的情況下精度達(dá)到了80.4%，無法體現(xiàn)航空公司的整體運(yùn)行狀況.由上，航班風(fēng)險(xiǎn)預(yù)測相關(guān)研究較少，至今缺乏精度較高的技術(shù)和方法，是民航安全領(lǐng)域尚未有效解決的關(guān)鍵問題之一.

深入分析其他領(lǐng)域的研究成果，發(fā)現(xiàn)多變量混沌時(shí)間序列預(yù)測方法是一種較好的短期預(yù)測技術(shù)，早期應(yīng)用于經(jīng)典的混沌序列中，證明了其比單變量預(yù)測更優(yōu)[14-16].近年來，該方法廣泛應(yīng)用于工程領(lǐng)域，典型研究成果包括：2013年，侯公羽等[17]對(duì)煤礦斜井的盾構(gòu)小時(shí)施工風(fēng)險(xiǎn)進(jìn)行預(yù)測，均方根誤差（Root mean square error, RMSE）為 0.4754；2016年，鐘儀華等[18]對(duì)油田日產(chǎn)油量進(jìn)行預(yù)測，均方根誤差為0.0061；2017年，杜柳青等[19]對(duì)數(shù)控機(jī)床的日運(yùn)動(dòng)精度進(jìn)行預(yù)測，均方根誤差低于0.02；2018年，張淑清等[20]對(duì)電力日負(fù)荷進(jìn)行預(yù)測，平均誤差為0.813%；2018年，黃發(fā)明等[21]對(duì)滑坡的月波動(dòng)項(xiàng)位移進(jìn)行預(yù)測，均方根誤差為23.71.上述工程應(yīng)用均表明多變量混沌時(shí)間序列預(yù)測方法具有良好的預(yù)測效果.

由上，本文首先對(duì)運(yùn)行數(shù)據(jù)時(shí)間序列進(jìn)行混沌識(shí)別和相空間重構(gòu)；然后，采用基于奇異值分解（Singular value decomposition, SVD）的主成分分析（Principal component analysis, PCA）法對(duì)相空間進(jìn)行降維處理；進(jìn)而，構(gòu)建極限學(xué)習(xí)機(jī)（Extreme learning machine, ELM）、RBF神經(jīng)網(wǎng)絡(luò)、回聲狀態(tài)網(wǎng)絡(luò)（Echo state network, ESN）和Elman神經(jīng)網(wǎng)絡(luò)4種風(fēng)險(xiǎn)預(yù)測模型，采取迭代方式進(jìn)行風(fēng)險(xiǎn)預(yù)測；最后計(jì)算并分析不同模型的預(yù)測精度，判斷其實(shí)踐操作價(jià)值，旨在構(gòu)建一種精度較高的航班運(yùn)行風(fēng)險(xiǎn)預(yù)測方法，

1 航班運(yùn)行風(fēng)險(xiǎn)時(shí)間序列及混沌特征識(shí)別

1.1 航班運(yùn)行風(fēng)險(xiǎn)

“風(fēng)險(xiǎn)”一詞最早僅代表客觀的危險(xiǎn)，狹義風(fēng)險(xiǎn)是指損失的不確定性，如工程建筑風(fēng)險(xiǎn)關(guān)注人員傷亡和設(shè)備損失；而廣義風(fēng)險(xiǎn)強(qiáng)調(diào)表現(xiàn)的不確定性，如經(jīng)營風(fēng)險(xiǎn)關(guān)注未來市場走向，如金融風(fēng)險(xiǎn)關(guān)注收益變化.“風(fēng)險(xiǎn)”最受認(rèn)可的定義是由Coleman與Marks于1999年提出的，風(fēng)險(xiǎn)是結(jié)果嚴(yán)重程度和發(fā)生概率的綜合[22].國際民航組織（ICAO）附件19中，延用了此定義，風(fēng)險(xiǎn)為某一危險(xiǎn)源預(yù)計(jì)產(chǎn)生的可能性與嚴(yán)重性相乘，并以風(fēng)險(xiǎn)矩陣劃分等級(jí).如表1，綠區(qū)為可接受風(fēng)險(xiǎn)，即航班可以正常執(zhí)行；紅區(qū)為不可接受風(fēng)險(xiǎn)，即航班不可執(zhí)行；而黃區(qū)需要進(jìn)行有效風(fēng)險(xiǎn)緩解后，航班才可執(zhí)行.

表1 ICAO風(fēng)險(xiǎn)矩陣與中國民航局（CAAC）風(fēng)險(xiǎn)值Table 1 ICAO risk matrix and CAAC risk value

但是在使用風(fēng)險(xiǎn)矩陣評(píng)判航班風(fēng)險(xiǎn)等級(jí)時(shí)，存在兩個(gè)問題：①對(duì)于中國民航，任一人員傷亡、飛機(jī)損傷報(bào)廢都無法接受，風(fēng)險(xiǎn)嚴(yán)重程度無法如金融行業(yè)等以金錢數(shù)字衡量；②根據(jù)《2019年民航機(jī)場生產(chǎn)統(tǒng)計(jì)公報(bào)》，我國持續(xù)安全記錄已達(dá)112個(gè)月、8068×104h，而2019年航班已達(dá)1166萬架次，事故、事故癥候、不安全事件等發(fā)生概率接近極小值.

因此，為了解決風(fēng)險(xiǎn)矩陣的半定量化，進(jìn)而發(fā)展風(fēng)險(xiǎn)評(píng)估與預(yù)測技術(shù)，民航局于2015年頒布《航空承運(yùn)人運(yùn)行控制風(fēng)險(xiǎn)管控系統(tǒng)實(shí)施指南（AC-121-FS-2015-125）》以政府指導(dǎo)的方式量化了風(fēng)險(xiǎn)，經(jīng)實(shí)踐多年已驗(yàn)證有效.其中，將風(fēng)險(xiǎn)矩陣中1E～5A轉(zhuǎn)化為數(shù)字1到10，即風(fēng)險(xiǎn)數(shù)值化，如表1.該體系具有一定科學(xué)性，應(yīng)用于數(shù)值化評(píng)估與預(yù)測時(shí)存在兩個(gè)問題：①尚有少量的半定性/半定量和定性指標(biāo)，需要依賴專家判斷；②部分指標(biāo)不具備可預(yù)測性，如人因和航空器類指標(biāo)是不含時(shí)間維度特性的.

1.2 運(yùn)行數(shù)據(jù)選取

由上，延續(xù)AC-121-FS-2015-125中對(duì)風(fēng)險(xiǎn)評(píng)定方式，使用可定量的數(shù)據(jù)來源.提取國內(nèi)某大型航空公司風(fēng)險(xiǎn)管控系統(tǒng)和航空安全系統(tǒng)中風(fēng)險(xiǎn)源指標(biāo)，統(tǒng)計(jì)2016—2018年連續(xù)以天為單位的航班數(shù)據(jù)，形成航班運(yùn)行風(fēng)險(xiǎn)數(shù)據(jù)項(xiàng)1096組，數(shù)據(jù)類型和數(shù)據(jù)樣本見表2和表3.

首先，處理15個(gè)風(fēng)險(xiǎn)項(xiàng)的數(shù)據(jù)，使用SPSS軟件，取置信區(qū)間95%，使用SNK法和Dunnet-t法對(duì)數(shù)據(jù)列做多個(gè)均數(shù)多重比較，檢驗(yàn)結(jié)果概率均小于0.05，說明數(shù)據(jù)間具有獨(dú)立性；然后，再經(jīng)飛行、運(yùn)控、機(jī)務(wù)、安監(jiān)等聯(lián)合專家組檢驗(yàn)，確認(rèn)數(shù)據(jù)可信，故作為后續(xù)計(jì)算的標(biāo)準(zhǔn)樣本.

1.3 混沌識(shí)別

根據(jù)航班運(yùn)行風(fēng)險(xiǎn)數(shù)據(jù)樣本，分別構(gòu)建A1～A15和總風(fēng)險(xiǎn)時(shí)間序列，其中總風(fēng)險(xiǎn)時(shí)間序列如圖1.采用Wolf法計(jì)算，結(jié)果表明總風(fēng)險(xiǎn)時(shí)間序列具有混沌特征，同理對(duì)15個(gè)風(fēng)險(xiǎn)源時(shí)間序列進(jìn)行混沌識(shí)別，計(jì)算結(jié)果見表4.

從表4可見，16個(gè)風(fēng)險(xiǎn)序列均滿足 λ＞0，具有混沌特征.因此，首先通過對(duì)總風(fēng)險(xiǎn)時(shí)間序列進(jìn)行單變量混沌預(yù)測，得到未來1 d的相對(duì)誤差均值（Mean absolute percentage error, MAPE）為 21.33%.發(fā)現(xiàn)單一序列輸入，不足以構(gòu)建精準(zhǔn)模型.故后續(xù)采用多變量時(shí)間序列增加輸入，通過15個(gè)風(fēng)險(xiǎn)源序列預(yù)測總風(fēng)險(xiǎn).

表2 航班運(yùn)行風(fēng)險(xiǎn)統(tǒng)計(jì)數(shù)據(jù)Table 2 Risk assessment statistical data

表3 航班運(yùn)行風(fēng)險(xiǎn)時(shí)間序列局部數(shù)據(jù)樣本Table 3 Time series sample data of flight operation risk

圖1 航班運(yùn)行總風(fēng)險(xiǎn)時(shí)間序列Fig.1 Time series example of flight operation total risk

2 多變量混沌預(yù)測模型構(gòu)建

2.1 多變量時(shí)間序列相空間重構(gòu)

多變量時(shí)間序列相空間重構(gòu)由單變量相空間重構(gòu)發(fā)展而成.以15個(gè)風(fēng)險(xiǎn)序列的多變量相空間重構(gòu)結(jié)果作為預(yù)測模型的輸入，可使模型輸入包含各個(gè)風(fēng)險(xiǎn)序列的歷史信息.

對(duì)于長度為N的n個(gè)時(shí)間序列，j=1,2,···,n，時(shí)間延遲為τj和嵌入維為mj，運(yùn)用坐標(biāo)延遲重構(gòu)法如式（1），則多變量時(shí)間序列重構(gòu)后的相空間如式（2）：

式中，ti=max{(mj-1)*τj}+1,···,N-1,N.

τj和mj選擇是相空間重構(gòu)的關(guān)鍵.根據(jù)Takens定理，設(shè)d為混沌吸引子的維數(shù)，若m≥2d+1，則可在m維空間中重構(gòu)出與風(fēng)險(xiǎn)序列混沌吸引子保持微分同胚的軌線[23].采用C-C方法計(jì)算風(fēng)險(xiǎn)序列的最佳τj和mj.以風(fēng)險(xiǎn)時(shí)間序列Xj(t)為例，計(jì)算過程如下：

（1）設(shè)N是風(fēng)險(xiǎn)序列數(shù)據(jù)的個(gè)數(shù)，σ是標(biāo)準(zhǔn)差，并定義為：

表4 時(shí)間序列的時(shí)間延遲、嵌入維和最大Lyapunov指數(shù)Table 4 Time delay, embedding dimension, and maximum Lyapunov exponent of the time series

式中，S(m,re,t)=C(m,re,t)-Cm(1,re,t)，ΔS(m,t)=max{S(m,re,t)-min{S(m,re,t)}.

τj和mj的合理選擇有助于對(duì)混沌系統(tǒng)的擬合和預(yù)測，表4為15個(gè)風(fēng)險(xiǎn)時(shí)間序列的最佳延遲時(shí)間和嵌入維，據(jù)此可進(jìn)行多變量時(shí)間序列相空間重構(gòu)，重構(gòu)后的維度為各個(gè)序列的總嵌入維

2.2 基于 PCA的相空間降維方法

由于航班運(yùn)行風(fēng)險(xiǎn)各指標(biāo)變量均代表不同含義、混沌重構(gòu)過程的τj和mj也不盡相同，所以在多變量相空間重構(gòu)后，不同變量在不同歷史時(shí)期的觀測數(shù)據(jù)混合在一起，可能造成了信息冗余或噪聲干擾，產(chǎn)生預(yù)測誤差.采用PCA抑制此問題.PCA基本思想是將一組具有一定相關(guān)性的變量重組為互不相關(guān)的主成分，通過保留方差貢獻(xiàn)率較大的主成分代替原變量，實(shí)現(xiàn)數(shù)據(jù)降維.采用PCA改進(jìn)航班運(yùn)行風(fēng)險(xiǎn)的預(yù)測過程，可以保證：①預(yù)測模型的不同輸入成分之間相互獨(dú)立、信息不重疊；②輸入的狀態(tài)空間在保證盡量多信息的情況下實(shí)現(xiàn)維度降低.

PCA常用的主成分提取方法有特征值分解（Eigen value decomposition, EVD）和奇異值分解（SVD），相對(duì)于EVD，SVD在數(shù)據(jù)重構(gòu)上誤差較小，故本文采用基于SVD的PCA進(jìn)行降維處理，設(shè)ξ1≥ ξ2≥ ···≥ ξq＞ 0為原相空間Y的奇異值，通過SVD進(jìn)行分解如式（7）：

式中，Q=[q1,q2,···,qN-τw]為Y的左 奇異矩陣，為Y的右奇異矩陣，D=diag(ξ1,ξ2,···,ξd).

假設(shè)前 γ個(gè)主成分的累計(jì)方差貢獻(xiàn)率為Gγ，見式（8），一般認(rèn)為，當(dāng)滿足Gγ≥85%時(shí)[24]，可將γ個(gè)主成分構(gòu)成的狀態(tài)空間，見式（9），代替原相空間，實(shí)現(xiàn)降維.

2.3 預(yù)測模型

混沌時(shí)間序列預(yù)測理論中常用方法有：全域法、局域法、基于Lyapunov指數(shù)預(yù)測法、基于Volterra自適應(yīng)濾波器預(yù)測法、基于神經(jīng)網(wǎng)絡(luò)預(yù)測法等.此處在選取方法時(shí)，一方面總結(jié)以往研究，發(fā)現(xiàn)神經(jīng)網(wǎng)絡(luò)類方法在混沌預(yù)測中表現(xiàn)較為突出；另一方面，通過查閱文獻(xiàn)，篩選出預(yù)測結(jié)果最好的幾種方法，包含了極限學(xué)習(xí)機(jī)[18]、RBF神經(jīng)網(wǎng)絡(luò)[25]、回聲狀態(tài)網(wǎng)絡(luò)[26]和Elman神經(jīng)網(wǎng)絡(luò)[21]等.而后，通過原理分析發(fā)現(xiàn)，上述4個(gè)模型均由傳統(tǒng)神經(jīng)網(wǎng)絡(luò)發(fā)展而來.兩方面相互印證，決定采用該4種方法，后續(xù)分別構(gòu)建短期預(yù)測模型，通過誤差分析，提出具有較好適應(yīng)性的預(yù)測方案.

對(duì)于給定的N個(gè)不同樣本{(u(i),O(i)),i=1,2,···N}，其中u(i)=[u1(i),u2(i),···,ul(i)]T∈Rl為輸入層，O(i)=[O1(i),O2(i),···,Ov(i)]T∈Rv為輸出層.統(tǒng)一設(shè)定隱含層神經(jīng)元個(gè)數(shù)為h，隱含層神經(jīng)元與輸入層的連接權(quán)值為、與輸出層的連接權(quán)值為，其中j=1,2,···,h.不同模型的原理如下：

（1）極限學(xué)習(xí)機(jī)ELM.

ELM是單隱層前饋神經(jīng)網(wǎng)絡(luò)，如圖2.設(shè)ELM隱含層神經(jīng)元的閾值為b，訓(xùn)練過程中，win和b隨機(jī)產(chǎn)生，并保持不變，加快了神經(jīng)網(wǎng)絡(luò)的訓(xùn)練速度并使之具有較好的泛化能力.

圖2 極限學(xué)習(xí)機(jī)結(jié)構(gòu)Fig.2 ELM structure

圖2中，從輸入到輸出的函數(shù)表達(dá)式為：

式中，i=1,2,···N；只需設(shè)定隱含層神經(jīng)元個(gè)數(shù)h，并確定激活函數(shù)，并通過求解式(11)即可得出wout，其中，激活函數(shù)多選用sigmoid函數(shù).

式中，H?為H的偽逆，H為隱含層輸出矩陣，即

（2）RBF 神經(jīng)網(wǎng)絡(luò).

RBF是從三層前向型網(wǎng)絡(luò)發(fā)展而來，不同之處在于其隱含層的神經(jīng)元變換函數(shù)使用徑向基函數(shù).RBF神經(jīng)網(wǎng)絡(luò)的優(yōu)點(diǎn)在于輸入層和隱含層之間不需通過權(quán)值連接，當(dāng)RBF中心點(diǎn)確定以后，即可確定二者的映射關(guān)系，從而使網(wǎng)絡(luò)結(jié)構(gòu)簡單、訓(xùn)練簡潔.

圖3中，{Rj,j=1,2,···,h}為激活函數(shù)，常使用高斯函數(shù)如式（12），由此，從輸入到輸出的函數(shù)表達(dá)式為：

圖3 RBF神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.3 RBF neural network structure

（3）回聲狀態(tài)網(wǎng)絡(luò).

回聲狀態(tài)網(wǎng)絡(luò)（ESN）是一種遞歸神經(jīng)網(wǎng)絡(luò)，結(jié)構(gòu)如圖4，與傳統(tǒng)神經(jīng)網(wǎng)絡(luò)相比，ESN使用大規(guī)模隨機(jī)連接的遞歸網(wǎng)絡(luò)形成一個(gè)儲(chǔ)備池，取代隱含層，通過預(yù)先設(shè)定儲(chǔ)備池權(quán)值矩陣的譜半徑可保證網(wǎng)絡(luò)的穩(wěn)定性.

ESN從輸入到輸出的表達(dá)式為：

圖4 回聲狀態(tài)網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.4 ESN structure

式中，φ(i)為儲(chǔ)備池的狀態(tài)向量，WDR為儲(chǔ)備池內(nèi)部連接權(quán)值矩陣；fin為激活函數(shù)，常取雙曲正切函數(shù)；fout為輸出函數(shù)，常取恒等函數(shù)，此時(shí)wout可通過下式求解：

式 中，B?為B={u(β+1),u(β+2),···,u(N)}的偽 逆，{u(1),u(2),···,u(β)}為舍棄的一段初始數(shù)據(jù).

（4）Elman神經(jīng)網(wǎng)絡(luò).

Elman是一種典型的局部回歸網(wǎng)絡(luò)，屬于反饋神經(jīng)網(wǎng)絡(luò)，如圖5，它與前向神經(jīng)網(wǎng)絡(luò)不同之處在于隱含層中增加了一個(gè)承接層，構(gòu)成局部的反饋.承接層用于接受隱含層前一時(shí)刻的輸出信號(hào)，并作為下一時(shí)刻輸入信號(hào)的一部分重新輸入到隱含層中，形成一個(gè)一步延時(shí)算子，使網(wǎng)絡(luò)具有動(dòng)態(tài)記憶功能.由于承接層具有可以記憶過去狀態(tài)的特性，非常適合時(shí)間序列預(yù)測問題.

圖5 Elman神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.5 Elman neural network structure

Elman從輸入到輸出的表達(dá)式為：

式中，χ(i)為隱含層神經(jīng)元向量，χc(i)為承接層的反饋向量，Wc為承接層到隱含層的連接權(quán)值，f為隱含層神經(jīng)元的傳遞函數(shù)，fout為輸出函數(shù).

2.4 預(yù)測流程

根據(jù)A1～A15風(fēng)險(xiǎn)時(shí)間序列，對(duì)未來pd的總風(fēng)險(xiǎn)值進(jìn)行短期預(yù)測，預(yù)測流程如圖6：

圖6 預(yù)測流程Fig.6 Prediction process

1）根據(jù)風(fēng)險(xiǎn)項(xiàng)A1～A15的混沌特征分析結(jié)果，確定延遲時(shí)間和嵌入維，進(jìn)行多變量相空間重構(gòu)；

2）對(duì)相空間進(jìn)行PCA分析，確定前若干個(gè)主成分構(gòu)成的狀態(tài)空間，實(shí)現(xiàn)相空間的降維；

3）根據(jù)迭代預(yù)測方法，分別構(gòu)建ELM、RBF、ESN和Elman風(fēng)險(xiǎn)短期預(yù)測.假設(shè)當(dāng)前時(shí)刻為ti，在迭代預(yù)測過程，首先以第ti天降維后的狀態(tài)空間作為預(yù)測模型的輸入，預(yù)測出第ti+1天A1～A15的風(fēng)險(xiǎn)值；然后，利用第ti+1天的預(yù)測值，重新進(jìn)行多變量相空間重構(gòu)并作降維處理，再將此時(shí)降維后的狀態(tài)空間輸入至預(yù)測模型中，預(yù)測出第ti+2天A1～A15的風(fēng)險(xiǎn)值.以此類推，直至完成第ti+p天的預(yù)測，得到A1～A15的pd預(yù)測結(jié)果.

4）采用支持向量機(jī)（SVM）進(jìn)行風(fēng)險(xiǎn)評(píng)估，得到總風(fēng)險(xiǎn)值的pd預(yù)測值，通過與總風(fēng)險(xiǎn)的真實(shí)值比較，計(jì)算預(yù)測誤差并分析不同模型的預(yù)測精度.

3 實(shí)例分析

3.1 主成分降維

將15個(gè)風(fēng)險(xiǎn)源時(shí)間序列數(shù)據(jù)作為預(yù)測樣本，其中，將前900 d作為訓(xùn)練集，剩余196 d作為測試集.經(jīng)過多變量相空間重構(gòu)后再使用PCA降維，各主成分的方差貢獻(xiàn)率如表5.

根據(jù)PCA原理，PCA處理后得到的主成分應(yīng)能解釋原始變量盡量多的信息，累計(jì)方差貢獻(xiàn)率達(dá)到85%以上為宜；同時(shí)，當(dāng)變量之間相關(guān)性較強(qiáng)的時(shí)候，少數(shù)幾個(gè)主成分的累積貢獻(xiàn)率就能達(dá)到較高水平，當(dāng)變量之間相關(guān)性較弱的時(shí)候，則需較多的主成分才能達(dá)到較高的累積貢獻(xiàn)率，但更為重要的是保證計(jì)算后達(dá)到有效的降維效果.因此，在較高的累計(jì)方差貢獻(xiàn)率基礎(chǔ)上結(jié)合碎石圖拐點(diǎn)出現(xiàn)位置，再綜合選擇.通過采用85%、90%和95%不同閾值，經(jīng)上文預(yù)測流程計(jì)算得到，當(dāng)保留前31個(gè)主成分累計(jì)方差貢獻(xiàn)率為90%時(shí)，此處預(yù)測結(jié)果最為理想，即維數(shù)由62維壓縮至31維.

3.2 參數(shù)優(yōu)化

根據(jù)訓(xùn)練集數(shù)據(jù)，采用交叉驗(yàn)證法、遺傳算法和粒子群算法對(duì)預(yù)測模型的參數(shù)進(jìn)行優(yōu)化.但由于后續(xù)預(yù)測過程采用多步迭代方式，過程復(fù)雜且易受預(yù)測步長的影響，因此假設(shè)多步迭代預(yù)測最優(yōu)的前提為單步預(yù)測結(jié)果最優(yōu)，綜合衡量不同參數(shù)優(yōu)化方法和多次試驗(yàn)的輸出，參數(shù)優(yōu)化結(jié)果如表6所示.

表5 主成分分析的方差貢獻(xiàn)率Table 5 Partial variance contribution rate of PCA

表6 預(yù)測模型的參數(shù)優(yōu)化結(jié)果Table 6 Parameter optimization results of the prediction models

3.3 短期預(yù)測的測試結(jié)果

結(jié)合航空公司實(shí)際需要，設(shè)定預(yù)測范圍為未來7 d，以第900天為起始點(diǎn)，向后預(yù)測得到第901～907天風(fēng)險(xiǎn)值，如圖7.但僅從一組數(shù)據(jù)無法說明模型預(yù)測精度，再從第901天至第1090天依次預(yù)測，總計(jì)190次，每次均得到7 d預(yù)測結(jié)果.統(tǒng)計(jì)所有預(yù)測的第1、3、5天結(jié)果如圖8；其中，圖8（a）均為第 1 天預(yù)測結(jié)果，其范圍為 [901，1090]；圖8（b）是第 3 天結(jié)果，范圍為 [903，1092]；圖8（c）是第 5 天結(jié)果，范圍為 [905，1094]，圖8（a）（b）（c）均包含190個(gè)數(shù)據(jù).

圖7 從第900天向后預(yù)測樣例Fig.7 Example of prediction results from the 900th day

3.4 結(jié)果分析

計(jì)算ELM、RBF、ESN和Elman等不同模型短期預(yù)測結(jié)果的相對(duì)誤差.統(tǒng)計(jì)相對(duì)誤差（Relative error, RE）的頻數(shù)分布，并以降維前的預(yù)測作對(duì)比，如表7.

從表8可見，RBF降維后的預(yù)測效果最佳，對(duì)于RE＜25%的預(yù)測結(jié)果，未來第1天占比達(dá)到82.62%，第3天為78.95%，第5天仍高于75%，隨預(yù)測范圍的增加，預(yù)測精度逐步遞減.4種預(yù)測模型在降維后，第1天預(yù)測結(jié)果RE＜50%的出現(xiàn)頻數(shù)平均為90.65%，第3天為 86.18%，第5天為 85.26%，表明RE＞50%的異常結(jié)果出現(xiàn)次數(shù)較少，4種預(yù)測模型均可行有效.進(jìn)一步對(duì)比降維前后結(jié)果發(fā)現(xiàn)，降維后RE＜50%占比相對(duì)于降維前提高了約6%，說明降維處理有助于提高RBF的預(yù)測精度.

圖8 預(yù)測結(jié)果示例.（a）第 1 天；（b）第 3 天；（c）第 5 天Fig.8 Example of prediction results: (a) day 1; (b) day 3; (c) day 5

3.5 誤差分析

從表7可見，不同模型預(yù)測相對(duì)誤差的頻數(shù)分布均呈右偏，若以均值反映相對(duì)誤差的集中趨勢，結(jié)果會(huì)偏大.因此，為使相對(duì)誤差均值（MAPE）更具代表性，采用修正的MAPE值，即去掉最大5%和最小5%的數(shù)據(jù)后所求得的MAPE值，結(jié)果見圖9～10.

從圖9可見，當(dāng)訓(xùn)練樣本為300時(shí)，誤差較大；隨著訓(xùn)練樣本增加，誤差開始下降；當(dāng)訓(xùn)練樣本從300增加至500時(shí)，7 d內(nèi)的修正MAPE值平均下降5.48%，預(yù)測效果明顯提升；當(dāng)訓(xùn)練樣本從700增加至900時(shí)，修正MAPE值平均下降1.27%，效果提升開始放緩；當(dāng)訓(xùn)練樣本從800增加至900時(shí)，修正MAPE值平均僅下降0.21%，可視為無明顯變化；因此，基于現(xiàn)有數(shù)據(jù)情況，后續(xù)使用訓(xùn)練樣本為900進(jìn)行計(jì)算.

圖10是基于降維后的RBF預(yù)測結(jié)果，修正MAPE值隨預(yù)測周期增長而逐步遞增，從第1天至第7天平均每天增加2.35%，呈現(xiàn)誤差累積現(xiàn)象.第1天修正MAPE值僅為11.32%，精度較高；第5天MAPE值仍為18.21%，其預(yù)測精度和預(yù)測周期均優(yōu)于文獻(xiàn)[12]和[13].根據(jù)國內(nèi)航空公司使用要求，當(dāng)預(yù)測相對(duì)誤差低于20%時(shí)，結(jié)果具有實(shí)踐操作價(jià)值，因此該方法所得未來5 d預(yù)測結(jié)果均可作為依據(jù)進(jìn)而指導(dǎo)生產(chǎn)實(shí)踐.

表7 各模型降維前后預(yù)測相對(duì)誤差頻數(shù)Table 7 RE frequency of each prediction model before and after dimension reduction

圖9 訓(xùn)練樣本數(shù)量與預(yù)測精度Fig.9 Number of training samples and prediction accuracy

圖10 RBF降維后的多變量預(yù)測精度Fig.10 Prediction accuracy of the RBF model after dimension reduction

4 結(jié)論

通過混沌識(shí)別和相空間重構(gòu)，構(gòu)建的4種風(fēng)險(xiǎn)預(yù)測模型在迭代預(yù)測后發(fā)現(xiàn)：在降維前后的ELM、RBF、ESN和Elman預(yù)測模型中，降維后的RBF預(yù)測模型效果最佳；其未來第1天預(yù)測結(jié)果RE＜25%的頻數(shù)可達(dá)到82.62%，修正MAPE值僅為11.32%；至第5天MAPE值仍為18.21%，說明未來5 d的預(yù)測結(jié)果均可滿足航空公司實(shí)踐操作要求.從預(yù)測精度和周期兩個(gè)角度綜合評(píng)判，多變量時(shí)間序列方法對(duì)航班運(yùn)行風(fēng)險(xiǎn)的預(yù)測可行且有效.

國內(nèi)航空公司以周為單位制定航班計(jì)劃，因此可在文中結(jié)果基礎(chǔ)上，最早于5 d前從飛機(jī)調(diào)配和人員調(diào)度等方面調(diào)整計(jì)劃編排，提前降低運(yùn)行風(fēng)險(xiǎn).該方案還可擴(kuò)展至其他領(lǐng)域，適用于結(jié)果受多變量綜合影響的、各變量為混沌時(shí)間序列的短期預(yù)測問題，其預(yù)測周期隨多變量的混沌時(shí)間序列特性而定.