導(dǎo)航衛(wèi)星軌道擬合方法與仿真研究

萬家歡，莊春華，陳秀萬，李智慧，張威奕

(1.北京大學(xué)遙感與地理信息系統(tǒng)研究所，北京100871;2.北京環(huán)球信息中心，北京100094;3.92292部隊，山東青島266405)

導(dǎo)航衛(wèi)星軌道擬合方法與仿真研究

萬家歡1，3，莊春華2，陳秀萬1，李智慧1，張威奕1

(1.北京大學(xué)遙感與地理信息系統(tǒng)研究所，北京100871;2.北京環(huán)球信息中心，北京100094;3.92292部隊，山東青島266405)

提出利用虛擬現(xiàn)實技術(shù)加強對導(dǎo)航衛(wèi)星星座的三維仿真模擬，用于系統(tǒng)的設(shè)計與監(jiān)控。考慮到GPS系統(tǒng)的成熟性和與北斗系統(tǒng)的相似性，以GPS精密星歷為基礎(chǔ)，深入研究利用GPS精密星歷擬合衛(wèi)星三維坐標(biāo)的方法，對兩種常用方法——拉格朗日多項式內(nèi)插和切比雪夫多項式擬合的處理效果進(jìn)行分析和比較，給出兩種擬合方法得到的衛(wèi)星位置誤差與階數(shù)的關(guān)系，探討短時間外推的效果，確定切比雪夫多項式擬合為最佳的衛(wèi)星軌道擬合方法，并據(jù)此求得任意時刻的衛(wèi)星坐標(biāo)，構(gòu)建了導(dǎo)航衛(wèi)星三維標(biāo)準(zhǔn)化軌道虛擬運行場景。

精密星歷;衛(wèi)星軌道;拉格朗日多項式;切比雪夫多項式;虛擬現(xiàn)實

一、引 言

在航空航天領(lǐng)域，衛(wèi)星信息顯示、軌道模擬、姿態(tài)異常等通常難以用數(shù)據(jù)形式表達(dá)，因此，國外紛紛建立三維仿真系統(tǒng)來研究并監(jiān)控在軌衛(wèi)星的運行狀態(tài)，國內(nèi)衛(wèi)星系統(tǒng)仿真還處于使用國外的可視化仿真工具來開發(fā)應(yīng)用程序的階段，如衛(wèi)星工具包STK(satellite tool kit)等，這些軟件價格昂貴，不公開源代碼，且在一些系統(tǒng)分析組件方面對我國用戶有很大的限制。鑒于GPS的成熟性和與北斗的相似性，利用虛擬現(xiàn)實技術(shù)加強對GPS系統(tǒng)星座的三維仿真模擬及相關(guān)分析，對加快我國新一代導(dǎo)航定位系統(tǒng)的發(fā)展具有很強的現(xiàn)實意義。

GPS精密星歷由國際衛(wèi)星導(dǎo)航服務(wù)(International GNSS Service，IGS)等國際組織發(fā)布，它給出了衛(wèi)星在ITRF中的三維坐標(biāo)及衛(wèi)星鐘差，其中三維坐標(biāo)的最終精度為±5 cm，可以滿足高精度定位的需要。但GPS精密星歷的數(shù)據(jù)采樣間隔一般為15 min，因此在GPS仿真系統(tǒng)三維軌道的生成中，用戶需選擇合理的數(shù)學(xué)模型對GPS精密星歷進(jìn)行插值或者擬合，才能獲得任意時刻衛(wèi)星的位置，這是生成虛擬衛(wèi)星軌道的基礎(chǔ)數(shù)據(jù)［1-5］。

基于GPS精密星歷的軌道逼近算法通常采用軌道擬合和坐標(biāo)內(nèi)插。其中，內(nèi)插多以拉格朗日內(nèi)插為主;擬合多以切比雪夫多項式為主，它們都能實現(xiàn)導(dǎo)航衛(wèi)星軌道擬合，但擬合精度和擬合效果不同，在實際工程中的可應(yīng)用性也不同。為了確定工程上可用的最佳衛(wèi)星軌道擬合方法和達(dá)到精度要求的擬合多項式階數(shù)，需構(gòu)建標(biāo)準(zhǔn)化的三維軌道虛擬運行場景。本文對這兩種常用方法的處理效果進(jìn)行了分析和比較，給出了兩種擬合方法得到的衛(wèi)星位置誤差與階數(shù)的關(guān)系，并進(jìn)一步探討了短時間外推的效果，最終確定了切比雪夫多項式擬合為最佳的衛(wèi)星軌道擬合方法。本文選用符合精度要求的多項式擬合階數(shù)得到了任意時刻的衛(wèi)星坐標(biāo)，在此基礎(chǔ)上開發(fā)了基于MFC和WTK的單視圖導(dǎo)航衛(wèi)星三維標(biāo)準(zhǔn)化軌道虛擬運行場景。

二、兩種擬合方法的數(shù)學(xué)模型

1.拉格朗日多項式內(nèi)插［6］

設(shè)y=f(x)是區(qū)間［a，b］上的一個實函數(shù)，xi(i=0，1，…，n)是［a，b］上n+1個互異實數(shù)，且y=f(x)在xi的值為yi=f(xi)，則區(qū)間［a，b］上任意一點x的n階拉格朗日多項式的代數(shù)表達(dá)式為

利用n階拉格朗日多項式內(nèi)插任意時刻的衛(wèi)星位置，必須至少有n+1個時刻的衛(wèi)星位置，并且內(nèi)插時刻要在已知時刻的區(qū)間內(nèi)。

2.切比雪夫多項式擬合

假設(shè)需要在時間間隔［t0，t0+Δt］計算n階切比雪夫多項式系數(shù)(其中，t0為起始?xì)v元時刻;Δt為擬合時間區(qū)間的長度)，先將變量t∈［t0，t0+Δt］變換成τ∈［－1，1］

則衛(wèi)星坐標(biāo)(X，Y，Z)的切比雪夫多項式為

式中，n為切比雪夫多項式的階數(shù);CXi、CYi、CZi分別為X坐標(biāo)分量、Y坐標(biāo)分量、Z坐標(biāo)分量的切比雪夫多項式系數(shù)，為待求量。Ti根據(jù)遞推公式(2)確定。根據(jù)所取時間段內(nèi)已知的GPS衛(wèi)星坐標(biāo)，利用最小二乘法可以計算出某一顆GPS衛(wèi)星在［t0，t0+Δt］時間區(qū)間內(nèi)的切比雪夫多項式擬合系數(shù)CXi、CYi、CZi。利用這些系數(shù)，根據(jù)式(3)可以計算出［t0，t0+ Δt］時間區(qū)間內(nèi)任意時刻的衛(wèi)星坐標(biāo)［7］。

三、精密星歷擬合方法分析和比較

在NASA官方網(wǎng)站上下載SP3精密星歷文件，本文選取起始?xì)v元在2008年5月1日0時0分0秒，終止歷元時刻在2008年5月1日23時45分0秒的SP3精密星歷文件作為試驗數(shù)據(jù)。

1.拉格朗日內(nèi)插效果分析

本文選用衛(wèi)星序號為11的歷元值作為試驗數(shù)據(jù)，取擬合的時間間隔為30 min，采用對稱內(nèi)插的方法，取30 min間隔的歷元作為內(nèi)插點計算各個內(nèi)插時間段中間時刻的衛(wèi)星位置，以便將內(nèi)插得到的衛(wèi)星位置同精密星歷給出的衛(wèi)星位置進(jìn)行比較。根據(jù)以上分析，利用VC6.0軟件編寫相應(yīng)程序，得到了11號衛(wèi)星5～18階的內(nèi)插結(jié)果(表1給出了8～18階的內(nèi)插結(jié)果)。

由表1可知，隨著階數(shù)逐漸升高，拉格朗日插值誤差逐漸減小。階數(shù)為12階時，X、Y、Z 3個坐標(biāo)的最大中誤差為±2.40 mm，最大誤差為6.05 mm，內(nèi)插的結(jié)果滿足精密定位的要求。當(dāng)階數(shù)高于12階時，插值精度基本不變。用同樣的方法內(nèi)插序號為4、21兩顆衛(wèi)星坐標(biāo)，與11號衛(wèi)星坐標(biāo)內(nèi)插結(jié)果進(jìn)行比較，并利用Matlab軟件繪制出中誤差變化趨勢圖(如圖1所示)。圖1給出了衛(wèi)星序號為4、11、21這3顆衛(wèi)星5～18階插值的點位中誤差變化趨勢圖，其中圖1(a)為5～18階變化全圖，圖1(b)為局部放大圖。由圖1可以看出，用30 min內(nèi)插計算的3顆衛(wèi)星軌道的精度變化趨勢是一致的，隨著內(nèi)插階數(shù)的增加，內(nèi)插精度逐漸升高，精度最終穩(wěn)定在±2 mm左右。

表1 拉格朗日內(nèi)插結(jié)果對比表 mm

圖1 內(nèi)插精度變化趨勢圖

2.切比雪夫多項式多項式擬合效果分析

在切比雪夫多項式擬合中，輸入的擬合參數(shù)有兩個:擬合階數(shù)和擬合時間段長度，下面分別探討這兩個參數(shù)對擬合結(jié)果的影響。

取擬合時間段為4 h，選用序號為4的衛(wèi)星歷元值作為試驗數(shù)據(jù)，不同階數(shù)的切比雪夫多項式擬合結(jié)果如表2所示。

表2 切比雪夫多項式擬合結(jié)果對比表 mm

由表2可以得出在一定時間間隔內(nèi)，切比雪夫多項式的階數(shù)取得越高，擬合的數(shù)據(jù)精度越高，階數(shù)取13時，可以達(dá)到毫米級的精度。但是在時間段取為4 h時，該時間段內(nèi)的精密星歷歷元數(shù)為16，根據(jù)切比雪夫多項式擬合的數(shù)學(xué)原理可知，在該時間段內(nèi)階數(shù)最多只能取到13階，這就說明在短時間弧段內(nèi)，擬合階數(shù)受到歷元數(shù)的限制，難以通過增加擬合階數(shù)來提高擬合的精度。

取不同的時間段，并且將該時間段內(nèi)所有的歷元數(shù)據(jù)全部作為已知數(shù)據(jù)參與擬合，擬合的結(jié)果如表3所示。

表3 不同時間段擬合結(jié)果對比表 mm

由表3可以看出，所取的擬合時間間隔越長，擬合的精度越差，在擬合時間段超過5 h以后，擬合的精度顯著下降。原因是IGS精密星歷提供的歷元時間間隔為15 min，處于15 min之間的時刻沒有衛(wèi)星數(shù)據(jù)，而時間段越長所選用的離散數(shù)據(jù)點相關(guān)性越差，影響到切比雪夫多項式的擬合精度。

由于IGS精密星歷給出的最后一個歷元時刻為23時45分0秒，而在實際應(yīng)用中，有時需要使用23時45分0秒至24時0分0秒之間的衛(wèi)星坐標(biāo)數(shù)據(jù)，因此需要利用該時段之前的歷元數(shù)據(jù)來進(jìn)行外推。根據(jù)切比雪夫多項式的數(shù)學(xué)原理可知，切比雪夫多項式可以對擬合時段之外的數(shù)據(jù)進(jìn)行外推，下面探討對IGS精密星歷外推的效果。

取時間段為4 h，階數(shù)為12對第4顆衛(wèi)星的軌道進(jìn)行擬合，得到一個時間段內(nèi)16個歷元的殘差分布表(如表4所示)。

由表4可以看出，16號歷元的殘差明顯大于其他歷元，它的殘差為分米級，而其他歷元的殘差為厘米級。因為16號歷元在擬合的時間段外，屬于外推范圍。由于IGS精密星歷的歷元時間間隔較大，外推的衛(wèi)星軌道與實際軌道相差較遠(yuǎn)，因此殘差也比較大。所以在利用切比雪夫多項式擬合精密星歷時，如果對擬合數(shù)據(jù)精度要求不高，可以利用切比雪夫多項式進(jìn)行有限的外推;如對擬合數(shù)據(jù)精度要求較高，如精密單點定位，則需要將所需時刻納入到擬合的時間段之內(nèi)。

表4 前16個歷元的殘差分布表 mm

3.兩種數(shù)學(xué)模型處理效果比較

本文選用11號衛(wèi)星的歷元值作為試驗數(shù)據(jù)，并比較了兩種數(shù)學(xué)模型處理的效果。圖2給出了階數(shù)為16的拉格朗日內(nèi)插后的X坐標(biāo)分量的殘差分布圖;圖3給出了階數(shù)為13、時間段長度為4 h的切比雪夫多項式擬合后X坐標(biāo)分量的殘差分布圖。

由圖2可知，16階拉格朗日內(nèi)插后的X坐標(biāo)分量的最大誤差為2.8 mm左右;由圖3可知，13階切比雪夫多項式擬合的X坐標(biāo)分量最大誤差為1.6 mm左右，且拉格朗日內(nèi)插的出的數(shù)據(jù)比切比雪夫多項式擬合出的數(shù)據(jù)誤差波動的范圍要大一些。

圖2 16階拉格朗日內(nèi)插殘差分布圖

圖3 13階切比雪夫多項式擬合殘差分布圖

另外，通過對兩種數(shù)學(xué)模型表達(dá)式(1)和式(3)分析比較后得出，利用切比雪夫多項式擬合衛(wèi)星軌道易計算，占用的計算機內(nèi)存小，并且可外推出擬合時間段以外時刻的衛(wèi)星坐標(biāo)，其精度完全滿足導(dǎo)航衛(wèi)星軌道擬合的要求。因此，確定切比雪夫多項式擬合為最佳的衛(wèi)星軌道擬合方法，并選用階數(shù)為13，時間段長度為4 h的切比雪夫多項式進(jìn)行坐標(biāo)擬合。

四、虛擬衛(wèi)星軌道的生成

SP3精密星歷是地心地固坐標(biāo)系下給出的衛(wèi)星位置軌道數(shù)據(jù)，而在進(jìn)行虛擬仿真時，為了表達(dá)地球自轉(zhuǎn)，衛(wèi)星軌道數(shù)據(jù)應(yīng)轉(zhuǎn)化到天固坐標(biāo)系中。SP3精密星歷的歷元間隔為15 min，隨著歷元的遞增，地球轉(zhuǎn)過的角度也遞增，將地固坐標(biāo)系下坐標(biāo)轉(zhuǎn)化到天固坐標(biāo)系時，相應(yīng)的旋轉(zhuǎn)角度也遞增，遞增的步長大小就是15 min內(nèi)地球轉(zhuǎn)過的角度，因此通過將不同歷元的衛(wèi)星坐標(biāo)按照其對應(yīng)的自轉(zhuǎn)角度進(jìn)行轉(zhuǎn)換，就可得到天固坐標(biāo)系下衛(wèi)星的運動軌跡。為了驗證此方法的正確性，本文通過Matlab，對所采用的設(shè)計思路進(jìn)行了仿真驗證。在不考慮地固坐標(biāo)系在天固坐標(biāo)系中初始定向的情況下，以IGS發(fā)布的2008年5月1日0時0分0秒的SP3星歷為例進(jìn)行坐標(biāo)變換，變換前后衛(wèi)星軌跡如圖4所示。

其中，圖4(a)表示的是地固坐標(biāo)系下的軌道坐標(biāo)直接在天固坐標(biāo)系下形成的衛(wèi)星軌跡，畫出的是一個封閉的馬鞍線，圖4(b)表示的是經(jīng)過坐標(biāo)轉(zhuǎn)換后天固坐標(biāo)系下的衛(wèi)星軌跡;圖4(c)表示地固坐標(biāo)系下的軌跡在天固坐標(biāo)系下畫出的所有GPS衛(wèi)星軌跡;圖4(d)表示經(jīng)坐標(biāo)轉(zhuǎn)換后的天固坐標(biāo)系下的軌跡分布。可以看出，衛(wèi)星大致分布在6個軌道平面內(nèi)，符合GPS實際空間分布情況。

圖4 變換前后衛(wèi)星軌跡

根據(jù)切比雪夫多項式擬合原理和分析比較結(jié)果，由式(3)對IGS發(fā)布的2008年5月1日SP3星歷進(jìn)行計算，得到任意時刻的衛(wèi)星位置數(shù)據(jù)。然后，在基于MFC和WTK的單文檔GPS衛(wèi)星仿真系統(tǒng)中通過調(diào)用WTK的API函數(shù)動態(tài)渲染生成虛擬運行場景，實現(xiàn)效果如圖5所示。

圖5 GPS衛(wèi)星仿真系統(tǒng)運行界面

五、結(jié)束語

本文深入研究了利用GPS精密星歷擬合衛(wèi)星三維坐標(biāo)的方法。對拉格朗日多項式內(nèi)插和切比雪夫多項式擬合的處理效果進(jìn)行了分析和比較，給出了兩種擬合方法衛(wèi)星位置誤差與階數(shù)的關(guān)系，確定了切比雪夫多項式擬合為最佳的衛(wèi)星軌道擬合方法，并在虛擬環(huán)境中進(jìn)行了驗證與實現(xiàn)，提高了衛(wèi)星導(dǎo)航系統(tǒng)模擬仿真的精度，為我國北斗導(dǎo)航定位系統(tǒng)建設(shè)中衛(wèi)星軌道的設(shè)計明確了計算方法，為我國衛(wèi)星軌道虛擬仿真提供了參考。

［1］ 蔡艷輝，程鵬飛，李夕銀.衛(wèi)星坐標(biāo)的內(nèi)插和擬合［J］.全球定位系統(tǒng)，2003，27(3):10-13.

［2］ 張守建，李建成，邢樂林，等.兩種IGS精密星歷插值方法的比較分析［J］.大地測量與地球動力學(xué)，2007，33(2):80-83.

［3］ 王義峰，文鴻雁，劉立龍.用切比雪夫多項式標(biāo)準(zhǔn)化GPS衛(wèi)星軌道［J］.測繪與空間地理信息，2008，31(1):21-23.

［4］ 陳兆林，張書畢，佟瑞菊.用拉格朗日多項式內(nèi)插計算GPS衛(wèi)星位置［J］.全球定位系統(tǒng)，2007，33(2): 33-35.

［5］ FENG Yanming，ZHENG Yi，BAI Zhengdong.Representing GPS Orbits and Corrections Efficiently for Precise Wide Area Positioning［C］∥Proceedings of the 17th International Technical Meeting of the Satellite Division of the Institute of Navigation ION GNSS2004.California:［s.n.］，2004.

［6］ 李征航，黃勁松.GPS測量與數(shù)據(jù)處理［M］.武漢:武漢大學(xué)出版社，2005:33-37.

［7］ 孔巧麗.用切貝雪夫多項式擬合 GPS衛(wèi)星精密坐標(biāo)［J］.測繪通報，2006(8):1-3.

［8］ 武漢大學(xué)測繪學(xué)院測量平差學(xué)科組.誤差理論與測量平差基礎(chǔ)［M］.武漢:武漢大學(xué)出版社，2003:102-104，118.

Research on the Orbit Fitting Methods of Navigation Satellite and Simulation Realization

WAN Jiahuan，ZHUANG Chunhua，CHEN Xiuwan，LI Zhihui，ZHANG Weiyi

0494-0911(2012)07-0001-05

P228.4

B

2011-07-07

國家國防科技工業(yè)局民用航天預(yù)研項目

萬家歡(1986—)，男，湖北浠水人，碩士，主要從事衛(wèi)星導(dǎo)航應(yīng)用方面的研究工作。