“天枢Ⅱ号”X射线脉冲星导航动态模拟系统及实验验证

2023-03-12盛立志郑伟苏桐张大鹏王奕迪杨向辉徐能李治泽

航空学报 2023年3期

盛立志，郑伟，苏桐，张大鹏，王奕迪，杨向辉，徐能，李治泽

1．中国科学院西安光学精密机械研究所，西安 710119 2．国防科技大学空天科学学院，长沙 410073

X射线脉冲星导航（X-ray Pulsar Navigation， XPNAV）是以高稳定的脉冲星辐射信号的时空信息为参考，从而实现航天器的位置、速度信息确定的新型导航方式，导航系统运行过程中不需要地面观测站的校准，具有完全自主导航的能力［1-6］。目前，国内外已围绕X射线脉冲星导航XPNAV进行了大量的理论及初步实验研究。美国国家航空航天局（NASA）于2017年发射了中子星内部结构探测器（Neutron-star Interior Com⁃position Explorer， NICER），NICER搭载了56个口径为10 cm的X射线聚焦镜，开展了X射线脉冲星计时研究及脉冲星导航验证研究，并通过对毫秒脉冲星的观测得出了导航精度可达10 km的结论［7］。中国于2016年11月发射了XPNAV-1，这是中国首颗脉冲星导航试验卫星。XPNAV-1中搭载了准直型和聚焦型两种X射线探测器，其中，聚焦型探测器有效探测能段为0.5～10 keV，有效面积为30 cm2［8］。卫星对PSR B0531+21脉冲星进行了长时间的观测并获得了有效的数据，实现了对X射线脉冲星“看得见”的目标。此外，中国于2016年9月发射的天宫二号上的伽马暴偏振探测仪（POLAR）［9］和2017年6月发射的用于空间天文研究的硬X射线调制望远镜（Hard Xray Modulation Telescope， HXMT）［10］均对X射线脉冲星进行了观测，并开展了导航试验空间验证。虽然上述空间观测和验证已经获得不少成果，但XPNAV工程化应用仍有一定的距离，完善的脉冲星数据库的建立、脉冲星信号降噪及高精度脉冲到达时间TOA（Timing of Arrival）估计、轻型化高效率X射线探测系统研制等关键技术仍需要深入研究。

X射线无法穿透大气层，导航工程化之前必须要开展X射线脉冲星地面模拟实验研究，为导航相关理论提供实验验证平台。NASA的戈达德空间飞行中心GSFC开发了X射线导航实验平台（X-ray Navigation Laboratory Testbed， GXLT），X射线模拟源采用调制X射线源（Modulated X-Ray Source， MXS），探测器为硅漂移探测器SDD（Sicicon Drift Detector）。MXS主要是通过紫外LED调制光电阴极，产生的光电子经过电子倍增器件后在阳极高压作用下轰击阳极靶材，从而产生调制的X射线。GXLT可以进行3个层次的仿真：纯软件层面、不模拟光子到达过程的仿真，纯软件模拟光子到达过程的软件层面闭环模拟实验，以及全物理闭环仿真实验［11］。中国也开展了脉冲星导航地面模拟实验验证技术的相关研究，并取得了一定的进展。国防科技大学于2011年研制了国内首套X射线脉冲星导航半实物仿真系统，可以利用探测器实时探测X射线光子信号，将光子到达时间通过数据采集传入导航计算机并进行处理，实现TOA估计和导航滤波算法，从而实现X射线脉冲星导航的半实物仿真［12］。中国航天科技集团五院建立了脉冲星导航探测器地面标定装置，用于探测器脉冲轮廓还原能力和探测效率、能量分辨率等指标的标定［13］。西安电子科技大学基于可见光搭建了脉冲星地面验证系统，模拟源采用线性较好的蓝色发光二极管LED（Light Emitting Diode），探测器为光电倍增管PMT（Photomultiplier Tubes），通过减光片将光子流量控制在单光子探测范畴，实现脉冲星信号地面模拟实验研究［14-15］。中国科学院西安光学精密机械研究所对X射线模拟源及探测器等关键技术进行了大量研究，提出了一种用于产生脉冲X射线的栅控X射线管，可以对任意脉冲星的轮廓进行模拟，开展了基于X射线脉冲星导航的静态模拟实验［16-18］。目前文献报道的基于X射线波段的半实物实验系统开展的大多局限于静态实验，并没有包含航天器轨道运动特性。

本文面向脉冲星导航应用研究的需求，采用栅控X射线管的方法开展基于脉冲星辐射特性和航天器轨道模型的X射线脉冲星动态信号模拟，构建地面实验系统，对X射线脉冲星静态辐射特性及包含各种动态效应的动态脉冲信号进行高精度模拟和效果评估。

1 “天枢Ⅱ号”X射线脉冲星导航动态模拟系统

“天枢Ⅱ号”X射线脉冲星导航地面实验系统的组成如图1所示，第1个模块主要实现对X射线脉冲星及背景噪声的模拟。PC软件控制终端、动态信号发生器、栅控信号调制器组成了脉冲星模拟控制器，用该控制器驱动栅控X射线管发射所要模拟的X射线脉冲信号，模拟不同的工作状态，利用GPS（Global Position System）驯服的原子钟来给出精确时标；直流X射线管和二次靶转台构成了单能X射线模拟器，在上位机软件控制终端及驱动电源的作用下可以实现对背景X射线噪声的模拟。第2个模块为真空实验平台，主要由真空计、真空管道、分子泵、机械干泵等构成，模拟太空中的真空环境。X射线聚焦光学（Nested X⁃ray Focusing Optics， NXFO） SDD探测器、制冷及电压控制模块、SDD电子学等构成了聚焦型探测系统，用于接收光子信息；接收到的光子信息在信号处理计算机中进行处理，完成相关实验验证工作。系统实物如图2所示。

图1 X射线脉冲星导航地面实验系统组成Fig. 1 Composition of X-ray pulsar navigation ground experiment system

图2 “天枢Ⅱ号”X射线脉冲星导航动态模拟系统Fig. 2 “TIANSHUⅡ”XPNAV ground experiment system

X射线脉冲星模拟源由脉冲星模拟控制终端和脉冲X射线发射源组成，实现对X射线脉冲星的模拟。脉冲星模拟控制终端用于实验参数的配置，包含航天器及星体的轨道信息和脉冲星辐射特征等信息，并根据设定的条件实时计算脉冲轮廓数据，用于脉冲X射线发射源的控制，实现导航的静态及动态实验仿真。

脉冲X射线发射源的核心器件为栅控X射线管，由灯丝、阴极罩、栅极、聚焦极、阳极靶、玻璃外壳组成，如图3所示。当灯丝加上电流，阳极加上高压，灯丝发射的电子就会在高压电场下加速运动并轰击阳极靶材产生X射线。如果栅极加上调制信号，轰击阳极靶的电子数目就随调制信号的强度的变化而变化，发射X射线的强度也随之变化。聚焦极的主要作用是改变球管内部电场分布实现对从阴极发射的电子汇聚，缩小电子束斑尺寸，从而降低电子束的弥散时间，提高频率响应能力。

图3 栅控X射线管结构图Fig. 3 Structure diagram of grid controlled X-ray tube

在模拟实验中，脉冲星模拟控制终端根据设定的实验条件实时计算脉冲轮廓幅度数据，然后经过数字信号传输接口发送给数据转换电子学，数据转换电子学将轮廓幅度信号转换为栅极电压信号，在栅极电压的作用下栅控X射线管发射与模拟脉冲星轮廓一致的X射线脉冲信号。

2 X射线脉冲星导航动态模拟信号

2.1 X射线脉冲星动态信号产生原理

脉冲星的标准脉冲轮廓是将长时间观测的光子到达时间序列转换到太阳系质心坐标SSB（Solar System Barycenter）处的太阳系质心坐标时，然后进行周期折叠得出的累积轮廓。在常规的静态模拟实验中，脉冲星模拟控制器依照标准脉冲轮廓模型计算X射线脉冲星轮廓数据。而在动态实验中，由于航天器及星体的运动会产生各种动态效应，航天器记录的光子到达时间序列周期折叠所得到的脉冲轮廓不再是标准脉冲轮廓，动态信号发生器依据航天器的运动状态和标准轮廓模型计算出动态脉冲轮廓数据，流程如图4所示。

图4 动态脉冲轮廓数据计算过程Fig. 4 Calculation process of dynamic pulse profile data

由于脉冲星模拟控制器为数字系统，模拟脉冲的相位并非连续的，因此只需根据计算周期转换等间隔的有限时间点。将航天器固有时τ按照等时间间隔Δτ0排列为τ0，τ1，…，τk，与之对应的坐标时为t0，t1，…，tk，有

式中：k∆τ=τk−τ0；k为固有时τ的等分数；rEAk为地心到航天器的距离矢量；rE为SSB到地心的距离矢量；VE为地球相对于SSB的速度矢量；VEAk为航天器相对地心的速度矢量。

根据式（1）就可将等间隔的固有时间序列转换为坐标系时间序列，坐标系时间序列并非等间隔。由于脉冲轮廓没有具体函数表达式，通常是等时间间隔存储有限的数据点。固有时间序列τ0，τ1，…，τk是等间隔排列的，而其对应的SSB坐标时时间序列tSSB0，tSSB1，…，tSSBk并非等间隔的。在地面模拟系统的脉冲信号发生器中，信号采样的频率是固定值，因而是等间隔的，因此需要通过插值的方法从等间隔的点轮廓数据点中恢复出任意时刻的轮廓数据。

航天器所在位置的TCB （Barycentric Coor⁃dinate Time）时间与SSB处TCB时间关系为

式中：tSSB为SSB处的TCB时；t为航天器处TCB时；n为太阳系到脉冲星视线方向的单位向量；r为SSB指向航天器的距离矢量；D0为太阳质心到脉冲星的距离；b为SSB到太阳质心的距离矢量；μs=GMs为太阳引力常数，G为引力常数，Ms为太阳质量。同样，对于选定的航天器轨道，r可以确定，其他量都为常量。对于航天器所在位置TCB时间tk，就能转化成SSB处TCB时间tSSBk。将式（2）表达为tSSBk=g(tk)，式（1）表达为tk=h(τk)，则

航天器接收的脉冲相位与所选坐标系无关，该相位值对应着SSB处tSSBk时刻的相位值，即

式中：ØSC(τk)为航天器处接收的脉冲相位；ØSSB(tSSBk)为SSB处的脉冲相位。

根据SSB处相位预测模型有

式中：f为脉冲星频率̈分别为其一阶导数和二阶导数，本文模拟的Crab脉冲星和PSR B1937+21周期稳定性较好，导数项可以略去，脉冲周期为P，有

2.2 生成脉冲星导航动态模拟轮廓数据

等间隔固有时间序列τ0，τ1，…，τk对应的SSB时间序列tSSB0，tSSB1，…，tSSBk并不是等间隔的，因此ØSC(τ1)，ØSC(τ2)，…，ØSC(τk)也不是等间隔的，不与标准轮廓的已知相位点重合。为得到这些相位点对应的脉冲强度，需要对数据序列进行插值。由于模拟源信号发生器的采样间隔为10 μs，远小于毫秒脉冲星的周期，因而采用辛格函数SINC插值方法可以获得较为平滑的数据序列。

设标准轮廓已知相位点的时间间隔为T0，时间点tSSBk可表示为tSSBk=mT0+δk，m为整数，δk表示与相邻已知相位点时间偏移量，则tSSBk对应脉冲强度为

为精确计算tSSBk处强度值，实际上不可能取遍无限个点。由于sinc函数在原点具有最大值，而随离开原点距离的增加其函数值衰减很快，所以式（7）中只需取n=m附近较少的点就可得到比较精确的脉冲强度值。

由式（4）有S(τk)=S(tSSBk)，于是时间序列τ0，τ1，…，τk对应的脉冲强度就可以确定了，它对应的是航天器在轨飞行时所接收到的脉冲星信号轮廓。

动态信号脉冲发生器采用模块化设计，主要包括顶层控制、指令通信、高速数据传输、航天器轨道积分、地球星历查询、相位计算、轮廓插值等模块。控制模块对其他各功能模块的运行状态进行管理，各功能模块之间通过控制模块实现数据的交换，独立并行运行。

3 X射线脉冲星导航地面模拟实验

3.1 脉冲星导航静态模拟实验

在X射线脉冲星导航地面静态模拟实验中，脉冲信号发生器根据脉冲星轮廓数据及相位预测模型计算脉冲轮廓数据，脉冲轮廓数据经过栅控信号调制器后作用到栅控X射线管，栅控X射线管发射所需模拟的X射线脉冲，栅控X射线管阳极高压设置为15 kV。X射线脉冲信号经过真空管道后被聚焦型探测系统接收［19］，经过信号处理后提取光子到达时间信息和能量信息，并存储在信号处理计算机中。为了消除脉冲星模拟发射源和探测系统的电子学之间的时钟误差造成的对导航模拟精度的影响，两部分电子学均采用经过GPS驯服的铷原子钟作为基准时钟，其长时间稳定性和短时间稳定性都可高达10−12s/s。同时为了降低模拟源和探测器的同步误差，在每组实验起始时刻栅控信号调制器都会向SDD信号处理电子学发送同步触发信号，使时间记录清零。

实验系统对静态脉冲信号的模拟与真实情况的符合程度主要由轮廓模拟相似度衡量。若实验中模拟的脉冲星为固定周期的，则对采集到的光子到达时间序列以真实周期进行周期折叠得到累积轮廓，比较累积轮廓与真实轮廓即标准轮廓的相似度及脉冲TOA差；若模拟的脉冲星周期是变化的，则通过对采集到的光子到达时间序列进行周期搜索比较测试周期值与理论周期值的误差［20-22］。

分别模拟了在SSB处接收PSR B0531+21、PSR B1937+21两颗脉冲星的X射线脉冲信号标准轮廓数据源自欧洲EPN（The European Pulsar Network）脉冲星数据库［23］，不计脉冲星相位预测模型中频率高阶导数项，即脉冲周期固定不变，分别为33730、1560 μs。每组实验时长45 min，探测器接收的光子计数率控制在300 count/s。

对采集的光子到达时间序列进行周期折叠，设采集时间为N个周期，每个周期分为Nb个bin区间，每个bin区间的大小为10 μs，得出的累积脉冲轮廓如图5所示。累积脉冲轮廓与标准脉冲轮廓之间的相似度用Pearson相关系数表示，两颗脉冲星累积脉冲轮廓与标准脉冲轮廓之间的相似度分别为99.5%、99.1%，均高于99%。

图5 XPVAV地面静态模拟实验结果Fig. 5 XPNAV ground simulation results in static experiment

3.2 脉冲星导航动态模拟实验

在动态模拟实验中，动态信号发生器根据轨道力学模型可以由一个初始轨道状态值推导出后续一定有效时间内的轨道状态，因此只需提供初始时刻的航天器位置和速度坐标。实验开始后，动态信号发生器实时计算包含各种动态效应的变化周期的动态脉冲轮廓幅值数据，再通过栅控信号调制器控制栅控X射线管输出变化周期的动态X射线脉冲信号。探测系统接收到X射线光子信息后存储在信号处理计算机中，通过对采集到的光子到达时间序列进行周期搜索，并消除动态效应后进行脉冲轮廓还原，然后与理论值进行对比从而对实验的质量和模拟相似性进行验证。

实验系统对动态脉冲信号的模拟与真实情况的符合程度主要由脉冲周期测试值与理论值误差、动态效应消除后累积轮廓相似度、脉冲TOA估计值与理论值的偏差衡量。对采集到的光子到达时间序列进行周期搜索，比较周期测试值与理论值的误差；此外，对采集到的光子到达时间序列，先将其转换到SSB处的TCB时间序列，然后根据真实周期进行周期折叠得到累积轮廓，比较累积轮廓与真实轮廓即标准轮廓的相似度和脉冲TOA差。

3.2.1 X射线脉冲星脉冲周期变化及轮廓还原

在近地圆轨道中航天器做匀速圆周运动，接收的X射线脉冲周期是变化的。模拟了航天器在近地轨道飞行一周的情形，初始状态假定为［x0，y0，z0，vx0，vy0，vz0］T=［514306.42 m， −1380476.65 m， 6723879.58 m， −6342.58 m/s， 3984.04 m/s， 1298.94 m/s］T，即航天器轨道半径6.88337×106m，速度7602 m/s，轨道周期5700 s。初始时间为MJD57723，对应地球位置状态为［Ex0，Ey0，Ez0］T=［53400782398 m， 126870439795 m， 54972975401 m］T。分别对航天器在轨飞行时探测器接收PSR B0531+21、PSR B1937+21脉冲星的脉冲轮廓情形进行了模拟，每组实验时长1.5 h。

对每组实验采集到的光子到达时间序列进行周期搜索，得到的脉冲周期变化趋势如图6所示。PSR B0531+21的周期搜索每150 s搜索一次周期，bin的长度为100 μs，周期测试值与理论值偏差的标准差为38451 ps。PSR B1937+21每30 s搜索一次周期，bin的长度为10 μs，周期测试值与理论值偏差的标准差为350 ps。分析两颗脉冲星地面模拟的周期测试值与理论值得偏差是由于信号模拟的栅控信号调制器的产生精度、X射线管长时间工作出射光子的强度稳定度等综合因素引起。在所有引起周期变化的动态效应中，由航天器的运动产生的多普勒效应占主导因素，因此，两组数据的周期变化趋势近似为正弦曲线，表征了航天器绕地球飞行一周时速度在脉冲星方向向量上投影的变化趋势。由于地球的运动同样会产生多普勒效应，会使中心周期（最大与最小周期的平均值）相对于静态周期有一定的偏移，两组数据的偏移量分别为844、−41 ns。

图6 动态模拟实验周期变化趋势测试Fig. 6 Pulse period changes in dynamic experiment

对每组实验数据，分别对观测时间为0～300、0～1800、0～3600、0～5400 s的数据进行周期搜索，然后使用搜索到的周期进行周期折叠得到的累积轮廓如图7所示。随着观测时间的增加，脉冲展宽越明显，轮廓畸变越严重。

图7 动态模拟实验不同时间段轮廓还原Fig. 7 Pulse profile at different time period in dynamic experiment

探测系统采集到的光子到达时间为航天器固有时间，根据实验设定的轨道信息及脉冲星相关参数，先将固有时间转换为TCB时间，再将其转换到SSB处的TCB时间，即可消除动态效应的影响。时间转换的步长为10 μs，即每10 μs计算一次轨道状态，并将该时间间隔内包含的光子事件进行时间补偿，得到其对应的SSB处的TCB时间值。对新的光子到达时间序列进行周期搜索，测试得出的周期值分别为3.373×1010、1.560×109ps，与理论值偏差可以忽略。将新的时间序列进行周期折叠，得出的累积脉冲轮廓如图8所示，累积轮廓与标准轮廓的相似度分别为99.86%、99.99%。

图8 动态实验推算到SSB处的还原轮廓Fig. 8 Pulse profile at the SSB in dynamic experiment

3.2.2 椭圆轨道超实时仿真实验

超实时仿真的主要目的是在保证仿真信息量不变的情况下缩短仿真时间，以满足长时间实验的需要。在本次实验中，模拟航天器在椭圆轨道飞行一周，假定轨道周期P为6 h，根据开普勒第三定律，轨道半长轴a的表达式为

式中：µ为开普勒常数，其值为3.9861×105km3/s2。则半长轴a为1.6763523×107m，为了使实验现象明显，将轨道偏心率设定较大为0.6，则半短轴b为1.3314309×107m，近地点为6.578000×106m。椭圆轨道航天器速度表达式为

式中：r为航天器到地球的距离。根据式（9）计算得近地点速度为9869.98 m/s，若将仿真起始点设置为近地点，则轨道初始状态设置为［x0，y0，z0，vx0，vy0，vz0］T= ［0 m， 6578000 m， 0 m， 9869.98 m/s， 0 m/s， 0 m/s］T。初始时间设定为MJD57723历元，对应地球位置状态为［Ex0，Ey0，Ez0］T=［53400782398 m， 126870439795 m， 54972975401 m］T。采用双倍速率的超实时仿真，即实际的实验时长为3 h。脉冲星选取PSR B1937+21，周期1.56 ms，探测器光子计数率为100 counts/s。

在超实时仿真模式下，探测系统时标会压缩到与脉冲X射线源相同的尺度。对采集到的光子到达时间序列进行周期搜索，bin的长度设置为10 μs，每120 s搜索一次，得出的脉冲周期变化趋势如图9所示。周期测试值与理论值的偏差的标准差为592 ps，中心周期与静态脉冲周期的偏移量主要由地球公转产生的多普勒效应引起，偏移量为−40.4 ns，与图6（b）一致。

图9 航天器为椭圆轨道条件下PSR B1937+21脉冲周期变化趋势Fig. 9 PSR B1937+21 pulse period changes with ellip⁃tical orbit

根据设定的轨道信息及脉冲星初始参数，将采集光子到达时间序列转换到SSB处的时间序列，转换步长为10 μs，得到的累积脉冲轮廓如图10所示，bin的大小为10 μs，与标准轮廓相似度为99.89%。由实验结果可知，地面实验系统可以有效对航天器在椭圆轨道飞行的情形进行模拟，超实时仿真对脉冲周期变化的模拟及轮廓模拟相似度与正常速率仿真模式处于同一水平，仍能保持较高的实验质量。

图10 椭圆轨道消除动态后PSR B1937+21还原轮廓Fig. 10 PSR B1937+21 pulse profile of elliptical orbit with dynamic effect eliminated

3.2.3 轨道机动模拟

在航天器的飞行任务中，轨道调整与机动是经常执行的环节。实验模拟了霍曼变轨的过程，初始轨道为高度是200 km的圆轨道，经过转移轨道将初始轨道高度增加1200 km，到达高度为1400 km的目标轨道，轨道机动示意图如图11所示。起始时刻，航天器在初始轨道的A点处，飞行半周后到达转移轨道的近地点B点处，经过轨道机动后沿着转移轨道飞行到远地点C点，再经过轨道机动后进入目标轨道，最后在目标轨道上飞行半周后到达D点。轨道机动模拟主要通过设置不同轨道状态的轨道根数，从而改变轨道积分模型的输入状态。A点处航天器状态为［x0，y0，z0，vx0，vy0，vz0］T= ［0 m，−6578000 m，0 m，−7784.43 m/s， 0 m/s， 0 m/s］T，机动前B点状态为［0 m，6578000 m，0 m，7784.43 m/s， 0 m/s，0 m/s］T，机动后B点的状态为［0 m，6578000 m，0 m，8103.25 m/s， 0 m/s， 0 m/s］T，其中椭圆轨道半长轴为7.178000×106m，半短轴为7.152879×106m，机动前C点状态为［0 m，−7778000 m，0 m，−6853.07 m/s， 0 m/s， 0 m/s］T，机动后C点的状态为［0 m，−7778000 m， 0 m，−7158.80 m/s，0 m/s， 0 m/s］T。脉冲星选取PSR B1937+21，探测器光子计数率为100 counts/s，A点处起始时刻为MJD57723历元，实验时长2.53 h。

图11 霍曼变轨示意图Fig. 11 Hohmann orbit change diagram

对采集到的光子到达时间序列进行周期搜索，bin的长度设置为10 μs，每60 s搜索一次，得出的脉冲周期变化趋势如图12所示，其中，周期测试值与理论值的偏差的标准差为637 ps。变化曲线中心周期为1.559959232×109ps，曲线与中心周期的交点分别在第0、2655、5681、9094 s，初始轨道、转移轨道、目标轨道的半周期分别为2655、3026、3413 s，因此，中心周期的取值分别在轨道的A、B、C、D点，与轨道模型一致。

图12 霍曼变轨过程中PSR B1937+21脉冲周期变化趋势Fig. 12 Pulse period variation of PSR B1937+21 in the Hohmann orbit transfer

根据以上动态模拟实验结果可知，“天枢Ⅱ号”地面实验系统可以实现对不同类型的航天器轨道模拟以及航天器在轨运动产生的各种动态效应模拟，其中，对动态效应中的短时间缓慢变化过程及长时间缓慢变化过程都能实现高相似度的模拟。因此，基于该地面实验系统，可以实现接近太空实验场景的模拟实验，开展更丰富的动态脉冲数据处理方法的研究以及系统性的导航滤波算法验证研究。基于该地面实验系统，可以为动态信号处理中各种影响因素的研究提供有力的实验数据支撑。