运用GPS对地壳进行监测，主要是在一系列固连在稳定基岩上的地面观测墩上进行观测，以此来监测地壳的运动，从而反映地球内部的构造运动。然而，季节性的水圈运动会引起岩石圈的周期性位移，GPS的观测数据中必然会包含非地球内部构造运动引起的地表位移。部分学者开始关注这一现象，并对区域或全球由陆地水(主要为土壤水和冰雪等)、大气、海洋和海洋非潮汐等负荷变化引起的GPS连续站季节性变化(尤其是垂直方向)进行了研究。如何改正GPS观测中非构造运动引起的地表位移，是目前的研究热点。

重力场恢复与气候实验卫星(Gravity Recovery and Climate Experiment, GRACE)以前所未有的精度和分辨率获取全球时变重力场信息，其提供的时变重力场模型可反映非大气、非海洋的质量随时间变化。GRACE重力卫星对陆地水负载敏感，可定量地分析陆地水变化，在季节性尺度上，对陆地区域而言，主要表现为陆地水储量的变化。因此，用GRACE改正GPS观测站的季节性变化是一条比较理想的途径。本文的思路就是利用GRACE观测到的季节性效应去改正GPS台站时间序列中的季节项并评估其有效性。

青藏高原南部地区受到夏季季风的影响，降水呈现显著的季节性效应，由此导致该区域的GPS台站存在季节性波动^[1-3]，该季节性波动是应用GPS进行青藏高原地区构造运动研究的一项重要误差。本文使用美国德克萨斯大学空间研究中心(Center for Space Research, CSR)提供的GRACE月重力数据，GPS数据采用中国大陆构造环境监测网络(Crustal Movement Observation Network of China，CMONOC)的GPS基准站和境外(主要分布在尼泊尔)的GPS连续站。LHAZ站为陆态网络1期站点，从1998年开始观测；国内其余站点均为陆态网络2期站点，从2010年开始记录数据；境外站点由加州理工大学布设，观测时间相对较长，跨度在1997年至2013年。图 1和表 1分别展示了站点分布和观测时间的详细信息。

图  1  青藏高原南部的GPS站点分布图

Figure  1.  Distribution of GPS Sites in the Southern Qinghai-Tibet Plateau

下载: 全尺寸图片 幻灯片

表  1  GPS站点信息

Table  1.  GPS Site Information

点名	纬度/(°)	经度/(°)	观测时间(年积日)	RWRMS	周年振幅/mm		半周年振幅/mm		相关系数
					GPS	GRACE	GPS	GRACE
BMCL	28.66	81.71	07077-13160	0.62	5.78	4.88	3.03	1.91	0.57
BYNA	29.47	81.20	09106-12179	0.67	5.87	4.18	3.07	1.76	0.85
CHLM	28.21	85.31	04091-13240	0.80	6.89	5.83	3.18	1.77	0.50
DAMA	27.61	85.11	04001-10099	0.73	6.23	6.21	3.90	1.89	0.86
DLPA	28.98	82.82	07131-13091	0.85	6.42	4.86	2.26	1.79	0.67
DNGD	28.75	80.58	08128-13241	0.61	3.58	4.55	2.93	1.92	0.20
DRCL	29.73	80.50	08073-13241	0.80	5.01	3.82	4.17	1.73	0.58
GNTW	29.18	80.63	08120-12182	0.44	9.68	4.26	4.40	1.84	0.72
GRHI	27.95	82.49	07128-13241	0.58	6.02	5.53	2.60	2.00	0.72
JMLA	29.28	82.19	07135-13241	0.70	6.95	4.53	4.76	1.76	0.83
JMSM	28.81	83.74	04124-12142	0.77	9.33	5.15	0.93	1.77	0.83
KKN4	27.80	85.28	04080-13241	0.50	10.48	6.11	3.70	1.84	0.75
KLDN	27.77	83.60	04102-12343	0.39	7.36	5.85	2.90	1.97	0.69
LHAZ	29.66	91.10	09001-13277	0.77	8.66	5.39	2.16	1.04	0.93
NPGJ	28.12	81.60	07139-13179	0.61	8.26	5.24	2.44	2.01	0.51
ODRE	26.87	87.39	04314-13153	0.70	10.01	7.11	3.13	1.77	0.62
RBIT	26.85	89.39	04001-06001	0.79	9.96	7.45	1.70	1.56	0.41
RMJT	27.31	86.55	09310-13215	0.65	11.17	6.67	5.79	1.80	0.92
RMTE	26.99	86.60	08266-13241	0.68	8.87	6.88	3.13	1.84	0.76
SIM4	27.17	84.99	04087-10098	0.63	10.19	6.48	4.03	1.96	0.82
SMKT	29.97	81.81	08139-13241	0.84	4.93	3.95	3.16	1.63	0.82
SRNK	28.26	83.94	05106-12157	0.72	4.56	5.57	3.69	1.86	0.82
SYBC	27.81	86.71	08277-09328	0.60	6.75	6.35	4.79	1.71	0.58
TIMP	27.47	89.63	04001-08095	0.77	10.28	7.05	2.75	1.47	0.79
TPLJ	27.35	87.71	04069-13091	0.96	6.56	6.85	1.30	1.68	0.46
XZAR	29.27	87.18	11117-13131	0.60	7.86	5.32	3.95	1.43	0.74
XZBG	30.84	81.43	11091-13131	0.49	7.08	3.26	3.45	1.45	0.68
XZDX	30.48	91.10	10263-13131	0.68	6.96	4.72	3.84	0.94	0.88
XZGE	32.52	80.11	10343-13131	0.25	2.98	1.97	4.63	1.12	0.44
XZNQ	31.49	92.11	10255-13131	0.50	6.42	4.03	4.93	0.75	0.82
XZRK	29.25	88.87	11160-13131	0.53	7.23	5.55	4.95	1.29	0.71
XZRT	33.39	79.72	12299-13131	*0.07	*30.83	1.56	*17.14	0.93	0.38
注：*表示该数据存在异常的情况，可能是数据观测周期太短，受到异常值的影响较大

下载: 导出CSV 

| 显示表格

1 大地测量数据及其处理

1.1   GPS数据处理

GPS数据处理采用GIPSY/OASIS软件(5.0版本)精密单点定位模式，利用测站每日记录的相位和伪距数据计算当日测站坐标三分量及其方差-协方差，每日形成一个测站坐标解算结果文件(单日解)。计算时固定模糊度，使用美国的喷气推进实验室(Jet Propulsion Laboratory, JPL)的精密轨道/钟差，绝对天线相位中心模型，对流层模型方面采用全球映射函数(global mapping function, GMF)及全球气压和温度(global pressure and temperature, GPT)模型提供的对流层干分量作为先验值。在海潮模型改正方面，为避免系统误差，并与JPL的精密轨道和钟差产品保持良好的自洽性^[4]，采用TPXO7.0海潮模型，并通过格林函数将其改正到地心(center Earth system, CM)参考框架。

采用上述解算策略和模型，获取测站的无基准单日解，通过相似变化得到ITRF2008坐标参考框架下各站点北(N)、东(E)、垂直(U)各坐标分量的时间序列。

1.2   GRACE数据处理

目前提供GRACE产品的机构有法国的空间大地测量学研究组(Space Geodesy Research Group, GRGS)、美国的CSR、德国地学研究中心(GeoForschungsZentrum，GFZ)和美国的JPL，其中除了GRGS的RL02产品为10日解，其他产品均为月解。

本文使用CSR提供的RL05(Level-2 Release-05)月重力场产品GSM，时间跨度为2002-2014年。数据处理中采用文献[5]的一阶项结果代替；采用高斯平滑滤波方法，选取平滑半径为400 km；GRACE轨道存在系统误差，使用高精度人卫激光测距仪(satellite laser ranging, SLR)计算的C20替换^[6]；利用文献[7]方法做去相关性处理。

1.3   GRACE数据反演地表垂直变形

根据文献[8]的理论，GRACE反演地表垂直位移Δh的数学表达为：

式中，r为地球半径；θ、φ为余纬和经度；P_lm为正则化的勒让德函数；h′ _l和k′ _l为负荷勒夫数，本文中采用文献[9]计算的数值；(ΔC_lm, ΔS_lm)是GRACE的球谐系数产品。另外，值得注意的是GRACE重力场的变化通过球谐系数的变化来体现，球谐系数的变化一般有3种方法进行计算：(1)减去研究时间区内的平均系数值；(2)减掉某一静态重力场模型的系数值；(3)减掉某一月的系数值。本文采取第一种计算方法，减去2002-2014年间的平均系数值。

2 GPS和GRACE反演垂直形变结果分析

2.1   周期项与振幅

GPS时间序列资料的分析，通常考虑长周期运动、同震变形以及震后变形等因素。GPS单站、单分量坐标序列y(t_i)一般用下列周期性模型来表示^[10]：

式中，n_g表示有n_g个地震(阶跌)发生；t_i是以年为单位；系数a、b分别代表地壳位置和线性变换率；系数c、d、e和f描述了周年和半周年运动振幅；系数g_j代表地震造成的同震偏移；H(t_i-T_gj)为阶跃函数；v_i为残差值，代表观测值与预测值间的差异。

时间序列中的周年项和半周年项的振幅(A_ann和A_semi-ann)、相位(Φ_ann和Φ_semi-ann)计算公式为：

GPS结果中包括由地理物理现象造成的长期线性项(例如板块运动和冰后回弹等)，而本文重点关注GPS站坐标时间序列的季节性波动，因此，首先扣除GPS时间序列中的趋势项部分，暂不对长期线性变化进行研究。根据式(2)~(4)，可获得GPS站坐标时间序列线性及周期项拟合结果，周年项及半年项参见表 1。从表 1中可以看出，GPS站的垂向位移时间序列周年项振幅最大可达11.17 mm，半周年项振幅最大为5.79 mm(均为位于尼泊尔境内的RMJT站)。

图 2所示为CHLM和KKN4 GPS站坐标时间序列的线性拟合与周期性拟合曲线，从图 2中可看出二者时间序列存在明显的周期性，GPS周年振幅分别高达6.89 mm和10.48 mm，这对地壳形变研究而言，都是不容忽视的误差。无论从相位还是振幅分析，CHLM站上GRACE计算结果和GPS观测结果都吻合很好。

图  2  GPS观测值及其拟合的线性部分和非线性部分与GRACE模型比较

Figure  2.  Seasonal Variations in Two Sites Compareing GPS Observations, GPS Nonlinear Fit, GPS Linear Fit with CSR RL05

下载: 全尺寸图片 幻灯片

2.2   相关性分析

两种不同空间技术获取的时间序列，其周期性若存在相关，则存在两种可能：

1) 两个周期项实为同一周期，即由同一个物理因素所引起，由于所采用分析方法的不同或者分辨率和定位精度所限，而显现出两个周期。

2) 引起两个周期项的物理因素之间具有相关性。

因此，可采用GPS与GRACE之间的相关性来判断两者是否为同一周期项，或者为由同一物理因素产生的周期振动。若相关性不显著，则需进一步探讨。

比较GPS和GRACE时间序列的相关系数(见表 1)，发现GPS和GRACE获取的垂向位移时间序列存在一定相关性，多数站点的相关系数在0.7左右，特别是沿喜马拉雅造山带区域内的站点相关性更高。结合GPS观测与GRACE反演模型时间序列相关性的空间展布(见图 3)，可以看出二者存在一定的空间相关性。

图  3  GPS观测与GRACE反演模型数据在垂直向相关性图

Figure  3.  Correlation Map Between GPS Observation and GRACE Inversion Model in Vertical Displacement

下载: 全尺寸图片 幻灯片

2.3   一致性分析

为了更进一步论证GPS与GRACE垂直分量的周期(季节项为主)一致性，定量评估GRACE改正GPS数据的有效性，可估算改正前后GPS测站时间序列加权均方根误差(weighted root mean square, WRMS)，公式为：

式中，n为单日解的个数；c_i(i=1, 2, 3…n)是测站坐标序列的垂直分量(去掉趋势项)；σ_i是GPS垂直分量标准误差；Δh_i^G是测站垂直分量的GRACE改正；σ_Gi是GRACE反演垂直分量标准误差；W_GPS和W_GPS-GRACE分别表示GPS测站时间序列WRMS和经过GRACE改正之后的WRMS。当W_GPS-GRACE小于W_GPS时说明附加GRACE季节项改正是行之有效的。同样也比照上面公式估算GPS坐标时间序列扣除自身拟合后的W_GPS-GPSfit。采用定量评估公式来计算GPS基准站的WRMS减小百分比R_WRMS，公式为^[11]：

R_WRMS可反映GPS与GRACE时间序列年周期项在振幅相位上的一致性，当R_WRMS值为1.0时，表示GPS与GRACE分别拟合得到的周期项垂直变形完全一致。理论上，GPS测站的WRMS减小比率不超过1。通过上述公式，分别计算了青藏高原南部的32个GPS测站的WRMS和R_WRMS值(见表 1)。

进一步统计发现，经GRACE改正之后，GPS的垂向位移时间序列WRMS减小，R_WRMS的平均值达到了0.64；各测站的R_WRMS值如图 1所示。

图 4显示了GRACE改正GPS垂向位移时间序列的效果。小圆点越靠近黄色方柱的底部，表示两者一致性越好，即R_WRMS数值也就也接近1。从图中看到W_GPS-GRACE比W_GPS小，这说明利用GRACE改正GPS观测值是有效的。

图  4  GRACE改正GPS之后WRMS的变化图

Figure  4.  WRMS Reduction when Removing GPS Seasonal Variation Using GRACE-Modeled Seasonal Displacements

下载: 全尺寸图片 幻灯片

3 讨论

3.1   季节项波动的时间与区域特征分析

从图 2可以明显看出GHLM与KKN4的时间序列存在明显的周期性，并且冬季上升，夏季降落，周期性曲线下降峰值出现在每年的7、8月，这与夏季季风带来强降水而引起地表质量负载增加的时间区间一致^[12]。由此可以推断，青藏高原南部地区所呈现的季节项波动很大程度受到降水而导致的水文负载影响，夏季时地表负荷变大而冬季时减弱。

从图 3中可以看出GPS与GRACE获得的时间序列存在相关性，但相关性在整个青藏高原南部地区并非完全一致，而是表现出一定的局部特征。图 3中的西北处(XZRT、XZGE)阿里地区比图 3中东北处(LHAZ、XZDX等)拉萨地区的相关性低，尼泊尔中部地区站点的相关性比东西两侧地区的站点高。XZRT和XZGE两个站点的相关性不显著，可能的原因是较之拉萨地区，此地的季节性的水文作用并不强烈，从而导致该地区站点的GPS与GRACE垂向负荷的相关性降低。

对比分析图 1和图 3可发现二者存在一定的空间一致性，按照GPS和GRACE的相关性和R_WRMS数值，大致可以将青藏高原南部地区划分成拉萨、尼泊尔中部和西北阿里3个区域。尼泊尔中部位于喜马拉雅山南麓，夏季季风带来丰富的降水；拉萨区域由于受到喜马拉雅山的阻挡，夏季季风影响有所减弱；西北阿里地区更远离喜马拉雅山，更接近青藏高原内陆，受夏季季风的影响更小。

3.2   季节项波动的成因分析

上述分析认为，青藏高原南部地区的季节项波动主要是由于夏季季风带来降雨并引起水文负载增加所致，但是GRACE反演结果与GPS观测的季节项并非完全一致，总体而言GRACE反演获得的垂向运动周年项和半周年项振幅均小于GPS的结果。从两种空间技术可能包含的误差分析，采用GRACE原始数据解算GSM重力场模型时，已使用相应背景重力场模型将非潮汐大气和高频海洋信号以及各种潮汐(包括海潮、固体潮、固体极潮)的影响扣除；GRACE计算的时间序列中主要包括由于水质量(地下水、土壤湿度、雪和冰，以及海洋质量)的重新分布造成的重力变化。而GPS观测的时间序列中包括由于地理物理现象造成的长期线性项(例如板块运动和冰后回弹等)和水文负荷、大气和海洋非潮汐因素影响。除此之外，GPS观测结果中还包括一些未模型化的误差(例如海潮改正残差、参考框架的周期变化、标墩的热胀冷缩、GPS轨道的交点年周期误差等)，该类误差也呈现出周期性特征。因此，即使从GPS测站时间序列中剔出了由于水文负载导致的季节项，还是会有其他季节项的残余，需要进一步研究。

4 结语

在GPS测站坐标时间序列中包含短周期(年-半年周期等)水文负荷造成的季节项，对研究构造形变而言是一种噪声，需要剔除。本文利用青藏高原南部GPS站点多年观测资料，结合GRACE数据反演得到该地区的季节项波动，对比GRACE对陆地水敏感的特性，并结合季节项出现的时间区间分析认为，青藏高原南部地区的季节项主要由夏季季风降水导致的水文负载变化引起。该地区GPS观测的垂向时间序列与GRACE反演的地表垂向位移存在相关性，多数站点相关性在0.7左右，使用GRACE观测获取的季节项对GPS测站时间序列进行改正，WRMS有明显减小。R_WRMS存在明显差异，最大和最小值分别是0.96(TPLJ)和0.24(XZGE)，结合相关性与R_WRMS，进一步研究发现不同区域季节项及其改正效果呈现较显著的局部特征，这也与夏季季风的影响范围和大小呈现一定的关联。除此之外，GPS测站时间序列中还存在一些未模型化的误差，也可能引起周期性变化，具体原因需要进一步分析。青藏高原南部地区内部不同地区构造运动强度、地表层的水负荷周期变化存在明显差异，其动力细节需要更进一步深入研究。