前言

向量叉乘

二维

a=[x1y1],b=[x2y2]⇓a×b=x1y2−y1x2a=\begin{bmatrix}x_1\\y_1\end{bmatrix},b=\begin{bmatrix}x_2\\y_2\end{bmatrix}\\\Downarrow\\a\times b=x_1y_2-y_1x_2

三维

a=[x1y1z1],b=[x2y2z2]⇓a×b=[y1z2−z1y2z1x2−x1z2x1y2−y1x2]a=\begin{bmatrix}x_1\\y_1\\z_1\end{bmatrix},b=\begin{bmatrix}x_2\\y_2\\z_2\end{bmatrix}\\\Downarrow\\a\times b=\begin{bmatrix}y_1z_2-z_1y_2\\z_1x_2-x_1z_2\\x_1y_2-y_1x_2\end{bmatrix}

欧拉角

是用三个角度表示物体的空间姿态，每个角度代表物体绕其一个固定轴的旋转量

  • pitch 俯仰角，一般来说向上转，也就是抬头为正
  • roll 翻滚角，一般是向右侧倾倒为正
  • yaw 航向角，一般是向右转为正

旋转矩阵

对于三轴角度的旋转矩阵，如下

  • 绕 x 轴旋转

    Rx(θ)=[1000cos⁡θ−sin⁡θ0sin⁡θcos⁡θ]R_x(\theta)=\begin{bmatrix}1&0&0\\0&\cos\theta&-\sin\theta\\0&\sin\theta&\cos\theta\end{bmatrix}

  • 绕 y 轴旋转

    Ry(γ)=[cos⁡γ0−sin⁡γ010sin⁡γ0cos⁡γ]R_y(\gamma)=\begin{bmatrix}\cos\gamma&0&-\sin\gamma\\0&1&0\\\sin\gamma&0&\cos\gamma\end{bmatrix}

  • 绕 z 轴旋转

    Rz(φ)=[cos⁡φ−sin⁡φ0sin⁡φcos⁡φ0001]R_z(\varphi)=\begin{bmatrix}\cos\varphi&-\sin\varphi&0\\\sin\varphi&\cos\varphi&0\\0&0&1\end{bmatrix}

为了表示一个物体绕着三个轴的旋转，可以简单的将三个单轴的旋转矩阵相乘

R=Rz(φ)Ry(γ)Rx(θ)R=R_z(\varphi)R_y(\gamma)R_x(\theta)

得到

R=[cos⁡γcos⁡φ−cosθsin⁡φ−cos⁡φsin⁡γsin⁡θsin⁡φsin⁡θ−cos⁡φcos⁡θsin⁡γcos⁡γsin⁡φcos⁡φcos⁡θ−sin⁡γsin⁡φsin⁡θ−cos⁡φsin⁡θ−cos⁡θsin⁡γsin⁡φ−sin⁡γcos⁡γsin⁡θcos⁡γcos⁡θ]R=\begin{bmatrix}\cos\gamma\cos\varphi&-cos\theta\sin\varphi-\cos\varphi\sin\gamma\sin\theta&\sin\varphi\sin\theta-\cos\varphi\cos\theta\sin\gamma\\\cos\gamma\sin\varphi&\cos\varphi\cos\theta-\sin\gamma\sin\varphi\sin\theta&-\cos\varphi\sin\theta-\cos\theta\sin\gamma\sin\varphi\\-\sin\gamma&\cos\gamma\sin\theta&\cos\gamma\cos\theta\end{bmatrix}

但是对于欧拉角的旋转，当任何一个坐标轴旋转角度为90度时，就会有两个轴的旋转动作起到对总体旋转结果相同的效果，这就被称为死锁

旋转矩阵与欧拉角

如上述公式，可以得到欧拉角与旋转矩阵的关系

γ=−arcsin⁡R31θ=arctan⁡R32R33φ=arctan⁡R21R11\gamma=-\arcsin R_{31}\\\theta=\arctan\frac{R_{32}}{R_{33}}\\\varphi=\arctan\frac{R_{21}}{R_{11}}

四元数

也是一种姿态的表示，相较于欧拉角，计算更加简便，通过绕空间中的任意一个轴，旋转一定角度，就可以得到姿态角的朝向。与欧拉角和旋转矩阵相比，四元数在某些计算中具有优势，如插值和归一化。此外，四元数避免了所谓的“万向节锁，仅需要存储 4 个浮点数，比矩阵更加轻量，无论是求逆，串联操作，都比矩阵更快

四元数的表达式为

q=q0+q1i+q2j+q3kq=q_0+q_1i+q_2j+q_3k

其中的虚向量满足

i2=j2=k2=ijk=−1i^2=j^2=k^2=ijk=-1

四元数可以理解为四维空间，三个虚部可理解为三维空间的三个正交基，而实部垂直于三维空间

单位四元数是四维空间中一个超球上面的点，满足 q02+q12+q22+q32=1q_0^2+q_1^2+q_2^2+q_3^2=1 ，而纯四元数是四维空间在 w=0w=0 时的一个子空间的点

单位四元数还有另一种表达公式

q=[cos⁡θ2sin⁡θ2v⃗]q=\begin{bmatrix}\cos\frac{\theta}{2}&\sin\frac{\theta}{2}\vec{v}\end{bmatrix}

共轭四元数指的是实部相同但是虚部相反的四元数，例如 qq 的共轭四元数为

q∗=[q0−q1−q2−q3]q^*=\begin{bmatrix}q_0&-q_1&-q_2&-q_3\end{bmatrix}

四元数微分运算

四元数微分通常用下式表达

dqdt=0.5⋅q⊗w\frac{dq}{dt}=0.5·q\otimes w

其中 ⊗\otimes 是四元数乘法

展开这个四元数乘法，写作矩阵形式

Q(q)=[q0−q1−q2−q3q1q0−q3q2q2q3q0−q1q3−q2q1q0]Q(q)=\begin{bmatrix}q_0&-q_1&-q_2&-q_3\\q_1&q_0&-q_3&q_2\\q_2&q_3&q_0&-q_1\\q_3&-q_2&q_1&q_0\end{bmatrix}

上述微分表达式可以写作

dqdt=0.5Q(q)[0wxwywz]\frac{dq}{dt}=0.5Q(q)\begin{bmatrix}0\\w_x\\w_y\\w_z\end{bmatrix}

上述四元数微分还有另一种表达方式

dqdt=0.5[0−wx−wy−wzwx0wz−wywy−wz0wxwzwy−wx0]q\frac{dq}{dt}=0.5\begin{bmatrix}0&-w_x&-w_y&-w_z\\w_x&0&w_z&-w_y\\w_y&-w_z&0&w_x\\w_z&w_y&-w_x&0\end{bmatrix}q

四元数与旋转矩阵

四元数与旋转矩阵之间可以相互转换

R=[1−2q22−2q322q1q2−2q0q32q1q3+2q0q22q1q2+2q0q31−2q12−2q322q2q3−2q0q12q1q3−2q0q22q2q3+2q0q11−2q12−2q22]R=\begin{bmatrix}1-2q_2^2-2q_3^2&2q_1q_2-2q_0q_3&2q_1q_3+2q_0q_2\\2q_1q_2+2q_0q_3&1-2q_1^2-2q_3^2&2q_2q_3-2q_0q_1\\2q_1q_3-2q_0q_2&2q_2q_3+2q_0q_1&1-2q_1^2-2q_2^2\end{bmatrix}

旋转矩阵需要九个参数，相比于四元数，参数量更多。但相比于四元数，旋转矩阵好处是更直观，代表着机体坐标系的三个轴的方向

四元数与欧拉角

通过上述旋转矩阵，可以得到欧拉角表示下的四元数

q=[cos⁡θ2cos⁡γ2cos⁡φ2+sin⁡θ2sin⁡γ2sin⁡φ2cos⁡θ2sin⁡γ2cos⁡φ2+sin⁡θ2cos⁡γ2sin⁡φ2cos⁡θ2cos⁡γ2sin⁡φ2−sin⁡θ2sin⁡γ2cos⁡φ2sin⁡θ2cos⁡γ2cos⁡φ2−cos⁡θ2sin⁡γ2sin⁡φ2]q=\begin{bmatrix}\cos\frac{\theta}{2}\cos\frac{\gamma}{2}\cos\frac{\varphi}{2}+\sin\frac{\theta}{2}\sin\frac{\gamma}{2}\sin\frac{\varphi}{2}\\\cos\frac{\theta}{2}\sin\frac{\gamma}{2}\cos\frac{\varphi}{2}+\sin\frac{\theta}{2}\cos\frac{\gamma}{2}\sin\frac{\varphi}{2}\\\cos\frac{\theta}{2}\cos\frac{\gamma}{2}\sin\frac{\varphi}{2}-\sin\frac{\theta}{2}\sin\frac{\gamma}{2}\cos\frac{\varphi}{2}\\\sin\frac{\theta}{2}\cos\frac{\gamma}{2}\cos\frac{\varphi}{2}-\cos\frac{\theta}{2}\sin\frac{\gamma}{2}\sin\frac{\varphi}{2}\end{bmatrix}

欧拉角与四元数关系图

同理可以得到四元数表示下的欧拉角

[θγφ]=[arctan⁡2q2q3+2q0q11−2q12−2q22−arcsin⁡(2q1q3−2q0q2)arctan⁡2q1q2+2q0q31−2q22−2q32]\begin{bmatrix}\theta\\\gamma\\\varphi\end{bmatrix}=\begin{bmatrix}\arctan\frac{2q_2q_3+2q_0q_1}{1-2q_1^2-2q_2^2}\\-\arcsin (2q_1q_3-2q_0q_2)\\\arctan\frac{2q_1q_2+2q_0q_3}{1-2q_2^2-2q_3^2}\end{bmatrix}

四元数运算

对于四元数的叉乘，在 matlab 中是不支持四元数的叉乘的，所以就有了如下两种定义

Q1⊗Q2⇓[Q1]L⋅Q2⇓[q0−q1−q2−q3q1q0−q3q2q2q3q0−q1q3−q2q1q0]⋅Q2Q_1\otimes Q_2\\\Downarrow\\\begin{bmatrix}Q_1\end{bmatrix}_L· Q_2\\\Downarrow\\\begin{bmatrix}q_0&-q_1&-q_2&-q_3\\q_1&q_0&-q_3&q_2\\q_2&q_3&q_0&-q_1\\q_3&-q_2&q_1&q_0\end{bmatrix}·Q_2

Q1⊗Q2⇓[Q2]R⋅Q1⇓[q0−q1−q2−q3q1q0q3−q2q2−q3q0q1q3q2−q1q0]⋅Q1Q_1\otimes Q_2\\\Downarrow\\\begin{bmatrix}Q_2\end{bmatrix}_R· Q_1\\\Downarrow\\\begin{bmatrix}q_0&-q_1&-q_2&-q_3\\q_1&q_0&q_3&-q_2\\q_2&-q_3&q_0&q_1\\q_3&q_2&-q_1&q_0\end{bmatrix}·Q_1

坐标系之间的映射

对于在 Earth 参考系下的任意一个向量

Ed=[0dxdydz]^Ed=\begin{bmatrix}0&d_x&d_y&d_z\end{bmatrix}

在 Sensor 坐标系下通过传感器读出的数据为

Sd=ESq∗×Ed×ESq^Sd=^S_Eq^*\times {^Ed}\times {^S_Eq}

Mahony滤波

这一部分参考于 陀螺仪姿态解算+mahony滤波算法（公式推导及其代码） - 知乎 (zhihu.com)

动机

只有角速度积分所得到的角度会漂移，需要有绝对真实的加速度去修正积分带来的误差，所以先将欧拉角的角速度得到四元数的角速度，然后积分得到当前四元数，再转化为欧拉角。

根据四元数得到重力方向

gx=2(q1q3−q0q2)gy=2(q0q1+q2q3)gz=q0q0−q1q1−q2q2+q3q3g_x=2(q_1q_3-q_0q_2)\\g_y=2(q_0q_1+q_2q_3)\\g_z=q_0q_0-q_1q_1-q_2q_2+q_3q_3

计算积分带来的误差

根据上述可以得到仅有重力加速度的加速度矢量

g=[gxgygz]g=\begin{bmatrix}g_x\\g_y\\g_z\end{bmatrix}

并且我们可以根据角速度积分得到的欧拉角和加速度所计算出的加速度矢量

a=[axayaz]a=\begin{bmatrix}a_x\\a_y\\a_z\end{bmatrix}

叉乘得到

∣a∣∣g∣sin⁡e=a×g|a||g|\sin e=a\times g

小角度近似 sin⁡ψ=ψ\sin\psi=\psi ，所以就是

e=a×g∣a∣∣g∣e=\frac{a\times g}{|a||g|}

这就能看出来积分得到的数据与真实数据之间的误差了

PI控制器修正误差

PI控制器可以看作是一个低通滤波器，它用于整合加速度计的低频信息并与陀螺仪读数互补，从而实现互补滤波的效果

[θ˙γ˙φ˙]=[θ˙γ˙φ˙]+kp[exeyez]+ki[exeyez]△t\begin{bmatrix}\dot\theta\\\dot\gamma\\\dot\varphi\end{bmatrix}=\begin{bmatrix}\dot\theta\\\dot\gamma\\\dot\varphi\end{bmatrix}+k_p\begin{bmatrix}e_x\\e_y\\e_z\end{bmatrix}+k_i\begin{bmatrix}e_x\\e_y\\e_z\end{bmatrix}\triangle t

得到当前四元数

由上述可知

dqdt=0.5[0−wx−wy−wzwx0wz−wywy−wz0wxwzwy−wx0]q\frac{dq}{dt}=0.5\begin{bmatrix}0&-w_x&-w_y&-w_z\\w_x&0&w_z&-w_y\\w_y&-w_z&0&w_x\\w_z&w_y&-w_x&0\end{bmatrix}q

则有

qk+1=0.5△t[0−wx−wy−wzwx0wz−wywy−wz0wxwzwy−wx0]qk+qkq_{k+1}=0.5\triangle t\begin{bmatrix}0&-w_x&-w_y&-w_z\\w_x&0&w_z&-w_y\\w_y&-w_z&0&w_x\\w_z&w_y&-w_x&0\end{bmatrix}q_k+q_k

代码

根据上述公式，可以得到如下 C 代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#define DELTAT 0.005f

typedef struct {
  float q0;
  float q1;
  float q2;
  float q3;
} Quaternion;

typedef struct {
  float roll;   // x
  float pitch;  // y
  float yaw;    // z
} Euler;

typedef struct {
  float x;
  float y;
  float z;
} Acceleration;

// 暂定，需要调试才能确定
Acceleration MahonyKp = { 0.1, 1, 1 };
Acceleration MahonyKi = { 0.1, 1, 1 };

void QuaternionToEuler(Quaternion _q, Euler* _e) {
  float rolly = 2 * _q.q2 * _q.q3 + 2 * _q.q0 * _q.q1;
  float rollx = 1 - 2 * _q.q1 * _q.q1 - 2 * _q.q2 * _q.q2;
  _e->roll    = atan2f(rolly, rollx);
  _e->pitch   = asin(2 * _q.q1 * _q.q3 - 2 * _q.q0 * _q.q2);
  float yawy  = 2 * _q.q1 * _q.q2 + 2 * _q.q0 * _q.q3;
  float yawx  = 1 - 2 * _q.q2 * _q.q2 - 2 * _q.q3 * _q.q3;
  _e->yaw     = atan2f(yawy, yawx);
}

void EulerToQuaternion(Euler _e, Quaternion* _q) {
  float cosRoll  = cosf(_e.roll / 2);
  float sinRoll  = sinf(_e.roll / 2);

  float cosPitch = cosf(_e.pitch / 2);
  float sinPitch = sinf(_e.pitch / 2);

  float cosYaw   = cosf(_e.yaw / 2);
  float sinYaw   = sinf(_e.yaw / 2);

  _q->q0         = cosRoll * cosPitch * cosYaw + sinRoll * sinPitch * sinYaw;
  _q->q1         = cosRoll * sinPitch * cosYaw + sinRoll * cosPitch * sinYaw;
  _q->q2         = cosRoll * cosPitch * sinYaw - sinRoll * sinPitch * cosYaw;
  _q->q3         = sinRoll * cosPitch * cosYaw - cosRoll * sinPitch * sinYaw;
}

void Mahony(Acceleration _a, Euler* _e, Quaternion* _q, Euler _dote) {
  float anorm     = sqrtf(_a.x * _a.x + _a.y * _a.y + _a.z * _a.z);
  _a.x           /= anorm;
  _a.y           /= anorm;
  _a.z           /= anorm;

  Acceleration g  = {
    2 * (_q->q1 * _q->q3 - _q->q0 * _q->q2),
    2 * (_q->q0 * _q->q1 + _q->q2 * _q->q3),
    _q->q0 * _q->q0 - _q->q1 * _q->q1 - _q->q2 * _q->q2 + _q->q3 * _q->q3
  };

  Acceleration err = {
    _a.y * g.z - _a.z * g.y,
    _a.z * g.x - _a.x * g.z,
    _a.x * g.y - _a.y * g.x,
  };

  _dote.roll  += MahonyKp.x * err.x + MahonyKi.x * err.x * DELTAT;
  _dote.pitch += MahonyKp.y * err.y + MahonyKi.y * err.y * DELTAT;
  _dote.yaw   += MahonyKp.z * err.z + MahonyKi.z * err.z * DELTAT;

  Quaternion qdot;
  qdot.q0  = -_dote.roll * _q->q1 - _dote.pitch * _q->q2 - _dote.yaw * _q->q3;
  qdot.q1  = _dote.roll * _q->q0 + _dote.yaw * _q->q2 - _dote.pitch * _q->q3;
  qdot.q2  = _dote.pitch * _q->q0 - _dote.yaw * _q->q1 + _dote.roll * _q->q3;
  qdot.q3  = _dote.yaw * _q->q0 + _dote.pitch * _q->q1 - _dote.roll * _q->q2;

  _q->q0  += 0.5 * DELTAT * qdot.q0;
  _q->q1  += 0.5 * DELTAT * qdot.q1;
  _q->q2  += 0.5 * DELTAT * qdot.q2;
  _q->q3  += 0.5 * DELTAT * qdot.q3;

  QuaternionToEuler(*_q, _e);
}

总结

上述的 Mahony 滤波算法实际上是基于就静止状态下的，只有在静止状态下才能使得机体的加速度 aa 是只有重力加速度反馈的，所以需要注意，当 ∣a∣|a| 的大小不等于重力加速度时，置信度应该降低，对这一部分可以采用自适应的扩展卡尔曼滤波来做权重的分配

Madgwick融合算法

v2-9ac058e63ce0e1bee5edf9069df0cb1d_720w.webp

通过陀螺仪的加速度初步计算传感器方向的偏移 ESqw,t^S_Eq_{w,t} 然后通过对重力加速度和地磁这两个恒定量在传感器坐标系下的投影计算传感器方向的偏移 ESq▽,t^S_Eq_{\triangledown,t}。最后通过互补滤波器融合 ESqw,t^S_Eq_{w,t} 和 ESq▽,t^S_Eq_{\triangledown,t} 来得到一个更加可信的结果，此外还要矫正地磁计数据畸变和陀螺仪零点漂移。其中 ESq^S_Eq 表示坐标系 Earth 相对于坐标系 Sensor 的单位旋转四元数

计算 ESqw,t^S_Eq_{w,t}

从陀螺仪中读取到的数据为

Swt=[0wxwywz]^Sw_t=\begin{bmatrix}0&w_x&w_y&w_z\end{bmatrix}

地面参考系相对于传感器参考系的旋转可以由其四元数的导数对时间的积分获得，由上述可知

ESqw,t=ESq^est,t−1+(12ESq^est,t−1⊗Swt)△t=ESq^est,t−1+ESq˙w,t△t{^S_Eq_{w,t}}={^S_E\hat{q}_{est,t-1}}+(\frac{1}{2}{^S_E\hat{q}_{est,t-1}}\otimes {^Sw_t})\triangle t={^S_E\hat{q}_{est,t-1}}+{^S_E\dot{q}_{w,t}}\triangle t

其中 ESq^est,t−1{^S_E\hat{q}_{est,t-1}} 是上一时刻的估计值，这个值就是根据角速度所 ESqw,t^S_Eq_{w,t}

计算 ESq▽,t^S_Eq_{\triangledown,t}

如果一个地球场的方向在地球框架内是已知的，那么在传感器框架内测量该场的方向就可以计算出传感器框架相对于地球框架的方向。然而，对于任何给定的测量，将不会有唯一的传感器方向解，相反，将有无限的解由所有这些方向表示，这些方向是通过将真实方向绕与场平行的轴旋转来实现的。

者是利用外部恒定的向量在各个角度的投影来估算传感器的位姿，在 IMU 中利用重力加速度，在 MARG 中同时利用重力加速度和地磁方向来计算。对于在 Earth 参考系下的任意一个向量

Ed=[0dxdydz]^Ed=\begin{bmatrix}0&d_x&d_y&d_z\end{bmatrix}

在 Sensor 坐标系下通过传感器直接读出的数据为 Ss^Ss ，还有就是通过传感器所测量出的四元数转换之后所得到的向量为

Sd=ESq∗×Ed×ESq^Sd=^S_Eq^*\times {^Ed}\times {^S_Eq}

其中 ESq∗^S_Eq^* 为 ESq{^S_Eq} 的共轭四元数。但是这两者是通过不同的方式得到的，也就有了估计 ESq^S_Eq 的一个方法，就是构建目标函数

min⁡f=ESq∗⊗Ed⊗ESq−Ss\min f=^S_Eq^*\otimes {^Ed}\otimes {^S_Eq}-{^Ss}

使得 ff 取得最小值时的 ESq^S_Eq 就是最终结果（最优结果），优化算法可以使用梯度下降法，通过迭代公式来计算

ESqk+1=ESqk−μ▽f∥▽f∥^S_Eq_{k+1}=^S_Eq_{k}-\mu \frac{\triangledown f}{\lVert\triangledown f\rVert}

其中 uu 就是学习率，而 ▽f\triangledown f 就是 ff 的梯度。计算由目标函数及其雅克比定义的解曲面的梯度

▽f=JT(ESqk,Ed)f\triangledown f=J^T({^S_Eq_k,{^Ed}})f

其中

J=∂f∂ESqkJ=\frac{\partial f}{\partial {^S_Eq_k}}

上述所得到的公式适用于空间中的任意一个方向的向量，太复杂了，不好求解，看其中的思路就是把叉乘的列向量扩展到了一个方阵，类似于上式中求解四元数的加速度的方程

1711678115108.png

1711678667170.png

IMU

如果使用的是 IMU ，那其中的 Ed^Ed 就是重力加速度向量 [0001]\begin{bmatrix}0&0&0&1\end{bmatrix}

并且 Ss^Ss 就是 IMU 传感器所测出来的在 Sensor 坐标下的重力加速度的向量 [0axayaz]\begin{bmatrix}0&a_x&a_y&a_z\end{bmatrix} ，带入上述公式可以得到

f=[2(q2q3−q1q3)−ax2(q1q2+q3q4)−ay2(12−q22−q32)−az]J=[−2q32q4−2q12q22q22q12q42q30−4q2−4q30]f=\begin{bmatrix}2(q_2q_3-q_1q_3)-a_x\\2(q_1q_2+q_3q_4)-a_y\\2(\frac{1}{2}-q_2^2-q_3^2)-a_z\end{bmatrix}\\J=\begin{bmatrix}-2q_3&2q_4&-2q_1&2q_2\\2q_2&2q_1&2q_4&2q_3\\0&-4q_2&-4q_3&0\end{bmatrix}

MARG

如果使用的是 MARG，在考虑重力加速度的基础上还需要将当地的地磁传感器的读数考虑进来。

Ed=[0bx0bz]^Ed=\begin{bmatrix}0&b_x&0&b_z\end{bmatrix}

这里是纯地磁传感器所得到的数据

Sm=[0mxmymz]\\{^Sm}=\begin{bmatrix}0&m_x&m_y&m_z\end{bmatrix}

再与上述的加速度数据融合，也就是

Sa=[0axayaz]{^Sa}=\begin{bmatrix}0&a_x&a_y&a_z\end{bmatrix}

写作一个矩阵，就是

f=[fg(ESq,Sa)fd(ESq,Ed,Sm)]f=\begin{bmatrix}f_g({^S_Eq,{^Sa}})\\f_d({^S_Eq,{^Ed},{^Sm}})\end{bmatrix}

对应的雅各比矩阵就是

J=[JgT(ESq)JdT(ESq,Ed)]J=\begin{bmatrix}J^T_g({^S_Eq})\\J_d^T({^S_Eq},{^Ed})\end{bmatrix}

根据迭代公式

ESq▽,t=ESq^est,t−1−μt▽f∥▽f∥^S_Eq_{\triangledown,t}={^S_E\hat{q}_{est,t-1}}-\mu_t \frac{\triangledown f}{\lVert\triangledown f\rVert}

其中

▽f={JgT(ESq^est,t−1)fg(ESq^est,t−1,Sat)JT(ESq^est,t−1,Eb^t)f(ESq^est,t−1,Sat,Ebt,Smt)\triangledown f=\left\{\begin{aligned}&J^T_g({^S_E\hat{q}_{est,t-1}})f_g({^S_E\hat{q}_{est,t-1}},{^Sa_t})\\&J^T({^S_E\hat{q}_{est,t-1}},{^E\hat{b}_t})f({^S_E\hat{q}_{est,t-1}},{^Sa_t},{^Eb_t},{^Sm_t})\end{aligned}\right.

学习率的选择需要满足

μt=α∥ESq˙w,t∥△t\mu_t=\alpha \lVert{^S_E\dot q_{w,t}}\rVert\triangle t

其中 α>1\alpha>1

融合滤波算法

将上述中计算出的方向四元数 ESqw,t^S_Eq_{w,t} 和 ESq▽,t^S_Eq_{\triangledown,t} 进行融合，可以估计得到传感器相对于地球坐标系的方向 ESqest,t^S_Eq_{est,t} ，如下

ESqest,t=γq▽,t+(1−γ)ESqw,t^S_Eq_{est,t}=\gamma q_{\triangledown,t}+(1-\gamma){^S_Eq_{w,t}}

其中 0≤γ≤10\leq \gamma\leq 1 ，其中 γ\gamma 的最优值可以定义为确保 ESqw,t^S_Eq_{w,t} 的加权发散度等于 ESq▽,t^S_Eq_{\triangledown,t} 的加权收敛度。可以将 γ\gamma 定义为如下公式

(1−γ)β=γμtΔγ=βμtΔ+β(1-\gamma)\beta=\gamma\frac{\mu_t}{\Delta}\\\gamma=\frac{\beta}{\frac{\mu_t}{\Delta}+\beta}

假设上述中由 α\alpha 控制的 ESq▽^S_Eq_{\triangledown} 的收敛速度等于或大于实际方向变化速率，上述公式就保证了 ESqw,t^S_Eq_{w,t} 和 ESq▽,t^S_Eq_{\triangledown,t} 的最优融合。

由上述可知 α\alpha 没有上界，如果 α\alpha 非常大时，那么上述公式中所定义的 μt\mu_t 就非常大，那这个滤波也就简化了，这意味着 ESq^est,t−1{^S_E\hat{q}_{est,t-1}} 可以忽略不计，所以公式可以写作

ESq▽,t=−μt▽f∥▽f∥^S_Eq_{\triangledown,t}=-\mu_t \frac{\triangledown f}{\lVert\triangledown f\rVert}

相应的公式中的 γ\gamma 也简化了，公式中的 β\beta 项可以忽略不计，因此可写作

γ=βΔtμt\gamma=\frac{\beta\Delta t}{\mu_t}

甚至可以假设 γ=0\gamma=0 。将上述得到的带入融合滤波公式，可以得到

ESqest,t=γq▽,t+(1−γ)ESqw,t=βΔtμt(−μt▽f∥▽f∥)+(1−0)(ESq^est,t−1+ESq˙w,t△t)=ESq^est,t−1+(ESq˙w,t−β▽f∥▽f∥)Δt^S_Eq_{est,t}=\gamma q_{\triangledown,t}+(1-\gamma){^S_Eq_{w,t}}\\=\frac{\beta\Delta t}{\mu_t}(-\mu_t \frac{\triangledown f}{\lVert\triangledown f\rVert})+(1-0)({^S_E\hat{q}_{est,t-1}}+{^S_E\dot{q}_{w,t}}\triangle t)\\={^S_E\hat{q}_{est,t-1}}+({^S_E\dot{q}_{w,t}}-\beta\frac{\triangledown f}{\lVert\triangledown f\rVert})\Delta t

上式就可以化简为

ESqest,t=ESq^est,t−1+ESq˙est,tΔt^S_Eq_{est,t}={^S_E\hat{q}_{est,t-1}}+{^S_E\dot{q}_{est,t}}\Delta t

其中

ESq˙est,t=ESq˙w,t−βESq^˙ε,tESq^˙ε,t=▽f∥▽f∥^S_E\dot{q}_{est,t}={^S_E}\dot{q}_{w,t}-\beta{^S_E\dot{\hat{q}}_{\varepsilon,t}}\\{^S_E\dot{\hat{q}}_{\varepsilon,t}}=\frac{\triangledown f}{\lVert\triangledown f\rVert}

可以看出，滤波器通过对估计的方向四元数导数 ESq˙est^S_E\dot{q}_{est} 进行数值积分得到方向 ESq˙est^S_E\dot{q}_{est} 。滤波器将其计算为陀螺仪测量方向变化率 ESq˙w^S_E\dot{q}_{w} 。有了陀螺仪测量误差 β\beta ，然后再估计误差方向上减去从加速度计和磁力计测量计算得到的 ESq^˙ε{^S_E\dot{\hat{q}}_{\varepsilon}}

代码

由于 matlab 中并不支持四元数的叉乘，所以写了一个函数来实现四元数的叉乘

首先是四元数转化为方阵形式，根据上述的说明，有两种转化方式

1
2
3
4
5
6
7
8
9
10
11
12
13
function res = transR(q1)
    res = [q1(1) -q1(2) -q1(3) -q1(4);
           q1(2) q1(1) q1(4) -q1(3);
           q1(3) -q1(4) q1(1) q1(2);
           q1(4) q1(3) -q1(2) q1(1);];
end
function res = transL(q1) 
    res = [q1(1) -q1(2) -q1(3) -q1(4);
           q1(2) q1(1) -q1(4) q1(3);
           q1(3) q1(4) q1(1) -q1(2);
           q1(4) -q1(3) q1(2) q1(1);];
    
end

所以此传感器的代码为（未加磁偏角）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
clc
clear

syms q1(t) q2(t) q3(t) q0(t) dx dy dz ax ay az bx bz mx my mz

q = [q0(t) q1(t) q2(t) q3(t)]';
q_conj = [q0(t) -q1(t) -q2(t) -q3(t)]';
g = [0 0 0 1]';
a = [0 ax ay az]';

fg = transR(q) * transL(q_conj) * g - a;
fg = subs(fg, [conj(ax), conj(ay), conj(az)], [ax, ay, az]);
fg = subs(fg, [conj(q0(t)), conj(q1(t)), conj(q2(t)), conj(q3(t))], [q0(t), q1(t), q2(t), q3(t)]);
fg = simplify(fg);
Jg = jacobian(fg, [q0; q1; q2; q3]);
Jg = simplify(Jg);

b = [0 bx 0 bz]';
m = [0 mx my mz]';

fb = transR(q) * transL(q_conj) * b - m;
fb = subs(fb, [conj(mx), conj(my), conj(mz)], [mx, my, mz]);
fb = subs(fb, [conj(bx), conj(bz)], [bx, bz]);
fb = subs(fb, [conj(q0(t)), conj(q1(t)), conj(q2(t)), conj(q3(t))], [q0(t), q1(t), q2(t), q3(t)]);
fb = simplify(fb);
Jb = jacobian(fb, [q0; q1; q2; q3]);
Jb = simplify(Jb);

Fg = Jg' * fg;
Fg = subs(Fg, [conj(ax), conj(ay), conj(az)], [ax, ay, az]);
Fg = subs(Fg, [conj(q0(t)), conj(q1(t)), conj(q2(t)), conj(q3(t))], [q0(t), q1(t), q2(t), q3(t)]);

Fb = [Jg', Jb'] * [fg; fb];
Fb = subs(Fb, [conj(mx), conj(my), conj(mz)], [mx, my, mz]);
Fb = subs(Fb, [conj(bx), conj(bz)], [bx, bz]);
Fb = subs(Fb, [conj(q0(t)), conj(q1(t)), conj(q2(t)), conj(q3(t))], [q0(t), q1(t), q2(t), q3(t)]);

Fgnorm = norm(Fg);
Fbnorm = norm(Fb);

beta = 1;
DeltaT = 0.001;

syms w0 w1 w2 w3
qw_dot = [w0 w1 w2 w3]';

q1 = q + (qw_dot  - beta * Fg / Fgnorm) * DeltaT
q2 = q + (qw_dot  - beta * Fb / Fbnorm) * DeltaT

其中 FgFg 和 FbFb 分别是在 IMU 和 MARG 条件下的 ▽f\triangledown f ，最终得到的 q1q_1 和 q2q_2 就是该条件下滤波的最终结果