切比雪夫多項式在GLONASS廣播星歷中的應(yīng)用

寇瑞雄，楊樹文，化希瑞

（蘭州交通大學(xué) 測繪與地理信息學(xué)院/地理國情監(jiān)測技術(shù)應(yīng)用國家地方聯(lián)合工程研究中心/甘肅省地理國情監(jiān)測工程實驗室，蘭州 730070）

0 引言

格洛納斯衛(wèi)星導(dǎo)航系統(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顆不同型號的導(dǎo)航衛(wèi)星[2]，并于2019年5月27日成功發(fā)射了一顆“GLONASS-M”導(dǎo)航衛(wèi)星[3]，因此，需要持續(xù)關(guān)注GLONASS的發(fā)展?fàn)顩r。實時導(dǎo)航定位需要用廣播星歷計算衛(wèi)星位置，但GLONASS計算衛(wèi)星位置時，需要以參考時刻為中心向前、向后積分并不斷迭代[4]，這樣會消耗較長時間和占用較多內(nèi)存空間，影響計算效率。為此，可用1個時間多項式來表示衛(wèi)星星歷，計算衛(wèi)星位置時只需調(diào)用多項式系數(shù)，這樣能提高運算效率[1]。

當(dāng)前研究全球定位系統(tǒng)（global positioning system, GPS）和北斗衛(wèi)星導(dǎo)航系統(tǒng)（BeiDou navigation satellite system, BDS）廣播星歷擬合的文獻(xiàn)較多[5-9]，且BDS和GPS使用廣播星歷計算衛(wèi)星坐標(biāo)都采用16參數(shù)模型[10]，而GLONASS廣播星歷計算衛(wèi)星坐標(biāo)是通過二次積分且不斷迭代的方法，與GPS和BDS所用方法完全不同，導(dǎo)致GLONASS廣播星歷擬合存在差異。文獻(xiàn)[11-12]只對GLONASS衛(wèi)星廣播星歷的精度進(jìn)行了分析；文獻(xiàn)[13]分析了GLONASS廣播星歷和GPS廣播星歷的精度，使用拉格朗日插值法對GLONASS廣播星歷進(jìn)行插值分析。但拉格朗日插值法在對離散數(shù)據(jù)進(jìn)行插值時，隨著階數(shù)的增高會出現(xiàn)明顯的“龍格”現(xiàn)象，影響插值精度。本文研究切比雪夫多項式擬合法對GLONASS衛(wèi)星坐標(biāo)擬合精度的影響。

1 GLONASS廣播星歷

1.1 GLONASS衛(wèi)星位置計算

GLONASS廣播星歷與GPS廣播星歷提供的參數(shù)完全不一樣：GPS給出衛(wèi)星的開普勒軌道數(shù)據(jù)和衛(wèi)星時鐘，每2 h廣播1次；而GLONASS廣播星歷提供在俄羅斯大地坐標(biāo)框架（parametry zelmy1990, PZ90）下參考時刻的衛(wèi)星位置（X,Y,Z）、衛(wèi)星3個方向的速度和3個方向的日月攝動加速度，每30 min廣播1次[1]。

GLONASS衛(wèi)星廣播星歷采用PZ90坐標(biāo)系，通常以衛(wèi)星參數(shù)為軌道初值，根據(jù)衛(wèi)星運動攝動力模型，使用數(shù)值積分法計算衛(wèi)星坐標(biāo)。文獻(xiàn)[14]給出詳細(xì)計算過程，其中加速度的簡略公式為

式中：GM為地球引力常量（GM =398 600.44 km3/s2）；為衛(wèi)星3維坐標(biāo)；為衛(wèi)星3個方向的速度；r為衛(wèi)星到地球質(zhì)心的距離，為地球的長半徑，為地球重力系數(shù)，為地球自轉(zhuǎn)加速度，ω=0.000 072 921 15 rad/s；（ax,ay,az）為日月攝動加速度。

將式（1）衛(wèi)星加速度表達(dá)成關(guān)于坐標(biāo)、速度和日月攝動加速度的函數(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ìn)行積分，經(jīng)過二次積分，會得到tb時刻衛(wèi)星位置。文獻(xiàn)[4]對使用定步長4階Runge-Kutte法的積分過程有詳細(xì)介紹。

1.2 GLONASS廣播星歷精度評定

廣播星歷對實時導(dǎo)航定位精度有著重要的作用，且GLONASS各系統(tǒng)一直處于不斷完善的狀態(tài)，所以評定GLONASS廣播星歷的精度有著重要意義。本文使用德國地學(xué)研究中心（Deutsches Geo Forschungs Zentrum, GFZ）提供的2019-08-29混合廣播星歷和混合精密星歷。用廣播星歷計算精密星歷歷元時刻的衛(wèi)星位置，由于混合星歷的時間系統(tǒng)不是GLONASS時，所以計算時，需要注意時間系統(tǒng)的統(tǒng)一[15]。

混合精密星歷一般提供24顆GLONASS衛(wèi)星的坐標(biāo)，但該天R04、R06和R11衛(wèi)星缺失數(shù)據(jù)，因此用剩下的21顆精密星歷坐標(biāo)為真值，驗證廣播星歷計算衛(wèi)星坐標(biāo)的精度。統(tǒng)計得到當(dā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。

2 切比雪夫多項式擬合原理

對GLONASS廣播星歷，在[t0,t0+Δt]時間間隔內(nèi)，使用n階切比雪夫多項式進(jìn)行擬合，其中t0為初始?xì)v元，Δt為擬合區(qū)間長度。首先將t∈ [t0,t0+Δt]處理成τ∈[ -1 , 1]的形式，即

因而衛(wèi)星的3維坐標(biāo)可表示為

式中：n為切比雪夫多項式階數(shù)；CXi、CYi和CZi分別為3個坐標(biāo)分量的切比雪夫多項式系數(shù)；切比雪夫多項式Ti(τ)用下式遞推得到，即

以X坐標(biāo)為例，使用間接平差的方法求解多項式系數(shù)，將式（6）的第1個式子列誤差方程并整理成矩陣形式為：

式中m為X坐標(biāo)個數(shù)。從而可以計算X坐標(biāo)的切比雪夫多項式系數(shù)C為

同理，使用上述方法可以求得Y和Z坐標(biāo)的切比雪夫多項式系數(shù)，根據(jù)計算得到的多項式可以計算擬合區(qū)間內(nèi)任意時刻的3維坐標(biāo)。

3 GLONASS廣播星歷擬合精度分析

本文使用2019-08-28—2019-09-06的總共10 d的GLONASS廣播星歷數(shù)據(jù)，由于GLONASS衛(wèi)星軌道都是圓形軌道，所以選取R01星的廣播星歷進(jìn)行實驗分析。GLONASS廣播星歷每30 min發(fā)布1次，每顆衛(wèi)星1 d內(nèi)共有48組軌道參數(shù)，R01星10 d共有480組軌道參數(shù)。

文獻(xiàn)[14]提出以參考?xì)v元前后15 min為有效的擬合區(qū)間，使用第2節(jié)的方法計算出所有參考?xì)v元前后15 min時間段內(nèi)，以30 s為時間間隔的衛(wèi)星3維坐標(biāo)，作為后續(xù)切比雪夫多項式擬合結(jié)果檢核的真值。將擬合出的衛(wèi)星坐標(biāo)與對應(yīng)時刻原有方法求得的衛(wèi)星坐標(biāo)求差，對殘差進(jìn)行統(tǒng)計，求得絕對值的最大值、平均值和均方差，進(jìn)行擬合精度評定。由于切比雪夫多項式擬合精度與時間間隔、節(jié)點個數(shù)和擬合階數(shù)都有關(guān)系，需要分兩個方面進(jìn)行討論：一方面是在固定時間間隔下，不同節(jié)點個數(shù)和擬合階數(shù)對精度的影響；另一方面是在不同時間間隔下，達(dá)到理想擬合精度時，節(jié)點個數(shù)和擬合階數(shù)最優(yōu)組合的變化情況。

3.1 不同節(jié)點數(shù)和擬合階數(shù)下的擬合精度

在1個參考?xì)v元前后15 min的時間段內(nèi)，設(shè)節(jié)點時間間隔為120 s，則30 min共有16個已知點；每30 s設(shè)置1個檢核點，則總共有61個檢核點，通過選取不同節(jié)點個數(shù)和擬合階數(shù)對廣播星歷進(jìn)行擬合，并統(tǒng)計在X、Y和Z方向上的擬合誤差。本節(jié)總共使用480個參考?xì)v元擬合區(qū)間，檢核點最多時達(dá)到29 280個，統(tǒng)計在不同節(jié)點個數(shù)和擬合階數(shù)下廣播星歷3個方向的擬合誤差，具體如表1至表3所示，需要注意的是，節(jié)點個數(shù)要大于等于擬合階數(shù)，所以表1至表3給出了表格的上三角部分。用誤差絕對值的最大值（Max）、誤差絕對值的平均值（Mean）和誤差的均方差（Std）作為擬合精度的評價指標(biāo)。

表1 X方向上節(jié)點個數(shù)與擬合階數(shù)不同組合的擬合誤差

（續(xù)表2）

表3 Z方向上節(jié)點個數(shù)與擬合階數(shù)不同組合的擬合誤差

從表1至表3可知，當(dāng)節(jié)點個數(shù)和擬合階數(shù)都為7時，X、Y和Z三個方向的誤差絕對值的最大值分別為1.613、1.583和0.359 mm；當(dāng)節(jié)點個數(shù)為10，擬合階數(shù)為8時，X、Y和Z三個方向的擬合精度指標(biāo)均小于1 mm；當(dāng)節(jié)點個數(shù)為13，擬合階數(shù)為8時，Z方向的誤差絕對值的最大值為0.246 mm、平均值為0.046 mm、均方差為0.04 mm，X和Y方向誤差絕對值的最大值小于2 mm。當(dāng)節(jié)點個數(shù)為16時，整體上隨著擬合階數(shù)的增加擬合精度迅速提高，擬合階數(shù)為9時，3個方向誤差絕對值的最大值均為小于1 mm；當(dāng)擬合階數(shù)取10～14的情況下，3個方向的誤差均在逐漸變大，3個方向的誤差絕對值的最大值為2.176 mm，遠(yuǎn)小于GLONASS廣播星歷米級的誤差；但當(dāng)擬合階數(shù)取15和16時，3個方向的擬合誤差較大且需要的消耗更多時間，因此在廣播星歷擬合過程要考慮擬合誤差和計算效率的問題。

3.2 不同時間間隔下的最優(yōu)擬合階數(shù)

考慮到在固定時間間隔下，擬合誤差不會隨擬合階數(shù)的增大而一直減小，因此，將在固定擬合區(qū)間里，不同時間間隔下，X、Y和Z三個方向的擬合精度指標(biāo)均達(dá)到亞毫米級且擬合階數(shù)最小的擬合階數(shù)定義為該時間間隔下的最優(yōu)擬合階數(shù)。繼續(xù)使用10 d的R01衛(wèi)星的廣播星歷，在每個參考?xì)v元前后15 min區(qū)間里，將時間間隔分別設(shè)為120、150、180和210 s，統(tǒng)計得到最優(yōu)組合擬合階數(shù)，擬合精度分析如表4所示。同時給出R01衛(wèi)星在時間間隔為210 s擬合最優(yōu)組合時的X、Y和Z三個方向是480個參考?xì)v元區(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，所以對廣播星歷的精度影響可以忽略。在實際應(yīng)用時，根據(jù)精度需求和運算效率，選擇上述不同的切比雪夫多項式擬合組合，可以很好地擬合GLONASS廣播星歷。

4 結(jié)束語

本文首先使用2019-08-28—2019-09-06共10 d的GLONASS廣播星歷數(shù)據(jù)，計算了衛(wèi)星坐標(biāo)，然后選取2019-8-29的衛(wèi)星坐標(biāo)數(shù)據(jù)，評定了廣播星歷的精度，最后利用切比雪夫多項式擬合法對10 d的衛(wèi)星坐標(biāo)進(jìn)行擬合，得出以下結(jié)論：

1）GLONASS廣播星歷計算出的坐標(biāo)衛(wèi)星（X、Y、Z），在3個方向上的均方差絕大多數(shù)不超過3 m，最大誤差小于6.5 m，誤差平均值基本小于2.5 m，總之廣播星歷的誤差處于米級。

2）當(dāng)節(jié)點時間間隔固定時，不同節(jié)點個數(shù)的最優(yōu)擬合階數(shù)不完全一樣。隨著節(jié)點個數(shù)的增加，擬合階數(shù)過低或者過高的擬合精度，均低于最優(yōu)擬合階數(shù)的精度，在時間間隔為120 s，節(jié)點個數(shù)為16，選擇擬合階數(shù)為9時，擬合誤差最大值為0.16 mm，遠(yuǎn)小于廣播星歷米級誤差，可以忽略不計。

3）在固定擬合區(qū)間的情況下，不同時間間隔的最優(yōu)擬合階數(shù)變化不大，當(dāng)擬合階數(shù)為9階時，擬合精度可以達(dá)到亞毫米級，滿足精度要求。盡管擬合時間間隔取值較密時，精度相對較高，但效率會有所下降。因此在實際使用過程中，須兼顧精度和效率，選擇合適的擬合時間間隔和階數(shù)。

4）切比雪夫多項式擬合法適用于GLONASS廣播星歷的衛(wèi)星坐標(biāo)擬合，在選取合適的擬合時間間隔和擬合階數(shù)時，擬合精度足以滿足要求。故切比雪夫多項式可以作為1種新的廣播星歷計算衛(wèi)星坐標(biāo)的形式。