基于EEMD-ELM的航班运行风险混沌预测

王岩韬，李景良，谷润平

(中国民航大学 空中交通管理学院，天津300300)

0 引 言

民航工作高质量发展依赖于航班运行风险防控能力的提升.欧美最早从1991年开始航班风险评价研究[1].研究理论和手段侧重于事故树等传统方法[2-3].实践方面，2007—2012年美国、加拿大民航局研发了Flight Risk Assessment Tool(FRAT)[4]、Recheck[5]等飞行风险评估工具，至今鲜见可信的使用报告和精度分析.国内，航班风险评估研究最早于1999年[6]，王岩韬等建立了一套航班运行风险评估体系[7]，运用多算法协作将风险分类正确率提高至95%[8].对于民航非动态风险评估，上述研究虽已取得较好成果，但将其应用于民航一线后发现，评估技术只能做到随航班要素变化而变化，只是对已有状态评估，而风险管控需要时间提前量，对此风险预测技术成为风险防控实质所亟需.

国外最早于2008年研究航空相关风险预测[9].2012年，欧盟启动航空运行主动安全管理项目PROSPERO(Proactive Safety Performance for Operations)预测航空系统风险，但未见详细方案和结果报告，最新消息表明该项目已于2015年戛然而止[10].Andrej Lalis 等[11]使用线性回归模型和ARMA模型对航空安全绩效进行预测，但以月为单位的宽幅时间序列严重削弱了其实用价值.Zhang等[12]融合支持向量机和深度神经网络算法预测航空事件严重程度，但其实质过程与文献[7]和[8]类似，对未来预测能力不足.王岩韬团队[13]于2019年提出基于动态贝叶斯网络的航班运行风险预测模型，主要针对单个航班的进近、着陆阶段进行风险预测，但无法体现航空公司整体运行状况.航班风险预测的专项研究很少，无论是小时、日、月等不同时间维度的预测方法，还是预测精度稳定性分析，亦或是预测精度提升技术，均缺少深入的研究，是民航安全领域尚未解决的关键问题之一.

本文旨在探索一种实用有效的航班运行风险短期预测方法.借鉴混沌时间序列预测理论，对数据进行相空间重构，实现混沌识别；基于集成经验模态分解(Ensemble Empirical Mode Decomposition,EEMD)进行阈值降噪处理，采用极限学习机(Extreme Learning Machine，ELM)构建改进的混沌短期预测模型，并对预测效果进行数值分析.

1 风险数据时间序列及其混沌识别

1.1 风险数据时间序列

航班运行总风险是指当日航班运行相关危险源导致的后果严重程度和发生概率的综合.中国民航制定咨询通告AC-121-FS-2015-125，使用数值1～10 量化评估了不同严重程度与发生概率组合.

表1 以天为单位，统计国内某大型航空公司2016—2018年连续1 096 d 的航班状况，提取风险管控系统和航空安全系统中所有差错指标，形成航班运行风险数据.使用SPSS 软件，取置信区间95%，使用SNK法和Dunnet-t法对数据列做多个均数多重比较，检验结果概率均小于0.05，说明数据间具有独立性；再经飞行、运控、机务、安监等联合专家组检验，确认数据可信.以航班运行总风险时间序列(简称风险序列)作为后续计算的标准样本，如图1所示.

表1 航班运行风险统计数据Table 1 Risk assessment statistical data

1.2 混沌识别

运用坐标延迟重构法进行相空间重构，设某一时刻ti的风险值为x(ti)，ti=1+τw,2+τw,3+τw,…,N，i为时间索引，样本量为N的风险序列为X(N)=[x(1),x(2),…,x(N)]，则ti时刻重构后的相空间为

式中：τw为延迟时间窗口，τw=(m-1) ·τ；τ为时间延迟，m为嵌入维，τ和m取值是相空间重构的关键，最佳τ和m采用C-C方法计算得出.

图1 风险序列数据样本Fig.1 Time series example of flight operation risk

重构相空间后，采用Wolf 方法计算最大Lyapunov 指数，识别风险序列的混沌特征.Lyapunov指数用于描述相邻相空间轨线随时间推移而偏离扩散程度，当最大Lyapunov 指数λ>0，意味着轨线出现指数扩散，认为风险序列具有混沌特征.

2 混沌短期预测模型构建

2.1 短期预测方法

短期预测方法分为直接预测和迭代预测.直接预测通过映射f(·)得到ti+p时刻的预测值，即，其中，p为预测周期.迭代预测：先确定映射g(·)，进行单步预测，得到ti+1 时刻预测值x(ti+1)=g[Y(ti)]；再根据x(ti+1) 重构相空间Y(ti+1) ，重新执行单步预测，得到预测值x(ti+2)；反复迭代p次，得到ti+p时刻预测值.直接预测不存在误差累积效应，但p较大时，难以构建满足精度要求的预测模型；迭代预测过程依赖于单步预测，误差随迭代次数增加而增大，当输入全为预测值时，误差会急剧增大.

多次实验发现，采用迭代预测方法优于直接预测，故后文结果均通过迭代预测得到.

2.2 ELM方法

ELM 是一种单隐层前馈神经网络，具有学习速度快，泛化性能好等特点，广泛应用于混沌时间序列预测[14].对于Q个样本其中，j为样本索引，Ij为输入层，，k为输入层神经元个数；Oj为期望输出值，，l为输出层神经元个数.设隐含层神经元个数为n，隐含层激活函数为g(·)，隐含层某一神经元与输入层、输出层的连接权值分别为和，其中为第i(i=1,2,…,n)个隐含层神经元与输入层第e(e=1,2,…,k)个神经元的连接权值为第i个隐含层神经元与输出层第u(u=1,2,…,l)个神经元的连接权值.采用ELM逼近具有Q个样本的目标输出时，可得

式中：H为隐含层输出矩阵 ，bi为神经网络的阈值，i=1,2,…,n；w为隐含层与输出层的连接权值，O为期望输出值，.

在ELM 训练过程，隐含层与输入层的连接权值win和阈值bi(i=1,2,…,n)随机产生，并在训练过程保持不变；只需设置隐含层神经元个数和激活函数，w求解公式为

其解为，H†为H的伪逆.

3 EEMD改进的风险短期预测模型

3.1 EEMD方法

经验模态分解 (Empirical Mode Decomposition, EMD)对序列逐级分解，产生多个具有不同特征尺度的本征模态函数IMF，即

使用集成经验模态分解(Ensemble Empirical Mode Decomposition,EEMD)是为克服EMD 的模态混叠现象.EEMD 改进之处为：原序列加入白噪声序列后，进行EMD 分解；反复加入不同的白噪声序列，反复分解；当重复次数足够多时，对IMF和残余项分别求平均，抵消噪声，得到EEMD分解结果.

3.2 端点延拓

端点效应是指分解过程中因序列两端点不一定是极值点，导致上、下包络线在序列两端出现发散的现象，该现象会逐渐对内部数据产生影响，使分解结果失真.对ti时刻的风险序列预测，模型输入向量属于序列X(ti)的端点数据，经过计算发现，端点效应对预测结果产生了明显的负面影响.

采用基于预测值的端点延拓方法，避免端点效应.根据第2 节混沌短期预测方法，计算未来一定时间范围的预测值，保证至少存在极大值点和极小值点；将原序列与预测序列连接，在分解前实现原序列的端点延拓.

3.3 阈值降噪

借鉴小波阈值降噪思路，选择合适阈值和阈值函数截断IMF，把处理后的IMF 合成新序列，生成降噪结果.

式中：为中噪声的方差，.

阈值函数分为软阈值函数和硬阈值函数.软阈值函数会使降噪结果失真，硬阈值函数会使降噪结果不连续.采用改进的软硬折中阈值函数[16]，即

式中：α为折中因子；为的降噪结果；sign(·) 为符号函数.

3.4 改进预测模型构建

根据EEMD阈值降噪方法对风险短期预测模型进行改进，预测流程如图2所示.

图2 改进混沌短期预测模型的预测流程Fig.2 Prediction process of improved chaotic short-term prediction model

图2中，噪声主导的IMF通过其自相关特征进行辨识[16].预测方式I 和方式II 的区别在于合并IMF的先后次序不同：方式I将合并后的风险序列输入至多步预测器中，计算预测值；方式II 将降噪后和其余各IMF 分别输入至多步预测器中，计算出各自的预测值后，合并得到.

4 实例分析

4.1 混沌识别

将1.1 节中1 096 d 风险序列数据作为混沌识别样本，采用C-C 方法计算相空间重构的最佳τ=4 、m=5，重构相空间为，其中，ti=17,18,19,…,1 096.采用Wolf 方法计算最大Lyapunov 指数λ=0.328 9>0，说明具有混沌特征，风险序列短期可预测.混沌最大可预测时间长度为t0=1λ=3，即在一般情况下，当预测时间超过3 d 时，误差将会放大.

4.2 经验模态分解

运用EMD 和EEMD 对风险序列进行经验模态分解，得到8个IMF分量和1个残余项r，如图3所示.可见：EMD 分解结果中，任一分量中均具有不同频率的成分，出现明显模态混叠现象；EEMD 分解结果中，从，频率从高到低递减，各IMF中仅含有频率相近的成分，说明模态混叠得到抑制.

图3 EMD 和EEMD 分解后部分结果Fig.3 Partial EMD and EEMD datas

4.3 阈值降噪处理

和频率高，噪声特征明显，故对两者进行阈值降噪.根据EEMD阈值降噪方法，频率较高，噪声成分较多，为使降噪后的IMF 平滑，阈值函数的折中因子α取值稍高；噪声特征比弱，有用成分较多，防止降噪后的失真，将α取值稍低；经过反复实验发现，，时，降噪后信噪比为7.275 9，效果较好.

为分析降噪后风险序列混沌特征，计算降噪后序列和IMF分量的τ、m和λ，如表2所示.经降噪处理，λ从0.328 9下降至0.237 7，仍具备混沌特征且最大可预测范围增为4 d.从计算结果可知，随着阶数增加，λ减小并接近于0，混沌特征变弱.

4.4 预测计算与结果

将风险序列中前900 d作为训练集，剩余196 d作为测试集，测试过程根据3.4 节的预测方式I 和方式II预测风险序列未来1～7 d的风险值.以未来1，3，5 d的预测结果为例，如图4所示.

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

图4 航班运行风险的未来1，3，5 d 预测结果Fig.4 Prediction results of flight operations risks in next 1,3,and 5 days

为定量说明预测效果，采用均方根误差(Root Mean Square Error,RMSE)和修正的相对误差均值(Mean Absolute Percentage Error, MAPE)评估预测精度.MAPE反映预测可信程度，计算后发现，结果相对误差值的频数分布呈高度右偏，为使MAPE更具代表性，去掉相对误差最大5%和最小5%的数据后再计算MAPE，定义为修正MAPE 值.计算降噪前后预测结果的RMSE 和修正MAPE 值如表3所示.

表3 不同方式预测结果的RMSE 及修正MAPE 值Table 3 RMSE and revised MAPE for prediction results in different ways

由表3 可知：未来1～7 d 预测结果的RMSE 和修正MAPE 值逐步增大，说明预测精度随预测周期延长逐步降低；经EEMD阈值降噪处理，未来1～7 d 预测结果的MAPE 值明显变小，预测方式II 的MAPE 值平均降低9%，表明降噪处理有助于提高预测精度.

根据该航空公司使用要求，当预测值相对误差低于20%时，结果具有实践操作价值；介于20%～25%时，具有实践参考价值.降噪后，方式II所得未来第1 天预测结果的修正MAPE 值仅为18.63%，满足该公司操作标准，可依照结果开展次日风险防控工作.

降噪后，方式II的预测结果精度变化从第4天开始放缓，至第7天修正MAPE值仍可保持在25%以下，具有实践参考价值.受限于航权、商务、空域等因素，国内航空公司制定航班计划以周为单位.该方案预测周期恰好对应，可于一周前通过飞机调配和人员调度等方式调整计划编排，提前降低日运行风险.

5 结 论

本文针对航班运行风险预测问题，在识别时间序列混沌特征的基础上，构建基于EEMD-ELM的航班运行风险混沌短期预测模型，并以某航空公司1 096 d数据进行验算.结果表明：航班运行风险时间序列λ=0.328 9>0，具有混沌特征，可预测范围为1～3 d；EEMD 方法可抑制序列的IMF 模态混叠现象，且在EEMD阈值降噪后，序列仍具备混沌特征且可预测范围增长至4 d；与降噪前相比，基于EEMD-ELM的7 d预测结果修正MAPE值平均降低了9%，说明阈值降噪方法可有效降低误差；使用EEMD-ELM 预测模型计算未来1 d 预测结果的MAPE 值仅为18.63%，对次日航班风险防控具有实践操作价值，未来7 d预测结果对周航班任务的优化编排具有实践参考价值.在此基础上，后续将对航班运行风险不同时间维度的预测及精度提升等开展研究.