基于多变量混沌时间序列的航班运行风险预测模型

工程科学学报 Chinese Journal of Engineering 基于多变量混沌时间序列的航班运行风险预测棋型 王岩韬李景良谷润平 Flight operation risk prediction model based on the multivariate chaotic time series WANG Yan-tao,LI Jing-liang.GU Run-ping 引用本文： 王岩韬，李景良，谷润平.基于多变量混沌时间序列的航班运行风险预测模型D.工程科学学报，2020,42(12：1664-1673. doi10.13374j.issn2095-9389.2019.12.09.002 WANG Yan-tao,LI Jing-liang.GU Run-ping.Flight operation risk prediction model based on the multivariate chaotic time series[J].Chinese Journal of Engineering,2020,42(12):1664-1673.doi:10.13374/j.issn2095-9389.2019.12.09.002 在线阅读View online::htps:/ldoi.org10.13374.issn2095-9389.2019.12.09.002 您可能感兴趣的其他文章 Articles you may be interested in BP神经网络F钢铝耗的预测模型 Prediction model of aluminum consumption with BP neural networks in IF steel production 工程科学学报.2017,394：511 https::/1doi.org10.13374.issn2095-9389.2017.04.005 基于多维时间序列形态特征的相似性动态聚类算法 Similarity dynamical clustering algorithm based on multidimensional shape features for time series 工程科学学报.2017.397)：1114htps:/doi.org10.13374.issn2095-9389.2017.07.019 基于机器学习的北京市PM2.5浓度预测模型及模拟分析 Machine-learning-based model and simulation analysis of PM2.5 concentration prediction in Beijing 工程科学学报.2019,41(3)401 https::/1doi.org/10.13374.issn2095-9389.2019.03.014 深度神经网络模型压缩综述 A survey of model compression for deep neural networks 工程科学学报.2019,41(10：1229htps:oi.org10.13374.issn2095-9389.2019.03.27.002 20 CrMnTit齿轮钢的点蚀敏感性及裂纹萌生风险 Pitting sensitivity and crack initiation risk of 20CrMnTi gear steel 工程科学学报.2017,395：731 https:/doi.org10.13374/j.issn2095-9389.2017.05.011 基于极限学习机ELM的连铸坯质量预测 Quality prediction of the continuous casting bloom based on the extreme learning machine 工程科学学报.2018.40(7)：815 https:1doi.org/10.13374issn2095-9389.2018.07.007

基于多变量混沌时间序列的航班运行风险预测模型 王岩韬 李景良 谷润平 Flight operation risk prediction model based on the multivariate chaotic time series WANG Yan-tao, LI Jing-liang, GU Run-ping 引用本文: 王岩韬, 李景良, 谷润平. 基于多变量混沌时间序列的航班运行风险预测模型[J]. 工程科学学报, 2020, 42(12): 1664-1673. doi: 10.13374/j.issn2095-9389.2019.12.09.002 WANG Yan-tao, LI Jing-liang, GU Run-ping. Flight operation risk prediction model based on the multivariate chaotic time series[J]. Chinese Journal of Engineering, 2020, 42(12): 1664-1673. doi: 10.13374/j.issn2095-9389.2019.12.09.002 在线阅读 View online: https://doi.org/10.13374/j.issn2095-9389.2019.12.09.002 您可能感兴趣的其他文章 Articles you may be interested in BP神经网络IF钢铝耗的预测模型 Prediction model of aluminum consumption with BP neural networks in IF steel production 工程科学学报. 2017, 39(4): 511 https://doi.org/10.13374/j.issn2095-9389.2017.04.005 基于多维时间序列形态特征的相似性动态聚类算法 Similarity dynamical clustering algorithm based on multidimensional shape features for time series 工程科学学报. 2017, 39(7): 1114 https://doi.org/10.13374/j.issn2095-9389.2017.07.019 基于机器学习的北京市PM2.5浓度预测模型及模拟分析 Machine-learning-based model and simulation analysis of PM2.5 concentration prediction in Beijing 工程科学学报. 2019, 41(3): 401 https://doi.org/10.13374/j.issn2095-9389.2019.03.014 深度神经网络模型压缩综述 A survey of model compression for deep neural networks 工程科学学报. 2019, 41(10): 1229 https://doi.org/10.13374/j.issn2095-9389.2019.03.27.002 20CrMnTi齿轮钢的点蚀敏感性及裂纹萌生风险 Pitting sensitivity and crack initiation risk of 20CrMnTi gear steel 工程科学学报. 2017, 39(5): 731 https://doi.org/10.13374/j.issn2095-9389.2017.05.011 基于极限学习机(ELM)的连铸坯质量预测 Quality prediction of the continuous casting bloom based on the extreme learning machine 工程科学学报. 2018, 40(7): 815 https://doi.org/10.13374/j.issn2095-9389.2018.07.007

工程科学学报.第42卷.第12期：1664-1673.2020年12月 Chinese Journal of Engineering,Vol.42,No.12:1664-1673,December 2020 https://doi.org/10.13374/j.issn2095-9389.2019.12.09.002;http://cje.ustb.edu.cn 基于多变量混沌时间序列的航班运行风险预测模型 王岩韬四，李景良，谷润平 中国民航大学空管学院，天津300300 通信作者，E-mail:CAUCwyt@126.com 摘要为了提升航班运行风险预测精度，基于某航空公司2016一2018年航班运行风险数据，在验证15个风险时间序列的 混沌特性后，构建基于多变量混沌时间序列的风险预测模型.首先，对15个风险时间序列进行多变量相空间重构，采用主成 分分析法(PCA)对相空间进行降维处理：然后，基于迭代预测的方式，分别采用极限学习机、RBF神经网络、回声状态网络和 Elm神经网络建立风险短期预测模型：最后，以降维后的相空间作为输入.计算并比较分析未来1~7d的风险预测结果.结 果表明：多变量相空间重构后总维数为62维.经PCA降维处理.降至31维：在不同的预测模型中，降维后RBF模型预测效果 最佳：其中，预测第1天结果相对误差<25%出现频数为82.62%.至第5天仍达75%以上：该模型第1天预测结果的修正平均 绝对百分比误差(MAPE)值为11.32%，且前5d均低于20%.满足航空公司使用要求，1~5d预测结果对航班风险管控具有实 践操作价值，证明基于多变量混沌时间序列的风险预测方案可行、有效. 关键词航班运行风险：风险预测：多变量混沌时间序列：相空间重构：神经网络 分类号N945.24;X949U8 Flight operation risk prediction model based on the multivariate chaotic time series WANG Yan-tao,LI Jing-liang.GU Run-ping School of Air Traffic Management,Civil Aviation University of China,Tianjin 300300,China Corresponding author,E-mail:CAUCwyt@126.com ABSTRACT With the development of civil aviation safety management,the flight operation risk of airlines is of increasing concern. Risk prediction technology extracts information from historical and current risk data and uses it to predict short-term trends in the future. thus helping identify emerging risks and providing more time for risk management.Compared with non-dynamic risk assessment,this technology is more substantial for the management and control of flight operation risk.To improve the accuracy of flight operation risk prediction,on the basis of the flight risk data of a certain airline in 2016-2018,the chaotic characteristics of 15 risk time series were verified and a short-term risk prediction model based on the multivariate chaotic time series was constructed.First,multivariate phase space reconstruction was performed on 15 risk time series,and the phase space was reduced by the principal component analysis(PCA) method.Then,four short-term risk prediction models,namely,extreme learning machine,radial basis function(RBF)neural network, echo state network,and Elman neural network,were built on the basis of iterative prediction.Finally,the phase space after dimension reduction was used as the model input,and the risk prediction results for 1-7 d were calculated and compared.Results show that the total number of dimensions after multivariable phase space reconstruction is 62,which is reduced to 31 by PCA dimension reduction.Of the four prediction models,the RBF neural network model after dimension reduction has the best prediction effect.The occurrence frequency of<25%relative error is 82.62%for the first day and 75%for the fifth day.The corrected mean absolute percentage error for the first day is 11.32%,and lower than 20%for the next 4 d.Thus,the calculation results meet the requirements of the airline.The prediction results within 1-5 d have practical value for flight risk management,proving that the risk prediction method based on the 收稿日期：2019-12-09 基金项目：国家自然科学基金资助项目(01933103)

基于多变量混沌时间序列的航班运行风险预测模型 王岩韬苣，李景良，谷润平 中国民航大学空管学院，天津 300300 苣通信作者，E-mail：CAUCwyt@126.com 摘    要    为了提升航班运行风险预测精度，基于某航空公司 2016—2018 年航班运行风险数据，在验证 15 个风险时间序列的 混沌特性后，构建基于多变量混沌时间序列的风险预测模型. 首先，对 15 个风险时间序列进行多变量相空间重构，采用主成 分分析法（PCA）对相空间进行降维处理；然后，基于迭代预测的方式，分别采用极限学习机、RBF 神经网络、回声状态网络和 Elman 神经网络建立风险短期预测模型；最后，以降维后的相空间作为输入，计算并比较分析未来 1～7 d 的风险预测结果. 结 果表明：多变量相空间重构后总维数为 62 维，经 PCA 降维处理，降至 31 维；在不同的预测模型中，降维后 RBF 模型预测效果 最佳；其中，预测第 1 天结果相对误差<25% 出现频数为 82.62%，至第 5 天仍达 75% 以上；该模型第 1 天预测结果的修正平均 绝对百分比误差（MAPE）值为 11.32%，且前 5 d 均低于 20%，满足航空公司使用要求. 1～5 d 预测结果对航班风险管控具有实 践操作价值，证明基于多变量混沌时间序列的风险预测方案可行、有效. 关键词    航班运行风险；风险预测；多变量混沌时间序列；相空间重构；神经网络 分类号    N945.24; X949; U8 Flight operation risk prediction model based on the multivariate chaotic time series WANG Yan-tao苣 ，LI Jing-liang，GU Run-ping School of Air Traffic Management, Civil Aviation University of China, Tianjin 300300, China 苣 Corresponding author, E-mail: CAUCwyt@126.com ABSTRACT    With the development of civil aviation safety management, the flight operation risk of airlines is of increasing concern. Risk prediction technology extracts information from historical and current risk data and uses it to predict short-term trends in the future, thus helping identify emerging risks and providing more time for risk management. Compared with non-dynamic risk assessment, this technology is more substantial for the management and control of flight operation risk. To improve the accuracy of flight operation risk prediction, on the basis of the flight risk data of a certain airline in 2016—2018, the chaotic characteristics of 15 risk time series were verified and a short-term risk prediction model based on the multivariate chaotic time series was constructed. First, multivariate phase space reconstruction was performed on 15 risk time series, and the phase space was reduced by the principal component analysis (PCA) method. Then, four short-term risk prediction models, namely, extreme learning machine, radial basis function (RBF) neural network, echo state network, and Elman neural network, were built on the basis of iterative prediction. Finally, the phase space after dimension reduction was used as the model input, and the risk prediction results for 1–7 d were calculated and compared. Results show that the total number of dimensions after multivariable phase space reconstruction is 62, which is reduced to 31 by PCA dimension reduction. Of the four prediction models, the RBF neural network model after dimension reduction has the best prediction effect. The occurrence frequency of <25% relative error is 82.62% for the first day and 75% for the fifth day. The corrected mean absolute percentage error for the first day is 11.32%, and lower than 20% for the next 4 d. Thus, the calculation results meet the requirements of the airline. The prediction results within 1–5 d have practical value for flight risk management, proving that the risk prediction method based on the 收稿日期: 2019−12−09 基金项目: 国家自然科学基金资助项目（U1933103） 工程科学学报，第 42 卷，第 12 期：1664−1673，2020 年 12 月 Chinese Journal of Engineering, Vol. 42, No. 12: 1664−1673, December 2020 https://doi.org/10.13374/j.issn2095-9389.2019.12.09.002; http://cje.ustb.edu.cn

王岩韬等：基于多变量混沌时间序列的航班运行风险预测模型 1665 multivariate chaotic time series is feasible and effective. KEY WORDS flight operation risk;risk prediction;multivariate chaotic time series;phase space reconstruction;artificial neural networks 近年不安全事件表明，民航工作高质量发展 术，早期应用于经典的混沌序列中，证明了其比单 依赖于航班运行风险防控能力的提升川针对航班 变量预测更优4-近年来，该方法广泛应用于工 运行风险影响因素繁杂、随时间变化、等级界定 程领域，典型研究成果包括：2013年，侯公羽等7 困难等问题，欧美国家最早从1991年开始航班风 对煤矿斜井的盾构小时施工风险进行预测，均方 险评价研究四从近年研究成果看，研究理论和手 根误差(Root mean square error,.RMSE)为0.4754； 段仍使用事故树等传统方法而在实践方面， 2016年，钟仪华等1对油田日产油量进行预测， 自2007一2012年间美国、加拿大民航局研发了 均方根误差为0.0061：2017年，杜柳青等对数控 Flight risk assessment tool(FRAT)、RA check!等 机床的日运动精度进行预测，均方根误差低于 飞行风险评估工具后，至今未有可信的使用报告 0.02:2018年，张淑清等20对电力日负荷进行预 和精度分析；除此之外，近年也未见其他突破.在 测，平均误差为0.813%：2018年，黄发明等2对滑 国内，航班风险评估研究最早从1999年起步)，而 坡的月波动项位移进行预测，均方根误差为23.71 后王岩韬等⑧]建立了具有指导性、实用性的航班 上述工程应用均表明多变量混沌时间序列预测方 运行风险评估体系，并运用多算法协作的方式将 法具有良好的预测效果 航班风险分类正确率提高至95%.对于民航非动 由上，本文首先对运行数据时间序列进行混 态风险评估，上述研究已取得较好成果.但将前续 沌识别和相空间重构：然后，采用基于奇异值分解 成果应用于民航一线运行后发现，航班评估技术 (Singular value decomposition,.SVD)的主成分分析 只能做到随航班要素快速变化而变化，评估结果 (Principal component analysis,PCA)法对相空间进 代表已有状态，而风险管控工作需要时间提前量， 行降维处理；进而，构建极限学习机(Extreme 因此风险预测技术才是风险防控实质所亟需 learning machine,.ELM)、RBF神经网络、回声状态 2012年，欧盟开展航空运行主动性安全管理 网络(Echo state network,ESN)和Elman神经网络 Proactive safety performance for operations, 4种风险预测模型，采取迭代方式进行风险预测： PROSPERO),旨在通过对航空运输系统的风险主 最后计算并分析不同模型的预测精度，判断其实 动预测，降低人为差错和空中交通事故率，但详细 方案和结果报告一直未见报道，并且最新消息表 践操作价值，旨在构建一种精度较高的航班运行 明该项目已于2015年戛然而止10.2018年Lalis等山 风险预测方法， 使用线性回归模型和ARMA模型对航空安全绩效 1航班运行风险时间序列及混沌特征识别 指标进行预测，但以月为单位的、本身就波动很小 的时间序列严重削弱了其实用价值.2019年， 1.1航班运行风险 Zhang与Mahadevans!2I融合支持向量机和深度神 “风险”一词最早仅代表客观的危险，狭义风 经网络两种算法应用于航空风险预测中，对单一 险是指损失的不确定性，如工程建筑风险关注人 事件预测精度达到了81%，但其过程实质上与风 员伤亡和设备损失；而广义风险强调表现的不确 险评估类似，预测周期过短，预测能力不足.在国 定性，如经营风险关注未来市场走向，如金融风险 内，王岩韬团队于2019年提出了基于动态贝叶斯 关注收益变化.“风险”最受认可的定义是由 网络的航班运行风险预测模型]，但该模型主要 Coleman与Marks于1999年提出的，风险是结果 针对单个航班的进近、着陆阶段的风险预测，在预 严重程度和发生概率的综合四国际民航组织 测周期很短的情况下精度达到了80.4%，无法体现 (ICAO)附件19中，延用了此定义，风险为某一危 航空公司的整体运行状况.由上，航班风险预测相 险源预计产生的可能性与严重性相乘，并以风险 关研究较少，至今缺乏精度较高的技术和方法，是 矩阵划分等级.如表1，绿区为可接受风险，即航 民航安全领域尚未有效解决的关键问题之一. 班可以正常执行：红区为不可接受风险，即航班不 深入分析其他领域的研究成果，发现多变量 可执行：而黄区需要进行有效风险缓解后，航班才 混沌时间序列预测方法是一种较好的短期预测技 可执行

multivariate chaotic time series is feasible and effective. KEY  WORDS    flight operation risk； risk prediction； multivariate chaotic time series； phase space reconstruction； artificial neural networks 近年不安全事件表明，民航工作高质量发展 依赖于航班运行风险防控能力的提升[1] . 针对航班 运行风险影响因素繁杂、随时间变化、等级界定 困难等问题，欧美国家最早从 1991 年开始航班风 险评价研究[2] . 从近年研究成果看，研究理论和手 段仍使用事故树等传统方法[3−4] . 而在实践方面， 自 2007—2012 年间美国、加拿大民航局研发了 Flight risk assessment tool (FRAT)[5]、 RA check[6] 等 飞行风险评估工具后，至今未有可信的使用报告 和精度分析；除此之外，近年也未见其他突破. 在 国内，航班风险评估研究最早从 1999 年起步[7] ，而 后王岩韬等[8−9] 建立了具有指导性、实用性的航班 运行风险评估体系，并运用多算法协作的方式将 航班风险分类正确率提高至 95%. 对于民航非动 态风险评估，上述研究已取得较好成果. 但将前续 成果应用于民航一线运行后发现，航班评估技术 只能做到随航班要素快速变化而变化，评估结果 代表已有状态，而风险管控工作需要时间提前量， 因此风险预测技术才是风险防控实质所亟需. 2012 年，欧盟开展航空运行主动性安全管理 项 目 （ Proactive safety performance for operations, PROSPERO），旨在通过对航空运输系统的风险主 动预测，降低人为差错和空中交通事故率，但详细 方案和结果报告一直未见报道，并且最新消息表 明该项目已于2015 年戛然而止[10] . 2018 年Lalis 等[11] 使用线性回归模型和 ARMA 模型对航空安全绩效 指标进行预测，但以月为单位的、本身就波动很小 的时间序列严重削弱了其实用价值 . 2019 年 ， Zhang 与 Mahadevans[12] 融合支持向量机和深度神 经网络两种算法应用于航空风险预测中，对单一 事件预测精度达到了 81%，但其过程实质上与风 险评估类似，预测周期过短，预测能力不足. 在国 内，王岩韬团队于 2019 年提出了基于动态贝叶斯 网络的航班运行风险预测模型[13] ，但该模型主要 针对单个航班的进近、着陆阶段的风险预测，在预 测周期很短的情况下精度达到了 80.4%，无法体现 航空公司的整体运行状况. 由上，航班风险预测相 关研究较少，至今缺乏精度较高的技术和方法，是 民航安全领域尚未有效解决的关键问题之一. 深入分析其他领域的研究成果，发现多变量 混沌时间序列预测方法是一种较好的短期预测技 术，早期应用于经典的混沌序列中，证明了其比单 变量预测更优[14−16] . 近年来，该方法广泛应用于工 程领域，典型研究成果包括：2013 年，侯公羽等[17] 对煤矿斜井的盾构小时施工风险进行预测，均方 根误差 （ Root mean square error, RMSE） 为 0.4754； 2016 年，钟仪华等[18] 对油田日产油量进行预测， 均方根误差为 0.0061；2017 年，杜柳青等[19] 对数控 机床的日运动精度进行预测 ，均方根误差低于 0.02； 2018 年，张淑清等[20] 对电力日负荷进行预 测，平均误差为 0.813%；2018 年，黄发明等[21] 对滑 坡的月波动项位移进行预测，均方根误差为 23.71. 上述工程应用均表明多变量混沌时间序列预测方 法具有良好的预测效果. 由上，本文首先对运行数据时间序列进行混 沌识别和相空间重构；然后，采用基于奇异值分解 （Singular value decomposition, SVD）的主成分分析 （Principal component analysis, PCA）法对相空间进 行降维处理 ；进而 ，构建极限学习机 （ Extreme learning machine, ELM）、RBF 神经网络、回声状态 网络（Echo state network, ESN）和 Elman 神经网络 4 种风险预测模型，采取迭代方式进行风险预测； 最后计算并分析不同模型的预测精度，判断其实 践操作价值，旨在构建一种精度较高的航班运行 风险预测方法， 1    航班运行风险时间序列及混沌特征识别 1.1    航班运行风险 “风险”一词最早仅代表客观的危险，狭义风 险是指损失的不确定性，如工程建筑风险关注人 员伤亡和设备损失；而广义风险强调表现的不确 定性，如经营风险关注未来市场走向，如金融风险 关注收益变化 . “ 风 险 ” 最受认可的定义是 由 Coleman 与 Marks 于 1999 年提出的，风险是结果 严重程度和发生概率的综合[22] . 国际民航组织 （ICAO）附件 19 中，延用了此定义，风险为某一危 险源预计产生的可能性与严重性相乘，并以风险 矩阵划分等级. 如表 1，绿区为可接受风险，即航 班可以正常执行；红区为不可接受风险，即航班不 可执行；而黄区需要进行有效风险缓解后，航班才 可执行. 王岩韬等： 基于多变量混沌时间序列的航班运行风险预测模型 · 1665 ·

1666 工程科学学报，第42卷，第12期 表1ICAO风险矩阵与中国民航局(CAAC)风险值 该体系具有一定科学性，应用于数值化评估与预 Table 1 ICAO risk matrix and CAAC risk value 测时存在两个问题：①尚有少量的半定性/半定量 Severity 和定性指标，需要依赖专家判断：②部分指标不具 Likelihood Catastrophic,Hazardous,Major,Minor,Negligible, 备可预测性，如人因和航空器类指标是不含时间 A B C D E 维度特性的 Frequent,5 5A/10 5B/0 5/10 5D/7 5E/6 1.2运行数据选取 Occasional,4 4A/10 48/8 4C/8 4D/6 4E/5 由上，延续AC-121-FS-2015-125中对风险评定 Remote,3 3A/9 3B/8 3C/7 3D15 3E/2 方式，使用可定量的数据来源.提取国内某大型航 Improbable,2 2A/8 2B/7 2C/6 2D/4 2E/1 空公司风险管控系统和航空安全系统中风险源指 Extremely improbable,1 1A/7 1B/4 1C/3 D/1 1E/1 标，统计2016一2018年连续以天为单位的航班数 据，形成航班运行风险数据项1096组，数据类型 但是在使用风险矩阵评判航班风险等级时， 和数据样本见表2和表3. 存在两个问题：①对于中国民航，任一人员伤亡、 首先，处理15个风险项的数据，使用SPSS软 飞机损伤报废都无法接受，风险严重程度无法如 件，取置信区间95%，使用SNK法和Dunnet-t法对 金融行业等以金钱数字衡量；②根据《2019年民航 数据列做多个均数多重比较，检验结果概率均小 机场生产统计公报》，我国持续安全记录已达 于0.05，说明数据间具有独立性：然后，再经飞行、 112个月、8068×10h,而2019年航班已达1166万 运控、机务、安监等联合专家组检验，确认数据可 架次，事故、事故症候、不安全事件等发生概率接 信，故作为后续计算的标准样本 近极小值 1.3混沌识别 因此，为了解决风险矩阵的半定量化，进而发 根据航班运行风险数据样本，分别构建A1~ 展风险评估与预测技术，民航局于2015年颁布《航 A5和总风险时间序列，其中总风险时间序列如 空承运人运行控制风险管控系统实施指南(AC 图1.采用Wolf法计算，结果表明总风险时间序列 121-FS-2015-125)》以政府指导的方式量化了风 具有混沌特征，同理对15个风险源时间序列进行 险，经实践多年已验证有效.其中，将风险矩阵中 混沌识别，计算结果见表4. 1E~5A转化为数字1到10，即风险数值化，如表1 从表4可见，16个风险序列均满足1>0，具有 表2航班运行风险统计数据 Table 2 Risk assessment statistical data Risk item Statistical meaning 4 Statistics of captain (non-instructor)and second co-pilot match times Statistics of crew duty time or flight time less than I h times 学 Statistics of quick access recorder(QAR)blue and yellow warnings times after flight 4 Statistics times of captain's experience less than 200 h and the total flight experience less than 3000 h 4 Statistics of the crew first match times to special airport in 12 calendar months 46 temporary airborne failures times Number of minimum equipment list(MEL)/configur cy list(CDL)reservations or MEL reservations affecting near landings As Statistics times cargo or dangerous goods As nd handling error times 410 Statistics of low fuel abnormal altitude,etc.in flight monitoring A Statistics of ti d to on the edge of the weather criteria An Statistics of navigation equipment degradations,operating standards and obstacles temporary changes Au Statistics of special cing,thunders rms,etc.,)which need to divert A14 Statistics of special operations(polar oper re-dispatch flights,extended-range operations(ETOPS). extended cross-water,overflying unmanned areas) Ais Statistics of route closure,flow control,height limit,etc. Total risk Comprehensive assessment from the above risk items

但是在使用风险矩阵评判航班风险等级时， 存在两个问题：①对于中国民航，任一人员伤亡、 飞机损伤报废都无法接受，风险严重程度无法如 金融行业等以金钱数字衡量；②根据《2019 年民航 机场生产统计公报 》 ，我国持续安全记录已 达 112 个月、8068×104 h，而 2019 年航班已达 1166 万 架次，事故、事故症候、不安全事件等发生概率接 近极小值. 因此，为了解决风险矩阵的半定量化，进而发 展风险评估与预测技术，民航局于 2015 年颁布《航 空承运人运行控制风险管控系统实施指南（AC- 121-FS-2015-125） 》以政府指导的方式量化了风 险，经实践多年已验证有效. 其中，将风险矩阵中 1E～5A 转化为数字 1 到 10，即风险数值化，如表 1. 该体系具有一定科学性，应用于数值化评估与预 测时存在两个问题：①尚有少量的半定性/半定量 和定性指标，需要依赖专家判断；②部分指标不具 备可预测性，如人因和航空器类指标是不含时间 维度特性的. 1.2    运行数据选取 由上，延续 AC-121-FS-2015-125 中对风险评定 方式，使用可定量的数据来源. 提取国内某大型航 空公司风险管控系统和航空安全系统中风险源指 标，统计 2016—2018 年连续以天为单位的航班数 据，形成航班运行风险数据项 1096 组，数据类型 和数据样本见表 2 和表 3. 首先，处理 15 个风险项的数据，使用 SPSS 软 件，取置信区间 95%，使用 SNK 法和 Dunnet-t 法对 数据列做多个均数多重比较，检验结果概率均小 于 0.05，说明数据间具有独立性；然后，再经飞行、 运控、机务、安监等联合专家组检验，确认数据可 信，故作为后续计算的标准样本. 1.3    混沌识别 根据航班运行风险数据样本，分别构建 A1～ A15 和总风险时间序列，其中总风险时间序列如 图 1. 采用 Wolf 法计算，结果表明总风险时间序列 具有混沌特征，同理对 15 个风险源时间序列进行 混沌识别，计算结果见表 4. 从表 4 可见，16 个风险序列均满足 λ > 0 ，具有 表 1 ICAO 风险矩阵与中国民航局（CAAC）风险值 Table 1 ICAO risk matrix and CAAC risk value Likelihood Severity Catastrophic, A Hazardous, B Major, C Minor, D Negligible, E Frequent, 5 5A / 10 5B / 10 5C / 10 5D / 7 5E / 6 Occasional, 4 4A / 10 4B / 8 4C / 8 4D / 6 4E / 5 Remote, 3 3A / 9 3B / 8 3C / 7 3D / 5 3E / 2 Improbable, 2 2A / 8 2B / 7 2C / 6 2D / 4 2E / 1 Extremely improbable, 1 1A / 7 1B / 4 1C / 3 1D / 1 1E / 1 表 2 航班运行风险统计数据 Table 2 Risk assessment statistical data Risk item Statistical meaning A1 Statistics of captain (non-instructor) and second co-pilot match times A2 Statistics of crew duty time or flight time less than 1 h times A3 Statistics of quick access recorder (QAR) blue and yellow warnings times after flight A4 Statistics times of captain’s experience less than 200 h and the total flight experience less than 3000 h A5 Statistics of the crew first match times to special airport in 12 calendar months A6 Statistics of temporary airborne failures times A7 Number of minimum equipment list (MEL)/configuration discrepancy list (CDL) reservations or MEL reservations affecting near landings A8 Statistics times of flights with special cargo or dangerous goods A9 Statistics of ground handling error times A10 Statistics of low fuel, yaw, abnormal altitude, etc. in flight monitoring A11 Statistics of times the airport failed to meet or on the edge of the weather criteria A12 Statistics of navigation equipment degradations, operating standards and obstacles temporary changes A13 Statistics of special weather (icing, thunderstorms, etc.,) which need to divert A14 Statistics of special operations (polar operations, re-dispatch flights, extended-range operations (ETOPS), extended cross-water, overflying unmanned areas) A15 Statistics of route closure, flow control, height limit, etc. Total risk Comprehensive assessment from the above risk items · 1666 · 工程科学学报，第 42 卷，第 12 期

王岩韬等：基于多变量混沌时间序列的航班运行风险预测模型 1667· 表3航班运行风险时间序列局部数据样本 Table 3 Time series sample data of flight operation risk Sequence number/d A As A13 414 A15 Total risk > o 10 8 8.5 2 6 9 6 9 -.P 6 9 6 3 6 7 3 6 3 5 6 5 2 … … 1096 6 6 5 6 3 5 0 运用坐标延迟重构法如式(1)，则多变量时间序列 重构后的相空间如式(2)： Yj()=[xj(ti).xj(ti-Tj).xj(ti-2r).....xj(i-(m-1))] (1) 200 400600800 1000 Sequence number/d Y(ti)={Y1(t),Y2(t),…,Yn(t)} (2) 图1航班运行总风险时间序列 式中，=max{(m-1)*tjl+1,…,N-1,N. Fig.1 Time series example of flight operation total risk t和m选择是相空间重构的关键.根据Takens 混沌特征.因此，首先通过对总风险时间序列进行 定理，设d为混沌吸引子的维数，若m≥2d+1,则可 单变量混沌预测，得到未来1d的相对误差均值 在m维空间中重构出与风险序列混沌吸引子保持 Mean absolute percentage error,MAPE)21.33%. 微分同胚的轨线1.采用C-C方法计算风险序列 发现单一序列输入，不足以构建精准模型.故后续 的最佳τ和m.以风险时间序列X)为例，计算过 采用多变量时间序列增加输人，通过15个风险源 程如下： 序列预测总风险 (1)设N是风险序列数据的个数，σ是标准差， 并定义为： 2多变量混沌预测模型构建 2.1多变量时间序列相空间重构 C(m.N.M(M-D 2 ∑8-IK0-KD (3) 多变量时间序列相空间重构由单变量相空间 式(3)中，M=N-(m-l)r,g()为heaviside函数. 重构发展而成.以15个风险序列的多变量相空间 重构结果作为预测模型的输入，可使模型输入包 (2》取m=2,345,e=罗e=123,4计算如 含各个风险序列的历史信息 式(4~6)： 对于长度为N的个时间序列x(0=x(1),x(2),, 4 S(m.re,t) (4) x(W,j=1,2.…,n,时间延迟为t和嵌入维为m, 表4时间序列的时间延迟、嵌人维和最大Lyapunov指数 Table 4 Time delay,embedding dimension,and maximum Lyapunov exponent of the time series Embedding Maximum Lyapunov Embedding Maximum Lyapunov Risk item Time delay,r Risk item Time delay,r dimension,m exponent, dimension,m exponent,A A 3 4 0.5440 Ag 9 0.2039 A 4 3 0.9583 A10 J 3 0.9572 A 3 4 0.5973 A 2 8 0.2010 Aa 6 0.9419 Ap 4 0.2723 4 4 2 0.8864 A3 3 4 0.7508 46 4 0.9689 A14 J 4 0.6747 Ar 4 1.1224 As 3 4 0.3707 As 3 6 0.2631 Total risk 4 5 0.3289

混沌特征. 因此，首先通过对总风险时间序列进行 单变量混沌预测，得到未来 1 d 的相对误差均值 （ Mean absolute percentage error, MAPE） 为 21.33%. 发现单一序列输入，不足以构建精准模型. 故后续 采用多变量时间序列增加输入，通过 15 个风险源 序列预测总风险. 2    多变量混沌预测模型构建 2.1    多变量时间序列相空间重构 多变量时间序列相空间重构由单变量相空间 重构发展而成. 以 15 个风险序列的多变量相空间 重构结果作为预测模型的输入，可使模型输入包 含各个风险序列的历史信息. N n Xj(t) = [ xj(1), xj(2),··· , xj(N) ] j = 1,2,··· ,n τj mj 对于长度为 的 个时间序列 ， ，时间延迟为 和嵌入维为 ， 运用坐标延迟重构法如式（1），则多变量时间序列 重构后的相空间如式（2）： Yj(ti) = [xj(ti), xj(ti −τj), x j(ti −2τ),··· , xj(ti −(m−1)τ)] （1） Y(ti) = {Y1(ti),Y2(ti),··· ,Yn(ti)} （2） t 式中， i = max{(mj −1) ∗ τj}+1,··· ,N −1,N. τj mj d m ⩾ 2d +1 m τj mj Xj(t) 和 选择是相空间重构的关键. 根据 Takens 定理，设 为混沌吸引子的维数，若 ，则可 在 维空间中重构出与风险序列混沌吸引子保持 微分同胚的轨线[23] . 采用 C-C 方法计算风险序列 的最佳 和 . 以风险时间序列 为例，计算过 程如下： （1）设 N 是风险序列数据的个数，σ是标准差， 并定义为： C(m,N,r,t) = 2 M(M −1) ∑ 1⩽i⩽j⩽M θ(r −∥X(i)− X(j)∥) （3） 式（3）中， M = N −(m−1)τ，θ(•) 为 heaviside 函数. m = 2,3,4,5 re = eσ 2 （ 2）取 ， , e = 1,2,3,4 ，计算如 式（4～6）： S¯(t) = 1 16 ∑ 5 m=2 ∑ 4 e=1 S (m,re,t) （4） 表 3 航班运行风险时间序列局部数据样本 Table 3 Time series sample data of flight operation risk Sequence number/d A1 A2 A3 A4 A5 …… A13 A14 A15 Total risk 1 7 7 9 7 10 …… 8 10 8 8.5 2 6 7 9 6 9 …… 6 9 6 8 3 4 4 6 4 7 …… 3 6 3 5 4 4 4 6 1 4 …… 2 5 2 3 …… …… …… …… …… …… …… …… …… …… …… 1096 6 6 8 3 6 5 6 3 5 5 表 4 时间序列的时间延迟、嵌入维和最大 Lyapunov 指数 Table 4 Time delay, embedding dimension, and maximum Lyapunov exponent of the time series Risk item Time delay, τ m Embedding dimension, λ Maximum Lyapunov exponent, Risk item Time delay, τ m Embedding dimension, λ Maximum Lyapunov exponent, A1 3 4 0.5440 A9 2 9 0.2039 A2 4 2 0.9583 A10 3 3 0.9572 A3 3 4 0.5973 A11 2 8 0.2010 A4 6 3 0.9419 A12 2 4 0.2723 A5 4 2 0.8864 A13 3 4 0.7508 A6 4 3 0.9689 A14 3 4 0.6747 A7 4 2 1.1224 A15 2 4 0.3707 A8 3 6 0.2631 Total risk 4 5 0.3289 200 400 600 800 1000 0 0 2 4 6 8 10 Sequence number/d Total risk 图 1 航班运行总风险时间序列 Fig.1 Time series example of flight operation total risk 王岩韬等： 基于多变量混沌时间序列的航班运行风险预测模型 · 1667 ·

1668 工程科学学报.第42卷.第12期 1了△Sm,t) △50=4 (5) m=2 Gy= (8) Scor(t)=AS(t)+5(t) (6) 式中，S(m,re,)=C(m,re,t)-Cm(l,re,),△S(m,)=max (S (m,re,t)-min(S (m,re,t)). (9) (3)依据上述计算结果绘制△S()和Sor()曲 线，对于△S()曲线，△S()取得第一极小值时的即 2.3预测模型 为最佳的时间延迟r:对于Scor(0)曲线，Scor()取得 混沌时间序列预测理论中常用方法有：全域 最小值的即为时间延迟窗口Tw,进而根据 法、局域法、基于Lyapunov指数预测法、基于 tw=(Om-1)*T,可计算出最佳的嵌入维m. Volterra自适应滤波器预测法、基于神经网络预测 t和m的合理选择有助于对混沌系统的拟合 法等.此处在选取方法时，一方面总结以往研究， 和预测，表4为15个风险时间序列的最佳延迟时 发现神经网络类方法在混沌预测中表现较为突 间和嵌人维，据此可进行多变量时间序列相空间 出：另一方面，通过查阅文献，筛选出预测结果最 重构，重构后的维度为各个序列的总嵌入维m= 好的几种方法，包含了极限学习机I、RBF神经网 络)、回声状态网络2和Elman神经网络2创等 而后，通过原理分析发现，上述4个模型均由传统 神经网络发展而来.两方面相互印证，决定采用 2.2基于PCA的相空间降维方法 该4种方法，后续分别构建短期预测模型，通过误 由于航班运行风险各指标变量均代表不同含 差分析，提出具有较好适应性的预测方案. 义、混沌重构过程的τ和m;也不尽相同，所以在多 对于给定的N个不同样本(u(,O(),i=1,2,…W, 变量相空间重构后，不同变量在不同历史时期的 其中u()=(,2(①，…，(T∈R为输入层，O()= 观测数据混合在一起，可能造成了信息冗余或噪 [O1(),O2(0,…,O(T∈R为输出层.统一设定隐 声干扰，产生预测误差.采用P℃A抑制此问题 含层神经元个数为h,隐含层神经元与输入层的连接权 PCA基本思想是将一组具有一定相关性的变量重 组为互不相关的主成分，通过保留方差贡献率较 值为w=费，w2,w、与输出层的连接权值 大的主成分代替原变量，实现数据降维，采用 为wg="(0,wz(0,…,w0(,其中j=1,2,…,h. 不同模型的原理如下： PCA改进航班运行风险的预测过程，可以保证： (1)极限学习机ELM ①预测模型的不同输入成分之间相互独立、信息 ELM是单隐层前馈神经网络，如图2.设ELM 不重叠：②输入的状态空间在保证尽量多信息的 隐含层神经元的阈值为b,训练过程中，wn和b随机 情况下实现维度降低. 产生，并保持不变，加快了神经网络的训练速度并 PCA常用的主成分提取方法有特征值分解 使之具有较好的泛化能力 （Eigen value decomposition,EVD)和奇异值分解 (SVD),相对于EVD,SVD在数据重构上误差较 一 g(u(i)),b 1() 小，故本文采用基于SVD的PCA进行降维处理， 2(i)1 g(u(i)),b2 设51≥2≥…≥5g>0为原相空间Y的奇异值，通过 SVD进行分解如式(7)： g(u(i)),ba Y=QDVT (7) Input layer Hidden layer Output layer 式中，Q=[q1,q2,…,qN-rw]为Y的左奇异矩阵，V= 图2极限学习机结构 [yT,T,…,]为r的右奇异矩阵，D=diag(51,52,…,a) Fig.2 ELM structure 假设前y个主成分的累计方差贡献率为Gy,见 图2中，从输入到输出的函数表达式为： 式(8)，一般认为，当满足G,≥85%时，可将y个主 成分构成的状态空间，见式(9)，代替原相空间.实 0(i0= wg"gw·u()+bj） (10) 现降维 l

∆S¯(t) = 1 4 ∑ 5 m=2 ∆S (m,t) （5） S cor(t) = ∆S¯(t)+   S¯(t)    （6） S (m,re,t) = C(m,re,t)−C m(1,re,t) ∆S (m,t) = max {S (m,re,t)−min{S (m,re,t)} 式中， ， . ∆S¯(t) S cor(t) ∆S¯(t) ∆S¯(t) t τj S cor(t) S cor(t) t τw τw = (m−1) ∗ τj mj （ 3）依据上述计算结果绘制 和 曲 线，对于 曲线， 取得第一极小值时的 即 为最佳的时间延迟 ；对于 曲线， 取得 最 小 值 的 即 为 时 间 延 迟 窗 口 ， 进 而 根 据 ，可计算出最佳的嵌入维 . τj mj m = ∑n j=1 mj 和 的合理选择有助于对混沌系统的拟合 和预测，表 4 为 15 个风险时间序列的最佳延迟时 间和嵌入维，据此可进行多变量时间序列相空间 重构，重构后的维度为各个序列的总嵌入维 . 2.2    基于 PCA 的相空间降维方法 τj mj 由于航班运行风险各指标变量均代表不同含 义、混沌重构过程的 和 也不尽相同，所以在多 变量相空间重构后，不同变量在不同历史时期的 观测数据混合在一起，可能造成了信息冗余或噪 声干扰，产生预测误差. 采用 PCA 抑制此问题. PCA 基本思想是将一组具有一定相关性的变量重 组为互不相关的主成分，通过保留方差贡献率较 大的主成分代替原变量 ，实现数据降维. 采用 PCA 改进航班运行风险的预测过程，可以保证： ①预测模型的不同输入成分之间相互独立、信息 不重叠；②输入的状态空间在保证尽量多信息的 情况下实现维度降低. ξ1 ⩾ ξ2 ⩾ ··· ⩾ ξq > 0 Y PCA 常用的主成分提取方法有特征值分解 （ Eigen value decomposition, EVD）和奇异值分解 （SVD），相对于 EVD，SVD 在数据重构上误差较 小，故本文采用基于 SVD 的 PCA 进行降维处理， 设 为原相空间 的奇异值，通过 SVD 进行分解如式（7）： Y = QDVT （7） Q = [q1, q2,··· , qN−τw ] Y V = [v T 1 , v T 1 ,··· , v T m] Y D = diag(ξ1, ξ2,··· , ξd) 式中， 为 的左奇异矩阵， 为 的右奇异矩阵， . γ Gγ Gγ ⩾ 85% γ 假设前 个主成分的累计方差贡献率为 ，见 式（8），一般认为，当满足 时[24] ，可将 个主 成分构成的状态空间，见式（9），代替原相空间，实 现降维. Gγ = ∑ γ i=1 ξi ∑ d i=1 ξi （8） Yˆ = ∑ γ i=1 ξiqiv T i （9） 2.3    预测模型 混沌时间序列预测理论中常用方法有：全域 法 、局域法 、基 于 Lyapunov 指数预测法 、基 于 Volterra 自适应滤波器预测法、基于神经网络预测 法等. 此处在选取方法时，一方面总结以往研究， 发现神经网络类方法在混沌预测中表现较为突 出；另一方面，通过查阅文献，筛选出预测结果最 好的几种方法，包含了极限学习机[18]、RBF 神经网 络[25]、回声状态网络[26] 和 Elman 神经网络[21] 等. 而后，通过原理分析发现，上述 4 个模型均由传统 神经网络发展而来. 两方面相互印证，决定采用 该 4 种方法，后续分别构建短期预测模型，通过误 差分析，提出具有较好适应性的预测方案. N {(u(i),O(i)),i = 1,2,···N} u(i) = [u1(i),u2(i),··· ,ul(i)]T ∈ R l O(i) = [O1(i),O2(i),··· ,Ov(i)]T ∈ R v h w in j = [w in j1 ,w in j2 ,··· ,w in jl] T w out j = [w out j1 (i),w out j2 (i),··· ,w out jv (i)]T j = 1,2,··· ,h 对于给定的 个不同样本 ， 其中 为输入层， 为输出层. 统一设定隐 含层神经元个数为 ，隐含层神经元与输入层的连接权 值为 、与输出层的连接权值 为 ，其中 . 不同模型的原理如下： （1）极限学习机 ELM. b w in b ELM 是单隐层前馈神经网络，如图 2. 设 ELM 隐含层神经元的阈值为 ，训练过程中， 和 随机 产生，并保持不变，加快了神经网络的训练速度并 使之具有较好的泛化能力. 图 2 中，从输入到输出的函数表达式为： O(i) = ∑ h j=1 w out j g(w in j ·u(i)+bj) （10） Input layer Hidden layer Output layer u1 (i) O(i) g(u(i)), b1 g(u(i)), b2 g(u(i)), bh w in w out u2 (i) … … ul (i) 图 2 极限学习机结构 Fig.2 ELM structure · 1668 · 工程科学学报，第 42 卷，第 12 期

王岩韬等：基于多变量混沌时间序列的航班运行风险预测模型 ·1669· 式中，i=1,2,…;只需设定隐含层神经元个数h, 40 并确定激活函数，并通过求解式(11)即可得出 wu,其中，激活函数多选用sigmoid函数 () W=HO (11) 式中，H为H的伪逆，H为隐含层输出矩阵，即H= {g(wn.u(间+b1）,gw·u()+b2),…,g(w7u()+bh)h Input layer Hidden layer Output layer (2)RBF神经网络. 图4回声状态网路结构图 RBF是从三层前向型网络发展而来，不同之 Fig.4 ESN structure 处在于其隐含层的神经元变换函数使用径向基函 连接权值矩阵；为激活函数，常取双曲正切函 数.RBF神经网络的优点在于输入层和隐含层之 数；fout为输出函数，常取恒等函数，此时wou可通 间不需通过权值连接，当RBF中心点确定以后，即 过下式求解： 可确定二者的映射关系，从而使网络结构简单、训 wout=Bo (16) 练简洁 式中，B为B={u(B+1),u(B+2),·,u(W的伪逆， 图3中，R,j=1,2,…,h为激活函数，常使用 高斯函数如式(12)，由此，从输入到输出的函数表 {u(1),u(2),…,(B)1为舍弃的一段初始数据 (4)Elman神经网络 达式为： Elman是一种典型的局部回归网络，属于反馈 R(u(间-c=eala0-c (12) 神经网络，如图5，它与前向神经网络不同之处在 于隐含层中增加了一个承接层，构成局部的反馈 Ok(i)= ∑wx"Ru(-c） (13) 承接层用于接受隐含层前一时刻的输出信号，并 作为下一时刻输入信号的一部分重新输入到隐含 层中，形成一个一步延时算子，使网络具有动态记 忆功能.由于承接层具有可以记忆过去状态的特 0i) 性，非常适合时间序列预测问题 0.(i) Output layer Input layer Hidden layer Output layer Hidden x(i 图3RBF神经网络结构图 laver Fig.3 RBF neural network structure Input layer Connection layer 式(13)中，u()-c为欧式范数，c为高斯函 1) 数的中心，σ为高斯函数的方差，h为隐含层神经元 图5 Elman神经网络结构图 个数，状为隐含层第个神经元到输出层第k个节 Fig.5 Elman neural network structure 点的连接权值，k=1,2,…,y,v为输出层节点个数 (3)回声状态网络 Elman从输人到输出的表达式为： 回声状态网络(ESN)是一种递归神经网络，结 0(i)=fout(woux(i)) (17) 构如图4，与传统神经网络相比，ESN使用大规模 X(i)=f(wiu(i-1)+wexe(i)) (18) 随机连接的递归网络形成一个储备池，取代隐含 X()=X(i-1) (19) 层，通过预先设定储备池权值矩阵的谱半径可保 式中，X(为隐含层神经元向量，X)为承接层的反 证网络的稳定性 馈向量，W心为承接层到隐含层的连接权值，f为隐 ESN从输入到输出的表达式为： 含层神经元的传递函数，fur为输出函数 (i)=fi(wiu(i-1)+WDR(i-1)) (14) 2.4预测流程 0()=fu(wu[u(0,p(]） (15) 根据A1~A1s风险时间序列，对未来pd的总风 式中，p(①为储备池的状态向量，WDR为储备池内部 险值进行短期预测，预测流程如图6：

i = 1,2,···N h w out 式中， ；只需设定隐含层神经元个数 ， 并确定激活函数 ，并通过求解式 (11) 即可得出 ，其中，激活函数多选用 sigmoid 函数. W = H †O （11） H† H H H = {g(w in 1 ·u(i)+b1), g(w in 2 ·u(i)+b2),··· ,g(w in h ·u(i)+bh)} 式中， 为 的伪逆， 为隐含层输出矩阵，即 . （2）RBF 神经网络. RBF 是从三层前向型网络发展而来，不同之 处在于其隐含层的神经元变换函数使用径向基函 数. RBF 神经网络的优点在于输入层和隐含层之 间不需通过权值连接，当 RBF 中心点确定以后，即 可确定二者的映射关系，从而使网络结构简单、训 练简洁. {Rj 图 3 中， , j = 1,2,··· ,h} 为激活函数，常使用 高斯函数如式（12），由此，从输入到输出的函数表 达式为： R(u(i)− cj)=e − 1 2σ2 ∥u(i)−c j∥ 2 （12） Ok(i) = ∑ h j=1 w out jk R(u(i)− cj) （13） u(i)− c j cj σ h w out jk j k k = 1,2,··· , v v 式（13）中， 为欧式范数， 为高斯函 数的中心， 为高斯函数的方差， 为隐含层神经元 个数， 为隐含层第 个神经元到输出层第 个节 点的连接权值， ， 为输出层节点个数. （3）回声状态网络. 回声状态网络（ESN）是一种递归神经网络，结 构如图 4，与传统神经网络相比，ESN 使用大规模 随机连接的递归网络形成一个储备池，取代隐含 层，通过预先设定储备池权值矩阵的谱半径可保 证网络的稳定性. ESN 从输入到输出的表达式为： φ(i) = f in(w inu(i−1)+WDRφ(i−1)) （14） O(i) = f out(w out[u(i),φ(i)]) （15） 式中， φ(i) 为储备池的状态向量， WDR 为储备池内部 f in f out w out 连接权值矩阵； 为激活函数，常取双曲正切函 数； 为输出函数，常取恒等函数，此时 可通 过下式求解： w out = B †O （16） B † B={u(β+1),u(β+2),··· ,u(N)} {u(1),u(2),··· ,u(β)} 式 中 ， 为 的 伪 逆 ， 为舍弃的一段初始数据. （4）Elman 神经网络. Elman 是一种典型的局部回归网络，属于反馈 神经网络，如图 5，它与前向神经网络不同之处在 于隐含层中增加了一个承接层，构成局部的反馈. 承接层用于接受隐含层前一时刻的输出信号，并 作为下一时刻输入信号的一部分重新输入到隐含 层中，形成一个一步延时算子，使网络具有动态记 忆功能. 由于承接层具有可以记忆过去状态的特 性，非常适合时间序列预测问题. Elman 从输入到输出的表达式为： O(i) = f out(w outχ(i)) （17） χ(i) = f(w inu(i−1)+Wcχ c (i)) （18） χ c (i) = χ(i−1) （19） χ(i) χ c (i) Wc f f out 式中， 为隐含层神经元向量， 为承接层的反 馈向量， 为承接层到隐含层的连接权值， 为隐 含层神经元的传递函数， 为输出函数. 2.4    预测流程 根据 A1～A15 风险时间序列，对未来 p d 的总风 险值进行短期预测，预测流程如图 6： u1 (i) O1 R (i) 1 R2 Rh O2 (i) Ov (i) u2 (i) … … … ul (i) Input layer Hidden layer Output layer 图 3 RBF 神经网络结构图 Fig.3 RBF neural network structure Input layer Output layer Hidden layer O(i) u1 (i) u2 (i) … ul (i) w in w out 图 4 回声状态网络结构图 Fig.4 ESN structure Output layer Hidden layer Input layer Connection layer O(i) χ(i) χ c (i) u(i−1) … … 图 5 Elman 神经网络结构图 Fig.5 Elman neural network structure 王岩韬等： 基于多变量混沌时间序列的航班运行风险预测模型 · 1669 ·

.1670 工程科学学报.第42卷，第12期 Multivariable chaotic time series of flight operation risks Iterative predictor Chaos Phase space reconstruction Phase space analysis reconstruction results PCA analysis Dimensionality reduction processing based on PCA Iterative predictor Prediction model Parameter ELM RBF optimization Prediction about t+1 methods ESN Elman F+1 N SVM-based risk assessment Iterative prediction results Short-term prediction results and accuracy analysis 图6预测流程 Fig.6 Prediction process 1)根据风险项A1~A15的混沌特征分析结果， 其中，将前900d作为训练集，剩余196d作为测试 确定延迟时间和嵌入维，进行多变量相空间重构； 集.经过多变量相空间重构后再使用PCA降维， 2)对相空间进行PCA分析，确定前若干个主 各主成分的方差贡献率如表5. 成分构成的状态空间，实现相空间的降维： 根据PCA原理，PCA处理后得到的主成分应 3)根据迭代预测方法，分别构建ELM、RBF、 能解释原始变量尽量多的信息，累计方差贡献率 ESN和Elman风险短期预测.假设当前时刻为t, 达到85%以上为宜：同时，当变量之间相关性较强 在迭代预测过程，首先以第天降维后的状态空间 的时候，少数几个主成分的累积贡献率就能达到 作为预测模型的输人，预测出第+1天A1~A15的 较高水平，当变量之间相关性较弱的时候，则需较 风险值；然后，利用第+1天的预测值，重新进行多 多的主成分才能达到较高的累积贡献率，但更为 变量相空间重构并作降维处理，再将此时降维后 重要的是保证计算后达到有效的降维效果.因此， 的状态空间输入至预测模型中，预测出第+2天 在较高的累计方差贡献率基础上结合碎石图拐点 A1~A1s的风险值.以此类推，直至完成第+p天的 出现位置，再综合选择.通过采用85%、90%和 预测，得到A1~A1s的pd预测结果 95%不同阈值，经上文预测流程计算得到，当保留 4)采用支持向量机(SVM)进行风险评估，得 前31个主成分累计方差贡献率为90%时，此处预 到总风险值的pd预测值，通过与总风险的真实值 测结果最为理想，即维数由62维压缩至31维 比较，计算预测误差并分析不同模型的预测精度. 3.2参数优化 3实例分析 根据训练集数据，采用交叉验证法、遗传算法 和粒子群算法对预测模型的参数进行优化.但由 3.1主成分降维 于后续预测过程采用多步迭代方式，过程复杂且 将15个风险源时间序列数据作为预测样本， 易受预测步长的影响，因此假设多步迭代预测最 表5主成分分析的方差贡献率 Table 5 Partial variance contribution rate of PCA Principal Variance Accumulative contribution Principal Variance Accumulative contribution component contribution/% rate/% component contribution/% rate/% 1 15.5907 15.5907 30 0.8018 89.8892 2 11.5223 27.1131 31 0.7893 90.6785 3 7.5821 34.6952 32 0.7054 91.3839 4 6.3632 41.0583 62 0.0092 100.0000

1）根据风险项 A1～A15 的混沌特征分析结果， 确定延迟时间和嵌入维，进行多变量相空间重构； 2）对相空间进行 PCA 分析，确定前若干个主 成分构成的状态空间，实现相空间的降维； ti ti ti +1 A1 A15 ti +1 ti +2 A1 A15 ti + p A1 A15 p 3）根据迭代预测方法，分别构建 ELM、RBF、 ESN 和 Elman 风险短期预测. 假设当前时刻为 ， 在迭代预测过程，首先以第 天降维后的状态空间 作为预测模型的输入，预测出第 天 ～ 的 风险值；然后，利用第 天的预测值，重新进行多 变量相空间重构并作降维处理，再将此时降维后 的状态空间输入至预测模型中，预测出第 天 ～ 的风险值. 以此类推，直至完成第 天的 预测，得到 ～ 的 d 预测结果. p 4）采用支持向量机（SVM）进行风险评估，得 到总风险值的 d 预测值，通过与总风险的真实值 比较，计算预测误差并分析不同模型的预测精度. 3    实例分析 3.1    主成分降维 将 15 个风险源时间序列数据作为预测样本， 其中，将前 900 d 作为训练集，剩余 196 d 作为测试 集. 经过多变量相空间重构后再使用 PCA 降维， 各主成分的方差贡献率如表 5. 根据 PCA 原理，PCA 处理后得到的主成分应 能解释原始变量尽量多的信息，累计方差贡献率 达到 85% 以上为宜；同时，当变量之间相关性较强 的时候，少数几个主成分的累积贡献率就能达到 较高水平，当变量之间相关性较弱的时候，则需较 多的主成分才能达到较高的累积贡献率，但更为 重要的是保证计算后达到有效的降维效果. 因此， 在较高的累计方差贡献率基础上结合碎石图拐点 出现位置，再综合选择. 通过采用 85%、90% 和 95% 不同阈值，经上文预测流程计算得到，当保留 前 31 个主成分累计方差贡献率为 90% 时，此处预 测结果最为理想，即维数由 62 维压缩至 31 维. 3.2    参数优化 根据训练集数据，采用交叉验证法、遗传算法 和粒子群算法对预测模型的参数进行优化. 但由 于后续预测过程采用多步迭代方式，过程复杂且 易受预测步长的影响，因此假设多步迭代预测最 表 5 主成分分析的方差贡献率 Table 5 Partial variance contribution rate of PCA Principal component Variance contribution/% Accumulative contribution rate/% Principal component Variance contribution/% Accumulative contribution rate/% 1 15.5907 15.5907 30 0.8018 89.8892 2 11.5223 27.1131 31 0.7893 90.6785 3 7.5821 34.6952 32 0.7054 91.3839 4 6.3632 41.0583 … … … … … … 62 0.0092 100.0000 Multivariable chaotic time series of flight operation risks Phase space reconstruction Chaos analysis results PCA analysis Iterative predictor ELM RBF ESN Elman SVM-based risk assessment Short-term prediction results and accuracy analysis Parameter optimization methods Iterative predictor Dimensionality reduction processing based on PCA Prediction model Iterative prediction results N Y Phase space reconstruction Prediction about ti+1 t t i>p i=ti+1 图 6 预测流程 Fig.6 Prediction process · 1670 · 工程科学学报，第 42 卷，第 12 期

王岩韬等：基于多变量混沌时间序列的航班运行风险预测模型 1671 优的前提为单步预测结果最优，综合衡量不同参 数优化方法和多次试验的输出，参数优化结果如 表6所示 (a) 901 940 980 1020 10601090 表6预测模型的参数优化结果 Sequence number/d Table 6 Parameter optimization results of the prediction models Parameter Value Parameter Value Hidden layer sizes of ELM 75 Leaking rate of ESN 0.80 b Spread of RBF Layer delays of Elman 4 903 940 980 1020 1060 1092 Sequence number/d Reservoir sizes of ESN 18 Hidden layer sizes of Elman 152 10 3.3短期预测的测试结果 6 结合航空公司实际需要，设定预测范围为未 (c) 来7d,以第900天为起始点，向后预测得到第 905 940 9801020 10601094 901~907天风险值，如图7.但仅从一组数据无法 Sequence number/d -Expected results-e-ELM results 说明模型预测精度，再从第901天至第1090天依 -RBF results ESN results-*-Elman results 次预测，总计190次，每次均得到7d预测结果.统 图8预测结果示例.(a)第1天：(b)第3天：(c)第5天 计所有预测的第1、3、5天结果如图8：其中，图8(a) Fig.8 Example of prediction results:(a)day 1;(b)day 3;(c)day 5 均为第1天预测结果，其范围为[901,1090：图8(b) 降维处理有助于提高RBF的预测精度 是第3天结果，范围为903,1092：图8(c)是第5天 3.5误差分析 结果，范围为[905,1094，图8(a)(b)(c)均包含 从表7可见，不同模型预测相对误差的频数分 190个数据 布均呈右偏，若以均值反映相对误差的集中趋势， 10 结果会偏大.因此，为使相对误差均值(MAPE)更 具代表性，采用修正的MAPE值，即去掉最大 6 -Expected results -e-ELM results -e-RBF results 5%和最小5%的数据后所求得的MAPE值，结果 ESN results 见图9~10. 0 -*-Elman results 94 896 898900902904 906 从图9可见，当训练样本为300时，误差较大； Sequence number/d 随着训练样本增加，误差开始下降：当训练样本从 因7从第900天向后预测样例 300增加至500时.7d内的修正MAPE值平均下 Fig.7 Example of prediction results from the 900th day 降5.48%，预测效果明显提升；当训练样本从 3.4结果分析 700增加至900时，修正MAPE值平均下降1.27%， 计算ELM、RBF、ESN和Elman等不同模型短 效果提升开始放缓：当训练样本从800增加至 期预测结果的相对误差.统计相对误差(Relative 900时，修正MAPE值平均仅下降0.21%，可视为 eror,RE)的频数分布，并以降维前的预测作对比， 无明显变化：因此，基于现有数据情况，后续使用 如表7. 训练样本为900进行计算. 从表8可见，RBF降维后的预测效果最佳，对 图10是基于降维后的RBF预测结果，修正 于RE50%的异常结果出现次数较少，4种预测模型 用要求，当预测相对误差低于20%时，结果具有实 均可行有效.进一步对比降维前后结果发现，降维 践操作价值，因此该方法所得未来5d预测结果均 后RE<50%占比相对于降维前提高了约6%，说明 可作为依据进而指导生产实践

优的前提为单步预测结果最优，综合衡量不同参 数优化方法和多次试验的输出，参数优化结果如 表 6 所示. 3.3    短期预测的测试结果 结合航空公司实际需要，设定预测范围为未 来 7 d， 以 第 900 天为起始点 ，向后预测得到 第 901～907 天风险值，如图 7. 但仅从一组数据无法 说明模型预测精度，再从第 901 天至第 1090 天依 次预测，总计 190 次，每次均得到 7 d 预测结果. 统 计所有预测的第 1、3、5 天结果如图 8；其中，图 8（a） 均为第 1 天预测结果，其范围为 [901，1090]；图 8（b） 是第 3 天结果，范围为 [903，1092]；图 8（c）是第 5 天 结果 ，范围为 [905， 1094]，图 8（ a） （ b） （ c）均包含 190 个数据. 3.4    结果分析 计算 ELM、RBF、ESN 和 Elman 等不同模型短 期预测结果的相对误差. 统计相对误差（Relative error, RE）的频数分布，并以降维前的预测作对比， 如表 7. 从表 8 可见，RBF 降维后的预测效果最佳，对 于RE50% 的异常结果出现次数较少，4 种预测模型 均可行有效. 进一步对比降维前后结果发现，降维 后 RE<50% 占比相对于降维前提高了约 6%，说明 降维处理有助于提高 RBF 的预测精度. 3.5    误差分析 从表 7 可见，不同模型预测相对误差的频数分 布均呈右偏，若以均值反映相对误差的集中趋势， 结果会偏大. 因此，为使相对误差均值（MAPE）更 具代表性 ，采用修正 的 MAPE 值 ，即去掉最 大 5% 和最小 5% 的数据后所求得的 MAPE 值，结果 见图 9～10. 从图 9 可见，当训练样本为 300 时，误差较大； 随着训练样本增加，误差开始下降；当训练样本从 300 增加至 500 时 ，7 d 内的修正 MAPE 值平均下 降 5.48%， 预 测 效 果 明 显 提 升 ； 当 训 练 样 本 从 700 增加至 900 时，修正 MAPE 值平均下降 1.27%， 效果提升开始放缓 ；当训练样本 从 800 增 加 至 900 时，修正 MAPE 值平均仅下降 0.21%，可视为 无明显变化；因此，基于现有数据情况，后续使用 训练样本为 900 进行计算. 图 10 是基于降维后的 RBF 预测结果 ，修正 MAPE 值随预测周期增长而逐步递增，从第 1 天至 第 7 天平均每天增加 2.35%，呈现误差累积现象. 第 1 天修正 MAPE 值仅为 11.32%，精度较高 ；第 5 天 MAPE 值仍为 18.21%，其预测精度和预测周 期均优于文献 [12] 和 [13]. 根据国内航空公司使 用要求，当预测相对误差低于 20% 时，结果具有实 践操作价值，因此该方法所得未来 5 d 预测结果均 可作为依据进而指导生产实践. 表 6 预测模型的参数优化结果 Table 6 Parameter optimization results of the prediction models Parameter Value Parameter Value Hidden layer sizes of ELM 75 Leaking rate of ESN 0.80 Spread of RBF 1 Layer delays of Elman 4 Reservoir sizes of ESN 18 Hidden layer sizes of Elman 152 894 896 898 900 902 904 906 Sequence number/d 0 2 4 6 8 10 Risk Expected results ELM results RBF results ESN results Elman results 图 7 从第 900 天向后预测样例 Fig.7 Example of prediction results from the 900th day 940 980 1020 1060 0 2 4 6 8 10 (a) (b) (c) 901 1090 Sequence number/d 903 940 980 1020 1060 1092 Sequence number/d 905 940 980 1020 1060 1094 Sequence number/d 0 2 4 6 8 10 0 2 4 6 8 10 Risk Risk Risk Expected results ELM results RBF results ESN results Elman results 图 8 预测结果示例. （a）第 1 天；（b）第 3 天；（c）第 5 天 Fig.8 Example of prediction results: (a) day 1; (b) day 3; (c) day 5 王岩韬等： 基于多变量混沌时间序列的航班运行风险预测模型 · 1671 ·

·1672 工程科学学报.第42卷，第12期 表7各模型降维前后预测相对误差频数 Table 7 RE frequency of each prediction model before and after dimension reduction RE frequency of prediction model before dimension reduction RE frequency of prediction model after dimension reduction Model Prediction range 50% 50% Day 1 131168.95% 41/21.58% 1819.47% 127166.84% 43/22.63% 20/10.53% ELM Day3 107156.32% 55/28.95% 28/14.74% 109/57.37% 46/24.21% 35/18.42% Day 5 115/60.53% 45/23.68% 30/15.79% 113/59.47% 46/24.21% 31/16.32% Day 1 157/82.63% 22/11.58% 11/5.79% 151/79.479% 27/14.21% 12/6.329% RBF Day 3 150/78.95% 20/10.53% 20/10.53% 138/72.63% 32/16.84% 20/10.53% Day 5 143/7526% 19/10.00% 21/11.05% 129/67.89% 28/14.74% 24/12.63% Day 1 111158.42% 53127.89% 26/13.68% 108/56.84% 56/29.47% 26/13.68% ESN Day 3 102/53.68% 58/30.53% 30115.79% 99152.11% 61/32.11% 30/15.79% Day 5 106155.79% 52/27.37% 32/16.84% 105/55.26% 46/24.21% 39/20.53% Day 1 142/74.74% 32/16.84% 16/8.42% 142/74.74% 35/18.42% 13/6.84% Elman Day3 116/61.05% 47/24.74% 27/14.21% 129/67.89% 31/16.329% 30/15.79% Day 5 109/57.37% 54/28.42% 27/14.21% 111/58.42% 47124.749% 32/16.84% 300 samples 500 samples 至第5天MAPE值仍为18.21%，说明未来5d的预 600 samples …700 samples --◆--800 samples -900 samples 测结果均可满足航空公司实践操作要求.从预测 35 30 精度和周期两个角度综合评判，多变量时间序列 方法对航班运行风险的预测可行且有效 国内航空公司以周为单位制定航班计划，因 10 此可在文中结果基础上，最早于5d前从飞机调配 5 和人员调度等方面调整计划编排，提前降低运行 0 Day 1 Day 2 Day 3 Day 4 Day 5 Day 6 Day 7 风险.该方案还可扩展至其他领域，适用于结果受 Prediction range 多变量综合影响的、各变量为混沌时间序列的短 图9训练样本数量与预测精度 期预测问题，其预测周期随多变量的混沌时间序 Fig.9 Number of training samples and prediction accuracy 列特性而定 30 25.44 25 2284 参考文献 15.1416.35 18.21 [Feng Z L.Promote the high-quality development of civil aviation 15 11.32 13.42 with the new development concept.People's Forum,2019(5):6 10 8 (冯正霖.以新发展理念引领新时代民航高质量发展.人民论坛 2019(5):6) 0 Day I Day 2 Day 3 Day 4 Day 5 Day 6 Day 7 [2]Bowen B D,Headley D E,Luedtke J R.A quantitative Prediction range methodology for measuring airline quality.J Aviation/Aerospace 图10RBF降维后的多变量预测精度 Educ Res,.1992,2(2):27 Fig.10 Prediction accuracy of the RBF model after dimension reduction [3]Lower M.Magott J,Skorupski J.Analysis of air traffic incidents using event trees with fuzzy probabilities.Fy Sets Syst,2016, 4结论 293(C:50 [4] Kokangu A,Polat U,Dagsuyu C.A new approximation for risk 通过混沌识别和相空间重构，构建的4种风险 assessment using the AHP and Fine Kinney methodologies.Saf 预测模型在迭代预测后发现：在降维前后的ELM、 Sci,2017,91:24 RBF、ESN和Elman预测模型中，降维后的RBF预 [5]FAA.Flight risk assessment tool[EB/OL].Federal Aviation Administration Online (2007-08-02)[2018-08-31].https://www.faa 测模型效果最佳；其未来第1天预测结果RE<25% gov/news/safety_briefing/2016/media/SE_Topic_16-12.pdf 的频数可达到82.62%，修正MAPE值仅为11.32%； [6]NATA.NATA 's risk assessment tool -RA Check[EB/OL]

4    结论 通过混沌识别和相空间重构，构建的 4 种风险 预测模型在迭代预测后发现：在降维前后的 ELM、 RBF、ESN 和 Elman 预测模型中，降维后的 RBF 预 测模型效果最佳；其未来第 1 天预测结果 RE50% 50% ELM Day 1 131 / 68.95% 41 / 21.58% 18 / 9.47% 127 / 66.84% 43 / 22.63% 20 / 10.53% Day 3 107 / 56.32% 55 / 28.95% 28 / 14.74% 109 / 57.37% 46 / 24.21% 35 / 18.42% Day 5 115 / 60.53% 45 / 23.68% 30 / 15.79% 113 / 59.47% 46 / 24.21% 31 / 16.32% RBF Day 1 157 / 82.63% 22 / 11.58% 11 / 5.79% 151 / 79.47% 27 / 14.21% 12 / 6.32% Day 3 150 / 78.95% 20 / 10.53% 20 / 10.53% 138 / 72.63% 32 / 16.84% 20 / 10.53% Day 5 143 / 75.26% 19 / 10.00% 21 / 11.05% 129 / 67.89% 28 / 14.74% 24 / 12.63% ESN Day 1 111 / 58.42% 53 / 27.89% 26 / 13.68% 108 / 56.84% 56 / 29.47% 26 / 13.68% Day 3 102 / 53.68% 58 / 30.53% 30 / 15.79% 99 / 52.11% 61 / 32.11% 30 / 15.79% Day 5 106 / 55.79% 52 / 27.37% 32 / 16.84% 105 / 55.26% 46 / 24.21% 39 / 20.53% Elman Day 1 142 / 74.74% 32 / 16.84% 16 / 8.42% 142 / 74.74% 35 / 18.42% 13 / 6.84% Day 3 116 / 61.05% 47 / 24.74% 27 / 14.21% 129 / 67.89% 31 / 16.32% 30 / 15.79% Day 5 109 / 57.37% 54 / 28.42% 27 / 14.21% 111 / 58.42% 47 / 24.74% 32 / 16.84% 0 5 10 15 20 25 30 35 Day 1 Day 2 Day 3 Day 4 Day 5 Day 6 Day 7 Modified MAPE/ % Prediction range 300 samples 500 samples 600 samples 700 samples 800 samples 900 samples 图 9 训练样本数量与预测精度 Fig.9 Number of training samples and prediction accuracy Prediction range 11.32 13.42 15.14 16.35 18.21 22.84 25.44 0 5 10 15 20 25 30 Day 1 Day 2 Day 3 Day 4 Day 5 Day 6 Day 7 Modified MAPE/ % 图 10 RBF 降维后的多变量预测精度 Fig.10 Prediction accuracy of the RBF model after dimension reduction · 1672 · 工程科学学报，第 42 卷，第 12 期