上海理工大学机械工程学院，上海

收稿日期：2022年6月13日；录用日期：2022年7月14日；发布日期：2022年7月20日

摘要

高效精确地对无人机飞行轨迹进行预测是无人机领域的关键技术之一。针对现在无人机飞行轨迹预测方法的预测精度不足的问题，本文在进行大地坐标与空间直角坐标转换的基础上，提出一种基于最小二乘的无人机飞行轨迹预测算法，在无人机的飞行过程中，该算法能够不间断地更新轨迹数据并拟合飞行轨迹方程，从而达到对飞行轨迹的精确预测。仿真实验结果表明，所提出的轨迹预测算法误差可控制在0.5 m以内，相比较于BP神经网络预测模型，该算法的可靠性和精度更优。

关键词

无人机，轨迹预测，BP神经网络，最小二乘法

Research on UAV Trajectory Prediction Based on Improved Least Squares Method

Qianqian Yuan, Song Huang, Ibulayimujiang Yasen

School of Mechanical Engineering, University of Shanghai for Science and Technology, Shanghai

Received: Jun. 13^th, 2022; accepted: Jul. 14^th, 2022; published: Jul. 20^th, 2022

ABSTRACT

Efficient and accurate prediction of UAV flight trajectory is one of the key technologies in the field of UAV. To address the problem of insufficient prediction accuracy of current UAV flight trajectory prediction methods, this paper proposes a least squares-based UAV flight trajectory prediction algorithm based on the conversion of geodetic coordinates and spatial Cartesian coordinates, which can uninterruptedly update the trajectory data and fit the flight trajectory equation during the flight of UAV, so as to achieve accurate prediction of flight trajectory. The simulation experimental results show that the error of the proposed trajectory prediction algorithm can be controlled within 0.5 m, and the reliability and accuracy of the algorithm is better compared with the BP neural network prediction model.

Keywords:UAV, Trajectory Prediction, BP Neural Network, Least Square Method

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

由于无人机的大量投入与广泛应用，导致无人机有更大风险闯入重要场所，对设备及人身安全构成了严重的威胁，所以对重要场所设置禁飞区且对禁飞区进行规避是必要的。如何通过已知历史位置数据预测无人机的轨迹，从而实现自主避让禁飞区，对保护禁飞区的安全具有重要意义。

目前，国内外研究学者对飞行轨迹的预测已展开了诸多的研究。这些研究主要从神经网络 [1]、模型滤波和数据挖掘 [2] [3] 这几个方面入手，进而进行轨迹预测算法的研究。文献 [4] 结合了机器学习的方法和RM-IMM的估计方法，从历史监控数据中训练出基于机器学习的轨迹预测模型，从中可以得到数据驱动的预测结果；文献 [5] 提出了一种用于飞行轨迹预测的受约束的长期内记忆网络，其模型能够与动态物理限制保持长期依赖性。Baklacioglu T等 [6] 利用遗传算法对运输机起飞上升和着陆下降段的轨迹预测模型进行改进，提升了模型预测精准度；文献 [7] 对Kalman滤波算法进行了改进，该算法对预测模型中的系统噪声进行了4D运动轨迹估计，并且相比较于传统的滤波算法，它的优势在于实时地对系统噪声进行评估，更有效地提高了轨迹预测的精度。且算法模型简单，可以稳定地按照环境变化调节自身对环境的适应能力，但该模型学习时间过长，在长期进行大量的数据集训练模型时，没有很好的性能并且预测差值偏大；文献 [8] 提出一种在数据挖掘和机器学习的基础上解决船舶轨迹预测问题。该算法通过对2000多条历史轨迹信息进行研究，利用数据挖掘的空间聚类算法对船舶的历史轨迹进行聚类，分析出这些历史轨迹的分布规律。但这种算法有一定缺陷，它响应时间较长，计算量大且需要占用大量内存对历史轨迹信息进行储存，因此对单架无人机某次飞行进行轨迹曲线预测时并不适用该算法。

以上这些方法虽然都能解决轨迹预测问题，但也各有不足。IMM算法是一种需要建立模型的估计预测算法，可以在移动目标产生多种不同机动时，通过切换到各自对应的运动模型来匹配移动目标的实际运动状态，但是也比较容易受到建模误差所带来的影响；基于Kalman滤波的轨迹预测算法虽然能够对目标进行轨迹预测，但算法需要建模，而且会产生滤波器发散的问题 [9]。而对无人机进行轨迹预测时，需要考虑到无人机是否会进入禁飞区对设备及人身安全造成威胁。本文分别针对BP神经网络预测模型 [10] [11] 和传统最小二乘法 [12] 的飞行轨迹预测进行了研究，同时在传统最小二乘法的基础上，对获得的离散轨迹点采取下一时刻舍去最末尾的点同时加入最新的轨迹点进行曲线拟合 [13]，进而预测出无人机的飞行轨迹。

2. 坐标转换

在进行曲线拟合之前，需要先获取无人机的轨迹点坐标。本文通过无人机上的机载GPS接收机，采集一段WGS-84坐标系下的飞行轨迹经纬度坐标点。为了保证后续轨迹预测的精确性，需要对采集到的经纬度进行处理，将其转变为CGCS2000坐标系下的平面直角坐标 [14]。

本文研究中，假定无人机定高飞行，采集的是二维离散飞行轨迹点，然后通过七参数法和高斯投影法将经纬度坐标(B, L)转换为高斯平面坐标(x, y)，此过程中未考虑高程，但这并不影响结果。通过MATLAB编写代码进行坐标转换，表1为通过MATLAB进行坐标转换前后的部分数据。

表1. 坐标转换前后的部分数据

3. 基于BP神经网络预测模型的飞行轨迹预测

由于BP神经网络模型为单输入单输出数据，在进行预测时，选取输入层和输出层的神经元个数均为1，根据经验公式 $l=\sqrt{m+n}+\alpha$ 可以计算出隐含层神经元个数在3~11个之间，通过不断改变神经元的数目观察样本均方根误差的变化，从图1可知选取隐含层神经元个数为7。

图1. 隐含层不同神经元个数训练结果

训练模型参数选择如表2所示。

表2. 训练模型参数

在进行BP神经网络预测时，使用上文中给出的参数进行训练，BP神经网络训练结果如图2所示。

图2. BP神经网络训练结果

图3为BP神经网络轨迹预测模型误差。从图中可以看出，BP神经网络预测的轨迹误差在5 m之内，对于执行任务中的无人机来说，该误差偏大。

图3. BP神经网络轨迹预测模型误差

4. 基于更新轨迹点的最小二乘法飞行轨迹研究

由于BP神经网络的预测误差偏大，本文在传统最小二乘法的基础上，对获得的离散轨迹点采取下一时刻舍去最末尾的点同时加入最新的轨迹点进行曲线拟合，从而达到提高预测精确度的效果 [15]。基于传统最小二乘法进行的无人机轨迹预测是移动目标单次运动时根据历史轨迹点数据对运动曲线进行预测。轨迹预测流程图如图4所示。

图4. 轨迹预测流程图

由GPS所测得的离散轨迹点只能表示出无人机当前的运动方向及运动位置。如果可以通过这些轨迹点找出精确地描述无人机的运动曲线函数表达式，就能够很好地预测出无人机的运动轨迹了。因此，本文通过无人机GPS接收机采集到实时的无人机飞行轨迹点，然后对数据进行处理后采取更新轨迹点的方法对无人机此时的运动轨迹进行最小二乘法拟合，不断利用此操作，反映出无人机下一时刻的运动轨迹。同时将其与传统的最小二乘法和BP神经网络预测模型进行比较。

4.1. 曲线拟合多项式阶次选择

首先通过无人机的机载GPS接收机获取无人机的离散轨迹点，本文截取了一组120个原始轨迹点的数据集。对该数据集根据经纬度坐标到空间直角坐标的转换关系，利用matlab软件对采集到的GPS轨迹点坐标进行数据转换，为了方便计算，本文将转换后的平面直角坐标系原点从(0, 0)平移至(41362000, 3461900)处。

然后再使这120个轨迹点作为拟合采样点，分别采用二阶拟合多项式、三阶拟合多项式、五阶拟合多项式来拟合无人机的飞行轨迹。仿真结果如图5所示。

图5. 传统最小二乘法曲线拟合结果

从图中可以明显的看出，五阶拟合效果更优，因此本文决定在五阶多项式拟合的基础上，对所采集的离散坐标点进行去掉此时最末尾的轨迹点加入下一时刻最新的轨迹点的拟合方法。

4.2. 基于传统最小二乘法的飞行轨迹预测

在得到一组离散数据后，通过选择合适的拟合表达式的待定参数，让数据点的拟合值与实际测量值的残差的加权平方和最小。对于给定一组离散点数据：

$\left(x_i,y_i\right)\left(i=1,2,\cdots ,m\right)$ (1)

可以用曲线多项式对这些数据进行描述，函数表达式为：

$y\left(x\right)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n\left(n<m\right)$ (2)

其中，n是多项式的阶数。为了确定拟合效果性能，本文通过确定参数a的值使拟合多项式与给定的离散数据在的残差 ${\delta }_i=f\left(x_i\right)-y_i$ 的加权平方和达到最小，即：

${\left|\delta \right|}^2={{\displaystyle {\sum }_i^m\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n\right)\right]}}^2$ (3)

采用最小二乘曲线拟合法，根据多元函数求极值，函数在某点有偏导数且为零，可得：

$\frac{\partial {\delta }^2}{\partial a_i}=0$ (4)

其矩阵形式如下：

$\left[\begin{array}{cccc}N& {\displaystyle {\sum }_{i=1}^N{x}_i}& \cdots & {\displaystyle {\sum }_{i=1}^N{x}_i^n}\\ {\displaystyle {\sum }_{i=1}^N{x}_i}& {\displaystyle {\sum }_{i=1}^N{x}_i^2}& \cdots & {\displaystyle {\sum }_{i=1}^N{x}_i^{n+1}}\\ ⋮& \text{}\text{}⋮& \ddots & ⋮\\ {\displaystyle {\sum }_{i=1}^N{x}_i^n}& {\displaystyle {\sum }_{i=1}^N{x}_i^{n+1}}& \cdots & {\displaystyle {\sum }_{i=1}^N{x}_i^{n+n}}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}{\displaystyle {\sum }_{i=1}^N{y}_i}\\ {\displaystyle {\sum }_{i=1}^N{y}_i}a\\ ⋮\\ {\displaystyle {\sum }_{i=1}^N{y}_i}{a}_n\end{array}\right]$ (5)

求解上述方程，得到系数矩阵 ${\left[\begin{array}{cccc}{a}_0& a_1& \cdots & a_n\end{array}\right]}^{\text{T}}$ 的唯一解，将 ${\left[\begin{array}{cccc}{a}_0& a_1& \cdots & a_n\end{array}\right]}^{\text{T}}$ 带入式(3)即可得到多项式的最优解。

根据以上步骤，首先将截取的120个原始轨迹点的数据集进行坐标转换后得到一组离散坐标点 $\left\{\begin{array}{cc}\left(x_i,y_i\right)& \left(i=0,1,2,\cdots ,119\right)\end{array}\right\}$

在进行公式(2)的求解时，本文选用5阶多项式进行曲线拟合，则所需求解的函数表达式为：

$y\left(x\right)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+a_4{x}^4+a_5{x}^5$ (6)

在得到无人机飞行轨迹表达式 $y\left(x\right)$ 函数之后，需要确定 $y\left(x\right)$ 的系数 $a_0,a_1,a_2,a_3,a_4,a_5$ 的值。根据公式(3)~(4)求误差平方和对系数得偏导，并且使偏导数为零。可得到下列方程组：

$\{\begin{array}{l}-2{\displaystyle \underset{i=0}{\overset{119}{\sum }}\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdots +a_5{x}_i^5\right)\right]}=0\hfill \\ -2{\displaystyle \underset{i=0}{\overset{119}{\sum }}\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdots +a_5{x}_i^5\right)\right]x_i=0}\hfill \\ ⋮\hfill \\ -2{\displaystyle \underset{i=0}{\overset{119}{\sum }}\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdots +a_5{x}_i^5\right)\right]x_i^5=0}\hfill \end{array}$ (7)

由于6个未知数 $a_0,a_1,a_2,a_3,a_4,a_5$，6个方程，所以此方程组可解。将其整理为正规方程组得：

$\{\begin{array}{l}119a_0+a_1{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i}+a_2{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i^2}+\cdots +a_5{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i^5}={\displaystyle \underset{i=0}{\overset{119}{\sum }}{y}_i}\hfill \\ a_0{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i+a_1{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i^2+a_2{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i^3+\cdots +a_5{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i^6={\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i{y}_i}}}}}\hfill \\ ⋮\hfill \\ a_0{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i^5+a_1{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i^6+a_2{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i^7+\cdots +a_5{\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i^{10}={\displaystyle \underset{i=0}{\overset{119}{\sum }}{x}_i^5{y}_i}}}}}\hfill \end{array}$ (8)

最后，通过MATLAB对所采集的120个离散坐标点进行传统的最小二乘法求解时，所得出的系数 $a_0,a_1,a_2,a_3,a_4,a_5$ 的最优解[−41.9924, 4.3074, −0.0275, 6.4738e−05, 6.2397e−08, 2.0729e−11]，获得的无人机飞行轨迹预测曲线函数表达式为：

$y\left(x\right)=-41.9924+4.3074x-0.0275x^2+\left(6.4738\text{e}-05\right)x^3+\left(6.2397\text{e}-08\right)x^3+\left(2.0729\text{e}-11\right)x^5$ (9)

五阶多项式的曲线拟合结果如图6所示。

图6. 五阶多项式曲线拟合结果

图7为传统最小二乘算法的曲线拟合误差。从图中可以看出，传统最小二乘算法的轨迹预测误差在30 m之内。

4.3. 更新轨迹点最小二乘法预测算法

采用传统的最小二乘法对无人机离散轨迹点进行轨迹拟合时可以看出，该算法可以反映出无人机大致的飞行轨迹曲线，也可对未来一段时间无人机的轨迹进行预测，但是预测误差偏差明显。这是由于未来短时间内的飞行轨迹只与此刻位置和最近一段时间位置信息有关，与更早的位置信息无关。基于这一思想，本文在最小二乘法的基础上，去掉此时最末尾的轨迹点加入下一时刻最新的轨迹点，即更新轨迹点的方法重新拟合新的轨迹曲线。以此不间断对无人机飞行轨迹进行预测，保证了更新轨迹点的实时性和准确性。更新轨迹点的最小二乘曲线拟合轨迹预测算法步骤如下：

1) 确定无人机飞行轨迹的离散坐标点 $\left\{\begin{array}{cc}\left(x_i,y_i\right)& i=0,1,2,\cdots ,m\end{array}\right\}$ ；

图7. 传统最小二乘法曲线拟合误差

2) 根据公式(3)~(5)，利用前r个离散坐标点对系数矩阵 ${\left[\begin{array}{cccc}{a}_0& a_1& \cdots & a_n\end{array}\right]}^{\text{T}}$ 求最优解，获得该曲线的函数表达式；

3) 在第i (r ≥ 2)个轨迹点采样时刻去掉此时最末尾的轨迹点 $\left(x_{i-r},y_{i-r}\right)$，加入下一时刻最新的轨迹点 $\left(x_{i+1},y_{i+1}\right)$ 重新拟合出新的曲线；

4) 不断重复执行上述(1)、(2)、(3)步骤，对无人机飞行轨迹进行预测(图8)。

图8. 无人机飞行轨迹预测算法流程图

在进行更新轨迹点重新拟合曲线之前，需要考虑的问题是在更新轨迹点的过程中，每一次的轨迹拟合选取几个离散轨迹点可以使运行时间更短，预测误差更小。上文中选择了5阶拟合多项式来拟合无人机的飞行轨迹，由公式(6)可知离散轨迹点数目需大于多项式阶数，因此本文分别对6个，8个，10个和12个轨迹点数目进行比较，选取最合适的轨迹点数目进行无人机飞行轨迹拟合。

本文对无人机轨迹预测进行研究的目的是更高效精确地对禁飞区进行规避，从而保护变电站的设备安全。因此为了更直观地分析四种不同轨迹点数目的优劣性，表3给出了运行时间和平均绝对误差的仿真结果。

表3. 不同轨迹点数目的运行时间和平均绝对误差

由上表可以看出，当轨迹点数目为6时，平均绝对误差最小为0.0011 m，但是运行时间最长，不利于无人机飞行时快速预测轨迹；当轨迹点数目为12时，运行时间较短为98.29 s，但平均绝对误差较大，使得无人机在飞行轨迹预测中得到的结果准确性低。综合考虑两者的影响，可以发现轨迹点数目选取8时，运行时间不算太长且平均绝对误差为0.0437 m，该误差结果在可接受的范围之内，对于后期无人机进行禁飞区的规避任务也不会有较大影响。

为了后续更准确地将各算法进行比较，更新轨迹点的最小二乘法的仿真实验依然使用这120个离散坐标点进行，且该算法对当前时刻之前的8个轨迹坐标点进行飞行曲线拟合，函数表达式设为：

$p_1\left(x\right)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+a_4{x}^4+a_5{x}^5$ (10)

在去掉此时最末尾的轨迹坐标点 $\left(x_0,y_0\right)$，加入下一时刻最新的轨迹坐标点 $\left(x_8,y_8\right)$ 重新拟合出 $\left(x_7,y_7\right)$ 与 $\left(x_8,y_8\right)$ 之间的轨迹曲线，函数表达式设为：

$p_2\left(x\right)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+a_4{x}^4+a_5{x}^5$ (11)

以此类推，不断重复，最终无人机的轨迹曲线表达式由以下分段函数表示：

$P\left(x\right)=\{\begin{array}{cc}{P}_1\left(x\right)& p_0\sim p_7\\ P_2\left(x\right)& p_7\sim p_8\\ ⋮& \\ P_{113}\left(x\right)& p_{118}\sim p_{119}\end{array}$ (12)

由于在无人机对禁飞区进行规避时，飞行轨迹预测效果与当前时刻和最近一段时间的轨迹点相关。因此在多段拟合的曲线中，本文只关注最新一段轨迹的参数，所求得的 $P_{113}\left(x\right)$ 函数的参数a的值为[78492565.9851, −463117.9450, −1.2886, 1092.6514, −1.2886, −1.7903e−07]。分段函数为：

$P\left(x\right)=\{\begin{array}{ll}-159.0305+33.9409x-2.5079x^2+0.0948x^3-0.0017x^4+1.1497x^5\hfill & p_0\sim p_7\hfill \\ -94.1679+18.8801x-1.2072x^2+0.0422x^3-0.0007x^4+\left(4.2828\text{e}-06\right)x^5\hfill & p_7\sim p_8\hfill \\ ⋮\hfill & \hfill \\ 78492565.9851-463117.9450x-1.2886x^2+1092.6514x^3-1.2886x^4+\left(1.7903\text{e}-07\right)x^5\hfill & p_{118}\sim p_{119}\hfill \end{array}$ (13)

为了验证该预测算法的有效性和可行性，依然采用Matlab仿真软件实现对更新轨迹点的最小二乘法代码的编写。拟合曲线如图9所示。

图9. 更新轨迹点的最小二乘曲线拟合结果

图10为更新轨迹点最小二乘轨迹预测误差。从图中可以直观地看出，不断更新轨迹点进行最小二乘曲线拟合算法的轨迹预测误差范围不超过0.5 m。

图10. 更新轨迹点最小二乘轨迹预测误差

5. 算法研究对比分析

本文所采用的最小二乘算法、BP神经网络预测模型和更新轨迹点的最小二乘算法误差如图11所示。将三种预测方法的误差进行比较，可以看出：更新轨迹点的最小二乘轨迹预测算法的误差在0.5 m以内，明显小于传统的最小二乘法和BP神经网络模型预测的误差。

图11. 三种预测方法的误差对比

6. 结束语

随着GPS定位在无人机领域的应用，大地坐标与空间直角坐标之间的转换也成为了非常重要的任务，本文实现WGS84坐标与CGCS2000坐标之间的换算，对无人机今后的研究具有重要意义。同时文章采用基于更新轨迹点的最小二乘曲线拟合的方法对无人机飞行过程中采集的离散坐标点的轨迹进行预测，且将其与传统的最小二乘算法及BP神经网络预测模型进行对比，通过仿真结果显示，用最小二乘曲线拟合的方法在利用舍去最末尾的坐标点加上最新的坐标点的情况下，具有较高的预测精度，对今后无人机轨迹预测方向具有一定的实际应用价值。

文章引用

原倩倩,黄 松,伊卜拉伊木江•亚森. 基于改进最小二乘法的无人机轨迹预测研究
Research on UAV Trajectory Prediction Based on Improved Least Squares Method[J]. 建模与仿真, 2022, 11(04): 1131-1142. https://doi.org/10.12677/MOS.2022.114103

参考文献