CN113962283A - 一种基于局部自适应动态时间规整的航空器轨迹聚类方法 - Google Patents

一种基于局部自适应动态时间规整的航空器轨迹聚类方法 Download PDF

Info

Publication number
CN113962283A
CN113962283A CN202111019773.0A CN202111019773A CN113962283A CN 113962283 A CN113962283 A CN 113962283A CN 202111019773 A CN202111019773 A CN 202111019773A CN 113962283 A CN113962283 A CN 113962283A
Authority
CN
China
Prior art keywords
track
distance
tracks
path
point
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202111019773.0A
Other languages
English (en)
Inventor
曾维理
褚晓
包杰
徐正凤
曹普文
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN202111019773.0A priority Critical patent/CN113962283A/zh
Publication of CN113962283A publication Critical patent/CN113962283A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/22Matching criteria, e.g. proximity measures
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/23Clustering techniques
    • G06F18/232Non-hierarchical techniques
    • G06F18/2321Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06QINFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00Administration; Management
    • G06Q10/04Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G06Q10/047Optimisation of routes or paths, e.g. travelling salesman problem

Landscapes

  • Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Human Resources & Organizations (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Economics (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Strategic Management (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Development Economics (AREA)
  • Game Theory and Decision Science (AREA)
  • Probability & Statistics with Applications (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开一种基于局部自适应动态时间规整的航空器轨迹聚类方法,首先,获取进离港航迹数据,并对航迹数据进行质量分析;并对航迹数据进行预处理;其次,对跑道进行区分,为每条轨迹分配跑道,并将在同一个跑道起降的轨迹再按照走廊口继续分类;然后构建局部自适应动态时间规整模型,找出最优路径以及该路径的总体距离;最后,基于划分方法对航空器轨迹进行聚类。本发明借鉴直观视觉区分轨迹的思想,通过人为施加较大权重的方式充分放大轨迹中相似度较低的部分,使得相似度较低部分的相似度更低,且不考虑相似部分对整个轨迹相似度的影响,使得聚类结果更加精确、具有解释性及更加具有针对性。

Description

一种基于局部自适应动态时间规整的航空器轨迹聚类方法
技术领域
本发明属于民航空管自动化与智能化技术领域,具体涉及一种基于局部自适应动态时间规整的航空器轨迹聚类方法。
背景技术
近年来,随着广播式自动相关监视***和二次雷达等技术的逐步成熟,海量航迹数据的获取变得愈发便捷。分析大量的航迹数据发现,许多航班在终端区空域由于航班量大、天气、军方管制或其他原因等并没有按照标准进离场程序飞行,这说明标准终端区进场和离场程序有时并不足以定义具有代表性的路线结构。因而从大量的历史航迹数据中挖掘出有代表性的实际飞行路线成为当下空中交通研究的热点。目前实现该目标最为主流且有效的方式是对轨迹数据进行聚类分析。轨迹聚类是将轨迹的集合划分为若干个具有相似特征子集的过程。其中,聚类中心能够在一定程度上代表大多数轨迹的盛行飞行模式,可以为评估标准飞行程序制定的合理性提供参考。通过观察轨迹聚类发现的代表轨迹和异常轨迹,能够窥探出当前空域结构的不足,对于空域扇区的科学划分、空域资源的高效利用和空域运行的安全性都有着至关重要的作用。另外,航迹聚类还是航迹预测的前提和基础,同时也是航迹生成和空域运行仿真模拟的一个途径。
其中如何定义对象之间的相似性对于聚类结果具有显著影响,因为它决定了聚类搜索的方向。相似性一般用对象之间的距离来定义。由于航迹长度不一致,目前对于轨迹长度的处理一般有两种方法,一种是对轨迹进行重采样、平滑和重建处理,使之具有相同的长度,进而采用一般的距离公式,如欧式距离等计算轨迹间的距离作为相似性;另一种则是对轨迹长度不做任何处理,直接采用能够度量不同长度轨迹间相似性的方法,如动态时间规整、弗雷歇距离、豪斯多夫距离等。
但在聚类过程中发现,有些轨迹由于大部分阶段距离比较近,但存在部分阶段的距离比较远,这种明显属于不同类别的轨迹在聚类时往往将其归为一类。通过调整参数将类别数变大,虽然有时能够解决该问题,但可能又会导致本来属于一类的轨迹变成两类。所以,从距离度量方式上入手,针对能够度量不同长度轨迹间相似性的动态时间规整方法,对其进行一定的改进来解决上述问题。
发明内容
发明目的:针对动态时间规整方法存在的上述问题,本发明提出一种基于局部自适应动态时间规整的航空器轨迹聚类方法,使得聚类结果更加精确、具有解释性及更加具有针对性。
技术方案:本发明所述的一种基于局部自适应动态时间规整的航空器轨迹聚类方法,包括以下步骤:
(1)获取进离港航迹数据,并对航迹数据进行质量分析;
(2)对步骤(1)获取的航迹数据进行预处理;
(3)对跑道进行区分,为每条轨迹分配跑道,并将在同一个跑道起降的轨迹再按照走廊口继续分类;
(4)构建局部自适应动态时间规整模型,找出最优路径以及该路径的总体距离;
(5)基于划分方法对航空器轨迹进行聚类。
进一步地,步骤(1)所述的航迹数据包括飞越、起飞和降落航迹数据。
进一步地,步骤(1)所述的对航迹数据进行质量分析包括:合并一架航班的一天内的所有轨迹;对每天的航班按照年月日加航迹序数的方式添加航迹编号,并对每条轨迹的所有轨迹点按照时间先后顺序排序;分析数据中有无缺失、重复和异常数据;分析有无轨迹点数过少和前后两点间间隔过大的轨迹。
进一步地,所述步骤(2)包括以下步骤:
(21)数据清洗:对于数据中重复的航迹点,保留第一个重复的航迹点,删除其他所有重复的航迹点数据;数据中存在的缺失值采用删除法和插补法;
(22)坐标转换:将二次雷达数据中的经纬度坐标由地理坐标系转换为以机场中心位置为原点的东北天直角坐标***的坐标值xeast、ynorth、zup
xeast=-sinθ0×(X-x0)+cosθ0×(Y-y0) (7)
Figure BDA0003240986090000021
Figure BDA0003240986090000022
其中,X、Y、Z为对应地心地固坐标***的坐标值;(x0,y0,z0)为终端区的中心的地心地固坐标,
Figure BDA0003240986090000031
表示原点纬度对应的弧度,θ0表示原点经度对应的弧度;
(23)特征规范化:对清洗后的航迹信息数据集进行最大-最小规范化,把数据映射到[-1,1]区间内;
(25)衍生属性构造:将所有进场和离场轨迹合并在一起;统计有多少个不同的时刻;统计每一时刻分别有多少个不同的进场航迹编号和离场航迹编号,并将其分别作为该时刻的当前进场航班量和离场航班量。
进一步地,所述步骤(3)实现过程如下:
在每条跑道附近划设一片包含跑道区域和***区域的矩形区域,以跑道中点处垂直于两跑道边缘画线,作为矩形区域的分界线,将该矩形区域按照面积二等分,并为每个区域添加标签以区分,分别对应各自所属跑道;分别统计轨迹中高度在跑道端高度至矩形区域边界处的高度的所有航迹点落在各区域的点数,落在哪个区域的点数最多,则该轨迹就属于该区域起飞或降落;对属于同一个跑道起飞的轨迹,计算它的最后一个轨迹点与各个走廊口的距离;对属于同一个跑道降落的轨迹,计算它的第一个轨迹点与各个走廊口的距离;与哪个走廊口最近,则属于哪个走廊口。
进一步地,所述步骤(4)包括以下步骤:
(41)创建初始距离矩阵:假设轨迹Ti的航迹点序列集合为:
Figure BDA0003240986090000032
轨迹Ti+1的航迹点序列集合为
Figure BDA0003240986090000033
创建一个m×n路径矩阵,其中矩阵的第(cth,dth)个元素是两点
Figure BDA0003240986090000034
Figure BDA0003240986090000035
的欧式距离再赋予一个权重,即
Figure BDA0003240986090000036
||·||2表示两点间的欧式距离;每一航迹点
Figure BDA0003240986090000037
都包含E坐标xeast、N坐标ynorth、U坐标zup、新升降率sjl、当前进场航班量fland、当前离场航班量ftake、航向D等七个属性值;
(42)规整路径服从约束如下:
1)边界点约束:规整路径的开始点和结束点必须是路径矩阵的第一个和最后一个点;
2)连续性约束:路径一次只能前进一步;
3)单调性约束:路径上面的点必须是随着时间单调进行;
(43)动态寻找最优路径:规整路径用L={u1,u2,…,uG}表示,G表示共有几对相似点,u代表规整路径某对相似点的索引,用(c,d)表示;规整路径相似点之间的欧式距离序列用DL={w1,w2,L,wG}表示,wi表示某对相似点之间的欧式距离;
结合连续性和单调性约束,每一个格点的路径只有三个方向:如果路径已经通过了格点(c,d),那么下一个通过的格点为(c+1,d)、(c,d+1)或(c+1,d+1);
使规整代价最小的路径:
Figure BDA0003240986090000041
定义一个累加距离,从(0,0)点开始匹配两个序列,每到一个点,之前所有的点计算的距离都会累加,到达终点(m,n)后,这累加距离就是最后的总的距离;累加距离为当前格点距离,即点
Figure BDA0003240986090000042
Figure BDA0003240986090000043
的欧式距离与可以到达该点的最小的邻近元素的累积距离之和,由此,找出最优路径:
Figure BDA0003240986090000044
(44)局部自适应加权计算相似度:
找出最优路径后,根据路径从原始距离矩阵中找出各个相似点之间的距离,按照从小到大的顺序排列。找出该距离序列的40%分位数pe,将小于pe的距离的权重设为0,即不考虑距离比较近的部分的距离,大于pe的距离的权重设为1.5,即放大轨迹中不相似部分的距离;然后,计算权重为1.5的加权距离的平均值,即为最终两轨迹间的局部自适应的DTW距离。
进一步地,所述步骤(5)包括以下步骤:
(51)从所有轨迹的集合中随机选择K条轨迹作为初始的代表轨迹;
(52)将每条剩余的轨迹分配到最近的代表轨迹,此处采用局部自适应的动态时间规整计算每个剩余轨迹与各个代表轨迹的相似度,每个剩余轨迹被分配到与其相似度最大(局部自适应动态时间规整距离最小)的代表轨迹;
(53)随机地选择一个非代表轨迹orandom
(54)计算用orandom代替代表轨迹oj的总代价S:
Figure BDA0003240986090000051
其中,dist(Ti,Ti+1)表示两条轨迹间的局部自适应动态时间规整距离,
Figure BDA0003240986090000052
是数据集中所有轨迹T与Cj的代表轨迹oj的绝对误差之和,
Figure BDA0003240986090000053
是数据集中所有轨迹T与Crandom的代表轨迹orandom的绝对误差之和;
(55)如果S<0,那么orandom替换oj,形成新的K个代表轨迹的集合;
(56)重复步骤(53)-(55),直到代表轨迹不发生变化;
(57)采用轮廓系数法确定K值;值在-1到1之间,K值越接近1,则轨迹聚类越合理;K值越接近-1,则轨迹应该聚类到其他簇;K值接近0,则样本在两个簇的边界上。
有益效果:与现有技术相比,本发明的有益效果:
1、对于轨迹部分过于相似、轨迹部分间距离过大的轨迹,传统的距离度量可能错误地将其归为一类,调整参数值都不能将其分开;本发明通过将轨迹相似的部分的距离权重设为0,即不考虑相似部分,并将不相似部分的距离放大,即设置较大的权重,从而将其分开,解决了上述问题;
2、对异常值不敏感:由于选用了基于划分的K-中心点聚类方法,它不同于K-均值聚类,前者选择实际对象作为代表轨迹,后者是计算类内所有对象的均值作为代表对,因此后者对异常值较敏感,而前者不会;
3、提出了当前进场航班量和当前离场航班量两个对聚类结果有重要影响的特征,能够将原来不加该特征属于一类的轨迹按照空域当前航班量状态分成多类轨迹,使得聚类结果更加精确和具有解释性;
4、按照跑道和走廊口先把轨迹分成几类,即把数据集首先分开,然后在此基础上再分别单独对每一类轨迹聚类,能够使聚类结果更加具有针对性,防止出现明显的错误,比如由于距离比较近,将本来不属于某跑道(或走廊口)的轨迹归到了该跑道(或走廊口)。
附图说明
图1为本发明的流程图;
图2为跑道分区示意图;
图3为规整路径示意图。
具体实施方式
下面结合附图对本发明做进一步详细说明。
本发明提出一种基于局部自适应动态时间规整的航空器轨迹聚类方法,如图1所示,具体包括以下步骤:
步骤1:获取进离港航迹数据,并对航迹数据进行质量分析。
根据已有某地区的所有航迹数据,包括飞越、起飞和降落航迹,提取属于某特定终端区的起飞和降落航迹数据,即以某机场为中心,80公里为半径的区域内在该机场起降的所有航迹,包括信息记录时间、航班号、航空器所处位置信息(经度、纬度和高度)、航向、升降率、航空器速度等,并划分进离港数据,此后的所有步骤都是针对进场和离场轨迹分开处理的。
将属于一天的所有轨迹合并。考虑到一架航班的起飞或降落可能分布在前一小时的后几分钟和后一小时的前几分钟,航班在终端区的起降时间不等,少则15分钟,多则半小时,所以在进行数据质量探索之前,首先合并一天24小时的航迹数据。
添加航迹编号区分所有轨迹。一天中可能会有一个航班号对应多条轨迹。所以,对这类轨迹,如果存在连续两点间的时间间隔超过15分钟,那么将其航班号人为区分开来,该间隔前是一条轨迹,该间隔后是另一条轨迹。即由于航空器不可能降落后马上起飞或起飞后马上降落,对于一个航班号对应多条轨迹的,肯定存在有连续两点间的时间间隔大于15分钟,另外,天与天之间、月与月之间、年与年之间更会存在相同的航班号,所以为了区分航迹,对每天的航班按照年月日加航迹序数的方式添加航迹编号,并对每条轨迹的所有轨迹点按照时间先后顺序排序。这样,对于一年以内的所有航迹数据,航迹之间都可以区分开来。
分析数据中有无缺失、重复和异常数据。对比航迹数据集中所有航迹点的所有属性,当两个航迹点之间的所有属性均相同时,则这两个航迹点为重复值;检查航迹点的属性值,当存在空值时,则该航迹信息为缺失值;异常值的判断采用聚类法,由于航迹信息数据是多变量数据,针对多变量异常值处理,采用快速聚类法,将数据对象分组成多个簇,从而挖掘出孤立点,判断出异常值数据。
分析有无轨迹点数过少和前后两点间间隔过大的轨迹。统计所有轨迹的点数和每条轨迹所有航迹点间的时间间隔,分析有无轨迹点数过少或者前后两点间时间间隔过大的航迹。
步骤2:对步骤1获取的航迹数据进行预处理。
数据清洗。为了保证后续聚类的质量和效率,对于数据中重复的航迹点,保留第一个重复的航迹点,删除其他所有重复的航迹点数据;数据中存在的缺失值采用删除法和插补法。当数据缺失属性占比大于85%时,删除该航迹点的信息数据,否则根据其他非缺失的属性值,运用回归的方法对缺失值进行插补;对于异常值,统一采取删除的方法;删除轨迹点数小于90和前后连续两点间时间间隔超过90秒的轨迹,因为前者不能体现整个进近过程,后者缺失航迹点过多,通过插值等方式补充会导致误差过大,影响后续聚类质量;
坐标转换。为了简化空间距离计算,更直观地呈现航空器航迹信息变化,需要将二次雷达数据中的经纬度坐标由地理坐标系转换为以机场中心位置为原点的东北天直角坐标***。首先将数据由地理坐标系转换为地心地固直角坐标系。其主要坐标变化过程如下:
1)计算空间球坐标***对应的椭球扁率F、偏心率e和曲率半径r:
Figure BDA0003240986090000071
e=F×(2-F) (2)
Figure BDA0003240986090000072
其中,
Figure BDA0003240986090000081
表示纬度对应的弧度,a表示地球长轴半径,b表示地球短轴半径。然后,计算对应地心地固坐标***的坐标值X、Y、Z:
Figure BDA0003240986090000082
Figure BDA0003240986090000083
Figure BDA0003240986090000084
其中,θ表示经度对应的弧度,h表示高度。
2)将地心地固坐标转化为以机场中心位置为原点的东北天直角坐标***,其Z轴正向与椭球的法线相重合,Y轴正向指向北向,X轴正向指向东向。
以终端区的中心的地心地固坐标(x0,y0,z0)作为原点位置,
Figure BDA0003240986090000085
表示原点纬度对应的弧度,θ0表示原点经度对应的弧度,h0表示原点高度。根据原点坐标计算ENU直角坐标***的坐标值xeast、ynorth、zup
xeast=-sinθ0×(X-x0)+cosθ0×(Y-y0) (7)
Figure BDA0003240986090000086
Figure BDA0003240986090000087
特征规范化。对清洗后的航迹信息数据集进行最大-最小规范化,把数据映射到[-1,1]区间内。
对规范化后的航向、E坐标、N坐标、U坐标属性分别乘2,使各特征拥有适当的影响权重。
衍生属性构造。当前进场航班量和离场航班量:航空器在实际运行过程中可能会受到当前空域内航班量的影响,当航班量较大时,航空器一般会严格按照五边进近或起飞,当航班量更大时管制员一般会对飞机实施雷达引导,那时的飞机轨迹更加具有不确定性;而当航班量较小时,飞机轨迹也会有不确定性,比如,为了节省不必要的时间,管制员可能会指挥飞机直接切到最后进场边。所以构造了当前进场航班量和当前离场航班量两个属性。具体构造步骤如下:
①将所有进场和离场轨迹合并在一起;
②针对所有航迹数据,统计时间属性中有多少个不同的值,即时间属性中共有多少个不同的时间值。
③统计每一时刻分别有多少个不同的进场航迹编号和离场航迹编号,并将其分别作为该时刻的当前进场航班量和离场航班量。
步骤3:对跑道进行区分,为每条轨迹分配跑道,并将在同一个跑道起降的轨迹再按照走廊口继续分类。
跑道匹配。首先,对跑道进行分区。为了判断每条轨迹在哪条跑道起飞或降落,在每条跑道附近划设一片平面矩形区域,该矩形区域不只包含跑道区域,还包含部分其***区域。以跑道中点处垂直于两跑道边缘画线,作为矩形区域的分界线,将该矩形区域按照面积二等分,并为每个区域添加标签以区分,分别对应各自所属跑道。比如,如果某机场有两条跑道,则将其分为4个区域,并用A、B、C、D四个标签来区分,如图2所示。其次,为轨迹分配跑道。查看轨迹数据在跑道端和矩形区域边界处的大致高度,大约分别是250米和1000米,然后分别统计轨迹中高度在250米至1000米的所有航迹点落在各区域的点数,落在哪个区域的点数最多,则该轨迹就属于该区域起飞或降落。不考虑轨迹中高度低于250米的航迹点是因为航空器在跑道速度较慢,雷达扫描周期不变的情况下轨迹在跑道上的点相对密集,且在滑行阶段的点可能不属于它起飞或降落的跑道区域,会对判断轨迹所属跑道造成干扰。
走廊口匹配。为了在聚类时更加具有针对性,需要将在同一个跑道起降的轨迹再按照走廊口继续分类。由于只选取了以机场为中心,80公里为半径区域内的轨迹,此区域可能不会包含部分走廊口,因此,对属于同一个跑道起飞的轨迹,计算它的最后一个轨迹点与各个走廊口的距离;对属于同一个跑道降落的轨迹,计算它的第一个轨迹点与各个走廊口的距离;与哪个走廊口最近,则属于哪个走廊口。
步骤4:构建局部自适应动态时间规整模型,找出最优路径以及该路径的总体距离。
创建初始距离矩阵。假设轨迹Ti的航迹点序列集合为
Figure BDA0003240986090000091
轨迹Ti+1的航迹点序列集合为
Figure BDA0003240986090000092
其中每一航迹点
Figure BDA0003240986090000101
都包含E坐标xeast、N坐标ynorth、U坐标zup、新升降率sjl、当前进场航班量fland、当前离场航班量ftake、航向D等七个属性值。创建一个m×n路径矩阵。其中矩阵的第(cth,dth)个元素是两点
Figure BDA0003240986090000102
Figure BDA0003240986090000103
的欧式距离,即
Figure BDA0003240986090000104
||·||2表示两点间的欧式距离。
动态寻找最优路径。规整路径用L={u1,u2,L,uG}表示,G表示共有几对相似点,u代表规整路径某对相似点的索引,用(c,d)表示,如图2所示。规整路径相似点之间的欧式距离序列用DL={w1,w2,…,wG}表示,wi表示某对相似点之间的欧式距离。
规整路径服从几个约束:
·边界点约束:规整路径的开始点和结束点必须是路径矩阵的第一个和最后一个点。
·连续性约束:路径一次只能前进一步。
·单调性约束:路径上面的点必须是随着时间单调进行的。
结合连续性和单调性约束,每一个格点的路径只有三个方向。例如,如果路径已经通过了格点(c,d),那么下一个通过的格点只可能是下列三种情况之一:(c+1,d),(c,d+1)或者(c+1,d+1)。
满足上述约束条件的路径有指数个,但要求的是使得下面的规整代价最小的路径:
Figure BDA0003240986090000105
定义一个累加距离,从(0,0)点开始匹配两个序列,每到一个点,之前所有的点计算的距离都会累加,到达终点(m,n)后,这个累加距离就是最后的总的距离。累加距离的计算方式如式(12)所示。其中,η(c,d)表示某格点的动态更新距离,
Figure BDA0003240986090000106
表示航迹Ti的航迹点
Figure BDA0003240986090000107
和航迹Ti+1的航迹点
Figure BDA0003240986090000108
间的欧式距离,累加距离为当前格点距离,也就是点
Figure BDA0003240986090000111
Figure BDA0003240986090000112
的欧式距离与可以到达该点的最小的邻近元素的累积距离之和。由此,找出最优路径。
Figure BDA0003240986090000113
局部自适应加权计算相似度:
找出最优路径后,根据路径从原始距离矩阵中找出各个相似点之间的距离,按照从小到大的顺序排列,找出该距离序列的40%分位数pe,将小于pe的距离的权重设为0,即不考虑距离比较近的部分的距离,大于pe的距离的权重设为1.5,即放大轨迹中不相似部分的距离。然后,计算权重为1.5的加权距离的平均值,即为最终两轨迹间的局部自适应的DTW距离。此处分位数和权重的设置要根据每一类轨迹的情况,如果轨迹相似的部分较多,分位数应选取的相对大一些,权重的设置则也要根据要聚类轨迹的情况来定,如果想要聚类的轨迹不相似部分间的距离较小,则权重应相对设置的大一些,以放大轨迹中不相似的部分,从而获得理想的聚类结果。
步骤5:基于划分方法对航空器轨迹进行聚类。
1)从所有轨迹的集合中随机选择K条轨迹作为初始的代表轨迹。
2)将每条剩余的轨迹分配到最近的代表轨迹,此处采用局部自适应的动态时间规整计算每个剩余轨迹与各个代表轨迹的相似度,每个剩余轨迹被分配到与其相似度最大(局部自适应动态时间规整距离最小)的代表轨迹。
3)随机地选择一个非代表轨迹orandom
4)计算用orandom代替代表轨迹oj的总代价S,如式(13)所示,
Figure BDA0003240986090000114
其中,dist(Ti,Ti+1)表示两条轨迹间的局部自适应动态时间规整距离,
Figure BDA0003240986090000115
是数据集中所有轨迹T与Cj的代表轨迹oj的绝对误差之和,
Figure BDA0003240986090000116
是数据集中所有轨迹T与Crandom的代表轨迹orandom的绝对误差之和。
5)如果S<0,那么orandom替换oj,形成新的K个代表轨迹的集合。
6)重复步骤3)至5),直到代表轨迹不发生变化。
7)采用轮廓系数法确定K值。
轮廓系数从簇内紧凑性和簇间分离性的角度出发来度量聚类的质量。它的值在-1到1之间,该值越接近1,说明轨迹聚类越合理;该值越接近-1,说明轨迹更应该聚类到其他簇;该值接近0,说明样本在两个簇的边界上。所有轨迹轮廓系数的均值作为总体的轮廓系数,并用其衡量聚类结果整体的好坏。
以南京禄口机场终端区为例,采用了2019年7月20日至2019年8月11日之间的飞行数据,总共包括10311条轨迹数据,其中降落轨迹5073条,起飞轨迹5238条。对于起飞和降落聚类的所有轨迹,采用轮廓系数法比较采用标准动态时间规整和局部自适应的动态时间规整计算相似度得出的效果,分别为0.6和0.8。改进的动态时间规整方法——局部自适应的动态时间规整能够使得整体的聚类效果比使用标准动态时间规整效果要好。

Claims (7)

1.一种基于局部自适应动态时间规整的航空器轨迹聚类方法,其特征在于,包括以下步骤:
(1)获取进离港航迹数据,并对航迹数据进行质量分析;
(2)对步骤(1)获取的航迹数据进行预处理;
(3)对跑道进行区分,为每条轨迹分配跑道,并将在同一个跑道起降的轨迹再按照走廊口继续分类;
(4)构建局部自适应动态时间规整模型,找出最优路径以及该路径的总体距离;
(5)基于划分方法对航空器轨迹进行聚类。
2.根据权利要求1所述的基于局部自适应动态时间规整的航空器轨迹聚类方法,其特征在于,步骤(1)所述的航迹数据包括飞越、起飞和降落航迹数据。
3.根据权利要求1所述的基于局部自适应动态时间规整的航空器轨迹聚类方法,其特征在于,步骤(1)所述的对航迹数据进行质量分析包括:合并一架航班的一天内的所有轨迹;对每天的航班按照年月日加航迹序数的方式添加航迹编号,并对每条轨迹的所有轨迹点按照时间先后顺序排序;分析数据中有无缺失、重复和异常数据;分析有无轨迹点数过少和前后两点间间隔过大的轨迹。
4.根据权利要求1所述的基于局部自适应动态时间规整的航空器轨迹聚类方法,其特征在于,所述步骤(2)包括以下步骤:
(21)数据清洗:对于数据中重复的航迹点,保留第一个重复的航迹点,删除其他所有重复的航迹点数据;数据中存在的缺失值采用删除法和插补法;
(22)坐标转换:将二次雷达数据中的经纬度坐标由地理坐标系转换为以机场中心位置为原点的东北天直角坐标***的坐标值xeast、ynorth、zup
xeast=-sinθ0×(X-x0)+cosθ0×(Y-y0) (7)
Figure FDA0003240986080000011
Figure FDA0003240986080000012
其中,X、Y、Z为对应地心地固坐标***的坐标值;(x0,y0,z0)为终端区的中心的地心地固坐标,
Figure FDA0003240986080000021
表示原点纬度对应的弧度,θ0表示原点经度对应的弧度;
(23)特征规范化:对清洗后的航迹信息数据集进行最大-最小规范化,把数据映射到[-1,1]区间内;
(24)衍生属性构造:将所有进场和离场轨迹合并在一起;统计有多少个不同的时刻;统计每一时刻分别有多少个不同的进场航迹编号和离场航迹编号,并将其分别作为该时刻的当前进场航班量和离场航班量。
5.根据权利要求1所述的基于局部自适应动态时间规整的航空器轨迹聚类方法,其特征在于,所述步骤(3)实现过程如下:
在每条跑道附近划设一片包含跑道区域和***区域的矩形区域,以跑道中点处垂直于两跑道边缘画线,作为矩形区域的分界线,将该矩形区域按照面积二等分,并为每个区域添加标签以区分,分别对应各自所属跑道;分别统计轨迹中高度在跑道端高度至矩形区域边界处的高度的所有航迹点落在各区域的点数,落在哪个区域的点数最多,则该轨迹就属于该区域起飞或降落;对属于同一个跑道起飞的轨迹,计算它的最后一个轨迹点与各个走廊口的距离;对属于同一个跑道降落的轨迹,计算它的第一个轨迹点与各个走廊口的距离;与哪个走廊口最近,则属于哪个走廊口。
6.根据权利要求1所述的基于局部自适应动态时间规整的航空器轨迹聚类方法,其特征在于,所述步骤(4)包括以下步骤:
(41)创建初始距离矩阵:假设轨迹Ti的航迹点序列集合为:
Figure FDA0003240986080000022
轨迹Ti+1的航迹点序列集合为
Figure FDA0003240986080000023
创建一个m×n路径矩阵,其中矩阵的第(cth,dth)个元素是两点
Figure FDA0003240986080000024
Figure FDA0003240986080000025
的欧式距离再赋予一个权重,即
Figure FDA0003240986080000026
||·||2表示两点间的欧式距离;每一航迹点
Figure FDA0003240986080000027
都包含E坐标xeast、N坐标ynorth、U坐标zup、新升降率sjl、当前进场航班量fland、当前离场航班量ftake、航向D等七个属性值;
(42)规整路径服从约束如下:
1)边界点约束:规整路径的开始点和结束点必须是路径矩阵的第一个和最后一个点;
2)连续性约束:路径一次只能前进一步;
3)单调性约束:路径上面的点必须是随着时间单调进行;
(43)动态寻找最优路径:规整路径用L={u1,u2,…,uG}表示,G表示共有几对相似点,u代表规整路径某对相似点的索引,用(c,d)表示;规整路径相似点之间的欧式距离序列用DL={w1,w2,L,wG}表示,wi表示某对相似点之间的欧式距离;
结合连续性和单调性约束,每一个格点的路径只有三个方向:如果路径已经通过了格点(c,d),那么下一个通过的格点为(c+1,d)、(c,d+1)或(c+1,d+1);
使规整代价最小的路径:
Figure FDA0003240986080000031
定义一个累加距离,从(0,0)点开始匹配两个序列,每到一个点,之前所有的点计算的距离都会累加,到达终点(m,n)后,这累加距离就是最后的总的距离;累加距离为当前格点距离,即点
Figure FDA0003240986080000032
Figure FDA0003240986080000033
的欧式距离与可以到达该点的最小的邻近元素的累积距离之和,由此,找出最优路径:
Figure FDA0003240986080000034
(44)局部自适应加权计算相似度:
找出最优路径后,根据路径从原始距离矩阵中找出各个相似点之间的距离,按照从小到大的顺序排列,找出该距离序列的40%分位数pe,将小于pe的距离的权重设为0,即不考虑距离比较近的部分的距离,大于pe的距离的权重设为1.5,即放大轨迹中不相似部分的距离;然后,计算权重为1.5的加权距离的平均值,即为最终两轨迹间的局部自适应的DTW距离。
7.根据权利要求1所述的基于局部自适应动态时间规整的航空器轨迹聚类方法,其特征在于,所述步骤(5)包括以下步骤:
(51)从所有轨迹的集合中随机选择K条轨迹作为初始的代表轨迹;
(52)将每条剩余的轨迹分配到最近的代表轨迹,此处采用局部自适应的动态时间规整计算每个剩余轨迹与各个代表轨迹的相似度,每个剩余轨迹被分配到与其相似度最大(局部自适应动态时间规整距离最小)的代表轨迹;
(53)随机地选择一个非代表轨迹orandom
(54)计算用orandom代替代表轨迹oj的总代价S:
Figure FDA0003240986080000041
其中,dist(Ti,Ti+1)表示两条轨迹间的局部自适应动态时间规整距离,
Figure FDA0003240986080000042
是数据集中所有轨迹T与Cj的代表轨迹oj的绝对误差之和,
Figure FDA0003240986080000043
是数据集中所有轨迹T与Crandom的代表轨迹orandom的绝对误差之和;
(55)如果S<0,那么orandom替换oj,形成新的K个代表轨迹的集合;
(56)重复步骤(53)-(55),直到代表轨迹不发生变化;
(57)采用轮廓系数法确定K值;值在-1到1之间,K值越接近1,则轨迹聚类越合理;K值越接近-1,则轨迹应该聚类到其他簇;K值接近0,则样本在两个簇的边界上。
CN202111019773.0A 2021-09-01 2021-09-01 一种基于局部自适应动态时间规整的航空器轨迹聚类方法 Pending CN113962283A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111019773.0A CN113962283A (zh) 2021-09-01 2021-09-01 一种基于局部自适应动态时间规整的航空器轨迹聚类方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111019773.0A CN113962283A (zh) 2021-09-01 2021-09-01 一种基于局部自适应动态时间规整的航空器轨迹聚类方法

Publications (1)

Publication Number Publication Date
CN113962283A true CN113962283A (zh) 2022-01-21

Family

ID=79460691

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111019773.0A Pending CN113962283A (zh) 2021-09-01 2021-09-01 一种基于局部自适应动态时间规整的航空器轨迹聚类方法

Country Status (1)

Country Link
CN (1) CN113962283A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115630131A (zh) * 2022-12-19 2023-01-20 北京码牛科技股份有限公司 一种轨迹数据时空切片方法、***及电子设备
CN116826977A (zh) * 2023-08-28 2023-09-29 青岛恒源高新电气有限公司 一种光储直柔微电网智能管理***

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115630131A (zh) * 2022-12-19 2023-01-20 北京码牛科技股份有限公司 一种轨迹数据时空切片方法、***及电子设备
CN115630131B (zh) * 2022-12-19 2023-04-07 北京码牛科技股份有限公司 一种轨迹数据时空切片方法、***及电子设备
CN116826977A (zh) * 2023-08-28 2023-09-29 青岛恒源高新电气有限公司 一种光储直柔微电网智能管理***
CN116826977B (zh) * 2023-08-28 2023-11-21 青岛恒源高新电气有限公司 一种光储直柔微电网智能管理***

Similar Documents

Publication Publication Date Title
CN110930770B (zh) 一种基于管制意图和飞机性能模型的四维航迹预测方法
CN109493644B (zh) 一种基于历史航迹数据挖掘的四维航迹推测方法
US20210383706A1 (en) System and methods for improving aircraft flight planning
CN113962283A (zh) 一种基于局部自适应动态时间规整的航空器轨迹聚类方法
CN110196962B (zh) 一种基于核密度估计的飞机速度异常识别方法
CN112905576B (zh) 一种基于农机作业轨迹确定农田和道路的方法及***
US11094206B2 (en) Vertical flightpath optimization
CN106933977B (zh) 一种基于大数据挖掘分类剔除飞行参数野值的方法
CN109215399B (zh) 一种终端区智能化流控策略生成方法
CN110060513A (zh) 基于历史轨迹数据的空中交通管制员工作负荷评估方法
CN110988880A (zh) 一种基于smr目标轨迹的地理信息提取及目标跟踪方法
CN112164247A (zh) 一种基于船舶轨迹聚类的船舶航线预测方法
CN109977546B (zh) 一种基于无监督学习的四维航迹在线异常检测方法
CN112215416B (zh) 智能规划巡检航线***及方法
CN107424441A (zh) 基于Hotelling′s T2统计量的航空器航迹变点检测与估计方法
CN110570693A (zh) 一种基于可靠性的航班运行时间预测方法
CN115064009A (zh) 一种终端区无人机与有人机冲突风险等级划分方法
Xuhao et al. Trajectory clustering for arrival aircraft via new trajectory representation
CN113284369B (zh) 一种基于ads-b实测航路数据的预测方法
Meijers et al. A data-driven approach to understanding runway occupancy time
CN109145415B (zh) 一种空中交通流时距分布分析方法
CN115862385A (zh) 基于标准飞行程序的机场终端区航空器飞行模式挖掘方法
CN111414000A (zh) 用于生成与飞行器在机场基础设施中的移动相关的操作数据的方法和***
CN110942675A (zh) 一种终端区b类空域划设的方法
CN114118578A (zh) 一种基于空中轨迹和大数据预测航班到达时间的计算方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination