1. 前言
本博客借鉴了崔华坤的《VINS论文推导及代码解析》和 VINS-Mono理论学习——IMU预积分 Pre-integration (Jacobian 协方差)的内容,因为确实写得太好了,然后有些地方加入自己一些理解。
VINS-MONO论文中的IV-B. IMU Pre-integration介绍了IMU预积分模型,Foster的倆篇论文对IMU预积分理论进行详细分析。
值得注意的是,本文使用了四元素中Hamilton(可参考《Quaternion kinematics for the error state Kalman filter》,右手系)的表示和JPL(可参考《Indirect Kalman Filter for 3D Attitude Estimation》,左手系)表示,因为Vins-mono里面也是用了俩种表示混用,俩者的具体区别可以看《Quaternion kinematics for the error state Kalman filter》的第17页,且论文中用到的JPL的右乘左乘定义还与《Indirect Kalman Filter for 3D Attitude Estimation》不同,这种问题估计让无数人心塞。/(ㄒoㄒ)/~~
为什么需要对IMU进行预积分?
传统捷联惯性导航算法,在已知 k 时刻下的IMU状态量(姿态、速度和位移)情况下,利用IMU测量的线加速度和角速度,通过积分预算得到k+1时刻下的状态量。
然而在非线性优化的VIO中,各个节点的状态量都是估计值,当算法对这些状态量优化时,每次调整都需要在它们之间重新积分,导致绝对位姿被优化时对状态量进行重复积分。IMU预积分的提出使得优化算法可对IMU的相对测量进行处理,使它与绝对位姿解耦,或者只要线性运算就可以进行矫正。
2. IMU模型
IMU测量值包括加速度计得到的 (测量值) 线加速度at^ 和陀螺仪得到的角加速度w^t 【论文式(1)】
a^t=at+bat+Rwtgw+na w^t=wt+bwt+nw其中t下标表示在IMU的体(body)坐标系下,at、wt分别表示IMU真实的线加速度和角速度,并受到加速度偏置(bias) bat、陀螺仪偏置bwt和附加噪声na、nnw的影响。计算得到的线加速度at^是重力加速度和物体加速度的合矢量。
假设附加噪声为高斯噪声:
na∼(0,σa2), na∼(0,σw2) 加速度计偏置和陀螺仪偏置被建模为随机游走,其导数为高斯分布:【论文式(2)】
b˙at=nba, b˙wt=nbw nba∼N(0,σba2), nbw∼N(0,σbw2)
3. 基于世界坐标系下的IMU运动模型
3.1 连续形式下的IMU运动模型
对于图像帧k和k+1, IMU body坐标系对应为bk和bk+1,位置、速度和姿态状态值PVQ(Pose、Velocity、Quaternion)可以根据[tk,tk+1]时间间隔内的IMU测量值,在世界坐标系下进行传递:【论文式(3)(4)】
pbk+1w=pbkw+vbkwΔtk+∫∫t∈[tk,tk+1](Rtw(a^t−bat−na)−gw)dt2 vbk+1w=vbkw+∫t∈[tk,tk+1](Rtw(a^t−bat−na)−gw)dt qbk+1w=qbkw⊗∫t∈[tk,tk+1]21qtbk⊗[(w^t−bwt−nw)0]dt=qbk⊗∫t∈[tk,tk+1]21Ω(w^t−bwt−nw)qtbkdt(1) Ω(w)=[−[w]×−wTw0], [w]×=⎣⎡0wz−wy−wz0wxwy−wx0⎦⎤其中 Δtk是[tk,tk+1]之间的时间间隔,Rtw为t时刻IMU body坐标系到世界坐标系的旋转矩阵,qtbk为用四元素表示的 t 时刻从IMU body坐标系到 k 时刻IMU body坐标系的旋转,这里的四元素实部在后,虚部在前, 为了与论文保持一致,即四元素JPL表达。这里的Ω(w)表示四元素右乘。
关于公式(1)的推导,这里首先引入四元素左乘右乘及导数定理:
这里与《Indirect Kalman Filter for 3D Attitude Estimation》中公式(10)不同,主要是因为左右手系不同,我们引入左乘和右乘符号如下:
qa⊗qb=R(qb)qa=⎣⎢⎢⎡sb−zbyb−xbzbsb−xb−yb−ybxbsb−zbxbybzbsb⎦⎥⎥⎤⎣⎢⎢⎡xayazasa⎦⎥⎥⎤=L(qa)qb=⎣⎢⎢⎡saza−ya−xa−zasaxa−yaya−xasa−zaxayazasa⎦⎥⎥⎤⎣⎢⎢⎡xbybzbsb⎦⎥⎥⎤
为了简化,令q=[x y z s]=[w s], 则有:
R(q)=Ω(w)+sI4×4=[−[w]×−wTw0]+sI4×4 L(q)=Ψ(w)+sI4×4=[[w]×−wTw0]+sI4×4
对于四元素的求导,我们定义qt为t时刻下的单位四元素(这里的四元素实部在后,虚部在前),w为qt确定的角速度,则关于qt的导数为: q˙t=21[−[w]×−wTw0]qt=21Ω(w)qt=21R([w0])qt=21qt⊗[w0]
当看到q˙=21qt⊗[0w] 说明该四元素实部在前,虚部在后(Hamilton表示),只是表示不同而已。
因此对于IMU连续形式下的旋转状态(用四元素表示)推导,我们有:
qbk+1w=qbkw⊗qbk+1bk=qbkw⊗∫t∈[tk,tk+1]qt˙bkdt=qbkw⊗∫t∈[tk,tk+1]21qtbk⊗[wtbk0]dt=qbkw⊗∫t∈[tk,tk+1]21qtbk⊗[wt^−bwt−nw0]dt==qbk⊗∫t∈[tk,tk+1]21Ω(w^t−bwt−nw)qtbkdt
3.2 离散形式下的IMU运动模型
3.2.1 欧拉法离散形式
使用欧拉法,即k+1时刻的位姿是用第k时刻的测量值a^bk, w^bk来计算的:
pbk+1w=pbkw+vbkwΔtk+21a^bkδt2 vbk+1w=vbkw+a^bkδt qbk+1w=qbkw⊗[121w^bkδt] 其中 a^bk=qbkw(abk−bak)−gww^bk=wbk−bwk
3.2.2 中值法离散形式
使用中值法,即k+1时刻的位姿是用俩个时刻k和k+1测量值a, w的平均值来计算的:
pbk+1w=pbkw+vbkwΔtk+21a^ˉtδt2 vbk+1w=vbkw+a^ˉtδt qbk+1w=qbkw⊗[121w^ˉtδt] 其中 a^ˉt=21[qbkw(abk−bak)−gw+qbk+1w(abk+1−bak+1)−gw]w^ˉt=21(wbk−bwk+wbk+1−bwk+1)=21(wbk+wbk+1)−bwk 假设在短时间内加速度计和陀螺仪的偏置不变,则有:bak=bak+1, bwk=bwk+1
4.IMU预积分 (基于第K帧IMU body坐标系下的运动模型)
通过公式(1)可以看到,IMU 的积分需要依赖与第 k 帧的 v 和 R (基于世界坐标系下的),当我们在后端进行非线性优化时,需要迭代更新第 k 帧的 v 和 R,这将导致我们需要根据每次迭代后的值重新进行积分,这将非常耗时。我们考虑将优化变量从第 k 帧到第 k+1 帧的 IMU 积分项中分离开来。
4.1 连续形式下的IMU运动模型
IMU预积分的思想就是将参考坐标系从世界坐标系w调整为第k帧的IMU body坐标系bk下,可通过在等式俩端同时乘以Rwbk得到:【论文式[5][6]】
Rwbkpbk+1w=Rwbk(pbkw+vbkwΔtk−21gwΔtk2)+αbk+1bk Rwbkvbk+1w=Rwbk(vbkw−gwΔtk)+βbk+1bk qwbk⊗qbk+1w=γbk+1bk(2)
其中 αbk+1bk=∫∫t∈[tk,tk+1]Rtbk(a^t−bat−na)dt2 βbk+1bk=∫t∈[tk,tk+1]Rtbk(a^t−bat−na)dt γbk+1=∫t∈[tk,tk+1]21Ω(w^t−bwt−nw)γtbkdt
此时的积分结果αbk+1bk、βbk+1bk、γbk+1bk可以理解为bk+1对bk的相对运动量,bk的状态并不会对其产生影响,因此将其作为非线性优化变量,可以避免状态的重复传递。
注意,这是在假设IMU偏置ba、bw已经确定的情况下,实际上偏置也是需要优化的变量,那么每次迭代时,ba、bw发生改变,得重新根据公式求得所有帧之间的IMU预积分。
当偏置变换很小时,可以将αbk+1bk、βbk+1bk、γbk+1bk按其偏置的一阶近似来调整,否则就进行重新传递。【论文式[12]】(这部分只是抛出一个概念,后面会讲为什么这样写)
αbk+1bk≈α^bk+1bk+Jbaαδba+Jbwαδbw βbk+1bk≈β^bk+1bk+Jbaβδba+Jbwβδbw γbk+1bk≈γ^bk+1bk[121Jbwγδbw]
4.2 离散形式下的IMU运动模型
4.2.1 两帧之间 PVQ 增量的欧拉法离散形式
欧拉法给出第i时刻与第i+1时刻的预积分量估计值的关系: 【论文式[7]】
α^i+1bk=α^ibk+β^ibkδt+21R(γ^ibk)(a^i−bai)δt2 β^i+1bk=β^ibk+21R(γ^ibk)(a^i−bai)δt γ^i+1bk=γ^ibk⊗γ^i+1i=γ^ibk⊗[121(w^i−bwi)δt]
其中i是[tk,tk+1]之间的离散时间
4.2.2 两帧之间 PVQ 增量的中值法离散形式
代码中采用的基于中值法的 IMU 预积分公式,这在Estimator::processIMU()函数中的IntegrationBase::push_back()函数中得以实现,注意这里积分出来的是前后两帧之间的 IMU 增量信息。
α^i+1bk=α^ibk+β^ibkδt+21a^ˉiδt2 β^i+1bk=β^ibk+a^ˉiδt γ^i+1bk=γ^ibk⊗γ^i+1i=γ^ibk⊗[121w^ˉiδt] 其中 a^ˉt=21[γ^ibk(a^i−bi)+γ^i+1bk(a^i+1−bai)] w^ˉi=21(w^i+w^i+1)−bwi
初始状态下αbkbk、βbkbk为0, γbkbk为单位四元素,na、nw被视为0,i为在[k,k+1]中IMU测量值的某一时刻,δt为IMU测量值i和i+1之间的时间间隔。
5. 非线性方程的误差递推方程
5.1 一段时间内多个IMU数据积分形成的预积分量的协方差计算
一个 IMU 数据作为测量值的噪声方差我们能够标定。现在,一段时间内多个 IMU 数据积分形成的预积分量的方差呢?
要推导预积分量的协方差,需要知道IMU噪声和预积分量之间的线性递推关系。
假设已知了相邻时刻误差的线性传递方程:
ηik+1=Fkηik+Gknk其中ηik为状态量误差且ηik=[δθik,δvik,δpik],nk为测量噪声且nk=[nkg,nkg]。
可以看出误差的传递由倆部分组成:当前时刻的误差传递给下一时刻,当前时刻测量噪声传递给下一时刻。
5.2 非线性方程的误差递推方程推导
通常对于状态量之间的递推关系是非线性的方程如xk+1=f(xk,uk),其中状态量x、u为系统的输入量。
可以用俩种方法来推导状态误差传递的线性递推关系:
- 一种是基于一阶泰勒展开的误差递推方程
- 一种是基于误差随时间变化的递推方程 (论文是基于误差随时间变化来推导的)
5.2.1 基于一阶泰勒展开的误差递推方程
令状态量为x^=x+δx,其中x^表示计算得到的值,带有误差,真值为x,误差为δx。另外,输入量u的噪声为n。
应用一阶泰勒展开(EKF的协方差预测也是用了泰勒展开),非线性系统x^k+1=f(x^k,u^k)的状态误差的线性递推关系为:
δxk+1=Fδxk+Gnk 其中,F是状态量xk+1对状态量xk的雅克比矩阵,G是状态量xk+1对输入量uk的雅克比矩阵。
证明:对非线性状态方程进行一阶泰勒展开有:
x^k+1=f(x^k,u^k) xk+1+δxk+1=f(xk+δxk,uk+nk) xk+1+δxk+1=f(xk,uk)+Fδxk+Gnk(3)
5.2.2 基于误差随时间变化(误差导数)的递推方程
如果我们能够推导状态误差随时间变化的导数关系,如:
δ˙x=Aδx+Bn 则误差状态的传递方程为:
δxk+1=δxk+δ˙xkΔt →δxk+1=(I+AΔt)δxk+BΔtnk(4)
由公式(3)(4)得 F=I+AΔt, G=BΔt
5.2.3 基于泰勒展开和误差随时间变化的方法对比
基于一阶泰勒展开的误差递推方程不香吗,为什么会用误差随时间的变化来进行误差递推方程的构建呢?
在VIO系统中,我们已经知道了状态的导数和状态之间的转移矩阵,如我们已经知道速度和状态量之间的关系:
v˙w=Rbwab+gw 那我们就可以推导速度的误差导数和状态误差之间的关系,在每一项上都加上各自的误差就有:
v˙w+δ˙vw=Rbwexp([δθ]×)(ab+δab)+g+δg v˙w+δ˙vw=Rbw(I+[δθ]×)(ab+δab)+g+δg v˙w+δ˙vw=Rbwab+g+Rbwδab+Rbw[δθ]×(ab+δab)+δgδ˙vw=Rbwδab+Rbw[δθ]×(ab+δab)+δg δ˙vw=Rbwδab−Rbw[ab]×δθ+δg(5) 由此就能以此类推,轻易写出整个 A 和 B 其他方程了。
考虑到公式的编辑篇幅,为了对一些求导公式进行简化,对求导公式进行了简化 (下同):
∂δθ∂aw=δθ→0limδθRbwexp([δθ]×)ab−Rbwab→∂δθ∂aw=∂δθ∂Rbwexp([δθ]×)ab
6. PVQ增量误差、协方差及雅克比矩阵的递推方程
6.1 连续形式下的误差、协方差及雅克比矩阵的递推方程
在第4部分我们已经完成了IMU预积分测量值的推导,而要将IMU预积分运用到非线性优化中,需要建立线性高斯误差状态递推方程 (可参考第5部分),由线性高斯系统的协方差,推到方程协方差矩阵,并纠结对应的雅克比矩阵。
首先由于四元素γtbk被参数化过,我们将其误差定义为围绕其均值的扰动:【论文式(8)】
γtbk≈γ^tbk⊗[121δθtbk]
由预积分的连续形式建立的运动模型,由公式(2)可得一段时间内IMU构建的预积分量作为测量值,对俩时刻之间的状态量进行约束,可得到其在k和k+1时刻下的误差项为:【论文公式(24)】
⎣⎢⎢⎢⎢⎢⎡δαbk+1bkδβbk+1bkδθbk+1bkδbaδbw⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡Rwbk(pbk+1w−pbkw−vkwΔt+21gwΔt2)−αbk+1bkRwbk(vbk+1w−vbkw+gwΔt)−βbk+1bk2[qwbk⊗qbk+1w⊗γbkbk+1]xyzbabk+1−babkbwbk+1−bwbk⎦⎥⎥⎥⎥⎥⎤
我们给出在t时刻误差项的线性化递推方程为:【论文式(9)】:
⎣⎢⎢⎢⎢⎡δα˙tbkδβ˙tbkδθ˙tbkδb˙atδb˙wt⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡00000I00000−Rtbk[a^t−bat]×−[w^−bwt]×000−Rtbk00000−I00⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡δαtbkδβtbkδθtbkδbatδbwt⎦⎥⎥⎥⎥⎤+⎣⎢⎢⎢⎢⎡0−Rtbk00000−I00000I00000I⎦⎥⎥⎥⎥⎤⎣⎢⎢⎡nanwnbanbw⎦⎥⎥⎤=Ftδztbk+Gtnt(6) 可以简写为:
δz˙tbk=Ftδztbk+Gtnt 其中Ft是15x15,Gt是12x12,δztbk是15x1,nt是12x1。
根据导数的定义有:
δz˙tbk=δt→0limδtδzt+δtbk−δztbk 即:δzt+δtbk=δztbk+δz˙tbkδt=δztbk+(Ftδztbk+Gtnt)δt=(I+Ftδt)δztbk+Gtδtnt 这个形式与公式(4)不谋而合,一般连续形式下的微小时间用δt表示,离散形式下的时间区间用Δt表示,其实都是一个意思。进一步的,令F=(I+Ftδt),V=Gtδt, 有:
δzt+δtbk=Fδztbk+Vnt(7) 由此我们已经知道了该非线性系统的线性化递推方程。
至于下一时刻的协方差,首先对于一个高斯分布的变量线性变换后得到新的变量的协方差计算,我们有以下推导:
已知一个变量y=Ax,x∼(0,Σx),则有 Σy=AΣxAT,即:
Σy=E((Ax)(Ax)T)=E(AxxTAT)=AΣxAT
得到了递推方程(7)之后,根据高斯分布线性变换协方差的传递规律,我们可以预测下一时刻的协方差:【论文式(10)】
Pt+δtbk=(I+Ftδt)Ptbk(I+Ftδt)T+(Gδt)Q(Gtδt)T
其中初始协方差 Pbkbk=0, Q代表噪声项的对角协方差矩阵:
Q12×12=⎣⎢⎢⎡σa20000σw20000σba20000σbw2⎦⎥⎥⎤
同时也可以根据递推方程得到对应的雅克比矩阵迭代公式:【论文式[11]】
Jt+δt=(I+Ftδt)Jt 其中初始雅克比矩阵为:Jbk=I。
关于公式(6)的推导:
先把之前的预积分再写一遍:
αbk+1bk=∫∫t∈[tk,tk+1]Rtbk(a^t−bat−na)dt2 βbk+1bk=∫t∈[tk,tk+1]Rtbk(a^t−bat−na)dt γbk+1=∫t∈[tk,tk+1]21Ω(w^t−bwt−nw)γtbkdt 1)首先对于δα˙tbk、δb˙at、δb˙wt, 显然根据定义有:
δα˙tbk=α˙^tbk−α˙tbk=β^tbk−βtbk=δβtbk δb˙at=b˙at−0=n