为此，本文以利用信标卫星进行脉冲星方位误差估计为背景，设计了两步卡尔曼滤波(Two-step Kalman filter，TSKF)算法解决了同时存在钟差和卫星位置误差影响的问题，从而提高了脉冲星方位误差估计的精度。
1 脉冲星方位误差估计的传统模型 利用卫星估计脉冲星方位误差是将卫星上X射线探测器测量的脉冲到达时间转换为到达太阳系质心坐标系(SSB)原点的时间t′SSB，并用该时间减去脉冲星相位时间模型计算的真实到达SSB原点的时间tSSB，而脉冲星方位误差就反映在两者之差中[3]，如图 1所示。
 |
| 图 1 脉冲星方位误差估计原理 Fig. 1 Principle of pulsar position error estimation |
| 图选项 |
转换过程如下[11]：
 | (1) |
式中：t′sat为信标卫星时钟测得的脉冲到达时间；n′为测量得到的带误差的单位方向矢量；rsat为卫星在SSB中的位置矢量；c为光速；o(t)为周年视差效应等引起的高阶项。由相位时间模型得到脉冲到达SSB的时间tSSB为
 | (2) |
式中：tsat为脉冲到达卫星的真实时间；n为脉冲星真实的单位方向矢量。设脉冲星的赤经为α，赤纬为β，则满足：
 | (3) |
设(Δα, Δβ)为脉冲星方位误差，则脉冲星方位信息的真实值与带误差的值之间满足：
 | (4) |
将式(4)代入式(3)，忽略二阶小项得
 | (5) |
脉冲星带误差的单位方向矢量为
 | (6) |
记单位方向矢量误差Δn为
 | (7) |
则可以将式(5)表示为
 | (8) |
若不考虑时钟钟差，则t′sat≈tsat，可得
 | (9) |
取状态变量为X=[Δα?? Δβ]T，则脉冲星方位误差估计的传统模型为
 | (10) |
 | (11) |
式中：Wk和ηk分别为系统噪声和量测噪声，k表示第k时刻; 状态转移矩阵
 |
rsat/x, rsat/y, rsat/z为在3个坐标轴方向的投影。
2 误差影响分析 2.1 钟差的影响 卫星上的时钟频率和相位会发生漂移，从而会造成脉冲到达卫星的真实时间和卫星时钟测得的脉冲到达时间存在偏差，这种偏差就是时钟钟差，设其为δt，则满足：
 | (12) |
此时
 | (13) |
卫星钟差模型为[12]
 | (14) |
式中：x1、x2分别为钟差、钟差漂移率；τ为步长；ω为白噪声，其方差为
 | (15) |
式中：q1、q2为噪声功率谱密度。其离散过程模型为
 | (16) |
式中: 状态转移矩阵为
取时钟误差漂移率为3.637 979×10－12，根据铷原子钟模型，取时钟噪声功率谱密度分别为q1=1.11×10－22 s和q2=2.22×10－32/s[9]。给定时钟初始时刻的钟差为0 s，取步长为1 s，则可得到钟差随时间的变化如图 2所示。
 |
| 图 2 钟差随时间变化 Fig. 2 Clock error changes with time |
| 图选项 |
由图 2可知，随着时间的推移，钟差达到5×10－6 s，将其与光速相乘理论上将造成1 500 m的误差，将会给系统造成不容忽略的影响。因此为验证其影响，仿真实验以脉冲星B0531+21作为观测脉冲星，其参数如表 1所示。
表 1 脉冲星B0531+21参数 Table 1 Parameters of pulsar B0531+21
| 参数 | 数值 |
| 赤经/(°) | 83.63 |
| 赤纬/(°) | 22.01 |
| 距离/kpc | 2.0 |
| P/ms | 33.4 |
| W/ms | 1.7 |
| Fx/(ph·cm－2·s－1) | 1.54 |
| pf/% | 70 |
| 注：kpc为秒差距，是天文学中使用的距离单位；ph·cm－2·s－1为宇宙背景噪声单位。 | |
表选项
其中，脉冲周期、宽度、辐射光子流量分别为P、W、Fx；pf为脉冲周期中辐射流量与平均值之比。脉冲星的观测噪声方差[13]可由式(17)计算得到：
 | (17) |
式中：设探测器有效面积A为1 m2；Bx=0.005 ph·cm－2·s－1为宇宙背景噪声；d为W与P之比；设观测时间tobs为1 000 s，则可计算得到[8]σ=(77.69 m)2。设方位误差初始值为(2, 2) mas(毫角秒)，并取初始状态为0。
表 2和图 3为卫星轨道参数和仿真结果。由图 3可知，在钟差的影响下赤经和赤纬误差估计结果均出现了较大的偏差，尤其是后者，甚至超过了40 mas，且估计结果都不能收敛到一个固定值。因此，为了使赤经和赤纬误差估计更加准确，有必要考虑修正钟差的影响。
表 2 卫星轨道参数 Table 2 Parameters of satellite orbit
| 参数 | 数值 |
| 半长轴/km | 7 460 |
| 离心率 | 4.55×10-16 |
| 轨道倾角/(°) | 25 |
| 近地点幅角/(°) | 45 |
| 升交点赤经/(°) | 0 |
| 初始真近地点/(°) | 30 |
| 起始时间 | 2015-07-01T12：00：00 |
表选项
 |
| 图 3 未修正钟差的估计结果 Fig. 3 Estimation result with uncorrected clock error |
| 图选项 |
2.2 卫星位置误差的影响 使用卫星估计方位误差的前提是卫星位置是精确已知的[3]，然而实际情况中，地面站得到的卫星位置会不可避免地存在偏差。因此当卫星位置存在误差Δr时，设rsat与偏差位置r′sat之间满足：
 | (18) |
忽略二阶小项，则观测模型变为
 | (19) |
此时若仍以式(9)作为观测模型，则会给系统引入一定的偏差。仿真时卫星轨道和脉冲星数据不变，设卫星位置误差在3个坐标方向的投影为100 m，估计结果如图 4所示。
 |
| 图 4 未修正卫星位置误差的估计结果 Fig. 4 Estimation results with uncorrected satellite position errors |
| 图选项 |
分析以上仿真结果可得，当信标卫星存在位置误差时，估计结果误差较大，且不收敛。可见，信标卫星位置误差也是不可忽略、需要被修正的。
2.3 钟差和卫星位置误差同时存在的影响 2.1节和2.2节是假设一种误差存在而另一种误差不存在而进行分析的，但是当利用信标卫星进行脉冲星方位误差估计时，时钟会发生漂移，同时卫星位置也不可避免地存在偏差。因此，当钟差和卫星位置误差同时存在时，可得系统新的观测模型为
 | (20) |
仿真条件不变，当2种误差都未修正时的仿真结果如图 5所示。
 |
| 图 5 2种误差都未修正的估计结果 Fig. 5 Estimation results with two types of error uncorrected |
| 图选项 |
可见，当同时存在钟差和卫星位置误差时，赤经和赤纬估计结果偏差较大，存在较大的发散。因此，为了得到更加精确的估计结果，要同时修正钟差和卫星位置误差的影响。
3 修正钟差和卫星位置误差的算法 3.1 增广状态模型 为了修正钟差的影响，对观测模型(20)，在传统模型状态变量为[Δα??Δβ]T的基础上，将钟差和钟差漂移率扩展为新的状态量。扩展后的状态量为X=[Δα??Δβ??x1??x2]T，进而得到新的状态和观测方程为
 | (21) |
 | (22) |
状态转移矩阵
 |
由于是线性时变系统，采用PWCS判据分析其可观性，将系统分为j个时间段，在各个时间段都将其视为线性定常系统[3]。
系统总的可观测性矩阵为
 | (23) |
第j个时间段的可观测矩阵为
 | (24) |
式中：
 | (25) |
其中任意相邻的3个时间段满足：
 | (26) |
对其进行初等行变换可得
 | (27) |
由于τ>0，且卫星运行时相邻3个位置一般不共面，可得
 | (28) |
 | (29) |
所以
 | (30) |
即系统满足可观测性条件。
3.2 TSKF算法 对于观测模型(22)，为修正卫星位置误差的影响，若仍采用增广状态法将卫星位置误差在3个坐标轴的分量增广为状态变量，则此时状态变量由4个增加到7个，不仅会增加计算量，还容易出现数值不稳定的问题[11]。因此为减小计算量并确保算法的可观测性，采用TSKF算法进行解算。
Friedland等[14-15]最先提出TSKF算法并将其用在处理系统常值偏差问题上，Ignagni[16]在此基础上将其应用到处理缓变偏差上。因此本文在新的状态和观测方程(21)和(22)基础上，根据两步卡尔曼滤波原理取第一步滤波的状态量为X=[Δα??Δβ??x1??x2]T，第二步滤波状态量b为卫星位置误差Δr，可将该算法的更新方程写为[17-18]
一步滤波方程
 | (31) |
二步滤波方程
 | (32) |
最终估计结果可写为
 | (33) |
式中：Qk为Wk的方差；Pk为状态量协方差；Vk-1为第二步状态变量对第一步状态变量的纠正矩阵；Rk为ηk的方差；Bk为Δr在式(21)中的驱动方程；Ck为Δr在式(22)中的驱动方程。由上述分析可得：Bk=0，Ck=n′。
仿真时，仍取脉冲星方位误差为(2，2)mas，时钟初始时刻的钟差为0 s，卫星轨道参数、脉冲星参数、钟差模型参数及其他参数不变，可得同时修正钟差和卫星位置误差的TSKF算法估计结果如图 6所示。
 |
| 图 6 TSKF算法估计结果 Fig. 6 Estimation results of TSKF algorithm |
| 图选项 |
为进一步验证TSKF算法的有效性，其他条件不变，设置不同的钟差和卫星位置误差，设置的条件和仿真结果分别如表 3和表 4所示。
表 3 条件设置 Table 3 Condition setup
| 钟差/(μs) | 卫星位置误差/km |
| 0 | (0.1, 0.1, 0.1) |
| 1 | (0, 0, 0) |
| 1 | (0.1，-0.1，0.1) |
| 2 | (0.1, 0.1, 0.1) |
| 2 | (0.2，-0.2, 0.2) |
| 2 | (0.5, 0.5, 0.5) |
| 5 | (0.1, 0.1, 0.1) |
| 5 | (0.2，-0.2, 0.2) |
表选项
表 4 仿真结果 Table 4 Simulation results
| 赤经估计偏差/mas | 赤纬估计偏差/mas |
| 0.050 | 0.043 |
| 0.061 | 0.047 |
| 0.150 | 0.065 |
| 0.150 | 0.072 |
| 0.172 | 0.110 |
| 0.221 | 0.150 |
| 0.182 | 0.143 |
表选项
由图 6可得，赤经误差估计结果收敛到2.2 mas以内，赤纬误差估计结果也在1.8 mas和2.1 mas之间，可见TSKF算法能使方位误差估计结果实现较好的收敛，使赤经和赤纬估计结果保持在约0.2 mas以内的精度。且由表 4可得，在不同的钟差和卫星位置误差下，TSKF算法均能使估计结果保持较高的精度。因此可得，TSKF算法能有效隔离钟差和卫星位置误差的影响，使赤经和赤纬误差估计结果保持在没有误差影响的水平下。
4 结论 1) 在利用信标卫星进行脉冲星方位误差估计时，时钟钟差和卫星位置误差是需要被修正的，否则会严重影响估计精度。
2) 在重新推得的包含钟差和卫星位置误差模型基础上，采用TSKF算法能有效克服钟差和卫星位置误差的影响，使脉冲星方位误差估计精度显著提高。
参考文献
| [1] | SHEIKH S I. The use of variable celestial X-ray sources for spacecraft navigation[D]. Maryland: University of Maryland, 2005. |
| [2] | 刘劲. 基于X射线脉冲星的航天器自主导航方法研究[D]. 武汉: 华中科技大学, 2011. LIU J. X-ray pulsar-based spacecraft autonomous navigation[D]. Wuhan: Huazhong University of Science and Technology, 2011(in Chinese). |
| [3] | 孙守明. 基于X射线脉冲星的航天器自主导航方法研究[D]. 长沙: 国防科学技术大学, 2011. SUN S M. Study on autonomous navigation method of spacecraft based on X-ray pulsars[D]. Changsha: National University of Defense Technology, 2011(in Chinese). |
| [4] | 孙守明, 郑伟, 汤国建. X射线脉冲星星表方位误差估计算法研究[J]. 飞行器测控学报, 2010, 29(2): 57-60. SUN S M, ZHENG W, TANG G J. A new estimation algorithm of the X-ray pulsar position error[J]. Journal of Spacecraft TT & C Technology, 2010, 29(2): 57-60. (in Chinese) |
| [5] | 熊凯, 魏春岭, 刘良栋. 鲁棒滤波技术在脉冲星导航中的应用[J]. 空间控制技术与应用, 2008, 34(6): 8-11. XIONG K, WEI C L, LIU L D. Application of robust filtering in pulsars-based navigation[J]. Aerospace Control and Application, 2008, 34(6): 8-11. (in Chinese) |
| [6] | XU Q, WANG H L, FENG L, et al. An improved augmented X-ray pulsar navigation algorithm based on the norm of pulsar direction error[J]. Advances in Space Research, 2018, 62(11): 3187-3198. DOI:10.1016/j.asr.2018.08.026 |
| [7] | NING X L, GUI M Z, FANG J C, et al. Differential X-ray pulsar aided celestial navigation for Mars exploration[J]. Aerospace Science and Technology, 2017, 62: 36-45. DOI:10.1016/j.ast.2016.10.032 |
| [8] | 王宏力, 许强, 由四海, 等. 考虑卫星位置误差的增广脉冲星方位误差估计算法[J]. 国防科技大学学报, 2018, 40(5): 177-182. WANG H L, XU Q, YOU S H, et al. Augmented estimation algorithm for pulsar position error with satellite position error[J]. Journal of National University of Defense Technology, 2018, 40(5): 177-182. (in Chinese) |
| [9] | 孙守明, 郑伟, 汤国建. X射线脉冲星/SINS组合导航中的钟差修正方法研究[J]. 国防科技大学学报, 2010, 32(6): 82-86. SUN S M, ZHENG W, TANG G J. A new clock error control algorithm of X-ray pulsars/SINS integrated navigation[J]. Journal of National University of Defense Technology, 2010, 32(6): 82-86. (in Chinese) |
| [10] | 王璐, 史晨曦, 李建勋, 等. 考虑钟差修正的脉冲星与多普勒差分组合导航[J]. 电波科学学报, 2018, 33(2): 217-224. WANG L, SHI C X, LI J X, et al. Pulsar and Doppler difference integrated navigation considering clock error control[J]. Chinese Journal of Radio Science, 2018, 33(2): 217-224. (in Chinese) |
| [11] | EMADZADEH A A, SPEYER J L. X-ray pulsar-based relative navigation using epoch folding[J]. IEEE Transactions on Aerospace and Electronic Systems, 2011, 47(4): 2317-2328. |
| [12] | 李孝辉, 吴海涛, 高海军, 等. 用Kalman滤波器对原子钟进行控制[J]. 控制理论与应用, 2003, 20(4): 551-554. LI X H, WU H T, GAO H J, et al. Clock disciplined method by using Kalman filter[J]. Control Theory & Applications, 2003, 20(4): 551-554. (in Chinese) |
| [13] | 帅平, 陈绍龙, 吴一帆, 等. X射线脉冲星导航技术研究进展[J]. 空间科学学报, 2007, 27(2): 169-176. SHUAI P, CHEN S L, WU Y F, et al. Advance in X-ray pulsar navigation technology[J]. Chinese Journal of Space Science, 2007, 27(2): 169-176. (in Chinese) |
| [14] | FRIEDLAND B. Treatment of bias in recursive filtering[J]. IEEE Transactions on Automatic Control, 1969, 14(4): 359-367. |
| [15] | HAESSIG D, FRIEDLAND B. Separate-bias estimation with reduced-order Kalman filters[J]. IEEE Transactions on Automatic Control, 1998, 43(7): 983-987. |
| [16] | IGNAGNI M B. Separate bias Kalman estimator with bias state noise[J]. IEEE Transactions on Automatic Control, 1990, 35(3): 338-341. |
| [17] | 王奕迪, 郑伟, 孙守明, 等. 考虑系统偏差的脉冲星守时算法研究[J]. 国防科技大学学报, 2013, 35(2): 12-16. WANG Y D, ZHENG W, SUN S M, et al. Algorithm for the pulsar timing system with the system bias[J]. Journal of National University of Defense Technology, 2013, 35(2): 12-16. (in Chinese) |
| [18] | 许强, 范小虎, 徐利国, 等. 脉冲星方位误差估计的TSKF算法[J]. 北京航空航天大学学报, 2020, 46(4): 761-768. XU Q, FAN X H, XU L G, et al. TSKF algorithm for pulsar position error estimation[J]. Journal of Beijing University of Aeronautics and Astronautics, 2020, 46(4): 761-768. (in Chinese) |
