一种曲线飞行航段及误差分布确立方法与流程

本发明涉及空中交通管理技术领域。

背景技术：

碰撞风险评估是民航保障安全底线的最重要工作之一，风险本质上是随机发生的，不可能根本的杜绝风险发生的可能性。为此国际民航组织对空中交通碰撞风险提出了安全目标水平的概念(每对飞机发生碰撞频率为5×10^-9次/飞行小时)，不同的学者也尝试着从各个方面度量碰撞风险，旨在合理的判断空中交通系统是否满足安全目标水平的要求、确定航空器间安全间隔等。碰撞风险评估的核心在于飞行误差的确立，通常涉及到纵向误差、侧向误差、竖向误差三个维度，现有方法侧重于宏观分析，对于某些微观场景精度不足。例如：

1、现有多采取航空器飞行轨迹直线模型评估碰撞风险，在航空器转弯的时候该模型将不再适用。

2、大量研究假设航空器的飞行误差在一段航段内满足固定的分布，虽然方便碰撞风险评估的计算，但是对风险评估精度造成一定的损失。

基于此，合理的航空器飞行轨迹曲线模型是提高碰撞风险评估精度的重大改进，特别是应用在航空器进离场阶段。同时，更为精确的，与航空器飞行轨迹中对应的空间位置相关的航空器飞行误差分布，在分析航空器爬升降落等微观场景时尤为重要。

技术实现要素：

本发明的目的是：针对现有航空器曲线飞行航段缺乏误差分析的技术缺陷，本发明利用大量样本数据通过统计分析技术拟合航段曲线、分析飞行误差，建立与航空器在航路中对应的空间位置匹配的航空器飞行误差分布，提高航空器碰撞评估的准确性，尤其适用于航空器曲线爬升阶段复杂情况或航班转弯情形。

本发明解决其技术问题所采用的技术方案是：

本发明的发明内容如下：

一种曲线飞行航段的确定方法，包括如下步骤：

s1建立航段曲线模型；

s2回归模拟航段；

s3估计航段侧向飞行误差分布；

s4估计航段垂向飞行误差分布。

优选的，所述s1建立航段曲线模型，包括：

航段曲线l(s)＝(x(s)，y(s))，包括长度为l0的直线、半径为r1圆心角为α1的圆弧、半径位r2圆心角位α2的圆弧、斜率k＝tan(α1-α2)的直线；

l(s)＝(x(s)，y(s))的数学表达式为：

其中，l0＝x0，l1＝x0+r1α1，l2＝x0+r1α1+r2α2，

x2＝x0+r1sinα1+r2(sinα1-sin(α1-α2))

y2＝y0+r1(1-cosα1)-r2(cosα1-cos(α1-α2))

l(s)由参数(x0，y0，r1，α1，r2，α2)唯一确定。

优选的，所述s2回归模拟航段，包括：

s201记di为样本点(xi，yi)到航路段l的距离，即：

其中inf代表实数集合的下确界

s202利用最小二乘法优化目标函数ls(x0，y0，r1，α1，r2，α2)，

得到的参数(x0，y0，r1，α1，r2，α2)为曲线l(s)的最佳拟合参数；

s203将最佳拟合参数带入式(1)、式(2)，获得l(s)数学表达式。

优选的，所述s3估计航段侧向飞行误差分布，包括：

s301取步长为δs，将航段曲线划分为子航段li＝l(s)，s∈[i×δs，(i+1)×δs]，假设每个子航段li内各点对应的侧向飞行误差εy(s)均满足标准差为σi的正态分布，即对任意的s∈[i×δs，(i+1)×δs]，εy(s)～n(0，σi)，其中σi为标准差；

s302航段曲线的法平面nplane(s)为：为竖直单位向量，为曲线l(s)的法向量；

利用nplane(i×δs)与nplane((i+1)×δs)之间的样本点，采用极大似然估计法，得到参数σi的参数估计值；

s303将参数σi的参数估计值带入εy(s)～n(0，σi)获得航段l(s)侧向飞行误差εy(s)的分布函数。

优选的，所述s4估计航段垂向飞行误差分布，包括：

s401对于指定s，航班的实际位置z(s)应该由其预计位置zintend(s)以及由航空器系统飞行误差组成；

假设航空器系统误差预计位置zintend(s)服从某一区间的均匀分布，即：zintend(s)～u[z1(s)，z2(s)]；取zintend(s)生成的中线为航段，那么由预计轨迹的不确定性产生的垂向飞行误差满足分布u[-z0(s)，z0(s)]，假设与系统飞行误差相互独立，那么随机向量的联合概率密度函数为

故εz(s)的概率密度函数p(x)为：

其中，函数erf(x)称为x的误差函数，定义为：

s402分别对每一段li，利用nplane(i×δs)与nplane((i+1)×δs)之间的样本点带入式(6)，样本数据采用极大似然估计法估计参数

s403将参数的参数估计值带式(6)获得航段l(s)垂向飞行误差εz(s)的概率密度函数p(x)。

本发明的有益效果：

1、本发明的曲线飞行航段确立方法适用于航空器曲线爬升阶段复杂情况或航班转弯的情形，是安全风险评估的重大改进。

2、本发明的曲线飞行航段确立方法克服现有技术在宏观上假设航空器的飞行误差总体满足固定的分布，对风险评估的精度有一定的损失，导致空域资源的利用得不到极大的经济化的缺陷。

3、本发明的方法采用的飞行误差分布与航空器在航路中对应的空间位置有关，在分析航空器爬升降落时是至关重要的。

附图说明

下面结合附图和实施例对本发明进一步说明。

图1为航段曲线模型示意图；

图2为航段碰撞风险区边界对应参数的示意图；

图3为航段曲线拟合效果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不限定本发明。

本发明所述航段均为民航航空器飞行航段，航段碰撞风险区是根据民航航空器飞行误差所确定的保障民航航空器和无人机碰撞概率小于指定目标安全水平e的区域。若无人机在该区域以外飞行，其与相应航段内飞行的民航客机发生碰撞的概率小于e，此时将被认为是可接受的安全风险。

本发明在假设获得了大量的样本数据的基础上，如何通过统计分析技术拟合航段曲线、分析飞行误差，进而判定碰撞风险区边界的技术。本发明处理的数据源来自于空中交通管制单位运行的空管自动化系统，其数据格式满足asterixcat62格式规范，假设基于离场程序dep，我们采集到足够多的样本点pi＝(xi，yi，zi)。

s1航段曲线建模

本发明重点针对起飞、降落阶段的航段曲线建模问题，通过前期大量的数据观察，不失一般性的，假设航段l(此刻考虑的是航段在水平面的投影)由直线、圆弧连续一阶导数的拼接而成。如图1所示，航段曲线可由如下几部分组成：

第一部分，长度为l0的直线，其终点为(x0，y0)；

第二部分，半径为r1，圆心角为α1的圆弧(拐弯弧度)；

第三部分，半径位r2，圆心角位α2的圆弧(矫正过度拐弯的弧度)

第四部分，斜率k＝tan(α1-α2)的直线。

航段曲线l(s)＝(x(s)，y(s))对应的模型的精确数学描述为：

其中，l0＝x0，l1＝x0+r1α1，l2＝x0+r1α1+r2α2，

x2＝x0+r1sinα1+r2(sinα1-sin(α1-α2))

y2＝y0+r1(1-cosα1)-r2(cosα1-cos(α1-α2))

航段的模型由参数x0，y0，r1，α1，r2，α2确定。该模型在参数特殊化的时候，涵盖了通常情况下的情况，当α2＝0时，对应的是拐弯精准的情况；当α1＝α2＝0时，对应的为直线飞行的情况。

s2航线段l回归拟合

记di表示为样本点(xi，yi)到航路段l的距离，即：

取目标优化目标函数：

得到的参数(x0，y0，r1，α1，r2，α2)为曲线的最佳拟合参数。故而通过前述分析，我们可以由实际采集的样本数据得到航段的具体形式l(s)。

s3侧向误差εy(s)分析

前述步骤已经通过采集的数据点拟合了航段l(s)，现在分析其飞行误差。取步长为δs，将航段l划分为子航段li＝l(s)，s∈[i×δs，(i+1)×δs]，假设每个子航段li内各点对应的侧向误差均满足标准差为σi的正态分布，即对任意的s∈[i×δs，(i+1)×δs]，有：

εy(s)～n(0，σi)

利用nplane(i×δs)与nplane((i+1)×δs)之间的样本点，采用极大似然估计法，得到参数σi的参数估计。

极大似然估计法原理就是固定样本观测值(x1，x2，...，xn)，挑选参数θ使l(x1，x2，...，xn；θ)＝maxl(x1，x2，...，xn；θ)，这样得到的θ与样本值有关，θ(x1，x2，...，xn)称为参数θ的极大似然估计值，其相应的统计量θ(x1，x2，...，xn)称为θ的极大似然估计量。利用lnl(θ)是l(θ)的增函数，故lnl(θ)与l(θ)在同一点处达到最大值，于是对似然函数l(θ)取对数，利用微分学知识转化为求解对数似然方程：

解此方程求得的驻点即为所求的最值点，就可以很到参数的极大似然估计值。

s4垂向误差εz(s)分析

由于起飞航空器在脱离跑道的位置并不相同，故假设航空器预计轨迹服从一定的均匀分布，对于指定s，航空器的实际位置z(s)应该由其预计位置zintend(s)以及由航空器系统飞行误差组成(航空器系统飞行误差是指航空器沿着预计轨迹运行时由各种原因产生的不可不免的实际位置与预计位置的偏差)，本发明假设航空器系统飞行误差zintend(s)服从某一区间的均匀分布，即：zintend(s)～u[z1(s)，z2(s)]。取zintend(s)生成的中线为航段，那么由预计轨迹的不确定性产生的误差满足分布u[-z0(s)，z0(s)]，假设与系统飞行误差相互独立，那么随机向量的联合概率密度函数为

故εz(s)的概率密度函数p(x)为：

其中，函数erf(x)称为x的误差函数，定义为：

按照前述方法，分别对每一段li，利用样本数据采用极大似然估计法估计参数进而得到εz(s)的分布函数。

s5碰撞风险区边界判定

如图2所示，航段s∈[0，s1]对应的碰撞风险区rz(s1)为一个包含s∈[0，s1]在内的三维区域，其边界由以下各个部分组成：

后侧竖直面bsf(s1)：即曲线l(s)在起点处的法平面

bsf(s1)＝nplane(0)

左侧竖直面lsf(s1)：位于曲线l左侧，且距l上每一点距离均相等的竖直面，即以曲线为底边的竖直面

右侧竖直面rsf(s1)：位于水平投影曲线l(s)右侧，且距l上每一点距离均相等的竖直面，即以曲线为底边的竖直面

前侧竖直面fsf(s1)：包含水平投影曲线l(s)在点l(s1)处法平面

fsf(s1)＝nplane(s0)

一个底面ssf(s1)：垂直于左侧向竖直面、右侧竖直面且与航迹曲线平行的曲面

碰撞风险区rz(s1)描述了由后侧竖直面bsf(s1)、左侧竖直面lsf(s1)、右侧竖直面rsf(s1)、前侧竖直面fsf(s1)、底面ssf(s1)围成的包含在内的三维区域，该区域是限制无人机运行的核心区域。

假设航空器实际运行航迹在点l(s)的侧向误差εy(s)，垂向误差εz(s)，fy(s，εy)是侧向误差εy(s)分布函数，fz(s，εz)是垂向误差εz(s)分布函数，e为可接受的碰撞概率，那么航段曲线对应的碰撞风险区为满足下述条件的最小区域：

条件一：假设该航段运行航空器飞机的最大侧宽为2λy，那么左、右侧竖直面参数dl、dr应满足：

其中p(x)是事件x发生的概率，fy(s，εy)是侧向误差分布函数，εy(s)是侧向误差；

条件二：假设该航段运行航空器的最大高度为2λz，那么底面参数应满足：

p(εz(s)＜-dz+λz)＝fz(s，-dz+λz)≤e

其中p(x)是事件x发生的概率，fz(s，εz)是垂向误差分布函数，εz(s)是垂向误差；

事实上满足上述两个条件的参数dl，dr，dz有无穷多个，理由如下，任意参数dl满足上述条件，那么可知dl+1依然满足条件；

本发明称满足以上两个条件的参数dl，dr，dz为e风险规避，将e风险规避的参数dl，dr，dz的集合分别表示为ral(e)，rar(e)，raz(e)，那么碰撞风险区边界对应的参数分别为：

dl＝infral(e)

dr＝infrar(e)

dz＝infraz(e)

其中：inf代表实数集合的下确界；

dl是航段曲线与左侧竖直面lsf(s1)之间的距离；

dz是航段曲线与底面ssf(s1)之间的距离；

dr是航段曲线与右侧竖直面rsf(s1)之间的距离。

实施案例一：

采集某机场运行数据，提取数据记录三维空间坐标位置(θ，φ，h)，其中(θ，φ)分别为对应的维度和经度，h为高度，单位为英尺。

该机场的基准点经纬度为(n294308，e1063829)，跑道对应的朝向为真向角17°，机场标高为h0＝1364英尺，该城市对应的地球半径大约为r＝6377830米，经计算机场基准点对应的经纬度的弧度制表示为：

θ0＝0.6186924

φ0＝1.8612433

那么可知沿着该跑道方向应作坐标变换：

x＝r(φ-φ0)sin(β)cos(θ0)+(θ-θ0)rcos(β)

y＝r(φ-φ0)cos(β)cos(θ0)-(θ-θ0)rsin(β)

根据微分几何以及微积分的知识可知，在机场周边50千米内，经度的误差δθ＜50000/r＜0.008，由泰勒展开

cosθ＝cos(θ0)-sin(θ0)δθ+o(δθ)²

故坐标变换在机场周边50km类的精度大于99％。

得到样本点pi的笛卡尔坐标(xi，yi)，利用样本(xi，yi)来确定模型参数(x0，y0，r1，α1，r2，α2)。样本点(xi，yi)与曲线距离di的计算过程如下：

最小化ls(x0，y0，r1，α1，r2，alpha2)得到参数结果为：

(3748.91714，2018.654，7830.656，0.3684，24117.5824，0.112529)

因此，得到的参数3748.91714，2018.654，7830.656，0.3684，24117.5824，0.112529)为曲线的最佳拟合参数，上述拟合效果如图3所示。

实施案例二：某机场进场程序下的碰撞风险区划设

部分跑道离场程序对应的航段为曲线(航班起飞的过程中有转弯的情况)。本例中对应的航段s∈[0，10000]为：

z(s)＝0.05s

该航段描述的是起飞2千米后沿着半径为6千米的圆弧转弯角度后继续直线飞行的情况。本例之中，依然假设航班侧向误差εy(s)与上述例子服从相同的分布，并假设εz(s)的概率密度函数p(x)为：

其中，函数erf(x)称为x的误差函数，定义为：

根据约束条件一、约束条件二，可取dl＝dr＝283.887，dz＝133.5。

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。尽管本发明已进行了一定程度的描述，明显地，在不脱离本发明的精神和范围的条件下，可进行各个条件的适当变化。可以理解，本发明不限于所述实施方案，而归于权利要求的范围，其包括所述每个因素的等同替换。