快速检索        
  武汉大学学报·信息科学版  2017, Vol. 42 Issue (8): 1035-1039

文章信息

魏二虎, 王凌轩, 张帅, 刘经南
WEI Erhu, WANG Lingxuan, ZHANG Shuai, LIU Jingnan
基于X射线脉冲星的环火探测器同步定位和授时仿真分析
Simulation Analysis of Mars Probe Autonomous Positioning and Timing Based on X-Ray Pulsars
武汉大学学报·信息科学版, 2017, 42(8): 1035-1039
Geomatics and Information Science of Wuhan University, 2017, 42(8): 1035-1039
http://dx.doi.org/10.13203/j.whugis20150437

文章历史

收稿日期: 2016-03-29
基于X射线脉冲星的环火探测器同步定位和授时仿真分析
魏二虎1, 王凌轩1, 张帅1, 刘经南2     
1. 武汉大学测绘学院, 湖北 武汉, 430079;
2. 武汉大学卫星导航定位中心, 湖北 武汉, 430079
摘要:利用X射线脉冲星观测量不仅能够为深空探测器提供稳定可靠的位置、速度和姿态信息,还能够对探测器上的原子钟进行长期误差校正。火星探测任务中,存在探测器飞行时间长,钟差随时间累积严重降低脉冲星自主定位系统的精度等问题。针对此,提出了一种环火探测器同步定位和授时方法,对钟差进行建模,将钟差模型参数纳入到滤波系统中,利用自适应扩展卡尔曼滤波器同时计算探测器的位置、速度和钟差模型参数。仿真实验表明,经滤波校正后的钟差保持在300 ns以内,探测器的位置和速度精度分别优于200 m和0.03 m/s,精度得到明显改善。
关键词环火探测器     X射线脉冲星     自适应扩展卡尔曼滤波     定位     授时    
Simulation Analysis of Mars Probe Autonomous Positioning and Timing Based on X-Ray Pulsars
WEI Erhu1, WANG Lingxuan1, ZHANG Shuai1, LIU Jingnan2     
1. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China;
2. GNSS Research Center, Wuhan University, Wuhan 430079, China
First author: WEI Erhu, PhD, professor, specializes in space geodesy and geodynamics. E-mail: ehwei@sgg.whu.edu.cn
Foundation support: The National Natural Science Foundation of China, No. 41374012
Abstract: Measurement based on X-Ray Pulsars is not only able to provide stable and reliable position, velocity and attitude information for deep space probe, but also can correct the long term bias of the atomic clock in the probe. Especially for the Mars exploring missions, the flight time of Mars probe is fairly long so that the clock bias is accumulating continuously, which will greatly damage the accuracy of the Mars probe autonomous positioning and timing system. In this paper, a novel synchronous Positioning/Timing method for Mars probe is proposed based on the principle of synchronous positioning and timing by X-Ray pulsars. And the clock correction is modeled accordingly. Pulsar optimization factor is used to conduct the optimal selection of pulsars. Then the Positioning/Timing model is presented and the clock corrections based on the new clock model are estimated together with the position and velocity of the Mars probe by using the adaptive extended Kalman filter. So that the position, velocity and clock model parameters of Mars probe could be calculated synchronously. Simulation experiments are carried out based on the proposed method. Results show that the clock bias can maintain in 300ns, while the position and velocity accuracy of Mars probe is better than 200m and 0.03m/s after corrected, respectively. Conclusions could be drawn that the proposed method of Mars Probe Autonomous Positioning and Timing Based on X-Ray Pulsars could overcome the initial state error effectively and obtain more accurate filter value of position, velocity and clock parameters for Mars probe.
Key words: X-ray pulsars     Mars probe     adaptive extended Kalman filter     positioning     timing    

火星探测任务中航天器的精确导航定位至关重要,而火星和地球之间的遥远距离和探测器所处的复杂太空环境对基于地面测控通讯站的传统导航方法而言是一个巨大挑战,航天器的自主导航将是突破该瓶颈的关键技术之一[1]

X射线脉冲星自主导航技术具有可靠性高、自主性强、精确性良好和使用范围广等优点,在火星探测器自主导航方面极具应用潜力。国内外一些学者在利用X射线脉冲星观测量为执行火星探测任务的飞行器提供自主导航服务方面开展了相关研究,包括相关时空基准的定义和转换、适用于火星探测器的脉冲星测量模型的建立、火星探测器轨道动力学研究和导航滤波算法等内容[2, 3]。本文在此基础上,针对初始状态误差较大会导致普通的滤波算法收敛较慢甚至发散,以及火星探测任务周期长,探测器上搭载的原子钟的钟差随时间累积,严重损害了自主定位系统的精度等问题,提出了利用X射线脉冲星观测量同步滤波估计探测器的位置、速度和钟差模型参数的同步定位和授时方法,并利用设计的自适应扩展卡尔曼滤波器进行了仿真实验。

1 脉冲星同步定位和授时模型

由X射线脉冲星自主导航原理可知,该方法是一种测定探测器相对于某一特定参考系的相对测量方法。针对火星轨道探测器,该方法选择火星质心作为基准点,建立J2000.0火星质心天球坐标系。此时,X射线脉冲星观测量即脉冲信号到达探测器和火星质心的时间延迟τ可以表示为探测器在J2000.0火星质心天球坐标系中的位置参数r和钟差ΔtS的函数。应用于工程计算的简化脉冲星轨道测量模型[1]如式(1) 所示:

$ \tau = \frac{1}{c}\mathit{\boldsymbol{n}} \cdot \mathit{\boldsymbol{r + }}\Delta {t_S} + v $ (1)

式中,n为脉冲星的方向矢量;v为测量噪声。若同时获得多于3颗脉冲星的测量信息,则可计算出探测器的位置参数r和钟差ΔtS

在航天器导航中,通常结合探测器的轨道动力学信息和观测信息,利用相应的卡尔曼滤波算法进行状态参数的预测。

1.1 轨道动力学模型

在导航滤波算法中,探测器的轨道动力学模型主要是在获取新的观测量之前,对探测器轨道进行预报。火星轨道探测器受到的主要作用力是火星的中心引力,对于中高轨道探测器而言,研究表明火星非球型摄动中的J2项为主要摄动项[2, 4]

在J2000.0火星质心天球坐标系中,火星轨道探测器的动力学模型可表示为:

$ \mathit{\boldsymbol{\ddot r = }} - \frac{{{\mu _M}\mathit{\boldsymbol{r}}}}{{{r^3}}} + {\mathit{\boldsymbol{a}}_{{\mathit{\boldsymbol{J}}_2}}} + \mathit{\boldsymbol{\omega }} $ (2)

式中,r=[x y z]T为探测器位置矢量;$\mathit{\boldsymbol{\ddot r}}$为其二阶导数,r=||r||; μM是火星引力常数;w为满足高斯分布的随机噪声;aJ2代表探测器受到的J2项摄动力加速度。

$ {\mathit{\boldsymbol{a}}_{{\mathit{\boldsymbol{J}}_2}}} = \left[ {\begin{array}{*{20}{c}} { - \frac{{{\mu _M}x}}{{{\mathit{\boldsymbol{r}}^3}}}{\mathit{\boldsymbol{J}}_2}{{\left( {\frac{{{R_M}}}{\mathit{\boldsymbol{r}}}} \right)}^2}\frac{3}{2}\left( {1 - 5\frac{{{z^2}}}{{{\mathit{\boldsymbol{r}}^2}}}} \right)}\\ { - \frac{{{\mu _M}y}}{{{\mathit{\boldsymbol{r}}^3}}}{\mathit{\boldsymbol{J}}_2}{{\left( {\frac{{{R_M}}}{\mathit{\boldsymbol{r}}}} \right)}^2}\frac{3}{2}\left( {1 - 5\frac{{{z^2}}}{{{\mathit{\boldsymbol{r}}^2}}}} \right)}\\ { - \frac{{{\mu _M}z}}{{{\mathit{\boldsymbol{r}}^3}}}{\mathit{\boldsymbol{J}}_2}{{\left( {\frac{{{R_M}}}{\mathit{\boldsymbol{r}}}} \right)}^2}\frac{3}{2}\left( {1 - 5\frac{{{z^2}}}{{{\mathit{\boldsymbol{r}}^2}}}} \right)} \end{array}} \right] $ (3)

式中,RM为火星平均半径。

在给定探测器初始轨道参数的情况下,即可利用轨道动力学模型进行数值积分计算下一时刻的探测器状态。

1.2 钟差模型

一般情况下原子钟的钟差可以用形如式(4) 的二阶多项式表示[5]

$ \begin{array}{l} {x_1}\left( {{t_{k + 1}}} \right) = {x_1}\left( {{t_k}} \right) + \tau \cdot {x_2}\left( {{t_k}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{{\tau ^2}}}{2}{x_3}\left( {{t_k}} \right) + \omega \left( k \right) \end{array} $ (4)

式中,x1(tk)、x2(tk)和x3(tk)分别表示tk时刻原子钟的钟差、频率漂移和频率漂移变化率,它们组成了钟差参数;ω(k)为相应的原子钟噪声。

以时钟钟差、时钟频率漂移和频率漂移变化率为状态量的钟差状态方程为:

$ {\left[ {\begin{array}{*{20}{c}} {{x_1}}\\ {{x_2}}\\ {{x_3}} \end{array}} \right]_{k + 1}} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k + 1,k}}{\left[ {\begin{array}{*{20}{c}} {{x_1}}\\ {{x_2}}\\ {{x_3}} \end{array}} \right]_k} = {\mathit{\boldsymbol{\omega }}_k} $ (5)

式中,ωk为钟差噪声向量;状态转移矩阵Φk+1, k为:

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k + 1,k}} = \left[ {\begin{array}{*{20}{c}} 1&\varepsilon &{{\varepsilon ^2}/2}\\ 0&1&\varepsilon \\ 0&0&1 \end{array}} \right] $ (6)

式中,ε为时间间隔;ωk满足以下统计学模型:

$ \begin{array}{*{20}{c}} {{Q_k}\left( \varepsilon \right) = E\left[ {{\mathit{\boldsymbol{\omega }}_k}\mathit{\boldsymbol{\omega }}_k^{\rm{T}}} \right] = }\\ {\left[ {\begin{array}{*{20}{c}} {{q_1}\varepsilon + \frac{1}{3}{q_2}{\varepsilon ^3} + \frac{1}{{20}}{q_3}{\varepsilon ^5}}&{\frac{1}{2}{q_2}{\varepsilon ^2} + \frac{1}{8}{q_3}{\varepsilon ^4}}&{\frac{1}{6}{q_3}{\varepsilon ^3}}\\ {\frac{1}{2}{q_2}{\varepsilon ^2} + \frac{1}{8}{q_3}{\varepsilon ^4}}&{{q_2}\varepsilon + \frac{1}{3}{q_3}{\varepsilon ^3}}&{\frac{1}{2}{q_3}{\varepsilon ^2}}\\ {\frac{1}{6}{q_3}{\varepsilon ^3}}&{\frac{1}{2}{q_3}{\varepsilon ^2}}&{{q_3}\varepsilon } \end{array}} \right]} \end{array} $ (7)

式(7) 中,q1q2q3为连续过程噪声功率谱密度。

1.3 状态方程

在J2000.0火星质心天球坐标系中,设探测器位置、速度以及3个钟差参数组成火星轨道探测器同步定位和授时系统状态向量为:

$ \mathit{\boldsymbol{x}}\left( t \right) = {\left[ {x\;y\;z\;{v_x}\;{v_y}\;{v_z}\;{x_1}\;{x_2}\;{x_3}} \right]^{\rm{T}}} $ (8)

结合火星轨道探测器动力学模型(2) 和钟差模型(4) 可得同步定位和授时系统状态方程为:

$ \mathit{\boldsymbol{\dot x}}\left( t \right) = \mathit{\boldsymbol{F}}\left( t \right) \cdot \mathit{\boldsymbol{x}}\left( t \right) + \mathit{\boldsymbol{W}}\left( t \right) $ (9)

式中,W(t)为系统状态噪声矩阵,状态方程系数矩阵F(t)为:

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{F}}\left( t \right) = }\\ {\left[ {\begin{array}{*{20}{c}} 0&0&0&1&0&0&0&0&0\\ 0&0&0&0&1&0&0&0&0\\ 0&0&0&0&0&1&0&0&0\\ {\frac{{\partial {a_x}}}{{\partial x}}}&{\frac{{\partial {a_x}}}{{\partial y}}}&{\frac{{\partial {a_x}}}{{\partial z}}}&0&0&0&0&0&0\\ {\frac{{\partial {a_x}}}{{\partial x}}}&{\frac{{\partial {a_y}}}{{\partial y}}}&{\frac{{\partial {a_z}}}{{\partial z}}}&0&0&0&0&0&0\\ {\frac{{\partial {a_z}}}{{\partial x}}}&{\frac{{\partial {a_z}}}{{\partial y}}}&{\frac{{\partial {a_z}}}{{\partial z}}}&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&1&0\\ 0&0&0&0&0&0&0&0&1\\ 0&0&0&0&0&0&0&0&0 \end{array}} \right]} \end{array} $ (10)
1.4 观测方程

若同一时段观测了p颗脉冲星,参考脉冲星测量模型(1),则环火探测器自主定位和授时系统的测量方程为:

$ \mathit{\boldsymbol{z = }}{\left[ {\Delta {t_1} \cdots \cdots \Delta {t_p}} \right]^{\rm{T}}} = \mathit{\boldsymbol{H}} \cdot \mathit{\boldsymbol{x + \eta }} $ (11)

式中,z表示脉冲星观测量;方程矩阵H为:

$ \mathit{\boldsymbol{H = }}\left[ \begin{array}{l} \begin{array}{*{20}{c}} {\frac{{{\mathit{\boldsymbol{n}}_1}}}{c}}&0&0&0&1&0&0 \end{array}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots \\ \begin{array}{*{20}{c}} {\frac{{{\mathit{\boldsymbol{n}}_p}}}{c}}&0&0&0&1&0&0 \end{array} \end{array} \right] $ (12)

完成上述同步定位和授时系统状态方程以及观测方程的建立后,便可以利用相应的自适应扩展卡尔曼滤波算法进行状态参数的解算[6, 7]

2 仿真实验与结果分析 2.1 导航脉冲星优选

在选择导航脉冲星时,除了应选择经过可见性条件检查的脉冲星,还要考虑脉冲星的测量误差及脉冲星组合的位置结构对定位精度的影响[8, 9]。实验中,利用脉冲星优选因子(POSI)来定量计算脉冲星测量误差和位置分布对定位精度的影响值:

$ \begin{array}{*{20}{c}} {{\rm{Cov}}\left( {\mathit{\boldsymbol{\hat x}}} \right) = E\left( {\mathit{\boldsymbol{\hat x}}{{\mathit{\boldsymbol{\hat x}}}^{\rm{T}}}} \right) = }\\ {{{\left( {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{H}}} \right)}^{ - 1}}{\mathit{\boldsymbol{H}}^{\rm{T}}}E\left( {\Delta \mathit{\boldsymbol{\rho }}\Delta {\mathit{\boldsymbol{\rho }}^{\rm{T}}}} \right)\mathit{\boldsymbol{H}}{{\left( {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{H}}} \right)}^{ - 1}}}\\ {{{\left( {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{H}}} \right)}^{ - 1}}{\mathit{\boldsymbol{H}}^{\rm{T}}}{\rm{diag}}\left( {\sigma _1^2,\sigma _2^2, \cdots \sigma _p^2} \right)\mathit{\boldsymbol{H}} \cdot }\\ {{{\left( {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{H}}} \right)}^{ - 1}} = {D_p}} \end{array} $ (13)

式中,$\mathit{\boldsymbol{\hat x}}$为状态向量x的修正量;Δρ为相应的脉冲星测量误差,其方差矩阵为diag(σ12, σ22, …σp2);H为简化的脉冲星观测方程矩阵:

$ \mathit{\boldsymbol{H = }}{\left[ {{\mathit{\boldsymbol{n}}_1}\;{\mathit{\boldsymbol{n}}_2}\; \cdots \;{\mathit{\boldsymbol{n}}_p}} \right]^{\rm{T}}} $ (14)

式中,n表示相应脉冲星的方向矢量。

此次环火探测器同步定位和授时实验中,我们从澳大利亚脉冲星数据库中挑选了14颗信噪比高、轮廓清晰的脉冲星,利用式(13) 进行了定量优选,最终确定B1937+21、B1821-24、B0531+21、B0540-69、B1957+20、B0614+091作为导航脉冲星,相应参数列于表 1中。

表 1 脉冲星测距精度表 Table 1 Pulsar Ranging Accuracy
脉冲星编号 测距精度σ/m
B1937+21 689.5
B1821-24 495.8
B0531+21 88.3
B0540-69 1 493.0
B1957+20 932.0
B0614+091 282.0
2.2 仿真环境设置

利用Satellite Tool Kit(STK)软件[10]仿真得到火星轨道探测器的“标称轨道”,采样间隔为3 600 s,总采样数为2 400个,初始轨道参数为:a=15 000 km, e=0.005, i=30°,Ω=30°,ω=30°,M0=0°。

在观测量模拟过程中,原子钟的钟差由钟差模型式(4) 模拟得到,初始钟差设为2×10-6 s,时钟频率漂移设为4.0×10-11 s/s, 时钟漂移变化率设为6.0×10-18 s/s2图 1为模拟的钟差累积示意图。从图 1中可知,若不进行定期校正,钟差将严重损害定位精度。

图 1 模拟钟差累积示意图 Figure 1 Schematic of Simulate Clock Correction Accumulation

其次,在观测量模拟过程中加入了-2σ~2σ的随机噪声,σ为相应脉冲星测距精度中误差。

导航滤波的初始误差设置为(10 km, 10 km, 10 km)、(2 m/s, 2 m/s, 2 m/s)和(2×10-7s, 0, 0),时钟频率漂移及其变化率同上。

2.3 实验结果和分析

完成上述模型建立及观测量模拟后便可以利用相应的自适应扩展卡尔曼滤波(AEKF)算法[11-13]进行导航滤波计算。

图 2为同步定位和授时算法滤波修正后的钟差示意图。从图 2中可以看出,经过15个历元左右滤波误差迅速收敛至400 ns以内,最终稳定在300 ns。

图 2 滤波修正后的钟差示意图 Figure 2 Clock Correction After Filtering

图 3为滤波修正后的探测器位置误差示意图,图 4为相应的速度误差示意图。从图 3~4中可知,探测器在3个坐标轴方向上的位置、速度误差经几个历元后便可以迅速收敛到1 000 m、0.1 m/s以内,经过十几个历元以后稳定在200 m、0.03 m/s以内。

图 3 滤波修正后的位置误差示意图 Figure 3 Position Error After Filtering
图 4 滤波修正后的速度误差示意图 Figure 4 Velocity Error After Filtering

图 3图 4可知,3个坐标轴方向上的滤波位置、速度精度水平相当,没有出现自主定位系统中X轴方向滤波精度差于Y、Z轴方向的精度情况[5]

统计结果表明,基于X射线脉冲星的环火探测器同步定位和授时滤波算法能够有效地克服初始状态误差的影响,获得较高精度的探测器位置、速度和钟差参数滤波值,验证了该算法在自主定位和授时方面的高效性。

3 结语

本文研究了基于X射线脉冲星的火星轨道探测器同步定位和授时原理,建立了相应的同步定位和授时模型,利用POSI进行导航脉冲星的优选,并开展了基于自适应扩展卡尔曼滤波的仿真实验。实验表明, 同步定位授时算法能够有效地校正原子钟钟差,削弱钟差对脉冲星观测量的影响,改善探测器位置和速度参数的精度。

参考文献
[1] Wu Weiren, Liu Wangwang, Tang Yuhua, et al. Deep Space Exploration and Several Key Technology Development Trends[C]. The 10th Annual Academic Conference of the China Aerospace Association Deep Space Exploration Technology Committee, Beijing, 2013 ( 吴伟仁, 刘旺旺, 唐玉华, 等. 深空探测及几项关键技术发展趋势[C]. 中国宇航协会深空探测技术专业委员会第十届学术年会, 北京, 2013 )
[2] Wei Erhu, Jin Shuanggen, Zhang Qi, et al. Autonomous Navigation of Mars Probe Using X-Ray Pulsars:Modeling and Results[J]. Advances in Space Research, 2013, 51(5): 849–857 DOI:10.1016/j.asr.2012.10.009
[3] Wu Weiren, Wang Dayi, Ning Xiaolin. Autonomous Navigation Principle and Technology of Deep Space Detector[M]. Beijing: China Space Press, 2011 ( 吴伟仁, 王大轶, 宁晓琳. 深空探测器自主导航原理与技术[M]. 北京: 中国宇航出版社, 2011 )
[4] Zhou Chuihong, Yu Shengxian, Liu Lin. Analytical Solutions of Coupled Perturbation of Tesseral Harmonic Terms of Mars Orbiters Under Nonspherical Gravitational Potential[J]. Acta Astronomica Sinica, 2012, 53(3): 205–212 ( 周垂红, 喻圣贤, 刘林. 火星非球形引力位田谐项联合摄动分析解[J]. 天文学报, 2012, 53(3): 205–212. )
[5] Shuai Ping. Principle and Method of X-Ray Pulsar Navigation System[M]. Beijing: China Space Press, 2009 ( 帅平. X射线脉冲星导航系统原理与方法[M]. 北京: 中国宇航出版社, 2009 )
[6] Wei Erhu, Zhang Shuai, Wei Jianan. Autonomous Navigation of Mars Probe Using X-Ray Pulsar Based on a Modified AEKF[J]. Geomatics and Information Science of Wuhan University, 2015, 40(5): 701–705 ( 魏二虎, 张帅, WeiJianan. 自适应滤波在环火探测器脉冲星自主导航中的应用[J]. 武汉大学学报·信息科学版, 2015, 40(5): 701–705. )
[7] Xia Qijun, Sun Youxian, Zhou Chunhui. An Optimal Adaptive Algorithm for Fading Kalman Filter and Its Application[J]. Acta Automatica Sinica, 1990, 16(3): 210–216 ( 夏启军, 孙优贤, 周春晖. 渐消卡尔曼滤波器的最佳自适应算法及其应用[J]. 自动化学报, 1990, 16(3): 210–216. )
[8] Mao Yue, Song Xiaoyong, Feng Laiping. Visibility Analysis of X-Ray Pulsar Navigation[J]. Geomatics and Information Science of Wuhan University, 2009, 34(2): 222–225 ( 毛悦, 宋小勇, 冯来平. X射线脉冲星导航可见性分析[J]. 武汉大学学报·信息科学版, 2009, 34(2): 222–225. )
[9] Chu Yonghui, Wang Dayi, Huang Xiangyu. Selection of the Beacon Pulsars in Pulsars Navigation[J]. Chinese Space Science and Technology, 2011, 31(5): 64–69 ( 褚永辉, 王大轶, 黄翔宇. 脉冲星导航中最优脉冲星组合选取方法[J]. 中国空间科学技术, 2011, 31(5): 64–69. )
[10] Yang Ying, Wang Qi. The Application of STK in Computer Simulation[M]. Beijing: National Defense Industry Press, 2005 ( 杨颖, 王琦. STK在计算机仿真中的应用[M]. 北京: 国防工业出版社, 2005 )
[11] Fu Mengyin, Deng Zhihong, Yan Liping. Theory and Application of Kalman Filter in Navigation System[M]. 2th Edition. Beijing: Science Press, 2010 ( 付梦印, 邓志红, 闫莉萍. Kalman滤波理论及其在导航系统中的应用[M]. 2版. 北京: 科学出版社, 2010 )
[12] Liu L, Tang J S. Orbit Variation Characteristics of the Mars' Orbiters[J]. Journal of Astronautics, 2008, 29(2): 461–466
[13] Wang Y, Zheng W, Sun S, et al. X-Ray Pulsar-Based Navigation Using Time-Differenced Measurement[J]. Aerospace Science & Technology, 2014, 36: 27–35