CN115577511A

CN115577511A - 基于无人机运动状态的短期航迹预测方法、装置及系统

Info

Publication number: CN115577511A
Application number: CN202211171519.7A
Authority: CN
Inventors: 钟罡; 周蒋颖; 张洪海; 费毓晗; 刘皞
Original assignee: Nanjing University of Aeronautics and Astronautics
Current assignee: Nanjing University of Aeronautics and Astronautics
Priority date: 2022-09-26
Filing date: 2022-09-26
Publication date: 2023-01-06
Anticipated expiration: 2042-09-26
Also published as: CN115577511B

Abstract

本发明公开了一种基于无人机运动状态的短期航迹预测方法、装置及系统，所述短期航迹预测方法包括对获取到的无人机轨迹数据进行预处理，获得若干个轨迹点数据；基于所述轨迹点数据，判断无人机的运动模式，所述运动模式包括悬停状态和运动状态；当判定无人机的运动模式为运动状态时，则将所述轨迹点数据输入至轨迹预测模型，进行短期轨迹预测。本发明将无人机的运动状态看作为一个二分类问题，考虑了运动状态的变化对于轨迹预测精度的影响，实现了无人机短期预测精度的提高。

Description

基于无人机运动状态的短期航迹预测方法、装置及系统

技术领域

本发明属于无人机轨迹预测技术领域，具体涉及一种基于无人机运动状态的短期航迹预测方法、装置及系统。

背景技术

无人机航迹预测是冲突检测的核心技术，是实现无人机自动化监管的基础。航迹预测是指通过计算预估航空器未来航迹的过程。其数据来源包括性能数据、基础数据与气象数据。航迹预测主要应用于冲突探测、短时流量管理、任务规划等领域。此前，针对航空器的四维航迹预测已有了许多的研究，航空器的四维航迹预测一般有三种模型：状态估计模型、动力学模型以及机器学习模型。在深度学习的算法中，已经有各种类型的模型应用于时间序列预测任务，其中的模型主要包括人工神经网络(ANN)、深度信念网络(BND)、受限玻尔兹曼机(RBM)、循环神经网络(RNN)、长短期记忆(LSTM)。这些算法能在一定程度上对航迹时间特性进行有效的预测，但往往忽略了无人机运动的空间特性，因此，如何充分捕捉无人机的时空特性是问题的关键所在。

无人机飞行轨迹具有连续性和时序性的特点。连续性是指无人机轨迹是连续变化的，而不是间断的。时序性是指轨迹数据是带有时间性质的，且后一个时刻的位置与前一个时刻的位置是有关的，所以无人机的轨迹数据本质上是一个时间序列数据。无人机的飞行轨迹属于四维航迹的一种，在传统的无人机真实的轨迹数据集中，存在大量符合真实运动模式的静止轨迹，用现有的预测方法处理静止轨迹的输入数据时，通常都会输出一条运动的轨迹，导致预测误差较大。因此通过对无人机两种运动状态的识别和分类，获得无人机实时的运行状态变化。此外，由于轨迹点都是由三维经纬度及高度组成，在空间范围内，许多数据之间差异非常微小，导致直接使用单个模型无法很好的对数据建模。因此，利用DBN和双向长短期记忆网络联合进行预测来提高预测精度，通过DBN网络捕获轨迹数据中的局部特征，通过双向长短期记忆网络挖掘轨迹数据中的时间依赖性。

发明内容

针对上述问题，本发明提出一种基于无人机运动状态的短期航迹预测方法、装置及系统，将无人机的运动状态看作为一个二分类问题，考虑了运动状态的变化对于轨迹预测精度的影响，实现了无人机短期预测精度的提高。

为了实现上述技术目的，达到上述技术效果，本发明通过以下技术方案实现：

第一方面，本发明提供了一种基于无人机运动状态的短期航迹预测方法，包括：

对获取到的无人机轨迹数据进行预处理，获得若干个轨迹点数据；

基于所述轨迹点数据，判断无人机的运动模式，所述运动模式包括悬停状态和运动状态；

当判定无人机的运动模式为运动状态时，则将所述轨迹点数据输入至轨迹预测模型，进行短期轨迹预测。

可选地，所述无人机的运动模式通过以下步骤判断：

计算出无人机的行为特征值属性向量l_(1*6)，l_(1*6)＝[x',y',z',x”,y”,z”]，其中x、y、z分别表示轨迹点的经度、纬度和高度，x',y',z'分别表示轨迹点的经度的一阶导数、纬度的一阶导数和高度的一阶导数，x”,y”,z”分别表示轨迹点的经度的二阶导数、纬度的二阶导数和高度的二阶导数；

设置无人机运动状态约束条件，所述约束条件包括：

约束1：

约束2：h≥H_min

若所述行为特征值属性向量l_(1*6)中的一阶导数值以及二阶导数值全部为0，则表明无人机的状态为静止状态；反之，无人机处于运动状态；

当判定无人机处于静止状态后，进一步判断无人机的飞行高度h与允许高度最小值H_min之间的关系，若满足约束2，则表明无人机的状态为悬停状态，否则无人机的状态为地面静止状态；

计算无人机掉头所需要的时间t_max：

其中，F_合为无人机运动所需的力，m为无人机的质量，v为无人机的运行速度，r为转弯半径；

若无人机完成一次掉头所需的时间为t_s秒，则在t_s<t_max秒之内无人机不存在掉头飞回原点的情况。

可选地，所述短期航迹预测方法还包括：

当判定无人机的运动模式为悬停状态时，则直接输出所述轨迹点数据。

可选地，所述短期轨迹预测包括以下步骤：

获取轨迹预测模型，所述轨迹预测模型以深度信念网络为编码器，以双向长短期记忆网络为解码器；所述深度信念网络由三层受限波尔兹曼机组成，且采用ReLU作为激活函数；

分别提取各轨迹点数据中的经度、纬度、高度数据，建立对应的无人机空间轨迹属性向量p_(1*3)，p_(1*3)＝[x,y,z]，其中x、y、z分别表示轨迹点的经度、纬度和高度；

将若干个无人机空间轨迹属性向量组成空间属性向量矩阵，并输入到轨迹预测模型中，进行短期轨迹预测。

可选地，所述轨迹点数据的获取方法包括：

提取无人机轨迹数据中的经度、纬度、高度数据，并将提取到的数据从Geodetic坐标系转换到ECEF坐标系，获得初始轨迹点数据；

对初始轨迹点数据进行平滑滤波；

当判定在一段时间内，初始轨迹点数据中轨迹点的数目小于预设的预测轨迹所需轨迹点数目，则认为初始轨迹点数据稀疏，对所述初始轨迹点数据进行补全处理，得到最终轨迹点数据。

可选地，所述初始轨迹点数据通过以下步骤获得：

获取到的原始无人机轨迹数据包含P_G＝[X,Y,Z]^T，其中，X表示Geodetic坐标系下的经度，Y表示Geodetic坐标系下的纬度，Z表示Geodetic坐标系下的高度；

将P_G＝[X,Y,Z]^T从Geodetic坐标系转换到ECEF坐标系，获得P_E，P_E＝[X_E,Y_E,Z_E]^T＝[(σ+Z)cos Y cos X,(σ+Z)cos Y sin X,[σ(1-e²)+Z]sin Y]^T，

其中，X_E表示ECEF坐标系下的经度，Y_E表示ECEF坐标系下的纬度，Z_E表示ECEF坐标系下的高度；σ是基准椭球体的卯酉圆曲率半径，e是地球离心率，a是参考椭球的赤道半径。

可选地，所述对初始轨迹点数据进行去噪，包括以下步骤：

采用Savitzky-Golay滤波器对初始轨迹点进行平滑滤波，Savitzky-Golay滤波器在同一段曲线上的任意位置可以任意选取不同的窗口宽度，对该曲线进行平滑滤波，具体包括：

假设一段无人机轨迹共有N个轨迹点，确定滤波的窗口宽度是N＝2w+1，其中，w表示需要拟合的单边轨迹点数，将这些轨迹点的X_E、Y_E、Z_E、v分别取出组成若干个与时间相关的点集，v表示无人机的速度；

采用平滑滤波公式分别对X_E、Y_E、Z_E、v的点集序列进行平滑滤波；

将平滑滤波后点集序列中相应位置的X_E、Y_E、Z_E、v还原至每个轨迹点，并以此为组合得到每个点平滑后的单个轨迹点数据。

可选地，所述对所述初始轨迹点数据进行补全处理为：采用反距离加权法对平滑滤波得到的轨迹数据进行插值，插值计算公式为：

式中，λ_j为权重系数，d_j为所插轨迹点到所有轨迹点的距离，j为第j个所插轨迹点，t(x_o,y_o,z_o)为插值后函数，t(x_j,y_j,z_j)为插值点函数，N'为预设的预测轨迹所需的插值数目。

第二方面，本发明提供了一种基于无人机运动状态的短期航迹预测装置，包括:

预处理模块，用于对获取到的无人机轨迹数据进行预处理，获得若干个轨迹点数据；

判断模块，用于基于所述轨迹点数据，判断无人机的运动模式，所述运动模式包括悬停状态或运动状态；

预测模块，用于当判定无人机的运动模式为运动状态时，则将所述轨迹点数据输入至轨迹预测模型，进行短期轨迹预测。

第三方面，本发明提供了一种基于无人机运动状态的短期航迹预测系统，其特征在于，包括存储介质和处理器；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行根据第一方面中任一项所述方法的步骤。

与现有技术相比，本发明的有益效果：

本发明提出的一种基于无人机运动状态的短期航迹预测方法、装置及系统，根据无人机的运动特性划分无人机的运动模式，综合考虑了无人机的时间特性和空间特性的约束，构建最小损失函数预测模型，从中得到最优的组合预测方法，提高无人机短期轨迹的预测精度。

附图说明

为了使本发明的内容更容易被清楚地理解，下面根据具体实施例并结合附图，对本发明作进一步详细的说明，其中：

图1为本发明一种实施例的基于无人机运动状态的短期航迹预测方法的流程示意图。

图2为本发明一种实施例的BDN神经网络示意图，其中，BDN由三层RBM组成。

图3为本发明一种实施例的Bi-LSTM示意图。

具体实施方式

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

下面结合附图对本发明的应用原理作详细的描述。

实施例1

本发明实施例中提供了一种基于无人机运动状态的短期航迹预测方法，所述无人机可以是城区物流无人机，还可以是适用于其他场景的无人机，具体根据实际需要进行场景设计，如图1所示，所述短期航迹预测方法具体包括以下步骤：

步骤1：通过GPS采集四旋翼无人机的运动状态相关信息，获取飞行中无人机的速度、经纬度、高度、接收日期等信息，形成无人机轨迹数据；

步骤1.1：确定无人机飞行区域，获取飞行中无人机的相关数据P₀，所述相关数据包括：

P₀＝[X,Y,Z,v,date]^T，

其中，X表示经度，Y表示纬度，Z表示高度，v为速度，date为接收日期。

步骤1.2：对所述相关数据进行筛选，解决数据异常、数据丢失和其他待清理的问题，形成无人机轨迹数据。

步骤2：将所述无人机轨迹数据转化为无人机控制需要的状态信息；具体为：先把所述无人机轨迹数据从大地坐标系(Geodetic)转换成笛卡尔坐标系(ECEF)下的数据，获得无人机在平面中的位置和速度。

步骤2.1：从步骤1中的无人机轨迹数据中提取经度、纬度、高度数据，再将提取到的数据从Geodetic坐标系转换到ECEF坐标系。

所述从步骤1中的无人机轨迹数据中提取到的经度、纬度、高度数据的表达式为：

P_G＝[X,Y,Z]^T

其中，X表示Geodetic坐标系下的经度，Y表示Geodetic坐标系下的纬度，Z表示Geodetic坐标系下的高度；

将提取到的数据从Geodetic坐标系转换到ECEF坐标系，所采用的转换公式为：

P_E＝[X_E,Y_E,Z_E]^T＝[(σ+Z)cos Y cos X,(σ+Z)cos Y sin X,[σ(1-e²)+Z]sin Y]^T

步骤2.2：采用Savitzky-Golay滤波器对轨迹点进行平滑滤波。

Savitzky-Golay滤波器在同一段曲线上的任意位置可以任意选取不同的窗口宽度，对该曲线进行平滑滤波。假设一段无人机轨迹共有N个轨迹点，确定滤波的窗口宽度是N＝2w+1，其中，w表示需要拟合的单边轨迹点数。将这些轨迹点的X_E、Y_E、Z_E、v分别取出组成若干个与时间相关的点集；

以轨迹点X_E坐标组成的点集为例，则在窗口内的采样轨迹点集为i＝(-w,-w+1,...,0,...,w-1,w)，采用k-1次多项式对窗口内的数据点进行拟合；

平滑滤波公式：

其中：h_i为平滑系数，由最小二乘法拟合多项式求得，i代表采样轨迹点集，k代表拟合多项式最高次系数。

采用以上公式分别对每个序列进行平滑滤波后，再将平滑后序列中相应位置的X_E、Y_E、Z_E、v还原至每个轨迹点，并以此为组合得到每个点平滑后的单个轨迹点数据格式为P'：

P'＝[X,Y,Z,v]^T

步骤2.3：由步骤2.2得知一段无人机轨迹共有N个轨迹点,设定预测轨迹所需的轨迹点数目N'，如果在一段时间内，轨迹点的数目小于一定阈值，则认为轨迹比较稀疏，需要进行插值等补全处理。无人机定位数据是一条时间序列，各轨迹点间以时间间隔为距离计算权重，采用反距离加权法对轨迹数据进行插值，各样本点以时间间隔为距离计算权重，以t_c＝0.1s为插值间隔，对轨迹数据量进行扩充。其中，t_c为插值时间。

插值计算公式为：

式中，λ_j为权重系数，d_j为所插轨迹点到所有轨迹点的距离，j为第j个所插轨迹点，t(x_o,y_o,z_o)为插值后函数，t(x_j,y_j,z_j)为插值点函数，N'为预设的预测轨迹所需的轨迹点数目。

将每个插值按时间顺序连接在一起得到的插值后的每条轨迹的轨迹点个数为：

(P₁,P₂,...,P_N,...,P_N')

为了后续说明方便，将每条轨迹表示为：

(P₁,P₂,...,P_N,...,P_N')＝(p₁,p₂,...,p_t)

步骤3：将无人机的GPS数据转化为可利用的数据后，引入一个无人机运动状态判断的约束条件；

步骤3.1：通过平滑插值后的单个轨迹点数据表示具有时间戳，经度，纬度，高度，速度等基础信息，如式所示：

P_t＝[t',x_t,y_t,z_t,v_t]

步骤3.2：本发明主要将无人机的运动模式分为两种：悬停(静止)和运动。捕捉无人机轨迹的空间分布特征，提取无人机轨迹数据中的经纬度、高度数据，建立无人机空间轨迹属性向量，作为轨迹预测的输入。空间轨迹属性向量p入下式：

p_(1*3)＝[x,y,z]

其中，x、y、z分别表示无人机轨迹的经度、纬度、高度。

步骤4：判断无人机的运动状态，若在设定的时间段内无人机的状态为悬停(静止)，则不通过轨迹预测模型进行预测，直接输出该轨迹作为预测的无人机轨迹；若无人机为运动状态，则将引入一个BDN-Bi-LSTM轨迹预测模型进行轨迹预测；

步骤4.1：建立运动状态下无人机的行为特征值属性向量。本发明以无人机轨迹点三个方向的一阶导数值以及二阶导数值和高度表征行为特征，定义行为特征值属性向量l如下式：

l_(1*6)＝[x',y',z',x”,y”,z”,h]

步骤4.2：设定无人机运动状态约束条件，当行为特征值属性向量l中一阶导数值以及二阶导数值均为0时，表明无人机速度和加速度均为0：

[x',y',z',x”,y”,z”]＝0

步骤4.3：由于无人机升限限制和城市低空空域管理政策影响，物流无人机的飞行高度h应高于允许高度最小值H_min，此时无人机的状态为悬停状态，否则无人机处于地面静止状态，需满足：

h≥H_min

步骤4.4：无人机的运动所需的力由螺旋桨和重力的合力F_合提供；

其中：m为无人机的质量，v为无人机的运行速度，r为转弯半径。

步骤4.5：假设无人机做掉头运动，相当于无人机做了一次自身的圆周运动，则：

2πr＝vt_max

步骤4.6：合并以上两个式子，计算无人机掉头所需要的时间t_max：

步骤4.7：取无人机质量为m＝5kg，运动速度v＝20m/s，螺旋桨和重力的合力提供向心力F_合＝200N，计算结果t_max：

通过以上假设可知无人机完成一次掉头所需的时间为π秒，近似将其认为掉头所需的时间约为3秒，考虑到在3秒之内无人机不存在掉头飞回原点的情况，且满足行为特征约束条件和飞行最小高度约束条件的情况下无人机处于悬停状态，因此基于本发明实施例中的方法，能够预测运动无人机未来3秒内的轨迹。

步骤5：通过步骤4将运动的轨迹和悬停的轨迹进行区分，直接输出与时间序列关系不大的悬停运动模式的轨迹，将运动的轨迹用来做预测，建立组合神经网络预测模型，将运动轨迹输入到组合神经网络中进行综合预测，提高轨迹预测的精度。

步骤5.1：建立轨迹预测模型(BDN-Bi-LSTM)，采用深度信念网络和双向长短期记忆网络联合解决。以深度信念网络为编码器，双向长短期记忆网络为解码器。采用了编码-解码模型对无人机轨迹进行预测。

步骤5.2：判断步骤4中是否满足约束。若满足，则无人机处于悬停(静止)状态，混合模型直接输出预测轨迹P，预测长度为L。

P＝(p′_t+1,p′_t+2,...,p′_t+L)，p′_i＝p_t

其中，p′_i表示P中的每个轨迹点，均等于步骤2.3中的p_t。

步骤5.3：若不满足约束，则无人机处于运动状态，建立空间特征属性向量，将向量数据输入到BDN-Bi-LSTM轨迹预测模型中进行预测。

(1)捕捉无人机轨迹的空间分布特征，提取无人机轨迹数据中的经纬度、高度数据，建立无人机空间轨迹属性向量。空间轨迹属性向量p入下式：

p_(1*3)＝[x,y,z]

(2)当无人机被判定为运动无人机时，将无人机轨迹点(p₁,p₂,...,p_t)组成的空间属性向量矩阵输入到BDN-Bi-LSTM轨迹预测模型中，输入数据为：

其中：M指由t个三维输入轨迹点组成的t*3维矩阵。

(3)采用编码-解码模型对输入的轨迹序列进行预测，先采用深度信念网络(DBN)中的受限玻尔兹曼机(RBM)，捕捉无人机轨迹的局部特性。如图2所示，本发明设置的深度信念网络由三层受限波尔兹曼机组成，其中每层中的神经元数目均设置为128。每层受限波尔兹曼机的隐藏层维度分别为3*128、128*128和128*128，激活函数为ReLU。RBM模块定义如下：

其中W_i为RBM中第i层的权重矩阵，b_i为第i层的偏置系数。

(4)本发明通过将RBM隐藏单元设置为128，利用以上公式对三维的输入轨迹M进行非线性转换，把三维向量从3×t转换为128×t。

(5)将Bi-LSTM的编码器表示为ENC(·)，编码器将的128×t的输入向量编码成一个128维的中间向量H_t。

(6)其中，编码器采用ReLU作为激活函数，它的第t个单元的编码公式为：

f_t＝ReLU(W_pfp_t+W_hfh_t-1+b_f)

i_t＝ReLU(W_pip_t+W_hih_t-1+b_i)

O_t＝ReLU(W_pop_t+W_hoh_t-1+b_o)

c_t＝f_tc_t-1+i_t·tanh(W_pcp_t+W_hch_t-1+b_c)

其中，W_pf，W_pi，W_po，W_pc,

是输入向量p_t连接到每个对应门结构的权重矩阵,W_hf,W_hi,W_ho,W_hc是每一层连接到前一个短期状态h_t-1的权重矩阵,tanh为tanh激活函数，b_i,b_f,b_o,b_c,b_y为对应层的偏置系数，g＝(e^2p-1)/(e^2p+1)，f_t为LSTM的遗忘门结构，i_t和O_t分别为输入门和输出门结构,c_t用于记忆输入样本的关键数值特征，y_t为输出序列，

是前向序列，

是

前一个时刻的序列，

是后向序列，

是

后一个时刻的序列；

(7)将Bi-LSTM的解码器表示为DEC(·)对中间编码向量H_t进行解码，解码器中的第一个LSTM单元的输入数据为h_t和p_t，输出数据为o_t+1和h_t+1，将输出数据输入到后续的DEC(·)单元中，得到向量

模型如下式：

(8)输出向量

通过全连接层生成无人机的预测位置(p′_t+1,p′_t+2,...,p′_t+L)，L为预测长度。

p′_t+L＝W_Ot+j+b,W∈R^3×128

(9)由步骤五(2)和步骤五(8)可以到无人机得到预测的无人机轨迹P。

实施例2

基于与实施例1相同的发明构思，本发明实施例中提供了一种基于无人机运动状态的短期航迹预测装置，包括:

实施例3

基于与实施例1相同的发明构思，本发明实施例中提供了一种基于无人机运动状态的短期航迹预测系统，包括存储介质和处理器；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行根据实施例1中任一项所述方法的步骤。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。