申少飛,雷偉偉,李振南,2,馬晨陽

(1. 河南理工大學(xué)測繪與國土信息工程學(xué)院,河南 焦作 454003; 2. 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院,湖北 武漢 430077)

北斗衛(wèi)星導(dǎo)航系統(tǒng)(BeiDou Navigation Satellite System,BDS)是由中國自行研制、面向全球的衛(wèi)星導(dǎo)航系統(tǒng)[1]。2020年6月23日,北斗三號最后一顆全球組網(wǎng)衛(wèi)星發(fā)射成功,標(biāo)志著北斗衛(wèi)星導(dǎo)航系統(tǒng)正式完成全球組網(wǎng),其在設(shè)計性能上優(yōu)于俄羅斯的GLONASS,且與GPS性能相當(dāng)[2]。北斗衛(wèi)星的定位精度可達分米級和厘米級,在定位過程中需要通過衛(wèi)星精密星歷獲取衛(wèi)星的坐標(biāo)信息[3]。GNSS服務(wù)組織提供的精密星歷采樣間隔一般為15或5 min,通常接收機的采樣間隔為30、15或1 s[4]。為獲取時間間隔更短的連續(xù)歷元的衛(wèi)星坐標(biāo),通常采用插值和擬合的方法[5]。

常用的插值和擬合方法有拉格朗日插值法、牛頓插值法、切比雪夫多項式擬合和勒讓德多項式擬合[6-10]。只要選取合適的時段和階數(shù),插值和擬合精度都能達到厘米甚至毫米級,完全能滿足精度要求。但國際全球衛(wèi)星導(dǎo)航系統(tǒng)服務(wù)(International GNSS Service,IGS)僅提供一天中00:00:00—23:45:00的星歷坐標(biāo)。因此,為了得到23:45:00—24:00:00的星歷坐標(biāo)有兩種方法:第1種是將當(dāng)天和第2天的精密星歷連接,采用插值或擬合的方法求取對應(yīng)的星歷坐標(biāo);第2種方法是依據(jù)當(dāng)天的星歷數(shù)據(jù)直接進行坐標(biāo)外推求得星歷坐標(biāo)[11]。采用第1種方法,當(dāng)兩天中某個時段的精密星歷中斷時,得出的結(jié)果精度不高,因此,通常選擇第2種外推方法。采用第2種方法要設(shè)法提高星歷外推精度。文獻[12]采用滑動式拉格朗日多項式插值和三角函數(shù)對IGS精密星歷進行外推,三角函數(shù)有效利用了衛(wèi)星在空間運動的特點,外推精度明顯優(yōu)于滑動式拉格朗日插值,且三角函數(shù)在6階時外推效果最優(yōu),小于滑動式拉格朗日插值的9階。文獻[13]運用拉格朗日插值、切比雪夫多項式擬合和三角函數(shù)對GPS精密星歷進行內(nèi)插和外推,結(jié)果表明3種方法內(nèi)插精度大致相同,但三角函數(shù)的外推精度明顯優(yōu)于另外兩種方法。文獻[14]利用廣義延拓法對GPS精密星歷進行內(nèi)插和外推,并與拉格朗日插值法的結(jié)果進行對比,結(jié)果表明廣義延拓逼近法插值精度小于5 cm,在外推30 min后依然能保持較高精度,優(yōu)于拉格朗日插值法。文獻[15]繪制11階多項式和9階三角函數(shù)在23:00—24:00時間段的外推圖像,結(jié)果表明三角函數(shù)在24:00外推誤差小于10 cm,精度遠高于多項式外推結(jié)果,進一步驗證了三角函數(shù)具有良好的外推能力。

以上只是針對GPS精密星歷的軌道坐標(biāo)外推研究,而關(guān)于北斗衛(wèi)星研究的較少。不同于GPS,BDS是由地球靜止軌道(geostationary earth orbit,GEO)、中圓地球軌道(medium earth orbit,MEO)和傾斜地球同步軌道(inclined geosynchronous satellite orbit,IGSO)3種不同軌道的衛(wèi)星星座組成[16]。本文采用三角函數(shù)插值法和廣義延拓逼近法,分別對采樣間隔為15 和5 min的北斗衛(wèi)星精密星歷進行5 min間隔外推,并與原始精密星歷作差后進行精度分析,對比兩種方法的外推能力。

1 兩種精密星歷軌道坐標(biāo)外推模型

1.1 三角函數(shù)插值多項式

BDS衛(wèi)星繞地球運動過程中是有規(guī)律性的,呈周期性變化。三角函數(shù)插值法能夠很好地反映衛(wèi)星的這種周期性運動規(guī)律,使插值精度更高。但衛(wèi)星在地固坐標(biāo)系下運動時,精密星歷引入了地球自身運動,影響衛(wèi)星的運動規(guī)律。而在慣性坐標(biāo)系下,衛(wèi)星運動更能反映衛(wèi)星本身的運動特點,更有規(guī)律。運用三角函數(shù)插值法對衛(wèi)星精密星歷進行外推,首先將精密星歷的坐標(biāo)由地固系轉(zhuǎn)換為慣性坐標(biāo)系,轉(zhuǎn)換方程[17]為

(1)

其中

M=R2(-xp)R1(-yp)

E=R3(GAST)

N=R1(-ε-Δε)R3(-ΔΨ)R1(ε)

P=R3(-90-ZA)R1(θA)R3(ξA)

式中,(x,y,z)為以CTP為指向的協(xié)議地球坐標(biāo);(X,Y,Z)為協(xié)議天球坐標(biāo)系;M、E、N、P分別為極移旋轉(zhuǎn)矩陣、旋轉(zhuǎn)春分點時角、章動旋轉(zhuǎn)矩陣、歲差旋轉(zhuǎn)矩陣;xp、yp為極移分量;GAST為格林尼治時角;Δε、ΔΨ、ε分別為交角章動、黃經(jīng)章動、黃赤交角;ξA、ZA、θA為歲差參數(shù);R1、R2、R3為旋轉(zhuǎn)矩陣,其表達式為

(2)

式中,α為旋轉(zhuǎn)角。精密星歷軌道坐標(biāo)轉(zhuǎn)換為慣性系坐標(biāo)后,運用三角函數(shù)插值法進行插值,得到慣性系坐標(biāo)插值結(jié)果,再轉(zhuǎn)換為地固系坐標(biāo),與已知坐標(biāo)作差后進行比較。

三角函數(shù)插值多項式的形式[12]為

X′t=a0+a1sin(ωt)+a2cos(ωt)+a3sin(2ωt)+

a4cos(2ωt)+…+ancos(nωt/2)

(3)

式中,X′t為對應(yīng)歷元時刻t下精密星歷軌道坐標(biāo)X方向分量;n為三角函數(shù)對應(yīng)的階數(shù);ω=2π/T(T為BDS衛(wèi)星運行周期,其中GEO和IGSO衛(wèi)星周期為23 h 56 min,MEO衛(wèi)星周期為12 h 50 min);a0、a1、…、an為求得的三角函數(shù)系數(shù)。

根據(jù)已獲得的星歷數(shù)據(jù),對式(3)進行展開為

(4)

對式(4)線性方程組求解,得到待求系數(shù)a0、a1、…、an。將待求系數(shù)代入式(3)中,即可得到任意時刻BDS衛(wèi)星軌道在X方向的坐標(biāo)分量,同理可求得BDS衛(wèi)星在Y方向和Z方向的坐標(biāo)分量。

1.2 廣義延拓逼近法外推模型

已知一組數(shù)據(jù)(x1,t1),(x2,t2),…,(xi,ti),…,(xn,tn) ,若求取tn+1時刻xn+1的值,則令tn為最新時刻,采用外推算法求取tn+1時刻的值xn+1,現(xiàn)以衛(wèi)星坐標(biāo)在X方向分量為例,建立廣義延拓外推模型[18]為

(5)

式中,n為已知衛(wèi)星軌道坐標(biāo)點數(shù)量;a1、a2、a3為待求系數(shù);tn為衛(wèi)星軌道坐標(biāo)對應(yīng)的時間點;xn為tn時刻對應(yīng)的衛(wèi)星軌道坐標(biāo)在X方向的分量。上述模型為有約束優(yōu)化模型,把上述模型展開求解,即按約束條件可得

(6)

把式(6)代入最小二乘公式可得

(7)

則按

(8)

可寫成矩陣形式

(9)

其中

由式(9)得

(10)

求得待求系數(shù)a2、a3后,代入式(6)中,即可得a1,則可得外推公式為

(11)

同理可得BDS衛(wèi)星在Y和Z方向軌道坐標(biāo)分量的外推公式。

2 兩種衛(wèi)星軌道坐標(biāo)外推方法精度分析

為驗證兩種衛(wèi)星軌道坐標(biāo)外推方法的精度,本文從IGS官網(wǎng)下載了2021年10月10日北斗衛(wèi)星軌道坐標(biāo)的精密星歷,選取PC16(GEO衛(wèi)星)、PC22(IGSO衛(wèi)星)及PC36(MEO衛(wèi)星)作為分析對象。歷元時刻為00:00:00—24:00:00,采樣時間間隔為15和5 min,選取的先驗數(shù)據(jù)點的個數(shù)為13、17、21和25。以先驗數(shù)據(jù)為已知數(shù)據(jù),運用兩種坐標(biāo)外推方法構(gòu)建數(shù)學(xué)模型,每隔5 min外推23:50:00—24:00:00時段的衛(wèi)星精密星歷。

假如采樣間隔為15 min,先驗數(shù)據(jù)點為13個,則在歷元時刻為20:45:00—23:45:00的時段中,每隔15 min選取衛(wèi)星精密星歷數(shù)據(jù),共選取13個作為先驗數(shù)據(jù),每隔5 min外推23:50:00—24:00:00時段的衛(wèi)星精密星歷,與外推點處歷元對應(yīng)的衛(wèi)星坐標(biāo)作差后進行比較,求出誤差絕對值的最大值和誤差中誤差,并求出相應(yīng)的點位中誤差。本文試驗運算環(huán)境為Matlab 9.8。IGS精密星歷數(shù)據(jù)誤差一般小于5 cm,因此兩種外推算法的精度至少小于5 cm才能滿足要求。

2.1 先驗數(shù)據(jù)個數(shù)為13個時兩種外推方法精度分析

2.1.1 GEO衛(wèi)星精密星歷外推精度分析

表1是先驗數(shù)據(jù)為13個時,GEO衛(wèi)星不同時間間隔X方向三角函數(shù)外推精度。當(dāng)時間間隔為15 min時,隨外推階數(shù)的增大,X方向差值和中誤差的精度不斷提高,在10階時分別為4.15和2.98 cm;當(dāng)時間間隔為5 min時,X方向差值和中誤差則隨著外推階數(shù)的增大而不斷增大,在6階時精度最高,分別為0.44和0.28 cm,比15 min時間間隔的最優(yōu)精度高出了1個數(shù)量級。表2是先驗數(shù)據(jù)為13時,GEO衛(wèi)星不同時間間隔X方向廣義延拓外推精度。當(dāng)時間間隔為15 min時,外推精度隨階數(shù)的增大呈現(xiàn)先增高后降低的趨勢,在第9階精度最高;當(dāng)時間間隔為5 min時,外推精度在7階時最高,差值和中誤差分別為1.23和0.77 cm,比15 min時間間隔的外推精度有了很大的提升,能夠滿足精度的要求。

表1 GEO衛(wèi)星X方向三角函數(shù)外推精度 cm

表2 GEO衛(wèi)星X方向廣義延拓外推精度 cm

2.1.2 IGSO衛(wèi)星精密星歷外推精度分析

表3是先驗數(shù)據(jù)為13個時,IGSO衛(wèi)星不同時間間隔X方向三角函數(shù)外推精度。當(dāng)時間間隔為15 min時,隨外推階數(shù)的增大,X方向三角函數(shù)外推的差值和中誤差精度總體上呈現(xiàn)不斷增高的趨勢,在10階時精度達到最高;當(dāng)時間間隔為5 min時,X方向三角函數(shù)外推的差值和中誤差隨階數(shù)的增大先減小后增大,在7階時達到最小,與15 min時間間隔最優(yōu)精度相比,精度相近。表4除了在5 min時間間隔的8階時IGSO衛(wèi)星X方向廣義延拓外推中誤差為4.88 cm,其余外推結(jié)果精度均較低,不能夠滿足衛(wèi)星精密星歷軌道坐標(biāo)的精度要求。

表3 IGSO衛(wèi)星X方向三角函數(shù)外推精度 cm

表4 IGSO衛(wèi)星X方向廣義延拓外推精度 cm

2.1.3 MEO衛(wèi)星精密星歷外推精度分析

表5是先驗數(shù)據(jù)為13個時, MEO衛(wèi)星X方向三角函數(shù)外推精度。當(dāng)時間間隔為15 min時,隨外推階數(shù)的增大不斷提高,在10階時精度達最高;5 min時間間隔的X方向誤差和中誤差在7階時精度最高,小于5 cm,符合精度要求。表6是先驗數(shù)據(jù)為13時,MEO衛(wèi)星不同時間間隔X方向廣義延拓外推精度。在15和5 min時間間隔時,X方向廣義延拓外推的差值和中誤差均隨階數(shù)增大先減小后增大,但無論是15 min,還是精度更高的5 min時間間隔的外推結(jié)果,其精度均不能滿足要求。

表5 MEO衛(wèi)星X方向三角函數(shù)外推精度 cm

表6 MEO衛(wèi)星X方向廣義延拓外推精度 cm

2.2 不同先驗數(shù)據(jù)個數(shù)時兩種外推方法最優(yōu)精度分析

2.2.1 GEO衛(wèi)星精密星歷外推最優(yōu)精度分析

表7為 GEO衛(wèi)星X方向三角函數(shù)外推精度的比較?？梢钥闯?隨著先驗數(shù)據(jù)個數(shù)的增加,15 min時間間隔的差值、中誤差和點位中誤差精度相近,且均小于5 cm;隨著先驗數(shù)據(jù)時間間隔的縮短,當(dāng)時間間隔為5 min時,精度有了較大的提高。先驗數(shù)據(jù)為13和17個時,X方向差值和中誤差精度都達到了毫米級。表8為GEO衛(wèi)星X方向廣義延拓外推精度的比較?？梢钥闯?相較于三角函數(shù)外推精度,15 min時間間隔的精度較低,且都大于5 cm,不能滿足精度的要求;相較于15 min時間間隔的外推精度,5 min時間間隔的外推精度達厘米級,先驗數(shù)據(jù)為13個,在7階時精度最高,差值、中誤差和點位中誤差分別為1.23、0.77和2.65 cm。

表7 GEO衛(wèi)星X方向三角函數(shù)外推最優(yōu)精度統(tǒng)計 cm

表8 GEO衛(wèi)星X方向廣義延拓外推最優(yōu)精度統(tǒng)計 cm

2.2.2 IGSO衛(wèi)星精密星歷外推最優(yōu)精度分析

表9為IGSO衛(wèi)星X方向三角函數(shù)外推精度的統(tǒng)計。當(dāng)時間間隔為15 min、先驗數(shù)據(jù)為13和25個的最優(yōu)階時,差值和中誤差的精度可以達到毫米級;當(dāng)時間間隔為5 min時,隨著先驗數(shù)據(jù)的增加,精度越來越高,先驗數(shù)據(jù)為25個時達到最優(yōu)精度。表10為IGSO衛(wèi)星X方向廣義延拓外推最優(yōu)精度的統(tǒng)計,只有在5 min時間間隔時,精度較高。在先驗數(shù)據(jù)為25個、最優(yōu)階數(shù)為9階時,差值、中誤差和點位中誤差分別為0.3、0.21和3.23 cm。

表9 IGSO衛(wèi)星X方向三角函數(shù)外推最優(yōu)精度統(tǒng)計 cm

表10 IGSO衛(wèi)星X方向廣義延拓外推最優(yōu)精度統(tǒng)計 cm

2.2.3 MEO衛(wèi)星精密星歷外推最優(yōu)精度分析

表11為MEO衛(wèi)星X方向三角函數(shù)外推最優(yōu)精度的結(jié)果統(tǒng)計。當(dāng)時間間隔為15 min、先驗數(shù)據(jù)為21個,最優(yōu)階為14時,差值為6.3 cm,精度最低;其余3個不同先驗數(shù)據(jù)個數(shù)最優(yōu)階精度相近,且都能符合精度要求。5 min時間間隔的外推精度隨著先驗數(shù)據(jù)的增加而提高,在先驗數(shù)據(jù)為25個時,差值、中誤差和點位中誤差分別為1.61、0.95和2.01 cm。表12為MEO衛(wèi)星X方向廣義延拓外推最優(yōu)精度的結(jié)果統(tǒng)計。當(dāng)時間間隔為15 min時,精度達到米級,不符合精度要求;隨著外推時間間隔的減小,時間間隔為5 min時,精度隨著先驗數(shù)據(jù)個數(shù)的增加先增高后降低,在先驗數(shù)據(jù)為17個時,差值和中誤差的精度能達到毫米級。

表11 MEO衛(wèi)星X方向三角函數(shù)外推最優(yōu)精度統(tǒng)計 cm

表12 MEO衛(wèi)星X方向廣義延拓外推最優(yōu)精度統(tǒng)計 cm

3 結(jié) 論

本文基于IGS官網(wǎng)發(fā)布的15和5 min時間間隔的北斗衛(wèi)星精密星歷軌道坐標(biāo),運用三角函數(shù)插值多項式和廣義延拓逼近法數(shù)學(xué)模型,對2021年10月10日PC16、PC22和PC36 3顆北斗衛(wèi)星精密星歷進行了5 min時間間隔的軌道坐標(biāo)外推計算和精度分析。結(jié)果表明:

(1)對不同類型北斗衛(wèi)星精密星歷進行軌道坐標(biāo)外推,即使是同一種方法,精度也不同。 GEO衛(wèi)星相較于IGSO衛(wèi)星和MEO衛(wèi)星,運用同一種方法進行軌道坐標(biāo)外推,無論是15 min還是5 min時間間隔,精度都相對更高。

(2)兩種外推方法在對同一顆北斗衛(wèi)星進行軌道坐標(biāo)外推時,三角函數(shù)插值多項式相比于廣義延拓逼近法,在15 min時間間隔時,三角函數(shù)外推精度比廣義延拓逼近法高1～2個數(shù)量級;在5 min時間間隔時,只要選擇合適的階數(shù),兩種外推方法精度都能達到厘米甚至毫米級,但是三角函數(shù)相較于廣義延拓逼近法精度更高。

(3)兩種外推方法對不同類型北斗衛(wèi)星精密軌道坐標(biāo)進行外推計算時,5 min時間間隔的計算精度相較于15 min時間間隔,精度更高,廣義延拓逼近法5 min時間間隔的計算精度相較于15 min時間間隔,精度能高出1～2個數(shù)量級。因此,可以優(yōu)先選用5 min時間間隔的衛(wèi)星精密星歷進行軌道坐標(biāo)外推計算。