CN113284369A - 一种基于ads-b实测航路数据的预测方法 - Google Patents

Abstract

本发明涉及一种基于ADS‑B实测航路数据的预测方法。ADS‑B航路数据通过连接ADS‑B获取实时的扇区航空器相关数据，利用MATLAB处理分析后得到。通过对ADS‑B数据的清洗，聚类，预测得到下一时刻的航路交通流量数据，预测航路拥塞发生的概率，从而降低运行风险，合理规划空域流量，也可以利用当前预测数据进行分析，为将来空域容量的评估提供支撑。

Description

一种基于ADS-B实测航路数据的预测方法

技术领域

本发明涉及一种空中交通流拥塞评估方法，特别涉及一种基于ADS-B实测航路数据的预 测方法，用于空中交叉航路交通流的预测。

背景技术

空中交通拥挤这一概念从上个世纪已经提出，到如今，空中交通工具数量愈来愈大，这 使得有限的空域资源无法再满足现代航空企业的需求。交通拥挤现象愈发严重，很大程度上 影响着空中交通运行的安全与效率。而要缓解空中交通拥挤则必须要先确定空中交通拥挤发 展态势，如此才能采取相应的空中交通管理措施。故而，空中交通拥挤的识别与预测成为了 国际民航学术界的一个重要研究内容。

当前空中交通流量拥塞预测方法主要分为4类：(1)基于飞行计划的确定性流量预测方 法，如基于飞行计划的4D航迹预测方法；(2)基于数理统计的预测方法，包括时间序列法、 回归分析法、卡尔曼滤波和指数平滑法等；(3)非线性预测方法，包括小波分析预测法、混 沌理论预测法和突变理论预测法等；(4)智能预测方法，包括神经网络法、支持向量机等。

然而，在实际运行中，当航班规模巨大时，计算复杂性对系统运算能力构成严峻考验。 空中交通实际运行时出现的随机特性一直是影响4D航迹精度的关键因素，在流量预测中不 可回避；在面对大量的航班数据时基于数理统计的预测方法运算效率较为低下，不宜在实际 运行中使用。

发明内容

鉴于现有技术的状况，而最近几年所研究的非线性预测方法，成果较为显著，本方法基 于非线性预测增加了时间序列的研究，使预测更为精准。

本发明采用ADS-B航路数据，预测航路拥塞发生的概率，从而降低运行风险，合理规划 空域流量；同时也可以利用当前预测数据进行分析，为将来空域容量的评估提供支撑。

本发明采用的技术方案是：一种基于ADS-B实测航路数据的预测方法，利用MATLAB作 为预测平台，包括以下步骤：

步骤一：提取航路ADS-B数据；

采用基于运行环境实测数据的分析方法，数据来源为ADS-B航迹数据，ADS-B为广播式 自动相关监视，它主要实现空对空监视，一般只需机载电子设备：GPS接收机、数据链收 发机及其天线、驾驶舱冲突信息显示器CDTI，不需要任何地面辅助设备即可获得航空器的 各种信息；

数据格式依次是：序号、注册号、ICAO 24位地址码、航班呼号、纬度、经度、高度、速度、航向、垂直速度、接收机IP地址和当前时间，根据ADS-B航迹数据，能够获得该时 间的航空器数据，本方法的研究对象为空中交通流，所以提取有关航空其轨迹的数据即可， 例如纬度，经度和时间，如何确定这些数据属于哪个航空器尤为重要，本方法在提取轨迹点数据时发现在数据接收过程中存在偶尔的数据丢失现象，因此不能以航班呼号作为轨迹唯一 主键，能够结合ICAO 24位地址码一起识别航空器的轨迹，且这一架飞机可能从两地之间往 返多次，所以也不能以ICAO 24位地址码作为轨迹唯一主键；

所以本方法以“ICAO 24位地址码+航班呼号”确定是否为一条轨迹数据，本方法所有 的数据处理使用的是MATLAB，导入轨迹数据后，将一条轨迹数据以double类型数据的存储 在一个cell中，存储信息包括每条轨迹数据的纬度、经度和时间，以下将介绍处理过程；

步骤二：利用Matlab进行航迹处理，包括数据读入、坐标转换、删除异常数据与轨迹 数据库构建四个步骤；

1.数据读入

将航迹数据txt格式的文件导入到MATLAB工作空间；

2.坐标转换

墨卡托投影是一种正轴等角切圆柱投影，投影经线和纬线是直线且相互垂直，墨卡托投 影的角度变形为0，是航空航海制图中常用的投影类型，利用投影的正解公式将经纬度坐标 转换为平面直角坐标，墨卡托投影坐标系以零度子午线或某一指定经线与赤道或某一标准纬 线的交点的投影为原点，零度子午线或某一标准经线的投影方向为x轴，赤道或某一标准纬 线的投影方向为y轴，单位为米，定义地球上某点的经纬度坐标为

投影后该点在墨卡 托投影坐标系内坐标为(X_N,Y_N)，则墨卡托投影的正解公式如方程组(1)所示：

方程组中a为参考椭球体长半轴，b为参考椭球体短半轴，e为第一偏心率，e’为第二 偏心率；λ₀为制定经线数，

为标准纬线数；

3.删除异常数据

根据原始数据利用Matlab绘图时会发现图中出现野点，这样异常的轨迹数据可以分为 零散轨迹、多连轨迹、错位轨迹和缺点轨迹，以下为处理异常数据流程：

1)识别轨迹数据中是否存在相邻两航迹点时间间隔大于100s的情况，如果存在则判定 为多连轨迹数据，将轨迹拆分成多条；

2)第一次计算每条轨迹中航迹点的个数并求平均值，以其十分之一为阈值，识别轨迹 数据中是否存在航迹点数低于的情况，如果存在则判定为零散轨迹数据删除；由于此时很多 航迹只有1-3个航迹点，这些零散数据影响平均航迹点数的取值，导致平均值结果偏小，所 以后续还会做一次平均值计算，清理一遍零散轨迹数据；

3)识别轨迹数据中是否存在锐角的情况，如果存在则判定为错位轨迹数据，取第一个 和第三个航迹点的空间坐标和时间坐标的中点，代替第二个航迹点，反复迭代多次，直到轨 迹中没有锐角为止；

4)识别轨迹数据中是否存在相邻两个航迹点之间的距离大于10km的情况，如果存在则 判定为缺点轨迹数据，当其时间间隔大于1s时，用线性插值方法以1s为步长插入1个等距 等时的航迹点，当其时间间隔为0或1s时，插入20个的航迹点，时间选前一点的时间；

5)第二次计算每条轨迹中航迹点的个数并求平均值，以其十分之一为阈值，识别轨迹 数据中是否存在航迹点数低于的情况，如果存在则判定为零散轨迹数据删除，考虑到第一次 设定的阈值偏小，并且之后又有轨迹数据的拆分和插值，所以需要再次检验结果；

6)再次检查是否存在相邻两航迹点时间间隔大于100s，相邻三个航迹点之间的夹角是 锐角和相邻两个航迹点之间的距离大于10km的情况；

至此，经过错误轨迹数据处理后的轨迹数据，数据准确，路径清晰，可作为空中交通流 研究对象；

步骤三：对航迹进行聚类得到所需航路的交通轨迹；

自动聚类模型分两阶段：第一阶段是计算结构相似度，第二阶段是调用谱聚类算法，采用 Normalized-cut，构建Normalized Laplacian矩阵，计算其特征值及特征向量，并求得谱 系值，自动确定终端区飞行轨迹类别数目，然后对其进行聚类分析；

1.基于飞行轨迹的欧氏距离求解航迹结构相似度

飞行轨迹由轨迹点组成，经过数据预处理的飞行轨迹包含相同数目的飞行轨迹点，且投 影后经纬坐标已转化为直角坐标系中的x值与y值，由此，顺序计算两条飞行轨迹中每对飞 行轨迹点的欧氏距离，然后取平均值、最值，定义两条轨迹的欧氏距离，如公式(2)是计 算两点间的欧氏距离，公式(3)是求解平均欧氏距离计算轨迹相似度的公式

式中：pi,pj代表两条飞行轨迹，均包含飞行轨迹点数N，pin是飞行轨迹pi的第n个点， pjn同理，xin，xjn，yjn，yjn各是飞行轨迹pi，pj第n个点的横、纵坐标，Eu_D(pi,pj) 是pi,pj的欧氏距离，dist(pi,pj)是两个飞行轨迹点的欧氏距离；

运用欧氏距离度量飞行轨迹相似性的优点在于操作步骤简单且易于理解，因而其应用范围 比较广泛；缺点在于受噪声影响较大，如果未预先处理噪声点则欧氏距离难以准确度量飞行 轨迹相似度，在模拟中若存在两条进场飞行轨迹pi与pj，各飞行轨迹点按顺序组成飞行轨 迹点对，计算点对的欧氏距离，然后取平均距离定义为最终轨迹距离；

2.谱聚类算法

谱聚类算法的核心思想与图论相仿，是一种具有高计算性能的现代聚类方法，它依据无向 图的多重分割理论对研究对象进行聚类分析，无向图G(V,E)的顶点V对应聚类问题中的数 据点，用距离函数求得两点间的距离对边E赋权值aij，A的定义是待处理点的亲和度矩阵， 也即G的邻接矩阵，图分割就是指把一个连通图的一些边切断，达到让原本连通的图打散成 一些独立连通的子图的效果，Min-cut是最早的划分准则，但该方法存在一个缺点：对极小 值很敏感，容易产生歪斜切割而导致结果不准确，因而，本方法选用Normalized-cut准则， 引入度限制，平衡划分子类间的大小；

以下是谱聚类的流程：

输入：终端区飞行轨迹集合PS

输出：飞行轨迹簇C{C1,C2，.，Ck}

步骤1：计算结构相似度并构建亲和度矩阵A；

步骤2：构造矩阵Lsym；

步骤3：求解Lsym的前k个最小特征值λ1，λ2，.，λk以及其特征向量v1,v2,.，vk；

步骤4：计算谱系值，求出最大谱系值对应的下标k并作为聚类簇数目构造特征子空间X；

步骤5：最后利用经典聚类算法分析，得到聚类结果C{C1，C2，.，Ck}；

步骤四：根据航迹聚类结果求解时间序列给出所需时间段的流量预测数据；

1.空中交通流时间序列的构建

提取一段时间内通过某航路点的航空器数量来表示流量值，以构建空中交通流的时间序 列，构建空中交通流的时间序列需要选取一个航路点a，取经过航路点a的一条交通流的轨 迹f，它是由很多个航迹点fi构成的，取每个航迹点经过墨卡托投影转换后的横坐标，纵 坐标和时间，然后计算每个航迹点与a的直线距离，距离最近的点fi即为航路点a的到达 航迹点；

通过以上方法得到交通流到达航迹点集合，将所有到达航迹点的时间值提取出来，并进行 升序排列，假设共有n个到达航迹点，得到单位为秒的到达时间数组，在此基础上，按照某 一时间尺度对航班数量进行统计，得到交通流量的时间序列T_f＝[t_f1，t_f2，...，t_fi，...,t_fn’]，序列长 度为n’，假设统计时间间隔Δt单位为分，总观察时间ts的单位为分，t_fj为第j个时间 间隔内通过交汇点的航班数量；

记到达时间数组T中相邻两时间值为t_i、t_i+1,计算t_gi＝t_i+1-t_i，得到航空器时距序列，即 微观交通流时间序列T_g＝[t_g1,t_g2,...，t_gi，...，t_g(n-1)]，序列长度为n-1，由此就得到了航迹的时间 序列；

2.基于最大Lyapunov指数的混沌时间序列预测

基于上述理论，在判定空中交通流量时间序列具备混沌特性的基础上，从Lyapunov指 数的定义出发，在得到时间序列后，重构相空间后经小数据量法求解出最大Lyapunov指数 λ1，最后预测N+1时刻时间序列为公式(4)，其中x_M+1(j)代表相点x_M+1的第j维分量， x_M(j)，x_k(j)，x_K+1(j)的定义类似，

通过求解λ1计算下一时刻的航路预测数据，与实际数据进行对比，识别空中交通拥挤 流量态势，判断航路的拥挤情况，从而实现了对交叉航路的拥挤概率的预测。

本发明的有益效果是：本方法兼顾空中交通流长期的周期性和短期的不确定性等特点， 运用基于时间序列法的航迹预测，对VYK台的航路流量进行大量的预测，经统计后小误差概 率P＝0.9700＞0.95，预测准确率比传统方法方法高出5％，预测结果相对于普通的预测算法 结果拟合性好，表明预测结果能够很好的反映交叉航路拥塞的趋势和规律，可用于潜在冲突 量以及拥塞的预测。

附图说明

图1为本发明航迹数据处理流程图；

图2为本发明错误轨迹数据处理算法流程图。

具体实施方式

下面对本发明进一步说明。

如图1图2所示，一种基于ADS-B实测航路数据的预测方法，利用MATLAB作为预测平台，包括以下步骤：

采用基于运行环境实测数据的分析方法，数据来源为ADS-B航迹数据，ADS-B是Automatic Dependent Surveillance–Broadcast的缩写，译为广播式自动相关监视，它 主要实现空对空监视，一般只需机载电子设备：GPS接收机、数据链收发机及其天线、驾 驶舱冲突信息显示器CDTI，不需要任何地面辅助设备即可获得航空器的各种信息；

数据格式依次是：序号、注册号、ICAO 24位地址码、航班呼号、纬度、经度、高度、速度、航向、垂直速度、接收机IP地址和当前时间，根据ADS-B航迹数据，可以获得该时 间的航空器数据，本方法的研究对象为空中交通流，所以提取有关航空其轨迹的数据即可， 例如纬度，经度和时间，如何确定这些数据属于哪个航空器尤为重要，本方法在提取轨迹点数据时发现在数据接收过程中存在偶尔的数据丢失现象，因此不能以航班呼号作为轨迹唯一 主键，可以结合ICAO 24位地址码一起识别航空器的轨迹，且这一架飞机可能从两地之间往 返多次，所以也不能以ICAO 24位地址码作为轨迹唯一主键；

所以本方法以“ICAO 24位地址码+航班呼号”确定是否为一条轨迹数据，本方法所有 的数据处理使用的是MATLAB，导入轨迹数据后，将一条轨迹数据以double类型数据的存储 在一个cell中，存储信息包括每条轨迹数据的纬度、经度和时间，以下将介绍处理过程；

步骤二：利用Matlab进行航迹处理，包括数据读入、坐标转换、删除异常数据与轨迹 数据库构建四个步骤；

1.数据读入

将航迹数据txt格式的文件导入到MATLAB工作空间；

2.坐标转换

墨卡托投影是一种正轴等角切圆柱投影，投影经线和纬线是直线且相互垂直，墨卡托投 影的角度变形为0，是航空航海制图中常用的投影类型，利用投影的正解公式将经纬度坐标 转换为平面直角坐标，墨卡托投影坐标系以零度子午线(经度＝0°)或某一指定经线(经度 ＝λ₀)与赤道(纬度＝0°)或某一标准纬线

的交点的投影为原点，零度子午线 或某一标准经线的投影方向为x轴，赤道或某一标准纬线的投影方向为y轴，单位为米，定 义地球上某点的经纬度坐标为

投影后该点在墨卡托投影坐标系内坐标为(X_N,Y_N)，则 墨卡托投影的正解公式如方程组(1)所示：

方程组中a为参考椭球体长半轴，b为参考椭球体短半轴，e为第一偏心率，e’为第二 偏心率；

3.删除异常数据

根据原始数据利用Matlab绘图时会发现图中出现野点，这样异常的轨迹数据可以分为 零散轨迹、多连轨迹、错位轨迹和缺点轨迹，以下为处理异常数据流程：

1)识别轨迹数据中是否存在相邻两航迹点时间间隔大于100s的情况，如果存在则判定 为多连轨迹数据，将轨迹拆分成多条；

2)第一次计算每条轨迹中航迹点的个数并求平均值，以其十分之一为阈值，识别轨迹 数据中是否存在航迹点数低于的情况，如果存在则判定为零散轨迹数据删除；(由于此时很 多航迹只有1-3个航迹点，这些零散数据影响平均航迹点数的取值，导致平均值结果偏小， 所以后续还会做一次平均值计算，清理一遍零散轨迹数据)；

3)识别轨迹数据中是否存在锐角的情况，如果存在则判定为错位轨迹数据，取第一个 和第三个航迹点的空间坐标和时间坐标的中点，代替第二个航迹点，反复迭代多次，直到轨 迹中没有锐角为止；

4)识别轨迹数据中是否存在相邻两个航迹点之间的距离大于10km的情况，如果存在则 判定为缺点轨迹数据，当其时间间隔大于1s时，用线性插值方法以1s为步长插入1个航迹 点(等距等时)，当其时间间隔为0或1s时，插入20个航迹点(等距)，时间选前一点的时间；

5)第二次计算每条轨迹中航迹点的个数并求平均值，以其十分之一为阈值，识别轨迹 数据中是否存在航迹点数低于的情况，如果存在则判定为零散轨迹数据删除，(考虑到第一 次设定的阈值偏小，并且之后又有轨迹数据的拆分和插值，所以需要再次检验结果)；

6)再次检查是否存在相邻两航迹点时间间隔大于100s，相邻三个航迹点之间的夹角是 锐角和相邻两个航迹点之间的距离大于10km的情况；

至此，经过错误轨迹数据处理后的轨迹数据，数据准确，路径清晰，可作为空中交通流 研究对象；

步骤三：对航迹进行聚类得到所需航路的交通轨迹；

自动聚类模型分两阶段：第一阶段是计算结构相似度，第二阶段是调用谱聚类算法，采用 Normalized-cut，构建Normalized Laplacian矩阵，计算其特征值及特征向量，并求得谱 系值，自动确定终端区飞行轨迹类别数目，然后对其进行聚类分析。

1.基于飞行轨迹的欧氏距离求解航迹结构相似度

飞行轨迹由轨迹点组成，经过数据预处理的飞行轨迹包含相同数目的飞行轨迹点，且投 影后经纬坐标已转化为直角坐标系中的x值与y值，由此，可以顺序计算两条飞行轨迹中每 对飞行轨迹点的欧氏距离，然后取平均值、最值等，定义两条轨迹的欧氏距离，如公式(2) 是计算两点间的欧氏距离，公式(3)是求解平均欧氏距离计算轨迹相似度的公式

式中：pi,pj代表两条飞行轨迹，均包含飞行轨迹点数N，pin是飞行轨迹pi的第n个点， pjn同理，xin，xjn，yjn，yjn各是飞行轨迹pi，pj第n个点的横、纵坐标，Eu_D(pi,pj) 是pi,pj的欧氏距离，dist(pi,pj)是两个飞行轨迹点的欧氏距离；

运用欧氏距离度量飞行轨迹相似性的优点在于操作步骤简单且易于修改，因而其应用范围 比较广泛；缺点在于受噪声影响较大，如果未预先处理噪声点则欧氏距离难以准确度量飞行 轨迹相似度，在模拟中若存在两条进场飞行轨迹pi与pj，各飞行轨迹点按顺序组成飞行轨 迹点对，计算点对的欧氏距离，然后取平均距离定义为最终轨迹距离；

2.谱聚类算法

谱聚类算法的核心思想与图论(Graph Theory)相仿，是一种具有高计算性能的现代聚类方 法，它依据无向图的多重分割理论对研究对象进行聚类分析，无向图G(V,E)的顶点V对应 聚类问题中的数据点，用距离函数求得两点间的距离对边E赋权值aij，A的定义是待处理 点的亲和度矩阵，也即G的邻接矩阵，图分割(Graph Cut)就是指把一个连通图的一些边切 断，达到让原本连通的图打散成一些独立连通的子图的效果，Min-cut是最早的划分准则， 但该方法存在一个缺点：对极小值很敏感，容易产生歪斜切割而导致结果不准确，因而，本 方法选用Normalized-cut准则，引入度限制，平衡划分子类间的大小；

以下是谱聚类的流程：

输入：终端区飞行轨迹集合PS；

输出：飞行轨迹簇C{C1,C2,.,Ck}；

步骤1：计算结构相似度并构建亲和度矩阵A；

步骤2：构造矩阵Lsym；

步骤3：求解Lsym的前k个最小特征值λ1,λ2,.,λk以及其特征向量v1,v2,.,vk；

步骤4：计算谱系值，求出最大谱系值对应的下标k并作为聚类簇数目构造特征子空间X；

步骤5：最后利用经典聚类算法分析，得到聚类结果C{C1,C2,.,Ck}；

步骤四：根据航迹聚类结果求解时间序列给出所需时间段的流量预测数据；

1.空中交通流时间序列的构建

提取一段时间内通过某航路点的航空器数量来表示流量值，以构建空中交通流的时间序 列，构建空中交通流的时间序列需要选取一个航路点a，取经过航路点a的一条交通流的轨 迹f，它是由很多个航迹点fi构成的，取每个航迹点经过墨卡托投影转换后的横坐标，纵 坐标和时间，然后计算每个航迹点与a的直线距离，距离最近的点fi即为航路点a的到达 航迹点；

通过以上方法得到交通流到达航迹点集合，将所有到达航迹点的时间值提取出来，并进行 升序排列，假设共有n个到达航迹点，得到单位为秒的到达时间数组，在此基础上，按照某 一时间尺度对航班数量进行统计，得到交通流量的时间序列T_f＝[t_f1，t_f2，...,t_fi,...,t_fn’]，序列长 度为n’，假设统计时间间隔Δt单位为分，总观察时间ts的单位为分，t_fj为第_j个时间间 隔内通过交汇点的航班数量；

记到达时间数组T中相邻两时间值为t_i、t_i+1，计算t_gi＝t_i+1-t_i，得到航空器时距序列，即 微观交通流时间序列T_g＝[t_g1，t_g2，...,t_gi,...,t_g(n-1)]，序列长度为n-1，由此就得到了航迹的时间 序列；

2.基于最大Lyapunov指数的混沌时间序列预测

基于上述理论，在判定空中交通流量时间序列具备混沌特性的基础上，从Lyapunov指 数的定义出发，在得到时间序列后，重构相空间后经小数据量法求解出最大Lyapunov指数 λ1，最后预测N+1时刻时间序列为公式(4)(其中x_M+1(j)其中代表相点x_M+1的第j维分量， x_M(j)，x_k(j)，x_K+1(j)的定义类似)，

通过求解λ1计算下一时刻的航路预测数据，与实际数据进行对比，识别空中交通拥挤 流量态势，判断航路的拥挤情况，从而实现了对交叉航路的拥挤概率的预测。

Claims (1)

1.一种基于ADS-B实测航路数据的预测方法，利用MATLAB作为预测平台，其特征在于，包括以下步骤：

步骤一：提取航路ADS-B数据；

采用基于运行环境实测数据的分析方法，数据来源为ADS-B航迹数据，ADS-B为广播式自动相关监视，它主要实现空对空监视，一般只需机载电子设备：GPS接收机、数据链收发机及其天线、驾驶舱冲突信息显示器CDTI，不需要任何地面辅助设备即可获得航空器的各种信息；

数据格式依次是：序号、注册号、ICAO 24位地址码、航班呼号、纬度、经度、高度、速度、航向、垂直速度、接收机IP地址和当前时间，根据ADS-B航迹数据，能够获得该时间的航空器数据，本方法的研究对象为空中交通流，所以提取有关航空其轨迹的数据即可，例如纬度，经度和时间，如何确定这些数据属于哪个航空器尤为重要，本方法在提取轨迹点数据时发现在数据接收过程中存在偶尔的数据丢失现象，因此不能以航班呼号作为轨迹唯一主键，能够结合ICAO 24位地址码一起识别航空器的轨迹，且这一架飞机可能从两地之间往返多次，所以也不能以ICAO 24位地址码作为轨迹唯一主键；

所以本方法以“ICAO 24位地址码+航班呼号”确定是否为一条轨迹数据，本方法所有的数据处理使用的是MATLAB，导入轨迹数据后，将一条轨迹数据以double类型数据的存储在一个cell中，存储信息包括每条轨迹数据的纬度、经度和时间，以下将介绍处理过程；

步骤二：利用Matlab进行航迹处理，包括数据读入、坐标转换、删除异常数据与轨迹数据库构建四个步骤；

1.数据读入

将航迹数据txt格式的文件导入到MATLAB工作空间；

2.坐标转换

墨卡托投影是一种正轴等角切圆柱投影，投影经线和纬线是直线且相互垂直，墨卡托投影的角度变形为0，是航空航海制图中常用的投影类型，利用投影的正解公式将经纬度坐标转换为平面直角坐标，墨卡托投影坐标系以零度子午线或某一指定经线与赤道或某一标准纬线的交点的投影为原点，零度子午线或某一标准经线的投影方向为x轴，赤道或某一标准纬线的投影方向为y轴，单位为米，定义地球上某点的经纬度坐标为

投影后该点在墨卡托投影坐标系内坐标为(X_N,Y_N)，则墨卡托投影的正解公式如方程组(1)所示：

方程组中a为参考椭球体长半轴，b为参考椭球体短半轴，e为第一偏心率，e’为第二偏心率；λ₀为制定经线数，

为标准纬线数；

3.删除异常数据

根据原始数据利用Matlab绘图时会发现图中出现野点，这样异常的轨迹数据可以分为零散轨迹、多连轨迹、错位轨迹和缺点轨迹，以下为处理异常数据流程：

1)识别轨迹数据中是否存在相邻两航迹点时间间隔大于100s的情况，如果存在则判定为多连轨迹数据，将轨迹拆分成多条；

2)第一次计算每条轨迹中航迹点的个数并求平均值，以其十分之一为阈值，识别轨迹数据中是否存在航迹点数低于的情况，如果存在则判定为零散轨迹数据删除；由于此时很多航迹只有1-3个航迹点，这些零散数据影响平均航迹点数的取值，导致平均值结果偏小，所以后续还会做一次平均值计算，清理一遍零散轨迹数据；

3)识别轨迹数据中是否存在锐角的情况，如果存在则判定为错位轨迹数据，取第一个和第三个航迹点的空间坐标和时间坐标的中点，代替第二个航迹点，反复迭代多次，直到轨迹中没有锐角为止；

4)识别轨迹数据中是否存在相邻两个航迹点之间的距离大于10km的情况，如果存在则判定为缺点轨迹数据，当其时间间隔大于1s时，用线性插值方法以1s为步长插入1个等距等时的航迹点，当其时间间隔为0或1s时，插入20个的航迹点，时间选前一点的时间；

5)第二次计算每条轨迹中航迹点的个数并求平均值，以其十分之一为阈值，识别轨迹数据中是否存在航迹点数低于的情况，如果存在则判定为零散轨迹数据删除，考虑到第一次设定的阈值偏小，并且之后又有轨迹数据的拆分和插值，所以需要再次检验结果；

6)再次检查是否存在相邻两航迹点时间间隔大于100s，相邻三个航迹点之间的夹角是锐角和相邻两个航迹点之间的距离大于10km的情况；

至此，经过错误轨迹数据处理后的轨迹数据，数据准确，路径清晰，可作为空中交通流研究对象；

步骤三：对航迹进行聚类得到所需航路的交通轨迹；

自动聚类模型分两阶段：第一阶段是计算结构相似度，第二阶段是调用谱聚类算法，采用Normalized-cut，构建Normalized Laplacian矩阵，计算其特征值及特征向量，并求得谱系值，自动确定终端区飞行轨迹类别数目，然后对其进行聚类分析；

1.基于飞行轨迹的欧氏距离求解航迹结构相似度

飞行轨迹由轨迹点组成，经过数据预处理的飞行轨迹包含相同数目的飞行轨迹点，且投影后经纬坐标已转化为直角坐标系中的x值与y值，由此，顺序计算两条飞行轨迹中每对飞行轨迹点的欧氏距离，然后取平均值、最值，定义两条轨迹的欧氏距离，如公式(2)是计算两点间的欧氏距离，公式(3)是求解平均欧氏距离计算轨迹相似度的公式

式中：pi，pj代表两条飞行轨迹，均包含飞行轨迹点数N，pin是是飞行轨迹pi的第n个点，pjn同理，xin，xjn，yjn，yjn各是飞行轨迹pi，pj第n个点的横、纵坐标，Eu_D(pi,pj)是pi，pj的平均欧氏距离，dist(pi,pj)是两个飞行轨迹点的欧氏距离；

运用欧氏距离度量飞行轨迹相似性的优点在于操作步骤简单且易于理解，因而其应用范围比较广泛；缺点在于受噪声影响较大，如果未预先处理噪声点则欧氏距离难以准确度量飞行轨迹相似度，在模拟中若存在两条进场飞行轨迹pi与pj，各飞行轨迹点按顺序组成飞行轨迹点对，计算点对的欧氏距离，然后取平均距离定义为最终轨迹距离；

2.谱聚类算法

谱聚类算法的核心思想与图论相仿，是一种具有高计算性能的现代聚类方法，它依据无向图的多重分割理论对研究对象进行聚类分析，无向图G(V，E)的顶点V对应聚类问题中的数据点，用距离函数求得两点间的距离对边E赋权值aij，A的定义是待处理点的亲和度矩阵，也即G的邻接矩阵，图分割就是指把一个连通图的一些边切断，达到让原本连通的图打散成一些独立连通的子图的效果，Min-cut是最早的划分准则，但该方法存在一个缺点：对极小值很敏感，容易产生歪斜切割而导致结果不准确，因而，本方法选用Normalized-cut准则，引入度限制，平衡划分子类间的大小；

以下是谱聚类的流程：

输入：终端区飞行轨迹集合PS

输出：飞行轨迹簇C{C1,C2,.,Ck}

步骤1：计算结构相似度并构建亲和度矩阵A；

步骤2：构造矩阵Lsym；

步骤3：求解Lsym的前k个最小特征值λ1,λ2,.,λk以及其特征向量v1,v2,.,vk；

步骤4：计算谱系值，求出最大谱系值对应的下标k并作为聚类簇数目构造特征子空间X；

步骤5：最后利用经典聚类算法分析，得到聚类结果C{C1,C2,.,Ck}；

步骤四：根据航迹聚类结果求解时间序列给出所需时间段的流量预测数据；

1.空中交通流时间序列的构建

提取一段时间内通过某航路点的航空器数量来表示流量值，以构建空中交通流的时间序列，构建空中交通流的时间序列需要选取一个航路点a，取经过航路点a的一条交通流的轨迹f，它是由很多个航迹点fi构成的，取每个航迹点经过墨卡托投影转换后的横坐标，纵坐标和时间，然后计算每个航迹点与a的直线距离，距离最近的点fi即为航路点a的到达航迹点；

通过以上方法得到交通流到达航迹点集合，将所有到达航迹点的时间值提取出来，并进行升序排列，假设共有n个到达航迹点，得到单位为秒的到达时间数组，在此基础上，按照某一时间尺度对航班数量进行统计，得到交通流量的时间序列T_f＝[t_f1,t_f2,...,t_fi,...,t_fn’]，序列长度为n’，假设统计时间间隔Δt单位为分，总观察时间ts的单位为分，t_fj为第j个时间间隔内通过交汇点的航班数量；

记到达时间数组T中相邻两时间值为t_i、t_i+1,计算t_gi＝t_i+1-t_i，得到航空器时距序列，即微观交通流时间序列T_g＝[t_g1,t_g2,...,t_gi,...,t_g(n-1)]，序列长度为n-1，由此就得到了航迹的时间序列；

2.基于最大Lyapunov指数的混沌时间序列预测

基于上述理论，在判定空中交通流量时间序列具备混沌特性的基础上，从Lyapunov指数的定义出发，在得到时间序列后，重构相空间后经小数据量法求解出最大Lyapunov指数λ1，最后预测N+1时刻时间序列为公式(4)，其中x_M+1(j)代表相点x_M+1的第j维分量，x_M(j)，x_κ(j)，x_K+1(j)的定义类似，

通过求解λ1计算下一时刻的航路预测数据，与实际数据进行对比，识别空中交通拥挤流量态势，判断航路的拥挤情况，从而实现了对交叉航路的拥挤概率的预测。

Images