低轨导航增强卫星广播星历拟合模型研究

2023-09-09杨开伟刘涵宇孙腾达孙秀聪徐明

航天器工程 2023年4期

杨开伟刘涵宇孙腾达孙秀聪徐明

(1 中国电子科技集团公司第五十四研究所,石家庄 050081)(2 北京航空航天大学宇航学院,北京 102206)

近年来,全球卫星导航系统(GNSS)迅速发展,在人类社会生活的方方面面发挥着不可或缺的作用,涉及到测绘、交通、救援等诸多方面。随着自动驾驶、物联网、5G等新兴产业的发展,社会生产生活对卫星导航服务提出了更高的要求,对时空信息的需求发展为精准、实时、动态和全球[1-2]。GNSS的不足之处也逐渐显露出来。GNSS在常规情况下只能提供10m左右的定位精度,而无法满足用户高精度的需求;此外,导航卫星信号到达地面时极其微弱,不但导致其极易受到干扰和欺骗,而且导致其难以穿透物理遮蔽,无法在室内、水下等复杂地形中提供高精度的定位服务[1,3-4]。

导航增强技术是解决上述问题、提升导航服务性能的一种可行方案,能利用低轨导航增强卫星产生测距信号并与中高轨GNSS系统协同提供导航服务[5]。传统的导航增强系统采用的是地球静止轨道(GEO)卫星,但GEO卫星轨道资源有限,而且用户接收到的信号弱,容易受到干扰,无法满足用户日益增长的高精度定位需求[6]。低轨卫星则没有上述问题,其轨道高度低,信号的自由空间损耗低,落地信号更强,抗干扰性能更好。此外,由于低轨导航增强卫星的飞行速度较快,相等时间的位置变化较大,从而使得历元间量测信息的相关性下降,能提高系统的可观测度,加快载波相位模糊度的收敛速度[4]。例如,美国的铱下一代(Iridium NEXT)星座提供的卫星授时与定位服务,已经能够脱离GPS单独提供定位服务,因而低轨导航增强卫星可以在GPS受到干扰的情况下提供备份的导航服务[4]。我国的珞珈一号卫星搭载有导航增强载荷,并开展了一系列低轨卫星增强相关试验[7]。

高精度确定用户的位置需要导航卫星高精度的位置信息,导航卫星的位置信息由广播星历计算,因此广播星历的高精度计算方法得到了广泛的研究。目前,应用广泛的广播星历参数模型有2种,主要是针对中高轨卫星,一种是以GPS为代表的基于轨道根数的广播星历参数模型;一种是以格洛纳斯(GLONASS)为代表的基于轨道状态的广播星历参数模型[8]。GPS广播星历参数模型的精度和GLONASS相当,但前者的预报能力更强,用户算法更简单;后者的优势是其参数更少,占用导航信号的资源少。卫星设备、通信链路和用户性能等多方面的限制,要求广播星历模型具有参数少、外推能力强、用户算法简单等特点,因此GPS广播星历参数模型应用更加广泛[9]。

传统的GPS系统采用的是包括6个开普勒轨道参数、3个轨道参数变化项、6个调和项系数和1个参考历元的16参数模型[10]。在原来16参数的基础上,GPS系统后续又发展了18参数模型。18参数模型中考虑卫星所受摄动力对卫星轨道半长轴、平均角速度及升交点赤经的影响,相较于16参数模型更加注重卫星轨道的瞬时变化,拟合精度更高[11]。16参数模型和18参数模型都是针对中高轨卫星设计的,而低轨卫星与中高轨卫星所受摄动力有所区别,如果将这种参数模型直接应用在低轨卫星上,将无法满足高精度的定位需求[12],因此需要根据低轨卫星的特点重新设计广播星历参数模型。文献[12]中利用轨道高度约1000km的低轨卫星在16参数模型的基础上设计了26参数的广播星历拟合算法,增加了轨道偏心率的变化率、升交角线性变化率和8个调和改正项振幅共10个参数,其2h轨道弧段的均方差小于10m,局部弧段拟合误差小于25m,但参数过多,因此占用的导航信号资源多。文献[6]中考虑了低轨卫星摄动力的短期变化快,在GLONASS广播星历参数模型的基础上设计了基于轨道状态的22参数模型。其研究显示:当低轨卫星轨道高度大于700km时,用户距离误差(FURE)精度在5cm以内;当轨道高度降低至500km时,FURE精度优于0.1m。虽然相对于26参数模型来说22参数模型参数少,但是用于确定卫星位置的用户算法外推效率较低,因此限制了其被广泛应用。针对26参数和22参数占用导航信号资源多、用户算法复杂的问题,本文提出了一种兼具拟合精度、计算效率和占用导航信号资源少的模型。考虑到大气阻力主要改变平面内的轨道要素,在18参数模型的基础上仅增加偏心率变化率这一参数,设计了19参数模型,并针对500～2000km高度的低轨卫星进行了19参数模型最小二乘拟合,检验了这种星历模型的有效性,验证了其拟合精度要优于18参数模型。本文提出的19参数模型对大倾角(不小于30°)卫星的长时弧段(4h)的位置拟合精度与26参数模型近似,但所需参数量大大减少,可为低轨导航增强卫星的广播星历提供一种兼顾精度与计算资源的模型。

1 19参数模型

1.1 广播星历参数

18参数模型采用轨道根数加摄动参数的形式表示,包括6个开普勒轨道参数、5个轨道根数变化项、6个调和项系数和1个参考历元。19参数模型在18参数模型的基础上添加偏心率变化率这一参数。表1给出了参考历元toe时刻19参数广播星历的符号说明。

表1 19参数广播星历参数定义Table 1 19-parameter broadcast ephemeris parameter definition

1.2 19参数广播星历用户算法

用户算法是指用户使用广播星历计算卫星位置的公式。在短时间内,使用广播星历可以较精确地描述卫星的位置。在之后的一段时间里,广播星历仍能描述卫星的位置,只是误差会随着外推时间变长而增大。对于一般轨道,19参数广播星历用户算法如下。

1)计算观测历元t到参考历元t0的时间差tk

tk=t-t0

(1)

式中:k为第k时刻。

2)计算观测历元半长轴A

(2)

3)计算卫星平均角速度

(3)

式中:μ为地心引力常数;Δn为观测历元相比参考历元的卫星平均角速度变化量;n为观测历元卫星平均角速度。

4)计算观测瞬间的平近点角Mk

Mk=M0+ntk

(4)

5)计算偏心率e

(5)

6)迭代计算偏近点角Ek

Ek=Mk+esinEk

(6)

7)计算真近点角Vk

(7)

8)计算摄动改正前的纬度辐角Φk

Φk=Vk+ω

(8)

9)计算摄动改正后的纬度辐角uk、卫星矢量半径rk和轨道倾角ik的摄动改正项δuk,δrk,δik

(9)

10)计算经过摄动改正的uk,rk,ik

(10)

11)计算卫星在轨道平面上的位置(xk,yk)

(11)

12)计算观测时刻的升交点经度(格林威治子午圈到卫星轨道升交点)

(12)

13)计算卫星在地心固连坐标系中的位置(X,Y,Z)

(13)

当直接利用常规的轨道根数方法对轨道倾角近似为0°的卫星进行广播星历参数拟合时,会导致法方程系数矩阵成为病态矩阵,使方程组的解误差变大,进而造成拟合精度降低。因此,为了解决小轨道倾角造成的卫星广播星历拟合精度降低甚至拟合失败的问题,本文采用轨道旋转法[11],先将地心惯性坐标系绕X轴(由地心指向春分点)旋转一个角度,在新坐标系计算广播星历,在用户预测卫星位置时将惯性坐标系旋转相反的角度,得到位置矢量在原惯性坐标系下的表达,进一步得到地心固连坐标系下的位置矢量。对于轨道倾角近似为0°的情况,坐标转换的过程见图1。

2 19参数广播星历拟合算法

19参数广播星历拟合算法流程如图2所示。首先,输入被拟合卫星一段时间的精密位置和初始时刻的精密速度,设置拟合时长,并对小倾角情况进行坐标系旋转。然后,设计19参数的迭代初值,除轨道六要素对应的参数外其他参数置为零。最后,利用最小二乘拟合算法迭代计算19参数,直至迭代改正量足够小时停止。其中,当倾角不大于ε时,认为属于小倾角情况,需要进行坐标系旋转,本文取ε为1°。通过判断改正量的模是否足够小(可以取阈值为1×10-2)或已经达到最大迭代次数(可以取10000次)作为迭代终止条件。下面对拟合算法进行详细介绍。

图2 广播星历拟合算法流程Fig.2 Broadcast ephemeris fitting algorithm flow

2.1 最小二乘拟合算法

广播星历参数的设计是实现高精度拟合的前提条件,而拟合算法的设计则是拟合精度和稳定性的重要保证。拟合算法的核心内容是参数估计,经典的参数估计法有最小二乘估计、极大似然估计、极大验后估计、贝叶斯估计等。其中,经典的最小二乘估计法是比较简单且容易理解的估计方法,在参数模型合理、观测值误差较小的情况下,采用最小二乘拟合算法就可以保证较高的拟合精度。

批处理最小二乘拟合算法的基本原理为是结合动力学模型与观测模型对状态量进行迭代平差,使得给定时间段内测量估计值与实际测量值之差的平方和最小。

对于一个非线性系统,有

z=h(x)+ν

(14)

式中:z为低轨卫星的精确位置信息;x为19参数广播星历;h(x)为通过用户算法计算出的低轨卫星的位置信息;ν为精确位置和计算位置的误差。

最小二乘拟合算法是使残差的平方和最小,即让评价函数J最小。

J=(z-h(x))T(z-h(x))

(15)

由于最小二乘拟合算法只能用于线性系统,因此需要对系统线性化,在x0处线性化后的残差可以表示为

Δz-HΔx0

(16)

可以使用线性最小二乘拟合算法求解式(14),得到

Δx0=(HTH)-1(HTΔz)

(17)

进而,原问题的最小二乘解为x=x0+Δx0。在线性化过程中舍弃了系统二阶及以上的信息,所以不可避免地引入误差,造成精度损失,即所得到的估计值并不是原问题的最优估计值。因此,还需要以所得估计值作为新的初值,重新迭代计算,直到算法收敛,即可得到状态量的最优估计值。

(18)

式中:rk为第k时刻的卫星位置;∂rk/∂e由文献[11]给出。

2.2 数值稳定修正

数值稳定修正是为了确保拟合算法能得到正确的结果,下面对修正方法进行具体介绍。

(2)升交点赤经、近地点幅角、平近点角。为保证数值稳定性,需要在每次迭代施加改正量之后对升交点赤经、近地点幅角、平近点角进行归一化,也就是让其落在[0,2π],之后再继续进行下一次迭代。例如,当Mk过大时,可能会因为数值问题无法通过迭代求出Ek。

(3)解线性方程组。最小二乘拟合算法问题需要求解线性方程组(HTH)Δx0=HTΔz,但HTH条件数很小,甚至接近病态,直接求逆或者直接使用高斯消元法求解误差较大。因此,本文使用下三角上三角(LU)分解法求解线性方程组。LU分解法可以将1个矩阵分解为1个下三角矩阵和1个上三角矩阵的乘积。如果使用杜利特尔(Doolittle)分解,得到的下三角矩阵为单位下三角矩阵;如果使用克劳特(Crout)分解,得到的上三角矩阵为单位上三角矩阵。首先,对(HTH)Δx0=HTΔz中的HTH进行LU分解,即LU=HTH,则(HTH)Δx0=HTΔz则可以写为LUΔx0=HTΔz,令y=UΔx0,可以先解线性方程组Ly=HTΔz,得到y,再解线性方程组y=UΔx0,得到Δx0,从而相比于直接解线性方程组能够降低系数矩阵的病态程度,误差更小。

3 仿真结果及分析

本文使用高精度轨道预报器预报的卫星位置作为精密星历,高精度轨道预报器中非球形引力摄动考虑70阶70次全球超高阶地球重力场模型(EGM2008),大气阻力模型使用质谱仪非相关散射模型(NRLMSISE-00),同时考虑日月三体引力;参考时刻为2010-09-19T04:00(UTC),参考时刻的偏心率取0,升交点赤经、近地点幅角和平近点角均取0°。

下面对不同轨道高度(500km,1000km,1500km,2000km)及不同轨道倾角(0°,30°,60°,90°)的轨道进行19参数广播星历最小二乘拟合。设置拟合初值如下。参考时刻半长轴取为参考时刻卫星的半长轴,参考时刻升交点赤经变化率取为0,19参数中的参考时刻、偏心率、升交点赤经、近地点幅角和平近点角按照卫星参考时刻的实际数据设置,其余参数设为0。卫星位置采样时间间隔为1min,即每分钟获取1组位置,拟合时长分别考虑20min,30min,60min,2h,3h,4h的情况。

绘制上述仿真场景下的位置均方根误差(RMSE),计算公式为

(19)

式中:Δx,Δy,Δz为地心固连坐标系下的位置误差。

表2～5中给出了18参数和19参数模型拟合的RMSE曲线对各个时刻的均方根,图3～6中给出了同一轨道倾角、同一拟合时长时不同轨道高度对拟合效果的影响。

图3 轨道倾角0°、拟合时长60min时不同轨道高度对18参数和19参数模型的拟合效果Fig.3 Fitting results of different orbit altitudes on 18-parameter and 19-parameter models when inclination angle is 0° and fitting time is 60 minutes

表2 轨道倾角为0°时不同轨道高度拟合结果的RMSETable 2 RMSEs of fitting results for different orbit altitudes at an inclination angle of 0° m

从表2和图3可以看出:对于任何拟合时长,拟合RMSE随轨道高度的升高而减小,原因是轨道高度越高,同等时间下轨道的非线性越弱,最小二乘法更能发挥优势;同时,随着轨道高度的升高及非球形引力摄动和大气阻力摄动的下降,轨道更接近二体轨道,因此拟合误差会有所下降。19参数模型拟合时长为20min时精度达到厘米级,30min时精度达到分米级,60min时精度达到米级。随着拟合时间的增长,RMSE有所增大,这是因为轨道根数预报是基于二体模型的,即使广播星历增加了摄动修正项,但其只有6个系数,远远不足以描述全部的摄动力,因此也只能解决短时间内的拟合问题,随着时间增长,误差仍然会越来越大。总体上来说,19参数模型的拟合精度要高于18参数模型。

从表3和图4可以看出:轨道倾角为30°时同样可以得到与轨道倾角为0°时相同的结论,即对于任何拟合时长,拟合RMSE随轨道高度升高而减小,整体上19参数模型的精度要高于18参数模型。对比轨道倾角为30°和0°的2种场景,可以发现:当轨道高度和拟合时长相等时,30°轨道的拟合误差略大于0°轨道的拟合误差。

图4 轨道倾角30°、拟合时长60min时不同轨道高度对18参数和19参数模型的拟合效果Fig.4 Fitting results of different orbit altitudes on 18-parameter and 19-parameter models when inclination angle is 30° and fitting time is 60 minutes

表3 轨道倾角为30°时不同轨道高度拟合结果的RMSETable 3 RMSEs of fitting results for different orbit altitudes at an inclination angle of 30° m

从表4和图5可以看出:轨道倾角为60°时同样可以得到与轨道倾角为0°和30°时相同的结论。轨道倾角60°时的拟合误差大于轨道倾角0°时的拟合误差,与30°倾角时近似。

图5 轨道倾角60°、拟合时长60min时不同轨道高度对18参数和19参数模型的拟合效果Fig.5 Fitting results of different orbit altitudes on 18-parameter and 19-parameter models when inclination angle is 60° and fitting time is 60 minutes

表4 轨道倾角为60°时不同轨道高度拟合结果的RMSETable 4 RMSEs of fitting results for different orbit altitudes at an inclination angle of 60° m

从表5和图6可以看出:轨道倾角为90°时同样可以得到与轨道倾角为0°、30°和60°时相同的结论。

图6 轨道倾角90°、拟合时长60min时不同轨道高度对18参数和19参数模型的拟合效果Fig.6 Fitting results of different orbit altitude on 18-parameter and 19-parameter models when inclination angle is 90° and fitting time is 60 minutes

表5 轨道倾角为90°时不同轨道高度拟合结果RMSETable 5 RMSEs of fitting results for different orbit altitudes at an inclination angle of 90° m

总的来看,19参数模型对不同轨道高度(500～1000km)、不同轨道倾角的卫星在不同拟合弧段时长下的拟合均具有很好的适用性;同时,19参数模型的精度要高于18参数模型。

进一步对比19参数模型和26参数模型的长弧段(4h)位置拟合精度,结果如表6所示。从表6可以看出:当轨道倾角不小于30°时,在12种不同轨道高度和倾角的情况中,有8种情况是19参数模型精度更高,另外4种情况是26参数模型精度更高。这证明了增加轨道倾角变化率这一参数对长弧段拟合精度的提升是有显著帮助的,使得19参数模型对于大倾角卫星的长弧段拟合达到了与26参数模型近似的精度。同时,26参数模型在有些情况下拟合精度低于18参数模型,这是因为26参数模型是由16参数模型发展而来的,其没有考虑18参数模型相较16参数模型增加的轨道高度变化率和角速度变化率。

表6 不同轨道高度、不同倾角的4h拟合结果RMSETable 6 RMSEs of fitting results for different orbit altitudes and inclination angles when fitting time is 4 hours m

本文在GPS的18参数模型基础上增加了偏心率变化率这一参数,仿真结果证明这一参数对于提高低轨导航增强卫星的位置拟合精度具有重要作用。如何在本文提出的19参数模型基础上,在不增加过多参数前提下进一步提高位置拟合精度,值得进一步研究。表6为此提供了重要方向,即:对于小倾角(0°)卫星的长弧段拟合,26参数模型的位置拟合精度远高于19参数模型;而对于大倾角(不小于30°)卫星,二者精度则近似。因此,导致26参数模型对小倾角卫星的位置拟合精度较高的主要参数十分值得研究,进而可以将这些起主要作用的参数与19参数模型结合,以提高后者对于小倾角卫星的位置拟合精度。