一种基于GNSS广播星历的高轨卫星高精度自主定轨方法与流程

本发明属于卫星导航技术领域，涉及一种基于gnss广播星历的高轨卫星定轨方法。

背景技术：

基于星载gnss接收机的定位与授时技术具有自主性强、定位精度高和体积功耗小等优点，已经被广泛应用于低轨卫星的自主导航与精密定轨任务。但是，对于高轨任务而言，由于卫星轨道高度高于gnss卫星，因此只能接收来自于地球另外一边的gnss信号。由于主瓣信号大部分被地球遮掩而需要接收旁瓣信号，这种远距离的信号传播方式会引入很大的空间衰减，使得高轨场景下gnss可用性急剧下降。极弱信号的捕获和跟踪能够大幅提升gnss的可用性，但由于在高轨场景下的gnss观测量噪声大，gnss卫星的几何分布差，导致导航定位误差较大。

现有的基于gnss的高轨卫星自主定轨方法，在定轨过程中并未对gnss广播星历引起的缓变系统误差(gnss卫星轨道误差和钟差)进行建模估计，而将其作为gnss测量误差的一部分，并容许其传递至定轨定位结果中去。在高轨场景下，由于gnss卫星的几何分布很差，几何因子pdop的数值大都在十几至几十量级，其与由gnss广播星历引起的用户测距误差ure相乘后，可造成数十米级的定轨定位误差。

文献《基于高轨航天器的gnss接收机技术》(《中国空间科学技术》，vol.37,no.3,jun.2017)针对高轨道航天器gnss信号微弱的特点，采用长时间积分处理的梳状滤波方法、差分相干累加比特同步算法和基于动力学模型补偿的扩展卡尔曼滤波自主定轨算法设计gnss接收机。文献《gnss接收机在探月三期中的应用研究》(《飞行器测控学报》，vol.34,no.2,apr.,2015)针对低信噪比环境下的高灵敏度接收机系统应用进行设计，采用弱信号捕获算法和弱信号跟踪算法实现-175dbw的灵敏度，最终采用基于轨道动力学模型的卡尔曼滤波方法实现了gnss自主定轨算法。

上述两篇文献均介绍了针对高轨和月球场景的高灵敏度gnss接收机，分析了gnss可用性，并为了提高自主导航的精度，引入了轨道动力学模型和扩展卡尔曼滤波进行了自主定轨。然而，对于由于gnss广播星历引入的系统误差，两篇文献皆未提及如何通过在轨处理予以消除或削弱，在其各自的仿真场景设计中，也未明确对此部分真实存在的误差进行模拟。

技术实现要素：

本发明解决的技术问题是：克服现有技术的不足，提供了一种基于gnss广播星历的高轨卫星高精度自主定轨方法，可以实现在没有实时gnss精密星历支持情况下，仅使用gnss广播星历和gnss观测量完成高轨卫星高精度自主定轨。

本发明的主要原理是：通过星载增强扩展卡尔曼滤波器(aekf)完成gnss观测量和轨道动力学模型的紧耦合，获取轨道动力学模型递推结果短期精度高和gnss观测量长期不发散的优点；同时，针对广播星历引入的gnss卫星轨道误差和时钟误差缓慢变化的特性，在滤波过程对其进行联合估计，并从gnss原始观测量中予以扣除，以削弱其对定轨结果的影响。

本发明的技术解决方案是：一种基于gnss广播星历的高轨卫星高精度自主定轨方法，包括如下步骤：

(1)获取gnss卫星跟踪列表prns，确定与增强状态xaug相对应的协方差矩阵paug；其中xaug＝[δ1,δ2,δ3…δn]n，δ1,δ2,δ3…δn为卫星跟踪列表prns中正在跟踪的n颗gnss卫星由于广播星历误差引入的系统偏差，paug为n×n维对角矩阵，初始化为对角阵，对角线元素为与n颗gnss卫星相对应的初始方差值；

(2)根据gnss观测量的噪声水平生成观测噪声矩阵rk；

(3)利用轨道动力学模型，将上一历元滤波估计的卫星位置和速度通过轨道递推，计算得到当前历元的预测值

(4)将接收机状态xrecv和增强状态xaug一起组成系统状态x＝[xrecv,xaug]n+8，采用增强扩展卡尔曼滤波器对所述系统状态进行滤波，得到更新后的卫星位置和速度；其中，rx,ry,rz为gnss接收机的位置矢量，b为gnss接收机的钟差，vx,vy,vz为gnss接收机的速度矢量，为gnss接收机的钟速。

所述的卫星跟踪列表prns中，若有gnss卫星跟踪状态发生变化导致prns发生变化，则对增强状态xaug及其对应的协方差矩阵paug进行调整，具体为：

若有一颗新的gnss卫星出现，则将增强状态xaug增加一维，并将对应于新gnss卫星的广播星历误差初始化为0，将增强状态的协方差矩阵paug从n×n维增加至(n+1)×(n+1)维，并将新加入的行和列初始化；

若有一颗gnss卫星停止跟踪，则将增强状态xaug减少一维，删除停止跟踪的gnss卫星的广播星历误差，同时将增强状态的协方差矩阵paug从n×n维减少至(n-1)×(n-1)维，删除停止跟踪的gnss卫星的广播星历误差对应的行和列；

若跟踪的gnss卫星保持不变，但有某颗gnss卫星的广播星历发生了更新，则将增强状态xaug中对应的元素初始化为0，将增强状态的协方差矩阵paug中对应的行和列初始化。

所述的rk矩阵为对角阵，对角元素由gnss伪距噪声σρ,k和伪速噪声组成，具体表示如下：

所述的轨道递推方法为龙格-库塔4阶积分法。

所述的增强扩展卡尔曼滤波器中，状态量的递推通过轨道力学模型积分得到，系统矩阵通过在状态量处线性化轨道力学模型得到，gnss原始观测作为观测量，观测矩阵为由gnss卫星的位置和估计的接收机位置计算得到，系统误差矩阵qk由轨道动力学模型的精度决定，观测误差矩阵rk由gnss观测的噪声水平决定。

所述的增强扩展卡尔曼滤波器中，系统转移矩阵为

其中，表示k时刻的状态预测值，表示k-1时刻的状态估计值，f为系统函数，表示在t时刻状态对时间的导数，τs为状态递推的时间间隔，

ax,ay,az表示gnss接收机的加速度，由轨道动力学模型计算得到。

所述的增强扩展卡尔曼滤波器中，观测矩阵hk＝[hrecvhaug]2n×(n+8)

表示第一颗gnss卫星和接收机间的视线矢量，角标为n的代表第n颗gnss卫星和接收机间的视线矢量。

所述的增强扩展卡尔曼滤波器中，协方差矩阵p表示为：

precv为接收机状态xrecv的协方差矩阵。

本发明与现有技术相比的优点在于：本发明方法可以在没有实时gnss精密星历支持情况下，仅使用gnss广播星历和gnss观测量完成高轨卫星高精度自主定轨。方法中对gnss广播星历的缓变误差进行了估计，在自主定轨过程中从gnss观测量中予以扣除，削弱了gnss广播星历误差对自主定轨的影响。本发明通过使用增强扩展卡尔曼滤波器(aekf)融合gnss观测量和轨道动力学模型进行星上实时滤波处理，在geo轨道可达到优于10m(3drms)的自主定轨精度。

附图说明

图1为本发明方法的原理框图；

图2为本发明实施例中的gnss广播星历轨道误差示意图；

图3为本发明实施例中的gnss广播星历星钟误差示意图；

图4为本发明实施例中现有自主定轨方法与本发明方法结果比对示意图。

具体实施方式

如图1所示，为本发明方法的基本原理框图，本发明方法主要包括以下四个主要的环节：

1、gnss系统偏差处理；

2、自适应观测噪声矩阵生成；

3、轨道动力学模型辅助；

4、增强扩展卡尔曼滤波器aekf滤波。

如图1所示，本发明方法的输入分别为k时刻可用gnss卫星列表prns，gnss卫星在惯性系下的位置和速度rsat,k,vsat,k，gnss伪距测量值和伪速测量值ρk,伪距误差估计值和伪速的误差估计值σρ,k,输出为高轨卫星在k时刻的位置估计值和速度估计值

一、gnss系统偏差编辑处理

这里的gnss系统偏差编辑，主要是根据gnss卫星跟踪列表prns和广播星历状态对增强状态增强状态的协方差矩阵的维数和数值进行调整。这里xaug＝[δ1,δ2,δ3…δn]n，δ1,δ2,δ3…δn为卫星跟踪列表prns中正在跟踪的n颗gnss卫星由于广播星历误差引入的系统偏差，paug为xaug的协方差矩阵，k-1表示k时刻的前一时刻，其为aekf滤波器的协方差矩阵p中与增强状态xaug相关的分块，表示如下：

对于aekf滤波器将在后续步骤中进行详细的说明，但是由于滤波器中涉及该参数初始化，此处先行说明。precv为接收机状态xrecv的协方差矩阵，rx,ry,rz为gnss接收机的位置矢量，b为gnss接收机的钟差，vx,vy,vz为gnss接收机的速度矢量，为gnss接收机的钟速。

具体的处理流程包括：

(1)定轨计算开始前，对进行初始化操作。

根据卫星跟踪列表prns中的gnss卫星数量n生成n×n维对角矩阵，对角元素置初始方差值，该初始方差值由gnss广播星历的误差量级和卡尔曼滤波最优调谐状态共同决定。本发明中，初始方差值取(10m)²。

(2)当前历元若有一颗新的gnss卫星出现，则将增强状态从n维增加至n+1维，并将对应于新gnss卫星的广播星历误差δn+1初始化为0，将增强状态的协方差矩阵从n×n维增加至(n+1)×(n+1)维，并将新加入的行和列初始化(除对角线以外的元素均置零，对角线元素置初始方差值)。若新增多颗gnss卫星，每一颗新增卫星的处理方法均和前述一致；

(3)当前历元相比上个历元若有gnss卫星停止跟踪，则将增强状态从n维减少至n-1维，删除停止跟踪的gnss卫星的广播星历误差；同时将增强状态的协方差矩阵从n×n维减少至(n-1)×(n-1)维，删除停止跟踪的gnss卫星的广播星历误差对应的行和列；

(4)当前历元若有某颗gnss卫星的广播星历发生了更新，则将增强状态中对应的元素初始化为0；将增强状态的协方差矩阵中对应的行和列初始化(除对角线以外的元素均置零，对角线元素置初始方差值)。

二、自适应观测噪声矩阵生成

这里的自适应观测噪声矩阵生成，是指根据gnss观测量的噪声水平生成aekf的观测噪声矩阵rk，使aekf可以根据gnss观测量的质量自适应的调节滤波参数，达到最优估计。

这里，rk矩阵为对角阵，对角元素由gnss伪距噪声σρ,k和伪速噪声组成，具体表示如下：

三、轨道动力学模型辅助

这里的轨道动力学模型，用于计算卫星的轨道位置和速度，它对作用于空间飞行器的各种轨道摄动因素进行建模，对轨道摄动力引起的加速度进行数值积分，计算飞行器的运动状态。

本发明针对高轨卫星的轨道力学特性，考虑的轨道摄动因素主要包括地球低阶引力场(6×6阶)，月球和太阳的三体引力，太阳辐射压力等，并使用龙格-库塔4阶积分法对轨道运动状态进行递推计算。将上一历元滤波估计的高轨卫星位置和速度通过轨道递推计算当前历元的预测值

轨道摄动因素的建模方法和龙格-库塔4阶级分方法的具体实现参考《卫星轨道—模型,方法和应用》。

使用轨道动力学模型是为了给增强扩展卡尔曼滤波器的系统矩阵计算提供必要的信息，同时轨道动力学递推的结果作为增强扩展卡尔曼滤波器系统状态历元间递推的结果。

四、增强扩展卡尔曼滤波器aekf滤波

这里的增强扩展卡尔曼滤波器aekf，对gnss观测量和轨道动力学模型进行距离域耦合(紧耦合)，aekf的待估系统状态x由接收机状态xrecv和增强状态xaug一起组成：

x＝[xrecv,xaug]n+8

xaug＝[δ1,δ1,δ1...δn]n

其中，rx,ry,rz为gnss接收机的位置矢量，b为gnss接收机的钟差，vx,vy,vz为gnss接收机的速度矢量，为gnss接收机的钟速，aekf状态量的递推通过轨道力学模型积分得到，系统矩阵通过在aekf状态量处线性化轨道力学模型得到。gnss原始观测量作为aekf的观测，观测矩阵为由gnss卫星的位置和估计的接收机位置计算得到；系统误差矩阵qk由轨道动力学模型的精度决定，观测误差矩阵rk由gnss观测的噪声水平决定。

增强扩展卡尔曼滤波器aekf具体的处理流程包括：

(1)状态预测：

通过上一时刻状态估计值和系统函数递推估计当前时刻的状态。其中表示k时刻的状态预测值，表示k-1时刻的状态估计值，f为系统函数，表示在t时刻状态对时间的导数，即：

ax,ay,az表示gnss接收机的加速度，由轨道动力学模型计算得到。

对于增强状态xaug，认为其在当前历元的预测值等于上个历元估计值：

(2)系统矩阵获取：

接收机状态的系统矩阵frecv,k-1可通过在上一时刻的接收机状态估计值处对系统函数f进行线性化得到。

(3)转移矩阵求解：

通过级数展开对转移矩阵做近似求解。其中τs为状态递推的时间间隔，级数可取n＝2。将frecv,k-1带入此公式得到接收机状态的转移矩阵φrecv,k-1；同时，增强状态的转移矩阵φaug,k-1为单位阵in×n，则总的系统转移矩阵为：

(4)系统协方差矩阵预测：pk^-＝φk-1pk-1⁺φ^tk-1+qk

其中，pk-1⁺表示在k-1时刻系统状态误差矩阵估计值，pk^-表示在k时刻系统协方差矩阵预测值，在自主定轨开始时刻，为对角阵，对角元素为相对应的系统状态的初始方差估计值；qk表示系统噪声矩阵，由轨道动力学模型的误差特性和卡尔曼滤波最优调谐状态共同决定。

(5)卡尔曼增益计算：

其中，hk表示观测矩阵，由接收机预测位置、gnss卫星位置及增强状态与观测量间的映射关系共同计算得到：hk＝[hrecvhaug]2n×(n+8)

表示第一颗gnss卫星和接收机间的视线矢量，依此类推，角标为n的代表第n颗gnss卫星和接收机间的视线矢量。

(6)状态更新：其中，zk表示gnss原始观测量，即gnss伪距测量值ρ和伪速测量值构成：系统协方差矩阵更新pk⁺＝(i-kkhk)pk^-。

实施例

对本发明方法进行测试验证，使用双频gps接收机和spirentgss9000信号模拟器进行硬件测试。

spirentgss9000信号模拟器产生真实的gps信号，接收场景模拟地球同步轨道；gps星座参数使用2018年年积日125天真实的gps星座状态，包含31颗卫星；信号落地功率与gpsicdis-gps-200j保持一致。

在双频gps接收机获取的原始伪距中加入了由于广播星历引入的误差。广播星历误差由igs提供的最终精密星历产品作为真值和广播星历求差进行模拟，如图2和图3所示。图2表示在2018年年积日125天时，以igs提供的最终精密轨道作为真值，计算得到的各gps卫星广播星历轨道单轴误差；图3表示在2018年年积日125天时，以igs提供的最终精密星钟作为真值，计算得到的各gps卫星广播星历星钟误差。

同时使用单点最小二乘定位方法，其他文献中描述的基于ekf的轨道滤波自主定轨方法以及本发明的自主定轨方法(aeke)，结果如图4所示，使用最小二乘方法得到的自主定轨误差最大，具有较强的随机误差(误差曲线较‘胖’，快速的随机波动)及系统偏差(误差曲线在不同时段具有不同系统偏置，且系统偏置在缓慢变化)。使用经典的ekf轨道滤波方法得到的自主定轨误差得到了改善，随机误差大幅减小(误差曲线较‘光滑’)，但是仍然存在系统偏差(误差曲线的均值发生大幅变化)。使用本发明的aekf轨道滤波方法得到的自主定轨结果最优，不但随机误差较小，且很大程度削弱了系统偏差(误差曲线一直较为平稳，没有出现明显的系统性变化)。综上所述，本发明的自主定轨方法比使其他方法在精度上得到了大幅度提升。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。