地球物理学报  2014, Vol. 57 Issue (6): 1792-1804   PDF    
利用GPS与GRACE监测陆地水负荷导致的季节性水平形变:以喜马拉雅山地区为例
王林松1,2, 陈超1, 邹蓉1, 杜劲松1, 陈晓东2    
1. 中国地质大学(武汉)地球物理与空间信息学院, 武汉 430074;
2. 中国科学院测量与地球物理研究所 大地测量与地球动力学国家重点实验室, 武汉 430077
摘要:地表陆地水负荷变化是引起重力场和地壳形变呈现季节性特征的主要因素,并且能够利用地表及空间大地测量技术对其进行有效的监测.本文通过对质量负荷形变效应的理论模拟,描述了水平分量的形变指向以及垂直与水平分量的幅值比可以提高对负荷区域的辨别程度,并且联合GPS坐标时间序列及GRACE模型对喜马拉雅山地区的季节性负荷形变进行了详细对比分析,研究结果显示两者垂直分量的季节性变化具有较好的一致性,且GPS周年项幅值要大于GRACE.而由GRACE解算得到的水平分量结果表明该地区季节性形变主要受东南亚及印度东北部地区的陆地水负荷控制,位于喜马拉雅山地区多数GPS台站的垂直分量及北向分量的初相位与GRACE模型解算结果相近,而部分GPS台站的东向分量与GRACE模型存在明显不同,由此导致GPS与GRACE监测到的形变指向存在差异.通过对GRACE估算精度以及GPS垂直与水平分量幅值比的深入分析,发现GPS对局部周边地区的河流、谷地及农田灌溉等负荷变化造成的形变效应较为敏感,而GRACE由于截断阶次及平滑滤波等影响因素,不仅造成在水平分量上的分辨率远低于垂直分量,而且整体估算精度要低于GPS观测得到的形变信息.
关键词GPS     GRACE     季节性水平形变     地表陆地水负荷     喜马拉雅山    
Using GPS and GRACE to detect seasonal horizontal deformation caused by loading of terrestrial water:A case study in the Himalayas
WANG Lin-Song1,2, CHEN Chao1, ZOU Rong1, DU Jin-Song1, CHEN Xiao-Dong2    
1. Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China;
2. State Key Laboratory of Geodesy and Earth's Dynamics, Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan 430077, China
Abstract: Seasonal variations in the gravitational field and crustal deformation are mainly caused by surface terrestrial water loading, which can be detected using ground- or space-based geodetic techniques. In this paper,in terms of theoretical modeling for the deformation effect of mass loading,we describe the direction of horizontal deformations and the ratio of the vertical displacement to the horizontal displacement which can improve discrimination of the surface loading region. And then,with the help of the GPS time series and the GRACE data we analyze the seasonal deformations of surface loads in the Himalayas. The results show the significant seasonal variations and a good correlation between GPS and GRACE in the vertical component,and the annual amplitude of GPS is greater than GRACE. However,the horizontal component variations derived from GRACE monthly gravity field models are dominated by the surface terrestrial water loads of Southeast Asia and Northeast India in the south region of the Himalayas. In addition, the phase of the vertical component and the northern component from the most GPS time series which is similar to estimation results derived from the GRACE in the Himalayas,but the eastern components from a few GPS stations and GRACE show a significant difference,which results in the difference between GPS and GRACE in direction of deformation. Through the analysis of the GRACE estimation accuracy and the ratio of GPS,further comparison indicates that the GPS is more sensitive to the deformation effect induced by the local load changes (i.e.,rivers,valleys and farmland irrigation,etc.). Furthermore,because of the influence factors such as the truncation error,and as well as the smoothing filtering effect,GRACE data shows that not only spatial resolution of the horizontal component is much lower than vertical component,but the accuracy of the GRACE deformation prediction is also less than GPS.
Key words: GPS     GRACE     Seasonal horizontal deformation     Surface terrestrial water loading     Himalayas    

 1 引言

利用GPS监测地壳运动,其坐标时间序列不仅包含构造运动引起的线性变化,并且也呈现出非常显著的非线性特征(通常为周年尺度的季节性变化),姜卫平等(2013)认为地球物理效应及GPS技术相关的系统误差是造成非线性变化的主要因素.这种地球物理效应一般是由环境变化(如大气、海洋及陆地水)主导的地表质量负荷响应,所导致的地壳形变为国内外学者利用GPS研究环境效应提供了广阔空间,取得的研究成果主要在大气负荷(van Dam et al., 1994Boehm et al., 2007)、潮汐海洋及非潮汐海洋负荷(van Dam et al., 1997Kusche and Schrama, 2005)及陆地水储量(Sauber et al., 2000Blewitt et al., 2001Bevis et al., 2005van Dam et al., 2007)等几个方面.除地表质量负荷影响之外,还有一些可能引起季节性形变的系统因素,如校正模型的适用性、基岩的热胀冷缩、相位中心模型误差、软件和解算策略的影响等(Dong et al., 2002张飞鹏等,2002Yan et al., 2009),而对此类因素的影响程度等问题的研究目前正处于起步阶段(姜卫平等,2013).

近年来,卫星重力观测技术的快速发展,特别是 GRACE(Gravity Recovery and Climate Experiment)卫星任务的发射与实施取得巨大成功,其主要任务是探测地表质量变迁导致的重力场变化(Tapley et al., 2004),相关的研究主要涵盖了地表水、大型流域与河流、冰雪覆盖等引起的时变重力响应(Chen et al., 2006Luthcke et al., 2008杨元德等,2009廖海华等,2010Bruinsma et al., 2010). 由于GRACE模型系数能够直接转换成因表层质量负荷引起的地壳形变的三个分量(U分量、N分量及E分量),因此联合GRACE及GPS两个观测手段,能够从质量负荷的研究角度定量地对比分析两者所获得的季节性(van Dam et al., 2007)或长期性(Khan et al., 2010)负荷形变特征.

喜马拉雅山地区位于青藏高原南巅边缘,作为地球上海拔最高的山脉,不仅是中国与印度、尼泊尔、不丹、巴基斯坦等国的天然国界,而且该地区又因印度板块与欧亚板块的碰撞,其形成的起因与演变一直是国内外研究的热点问题,例如青藏高原的隆升与增厚(Sun et al., 2009)及板块推挤速度(Ader et al., 2012)等.同时,该地区巨大的冰川积雪覆盖是亚洲地区大型河流与人类生活用水的主要来源,并且积雪量呈现出明显的季节性变化,其周年变化的积雪与融化是喜马拉雅山南端恒河流域水量季节性变化的主要贡献者(Fu and Freymueller, 2012),也是造成该地区地震活动呈季节性特征的主要影响因素(Bettinelli et al., 2008).

以往利用GPS或GRACE监测因地表质量负荷引起的季节性形变的研究中,重点主要集中在地壳运动的垂直分量(van Dam et al., 2007Tesmer et al., 2011Fu and Freymueller, 2012),由于缺少水平分量的综合解释,因此不能详细提供单一或主要负荷区域的位置信息.而Wahr等(2013)研究发现水平位移对约束负荷区域的位置以及验证垂直位移是非常有价值的. Fu等(2013)也证实了利用GPS和GRACE的季节性水平位移能够确定因雨季造成的地球弹性响应的主要负荷区域. 本文基于Wahr等(2013)的研究思路,主要目的在于如何利用地壳形变的水平分量研究季节性变化特征,通过简单理论模型,分析主要地表质量负荷区域的形变指向,并且结合形变垂直与水平分量周年项的幅值及初相位的关系,区分单一及复杂负荷条件下的形变特点;本文以喜马拉雅山地区为例,探讨利用GPS和GRACE所获得的水平位移分析影响该地区季节性形变的主要陆地水负荷区域,并且对GPS与GRACE研究因陆地水效应引起的季节性变化的适用性与影响因素等方面进行了讨论. 2 基本理论 2.1 地表负荷形变

假设在球对称、自引力、弹性和各向同性的地球表面 θ 位置施加一个角半径为 α 的均匀圆盘质量负荷,则在地表(地球半径为 a)等效质量分布能够表示为(Farrell,1972)

而地表质量负荷引起的重力位扰动为(Farrell,1972Wahr et al., 2013)

其中,G为牛顿万有引力常数. Pn(cosθ)为n阶的勒让德多项式;结合负荷勒夫数,由负荷引起的地表垂直位移(向上为正)及水平位移(远离负荷为正)的格林函数可以分别表示为(Farrell,1972Wang,2000Wang et al., 2002)

式中,me为地球质量,hn与ln为n阶的负荷勒夫数,本文采用Jentzsch(1997)文献中利用PREM模型(Dziewonski and Anderson, 1981)所给出的负荷勒夫数.

2.2 形变特征描述

由上述公式(1)—(5)得知,地表任意一点的质量负荷变化,能够通过GPS的观测给出该负荷特征信息,例如垂直位移的正负变化可以确定负荷的移除或增加,但无法得知具体的负荷区域的位置信息.如果结合水平位移信息,我们不仅能够了解负荷质量的增加或减少,而且可以分辨出观测位置相对负荷位置的运动特征(背离或靠近).因此,利用垂直与水平形变信息可确定简单负荷模型的强度与空间位置特征.

根据负荷形变理论,假设一个半径为500 km,等效水柱高为50 cm的圆盘水文模型加载到地球表面,由此导致的周边离负荷中心点1250 km以内的垂直形变在0 ~-1.63 cm(向上为正)范围内变化(如图 1).而结合水平形变的正负变化,我们能够定 义观测位置对负荷区域的形变指向,例如水平形变的N分量(向北为正)、E分量(向东为正)分别与U分量(向上为正)变化方向相同时定义为正,相反时为负(图 1),即当负荷加载时,负荷区域的西北方向水平分量的变化方向为(+,-)、东北方向为(+,+)、东南方向为(-,+)及西南方向为(-,-),具体指向角度则由E分量与N分量幅值大小决定(angle=arctan(E/N)).

图 1 负荷形变特征
图中黑点代表负荷质量为半径500 km且等效水柱高为50 cm的圆盘模型,箭头表示利用水平分量与垂直分量变化方向的关 系定义的形变指向,例如(+,-)表示N分量变化方向与U分 量相同,E分量与U分量相反.
Fig. 1 Deformation characteristics caused by loading
Black dots represent the surface mass model,which has a uniform disc load of radius 500 km and 50 cm equivalent water thickness. The arrows represent the direction of motion using the relationship between horizontal component and vertical component, and sign “+” shows that the direction of N component is same as the U component, and “-” indicates E and U component are in the opposite directions.

对于集中的质量负荷,垂直形变量往往大于水平形变量.在理论模拟中,采用不同的地球模型计算得到的垂直与水平分量的比值关系稍有不同(Pinel et al., 2007Wahr et al., 2013);对于同一地球模型,影响其比值的因素主要与负荷模型性质及观测点与负荷区域的距离有关(图 2).图 2a图 2b显示了不同半径的单一圆盘负荷模型(加载50 cm的等效水柱高)引起的垂直与水平的比值在2~3之间,随着观测距离的增加而趋于稳定(约2.3),其中图 2a中等效水柱高为50 cm不同半径模型加载的理论模拟结果与Wahr等(2013)采用的等效水柱高为100 cm的预测结果非常相近,说明单一负荷模式下的垂直与水平位移之间有一定的规律可循;而在非单一负荷模式下,即外部区域周边也存在负荷质量变化(单一圆盘负荷模型不变)的垂直与水平位移的比值,则依赖于观测点周边负荷变化特征及观测距离,例如负荷外部区域周围同样存在负荷质量加载5 cm的等效水柱高时比值较大(图 2c),相反负荷外部区域周围存在负荷质量移除5 cm的等效水柱高且远离负荷中心一定距离内(100~200 km)时比值减小(图 2d).由此可见,多个负荷源造成的形变叠加后垂直形变与水平形变的相对关系是与负荷大小、作用形式(加载或移除)及距离等因素的综合反应,很难有规律可循,因此图 2c和2d是在特定的负荷叠加条件下的结果,并不具有普遍意义.但对于实测资料而言,根据监测到的垂直与水平位移不同的比值关系,结合理论模拟结果,有助于我们对观测点局部或区域负荷形式进行判别.

图 2 相同地球模型条件下的垂直与水平分量比值的模拟
(a)为单一不同半径(10~100 km),加载50 cm的等效水柱高负荷模型下的比值变化;(b)为单一不同半径(100~500 km),加载50 cm的 等效水柱高负荷模型下的比值变化;(c)与(d)负荷模式与(b)相同,但负荷外部区域存在加载5 cm的等效水柱高的负荷模式(c)与移除 5 cm的等效水柱高时的负荷模式(d).
Fig. 2 Simulated ratios of vertical to horizontal displacements in the same Earth model
(a)Ratios in the model loaded by single 50 cm equivalent water column with different radii from 10 km to 100 km;(b)Same as(a)but with different radii from 100 km to 500 km;(c)Same as(b)but with additional 5 cm equivalent water thickness in outside area;(d)Same as(b)but with removing 5 cm equivalent water thickness.
3 观测资料 3.1 GPS坐标时间序列

本文采用的GPS坐标时间序列主要来自于Caltech研究小组安装在尼泊尔地区观测时间大于2年的连续台站(图 3),经由Caltech′s Tectonics Observatory 网站下载得到,数据经过了GAMIT/GLOBK软件包(Herring et al., 2009)的处理,且基于ITRF2005(Altamimi,2009)框架下的全球IGS观测网络提供的5个台站分析的综合解,其处理过程中扣除了固体潮、海潮和极移等影响,但校正过程未包含大气、非潮汐海洋及陆地水等负荷的影响.对于大气负荷的影响,我们采用由GGFC(Global Geophysical Fluid Center)提供的大气负荷形变效应计算程序(van Dam,2010. http://geophy.uni.lu/ncep-loading.html),其中,大气压力来自NCEP(National Centers for Environmental Prediction)提供的采样间隔为6小时,空间分辨率为2.5°×2.5°的数据.而对于非潮汐海洋负荷引起的非线性变化,由于其对非近海台站引起的形变幅值较小(王敏等,2005姜卫平等,2013),因此该影响包含在本文的GPS坐标时间序列中.

图 3 位于喜马拉雅地区GPS连续台站的位置示意图
研究区域内(黑色虚线方框)18个台站位于尼泊尔境内,1个台站位于不丹境内,2个台站位于拉萨地区,蓝色菱形为本文列举的GPS坐标时间序列实例的台站位置.
Fig. 3 Locations of continuous GPS stations in the Himalayas
In the study area(black dotted box),18 stations lies in Nepal,1 station in Bhutan and 2 stations in Lhasa, and locations of example GPS coordinate time series in this paper are shown by blue diamonds.

由于GPS数据持续时间长短不一,在时间尺度上主要反映的是季节性和年际变化,因此我们可以采用如下形式对GPS坐标时间序列的各个分量进行最小二乘拟合:

其中a为常数项,b为线性项系数,ω为频率(ω=1/T,本文以年为周期时T=1),c为年变化振幅,d为初始相位,t为时间(以年为单位),初始历元t0=2004.0;同时考虑到与GRACE模型数据的对比分析,下文中利用公式(6)对GRACE估算得出的季节性变化进行相同的最小二乘拟合.

表 1为拟合得到的GPS三个分量的周年振幅与初相位,同时也给出了U分量与H分量(H= )振幅比值,根据与理论模拟的比值对比,可将U与H分量周年振幅比值分为三个范围(R<2.0、2.0<R<4.0及R>4.0).表 1中的正负表示为观测位置对负荷区域的形变指向,其正负由GPS坐标时间序列各个分量的初相位之差决定.当初相位之差小于90°时,我们近似定义为变化方向相同(即为正号);当初相位大于90°且小于180°时,近似定义为变化方向相反(即为负号).严格来讲,由U分量与N分量初相位的对比结果表明,位于喜马拉雅山地区的多数GPS连续台站的初相位相差较小(<20°),说明引起喜马拉雅山地区季节性形变的负荷区域主要位于喜马拉雅山脉的南端;但部分观测台站(DNGD、NPGJ及GRHI)的垂直与水平分量的初相位也存在较大差异,表明这些台站GPS坐标时间序列的季节性变化受局部及区域性多种负荷因素的影响.

表 1 各台站GPS与GRACE季节性变化周年振幅、初相位、幅值比值及水平形变方向 Table 1 Amplitudes,initial phases,amplitude ratios and deformation directions derived from GPS and GRACE
3.2 GRACE模型

利用GRACE模型给出的地球重力场球谐系数以及负荷勒夫数(Jentzsch,1997),不仅能够估算出负荷质量的变化(Wahr et al., 1998),而且可以根据弹性负荷形变理论(Farrell,1972)获得负荷引起的形变信息(Kusche and Schrama, 2005van Dam et al., 2007).目前GRACE月平均重力场模型的处理 与发布主要由CSR(Center for Space Research,USA)、 GFZ(GeoForschungsZentrum,Germany)及JPL(Jet Propulsion Laboratory,USA)三个研究中心 提供;此外,GFZ与GRGS(Space Geodesy Research Group,France)分别提供了时间分辨率为7天与10天的重力场模型.研究表明,除个别发布的版本之外(例如JPL RL04),不同的GRACE模型估算结果显示出非常好的一致性(杨元德等,2009Fu and Freymueller, 2012).本文采用CSR提供的GRACE RL05模型,时间跨度从2003年1月到2012年12月,该系数最高阶次为60阶,并且扣除了潮汐影响和非潮汐的大气和海洋影响.由于C20项误差较大,本文将此项替换为利用SLR观测的结果(Cheng and Tapley, 2004),而模型中的一阶项我们则采用Swenson等(2008)提供的一阶项系数.由于模型受卫星的轨道误差和球谐系数的截断误差等影响,此处采用多项式去条带滤波(Swenson and Wahr, 2006)和Fan平滑滤波(Zhang et al., 2009)相结合的方法对模型系数进行滤波,其中Fan滤波的径向与纬向高斯滤波半径分别取400 km.

根据文献(Wahr et al., 1998;Kusche and Schrama, 2005van Dam et al., 2007)给出的公式,推导由GRACE重力场模型解算的地表三个形变分量(U分量、N分量及E分量)分别为

式中R为地球平均半径,θФ 分别为余纬度和经度;Nmax为月平均重力场的最大阶数,取值为60, l,m 为完全规格化的缔合勒让德多项式,对应的阶次分别为l及mClm与Slm为每月重力场模型的球谐系数与整个研究时间范围内的系数平均值之差; ll,hl及kll阶的负荷勒夫数,同样采用Jentzsch(1997)文献中利用PREM模型(Dziewonski and Anderson, 1981)所给出的结果.

根据CSR提供的RL05模型数据,经过解算得到的陆地水储量具有明显的季节性变化特征,图 4上部给出了变化幅值最大与最小(峰-峰值)两个月的解算结果,其中本文研究区域等效水柱高的变化范围能够达到约400 mm,主要陆地水负荷区域集中在位于喜马拉雅山地区南端的恒河流域、印度东北部及东南亚地区,由于在解算过程中受到系数截断误差及滤波等影响,整个区域的空间分辨率较低,在研究区内的等效水柱高变化较为平缓,但变化幅值南高北低的特点依然清晰可见.图 4下部的(a)、(b)及(c)给出了分别位于印度境内(HYDE,IISC)及中国境内(KUNM)三个观测台站垂直分量的时间序列,包括已去除趋势项的GPS时间序列、去除趋势的GPS拟合结果以及GRACE解算的垂直位移.由GPS与GRACE垂直位移幅值对比分析可发 现,两者的幅值在靠近负荷区域的HYDE与KUNM 台站(图 4上部)较为相近,而位于印度南部远离负荷区域的IISC台站监测结果显示GPS幅值要大于GRACE监测到的垂直位移,说明IISC台站位置除受GRACE监测到的负荷区域影响外,其台站也受局部一定范围内负荷的影响.利用前文中所描述的关于水平形变指向的定义,根据拟合得到的GPS和GRACE水平分量与垂直分量的幅值及初相位之差,结果显示出三个台站GPS水平形变都能较好地指向GRACE估算的主要质量负荷区域(图 4上部),而本文研究的喜马拉雅山地区位于GRACE估算陆地水负荷变化区域的西北部.因此,联合GPS和GRACE研究该地区的季节性形变具有一定的可靠性.此外,IISC台站GPS和GRACE水平分量指向的偏差也验证了两者在垂直分量幅值上的不同;并且GPS水平分量幅值要远大于GRACE解算得到的水平形变,但两者垂直分量相差相对较小(图 4下部),说明GRACE在空间上的水平分辨率要低于 垂直分辨率以及地表观测的GPS时间序列,其中关 于GRACE估算精度的影响因素将在下文中详细讨论.

图 4 喜马拉雅山及邻区的等效水柱高分布图及观测台站垂直分量的时间序列
上部为GRACE解算得到的两个月(峰-峰值)的等效水柱高的差值,图中箭头为负荷异常引起的水平形变指向(蓝色为GPS,灰色为GRACE);下部(a)、(b)及(c)三个垂直分量的时间序列分别为经过处理与拟合的GPS观测结果及GRACE解算结果.
Fig. 4 Map showing equivalent water columns and time series of vertical component in the Himalayas and adjacent areas
Upper: Differential equivalent water columns for two months(peak-to-peak)derived from GRACE,the arrows represent directions of deformation caused by abnormal loads(blue for GPS and gray for GRACE). Lower: Time series of vertical displacement derived from GPS and GRACE after processing and fitting for three GPS stations.
4 结果与分析 4.1 GPS与GRACE的对比分析

通过GRACE模型解算可得到的各个台站形变分量的时间序列,利用处理GPS数据时的公式(6)及最小二乘方法对其进行拟合,各台站所在位置拟合得到的GRACE周年项振幅、初相位及幅值比值如表 1所示,各台站水平分量幅值相对垂直分量要小很多,即垂直与水平分量的比值在3~8之间;因负荷区域位于台站南部,因此水平分量中N分量幅值要大于E分量.结合该地区各台站水平与垂直分量初相位之间的关系(图 1中的水平形变指向的理论模拟结合表 1中的(N,E)),显示出由GRACE模型得到各台站水平形变指向大致相同,即负荷区域位于观测点位置东南方向的陆地水季节性负荷区域.

而相比GRACE解算结果较小的差异性,各 GPS台站之间的坐标时间序列无论在周年项以及初相位上都或多或少存在差异,表明GPS对短波长的周边地表负荷有强烈的敏感性.在图 5中,我们将6个台站(图 3b中蓝色菱形观测点)的GPS坐标时间序列与GRACE解算结果进行了比较,两者在垂直分量的时间序列中显示出较好的一致性,说明季节性陆地水负荷是造成形变垂直分量的主要影响因素(Tesmer et al., 2011).对于水平分量,GPS与GRACE的结果虽然显示出两者能够识别大尺度负荷引起的水平形变,但也因台站位置的不同而存在差异,而由表 1拟合结果得知各台站GRACE水平分量有较好的一致性,因此GPS与GRACE不同的主要原因在于各台站GPS观测结果的差异,说明位于研究区域内因局部负荷效应引起的地壳形变也是影响GPS结果的因素之一;相比之下,空间分辨率较低的GRACE模型主要监测范围为长波长的大尺度区域,不能较好的解释各台站水平分量的负荷成因,而GPS观测得到的水平形变受观测点周边一定区域内负荷的影响较大,例如大多数台站GPS水平分量幅值要高于GRACE的估算结果,并且多种不同负荷区域导致的地壳形变也是造成GPS与GRACE水平形变初相位存在不同的原因之一.

图 5 GPS与GRACE三个分量的时间序列对比
GPS时间序列已经过平滑及去除趋势项等处理,拟合结果基于处理后的GPS结果.
Fig. 5 Comparison between GPS time series and GRACE-derived in the three components
GPS time series has been processed by smoothing and removing trends. Fitting results are based on GPS data after processing.
4.2 GRACE估算精度

针对GRACE估算精度,我们可以由GRACE模型的解算理论出发,对于公式(7)—(9)中 Clm与Slm,与地表负荷质量的贡献 σ(θ,Ф) 的关系可表示为

因此,将式(10)代入式(7)可得:

其中,垂直位移的格林函数为

式(12)中的 α 为观测点(θ,Ф)与负荷点(θ′,Ф′)之间的夹角.

由于GRACE模型的最大截取阶数为60阶,因此,对于公式(7)的计算可等价为利用负荷格林函数估算到 lmax=60 的位移变化,即

因此,对于地表质量负荷引起的形变效应,GPS与 GRACE观测结果的差异可视为不同截取阶次(lmax=∞与lmax=60)的格林函数计算结果的不同.

因此,两种方法监测喜马拉雅山地区季节性变 化得到的周年项比值,从另一种角度反映了GRACE 估计季节性形变效应的准确度(如图 6),其中GPS与GRACE的U、N及E分量周年项的平均比值分别为:1.31±0.24 mm、2.59±1.80 mm及3.6±2.46 mm,两者比值大于1表明GRACE低估了实际变化的位移量;相对于垂直分量,GRACE在水平分量上的估算准确度要低,其中N分量比值低于E分量表明该地区的水平位移主要受南北方向上的质量负荷控制.另一方面,对于观测点周边非单一负荷区域,多种负荷(初相位相近)通常会增加垂直分量的幅值,但在水平方向上则有可能造成幅值上的抵消,而这种变化也是GRACE在精度与分辨率方面所不能达到的.

图 6 GPS与GRACE周年项幅值比 Fig. 6 GPS-to-GRACE ratios of seasonal amplitudes
4.3 负荷区域的形变指向

在联合GPS和GRACE对陆地水负荷产生的形变效应进行研究时,许多学者(van Dam et al., 2007Tregoning et al., 2009Tesmer et al., 2011; Fu and Freymueller, 2012)通常利用GPS减去GRACE前后时间序列的WRMS(weighted root-mean-squares)来评价两者之间的一致性,大多数台站在扣除GRACE时间序列后,GPS垂直分量的WRMS都有明显的减小.一些季节性水平形变的研究(Tregoning et al., 2009)表明WRMS也同样会减小,尤其研究地区靠近较大的陆地水储量变化(例如亚马逊流域)的区域.然而,造成一些观测实例中WRMS并没有减小甚至增加的原因主要由GPS与GRACE两者存在较大的相位差引起.同样,在本文对负荷区域的形变指向中,水平与垂直分量之间的相位差是影响该指向不同的关键因素,由GRACE模型解算得到的同一区域水平形变指向基本相同(如图 7),说明东南亚及印度东北部地区大尺度范围内的陆地水负荷是控制GRACE估算水平形变的主要因素.但GPS的形变指向各异,尤其是位于尼泊尔西南部地区的DNGD、NPGJ及GRHI三个台站的差异最大,DNGD台站GPS与GRACE的形变指向相反主要由于该站点周边为大面积的农田,GPS可能受该地区利用地下水灌溉的影响(Shah et al., 2006),而NPGJ与GRHI位于尼泊尔南部的加德满都谷地的Dang地区,该两个台站的水平形变可能由谷地局部地区的负荷效应主导(Fu et al., 2013).因此,这些地区GPS减去GRACE前后时间序列的WRMS不会有明显的改善,观测点周边局部陆地水变化可能会引起GPS明显的位移变化,但由于GRACE估算精度,这种局部的变化将不会被GRACE有效地分辨.

图 7 GPS坐标时间序列水平分量的形变指向及垂直与水平分量周年项幅值比 Fig. 7 Motion directions and vertical-to-horizontal ratios derived from GPS coordinate time series

局部质量负荷变化除能够导致GPS观测指向不同以外,也是影响GPS坐标时间序列垂直与水平分量周年项振幅比的主要因素,而GRACE解算得到的负荷区域(图 4)可近似为由图 2理论模拟中的单一负荷区域,根据观测点局部负荷变化,图 7中垂直与水平的比值大于4的GPS观测台站周边可能存在幅度不一的陆地水进程,且初相位与主要负荷地区相近,结合前文中理论模拟结果(图 2c),说明这些台站在受东南亚及印度东北部地区陆地水负荷影响的同时,其局部也存在与东南亚及印度东北部地区周期相近的陆地水变化,因此引起垂直与水平分量存在较大差异;比值小于4的台站(表 1)除个别台站存在异常比值小之外(DNGD与DRCL),多数台站的比值在2~3之间,部分台站(NPGJ、GRHI及GUMB)的GPS初相位与主要负荷地区相差较大,影响这些台站的主要陆地水负荷可能位于台站周边短波长的一定范围,例如台站几百公里内的恒河流域及各支流的影响,由此说明这些台站局部存在多种的陆地水负荷进程,增加了水平分量的复杂程度,由于本文无地表实测陆地水数据,因此无法从实际负荷的角度解释GPS形成的形变指向差异,这也是今后所需的进一步验证工作. 5 结论

在本文联合GPS和GRACE对喜马拉雅山地区地表质量负荷引起的地壳形变的研究中,我们首先利用单一且集中的负荷模型定性描述了水平分量与垂直分量的关系,包括根据形变方向确定负荷区域的位置,垂直与水平分量的幅值比区分观测点局部与主要负荷区域的差异,这些简单概念的给出能够丰富目前单一垂直分量的解释手段,特别是针对简单且集中(大型湖泊、水库及冰川等)的负荷形式,水平形变结合垂直形变,更能全面分析因负荷形式改变而引起的地壳中长期形变.

通过GPS和GRACE的对比分析,结果显示出在垂直分量时间序列中两者有较高的一致性,而水平分量因观测位置的不同而各有差异.位于喜马拉雅山脉上的GPS垂直分量与北向分量初相位相近,表明该地区水平分量主要受南端的东南亚及印度东北部长波长的负荷区域控制,此结论也与GRACE解算结果相符合,但也存在一些GPS台站与GRACE模型在形变指向以及垂直与水平分量幅值比有较明显的不同.本文结合对影响GRACE分辨率的估算精度的深入分析,发现GRACE由于截断阶次的影响,造成对水平形变监测的分辨率远低于GPS的观测结果,较难识别出局部负荷变化造成的形变效应;而GPS则对此项影响较为敏感,除了受长波长负荷区域(东南亚及印度东北部)的影响外,局部周边地区的河流、谷地及农田灌溉等因素也会导致喜马拉雅山地区的弹性形变,因缺少地表实际质量负荷信息,因此本文未对局部负荷的影响程度进行深入分析.

此外,针对引起季节性形变影响因素的讨论,本文仅从地表质量负荷的角度出发,重点探讨陆地水储量变化对地壳形变的响应,然而GPS观测台站非线性变化由多种因素控制,大气模型的适用性、非潮汐海洋负载等负荷因素都会导致不同程度位移变化.除此之外,与GPS及GRACE观测技术本身相关的潜在因素也不能忽视,例如电离层与对流层延迟,观测点基岩热膨胀效应、映射函数、天线相位中心改正模型以及GRACE轨道与截断误差等. 致谢 感谢审稿人在本文修改过程中给予的宝贵建议.感谢美国科罗拉多大学的John Wahr教授在本文撰写过程中给予负荷形变及GRACE估算精度方面的建议与帮助.文中所用GPS台站数据、GRACE RL05模型数据、计算大气负荷所用数据和程序分别由Caltech Tectonics Observatory、CSR研究中心、Tonie van Dam与NCEP提供,文中部分图件由GMT软件绘制,在此一并表示感谢.

参考文献
[1] Ader T, Avouac J P, Zeng J L, et al. 2012. Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: Implications for seismic hazard. J. Geophys. Res., 117(B4): B04403, doi: 10.1029/2011JB009071.
[2] Altamimi Z. 2009. The International Terrestrial Reference Frame (ITRF2005).//Geodetic Reference Frames, International Association of Geodesy Symposia. Berlin Heidelberg: Springer-Verlag, 134: 81-82.
[3] Bettinelli P, Avouac J P, Flouzat M, et al. 2008. Seasonal variations of seismicity and geodetic strain in the Himalaya induced by surface hydrology. Earth and Planet. Sci. Lett., 266(3-4): 332-344.
[4] Bevis M, Alsdorf D, Kendrick E, et al. 2005. Seasonal fluctuations in the mass of the Amazon River system and Earth's elastic response. Geophys. Res. Lett., 32(16): L16308, doi:10.1029/2005GL023491.
[5] Blewitt G, Lavallée D, Clarke P, et al. 2001. A new global model of earth deformation: Seasonal cycle detected. Science, 294(5550): 2342-2345.
[6] Boehm J, Heinkelmann R, Schuh H. 2007. Short note: A global model of pressure and temperature for geodetic applications. J. Geod., 81(10): 679-683.
[7] Bruinsma S, Lemoine J M, Biancale R, et al. 2010. CNES/GRGS 10-day gravity field models (release 2) and their evaluation. Adv. Space Res., 45(4): 587-601.
[8] Chen J L, Wilson C R, Tapley B D. 2006. Satellite gravity measurements confirm accelerated melting of Greenland Ice Sheet. Science, 313(5795): 1958-1960.
[9] Cheng M K, Tapley B D. 2004. Variations in the Earth's oblateness during the past 28 years. J. Geophys. Res., 109(B9): B09402, doi: 10.1029/2004JB003028.
[10] Dong D N, Fang P, Boek Y, et al. 2002. Anatomy of apparent seasonal variations from GPS-derived site position time series. J. Geophys. Res., 107(B4): ETG9-1-ETG9-16.
[11] Dziewonski A M, Anderson D L. 1981. Preliminary reference Earth model. Phys. Earth Planet. Inter., 25(4): 297-356.
[12] Farrell W E. 1972. Deformation of the Earth by surface loads. Rev. Geophys., 10(3): 761-797.
[13] Fu Y N, Freymueller J T. 2012. Seasonal and long-term vertical deformation in the Nepal Himalaya constrained by GPS and GRACE measurements. J. Geophys. Res., 117(B3): B03407, doi: 10.1029/2011JB008925.
[14] Fu Y N, Argus D F, Freymueller J T, et al. 2013. Horizontal motion in elastic response to seasonal loading of rain water in the Amazon Basin and monsoon water in Southeast Asia observed by GPS and inferred from GRACE. Geophys. Res. Lett., 40(23): 6048-6053.
[15] Herring T, King R, McClusky S. 2009. GAMIT reference manual and GLOBK reference manual, release 10.3. Cambridge: Massachusetts Institute of Technology.
[16] Jentzsch G. 1997. Earth tides and ocean tidal loading. Lecture Notes in Earth Sciences, 66: 145-171.
[17] Jiang W P, Li Z, Liu H F, et al. 2013. Cause analysis of the non-linear variation of the IGS reference station coordinate time series inside China. Chinese J. Geophys. (in Chinese), 56(7): 2228-2237,doi:10.6038/cjg20130710.
[18] Khan S A, Wahr J, Bevis M, et al. 2010. Spread of ice mass loss into northwest Greenland observed by GRACE and GPS. Geophys. Res. Lett., 37(6): L06501, doi: 10.1029/2010GL042460.
[19] Kusche J, Schrama E J O. 2005. Surface mass redistribution inversion from global GPS deformation and Gravity Recovery and Climate Experiment (GRACE) gravity data. J. Geophys. Res., 110(B9): B09409, doi: 10.1029/2004JB003556.
[20] Liao H H, Zhong M, Zhou X H. 2010. Climate-driven annual vertical deformation of the solid Earth calculated from GRACE. Chinese J. Geophys. (in Chinese), 53(5): 1091-1098.
[21] Luthcke S B, Arendt A A, Rowlands D D, et al. 2008. Recent glacier mass changes in the Gulf of Alaska region from GRACE mascon solutions. J. Glaciol., 54(188): 767-777.
[22] Pinel V, Sigmundsson F, Sturkell E, et al. 2007. Discriminating volcano deformation due to magma movements and variable surface loads: application to Katla subglacial volcano, Iceland. Geophys. J. Int., 169(1): 325-338.
[23] Sauber J, Plafker G, Molnia B F, et al. 2000. Crustal deformation associated with glacial fluctuations in the eastern Chugach Mountains, Alaska. J. Geophys. Res., 105(B4): 8055-8077.
[24] Shah T, Singh O P, Mukherji A. 2006. Some aspects of South Asia's groundwater irrigation economy: analyses from a survey in India, Pakistan, Nepal Terai and Bangladesh. Hydrogeol. J., 14(3): 286-309.
[25] Sun W K, Wang Q, Li H, et al. 2009. Gravity and GPS measurements reveal mass loss beneath the Tibetan Plateau: Geodetic evidence of increasing crustal thickness. Geophys. Res. Lett., 36(2): L02303, doi: 10.1029/2008GL036512.
[26] Swenson S, Wahr J. 2006. Post-processing removal of correlated errors in GRACE data. Geophys. Res. Lett., 31: L23402, doi: 10.1029/2005GL025285.
[27] Swenson S, Chambers D, Wahr J. 2008. Estimating geocenter variations from a combination of GRACE and ocean model output. J. Geophys. Res., 113(B8): B08410, doi: 10.1029/2007JB005338.
[28] Tapley B D, Bettadpur S, Watkins M, et al. 2004. The gravity recovery and climate experiment: Mission overview and early results. Geophys. Res. Lett., 31(9): L09607, doi: 10.1029/2004GL019920.
[29] Tesmer V, Steigenberger P, van Dam T, et al. 2011. Vertical deformations from homogeneously processed GRACE and global GPS long-term series. J. Geod., 85(5): 291-310.
[30] Tregoning P, Watson C, Ramillien G, et al. 2009. Detecting hydrologic deformation using GRACE and GPS. Geophys. Res. Lett., 36(15): L15401, doi: 10.1029/2009GL038718.
[31] Van Dam T M, Blewitt G, Heflin M B. 1994. Atmospheric pressure loading effects on Global Positioning System coordinate determinations. J. Geophys. Res., 99(B12): 23939-23950.
[32] Van Dam T M, Wahr J, Lavallée D. 2007. A comparison of annual vertical crustal displacements from GPS and Gravity Recovery and Climate Experiment (GRACE) over Europe. J. Geophys. Res., 112(B3): B03404, doi: 10.1029/2006JB004335.
[33] Van Dam T M, Wahr J, Chao Y, et al. 1997. Predictions of crustal deformation and of geoid and sea-level variability caused by oceanic and atmospheric loading. Geophys. J. Int., 129(3): 507-517.
[34] Wahr J, Khan S A, van Dam T M, et al. 2013. The use of GPS horizontals for loading studies, with applications to northern California and southeast Greenland. J. Geophys. Res., 118(4): 1795-1806.
[35] Wahr J, Molenaar M, Bryan F. 1998. Time variability of the Earth's gravity field: Hydrological and oceanic effects and their possible detection using GRACE. J. Geophys. Res., 103(B12): 30205-30229.
[36] Wang H S, Hsu H T, Zhu Y Z. 2002. Prediction of surface horizontal displacements, and gravity and tilt changes caused by filling the Three Gorges Reservoir. J. Geod., 76(2): 105-114.
[37] Wang H S. 2000. Surface vertical displacements and level plane changes in the front reservoir area caused by filling the Three Gorges Reservoir. J. Geophys. Res., 105(B6): 13211-13220.
[38] Wang M, Shen Z K, Dong D N. 2005. Effects of non-tectonic crustal deformation on continuous GPS position time series and correction to them. Chinese J. Geophys. (in Chinese), 48(5): 1045-1052.
[39] Yan H M, Chen W, Zhu Y Z, et al. 2009. Contributions of thermal expansion of monuments and nearby bedrock to observed GPS height changes. Geophys. Res. Lett., 36(13): L13301, doi: 10.1029/2009GL038152.
[40] Yang Y D, E D C, Chao D B, et al. 2009. Seasonal and inter-annual change in land water storage from GRACE. Chinese J. Geophys. (in Chinese), 52(12): 2987-2992.
[41] Zhang F P, Dong D N, Cheng Z Y, et al. 2002. Seasonal vertical crustal motions in China detected by GPS. Chinese Sci. Bull., 47(21): 1772-1780.
[42] Zhang Z Z, Chao B F, Lu Y, et al. 2009. An effective filtering for GRACE time-variable gravity: Fan filter. Geophys. Res. Lett., 36(17): L17311, doi: 10.1029/2009GL039459.
[43] 姜卫平, 李昭, 刘鸿飞等. 2013. 中国区域IGS基准站坐标时间序列非线性变化的成因分析. 地球物理学报, 56(7): 2228-2237,doi:10.6038/cjg20130710. 
[44] 廖海华, 钟敏, 周旭华. 2010. 利用GRACE卫星重力资料解算气候驱动的地表周年垂直形变. 地球物理学报, 53(5): 1091-1098. 
[45] 王敏, 沈正康, 董大南. 2005. 非构造形变对GPS连续站位置时间序列的影响和修正. 地球物理学报, 48(5): 1045-1052. 
[46] 杨元德, 鄂栋臣, 晁定波等. 2009. GRACE估算陆地水储量季节和年际变化. 地球物理学报, 52(12): 2987-2992. 
[47] 张飞鹏, 董大南, 程宗颐等. 2002. 利用GPS监测中国地壳的垂向季节性变化. 科学通报, 47(18): 1370-1377.