多传感器融合定位十一-基于滤波的融合方法Ⅱ

news2024/11/27 12:52:54

多传感器融合定位十一-基于滤波的融合方法Ⅱ

  • 1. 编码器运动模型及标定
    • 1.1 编码器基础知识
    • 1.2 编码器运动模型
      • 1.2.1 旋转半径求解
      • 1.2.2 角速度求解
      • 1.2.3 线速度求解
      • 1.2.4 位姿求解
    • 1.3 编码器的标定
      • 1.3.1 轮子半径标定
      • 1.3.2 轮子与底盘中心距离标定
  • 2. 融合编码器的滤波方法
    • 2.1 核心思路
    • 2.2 观测量定义
    • 2.3 观测方程推导
  • 3. 融合运动约束的滤波方法
  • 4. 融合点云特征的滤波方法
    • 4.1 整体思路
    • 4.2 滤波模型
    • 4.3 位姿更新
    • 4.4 相似工作

Reference:

  1. 深蓝学院-多传感器融合
  2. 多传感器融合定位理论基础

文章跳转:

  1. 多传感器融合定位一-3D激光里程计其一:ICP
  2. 多传感器融合定位二-3D激光里程计其二:NDT
  3. 多传感器融合定位三-3D激光里程计其三:点云畸变补偿
  4. 多传感器融合定位四-3D激光里程计其四:点云线面特征提取
  5. 多传感器融合定位五-点云地图构建及定位
  6. 多传感器融合定位六-惯性导航原理及误差分析
  7. 多传感器融合定位七-惯性导航解算及误差分析其一
  8. 多传感器融合定位八-惯性导航解算及误差分析其二
  9. 多传感器融合定位九-基于滤波的融合方法Ⅰ其一
  10. 多传感器融合定位十-基于滤波的融合方法Ⅰ其二
  11. 多传感器融合定位十一-基于滤波的融合方法Ⅱ
  12. 多传感器融合定位十二-基于图优化的建图方法其一
  13. 多传感器融合定位十三-基于图优化的建图方法其二
  14. 多传感器融合定位十四-基于图优化的定位方法
  15. 多传感器融合定位十五-多传感器时空标定(综述)

1. 编码器运动模型及标定

1.1 编码器基础知识

在这里插入图片描述
编码器感应轮子的旋转,并在旋转时输出脉冲,脉冲数与转过的角度呈线性比例关系

脉冲对应的是角度增量,有时也用增量除以时间,转成轮子转动的角速度输出。

需要注意的是,编码器只是各种转角测量方式中的一种,其他还有轮速计、霍尔传感器等,本课程以编码器为例子讲解模型,但同样适用于其他形式的传感器。

编码器安装方式有单轮、双轮、三轮,本课程推导仅围绕双轮差分模型进行展开。
该模型中,需要用到以下变量:

  • r L \boldsymbol{r}_L rL:左轮半径
  • r R \boldsymbol{r}_R rR:右轮半径
  • d \boldsymbol{d} d:轮子离底盘中心的距离
  • ω L \boldsymbol{\omega}_L ωL:左轮自转角速度
  • ω R \boldsymbol{\omega_R} ωR:右轮自转角速度
  • v L \boldsymbol{v}_L vL:左轮线速度
  • v R \boldsymbol{v}_R vR:右轮线速度

实际使用时,标定完成后, r L \boldsymbol{r}_L rL r R \boldsymbol{r}_R rR d \boldsymbol{d} d 为固定参数, ω L \boldsymbol{\omega}_L ωL ω R \boldsymbol{\omega}_R ωR 为测量值,而 v L \boldsymbol{v}_L vL v R \boldsymbol{v}_R vR 可以通过下式计算得到:
v L = ω L r L v R = ω R r R \begin{aligned} \boldsymbol{v}_L & =\boldsymbol{\omega}_L \boldsymbol{r}_L \\ \boldsymbol{v}_R & =\boldsymbol{\omega}_R \boldsymbol{r}_R \end{aligned} vLvR=ωLrL=ωRrR

1.2 编码器运动模型

在这里插入图片描述

运动模型的作用是,使用前述已知量,求解以下变量:
ω \boldsymbol{\omega} ω : 底盘中心的角速度
v \boldsymbol{v} v : 底盘中心的线速度
r \boldsymbol{r} r : 底盘中心圆弧运动旋转半径

1.2.1 旋转半径求解

双轮差分模型下,左右轮圆弧运动的角速度相等,且等于底盘中心圆弧运动的角速度(两个轮子的自转角速度是相同的),因此有:
ω = v L r − d = v R r + d \boldsymbol{\omega}=\frac{\boldsymbol{v}_L}{\boldsymbol{r}-\boldsymbol{d}}=\frac{\boldsymbol{v}_R}{\boldsymbol{r}+\boldsymbol{d}} ω=rdvL=r+dvR由此可以得出:
v L ( r + d ) = v R ( r − d ) \boldsymbol{v}_L(\boldsymbol{r}+\boldsymbol{d})=\boldsymbol{v}_R(\boldsymbol{r}-\boldsymbol{d}) vL(r+d)=vR(rd)移项可得:
( v R − v L ) r = ( v R + v L ) d \left(\boldsymbol{v}_R-\boldsymbol{v}_L\right) \boldsymbol{r}=\left(\boldsymbol{v}_R+\boldsymbol{v}_L\right) \boldsymbol{d} (vRvL)r=(vR+vL)d从而可以得到:
r = v R + v L v R − v L d \boldsymbol{r}=\frac{\boldsymbol{v}_R+\boldsymbol{v}_L}{\boldsymbol{v}_R-\boldsymbol{v}_L} \boldsymbol{d} r=vRvLvR+vLd

1.2.2 角速度求解

把旋转半径的求解结果,代入角速度公式,即可得到:
ω = v L v R + v L v R − v L d − d = v R − v L 2 d \boldsymbol{\omega}=\frac{\boldsymbol{v}_L}{\frac{\boldsymbol{v}_R+\boldsymbol{v}_L}{\boldsymbol{v}_R-\boldsymbol{v}_L} \boldsymbol{d}-\boldsymbol{d}}=\frac{\boldsymbol{v}_R-\boldsymbol{v}_L}{2 \boldsymbol{d}} ω=vRvLvR+vLddvL=2dvRvL

1.2.3 线速度求解

利用旋转角速度和旋转半径的结果,可以直接得到线速度:
v = ω r = v R − v L 2 d v R + v L v R − v L d = v R + v L 2 \boldsymbol{v}=\boldsymbol{\omega} \boldsymbol{r}=\frac{\boldsymbol{v}_R-\boldsymbol{v}_L}{2 \boldsymbol{d}} \frac{\boldsymbol{v}_R+\boldsymbol{v}_L}{\boldsymbol{v}_R-\boldsymbol{v}_L} \boldsymbol{d}=\frac{\boldsymbol{v}_R+\boldsymbol{v}_L}{2} v=ωr=2dvRvLvRvLvR+vLd=2vR+vL

1.2.4 位姿求解

在这里插入图片描述
假设 x k , y k , θ k \boldsymbol{x}_k, \boldsymbol{y}_k, \boldsymbol{\theta}_k xk,yk,θk 为当前时刻位姿, x k − 1 , y k − 1 , θ k − 1 \boldsymbol{x}_{k-1}, \boldsymbol{y}_{k-1}, \boldsymbol{\theta}_{k-1} xk1,yk1,θk1 为上一时刻的位姿,则有:
θ k = θ k − 1 + ω Δ t x k = x k − 1 + v Δ t cos ⁡ ( θ k − 1 ) y k = y k − 1 + v Δ t sin ⁡ ( θ k − 1 ) \begin{aligned} & \boldsymbol{\theta}_k=\boldsymbol{\theta}_{k-1}+\boldsymbol{\omega} \Delta t \\ & \boldsymbol{x}_k=\boldsymbol{x}_{k-1}+\boldsymbol{v} \Delta t \cos \left(\boldsymbol{\theta}_{\boldsymbol{k - 1}}\right) \\ & \boldsymbol{y}_k=\boldsymbol{y}_{k-1}+\boldsymbol{v} \Delta t \sin \left(\boldsymbol{\theta}_{\boldsymbol{k}-1}\right) \end{aligned} θk=θk1+ωΔtxk=xk1+vΔtcos(θk1)yk=yk1+vΔtsin(θk1)其中:
Δ t = t k − t k − 1 \Delta t=t_k-t_{k-1} Δt=tktk1

1.3 编码器的标定

标定可以理解为运动模型求解过程的反向过程,具体是指在已知底盘中心线速度、角速度的情况下,求解轮子半径、轮子离底盘中心距离等。
已知量:
v v v:底盘中心的线速度
ω \omega ω:底盘中心的角速度
待求解量:
r L \boldsymbol{r}_L rL:左轮半径
r R \boldsymbol{r}_R rR:右轮半径
d \boldsymbol{d} d:轮子离底盘中心的距离
实际标定时,线速度、角速度由其他传感器提供(比如雷达点云和地图匹配),且为了简化模型,认为雷达装在底盘中心正上方。(这里雷达最好的方法是先建好一个点云地图,然后在点云地图里面匹配,然后拿它做观测,去提供线速度和角速度的结果,而不要用激光里程计-----里程计本身就是有累计误差的)

1.3.1 轮子半径标定

由于速度的求解公式为:
v = v R + v L 2 = ω R r R + ω L r L 2 \boldsymbol{v}=\frac{\boldsymbol{v}_R+\boldsymbol{v}_L}{2}=\frac{\boldsymbol{\omega}_R \boldsymbol{r}_R+\boldsymbol{\omega}_L \boldsymbol{r}_L}{2} v=2vR+vL=2ωRrR+ωLrL它可以重新写为:
[ ω R ω L ] [ r R r L ] = 2 v \left[\begin{array}{ll} \boldsymbol{\omega}_R & \boldsymbol{\omega}_L \end{array}\right]\left[\begin{array}{l} \boldsymbol{r}_R \\ \boldsymbol{r}_L \end{array}\right]=2 \boldsymbol{v} [ωRωL][rRrL]=2v当有多组测量值时,可以构成如下方程组:
[ ω R 0 ω L 0 ω R 1 ω L 1 ⋮ ⋮ ω R N ω L N ] [ r R r L ] = [ 2 v 0 2 v 1 ⋮ 2 v N ] \left[\begin{array}{cc} \boldsymbol{\omega}_{R 0} & \boldsymbol{\omega}_{L 0} \\ \boldsymbol{\omega}_{R 1} & \boldsymbol{\omega}_{L 1} \\ \vdots & \vdots \\ \boldsymbol{\omega}_{R N} & \boldsymbol{\omega}_{L N} \end{array}\right]\left[\begin{array}{l} \boldsymbol{r}_R \\ \boldsymbol{r}_L \end{array}\right]=\left[\begin{array}{c} 2 \boldsymbol{v}_0 \\ 2 \boldsymbol{v}_1 \\ \vdots \\ 2 \boldsymbol{v}_N \end{array}\right] ωR0ωR1ωRNωL0ωL1ωLN [rRrL]= 2v02v12vN 这是典型的最小二乘问题,可用最小二乘标准形式计算。

1.3.2 轮子与底盘中心距离标定

由于角速度的求解公式为:
ω = v R − v L 2 d \boldsymbol{\omega}=\frac{\boldsymbol{v}_R-\boldsymbol{v_L}}{2 \boldsymbol{d}} ω=2dvRvL在经过轮子半径标定之后,分子上的两项可认为是已知量,因此可以得到:
d = v R − v L 2 ω \boldsymbol{d}=\frac{\boldsymbol{v}_R-\boldsymbol{v}_L}{2 \boldsymbol{\omega}} d=2ωvRvL虽然可直接求解,但是为了抑制噪声带来的影响,因多次采样计算取平均。

2. 融合编码器的滤波方法

2.1 核心思路

在上一节课滤波模型的基础上增加编码器进行融合,有一种非常简单的方法,即使用编码器解算的速度作为观测量,加入原来模型的观测方程中,而其他环节保持不变。

2.2 观测量定义

编码器提供的是载体系下的速度观测,在前 ( x ) (\mathrm{x}) (x)-左 ( y ) (\mathrm{y}) (y)-上(z)坐标系的定义下, x \mathrm{x} x 方向的速度分量是已知的 v x b = v m \boldsymbol{v}_x^b=\boldsymbol{v}_m vxb=vm。 另外,在以车作为载体的情况下,由于车的侧向和天向没有运动,因此又有 v y b = 0 \boldsymbol{v}_y^b=0 vyb=0 v z b = 0 \boldsymbol{v}_z^b=0 vzb=0
基于此,我们可以认为 b b b 系 3 个维度的速度分量都是可观测的。

2.3 观测方程推导

由于导航解算得到的是 w w w 系下得速度,而速度观测是 b b b 系下得,因此需要推导二者之间的误差关系,才能得到相应的观测方程。
推导方法仍按照第6讲的固定套路进行。

  1. 写出不考虑误差时的方程:
    v b = R b w v w \boldsymbol{v}^b=\boldsymbol{R}_{b w} \boldsymbol{v}^w vb=Rbwvw
  2. 写出考虑误差时的方程:
    v ~ b = R ~ b w v ~ w \tilde{\boldsymbol{v}}^b=\tilde{\boldsymbol{R}}_{b w} \tilde{\boldsymbol{v}}^w v~b=R~bwv~w
  3. 写出真实值与理想值之间的关系:
    v ~ b = v b + δ v b v ~ w = v w + δ v w R ~ b w = R ~ w b T = ( R w b ( I + [ δ θ ] × ) ) T = ( I − [ δ θ ] × ) R b w \begin{aligned} & \tilde{\boldsymbol{v}}^b=\boldsymbol{v}^b+\delta \boldsymbol{v}^b \\ & \tilde{\boldsymbol{v}}^w=\boldsymbol{v}^w+\delta \boldsymbol{v}^w \\ & \tilde{\boldsymbol{R}}_{b w}=\tilde{\boldsymbol{R}}_{w b}^T=\left(\boldsymbol{R}_{w b}\left(\boldsymbol{I}+[\delta \boldsymbol{\theta}]_{\times}\right)\right)^T \\ & =\left(\boldsymbol{I}-[\delta \boldsymbol{\theta}]_{\times}\right) \boldsymbol{R}_{b w} \end{aligned} v~b=vb+δvbv~w=vw+δvwR~bw=R~wbT=(Rwb(I+[δθ]×))T=(I[δθ]×)Rbw
  4. 把3)中的关系带入2)式:
    v b + δ v b = ( I − [ δ θ ] × ) R b w ( v w + δ v w ) \boldsymbol{v}^b+\delta \boldsymbol{v}^b=\left(\boldsymbol{I}-[\delta \boldsymbol{\theta}]_{\times}\right) \boldsymbol{R}_{b w}\left(\boldsymbol{v}^w+\delta \boldsymbol{v}^w\right) vb+δvb=(I[δθ]×)Rbw(vw+δvw)
  5. 把1)中的关系带入4)式:
    R b w v w + δ v b = ( I − [ δ θ ] × ) R b w ( v w + δ v w ) \boldsymbol{R}_{b w} \boldsymbol{v}^w+\delta \boldsymbol{v}^b=\left(\boldsymbol{I}-[\delta \boldsymbol{\theta}]_{\times}\right) \boldsymbol{R}_{b w}\left(\boldsymbol{v}^w+\delta \boldsymbol{v}^w\right) Rbwvw+δvb=(I[δθ]×)Rbw(vw+δvw)
  6. 化简方程:
    δ v b = R b w δ v w − [ δ θ ] × R b w v w = R b w δ v w − [ δ θ ] × v b = R b w δ v w + [ v b ] × δ θ \begin{aligned} \delta \boldsymbol{v}^b & =\boldsymbol{R}_{b w} \delta \boldsymbol{v}^w-[\delta \boldsymbol{\theta}]_{\times} \boldsymbol{R}_{b w} \boldsymbol{v}^w \\ & =\boldsymbol{R}_{b w} \delta \boldsymbol{v}^w-[\delta \boldsymbol{\theta}]_{\times} \boldsymbol{v}^b \\ & =\boldsymbol{R}_{b w} \delta \boldsymbol{v}^w+\left[\boldsymbol{v}^b\right]_{\times} \delta \boldsymbol{\theta} \end{aligned} δvb=Rbwδvw[δθ]×Rbwvw=Rbwδvw[δθ]×vb=Rbwδvw+[vb]×δθ

根据前一章内容,状态量为
δ x = [ δ p δ v δ θ δ b a δ b ω ] \delta \boldsymbol{x}=\left[\begin{array}{c} \delta \boldsymbol{p} \\ \delta \boldsymbol{v} \\ \delta \boldsymbol{\theta} \\ \delta \boldsymbol{b}_a \\ \delta \boldsymbol{b}_\omega \end{array}\right] δx= δpδvδθδbaδbω 而融合编码器以后,观测量变为
y = [ δ p ‾ δ v ‾ b δ θ ‾ ] \boldsymbol{y}=\left[\begin{array}{c} \delta \overline{\boldsymbol{p}} \\ \delta \overline{\boldsymbol{v}}^b \\ \delta \overline{\boldsymbol{\theta}} \end{array}\right] y= δpδvbδθ 其中 δ v ‾ b \delta \overline{\boldsymbol{v}}^b δvb 的观测值可以通过下式获得
δ v ‾ b = v ~ b − v b = R ~ b w v ~ w − [ v m 0 0 ] \delta \overline{\boldsymbol{v}}_b=\tilde{\boldsymbol{v}}^b-\boldsymbol{v}^b=\tilde{\boldsymbol{R}}_{b w} \tilde{\boldsymbol{v}}^w-\left[\begin{array}{c} \boldsymbol{v}_m \\ 0 \\ 0 \end{array}\right] δvb=v~bvb=R~bwv~w vm00 此时的观测方程 y = G t δ x + C t n \boldsymbol{y}=\boldsymbol{G}_t \delta \boldsymbol{x}+\boldsymbol{C}_t \boldsymbol{n} y=Gtδx+Ctn 中的各变量应重新写为
G t = [ I 3 0 0 0 0 0 R b w [ v b ] × 0 0 0 0 I 3 0 0 ] C t = [ I 3 0 0 0 I 3 0 0 0 I 3 ] n = [ n δ p ˉ x n δ p ˉ y n δ p ˉ z n δ v ˉ x b n δ v ˉ y b n δ v ˉ z b n δ θ ˉ x n δ θ ˉ y n δ θ ˉ z ] T \begin{aligned} & \boldsymbol{G}_t=\left[\begin{array}{ccccc} \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{R}_{b w} & {\left[\boldsymbol{v}^b\right]_{\times}} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} \end{array}\right] \\ & \boldsymbol{C}_t=\left[\begin{array}{ccc} \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{I}_3 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \boldsymbol{I}_3 \end{array}\right] \\ & \boldsymbol{n}=\left[\begin{array}{lllllllll} n_{\delta \bar{p}_x} & n_{\delta \bar{p}_y} & n_{\delta \bar{p}_z} & n_{\delta \bar{v}_x^b} & n_{\delta \bar{v}_y^b} & n_{\delta \bar{v}_z^b} & n_{\delta \bar{\theta}_x} & n_{\delta \bar{\theta}_y} & n_{\delta \bar{\theta}_z} \end{array}\right]^T \end{aligned} Gt= I3000Rbw00[vb]×I3000000 Ct= I3000I3000I3 n=[nδpˉxnδpˉynδpˉznδvˉxbnδvˉybnδvˉzbnδθˉxnδθˉynδθˉz]T随后,便可以使用新的观测方程,不改变其他方程,直接按照原有流程进行Kalman滤波融合。

3. 融合运动约束的滤波方法

很多时候,硬件平台并没有编码器,不能直接使用上一小节的模型,但是车本身的运动特性(即侧向速度和天向速度为0)仍然可以使用。

它对观测量带来的改变仅仅是少了一个维度( x x x 方向),而推导方法并没有改变,因此此处直接给出该融合模式下的推导结果。

新的观测量为
y = [ δ p ‾ [ δ v ‾ b ] y z δ θ ‾ ] \boldsymbol{y}=\left[\begin{array}{c} \delta \overline{\boldsymbol{p}} \\ {\left[\delta \overline{\boldsymbol{v}}^b\right]_{y z}} \\ \delta \overline{\boldsymbol{\theta}} \end{array}\right] y= δp[δvb]yzδθ [ ∙ ] y z [\bullet]_{y z} []yz 表示只取三维向量或矩阵的后2行

此时的观测方程 y = G t δ x + C t n \boldsymbol{y}=\boldsymbol{G}_t \delta \boldsymbol{x}+\boldsymbol{C}_t \boldsymbol{n} y=Gtδx+Ctn 中的各变量应重新写为
G t = [ I 3 0 0 0 0 0 [ R b w ] y z [ [ v b ] × ] y z 0 0 0 0 I 3 0 0 ] C t = [ I 3 0 0 0 I 2 0 0 0 I 3 ] n = [ n δ p ˉ x n δ p ˉ y n δ p ˉ z n δ v ˉ y b n δ v ˉ z b n δ θ ˉ x n δ θ ˉ y n δ θ ˉ z ] T \begin{aligned} & \boldsymbol{G}_t=\left[\begin{array}{ccccc} \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & {\left[\boldsymbol{R}_{b w}\right]_{y z}} & {\left[\left[\boldsymbol{v}^b\right]_{\times}\right]_{y z}} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} \end{array}\right] \\ & \boldsymbol{C}_t=\left[\begin{array}{ccc} \boldsymbol{I}_3 & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{I}_2 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \boldsymbol{I}_3 \end{array}\right] \\ & \boldsymbol{n}=\left[\begin{array}{llllllll} n_{\delta \bar{p}_x} & n_{\delta \bar{p}_y} & n_{\delta \bar{p}_z} & n_{\delta \bar{v}_y^b} & n_{\delta \bar{v}_z^b} & n_{\delta \bar{\theta}_x} & n_{\delta \bar{\theta}_y} & n_{\delta \bar{\theta}_z} \end{array}\right]^T \end{aligned} Gt= I3000[Rbw]yz00[[vb]×]yzI3000000 Ct= I3000I2000I3 n=[nδpˉxnδpˉynδpˉznδvˉybnδvˉzbnδθˉxnδθˉynδθˉz]T随后的Kalman流程仍然与之前保持一致。

4. 融合点云特征的滤波方法

4.1 整体思路

以IMU做状态预测,以特征中的点-面距离、点-线距离为约束(观测),修正误差。
在这里插入图片描述
论文题目:LINS: A Lidar-Inertial State Estimator for Robust and Efficient Navigation

4.2 滤波模型

  1. 状态定义
    位姿定义:
    x w b k : = [ p w b k , q w b k ] \mathbf{x}_w^{b_k}:=\left[\mathbf{p}_w^{b_k}, \mathbf{q}_w^{b_k}\right] xwbk:=[pwbk,qwbk]相对位姿相关:
    x b k + 1 b k : = [ p b k + 1 b k , v b k + 1 b k , q b k + 1 b k , b a , b g , g b k ] \mathbf{x}_{b_{k+1}}^{b_k}:=\left[\mathbf{p}_{b_{k+1}}^{b_k}, \mathbf{v}_{b_{k+1}}^{b_k}, \mathbf{q}_{b_{k+1}}^{b_k}, \mathbf{b}_a, \mathbf{b}_g, \mathbf{g}^{b_k}\right] xbk+1bk:=[pbk+1bk,vbk+1bk,qbk+1bk,ba,bg,gbk]状态量:
    δ x : = [ δ p , δ v , δ θ , δ b a , δ b g , δ g ] \delta \mathbf{x}:=\left[\delta \mathbf{p}, \delta \mathbf{v}, \delta \boldsymbol{\theta}, \delta \mathbf{b}_a, \delta \mathbf{b}_g, \delta \mathbf{g}\right] δx:=[δp,δv,δθ,δba,δbg,δg]状态量修正:
    x b k + 1 b k = − x b k + 1 b k ⊞ δ x = [ − p b k + 1 b k + δ p − v k b k + δ v − q b k + 1 b k ⊗ exp ⁡ ( δ θ ) − b a + δ b a − b g + δ b g − g b k + δ g ] \mathbf{x}_{b_{k+1}}^{b_k}={ }^{-} \mathbf{x}_{b_{k+1}}^{b_k} \boxplus \boldsymbol{\delta} \mathbf{x}=\left[\begin{array}{c} -\mathbf{p}_{b_{k+1}}^{b_k}+\boldsymbol{\delta} \mathbf{p} \\ -\mathbf{v}_k^{b_k}+\boldsymbol{\delta} \mathbf{v} \\ -\mathbf{q}_{b_{k+1}}^{b_k} \otimes \exp (\boldsymbol{\delta} \boldsymbol{\theta}) \\ -\mathbf{b}_a+\boldsymbol{\delta} \mathbf{b}_a \\ -\mathbf{b}_g+\delta \mathbf{b}_g \\ -\mathbf{g}^{b_k}+\boldsymbol{\delta} \mathbf{g} \end{array}\right] xbk+1bk=xbk+1bkδx= pbk+1bk+δpvkbk+δvqbk+1bkexp(δθ)ba+δbabg+δbggbk+δg 其中 − x b k + 1 b k { }^{-} \mathbf{x}_{b_{k+1}}^{b_k} xbk+1bk 表示预测值。

  2. 状态方程
    δ x ˙ ( t ) = F t δ x ( t ) + G t w \delta \dot{\mathbf{x}}(t)=\mathbf{F}_t \delta \mathbf{x}(t)+\mathbf{G}_t \mathbf{w} δx˙(t)=Ftδx(t)+Gtw其中
    F t = [ 0 I 0 0 0 0 0 0 − R t b k [ a ^ t ] × − R t b k 0 − I 3 0 0 − [ ω ^ t ] × 0 − I 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ] G t = [ 0 0 0 0 − R t b k 0 0 0 0 − I 3 0 0 0 0 I 3 0 0 0 0 I 3 0 0 0 0 ] w = [ n a T , n g T , n b a T , n b g T ] T \begin{aligned} \mathbf{F}_t & =\left[\begin{array}{cccccc} 0 & \mathbf{I} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & -\mathbf{R}_t^{b_k}\left[\hat{\mathbf{a}}_t\right]_{\times} & -\mathbf{R}_t^{b_k} & \mathbf{0} & -\mathbf{I}_3 \\ 0 & 0 & -\left[\hat{\omega}_t\right]_{\times} & \mathbf{0} & -\mathbf{I}_3 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \end{array}\right] \quad \mathbf{G}_t=\left[\begin{array}{cccc} \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ -\mathbf{R}_t^{b_k} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & -\mathbf{I}_3 & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{I}_3 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I}_3 \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \end{array}\right] \\ \mathbf{w} & =\left[\mathbf{n}_a^T, \mathbf{n}_g^T, \mathbf{n}_{b_a}^T, \mathbf{n}_{b_g}^T\right]^T \end{aligned} Ftw= 000000I000000Rtbk[a^t]×[ω^t]×0000Rtbk000000I30000I30000 Gt= 0Rtbk000000I3000000I3000000I30 =[naT,ngT,nbaT,nbgT]T

  3. 观测方程
    观测的计算与loam中前后帧匹配的思想一致,都是计算 点-面、点-线的残差
    其中
    p ^ i l k = ( R l b ) T ( R b k + 1 b k ( R l b p i l k + 1 + p l b ) + p b k + 1 b k − p l b ) \hat{\mathbf{p}}_i^{l_k}=\left(\mathbf{R}_l^b\right)^T\left(\mathbf{R}_{b_{k+1}}^{b_k}\left(\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right)+\mathbf{p}_{b_{k+1}}^{b_k}-\mathbf{p}_l^b\right) p^ilk=(Rlb)T(Rbk+1bk(Rlbpilk+1+plb)+pbk+1bkplb)为了计算观测方程,需要计算残差对状态量的雅可比
    H k = ∂ f ∂ p ^ i l k ⋅ ∂ p ^ i l k ∂ δ x \mathbf{H}_k=\frac{\partial \mathbf{f}}{\partial \hat{\mathbf{p}}_i^{l_k}} \cdot \frac{\partial \hat{\mathbf{p}}_i^{l_k}}{\partial \delta \mathbf{x}} Hk=p^ilkfδxp^ilk上式中包含两部分,第一部分已经讲过,此处只推导第 二部分。
    除了旋转和平移外,残差项对其它量的导数均为 0 。
    a.对平移求导
    ∂ p ^ i l k ∂ δ p = ∂ [ ( R l b ) ⊤ ( R b k + 1 b k ( R l b p i l k + 1 + p l b ) + p b k + 1 b k − p l b ) ] ∂ δ p = ∂ [ ( R l b ) ⊤ p b k + 1 b k ] ∂ δ p = ( R l b ) ⊤ \begin{aligned} & \frac{\partial \hat{\mathbf{p}}_i^{l_k}}{\partial \delta \mathbf{p}} \\ = & \frac{\partial\left[\left(\mathbf{R}_l^b\right)^{\top}\left(\mathbf{R}_{b_{k+1}}^{b_k}\left(\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right)+\mathbf{p}_{b_{k+1}}^{b_k}-\mathbf{p}_l^b\right)\right]}{\partial \delta \mathbf{p}} \\ = & \frac{\partial\left[\left(\mathbf{R}_l^b\right)^{\top} \mathbf{p}_{b_{k+1}}^{b_k}\right]}{\partial \delta \mathbf{p}} \\ = & \left(\mathbf{R}_l^b\right)^{\top} \end{aligned} ===δpp^ilkδp[(Rlb)(Rbk+1bk(Rlbpilk+1+plb)+pbk+1bkplb)]δp[(Rlb)pbk+1bk](Rlb)b.对旋转求导
    ∂ p ^ i l k ∂ δ x = ∂ [ ( R l b ) T ( R b k + 1 b k ( R l b p i l k + 1 + p l b ) + p b k + 1 b k − p l b ) ] ∂ δ θ = ( R l b ) T ∂ [ R b k + 1 b k ( R l b p i l k + 1 + p l b ) ] ∂ R b k + 1 b k ∂ R b k + 1 b k ∂ δ θ = − ( R l b ) T [ R b k + 1 b k ( R l b p i l k + 1 + p l b ) ] × R b k + 1 b k J r − 1 ( θ ) = − ( R l b ) T R b k + 1 b k [ R l b p i l k + 1 + p l b ] × R b k b k + 1 R b k + 1 b k J r − 1 ( θ ) = − ( R l b ) T R b k + 1 b k [ R l b p i l k + 1 + p l b ] × J r − 1 ( θ ) \begin{aligned} & \frac{\partial \hat{\mathbf{p}}_i^{l_k}}{\partial \delta \mathbf{x}} \\ = & \frac{\partial\left[\left(\mathbf{R}_l^b\right)^T\left(\mathbf{R}_{b_{k+1}}^{b_k}\left(\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right)+\mathbf{p}_{b_{k+1}}^{b_k}-\mathbf{p}_l^b\right)\right]}{\partial \delta \boldsymbol{\theta}} \\ = & \left(\mathbf{R}_l^b\right)^T \frac{\partial\left[\mathbf{R}_{b_{k+1}}^{b_k}\left(\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right)\right]}{\partial \mathbf{R}_{b_{k+1}}^{b_k}} \frac{\partial \mathbf{R}_{b_{k+1}}^{b_k}}{\partial \delta \boldsymbol{\theta}} \\ = & -\left(\mathbf{R}_l^b\right)^T\left[\mathbf{R}_{b_{k+1}}^{b_k}\left(\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right)\right]_{\times} \mathbf{R}_{b_{k+1}}^{b_k} \boldsymbol{J}_r^{-1}(\boldsymbol{\theta}) \\ = & -\left(\mathbf{R}_l^b\right)^T \mathbf{R}_{b_{k+1}}^{b_k}\left[\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right]_{\times} \mathbf{R}_{b_k}^{b_{k+1}} \mathbf{R}_{b_{k+1}}^{b_k} \boldsymbol{J}_r^{-1}(\boldsymbol{\theta}) \\ = & -\left(\mathbf{R}_l^b\right)^T \mathbf{R}_{b_{k+1}}^{b_k}\left[\mathbf{R}_l^b \mathbf{p}_i^{l_{k+1}}+\mathbf{p}_l^b\right]_{\times} \boldsymbol{J}_r^{-1}(\boldsymbol{\theta}) \end{aligned} =====δxp^ilkδθ[(Rlb)T(Rbk+1bk(Rlbpilk+1+plb)+pbk+1bkplb)](Rlb)TRbk+1bk[Rbk+1bk(Rlbpilk+1+plb)]δθRbk+1bk(Rlb)T[Rbk+1bk(Rlbpilk+1+plb)]×Rbk+1bkJr1(θ)(Rlb)TRbk+1bk[Rlbpilk+1+plb]×Rbkbk+1Rbk+1bkJr1(θ)(Rlb)TRbk+1bk[Rlbpilk+1+plb]×Jr1(θ)预测 :
    δ x t τ = ( I + F t τ Δ t ) δ x t τ − 1 P t τ = ( I + F t τ Δ t ) P t τ − 1 ( I + F t τ Δ t ) T + ( G i τ Δ t ) Q ( G t t Δ t ) T K k , j = P k H k , j T ( H k , j P k H k , j T + J k , j M k J k , j T ) − 1 \begin{aligned} & \delta \mathbf{x}_{t_\tau}=\left(\mathbf{I}+\mathbf{F}_{t_\tau} \Delta t\right) \delta \mathbf{x}_{t_{\tau-1}} \\ & \mathbf{P}_{t_\tau}=\left(\mathbf{I}+\mathbf{F}_{t_\tau} \Delta t\right) \mathbf{P}_{t_{\tau-1}}\left(\mathbf{I}+\mathbf{F}_{t_\tau} \Delta t\right)^T+\left(\mathbf{G}_{i_\tau} \Delta t\right) \mathbf{Q}\left(\mathbf{G}_{t_t} \Delta t\right)^T \\ & \mathbf{K}_{k, j}=\mathbf{P}_k \mathbf{H}_{k, j}^T\left(\mathbf{H}_{k, j} \mathbf{P}_k \mathbf{H}_{k, j}^T+\mathbf{J}_{k, j} \mathbf{M}_k \mathbf{J}_{k, j}^T\right)^{-1} \end{aligned} δxtτ=(I+FtτΔt)δxtτ1Ptτ=(I+FtτΔt)Ptτ1(I+FtτΔt)T+(GiτΔt)Q(GttΔt)TKk,j=PkHk,jT(Hk,jPkHk,jT+Jk,jMkJk,jT)1迭代观测: Δ x j = K k , j ( H k , j δ x j − f ( − x b k + 1 b k ⊞ δ x j ) ) \Delta \mathbf{x}_j=\mathbf{K}_{k, j}\left(\mathbf{H}_{k, j} \delta \mathbf{x}_j-f\left({ }^{-} \mathbf{x}_{b_{k+1}}^{b_k} \boxplus \boldsymbol{\delta} \mathbf{x}_j\right)\right) Δxj=Kk,j(Hk,jδxjf(xbk+1bkδxj))
    δ x j + 1 = δ x j + Δ x j \delta \mathbf{x}_{j+1}=\delta \mathbf{x}_j+\Delta \mathbf{x}_j δxj+1=δxj+Δxj后验方差: P k + 1 = ( I − K k , n H k , n ) P k ( I − K k , n H k , n ) T + K k , n M k K k , n T \mathbf{P}_{k+1}=\left(\mathbf{I}-\mathbf{K}_{k, n} \mathbf{H}_{k, n}\right) \mathbf{P}_k\left(\mathbf{I}-\mathbf{K}_{k, n} \mathbf{H}_{k, n}\right)^T+\mathbf{K}_{k, n} \mathbf{M}_k \mathbf{K}_{k, n}^T Pk+1=(IKk,nHk,n)Pk(IKk,nHk,n)T+Kk,nMkKk,nT

4.3 位姿更新

x w b k + 1 = [ p w b k + 1 q w b k + 1 ] = [ R b k b k + 1 ( p w b k − p b k + 1 b k ) q b k b k + 1 ⊗ q w b k ] \mathbf{x}_w^{b_{k+1}}=\left[\begin{array}{c} \mathbf{p}_w^{b_{k+1}} \\ \mathbf{q}_w^{b_{k+1}} \end{array}\right]=\left[\begin{array}{c} \mathbf{R}_{b_k}^{b_{k+1}}\left(\mathbf{p}_w^{b_k}-\mathbf{p}_{b_{k+1}}^{b_k}\right) \\ \mathbf{q}_{b_k}^{b_{k+1}} \otimes \mathbf{q}_w^{b_k} \end{array}\right] xwbk+1=[pwbk+1qwbk+1]=[Rbkbk+1(pwbkpbk+1bk)qbkbk+1qwbk]

4.4 相似工作

论文题目:FAST-LIO: A Fast, Robust LiDAR-inertial Odometry Package by Tightly-Coupled Iterated Kalman Filter

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/335998.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

调用chatgpt的api, 必须知道的三件事

牙叔教程 简单易懂 调用api的代码 let url "https://api.openai.com/v1/completions"; let answer await axios // 使用axios发送post请求.post(url, data, { headers: headers }).then((res) > {return res.data.choices[0].text.trim();}).catch((err) >…

谈谈会话管理

客户端和服务器之间进行数据传输遵循的是HTTP协议, 此协议属于无状态协议(一次请求对应一次响应, 响应完之后断开连接), 服务器是无法跟踪客户端的请求, 通过cookie技术可以给客户端添加一个标识, 客户端之后发出的每次请求都会带着这个标识从而让服务器识别此客户端, 但由于co…

PostgreSQL入门

PostgreSQL入门 简介 PostgreSQL是以加州大学伯克利分校计算机系开发的POSTGRES, 版本 4.2为基础的对象关系型数据库管理系统(ORDBMS) 支持大部分SQL标准并且提供了许多现代特性 复杂查询外键触发器可更新视图事务完整性多版本并发控制 …

引导滤波code

文章目录1. 原理概述2. 实验环节2.1 验证与opencv 库函数的结果一致2.2 与 双边滤波比较2.3 引导滤波应用,fathering2.3 引导滤波应用,图像增强2.4 灰度图引导,和各自通道引导的效果差异2.5 不同参数设置影响3. 参考引导滤波1. 原理概述 引导…

VHDL语言基础-状态机设计-ASM图法状态机设计

目录 有限状态机的描述方法: ASM图: 状态转移图: 状态转移列表: MDS图: ASM图法状态机设计: ASM图的组成: 状态框: 判断框: 条件框: 状态框与条件框…

Python之FileNotFoundError: [Errno 2] No such file or directory问题处理

错误信息:FileNotFoundError: [Errno 2] No such file or directory: ../AutoFrame/temp/report.xlsx相对于当前文件夹的路径,其实就是你写的py文件所在的文件夹路径!python在对文件的操作时,需要特别注意文件地址的书写。文件的路…

上海亚商投顾:三大指数集体调整 消费板块逆市活跃

上海亚商投顾前言:无惧大盘涨跌,解密龙虎榜资金,跟踪一线游资和机构资金动向,识别短期热点和强势个股。市场情绪三大指数今日集体调整,沪指全天弱势震荡,创业板指盘中跌超1%。旅游、食品、乳业等大消费板块…

渗透测试 -- 网站信息收集

数据来源 01 网站指纹识别 网站的最基本组成:服务器(操作系统)、中间件(web容器)、脚本语言(php、java、...)、数据库(Mysql、...)为什么要了解这些? 举个例子:发现了一…

vue3组件库项目学习笔记(八):Git 使用总结

目前组件库的开发已经接近尾声,因为这次是使用 git 进行协作的开发模式,在团队协作的时候遇到很多的问题,开发过程中发现小伙伴们对于 git 的使用还不是很熟练,这里就简单总结一下常用的 git 的操作,大致有&#xff1a…

Revit快速材质切换:同一墙面赋予不同材质的方法

一、Revit中对同一墙面赋予不同材质的方法 方法1:零件法 重点:通过工作平面面板上的设置工作平面命令选取正确的面取消勾选通过原始分类的材质,如图1所示 方法2:拆分构造层绘制一道墙体,选择创建的墙体,单击…

判断元素是否在可视区域

前言 在日常开发中,我们经常需要判断目标元素是否在视窗之内或者和视窗的距离小于一个值(例如 100 px),从而实现一些常用的功能,例如: 图片的懒加载列表的无限滚动计算广告元素的曝光情况可点击链接的预加…

关于TL431和光耦PC817反馈控制部分电阻取值计算

计算R4和R3。TL431的R端流入电流2uA,为了保证取样精度,即不让TL431的R端吸取电流参与R3和R4的分压,可以设置R4的电流大于TL431的R端吸取电流的100倍,此时工程实践上基本可以忽略掉TL431的R端吸取电流的影响了。R4电流最小为2uA*10…

CSS基础:选择器和声明样式

CSS概念 CSS(Cascading Style Sheets)层叠样式表,又叫级联样式表,简称样式表 CSS用于HTML文档中元素样式的定义 使用css让网页具有美观一致的页面 语法 CSS 规则由两个主要的部分构成:选择器和声明样式 选择器通常…

DaVinci:色度 - 亮度网格应用

调色页面:色彩扭曲器 Color:Color Warper色彩扭曲器中的色度 - 亮度 Chroma - Luma网格提供了强大且直观的调色功能,相对于色相 - 饱和度网格,色度 - 亮度网格在颜色的亮度控制上更具优势。在使用网格调色之前,最好先确…

DOM编程-复选框的全选和取消全选

<!DOCTYPE html> <html> <head> <meta charset"utf-8"> <title>复选框的全选和取消全选</title> </head> <body> <script type"text/javascript"> …

LoadRunner

目录 为什么需要性能测试 性能测试实施流程 常见的性能测试指标 性能测试分类 1、一般性能测试 2、负载测试 3、压力测试 LoadRunner LoadRunner包括三个组件 VUG Controller Analysis 一个网站或者app的性能差&#xff0c;用户的使用体验就会很差 常见的性能问题&a…

LabVIEW中使用.NET方法时出现错误1316

LabVIEW中使用.NET方法时出现错误1316为什么不能调用带有泛型参数的方法&#xff1f;LabVIEW不支持哪些.NET功能&#xff1f;为什么会收到以下错误&#xff1a;发生此错误的原因是正在调用LabVIEW中不支持的.NET功能。有关解决方法&#xff0c;请参阅“其他信息”部分。可以在下…

04- Matplotlib数据可视化详解 (数据库)

Matplotlib的亮点: import matplotlib.pyplot as plt # 导包plt.figure(figsize (9, 6) , 设置图片大小plt. plot(x, y), 画图绘制网格线: 线型, 颜色, 透明度plt.grid(linestyle --, color green, alpha0.75) # linestyle: 样式, color: 颜色, alpha: 透明度plt.axis(…

【堆】数据结构堆的实现(万字详解)

前言&#xff1a; 在上一期中我们讲到了树以及二叉树的基本的概念&#xff0c;有了之前的认识&#xff0c;今天我们将来具体实现一种二叉树的存储结构“堆”&#xff01;&#xff01;&#xff01; 目录1.二叉树顺序结构介绍2.堆的概念及结构3.调整算法3.1向上调整算法3.1.1算法…

消息中间件-RocketMQ入门 消息发送的三种方式

消息中间件-RocketMQ入门 消息发送的三种方式消息中间件简介应用场景常用消息中间件RocketMQ核心概念入门案例-生产者和消费者代码逻辑消息发送的三种方式同步发送异步发送一次性消息消息中间件简介 应用场景 假设现在有订单微服务和积分微服务,正常请求流程之后是不是一个订…