目前面向空中交通管理(ATM)研究的轨迹预测方法主要分为短期轨迹预测和中长期轨迹预测。短期轨迹预测^[1-2]主要用于冲突探测与解脱。这类算法对飞机运动过程的建模非常精细，多数基于概率和最优滤波^[2]。由于其对预测精度的严格要求以及受限的运算速度，往往仅用于预测未来短时间内的飞行轨迹。中长期轨迹预测^[3-9](如未来1小时或者整个飞行过程)一般用于为ATM系统级模拟运行提供数据支持。中长期轨迹预测对单架飞机轨迹的精度要求不如短期轨迹预测法，但由于ATM系统往往包含多架不同型号飞机在广阔空域内同时运行，中长期轨迹预测需要具备以下几个特点^[5]： 快速、通用、覆盖范围广。

近年来，研究人员相继提出了一系列中长期轨迹预测方法^[5-9]，并建立了一些基于这些方法的飞行仿真系统 ^{[3-4, 7-8]}。其中，部分针对特定的机型或者特定的空域(如CTAS^[4]主要仿真终端区及临近空域运行)； 部分对飞行程序的处理比较简单，未考虑实际管制约束或者仅以直线段的前后连接来表示飞行轨迹。此外，文^[6-9]中的方法基于近似的二维平面坐标体系，这在进行大范围空间计算时会引进较大误差。

本文提出一种基于导航数据库^[10]和飞行计划的飞行全过程快速4-D轨迹预测方法。按照实际飞行计划，预测从起飞跑道到降落跑道的飞行全过程飞行轨迹。采用质点飞机模型保证方法的快速性和通用性。该方法考虑了专为空管仿真而提出的BADA^[11]模型参数与规则、导航数据库^[10]和空管约束，并在WGS-84地球坐标体系下进行空间计算，保证了预测轨迹具有较好的可信度和精度。

1 名词定义与假设

航段或程序段(leg)： 特指导航数据库中一段航段/程序段数据记录，其中存储了航段/程序段的类型、结构、约束等信息。

水平预测段(horizontal segment，HS)： 特指预测轨迹在水平方向上的一段投影，将所有HS前后顺次相连则得到完整预测轨迹在水平方向的投影，也可称为水平飞行剖面。

垂直预测段(vertical segment，VS)： 特指预测轨迹在垂直方向的一段投影，将所有VS前后顺次相连则得到完整预测轨迹在垂直方向的投影，也可称为垂直飞行剖面。由于速度计划与垂直剖面密切相关，本文将速度/时间信息整合在VS中。

预测轨迹段(segment)： 经过水平、垂直剖面整合以后得到的一段4-D预测轨迹，将所有预测轨迹段前后顺次连接则得到完整的4-D预测轨迹。

由于中长期轨迹预测属于开环预测，并不考虑飞行控制系统的作用。因此，如果在计算轨迹时考虑风速，仅需要将计算所得真空速与风速进行矢量求和得到地速即可。而实际上，多数中长期轨迹预测方法或者假设无风环境^{[6, 9]}，或者假定恒定风速^[7]，或者给定风场模型^[8]。这些假设对预测方法本身均无本质影响。本文限于篇幅，不再介绍有关风速的内容。

2 水平剖面预测

现代典型民航飞行全过程由标准离场程序(SID)、 途中航路(En-Route)、 标准进场程序(STAR)、 仪表进近程序(IAP)四部分组成，每一部分均由多段航段或程序段组成。

ARINC 424导航数据库规范^[10]中定义了23种航段类型，分为3类： 可完全预测类包括TF、 DF、 IF、 FD、 FC、 RF、 AF、 CF、 CD、 CR、 CI、 HF，可部分预测类包括CA、 FA、 HA，不可预测类包括VI、 VR、 VA、 VD、 FM、 VM、 HM、 PI。可部分预测类是由于其航段类型均为高度截止型，轨迹预测时只能大致估计飞机到达该高度的位置，在进行仿真飞行时，当检测到飞机到达截止高度后，需在线重预测剩余轨迹。不可预测类中以“V”开头的为航向保持类，该模式下轨迹受多种因素影响，不能提前预测； 以“M”结尾的为人工终止类，这一类由于需要管制员人工终止，也不可提前预测； PI航段由于其对飞机在航段内转弯时的定义并不精确，只有范围建议，同样不可提前预测。

水平剖面预测不涉及不可预测类以及等待模式下的3种航段(HF、HA、HM)，也忽略点类型航段IF。剩余12种航段根据其基本结构特征可进一步分为直线型(包括TF、DF、FD、FC、RF、AF、CF、CD、CR、CI、CA、FA)和圆弧型(包括AF、RF)。

水平剖面预测是一个循环过程，从SID程序的第一航段开始，每一步预测对应两相邻航段，分别设为前序航段和后续航段，完成一步预测之后，把当前后续航段设为下一步预测中的前序航段，并读取飞行计划中的下一航段，设为后续航段，开始下一步预测。当后续航段为IAP程序的最后航段时，预测终止，每一步预测中，根据前后两相邻航段的不同类型，可分为直线—直线型预测以及其他类型预测两大类。

2.1 直线—直线型预测

在一步水平预测中，如果前后两航段均为直线型航段，则该步预测为直线—直线型预测。如果两航段朝向相等或者接近(可设为低于某一阀值)，则不需要在航段之间插入衔接过渡段，只需生成一段与前序航段完全重合的水平预测段即可。如果前后两航段朝向有差异，需要根据导航数据库中规定的转弯规则设计转弯策略并预测转弯过渡轨迹。

2.1.1 切点转弯(fly-by)

当前序航段的终点和后续航段的起点重合且没有特别要求过点转弯时，使用切点转弯，以保证平滑过渡。典型的切点转弯如图 1所示。此时，在一步预测中需要构建2段HS； 一段如线段AB所示，与前序航段完全重合； 另一段为圆弧BC段。

计算圆弧段BC需要包含以下几个步骤：

步骤 1 估算本次转弯时的速度V_turn。当该点有速度限制V_lim时，令V_turn=V_lim； 否则根据该点高度限制查询速度计划表(见表 1)获取该点计划速度V_sche，令V_turn=V_sche。如果高度限制和速度限制均没有规定，则从起飞跑道开始，沿着已经预测的水平轨迹，假设固定爬升或下降坡度，估计转弯时飞机所处高度，并以此高度查询速度计划表，获取估计速度V_etm，令V_turn=V_etm。

表 1中，V_base=C_Vmin(V_stall)_To，V_cl,1、 V_cl,2、 Vd_CL,i(i=1,2,3,4,5)、 C_Vmin、 (V_stall)_To和M_cl均为飞机模型参数；Alt_trans和Alt_cruise分别代表等校正空速(CAS)等马赫数转换高度和巡航高度。

步骤 2 使用步骤1中估算的高度，根据“高度—横倾角”对应关系^[11]获得估计横倾角Φ。计算转弯半径：

$R=\frac{{{V}^{2}}_{turrn}}{g~tan\phi }.$

步骤 3 设转弯圆心经纬度坐标为(LON,LAT)，利用式(1)将其转换到“北—东—地”坐标系中，设坐标为P_c(x,y,z)

$\begin{align} & x=cosLAT~cosLON, \\ & y=cosLAT~sinLON, \\ & z=sinLON. \\ \end{align}$

(1)

求解式(2)计算P_c(x,y,z)：

$\begin{align} & {{N}_{1}}{{P}_{c}}=\frac{{{R}_{n}}}{{{R}_{e}}}, \\ & {{N}_{2}}{{P}_{c}}=\frac{{{R}_{n}}}{{{R}_{e}}}, \\ & \left\| {{P}_{c}} \right\|=1. \\ \end{align}$

(2)

其中: N₁和N₂分别代表前和后航段所在大圆面在“北—东—地”坐标系中单位法向量； R_n为航段衔接点所在处地球半径，R_n=R_e[1－sin²LAT]，R_e代表地球标准椭球模型长半轴半径。

步骤 4 按式(3)计算转弯起点坐标P_start和转弯终点坐标P_end，再利用式(1)计算两点经纬度。

$\begin{align} & {{P}_{start}}={{N}_{1}}\times ({{P}_{c}}\times {{N}_{1}}), \\ & {{P}_{end}}={{N}_{2}}\times ({{P}_{c}}\times {{N}_{2}}). \\ \end{align}$

(3)

之后将点C设为下一步预测的起点。

2.1.2 过点转弯(fly-over)

过点转弯发生在前后航段未相连(即前序航段终点与后续航段起点不重合或者后续航段起点不确定)或者当导航数据库特别要求过点转弯时。执行过点转弯时，飞机必须飞过前序航段的终点，在此以前不允许转弯。本文采用三步法来预测过点转弯要求下的过渡轨迹。

步骤 1 判断前后航段是否空间相连，如果是则转步骤3，否则继续步骤2。

步骤 2 将前后航段按照图 2中的方式延长，获得虚拟连接点。之后按照切点转弯中的方法设计一段转弯圆弧。如果求出的转弯起点或终点不满足实际需求(如转弯起点在前序航段终点以前，转弯终点离后续航段终点太近甚至在后续航段终点之后等情况)，转步骤3，否则结束本步水平预测。

步骤 3 图 3中的过点转弯过程包含2段转弯(弧BDE与弧FH)及一段直线(线段EF)。转弯开始于前序航段终点(点B)，之后引导飞机切入后续航段，切入角度与后续航段夹角为φ，最后通过圆弧FH完成切入过渡。

按照节2.1.1 中介绍的方法估算转弯半径，并预测过渡轨迹如下：

1) 计算线段BC₁在点B切向方向θ_BC1=θ_pre－90°,θ_pre为前序航段在点B切向方向，根据点B坐标P_B、 R和θ_BC1，采用文^[12]中已知参考点位置，求固定方向上一定距离外目标点坐标的方法，计算点C₁坐标P_C1。

2) 计算大圆线₁E在点₁的切向方向θ_C1E=θ_post+90°-Φ,θ_post为后续航段朝向。根据R、 θ_C1E和P_C1，采用文^[12]中已知参考点位置，求固定方向上一定距离外目标点坐标的方法，求得点E坐标P_E。

3) 设后续航段终点坐标P_J，线段EF方向θ_EF=θ_post－Φ。根据P_E、 P_J、 θ_EF和θ_post，采用文^[12]中求两大圆线交点直角坐标的方法，求得点G坐标P_G。

4) 计算线段GC₂在点G的方向θ_GC2=θ_post+90°－Φ/2，根据R、 θ_GC2和P_G，采用文^[12]中已知参考点位置，求固定方向上一定距离外目标点坐标的方法，求得点C₂坐标P_C2。

5) 设直线EF所在大圆面法向量N_EF=P_E×P_G，后续航段所在大圆面法向量N_GJ=P_G×P_J，利用式(3)计算点C₂在两线上的投影点F、 H的坐标P_F、 P_H，将点F设为第二次转弯的起点，点H设为第二次转弯的终点。

上述过程共生成4段水平预测段，分别为AB(直线段)、 BDE(圆弧段)、 EF(直线段)和FH(圆弧段)。然后，将点H设为下一步预测起点。

2.2 其余类型预测

除节2.1中介绍的直线—直线型预测以外，其余3种前后连接方式(包括直线—圆弧型、圆弧—直线型、和圆弧—圆弧型)均不需要设计专门的过渡段，在一步水平预测中只生成一段与前序航段完全重合的水平预测段即可。

3 垂直飞行剖面预测

典型民航飞机垂直剖面分为爬升、巡航、下降3个阶段。其中，爬升阶段和下降阶段存在很大程度的对称性。另一方面，本文与文^[3-9]相同，假设巡航过程为定高定速巡航。因此，本节只介绍爬升阶段的垂直剖面预测方法。

爬升过程中，飞机真空速和飞行航线角(FPA)都在不断变化。即使更为精确的机载飞管系统(FMS)^[13-14]中的航迹预测功能也需要以数值积分的形式做近似计算。本文为了保证预测方法的通用性，并参照文^[11]中的规则，在FMS预测方法的假设基础上做出以下近似：

1) 假设飞机爬升速度计划表如表 1所示，即按照高度分段等CAS爬升。

2) 在假设1的基础上进一步假设每段等CAS高度区间内，飞机按照等真空速(TAS)爬升，其 TAS值取起始高度和终止高度处TAS值的平均值。

3) 每段等速区间内的爬升率相等。

爬升垂直剖面预测如图 4所示，由以下几步构成：

步骤 1 初始剖面构建，从起飞跑道出口高度开始，直到设定的巡航高度，按照表 2中的速度计划，在每个等CAS高度区间内根据对应的CAS值和起点、终点高度值计算TAS_start、 TAS_end，和TAS_aver=(TAS_start+TAS_end)/2。再利用TAS_aver按照文^[11]中节3.2中的方法计算该段FPA。

步骤 2 加入平飞加速段，根据表 1中的速度计划，1 829 ~3 048 m段与下一段速度相差较大，如果与其他各段一样直接假设速度突变可能引起较大误差。因此，在3 048 m处加入平飞加速段表征V₆到V_cl,2之间的加速过程(见图 4b中线段AB)，加速度按照文^[11]中式(3.2-1)计算。

步骤 3 速度限制调整，如果步骤2预测剖面中某段速度违反该段中某定位点(速度限制违反点，见图 4c中五角星)的速度限制，则对该点以前所有TAS值高于该限制值的垂直计划段修改其原定TAS至该限制值。修改过后，受影响垂直预测段需要重新计算FPA。如果由此造成速度限制违反点处预测高度变化，还需按照步骤1中的方法逐段修改其后续垂直预测段。

步骤 4 高度限制调整，如果步骤3预测剖面违反某定位点的高度限制，则以限制高度值做水平线，将其与垂直剖面交点(见图 4d中点C)与高度限制违反点(见图 4d中点F)之间的轨迹由爬升改为平飞，之后修改高度限制违反点之后垂直剖面。如此顺次检查爬升剖面中各点，直到所有定位点高度限制都被满足。

4 水平垂直剖面整合

图 5中，将各水平预测段(h₁h₂、h₂h₃、…、h₈h₉)按照其水平方向长度投影到水平距离坐标轴上，例如h₁h₂投影到hd₁hd₂。各段端点都对应水平距离坐标轴上的一个投影点，如h₁对应hd₁，h₂ 对应hd₂。同理将各垂直预测段(v₁v₂、v₂v₃、…、v₅v₆)按照其水平方向距离投影到水平距离坐标轴上，例如v₁v₂投影到vd₁vd₂。

最后根据水平距离轴上所有投影点到零距离点的距离，在相邻两点之间构成一段预测段。例如： vd₂hd₂段在水平方向特征取自h₁h₂段，在垂直方向上特征取自v₂v₃段，因此该段为直线平飞段。

5 仿真实验

本文以北京首都(ZBAA)机场到上海虹桥(ZSSS)机场飞行全过程为例，检验所提轨迹预测方法的可行性和适用性。实验所用飞行计划信息如表 2所示。

试验中采用BADA 3.10^[10]版本中的B737-300飞机模型，设定巡航高度9 878 m，爬升关键速度采用250节/290节/0.79马赫模式，即V_cl,2=290节，M_cl=0.79马赫，min(V_cl,1,250节)=250节。

试验在配置Intel Core i3-2120，3.3 GHz处理器，4 GB内存的PC平台上运行，预测全程航迹仅耗时0.5 s。试验中一共产生31段水平预测段，26段垂直预测段，整合之后一共50段预测段。预测总飞行时间1 h 34 min 33 s，与实际统计数据高度吻合。表 3列出本文预测的飞行过程中部分关键点到达时间与实际记录中到达时间的对比。可以看出3个关键点的到达时间差均小于10 min，考虑到实际样本数据受到风速，流量管制，偏航等因素影响，本就与仿真数据存在一定差距，本方法预测轨迹精度完全能够满足ATM中系统级流量调度和空域规划研究的需求。

限于篇幅，本文不对全过程预测结果进行详细分析，只列出爬升阶段的预测结果。图 6a中白色曲线为整个SID程序和开始巡航段的水平剖面，标识“TOC”处为爬升顶点，飞机到达此处后结束爬升过程，开始进入巡航阶段。可以看到TOC所处区域仍在SID程序范围内。图 6b中所示只对应起飞跑道到TOC之间的垂直剖面。表 4详细列出图 6中2个剖面整合后的预测飞行段信息。其中，A/C代表圆弧爬升型预测段； S/C代表直线爬升型预测段； S/P代表直线平飞型预测段。

6 结 论

本文提出一种基于飞行计划和导航数据库^[13]的飞行全过程快速4-D轨迹预测方法。相比机载飞行管理系统的轨迹预测功能，本方法为了提高运算速度和通用性，在飞机模型、速度计划和FPA等方面做了一定的简化。但由于在计算中引入了导航数据库^[10]，考虑了真实的飞行程序约束，采用了BADA模型^[11]，并且在空间计算中采用符合WGS-84坐标体系的运算方法，使得预测的飞行轨迹可以满足空管研究的精度要求。下一步工作将继续完善垂直剖面预测的功能，使其适应于更多更复杂的终端区飞行情况。