为了评估各系统近期广播星历的精度，考虑到Galileo自2019-07-12全星座暂停向用户发播信号^[4]，直至2019-07-18恢复正常^[5]，所选时段需要满足各系统均正常播发信号的要求，因此，评估时段选定为2019-09-01-2019-11-01，即年积日（day of year，DOY）第244~305天，共计62 d。卫星系统的星座状态包括卫星发射时间、开始服役的时间、已服役时长、轨道类型、搭载卫星钟的种类等，GPS、GLONASS、Galileo、BDS的星座状态均可以在GPS World的THE ALMANAC栏目中查得^[3]，但是该网址不包含BDS搭载的卫星钟情况，BDS详细星座状态可以在中国卫星导航办的官网上查得^[6]。需要说明的是，本文均采用伪随机码（pseudo-random noise，PRN）编号识别卫星，因为所选评估时段内各系统卫星对应的飞行器编号不变，不涉及到卫星型号的更替问题，PRN编号前的大写字母G、R、E、C分别代表GPS、GLONASS、Galileo和BDS。

本文采用多模GNSS实验（multi-GNSS experiment，MGEX）发布的多系统混合广播星历（ftp://cddis.gsfc.nasa.gov/pub/gps/data/campaign/mgex/daily/rinex3/2019/brdm/）与武汉大学（WUM）、中国科学院上海天文台（SHA）分析中心发布的15 min事后精密星历（ftp://igs.ign.fr/pub/igs/products/mgex/）进行星历比较，选用WUM的精密星历产品是因为只有WUM有较完整的BDS-3数据，考虑到2019年DOY 299的WUM精密星历文件缺失，故使用SHA的精密星历替代，该天SHA精密星历对于BDS仅提供BDS-2的数据。

对于所选用的精密星历产品，WUM使用IGb08参考框架，SHA使用IGS08参考框架，两者均每隔15 min（900 s）提供一组各PRN编号的卫星质心在各自参考框架下以km为单位的三维坐标和以μs为单位的钟差。WUM精密星历轨道与钟差的精度与稳定性和MGEX各分析中心持平，对于BDS的轨道精度，在地球同步轨道（geostationary orbit，GEO）、倾斜地球同步轨道（inclined geosynchronous orbit，IGSO）、地球中轨（medium earth orbit，MEO）的精度分别达到分米级、分米级、厘米级^[7]，钟差精度约为0.1 ns^[8]；而SHA精密星历在轨道与钟差精度上与国际GNSS服务组织（international GNSS service，IGS）各分析中心相当^[9]。在评估广播星历精度时，二者提供的坐标与钟差均可作为真值使用。

为评估广播星历精度，需要考虑轨道误差、钟差以及二者的综合影响。在星固坐标系中，利用广播星历与精密星历得到的卫星位置在径向、法向、切向上的坐标差可得到轨道精度；利用广播星历计算得到的钟差与精密星历钟差之差$\mathrm{\delta }{d}_C$可以衡量钟差精度；利用SISRE来衡量包含轨道误差和钟差影响的总体测距精度，计算式为:

式中，R、A、C分别为径向、法向、切向的坐标差；${\omega }_R$、${\omega }_{A,C}^{\mathrm{2}}$是权系数，不同卫星系统和轨道类型的卫星取值不同，详见表 1。

若仅衡量整体轨道精度，可以不考虑钟差项，得到仅含卫星轨道误差的测距差SISRE_O，计算式为:

对上述指标求取均值和均方根（root mean square，RMS），可得到各项指标的统计特性。

在所选评估时段内，GPS、GLONASS、Galileo、BDS的广播星历分别使用WGS-84^[11]、PZ-90^[12]、GTRF^[13]、BDCS^[14]参考框架，与精密星历参考框架间互差最大为厘米级^[15-18]，因此在评估广播星历精度时忽略参考框架不同带来的影响。

通过精密星历头文件可知，所用精密星历均采用GPST（GPS time），因此需要将各系统时转为GPST: （1）对于GLONASS，RENIX 3.03给出的是UTC时间系统下的周内秒，仅需加上2019年对应的闰秒数（18 s）即可转化为GPST，需要注意的是，由于并非所有MGEX混合广播星历头文件中均含有LEAP SECONDS信息，所以不能采取读头文件的方式读取闰秒数；（2）对于Galileo，由于GST（Galileo system time）中Galileo周和GPS周相同，两个时间系统仅差几十ns^[19]，故视为一致；（3）对于BDS，由于BDS周比GPS周晚1 356周（BDT（BeiDou time）开始于2006-01-01 00:00:00，GPST开始于1999-08-22 00:00:00），因此BDT比GPST晚14 s。利用导航星历计算卫星位置和钟差的具体算法参见各系统最新版本的接口控制文件（interface control document，ICD）^[11-13，20]。

广播星历计算得到的卫星位置是卫星天线相位中心的位置，而精密星历提供的是卫星质心的位置，在进行轨道比较时，应把精密星历提供的卫星位置改化至天线相位中心。需要注意的是，此时对精密星历卫星位置加入的天线相位中心改正值不能采用IGS发布的ANTEX文件提供的改正值，虽然IGS发布的天线模型是经平差处理后与真实天线相位中心位置最接近的结果，但是该模型是与IGS发布的轨道及钟差产品配套使用的，和广播星历生成时使用的各GNSS卫星生产商提供的天线相位中心偏移量不能混淆^[21]，广播星历在生成时是将生产商提供的偏移量作为改正值将卫星质心坐标转化为天线相位中心坐标提供给用户使用，因此，在利用精密星历评估广播星历轨道精度时，应使用与广播星历生成相匹配的天线相位中心改正值。

目前，四大GNSS并没有官方公布各自广播星历生成时使用的具体天线相位中心改正值。因此，对于GPS和GLONASS，采取Montenbruck等于2015年^[10]和2017年^[22]通过对比精密星历和广播星历的差值计算得到的天线相位中心偏差；对于Galileo，采取欧洲卫星导航服务中心（European GNSS service centre，ESC）发布的卫星天线相位中心与某一确定物理参考点在星固坐标系中三方向的差值^[23]，同一型号的卫星取均值作为该型号卫星天线相位中心偏移量；对于BDS，早期一些星历质量评估文献由于BDS官方未公布改正参数，加之文献[10]发现BDS二代卫星质心与相心较为接近，于是忽略此项改正^[10，24-25]。随着星座构建日益完善，部分文献采用的是MGEX推荐值^[26]或早期硬件制造商公布的一组改正值^[27]，经验算，以上方法对于BDS广播星历轨道精度的评估仍存在系统性误差，注意到BDS官方于2019-12-30发布的北斗卫星参数文件中包含制造商公布的每一颗BDS卫星的天线相位中心相对于质心的改正值^[28]，采取对各型号卫星三方向改正值求均值的方法，计算得到本文后续采用的天线相位中心改正值。综上，本文对各系统不同型号卫星星固坐标系下x、y、z 3个方向进行天线相位中心改正值的取值，结果见表 2。

广播星历中，不同GNSS系统的钟差参数基准不同，通常定义为某个指定观测量，该观测量可以是单频观测量，如BDS的B3，也可以是双频消电离层组合观测量^[29]，如GPS的L1/L2、GLONASS的G1/G2、Galileo的E1/E5a（对应INAV）或E1/E5b（对应FNAV）^[19]。以双频消电离层组合观测量为例，假设$P_i$（i=1，2）代表某GNSS系统两个频率上的码观测量，其非差观测方程可以简写为:

式中，$\rho$代表接收机至卫星天线相位中心的距离；$c$为光速；$\mathrm{d}{t}_{}^S$为卫星钟差；$K_i^S$是卫星的硬件延迟偏差；${\epsilon }_i$包含接收机钟差和硬件延迟偏差、相对论效应、对流层延迟、电离层延迟和观测噪声。需要注意的是，对于GLONASS，与GPS、Galileo和BDS不同，在广播星历中已经加入相对论改正，用户在使用广播星历进行位置解算时不必考虑该项^[12]，但在星历比较时，为了与精密星历对应，需要将钟差的相对论改正项$\mathrm{\Delta }{\rho }_{\mathrm{r}\mathrm{e}\mathrm{l}}^{}$从GLONASS广播星历计算的钟差中减去，$\mathrm{\Delta }{\rho }_{\mathrm{r}\mathrm{e}\mathrm{l}}^{}$计算式为:

式中，$\mathit{r}$、$\mathit{v}$分别是利用广播星历计算出的卫星位置向量和速度向量。

设$f_i$为频率，令${\tilde{\alpha }}_{\mathrm{1}}=f_{\mathrm{2}}^{\mathrm{2}}/(f_{\mathrm{1}}^{\mathrm{2}}-f_{\mathrm{2}}^{\mathrm{2}})$，${\tilde{\alpha }}_{\mathrm{2}}=f_{\mathrm{1}}^{\mathrm{2}}/(f_{\mathrm{1}}^{\mathrm{2}}-f_{\mathrm{2}}^{\mathrm{2}})$，由$P_i$进行消电离层组合得到的组合观测量$P_C$可以表示为:

式中，${\epsilon }_C$为各模型误差的消电离层组合量；$\mathrm{\delta }{t}^S$即为广播星历所给钟差参数计算得到的卫星钟差。观察式（3）和式（5），得到关系式:

由式（6）可知，$c\mathrm{\delta }{t}^S$与$c\mathrm{d}{t}_{}^S$的差值在星历比较中是一个系统性偏差，对于某PRN编号的卫星来说硬件延迟偏差在一天之内可以看作一个常数。考虑到在评估广播星历精度时，SISRE指的是用户接收机至卫星天线相位中心的测距误差，不应含有系统性偏差，因此，广播星历精度评估采用的应该是$\mathrm{d}{t}_{}^S$所表示的卫星钟差。

由于GPS和Galileo广播星历中给出的钟差基于双频消电离层组合观测量，无需考虑群延迟（time group delay，TGD）改正，GLONASS则是在广播星历中直接给出包含TGD改正的基于G1的钟差，在与精密星历比较时也无需考虑该项改正^[21]。BDS广播星历钟差参数基准为B3，由于WUM分析中心发布的精密星历最终产品自2019年DOY 001开始从使用B1I和B2I消电离层组合改为使用B1I和B3I消电离层组合^[30]，在与精密星历钟差比较时，还应加入TGD改正，使其与精密星历消电离层组合基准一致。TGD改正公式为:

式中，$f_{\mathrm{B}\mathrm{1}}^{\mathrm{2}}$、$f_{\mathrm{B}\mathrm{3}}^{\mathrm{2}}$为B1、B3频点的频率值；$\mathrm{T}\mathrm{G}{\mathrm{D}}_{\mathrm{1}}$为BDS广播星历中B1频点的群延迟改正数。

由于DOY 299采用SHA精密星历来评估当天广播星历的精度，虽然精密星历的互差在广播星历精度评估时可忽略不计，但目前尚无SHA在产生精密星历时使用的组合观测量有关的参考文献，因此，式（7）的适用性无法确定。根据TGD的数值由载波的频率和该频率上卫星的硬件延迟偏差决定的性质，可以把TGD对单颗星一日内的影响视为系统性偏差处理。

消除广播星历钟差系统性偏差的方法有很多: （1）选取一颗卫星为基准星，其余卫星与其作二次差^[31-33]，但是该方法对基准星要求较高，基准星选取不当会将基准星误差引入其他卫星钟差；（2）求取一个GNSS系统内所有可用卫星钟差的均值，其余卫星与其作二次差^{[10，24-25，34]}，但是该方法仅消除了星座整体时标偏差，每颗卫星的钟差仍然包含一个均值偏差^[35]；（3）利用加权最小二乘进行迭代，以残差最小为目标函数，计算所有可用卫星的钟差一次差来求解公共偏差^[36]，但是该方法计算量大，过程较为繁琐。为了最大限度地消除系统性偏差，本文采用一种对单颗星一日内广播星历与精密星历钟差之差求均值，在逐历元与精密星历比较时，不同历元的钟差再与该均值作差的方法，消除该星钟差的系统性偏差。经验证，采用该方法作二次差，在计算广播星历钟差时是否加入TGD改正对最终钟差精度没有影响，也证明了TGD和其他系统误差通过上述二次差方法均可得到有效消除。

以精密星历的历元间隔15 min（900 s）为比较间隔，每隔15 min选择健康卫星的数据计算一次广播星历各系统各PRN编号卫星的位置与钟差，GPS、GLONASS、Galileo、BDS与精密星历的比较结果分别如图 1~图 4所示。

图  1  GPS广播星历与精密星历对比结果

Figure 1.  Comparison of GPS Broadcast Ephemeris and Precise Ephemeris

图  2  GLONASS广播星历与精密星历对比结果

Figure 2.  Comparison of GLONASS Broadcast Ephemeris and Precise Ephemeris

图  3  Galileo广播星历与精密星历对比结果

Figure 3.  Comparison of Galileo Broadcast Ephemeris and Precise Ephemeris

图  4  BDS广播星历与精密星历对比结果

Figure 4.  Comparison of BDS Broadcast Ephemeris and Precise Ephemeris

需要指出的是，在画图过程中保留了所有数据，不剔除粗差，以便观察各误差变化规律，各图块纵轴范围内显示的数据均不包含误差值过大的个别数据点。各系统不同类型卫星的轨道精度、钟差精度、SISRE的RMS等统计结果详见表 3。在统计精度时需要进行粗差的剔除，本文设置的剔除条件为: 若某颗星在某历元SISRE的RMS大于10 m，则判定为粗差，并剔除该星在该历元所有指标数据。对于BDS，所使用广播星历BDS在GEO、IGSO的卫星钟龄期（age of data clock，AODC）、星历龄期（age of data ephemeris，AODE）数值均小于25，MEO的AODC、AODE 24 h以内的各占98.32%、98.51%。需要指出的是，BDS广播星历精度与数据龄期并非简单的线性相关的关系，其中还存在周期性变化^[37]，这里提供的龄期状态仅为后续比较结果提供参考。

结合图 1~图 4与表 3，具体分析各系统各指标精度如下:

1）对于GPS，由图 1可以看出，径向、法向、切向和钟差的误差值均围绕0刻度上下波动，说明计算结果不含系统性偏差，3个方向中径向精度最高，且数据质量较为稳定，没有过多的散点。结合表 3可知，整体轨道精度的RMS优于0.3 m，SISRE的RMS优于0.5 m。其中，Block IIA由于型号老旧，钟差的RMS大于同样将Rb钟作为工作钟的Block IIR/IIR-M和Block IIF型号的卫星，精度较差；对于Block IIF卫星，使用Cs钟的两颗卫星钟差的RMS大于使用Rb钟的其他同型号卫星，由于轨道精度相近，且根据式（1），其SISRE的RMS也大于其他同型号卫星。对于使用Rb钟的GPS卫星，从表 3中可以看出，Block IIA至Block IIF，钟差与SISRE的精度均逐渐提高。

2）对于GLONASS，由图 2可以看出，轨道误差中3个方向以及钟差的误差值在四大导航系统MEO轨道卫星中是最大的，但是通过图 2中数据点分布可以看出，其数据不存在明显的分散以及系统性偏差，由文献[3]可知，GLONASS全星座均使用Cs钟，结合表 3，该系统整体钟差的RMS为1.5 m左右，轨道精度的RMS优于0.8 m，SISRE的RMS在2 m以内。由于统计时段内有数据参与统计的GLONASS-K1卫星仅有一颗，且由于该统计仅反映近期精度情况，因此K1型号卫星在该统计时段内整体精度稍差。

3）对于Galileo，由图 3可以看出，除数据没有明显的系统性偏差外，3个方向的误差较小且精度稳定，但钟差存在一条稀疏的离散点条带，且个别日期的数据存在整体的线性变化，由式（1）知，SISRE数值也随之变化，结合表 3，该系统整体轨道精度的RMS优于0.2 m，钟差和SISRE的RMS均优于0.4 m。使用Rb钟的E11卫星广播星历预报钟差的RMS稍大于使用被动型氢钟的其他卫星，对于全部使用氢钟的FOC卫星，整体钟差精度的RMS要稍好于测试用的IOV卫星。

4）对于BDS，除各项精度指标不含系统性偏差外，由图 4（a）可知，由于轨道高度原因，BDS-2的GEO卫星在法向和切向上精度明显较差，且数值变化周期长，有较多的数据分散点；由图 4（b）可知，所有指标在DOY 299发生数据缺失，这是因为该日使用的SHA精密星历中没有BDS-3的数据，BDS-3轨道精度明显优于BDS-2，钟差精度较BDS-2也更加接近0，但相比于其他3个系统，BDS总体的钟差精度稳定性较差，图 4中表现为钟差图块中离散的点数较多。结合表 3，BDS-2钟差精度由高到低按轨道高度排列依次为MEO、IGSO、GEO，BDS-2的整体轨道精度优于0.7 m，SISRE的RMS优于1 m。相比BDS-2，BDS-3的轨道精度和钟差精度均有较大提高，SECM制造的BDS-3卫星使用氢钟，该批卫星整体广播星历预报钟差精度要优于CAST制造的使用Rb钟的BDS-3卫星，BDS-3整体轨道精度的RMS在0.2 m以内，SISRE的RMS在0.4 m以内。同时，由于BDS-3可用观测值数量低于BDS-2，因此WUM精密星历BDS-3的轨道精度和钟差精度较低^[30]，对BDS-3广播星历精度评估有一定的影响。

相比于已有的评估工作，本文在评估内容与策略上的改进主要体现在以下几个方面: （1）评估时段较新，对于BDS，采用了最新发布的制造商对不同型号卫星天线相位中心相对于质心偏移量的标定值进行天线相位中心改正。（2）对于钟差的系统性偏差的处理，在广播星历质量评估类文章中，首次利用非差观测值方程，推导出广播星历所提供钟差的具体表达形式，为充分消除系统性偏差提供理论依据，并提出一种不依赖于精密星历所使用的消电离层的频点的系统性偏差消除方法。（3）涵盖四大导航系统，目前，比较时段较近的涵盖四大导航系统广播星历质量评估的工作是2018年Montenbruck等^[21]开展的，随着2019年WUM改变精密星历消电离层使用的频点以及BDS官方发布北斗卫星参数文件，原有的四大系统广播星历质量评估与比较的结果不再具有时效性。

从本文评估结果来看，四大导航系统广播星历整体精度从高到低依次为Galileo、GPS、BDS、GLONASS。其中，Galileo整体轨道精度优于0.2 m，钟差和SISRE的RMS均优于0.4 m；GPS整体轨道精度优于0.3 m，钟差的RMS优于0.4 m，SISRE的RMS优于0.5 m；对于BDS，BDS-2整体轨道精度优于0.7 m，钟差的RMS优于0.6 m，SISRE的RMS优于1 m；BDS-3的轨道精度优于0.2 m，钟差的RMS优于0.4 m，SISRE的RMS优于0.4 m，且BDS-3广播星历的整体轨道精度、钟差精度和SISRE的RMS和四大系统整体结果相比均最优；由于采用了最新的天线相位中心改正值以及新的钟差系统性偏差处理方法，相比于刘凡等^[34]的评估结果，BDS-2与BDS-3 SISRE的RMS精度均提高0.2 m左右；BDS总体轨道精度和钟差精度均在0.4 m左右，SISRE优于0.6 m；GLONASS整体轨道精度在0.5 m左右，钟差的RMS约为1.5 m，SISRE的RMS约为1.6 m。此外，广播星历钟差精度不完全取决于钟的稳定度，更取决于预报时长的差异: 以BDS为例，根据毛亚等^[38]的评估结果，氢原子钟千秒稳定性、万秒稳定性和日稳定性分别高于铷原子钟40.16%、47.10%、65.60%，且氢原子钟的长期稳定性要明显优于铷原子钟，但是运控系统监测的各卫星导航星历钟差精度差别并不明显。造成各系统广播星历精度评估结果差异的原因也是多方面的，如预报时长、原子钟类型、信号调制方式、信号频率、地面控制中心的布址分散程度、空间星座的轨道设计、不同生产厂家生产卫星的工艺，以及生成广播星历与精密星历采用的算法等。随着各系统卫星的更新换代，各GNSS广播星历的精度还会有所提高，尤其是BDS，随着BDS-3的GEO和IGSO卫星逐步替代原有的BDS-2相应轨道高度的卫星，BDS整体广播星历精度还会有较大的提升。