基于龙格库塔法的GLONASS轨道仿真研究*

卢 祥 袁运彬 蒋振伟

1)大地测量与地球动力学国家重点实验室，武汉 430077

2)中国科学院大学，北京 100049

基于龙格库塔法的GLONASS轨道仿真研究*

卢 祥1，2)袁运彬1)蒋振伟1，2)

1)大地测量与地球动力学国家重点实验室，武汉 430077

2)中国科学院大学，北京 100049

将龙格库塔法引入GLONASS轨道仿真，研究步长参数设置对结果的影响，定量分析仿真结果与广播星历的差异和算法效率，并进行长时间仿真试验分析。结果表明，步长参数设置为1～300 s时能兼顾算法效率与结果精度;仿真结果与广播星历的差异随时间延长而增大，仿真300 min时X、Y、Z方向差异分别为266.43、246.13、336.06 m，每仿真一个历元24颗星数据耗时约2 ms。龙格库塔法适用于GLONASS轨道仿真。

龙格库塔法;GLONASS;轨道仿真;积分步长;积分时长

GNSS(GPS、BD、GALLEO、GLONASS)仿真系统在接收机测试、算法验证等工作中有着重要作用［1］。与 GPS、BD、GALILEO 系统仿真情况不同，国内外对GLONASS仿真研究较少。Spirent公司于2007年实现了GLONASS的系统仿真，但其方法及代码均未公开。国内开展GNSS系统仿真较早的单位实现了GPS、GALILEO、BD 3大系统的仿真，但鲜见公开报道GLONASS系统仿真研究。赵春梅等运用动力学定轨法进行了GALILEO轨道仿真，但由于GLONASS轨道采用卫星坐标、速度、加速度方式发布，不同于其他3大系统，传统 GPS、BD、GALILEO中的动力学定轨方法不能直接应用于GLONASS轨道仿真。传统动力学定轨方法公式繁杂、计算量大，通常难以满足 GLONASS高性能仿真的要求［2－6］。本文首次将龙格库塔法引入到GLONASS轨道仿真，分析步长、仿真时间对仿真结果的影响，定量分析算法效率，并测试算法的稳定性。

1 基于龙格库塔法的GLONASS轨道仿真方法

采用简化GLONASS轨道力学模型进行龙格库塔积分，简化力学模型如下:

根据以上微分方程，采用四阶龙格库塔积分对GLONASS轨道进行积分，积分模型为［7］:

根据式(7)～(11)积分模型，可按以下步骤进行积分:

由初始位置X0=X(t0=tb)计算时刻t的卫星位置，必须以积分步长h重复积分，t1=t0+h，t2=t1+h，…，tm=tm－1+h，其中 t－ h≤tm＜ t0如果 tm=th，最后一次积分步长保持h不变;如果tm＞t－h，则最后一次积分步长h=t－tm。本文在计算过程中，将x··LS、y··LS、z··LS固定为常数。

2 仿真试验及分析

2.1 步长参数选择

利用龙格库塔积分法仿真GLONASS轨道时，步长参数的选择直接影响到仿真效率与仿真精度。较短的步长仿真精度较高，但积分效率较低;较长的步长仿真精度较差，但积分效率较高。选择合适的积分步长可以在保证积分效率的同时达到较好的精度。为分析龙格库塔积分方法中步长选择对仿真结果产生的影响，选取2012-02-20由IGS提供的GLONASS广播星历进行试验。

以广播星历中1号星为例进行试验，以0:15:00 数据为起始数据(图 1)，分别以 1、10、30、60、90、120、300、900 s为积分步长，积分 30 min，得到卫星位置 Xn、速度值(即 0:45:00 的数据)，将 Xn、与广播星历中0:45:00的卫星位置、速度作差，计算得到不同积分步长下的仿真结果与广播星历的差异如图2所示。

图1 试验示意图Fig.1 Experimental sketch

图2显示，积分步长在1～300 s变化时，仿真30 min轨道数据，与真实轨道差异均保持在1 m左右;但当积分步长升至900 s时，误差显著增长，X、Y、Z 方向分别为 29.31、26.48、18.92 m。因此，在GLONASS轨道仿真中，宜将积分步长选择在1～300 s之间，具体选择还要结合轨道仿真要求的数据频率确定。

图2 不同步长30 min积分结果与广播星历差异Fig.2 Difference between integration results and broadcast ephemeris with different step length

2.2 仿真结果与广播星历的差异

为分析利用龙格库塔法仿真GLONASS轨道的效果，设计试验定量研究仿真结果与广播星历的差异。

图3 试验示意图Fig.3 Experimental sketch

图4 不同积分时长积分结果与广播星历差异Fig.4 Difference between integration results and broadcast ephemeris with different time duration

在2012-02-20 IGS发布的GLONASS广播星历中以1号星0:15:00为起点(图3)，以60 s为积分步长分别积分 30 min、60 min、90 min、…、300 min，以广播星历中对应时刻的轨道位置作为参考值，计算得到不同积分时长下的轨道差异如图4所示。

随着时间的延长，仿真结果与广播星历差异逐渐增大。30 min时，X、Y、Z三个方向的差异分别为1.25、0.75、1.21 m，优于文献［3］采用 4、5、6、7 阶多项式仿真GPS系统轨道30 min的精度(10 m级)，与文献［3］采用20阶多项式拟合法仿真30 min GPS系统的精度相当(m级)。90 min仿真结果与广播星历的差异增至10 m级，到300 min时X、Y、Z 三个分量差异分别为 266.4、246.1、336.0 m，相比GLONASS卫星轨道高度19 100 km，这些差异非常微小，可认为与广播星历基本吻合，表明龙格库塔法适用于GLONASS轨道仿真。

2.3 算法效率

接收机进行野外测试前，需进行大量的仿真测试，特别是检测50 Hz甚至更高频的数据时，对轨道仿真算法效率有很高的要求。定量研究龙格库塔法的效率，选取不同时长，积分步长设置为1 s，仿真GLONASS星座(24颗星)轨道数据(表1)。图5为算法效率图，龙格库塔法每仿真一个历元24颗星数据耗时约2 ms，可满足50 Hz甚至更高频率GLONASS的仿真要求。

2.4 长时间仿真稳定性

使用龙格库塔法进行长时间轨道仿真时，仿真结果发散性、算法稳定性都是关键问题。为分析龙格库塔法进行GLONASS轨道长时间仿真的稳定性，以2012-02-20 00:15:00数据为起点进行GLONASS星座(24颗星)轨道仿真，得到 1 d、2 d、20 d的仿真轨道。选取1号星为例，绘制广播星历、仿真轨道1 d、2 d、20 d轨道轨迹图(图6)，将20天1号星的仿真轨道与广播星历作差得到20 d仿真轨道与广播星历X、Y、Z方向差异(图7)。

由图6与7可知，利用龙格库塔法仿真所得的GLONASS轨道轨迹与广播星历轨迹高度相似;使用龙格库塔法进行长时间GLONASS轨道仿真时，尽管仿真轨道形状与广播星历类似，但仿真轨道与广播星历的差值随时间延长逐渐增大。仿真时间到达480 h(20 d)时，X、Y、Z三个方向与广播星历的差异分别为 44 906.55、43 388.11、110 309.45 m。

表1 算法效率Tab.1 Efficiency of the algorithm

图5 算法效率Fig.5 Efficiency of the algorithm

3 结论

系统研究了利用龙格库塔法实现GLONASS轨道仿真的性能与涉及的关键技术，给出龙格库塔法的数学模型。研究结果显示:

1)步长参数设置为1～300 s时能兼顾算法效率与结果精度;

2)以60 s为积分步长，30 min仿真轨道与广播星历 X、Y、Z 三个方向差异分别为 1.25、0.75、1.21 m，300 min仿真轨道与广播星历X、Y、Z分量差异分别为 266.4、246.1、336.0 m;

3)利用龙格库塔法每仿真一个历元24颗星数据耗时约为2 ms，可满足接收机仿真测试工作中常用的50 Hz高性能数据仿真要求;

4)使用龙格库塔法仿真20 d GLONASS轨道，仿真轨道形状与广播星历类似，但仿真轨道与广播星历的差值随时间的延长逐渐增大，20 d时X、Y、Z三个方向与广播星历的差异分别为44 906.55、43 388.11、110 309.45 m。

1 赵春梅，欧吉坤，文援兰.基于GALILEO及GPS-GALILEO组合系统的仿真分析［J］.系统仿真学报，2005，17(4):1 008 －1 011.(Zhao Chunmei，Ou Jikun，Wen Yuanlan.Simulation analysis based on the combination of GALILEO and integrated GPS-GALILEO［J］.Journal of System Simulation，2005，17(4):1 008 －1 011)

2 李得海，等.基于GALELIO/GPS仿真系统的单点定位研究［J］.系统仿真学报，2009，21(10):2 828 －2 831.(Li Dehai，et al.Research on point positioning based on simulated GALILEO and GPS system［J］.Journal of System Simulation，2009，21(10):2 828 －2 831)

3 白文娟.基于OpenGL的GPS卫星轨道仿真与可视化实现［D］.哈尔滨:哈尔滨工程大学，2009.(Bai Wenjuan.Simulation and visualization realization of GPS orbit based on PpenGL［D］.Harbin:Harbin Engineering University，2009)

4 戴冲.卫星导航系统空间段仿真关键技术研究［D］.长沙:国防科技大学，2011.(Dai Chong.Study on key technology of space segment of satellite navigation system［D］.Changsha:National University of Defense Technology，2011)

5 谢杰，等.GPS信号仿真器方案设计与实现［J］.计算机仿真，2012，29(2):36 －39.(Xie Jie，et al.Design and implementation of GPS signal simulator based on ARM［J］.Computer Simulation，2012，29(2):36 －39)

6 于清德，等.GNSS几何性能仿真分析［J］.计算机仿真，2010，27(6):87 －91.(Yu Qingde，et al.simulation analysis of GNSS geometric performance［J］.Computer Simulation，2010，27(6):87 －91)

7 李庆杨，王能超，易大义.数值分析［M］.北京:清华大学出版社，2008.(Li Qingyang，Wang Nengchao，Yi Dayi.Numerical analysis［M］.Beijing:Tsinghua University Press，2008)

RESEARCH ON GLONASS ORBIT SIMULATION BASED ON RUNGE-KUTTA METHOD

Lu Xiang1，2)，Yuan Yunbin1)and Jiang Zhenwei1，2)
1)State Key Laboratory of Geodesy and Earth’s Dynamics，Wuhan 430077
2)University of Chinese Academy of Sciences，Beijing100049

Runge-Kutta method was introduced to GLONASS orbit simulation to analyze how the step length parameter effects on simulated results，the efficiency of Runge-Kutta method and the difference between simulated results and broadcast orbit.The result shows that step the balance between algorithm efficiency and accuracy can be keeped though setting step length as 1 s－300 s.The difference between simulated results and broadcast ephemeris increases with time，and reaches about 250m in X，Y，Z directions for 300 minutes simulation time duration.It takes 2 ms to simulate 1 epoch 24 satellites orbit data.According the research，Runge-Kutta method can be used for GLONASS orbit simulation.

Runge-Kutta;GLONASS;orbit simulation;integration step length;integration time duration

P228

A

1671-5942(2014)03-0137-05

2013-10-13

国家自然科学基金重点项目(41231064);国家863计划项目(2012AA121803)。

卢祥，男，1989年生，硕士生，主要研究方向为GNSS精密定位。E－mail:luxiang.whigg@gmail.com。