基于飞行状态划分的非合作无人机轨迹分布预测方法

1.本发明属于无人机交通管理技术领域，具体涉及一种基于飞行状态划分的非合作无人机轨迹分布预测方法。


背景技术：

2.随着无人机在各行业的迅速普及，低空无人机发展迅速且飞行需求与日俱增，空中交通活动日趋频繁。对无人机轨迹预测技术的研究，可以推动无人机在有限的低空空域安全、有序、高效地运行，计算无人机在空域中某一时刻的轨迹分布是构建危险行为识别模型的基础。
3.现有轨迹预测方法大多得到无人机未来的位置估计点，而非无人机所在的位置区间，在进行无人机运行安全风险分析时(入侵固定区域或冲突风险判断)，有可能会因为该位置预测结果存在的误差而忽略原本的危险行为，或产生虚警。此外，现有无人机的轨迹分布预测大多依据轨迹可达最大范围或部分已知条件约束下的轨迹可达范围，而在实际需要实时轨迹预测的任务场景中，可能面向的是首次出现的非合作类目标，其身份信息在实施监控管理时往往是未知的，对该类无人机操纵者的意图进行建模预估十分困难。


技术实现要素：

4.针对于上述现有技术的不足，本发明的目的在于提供一种基于飞行状态划分的非合作无人机轨迹分布预测方法，以解决现有技术中的轨迹预测方法在面向非合作类目标时无人机操纵者的意图难以预估、未来轨迹点预测精度不高导致漏警或虚警的问题；本发明方法通过无人机初始飞行状态划分，结合概率分布预测替代未来轨迹点位预测，生成轨迹分布空间，减少漏警、虚警的发生，同时缩小了非零概率的轨迹空间范围，实现更加精准的轨迹分布预测。
5.为达到上述目的，本发明采用的技术方案如下：
6.本发明的一种基于飞行状态划分的非合作无人机轨迹分布预测方法，步骤如下：
7.(1)针对所选监控空域的范围，构建栅格化空域；
8.(2)划分无人机飞行状态：依据无人机在初始时刻的飞行状态，将无人机的飞行状态划分为悬停状态和运动状态，设定无人机轨迹分布预测参数；
9.(3)筛选相似轨迹数据集，构建基于数据迁移的轨迹预测模型；
10.(4)生成轨迹可达空间：根据步骤(2)中划分的无人机飞行状态，求解无人机在预测时刻所有可能的轨迹位置的集合，作为具有闭合边界的区域，生成无人机在预测时刻的轨迹可达空间；
11.(5)索引轨迹覆盖栅格坐标：将步骤(4)生成的轨迹可达空间在栅格化空域中近似表达为其覆盖的栅格集合，求解所覆盖的栅格在空域笛卡尔坐标系中的坐标；
12.(6)生成轨迹概率分布：求解监测时刻无人机位于空域笛卡尔坐标系中各栅格的概率，得到栅格集合相对应的概率分布。
13.进一步地，所述步骤(1)中的监控空域范围以左下角为原点建立空域笛卡尔坐标系ll，通过设定x轴和y轴的栅格数量，将空域划分为带有唯一坐标的栅格集合，栅格坐标为(i,j)表示该栅格是空域x轴方向第i个，y轴方向第j个栅格，每个栅格保存的信息包括：栅格中心点和四个顶点在空域笛卡尔坐标系ll中的坐标。
14.进一步地，所述步骤(2)的具体过程如下：
15.(21)划分无人机飞行状态；
16.依据无人机在初始时刻的飞行状态，将无人机的飞行状态划分为悬停状态和运动状态；
17.(22)设定无人机轨迹分布预测参数；
18.(221)设定无人机飞行初始速度；
19.针对无人机在轨迹预测过程中的初始速度，设定悬停状态初始速度绝对值为零，运动状态初始速度绝对值大于零；
20.(222)设定无人机最大飞行地速及最大水平加速度；
21.根据无人机的出厂性能参数，预估无人机的静风下最大允许飞行速度均小于30m/s，设定最大飞行地速为30m/s，最大水平加速度为6m/s2；
22.(223)设定无人机轨迹预测时长；
23.根据步骤(222)中设定的无人机的最大飞行地速、最大水平加速度、最大静风速度和无人机危险行为识别的预先时间需求，设定预测时长大于无人机从悬停状态以最大水平加速度加速到最大飞行地速需要的时间。
24.进一步地，所述步骤(3)的具体过程如下：
25.(31)筛选相似轨迹数据集；
26.采集非合作类无人机的轨迹片段，根据已有的合作类无人机轨迹数据集，通过飞行速度方差和累积航向改变量来表达无人机轨迹的潜在特征，进而描述每段轨迹飞行状态的稳定性，基于飞行速度方差和累积航向改变量数值的相似，衡量轨迹之间的距离，筛选出与非合作无人机轨迹相似的轨迹数据集；
27.(32)构建基于数据迁移的轨迹预测模型；
28.利用步骤(31)中筛选出的相似轨迹数据集作为d-gru轨迹预测模型的训练样本，将训练后的d-gru轨迹预测模型视为基于数据迁移的轨迹预测模型。
29.进一步地，所述步骤(4)的具体过程如下：
30.(41)时空约束布朗桥；
31.针对步骤(21)中划分的不同飞行状态下无人机的运动特点，依据截断布朗桥表示出t时刻无人机的位置为(x(t)，y(t))时，x(t)服从截断概率密度函数：
[0032][0033]
t时刻无人机在轨迹可达空间内到达坐标对任意的y(t)服从截断概率密度函数：
[0034][0035]
此外，无人机在轨迹可达空间内到达坐标的概率为：
[0036][0037]
式中，(x(t)，y(t))为t时刻无人机的位置，对于正态分布x～n(μ，σ2)，为概率密度函数，φ(x)为分布函数，u
x
(t)、uy(t)是x，y方向的上边界，l
x
(t)、ly(t)是x，y方向的下边界；
[0038]
(42)初始悬停状态下无人机轨迹可达空间；
[0039]
对于初始时刻处于悬停状态的无人机，在预测时刻的轨迹可达空间为一个圆，圆心为初始时刻的位置，半径为到预测时刻的最大飞行距离，最大飞行距离包含加速飞行距离和以最大飞行速度匀速飞行距离两部分；
[0040]
(43)初始运动状态下无人机轨迹可达空间；
[0041]
(431)确定最小转弯半径；
[0042]
依据无人机在设计时允许的最大过载系数得到其最小转弯半径：
[0043][0044]
式中，r
min
为无人机的最小转弯半径，g表示重力加速度，|v0|表示开始转弯时无人机的飞行速度大小，n
max
为最大过载系数，设定具有初始速度的无人机在初始时刻开始转弯，转弯时初始速度保持大小不变，只改变方向；
[0045]
(432)设定转弯角度的方向；
[0046]
将转弯前后无人机速度方向之差定义为转弯角度，右转弯时，转弯角度为正；左转弯时，转弯角度为负；保持初始速度方向不变时，转弯角度为0；
[0047]
(433)建立无人机笛卡尔坐标系；
[0048]
以初始时刻无人机位置为坐标原点，无人机的初始速度的方向为纵轴y的正向，机身右侧为横轴x的正向，建立无人机笛卡尔坐标系oxy；
[0049]
(434)按照不同运动阶段分析轨迹可达空间；
[0050]
基于步骤(433)中建立的坐标系，对无人机转弯阶段按照步骤(432)中转弯角度的方向求解无人机转弯结束后的坐标；对于无人机加速阶段，设定由初始速度以最大水平加速度加速到最大飞行速度；对于无人机匀速飞行阶段，则以最大飞行速度匀速飞行；求解无人机在预测时刻即匀速飞行阶段最后的位置；
[0051]
(435)离散化轨迹可达空间的边界；
[0052]
预测时刻无人机的轨迹可达空间边界由转弯角度在区间(-360，360)下无人机在该时刻的边界坐标集合生成，离散化转弯角度，将转弯角度以等间隔度数进行采样，前一飞行阶段的终点为后一阶段的起点生成轨迹可达空间；
[0053]
(436)清除轨迹可达空间重叠区域；
[0054]
求解无人机在转弯和加速运行的时间长度内，其飞行轨迹首次出现交点时加速轨迹的转弯角度大小：
[0055][0056]
式中，sa为无人机加速阶段的飞行距离，r
min
为无人机的最小转弯半径；清除过大角度转弯时出现重叠的轨迹可达空间；
[0057]
(437)重新生成不同运动阶段的轨迹可达空间；
[0058]
将步骤(436)中得到的转弯角度区间以等间隔度数进行采样，生成初始运动在不同阶段覆盖的轨迹可达空间。
[0059]
进一步地，所述步骤(5)的具体过程如下：
[0060]
(51)初始悬停状态轨迹覆盖栅格坐标索引；
[0061]
(511)生成初始覆盖空间；
[0062]
依据步骤(42)中圆的外接正方形生成初始覆盖空间，保存外接正方形左上，左下，右上，右下四个顶点坐标集合；
[0063]
(512)扩展初始覆盖空间；
[0064]
基于步骤(511)中初始覆盖空间边界可能出现栅格覆盖不完整的情况，对于x，y方向下边界的栅格坐标向下求整，对于上边界的栅格坐标向上求整得到扩展后的覆盖空间；
[0065]
(513)求解栅格集合在栅格化空域笛卡尔坐标系ll中的坐标集合；
[0066]
遍历步骤(512)中覆盖空间包含的栅格，保留栅格中心到无人机初始位置的距离小于或等于最大飞行距离的栅格，得到栅格集合在栅格化空域笛卡尔坐标系ll中的坐标集合；
[0067]
(52)初始运动状态轨迹覆盖栅格坐标索引；
[0068]
(521)构建离散边界坐标集合；
[0069]
基于无人机初始时刻运动状态，得到预测时刻其轨迹可达空间基于步骤(433)中无人机笛卡尔坐标系oxy中的离散边界坐标集合；
[0070]
(522)坐标系转换；
[0071]
针对步骤(521)无人机笛卡尔坐标系oxy中的每一个离散边界坐标(x
c，β
，y
c，β
)，转换到步骤(1)中栅格所在的空域笛卡尔坐标系ll中的坐标(x
lc，β
，y
lc，β
)：
[0072][0073]
式中，α
l，o
表示由空域笛卡尔坐标系ll逆时针旋转到无人机笛卡尔坐标系oxy的角度，初始时刻无人机位置坐标为(xo，yo)，得到空域笛卡尔坐标系ll上的边界坐标集合；
[0074]
(523)求得覆盖空间；
[0075]
按照步骤(511)生成初始覆盖空间，然后根据步骤(512)扩展初始覆盖空间，即得到轨迹可达空间最小整数化栅格包含的覆盖空间；
[0076]
(524)曲线拟合圆方程；
[0077]
基于步骤(43)得到的轨迹可达空间，将栅格中心坐标(x
lck
，y
lck
)转移到无人机笛
卡尔坐标系oxy中的(x
ck
，y
ck
)：
[0078][0079]
式中，α
l，o
表示由空域笛卡尔坐标系ll逆时针旋转到无人机笛卡尔坐标系oxy的角度，初始时刻无人机位置坐标为(xo，yo)，进而计算无人机笛卡尔坐标系oxy中(x
ck
，y
ck
)与无人机初始位置o(0，0)点的连线与y正轴的夹角β
r，o
：
[0080][0081]
式中，(x
ck
，y
ck
)为任意栅格ck的中心在无人机笛卡尔坐标系oxy中的坐标，β
r，o
∈[-180，180]，进而通过最小二乘法构建圆弧拟合模型fit：
[0082]
(x
ck，c
，y
ck，c
)＝fit({x
c，β
，y
c，β
})，β∈[β
r，o-ε，β
r，o
ε]
ꢀꢀꢀꢀ
(9)
[0083]
模型输入为无人机预测时刻的边界坐标集合x
k，r
＝{(x
c，β
，y
c，β
)}，β∈[β
r，o-ε，β
r，o
ε]，ε用于定义角度范围，将边界坐标集合x
k，r
＝{(x
c，β
，y
c，β
)}中小的角度范围内的点视为圆弧，通过最小二乘法进行曲线拟合，输出为圆弧的圆心(x
ck，c
，y
ck，c
)，以及半径r
ck，c
，得到误差最小化对应的拟合圆方程；
[0084]
(525)判断栅格中心坐标与拟合圆的位置关系；
[0085]
计算栅格中心坐标与步骤(524)中拟合圆心之间的欧式距离，通过比较该欧式距离与拟合圆半径的大小，判断栅格中心坐标是否位于拟合圆弧内侧；将位于圆弧内侧的栅格视为在轨迹可达空间内部；
[0086]
(526)求解栅格集合在栅格化空域笛卡尔坐标系ll中的坐标集合；
[0087]
遍历步骤(523)中覆盖空间所含栅格，通过步骤(525)判断栅格中心坐标与拟合圆的位置关系，得到栅格集合在栅格化空域笛卡尔坐标系ll中的坐标集合；
[0088]
(53)表示索引的栅格坐标；
[0089]
将不同飞行状态下的栅格集合通过索引得到的坐标集合表示为集合{ck}，ck＝(ik，jk)表示栅格集合中栅格ck是空域笛卡尔坐标系ll中x方向第ik，y方向第jk个栅格。
[0090]
进一步地，所述步骤(6)的具体过程如下：
[0091]
(61)初始悬停状态轨迹概率分布求解；
[0092]
(611)坐标系平移；
[0093]
初始时刻为悬停状态的无人机的后续运动视为无向布朗运动，将轨迹可达空间从原有空域笛卡尔坐标系ll中的(x，y)平移到以初始时刻无人机位置为原点的新坐标系(x
′
，y
′
)：
[0094][0095]
式中，(x0，y0)为初始时刻无人机的位置坐标；
[0096]
(612)求解轨迹概率分布；
[0097]
在步骤(611)平移后的坐标系中，通过公式(3)求得预测时刻无人机位置坐标的概率密度函数，通过公式(10)将步骤(513)栅格集合中的栅格坐标进行转换，栅格ck＝(ik，jk)的栅格中心，左下，左上，右上，右下顶点在平移坐标系中的坐标依次为{(x
ck
′
(t)，y
ck
′
(t))，(x
ck，l
′
(t)，y
ck，l
′
(t))，(x
ck，l
′
(t)，y
ck，u
′
(t))，(x
ck，r
′
(t)，y
ck，u
′
(t))，(x
ck，r
′
(t)，y
ck，l
′
(t))}，求解预测时刻即t时刻无人机在轨迹可达空间内的位置位于栅格ck内的概率p
1，ck
(t)：
[0098][0099]
式中，x
ck，l
′
(t)，x
ck，r
′
(t)分别为栅格左上和左下，右上和右下顶点的横坐标，y
ck，l
′
(t)，y
ck，u
′
(t)分别为栅格左下和右下，左上和右上顶点的纵坐标，为初始悬停状态t时刻坐标的概率密度函数；
[0100]
(613)更新轨迹概率分布；
[0101]
依据步骤(612)得到的每个栅格内的概率值进行求和得到总概率p1(t)，p1(t)在栅格化足够精细时极其接近于1：
[0102][0103]
式中，p
1，ck
(t)为初始悬停状态t时刻无人机在轨迹可达空间内的位置位于栅格ck内的概率，为使模型扩展到栅格粒度较大的场景，通过公式(12)对p
1，ck
(t)进行更新：
[0104][0105]
得到无人机初始悬停状态时，t时刻无人机位于栅格集合中每个栅格的概率集合；
[0106]
(62)初始运动状态轨迹概率分布求解；
[0107]
(621)求解预测时刻轨迹可达空间的期望；
[0108]
初始时刻为运动状态的无人机的后续运动视为无向布朗运动和有向布朗运动的结合，利用步骤(32)中基于数据迁移的轨迹预测模型求解无人机在预测时刻的轨迹预测值，将其视为预测时刻轨迹可达空间的期望e＝(μ
x
(t)，μy(t))；
[0109]
(622)平移旋转坐标系；
[0110]
将轨迹可达空间从栅格化空域笛卡尔坐标系ll平移到以步骤(621)中得到的预测时刻无人机位置期望值e为原点，且旋转到以初始时刻的坐标与e的连线为新坐标系exy的x轴正向，进行原坐标系在新坐标系下的坐标转换：
[0111]
[0112]
式中，α
l，e
表示空域笛卡尔坐标系ll与exy坐标系x轴正向之间的夹角；
[0113]
(623)确定预测时刻无人机在轨迹可达空间内横坐标的上下限；
[0114]
将空域笛卡尔坐标系ll上的边界坐标集合转换到步骤(522)中的坐标系exy中，表示出预测时刻无人机在轨迹可达空间内横坐标的上下限；
[0115]
(624)基于横坐标确定纵坐标的上下限；
[0116]
计算步骤(623)中横坐标与转换后的边界坐标集合的差值绝对值，求解该集合中最小值对应的索引i
x，max
及该索引对应的yi，若yi≥0，设定纵坐标上限为yi，执行循环，设置该集合索引i
x，max
处为 ∞，求解该集合最小值对应的索引i
x，max
，直到索引i
x，max
对应的yi＜0，停止循环，设定纵坐标下限为yi；若yi＜0，设定纵坐标下限为yi，执行循环，设置该集合索引i
x，max
处为 ∞，求解该集合最小值对应的索引i
x，max
，直到索引i
x，max
对应的yi≥0，停止循环，设定纵坐标上限为yi；
[0117]
(625)求解预测时刻概率密度函数；
[0118]
根据步骤(623)中确定的横坐标的上下限和步骤(624)中确定的纵坐标的上下限，通过公式(3)求得初始状态为运动时预测时刻无人机位置坐标在轨迹可达空间内的概率密度函数；
[0119]
(626)选择概率密度函数的积分范围；
[0120]
基于步骤(622)旋转后栅格的顶点坐标将导致二元积分域变大，计算真实栅格积分域内部基于坐标系旋转角度得到的面积s1，建立真实栅格积分域s
ck
与旋转后栅格二元积分域面积s2之间的关系：
[0121][0122]
(627)求解轨迹概率分布；
[0123]
在步骤(622)平移旋转后的坐标系中，通过公式(14)将步骤(526)栅格集合中的栅格坐标进行转换，栅格ck＝(ik，jk)的栅格中心，左下，左上，右上，右下顶点在平移坐标系中的坐标依次为{(x
ck
′
(t)，y
ck
′
(t))，(x
ck，ll
′
(t)，y
ck，ll
′
(t))，(x
ck，lu
′
(t)，y
ck，lu
′
(t))，(x
ck，ru
′
(t)，y
ck，ru
′
(t))，(x
ck，rl
′
(t)，y
ck，rl
′
(t))}，依据步骤(626)中得到的关系，分别建立面积s1和s2上的二元积分p
2，ck1
(t)和p
2，ck2
(t)：
[0124][0125][0126]
式中，(x
ck
′
(t)，y
ck
′
(t))为栅格ck的中心坐标，a和b分别为面积s1中两个相邻边长的二分之一，u
ck，x
′
(t)和l
ck，x
′
(t)分别为栅格ck左下，左上，右上，右下四个顶点横坐标中的最大值和最小值，u
ck，y
′
(t)和l
ck，y
′
(t))分别为栅格ck左下，左上，右上，右下四个顶点纵坐标中的最大值和最小值，为初始运动状态t时刻坐标的概率密度函数，根据公式(15)得到初始运动状态t时刻无人机在轨迹可达空间内的位置位于栅格ck内的概率p
2，ck
(t)：
[0127][0128]
(628)求解概率分布集合；
[0129]
按照步骤(613)中的概率更新公式(13)对步骤(627)中的概率p
2，ck
(t)进行更新，得到无人机初始运动状态时，无人机位于栅格集合中每个栅格的概率集合。
[0130]
本发明的有益效果：
[0131]
本发明的方法提高了无人机轨迹信息不可知时其轨迹分布预测的精度，为无人机飞行过程中的危险行为识别和冲突风险评估提供依据，为低空空中交通安全管理提供理论支持。
[0132]
本发明方法通过将无人机的飞行状态细分，结合基于数据迁移的轨迹预测方法，考虑无人机操纵者意图不确定性的运动建模为布朗运动，采用截断布朗桥方法对非合作类无人机的位置分布进行建模，实现更加精准的轨迹分布预测，缩小非零概率的轨迹空间范围。
附图说明
[0133]
图1为本发明方法的流程图。
[0134]
图2a为实施例中初始悬停状态在预测时长为10s时的轨迹可达空间示意图。
[0135]
图2b为实施例中初始悬停状态在预测时长为20s时的轨迹可达空间示意图。
[0136]
图2c为实施例中初始悬停状态在预测时长为30s时的轨迹可达空间示意图。
[0137]
图3a为本发明实施例中初始运动状态在预测时长为5s时转弯阶段的轨迹可达空间示意图。
[0138]
图3b为本发明实施例中初始运动状态在预测时长为5s时加速阶段的轨迹可达空间示意图。
[0139]
图3c为本发明实施例中初始运动状态在预测时长为5s时匀速阶段的轨迹可达空间示意图。
[0140]
图4a为本发明实施例中初始运动状态在预测时长为30s时转弯阶段的轨迹可达空间示意图。
[0141]
图4b为本发明实施例中初始运动状态在预测时长为30s时加速阶段的轨迹可达空间示意图。
[0142]
图4c为本发明实施例中初始运动状态在预测时长为30s时匀速阶段的轨迹可达空间示意图。
[0143]
图5为本发明实施例中轨迹可达空间s
p,t
,、初始覆盖空间s
p,t,1
及最小整数化栅格s
p,t,2
关系示意图。
[0144]
图6a为本发明实施例中坐标系旋转角为12
°
时栅格积分区间示意图。
[0145]
图6b为本发明实施例中坐标系旋转角为50
°
时栅格积分区间示意图。
具体实施方式
[0146]
为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。
[0147]
参照图1所示，本发明的一种基于飞行状态划分的非合作无人机轨迹分布预测方法，步骤如下：
[0148]
(1)针对所选监控空域的范围，构建栅格化空域；
[0149]
其中，所述监控空域范围以左下角为原点建立空域笛卡尔坐标系ll，通过设定x轴和y轴的栅格数量，将空域划分为带有唯一坐标的栅格集合，栅格坐标为(i,j)表示该栅格是空域x轴方向第i个，y轴方向第j个栅格，每个栅格保存的信息包括：栅格中心点和四个顶点在空域笛卡尔坐标系ll中的坐标。
[0150]
(2)划分无人机飞行状态：依据无人机在初始时刻的飞行状态，将无人机的飞行状态划分为悬停状态和运动状态，设定无人机轨迹分布预测参数；
[0151]
(21)划分无人机飞行状态；
[0152]
依据无人机在初始时刻的飞行状态，将无人机的飞行状态划分为悬停状态和运动状态；
[0153]
(22)设定无人机轨迹分布预测参数；
[0154]
(221)设定无人机飞行初始速度；
[0155]
针对无人机在轨迹预测过程中的初始速度，设定悬停状态初始速度绝对值为零，运动状态初始速度绝对值大于零；
[0156]
(222)设定无人机最大飞行地速及最大水平加速度；
[0157]
根据无人机的出厂性能参数，预估无人机的静风下最大允许飞行速度均小于30m/s，设定最大飞行地速为30m/s，最大水平加速度为6m/s2；
[0158]
(223)设定无人机轨迹预测时长；
[0159]
根据步骤(222)中设定的无人机的最大飞行地速、最大水平加速度、最大静风速度和无人机危险行为识别的预先时间需求，设定预测时长大于无人机从悬停状态以最大水平加速度加速到最大飞行地速需要的时间。
[0160]
(3)筛选相似轨迹数据集，构建基于数据迁移的轨迹预测模型；
[0161]
(31)筛选相似轨迹数据集；
[0162]
采集非合作类无人机的轨迹片段，根据已有的合作类无人机轨迹数据集，通过飞行速度方差和累积航向改变量来表达无人机轨迹的潜在特征，进而描述每段轨迹飞行状态的稳定性，基于飞行速度方差和累积航向改变量数值的相似，衡量轨迹之间的距离，筛选出与非合作无人机轨迹相似的轨迹数据集；
[0163]
(32)构建基于数据迁移的轨迹预测模型；
[0164]
利用步骤(31)中筛选出的相似轨迹数据集作为d-gru轨迹预测模型的训练样本，将训练后的d-gru轨迹预测模型视为基于数据迁移的轨迹预测模型。
[0165]
(4)生成轨迹可达空间：根据步骤(2)中划分的无人机飞行状态，求解无人机在预测时刻所有可能的轨迹位置的集合，作为具有闭合边界的区域，生成无人机在预测时刻的轨迹可达空间；
[0166]
(41)时空约束布朗桥；
[0167]
针对步骤(21)中划分的不同飞行状态下无人机的运动特点，依据截断布朗桥表示出t时刻无人机的位置为(x(t)，y(t))时，x(t)服从截断概率密度函数：
[0168][0169]
t时刻无人机在轨迹可达空间内到达坐标对任意的y(t)服从截断概率密度函数：
[0170][0171]
此外，无人机在轨迹可达空间内到达坐标的概率为：
[0172][0173]
式中，(x(t)，y(t))为t时刻无人机的位置，对于正态分布x～n(μ，σ2)，为概率密度函数，φ(x)为分布函数，u
x
(t)、uy(t)是x，y方向的上边界，l
x
(t)、ly(t)是x，y方向的下边界；
[0174]
(42)初始悬停状态下无人机轨迹可达空间；
[0175]
对于初始时刻处于悬停状态的无人机，在预测时刻的轨迹可达空间为一个圆，圆心为初始时刻的位置，半径为到预测时刻的最大飞行距离，最大飞行距离包含加速飞行距离和以最大飞行速度匀速飞行距离两部分；预测时长取10s、20s、30s，以60度为间隔采样的无人机的轨迹可达空间如图2a-图2c，图中所示线条为以60度为间隔采样的6个无人机航向，靠近圆心较细的线条为无人机加速飞行阶段，靠近外侧较粗的线条为无人机以最大飞行速度匀速飞行的阶段，随着危险行为识别预测时长的增加，以最大飞行速度匀速飞行的距离占半径的比重会显著增大；
[0176]
(43)初始运动状态下无人机轨迹可达空间；
[0177]
(431)确定最小转弯半径；
[0178]
依据无人机在设计时允许的最大过载系数得到其最小转弯半径：
[0179][0180]
式中，r
min
为无人机的最小转弯半径，g表示重力加速度，|v0|＝10m/s表示开始转弯时无人机的飞行速度大小，n
max
＝3.5为最大过载系数，设定具有初始速度的无人机在初始时刻开始转弯，转弯时初始速度保持大小不变，只改变方向；
[0181]
(432)设定转弯角度的方向；
[0182]
将转弯前后无人机速度方向之差定义为转弯角度，右转弯时，转弯角度为正；左转弯时，转弯角度为负；保持初始速度方向不变时，转弯角度为0；
[0183]
(433)建立无人机笛卡尔坐标系；
[0184]
以初始时刻无人机位置为坐标原点，无人机的初始速度的方向为纵轴y的正向，机身右侧为横轴x的正向，建立无人机笛卡尔坐标系oxy；
[0185]
(434)按照不同运动阶段分析轨迹可达空间；
[0186]
基于步骤(433)中建立的坐标系，对无人机转弯阶段按照步骤(432)中转弯角度的方向求解无人机转弯结束后的坐标；对于无人机加速阶段，设定由10m/s以6m/s2加速到30m/s；对于无人机匀速飞行阶段，则以最大飞行速度匀速飞行；求解无人机在预测时刻即匀速飞行阶段最后的位置；
[0187]
(435)离散化轨迹可达空间的边界；
[0188]
预测时刻无人机的轨迹可达空间边界由转弯角度在区间(-360，360)下无人机在该时刻的边界坐标集合生成，离散化转弯角度，将转弯角度以5度为间隔度数进行采样，前一飞行阶段的终点为后一阶段的起点生成轨迹可达空间；
[0189]
(436)清除轨迹可达空间重叠区域；
[0190]
求解无人机在转弯和加速运行的时间长度内，其飞行轨迹首次出现交点时加速轨迹的转弯角度大小：
[0191][0192]
式中，sa为无人机加速阶段的飞行距离，r
min
为无人机的最小转弯半径；清除过大角度转弯时出现重叠的轨迹可达空间，重叠区域的出现是由于转弯角度绝对值大于β0的边界坐标集合xd＝{(x
c，β
，y
c，β
)}，β∈(-360，-β0)∪(β0，360)被包含在转弯角度绝对值小于等于β0的边界坐标集合xk＝{(x
c，β
，y
c，β
)}，β∈[-β0，β0]生成的闭合空间内，因此需要清除转弯角度绝对值大于β0时的边界坐标集合。
[0193]
(437)重新生成不同运动阶段的轨迹可达空间；
[0194]
基于步骤(436)得到的转弯角度区间，基于此区间以5度为间隔进行采样，生成初始运动预测时长为5s、30s时在不同阶段覆盖的轨迹可达空间如图3a-图4c，前一飞行阶段的终点为后一阶段的起点，对于旋翼无人机，转弯阶段的移动范围相比于加速和匀速更小，随着预测时长的增加，转弯时间导致的无人机到达最大飞行速度前的位移损失相比匀速阶段的总位移占比将越少，即预测时间跨度越大，无人机的轨迹可达空间的边界越接近一个圆。
[0195]
(5)索引轨迹覆盖栅格坐标：将步骤(4)生成的轨迹可达空间在栅格化空域中近似表达为其覆盖的栅格集合，求解所覆盖的栅格在空域笛卡尔坐标系中的坐标；
[0196]
(51)初始悬停状态轨迹覆盖栅格坐标索引；
[0197]
(511)生成初始覆盖空间；
[0198]
依据步骤(42)中圆的外接正方形生成初始覆盖空间，保存外接正方形左上，左下，右上，右下四个顶点坐标集合；
[0199]
(512)扩展初始覆盖空间；
[0200]
基于步骤(511)中初始覆盖空间边界可能出现栅格覆盖不完整的情况，对于x，y方向下边界的栅格坐标向下求整，对于上边界的栅格坐标向上求整得到扩展后的覆盖空间；如图5所示，轨迹可达空间s
p，t
、初始覆盖空间s
p，t，1
以及最小整数化栅格s
p，t，2
之间呈逐层包含的关系；
[0201]
(513)求解栅格集合在栅格化空域笛卡尔坐标系ll中的坐标集合；
[0202]
遍历步骤(512)中覆盖空间包含的栅格，保留栅格中心到无人机初始位置的距离小于或等于最大飞行距离的栅格，得到栅格集合在栅格化空域笛卡尔坐标系ll中的坐标集合；
[0203]
(52)初始运动状态轨迹覆盖栅格坐标索引；
[0204]
(521)构建离散边界坐标集合；
[0205]
基于无人机初始时刻运动状态，得到预测时刻其轨迹可达空间基于步骤(433)中无人机笛卡尔坐标系oxy中的离散边界坐标集合；
[0206]
(522)坐标系转换；
[0207]
针对步骤(521)无人机笛卡尔坐标系oxy中的每一个离散边界坐标(x
c，β
，y
c，β
)，转换到步骤(1)中栅格所在的空域笛卡尔坐标系ll中的坐标(x
lc，β
，y
lc，β
)：
[0208][0209]
式中，α
l，o
表示由空域笛卡尔坐标系ll逆时针旋转到无人机笛卡尔坐标系oxy的角度，初始时刻无人机位置坐标为(xo，yo)，得到空域笛卡尔坐标系ll上的边界坐标集合；
[0210]
(523)求得覆盖空间；
[0211]
按照步骤(511)生成初始覆盖空间，然后根据步骤(512)扩展初始覆盖空间，即得到轨迹可达空间最小整数化栅格包含的覆盖空间；
[0212]
(524)曲线拟合圆方程；
[0213]
基于步骤(43)得到的轨迹可达空间，将栅格中心坐标(x
lck
，y
lck
)转移到无人机笛卡尔坐标系oxy中的(x
ck
，y
ck
)：
[0214][0215]
式中，α
l，o
表示由空域笛卡尔坐标系ll逆时针旋转到无人机笛卡尔坐标系oxy的角度，初始时刻无人机位置坐标为(xo，yo)，进而计算无人机笛卡尔坐标系oxy中(x
ck
，y
ck
)与无人机初始位置o(0，0)点的连线与y正轴的夹角β
r，o
：
[0216][0217]
式中，(x
ck
，y
ck
)为任意栅格ck的中心在无人机笛卡尔坐标系oxy中的坐标，β
r，o
∈[-180，180]，进而通过最小二乘法构建圆弧拟合模型fit：
[0218]
(x
ck，c
，y
ck，c
)＝fit({x
c，β
，y
c，β
})，β∈[β
r，o-ε，β
r，o
ε]
ꢀꢀꢀ
(9)
[0219]
模型输入为无人机预测时刻的边界坐标集合x
k，r
＝{(x
c，β
，y
c，β
)}，β∈[β
r，o-ε，β
r，o
ε]，ε用于定义角度范围，将边界坐标集合x
k，r
＝{(x
c，β
，y
c，β
)}中小的角度范围内的点视为圆弧，通过最小二乘法进行曲线拟合，输出为圆弧的圆心(x
ck，c
，y
ck，c
)，以及半径r
ck，c
，得到误
差最小化对应的拟合圆方程；
[0220]
(525)判断栅格中心坐标与拟合圆的位置关系；
[0221]
计算栅格中心坐标与步骤(524)中拟合圆心之间的欧式距离，通过比较该欧式距离与拟合圆半径的大小，判断栅格中心坐标是否位于拟合圆弧内侧；将位于圆弧内侧的栅格视为在轨迹可达空间内部；
[0222]
(526)求解栅格集合在栅格化空域笛卡尔坐标系ll中的坐标集合；
[0223]
遍历步骤(523)中覆盖空间所含栅格，通过步骤(525)判断栅格中心坐标与拟合圆的位置关系，得到栅格集合在栅格化空域笛卡尔坐标系ll中的坐标集合；
[0224]
(53)表示索引的栅格坐标；
[0225]
将不同飞行状态下的栅格集合通过索引得到的坐标集合表示为集合{ck}，ck＝(ik，jk)表示栅格集合中栅格ck是空域笛卡尔坐标系ll中x方向第ik，y方向第jk个栅格。
[0226]
(6)生成轨迹概率分布：求解监测时刻无人机位于空域笛卡尔坐标系中各栅格的概率，得到栅格集合相对应的概率分布；
[0227]
(61)初始悬停状态轨迹概率分布求解；
[0228]
(611)坐标系平移；
[0229]
初始时刻为悬停状态的无人机的后续运动视为无向布朗运动，将轨迹可达空间从原有空域笛卡尔坐标系ll中的(x，y)平移到以初始时刻无人机位置为原点的新坐标系(x
′
，y
′
)：
[0230][0231]
式中，(x0，y0)为初始时刻无人机的位置坐标；
[0232]
(612)求解轨迹概率分布；
[0233]
在步骤(611)平移后的坐标系中，通过公式(3)求得预测时刻无人机位置坐标的概率密度函数，通过公式(10)将步骤(513)栅格集合中的栅格坐标进行转换，栅格ck＝(ik，jk)的栅格中心，左下，左上，右上，右下顶点在平移坐标系中的坐标依次为{(x
ck
′
(t)，y
ck
′
(t))，(x
ck，l
′
(t)，y
ck，l
′
(t))，(x
ck，l
′
(t)，y
ck，u
′
(t))，(x
ck，r
′
(t)，y
ck，u
′
(t))，(x
ck，r
′
(t)，y
ck，l
′
(t))}，求解预测时刻即t时刻无人机在轨迹可达空间内的位置位于栅格ck内的概率p
1，ck
(t)：
[0234][0235]
式中，x
ck，l
′
(t)，x
ck，r
′
(t)分别为栅格左上和左下，右上和右下顶点的横坐标，y
ck，l
′
(t)，y
ck，u
′
(t)分别为栅格左下和右下，左上和右上顶点的纵坐标，为初始悬停状态t时刻坐标的概率密度函数；
[0236]
(613)更新轨迹概率分布；
[0237]
依据步骤(612)得到的每个栅格内的概率值进行求和得到总概率p1(t)，p1(t)在栅格化足够精细时极其接近于1：
[0238][0239]
式中，p
1，ck
(t)为初始悬停状态t时刻无人机在轨迹可达空间内的位置位于栅格ck内的概率，为使模型扩展到栅格粒度较大的场景，通过公式(12)对p
1，ck
(t)进行更新：
[0240][0241]
得到无人机初始悬停状态时，t时刻无人机位于栅格集合中每个栅格的概率集合；
[0242]
(62)初始运动状态轨迹概率分布求解；
[0243]
(621)求解预测时刻轨迹可达空间的期望；
[0244]
初始时刻为运动状态的无人机的后续运动视为无向布朗运动和有向布朗运动的结合，利用步骤(32)中基于数据迁移的轨迹预测模型求解无人机在预测时刻的轨迹预测值，将其视为预测时刻轨迹可达空间的期望e＝(μ
x
(t)，μy(t))；
[0245]
(622)平移旋转坐标系；
[0246]
将轨迹可达空间从栅格化空域笛卡尔坐标系ll平移到以步骤(621)中得到的预测时刻无人机位置期望值e为原点，且旋转到以初始时刻的坐标与e的连线为新坐标系exy的x轴正向，进行原坐标系在新坐标系下的坐标转换：
[0247][0248]
式中，α
l，e
表示空域笛卡尔坐标系ll与exy坐标系x轴正向之间的夹角；
[0249]
(623)确定预测时刻无人机在轨迹可达空间内横坐标的上下限；
[0250]
将空域笛卡尔坐标系ll上的边界坐标集合转换到步骤(522)中的坐标系exy中，表示出预测时刻无人机在轨迹可达空间内横坐标的上下限；
[0251]
(624)基于横坐标确定纵坐标的上下限；
[0252]
计算步骤(623)中横坐标与转换后的边界坐标集合的差值绝对值，求解该集合中最小值对应的索引i
x，max
及该索引对应的yi，若yi≥0，设定纵坐标上限为yi，执行循环，设置该集合索引i
x，max
处为 ∞，求解该集合最小值对应的索引i
x，max
，直到索引i
x，max
对应的yi＜0，停止循环，设定纵坐标下限为yi；若yi＜0，设定纵坐标下限为yi，执行循环，设置该集合索引i
x，max
处为 ∞，求解该集合最小值对应的索引i
x，max
，直到索引i
x，max
对应的yi≥0，停止循环，设定纵坐标上限为yi；
[0253]
(625)求解预测时刻概率密度函数；
[0254]
根据步骤(623)中确定的横坐标的上下限和步骤(624)中确定的纵坐标的上下限，通过公式(3)求得初始状态为运动时预测时刻无人机位置坐标在轨迹可达空间内的概率密度函数；
[0255]
(626)选择概率密度函数的积分范围；
[0256]
基于步骤(622)旋转后栅格的顶点坐标将导致二元积分域变大，计算真实栅格积
分域内部基于坐标系旋转角度得到的面积s1如图6a和图6b，坐标系旋转角α
l，e
分别为12
°
和50
°
，对建立在exy坐标系上的概率密度函数进行积分时，选择旋转后栅格ck的顶点坐标将导致二元积分域变大，即图6a和图6b中的面积s2，真实的栅格ck的积分域应当为图6a和图6b中的面积s
c，k
，建立真实栅格积分域s
ck
与旋转后栅格二元积分域面积s2之间的关系：
[0257][0258]
(627)求解轨迹概率分布；
[0259]
在步骤(622)平移旋转后的坐标系中，通过公式(14)将步骤(526)栅格集合中的栅格坐标进行转换，栅格ck＝(ik，jk)的栅格中心，左下，左上，右上，右下顶点在平移坐标系中的坐标依次为{(x
ck
′
(t)，y
ck
′
(t))，(x
ck，ll
′
(t)，y
ck，ll
′
(t))，(x
ck，lu
′
(t)，y
ck，lu
′
(t))，(x
ck，ru
′
(t)，y
ck，ru
′
(t))，(x
ck，rl
′
(t)，y
ck，rl
′
(t))}，依据步骤(626)中得到的关系，分别建立面积s1和s2上的二元积分p
2，ck1
(t)和p
2，ck2
(t)：
[0260][0261][0262]
式中，(x
ck
′
(t)，y
ck
′
(t))为栅格ck的中心坐标，a和b分别为面积s1中两个相邻边长的二分之一，u
ck，x
′
(t)和l
ck，x
′
(t)分别为栅格ck左下，左上，右上，右下四个顶点横坐标中的最大值和最小值，u
ck，y
′
(t)和l
ck，y
′
(t))分别为栅格ck左下，左上，右上，右下四个顶点纵坐标中的最大值和最小值，为初始运动状态t时刻坐标的概率密度函数，根据公式(15)得到初始运动状态t时刻无人机在轨迹可达空间内的位置位于栅格ck内的概率p
2，ck
(t)：
[0263][0264]
(628)求解概率分布集合；
[0265]
按照步骤(613)中的概率更新公式(13)对步骤(627)中的概率p
2，ck
(t)进行更新，得到无人机初始运动状态时，无人机位于栅格集合中每个栅格的概率集合。
[0266]
本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。