CN112116830B - 一种基于空域网格化的无人机动态地理围栏规划方法 - Google Patents

Abstract

本发明公开了一种基于低空空域网格化的无人机动态地理围栏规划方法，包括考虑无人机实际航迹随机性特性，对无人机飞行计划进行柔性航迹预测；结合无人机柔性航迹预测结果，进行动态地理围栏纵向长度尺寸设计；对目标空域范围内的历史航迹数据进行预处理，并对空域进行均匀网格化处理与编号；定义影响低空空域复杂性的关键指标，采用核K‑均值聚类方法对空域网格进行聚类分析实现低空空域复杂性分级度量结果；结合低空空域网格复杂度，进行动态地理围栏动态地理围栏横向范围规划。本方法首次提出一套完整的无人机动态地理围栏规划方法并将其运用到低空空域空中交通管理中，有效提高低空空域的利用率，有利于缓解低空空域资源紧张问题。

Description

一种基于空域网格化的无人机动态地理围栏规划方法

技术领域

本发明涉及空中交通规划管理技术，特别是涉及一种基于空域网格化的无人机动态地理围栏规划方法。

背景技术

随着民用无人机应用的不断深入和无人机市场的快速扩大，未来低空空域的运行管理将面临巨大的挑战。低空空域资源紧缺和分配不当将带来许多问题，一方面，关于低空空域的空中交通管理方式，采用传统管理模式下划分隔离空域方式或规划静态地理围栏方式，都是较为粗放型的管理模式，在部分因经济或地理因素所导致的较繁忙空域资源将无法满足需求。另一方面，由于低空空域的航空器种类、机型大小以及飞行方式各种各样，导致低空空域的管理模式无法照搬运输航空空域的扇区划分管制模式。因此关于低空空域的空中交通管理模式应该结合其自身特点，针对空域资源紧缺，无人机之间潜在冲突等问题进行合理规划。

发明内容

发明目的：本发明的目的是提供一种基于空域网格化的无人机动态地理围栏规划方法，解决现有的划设隔离空域或静态地理围栏规划等低空空中交通管理方式所带来的低空空域资源紧张以及利用率低下等实际问题。

技术方案：本发明的基于低空空域网格化的无人机动态地理围栏规划方法，包括以下步骤：

S1、考虑无人机实际航迹随机性特性，对无人机飞行计划进行柔性航迹预测；

S2、结合无人机柔性航迹预测结果，进行动态地理围栏纵向长度尺寸设计；

S3、对目标空域范围内的历史航迹数据进行预处理，并对空域进行均匀网格化处理与编号；

S4、定义影响低空空域复杂性的关键指标，采用核K-均值聚类方法对空域网格进行聚类分析实现低空空域复杂性分级度量结果；

S5、结合低空空域网格复杂度，进行动态地理围栏横向范围规划。

进一步的，步骤S1具体为：

首先考虑无人机性能及飞行计划数据，综合推算出水平轨迹、高度剖面以及速度剖面，从而获得无人机初始4D航迹；在此基础上，考虑无人机实际航迹随机特性，即风速的不确定性和指示空速的不稳定性导致实际的平均地速与计算得到的结果产生的差异；

假设无人机在巡航阶段的指示空速V_TAS和风速V_WS均服从正态分布，其中指示空速

风速

其概率密度函数分别为：

结合正态分布的特点和平均地速的表达式可知平均地速V_GS服从正态分布，即

且其概率密度函数为：

由此推算出无人机在某航段上的平均空速实际取值在μ±3σ的范围内，即

进一步的，步骤S2具体为：

在飞行习惯和风速的综合影响下，无人机在t时刻的平均地速在[μ_GS-3σ_GS,μ_GS+3σ_GS]范围内，则对应的动态地理围栏的纵向长度

为：

因此无人机动态地理围栏的纵向长度为6tσ_GS，同时会随着规划时间步长增加而增加，符合长期规划时不确定性增加的基本事实。

进一步的，步骤S3具体为：

首先选取目标空域内的历史航迹数据，然后对历史航迹数据进行预处理，剔除异常数据、非目标空域范围内的数据和非低空空域内的航迹点数据；然后，依据墨卡托投影将目标空域转为平面直角坐标系，然后在平面坐标系上对目标空域进行均匀离散化，得到二维平面上的均匀空域网格并进行编号；然后对广播式自动相关监视数据(ADS-B：AutomaticDependent Surveillance-Broadcast)格式的历史航迹数据进行解析处理，获得每个空域网格的原始指标数据；并进行编号；以坐标原点为起点，主方向为以x轴方向顺序递增，次方向为以y轴方向递增；同时，已知每个网格空域的单位长度，取每个网格的左下角坐标表示其网格角标。

进一步的，步骤S4具体为：

S41、在空域网格化的基础上引入航迹历史数据，将航迹密度，航迹点速度分布，航迹点高度分布，机型混合程度，潜在冲突程度，超障高度以及航迹时间密度几个方面作为影响低空空域复杂性的关键指标；采用线性加权法对关键指标进行处理，得到网格空域复杂度的综合评价值，即为复杂度综合评价模型；线性加权法为：

S42、结合单位网格多个复杂性影响指标的高维特征和综合评价结果，引入核函数理论并结合K-均值聚类算法对空域网格复杂性进行分级量化研究。

首先对多指标综合量化结果进行处理，根据实际问题确定类簇个数K，并根据样本总数m将样本均分为K份，标记为C＝{c₁,c₂,…,c_K}，取每个样本中综合量化值中位数作为初始聚类中心；然后将原始数据经过非线性变换，映射到高维特征空间，再在特征空间中进行K-均值聚类；

已知选取的空域网格样本总数为m，则样本集为{g⁽¹⁾,g⁽²⁾,…,g^(m)}，其中g⁽ⁱ⁾∈G^N，i＝1,2,3,...m，依据Mercer定理，若已知映射Φ:

使得K(g⁽ⁱ⁾,g^(j))＝Φ(g⁽ⁱ⁾)^TΦ(g^(j))，K(g⁽ⁱ⁾,g^(j))则能够用来表示核函数，其中G^N是选取的空域网格样本g⁽ⁱ⁾所在的空间，G^F是选取的空域网格样本g⁽ⁱ⁾经过映射Φ后得到的Φ(g⁽ⁱ⁾)所在的空间；利用核函数改变样本的分布形式，而不必知道其在映射高维空间中的具体表现形式；核K-均值聚类就是讨论原始样本集{g⁽¹⁾,g⁽²⁾,…,g^(m)}经过映射Φ后得到的数据集{Φ(g⁽¹⁾),Φ(g⁽²⁾),…,Φ(g^(m))}在G^F空间中的聚类情况；样本空间中第i个样本到第k类中心的距离为：

以综合量化值对样本进行顺序排序，按照设置的复杂性等级数量、各等级样本比例，初步划分样本的复杂性等级；然后将样本映射到高维核空间，取每个复杂性等级样本的综合量化值的均值作为初始聚类中心。

进一步的，步骤S5具体为：

S51、求解目标预测航迹点；

基于S1中提出的无人机柔性航迹预测技术，预测无人机目标的航迹点；

S52、求解基本网格空域集合；

两个航迹点的坐标P_start(x₀,y₀)和P_end(x₁,y₁)，根据直线方程求出两点连线的表达式y＝rx+b；首先通过计算Δx＝|x₀-x₁|和Δy＝|y₀-y₁|来确定移动的主方向：

如果Δx>Δy，说明x轴的最大差值大于y轴的最大差值，x轴方向为移动的主方向；如果Δx<Δy，说明y轴的最大差值大于x轴的最大差值，y轴方向为移动的主方向；将起始点坐标分别对网格边长向下取整，即令

其中a是网格边长；

当x轴为移动的主方向时，设当前点为P′_start(x′₀,y′₀)，其中x′₀＝ai，y′₀＝aj，当x＝x′₀+a时，根据直线公式算出y＝k(x′₀+a)+b，对y以网格边长向下取整得，

令σ＝y′-y′₀，则：

计算

步后，当x＝x′₁，即终点所在网格时，求出起点和终点连线所在的外包网格空域集合G_base，当y轴为移动的主方向时同理；

S53、滑动邻域法求解扩展空域；

设S为初始解，即基本空域集合，邻域的个数为M，N_l表示第l个邻域。在进行邻域搜索时，只要当前邻域搜索结束后优化了S，即表现为效益函数增大时，则将当前邻域加入最优解，并进行下一个邻域结构进行运算，否则停止运算，输出当前解；

根据取整法得到航迹点连线穿过的网格空域集合G_base，以G_base作为初始解S；将初始解向两侧依次平移M次，得到M个递增的邻域，作为邻域结构集，分别记为N₁,N₂,…,N_M；

邻域结构间的搜索顺序可以通过改变邻域结构之间的次序，一般来说会按照邻域的大小顺序进行排序，从小到大，即：N₁(x)≤N₂(x)≤N₃(x)≤…≤N_max(x)；根据实际情况，邻域结构之间的次序为以初始解为基准向两侧依次交叉排列。

S54、定义效益函数F_obj表示某规划阶段所分配的网格空域的最大效益，其中该空域集合由基本空域和扩展空域两部分组成，分别为航迹中心线所在网格空域和与其相邻的部分网格空域；用f_dis表示增加当前扩展空域带来的自由度效益，f_c表示该网格空域产生的复杂性负效益，g_ij为决策变量，表示当增加某网格带来的总效益为正时，则将其划分为扩展空域，效益函数表达式为：

其中，由于目标为航线运行规划，考虑到距离航迹点连线过于远的网格对于地理围栏的规划没有实用意义，因此趋向于在预测航迹附近寻找复杂性等级较低的连续网格空域集合，即采用当前空域网格与预测航迹所在直线的距离的倒数来表示自由度效益：

其中，μ为自由度可调参数，d为当前空域网格与预测航迹所在直线的距离，s为复杂度值，λ用来保证规划空域网格的时候，距离航线太远的网格即使它复杂度很低也不把它划进规划队列，λ的值根据实际情况设置，σ根据实际需求用来调节复杂度带来负效益的影响，σ的值根据实际情况设置；

由于自由度效益函数为双曲线中的第一象限中的单支，因此随着距离d的增大，自由度效益f_dis会逐渐减小，设置μ为自由度可调参数，当d大于μ时，自由度效益小于1，同时，由于复杂度值由经过标准化后的指标数据加权平均值表示，考虑所有指标对复杂性的影响程度相同，则复杂度值的取值范围为s∈[0,1]，令λ＝1.5μ,σ＝1，则f_c∈[1,1.5μ+1]，当d大于μ时，效益函数F_obj≤0；

模型约束条件为：

决策变量为0-1变量，表示当前网格是否在最优解集中，当g_ij为1时，表示最优解集包含当前网格，当g_ij为0时，则不包括当前网格；

g_ij＝0,1；

作为无人机完成飞行任务的前提条件，必须有一个基本的网格空域满足当前飞行任务的需要，在模型中则要求起始点和终点连线所在网格必须为最优解集中；

g_ij＝1(若g_ij∈G_base)；

可行域范围约束，若起始点和终点所在网格分别为g_i0j0和g_i1j1，则当起点和终点所在直线斜率r<0时：

斜率r>0时：

斜率r＝0时：

i₀≤i≤i₁,i,j∈N；

斜率r不存在时：

j₀≤j≤j₁,i,j∈N；

连通性约束，已知起点和终点连线的直线方程为Ax+By+C＝0，其中A,B不同时为0，则当Ai+Bj+C≤0时，即网格位于连线下方时：

当Ai+Bj+C≥0时，即网格位于连线下方时：

S55、利用S4提出的空域复杂分级度量方法和S54提出的效益函数，计算基本空域网格的效益函数是否为正，若为正，则表示将基本空域网格集合向两侧依次滑动，若为负，则输出所有扩展空域网格集合。

有益效果：与现有技术相比，本发明首次提出一套完整的基于空域网格化的无人机动态地理围栏规划方法，弥补了规划静态地理围栏以及划设隔离空域等空域资源分配方式的不足，通过时隙分配和低空空域复杂性自适应地调整地理围栏大小，使得空域复杂性较高的繁忙低空空域资源得到更加细致、精准和高效地利用，同时利用复杂性较低的空闲低空空域增加无人机飞行活动的空间自由度以及安全性。因此，对无人机进行动态地理围栏的规划在一定程度上能实现合理分配低空空域资源的目的，对推进低空空域的精细化管理模式具有十分重要的意义。

附图说明

图1是无人机动态地理围栏规划流程示意图；

图2是本发名明中动态地理围栏纵向尺寸设计流程示意图；

图3是低空空域复杂度分级度量模型示意图；

图4是本发明中动态地理围栏横向范围规划流程示意图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

图1为本发明的流程图，具体可包括以下步骤：

S1、考虑无人机实际航迹随机性特性，对无人机飞行计划进行柔性航迹预测；

对于目标飞行计划，首先结合无人机性能参数建立了无人机初始静态4D航迹模型，然后分析无人机在实际飞行中的航迹随机性特征，根据航迹随机性特征对无人机初始4D航迹模型进行改进，最后生成柔性4D航迹预测模型。

无人机初始静态4D航迹预测模型具体为：4D航迹是指基于四维从当前位置到结束位置连续的对于航空器飞行轨迹的精确描述，在4D航迹上的每个点都精确的与时间关联。运输航空4D航迹研究是为了提高空中交通安全和空中交通效率，主要是以航空器航迹管理为中心，通过准确安排航空器运行航迹和时间，提高空域和机场资源的利用率和安全性；4D航迹预测模型主要包括a)水平轨迹生成；b)飞行高度剖面；c)飞行速度剖面和d)初始航迹生成。

a)水平轨迹生成：无人机的水平航迹是指无人机在整个飞行任务中的航迹在水平面的投影。由于固定翼无人机的飞行原理可知，无人机在调整航向时无法悬停转弯，需要进行水平转弯机动，因此无人机实际航迹是由曲线与直线段共同构成。当无人机的航程较长，直线段的长度远远大于圆弧时，可以适当忽略转弯时的曲线，则水平航迹为纯直线段组成。当航迹的精确度要求较高时，需要考虑无人机在转弯时的航迹曲线，

b)飞行高度剖面：无人机的飞行高度剖面是指空间直角坐标系中，无人机从起飞至降落的航迹在铅锤面上的投影。爬升阶段的飞行高度剖面可以按照无人机的爬升率或爬升梯度可以算出以时间间隔或航程间隔下的各个航路点的飞行高度；巡航阶段飞行高度则保持不变；下降阶段的飞行高度剖面以巡航阶段高度为起始高度，根据下降率和下降梯度计算得出。

c)飞行速度剖面：无人机的飞行速度剖面是指无人机速度随航程变化的规律。由于本文的研究对象为在低空空域进行飞行作业的固定翼无人机，因此在飞行高度变化不大的前提下，忽略飞行高度的影响，近似认为指示空速等于真空速。典型的无人机飞行速度剖面一般分为三个阶段，即加速段、匀速段以及减速段，在无人机达到巡航指示空速之前需要在其自身性能允许的情况下以一定的加速度加速飞行，然后保持匀速进行平飞直到需要进行减速时。

d)初始航迹生成：根据水平轨迹生成得到每个航路点的地理坐标，再根据速度剖面信息、航路点之间的距离和风场信息，即可计算从当前航路点到下一个航路点所需要的时间，从第一个航路点算起，即可得到无人机经过每一个航路点的时间，结合无人机飞行高度剖面，即可实现无人机初始静态4D航迹预测

无人机柔性4D航迹预测模型为：实际飞行过程中，由于既无法保证无人机严格按照最大巡航速度飞行，同时由于风对无人机会产生一定的干扰，导致实际飞行与理想情况存在一些误差。因此综合考虑平均地速的正态分布特征，根据“3σ原则”，得到无人机在每一航段内平均地速的取值范围，进而得到一个更合理的无人机当前时刻所在的位置范围。

如图2所示，首先考虑无人机性能及飞行计划数据，综合推算出水平轨迹、高度剖面以及速度剖面，从而获得无人机初始4D航迹。在此基础上，考虑无人机实际航迹随机特性，即风速的不确定性和指示空速的不稳定性导致实际的平均地速与计算得到的结果产生的差异。

假设无人机在巡航阶段的指示空速V_TAS和风速V_WS均服从正态分布，其中指示空速

风速

其概率密度函数分别为：

其中，μ_TAS表示平均指示空速，σ_TAS表示指示空速标准差，μ_WS表示风速平均值，σ_WS表示风速标准差。结合正态分布的特点和平均地速的表达式可知平均地速V_GS服从正态分布，即

且其概率密度函数为：

由此可推算出无人机在某航段上的平均空速实际取值在μ±3σ的范围内，即

其中DA为偏流角，表示航迹线与航向线的偏离程度，WA为风角，表示风向线同航迹线之间的夹角。

S2、结合无人机柔性航迹预测结果，进行动态地理围栏纵向长度尺寸设计；

如图2所示，综合考虑平均地速的正态分布特征，从而对动态地理围栏的纵向长度计算方法进行优化调整，得到一个更合理的无人机当前时刻所在的位置范围，根据柔性航迹预测结果，已知无人机在某时刻的平均地速范围和航段平均飞行时间，可以得到对应的动态地理围栏的纵向长度。在飞行习惯和风速的综合影响下，无人机在t时刻的平均地速在[μ_GS-3σ_GS,μ_GS+3σ_GS]范围内，则对应的动态地理围栏的纵向长度

为：

因此无人机动态地理围栏的纵向长度为6tσ_GS，同时会随着规划时间步长增加而增加，符合长期规划时不确定性增加的基本事实。其中，μ_GS表示平均地速，σ_GS表示地速的标准差。

选用某型的固定翼无人机在货物运输活动过程中进行案例仿真，选取某个已处理的飞行计划为具体分析对象，其信息包括计划编号P-20191201-1、计划时刻12:00、航路点1、2，其经纬度坐标分别为(31.5,119)、(31.8,119.2)。现对两航路点的航段进行分析，依据初始航迹生成方法得到无人机在该航段上的预测航迹及航段飞行时间。

当无人机在该航段上单向飞行时，假设已知当前航段上的指示空速和风速的分布律，根据蒙特卡洛仿真可以得到平均地速的分布参数，同时根据两个航路点的经纬度可以计算出航程，从而计算该航段的平均飞行时间为0.322h。无人机的地速服从正态分布，均值为119.792，方差为16.752，最大值为135.034，最小值为101.621。

定义动态地理围栏纵长L_CD，即无人机围栏在航向上的长度，表示为动态地理围栏的最前端和目的航路点两者之间的较小值与动态地理围栏后端和起始航路点之间的最大值之差。取规划步长τ＝5s，计算出采用速度范围的方法时该航段上每个规划时刻的L_CD，即以无人机性能的速度范围计算围栏的纵向长度，然后考虑航迹随机性的情况下的动态地理围栏纵长

对比结果如表1所示。

表1

从表中可以大致看出

小于L_CD，即采用的动态地理围栏纵向长度规划方法更加节省空域资源，在实际规划中，可以按需合理调整规划步长τ，尤其是在规划初期

较短，可适当将步长增加以保证无人机安全飞行。

S3、对目标空域范围内的历史航迹数据进行预处理，并对空域进行均匀网格化处理与编号；

如图3所示，首先选取目标空域内的历史航迹数据，然后对历史航迹数据进行预处理，剔除异常数据、非目标空域范围内的数据和非低空空域内的航迹点数据。然后，依据墨卡托投影将目标低空空域地理坐标系转为平面直角坐标系，然后在平面直角坐标系上以一定的划分级别对目标低空空域进行均匀离散化，得到二维平面上的均匀空域网格，然后对广播式自动相关监视数据(ADS-B：Automatic Dependent Surveillance-Broadcast)格式的历史航迹数据进行解析处理，获得每个空域网格的原始指标数据；并进行编号。以坐标原点为起点，主方向为以x轴方向顺序递增，次方向为以y轴方向递增。同时，已知每个网格空域的单位长度，取每个网格的左下角坐标表示其网格角标。

S4、定义影响低空空域复杂性的关键指标，采用核K-均值聚类方法对空域网格进行聚类分析实现低空空域复杂性分级度量结果；

如图3所示，包括以下步骤：

S41、在空域网格化的基础上引入航迹历史数据，将航迹密度，航迹点速度分布，航迹点高度分布，机型混合程度，潜在冲突程度，超障高度以及航迹时间密度几个方面作为影响低空空域复杂性的关键指标。采用线性加权法对关键指标进行处理，得到网格空域复杂度的综合评价值，即为复杂度综合评价模型。线性加权法为：

其中，CR_i是第i个网格空域复杂度的综合评价值，ω_j是第j个指标的权重，F′_j是第j个指标的评价值，n表示选取的影响低空空域复杂性的关键指标的个数。

根据指标归一化后的结果，可以计算第i个网格空域的复杂性相对值CR_i，设向量中有n个指标，当所有指标被认为有同等权重时，则ω_j都取

否则ω_j可根据当前指标的重要程度取不同的值。

S42、结合单位网格多个复杂性影响指标的高维特征和综合评价结果，引入核函数理论并结合K-均值聚类算法对空域网格复杂性进行分级量化研究。

首先对多指标综合量化结果进行处理，根据实际问题确定类簇个数K，并根据样本总数m将样本均分为K份，标记为C＝{c₁,c₂,…,c_K}，取每个样本中综合量化值中位数作为初始聚类中心。以此来保证聚类结果有较好的综合量化值分级特征。然后将原始数据经过非线性变换，映射到高维特征空间，再在特征空间中进行K-均值聚类。

已知选取的空域网格样本总数为m，则样本集为{g⁽¹⁾,g⁽²⁾,…,g^(m)}。其中g⁽ⁱ⁾∈G^N，i＝1,2,3,...m。依据Mercer定理，若已知映射Φ:

使得K(g⁽ⁱ⁾,g^(j))＝Φ(g⁽ⁱ⁾)^TΦ(g^(j))，K(g⁽ⁱ⁾,g^(j))则可以用来表示核函数，其中G^N是选取的空域网格样本g⁽ⁱ⁾所在的空间，G^F是选取的空域网格样本g⁽ⁱ⁾经过映射Φ后得到的Φ(g⁽ⁱ⁾)所在的空间。利用核函数改变样本的分布形式，而不必知道其在映射高维空间中的具体表现形式。核K-均值聚类就是讨论原始样本集{g⁽¹⁾,g⁽²⁾,…,g^(m)}经过映射Φ后得到的数据集{Φ(g⁽¹⁾),Φ(g⁽²⁾),…,Φ(g^(m))}在G^F空间中的聚类情况。样本空间中第i个样本到第k类中心的距离为：

以综合量化值对样本进行顺序排序，按照设置的复杂性等级数量、各等级样本比例，初步划分样本的复杂性等级。然后将样本映射到高维核空间，取每个复杂性等级样本的综合量化值的均值作为初始聚类中心。根据此方法得到的聚类结果会使每个类簇的综合量化值在一定的范围内，并且每个类簇的综合量化值基本呈现梯度分布，只需对其交叉重复部分和离群点进行处理即可得到较好的复杂性分级效果，实现了复杂度分级度量模型的建立。该结果兼顾了综合量化值和原始指标数据的高维结构特征，具有较好的实际意义。

选取某机场2019年7月23日的ADS-B数据进行案例仿真，首先对数据进行预处理。选取机场周围50km低空空域范围，纬度大约在34.9506292053至35.8544246137，经度在118.789668151至119.870509340的空域，然后将以WGS_84标准的经纬度表示的目标空域(同高度层，忽略同高度层上下界带的偏差)的地理坐标系转换为墨卡托投影平面直角坐标系，为了便于计算，将坐标原点平移到O′，即将令x′＝x-O(x)，y′＝y-O(y)，得到以目标空域范围左下角坐标为原点的平面直角坐标系，其中，x、y是原始平面直角坐标系中点在x轴和y轴上的坐标，x′、y′是平移后平面直角坐标系中点在x轴和y轴上的坐标，O(x)、O(y)分别是沿x轴和y轴上平移的长度。将直角坐标系下的目标空域范围均匀划分为30×30的单位网格空域，并对每个空域网格进行编号，然后对该机场三个自然日的ADS-B历史航迹数据进行处理，剔除异常数据、非目标空域范围内的数据和非低空空域内的航迹点数据，并根据位置信息将航迹点与网格一一对应。经过初步处理得到900个空域网格的航迹点数据，然后逐步对每个网格的原始数据进行综合量化值计算，得到900组8个复杂性指标的数据列表，部分结果如表2所示。

表2

为了呈现目标空域的复杂性等级分层情况，考虑对得到的标准化指标数据进行聚类分析，不考虑指标的综合量化结果，首先仅采用传统的K-means聚类方法进行分析。使用Python进行K-means方法聚类分析，以每个网格空域复杂性的影响指标间的欧氏距离为聚类依据，结果得到6个非均匀的类簇。将聚类结果在实际目标空域网格中对应后，得到该机场附近低空空域复杂程度的分布。但是这种聚类方法只能体现以多指标蕴含的复杂性相似的样本特征，如果需要对不同类之间进行复杂性等级划分，则需要进一步分析作为聚类中心的样本原始指标数据。为了弥补采用综合评价量化值表达复杂性的单一性缺陷，以及传统聚类分析法无法较好地进行量化分级的不足，采用核K-均值聚类方法对指标数据进行处理，使聚类结果包含综合评价量化结果。首先对样本进行初步处理，将900个样本按照指标综合评价值按从小到大的次序排序，将排序后的样本均匀分为6组，取每组样本的综合评价指标中位数所属的网格样本点作为初始聚类中心，结果为表3所示。

表3

然后选择高斯核函数表示两个样本点的高维空间距离，最后求出聚类分析结果。经分析，采用核K-均值聚类方法得到的结果在复杂性分布上产生了较好的分级效果。当然复杂性度在各类之间还是存在小部分重叠，这是由于聚类算法的本质是优先考虑高维核空间的距离，尽管这些样本点的综合评价指标比较接近，但其原始指标数据的高维空间结构却存在一定的差异，这个结果也从侧面说明了单纯考虑多指标的加权平均值不能全面地反映网格空域的复杂程度，因此，采用综合量化值与高维空间聚类相结合的方法对复杂性的描述更加科学和完善。

S5、结合低空空域网格复杂度，进行动态地理围栏动态地理围栏横向范围规划；

根据动态地理围栏纵向长度规划结果，以当前航迹点和下一航迹点为规划基础，将两个航迹点连线经过的所有单位网格空域标记为基本空域，结合空域复杂性度量结果，将基本空域相邻的若干空域复杂性较低的单位网格空域合理标记为扩展空域，与基本空域共同构成该时段无人机动态地理围栏的横向空域范围。

如图4所示，所述的动态地理围栏横向范围规划方法，具体包括以下步骤：

S51、求解目标预测航迹点；

基于S1中提出的无人机柔性航迹预测技术，预测无人机目标的航迹点；

S52、求解基本网格空域集合；

无人机的起点和终点坐标分别为P_start(x₀,y₀)和P_end(x₁,y₁)，根据直线方程可以求出两点连线的表达式y＝rx+b。首先通过计算Δx＝|x₀-x₁|和Δy＝|y₀-y₁|来确定移动的主方向：

如果Δx>Δy，说明x轴的最大差值大于y轴的最大差值，x轴方向为移动的主方向；如果Δx<Δy，说明y轴的最大差值大于x轴的最大差值，y轴方向为移动的主方向。将起始点坐标分别对网格边长向下取整，即令

其中a是网格边长，则P_start(x₀,y₀)位于g_ij所代表的空域网格。

当x轴为移动的主方向时，设当前点为P′_start(x′₀,y′₀)，其中x′₀＝ai，y′₀＝aj，当x＝x′₀+a时，根据直线公式可以算出y＝k(x′₀+a)+b，对y以网格边长向下取整可得，

令σ＝y′-y′₀，则：

计算

步后，当x＝x′₁，即终点所在网格时，可求出起点和终点连线所在的外包网格空域集合G_base，当y轴为移动的主方向时同理。

S53、滑动邻域法求解扩展空域；

设S为初始解，即基本空域集合，邻域的个数为M，N₁表示第l个邻域。在进行邻域搜索时，只要当前邻域搜索结束后优化了S，即表现为效益函数增大时，则将当前邻域加入最优解，并进行下一个邻域结构进行运算，否则停止运算，输出当前解。

根据取整法得到航迹点连线穿过的网格空域集合G_base，以G_base作为初始解S。将初始解向两侧依次平移M次，得到M个递增的邻域，作为邻域结构集，分别记为N₁,N₂,…,N_M；

邻域结构间的搜索顺序可以通过改变邻域结构之间的次序，一般来说会按照邻域的大小顺序进行排序，从小到大，即：N₁(x)≤N₂(x)≤N₃(x)≤…≤N_max(x)。根据实际情况，邻域结构之间的次序为以初始解为基准向两侧依次交叉排列。停止准则的选择对算法全局收敛性和时效性都有直接并且重要的影响。

S54、定义效益函数F_obj表示某规划阶段所分配的网格空域的最大总效益，其中该空域集合由基本空域和扩展空域两部分组成，分别为航迹中心线所在网格空域和与其相邻的部分网格空域。用f_dis表示增加当前扩展空域带来的自由度效益，f_c表示该网格空域产生的复杂性负效益，g_ij为决策变量，表示当增加某网格带来的总效益为正时，则将其划分为扩展空域，效益函数表达式为：

其中，N为x轴方向和y轴方向可行域网格的数量，由于目标为航线运行规划，考虑到距离航迹点连线过于远的网格对于地理围栏的规划没有实用意义，因此趋向于在预测航迹附近寻找复杂性等级较低的连续网格空域集合，即采用当前空域网格与预测航迹所在直线的距离的倒数来表示自由度效益：

其中，μ为自由度可调参数，d为当前空域网格与预测航迹所在直线的距离，s为复杂度值，λ用来保证规划空域网格的时候，距离航线太远的网格即使它复杂度很低也不把它划进规划队列，λ的值根据实际情况设置，σ用来根据实际需求用来调节复杂度带来负效益的影响，σ的值根据实际情况设置。

由于自由度效益函数为双曲线中的第一象限中的单支，因此随着距离d的增大，自由度效益f_dis会逐渐减小，设置μ为自由度可调参数，当d大于μ时，自由度效益小于1，其意义为距离预测航线过于远的空域对其飞行没有实际帮助，即自由度效益对总效益值贡献趋于0。同时，由于复杂度值由经过标准化后的指标数据加权平均值表示，考虑所有指标对复杂性的影响程度相同，则复杂度值的取值范围为s∈[0,1]，令λ＝1.5μ,σ＝1，则f_c∈[1,1.5μ+1]，当d大于μ时，效益函数F_obj≤0。

同时，效益函数的实际意义也表示随着扩展空域网格与预测航迹的距离逐渐增大，对其复杂性的要求也越来越低，而预测航线附近的网格空域可适当放宽复杂性等级要求，以保证为其分配的空域能满足安全飞行的要求。结合实际情况，效益函数模型约束条件为：

决策变量为0-1变量，表示当前网格是否在最优解集中，当g_ij为1时，表示最优解集包含当前网格，当g_ij为0时，则不包括当前网格。

g_ij＝0,1；

作为无人机完成飞行任务的前提条件，必须有一个基本的网格空域满足当前飞行任务的需要，在效益函数模型中则要求起始点和终点连线所在网格必须属于最优解集。

g_ij＝1(若g_ij∈G_base)；

其中，G_base是起点和终点连线所在的外包网格空域集合。

可行域范围约束，若起始点和终点分别位于g_i0j0和g_i1j1代表的空域网格，i₀和j₀表示起始点所在网格编号的下标，i₁和j₁表示终点所在网格编号的下标，i和j表示可行域网格编号的下标，则当起点和终点所在直线

斜率r<0时：

斜率r>0时：

斜率r＝0时：

i₀≤i≤i₁,i,j∈N；

斜率r不存在时：

j₀≤j≤j₁,i,j∈N；

连通性约束，已知起点和终点连线的直线方程为Ax+By+C＝0(A,B不同时为0),则当Ai+Bj+C≤0时，即网格位于连线下方时：

当Ai+Bj+C≥0时，即网格位于连线下方时：

其中，g_i,j+1是下标为(i,j+1)的空域网格对应的决策变量，g_i+1,j是下标为(i+1,j)的空域网格对应的决策变量，g_i-1,j是下标为(i-1,j)的空域网格对应的决策变量，g_i,j-1是下标为(i,j-1)的空域网格对应的决策变量，

S55、利用S4提出的空域复杂分级度量方法和S54提出的效益函数，计算基本空域网格的效益函数是否为正，若为正，则表示将基本空域网格集合向两侧依次滑动，若为负，则输出所有扩展空域网格集合。

选取某机场方圆50km范围的低空空域作为规划目标，根据坐标转换方法将两个航迹点(35.6078399956,119.4278243855)与(35.6322514834,119.5232603629)的地理转化为直角坐标，结果为P_start(71039.227,89620.1935)和P_end(81663.1114,92963.1438)，取单位网格空域长度为

对a向下取整得到其所在的网格编号为g_17,22和g_20,23，然后计算得到基本空域网格集合G_base。首先根据两个点的坐标得到两点间的直线方程为f(x)＝0.3146x+67266.7306，由于|k|<1，因此可确定x为主移动方向，即只需计算20-17+1＝4步即可得到G_base。

当m＝0时，

当m＝1时，将起点坐标经过规整得到P′_start(68000,88000)，令x增加一个网格单位长度，得到x₁＝72000，将x₁代入直线方程得到y₁＝89922.5138，向下取整得到y₁′＝88000，计算σ_0,1＝y′₁-y′₀＝0，根据基本网格空域集合求解公式可知，

以此类推，可以算出当m＝4时，

然后根据变邻域法求解扩展空域网格集合G_exp，本例中，主移动方向为x轴，且k>0，可以算出上邻域

下邻域

然后根据每个网格的复杂度值和到预测航迹的距离，计算每个邻域的效益函数F_obj，得到最终分配的空域网格集合为

具体网格编码结果如表4所示。

表4

本发明的一种基于低空空域网格复杂性的无人机动态地理围栏规划方法。该方法在低空空域网格化的基础上，建立了以航迹密度、速度波动、高度变化等关键影响指标的低空空域复杂性度量模型，结合无人机性能参数，考虑无人机在实际飞行中的航迹随机性特征，生成柔性4D航迹预测模型。通过柔性4D航迹预测结果规划动态地理围栏的纵向长度，结合空域网格复杂性的度量结果，以效益函数最大化为目标，确定动态地理围栏的横向范围。本方法首次提出一套完整的无人机动态地理围栏规划方法并将其运用到低空空域空中交通管理中，有效提高低空空域的利用率，有利于缓解低空空域资源紧张问题。

Claims (5)

1.一种基于低空空域网格化的无人机动态地理围栏规划方法，其特征在于，包括以下步骤：

S1、考虑无人机实际航迹随机性特性，对无人机飞行计划进行柔性航迹预测；

S2、结合无人机柔性航迹预测结果，进行动态地理围栏纵向长度尺寸设计；

S3、对目标空域范围内的历史航迹数据进行预处理，并对空域进行均匀网格化处理与编号；

S4、定义影响低空空域复杂性的关键指标，采用核K-均值聚类方法对空域网格进行聚类分析实现低空空域复杂性分级度量结果；

S5、结合低空空域网格复杂度，进行动态地理围栏横向范围规划；

具体为：

S51、求解目标预测航迹点；

基于S1中提出的无人机柔性航迹预测技术，预测无人机目标的航迹点；

S52、求解基本网格空域集合；

两个航迹点的坐标P_start(x₀,y₀)和P_end(x₁,y₁)，根据直线方程求出两点连线的表达式y＝rx+b；首先通过计算Δx＝|x₀-x₁|和Δy＝|y₀-y₁|来确定移动的主方向：

如果Δx>Δy，说明x轴的最大差值大于y轴的最大差值，x轴方向为移动的主方向；如果Δx<Δy，说明y轴的最大差值大于x轴的最大差值，y轴方向为移动的主方向；将起始点坐标分别对网格边长向下取整，即令

其中a是网格边长；

当x轴为移动的主方向时，设当前点为P′_start(x′₀,y′₀)，其中x′₀＝ai，y′₀＝aj，当x＝x′₀+a时，根据直线公式算出y＝k(x′₀+a)+b，对y以网格边长向下取整得，

令σ＝y′-y′₀，则：

计算

步后，当x＝x′₁，即终点所在网格时，求出起点和终点连线所在的外包网格空域集合G_base，当y轴为移动的主方向时同理；

S53、滑动邻域法求解扩展空域；

设S为初始解，即基本空域集合，邻域的个数为M，N_l表示第l个邻域；在进行邻域搜索时，只要当前邻域搜索结束后优化了S，即表现为效益函数增大时，则将当前邻域加入最优解，并进行下一个邻域结构进行运算，否则停止运算，输出当前解；

根据取整法得到航迹点连线穿过的网格空域集合G_base，以G_base作为初始解S；将初始解向两侧依次平移M次，得到M个递增的邻域，作为邻域结构集，分别记为N₁,N₂,…,N_M；

邻域结构间的搜索顺序可以通过改变邻域结构之间的次序，一般来说会按照邻域的大小顺序进行排序，从小到大，即：N₁(x)≤N₂(x)≤N₃(x)≤…≤N_max(x)；根据实际情况，邻域结构之间的次序为以初始解为基准向两侧依次交叉排列；

S54、定义效益函数F_obj表示某规划阶段所分配的网格空域的最大效益，其中该空域集合由基本空域和扩展空域两部分组成，分别为航迹中心线所在网格空域和与其相邻的部分网格空域；用f_dis表示增加当前扩展空域带来的自由度效益，f_c表示该网格空域产生的复杂性负效益，g_ij为决策变量，表示当增加某网格带来的总效益为正时，则将其划分为扩展空域，效益函数表达式为：

其中，由于目标为航线运行规划，考虑到距离航迹点连线过于远的网格对于地理围栏的规划没有实用意义，因此趋向于在预测航迹附近寻找复杂性等级较低的连续网格空域集合，即采用当前空域网格与预测航迹所在直线的距离的倒数来表示自由度效益：

其中，μ为自由度可调参数，d为当前空域网格与预测航迹所在直线的距离，s为复杂度值，λ用来保证规划空域网格的时候，距离航线太远的网格即使它复杂度很低也不把它划进规划队列，λ的值根据实际情况设置，σ根据实际需求用来调节复杂度带来负效益的影响，σ的值根据实际情况设置；

由于自由度效益函数为双曲线中的第一象限中的单支，因此随着距离d的增大，自由度效益f_dis会逐渐减小，设置μ为自由度可调参数，当d大于μ时，自由度效益小于1，同时，由于复杂度值由经过标准化后的指标数据加权平均值表示，考虑所有指标对复杂性的影响程度相同，则复杂度值的取值范围为s∈[0,1]，令λ＝1.5μ,σ＝1，则f_c∈[1,1.5μ+1]，当d大于μ时，效益函数F_obj≤0；

模型约束条件为：

决策变量为0-1变量，表示当前网格是否在最优解集中，当g_ij为1时，表示最优解集包含当前网格，当g_ij为0时，则不包括当前网格；

g_ij＝0,1；

作为无人机完成飞行任务的前提条件，必须有一个基本的网格空域满足当前飞行任务的需要，在模型中则要求起始点和终点连线所在网格必须为最优解集中；

g_ij＝1(若g_ij∈G_base)；

其中，G_base是起点和终点连线所在的外包网格空域集合；

可行域范围约束，若起始点和终点所在网格分别为

和

i₀和j₀表示起始点所在网格编号的下标，i₁和j₁表示终点所在网格编号的下标，i和j表示可行域网格编号的下标；则当起点和终点所在直线斜率r<0时：

斜率r>0时：

斜率r＝0时：

i₀≤i≤i₁,i,j∈N；

斜率r不存在时：

j₀≤j≤j₁,i,j∈N；

连通性约束，已知起点和终点连线的直线方程为Ax+By+C＝0，其中A,B不同时为0，则当Ai+Bj+C≤0时，即网格位于连线下方时：

当Ai+Bj+C≥0时，即网格位于连线下方时：

其中，g_i,j+1是下标为(i,j+1)的空域网格对应的决策变量，g_i+1,j是下标为(i+1,j)的空域网格对应的决策变量，g_i-1,j是下标为(i-1,j)的空域网格对应的决策变量，g_i,j-1是下标为(i,j-1)的空域网格对应的决策变量；

S55、利用S4提出的空域复杂分级度量方法和S54提出的效益函数，计算基本空域网格的效益函数是否为正，若为正，则表示将基本空域网格集合向两侧依次滑动，若为负，则输出所有扩展空域网格集合。

2.根据权利要求1所述的基于低空空域网格化的无人机动态地理围栏规划方法，其特征在于，步骤S1具体为：

首先考虑无人机性能及飞行计划数据，综合推算出水平轨迹、高度剖面以及速度剖面，从而获得无人机初始4D航迹；在此基础上，考虑无人机实际航迹随机特性，即风速的不确定性和指示空速的不稳定性导致实际的平均地速与计算得到的结果产生的差异；

假设无人机在巡航阶段的指示空速V_TAS和风速V_WS均服从正态分布，其中指示空速

风速

其概率密度函数分别为：

其中，μ_TAS表示平均指示空速，σ_TAS表示指示空速标准差，μ_WS表示风速平均值，σ_WS表示风速标准差；

结合正态分布的特点和平均地速的表达式可知平均地速V_GS服从正态分布，即

且其概率密度函数为：

由此推算出无人机在某航段上的平均空速实际取值在μ±3σ的范围内，即

其中DA为偏流角，表示航迹线与航向线的偏离程度，WA为风角，表示风向线同航迹线之间的夹角。

3.根据权利要求1所述的基于低空空域网格化的无人机动态地理围栏规划方法，其特征在于，步骤S2具体为：

在飞行习惯和风速的综合影响下，无人机在t时刻的平均地速在[μ_GS-3σ_GS,μ_GS+3σ_GS]范围内，则对应的动态地理围栏的纵向长度

为：

因此无人机动态地理围栏的纵向长度为6tσ_GS，同时会随着规划时间步长增加而增加，符合长期规划时不确定性增加的基本事实，其中，μ_GS表示平均地速，σ_GS表示地速的标准差。

4.根据权利要求1所述的基于低空空域网格化的无人机动态地理围栏规划方法，其特征在于，步骤S3具体为：

首先选取目标空域内的历史航迹数据，然后对历史航迹数据进行预处理，剔除异常数据、非目标空域范围内的数据和非低空空域内的航迹点数据；然后，依据墨卡托投影将目标空域转为平面直角坐标系，然后在平面坐标系上对目标空域进行均匀离散化，得到二维平面上的均匀空域网格并进行编号；然后对广播式自动相关监视数据格式的历史航迹数据进行解析处理，获得每个空域网格的原始指标数据；并进行编号；以坐标原点为起点，主方向为以x轴方向顺序递增，次方向为以y轴方向递增；同时，已知每个网格空域的单位长度，取每个网格的左下角坐标表示其网格角标。

5.根据权利要求1所述的基于低空空域网格化的无人机动态地理围栏规划方法，其特征在于，步骤S4具体为：

S41、在空域网格化的基础上引入航迹历史数据，将航迹密度，航迹点速度分布，航迹点高度分布，机型混合程度，潜在冲突程度，超障高度以及航迹时间密度几个方面作为影响低空空域复杂性的关键指标；采用线性加权法对关键指标进行处理，得到网格空域复杂度的综合评价值，即为复杂度综合评价模型；线性加权法为：

其中，CR_i是第i个网格空域复杂度的综合评价值，ω_j是第j个指标的权重，F′_j是第j个指标的评价值，n表示选取的影响低空空域复杂性的关键指标的个数；

S42、结合单位网格多个复杂性影响指标的高维特征和综合评价结果，引入核函数理论并结合K-均值聚类算法对空域网格复杂性进行分级量化研究；

首先对多指标综合量化结果进行处理，根据实际问题确定类簇个数K，并根据样本总数m将样本均分为K份，标记为C＝{c₁,c₂,…,c_K}，取每个样本中综合量化值中位数作为初始聚类中心；然后将原始数据经过非线性变换，映射到高维特征空间，再在特征空间中进行K-均值聚类；

已知选取的空域网格样本总数为m，则样本集为{g⁽¹⁾,g⁽²⁾,…,g^(m)}，其中g⁽ⁱ⁾∈G^N，i＝1,2,3,...m，依据Mercer定理，若已知映射Φ:

使得K(g⁽ⁱ⁾,g^(j))＝Φ(g⁽ⁱ⁾)^TΦ(g^(j))，K(g⁽ⁱ⁾,g^(j))则能够用来表示核函数，其中G^N是选取的空域网格样本g⁽ⁱ⁾所在的空间，G^F是选取的空域网格样本g⁽ⁱ⁾经过映射Φ后得到的Φ(g⁽ⁱ⁾)所在的空间；利用核函数改变样本的分布形式，而不必知道其在映射高维空间中的具体表现形式；核K-均值聚类就是讨论原始样本集{g⁽¹⁾,g⁽²⁾,…,g^(m)}经过映射Φ后得到的数据集{Φ(g⁽¹⁾),Φ(g⁽²⁾),…,Φ(g^(m))}在G^F空间中的聚类情况；样本空间中第i个样本到第k类中心的距离为：

以综合量化值对样本进行顺序排序，按照设置的复杂性等级数量、各等级样本比例，初步划分样本的复杂性等级；然后将样本映射到高维核空间，取每个复杂性等级样本的综合量化值的均值作为初始聚类中心。

Images