一种基于星载GNSS接收机的高轨卫星自主导航滤波算法①

张 蓬，杨克元，王延光，蒙艳松

(中国空间技术研究院西安分院，西安 710000)

0 引言

全球导航系统(Global Navigation Satellite System,GNSS)设计之初主要针对地面用户和低空用户提供服务。由于GNSS在卫星导航与定轨领域显示出了明显的优势：极大的减轻了地面测控压力，提升了自主性、实时性与导航定轨精度，针对高轨卫星应用的特种GNSS接收机研究逐渐成为了热点。2016年，国际GNSS组织工作组B(ICG WG-B)在ICG工作会上达成了GNSS提供空间服务(Space Volume Service,SSV)的共识,明确规定了在GEO轨道的系统指标[1]。美国哥达德航天中心研制的高轨GNSS接收机Navigator，于2015年在MMS(Magnetospheric Multiscale)任务中成功服役，其上搭载的GEONS导航滤波器在Phase 1阶段中12倍地球半径的远地点可达到12米的导航精度[2]。美国的地球同步轨道(GEO)气象卫星GOCE-R上搭载了另一款高轨接收机Viceroy，其在仿真条件下在GEO轨道的导航精度达到了18米[3]。国内对高轨GNSS接收机的研究起步较晚，文献[4]、[5]分别针对GEO轨道和60000公里以下的地月转移轨道进行了基于GNSS的自主导航算法设计，达到了优于50米和100米的位置精度。

在高轨场景下基于GNSS的自主导航，由GNSS广播星历引入的用户测距误差(User Range Error，URE)会引起比低轨场景更加严重的位置精度损失。广播星历引入的GNSS卫星轨道误差和时钟误差不是随机分布的白噪声序列，在同一广播星历龄期内一般发生缓慢连续的变化。经典的扩展卡尔曼滤波器(Extended Kalman Filter，EKF)只有在观测噪声满足白噪声的假设条件时才能达到最优估计6]。因此，基于EKF的导航滤波器难以削弱广播星历引入的误差，该误差会直接由几何分布因子放大后反映在最终的定位结果中，造成数十米的位置误差。文章提出的自主导航滤波算法，使用增强扩展卡尔曼滤波器(Augmented Extended Kalman Filter，AEKF)对GNSS观测和轨道动力学模型进行紧耦合，同步估计高轨卫星的运动状态和广播星历的缓变误差，以期从GNSS原始观测中削弱和消除GNSS广播星历不准确的影响。文章使用STK 11.0软件和Spirent GNSS GSS9000信号模拟器共同搭造算法仿真平台，对GEO卫星星载GNSS接收机的标准最小二乘方法、基于经典EKF的导航滤波方法以及本文提出的基于AEKF的导航滤波方法的性能进行了比对，验证了本文算法的有效性。

1 仿真平台建立

使用STK11.0软件中的高精度轨道递推器(HPOP)生成高轨卫星的高保真动力学轨道，将此卫星轨道数据作为接收机运动轨迹导入Spirent信号模拟器。在Spirent模拟器中设置GNSS星座参数和接收机天线参数后进行共24小时场景仿真，生成GNSS卫星、接收机的位置/速度以及信号接收功率信息。之后根据GNSS卫星和高轨卫星的相对位置和速度精准计算真实距离及距离变化率，在其上叠加模拟的伪距误差和伪速误差即可得到伪距和伪速模拟量。最后，可使用生成的GNSS模拟观测量在Matlab中进行导航滤波算法仿真验证。

1.1 高轨卫星运动参数

高轨卫星设置为地球同步轨道(GEO)卫星，初始轨道6根数设置如表 1 所示，该GEO卫星的轨道示意图如图1所示。

1.2 GNSS系统参数

在本文的研究中，使用单GPS星座，从CelesTrak官网(www.celestrak.com)获取2018年年积日125天GPS星座的历书，并导入Spirent信号模拟器进行GPS星座模拟，共包含31颗有效GPS卫星。使用的信号为民用信号GPS L1 C/A。该信号的落地功率按照ICD文件[7]设置为-128.5dBm。对于高轨应用场景， GNSS卫星发射天线的旁瓣辐射性能至关重要。GPS IIM、IIR系列卫星的天线方向图使用洛克希德马丁公司的公开数据[8]，对于IIF系列卫星，由于缺乏公开数据，可使用其他系列卫星方向图数据的平均值进行代替[9]，如图2所示。

图2 GPS IIF系列卫星发射天线L1频段增益方向图[9]Fig.2 Transmission Antenna pattern at L1 frequency band of GPS IIF satellites[9]

1.3 接收机参数

使用类似文献[10]中设计的高轨高灵敏度GPS接收机。捕获灵敏度门限设定为-149dBm，跟踪灵敏度设定为-151dBm，具体参数见表2。

表2 高轨GNSS接收机参数Tab.2 HEO GNSS receiver parameters

高增益指向型接收天线会带来天线尺寸和重量的增加，会给高轨卫星平台提出更苛刻的要求。故在GPS可用性足够的条件下使用普通增益的接收天线，其幅度方向图如图3所示。

图3 高轨GPS接收机L1频段接收天线增益方向图Fig.3 Reception Antenna pattern at L1 frequency band of HEO GPS receiver

1.4 GNSS观测模拟

伪距模拟通过计算GNSS卫星和接收机间的真实距离ρ以及各种误差模拟量完成，公式如式(1)所示：

P=ρ-cδts+cδtr+δρbrdc+δρion+δρtrop+δρmul+εP

(1)

其中，cδts代表卫星时钟偏差，按照2018年年积日125天的历书中星钟参数进行模拟；cδts代表接收机时间偏差，使用一阶斜升模型机型模拟；δρbrdc代表由广播星历误差引入的距离误差，使用2018年年积日125天的GPS广播星历和国际GNSS服务组织IGS提供的GPS最终精密星历，以最终精密星历作为真值，获取GPS系统段轨道误差和时钟误差，并将其投影至GNSS卫星和接收机连线进行模拟。需要说明的是，广播轨道是以GNSS发射天线相位中心为参考点，而精密轨道则以卫星质心为参考点，需要在精密轨道上进行GNSS发射天线相位中心偏差(Phase Center Offsite，PCO)修正才可和广播轨道进行比对。由广播星历引入的距离误差见图 4，可明显看到其成缓慢变化的系统偏差特征；δρion表示由电离层引入的测距误差，δρtrop表示由对流层引入的测距误差，根据文献[11]的分析，高轨场景下直接丢弃穿过电离层的观测数据比使用电离层消除组合可以带来更加优异的定位性能，因此，本文的导航滤波方法采取只使用没有经过电离层的观测数据，δρion和δρtrop可建模为0；δρmul表示多径误差，建立为0.2m的随机误差[12]；εP表示热噪声误差，由信号接收功率和接收机的环路参数决定，利用文献[12]中介绍的码环跟踪热噪声经验公式进行计算。还有一些可以被模型精确改正的误差(例如相对论误差、接收天线PCO误差和硬件延迟误差等)未进行模拟。

图4 GPS 2018年年积日125天广播星历引起的伪距误差Fig.4 User range error caused by GPS broadcast ephemeris on 2018 doy 125

伪速模拟通过计算GNSS卫星和接收机间的真实距离变化率以及热噪声误差模拟量完成，如式(2)所示：

(2)

其中，εD表示热噪声误差，由信号接收功率和接收机的环路参数决定，利用文献[12]中介绍的载波环跟踪热噪声经验公式进行计算。

2 GEO轨道GNSS可用性

在建立的Spirent仿真场景下收集生成的仿真数据：接收GNSS信号功率、GNSS卫星位置/速度、接收机位置/速度等。首先根据接收机的捕获时间、捕获灵敏度、跟踪灵敏度等参数筛选可见GNSS卫星；然后根据GNSS卫星和接收机的位置进行过电离层检测并予以剔除，得到最终可用GNSS星表，如图5所示。不考虑电离层影响的可用星数均值约为8.40颗，删除过电离层卫星后可用星数约8.11颗，如文献[11]所分析，可用性并不会明显改变。

图5 GEO轨道 GPS可见星数示意图Fig.5 Visible GPS PRN numbers in GEO orbit

使用剔除电离层后的GNSS卫星和接收机位置计算GNSS几何分布因子，计算方法见文献[12]。结果如图6所示。可以看到，GEO场景下，高灵敏度GPS接收机的GDOP均值约为19.1(剔除过电离层的GNSS数据后)，在几何构型较差时可出现大于100的峰值。

图6 GEO轨道GPS卫星几何分布因子示意图Fig.6 GDOP values of GPS in GEO orbit

3 导航滤波算法设计

表3 经典EKF计算流程Tab.3 Calculation procedures of classical EKF

经典EKF的推导过程建立在系统噪声Q和观测噪声R为高斯白噪声的前提下。当Q或R呈现某种系统性的偏差特性时时，EKF无法达到最优估计。AEKF在系统状态向量中加入增强状态xaug，表示缓变误差分量，与原系统状态组成新的系统状态x′=[xxaug]，在滤波过程中使xaug吸收非白噪声误差以优化估计结果。

增强扩展卡尔曼滤波器AEKF；

轨道动力学模型；

GNSS系统偏差估计量编辑；

图7 AEKF导航滤波算法框图Fig.7 Diagram of AEKF based navigation filtering algorithm

3.1 AEKF滤波器

AEKF对GNSS观测量和轨道动力学模型进行距离域耦合，AEKF的系统状态x由接收机状态xrecv和增强状态xaug组成：

x=[xrecv,xaug]N+8

xaug=[δ1,δ2,δ3…δN]N

(3)

(4)

其中，f(xrecv,t)接收机状态的系统函数，表示在t时刻状态xrecv对时间的导数，含义如下：

(5)

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

转移矩阵求解：

Φaug,k-1=IN×N

其中，Φrecv,k-1为接收机状态的转移矩阵，通过Frecv,k-1近似计算(Frecv,k-1通过f(xrecv,t)线性化展开求得)，τs为状态推的时间间隔。增强状态的转移矩阵Φaug,k-1为单位阵I(N×N)。

观测矩阵计算：

Hk=[HrecvHaug]2N×(N+8)

(6)

Haug=I

02N×N

其中，[-uia1,x-uia1,y-uia1,z]表示第一颗GNSS卫星和接收机间的视线矢量。

3.2 轨道力学模型

针对高轨卫星的轨道力学特性，为了简化在轨运算负荷，只考虑主要的轨道摄动因素。轨道力学模型主要包括地球低阶引力场(6×6阶)，月球和太阳的三体引力，太阳辐射压力等。

3.3 系统状态及其协方差矩阵编辑

根据GNSS卫星跟踪列表和广播星历状态对增强状态xaug、增强状态的协方差矩阵Paug的维数和数值进行调整：

当前历元若有新的GNSS卫星出现，则将增强状态xaug从N维增加至N+1维，并将对应于新GNSS卫星的广播星历偏差初始化为0；将增强状态的协方差矩阵Paug从N×N维增加至(N+1)×(N+1)维，并将新加入的行和列初始化(除对角线置零，对角线元素置初始方差值)；

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

当前历元若有某颗GNSS卫星的广播星历发生了更新，则将增强状态xaug中对应的元素初始化为0；将增强状态的协方差矩阵Paug对应的行和列初始化。

4 仿真验证

分别采用最小二乘方法，经典的EKF导航滤波方法以及本文提出的AEKF导航滤波方法对GEO卫星进行自主导航仿真，结果如图 8及表 4所示。可以看到，即使使用高灵敏度接收机，最小二乘方法的3D RMS位置误差也达到了31.517米，使用经典的EKF导航滤波法可以将位置精度提升至13.676米，若使用本文提出的AEKF方法，则可进一步提升至8.874米，相对于经典EKF方法，提升幅度约33.3%，相对于最小二乘方法，提升了70.9%；使用最小二乘方法时，在几何分布因子较差时可达到726.124米的3D最大位置误差，但若使用滤波方法，则可借助于动力学轨道递推精度高的特点，使导航定位结果不出现大的跳变，因此，EKF方法和AEKF方法的3D最大位置误差仅分别为45.978米和32.707米，AEKF方法相对于EKF方法，性能提升约28.8%，相对于最小二乘方法，提升了约95.5%；表5给出了3种方法的1D位置误差，观察可发现，径向误差相对于切向和法向误差大了一个量级，这是因为在高轨场景下，GNSS卫星集中分布在天顶区域，造成径向距离和接收机钟差高度相关，难以解耦，导致径向误差变大。AEKF方法和EKF方法的性能差别也主要体现在径向误差，AEKF方法相对于EKF方法有约35.1%的性能提升，相对于最小二乘方法有约71.8%的性能提升。

图8 GEO轨道自主导航定位性能示意图Fig.8 Autonomous navigation result in GEO orbit

表4 GEO轨道自主导航定位3D位置精度Tab. 4 3D position accuracy of autonomous navigation in GEO orbit

表5 GEO轨道自主导航定位1D位置精度Tab.5 1D position accuracy of autonomous navigation in GEO orbit

5 结论

本文提出了一种基于GNSS接收机的高轨卫星自主导航滤波算法，针对高轨场景下GNSS几何分布因子差、用户等效测距误差大且包含GNSS广播星历引入的缓变误差等特点，利用AEKF在非白噪声环境下更加突出的性能，融合星载GNSS数据与轨道动力学模型，估计并削弱了GNSS广播星历误差影响，提升了高轨卫星自主导航定位精度。并通过STK11.0和Spirent GSS9000模拟器搭造仿真平台进行了算法的有效性验证。仿真结果显示，AEKF方法在GEO轨道进行自主导航可以达到约9.32米的3D RMS 位置精度；相对于经典EKF滤波方法，3D RMS位置误差提升了约33.3%，3D最大位置误差提升了约28.8%，相对于最小二乘方法，3D RMS位置误差提升了约70.9%，3D最大位置误差提升了约95.5%。