Corresponding author: WU Meiping E-mail：meipingwu@263.net
1 引 言
进入21世纪以来，卫星导航定位技术得以迅速发展。美国的GPS、俄罗斯的GLONASS、欧盟的Galileo和中国的BeiDou等作为GNSS的主要成员，在导航定位领域发挥着越来越重要的作用，就单独系统而言，GPS、Galileo、GLONASS属于全球定位系统，为保证全球对地覆盖性能，其导航卫星被均匀地布设在几个近圆的中高轨道(medium earth orbit,MEO)面内。而北斗系统二期将在2012年建成区域卫星导航系统，其空间星座包括5颗对地静止轨道(geostationary earth orbit,GEO)卫星、4颗MEO卫星和5颗倾斜地球同步轨道(inclined geosynchronous orbit,IGSO)卫星[1,2]，因此也被称为混合星座卫星导航系统。
除了确定用户位置以外，利用GNSS确定用户的速度和加速度也有重要的应用需求。例如在航空重力测量中，为确定航空载体所在空间位置的重力异常信息，需利用GNSS以mGal(1 mGal=1×10－5 m/s2)级精度确定载体的运动加速度[3,4,5]；在强实时应用领域中，如运动状态监测[6]、自动刹车系统[7]等，则需要以高采样率计算载体的运动状态信息，因此对计算的实时性有较高的要求。这些应用均需求解导航卫星的在轨速度和加速度。作为GNSS的空间基准，导航卫星速度和加速度的计算精度直接影响最终解算结果，因此需对其计算方法和误差特性进行深入地分析。
导航卫星的在轨状态可由广播星历或精密星历进行描述。广播星历事先计算，用户可实时获取，实时性强，但目前的精度约为3 m左右[8]；精密星历由IGS (International GNSS Service)等机构滞后发布，但可达优于5 cm的精度[9]。目前常见的广播星历类型包括Kepler根数型、GEO型和位置速度型。在卫星导航系统的接口控制文件(interface control document,ICD)中，一般仅给出导航卫星的位置、速度计算方法。本文详细推导了基于这3种广播星历的卫星加速度的计算公式，并对计算精度进行分析。
导航卫星精密星历一般以一定采样率给出卫星在地心地固坐标系(earth-centred-earth-fixed,ECEF)下的位置，其他时刻的位置通过插值得到，并通过对计算的卫星位置时间序列进行差分运算以获得卫星的速度和加速度。
文献[10]对基于高精度定位结果、原始多普勒频移观测量、载波相位差分导出的多普勒频移观测值等3种速度确定方法进行了对比分析。国内外有较多文献对航空重力测量中如何利用GPS和数字滤波技术确定载体垂直加速度的方法进行了详细的阐述[11,12,13]。但这些文献均未对导航卫星速度和加速度的计算方法做系统的总结，特别是没有对各种方法的计算精度特性和算法适用性等问题进行深入探讨。
本文首先给出几种导航卫星速度和加速度的计算方法，包括：① 基于广播星历参数的公式法；② 基于导航卫星位置序列的数值差分方法；③ 基于导航卫星位置序列的解析差分方法。对算法的计算精度和适用性进行了比较分析，并给出了相应的结论，最后通过在坐标精确已知的CORS站实测数据验证了上述结论的正确性。
2 基于广播星历的卫星速度和加速度计算目前广泛使用的3类广播星历的参数类型及示范系统如表 1所示。
| 类型 | 参数内容 | 示范系统 |
| Kepler 根数 | 对应参考历元的6个Kepler椭圆参数；3个长期项摄动参数和6个周期性摄动参数 | GPS Galileo BeiDou |
| GEO | 同Kepler型，但采用了不同的星历拟合方法以解决GEO卫星的小倾角问题 | BeiDou |
| 位置- 速度 | 对应参考历元的卫星瞬时位置、速度；日月引力摄动加速度 | GLONASS |
根据ICD文件描述的广播星历参数及对应的用户算法[14]，卫星在ECEF坐标系统下的位置可计算为
式中，rECEF、rorb分别为导航卫星在ECEF坐标系和轨道坐标系下的位置矢量(轨道坐标系定义为以卫星轨道面为参考平面，以轨道面的法线方向为Z轴，以升交点矢径方向为X轴，Y轴与XZ轴构成右手坐标系)
为从轨道坐标系到地固坐标系的转换矩阵，并且
式中，ik、Lk分别为计算时刻tk对应的轨道倾角和升交点大地经度，并有
式中，rk、uk分别为计算时刻卫星的矢径大小和经摄动改正的升交角距，均可由导航电文直接求得。R3、R1为坐标旋转矩阵，且有
ik、Lk可由下式计算
式中，Φk为升交距角，可由广播星历直接计算；ωe为地球自转角速度；Cic、Cis为导航电文给出的轨道倾角周期项改正系数；
为升交点赤经的变率；i0、
分别为参考时刻的卫星轨道倾角及变率；te、tk分别对应星历历元时刻和计算时刻。
对式(1)进行一阶微分，可得到导航卫星的速度计算公式，即
由式(5)知
联合式(6)，有
式中，旋转矩阵的导数可求为
卫星在轨道坐标系中的速度
由式(3)微分得
式中的导数项按下列公式计算
式中，fk、Ek为导航卫星轨道的真近点角和偏近点角；n0、Δn分别为卫星的平均角速度及改正；a、e分别为轨道的半长轴和偏心率；Crc、Crs、Cuc、Cus为摄动力的调和改正系数，这些参数均由广播星历得到。对式(6)再次微分，可得导航卫星加速度的计算公式
式(12)中旋转矩阵的二阶导数项
的计算需要计算
、
。这两项参数在导航星历中并未给出，但由于其量级相对较小，可忽略其影响[15]，即
、
，则旋转矩阵的二阶导数可求为
另外根据牛顿定律，
可计算为
式中，μ=3.986 004 418×1014 m3/s2为地球引力常数。
2.2 GEO型广播星历的卫星速度加速度
正如前文所述，北斗卫星导航系统星座中包括GEO、IGSO、MEO等轨道类型的导航卫星，其中GEO卫星的倾角i接近于0，将导致升交点赤经Ω和近地点角ω具有奇异性。如果直接使用经典的最小二乘法进行广播星历参数拟合[16]，就会出现拟合精度较差或拟合失败的情况。一般通过旋转参考坐标系来解决这一问题，即将Kepler星历拟合中参考的地固坐标系绕X轴旋转θ得到新的地固坐标系，并在新的坐标系中进行星历参数的拟合[17,18]。以旋转坐标系方法拟合得到的卫星广播星历称为GEO型星历。下面推导基于旋转方案的GEO星历卫星速度和加速度的计算方法。
对GEO卫星，利用式(1)计算的卫星位置是GEO卫星在旋转坐标系下的位置，需进一步进行旋转改正以获得ECEF坐标系下的位置。改正方法如下
式中，rR为GEO卫星在旋转坐标系下的位置。对式(15)进行微分运算可得GEO卫星的速度解算公式为 式中 再次微分可得加速度计算公式为 式(16)、(18)中，
、
的计算方法同Kepler根数型星历中<
和
的计算方法，并且
2.3 位置-速度型广播星历的卫星速度加速度
GLONASS采用了位置-速度型广播星历形式[19]，以每半小时为更新周期给出卫星在地固坐标系下的位置矢量、速度矢量和日月引力摄动加速度。用户采用积分方式获取即时位置和速度，积分方程可表示为
式中,矢量(x,y,z)、(Vx,Vy,Vz)和(ax,ay,az)分别表示卫星的位置、速度和加速度矢量；
为日月摄动引力加速度矢量。
由此知，位置-速度型广播星历直接给出了卫星加速度计算的解析公式，而卫星的速度计算应根据式进行数值积分获得，常用的积分算法是四阶龙格-库塔积分[20]。
3 基于数值差分法计算卫星速度和加速度在利用卫星广播星历或精密星历得到了导航卫星的位置序列后，可通过数值差分法得到卫星的速度和加速度。理想数值差分器的频率响应为[21]
式中，H表示差分器；ω、ωs分别为信号频率和采样频率；T为ωs对应的采样时间。在对离散信号进行差分处理时，常选用有限冲激响应(finite impulse response,FIR)差分滤波器，即 式中，
为输入离散信号Φ(t)对时间t的一阶导数；N为滤波器长度。再次FIR滤波可得离散信号的二阶导数
h为N－1阶脉冲响应，其定义为
在数值差分器的应用时需综合考虑信号和差分器的频率特性，实际FIR差分滤波器仅是对理想滤波器的一种近似，因此存在一个带宽限制，当信号频率位于滤波器带宽之外时，差分滤波器对信号将产生抑制作用，应用中可利用此性质抑制信号中的噪声信号。
由于导航卫星的低动态特性，在计算其速度、加速度时，低阶的差分器即可满足要求。如一阶中心差分器，滤波器系数可取为
由此得到卫星速度的计算公式 式中，Δt为卫星位置序列的时间间隔；x、
为卫星t时刻的位置、速度矢量。同理可得加速度的计算公式
式中，
为卫星加速度矢量。根据误差传播定律得到速度、加速度的计算误差为
当采样频率为1 Hz时，一阶中心差分器的幅频响应特性如图 1所示，可见频率为0~0.5 Hz的信号能通过一阶差分滤波器，但仅有频率为0.25 Hz的信号可完全无衰减通过该差分器，滤波器对其余频率信号均有抑制作用。同时由图 1所示的滤波器低频特性可知，一阶差分滤波器可完全消除信号中的常值分量和大部分低频分量，因此对慢变型误差不敏感。
|
| 图 1 一阶中心差分器幅频响应特性 Fig. 1 Frequency response of the first-order central differentor |
在基于GNSS的载体速度和加速度确定的相关数据处理中，卫星的速度、加速度值在每一次迭代运算中一般都需要重新计算，但每次迭代运算一般对应不同时刻，数值差分法效率较低，此时可采用解析差分法。该方法的思路是首先利用短时段内的卫星位置序列建立卫星的轨道模型，并通过对模型进行差分运算获得速度和加速度计算的解析公式，利用该解析公式可直接计算出该时段内任意时刻的卫星速度、加速度。例如可用Lagrange多项式对卫星位置序列进行建模，即
式中，x(tk)对应卫星tk时刻的位置分量；Xi为已知的ti时刻的卫星位置；n为多项式的阶数。为保证插值精度，一般要求tn≤tk≤t－n。将式(29)对时间分别求一阶、二阶导数可获得计算卫星速度、加速度的解析计算公式，即 和 也可采用多项式插值公式，即 则对应的速度和加速度计算公式为
值得注意的是，由于精密星历的历元间间隔较长(一般为15 min)，导致解析差分法计算的卫星速度、加速度的精度无法保证，可首先适当对位置序列进行加密以获得更为准确的解析模型。
5 计算精度比较分析为对上述3种算法求解的卫星速度、加速度进行比较分析，以精密星历并运用一阶数值差分法计算的卫星速度、加速度作为参照标准，计算各算法结果与该标准值之间的偏差，并通过在坐标已知点上的速度、加速度确定实验对各方法计算的精度进行评估。
5.1 广播星历法精度分析以北斗卫星中的GEO、IGSO和MEO卫星为例，其中GEO卫星采用GEO型广播星历，IGSO和MEO采用Kepler根数型广播星历。计算24 h的卫星速度、加速度与标准值之间的偏差分别如图 2、图 3所示。
|
| 图 2 卫星速度偏差 Fig. 2 Satellite velocity bias based on Kepler elements |
|
| 图 3 卫星加速度偏差 Fig. 3 Satellite acceleration bias based on Kepler elements |
由图 2、图 3知，3类卫星的速度计算偏差均小于1 mm/s，相比较而言，高轨道的GEO、IGSO卫星的计算精度优于MEO卫星。由图 3知，由于在式(12)转换矩阵
的计算中忽略了
、
的影响，这一简化对加速度的计算精度有较大的损失，只能达到mGal(1 mGal=10-3 cm/s2)级的计算精度，其中对MEO卫星的影响最大，计算误差达到10 mGal的量级。同时，由于在GEO星历拟合过程中，坐标系的旋转致了Lk计算的不连续，进而导致计算的加速度也呈现不连续的特点。
GLONASS星历采用位置速度序列表征导航卫星轨道，并直接给出了摄动影响较大的日月引力加速度。根据前面的分析可知，GLONASS卫星的加速度可根据解析公式直接求得，而其速度则需通过数值积分法获得，计算结果与标准值之间的偏差如图 4所示。
|
| 图 4 GLONASS卫星位置-速度型广播星历计算的计算结果 Fig. 4 Velocity and acceleration biases of position-velocity broadcast ephemeris of GLONASS |
由于GLONASS星历求取加速度时未经过数学拟合，所表达的物理含义更为明确，因此方法简单，精度较高，两个轨道周期内的计算结果表明，GLONASS星历计算的卫星加速度误差优于0.5 mGal，但由于舍弃了其他摄动力的影响，导致通过积分求取的卫星速度发散较快、精度较差，和Kepler根数型广播星历相比，精度约低一个量级，最大偏差可达约1 cm/s。
5.2 数值差分法精度分析通过一阶数值差分器的频谱特性可知，差分器对常值分量可完全滤除，对低频和高频噪声则有抑制作用，广播星历的偏差属于慢变型误差，这种慢变型误差运用数字差分时将被消除。分别用广播星历和精密星历计算出卫星的位置序列，再采用一阶中心差分器分别求取卫星的速度和加速度，GEO、IGSO和MEO卫星的计算结果偏差如图 5所示。
|
| 图 5 广播星历和精密星历的数值差分法加速度计算偏差 Fig. 5 Acceleration bias using numerical differentor between broadcast and precision ephemeris |
图 5反映了由广播星历和精密星历计算的卫星位置经过一阶中心差分计算的卫星加速度之间的差异不大，均优于0.2 mGal。计算结果进一步说明了广播星历偏差的低频特性，因此在强实时性需求的应用中，采用数值差分法时，可利用广播星历替代精密星历。
5.3 解析差分法精度分析以GPS卫星中的PRN13在24 h内的实际数据为例，解析差分计算的卫星速度加速度偏差如图 6所示。
|
| 图 6 PRN13的解析差分计算结果 Fig. 6 Velocity and acceleration bias using analytical differencing for PRN13 |
从图 6所示的计算结果可以看出解析差分计算的速度和标准值之间的偏差较大，达到1 mm/s；但计算的加速度和标准值之间的偏差不大，小于0.5 mGal。主要原因是根据短时段内的卫星位置序列所建立的轨道模型并不符合轨道的实际运行规律，这种模型偏差将对速度计算结果产生影响。
影响解析差分法计算精度的另一个原因是用于解析模型建立的数据序列的位置序列之间的时间间隔。表 2反映了不同的数据采样率与模型精度之间的关系。
| ωs/Hz | 速度偏差/(mm/s) | 加速度偏差/mGal | ||||
| vx | vy | vz | ax | ay | az | |
| 100 | 0.007 | 0.007 | 0.008 | 2.004 | 1.944 | 1.664 |
| 10 | 0.007 | 0.007 | 0.008 | 0.116 | 0.108 | 0.053 |
| 2 | 0.005 | 0.005 | 0.006 | 0.020 | 0.018 | 0.009 |
| 0.5 | 0.021 | 0.021 | 0.024 | 0.015 | 0.017 | 0.009 |
| 0.1 | 0.698 | 0.700 | 0.790 | 0.053 | 0.054 | 0.047 |
| 0.01 | 70.49 | 70.73 | 79.74 | 3.069 | 3.069 | 2.323 |
表 2中，均采用九阶的Lagrange多项式对卫星轨道进行建模，表中第1列ωs代表用于多项式建模的位置序列的采样率，均以ωs=1 Hz为参照准。结果表明过低采样率对速度和加速度的计算均会产生影响，另外当过高的采样率也会导致加速度的计算精度将会降低，从式(28)所示的误差表达式可知，加速度计算误差与采样间隔的四次方成反比，过低的采样率会放大位置序列中的高频噪声，对计算结果产生不利影响。
5.4 精度对比为对上述各方法计算的精度卫星速度和加速度进行评估，采用在坐标已知点的静止地面试验进行评估。试验选择某CORS网的两个静态站点，其接收机类型均为Leica GRX1200GGPRO高精度测地型接收机，基线距离约为20 km。试验时间为2011年5月10日0时至2时共2 h(均为GPST)。试验仅接收GPS卫星信号，卫星截止仰角设定为15°，基线双差模糊度通过BERNESE求解[22]，采用双差算法解算接收机的速度和加速度，计算结果如表 3所示。
| 算法类型 | 速度偏差/(mm/s) | 加速度偏差/mGal | ||||
| North | East | Up | North | East | Up | |
| 解析公式 | 1.11 | 1.58 | 2.52 | 88.36 | 119.25 | 201.65 |
| 数值差分 | 1.01 | 1.38 | 2.29 | 83.17 | 114.91 | 191.57 |
| 解析差分 | 2.53 | 2.09 | 4.56 | 83.32 | 115.34 | 191.96 |
| 注:表中,速度、加速度解算结果已转换到当地地理坐标系,North、East和Up即分别为该坐标系的北、东、天方向 | ||||||
从表 3可以看出，基于精密星历数值差分法解算得到的静态接收机的速度和加速度具有最高的精度，因此前文中其余方法与此方法计算结果的偏差也可在很大程度上反映这些计算方法的精度；对速度计算精度而言，广播星历解析公式也可达到mm/s的精度，因此可在实时测速应用中直接采用广播星历法，解析差分法较数值差分法的测速精度低，但加速度的计算精度和数值差分法接近。
6 结 论导航卫星自身的速度和加速度计算是用户进行速度和加速度确定时需要解决的关键问题之一，其计算精度也直接影响最终的解算精度。针对这一问题，结合前面的分析，可以得到以下结论。
(1) 通过合理的简化和数值计算，从常用的3类广播星历，即Kepler根数型、GEO型和位置-速度型均可得到卫星速度和加速度的解析计算公式，但计算精度低，难以满足高精度测定速度和加速度应用的需求。相比较而言，位置-速度型星历直接给出了加速度计算的解析公式，因此加速度计算方法简单、精度最高，但速度计算精度较差；Kepler根数型星历和GEO型星历的精度受到所忽略的二阶导数项的影响，加速度计算精度较低；另外高轨道卫星(IGSO、GEO)的计算精度优于中轨卫星(MEO)。
(2) 数值差分法具有最高的速度、加速度计算精度；解析差分法加速度计算精度接近数值差分法，但速度计算精度较差。
(3) 数值差分器对信号的低频分量具有抑制作用，因此对由广播星历计算的卫星位置序列进行数值差分获得的速度和加速度也具有和精密星历相当的计算精度。
| [1] | YANG Yuanxi. Progress, Contribution and Challenges of Compass/Beidou Satellite Navigation System[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(1):1-6.(杨元喜.北斗卫星导航系统的进展、贡献与挑战[J].测绘学报,2010, 39(1):1-6.) |
| [2] | Administrator Office of China Satellite Navigation. Interface Control Document of the Space Signal of Beidou Satellite Navigation System [EB/OL]. [2011-12-27]. (中国卫星导航系统管理办公室.北斗卫星导航系统空间信号接口控制文件[EB/OL]. [2011-12-27]. |
| [3] | JEKELI C. On Precision Kinematic Accelerations for Airborne Gravimetry[J]. Journal of Geodetic Science, 2011, 1(4): 367-378. |
| [4] | SUN Zhongmiao, ZHAI Zhenhe, LI Yingchun. Analysis of the Resolution and Accuracy of Airborne Gravity Survey[J]. Progress in Geophys, 2010, 25(3) :795-798. (孙中苗,翟振和,李迎春. 航空重力测量的分辨率和精度分析[J]. 地球物理学进展, 2010, 25(3):795-798.) |
| [5] | JEKELI C, GARCIA R. GPS Phase Acceleration for Moving-base Vector Gravimetry[J]. Journal of Geodesy, 1997, 71(10): 630-639. |
| [6] | XU Guochang. A Concept of Precise Kinematic Positioning and Flight-state Monitoring from the AGMASCO Practice[J]. Earth Planets Space, 2000, 52(10): 831-835. |
| [7] | WANG Liang. Experiment and Result Analysis of the Performance of ABS Based on VBOX Facilities [J]. Light Vehicles, 2010(5-6): 25-35. (王亮. 基于VBOX设备的ABS性能实验和结果分析[J]. 轻型汽车技术, 2010(5-6): 25-35.) |
| [8] | SHUAI Ping, CHEN Dingchang, JIANG Yong. Error of GPS Broadcast Ephemeris and Their Effects on Navigation and Positioning Accuracy [J]. Journal of Data Acquisition and Processing, 2004, 19(1): 107-110. (帅平,陈定昌,江涌. GPS广播星历误差及其对导航定位精度的影响[J]. 数据采集与处理, 2004, 19(1): 107-110.) |
| [9] | IGS. IGS Product Table: International GNSS Service[EB/OL]. [2012-01-22]. |
| [10] | HE Haibo, YANG Yuanxi, SUN Zhongmiao. A Comparison of Several Approaches for Velocity Determination with GPS[J]. Acta Geodaetica et Cartographica Sinica, 2002, 31(3): 217-221. (何海波,杨元喜,孙中苗. 几种GPS测速方法的对比分析[J]. 测绘学报, 2002, 31(3): 217-221.) |
| [11] | SUN Zhongmiao, SHI Pan, XIA Zheren, et al. Determination of the Vertical Acceleration for the Airborne Gravimetry Using GPS and Digital Filtering[J]. Acta Geodaetica et Cartographica Sinica, 2004, 33(2): 110-115. (孙中苗, 石磐, 夏哲仁, 等. 利用GPS和数值滤波技术确定航空重力测量中的垂直加速度[J]. 测绘学报, 2004, 33(2): 110-115.) |
| [12] | JEKELI C. On the Computation of Vehicle Accelerations Using GPS Phase Acceleration[C]// Proceedings of the International Symposium on Kinematic Systems in Geodesy, Geomatics, and Navigation (KIS94). Bunff: [s.n.], 1994: 473-481 |
| [13] | BRUTON A M, KERN M, SCHWARZ K P. et al. On the Accuracy of Kinematic Carrier Phase DGPS for Airborne Mapping[J]. Geomatica, 2001, 55(4): 491-507. |
| [14] | ARINC. Interface Control Document, Navstar GPS Space Segment/Navigation User Interface, IS-GPS-200[EB/OL]. [2006-05-15]. |
| [15] | JASON Z, ZHANG K F, RON G, et al. GPS Satellite Velocity and Acceleration Determination Using the Broadcast Ephemeris [J]. The Journal of Navigation, 2006, 59(2): 293-305. |
| [16] | CUI Xianqing, JIAO Wenhai, JIA Xiaonlin, et al. The Fitting Algorithm of GPS Broadcast Ephemeris Parameters[J]. Journal of Institute of Surveying and Mapping, 2004, 21(4): 244-246. (崔先强,焦文海,贾小林, 等. GPS广播星历参数拟合算法[J]. 测绘学院学报, 2004, 21(4): 244-246.) |
| [17] | GAO Yudong, XI Xiaoning, WANG Wei. An Improved Fitting Algorithm Design for Broadcast Ephemeris for GEO Satellite[J]. Journal of National University of Defense Technology, 2007, 29(5): 18-22. (高玉东,郗晓宁,王威. GEO导航星广播星历拟合改进算法设计[J]. 国防科技大学学报, 2007, 29(5): 18-22.) |
| [18] | CHEN Liucheng, TANG Bo. Impact of Coordinate Transformation on Fitting Accuracy of Kepler Broadcast Ephemeris[J]. Journal of Spacecraft TT&C Technology, 2006, 25(4): 19-25. (陈刘成,唐波. 参考系选择对Kepler广播星历参数拟合精度的影响[J]. 飞行器测控学报, 2006, 25(4): 19-25.) |
| [19] | Russian Institute of Space Device Engineering. Interface Control Document of Global Navigation Satellite System[M]. 5th ed. Moscow: [s.n.], 2008: 25-43. |
| [20] | GE Kui, WANG Jiexian. Calculation of GLONASS Satellite Station and the Realization of Program[J]. Geomatics and Spatial Information Technology, 2009, 32(2): 137-140. (葛奎,王解先. GLONASS卫星位置计算与程序实现[J]. 测绘与空间地理信息, 2009, 32(2): 137-140.) |
| [21] | KENNEDY S L. Precise Acceleration Determination from Carrier-Phase Measurements [J]. Journal of the Institute of Navigation, 2003, 50(1): 9-19. |
| [22] | ROLF D, URS H, PIERRE F, et al. Bernese GPS Software Version 5.0[R]. Berne: University of Berne, 2007. |