飞行冲突是指飞行过程中飞机之间的距离小于规定的最小间隔，一般分为水平间隔和垂直间隔，相当于在飞机外面形成了一个圆柱形保护区，如果有其他的飞机进入了这个保护区，则形成飞行冲突。低空飞行环境中的空域形态复杂，飞行器密度高，监视系统难以有效覆盖该空域，再加上低空空域内活动的飞行器航线自由度大，导致低空飞行冲突发生的可能性大大增加，飞行安全受到极大威胁。为此要有更加高效的冲突探测和避免手段，以保证低空飞行的安全性。
飞行冲突探测是在航迹预测基础上，建立2架或多架飞行器航迹误差模型和飞行碰撞危险模型，计算飞行冲突概率。冲突概率的计算精度决定了预警的准确程度，精度越高准确度越大，飞行冲突探测目前主要采用概率的方法^[1]。概率冲突探测方法考虑了飞行器行为的不确定性，采用不确定性模型描述潜在的航迹位置偏差，建立飞行安全空间模型，计算冲突发生的概率。该方法冲突探测准确率较高，漏警率和虚警率较低，因此应用较为广泛。通常航迹预估是在预定航迹上加一个扰动项，对扰动的建模精度决定了航迹预估的精度，因此航迹预估的难点主要在于对飞行器飞行的干扰因素进行建模。对于低空自由飞行而言，飞行器速度相对较低，且没有航路历史飞行统计数据，因此需要对预测误差进行精细建模，本文把低空自由飞行的误差来源具体地分为3项：导航误差、控制误差、风扰动，对每一项误差来源进行单独建模和分析。
经典的冲突探测概率算法如Paielli和Erzberger (P&E)^[2-3]的近似算法，假设了飞行速度大小和方向均不改变，只适用于直线飞行，不适用于机动飞行的情况。在P&E近似算法的基础上，很多人做了一些扩展性的研究^[4-6]。Prandini等^[7-8]提出的基于布朗运动的算法，把风等不确定因素对飞行的影响计入动力学方程，来计算中短期的冲突概率。文献[9-11]针对不同的情况下改进了基于布朗运动的算法。Jilkov等^[12]在文献[13]的基础上从航迹预测方面采用多模型来提高预测精度，从而提高了冲突概率的估算精度。Jacquemart和Morio^[14]采用IPS(Interacting Particle System)算法来提高冲突概率的估算精度，改进了算法的收敛速度。对于机动飞行情况下的冲突探测，Monto Carlo^[15-18]算法虽然能得到较好的结果，但计算量大，计算过程复杂，无法满足实时性要求。
本文基于概率型冲突探测的思想，假设可能处于飞行冲突中的飞行器之间能够互相交换飞行状态信息，对2架飞行器的航迹进行预估，分析了飞行器航迹预测的3种误差来源并分别建模。考虑到低空飞行中飞行器种类的多样性以及相对于民航的低速性，把飞行器的保护区建模为椭球。采用坐标变换以及扩展冲突区域的方法，提出了在机动飞行情况下，计算冲突概率的近似解析算法。可以计算短期内任意时刻飞行器之间的瞬时冲突概率，以及冲突概率随时间的变化，适用于低空飞行机动目标的冲突探测。最后通过仿真实验对算法进行了验证和比较，说明了算法的准确性，快速性和实时性。
1 飞行冲突保护区 对于低空自由飞行而言，主要活动的飞行器为直升机，飞行速度相对较低，而且低空航线自由，经常会有垂直方向的机动，也容易遭遇起飞和降落的其他飞行器，故而低空飞行的垂直间隔和水平间隔应当具有同样量级的要求，本文以直升机为对象，采用椭球作为低空飞行器的保护区来进行冲突探测比较合适。即

(1)

式中：S_x(t)为t时刻飞行器的保护区; x(t)为t时刻飞行器的位置; U=diag{1/σ₁², 1/σ₂², 1/σ₃²}，σ₁、σ₂和σ₃分别为飞行前后、左右和垂直间隔，也是保护区的3个半轴。通常情况下，民航飞行中水平间隔为6 km，垂直间隔为0.3 km，即σ₁=σ₂=6 km, σ₃=0.3 km如果一架飞行器B进入另一架飞行器A的保护区即发生飞行冲突(如图 1所示)，因此在给定的前瞻时间[0, T_f]内任意时刻的冲突概率P_C(t)定义如下：

(2)

图 1 低空自由飞行的保护区 Fig. 1 Protection zone of aircraft for free flight in low-altitude airspace

图选项

式中:x^B(t)为t时刻飞行器B的位置。
2 航迹预估模型 冲突概率的计算依赖于飞行器的航迹预估模型，在飞行过程中飞行器受各种随机扰动的影响，需要对扰动因素进行分析建模，建立适当的误差模型是冲突探测的基础，直接影响冲突探测的准确性。
位置预测误差的来源主要有以下几个方面：① 导航误差(定位误差、速度测量误差)会导致位置预测误差；② 控制误差，飞行器做机动飞行时需要保持预期的过载，但由于油门扰动、控制面扰动、反馈系统噪声等影响会产生控制误差并随时间积累，进而导致位置预测误差；③ 风扰动对位置预测也有影响，通常很难建立一个确定的风扰动模型。位置预测误差是由这些因素综合作用的结果，通常可认为它服从期望值为零的高斯分布。
假设飞行器之间可以实时交换状态信息，包括当前时刻的位置、速度等信息。数据链系统以较短的周期更新状态信息，当飞行器获得新的状态信息时，记当前时刻为零时刻。xⁱ(t)、、分别表示地面坐标系中，飞行器i=(A, B)在t时刻后的真实位置、预测位置及预测位置误差，三者之间的关系如下：

(3)

t时刻预测位置可由当前时刻的位置信息、运动信息、控制信息外推得到：

(4)

(5)

式中：和分别为当前时刻导航仪器测得的位置和运动速度; 为预期控制输入; 、、和分别为定位误差、速度测量误差、控制误差、随机风扰动引起的位置预测误差分量；Rⁱ为机体坐标系到地面坐标系的旋转矩阵，i=(A, B)。

(6)

式中：ψ、θ和?分別为飞行器的偏航角、俯仰角和滚转角。
位置预测误差的协方差矩阵为

(7)

位置预测误差服从期望值为零的高斯分布，即。
卫星定位系统在测量飞行器位置时存在误差，该误差通常可描述为期望值为零的高斯分布，定位误差协方差矩阵为

(8)

式中：为功率谱矩阵，下文中所有Σ皆为相应上下标所指对象的功率谱矩阵。
导航仪器在测量飞行器飞行速度时存在误差，由式(4) 可知，速度测量误差会引起位置预测误差，并随时间累积。速度测量误差协方差距阵为

(9)

飞行器做机动飞行时需要保持预期的过载，但由于油门扰动、控制面扰动、反馈系统噪声等影响会产生控制误差，在短时间内，控制误差可以描述为作用于加速度的布朗运动。布朗运动的线性变换仍然为布朗运动^[19]，可得如下表达式：

(10)

式中：B_aⁱ(t)为三维标准布朗运动。布朗运动的积分为高斯过程，任意时刻的协方差矩阵为^[19]

(11)

文献[20]提出了空间相关风力场，描述风对位置预测误差的影响。认为风由确定部分和随机部分组成，确定部分是均匀、各向同性的，可以通过测量补偿抵消，因此，只关心随机部分对位置预测误差的影响。随机风扰动作用在飞行器上，使飞行器的空气动力产生扰动，进而引起加速度变化。将随机风对加速度的扰动描述为布朗运动，仿照式(11) 的推导可得

(12)

式中：为功率谱矩阵。
根据高斯分布的性质：有限个服从高斯分布的随机变量的线性组合仍然服从高斯分布。位置预测误差是由导航误差、随机风扰动、控制误差等因素引起的，可认为三者是相互独立的，因此，飞行器在地面坐标系内的位置预测误差协方差矩阵为

(13)

每个飞行器的位置预测误差描述为高斯分布，因此可以将2个飞行器各自的协方差矩阵合并为2个飞行器之间相对位置误差的联合协方差矩阵。将联合协方差矩阵施加在其中的一个飞行器上，称为“随机”飞行器，则另一架飞行器可以看作是没有误差的，称为“参考”飞行器。将A作为随机飞行器，相对位置固定在原点；将B看作参考飞行器，相对于随机飞行器运动。
两飞行器在地面坐标系内的相对位置真实值为

(14)

两飞行器在地面坐标系内的相对位置测量值为

(15)

两飞行器在地面坐标系的相对位置误差为

(16)

联合协方差矩阵为

(17)

其中:Q^AB(t)为交叉相关项，定义如下:

(18)

导航误差与控制误差来源于两飞行器上两套互相独立的导航设备与控制系统，因此两者引起的交叉相关项为零。但随机风扰动引起的误差交叉相关项在近距离内不可以忽略^[20]。

(19)

Blin等^[21]推导了二维平面内随机风扰动引起的交叉相关项。将其扩展为三维情况:

(20)

式中：eⁱ=(xⁱ, yⁱ, zⁱ)^T为飞行器i=(A, B)的单位方向向量。
3 冲突探测算法 由于预测位置相对于真实位置的不确定性带来的误差符合特定的概率分布，故可依据冲突探测算法预测前瞻时间[0, T_f]内任意时刻发生冲突的概率，并以此时间段内冲突概率的最大值作为衡量冲突发生可能性的依据。
Paielli和Erzberger^[2-3]基于常速度和直线飞行假设，使用平行于相对速度方向的无限延长的条形扩展区域计算冲突概率，这种方法得到直线航路内冲突概率的上限值，高估了冲突概率且无法用于机动飞行情况。其他****使用Monto Carlo算法计算冲突概率，虽然适用于计算机动情况下的冲突概率，但计算量大，计算过程复杂，无法满足实时性要求，本文仅使用Monto Carlo算法作为对比算法验证。
由于联合误差分布不规则，坐标之间是耦合的，积分运算十分困难，因此需要进行适当的变换，使变换后的误差分布在新坐标系下服从独立的高斯分布。
设εⁱ(t)为新坐标系O-ξηγ中的位置向量，T为从地面坐标系O-xyz到新坐标系O-ξηγ的旋转矩阵：

(21)

新坐标系中相对位置向量为

(22)

新坐标系中联合协方差矩阵为

(23)

联合协方差矩阵为实对称阵，通过矩阵正交分解可得正交矩阵T，使得联合协方差矩阵为

(24)

式中：σ_ξ²、σ_η²、σ_γ²为方差在坐标轴方向的特征分量。由此可得到位置预测误差在坐标系O-ξηγ内的联合概率密度函数：

(25)

在新坐标系O-ξηγ内，飞行冲突保护区为

(26)

式中：W=T^－1。
任意时刻的冲突概率等于联合概率密度函数在变换后的保护区内的积分

(27)

尽管概率密度函数是关于坐标解耦的，但求解该三重积分表达式的解析解十分困难，因此寻找一种精度高且满足实时性要求的算法是问题的关键所在。首先考虑二维冲突情况，即假设两飞行器在同一高度飞行，然后将二维冲突探测算法拓展到三维空间。
设经坐标变换后的保护区二次型矩阵为

(28)

式中:W、U为退化的二阶方阵。
图 2为二维平面冲突概率计算的示意图，其中保护区退化为平面椭圆形，其边界为

(29)

图 2 二维平面冲突概率计算示意图 Fig. 2 Schematic diagram of calculation of conflict probability in two-dimensional plane

图选项

其外切矩形的边长分别为

(30)

(ξ_i, η_i)、δξ_i、δη_i可由式(31) 和式(32) 确定:

(31)

(32)

式中：。
由于概率密度函数关于坐标解耦，且扩展区域的边界平行(垂直)于坐标轴，假设联合概率密度函数在扩展区域内的积分为

(33)

令

(34)

式中：为标准正态分布的累积分布函数。
则式(27) 的冲突概率为可表示为

(35)

由于P_extend(t)无法求得解析表达式，本文用P(t)来近似代替P_C(t)，造成的误差为P_extend(t)。
现在将二维冲突探测算法拓展到三维空间，如图 3所示。

图 3 三维空间冲突概率计算示意图 Fig. 3 Schematic diagram of calculation of conflict probability in three-dimensional airspace

图选项

设经坐标变换后保护区的二次型矩阵为

(36)

保护区的边界为

(37)

其外切长方体的边长分别为

(38)

(ξ_i, η_i, γ_i)、δξ_i、δη_i、δγ_i可由式(39) 确定:

(39)

(40)

式中: 。
三维冲突概率的近似表达式可仿照式(34) 推导：

(41)

误差为。
4 算例仿真 4.1 直线相向飞行 考虑如图 4(a)所示的相向飞行情况。由于只需考虑2架飞行器之间的相对运动，假设飞行器A在相对坐标系中速度为0，即固定于原点，飞行器B将以常相对速度接近A。飞行器B在相对坐标系中的初始位置为x^B=[－9, －5, 0]^T km，相对速度v^B=[0.3, 0, 0]^T km/s。考虑到短期冲突预测的实时性要求，现选取[0, 60]s作为冲突预测的前瞻时间，分别使用P&E近似算法、Monte Carlo仿真算法(抽样10⁵次)及本文算法计算[0, 60]s内的冲突概率，计算结果如图 4(b)所示。当t较小(t < 10 s)时，由于A、B两飞行器距离较远，Monte Carlo仿真算法计算的瞬时冲突概率较小。t较大(t > 50 s)时，同样如此。t=30 s时两飞行器之间相对距离最小，冲突概率达到最大值0.27，P&E近似算法与本文算法计算结果分别为0.35、0.29。本文算法与Monte Carlo仿真算法的最大误差为0.02，相对误差为7.4%。与P&E近似算法相比，本文算法更符合Monte Carlo仿真算法结果的变化趋势。

图 4 直线相向飞行结果比较 Fig. 4 Result comparison of straight toward flight

图选项

4.2 转弯机动飞行 考虑如图 5(a)所示的定速转弯机动情况。假设飞行器A在相对坐标系中速度为0，即固定于原点，飞行器B将做机动飞行。由于P&E近似算法是基于常速度假设的，因此在使用此算法时需将飞行器B的速度在时间段[0, 60]s内做一个分段线性假设，即假设在每一个小时间段内飞行器B的速度大小是不变的，本例中共分为7个小时间段[0, 11]，[11, 15]，[15, 21]，[21, 24]，[24, 28]，[28, 31]，[31, 60] s。飞行器B在相对坐标系中的初始位置为x^B=[－13, －3.8, 0]^T km，以水平速度v^B=[0.3, 0, 0]^T km/s飞行11 s后进行逆时针定速转弯，转弯速率为 rad/s，20 s后以速度v^B=[0, 0.3, 0]^T km/s继续水平飞行。分别使用3种算法计算[0, 60]s内的冲突概率，计算结果如图 5(b)所示。当t较小时，由于A、B两飞行器距离较远，Monte Carlo仿真算法计算的瞬时冲突概率较小，但P&E近似算法计算的冲突概率却非常大，这是因为P&E近似算法是在相对速度方向无限延长冲突扩展区域，在t=16 s左右，相对速度方向指向飞行器A。t=31 s时，两飞行器之间相对距离最小，冲突概率达到最大值，Monte Carlo仿真算法与本文算法计算结果分别为0.25、0.28。本文算法与Monte Carlo仿真算法的最大误差为0.03，相对误差为12%。本文算法的计算结果与Monte Carlo仿真算法结果的变化趋势相同。

图 5 二维情况转弯机动飞行结果比较 Fig. 5 Result comparison of turning maneuver flight in two-dimensional case

图选项

4.3 三维机动飞行 考虑如图 6(a)所示的三维机动情况。假设飞行器A在相对坐标系中速度为0，即固定于原点，飞行器B做如图 6(a)所示的爬升转弯机动飞行。飞行器B在相对坐标系中的初始位置为x^B=[－15, －3.8, －0.1]^T km，飞行速度v^B=[0.3, 0, 0]^T km/s。以垂直爬升率20 m/s爬升5 s后继续水平飞行，11 s后进行逆时针定速转弯，转弯速率为 rad/s，20 s后以速度v^B=[0, 0.3, 0]^T km/s继续水平飞行。由于P&E近似算法是基于常速度的直线飞行假设，故不适用于本例分析，现分别使用Monto Carlo仿真算法与本文算法计算时间[0, 60]s内的冲突概率，计算结果如图 6(b)所示。Monte Carlo仿真算法与本文算法计算的最大概率分别为0.36、0.41，最大误差为0.05，相对误差为13.9%。本文算法的计算结果与Monte Carlo仿真算法结果的变化趋势相同，且基本吻合。

图 6 三维情况爬升转弯机动飞行结果比较 Fig. 6 Result comparison of climbing turning maneuver flight in three-dimensional case

图选项

在同样的条件下对3种算法计算耗时的比较结果见表 1。P&E近似算法计算简单，耗时少，但高估了冲突概率，也不适用于机动飞行。本文算法的计算准确性大于P&E近似算法，且计算时间远远小于Monte Carlo仿真算法，满足实时性要求。
表 1 3种算法计算时间比较 Table 1 Computation time comparison among three algorithms

算法	计算时间/ms
	直线相向飞行	转弯机动飞行	三维机动飞行
Monte Carlo仿真算法	7.039×103	7.039×103	1.047×104
P&E近似算法	0.06	0.48
本文算法	2.58	2.58	3.72

表选项






5 结论 1) 分析了航迹预估中导航误差、控制误差以风扰动3种误差来源并分别进行了建模，建立了适当的误差模型，提高了航迹预估的精度。
2) 提出了低空自由飞行下冲突概率探测方法的近似解析算法，可用于计算三维空间内的冲突概率。仿真结果表明，本文算法适用于机动飞行的情况，可实现复杂环境下的冲突探测，具有足够的计算精度并满足实时性要求。

参考文献

[1]	KUCHAR J, YANG L. A review of conflict detection and resolution modeling methods[J].IEEE Transactions on Intelligent Transportation Systems, 2000, 1(4): 179–189.DOI:10.1109/6979.898217
[2]	PAIELLI R, ERZBERGER H. Conflict probability estimation for free flight[J].Journal of Guidance, Control, and Dynamics, 1997, 20(3): 588–596.DOI:10.2514/2.4081
[3]	PAIELLI R, ERZBERGER H. Conflict probability estimation generalized to non-level flight[J].Air Traffic Control Quarterly, 1999, 7(3): 195–222.DOI:10.2514/atcq.7.3.195
[4]	邓炜, 张军, 吴限, 等. 一种适用于航路改变情况的冲突概率探测算法[J].北京航空航天大学学报, 2005, 31(12): 1327–1332. DENG W, ZHANG J, WU X, et al. Algorithm of conflict probability prediction for the case of trajectory change[J].Journal of Beijing University of Aeronautics and Astronautics, 2005, 31(12): 1327–1332.DOI:10.3969/j.issn.1001-5965.2005.12.013(in Chinese)
[5]	高杨, 徐浩军, 郑海峰. 计算航路上飞行冲突概率的一种方法[J].飞行力学, 2009, 27(2): 50–53. GAO Y, XU H J, ZHENG H F. Algorithm to estimate the conflict probability on the airway[J].Flight Dynamics, 2009, 27(2): 50–53.(in Chinese)
[6]	沈笑云, 周波, 曹博, 等. 基于冲突概率的低空自由飞行冲突检测算法[J].电光与控制, 2014, 21(6): 43–47. SHEN X Y, ZHOU B, CAO B, et al. A free flight conflict detection algorithm of low-altitude airspace based on conflict probability[J].Electronics Optics & Control, 2014, 21(6): 43–47.(in Chinese)
[7]	PRANDINI M, HU J. A probabilistic approach to aircraft conflict detection[J].IEEE Transactions on Intelligent Transportation Systems, 2001, 1(4): 199–220.
[8]	PRANDINI M, PUTTA V. A probabilistic measure of air traffic complexity in three-dimensional airspace[J].International Journal of Adaptive Control and Signal Processing, 2010, 24(10): 813–829.DOI:10.1002/acs.1192
[9]	李丹, 崔德光. 基于布朗运动的空中交通短期冲突探测[J].清华大学学报(自然科学版), 2008, 48(4): 477–481. LI D, CUI D G. Air traffic control conflict detection algorithm based on Brownian motion[J].Journal of Tsinghua University(Science and Technology), 2008, 48(4): 477–481.(in Chinese)
[10]	刘小龙, 罗以宁, 赵喜求, 等. 一种改进的Prandini概率型中期冲突探测方法[J].计算机技术与发展, 2013, 23(1): 214–216. LIU X L, LUO Y N, ZHAO X Q, et al. An improved medium-term conflict detection methods of Prandini probability type[J].Computer Technology and Development, 2013, 23(1): 214–216.(in Chinese)
[11]	石磊, 吴仁彪, 黄晓晓. 基于总体冲突概率和三维布朗运动的冲突探测算法[J].电子与信息学报, 2015, 37(2): 360–366. SHI L, WU R B, HUANG X X. Conflict detection algorithm based on overall conflict probability and three dimensional Brownian motion[J].Journal of Electronics & Information Technology, 2015, 37(2): 360–366.DOI:10.11999/JEIT140363(in Chinese)
[12]	JILKOV V P, LI X R, LEDET J H.Improved estimation of conflict probability for aircraft collision avoidance[C]//17th International Conference on Information Fusion.Piscataway, NJ:IEEE Press, 2014:1-7.
[13]	BLOM H, BAKKER G.Conflict probability and incrossing probability in air traffic management[C]//Proceedings of the 41st IEEE Conference on Decision and Control.Piscataway, NJ:IEEE Press, 2002, 3:2421-2426.
[14]	JACQUEMART D, MORIO J. Adaptive interacting particle system algorithm for aircraft conflict probability estimation[J].Aerospace Science and Technology, 2016, 55(5): 431–438.
[15]	VISINTINI A, GLOVER W, LYGEROS J. Monte Carlo optimization for conflict resolution in air traffic control[J].IEEE Transactions on Intelligent Transportation Systems, 2006, 7(4): 470–482.DOI:10.1109/TITS.2006.883108
[16]	WOLF T B, KOCHENDERFER M J. Aircraft collision avoidance using Monte Carlo real-time belief space search[J].Journal of Intelligent & Robotic Systems, 2011, 64(2): 277–298.
[17]	BELKHOUCHE F. Modeling and calculating the collision risk for air vehicles[J].IEEE Transactions on Vehicular Technology, 2013, 62(5): 2031–2041.DOI:10.1109/TVT.2013.2238265
[18]	LAMBERT A, GRUYER D.A fast Monte Carlo algorithm for collision probability estimation[C]//10th International Conference on Control, Automation, Robotics and Vision.Piscataway, NJ:IEEE Press, 2008:406-411.
[19]	KARLIN S, TAYLOR H M. 随机过程初级教程[M]. 庄无兴, 陈宗洵, 陈庆华, 译. 北京: 人民邮电出版社, 2007: 306-351. KARLIN S, TAYLOR H M.A first course in stochastic processess[M].ZHUANG W X, CHEN Z X, CHEN Q H, translated.Beijing:Posts & Telecom Press, 2007:306-351(in Chinese).http://210.37.2.189/opac/item.php?marc_no=0000326501
[20]	HU J H, PRANDINI M, SASTRY S. Aircraft conflict prediction in the presence of a spatially correlated wind field[J].IEEE Transactions on Intelligent Transportation Systems, 2005, 6(3): 326–340.DOI:10.1109/TITS.2005.853699
[21]	BLIN K, AKIAN M, BONNANS F, et al.A stochastic conflict detection method integrating planned heading and velocity changes[C]//Proceedings of 39th IEEE Conference on Decision and Control.Piscataway, NJ:IEEE Press, 2000, 5:4717-4722.