王禹淞,王奕迪,郑 伟

(国防科技大学空天科学学院,长沙 410073)

0 引言

卫星导航系统是以时间测量为基础的系统,精确的位置测量本质上依赖于精确的时间测量[1]。导航卫星的时间参考由其携带的星载原子钟提供,因此,星载原子钟的钟差直接影响了卫星导航系统的定位与导航精度[2-4]。当前,星载原子钟的钟差修正主要依靠地面测控系统测量并将修正量上传至卫星[5]。该方法精度较高且十分稳定,但也存在地面测控系统的负担大、卫星导航系统对地面测控系统依赖程度高的问题[5]。

脉冲星是一类高速自转的中子星,其自转的长期稳定性极佳,被誉为自然界最稳定的天然时钟[6]。由于脉冲星具有良好的长期稳定性,众多学者在脉冲星长期稳定度分析以及脉冲星计时等方面开展了许多研究工作[7-8]。Hobbs等[9]利用国际脉冲星计时阵(international pulsar timing array,IPTA)第一次释放的数据,综合利用48颗毫秒脉冲星长期计时资料,构建了综合脉冲星时。Li等[10]使用昆明40 m 射电望远镜的脉冲星计时观测数据,实现了80 ns量级的本地时钟钟跳检测。Deneva等[11]利用NICER 探测器1 年左右的PSR B1937+21 和9 个月左右的PSR B1821-24观测数据,开展了X 射线脉冲星计时分析。Sun等[12]利用NICER 探测器5年的PSR B1937+21脉冲星的X 射线观测数据进行计时分析,并与利用射电波段观测数据的计时结果进行了对比。除了脉冲星计时外,还可利用脉冲星观测数据进行星载原子钟钟差的估计,从而实现导航卫星自主守时[13]。Sun等[13]分析了脉冲星守时可行性,提出了基于卡尔曼滤波的原子钟钟差估计方法。王奕迪等[14]在考虑脉冲星到达时间测量存在系统误差的情况下,提出了二步卡尔曼滤波算法,有效降低了系统误差对守时性能的影响。童明雷等[15]利用脉冲星实测数据,提出了利用脉冲星校准星载原子钟频率的方法,为脉冲星守时提供了新的思路。

在上述研究中,星载原子钟的钟差、频率等数据均为仿真生成,并不能完全反映脉冲星守时系统的性能。本文利用武汉大学数据中心发布的星载原子钟钟差数据,研究和分析了利用脉冲星对星载原子钟钟差的估计方法,为验证脉冲星守时系统的可行性提供了参考。

1 星载原子钟数据处理与分析

1.1 数据预处理

由于星载原子钟在运行过程中受到环境变化、设备老化等因素的影响,其钟差和频率数据中可能存在粗差或中断[1-2]。因此,在进行星载原子钟性能评估前,需要对原始钟差数据进行预处理。对于数据中断的情况,可采用线性插值进行拟合[2]。对于数据粗差,需要对其进行检测和剔除[1]。通常,粗差的检测一般采用中位数法(median absolute deviation,MAD),检测对象通常为由原始钟差数据转化得到的频率数据。钟差数据转化为频率数据的方法为

式中,xi表示第i个钟差数据;yi表示i个频率数据;Δt为两个钟差数据之间的时间间隔。

将钟差数据转化为频率数据后,将满足以下条件的数据认为是粗差并加以剔除

式中,m为y序列的中位数;n为经验取值,一般为3～5;M=Median)/0.674 5。 此外,对于存在较大频率漂移的数据,需要先去除频率漂移的趋势项,再进行粗差检测。图1为G07卫星(铷钟)剔除粗差后的频率数据和钟差数据。

图1 G07的钟差和频差数据Fig.1 Clock deviation and relative frequency deviation data of G07

1.2 原子钟钟差模型

多项式模型是目前最常用的钟差模型之一,对于一组钟差序列xi(i=1,2,…,N),xi可用多项式模型表示为[16]

式中,a0,a1,a2,…,an为n个模型参数;n为多项式阶数;ei为模型误差。通常,铷钟选用2阶模型,铯钟和氢钟选用1阶模型,本文的研究主要考虑2阶模型的情况。

根据式(3),铷钟的状态方程可表示为[16]

式中,τ为采样时间;x(t),y(t),z(t)分别为原子钟的钟差、钟速和钟飘;εx,εy,εz为均值为0的模型随机误差。

式(4)可简化为

式中,Xk=[x(t+τ)y(t+τ)z(t+τ)]为tk时刻的系统状态;Wk为系统模型误差,其协方差阵Qk可表示为[16]

式中,q1为对应εx的过程噪声参数,表现为调相随机游走噪声;q2为对应于εy的过程噪声参数,表现为调频随机游走噪声;q3为对应于εz的过程噪声参数,表现为调频随机奔跑噪声。

1.3 频率稳定度分析

频率稳定度是衡量原子钟性能的重要指标之一,反映了原子钟的输出频率受噪声的影响而产生的随机起伏的情况[1]。阿伦方差(Allan variance)是一种常用的在时域评估稳定度的方法。对于采样间隔为τ0的频率数据序列{yn,n=1,2,…,L},该序列的阿伦方差表达式为[16]

式中,τ=lτ0表示平滑时间;¯yi(l)=表示第i个平滑时间内的l个频率数据的均值;L'为的个数。本文利用武汉大学分析中心发布的钟差数据(采样间隔为900 s)根据G07卫星自2016年7月1日至2020年6月26日的数据,计算了其阿伦标准差,结果如图2所示,该图为平滑时间和阿伦标准差的双对数坐标图。从图中可以看出,随着平滑时间的增长,阿伦标准差先逐渐减小,后逐渐增大。

图2 基于Allan标准差的频率稳定度评估结果Fig.2 Results of frequency stability based on Allan deviation

阿伦方差还可用于原子钟过程噪声参数的估计,具体方法为[16]

根据式(8),采用最小二乘法,即可计算出过程噪声参数q1,q2,q3。

1.4 原子钟误差特性分析

若真实的原子钟钟差系数已知,则原子钟的钟差仅受随机性误差的影响。因此,首先分析钟差的随机性误差。本文以钟差的功率谱表征钟差的随机性误差水平。根据原子钟的钟差数据,计算了其功率谱,如图3所示,该图为频率和功率谱的双对数坐标图。图中,蓝色线为根据钟差数据计算的功率谱,红色直线为对功率谱进行拟合的曲线,黄色和紫色曲线分别表示其功率谱的上、下包络。

图3 G07的功率谱Fig.3 Power spectrum of G07

为对比脉冲星与原子钟的性能差异,利用NICER 探测器的PSR B1937+21 脉冲星在2017年至2019年的实测数据,计算脉冲星的功率谱,并将其与原子钟的功率谱进行对比,结果如图4所示,该图为频率和功率谱的双对数坐标图。图中,蓝色线表示脉冲星的功率谱,红色线为脉冲星功率谱的拟合曲线,黄色线为原子钟功率谱的拟合曲线,紫色和绿色的虚线为原子钟功率谱的上、下包络。

图4 G07与脉冲星的功率谱对比Fig.4 Compare of the power spectrums of G07 and PSR B1937+21

如图4所示,脉冲星的拟合曲线(红色)与铷钟的拟合曲线(黄色)在约5.3天相交,之后脉冲星的功率谱拟合曲线(红色)低于铷钟的功率谱拟合曲线(黄色)。在约33.2天后,脉冲星实测功率谱(蓝色线)低于铷钟功率谱的下包络。因此,脉冲星的长期性能优于原子钟,可利用脉冲星的观测数据校准原子钟。

图1 给出了G07(铷钟)的钟差和频差数据。可以看出,该铷钟在约3年的时间内,钟差变化了接近800μs。因此,在含系统性误差影响的情况下,星载原子钟的钟差变化很大。

图5给出了PSR B1937+21 的计时残差变化情况,可以看出,在接近3年的时间内,脉冲星的计时残差不存在系统性误差,始终不超过15μs。因此,虽然PSR B1937+21的随机性误差在33天才优于G07号卫星的原子钟,但是由于星载原子钟存在系统性误差,即使在较短的时间内,也可以通过脉冲星来估计原子钟的钟差。

图5 PSR B1937+21的计时残差Fig.5 Residual of PSR B1937+21

2 脉冲星守时方法

2.1 脉冲星守时方案设计

脉冲星守时方案如图6所示。首先,通过星载原子钟测量脉冲星的脉冲到达时间。而后,利用该脉冲到达时间,对星载原子钟的钟差进行估计。从图中可以看出,脉冲星守时的关键,在于如何利用脉冲到达时间,实现原子钟钟差的准确估计。

图6 脉冲星守时总体方案Fig.6 The overall plan for pulsar timing

图7给出了利用脉冲星估计原子钟钟差的方法示意图。首先,利用星载原子钟,可以测量脉冲到达时间。此外,对原子钟的钟差模型进行建模。结合脉冲到达时间和原子钟的钟差模型,利用卡尔曼滤波算法,即可实现原子钟的钟差估计。

图7 利用脉冲星估计原子钟钟差方法示意图Fig.7 Schematic diagram of using pulsar to estimate atomic clock deviation method

2.2 脉冲星测量模型

脉冲星的时间相位模型是实现脉冲星守时的基础,时间相位模型通常建立在太阳系质心(solar system barycenter,SSB)。假设由星载原子钟记录的脉冲到达卫星的时间为tSC,则该脉冲到达SSB的时间为[17]

式中,n表示脉冲星方向矢量;r(tSC)表示卫星相对于地球的位置矢量;rE(tSC)表示地球相对于SSB的位置矢量;c表示光速;Mk表示第k个天体的质量;rk(tSC)表示第k个天体相对于卫星的位置矢量;H.O.T 表示可忽略的高阶项。

在SSB处的脉冲星时间相位模型为[17]

式中,φ表示脉冲相位;φ0表示tSSB,0=g(t0)时刻的脉冲相位;f(n)为tSSB,0时刻的脉冲星自转频率的n阶导数。利用式(10)所示的时间相位模型,即可预测脉冲星信号到达SSB的时间SSB。

令

因此,脉冲到达时间的预测值与测量值之差即反映了原子钟的钟差。因此,脉冲星测量方程可表示为

2.3 脉冲到达时间精度分析

脉冲到达时间的精度可通过以下公式表示[18]

式中,T50是半流量密度持续时间,Tb是敏感器的时间分辨率,B是背景噪声的光子数,S是来自脉冲星的光子数。

式中,Δt是观测持续时间,λp和λn分别是脉冲信号和背景噪声的平均流量密度,A是探测器有效面积。

将式(15)(16)代入式(14),最后可得脉冲TOA测量精度的表达式为

3 守时性能分析

本章以武汉大学发布的G07卫星自2016年7月1日至2020年6月26日的钟差数据为基础,分析脉冲星守时系统性能。假设脉冲到达时间的测量精度为1μs,观测周期为30天。原子钟钟差的过程噪声参数通过式(8)进行拟合。在进行脉冲星守时计算时,首先使用前3个测量值,利用最小二乘法计算原子钟的状态及其状态协方差矩阵。而后,将最小二乘计算的原子钟状态和状态协方差矩阵作为卡尔曼滤波的初值进行解算。

图8(a)给出了脉冲星守时系统的仿真计算结果,可以看出,最初用最小二乘法解算的结果钟差较大,之后脉冲星守时系统精度逐渐收敛至1μs左右。图8(b)给出了利用批量最小二乘的脉冲星守时系统的钟差,可以看出,该守时系统的钟差保持在6μs以内,守时精度低于卡尔曼滤波方法。这是因为,在批量最小二乘算法中未考虑钟差模型的过程噪声,导致基于批量最小二乘的脉冲星守时系统钟差大于基于卡尔曼滤波的脉冲星守时系统钟差。

图8 脉冲星守时系统仿真结果Fig.8 Simulation results of pulsar timing system

对基于卡尔曼滤波和批量最小二乘的脉冲星守时系统进行了500次蒙特卡罗仿真计算,得到基于卡尔曼滤波的守时系统精度为0.78μs,基于批量最小二乘的守时系统精度为1.86μs。可以看出,基于卡尔曼滤波的脉冲星守时系统精度优于基于批量最小二乘的脉冲星守时系统。因此,分析并在钟差估计过程中考虑星载原子钟模型的过程噪声,对提高守时系统的性能是有意义的。

分析脉冲星观测策略和探测器有效面积对脉冲星守时系统性能的影响。表1给出了不同观测策略下的脉冲星守时系统精度和对应需要的探测器有效面积。根据表中前3列可以看出,相同的观测周期下,脉冲到达时间解算精度越高,则脉冲星守时系统精度越高。脉冲到达时间解算精度相同的情况下,观测周期越短则脉冲星守时系统精度越高。在计算探测器有效面积时,假设观测的脉冲星为PSR B1937+21,其半流量密度持续时间为21μs,λp=0.145 photon/(cm2·s),λn=12 photon/(cm2·s),敏感器的时间分辨率为0.1μs。根据表中第3列和第4列可以看出,对于观测同一颗脉冲星的情况,探测器有效面积越大,则脉冲星守时系统的性能越好。增大探测器的有效面积对提高守时系统的性能具有积极作用。然而,探测器的有效面积增大意味着航天器的搭载成本提高。因此,需要根据实际的工程应用场景和任务需求,对探测器有效面积进行设计。

表1 不同观测策略下的脉冲星守时系统性能及探测器有效面积Tab.1 Performance of pulsar timing systems and effective detector areas under different observation strategies

此外,需要说明的是,在本文的仿真分析中,脉冲星观测策略及需要的探测器有效面积均是根据理论公式计算的,与实际空间环境下真实的脉冲星观测性能可能存在一定差距。在后续的脉冲星守时方法研究中,应该结合X 射线探测器的实际观测能力或国内外脉冲星实测数据开展分析与研究。

4 结论

本文利用实测星载原子钟钟差数据和脉冲星观测数据,开展了脉冲星守时方法研究,得到如下结论:

1)通过对比脉冲星计时残差和星载原子钟钟差的功率谱,可以看出脉冲星的长期性能优于星载原子钟。

2)由于包含系统误差的影响,星载原子钟的钟差变化剧烈,G07号卫星的原子钟在约3年的时间内,钟差变化了接近800μs。而脉冲星的计时残差几乎不存在系统误差,PSR B1937+21在接近3年的时间内计时残差始终不超过15μs。因此,虽然PSR B1937+21 的随机性误差在33 天后才优于G07号卫星的原子钟,但是由于星载原子钟存在系统性误差,即使在较短的时间内,也可以通过脉冲星来估计原子钟的钟差。

3)基于实测星载原子钟钟差数据,开展了脉冲星守时性能分析。计算结果表明,若脉冲星的脉冲到达时间解算精度为1μs/30 d,则原子钟钟差估计精度可达到优于1μs的水平,初步验证了脉冲星守时系统的可行性。