郑 磊，池 宏，邵雪焱

(1．中国科学院科技战略咨询研究院，北京 100190；2．中国科学院大学，北京 100049)

基于QAR数据的飞行操作模式及其风险分析

郑 磊1,2，池 宏1，邵雪焱1

(1．中国科学院科技战略咨询研究院，北京 100190；2．中国科学院大学，北京 100049)

QAR海量数据分析是航空安全管理的一个重要课题，准确而高效地挖掘其中蕴含的飞行操作模式、评估其风险水平，是切实提高飞行操作品质的关键。然而，飞行过程中飞行员需要根据动态变化的环境条件、飞机状态，对油门、杆、盘、舵等设备进行操作，是一个需要不断调整、相互配合的过程，这增加了对飞行操作模式及其风险进行分析的难度。本文从飞行参数多元时间序列数据的特征提取及相似性度量入手，通过聚类挖掘蕴含的飞行操作模式，进而，通过分析飞行操作模式与QAR监控指标值之间的关联关系，量化不同飞行操作模式的风险水平。最后，本文以首都机场特定机型着陆阶段QAR实际数据为例，验证了模型的有效性。

QAR；多元时间序列；操作模式；风险分析

1 引言

伴随飞行数据记录设备的不断更新升级，飞行数据的采集数量与质量都有很大提高。快速存取记录器(Quick access recorder，QAR)记录了飞行全过程发动机等设备状况、外部环境状况，以及飞行员的操作状况，对QAR海量数据的深入挖掘分析日益受到航空监管机构及航空运营商的高度重视，成为事故调查、机务维修、飞行品质监控等多项工作的有力支撑。NASA支持研发的航空性能测量系统(Aviation Performance Measuring System，APMS)主要围绕飞行数据分析的先进理念和软件原型开展工作。

运用QAR数据进行飞行品质分析也受到国内外学者的关注，曹海鹏等[1]利用神经网络对QAR超限事件进行诊断分析，并以某航空公司B737机型实际数据对模型进行验证，结果显示基于神经网络的模型不仅能够判断航空器是否重着陆，还能够分析导致重着陆的原因，同时训练好的网络执行速度快，可以满足实时诊断的需求。祁明亮等[2]和Shao Xueyan等[3]以着陆阶段的一类QAR超限事件为研究对象，建立数学规划模型，对操作参数的风险区域进行分析。曹海鹏[4]研究了使用灰色预测模型对某超限事件进行时间序列预测，此外，还与运用指数平滑法的计算结果进行了对比。结果显示，灰色预测的效果好于指数平滑。张荣家等[5]利用QAR数据研究了A321型飞机在着陆阶段易引发超限事件的原因，并提出安全操作建议。孙瑞山和韩文律[6]利用差异检验对QAR超限事件及其飞行记录数据进行分析，通过对爬升速度大(35-1000ft)，下降率大(2000-1000ft)及接地点远三类超限事件的分析，得到K-W检验能够用于分析飞行参数之间的分布差异。孙瑞山和肖亚兵[7]结合QAR飞行记录，从人的行为特征出发，确定了描述飞行员操作的完整性、平稳性、准确性、及时性的操作特征向量指标。Li Lishuai[8]运用主成分分析对起飞和降落阶段的多元时间序列飞行参数先进行降维，再进行分析。文章将时间序列中的变量视为随机变量，是一种统计意义上的度量方法，忽略了观察值的时序关系。综上，目前利用QAR数据进行飞行品质研究主要使用的是截面数据，即飞行参数的某些状态值，对于飞行参数多元时序数据隐含规律的挖掘亟待加强。

多元时序数据挖掘也是数据挖掘领域的研究热点，其中的多元时间序列相似性研究方面积累了相当数量的成果，包括基于统计特征的方法、基于符号的表示方法、动态时间规整方法(Dynamic Time Warping，DTW)、奇异值分解方法(Singular Value Decomposition，SVD)、基于点分布特征的方法，以及集成方法等。闫相斌等[9]和Cai Qinglin等[10]提出使用时间序列的统计特征如趋势、平均值、方差等作为聚类特征向量，该方法对序列的间隔与长度无特殊的要求，具有较好的适应性，但由于采用的都是全局特征描述变量，难以细致刻画时间序列的局部特征。徐梅和黄超[11]引入符号时间序列分析方法，将收益时间序列转化为符号序列，由符号序列中不同的字代表不同的收益变化模式，该方法可以极大的减小数据的规模，但是产生的序列值比较离散，无法精确比较序列的相似性。李正欣等[12]、吴虎胜等[13]和Keogh等[14]都是基于动态时间规整提出的时间序列数据分析模型，都采用了对时间序列数据进行分段拟合，并基于拟合后的分段计算DTW距离进行相似性分析。李正欣等[12]采用自下而上的分段方式，将长度为n的多元时间序列分成n/2段，递归地合并拟合误差最小的相邻段，直至任意两相邻段合并的拟合误差都大于给定的最大误差阈值。在定义相似性分析的距离时，时间跨度和角度的差异均作为权重进行计算。吴虎胜等[13]采用滑窗的分段方式，从初始时刻开始，每次加入一个时间单位跨度的数据点，计算窗口内数据点的直线拟合误差，当误差大于某个阈值时进行切分。在定义相似性分析的距离时，以时间跨度作为权重进行计算。Keogh等[14]较早的提出使用分段拟合后的时间序列来计算DTW距离。毛红保等[15]引入奇异值分解(Singular Value Decomposition，SVD)方法，采用相关系数矩阵作为特征提取的基础，依据线性空间中坐标变换原理构造模式匹配模型。该方法仅使用提取的奇异值信息，未使用特征矩阵信息，存在误判风险且计算复杂度较高。管河山等[16]使用基于点分布(Point Distribution，PD)特征的方法抽取局部重要点刻画多元时间序列特征，并根据重要点的分布特征构建模式匹配模型，该方法处理小规模多元时间序列性能较好。ten Holt等[17]提出MD-DTW距离，将动态时间规整方法扩展到多元时间序列的相似性研究中。DTW方法适用于时间轴的伸缩，定义了一种允许对时间轴进行伸缩下的最佳相似关系，但其计算复杂度高，在用于海量时间序列数据分析时存在困难。Bankó等[18]提出将基于均值的DTW距离和基于角度的DTW距离结合起来的混合DTW距离，可以较好的提升DTW算法的效率，减少运算时间，同时可以获得满意的运算结果。综上，围绕具体分析问题，在多元时序数据的特征提取，以及相似性度量方面，需要研究有针对性的方法。

2 问题描述

在飞行过程中，飞行员需要根据动态的环境条件、飞机状态，综合操作油门、杆位等设备，实现对飞机高度、速度、姿态等的控制。QAR记录了飞行全过程的环境条件参数、飞机状态参数，以及飞行员的操作参数，是一类多元时间序列，可以表示为时间序列X=[x1,x2,…,xn],每个元素xi是由m个飞行参数组成的向量，长度为n的时间序列记为矩阵X(X∈Rm×n)。

通过监控QAR数据，可以发现飞行员何时加油门、何时收襟翼、何时收轮，进而与规章、手册的要求进行比对。对于飞行员操作中存在的共性特征，譬如，“抬头慢”、“人工干预自动油门”，还需要进一步对此时的重点安全监控参数，如飞机的着陆垂直载荷等数值进行分析，从而量化飞行操作模式的风险大小。

综上，本文主要研究两方面问题：一是，飞行参数多元时间序列数据的相似性度量，即从多元时序数据中挖掘相似的飞行操作模式；二是，分析飞行操作模式与QAR监控事件之间的关联关系，量化不同操作特征的风险水平。

3 飞行操作模式分析

研究飞行操作模式，目的是发现具有代表性的飞行员操作手法，即对油门、杆位、盘、舵等设备的操作，是进一步针对该类操作手法进行安全风险分析的基础，这里，有代表性的操作手法是指具有一定数量的飞行过程中出现的相似操作，因此，我们将飞行操作模式的分析归结为对飞行参数多元时间序列数据相似性的研究。分析多元时间序列的相似性需要解决两个关键问题：一是，如何抽象和概括时间序列的特征；二是，如何度量时间序列的相似性。

3.1飞行参数多元时序数据的特征提取

QAR可记录几十甚至上百飞行小时的数据，以波音747飞机ARINC717(128WPS)为例，其QAR原始代码的帧结构为1帧4秒，1帧包括4副帧，每副帧1秒，包括128个数据字[19]，因此，QAR的数据量是非常庞大的，对QAR数据的挖掘技术提出了很高的要求。时间序列一般具有趋势性、季节性、周期性等特征，针对本文所要研究的飞行操作模式，重点关注QAR数据的趋势性特征。由于QAR数据的海量性、复杂性和噪音干扰等特点，直接在原始数据上进行分析计算量大，还可能影响算法的准确性和可靠性，因此，可采用时间序列模式表示方法对飞行参数进行处理，不仅可大大压缩数据量，还可有效剔除原始参数序列中的噪声干扰[20]。

常用的序列模式表示方法有离散傅立叶变换、离散小波变换、奇异值分解，以及分段线性，考虑本文在研究飞行操作模式时，希望直观反映时间轴上飞行参数的趋势特征，故采用分段线性的方法，研究飞行参数时序数据的特征提取。在对飞行参数多元时序数据进行分段线性拟合时，选取滑动时间窗算法，该算法的基本思想是用一段直线拟合同一窗口内的数据点，当窗口内数据的拟合误差小于给定阈值时，将新数据加入，增加窗口宽度，直到累积的拟合误差大于给定阈值时，开启新的窗口。

在分段基础上，采用拟合线段的倾斜角与窗口的时间跨度刻画时序数据的趋势特征，定义Aij表示飞行参数i在窗口j内拟合线段的倾斜角，tj表示窗口j的时间跨度，则对于包含m个飞行参数，时间长度为n的多元时间序列，可采用如下矩阵表示：

(1)

其中，i=1,2,…m；j=1,2,…,Q，Q为分段后的时间窗数目。

3.2飞行参数多元时序数据的相似性度量

如前所述，在多元时间序列相似性研究方面形成了动态时间规整、奇异值分解、主成分分析，以及基于点分布特征等方法。本文所要研究的飞行操作模式，重点关注QAR数据的趋势特征，故采用动态时间规整(DTW)方法进行研究。该方法被广泛应用于在语音识别、手势识别等问题。设有时间序列C和Q，如图1(a)，DTW方法的处理结果如图1(b)所示，即找到两条时间序列曲线上每个点的对应点，使得C和Q上对应点对之间距离的和最小。

图1 动态时间规整方法原理示意

设时间序列X=[x1，x2,…,xn]和Y=[y1，y2,…,yn],则两个序列的DTW距离Ddtw(X,Y)定义如下式，其中，Dbase(xi,yj)表示向量点xi和yj之间的基距离，可以根据情况选择不同的距离度量。

Ddtw(X,Y)=Dbase(x1,y1)+

(2)

本文中对于飞行参数时序数据进行分段线性拟合，并提取拟合线段的倾斜角角度及时间跨度刻画时序数据特征，因此，在采用DTW方法进行飞行参数多元时序数据相似性分析时借鉴文章[13]，定义基距离如下：

(3)

其中，Dbase(xp,yq)表示列向量xp和yq之间的基距离,tX和tY分别为多元时间序列X和Y的时间跨度。通过定义Dbase(xp,yq)将多维时间序列降为一维时间序列。

3.3飞行操作多元时序数据的聚类分析

在飞行参数多元时序数据的相似性度量基础上，本文采用聚类的方法挖掘有代表性的飞行操作手法，即有代表性的飞行操作模式。由于本文对于飞行员操作模式的分析，并没有关于模式个数的经验信息，故采用系统聚类进行研究，类间距的计算选取类平均法。

(4)

3.4案例分析

根据波音公司2006-2015年的统计数据[21]，着陆阶段虽然在时间上只占整个航程的1%左右，但发生的重大事故占全部航段重大事故的20%左右，因此，本文采集着陆前20秒的飞行参数数据，并以油门(单位，%RPM)、杆位(单位，DEG)两项参数QAR记录为例进行研究。

图变化情况

表1 系统聚类结果

考虑到拟合后的轨迹较好地展现了整体趋势，本文采用拟合后的轨迹对类操作特征进行研究，如图3。本文主要以油门和杆的操作进行分析，在数据筛选时尽可能排除风速、风向对飞机着陆操作的影响。飞机着陆过程需要收油门减速，飞机的高度和速度不断下降，飞行员根据飞机速度和高度的变化情况控制操纵杆。当飞机下降到特定高度标准时，就要开始拉杆，即开始拉平，使飞机下降的速度逐渐变缓，再达到特定高度标准后拉平，飞机平飘，之后速度降低飞机慢慢着陆。

在本文样本机型的《飞行训练手册》中，对于着陆阶段的操作有如下描述：1)3度进近下滑轨迹；2)拉平距离跑道头大约1000-2000英尺；3)典型的着陆拉平时间范围在4-8秒，并且随进近速度而变化；4)飞机机身姿态基于典型的着陆重量，襟翼30，VREF30+5节(进近)和VREF30+0节(接地)，每高于此速度5节应减小1度；5)主起落架的跑道入口高度和飞行员的目视高度参考两排杆/三排杆VASI着陆几何图。

在图3中，(a)表示在落地前约7-8秒开始收油门，先缓慢拉杆，在落地前约6秒开始推杆；(b)表示在落地前约4-9秒开始收油门，先缓慢拉杆，在落地前10秒左右加速拉杆，在落地前1秒左右快速推杆；(c)表示在落地前约5秒开始收油门，先缓慢拉杆，在落地前13秒左右加速拉杆，在落地前8秒左右快速推杆，最后在落地前3秒左右再次拉杆；(d)表示在落地前约7秒开始收油门，先缓慢拉杆，在落地前3秒左右推杆；(e)表示在落地前约9秒开始收油门，先缓慢拉杆，在落地前12秒左右加速拉杆，在落地前4秒左右缓慢推杆；(f)表示在落地前约7秒开始收油门，先缓慢拉杆，在落地前1秒左右快速推杆；(g)表示在落地前约6秒开始收油门，持续缓慢拉杆；(h)表示在落地前约6秒开始收油门，持续缓慢拉杆；(i)表示缓慢收油门，持续缓慢拉杆。综上，第七、八、九类关于杆的操作和训练手册描述较为吻合，在着陆阶段持续拉平，但是第九类的油门操作与第七、八类有较大差异，收油门持续时间较长。

4 飞行操作模式的风险分析

4.1分析方法

在本文研究工作中，以QAR监控标准作为飞行操作模式风险高低的判断依据，采用Kolmogorov-Smirnov检验(K-S检验)，对不同飞行操作模式的监控指标分布情况的差异性进行检验，进一步做出风险评判。K-S检验用于验证一个经验分布是否符合某种理论分布，或用于比较两个经验分布是否有显著性差异。两样本K-S检验对于两样本经验分布函数的位置和形状参数的差异都较为敏感，且不依赖于累积分布函数。

4.2案例分析

针对如前所述特定月份首都机场1210次降落的QAR记录，统计落地G值的分布情况，G值范围在[1.1445，1.8867]。对于样本机型，着陆垂直载荷(G值)的监控标准为1.6，飞行操作模式G值的分布如图4所示。

飞行操作模式G值描述性统计见表2。

为验证以G值超过监控标准样本的占比为依据进行的风险划分是否合理，采用K-S检验对不同飞行操作模式下G值分布的差异性进行检验，结果见表3。

图4 G值分布图

类别类样本数类均值类标准差G值<1.6的样本数G值≥1.6的样本数G值≥1.6样本占比第一类281.270.082800.00%第二类421.320.114112.38%第三类321.300.093200.00%第四类541.280.095400.00%第五类351.290.133412.86%第六类461.280.104512.17%第七类1971.290.0919431.52%第八类2071.290.0920520.97%第九类481.310.094712.08%其他521--51110-总计12101.300.101191191.57%

4.3方法对比

接下来采用神经网络和支持向量机(SVM)来分析QAR监控事件着陆垂直载荷大。本文采用BP神经网络，网络包含输入层、隐藏层和输出层共三层，输入层2个单元、隐藏层20个单元、输出层2个单元，后向传播学习率2.5；支持向量机采用LIBSVM，使用默认的C-SVC分类机，采用RBF核函数，参数设置c=1000，gamma=10。

表3 K-S检验结果

采用的数据同上，使用落地时刻的油门(单位，%RPM)、杆位(单位，DEG)两项参数的QAR记录作为模型的输入，将落地时刻的G值是否不小于1.6作为模型的输出。训练集、测试集数据情况如表4。

表4 数据情况

分类结果如下：

表5 BP神经网络分类的混淆矩阵

表6 LIBSVM分类的混淆矩阵

以G值≥1.6的样本为研究对象(正例)，由表5、6可知，神经网络的分类准确率为98.4%，支持向量机的分类准确率为96.9%，准确率都很高；考虑到存在类不平衡的情况，考察两个分类器的灵敏度(真正例率)和特效性(真负例率)，神经网络和支持向量机的灵敏度均为0，神经网络特效性为100%，支持向量机特效性为98.4%。着陆垂直载荷大这一事件本身发生概率很低，(1210个样本中G值≥1.6的样本仅有19个，占比1.57%)，即使将所有样本均判定为G值<1.6，最终的准确率也很高。在这种类不平衡的情况下，应更关注灵敏度和特效性指标，但两个分类器灵敏度表现都很差，因此两个分类器正确标记稀有类(正例)的能力较差。此外，两种分类器除了最终判断结果外，提供的其他有用信息较少。神经网络用于分类预测时，解释性较弱；支持向量机由于需要空间维度转换处理，同样存在解释性较弱的问题。

本文研究落地前一段时间内关键参数的变化特征，不仅可以根据关键参数特征变化的相似性确定操作模式，而且可以进一步分析操作模式的风险高低，本文从参数曲线的相似性出发研究飞行操作模式，直观、解释性强。

5 结语

本文围绕QAR记录的飞行员对油门、杆位等设备的操作参数数据，从操作参数多元时间序列的特征提取入手，进而分析飞行操作的相似性，挖掘蕴含的飞行操作模式。针对不同的飞行操作模式，研究相关QAR监控指标值的分布情况，量化不同飞行操作模式的风险水平。通过对首都机场特定机型在特定季节的着陆阶段QAR记录的分析，对模型的有效性进行了验证。

在本文的研究工作中，对于多元时序数据进行特征提取时主要考虑了其中的趋势信息，后续将结合实际问题，研究时序数据所含其他信息的提取方法对飞行操作进行聚类分析的终止条件也有待进一步研究。此外，飞行过程是对飞机高度、速度、姿态等的复杂控制过程，并受到地理条件、气象条件等不确定性因素的影响，对飞行操作模式及其风险的分析离不开对专家经验深入运用。

本文的研究工作可拓展至对其他交通工具操作模式的研究，比如汽车驾驶员根据传感器记录的驾驶员对方向盘、刹车和油门等设备的操作数据，采用本文方法，可对驾驶员的操作模式进行挖掘，并对风险率水平进行评估，有助于提高安全驾驶水平。

[1] 曹海鹏, 舒平, 黄圣国. 基于神经网络的民用飞机重着陆诊断技术研究[J]. 计算机测量与控制, 2008, 16(7): 906-908.

[2] 祁明亮, 邵雪焱, 池宏. QAR超限事件飞行操作风险诊断方法[J]. 北京航空航天大学学报, 2011, 37(10): 1207-1210.

[3] Shao Xueyan, Qi Mingliang, Gao Mingang. A risk analysis model of flight operations based on region partition. Kybernetes, 2012, 41(10): 1497-1508.

[4] 曹海鹏. B737重着陆诊断技术研究[D]. 南京:南京航空航天大学, 2008.

[5] 张荣家, 钟力生. A321机型大过载着陆的QAR数据分析和探讨[J]. 中国民航飞行学院学报, 2012, 23(1): 49-53.

[6] 孙瑞山, 韩文律. 基于差异检验的飞行超限事件参数特征分析[J]. 中国安全生产科学技术, 2011, 7(2): 22-27.

[7] 孙瑞山, 肖亚兵. 基于QAR记录数据的民航飞行员操作特征指标结构研究[J]. 中国安全生产科学技术, 2012, 8(11): 49-54.

[8] Li Lishuai. Anomaly detection in airline routine operations using flight data recorder data [D]. Cambridge: MIT, 2013.

[9] 闫相斌, 李一军, 邹鹏, 等. 动静态属性数据相结合的客户分类方法研究[J]. 中国管理科学, 2005, 13(2): 95-100.

[10] Cai Qinglin, Chen Ling, Sun Jianling. Piecewise statistic approximation based similarity measure for time series[J] Knowledge-Based Systems,2015,85:181-195.

[11] 徐梅, 黄超. 基于符号时间序列方法的金融收益分析与预测[J]. 中国管理科学, 2011, 19(5): 1-9.

[12] 李正欣，张凤鸣，李克武. 多元时间序列模式匹配方法研究[J]. 控制与决策，2011，26(4):565-570.

[13] 吴虎胜, 张凤鸣, 张超, 等. 多元时间序列的相似性匹配[J]. 应用科学学报, 2013, 31(6): 643-649.

[14] Keogh E J, Pazzani M J. Scaling up dynamic time warping to massive datasets[C]//Proceedings of the Third European Conference on Principles and Practice of Knowledge Discovery in Databases，Prague,Czech Republic,September 15-18,1999.

[15] 毛红保，张凤鸣，冯卉. 基于奇异值分解的飞行动作评价方法研究[J].计算机工程与应用，2008，44(32):240-245.

[16] 管河山，姜青山，王声瑞. 基于点分布特征的多元时间序列模式匹配方法[J]. 软件学报，2009, 20(1): 67-79.

[17] ten Holt G A, Reinders M J T, Hendriksa E A. Multi-dimensional dynamic time warping for gesture recognition[C]//Proceedings of the 13th Annual Conference of the Advanced School for Computing and Imaging，2007.

[18] Bankó Z, Abonyi J. Mixed dissimilarity measure for piecewise linear approximation based time series applications[J].Expert Systems with Applications,2015,42(21): 7664-7675.

[19] 曹惠玲, 张浩, 王立鑫. QAR译码技术研究与自动译码系统的建立[J]. 航空维修与工程, 2016,(5): 74-77.

[20] 刘陶, 李晓白, 郞荣玲, 等. 飞行参数的分段线性表示及其应用[J]. 微计算机信息, 2010, 26(5-1): 199-201.

[21] Boeing Commercial Airplanes. Statistical summary of commercial jet airplane accidents worldwide operations 1959-2015[R]. Seattle, WA: Boeing Commercial Airplanes, 2015.

PatternRecognitionandRiskAnalysisforFlightOperations

ZHENGLei1,2,CHIHong1,SHAOXue-yan1

(1.Institutes of Science and Development, Chinese Academy of Sciences, Beijing 100190, China；2.University of Chinese Academy of Sciences, Beijing 100049, China)

The analysis of QAR data is important to continuously improving the quality of flight operations. During the flight, the pilot controls the equipment, such as the rod, the plate, the rudder according to the dynamic changes of environmental conditions and the state of the aircraft. It is a process of constant adjustment and coordination, which increases the difficulty of data analysis. So whether pilots have similar operation patterns and what effects these operations hase on the QAR monitoring indexare of great interest to us.

1003-207(2017)10-0109-10

10.16381/j.cnki.issn1003-207x.2017.10.012

V328.1;C931

A

2016-06-30；

2017-03-27

邵雪焱(1978-)，女(汉族)，山东巨野人，中国科学院科技战略咨询研究院助理研究员，研究方向：风险管理、数据挖掘；E-mail: xyshao@casipm.ac.cn.

In this paper, by studying the feature extraction method of the multivariate time series data of flight parameters, the definition of similarities of flight operations is analyzed.The piecewise linear fitting based multivariate Dynamic Time Warping distance is employed to depict the similarities of flight operations.The hierarchical clustering analysis is used to recognize the similar patterns of flight operations. And then, the descriptive statistics and the Kolmogorov-Smirnov test is adopted to quantify the relations between flight patterns and the QAR monitoring index. The judgement of risk levels is obtained. Finally, the validity of the model is verified by using the actual QAR data recorded during the landing stage of a specific aircraft.Other classifiers like BP Neural Networks and Support Vector Machine are used to compare with the proposed method.It turns out that the raised method provides an effective way to analyze flight operationsand the relationship between flight patterns and the QAR monitoring index.

In the future studies,focus will be put on the better description of multivariate time series and clustering methods for multivariate time series. The proposed approach could also be applied in the analysis of other vehicle driving, for example the monitoring of car driving. The method advocated couldhelp tofind the recurring patterns of drivers and how they affectsafety.

Keywords： QAR; multivariate time series; operation pattern; risk analysis