应用星载GNSS接收机的Molniya轨道卫星测定轨方法

董召勇 尤超蓝 李文峰

(上海卫星工程研究所，上海 201109)



应用星载GNSS接收机的Molniya轨道卫星测定轨方法

董召勇 尤超蓝 李文峰

(上海卫星工程研究所，上海 201109)

针对传统地基测定轨技术应用于Molniya轨道卫星时跟踪弧段不足和定位精度低的问题，提出了应用GNSS信息的测定轨方法，以满足任务对定轨和预报精度要求。分析了Molniya轨道对GNSS星座的导航信息可用性，研究了应用GNSS导航解的轨道确定及预报方法，仿真计算了两种数据精度模式下导航解作为测量数据参与轨道估算的定轨预报误差，并分析了测量轨道圈数和跟踪模型偏差等因素对定轨预报精度的影响。仿真结果表明：利用2圈的导航解参与运算的定轨误差在15m以内，预报6天的位置误差在60m以内。此方法具有精度高、稳定性好和数据需求量少等优点，对未来大椭圆轨道卫星测定轨工程实践具有借鉴意义。

Molniya轨道；全球卫星导航系统；导航解；精密定轨；最小二乘估计

1 引言

Molniya轨道是一种临界倾角的大椭圆轨道(HEO)，其远地点位于北半球高纬度地区，远地点高度和星下点位置较为稳定，具有长期保持在北半球上空的“逗留”特征[1]，1个周期内90%以上的弧段位于北半球上空。Molniya轨道周期通常为12 h，地面轨迹2圈后回归，在高轨段卫星载荷视场可以完全覆盖北极地区，广泛应用于通信、电子侦察和空间物理探测等任务。

目前，针对HEO卫星的测定轨，比较成熟的是采用地基多站测距与测角结合的方式，将地基测量的距离、方位角、俯仰角信息进行批处理后得到初始轨道状态，进行轨道预报。与天基测量相比，地基测量对卫星轨道的覆盖率等受地面测控网布设不足的影响十分严重，考虑到我国境外测控站部署较少，地基可跟踪弧段和跟踪时间受到极大限制。因此，地基测定轨技术应用于HEO卫星的测定轨时，一般定轨精度不高，卫星位置确定精度仅为几百米。采用星载全球卫星导航系统(GNSS)测定轨技术可克服此不足，并有望实现定轨误差50 m以内、预报6天误差100 m以内的能力。

星载GNSS测定轨技术即在卫星上配置GNSS接收机，以导航卫星作为空间基准点测量卫星的位置、速度等轨道状态信息[2]。出于全球覆盖的考虑，导航卫星空间部署数量较多，能够满足Molniya轨道大弧段覆盖，在提升测定轨精度的同时，可以实现卫星自主定轨的操作。尽管星载GNSS测定轨技术大多在中低轨道应用，应用于高轨和大椭圆轨道时存在信号“稀疏”等问题，但国内外对其在高轨的应用已进行了积极探索，取得了一定成果。

1997年末，“团队卫星”(TEAMSAT)/“年轻工程师卫星”(YES)和“赤道卫星”(EQUATOR-S)[3]上的GPS试验，通过人为地命令接收机去跟踪特定GPS卫星，成功地在高轨段捕获到信号；美国空军研究院主持的Falcon Gold卫星，使用采样接收机在地球转移轨道(GTO)上记录了稀疏的GPS频谱样本[4]。这些试验是卫星在高轨上跟踪GPS信号的最早案例。

无线电爱好者卫星公司(AMSAT)研制的“AMSAT奥斯卡”卫星(AO-40)[5]运行于HEO轨道，卫星上实施的GPS试验，表明在高于GPS星座的58 000 km轨道段也可以追踪到GPS信号，实现无地面介入情况下的全轨段GPS信号自治获取和跟踪。在卫星运行速度达10 km/s的低轨段，单个6通道GPS接收机采用简化的信号获取方式即可追踪到4个及以上的GPS卫星。伪距导出的导航解与伪距直接作为测量数据所确定的轨道，两者之间在3个方向的偏差均值均小于3 km，标准差在几百米量级。

针对高轨和HEO卫星应用中可见GPS卫星不足而难以获得导航解的情况，文献[6]提出了一种HEO轨道确定初始化算法。此算法事先假定标称轨道根数信息，应用收集的伪距和多普勒测量信息搜索出经过近地点的时间，得到历元时刻的状态估计作为滤波初始值，同时提出了一种应用多普勒信息的批处理滤波方法[6]，以计算AO-40卫星的轨道，并成功对无线电爱好者卫星公司的实际飞行数据和戈达德航天飞行中心的飞行仿真数据进行了处理。采用此方法，AO-40卫星沿运行轨迹和垂直运行轨迹方向的位置确定精度分别达到了35 m和87 m。

文献[7]介绍了一种基于GNSS星间链路的中高轨卫星定轨方法，假定用户卫星也能像导航星座内部中高轨卫星一样接收GNSS星间链路信号，此方法对HEO卫星在径向的定轨精度可达6 m以内，切向和法向在2 m以内。目前世界主要卫星导航系统中，只有GPS和GLONASS建设有星间链路，北斗二号卫星导航系统星间链路的建设正在进行。短期内，基于GNSS星间链路的测定轨技术难以应用于我国卫星测定轨的工程实践。

本文提出一种应用GNSS导航信息的Molniya轨道卫星测定轨方法，即将星载GNSS接收机解算出的导航解信息经过批处理滤波改进历元时刻轨道根数而确定卫星轨道，并实施轨道预报。首先，给出了典型的混合GNSS星座和Molniya卫星轨道的主要参数，并建立了Molniya卫星可利用GNSS信息的限制条件的分析模型，据此计算出用户星可用GNSS信息的轨道弧段，然后，对Molniya轨道卫星获取的导航解信息以及轨道确定和预报方法进行了研究，最后，通过仿真算例验证了方法的可行性。

2 Molniya轨道描述及GNSS星座定义

用户星运行于Molniya轨道，其轨道根数：半长轴a=26 561.72 km，偏心率e=0.703 403，倾角i=63.4°，近地点幅角ω=270°，升交点赤经Ω=126.091°，平近点角M=0°，初始历元为2018年7月1日00:00:00(J2000.0)，星下点轨迹如图1所示。

图1 Molniya轨道星下点轨迹Fig.1 Sub-track of Molniya orbit satellite

导航星运行于混合GNSS星座，星座包含地球静止轨道(GEO)卫星、倾斜地球同步轨道(IGSO)卫星和中圆地球轨道(MEO)卫星三类卫星。

3 GNSS星座导航信息可用性分析

GNSS接收机在某时刻同时接收到不少于4颗GNSS卫星的信号时，才能够由伪距测量值求出和用户星位置相关的导航解，进而认为用户星在此时刻有“观测值”(即导航解)。需要开展Molniya轨道

卫星对GNSS星座的导航信息可用性分析，得到用户星可解算出导航解的轨道弧段。

GNSS系统设计之初主要是针对地面用户，为其提供连续、实时、高精度的定位导航服务，故GNSS卫星导航信号发射天线主波束指向地球。影响用户星使用GNSS信息的限制条件[8]有如下4个：

(1)导航星和用户星之间的几何视线不被遮挡。

(2)用户星接收天线在传输的GNSS信号主波束范围内。

Molniya轨道用户星近地点高度为1500 km，远地点高度为38 867.2 km，远地点超过各GNSS卫星的轨道高度，导航星和用户星的几何配置如图2所示，地球遮挡GNSS信号的半锥角为

(1)

式中：Re为地球半径，rGNSS为地心指向导航星的矢量。

导航星的L1频段发射天线主瓣波束宽度半角为γ，则约束条件(1)和(2)可抽象为

(2)

式中：ra=-rGNSS，rb=rHEO-rGNSS，rHEO为地心指向用户星的矢量。

图2 导航星和用户星的几何配置Fig.2 Geometric relationship between navigation satellites and user satellites

(3)用户星接收到的信号功率超过接收机灵敏度阈值。GNSS接收机接收信号的功率Pr与GNSS卫星发射功率、发射天线增益、链路损耗、接收天线增益等诸多因素有关，其计算公式为

(3)

式中：PEIRP为等效各向同性辐射功率，它的定义是卫星天线发送出的功率Pt和该天线增益Gt的乘积，对于GNSS卫星，L1频段上发射C/A码的等效各向同性辐射功率为+28.8 dBW；Ld为卫星通信的天波在传播中受到的自由空间传播损耗，传播距离越大，自由空间损耗越大；Gr为GNSS接收天线的增益，一般为3.0 dB；λ为发射信号的波长；R为GNSS卫星到接收机的距离。

接收机灵敏度阈值Pts，目前星载GNSS 接收机中，采用快速傅里叶变换(FFT)的快速捕获方法，实现了-162 dBW电平下的快速捕获；在高轨卫星上，GNSS接收机接收到的信号非常微弱，其灵敏度一般较高，为-173 dBW(无星历)。

则约束条件(3)可抽象为Pr≥Pts。

(4) 剔除卫星数量不够和相对几何关系不利弧段，精度衰减因子(DOP)值满足要求。当可见GNSS卫星数目过少时，用户星无法通过GNSS星座进行导航定位；当可见GNSS卫星数目足够多但相对几何关系较差时，用户星使用GNSS星座进行导航定位的精度较差。由卫星几何关系造成的伪距测量与用户定位误差间的比例关系，用几何精度衰减因子DOP来表征，在伪距观测精度确定的情况下，需要尽可能减小精度因子的数值，以提升对用户星的定位精度。星座DOP值是影响测量精度的一个重要因素，本文作为仿真研究仅考虑可见GNSS卫星数目的影响，而暂不考虑卫星相对几何的影响。

4 GNSS接收机获取的信息及轨道确定与预报方法

根据接收到的4颗或者更多导航星给出的原始伪距测量值，Molniya轨道卫星搭载的GNSS接收机可解算出用户星的空间位置ru和时钟偏移量tu，具体解算过程见文献[9]。导航解即为GNSS接收机解算出的用户星空间位置，文献[10]指出，导航解实际上可用作伪测量值进行卫星轨道确定，采用导航解而不是伪距的内在优势是其简单的测量模型。针对文章提出的方法，建立了摄动卫星轨道动力学模型和导航解作为测量值的观测模型，并给出了基于GNSS导航解的批处理算法。

应用星载GNSS接收机的轨道确定与预报方法将涉及诸多时空参考系，时间参考系[11]包括TAI，UTC，GMST，UT1，TT和GPS时间；所涉及的空间参考系[11]包括国际天球参考系(ICRS)、国际地球参考系(ITRS)、WGS84大地坐标系、轨道平面坐标系和RTN坐标系。

4.1 测量模型

在轨道估计问题中将导航解作为观测值处理，导航解ru在WGS84参考系表示，与ICRS下卫星位置矢量r的关系为

ru(t)=MU(t)r(t)

(4)

式中：U(t)是考虑岁差、章动、地球自转和极移的从ICRS到ITRS的转换矩阵；M为ITRS大地坐标系至WGS84大地坐标系的转换矩阵。经计算，ITRS和WGS84(G873)只有厘米量级上的误差，两大地坐标系之间的差异可忽略不计。

引入的测量误差通常在RTN坐标系下表述，RTN坐标系坐标原点为卫星质心，R轴为径向(Radial)，与地心到卫星质心连线方向一致；T轴为横向(Transverse)，在轨道平面内与R轴垂直，指向卫星运动方向；N轴为轨道正法向(Normal)，和瞬时角动量方向一致，与R轴、T轴成右手系。假设在任一时刻，RTN坐标系原点都固定在密切轨道各点上，根据RTN坐标系定义，得到以下关系：

(5)

式中：uR，uT，uN是RTN坐标系统的坐标轴方向矢量；r(t)和v(t)是卫星在ICRS下的瞬时位置和速度矢量；i，j，k是ICRS的坐标轴方向单位矢量。

因此，从ICRS到RTN坐标系转换矩阵为

(6)

则导航解所在WGS84参考系到RTN坐标系的转换矩阵为R[MU(t)]-1。

导航解ru对瞬时位置矢量r的偏导数为

∂ru/∂r=MU(t)

(7)

导航解ru对瞬时速度矢量v的偏导数为0。

4.2 轨道动力学模型

摄动卫星问题的运动微分方程在ICRS坐标系下描述为

(8)

式中：μ为地球引力常数；r为惯性系下卫星相对地球质心的位置矢量；f为单位质量的摄动力，包括地球非球形引力摄动、日月引力摄动、太阳光压摄动和大气阻力。Molniya轨道卫星大部分时间处于轨道高度12 000km以上的中高轨段，对卫星轨道影响量级大于10-10km/s2的摄动力仅包括：地球引力场J2,0项摄动、月球引力、太阳引力、地球引力场J2,2项摄动和太阳辐射压(按量级从大到小排序)[11]。尽管卫星近地点高度达到1500 km，大气阻力量级较小，但是对于Molniya这种大偏心率轨道，大气阻力的影响不可忽视。由于太阳光压系数Cr通常无法由材料常数和卫星几何来精确计算，需要作为状态参数在轨道运动方程中进行解算。

4.3 数据批处理算法

状态方程中除地球中心引力项外，还考虑了J2,0、J2,2、日月引力、太阳辐射压和大气阻力等摄动，状态方程表述为

(9)

式中：w(t)为状态模型误差，设为独立零均值的Gauss白噪声。

观测方程为

Yi=G(Xi,ti)+εii=1,…,l

(10)

式中：Yi是各离散的观测时刻ti对应的导航解测量值；v(t)为观测噪声，设为独立零均值的Gauss白噪声。

待估历元记为t0，相应的待估状态量记为X0，x0待估状态量实际值与参考值之差，借助最小二乘估计的批处理算法进行定轨，可得到历元时刻x0的估计值为

(11)

通过轨道确定，一方面得到历元时刻卫星的精确运动状态，另一方面将轨道动力学模型描述得更加精准，可以将其用作预报公式，通过数值积分方法推算出任意时刻卫星的运动状态，以解决GNSS不可用弧段的测轨问题。

5 算例及计算结果分析

假设一种典型的混合GNSS星座进行仿真分析，其组成为：3颗GEO卫星、3颗倾斜地球同步轨道(IGSO)卫星和24颗MEO卫星，合计30颗导航卫星。其中，3颗GEO卫星轨道高度35 788km，分别定点于80°E、110.5°E和140°E；3颗IGSO卫星轨道高度35 788km，倾角55°，星下点轨迹交点分别位于95°E、118°E和147°E附近；24颗MEO卫星轨道高度为20 280km，倾角约55°，近似按照Walker(24/6/1)星座部署于6个轨道面上。各导航卫星星下点轨迹如图3和图4所示。

图3 GEO和IGSO导航卫星星下点轨迹图Fig.3 Sub-track of GEO and IGSO navigation satellites

图4 MEO导航卫星星下点轨迹图Fig.4 Sub-track of MEO navigation satellites

5.1 导航信息可用性计算

当用户星可同时获取4颗及以上的GNSS星信号时，视为在此时刻GNSS信息“可用”。根据轨道特性，各GNSS星和用户星之间的几何关系每隔24 h复现一次，只需分析24 h内的情况，就可以推广到整个用户寿命周期。图5是接收机灵敏度分别定为-162 dBW，-173 dBW时GNSS可用卫星数量随时间的变化情况。

根据图5可知:

(1)对于所研究的Molniya轨道，可见GNSS星数量≥4颗的弧段集中在轨道高度小于10 000 km的低轨部分；

(2)提高接收机灵敏度可同时增加中高轨及低轨弧段的GNSS可见星的数量；

(3)去除零星的短暂可见部分，接收机灵敏度为-162 dBW时，1天内GNSS信息可用时间为3.03 h，可用弧段如图6中加粗部分所示；接收机灵敏度提高到-173 dBW时，1天内GNSS信息可用时间增加至4.24 h。

图5 接收机灵敏度为-162 dBW和-173 dBW对应的GNSS可见星数量Fig.5 Number of visible GNSS stars setting the sensitivity of the receiver as -162 dBW and -173 dBW

图6 接收机灵敏度为-162 dBW时对应的GNSS信息可用弧段Fig.6 GNSS information available arcs setting the sensitivity of the receiver as -162 dBW

5.2 定轨预报精度计算

根据GNSS导航信息可用性分析结果，将每圈测量时段设定为：00:00—00:40，11:20—12:00。对于多圈观测，当前轨道周期的测量时段加上12 h可得下个轨道周期的测量时段，依次类推。每隔60 s取一组数据，一圈可得到80组导航解，利用导航解进行精密定轨计算。

给定两种导航解精度(WGS84下)：

(1) 对于所有轨道高度的弹道点，3个方向(即RTN坐标系下径向、横向和法向，下同)分别添加10 m(3σ)的随机误差，简称10 m/10 m误差模式；

(2) 对于距离地面3000 km以下的弹道点，3个方向分别添加10 m(3σ)的随机误差；对于3000 km以上的弹道点，3个方向分别添加50 m(3σ)的随机误差，简称10 m/ 50 m误差模式。

卫星在轨运行过程中，难以准确获知其卫星结构、方向和反射系数等信息，光压摄动的建模会产生偏差。考虑一个存在20%光压模型偏差(面质比为0.005，光压系数偏差20%)的动力学模型偏差，参与批处理滤波计算。分别用1圈、2圈、3圈、4圈和5圈的导航解参与批处理滤波计算并预报至第7天，得到表1和表2所示的10 m/10 m和10 m/50 m误差模式下的定轨预报误差。

由表1，10 m/10 m误差模式下，只考虑随机差而不考虑模型差时，测量圈数超过1时，定轨误差和预报误差均可保持2 m以内；同时考虑随机差和模型差时，随着测量圈数的增加，定轨误差从6.87 m逐渐增加至27.61 m，而预报误差从72.68 m减小至59.68 m后，又逐渐增加至66.89 m。由表2，10 m/50 m误差模式下，只考虑随机差而不考虑模型差时，随着测量圈数的增加，定轨误差和预报误差均先减小后增大，保持10 m以内；同时考虑随机差和模型差时，随着测量圈数的增加，定轨误差从13.29 m逐渐增加至27.05 m，而预报误差基本不变，维持60 m左右，该方法对测量圈数的依赖性不强。需要说明的是，由于测量数据生成时引入的误差值具有随机性，本文针对同一种工况多次模拟了数据生成和批处理滤波过程，通过蒙特卡罗仿真分析定轨误差分布规律。然而，可能因为同一工况模拟仿真次数不足，未能完全规避误差值的偶然性因素，出现了应用2圈测量数据的定轨和预报精度都略优于3圈的情况。

表1 10 m/10 m测量误差模式下的定轨预报误差

表2 10 m/50 m误差模式下的定轨预报误差

比较表1和表2数据，只考虑随机差而不考虑模型差时，10 m/10 m误差模式对应的定轨误差和预报误差都远小于10 m/50 m误差模式对应的值；同时考虑随机差和模型差时，10 m/10 m和10 m/50 m误差模式对应的定轨误差和预报误差基本相等。测量随机误差对定轨预报精度的影响不明显，模型差的影响比测量误差显著。

图7给出了利用2圈的导航解进行轨道估计计算，两种误差模式下的定轨预报误差曲线。其中，定轨段为0～24 h，预报段为24～168 h。

图7中，实线表示10 m/10 m误差模式对应的定轨预报误差曲线，虚线表示10 m/50 m误差模式对应的定轨预报误差曲线。从图7看出，利用2圈的导航解进行轨道估计计算，两种误差模式下的定轨预报误差在R方向和T方向的变化曲线基本重合，只是在N方向略有差异。

综上所述，应用星载GNSS的测定轨方法的定轨预报精度对测量圈数的依赖性不强，2圈的GNSS导航解数据即可确保较高的定轨预报精度。引入系统误差和模型偏差后，定轨预报精度虽略有下降，但不会急剧发散，说明此方法对测量误差的影响不敏感，稳定性好。利用2圈的导航解参与运算的定轨误差在15 m以内，预报6天的位置预报误差在60 m以内，优于传统地基测定轨方法，满足任务需求。

图7 定轨预报误差变化曲线Fig.7 Orbit determination and prediction error curve

6 结束语

作为一种大椭圆轨道，Molniya轨道偏心率很大，传统地基测定轨方法的卫星定轨误差对测量误差是极其敏感的，测量数据的微小误差就能引起定轨预报的较大误差[12]。应用星载GNSS接收机的测定轨方法将接收机解算出的导航解信息直接作为测量数据进行轨道确定和轨道预报，具有精度高、稳定性好和数据需求量少等优点，实现了卫星自主定轨的操作。与传统地基测定轨方法相比，定轨预报精度从百米量级提升到了几十米量级，对测量误差和模型误差的敏感性也显著降低，此方法具有一定的工程参考价值。

References)

[1]Ilcev S D. Highly elliptical orbits (HEO) for high latitudes and polar coverage[C]// 20th International Crimean Conference on Microwave & Telecommunication Technology.New York: IEEE ,2010

[2]樊士伟,孟轶男,高为广,等. 航天器测定轨技术发展综述[J]. 测绘科学技术学报,2013,30(6):549-554

Fan Shiwei,Meng Yinan,Gao Weiguang,et al. Summarizing on the development of spacecraft orbit determination technology[J]. Journal of Geometrics Science and Technology,2013,30(6):549-554 (in Chinese)

[3]Balbach O,Eissfeller B,Hein G W,et al. Tracking GPS above GPS satellite altitude: first results of the GPS experiment on the HEO mission equator-S[C]//Position Location and Navigation Symposium. New York: IEEE,1998:243-249

[4]Silva R,Chao C C,Sedlacek S B,et al. GPS signals in a geosynchronous transfer orbit: Falcon Gold data processing[R]. Colorado: United States Air Force Academy,1999

[5]Davis G,Moreau M,Carpenter R,et al. GPS-based navigation and orbit determination for the AMSAT AO-40 satellite[R]. Maryland: Goddard Space Flight Center,2002

[6]Axelrad P,Speed E. Algorithms for autonomous GPS orbit determination and formation flying: investigation of initialization approaches and orbit determination for HEO[R]. Colorado: University of Colorado,2002

[7]孟轶男,樊士伟,李罡,等. 利用GNSS星间链路对中高轨航天器测定轨的可行性分析[J]. 武汉大学学报信息科学版,2014,39(4): 445-449

Meng Yinan,Fan Shiwei,Li Gang,et al. Orbit determination of medium-high earth orbital satellite using GNSS crosslink ranging observations[J]. Geomatics and Information Science of Wuhan University,2014,39(4):445-449 (in Chinese)

[8]Li T,Liu J,Huang Z,et al. Observability of HEO satellite autonomous navigation system using GPS[C]// 2010 International Conference on Multimedia Technology. New York: IEEE,2010

[9] Kaplan E D,Hegarty C J. Understanding GPS:principles and applications,2nd edition[M]. London: Artech House,Inc.,2006

[10]Carter S S,Cefolo P J,Proulx R J. The determination of precision mean element sets from GPS receiver on-board navigation solution [C]//Proceedings of the AAS/AIAA Astrodynamics Conference.California: Univelt,Inc.,1995

[11]Montenbruck O,Gill E. Satellite orbits: models,methods,and applications[M]. Berlin: Springer,2000

[12]Antreasian P G,Baird D T,Border J S,et al. 2001 Mars Odyssey orbit determination during interplanetary cruise[J]. Journal of Spacecraft and Rockets,2005,42(2):394-405

(编辑：李多)

Precise Orbit Determination Method of Molniya Orbit Satellites Using Onboard GNSS Receiver

DONG Zhaoyong YOU Chaolan LI Wenfeng

(Shanghai Institute of Satellite Engineering,Shanghai 201109,China)

Considering the problems of lacking visible arc and the low positioning accuracy using the traditional ground-based orbit determination technology,the method using onboard GNSS information is proposed in this paper in order to meet the orbit determination and prediction accuracy from the space mission requirements. The navigation information availability of GNSS constellations is analyzed and the algorithm of orbit determination based on GNSS navigation solutions is studied,then the orbit determination and prediction errors are calculated taking navigation solutions of different accuracy as measurement data. The influence of measured orbital circles and the accuracy of tracking model to orbit determination and prediction accuracy is analyzed. The simulation results show that the orbit determination errors can be restricted below 15 meters and after six days’ prediction the position errors can be restricted below 60 meters using navigation solutions of two orbit circles,when taking the random errors of measurement data and the model deviation into consideration. The method has advantages of high precision,strong stability and low data demand,which may be used for future HEO satellite orbit determination engineering.

Molniya orbit; GNSS; navigation solution; precise orbit determination; least square estimation

2017-04-18；

2017-05-10

董召勇，男，硕士，工程师，研究方向为卫星姿轨控设计。Email：dongzhaoyong1990@163.com。

V

A

10.3969/j.issn.1673-8748.2017.03.004