IT数码 购物 网址 头条 软件 日历 阅读 图书馆
TxT小说阅读器
↓语音阅读,小说下载,古典文学↓
图片批量下载器
↓批量下载图片,美女图库↓
图片自动播放器
↓图片自动播放器↓
一键清除垃圾
↓轻轻一点,清除系统垃圾↓
开发: C++知识库 Java知识库 JavaScript Python PHP知识库 人工智能 区块链 大数据 移动开发 嵌入式 开发工具 数据结构与算法 开发测试 游戏开发 网络协议 系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑 笔记本 显卡 显示器 固态硬盘 硬盘 耳机 手机 iphone vivo oppo 小米 华为 单反 装机 图拉丁
 
   -> 数据结构与算法 -> VIO(notes) —— (3)VIO残差构建与IMU预积分 -> 正文阅读

[数据结构与算法]VIO(notes) —— (3)VIO残差构建与IMU预积分

一、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?=[zx??uzy??v?]
其中,待估计的状态量为特征点的三维空间坐标 ( 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????=λ1????uv1????
其中 λ = 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] ?????xcj??ycj??zcj??1??????=Tbc?1?Twbj??1?Twbi??Tbc??????λ1?uci??λ1?vci??λ1?1??????
其中:
T = [ R t 0 T 1 ] \begin{aligned} T=\left[\begin{array}{ll} R & t \\ 0^{T} & 1 \end{array}\right] \end{aligned} T=[R0T?t1?]?

  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?=???zcj??xcj????ucj??zcj??ycj????vcj???????

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??=vtw?v˙tw?=atw?q˙?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]?(qwbt??abt??gw)δt2vjw?=viw?+t[i,j]?(qwbt??abt??gw)δ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??=qwbi???qbi?bt??
那么,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?Δt?21?gwΔt2+qwbi???t[i,j]?(qbi?bt??abt?)δt2vjw?=viw??gwΔt+qwb???t[i,j]?(qbz?bt??abt?)δtqwbj??=qwbi??t[i,j]?qbi?bt???[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} αbi?bj??βbi?bj??qbi?bj???=?t[i,j]?(qbi?bt??abt?)δt2=t[i,j]?(qbi?bt??abt?)δt=t[i,j]?qbi?bt???[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] ???????pwbj??vjw?qwbj??bja?bjg?????????=???????pwbi??+viw?Δt?21?gwΔt2+qwbi??αbi?bj??viw??gwΔt+qwbi??βbi?bj??qwbq??qbi?bj??bia?big?????????

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] ???????rp?rq?rv?rba?rbg?????????15×1?=???????qbi?w?(pwbj???pwbi???viw?Δt+21?gwΔt2)?αbi?bj??2[qbj?bi???(qbi?w??qwbj??)]xyz?qbi?w?(vjw??viw?+gwΔt)?βbi?bj??bja??bia?bjg??big?????????
上面误差中位移,速度,偏置都是直接相减得到。第二顶是关于四元数的旋转误差,其中 [.] 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} ωqbi?bk+1??aαbi?bk+1??βbi?bk+1??bk+1a?bk+1g??=21?((ωbk??bkg?)+(ωbk+1??bkg?))=qbi?bk???[121?ωδt?]=21?(qbi?bk??(abk??bka?)+qbi?bk+1??(abk+1??bka?))=αbi?bk??+βbi?bk??δt+21?aδt2=βbi?bk??+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Σx?A?
Σ 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(Axx?A?)=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?=Fk?1?ηik?1?+Gk?1?nk?1?
比如: 状态量误差为 η 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?=Fk?1?Σik?1?Fk?1??+Gk?1?Σn?Gk?1T?
其中, Σ 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(xk?1?,uk?1?), 其中状态量为 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(xk?1?,uk?1?) 的状态误差的线性递推关系如下:
δ x k = F δ x k ? 1 + G n k ? 1 \delta x_{k}=F \delta x_{k-1}+G n_{k-1} δxk?=Fδxk?1?+Gnk?1?
其中, F \mathrm{F} F 是状态量 x k x_{k} xk? 对状态量 x k ? 1 x_{k-1} xk?1? 的雅克比矩阵, G \mathrm{G} G 是状态量 x k x_{k} xk? 对输入量 u k ? 1 u_{k-1} uk?1? 的雅克比矩阵。
证明:对非线性状态方程进行一阶泰勒展开有:
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(xk?1?,uk?1?)x^k?+δxk?=f(x^k?1?+δxk?1?,u^k?1?+nk?1?)x^k?+δxk?=f(x^k?1?,u^k?1?)+Fδxk?1?+Gnk?1??

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?=δxk?1?+δxk?1?Δtδxk?=(I+AΔt)δxk?1?+BΔtnk?1??
这两种推导方式的可以看出有:
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δab?R[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} ωqbi?bk+1??aαbi?bk+1??βbi?bk+1??bk+1a?bk+1g??=21?((ωbk?+nkg??bkg?)+(ωbk+1?+nk+1g??bkg?))=qbi?bk???[121?ωδt?]=21?(qbi?bk??(abk?+nka??bka?)+qbi?bk+1??(abk+1?+nk+1a??bka?))=αbi?bk??+βbi?bk??δt+21?aδt2=βbi?bk??+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δxk?1?+Gnk?1?, 我们?望能推导出如下的形式:
[ δ α 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+1?bk+1??δθbk+1?bk+1??δβbk+1?bk+1??δbk+1a?δbk+1g?????????=F???????δαbk?bk??δθbk?bk??δβbk?bk??δbka?δbkg?????????+G?????????nka?nkg?nk+1a?nk+1g?nbka??nbkg????????????
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=???????I0000?f12?I?[ω]×?f32?00?Iδt0I00??41?(qbi?bk??+qbi?bk+1??)δt20?21?(qbi?bk??+qbi?bk+1??)δtI0?f15??Iδtf35?0I????????
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=???????41?qbi?bk??δt2021?qbi?bk??δt00?g12?21?Iδtg32?00?41?qbi?bk+1??δt2021?qbi?bk+1??δt00?g14?21?Iδtg34?00?000Iδt0?0000Iδ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] ?????xcj??ycj??zcj??1??????=Tbc?1?Twbj??1?Twbi??Tbc??????λ1?uci??λ1?vci??λ1?1??????
拆成三维坐标形式为:
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??=???xcj??ycj??zcj???????=Rbc??Rwbj???Rwbi??Rbc?λ1????uci??vci??1????+Rbc??(Rwbj???((Rwbi??pbc?+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} fbi??fw?fbj???=Rbc?fci??+pbc?=Rwbi??fbi??+pwbi??=Rwbj???(fw??pwbj??)?
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=[?[δpbi?bi??δθbi?bi???]?rc???[pbj?bj??δθbj?bj???]?rc???[δpcc?δθcc??]?rc???δλ?rc???]

参考文献

VIO 中残差雅可比的推导

  数据结构与算法 最新文章
【力扣106】 从中序与后续遍历序列构造二叉
leetcode 322 零钱兑换
哈希的应用:海量数据处理
动态规划|最短Hamilton路径
华为机试_HJ41 称砝码【中等】【menset】【
【C与数据结构】——寒假提高每日练习Day1
基础算法——堆排序
2023王道数据结构线性表--单链表课后习题部
LeetCode 之 反转链表的一部分
【题解】lintcode必刷50题<有效的括号序列
上一篇文章      下一篇文章      查看所有文章
加:2022-04-22 19:01:27  更:2022-04-22 19:02:59 
 
开发: C++知识库 Java知识库 JavaScript Python PHP知识库 人工智能 区块链 大数据 移动开发 嵌入式 开发工具 数据结构与算法 开发测试 游戏开发 网络协议 系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑 笔记本 显卡 显示器 固态硬盘 硬盘 耳机 手机 iphone vivo oppo 小米 华为 单反 装机 图拉丁

360图书馆 购物 三丰科技 阅读网 日历 万年历 2024年11日历 -2024/11/26 7:31:00-

图片自动播放器
↓图片自动播放器↓
TxT小说阅读器
↓语音阅读,小说下载,古典文学↓
一键清除垃圾
↓轻轻一点,清除系统垃圾↓
图片批量下载器
↓批量下载图片,美女图库↓
  网站联系: qq:121756557 email:121756557@qq.com  IT数码