VIO(notes) —— (3)VIO残差构建与IMU预积分

article/2025/9/17 1:43:00

VIO残差构建与IMU预积分

  • 一、VIO残差函数的构建
    • 1. 系统所需的状态变量
    • 2. 视觉重投影误差
      • 2.1 视觉重投影误差
      • 2.2 逆深度参数化
      • 2.3 VIO 中基于逆深度的重投影误差
    • 3. 预积分模型由来及意义
      • 3.1 为什么需要预积分?
      • 3.2 怎么预积分?
      • 3.3 预积分是什么?
    • 4. 预积分量方差的计算
    • 5. 预积分离散方法
    • 6. 预积分量的方差
    • 7. 状态误差线性递推公式的推导
      • 7. 1基于泰勒展开的误差传送(应用于 `EKF 的协方差`预测)
      • 7.2 基于误差随时间变化的递推方程
    • 8. 预积分的误差递推公式推导
  • 二、残差Jacobian的推导
    • 1. 视觉重投影残差的Jacobian
  • 参考文献

一、VIO残差函数的构建

1. 系统所需的状态变量

为了节约计算量采用滑动窗口形式的 Bundle Adjustment, 在 i i i 时刻, 滑动窗口内待优化的系统状态量定义如下:
X = [ x n , x n + 1 , … , x n + N , λ m , λ m + 1 , … , λ m + M ] x i = [ p w b i , q w b i , v i w , b a b i , b g b i ] ⊤ , i ∈ [ n , n + N ] \begin{aligned} &\mathcal{X}=\left[\mathbf{x}_{n}, \mathbf{x}_{n+1}, \ldots, \mathbf{x}_{n+N}, \lambda_{m}, \lambda_{m+1}, \ldots, \lambda_{m+M}\right] \\ &\mathbf{x}_{i}=\left[\mathbf{p}_{w b_{i}}, \mathbf{q}_{w b_{i}}, \mathbf{v}_{i}^{w}, \mathbf{b}_{a}^{b_{i}}, \mathbf{b}_{g}^{b_{i}}\right]^{\top}, i \in[n, n+N] \end{aligned} X=[xn,xn+1,,xn+N,λm,λm+1,,λm+M]xi=[pwbi,qwbi,viw,babi,bgbi],i[n,n+N]
其中:

  • x i \mathbf{x}_{i} xi 包含 i i i 时刻 IMU 机体的在惯性坐标系中的位置,速度,姿态, 以及 IMU 机体坐标系中的加速度和角速度的偏置量估计。
  • n , m n, m n,m 分别是机体状态量, 路标在滑动窗口里的起始时刻。
  • N N N 滑动窗口中关键帧数量。
  • M M M 是被滑动窗口内所有关键帧观测到的路标数量。

2. 视觉重投影误差

2.1 视觉重投影误差

定义: 一个特征点在归一化相机坐标系下的估计值与观测值的差。
r c = [ x z − u y z − v ] \mathbf{r}_{c}=\left[\begin{array}{l} \frac{x}{z}-u \\ \frac{y}{z}-v \end{array}\right] rc=[zxuzyv]
其中,待估计的状态量为特征点的三维空间坐标 ( x , y , z ) ⊤ (x, y, z)^{\top} (x,y,z), 观测值 ( u , v ) ⊤ (u, v)^{\top} (u,v) 为特征在相机归一化平面的坐标。

2.2 逆深度参数化

特征点在归一化相机坐标系与在相机坐标系下的坐标关系为:
[ x y z ] = 1 λ [ u v 1 ] \left[\begin{array}{l} x \\ y \\ z \end{array}\right]=\frac{1}{\lambda}\left[\begin{array}{l} u \\ v \\ 1 \end{array}\right] xyz=λ1uv1
其中 λ = 1 / z \lambda=1 / z λ=1/z 称为逆深度。

  • 为什么使用逆深度表示而不是深度值呢?
  1. 使用逆深度的缘故是相较于深度值其更符合高斯分布的特性;
  2. 深度值一般很大,会影响优化的节奏,使用倒数更能保证数值的稳定性

2.3 VIO 中基于逆深度的重投影误差

特征点逆深度在第 i i i 帧中初始化得
到,在第 j j j 帧又被观测到,预测其 在第 j j j 中的坐标为:
[ x c j y c j z c j 1 ] = T b c − 1 T w b j − 1 T w b i T b c [ 1 λ u c i 1 λ v c i 1 λ 1 ] \left[\begin{array}{c}x_{c_{j}} \\ y_{c_{j}} \\ z_{c_{j}} \\ 1\end{array}\right]=\mathbf{T}_{b c}^{-1} \mathbf{T}_{w b_{j}}^{-1} \mathbf{T}_{w b_{i}} \mathbf{T}_{b c}\left[\begin{array}{c}\frac{1}{\lambda} u_{c_{i}} \\ \frac{1}{\lambda} v_{c_{i}} \\ \frac{1}{\lambda} \\ 1\end{array}\right] xcjycjzcj1=Tbc1Twbj1TwbiTbcλ1uciλ1vciλ11
其中:
T = [ R t 0 T 1 ] \begin{aligned} T=\left[\begin{array}{ll} R & t \\ 0^{T} & 1 \end{array}\right] \end{aligned} T=[R0Tt1]

  1. 将i帧中观测到的数据变换到相机坐标系,将相机坐标系变换到body坐标系
  2. 将第i个body坐标系变换到世界坐标系;
  3. 将世界坐标系变换到第j个body坐标系;
  4. body坐标系变换到相机坐标系,得到第j帧的预测值。

这期间相对于纯视觉多了相机坐标系变换到body坐标系,然后再由body坐标系变换回相机坐标系的过程。

视觉重投影误差为:
r c = [ x c j z c j − u c j y c j z c j − v c j ] \begin{aligned} &\mathbf{r}_{c}=\left[\begin{array}{l} \frac{x_{c_{j}}}{z_{c_{j}}}-u_{c_{j}} \\ \frac{y_{c_{j}}}{z_{c_{j}}}-v_{c_{j}} \end{array}\right] \\ \end{aligned} rc=zcjxcjucjzcjycjvcj

3. 预积分模型由来及意义

3.1 为什么需要预积分?

IMU的测量值为 ω , b \omega, b ω,b ,则有
ω ~ b = ω b + b g + n g a ~ b = q b w ( a w + g w ) + b a + n a \begin{gathered} \tilde{\omega}^{b}=\omega^{b}+b^{g}+n^{g} \\ \tilde{a}^{b}=q_{b w}\left(a^{w}+g^{w}\right)+b^{a}+n^{a} \end{gathered} ω~b=ωb+bg+nga~b=qbw(aw+gw)+ba+na
PVQ对时间的导数可写成
p ˙ w b t = v t w v ˙ t w = a t w q ˙ w b t = q w b t ⊗ [ 0 1 2 ω b t ] \begin{gathered} \dot{p}_{w b_{t}}=v_{t}^{w} \\ \dot{v}_{t}^{w}=a_{t}^{w} \\ \dot{q}_{w b_{t}}=q_{w b_{t}} \otimes\left[\begin{array}{c} 0 \\ \frac{1}{2} \omega^{b_{t}} \end{array}\right] \end{gathered} p˙wbt=vtwv˙tw=atwq˙wbt=qwbt[021ωbt]
根据上面的导数关系,可以从第 i \mathrm{i} i 时刻的 P V Q \mathrm{PVQ} PVQ ,通过对 I M U \mathrm{IMU} IMU 的测量值进行积分,得到第时刻的 PVQ:
p w b j = p w b i + v i w Δ t + ∬ t ∈ [ i , j ] ( q w b t a b t − g w ) δ t 2 v j w = v i w + ∫ t ∈ [ i , j ] ( q w b t a b t − g w ) δ t q w b j = ∫ t ∈ [ i , j ] q w b t ⊗ [ 0 1 2 ω b t ] δ t \begin{aligned} &\mathbf{p}_{w b_{j}}=\mathbf{p}_{w b_{i}}+\mathbf{v}_{i}^{w} \Delta t+\iint_{t \in[i, j]}\left(\mathbf{q}_{w b_{t}} \mathbf{a}^{b_{t}}-\mathbf{g}^{w}\right) \delta t^{2} \\ &\mathbf{v}_{j}^{w}=\mathbf{v}_{i}^{w}+\int_{t \in[i, j]}\left(\mathbf{q}_{w b_{t}} \mathbf{a}^{b_{t}}-\mathbf{g}^{w}\right) \delta t \\ &\mathbf{q}_{w b_{j}}=\int_{t \in[i, j]} \mathbf{q}_{w b_{t}} \otimes\left[\begin{array}{c} 0 \\ \frac{1}{2} \boldsymbol{\omega}^{b_{t}} \end{array}\right] \delta t \end{aligned} pwbj=pwbi+viwΔt+t[i,j](qwbtabtgw)δt2vjw=viw+t[i,j](qwbtabtgw)δtqwbj=t[i,j]qwbt[021ωbt]δt

  • 问题: 为什么需要预积分?
    每次 q ω b t q_{\omega b_{t}} qωbt 更新后,都需要重新进行积分,运算量大;我们知道在紧耦合中待优化的状态变量是 p p p, q q q, v v v在优化的一次次迭代中这些变量是在不停变化的。由公式可以看到当 b i b_i bi时刻的状态发生改变的时候,我们需要一次次的积分来求取 b j b_j bj时刻的状态,这个是很耗费计算资源的,为了节省计算资源和时间,我们希望能够避免积分部分的重复计算。

3.2 怎么预积分?

一个很简单的公式转换,就可以将积分模型转为预积分模型:
q w b t = q w b i ⊗ q b i b t \mathbf{q}_{w b_{t}}=\mathbf{q}_{w b_{i}} \otimes \mathbf{q}_{b_{i} b_{t}} qwbt=qwbiqbibt
那么,PVQ 积分公式中的积分项则变成相对于第i 时刻的姿态,而不是相对于世界坐标系的姿态(注意,积分的变化):
p w b j = p w b i + v i w Δ t − 1 2 g w Δ t 2 + q w b i ∬ t ∈ [ i , j ] ( q b i b t a b t ) δ t 2 v j w = v i w − g w Δ t + q w b ı ∫ t ∈ [ i , j ] ( q b z b t a b t ) δ t q w b j = q w b i ∫ t ∈ [ i , j ] q b i b t ⊗ [ 0 1 2 ω b t ] δ t \begin{aligned} &\mathbf{p}_{w b_{j}}=\mathbf{p}_{w b_{i}}+\mathbf{v}_{i}^{w} \Delta t-\frac{1}{2} \mathbf{g}^{w} \Delta t^{2}+\mathbf{q}_{w b_{i}} \iint_{t \in[i, j]}\left(\mathbf{q}_{b_{i} b_{t}} \mathbf{a}^{b_{t}}\right) \delta t^{2} \\ &\mathbf{v}_{j}^{w}=\mathbf{v}_{i}^{w}-\mathbf{g}^{w} \Delta t+\mathbf{q}_{w b_{\imath}} \int_{t \in[i, j]}\left(\mathbf{q}_{b_{z} b_{t}} \mathbf{a}^{b_{t}}\right) \delta t \\ &\mathbf{q}_{w b_{j}}=\mathbf{q}_{w b_{i}} \int_{t \in[i, j]} \mathbf{q}_{b_{i} b_{t}} \otimes\left[\begin{array}{c} 0 \\ \frac{1}{2} \omega^{b_{t}} \end{array}\right] \delta t \end{aligned} pwbj=pwbi+viwΔt21gwΔt2+qwbit[i,j](qbibtabt)δt2vjw=viwgwΔt+qwbıt[i,j](qbzbtabt)δtqwbj=qwbit[i,j]qbibt[021ωbt]δt
这里的预积分量仅仅跟IMU 测量值有关,它将一段时间内的IMU 数据直接积分起来就得到了预积分量:
α b i b j = ∬ t ∈ [ i , j ] ( q b i b t a b t ) δ t 2 β b i b j = ∫ t ∈ [ i , j ] ( q b i b t a b t ) δ t q b i b j = ∫ t ∈ [ i , j ] q b i b t ⊗ [ 0 1 2 ω b t ] δ t \begin{aligned} \boldsymbol{\alpha}_{b_{i} b_{j}} &=\iint_{t \in[i, j]}\left(\mathbf{q}_{b_{i} b_{t}} \mathbf{a}^{b_{t}}\right) \delta t^{2} \\ \boldsymbol{\beta}_{b_{i} b_{j}} &=\int_{t \in[i, j]}\left(\mathbf{q}_{b_{i} b_{t}} \mathbf{a}^{b_{t}}\right) \delta t \\ \mathbf{q}_{b_{i} b_{j}} &=\int_{t \in[i, j]} \mathbf{q}_{b_{i} b_{t}} \otimes\left[\begin{array}{c} 0 \\ \frac{1}{2} \boldsymbol{\omega}^{b_{t}} \end{array}\right] \delta t \end{aligned} αbibjβbibjqbibj=t[i,j](qbibtabt)δt2=t[i,j](qbibtabt)δt=t[i,j]qbibt[021ωbt]δt
重新整理下PVQ 的积分公式,有:
[ p w b j v j w q w b j b j a b j g ] = [ p w b i + v i w Δ t − 1 2 g w Δ t 2 + q w b i α b i b j v i w − g w Δ t + q w b i β b i b j q w b q q b i b j b i a b i g ] \left[\begin{array}{c} \mathbf{p}_{w b_{j}} \\ \mathbf{v}_{j}^{w} \\ \mathbf{q}_{w b_{j}} \\ \mathbf{b}_{j}^{a} \\ \mathbf{b}_{j}^{g} \end{array}\right]=\left[\begin{array}{c} \mathbf{p}_{w b_{i}}+\mathbf{v}_{i}^{w} \Delta t-\frac{1}{2} \mathbf{g}^{w} \Delta t^{2}+\mathbf{q}_{w b_{i}} \alpha_{b_{i} b_{j}} \\ \mathbf{v}_{i}^{w}-\mathbf{g}^{w} \Delta t+\mathbf{q}_{w b_{i}} \boldsymbol{\beta}_{b_{i} b_{j}} \\ \mathbf{q}_{w b_{\mathbf{q}}} \mathbf{q}_{b_{i} b_{j}} \\ \mathbf{b}_{i}^{a} \\ \mathbf{b}_{i}^{g} \end{array}\right] pwbjvjwqwbjbjabjg=pwbi+viwΔt21gwΔt2+qwbiαbibjviwgwΔt+qwbiβbibjqwbqqbibjbiabig

3.3 预积分是什么?

预积分其实是加速度和角速度引起的 b j b_j bj时刻的状态量关于 b i b_i bi时刻的增量
所谓的预积分,就是预先积分好,之后需要的时候拿来用。我们知道IMU测量的是线性加速度和角速度。状态量中速度是关于加速度的积分,位置是关于加速度的二次积分,姿态是关于角速度的积分。因此,在当前坐标系下(因为IMU的测量值就是在当前坐标系下的)把IMU测量值积分好,当需要的是,只需要根据情况进行坐标系转换和加减了。
可以看到预积分与待优化的状态量(忽略零偏)无关,它是一个固定的值。

4. 预积分量方差的计算

预积分误差定义:一段时间内IMU 构建的预积分量作为测量值,对两时刻之间的状态量进行约束
[ r p r q r v r b a r b g ] 15 × 1 = [ q b i w ( p w b j − p w b i − v i w Δ t + 1 2 g w Δ t 2 ) − α b i b j 2 [ q b j b i ⊗ ( q b i w ⊗ q w b j ) ] x y z q b i w ( v j w − v i w + g w Δ t ) − β b i b j b j a − b i a b j g − b i g ] \left[\begin{array}{c} \mathbf{r}_{p} \\ \mathbf{r}_{q} \\ \mathbf{r}_{v} \\ \mathbf{r}_{b a} \\ \mathbf{r}_{b g} \end{array}\right]_{15 \times 1}=\left[\begin{array}{c} \mathbf{q}_{b_{i} w}\left(\mathbf{p}_{w b_{j}}-\mathbf{p}_{w b_{i}}-\mathbf{v}_{i}^{w} \Delta t+\frac{1}{2} \mathbf{g}^{w} \Delta t^{2}\right)-\alpha_{b_{i} b_{j}} \\ 2\left[\mathbf{q}_{b_{j} b_{i}} \otimes\left(\mathbf{q}_{b_{i} w} \otimes \mathbf{q}_{w b_{j}}\right)\right]_{x y z} \\ \mathbf{q}_{b_{i} w}\left(\mathbf{v}_{j}^{w}-\mathbf{v}_{i}^{w}+\mathbf{g}^{w} \Delta t\right)-\beta_{b_{i} b_{j}} \\ \mathbf{b}_{j}^{a}-\mathbf{b}_{i}^{a} \\ \mathbf{b}_{j}^{g}-\mathbf{b}_{i}^{g} \end{array}\right] rprqrvrbarbg15×1=qbiw(pwbjpwbiviwΔt+21gwΔt2)αbibj2[qbjbi(qbiwqwbj)]xyzqbiw(vjwviw+gwΔt)βbibjbjabiabjgbig
上面误差中位移,速度,偏置都是直接相减得到。第二顶是关于四元数的旋转误差,其中 [.] xyz 表示只取四元数的虚部 ( x , y , z ) (x, y, z) (x,y,z) 组成的三维向量。

5. 预积分离散方法

关于预积分的计算,前面提到过有欧拉法与中值法。这里使用mid-point 方法,即两个相邻时刻k 到 k + 1 k+1 k+1 的位姿是用两个时刻的测量值 a , w a, w a,w 的平均值来计算:
ω = 1 2 ( ( ω b k − b k g ) + ( ω b k + 1 − b k g ) ) q b i b k + 1 = q b i b k ⊗ [ 1 1 2 ω δ t ] a = 1 2 ( q b i b k ( a b k − b k a ) + q b i b k + 1 ( a b k + 1 − b k a ) ) α b i b k + 1 = α b i b k + β b i b k δ t + 1 2 a δ t 2 β b i b k + 1 = β b i b k + a δ t b k + 1 a = b k a + n b k a δ t b k + 1 g = b k g + n b k g δ t \begin{aligned} \boldsymbol{\omega} &=\frac{1}{2}\left(\left(\boldsymbol{\omega}^{b_{k}}-\mathbf{b}_{k}^{g}\right)+\left(\boldsymbol{\omega}^{b_{k+1}}-\mathbf{b}_{k}^{g}\right)\right) \\ \mathbf{q}_{b_{i} b_{k+1}} &=\mathbf{q}_{b_{i} b_{k}} \otimes\left[\begin{array}{c} 1 \\ \frac{1}{2} \boldsymbol{\omega} \delta t \end{array}\right] \\ \mathbf{a} &=\frac{1}{2}\left(\mathbf{q}_{b_{i} b_{k}}\left(\mathbf{a}^{b_{k}}-\mathbf{b}_{k}^{a}\right)+\mathbf{q}_{b_{i} b_{k+1}}\left(\mathbf{a}^{b_{k+1}}-\mathbf{b}_{k}^{a}\right)\right) \\ \alpha_{b_{i} b_{k+1}} &=\alpha_{b_{i} b_{k}}+\beta_{b_{i} b_{k}} \delta t+\frac{1}{2} \mathbf{a} \delta t^{2} \\ \beta_{b_{i} b_{k+1}} &=\beta_{b_{i} b_{k}}+\mathbf{a} \delta t \\ \mathbf{b}_{k+1}^{a} &=\mathbf{b}_{k}^{a}+\mathbf{n}_{\mathbf{b}_{k}^{a} \delta t} \\ \mathbf{b}_{k+1}^{g} &=\mathbf{b}_{k}^{g}+\mathbf{n}_{\mathbf{b}_{k}^{g} \delta t} \end{aligned} ωqbibk+1aαbibk+1βbibk+1bk+1abk+1g=21((ωbkbkg)+(ωbk+1bkg))=qbibk[121ωδt]=21(qbibk(abkbka)+qbibk+1(abk+1bka))=αbibk+βbibkδt+21aδt2=βbibk+aδt=bka+nbkaδt=bkg+nbkgδt

6. 预积分量的方差

疑问?:
一个 IMU 数据作为测量值的噪声方差我们能够标定。现在,一段时间内多个 IMU 数据积分形成的预积分量的方差呢?
Covariance Propagation(协方差传播)
举例说明:
已知一个变量 y = A x , x ∈ N ( 0 , Σ x ) y=A x, x \in N\left(0, \Sigma_{x}\right) y=Ax,xN(0,Σx), 则有 Σ y = A Σ x A ⊤ \Sigma_{y}=A \Sigma_{x} A^{\top} Σy=AΣxA
Σ y = E ( ( A x ) ( A x ) ⊤ ) = E ( A x x ⊤ A ⊤ ) = A Σ x A ⊤ \begin{aligned} \Sigma_{y}=& E\left((A x)(A x)^{\top}\right) \\ =& E\left(A x x^{\top} A^{\top}\right) \\ &=A \Sigma x A^{\top} \end{aligned} Σy==E((Ax)(Ax))E(AxxA)=AΣxA
所以,要推导预积分量的协方差,我们需要知道 imu 噪声和预积分量之间的线性递推关系
假设已知了相邻时刻误差的线性传递方程:
η i k = F k − 1 η i k − 1 + G k − 1 n k − 1 \eta_{i k}=F_{k-1} \eta_{i k-1}+G_{k-1} n_{k-1} ηik=Fk1ηik1+Gk1nk1
比如: 状态量误差为 η i k = [ δ θ i k , δ v i k , δ p i k ] \eta_{i k}=\left[\delta \theta_{i k}, \delta v_{i k}, \delta p_{i k}\right] ηik=[δθik,δvik,δpik], 测量唱声为 n k = [ n k g , n k a ] n_{k}=\left[n_{k}^{g}, n_{k}^{a}\right] nk=[nkg,nka]
误差的传递由两部分组成当前时刻的误差 传递给下一时刻,当前时刻测量噪声传递给下一时刻。
一个有趣的例子:
综艺节目中常有传递信息的节目,前一个人根据上一个人的信息 + 自己的理解(测量)传递给下一个人,导致这个信息越传越错。 协方差矩阵可以通过递推计算得到:
Σ i k = F k − 1 Σ i k − 1 F k − 1 ⊤ + G k − 1 Σ n G k − 1 T \Sigma_{i k}=F_{k-1} \Sigma_{i k-1} F_{k-1}^{\top}+G_{k-1} \Sigma_{n} G_{k-1}^{T} Σik=Fk1Σik1Fk1+Gk1ΣnGk1T
其中, Σ n \Sigma_{n} Σn 是测量噪声的协方差矩阵, 方差从 i 时刻开始进行递推, Σ i i = 0 \Sigma_{i i}=0 Σii=0

7. 状态误差线性递推公式的推导

通常对于状态量之间的递推关系是非线性的方程如 x k = f ( x k − 1 , u k − 1 ) x_{k}=f\left(x_{k-1}, u_{k-1}\right) xk=f(xk1,uk1), 其中状态量为 x , u x, u x,u 为系统的输入量。

  • 我们可以用两种方法来推导状态误差传递的线性递推关系:
  1. 一种是基于一阶泰勒展开的误差递推方程。
  2. 一种是基于误差随时间变化的递推方程。

7. 1基于泰勒展开的误差传送(应用于 EKF 的协方差预测)

令状态量为 x = x ^ + δ x x=\hat{x}+\delta x x=x^+δx, 其中, 真值为 x ^ \hat{x} x^, 误差为 δ x \delta x δx 。另外, 输入量 u u u 的噪声为 n n n
非线性系统 x k = f ( x k − 1 , u k − 1 ) x_{k}=f\left(x_{k-1}, u_{k-1}\right) xk=f(xk1,uk1) 的状态误差的线性递推关系如下:
δ x k = F δ x k − 1 + G n k − 1 \delta x_{k}=F \delta x_{k-1}+G n_{k-1} δxk=Fδxk1+Gnk1
其中, F \mathrm{F} F 是状态量 x k x_{k} xk 对状态量 x k − 1 x_{k-1} xk1 的雅克比矩阵, G \mathrm{G} G 是状态量 x k x_{k} xk 对输入量 u k − 1 u_{k-1} uk1 的雅克比矩阵。
证明:对非线性状态方程进行一阶泰勒展开有:
x k = f ( x k − 1 , u k − 1 ) x ^ k + δ x k = f ( x ^ k − 1 + δ x k − 1 , u ^ k − 1 + n k − 1 ) x ^ k + δ x k = f ( x ^ k − 1 , u ^ k − 1 ) + F δ x k − 1 + G n k − 1 \begin{gathered} x_{k}=f\left(x_{k-1}, u_{k-1}\right) \\ \hat{x}_{k}+\delta x_{k}=f\left(\hat{x}_{k-1}+\delta x_{k-1}, \hat{u}_{k-1}+n_{k-1}\right) \\ \hat{x}_{k}+\delta x_{k}=f\left(\hat{x}_{k-1}, \hat{u}_{k-1}\right)+F \delta x_{k-1}+G n_{k-1} \end{gathered} xk=f(xk1,uk1)x^k+δxk=f(x^k1+δxk1,u^k1+nk1)x^k+δxk=f(x^k1,u^k1)+Fδxk1+Gnk1

7.2 基于误差随时间变化的递推方程

如果我们能够推导状态误差随时间变化的导数关系, 比如:
δ x ′ = A δ x + B n \delta x^{\prime}=A \delta x+B n δx=Aδx+Bn
则误差状态的传递方程为:
δ x k = δ x k − 1 + δ x k − 1 ′ Δ t → δ x k = ( I + A Δ t ) δ x k − 1 + B Δ t n k − 1 \begin{gathered} \delta x_{k}=\delta x_{k-1}+\delta x_{k-1}^{\prime} \Delta t \\ \rightarrow \delta x_{k}=(I+A \Delta t) \delta x_{k-1}+B \Delta t n_{k-1} \end{gathered} δxk=δxk1+δxk1Δtδxk=(I+AΔt)δxk1+BΔtnk1
这两种推导方式的可以看出有:
F = I + A Δ t , G = B Δ t F=I+A \Delta t, G=B \Delta t F=I+AΔt,G=BΔt

  • 第一种方法不是很好么,为什么会随着去弄误差随时间的变化呢?

这是因为 VIO 系统中已经知道了状态的导数和状态之间的转移矩阵。如:我们已经知道速度和状态量之间的关系:
v ˙ = R a b + g \dot{v}=R a^{b}+g v˙=Rab+g
那我们就可以推导速度的误差和状态误差之间的关系,再每一项上都加上各自的误差 就有:
v ˙ + δ v ˙ = R ( I + [ δ θ ] × ) ( a b + δ a b ) + ( g + δ g ) δ v ˙ = R δ a b + R [ δ θ ] × ( a b + δ a b ) + δ g δ v ˙ = R δ a b − R [ a b ] × δ θ + δ g \begin{gathered} \dot{v}+\delta \dot{v}=R\left(I+[\delta \theta]_{\times}\right)\left(a^{b}+\delta a^{b}\right)+(g+\delta g) \\ \delta \dot{v}=R \delta a^{b}+R[\delta \theta]_{\times}\left(a^{b}+\delta a^{b}\right)+\delta g \\ \delta \dot{v}=R \delta a^{b}-R\left[a^{b}\right]_{\times} \delta \theta+\delta g \end{gathered} v˙+δv˙=R(I+[δθ]×)(ab+δab)+(g+δg)δv˙=Rδab+R[δθ]×(ab+δab)+δgδv˙=RδabR[ab]×δθ+δg
由此就能依次类推,轻易写出整个 A A A B B B 其他方程了。

8. 预积分的误差递推公式推导

首先回顾预积分的误差递推公式, 将测量噪声也考虑进模型:
ω = 1 2 ( ( ω b k + n k g − b k g ) + ( ω b k + 1 + n k + 1 g − b k g ) ) q b i b k + 1 = q b i b k ⊗ [ 1 1 2 ω δ t ] a = 1 2 ( q b i b k ( a b k + n k a − b k a ) + q b i b k + 1 ( a b k + 1 + n k + 1 a − b k a ) ) α b i b k + 1 = α b i b k + β b i b k δ t + 1 2 a δ t 2 β b i b k + 1 = β b i b k + a δ t b k + 1 a = b k a + n b k a δ t b k + 1 g = b k g + n b k g δ t \begin{aligned} \boldsymbol{\omega} &=\frac{1}{2}\left(\left(\boldsymbol{\omega}^{b_{k}}+\mathbf{n}_{k}^{g}-\mathbf{b}_{k}^{g}\right)+\left(\boldsymbol{\omega}^{b_{k+1}}+\mathbf{n}_{k+1}^{g}-\mathbf{b}_{k}^{g}\right)\right) \\ \mathbf{q}_{b_{i} b_{k+1}} &=\mathbf{q}_{b_{i} b_{k}} \otimes\left[\begin{array}{c} 1 \\ \frac{1}{2} \boldsymbol{\omega} \delta t \end{array}\right] \\ \mathbf{a} &=\frac{1}{2}\left(\mathbf{q}_{b_{i} b_{k}}\left(\mathbf{a}^{b_{k}}+\mathbf{n}_{k}^{a}-\mathbf{b}_{k}^{a}\right)+\mathbf{q}_{b_{i} b_{k+1}}\left(\mathbf{a}^{b_{k+1}}+\mathbf{n}_{k+1}^{a}-\mathbf{b}_{k}^{a}\right)\right) \\ \boldsymbol{\alpha}_{b_{i} b_{k+1}} &=\boldsymbol{\alpha}_{b_{i} b_{k}}+\boldsymbol{\beta}_{b_{i} b_{k}} \delta t+\frac{1}{2} \mathbf{a} \delta t^{2} \\ \boldsymbol{\beta}_{b_{i} b_{k+1}} &=\boldsymbol{\beta}_{b_{i} b_{k}}+\mathbf{a} \delta t \\ \mathbf{b}_{k+1}^{a} &=\mathbf{b}_{k}^{a}+\mathbf{n}_{\mathbf{b}_{k}^{a}} \delta t \\ \mathbf{b}_{k+1}^{g} &=\mathbf{b}_{k}^{g}+\mathbf{n}_{\mathbf{b}_{k}^{g}} \delta t \end{aligned} ωqbibk+1aαbibk+1βbibk+1bk+1abk+1g=21((ωbk+nkgbkg)+(ωbk+1+nk+1gbkg))=qbibk[121ωδt]=21(qbibk(abk+nkabka)+qbibk+1(abk+1+nk+1abka))=αbibk+βbibkδt+21aδt2=βbibk+aδt=bka+nbkaδt=bkg+nbkgδt
确定误差传递的状态量,噪声量,然后开始构建传递方程。
用前面一阶泰勒展开的推导方式 δ x k = F δ x k − 1 + G n k − 1 \delta x_{k}=F \delta x_{k-1}+G n_{k-1} δxk=Fδxk1+Gnk1, 我们㠻望能推导出如下的形式:
[ δ α b k + 1 b k + 1 ′ δ θ b k + 1 b k + 1 ′ δ β b k + 1 b k + 1 ′ δ b k + 1 a δ b k + 1 g ] = F [ δ α b k b k ′ δ θ b k b k ′ δ β b k b k ′ δ b k a δ b k g ] + G [ n k a n k g n k + 1 a n k + 1 g n b k a n b k g ] \left[\begin{array}{c} \delta \boldsymbol{\alpha}_{b_{k+1} b_{k+1}^{\prime}} \\ \delta \boldsymbol{\theta}_{b_{k+1} b_{k+1}^{\prime}} \\ \delta \boldsymbol{\beta}_{b_{k+1} b_{k+1}^{\prime}} \\ \delta \mathbf{b}_{k+1}^{a} \\ \delta \mathbf{b}_{k+1}^{g} \end{array}\right]=\mathbf{F}\left[\begin{array}{c} \delta \boldsymbol{\alpha}_{b_{k} b_{k}^{\prime}} \\ \delta \boldsymbol{\theta}_{b_{k} b_{k}^{\prime}} \\ \delta \boldsymbol{\beta}_{b_{k} b_{k}^{\prime}} \\ \delta \mathbf{b}_{k}^{a} \\ \delta \mathbf{b}_{k}^{g} \end{array}\right]+\mathbf{G}\left[\begin{array}{c} \mathbf{n}_{k}^{a} \\ \mathbf{n}_{k}^{g} \\ \mathbf{n}_{k+1}^{a} \\ \mathbf{n}_{k+1}^{g} \\ \mathbf{n}_{\mathbf{b}_{k}^{a}} \\ \mathbf{n}_{\mathbf{b}_{k}^{g}} \end{array}\right] δαbk+1bk+1δθbk+1bk+1δβbk+1bk+1δbk+1aδbk+1g=Fδαbkbkδθbkbkδβbkbkδbkaδbkg+Gnkankgnk+1ank+1gnbkanbkg
F , G \mathrm{F}, \mathrm{G} F,G 为两个时刻间的协方差传递矩阵。
这里我们直接给出 F , G F, G F,G 的最终形式,后面会对部分项进行详细推导:
F = [ I f 12 I δ t − 1 4 ( q b i b k + q b i b k + 1 ) δ t 2 f 15 0 I − [ ω ] × 0 0 − I δ t 0 f 32 I − 1 2 ( q b i b k + q b i b k + 1 ) δ t f 35 0 0 0 I 0 0 0 0 0 I ] \mathbf{F}=\left[\begin{array}{ccccc} \mathbf{I} & \mathbf{f}_{12} & \mathbf{I} \delta t & -\frac{1}{4}\left(\mathbf{q}_{b_{i} b_{k}}+\mathbf{q}_{b_{i} b_{k+1}}\right) \delta t^{2} & \mathbf{f}_{15} \\ \mathbf{0} & \mathbf{I}-[\boldsymbol{\omega}]_{\times} & \mathbf{0} & \mathbf{0} & -\mathbf{I} \delta t \\ \mathbf{0} & \mathbf{f}_{32} & \mathbf{I} & -\frac{1}{2}\left(\mathbf{q}_{b_{i} b_{k}}+\mathbf{q}_{b_{i} b_{k+1}}\right) \delta t & \mathbf{f}_{35} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I} \end{array}\right] F=I0000f12I[ω]×f3200Iδt0I0041(qbibk+qbibk+1)δt2021(qbibk+qbibk+1)δtI0f15Iδtf350I
G = [ 1 4 q b i b k δ t 2 g 12 1 4 q b i b k + 1 δ t 2 g 14 0 0 0 1 2 I δ t 0 1 2 I δ t 0 0 1 2 q b i b k δ t g 32 1 2 q b i b k + 1 δ t g 34 0 0 0 0 0 0 I δ t 0 0 0 0 0 0 I δ t ] \mathbf{G}=\left[\begin{array}{cccccc} \frac{1}{4} \mathbf{q}_{b_{i} b_{k}} \delta t^{2} & \mathbf{g}_{12} & \frac{1}{4} \mathbf{q}_{b_{i} b_{k+1}} \delta t^{2} & \mathbf{g}_{14} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \frac{1}{2} \mathbf{I} \delta t & \mathbf{0} & \frac{1}{2} \mathbf{I} \delta t & \mathbf{0} & \mathbf{0} \\ \frac{1}{2} \mathbf{q}_{b_{i} b_{k}} \delta t & \mathbf{g}_{32} & \frac{1}{2} \mathbf{q}_{b_{i} b_{k+1}} \delta t & \mathbf{g}_{34} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I} \delta t & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I} \delta t \end{array}\right] G=41qbibkδt2021qbibkδt00g1221Iδtg320041qbibk+1δt2021qbibk+1δt00g1421Iδtg3400000Iδt00000Iδt
在这里插入图片描述

二、残差Jacobian的推导

1. 视觉重投影残差的Jacobian

视觉残差为:
对于第 i \mathrm{i} i 帧中的特征点, 它投影到第 j \mathrm{j} j 帧相机坐标系下的值为:
[ x c j y c j z c j 1 ] = T b c − 1 T w b j − 1 T w b i T b c [ 1 λ u c i 1 λ v c i 1 λ 1 ] \left[\begin{array}{c} x_{c_{j}} \\ y_{c_{j}} \\ z_{c_{j}} \\ 1 \end{array}\right]=\mathbf{T}_{b c}^{-1} \mathbf{T}_{w b_{j}}^{-1} \mathbf{T}_{w b_{i}} \mathbf{T}_{b c}\left[\begin{array}{c} \frac{1}{\lambda} u_{c_{i}} \\ \frac{1}{\lambda} v_{c_{i}} \\ \frac{1}{\lambda} \\ 1 \end{array}\right] xcjycjzcj1=Tbc1Twbj1TwbiTbcλ1uciλ1vciλ11
拆成三维坐标形式为:
f c j = [ x c j y c j z c j ] = R b c ⊤ R w b j ⊤ R w b i R b c 1 λ [ u c i v c i 1 ] + R b c ⊤ ( R w b j ⊤ ( ( R w b i p b c + p w b i ) − p w b j ) − p b c ) \begin{aligned} \mathbf{f}_{c_{j}}=\left[\begin{array}{c} x_{c_{j}} \\ y_{c_{j}} \\ z_{c_{j}} \end{array}\right] &=\mathbf{R}_{b c}^{\top} \mathbf{R}_{w b_{j}}^{\top} \mathbf{R}_{w b_{i}} \mathbf{R}_{b c} \frac{1}{\lambda}\left[\begin{array}{c} u_{c_{i}} \\ v_{c_{i}} \\ 1 \end{array}\right] \\ &+\mathbf{R}_{b c}^{\top}\left(\mathbf{R}_{w b_{j}}^{\top}\left(\left(\mathbf{R}_{w b_{i}} \mathbf{p}_{b c}+\mathbf{p}_{w b_{i}}\right)-\mathbf{p}_{w b_{j}}\right)-\mathbf{p}_{b c}\right) \end{aligned} fcj=xcjycjzcj=RbcRwbjRwbiRbcλ1ucivci1+Rbc(Rwbj((Rwbipbc+pwbi)pwbj)pbc)
再推导各类 Jacobian 之前, 为了简化公式, 先定义如下变量:
f b i = R b c f c i + p b c f w = R w b i f b i + p w b i f b j = R w b j ⊤ ( f w − p w b j ) \begin{aligned} \mathbf{f}_{b_{i}} &=\mathbf{R}_{b c} \mathbf{f}_{c_{i}}+\mathbf{p}_{b c} \\ \mathbf{f}_{w} &=\mathbf{R}_{w b_{i}} \mathbf{f}_{b_{i}}+\mathbf{p}_{w b_{i}} \\ \mathbf{f}_{b_{j}} &=\mathbf{R}_{w b_{j}}^{\top}\left(\mathbf{f}_{w}-\mathbf{p}_{w b_{j}}\right) \end{aligned} fbifwfbj=Rbcfci+pbc=Rwbifbi+pwbi=Rwbj(fwpwbj)
Jacobian 为视觉误差对两个时刻的状态量,外参,以及逆深度求导:
J = [ ∂ r c ∂ [ δ p b i b i ′ δ θ b i b i ′ ] ∂ r c ∂ [ p b j b j ′ δ θ b j b j ′ ] ∂ r c ∂ [ δ p c c ′ δ θ c c ′ ] ∂ r c ∂ δ λ ] \mathbf{J}=\left[\begin{array}{lll}\frac{\partial \mathbf{r}_{c}}{\partial\left[\begin{array}{l}\delta \mathbf{p}_{b_{i} b_{i}^{\prime}} \\ \delta \boldsymbol{\theta}_{b_{i} b_{i}^{\prime}}\end{array}\right]} \quad \frac{\partial \mathbf{r}_{c}}{\partial\left[\begin{array}{l}\mathbf{p}_{b_{j} b_{j}^{\prime}} \\ \delta \boldsymbol{\theta}_{b_{j} b_{j}^{\prime}}\end{array}\right]} \quad \frac{\partial \mathbf{r}_{c}}{\partial\left[\begin{array}{l}\delta \mathbf{p}_{c c^{\prime}} \\ \delta \boldsymbol{\theta}_{c c^{\prime}}\end{array}\right]} \quad \frac{\partial \mathbf{r}_{c}}{\partial \delta \lambda}\end{array}\right] J=[[δpbibiδθbibi]rc[pbjbjδθbjbj]rc[δpccδθcc]rcδλrc]

参考文献

VIO 中残差雅可比的推导


http://chatgpt.dhexx.cn/article/LeWbHf9b.shtml

相关文章

深蓝学院VIO课程学习笔记 VIO概述

VIO概述 1. VIO整体概述 松耦合:各部分自己算自己的,最后单独把数据来算 紧耦合:同时考虑这两个问题(效果更好) IMUGPS精度可以达到cm级,但是受环境影响比较大 融合方案 采用卡尔曼滤波,当一边…

海思3559 sample解析:vio

前言 拿到开发板,编完了平台sample,自然按捺不住要去简单学习测试了。打开最直观相对也比较简单的vio例程做个到手分析和流程梳理吧 测试 一开始自然是最磕磕绊绊的,连上HDMI线,串口登录后运行,屏幕乌漆嘛黑&#xff…

从零手写VIO(7)

从零手写VIO(7) 文章目录 从零手写VIO(7)前言一、VINS-Course代码解析二、作业(7)1.simulation-test.cpp修改2.Sysytem.cpp修改3.config文件夹下euroc_config.yaml参数修改4.param.h修改4.1无噪声4.2小噪声4.3大噪声 总结 前言 一、VINS-Course代码解析…

运行msckf_vio

1、编译 cd ~/msckf catkin_make --pkg msckf_vio --cmake-args -DCMAKE_BUILD_TYPERelease2、运行(运行euroc数据集) 首先roscore开启ros节点 cd ~/msckf source ~/msckf/devel/setup.bash roslaunch msckf_vio msckf_vio_euroc.launchcd ~/msckf sou…

深蓝学院-手写VIO作业-第二章

文章目录 一、基础作业,必做环境配置说明a. ROS环境搭建b. Matlab安装 1、设置IMU 仿真代码中的不同的参数,生成Allen 方差标定曲线。a. 对于非ROS:生成运动imu数据b.对于ROS: 专门生成静止 imu 数据,用于 allan 方差标定 2、将IM…

【代码阅读】PL-VIO

〇、写在前面 PL-VIO采用的通信是ROS,所以并不能像ORBSLAM那样按照执行顺序来理顺,因为ORBSLAM是有一个真正意义上的主函数的,经过CMakeList的编辑产生的可执行文件会有一个开始,但是PL-VIO用的是ROS,其内部通信是节点…

VIO:飞行机器人单目VIO算法测评

转:https://blog.csdn.net/u012348774/article/details/81414264 泡泡图灵智库,带你精读机器人顶级会议文章 标题:A Benchmark Comparison of Monocular Visual-Inertial Odometry Algorithms for Flying Robots 作者:Jeffrey De…

VIO学习笔记一

1. IMU(Inertial Measurement Unit,惯性测量单元) 典型6轴IMU以较高频率(≥100Hz)返回被测量物体的角速度与加速度。受自身温度、零偏、振动等因素干扰,积分得到的平移和旋转容易漂移。IMU本身由一个陀螺仪…

VIO仿真

使用turtelbot3仿真,发现gazebo的imu没有重力加速度。放弃。还是使用公开数据集。 使用这个仿真​​​​​​vio_data_simulation/src at ros_version HeYijia/vio_data_simulation GitHub 看一下如何用这个仿真跑起来vio。 -- 将特征点反投回图像,…

3.4.1 VIO虚拟以太网原理

最后更新2021/08/12 VIO Server在此完全实现了一个标准的以太网交换机的功能,现在业界都有了高大上的名称:SDN(Software Defined Network),如果没有足够的网络知识(幸好只是网络链路层)&#x…

Vivado调用VIO核

文章目录 前言一、IP核的介绍二、VIO核1.作用2.调用方法 总结 前言 提示:本篇文章所使用的软件为Vivado2018.3: 以四选一数据选择器为例,使用verilog hdl语言以及Vivado自带的VIO,IP来实现功能 提示:以下是本篇文章正文内容&…

海思3519 VIO Sample例程讲解

海思VIO Sample例程讲解 海思SDK解压出来后,Sample包含各个功能模块的历程,本篇讲解VIO Sample历程。 进入VIO模块可以看到,VIO的main函数文件,先从main函数执行程序。 进入文件后首先看下VIO实现的功能,可以看到VIO…

PL-VIO论文阅读

PL-VIO: Tightly-Coupled Monocular Visual–Inertial Odometry Using Point and Line Features Yijia He 1,2,* , Ji Zhao 3, Yue Guo 1,2, Wenhao He 1 and Kui Yuan 1 2018 摘要 To address the problem of estimating camera trajectory and to build a structural 3D m…

DM-VIO简析

今天主要是针对DMVIO/DM-VIO的简析,中文网上有的东西都太少了,只能靠看完论文和组员们一起改代码。Lukas组这个东西在中文网被称为有史以来最好的VIO,但是实际过程中我们还是发现了许多不完美的地方。。。(比如ZUPT更新改造中该有的问题仍然在…

VIOSLAM 综述

文章目录 1.VIO 松耦合/紧耦合。2. 相机和IMU的缺点及互补性3. VIO融合算法流程及其模块分解:4. VIO 算法核心:5. 实验结果与总结:6. 参考文献: 1.VIO 松耦合/紧耦合。 Visual-Inertial Odometry(VIO)即视觉惯性里程计,有时也叫视觉惯性系统…

VIO系统介绍

VIO(visual-inertial odometry)即视觉惯性里程计,有时也叫视觉惯性系统(VINS,visual-inertial system),是融合相机和IMU数据实现SLAM的算法,根据融合框架的区别又分为紧耦合和松耦合…

vivado VIO (virtual input output)虚拟IO的使用

转自:https://blog.csdn.net/wordwarwordwar/article/details/77150930 一般情况下ILA和VIO都是用在chipscope上使用,VIO可以作为在chipscope时模拟IO。 譬如: 在使用chipscope时需要使用按键出发,但是没有设计按键或者板子不再身…

【Vivado那些事儿】-VIO原理及应用

虚拟输入输出(Virtual Input Output,VIO)核是一个可定制的IP核,它可用于实时监视和驱动内部FPGA的信号,如图所示。 可以定制VIO的输入和输出端口的数量与宽度,用于和FPGA设计进行连接。由于VIO核与被监视和驱动的设计同步&#xf…

python logger.exception_Python logging设置和logger解析

一、logging模块讲解 1.函数:logging.basicConfig() 参数讲解: (1)level代表高于或者等于这个值时,那么我们才会记录这条日志 (2)filename代表日志会写在这个文件之中,如果没有这个字段则会显示在控制台上 (3)format代表我们的日志显示的格式自定义,如果字段为空,那么默认…

Logger 基本用法

Logger 基本用法 简介 Simple, pretty and powerful logger for android 为Android提供的,简单、强大而且格式美观的工具 本质就是封装系统提供的Log类,加上一些分割线易于查找不同的Log;logcat中显示的信息可配置。最初的样子如下图 包含…