CN111505679B - 一种基于星载gnss的leo初轨确定方法 - Google Patents

Abstract

本发明公开一种基于星载GNSS的LEO初轨确定方法，包括：步骤1)利用GNSS伪距观测值进行SPP解算，得到地心地固系下所有有效历元的接收机天线位置序列，所有有效历元组成的时段为begT1‑endT1；步骤2)计算惯性系下的接收机天线位置序列；步骤3)利用begT1和其下一个有效历元的惯性系下的接收机天线位置计算轨道根数，从而计算LEO卫星质心位置序列；步骤4)使用动力学模型进行轨道积分并对LEO卫星质心位置序列进行轨道拟合，若轨道拟合异常，则将begT1增加一个时间间隔，转入步骤3)；否则，进入步骤5)；步骤5)轨道拟合成功的时段为begT2‑endT1，获取该时段拟合后的惯性系下的LEO卫星质心位置和速度序列；步骤6)计算begT0‑begT2和endT0‑endT1内的惯性系下的LEO卫星质心位置和速度序列。

Description

一种基于星载GNSS的LEO初轨确定方法

技术领域

本发明属于LEO精密定轨领域，尤其涉及一种基于星载GNSS的LEO初轨确定方法。

背景技术

基于星载GNSS的LEO精密定轨，是指在LEO卫星上搭载测地型GNSS接收机，接收LEO卫星天顶上空的GNSS信号，得到伪距、载波相位和载噪比等观测值，利用上述观测值结合卫星动力学模型进行数据处理解算得到LEO卫星质心的精密轨道。在基于星载GNSS的LEO精密定轨解算中，通常是采用最小二乘等方法进行参数估计，最终得到待估参数(包括：LEO卫星轨道根数或对应的位置和速度、动力学参数以及GNSS载波相位模糊度、接收机钟差等)。由于GNSS观测方程和卫星动力学模型联立得到的LEO精密定轨解算方程是非线性的，使用小二乘等方法进行解算的时候首先需要把解算方程进行线性化。在解算方程线性化的时候，需要有待估参数的近似初始值。这些待估参数的初始近似值中LEO轨道根数(或对应的位置和速度)的近似初始值简称初轨。基于星载GNSS的LEO初轨确定所采用的方法，通常是使用GNSS伪距观测值进行标准单点定位(SPP)计算，然后使用起始几个历元的SPP定位结果转换成密切轨道根数(或转换成对应的位置和速度)。但是上述传统的初轨确定方式，如果起始几个历元的GNSS伪距观测值质量不佳或者观测值缺失(接收机异常或者某些GNSS***星座建设不完全等原因)，造成起始几个历元的SPP定位结果存在误差很大值或者序列不连续，则会导致转换成的密切轨道根数(或对应的位置和速度)错误甚至无解，进而引起后续整个时段的LEO精密定轨解算都无法正常进行，因此初轨的可靠性确定对LEO精密定轨的业务化稳定运行具有举足轻重的作用。

发明内容

鉴于传统的基于GNSS的LEO初轨确定方法可靠性方面的不足，本发明给出了一种基于星载GNSS的LEO初轨确定方法，该方法利用GNSS观测数据连续采样的特点并充分结合卫星动力学模型进行轨道积分及轨道拟合，轨道拟合过程中进行异常判断，并按历元时间顺序进行递进迭代处理，最终更可靠的确定原始数据时段的LEO初轨。用这种方法确定LEO初轨，具有不需依赖外部数据条件、算法简易、结果稳定可靠等特点。

为实现上述目的，本发明提供了一种基于星载GNSS的LEO初轨确定方法，所述方法包括：

步骤1)在原始数据时段begT0-endT0内，利用GNSS伪距观测值进行SPP解算，得到地心地固系下所有有效历元的接收机天线位置序列，所有有效历元组成的时段为begT1-endT1；所述有效历元为能够获得SPP解算结果的历元；

步骤2)在时段begT1-endT1内，基于地心地固系下所有有效历元的接收机天线位置序列，计算相应有效历元的惯性系下的接收机天线位置序列；

步骤3)利用begT1和其下一个有效历元的惯性系下的接收机天线位置计算轨道根数；从而计算近似LEO卫星质心位置和速度序列，以此作为卫星姿态矢量将惯性系下接收机天线位置序列转换成LEO卫星质心位置序列；

步骤4)使用动力学模型进行轨道积分并对LEO卫星质心位置序列进行轨道拟合，若轨道拟合异常，则将begT1增加一个时间间隔，若新的begT1为有效历元，转入步骤3)；否则，进入步骤5)；

步骤5)记轨道拟合成功的时段为begT2-endT1，begT2≥begT1；获取该时段拟合后的惯性系下的LEO卫星质心位置和速度序列；

步骤6)使用步骤5)得到的惯性系下的LEO卫星位置和速度序列的第一个向量和最后一个向量，利用动力学模型计算加速度并分别向前和向后进行两次积分，分别得到begT0-begT2和endT0-endT1内的惯性系下的LEO卫星质心位置和速度序列；

步骤7)按时间顺序合成时段begT0-endT0内连续等采样间隔的惯性系下LEO卫星质心位置和速度序列；并转换成地心地固系下LEO卫星质心位置和速度序列。

作为上述方法的一种改进，所述步骤1)具体包括：

在原始数据时段begT0-endT0内，基于伪距观测方程，计算地心地固坐标系下接收机天线位置；在计算过程中，当某个历元的有效观测卫星数小于4颗或者PDOP>5，则不予进行SPP解算，该历元为无效历元；

在时间区间begT0～endT0剔除无效历元，将第一个有效历元时刻记为begT1，begT1≥begT0，最后一个有效历元时刻记为endT1，endT1≤endT0，则begT1～endT1为SPP解算成功的时段，该时段内地心地固坐标系下的接收机天线位置序列为

k为有效历元的个数。

作为上述方法的一种改进，所述步骤3)具体包括：

在时段begT1-endT1内，惯性系下接收机天线位置序列为

利用begT1和其下一个有效历元的惯性系下的接收机天线位置，计算得begT1时刻处的密切轨道根数；

利用begT1时刻处的密切轨道根数，计算时段begT1-endT1内近似LEO卫星质心位置和速度序列，该序列与有效历元的时刻相对应；

根据近似LEO卫星质心位置和速度序列，将惯性系下接收机天线位置序列

转换成LEO卫星质心位置序列

作为上述方法的一种改进，所述卫星动力学模型采用JGM3地球重力场模型、日月引力摄动、光压摄动和大气阻力摄动；轨道积分采用四阶龙格库塔方法；在轨道积分时，计算每个历元积分的卫星位置与地心位置(0，0，0)之间的距离，一旦距离小于地球平均半径，则判断轨道积分计算失败，同时表明该次轨道拟合异常；

所述轨道拟合使用扩展卡尔曼滤波器进行，未知参数只估计位置和速度向量，扩展卡尔曼滤波器的观测方程的矩阵A为：

扩展卡尔曼滤波器从begT1开始进行前向滤波，到endT1结束后，然后从endT1开始往后进行后向平滑，最终实现begT1～endT1时段精度均匀的轨道拟合值；

若SPP解算结果和扩展卡尔曼滤波器一步预测值相差大于3δ，δ为SPP标称定位精度，则SPP解算结果作为粗差剔除；若从某个历元开始SPP解算结果和卡尔曼滤波器一步预测值相差一直大于3δ，则认为卡尔曼滤波器发散，表明该次轨道拟合异常。

作为上述方法的一种改进，所述步骤6)具体包括：

若begT2>begT0，使用

的begT2时刻的位置和速度向量

利用卫星动力学模型计算卫星加速度，从begT2时刻开始向前进行一次和二次轨道积分，得到begT0～begT2时段的性坐标系下的LEO卫星质心位置和速度序列

若endT0>endT1，使用

最后一个时刻endT1的坐标和速度向量

利用卫星动力学模型计算卫星加速度，从endT1时刻开始向后进行一次和二次轨道积分，得到endT1～endT0时段的惯性坐标系下的LEO卫星质心位置和速度序列

作为上述方法的一种改进，所述步骤7)具体包括：

把惯性系下的LEO卫星质心位置和速度序列

和

按时间顺序合并，最终得到begT0和endT0时间段内惯性系下的连续等间隔采样的LEO卫星质心位置和速度序列，记为

把惯性系下LEO卫星位置和速度向量序列

转换成地心地固系下连续等间隔采样的LEO卫星位置和速度向量序列

作为上述方法的一种改进，所述方法还包括：将地心地固系下连续等间隔采样的LEO卫星位置和速度向量序列

存入初轨标准交换文件，供后续精密定轨作为初轨输入文件使用。

本发明的优势在于：

1、本发明的方法无需额外的地基监测数据，只靠自身星载GNSS观测值，可实现自主获取稳定可靠的LEO卫星初轨；

2、本发明的方法不需要依赖LEO星载GNSS的多普勒观测值，仅使用GNSS单频或双频伪距观测值；

3、本发明的方法以历元为单位按时间先后顺序进行递进计算及迭代判断，过程简单、结果稳定可靠；

4、本发明的方法确定的初轨对应的位置精度达米级，完全满足后续LEO卫星精密定轨的需求。

附图说明

图1为本发明的基于星载GNSS的LEO初轨确定方法的流程图。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明。

鉴于传统基于星载GNSS的LEO初轨确定在可靠性方面的漏洞，本发明设计了一种遍历递进方法，并充分利用轨道积分及轨道拟合，最终得到更加可靠的初轨。本发明的技术方案为：首先进行伪距单点定位(SPP)解算，得到有效历元时段的接收机天线位置序列；然后根据开始两个相邻有效历元的SPP位置把接收机天线坐标序列转换成LEO质心位置序列(偏向改正)；把SPP得到的LEO质心位置序列作为观测值，使用动力学模型和轨道积分对LEO质心位置序列进行轨道拟合，得到优化后的有效时段的卫星位置/速度序列。如果轨道拟合失败，则有效开始时刻向前递进，然后从偏向改正开始处理。不断循环迭代，直到处理成功。最后进行坐标转换并写入初轨标准交换文件，供后续精密定轨(POD)使用。

如图1所示，本发明提供了一种基于星载GNSS的LEO初轨确定方法，在本实施例中，采用LEO星载GPS的单频定位观测数据，具体指L1频点上的C/A码伪距观测值，采样间隔为30秒。本申请还可以采用北斗定位观测数据或其它导航卫星的定位观测数据。

该方法的具体实施步骤如下：

步骤一、对待解算数据时段内所有观测历元全部进行SPP解算，每个历元时刻处的SPP解算方法如下：

基于GPS的单频伪距观测方程(忽略电离层影响)为：

式中，

为L1频点上的C/A码伪距的观测值，上标j表示卫星序号。

ρ^j为几何距离，上标j表示卫星序号。其具体表达式为：

式中，(X_s,Y_s,Z_s)为GPS卫星天线相位中心位置，通过广播星历文件计算得到，(x,y,z)为接收机天线相位中心位置，为待估参数。

dT为接收机钟差。

dt^j为GPS卫星钟差，上标j表示卫星序号。

为C/A码伪距的的测量误差，包括多路径和热噪声，上标j表示卫星序号。

将式(1)线性化后的误差方程矩阵形式可表示为：

y＝Aδ-W，P_l (2)

式中，

其中，

x＝[x y z dT]^T

对每一个历元的观测方程，由式(2)进行参数估计的步骤为：

(1)第1次计算

令X⁰中各参数的值均为0，由式(2)依据最小二乘原则可求得未知参数改正量δ为：

δ＝(A^TP_lA)^-1A^TP_lW (3)

由此可得到估计的未知参数为：

(2)第i次迭代(i>1)

用i表示迭代次数，将上一次的参数估计值作为当前迭代的初值，即

有第i次迭代的未知参数改正量为：

δ_i＝(A^TP_lA)^-1A^TP_lW_i (5)

由此可得到估计的未知参数为：

其中，判断迭代收敛的条件为：

式(6)中的

为GPS接收机天线相位中心处的地心地固坐标系(WGS84)下的位置。在计算过程中，有效观测卫星数小于4颗或者位置精度因子(PDOP)>5，则不予进行SPP解算，该历元SPP位置缺失。

令原始数据时段的开始时刻和结束时刻分别记为begT0和endT0。搜索时段begT0～endT0运行SPP成功的历元，第一个成功的历元时刻记为begT1(begT1≥begT0)，最后一个成功的历元时刻记为endT1(endT1≤endT0)，时段begT1～endT1也称为SPP解算有效历元序列时段。该时段内地心地固坐标系下的接收机天线位置序列为

k为有效历元的个数。

在时段begT1-endT1内，基于地心地固系下所有有效历元的接收机天线位置序列，计算相应有效历元的惯性系下的接收机天线位置序列

步骤二、获取第一个有效历元时刻begT1以及下一个相邻有效历元时刻对应的位置

和

通过目前标准转换公式计算得到begT1时刻处的密切轨道根数。利用begT1时刻处的密切轨道根数计算begT1到endT1时段内SPP有效历元时刻对应的J2000惯性系下LEO卫星质心位置和速度序列，以此作为卫星姿态矢量将惯性系下接收机天线位置序列转换成LEO卫星质心位置序列(该过程称为偏心改正)，记为

步骤三、把SPP的定位结果转换得到的位置序列

作为观测值，使用动力学模型进行轨道积分并对位置序列

进行轨道拟合，得到J2000惯性系下begT1到endT1时间段的坐标和速度序列

轨道动力学拟合使用扩展卡尔曼滤波器进行，未知参数只估计位置和速度向量，不对动力学参数进行估计。其中扩展卡尔曼滤波器的观测方程设计矩阵为：

卡尔曼滤波器的协方差更新使用状态转移矩阵通过误差传播率计算，其状态转移矩阵仅考虑二体问题进行计算；卡尔曼滤波器的状态更新采用上一个历元的卡尔曼滤波的状态和动力学模型计算的卫星加速度进行二次积分实现；位置的过程噪声分别设为：径向R方向0.01m/s，切向T方向0.02m/s，法向N方向0.16m/s；速度的过程噪声设为0。

卫星动力学模型考虑20阶JGM3地球重力场模型、日月引力摄动、光压摄动(光压系数取近似值，默认取1.0)和大气阻力摄动(阻尼系数取近似值，默认取2.0)。轨道积分采用四阶龙格库塔方法，积分步长选30秒。其他步骤中涉及的卫星动力学模型和轨道积分方法均与此相同。

在轨道积分的时候，每个历元积分的卫星位置都与地心位置(0，0，0)之间计算距离，一旦距离小于地球平均半径6371km，则轨道积分计算失败，表明该次轨道拟合计算失败。

卡尔曼滤波器从begT1开始进行前向滤波，到endT1结束后，然后从endT1开始往后进行后向平滑，最终实现begT1～endT1时段精度均匀的轨道拟合值。

在卡尔曼滤波计算过程中，若SPP结果和卡尔曼滤波器一步预测值相差大于3δ(δ为SPP标称定位精度，GPS取5米)，则SPP结果作为粗差剔除。若从某个历元开始SPP结果和卡尔曼滤波器一步预测值相差一直大于3δ，则认为卡尔曼滤波器发散，表明该次轨道拟合计算失败。

在轨道拟合执行过程中，若道拟合计算失败，则begT1的时间增加5分钟(即begT1＝begT1+5min，endT1保持不变，相应的持续时间变更为lastSec＝lastSec-5min)，如果新的begT1不是有效历元，则begT1＝begT1+5min，直至begT1是有效历元。然后重新从步骤二开始执行，直到轨道拟合执行成功，则跳出开始历元的遍历循环，继续向下执行步骤四和步骤五；否则，若轨道拟合一直失败，直到历元遍历执行到begT1>endT1-5min时，轨道拟合仍然没有执行成功，则表明使用星载GNSS数据进行初轨确定失败(表明所有历元观测值均质量不佳，这种情况即使有初轨，后续的基于星载GNSS的LEO精密定轨也无法进行)，不再向下执行步骤四和步骤五。

步骤四、若begT1>begT0，使用

的第一个时刻begT1的坐标和速度向量

利用卫星动力学模型计算卫星加速度，从begT1时刻开始向前进行一次和二次轨道积分，得到begT0～begT1时段的J2000惯性坐标系下的LEO卫星质心位置和速度序列

若endT0>endT1，使用

最后一个时刻endT1的坐标和速度向量

利用卫星动力学模型计算卫星加速度，从endT1时刻开始向后进行一次和二次轨道积分，得到endT1～endT0时段的J2000惯性坐标系下的LEO卫星质心位置和速度序列

把J2000惯性坐标系下的LEO卫星质心位置和速度序列

和

按时间顺序合并，最终得到begT0和endT0时间段内J2000惯性坐标系下的LEO卫星质心位置和速度序列，记为

步骤五、把步骤四得到的J2000惯性坐标系下卫星位置和速度向量序列

转换成WGS84坐标系下30秒等间隔采样的卫星位置和速度向量序列

并存入SP3格式的标准交换文件中。该SP3文件作为后续LEO卫星精密定轨解算的初轨输入文件使用。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (6)

1.一种基于星载GNSS的LEO初轨确定方法，所述方法包括：

步骤1)在原始数据时段begT0-endT0内，利用GNSS伪距观测值进行SPP解算，得到地心地固系下所有有效历元的接收机天线位置序列，所有有效历元组成的时段为begT1-endT1；所述有效历元为能够获得SPP解算结果的历元；

步骤2)在时段begT1-endT1内，基于地心地固系下所有有效历元的接收机天线位置序列，计算相应有效历元的惯性系下的接收机天线位置序列；

步骤3)利用begT1和其下一个有效历元的惯性系下的接收机天线位置计算轨道根数；从而计算近似LEO卫星质心位置和速度序列，以此作为卫星姿态矢量将惯性系下接收机天线位置序列转换成LEO卫星质心位置序列；

步骤4)使用动力学模型进行轨道积分并对LEO卫星质心位置序列进行轨道拟合，若轨道拟合异常，则将begT1增加一个时间间隔，若新的begT1为有效历元，转入步骤3)；否则，进入步骤5)；

步骤5)记轨道拟合成功的时段为begT2-endT1，begT2≥begT1；获取该时段拟合后的惯性系下的LEO卫星质心位置和速度序列；

步骤6)使用步骤5)得到的惯性系下的LEO卫星位置和速度序列的第一个向量和最后一个向量，利用动力学模型计算加速度并分别向前和向后进行两次积分，分别得到begT0-begT2和endT0-endT1内的惯性系下的LEO卫星质心位置和速度序列；

步骤7)按时间顺序合成时段begT0-endT0内连续等采样间隔的惯性系下LEO卫星质心位置和速度序列；并转换成地心地固系下LEO卫星质心位置和速度序列；

所述卫星动力学模型采用JGM3地球重力场模型、日月引力摄动、光压摄动和大气阻力摄动；轨道积分采用四阶龙格库塔方法；在轨道积分时，计算每个历元积分的卫星位置与地心位置(0，0，0)之间的距离，一旦距离小于地球平均半径，则判断轨道积分计算失败，同时表明该次轨道拟合异常；

所述轨道拟合使用扩展卡尔曼滤波器进行，未知参数只估计位置和速度向量，扩展卡尔曼滤波器的观测方程的矩阵A为：

扩展卡尔曼滤波器从begT1开始进行前向滤波，到endT1结束后，然后从endT1开始往后进行后向平滑，最终实现begT1～endT1时段精度均匀的轨道拟合值；

若SPP解算结果和扩展卡尔曼滤波器一步预测值相差大于3δ，δ为SPP标称定位精度，则SPP解算结果作为粗差剔除；若从某个历元开始SPP解算结果和卡尔曼滤波器一步预测值相差一直大于3δ，则认为卡尔曼滤波器发散，表明该次轨道拟合异常；

初轨为LEO轨道根数的近似初始值或者LEO轨道对应的位置和速度的近似初始值。

2.根据权利要求1所述的基于星载GNSS的LEO初轨确定方法，其特征在于，所述步骤1)具体包括：

在原始数据时段begT0-endT0内，基于伪距观测方程，计算地心地固坐标系下接收机天线位置；在计算过程中，当某个历元的有效观测卫星数小于4颗或者PDOP>5，则不予进行SPP解算，该历元为无效历元；

在时间区间begT0～endT0剔除无效历元，将第一个有效历元时刻记为begT1，begT1≥begT0，最后一个有效历元时刻记为endT1，endT1≤endT0，则begT1～endT1为SPP解算成功的时段，该时段内地心地固坐标系下的接收机天线位置序列为

k为有效历元的个数。

3.根据权利要求2所述的基于星载GNSS的LEO初轨确定方法，其特征在于，所述步骤3)具体包括：

在时段begT1-endT1内，惯性系下接收机天线位置序列为

利用begT1和其下一个有效历元的惯性系下的接收机天线位置，计算得begT1时刻处的密切轨道根数；

利用begT1时刻处的密切轨道根数，计算时段begT1-endT1内近似LEO卫星质心位置和速度序列，该序列与有效历元的时刻相对应；

根据近似LEO卫星质心位置和速度序列，将惯性系下接收机天线位置序列

转换成LEO卫星质心位置序列

4.根据权利要求1所述的基于星载GNSS的LEO初轨确定方法，其特征在于，所述步骤6)具体包括：

若begT2>begT0，使用

的begT2时刻的位置和速度向量

利用卫星动力学模型计算卫星加速度，从begT2时刻开始向前进行一次和二次轨道积分，得到begT0～begT2时段的性坐标系下的LEO卫星质心位置和速度序列

若endT0>endT1，使用

最后一个时刻endT1的坐标和速度向量

利用卫星动力学模型计算卫星加速度，从endT1时刻开始向后进行一次和二次轨道积分，得到endT1～endT0时段的惯性坐标系下的LEO卫星质心位置和速度序列

5.根据权利要求4所述的基于星载GNSS的LEO初轨确定方法，其特征在于，所述步骤7)具体包括：

把惯性系下的LEO卫星质心位置和速度序列

和

按时间顺序合并，最终得到begT0和endT0时间段内惯性系下的连续等间隔采样的LEO卫星质心位置和速度序列，记为

把惯性系下LEO卫星位置和速度向量序列

转换成地心地固系下连续等间隔采样的LEO卫星位置和速度向量序列

6.根据权利要求5所述的基于星载GNSS的LEO初轨确定方法，其特征在于，所述方法还包括：将地心地固系下连续等间隔采样的LEO卫星位置和速度向量序列

存入初轨标准交换文件，供后续精密定轨作为初轨输入文件使用。

Images