盧 祥 袁運(yùn)彬 蔣振偉

1)大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430077

2)中國(guó)科學(xué)院大學(xué)，北京 100049

基于龍格庫(kù)塔法的GLONASS軌道仿真研究*

盧 祥1，2)袁運(yùn)彬1)蔣振偉1，2)

1)大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430077

2)中國(guó)科學(xué)院大學(xué)，北京 100049

將龍格庫(kù)塔法引入GLONASS軌道仿真，研究步長(zhǎng)參數(shù)設(shè)置對(duì)結(jié)果的影響，定量分析仿真結(jié)果與廣播星歷的差異和算法效率，并進(jìn)行長(zhǎng)時(shí)間仿真試驗(yàn)分析。結(jié)果表明，步長(zhǎng)參數(shù)設(shè)置為1～300 s時(shí)能兼顧算法效率與結(jié)果精度;仿真結(jié)果與廣播星歷的差異隨時(shí)間延長(zhǎng)而增大，仿真300 min時(shí)X、Y、Z方向差異分別為266.43、246.13、336.06 m，每仿真一個(gè)歷元24顆星數(shù)據(jù)耗時(shí)約2 ms。龍格庫(kù)塔法適用于GLONASS軌道仿真。

龍格庫(kù)塔法;GLONASS;軌道仿真;積分步長(zhǎng);積分時(shí)長(zhǎng)

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

1 基于龍格庫(kù)塔法的GLONASS軌道仿真方法

采用簡(jiǎn)化GLONASS軌道力學(xué)模型進(jìn)行龍格庫(kù)塔積分，簡(jiǎn)化力學(xué)模型如下:

根據(jù)以上微分方程，采用四階龍格庫(kù)塔積分對(duì)GLONASS軌道進(jìn)行積分，積分模型為［7］:

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

由初始位置X0=X(t0=tb)計(jì)算時(shí)刻t的衛(wèi)星位置，必須以積分步長(zhǎng)h重復(fù)積分，t1=t0+h，t2=t1+h，…，tm=tm－1+h，其中 t－ h≤tm＜ t0如果 tm=th，最后一次積分步長(zhǎng)保持h不變;如果tm＞t－h(huán)，則最后一次積分步長(zhǎng)h=t－tm。本文在計(jì)算過(guò)程中，將x··LS、y··LS、z··LS固定為常數(shù)。

2 仿真試驗(yàn)及分析

2.1 步長(zhǎng)參數(shù)選擇

利用龍格庫(kù)塔積分法仿真GLONASS軌道時(shí)，步長(zhǎng)參數(shù)的選擇直接影響到仿真效率與仿真精度。較短的步長(zhǎng)仿真精度較高，但積分效率較低;較長(zhǎng)的步長(zhǎng)仿真精度較差，但積分效率較高。選擇合適的積分步長(zhǎng)可以在保證積分效率的同時(shí)達(dá)到較好的精度。為分析龍格庫(kù)塔積分方法中步長(zhǎng)選擇對(duì)仿真結(jié)果產(chǎn)生的影響，選取2012-02-20由IGS提供的GLONASS廣播星歷進(jìn)行試驗(yàn)。

以廣播星歷中1號(hào)星為例進(jìn)行試驗(yàn)，以0:15:00 數(shù)據(jù)為起始數(shù)據(jù)(圖 1)，分別以 1、10、30、60、90、120、300、900 s為積分步長(zhǎng)，積分 30 min，得到衛(wèi)星位置 Xn、速度值(即 0:45:00 的數(shù)據(jù))，將 Xn、與廣播星歷中0:45:00的衛(wèi)星位置、速度作差，計(jì)算得到不同積分步長(zhǎng)下的仿真結(jié)果與廣播星歷的差異如圖2所示。

圖1 試驗(yàn)示意圖Fig.1 Experimental sketch

圖2顯示，積分步長(zhǎng)在1～300 s變化時(shí)，仿真30 min軌道數(shù)據(jù)，與真實(shí)軌道差異均保持在1 m左右;但當(dāng)積分步長(zhǎng)升至900 s時(shí)，誤差顯著增長(zhǎng)，X、Y、Z 方向分別為 29.31、26.48、18.92 m。因此，在GLONASS軌道仿真中，宜將積分步長(zhǎng)選擇在1～300 s之間，具體選擇還要結(jié)合軌道仿真要求的數(shù)據(jù)頻率確定。

圖2 不同步長(zhǎng)30 min積分結(jié)果與廣播星歷差異Fig.2 Difference between integration results and broadcast ephemeris with different step length

2.2 仿真結(jié)果與廣播星歷的差異

為分析利用龍格庫(kù)塔法仿真GLONASS軌道的效果，設(shè)計(jì)試驗(yàn)定量研究仿真結(jié)果與廣播星歷的差異。

圖3 試驗(yàn)示意圖Fig.3 Experimental sketch

圖4 不同積分時(shí)長(zhǎng)積分結(jié)果與廣播星歷差異Fig.4 Difference between integration results and broadcast ephemeris with different time duration

在2012-02-20 IGS發(fā)布的GLONASS廣播星歷中以1號(hào)星0:15:00為起點(diǎn)(圖3)，以60 s為積分步長(zhǎng)分別積分 30 min、60 min、90 min、…、300 min，以廣播星歷中對(duì)應(yīng)時(shí)刻的軌道位置作為參考值，計(jì)算得到不同積分時(shí)長(zhǎng)下的軌道差異如圖4所示。

隨著時(shí)間的延長(zhǎng)，仿真結(jié)果與廣播星歷差異逐漸增大。30 min時(shí)，X、Y、Z三個(gè)方向的差異分別為1.25、0.75、1.21 m，優(yōu)于文獻(xiàn)［3］采用 4、5、6、7 階多項(xiàng)式仿真GPS系統(tǒng)軌道30 min的精度(10 m級(jí))，與文獻(xiàn)［3］采用20階多項(xiàng)式擬合法仿真30 min GPS系統(tǒng)的精度相當(dāng)(m級(jí))。90 min仿真結(jié)果與廣播星歷的差異增至10 m級(jí)，到300 min時(shí)X、Y、Z 三個(gè)分量差異分別為 266.4、246.1、336.0 m，相比GLONASS衛(wèi)星軌道高度19 100 km，這些差異非常微小，可認(rèn)為與廣播星歷基本吻合，表明龍格庫(kù)塔法適用于GLONASS軌道仿真。

2.3 算法效率

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

2.4 長(zhǎng)時(shí)間仿真穩(wěn)定性

使用龍格庫(kù)塔法進(jìn)行長(zhǎng)時(shí)間軌道仿真時(shí)，仿真結(jié)果發(fā)散性、算法穩(wěn)定性都是關(guān)鍵問(wèn)題。為分析龍格庫(kù)塔法進(jìn)行GLONASS軌道長(zhǎng)時(shí)間仿真的穩(wěn)定性，以2012-02-20 00:15:00數(shù)據(jù)為起點(diǎn)進(jìn)行GLONASS星座(24顆星)軌道仿真，得到 1 d、2 d、20 d的仿真軌道。選取1號(hào)星為例，繪制廣播星歷、仿真軌道1 d、2 d、20 d軌道軌跡圖(圖6)，將20天1號(hào)星的仿真軌道與廣播星歷作差得到20 d仿真軌道與廣播星歷X、Y、Z方向差異(圖7)。

由圖6與7可知，利用龍格庫(kù)塔法仿真所得的GLONASS軌道軌跡與廣播星歷軌跡高度相似;使用龍格庫(kù)塔法進(jìn)行長(zhǎng)時(shí)間GLONASS軌道仿真時(shí)，盡管仿真軌道形狀與廣播星歷類(lèi)似，但仿真軌道與廣播星歷的差值隨時(shí)間延長(zhǎng)逐漸增大。仿真時(shí)間到達(dá)480 h(20 d)時(shí)，X、Y、Z三個(gè)方向與廣播星歷的差異分別為 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 結(jié)論

系統(tǒng)研究了利用龍格庫(kù)塔法實(shí)現(xiàn)GLONASS軌道仿真的性能與涉及的關(guān)鍵技術(shù)，給出龍格庫(kù)塔法的數(shù)學(xué)模型。研究結(jié)果顯示:

1)步長(zhǎng)參數(shù)設(shè)置為1～300 s時(shí)能兼顧算法效率與結(jié)果精度;

2)以60 s為積分步長(zhǎng)，30 min仿真軌道與廣播星歷 X、Y、Z 三個(gè)方向差異分別為 1.25、0.75、1.21 m，300 min仿真軌道與廣播星歷X、Y、Z分量差異分別為 266.4、246.1、336.0 m;

3)利用龍格庫(kù)塔法每仿真一個(gè)歷元24顆星數(shù)據(jù)耗時(shí)約為2 ms，可滿足接收機(jī)仿真測(cè)試工作中常用的50 Hz高性能數(shù)據(jù)仿真要求;

4)使用龍格庫(kù)塔法仿真20 d GLONASS軌道，仿真軌道形狀與廣播星歷類(lèi)似，但仿真軌道與廣播星歷的差值隨時(shí)間的延長(zhǎng)逐漸增大，20 d時(shí)X、Y、Z三個(gè)方向與廣播星歷的差異分別為44 906.55、43 388.11、110 309.45 m。

1 趙春梅，歐吉坤，文援蘭.基于GALILEO及GPS-GALILEO組合系統(tǒng)的仿真分析［J］.系統(tǒng)仿真學(xué)報(bào)，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)的單點(diǎn)定位研究［J］.系統(tǒng)仿真學(xué)報(bào)，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)星軌道仿真與可視化實(shí)現(xiàn)［D］.哈爾濱:哈爾濱工程大學(xué)，2009.(Bai Wenjuan.Simulation and visualization realization of GPS orbit based on PpenGL［D］.Harbin:Harbin Engineering University，2009)

4 戴沖.衛(wèi)星導(dǎo)航系統(tǒng)空間段仿真關(guān)鍵技術(shù)研究［D］.長(zhǎng)沙:國(guó)防科技大學(xué)，2011.(Dai Chong.Study on key technology of space segment of satellite navigation system［D］.Changsha:National University of Defense Technology，2011)

5 謝杰，等.GPS信號(hào)仿真器方案設(shè)計(jì)與實(shí)現(xiàn)［J］.計(jì)算機(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］.計(jì)算機(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 李慶楊，王能超，易大義.數(shù)值分析［M］.北京:清華大學(xué)出版社，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

國(guó)家自然科學(xué)基金重點(diǎn)項(xiàng)目(41231064);國(guó)家863計(jì)劃項(xiàng)目(2012AA121803)。

盧祥，男，1989年生，碩士生，主要研究方向?yàn)镚NSS精密定位。E－mail:luxiang.whigg@gmail.com。