无人机轨迹的高斯变分推断精确学习方法

1.本发明涉及非线性、批量化的状态估计领域，主要涉及无人机轨迹精确预测的高斯变分推断状态估计方法。


背景技术：

2.状态估计是无人机系统导航的核心部分。状态估计的目标是通过已知的能观测到的状态估计出系统中无法观测到的状态。当无人机系统模型参数较复杂、系统参数随着传感器变化而变化，且状态具有非线性、批量化的特点时，传统的状态估计方法存在误差大，算法迭代次数多等缺陷。
3.状态估计可以看作是对非线性系统的识别，故可以使用机器学习技术来进行处理。机器学习方法的优点在于不用对估计对象进行精确完整地建模，不考虑对象的运动细节，单纯从数据的角度出发，在估计对象模型以及运动简单时，这种方法能够精确地估计出运动对象的状态，但是当运动对象模型以及运动复杂时，无法进行长期稳定有效的状态估计。
4.另一种方法是对运动对象进行概率状态估计，运用卡尔曼滤波、贝叶斯滤波、sigmapoint滤波、粒子滤波、高斯过程回归、最大后验估计(maximum aposteriori,map)和高斯变分推断(gaussian variational inference,gvi)等方法估计运动对象在空间中移动时的状态。这些方法在点云对齐、位姿图松弛、光束平差法以及同时定位与地图构建等方面得到了重要应用。传统的gvi方法的特点是初始就确定优化目标的最佳均值和协方差，如果均值维数为n，协方差矩阵将变为n
×
n的矩阵，优化计算复杂度将会大幅提高。因此，gvi框架虽然具有计算效率高、误差小的优势，但是在针对无人机模型以及系统参数非线性、批量化的特点时，其无法进行此类的大批量状态估计。


技术实现要素：

5.本发明提出了无人机估计精确预测的高斯变分推断状态估计方法，改进gvi算法使其适用于批量状态估计，将map方法的算法结构结合到gvi方法中以构建完整的优化迭代算法。改进后的gvi方法将以小幅增加计算成本为代价，提高批量状态估计的精确性，进而进一步近似完整的贝叶斯后验。改进后的gvi方法通过kl散度找到最接近完整贝叶斯后验的高斯近似。与传统的gvi方法不同，本文将在算法中使用自然梯度下降法对均值以及协方差进行迭代更新，使改进后的gvi方法能够进行批量状态估计。相比于传统方法，改进后的gvi方法因为其在精确性方面有更好的表现，所以将该批量状态估计与轨迹预测方法应用到无人机轨迹预测中。
6.于此，本发明的目的在于改进一种无人机轨迹的高斯变分推断精确学习方法。本发明能够提升在无人机轨迹精确轨迹过程中大批量状态估计的准确性。
7.为实现上述目的，本发明提供的无人机轨迹的高斯变分推断精确学习方法包括如下步骤：
8.步骤s1：根据无人机在实际飞行中的空气动力学原理以及运动学原理构建无人机的数学模型；
9.步骤s2：基于高斯变分推断，通过kl散度将真实贝叶斯后验和近似后验联系起来，将状态的估计问题转化为后验的近似问题，再通过定义损失函数，将kl散度的计算转换为对损失函数的优化问题；
10.步骤s3：采用牛顿法更新损失函数以及采用自然梯度下降法对均值和协方差逆矩阵进行迭代优化，使损失函数v(q)降到最低以达到有效近似后验的目的。进一步的，步骤1中构建的无人机数学模型如下：
11.在惯性坐标系oxy中，假设无人机位于空间内任意一点，坐标为ξ＝[x，y，z]，欧拉角为其动力学方程式为：
[0012][0013]
式中，v表示无人机速度；n
x
，nz分别表示无人机切向过载以及法向过载；g表示重力加速度，取9.8m/s2。
[0014]
进一步的，步骤2中基于高斯变分推断，通过kl散度将真实贝叶斯后验和近似后验联系起来，得到近似后验和真实贝叶斯后验的kl散度展开式如下：
[0015][0016]
其中，是高斯分布的信息熵，ln p(z)是常数；
[0017]
将kl散度转化为损失函数，定义如下：
[0018][0019]
其中，φ(x)＝-ln p(x|z)；v(q)作为损失函数，后续步骤将对损失函数v(q)进行优化，求得损失函数的最小值。
[0020]
进一步的，步骤3具体包括如下步骤：
[0021]
步骤s301：将无人机的状态/数据的联合似然函数带到损失函数v(q)中
[0022]
无人机的状态/数据的联合似然函数如下：
[0023]
[0024]
其中，k表示子集中数据数量，k表示第k组数据，const表示常数值将此公式替换到损失函数v(q)中的φ(x)；
[0025]
步骤s302：利用牛顿法对损失函数v(q)进行更新；
[0026]
对损失函数v(q)关于高斯参数μ和∑-1
的导数如下：
[0027][0028][0029]
损失函数v(q)的二阶导如下：
[0030][0031]
比较后两式可得
[0032][0033]
v(q)的泰勒级数展开式如下：
[0034][0035]
其中：δμ＝μ
(i 1)-μ(i)，δ∑-1
＝(∑-1
)
(i 1)-(∑-1
)i，i是优化方案的迭代指数；
[0036]
为了将损失函数v(q)最小化，需要选择合适的δμ以及δ∑-1
；
[0037]
对于协方差∑-1
，令
[0038]
得到：
[0039][0040]
使用上式中的迭代方案后，逆协方差变化量为：
[0041][0042]
使用最大后验估计方法中的高斯非线性批量状态估计框架以及牛顿更新法求均值μ；
[0043]
为了确定最小值，为δμ取导数并令其为零，使得δμ的线性方程组转化为如下形式：
[0044][0045]
将δμ和δ∑-1
的优化方案代入损失函数的损耗近似中，得到：
[0046][0047]
当δμ和δ∑-1
不都为零，且关于μ和∑-1
的导数不都为零时，v(q)的损耗将在迭代过程中减少。
[0048]
进一步的，步骤3具体还包括如下步骤：
[0049]
步骤s303：通过自然梯度下降法更新δμ和δ∑-1
[0050]
自然梯度下降法利用信息的效率比常规的梯度下降法更高；vec(
·
)运算符将变分参数集中到单个列向量α中，将矩阵转化为向量：
[0051][0052][0053]
上式是损失函数v(q)关于α的梯度；将自然梯度下降法更新方案定义为：
[0054][0055]
γa是变分参数α的fisher信息矩阵；结合上述公式，得到：
[0056][0057]
是克罗内克内积；将单个更新提取：
[0058][0059][0060]
根据得到更新：
[0061][0062][0063]
最终，将迭代更新法案表示为：
[0064]
μ
(i 1)
←
μ(i) αbδμ
[0065]
∑-1(i 1)
←
∑-1(i)
αbδ∑-1
[0066]
其中，α＝0.95，b＝0，1，2，
…
一直增加，直到损失函数相比于前一次迭代的值下降为止。
[0067]
本发明的有益效果是：本发明能够对无人机的飞行系统采用自然梯度下降法对高斯变分推断进行改进，实现了对非线性系统进行批量状态估计。引入高斯变分推断况下的状态估计，相比较最大后验估计方法具有更高的精确性，能够更准确的估计无人机的飞行轨迹。
附图说明
[0068]
图1是map方法与gvi方法估计出的轨迹与真实值比较结果图。
[0069]
图2是无人机轨迹估计实验中map与gvi的损失函数值以及迭代次数图。
[0070]
具体实施方法
[0071]
为了使本发明实施例的目的、技术方案和优点更加明确，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0072]
针对无人机轨迹的高斯变分推断精确学习方法，首先必须要介绍本算法涉及的主要的算法：
[0073]
1、高斯变分推断
[0074]
变分推断方法希望将真实贝叶斯后验p(x|z)和后验的估计q(x)之间的kl散度最小化。x是希望估计的状态，z是已知观测到的状态。高斯变分推断则是假定q(x)中的所有状态x满足高斯分布，则q(x)可表述为如下形式：
[0075][0076]
散度的不对称性，首先需要选择使用kl(p||q)或者kl(q||p)。kl(p||q)前向kl散度表示如下：
[0077][0078]
kl(p||q)后向kl散度表示如下：
[0079][0080]
其中x∈rn，是试图从数据z∈rd推断的不可观测到的状态，e|
·
|代表期望。选择kl(q||p)的原因是式中的期望是后验的估计值q(x)的期望，利用该性质为q(x)设计的迭代方案能够更加有效地近似后验。
[0081]
对kl(q||p)表达式展开得：
[0082][0083]
是高斯分布的信息熵，ln p(z)是常数。将kl散度转化为损失函数，定义如下：
[0084][0085]
其中，φ(x)＝-ln p(x|z)。
[0086]
实施例1
[0087]
本实施例提出了一种无人机轨迹的高斯变分推断精确学习方法，实施步骤如下：
[0088]
步骤s1、根据无人机在实际飞行中的空气动力学原理以及运动学原理构建无人机
的数学模型。在实际飞行中，无人机的运动模型具有连续性以及非线性的特点，对其进行轨迹估计的前提是对无人机的运动参数进行状态估计即对无人机的位置、速度以及对应的加速度进行状态估计。
[0089]
为提升状态估计与轨迹估计的精度，首先对无人机进行建模，假设无人机模型满足如下条件：
[0090]
①
将无人机视作刚体，忽略旋翼旋叶的振动以及产生的形变；
[0091]
②
将无人机机身视作绝对对称，质量均匀的物体，同时重心与坐标原点重合；
[0092]
③
忽略地面效应作用在无人机上所产生的影响；
[0093]
④
忽略地球的自转，重力加速度保持恒定，将地面坐标系视作惯性坐标系。
[0094]
在惯性坐标系oxy中，假设无人机位于空间内任意一点，坐标为ξ＝[x，y，z]，欧拉角为其动力学方程式为：
[0095][0096]
式中，v表示无人机速度；n
x
，nz分别表示无人机切向过载以及法向过载；g表示重力加速度，取9.8m/s2。
[0097]
无人机结构参数与力学常量如表1所示，包括无人机的质量、轴半径、x轴、y轴以及z轴的转动惯量。
[0098]
表1 无人机结构参数与力学常量
[0099][0100]
为了评估gvi方法在轨迹估计方面的性能表现，将无人机作为研究对象对其进行轨迹预测。本实验将无人机飞行高度进行锁定，可以看作定高飞行。
[0101]
在轨迹估计实验中，需要测量无人机的位置信息以作为真实值来与估计值进行对比，主要包括无人机的位置、速度、角度、角速度。
[0102]
整个数据集包含10,000个时间步长，将其分为5个子集，每个子集包含2,000个时
间步长。
[0103]
需要估计的状态为：
[0104][0105]
其中，xk是无人机状态，对5个子集都有k＝2000.xk,yk,ψk分别是无人机在x轴，y轴的位置以及无人机的偏航角；分别是对应的速度以及角速度。
[0106]
对于无人机状态的线性先验因子，有以下形式：
[0107][0108]
其中
[0109][0110][0111]
qc＝diag(q
c，1
，q
c，2
，q
c，3
)
[0112]
t是离散时间采样周期，q
c，i
是功率谱密度，是各个量对应的方差。
[0113]
无人机非线性里程因子是：
[0114][0115]
其中：
[0116][0117][0118]
υk包括无人机的前向速度、侧向速度以及偏航角速度；是测量噪声方差。
[0119]
状态/数据的联合似然函数：
[0120][0121]
步骤s2：基于高斯变分推断，通过kl散度将真实贝叶斯后验和近似后验联系起来，
将状态的估计问题转化为后验的近似问题。再通过定义损失函数，将kl散度的计算转换为对损失函数的优化问题。
[0122]
近似后验和真实贝叶斯后验的kl散度展开式如上述高斯变分推断推出结果如下：
[0123][0124]
是高斯分布的信息熵，ln p(z)是常数，将kl散度转化为损失函数，定义如下：
[0125][0126]
其中，φ(x)＝-ln p(x|z)。v(q)作为损失函数，后续步骤将对v(q)进行优化，求得损失函数的最小值。
[0127]
步骤s3：采用牛顿法更新损失函数以及采用自然梯度下降法对均值和协方差逆矩阵进行迭代优化，使损失函数v(q)降到最低以达到有效近似后验的目的。
[0128]
步骤s301：将无人机的状态/数据的联合似然函数带到损失函数v(q)中
[0129]
无人机的状态/数据的联合似然函数带入到v(q)的损失函数中的φ(x)中，令。
[0130]
步骤s302：利用牛顿法对损失函数v(q)进行更新。
[0131]
对损失函数v(q)关于高斯参数μ和∑-1
的导数如下：
[0132][0133][0134]
损失函数v(q)的二阶导如下：
[0135][0136]
比较后两式可得
[0137][0138]
v(q)的泰勒级数展开式如下：
[0139][0140]
其中：δμ＝μ
(i 1)-μ(i)，δ∑-1
＝(∑-1
)
(i 1)-(∑-1
)i，i是优化方案的迭代指数。为了将损失函数v(q)最小化，需要选择合适的δμ以及δ∑-1
。
[0141]
对于协方差∑-1
，令
[0142][0143]
得到
[0144]
[0145]
联立公式，使用公式9中的迭代方案后，逆协方差变化量为：
[0146][0147]
使用最大后验估计方法中的高斯非线性批量状态估计框架以及牛顿更新法求均值μ。为了确定最小值，为δμ取导数并令其为零，使得δμ的线性方程组转化为如下形式：
[0148][0149]
将δμ和δ∑-1
的优化方案代入损失函数的损耗近似中，得到：
[0150][0151]
将δμ和δ∑-1
不都为零，且关于μ和∑-1
的导数不都为零时，v(q)的损耗将在迭代过程中减少。
[0152]
步骤s303：通过自然梯度下降法更新δμ和δ∑-1
[0153]
自然梯度下降法利用信息的效率比常规的梯度下降法更高。vec(
·
)运算符将变分参数集中到单个列向量α中，将矩阵转化为向量：
[0154][0155][0156]
上式是损失函数v(q)关于α的梯度。将自然梯度下降法更新方案定义为：
[0157][0158]
γa是变分参数α的fisher信息矩阵。结合上述公式，得到：
[0159][0160]
是克罗内克内积。将单个更新提取：
[0161][0162][0163]
根据得到更新：
[0164]
[0165][0166]
最终，将迭代更新法案表示为：
[0167]
μ
(i 1)
-μ(i) αbδμ
[0168]
∑-1(i 1)
←
∑-1(i)
αbδ∑-1
其中，α＝0.95，b＝0，1，2，
…
一直增加，直到损失函数相比于前一次迭代的值下降为止。
[0169]
图1表示了实验中map方法与gvi方法估计的轨迹与真实轨迹，实线为真实轨迹，虚线为map方法估计出的轨迹，点线为gvi方法估计出的轨迹。gvi方法估计的轨迹更加接近真实轨迹，map方法估计的轨迹偏离真实轨迹更多。在真实轨迹出现较大弯曲处，map方法估计的轨迹偏离真实值更多。
[0170]
图2中实线箱线与虚线箱线分别为无人机轨迹估计实验中map方法与gvi方法的损失函数值。图2中map方法与gvi方法损失函数上边缘值分别为5.923
×
103和5.918
×
103，上四分位值分别为5.913
×
103和5.908
×
103，中位数分别为5.903
×
103和5.900
×
103，下四分位值分别为5.882
×
103和5.880
×
103，下边缘值分别为5.875
×
103和5.871
×
103。对比数据得到gvi方法的损失函数值小于map方法，gvi方法能使损失函数收敛于更小的值。
[0171]
需要指出的是，以上，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化和替换，都应涵盖在本发明的保护范围内。因此，本发明的保护范围应以权利要求的保护范围为准。