X射线脉冲星是宇宙中高速旋转并周期性地辐射电磁脉冲的中子星。由于具有高度稳定的旋转周期, 使得脉冲星所辐射脉冲的周期及轮廓具有较强的不变性和较高的可识别性。X射线脉冲星导航系统就是利用其稳定的X射线脉冲信号作为观测量进行航天器导航的系统, 它具有隐蔽性高、自主性强、抗干扰性好等优点。近年来关于X射线脉冲星导航的相关研究也日益深入。美国国防高等研究计划局早在2004年就已经正式提了X射线脉冲星导航计划—XNAV计划, 并于2005年开始研制[1], 欧洲航空局也在2005年制定了类似的“深空探测器脉冲星导航研究(deep space vessel pulsar navigation study)”计划。后续俄罗斯、日本等也都相继开展了相关的研究项目。中国在2016年11月10日于酒泉卫星发射中心成功发射了第一颗脉冲星导航试验卫星XPNAV-1[2]。
X射线脉冲星导航的测量原理是根据航天器与太阳系质心(SSB)在脉冲星方向的多普勒延迟得到航天器在SSB坐标系中的状态[3]。观测量是脉冲到达航天器与SSB的时间差或其与光速的乘积。但在宇宙中, 脉冲信号的传播会受天体运动、相对论效应等的影响而产生变化, 也就在观测方程之中引入了高阶项[4]。在当前研究中, 扩展卡尔曼滤波(EKF)算法因具有工程实现简单、算法成熟、实时性好等特点而得到广泛应用。但是在X射线脉冲星导航中应用EKF算法时却常常因为线性化需要而忽略观测模型高阶项的影响[5-8], 这实际会给导航结果带来较大截断误差甚至导致发散。这些高阶项的变化大多是慢时变的, 短时间内可以认为是常值并进行补偿, 长时间来看却有着较大的变化波动。对于X射线脉冲星导航在深空探测器、长期在轨卫星上的应用影响较大, 有必要通过在线估计的方法予以补偿。
针对以上问题, 本文以地球长期在轨卫星为背景, 通过对高阶项的产生机理及数值量级进行仿真分析, 提出一种改进的线性观测方程。该方程既不需要额外增加状态量, 也满足以EKF为基础的各种拓展滤波算法应用于脉冲星导航时对观测方程线性化的要求。并通过对比不同时段X射线脉冲星导航系统的仿真结果, 证明了该方程的有效性。
1 典型脉冲星导航空间方程 1.1 状态方程X射线脉冲星导航的状态方程是基于航天器在地球质心惯性系中的动力学方程建立起来的。状态量选为航天器在地球质心惯性系中的位置和速度。在考虑二阶带谐项时, 其状态方程为[9]
|
(1) |
式中:矢径
则可将其状态方程简写为
|
(2) |
式中:w(t)为系统噪声矩阵, 且
|
(3) |
本文X射线脉冲星导航中, 观测量选用X射线脉冲信号到达航天器和SSB的时间差Δt与光速c的乘积。如果不考虑其他物理效应的影响, 其简化的观测模型如下[10]:
|
(4) |
式中:cΔt(i)为第i颗脉冲星对应的观测量; ni为第i颗脉冲星的方向矢量, 下标i为脉冲星的编号; rsat为卫星在SSB坐标系中的位置矢量。由于状态方程中使用的状态量是航天器在地球质心惯性坐标系中的位置和速度, 所以可将rsat表示为rsat=rE+r, rE为SSB坐标系中地球的位置, 可以通过地球星历预报得到, r为卫星在地球质心坐标系的位置矢量。所以, X射线脉冲星的观测方程可以简写为
|
(5) |
式中:yk为观测量; Vk为系统的量测噪声; hk=[hk(1) hk(2) … hk(i)], 且hk(i)=[niT 01×3]。
2 观测模型高阶项分析如图 1所示, 简化的线性观测模型中仅有几何延迟项, 即nrsat=cΔt, n为脉冲星的方向矢量。脉冲信号在宇宙中传播时理论上能够对传播时间产生影响的还有相对论效应、色散延迟等[4]。但是考虑到当前光子到达时间的测量精度在微秒量级, 理论精度应达到10-8 s, 所以能够产生截断误差的因素主要考虑地球周年运动导致的周年视差效应和恒星引力场作用下光线弯曲产生的引力延迟效应[11]。
|
| OS-XSYSZS—太阳质心坐标系; OSSB-XSSBYSSBZSSB—太阳系质心坐标系; b—太阳质心到SSB的矢量; rE—地球在太阳系质心坐标系中的位置。 图 1 脉冲星导航基本原理 Fig. 1 Basic principle of pulsar navigation |
假设观测者在地球上对恒星进行观测, 由于地球的公转, 恒星在天球上的位置也会发生改变。以太阳上观测的恒星在天球上位置作为其平均位置, 从地球上观测的位置为其实际位置。由于地球围绕太阳公转, 2个位置会存在一定的偏差, 也就是所谓的恒星周年视差。周年视差在日地连线(太阳与地球连线)同星日连线(脉冲星与太阳连线)垂直时达到极大值, 如图 2所示。
|
| 图 2 周年视差效应原理 Fig. 2 Principle of annual parallax effect |
理论上而言, 被观测恒星距离太阳的距离越远, 周年视差效应会越不明显, 所以脉冲星导航中此项的时间量也是一个很小的数值。但是在X射线脉冲星导航的解算过程中, 该项会直接与光速相乘, 而且由几何原理可知, 该项产生的误差投影到3个坐标方向上之后, 各方向误差的代数和会变大, 如式(6)所示:
|
(6) |
式中:cΔt1(ij)为第i颗脉冲星在j方向上造成的位置估计误差; Δt1(i)为第i颗脉冲星由于周年视差效应造成的时间差; x、y和z分别为坐标轴的3个方向。
在短时间内周年视差效应引起的偏差变化不明显, 但是一年内的变化幅值却比较大。如果X射线脉冲星导航系统长时间的运行, 那该项的影响必须予以考虑。关于周年视差效应在X射线脉冲星观测中的影响的具体推导属于天文学范畴, 在此直接由式(7)给出[11-12]:
|
(7) |
式中:D0为脉冲星距离SSB的距离。
2.2 引力延迟效应引力延迟效应也称Shapiro效应, 是在太阳系内验证广义相对论的4个经典试验之一。它是指当射电或者光信号经过大质量的天体时, 受天体重力场的影响其速度方向会发生变化, 传播的行程增加, 传播到目的地或往返的时间也会因此增加, 如图 3所示[13]。所以如果观测者同样在地球上, 当地球位于星日连线上时, 引力延迟效应最明显; 当地球位于星日连线的延长线上时, 引力延迟效应相互抵消, 影响效果最小。
引力延迟效应对导航结果的影响同周年视差效应相同, 在此不再赘述。根据天文学知识, 可以得到太阳系内天体引力延迟效应的总和为[14]
|
(8) |
式中:Ml为对应天体的质量; bl为SSB相对于第l颗天体的位置; dlp为天体l中心到脉冲星的位置矢量; pl为航天器相对天体l中心的位置矢量。
由于太阳系内太阳质量最大, 产生的引力延迟也就最大, 约为微秒量级[15]。而其他天体的引力延迟均不大于10-8 s, 所以受光子测量精度影响可以不予考虑。若仅考虑太阳引力场的作用, 则式(8)可以简化为
|
(9) |
式中:Msun为太阳质量。
2.3 数值分析通过以上分析, 能够产生截断误差的主要为地球运动的周年视差效应和太阳的引力延迟效应。利用STK9.0版本软件按照表 1中的轨道参数分别仿真产生一天及一年内的高精度HPOP卫星运行轨道, 并以表 2[16]中B1821-24脉冲星的参数为例计算式(7)和式(9)的结果, 太阳质量Msun=1.989 1×1030 kg, 万有引力常数G=6.67×10-11 m3·kg-1·s-2。一天和一年内高阶项变化情况结果分别如图 4和图 5所示。
| 参数 | 数值 |
| 半长轴/km | 7350 |
| 离心率 | 0 |
| 轨道倾角/(°) | 45 |
| 升交点赤经/(°) | 0 |
| 近地点幅角/(°) | 30 |
| 初始真近点/(°) | 0 |
| 起始时刻 | 2015-10-17 T04:00:00 |
| 结束时刻 | 2016-10-16 T04:00:00 |
|
| 图 4 一天内高阶项变化情况 Fig. 4 Variation of higher order terms in a day |
|
| 图 5 一年内高阶项变化情况 Fig. 5 Variation of higher order terms in a year |
通过分析可以发现:周年视差效应和太阳引力延迟效应的变化幅度明显, 最大与最小值能相差10倍之多, 一年内周期变化明显。最大值最小点计算的结果与第2.1节及2.2节的理论分析相吻合。
3 线性观测方程的改进由于滤波周期内卫星绕地球的运动相对于卫星到SSB的距离|rsat|来说非常小, 所以结合式(7), 在k时刻可对脉冲星导航中的周年视差效应做如下化简:
|
(10) |
式中:rsat(k-1)为k-1时刻卫星在SSB中的位置; rsat(k)为k时刻卫星在SSB中的位置。
引力延迟效应的影响量级最大, 同理也可对k时刻的式(9)进行如下简化处理:
|
(11) |
令
|
(12) |
则将式(11)在k-1时刻对rsat进行一阶泰勒展开:
|
(13) |
式中:Δ为泰勒展开中二阶以上高阶项。
省略二阶以上高阶项后可得
|
(14) |
综上所述, 结合式(5)、式(10)、式(12)、式(14), 可将改进的带有高阶项的线性观测方程写为
|
(15) |
式中: 
|
|
为验证本文所提出的改进观测方程对导航结果的提高效果, 本文利用表 1中的数据仿真产生卫星运行轨道。所用的X射线背景流量BX=0.005 ph·cm-2·s-1, 探测器面积A=1 m2, 脉冲星导航周期为60 s, 导航采用B0531+21、B1821-24及B1937+21三颗脉冲星, 具体参数见表 3[16]。地球引力常数为3.986 004 418×1014 N·m2/kg, 光速为3×108 m/s, 重力二阶带谐系数为0.001 082 63, 地球赤道半径为6.378 137×106 m。本文中脉冲星导航的其他初始条件设置如下:
1) 初始误差为
δx=(1 000 m, 1 000 m, 1 000 m, 20 m/s, 20 m/s, 20 m/s)
2) 初始误差协方差为
P(0)=diag[1 0002 m2, 1 0002 m2, 1 0002 m2, 202 m2/s2, 202 m2/s2, 202 m2/s2]
3) 系统噪声协方差为
Qk=diag[q12, q12, q12, q22, q22, q22]
式中:q1=8 cm; q2=0.05 mm/s。
4) 量测噪声协方差为
Rk=diag[1092 m2, 3252 m2, 3442 m2]
通过第2节对高阶项的分析可以得知, 假如从卫星上对脉冲星进行观测, 在一年中有4个点比较特殊, 都是在日卫连线(太阳与卫星连线)与星日连线相共线或垂直时。由于导航的截断误差同时受3颗脉冲星的高阶项变化影响, 所以3颗脉冲星高阶项总和的周年变化在脉冲星导航中是有实际影响意义的。
3颗脉冲星高阶项总和的周年变化情况如图 6所示。为了验证改进的线性观测方程的有效性, 本文选取其最大值和最小值点, 即2015年10月17日与2015年12月29日进行仿真运算, 仿真结果见图 7~图 10。同时为了更好地评价算法的估计性能, 本文采用均方根误差(RMSE)作为导航误差的计算公式[17], 统计区间为一天中的(900, 1 440) min区间, RMSE的表达式为
|
(16) |
|
| 图 6 全年3颗脉冲星高阶项总和变化 Fig. 6 Whole year variation of summation of three pulsars' higher order terms |
|
| 图 7 2015年10月17日简化线性观测方程估计结果 Fig. 7 Estimation results of simplified linear measurement equation on October 17, 2015 |
|
| 图 8 2015年10月17日改进线性观测方程估计结果 Fig. 8 Estimation results of improved linear measurement equation on October 17, 2015 |
|
| 图 9 2015年12月29日简化线性观测方程估计结果 Fig. 9 Estimation results of simplified linear measurement equation on December 29, 2015 |
|
| 图 10 2015年12月29日改进线性观测方程估计结果 Fig. 10 Estimation results of improved linear measurement equation on December 29, 2015 |
式中:Δrk为第k时刻真实的轨道位置与滤波器估计位置之间的距离。
为了降低随机因素对结果的影响, 将2个线性观测方程分别独立运算50次并对导航结果求平均值, 结果信息统计如表 4所示。通过运行结果的图表来看, 在考虑高阶项的情况下使用简化线性观测方程进行EKF解算时, 忽略高阶项影响而产生的较大截断误差已经导致滤波发散, 而改进的线性观测方程却仍然能够较好较快地将位置估计误差收敛到250 m以内, 速度估计误差2 m/s以内, 并且在高阶项的最大值点与最小值点的位置估计误差仅相差2.6 m, 速度误差仅相差0.038 2 m/s。方程对高阶项的变化表现出较好的鲁棒性。通过运算时间反应运算量, 在60 s更新周期的条件下, 用配置i3处理器的台式计算机运行2012A版本MATLAB进行导航解算, 2种方程进行一天内的导航解算, 总运行时间分别为2.054 4 s和1.749 0 s, 平均每次计算仅相差0.000 2 s。
| 日期 | 简化观测方程 | 改进观测方程 | |||
| 位置误差/km | 速度误差/(m·s-1) | 位置误差/km | 速度误差/(m·s-1) | ||
| 2015-10-17 (最大值点) |
51.631 1 | 120.003 9 | 0.244 5 | 1.871 5 | |
| 2015-12-29 (最小值点) |
40.132 7 | 74.361 0 | 0.241 9 | 1.909 7 | |
5 结论
本文提出一种适用于脉冲星导航的改进线性观测方程。该方程包含了脉冲星观测模型中会对观测量造成较大高阶截断误差的周年视差效应和引力延迟效应。仿真结果表明:
1) 在考虑周年视差效应和引力延迟效应影响的情况下, 简化的线性观测方程在EKF算法中会产生发散, 而改进的线性观测方程却能较快的得到250 m和2 m/s以内的位置和速度估计误差。
2) 改进的线性观测方程对高阶项的变化具有一定的鲁棒性。在高阶项变化的最大值和最小值点位置估计误差仅相差2.6 m, 速度估计误差仅相差0.038 2 m/s, 适合长时在轨航天器的导航解算。
| [1] | PINES D J. ARPA/DARPA space programs[M]. Arlington: XNAV Industry Day, 2004: 1-5. |
| [2] |
帅平, 张新源, 黄良伟, 等. 脉冲星导航试验卫星科学观测数据分析[J].
空间控制技术与应用, 2017, 43 (2): 1–6.
SHUAI P, ZHANG X Y, HUANG L W, et al. X-ray pulsar navigation test satellite science data analysis[J]. Aerospace Control and Application, 2017, 43 (2): 1–6. DOI:10.3969/j.issn.1674-1579.2017.02.001 (in Chinese) |
| [3] |
帅平, 陈绍龙, 吴一凡, 等. X射线脉冲星导航原理[J].
宇航学报, 2007, 28 (6): 1538–1543.
SHUAI P, CHEN S L, WU Y F, et al. Navigation principles using X-ray pulsars[J]. Journal of Astronautics, 2007, 28 (6): 1538–1543. DOI:10.3321/j.issn:1000-1328.2007.06.020 (in Chinese) |
| [4] |
毛悦, 宋小勇. 脉冲星时间模型精化及延迟修正分析[J].
武汉大学学报(信息科学版), 2009, 34 (5): 581–584.
MAO Y, SONG X Y. Accurating and delay correction analysis of pulsar timing model[J]. Geomatics and Information Science of Wuhan University, 2009, 34 (5): 581–584. (in Chinese) |
| [5] |
刘劲.基于X射线脉冲星的航天器自主导航方法研究[D].武汉: 华中科技大学, 2008: 55-56.
LIU J.X-ray pulsar-based spacecraft autonomous navigation[D].Wuhan: Huazhong University of Science and Technology, 2008: 55-56(in Chinese). |
| [6] |
乔黎.X射线脉冲星高轨道卫星自主导航及其应用技术研究[D].南京: 南京航空航天大学, 2010: 23-27.
QIAO L.X-ray pulsar-based autonomous navigation and its application to high earth orbits satellites[D].Nanjing: Nanjing University of Aeronautics and Astronautics, 2010: 23-27(in Chinese). |
| [7] |
王敏.基于X射线脉冲星的航天器自主导航滤波算法研究[D].哈尔滨: 哈尔滨工业大学, 2015: 14-18.
WANG M.The research of spacecraft autonomous navigation algorithm based on the X-ray pulsars[D].Harbin: Harbin Institute of Technology, 2015: 14-18(in Chinese). |
| [8] | CHEN P T, SPEYER J L, BAYARD D S, et al. Autonomous navigation using X-ray pulsars and multirate processing[J]. Journal of Guidance, Control, and Dynamics, 2017, 40 (9): 2237–2249. DOI:10.2514/1.G002705 |
| [9] | NING X L, GUI M Z, ZHANG J, et al. Impact of the pulsar's direction on CNS/XNAV integrated navigation[J]. Transactions on Aerospace and Electronic Systems, 2017, 53 (6): 3043–3055. DOI:10.1109/TAES.2017.2725518 |
| [10] |
宁晓琳, 马辛, 张学亮, 等. 基于ASUKF的火星探测器脉冲星自主导航方法[J].
北京航空航天大学学报, 2012, 38 (1): 22–27.
NING X L, MA X, ZHANG X L, et al. Autonomous pulsars navigation method based on ASUKF for Mars probe[J]. Journal of Beijing University of Aeronautics and Astronautics, 2012, 38 (1): 22–27. (in Chinese) |
| [11] |
费保俊, 孙维瑾, 潘高田, 等. X射线脉冲星自主导航的光子到达时间转换[J].
空间科学学报, 2010, 30 (1): 85–90.
FEI B J, SUN W J, PAN G T, et al. Transformation of photon time of arrival in XNAV[J]. Chinese Journal of Space Science, 2010, 30 (1): 85–90. (in Chinese) |
| [12] |
孙海峰.X射线脉冲星导航信号特性分析及具有多物理特性的仿真系统研究[D].西安: 西安电子科技大学, 2015: 25-28
SUN H F.Study on signal characteristics of X-ray pulsar based navigation and simulation experiment system with multi-physical features[D].Xi'an: Xidian University, 2015: 25-28(in Chinese). |
| [13] |
李小平, 方海燕, 孙海峰, 等. X射线脉冲星大尺度时间转换模型研究[J].
载人航天, 2015, 21 (6): 628–634.
LI X P, FANG H Y, SUN H F, et al. Research on large-scale time transform model of X-ray pulsar[J]. Manned Spaceflight, 2015, 21 (6): 628–634. DOI:10.3969/j.issn.1674-5825.2015.06.016 (in Chinese) |
| [14] | SHEIKH S I, HELLINGS R W, MATZNER R A.High-order pulsar timing for navigation[C]//Proceeding of the 63rd Annual Meeting of the Institute of Navigation.Cambridge: Institute of Navigation, 2007: 432-443. |
| [15] |
杨廷高. 关于脉冲星脉冲到达时间转换方程[J].
时间频率学报, 2009, 32 (2): 154–159.
YANG T G. On transfer equation of pulsar pulse arrival time[J]. Journal of Time and Frequency, 2009, 32 (2): 154–159. DOI:10.3969/j.issn.1674-0637.2009.02.012 (in Chinese) |
| [16] |
李志豪.基于X射线脉冲星的航天器导航滤波算法仿真分析[D].长沙: 国防科学技术大学, 2008: 52.
LI Z H.Simulation and analysis of the filters for spacecraft navigation based on X-ray pulsars[D].Changsha: National University of Defense Technology, 2008: 52(in Chinese). |
| [17] |
李晓宇, 姜宇, 金晶, 等. 脉冲星导航系统的星历表误差RKF校正算法[J].
宇航学报, 2017, 38 (1): 26–33.
LI X Y, JIANG Y, JIN J, et al. RKF method for pulsar based navigation with emphasis error correction[J]. Journal of Astronautics, 2017, 38 (1): 26–33. DOI:10.3873/j.issn.1000-1328.2017.01.004 (in Chinese) |