CN104764449A - 一种基于星历修正的捕获段深空探测器自主天文导航方法 - Google Patents

Abstract

本发明公开了一种基于星历修正的捕获段深空探测器自主天文导航方法，首先建立目标天体星历误差状态模型和量测模型，并根据预测与实际天体图像的位置差获取星历误差量测量；其次将星历误差状态模型和轨道动力学模型联立作为天文导航系统状态模型，并将星历误差量测模型和星光角距量测模型作为天文导航系统的量测模型，采用Unscented卡尔曼滤波方法，估计探测器位置、速度和目标天体星历误差，并将所估计的星历误差反馈至状态模型中，修正目标天体星历数据，获得自主校正星历误差后的相对于目标天体和相对于日心的探测器位置和速度。本发明属于航天导航技术领域，可在线估计天体星历误差，修正导航系统的模型误差，适用于探测器捕获段。

Description

一种基于星历修正的捕获段深空探测器自主天文导航方法

技术领域

本发明涉及在深空探测器处于捕获段时，基于目标天体图像和预测天体图像进行星历修正的自主天文导航方法，是一种非常适用于深空探测器捕获段的导航方法。

背景技术

捕获段是指从深空探测器进入目标影响球开始到点火制动的全过程，处于该阶段的深空探测器飞行速度快，飞行弧段短，控制精度要求高，且机会唯一。深空探测器减速制动的入轨点距离目标行星表面非常近(近地点)，被捕获是整个深空探测任务的一个关键时间节点，捕获阶段相对目标天体的相对导航精度和相对于太阳或者地球的绝对导航精度对后续探测任务有直接影响。然而深空探测器在捕获段航行速度快、空间电离环境未知、大气环境复杂，这些因素都对深空探测器的入轨精度有很大的影响，也制约着深空探测器捕获后绕飞、着陆等阶段的导航精度，当入轨精度不能达到指标要求时，甚至无法完成科学探测任务，导致整个任务的失败。

目标天体星历数据是影响深空探测器捕获段导航性能的主要因素之一。目标天体星历数据是描述目标天体位置、速度等特征的一类天体数据库，由长期天文观测拟合而成的。若一段时间没有天文观测信息进行修正，目标天体星历数据误差会随时间的递推而增加。太阳系内行星(除地球外)目前的星历误差约为200m～100km。

处于捕获段的深空探测器距离地球较远，飞行速度较快，飞行距离远，时间长，环境未知性较强，传统上依靠地面测控的深空探测器导航与控制方法在实时性、成本和资源上受到种种限制，存在很多不足，很难满足深空探测一些特殊任务对高精度导航与控制的需要，因而对自主导航的要求更为迫切。

传统基于天体星历修正的导航方法主要有两种，一种方法为将无线电导航信息与天文导航信息组合进行几何解析，求解目标天体星历及星历误差，此后将无线电导航导航结果转换至目标天体为中心的惯性坐标系中，获得较高精度的导航信息。但是这种方法由于没有实时修正状态模型的模型误差，仍然无法消除目标天体星历误差对探测器状态模型精度的影响，导航精度仍较低。另一种方法(专利号：201310006575.X一种基于星历修正的深空探测器天文/无线电组合导航方法)为利用天文和无线电两个导航系统所估计的探测器位置，结合目标天体、探测器、地球的几何关系得到目标天体星历误差量测量，并获得目标天体星历误差量测量，利用卡尔曼滤波方法估计目标天体星历误差，并反馈至导航系统模型中，再进行天文/无线电测控导航信息融合。但是这两种方法都必须依赖无线电测控信息的协助，才可确定目标天体星历误差，仍然无法自主实时修正目标天体星历误差及其对状态模型误差和导航结果的影响。

发明内容

本发明要解决的技术问题是：克服目标天体星历误差对相对地球导航精度的影响，弥补现有天文/无线电测控导航方法难以自主消除目标天体星历误差对探测器状态模型的影响这一不足，合理有效利用天文导航所提供的导航信息，为深空探测器捕获段提供一种高精度自主修正目标天体星历的天文导航方法。

本发明解决其技术问题所采用的技术方案为：首先建立目标天体星历误差状态模型和量测模型，并根据预测天体图像与实际天体图像的位置差获取目标天体星历误差量测量；其次将目标天体星历误差状态模型和基于太阳和八大行星引力的轨道动力学模型联立作为天文导航系统状态模型，并将目标天体星历误差量测模型和天文星光角距量测模型作为天文导航系统的量测模型，根据探测器的状态模型和量测模型都呈现非线性特性，且系统噪声为非高斯噪声的特点，采用Unscented卡尔曼滤波方法，估计探测器位置、速度和目标天体星历误差，并将所估计的目标天体星历误差反馈至状态模型中，修正目标天体星历数据，提高状态模型的模型精度，并最终提高相对于目标天体和相对于日心的探测器导航信息的估计精度。

具体包括以下步骤：

1.建立目标天体星历误差状态模型

根据目标天体星历误差在捕获段时间段内慢变这一特点，建立目标天体星历误差状态模型为：

$\left\{\begin{array}{c}{\stackrel{\·}{b}}_x=w_{b_x}\\ {\stackrel{\·}{b}}_y=w_{b_y}\\ {\stackrel{\·}{b}}_z=w_{b_z}\end{array}\right.---\left(1\right)$

式中，b_x,b_y,b_z为目标天体星历三轴位置误差，为目标天体星历三轴位置误差的微分，为目标天体星历误差状态模型误差；

2.建立目标天体星历误差量测模型

B_m＝[b_x b_y b_z]^T+w_m    (2)

式中，w_m为目标天体星历误差的测量噪声；B_m为目标天体星历误差量测量。

3.量测量目标天体星历误差的获取

图2给出了天文导航系统所用目标天体敏感器的成像原理图，目标天体卫星敏感器成像过程与之类似。目标天体敏感器主要由光学透镜和二维成像面阵组成，在敏感器测量坐标系O′X_cY_cZ_c中沿目标天体到探测器的矢量方向目标天体反射太阳光线射向目标天体敏感器，此时，目标天体在敏感器测量坐标系中的坐标为(x_c,y_c,z_c)；目标天体敏感器的光学透镜以焦距f将目标天体的光线折射后成像在二维成像面阵上，二维成像面阵将照在每个成像单元上的图像亮度信号储存。

(1)目标天体矢量方向的获取

根据目标天体敏感器的成像原理，目标天体矢量方向的获取步骤如下所示：

A.天体图像的处理

由于目标天体在二维成像面阵上的图像并不是一个点，而是一个圆，通过质心识别等图像处理技术确定目标天体图像在二维成像平面坐标系OX_2dY_2d的质心(x_2d,y_2d)，这一中心可以用像元像线坐标系O_plX_plY_pl中的像元像线(p,l)表示。

B.质心坐标从像元像线坐标系转换至二维成像平面坐标系的转换

根据像元像线坐标系与二维成像平面坐标系之间的关系，将目标天体质心坐标从像元像线坐标系转换至二维成像平面坐标系：

$\left[\begin{array}{c}{x}_{2d}\\ y_{2d}\end{array}\right]=K^{-1}(\left[\begin{array}{c}p\\ l\end{array}\right]-\left[\begin{array}{c}{p}_0\\ l_0\end{array}\right])---\left(3\right)$

式中，(x_2d,y_2d)为目标天体在二维成像平面OX_2dY_2d中的坐标，p和l分别为目标天体在敏感器二维成像平面的像元和像线， $K=\left[\begin{array}{cc}{K}_x& K_{\mathrm{yx}}\\ K_{\mathrm{xy}}& K_y\end{array}\right]$ 为由毫米转为像素的敏感器转换矩阵,p₀和l₀分别为目标天体敏感器中心在像元像线坐标系OX_plY_pl中的像元和像线。

C.质心坐标从二维成像平面坐标系转换至敏感器测量坐标系的转换

根据透镜成像原理，将目标天体质心坐标从二维成像平面坐标系转换至敏感器测量坐标系：

${\stackrel{\→}{l}}_{\mathrm{pm}}^s=\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]=\frac{1}{\sqrt{x_{2d}^2+y_{2d}^2+f^2}}\left[\begin{array}{c}{x}_{2d}\\ y_{2d}\\ -f\end{array}\right]---\left(4\right)$

式中，f为目标天体敏感器的焦距，为在目标天体敏感器测量坐标系中从目标天体至探测器的单位矢量。同理可获得目标天体图像中第i颗背景恒星在敏感器测量坐标系中的单位矢量i＝1,2,3。

因此，可以根据预测图像获得预测目标天体质心所对应的单位矢量实际图像中目标天体质心所对应的单位矢量

(2)目标天体图像半径的获取

由于目标天体在二维成像面阵上的图像并不是一个点，而是一个圆，利用边缘检测算法，提取火星图像轮廓，计算目标天体图像轮廓距离目标天体质心的距离，获得在二维成像平面坐标系中目标天体实际图像的半径和在二维成像平面坐标系中目标天体预测图像的半径

(3)探测器到目标天体质心距离的解算

探测器到实际目标天体质心之间的距离可以表示为：

$r_t^a=l_{\mathrm{pt}}^a{R}_{t\arg \mathrm{et}}/R_{\mathrm{image}}^a---\left(5\right)$

式中，R_target为目标天体半径，可由地面天文数据库获得，f为目标天体敏感器的焦距，为实际目标天体质心在二维成像平面坐标系中的坐标；

探测器到星历所预测的目标天体质心之间的距离可以表示为：

$r_t^p=l_{\mathrm{pt}}^p{R}_{t\arg \mathrm{et}}/R_{\mathrm{image}}^p---\left(6\right)$

式中，为根据星历所预测的目标天体质心在二维成像平面坐标系中的坐标。

(4)目标天体星历误差的解算

在敏感器测量坐标系中，目标天体星历误差矢量可以根据探测器到实际目标天体质心的矢量和探测器到星历所预测目标天体质心的矢量计算得出：

${\mathrm{\δl}}_{2d}^c=r_t^p{l}_t^p-r_t^a{l}_t^a---\left(7\right)$

由此可得，在惯性坐标系中目标天体星历误差可以表示为：

$B=A_{\mathrm{ib}}{A}_{\mathrm{bc}}{\mathrm{\δl}}_{2d}^c---\left(8\right)$

式中，A_ib为探测器本体系至惯性系的姿态转移矩阵，A_bc为探测器敏感器坐标系至本体系的转移矩阵；

4.建立深空探测器天文导航系统的状态模型

在目标天体为中心的惯性坐标系中，建立深空探测器基于太阳和八大行星引力轨道动力学模型，并与步骤1所建立的目标天体星历误差模型联立，构成天文导航系统的状态模型；

$\left\{\begin{array}{c}\stackrel{\·}{x}=v_x+w_x\\ \stackrel{\·}{y}=v_y+w_y\\ \stackrel{\·}{z}=v_z+w_z\\ {\stackrel{\·}{v}}_x=-{\mathrm{\μ}}_m\frac{x}{r_{\mathrm{pm}}^3}-{\mathrm{\μ}}_s[\frac{x-x_s}{r_{\mathrm{ps}}^3}+\frac{x_s}{r_{\mathrm{ms}}^3}]-\underset{i_c}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i_c}[\frac{x-x_{i_c}}{r_{{\mathrm{pi}}_c}^3}+\frac{x_{i_c}}{r_{{\mathrm{mi}}_c}^3}]+w_{v_x}\\ {\stackrel{\·}{v}}_y=-{\mathrm{\μ}}_m\frac{y}{r_{\mathrm{pm}}^3}-{\mathrm{\μ}}_s[\frac{y-y_s}{r_{\mathrm{ps}}^3}+\frac{y_s}{r_{\mathrm{ms}}^3}]-\underset{i_c}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i_c}[\frac{y-y_{i_c}}{r_{{\mathrm{pi}}_c}^3}+\frac{y_{i_c}}{r_{{\mathrm{mi}}_c}^3}]+w_{v_y}\\ {\stackrel{\·}{v}}_z=-{\mathrm{\μ}}_m\frac{z}{r_{\mathrm{pm}}^3}-{\mathrm{\μ}}_s[\frac{z-z_s}{r_{\mathrm{ps}}^3}+\frac{z_s}{r_{\mathrm{ms}}^3}]-\underset{i_c}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i_c}[\frac{z-z_{i_c}}{r_{{\mathrm{pi}}_c}^3}+\frac{z_{i_c}}{r_{{\mathrm{mi}}_c}^3}]+w_{v_z}\\ {\stackrel{\·}{b}}_x=w_{b_x}\\ {\stackrel{\·}{b}}_y=w_{b_y}\\ {\stackrel{\·}{b}}_z=w_{b_z}\end{array}\right.---\left(9\right)$

式中，x,y,z为探测器在目标天体质心惯性坐标系中三轴位置，v_x,v_y,v_z为探测器在目标天体质心惯性坐标系中三轴速度，b_x,b_y,b_z为目标天体星历三轴位置误差，为探测器在目标天体质心惯性坐标系中三轴位置的微分，为探测器在目标天体质心惯性坐标系中三轴速度的微分，为目标天体星历三轴位置误差的微分，μ_s、μ_m和μ_ic分别为太阳、目标天体和第i_c颗行星的引力常数；r_ps为日心到探测器的距离；r_pm为目标天体质心到探测器的距离；r_ms为日心到目标天体质心的距离；为第i_c颗行星到探测器的距离；为第i_c颗行星质心到火心的距离；(x_s,y_s,z_s)，分别为太阳、第i_c颗行星在目标天体质心惯性坐标系中的三轴位置坐标，可根据时间由行星星历表获得，w_x,w_y,w_z,分别为状态模型中探测器三轴位置和速度状态模型误差；i_c表示太阳和八大行星中从内至外的第i_c颗行星，如目标天体为火星，则i_c＝1,2,3...,N(i_c≠4)，N＝8，由于i_c＝4表示第4颗行星(火星)，因此不包含在求和表达式中。

式中的各变量都是与时间t有关的变量，可简写为：

$\stackrel{\·}{X}\left(t\right)=f(X\left(t\right),t)+w\left(t\right)---\left(10\right)$

式中，X(t)＝[x,y,z,v_x,v_y,v_z,b_x,b_y,b_z]^T为天文导航系统状态模型的状态向量，为X(t)的微分，f(X(t),t)为天文导航系统状态模型的系统非线性连续状态转移函数， $w\left(t\right)={\left[\begin{array}{ccccccccc}{w}_x& w_y& w_z& w_{v_x}& w_{v_y}& w_{v_z}& w_{b_x}& w_{b_y}& w_{b_z}\end{array}\right]}^T$ 为天文导航系统状态模型误差。

5.建立天文导航系统量测模型

天文导航系统量测模型选取星光角距模型和步骤2中所建立的星历误差模型作为量测模型，其中目标天体和两个卫星与三颗背景恒星之间的星光角距θ_1i、θ_2i和θ_3i(i＝1,2,3)表达式为：

$\left\{\begin{array}{c}{\mathrm{\θ}}_{1i}=\arccos (-{\stackrel{\→}{l}}_{p1}^I\·{\stackrel{\→}{s}}_{1i}^I)+v_{1i}\\ {\mathrm{\θ}}_{2i}=\arccos (-{\stackrel{\→}{l}}_{p2}^I\·{\stackrel{\→}{s}}_{2i}^I)+v_{2i}\\ {\mathrm{\θ}}_{3i}=\arccos (-{\stackrel{\→}{l}}_{p3}^I\·{\stackrel{\→}{s}}_{3i}^I)+v_{3i}\end{array}\right.---\left(11\right)$

式中，为在惯性坐标系中目标天体到探测器的单位矢量，为目标天体图像中第i个恒星在惯性坐标系中的单位矢量；为在惯性坐标系中目标天体的第一个卫星到探测器的单位矢量，为目标天体第一个卫星图像中第i个恒星在惯性坐标系中的单位矢量；为在惯性坐标系中第二个目标天体的卫星到探测器的单位矢量，为目标天体第二个卫星图像中第i个恒星在惯性坐标系中的单位矢量。

与目标天体星历误差量测模型式联立，得天文导航系统量测模型为：

$\left\{\begin{array}{c}{\mathrm{\θ}}_{1i}=\arccos (-{\stackrel{\→}{l}}_{p1}^I\·{\stackrel{\→}{s}}_{1i}^I)+v_{1i}\\ {\mathrm{\θ}}_{2i}=\arccos (-{\stackrel{\→}{l}}_{p2}^I\·{\stackrel{\→}{s}}_{2i}^I)+v_{2i}\\ {\mathrm{\θ}}_{3i}=\arccos (-{\stackrel{\→}{l}}_{p3}^I\·{\stackrel{\→}{s}}_{3i}^I)+v_{3i}\\ B_m={\left[\begin{array}{ccc}{b}_x& b_y& b_z\end{array}\right]}^T+v_m\end{array}\right.---\left(12\right)$

设天文导航系统量测量Z＝[θ₁₁,θ₁₂,θ₁₃,θ₂₁,θ₂₂,θ₂₃,θ₃₁,θ₃₂,θ₃₃,b_x,b_y,b_z]^T，天文导航系统量测噪声 $v={[v_{{\mathrm{\θ}}_{11}},v_{{\mathrm{\θ}}_{12}},v_{{\mathrm{\θ}}_{13}},v_{{\mathrm{\θ}}_{21}},v_{{\mathrm{\θ}}_{22}},v_{{\mathrm{\θ}}_{23}},v_{{\mathrm{\θ}}_{31}},v_{{\mathrm{\θ}}_{32}},v_{{\mathrm{\θ}}_{33}},v_{b_x},v_{b_y},v_{b_z}]}^T,$ 分别为测量θ₁₁,θ₁₂,θ₁₃,θ₂₁,θ₂₂,θ₂₃,θ₃₁,θ₃₂,θ₃₃,b_x,b_y,b_z的观测误差，由于各变量都是与时间t有关的变量，则可根据式(12)建立天文导航系统量测模型的表达式为：

Z(t)＝h[X(t),t]+v(t)    (13)

式中，h[X(t),t]为天文导航系统非线性连续量测函数。

6.对步骤4和步骤5中的状态模型和量测模型进行离散化

X(k)＝F(X(k-1),k-1)+W(k-1)    (14)

Z(k)＝H(X(k),k)+V(k)    (15)

式中，k＝1,2,…，F(X(k-1),k-1)为f(X(t),t)离散后从第k-1时刻到第k时刻的非线性状态转移函数，H(X(k),k)为h(X(t),t)离散后第k时刻的非线性量测函数，W(k)和V(k)为w(t)和v(t)离散后第k时刻的等效噪声，且W(k)和V(k)互不相关。

7.天文导航量测量的获取及处理

天文导航系统的量测量Z＝[θ₁₁,θ₁₂,θ₁₃,θ₂₁,θ₂₂,θ₂₃,θ₃₁,θ₃₂,θ₃₃,b_x,b_y,b_z]^T可分为两部分：(1)星光角距；(2)目标天体星历误差。其中量测量星光角距的获取可根据在火星敏感器测量坐标系中目标天体至探测器的单位矢量与目标天体图像中第i颗背景恒星的单位矢量计算两个矢量之间的星光角距θ_1i：

${\mathrm{\θ}}_{1i}=\arccos [(-{\stackrel{\→}{l}}_{\mathrm{pm}}^S)\·{\stackrel{\→}{s}}_{1i}^S]---\left(16\right)$

同理可获得目标天体卫星与其背景恒星的星光角距θ_2i,θ_3i。而量测量目标天体星历误差的获取可根据步骤3中的式(3)～式(8)获得；

8.对天文导航系统进行Unscented卡尔曼滤波

根据天文导航系统状态模型、天文导航系统量测模型、天文导航敏感器获得的量测量，进行天文导航系统Unscented卡尔曼滤波，获得在目标天体惯性坐标系中表示深空探测器位置、速度和目标天体星历误差的估计状态向量 ${\hat{x}}_k={\left[\begin{array}{ccccccccc}{\hat{x}}_k& {\hat{y}}_k& {\hat{z}}_k& {\hat{v}}_{\mathrm{xk}}& {\hat{v}}_{\mathrm{yk}}& {\hat{v}}_{\mathrm{zk}}& {\hat{b}}_{\mathrm{xk}}& {\hat{b}}_{\mathrm{yk}}& {\hat{b}}_{\mathrm{zk}}\end{array}\right]}^T$ 和估计均方误差阵P_k，其中为第k时刻探测器在目标天体质心惯性坐标系中的三轴位置估计、三轴速度估计、目标天体星历三轴估计误差。

9.天文导航系统对目标天体星历误差进行反馈校正

将第8步中获得的目标天体星历误差和目标天体星历估计均方误差阵反馈回深空探测器的状态模型中，并重新确定状态模型的模型误差协方差阵，最后将校正后的模型误差协方差阵输入至天文导航系统Unscented卡尔曼滤波中，修正下一时刻的导航结果。

10.计算日心惯性坐标系中的探测器位置、速度估计状态向量

将第8步获得的第k时刻探测器在目标天体质心惯性坐标系中的位置、速度估计状态向量转换至日心惯性坐标系中，得到日心惯性坐标系中的估计状态向量

${\hat{X}}^s\left(k\right)={\hat{X}}^m\left(k\right)+X_{t\arg \mathrm{et}}^s\left(k\right)+\mathrm{\δ}{\hat{X}}_{t\arg \mathrm{et}}^s\left(k\right)---\left(17\right)$

式中，为探测器第k时刻在火心惯性坐标系中的位置、速度估计状态向量 ${\hat{X}}^m\left(k\right)={\left[\begin{array}{cccccc}{\hat{x}}_k& {\hat{y}}_k& {\hat{z}}_k& {\hat{v}}_{\mathrm{xk}}& {\hat{v}}_{\mathrm{yk}}& {\hat{v}}_{\mathrm{zk}}\end{array}\right]}^T;$ 为目标天体在日心惯性坐标系中第k时刻的位置、速度，可由目标天体星历数据库获得；为目标天体第k时刻的位置、速度误差，且 $\mathrm{\δ}{\hat{X}}_{t\arg \mathrm{et}}^s\left(k\right)={\left[\begin{array}{cccccc}{\hat{b}}_{\mathrm{xk}}& {\hat{b}}_{\mathrm{yk}}& {\hat{b}}_{\mathrm{zk}}& 0& 0& 0\end{array}\right]}^T.$

最终将第8步和第10步获得的第k时刻在目标天体质心惯性坐标系中和在日心惯性坐标系中的估计状态向量和估计均方误差阵P_k输出，估计状态向量分别表示在目标天体质心惯性坐标系中和日心惯性坐标系中探测器的位置、速度信息，输出的估计均方误差阵P_k表示了滤波估计的性能，并将这些导航信息分别返回天文导航系中，用于k+1时刻的位置、速度导航信息，k＝1,2,...。

本发明的原理是：由于天文导航敏感器可以直接测量目标天体图像，而利用天体星历数据也可预测天文导航系统所获得的目标天体图像，但预测天体图像与实际天体图像在敏感器成像平面的位置不同，位置的变化反映了目标天体星历误差的大小，本发明利用该信息自主修正目标天体星历误差。首先结合捕获段持续时间短、目标天体星历误差在捕获段变化小的特点，建立目标天体星历误差状态模型和量测模型，并根据预测天体图像与实际天体图像的位置差获取目标天体星历误差量测量；其次将目标天体星历误差状态模型和基于太阳和八大行星引力的轨道动力学模型联立作为天文导航系统状态模型，并将目标天体星历误差量测模型和天文星光角距量测模型作为天文导航系统的量测模型，根据探测器的状态模型和量测模型都呈现非线性特性，且系统噪声为非高斯噪声的特点，采用Unscented卡尔曼滤波方法，估计探测器位置、速度和目标天体星历误差，并将所估计的目标天体星历误差反馈至状态模型中，修正目标天体星历数据，提高状态模型的模型精度；提高对相对于目标天体和相对于日心的探测器导航信息的估计精度。

本发明与现有技术相比的优点在于：一方面实现了不再需要无线电测控信息就可以自主获得目标天体的星历误差，实现了对目标天体星历误差的高精度估计；另一方面天文导航系统可自主获得相对目标天体和相对日心的探测器高精度导航信息，并自主修正了导航系统的状态模型，减小了目标天体星历误差对状态模型精度的影响。

附图说明

图1为本发明基于自主星历修正的天文导航流程图。

图2为本发明天文导航系统所用天文导航敏感器的成像原理图。

具体实施方式

如图1所示，前述技术方案中所涉及的目标天体可以为火星、金星、木星、土星等太阳系内的天体，以下以火星作为实施例，说明本发明的具体实施过程：

1.建立火星星历误差状态模型

根据火星星历误差在捕获段时间段内慢变这一特点，建立火星星历误差状态模型为：

$\left\{\begin{array}{c}{\stackrel{\·}{b}}_x=w_{b_x}\\ {\stackrel{\·}{b}}_y=w_{b_y}\\ {\stackrel{\·}{b}}_z=w_{b_z}\end{array}\right.---\left(1\right)$

式中，b_x,b_y,b_z为火星星历三轴位置误差，为火星星历三轴位置误差的微分，为火星星历误差状态模型误差；

2.建立火星星历误差量测模型

B_m＝[b_x b_y b_z]^T+w_m    (2)

式中，w_m为火星星历误差的测量噪声；B_m为火星星历误差量测量。

3.量测量火星星历误差的获取

图2给出了天文导航系统所用火星敏感器的成像原理图，火卫一敏感器、火卫二敏感器成像过程与之类似。火星敏感器主要由光学透镜和二维成像面阵组成，在敏感器测量坐标系O′X_cY_cZ_c中沿火星到探测器的矢量方向火星反射太阳光线射向火星敏感器，此时，火星在火星敏感器测量坐标系中的坐标为(x_c,y_c,z_c)；火星敏感器的光学透镜以焦距f将火星的光线折射后成像在二维成像面阵上，二维成像面阵将照在每个成像单元上的图像亮度信号储存。

(1)火星矢量方向的获取

根据敏感器的成像原理，天体矢量方向的获取步骤如下所示：

A.天体图像的处理

由于火星在二维成像面阵上的图像并不是一个点，而是一个圆，通过质心识别等图像处理技术确定火星图像在二维成像平面坐标系OX_2dY_2d的质心(x_2d,y_2d)，这一中心可以用像元像线坐标系O_plX_plY_pl中的像元像线(p,l)表示。

B.质心坐标从像元像线坐标系转换至二维成像平面坐标系的转换

根据像元像线坐标系与二维成像平面坐标系之间的关系，将火星质心坐标从像元像线坐标系转换至二维成像平面坐标系：

$\left[\begin{array}{c}{x}_{2d}\\ y_{2d}\end{array}\right]=K^{-1}(\left[\begin{array}{c}p\\ l\end{array}\right]-\left[\begin{array}{c}{p}_0\\ l_0\end{array}\right])---\left(3\right)$

式中，(x_2d,y_2d)为火星在二维成像平面OX_2dY_2d中的坐标，p和l分别为火星在火星敏感器二维成像平面的像元和像线， $K=\left[\begin{array}{cc}{K}_x& K_{\mathrm{yx}}\\ K_{\mathrm{xy}}& K_y\end{array}\right]$ 为由毫米转为像素的敏感器转换矩阵,p₀和l₀分别为火星敏感器中心在像元像线坐标系OX_plY_pl中的像元和像线。

C.质心坐标从二维成像平面坐标系转换至敏感器测量坐标系的转换

根据透镜成像原理，将火星质心坐标从二维成像平面坐标系转换至敏感器测量坐标系：

${\stackrel{\→}{l}}_{\mathrm{pm}}^s=\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]=\frac{1}{\sqrt{x_{2d}^2+y_{2d}^2+f^2}}\left[\begin{array}{c}{x}_{2d}\\ y_{2d}\\ -f\end{array}\right]---\left(4\right)$

式中，f为火星敏感器的焦距，为在火星敏感器测量坐标系中从火星至探测器的单位矢量。同理可获得火星图像中第i颗背景恒星在火星敏感器测量坐标系中的单位矢量i＝1,2,3。

因此，可以根据预测火星图像获得预测火星图像中火星质心所对应的单位矢量实际火星图像中火星质心所对应的单位矢量

(2)火星图像半径的获取

由于火星在二维成像面阵上的图像并不是一个点，而是一个圆，利用边缘检测算法，提取火星图像轮廓，计算火星图像轮廓距离火星质心的距离，获得在二维成像平面坐标系中火星实际图像的半径和在二维成像平面坐标系中火星预测图像的半径

(3)火星探测器到火星质心距离的解算

火星探测器到实际火星质心之间的距离可以表示为：

$r_t^a=l_{\mathrm{pt}}^a{R}_{t\arg \mathrm{et}}/R_{\mathrm{image}}^a---\left(5\right)$

式中，R_target＝3397km为火星半径，f为火星敏感器的焦距，为实际火星质心在二维成像平面坐标系中的坐标；

火星探测器到火星星历所预测的火星质心之间的距离可以表示为：

$r_t^p=l_{\mathrm{pt}}^p{R}_{t\arg \mathrm{et}}/R_{\mathrm{image}}^p---\left(6\right)$

式中，为根据火星星历所预测的火星质心在二维成像平面坐标系中的坐标。

(4)火星星历误差的解算

在敏感器测量坐标系中，火星星历误差矢量可以根据探测器到实际火星质心的矢量和探测器到火星星历所预测火星质心的矢量计算得出：

${\mathrm{\δl}}_{2d}^c=r_t^p{l}_t^p-r_t^a{l}_t^a---\left(7\right)$

由此可得，在惯性坐标系中火星星历误差可以表示为：

$B=A_{\mathrm{ib}}{A}_{\mathrm{bc}}{\mathrm{\δl}}_{2d}^c---\left(8\right)$

式中，A_ib为探测器本体系至惯性系的姿态转移矩阵，A_bc为探测器敏感器坐标系至本体系的转移矩阵；

4.建立深空探测器天文导航系统的状态模型

首先初始化位置、速度、火星星历误差，设X＝[x y z v_x v_y v_z b_x b_y b_z]^T为在火心惯性坐标系中深空探测器天文导航系统的状态向量，x,y,z,v_x,v_y,v_z分别为探测器在火心惯性坐标系中三轴的位置和速度，b_x,b_y,b_z为火星星历三轴位置误差，上述变量都是与t有关的函数，根据探测器的轨道设计和火星星历数据库轨道，选取探测器的位置、速度初值为x₀,y₀,z₀,火星星历三轴位置误差初值为则深空探测器天文导航系统的状态向量初值为 $X\left(0\right)=\left[\begin{array}{ccccccccc}{x}_0& y_0& z_0& v_{x_0}& v_{y_0}& v_{z_0}& b_{x_0}& b_{y_0}& b_{z_0}];\end{array}\right]$

考虑太阳和火星、地球等八大行星对探测器的引力作用，选取火心惯性坐标系，联立步骤1获得的火星星历误差状态模型，可得深空探测器在火心惯性坐标系中状态模型：

$\left\{\begin{array}{c}\stackrel{\·}{x}=v_x+w_x\\ \stackrel{\·}{y}=v_y+w_y\\ \stackrel{\·}{z}=v_z+w_z\\ {\stackrel{\·}{v}}_x=-{\mathrm{\μ}}_m\frac{x}{r_{\mathrm{pm}}^3}-{\mathrm{\μ}}_s[\frac{x-x_s}{r_{\mathrm{ps}}^3}+\frac{x_s}{r_{\mathrm{ms}}^3}]-\underset{i_c}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i_c}[\frac{x-x_{i_c}}{r_{{\mathrm{pi}}_c}^3}+\frac{x_{i_c}}{r_{{\mathrm{mi}}_c}^3}]+w_{v_x}\\ {\stackrel{\·}{v}}_y=-{\mathrm{\μ}}_m\frac{y}{r_{\mathrm{pm}}^3}-{\mathrm{\μ}}_s[\frac{y-y_s}{r_{\mathrm{ps}}^3}+\frac{y_s}{r_{\mathrm{ms}}^3}]-\underset{i_c}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i_c}[\frac{y-y_{i_c}}{r_{{\mathrm{pi}}_c}^3}+\frac{y_{i_c}}{r_{{\mathrm{mi}}_c}^3}]+w_{v_y}\\ {\stackrel{\·}{v}}_z=-{\mathrm{\μ}}_m\frac{z}{r_{\mathrm{pm}}^3}-{\mathrm{\μ}}_s[\frac{z-z_s}{r_{\mathrm{ps}}^3}+\frac{z_s}{r_{\mathrm{ms}}^3}]-\underset{i_c}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i_c}[\frac{z-z_{i_c}}{r_{{\mathrm{pi}}_c}^3}+\frac{z_{i_c}}{r_{{\mathrm{mi}}_c}^3}]+w_{v_z}\\ {\stackrel{\·}{b}}_x=w_{b_x}\\ {\stackrel{\·}{b}}_y=w_{b_y}\\ {\stackrel{\·}{b}}_z=w_{b_z}\end{array}\right.---\left(9\right)$

式中，x,y,z为探测器在火心惯性坐标系中三轴位置，v_x,v_y,v_z为探测器在火心惯性坐标系中三轴速度，b_x,b_y,b_z为火星星历三轴位置误差，为探测器在火心惯性坐标系中三轴位置的微分，为探测器在火心惯性坐标系中三轴速度的微分，为火星星历三轴位置误差的微分，μ_s、μ_m和分别为太阳、火星和第i_c颗行星的引力常数；r_ps为日心到探测器的距离；r_pm为火星到探测器的距离；r_ms为日心到火心的距离；为第i_c颗行星到探测器的距离；为第i_c颗行星质心到火心的距离；(x_s,y_s,z_s)，分别为太阳、第i_c颗行星在火心惯性坐标系中的三轴位置坐标，可根据时间由行星星历表获得，w_x，w_y，w_z，分别为状态模型中探测器三轴位置和速度的状态模型误差；i_c表示太阳和八大行星中从内至外的第i_c颗行星，i_c＝1,2,3...,N(i_c≠4)，N＝8，由于i_c＝4表示第4颗行星(火星)，因此不包含在求和表达式中。

式(9)中的各变量都是与时间t有关的变量，可简写为：

$\stackrel{\·}{X}\left(t\right)=f(X\left(t\right),t)+w\left(t\right)---\left(10\right)$

式中，X(t)＝[x,y,z,v_x,v_y,v_z,b_x,b_y,b_z]^T为天文导航系统状态模型的状态向量，为X(t)的微分，f(X(t),t)为天文导航系统状态模型的系统非线性连续状态转移函数， $w\left(t\right)={\left[\begin{array}{ccccccccc}{w}_x& w_y& w_z& w_{v_x}& w_{v_y}& w_{v_z}& w_{b_x}& w_{b_y}& w_{b_z}\end{array}\right]}^T$ 为天文导航系统状态模型误差。

5.建立天文导航系统量测模型；

天文导航系统量测模型选取火星、火卫一、火卫二与各自背景恒星的星光角距模型和星历误差模型作为量测模型。

火星与第i颗背景恒星之间星光角距θ_mi的大小在不同坐标系中是相同的，因此其表达式为：

${\mathrm{\θ}}_{\mathrm{mi}}=\arccos (-{\stackrel{\→}{l}}_{\mathrm{pm}}^S\·{\stackrel{\→}{s}}_{1i}^S)+v_{\mathrm{mi}}=\mathrm{arcoos}(-{\stackrel{\→}{l}}_{\mathrm{pm}}^I\·{\stackrel{\→}{s}}_{1i}^I)+-v_{\mathrm{mi}}--\left(11\right)$

式中，v_mi为θ_mi的测量噪声，为在火星敏感器测量坐标系中从火星至探测器的单位矢量，为在惯性坐标系中火星到探测器的单位矢量，可表示为：

${\stackrel{\→}{l}}_{\mathrm{pm}}^I=\frac{1}{\sqrt{x^2+y^2+z^2}}\left[\begin{array}{c}x\\ y\\ z\end{array}\right]---\left(12\right)$

式中，(x,y,z)为探测器在火心惯性坐标系中三轴位置坐标，为火星图像中第i颗背景恒星在火星敏感器测量坐标系中的单位矢量，为在火星图像中第i颗背景恒星星光在惯性坐标系中的单位矢量，i＝1,2,3，可由恒星星历数据库直接得到；

同理可得火卫一和火卫二与其第i颗背景恒星之间星光角距θ_pi和θ_di的表达式为：

${\mathrm{\θ}}_{\mathrm{pi}}=\arccos (-{\stackrel{\→}{l}}_{\mathrm{pp}}^I\·{\stackrel{\→}{s}}_{2i}^I)+v_{\mathrm{pi}}---\left(13\right)$

${\mathrm{\θ}}_{\mathrm{di}}=\arccos (-{\stackrel{\→}{l}}_{\mathrm{pd}}^I\·{\stackrel{\→}{s}}_{3i}^I)+v_{\mathrm{di}}---\left(14\right)$

式中，v_pi,v_di分别为θ_pi,θ_di的测量噪声。

与火星星历误差量测模型式联立，得天文导航系统量测模型为：

$\left\{\begin{array}{c}{\mathrm{\θ}}_{\mathrm{mi}}=\arccos (-{\stackrel{\→}{l}}_{\mathrm{pm}}^I\·{\stackrel{\→}{s}}_{1i}^I)+v_{\mathrm{mi}}\\ {\mathrm{\θ}}_{\mathrm{pi}}=\arccos (-{\stackrel{\→}{l}}_{\mathrm{pp}}^I\·{\stackrel{\→}{s}}_{2i}^I)+v_{\mathrm{pi}}\\ {\mathrm{\θ}}_{\mathrm{di}}=\arccos (-{\stackrel{\→}{l}}_{\mathrm{pd}}^I\·{\stackrel{\→}{s}}_{3i}^I)+v_{\mathrm{di}}\\ B_m={\left[\begin{array}{ccc}{b}_x& b_y& b_z\end{array}\right]}^T+v_m\end{array}\right.---\left(15\right)$

设天文导航系统量测量Z＝[θ_m1,θ_m2,θ_m3,θ_p1,θ_p2,θ_p3,θ_d1,θ_d2,θ_d3,b_x,b_y,b_z]^T，天文导航系统量测噪声 $v={[v_{{\mathrm{\θ}}_{m1}},v_{{\mathrm{\θ}}_{m2}},v_{{\mathrm{\θ}}_{m3}},v_{{\mathrm{\θ}}_{p1}},v_{{\mathrm{\θ}}_{p2}},v_{{\mathrm{\θ}}_{p3}},v_{{\mathrm{\θ}}_{d1}},v_{{\mathrm{\θ}}_{d2}},v_{{\mathrm{\θ}}_{d3}},v_{b_x},v_{b_y},v_{b_z}]}^T,$ 分别为测量θ_m1,θ_m2,θ_m3,θ_p1,θ_p2,θ_p3,θ_d1,θ_d2,θ_d3,b_x,b_y,b_z的观测误差，由于各变量都是与时间t有关的变量，则可根据式(15)建立天文导航系统量测模型的表达式为：

Z(t)＝h[X(t),t]+v(t)    (16)

式中，h[X(t),t]为天文导航系统非线性连续量测函数。

6.对步骤4和步骤5中的状态模型和量测模型进行离散化

X(k)＝F(X(k-1),k-1)+W(k-1)    (17)

Z(k)＝H(X(k),k)+V(k)    (18)

式中，k＝1,2,…，F(X(k-1),k-1)为f(X(t),t)离散后从第k-1时刻到第k时刻的非线性状态转移函数，H(X(k),k)为h(X(t),t)离散后第^k时刻的非线性量测函数，W(k)和V(k)为w(t)和v(t)离散后第k时刻的等效噪声，且W(k)和V(k)互不相关。

7.天文导航系统量测量的获取及处理

天文导航系统的量测量Z＝[θ_m1,θ_m2,θ_m3,θ_p1,θ_p2,θ_p3,θ_d1,θ_d2,θ_d3,b_x,b_y,b_z]^T可分为两部分：(1)星光角距；(2)火星星历误差。其中量测量星光角距的获取可根据在火星敏感器测量坐标系中火星至探测器的单位矢量与火星图像中第i颗背景恒星的单位矢量计算两个矢量之间的星光角距θ_mi

${\mathrm{\θ}}_{\mathrm{mi}}=\arccos [(-{\stackrel{\→}{l}}_{\mathrm{pm}}^S)\·{\stackrel{\→}{s}}_{1i}^S]---\left(19\right)$

同理可获得火卫一与其背景恒星、火卫二与其背景恒星之间的星光角距θ_pi,θ_di。

而量测量火星星历误差的获取可根据步骤3中的式(3)～式(8)获得；

8.对天文导航系统进行Unscented卡尔曼滤波

根据状态模型式(17)、天文导航系统量测模型式(18)、由天文导航敏感器获得的量测量式(8)和式(19)，进行天文导航系统Unscented卡尔曼滤波。具体步骤如下：

A.初始化

${\hat{x}}_0=E\left[x_0\right],P_0=E\left[(x_0-{\hat{x}}_0){(x_0-{\hat{x}}_0)}^T\right]---\left(20\right)$

式中，为第0时刻(初始时刻)在火心惯性坐标系中探测器的三轴位置和速度估计值，x₀为第0时刻(初始时刻)在火心惯性坐标系中探测器的三轴位置和速度真实值，P₀为状态向量的初始均方误差阵。

B.计算采样点

在天文导航系统第k-1时刻状态向量附近选取一系列样本点，这些样本点的均值和均方误差阵分别为和P_k-1。设状态向量为6×1维，那么13个天文导航系统的样本点χ_0,k,...,χ_12,k及其权重W₀…W₁₂分别如下：

${\mathrm{\χ}}_{0,k-1}={\hat{x}}_{k-1},$ W₀＝-1

${\mathrm{\χ}}_{j,k-1}={\hat{x}}_{k-1}+\sqrt{3}{\left(\sqrt{P_{k-1}}\right)}_j,W_j^{\′}=1/6---\left(21\right)$

${\mathrm{\χ}}_{j+6,k-1}={\hat{x}}_{k-1}-\sqrt{3}{\left(\sqrt{P_{k-1}}\right)}_j,$ W′_j+6＝1/6

式中,当P_k-1＝A^TA时，取A的第j行，当P_k-1＝AA^T时，取A的第j列，得第k-1时刻采样点χ_k-1的统一表达式为：

C.时间更新

天文导航系统状态向量的一步预测χ_k|k-1为：

χ_k|k-1＝F₁(χ_k-1,k-1)    (23)

天文导航系统所有采样点状态向量的一步预测加权后结果为：

${\hat{x}}_k^{-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j{\mathrm{\χ}}_{j,k|k-1}---\left(24\right)$

式中，W_j为第j个采样点的权值；

天文导航系统状态向量的估计均方误差阵一步预测为：

$P_k^{-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j[{\mathrm{\χ}}_{j,k|k-1}-{\hat{x}}_k^{-}]{[{\mathrm{\χ}}_{j,k|k-1}-{\hat{x}}_k^{-}]}^T+Q_k---\left(25\right)$

式中，Q′_k为k时刻天文导航系统的状态模型误差协方差阵；

天文导航系统采样点对应的量测估计向量Z_k|k-1

Z_k|k-1＝H₁(χ_k|k-1,k)    (26)

天文导航系统所有采样点量测估计加权向量

${\hat{z}}_k^{-}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j{Z}_{j,k|k-1}---\left(27\right)$

D.量测更新

天文导航系统量测均方误差阵为：

$P_{{\hat{z}}_k{\hat{z}}_k}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j[Z_{j,k|k-1}-{\hat{z}}_k^{-}]{[Z_{j,k|k-1}-{\hat{z}}_k^{-}]}^T+R_k---\left(28\right)$

式中，R_k为k时刻天文导航系统的量测噪声协方差阵；

天文导航系统状态向量量测量均方误差阵

$P_{{\hat{x}}_k{\hat{z}}_k}=\underset{j=0}{\overset{12}{\mathrm{\Σ}}}{W}_j[{\mathrm{\χ}}_{j,k|k-1}-{\hat{x}}_k^{-}]{[Z_{j,k|k-1}-{\hat{z}}_k^{-}]}^T---\left(29\right)$

天文导航系统滤波增益K_k为：

$K_k=P_{{\hat{x}}_k{\hat{z}}_k}{P}_{{\hat{z}}_k{\hat{z}}_k}^{-1}---\left(30\right)$

天文导航系统的估计状态向量和估计均方误差阵P_k为：

${\hat{x}}_k={\hat{x}}_k^{-}+K_k(Z_k-{\hat{z}}_k^{-})---\left(31\right)$

$P_k=P_k^{-}-K_k{P}_{{\hat{z}}_k{\hat{z}}_k}{K}_k^T---\left(32\right)$

式中，天文导航系统的估计状态向量 ${\hat{x}}_k={\left[\begin{array}{ccccccccc}{\hat{x}}_k& {\hat{y}}_k& {\hat{z}}_k& {\hat{v}}_{\mathrm{xk}}& {\hat{v}}_{\mathrm{yk}}& {\hat{v}}_{\mathrm{zk}}& {\hat{b}}_{\mathrm{xk}}& {\hat{b}}_{\mathrm{yk}}& {\hat{b}}_{\mathrm{zk}}\end{array}\right]}^T,$ 为第k时刻探测器在火心惯性坐标系中的三轴位置估计、三轴速度估计、火星星历三轴估计误差。

9.天文导航系统对火星星历误差进行反馈校正

将第8步的D中获得的火星星历误差和火星星历估计均方误差阵反馈回深空探测器的状态模型中，并重新确定状态模型的模型误差协方差阵Q_k，最后将校正后的模型误差协方差阵Q_k输入至天文导航系统Unscented卡尔曼滤波中，修正下一时刻的导航结果。

10.计算日心惯性坐标系中的探测器位置、速度估计状态向量

将式(31)获得的探测器第k时刻在火心惯性坐标系中的位置、速度估计状态向量转换至日心惯性坐标系中，得到日心惯性坐标系中的估计状态向量

${\hat{X}}^s\left(k\right)={\hat{X}}^m\left(k\right)+X_{\mathrm{tmars}}^s\left(k\right)+\mathrm{\δ}{\hat{X}}_{\mathrm{tmars}}^s\left(k\right)---\left(33\right)$

式中，为探测器第k时刻在火心惯性坐标系中的位置、速度估计状态向量 ${\hat{X}}^m\left(k\right)={\left[\begin{array}{cccccc}{\hat{x}}_k& {\hat{y}}_k& {\hat{z}}_k& {\hat{v}}_{\mathrm{xk}}& {\hat{v}}_{\mathrm{yk}}& {\hat{v}}_{\mathrm{zk}}\end{array}\right]}^T,$ 可由式(31)获得；为火星在日心惯性坐标系中第k时刻的位置、速度，可由火星星历数据库获得；为火星第k时刻的位置、速度误差，且 $\mathrm{\δ}{\hat{X}}_{\mathrm{mars}}^s\left(k\right)={\left[\begin{array}{cccccc}{\hat{b}}_{\mathrm{xk}}& {\hat{b}}_{\mathrm{yk}}& {\hat{b}}_{\mathrm{zk}}& 0& 0& 0\end{array}\right]}^T.$

最终将式(31)～式(33)获得的第k时刻在火心惯性坐标系中和在日心惯性坐标系中的估计状态向量和估计均方误差阵P_k输出，估计状态向量分别表示在火心惯性坐标系中和日心惯性坐标系中探测器的位置、速度信息，输出的估计均方误差阵P_k表示了滤波估计的性能，并将这些导航信息分别返回天文导航系中，用于k+1时刻的位置、速度导航信息，k＝1,2,...。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (1)

1.一种基于星历修正的捕获段深空探测器自主天文导航方法，其特征在于：建立目标天体星历误差状态模型和量测模型，并根据预测天体图像与实际天体图像的位置差获取目标天体星历误差量测量；其次将目标天体星历误差状态模型和基于太阳和八大行星引力的轨道动力学模型联立作为天文导航系统状态模型，并将目标天体星历误差量测模型和星光角距量测模型作为天文导航系统的量测模型，采用Unscented卡尔曼滤波方法，估计探测器位置、速度和目标天体星历误差，并将所估计的目标天体星历误差反馈至状态模型中，修正目标天体星历数据，提高状态模型的模型精度，获得自主校正星历误差后的相对于目标天体和相对于日心的探测器位置和速度；具体包括以下步骤：

①建立目标天体星历误差状态模型

$\left\{\begin{array}{c}{\stackrel{.}{b}}_x=w_{b_x}\\ {\stackrel{.}{b}}_y=w_{b_y}\\ {\stackrel{.}{b}}_z=w_{b_z}\end{array}\right.$

式中，b_x,b_y,b_z为目标天体星历三轴位置误差，为目标天体星历三轴位置误差的微分，w_bx,w_by,w_bz为目标天体星历误差状态模型误差；

②建立目标天体星历误差量测模型

B_m＝[b_x b_y b_z]^T+w_m

式中，w_m为目标天体星历误差的测量噪声；B_m为目标天体星历误差量测量；

③量测量目标天体星历误差的获取

A.目标天体矢量方向的获取

根据目标天体敏感器成像原理，获得预测目标天体质心所对应的单位矢量实际图像中目标天体质心所对应的单位矢量

B.目标天体图像半径的获取

由于目标天体在二维成像面阵上的图像并不是一个点，而是一个圆，利用边缘检测算法，提取火星图像轮廓，计算目标天体图像轮廓距离目标天体质心的距离，获得在二维成像平面坐标系中目标天体实际图像的半径和在二维成像平面坐标系中目标天体预测图像的半径

C.探测器到目标天体质心距离的解算

探测器到实际目标天体质心之间的距离可以表示为：

$r_t^a=l_{\mathrm{pt}}^a{R}_{t\arg \mathrm{et}}/R_{\mathrm{image}}^a$

式中，R_target为目标天体半径，可由地面天文数据库获得，f为目标天体敏感器的焦距，为实际目标天体质心在二维成像平面坐标系中的坐标；

探测器到星历所预测的目标天体质心之间的距离可以表示为：

$r_t^p=l_{\mathrm{pt}}^p{R}_{t\arg \mathrm{et}}/R_{\mathrm{image}}^p$

式中， 为根据星历所预测的目标天体质心在二维成像平面坐标系中的坐标；

D.目标天体星历误差的解算

在敏感器测量坐标系中，目标天体星历误差矢量可以根据探测器到实际目标天体质心的矢量和探测器到星历所预测目标天体质心的矢量计算得出：

${\mathrm{\δl}}_{2d}^c=r_t^p{l}_t^p-r_t^a{l}_t^a$

由此可得，在惯性坐标系中目标天体星历误差可以表示为：

$B=A_{\mathrm{ib}}{A}_{\mathrm{bc}}{\mathrm{\δl}}_{2d}^c$

式中，A_ib为探测器本体系至惯性系的姿态转移矩阵，A_bc为探测器敏感器坐标系至本体系的转移矩阵；

④建立深空探测器天文导航系统的星历误差扩维状态模型

在目标天体为中心的惯性坐标系中，建立深空探测器基于太阳和八大行星引力轨道动力学模型，，并与步骤①所建立的目标天体星历误差模型联立，构成天文导航系统的状态模型；

$\stackrel{.}{X}\left(t\right)=f(X\left(t\right),t)+w\left(t\right)$

式中，X(t)＝[x,y,z,v_x,v_y,v_z,b_x,b_y,b_z]^T为天文导航系统状态模型的状态向量，x,y,z,v_x,v_y,v_z分别为探测器在目标天体质心惯性坐标系中三轴的位置和速度，b_x,b_y,b_z为目标天体星历三轴位置误差，为X(t)的微分，f(X(t),t)为天文导航系统状态模型的系统非线性连续状态转移函数， $w\left(t\right)={\left[\begin{array}{ccccccccc}{w}_x& w_y& w_z& w_{v_x}& w_{v_y}& w_{v_z}& w_{b_x}& w_{b_y}& w_{b_z}\end{array}\right]}^T$ 为天文导航系统状态模型误差；

⑤建立天文导航系统量测模型

天文导航系统量测模型选取星光角距模型和步骤②中所建立的星历误差模型作为量测模型：

$\left\{\begin{array}{c}{\mathrm{\θ}}_{1i}=\arccos (-{\stackrel{\→}{l}}_{p1}^I\·{\stackrel{\→}{s}}_{1i}^I)+v_{1i}\\ {\mathrm{\θ}}_{2i}=\arccos (-{\stackrel{\→}{l}}_{p2}^I\·{\stackrel{\→}{s}}_{3i}^I)+v_{3i}\\ {\mathrm{\θ}}_{3i}=\arccos (-{\stackrel{\→}{l}}_{p3}^I\·{\stackrel{\→}{s}}_{3i}^I)+v_{3i}\\ B_m={\left[\begin{array}{ccc}{b}_x& b_y& b_z\end{array}\right]}^T+v_m\end{array}\right.$

式中，θ_1i、θ_2i和θ_3i(i＝1,2,3)为目标天体和两个卫星与三颗背景恒星之间的星光角距，为在惯性坐标系中目标天体到探测器的单位矢量，为目标天体图像中第i个恒星在惯性坐标系中的单位矢量；为在惯性坐标系中目标天体第一个卫星到探测器的单位矢量，为目标天体第一个卫星图像中第i个恒星在惯性坐标系中的单位矢量；为在惯性坐标系中目标天体第二个卫星到探测器的单位矢量，为目标天体第二个卫星图像中第i个恒星在惯性坐标系中的单位矢量，v_1i,v_2i,v_3i分别为θ_1i,θ_2i,θ_3i的测量误差，B_m为量测量星历误差，b_x,b_y,b_z为火星星历三轴位置误差，为星历误差模型误差；

设天文导航系统量测量Z＝[θ₁₁,θ₁₂,θ₁₃,θ₂₁,θ₂₂,θ₂₃,θ₃₁,θ₃₂,θ₃₃,b_x,b_y,b_z]^T，天文导航系统量测噪声 $v={[v_{{\mathrm{\θ}}_{11}},v_{{\mathrm{\θ}}_{12}},v_{{\mathrm{\θ}}_{13}},v_{{\mathrm{\θ}}_{21}},v_{{\mathrm{\θ}}_{22}},v_{{\mathrm{\θ}}_{23}},v_{{\mathrm{\θ}}_{31}},v_{{\mathrm{\θ}}_{32}},v_{{\mathrm{\θ}}_{33}},v_{b_x},v_{b_y},v_{b_z}]}^T$ 分别为测量θ₁₁,θ₁₂,θ₁₃,θ₂₁,θ₂₂,θ₂₃,θ₃₁,θ₃₂,θ₃₃,b_x,b_y,b_z的观测误差，由于各变量都是与时间t有关的变量，则可建立天文导航系统量测模型的表达式为：

Z(t)＝h[X(t),t]+v(t)

式中，h[X(t),t]为天文导航系统非线性连续量测函数；

⑥对步骤④和步骤⑤中的状态模型和量测模型进行离散化；

⑦天文导航量测量的获取及处理；

目标天体星历误差量测量的获取可根据步骤③获得；星光角距量测量由天文导航敏感器获得目标天体的图像，利用图像处理技术，确定天体质心的像元像线坐标；经过从像元像线坐标系到二维成像平面坐标系、从二维成像平面坐标系到敏感器测量坐标系的三次转换，确定天体及其背景恒星在敏感器坐标系中的单位矢量；最后计算天体单位矢量与背景恒星单位矢量间的星光角距；

⑧对天文导航系统进行Unscented卡尔曼滤波

根据天文导航系统状态模型、天文导航系统量测模型、天文导航敏感器获得的量测量，进行天文导航系统Unscented卡尔曼滤波，获得在目标天体惯性坐标系中表示深空探测器位置、速度和目标天体星历误差的估计状态向量和估计均方误差阵P_k；

⑨天文导航系统对目标天体星历误差进行反馈校正

将步骤⑧中获得的目标天体星历误差和目标天体星历估计均方误差阵反馈回深空探测器的状态模型中，并重新确定状态模型的模型误差协方差阵，最后将校正后的模型误差协方差阵输入至天文导航系统Unscented卡尔曼滤波中，修正下一时刻的导航结果；

⑩计算日心惯性坐标系中的探测器位置、速度估计状态向量

将步骤⑧获得的第k时刻探测器在目标天体质心惯性坐标系中的位置、速度估计状态向量转换至日心惯性坐标系中，得到日心惯性坐标系中的估计状态向量

最终输出在目标天体质心惯性坐标系中和在日心惯性坐标系中的估计状态向量和估计均方误差阵P_k，估计状态向量分别表示在目标天体质心惯性坐标系中和日心惯性坐标系中探测器的位置、速度信息，输出的估计均方误差阵P_k表示了滤波估计的性能，并将这些导航信息分别返回天文导航系中，用于k+1时刻的位置、速度导航信息，k＝1,2,...。