随着卫星技术和数据处理方法的进步，空间大地测量技术已经具备了以高时间分辨率和高空间分辨率对由地表质量重新分布导致的地球形状变化和重力场变化进行高精度监测的能力(Blewitt et al.，2001; Tapley et al.，2004). IGS(International GNSS Service)自1994年成立以来，已提供了全球400多站近20年的地表形变时间序列.GRACE(Gravity Recovery and Climate Experiment)计划实施以来，利用其提供的时变重力场研究全球质量变化受到了广泛关注，成为国内外地学界的一个研究热点(Tapley et al.，2004; 胡小工等，2006).

现有研究成果表明，由GPS和GRACE数据获得的季节性地表形变存在差异，但随着数据处理方法的不断精化，这些差异在不断缩小(Van Dam et al.，2007; Tregoning et al.，2009; Tesmer et al.，2011; 张良镜等，2012).GPS通过精确确定地面点的位置监测地表形变，是一种点测量模式，容易受局部负载的影响(Nahmani et al.，2012)；此外技术相关误差(Ray et al.，2008)和误差改正模型的不完善(Penna et al.，2007)等也会产生一定的影响.而GRACE监测到的是一定区域内的，经过平滑的地表质量变化的平均值，是一种面测量模式，空间分辨率比GPS低，受滤波半径等的影响(Schrama et al.，2007)；而且由GRACE重力场球谐系数计算得到的地表形变不受温度影响，因为GRACE仅对地表质量变化敏感.另外，目前对GPS和GRACE季节性地表形变的一致性研究主要集中在高程方向，对水平方向的关注较少，仅有部分学者针对特定区域比较了GPS/GRACE水平方向的一致性(Fu et al.，2013; 王林松等，2014).水平方向的负载形变大致是高程方向的1/3，信号相对较弱，容易受观测噪声影响(Van Dam et al.，2002)；而且水平方向的构造信号显著，对周期性信号的提取也有影响.但考虑到GPS水平方向的定位精度(~3 mm)明显优于高程方向的(~6 mm)，因此如何有效地利用这一特征对全球范围内的基于GPS和GRACE数据的三维地表形变进行全面的比较具有十分重要的意义.

基于以上背景，本文旨在分析GPS/GRACE的三维季节性地表形变的差异，进而探讨其原因，明确两种技术监测地表形变的能力，以期从空间大地测量数据中获取更加准确可靠的季节性地表质量变化信息.

GPS站地表形变采用IGN(Institut Géographique National)提供的ITRF2008-GPS站坐标残差时间序列(GPS^ITRF08)，由CATREF软件处理得到(Altamimi et al.，2011).它是ITRF2008数据处理的副产品，已消除了线性速度和测站跳变.ITRF2008是目前最新版本的地球参考框架，输入数据是第一次重新处理的IG1周解.IG1周解是对各分析中心基于最新的数据处理策略和改正模型进行统一处理的解进行组合得到的(Steigenberger et al.，2006).选取的测站为IGb08的232个框架站.

采用CSR(Center for Space Research)提供的RL05版本的重力场球谐系数计算地表形变，截止至60阶(GRC^CSR).首先利用GAC系数恢复大气和海洋的高频非潮汐变化.由于GRACE数据处理采用的是CM(Center-of-Mass of the Earth plus load)框架，需要利用Swenson等提供的一阶系数获得CF(Center-of-Figure of the Earth)框架下的1~60阶的重力场球谐系数(Swenson et al.，2008).二阶带谐项用CSR提供基于SLR的ΔC₂₀进行替换.最后根据Farrell的负荷理论(Farrell，1972)，可由CF框架下的负荷Love数计算测站的地表负载形变.

基于地球物理模型的地表形变由大气、海洋和陆地水模型计算得到的负荷形变之和得到(AOC)，采用QOCA软件计算.大气数据由NCEP(National Centers for Environmental Prediction)提供，全球格网的分辨率为2.5°×2.5°，采用反变大气压假设计算海洋上方大气变化导致的地表形变(Van Dam and Wahr，1987).海底压强数据由ECCO(Estimating the Circulation & Climate of the Ocean)提供，该模型基于MIT的全球环流模型.覆盖区域从79.5°S到78.5°N，除了赤道附近(20°S到20°N)纬度方向的分辨率为0.3°外，全球格网的分辨率为1°×1°.陆地水模型由GLDAS(Global L and Data Assimilation System)模型提供，该模型以全球1°×1°分辨率提供雪、土壤湿度和植物树冠的质量变化.

利用公式(1)计算由GPS、GRACE和地球物理模型获得的三维地表形变的季节性变化：

其中w=2π/365.25，A_a和A_sa分别是周年和半周年振幅，Φ_a和Φ_sa分别是周年和半周年相位，单位为年积日doy.

图 1给出的是由GPS^ITRF08(红色)、GRC^CSR(蓝色)和AOC(绿色)计算得到的高程方向的周年变化.由图可知，除南极和南美南部外，GPS^ITRF08与GRC^CSR的周年相位符合得比较好.有65%的测站GPS^ITRF08与GRC^CSR的U方向的周年相位差小于30天.GPS^ITRF08的周年振幅总体大于GRC^CSR的，特别是在南极、南美南部和部分岛屿.GPS^ITRF08与GRC^CSR周年振幅差值均值为0.44±1.65 mm.由图 1还可看出，与GPS^ITRF08相比，GRC^CSR与AOC的周年振幅和周年相位的符合都更好.这一现象是合理的，因为GRACE和地球物理模型反映的是大尺度的质量迁移，而GPS站地表形变还对区域性质量变化等因素敏感.

图 1

Fig. 1

图 1 由GPS(红色)、GRACE(蓝色)和地球物理模型(绿色)计算得到的高程方向的周年变化 (箭头长度代表周年振幅(mm)，箭头方向代表周年相位(从正东方向逆时针，相位起点为一月份)) Fig. 1 Annual vertical variations from GPS (red), GRACE (blue) and loading models (green). The length of arrows shows the annual amplitude (mm) and the direction of arrows shows the annual phase (anti-clockwise from the east)

图 2和图 3给出的是由GPS^ITRF08(红色)、GRC^CSR(蓝色)和AOC(绿色)计算得到的水平方向的周年变化.总体而言，GPS^ITRF08和GRC^CSR水平方向的周年变化一致性不如高程方向，但GRC^CSR和AOC的水平方向周年振幅和相位仍然符合得比较好.GPS^ITRF08的EN方向的周年振幅比GRC^CSR大很多：GPS^ITRF08的EN方向周年振幅均值分别为1 mm和1.3 mm，而GRC^CSR分别为0.4 mm和0.4 mm.GPS^ITRF08与GRC^CSR的EN方向的周年相位差小于30天的测站比例分别只有32%和22%，近一半测站EN方向的周年相位差大于60天(E方向为48%，N方向为51%).除欧洲区域和部分高纬区域外，GPS^ITRF08的N方向的最大振幅都出现在一季度.在欧洲区域，GRC^CSR和AOC的E方向呈现出比其他区域更加明显的周年变化，而GPS^ITRF08的E方向并没有呈现出这一特征.相对于其他区域，欧洲区域的GPS^ITRF08在EN方向的周年相位显得杂乱无章.

图 2

Fig. 2

图 2 由GPS(红色)、GRACE(蓝色)和地球物理模型(绿色)计算得到的East方向周年变化.(箭头大小和方向的含义同图1) Fig. 2 Annual east variations from GPS (red), GRACE (blue) and loading models (green). The length and direction of the arrows have the same meaning as Fig.1

图 3

Fig. 3

图 3 由GPS(红色)、GRACE(蓝色)和地球物理模型(绿色)计算得到的North方向周年变化.(箭头大小和方向的含义同图1) Fig. 3 Annual north variations from GPS (red), GRACE (blue) and loading models (green). The length and direction of the arrows have the same meaning as Fig.1

综上所述，GPS^ITRF08和GRC^CSR的周年变化在高程方向上具有较好的一致性，但水平方向的差异明显，GPS^ITRF08水平方向的周年变化明显比GRC^CSR大.无论是水平方向还是高程方向，GRC^CSR和AOC在反映大尺度负载形变上具有较好的一致性.

为了进一步分析GPS^ITRF08和GRC^CSR季节性地表形变差异的原因，分别对GPS和GRACE技术确定地表形变的不确定度进行了分析.首先将JPL提供的GPS站时间序列(用GPS^JPL表示)和GPS^ITRF08进行对比，分析不同GPS数据处理策略和处理软件计算的地表形变的不确定度.GPS^JPL的数据处理流程与GPS^ITRF08不同，具体为：采用GIPSY软件，利用JPL无基准约束的卫星轨道和钟差，对所有测站进行PPP定位，然后再利用七参数转换，转换到IGS框架下.比较结果表明，GPS^ITRF08与GPS^JPL在ENU方向的周年相位差小于30天的测站比例分别为60%，67%和84%.54%测站E方向的周年振幅差异小于0.3 mm，77%测站E方向的周年振幅差异小于0.5 mm.65%测站N方向的周年振幅差异在0~1 mm之间，差值均值为0.5 mm，如图 4所示.具体原因有待进一步研究.78%测站U方向的周年振幅差异小于1 mm.

图 4

Fig. 4

图 4 GPSITRF08与GPSJPL 在N方向周年振幅差异的直方图 Fig. 4 Histogram of annual amplitude difference between GPSITRF08 and GPSJPLin the north direction

同样地，将基于JPL(Jet Propulsion Laboratory)和GFZ(GeoForschungsZentrum)提供的GRACE重力场系数计算的地表形变(分别用GRC^JPL和GRC^GFZ表示)与GRC^CSR进行比较.除GFZ的ΔC₂₀不需要进行替换外，由JPL和GFZ重力场系数计算地表形变的处理流程参见2.2节.比较结果表明，三种GRACE形变序列EN方向的周年振幅和周年相位都符合得非常好.U方向的差异相对水平方向稍大，主要是在高纬度地区.以GRC^CSR和GRC^JPL为例，94%的测站U方向的周年相位差小于20天，85%的测站周年振幅差小于0.3 mm.

表1(Table 1)

表 1 基于GPS与GRACE数据的地表形变ENU方向的周年振幅和周年相位差值均值及标准差 Table 1 Average and STD of the difference of annual variations between GPS and GRACE in ENU direction Amplitude (mm) Phase (doy) E N U E N U GPSITRF08/GRCCSR 0.61±0.69 0.88±0.69 0.44±1.65 -5±89 -7±84 -5±49 GPSITRF08/GPSJPL 0.09±0.52 0.40±0.55 0.01±0.93 -10±46 -12±49 -6±29 GRCCSR/GRCJPL 0.00±0.02 0.02±0.03 0.12±0.38 1±5 1±5 -3±10 GRCCSR/GRCGFZ 0.00±0.03 0.01±0.04 -0.08±0.28 -1±7 -2±6 -6±12

表 1 基于GPS与GRACE数据的地表形变ENU方向的周年振幅和周年相位差值均值及标准差 Table 1 Average and STD of the difference of annual variations between GPS and GRACE in ENU direction

综上所述，不同GRACE形变序列间的一致性很好.不同GPS形变序列间的不确定性比不同GRACE形变序列间的不确定性要大的多.不同GPS形变序列间水平方向周年相位的不确定度比高程方向的大.因此我们推断GPS形变本身的不确定度是导致GPS/GRACE季节性地表形变差异的一个重要原因，特别是水平方向.

GPS站的季节性地表形变除受地表质量重新分布的影响外，还受与地表温度变化相关的热弹性形变的影响.而GRACE重力场系数仅对质量变化敏感，不受温度影响.Fang等建立了由全球温度变化导致的地表形变的模型，包括高程方向和水平方向(Fang et al.，2014).而之前学者大多基于半空间模型假设研究热弹性形变，主要研究的是高程方向或区域内的水平变化，并没有给出全球三维热弹性形变的计算模型(Berger，1975; Dong et al.，2002; Tsai，2011).本节采用Fang等的研究结果来分析热弹性形变对GPS/GRACE地表形变差异的影响.Fang等指出，与温度相关的热弹性形变最大可导致高程方向~2 mm，水平方向~1 mm的形变.E方向受热弹性形变的影响总体较小，只有欧洲和北美西海岸较为显著，在此暂不详细讨论E方向.全球地表质量迁移与温度变化是密切相关的，夏季北半球质量负载最小，地表向上膨胀；与此同时，北半球温度最高，地表也是向上膨胀的.因此质量负载导致的形变和热弹性形变方向在高程方向是基本一致的，由此可以解释GPS站观测到的形变(包括负载形变和热弹性形变)比GRACE的大.

下面主要分析热弹性形变对GPS/GRACE的N方向差异的影响.图 6给出了GPS^ITRF08和GRC^CSR之差(用GPS^ITRF08-GRC^CSR表示)的N方向的周年变化.除欧洲区域和部分海岛外，两者间差异的变化规律与Fang等给出的热弹性形变分布(参见Fang等图 3)具有较好的空间相似性，初步推测热弹性形变是导致GPS^ITRF08-GRC^CSR的N方向形变差异的重要原因之一.具体表现为：除南北两极外，全球范围内GPS^ITRF08-GRC^CSR的N方向周年相位基本都在冬季，与Fang等给出的N方向的热弹性形变的运动方向是一致的.因为冬季北半球温度低，南半球温度高，受温度梯度影响，地表由南向北运动；夏季N-S方向的运动正好相反.GPS^ITRF08-GRC^CSR的N方向周年振幅在青藏高原区域较大，澳大利亚由北至南呈条带状，这些都与Fang等的结果一致.但是GPS^ITRF08-GRC^CSR的N方向周年振幅在海洋上和欧洲区域与Fang等的结果有较大的差异.究其原因，Fang等是利用全球温度变化场的球谐系数计算热弹性形变，但在计算全球温度变化球谐系数时，并未考虑海洋的温度变化的影响.也就是说部分海岛上测站受温度变化导致的热弹性形变并未考虑在内.另外，GPS^ITRF08-GRC^CSR在欧洲区域的变化也与Fang等的结果不一致，其原因将在4.3节予以讨论.

图 5

Fig. 5

图 5 GPSITRF08-GRCCSR的N方向的周年相位 Fig. 5 Annual phase variations of GPSITRF08-GRCCSR in the north direction

图 6

Fig. 6

图 6 GPSITRF08-GRCCSR的N方向的周年振幅 Fig. 6 Annual amplitude variations of GPSITRF08-GRCCSR in the north direction

另外需要注意的是，GPS^ITRF08-GRC^CSR的 N方向的周年振幅总体比Fang等给出的量级要大，可能是由GPS^ITRF08周年振幅略大(见图 4)引起的，具体原因将另文研究.

Wahr等提出综合利用GPS站水平和高程方向的地表形变大小和方向来判断负载所在的方向(Wahr et al.，2013).本节采用这种方法来区分大尺度负载和区域性负载引起的形变，进而分析区域负载对GPS/GRACE地表形变差异的影响，重点分析欧洲区域GPS^ITRF08-GRC^CSR与Fang等结果的差异的原因.图 7给出的是GPS^ITRF08-GRC^CSR的水平方向的周年振幅及指向(水平向右为东方向，垂直向上为北方向)，其中红色箭头与GRC^CSR高程方向的周年相位对应(与温度变化对应)；蓝色箭头与GPS^ITRF08-GRC^CSR高程方向的周年相位对应，具体计算方法为：首先计算出GPS^ITRF08-GRC^CSR的ENU方向的周年振幅和相位，然后根据指定的高程方向周年相位的时间点计算对应的EN方向的波峰和波谷的差值作为振幅，箭头大小和方向分别为EN方向的振幅大小之和与方向之和.

图 7

Fig. 7

图 7 GPSITRF08-GRCCSR水平方向周年振幅及指向 红色：与GRCCSR高程方向的周年相位对应；蓝色：与GPSITRF08-GRCCSR高程方向的周年相位对应. Fig. 7 Annual amplitude and direction of GPSITRF08-GRCCSR in the horizontals Red: w.r.t. the annual phase of GRCCSR; Blue: w.r.t. the annual phase of GPSITRF08-GRCCSR.

从图 7给出的与GRC^CSR高程周年相位对应的GPS^ITRF08-GRC^CSR的水平方向的振幅及指向(红色箭头)可以看出，欧洲区域的指向具有一定的一致性.当GRC^CSR高程振幅最大时(与欧洲区域的最高气温夏季对应)，GPS^ITRF08-GRC^CSR存在一个向东的系统性运动，与Fang等图 3给出的GPS站的热弹性形变运动规律一致.另外，从图 7给出的与GPS^ITRF08-GRC^CSR高程周年相位对应的GPS^ITRF08-GRC^CSR的水平方向的指向(蓝色箭头)可以看出，此时欧洲区域的指向是没有规律的，初步推测是由区域负载引起的，后续还需进行深入分析.除欧洲区域外，图 7的蓝色箭头在其他区域指向大致是一致的，即没有明显的区域负载的影响.因此相对于其他区域，欧洲区域的区域性形变更为显著，导致了图 3所示的GPS^ITRF08在EN方向的周年相位显得杂乱无章.

GPS^ITRF08和GRC^CSR地表形变的周年变化在高程方向上具有较好的一致性，但水平方向的差异明显.一方面是因为不同GPS站形变时间序列水平方向间的不确定度比高程方向的大：虽然GPS站水平方向的重复性比高程方向的好，但由于水平方向的线性运动显著、季节性信号比高程方向的小，计算得到的水平方向的周年变化受不同数据处理软件和模型的影响大，外符合精度差.

另一方面，GPS站除了能够监测GRACE反映的大尺度地表负载形变，还受热弹性形变和区域负载的影响.如前所述，由于全球温度变化与地表质量变化具有一定的同步性，这两种因素引起的高程方向的周年相位是基本相同的(例如夏季北半球由于温度高、地表负载最小，都是向上膨胀形变的).而且地表质量负载引起的高程方向的最大振幅~9 mm，平均~3 mm；而热弹性形变导致的高程方向最大振幅~2 mm.因此热弹性形变并未明显改变GPS站高程方向的周年相位(与GRACE相比)，但可使高程方向的周年振幅略微增加.水平方向则不同，热弹性形变导致的水平方向的形变大小与负载形变(平均~0.4 mm)相当，形变方向也有所差异，导致GPS/GRACE水平方向间的差异明显.

区域负载导致的水平方向的形变比高程方向更加复杂.高程方向的负载形变的方向与负载的增减有关(负载增加即向下移动，负载减少即向上移动)；而水平方向的负载形变的方向不仅与负载的增减有关，还与负载所在的方向有关.一定区域内负载的增减是规律的，但负载的分布是随机的，也就导致了在区域负载显著的地方，水平方向的负载形变杂乱无章，尤其是欧洲地区.

与地球物理模型对比发现，GRACE技术能对大尺度的质量迁移进行高质量的监测.由于测量模式的不同，GPS监测的地表形变与GRACE还存在一定的差异.本文对GPS/GRACE的季节性地表差异进行了详细的比较，并初步探讨了差异的原因.结果表明，GPS^ITRF08和GRC^CSR地表形变的周年变化在高程方向上具有较好的一致性：有65%的测站U方向的周年相位差小于30天.GPS的周年振幅总体略大于GRACE；但水平方向的差异明显：近一半测站EN方向的周年相位差大于60天，GPS的EN方向周年振幅为GRACE的2~3倍.

在此基础上研究了三个可能影响GPS/GRACE地表形变差异的因素：不同GPS站地表形变时间序列间的不确定性，热弹性形变和区域形变.结合Fang等的理论指出热弹性形变是导致GPS和GRACE的N方向差异的主要原因.区域负载对GPS/GRACE水平方向差异的影响也是不可忽略的，特别是欧洲区域.另外，GPS形变本身的不确定度也是导致GPS/GRACE差异的一个重要原因，特别是水平方向，但GPS/GRACE之间的季节性地表形变的差异并不主要由GPS技术误差引起.

本文仅从定性的角度分析了热弹性形变对GPS/GRACE形变差异影响；后续将利用Fang等建立的热弹性形变模型计算GPS测站的热弹性形变，并分析消除热弹性形变后的GPS形变序列与GRACE形变序列的差异.另外本文尚未考虑某些GPS和GRACE技术相关误差的影响，如周期为1.04cpy(cycles per year)的虚假信号和多路径误差等等.

利用GPS数据研究地表质量迁移时，有必要消除热弹性形变的影响.另外，与GRACE相比，GPS还能反映更小尺度上的质量负载.一旦GPS测站达到一定的密度，基于GPS站监测地表质量变化，进而监测全球和区域气候变化和水循环具有十分重要的意义.