不同軌道類型LEO衛(wèi)星軌道擬合及預(yù)報(bào)精度研究

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

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

0 引 言

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

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

1 算法原理

1.1 滑動切比雪夫多項(xiàng)式算法

采用切比雪夫多項(xiàng)式擬合全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)衛(wèi)星軌道時(shí)，待擬合點(diǎn)位于擬合弧段中間部分可以獲得高穩(wěn)定、高精度的位置坐標(biāo)[11-12]. 滑動切比雪夫多項(xiàng)式算法的實(shí)質(zhì)及流程可概述為：根據(jù)軌道待擬合點(diǎn)的時(shí)間，選擇合適的擬合軌道弧段，使得待擬合點(diǎn)位于弧段中間，再采用切比雪夫多項(xiàng)式計(jì)算待擬合點(diǎn)的位置坐標(biāo)[13-15]，其計(jì)算原理如下：

1) 切比雪夫多項(xiàng)式的自變量區(qū)間為[-1，1]，因此需將擬合軌道弧段內(nèi)各點(diǎn)的時(shí)間t歸化到[-1，1].設(shè)擬合軌道弧段對應(yīng)的時(shí)間段為[t1，t2]，則t對應(yīng)的變量τ可表示為：

2) 根據(jù)切比雪夫多項(xiàng)式，衛(wèi)星各個(gè)時(shí)刻t對應(yīng)的坐標(biāo)可以表示為

式中：QX,i、QY,i、QZ,i分別表示衛(wèi)星參與擬合各歷元X、Y、Z坐標(biāo)分量對應(yīng)的切比雪夫系數(shù)；n為切比雪夫多項(xiàng)式的階數(shù)；Ti(τ) 為切比雪夫多項(xiàng)式，其計(jì)算方法為

3)根據(jù)式(2)、(3)分別構(gòu)建X、Y、Z方向方程及矩陣，求解擬合弧段內(nèi)各參與擬合歷元對應(yīng)的系數(shù).此處以X方向?yàn)槔?，假設(shè)參與擬合的歷元數(shù)目為m(m≥n)，則切比雪夫多項(xiàng)式矩陣T可表示為

同時(shí)切比雪夫系數(shù)矩陣、及X坐標(biāo)矩陣可分別表示為QX、Xf：

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

同理可求解出Y、Z方向?qū)?yīng)的切比雪夫系數(shù).

4)根據(jù)原理3)中求得的切比雪夫系數(shù)，將待擬合點(diǎn)的時(shí)間t帶入式(1)、(2)、(3)，求出待擬合點(diǎn)的位置坐標(biāo).

1.2 克里金插值算法

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

1)根據(jù)待插值點(diǎn)的時(shí)間、以及參與擬合點(diǎn)的數(shù)目n，提取相應(yīng)參與擬合點(diǎn)的歷元時(shí)間及空間坐標(biāo).

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

式中：h為擬合點(diǎn)對應(yīng)的時(shí)間間隔；Z(Ti) 、Z(Ti+h) 分別表示Ti、Ti+h時(shí) 刻對應(yīng)軌道的X坐標(biāo)；N(h) 表示時(shí)間間隔為h的樣本點(diǎn)對總數(shù).

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

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

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

式中， λi為第i個(gè)參與擬合點(diǎn)對應(yīng)的權(quán)重系數(shù). 由無偏估計(jì)性質(zhì)可以得

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

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

式中：

5)根據(jù)式(12)求解得到各個(gè)擬合點(diǎn)對應(yīng)的權(quán)重系數(shù)，再利用式 (9) 得到待插值點(diǎn)的插值結(jié)果.

2 數(shù)據(jù)來源

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

表1 實(shí)驗(yàn)選用LEO衛(wèi)星軌道信息

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

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

3 擬合精度分析

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

式中： Δxi、 Δyi、 Δzi表示待求點(diǎn)擬合插值結(jié)果與對應(yīng)精密星歷的坐標(biāo)差；k為待求點(diǎn)的總數(shù).

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

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

表2 各顆LEO衛(wèi)星軌道的擬合點(diǎn)位空間誤差RMS cm

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

圖2 各顆LEO衛(wèi)星軌道的擬合點(diǎn)位空間誤差RMS

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

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

3.2 克里金擬合精度分析

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

表3 采用6、8個(gè)歷元進(jìn)行克里金擬合軌道的空間誤差RMS m

圖5 采用10、12、14、16個(gè)歷元進(jìn)行克里金擬合軌道的空間誤差RMS

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

圖6 Sentinel-1B、Jason-3采用10個(gè)歷元進(jìn)行克里金擬合的誤差序列

4 預(yù)報(bào)精度分析

LEO衛(wèi)星軌道高度低、運(yùn)動速度極快、且在大氣層內(nèi)，同時(shí)軌道所受攝動復(fù)雜，故LEO衛(wèi)星軌道不可外推時(shí)間太長. 此節(jié)分析了切比雪夫多項(xiàng)式、克里金算法在不同LEO衛(wèi)星軌道上短時(shí)間內(nèi)的預(yù)報(bào)精度.

4.1 切比雪夫多項(xiàng)式預(yù)報(bào)精度分析

圖7為各顆LEO衛(wèi)星采用不同切比雪夫多項(xiàng)式策略進(jìn)行外推1歷元的軌道空間誤差. 由圖7可知，采用不同的擬合階數(shù)、擬合點(diǎn)數(shù)目對LEO衛(wèi)星軌道的預(yù)報(bào)精度影響很大；對于各顆LEO衛(wèi)星，采用策略[8, 10]進(jìn)行外推的軌道精度最高. 同時(shí)，結(jié)合表1可以發(fā)現(xiàn)，在采取相同策略進(jìn)行預(yù)報(bào)時(shí)，軌道高度越高，預(yù)報(bào)的精度也越高.

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

同時(shí)根據(jù)圖7也可發(fā)現(xiàn)，在采用8階及以上階數(shù)切比雪夫多項(xiàng)式進(jìn)行預(yù)報(bào)時(shí)，在階數(shù)固定的情況下，隨著參與擬合點(diǎn)數(shù)目的增加，預(yù)報(bào)精度逐漸提高，但當(dāng)擬合點(diǎn)數(shù)目增加到一定程度時(shí)，預(yù)報(bào)精度會隨著擬合點(diǎn)數(shù)目的增加而逐漸降低；在擬合點(diǎn)數(shù)目固定的情況下，采用低階切比雪夫多項(xiàng)式的預(yù)報(bào)精度更高.

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

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

4.2 克里金算法預(yù)報(bào)精度分析

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

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

圖9為采用20個(gè)歷元進(jìn)行克里金預(yù)報(bào)1歷元(60 s)時(shí)，各顆LEO衛(wèi)星預(yù)報(bào)軌道的空間誤差分布序列. 由圖9可知，軌道高度越高，預(yù)報(bào)軌道的誤差越小且分布更集中.

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

表5為采用20個(gè)歷元進(jìn)行克里金預(yù)報(bào)時(shí)各顆LEO衛(wèi)星對應(yīng)外推各歷元的軌道精度，由表5可知采用克里金算法外推時(shí)，隨著外推時(shí)間的增加，外推軌道的精度急劇下降；外推240 s的軌道精度為百米級別，只能滿足LEO衛(wèi)星低精度應(yīng)用的需求. 對比表4～5可知：克里金算法外推LEO軌道的精度低于滑動切比雪夫算法.

表5 采用20個(gè)歷元進(jìn)行克里金預(yù)報(bào)時(shí)外推軌道的空間誤差

5 結(jié) 語

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

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

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

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