一种小天体绕飞轨道及目标天体参数确定方法

王卫东，张泽旭，朱圣英，崔平远

(1．哈尔滨工业大学 深空探测基础研究中心，150001哈尔滨，zexuzhang@hit．edu．cn;2．北京理工大学宇航学院，100081北京)

一种小天体绕飞轨道及目标天体参数确定方法

王卫东1，张泽旭1，朱圣英2，崔平远2

(1．哈尔滨工业大学 深空探测基础研究中心，150001哈尔滨，zexuzhang@hit．edu．cn;2．北京理工大学宇航学院，100081北京)

提出一种基于SRIF滤波器的小天体探测器绕飞轨道及目标天体参数确定方法．考虑绕飞小天体初期待估参数多、动力学环境复杂等特点，利用SRIF滤波器对组合导航数据进行处理，对探测器绕飞轨道以及目标天体引力场模型、自旋状态、星历信息等动力学参数进行估计．利用Eros433的观测数据进行数学仿真表明，本文给出的方法能够有效地计算绕飞轨道和确定目标天体的物理参数．

小天体探测;绕飞轨道;光学导航;均方根信息滤波

当探测器最终捕获目标小天体并形成稳定的绕飞轨道以后，需要确定探测器相对小天体的位置、速度以及姿态等参数，同时在此阶段，对目标天体参数的评估也是对探测器的1个重要需求，这些参数包括影响探测器轨道的小天体的引力场系数、自旋角速度、自旋轴指向以及小天体形状和星历信息等．NEAR、Rosetta 任务［1－2］证实将小天体表面存在的弹坑地表特征作为导航路标，比利用小天体形状特征具有更好的轨道确定性能．

在小天体探测器绕飞轨道确定方面，日本ISAS的Kawaguchi等研究了绕飞小行星等小天体的自主光学制导与导航技术［3－4］．美国 JPL 的Bhaskaran等研究了低成本的绕飞小天体轨道自主确定技术，提出了一种利用宽视场相机和预处理的目标小天体模型自主确定绕飞小天体轨道的算法［5－6］．J．J．Bordi利用 NEAR 任务中测量数据研究了激光测距仪的测量对Eros形状评估和探测器轨道确定的影响，并对只借助于激光测据仪导航方式的导航精度进行了分析，结果显示此时的导航精度在100 m左右［7］．

利用探测器绕飞小天体的观测数据来确定小天体模型参数的研究，主要是JPL实验室和ISAS在进行这方面的研究．JPL实验室一般利用下面两种技术，第一种技术是使用PCODP(PC Orbit Determination Program)软件处理大约30 d的探测器斜距和多普勒数据及小天体表面的光学观测数据［8－9］，第二种技术是利用另 1 个独立的软件ODP(Orbit Determination Program)处理整个轨道的斜距和多普勒数据［10］．Kominato 等［11］研究了利用在轨光学导航参数确定小天体相关物理参数的问题;Catherine等［12］利用飞越小天体的观测数据研究了小天体的质量确定问题．Johnson等［13］人研究了基于视觉分析的火星软着陆导航方法．

均方根信息滤波SRIF算法在处理包含了光学数据的观测数据时能有效地克服滤波器的发散，具有较高的数值稳健型和计算高效性［14］．本文将利用SRIF对绕飞小天体的轨道确定滤波器进行设计，考虑未建模干扰加速度的影响，通过测距数据、测速数据、VLBI数据及导航路标光学数据的处理，对探测器绕飞轨道、小天体物理参数进行确定．

1 变分方程与观测方程

定义小天体固联坐标系Σa:坐标原点oa位于小天体的质量中心，za轴沿小天体自旋轴方向，xa轴沿小天体最小惯量轴方向，ya轴定义满足右手系法则.在小天体固联坐标系下，绕飞探测器的轨道动力学方程可以表示为

式中:x、y、z为探测器的三轴位置;ω为小天体的自旋角速度;Tx、Ty、Tz为建模加速度包括控制加速度、太阳光压加速度等;nx、ny、nz为未建模的干扰加速度;V为小天体的引力位函数，利用球体调和函数表示如下:

其中:a为小天体名义半径;θ、φ为探测器所处的经纬度;Pnm为缔结勒让德多项式.

定义小天体惯性坐标系ΣA:坐标原点oA位于小天体的质量中心，zA轴沿小天体自旋轴方向，xA轴位于黄道平面内垂直zA轴，yA轴定义满足右手系法则.令小天体自旋轴在J2000坐标系下的赤经、赤纬分别为θa、φa，则小天体惯性坐标系相对J2000坐标系的坐标转换矩阵为小天体固联坐标系相对于小天体惯性系坐标转换矩阵为

其中λ为某历元时刻xa与xA之间夹角，t为当前时刻距历元时刻的时间.

探测器绕飞轨道定轨数据有地面站测轨数据和星载相机拍摄的路标图像数据，其中地面站测轨数据提供探测器绝对位置、速度信息，其本质包含了探测器运动距离信息．路标图像数据提供了探测器相对目标天体的位置信息，其本质是角度测量信息，结合地面测轨数据包含的运动距离信息和星敏感器提供的姿态信息，能够确定探测器在小天体固联坐标系下的三维位置、运动状态，同时建立小天体固联坐标系与J2000坐标系之间的转换关系．路标图像数据为小天体表面弹坑相对应的像素p、像线l坐标值，其表达式为

其中xci、yci、zci是第i个导航路标在相机坐标系下相对探测器位置矢量rci的三轴分量，f为导航相机焦距.其中

式中r，ρi分别为小天体固联坐标系下探测器和第i个导航路标的位置矢量;Cca是相机坐标系相对小天体固联坐标系的坐标转换矩阵，有Cca=CcbCbICIa成立，Ccb为导航相机坐标系在探测器本体坐标系下的安装矩阵，CbI为探测器本体坐标系相对惯性空间的姿态转换矩阵，CIa为J2000坐标系相对小天体固联坐标系姿态转换矩阵．

基于以上给出的轨道动力学方程与坐标转换关系，建立变分方程、量测方程，可以求取状态转移偏导数与各观测数据对应的观测偏导数，进而确定量测偏导数，以对待估参数进行确定．

2 待估参数与SRIF滤波器设计

2.1 待估参数选取

参考小天体绕飞段轨道变分方程和量测方程，可见直接影响观测量的因素有探测器相对目标天体的位置、速度和姿态、目标天体的运行轨道、测控站的三维位置、导航路标在目标天体固联坐标系下的位置，间接影响观测量的因素有太阳光压系数、目标天体引力场系数、目标天体自旋状态、随机加速度等动力学参数．

在探测器绕飞目标天体阶段，星敏感器与陀螺联合定姿系统能够提供探测器相对惯性空间的姿态信息，而导航路标所对应像点坐标与探测器相对目标天体的姿态相关，因此光学观测数据对目标天体相对惯性空间姿态可观．这里采用小天体固联坐标系相对J2000姿态转换矩阵表示目标天体的姿态，该姿态可以利用小天体自旋轴指向、子午线历元时刻指向、小天体旋转角速度表示．光学观测数据同时与探测器相对目标天体的位置相关，地面测控数据则与探测器相对日心惯性空间轨道状态有关，因此，结合光学观测数据与地面测控数据能够对目标天体的星历信息进行估计，这里采用历元时刻目标天体相对日心的位置、速度表示其星历状态．太阳光压力、小天体引力、向心加速度、苛氏加速度等动力学参数影响探测器的运行轨道，间接影响各观测量数据，其隐含在变分方程中，可以通过状态转移矩阵对这些动力学参数进行估计．综合如上因素，考虑科学任务和工程任务的需求，本文选取SRIF滤波器待估参数如表1所示．

表1 轨道确定滤波器估计参数

在小天体表面存在大量弹坑地形，后续自主光学导航需要选取弹坑为导航路标，同时对其三维位置信息进行确定．待估参数中的随机加速度项为探测器所受到的未建模加速度部分，主要是由于太阳光压力变化、燃料泄漏引起的，参考国外测控经验，用高斯－马尔科夫过程对其建模，该加速度包括时间相关的部分和纯粹随机的部分，其在一定时间内变化不大．SRIF滤波器采用批处理与贯序处理相结合的滤波技术，在滤波过程中，每个批处理过程中均假设未建模加速度为常值，在批与批之间的贯序处理过程中时间相关变化．

基于SRIF滤波的基本原理，考虑到实现绕飞小天体轨道确定的复杂性，将待估参数分为三类进行处理:1)P为相关过程噪声参数，即为上述中的未建模加速度;2)X为随时间变化的状态参数，但并不完全依赖于白噪声的量，如探测器的相对位置和速度、小天体星历状态等;3)Y为不随时间变化的状态参数，如探测器太阳光压系数、小天体引力场系数、小天体姿态状态、弹坑三维位置等常值参数.此时，滤波器的状态方程可表示为

式中:下标j表示第j批参数;Wj+1为滤波器的过程噪声;VP(j)为未建模加速度对状态参数的作用矩阵，可以表示为

其中ΦX和ΦXP代表了相应的状态转移矩阵.由于假设探测器所受未建模加速度满足高斯马尔科夫过程，因此式(1)中转移矩阵Mj可以表示为

其中Δtj表示第j批与第j+1批间隔时间，τi表示噪声的相关时间.

选取每步待估参数的更新值作为递推状态，令初始待估参数服从的误差协方差阵为P0，对P0阵进行Cholesky变换获取初始均方差阵S0，进而求得初始信息矩阵R0，满足

式中R0为上三角阵，对应各待估参数可以表示为

式中 zX、zP、zY为初始虚拟观测值，vX、vP、vY为初始观测噪声，由于初始改变值为零值，所以zX、zP、zY和 vX、vP、vY也均为零值.

2.2 时间更新方程

方程(1)可以等价表示为

将上式代入式(2)中，有

由于Wj+1可以认为是一独立的随机过程噪声，并服从N( 0，σ2)分布，可以用下式描述:式中rw为相应过程噪声的标准差倒数，考虑到Wj+1为零均值，将zW取为零值.综合式(4)和式(5)并代入式(3)，有

利用Householder变换［15］构造正交阵T左乘上式，使左端动态矩阵转化为上三角矩阵，上式变为

将式(6)第2行到第4行的子矩阵保留，即为时间更新到第j+1步的信息矩阵等式，上述过程为SRIF滤波器的时间更新过程，主要考虑了未建模加速度对系统估计参数的影响．

2.3 测量更新方程

利用地面测轨数据与光学数据的线性化量测偏导数，可以将测量残差表示成相应待估参数的形式

考虑测量数据的噪声特性，将上式与信息矩阵等式合并，有

其中N为观测噪声服从的误差标准差，即R=NNT，R为观测噪声服从的误差协方差阵.

利用Householder变换构造正交阵T左乘上式，使左端动态矩阵转化为上三角矩阵，则上式变为

求解式(8)，可得到第j+1步的状态参数以及相应的协方差阵，利用式(7)中量测更新的信息矩阵等式，继续进行时间更新，以便数据的贯序批处理，完成探测器轨道状态和目标天体物理参数的确定.

3 数学仿真与分析

为了验证轨道确定方案的可行性以及SRIF滤波器的性能，以Eros433为目标星，利用从2012年12月28日起的14天观测数据对其绕飞轨道与目标天体参数进行估计．小行星Eros433自旋角速度1639.4(°/d)，名义半径16 km，引力常数4.462 1×105m3/s2．设探测器绕飞轨道高度为80 km，在该轨道高度探测器受到小天体引力场高阶项的摄动已较为明显，如图1所示．在两周的引力摄动作用下，探测器轨道漂移将会达到3 km以上，因此在该轨道高度上能够对目标天体的引力场进行较高精度的评估．这里选取待估目标天体引力场6阶系数，绕飞轨道偏心率 0，轨道倾角176°，探测器面质比0. 01，光压系数为0.6．

参考国外测控经验，SRIF滤波器采用一阶高斯马尔科夫过程模拟未建模加速度，其相关时间选为2 d，量级为10－12km/s．探测器在80 km 轨道上绕飞阶段，地面站每天跟踪测轨1次，每次2 h，导航相机每天拍照60幅．斜距测量系统差、随机差均为10 m，斜距变化率测量系统差、随机差为1 mm/s，VLBI测量测角精度为5 nrad，导航相机视场角为3°×3°，分辨率为1 024 ×1 024，图像处理精度0.1 piexl，导航路标选取100个，星敏感器确定精度为50 nrad．待估参数初始不确定度及经SRIF滤波器迭代30次得到的最终确定结果如表2所示．

图1 探测器绕飞小天体Eros433轨道

表2 待估参数初始条件及结果

初始探测器相对目标天体的不确定位置精度为三轴各300 m，速度0.01 m/s，经观测处理后，位置、速度状态确定误差均得到减小，由结果可见，探测器在小天体固联坐标系下z轴的导航精度明显地要高于其它两轴，这与探测器绕飞轨道倾角有关．在本仿真中，绕飞轨道近似为赤道面轨道，导航路标观测几何关系使自旋轴方向的可观度大于其它两轴．由目标天体自旋轴、子午线与旋转角速度的估计结果可以看出，采用光学路标导航辅助的方法，可以对目标天体在绝对空间下的姿态指向进行较高精度的估计，其估计精度与星敏感器同等量级．小天体引力常数由初始2 000 m3/s2的不确定度，下降到100 m3/s2以内，估计精度达到了0.02%以内，对小天体各阶引力系数估计精度也提高了2个数量级，说明滤波器对小天体质量精确确定的同时，也能够对小天体的质量分布进行很好地估计，结合小天体形状信息，通过对小天体密度的分析，可以对其内部组成成分及物质构成提供参考．另一项科学考察数据，小天体的星历状态可以达到10 km以内，速度在0.1 m/s以内，达到与基于地面站对深空探测器进行测轨的同等精度量级．SRIF滤波器对100个导航路标的三维位置进行估计，表2中误差表示形式为位置误差平均值，其表达式为

同时导航相机可以拍摄到其它路标的大量图像数据，这些路标虽然没有在SRIF滤波器中对其位置进行估计，但是结合拍照时刻探测器的位置、姿态等信息，利用多帧图像数据能够对这些路标的位置进行估计，这些路标与SRIF滤波器中采用的导航路标联合组成了导航路标库，以支持后续探测任务的自主光学导航的实现．

4 结论

本文对探测器绕飞轨道及目标天体参数确定方法进行了研究，针对绕飞小天体初期待估参数多，动力学不确定环境复杂等特点，同时考虑未建模加速度的影响，提出了一种基于SRIF滤波器的轨道确定和目标天体参数确定方法．该方法利用地面测控数据、导航路标的光学数据对探测器绕飞轨道确定，同时对目标天体引力场模型、自旋状态、星历信息等动力学参数进行估计．仿真实验表明该方法能够有效地对绕飞轨道和目标天体参数进行计算．

［1］KONOPLIV A S，MILLER J K，OWEN W M，et al.A global solution for the gravity field，rotation，landmarks，and ephemeris of eros［J］．Icarus， 2001，160(2):289－299．

［2］ BIELE J，ULAMEC S．Capabilities of philae，the rosetta lander［J］．Space Science Reviews， 2008，138(1/2/3/4):275－289．

［3］KAWAGUCHI J，HASHIMOTO T，MISU T，et al.An autonomous optical guidance and navigation around the asteroid［C］//Proceedings of the 47th International Astronautical Congress．Beijing:［s．n．］，1996．

［4］KAWAGUCHI J，HASHIMOTO T，KUBOTA T，et al.Autonomous optical guidance and navigation strategy around a small body［J］．Journal of Guidance，Control and Dynamics， 1997，20(5):1210 －1220．

［5］BHASKARAN S，RIEDEL J E，SYNNOTT S P．Demonstration of autonomous orbit determination around small bodies［J］．Advances in the Astronautical Sciences Series， 1996，90(2):1297 －1308．

［6］TIMOTHY C D．NEAR Laser rangefinder:a tool for the mapping and topologic study of asteroid 433 Eros［J］．Johns Hopkins APL Technical Digest， 1998，19(2):142－157．

［7］BORDI J J，MILLER J K，WILLIAMS B G．Near Earth Asteroid Rendezvous(NEAR)navigation using altimeter range observations［R］．The Interplanetary Network Progress Report， 2001，IPN PR －42－146．

［8］MILLER J K，WEEKS C J，WOOD L J．Orbit determination strategy and accuracy for a comet rendezvous mission［J］．Journal of Guidance，Control and Dynamics， 1990，13:775 －784．

［9］CHESLEY S R，YEOMANS D K．Comet 9P/Tempel 1 ephemeris development for the deep impact mission［J］．Advances in the Astronautical Sciences， 2006， 2:1271－1282．

［10］MILLER J K，KONOPLIV A S，ANTREASIAN P G，et al.Determination of shape，gravity，and rotational state of asteroid 433 Eros［J］．Icarus， 2002，155(1):3 －17．

［11］KOMINATO T，MATSUOKA M，UO M．Optical hybrid navigation in Hayabusa-Approach，station keeping＆hovering［J］．Advances in the Astronautical Sciences， 2006，124(2):1753－1772．

［12］CATHERINE B S，JEREMY B J，JOHN D A．Small body mass determination from a flyby［J］．Advances in the Astronautical Sciences， 1991，75(2):885 －896．

［13］JOHNSON A E，WILLSON R，CHENG Y．Design through operation of an image-based velocity estimation system for mars landing［J］．International Journal of Computer Vision， 2007，74(3):319－341．

［14］BIERMAN G J．Factorization Methods for Discrete Sequential Estimation．New York:Academic Press，1976．

［15］ALSTON S．Householder，Unitary Triangularization of a Nonsymmetric Matrix［J］．Journal ACM， 1958，5(4)，339－342．

An algorithm of around orbit and parameters determination for small body

WANG Wei-dong1，ZHANG Ze-xu1，ZHU Sheng-ying2，CUI Ping-yuan2

(1．Deep Space Exploration Research Center，Harbin Institute of Technology，150001 Harbin，China，zexuzhang@hit．edu．cn;2．Aerospace College，Beijing Institute of Technology，100081 Beijing，China)

The around orbit parameters and small body physical parameters determination methods are studied based on Square Root Information Filter(SRIF)．Due to a great deal of parameters to be estimated and the complexity of the dynamics，the SRIF is designed to process the integrated navigation data and estimate the gravity field，rotation，and ephemeris of small body and orbit parameters of probe．The mathematic simulation experiments using the measuring data of Eros433 small body have validated the algorithm．

small body exploration;around orbit;optical navigation;Square Root Information Filter(SRIF)

文献标志码:A 文章编号:0367－6234(2011)09－0019－06

2010－04－09．

航天创新基金资助项目(CASC200902－4)，深空探测着陆与返回控制技术国防重点学科实验室，开放基金资助项目(HIT．KLOF．2009070)．

王卫东(1969—)，男，博士研究生;

崔平远(1961—)，男，教授，博士生导师．

(编辑 张 宏)