一种深空探测巡视器惯性/视觉组合导航方法

摘要

本发明涉及一种深空探测巡视器惯性/视觉组合导航方法，根据行星固连坐标系下捷联惯导的力学编排，建立巡视器惯性/视觉组合导航的状态模型；利用双目视觉相机对巡视器周围环境进行拍摄，获得立体图像序列，通过图像特征提取算法得到特征点的像素坐标，建立以特征点像素坐标为量测向量的量测模型；并使用Unscented卡尔曼滤波估计深空探测巡视器的位置、速度和姿态。本发明能够有效修正惯导的位置误差，提高导航精度，非常适用于深空探测巡视器的自主导航，本发明属于航天导航技术领域，不仅可以为深空探测巡视器提供高精度导航参数，而且可为其自主导航系统设计提供参考。

主权项

1.一种深空探测巡视器的惯性/视觉组合导航方法，其特征在于：首先根据行星固连坐标系下的捷联惯导力学编排方程，建立惯性/视觉组合导航的状态模型，再利用双目视觉相机对周围环境进行拍摄，获得立体图像序列，通过图像特征提取算法得到特征点的像素坐标，建立以特征点像素坐标为量测量的量测模型，并使用Unscented卡尔曼滤波估计深空探测巡视器的位置、速度和姿态；具体包括以下步骤：①建立深空探测巡视器惯性/视觉组合导航系统状态模型；选取行星固连坐标系为全局参考坐标系，可得组合导航的状态模型：$\left\{\begin{array}{c}{\stackrel{·}{q}}_p^b=\frac{1}{2}\Omega ({\omega }_{\mathrm{ib}}^b-b_g-C_p^b{\omega }_{\mathrm{ip}}^p)q_p^b+w_g\\ {\stackrel{·}{p}}_b^p=v_b^p\\ {\stackrel{·}{v}}_b^p=C_b^p(f_b-b_a)-2{\omega }_{\mathrm{ip}}^p\times v_b^p+g^p+w_a\\ {\stackrel{·}{b}}_g=0\\ {\stackrel{·}{b}}_a=0\end{array}\right.---\left(1\right)$式中，为从行星固连坐标系到巡视器本体坐标系旋转的单位四元数，为4×1维的向量，和分别为巡视器在行星固连坐标系下的位置和速度，均为3×1维的向量，b_g和b_a分别为巡视器本体坐标系下陀螺仪和加速度计的常值偏置，均为3×1维的向量，为巡视器本体坐标系相对惯性坐标系的旋转角速率在巡视器本体坐标系中的表示，由陀螺仪直接测得，f_b为比力信息，由加速度计直接测得，为行星自转角速率在行星固连坐标系下的矢量表示，为巡视器本体坐标系相对行星固连坐标系的方向余弦矩阵，g^p为重力矢量在行星固连坐标系下的表示，Ω(·)为四元数矩阵，w_g和w_a分别为陀螺仪和加速度计的随机误差，均为零均值的高斯白噪声；式(1)中变量都是与t有关的函数，可简写为：$\stackrel{·}{X}\left(t\right)=f(X\left(t\right),t)+w\left(t\right)---\left(2\right)$状态变量为f(X(t),t)为系统非线性连续状态转移函数，状态噪声为w(t)＝[w_g,w_a]^T；②建立基于特征点像素坐标为量测量的量测模型；巡视器本体坐标系与IMU坐标系重合，摄像机坐标系和IMU坐标系之间的转换关系已知，t_k时刻和t_k+1时刻摄像机坐标系之间的旋转矩阵和平移矢量为其中为3×3维的矩阵，为3×1维的向量，计算公式为：$R_k^{k+1}=C_b^c{C}_{p,k+1}^b{\left(C_b^c{C}_{p,k}^b\right)}^T$                                           (3)$T_k^{k+1}=-C_b^c{C}_{p,k+1}^b(p_{b,k+1}^p+{C_{p,k+1}^b}^T{p}_c^b)+C_b^c{C}_{p,k+1}^b(p_{b,k}^p+{C_{p,k}^b}^T{p}_c^b)$式中，k表示时间序号，k＝1,2,…，为摄像机相对巡视器本体坐标系的安装矩阵，为摄像机在巡视器本体坐标系中的位置，为3×1维的向量，由标定可知，和分别为t_k时刻和t_k+1时刻巡视器在行星固连坐标系中的位置，均为3×1维的向量，和分别为t_k时刻和t_k+1时刻巡视器本体坐标系相对行星固连坐标系的姿态矩阵，均为3×3维，和根据t_k时刻组合导航的滤波结果获得，和由t_k+1时刻惯导的状态更新结果获得；跟踪得到的第i个特征点在相邻时刻摄像机坐标系下的三维坐标关系式为：${\hat{p}}_{i,k+1}^{c_l}=R_k^{k+1}{p}_{i,k}^{c_l}+T_k^{k+1}---\left(4\right)$式中，k＝1,2,…，为第i个特征点在t_k时刻左摄像机坐标系下的三维坐标，为第i个特征点在t_k+1时刻左摄像机坐标系下三维坐标的估计值，分别为从t_k时刻到t_k+1时刻摄像机坐标系之间的旋转矩阵和平移矢量，由式(3)获得；根据上述求出的和摄像机成像的针孔模型，可以求出t_k+1时刻特征点在左摄像机成像平面内的像素坐标预测值，对所有M个特征点的表达式为：${\hat{Z}}_{i,k+1}^l=\left[\begin{array}{c}{\hat{u}}_{i,k+1}^l\\ {\hat{v}}_{i,k+1}^l\end{array}\right]=\left[\begin{array}{cc}\frac{1}{\mathrm{dX}}& 0\\ 0& \frac{1}{\mathrm{dY}}\end{array}\right]\left[\begin{array}{c}\frac{f_{c_l}{\hat{p}}_{i,k+1,x}^{c_l}}{{\hat{p}}_{i,k+1,z}^{c_l}}\\ \frac{f_{c_l}{\hat{p}}_{i,k+1,y}^{c_l}}{{\hat{p}}_{i,k+1,z}^{c_l}}\end{array}\right]+\left[\begin{array}{c}{u}_{l0}\\ v_{l0}\end{array}\right]+{\tilde{v}}_{i,k+1}^l,i=\mathrm{1,2},...,M---\left(5\right)$式中，为第i个特征点在左摄像机中对应的量测信息，为第i个特征点t_k+1时刻在左摄像机成像平面中的像素坐标预估值，dX,dY分别为每个像素在左摄像机成像平面横轴和纵轴方向的长和宽，(u_l0,v_l0)为左摄像机光轴与成像平面交点的像素坐标，为左摄像机焦距，为第i个特征点在t_k+1时刻左摄像机坐标系下的三维坐标预估值，M为t_k时刻到t_k+1时刻之间跟踪到的特征点个数，为左图像中特征点提取误差，为2×1维列向量；根据平行双目立体视觉系统的外部参数，可以得到t_k+1时刻第i个特征点在右摄像机坐标系下三维坐标的估计值表达式为：${\hat{p}}_{i,k+1}^{c_r}=R({\hat{p}}_{i,k+1}^{c_l}-T)---\left(6\right)$式中，R,T分别为右摄像机坐标系相对于左摄像机坐标系的旋转矩阵和平移矢量；根据上述求出的和摄像机成像的针孔模型，可以求出t_k+1时刻特征点在右摄像机成像平面内的像素坐标预测值，对所有M个特征点的表达式为：${\hat{Z}}_{i,k+1}^r=\left[\begin{array}{c}{\hat{u}}_{i,k+1}^r\\ {\hat{v}}_{i,k+1}^r\end{array}\right]=\left[\begin{array}{cc}\frac{1}{\mathrm{dX}}& 0\\ 0& \frac{1}{\mathrm{dY}}\end{array}\right]\left[\begin{array}{c}\frac{f_{c_r}{\hat{p}}_{i,k+1,x}^{c_r}}{{\hat{p}}_{i,k+1,z}^{c_r}}\\ \frac{f_{c_r}{\hat{p}}_{i,k+1,y}^{c_r}}{{\hat{p}}_{i,k+1,z}^{c_r}}\end{array}\right]+\left[\begin{array}{c}{u}_{r0}\\ v_{r0}\end{array}\right]+{\tilde{v}}_{i,k+1}^r,$i＝1,2,…,M             (7)式中，为第i个特征点在右摄像机中对应的量测信息，为第i个特征点t_k+1时刻在右摄像机成像平面中的像素坐标预估值，dX,dY分别为每个像素在右摄像机成像平面横轴和纵轴方向的长和宽，其值与左摄像机相同，(u_r0，v_r0)为右摄像机光轴与成像平面交点的像素坐标，为右摄像机焦距，为第i个特征点在t_k+1时刻右摄像机坐标系下的三维坐标预估值，M为t_k时刻到t_k+1时刻之间跟踪到的特征点个数，为右图像中特征点提取误差，为2×1维列向量；将全部特征点在左、右摄像机成像平面中的像素坐标构成一个列向量$Z={[{\hat{Z}}_{1,k+1}^l,{\hat{Z}}_{1,k+1}^r,{\hat{Z}}_{2,k+1}^l,{\hat{Z}}_{2,k+1}^r,...,{\hat{Z}}_{M,k+1}^l,{\hat{Z}}_{M,k+1}^r]}^T,$即系统量测量，共4M维，设量测噪声$v={[{\tilde{v}}_1,{\tilde{v}}_2,{\tilde{v}}_3,...,{\tilde{v}}_{4M}]}^T,$${\tilde{v}}_1,{\tilde{v}}_2,{\tilde{v}}_3,...,{\tilde{v}}_{4M}$分别为$u_1^l,v_1^l,u_1^r,v_1^r,u_2^l,v_2^l,u_2^r,v_2^r,...,u_M^l,v_M^l,u_M^r{,v}_M^r$的观测误差，由于各变量都是与时间t有关的变量，则可建立基于像素坐标量测信息的量测模型表达式为：Z(t)＝h(X(t),t)+v(t)                    (8)③对步骤①中式(2)所示的状态模型及步骤②中式(8)所示的量测模型进行离散化：X(k+1)＝F(X(k),k)+w(k)                   (9)Z(k)＝H(X(k),k)+v(k)                        (10)式中，k＝1,2,…，F(X(k),k)为f(X(t),t)离散后的非线性状态转移函数，H(X(k),k)为h(X(t),t)离散后的非线性量测函数，w(k)、v(k)分别为离散后的系统噪声和量测噪声，它们互不相关；④利用Unscented卡尔曼滤波算法，结合步骤①和步骤②所述的状态模型和量测模型进行滤波，利用双目视觉相机拍摄立体图像序列，再进行特征提取获得特征点的像素坐标信息，通过量测量与量测模型得到的量测预测值相减得到系统量测残差，用系统这一残差校正状态一步预测的误差；根据状态向量可得相应的Unscented采样点，利用系统状态模型，对采样点进行一步预测，并得出与上一时刻滤波得到的迭代状态值之间的协方差阵，以消除状态模型中模型误差的影响；⑤输出导航信息；将步骤④中获得的状态量的估计值和状态估计方差用于下一时刻的滤波，并输出状态估计值和状态估计方差，同时，将步骤④中获得的行星固连坐标系中的位置、速度和姿态信息通过矩阵变换转换到以东北天坐标系表示的导航坐标系中并输出。