在漫长的时间里，地球表面液体层的质量是守恒的，其在地表的再分布周期从1 d~1 a不等，与地球的自转(来自太阳的热量驱使大气)和公转(每年的全球水循环)相关^[1]。对于一个弹性地球来说，地表液体层的质量再分布势必会引起固体地球的形变^[2]。从全球范围来看，如果地表液体层由于季节性循环主要集中在北半球，导致北半球的负载增加，则北半球被挤压，南半球扩张，反之亦然，这就是季节性地球形变的原因。

目前,对于弹性地球形变的研究，传统上习惯使用Farrell的表面负荷理论即负载勒夫数理论。该理论广泛应用于大气层和水体等负载模型的理论建模中。1972年，Farrell首次对弹性地球进行建模，阐述了地球的形变响应是由地表负载所引起的引力位变化的球谐分量的函数^[3]。随着地球模型研究的不断深入和完善，SNREI模型逐渐成为人们普遍接受的“标准地球模型”。在这个模型下，地球是对称、非旋转、完全弹性和各向同性的，当其表面负载发生变化的时候，固体地球会产生静态形变^{[4, 5]}。Blewitt首次利用地表引力位一阶球谐分量的变化和负载勒夫数理论，通过GPS数据反演得到了季节性的地球形变^[6]。孙付平等利用GPS和VLBI组合的站坐标、站速度，采用Delaunay算法生成的三角形来逼近地形表面，研究了地球的表面积、体积及它们的变化趋势。结果表明，目前地球半径每年大约缩短3~4 mm，地球整体上处于压缩变形中^{[7, 8]}。申文斌等基于ITRF2000框架的671个台站，利用同样的方法研究得到地球目前正处于南胀北缩的状态中，但地球整体处于膨胀状态^[9]。此外，夏一飞等研究了日月引潮力以及变化着的地球自转惯性离心力的作用下所产生的地球周期性弹性形变^[10]；徐建桥等研究了地球固体部分对液核动力学效应引发的核幔边界和内核边界上压力和引力扰动的形变响应^[11]。 1 弹性地球形变与地球参考框架原点

国际地球参考框架的原点定义在整个地球系统的质量中心(CM)，由于地球整体在惯性空间下不受外力作用，所以CM的位置并不会因为地球表面的质量再分布而产生变化，变化的是固体地球质量中心(CE)。虽然参考框架的原点定义在CM，但实际上却是通过与CE相差只有2%的固体地球外表面的几何中心(CF)来实现的^{[1, 2]}地表质量再分布，CF框架点是由地球表面覆盖的统一的、无限多的点构成的，这些点的运动也需要考虑。实际上可以通过在全球建立合适密度的观测站来实现，使用GPS易于实现。

地表液体层质量再分布的过程会通过两种方式影响到大地测量网络：1)伴随地表质量再平衡的过程所产生的CE相对于CM的位移，即地心运动，其对ms级地球参考框架所造成的影响是不容忽视的^[12]；2)地表跟踪站本身伴随着固体地球形变而产生的位移^{[1, 12]}。基于此，如果能够得到地表跟踪站的精确位置变化，应用相应的弹性地球负载模型，就可以探索固体地球的弹性形变机制。目前的研究表明，地表负载的一阶变化项与地表负载的质量中心和地心直接相关，而高阶项与负载的质量中心之间的关系尚不明确，所以本文只讨论一阶负载下的地球弹性形变。 2 负载力矩理论及其与弹性地球形变的关系 2.1 负载力矩的基本理论

由于整个地球系统在惯性空间下的运动是不受外力的，即满足动量守恒的条件，因此伴随地表质量再分布的过程有：

式中，Δ_L是地表液体层的质量中心在CE框架下的位置变化；M_L为地表液体层的质量；M_E为固体地球的质量；则CE相对于CM的变化为：

负载力矩矢量被定义为：

应用“负载力矩”理论可以更加方便地研究地表负载变化驱使固体地球形变的机制。 2.2 负载力矩理论反演弹性地球形变的方法

由于地表质量再分布所引起的CE相对于CM的位置变化会导致地球表面引力位发生变化，其一阶分量的球谐函数可以表示为：

式中，Ω为地理位置(纬度φ，经度λ)；(Ω)是径向单位向量； g 为地表重力加速度。式(2)是一个纯粹的球谐函数，可以直接应用负载勒夫数理论。根据负载勒夫数理论，CE框架下地球表面位移ΔS_h(Ω)(径向)，ΔS_l(Ω)(切向)可以表示为： <

式中，是切向单位向量；地表梯度算子；此处使用Farrell的负载勒夫数模型^[4]：

理论上，CF框架原点满足地球表面位移场的积分为0的条件，这与沿着各个轴向的投影的无网络变换是一致的。CE框架的最大不足之处在于它不可以被直接观测。而CF框架与CE框架特别近似，相比之下，CF框架更方便于实际观测^[13],CE框架更适合于理论研究。若将CF框架与CE框架进行严格比较，则需要一个从CE框架到CF框架的坐标转换。根据Blewitt的负载勒夫数转换理论，CF框架下的负载勒夫数与CE框架下的负载勒夫数之间的关系为^[2]：

CF框架下地表某处点位位移(为了看起来更简洁，这里省略了Ω)：

因此，CF框架下地表GPS测站的位移 [Δs]_CF=(Δs_x,Δs_y,Δs_z)^T可以表示为：

式中，G 是3×3维的将地心坐标转换为地面(东西、北南、垂直)坐标的雅可比矩阵：

事实上，由于框架点全球分布的不均匀性，地球参考框架从本质上来说只是 CF框架的一个近似，实际上为CN(the center of network)框架，即由地面网实现的参考框架。利用GPS数据得到的[Δs]_CF是CN框架下测站的形变，而式(6)直接采用的是CF框架下的负载勒夫数，这个误差可以通过附加平移参数t来减弱或消除^[14]：

因此，通过计算全球范围内的GPS跟踪站的位移[Δ s]_OBS，就可以通过式(7)计算得到“地表负载力矩” m ，进而探索固体地球形变规律；相反，如果地表负载变化能够被建模，我们也可以通过式(7)计算得到由于弹性地球形变而造成的测站位移。

2.3 最小二乘法的设计矩阵

根据式(7)，其未知参数矩阵为：

设计矩阵为：

国际 GNSS服务(IGS)每星期为用户提供精确的IGS跟踪站坐标及变化率，这为利用IGS数据探索固体地球形变提供可能。截至目前，IGS共有11个分析中心，GPS卫星跟踪网中共有856个跟踪站^[15]。利用这些可靠的数据源，可以精确解算地表负载力矩，进而研究固体地球形变机制。本文选择IGS05框架的132个框架点(图 1)，利用IGS重新处理得到的近10 a的igl周解数据^[16](共522周)，反演得到了近10 a地球表面负载力矩变化的时间序列(图 2)。图 2中m_x、m_y、m_z的方向与WGS-84定义的XYZ三轴方向一致，igl解在数据处理时根据阈值剔除了质量不好的测站，较为完善地考虑了各种突变，同时，igl解准确合理地消除了由板块运动、冰后回弹等引起的测站位置的长期变化，剔除了仪器变动和同震位移等引起的测站位置突变。

图 1 IGS05框架的132个框架点Fig. 1 132 Sites of the IGS05 Reference Frame

图选项

图 2 近十年地表负载力矩变化 Fig. 2 The Variations of Load Moment in Recent Ten Years

图选项

为了计算负载力矩时间序列能量谱密度并确定其主要周期，对m_x、m_y、m_z序列分别进行傅里叶变换，得到其功率谱(图 3)。功率谱密度越大，说明该频率对负载力矩分量的贡献就越大，因此可由图 3确定负载力矩各分量的主频率。

图 3 地表负载力矩时间序列的功率谱分析Fig. 3 The Power Spectrum Analysis of Load Moment Time Series

图选项

计算分析表明：负载力矩各分量都有明显的年周期变化和相对较弱的长周期变化，m_x方向有一个1.3 a和一个2.1 a的周期，m_y和m_z均有一个2 a周期。由于反演的数据源采用的是周解，相当于采样率为一周，因此还不足以发现更短的周期。 4.2 负载力矩时间序列的回归分析

本文负载力矩时间序列m(t)的回归模型定义为：

式中，a₀为偏差；a′₀为变化趋势项；ω_n为不同周期变化的角频率；a_n和φ_n为相应周期变化的振幅和相位；Δt=t_i－t₀，本文中t₀=2 000.0。

根据式(8)，计算得到三个方向负载力矩时间序列的周年项振幅和相位，并将本文估计结果和Blewitt的实验结果分别列于表 1。

根据负载力矩各方向的周年相位，我们可以大致推断出负载力矩每年在地表的运动轨迹(图 4)。分析表 1中m_z方向的相位可知，北半球地表负载力矩的最大值出现在每年的2、3月份，而南半球出现在8、9月份；结合表 1和图 4可知，地表负载力矩的最小值分别于5月和11月出现在古巴和印度尼西亚附近；每年的6~11月地表负载力矩主要集中在南半球，而12~5月主要集中在北半球，其中12月份出现在中国境内，这个时候地表负载力矩对中国范围的影响相对来说最为明显；此外，无论南北极点，其地表负载力矩最大值均出现在各自的冬季，中国也是如此。

图 4 每年负载力矩在地球表面的运动轨迹Fig. 4 Annual Trajectory of the Load Moment on the Earth Surface

图选项

为了更加深入研究地表负载力矩对中国范围 内测站的影响，选择09年12月的地表负载力矩时间序列，利用式(7)反算中国境内6个IGS站(BJFS、SHAO、WUHN、LHAZ、URUM、KUNM)的位移，结果列于表 2。

由表 2可知,由于地表负载力矩的作用，中国范围内的IGS站在12月份会下降2 mm左右，同 时伴有1 mm左右的水平位移。 4.4 分析与讨论

与Blewitt实验结果相比，本文估计的m_x方向相位比Blewitt的结果要小，振幅与Blewitt的结果基本一致； m_y和m_z方向相位与Blewitt的结果符合较好，振幅均比Blewitt的结果小得多。

分析认为，造成与Blewitt结果振幅相差较大的原因是Blewitt反演的数据源时间跨度为1995.642~ 2001.046，要早于本文数据源2000.0~2010.0，且时间跨度较短；本文所选用的测站坐标统一归算到ITRF05框架，其数据处理所采用模型的精度及接收机质量较Blewitt所使用的ITRF97框架的站点均有显著提高；此外，Blewitt选取的站点是3个IGS分析中心的66个GPS跟踪站(图 5)，而本文选取的站点是IGS05框架的132个框架点(图 1)，测站分布更广，更均匀，尤其在海洋上也有一定的观测资料，较之Blewitt的数据源更加稳定可靠。最后，Blewitt的模型中没有加入平移参数来削弱CN框架和CF框架之间的误差。

图 5 Blewitt实验数据源的测站分布Fig. 5 The Distribution of Blewitt’s Experimental Data Source

图选项

本文利用IGS数据反演得到的季节性固体地球形变机制与已知的地球物理知识基本相符。研究结果显示，每年的2、3月份北半球地表负载最大，而8、9月份南半球地表负载最大；12月份地表负载力矩出现在中国境内，这将会造成中国范围的测站在水平和垂直方向分别产生1 mm和 2 mm左右的位移。因此，季节性固体地球形变对我国mm级CGCS2000框架的影响是不容忽视的。

地表液体层的质量再分布与地心运动密切相关，因此深入研究季节性固体地球形变对于季节性地心变化的研究具有十分重要的借鉴意义。

随着IGS服务的不断完善和发展，GPS对固体地球物理过程研究的意义越来越重大，然而GPS也有自身的缺陷，即IGS其跟踪站主要集中在陆地，占地球面积70%的海洋却缺少足够的观测资料；此外GPS跟踪站在陆地上的分布南北半球也不均匀，这会对单独使用GPS数据反演地表负载力矩的精度带来影响。如何将海洋数据和GPS数据融合来研究地表负载力矩，是下一步本文需要研究的内容。