基于多變量混沌時間序列的航班運行風險預(yù)測模型

王巖韜，李景良，谷潤平

中國民航大學空管學院，天津 300300

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

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

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

由上，本文首先對運行數(shù)據(jù)時間序列進行混沌識別和相空間重構(gòu)；然后，采用基于奇異值分解（Singular value decomposition, SVD）的主成分分析（Principal component analysis, PCA）法對相空間進行降維處理；進而，構(gòu)建極限學習機（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種風險預(yù)測模型，采取迭代方式進行風險預(yù)測；最后計算并分析不同模型的預(yù)測精度，判斷其實踐操作價值，旨在構(gòu)建一種精度較高的航班運行風險預(yù)測方法，

1 航班運行風險時間序列及混沌特征識別

1.1 航班運行風險

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

表1 ICAO風險矩陣與中國民航局（CAAC）風險值Table 1 ICAO risk matrix and CAAC risk value

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

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

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

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

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

1.3 混沌識別

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

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

表2 航班運行風險統(tǒng)計數(shù)據(jù)Table 2 Risk assessment statistical data

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

圖1 航班運行總風險時間序列Fig.1 Time series example of flight operation total risk

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

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

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

對于長度為N的n個時間序列，j=1,2,···,n，時間延遲為τj和嵌入維為mj，運用坐標延遲重構(gòu)法如式（1），則多變量時間序列重構(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)出與風險序列混沌吸引子保持微分同胚的軌線[23].采用C-C方法計算風險序列的最佳τj和mj.以風險時間序列Xj(t)為例，計算過程如下：

（1）設(shè)N是風險序列數(shù)據(jù)的個數(shù)，σ是標準差，并定義為：

表4 時間序列的時間延遲、嵌入維和最大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的合理選擇有助于對混沌系統(tǒng)的擬合和預(yù)測，表4為15個風險時間序列的最佳延遲時間和嵌入維，據(jù)此可進行多變量時間序列相空間重構(gòu)，重構(gòu)后的維度為各個序列的總嵌入維

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

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

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

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

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

2.3 預(yù)測模型

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

對于給定的N個不同樣本{(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)元個數(shù)為h，隱含層神經(jīng)元與輸入層的連接權(quán)值為、與輸出層的連接權(quán)值為，其中j=1,2,···,h.不同模型的原理如下：

（1）極限學習機ELM.

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

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

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

式中，i=1,2,···N；只需設(shè)定隱含層神經(jīng)元個數(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)點在于輸入層和隱含層之間不需通過權(quán)值連接，當RBF中心點確定以后，即可確定二者的映射關(guān)系，從而使網(wǎng)絡(luò)結(jié)構(gòu)簡單、訓練簡潔.

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

圖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ī)模隨機連接的遞歸網(wǎng)絡(luò)形成一個儲備池，取代隱含層，通過預(yù)先設(shè)定儲備池權(quán)值矩陣的譜半徑可保證網(wǎng)絡(luò)的穩(wěn)定性.

ESN從輸入到輸出的表達式為：

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

式中，φ(i)為儲備池的狀態(tài)向量，WDR為儲備池內(nèi)部連接權(quán)值矩陣；fin為激活函數(shù)，常取雙曲正切函數(shù)；fout為輸出函數(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òu)成局部的反饋.承接層用于接受隱含層前一時刻的輸出信號，并作為下一時刻輸入信號的一部分重新輸入到隱含層中，形成一個一步延時算子，使網(wǎng)絡(luò)具有動態(tài)記憶功能.由于承接層具有可以記憶過去狀態(tài)的特性，非常適合時間序列預(yù)測問題.

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

Elman從輸入到輸出的表達式為：

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

2.4 預(yù)測流程

根據(jù)A1～A15風險時間序列，對未來pd的總風險值進行短期預(yù)測，預(yù)測流程如圖6：

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

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

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

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

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

3 實例分析

3.1 主成分降維

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

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

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

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

表5 主成分分析的方差貢獻率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è)定預(yù)測范圍為未來7 d，以第900天為起始點，向后預(yù)測得到第901～907天風險值，如圖7.但僅從一組數(shù)據(jù)無法說明模型預(yù)測精度，再從第901天至第1090天依次預(yù)測，總計190次，每次均得到7 d預(yù)測結(jié)果.統(tǒng)計所有預(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個數(shù)據(jù).

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

3.4 結(jié)果分析

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

從表8可見，RBF降維后的預(yù)測效果最佳，對于RE＜25%的預(yù)測結(jié)果，未來第1天占比達到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ù)測模型均可行有效.進一步對比降維前后結(jié)果發(fā)現(xiàn)，降維后RE＜50%占比相對于降維前提高了約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ù)測相對誤差的頻數(shù)分布均呈右偏，若以均值反映相對誤差的集中趨勢，結(jié)果會偏大.因此，為使相對誤差均值（MAPE）更具代表性，采用修正的MAPE值，即去掉最大5%和最小5%的數(shù)據(jù)后所求得的MAPE值，結(jié)果見圖9～10.

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

圖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)于文獻[12]和[13].根據(jù)國內(nèi)航空公司使用要求，當預(yù)測相對誤差低于20%時，結(jié)果具有實踐操作價值，因此該方法所得未來5 d預(yù)測結(jié)果均可作為依據(jù)進而指導(dǎo)生產(chǎn)實踐.

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

圖9 訓練樣本數(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é)論

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

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