不同軌道類型LEO衛星軌道擬合及預報精度研究

譚理慶，彭琦，曹陽，楊鑫，唐帥，劉俊

( 重慶兩江衛星移動通信有限公司, 重慶 401120 )

0 引 言

近年來，世界許多大國及公司陸續提出并開始建設服務于衛星互聯網、物聯網的低軌衛星星座，其中部分星座還具備導航增強服務功能. 未來，低軌衛星系統不僅可以促進衛星導航增強系統向星/地基增強一體化方向快速發展[1-2]，更能助力構建集通信、導航、遙感于一體的天基信息實時服務系統[3]，而高精度、高可靠的低軌衛星軌道是實現上述功能的前提保障條件之一. 高精度最終/預報精密星歷通常都只給出具有一定時間間隔的衛星位置，因此必須通過擬合插值獲得需要時刻的衛星位置. 同時為了保障星座的安全運行，當低軌道地球(LEO)衛星無法獲得實時定軌結果時，還需要利用前面歷元的定軌結果進行一定時間的外推預報.

目前通常以一定時間間隔給出衛星軌道根數或空間位置，可通過擬合/插值獲得任意所需時刻的衛星位置；對于LEO衛星軌道擬合/插值的研究大多基于切比雪夫多項式、拉格朗日多項式、牛頓多項式[4-6]，上述三種方法在實際應用中插值點在插值弧段中間部分可以取得較高精度，但在靠近插值弧段兩端部分會隨著多項式階數的增加而出現龍格現象. 此外，最佳平方逼近多項式[7]、最小二乘曲線擬合[8]等方法也被用于LEO衛星軌道擬合插值，均取得了厘米級精度. 在LEO衛星軌道短期預報方面，基于動力學的方法可以實現對低軌衛星長時間、高精度的預報[9-10]，但計算復雜，實時高精度的預報難以實現；基于先驗衛星軌道位置，采用多項式的方法可以進行實時、短期高精度的預報，目前切比雪夫多項式、最佳平方逼近多項式等方法均取得了較高的短期預報精度[4-8].但目前關于克里金算法在LEO軌道擬合與預報的研究尚未發現，本文詳細研究了滑動切比雪夫多項式、克里金算法在不同類型LEO軌道的插值擬合精度，以及兩種算法的短期軌道預報效果，以其為未來LEO衛星相關應用做出些許貢獻.

1 算法原理

1.1 滑動切比雪夫多項式算法

采用切比雪夫多項式擬合全球衛星導航系統(GNSS)衛星軌道時，待擬合點位于擬合弧段中間部分可以獲得高穩定、高精度的位置坐標[11-12]. 滑動切比雪夫多項式算法的實質及流程可概述為：根據軌道待擬合點的時間，選擇合適的擬合軌道弧段，使得待擬合點位于弧段中間，再采用切比雪夫多項式計算待擬合點的位置坐標[13-15]，其計算原理如下：

1) 切比雪夫多項式的自變量區間為[-1，1]，因此需將擬合軌道弧段內各點的時間t歸化到[-1，1].設擬合軌道弧段對應的時間段為[t1，t2]，則t對應的變量τ可表示為：

2) 根據切比雪夫多項式，衛星各個時刻t對應的坐標可以表示為

式中：QX,i、QY,i、QZ,i分別表示衛星參與擬合各歷元X、Y、Z坐標分量對應的切比雪夫系數；n為切比雪夫多項式的階數；Ti(τ) 為切比雪夫多項式，其計算方法為

3)根據式(2)、(3)分別構建X、Y、Z方向方程及矩陣，求解擬合弧段內各參與擬合歷元對應的系數.此處以X方向為例，假設參與擬合的歷元數目為m(m≥n)，則切比雪夫多項式矩陣T可表示為

同時切比雪夫系數矩陣、及X坐標矩陣可分別表示為QX、Xf：

根據最小二乘原理，切比雪夫系數可表示為式(6)，詳細解算原理可參考[10]

同理可求解出Y、Z方向對應的切比雪夫系數.

4)根據原理3)中求得的切比雪夫系數，將待擬合點的時間t帶入式(1)、(2)、(3)，求出待擬合點的位置坐標.

1.2 克里金插值算法

克里金算法由南非采礦工程師 D.G.Krige于1951年首次提出，是一種求最優、線形、無偏的空間內插方法[16]，該方法充分考慮到參與擬合各點之間、待插值點與各參與擬合點之間的空間相互關系，對每一個參與擬合的點賦予一定的權重系數，加權得到待插值點值. 同時，已有的研究表明當待插值點位于擬合弧段中間部分時，克里金插值算法獲得的精度較高[17]，本文在利用克里金插值時讓待插值點位于所使用弧段中間部分. 以X方向為例，利用克里金插值算法內插精密軌道的原理及流程如下：

1)根據待插值點的時間、以及參與擬合點的數目n，提取相應參與擬合點的歷元時間及空間坐標.

2)根據參與擬合點的歷元時間Ti，計算半變異函數值r*(h) ，其計算公式為

式中：h為擬合點對應的時間間隔；Z(Ti) 、Z(Ti+h) 分別表示Ti、Ti+h時 刻對應軌道的X坐標；N(h) 表示時間間隔為h的樣本點對總數.

3)根據半變異函數值r*(h)擬合理論變差函數，克里金算法中常用的理論變差函數模型有球狀模型、高斯模型、冪函數模型等. 本文采用高斯模型，其函數模型為

式中：r(h) 表 示半變異函數值；h含義同式(7)；a、c分別表示變程、基臺值.

4)克里金模型求待插值點的方程可表示為

式中， λi為第i個參與擬合點對應的權重系數. 由無偏估計性質可以得

根據步驟3)中擬合得到的理論變差函數模型，

計算參與擬合點中第i點和第k點對間距離對應的半變異函數值rij，以及參與擬合點與待插值點的半變異函數值riX，根據式(10)、rij、riX構建克里金方程組，求解系數 λi. 其中克里金方程組可以表示為

式中：

5)根據式(12)求解得到各個擬合點對應的權重系數，再利用式 (9) 得到待插值點的插值結果.

2 數據來源

為研究滑動切比雪夫多項式、克里金算法在不同類型LEO軌道的插值擬合及預報的精度，本文選取了GRACE-B、Jason-3、Sentinel-1B、HY-2A衛星軌道作為研究對象. 其中GRACE-B、Sentinel-1B衛星的精密軌道來自瑞士CODE中心，HY-2A和Jason-3衛星的精密軌道數據采用法國國家空間中心(CENS)發布的事后精密軌道，本文在進行上述研究時軌道采樣間隔均采用60 s，擬合弧段為24 h，各顆衛星的軌道信息如表1所示.

表1 實驗選用LEO衛星軌道信息

圖1為GRACE-B、Jason-3衛星24 h的精密軌道時間序列，從圖中可以發現無論軌道高度及軌道類型，LEO衛星的軌道整體平滑且有明顯的周期性.

圖1 GRACE-B、Jason-3精密軌道時間序列

3 擬合精度分析

本文在進行滑動切比雪夫擬合、克里金插值時，選取待求點前后相同數目的點進行擬合，精密星歷給出的待求點坐標作為真值，將擬合插值求得的待求點結果與真值比較得到誤差結果. 其中擬合軌道的空間點位誤差均方根(RMS)計算公式為

式中： Δxi、 Δyi、 Δzi表示待求點擬合插值結果與對應精密星歷的坐標差；k為待求點的總數.

3.1 滑動切比雪夫擬合精度分析

根據滑動切比雪夫插值策略，本文分析了6、8、10、12、14階切比雪夫多項式在4顆LEO衛星軌道的擬合精度. 表2為采用6階切比雪夫多項式擬合時，各顆LEO衛星軌道的擬合點位空間誤差；圖2為各顆LEO衛星軌道的擬合點位空間誤差；表2中[6, 8]表示擬合階數為6，所用擬合點數目為8，下文含義相同.

表2 各顆LEO衛星軌道的擬合點位空間誤差RMS cm

由表2、圖2分析可知，對于各顆LEO衛星均存在如下規律：軌道高度越高，擬合誤差越小. 在擬合階數確定時，參與擬合點的數目與階數相近時，擬合誤差較小，隨著擬合點數目的增加，擬合的誤差逐漸變大；在給定擬合點數目時，當階數小于擬合點數目時，隨著擬合階數的增加，擬合誤差逐漸變小. 同時也可發現，使用6階切比雪夫多項式進行擬合時，各顆衛星的擬合精度可以達到厘米級；各顆LEO衛星采用[8, 8]、[10, 10]、[12, 12]、[14, 14]策略擬合軌道的精度逐漸提高，但提升幅度極小.

圖2 各顆LEO衛星軌道的擬合點位空間誤差RMS

圖3～4分別為HY-2A、GRACE-B采用[8, 8]、[12, 12]策略擬合時的誤差序列. 從圖中可以發現，采用相同的擬合策略時，軌道高度越高，擬合誤差越小且分布更集中；對于同一顆衛星，擬合階數及擬合點數目較多時，誤差分布更集中.

圖4 GRACE-B采用[8, 8](左)、[12, 12](右)策略擬合時的誤差序列

3.2 克里金擬合精度分析

本文分別選取6、8、10、12、14、16個歷元來進行克里金擬合，分析各顆LEO衛星軌道擬合精度.表3、圖5為4顆LEO衛星采用不同歷元數目進行克里金擬合的軌道空間誤差，結合表3、圖5可知，采用6、8個歷元進行克里金擬合時，各顆LEO衛星的軌道精度在20～50 m，只能滿足LEO衛星測運控等低精度應用對軌道精度的要求. 同時也可發現，克里金擬合的精度對參與擬合的歷元數目特別敏感，采用6、8個歷元進行克里金擬合只能到達數十米級精度，當增加到10個歷元時卻可以獲得毫米級精度且精度達到最高；同時也需注意，在參與擬合的歷元數目大于10時，隨著歷元數目的增加，克里金擬合的精度逐漸降低.

表3 采用6、8個歷元進行克里金擬合軌道的空間誤差RMS m

圖5 采用10、12、14、16個歷元進行克里金擬合軌道的空間誤差RMS

圖6為Sentinel-1B(a)、Jason-3(b)采用10個歷元進行克里金擬合時，軌道在X、Y、Z方向及空間的擬合誤差序列. 結合圖1分析發現，采用克里金擬合時，各個插值點在X、Y、Z方向的擬合誤差與插值點在軌道中的位置表現出明顯的相關性；同時也可發現，相較于滑動切比雪夫擬合算法，克里金算法擬合軌道的空間誤差分布更為集中，未隨著歷元變化出現大幅波動.

圖6 Sentinel-1B、Jason-3采用10個歷元進行克里金擬合的誤差序列

4 預報精度分析

LEO衛星軌道高度低、運動速度極快、且在大氣層內，同時軌道所受攝動復雜，故LEO衛星軌道不可外推時間太長. 此節分析了切比雪夫多項式、克里金算法在不同LEO衛星軌道上短時間內的預報精度.

4.1 切比雪夫多項式預報精度分析

圖7為各顆LEO衛星采用不同切比雪夫多項式策略進行外推1歷元的軌道空間誤差. 由圖7可知，采用不同的擬合階數、擬合點數目對LEO衛星軌道的預報精度影響很大；對于各顆LEO衛星，采用策略[8, 10]進行外推的軌道精度最高. 同時，結合表1可以發現，在采取相同策略進行預報時，軌道高度越高，預報的精度也越高.

圖7 采用不同切比雪夫策略外推1歷元的軌道空間誤差

同時根據圖7也可發現，在采用8階及以上階數切比雪夫多項式進行預報時，在階數固定的情況下，隨著參與擬合點數目的增加，預報精度逐漸提高，但當擬合點數目增加到一定程度時，預報精度會隨著擬合點數目的增加而逐漸降低；在擬合點數目固定的情況下，采用低階切比雪夫多項式的預報精度更高.

表4為采用[8，10]策略進行預報時各顆LEO衛星對應歷元的軌道預報精度. 由表4可知，對于軌道高度較高的HY-2A、Jason-3衛星，120 s內預報的軌道精度在厘米級，外推240 s的軌道精度可以保持在5 m以內；對于軌道高度較低的GRACE-B、Sentinel-1B衛星，120 s內預報的軌道精度在5 m以內，外推360 s的軌道精度可以保持在100 m以內；值得注意隨著外推時間的增加，外推的精度急劇下降，很難用切比雪夫多項式模型對軌道進行逼近.

表4 采用[8,10]策略切比雪夫多項式外推軌道的空間誤差

4.2 克里金算法預報精度分析

本節選取6、8、10、12、14、16、20、30個連續歷元來分別進行克里金擬合外推各顆LEO衛星軌道. 分析發現采用6、8個歷元來進行預報的精度很差，外推1個歷元的誤差在100 m以上，已無法滿足大部分應用對LEO衛星軌道精度的要求. 圖8為各顆LEO衛星采用10個及以上歷元進行克里金算法外推1歷元的軌道空間誤差. 從圖中可以發現，采用相同歷元數目進行預報時，軌道高度越高，預報精度越高；隨著參與擬合歷元數目的增加，克里金預報精度逐漸提高，但當擬合點數目增加到一定程度時，預報精度會隨著擬合歷元數目的增加而逐漸降低；采用20個歷元進行預報1歷元(60 s)時，各顆LEO衛星的軌道預報精度在1～2.5 m.

圖8 采用不同擬合點數目克里金算法外推1歷元的軌道空間誤差

圖9為采用20個歷元進行克里金預報1歷元(60 s)時，各顆LEO衛星預報軌道的空間誤差分布序列. 由圖9可知，軌道高度越高，預報軌道的誤差越小且分布更集中.

圖9 采用20歷元進行克里金外推1歷元的軌道空間誤差序列

表5為采用20個歷元進行克里金預報時各顆LEO衛星對應外推各歷元的軌道精度，由表5可知采用克里金算法外推時，隨著外推時間的增加，外推軌道的精度急劇下降；外推240 s的軌道精度為百米級別，只能滿足LEO衛星低精度應用的需求. 對比表4～5可知：克里金算法外推LEO軌道的精度低于滑動切比雪夫算法.

表5 采用20個歷元進行克里金預報時外推軌道的空間誤差

5 結 語

本文利用60 s采樣間隔LEO精密星歷數據詳細研究了滑動切比雪夫多項式、克里金算法在不同類型LEO軌道的插值擬合精度，以及兩種算法的短期外推軌道精度，研究結果表明：

1)運用滑動切比雪夫算法進行LEO軌道擬合時，擬合階數接近擬合點數目取得的插值精度相對較高；在同時考慮計算量、插值精度的情況下，推薦使用[8,8]策略進行擬合，插值精度優于4 mm；采用[6,6]策略進行擬合時，各顆衛星的軌道插值可以獲得優于5 cm的精度，仍可以滿足大部分應用的需求.同時也需注意：滑動切比雪夫多項式插值精度與插值點的空間位置密切相關.

2)克里金擬合的精度對參與擬合的歷元數目特別敏感，采用6、8個歷元進行克里金擬合的精度在數十米級別；采用在擬合歷元數為10時，插值精度最高并優于6 mm. 相較于滑動切比雪夫多項式，克里金算法擬合軌道的空間誤差分布更為集中，未隨著歷元變化出現大幅波動.

3)總體上，克里金算法外推LEO軌道的精度低于滑動切比雪夫算法；在采取相同策略進行預報時，衛星軌道高度越高，預報的精度也越高；隨著外推時間的增加，兩種算法外推軌道的精度急劇下降. 采用克里金算法預報60 s，各顆LEO軌道預報的精度在1～2.5 m；采用滑動切比雪夫多項式預報120 s，可獲得優于5 m的軌道精度.