基于k-SC聚类的飞行操作模式及危险性分析

2021-10-13孙瑞山李重锋

中国安全生产科学技术 2021年9期

孙瑞山，李重锋

(中国民航大学安全科学与工程学院，天津 300300)

0 引言

2020年《中国民航航空安全报告》表明：2010—2019年的10 a间，我国民航由机组原因导致的飞行事故占到67.90%[1]，其中飞行员操作偏差是导致不安全事件发生的重要原因。快速存取记录器(QAR)可全面记录飞行员手动操作过程中产生的各类参数，包括飞行员操作参数，飞机性能参数与运行环境参数。飞行操作模式是指飞行员具有“代表性”的操作手法，基于QAR数据对其进行挖掘并建立其与不安全事件之间的映射关系，对于有针对性地提升飞行员操作能力、改善航司飞行训练水平、保障飞机的运行安全等具有重要意义。

当前已有学者基于QAR数据对飞行员操作特征展开研究。首先是利用QAR截面数据对飞行员操作特性进行研究。例如，孙瑞山等[2]结合QAR数据及相关监控标准，确定了描述飞行员操作的完整性、平稳性、准确性和及时性的操作特征向量指标；董传亭[3]利用功率谱密度方法获得飞行操作的平稳性特征指标，进一步基于QAR数据构建出1套飞行操作评价体系；祁明亮等[4]以着陆阶段的1类QAR超限事件为研究对象建立数学规划模型，并针对该类事件的操作参数构成的“高风险子空间”进行分析。以上利用截面数据，即QAR数据的状态值进行的研究忽略了数据的时序规律，特别是对于连续的飞行员操作而言，更需要从时间维度上挖掘相关数据的信息。针对这个问题，学者们将时间序列的相关算法应用于飞行员操作的相关研究中。例如，汪磊等[5]选取无线电高度、驾驶杆位以及油门杆位作为评价飞行着陆操作的关键指标，构建基于曲线相似度的飞行着陆操作评价模型；郑磊等[6]使用基于分段线性表示和动态时间规整(DTW)的系统聚类方法挖掘QAR数据中蕴含的飞行操作模式，并在此基础上利用卷积神经网络模型构建精度较高的重着陆预警模型[7]，结果表明在已知飞行操作模式的前提下构建重着陆预警模型的召回率更高。

挖掘飞行数据中飞行操作特征的研究相对较少，如Wang等[8-9]基于QAR数据应用统计建模的方法对着陆阶段的事件及操作风险进行分析，指出拉平操作对着陆绩效的重要性；Matthews等[10]基于飞行品质监控数据提出1种新的多元时间序列搜索算法，用于识别包括人为因素在内的各类因素导致的重大运行事件；Kraemer等[11]基于飞行模拟器数据利用统计和聚类分析飞行员在正常、发动机故障和襟翼故障3种情况下的起飞操作并得出差异化结果。总体来说，当前基于飞行数据的研究更多注重飞行品质监控和异常检测，对于飞行数据中飞行员操作特征的挖掘相对较少。

综上所述，本文引入k-SC聚类算法挖掘QAR数据中蕴含的飞行员操作特征，该算法基于时间序列形状相似度的度量完成飞行员操作数据的时间序列的聚类，比传统的基于欧式距离的算法优越，同时比基于动态时间规整的算法更高效[12]。同时，结合不安全事件对k-SC聚类结果所反映的飞行操作模式的危险性进行分析，从而合理量化飞行操作模式的危险性，以期为和飞行员操作特征相关的研究提供理论参考。

1 时间序列相关理论

1.1 时间序列相似度不变性

分析时间序列的相似度是时间序列聚类的基础，一般采用能识别其多种变化的距离度量方式来衡量[13]。即距离度量方式应满足：

1)尺度变化：时间序列在尺度(缩放)中存在差异。如序列X变换为Y=aX+b，其中a和b为常数，X与Y的相似度不变；

2)位移变化：2个时间序列相位具有一定偏差，如Y(t)=X(t-t0)，X与Y的相似度不变；

3)噪声变化：当2个时间序列具有相似的形态，但受到不同程度噪声干扰时，2者相似度不变。

1.2 常用时间序列距离度量方法

在进行时间序列数据聚类时，常采用以下2种方法来进行相似度计算。

1)欧氏距离：对于同为n维的时间序列X=[x1,…,xn]和Y=[y1,…,yn]，2者的欧氏距离ED(X,Y)如式(1)所示：

(1)

式中：X与Y为时间序列；n为X的维度；ED(X,Y)为X与Y的欧式距离；xi和yi分别为X与Y中某时刻的值。欧氏距离在衡量时间序列相似性时对噪声和异常点较敏感[14]，因此该方法虽然简单高效却很难精确度量时间序列之间的相似性。

2)DTW距离：处理时间序列分类时，序列长度可能不同，同时序列之间可能存在局部的扭曲。DTW允许时间序列弯曲时间轴，灵活地对时间序列进行匹配，实现局部范围内匹配最优。对于n维X=[x1,…,xn]和m维Y=[y1,…,ym]，可利用动态规整方法寻找到规整成本最小的路径，获得最小动态时间弯曲距离，如式(2)所示：

DTW(X,Y)=τ(n,m)

(2)

式中：DTW(X,Y)为X与Y的DTW距离；τ(n,m)为累计距离；m为Y的维度。

k-SC聚类算法基于形状相似性度量处理时间序列在尺度和位移上的扭曲，满足尺度与位移不变性且计算效率高；对比基于ED度量时间序列相似性的k-means聚类算法，k-SC算法具有更高的分类准确性[12]。因此，k-SC聚类算法适用于挖掘海量QAR数据中蕴含的飞行操作模式。

2 飞行操作模式挖掘及危险性分析方法

2.1 问题定义

飞行过程中，飞行员主要通过对杆、油门、舵等的操纵实现对飞机高度、速度、姿态等的控制。QAR可记录飞行操作数据(如杆位、油门位参数等)随时间变化的情况，即操作参数i对应1个长度为n的时间序列Xi=[xi1,…,xin]。通过分析飞行员的QAR时间序列数据，可以发现飞行员控制驾驶杆、油门杆、升降舵等的时机和变化情况，进而与规章、手册的要求或建议进行对比并发现问题；对于飞行员操作中存在的共性特征，即本文定义的飞行操作模式，如“抬头慢”、“拉平早”等还需要进一步对重点安全监控参数，如飞机的着陆垂直载荷、15.24 m至接地平飘距离等进行分析并量化飞行操作模式的危险性，进而有针对性地排除飞行员的操作隐患。

综上所述，本文主要研究2方面问题：

1)飞行操作参数的时间序列聚类，即从某一参数的时序数据中挖掘相似的飞行操作模式。

2)分析飞行操作模式与QAR监控事件之间的关联关系，量化不同飞行操作模式的危险性。

2.2 飞行操作模式分析

针对问题的定义，飞行操作模式的挖掘可转化为对飞行操作参数序列的相似性研究。k-SC聚类基于某飞行操作参数序列的相似性分析完成聚类。

具体地，对于某操作参数序列集合C的任意2个序列X=[x1,…,xn],Y=[y1,…,yn]，k-SC聚类算法的序列相似度度量方法如式(3)所示：

(3)

该方法找到最优平移q和比例系数α来匹配2个时间序列的形状。具体为，首先找到q′使X,Y在同一时间达到峰值，然后在q′周围寻找最优q。固定q后，式(3)是1个关于α的凸问题，通过设置梯度为0计算出α最优值如式(4)所示：

(4)

式中：α*为α的最优值。

在飞行操作参数序列的相似性度量基础上进行k-SC聚类，定义方程F如式(5)所示：

(5)

式中：F为度量聚类效果评价值，F越小，k-SC聚类效果越好；K为聚类数；Xi为某个类中的时间序列；Ck为第k类的序列集合；μk为第k类的聚类中心。

问题转化为寻找最优聚类中心使得F最小，如式(6)所示：

(6)

由式(3)与式(6)得式(7)：

(7)

式中：αi和qi为对Xi的比例系数和平移量。

考虑(不失一般性)Xi已经被qi移位并将α*代入式(7)，结果如式(8)所示：

(8)

式(8)化简后如式(9)所示：

(9)

此外，采取经典的肘部法则来确定最佳聚类数[16]。不同K值对应不同的F，随着聚类数K值的增加，每类样本数的减少，样本距其聚类中心的距离减小，平均畸变程度降低。K值增大过程中，F值下降幅度最大的位置对应的K′值就是肘部，K′即为最佳聚类数。根据k-SC聚类的特点对肘部法则进行改进，采用改进的误差平方和(SSE)确定，结合式(3)得出改进的SSE计算方法如式(10)所示：

(10)

式中：SSE′为改进的误差平方和。

设定初始聚类数Kmin和最大聚类数Kmax后，具体的飞行操作模式聚类流程如图1所示。

图1 飞行操作模式挖掘流程Fig.1 Flow chart of flight operation patterns mining

2.3 飞行操作模式危险性分析

以QAR不安全事件的监控标准作为飞行操作模式危险性高低的判断依据，采用Kruskal-Wallis检验(K-W检验)，对飞行操作模式聚类结果监控参数分布的差异性进行检验[17]，进一步做出危险性评判。K-W检验是1种推广的多独立样本平均秩检验，属于非参数检验。K-W检验的原假设H0为：样本来自的多个独立总体的分布无显著差异。

基于K-W检验的飞行操作模式危险性分析步骤为：

1)确定机队飞行操作模式危险性分析的不安全事件及对应的监控参数和监控标准。

2)依照k-SC飞行操作模式聚类结果对数据进行分组并对各组不安全事件发生的占比进行统计，选取平均水平作为划分危险性高低的依据。

3)针对所选取的不安全事件监控参数对各组聚类结果进行K-W检验。

4)基于K-W检验结果和各组不安全事件占比分析飞行操作模式的危险性。

3 案例分析

3.1 数据处理

根据波音公司2009至2018年的统计数据[18]，着陆阶段虽然在时间上只占整个航程的1%左右，但发生的重大事故占全部航段重大事故的24%左右；波音公司的机组训练手册表明，在着陆阶段飞行员需要操纵驾驶杆完成正常的运动轨迹、速度和正确的配平操作，并在主起落架高于跑道约6 m时操纵驾驶杆逐渐增加俯仰姿态约2°～3°开始拉平，以减缓飞机的下降率；同时，相关学者基于QAR数据应用统计建模的方法对着陆阶段的事件及操作风险进行分析，指出飞行员对驾驶杆的操作，特别是对拉平环节的操作影响着飞机的着陆绩效[8-9]。综上所述，本文采集某机队飞机着陆接地前20 s的飞行参数数据，并以驾驶杆操作的时间序列为例进行研究。

由于机场条件、气象条件等外部因素均会影响到飞行操作，不同机场海拔、跑道条件等环境因素下的飞行操作要求不尽相同，因此，本文依据以下条件对数据进行筛选：落地机场固定；特定机型；着陆重量小于60 000 kg；飞机着陆形态为襟翼30°；机场温度30～35 ℃；顶风风量小于10 m/s，顺风分量小于5 m/s的样本记录进行分析。最终，共收集422次QAR数据记录并利用Python进行批量化预处理，包括采用极值归一化将序列数据压缩在区间[0,1]之间。

3.2 驾驶杆操作模式分析

为说明实验方法的可行性，选取Kmin为2，Kmax为10。利用式(10)计算不同聚类数K对应的SSE′值，其对应关系如图2所示。

图2 K与SSE′之间的关系Fig.2 Relationship between K and SSE′

可以看出SSE′在K=5处变化最大，选取最佳聚类数为5。在此基础上由式(10)得出最终的着陆前20 s驾驶杆位时间序列的各个聚类中心如图3所示。

图3中，横坐标起始时间为飞机着陆前20 s的第1 s。

图3 k-SC聚类结果Fig.3 Results of k-SC clustering

图3所示的5类操作模式中，第1类表示飞行员在第10 s左右开始拉杆，在出现2次幅度较大的推拉杆动作后小幅拉杆接地；第2类表示飞行员初始操作较不稳定，在第11 s左右拉杆，在第18 s左右出现持续1 s左右的小幅推杆后继续拉杆并接地；第3类表示飞行员在第11 s左右开始拉杆，在第19 s左右推杆接地；第4类表示飞行员在第10 s左右拉杆，做出3次小幅推拉杆动作后推杆接地；第5类表示飞行员在第11 s左右大幅拉杆，做出2次幅度较大的推拉杆动作后拉杆接地。

可以看出，5类操作模式对应驾驶杆位时间序列的形状存在一定差异，主要体现在初始拉杆时机、拉杆幅度及拉平稳定性等的控制上。

3.3 驾驶杆操作模式危险性分析

为进一步分析不同操作模式的危险性，选取长着陆事件为分析目标。长着陆是1种着陆不安全事件，会降低跑道的使用效率，增加飞机冲出跑道的概率[19]。长着陆对应飞行品质监控项目为15.24 m至接地距离远，监控参数为15.24 m至接地段的地速积分距离D，具体的监控标准参考波音飞行品质监控轻度超限取值为750 m。按照k-SC聚类结果将D分为5组进行K-W检验，并设检验显著性水平α为0.01，K-W检验结果如表1所示。

表1 K-W检验结果Table 1 Results of K-W test

K-W检验结果表明，不同驾驶杆操作模式的地速积分距离D的分布不同，第1类操作模式的秩平均值最高，第3类操作模式的秩平均值最低。将各组对应长着陆事件数的占比与平均占比水平(55.68%)对比得出高危险性操作模式为第1类和第4类，第4类驾驶杆操作模式对应发生长着陆的危险性最高。

总体来看，特定条件下该机队飞行员的5类驾驶杆操作模式起始拉平时间近似相同，其中长着陆危险性较高的驾驶杆操作模式具有推拉杆次数较多的特点，即拉平至接地段对驾驶杆操作平稳性较差。对比危险性较低的操作模式，可以看出平稳连续地拉杆可以降低飞机的平飘着陆距离，从而减小飞机发生长着陆的危险性。因此，航空公司可从以上角度对机队现存的长着陆危险性较高的飞行操作模式进行改进训练，以持续提升机队安全运行水平。