青藏高原南缘作为印度板块与欧亚板块碰撞的前缘地带，其地壳形变问题一直是地学界的研究热点，并且取得了显著成果^[1-4]。传统的大地测量手段(如全站仪、水准仪、GPS等)观测覆盖范围较小、人力成本高、测站空间分布不均匀，从而使得中长空间尺度重力变化估计的不确定性较大。GRACE卫星的发射开创了高精度全球重力场观测和气候变化实验的新纪元，被认为是当前中长空间尺度高精度地球重力场恢复、地球表层陆地水和海水物质质量空间分布和变化探测研究中极具发展潜力的方法之一^[5-7]。

张良镜等^[8]基于GPS和卫星重力观测对地表形变的季节性变化进行研究，结果表明两者周年变化的相位和振幅在大部分地区具有较好的一致性；Steffen等^[9]基于GRACE时变资料和地面绝对重力数据对Fennoscandia地区的隆升速率进行分析，认为绝对重力结果能较好地验证GRACE结果。Fu等^[10]和王林松等^[11]分别在喜马拉雅地区利用GPS与GRACE技术监测地壳垂直运动情况和陆地水负荷导致的季节性水平形变，结果显示两者垂直分量的季节性变化具有较好的一致性。研究表明，卫星重力结合GPS、绝对重力等传统大地观测手段可有效获取大范围的地表质量负荷迁移及地壳形变信息。本文结合GRACE与GPS以及地面绝对重力技术对青藏高原南部地区地壳形变进行探讨，以期为该地区重力场时空演化特征提供相关参考。

1 数据和方法 1.1 GRACE数据处理

本文同时采用CSR、GFZ和JPL三家机构提供的Level(RL05)数据，并将平均结果作为最终的估计。GRACE卫星对低阶项不敏感，本文采用Swenson等^[12]提供的一阶项和卫星激光测距(SLR)得到的二阶项进行替换。对2003~2013年所有月份的球谐系数求取平均值，然后将每个月的球谐系数与该平均值作差得到异常场。基于P5M15策略对高于15阶的奇偶项分别采用5阶多项式进行拟合并予以扣除来削弱南北条带误差的影响；采用200 km扇形滤波^[13]削弱高阶球谐系数的噪声影响。球谐系数可转换为分辨率为1°×1°全球分布的重力异常、等效水柱高或垂向位移^[14]：

$ \begin{align} & \ \ \ \ \ \ \ \ \text{d}{{r}_{h}}(\theta , \lambda )=R\sum\limits_{n=1}^{{{n}_{\max }}}{\frac{{{{{h}{'}}}_{n}}}{1+{{{{k}{'}}}_{n}}}\sum\limits_{m=0}^{n}{{}}} \\ & [\Delta {{C}_{nm}}\cos (m\lambda )+\Delta {{S}_{nm}}\sin (m\lambda )]{{{\tilde{P}}}_{nm}}(\cos \theta ) \\ \end{align} $

(1)

式中，dr_h(θ, λ) 为地球表面垂向位移变化，h′_n和k′_n为n阶负荷勒夫数，R为地球平均半径，${{{\tilde{P}}}_{nm}}$表示阶数为n、次数为m的归一化缔合勒让德函数，ΔC_nm和ΔS_nm为Stokes系数变化值。勒夫数前3阶分别取：h₁=-0.286, l₁=0.104, k₁=0.027；h₂=-0.991, l₂=0.233, k₂=-0.305；h₃=-0.105, l₃=0.699, k₃=-0.196。为了与GPS观测的整体质量负荷影响进行比较，在GSM资料基础上需要加上反映大气和海洋质量变化引起时变重力场变化的GAC产品系数。

1.2 GPS数据

GPS数据由南加州综合网络(southern California integrated GPS network, SCIGN)构建的远程连续GPS参考站(remote continuous GPS station package，RCGSP)提供，数据处理基于ITRF2005框架，由GAMIT/GLOBK软件分析得到单日解时间序列，本文共选取LHAZ(91.1°E, 29.7°N)、TPLJ(87.7°E, 27.4°N)、CHLM(85.3°E, 28.2°N)、JMSM(83.7°E, 28.8°N) 四个测站参与研究。

1.3 绝对重力数据

FG5绝对重力仪的标称精度优于5 μGal，不同仪器之间的较差约为1~2 μGal，无明显系统误差。假定地球内部无质量迁移，1 cm地面垂直运动引起的重力变化约为3 μGal^[15]。为得到高精度的地面绝对重力值，对各期观测数据进行处理^[16]，主要包括大气压力、极移、固体潮、海潮等改正。本文利用LHAZ测站1992~2014年的绝对重力资料进行分析比较。

2 结果与分析 2.1 GRACE和GPS结果比较

图 1为GPS得到的LHAZ测站长期性时间变化序列，从图中可以看出，LHAZ测站1999~2015年除具有明显的周期性波动变化外，还存在明显的线性变化趋势，主要反映出相对稳定的垂向构造运动。周期性变化振幅约为15 mm，与季节变化密切相关。对于线性变化趋势，LHAZ测站垂向速率约为1.5 mm/a，显示出LHAZ测站持续隆升的特点。

图 2为由GRACE和GPS得到的2003~2013年垂向位移长期性时间变化序列，其中部分GPS资料因数据采集等因素造成缺失。从图 2可以看出，4个测站由GPS和GRACE得到的垂向位移时间序列变化中均具有明显的周年变化特征，与季节性变化密切相关，因此认为由相同的物理因素(同源)所激发。部分学者^{[10, 17-18]}从加权均方根WRMS角度评价两个时间序列之间的一致性，结果表明，大多数台站在扣除GRACE时间序列后，GPS垂直分量明显减小，这在某种程度上表明季节性陆地水负荷是造成形变垂直分量的主要因素。

对于GPS数据，选择1σ标准偏差作为误差估计，通过最小二乘拟合得到CHLM、JMSM、TPLJ测站的上升速率分别为3.8±0.09 mm/a、2.6±0.14 mm/a和1.4±0.08 mm/a。其中CHLM和TPLJ测站的建站时间分别为2006年和2007年，JMSM测站部分年份数据丢失，3个测站均可以看出明显的地壳上升趋势。对4个测站的上升速率进行平均，得到研究区平均隆升速率为2.3±0.11 mm/a。

对于GRACE资料，考虑到其误差为测站所在经度加90°换算至太平洋上相应纬度地区的数值，因此图 3中海洋地区径向位移值(红色空心圆)可作为相应测站的误差估计。从GRACE得到的时间序列变化中可以看出明显的周年变化趋势，基于最小二乘原理得到4个测站的上升速率分别为：0.29±0.05 mm/a(LHAZ)、0.34±0.04 mm/a(CHLM)、0.37±0.05 mm/a(JMSM)、0.39±0.04 mm/a(TPLJ)，平均上升速率为0.35 mm/a。

地面测站的长期变化主要受大尺度板块运动、局部构造形变、负荷质量迁移及冰川均衡调整(glacial isostatic adjustment, GIA)等因素的综合影响，空间大地测量得到的测站速率除观测误差外，主要反映这些因素综合引起的测站运动。如果测站位于构造稳定地区，则测站垂直分量主要反映构造运动、负荷迁移引起的形变及GIA引起的垂直运动。GRACE结果为全球质量和地表弹性负荷变化导致固体地球重新产生负荷形变的综合影响，可反映一定尺度范围内总体位移的变化情况，本研究为200 km扇形滤波空间尺度范围的平滑结果。GPS结果可反映测站及周边地区的变化情况，其更容易受到局部因素影响而表现出高频变化特征；基于GRACE估算的负荷位移主要反映流体质量迁移的影响，其长期变化机理较为复杂。

2.2 绝对重力结果

高精度绝对重力测量(FG5)一方面可反映地壳上某点的垂直运动信息，即地壳内部质量的运动；另一方面可测量地表相对于地心的垂直变化，可作为一种附加和独立的与地心有关的垂直参考基准。而GPS、水准测量等常规大地测量技术得到的仅为纯几何高程变化，无法得到地壳内部物质的运动信息。张为民等^[19]通过FG5绝对重力仪监测青藏高原的隆起，对比分析1993年和1999年LHAZ测站的重力观测资料，结果表明，6 a内重力值减小12 μGal，证实LHAZ处于隆升阶段，并初步估计上升速率为10 mm/a。王勇等^[15]对滇西地区和西藏拉萨的绝对重力重复观测资料进行分析，结果显示LHAZ测站的绝对重力以-1.82±0.91 μGal/a的速率下降，从重力学角度反映出青藏高原的隆升，并通过俯冲模型计算重力变化，认为其主要原因为印度板块向欧亚大陆俯冲。邢乐林等^[20]利用绝对重力资料和GPS结果定量分析认为，青藏高原拉萨测站地下地壳具有增厚特征。Sun等^[21]对拉萨、昆明、大理3个测站的GPS与绝对重力数据进行分析，以大地测量技术揭示了青藏高原底部质量损失，显示出该地区地壳增厚的特点。Sun等^[22]加入卫星数据重新对拉萨、昆明、大理3个测站进行研究，得到3个测站的平均速率为-0.66±0.49 μGal/a，地壳增厚速率为1.9±1.4 cm/a。本文按照标准绝对重力数据的处理步骤^[19]，将重力观测值均归算至测站标石表面高度，得到拉萨测站绝对重力变化时间序列(图 4)。

图 4中纵轴已经从重力观测值中扣除一个均值，因此反映的是LHAZ测站1993~2014年重力值的相对变化情况。关于结果的误差分析具体见王勇等^[15]和张为民等^[19]的研究。从图 4中可以看出，除2010-10采用A10绝对重力仪外，其余年份均采用FG5仪器。何志堂等^[23]对A10与FG5绝对重力仪进行2次比对测量实验，结果显示两者之间无明显差异，一致性较好。

3 讨论与分析

前文通过拉萨测站(LHAZ)和尼泊尔地区RCGSP网络中3个测站(TPLJ、CHLM、JMSM)的GPS资料得到4个测站近10 a的平均上升速率为2.3±0.11 mm/a，显示青藏高原南部地区在持续隆升。通过由GRACE得到的垂向位移时间序列可知，2003~2013年除季节性变化外还具有轻微上升趋势，反映出负荷质量迁移对测站周边地区的综合影响。

基于LHAZ测站1993~2014年的绝对重力资料，在考虑大气压力、极移、固体潮、海潮、仪器高归算等因素后，通过最小二乘拟合得到LHAZ测站的重力变化速率约为-1.76 μGal/a，与王勇等^[15]得到的-1.82±0.9 μGal/a和邢乐林等^[20]得到的-1.85±0.79 μGal/a的结果较为一致，表明LHAZ测站绝对重力值在持续减小，且速率较稳定。分析表明，LHAZ测站1998-09~2000-09期间下降趋势最为显著，约为6 μGal/a；2004-08~2013-09下降趋势较为平缓，约为1 μGal/a。整体来说，LHAZ测站重力值在近20 a下降了35 μGal，反映出拉萨及周边地区地壳的垂直运动(隆升)和内部质量迁移(亏损)情况，本文认为其主要原因为印度板块与欧亚板块挤压碰撞进而引起青藏高原地壳形变，物质分布发生变化。王勇等^[15]将重力变化全部归因于地壳隆升，得到拉萨测站的隆升速率为8.7 mm/a。Sun等^[21]认为印度板块向欧亚板块俯冲不仅造成青藏高原隆升，同时也会造成地壳底部增厚，且为持续平稳的过程。强调地壳的双倍增厚是高原隆升的主导机制，也是目前比较合理和普遍接受的隆升机制^[24]，低密度的山根挤出高密度的地幔物质导致区域内物质亏损，进而造成重力减小。对于高原隆升和底部地壳增厚，根据物理大地测量基本理论，存在如下重力变化关系^[25]：

$ \begin{align} & 2\pi G({{\rho }_{m}}-{{\rho }_{c}}){{h}_{m}}={{\delta }_{g}}-3.086{{h}_{t}}+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ 2\pi G{{\rho }_{c}}{{h}_{t}}-{{\delta }_{{{g}_{d}}}}-{{\delta }_{\text{GIA}}} \\ \end{align} $

(2)

式中，ρ_m为地幔平均密度，取值3.3 g/cm³；ρ_c为地壳平均密度，取值2.78 g/cm³；δ_g为某段时间内的重力变化；δ_gd为侵蚀引起高原水土流失造成的重力变化；h_t为地表几何垂直位移变化量；h_m为地壳底部增厚量。

将式(2)除以Δt得到地壳厚度的变化率为：

$ \begin{align} & \overset{\centerdot }{\mathop{{{h}_{m}}}}\, \text{=}(\overset{\centerdot }{\mathop{{{\delta }_{g}}}}\, -3.086{{\overset{\centerdot }{\mathop{h}}\, }_{t}}+2\pi G{{\rho }_{c}}{{\overset{\centerdot }{\mathop{h}}\, }_{t}}-{{\overset{\centerdot }{\mathop{\delta }}\, }_{{{g}_{d}}}}-{{\overset{\centerdot }{\mathop{\delta }}\, }_{\text{GIA}}})/ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2\pi G({{\rho }_{m}}-{{\rho }_{c}}) \\ \end{align} $

(3)

Sun等^[22]研究认为，拉萨地区由于侵蚀引起水土流失造成的重力影响根据Bougour层改正约为-0.25±0.1 μGal/a，GIA改正约为-0.25±0.1 μGal/a。对于拉萨测站，由绝对重力得到的重力变化率为-1.76 μGal/a，由GPS结果得到的垂直位移增长率为1.5 mm/a，引起的重力变化为-0.29 μGal，将剩余未解释部分(-0.97 μGal/a)归因为地壳底部增厚，代入式(3)可得拉萨测站地壳底部增厚速率约为4.4 cm/a。由于仅有拉萨测站的绝对重力资料，其结果并不能代表整个青藏高原南部的地壳形变。相关研究显示，高原地壳并不是印度地壳与西藏地壳的简单叠加，地壳厚度存在约40 km的起伏，整个青藏高原地壳平均厚度约为70 km，拉萨块体与羌塘块体地壳厚度相对平坦，与高原地壳厚度接近。因此，拉萨测站的结果基本可反映青藏高原整体的地壳变化趋势，即地壳隆升与底部增厚^[20]。

4 结语

本文联合GRACE卫星重力、GPS和绝对重力资料对青藏高原南缘地壳形变进行分析，结果表明，GPS与GRACE得到的结果趋势较为一致，均具有明显的周年变化特征，并可从趋势项结果中看出青藏高原南部区域在持续隆升。结合拉萨测站绝对重力数据，考虑土壤剥蚀及GIA影响，估计该测站地壳增厚速率为4.4 cm/a，反映出青藏高原南部地壳形变特征——地壳隆升与底部增厚。因搜集到的绝对重力资料较少，今后可结合更多的重力数据分析青藏高原南缘地壳形变问题，以增强结果的可靠性。