Python写基于非线性优化的2D-SLAM系统（已开源）

Github代码地址:
https://github.com/liuzhenboo/2D-SLAM-By-Nonlinear-Optimization

1. SLAM问题概率模型

1.1 最大后验到最小二乘

SLAM问题其实就是一个状态估计的问题，就是要根据一系列观测来推测状态量；一般情况下，我们把SLAM问题建立在概率论的框架下。

说白了，SLAM就是要是解决这样一类后验概率问题：

$$p(x|z)$$

其中$x$是系统当前状态量，$z$是与状态量相关的观测量；在求出这个条件概率之后，把观测值$Z$带入，就可以获得$p(x|z=Z)$，这也就是我们需要的分布。

根据概率公式展开：

$$p(x|z=Z)=\frac{p(z=Z|x)p(x)}{p(z=Z)}$$

其中分母是常量，当作是归一化因子$\eta$即可，所以也可写为：

$$p(x|z=Z)=\eta p(z=Z|x)p(x)$$

$p(z|x)$表示在$x$已知的情况下，$z$的概率分布；在求出$p(z|x)$之后，把$z=Z$带入，得到$p(z=Z|x)$，它也被称为似然项，它表示在$x$取不同值的情况下，$z=Z$的概率。

$p(x)$也被称为先验，也就是在观测之前$x$服从什么分布。这可以通过其他信息源得到，比如GPS，惯导等。如果事前不知道任何先验分布信息，那么可以将其设为1，表示不相信先验信息，只根据系统使用的量测来估计。

$$x_{MAP}=\underset{x}{\mathrm{argmax}}p(z=Z|x)p(x)$$

因为每次观测相互独立，所以上式可以写为：

$$x_{MAP}=\underset{x}{\mathrm{argmax}}\prod _{i}p(z_i=Z_i|x)p(x)$$

取负对数得到：

$$x_{MAP}=\underset{x}{\mathrm{argmin}}\left[-\sum _{i}logp(z_i=Z_i|x)-logp(x)\right]$$
如果假设观测$p(z_i|x)$服从高斯分布$N(h_i(x)，{\sum }_i)$；先验$p(x)$服从$N(x_p，{\sum }_x)$那么有：

$$x_{MAP}=\underset{x}{\mathrm{argmin}}\sum _i||z_i-h_i(x)|{|}_{{\sum }_i}^2+||x-x_p|{|}_{{\sum }_x}^2$$

如果没有先验知识，那么可以写为：

$$x_{MAP}=\underset{x}{\mathrm{argmin}}\sum _i||z_i-h_i(x)|{|}_{{\sum }_i}^2$$
显而易见，在SLAM的过程中，不断有新的残差项加入到系统，上述最小乘问题不断扩大。

1.2 非线性优化求解

设当前时刻状态向量为x，维度为n，我们将残差表示成向量的形式：

$$e(x)=\left[\begin{array}{c}{e}_1(x)\\ ...\\ e_m(x)\end{array}\right]$$
这里我们使用马氏范数，损失函数定义为：

$$E(x)=\frac{1}{2}||e(x)|{|}_{{\Sigma }_e}^2=\frac{1}{2}{e(x)}^T{\Sigma }_e^{-1}e(x)$$

记$J_i(x)=\frac{\partial e_i(x)}{\partial x}$，那么可以得到：

$$\frac{\partial e(x)}{\partial x}=J=\left[\begin{array}{c}{J}_1(x)\\ ...\\ J_m(x)\end{array}\right]$$

$$J_i(x)=\left[\begin{array}{c}\frac{\partial e_i(x)}{\partial x_1}...\frac{\partial e_i(x)}{\partial x_n}\end{array}\right]$$

对残差函数泰勒展开，得到一阶近似：

$$e(x+\delta x)\approx e(x)+J\delta x$$
那么带入损失函数可以近似得到：

$$\begin{gathered} E(x+\delta x)=\frac{1}{2}{e(x+\delta x)}^T{\Sigma }_e^{-1}e(x+\delta x) \\ \approx \frac{1}{2}{e(x)}^T{\Sigma }_e^{-1}e(x)+{\delta x}^T{J}^T{\Sigma }_e^{-1}e(x)+\frac{1}{2}{\delta x}^T{J}^T{\Sigma }_e^{-1}J\delta x \end{gathered}$$

所以损失函数就转换为关于$\delta x$的二次函数，并且如果$J^T{\Sigma }_e^{-1}J$正定，那么损失函数有最小值。

对上式关于$\delta x$求一阶导数，并另其为零，得：

$$J^T{\Sigma }_e^{-1}J\delta x=-J^T{\Sigma }_e^{-1}e(x)$$

也常常写为

$$H\delta x=b$$

这样就可以得到增量$\delta x$，这样一直迭代，就可以不断优化当前解，直到误差小于阈值。这种方法被称为高斯牛顿法。

不过在实际中，求解上述方程需要$H$矩阵可逆，而我在编程中，经常遇到$H$不可逆的情况，所以实际代码中我使用了Levenberg-Marquardt法。其是对高斯牛顿法进行了改进，增加了阻尼因子一项，如下所示：

( H + μ I ) δ x = b (H+\mu I)\delta x = b