基于EEMD-ELM的航班運行風(fēng)險混沌預(yù)測

王巖韜，李景良，谷潤平

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

0 引 言

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

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

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

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

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

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

表1 以天為單位，統(tǒng)計國內(nèi)某大型航空公司2016—2018年連續(xù)1 096 d 的航班狀況，提取風(fēng)險管控系統(tǒng)和航空安全系統(tǒng)中所有差錯指標(biāo)，形成航班運行風(fēng)險數(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ù)可信.以航班運行總風(fēng)險時間序列(簡稱風(fēng)險序列)作為后續(xù)計算的標(biāo)準(zhǔn)樣本，如圖1所示.

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

1.2 混沌識別

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

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

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

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

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

2.1 短期預(yù)測方法

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

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

2.2 ELM方法

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

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

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

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

3 EEMD改進的風(fēng)險短期預(yù)測模型

3.1 EEMD方法

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

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

3.2 端點延拓

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

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

3.3 閾值降噪

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

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

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

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

3.4 改進預(yù)測模型構(gòu)建

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

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

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

4 實例分析

4.1 混沌識別

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

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

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

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

4.3 閾值降噪處理

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

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

4.4 預(yù)測計算與結(jié)果

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

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

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

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

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

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

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

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

5 結(jié) 論

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