李 笛，孙腾达

(卫星导航系统与装备技术国家重点实验室，河北 石家庄 050081)

0 引言

卫星导航模拟器[1]可以通过数据仿真产生在接收机各种运动状态下的导航信号，可以直观地呈现卫星导航系统卫星的空间分布与可见性信息，用户的位置、速度等轨迹信息和接收机接收观测量信息等。因此，卫星导航模拟器广泛用于各种应用场景中接收机的调试与测试[2]。低轨卫星广泛用于探测、通信和导航增强等[3]，发挥着越来越重要的作用，在低轨卫星任务的地面仿真中，卫星导航模拟器发挥着不可替代的作用。其中，低轨用户轨迹长时间、高精度仿真是模拟器正确仿真观测量的基础，其仿真精度直接影响低轨用户数据仿真的逼真度。由于低轨卫星受力情况非常复杂，实现其高精度轨迹仿真特别困难，同时在模拟器中运行还需考虑工程可实现性，针对上述问题，开展低轨用户高精度轨迹仿真方法研究，并兼顾精度与计算量，提升工程的可实现性。

1 低轨用户动力学模型

低轨用户轨迹分布取决于空间复杂的力学环境[4]，故建立相对精确的动力学模型[5]对低轨用户轨迹仿真至关重要。低轨用户在轨运行期间，除了受地球中心引力[6]作用外，还将受到地球非球形引力[7]，太阳、月球及其他天体引力的影响，以及太阳光压、大气阻力和地球潮汐力等因素的影响。

1.1 重力场模型

在初步分析航天器运动规律时，把地球看作理想球体，其对航天器产生指向地心的、大小与航天器质量m成正比而与地心与航天器距离r的平方成反比的引力F[8]：

(1)

由此得到地球引力常数：

μ=GM=3.986 000 441 5×1014m3/s2。

(2)

于是引力加速度为g=μ/r2，此情况下地球引力场称为中心引力场。引力加速度为引力势Uc的梯度，即Uc=μ/r。

地球并不是完美球体，具有扁平度、梨形和赤道变形等，因此，引力加速度是地心经度φ、纬度λ和r的函数。则地球引力加速度矢量g为地球引力势函数U的梯度：

g=gradU(r，φ，λ)。

(3)

目前引力势函数的普遍表达[8]为：

([Cn，mcos(mλ)+Sn，msin(mλ)]}，

(4)

式中，RE=6 378.136 3 km为地球赤道半径；Pn，m(x)为n阶m级的Legendre函数。为克服各阶数系数相差太大而引起的计算误差，采用的范化式[6]为：

(5)

(6)

式中，

(7)

地球引力势模型系数Jn(即-Cn，0)，Cn，m，Sn，m由测量得到。许多国际组织通过不同手段测量重力场系数并创建对应的重力场模型，比较著名的重力场模型包括WGS系列、GEM系列和GRIM系列模型。WGS和GEM模型已经合并为EGM96模型，IER2003推荐采用EGM96模型。

1.2 大气阻力模型

大气对低轨用户的作用包括三部分[6]，与用户对大气流相对速度方向相反的大气阻力、升力和副法线作用力，其中大气阻力占支配地位，其余2种可以忽略。对低轨用户而言，大气阻力为其受到的最大的非引力性摄动力。

低轨用户的大气阻力加速度[6]为：

(8)

式中，r为低轨用户位置矢量；CD为大气阻力系数；A为低轨用户迎风面；ρ为低轨用户处的大气密度；vr为用户对大气流的相对速度的大小；ev为相对速度的单位方向向量。

在大气密度[9]计算方面，可采用目前最常用的3种经验大气密度模型，即Jacchia系列模型、DTM系列模型和MSIS系列模型。由于大气变化机制复杂，任何经验大气密度模型都很难精确地刻画大气密度的实际变化。因此，为进一步提高大气模型精度，可利用实际数据对经验模型进行进一步校准。同时，由于不同轨道高度大气密度不同，所以针对不同轨道类型的低轨卫星，还需研究不同经验大气密度模型下的轨道预报精度，探索选取最优的大气密度计算方法。

在大气阻力系数[10]计算方面，一般而言，在低轨卫星精密定轨中，常常将大气阻力系数作为未知量与卫星的运动状态矢量一起解算，解算得到的大气阻力系数会吸收模型误差，在轨道预报中利用解算得到的大气阻力系数来计算大气阻力加速度，但大气阻力系数的较难准确地确定是制约解算大气阻力加速度精度的一大原因。基于目前已有研究，可采用线性回归分析建立大气阻力系数的补偿算法，还可基于空间环境数据和神经网络模型对空间目标大气阻力系数进行修正。同样，针对不同轨道类型的低轨卫星还需研究不同大气阻力系数算法下的轨道预报精度，探索选取最优的大气阻力系数计算方法。

在低轨卫星迎风面积[11]计算方面，需对低轨卫星的几何形状进行描述，可采用Cannon-Ball模型和Macro模型。Canno-Ball模型是假设低轨卫星为一个球体，具有固定的形状及单一的表面性质，对几何形状较为简单的低轨卫星十分适用。Macro模型也被称为Box-Wing模型，它会对低轨卫星每个面的几何形状都进行描述，包括面积、单位外法向量以及该面的光学特性(对可见光和红外光的反射和吸收)，特别当低轨卫星在飞行过程中太阳帆板发生翻转时，Macro模型亦能做出几何形状的描述并反映到阻力计算中。

1.3 太阳光压模型

照射到卫星上的太阳光对卫星产生入射作用力，卫星吸收一部分太阳光转变成热能和电能，另一部分被反射回去[12]。入射力和反射力统称为太阳辐射压力，也称太阳光压力。太阳辐射压力产生的加速度与太阳光强度、垂直于入射方向的有效面积、表面反射率、与到太阳的距离和光速有关。由于卫星表面材料的老化、太阳能量随太阳活动的变化以及卫星姿态控制的误差等因素的影响，使得太阳辐射压摄动也是难以精确模型的摄动力。

直接的太阳光压力[6]可以表示为：

(9)

式中，P⊙为太阳常数，等于4.560 4×10-6N/m2，其物理意义是完全吸收的物体在距太阳一个天文单位(AU)处，单位面积上所受的辐射压力；CR为太阳光压系数；A为垂直于卫星与太阳方向的截面积；ν为蚀因子(地影区ν=0，卫星光照区ν=1，阴影区0<ν<1)；r⊙为太阳的位置矢量。

1.4 N体摄动模型

太阳、月球和行星的引力摄动可近似用点质量来描述，在地心惯性系中，天体对卫星的N体引力加速度为：

(10)

式中，GMi为质点的引力常数；ri为质点在地心惯性系中的位置矢量。

2 坐标系及其转换

为了减少牵连运动所引起的附加速度，使卫星的运动方程相对简化，一般在惯性坐标系[5]中建立和求解卫星的运动方程。而地球对卫星的引力是在地固坐标系[13]中定义的。此外，卫星摄动力分析计算要涉及到卫星轨道坐标系，把加速度转换到惯性坐标系或地固坐标系要涉及卫星固定坐标系。因此要建立卫星运动方程和参数估计的观测方程，必须给出这些坐标系统的定义及其相互转换关系。

2.1 坐标系统定义

2.1.1 协议惯性坐标系

具有绝对意义的惯性坐标系统是无法获得的，人类只能根据科学任务的需要和一些约定建立相应精度的惯性坐标系统。通常使用的惯性坐标系统是协议惯性系[14](CIS)，其具体实现就是国际协议天球参考框架(ICRF)。在CIS下，物体的运动可以以很高的精度满足Newton力学定理，这对于动力学低轨轨迹仿真的相关研究十分重要。

2.1.2 协议地球坐标系

协议地球坐标系地固坐标系统[14](CTS)是为了描述地面观测站位置所建立的与地球体联结的坐标系统，依表达方式又可分为直角坐标系及大地坐标系。CTS是使用起来最直观的坐标系统，它以特定的协议与地球固联，并由一组全球参考站维持。

2.1.3 卫星轨道坐标系

为了便于对特定摄动力进行分析，将作用于低轨用户的力分解到卫星轨道坐标系统[5](RTN)下进行分析。此外，在RTN系统下，还便于对各项误差进行分析，如利用HL-SST模式研究恢复地球重力场时，径向的轨道误差占非常重要地位。

2.1.4 卫星固定坐标系

卫星固定坐标系[4](SBF)的主要作用是定义卫星在惯性空间中的姿态，同时建立各相关传感器坐标系与惯性系的关系。如加速度计是固联在卫星质心处，其轴系相应地平行于卫星固联坐标系，在转换加速度计观测量到CIS或CTS时都需要使用SBF 作为中间过渡坐标系统。

2.2 坐标系统转换方法

2.2.1 协议地球坐标系到协议惯性坐标系的转换

如果以RCIS表示某点在历元J2000.0对应的CIS中的坐标，以RCTS表示CTS中的坐标，则有[8]：

RCIS=P(t)N(t)S(t)Pm(t)RCTS，

(11)

式中，P(t)，N(t)，S(t)，Pm(t)分别为岁差矩阵、章动矩阵、周日自转矩阵和极移矩阵。极移矩阵将历元平地球坐标系，即CTS转换到瞬时极地球坐标系。地球自转矩阵将瞬时极地球坐标系转换到真天球坐标系。章动矩阵将真平天球坐标系转换到瞬平天球坐标系。岁差矩阵将瞬平天球坐标系转换到CIS[J2000.0]。根据IERS规范(McCarthy，1996)，有[8]：

(12)

(13)

(14)

(15)

2.2.2 协议惯性坐标系到卫星轨道坐标系的转换

RTN的3个坐标轴在CIS下的单位向量[14]为：

(16)

即：

RCIS=MRRTN

(17)

式中，RRTN为质点RTN坐标系中的三维位置；M为从RTN到CIS的转换矩阵：

M=

(18)

2.2.3 卫星固定坐标系到协议惯性坐标系的转换

从SBF到CIS的转换可以通过太阳和卫星在协议惯性坐标系中的位置矢量来实现。更高精度转换一般需要用到卫星的姿态数据。姿态数据给出的是由星体固定坐标系到协议惯性系转换的四元数(或称欧拉参数)q1，q2，q3，q4。由四元数表示的坐标转换矩阵[14]为：

Q=

(19)

所以，将SBF下的向量RSBF转换到CIS下的向量RCIS的公式为：

RCIS=QRSBF。

(20)

3 算法流程设计

在利用低轨卫星力学模型和坐标系统转换建立低轨卫星轨迹信息(位置+速度)状态方程的基础上，兼顾轨迹准确度与工程计算量的要求，利用Runge-Kutta算法完成数值积分，得到各历元时刻的低轨卫星轨迹信息。

3.1 Runge-Kutta数值积分算法

Runge-Kutta算法[6]是一种工程上广泛应用的高精度单步算法。已知方程导数与初值信息：

(21)

式中，y0为矢量y的初值；f(t，y)为矢量y关于变量t的一阶偏导数。则根据Runge-Kutta算法：

y(tn+1)≈y(tn)+h·Φ=η(tn+1) ,

(22)

式中，η(tn+1)为y(tn+1)的近似解；Φ为增量函数，h=tn+1-tn。根据经典RK4 Runge-Kutta算法[6]，增量方程为：

(23)

式中，

(24)

根据Oliver Montenbruck和Eberhard Gill的研究，RK4解的局部截断误差[6]有以下特征：

eRK4=|y(tn+1)-η(tn+1))≤const·h5，

(25)

且RK4解与四阶泰勒展开多项式精度相当：

(26)

式中，y(tn)(i)，(i=2，3，4)为矢量y关于变量t的i阶偏导数。

3.2 低轨用户轨迹仿真算法设计

低轨用户轨迹状态量为其位置与速度[15]，即

(27)

式中，R为低轨用户三维位置矢量；V为其三维速度矢量。对于任意已知R，V，可低轨卫星受力模型与坐标转化方法得到其对应的加速度α，则轨迹状态量关于时间t的一阶偏导数为：

(28)

在此基础上，在给定初值R0，V0的前提下，利用RK4 Runge-Kutta算法即可计算得到各个历元时刻的轨迹信息。

低轨轨迹仿真算法流程如图1所示。

图1 低轨轨迹仿真算法流程Fig.1 Flow chart of LEO track simulation algorithm

4 算法仿真及验证

为了测试低轨用户轨迹仿真算法的逼真性与可靠性[16-17]，以STK10.0软件仿真的低轨卫星轨道数据作为理论参考值进行测试验证。

4.1 轨迹参数设定

轨迹仿真开始于UTC时间2010年1月1日0时0分0秒，持续86 400 s，仿真步长为10 s，输出地心地固坐标系的位置速度。低轨用户的六参数设置如表1所示。

表1 低轨用户轨迹初值设置Tab.1 Initial value setting of LEO user track

地球重力场选用70阶EGM96模型，不考虑固体潮影响；大气密度计算选用NRLMSISE2000模型；太阳光压选用圆柱地影模型；考虑太阳和月球的三体引力影响，选用JPL的DE405星历计算日月位置；考虑相对论效应。

4.2 仿真结果及分析

本次仿真的轨迹与STK10.0软件仿真轨迹比较结果显示，经过统计24 h后其三维位置误差保持在10.0 m之内，速度误差保持在0.012 m/s之内。

三维位置误差分布如图2所示，三维速度误差分布如图3所示。

图2 三维位置误差分布Fig.2 Error distribution graph of three-dimensional position

图3 三维速度误差分布Fig.3 Error distribution graph of three-dimensional velocity

通过上述测试结果分析，低轨用户轨迹仿真精度能到达较高水平，典型摄动力模型下逼近国外STK10.0软件轨道仿真精度水平。因此，该仿真算法在逼真性和可靠性上完全可以满足卫星导航模拟器数学仿真的需要。

5 结束语

在调研国内外已有轨迹仿真技术方案的基础上，结合卫星导航模拟器应用实际需求，深入研究了低轨用户轨迹仿真的数学模型，确定了低轨用户轨迹仿真的算法，完成了低轨用户轨迹仿真的流程设计与软件实现。经过仿真与分析，低轨用户轨迹仿真软件与STK10.0软件的仿真结果偏差小于10 m，达到了很高的精度，提升了卫星导航模拟器的性能，有力支持了低轨卫星地面验证工作。