脉冲星方位误差估计的TSKF算法

2020-06-01许强范小虎徐利国王宏力冯磊

北京航空航天大学学报 2020年4期

许强，范小虎，徐利国，王宏力，冯磊，*

（1.青州高新技术研究所测试控制系，潍坊262500； 2.火箭军工程大学导弹工程学院，西安710025）

X射线脉冲星导航作为一种新兴的天文自主导航，在深空探测、卫星授时等领域引起了研究人员的极大关注并掀起了广泛的研究热潮［1-3］。由于使用脉冲星发出的辐射信号作为信号源，X射线脉冲星导航相对于其他导航方式具有信号稳定、抗干扰性强及自主程度高等特点，能够较好满足太空卫星、深空探测器等航天器的任务需求。虽然目前相关国家针对X射线脉冲星导航的研究已经进入了实验验证阶段，但距离真正的工程应用仍然有一定的距离。其中一个较为重要的影响因素即为脉冲星存在的方位误差问题。

由于大气层对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 Propermotion distribution of pulsars

图2 脉冲星方位自行速度不确定度分布Fig.2 Propermotion uncertainty distribution of pulsars

通过图1和图2可以看出，大多数脉冲星的赤经、赤纬自行速度在50mas／a以内，不确定度基本在10 mas／a以内。如果在方位误差估计中不考虑方位自行速度的影响，势必会降低估计精度甚至产生发散。

以脉冲星B1821-24为例，其基本参数如表1所示。当使用文献［9］中的增广算法进行方位误差估计时，假设探测器面积为1 m2，观测周期为1 000 s，宇宙X 射线背景流量Bx＝0.005 ph／（cm2·s）（ph表示通过的光子个数），则观测噪声方差可计算得［6］σR＝（230.01m）2。使用同一卫星轨道，分别在有方位自行速度和无方位自行速度的条件下进行导航计算。2种条件下卫星的位置误差均为［100，100，100］m，脉冲星的方位误差为［2，2］mas。设定存在的脉冲星方位自行速度为［10，10］mas／a。具体仿真结果如图3所示。

通过对图3的分析发现，在无方位自行速度的情况下，增广算法可有效隔离卫星位置误差的影响，进而实现脉冲星方位误差的精准估计。但是存在方位自行速度时，估计算法会产生较大的发散，无法正常工作。

以上数值仿真说明，在对脉冲星方位误差估计时有必要考虑方位自行速度的影响。

表1 脉冲星B1821-24参数Tab le 1 Param eters of pulsar B1821-24

1.2 卫星位置误差的影响

孙守明等［10］在后续研究中提出了基于CV模型的X射线脉冲星方位误差估计算法。该算法将方位自行速度作为状态量单独估计，有效解决了方位自行速度的影响问题，但该算法并没有考虑到卫星位置误差带来的影响。

卫星位置误差对脉冲星方位估计的影响原理如图4所示。当卫星位置不存在误差时，根据脉冲星方位误差的观测模型，可以认为脉冲信号到达太阳系质心（Solar System Barycenter，SSB）处的时间延迟仅与脉冲星方位误差有关。假设n为真实方向，～n为带误差的脉冲星方向，其满足关系：

式中：Δn为方位误差。

图3 不同条件下增广算法仿真结果Fig.3 Simulation results of augmented algorithm under different conditions

则观测模型为［14］

式中：c为光速；tSSB为真实SSB处脉冲到达时间（Time-of-Arrival，TOA）；～tSSB为根据观测数据推算得到的SSB处到达时间；rsat为真实卫星位置。

但若卫星存在位置误差Δr，在误差影响下实际的观测模型就发生了改变。存在Δr时，推算得到的到达时间＾tSSB与实际的到达时间tSSB之差就不再是单纯的由脉冲星方位误差导致。假设卫星导航位置

＾rsat与真实位置rsat满足：

图4 卫星位置误差对估计算法的影响原理Fig.4 Principle of influence of satellite position error on estimation algorithm

则此时观测模型应当变为

如果此时仍然以式（2）作为观测模型，不仅会引入一定的系统偏差，还有可能在地球自转的影响下产生发散。

同样以脉冲星B1821-24为例，假设脉冲星方位误差为［2，2］mas，方位自行速度为［10，10］mas／a。其他条件不变，采用基于CV模型的X射线脉冲星方位误差估计算法分别在无卫星位置误差和有卫星位置误差情况下进行导航解算，卫星位置误差设为［100，100，100］m，具体仿真结果如图5所示。

通过对图5的分析可见，在加入卫星位置误差之前，基于CV模型的脉冲星方位误差估计算法可以较为精准地估计出当前的方位误差和方位自行速度。但是在引入卫星位置误差之后，由于地球的自转，估计结果无法收敛在某一固定值。结合文献［9］的分析，说明无论是否存在脉冲星方位自行速度，信标卫星的位置误差都应当成为估计算法重点解决的工程问题之一。

图5 不同条件下基于CV模型的估计算法仿真结果Fig.5 Simulation results of estimation algorithm based on CV model under different conditions

2 TSKF算法

TSKF算法最早由Friedland［15］提出，用于解决线性系统中的定常偏差问题。Hsieh 和Chen［16］将两级滤波思想用于标准卡尔曼滤波算法，证明最高可以将计算量降低59%。

本文在基于CV模型的方位误差估计算法基础上，采用两级滤波的方法，将脉冲星方位误差和方位自行速度作为第一级滤波状态量，卫星位置误差作为第二级滤波的状态量，在不增加状态维数的前提下实现同步估计，有效隔离卫星位置误差对估计算法的影响。

结合CV模型和第1节的分析，可将TSKF算法的离散空间状态方程写为

式中：Xk＝［ΔαkΔ˙αkΔδkΔ˙δk］T为第一级滤波的状态量，分别代表赤经、赤经自行速度、赤纬、赤纬自行速度；T为计算步长；Ak－1为状态转移矩阵；Wk为系统噪声；Zk为观测量；Δt为脉冲到达航天器与SSB的时间差；ηk为观测噪声；Hk为观测矩阵且满足：

第二级滤波状态更新：

分析以上过程还可以发现，TSKF算法的第一级滤波与常规的卡尔曼滤波算法完全相同，只是第二级滤波与一般滤波过程不同。实际上，在第二级滤波过程中公式Mk的作用是描述状态量bk的估计方差，其对第二级滤波增益¯Kk的计算具有重要的影响作用。因此，式（14）相当于第一级滤波中的式（8），是估计方差阵的更新方程。在得到第二级滤波的增益矩阵¯Kk后，式（17）便相当于状态量的更新估计，与第1节滤波中的式（10）作用相同，其中（I－¯KkSk）bk－1为状态量bk的一步预测，相当于第一级滤波中的式（7）。

而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。

其他条件不变，将TSKF算法在不同卫星位置误差及方位自行速度条件下分别运行50次，其中每次运行的结果取为最后一天所有计算值的平均值。将每个算法运行50次的结果再取平均值作为此时该算法的最终精度。具体条件及误差统计如表2和表3所示。

表2 仿真条件设置Tab le 2 Sim u lation condition setup

通过分析图7、表2和表3可见，在不同误差条件下，TSKF算法均可较好地完成收敛，并实现方位误差最大约0.1mas及方位自行速度估计误差最大约1.1mas／a的精度。同时，其收敛速度也明显快于第1节提到的增广算法和基于CV模型的估计算法。

考虑到实际卫星在轨运行时，位置误差可能根据轨道周期变化，故将卫星位置误差设置为随卫星轨道呈三角函数变化的形式。为体现普遍性，其具体关系式满足：

表3 仿真结果统计Tab le 3 Sim ulation result statistics

式中：Ts为卫星轨道周期；t为卫星运行时间；Lx、Ly、Lz为对应的幅值。

假设脉冲星方位误差及方位自行速度同样为［2，2］mas和［10，10］mas／a，则当Lx、Ly、Lz均为100m时，TSKF算法的具体运行过程如图8所示。

采用同样的统计方法，将不同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模型的基础上继续采用状态增广的方法带来的计算负担小。