邵 璠124，王小亚1234

(1.中国科学院 上海天文台，上海 200030；2.中国科学院大学，北京 100049；3.中国科学院 上海天文台上海市空间导航与定位技术重点实验室，上海 200030；4.宇航动力学国家重点实验室，西安 710043)

1 引言

与地球相关的参考架原点主要有3个：包括海洋和大气的整个地球的质量中心(center of mass,CM)、不包括质量负荷的固体地球的质量中心(center of Earth,CE)和固体地球外表面的形状中心(center of figure,CF)。CM通常被空间大地测量所采用，因为卫星动力学是针对CM的，理想地球参考架定义要求地球参考架原点定义在CM。CE应用于某些地球物理学理论研究，如勒夫数的计算等。CF通常应用于与地面测量相关的学科，这里地面点间的几何形状是唯一可测量的量，所以一般应用于与地面大地测量和形变测量等有关的研究[1,2]。国际地球自转服务(International Earth Rotation and Reference System Service,IERS)发布的国际地球参考架(International Terrestrial Reference Frame,ITRF)是通过SLR、甚长基线干涉测量(Very Long Baseline Interferometry,VLBI)、全球定位系统(Global Positioning System,GPS)和星载多普勒定轨定位系统(Doppler orbitography and radio-positioning integrated by satellite,DORIS)四种空间大地测量技术建立和维持的，由各技术的全球观测网数据得到的技术内综合解，经由技术间再综合得到台站坐标和速度场来实现，其最初目标是要实现理想的地球参考架。由于四种空间大地测量技术台站都固定在地壳上，由这些台站坐标和速度场确定的ITRF并不能直接反映地球整体质量分布和地心运动，而是反映了由这些台站构成的多面体中心，即观测站网的中心(center of network,CN)，其近似于地球形状中心CF[1−4]。由于地球上海平面变化、冰川融化、大气环流、地幔对流和液核振荡等原因，地球质量会发生迁移，导致地球质心CM相对于地球形状中心(即参考架原点)发生位移，这个位移被称为地心运动[2]。地心运动必然会影响卫星精密定轨、地面点的精密定位、地球定向参数的确定、地球低阶引力场的确定以及所有以ITRF作为参考架的空间大地测量结果。因此，对地心运动的研究十分重要。

Cheng等人[5]利用5颗激光卫星的SLR数据，按照IERS2003规范，通过估算地球引力位一阶系数(采用的参考架为LPOD2005，在SLRF2005的基础上更新了一些测站坐标和速度)，研究了地心的周年运动。他们与GRACE得到的结果进行了比较，发现其X方向(2.5 mm)与ILRSX方向(2.7 mm)和GRACEX方向(3 mm)符合得很好，而Y和Z方向与ILRS和GRACE结果符合较差。Collilieux等人[6]以全球SLR数据解出的测站坐标作为类观测，进行参考架转换，得到平移参数，但是，认为该平移参数因为测站分布的不均匀性和稀疏性而不能被严格地当做地心运动，它们之间的差别定义为网形效应。国内也有一些学者对地心运动进行了研究，但地心运动的振幅与周期特性不完全一致。吴斌等人[7]认为地心变化幅度不超过5 cm，秦显平和杨元喜[8]发现地心运动存在近周年和半年的周期分量，但地心运动的长期变化并不明显；朱元兰和冯初刚[9]研究表明地心X、Y方向变化幅度为1 cm，Z方向为3∼4 cm；郭金运等人[10]发现地心运动存在长期和周期变化；李九龙[11]研究表明地心运动主要呈现出季节性周期变化，其中以周年项变化为主，半周年项的变化不够明显。究其原因，很大可能是基于老的IERS规范和针对单一卫星进行解算得到的地心运动并不可靠。我们利用上海天文台SHORD II软件，处理了Lageos-1,Lageos-2,Etalon-1和Etalon-2共4颗卫星从1993年1月―2017年12月的25年SLR数据，得到上海天文台(Shanghai Astronomical Observatory,SHAO)周解；并与ILRS各分析中心ACs的周解进行技术内综合，获得几何法的地心运动序列；又利用SHORD II软件，直接求解地心运动序列，并将这两类解的结果与CSR动力法结果进行比较；最后分析和探讨了不同地心运动确定方法结果的特征、可靠性和影响因素。

2 地心运动的计算方法

目前主要有四种空间大地测量技术，分别为VLBI,SLR,GPS和DORIS。其中，VLBI技术是一种相对测量，其观测目标是遥远的河外射电源，它们并不围绕地球转动，且其数据处理采用的是几何方法，所以VLBI技术自身无法确定地球参考架的原点。GPS技术由于受轨道模型误差及其他因素如卫星钟差、整周模糊度、天线相位中心改正等相关性的影响，所以无法给出精确的地心运动变化，同时，GPS卫星轨道较高且其数据处理经常采用的差分法也降低了对地心变化的敏感性。DORIS技术由于受自身系统误差较严重和在轨观测卫星数变化影响，不能精确地解算出地球质心的位置[11]。SLR技术是目前最精密的地球参考架原点确定技术，其观测的目标大都是面质比小、结构简单的地球动力学卫星，且大部分激光卫星离地球比较近，对地心运动很敏感，这为监测地心运动提供了有利条件。

目前，采用SLR技术确定地心运动的方法有三种：几何法、直接法和动力学法。

2.1 几何法

几何法就是先求解测站坐标序列，然后再求解地心运动序列[8]。该方法先利用卫星精密定轨得到测站在地心坐标系下的坐标序列，然后选择国际地球参考框架ITRF(或其他较权威的地球参考框架)为参考基准，由HELMERT 7参数转换到这个ITRF地球参考框架，得出测站坐标相对于该ITRF地球参考架的平移量。该平移量即为地心运动时间序列。

其中，T1,T2,T3是3个平移参数，也是所要求取的地心运动3个分量；R1,R2,R3为3个旋转参数；D为尺度因子。

几何法受到网形的影响，因此，为了使计算的地心运动序列更加真实可靠，我们应均匀选取测站。我们在Pavlis和Luceri[12]给出的SLR核心站基础上，选择性地剔除数据不好的站点(转换到新的参考架后，与框架结果的差大于0.3 m)作为SLR核心站，见表1。通过这些核心站进行7参数转换，得到的其中3个平移参数即是地心运动序列的3个分量。由于该方法与核心站的选取、数目和测站分布有关，所以其可靠性会受这些因素的影响。

表1 SLR核心站列表

我们采用了几何法来确定地心运动序列。我们先通过对Lageos-1,Lageos-2,Etalon-1和Etalon-2这4颗卫星进行弧长为7 d的松约束(站坐标先验标准差为1 m，极移先验标准差为20 mas，LOD先验标准差为2 ms)定轨，得到SHAO周解，然后对SHAO周解与ILRS各分析中心的周解进行加权综合，得到测站坐标，具体综合过程如图1所示。对于SINEX解中存在的松约束不予消除，我们进行直接综合，然后利用方差分量估算法对其进行加权综合，这样得到的综合周解还是松约束解[13,14]。我们再对混合后的周解相对于SLRF2008和SLRF2014分别进行7参数转换求平移参数，得到了两个地心运动序列。为了在下文图表中表述得清晰简洁，我们将其分别命名为SLRF2008几何法和SLRF2014几何法。

图1 SLR技术内综合流程图

2.2 动力学法

动力学法是通过处理SLR数据，利用动力学法进行卫星精密定轨来估算地球引力位一阶系数，从而确定地心运动序列[9]。地球引力场是在地固坐标系中以球谐函数展开的，如果坐标原点就是地球质心，那么引力位的一阶项系数C11,S11,C10将全部为0。但实际上，地球质心一直在运动，也就是说C11,S11,C10并不为0，动力学法就是在卫星精密定轨过程中，将地球引力位一阶系数作为待估参数，通过解算引力位一阶系数来确定地球的质心运动。地球引力位的球谐系数的级数展开形式一般为：

式中，G为地球引力常数，M为地球质量，r为地心向径，R为地球半径，ϕ为地心纬度，λ为地心经度，Pnm为缔合勒让德函数，Cnm和Snm为位系数，n和m为位系数的阶和次，Cnm和Snm的大小反映了地球内部的物质分布情况，其定义为：

式中，当m=0时δnm=1，否则δnm=0；dM表示坐标为(r,ϕ,λ)的空间点处的质量元。

在大多数情况下，大气、海洋、陆地水等都近似处理为地球表面质量负荷变化，这样式(4)和式(5)的时变部分可以简化为[11]：

式中，L(ϕ,λ,t)为地球表面上随时间t变化的一点，dS为该处的质量元，式(3)中没有一次项，这是假定CM在地球参考架原点，当存在地心运动时，其定义为：

式中，x,y,z为地球表面一点的三维坐标。如只考虑地球表面质量负荷L(ϕ,λ,t)变化引起的地球质心变化，则上述方程可重写为：

式中，∆x,∆y,∆z为地球质心运动的3个分量，C11,S11,C10为地球引力位一阶系数。

2.3 直接法

直接法就是在卫星精密定轨过程中将地心运动作为待估参数，在测量模型中加入地心运动改正项，通过循环迭代求得地心的位移量[9]。台站坐标定义在形心坐标系下，卫星轨道定义在质心坐标下，由于地球质心相对于形状中心存在地心运动，因此激光测距观测值须增加这一项改正。地心运动的具体改正方程为：

该方法的一个本质缺陷就是其地心运动估值与SLR测量误差中的偏心改正强相关，SLR测站偏心改正公式为：

其中，∆ρRO为测站偏心改正；ρ为卫星到测站的距离，ρx,ρy,ρz分别为ρ在测站坐标系中的3个分量；ECE为测量设备的测量中心相对于测站标定坐标在东向的偏心量；ECN为测量设备的测量中心相对于测站标定坐标在北向的偏心量；ECU为测量设备的测量中心相对于测站标定坐标在高程方向的偏心量。

在SLR数据处理中，通常不计算测站偏心改正，从式(11)和(12)可以看出，其公式非常相似，即地心改正实际估算了所有测站的共同偏心改正，测站的偏心改正误差也在其中。

本文分别利用Lageos-1和Lageos-2的SLR数据，采用直接法进行了地心运动解算，为了表述清晰简洁，其结果分别命名为Lageos-1直接法和Lageos-2直接法。直接法计算地心运动所采用的模型和解算策略如表2所示。

表2 SLR直接法计算地心运动所采用的模型和策略

3 地心运动特征分析

对获得的地心序列首先进行了傅里叶频谱分析，分析其中所含周期，然后假设地心运动同时具有长期项和周期项，则地心运动的时间序列可以表示为：

其中，c为长期变化中的常数项，d为长期变化中的系数项，k为周期函数的个数，Ah为对应第h个周期项的振幅，fh为对应的频率，θh为对应相位。将信号中的周期项展开可得：

将求解的地心运动∆Rj当作已知观测量，代入式(13)中，则可得到误差方程：

式中，V为地心运动∆Rj的残差向量，X为未知数向量(即所求的周期项和长期项系数)，B为系数矩阵，L为由∆Rj组成的向量。通过最小二乘法可得未知数向量的估值为：

其中，P为地心运动的权矩阵，此处取单位权。当求得参数向量X后，可通过式(15)和式(16)得到周期项的振幅和相位，这样就可以获得地心运动的长期变化和周期变化特征。

4 结果分析

本文利用上海天文台SHORD II软件处理了Lageos-1,Lageos-2,Etalon-1和Etalon-2这4颗卫星从1993年1月―2017年12月期间的SLR数据，将得到的SHAO周解与ILRS各分析中心(包括ASI,BKG,DGFI,ESA,GRGS,JCET和NSGF)提供的周解进行加权混合，由混合后的解进行7参数转换，分别得到以SLRF2008和SLRF2014作为参考的两个地心运动序列，并进行了傅里叶变换及长期和周期性特征分析，其结果见图2―4。图2显示了SLRF2008几何法确定的地心序列傅里叶频谱分析，证实了地心在X,Y,Z三个方向上有明显的周年项和半年项，周期分别为0.9996 a和0.4996 a，然后对地心序列中的周年项、半年项的振幅和相位以及长期项按照第3章方法进行最小二乘法拟合，这里只给出了基于SLRF2008的几何法傅里叶频谱分析结果，其他结果类似，并且周期也相同。图3和图4分别给出了SLRF2008几何法和SLRF2014几何法确定的地心运动及其特征拟合序列，从图中可以看出周期拟合序列较好地模拟了原始的地心运动序列，但是在振幅上的拟合结果还不是很好，特别是在Z方向。

利用Lageos-1和Lageos-2的SLR数据并通过直接法确定了地心运动，结果如图5和图6所示。图7给出了CSR利用动力学法确定的地心运动。图中，黑色曲线代表最小二乘法拟合的线性项和周期项之和。可以发现，地心运动在X,Y方向变化幅度均比较小，而在Z方向变化幅度较大，且通过直接法求得的地心序列比几何法和动力学法求得的变化幅度更大，更不规则，特别是在Z方向。这可能是因为直接法求解的质心运动与SLR测站偏心改正强相关，其可靠性相对较低。

图2 SLRF2008几何法确定的地心运动序列傅里叶频谱分析图

图3 SLRF2008几何法确定的地心运动序列及其特征拟合序列

图4 SLRF2014几何法确定的地心运动序列及其特征拟合序列

图5 Lageos-1直接法确定的地心运动序列

图6 Lageos-2直接法确定的地心运动序列

图7 CSR动力学法确定的地心运动序列

我们采用几何法和直接法计算了地心运动与CSR提供的动力学法地心运动序列的平均统计结果(如表3所示)，其中CSR的结果来自于ftp.csr.utexas.edu/pub/slr/geocenter发布的地心运动序列，其时间跨度为2002―2017年，解算频率为每30 d提供一组解，选择的参考架原点为SLRF2014原点。从表3中可以看出，利用直接法求得的地心运动平均值比几何法和动力学法的结果偏大，几何法求得的地心运动与CSR提供的动力法结果较接近。

表3 SLR不同方法监测地心运动的平均统计结果 mm

表4给出了本文和基于CSR提供的动力学地心运动序列分别拟合得到的地心运动长期变化率。从表中可以看出，利用参考SLRF2014的几何法确定的地心长期运动速率比利用参考SLRF2008几何法确定的结果明显小，特别是Z分量，说明地球参考架SLRF2014相对于SLRF2008定义的原点更加稳定。利用几何法和直接法求得的地心运动速率Z分量与CSR结果存在差异，其原因可能是几何法求得结果受网形影响，由于核心站南北分布不均匀，导致解算的地心运动Z分量可靠性降低。而直接法使用的是单星解，卫星观测几何强度不够，解算的地心运动Z分量存在误差，但从整体结果来看，本文计算的地心运动速率与CSR提供的速率较一致，都比较小，说明SLR确定的地心比较稳定，其长期变化并不明显。其中，基于SLRF2014几何法确定的地心稳定性好于0.1 mm·a−1，其他方法确定的地心稳定性差于 0.1 mm·a−1。

表 4 地心运动长期变化率 mm·a−1

本文对地心运动在X,Y,Z三个方向上的周年及半年项振幅和相位(相位的参考历元为1月1日)进行了解算，并与CSR的结果进行对比(如表5所示)，CSR结果来自对其提供的地心运动序列进行类似的周年项及半年项最小二乘法拟合。从表5可以看出，使用Lageos-1直接法与Lageos-2直接法确定的地心运动序列，其半年项比较接近，但是周年项相位相差较大，这一方面可能与该方法参数之间的相关性较强有关，也可能与单颗卫星解算地心运动时的观测几何强度不够有关。通过比较还发现SLRF2014几何法与SLRF2008几何法得到的地心运动较接近，其中SLRF2014几何法得到的周年项和半年项与CSR动力学法结果更接近，这说明通过对SHAO提供的SLR周解与ILRS各分析中心提供的周解进行技术内加权综合后，可以得到较可靠的地心运动序列。

表5 地心运动周年项的振幅和相位与半年项的比较

5 结论

本文处理了Lageos-1,Lageos-2,Etalon-1和Etalon-2这4颗卫星在1993―2017年期间的全球SLR观测数据，分别通过几何法(参考SLRF2008和SLRF2014)和直接法(针对Lageos-1和Lageos-2)确定了SLR地心运动序列，利用傅里叶变换和最小二乘法分析了两种方法解算的地心运动的长期变化率、周年项和半年项振幅和相位，并与CSR提供的动力法结果进行了比较，发现：(1)地心运动的长期项比周期项小1个量级；(2)与SLRF2008几何法相比，SLRF2014几何法得到的地心运动长期项明显减小，特别是在Z方向，且其得到的地心运动周年项和半年项与CSR动力学法得到的结果很相似，说明综合SHAO与ILRS ACs周解来确定地心运动是可行的；(3)将Lageos-1直接法与Lageos-2直接法确定的地心运动序列相比，其半年项有较好的一致性，但周年项相位相差较大，与几何法及CSR提供的动力学结果在Z方向差异也很大，这可能与该方法解算的地心运动与测站偏心改正强相关有关，也可能与单颗卫星解算地心运动几何结构不够好有关，该方法的可靠性有待进一步提高。