SLAM学习笔记(五)——后端优化(一)

SLAM学习笔记(五)——后端优化(一)

先说明一下,书里面第八讲视觉里程计二,我写的文档丢失了,暂时先不去管他了,先继续往后看吧,之前是准备先做一个VO框架的,但是仔细思考了一下,准备先来实现一下这边,然后再去做那个VO框架。

我个人认为这里在我未来的项目中,不会是我做的重点,所以我在这块不做太过详细的学习。

EKF后端

我们知道,实际上,仅有前端的SLAM就可以支持简单的定位、建图的工作了,但是前端的结果往往是粗糙的,后端的存在就是为了继续优化前端的结果。它的本质是一种状态估计问题,主流的后端可以分为以下两类:

  • 渐进式: 保持当前状态估计, 加入新信息时,更新已有的估计
  • 批量式: 给定一定规模的数据, 计算该数据下的最优估计

接下来回到我们看了无数遍的方程:

image-20200830224611941

这里的x是机器人的位姿,z是观测数据,w与v是噪声,我们如果做一个不太严谨的定义:

image-20200830225535458

这里的 x k x_k xk代表k时刻所有的未知量,那么我们就可以重写之前的方程:

那么这个时候,我们就可以用数学语言来表示我们的问题了:

1

我们要做的事情实质上是:已知初始时刻的状态,以及中间过程的观测数据与误差(这里的误差或者说是噪声是可以人为假设出来的,比如高斯):

这里就是把原来的拆成了一个似然和一个先验的乘积,这个先验非常眼熟,上面的方程中 z k z_k zk的计算方法蕴含了其计算形式,后面的先验可以继续展开:

这个可以拆成两部分:

  • k时刻是如何受到之前状态的影响的
  • k-1时刻的状态估计

在这里时产生了一定的分歧:

  • 可以假设k时刻的状态只和k-1时刻有关
  • 假设k时刻与之前的所有状态有关

如果按照第一种方法来,那么第一项可以被写为:

第二项可以被写为:

这里的第一项表示:已知上一时刻的状态时与当前误差时,产生当前情况的概率,也就是我们所说的当前状态只与上一时刻状态相关。

这里的第二项与之前的差别在于将k时刻的误差去除了,因为想要预估k-1时刻的状态,无需知道所有时刻的误差,只需要知道到k-1时刻即可。

在开始下一步之前,我们先来明确一个问题,啥是卡尔曼滤波?

卡尔曼滤波(KF)

假设我们的运动方程和观测方程可以用线性方程来描述:

image-20200831112404112

这里的 w k , v k w_k, v_k wk,vk代表的是噪声,同时我们也假设这两个噪声服从高斯分布:

image-20200831112440905

首先我们通过运动方程来确定 x k x_k xk的先验分布,后面的高斯是通过上面的线性方程得到的:

image-20200831114803713

这个式子显示了如何通过上一时刻的变换来计算下一时刻的状态分布。这个分布也就是先验。

先验和后验(补充一下)
  • 可以说先验就是不基于经验的
  • 与之相对的是后验是基于经验的,后验需考虑当前的情况来对估计进行调整(贝叶斯)

那么我们这里的先验就记作:

image-20200831115701029

这里我们找到了每个时刻的先验的状态量与协方差与上一时刻的状态与协方差的关系,这个过程实质上就是直接将上面高斯的均值与方差拿出来了。

接下来,我们可以计算在某个状态下会产生怎么样的观测数据(这里是通过上面线性方程的第二个式子得来的)

image-20200831115954342

这个时候我们回看一下上一小节中我们想要求解的公式:

image-20200831120246627

我们发现,我们已经将两项的分布都求出来了,接下来我们对他们两个相乘就可以求出我们目标的概率,但是在此之前,我们需要假设我们乘出来的结果也是一个高斯分布,那么我们就得到了下面的结果:

image-20200831120013960

等式后面的第二项的均值与方差来自于先验,而第一项来自于我们刚才的推导,我们这里先来复习一下高斯(这里开始参考这篇博客):
f ( x ) = 1 2 π σ e − ( x − μ ) 2 2 σ 2 f(x) = \frac{1}{\sqrt{2π}σ}e^{-\frac{(x-μ)^2}{2σ^{2}}} f(x)=2π σ1e2σ2(xμ)2
上面我们看到的是一个普通一元的高斯,其中 μ \mu μ是均值, σ σ σ是方差。显然,我们的变量(状态变量)是非常多的变量构成的,也就是说这里的高斯不再是简单的一元高斯,而是多元高斯,那么多元高斯是怎么写的呢?我们假设所有的变量互不相关,就可以有下面的关系:
f ( x ) = p ( x 1 , x 2 . . . . x n ) = p ( x 1 ) p ( x 2 ) . . . . p ( x n ) = 1 ( 2 π ) n σ 1 σ 2 ⋯ σ n e − ( x 1 − μ 1 ) 2 2 σ 1 2 − ( x 2 − μ 2 ) 2 2 σ 2 2 ⋯ − ( x n − μ n ) 2 2 σ n 2 f(x) = p(x_{1},x_{2}....x_{n}) = p(x_{1})p(x_{2})....p(x_{n}) = \frac{1}{(\sqrt{2π})^nσ_{1}σ_{2}\cdotsσ_{n}}e^{-\frac{(x_{1}-μ_{1})^2}{2σ_{1}^2}-\frac{(x_{2}-μ_{2})^2}{2σ_{2}^2}\cdots-\frac{(x_{n}-μ_{n})^2}{2σ_{n}^2}} f(x)=p(x1,x2....xn)=p(x1)p(x2)....p(xn)=(2π )nσ1σ2σn1e2σ12(x1μ1)22σ22(x2μ2)22σn2(xnμn)2
假如我们有: z 2 = ( x 1 − μ 1 ) 2 σ 1 2 + ( x 2 − μ 2 ) 2 σ 2 2 ⋯ + ( x n − μ n ) 2 σ n 2 z^{2} = \frac{(x_{1}-μ_{1})^2}{σ_{1}^2}+\frac{(x_{2}-μ_{2})^2}{σ_{2}^2}\cdots+\frac{(x_{n}-μ_{n})^2}{σ_{n}^2} z2=σ12(x1μ1)2+σ22(x2μ2)2+σn2(xnμn)2,那么我们能够把原来的式子简写成:
f ( z ) = 1 ( 2 π ) n σ z e − z 2 2 f(z) = \frac{1}{(\sqrt{2π})^nσ_{z}}e^{-\frac{z^2}{2}} f(z)=(2π )nσz1e2z2
我们可以将 z 2 z^2 z2展开为矩阵的形式:
z 2 = z T z = [ x 1 − μ 1 , x 2 − μ 2 , ⋯   , x n − μ n ] [ 1 σ 1 2 0 ⋯ 0 0 1 σ 2 2 ⋯ 0 ⋮ ⋯ ⋯ ⋮ 0 0 ⋯ 1 σ n 2 ] [ x 1 − μ 1 , x 2 − μ 2 , ⋯   , x n − μ n ] T z^2 = z^\mathrm{T}z = \left[ \begin{matrix} x_{1} - μ_{1}, x_{2} - μ_{2}, \cdots,x_{n} - μ_{n}\end{matrix}\right] \left[ \begin{matrix} \frac{1}{σ_{1}^2}&0&\cdots&0\\ 0&\frac{1}{σ_{2}^2}&\cdots&0\\ \vdots&\cdots&\cdots&\vdots\\ 0&0&\cdots&\frac{1}{σ_{n}^2} \end{matrix}\right]\left[ \begin{matrix} x_{1} - μ_{1}, x_{2} - μ_{2}, \cdots,x_{n} - μ_{n}\end{matrix}\right]^\mathrm{T} z2=zTz=[x1μ1,x2μ2,,xnμn]σ121000σ221000σn21[x1μ1,x2μ2,,xnμn]T
我们可以把左边的矩阵转换成两个矩阵差的形式,即: x − μ x = [ x 1 − μ 1 , x 2 − μ 2 , ⋯   , x n − μ n ] T x - μ_{x} = \left[ \begin{matrix} x_{1} - μ_{1}, x_{2} - μ_{2}, \cdots,x_{n} - μ_{n}\end{matrix}\right]^\mathrm{T} xμx=[x1μ1,x2μ2,,xnμn]T

再把中间的协方差矩阵做一个符号的约定:
∑ = [ σ 1 2 0 ⋯ 0 0 σ 2 2 ⋯ 0 ⋮ ⋯ ⋯ ⋮ 0 0 ⋯ σ n 2 ] ∑_{}^{} = \left[ \begin{matrix} σ_{1}^2&0&\cdots&0\\ 0&σ_{2}^2&\cdots&0\\ \vdots&\cdots&\cdots&\vdots\\ 0&0&\cdots&σ_{n}^2 \end{matrix}\right] =σ12000σ22000σn2
因为他是一个对角矩阵,所以他的逆就是全都取倒数就可以了:
( ( ∑ ) − 1 = [ 1 σ 1 2 0 ⋯ 0 0 1 σ 2 2 ⋯ 0 ⋮ ⋯ ⋯ ⋮ 0 0 ⋯ 1 σ n 2 ] ( (∑_{}^{})^{-1} = \left[ \begin{matrix} \frac{1}{σ_{1}^2}&0&\cdots&0\\ 0&\frac{1}{σ_{2}^2}&\cdots&0\\ \vdots&\cdots&\cdots&\vdots\\ 0&0&\cdots&\frac{1}{σ_{n}^2} \end{matrix}\right] (()1=σ121000σ221000σn21
当然了,我们在分母项里面需要做一个累乘,所以我们还是做下面的约定:
σ z = ∣ ∑ ∣ 1 2 = σ 1 σ 2 . . . . . σ n σ_{z}= \left|∑_{}^{}\right|^\frac{1}{2} =σ_{1}σ_{2}.....σ_{n} σz=21=σ1σ2.....σn
最后,我们把说的这些东西综合起来, z 2 z^2 z2就可以写成下面的形式:
z T z = ( x − μ x ) T ∑ − 1 ( x − μ x ) z^\mathrm{T}z = (x - μ_{x})^\mathrm{T} \sum_{}{}^{-1} (x - μ_{x}) zTz=(xμx)T1(xμx)
最后,我们来看部分展开的多元高斯:
f ( z ) = 1 ( 2 π ) n σ z e − z 2 2 = 1 ( 2 π ) n ∣ ∑ ∣ 1 2 e − ( x   −   μ x ) T   ( ∑ ) − 1   ( x   −   μ x ) 2 f(z) = \frac{1}{(\sqrt{2π})^nσ_{z}}e^{-\frac{z^2}{2}} = \frac{1}{(\sqrt{2π})^{n}\left|∑_{}^{}\right|^\frac{1}{2}}e^{-\frac{ (x\ -\ μ_{x})^\mathrm{T}\ (\sum_{}{})^{-1}\ (x\ -\ μ_{x})}{2}} f(z)=(2π )nσz1e2z2=(2π )n211e2(x  μx)T ()1 (x  μx)
由于这里的两侧都是高斯,所以我们就只考虑指数部分就可以了,我们就有:

image-20200831153731206

可以将两侧展开,分别比较x的二次和一次系数,对于二次系数,有:

image-20200831153828975

这里的二次项就是展开后含有 x k 2 x_k^2 xk2的项,这里为了简化计算,我们定义一个中间变量:

image-20200831154546755

接下来,我们将上面式子的左边通过乘上他的逆的形式转换为一个单位矩阵,就可以得到下面的式子:

image-20200831154754733

这里转换的过程中,同时左乘了一个 P k P_k Pk帽,然后又使用刚才定义的 K K K,简化了这个式子。然后通过简单变换,可以写出:

image-20200831154940187

这里是二次项,接下来我们来比较一次项:

image-20200831155043748

经过一波化简,(先取转置,然后再消元)就可以得到下面的形式:

image-20200831161548477

最后,两侧同时乘以 P k P_k Pk帽,并且带入上面的中间变量 K K K P k P_k Pk帽的表达就可以得到下面的式子:

image-20200831161832674

于是呢,我们又找到了 x k x_k xk帽与其他之间的关系,同时我们根据二次项也可以知道 P k P_k Pk帽与其他之间的关系,于是我们就可以根据之前的信息对当前的信息进行更新,这一步我们叫做预测

image-20200831163451457

那么预测完成后,我们要计算非常重要的K:

image-20200831163549235

最终,计算后验分布的概率:

image-20200831163614442

这就是卡尔曼滤波的全过程,说实话我感觉我顺一遍还是有点小问题,难顶。

回到EKF后端

我们上面说的是KF,也就是标准的卡尔曼滤波,他只能运用到线性系统中,那么我们需要用EKF来做实际应用,接下来我们看具体是怎么做的:

首先,我们用泰勒做一个一阶展开:

image-20200831164541206

我们将这里的偏导数记作:

image-20200831164606971

对于观测方程也是如此:

image-20200831164635364

一样的道理,我们用符号来表达一下:

image-20200831164653903

这里与上面的卡尔曼滤波的过程相似,最终得到的结果是:

image-20200831165036774

我们目前对卡尔曼滤波仅作简单掌握吧,我们可以看到,卡尔曼滤波存在一个较大的局限:我们展开的过程是做的一阶展开,实质上是在用线性函数来近似非线性函数,这样的做法有非常大的局限性,并且不够精准,接下来我们来看下一个方法——BA与图优化。

BA与图优化

在看这个之前,我们需要复习一下我们的相机模型

  • 首先从世界坐标到相机坐标:image-20200831172305492
  • 从相机坐标到归一化坐标: image-20200831172341318
  • 进行畸变处理: image-20200831172411110
  • 根据内参模型,来计算图像坐标:image-20200831172433324

我们将上面四个步骤使用一个函数来概括:
z = h ( x , y ) z = h(x,y) z=h(x,y)
我们可以用流程图来概括 h ( x , y ) h(x,y) h(x,y)的作用:

image-20200831173018908

这里的x指代当前相机的位姿,也就是外参 R , t R,t R,t,它对应的李代数是 ξ \xi ξ。路标y是这里的三维点,而观测数据是: z Δ = [ u s , v s ] T z^{\Delta} = \left[ u_s, v_s\right]^T zΔ=[us,vs]T,我们可以构建残差:
e = z − h ( x , y ) e = z - h(x,y) e=zh(x,y)
他的实质含义与我们以前经常说的一样:实际观测到的数据,与我们预测出的数据的差异。接下来我们将总体(也就是不同时刻)的残差一并考虑进来,即可构建代价函数:

image-20200831173535197

我们对这个目标方程来做一个优化就是传说中的BA了。

BA求解

我们可以判断,上面的目标函数并不是一个线性函数,所以我们需要使用一些非线性优化的方法来优化他,具体来说就是通过从某一个初值点开始,不断地寻找下降方向 Δ x \Delta x Δx来找到目标函数的最优解。对于整体的目标函数,我们必须把自变量定义成所有待优化的变量:
x = [ ξ 1 , . . . , ξ m , p 1 , . . . , p n ] T x = \left[ \xi_1, ..., \xi_m, p_1, ..., p_n \right]^T x=[ξ1,...,ξm,p1,...,pn]T
相应的,当我们对方程中的x增加一定的增量时,其误差变化如下:

image-20200831202853038

我们可以将原来的x拆分成两个部分:
x ξ = [ ξ 1 , ξ 2 , . . . , ξ m ] , x p = [ p 1 , . . . , p n ] x_{\xi} = \left[ \xi_1, \xi_2 , .. ., \xi_m \right], x_{p} = [p_1, ..., p_n] xξ=[ξ1,ξ2,...,ξm],xp=[p1,...,pn]
于是这个式子就可以化简成这个形式:

image-20200831203357255

我们面对的是这样的一个增量线性方程:
H Δ x = g H\Delta x = g HΔx=g
其中,由于我们吧变量分为了位姿和空间点两种,所以雅可比矩阵也可以分块:
J = [ F E ] J = \left[ \begin{matrix} F & E \end{matrix}\right] J=[FE]
以高斯牛顿为例:
H = J T J H = J^TJ H=JTJ
这样就可以进行求解了,后面好像还讲了一些其他的,关于这个矩阵H的逆的求法的问题,但是由于篇幅太长,而且和我做的事情的关系不是很大的原因,我决定先跳过这一部分,日后再看。

小结

这节主要讲了两种后端优化方法:

  • 扩展卡尔曼滤波:实质上还是线性的卡尔曼滤波,只不过是用了一阶展开来做了线性近似,有一定的局限性,但是效果不错(据说)
  • BA:使用最优化方法来解的全局能量函数最小化问题,比较灵活,看作者的意思是,比较推荐用BA
  • 2
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值