由于大气层对X射线信号的衰减作用，地面只能收到脉冲星在射电频段的极其微弱信号，所以为了得到较为准确的观测数据，需要使用甚长基线干涉测量(Very Long Baseline Interferometry, VLBI)技术。这就导致目前对脉冲星的观测不仅成本高，而且不同频段的观测精度只能保证在毫角秒(milliarcsecond, mas)量级[4-5]。但对于脉冲星导航而言，几毫角秒的方位误差都有可能带来数百米的导航偏差。为此，国内许多****从组合导航[6]、鲁棒滤波算法[7]等方面开展了相关研究。同时，孙守明等[8]也提出了基于信标卫星的脉冲星方位误差估计算法，以提高脉冲星的方位精度。文献[9]采用增广扩展卡尔曼滤波的方法对该算法进行了改进，考虑了实际存在的信标卫星位置误差的影响。但是，以上2种算法都未考虑脉冲星方位的自行速度。虽然孙守明等[10]在后续的研究中采用匀速(Constant Velocity, CV)模型，在算法中加入了脉冲星的方位自行速度，但仍未消除信标卫星位置误差的干扰。然而这部分干扰是无法避免且较为致命的。
若在CV模型基础上继续采用文献[9]中的增广算法解决信标卫星位置误差的影响，将会使得矩阵运算由4维变为7维。当对多个脉冲星同时观测时会大幅增加运算负担，还容易出现数值病态的问题。为此，本文设计了两级卡尔曼滤波(Two-Stage Kalman Filter, TSKF)算法。该算法可以在不增加维数的情况下实现对系统偏差的在线估计，兼顾卫星位置误差和方位自行速度影响的同时，实现较高精度的方位误差估计。
1 误差影响 1.1 方位自行速度的影响 对于宇宙中的绝大多数脉冲星而言，其方位并非一成不变，而是存在一定自行速度的。研究人员认为产生自行的因素很可能是超新星爆发不是各向同性的[11]。这种自行短时间内受观测技术限制难以准确测出[12]，但长时间看又存在跟随脉冲星自身运动产生突变的可能性，影响脉冲星导航的使用。
以澳大利亚国家天文台(Australia Telescope National Facility，ATNF)提供的脉冲星数据库为例，所有2 702颗脉冲星中，已知自行数据的有306颗[13]，且均具有一定的不确定度。其方位自行速度及其满足1个σ标准差分布时不确定度情况分别如图 1和图 2所示。
 |
| 图 1 脉冲星方位自行速度分布 Fig. 1 Proper motion distribution of pulsars |
| 图选项 |
 |
| 图 2 脉冲星方位自行速度不确定度分布 Fig. 2 Proper motion uncertainty distribution of pulsars |
| 图选项 |
通过图 1和图 2可以看出，大多数脉冲星的赤经、赤纬自行速度在50 mas/a以内，不确定度基本在10 mas/a以内。如果在方位误差估计中不考虑方位自行速度的影响，势必会降低估计精度甚至产生发散。
以脉冲星B1821-24为例，其基本参数如表 1所示。当使用文献[9]中的增广算法进行方位误差估计时，假设探测器面积为1 m2，观测周期为1 000 s，宇宙X射线背景流量Bx=0.005 ph/(cm2·s)(ph表示通过的光子个数)，则观测噪声方差可计算得[6]σR=(230.01 m)2。使用同一卫星轨道，分别在有方位自行速度和无方位自行速度的条件下进行导航计算。2种条件下卫星的位置误差均为[100,100, 100]m，脉冲星的方位误差为[2, 2]mas。设定存在的脉冲星方位自行速度为[10, 10]mas/a。具体仿真结果如图 3所示。
表 1 脉冲星B1821-24参数 Table 1 Parameters of pulsar B1821-24
| 参数 | 数值 |
| 赤经/(°) | 276.13 |
| 赤纬/(°) | -24.87 |
| 距离/(1020m) | 1.694 |
| 脉冲周期P/(10-3s) | 3.045 |
| 脉冲宽度W/(10-5s) | 5.5 |
| 光子流量Fx/(10-4ph·(cm2·s)-1) | 1.93 |
| 脉冲辐射流量与平均辐射流量之比Pf/% | 98 |
表选项
 |
| 图 3 不同条件下增广算法仿真结果 Fig. 3 Simulation results of augmented algorithm under different conditions |
| 图选项 |
通过对图 3的分析发现，在无方位自行速度的情况下，增广算法可有效隔离卫星位置误差的影响，进而实现脉冲星方位误差的精准估计。但是存在方位自行速度时，估计算法会产生较大的发散，无法正常工作。
以上数值仿真说明，在对脉冲星方位误差估计时有必要考虑方位自行速度的影响。
1.2 卫星位置误差的影响 孙守明等[10]在后续研究中提出了基于CV模型的X射线脉冲星方位误差估计算法。该算法将方位自行速度作为状态量单独估计，有效解决了方位自行速度的影响问题，但该算法并没有考虑到卫星位置误差带来的影响。
卫星位置误差对脉冲星方位估计的影响原理如图 4所示。当卫星位置不存在误差时，根据脉冲星方位误差的观测模型，可以认为脉冲信号到达太阳系质心(Solar System Barycenter, SSB)处的时间延迟仅与脉冲星方位误差有关。假设n为真实方向，
 |
| 图 4 卫星位置误差对估计算法的影响原理 Fig. 4 Principle of influence of satellite position error on estimation algorithm |
| 图选项 |
 | (1) |
式中:Δn为方位误差。
则观测模型为[14]
 | (2) |
式中：c为光速；tSSB为真实SSB处脉冲到达时间(Time-of-Arrival, TOA)；
但若卫星存在位置误差Δr，在误差影响下实际的观测模型就发生了改变。存在Δr时，推算得到的到达时间SSB与实际的到达时间
 | (3) |
则此时观测模型应当变为
 | (4) |
如果此时仍然以式(2)作为观测模型，不仅会引入一定的系统偏差，还有可能在地球自转的影响下产生发散。
同样以脉冲星B1821-24为例，假设脉冲星方位误差为[2, 2]mas，方位自行速度为[10, 10]mas/a。其他条件不变，采用基于CV模型的X射线脉冲星方位误差估计算法分别在无卫星位置误差和有卫星位置误差情况下进行导航解算，卫星位置误差设为[100, 100, 100]m，具体仿真结果如图 5所示。
 |
| 图 5 不同条件下基于CV模型的估计算法仿真结果 Fig. 5 Simulation results of estimation algorithm based on CV model under different conditions |
| 图选项 |
通过对图 5的分析可见，在加入卫星位置误差之前，基于CV模型的脉冲星方位误差估计算法可以较为精准地估计出当前的方位误差和方位自行速度。但是在引入卫星位置误差之后，由于地球的自转，估计结果无法收敛在某一固定值。结合文献[9]的分析，说明无论是否存在脉冲星方位自行速度，信标卫星的位置误差都应当成为估计算法重点解决的工程问题之一。
2 TSKF算法 TSKF算法最早由Friedland[15]提出，用于解决线性系统中的定常偏差问题。Hsieh和Chen[16]将两级滤波思想用于标准卡尔曼滤波算法，证明最高可以将计算量降低59%。
本文在基于CV模型的方位误差估计算法基础上，采用两级滤波的方法，将脉冲星方位误差和方位自行速度作为第一级滤波状态量，卫星位置误差作为第二级滤波的状态量，在不增加状态维数的前提下实现同步估计，有效隔离卫星位置误差对估计算法的影响。
结合CV模型和第1节的分析，可将TSKF算法的离散空间状态方程写为
 | (5) |
 | (6) |
式中：
 |
其中：
分析以上模型，可见此时算法中出现的常值偏差仅存在于观测方程，与状态方程无关。取第二级滤波状态量b为卫星位置误差Δr，结合算法流程，可将方位误差估计的TSKF算法更新方程表达如下。
第一级滤波时间更新：
 | (7) |
 | (8) |
第一级滤波状态更新：
 | (9) |
 | (10) |
 | (11) |
第二级滤波时间更新：
 | (12) |
 | (13) |
 | (14) |
第二级滤波状态更新：
 | (15) |
 | (16) |
 | (17) |
最终估计结果为
 | (18) |
式中：Qk为系统噪声Wk的方差；Pk为状态量的协方差；Pk k-1 为状态量的一步预测协方差；Rk为观测噪声ηk的方差；Kk为第一级滤波器的滤波增益；Vk为第二级滤波状态量bk对第一级滤波状态量Xk的纠正矩阵；Sk和Uk为迭代计算过程的中间量；Mk为第二级滤波中状态量的协方差；Kk为第二级滤波增益矩阵；Bk-1为Δr在状态方程式(5)中的驱动矩阵；Ck为Δr在观测方程式(6)中的驱动方程；I为单位阵。结合前述分析，已知
结合式(7)~式(18)可得，算法运行过程中除最后结果的整合外，第二级滤波仅用到第一级的增益值Kk和残差Zk－HkXk k－1 。故可将两级滤波并行计算以提高运算效率，其并行计算的流程如图 6所示。
 |
| 图 6 并行计算流程图 Fig. 6 Parallel computing flowchart |
| 图选项 |
分析以上过程还可以发现，TSKF算法的第一级滤波与常规的卡尔曼滤波算法完全相同，只是第二级滤波与一般滤波过程不同。实际上，在第二级滤波过程中公式Mk的作用是描述状态量bk的估计方差，其对第二级滤波增益Kk的计算具有重要的影响作用。因此，式(14)相当于第一级滤波中的式(8)，是估计方差阵的更新方程。在得到第二级滤波的增益矩阵Kk后，式(17)便相当于状态量的更新估计，与第1节滤波中的式(10)作用相同，其中
而Vk的作用是纠正常值偏差bk对第一级滤波估计值Xk的传递影响作用，故可将其命名为纠正矩阵。由于第一级滤波估计中完全没有涉及到常值偏差bk的计算，所以得到的估计值Xk必然是带有一定误差的。这部分误差会随着时间的推移而不断变化。为此，在第二级滤波估计中，不仅要估计出常值偏差bk的值，还要利用式(15)实时求解当前的纠正矩阵。最后，通过式(18)将纠正矩阵Vk与第二级滤波估计值bk的乘积加到第一级滤波的估计结果中便实现了常值偏差与第一级滤波状态量之间的隔离。第二级滤波时间更新环节中对矩阵Uk和Sk的计算均为计算纠正矩阵Vk的中间过程。
3 仿真分析 为证明TSKF算法的有效性，在方位自行速度及方位误差都存在的情况下进行仿真验证。所选用的脉冲星及其他相关参数与第1节相同。脉冲星方位误差为[2, 2]mas，方位自行速度为[10, 10]mas/a，卫星位置误差为[100, 100, 100]m。
系统噪声的方差Qk=diag[q12 q22 q12 q22], 其中q1=10－9arc sec≈2.78×10－13(°)，q2=10－12arc sec/s≈2.78×10－16(°)/s。M0=diag[m02 m02 m02]，其中m0=0.1 km。方位自行速度噪声的标准差为[10-18, 10-18](°)/s。具体仿真过程如图 7所示。
 |
| 图 7 TSKF算法仿真过程 Fig. 7 Simulation process of TSKF algorithm |
| 图选项 |
其他条件不变，将TSKF算法在不同卫星位置误差及方位自行速度条件下分别运行50次，其中每次运行的结果取为最后一天所有计算值的平均值。将每个算法运行50次的结果再取平均值作为此时该算法的最终精度。具体条件及误差统计如表 2和表 3所示。
表 2 仿真条件设置 Table 2 Simulation condition setup
| 实验编号 | 方位自行速度/ (mas·a-1) | 卫星位置误差/m |
| 1 | (0, 0) | (100, 100, 100) |
| 2 | (10, 10) | (0, 0, 0) |
| 3 | (-10, 10) | (100, 100, 100) |
| 4 | (20, 20) | (100, 100, 100) |
| 5 | (10, 10) | (100, -100, 100) |
| 6 | (10, 10) | (500, 500, 500) |
表选项
表 3 仿真结果统计 Table 3 Simulation result statistics
| 实验编号 | 方位误差/mas | 速度误差/(mas·a-1) | |||
| 赤经 | 赤纬 | 赤经 | 赤纬 | ||
| 1 | -0.003 9 | -0.005 8 | 0.016 9 | 0.201 2 | |
| 2 | 0.004 9 | -0.024 6 | -0.003 6 | -0.542 9 | |
| 3 | 0.024 1 | -0.106 4 | 0.254 6 | -1.013 5 | |
| 4 | 0.004 9 | 0.014 3 | -0.159 6 | 0.116 7 | |
| 5 | -0.005 6 | 0.092 9 | -0.103 4 | 1.081 4 | |
| 6 | -0.001 6 | -0.113 8 | 0.020 8 | 0.142 7 | |
表选项
通过分析图 7、表 2和表 3可见，在不同误差条件下，TSKF算法均可较好地完成收敛，并实现方位误差最大约0.1 mas及方位自行速度估计误差最大约1.1 mas/a的精度。同时，其收敛速度也明显快于第1节提到的增广算法和基于CV模型的估计算法。
考虑到实际卫星在轨运行时，位置误差可能根据轨道周期变化，故将卫星位置误差设置为随卫星轨道呈三角函数变化的形式。为体现普遍性，其具体关系式满足：
 | (19) |
式中：Ts为卫星轨道周期；t为卫星运行时间；Lx、Ly、Lz为对应的幅值。
假设脉冲星方位误差及方位自行速度同样为[2, 2] mas和[10, 10]mas/a，则当Lx、Ly、Lz均为100 m时，TSKF算法的具体运行过程如图 8所示。
 |
| 图 8 卫星位置误差周期变化时TSKF算法仿真结果 Fig. 8 Simulation results of TSKF algorithm when satellite position error period changes |
| 图选项 |
采用同样的统计方法，将不同Lx、Ly、Lz取值时的仿真结果统计如表 4所示。
通过分析图 8可得，当卫星位置误差出现周期性的变化时只会导致曲线的“毛刺”愈加明显。这是因为算法在达到稳态后，位置误差的周期变化相当于系统噪声有所增加，所以TSKF算法的估计结果会出现“毛刺”增加现象。但从表 4的结果来看，这对估计结果的影响非常的小。这是因为本文将一段时间内估计结果的平均值作为最终的结果，消除了“毛刺”的影响。因此，若采取本文类似处理措施或在算法中添加相应的平滑处理环节，那么便可消除位置误差周期变化的影响。
最后，为证明TSKF算法的高效性，将其与文献[10]中4状态量的基于CV模型的方位误差估计算法进行对比。将仿真运行时间设为一个自然年，分别统计2个算法MATLAB程序中除参数初始化及观测数据模拟部分的浮点运算次数。其中TSKF算法为30 045 202 038次，基于CV模型的估计算法为30 030 750 327次，前者仅比后者增加了0.048%，显然比在CV模型的基础上继续采用状态增广的方法带来的计算负担小。
表 4 卫星位置误差周期变化时仿真结果统计 Table 4 Simulation result statistics when satellite position error period changes
| 实验条件/m | 方位误差/mas | 速度误差/(mas·a-1) | ||||||
| Lx | Ly | Lz | 赤经 | 赤纬 | 赤经 | 赤纬 | ||
| 100 | 100 | 100 | 0.015 1 | -0.077 6 | -0.085 6 | 0.153 0 | ||
| 200 | 200 | 200 | -0.016 6 | -0.121 9 | -0.025 2 | 0.059 7 | ||
| 300 | 300 | 300 | -0.020 5 | -0.132 4 | -0.047 3 | 0.103 7 | ||
| 100 | 200 | 300 | -0.008 3 | -0.138 5 | -0.059 4 | 0.126 7 | ||
| 300 | 200 | 100 | 0.015 4 | -0.075 5 | -0.051 4 | 0.092 8 | ||
表选项
4 结论 1) 在脉冲星方位误差估计中，方位自行速度和卫星位置误差均有必要考虑在内，否则会严重影响算法的精度。
2) TSKF算法可以在方位自行速度和卫星位置误差都存在的情况下正常工作，并在仿真实验中基本达到了0.1 mas的方位估计精度和1.1 mas/a的方位自行速度估计精度。
3) TSKF算法的最大状态量维数与基于CV模型的估计算法相同，且浮点运算数也仅增加了0.048%。相对于增广的方法而言，两级滤波的方法计算效率更高，更不容易出现数值病态等问题。
参考文献
| [1] | GRAVEN P, COLLINS J, SHEIKH S, et al. XNAV for deep space navigation[J]. Advances in the Astronautical Sciences, 2008, 131: 349-364. |
| [2] | SHEIKH S I, HELLINGS R W, MATZNER R A.High-order pulsar timing for navigation[C]//Proceedings of Annual Meeting of the Institute of Navigation.Manassas: Institute of Navigation, 2007: 432-443. |
| [3] | YOU S, WANG H, HE Y, et al. Pulsar profile construction based on double-redundant-dictionary and same-scale L1-norm compressed sensing[J]. Optik, 2018, 164: 617-623. DOI:10.1016/j.ijleo.2018.03.066 |
| [4] | TAO A, HONG X Y, ZHENG W M, et al. Space very long baseline interferometry in China[J]. Advances in Space Research, 2020, 65(2): 850-855. DOI:10.1016/j.asr.2019.03.030 |
| [5] | WAJIMA K, HAGIWARA Y, TAO A, et al. The East-Asian VLBI network[J]. Physics, 2015, 502: 403-413. |
| [6] | NING X, GUI M, FANG J, et al. Differential X-ray pulsar aided celestial navigation for Mars exploration[J]. Aerospace Science & Technology, 2017, 62: 36-45. |
| [7] | XU Q, WANG H, 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 |
| [8] | 孙守明, 郑伟, 汤国建. 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) |
| [9] | 王宏力, 许强, 由四海, 等. 考虑卫星位置误差的增广脉冲星方位误差估计算法[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) |
| [10] | 孙守明, 郑伟, 汤国建. 基于CV模型的X射线脉冲星位置误差估计[J]. 系统仿真学报, 2010, 22(11): 2712-2714. SUN S M, ZHENG W, TANG G J. Position error estimation of X-ray pulsar based on CV model[J]. Journal of System Simulation, 2010, 22(11): 2712-2714. (in Chinese) |
| [11] | 孙守明.基于X射线脉冲星的航天器自主导航方法研究[D].长沙: 国防科学技术大学, 2011. SUN S M.Study on autonomous navigation method of spacecraft based on X-ray pulsar[D].Changsha: National University of Defense Technology, 2011(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-90002-2007140905.htm |
| [12] | BOBOLTZ D A, FEY A L, JOHNSTON K J, et al. Astrometric positions and proper motions of 19 radio stars[J]. Astronomical Journal, 2003, 126(1): 484-493. DOI:10.1086/375462 |
| [13] | MANCHESTER R N, HOBBS G B, TEOH A, et al. The australia telescope national facility pulsar catalogue[J]. Astronomical Journal, 2005, 129(4): 1993-2006. DOI:10.1086/428488 |
| [14] | 许强, 王宏力, 何贻洋, 等. 基于截断误差的改进脉冲星导航观测方程[J]. 北京航空航天大学学报, 2018, 44(9): 1974-1981. XU Q, WANG H L, HE Y Y, et al. Improved pulsar navigation measurement equation based on truncation errors[J]. Journal of Beijing University of Aeronautics and Astronautics, 2018, 44(9): 1974-1981. (in Chinese) |
| [15] | FRIEDLAND B. Treatment of bias in recursive filtering[J]. IEEE Transactions on Automatic Control, 1969, AC-14(4): 359-367. |
| [16] | HSIEH C S, CHEN F C. General two-stage Kalman filters[J]. IEEE Transactions on Automatic Control, 2000, 45(4): 819-824. |
