盧 祥 袁運(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)定性。
采用簡(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ù)。
利用龍格庫(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
為分析利用龍格庫(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軌道仿真。
接收機(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的仿真要求。
使用龍格庫(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
系統(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。