基于蒙特卡罗法的航路区与终端区冲突检测

王献锋，甘旭升，刘 飞，孙静娟

（1.西京学院理学院，西安 710123；2.空军工程大学空管领航学院，西安 710051）

0 引言

1966 年，Reich 提出一种内外两层的立方体飞行冲突模型，并充分考虑了管制的具体要求，引发了对冲突检测研究的小高潮［1］。1997 年，Erzberger等针对冲突概率估算问题提出了一种二维解析算法，并扩展到三维空间中［2］。同年，Fulton 等基于计算几何构建了Volonoi 多边形，并采用Delaunay 三角来降低多机冲突情况下的计算复杂性，但该方法受限于飞机速度的改变［3］。2000 年，KarineBlin 等在研究误差概率模型中引入位置误差，其能够在冲突检测中定位误差的影响，并最终确定飞机的动态位置［4］。国内相关研究虽然起步较晚，但也取得了诸多成果，李春锦等采用Reich 模型研究了平行航路的冲突检测问题，但没有在横向、纵向及垂直方向深入分解［5］。陈晨等研究了基于布朗运动的冲突概率算法，并分析了多种影响因素概率［6］。石磊等提出采用预测位置空间离散化的新算法来计算航路片段的瞬时冲突概率［7］。王凯等结合空管中的航路冲突标准，构建了双机基于管制规则库和决策树的航路冲突探测模型［8］。张中广等提出了基于ADS-B信息的终端区飞行冲突探测与告警算法［9］。纵观目前国内外文献可发现：它们通常集中讨论两架飞机之间或者一架飞机与其他多架飞机之间的冲突检测，且没有考虑不确定因素的影响，这一点在长期冲突探测中尤为重要，而这往往正是冲突检测成果的实用性体现。

有鉴于此，本文在随机微分方程位置预测基础上，提出采用蒙特卡罗法解决航路区与终端区的飞行冲突检测问题，并进行仿真验证。

1 飞行冲突探测模型

从飞行程序角度，飞机从机场起飞到机场降落，需要依次经过终端区-航路区-终端区。在此过程中飞机的飞行状态是不一样的，飞机在航路区飞行时，飞行高度和速度通常不会发生变化；而在终端区上，飞机由于要实施进近程序则会有较多飞行状态的改变。因此，本文主要研究飞机航路区飞行和终端区飞行的冲突探测。

飞行冲突探测可以分为3 步：第1 步是对信息的整合处理，这些信息包括飞行数据、航线信息、误差影响等等；第2 步是冲突探测过程，包括对飞机未来位置的预测，模拟计算是否会发生飞行冲突；第3 步是输出探测结果，用发生飞行冲突的实验次数除以总实验次数得到发生飞行冲突的概率值。具体流程图如图1 所示。

图1 飞行冲突探测流程

冲突探测模型即通常所说的保护区模型，根据安全间隔，每个飞机都会有自己的保护区域，其他飞机如果进入这个区域，则视为其他飞机入侵了我机的保护区，并将这种情况定义为发生了飞行冲突。

现在常用的冲突探测模型有圆柱型模型、椭球形模型（E 模型）和球形模型，如图2 所示。

假设两飞机的位置矢量为（x，y，z）和（x0，y0，z0）。圆柱形保护区水平截面圆半径为r，高度为h；椭球形保护区水平半长轴为a，垂直半短轴为c；球形保护区半径为r。则两飞机冲突的条件分别是：

图2 保护区模型

根据最小安全间隔距离的设置，发现最小水平间隔远大于最小垂直间隔，对比这3 种模型，球形保护区上任意一点到飞机（圆心）距离是一样的，很明显不利于飞行冲突的判定；圆柱形保护区和椭球形保护区都满足最小安全间隔的设定，但二者相比，椭球形保护区更加符合实际，运用该模型（E 模型）得到的冲突概率更加接近真实值，所以本次研究我们决定采用E 模型。

当在任意时刻两机位置关系满足冲突条件时，即潜在冲突机进入到目标机保护区内，则认为两机发生飞行冲突。为符合ATC 标准，取椭球体的长焦距a=5 n mile，短焦距c=2 000 ft。

一次冲突探测的主要判定流程为：两机沿预定航线飞行，其预估航迹的期望值受不确定性误差的影响。在每一单位时间长内，判断单位时间内潜在冲突是否发生在目标飞机的保护区内，若满足条件，则认为整个过程存在飞行冲突。

2 基于随机微分方程的位置预测

飞机在航路区的飞行高度普遍较高，容易受到高空风的影响，风往往被认为是对飞机影响最大的因素，它和其他的误差因素作用于飞机上，使飞机的飞行速度发生变化。假设在一段近似为零的时间内，由于外界因素影响，飞机的飞行速度变化符合高斯分布。构建一个xyz 坐标系，在t 时刻时（t 满足t∈T），飞机的位置方程为X（t）∈R3。则飞机的飞行速度变化可表述为如下数学表达式

式中，X（t）为t 时刻飞机的位置，μ（t）为t 时刻飞机的空速，ω（t）为均值为零的高斯分布的随机变量。

对该公式进行积分，则可以得到

其中，b（t）为布朗运动，由上式可知飞机位置是速度在时间上的积分附加一个布朗运动。

布朗运动应用于概率论的研究。物理学中对布朗运动有这样的描述：一个完全浸没于一种液体或者气体中的小粒子显示出的运动叫作布朗运动。假设一个随机过程{x（t），t≥0}为布朗运动，则该过程满足：①x（0）=0；②{x（t），t≥0}有平稳独立增量；③对于每个t＞0，x（t）服从高斯分布，且均值为0，方差为c2t。其中，c 为布朗运动方差的系数，当c=1 的时候，x（t）为标准布朗运动。

该模型假定飞机在三维平面飞行，则t 时刻飞行方向在水平面的投影与惯性坐标系x 轴正向夹角为θ（t），如图3 所示，飞机位置x（t）可以用随机微分方程表示［10］，即

其中，B（t）为标准布朗运动：

其中，σ2x、σ2y、σ2z分别是飞机水平航向、水平侧向和垂直方向速度扰动的功率谱密度，R（θ（t））为t 时刻的旋转矩阵，由于风对飞机水平速度和垂直速度的影响可以认为是独立的，因此，R（θ（t））为

3 冲突计算方法

在进行飞机位置预测时，会因为各种因素出现随机干扰，这会导致结果出现一定的偏差，所以在进行冲突探测时，一次实验得到的结果不能代表飞机发生冲突的概率，需要增加实验次数，测算一个接近真实概率的概率值。目前关于冲突探测的研究主要是对两种概率进行测算，一种是总体冲突概率，另一种是最大瞬时冲突概率。总体冲突概率因为考虑到了航迹上的每一个航路片段，所以相对最大瞬时冲突概率更全面，更具有说服力。在飞机飞行过程中，采样数量是有限的，不一定能够得到最优解，本文主要采用蒙特卡罗法计算总体冲突概率［11］。

蒙特卡罗法是一类随机方法的统称，而并不是一种算法的名称。它的特点是可以在随机采样上计算得到近似结果，随着采样的增多，得到的结果是正确结果的概率逐渐加大，但在（放弃随机采样，而采用类似全采样这样的确定性方法）获得真正的结果之前，无法知道目前得到的结果是不是真正的结果。也就是在采样不全时，通常不能保证找到最优解，这时蒙特卡罗方法就会认为采样的数目越多，也就越接近最优解。

4 仿真分析

根据航路区和终端区的不同特点，模拟不同场景，通过仿真实验预测飞行冲突概率。通过改变误差参数，与相同场景下不同误差条件下的冲突概率对比、分析和研究，找出其中相同之处以及区别所在，以验证所提出方法。

4.1 前期准备工作

4.1.1 飞机保护区模拟

在前面已经通过比较选定了适合的保护区，即椭球形模型（E 模型），根据ATC 标准，取椭球体的长焦距a=5 n mile，短焦距c=2 000 ft，将单位换算成m，分别是a=9 260 m，c=609.6 m。在MATLAB 上模拟该保护区模型，模拟椭球形模型如图3 所示。

图3 椭球形模型

4.1.2 飞机运动轨迹与保护区结合

由于研究对象为两架飞行航线存在冲突的飞机，不包括两架以上，故两架飞机的运动情况可以分两种情况：相向而行与同向而行。在发生的飞机空中相撞事件中基本都为相向而行，因为两架飞机相向而行，相对速度都比较大，飞机调整非常困难，如果没能提前预警进行规避，飞机相撞的可能性非常高，所以，我们的研究只考虑相向而行的情况。相向而行即两架飞机飞行方向相反，互相朝着对方飞行，该情况下两架飞机飞行速度无限制，高度需在一个高度层才会产生飞行冲突，使用MATLAB 软件模拟如下页图4 所示。

图4 相向飞行轨迹模拟

4.2 航路区飞行冲突探测

飞机在航路区飞行，可视为匀速直线运动，且飞机状态几乎不改变。误差因素只考虑风的影响，这里选择将飞机在未来的运动轨迹作为布朗运动，以其中一架飞机为参考系，将两架飞机的航迹分别看成是一条直线和一条不规则的曲线，做随机运动的飞机进入另一架飞机的保护求得概率即为飞行冲突概率。

选取乌伯林根空难作为背景，事发当日BTC2937班机和DHX611 班机的飞行高度都为11 000 m，且航线存在交点，即两机可能会相撞。空难发生前1 min，管制员发现情况，首先同BTC2937 班机取得联系，通知其飞行员降低高度300 m，以避免同DHX611 班机相撞。机组依照指挥下降高度，但飞机的空中防撞系统（TCAS）却提示拉高飞机，同时DHX611 班机按照TCAS 提示也开始下降高度。之后由于多方面因素，两架飞机及管制员都没意识到向着冲突的方向飞行，最终两机在10 068 m 高空相撞。

实验设定预警时间1 min，即管制员发现情况，两架飞机根据各自提示规避之后，两机位于同一高度层，且下降速率相同，两机飞行速度均为巡航速度900 km/h。设定风等因素影响之后，在MATLAB上进行5 000 次冲突模拟。

经模拟计算，航路区内两机飞行轨迹如图5 所示，其总体冲突概率与时刻关系（P-T）曲线如图6所示。调节风的误差系数再进行模拟，得到两机飞行轨迹如图7 所示，其P-T 曲线如图8 所示，也就是说，若存在风对飞机的误差影响，两机发生飞行冲突的概率会随着时间持续增至100 %。之前的背景是两架飞机最终相遇于一点，以该点为最终点往前推算两架飞机的起始点，为便于对比，还可再模拟一个两机未相撞但入侵了各自保护区的背景。此时，两机都处于同一高度层。速度设置为900 km/h，第1 架飞机航向为90，第2 架飞机航向为180，两机航向交叉角为90°。经过模拟后两机轨迹下页如图9 所示，从航迹预推图中可看出，并非所有模拟都发生了飞行冲突，其P-T 曲线如图10 所示。

图5 飞行轨迹

图6 P-T 曲线

从图中可以看出，两架飞机在1 min 之后发生飞行冲突的概率不足4%。修改风的影响系数，再进行模拟计算，得到的轨迹模拟如图11 所示，其P-T曲线如图12 所示，这种情况下，两架飞机1 min 之后发生飞行冲突的概率为12%。通过对比两种背景下的总体冲突概率可以发现，风对飞机的影响不可忽视，同时，也表明用蒙特卡罗法计算航路区的飞行冲突概率是可行的。

图7 修改风影响下的飞行轨迹

图8 修改风影响下的P-T 曲线

图9 两机90°交叉角下的飞行轨迹

图10 两机90°交叉角下的P-T 曲线

图11 两机90°交叉角下修改风影响的飞行轨迹

图12 两机90°交叉角下修改风影响的P-T 曲线

4.3 终端区飞行冲突探测

飞机在终端区的飞行状态多变，在这里只对处于同一边（上升或下降）的两架飞机进行冲突探测。

在2016 年上映的电影《火海凌云》中有一个场景，由科兹洛夫斯基饰演的实习飞行员古辛在降落过程中，由于一架航班飞错梯队而降落到古辛的班机本应降落的跑道上，并且该航班已经起飞，当时机场云量大且伴有浓雾，两架飞机已处于飞行冲突状态下，但是飞机的告警设备并没有提醒，若不避让就会发生相撞。古辛在降落中感知及时预判到该情况，及时操作飞机避让，继续绕飞一圈，避免了事故的发生。

将该事件作为终端区冲突探测的背景，假设两架飞机处于同一边上，第1 架飞机正在按照进近程序实施下降，第2 架飞机由于操作失误在同一边实施起飞程序。两机航向相对，即航向相差180°，由于终端区的飞机飞行状态多变，所以探测时间缩减为1 min，假设这1 min 内两架飞机的速度改变量很小，即表速相同，都为300 km/h。飞机在终端区的飞行间隔标准为垂直间隔100 m，纵向间隔3 000 m，横向间隔1 000 m，为提高预警概率，把横向的间隔也扩大为3 000 m。

两机飞行轨迹如图13 所示，其P-T 曲线如图14 所示。从图中可以看出，这两架飞机最初发生飞行冲突的概率是很低的，但是在最后阶段却突然升高，且发生冲突的可能性高达95%，此时管制员应及时下令让两机进行规避，以免发生飞行冲突甚至空中相撞。改变两机初始位置，再进行模拟实验，检验是否还会发生飞行冲突，飞机轨迹预推如下页图15 所示。由图可看出，1 min 后，两机预推点的距离仍比较大，在1 min 之内两机发生冲突的概率始终为0。可以看出，此情况下两机可正常进行上升或下降程序，不会发生冲突。

图13 终端区飞行轨迹

图14 终端区P-T 曲线

图15 终端区内改变初始位置的飞行轨迹

4.4 启示与意义

从以上两个不同区域的冲突探测研究中可以看出，蒙特卡罗法应用在冲突探测概率估算中比较具有说服力，并且探测时间范围比较广，可以帮助管制指挥人员预知发生飞行冲突的可能性。对于民航而言，可以将这一研究方法与机载预警设备（TCAS）相结合，通过将两种告警方法结合，以避免管制人员或者飞行人员因错误操作而造成飞行事故。同时，这两种方法也能互相辅助，大大地增加了预警的成功率，有效降低虚警率。对于军航而言，由于目前军机上没有安装机载预警设备，那么，这种方法可以成为后期军机改装的一项课题，从而让军机也具备了自动预警功能。而且这种方法也能够减轻指挥员的工作压力，使得指挥员能够提前对飞行员进行预警，能够有更多的时间来进行动作规避，提高了飞机的存活率。

5 结论

1）出于飞行冲突探测目的，选择采用随机微分方程构建适合的飞机位置预测模型，并通过飞行冲突探测仿真实例分析，证明了预测模型的有效性。

2）从真实案例中分别构建出航路区与终端区的飞行模拟场景，基于所选择的飞行冲突探测模型，通过蒙特卡罗法模拟计算飞行冲突的总体概率。仿真结果得出了较为准确的总体飞行冲突概，也说明了蒙特卡罗法应用于冲突探测研究的可行性，并获得了启示。

3）研究中只简单地讨论了风的影响，下一步将考虑更多因素，从而使研究背景更加真实，更加贴近实际。