Madgwick算法的数学推演：基于惯性导航元件的位姿估计

《An efficient orientation filter for inertial and inertial/magnetic sensor arrays》 Sebastian O.H. Madgwick

• 文章简介：这篇文章详细的介绍了基于惯性导航元件位姿估计的Madgwick算法的推导过程并将实验结果和基于卡尔曼滤波算法的结果比对，证明了Madgwick算法的优秀性能。在摘要中主要提到了三个亮点：
  1) Wadgwick算法对于IMU（Inertial Measurement Unit）每次迭代仅使用109次运算，对于MARG（IMU + magntetometers）传感器使用277次运算，非常适合开发轻量级惯导设备。
  2) Wadgwick算法在低采样频率的情况下仍保持良好性能。
  3) 对于IMU传感器有一个可调参数，对于MARG传感器有仅有两个可调参数。下面对文章中的数学过程进行逐一翻译和梳理。
   
• 需要提前掌握的数学内容：
  四元数的概念，四元数的乘法以及含义，四元数的共轭，用四元数计算向量的旋转，四元数与旋转矩阵的转换，四元数与欧拉角的转换。文章中的四元数$_B^A\bm{\hat q}_{}$BA​q^​​含义为：坐标系B相对于坐标系A的单位旋转四元数。
   
• 算法核心思想：
  Wadgwick算法的核心思想如下：首先，通过陀螺仪的角速度初步计算传感器方向偏移 $_E^S\bm{q}_{\omega,t}$ES​qω,t​；然后，通过对重力加速度和地磁这两个恒定量在传感器坐标系下的投影计算传感器方向偏移 $_E^S\bm{q}_{\nabla,t}$ES​q∇,t​；最后通过互补滤波器融合 $_E^S\bm{q}_{\omega,t}$ES​qω,t​ 和 $_E^S\bm{q}_{\nabla,t}$ES​q∇,t​ 来得到一个更加可信的结果。除此之外还需要校正地磁计数据畸变和陀螺仪零点漂移。我们下面会一步一步看各部分的数学推演。
   
• 通过陀螺仪的角速度对$_E^S\bm{q}_{\omega,t}$ES​qω,t​的计算：
  从陀螺仪中直接读取的数据为$^S\bm{\omega}=[0,{\omega_x},{\omega_y},{\omega_z}]$Sω=[0,ωx​,ωy​,ωz​]。地面参考系（Earth）相对于传感器参考系（Sensor）的旋转可以由其四元数的导数对时间的积分获得。
   
  $_E^S\bm{q}_{\omega,t}=_E^S\bm{\hat q}_{est,t-1}+_E^S\bm{\dot q}_{\omega,t}\Delta t$ES​qω,t​=ES​q^​est,t−1​+ES​q˙​ω,t​Δt
   
  其中$_E^S\bm{\hat q}_{est,t-1}$ES​q^​est,t−1​是在t-1时刻得到的最佳单位位姿四元数。而$_E^S\bm{\dot q}_{\omega,t}$ES​q˙​ω,t​可以用David et al.提出的Npsnet中的数学模型进行运算，公式如下：
   
  $_E^S\bm{\dot q}_{\omega,t}=\frac{1}{2}_E^S\bm{\hat q}_{est,t-1}×^S\bm{\omega}_t$ES​q˙​ω,t​=21​ES​q^​est,t−1​×Sωt​
   
  公式中的角速度可以直接从传感器读出，上一时刻的位姿也是已知，时间间隔也可以从传感器芯片中的获得。可以轻松的计算出$_E^S\bm{q}_{\omega,t}$ES​qω,t​的数值结果。
   
• 通过恒定向量对$_E^S\bm{q}_{\nabla,t}$ES​q∇,t​的计算：
  这一部分是利用外部恒定的向量在各个角度的投影来估算传感器的位姿。在IMU中将利用重力加速度，而在MARG中将同时利用重力加速度和地磁方向来计算。对于在Earth参考系下的恒定向量$^E\bm{\hat d}=[0,d_x,d_y,d_z]$Ed^=[0,dx​,dy​,dz​]，在Sensor坐标系下通过传感器读出的数据为$^S\bm{\hat s}=[0,s_x,s_y,s_z]$Ss^=[0,sx​,sy​,sz​]。在理论上，这二者存在如下关系：
   
  $^S\bm{\hat s}=_E^S\bm{\hat q}^*×^E\bm{\hat d}×_E^S\bm{\hat q}$Ss^=ES​q^​∗×Ed^×ES​q^​
   
  但是这二者是通过不同的方式获得的，其间的等式关系基本不存在，这也恰好提供了一个估计 $_E^S\bm{\hat q}$ES​q^​ 的一个好方法，就是：构建优化目标函数 $f=_E^S\bm{\hat q}^*×^E\bm{\hat d}×_E^S\bm{\hat q}-^S\bm{\hat s}$f=ES​q^​∗×Ed^×ES​q^​−Ss^，在 $f$f 取得最小值时的 $_E^S\bm{\hat q}$ES​q^​ 就是我们想要的结果。优化算法使用的是大名鼎鼎的梯度下降法，迭代公式如下：
   
  $_E^S\bm{q}_{k+1}=_E^S\bm{q}_{k}-\mu{\frac{\nabla \bm f}{\|{\nabla \bm f}\|}}$ES​qk+1​=ES​qk​−μ∥∇f∥∇f​
   
  其中$\mu$μ是学习率，${\nabla \bm f}$∇f在论文中也给出了解析表达式，可以方便查阅。如果使用的是IMU，仅需要将其中的$^E\bm{\hat d}$Ed^看成重力加速度$[0,0,0,1]$[0,0,0,1]，$^S\bm{\hat s}$Ss^就是加速度计的单位四元数。若使用的是MARG，在考虑重力加速度的基础上还需要将当地的准确地磁和地磁传感器的读出考虑进来。最终通过梯度下降的方法就能够得到$_E^S\bm{q}_{\nabla,t}$ES​q∇,t​的结果。值得一提的是，学习率的选择上需要满足条件：
   
  $\mu_t=\alpha\|_E^S\bm{\dot q}_{\omega,t}\|\Delta t$μt​=α∥ES​q˙​ω,t​∥Δt，其中$\alpha>1$α>1
   
  上述式子的含义是$_E^S\bm{q}_{\nabla,t}$ES​q∇,t​的收敛速度需要大于$_E^S\bm{q}_{\omega,t}$ES​qω,t​的发散速度。其中$_E^S\bm{\dot q}_{\omega,t}$ES​q˙​ω,t​就是$_E^S\bm{q}_{\omega,t}$ES​qω,t​的发散速度，而$\frac {\mu_t}{\Delta t}$Δtμt​​则是通过梯度下降运算中的$_E^S\bm{q}_{\nabla,t}$ES​q∇,t​的收敛速度。
   
• 通过互补滤波方法融合 $_E^S\bm{q}_{w,t}$ES​qw,t​ 和 $_E^S\bm{q}_{\nabla,t}$ES​q∇,t​ 的结果：
  互补滤波的方法非常简单，就是这两个部分的加权平均而已：
   
  $_E^S\bm{q}_{est,t}={\gamma_t}_E^S\bm{q}_{\nabla,t}+(1-\gamma_t)_E^S\bm{q}_{w,t}$ES​qest,t​=γt​ES​q∇,t​+(1−γt​)ES​qw,t​，其中$0&lt;\gamma_t&lt;1$0<γt​<1
   
  如果假设 $\beta$β 是 $_E^S\bm{q}_{\omega,t}$ES​qω,t​ 的发散速度，一般用陀螺仪的误差四元数导数的模来表示。如果 $_E^S\bm{q}_{\omega,t}$ES​qω,t​ 的发散速度和 $_E^S\bm{q}_{\nabla,t}$ES​q∇,t​ 的收敛速度相等，则$\gamma_t$γt​有：
   
  $(1-\gamma_t)\beta={\gamma_t}{\frac{\mu_t}{\Delta t}}$(1−γt​)β=γt​Δtμt​​
   
  进而得到 $\gamma_t={\frac{\beta}{\frac{\mu_t}{\Delta t}+\beta}}$γt​=Δtμt​​+ββ​
   
  假设 $_E^S\bm{q}_{\nabla}$ES​q∇​ 的收敛速度突破了物理极限，则 $\alpha$α 会远远大于 $\mu_t$μt​，这会带来一些近似的化简（很简单，这里不详细列出）。这些化简步骤将会进一步减少运算。最终我们得到的迭代公式如下：
   
  $_E^S\bm q_{est,t}=_E^S\bm {\hat q}_{est,t}+(_E^S\bm {\dot q}_{\omega,t}-\beta\frac{\nabla \bm f}{\|\nabla \bm f\|})\Delta t$ES​qest,t​=ES​q^​est,t​+(ES​q˙​ω,t​−β∥∇f∥∇f​)Δt
   
• 地磁计数据畸变校正：
  地磁计的畸变产生是由于周围的电磁波的干扰。在传感器坐标系下静止的外界电磁干涉源可以通过地磁计的校准进行校正剔除。而在Earth坐标系下的干扰源将会主要产生地磁的测量偏差。在t时刻的地磁方向$^E\bm {\hat h}_t$Eh^t​可以通过地磁计的读数$^S\bm {\hat m}_t$Sm^t​计算：
   
  $^E\bm {\hat h}_t=[0,h_x,h_y,h_z]=_E^S\bm{\hat q}_{est,t-1}×^S\bm {\hat m}_t×_E^S\bm{\hat q}^*_{est,t-1}$Eh^t​=[0,hx​,hy​,hz​]=ES​q^​est,t−1​×Sm^t​×ES​q^​est,t−1∗​
   
  最终参与位姿估计运算的地磁向量$^E\bm {\hat b}_t$Eb^t​为：
   
  $^E\bm {\hat b}_t=[0,\sqrt{h_x^2+h_y^2},0,h_z]$Eb^t​=[0,hx2​+hy2​​,0,hz​]
   
  这样的校准计算可以将外界的电磁干扰影响约束在对传感器指向的方向估计中。而且不需要预先给定地磁方向。
• 陀螺仪零点漂移校正：
  $_E^S\bm{\dot{\hat q}}_\epsilon$ES​q^​˙​ϵ​代表了陀螺仪在各个轴的角误差，使用前文提到的$_E^S\bm{\dot q}_{\omega,t}=\frac{1}{2}_E^S\bm{\hat q}_{est,t-1}×^S\bm{\omega}_t$ES​q˙​ω,t​=21​ES​q^​est,t−1​×Sωt​公式的逆运算可以得出陀螺仪当前度数的可能偏差：
   
  $^S\bm {\omega}_{\epsilon,t}=2_E^S\bm{\hat q}_{est,t-1}^*×_E^S\bm{\dot{\hat q}}_{\epsilon,t}$Sωϵ,t​=2ES​q^​est,t−1∗​×ES​q^​˙​ϵ,t​
   
  那么陀螺仪随着时间的漂移量为：
   
  $^S\bm {\omega}_{b,t}=\zeta \sum{^S\bm {\omega}_{\epsilon,t}}\Delta t$Sωb,t​=ζ∑Sωϵ,t​Δt
   
  最后得到校正后的陀螺仪读数为：
   
  $^S\bm {\omega}_{c,t}=^S\bm {\omega}_{t}-^S\bm {\omega}_{b,t}$Sωc,t​=Sωt​−Sωb,t​
   
• 算法的整体流程图：
  
• 论文结尾附有该算法在C中的实现，详情请查阅论文。