丁陶伟，帅平，黄良伟，张新源

中国空间技术研究院 钱学森空间技术实验室，北京 100094

X射线脉冲星属于高速旋转的中子星，具有极稳定的周期性。利用X射线脉冲星导航能够为航天器提供包括时间、位置、速度和姿态在内的全面的导航参考信息，为解决近地轨道和深空探测航天器的持续高精度自主导航难题提供了一种新思路。近十余年来，利用X射线脉冲星的航天器自主导航一直是国内外航天前沿技术研究的热点领域[1-2]。2004年，美国国防部国防预先研究计划局(DARPA)提出“X射线导航与自主定位验证(XNAV)”计划[3]，其目标是建立一个脉冲星网络，利用脉冲星辐射的X射线信号实现航天器长时间自主导航。自2005年以来，中国持续开展X射线脉冲星导航系统理论与方法研究，研制了多种类型的X射线探测器产品、脉冲星导航地面试验系统和脉冲星导航专用试验卫星[4]。2016年11月10日，中国成功发射全球首颗脉冲星导航专用试验卫星——X射线脉冲星导航1号(XPNAV-1)，率先开展脉冲星导航空间飞行试验[5-6]。2017年6月3日，美国宇航局(NASA)将“空间站X射线计时与导航技术试验(SEXTANT)”设备发射到国际空间站(ISS)上，开展脉冲星导航在轨演示验证试验[7]。该试验项目是“中子星内部结构探测器(NICER)”任务之一，其X射线计时仪的有效探测面积大于1 800 cm2，探测能谱范围为0.2～1.2 keV。

在利用X射线脉冲星的航天器自主轨道确定算法研究方面，前人已经进行了大量研究工作，主要包括针对模型不确定性的鲁棒滤波算法[8]、为降低非线性影响的非线性预测滤波(NPF)算法[9]、针对脉冲星方向误差的扩展状态无迹卡尔曼滤波(ASUKF)算法[10]、基于多模自适应估计的无迹卡尔曼滤波(UKF)算法[11]等，但这些算法都是利用仿真数据来实现的。目前，美国戈达德太空飞行中心(GSFC)利用SEXTANT设备的观测数据开展轨道确定试验，但尚未检索到其具体定轨算法。XPNAV-1卫星在轨运行4年来，获取了大量观测数据，按计划完成了X射线探测器性能测试、典型目标脉冲星观测及脉冲星导航系统体制验证等试验任务[12-14]。利用XPNAV-1卫星对蟹状星云(Crab)脉冲星的完整观测数据，采用基于轨道法平面几何约束方法，验证了该卫星轨道改进的有效性和脉冲星导航系统体制[15]，但利用单颗脉冲星的观测数据长时间定轨过程存在发散问题。本文针对XPNAV-1卫星拓展任务及后续发展需求，建立其轨道力学模型和观测方程，并论述扩展卡尔曼滤波(EKF)算法和分段式定常系统(PWCS)的可观测性分析方法。

1 XPNAV-1卫星任务

XPNAV-1卫星是一颗整星质量为243 kg的小卫星。其科学试验目标是：1)在空间环境下实测验证X射线探测器性能，研究宇宙背景噪声对探测器作用机理；2)探测Crab脉冲星辐射的X射线光子，提取脉冲轮廓曲线；3)长时间累积探测脉冲星，积累观测数据，探索验证脉冲星导航体制。

XPNAV-1卫星采用降交点地方时为6:00时的太阳同步晨昏轨道，其轨道半长轴为6 878.137 km，轨道倾角为97.4°，偏心率为0.010 3。该卫星的有效载荷为两种不同类型的X射线探测器。其中一种是掠入射式聚焦型探测器，其有效面积为30 cm2，探测能谱范围为0.5～10 keV，时间分辨率为1.5 μs，视场为15′；另一种是微通道板准直型探测器，其有效探测面积为1 200 cm2，探测能谱范围为1～10 keV，时间分辨率为100 ns，视场为2°。针对核心试验任务需求，优选了8颗目标脉冲星进行重点观测研究，其脉冲星编号、J2000.0国际天球参考系下的赤经和赤纬、脉冲周期及光子流量如表1所示。表1中，序号1～4为独立旋转供能脉冲星(IRP)，其中第1颗脉冲星B0531+21为Crab脉冲星，序号5～8为X射线双星(XB)。

XPNAV-1卫星在轨运行近4年来，对8颗目标脉冲星、超新星遗迹以及银河系X射线等天体进行了观测研究。由表1可见，XB的光子流量比IRP高2～3个数量级，观测XB的主要目的是进行探测器性能测试。考虑到XB系统本身轨道动力学的复杂性，本文主要利用上述4颗IRP开展自主定轨算法研究。

表1 8个候选观测对象参数列表

2 XPNAV-1卫星自主定轨数学模型

XPNAV-1卫星轨道为太阳同步轨道，在轨道可见性分析过程中要求脉冲星方向与太阳夹角大于45°，与月球夹角大于5°，并且要考虑地球遮挡、规避南大西洋异常区和高纬度极区辐射等约束条件。针对该卫星轨道特点，建立轨道力学模型和自主定轨观测方程。

2.1 轨道力学模型

XPNAV-1卫星轨道力学方程可表示为:

(1)

式中:X(t)为状态变量；f[X(t),t]为轨道动力学方程；w(t)为系统噪声项。

在该卫星自主定轨研究中，采用只考虑地球中心引力加速度和J2项非球形摄动加速度的卫星轨道力学模型，代入式(1)得到卫星轨道力学方程:

(2)

式中:r=[rx,ry,rz]T为卫星位置矢量；v=[vx,vy,vz]T为卫星速度矢量；X=[rx,ry,rz,vx,vy,vz]T为卫星状态矢量；w=[wx,wy,wz,wvx,wvy,wvz]T为系统噪声；μ为地球引力常数；J2为带谐项系数；Re为地球半径;r为卫星到地心的距离。

2.2 自主定轨观测方程

在太阳系质心(solar system barycenter,SSB)坐标系下，根据脉冲星导航的基本原理，可以得到脉冲星导航的基本测量方程[16]：

(3)

式中：tSSB和tsat为脉冲星同一脉冲到达SSB和卫星的时间；n为观测脉冲星方向向量；RSC为卫星相对于SSB的位置矢量；c为光速；δt为星载原子时钟偏差。

相较于脉冲星测距精度和观测数目等因素，卫星钟差对定轨精度的影响较小。不失一般性，针对本文开展的自主定轨研究，卫星钟差忽略不计，并考虑地球、卫星和SSB三者位置关系，将式(3)改写如下：

c·(tSSB-tsat)-n·RE=n·R

式中：RE为地球相对于SSB的位置矢量；R为卫星相对于地球质心的位置矢量。

在tk时刻卫星观测1颗脉冲星时，测量方程可以表示为：

Zk=HkXk+Vk

(4)

在X射线脉冲星定轨中决定定轨精度的主要因素是脉冲到达时间(TOA)的精度，TOA观测量的测量噪声方程可以估算为[17]：

(5)

式中：FB为X射线背景辐射流强；Fx为X射线辐射流强；tobs为信号观测时间；P为脉冲周期；W为脉冲宽度；pf为脉冲调制度；Ae为探测器有效面积。

3 EKF算法及PWCS分析方法

标准卡尔曼滤波主要针对随机离散线性系统模型，针对轨道动力学非线性的特性，采用EKF来解决随机连续非线性系统问题。并在对系统进行滤波处理之前，用PWCS方法定性地分析系统状态的可观测性，这是反映滤波收敛精度和速度的重要指标。

3.1 EKF导航滤波算法

EKF的基本思想是将导航系统的状态方程进行泰勒展开，忽视二阶及二阶以上的高阶项，将展开式的一阶线性项作为原状态方程的近似表达式，再采用线性问题的标准卡尔曼滤波方法对近似的状态量进行估计，从而得到最终状态的次优近似值。

对式(2)和式(4)相应简化后可以得到如下系统模型：

(6)

式中：t为时间参数；w(t)和v(t)分别为系统噪声和观测噪声，且彼此不相关。

式中：δX(t)为两者的状态估计误差。

进而得到误差状态的动力学方程为：

式中：F(t)为雅可比矩阵。进一步将状态方程进行离散化处理，得到：

δXk=Φk,k-1δXk-1+Wk-1

(7)

式中：Φk,k-1=I+F(tk-1)·T，为k-1时刻到k时刻的状态转移阵；T为滤波周期。

由式(4)和式(7)可以建立系统的离散型线性干扰方程：

扩展卡尔曼滤波流程如下：

(1)滤波初始化

根据估计状态信息，给出状态变量及其相应误差方差阵初始值：

(2)状态更新

(3)量测更新

Pk=(I-KkHk)Pk/k-1(I-KkHk)T+

3.2 PWCS可观测性分析方法

不论采用卡尔曼滤波或其他滤波方式，系统状态变量的可观测性是滤波收敛的前提条件。针对X射线脉冲星自主定轨系统是典型的非线性时变系统，提出基于分段式线性定常系统(PWCS)的可观测度分析方法。在足够小的时间区间内，可以忽略线性时变系统的系数矩阵变化，在该时间区间内就可以把时变系统近似为一个分段线性定常系统进行处理，这往往可以大大简化分析过程。

离散型PWCS系统的一般形式为：

式中：j=1,2,…,q，为时变系统分段间隔序号。于是，系统总的可观测性矩阵(TOM)可以表示为：

式中：

根据PWCS分析定理，如果系统满足

null(Qj)⊂null(Φj), 1≤j≤q

则有

因此，可以用QS来代替QT分析系统的可观测性，使问题得到简化。又因系统状态呈递推关系，每个时间段的系统状态初值就是前一时间段的系统状态终值，则可观测性矩阵仍具有累积继承的性质。这样可以用某时段j的可观测性矩阵Qj进一步代替QS，使系统的可观测性分析进一步简化。

系统的可观测矩阵Qj可以表示为：

(8)

式中：

(9)

S中每个元素的量级为μ/r3，对于轨道高度为500 km的卫星，其数量级为10-6，近似为零，将式(9)其代入式(8)，并进行初等行变换，得到：

通过提取可观测矩阵，进一步进行可观测度的计算，用可观测矩阵的条件数来对可观测性进行量化，对Qj进行奇异值分解：

Qj=UΛVT

定义可观测矩阵的条件数如下：

式中：σmax和σmin分别为M的最大奇异值和最小奇异值。矩阵的条件数η越大，说明矩阵越接近奇异。

系统可观测性矩阵的条件数能够反映系统状态的可观测程度，即系统可观测性矩阵的条件数越大，系统的可观测度越差。如果系统可观测性矩阵的条件数是无穷大，则系统不可观测；如果系统可观测性矩阵的条件数等于1，系统的可观测性最好。

4 数值分析试验

4.1 XPNAV-1卫星的观测数据分析

目前，XPNAV-1卫星仍正常在轨运行，对Crab脉冲星和3颗低流量IRP进行了详细观测，积累了大量观测数据。从实测数据分析来看，掠入射式聚焦型探测器观测的数据有较高的信噪比，因此在本文定轨算法研究中，主要利用该探测器采集的数据进行数值分析试验。

自2016年11月18日以来，XPNAV-1卫星一直对Crab脉冲星进行重点观测，除每年4月30日至8月1日期间Crab方向和太阳夹角小于45°，以及考虑与月球夹角小于5°不可观测外，其余每天进行观测，每个观测弧段持续观测时间为30～50 min。目前，通过累积观测获取了比较完整的Crab脉冲星观测数据，并提取得到精化的脉冲轮廓，通过计时拟合建立了精度为55.14 μs的计时模型。

此外，该卫星在2017年2月21日至2017年8月31日，分别对3颗低流量IRP进行了长期观测，每个轨道周期持续开展5～15 min低流量脉冲星探测，累积探测数据。经统计分析得到了3颗低流量IRP的流量统计以及能谱分析结果。但由于探测器有效面积较小，脉冲星光子流量较低，目前未能提取到低流量脉冲星完整的脉冲轮廓。

XPNAV-1卫星对Crab脉冲星的进行了完整的观测。其中，时间跨度从2016年11月18日至2017年3月26日的观测效果较好，共179段观测数据，每个观测弧段平均观测时间约45 min。对该段实测数据做如下处理：1)在0.5～9.0 keV的能带内选择光子；2)筛除有效观测时间短于40 min的观测数据；3)筛除通量小于14个光子/s的观测数据。最后从中选取时间跨度为30 d的数据用作数值试验分析，共10组观测值，其TOA误差及相应的测距误差如表2所示，实测数据的起止时间为2017年1月5日至2017年2月5日。

但是由于实测数据的信号积分时间较长，观测量不充足以及测量序列不规则,所以仅仅只利用以上Crab的实测数据会导致滤波结果的发散。因此,基于Crab脉冲星实测数据的测距精度来模拟生成额外的观测值，对30 d内的实测数据进行补充，对观测周期进行压缩，每500 s生成一个观测值。

考虑到单颗脉冲星定轨系统的可观测性差，再结合B0540-69，B1509-58，J1846-0258这3颗低流量IRP，利用1颗、2颗、3颗、4颗脉冲星的观测数据组成数据组，分别进行轨道确定，分析观测脉冲星数目对定轨精度的影响。利用Crab脉冲星实测数据，由Crab脉冲星观测TOA与Crab脉冲星模型预报TOA的差，可以得到每个观测时间段的TOA误差(见表2)，计算得到TOA均方根误差为99.05 μs，对应29.70 km的测距均方根误差。由于3颗低流脉冲星光子流量较低，探测器有效面积较小，未能得到其完整的脉冲轮廓。由Crab脉冲星的完整观测数据反算出背景噪声，再根据3颗低流量脉冲星的实测光子流量统计以及探测器实际有效面积，利用式(5)的估计方程推算得到其测距精度，并由测距精度模拟生成观测数据。表3给出了4颗脉冲星的赤经、赤纬和测距精度。

表2 基于Crab脉冲星观测结果的TOA误差及测距误差

表3 4颗IRP的测距精度估计

4.2 颗脉冲星的可见性分析

在XPNAV-1卫星轨道覆盖性分析过程中，影响其对目标脉冲星观测的约束条件是：1)地球处于观测目标与卫星之间，遮挡观测信号；2)观测目标与太阳之间的夹角小于45°，或者观测目标与月球之间的夹角小于5°，高能粒子对探测器造成影响。XPNAV-1卫星的轨道周期约为94 min，图1和图2给出了4颗脉冲星在一个轨道周期内的可见性状况。

图1 4颗脉冲星的可见性分析Fig.1 Visibility analysis of four pulsars

图2 一个轨道周期内可见星个数Fig.2 Number of visible pulsars in an orbit period

由图2可以看出，一个轨道周期内的任意时刻均最少有2颗脉冲星可被观测到，脉冲星可见的时间百分比为100%，依靠4颗目标脉冲星的观测信息，可以实现全时段定轨，从每个观测弧段的可见星中选取单星进行导航，可以保证观测信息的完备性。

4.3 滤波系统状态可观测性分析

系统可观测度的数值试验结果如图3所示，纵坐标表示可观测矩阵条件数，且取底数为10的对数。

表4给出了观测不同脉冲星数目时，系统可观测矩阵的条件数。当全程只观测1颗或2颗脉冲星时，可观测矩阵Qj条件数的数量级较大，矩阵接近奇异，定轨系统可观测性差，无法实现滤波过程收敛。当观测3颗脉冲星时，可观测矩阵条件数明显变小，系统可观测度提高。当观测4颗脉冲星时，可观测矩阵Qj的条件数均值约为57，系统可观测度较好，明显优于前3组结果。可以总结出当观测脉冲星数目增加时，可观测矩阵条件数变小，系统可观测度提高，数值试验结果与上节理论分析结果一致。

图3 基于4颗脉冲星的可观测矩阵条件数Fig.3 Conditional number of observable matrix based on observing four pulsars

表4 脉冲星数量对系统可观测度影响

4.4 自主定轨结果分析

利用第3组脉冲星数据组进行定轨计算，得到轨道误差和速度误差序列分别如图4和图5所示，进一步统计分析得到位置误差均方根为56.65 km，速度误差均方根为0.055 1 km/s。利用第4组脉冲星数据进行定轨计算的轨道误差和速度误差序列分别如图6和图7所示，进一步统计分析得到位置误差均方根为40.24 km，速度误差均方根为0.043 3 km/s。

图4 基于观测3颗脉冲星的轨道误差Fig.4 Position error based on observing three pulsars

图5 基于观测3颗脉冲星的速度误差Fig.5 Speed error based on observing three pulsars

图6 基于观测4颗脉冲星的轨道误差Fig.6 Position error based on observing four pulsars

图7 基于观测4颗脉冲星的速度误差Fig.7 Speed error based on observing four pulsars

由图4～图7可见，当观测3颗或4颗脉冲星时，滤波过程收敛，并且随着观测脉冲星数目的增加，轨道误差减小，定轨精度提高。

5 结束语

针对XPNAV-1卫星拓展试验任务及脉冲星导航后续发展需求，利用多颗脉冲星实测数据，对EKF自主定轨算法进行数值试验验证。研究结果表明：当仅观测1颗脉冲星时，系统的可观测矩阵条件数极大，滤波结果很快发散；当观测2颗脉冲星时，系统可观测度提高，但条件数依然较大，滤波结果经过长时间依旧发散；当观测3颗脉冲星时，系统可观测度明显提高，脉冲星对轨道完全覆盖，滤波过程收敛，轨道确定精度为56.65 km；当观测4颗脉冲星时，系统可观测度进一步提高，脉冲星对轨道覆盖性能好，滤波过程收敛，轨道确定精度为40.24 km；在所选脉冲星满足该卫星整个轨道周期覆盖的条件下，每个弧段仅需要观测1颗脉冲星就可以进行定轨计算，其定轨精度取决于所选脉冲星的空间分布和数量。本文提出的EKF自主定轨算法过程收敛，试验结果验证了算法的合理性和有效性。EKF算法适用于线性高斯白噪声的系统状态估计，而实际的航天器轨道确定数学模型往往是非线性或非高斯的，存在线性化误差、建模误差和有色噪声等问题。下一步将利用脉冲星实测数据，开展基于UKF、H∞滤波和粒子滤波等算法研究，对比分析各种算法性能，提出最优的定轨算法解决方案，为脉冲星导航算法研究和后续空间飞行试验提供技术储备。