0.   引　言

地球定向参数(earth orientation parameter，EOP)是国际天球参考系统(International Celestial Reference System，ICRS)与国际地球参考系统(International Terrestrial Reference System，ITRS)转换的关键参数，其包含丰富的地球动力学信息^[1-3]，其精度直接影响卫星^[4]或深空探测器的精密定轨精确度^[5]，在大地测量领域具有极其重要的作用. EOP由表示地球自转轴指向信息的一系列可观测的参数表示，包含岁差和章动参数，极移(polar motion，PM)参数及日长变化(length of day，LOD)^[6]. 其中，岁差和章动参数描述地球自转轴相对于天球的运动，它们的变化可以通过现有的IAU2000和IAU2006模型进行精确预测，预测精度优于1 mAs^[7]. 在精度要求不严格的实际应用中，可以忽略其预测精度不足带来的影响，故后续不再讨论. 此外，极移参数描述地球自转轴相对于地球本体之间的运动，LOD参数表征地球自转角速率的变化，直接观测量为考虑极移改正的世界时(universal time，UT1)与协调世界时(coordinated universal time，UTC)之间的差值，它们统称为地球自转参数(earth rotation parameters，ERP).

为推进导航信息文件现代化，以适应所有GNSS星座的导航信息与系统数据信息，接收机独立交换(Receiver Independent Exchange，RINEX)4.00协议在2021年12月举行的第59届理事会会议上正式获批. 该格式的导航星历文件首次包含ERP数据记录：极移$ {x_{{p}}} $、$ {y_{{p}}} $，$ {\Delta _{{\text{UT1}}}} $及其变化速率^[8-9]，以便开展低地球轨道(low earth orbit，LEO)卫星和高地球轨道(high earth orbit，HEO)卫星的星载实时定轨研究^[10]. 目前有美国的GPS、中国的北斗卫星导航系统(BeiDou Navigation Satellite System，BDS)2个全球卫星导航系统与日本的准天顶卫星系统(Quasi-Zenith Satellite System，QZSS) 、印度区域导航卫星系统(Indian Regional Navigation Satellite System，IRNSS)2个区域卫星导航系统在导航电文中播发ERP.

由于ERP对卫星轨道精度具有重要意义^[11]，其预测精度将直接影响广播星历精度从而造成定位结果的偏差^[12]. 不同卫星导航系统的ERP预报精度存在差异，从而导致广播星历实现的地球框架不一致^[13]，进而影响用户端的定位、导航与授时服务^[14-15]. 在Steigenberger等^[8]和Chen等^[16]的研究中，发现广播星历和精密轨道之间的旋转参数与ERP误差之间存在非常一致的变化趋势，表明卫星广播星历中播发的ERP与地面生成广播轨道时使用的ERP很可能是相同的. 因此，广播ERP的精度可以提供星历误差相关的信息. 本文针对RINEX4.00格式的广播星历数据，着力于GPS和BDS-3播发的ERP参数分析评估，并设计了基于广播ERP的星历误差实时补偿方案.

1.   数据与方法

1.1   数据集

本研究共收集德国宇航中心(German Aerospace Center，DLR)发布的自2022年1月到2023年12月共24个月的RINEX4.00广播星历文件(后文简称BRD4)^[9]. 另外，为评估广播星历ERP参数，将 ERP与IERS 发布的最终解 C04 序列进行比较，通过求取偏差进行精度评估. C04产品由多种空间大地测量技术联合观测获得^[17]，具有较高的精度，已在相关研究中取得广泛应用.

为估计广播星历坐标可能存在的旋转误差，本文采用国际GNSS服务(International GNSS Service，IGS)分析中心德国地学研究中心(German Research Centre for Geosciences，GFZ)发布的精密轨道产品作为参考，通过7参数赫尔默特变换方法估计旋转参数. 由于本文的研究针对GPS与BDS-3，最好选择在研究期内生成策略一致的轨道产品，GFZ发布的精密星历产品包含多星座系统且数据质量稳定^[18-19]，被各地学者广泛采用，故将其作为本研究的参考星历.

1.2   ERP参数计算

BRD4文件中ERP参数记录包含参考历元${t_{{r}}}$的$ {x_{{p}}} $、$ {y_{{p}}} $、$ {\Delta _{{\text{UT1}}}} $及其一、二阶的变化速率$ {\dot x_{{p}}} $、$ {\dot y_{{p}}} $、$ {\dot \Delta _{{\text{UT1}}}} $、$ {\ddot x_{{p}}} $、$ {\ddot y_{{p}}} $、$ {\ddot \Delta _{{\text{UT1}}}} $. 另外，BRD4中还引入1个额外参数${t_{{m}}}$，用于表示接收到该条ERP数据记录的第1个传输历元. 由于实验阶段内BRD4 播发的ERP参数的二阶变化速率均为0，因此在计算中只考虑一阶变化速率，某个特定历元$t$的ERP参数可以通过以下方式进行预测：

式中，${t_{{{r,i}}}} < t < {t_{{{r,i + 1}}}}$，下标$i$表示BRD4文件中第$i$条数据记录.

1.3   轨道旋转参数估计

由于BRD4中的ERP参数是地面站预测得到的，并且合理推测在广播星历的生成中也采用同样的预测ERP参数^[20]. 而该预测ERP与各分析中心采用的PM解之间存在偏差，从而造成广播星历与精密星历存在显著的不一致性. 为评估该不一致性与ERP预测误差的关系，本研究采用7参数赫尔默特变换方法估计广播星历与精密星历之间的旋转参数：

式中：$ {X_{{\text{pre}}}} $、$ {Y_{{\text{pre}}}} $、$ {Z_{{\text{pre}}}} $与$ {X_{{\text{brd}}}} $、$ {Y_{{\text{brd}}}} $、$ {Z_{{\text{brd}}}} $分别为精密星历、广播星历卫星在地心地固坐标系中的位置；${R_{{x}}}$、${R_{{y}}}$、${R_{{z}}}$分别为待估的3轴旋转参数；$ {\Delta _{\,{{X}}}} $、$ {\Delta _{\,{{Y}}}} $、$ {\Delta _{\,{{Z}}}} $分别为3轴平移分量，$m$为比例因子. 为了构造最小二乘法方程，令

$ {a_1} $~$ {a_4} $为中间变量，代入式(2)中，可以得到

式中变量含义与式(2)一致. 在采用精密星历和广播星历进行参数估计时，要注意精密星历的参考点为卫星质心，而广播星历的参考点为天线相位中心（antenna phase center，APC），需要将精密星历的参考点转换到APC才能够进行对比. 但是不能采用IGS发布的ATX产品进行APC改正，因为该产品中的改正值是事后经过平差、考虑方位角及高度角对改正值变化影响的结果^[21]，需要与其他轨道产品搭配使用，因此本文采用各个卫星制造商提供的在星历生产时应用的改正模型进行相位中心改正. 然而四大GNSS官方并没有完全公布上述提到的相位中心改正值，根据张熙等^[21]的研究，采用其收集计算的改正值作为本次研究中的改正取值，具体改正数如表1所示.

表  1  GPS、BDS卫星APC改正值

卫星系统	卫星型号	PRN编号	Δ$x $/m	Δ$y $/m	Δ$z $/m
GPS	Block ⅡA	G18	+0.279	0.000	+0.920
	Block ⅡR-A	G11, G13, G14, G16, G20, G21, G28	0.000	0.000	+1.610
	Block ⅡR-B	G02, G19, G22, G23	0.000	0.000	−0.040
	Block ⅡR-M	G05, G07, G12, G15, G17, G29, G31
	Block ⅡF	G01, G03, G06, G08~G10, G24~G27, G30, G32	+0.394	0.000	+1.160
BDS	2G	C01~C05	−0.615	−0.570	+1.100
	2I	C06~C10, C13, C16	−0.419	−0.573	+1.100
	2M	C11, C12, C14	−0.428	−0.179	+1.100
	3M	C19~C30, C32~C37	+0.613	+0.128	+0.953

下载: 导出CSV 

| 显示表格

2.   广播星历ERP参数评估

2.1   数据编辑方法

本节首先提取2022—2023年间播发的BRD4中包含的ERP参数，结果发现有效ERP记录共3845条，其中，GPS包含2640条数据记录，每天约有不同卫星播发的1~4条，并且，通过观察传输历元${t_{{m}}}$发现，GPS大致在午夜22:30—01:00之间播发ERP. BDS-3包含1205条数据记录，每天1~2条，其传输时间大致在凌晨00:00—01:00，偶尔出现14:00—16:00传输的记录. 通过对数据的筛选发现存在大量重复数据，其具体表现为不同传输历元接收到的数据记录完全相同. 通过筛查发现，较早播发的重复ERP记录是正式播发数据的一个预播发，含有较大粗差的数据记录均为重复记录中传输历元靠前的记录. 因此，本文采用一种简明有效的数据编辑策略：对重复数据记录保留与参考历元间隔最近的一组，即传输历元靠后的数据记录.

除完全一致的重复播发数据记录，在BDS-3数据中还包含参考历元与PM参数相同而$ {\Delta _{{\text{UT1}}}} $不同的情况，对照C04产品，将偏差较大的数据记录视为粗差进行剔除. 另外，GPS的ERP参考历元均在16:00—20:00，而BDS-3大部分在00:00时刻附近，但少量数据参考历元在14:00之后，通过检验发现参考历元异常的数据均含有较大粗差，本文将其作为错误播发数据予以剔除. 此外还发现除参考历元外其他参数均为0的异常记录，直接剔除.

通过上述数据编辑，GPS与BDS-3广播ERP基本上都保持每天1~2条的播发频率，最终参与实验对比分析的有效ERP数据GPS有732条，剔除率高达72.27%(1 908条/2640条)，BDS-3有730条，剔除率达39.42%(475条/1205条). 可以发现，广播ERP参数的播发中，大部分数据为重复播发或错误播发，为保证后续评估的可靠性必须进行细致的数据筛查与编辑.

2.2   数据精度分析

本文通过与精密后处理ERP产品C04作差来评估GPS与BDS-3广播ERP的精度. 根据GPS、BDS接口控制文件以及IERS公约，广播ERP参数包含高频PM变化(周期小于1 d)，而在C04中该变化已被修正. 为了更为严谨的研究结果，在前文数据编辑的基础上，需要对广播ERP参数进行该高频项的修正. 图1展示了从2022—2023年共730 d的GPS、BDS-3广播星历ERP相对于C04的差值序列. 图中，左右图像分别对应GPS，BDS-3，横坐标表示自从2022年1月1日开始计数的天数. 可以看到，GPS与BDS-3广播ERP差值整体在0值附近，伴有小周期的上下波动，这是因为ERP的线性预测模型误差并不能完全表征ERP的变化，以及BRD4生产时的ERP低频观测并不能表征其高频极移变化. 对结果进行量化统计，如表2所示.

图  1  广播ERP参数相对于C04的差值序列

下载: 全尺寸图片 幻灯片

表  2  广播ERP相对于C04差值的RMS统计

参数	GPS	BDS-3
${x_{{p}}}$/mAs	2.309	2.769
${y_{{p}}}$/mAs	1.433	1.748
$ {\Delta _{{\text{UT}}\,{\text{1}}}} $/ms	0.303	0.372

下载: 导出CSV 

| 显示表格

其中，GPS广播${x_{{p}}}$差值波动在$ \pm7\; \mathrm{mAs} $以内，${y_{{p}}}$差值波动在$ \pm 5\;{\mathrm{mAs}}$以内，$ {\Delta _{{\text{UT1}}}} $差值波动在$ \pm 2\;{\mathrm{ms}}$，3个参数相对于C04差值的RMS分别为：2.309 mAs、1.433 mAs、0.303 ms. 而BDS-3广播${x_{{p}}}$差值波动在$ \pm 11\;{\mathrm{mAs}}$以内，${y_{{p}}}$差值波动在$ \pm 8\;{\mathrm{mAs}}$以内，$ {\Delta _{{\text{UT1}}}} $差值波动在$ \pm 3\;{\mathrm{ms}}$以内，3个参数相对于C04差值的RMS分别为：2.769 mAs、1.748 mAs、0.372 ms. 通过图1和表2可以看到，BDS-3广播ERP误差相对于GPS的波动幅度较大，稳定性较差， BDS-3 ERP参数相对于C04差值的RMS相较于GPS增大19.9%、22.0%与22.8%.

本文将实验期内的数据划分为不同的BDS-3周内日进行精度指标统计. 图2中，(a)、(b)2个子图分别表示GPS与BDS-3的广播ERP相对于C04的差值平均值与标准差(standard deviation，STD)统计，不同颜色的字体代表对应颜色ERP参数差值的RMS值，柱状图表示其平均值. 可以看到，GPS广播ERP差值统计值在不同周内日的结果比较一致，而BDS-3广播ERP在序号为5的周内日平均值与标准差最小，整体依周内日呈现出逐渐递增的趋势. 可以推断由于GPS与BDS-3地面站使用的ERP预测模型、更新频率不同，BDS-3的ERP预测更新频率较低，误差随时间累积，造成不同BDS-3周内日误差的递增趋势，且预报模型预报误差较大，从而导致BDS-3的广播ERP精度略低.

图  2  不同系统北斗周内日ERP差值统计

下载: 全尺寸图片 幻灯片

3.   BDS-3星历误差补偿

3.1   可行性分析

依前文，广播星历与精密星历之间的旋转参数与广播ERP参数误差之间存在一致的变化趋势，而BDS-3相较于GPS的ERP参数精度较低，所以推测BDS-3的星历误差可以通过广播ERP进行补偿. 为验证其可行性，本节进一步研究星历误差与广播ERP参数相对于C04的差值之间的相关性.

为统一轨道旋转参数与广播ERP参数中的$ {\Delta _{{\text{UT1}}}} $，本研究通过线性变换，导出地球旋转角(earth rotation angle，ERA)：

式中，$ {T_{{\text{JD,UT1}}}} $表示待求UT1时刻的儒略日，$ {\theta _{{\text{ERA,UT1}}}} $表示该时刻的ERA结果.

图3展示了广播ERP相比于C04的差值${d_{{{xp}}}}$、${d_{{{yp}}}}$、${d_{{\text{ERA}}}}$与轨道旋转参数${R_{{x}}}$、${R_{{y}}}$、$ {R_z} $的变化趋势，图中用对应曲线的颜色文字标注出了每个变量的平均值与STD. 可以看到，对于轨道旋转参数的序列，无论是GPS还是BDS-3，其与广播ERP PM参数差值序列具有相似的变化幅度和波动趋势，整体呈现出高度相关性. 具体来说，ERP PM参数${x_{{p}}}$、${y_{{p}}}$差值分别与轨道旋转参数${R_{{y}}}$、${R_{{x}}}$具有一致性. 但这一结论对$ {\Delta _{{\text{UT1}}}} $差值序列并不适用，其与${R_{{z}}}$旋转参数并没有呈现出与前者类似的相关性. 进一步地，绘制BDS-3广播ERP差值与轨道旋转参数的相关性图.

图  3  广播ERP参数相对于C04的差值与轨道旋转参数序列

下载: 全尺寸图片 幻灯片

图4中黑色虚线为拟合直线，拟合的回归方程与相关系数也在图中进行标注. 可以看到， BDS-3的广播${x_{{p}}}$差值与${R_{{y}}}$参数表现出较强的相关性，相关系数达到0.76，${y_{{p}}}$差值与${R_{{x}}}$参数相关系数达到0.66，而$ {\Delta _{{\text{UT1}}}} $差值与${R_{{z}}}$参数的相关系数仅为0.16，相关性较弱，与前文结论一致.

图  4  BDS-3广播ERP参数相对于C04的差值与轨道旋转参数相关性

下载: 全尺寸图片 幻灯片

3.2   实验验证

本节尝试采用广播ERP参数补偿星历误差，一种简单的方法见下式

式中：${\varPhi _{{X}}}$、${\varPhi _{{Y}}}$、${\varPhi _{{Z}}}$为从广播ERP中导出的旋转矫正参数；${X_{{{{t}}_{{i}}}}}$、${Y_{{{{t}}_{{i}}}}}$、${Z_{{{{t}}_{{i}}}}}$为广播星历计算的卫星${t_{{i}}}$时刻在地固系的位置；$\Delta {X_{{{{t}}_{{i}}}}}$、$\Delta {Y_{{{{t}}_{{i}}}}}$、$\Delta {Z_{{{{t}}_{{i}}}}}$为该时刻卫星位置的旋转补偿值. 本文中用于轨道误差补偿的${\varPhi _{{X}}}$、${\varPhi _{{Y}}}$、${\varPhi _{{Z}}}$参数可以用下式计算

式中：$ {x_{p,{\text{BDS-3}}}} $，$ {y_{p,{\text{BDS-3}}}} $参数分别为BDS-3广播ERP PM参数；$ {y_{p,{\text{GPS}}}} $、$ {x_{p,{\text{GPS}}}} $为用于补偿的GPS广播ERP极移参数. 由于$ {\Delta _{{\text{UT1}}}} $差值与${R_{{z}}}$参数相关性弱，所以不利用其进行改正，将$ {\varPhi _{{Z}}} $参数直接设为0. 为验证GPS广播ERP的实时补偿效果，本研究采用了2种改正策略对BDS-3 MEO卫星星历误差进行补偿：一种采用GPS广播ERP参数计算BDS-3广播ERP参数偏差，标记为COR_GPS方法；另一种采用C04事后产品计算BDS-3广播ERP偏差用于星历旋转补偿，标记为COR_C04方法.

图5展示了2022—2023年期间不同策略下的轨道改正结果. 从上到下分别表示BDS-3系统的地心地固坐标系下每颗卫星的轨道误差(与精密星历产品对比)在卫星轨道坐标系下3轴分量星历误差的平均值. 蓝色、黄色、橘红色分别代表未经过补偿的原始星历误差、经过GPS广播ERP、C04精密ERP产品补偿的星历误差. 可以看到，2种补偿策略均有效减小了广播星历误差，采用GPS广播ERP参数进行补偿后，星历误差最高减小3.1 cm，而采用C04产品进行补偿，误差最高减小5.2 cm，对星历误差的补偿效果明显. 但是可以观察到个别卫星补偿后特定方向的误差更大，若采用三维误差作为评价指标，实验中每颗卫星的星历误差均得到有效补偿. 无论使用何种策略，径向的星历误差没有得到补偿，这是因为旋转改正不能作用于径向. 对所有卫星的改正结果进行汇总，统计不同改正策略下各轴改正值的RMS，可以得到表3.

图  5  采用不同补偿方式对BDS-3 MEO卫星进行星历误差改正结果

下载: 全尺寸图片 幻灯片

表  3  不同策略的星历误差修正值RMS统计

方向	ERP_GPS/cm	C04/cm	占比/%
切向	1.56	2.84	54.93
法向	0.39	0.88	44.32
径向	0.00	0.00	-

下载: 导出CSV 

| 显示表格

可以看到，采用GPS广播ERP参数进行补偿后，星历误差修正值RMS为1.56 cm与0.39 cm，而采用C04产品进行补偿，误差修正值的RMS为2.84 cm、0.88 cm，C04产品达到了更好的误差改正效果，进一步验证了采用ERP参数改正星历误差的有效性. 表格最后一列为2种补偿策略的改正值RMS的比值，表示采用实时播发的GPS广播ERP进行星历误差补偿，在切向、法向分量上分别可以达到采用事后精密产品补偿效果的54.93%和44.32%，这一实验结果也验证了采用实时的GPS广播ERP参数进行BDS-3的星历误差补偿的可行性.

4.   结　论

本文设计了一套广播星历ERP重复记录与粗差的剔除方法，对2022—2023年的GPS与BDS-3广播ERP数据记录进行筛选，GPS与BDS-3分别剔除72.27%、39.42%的ERP记录，然后采用C04评估GPS与BDS-3广播ERP参数的精度. 结果表明：BDS-3广播ERP参数相对于C04的差值相较于GPS广播${x_{{p}}}$，${y_{{p}}}$以及$ {\Delta _{{\text{UT1}}}} $差值的RMS增大19.9%、22.0%与22.8%，精度较低. 本文还研究了广播ERP参数与星历误差的相关性，进行星历误差补偿的实验验证，发现采用实时GPS广播ERP进行改正在切向、法向上分别能达到采用C04产品改正效果的54.93%和44.32%. 上述结论有助于GPS、BDS-3广播ERP参数的筛选、处理及其精度的度量. 并且，证实了采用实时的GPS广播ERP补偿BDS-3星历误差的可行性，对于实时精密定轨有一定的参考价值.

致谢：感谢德国宇航中心(DLR)提供的RINEX4.00格式的广播星历、地学研究中心(GFZ)提供的精密星历产品，以及国际地球自转服务组织(IERS)提供的C04产品.