兩種北斗衛星精密星歷外推方法的精度分析

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

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

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

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

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

1 兩種精密星歷軌道坐標外推模型

1.1 三角函數插值多項式

BDS衛星繞地球運動過程中是有規律性的,呈周期性變化。三角函數插值法能夠很好地反映衛星的這種周期性運動規律,使插值精度更高。但衛星在地固坐標系下運動時,精密星歷引入了地球自身運動,影響衛星的運動規律。而在慣性坐標系下,衛星運動更能反映衛星本身的運動特點,更有規律。運用三角函數插值法對衛星精密星歷進行外推,首先將精密星歷的坐標由地固系轉換為慣性坐標系,轉換方程[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為指向的協議地球坐標;(X,Y,Z)為協議天球坐標系;M、E、N、P分別為極移旋轉矩陣、旋轉春分點時角、章動旋轉矩陣、歲差旋轉矩陣;xp、yp為極移分量;GAST為格林尼治時角;Δε、ΔΨ、ε分別為交角章動、黃經章動、黃赤交角;ξA、ZA、θA為歲差參數;R1、R2、R3為旋轉矩陣,其表達式為

(2)

式中,α為旋轉角。精密星歷軌道坐標轉換為慣性系坐標后,運用三角函數插值法進行插值,得到慣性系坐標插值結果,再轉換為地固系坐標,與已知坐標作差后進行比較。

三角函數插值多項式的形式[12]為

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

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

(3)

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

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

(4)

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

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

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

(5)

式中,n為已知衛星軌道坐標點數量;a1、a2、a3為待求系數;tn為衛星軌道坐標對應的時間點;xn為tn時刻對應的衛星軌道坐標在X方向的分量。上述模型為有約束優化模型,把上述模型展開求解,即按約束條件可得

(6)

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

(7)

則按

(8)

可寫成矩陣形式

(9)

其中

由式(9)得

(10)

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

(11)

同理可得BDS衛星在Y和Z方向軌道坐標分量的外推公式。

2 兩種衛星軌道坐標外推方法精度分析

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

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

2.1 先驗數據個數為13個時兩種外推方法精度分析

2.1.1 GEO衛星精密星歷外推精度分析

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

表1 GEO衛星X方向三角函數外推精度 cm

表2 GEO衛星X方向廣義延拓外推精度 cm

2.1.2 IGSO衛星精密星歷外推精度分析

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

表3 IGSO衛星X方向三角函數外推精度 cm

表4 IGSO衛星X方向廣義延拓外推精度 cm

2.1.3 MEO衛星精密星歷外推精度分析

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

表5 MEO衛星X方向三角函數外推精度 cm

表6 MEO衛星X方向廣義延拓外推精度 cm

2.2 不同先驗數據個數時兩種外推方法最優精度分析

2.2.1 GEO衛星精密星歷外推最優精度分析

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

表7 GEO衛星X方向三角函數外推最優精度統計 cm

表8 GEO衛星X方向廣義延拓外推最優精度統計 cm

2.2.2 IGSO衛星精密星歷外推最優精度分析

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

表9 IGSO衛星X方向三角函數外推最優精度統計 cm

表10 IGSO衛星X方向廣義延拓外推最優精度統計 cm

2.2.3 MEO衛星精密星歷外推最優精度分析

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

表11 MEO衛星X方向三角函數外推最優精度統計 cm

表12 MEO衛星X方向廣義延拓外推最優精度統計 cm

3 結 論

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

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

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

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