CN105203110B - 一种基于大气阻力模型补偿的低轨卫星轨道预报方法 - Google Patents

Abstract

本发明公开了一种基于大气阻力模型补偿的低轨卫星轨道预报方法，涉及航空航天领域。该方法：S1，利用星载GNSS接收机伪距观测量和IGS发布的超快星历对卫星进行定轨，得到卫星的动力学参数，所述动力学参数包括卫星位置、速度、C_D、C_R和逐段常量经验加速度参数；S2，对经验加速度建模，然后该模型与确定性动力学模型组成增强型动力学模型，最后利用定轨最终历元的卫星的位置和速度作为初值进行轨道预报。本发明所述方法只要求存在星载接收机伪距观测量即可；并且易于实现，对定轨后的经验加速度拟合，即可应用于轨道预报；同时，所述方法还具有一定普适性和可操作性，在一定程度上补偿了大气阻力模型误差，提高预报精度，应用于在轨高精度导航。

Description

一种基于大气阻力模型补偿的低轨卫星轨道预报方法

技术领域

本发明涉及航空航天领域，尤其涉及一种基于大气阻力模型补偿的低轨卫星轨道预报方法。

背景技术

目前低轨卫星的轨道预报精度不高，尤其是轨道迹向(切向)预报误差随时间增加迅速变差，主要原因在于很难建立高精度大气阻力模型。现有技术都集中在对大气阻力模型本身的改进上，通过调节大气阻力模型中一、两个系数或是通过改进大气密度模型期望提高预报精度。然而现有方法扭曲了大气模型中动力学参数的物理意义而且未明显提高预报精度，不具有普适性。

发明内容

本发明的目的在于提供一种基于大气阻力模型补偿的低轨卫星轨道预报方法，在已有GNSS观测量的时段对未模型化力进行估计，并经建模后用于轨道预报，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明所述基于大气阻力模型补偿的低轨卫星轨道预报方法，通过对定轨过程中的经验加速度建模，与预先定义的确定性动力学模型组成增强型动力学模型，在一定程度上补偿预报期间的大气阻力模型误差，具体包括以下步骤：

S1，利用星载接收机伪距观测量和IGS发布的超快星历对卫星进行定轨，得到卫星的动力学参数，所述动力学参数包括卫星位置、速度、C_D、C_R和经验加速度；

S2，对经验加速度建模，然后与确定性动力学模型组成增强型动力学模型，最后利用定轨最终历元的卫星的位置和速度作为初值进行轨道预报。

优选地，步骤S1，具体按照下述步骤实现：

S11，利用IGS发布的超快星历，对已有星载接收机观测数据时段的伪距值做单点定位，得到观测历元的单点位置值；

S12，采用简化动力学最小二乘批处理方法对单点位置值进行平滑，完成定轨，得到平滑后的卫星位置和速度，还得到C_D、C_R和经验加速度。

更优选地，步骤S12，采用简化动力学最小二乘批处理方法对单点位置值进行平滑，完成定轨，具体按照下述步骤实现：

A1，设定子区间长度为τ，将整个定轨时间区间[t₀,t_n]等间隔划分为n段；

A2，经验加速度a_i＝(a_iR,a_iT,a_iN)^T在子区间[t_i,t_i+1)上保持常量，其中t_i＝t₀+iτ，i＝0,1,…,n；经验加速度a_i的方向分别为卫星轨道径向e_R(t)、切向e_T(t)和法向e_N(t)，其中，轨道径向e_R(t)、切向e_T(t)和法向e_N(t)的计算公式为公式组(1)；

其中，r和v为卫星在地心惯性系中的位置和速度向量；

A3，在确定性动力学模型的基础上增加经验加速度项，得到卫星动力学模型(2)，

其中a为确定性动力学模型加速度矢量和；a_iR表示径向经验加速度；a_iT表示切向经验加速度；a_iN表示法向经验加速度；ξ_i(t)定义为公式(3)：

A4，利用单点位置值作为伪观测量，待估参数Y包括卫星初值位置和初值速度y₀、C_D、C_R以及经验加速度a_i，i＝0,1,…,n-1，见公式(4)；

其中，表示卫星初始状态；

A5，最小二乘法方程(5)为计算待估参数初值Y₀的修正量ΔY；

H^TH·ΔY＝H^T(z-h(Y₀)), (5)

其中，h表示轨道单点位置值模型，h(Y₀)表示加入初值的模型值，为设计矩阵；z为单点位置值由单点定位过程得到，假设共有n_T个历元数据，任意历元单点位置值z_i＝(x_i,y_i,z_i)中包括卫星位置三分量，i＝0,1,…,n_T-1；

A6，轨道位置模型值h的计算利用初值状态y₀和动力学参数作为输入，根据卫星动力学模型(2)进行数值积分得到；

最小二乘批处理过程是迭代过程，每次迭代完成后得到待估参数初值Y₀的修正量ΔY，并根据修正量ΔY对待估参数进行更新，把更新后的参数值Y₀+ΔY作为初始值再进行迭代计算，当待估参数收敛时，完成定轨。

更优选地，在待估参数收敛前，经过3～4次迭代计算。

更优选地，所述确定性动力学模型包括：

静态重力场：GGM01C 120×120；

潮汐力：固体潮IERS 2003、海潮UT/CSR3.0和极潮IERS 2003；

三体引力：日月引力解析表达式；

大气阻力：NRLMSISE 2000、NOAA solar flux daily、geomagnetic activity3hourly；

太阳辐射光压：双锥体地球阴影模型；

地球反照光压：短波可见光辐射和长波红外辐射。

优选地，步骤S2具体按照下述步骤实现：

S21，采用高阶傅立叶级数对切向经验加速度时间序列{a_iT|i＝0,1,…,n}进行拟合，拟合模型(6)为：

其中，a₀为常偏系数，a_j和b_j为频率系数，ω为基频，t₀≤t≤t_n；

S22，所述拟合模型与确定性动力学模型组成增强型动力学模型，并利用定轨得到的卫星初值状态y₀作为输入，通过轨道外推完成轨道预报。

本发明的有益效果是：

本发明所述方法只要求存在星载接收机伪距观测量即可；并且易于实现，定轨后的经验加速度经拟合，即可应用于轨道预报；同时，所述方法还具有一定普适性和可操作性，在一定程度上补偿了大气阻力模型误差，提高预报精度，应用于在轨高精度导航。

附图说明

图1是基于大气阻力模型补偿的低轨卫星轨道预报流程图；

图2是利用2006年1月3日的GRACE-A星载伪距观测量对切向经验加速度进行估计的结果(离散点)以及傅立叶级数拟合曲线；

图3是在预报时间起点为2006年1月7日零时，预报3天，考虑和不考虑切向经验加速度傅立叶级数拟合模型的轨道切向预报精度对比；三角形组成的曲线表示不考虑拟合模型的轨道切向预报精度，正方形组成的曲线表示考虑拟合模型的轨道切向预报精度；

图4是在预报时间起点为2006年1月6日零时，预报3天，考虑和不考虑切向经验加速度傅立叶级数拟合模型的轨道切向预报精度对比；三角形组成的曲线表示不考虑拟合模型的轨道切向预报精度，正方形组成的曲线表示考虑拟合模型的轨道切向预报精度；

图5是在预报时间起点为2006年1月5日零时，预报3天，考虑和不考虑切向经验加速度傅立叶级数拟合模型的轨道切向预报精度对比；三角形组成的曲线表示不考虑拟合模型的轨道切向预报精度，正方形组成的曲线表示考虑拟合模型的轨道切向预报精度；

图6是在预报时间起点为2006年1月4日零时，预报3天，考虑和不考虑切向经验加速度傅立叶级数拟合模型的轨道切向预报精度对比；三角形组成的曲线表示不考虑拟合模型的轨道切向预报精度，正方形组成的曲线表示考虑拟合模型的轨道切向预报精度；

图7是在3天的预报时间内，图3至图6中两种轨道预报策略的最大切向预报偏差对比示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

在基于GNSS星载接收机观测量(伪距)的简化动力学定轨过程中，预先定义的逐段常量经验加速度使卫星的确定性动力学模型能够匹配GNSS观测量的精度(伪距测量精度为分米级)的同时，而且保持卫星轨道的鲁棒性。经验加速度作为确定性动力学模型的一种补偿方式，其大小事先并不知道，而是随其他动力学参数在定轨过程中进行估计。定轨后的经验加速度参数具有伪周期、余弦曲线特点，幅值较为稳定，若将其特征模型化并应用于接下来几天的预报中，将在一定程度上补偿确定性动力学模型精度，提高低轨卫星预报精度。所以本发明的原理来源于精密定轨过程中用于补偿确定性动力学模型误差的经验加速度可模型化。因预报误差主要集中在切向上(其与径向和法向有数量级上的差别)，原因是大气阻力模型误差对于低轨卫星来说较大。需要说明的是本发明并不对大气阻力模型本身进行改进，而是利用切向经验加速度模型对大气阻力模型进行补偿。关于本申请更详细的解释如下：

本发明所述方法通过对定轨过程中的经验加速度建模，与预先定义的确定性动力学模型组成增强型动力学模型，在一定程度上补偿预报期间的大气阻力模型误差，具体包括以下步骤：

S1，利用星载接收机伪距观测量和IGS发布的超快星历对卫星进行定轨，得到卫星的动力学参数，所述动力学参数包括卫星位置、速度、C_D、C_R和经验加速度；

S2，对经验加速度建模，然后与确定性动力学模型组成增强型动力学模型，最后利用定轨最终历元的卫星的位置和速度作为初值进行轨道预报。

(一)首先是定轨过程

只要存在卫星星载接收机的观测数据就可以完成定轨，具体按照下述步骤实现：

S11，利用IGS发布的超快星历，对已有星载接收机观测数据时段的伪距值做单点定位，得到观测历元的单点位置值；在该过程中，星载接收机伪距观测量的测量精度已满足需求，同时采用IGS发布的超快星历是合适的；

S12，采用简化动力学最小二乘批处理方法对单点位置值进行平滑(单点定位过程包括在该批处理方法中)，完成定轨，得到平滑后的卫星位置和速度，还得到C_D、C_R和经验加速度。该步骤目的为得到较高精度的卫星预报初始时刻的位置和速度，以及相关的动力学参数。本申请中所述最小二乘批处理即同时处理所有观测数据，实现观测值与模型值在最小二乘意义下的最优拟合，该方法是相对于序贯处理方法而言的。

因单点定位原理与方法是本领域公知的尝试已比较成熟，本申请中不做介绍了，下面详细介绍平滑过程：

A0，在简化动力学最小二乘批处理方法(Batch Least-Squares,Batch LSQ)中引入经验加速度模型；

A1，设定子区间长度为τ，将整个定轨时间区间[t₀,t_n]等间隔划分为n段；

A2，经验加速度a_i＝(a_iR,a_iT,a_iN)^T在子区间[t_i,t_i+1)上保持常量，其中t_i＝t₀+iτ，i＝0,1,…,n；经验加速度a_i的方向分别为卫星轨道径向e_R(t)、切向e_T(t)和法向e_N(t)，其中，轨道径向e_R(t)、切向e_T(t)和法向e_N(t)的计算公式为公式组(1)；

其中，r和v为卫星在地心惯性系中的位置和速度向量；经验加速度a_i的取值预先未知，需要在定轨过程中作为参数估计；

A3，在确定性动力学模型的基础上增加经验加速度项，得到卫星动力学模型(2)，

其中a为确定性动力学模型加速度矢量和；a_iR表示径向经验加速度；a_iT表示切向经验加速度；a_iN表示法向经验加速度；ξ_i(t)定义为公式(3)：

A4，利用单点位置值作为伪观测量，待估参数Y包括卫星初值位置和初值速度y₀、C_D、C_R以及经验加速度a_i，i＝0,1,…,n-1，见公式(4)；

其中，表示卫星初始状态；

A5，最小二乘法方程(5)为计算待估参数初值Y₀的修正量ΔY；

H^TH·ΔY＝H^T(z-h(Y₀)), (5)

其中，h表示轨道单点位置值模型，h(Y₀)表示加入初值的模型值，为设计矩阵；z为单点位置值由单点定位过程得到，假设共有n_T个历元数据，任意历元单点位置值z_i＝(x_i,y_i,z_i)中包括卫星位置三分量，i＝0,1,…,n_T-1；

A6，轨道位置模型值h的计算利用初值状态y₀和动力学参数作为输入，根据卫星动力学模型(2)进行数值积分得到；

最小二乘批处理过程是迭代过程，每次迭代完成后得到待估参数初值Y₀的修正量ΔY，并根据修正量ΔY对待估参数进行更新，把更新后的参数值Y₀+ΔY作为初始值再进行迭代计算，当待估参数收敛时，完成定轨。

步骤A3中，所述确定性力学模型为表1；

表1：确定性力学模型

(二)预报过程

因定轨后的经验加速度呈伪周期、余弦曲线特点，而高阶傅立叶级数能将经验加速度参数的周期性特征分解为具有不同周期的简单函数的和，包括正弦和余弦函数。前面讲到，本申请主要关心用以补偿大气阻力模型误差的切向经验加速度，所以只需对切向经验加速度拟合。采用高阶傅立叶级数对切向经验加速度时间序列{a_iT|i＝0,1,…,n}进行拟合，拟合模型(6)为：

其中，a₀为常偏系数，a_j和b_j为频率系数，ω为基频，t₀≤t≤t_n，以上系数都为待拟合参数；

S22，所述拟合模型与确定性动力学模型组成增强型动力学模型，并利用定轨得到的卫星初值状态y₀作为输入，通过轨道外推完成轨道预报。

在本申请中经验加速度参数在定轨和预报这两个过程中都发挥了至关重要的作用，定轨中的经验加速度是为了使动力学模型能够匹配星载接收机观测量的精度，而拟合后经验加速度也可应用在接下来的轨道预报中补偿动力学模型误差。

实施例

利用GRACE-A卫星实际在轨观测数据对本方法的轨道预报效果做测试。GRACE-A卫星平均轨道高度460km，为低轨近圆轨道卫星，配备双频GPS接收机。数据来源为：星载接收机观测数据和参考轨道来自于美国宇航局JPL实验室物理海洋学数据分发存档中心(PO.DAAC)，选取4天数据(2006年1月3日～6日)。为满足轨道预报的需求，GPS卫星星历采用IGS发布的超快星历(观测部分)，其延迟3-9小时，位置精度优于0.1m。按照图1，实现定轨和预报，并利用JPL发布的参考轨道(精度为厘米级)对预报精度进行评估。

定轨过程：GRACE-A星载GPS伪距观测数据经单点定位计算后，再应用简化动力学最小二乘批处理定轨。定轨参数设置如下：设置GRACE-A卫星质量487.2kg，大气阻力面积0.97m²，辐射光压面积3.41m²；设置观测历元间隔30s；设置经验加速度相关时间τ为600s。确定性动力学模型如表1所示，动力学轨道积分采用变步长RKF78阶方法，可以满足高精度数值积分要求。对时长为一天的观测数据为进行批处理，共完成4组测试结果。单点定位的位置精度为2m(1σ，下同)，最小二乘批处理后的切向位置精度为0.2m，切向速度精度优于1.0×10^-4m/s。C_D和C_R的平均值分别为3.3和1.5。以2006年1月3日为例，给出切向经验加速度的估计值，如图2离散点所示。经验加速度为逐段常量，呈现准周期、余弦曲线的特点，均值为零，幅值不超过30nm/s²。拟合过程以实现周期匹配为准，不用过度关注幅值的拟合，因为未模型化力的幅值具有一定随机性。其他日期的切向经验加速度结果具有类似结果。

预报过程：利用高阶傅立叶级数对切向经验加速度拟合，式(6)为拟合式，表示为时间的函数：

仍以2006年1月3日为例给出拟合结果，如图2中连续曲线所示。其中阶数n为8，a₀为-0.693，ω为0.5855，a_i和b_i如表2所示。傅立叶级数拟合曲线只需在周期上实现较好的拟合即可，无需在幅值上准确拟合，因未模型化的大气阻力误差具有一定随机性。其他日期的切向经验加速度采用相同方法拟合。

表2切向经验加速度的傅立叶级数拟合参数(2006年1月3日)

i	1	2	3	4	5	6	7	8
									ai	-1.706	-2.438	0.5277	-2.315	0.3038	1.268	4.307	0.7206
bi	1.011	2.339	-0.7873	-3.724	-0.8137	1.851	-3.704	-0.4407

每组定轨过程最后历元的位置和速度值作为预报的初值。动力学模型由确定性动力学模型(见表1)和切向经验加速度级数拟合模型(6)组成，见公式(7)

其中a为确定性动力学模型加速度矢量和，注意：每组数据对应的拟合模型的参数是不尽相同的。数值积分方法仍采用与定轨相同的RKF78阶方法。本申请做了4组轨道预报，2006年1月3日至6日每天观测数据完成定轨后立刻预报3天，如利用1月3日伪距观测量定轨，轨道预报时间从1月4日零时至7日零时。

对比两种轨道预报策略——考虑切向经验加速度级数和不考虑切向经验加速度级数拟合模型，预报精度如图3所示，用于验证的GRACE-A参考星历来自JPL。从图中看出，切向误差随预报时间增加而迅速增大，未考虑级数拟合模型策略预报3天的最大误差从50到150米不等，而考虑级数拟合模型策略的最大误差都在60米以内。级数拟合模型或延缓预报误差的偏移速度(如图3、图4和图5)，或逆转预报误差的偏移方向(如图6)。图7为在3天的预报时间内，两种轨道预报策略的最大切向偏差对比，考虑级数拟合模型策略将预报精度平均提高2.3倍，表明级数拟合模型补偿了一定大气阻力模型误差。由于径向和法向的预报误差比较小，所以并未考虑径向和法向经验加速度的补偿。经验加速度与当时的空间环境有密切关系，所以只要存在星载接收机观测数据就进行定轨，不断更新切向经验加速度的拟合模型，更新轨道预报结果。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：本发明所述预报新方法只要求存在星载接收机伪距观测量即可；并且易于实现，对定轨的经验加速度拟合后，即可应用于轨道预报；同时本发明方法也可应用于在轨高精度导航。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (5)

1.一种基于大气阻力模型补偿的低轨卫星轨道预报方法，其特征在于，通过对定轨过程中的经验加速度建模，与预先定义的确定性动力学模型组成增强型动力学模型，补偿预报期间的大气阻力模型误差，具体包括以下步骤：

S1，利用星载接收机伪距观测量和IGS发布的超快星历对卫星进行定轨，得到卫星的动力学参数，所述动力学参数包括卫星位置、速度、C_D、C_R和经验加速度；

S2，对经验加速度建模，然后与确定性动力学模型组成增强型动力学模型，最后利用定轨最终历元的卫星的位置和速度作为初值进行轨道预报；

其中，步骤S2具体按照下述步骤实现：

S21，采用高阶傅立叶级数对切向经验加速度时间序列{a_iT|i＝0,1,…,n}进行拟合，拟合模型为：

其中，a₀为常偏系数，a_j和b_j为频率系数，ω为基频，t₀≤t≤t_n；n表示定轨时间区间[t₀,t_n]等间隔划分的段数；j＝1,2,...,m表示傅里叶级数的阶数；

S22，所述拟合模型与确定性动力学模型组成增强型动力学模型，并利用定轨得到的卫星初值状态y₀作为输入，通过轨道外推完成轨道预报。

2.根据权利要求1所述方法，其特征在于，步骤S1，具体按照下述步骤实现：

S11，利用IGS发布的超快星历，对已有星载接收机观测数据时段的伪距值做单点定位，得到观测历元的单点位置值；

S12，采用简化动力学最小二乘批处理方法对单点位置值进行平滑，完成定轨，得到平滑后的卫星位置和速度，还得到C_D、C_R和经验加速度。

3.根据权利要求2所述方法，其特征在于，步骤S12，采用简化动力学最小二乘批处理方法对单点位置值进行平滑，完成定轨，具体按照下述步骤实现：

A1，设定子区间长度为τ，将整个定轨时间区间[t₀,t_n]等间隔划分为n段；

A2，经验加速度a_i＝(a_iR,a_iT,a_iN)^T在子区间[t_i,t_i+1)上保持常量，其中t_i＝t₀+iτ，i＝0,1,…,n；经验加速度a_i的方向分别为卫星轨道径向e_R(t)、切向e_T(t)和法向e_N(t)，其中，轨道径向e_R(t)、切向e_T(t)和法向e_N(t)的计算公式为公式组(1)；

其中，r和v为卫星在地心惯性系中的位置和速度向量；

A3，在确定性的动力学模型的基础上增加经验加速度项，得到卫星动力学模型(2)，

其中，a为确定性动力学模型加速度矢量和；a_iR表示径向经验加速度；a_iT表示切向经验加速度；a_iN表示法向经验加速度；ξ_i(t)定义为公式(3)：

A4，利用单点位置值作为伪观测量，待估参数Y包括卫星初值位置和初值速度y₀、C_D、C_R以及经验加速度a_i，i＝0,1,…,n-1，见公式(4)；

其中，表示卫星初始状态；

A5，最小二乘法方程(5)为计算待估参数初值Y₀的修正量△Y；

H^TH·△Y＝H^T(z-h(Y₀)), (5)

其中，h表示轨道单点位置值模型，h(Y₀)表示加入初值的模型值，为设计矩阵；z为单点位置值由单点定位过程得到，假设共有n_T个历元数据，任意历元单点位置值z_i＝(x_i,y_i,z_i)中包括卫星位置三分量，i＝0,1,…,n_T-1；

A6，轨道位置模型值h的计算利用初值状态y₀和动力学参数作为输入，根据卫星动力学模型(2)进行数值积分得到；

最小二乘批处理过程是迭代过程，每次迭代完成后得到待估参数初值Y₀的修正量△Y，并根据修正量△Y对待估参数进行更新，把更新后的参数值Y₀+△Y作为初始值再进行迭代计算，当待估参数收敛时，完成定轨。

4.根据权利要求3所述方法，其特征在于，在待估参数收敛前，经过3～4次迭代计算。

5.根据权利要求3所述方法，其特征在于，所述确定性动力学模型包括：

静态重力场：GGM01C 120×120；

潮汐力：固体潮IERS 2003、海潮UT/CSR3.0和极潮IERS 2003；

三体引力：日月引力解析表达式；

大气阻力：NRLMSISE 2000、NOAA solar flux daily、geomagnetic activity 3hourly；

太阳辐射光压：双锥体地球阴影模型；

地球反照光压：短波可见光辐射和长波红外辐射。