利用GECAM卫星Crab脉冲星观测数据的定轨分析

韩大炜，郑世界，庹攸隶，葛明玉，宋黎明，李新乔，文向阳，熊少林

中国科学院 高能物理研究所 粒子天体物理重点实验室，北京 100049

脉冲星被称为天然的GPS卫星，其脉冲信号的长期时间稳定度很高［1］。空间飞行器利用其脉冲信息可实现自主导航，从而为太阳系内距离地球较远的行星乃至星际旅行的航天器进行导航定轨、守时等提供帮助。自1967年人类发现第一颗脉冲星［2］、1974年Downs提出射电脉冲星导航概念［3］、1981年Chester和Butman提出X射线脉冲星星际导航［4］以来，X射线脉冲星导航技术一直被世界各国争相研究。

国际上脉冲星导航的研究经历了从理论提出到算法研究、空间试验等阶段。目前开展在轨试验仍是脉冲星导航研究的主要手段。1999年美国在ARGOS卫星上开展了非常规恒星定位试验（Unconventional Stellar Aspect，USA）进行脉冲星导航研究［5-7］。2017年NASA实施了国际空间站上的X射线计时导航技术探索（The Sta⁃tion Explorer for X-ray Timing and Navigation Technology，SEXTANT）项目，使用的载荷称为中子星内部组成探测器（Neutron Star Interior Composition Explorer，NICER），从而开展了脉冲星自主导航试验；NICER由56个相同的X射线望远镜组成，总有效面积约为1800 cm2，每个望远镜都由一个集中的X射线光学元件和一个像素硅漂移探测器（SDD）组成；SEXTANT的试验结果显示通过观测4颗毫秒脉冲星获得的定位精度可达5 km内（1σ，其中σ为标准差）［8］。

21世纪以来中国的X射线脉冲星导航研究蓬勃发展，涉及时空坐标体系、导航脉冲星源［9］，导航理论、算法与误差分析及空间X射线探测器研制等各个方面，研究者们利用空间高能探测设备积极开展导航理论研究、脉冲星导航体制探索［10］及空间在轨验证。迄今为止开展的空间试验包括天宫二号空间实验室上的“天极望远镜”——伽马射线暴偏振探测器（Gamma-Ray Bursts Polarimeter，POLAR），主要科学目标是用于伽马暴的偏振测量，探测能段为15～500 keV。其具有很大的视场（超过2π立体角）且在该能段具有较大的有效面积（约200 cm2）、较高的时间分辨率（80 ns），因此也能探测到Crab脉冲星的X射线辐射信号，从而成功进行了脉冲星导航定轨试验，定轨精度为10 km量级［11］。

脉冲星导航试验卫星01星（X-Ray Pulsar Navigation Test Satellite No.1，XPNAV-1）使用Wolter-I型望远镜［12-13］和SDD探测器观测Crab给出了初步的导航精度，利用25次观测的数据计算出脉冲周期的均方根误差（RMS）为0.00264，对应的脉冲到达时间（Time of Arrival，TOA）误 差 为89.2 μs，距 离 误 差 为26.8 km［14］。

中国首颗X射线天文卫星——“慧眼”卫星（Insight， Hard X-ray Modulation Telescope）使用3台望远镜、通过对蟹状脉冲星的观测成功确定了轨道。试验分别使用“慧眼”卫星3台望远镜（高能X射线望远镜（High Energy X-Ray Tele⁃scope，HE）、中能X射线望远镜（Medium Energy X-Ray Telescope，ME）和低能X射线望远镜（Low Energy X-Ray Telescope，LE））各5天的观测数据，均成功确定了轨道参数，综合所有数据可知定位精度在10 km（3σ）［15］。

本文通过分析GECAM卫星对Crab脉冲星的观测数据说明基于轨道动力学的脉冲轮廓波形典型显著性分析方法进行脉冲星定轨的原理和步骤，并给出计算结果和与其他卫星的对比。

1 GECAM卫星

怀柔一号引力波暴高能电磁对应体全天监测器（Gravitational Wave High-Energy Elec⁃tromagnetic Counterpart All-Sky Monitor，GECAM）卫星是中国科学院“空间科学”（二期）战略性先导科技专项的天文卫星，其科学目标是对引力波伽马暴、快速射电暴高能辐射、特殊伽马暴和磁星爆发等全天随机爆发的高能暂现源的监测研究。GECAM卫星采用600 km高度低轨和29°低倾角的天文观测轨道［16］，包含2颗相同的运行于同一轨道面的微小卫星，分别为GECAM-A和GECAM-B。伽马射线探测器（Gamma-Ray Detector，GRD）是GECAM的主要探测器载荷，其主要任务和功能是实现对空间伽马射线暴发源的宽能段能谱和时变观测，并通过多个探头实现对全天的监测及对伽马暴的定位。

根据GECAM卫星的构型，GRD（图1中的圆形探测器）近于均匀地分布在椭球形的穹顶结构上，每颗卫星共有25个GRD模块，均安装在卫星的穹顶舱，每个GRD模块均能探测从前部入射的伽马射线光子，视场约为2π。25个GRD分别指向不同方向，在立体角上基本均匀分布［17］。平均情况下，约12个探头可观测到Crab脉冲星。某个方向上所有探头理论上均可接收X射线光子，有效面积合计约250 cm2［18］。GRD的主要技术指标如表1［17］所示。

图1 GECAM有效载荷的构型设计Fig.1 Configuration design of GECAM payload

表1 GECAM-GRD的主要技术指标［17］Table 1 Main technical parameters of GECAM-GRD［17］

总之GECAM卫星的探测能区宽、时间分辨率和时间精度都很高［19］，观测天区覆盖广。因此GECAM能长时间接收到Crab脉冲星发射出的X射线光子，可用于开展脉冲星导航试验。

2 SEPO算法

基于轨道动力学的脉冲轮廓波形典型显著性分析方法（Significance Enhancement of Pulseprofile with Orbit-dynamics，SEPO）是结合脉冲轮廓和轨道动力学显著性分析的方法。经典导航领域常用的卡尔曼滤波方法已在中国广泛应用于X射线脉冲星的仿真分析和在轨导航试验中［20-21］。笔者采用与之不同的SEPO算法，此方法已在天极望远镜（POLAR）和“慧眼”卫星（Insight-HXMT）中进行了试验［11，15］。

SEPO算法的基本原理是轨道动力学模型能在短时间内得到高精度的轨道预报，但具有长期性漂移，特别是低轨卫星由于大气模型的变化导致轨道的长期预报偏差很大。另外当初始轨道存在偏差时预报轨道与真实轨道自然会出现偏离，且该偏离随着时间的延长越来越大。根据真实轨道能计算出脉冲轮廓，而轨道偏离会导致脉冲轮廓发生变形，导致轮廓信号的显著性降低。与“慧眼”卫星类似，GECAM卫星工作在约600 km高度的低地球圆轨道，其轨道预报的变化效应和对轮廓的影响比较类似，同样可采用SEPO算法进行定轨。

SEPO算法的整个计算过程包括利用轨道动力学外推、计算脉冲波形的显著性和基于网格搜索的轨道拟合3大步骤，详见文献［11，15］。

SEPO算法的优势是仅使用一颗脉冲星即可实现卫星定轨。目前的试验结果表明这种方法更适用于轨道根数变化剧烈且不存在完全“平行”的轨道，这样这些变化可很好地显示在脉冲轮廓中。另外由于脉冲星在不同能段的轮廓形状会有所差异，所以不同类型的空间探测器由于工作原理和能段不同无法套用“脉冲标准轮廓”，SEPO算法使用自身探测的轮廓突破了该限制。

相较于最初在POLAR试验中的应用，使用的SEPO算法做了进一步的完善，体现在利用轨道动力学加上偏差计算轨道时由POLAR时的向后外推更改为往前推算。因此计算结果可得到“现在”的轨道根数，而非POLAR计算时得到的“过去”的某个起始点的轨道根数，因而本文的计算具有了更大的工程应用价值。

3 数据筛选和预处理

通常情况下使用X射线脉冲星的观测数据进行导航定轨的处理步骤分为4步：事例挑选、太阳系质心修正、自转频率搜索和脉冲轮廓叠加、脉冲TOA的计算。

对GECAM卫星的观测数据处理利用了前3个步骤，使用SEPO算法并未用到TOA的计算。具体步骤如下。

3.1 事例挑选

1） 剔除地球遮挡数据

根据Jodrell Bank射电天文台提供的Crab坐标值可知其在J2000坐标系下为赤经RA=05d34m31.972s、赤纬DEC=22h00m52.07s，可将该RA、DEC转换到直角坐标系（xcrab，ycrab，zcrab）中。卫星的姿态数据中提供了卫星在J2000的直角坐标系中每秒的位置（x，y，z），这样可先计算卫星的位置向量与Crab方向向量之间的夹角，从而可计算GECAM卫星与Crab的连线与地心之间的夹角。若Crab完全被地球遮挡（夹角小于67°即被地球遮挡），则剔除该部分的光子事例。

2） 探测器筛选

根据卫星的姿态及各GRD探头的位置可计算Crab脉冲信号分别入射到25个GRD探测器的夹角。在卫星的姿态文件中给出了逐秒的卫星姿态四元数。根据姿态四元数及卫星载荷信息中各个GRD探头与卫星主轴的位置关系可逐秒计算得到Crab脉冲星的位置矢量与各个GRD探头的夹角。对于各个GRD探测器探测到的光子事例，选择Crab与各探头夹角小于70°时间段内的光子。

3） 能段筛选

选择GECAM的能道不大于100的高增益事例，对应的是约10～100 keV的低能段事例。这样能使来自于非目标源Crab的光子事例占比较小，增加Crab脉冲信号的显著性。

选取GECAM-B星的数据，时间段为2021年7月23日—8月31日，总光子数为20.98亿个。对所有光子进行折叠并扣除本底后统计脉冲轮廓的光子数，脉冲轮廓的光子数约为800万个，约占总光子数的0.36%。

3.2 太阳系质心修正

太阳系质心坐标系是一个较好的惯性参考系，为得到准确的脉冲星脉冲到达时间，需将探测到的脉冲光子到达时间转换到太阳系质心坐标系中。这主要包括真空延迟修正、相对论修正及不同时间系统之间的转换等。使用太阳系星历DE405进行太阳系质心修正。

3.3 脉冲轮廓折叠

由于直接利用了Fermi卫星的Crab脉冲星星历，因此未作自转频率搜索。根据星历计算出所有观测的X射线脉冲光子的相位，统计所有光子的相位信息得到Crab脉冲星的X射线脉冲轮廓。

4 计算结果与分析

利用SEPO算法计算时首先给出估计的初始轨道参数和参数空间，利用轨道动力学模型进行30天的轨道预报，计算得到基于该预报轨道的观测脉冲轮廓；然后通过对此观测脉冲轮廓的显著性分析搜索最优轨道参数。具体流程如图2所示。

图2 脉冲星定轨算法的处理流程Fig.2 Processing flow of pulsar orbit determination

进行轨道外推时使用的积分器为Runge-Kutta-Fehlberg7（8），积分步长为60 s，在此情况下可满足轨道外推的需要，考虑的主要动力学模型参数如表2所示。

表2 轨道动力学模型中的参数与设置Table 2 Parameters and settings in orbital dynamic models

根据轨道动力学模型通过起始点向前推算了30天轨道位置与GECAM卫星真实测量的位置存在偏差，计算结果如图3所示。可见偏差的均值为−5.10、−3.17、−0.77 km，1σ标准差分别为7.56、9.28、5.22 km。在计算的30天内偏差总的趋势是逐渐放大，呈现比较复杂的变化趋势，平均在km级。计算Crab脉冲星的轮廓分为500个相位，每个相位的周期约为66 μs，折算到距离为20 km。可见轨道计算的偏差对脉冲星轮廓的影响在一个脉冲轮廓相位的范围内，未影响计算结果。

图3 GECAM轨道起点向前外推30天的位置与真实测量值的偏差Fig.3 Deviations between positions extrapolated for⁃ward for 30 days from orbit starting point of GECAM and true measured values

GECAM卫星的轨道状态用J2000坐标系下的6个轨道根数表示，分别是半长轴（a）、偏心率（e）、轨道倾角（i）、升交点赤经（ω）、近地点幅角（Ω）和平近地点角（w）。

首先选取GECAM卫星在观测结束时刻的位置作为初始值，兼顾轨道参数覆盖范围和计算效率，选用轨道根数初始值的偏移量如表3所示。

表3 轨道参数初始值的偏移量Table 3 Bias of initial value of orbital elements

在此参数空间内的每个格点处分别计算得观测脉冲轮廓。图4为观测Crab的原始轮廓和半长轴分别改变−520 m和+480 m的模拟轮廓，可见其发生了畸变。

图4 GECAM观测Crab脉冲星轮廓和半长轴变化的仿真轮廓Fig.4 Profile of Crab pulsar observed with GECAM and profile simulated with semimajor axis

直接对观测脉冲轮廓进行处理。定义脉冲轮廓的显著性χ2为

式中：p（i）为经过归一化的脉冲轮廓；pˉ为p（i）的平均值；k为脉冲轮廓的相位数目。

计算每个轨道参数格点对应脉冲轮廓的显著性，结果如图5所示，图中计算结果为不同参数的带误差卡方计算结果，拟合结果为最优值拟合结果（高斯拟合），0值点为实际轨道根数值。可见脉冲轮廓的显著性χ2随轨道参数的偏离有显著变化，在真实值（0值）附近χ2具有极大值。通过高斯函数曲线拟合得到的最大值可分别给出6个轨道参数的最优解。以半长轴为例，当起始时刻的半长轴偏离±50 m时计算得到的脉冲轮廓会发生畸变，其显著性χ2也从3.10×104下降到2.85×104（偏离±50 m）。

图5 χ2 随轨道参数偏移量的变化结果Fig.5 χ2 variation results with orbit elements bias

采用Sheikh等计算脉冲TOA误差的方法［6，22］计算轨道参数误差，每个参数的误差σ为

式中：W为参数最优值的高速拟合过程中的半高全宽（FWHM）；S/N为信噪比。采用相同的计算方法［6］：

式中：S/Nfiltered为加了滤波限幅的信噪比。为限制S/N可能出现极大的异常值，采用限幅公式对式（3）进行限制，最大为1000。从而得到每个轨道根数的最优值和3σ误差（99.7%置信度），如表4［11，15］所示。

表4 轨道参数最优值和误差［11，15］Table 4 Best estimated values of elements and errors［11，15］

通过比较表4［11，15］可见3次试验中半长轴的误差比较接近；分析原因是半长轴的轻微改变就会导致卫星运动的轨道周期发生变化，使卫星位置的变化随时间不断累积，从而使脉冲轮廓发生较大的改变，因此SEPO方法对半长轴有较好的约束。其余5个轨道参数的误差中“慧眼”卫星是最优的，GECAM卫星与POLAR望远镜相比偏心率、倾角和升交点赤经的误差较小，而近地点幅角和平近地点角的误差较大。

5 总结与展望

利用“怀柔一号”GECAM卫星对Crab脉冲星的在轨观测数据进行了脉冲星导航定轨的原理验证。通过建立基于轨道动力学模型和轮廓分析的导航算法初步实现了利用脉冲星探测进行空间飞行器的定轨。针对GECAM-B星40天的在轨观测数据得到了初始历元时刻6个轨道参数的最优值与误差，其最优值均与真实值在3倍误差范围内吻合。3σ的定轨误差与“慧眼”卫星（Insight-HXMT）脉冲星导航的试验结果［15］相比除半长轴和倾角相近外，其他4个轨道根数的误差是“慧眼”卫星的2～4倍。证明了GECAM卫星同样可用于确定轨道根数，但由于本底误差远高于“慧眼”卫星，故导航误差也大于“慧眼”卫星的10 km。

另外值得注意的是由于GECAM作为具有极宽视场的全天监测的天文卫星，尽管选取了70°的视场本底仍然远高于“慧眼”卫星，因此选取40天的数据后轨道根数变化显著性的拟合才有比较明显的概率分布，且最优值偏离在3倍标准差之内。可见GECAM用于定轨所需的时间远超过“慧眼”卫星的5天，而与POLAR试验中的31天接近。这个结果也说明本底会对计算造成较大的影响。

需要说明的是，截至目前的3次试验均未考虑航天器上的时钟时间偏移及脉冲星内在的红噪声和跃变属性等的影响。

从望远镜本身的研制而言，GECAM卫星的特点在于载荷对卫星平台的资源需求较低。相较于“慧眼”卫星观测脉冲星时指向特定天区的方式，GECAM采用大视场的设计，指向精度和姿轨控制的要求低。定轨数据处理过程中通过卫星姿态和观测角度的方法排除一部分非Crab的本底事例。另外单颗GECAM卫星的质量只有180 kg，单星发射包络尺寸不大于580 mm（X）×1050 mm（Y）×1364 mm（Z）［16］，与“慧眼”卫星的约2500 kg相比GECAM卫星的质量大幅度减少。

而安装在中国天宫二号上的POLAR载荷尽管其本身质量较小，包含电控箱只有30 kg，但依赖于空间实验室提供能源、姿轨控、测控等支持。GECAM卫星实现了载荷卫星一体化研制。

总之，本文在“慧眼”等传统的天文卫星和POLAR等空间实验室的研究基础上，采用卫星小型化和载荷轻型化之后的GECAM卫星成功实现了定轨，为中国未来脉冲星定轨和导航空间试验的开展提供了新思路。