MSCKF原理解析(参照论文)

参考论文地址 https://www-users.cs.umn.edu/~stergios/papers/ICRA07-MSCKF.pdf

算法特点

• 基于EKF框架（个人觉得叫ESKF更妥，文章中进行更新的状态都是误差状态）
• 提出了一种不需要将3D特征点包含到状态向量中的测量模型，使得计算规模小，复杂度低．

算法组成

• 状态的定义
• IMU预测(包括状态预测、协方差的更新)
• 状态增广
• 测量模型的构建
• EKF更新

下面将详细描述上述每个过程

状态定义

• IMU状态:
  $$X_{IMU}={\left[\begin{array}{ccccc}{}_G^I{\bar{q}}^T& {\mathbf{b}}_g^T& {}^G{\mathbf{v}}_I^T& {\mathbf{b}}_a^T& {}^G{\mathbf{p}}_I^T\end{array}\right]}^T$$
  式中$\mathbf{p},\mathbf{v},\bar{\mathbf{q}}$分别为IMU相对于世界坐标系的位置、速度和归一化四元素姿态，${\mathbf{b}}_g$,${\mathbf{b}}_a$分别为IMU陀螺仪和加速度计偏置．
• IMU误差状态:
  $${\tilde{X}}_{IMU}={\left[\begin{array}{ccccc}\delta {\mathbf{\theta }}_I^T& \widetilde{{\mathbf{b}}_g^T}& \widetilde{{}^G{\mathbf{v}}_I^T}& \widetilde{{\mathbf{b}}_a^T}& \widetilde{{}^G{\mathbf{p}}_I^T}\end{array}\right]}^T$$
  对于位置、速度和偏置状态，其误差状态直接通过标准加性误差$\tilde{x}=x-\hat{x}$定义．而对于姿态状态，其满足的关系是$\bar{q}=\delta q\otimes \hat{\bar{q}}$，其中$\delta q\simeq {\left[\begin{array}{cc}\frac{1}{2}\delta \mathbf{\theta }& 1\end{array}\right]}^T$．因此这里的姿态误差状态可以用$\delta q$或者$\delta \mathbf{\theta }$来表示，但$\delta q$是超参数表示，不利于后期的优化过程，所以使用$\delta \mathbf{\theta }$这种最小表达形式．
• EKF状态向量:
  $$\widehat{{\mathbf{X}}_k}=\left[\begin{array}{cccccc}{\hat{\mathbf{X}}}_{IM{U}_K}^T& {}_G^{C_1}{\hat{\bar{q}}}^T& {}^G{\hat{\mathbf{p}}}_{C_1}^T& ...& {}_G^{C_N}{\hat{\bar{q}}}^T& {}^G{\hat{\mathbf{p}}}_{C_N}^T\end{array}\right]$$
  该算法框架中长期维护的状态向量包括k时刻(当前时刻)IMU状态${\hat{\mathbf{X}}}_{IM{U}_K}$和N个相机状态$\left[\begin{array}{cc}{}_G^C{\hat{\bar{q}}}^T& {}^G{\hat{\mathbf{p}}}_C^T\end{array}\right]$
• EKF误差状态向量:
  $$\widetilde{{\mathbf{X}}_k}=\left[\begin{array}{cccccc}{\tilde{\mathbf{X}}}_{IM{U}_K}^T& \delta {\mathbf{\theta }}_{C_1}^T& {}^G{\tilde{\mathbf{p}}}_{C_1}^T& ...& \delta {\mathbf{\theta }}_{C_N}^T& {}^G{\tilde{\mathbf{p}}}_{C_N}^T\end{array}\right]$$
  该算法框架中利用的误差状态向量包括k时刻IMU误差状态${\tilde{\mathbf{X}}}_{IM{U}_K}$和N个相机误差状态$\left[\begin{array}{cc}\delta {\mathbf{\theta }}_C^T& {}^G{\tilde{\mathbf{p}}}_C^T\end{array}\right]$

IMU预测(当接收到IMU信息时)

1. 状态预测

• 真实状态系统模型
  $$\begin{array}{c}{}_G^I\dot{\bar{q}}(t)=\frac{1}{2}\Omega (\mathbf{\omega }(t)){}_G^I\bar{q}(t)\\ \dot{{\mathbf{b}}_g}(t)={\mathbf{n}}_{wg}(t)\\ {}^G\dot{{\mathbf{v}}_I}(t)={}^G\mathbf{a}(t)\\ \dot{{\mathbf{b}}_a}(t)={\mathbf{n}}_{wa}(t)\\ {}^G\dot{{\mathbf{p}}_I}(t)={}^G{\mathbf{v}}_I^T\end{array}$$
  式中$\mathbf{\omega }(t)=\left[\begin{array}{ccc}{\omega }_x& {\omega }_y& {\omega }_z\end{array}\right]$为实际角速度值，${}^G\mathbf{a}(t)$为实际加速度值，$\Omega (\mathbf{\omega })=\left[\begin{array}{cc}-\lfloor \mathbf{\omega }\times \rfloor & \mathbf{\omega }\\ -{\mathbf{\omega }}^T& 0\end{array}\right]$，$\lfloor \mathbf{\omega }\times \rfloor =\left[\begin{array}{ccc}0& -{\omega }_z& {\omega }_y\\ {\omega }_z& 0& -{\omega }_x\\ -{\omega }_y& {\omega }_x& 0\end{array}\right]$
• 名义状态系统模型
  $$\begin{array}{c}{}_G^I\dot{\hat{\bar{q}}}=\frac{1}{2}\Omega (\hat{\mathbf{\omega }}){}_G^I\hat{\bar{q}}\\ \dot{\widehat{{\mathbf{b}}_g}}=0_{3\times 1}\\ {}^G\dot{\widehat{{\mathbf{v}}_I}}={\mathbf{C}}_{\hat{q}}^T\hat{\mathbf{a}}+{}^G\mathbf{g}\\ {\dot{\hat{\mathbf{b}}}}_a=0_{3\times 1}\\ {}^G{\dot{\hat{\mathbf{p}}}}_I={}^G{\hat{\mathbf{v}}}_I\end{array}$$
  相比与论文中的模型，该模型忽略了地球自转影响，对该算法影响不大．式中${\mathbf{C}}_{\hat{q}}=\mathbf{C}({}_G^I\hat{\bar{q}})$为通过姿态估计值转换成的旋转矩阵，式中$\hat{\mathbf{a}}={\mathbf{a}}_m-{\hat{\mathbf{b}}}_a$，$\hat{\mathbf{\omega }}={\mathbf{\omega }}_m-{\hat{\mathbf{b}}}_g-{\mathbf{C}}_{\hat{q}}{\mathbf{\omega }}_G$，$({\mathbf{a}}_m,{\mathbf{\omega }}_m)$分别为加速度计和陀螺仪测量值．
• 误差状态系统模型(即EKF的预测模型)
  根据误差状态＝真实状态－名义状态定义，可得
  $${\dot{\tilde{\mathbf{X}}}}_{IMU}=\mathbf{F}{\tilde{\mathbf{X}}}_{IMU}+\mathbf{G}{\mathbf{n}}_{IMU}$$
  其中:
  $\mathbf{F}=\left[\begin{array}{ccccc}-\lfloor \hat{\mathbf{\omega }}\times \rfloor & -{\mathbf{I}}_3& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ -{\mathbf{C}}_{\hat{q}}^T\lfloor \hat{\mathbf{a}}\times \rfloor & 0_{3\times 3}& 0_{3\times 3}& -{\mathbf{C}}_{\hat{q}}^T& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& {\mathbf{I}}_3& 0_{3\times 3}& 0_{3\times 3}\end{array}\right],\mathbf{G}=\left[\begin{array}{cccc}-{\mathbf{I}}_3& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& {\mathbf{I}}_3& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& -{\mathbf{C}}_{\hat{q}}^T& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& {\mathbf{I}}_3\\ 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\end{array}\right]$
  $${\mathbf{n}}_{IMU}=\left[\begin{array}{cccc}{\mathbf{n}}_g^T& {\mathbf{n}}_{wg}^T& {\mathbf{n}}_a^T& {\mathbf{n}}_{wa}^T\end{array}\right]为IMU误差模型参数$$
  详细推导可查看我博客Quaternion kinematics for ESKF－预测中"连续时间内的系统运动方程"部分内容．

2. 协方差预测

• k时刻协方差
  $${\mathbf{P}}_{k|k}=\left[\begin{array}{cc}{\mathbf{P}}_{I{I}_{k|k}}& {\mathbf{P}}_{I{C}_{k|k}}\\ {\mathbf{P}}_{I{C}_{k|k}}^T& {\mathbf{P}}_{C{C}_{k|k}}\end{array}\right]$$
• 预测的k+1时刻协方差
  $${\mathbf{P}}_{k+1|k}=\left[\begin{array}{cc}{\mathbf{P}}_{I{I}_{k+1|k}}& \mathbf{\Phi }(t_k+T,t_k){\mathbf{P}}_{I{C}_{k|k}}\\ {\mathbf{P}}_{I{C}_{k|k}}^T\mathbf{\Phi }(t_k+T,t_k{)}^T& {\mathbf{P}}_{C{C}_{k|k}}\end{array}\right]$$
  其中${\mathbf{P}}_{I{I}_{k+1|k}}$通过李雅普诺夫方程计算:
  $${\dot{\mathbf{P}}}_{II}={\mathbf{F}\mathbf{P}}_{II}+{\mathbf{P}}_{II}{\mathbf{F}}^T+{\mathbf{G}\mathbf{Q}}_{IMU}{\mathbf{G}}^T$$
  $\mathbf{\Phi }(t_k+T,t_k)$通过下面方程进行数值积分可得(初始值$\mathbf{\Phi }(t_k,t_k)={\mathbf{I}}_{15\times 15}$):
  $$\dot{\mathbf{\Phi }}(t_k+\tau ,t_k)=\mathbf{F}\mathbf{\Phi }(t_k+\tau ,t_k),\tau \in [0,T]$$

状态增广(当接收到新图片信息时)

• 状态增广
  当接收到新图片时，需要将新的相机位姿添加到滤波的状态向量中，添加的状态按如下方法获取:
  $${}_G^C\hat{\bar{q}}={}_I^C\bar{q}\otimes {}_G^I\hat{\bar{q}},{}^G{\hat{\mathbf{p}}}_C={}^G{\hat{\mathbf{p}}}_I+{\mathbf{C}}_{\hat{q}}^T{}^I{\mathbf{p}}_C$$
• 协方差增广
  由于状态的增广，使得协方差矩阵应对应增广，方法如下:
  $${\mathbf{P}}_{k|k}=\left[\begin{array}{c}{\mathbf{I}}_{6N+15}\\ \mathbf{J}\end{array}\right]{\mathbf{P}}_{k|k}{\left[\begin{array}{c}{\mathbf{I}}_{6N+15}\\ \mathbf{J}\end{array}\right]}^T$$
  其中:
  $$\mathbf{J}=\left[\begin{array}{cccc}\mathbf{C}({}_I^C\bar{\mathbf{q}})& 0_{3\times 9}& 0_{3\times 3}& 0_{3\times 6N}\\ \lfloor {\mathbf{C}}_{\hat{q}}^T{}^I{\mathbf{p}}_C\times \rfloor & 0_{3\times 9}& {\mathbf{I}}_3& 0_{3\times 6N}\end{array}\right]$$

测量模型的构建

这部分是该算法的核心，保证了算法的高效运行．传统的基于EKF框架的vio算法都会将特征点作为状态向量的一部分在整个算法运行过程中进行更新矫正，这就导致了状态向量的维度急速提高，要知道在进行滤波算法的过程中，是需要计算观测方程的jacobi矩阵的，维度一旦太大，必然导致计算复杂性快速提高，导致算法的效率降低，不利于大规模场景的使用．
该算法不维护特征点信息，而是当一个特征点未被继续观察到时，通过观察到该特征点的相机状态利用最小二乘方法获取特征点的位置，然后利用该位置信息和观测到的特征信息来构建观测模型．

• 实际测量模型
  假设单个特征$f_j$位置${\mathbf{p}}_{f_j}$被$S_j$个相机状态$({}_G^{C_i}\bar{q},{}^G{\mathbf{p}}_{C_i}),i\in S_j$观测到，则可建立如下实际测量模型:
  $${\mathbf{z}}_i^{(j)}=\frac{1}{{}^{C_i}{Z}_j}\left[\begin{array}{c}{}^{C_i}{X}_j\\ {}^{C_i}{Y}_j\end{array}\right]+{\mathbf{n}}_i^{(j)},i\in S_j$$
  其中$[{}^{C_i}{X}_j,{}^{C_i}{Y}_j,{}^{C_i}{Z}_j]$为特征点$f_j$在相机状态$C_j$坐标系下的坐标，可通过下式计算:
  $${}^{C_i}{\mathbf{p}}_{f_j}=[{}^{C_i}{X}_j,{}^{C_i}{Y}_j,{}^{C_i}{Z}_j{]}^T=\mathbf{C}({}_G^{C_j}\bar{q})({}^G{\mathbf{p}}_{f_j}-{}^G{\mathbf{p}}_{c_i})$$
  另外${\mathbf{n}}_i^{(j)}$是图像的噪声向量，满足${\mathbf{n}}_i^{(j)}\sim N(0,{\sigma }_{im}^2\mathbf{I})$，为高斯白噪声．
• 名义测量模型
  $${\hat{\mathbf{z}}}_i^{(j)}=\frac{1}{{}^{C_i}\widehat{Z_j}}\left[\begin{array}{c}{}^{C_i}\widehat{X_j}\\ {}^{C_i}\widehat{Y_j}\end{array}\right],{}^{C_i}{\hat{\mathbf{p}}}_{f_j}=[{}^{C_i}\widehat{X_j},{}^{C_i}\widehat{Y_j},{}^{C_i}\widehat{Z_j}{]}^T=\mathbf{C}({}_G^{C_j}\hat{\bar{q}})({}^G{\hat{\mathbf{p}}}_{f_j}-{}^G{\hat{\mathbf{p}}}_{c_i})$$
  其中的名义状态${}_G^{C_j}\hat{\bar{q}},{}^G{\hat{\mathbf{p}}}_{f_j},{}^G{\hat{\mathbf{p}}}_{c_i}$均为状态增广过程添加的相机状态．
• 观测方程
  建立观测方程如下:$${\mathbf{r}}_i^{(j)}={\mathbf{z}}_i^{(j)}-{\hat{\mathbf{z}}}_i^{(j)}$$
  由于该方程为非线性方程，需要进行泰勒展开是线性方程，才能进行滤波，展开后方程:
  $${\mathbf{r}}_i^{(j)}\simeq {\mathbf{H}}_{{\mathbf{X}}_i}^j\tilde{X}+{\mathbf{H}}_{f_j}^j{}^G{\tilde{\mathbf{p}}}_{f_j}+{\mathbf{n}}_i^j$$

其中${\mathbf{H}}_{{\mathbf{X}}_i}^j,{\mathbf{H}}_{f_j}^j$分别为测量值${\mathbf{z}}_i^{(j)}$相对于状态和特征位置的jacobia矩阵．
将所有观测到特征$f_j$的相机状态构建的观测方程累加起来如下:
$${\mathbf{r}}^{(j)}\simeq {\mathbf{H}}_{\mathbf{X}}^j\tilde{X}+{\mathbf{H}}_f^j{}^G{\tilde{\mathbf{p}}}_{f_j}+{\mathbf{n}}^j$$
从上面的展开后的观测方程可以看出，观测过程不仅与误差状态有关还与特征位置有关，因此不能直接将上式作为最终的观测方程，否则达不到该算法的特点:不将特征加入状态向量中．为了达到算法效果，将${\mathbf{r}}^{(j)}$投影到${\mathbf{H}}_f^{(j)}$的左零空间中，可获得如下改进后的观测方程:
$${\mathbf{r}}_o^{(j)}={\mathbf{A}}^T({\mathbf{z}}^{(j)}-{\hat{\mathbf{z}}}^{(j)})\simeq {\mathbf{A}}^T{\mathbf{H}}_{\mathbf{X}}^{(j)}\tilde{\mathbf{X}}+{\mathbf{A}}^T{\mathbf{n}}^{(j)}={\mathbf{H}}_o^{(j)}{\tilde{X}}^{(j)}+{\mathbf{n}}_o^{(j)}$$
其中A为半正交矩阵，它的列构成了$H_f$的左零空间，${\mathbf{r}}_o$为投影后的残差项．

EKF更新

当某个特征点不被继续观测到或者相机状态太多需要消除老状态时，进行EKF的更新．由于${\mathbf{H}}_o$矩阵比较大，为了降低计算的复杂性，对${\mathbf{H}}_o$矩阵进行QR分解:
$${\mathbf{H}}_o=\left[\begin{array}{cc}{\mathbf{Q}}_1& {\mathbf{Q}}_2\end{array}\right]\left[\begin{array}{c}{\mathbf{T}}_H\\ 0\end{array}\right]$$
其中${\mathbf{Q}}_1,{\mathbf{Q}}_2,{\mathbf{T}}_H$分别为酉矩阵、酉矩阵和上三角矩阵．
于是原观测方程可表达为:
$${\mathbf{r}}_o=\left[\begin{array}{cc}{\mathbf{Q}}_1& {\mathbf{Q}}_2\end{array}\right]\left[\begin{array}{c}{\mathbf{T}}_H\\ 0\end{array}\right]\tilde{X}+{\mathbf{n}}_o\Rightarrow \left[\begin{array}{c}{\mathbf{Q}}_1^T{\mathbf{r}}_o\\ {\mathbf{Q}}_2^T{\mathbf{r}}_o\end{array}\right]=\left[\begin{array}{c}{\mathbf{T}}_H\\ 0\end{array}\right]\tilde{X}+\left[\begin{array}{c}{\mathbf{Q}}_1^T{\mathbf{n}}_o\\ {\mathbf{Q}}_2^T{\mathbf{n}}_o\end{array}\right]$$
由于${\mathbf{Q}}_2^T{\mathbf{r}}_o$只是噪声项，可以忽略，因此可以得到最终的观测方程如下:
$${\mathbf{r}}_n={\mathbf{Q}}_1^T{\mathbf{r}}_o={\mathbf{T}}_H\tilde{X}+{\mathbf{n}}_n$$
卡尔曼增益:
$$\mathbf{K}=\mathbf{P}{\mathbf{T}}_H^T({\mathbf{T}}_H\mathbf{P}{\mathbf{T}}_H^T+{\mathbf{R}}_n{)}^{-1}$$
误差状态向量增量:
$$\Delta \mathbf{X}=\mathbf{K}{\mathbf{r}}_n$$
协方差矩阵更新:
$${\mathbf{P}}_{k+1|k+1}=({\mathbf{I}}_{\xi }-{\mathbf{K}\mathbf{T}}_H){\mathbf{P}}_{k+1|k}({\mathbf{I}}_{\xi }-{\mathbf{K}\mathbf{T}}_H{)}^T+{\mathbf{K}\mathbf{R}}_n{\mathbf{K}}^T$$