刘 岩,饶才杰,吕 达,白志强
(中国航天科工信息技术研究院,北京100070)
组合GPS/GLONASS增加了观测卫星的个数,提高了卫星的可用性,改善了卫星定位的几何结构,通过融合可以提高定位精度和可用性,解决GPS网络RTK定位中困难地区的定位问题,具有重要应用价值[1-2]。因此,研究 GLONASS系统定位尤其是组合定位有重要实用价值。如何根据GLONASS广播星历进行卫星位置计算是组合定位的重要一环。GLONASS广播导航电文没有给出开普勒轨道参数,对于一个给定的时刻,不能通过最新的导航电文参考时刻的开普勒轨道参数直接计算卫星位置。GLONASS导航电文给出的是PZ-90坐标系下参考时刻卫星的运动状态向量,包括参考时刻卫星的位置、速度、日月摄动加速度等,每30min更新一次,这决定了对于某一时刻的卫星位置,需要以上述的卫星运动状态向量为初值,采用受力模型通过积分得到[3]。
大量的国内外文献对GPS卫星位置的计算方法进行了阐述,但对GLONASS卫星位置计算方法的研究较少。目前GLONASS卫星轨道数值积分方法主要有龙格-库塔法(Runge-Kutta),其中以四阶龙格-库塔法最为常用。文献[4]等推导了GLONASS卫星在地固坐标系中的运动方程,并简要介绍了利用龙格-库塔法进行轨道积分的数值结果;文献[5]对龙格-库塔法进行了分析;文献[6]利用龙格-库塔法实现了GLONASS卫星轨道求解。文献[7]完成了定步长龙格-库塔法的计算程序。但是,龙格-库塔法在计算某一时刻卫星位置和速度时只用到前一步的信息,为了保证定位精度,需要多次重新计算多个点处的函数值[8],计算量较大。采用Adams线性多步预测校正法来实现GLONASS卫星位置解算,该方法通过较多地利用前面几步的已知信息,在保证较高精度的前提下,运算量降低为龙格-库塔法的50%,同时稳定性更好。实验结果表明,Adams线性多步预测校正法是一种求解GLONASS位置的有效方法。
GLONASS卫星在PZ-90坐标系下的运动方程为[3]
式中,x,y,z分别表示卫星在PZ-90坐标系下的坐标;vx,vy,vz分别表示三个方向的卫星速度,其通式为而V={v,v,v}为卫星的速度向xyz量,dX={dx,dy,dz).
式中:地球引力常量μ=398 600.44km3/s2;卫星到地球质心的距离地球的长半轴ae=6 378.136km;地球重力位第二带谐系数J20=0.001 082 63;地球自转速度ω=0.000 072
为日月摄动加速度。t时刻式(2)的函数形式可以表示为
其通式为
式中:a为加速度向量;X={x,y,z}为卫星的坐标向量;为日月的摄动加速度向量。
GLONASS每30min更新一组新的星历参数,某一时刻ti的卫星位置和速度需要采用数值积分的方法求解,要以最新星历的参考时刻tb为参考时刻,并以该时刻卫星的速度和坐标为起算数据,对加速度和速度进行积分,其积分通式为
P阶龙格-库塔法原始公式如下:
使用龙格-库塔法由yn计算un+1时,只用到前一步的信息yn,为了提高运算精度,需要重新计算多次函数值f(x),由于解决实际问题所花费的代价主要取决于求解函数值的个数,所以运算量较大。
经过一个步长后,使用四阶龙格-库塔积分法(以下简称R-K法)计算GLONASS卫星的速度和位置坐标为
R-K法属于单步法,在计算时刻GLONASS卫星的速度Vi+1和位置坐标Xi+1时,只用到前一步的信息Vi和Xi,为了提高运算精度,重新计算了4次加速度值a1,a2,a3,a4,由于进行轨道积分所花费的代价主要在于求解加速度,如式(2)所示,因此,使用R-K方法进行GLONASS卫星位置解算运算量较大。
针对龙格-库塔方法进行轨道积分运算量大的缺点,研究应用Adams线性多步法来进行GLONASS卫星位置解算,通过较多地利用前面几步的已知信息,构造计算量小、精度高的GLONASS卫星位置解算的新算法。Adams线性多步公式包括隐式公式和显式公式,同阶的隐式公式比显式公式精确,而且稳定性也好。Adams线性多步法公式为
当βk=0为显式公式,当βk≠0为隐性公式。通常将四步四阶显式公式作为预测公式,将三步四阶隐式公式作为校正公式联合使用,如式(10)所示。选用四步五阶隐式公式作为校正公式,精度提高一阶,本文选用的预测校正公式如式(11)所示:
采用Adams线性多步预测校正法(以下简称Adams法)计算tt+1时刻的速度Vi+1和位置坐标Xi+1公式如下:
通过显式公式对ti+1时刻GLONASS卫星位置和速度值进行预测,如式(12)所示:
通过隐式公式对ti+1时刻GLONASS卫星位置和速度预测值进行校正,从而得到准确值,如式(13)所示:
采用Adams法计算ti+1时刻的速度Vi+1和位置坐标 Xi+1时,需要用到前4步的信息Vi,Vi-1,Vi-2,Vi-3和 Xi,Xi-1,Xi-2,Xi-3,通过较多地利用已知信息,实现了运算量小,高精度的解算。Adams线性多步预测校正法解算GLONASS卫星位置的特点如下:
1)Adams算法的计算量大约是R-K方法的50%.因为采用 Adams法计算ti+1时刻的速度Vi+1和位置坐标Xi+1时只需计算两次加速度值ai,ai+1,而采用R-K法时需计算4次加速度值a1,a2,a3,a4,由于计算加速度是GLONASS轨道积分过程中运算量较大的部分,由式(2)可见,所以Adams算法可以在很大程度上降低运算量。
2)Adams算法比R-K算法截断误差高一阶。R-K法具有4阶精度,截断误差为O(h2);Adams方法显示公式具有4阶精度,截断误差O(h5)为;本文采用的隐式公式具有5阶精度,截断误差为O(h6)[8].
选取2009年8月29日23∶45至2009年8月31日16∶45每间隔半小时共82个时刻,GLONASS 2号卫星的广播星历作为实验数据。
分别使用Adams方法与R-K方法,以2009年8月29日23∶45至2009年8月31日15∶45每半小时间隔的参考时间向前轨道积分30min,得到2009年8月30日00∶15至2009年8月31日16∶15半小时间隔的共81个时刻的,GLONASS 2号卫星位置的X、Y、Z坐标,再将计算结果与广播星历里给出的卫星坐标进行比较,得到两者差值的统计结果如表1所示。
分别使用Adams方法与R-K方法,以2009年8月29日23∶45至2009年8月31日15∶15每半小时间隔的参考时间向前轨道积分60min,得到2009年8月30日00∶45至2009年8月31日16∶15半小时间隔的共80个时刻的GLONASS 2号卫星位置的X、Y、Z坐标,再将计算结果与广播星历里给出的卫星坐标进行比较,得到两 者差值的统计结果如表2所示。
表1 30min轨道积分结果与星历所给坐标之差的绝对值/m
表2 60min轨道积分结果与星历所给坐标之差的绝对值/m
对试验数据分析如下:
1)向前积分30min,选取不同步长对两种算法的误差影响不大,两种算法误差都比较小,但Adams算法比R-K算法误差的均值要小。向前积分60min,Adams算法优势更加明显,如以200为步长积分时,R-K算法距离误差均值为8.971 7,而Adams法为8.689 3.
2)随着积分时间延长,两种算法的误差随之增大,但Adams算法的均值误差小于R-K算法,且随着积分时间越长越明显,说明Adams算法稳定性[8]优于 R-K 算法。
3)采用R-K方法进行一步积分,需要计算4次三维方向的加速度;采用Adams进行一步运算,只需要计算2次;由于轨道积分运算的代价主要取决于计算加速度的次数,所以Adams算法的运算量约为R-K算法的50%.
引入了Adams线性多步预测校正法进行GLONASS卫星位置解算,并采用GLONASS的广播星历对Adams算法和R-K算法的积分结果和运算效率进行分析,证明了所提方法是一种有效的GLONASS卫星位置解算方法。与四阶R-K算法相比,Adams算法精度高、稳定性好,而且运算量降低了50%.采用Adams法进行GLONASS位置解算,对确保多系统兼容GNSS接收机的精度、实时性等方面有重要意义。
[1]BRUYNINX C.Comparing GPS-only with GPS+GLONASS positioning in a regional permanent GNSS network[J].GPS Solutions,2007,11(2):97-106.
[2]ZINOVIEV A E.Using GLONASS in combined GNSS receivers:current status[C]//ION GNSS 18th International Technical Meeting of the Satellite Division.Long Beach,CA,USA,2005:1046-1057.
[3]ICD-GLONASS.GLONASS interface control document version 5.1,ICD 2008E[R].Moscow:Coordination Scientific Information Center,2008.
[4]葛茂荣,过静珺,葛胜杰.GLONASS卫星坐标的计算方法[J].测绘通报,1999(2):2-4.
[5]杨 剑,王泽民,孟 泱,等.GLONASS卫星轨道积分算法分析[J].武汉大学学报·信息科学版,2006,31(7):613-615.
[6]李建文.GLONASS卫星导航系统及 GPS/GLONASS组合应用研究[D].郑州:信息工程大学,2001.
[7]葛 奎,王解先.GLONASS卫星位置计算与程序实现[J].测绘与空间地理信息,2009,32(2):137-140.
[8]丁丽娟,程杞元.数值计算方法[M].北京:北京理工大学出版社,2008.