2. 北京环球信息应用开发中心，北京 100094
2. Beijing Global Information Application and Development Center,Beijing 100094,China
1 引 言
为了提供导航、定位和定时服务，精密的预报轨道和星钟钟差序列需要参数化即拟合出星历参数和钟参数并形成导航电文发播给用户[1, 2, 3]。在目前的几种GNSS星座中，北斗二期采用了GEO/IGSO/MEO的混合星座，因此从方便用户的角度出发，应尽量发播一致的星历和星钟格式。由于钟参数统一采用低阶多项式模型，更多的设计工作则集中在星历模型上[4, 5]。
北斗采用类GPS的广播星历形式。这种星历能够描述3类轨道的短期变化特性，但是GEO卫星的直接星历拟合存在奇点问题，影响拟合精度甚至导致拟合失败[6]。目前GEO卫星的星历拟合采用人为旋转5°倾角的策略，即拟合得到的是接近5°倾角的IGSO卫星结果，一旦拟合出这种“伪GEO”的星历参数并发播给用户，用户计算卫星位置时还需要反向旋转5°倾角以恢复GEO轨道的正确定向[7, 8, 9, 10]。因此从完整流程来看，GEO的星历模型与IGSO和MEO并非完全一致。此外，研究表明，这种伪GEO的星历拟合仍具有一定程度的奇点问题，且拟合时段不能少于4小时[11]。
针对小偏心率和小倾角的GEO卫星轨道，通常有两种专用的轨道根数描述方法，即同步根数和第二类无奇点根数[12, 13]。文献[14]给出了基于同步根数的GEO卫星星历拟合参数设计，由于参数设计中顾及了卫星星下点东西漂移律的影响，拟合算法的稳定性和拟合精度均能够得到大幅提高。但是，同步根数通常描述GEO在地固系下的运动特性，且设计的17参数星历与现行的类GPS的广播星历个数不同，因而影响了其通用性。
本文提出一组针对GEO卫星的16参数广播星历。首先，星历参数基于第二类无奇点根数，能够同时消除小偏心率和小倾角引起的数学奇异；其次，参数拟合的数据时段可缩短至2 h，从而减弱星下点东西漂移对半长轴的长期累积影响。采用定点在84°E附近的某GEO卫星的模拟星历数据，通过星历拟合和卫星位置计算验证了新的广播星历模型，分析了包含部分地影期的10 d共计258组星历参数。非地影期的拟合结果表明，2 h的拟合用户距离误差(URE)精度优于0.05 m。
2 新的GEO广播星历卫星在惯性系下的轨道运动描述常采用3种轨道根数，即开普勒根数、第一和第二类无奇点根数[15]。其中，第二类无奇点根数是专用于描述小偏心率和小倾角轨道及其轨道变化的无奇点根数，并已长期用于GEO卫星的轨道描述，如计算GEO定点机动的控制参数等[12]。因此，第二类无奇点根数同样适用于GEO卫星广播星历的参数设计。
2.1 第二类无奇点根数经典开普勒根数(a,e,i,Ω,ω,M)分别表示轨道半长轴、偏心率、轨道倾角、升交点赤经、近点角距和平近点角。
令第二类无奇点根数表示为(a,e,i,M*)，其与经典开普勒根数的关系为[15]
式中，ω=Ω+ω；e和i分别称为偏心率向量和倾角向量；ex、ey和ix、iy分别称为偏心率向量的二维分量和倾角向量的二维分量；M*为平经度，从春分点开始计量。
与3类近点角类似，可分别定义与平经度M*相对应的真经度f*和偏经度E*
式中，f和E分别为真近点角和偏近点角。 2.2 星历参数设计
在各种摄动力作用下，轨道的变化对应了轨道根数的摄动。通常轨道根数的变化依据时变特性分为3类：长期项、长周期项和短周期项。经典的轨道分析能够给出轨道根数的一阶解析解和部分二阶解析解[15]，但是其公式过于复杂，并不适用于描述数小时内的轨道变化。
GPS标准星历利用16个参数对短期轨道变化进行参数化表达(见表 1)。针对数小时内轨道根数的3类变化，其处理方法为：①将长期项和长周期项合并为长期项，且保留为根数的变率形式；②将根数的各种短周期项影响分别投影到径向、沿迹向和轨道外法向上(RTN)，且保留综合后的主项。如此，既可以用有限数量的参数表征主要的轨道变化，同时用户仍能够采用简便的解析公式计算卫星的受摄位置向量[16]。
| 类GPS标准星历 | GEO设计星历 | |
| 参考时刻 | toe | toe |
| 轨道根数 | ( ,e,i0,Ω0e,ω0,M0) | (,ex,ey,ix0,iy0,M0*)
|
| 长期改正参数 | (Δn,idot, ) | (Δn,ixdot,iydot) |
| 短周期改正参数 | (Crc,Crs,Cuc,Cus,Cic,Cis) | (Crc,Crs,Cλc,Cλs,CNc,CNs) |
| 拟合URE精度 | ＜0.1 m RMS (root-mean-square) | ＜0.05 m RMS |
| 拟合时段 | ≥4 h | ≥2 h |
| 有效时段 | 2 h | 2 h |
在设计GEO星历参数时，针对第二类无奇点根数表示法进行相应修改，可以得到16参数的GEO广播星历表示法，即
为轨道半长轴的均方根；ex、ey为偏心率向量的二维分量；ix0、iy0为倾角向量的二维分量；M0*为平经度；Δn为平均角速度与计算值的偏差；ixdot、iydot为倾角向量的二维变率；Crc、Crs为地心距的短周期改正项系数；Cλc、Cλs为真经度的短周期改正项系数；CNc、CNs为法向距离的短周期改正项系数。
表 1列出了新设计的GEO星历与现用的类GPS标准星历的对比关系。两者最大的区别即为轨道根数的选取不同，前者采用第二类无奇点根数，后者为开普勒根数。因此，在设计参数中：①轨道面在地心惯性系(ECI)下的空间定向及其变率，取代常规的倾角和升交点赤经(i,Ω)及其变率，采用倾角向量i=[ix iy]T及其变率表征；②沿迹方向的轨道运动快变量，则是取代了幅角u=ω+f，而采用真经度f*=Ω+u表征。
2.3 GEO卫星的星历计算在轨道上新建一个拟春分点，使得拟春分点至升交点的角距为Ω。定义轨道坐标系O-PQW：原点O为地球质心，OP轴指向拟春分点方向，OW轴沿轨道角动量方向，OQ轴由右手法则确定。
令θg0为历元时刻toe的Greenwich平恒星时，由历元时刻的星历参数可以计算t时刻卫星的地心地固坐标系(ECEF)下的位置向量。其流程如下：
(1) 计算半长轴和平均运动速度
式中，μ是地球引力常数。
(2) 计算t时刻的平经度、偏经度和真经度
由广义开普勒方程迭代计算偏经度 以及偏经度E*和真经度f0*的关系式 式中
(3) 计算包含短周期项改正的RTN三分量
(4) 计算卫星在轨道坐标系下的位置向量
(5)计算轨道坐标系O-PQW至惯性系(ECI)的旋转矩阵
式中
此外，地心惯性系至地心地固系的关系为旋转矩阵Rz(θg)，且
式中，ωe是地球自转角速度。
(6) 计算卫星在ECEF下的位置向量
3 拟合试验和分析
令GEO卫星定点在84°E附近，生成星历数据的力模型为10阶次地球引力场、日月引力摄动、光压、地球固体潮和海潮，采用IERS发布的地球定向参数。由于GEO卫星有频繁的定点机动控制，这里仅给出10 d的地固系位置向量序列。此外，GEO卫星每年有2次地影季，分别在春分日和秋分日附近，每次持续46 d[13]；为了同时分析GEO卫星在地影季的星历拟合，仿真从2012-04-06-12:00开始，卫星在前7 d处于地影季，期间每天有约1 h弧段位于地影内(称为地影期)。当拟合时段同时包含有进出地影数据时，由于光压的间断性影响了轨道数据的平滑性，会降低拟合效果。
为适应每小时的星历更新模式，将10 d数据逐小时分割为258组，每组分别包含2/3/4 h的时段，时段内数据点间隔300 s。星历参数解算采用最小二乘拟合，当前后两次拟合的残差均方差(RMS)的相对偏差小于0.001时，迭代结束[17, 18]。
为方便参数迭代初值的选取，将星历参数中的(，M0*)转换为两个新的拟合参数，即(Δa，Δλe)，分别表示轨道半长轴与参考值aref的偏差、星下点经度与定点经度λe的偏差，其中半长轴的参考值取为42 165 760 m[12]，定点经度取为84°E。迭代结束后，可由拟合参数反算相应的星历参数，即
此时，第一组拟合参数的初值可全部取为0，后续组的参数初值则分别采用其前一组的拟合结果。这里，星历数据将作为轨道位置真值，用于评定拟合星历和拟合URE的精度。
2/3/4 h的星历拟合的成功率均为100%，除首组拟合迭代次数为5，其余组均稳定为4次迭代。
3.1 轨道拟合精度令卫星径向、沿迹和外法向位置分量误差分别为ΔR、ΔT和ΔN，图 1给出了3种拟合时段条件下拟合残差的RMS统计结果。
|
| 图 1 拟合时长对拟合残差RMS的影响 Fig. 1 RMS of fit residuals vs. the length of fit arcs |
(1) 在常规的非地影期(图 1(a))，N分量的拟合误差RMS的平均值最小，且几乎不受拟合时长影响，保持在毫米级精度；T和R位置分量的拟合误差随拟合时段的增加快速增大，4 h拟合的分量误差分别达到0.7 m和0.2 m，因此，拟合时段应控制在2～3 h为宜。
(2) 在特殊的地影期(图 1(b))，三分量的拟合误差RMS的最大值随拟合时段的变化与非地影期较一致。但是，进出地影导致光压间断，不连续的力模型影响到星历数据的平滑性，从而引起较大的拟合误差。3 h拟合的T和R位置分量误差最大值分别超过0.5 m和0.1 m，此时，拟合时段应进一步控制在2 h以内。
3.2 拟合参数的轨道表征图 2列出了2 h拟合的部分轨道参数的时间序列。图 2(a)中的偏心率和轨道倾角由相关拟合参数计算得到。令T为卫星轨道周期，其中轨道倾角序列变化平稳，而偏心率由于没有设计相应的半轨道周期T/2的短周期项改正参数，其序列表现出T/2(近12 h)短周期项改正参数，其序列表现出T/2的主周期变化。
|
| 图 2 2 h拟合得到的部分轨道参数时序 Fig. 2 Time series of some orbital elements from 2 h fit |
图 2(b)分别列出了半长轴、平均角速度和定点经度与其参考值或计算值的拟合偏差序列。由于半长轴的T/2主周期项已由向径方向的短周期参数吸收，半长轴偏差Δa时序表现出T/3(近8 h)的次周期变化[11]，且振幅达到500 m。根据开普勒第三定律，Δn与Δa具有线性负相关性，两者的时序变化也呈现出该规律性。此外，注意到定点经度在10 d内向西漂移了约0.1°，同时轨道半长轴表现了一致的增大趋势，即轨道逐渐升高。事实上，在16参数的星历设计中，半长轴仅有短周期摄动变化，并未考虑其长期变化[14, 19]，因此，缩短拟合时段可以控制半长轴长期变化的累积效应，改善轨道的短期变化表征能力，这也是图 1中2 h拟合优于3 h和4 h拟合的主要原因。
图 2(c)给出了二维倾角向量的变率序列。注意到前7 d均有1次跳跃现象，且相隔1 d，这是进出地影的数据对相关时段拟合结果的影响。需要指出的是，卫星地影导致的光压间断，原本应该主要影响轨道面内的拟合参数及其变率，但是由于摄动参数个数限制，不恰当地反映在轨道面的空间定向元素变率中，这也表明现有的星历参数选取仍有改进的余地。
3.3 拟合URE精度URE是计算星历和钟差误差对定位影响的一个重要评价指标[1, 20]。URE的来源主要是外推的轨道和钟差误差，此外还包括星历参数拟合和钟参数拟合的拟合误差。这里仅分析星历参数拟合的单项影响，故称为拟合URE。计算公式为
URE=(0.96ΔR2+0.04ΔT2+0.04ΔN2)1/2
显然，径向拟合误差ΔR对URE贡献最大，其他两方向分量误差的影响因子要小得多。为便于比较，对于同组的 2/3/4 h拟合结果，URE的统计计算均采用时段中间2h，即同一参考历元前后各1 h的拟合残差。
表 2为3种拟合时段条件下拟合URE的RMS统计。显然，拟合时段越短，拟合URE精度越高；2/3/4 h拟合URE的RMS均值分别优于1 cm、5 cm和0.1 m。此外，包含进/出地影时段的拟合URE的精度有明显降低。与非地影期的最大RMS值相比较，地影期2/3/4 h拟合的相应结果分别放大约4倍、2倍和1.5倍。
| m | |||||
| 拟合时段 | 非地影期 | 地影期 | |||
| max | min | mean | max | ||
| 4 h | 0.169 | 0.024 | 0.099 | 0.278 | |
| 3 h | 0.063 | 0.005 | 0.036 | 0.111 | |
| 2 h | 0.015 | 0.001 | 0.009 | 0.055 | |
地影季期间，GEO卫星每天进出地影的时刻变化不大，且在影持续时间最长可达72 min[19]，因此有可能影响到相邻的1～3组拟合时段。图 3是2 h拟合URE的RMS时序。在拟合结果的前7 d，一致出现了1～2次较大拟合偏差，且相隔1 d。因此，对于地影期星历拟合，应适当减小拟合时长以满足精度需求。
|
| 图 3 2 h拟合URE的RMS时序 Fig. 3 RMS time series of the fitting URE of 2 h fit |
现有的GEO卫星广播星历，保持了与MEO/IGSO一致的发播参数格式，但是其星历参数的拟合算法和用户卫星位置计算，均需要旋转一定的轨道倾角。
针对GEO轨道特性，提出基于第二类无奇点根数的16参数GEO广播星历及用户卫星位置计算方法。与现有星历模型相比，参数个数相同，卫星的ECEF位置向量计算量相当，但是具有无奇点和拟合精度高的优点。此外，拟合初值的设置简单，拟合算法与用户算法均无需引入轨道倾角的人为旋转。
拟合试验结果表明，新的16参数设计更适合描述2 h内的定点运动。在非地影期间，2 h和3 h星历拟合的拟合URE精度分别优于0.05 m和0.1 m。此外，为降低跨地影数据的影响，也应适当缩短拟合时段。
为进一步提高星历拟合精度，并且与GPS的增强型18参数广播星历相匹配，也已设计了配套的基于第二类无奇点根数的18参数GEO广播星历及用户卫星位置计算方法将另文讨论。
| [1] | ZHOU S S, CAO Y L, ZHOU J H, et al. Positioning Accuracy Assessment for the 4GEO/5IGSO/2MEO Constellation of COMPASS[J]. Science China: Physics, Mechanics and Astronomy, 2012, 55(12): 2290-2299. |
| [2] | YANG Yuanxi, LI Jinlong, XU Junyi, et al. Contribution of the Compass Satellite Navigation System to Global PNT Users[J]. Chinese Science Bulletin, 2011, 56(26): 2813-2819. |
| [3] | MONTENBRUCK O, HAUSCHILD A, STEIGENBERGER P, et al. Initial Assessment of the COMPASS/BeiDou-2 Regional Navigation Satellite System[J]. GPS Solutions, 2013, 17(2): 211-222. |
| [4] | DIESPOSTI R, DILELLIO J, KELLEY C, et al. The Proposed State Vector Representation of Broadcast Navigation Message for User Equipment Implementation of GPS Satellite Ephemeris Propagation[C]//Proceedings of Institute of Navigation National Technology Meeting 2004. San Diego: ION, 2004: 92-134. |
| [5] | FU X F, WU M P. Optimal Design of Broadcast Ephemeris Parameters for a Navigation Satellite System [J]. GPS Solutions. 2012, 16(4): 439-448. |
| [6] | HAN Xingyuan, XIANG Kaiheng, WANG Haihong. Research on Broadcast Ephemeris Parameters Fitting Algorithm Based on the First Class of No Singularity Variables[J]. Spacecraft Engineering, 2011, 20(4): 54-59. (韩星远,向开恒,王海红. 第一类无奇点变量的广播星历参数拟合算法[J].航天器工程, 2011, 20(4): 54-59.) |
| [7] | HE Feng, WANG Gang, LIU Li, et al. Ephemeris Fitting and Experiments Analysis of GEO Satellite[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(Sup.): 52-58.(何峰,王刚,刘利,等. 地球静止轨道卫星广播星历参数拟合与试验分析[J]. 测绘学报,2011, 40(增刊): 52-58.) |
| [8] | CHEN Liucheng, HAN Chunhao, CHEN Jingping. The Research of Satellites Broadcast Ephemeris Parameters Fitting Arithmetic[J]. Science of Surveying and Mapping, 2007, 32(3):12-14. (陈刘成,韩春好,陈金平. 广播星历参数拟合算法研究[J]. 测绘科学,2007, 32(3):12-14.) |
| [9] | HUANG Yong, HU Xiaogong, WANG Xiaoya, et al. Precision Analysis of Broadcast Ephemeris for Medium and High Orbit Satellite[J]. Progress in Astronomy, 2006, 24(1):81-87. (黄勇,胡小工,王小亚,等. 中高轨卫星广播星历精度分析[J]. 天文学进展,2006, 24(1): 81-87.) |
| [10] | RUAN Rengui, JIA Xiaolin, WU Xianbing, et al. Broadcast Ephemeris Parameters Fitting for GEO Satellites Based on Coordinate Transformation[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(Sup.): 145-150. (阮仁桂,贾小林,吴显兵,等. 关于坐标旋转法进行地球静止轨道导航卫星广播星历拟合的探讨[J]. 测绘学报, 2011, 40(增刊): 145-150.) |
| [11] | HUANG Hua, LIU Lin, ZHOU Jianhua, et al. Research on 18 Elements Broadcast Ephemeris Model[J]. Journal of Spacecraft TT＆C Technology, 2012, 31(3): 80-84. (黄华, 刘林, 周建华, 等. 18参数广播星历分析研究[J]. 飞行器测控学报, 2012, 31(3): 80-84.) |
| [12] | LI Hengnian. Geostationary Satellite Orbital Analysis and Collocation Strategies[M]. Beijing: National Defense Industry Press, 2010. (李恒年. 地球静止卫星轨道与共位控制技术[M]. 北京: 国防工业出版社, 2010.) |
| [13] | SOOP E M. Handbook of Geostationary Orbits[M]. Beijing: National Defense Industry Press, 1999. (SOOP E M. 地球静止轨道手册[M].北京: 国防工业出版社, 1999.) |
| [14] | DU Lan, ZHANG Zhongkai, LIU Li, et al. Design of Fitting Parameters for GEO Broadcast Ephemeris[J]. Chinese Space Science and Technology, 2013, 33(3): 46-51. (杜兰, 张中凯, 刘利, 等. GEO卫星广播星历的拟合参数设计[J]. 中国空间科学技术, 2013, 33(3): 46-51.) |
| [15] | LIU Lin, HU Songjie, WANG Xin. An Introduction of Astrodynamics[M]. Nanjing: Nanjing University Press, 2006. (刘林, 胡松杰, 王歆.航天动力学引论[M].南京: 南京大学出版社, 2006.) |
| [16] | XU G C. GPS: Theory, Algorithms, and Applications[M]. Berlin: Springer-Verlag, 2007. |
| [17] | MONTENBRUCK O, EBERHARD G. Satellite Orbits: Models, Methods, and Applications[M]. Berlin: Springer-Verlag, 2000. |
| [18] | TAPLEY B D, SCHUTZ B E, BORN G H. Statistical Orbit Determination [M]. New York: Elsevier Academic Press, 2004. |
| [19] | MULLINS L D.Calculating Satellite Umbra/Penumbra Entry and Exit Positions and Times[J]. Journal of the Astronautical Sciences, 1991, 39: 411-422. |
| [20] | GUO Rui, ZHOU Jianhua, HU Xiaogong, et al. A New Strategy of Rapid Orbit Recovery for the Geostationary Satellite[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(Sup.): 19-25. (郭睿, 周建华, 胡小工, 等. 一种地球静止轨道卫星的快速恢复定轨方法[J]. 测绘学报, 2011, 40(增刊): 19-25.) |
| [21] | KUZNETSOV E D. The Effect of the Radiation Pressure on the Orbital Evolution of Geosynchronous Objects[J]. Solar System Research, 2011, 45(5): 433-446. |