CN115862385A - 基于标准飞行程序的机场终端区航空器飞行模式挖掘方法 - Google Patents

Abstract

本发明公开一种基于标准飞行程序的机场终端区航空器飞行模式挖掘方法，首先，对航空器的进离港数据以及标准飞行程序进行质量分析，并经过预处理获得需要的特征向量；其次，采用基于标准程序引导的方法对航迹进行飞行模式挖掘；在聚类时确定中心轨迹的过程中，使用标准轨迹作为部分中心聚类轨迹，并通过动态时间规整计算相似性度量的方法实现轨迹的距离计算。本发明在原本k‑均值挖掘飞行模式的基础上对中心轨迹的选取方式做出改进，以标准程序作为主要聚类中心对轨迹进行分类，减少迭代次数的同时可以提高聚类精确程度、使得结果更加具有解释性。

Description

基于标准飞行程序的机场终端区航空器飞行模式挖掘方法

技术领域

本发明属于民航空管自动化与智能化技术领域，具体涉及一种基于标准飞行程序的机场终端区航空器飞行模式挖掘方法。

背景技术

随着民航的飞速发展，空中交通流量急剧增加，导致了空中交通情况相较于之前更加复杂。在进离场飞行程序不能满足当前空域状态的同时，如何从大量的航迹数据中挖掘出主流飞行模式也成为了如今需要解决的问题。对于机场终端区来说，分析航迹数据、识别机场终端区的飞行模式对航空公司和空中交通管制来说都有着至关重要的作用。对于航空公司而言，航空器的盛行飞行模式可以直观地反应航班的整体运行状况，从而可以对航班进行调整。而空中交通管制员可以根据飞行模式更便捷地判断航迹异常情况，同时可以利用相关数据进行飞行程序的优化。

目前飞行模式挖掘采用的主流方法是基于聚类的方法。聚类是一种广泛使用的无监督和半监督技术，它根据成对距离或相似函数的定义将相似的数据实例分组成簇。它具有结构简单，算法可靠，计算效率高，适用于较大数据集等优点，因此被广泛用于识别相关的飞行操作事件。通过对航迹的聚类可以从海量航迹数据中提取出代表主流交通模式的中心航迹，从而获取终端区交通流的分布状况。

目前使用聚类方法挖掘飞行模式存在三个需要解决的问题：首先，航迹聚类的簇数对于聚类准确程度来说至关重要。而随着航迹数量的不断增加，空中交通流的分布状况更加复杂，簇数将会更加难以确定。其次，当前大体量的航迹数据会大幅增加聚类的迭代次数，因此需要寻找一种合适的方法提高聚类的效率和精确程度。同时，由于目前获取的航迹数据由航迹点构成，计划航迹点的时间序列不匹配将会影响聚类的精度。

飞行模式挖挖掘是四维航迹预测技术、进离场排序技术、飞行模拟仿真技术的一项基础技术。本发明提出了一种基于标准飞行程序引导的飞行模式挖掘聚类方法。首先对获取的航迹数据进行分析，获得航迹点坐标以及升降率、速度、航向、角度等航迹点数据作为聚类的特征数据。在聚类过程中标准飞行程序为中心轨迹的选取提供方案，从而提高飞行模式挖掘的精度。该方法通过标准飞行程序决定初始飞行模式的聚类中心，并在此基础上采用动态时间规整方法实现航迹点的匹配和分类。在迭代的过程中，每一次迭代结果都与标准飞行程序相关联，从而实现与标准程序相关的盛行飞行模式挖掘分类。本发明解决了航迹点时间序列难以匹配以及聚类簇数的确定问题，提高了飞行模式的挖掘精度。

发明内容

发明目的：为了使用标准飞行程序减少聚类迭代次数、提升聚类结果的精度，本发明提出一种基于标准飞行程序的机场终端区航空器飞行模式挖掘方法；采用标准飞行程序轨迹作为部分初始中心轨迹，并且在每次迭代过程，始终保持标准飞行程序轨迹作为每类飞行模式的中心，从而挖掘出的飞行模式更加符合实际。

技术方案：本发明提供了一种基于标准飞行程序的机场终端区航空器飞行模式挖掘方法，具体包括以下步骤：

步骤1：获取进离港起飞航迹、着陆航迹以及飞越航迹数据及标准飞行程序；

步骤2：对获取的数据进行质量分析，并进行预处理；

步骤3：对航迹点数据进行坐标转换及构造多维航迹特征向量；

步骤4：确定每类飞行模式的中心轨迹，并获得对应的飞行模式特征。

进一步地，步骤1所述起飞航迹、着陆航迹以及飞越航迹数据包括信息记录时间、航班号、航空器所处位置信息、航向、升降率、航空器速度。

进一步地，所述步骤2包括以下步骤：

步骤2.1：航空器轨迹合并：获取的航迹数据按照小时进行归类，每个小时内的轨迹点信息为一组；

步骤2.2：给航迹编号：当同一航班号的相邻航迹点时间间隔超过15分钟时，界定这两个航迹点来自于不同的航班；对每天的航班采用年月日加航迹序数的编号方式进行轨迹的编排区分；

步骤2.3：对于重复值、缺失值和异常值：对于重复值采取删除的方法，保留第一组重复数据并删除其余重复数据；对于缺失值采用插补或删除的方式进行处理；缺失率过小的数据采用回归的方法对缺失数据进行填补，缺失率过大的数据则将该航迹点全部信息进行删除；对于异常值，统一采用删除的方法；对于轨迹点数过少和前后两点间间隔过大的轨迹进行删除。

进一步地，所述步骤3包括以下步骤：

步骤3.1：对航迹点位置信息进行坐标转换：

将二次雷达返回的经度数据和纬度数据，转换成东北天坐标，体现与机场之间的距离关系，反映航空器的实际轨迹位置信息；以机场的导航台为坐标原点，创建直角坐标系，得到的东北天坐标的e、n、u坐标数值；以用户所在坐标原点(x₀,y₀,z₀)为坐标轴原点，设定坐标点(x,y,z)在东北天坐标系中的位置为(e,n,u)；同时设坐标原点的经纬高坐标点为LLA₀＝(lon₀,lat₀,alt₀)，则相应的转换公式为：

步骤3.2：构造多维航迹特征向量：

获取的航迹数据包括经度、纬度、高度、速度、航向等多维特征，通过构造升降率a和角度θ对航迹的形状进行完整描述；

构造升降率a和角度θ两个特征相应计算方法如下：

/>

其中，v_j表示t＝j时刻该点的航空器速度；dist指计算两点间欧几里得距离；Px_j、Py_j分别表示了该点对应的横纵坐标；

根据转换后坐标e、n、u，及升降率a和角度θ，以及已知的信息速度v、航向γ等特征，获得每个轨迹点的特征数据P_j：

P_j＝(e,n,u,a,v,γ,θ) (5)

其中，P_j为轨迹在t＝j时刻对应轨迹点，公式(5)表示轨迹在t＝j时刻对应多维轨迹点的e、n、u坐标，及升降率、速度、航向、角度这七个特征参数；

轨迹TR为1组轨迹点序列，j∈[1,m]为轨迹点编号，m为轨迹点总数，则轨迹点随着时间变化的集合为：

TR＝{P₁,P₂,...,P_j,...,P_m} (6)

假设一共有h条轨迹，TD为航空器轨迹的集合，i∈[1,h]为航迹编号，h为轨迹总数；航空器的轨迹集合可用如下集合的形式表示：

TD＝{TR₁,TR₂,...,TR_i,...,TR_h}。 (7)

进一步地，所述步骤4包括以下步骤：

步骤4.1：从历史轨迹数据中搜索与标准飞行程序距离最近的航迹：

寻找与标准飞行程序最接近的实际轨迹：从所有实际轨迹TD_B＝{TR_b1,TR_b2,...,TR_bQ}中使用欧几里得距离方法计算与全部标准程序TD_A＝{TR_a1,TR_a2,...,TR_aR}轨迹的距离，并从TD_B中找到离每一条标准程序轨迹最近的实际轨迹，用TD_cen表示这些轨迹的集合；

步骤4.2：将实际轨迹与TD_cen进行匹配：

分别计算集合TD_B中其他轨迹与这些TD_cen轨迹的动态时间规整距离：

δ(P_ay,P_bz)＝(P_ay-P_bz)² (8)

其中，P_ay定义为TR_a＝{P_a1,P_a2,...,P_ay,...,P_aY}的中心轨迹序列，P_bz定义为TR_b＝{P_b1,P_b2,...,P_bz,...,P_bZ}的实际航迹序列，d(P_ay,P_bz)即为两点间DTW距离，以航迹点DTW距离之和作为该组航迹的相关度；在TD_cen集合将中每条轨迹附近选择若干条实际轨迹作为其TD_C∈TD_B对应簇的匹配轨迹，完成基于标准飞行程序的轨迹初始聚类；实际轨迹TD_B分为了TD_cen、TD_cen中每条轨迹对应的簇轨迹TD_C、剩余轨迹TD_D三类轨迹，在剩余轨迹中将选择新的聚类中心轨迹；

步骤4.3：随机选择剩余的初始中心轨迹：

从TD_B中的剩余轨迹TD_D中随机选择k-r条轨迹作为其余的初始中心轨迹，记为TD’_cen，为TD_D中的其他轨迹进行匹配；集合TD_cen和集合TD’_cen组成了所有的初始中心轨迹；其余的轨迹使用DTW匹配原理分配到TD’_cen集合中轨迹对应的簇中；

步骤4.4：匹配度分析：根据k条中心轨迹，计算实际轨迹与每一条中心轨迹的DTW距离d(P_ay,P_bz)；将d(P_ay,P_bz)计算的子式结果存储在累积距离矩阵D(Y,Z)中，其中D(Y,Z)为Y行Z列的矩阵，便于保存d(P_ay,P_bz)的全部结果，d(Y,Z)作为最优全局距离，表示了两条轨迹之间的相似度；

计算第l条中心轨迹序列TR_l与给定的实际航迹序列的相似度距离，记为d_l(Y,Z)；引入w作匹配度分析：

w＝min(d_l(X,Y)) (10)

每一条实际轨迹TR_s都对应了所匹配的中心轨迹匹配度w_s；采用最小-最大规范化的归一化方法将数据线性映射到[0,1]区间里：

其中，min(w_s)和max(w_s)分别对应所有实际轨迹得到的w值中的最小和最大值，最终得到的w’_s表示了第s条实际轨迹TR_s的归一化匹配度；

步骤4.5：更新的聚类中心：对于中心轨迹为TD_cen所在的类，其中心轨迹的每一个点坐标均不作变化，中心轨迹为TD’_cen所在的类计算每一类坐标的平均值，并选取最接近的轨迹作为新的聚类中心，更新原来的TD’_cen轨迹对应坐标位置，至此完成了第一次迭代；

步骤4.6：重复步骤4.4和步骤4.5，直到聚类中心不再发生变化或者聚类次数达到要求，即聚类结果达到稳定状态后输出结果。

有益效果：与现有技术相比，本发明的有益效果：

1、本发明使用标准飞行程序协助获取初始中心轨迹，使得航迹的分离以空中交通管制规定的标准飞行程序相关联；采用此方法聚类能够得到与实际飞行模式相对应的聚类结果，更加符合航空器的实际飞行规律；

2、使用标准飞行程序轨迹替代随机选取的k-均值初始中心轨迹，在航迹分布复杂且密集的情况下这种改进能够大幅度减少聚类的迭代次数，同时可以在每一次迭代时减少寻找新的聚类中心的步骤，提高计算效率；

3、在采用空间特征进行聚类的同时引入航空器的速度、航向、升降率等特征数据，使得航迹能够根据不同速度、不同升降率等条件得以区分，在最大化利用历史航迹数据信息的同时使得聚类结果更加具有解释性。

附图说明

图1为本发明的流程图；

图2为角度计算示意图；

图3为DTW路径匹配原理图；

图4为进场(左)和离场(右)聚类后航迹分布图；

图5为进场(左)和离场(右)聚类中心分布图。

具体实施方式

下面结合附图对本发明做进一步详细说明。

本发明提供了一种基于标准飞行程序的机场终端区航空器飞行模式挖掘方法，如图1所示，包括以下步骤：

步骤1：获取进离港起飞航迹、着陆航迹以及飞越航迹数据及标准飞行程序。

根据某地区的航迹数据，包含起飞航迹、着陆航机以及飞越航迹，从中该机场特定终端区区域内的航迹信息。实验选取的二次雷达数据需要包含以机场终端区域内的全部起降轨迹，即以机场为中心半径80公里范围内的轨迹数据，包括信息记录时间、航班号、航空器所处位置信息(经度、纬度和高度)、航向、升降率、航空器速度等信息。根据起降高度将所有航线分为进场轨迹和离场轨迹。进离场轨迹将在后续操作中分别进行聚类处理。

步骤2：对获取的数据进行质量分析，并进行预处理。

步骤2.1：航空器轨迹合并。

获取的航迹数据按照小时进行归类，每个小时内的轨迹点信息为一组。为了使得一条航迹的所有航迹点能够集中在同一个组中，将一天的航迹数据进行合并，将24小时的航迹数据归到一个组。

步骤2.2：给航迹编号。

为了保证航班的区分，相同航班号的航班运行间隔会超过一次航班的运行时间。因此，当同一航班号的相邻航迹点时间间隔超过15分钟时，可以界定这两个航迹点来自于不同的航班。为了避免不同的天数也会有相同航班号的情况发生，可以对每天的航班采用年月日加航迹序数的编号方式进行轨迹的编排区分，以此做到一年内每一条轨迹都有其不会重复的编号。

步骤2.3：重复值、缺失值和异常值。

将航迹点的全部属性进行对比分析，对于其中属性完全相同、属性中存在空值的情况，需要通过删除、插补等方式使数据更加完整。对于数据中的异常值，可以通过聚类的方法进行判别，将数据分为多个簇，并找出明显的离群点，进行数据的修正。雷达数据的记录中，可能会产生航迹较短、可分析航迹点过少的航迹。因此需要找出这一类对聚类可能会产生一定干扰的航迹。对于重复值采取删除的方法，保留第一组重复数据并删除其余重复数据；对于缺失值采用插补或删除的方式进行处理。缺失率过小的数据采用回归的方法对缺失数据进行填补，缺失率过大的数据则将该航迹点全部信息进行删除；对于异常值，统一采用删除的方法。

步骤2.4：轨迹点数过少和前后两点间间隔过大的轨迹。

同时，对于轨迹点数过少和前后两点间间隔过大的轨迹，这些都将影响聚类的结果，应当将相应轨迹进行删除

步骤3：对航迹点数据进行坐标转换及构造多维航迹特征向量。

步骤3.1：坐标转换。

二次雷达返回了的坐标数据是经度数据和纬度数据，将其转换成东北天坐标，体现与机场之间的距离关系，才能更好地反映航空器的实际轨迹位置信息。以机场的导航台为坐标原点，创建直角坐标系，得到的东北天坐标的e、n、u坐标数值。以用户所在坐标原点(x₀,y₀,z₀)为坐标轴原点，设定坐标点(x,y,z)在东北天坐标系中的位置为(e,n,u)。同时设坐标原点的经纬高坐标点为LLA₀＝(lon₀,lat₀,alt₀)。则相应的转换公式为：

/>

步骤3.2：构造特征向量。

二次雷达获取的航迹数据包括经度、纬度、高度、速度、航向等多维特征，同时航迹的形状特征仅采用航迹的三位空间数据难以完整表示，可以通过构造升降率a和角度θ对航迹的形状进行更加完整的描述。由此获得的多维航迹特征向量可以更好地描述不同飞行模式，提高聚类结果的解释性。

构造升降率a和角度θ两个特征相应计算方法如下：

其中，v_j表示t＝j时刻该点的航空器速度；dist指计算两点间欧几里得距离；Px_j、Py_j分别表示了该点对应的横纵坐标。其中，计算获得的角度在0°到360°之间，图2中表示了不同航迹点计算的夹角有所不同。并且由于第一个点P₁和最后一个点P_m无法形成夹角，该航迹点的角度由P₂或P_m-1的角度进行替代。

获得了转换后坐标e、n、u，同时根据计算得到的升降率a和角度θ，以及已知的信息速度v、航向γ等特征，最终获得每个轨迹点的特征数据P_j：

P_j＝(e,n,u,a,v,γ,θ) (5)

其中，P_j为轨迹在t＝j时刻对应轨迹点。公式(5)表示第i条轨迹在t＝j时刻对应多维轨迹点的e、n、u坐标，及升降率、速度、航向、角度这7个特征参数。

轨迹TR_i为1组轨迹点序列，j∈[1,m]为轨迹点编号，m为轨迹点总数。则轨迹点随着时间变化的集合为：

TR_i＝{P_i1,P_i2,...,P_ij,...,P_im} (6)

假设一共有h条轨迹，TD为航空器轨迹的集合，i∈[1,h]为航迹编号，h为轨迹总数。航空器的轨迹集合可用如下集合的形式表示：

TD＝{TR₁,TR₂,...,TR_i,...,TR_h} (7)

步骤4：确定每类飞行模式的中心轨迹，并获得对应的飞行模式特征。

步骤4.1：从历史轨迹数据中搜索与标准飞行程序距离最近的航迹。

寻找与标准飞行程序最接近的实际轨迹：从所有实际轨迹TD_B＝{TR_b1,TR_b2,...,TR_bQ}中使用欧几里得距离方法计算与全部标准程序TD_A＝{TR_a1,TR_a2,...,TR_aR}轨迹的距离，并从TD_B中找到离每一条标准程序轨迹最近的实际轨迹，用TD_cen表示这些轨迹的集合。

步骤4.2：将实际轨迹与TD_cen进行匹配。

分别计算集合TR_B中其他轨迹与这些TD_cen轨迹的动态时间规整(DTW)距离。动态时间规整方法被广泛用于衡量不等长时间序列的相似程度。由于时间序列中，存在时间序列长度不相等或者不同序列在时间轴上存在位移的问题，导致了传统的欧几里得距离无法准确地计算两个时间序列之间的相似性。因此可以通过将时间序列进行有效的缩短或延伸，从而获得航路的全局相似度。如图3所示，两条时间序列分别表示飞行程序航路点序列以及实际航路点序列。当两条航路点序列不能做到直接地明显匹配时，通过DTW计算对应相似点的距离之和，从而得到两条轨迹之间的相似距离。即规整路径距离。

DTW的计算方法具体如下所示：

δ(P_ay,P_bz)＝(P_ay-P_bz)² (8)

其中，P_ay定义为TR_a＝{P_a1,P_a2,...,P_ay,...,P_aY}的中心轨迹序列，P_bz定义为TR_b＝{P_b1,P_b2,...,P_bz,...,P_bZ}的实际航迹序列。

为了得到最优的成本之和，即DTW距离，可以使用以下方法进行递归计算：

为了找到每一条TD_cen中的轨迹匹配的簇，在TD_cen集合中将每条轨迹附近选择若干条实际轨迹作为其TD_C∈TD_B对应簇的匹配轨迹，完成基于标准飞行程序的轨迹初始聚类。至此，实际轨迹B分为了TD_cen、TD_cen中每条轨迹对应的簇轨迹TD_C、剩余轨迹TD_D三类轨迹，在剩余轨迹中将选择新的聚类中心轨迹。

步骤4.3：随机选择剩余的初始中心轨迹。

从TD_B中的剩余轨迹TD_D中随机选择k-r条轨迹作为其余的初始中心轨迹，记为TD’_cen，为TR_D中的其他轨迹进行匹配。集合TD_cen和集合TD’_cen组成了所有的初始中心轨迹。其余的轨迹使用DTW匹配原理分配到TD’_cen集合中轨迹对应的簇中。

步骤4.4：匹配度分析。

根据k条中心轨迹，计算实际轨迹与每一条中心轨迹的DTW距离d(P_ay,P_bz)。将d(P_ay,P_bz)计算的子式结果存储在累积距离矩阵D(Y,Z)中，其中D(Y,Z)为Y行Z列的矩阵，便于保存d(P_ay,P_bz)的全部结果，d(Y,Z)作为最优全局距离，表示了两条轨迹之间的相似度。

计算第l条中心轨迹序列TR_l与给定的实际航迹序列的相似度距离，记为d_l(Y,Z)；引入w作匹配度分析：

w＝min(d_l(X,Y)) (10)

每一条实际轨迹TR_s都对应了所匹配的中心轨迹匹配度w_s；采用最小-最大规范化的归一化方法将数据线性映射到[0,1]区间里：

其中，min(w_s)和max(w_s)分别对应所有实际轨迹得到的w值中的最小和最大值，最终得到的w’_s表示了第s条实际轨迹TR_s的归一化匹配度。

步骤4.5：更新的聚类中心。

在对轨迹进行了匹配个归类后，对于中心轨迹为TD_cen所在的类，其中心轨迹的每一个点坐标均不做变化，中心轨迹为TD’_cen所在的类计算每一类坐标的平均值，并选取最接近的轨迹作为新的聚类中心，更新原来的TD’_cen轨迹对应坐标位置。至此完成了第一次迭代。

步骤4.6：迭代至获得最终结果。

重复步骤4.4和步骤4.5，直到聚类中心不再发生变化或者聚类次数达到要求，即聚类结果达到稳定状态后输出结果。

以南京禄口机场为例，采用了南京禄口机场2019年7月20日至8月11日的终端区区域二次雷达数据。选取轨迹14337条，经过数据清洗筛选后得到进场轨迹4891条,离场轨迹5387条。通过误差平方和计算获得进场和离场轨迹所需要聚类的最佳聚类数k值均为22。

对于起飞和降落的轨迹，通过采用轮廓系数对比的方法来比较标准程序引导和随机选取初始聚类中心的聚类结果优劣。从表1得出的结果可以看出，通过标准飞行程序引导改进后的k-均值聚类方法得到的聚类效果优于原始聚类效果。

表1结果对照图

对南京禄口机场终端区进离场航迹的聚类结果进行分析，获得的聚类结果可见图4。进离场航迹聚类后对应的聚类中心如图5所示。

Claims (5)

1.一种基于标准飞行程序的机场终端区航空器飞行模式挖掘方法，其特征在于，包括以下步骤：

步骤1：获取进离港起飞航迹、着陆航迹以及飞越航迹数据及标准飞行程序；

步骤2：对获取的数据进行质量分析，并进行预处理；

步骤3：对航迹点数据进行坐标转换及构造多维航迹特征向量；

步骤4：确定每类飞行模式的中心轨迹，并获得对应的飞行模式特征。

2.根据权利要求1所述的基于标准飞行程序的机场终端区航空器飞行模式挖掘方法，其特征在于，步骤1所述起飞航迹、着陆航迹以及飞越航迹数据包括信息记录时间、航班号、航空器所处位置信息、航向、升降率、航空器速度。

3.根据权利要求1所述的基于标准飞行程序的机场终端区航空器飞行模式挖掘方法，其特征在于，所述步骤2包括以下步骤：

步骤2.1：航空器轨迹合并：获取的航迹数据按照小时进行归类，每个小时内的轨迹点信息为一组；

步骤2.2：给航迹编号：当同一航班号的相邻航迹点时间间隔超过15分钟时，界定这两个航迹点来自于不同的航班；对每天的航班采用年月日加航迹序数的编号方式进行轨迹的编排区分；

步骤2.3：对于重复值、缺失值和异常值：对于重复值采取删除的方法，保留第一组重复数据并删除其余重复数据；对于缺失值采用插补或删除的方式进行处理；缺失率过小的数据采用回归的方法对缺失数据进行填补，缺失率过大的数据则将该航迹点全部信息进行删除；对于异常值，统一采用删除的方法；对于轨迹点数过少和前后两点间间隔过大的轨迹进行删除。

4.根据权利要求1所述的基于标准飞行程序的机场终端区航空器飞行模式挖掘方法，其特征在于，所述步骤3包括以下步骤：

步骤3.1：对航迹点位置信息进行坐标转换：

将二次雷达返回的经度数据和纬度数据，转换成东北天坐标，体现与机场之间的距离关系，反映航空器的实际轨迹位置信息；以机场的导航台为坐标原点，创建直角坐标系，得到的东北天坐标的e、n、u坐标数值；以用户所在坐标原点(x₀,y₀,z₀)为坐标轴原点，设定坐标点(x,y,z)在东北天坐标系中的位置为(e,n,u)；同时设坐标原点的经纬高坐标点为LLA₀＝(lon₀,lat₀,alt₀)，则相应的转换公式为：

步骤3.2：构造多维航迹特征向量：

获取的航迹数据包括经度、纬度、高度、速度、航向等多维特征，通过构造升降率a和角度θ对航迹的形状进行完整描述；

构造升降率a和角度θ两个特征相应计算方法如下：

其中，v_j表示t＝j时刻该点的航空器速度；dist指计算两点间欧几里得距离；Px_j、Py_j分别表示了该点对应的横纵坐标；

根据转换后坐标e、n、u，及升降率a和角度θ，以及已知的信息速度v、航向γ等特征，获得每个轨迹点的特征数据P_j：

P_j＝(e,n,u,a,v,γ,θ) (5)

其中，P_j为轨迹在t＝j时刻对应轨迹点，公式(5)表示轨迹在t＝j时刻对应多维轨迹点的e、n、u坐标，及升降率、速度、航向、角度这七个特征参数；

轨迹TR为1组轨迹点序列，j∈[1,m]为轨迹点编号，m为轨迹点总数，则轨迹点随着时间变化的集合为：

TR＝{P₁,P₂,...,P_j,...,P_m} (6)

假设一共有h条轨迹，TD为航空器轨迹的集合，i∈[1,h]为航迹编号，h为轨迹总数；航空器的轨迹集合可用如下集合的形式表示：

TD＝{TR₁,TR₂,...,TR_i,...,TR_h}。 (7)

5.根据权利要求1所述的基于标准飞行程序的机场终端区航空器飞行模式挖掘方法，其特征在于，所述步骤4包括以下步骤：

步骤4.1：从历史轨迹数据中搜索与标准飞行程序距离最近的航迹：

寻找与标准飞行程序最接近的实际轨迹：从所有实际轨迹TD_B＝{TR_b1,TR_b2,...,TR_bQ}中使用欧几里得距离方法计算与全部标准程序TD_A＝{TR_a1,TR_a2,...,TR_aR}轨迹的距离，并从TD_B中找到离每一条标准程序轨迹最近的实际轨迹，用TD_cen表示这些轨迹的集合；

步骤4.2：将实际轨迹与TD_cen进行匹配：

分别计算集合TD_B中其他轨迹与这些TD_cen轨迹的动态时间规整距离：

δ(P_ay,P_bz)＝(P_ay-P_bz)² (8)

其中，P_ay定义为TR_a＝{P_a1,P_a2,...,P_ay,...,P_aY}的中心轨迹序列，P_bz定义为TR_b＝{P_b1,P_b2,...,P_bz,...,P_bZ}的实际航迹序列，d(P_ay,P_bz)即为两点间DTW距离，以航迹点DTW距离之和作为该组航迹的相关度；在TD_cen集合将中每条轨迹附近选择若干条实际轨迹作为其TD_C∈TD_B对应簇的匹配轨迹，完成基于标准飞行程序的轨迹初始聚类；实际轨迹TD_B分为了TD_cen、TD_cen中每条轨迹对应的簇轨迹TD_C、剩余轨迹TD_D三类轨迹，在剩余轨迹中将选择新的聚类中心轨迹；

步骤4.3：随机选择剩余的初始中心轨迹：

从TD_B中的剩余轨迹TD_D中随机选择k-r条轨迹作为其余的初始中心轨迹，记为TD’_cen，为TD_D中的其他轨迹进行匹配；集合TD_cen和集合TD’_cen组成了所有的初始中心轨迹；其余的轨迹使用DTW匹配原理分配到TD’_cen集合中轨迹对应的簇中；

步骤4.4：匹配度分析：根据k条中心轨迹，计算实际轨迹与每一条中心轨迹的DTW距离d(P_ay,P_bz)；将d(P_ay,P_bz)计算的子式结果存储在累积距离矩阵D(Y,Z)中，其中D(Y,Z)为Y行Z列的矩阵，便于保存d(P_ay,P_bz)的全部结果，d(Y,Z)作为最优全局距离，表示了两条轨迹之间的相似度；

计算第l条中心轨迹序列TR_l与给定的实际航迹序列的相似度距离，记为d_l(Y,Z)；引入w作匹配度分析：

w＝min(d_l(X,Y)) (10)

每一条实际轨迹TR_s都对应了所匹配的中心轨迹匹配度w_s；采用最小-最大规范化的归一化方法将数据线性映射到[0,1]区间里：

其中，min(w_s)和max(w_s)分别对应所有实际轨迹得到的w值中的最小和最大值，最终得到的w’_s表示了第s条实际轨迹TR_s的归一化匹配度；

步骤4.5：更新的聚类中心：对于中心轨迹为TD_cen所在的类，其中心轨迹的每一个点坐标均不作变化，中心轨迹为TD’_cen所在的类计算每一类坐标的平均值，并选取最接近的轨迹作为新的聚类中心，更新原来的TD’_cen轨迹对应坐标位置，至此完成了第一次迭代；

步骤4.6：重复步骤4.4和步骤4.5，直到聚类中心不再发生变化或者聚类次数达到要求，即聚类结果达到稳定状态后输出结果。

Images