盧 祥 袁運彬 蔣振偉

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)仿真系統(tǒng)在接收機測試、算法驗證等工作中有著重要作用［1］。與 GPS、BD、GALILEO 系統(tǒng)仿真情況不同，國內外對GLONASS仿真研究較少。Spirent公司于2007年實現了GLONASS的系統(tǒng)仿真，但其方法及代碼均未公開。國內開展GNSS系統(tǒng)仿真較早的單位實現了GPS、GALILEO、BD 3大系統(tǒng)的仿真，但鮮見公開報道GLONASS系統(tǒng)仿真研究。趙春梅等運用動力學定軌法進行了GALILEO軌道仿真，但由于GLONASS軌道采用衛(wèi)星坐標、速度、加速度方式發(fā)布，不同于其他3大系統(tǒng)，傳統(tǒng) GPS、BD、GALILEO中的動力學定軌方法不能直接應用于GLONASS軌道仿真。傳統(tǒng)動力學定軌方法公式繁雜、計算量大，通常難以滿足 GLONASS高性能仿真的要求［2－6］。本文首次將龍格庫塔法引入到GLONASS軌道仿真，分析步長、仿真時間對仿真結果的影響，定量分析算法效率，并測試算法的穩(wěn)定性。

1 基于龍格庫塔法的GLONASS軌道仿真方法

采用簡化GLONASS軌道力學模型進行龍格庫塔積分，簡化力學模型如下:

根據以上微分方程，采用四階龍格庫塔積分對GLONASS軌道進行積分，積分模型為［7］:

根據式(7)～(11)積分模型，可按以下步驟進行積分:

由初始位置X0=X(t0=tb)計算時刻t的衛(wèi)星位置，必須以積分步長h重復積分，t1=t0+h，t2=t1+h，…，tm=tm－1+h，其中 t－ h≤tm＜ t0如果 tm=th，最后一次積分步長保持h不變;如果tm＞t－h(huán)，則最后一次積分步長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，得到衛(wèi)星位置 Xn、速度值(即 0:45:00 的數據)，將 Xn、與廣播星歷中0:45:00的衛(wèi)星位置、速度作差，計算得到不同積分步長下的仿真結果與廣播星歷的差異如圖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發(fā)布的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，優(yōu)于文獻［3］采用 4、5、6、7 階多項式仿真GPS系統(tǒng)軌道30 min的精度(10 m級)，與文獻［3］采用20階多項式擬合法仿真30 min GPS系統(tǒng)的精度相當(m級)。90 min仿真結果與廣播星歷的差異增至10 m級，到300 min時X、Y、Z 三個分量差異分別為 266.4、246.1、336.0 m，相比GLONASS衛(wèi)星軌道高度19 100 km，這些差異非常微小，可認為與廣播星歷基本吻合，表明龍格庫塔法適用于GLONASS軌道仿真。

2.3 算法效率

接收機進行野外測試前，需進行大量的仿真測試，特別是檢測50 Hz甚至更高頻的數據時，對軌道仿真算法效率有很高的要求。定量研究龍格庫塔法的效率，選取不同時長，積分步長設置為1 s，仿真GLONASS星座(24顆星)軌道數據(表1)。圖5為算法效率圖，龍格庫塔法每仿真一個歷元24顆星數據耗時約2 ms，可滿足50 Hz甚至更高頻率GLONASS的仿真要求。

2.4 長時間仿真穩(wěn)定性

使用龍格庫塔法進行長時間軌道仿真時，仿真結果發(fā)散性、算法穩(wěn)定性都是關鍵問題。為分析龍格庫塔法進行GLONASS軌道長時間仿真的穩(wěn)定性，以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 結論

系統(tǒng)研究了利用龍格庫塔法實現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組合系統(tǒng)的仿真分析［J］.系統(tǒng)仿真學報，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仿真系統(tǒng)的單點定位研究［J］.系統(tǒng)仿真學報，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衛(wèi)星軌道仿真與可視化實現［D］.哈爾濱:哈爾濱工程大學，2009.(Bai Wenjuan.Simulation and visualization realization of GPS orbit based on PpenGL［D］.Harbin:Harbin Engineering University，2009)

4 戴沖.衛(wèi)星導航系統(tǒng)空間段仿真關鍵技術研究［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。