寇瑞雄,楊樹文,化希瑞
(蘭州交通大學 測繪與地理信息學院/地理國情監(jiān)測技術應用國家地方聯(lián)合工程研究中心/甘肅省地理國情監(jiān)測工程實驗室,蘭州 730070)
格洛納斯衛(wèi)星導航系統(tǒng)(global navigation satellite system, GLONASS)于1996年初衛(wèi)星組網(wǎng)成功并正式投入運行,2003年,俄羅斯政府啟動GLONASS現(xiàn)代化工作,2010年重新建成了24顆衛(wèi)星組成的星座,2015年發(fā)射新衛(wèi)星并改善了各系統(tǒng)的性能[1]。俄羅斯航天局計劃在2019—2033年間,發(fā)射46顆不同型號的導航衛(wèi)星[2],并于2019年5月27日成功發(fā)射了一顆“GLONASS-M”導航衛(wèi)星[3],因此,需要持續(xù)關注GLONASS的發(fā)展狀況。實時導航定位需要用廣播星歷計算衛(wèi)星位置,但GLONASS計算衛(wèi)星位置時,需要以參考時刻為中心向前、向后積分并不斷迭代[4],這樣會消耗較長時間和占用較多內(nèi)存空間,影響計算效率。為此,可用1個時間多項式來表示衛(wèi)星星歷,計算衛(wèi)星位置時只需調(diào)用多項式系數(shù),這樣能提高運算效率[1]。
當前研究全球定位系統(tǒng)(global positioning system, GPS)和北斗衛(wèi)星導航系統(tǒng)(BeiDou navigation satellite system, BDS)廣播星歷擬合的文獻較多[5-9],且BDS和GPS使用廣播星歷計算衛(wèi)星坐標都采用16參數(shù)模型[10],而GLONASS廣播星歷計算衛(wèi)星坐標是通過二次積分且不斷迭代的方法,與GPS和BDS所用方法完全不同,導致GLONASS廣播星歷擬合存在差異。文獻[11-12]只對GLONASS衛(wèi)星廣播星歷的精度進行了分析;文獻[13]分析了GLONASS廣播星歷和GPS廣播星歷的精度,使用拉格朗日插值法對GLONASS廣播星歷進行插值分析。但拉格朗日插值法在對離散數(shù)據(jù)進行插值時,隨著階數(shù)的增高會出現(xiàn)明顯的“龍格”現(xiàn)象,影響插值精度。本文研究切比雪夫多項式擬合法對GLONASS衛(wèi)星坐標擬合精度的影響。
GLONASS廣播星歷與GPS廣播星歷提供的參數(shù)完全不一樣:GPS給出衛(wèi)星的開普勒軌道數(shù)據(jù)和衛(wèi)星時鐘,每2 h廣播1次;而GLONASS廣播星歷提供在俄羅斯大地坐標框架(parametry zelmy1990, PZ90)下參考時刻的衛(wèi)星位置(X,Y,Z)、衛(wèi)星3個方向的速度和3個方向的日月攝動加速度,每30 min廣播1次[1]。
GLONASS衛(wèi)星廣播星歷采用PZ90坐標系,通常以衛(wèi)星參數(shù)為軌道初值,根據(jù)衛(wèi)星運動攝動力模型,使用數(shù)值積分法計算衛(wèi)星坐標。文獻[14]給出詳細計算過程,其中加速度的簡略公式為
式中:GM為地球引力常量(GM =398 600.44 km3/s2);為衛(wèi)星3維坐標;為衛(wèi)星3個方向的速度;r為衛(wèi)星到地球質(zhì)心的距離,為地球的長半徑,為地球重力系數(shù),為地球自轉(zhuǎn)加速度,ω=0.000 072 921 15 rad/s;(ax,ay,az)為日月攝動加速度。
將式(1)衛(wèi)星加速度表達成關于坐標、速度和日月攝動加速度的函數(shù)形式,可得ti時刻的加速度函數(shù)式為
式中:(u i,vi,wi)為ti時刻(xi,yi,zi)方向的速度;(axi,ayi,azi)為ti時刻時3個方向上的加速度;(ax,ay,az)為日月攝動加速度。
以參考時刻t0為初始狀態(tài),則tb時刻衛(wèi)星位置的積分方程表示為:
式中:s(t0)為t0時刻的衛(wèi)星位置;a為衛(wèi)星加速度;v0為t0時刻的衛(wèi)星速度。使用定步長4階朗格-庫特(Runge-Kutte)法對衛(wèi)星軌道進行積分,經(jīng)過二次積分,會得到tb時刻衛(wèi)星位置。文獻[4]對使用定步長4階Runge-Kutte法的積分過程有詳細介紹。
廣播星歷對實時導航定位精度有著重要的作用,且GLONASS各系統(tǒng)一直處于不斷完善的狀態(tài),所以評定GLONASS廣播星歷的精度有著重要意義。本文使用德國地學研究中心(Deutsches Geo Forschungs Zentrum, GFZ)提供的2019-08-29混合廣播星歷和混合精密星歷。用廣播星歷計算精密星歷歷元時刻的衛(wèi)星位置,由于混合星歷的時間系統(tǒng)不是GLONASS時,所以計算時,需要注意時間系統(tǒng)的統(tǒng)一[15]。
混合精密星歷一般提供24顆GLONASS衛(wèi)星的坐標,但該天R04、R06和R11衛(wèi)星缺失數(shù)據(jù),因此用剩下的21顆精密星歷坐標為真值,驗證廣播星歷計算衛(wèi)星坐標的精度。統(tǒng)計得到當天有數(shù)據(jù)衛(wèi)星的廣播星歷相對于精密星歷的誤差絕對值的最大值、絕對值的平均誤差和均方差,如圖1、圖2和圖3所示。
圖1 GLONASS廣播星歷最大誤差
圖2 GLONASS廣播星歷平均誤差
從圖1至圖3可得,該天GLONASS廣播星歷的最大誤差出現(xiàn)在R07星的X方向上,為6.38 m,大多數(shù)衛(wèi)星在X、Y和Z方向上的最大誤差不超過5 m,最大平均誤差出現(xiàn)在R13星的Z方向上,為3.24 m,只有R01星、R05星、R07星和R13星在部分方向上的平均誤差大于2.5 m,R11星在Y方向上的平均誤差最小,為0.72 m。R01星、R02星、R03星、R05星、R07星和R08星在每個方向上的均方差為2.5~3.5 m之間;R13星在Z方向上的均方差為3.6 m,是該天均方差的最大值,除了R13星之外,在R11星至R24星的均方差基本都小于2 m。從整體上看,GLONASS廣播星歷給出的大多數(shù)衛(wèi)星在每個方向上的均方差都小于3 m。
對GLONASS廣播星歷,在[t0,t0+Δt]時間間隔內(nèi),使用n階切比雪夫多項式進行擬合,其中t0為初始歷元,Δt為擬合區(qū)間長度。首先將t∈ [t0,t0+Δt]處理成τ∈[ -1 , 1]的形式,即
因而衛(wèi)星的3維坐標可表示為
式中:n為切比雪夫多項式階數(shù);CXi、CYi和CZi分別為3個坐標分量的切比雪夫多項式系數(shù);切比雪夫多項式Ti(τ)用下式遞推得到,即
以X坐標為例,使用間接平差的方法求解多項式系數(shù),將式(6)的第1個式子列誤差方程并整理成矩陣形式為:
式中m為X坐標個數(shù)。從而可以計算X坐標的切比雪夫多項式系數(shù)C為
同理,使用上述方法可以求得Y和Z坐標的切比雪夫多項式系數(shù),根據(jù)計算得到的多項式可以計算擬合區(qū)間內(nèi)任意時刻的3維坐標。
本文使用2019-08-28—2019-09-06的總共10 d的GLONASS廣播星歷數(shù)據(jù),由于GLONASS衛(wèi)星軌道都是圓形軌道,所以選取R01星的廣播星歷進行實驗分析。GLONASS廣播星歷每30 min發(fā)布1次,每顆衛(wèi)星1 d內(nèi)共有48組軌道參數(shù),R01星10 d共有480組軌道參數(shù)。
文獻[14]提出以參考歷元前后15 min為有效的擬合區(qū)間,使用第2節(jié)的方法計算出所有參考歷元前后15 min時間段內(nèi),以30 s為時間間隔的衛(wèi)星3維坐標,作為后續(xù)切比雪夫多項式擬合結(jié)果檢核的真值。將擬合出的衛(wèi)星坐標與對應時刻原有方法求得的衛(wèi)星坐標求差,對殘差進行統(tǒng)計,求得絕對值的最大值、平均值和均方差,進行擬合精度評定。由于切比雪夫多項式擬合精度與時間間隔、節(jié)點個數(shù)和擬合階數(shù)都有關系,需要分兩個方面進行討論:一方面是在固定時間間隔下,不同節(jié)點個數(shù)和擬合階數(shù)對精度的影響;另一方面是在不同時間間隔下,達到理想擬合精度時,節(jié)點個數(shù)和擬合階數(shù)最優(yōu)組合的變化情況。
在1個參考歷元前后15 min的時間段內(nèi),設節(jié)點時間間隔為120 s,則30 min共有16個已知點;每30 s設置1個檢核點,則總共有61個檢核點,通過選取不同節(jié)點個數(shù)和擬合階數(shù)對廣播星歷進行擬合,并統(tǒng)計在X、Y和Z方向上的擬合誤差。本節(jié)總共使用480個參考歷元擬合區(qū)間,檢核點最多時達到29 280個,統(tǒng)計在不同節(jié)點個數(shù)和擬合階數(shù)下廣播星歷3個方向的擬合誤差,具體如表1至表3所示,需要注意的是,節(jié)點個數(shù)要大于等于擬合階數(shù),所以表1至表3給出了表格的上三角部分。用誤差絕對值的最大值(Max)、誤差絕對值的平均值(Mean)和誤差的均方差(Std)作為擬合精度的評價指標。
表1 X方向上節(jié)點個數(shù)與擬合階數(shù)不同組合的擬合誤差
(續(xù)表2)
表3 Z方向上節(jié)點個數(shù)與擬合階數(shù)不同組合的擬合誤差
從表1至表3可知,當節(jié)點個數(shù)和擬合階數(shù)都為7時,X、Y和Z三個方向的誤差絕對值的最大值分別為1.613、1.583和0.359 mm;當節(jié)點個數(shù)為10,擬合階數(shù)為8時,X、Y和Z三個方向的擬合精度指標均小于1 mm;當節(jié)點個數(shù)為13,擬合階數(shù)為8時,Z方向的誤差絕對值的最大值為0.246 mm、平均值為0.046 mm、均方差為0.04 mm,X和Y方向誤差絕對值的最大值小于2 mm。當節(jié)點個數(shù)為16時,整體上隨著擬合階數(shù)的增加擬合精度迅速提高,擬合階數(shù)為9時,3個方向誤差絕對值的最大值均為小于1 mm;當擬合階數(shù)取10~14的情況下,3個方向的誤差均在逐漸變大,3個方向的誤差絕對值的最大值為2.176 mm,遠小于GLONASS廣播星歷米級的誤差;但當擬合階數(shù)取15和16時,3個方向的擬合誤差較大且需要的消耗更多時間,因此在廣播星歷擬合過程要考慮擬合誤差和計算效率的問題。
考慮到在固定時間間隔下,擬合誤差不會隨擬合階數(shù)的增大而一直減小,因此,將在固定擬合區(qū)間里,不同時間間隔下,X、Y和Z三個方向的擬合精度指標均達到亞毫米級且擬合階數(shù)最小的擬合階數(shù)定義為該時間間隔下的最優(yōu)擬合階數(shù)。繼續(xù)使用10 d的R01衛(wèi)星的廣播星歷,在每個參考歷元前后15 min區(qū)間里,將時間間隔分別設為120、150、180和210 s,統(tǒng)計得到最優(yōu)組合擬合階數(shù),擬合精度分析如表4所示。同時給出R01衛(wèi)星在時間間隔為210 s擬合最優(yōu)組合時的X、Y和Z三個方向是480個參考歷元區(qū)間內(nèi)的所有檢核點擬合誤差,如圖4至圖6所示。
表4 不同時間間隔下的節(jié)點個數(shù)和擬合階數(shù)的最佳組合
圖4 210 s時間間隔最佳擬合階數(shù)下X方向的擬合誤差
圖5 210 s時間間隔最佳擬合階數(shù)下Y方向的擬合誤差
圖6 210 s時間間隔最佳擬合階數(shù)下Z方向的擬合誤差
分析表4可得:在固定擬合區(qū)間內(nèi),不同時間間隔下,節(jié)點個數(shù)不同,擬合階數(shù)取9階的情況下,最大誤差、平均誤差和均方差都小于1 mm,滿足精度需求;在相同擬合條件下,Z方向的擬合精度高于X和Y方向的擬合精度。同樣從圖4至圖6明顯看出,Z方向的擬合誤差小于X和Y方向的擬合誤差,但3個方向的擬合誤差均小于1 mm,所以對廣播星歷的精度影響可以忽略。在實際應用時,根據(jù)精度需求和運算效率,選擇上述不同的切比雪夫多項式擬合組合,可以很好地擬合GLONASS廣播星歷。
本文首先使用2019-08-28—2019-09-06共10 d的GLONASS廣播星歷數(shù)據(jù),計算了衛(wèi)星坐標,然后選取2019-8-29的衛(wèi)星坐標數(shù)據(jù),評定了廣播星歷的精度,最后利用切比雪夫多項式擬合法對10 d的衛(wèi)星坐標進行擬合,得出以下結(jié)論:
1)GLONASS廣播星歷計算出的坐標衛(wèi)星(X、Y、Z),在3個方向上的均方差絕大多數(shù)不超過3 m,最大誤差小于6.5 m,誤差平均值基本小于2.5 m,總之廣播星歷的誤差處于米級。
2)當節(jié)點時間間隔固定時,不同節(jié)點個數(shù)的最優(yōu)擬合階數(shù)不完全一樣。隨著節(jié)點個數(shù)的增加,擬合階數(shù)過低或者過高的擬合精度,均低于最優(yōu)擬合階數(shù)的精度,在時間間隔為120 s,節(jié)點個數(shù)為16,選擇擬合階數(shù)為9時,擬合誤差最大值為0.16 mm,遠小于廣播星歷米級誤差,可以忽略不計。
3)在固定擬合區(qū)間的情況下,不同時間間隔的最優(yōu)擬合階數(shù)變化不大,當擬合階數(shù)為9階時,擬合精度可以達到亞毫米級,滿足精度要求。盡管擬合時間間隔取值較密時,精度相對較高,但效率會有所下降。因此在實際使用過程中,須兼顧精度和效率,選擇合適的擬合時間間隔和階數(shù)。
4)切比雪夫多項式擬合法適用于GLONASS廣播星歷的衛(wèi)星坐標擬合,在選取合適的擬合時間間隔和擬合階數(shù)時,擬合精度足以滿足要求。故切比雪夫多項式可以作為1種新的廣播星歷計算衛(wèi)星坐標的形式。