切比雪夫多項式在GLONASS廣播星歷中的應用

寇瑞雄，楊樹文，化希瑞

（蘭州交通大學 測繪與地理信息學院/地理國情監測技術應用國家地方聯合工程研究中心/甘肅省地理國情監測工程實驗室，蘭州 730070）

0 引言

格洛納斯衛星導航系統（global navigation satellite system, GLONASS）于1996年初衛星組網成功并正式投入運行，2003年，俄羅斯政府啟動GLONASS現代化工作，2010年重新建成了24顆衛星組成的星座，2015年發射新衛星并改善了各系統的性能[1]。俄羅斯航天局計劃在2019—2033年間，發射46顆不同型號的導航衛星[2]，并于2019年5月27日成功發射了一顆“GLONASS-M”導航衛星[3]，因此，需要持續關注GLONASS的發展狀況。實時導航定位需要用廣播星歷計算衛星位置，但GLONASS計算衛星位置時，需要以參考時刻為中心向前、向后積分并不斷迭代[4]，這樣會消耗較長時間和占用較多內存空間，影響計算效率。為此，可用1個時間多項式來表示衛星星歷，計算衛星位置時只需調用多項式系數，這樣能提高運算效率[1]。

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

1 GLONASS廣播星歷

1.1 GLONASS衛星位置計算

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

GLONASS衛星廣播星歷采用PZ90坐標系，通常以衛星參數為軌道初值，根據衛星運動攝動力模型，使用數值積分法計算衛星坐標。文獻[14]給出詳細計算過程，其中加速度的簡略公式為

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

將式（1）衛星加速度表達成關于坐標、速度和日月攝動加速度的函數形式，可得ti時刻的加速度函數式為

式中：(u i,vi,wi)為ti時刻(xi,yi,zi)方向的速度；（axi,ayi,azi）為ti時刻時3個方向上的加速度；（ax,ay,az）為日月攝動加速度。

以參考時刻t0為初始狀態，則tb時刻衛星位置的積分方程表示為：

式中：s(t0)為t0時刻的衛星位置；a為衛星加速度；v0為t0時刻的衛星速度。使用定步長4階朗格-庫特（Runge-Kutte）法對衛星軌道進行積分，經過二次積分，會得到tb時刻衛星位置。文獻[4]對使用定步長4階Runge-Kutte法的積分過程有詳細介紹。

1.2 GLONASS廣播星歷精度評定

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

混合精密星歷一般提供24顆GLONASS衛星的坐標，但該天R04、R06和R11衛星缺失數據，因此用剩下的21顆精密星歷坐標為真值，驗證廣播星歷計算衛星坐標的精度。統計得到當天有數據衛星的廣播星歷相對于精密星歷的誤差絕對值的最大值、絕對值的平均誤差和均方差，如圖1、圖2和圖3所示。

圖1 GLONASS廣播星歷最大誤差

圖2 GLONASS廣播星歷平均誤差

從圖1至圖3可得，該天GLONASS廣播星歷的最大誤差出現在R07星的X方向上，為6.38 m，大多數衛星在X、Y和Z方向上的最大誤差不超過5 m，最大平均誤差出現在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廣播星歷給出的大多數衛星在每個方向上的均方差都小于3 m。

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

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

因而衛星的3維坐標可表示為

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

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

式中m為X坐標個數。從而可以計算X坐標的切比雪夫多項式系數C為

同理，使用上述方法可以求得Y和Z坐標的切比雪夫多項式系數，根據計算得到的多項式可以計算擬合區間內任意時刻的3維坐標。

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

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

文獻[14]提出以參考歷元前后15 min為有效的擬合區間，使用第2節的方法計算出所有參考歷元前后15 min時間段內，以30 s為時間間隔的衛星3維坐標，作為后續切比雪夫多項式擬合結果檢核的真值。將擬合出的衛星坐標與對應時刻原有方法求得的衛星坐標求差，對殘差進行統計，求得絕對值的最大值、平均值和均方差，進行擬合精度評定。由于切比雪夫多項式擬合精度與時間間隔、節點個數和擬合階數都有關系，需要分兩個方面進行討論：一方面是在固定時間間隔下，不同節點個數和擬合階數對精度的影響；另一方面是在不同時間間隔下，達到理想擬合精度時，節點個數和擬合階數最優組合的變化情況。

3.1 不同節點數和擬合階數下的擬合精度

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

表1 X方向上節點個數與擬合階數不同組合的擬合誤差

（續表2）

表3 Z方向上節點個數與擬合階數不同組合的擬合誤差

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

3.2 不同時間間隔下的最優擬合階數

考慮到在固定時間間隔下，擬合誤差不會隨擬合階數的增大而一直減小，因此，將在固定擬合區間里，不同時間間隔下，X、Y和Z三個方向的擬合精度指標均達到亞毫米級且擬合階數最小的擬合階數定義為該時間間隔下的最優擬合階數。繼續使用10 d的R01衛星的廣播星歷，在每個參考歷元前后15 min區間里，將時間間隔分別設為120、150、180和210 s，統計得到最優組合擬合階數，擬合精度分析如表4所示。同時給出R01衛星在時間間隔為210 s擬合最優組合時的X、Y和Z三個方向是480個參考歷元區間內的所有檢核點擬合誤差，如圖4至圖6所示。

表4 不同時間間隔下的節點個數和擬合階數的最佳組合

圖4 210 s時間間隔最佳擬合階數下X方向的擬合誤差

圖5 210 s時間間隔最佳擬合階數下Y方向的擬合誤差

圖6 210 s時間間隔最佳擬合階數下Z方向的擬合誤差

分析表4可得：在固定擬合區間內，不同時間間隔下，節點個數不同，擬合階數取9階的情況下，最大誤差、平均誤差和均方差都小于1 mm，滿足精度需求；在相同擬合條件下，Z方向的擬合精度高于X和Y方向的擬合精度。同樣從圖4至圖6明顯看出，Z方向的擬合誤差小于X和Y方向的擬合誤差，但3個方向的擬合誤差均小于1 mm，所以對廣播星歷的精度影響可以忽略。在實際應用時，根據精度需求和運算效率，選擇上述不同的切比雪夫多項式擬合組合，可以很好地擬合GLONASS廣播星歷。

4 結束語

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

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

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

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

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