本小节介绍了传统的旋转平移矩阵、四元数、对偶四元数以及鼎立四元数的定义以及运动学表示。
旋转平移矩阵
定义与性质
旋转矩阵 R \mathbf{R} R 是一个用于表示三维空间中绕原点的旋转的正交矩阵。它有以下性质:
-
行列式: det ( R ) = 1 \det(\mathbf{R}) = 1 det(R)=1,意味着它是特殊正交矩阵,属于旋转群 S O ( 3 ) SO(3) SO(3)。
-
逆矩阵:旋转矩阵的逆等于它的转置,即 R − 1 = R ⊤ \mathbf{R}^{-1} = \mathbf{R}^\top R−1=R⊤。
-
正交性:所有行向量和列向量都是单位向量,并且相互正交。
平移向量 t \mathbf{t} t 是一个三维向量,用于表示空间中的直线平移。它可以定义为:
t = [ t x t y t z ] \mathbf{t} = \begin{bmatrix} t_x \\ t_y \\ t_z \end{bmatrix} t= txtytz
其中 t x , t y , t z t_x, t_y, t_z tx,ty,tz 分别是沿着 x x x, y y y, 和 z z z 轴的平移量。
刚体变换包括一个旋转和一个平移,可以用齐次坐标表示为 4 × 4 4 \times 4 4×4 矩阵:
T = [ R t 0 1 ] \mathbf{T} = \begin{bmatrix} \mathbf{R} & \mathbf{t} \\ \mathbf{0} & 1 \end{bmatrix} T=[R0t1]
其中 R \mathbf{R} R 是旋转矩阵, t \mathbf{t} t 是平移向量, 0 \mathbf{0} 0 是 1 × 3 1 \times 3 1×3 的零向量。这个变换矩阵可以用来将一个点 p \mathbf{p} p 从一个坐标系变换到另一个坐标系:
p ′ = T p \mathbf{p}' = \mathbf{T} \mathbf{p} p′=Tp
这里 p \mathbf{p} p 和 p ′ \mathbf{p}' p′ 都是齐次坐标下的点。
运动学表示
在一个SINS系统中,传统方法对姿态的描述包括旋转、速度、位置,并通过IMU所采集的数据为载体系的角速度与加速度。对于传统方法,我们的系统状态量为 χ i = [ R b i e , v i e , t i e ] , \mathbf{\chi}_i = [\mathbf{R}_{b_i}^e,\mathbf{v}_i^e,\mathbf{t}_i^e], χi=[Rbie,vie,tie], 其中 χ i {\mathbf{\chi}_i} χi为 i i i时刻姿态表示的状态, R b i e \mathbf{R}_{b_i}^e Rbie为载体相对于参考系的旋转矩阵, v i e \mathbf{v}_i^e vie为载体相对参考系的速度, t i e \mathbf{t}_i^e tie为相对参考系的位移。设在 i i i时刻时IMU所测得的角速度为 ω i b \mathbf{\omega}_i^b ωib,加速度为 a i b \mathbf{a}_i^b aib,则从时刻 i i i到 i + 1 i+1 i+1的状态推算需要计算状态随时间变化的导数,其中, v i e \mathbf{v}_i^e vie与 t i e \mathbf{t}_i^e tie的计算易得: v ˙ i e = R b i e a i b − g , \dot{\mathbf{v}}_i^e = \mathbf{R}_{b_i}^e\mathbf{a}_i^b - \mathbf{g}, v˙ie=Rbieaib−g, t ˙ i e = v i e , \dot{\mathbf{t}}_i^e = {\mathbf{v}}_i^e, t˙ie=vie, 其中 g \mathbf{g} g为重力加速度。在SLAM系统中,地球的自转通常可以被忽略,因此我对重力模型相较于论文中进行了简化。旋转矩阵的求导相对复杂,首先这里先介绍罗德里格斯公式(Rodrigues’ rotation formula)。
罗德里格斯公式提供了一个绕单位向量 u = ( u x , u y , u z ) \mathbf{u} = (u_x, u_y, u_z) u=(ux,uy,uz) 旋转角度 θ \theta θ 的旋转矩阵 R \mathbf{R} R。公式如下:
R = ( cos θ ) I + ( sin θ ) K + ( 1 − cos θ ) K 2 , \mathbf{R} = (\cos \theta) \mathbf{I} + (\sin \theta) \mathbf{K} + (1 - \cos \theta) \mathbf{K}^2, R=(cosθ)I+(sinθ)K+(1−cosθ)K2, 其中 I \mathbf{I} I 是单位矩阵, K \mathbf{K} K 是关于 u \mathbf{u} u 的反对称矩阵,定义为: K = ( 0 − u z u y u z 0 − u x − u y u x 0 ) = u × . \mathbf{K} = \begin{pmatrix} 0 & -u_z & u_y \\ u_z & 0 & -u_x \\ -u_y & u_x & 0 \end{pmatrix} = \mathbf{u} \times. K= 0uz−uy−uz0uxuy−ux0 =u×. 对于旋转矩阵求导,根据导数定义: R ˙ b i e = lim Δ t → 0 R b i e ( t + Δ t ) − R b i e ( t ) Δ t . \dot{\mathbf{R}}_{b_i}^e = \lim_{\Delta t \to 0} \frac{\mathbf{R}_{b_i}^e(t + \Delta t) - \mathbf{R}_{b_i}^e(t)}{\Delta t}. R˙bie=Δt→0limΔtRbie(t+Δt)−Rbie(t). 在极短的时间内,假设从 i i i到 i + 1 i+1 i+1时刻,可以得到 R b i e ( t + Δ t ) = R b i + 1 e = R b i e ( t ) R b i + 1 b i , \mathbf{R}_{b_i}^e(t + \Delta t) = \mathbf{R}_{b_{i+1}}^e = \mathbf{R}_{b_i}^e(t) \mathbf{R}_{b_{i+1}}^{b_i}, Rbie(t+Δt)=Rbi+1e=Rbie(t)Rbi+1bi, 其中, R b i + 1 b i \mathbf{R}_{b_{i+1}}^{b_i} Rbi+1bi为载体在 i + 1 i+1 i+1时刻相对于 i i i时刻的相对旋转矩阵。当旋转在极短时间内为小量时,利用罗德里格斯公式,对 R b i + 1 b i \mathbf{R}_{b_{i+1}}^{b_i} Rbi+1bi进行近似,得到: R b i + 1 b i ≈ I + θ b i + 1 b i K b i + 1 b i . \mathbf{R}_{b_{i+1}}^{b_i} \approx \mathbf{I} + \theta_{b_{i+1}}^{b_i} \mathbf{K}_{b_{i+1}}^{b_i}. Rbi+1bi≈I+θbi+1biKbi+1bi. 将(8)代入(6)中得到: R ˙ b i e = lim Δ t → 0 R b i e ( t ) ( I + θ b i + 1 b i K b i + 1 b i ) − R b i e ( t ) Δ t = lim Δ t → 0 R b i e ( t ) θ b i + 1 b i K b i + 1 b i Δ t , \dot{\mathbf{R}}_{b_i}^e = \lim_{\Delta t \to 0} \frac{\mathbf{R}_{b_i}^e(t)(\mathbf{I} + \theta_{b_{i+1}}^{b_i} \mathbf{K}_{b_{i+1}}^{b_i}) - \mathbf{R}_{b_i}^e(t)}{\Delta t} = \lim_{\Delta t \to 0} \frac{\mathbf{R}_{b_i}^e(t)\theta_{b_{i+1}}^{b_i} \mathbf{K}_{b_{i+1}}^{b_i}} {\Delta t}, R˙bie=Δt→0limΔtRbie(t)(I+θbi+1biKbi+1bi)−Rbie(t)=Δt→0limΔtRbie(t)θbi+1biKbi+1bi, 显然, lim Δ t → 0 θ b i + 1 b i K b i + 1 b i Δ t = ω i b × . \lim_{\Delta t \to 0} \frac{\theta_{b_{i+1}}^{b_i} \mathbf{K}_{b_{i+1}}^{b_i}} {\Delta t} = \mathbf{\omega}_i^b \times. Δt→0limΔtθbi+1biKbi+1bi=ωib×. 因此,旋转矩阵对时间的导数为: R ˙ b i e = R b i e ω i b × . \dot{\mathbf{R}}_{b_i}^e = \mathbf{R}_{b_i}^e \mathbf{\omega}_i^b \times. R˙bie=Rbieωib×.
综上,传统的姿态状态表示方式需要 3 × 3 + 3 × 1 + 3 × 1 = 15 3 \times 3 + 3 \times 1 + 3 \times 1 = 15 3×3+3×1+3×1=15维来描述。在旋转、平移、速度这三个量上的表示相对独立。
状态空间 | 运动学表示 |
---|---|
χ i = [ R b i e , v i e , t i e ] \mathbf{\chi}_i = [\mathbf{R}_{b_i}^e,\mathbf{v}_i^e,\mathbf{t}_i^e] χi=[Rbie,vie,tie] | R ˙ b i e = R b i e ω i b × \dot{\mathbf{R}}_{b_i}^e = \mathbf{R}_{b_i}^e \mathbf{\omega}_i^b \times R˙bie=Rbieωib× ; v ˙ i e = R b i e a i b − g \dot{\mathbf{v}}_i^e = \mathbf{R}_{b_i}^e\mathbf{a}_i^b - \mathbf{g} v˙ie=Rbieaib−g ; t ˙ i e = v i e \dot{\mathbf{t}}_i^e = {\mathbf{v}}_i^e t˙ie=vie |
传统运动学表示
传统四元数
定义与性质
四元数是复数的一个扩展,用于在三维空间中表示旋转。它由一个实部和三个虚部组成,可以表示为 q = a + b i + c j + d k q = a + bi + cj + dk q=a+bi+cj+dk,其中 a , b , c , d ∈ R a, b, c, d \in \mathbb{R} a,b,c,d∈R,并且 i , j , k i, j, k i,j,k 是四元数的虚单位。
四元数的基本性质包括:
-
四元数的加法和实数的加法类似,遵循交换律和结合律。
-
四元数的乘法不遵循交换律,但遵循结合律。
-
四元数的虚单位满足以下关系: i 2 = j 2 = k 2 = i j k = − 1 , i j = k , j i = − k , j k = i , k j = − i , k i = j , i k = − j . \begin{aligned} i^2 &= j^2 = k^2 = ijk = -1, \\ ij &= k, \quad ji = -k, \\ jk &= i, \quad kj = -i, \\ ki &= j, \quad ik = -j. \end{aligned} i2ijjkki=j2=k2=ijk=−1,=k,ji=−k,=i,kj=−i,=j,ik=−j.
-
四元数的共轭 q ∗ q^* q∗ 定义为 q ∗ = a − b i − c j − d k q^* = a - bi - cj - dk q∗=a−bi−cj−dk。
-
四元数的模(或范数) ∥ q ∥ \|q\| ∥q∥ 定义为 ∥ q ∥ = a 2 + b 2 + c 2 + d 2 \|q\| = \sqrt{a^2 + b^2 + c^2 + d^2} ∥q∥=a2+b2+c2+d2。
-
四元数的逆 q − 1 q^{-1} q−1(如果 q q q 不是零四元数)定义为 q − 1 = q ∗ ∥ q ∥ 2 q^{-1} = \frac{q^*}{\|q\|^2} q−1=∥q∥2q∗。
-
一个单位四元数满足 ∥ q ∥ = 1 \|q\| = 1 ∥q∥=1。
两个四元数 q 1 = a 1 + b 1 i + c 1 j + d 1 k q_1 = a_1 + b_1i + c_1j + d_1k q1=a1+b1i+c1j+d1k 和 q 2 = a 2 + b 2 i + c 2 j + d 2 k q_2 = a_2 + b_2i + c_2j + d_2k q2=a2+b2i+c2j+d2k 的加法定义为:
q 1 + q 2 = ( a 1 + a 2 ) + ( b 1 + b 2 ) i + ( c 1 + c 2 ) j + ( d 1 + d 2 ) k q_1 + q_2 = (a_1 + a_2) + (b_1 + b_2)i + (c_1 + c_2)j + (d_1 + d_2)k q1+q2=(a1+a2)+(b1+b2)i+(c1+c2)j+(d1+d2)k
四元数的乘法则较为复杂,其定义为:
q
1
⊗
q
2
=
(
a
1
a
2
−
b
1
b
2
−
c
1
c
2
−
d
1
d
2
)
+
(
a
1
b
2
+
b
1
a
2
+
c
1
d
2
−
d
1
c
2
)
i
+
(
a
1
c
2
−
b
1
d
2
+
c
1
a
2
+
d
1
b
2
)
j
+
(
a
1
d
2
+
b
1
c
2
−
c
1
b
2
+
d
1
a
2
)
k
\begin{aligned} q_1 \otimes q_2 = \, & (a_1a_2 - b_1b_2 - c_1c_2 - d_1d_2) \\ &+ (a_1b_2 + b_1a_2 + c_1d_2 - d_1c_2)i \\ &+ (a_1c_2 - b_1d_2 + c_1a_2 + d_1b_2)j \\ &+ (a_1d_2 + b_1c_2 - c_1b_2 + d_1a_2)k \end{aligned}
q1⊗q2=(a1a2−b1b2−c1c2−d1d2)+(a1b2+b1a2+c1d2−d1c2)i+(a1c2−b1d2+c1a2+d1b2)j+(a1d2+b1c2−c1b2+d1a2)k
四元数在三维空间中表示旋转时在计算上比旋转矩阵更高效,减少了旋转矩阵的维数冗余。
运动学表示
传统四元数在SLAM中一般与平移向量和速度向量联合使用。平移向量和速度向量的时间导数在上一小节已经讨论过,下面讨论传统四元数对时间的导数。
设 q b i e q_{b_i}^e qbie为载体相对于参考系的旋转四元数。 对于旋转四元数求导,根据导数定义: q ˙ b i e = lim Δ t → 0 q b i e ( t + Δ t ) − q b i e ( t ) Δ t . \dot{q}_{b_i}^e = \lim_{\Delta t \to 0} \frac{q_{b_i}^e(t + \Delta t) - q_{b_i}^e(t)}{\Delta t}. q˙bie=Δt→0limΔtqbie(t+Δt)−qbie(t). 在极短的时间内,假设从 i i i到 i + 1 i+1 i+1时刻,可以得到 q b i e ( t + Δ t ) = q b i + 1 e = q b i e ( t ) q b i + 1 b i , q_{b_i}^e(t + \Delta t) = q_{b_{i+1}}^e = q_{b_i}^e(t) q_{b_{i+1}}^{b_i}, qbie(t+Δt)=qbi+1e=qbie(t)qbi+1bi, 其中, q b i + 1 b i q_{b_{i+1}}^{b_i} qbi+1bi为载体在 i + 1 i+1 i+1时刻相对于 i i i时刻的相对旋转四元数。与旋转矩阵相似,给定一个单位旋转轴 u = ( u x , u y , u z ) \mathbf{u} = (u_x, u_y, u_z) u=(ux,uy,uz) 和一个旋转角 θ \theta θ,相应的单位四元数 q q q 可以表示为:
q = cos ( θ 2 ) + ( u x i + u y j + u z k ) sin ( θ 2 ) , q = \cos\left(\frac{\theta}{2}\right) + (u_xi + u_yj + u_zk) \sin\left(\frac{\theta}{2}\right), q=cos(2θ)+(uxi+uyj+uzk)sin(2θ),
其中 θ \theta θ 是旋转角, u \mathbf{u} u 是旋转轴,且 u x , u y , u z u_x, u_y, u_z ux,uy,uz 是 u \mathbf{u} u 的坐标。这里 i , j , k i, j, k i,j,k 是四元数的虚单位。该四元数 q q q 表示绕轴 u \mathbf{u} u 旋转 θ \theta θ 弧度的旋转。 当在极短时间内旋转一个小角度时, q b i + 1 b i q_{b_{i+1}}^{b_i} qbi+1bi可近似表示为 q b i + 1 b i ≈ 1 + ( u x i + u y j + u z k ) θ b i + 1 b i 2 q_{b_{i+1}}^{b_i} \approx 1 + (u_xi + u_yj + u_zk)\frac{\theta_{b_{i+1}}^{b_i}}{2} qbi+1bi≈1+(uxi+uyj+uzk)2θbi+1bi 将(14)代入(12)中得到: q ˙ b i e = lim Δ t → 0 q b i e ( t ) ( 1 + ( u x i + u y j + u z k ) θ b i + 1 b i 2 ) − q b i e ( t ) Δ t = lim Δ t → 0 q b i e ( t ) ⊗ [ 0 ; u θ b i + 1 b i 2 ] Δ t , \dot{q}_{b_i}^e = \lim_{\Delta t \to 0} \frac{q_{b_i}^e(t)(1 + (u_xi + u_yj + u_zk)\frac{\theta_{b_{i+1}}^{b_i}}{2}) - q_{b_i}^e(t)}{\Delta t} = \lim_{\Delta t \to 0} \frac{q_{b_i}^e(t) \otimes [0;\mathbf{u}\frac{\theta_{b_{i+1}}^{b_i}}{2}] } {\Delta t}, q˙bie=Δt→0limΔtqbie(t)(1+(uxi+uyj+uzk)2θbi+1bi)−qbie(t)=Δt→0limΔtqbie(t)⊗[0;u2θbi+1bi], 显然, lim Δ t → 0 [ 0 ; u θ b i + 1 b i 2 ] Δ t = [ 0 ; ω i b 2 ] . \lim_{\Delta t \to 0} \frac{[0;\mathbf{u}\frac{\theta_{b_{i+1}}^{b_i}}{2}]} {\Delta t} = [0;\frac{\mathbf{\omega}_i^b}{2}]. Δt→0limΔt[0;u2θbi+1bi]=[0;2ωib]. 因此, q ˙ b i e = q b i e ⊗ [ 0 ; ω i b 2 ] . \dot{q}_{b_i}^e = q_{b_i}^e \otimes [0;\frac{\mathbf{\omega}_i^b}{2}]. q˙bie=qbie⊗[0;2ωib].
综上,传统四元数可以高效表示旋转,用4维代替了旋转矩阵的9维对旋转进行了表示。
对偶四元数
定义与性质
对偶四元数是一种数学工具,用于在三维空间中同时表示旋转和平移。一个对偶四元数可以写为 q ^ = q + ϵ q ′ \hat{q} = q + \epsilon q' q^=q+ϵq′,其中 q q q 是标准四元数部分,用于表示旋转,而 q ′ q' q′ 是对偶部分,用于表示平移,且 ϵ \epsilon ϵ 是对偶单位,满足 ϵ 2 = 0 \epsilon^2 = 0 ϵ2=0 但 ϵ ≠ 0 \epsilon \neq 0 ϵ=0。 一个对偶四元数 q ^ \hat{q} q^ 由一个实四元数 q q q 和一个对偶四元数 ϵ q ′ \epsilon q' ϵq′ 组成:
q ^ = q + ϵ q ′ \hat{q} = q + \epsilon q' q^=q+ϵq′
其中 q = a + b i + c j + d k q = a + bi + cj + dk q=a+bi+cj+dk 是旋转部分, q ′ = a ′ + b ′ i + c ′ j + d ′ k q' = a' + b'i + c'j + d'k q′=a′+b′i+c′j+d′k 是平移部分, a , b , c , d , a ′ , b ′ , c ′ , d ′ a, b, c, d, a', b', c', d' a,b,c,d,a′,b′,c′,d′ 是实数。 对偶四元数的基本性质包括:
-
加法:两个对偶四元数的加法遵循逐元素加法规则。
-
乘法:两个对偶四元数的乘法不遵循交换律,但遵循结合律和分配律。
-
共轭:对偶四元数的共轭定义为 q ^ ∗ = q ∗ + ϵ q ′ ∗ \hat{q}^* = q^* + \epsilon q'^* q^∗=q∗+ϵq′∗。
-
模:对偶四元数的模定义为 ∣ q ^ ∣ = ∣ q ∣ + ϵ q ⋅ q ′ ∣ q ∣ |\hat{q}| = |q| + \epsilon \frac{q \cdot q'}{|q|} ∣q^∣=∣q∣+ϵ∣q∣q⋅q′。
对偶四元数能够紧凑地表示三维空间中的刚体变换(旋转和平移)。
运动学表示
在一个SINS系统中,对姿态的描述包括旋转、速度、位置,并通过IMU所采集的数据为载体系的角速度与加速度。对于对偶四元数方法,我们的系统状态量为 χ i = [ q ^ b i e , v i e ] , \mathbf{\chi}_i = [\hat{q}_{b_i}^e,\mathbf{v}_i^e], χi=[q^bie,vie], 其中 χ i {\mathbf{\chi}_i} χi为 i i i时刻姿态表示的状态, q ^ b i e \hat{q}_{b_i}^e q^bie为载体相对于参考系的对偶四元数,其中包含了旋转四元数 q b i e {q}_{b_i}^e qbie,及其对偶部分,其对偶部分可以抽象理解为平移的表示,可以表示为 1 2 t i e ⊗ q b i e \frac{1}{2} \mathbf{t}_i^e\otimes {q}_{b_i}^e 21tie⊗qbie, t i e \mathbf{t}_i^e tie为相对参考系的位移。对偶四元数形式如下: q ^ b i e = q b i e + ϵ 2 t i e ⊗ q b i e , \hat{q}_{b_i}^e = {q}_{b_i}^e + \frac{\epsilon}{2} \mathbf{t}_i^e \otimes {q}_{b_i}^e, q^bie=qbie+2ϵtie⊗qbie, v i e \mathbf{v}_i^e vie为载体相对参考系的速度。设在 i i i时刻时IMU所测得的角速度为 ω i b \mathbf{\omega}_i^b ωib,加速度为 a i b \mathbf{a}_i^b aib,则从时刻 i i i到 i + 1 i+1 i+1的状态推算需要计算状态随时间变化的导数,其中, v i e \mathbf{v}_i^e vie与 t i e \mathbf{t}_i^e tie的计算易得: v ˙ i e = q b i e ⊗ a i b ⊗ q b i e ∗ − g , \dot{\mathbf{v}}_i^e = {q}_{b_i}^e \otimes \mathbf{a}_i^b \otimes {{q}_{b_i}^e}^* - \mathbf{g}, v˙ie=qbie⊗aib⊗qbie∗−g, t ˙ i e = v i e , \dot{\mathbf{t}}_i^e = {\mathbf{v}}_i^e, t˙ie=vie,
对于对偶四元数的时间导数可以表示为:
q ^ ˙ b i e = q ˙ b i e + ϵ 2 ( t i e ⊗ q b i e ) ′ , \dot{\hat{q}}_{b_i}^e = \dot{q}_{b_i}^e + \frac{\epsilon}{2} (\mathbf{t}_i^e \otimes {q}_{b_i}^e)', q^˙bie=q˙bie+2ϵ(tie⊗qbie)′, 利用复合函数求导规则可以得到 q ^ ˙ b i e = q ˙ b i e + ϵ 2 ( t ˙ i e ⊗ q b i e + t i e ⊗ q ˙ b i e ) . \dot{\hat{q}}_{b_i}^e = \dot{q}_{b_i}^e + \frac{\epsilon}{2} (\dot{\mathbf{t}}_i^e \otimes {q}_{b_i}^e + \mathbf{t}_i^e \otimes \dot{q}_{b_i}^e). q^˙bie=q˙bie+2ϵ(t˙ie⊗qbie+tie⊗q˙bie). 将式(18)与式(22)代入(24)可以得到 q ^ ˙ b i e = q b i e ⊗ [ 0 ; ω i b 2 ] + ϵ 2 ( [ 0 ; v i e ] ⊗ q b i e + t i e ⊗ q b i e ⊗ [ 0 ; ω i b 2 ] ) . \dot{\hat{q}}_{b_i}^e = q_{b_i}^e \otimes [0;\frac{\mathbf{\omega}_i^b}{2}] + \frac{\epsilon}{2} ({[0;\mathbf{v}}_i^e] \otimes {q}_{b_i}^e + \mathbf{t}_i^e \otimes q_{b_i}^e \otimes [0;\frac{\mathbf{\omega}_i^b}{2}]). q^˙bie=qbie⊗[0;2ωib]+2ϵ([0;vie]⊗qbie+tie⊗qbie⊗[0;2ωib]).
整理并化简上式可得 q ^ ˙ b i e = q ^ b i e ⊗ 1 2 ( ω i b + ϵ v i e ) , \dot{\hat{q}}_{b_i}^e = \hat{q}_{b_i}^e \otimes \frac{1}{2}(\boldsymbol{\omega}_i^b + \epsilon \mathbf{v}_i^e), q^˙bie=q^bie⊗21(ωib+ϵvie),
这里 ω i b \boldsymbol{\omega}_i^b ωib 和 v i e \mathbf{v}_i^e vie 分别被表示为四元数形式,即实部为零。若记 ω ^ i b = ( ω i b + ϵ v i e ) \hat{\omega}_i^b=(\boldsymbol{\omega}_i^b + \epsilon \mathbf{v}_i^e) ω^ib=(ωib+ϵvie),那么对偶四元数对时间的求导形式可与传统四元数相同: q ^ ˙ b i e = q ^ b i e ⊗ 1 2 ω ^ i b . \dot{\hat{q}}_{b_i}^e = \hat{q}_{b_i}^e \otimes \frac{1}{2}\hat{\omega}_i^b. q^˙bie=q^bie⊗21ω^ib.
综上,对偶四元数将旋转和平移统一表示在一个对偶四元数中,但是对于速度状态仍然单独考虑。
状态空间 | 运动学表示 |
---|---|
χ i = [ q ^ b i e , v i e ] \mathbf{\chi}_i = [\hat{q}_{b_i}^e,\mathbf{v}_i^e] χi=[q^bie,vie] | q ^ ˙ b i e = q ^ b i e ⊗ 1 2 ω ^ i b \dot{\hat{q}}_{b_i}^e = \hat{q}_{b_i}^e \otimes \frac{1}{2}\hat{\omega}_i^b q^˙bie=q^bie⊗21ω^ib ; v ˙ i e = q b i e ⊗ a i b − g \dot{\mathbf{v}}_i^e = {q}_{b_i}^e \otimes \mathbf{a}_i^b - \mathbf{g} v˙ie=qbie⊗aib−g |
对偶四元数运动学表示
鼎立四元数
定义与性质
鼎立四元数(TridentQs)在对偶四元数的基础上扩展了一维虚部,定义为 q ˘ = q + ϵ 1 q ′ + ϵ 2 q ′ ′ \breve{q} = q + \epsilon_1 q' + \epsilon_2 q'' q˘=q+ϵ1q′+ϵ2q′′,其中 q , q ′ , q ′ ′ q, q', q'' q,q′,q′′ 是常规的四元数,对偶单位 ϵ 1 , ϵ 2 \epsilon_1, \epsilon_2 ϵ1,ϵ2 满足 ϵ 1 2 = ϵ 2 2 = ϵ 1 ϵ 2 = 0 \epsilon_1^2 = \epsilon_2^2 = \epsilon_1 \epsilon_2 = 0 ϵ12=ϵ22=ϵ1ϵ2=0,且 ϵ 1 ≠ 0 \epsilon_1 \neq 0 ϵ1=0 和 ϵ 2 ≠ 0 \epsilon_2 \neq 0 ϵ2=0。
对于任意两个鼎立四元数 q ˘ 1 = q 1 + ϵ 1 q 1 ′ + ϵ 2 q 1 ′ ′ \breve{q}_1 = q_1 + \epsilon_1 q'_1 + \epsilon_2 q''_1 q˘1=q1+ϵ1q1′+ϵ2q1′′ 和 q ˘ 2 = q 2 + ϵ 1 q 2 ′ + ϵ 2 q 2 ′ ′ \breve{q}_2 = q_2 + \epsilon_1 q'_2 + \epsilon_2 q''_2 q˘2=q2+ϵ1q2′+ϵ2q2′′,算术运算如下:
-
加法: q ˘ 1 + q ˘ 2 = q 1 + q 2 + ϵ 1 ( q 1 ′ + q 2 ′ ) + ϵ 2 ( q 1 ′ ′ + q 2 ′ ′ ) \breve{q}_1 + \breve{q}_2 = q_1 + q_2 + \epsilon_1 (q'_1 + q'_2) + \epsilon_2 (q''_1 + q''_2) q˘1+q˘2=q1+q2+ϵ1(q1′+q2′)+ϵ2(q1′′+q2′′)
-
乘法: q ˘ 1 ⊗ q ˘ 2 = q 1 ⊗ q 2 + ϵ 1 ( q 1 ⊗ q 2 ′ + q 1 ′ ⊗ q 2 ) + ϵ 2 ( q 1 ⊗ q 2 ′ ′ + q 1 ′ ′ ⊗ q 2 ) \breve{q}_1 \otimes \breve{q}_2 = q_1 \otimes q_2 + \epsilon_1 (q_1 \otimes q'_2 + q'_1 \otimes q_2) + \epsilon_2 (q_1 \otimes q''_2 + q''_1 \otimes q_2) q˘1⊗q˘2=q1⊗q2+ϵ1(q1⊗q2′+q1′⊗q2)+ϵ2(q1⊗q2′′+q1′′⊗q2)
运动学表示
在一个SINS系统中,对姿态的描述包括旋转、速度、位置,并通过IMU所采集的数据为载体系的角速度与加速度。对于对偶四元数方法,我们的系统状态量为 χ i = [ q ˘ b i e ] , \mathbf{\chi}_i = [\breve{q}_{b_i}^e], χi=[q˘bie], 其中 χ i {\mathbf{\chi}_i} χi为 i i i时刻姿态表示的状态, q ^ b i e \hat{q}_{b_i}^e q^bie为载体相对于参考系的对偶四元数,其中包含了旋转四元数 q b i e {q}_{b_i}^e qbie,及其两个对偶部分,第一个对偶部分抽象表示了载体速度 1 2 v i e ⊗ q b i e \frac{1}{2} \mathbf{v}_i^e\otimes {q}_{b_i}^e 21vie⊗qbie, t i e \mathbf{t}_i^e tie为相对参考系的载体速度,第二个对偶部分抽象表示了载体平移,可以表示为 1 2 t i e ⊗ q b i e \frac{1}{2} \mathbf{t}_i^e\otimes {q}_{b_i}^e 21tie⊗qbie, t i e \mathbf{t}_i^e tie为相对参考系的载体位移。鼎立四元数形式如下: q ˘ b i e = q b i e + ϵ 1 2 v i e ⊗ q b i e + ϵ 2 2 t i e ⊗ q b i e , \breve{q}_{b_i}^e = {q}_{b_i}^e + \frac{\epsilon_1}{2} \mathbf{v}_i^e \otimes {q}_{b_i}^e + \frac{\epsilon_2}{2} \mathbf{t}_i^e \otimes {q}_{b_i}^e, q˘bie=qbie+2ϵ1vie⊗qbie+2ϵ2tie⊗qbie,
对于对偶四元数的时间导数可以表示为:
q ˘ ˙ b i e = q ˙ b i e + ϵ 1 2 ( v i e ⊗ q b i e ) ′ + ϵ 2 2 ( t i e ⊗ q b i e ) ′ , \dot{\breve{q}}_{b_i}^e = \dot{q}_{b_i}^e + \frac{\epsilon_1}{2} (\mathbf{v}_i^e \otimes {q}_{b_i}^e)' + \frac{\epsilon_2}{2} (\mathbf{t}_i^e \otimes {q}_{b_i}^e)', q˘˙bie=q˙bie+2ϵ1(vie⊗qbie)′+2ϵ2(tie⊗qbie)′, 利用复合函数求导规则可以得到 q ˘ ˙ b i e = q ˙ b i e + ϵ 1 2 ( v ˙ i e ⊗ q b i e + v i e ⊗ q ˙ b i e ) + ϵ 2 2 ( t ˙ i e ⊗ q b i e + t i e ⊗ q ˙ b i e ) . \dot{\breve{q}}_{b_i}^e = \dot{q}_{b_i}^e + \frac{\epsilon_1}{2} (\dot{\mathbf{v}}_i^e \otimes {q}_{b_i}^e + \mathbf{v}_i^e \otimes \dot{q}_{b_i}^e) + \frac{\epsilon_2}{2} (\dot{\mathbf{t}}_i^e \otimes {q}_{b_i}^e + \mathbf{t}_i^e \otimes \dot{q}_{b_i}^e). q˘˙bie=q˙bie+2ϵ1(v˙ie⊗qbie+vie⊗q˙bie)+2ϵ2(t˙ie⊗qbie+tie⊗q˙bie). 设在 i i i时刻时IMU所测得的角速度为 ω i b \mathbf{\omega}_i^b ωib。在鼎立四元数的推导过程中,为了形式更加优美,设加速度为载体的纯加速度,不考虑重力的影响,可表示为 a i b \mathbf{a}_i^b aib,则从时刻 i i i到 i + 1 i+1 i+1的状态推算需要计算状态随时间变化的导数,其中, v i e \mathbf{v}_i^e vie与 t i e \mathbf{t}_i^e tie的计算易得: v ˙ i e = q b i e ⊗ a i b ⊗ q b i e ∗ , \dot{\mathbf{v}}_i^e = {q}_{b_i}^e \otimes \mathbf{a}_i^b \otimes {{q}_{b_i}^e}^*, v˙ie=qbie⊗aib⊗qbie∗, t ˙ i e = v i e , \dot{\mathbf{t}}_i^e = {\mathbf{v}}_i^e, t˙ie=vie, 将式(32)及(33)代入(31)得到 q ˘ ˙ b i e = q b i e ⊗ [ 0 ; ω i b 2 ] + ϵ 1 2 ( q b i e ⊗ a i b ⊗ q b i e ∗ ⊗ q b i e + v i e ⊗ q b i e ( t ) ⊗ [ 0 ; ω i b 2 ] ) + ϵ 2 2 ( [ 0 ; v i e ] ⊗ q b i e + t i e ⊗ q b i e ( t ) ⊗ [ 0 ; ω i b 2 ] ) . \begin{aligned} \dot{\breve{q}}_{b_i}^e = q_{b_i}^e \otimes [0;\frac{\mathbf{\omega}_i^b}{2}] + \frac{\epsilon_1}{2} ({q}_{b_i}^e \otimes \mathbf{a}_i^b \otimes {{q}_{b_i}^e}^* \otimes {q}_{b_i}^e + \mathbf{v}_i^e \otimes q_{b_i}^e(t) \otimes [0;\frac{\mathbf{\omega}_i^b}{2}]) + \\ \frac{\epsilon_2}{2} ({[0;\mathbf{v}}_i^e] \otimes {q}_{b_i}^e + \mathbf{t}_i^e \otimes q_{b_i}^e(t) \otimes [0;\frac{\mathbf{\omega}_i^b}{2}]). \end{aligned} q˘˙bie=qbie⊗[0;2ωib]+2ϵ1(qbie⊗aib⊗qbie∗⊗qbie+vie⊗qbie(t)⊗[0;2ωib])+2ϵ2([0;vie]⊗qbie+tie⊗qbie(t)⊗[0;2ωib]).
通过四元数结合律可得 q ˘ ˙ b i e = q b i e ⊗ [ 0 ; ω i b 2 ] + ϵ 1 2 ( q b i e ⊗ a i b + v i e ⊗ q b i e ( t ) ⊗ [ 0 ; ω i b 2 ] ) + ϵ 2 2 ( [ 0 ; v i e ] ⊗ q b i e + t i e ⊗ q b i e ( t ) ⊗ [ 0 ; ω i b 2 ] ) . \begin{aligned} \dot{\breve{q}}_{b_i}^e = q_{b_i}^e \otimes [0;\frac{\mathbf{\omega}_i^b}{2}] + \frac{\epsilon_1}{2} ({q}_{b_i}^e \otimes \mathbf{a}_i^b + \mathbf{v}_i^e \otimes q_{b_i}^e(t) \otimes [0;\frac{\mathbf{\omega}_i^b}{2}]) + \\ \frac{\epsilon_2}{2} ({[0;\mathbf{v}}_i^e] \otimes {q}_{b_i}^e + \mathbf{t}_i^e \otimes q_{b_i}^e(t) \otimes [0;\frac{\mathbf{\omega}_i^b}{2}]). \end{aligned} q˘˙bie=qbie⊗[0;2ωib]+2ϵ1(qbie⊗aib+vie⊗qbie(t)⊗[0;2ωib])+2ϵ2([0;vie]⊗qbie+tie⊗qbie(t)⊗[0;2ωib]).
整理并化简上式可得 q ˘ ˙ b i e = q ˘ b i e ⊗ 1 2 ( ω i b + ϵ 1 a i b + ϵ 2 v i e ) , \dot{\breve{q}}_{b_i}^e = \breve{q}_{b_i}^e \otimes \frac{1}{2}(\boldsymbol{\omega}_i^b + \epsilon_1 \mathbf{a}_i^b + \epsilon_2 \mathbf{v}_i^e), q˘˙bie=q˘bie⊗21(ωib+ϵ1aib+ϵ2vie),
这里 ω i b \boldsymbol{\omega}_i^b ωib 、 a i b \mathbf{a}_i^b aib和 v i e \mathbf{v}_i^e vie 分别被表示为四元数形式,即实部为零。若记 ω ˘ i b = ( ω i b + ϵ 1 a i b + ϵ 2 v i e ) \breve{\omega }_i^b=(\boldsymbol{\omega}_i^b + \epsilon_1 \mathbf{a}_i^b + \epsilon_2 \mathbf{v}_i^e) ω˘ib=(ωib+ϵ1aib+ϵ2vie),那么鼎立四元数对时间的求导形式可与传统四元数相同: q ˘ ˙ b i e = q ˘ b i e ⊗ 1 2 ω ˘ i b . \dot{\breve{q}}_{b_i}^e = \breve{q}_{b_i}^e \otimes \frac{1}{2}\breve{\omega }_i^b. q˘˙bie=q˘bie⊗21ω˘ib.
综上,鼎立四元数将旋转、速度和平移统一表示在一个鼎立四元数中,是目前位置最简洁的姿态表示方式。
状态空间 | 运动学表示 |
---|---|
χ i = [ q ˘ b i e ] \mathbf{\chi}_i = [\breve{q}_{b_i}^e] χi=[q˘bie] | q ˘ b i e ⊗ 1 2 ω ˘ i b \breve{q}_{b_i}^e \otimes \frac{1}{2}\breve{\omega }_i^b q˘bie⊗21ω˘ib |
鼎立四元数运动学表示
TQ切比雪夫多项式迭代积分
在实际系统中,由于IMU采集的频率有限,因此数值积分存在计算误差。在常见的SLAM系统中,通常使用中点积分的方式来缓解积分数值误差。然而,在高速场景中,中点积分的性能也会迅速下降。论文中使用的切比雪夫多项式迭代积分方法可以实现更加准确的数值积分。接下来针对鼎立四元数的积分求解位姿进行讨论。
对于理想的位姿积分,积分过程为下式: q ˘ ( t ) = q ˘ ( 0 ) + ∫ 0 t q ˘ ˙ ( t ) d t . \breve{q}(t) = \breve{q}(0) + \int_{0}^{t} \dot{\breve{q}}(t)dt. q˘(t)=q˘(0)+∫0tq˘˙(t)dt. 根据2.4小节讨论的鼎立四元数动力学,积分过程可以近似为: q ˘ ( t ) = q ˘ ( 0 ) + 1 2 ∫ 0 t q ˘ ( t ) ⊗ ω ˘ b ( t ) d t . \breve{q}(t) = \breve{q}(0) + \frac{1}{2} \int_{0}^{t} {\breve{q}}(t) \otimes \breve{\omega }^b(t) dt. q˘(t)=q˘(0)+21∫0tq˘(t)⊗ω˘b(t)dt.
传统中点积分法
对于传统中点法,采用中点处的斜率直接近似两个采样点间的变化。其离散化推导如下: q ˘ i + 1 = q ˘ i + 1 4 q ˘ i ⊗ ( ω ˘ i + 1 b + ω ˘ i b ) . \breve{q}_{i+1} = \breve{q}_i + \frac{1}{4} {\breve{q}}_i \otimes (\breve{\omega }^b_{i+1} + \breve{\omega }^b_{i}). q˘i+1=q˘i+41q˘i⊗(ω˘i+1b+ω˘ib). 这种计算方式效率很高,但是只用两个时刻前后的采样值进行近似误差较大。适合用于短时间低速运动又需要快速推算的场景下。
切比雪夫逼近迭代法
切比雪夫迭代法的思路是通过切比雪夫多项式逼近真值运动。由于需要应用切比雪夫多项式逼近,首先需要将时间归一化至[-1,1]之间。假设采样持续时间为 t N t_N tN,每个采样时刻为 t k t_k tk,令 t = t N ( 1 + τ ) / 2 t = t_N(1+\tau)/2 t=tN(1+τ)/2,则式(37)可变换为 q ˘ ( τ ) = q ˘ ( 0 ) + t N 4 ∫ − 1 τ q ˘ ( τ ) ⊗ ω ˘ b ( τ ) d τ . \breve{q}(\tau) = \breve{q}(0) + \frac{t_N}{4} \int_{-1}^{\tau} {\breve{q}}(\tau) \otimes \breve{\omega }^b(\tau) d\tau. q˘(τ)=q˘(0)+4tN∫−1τq˘(τ)⊗ω˘b(τ)dτ. 首先,对于由角速度、加速度和速度构成的 ω ˘ b ( τ ) \breve{\omega }^b(\tau) ω˘b(τ)进行切比雪夫多项式逼近,得到 ω ˘ b ( τ ) = ∑ j = 0 n ω c ˘ j F j ( τ ) , \breve{\omega }^b(\tau) = \sum_{j=0}^{n_{\omega}} \breve{c}_j F_j(\tau), ω˘b(τ)=j=0∑nωc˘jFj(τ), 其中, n ω n_{\omega} nω为多项式阶数, c ˘ j \breve{c}_j c˘j为切比雪夫系数, F j ( τ ) F_j(\tau) Fj(τ)为 j j j阶切比雪夫多项式。
接着,对鼎立四元数所表示的状态进行切比雪夫逼近初始化,这里的初始化切比雪夫系数可以由传统中点法获得,得到 q ˘ ( τ ) = ∑ i = 0 n q b ˘ l , i F i ( τ ) , \breve{q }(\tau) = \sum_{i=0}^{n_{q}} \breve{b}_{l,i} F_i(\tau), q˘(τ)=i=0∑nqb˘l,iFi(τ), 其中, n q n_{q} nq为多项式阶数, b ˘ l , i \breve{b}_{l,i} b˘l,i为第 l l l次迭代的切比雪夫系数, F i ( τ ) F_i(\tau) Fi(τ)为 j j j阶切比雪夫多项式。将式(39)和式(40)代入式(41)可以得到 q ˘ l + 1 ( τ ) = q ˘ ( 0 ) + t N 4 ∫ − 1 τ ∑ i = 0 n q b ˘ l , i F i ( τ ) ⊗ ∑ j = 0 n ω c ˘ j F j ( τ ) d τ . \breve{q}_{l+1}(\tau) = \breve{q}(0) + \frac{t_N}{4} \int_{-1}^{\tau} \sum_{i=0}^{n_{q}} \breve{b}_{l,i} F_i(\tau) \otimes \sum_{j=0}^{n_{\omega}} \breve{c}_j F_j(\tau) d\tau. q˘l+1(τ)=q˘(0)+4tN∫−1τi=0∑nqb˘l,iFi(τ)⊗j=0∑nωc˘jFj(τ)dτ. 将积分和求和顺序进行变换后得到 q ˘ l + 1 ( τ ) = q ˘ ( 0 ) + t N 4 ∑ i = 0 n q b ˘ l , i ⊗ ∑ j = 0 n ω c ˘ j ∫ − 1 τ F i ( τ ) F j ( τ ) d τ . \breve{q}_{l+1}(\tau) = \breve{q}(0) + \frac{t_N}{4} \sum_{i=0}^{n_{q}} \breve{b}_{l,i} \otimes \sum_{j=0}^{n_{\omega}} \breve{c}_j \int_{-1}^{\tau} F_i(\tau) F_j(\tau) d\tau. q˘l+1(τ)=q˘(0)+4tNi=0∑nqb˘l,i⊗j=0∑nωc˘j∫−1τFi(τ)Fj(τ)dτ. 通过式(43)的计算得到 q l + 1 ( τ ) q_{l+1}(\tau) ql+1(τ)的值后,将结果代回到式(41),再次对状态进行切比雪夫逼近,重复此操作不断迭代直到收敛。
由于鼎立四元数统一表示了旋转、位置和速度,因此在切比雪夫迭代过程中,只需要进行一次迭代运算,而传统的表示方法需要分别对旋转、速度和平移进行3次迭代,提升了计算效率。