2015, Vol. 58
2. 中国科学院大学, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
A reduction of hydrological signals from the detected integral secular trends using global hydrological models (Global Land Data Assimilation System, GLDAS and Climate Prediction Center, CPC, averaged) is attempted. The glacier model result provided by Paulson is used to reduce the GIA (Glacial Isostatic Adjustment) effect. Also, we use the scaling factor method to weaken the GRACE post-process errors. At last we obtain the secular trends of gravity anomalies in China and adjacent areas. The results are largely similar to the previous research. Also, we find a remarkable increase signal of gravity in the Tibetan Plateau. Such a change cannot be attributed to the lake expanding and precipitation increase in the Tibetan Plateau. Instead it is possibly caused by the crustal thickening due to the subduction of Indian plate beneath the Eurasian plate.
The continuing collision and compression between the India and Eurasia plates result in intense uplift and orogenic extension of the Tibetan Plateau. We built a forward model with 3D rectangular prisms according to the GPS results and Airy isostatic hypothesis. Using the Bangong-Nujiang suture zone as the boundary, the plateau was divided into southern and northern parts. The modeling results show the gravity anomalies are increasing at a rate of 0.2 μGal·a-1 in the north, which is a little less than the result from GRACE solution 0.3±0.08 μGal·a-1 (corresponding to a crustal thickening rate of 3 mm·a-1). The rest unexplained part is possibly related to lakes, glacier, permafrost and others. Meanwhile, our results show the southern part has a rate of -0.62±0.06 μGal·a-1, which is presumably the combined influence of tectonics and groundwater and so on. It is a challenging to explain the positive signal in the Tibetan Plateau clearly with current knowledge. The new GRACE Follow-On satellites may provide us more means and evidence. With the continuing accumulation of observations, we can combine a variety of modern geodetic techniques (GRACE, GPS, absolute gravity, remote sensing, and so on) to study the geodynamics of the Tibetan Plateau more accurately.
青藏高原(25°N—40°N，74°E—104°E)平均海拔高度在4000 m以上，有“世界屋脊”和“第三极”之称.印度板块的持续向北移动通过折叠和各种断层错动机制形成了现今高原地壳的形状(孙文科等，2011)，沿着喜马拉雅山的高频率地震活动反映出在印度和亚洲之间的板块碰撞仍在推着青藏高原向上抬升.Liang等(2013)根据青藏高原及周边地区750个GPS测站十多年的观测数据获取了现今青藏高原地壳运动三维速度场，垂直速度场显示整个青藏高原作为整体相对于北部稳定块体一直在持续隆升.众多证据(Tapponnier et al.，2001; Gan et al.，2007; Royden et al.，2008)均表明青藏高原仍在继续变形.近些年，空间大地测量技术在研究青藏高原动力学过程中发挥了重要的作用，特别是于2002年3月升空由美德联合开发旨在获取地球重力场中长波部分及全球重力场的时变特征的 GRACE(Gravity Recovery and Climate Experiment)卫星，其提供的高精度观测数据为研究全球物质分布及季节性变化提供了重要参考依据(Cazenave and Chen，2010)，GRACE卫星重力技术已经广泛应用于地壳形变(王林松等，2014)，陆地水变化(Hu et al.，2006)，极地冰川变化(Chen et al.，2006; Velicogna and Wahr，2005)，海平面变化等方面(Woodworth and Gregory，2003; Lombard et al.，2007)，在研究青藏高原质量迁移变化及亚洲高山冰雪融化等方面也 取得了重要进展(Matsuo and Heki，2010). GRACE重力卫星发射至今已超过11个年头，过去GRACE卫星重力资料的周期比较短，主要利用GRACE研究大尺度季节性信号变化，随着资料的累积和数据处理方法的改进与成熟，用GRACE卫星研究长期变化成为新的热点.
关于青藏高原重力长期变化国内外学者开展了大量工作.对于青藏高原南部重力下降归结为最近几十年喜马拉雅冰川消融(Matsuo and Heki，2010)和印度北部地区的地下水抽采(Tiwari et al.，2009).然而不同的研究结果存在分歧.Matsuo和Heki(2010)最早采用卫星重力技术研究亚洲高山时变冰川消融情况，结果显示2003—2009冰消融率为-47±12 Gt/a，等效于0.13±0.04 mm·a-1的海平面升高；Jacob等(2012)研究了2003—2010全球冰川冰盖对海平面变化的贡献，显示亚洲高山冰川消融率约为-4±20 Gt/a，对于研究区域的划分及对地下水考虑的不同是造成结果与Matsuo和Heki(2010)差异的主要原因(Yi and Sun，2014).Ye等(2009)基于遥感卫星技术采用多时相网格方法定量分析青藏高原地区的冰川变化情况，结果显示最近几十年喜马拉雅冰川既有消融同时也有增长，且消融占主导地位.Sun等(2009)基于绝对重力和GPS技术通过对拉萨、昆明、大理三个台站进行研究，用大地测量的证据揭示青藏高原地下质量损失，得到三个测站的平均重力变化率为-0.78±0.47 μGal·a-1，认为青藏高原的主要构造特征是连续变形，且该连续变形基本上限定在高原本身并被地壳增厚吸收，进而引起地壳垂直运动；Sun等(2011)通过对拉萨测站不同高斯滤波半径的模拟研究，计算了青藏高原的空间重力分布，结果显示重力变化率依赖于高斯半径，同时也讨论了引起拉萨重力变化的物质来源，认为该地区重力减少的主要原因在于印度板块碰撞构造形变而非冰雪融化.
目前，关于青藏高原中北部重力上升信号开展了一些工作，却还没有合理的解释.Jacob等(2012)结果中显示青藏高原及祁连地区大致有7±7 Gt/a的质量增加趋势；Zhang等(2013)认为该青藏高原的质量增加信号可能与青藏高原湖水位的升高有关而不仅仅是冰川物质平衡，并通过ICEsat卫星数据予以佐证； Yi和Sun(2014)在对亚洲高山冰川变化进行评估时也发现这一明显的质量增加信号，认为Jacob等(2012)及Zhang等(2013)给出的量级(+7 Gt/a，+8.06 Gt/a)不足以用来解释高原明显的质量增加(+30 Gt/a)，且认为目前要解释清楚该质量增加信号是有些困难的.Liang等(2013)得到的现今青藏高原地壳运动三维速度场显示印度板块向北挤压速率为4 cm·a-1，青藏高原整体大约以2 mm·a-1的速率隆升，且南部上升速率约为3~4 mm·a-1.印度板块的挤压使得青藏高原地壳形变、物质分布变化可能是导致青藏高原重力上升的主要因素，而冰雪消融不是青藏高原北部重力上升主要原因.
本文对最新的GRACE(RL05)数据采用尺度因子的方法减少数据处理过程中误差的影响，最终基于最小二乘计算方法得到2003—2013中国及周边地区长期重力异常变化.依据青藏高原地壳结构和GPS观测结果，考虑地壳均衡和构造因素构建了地壳形变模型，并进行正演模拟分析，计算结果对于认识青藏高原隆升动力学有一定参考意义.
2 数据处理为增强结果的可靠性，这里同时采用CSR、GFZ和JPL这三家机构提供的Level-2(Rlease-05)数据，并对三家结果进行平均作为最终的估计.其中CSR提供的最新RL05模型为96阶次，将其截断至与GFZ一致的90阶，将JPL提供的RL05模型统一截断至60阶次，共计124个月份，其他月份因加速度数据或者K波段测距数据丢失.最新的RL05数据已扣除各种潮汐、固体潮、非潮汐大气、海洋大气等的影响，进一步降低了南北条带误差，东西方向上的误差有了一定程度的提高(Bettadpur，2012)，因而该数据为本文的研究提供了很可靠的保障.
2.1 GRACE卫星数据处理GRACE卫星对于低阶项不敏感，这里分别采用Swenson等(2008)提供的一阶项，和卫星激光测距SLR(Cheng and Tapley，2004)得到的二阶项来代替球谐系数中的相应项；接着对所有月份的球谐系数求取平均值，各个月份的球谐系数与该平均值做差得到异常场；虽然最新的RL05数据已经进一步考虑了南北向条带的影响，通过研究发现大部分地区仍存在轻微的条带影响，故这里选择P5M15策略对高于15阶次的奇偶项分别采用5阶多项式进行拟合并予以扣除以此来削弱南北条带误差的影响；采用200 km的扇形滤波(Zhang et al.，2009)来削弱高阶球谐系数的噪声影响；球谐系数最终可以转换为分辨率为1°×1°的全球分布的重力异常、等效水柱高或大地水准面起伏等变化(Wahr et al.，1998).
2.2 GRACE卫星得到长期性重力异常变化基于上面提到的方法理论，我们通过对2003—2013这11年间GRACE球谐系数产品进行分析，利用下面公式可得到每个格网点的重力变化值(Ray et al.，2003)：
式中θ和λ分别表示地球上某一点的余纬和经度，Δgk为地球上某点的重力异常变化值，GM为地球引力常数，r为地球半径，扇形滤波Wn和Wm表示阶n和次m方向上的高斯系数，ΔCnm和ΔSnm为每个月重力场球谐系数与平均月重力场模型球谐系数均值之差，
nm表示阶数为n，次数为m的归一化缔合勒让德函数.
关于长期性变化趋势，Ray等(2003)通过对卫星重力潮汐模型的研究发现S2，K2，K1潮汐的影响造成的混叠分别会造成161天，3.7年和7.4年周期的影响，因为本文研究的时间尺度为11年，故这里暂不考虑K1潮汐引起的7.4年周期.对每个格网点得到的重力异常值通过最小二乘算法进行拟合，参考下面的公式(Steffen et al.，2009)：
式中A代表常数项，B代表长期项年变化率，t为重力场模型的时间，周期项ωi的振幅为
，其中i=1和i=2分别表示年周期项和半年周期项，i=3为与S2半日波相关的161天周期项，i=4表示与K2相关的3.7年周期项，ε表示拟合残差.
GRACE卫星数据广泛用来监测地球上的时变质量变化，但获得的数据不可避免的存在误差，通过相关滤波及去相关(Han et al.，2005)的方法可以削弱相应误差，但是在削弱误差的同时也会造成真实信号的衰减，这就需要在数据处理过程中选择一个折中：信号损失最小同时噪声减弱达到最大，这里通过尺度因子的方法来考虑.Landerer和Swenson(2012)对GRACE获得的陆地水估计的精度进行了分析，重点考虑了基于流域平均的单一尺度因子，基于单个格网点的多尺度因子，和基于格网点时间频率函数的尺度因子，并进行了详细的讨论，这里采用较为常用的第一种策略得到研究区域的单一尺度因子，具体做法为：对水文模型时间序列ΔST数据进行球谐展开，并截断至90阶，采用与GRACE数据一样的数据处理，去相关、扇形滤波，得到研究区域的时间变化序列ΔSF，最后再基于最小二乘根据下面公式得到尺度因子k.
由GRACE得到的时间序列乘以该尺度因子k作为研究区域最终的估计.
2.4 GRACE误差考虑由GRACE来估计地球表面质量变化的误差(Wahr et al.，2006)大致可以分为两大类：对球谐系数的处理等误差引起的和由其他信号泄露引起的.对于第一类误差主要包括：仪器误差、数据处理误差、混淆误差.混淆误差可以理解为重力场中那些次月级的短周期变化混入月平均重力场而引起，通常利用相应的模型来扣除，GRACE数据处理中，即利用模型化的数据移去固体地球、海洋潮汐、大气质量变化等影响，这些模型误差都会与重力场模型存在混淆误差.当研究区域大于目标区域时，得到的结果不仅受到所选区域内质量变化的影响，同时也受到此区域外信号的影响，即所谓信号的泄露.关于信号的泄露，Baur等(2009)通过对格陵兰岛的冰川质量变化进行研究构建了一套扣除泄露效应的方法，并给出三个不同机构(CSR，GFZ，JPL)得到的冰川质量年变化率.
这里对GRACE误差的考虑主要从两方面，其一，在研究长期性变化中对每个格网点的重力异常值进行了最小二乘拟合，将拟合的残差值作为每个格网点的误差，得到全球或者研究区域的RMSE值；其二，GRACE得到的球谐系数已经扣除了海洋大气的影响，理论上海洋地区质量变化值为0，但是由于误差的影响海洋地区仍旧有一些残余信号存在，可以将研究区域转化为同纬度的海洋地区，得到其变化值作为研究区域的误差.
2.5 水文模型本文采用两个水文模型GLDAS(Rodell et al.，2004)与CPC(Fan and van den Dool，2004)分别来估计青藏高原地区土壤水及雪水的变化.GLDAS水文模型由美国宇航局哥达航空中心(Goddard Space Flight Center，NASA)和美国国家 环境预报中心(NCEP，National Centers of Environmental Prediction)共同构建的全球水文模型，其主要通过近实时的地面和空间数据来约束其模型来获取陆地表面变化的近实时信息.这里选择其提供的1°×1°格网间隔的产品，采用GLDAS_Noah模型提供的土壤湿度(0~0.1 m，0.1~0.4 m，0.4~1 m和1~2 m)及雪水数据来计算青藏高原地区陆地水变化.CPC水文模型是由美国国家海洋和大气管理局气象预报中心同化系统提供的数据集，时间间隔为1 个月，空间分辨率为0.5°×0.5°，土壤水深度至1.6 m. 为增强结果的可靠性，这里分别对这两个水文模型(GLDAS，CPC)进行与GRACE一样的数据处理，即先将格网点数据截断转化为90阶的球谐系数，然后进行去平均和扇形滤波处理，再将得到的两个水文模型进行平均作为该地区土壤水变化的估计，将两个水文模型的差值作为误差估计.
2.6 关于降雨降雨资料数据集是基于国家气象信息中心基础资料专项最新整编的中国地面约2400台站降水资料，利用ANUSPLIN软件的薄盘样条法(TPS，Thin Plate Spline)进行空间插值，生成1961年至最新的水平分辨率0.5°×0.5°的中国降水月值格点数据.为了更清晰地比较降雨较之前是否有明显的增多趋势，我们选择1990—2013年共24年降雨资料，分析了青藏高原地区的月降雨量变化和年降雨量变化情况，见图 1.
在图 1中并未看出青藏高原地区降雨有明显的增加趋势.同时也采用国家气候中心提供的160测站数据，分析了几个测站(拉萨，玉树，昌都等)的降雨变化，通过近63年的降雨资料，我们均未发现这些年明显的降雨增多趋势，故在后面结果讨论中并没考虑降雨因素的影响.
 | 图 1 青藏高原地区降雨量变化图 (a) 月降雨量变化图； (b) 年降雨量变化图. Fig. 1 The precipitation variations of Tibet Plateau |
冰川均衡调整(GIA)是黏弹地球对末次冰进期和冰退期地表冰和海水负荷改变的影响，是一种重要的地球动力学现象(汪汉胜等，2009)，GIA过程主要体现在地幔物质的流动、地壳的运动和地球重力场的变化.关于GIA的考虑，这里采用Paulson等(2007)提供的冰川模型结果予以扣除.
3 结果与讨论通过上文提到的理论基于GRACE资料得到2003—2013年间中国及周边地区重力长期变化趋势，如图 2.
 | 图 2 由GRACE得到的中国及周边地区重力异常长期变化趋势 Fig. 2 Secular trends of gravity anomaly in China and adjacent area from GRACE solutions |
图 2是对三家机构(CSR，GFZ，JPL)平均之后得到的结果，通过三家机构与平均值做差发现研究区域异常值最大仅为0.02 μGal·a-1，这说明三家机构在该地区得到的结果比较相似，结论比较可靠.在上面结果中我们可以看出在华北平原、天山附近、印度北部及喜马拉雅山区等地方有显著的重力减少趋势，同时在青藏高原祁连山周围存在明显的重力增加信号.
如前文所述，GRACE在数据处理中削弱误差的同时也造成了自身信号的衰减，这里通过尺度因子的方法来对结果进行恢复.首先将对GLDAS_Noah水文模型格网数据展开至90阶次，通过与GRACE一样的数据处理，基于最小二乘得到该地区的尺度因子为1.25，故GRACE得到的结果都乘以该因子来作为研究区域最终的估计.接着从GRACE得到的结果中扣除土壤水及雪水的影响，为增加结果的可靠性，本文对水文模型GLDAS_Noah和CPC得到的结果进行平均作为该地区土壤水变化的估计；关于降雨因素，通过对青藏高原地区月降雨量和年降雨量的分析并未发现这些年明显的降雨增加趋势，故这里未考虑降雨的影响；对于GIA的考虑，我们选择Paulson等(2007)提供的冰川模型结果予以扣除，该模型主要考虑在北欧、北美及南极等地区，结果显示青藏高原地区由于GIA的影响造成重力变化约为0.05 μGal·a-1.在扣除土壤水及GIA因素之后，可以得到青藏高原及周边地区2003—2013年重力异常长期性变化情况如图 3.
 | 图 3 GRACE结果中扣除土壤水及GIA影响之后得到的中国及周边地区重力异常长期变化趋势 Fig. 3 Secular trends of gravity anomaly in China and adjacent area after the deduction of soil moisture and GIA effects from GRACE solutions |
在图 3中，我们依旧可以看出中国及周边地区重力异常长期性变化一些明显的增加及减少趋势.对于重力增加区域如青藏高原和重力减少区域如华北平原、天山附近、印度北部等地区，首先提取出研究区域所有格网点坐标，接着对所有格网点采用纬度加权的方法求取平均值，最终得到该区域的估值.对于估值的不确定性，这里参考Chen等(2009)选择与研究区域同纬度海洋地区的相应变化值.
关于那些质量减少区域，许多学者已进行了相关研究和探讨.通过计算，我们得到华北地区 2003—2013重力异常大致以-0.54±0.09 μGal·a-1(-8±1.1 Gt/a)的速率递减，与Feng等(2013)得到2003—2010华北平原地下水损耗的-8.3±1.1 Gt/a较为一致.天山地区2003—2013重力异常以-0.6±0.1 μGal·a-1(-3.5±0.6 Gt/a)的速率递减，略小于Jacob等(2012)得到的结果-5±6 Gt/a，结果与研究时间范围及研究区域的面积有较大的关系.印度北部地区重力异常以-1.17±0.18 μGal·a-1(-28±3 Gt/a)的速率递减，与Yi和Sun(2014)最新的结果为-30.6±5 Gt/a较为一致.
在图 3中还可以看出青藏高原明显的质量增加信号，Jacob等(2012)在研究全球冰川冰盖对海平面变化贡献时发现该地区约以7±7 Gt/a的速率在增加.Zhang等(2013)认为该质量增加信号可能与青藏高原湖水位的上升有关，通过对占青藏高原53%的湖区近4~7年的ICESat数据分析，结果显示湖水面平均上升速率为0.14 m·a-1，等效为4.95 Gt/a.然而Yi和Sun(2014)的结果则显示青藏高原的重力增加为+30 Gt/a，远高于Jacob等(2012)得到的7 Gt/a，且认为目前要解释清楚该正信号是有些困难的.对于Yi和Sun(2014)得到30 Gt/a的质量增加信号，我们认为其结果与研究区域的面积有较大的关系，由于其并未给出该地区详细的划分，依照图 3所示勾勒出青藏高原(约160万km2)，本文的结果显示质量增加约为12.5±3 Gt/a(0.3±0.08 μGal·a-1)，故Yi和Sun(2014)可能高估了青藏高原的质量增加.Zhang等(2013)的结果显示由于湖泊水位升高造成青藏高原质量增加约为8 Gt/a，我们通过对青藏高原地区湖泊的分布区域及形态研究发现其并不能与GRACE得到的质量增加信号区域很好的拟合.结合这段时间降雨资料也并未发现该地区明显的降雨增多趋势，对于湖泊水位升高的来源也不能很好地解释.青藏高原地区湖区面积(马荣华等，2011)约为41831.7 m2，仅占整个青藏高原地区的1.6%，我们认为由湖泊水位上升引起的质量增加(8 Gt/a)不能很好地解释青藏高原的质量增加信号(12.5±3 Gt/a)，故接下来尝试通过构建地壳形变模型对青藏高原的质量增加信号进行相关探讨，给出青藏高原每年重力异常变化量，并对重力增加的原因进行解释.
4 正演模型的构建印度板块与欧亚板块的持续碰撞与挤压，使得青藏高原以强烈的隆升和造山伸展变形方式活动(邢乐林等，2011)，且该活动是个持续平稳的过程，强调地壳的双倍增厚是高原隆升的主导机制，也是目前认为比较合理和被普遍接受的隆升机制(刘燊等，2001).在挤压过程中，青藏高原地壳隆升且不断地增厚分别形成了构造山脉和山根，低密度的山根挤出了高密度的地幔物质导致区域内物质亏损.在前面GRACE结果的分析中，可以清晰地看出在青藏高原内陆地区有明显的质量增加信号(12.5±3 Gt/a)，对于该正信号目前还存在较大的争议，我们认为由湖泊水位上升引起的质量增加不能很好地解释该正信号，尝试从构造的角度通过正演模型来进行分析讨论，关于挤压碰撞过程可以通过简单的动力学模型示意图予以说明.
如图 4，印度板块碰撞青藏高原会引起青藏高原的隆升，在抬升的同时也会引起底部地壳的增厚(Tapponnier et al.，2001).这里假设青藏高原平均厚度为70 km，Liang等(2013)的结果显示青藏高原整体大约以2 mm·a-1的垂直速率一直在隆升，且南部上升速率约为3~4 mm·a-1，显示青藏高原地壳上升速率不同.张燕等(2013)基于全新的 1 ∶ 100万区域重力工作成果认为青藏高原重力异常形态是具有多条结合带的拼合体，班公湖—怒江结合带是高原内最主要的重力高异常带，将不同深度层次的重力场分为截然不同的南北两大区块.丁志峰等(2001)，吴庆举和曾融生(1998)，苏伟等(2002)通过地震波等手段研究青藏高原的地壳结构，结果显示青藏高原以班公湖—怒江缝合带为界，南北分带特征明显，南部是高速、冷的岩石圈上地幔；北部是低速、热的岩石圈上地幔.青藏高原的Moho界面在班公湖—怒江缝合带附近存在明显的深度错断，并认为该错断现象一个可能的解释是班公湖—怒江缝合带是印度地壳向欧亚下地壳挤入的前沿.这些结果显示青藏高原地壳以班公湖—怒江缝合带为界，南、北存在着明显的差异使得南、北地壳形变速率也不相同.南北速度结构上的不同可能与青藏高原的形成过程有关，现今地壳形变在青藏高原南部隆升的同时地壳增厚明显，而在北部区域只有微量的隆升，地壳增厚幅度很小.为简化地壳形变模型，依据青藏高原的地壳结构以班公湖—怒江断裂带为界，将青藏高原分为南、北两部分，青藏高原南部地区地壳隆升的同时又有增厚，而北部地区地壳以隆升为主.
 | 图 4 印度板块碰撞欧亚板块动力学模型示意图 Fig. 4 Dynamical model of the collision between India Plate and Eurasia Plate |
Hetényi等(2007)讨论了青藏高原下部印度板 块的密度分布，显示该地区地下0~20 km密度约为2.67 g·cm-3，20~35 km密度约为2.74 g·cm-3，35~60 km密度约为2.9 g·cm-3，地幔密度为3.3 g·cm-3.我们按各层所占比例加权平均得到青藏高原地壳平均密度为2.78 g·cm-3.再根据艾黎均衡模型原理依照公式(4)，得到山根与地壳增厚的比例r/h=ρc/Δρ=5.35.
基于GPS结果假设高原南部隆升速率为 4 mm·a-1，那么底部增厚的速率约为2.1 cm·a-1； 北部主要表现为微量的隆升，隆升速率为2 mm·a-1，底部增厚速率为0 mm·a-1.地壳的上升和增厚均以长方体薄板来表示，通过直立长方体模型予以正演模拟分析.长方体是位场正演计算中最为常用的三度体模型，可以用平行于直角坐标面的平面对任意形状的三维地质体进行分割，将三维复杂地质形体重力异常的正演计算转化为一系列长方体元重力异常的叠加.考虑如下：以班公湖—怒江为界，提取出整个青藏高原地区所有1°×1°格网点的经纬度坐标：其中南部共有160个格网点，北部共179个格网点，共计339个格网点.以格网点的经纬度作为长方体的中心，以该中心点±0.5°作为长方体的长宽边界，以青藏高原南部地壳上升4 mm·a-1，地壳增厚21 mm·a-1；北部地壳上升2 mm·a-1分别作为每个长方体单元的高.对于地壳上升参考面选择大地水准面为基准，对于地壳底部增厚参考面选择地壳平均厚度70 km为基准.地壳上升的单元密度为2.78 g·cm-3，增厚单元密度为-0.52 g·cm-3(地幔密度为3.3 g·cm-3).采用下面公式得到该长方体对任意格网点的重力异常变化(王谦身等，2003).
公式(5)中(a1，a2)，(b1，b2)，(h1，h2)分别为长方体在三个方向的边界，p(x，y，z)为长方体外任一点，ρ为长方体密度，q(ε，η，τ)为长方体内一点元.我们对所有小长方体进行正演计算，然后对重力效应叠加进而得到印度板块碰撞引起青藏高原地壳形变产生的重力变化，见图 5.
 | 图 5 (a) 直立长方体正演模型得到的青藏高原重力长期性变化； (b) GRACE扣除土壤水、GIA影响得到的结果 Fig. 5 (a) Secular trends of gravity anomaly in Tibetan Plateau from the forwarding model；(b) Results after the deduction of soil moisture and GIA effects from GRACE solutions |
图 5中白色实线为班公湖—怒江断裂带，其中图 5a为直立长方体正演模型得到的重力变化结果，显示班公湖—怒江断裂带以北地区重力异常大致约以+0.2 μGal·a-1的速率在递增，断裂带以南地区重力异常约以0.09 μGal·a-1的速率在增加.该结 果是对我们模拟地壳隆升与增厚的综合反映，与所选直立长方体正演参数有关，针对青藏高原，地壳上升2 mm·a-1会引起重力异常约0.2 μGal·a-1 的增加.图 5b为GRACE在扣除土壤水、GIA影响之后得到的结果，显示青藏高原(图中折线所围)重力异常约以0.3±0.08 μGal·a-1的速率变化，与图 5a正演结果比较可以看出，正演模拟结果只占到GRACE结果的约2/3左右，对于其余未能解释部 分，可能与湖水、冰川因素(Lei et al.，2014)，冻土因素(罗栋梁等，2012)等有关，也可能是这些因素的综合影响.如果GRACE观测到的质量变化(0.3 μGal·a-1)全部由地壳隆升引起，那么对应地壳隆升的速率约为3 mm·a-1.
同时在图 5b可以看出，由GRACE扣除GIA效应和水文模型的结果中显示青藏高原南部地区重力异常大致以-0.62±0.06 μGal·a-1的速率在递减，该结果是包括构造及地下水等因素的综合影响，
而由正演模型(图 5a)得到青藏高原南部地区重力异常仅以0.1 μGal·a-1的速率在递增，说明构造因素造成的重力变化较地下水冰雪等因素的影响处于次要地位.Sun等(2009)在用大地测量证据揭示青藏高原内部质量损失，进行了相关模型研究，结果显示青藏高原底部增厚速率约2.3±1.3 cm·a-1，与本文得到的2.1 cm·a-1较为一致.邢乐林等(2011)采用了14年绝对重力资料及GPS研究成果，从大地测量的角度检测了青藏高原拉萨点地壳的增厚，速率为3.9±0.8 cm·a-1，并认为随着更多监测系统的建立，越来越多的观测资料可用于研究青藏高原的地球动力学特征.
5 结论与讨论本文通过对三家机构(CSR，GFZ，JPL)提供的 2003—2013最新版本RL05数据进行处理分析，基于尺度因子的方法减少了数据处理中的误差，采用水文模型(GLDAS和CPC的平均值)扣除土壤水的影响，采用Paulson等(2007)提供的冰川模型结果扣除GIA影响，最终得到中国及周边地区长期性重力变化.在结果中可以看到一些明显的重力减少和增加区域，对于重力减少的区域(Matsuo and Heki，2010; Sun et al.，2009; Feng et al.，2013; Jacob et al.，2012)已经进行过较多的分析.在青藏高原有较为明显的质量增加信号，Zhang等(2013)认为此增加可能与该地区湖水面的升高有关.我们通过对同时间范围内降水资料的分析并没有发现明显的降雨增加趋势，且由湖泊水位升高引起的质量增加(8 Gt/a)不能很好地解释GRACE卫星得到的青藏高原内陆地区的质量变化(12.5±3 Gt/a).故认为湖水位的上升引起的重力增加不能很好解释该区域内重力上升.基于青藏高原一直在继续隆升变形这一事实，通过直立长方体模型进行正演模拟研究，以班公湖—怒江断裂带为界将青藏高原划分为南北两大区块，得到青藏高原地区的重力异常变化情况.结果显示青藏高原大致以0.2 μGal·a-1的速率在递增，小于GRACE得到的0.3±0.08 μGal·a-1的速率，我们认为剩余未解释部分可能与湖水、冰川因素、冻土因素等有关，也可能是这些因素的综合影响.目前关于青藏高原卫星重力上升信号的解释存在较大争议，且纯粹依据水质量增加是不能解释的.本文这里提供了一个思路，从构造的角度尝试对青藏高原质量增加信号进行解释.要想具体分离出造成该质量增加信号的原因，还需要更多高新技术的发展，如新一代重力卫星的发射可能会提供更多的手段和证据.
论文对青藏高原质量增加信号是在扣除水文影响与GIA影响之后归结为构造形变因素的一种尝试性考虑.由于资料有限我们的地壳形变模型较为粗略，模拟结果能够显示整体重力变化状态，对细节的反映不够.青藏高原地区地广人稀，要想获得整个地区的水文资料还是比较困难的，这里对通用的两个水文模型GLDAS与CPC进行平均作为对整个地区的估计，难免存在误差.如果有更多的实测资料加入，结果的可靠性会大大增强.冻土是一个较复杂的过程，温度升高会导致冻土融化，但是冻土层又具备一定的蓄水能力(Yi and Sun，2014)，关于冻土的研究还有较大的不确定性本文未予以考虑.随着观测资料的持续累积，我们可以结合多种现代大地测量手段(如GRACE，GPS，绝对重力，超导重力，遥感，测高卫星等)采用数据融合技术，更加精确地分析青藏高原的动力学问题.
致谢 感谢两位匿名评审专家和编辑提出的宝贵建议和帮助；感谢CSR，GFZ，JPL提供的GRACE RL05时变重力场模型数据；文中所涉图件大多由GMT软件绘制.
| [1] | Baur O, Kuhn M, Featherstone W E. 2009. GRACE-derived ice-mass variations over Greenland by accounting for leakage effects. J. Geophys. Res., 114(B6): B06407. |
| [2] | Bettadpur S. 2012. UTCSR Level-2 Processing Standards Document for Product Release 05, GRACE 327-742, Center for Space Research, The University of Texas at Austin. |
| [3] | Cazenave A, Chen J L. 2010. Time-variable gravity from space and present-day mass redistribution in the Earth system. Earth and Planetary Science Letters, 298(3-4): 263-274. |
| [4] | Chen J L, Wilson C R, Blankenship D D, et al. 2006. Antarctic mass rates from GRACE. Geophys. Res. Lett., 33(11): L11502. |
| [5] | Chen J L, Wilson C R, Tapley B D, et al. 2009. 2005 drought event in the Amazon River basin as measured by GRACE and estimated by climate models. J. Geophys. Res., 114, B05404, doi:10.1029/2008JB006056. |
| [6] | Cheng M K, Tapley B D. 2004. Variations in the Earth's oblateness during the past 28 years. J. Geophys. Res., 109(B9): B09402. |
| [7] | Ding Z F, He Z Q, Wu J P, et al. 2001. Research on the 3-D seismic velocity structures in Qinghai-Xizang Plateau. Earthquake Research in China (in Chinese), 17(2): 202-209. |
| [8] | Fan Y, van den Dool H. 2004. Climate Prediction Center global monthly soil moisture data set at 0.5° resolution for 1948 to present. J. Geophys. Res.,109, D10102, doi:10.1029/2003JD004345. |
| [9] | Feng W, Zhong M, Lemoine J M, et al. 2013. Evaluation of groundwater depletion in North China using the gravity recovery and climate experiment (GRACE) data and ground-based measurements. Water Resour. Res., 49(4): 2110-2118. |
| [10] | Gan W J, Zhang P Z, Shen Z K, et al. 2007. Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements. J. Geophys. Res.: Solid Earth (1978—2012), 112, B08416, doi:10.1029/2005JB004120. |
| [11] | Han S C, Shum C K, Jekeli C, et al. 2005. Non-isotropic filtering of GRACE temporal gravity for geophysical signal enhancement. Geophys. J. Int., 163(1): 18-25. |
| [12] | Hetényi G, Cattin R, Brunet F, et al. 2007. Density distribution of the India plate beneath the Tibetan plateau: geophysical and petrological constraints on the kinetics of lower-crustal eclogitization. Earth and Planetary Science Letters, 264(1-2): 226-244. |
| [13] | Hu X G, Chen J L, Zhou Y H, et al. 2006. Seasonal water storage change of the Yangtze River basin detected by GRACE. Science in China Series D, 49(5): 483-491. |
| [14] | Jacob T, Wahr J, Pfeffer W T, et al. 2012. Recent contributions of glaciers and ice caps to sea level rise. Nature, 482(7386): 514-518. |
| [15] | Landerer F W, Swenson S C. 2012. Accuracy of scaled GRACE terrestrial water storage estimates. Water Resources Research, 48(4): W04531. |
| [16] | Lei Y B, Yang K, Wang B, et al. 2014. Response of inland lake dynamics over the Tibetan Plateau to climate change. Climatic Change, 125(2): 281-290. |
| [17] | Liang S M, Gan W J, Shen C Z, et al. 2013. Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements. J. Geophys. Res., 118(10): 5722-5732. |
| [18] | Liu S, Chi X G, Li C, et al. 2001. The summarizing for the forming and uplifted mechanism of Qinghai-Tibet Plateau. World Geology (in Chinese), 20(2): 105-112. |
| [19] | Lombard A, Garcia D, Ramillien G, et al. 2007. Estimation of steric sea level variations from combined GRACE and Jason-1 data. Earth and Planetary Science Letters, 254(1-2): 194-202. |
| [20] | Luo D L, Jin H J, Lin L, et al. 2012. Degradation of permafrost and cold-environments on the interior and eastern Qinghai Plateau. Journal of Glaciology and Geocryology (in Chinese), 34(3): 538-546. |
| [21] | Ma R H, Yang G S, Duan H T, et al. 2011. China's lakes at present: Number, area and spatial distribution. Sci. China: Earth Sci., 54(2): 283-289. |
| [22] | Matsuo K, Heki K. 2010. Time-variable ice loss in Asian high mountains from satellite gravimetry. Earth and Planetary Science Letters, 290(1-2): 30-36. |
| [23] | Paulson A, Zhong S J, Wahr J. 2007. Inference of mantle viscosity from GRACE and relative sea level data. Geophys. J. Int., 171(2): 497-508. |
| [24] | Ray R D, Rowlands D D, Egbert G D. 2003. Tidal models in a new era of satellite gravimetry. Space Science Reviews, 108(1-2): 271-282. |
| [25] | Rodell M, Houser P R, Jambor U, et al. 2004. The global land data assimilation system. Bulletin of the American Meteorological Society, 85(3): 381-394. |
| [26] | Royden L H, Burchfiel B C, van der Hilst R D, et al. 2008. The geological evolution of the Tibetan Plateau. Science, 321(5892): 1054-1058. |
| [27] | Steffen H, Petrovic S, Müller J, et al. 2009. Significance of secular trends of mass variations determined from GRACE solutions. Journal of Geodynamics, 48(3-5): 157-165. |
| [28] | Su W, Peng Y J, Zheng Y J, et al. 2002. Crust and upper mantle shear velocity structure beneath the Tibetan Plateau and adjacent areas. Acta Geoscientia Sinica (in Chinese), 23(3): 193-200. |
| [29] | 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. |
| [30] | Sun W K, Hasegawa T, Zhang X L, et al. 2011. Effects of Gaussian filter in processing GRACE data: Gravity rate of change at Lhasa, southern Tibet. Sci. China: Earth Sci., 54(9): 1378-1385. |
| [31] | 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. |
| [32] | Tapponnier P, Xu Z Q, Roger F, et al. 2001. Oblique stepwise rise and growth of the Tibet Plateau. Science, 294(5547): 1671-1677. |
| [33] | Tiwari V M, Wahr J, Swenson S. 2009. Dwindling groundwater resources in northern India, from satellite gravity observations. Geophys. Res. Lett., 36(18): L18401. |
| [34] | Velicogna I, Wahr J. 2005. Greenland mass balance from GRACE. Geophys. Res. Lett., 32(18): L18505. |
| [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-30225. |
| [36] | Wahr J, Swenson S, Velicogna I. 2006. Accuracy of GRACE mass estimates. Geophys. Res. Lett., 33(6):L06401, doi:10.1029/2005GL025305. |
| [37] | Wang H S, Wu P, Xu H Z. 2009. A review of research in glacial isostatic adjustment. Progress in Geophys. (in Chinese), 24(6): 1958-1967, doi: 10.3969/j.issn.1004-2903.2009.06.005. |
| [38] | Wang L S, Chen C, Zou R, et al. 2014. Using GPS and GRACE to detect seasonal horizontal deformation caused by loading of terrestrial water: A case study in the Himalayas. Chinese J. Geophys. (in Chinese), 57(6): 1792-1804, doi: 10.6038/cjg20140611. |
| [39] | Wang Q S, An Y L, Zhang C J,et al. 2003. Gravitology (in Chinese). Beijing: Seismological Press, 163-165. |
| [40] | Woodworth P L, Gregory J M. 2003. V: sea level: Benefits of GRACE and GOCE to sea level studies. Space Science Reviews, 108(1): 307-317. |
| [41] | Wu Q J, Zeng R S. 1998. The crustal structure of Qinghai-Xizang Plateau inferred from broadband teleseismic waveform. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese), 41(5): 669-679. |
| [42] | Xing L L, Sun W K, Li H, et al. 2011. Present-day crust thickness increasing beneath the Qinghai-Tibetan Plateau by using geodetic data at Lhasa station. Acta Geodaetica et Cartographica Sinica (in Chinese), 40(1): 41-44. |
| [43] | Ye Q H, Chen F, Stein A, et al. 2009. Use of a multi-temporal grid method to analyze changes in glacier coverage in the Tibetan Plateau. Progress in Natural Science, 19(7): 861-872. |
| [44] | Yi S, Sun W K. 2014. Evaluation of glacier changes in high-mountain Asia based on 10 year GRACE RL05 models. J. Geophys. Res., 119(3): 2504-2517. |
| [45] | Zhang G Q, Yao T D, Xie H J, et al. 2013. Increased mass over the Tibetan Plateau: From lakes or glaciers? Geophys. Res. Lett., 40(10): 2125-2130. |
| [46] | Zhang Y, Cheng S Y, Zhao B K, et al. 2013. The feature of tectonics in the Tibet Plateau from new regional gravity signals. Chinese J. Geophys. (in Chinese), 56(4): 1369-1380, doi: 10.6038/cjg20130431. |
| [47] | 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, L17311, doi:10.1029/2009GL039459. |
| [48] | 丁志峰, 何正勤, 吴建平等. 2001.青藏高原地震波三维速度结构的研究. 中国地震, 17(2):202-209. |
| [49] | 刘遷, 迟效国, 李才等. 2001.青藏高原的形成和隆升机制综述. 世界地质, 20(2):105-112. |
| [50] | 罗栋梁, 金会军, 林琳等. 2012.青海高原中、东部多年冻土及寒区环境退化. 冰川冻土, 34(3):538-546. |
| [51] | 马荣华, 杨桂山, 段洪涛等. 2011.中国湖泊的数量、面积与空间分布. 中国科学:地球科学, 41(3):394-401. |
| [52] | 苏伟, 彭艳菊, 郑月军等. 2002.青藏高原及其邻区地壳上地幔S波速度结构. 地球学报, 23(3):193-200. |
| [53] | 汪汉胜, WuP, 许厚泽. 2009.冰川均衡调整(GIA)的研究. 地球物理学进展, 24(6):1958-1967, doi:10.3969/j.issn.1004-2903.2009.06.005. |
| [54] | 王林松, 陈超, 邹蓉等. 2014.利用GPS与GRACE监测陆地水负荷导致的季节性水平形变:以喜马拉雅山地区为例. 地球物理学报, 57(6):1792-1804, doi:10. 6038/cjg20140611. . |
| [55] | 王谦身, 安玉林, 张赤军等. 2003.地震学. 北京:地震出版社, 163-165. |
| [56] | 吴庆举, 曾融生. 1998.用宽频带远震接收函数研究青藏高原的地壳结构. 地球物理学报, 41(5):669-679. |
| [57] | 邢乐林, 孙文科, 李辉等. 2011.用拉萨点大地测量资料检测青藏高原地壳的增厚. 测绘学报, 40(1):41-44. |
| [58] | 张燕, 程顺有, 赵炳坤等. 2013.青藏高原构造结构特点:新重力异常成果的启示. 地球物理学报, 56(4):1369-1380, doi:10. 6038/cjg20130431 |