相对论体系下卫星原时及星间钟差机理建模

贾前，李青，刘磊

（1.西北工业大学航天学院，西安 710072；2.陕西省空天飞行器设计重点实验室，西安 710072；3.清华大学航天航空学院，北京 100084）

0 引言

星载原子钟是卫星导航系统的核心部件，其性能决定了导航系统时间基准的精度［1-2］。目前，地面光学原子钟的稳定度即将突破10-19量级［3-4］，JPL 的星载汞离子钟已通过在轨验证，比目前应用的星载原子钟精度提高了1 个量级以上，有望应用于新一代卫星导航系统［5］。然而，空间中变化的引力场以及卫星的运动会导致星载原子钟计量的原时产生偏差［6］，不同卫星的原时偏差导致了星间钟差，影响对原子钟性能极限的利用［7］。新一代导航卫星均计划或已经通过确定星间钟差来建立与维持星座自主时间基准［8］，以实现导航系统的自主运行，但卫星原时及星间皮秒钟差机理尚不清晰。为满足新一代导航系统对高精度自主时间基准的需求，亟需开展卫星原时偏差以及星间皮秒钟差模型研究。

在卫星原时偏差方面，国内外学者基于卫星钟的物理模型以及相对论效应研究了多种模型。基于卫星钟的物理模型主要利用原子钟的数据进行模型参数求解，进而得到钟差数据，精度在ns 量级［9］。常用的模型有多项式模型［10］、灰色模型［11］、谱分析模型［12］、时间序列模型［13］以及组合模型［14］等。文献［15］提出了二次多项式拟合预报模型，基于卫星钟的物理特性拟合得到卫星钟差，文献［9］在二次多项式模型基础上增加了显著周期项，并结合BP神经网络设计了超快速钟差预报算法实现流程。然而上述模型需要卫星的精密钟差数据与原子钟的随机先验信息，且精度无法满足新一代导航系统分米级导航定位需求［9］。基于相对论效应的模型通过研究卫星原时速率与引力场及卫星运动的关系得到卫星原时偏差。文献［6］与文献［16］在开普勒轨道假设下研究了GPS 和北斗卫星导航系统（BDS）卫星原子钟的相对论效应，并指出周期性累积钟差需要根据卫星星历实时进行计算和改正，但未考虑摄动对星载原子钟的影响。文献［17］以地球同步轨道卫星为研究对象分析了轨道摄动引起的钟差，文献［18］估计了地球J2 项摄动引起Galileo 卫星的原时偏差幅值，文献［19］建立了卫星原时与坐标时转换模型，并指出J2 项摄动引起的原时偏差对于中轨道卫星达到了纳秒量级。根据文献［6］和文献［16-19］的研究可知，基于广义相对论的钟差模型可实现的精度更高，但仍需要在开普勒轨道假设的基础上引入J2 项修正以达到皮秒级精度。星间钟差主要基于星间链路获得的观测数据进行计算，目前双向测量工作体制的星间链路已经用于高精度的钟差测量［20］。文献［21］提出了一种利用星间伪距拟合多项式和钟差拟合多项式联合求解星间钟差算法，精度在5ns 以内。文献［22］针对北斗2 导航系统的星间链路建立了同步时分双工（STDD）体制下星间钟差测量模型，精度优于1ns；文献［23］研究了STDD体制下星间钟差的测量误差模型，并与GPS 采用的轮询时分双工体制（PTDD）进行了对比，验证了其性能。文献［21-23］所采用的方法主要以激光或者微波为信号载体，通过计算两终端发射接收信号的时间差来确定钟差，精度取决于不同路径中延时修正模型，且未从原时模型角度研究星间皮秒钟差机理。文献［19］在开普勒轨道的假设下利用卫星原时模型建立了星间钟差模型，能够有效计算星间钟差，但该模型并未考虑摄动加速度，误差会随时间增加而增大。

本文从广义相对论出发，提出一种考虑J2 项摄动的卫星原时和星间钟差模型，根据卫星的位置速度信息计算得到卫星原时偏差以及星间钟差，并通过仿真对模型进行验证，为星载原子钟的校正和星间钟差确定提供了一种可行的建模方法。

1 坐标系与卫星运动方程

1.1 坐标系描述

为建立卫星运动与时间的联系，首先借助广义相对论建立卫星的运动方程。运动的描述离不开相应的坐标系，本文选择地心参考坐标系（Geocentric coordinate reference system，GCRS）建立地球引力场中卫星的运动模型。GCRS的原点在地球质量中心，并且由四维坐标X≡(ct，x，y，z)T所定义，其中c为真空中的光速，值为2.997 924 58×108m/s，t为坐标系对应的坐标时，x，y，z为卫星的空间位置坐标。在GCRS中，两个确定事件的间隔微分ds可以表示为：

式中：gμυ为GCRS 的度规张量；xμ，υ(μ，υ=0，1，2，3)为GCRS中的时空坐标分量，μ，υ=0时，x0=ct，μ，υ=1，2，3分别对应卫星的三轴坐标分量x，y，z。

度规张量可以由一组谐波势表示，谐波势包含标量势w与矢量势wE，保留至c-4项，选取度规张量为［24］：

式中：α，β和λ为索引标号，从1取至3；γ=diag(-1，-1，-1) ;为wE的第λ个分量。

矢量势wE的表达式为：

式中：G为引力常数，值为6.673 5×10-11N·m2/kg2；ME为地球的质量，为5.974 2×1024kg；r=‖r‖2为卫星地心距；r=[x，y，z]T为卫星在GCRS 中的位置矢量；SE=[0，0，9.8 × 108]Tm2s为地球单位质量的角动量；“(·)×”符号表示矢量叉乘运算的反对称矩阵：

地球引力场中，标量势w由地球引起的引力势UE和其他天体引起的潮汐势组成，表达式为：

考虑地球引力势至J4项，UE的表达式为：

式中：RE为地球平均赤道半径，取值为6 378 km；J2，J3和J4为地球的二至四阶带谐系数，值分别为1.082 6×10-3，-2.532 7×10-6和-1.619 6×10-6；为简化描述，定义UJ为地球的二至四阶带谐项引起的引力势变化的总和。

式中：Mb为天体b 的质量，太阳质量MS为1.980 4×1030kg，月球质量MM为7.336 9×1022kg；rbE为天体b与地球之间的距离；nbE=rbErbE为rbE的单位向量；符号S和M分别代表太阳和月球。

1.2 卫星运动方程

考虑卫星运动的相对论效应时，可以在牛顿体系下卫星的运动方程中添加相应的相对论修正项，也可以利用度规张量直接在广义相对论体系下建立卫星的运动学方程，本文选择第2种方式。

根据式（1）以及选择的度规张量，可以得到两事件的时空间隔微分为：

卫星的运动方程以坐标时t为变量，因此引入拉格朗日乘子L=dsdt，根据式（8）可得拉格朗日乘子的平方表达式为：

式中：v=为卫星运动速度大小；vx=dxdt，vy=dydt，vz=dzdt分别为卫星各轴向的速度大小。

对式（9）进行开方并保留至c-4项，可得拉格朗日乘子L的表达式为：

拉格朗日乘子L应满足Euler-Lagrange方程：

将式（10）代入式（11），并引入非引力项以及日地相对运动引起的加速度项［24］，推导可得广义相对论体系下卫星的运动方程为：

式中：n=r r为卫星位置r的单位向量；rSE为太阳相对地球的位置向量，且为非引力项产生的卫星加速度。

2 考虑J2项修正的卫星原时以及星间钟差模型

2.1 卫星原时模型

卫星在轨运行时，星载原子钟记录的时间称为原时τ，原时速率会受到卫星运动与引力场的影响发生改变，而坐标时仅与选定的坐标系相关［16］，因此可通过建立坐标时与原时的关系得到卫星的原时模型。在卫星附近的无穷小区域内，度规张量为常值，且原子钟相对卫星静止，即dx=dy=dz=0，因此可得：

联立式（8）与（13）可得卫星原时模型为：

式中：δ(r，v)为由于卫星运动与引力势的存在而产生的原时速率偏差。假设坐标时为导航系统的参考时间，此时原时偏差即为对应的卫星钟差，如图1所示。

图1 引力场中原时偏差示意图Fig.1 Schematic of proper time deviation in the gravitational field

仅考虑地球引力势中的牛顿主项，忽略其他天体引起的潮汐势，即标量势w简化为：

在地球引力场内，卫星所在位置处的引力势主要为地球产生的影响，并且J2 项的量级远大于高阶项，因此在引力势中引入地球非球形J2 项摄动，得到考虑J2项修正的卫星原时模型为：

卫星的状态也可以通过轨道六要素来表示，将式（17）中对应的变量转换为轨道六要素，即可得到基于轨道六要素的卫星原时模型。

卫星地心距r的表达式为：

式中：a为卫星所在轨道的半长轴；e为轨道偏心率；E为偏近点角。

卫星的z轴位置分量可表示为：

式中：i为轨道倾角；f为真近点角；ω为近地点幅角。

卫星的地心距与速度大小满足：

联立式（18）和（20）可得：

在轨道偏心率较小的情况下，根据轨道摄动理论，考虑J2项摄动，相关轨道参数可以表示为：

式中：a0，i0，e0，ω0和E0分别为不受摄动影响的卫星轨道对应的半长轴，轨道倾角，偏心率，近地点幅角以及偏近点角。

将式（18），（19），（21）以及式（22）联立，忽略高阶小量，代入式（17）可得到基于轨道六要素的卫星原时模型为：

根据式（23）可以看到，卫星运行过程中原时与坐标时产生的差异由四部分组成，前两项表示仅考虑牛顿主项加速度时卫星在轨运动引起的相对论效应，后两项为J2项摄动引起的相对论效应。

2.2 星间钟差模型

根据式（17）可以看到，当卫星的状态不同时，卫星原时相对参考时间的差异不同。而导航系统通常由多颗卫星组成，这些卫星分布在不同轨道面或者同一轨道面的不同位置处，因此导航卫星的状态间存在差异，在轨运行时不同的卫星原时会产生对应的星间钟差，如图2所示。

图2 星间钟差示意图Fig.2 Schematic of inter-satellite clock difference

考虑导航系统中两颗卫星A和B，假设在坐标时t0时刻，两颗卫星的星载原子钟读数即原时分别为τA0和τB0，初始钟差为Δτ0=τB0-τA0，则在t1时刻两颗卫星的原时可以表示为：

进一步对两颗卫星的原时作差，可以得到卫星A和B的星间钟差为：

式（25）即为考虑J2 项修正的星间钟差模型，该模型分为四部分：（1）卫星初始钟差；（2）卫星地心距差异引起的星间钟差；（3）卫星速度差异引起的星间钟差；（4）J2项摄动引起的星间钟差。

3 仿真分析

3.1 卫星轨道误差仿真

北斗卫星导航系统（BDS）是我国自主建立的全球卫星导航系统，可在全球范围内全天候、全天时为各类用户提供高精度、高可靠定位、导航、授时服务。本节以BDS 中北斗三号MEO-01 卫星为研究对象，对建立的卫星运动方程进行仿真分析，非引力项加速度考虑为太阳光压引起的加速度，计算方式参考文献［19］。仿真中卫星的初始轨道要素如表1所示，仿真时长为两个轨道周期，仿真步长为1 s，仿真时为坐标时。

表1 北斗三号-MEO-01卫星初始轨道要素Table 1 Initial orbital elements of BeiDou-3 MEO-01 satellite

考虑式（12）右端的全部加速度项时，摄动加速度引起的卫星位置误差以及速度误差分别如图3和图4 所示，两者均随时间的增加而增大。对于MEO-01 卫星而言，在两个轨道周期时，位置误差已经达到了25.48 km，速度误差达到了3.78 m/s，根据卫星原时模型可知，这些误差均会导致原时与坐标时之间产生偏差。

图4 摄动加速度引起的速度总误差Fig.4 Total velocity error induced by perturbation acceleration

仅考虑J2 项摄动加速度时，MEO-01 卫星的位置误差以及速度误差如图5 和图6 所示。根据仿真结果可知，J2 项摄动加速度引起的位置误差与速度误差与模型考虑的全部摄动加速度引起的误差变化趋势相同。两个轨道周期时，J2 项摄动加速度引起的位置误差为21.11 km，速度误差为3.22 m/s。与图3和图4仿真结果对比可知，在考虑的摄动加速度中，J2 项摄动加速度对MEO-01 卫星的影响最大，并且对于地球引力场中的卫星，J2 项的影响会随着轨道高度的减小而增大，因此，对于高精度时间基准的建立，J2项摄动加速度引起的原时偏差不可忽略。

图5 仅考虑J2项摄动加速度时位置误差Fig.5 Position error only considering J2 perturbation acceleration

图6 仅考虑J2项摄动加速度时速度误差Fig.6 Velocity error only considering J2 perturbation acceleration

3.2 卫星钟差仿真

本节根据建立的卫星原时模型对卫星钟差进行仿真分析。仿真中利用卫星运动方程（12）得到卫星状态，将状态信息代入式（14）计算得到的原时偏差并将其作为实际钟差，与考虑J2 项修正的卫星原时模型（17）得到的钟差进行对比，进而分析模型精度，仿真参数与3.1节相同。

实际钟差、考虑J2 项修正的模型（17）以及简化模型（16）的计算结果如图7 所示。根据仿真结果可知，卫星原时与坐标时的偏差即卫星钟差随着时间的增加而增大，在偏心率接近0时，卫星钟差与坐标时之间近似为线性关系，与文献［19］结论一致，在两个轨道周期时，MEO-01 卫星的钟差已经达到了2.21×10-5s。

图7 卫星钟差计算结果Fig.7 Satellite clock difference calculation results

两种模型的误差如图8所示，随着时间的增加，两种模型的计算误差都逐渐增大。在两个轨道周期内，考虑J2项修正的模型误差峰值为26.97 ps，简化模型的误均差峰值为47.95 ps，两者比实际钟差小6个数量级，即两种模型都能用于计算卫星钟差。但由于简化模型未考虑J2 项摄动对卫星原时速率的影响，计算得到的钟差波动大，误差的均方根为27.23 ps，而考虑J2 项修正的模型计算结果波动较小，误差的均方根为16.27 ps，提高了40.26%。仿真结果表明考虑J2 项修正的卫星原时模型可以有效提高卫星钟差的计算精度。

图8 原时模型误差Fig.8 Proper time model error

模型（23）为卫星钟差的轨道要素表示形式，该形式表明J2 项引起的钟差中包含常值项与周期项，轨道倾角为0°时，常值项达到最大值，此最大值与轨道半长轴成反比；当轨道倾角为54.73°时，常值项为0。周期项的周期为卫星轨道周期的一半，幅值与轨道半长轴及轨道倾角有关。在小偏心率情况下，真近点角可通过轨道平均角速度n近似为：

对周期项积分可得：

周期项幅值与轨道半长轴和轨道倾角的关系如图9所示，随着轨道倾角的增加，J2项摄动加速度引起的周期项钟差逐渐增大，并随着轨道半长轴的增加而快速减小。

图9 J2项摄动加速度引起的周期项钟差幅值与轨道倾角及半长轴关系Fig.9 Relationship between the amplitude of periodic proper time deviation caused by J2 perturbation acceleration，orbital inclination，and semi-major axis

3.3 星间钟差仿真

3.2 节的卫星原时仿真结果表明，卫星在轨运行期间星载原子钟会产生钟差，并且钟差与卫星轨道相关，进而导致不同卫星间产生星间钟差。本节对两卫星间的钟差进行仿真研究。文献［19］的研究表明两卫星轨道半长轴不同时，星间钟差分为线性与周期项两部分，线性部分会随着半长轴差异的增加而迅速累积，但通常卫星发射前都会根据目标轨道对原子钟进行频率校正［16］，因此本节只研究卫星轨道半长轴相同的情况下星间钟差的变化规律，且假设初始时刻星载原子钟均经过校准，即Δτ0=0。与卫星钟差仿真相同，仿真中两颗卫星的状态由运动方程（12）得到，卫星原时通过式（14）计算，作差后得到实际星间钟差，与式（25）计算得到的星间钟差进行对比，分析模型误差。

仿真中一颗卫星轨道参数与表1 相同，另一颗卫星初始真近点角为180°，即两卫星轨道面相同，在轨道平面内初始位置相差180°。两卫星间的钟差以及模型误差如图10 和图11 所示。仿真结果表明，两个轨道周期内，模型的计算结果与理论钟差基本吻合，星间钟差呈周期性变化，周期与卫星轨道周期相同，仿真中两颗卫星的钟差峰值为5.89 ns；模型误差也近似为周期变化，最大误差为8.86×10-2ps。

图10 星间钟差Fig.10 Inter-satellite clock difference

图11 星间钟差模型误差Fig.11 Inter-satellite clock difference model error

为说明J2 项修正的有效性，采用文献［19］的模型对相同参数的卫星进行星间钟差仿真，模型误差如图12所示。仿真结果表明，无J2项修正的星间钟差模型误差随时间逐渐增大，在两个轨道周期时，误差为16.37 ps，远大于本文提出的考虑J2 项修正的星间钟差模型的误差。

除了中轨（MEO）卫星，BDS 星座还包含高轨同步（GEO）卫星及高轨倾斜（IGSO）卫星。美国的GPS星座由近圆形轨道（e=0.02）的中轨卫星组成［6］，选取对应轨道上的一颗卫星，另一颗卫星与其在同一轨道平面上，初始位置相差180°，根据星间钟差模型计算这两颗卫星间的钟差，卫星的轨道参数以及计算结果如表2 所示。根据表2 的仿真结果可以看到，四种轨道平面中星间钟差峰值从小到大依次为BDS-GEO，BDS-MEO，BDS-IGSO 和GPS-MEO。这是因为随着轨道偏心率的增加，卫星速度与位置的变化引起的原时偏差增大，进而引起星间钟差峰值变大。对于BDS-IGSO 和GPS-MEO 这两个轨道平面，由于其半长轴、偏心率和轨道倾角的影响，两个轨道平面上卫星受到的太阳及月球的影响较大，而本文建立的星间钟差模型（25）忽略了其余天体的影响，因此模型误差的峰值大于其余两个轨道平面的误差峰值，但仍在皮秒量级。这些算例的计算结果表明考虑J2 项修正的星间钟差模型误差最大在皮秒量级，能够准确计算星间钟差。

表2 不同轨道的星间钟差计算结果Table 2 Calculation results of inter-satellite clock difference in different orbits

4 结论

本文针对新一代导航卫星对高精度时间基准的需求，在广义相对论体系下建立了卫星的运动方程，考虑地球引力场内卫星受到的J2 项摄动影响，提出了包含J2 项修正的卫星原时模型以及星间钟差模型，模型能够有效地计算出卫星在引力场中运动时产生的原时偏差以及星间钟差。相比基于开普勒轨道假设得到的模型，考虑J2 项修正可以提高卫星原时以及星间钟差的计算精度。仿真结果表明，卫星运动过程中产生的原时偏差与卫星轨道要素紧密关联，并随着时间增加而增大，考虑J2 项修正计算得到的原时偏差与实际钟差更为接近。同时本文提出的星间钟差模型能够在皮秒量级的精度上计算星间钟差。本文工作可为导航系统的时间基准建立与自主维持提供理论支持。