董志華 陳俊平 周旭華

(1 中國科學(xué)院上海天文臺 上海 200030)

(2 中國科學(xué)院大學(xué) 北京 100049)

1 引言

國際GNSS (Global Navigation Satellite System)服務(wù)(IGS)[1]旨在為導(dǎo)航定位用戶提供高精度GNSS數(shù)據(jù)、產(chǎn)品和服務(wù)[2],其產(chǎn)品精度、一致性與連續(xù)性是高精度定位和導(dǎo)航的基礎(chǔ)和前提.按照目前的國際規(guī)范,高精度GNSS軌道產(chǎn)品一般以天為周期進行發(fā)布,提供給用戶使用,但連續(xù)多天的產(chǎn)品存在不同天之間的跳變問題.精密產(chǎn)品的跳變,將造成用戶精密位置計算的不精確,從而影響針對地震監(jiān)測、海平面變化等毫米級位置精度需求的應(yīng)用.此外,由于IGS各個分析中心(Analysis Centers,AC)采用的數(shù)據(jù)處理策略和力學(xué)模型不同,其產(chǎn)品連續(xù)性也存在一定差異.Griffiths等人于2009年首次提出日邊界不連續(xù)性(Day Boundary Discontinuities,DBD)分析法[3]進行軌道性能分析,之后許多學(xué)者也采用DBD進行GNSS精密產(chǎn)品精度的評定[4–7].DBD產(chǎn)生的原因在于動力學(xué)模型不完善等引起的軌道確定和預(yù)報的誤差.因此,通過分析DBD的長期序列,能夠評估目前GNSS軌道的精度情況,并能由其特性分析研究軌道確定過程中存在的誤差.論文利用德國地學(xué)研究中心(GFZ)、歐洲定軌中心(COD)、歐空局(ESA)、美國噴氣試驗室(JPL)以及上海天文臺(SHA)[8–9]5個GNSS分析中心2013—2017年的軌道產(chǎn)品,分析了軌道跳變的特性,發(fā)現(xiàn)跳變序列存在顯著的周期項,并驗證了地球反照輻射壓和太陽光壓模型等動力學(xué)模型對軌道的影響.

2 軌道連續(xù)性分析

2.1 軌道日邊界不連續(xù)

DBD是連續(xù)兩天的精密星歷軌道在某一相同歷元的3維坐標(biāo)差異.A、B兩天相同歷元軌道在徑向(R)、切向(T)和法向(N) 3個方向上的DBD,PDR、PDT、PDN可分別表示為:

其中, RA、TA、NA和RB、TB、NB分別是A、B兩天的坐標(biāo).

目前國際各GNSS分析中心提供的單天精密衛(wèi)星軌道產(chǎn)品所給的衛(wèi)星坐標(biāo)范圍為00:00:00—23:45:00,兩天數(shù)據(jù)間沒有重疊時段,DBD的計算需要將第1天(Day A)的軌道外推一個歷元至下一天(Day B)的零點,并與下一天SP3文件的零點軌道作差.

在利用圖1所示方案進行DBD分析時,需要將軌道外推15 min.軌道外推采用的方法為軌道動力學(xué)擬合及外推.具體分析方法為: (1)利用精密星歷數(shù)據(jù)中起點Tb至終點Te的精密軌道數(shù)據(jù),進行精密軌道動力學(xué)擬合; (2)根據(jù)軌道擬合得到的衛(wèi)星初始軌道根數(shù)以及太陽光壓等動力學(xué)參數(shù),利用軌道積分獲取指定弧段的衛(wèi)星軌道.

圖1 軌道DBD的示意圖,PD代表A、B連續(xù)兩天SP3文件的差值, Tb代表起始時刻, Te代表終止時刻.Fig.1 Illustration of the day boundary discontinuity in the satellite orbits,PD is the position difference between A and B successive SP3 daily files, Tb is the beginning time, Te is the end time.

2.2 精度評估方法

本文采用均方根(root mean square,RMS)誤差作為每顆衛(wèi)星DBD評估的指標(biāo)來分析與評估精密軌道產(chǎn)品的不連續(xù)性,對衛(wèi)星s,其公式為

式中, m為誤差分析的數(shù)據(jù)總天數(shù); xs為第i天該GNSS衛(wèi)星精密軌道產(chǎn)品的真值;為第i天該衛(wèi)星精密軌道的外推值或擬合值.對于軌道外推和軌道不連續(xù)性,上述公式即其m天內(nèi)衛(wèi)星s的誤差值; 對于軌道擬合,先用上述公式計算出1 d內(nèi)每顆衛(wèi)星軌道擬合的RMS值,再根據(jù)上述公式計算每顆衛(wèi)星在m天內(nèi)的RMS值.

3 數(shù)據(jù)分析

3.1 軌道外推精度分析

軌道外推的精度將影響DBD評估,首先進行軌道外推精度的分析.在以上軌道外推策略中,取Tb為A天00:00:00, Te為A天23:30:00,利用軌道外推方法獲取A天23:45:00時刻的衛(wèi)星軌道,并將該時刻的外推軌道與精密軌道的真實3維坐標(biāo)進行比較,得到兩者差值.使用BERNESE 5.2軟件[10]進行軌道外推處理,軌道積分采用的策略如表1所示.

表1 軌道積分策略Table 1 The strategy of orbits integration

采用GFZ、COD、ESA、SHA和JPL 5個分析中心2013—2017年5 yr的精密軌道,進行軌道外推,分析所有衛(wèi)星的外推精度情況.表2列出了所有GPS/GLONASS (GLOBAL NAVIGATION SATELLITE SYSTEM)衛(wèi)星軌道外推3維坐標(biāo)誤差5 yr的統(tǒng)計的平均值,表中,JPL只提供GPS的精密軌道,Exp (Extrapolation)為軌道外推精度,Fit為軌道擬合精度.

由表2可以看出,各個分析中心的3維軌道擬合精度平均為2.6 mm; 其中GFZ最大,為2.93 mm; COD最小,為2.28 mm; SHA和JPL相當(dāng),分別為2.65 mm和2.63 mm;ESA為2.44 mm.軌道外推精度平均為4.6 mm,其中JPL最大,為6.04 mm; ESA最小,為3.20 mm; GFZ、COD和SHA依次為5.62 mm、3.71 mm和4.49 mm.軌道擬合和軌道外推的平均精度都小于5 mm,相比軌道擬合,軌道外推增加了大約2 mm的預(yù)報誤差.

表2 各GNSS分析中心外推軌道3維誤差(單位: 毫米)Table 2 Three-dimensional error of the orbits extrapolation of different GNSS analysis centers (unit: mm)

3.2 軌道日邊界不連續(xù)性分析

以上分析表明,軌道15 min的預(yù)報誤差平均為4.6 mm.進一步分析GNSS軌道日邊界不連續(xù)性,根據(jù)圖1所示方案,取Tb為A天00:00:00, Te為A天23:45:00,利用軌道外推方法獲取該天24:00:00時刻的衛(wèi)星軌道,并將外推獲取的24:00:00時刻軌道與B天00:00:00時刻精密軌道的真實3維坐標(biāo)進行比較,得到兩者差值.利用該方法,對GFZ、COD、ESA、JPL以及SHA 5個分析中心2013—2017年5 yr的數(shù)據(jù)進行DBD分析.部分GPS/GLONASS衛(wèi)星5 yr DBD值序列的RMS情況如圖2–4所示.

圖2 GPS/GLONASS衛(wèi)星軌道徑向跳變值的RMS值Fig.2 RMS values of the GPS/GLONASS satellites’ radial orbit jump

由圖2–4可以看出,各個分析中心的軌道跳變均在厘米量級,遠(yuǎn)大于表2所統(tǒng)計出的平均為4.6 mm的軌道擬合和軌道外推殘差,所以軌道擬合和外推對于分析軌道不連續(xù)性帶來的誤差可以忽略不計.對所有GPS衛(wèi)星軌道跳變進行統(tǒng)計,各個分析中心R、T、N 3個方向的軌道跳變RMS平均值分別如下: GFZ為2.2 cm、2.7 cm、2.7 cm,ESA為2.1 cm、2.5 cm、2.1 cm,COD為0.7 cm、0.8 cm、0.7 cm,SHA為3.1 cm、4.5 cm、3.9 cm,JPL為1.5 cm、1.3 cm、1.5 cm.對所有GLONASS衛(wèi)星軌道跳變進行統(tǒng)計,各個分析中心R、T、N 3個方向的軌道跳變RMS平均值分別如下: GFZ為5.3 cm、7.6 cm、5.9 cm,ESA為2.5 cm、9.1 cm、6.8 cm,COD為0.7 cm、1 cm、1.2 cm,SHA為7.9 cm、11.9 cm、8.7 cm.以上差異表明各GNSS分析中心在數(shù)據(jù)處理及軟件策略方面仍存在一定差異.此外,徑向軌道跳變量小于其他兩個方向,這與軌道徑向精度高于其他兩個方向的情況一致.

圖3 GPS/GLONASS衛(wèi)星軌道切向跳變值的RMS值Fig.3 RMS values of the GPS/GLONASS satellites’ tangential orbit jump

圖4 GPS/GLONASS衛(wèi)星軌道法向跳變值的RMS值Fig.4 RMS values of the GPS/GLONASS satellites’ normal orbit jump

將5個分析中心各顆衛(wèi)星的3維軌道DBD進行RMS統(tǒng)計,如表3所示.對于GPS衛(wèi)星軌道3維跳變值的RMS統(tǒng)計,GFZ為4.48 cm、ESA為3.90 cm、COD為1.27 cm、SHA為6.77 cm、JPL為2.51 cm.對于GLONASS衛(wèi)星軌道跳變值的RMS統(tǒng)計,GFZ為11.11 cm、ESA為11.63 cm、COD為1.74 cm、SHA為16.72 cm.COD在生成精密軌道產(chǎn)品采用的軟件為本文使用的BERNESE軟件,這可能是COD的DBD最小的原因之一.這也表明目前國際上各GNSS分析中心在數(shù)據(jù)處理上還存在模型的差異.此外,GLONASS衛(wèi)星的軌道跳變均大于GPS衛(wèi)星,這與目前GLONASS衛(wèi)星精密產(chǎn)品精度相對較低一致.

表3 GPS/GLONASS衛(wèi)星軌道3維跳變值的平均RMS值(單位: 厘米)Table 3 Mean RMS values of the GPS/GLONASS satellites’ three-dimensional orbit jump (unit: cm)

圖5為COD、GFZ分析中心G29和G05衛(wèi)星2013—2017年5 yr的軌道跳變序列,其中STDR、STDT、STDN分別為R、T、N 3個方向軌道跳變的標(biāo)準(zhǔn)差(STD),RMSR、RMST、RMSN分別為R、T、N 3個方向軌道跳變的均方根.從中可以看出,兩個分析中心的軌道跳變在各個分量方向上均存在明顯的周期特性.采用最小二乘頻譜分析(least square spectrum analysis,LSSA)對軌道跳變的周期特性進行分析,并在頻譜圖峰值附近搜索主要周期[12].圖6給出了GFZ分析中心G29衛(wèi)星和COD分析中心G05衛(wèi)星精密軌道在R、T、N 3個分量上跳變的頻譜圖.從圖6可以看出,較顯著的周期約為90 d、120 d、175 d、340 d、352 d; 其中90 d、120 d和340 d的周期項對應(yīng)海潮對地球自轉(zhuǎn)的影響,其振幅為數(shù)毫米至1 cm左右,表明軌道確定中的潮汐改正模型需進一步改進; 而175 d和352 d則是與衛(wèi)星星座相關(guān)的1/2和1個交點年周期項[13–14].對其他衛(wèi)星分析,都存在與這兩顆星相近的周期特性.

圖5 GFZ/COD精密產(chǎn)品衛(wèi)星G29 (左),G05 (右)軌道跳變Fig.5 The orbit jump in GFZ/COD precise products of satellites G29 (left),G05 (right)

圖6 GFZ/COD精密產(chǎn)品衛(wèi)星G29 (左),G05 (右)軌道跳變頻譜分析Fig.6 The spectrum of the orbit jump in GFZ/COD precise products of satellites G29 (left),G05 (right)

4 力學(xué)模型對軌道外推精度的影響

分析DBD時采用的數(shù)據(jù)段2013—2017年5 yr,在此期間COD分析中心發(fā)生了兩次比較大的力學(xué)模型的更新,分別是2013年7月14日(年積日195 d)地球反照輻射模型[15]的引入以及2015年1月4日(年積日4 d)太陽光壓模型[16]的改變.其中,地球反照輻射是指太陽光照射至地表后部分反射回衛(wèi)星處對衛(wèi)星所引起的攝動.而導(dǎo)航衛(wèi)星的太陽光壓力是對軌道影響最大的非保守力,本文選取的數(shù)據(jù)時間段2013—2017年期間,COD采用的光壓模型從2015年1月4日開始由ECOM1 (reduced ECOM)變更為ECOM2 (extended ECOM)[16–17].其中,ECOM1模型[17–18]為:

式中, D0、Y0、B0、B1,C和B1,S代表ECOM1模型的5個光壓參數(shù), u代表衛(wèi)星在軌道平面上距升交點的角度.ECOM2模型[11,18]為:

式中, D2,C、D2,S、D4,C和D4,S代表ECOM2模型中相較ECOM1模型多出的4個光壓參數(shù).兩種光壓模型的光壓參數(shù)分別如表4所示.

表4 光壓模型及參數(shù)[16?17]Table 4 Solar radiation pressure model and parameters[16?17]

可以看到ECOM2模型在ECOM1的基礎(chǔ)上引入了包含2u和4u的周期項.為分析以上力學(xué)模型更新對軌道跳變分析的影響,對COD分析中心精密軌道的擬合和預(yù)報進行進一步分析.將表1中的太陽光壓模型變?yōu)镋COM1重新進行軌道擬合及預(yù)報,重新分析COD中心2013—2017年軌道外推精度; 并與上述ECOM2模型的結(jié)果進行對比,結(jié)果如圖7–8所示.

圖7 COD分析中心G07衛(wèi)星軌道外推精度,左圖采用ECOM1,右圖采用ECOM2.Fig.7 The precision of the orbits extrapolation of G07 satellite using ECOM1 (left) and ECOM2 (right)models.

圖8 COD分析中心R07衛(wèi)星軌道外推精度,左圖采用ECOM1,右圖采用ECOM2.Fig.8 The precision of the orbits extrapolation of R07 satellite using ECOM1 (left) and ECOM2 (right)models.

圖7和圖8結(jié)果表明,采用ECOM1光壓模型時,COD分析中心GPS/GLONASS衛(wèi)星的外推軌道從第734日(即2015年1月4日)精度急劇降低,這與COD更新太陽光壓模型的日期一致; 表明太陽光壓模型的不一致,將使軌道外推精度劣化至厘米級.采用ECOM2光壓模型時,COD分析中心GPS/GLONASS衛(wèi)星的外推軌道在第195日,即2013年第195日精度迅速下降,這與COD引入地球反照輻射壓的日期一致; 這表明,地球反照輻射壓攝動同樣會造成軌道外推的惡化; 并且可以看出衛(wèi)星徑向軌道存在系統(tǒng)偏差,表明地球反照輻射會引起GPS衛(wèi)星徑向軌道誤差[15].

以上對比說明,地球反照輻射壓、太陽光壓模型等的差異會對軌道連續(xù)性產(chǎn)生較大影響.對ECOM1模型下G07衛(wèi)星第734日后的3維外推軌道誤差進行頻譜分析,結(jié)果如圖9所示.可以看到,光壓模型的不一致,將為外推軌道引入1/4個交點年(即87 d左右)的周期誤差,這與(3)–(4)式中4u項的周期一致; 而公式中2u項對軌道的影響不夠顯著,表明該項對軌道的影響有限.

圖9 COD分析中心G07衛(wèi)星外推軌道誤差的頻譜分析Fig.9 The spectrum of the 3D orbits extrapolation of G07 satellite

5 結(jié)論

本文計算了GFZ、COD、ESA、SHA和JPL這5個分析中心在2013—2017年5 yr間精密軌道產(chǎn)品天與天交界處的跳變,分析了其產(chǎn)品的連續(xù)性.計算結(jié)果表明:GFZ、COD、ESA、SHA和JPL的平均3維軌道跳變分別為7.79 cm、1.51 cm、7.77 cm、11.75 cm和2.51 cm,這些差異說明目前國際上各分析中心對這兩個系統(tǒng)衛(wèi)星的數(shù)據(jù)處理模型還存在一定的差異; 在此基礎(chǔ)上,并對連續(xù)5 yr的軌道跳變序列進行最小二乘頻譜分析,發(fā)現(xiàn)了與潮汐改正模型周期相近的周期為90 d、120 d、340 d左右的顯著周期項誤差,表明潮汐模型可能是影響GNSS軌道確定精度的因素; 此外,GPS軌道跳變序列呈現(xiàn)出與衛(wèi)星星座相關(guān)的175 d和352 d左右的交點年顯著周期項; 而地球反照輻射壓和太陽光壓模型等動力學(xué)模型不一致,也會引起軌道外推精度的下降,為外推軌道引入1/4個交點年(即87 d左右)的周期誤差; 地球反照輻射壓攝動還將引起衛(wèi)星徑向軌道的誤差.