参考这篇论文第四章
1. 字符含义
b
b
b: 载体坐标系
i
i
i: 惯性坐标系
D
D
D:odometer
n
n
n:导航坐标系(东北天坐标系)
m
m
m:IMU坐标系
C
b
n
C_b^n
Cbn:表示b到n旋转矩阵或者是姿态矩阵
q
b
n
q_b^n
qbn:表示b到n旋转四元素
w
n
b
b
w_{nb}^b
wnbb:表示b系相对n系的旋转在b系中的表达
2. 系统建模
状态量:位置(东北天)误差、速度误差、姿态误差(失准角)、陀螺仪零偏误差、加速计零偏误差、安装误差角误差(heading,pitch)、轮速计尺度因子误差,共18维。
[
δ
p
x
,
δ
p
y
,
δ
p
z
,
δ
v
x
,
δ
v
y
,
δ
v
z
,
ϕ
x
,
ϕ
y
,
ϕ
z
,
δ
∇
x
,
δ
∇
y
,
δ
∇
z
,
δ
ϵ
x
,
δ
ϵ
y
,
δ
ϵ
z
,
δ
α
y
,
δ
α
z
,
δ
s
]
[\delta p_x, \delta p_y, \delta p_z, \delta v_x, \delta v_y, \delta v_z, \phi_x,\phi_y,\phi_z, \delta\nabla_x, \delta\nabla_y , \delta \nabla_z, \delta\epsilon_x, \delta\epsilon_y, \delta \epsilon_z, \delta \alpha_y, \delta \alpha_z, \delta s]
[δpx,δpy,δpz,δvx,δvy,δvz,ϕx,ϕy,ϕz,δ∇x,δ∇y,δ∇z,δϵx,δϵy,δϵz,δαy,δαz,δs]
观测量1:GPS位置、速度和航向。
观测量2:[轮速,0,0]。采用非完整性约束。
特别注意1:此处失准角的定义为左扰动模型,定义为真值n系旋转到计算机推算得到n’系之间的旋转角度。必须与后面的correct相对应。因为姿态角计算不是简单的加减法。
C
n
′
n
=
I
+
[
ϕ
]
×
C_{n'}^n=I+[\phi]_\times
Cn′n=I+[ϕ]×
因此加上右扰动的姿态求导如下式,这个式子是后面求误差传播方程的基础。
δ
C
b
n
=
C
n
n
′
∗
C
b
n
−
C
b
n
=
−
[
ϕ
]
×
C
b
n
\delta C_b^n=C_n^{n'}*C_b^n-C_b^n=-[\phi]_\times C_b^n
δCbn=Cnn′∗Cbn−Cbn=−[ϕ]×Cbn
特别注意2:此处安装角的失准角使用了右扰动模型,定义为真值m系旋转到计算机推算得到的m’系之间的旋转角。必须与后面的correct相对应。因为姿态角计算不是简单的加减法。
C
m
′
m
=
I
+
[
δ
α
]
×
C_{m'}^m=I+[\delta \alpha]_\times
Cm′m=I+[δα]×
因此加上左扰动的姿态求导如下式,这个式子是后面求误差传播方程的基础。
δ
C
m
b
=
C
m
b
∗
C
m
′
m
−
C
m
b
=
C
m
b
∗
[
δ
α
]
×
\delta C_m^b=C_m^b*C_{m'}^m-C_m^b=C_m^b*[\delta \alpha]_\times
δCmb=Cmb∗Cm′m−Cmb=Cmb∗[δα]×
特别注意3
使用eigen求等效旋转向量时
- IMU测量了一个 b k b_k bk系到 b k + 1 b_{k+1} bk+1系的角度旋转量 Δ θ \Delta \theta Δθ,则将 b k + 1 b_{k+1} bk+1系的坐标转到 b k b_k bk系的坐标的等效旋转矢量可以Eigen::AngleAxisd(theta.norm(), theta.normalized()).toRotationMatrix()求出。
- 1中所述也可以和失准角联系起来。
3 状态微分方程
C ˙ b n = C b n ∗ ω n b b v ˙ n = C b n ∗ f b + g n p ˙ n = v n \dot C_b^n = C_b^n *\omega_{nb}^b \\ \dot v^n=C_b^n*f^b+g^n \\ \dot p^n = v^n C˙bn=Cbn∗ωnbbv˙n=Cbn∗fb+gnp˙n=vn
4. 误差传播方程
4.1 姿态误差传播方程
推导过程从教科书中15维姿态误差微分方程推出。
ϕ
˙
n
=
−
C
b
n
∗
δ
ω
n
b
b
ω
n
b
b
=
C
m
b
∗
(
ω
n
m
m
−
ϵ
n
m
m
)
δ
ω
n
b
b
=
δ
C
m
b
∗
(
ω
n
m
m
−
ϵ
n
m
m
)
−
C
m
b
∗
δ
ϵ
n
m
m
δ
ω
n
b
b
=
C
m
b
∗
[
δ
α
]
∗
(
ω
n
m
m
−
ϵ
n
m
m
)
−
C
m
b
∗
δ
ϵ
n
m
m
\dot \phi^n = -C_b^n *\delta \omega_{nb}^b \\ \omega_{nb}^b = C_m^b*(\omega_{nm}^m-\epsilon_{nm}^m) \\ \delta\omega_{nb}^b = \delta C_m^b* (\omega_{nm}^m-\epsilon_{nm}^m)-C_m^b* \delta\epsilon_{nm}^m\\ \delta\omega_{nb}^b = C_m^b*[\delta \alpha]*( \omega_{nm}^m-\epsilon_{nm}^m)-C_m^b* \delta\epsilon_{nm}^m
ϕ˙n=−Cbn∗δωnbbωnbb=Cmb∗(ωnmm−ϵnmm)δωnbb=δCmb∗(ωnmm−ϵnmm)−Cmb∗δϵnmmδωnbb=Cmb∗[δα]∗(ωnmm−ϵnmm)−Cmb∗δϵnmm
因此:
ϕ
˙
n
=
C
b
n
∗
C
m
b
∗
[
ω
n
m
m
−
ϵ
n
m
m
]
∗
δ
α
+
C
b
n
∗
C
m
b
∗
δ
ϵ
n
m
m
\dot \phi^n = C_b^n *C_m^b* [\omega_{nm}^m-\epsilon_{nm}^m]* \delta \alpha +C_b^n *C_m^b*\delta\epsilon_{nm}^m \\
ϕ˙n=Cbn∗Cmb∗[ωnmm−ϵnmm]∗δα+Cbn∗Cmb∗δϵnmm
4.2 速度误差传播方程
推导过程从教科书中15维姿态误差微分方程推出。
δ
v
˙
n
=
C
b
n
∗
f
b
×
ϕ
n
+
C
b
n
∗
δ
f
b
f
b
=
C
m
b
∗
(
f
m
−
∇
m
)
δ
f
b
=
δ
C
m
b
∗
(
f
m
−
∇
m
)
−
C
m
b
∗
δ
∇
m
δ
f
b
=
C
m
b
∗
[
δ
α
]
∗
(
f
m
−
∇
m
)
−
C
m
b
∗
δ
∇
m
\delta \dot v^n=C_b^n*f^b \times \phi^n + C_b^n *\delta f^b\\ f^b=C_m^b*( f^m-\nabla^m) \\ \delta f^b=\delta C_m^b* (f^m-\nabla^m)-C_m^b* \delta \nabla^m\\ \delta f^b= C_m^b*[\delta \alpha] * (f^m-\nabla^m)-C_m^b* \delta \nabla^m
δv˙n=Cbn∗fb×ϕn+Cbn∗δfbfb=Cmb∗(fm−∇m)δfb=δCmb∗(fm−∇m)−Cmb∗δ∇mδfb=Cmb∗[δα]∗(fm−∇m)−Cmb∗δ∇m
因此:
δ
v
˙
n
=
C
b
n
∗
C
m
b
∗
(
f
m
−
∇
m
)
×
ϕ
n
−
C
b
n
∗
C
m
b
∗
[
f
m
−
∇
m
]
∗
δ
α
−
C
b
n
∗
C
m
b
∗
δ
∇
n
m
m
\delta \dot v^n=C_b^n*C_m^b*(f^m-\nabla^m) \times \phi^n - C_b^n * C_m^b* [f^m-\nabla^m]*\delta \alpha - C_b^n *C_m^b* \delta\nabla_{nm}^m \\
δv˙n=Cbn∗Cmb∗(fm−∇m)×ϕn−Cbn∗Cmb∗[fm−∇m]∗δα−Cbn∗Cmb∗δ∇nmm
4.3 位置误差传播方程
推导过程从教科书中15维姿态误差微分方程推出。
δ
p
˙
n
=
δ
v
n
\delta \dot p^n = \delta v^n
δp˙n=δvn
5.误差传播矩阵
X
˙
=
F
X
+
G
W
\dot X=FX+GW
X˙=FX+GW
令
F
ϕ
=
C
b
n
∗
C
m
b
∗
[
ω
n
m
m
−
ϵ
n
m
m
]
F
ϕ
=
[
F
ϕ
(
:
,
2
:
3
)
;
0
3
×
1
]
F
v
=
−
C
b
n
∗
C
m
b
∗
[
f
m
−
∇
n
m
m
]
F
v
=
[
F
v
(
:
,
2
:
3
)
;
0
3
×
1
]
F_\phi=C_b^n * C_m^b *[ \omega_{nm}^m-\epsilon_{nm}^m]\\ F_\phi=[F_\phi(:,2:3);0_{3\times1}]\\ F_v=- C_b^n * C_m^b * [f^m-\nabla_{nm}^m]\\ F_v=[F_v(:,2:3);0_{3\times1}]\\
Fϕ=Cbn∗Cmb∗[ωnmm−ϵnmm]Fϕ=[Fϕ(:,2:3);03×1]Fv=−Cbn∗Cmb∗[fm−∇nmm]Fv=[Fv(:,2:3);03×1]
5.1 误差传播矩阵F
F = [ 0 3 × 3 I 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 [ C b n C m b f m ] C b n C m b 0 3 × 3 F v 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 − C b n C m b F ϕ 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 0 3 × 3 ] F = \begin{bmatrix} 0_{3\times3} & I_{3 \times 3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} \\ 0_{3\times3} & 0_{3\times3} &[C_b^n C_m^bf^m] & C_b^n C_m^b & 0_{3\times3} & F_v\\ 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & -C_b^n C_m^b & F_\phi\\ 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} \\ 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} \\ 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} \\ \end{bmatrix} F= 03×303×303×303×303×303×3I3×303×303×303×303×303×303×3[CbnCmbfm]03×303×303×303×303×3CbnCmb03×303×303×303×303×303×3−CbnCmb03×303×303×303×3FvFϕ03×303×303×3
5.2 噪声分配矩阵G
系统噪声输入
ω
n
o
i
s
e
=
[
w
a
w
a
l
k
_
x
w
a
w
a
l
k
_
y
w
a
w
a
l
k
_
z
w
g
w
a
l
k
_
x
w
g
w
a
l
k
_
y
w
g
w
a
l
k
_
z
w
a
o
f
f
s
e
t
_
x
w
a
o
f
f
s
e
t
_
y
w
o
f
f
s
e
t
_
z
w
g
o
f
f
s
e
t
_
x
w
g
o
f
f
s
e
t
_
y
w
g
o
f
f
s
e
t
_
z
]
\omega_{noise}= \begin{bmatrix} wa_{walk\_x} &wa_{walk\_y} &wa_{walk\_z} & wg_{walk\_x} &wg_{walk\_y} &wg_{walk\_z} & wa_{offset\_x} &wa_{offset\_y} &w_{offset\_z}& wg_{offset\_x} &wg_{offset\_y} &wg_{offset\_z} \end{bmatrix}
ωnoise=[wawalk_xwawalk_ywawalk_zwgwalk_xwgwalk_ywgwalk_zwaoffset_xwaoffset_ywoffset_zwgoffset_xwgoffset_ywgoffset_z]
噪声分配矩阵
G
=
[
0
3
×
3
0
3
×
3
0
3
×
3
0
3
×
3
0
3
×
1
C
b
n
0
3
×
3
0
3
×
3
0
3
×
3
0
3
×
1
0
3
×
3
C
b
n
0
3
×
3
0
3
×
3
0
3
×
1
0
3
×
3
0
3
×
3
E
3
×
3
0
3
×
3
0
3
×
1
0
3
×
3
0
3
×
3
0
3
×
3
E
3
×
3
0
3
×
1
0
3
×
3
0
3
×
3
0
3
×
3
0
3
×
3
0
3
×
1
]
G= \begin{bmatrix} 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times1} \\ C_b^n & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times1} \\ 0_{3\times3} & C_b^n & 0_{3\times3} & 0_{3\times3} & 0_{3\times1} \\ 0_{3\times3} & 0_{3\times3} & E_{3\times3} & 0_{3\times3} & 0_{3\times1} \\ 0_{3\times3} & 0_{3\times3} &0_{3\times3} & E_{3\times3} & 0_{3\times1} \\ 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times3} & 0_{3\times1} \end{bmatrix}
G=
03×3Cbn03×303×303×303×303×303×3Cbn03×303×303×303×303×303×3E3×303×303×303×303×303×303×3E3×303×303×103×103×103×103×103×1
6.误差传播离散化
X
k
+
1
=
Φ
k
+
1
,
k
X
k
+
Γ
k
+
1
,
k
W
X_{k+1}=\Phi_{k+1,k}X_{k}+\Gamma_{k+1,k}W
Xk+1=Φk+1,kXk+Γk+1,kW
其中
Φ
k
+
1
,
k
=
1
+
F
∗
T
Γ
k
+
1
,
k
=
G
∗
T
\Phi_{k+1,k} = 1+F*T \\ \Gamma_{k+1,k} = G*T
Φk+1,k=1+F∗TΓk+1,k=G∗T
7.状态更新方程
第一步对IMU测量去零偏
f
u
n
b
i
a
s
m
=
f
m
−
f
b
i
a
s
g
u
n
b
i
a
s
m
=
g
m
−
g
b
i
a
s
f_{unbias}^m=f^m-f_{bias}\\ g_{unbias}^m=g^m-g_{bias}
funbiasm=fm−fbiasgunbiasm=gm−gbias
第二步对校正安装误差角
f
b
=
C
m
b
∗
f
u
n
b
i
a
s
m
g
b
=
C
m
b
∗
g
u
n
b
i
a
s
m
f^b = C_m^b*f_{unbias}^m\\ g^b = C_m^b*g_{unbias}^m
fb=Cmb∗funbiasmgb=Cmb∗gunbiasm
第三步进行状态量更新
p
(
k
+
1
)
=
p
(
k
)
+
v
(
k
)
∗
Δ
T
+
1
2
∗
(
C
b
n
∗
f
b
+
g
n
)
T
2
v
(
k
+
1
)
=
v
(
k
)
+
(
C
b
n
∗
f
b
+
g
n
)
∗
T
C
b
(
k
+
1
)
n
=
C
b
(
k
)
n
∗
C
b
(
k
+
1
)
b
(
k
)
=
C
b
(
k
)
n
∗
C
(
Δ
θ
)
a
o
f
f
s
e
t
_
k
+
1
=
a
o
f
f
s
e
t
_
k
w
o
f
f
s
e
t
_
k
+
1
=
w
o
f
f
s
e
t
_
k
α
k
+
1
=
α
k
s
k
+
1
=
s
k
p(k+1) = p(k) + v(k)*\Delta T + \frac {1}{2}*(C_b^n*f^b+g^n) T^2\\ v(k+1) = v(k) + (C_b^n*f^b+g^n)* T\\ C_{b(k+1)}^n=C_{b(k)}^n*C_{b(k+1)}^{b(k)}=C_{b(k)}^n*C(\Delta \theta)\\ a_{offset\_k+1}=a_{offset\_k}\\ w_{offset\_k+1}=w_{offset\_k}\\ \alpha_{k+1}=\alpha_{k}\\ s_{k+1}=s_{k}
p(k+1)=p(k)+v(k)∗ΔT+21∗(Cbn∗fb+gn)T2v(k+1)=v(k)+(Cbn∗fb+gn)∗TCb(k+1)n=Cb(k)n∗Cb(k+1)b(k)=Cb(k)n∗C(Δθ)aoffset_k+1=aoffset_kwoffset_k+1=woffset_kαk+1=αksk+1=sk
8. 观测方程
观测量采用序贯处理的办法
8.1 GPS 的观测如下
δ p o b s = δ ( p g p s − p ( k ) ) = − E 3 ∗ 3 ∗ δ p δ v o b s = δ ( v g p s − v ( k ) ) = − E 3 ∗ 3 ∗ δ v δ y a w o b s = δ ( y a w g p s − y a w ( k ) ) = − E 1 ∗ 1 ∗ ϕ z \delta p_{obs} = \delta (p_{gps}-p(k))=-E_{3*3} * \delta p \\ \delta v_{obs} =\delta (v_{gps}-v(k))=-E_{3*3} * \delta v \\ \delta yaw_{obs} =\delta (yaw_{gps}-yaw(k)) = -E_{1*1} * \phi_z δpobs=δ(pgps−p(k))=−E3∗3∗δpδvobs=δ(vgps−v(k))=−E3∗3∗δvδyawobs=δ(yawgps−yaw(k))=−E1∗1∗ϕz
8.2 轮速的观测如下
v
D
b
=
[
s
∗
v
w
h
e
e
l
0
0
]
′
v
d
b
=
[
v
w
h
e
e
l
0
0
]
′
v
D
n
=
C
b
n
∗
v
D
b
v_D^b = [s*v_{wheel} \quad 0 \quad 0]' \\ v_d^b = [v_{wheel} \quad 0 \quad 0]' \\ v_D^n = C_b^n*v_D^b
vDb=[s∗vwheel00]′vdb=[vwheel00]′vDn=Cbn∗vDb
因此,可以求得里程计速度的误差表达式。
δ
v
D
n
=
δ
C
b
n
∗
v
D
b
+
C
b
n
∗
δ
v
D
b
δ
v
D
n
=
C
b
n
∗
v
D
b
×
ϕ
+
C
b
n
∗
δ
v
D
b
δ
v
D
b
=
[
v
w
h
e
e
l
0
0
]
′
∗
δ
s
=
v
d
∗
δ
s
\delta v_D^n =\delta C_b^n*v_D^b + C_b^n*\delta v_D^b \\ \delta v_D^n = C_b^n*v_D^b \times \phi + C_b^n*\delta v_D^b \\ \delta v_D^b = [v_{wheel} \quad 0 \quad 0]'*\delta s=v_d*\delta s
δvDn=δCbn∗vDb+Cbn∗δvDbδvDn=Cbn∗vDb×ϕ+Cbn∗δvDbδvDb=[vwheel00]′∗δs=vd∗δs
于是得到轮速的观测误差方程如下式。
δ
v
o
b
s
=
δ
(
v
D
n
−
v
(
k
)
)
=
δ
v
D
n
−
δ
v
(
k
)
=
C
b
n
∗
v
D
b
×
ϕ
+
C
b
n
∗
v
d
∗
δ
s
−
δ
v
(
k
)
\delta v_{obs}= \delta (v_D^n - v(k))= \delta v_D^n - \delta v(k)\\ =C_b^n*v_D^b \times \phi + C_b^n*v_d*\delta s-\delta v(k)
δvobs=δ(vDn−v(k))=δvDn−δv(k)=Cbn∗vDb×ϕ+Cbn∗vd∗δs−δv(k)
因此,观测矩阵可以表达为下式。
H
=
[
0
3
×
3
−
E
3
×
3
[
C
b
n
∗
v
D
b
]
0
3
×
3
0
3
×
3
0
3
×
2
C
b
n
∗
v
d
]
H = \begin{bmatrix} 0_{3\times3} & -E_{3\times3} & [C_b^n*v_D^b] & 0_{3\times3} & 0_{3\times3} & 0_{3\times2} & C_b^n*v_d \end{bmatrix}
H=[03×3−E3×3[Cbn∗vDb]03×303×303×2Cbn∗vd]
注1:按照上面观测的方程写法,求状态量更新的时候是X=X_pre-delta_X。如果想使用加法需要加个负号。
注2:由于姿态correct的时候不是简单的加减法,所以此处不要改变失准角对应的观测矩阵项的符号,按照失准角的定义进行correct。
如果想用X=X_pre+delta_X,GPS观测矩阵是单位矩阵E,轮速的观测矩阵如下。
H
=
[
0
3
×
3
E
3
×
3
[
C
b
n
∗
v
D
b
]
0
3
×
3
0
3
×
3
0
3
×
2
-
C
b
n
∗
v
d
]
H = \begin{bmatrix} 0_{3\times3} & E_{3\times3} & [C_b^n*v_D^b] & 0_{3\times3} & 0_{3\times3} & 0_{3\times2} & -C_b^n*v_d \end{bmatrix}
H=[03×3E3×3[Cbn∗vDb]03×303×303×2-Cbn∗vd]
姿态correct都是下式,必须和模型建立时的右扰动定义和失准角定义相同。
C
b
n
=
C
n
′
n
∗
C
b
n
′
=
(
I
+
[
ϕ
]
)
C
b
n
′
C
m
b
=
C
m
′
b
∗
C
m
m
′
=
C
m
′
b
∗
(
I
−
[
δ
α
]
)
C_b^n =C_{n'}^n *C_b^{n'}=(I + [\phi])C_b^{n'}\\ C_m^b =C_{m'}^b *C_m^{m'}=C_{m'}^b*(I-[\delta \alpha])
Cbn=Cn′n∗Cbn′=(I+[ϕ])Cbn′Cmb=Cm′b∗Cmm′=Cm′b∗(I−[δα])
7. 杆臂补偿
假设IMU安装在后轴中心,是车体的原点位置。
P
g
P_g
Pg:代表gnss观测在导航系的位置
v
g
v_g
vg:代表gnss观测在导航系的速度
C
b
n
C_b^n
Cbn: 代表当前车体在导航系的姿态
L
L
L:代表杆臂向量,gnss天线在车体系下的位置
P
c
a
r
P_{car}
Pcar:待求的GNSS观测的车体在导航系得位置。
因此位置杆臂补偿如下式。
P
g
=
C
b
n
∗
L
+
P
c
a
r
P
c
a
r
=
P
g
−
C
b
n
∗
L
P_g=C_b^n*L+P_{car}\\ P_{car} = P_g - C_b^n*L
Pg=Cbn∗L+PcarPcar=Pg−Cbn∗L
车体为刚体,天线和车体的转动角速度是一致的,因此速度的杆臂补偿公式如下
V
c
a
r
=
V
g
−
C
b
n
∗
(
ω
×
L
)
V_{car} = V_g-C_b^n*(\omega \times L)
Vcar=Vg−Cbn∗(ω×L)
ω
\omega
ω为IM U测量的角速度。
8. 完整性约束
轮速观测中已经带入了完整性约束。
9. 零速修正&零航向修正
10. 时间同步
11. 误差标定处理
系统建模时已经带入了安装误差和轮速误差建模。