王巖韜，李景良，谷潤(rùn)平

(中國(guó)民航大學(xué) 空中交通管理學(xué)院，天津300300)

0 引 言

民航工作高質(zhì)量發(fā)展依賴于航班運(yùn)行風(fēng)險(xiǎn)防控能力的提升.歐美最早從1991年開始航班風(fēng)險(xiǎn)評(píng)價(jià)研究[1].研究理論和手段側(cè)重于事故樹等傳統(tǒng)方法[2-3].實(shí)踐方面，2007—2012年美國(guó)、加拿大民航局研發(fā)了Flight Risk Assessment Tool(FRAT)[4]、Recheck[5]等飛行風(fēng)險(xiǎn)評(píng)估工具，至今鮮見(jiàn)可信的使用報(bào)告和精度分析.國(guó)內(nèi)，航班風(fēng)險(xiǎn)評(píng)估研究最早于1999年[6]，王巖韜等建立了一套航班運(yùn)行風(fēng)險(xiǎn)評(píng)估體系[7]，運(yùn)用多算法協(xié)作將風(fēng)險(xiǎn)分類正確率提高至95%[8].對(duì)于民航非動(dòng)態(tài)風(fēng)險(xiǎn)評(píng)估，上述研究雖已取得較好成果，但將其應(yīng)用于民航一線后發(fā)現(xiàn)，評(píng)估技術(shù)只能做到隨航班要素變化而變化，只是對(duì)已有狀態(tài)評(píng)估，而風(fēng)險(xiǎn)管控需要時(shí)間提前量，對(duì)此風(fēng)險(xiǎn)預(yù)測(cè)技術(shù)成為風(fēng)險(xiǎn)防控實(shí)質(zhì)所亟需.

國(guó)外最早于2008年研究航空相關(guān)風(fēng)險(xiǎn)預(yù)測(cè)[9].2012年，歐盟啟動(dòng)航空運(yùn)行主動(dòng)安全管理項(xiàng)目PROSPERO(Proactive Safety Performance for Operations)預(yù)測(cè)航空系統(tǒng)風(fēng)險(xiǎn)，但未見(jiàn)詳細(xì)方案和結(jié)果報(bào)告，最新消息表明該項(xiàng)目已于2015年戛然而止[10].Andrej Lalis 等[11]使用線性回歸模型和ARMA模型對(duì)航空安全績(jī)效進(jìn)行預(yù)測(cè)，但以月為單位的寬幅時(shí)間序列嚴(yán)重削弱了其實(shí)用價(jià)值.Zhang等[12]融合支持向量機(jī)和深度神經(jīng)網(wǎng)絡(luò)算法預(yù)測(cè)航空事件嚴(yán)重程度，但其實(shí)質(zhì)過(guò)程與文獻(xiàn)[7]和[8]類似，對(duì)未來(lái)預(yù)測(cè)能力不足.王巖韜團(tuán)隊(duì)[13]于2019年提出基于動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)的航班運(yùn)行風(fēng)險(xiǎn)預(yù)測(cè)模型，主要針對(duì)單個(gè)航班的進(jìn)近、著陸階段進(jìn)行風(fēng)險(xiǎn)預(yù)測(cè)，但無(wú)法體現(xiàn)航空公司整體運(yùn)行狀況.航班風(fēng)險(xiǎn)預(yù)測(cè)的專項(xiàng)研究很少，無(wú)論是小時(shí)、日、月等不同時(shí)間維度的預(yù)測(cè)方法，還是預(yù)測(cè)精度穩(wěn)定性分析，亦或是預(yù)測(cè)精度提升技術(shù)，均缺少深入的研究，是民航安全領(lǐng)域尚未解決的關(guān)鍵問(wèn)題之一.

本文旨在探索一種實(shí)用有效的航班運(yùn)行風(fēng)險(xiǎn)短期預(yù)測(cè)方法.借鑒混沌時(shí)間序列預(yù)測(cè)理論，對(duì)數(shù)據(jù)進(jìn)行相空間重構(gòu)，實(shí)現(xiàn)混沌識(shí)別；基于集成經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition,EEMD)進(jìn)行閾值降噪處理，采用極限學(xué)習(xí)機(jī)(Extreme Learning Machine，ELM)構(gòu)建改進(jìn)的混沌短期預(yù)測(cè)模型，并對(duì)預(yù)測(cè)效果進(jìn)行數(shù)值分析.

1 風(fēng)險(xiǎn)數(shù)據(jù)時(shí)間序列及其混沌識(shí)別

1.1 風(fēng)險(xiǎn)數(shù)據(jù)時(shí)間序列

航班運(yùn)行總風(fēng)險(xiǎn)是指當(dāng)日航班運(yùn)行相關(guān)危險(xiǎn)源導(dǎo)致的后果嚴(yán)重程度和發(fā)生概率的綜合.中國(guó)民航制定咨詢通告AC-121-FS-2015-125，使用數(shù)值1～10 量化評(píng)估了不同嚴(yán)重程度與發(fā)生概率組合.

表1 以天為單位，統(tǒng)計(jì)國(guó)內(nèi)某大型航空公司2016—2018年連續(xù)1 096 d 的航班狀況，提取風(fēng)險(xiǎn)管控系統(tǒng)和航空安全系統(tǒng)中所有差錯(cuò)指標(biāo)，形成航班運(yùn)行風(fēng)險(xiǎn)數(shù)據(jù).使用SPSS 軟件，取置信區(qū)間95%，使用SNK法和Dunnet-t法對(duì)數(shù)據(jù)列做多個(gè)均數(shù)多重比較，檢驗(yàn)結(jié)果概率均小于0.05，說(shuō)明數(shù)據(jù)間具有獨(dú)立性；再經(jīng)飛行、運(yùn)控、機(jī)務(wù)、安監(jiān)等聯(lián)合專家組檢驗(yàn)，確認(rèn)數(shù)據(jù)可信.以航班運(yùn)行總風(fēng)險(xiǎn)時(shí)間序列(簡(jiǎn)稱風(fēng)險(xiǎn)序列)作為后續(xù)計(jì)算的標(biāo)準(zhǔn)樣本，如圖1所示.

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

1.2 混沌識(shí)別

運(yùn)用坐標(biāo)延遲重構(gòu)法進(jìn)行相空間重構(gòu)，設(shè)某一時(shí)刻ti的風(fēng)險(xiǎn)值為x(ti)，ti=1+τw,2+τw,3+τw,…,N，i為時(shí)間索引，樣本量為N的風(fēng)險(xiǎn)序列為X(N)=[x(1),x(2),…,x(N)]，則ti時(shí)刻重構(gòu)后的相空間為

式中：τw為延遲時(shí)間窗口，τw=(m-1) ·τ；τ為時(shí)間延遲，m為嵌入維，τ和m取值是相空間重構(gòu)的關(guān)鍵，最佳τ和m采用C-C方法計(jì)算得出.

圖1 風(fēng)險(xiǎn)序列數(shù)據(jù)樣本Fig.1 Time series example of flight operation risk

重構(gòu)相空間后，采用Wolf 方法計(jì)算最大Lyapunov 指數(shù)，識(shí)別風(fēng)險(xiǎn)序列的混沌特征.Lyapunov指數(shù)用于描述相鄰相空間軌線隨時(shí)間推移而偏離擴(kuò)散程度，當(dāng)最大Lyapunov 指數(shù)λ>0，意味著軌線出現(xiàn)指數(shù)擴(kuò)散，認(rèn)為風(fēng)險(xiǎn)序列具有混沌特征.

2 混沌短期預(yù)測(cè)模型構(gòu)建

2.1 短期預(yù)測(cè)方法

短期預(yù)測(cè)方法分為直接預(yù)測(cè)和迭代預(yù)測(cè).直接預(yù)測(cè)通過(guò)映射f(·)得到ti+p時(shí)刻的預(yù)測(cè)值，即，其中，p為預(yù)測(cè)周期.迭代預(yù)測(cè)：先確定映射g(·)，進(jìn)行單步預(yù)測(cè)，得到ti+1 時(shí)刻預(yù)測(cè)值x(ti+1)=g[Y(ti)]；再根據(jù)x(ti+1) 重構(gòu)相空間Y(ti+1) ，重新執(zhí)行單步預(yù)測(cè)，得到預(yù)測(cè)值x(ti+2)；反復(fù)迭代p次，得到ti+p時(shí)刻預(yù)測(cè)值.直接預(yù)測(cè)不存在誤差累積效應(yīng)，但p較大時(shí)，難以構(gòu)建滿足精度要求的預(yù)測(cè)模型；迭代預(yù)測(cè)過(guò)程依賴于單步預(yù)測(cè)，誤差隨迭代次數(shù)增加而增大，當(dāng)輸入全為預(yù)測(cè)值時(shí)，誤差會(huì)急劇增大.

多次實(shí)驗(yàn)發(fā)現(xiàn)，采用迭代預(yù)測(cè)方法優(yōu)于直接預(yù)測(cè)，故后文結(jié)果均通過(guò)迭代預(yù)測(cè)得到.

2.2 ELM方法

ELM 是一種單隱層前饋神經(jīng)網(wǎng)絡(luò)，具有學(xué)習(xí)速度快，泛化性能好等特點(diǎn)，廣泛應(yīng)用于混沌時(shí)間序列預(yù)測(cè)[14].對(duì)于Q個(gè)樣本其中，j為樣本索引，Ij為輸入層，，k為輸入層神經(jīng)元個(gè)數(shù)；Oj為期望輸出值，，l為輸出層神經(jīng)元個(gè)數(shù).設(shè)隱含層神經(jīng)元個(gè)數(shù)為n，隱含層激活函數(shù)為g(·)，隱含層某一神經(jīng)元與輸入層、輸出層的連接權(quán)值分別為和，其中為第i(i=1,2,…,n)個(gè)隱含層神經(jīng)元與輸入層第e(e=1,2,…,k)個(gè)神經(jīng)元的連接權(quán)值為第i個(gè)隱含層神經(jīng)元與輸出層第u(u=1,2,…,l)個(gè)神經(jīng)元的連接權(quán)值.采用ELM逼近具有Q個(gè)樣本的目標(biāo)輸出時(shí)，可得

式中：H為隱含層輸出矩陣 ，bi為神經(jīng)網(wǎng)絡(luò)的閾值，i=1,2,…,n；w為隱含層與輸出層的連接權(quán)值，O為期望輸出值，.

在ELM 訓(xùn)練過(guò)程，隱含層與輸入層的連接權(quán)值win和閾值bi(i=1,2,…,n)隨機(jī)產(chǎn)生，并在訓(xùn)練過(guò)程保持不變；只需設(shè)置隱含層神經(jīng)元個(gè)數(shù)和激活函數(shù)，w求解公式為

其解為，H?為H的偽逆.

3 EEMD改進(jìn)的風(fēng)險(xiǎn)短期預(yù)測(cè)模型

3.1 EEMD方法

經(jīng)驗(yàn)?zāi)B(tài)分解 (Empirical Mode Decomposition, EMD)對(duì)序列逐級(jí)分解，產(chǎn)生多個(gè)具有不同特征尺度的本征模態(tài)函數(shù)IMF，即

使用集成經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition,EEMD)是為克服EMD 的模態(tài)混疊現(xiàn)象.EEMD 改進(jìn)之處為：原序列加入白噪聲序列后，進(jìn)行EMD 分解；反復(fù)加入不同的白噪聲序列，反復(fù)分解；當(dāng)重復(fù)次數(shù)足夠多時(shí)，對(duì)IMF和殘余項(xiàng)分別求平均，抵消噪聲，得到EEMD分解結(jié)果.

3.2 端點(diǎn)延拓

端點(diǎn)效應(yīng)是指分解過(guò)程中因序列兩端點(diǎn)不一定是極值點(diǎn)，導(dǎo)致上、下包絡(luò)線在序列兩端出現(xiàn)發(fā)散的現(xiàn)象，該現(xiàn)象會(huì)逐漸對(duì)內(nèi)部數(shù)據(jù)產(chǎn)生影響，使分解結(jié)果失真.對(duì)ti時(shí)刻的風(fēng)險(xiǎn)序列預(yù)測(cè)，模型輸入向量屬于序列X(ti)的端點(diǎn)數(shù)據(jù)，經(jīng)過(guò)計(jì)算發(fā)現(xiàn)，端點(diǎn)效應(yīng)對(duì)預(yù)測(cè)結(jié)果產(chǎn)生了明顯的負(fù)面影響.

采用基于預(yù)測(cè)值的端點(diǎn)延拓方法，避免端點(diǎn)效應(yīng).根據(jù)第2 節(jié)混沌短期預(yù)測(cè)方法，計(jì)算未來(lái)一定時(shí)間范圍的預(yù)測(cè)值，保證至少存在極大值點(diǎn)和極小值點(diǎn)；將原序列與預(yù)測(cè)序列連接，在分解前實(shí)現(xiàn)原序列的端點(diǎn)延拓.

3.3 閾值降噪

借鑒小波閾值降噪思路，選擇合適閾值和閾值函數(shù)截?cái)郔MF，把處理后的IMF 合成新序列，生成降噪結(jié)果.

式中：為中噪聲的方差，.

閾值函數(shù)分為軟閾值函數(shù)和硬閾值函數(shù).軟閾值函數(shù)會(huì)使降噪結(jié)果失真，硬閾值函數(shù)會(huì)使降噪結(jié)果不連續(xù).采用改進(jìn)的軟硬折中閾值函數(shù)[16]，即

式中：α為折中因子；為的降噪結(jié)果；sign(·) 為符號(hào)函數(shù).

3.4 改進(jìn)預(yù)測(cè)模型構(gòu)建

根據(jù)EEMD閾值降噪方法對(duì)風(fēng)險(xiǎn)短期預(yù)測(cè)模型進(jìn)行改進(jìn)，預(yù)測(cè)流程如圖2所示.

圖2 改進(jìn)混沌短期預(yù)測(cè)模型的預(yù)測(cè)流程Fig.2 Prediction process of improved chaotic short-term prediction model

圖2中，噪聲主導(dǎo)的IMF通過(guò)其自相關(guān)特征進(jìn)行辨識(shí)[16].預(yù)測(cè)方式I 和方式II 的區(qū)別在于合并IMF的先后次序不同：方式I將合并后的風(fēng)險(xiǎn)序列輸入至多步預(yù)測(cè)器中，計(jì)算預(yù)測(cè)值；方式II 將降噪后和其余各IMF 分別輸入至多步預(yù)測(cè)器中，計(jì)算出各自的預(yù)測(cè)值后，合并得到.

4 實(shí)例分析

4.1 混沌識(shí)別

將1.1 節(jié)中1 096 d 風(fēng)險(xiǎn)序列數(shù)據(jù)作為混沌識(shí)別樣本，采用C-C 方法計(jì)算相空間重構(gòu)的最佳τ=4 、m=5，重構(gòu)相空間為，其中，ti=17,18,19,…,1 096.采用Wolf 方法計(jì)算最大Lyapunov 指數(shù)λ=0.328 9>0，說(shuō)明具有混沌特征，風(fēng)險(xiǎn)序列短期可預(yù)測(cè).混沌最大可預(yù)測(cè)時(shí)間長(zhǎng)度為t0=1λ=3，即在一般情況下，當(dāng)預(yù)測(cè)時(shí)間超過(guò)3 d 時(shí)，誤差將會(huì)放大.

4.2 經(jīng)驗(yàn)?zāi)B(tài)分解

運(yùn)用EMD 和EEMD 對(duì)風(fēng)險(xiǎn)序列進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，得到8個(gè)IMF分量和1個(gè)殘余項(xiàng)r，如圖3所示.可見(jiàn)：EMD 分解結(jié)果中，任一分量中均具有不同頻率的成分，出現(xiàn)明顯模態(tài)混疊現(xiàn)象；EEMD 分解結(jié)果中，從，頻率從高到低遞減，各IMF中僅含有頻率相近的成分，說(shuō)明模態(tài)混疊得到抑制.

圖3 EMD 和EEMD 分解后部分結(jié)果Fig.3 Partial EMD and EEMD datas

4.3 閾值降噪處理

和頻率高，噪聲特征明顯，故對(duì)兩者進(jìn)行閾值降噪.根據(jù)EEMD閾值降噪方法，頻率較高，噪聲成分較多，為使降噪后的IMF 平滑，閾值函數(shù)的折中因子α取值稍高；噪聲特征比弱，有用成分較多，防止降噪后的失真，將α取值稍低；經(jīng)過(guò)反復(fù)實(shí)驗(yàn)發(fā)現(xiàn)，，時(shí)，降噪后信噪比為7.275 9，效果較好.

為分析降噪后風(fēng)險(xiǎn)序列混沌特征，計(jì)算降噪后序列和IMF分量的τ、m和λ，如表2所示.經(jīng)降噪處理，λ從0.328 9下降至0.237 7，仍具備混沌特征且最大可預(yù)測(cè)范圍增為4 d.從計(jì)算結(jié)果可知，隨著階數(shù)增加，λ減小并接近于0，混沌特征變?nèi)?

4.4 預(yù)測(cè)計(jì)算與結(jié)果

將風(fēng)險(xiǎn)序列中前900 d作為訓(xùn)練集，剩余196 d作為測(cè)試集，測(cè)試過(guò)程根據(jù)3.4 節(jié)的預(yù)測(cè)方式I 和方式II預(yù)測(cè)風(fēng)險(xiǎn)序列未來(lái)1～7 d的風(fēng)險(xiǎn)值.以未來(lái)1，3，5 d的預(yù)測(cè)結(jié)果為例，如圖4所示.

表2 IMF 分量和降噪前后的混沌特性Table 2 Chaos characteristics of IMF and before and after de-noising

圖4 航班運(yùn)行風(fēng)險(xiǎn)的未來(lái)1，3，5 d 預(yù)測(cè)結(jié)果Fig.4 Prediction results of flight operations risks in next 1,3,and 5 days

為定量說(shuō)明預(yù)測(cè)效果，采用均方根誤差(Root Mean Square Error,RMSE)和修正的相對(duì)誤差均值(Mean Absolute Percentage Error, MAPE)評(píng)估預(yù)測(cè)精度.MAPE反映預(yù)測(cè)可信程度，計(jì)算后發(fā)現(xiàn)，結(jié)果相對(duì)誤差值的頻數(shù)分布呈高度右偏，為使MAPE更具代表性，去掉相對(duì)誤差最大5%和最小5%的數(shù)據(jù)后再計(jì)算MAPE，定義為修正MAPE 值.計(jì)算降噪前后預(yù)測(cè)結(jié)果的RMSE 和修正MAPE 值如表3所示.

表3 不同方式預(yù)測(cè)結(jié)果的RMSE 及修正MAPE 值Table 3 RMSE and revised MAPE for prediction results in different ways

由表3 可知：未來(lái)1～7 d 預(yù)測(cè)結(jié)果的RMSE 和修正MAPE 值逐步增大，說(shuō)明預(yù)測(cè)精度隨預(yù)測(cè)周期延長(zhǎng)逐步降低；經(jīng)EEMD閾值降噪處理，未來(lái)1～7 d 預(yù)測(cè)結(jié)果的MAPE 值明顯變小，預(yù)測(cè)方式II 的MAPE 值平均降低9%，表明降噪處理有助于提高預(yù)測(cè)精度.

根據(jù)該航空公司使用要求，當(dāng)預(yù)測(cè)值相對(duì)誤差低于20%時(shí)，結(jié)果具有實(shí)踐操作價(jià)值；介于20%～25%時(shí)，具有實(shí)踐參考價(jià)值.降噪后，方式II所得未來(lái)第1 天預(yù)測(cè)結(jié)果的修正MAPE 值僅為18.63%，滿足該公司操作標(biāo)準(zhǔn)，可依照結(jié)果開展次日風(fēng)險(xiǎn)防控工作.

降噪后，方式II的預(yù)測(cè)結(jié)果精度變化從第4天開始放緩，至第7天修正MAPE值仍可保持在25%以下，具有實(shí)踐參考價(jià)值.受限于航權(quán)、商務(wù)、空域等因素，國(guó)內(nèi)航空公司制定航班計(jì)劃以周為單位.該方案預(yù)測(cè)周期恰好對(duì)應(yīng)，可于一周前通過(guò)飛機(jī)調(diào)配和人員調(diào)度等方式調(diào)整計(jì)劃編排，提前降低日運(yùn)行風(fēng)險(xiǎn).

5 結(jié) 論

本文針對(duì)航班運(yùn)行風(fēng)險(xiǎn)預(yù)測(cè)問(wèn)題，在識(shí)別時(shí)間序列混沌特征的基礎(chǔ)上，構(gòu)建基于EEMD-ELM的航班運(yùn)行風(fēng)險(xiǎn)混沌短期預(yù)測(cè)模型，并以某航空公司1 096 d數(shù)據(jù)進(jìn)行驗(yàn)算.結(jié)果表明：航班運(yùn)行風(fēng)險(xiǎn)時(shí)間序列λ=0.328 9>0，具有混沌特征，可預(yù)測(cè)范圍為1～3 d；EEMD 方法可抑制序列的IMF 模態(tài)混疊現(xiàn)象，且在EEMD閾值降噪后，序列仍具備混沌特征且可預(yù)測(cè)范圍增長(zhǎng)至4 d；與降噪前相比，基于EEMD-ELM的7 d預(yù)測(cè)結(jié)果修正MAPE值平均降低了9%，說(shuō)明閾值降噪方法可有效降低誤差；使用EEMD-ELM 預(yù)測(cè)模型計(jì)算未來(lái)1 d 預(yù)測(cè)結(jié)果的MAPE 值僅為18.63%，對(duì)次日航班風(fēng)險(xiǎn)防控具有實(shí)踐操作價(jià)值，未來(lái)7 d預(yù)測(cè)結(jié)果對(duì)周航班任務(wù)的優(yōu)化編排具有實(shí)踐參考價(jià)值.在此基礎(chǔ)上，后續(xù)將對(duì)航班運(yùn)行風(fēng)險(xiǎn)不同時(shí)間維度的預(yù)測(cè)及精度提升等開展研究.