路 余，杜 蘭，張中凱，周佩元

(信息工程大學(xué)導(dǎo)航與空天目標(biāo)工程學(xué)院，鄭州450001)

0 引言

我國自主發(fā)展的北斗衛(wèi)星導(dǎo)航系統(tǒng)已于2012年底實(shí)現(xiàn)了亞太地區(qū)的區(qū)域?qū)Ш蕉ㄎ环?wù)［1］。導(dǎo)航星座由中高度圓軌道(MEO)、傾斜地球同步軌道(IGSO)和地球靜止軌道(GEO)三類衛(wèi)星組成。衛(wèi)星的廣播星歷提供了導(dǎo)航服務(wù)的空間基準(zhǔn)，因此是影響衛(wèi)星導(dǎo)航定位精度和可靠性的重要因素。目前廣播星歷參數(shù)有兩大類，一是基于軌道根數(shù)的GPS型廣播星歷［2］，另一類是基于衛(wèi)星狀態(tài)向量的GLONASS型廣播星歷［3］。后者由于模型參數(shù)少，參數(shù)的有效期較短，因而更新更加頻繁。

北斗混合星座的廣播星歷采用類似于GPS的標(biāo)準(zhǔn)廣播星歷格式，其中GEO衛(wèi)星采用了5度傾角旋轉(zhuǎn)法以避免軌道根數(shù)的奇點(diǎn)問題［4］。但是，用戶計(jì)算GEO衛(wèi)星位置還需要反向旋轉(zhuǎn)以恢復(fù)正確的軌道定向［5］。如此以來，破壞了用戶使用北斗混合星座廣播星歷的算法一致性。

本文提出了一種增強(qiáng)型軌道列表星歷參數(shù)模型。首先，仍采用基于衛(wèi)星狀態(tài)向量的參數(shù)表示法，保留了GLONASS型廣播星歷的無奇點(diǎn)特性;其次，改進(jìn)了衛(wèi)星的攝動(dòng)力模型，星歷有效期從原來的30分鐘擴(kuò)展到數(shù)小時(shí)，同時(shí)星歷參數(shù)也從原來的9個(gè)增加為15個(gè)。與目前采用的GPS型廣播星歷相比，新模型的參數(shù)個(gè)數(shù)、參數(shù)有效期和更新頻率均保持一致，且能夠同時(shí)適用于MEO、IGSO、GEO三類衛(wèi)星軌道的短期運(yùn)動(dòng)特性描述，無需對(duì)GEO衛(wèi)星進(jìn)行人為旋轉(zhuǎn)。對(duì)北斗混合星座的星歷參數(shù)擬合實(shí)驗(yàn)表明，三類衛(wèi)星的2小時(shí)擬合精度在非地影期均優(yōu)于1 cm。

1 GLONASS標(biāo)準(zhǔn)廣播星歷模型

GLONASS型廣播星歷是通過對(duì)地心地固坐標(biāo)系(ECEF)下的力模型進(jìn)行數(shù)值積分來計(jì)算衛(wèi)星位置的［6］。由于星歷更新較快，所需的軌道外推時(shí)間短，攝動(dòng)力僅考慮了對(duì)軌道影響較大的攝動(dòng)主項(xiàng)。

令衛(wèi)星的狀態(tài)向量為ECEF下的位置和速度向量，即XT=［RT˙RT］，則狀態(tài)微分方程為

式中:牛頓加速度主項(xiàng)a0、aJ2和a3b分別為ECEF下的地球中心引力加速度、地球扁率攝動(dòng)加速度和日月引力攝動(dòng)加速度向量;隨后的兩項(xiàng)是由于在非慣性系下描述二階時(shí)間導(dǎo)數(shù)而引進(jìn)的牽連加速度主項(xiàng)，分別為離心加速度和科里奧利加速度向量，其中示平均地球自轉(zhuǎn)角速度。

1. 1 9參數(shù)星歷模型

在GLONASS標(biāo)準(zhǔn)廣播星歷模型中，由于廣播星歷作用時(shí)間短(每30分鐘更新)，日月引力攝動(dòng)加速度a3b不再用其具體力模型計(jì)算，而是簡(jiǎn)單地令其在有效期內(nèi)保持為待求解的常向量。此時(shí)，擴(kuò)展的9參數(shù)狀態(tài)向量為XT=［RT˙RTaT3b］，且?guī)в谐踔禇l件的擴(kuò)展?fàn)顟B(tài)微分方程為

若已知?dú)v元t0時(shí)的星歷參數(shù)X0，利用積分器進(jìn)行數(shù)值積分，即可得到有效期內(nèi)任意時(shí)刻的ECEF下的衛(wèi)星位置和速度向量。為了和GLONASS星歷擬合方法一致，文中同樣采用的是Runge-Kutta四階積分器，文獻(xiàn)［7］分析了積分器誤差對(duì)GLONASS廣播星歷擬合的影響，可以通過縮短積分步長提高擬合精度。

1. 2 星歷擬合誤差

由擬合時(shí)段內(nèi)精密外推軌道的位置向量序列［Ri］(i=1，2，…n)，通過最小二乘平差方法可以擬合出式(2)中的星歷參數(shù)X0的最佳估值。具體解算流程與軌道確定相同，參見文獻(xiàn)［8］。

采用GLONASS標(biāo)準(zhǔn)9參數(shù)星歷模型，對(duì)北斗混合星座的三類衛(wèi)星進(jìn)行了星歷擬合仿真實(shí)驗(yàn)。擬合數(shù)據(jù)［Ri］由復(fù)雜力模型和 Runge-Kutta7(8)階積分器計(jì)算得到，并作為軌道位置真值。擬合時(shí)段分別為30分鐘和1小時(shí)，數(shù)據(jù)間隔均取300秒。擬合參數(shù)的初值設(shè)定中，位置速度向量取為歷元時(shí)刻t0外推值，日月引力攝動(dòng)加速度a3b取為零向量。

對(duì)北斗星座的14顆衛(wèi)星進(jìn)行了大量擬合實(shí)驗(yàn)，圖1給出了三類衛(wèi)星的一組典型的擬合誤差分布圖。可以看出:

(1)標(biāo)準(zhǔn)30分鐘擬合條件下，三類軌道都能達(dá)到2 cm的擬合精度，說明GLONASS的9參數(shù)模型可以直接用于北斗星座中高軌衛(wèi)星的30分鐘超短期星歷。

(2)當(dāng)擴(kuò)展到1小時(shí)時(shí)段的星歷擬合時(shí)，位置三分量的擬合誤差迅速增大了一個(gè)數(shù)量級(jí)(達(dá)到0.2 m)，遠(yuǎn)遠(yuǎn)超出擬合精度要求。

顯然，對(duì)于數(shù)小時(shí)擬合時(shí)段，由于考慮的力模型過于簡(jiǎn)單，擬合誤差還將進(jìn)一步增加。此時(shí)，改進(jìn)力模型將是有效方法。

圖1 M2/G4/I2衛(wèi)星30 min和1 h的擬合誤差Fig.1 Fitting errors of M2/G4/I2 satellites for 30 min and 1 h

2 增強(qiáng)型軌道列表星歷模型

目前北斗衛(wèi)星廣播星歷的有效期為2小時(shí)，為此首先分析了日月引力攝動(dòng)和太陽光壓攝動(dòng)在數(shù)小時(shí)內(nèi)的變化特性，并據(jù)此建立了針對(duì)日月引力攝動(dòng)的二次多項(xiàng)式逼近模型，從而提出了適用期更長的15參數(shù)星歷模型。

2. 1 日月引力攝動(dòng)及太陽光壓攝動(dòng)

日月引力攝動(dòng)和光壓是中高軌衛(wèi)星的兩個(gè)主要攝動(dòng)力。通常力模型是在地心慣性系(ECI)中建立，因此考察兩者在地心地固坐標(biāo)系下的短期波動(dòng)特性，還需要逐點(diǎn)進(jìn)行地球自轉(zhuǎn)變換。

在地心慣性系下，日月引力攝動(dòng)加速度為［9］:

式中:G為引力常數(shù)，Mi為攝動(dòng)天體的質(zhì)量，ri為攝動(dòng)天體相對(duì)于地心的位置矢量，ri為攝動(dòng)天體到地心的距離，Δi為衛(wèi)星相對(duì)于攝動(dòng)天體的位置矢量，Δi衛(wèi)星到攝動(dòng)天體的距離。

采用柱形地影模型的太陽光壓標(biāo)準(zhǔn)模型為［10］:

式中:ν=0時(shí)衛(wèi)星在地影期，ν=1時(shí)衛(wèi)星在非地影期。P⊙為衛(wèi)星所在區(qū)域的太陽輻射流量，CR為光壓系數(shù)，A為衛(wèi)星有效截面積，m為衛(wèi)星質(zhì)量，r⊙為太陽相對(duì)于地心位置，AU為一個(gè)天文單位，其值約為1.496×108km。

需要注意的是，衛(wèi)星進(jìn)出地影使得力模型出現(xiàn)階躍。地影是指太陽位于地球衛(wèi)星軌道面附近且日、地、衛(wèi)星近似在一條直線上時(shí)，受地球遮擋，衛(wèi)星不能受到太陽光線照射并感知太陽位置的現(xiàn)象［11］。地影會(huì)造成太陽光壓的間斷，而GEO、IGSO與MEO衛(wèi)星在入地影頻度、持續(xù)時(shí)間、最長日期變化等方面均有所不同，在地影期間衛(wèi)星受到的影響程度也不同。其中，GEO軌道的地影每年出現(xiàn)兩次，分別在春分日和秋分日附近，每次持續(xù)46天;在地影季期間，每天衛(wèi)星進(jìn)入地影的時(shí)段不同，春/秋分日最長，約為72分鐘。

2. 2 實(shí)驗(yàn)分析

考察中高軌衛(wèi)星的兩攝動(dòng)力的量級(jí)及其在ECEF坐標(biāo)系中數(shù)小時(shí)內(nèi)的變化特性。仿真時(shí)間選在春分日(3月21日)，此時(shí)太陽位于赤道面，光壓對(duì)GEO的攝動(dòng)影響最大。令MEO和IGSO衛(wèi)星軌道的赤道交點(diǎn)經(jīng)度與GEO衛(wèi)星定點(diǎn)經(jīng)度相同，光壓系數(shù)為 1.0，衛(wèi)星的有效面質(zhì)比為0.0003 m2/kg。

圖2給出了非地影期的4 h內(nèi)某MEO衛(wèi)星的日月引力攝動(dòng)、太陽光壓攝動(dòng)的加速度分量。首先，兩種攝動(dòng)力在30 min內(nèi)變化較小，可以近似用平均值表示，但隨著時(shí)間增加到2～4 h，用常值進(jìn)行逼近的誤差過大。其次，光壓攝動(dòng)的大小與星體參數(shù)即有效面質(zhì)比成正比，且比日月引力攝動(dòng)小2個(gè)數(shù)量級(jí)，可以忽略或由其它擬合參數(shù)吸收其影響，但是其造成的力模型間斷在短時(shí)段內(nèi)可能影響參數(shù)擬合精度。

對(duì)GEO和IGSO衛(wèi)星的仿真有類似結(jié)論。因此，隨著擬合時(shí)間的增加，太陽光壓和日月攝動(dòng)加速度可以采用多項(xiàng)式和(或)三角函數(shù)等進(jìn)行數(shù)學(xué)逼近表示。

2. 3 15參數(shù)星歷模型

為了控制星歷參數(shù)的個(gè)數(shù)，在數(shù)小時(shí)的擬合時(shí)段內(nèi)，選擇采用二次多項(xiàng)式表示日月引力攝動(dòng)加速度，即

則有擴(kuò)展的15參數(shù)狀態(tài)向量為

且?guī)в谐踔禇l件的擴(kuò)展?fàn)顟B(tài)微分方程為

圖2 MEO衛(wèi)星的太陽光壓(a)和日月攝動(dòng)的短期波動(dòng)變化(b)Fig.2 Acceleration change of the solar radiation pressure(a)and lunisolar gravitational perturbation for MEO satellite(b)

式中:參數(shù)初值X0的設(shè)定為［R0˙R0］取擬合時(shí)段中間點(diǎn)的位置和速度，力模型參數(shù)［a0a1a2］的初值全部取為零向量。

同理，若由擬合時(shí)段內(nèi)精密外推軌道的位置向量序列［Ri］，擬合出星歷參數(shù)X0的最佳估值，則可用Runge-Kutta四階積分器對(duì)式(6)積分得到有效期內(nèi)任意時(shí)刻ECEF下的衛(wèi)星位置和速度向量。

3 擬合實(shí)驗(yàn)和分析

將15參數(shù)的增強(qiáng)型軌道列表星歷模型用于目前北斗混合星座中14顆衛(wèi)星的星歷擬合(5個(gè)GEO衛(wèi)星G1～G5、5個(gè)IGSO衛(wèi)星I1～I(xiàn)5和4個(gè)MEO衛(wèi)星M1～M4)。首先分別對(duì)GEO/IGSO/MEO三類衛(wèi)星在地心慣性系(ECI)下進(jìn)行軌道外推，采用力模型包括15×15階的地球引力場(chǎng)、日月引力攝動(dòng)、太陽光壓、海洋潮汐及地球固體潮汐等，并采用IERS發(fā)布的地球定向參數(shù)將外推軌道轉(zhuǎn)換到ECEF坐標(biāo)系下。仿真數(shù)據(jù)從2013年3月23日開始共計(jì)11天，數(shù)據(jù)點(diǎn)間隔為300 s;按照擬合時(shí)段和每小時(shí)更新頻率進(jìn)行分組，每顆衛(wèi)星均包括264組星歷擬合數(shù)據(jù)，歷元時(shí)刻取為每組的中間時(shí)刻。需要指出的是，該實(shí)驗(yàn)時(shí)段內(nèi)，5個(gè)GEO衛(wèi)星均處于地影季，每天在地影內(nèi)的時(shí)段在55～70分鐘之間變動(dòng)(見圖3)。

圖3 GEO的地影季和試驗(yàn)數(shù)據(jù)段的地影持續(xù)時(shí)間Fig.3 Daily variations of eclipse duration in the GEO satellite eclipse season centered around vernal equinox day

為考察15參數(shù)增強(qiáng)型軌道列表星歷模型的適用性，擬合時(shí)段分別選取1/2/3/4 h進(jìn)行實(shí)驗(yàn)。參數(shù)解算采用最小二乘平差，迭代終止條件為迭代前后兩次單位權(quán)方差的相對(duì)變化小于0.001。參數(shù)初值中，位置和速度向量取擬合時(shí)段中間時(shí)刻的衛(wèi)星位置和速度，其它參數(shù)均設(shè)為0;采用Runge-Kutta四階積分器，分別向前和向后數(shù)值積分。

星歷擬合的擬合精度采用用戶距離誤差(URE)的統(tǒng)計(jì)值進(jìn)行評(píng)判［12］。MEO衛(wèi)星的 URE計(jì)算公式:

式中:ΔR、ΔT和ΔN分別為衛(wèi)星徑向、沿跡和外法向位置分量的擬合誤差。需要注意的是，這里僅分析星歷參數(shù)的數(shù)學(xué)擬合精度，不包含軌道的外推誤差。

圖4為14顆衛(wèi)星1/2/3/4 h的擬合URE的RMS最大值，初步表明增強(qiáng)型軌道列表星歷模型設(shè)計(jì)能夠滿足1～3 h內(nèi)北斗混合星座的廣播星歷參數(shù)擬合。

(1)在非地影期(左圖)，三類衛(wèi)星的星歷擬合誤差隨擬合時(shí)段的增加逐漸增大。擬合URE的RMS，在1～2 h擬合時(shí)段內(nèi)優(yōu)于1 cm，3 h擬合時(shí)段優(yōu)于5 cm，但是4 h擬合的最大RMS已接近0.2 m。

(2)對(duì)比三類衛(wèi)星可以看出，IGSO衛(wèi)星的擬合誤差最小，GEO衛(wèi)星的擬合誤差次之，MEO衛(wèi)星的擬合誤差最大?？赡苡捎谠诘毓滔迪拢琈EO衛(wèi)星的軌道較低，其速度變化較快，擬合精度較差。

(3)在地影期(右圖)，5個(gè)GEO衛(wèi)星廣播星歷的在2 h擬合時(shí)段內(nèi)的 URE的 RMS最大值為1 cm，3 h擬合時(shí)段達(dá)到8 cm，4 h擬合時(shí)段的不能滿足擬合要求。

圖4 14個(gè)衛(wèi)星擬合URE的RMS的最大值Fig.4 Max RMSvalue of 14 satellites fitting URE

圖5 是G4衛(wèi)星分別采用1/2/3/4小時(shí)的擬合URE的RMS分布圖。每當(dāng)衛(wèi)星進(jìn)出地影時(shí)，太陽光壓導(dǎo)致衛(wèi)星受力模型出現(xiàn)跳變，影響到外推軌道的局部平滑性［13］，進(jìn)而導(dǎo)致星歷擬合的較大誤差。由于每天的地影期持續(xù)1小時(shí)左右，因此僅能影響其附近的幾組數(shù)據(jù)，導(dǎo)致擬合誤差分群，跳變間隔約為1天，圖中共有5組跳變現(xiàn)象(注意到，圖5中左上1 h擬合時(shí)，若地影數(shù)據(jù)本身恰好組成一組，則不影響擬合精度，因此圖中的跳變并非1天)。此外，隨著擬合時(shí)段的增加，受地影數(shù)據(jù)影響的比例逐漸減小，擬合精度趨于平穩(wěn)，但是擬合精度自身隨時(shí)段延長衰減更快。因此，控制地影對(duì)星歷擬合的影響，應(yīng)在地影期間適當(dāng)縮短擬合時(shí)段。

圖6給出了G4星分別利用1/2/3/4 h擬合的星歷參數(shù)積分計(jì)算的ECEF速度精度。隨著擬合時(shí)長的增加，R和T方向的速度誤差快速增大，且幅度較為一致，N方向速度誤差最小，且對(duì)擬合時(shí)長不敏感。其中，2 h擬合結(jié)果計(jì)算的衛(wèi)星速度在R、T、N方向的誤差均優(yōu)于0.005 mm/s，即使在地影期也能控制在 0.05 mm/s。

4 結(jié)束語

本文分析了北斗混合星座的中高軌衛(wèi)星在數(shù)小時(shí)內(nèi)的攝動(dòng)運(yùn)動(dòng)變化特性，提出了基于擴(kuò)展?fàn)顟B(tài)向量的15參數(shù)增強(qiáng)型軌道列表星歷模型。通過對(duì)GEO、IGSO、MEO三類衛(wèi)星的星歷擬合實(shí)驗(yàn)表明，增強(qiáng)型廣播星歷能夠滿足1～3 h擬合時(shí)段的擬合URE精度要求。與現(xiàn)有的基于軌道根數(shù)的GPS型星歷模型相比，避免了奇點(diǎn)問題，數(shù)值積分可方便地同時(shí)得到衛(wèi)星位置和速度信息，因此可考慮在廣播星歷或者星基增強(qiáng)系統(tǒng)中采用和備份發(fā)播。

圖5 G4衛(wèi)星擬合URE的RMS分布Fig.5 RMSdistribution of the G4 satellite fitting URE

圖6 G4衛(wèi)星速度擬合誤差RMSFig.6 RMSerrors of the G4 satellite fitting velocity

［1］ Shi C，Zhao Q L，Li M，et al.Precise orbit determination of Beidou satellites with precise positioning［J］.Sci China Earth Sci，2012，55(7):1079 －1086.

［2］ 崔先強(qiáng)，焦文海，賈小林，等.GPS廣播星歷參數(shù)擬合算法［J］.測(cè)繪學(xué)院學(xué)報(bào)，2004，21(4):244－246.［Cui Xianqiang，Jiao Wen-hai，Jia Xiao-lin，et al.The fitting algorithm of GPSbroadcast ephemeris parameters［J］.Journal of Institute of Surveying and Mapping，2004，21(4):244 －246.］

［3］ 郭際明，孟祥廣，李宗華，等.GLONASS衛(wèi)星廣播星歷精度分析［J］.大地測(cè)量與地球動(dòng)力學(xué)，2011，31(1):68－71.［Guo Ji-ming，Meng Xiang-guang，Li Zong-hua，et al.Accuracy analysis of GLONASSsatellites broadcast ephemeris［J］.Journal of Geodesy and Geodynamics，2011，31(1):68 －71.］

［4］ 何峰，王剛，劉利，等.地球靜止軌道衛(wèi)星廣播星歷參數(shù)擬合與試驗(yàn)分析［J］.測(cè)繪學(xué)報(bào)，2011，40(增刊):52-58.［He Feng，Wang Gang，Liu Li，et al.The ephemeris fitting and experiments analysis of GEO satellite［J］.Acta Geodaetica et Cartographica Sinica，2011，40(Supplement):52 －58.］

［5］ 崔先強(qiáng)，楊元喜，吳顯兵.軌道面旋轉(zhuǎn)角對(duì)GEO衛(wèi)星廣播星歷參數(shù)擬合的影響［J］.宇航學(xué)報(bào)，2012，33(5):590－596.［Cui Xian-qiang，Yang Yuan-xi，Wu Xian-bing.Influence of the orbital plane rotation angle on GEO satellite broadcast ephemeris parameter fitting［J］.Journal of Astronautics，2012，33(5):590 －596.］

［6］ 葛奎，王解先.GLONASS衛(wèi)星位置計(jì)算與程序?qū)崿F(xiàn)［J］.測(cè)繪與空間地理信息，2009，32(2):19－22.［Ge Kui，Wang Jie-xian.Calculation of GLONASS satellite station and the realization of program［J］.Geomatics＆ Spatial Information Technology，2009，32(2):19 －22.］

［7］ 賈小林，焦文海，吳顯兵.GLONASS廣播星歷用戶算法精度分析［J］.測(cè)繪學(xué)院學(xué)報(bào)，2005，22(3):160－162.［Jia Xiaolin，Jiao Wen-hai，Wu Xian-bing.The accuracy analysis of the user arithmetic GLONASS broadcast ephemeris［J］.Journal of Institute of Surveying and Mapping，2005，22(3):160 －162.］

［8］ 杜蘭.GEO衛(wèi)星精密定軌技術(shù)研究［D］.鄭州:信息工程大學(xué)，2006.［Du Lan.Research on GEO satellite precision orbit determination［D］. Zhengzhou:Information Engineering University，2006.］

［9］ 李濟(jì)生.人造衛(wèi)星精密軌道確定［M］.北京:解放軍出版社，1995.

［10］ Montenbruck O，Gill E.Satellite orbits:models，methods and applications［M］.Berlin:Springer，2000.

［11］ 張中凱，杜蘭，王若璞，等.混合星座衛(wèi)星軌道陰影期分析［C］.第四屆中國衛(wèi)星導(dǎo)航學(xué)術(shù)年會(huì)，武漢，2013年5月15－ 17 日.［Zhang Zhong-kai，Du Lan，Wang Ruo-pu，et al.Analysis on eclipses of satellite orbits in mixed constellation［C］.China Satellite Navigation Conference，Wuhan，May 15 － 17，2013.］

［12］ Zhou S S， Cao Y L，Zhou J H，et al.Positioning accuracy assessment for the4GEO/5IGSO/2MEOconstellation of Compass［J］.Science China Physics，Mechanics and Astronomy，2012，55(12):2290－2299.

［13］ 杜蘭，張中凱，劉利，等.GEO衛(wèi)星廣播星歷的擬合參數(shù)設(shè)計(jì)［J］.中國空間科學(xué)技術(shù)，2013，33(3):46－51.［Du Lan，Zhang Zhong-kai，Liu Li，et al.Fitting parameters design for GEO broadcast ephemeris［J］.Chinese Space Science and Technology，2013，33(3):46 －51.］

［14］ 譚述森.北斗衛(wèi)星導(dǎo)航系統(tǒng)的發(fā)展與思考［J］.宇航學(xué)報(bào)，2008，29(2):391 －396.［Tan Shu-sen.Development and thought of Compass navigation satellite system ［J］.Journal of Astronautics，2008，29(2):391 －396.］

［15］ Mullins L D.Calculating satellite umbra/penumbra entry and exit positions and times［J］.Journal of the Astronautical Sciences，1991，39(4):411－422.

［16］ Zhou SS，Hu X G，Wu B，et al.Orbit determination and time synchronization for a GEO/IGSO satellite navigation constellation with regional tracking network ［J］.Science China Physics，Mechanics and Astronomy，2011，54(6):1089－1097.