高精度地球定向参数(Earth orientation parameter,EOP)是地球参考框架和天球参考框架之间相互转换的重要参数，也影响时间系统( 协调世界时UTC与世界时UT1) 的相互转换，在飞行器的精密定位和导航方面^{[1, 2]}应用广泛并且包含了丰富的地球动力学信息。目前主要利用甚长基线干涉测量(very long baseline intergerometrh,VLBI)观测数据解算EOP，精度高且稳定性强^[3]，但是测站较少且不连续；而全球定位系统(GPS)数据能够解算地球自转参数，在GPS测站很多且分布较为理想的情况下，其解算精度可以与VLBI解算精度相当。GPS技术能够提供高精度、高分辨率的地球自转参数(Earth rotation parameter,ERP) 序列值，有效地弥补了VLBI技术观测密度不够的缺陷，但是，GPS解算地球自转参数的稳定性有待提高。联合GPS和VLBI技术解算地球自转参数能较好地利用两者的优势，弥补互相的不足，有望获得精度高和稳定性好的ERP和日长变化解算结果。

采用双差模型解算地球自转参数^{[4, 5]}，待估参数为站坐标、卫星坐标、地球自转参数和中性大气延迟。GPS载波相位观测值表示为待估参数的函数模型：

式中，M表示观测量与待估参数的函数模型；t为时间；X_SP为初始时刻的轨道根数和摄动参数(辐射压模型)；X_T为测站坐标；X_N为相位模糊度；X_erp为地球自转参数；X_atm为大气延迟；ε为观测噪声。

将式(1)线性化为：

式中，C₀表示由近似参数计算出的理论观测值。地球自转参数偏导数可表示为：

式中,r和R₁分别表示为卫星和测站在惯性坐标系中的位置矢量;ρ为站星之间的距离;分别包括和3部分,x_p和y_p表示极移在x、y方向的分量，D_R表示UT1-TAI的一阶变化率。

计算星间距离ρ时，卫星及接收机的位置矢量要换算到同一坐标系下，通常为天球参考系统ICRS。但是，接收机的坐标一般都是在国际地球参考系统ITRS下给出的。所以需要将ITRS转换到ICRS下，其转换关系为R(t_i)=PNSWR₀ⁱ。其中,P、N为岁差、章动旋转矩阵，S为地球自转矩阵，W为极移矩阵。

其中,站坐标与卫星轨道参数来自IGS结果。通过以上解算，可求出地球自转参数信息。

通常用日长(length of day，LOD)的观测值相对于标准日长(86 400国际时秒)的变化ΔLOD来表示地球自转速率的变化^[6]。通过对GPS观测量进行解算便可得到相应的UT1-TAI。ΔLOD的计算方法为:

设VLBI的基本观测方程为：

式中，O_t是在t时刻获得的时间延迟和延迟率观测量；X是由与观测量有关的参数组成的向量；C(X,t)为参数X对观测量O_t影响的数学模型，也是观测量的理论值；V_t是观测量的噪声残差向量。

设参数X由先验值X₀和改正值x两部分组成，则将观测方程线性化有：

由式(6)整理得线性化观测方程为：

式中,y_t=O_t-C(X,t)，是t时刻观测值与理论值之差；A_t为偏导系数矩阵，它与参数值相对于观测值的变化有关。

若将观测方程用矩阵形式表示，则有:

并设观测权阵为P，则根据最小二乘理论可求接触参数的修正值为：

利用卡尔曼滤波进行VLBI大地测量的参数解算^[8]时，将用随机模型取代钟和大气延迟的多项式参数模型，对于非随机参数随机扰动矢量为零，这些参数包括台站位置、射电源坐标和地球定向参数等，而随机参数则包括用于表示钟和大气延迟扰动的随机过程分量。

对GPS和VLBI的解算结果进行联合，需要选取GPS和VLBI同一观测时间段，经过对比，本文选取了2005年9月和2006年1、2月的数据。在GPS方面，为了获得较好的精度而采用了24 h观测资料来进行求解。本文选取了全球均匀分布、几何结构良好、相对稳定的17个IGS站(具体测站名见表 1)的观测数据来解算地球自转参数。VLBI方面，本文选择国际上通用的X波段(波长约为3.6 cm)中观测精度较高的XA和XE两种类型观测数据进行分析。

GAMIT解算地球自转参数时，大部分的参数设置为默认值，详见文献[9]。

用标准均方根NRMS的大小来评定解算结果，一般为0.25左右。若大于0.5，表明处理结果有问题，需要重新处理基线;若NRMS小于0.5，则认为解算成功。本文解算的NRMS均在0.18~0.25范围内。最后得到的ERP和ΔLOD是每天UTC12:00:00的参数值，将其与IGS结果进行对比来评定外符合精度。差值见图1。

图 1 GPS解算值与IGS结果比较 Fig. 1 Comparation of GPS-Estimated ERP and ΔLOD to IGS Solutions

图选项

计算得出极移在X与Y方向上的分量与IGS结果比较的均方根分别为0.32 mas与0.28 mas，UT1-UTC和ΔLOD的均方根分别为0.029 ms和0.044 ms/d。

对2006年的解算结果先去除趋势项，然后进行频谱分析，发现极移在X方向与Y方向上的分量、UT1-UTC和ΔLOD存在周期项，分别在14 d和7 d左右。考虑到IGS是采用周解的方式得到结果，故与其做差后的差值序列会有相应的周期项。

本文用OCCAM 6.2处理IVS服务器上相同时段的XA、XE类型的NGS数据，为了与GAMIT保持一致，选择MHB2000章动模型。将解算结果与IVS结果进行比较，差值见图2。用差值的均方根来评定外符和精度，计算得VLBI数据解算结果与IVS结果的差值在极移X方向与Y方向上的均方根分别为0.45 mas与0.58 mas，UT1-UTC和ΔLOD的均方根分别为0.020 ms和0.059 ms/d。

先用三次拉格朗日插值对GPS和VLBI的解算结果进行插值，得到与IERS 08C04相同历元的ERP和ΔLOD序列；再根据误差传播定律，利用解算结果的中误差分别求出VLBI和GPS解算结果插值后序列的中误差；再以此中误差平方的倒数为权做加权平均，公见如下：

式中，l_i 为拉格朗日插值多项式的系数；t表示解算结果对应的历元；ERP_t 表示解算结果序列；σ_ERPt表示ERP_t的中误差；T表示IERS08C04公布值的历元；ERP(T)表示插值后的序列；σ_ERP(T)为插值后的序列的中误差；ERP_G和ERP_V分别表示利用GPS和VLBI结果插值后的序列；ERP_G+V表示加权平均后的序列；ERP_IERS表示IERS 08C04序列；RMS_G+V表示加权平均后序列与IERS 08C04公布值的均方差。计算结果见表2。

加权平均后的均方根均位于两者之间，UT1-UTC和ΔLOD的均方根比较靠近VLBI结果的均方根，是由于VLBI的解算结果精度高、权重大。

将2005年9月以及2006年1月、2月IERS 08C04的结果和VLBI的解算结果采用分段三次厄米多项式插值得到每天UTC12:00:00的参数值，再计算GPS解算结果和VLBI解算结果与IERS 08C04结果相比较的均方根，以均方根的平方的倒数作为权进行加权平均，公式如下：

式中，ERP_G、ERP_V和ERP_G+V分别表示利用GPS和VLBI结果插值后的序列及加权平均后的序列；RMS_V、RMS_G和RMS_G+V分别表示以上序列与ERP_IERS的均方根。计算后具体的差值见图3。

图 3 单独解算与加权平均后ERP和ΔLOD与IERS 08C04比较 Fig. 3 Comparison of ERP and ΔLOD Estimated from VLBI and GPS and GPS + VLBI Combined Solutions to IERS 08C04

图选项

最终计算得极移在X方向与Y方向上的分量与IERS 08C04公布值的均方根(见表3)。由表3的数据可以看出，经过加权平均后，极移在X方向上的分量X_p、极移在Y方向上的分量Y_p、日长变化ΔLOD各参数的均方根(RMS_X_p、RMS_Y_p、RMS_ΔLOD)有所下降，UT1-UTC的均方根(RMS_UT1-UTC)相当，总的来说， 两者联合提高了解算结果的可靠性，而且与IERS 08C04公布值的差值的波动也更加平稳，即解算结果的稳定性也有所提高。

本文利用GAMIT软件处理了2005年9月和2006年1月、2月的GPS数据，并用OCCAM 6.2软件处理了同时期的VLBI观测数据，解算的精度较高。 在解算出ERP和ΔLOD的基础上，本文运用分段三次厄米多项式的插值，将IERS 08C04的公布值和VLBI的解算结果插值到与GPS结果相同的时刻，IERS 08C04的数据密度大、精度高，所以插值效果好。VLBI数据的解算结果插值后精度略微有所下降，拉格朗日插值的结果和分段三次厄米多项式插值的结果的均方根差别较大，表明插值方法对其结果有较大影响。

对GPS和VLBI数据的解算结果进行加权平均综合后，ERP和ΔLOD的解算结果与IERS 08C04的公布值对比的均方根变小，差值的波动幅度变小，可靠性和稳定性都有所提高。