IMU预积分,耳熟能详,但是一直没有仔细了解。这里详细介绍一下。
## 运动模型
首先做一个简单规定,因为IMU预积分离不开IMU坐标系与世界坐标系的转换,因此我们将世界坐标下的值标记为,将IMU坐标系下的值标记为。而表示了从世界坐标到IMU坐标的转换。
一般来说,IMU包含了两部分——陀螺仪与加速计,因此IMU的raw
data有两类:
$$
\tilde{\boldsymbol{\omega}}^{b}(t)=\boldsymbol{\omega}^{b}(t)+\mathbf{b}_
{g}(t)+\boldsymbol{\eta}_ {g}(t)\\\tilde{\mathbf{a}}^{b}(t)=\mathbf{R}_
{b(t)}^{w T}\left(\mathbf{a}^{w}(t)-\mathbf{g}^{w}\right)+\mathbf{b}_
{a}(t)+\boldsymbol{\eta}_ {a}(t)
$$
其中与分别是读取到的角速度与加速度,也就是观察值,属于IMU坐标系下的读数,且包含噪声与偏差(分别是加速度accelerometer与陀螺仪gyro的噪声与偏差)。同时,IMU的加速度测量的是IMU坐标下的一个力,但是由于测量原理,它测出来的加速度去除了重力,例如一个三轴加速度计放在水平面上,那么它的输出将是铅垂线反向的;如果一个加速度计做自由落体运动,那么它的输出会是0。而我们关注的物体真正的加速度是世界坐标系下,当然也包含了重力。所以真实加速度需要去除重力加速度再转换到IMU坐标系下才是IMU的加速度读数。
同时,根据运动模型,我们可以得到旋转矩阵,速度以及位置的微分模型如下:
有意思的一件事是,虽然旋转矩阵是一个从世界坐标到IMU坐标的转换(也就世界坐标系下姿态的逆),但是微分里的角速度是IMU坐标系的,这意味着我们可以直接使用IMU去噪后的角速度读数。使用欧拉积分(?),可以得到下一刻的状态(离散形式):
上式中速度和位置都很容易理解,对于旋转矩阵,很显然上面给出的结果也是很符合直观的,就是将角速度与间隔时间相乘后得到一个旋转向量(角速度乘上时间也就得到了绕各个轴旋转的角度,可以直接得到旋转向量,所以一个旋转向量其实可以理解成绕坐标轴旋转的角度?),再转化到上进行乘积。而一个严格的解释如下:
为了简化符号,规定:
$$
\mathbf{R}(t) \doteq \mathbf{R}_ {b(t)}^{w} ; \quad
\boldsymbol{\omega}(t) \doteq \boldsymbol{\omega}^{b}(t) ; \quad
\mathbf{a}(t)=\mathbf{a}^{b}(t) ; \\\quad \mathbf{v}(t) \doteq
\mathbf{v}^{w}(t) ; \quad \mathbf{p}(t) \doteq \mathbf{p}^{w}(t) ; \quad
\mathbf{g} \doteq \mathbf{g}^{w}
$$
如果将测量模型带入到运动模型中,可以得到:
首先要注意的是,上面的噪声也从连续形式变成了离散形式,而二者的区别就在于他们有不同的协方差:
如果假设采样频率一定,时间间隔为,那么每个时刻可以表示为离散的,因此,进一步简化可以将下一个时刻的状态写为:
预积分
根据之前的内容,如果有两个时刻,且我们知道时刻的机器人状态,以及时刻之间的IMU状态,那么我们可以直接得到时刻的机器人状态,通过对时刻间的IMU状态进行上述的积分:
其中,也就是时刻之间的时间间隔。可以看到,上述估计中相邻两个IMU数据之间的状态用先前的IMU状态来代替。后面的速度与位置的更新,都出现了新的变量,也就意味着当有变化,所有的都会跟着改变,所有的量都需重新计算。为了避免这样的情况,很直接的方法就是记录着时刻之间的变化值,而这个变化值不应该随着
的改变而改变。因此,需要引入预积分。预积分的目的自然是为了消除IMU积分式中的以及与其相关的,这样之后改变后,可以直接利用先前计算的预积分得到更新的结果。因此我们定义下面为预积分的内容:
预积分的这三个值,第一个容易理解以及第二个是容易理解的,就是在IMU积分项中把,,提取出。而对于位置的预积分,更复杂一点,证明如下:
令,有
有了预积分项,想要得到积分后的结果就非常容易了。预积分的意义是找出时刻之间IMU状态贡献的部分,即使机器人状态有变化,固定时刻之间IMU状态测量值不变,这个贡献也不应该改变,因此IMU预积分是可以当作两个机器人状态之间的一个约束。
预积分测量值以及测量噪声
这一部分是为了分离出噪声,直接使用测量值的预积分。首先假设时刻之间的bias是一定的,那么
(1)对于有:
这一步是把看作是一个微小量,则有。
从这一步继续往下推导会稍微复杂一点,先给出结果:
为了符号简便,令:
$$
\mathbf{J}_ {r}^{k}=\mathbf{J}_
{r}\left(\left(\tilde{\boldsymbol{\omega}}_ {k}-\mathbf{b}_
{i}^{g}\right) \Delta t\right)\\\Delta \tilde{\mathbf{R}}_ {i j}=\prod_
{k=i}^{j-1} \operatorname{Exp}\left(\left(\tilde{\boldsymbol{\omega}}_
{k}-\mathbf{b}_ {i}^{g}\right) \Delta
t\right)\\\operatorname{Exp}\left(-\delta \vec{\phi}_ {i
j}\right)=\prod_ {k=i}^{j-1} \operatorname{Exp}\left(-\Delta
\tilde{\mathbf{R}}_ {k+1,j}^{T} \cdot \mathbf{J}_ {r}^{k} \cdot
\boldsymbol{\eta}_ {k}^{g d} \Delta t\right)
$$
则有:
其中就是旋转量的预积分测量值,后者为噪声。
下面我们证明。首先,我们利用的伴随性质可以得到:
这里我们对上式进行一个推广:
也就是,伴随性质是链式传播的,从最右侧传到最左侧,中间每一个量都会受到影响。为了符号简便,我们可以令,也就是:
则有:
证毕,需要注意的是。
(2)对于,将(1)中的结论带入有:
令:
则有:
其中前一项为测量项,而后一项为噪声值。
(3)对于,带入前面的结果,可得:
令:
则
前者为预积分测量值,后者为噪声值。
最终,将分离后的噪声与测量值带入理想状态的预积分形式下,可以得到:
类似于测量值等于理想值+噪声的形式。
在实现的时候需要注意,会用到的是,也就是不会用到,因此实现的时候需要注意计算顺序:。
噪声分布
令预积分的噪声为:
我们希望它满足高斯分布,也就是。这里是三种噪声的线性组合,因此我们来分别分析。
(1)首先是旋转的噪声分析。由之前的推导可以知道:
对等式两边同时取对数,可以得到:
令,由于是小量,因此也是小量,因此有(注意这里是雅可比矩阵,而不是之前的缩写)。由于本身是小量,因此任意是小量,因为有,我们可以得到
也就是
由于上式中都是已知量,而假设是零均值高斯噪声,因此也是零均值高斯噪声。
(2)由于我们分析得到了是零均值高斯噪声,根据表达式:
我们可以知道也是高斯分布的。
(3)类似的,根据的表达式:
可以知道也是高斯分布。但是,不是零均值。
噪声的递推形式
下面分析一下噪声的递推形式,,以及其协方差的递推形式,。
(1)旋转项:
(2)速度项:
(3)位置项:
综上所述,可以得到噪声预积分的递推形象如下(令):
其中:
而噪声的协方差也就有了下面的递推形式:
bias更新
目前为止做的计算都是假设了时刻之间的bias是不变化的。如果bias发生了更新,需要将预计分测量值整个重新计算(bias在这段时刻里依然是恒定的,但是bias的值可能会被优化等方式更新)。目前的做法是将这个过程线性化,来得到一阶近似结果。现在先规定一些符号,首先我们将旧的bias标记为,而新的bias是由旧的bias加上更新量得到的,即:
则有一阶近似如下(类似于性质):
如果做符号简化如下:
就可以写成:
下面我们需要推导的是式子中的各个偏导数。
(1)旋转项:
到了这一步后,使用的变形手法与之前分离测量值与噪声的时候类似,我们依然利用伴随性质来处理。先做一些符号简化,令:
$$
\mathbf{M}_
{k}=\operatorname{Exp}\left(\left(\tilde{\boldsymbol{\omega}}_
{k}-\overline{\mathbf{b}}_ {i}^{g}\right) \Delta
t\right)\\\operatorname{Exp}\left(\mathbf{d}_
{k}\right)=\operatorname{Exp}\left(-\mathbf{J}_ {r}^{k} \delta
\mathbf{b}_ {i}^{g} \Delta t\right)
$$
则有:
根据伴随性质,我们可以得到:
由于,则上式可以写为:
也就是:
令,由于很小,所以很小,可以利用性质以及,则有:
因此有:
所以可得:
即:
在实现时,可以进一步引出递推式:
(2)速度项():
上述中分别用到了一阶近似,以及,且忽略了高阶小项。所以有:
(3)位置项():
下面我们分别推导和部分。
的推导与速度项是非常类似的,也是忽略了高阶小项。将合起来,可以得到:
所以有:
残差
之所以需要残差,是因为IMU预积分的测量值实际上也是优化中的一个约束量,pose
graph中的边。而一般来说,预积分的残差是由计算值和测量值决定的,这里的计算值可能是通过视觉约束等等方法得到的。我们知道IMU预积分的理想值是如下:
其三项残差定义如下:
从定义上来看,残差也就是用来度量计算值与测量值的差别,这也是很合理的。从前面的推导可以知道,这个测量值是由噪声加上理想值得到的(实际上我们是从理想值中分离出了噪声)。
一般来说,在SLAM系统中,我们需要优化的是,不过在带有IMU的系统中,bias的影响往往不可忽视,因此一般也会加上对bias的优化。不过从残差的定义上来看,实际上优化的是bias的变化量。所以在IMU预积分中,完整的一个导航状态为:。
对各个状态的更新操作如下:
需要注意的是,由于我们优化的是bias的增量,因此这里的实际上是增量的增量,也就是通过增量的Jacobian来计算的。但是实际上增量的Jacobian与原来的Jacobian计算方法上是一样的,只不过它的初始状态是从0开始(实际上之前也运用过这样的方式进行优化过,每次优化的变量实际上是更新量)。
为何位置没有采用速度一样的方式去更新()?这个问题实际上如果是在上,也就是直接用transformation
matrix,得到的结果是:
也就是上述更新的操作,从物理意义上理解,就是先进行了旋转再进行了平移的变换。不过,毕竟我们一麻溜推导下来,用的都是上的操作,平移与旋转是分开操作的(平移的残差中已经包含了旋转的计算),我个人认为还是理论上更合理一些。在实际中,可能二者使用起来没什么太大的区别。
由于噪声不是的零均值正态分布,因此不同的协方差对应着不同约束项的权重(实际上权重是协方差的逆),这样给出来的结果是一个无偏估计。最终我们优化目标就是能使得加权残差的平方和最小。另外,在优化时候,实际上优化变量是由Jacobian矩阵来决定的(高斯牛顿的海森矩阵也是由Jacobian矩阵来近似的),因此我们的下一个任务就是求出Jacobian矩阵。
Jacobian
接下来是最后一步,计算出Jacobian矩阵。我们将Jacobian矩阵分成三类:(1)”0“矩阵,也就是残差和某些变量是完全无关的;(2)线性的,也就是残差和某些变量是线性相关的,这样的Jacobian也可以由系数直接得到;(3)其他的,也就是残差和某些变量的关系较为复杂,需要进行相应的变形才能求得Jacobian矩阵。
的Jacobian
(1)“0”矩阵。由于中不含,因此有:
$$
\frac{\partial \mathbf{r}_ {\Delta \mathbf{R}_ {ij}}}{\partial
\mathbf{p}_ {i}}=\mathbf{0},\frac{\partial \mathbf{r}_ {\Delta
\mathbf{R}_ {ij}}}{\partial \mathbf{p}_ {j}}=\mathbf{0},\\\frac{\partial
\mathbf{r}_ {\Delta \mathbf{R}_ {i j}}}{\partial \mathbf{v}_
{i}}=\mathbf{0},\frac{\partial \mathbf{r}_ {\Delta \mathbf{R}_ {i
j}}}{\partial \mathbf{v}_ {j}}=\mathbf{0},\\\frac{\partial \mathbf{r}_
{\Delta \mathbf{R}_ {i j}}}{\partial \delta \mathbf{b}_
{i}^{a}}=\mathbf{0}
$$
(2)线性类:无。
(3)其他类:
- 下面求关于的Jacobian:
上面的推导中,用到的性质有:,伴随性质,。综上所述,有:
- 使用类似的方法,可以得到的Jacobian:
因此有:
- 关于的Jacobian(都是与残差项有关,而则是与残差项有关):
令,则:
上面的推导中,用到的性质有:,,伴随性质,(可以看出来推导的方法是很类似的)。综上所述,可以得到:
这里需要注意的这儿我们是对求导,而也是新的bias()的测量预积分值。至于为什么要这么做,也许之后会找到答案。
的Jacobian
(1)“0”类:由于中不包含,因此:
(2)线性类:由于与是线性关系,因此有:
(3)其他类:
- 关于的Jacobian
因此有:
- 关于的Jacobian
因此有:
上述推导用到的规则有:,,。因此有:
的Jacobian
(1)“0”类:由于中不包含
(2)线性类:由于与是线性关系,因此有:
(3)其他类:
- 关于的Jacobian:
因此可以得到:
- 关于的Jacobian:
因此有:
- 关于的Jacobian:
因此有:
- 关于的Jacobian:
上述推导用到的规则有:,,。因此有:
注意,由于位置的更新方式是,因此我们在推导的时候也是加上了这个来求Jacobian。如果换了的更新方式,实际上结果也一样根据上面的方法很容易求出来。
至此,全文将 Foster 的 IMU
预积分理论中的公式推导结束,感谢邱笑晨提供的PDF。
Note: 在Vins mono中,IMU的重力也会被优化。