用cGPS研究青藏高原南缘现今垂向变动

方智伟， 邹蓉， 李志才， 王敏， 谭凯， 杨少敏， 王琪

1 中国地震局地震研究所， 地震大地测量重点实验室， 武汉 430071 2 中国地质大学地球物理与空间信息学院， 武汉 430074 3 国家基础地理信息中心， 北京 100044 4 中国地震局地质研究所， 北京 100029 5 武汉市应急管理局， 武汉 430010

0 引言

青藏高原隆升是新生代大陆活动构造最引人注目的地质事件(Molnar and Tapponnier，1975)，有关此次陆内造山的起始、过程、机制及力源一直是大陆动力学研究的热点问题(Tapponnier et al.，2001).青藏高原隆升研究涉及对其历史高程的推演及不同阶段隆升速率的厘定，目前前者主要借助高原内第三纪地层湖相沉积物的氧同位素分析(Rowley and Currie，2006)，后者包括从中新世以来出露变质岩的低温年代记录推算中下地壳折返率(Herman et al.，2010)，或利用山前第四纪冲积盆地测定抬升速率(Lavé and Avouac，2000)，以及利用大地测量(精密水准、GPS、InSAR)测定现今的地表垂直运动速率等(张青松等，1991； Jackson and Bilham，1994；Xu et al.，2000；张全德等，2001；姜卫平等，2008；Grandin et al.，2012)，观测手段虽多，但迄今获取的可靠资料却严重不足，限制了青藏高原构造演化研究的发展.

与地质学手段相比，用大地测量监测垂向变形具有量化清晰、物理意义明确、受客观因素制约较少的特点，对判断青藏高原隆升状态、动力作用具有指示意义.近期的连续GPS观测揭示了青藏部分地区的三维形变场，初显用大地测量方法解决大陆构造演化中动力学问题的良好发展势头(Sun et al.，2009；Fu and Freymueller，2012；Ader et al.，2012；Liang et al.，2013).尽管如此，现有观测结果仍有较大的不一致，例如早期跨青藏高原精密水准网复测资料显示其最大隆升位于藏南一带，隆升幅度达到10 mm·a-1(张青松等，1991)，这与尼泊尔境内水准资料推算的结果并不一致(Jackson and Bilham，1994).早期利用4年的GPS资料推测的高原中部唐古拉山一带隆升速率高达16 mm·a-1(Xu et al., 2000)，而用中国大陆构造环境监测网络(陆态网络，下同)GPS区域站10年资料解算的高原中部的构造隆升不过2 mm·a-1(Liang et al.，2013)，差异十分明显.

相对于精密水准和InSAR而言，GPS观测具有基准统一、精度均匀、观测简便，受地形、地貌限制较小等优势，但GPS站点的高程(24 h)观测精度不高于10 mm，而大陆内部地壳升降速率低于10 mm·a-1，垂向位移监测的信噪比低，定期、流动的GPS监测往往难以胜任.青藏地区连续GPS站点的增加，使得短时间内高精度测定垂直形变场逐渐成为可能(Sun et al.，2009；Ader et al.，2012；Hao et al.，2016).本文基于青藏高原南缘GPS连续站1995—2020年间观测数据，系统解算测站坐标时间序列(时序，下同)，求取站点垂直运动速率，合成具有信噪比、点密度较高的区域垂向位移场，为分析区域构造变形模式提供观测约束.

1 资料

本文研究汇集了印度次大陆、喜马拉雅及青藏南部地区固定观测98个GPS站点资料(图1)，这些站点大多设立在稳定岩石上，采用扼径圈天线，观测持续、质量稳定，可在相对较短的(5年)时间内，以优于1 mm·a-1的精度监测垂向位移(Lau et al.，2020).而区域内不定期观测的流动站的数据质量参差不齐，且观测10年左右的陆态网络区域站，其垂向位移观测精度一般在1～3 mm·a-1(Liang et al.，2013)，本文不予采用.

1.1 cGPS观测

从1997年以来，美国科罗拉多大学、加州理工学院与尼泊尔有关机构合作，布设基本覆盖尼泊尔全境的GPS连续台网(表1)，该网观测情况复杂，部分站点已不复存在或不再观测，有些则经历了2004年苏门答腊地震和2015年尼泊尔地震，维持至今.加州理工学院分析了该网早期站点的三维运动状况(Ader et al.，2012)，国内也有相关研究(Liang et al.，2013；Pan et al.，2018).本文系统处理该网从早期建站到今年春季数据，包括美方在加德满都等地为监测尼泊尔地震震后形变增补的新站(Zhao et al.，2017).

图1 cGPS站点分布及观测时长Fig.1 Distribution and observation time span of cGPS stations

表1 本文所用cGPS站观测情况Table 1 Observations of continuous GPS stations used in this paper

续表1

本文将周边国际GNSS服务组织(International GNSS Services，IGS)站点观测数据纳入处理，包括位于印度北方邦的LCK1-LCK4、印度中部的HYDE和IISC、不丹王国的TIMP、RBIT及位于拉萨的LHAZ(与陆态网络站LHAS并址)，这些测站中最长持续观测了25年.陆态网络12个连续观测基准站分布于青藏高原南部地区，此前一些研究多集中这些站点水平运动状况(Liang et al.，2013；Wang and Shen，2020)，与此前研究比，本文增加了最近几年观测资料.为增加西藏境内位移场的空间分辨率，本文还收集中国气象局于藏南布设的3个准连续站观测资料(表1)，这些站点观测时长约5年，尚未用于垂向位移的研究(Wang and Shen，2020).此外，本文引用了中国科学院青藏高原研究所于藏南布设的4个准连续站2007—2012年观测的处理成果(Liang et al.，2013)，以及印度国家地球物理研究所在喜马拉雅西段22个连续GPS站2009—2015年观测的处理成果(Yadav，2017；Yadav et al.，2019).

1.2 坐标解算

本文大部分GPS测站的原始数据处理基于喷气推进实验室(Jet Propulsion Laboratory，JPL)研制的GIPSY/OASIS Ⅱ-5.0软件.GIPSY处理普遍采用精密单点定位(Precise Point Positioning, PPP)模式(Zumberge et al.，1997)，获得ITRF2014框架下各站点的单日(UTC：0～24 h)时段坐标(Altamimi et al.，2016).本文具体处理策略包括：载波相位和伪距数据来自RINEX格式原始观测文件(采样率30 s)，L1和L2双频载波相位经过基于精码伪距辅助的周跳修复后，转换为无电离层折射L3组合相位(其截止高度角为10°)，并用IGS绝对天线相位中心模型修正测站天线中心相位偏差.对L3相位观测值的处理采用卡尔曼滤波/平滑算法，同步估算测站坐标、卫星单差的组合相位模糊度，以及作为随机游走变量的测站天顶方向大气延迟湿分量.受软件版本限制，没有借用JPL提供的卫星初始相位改正文件以辅助固定卫星单差组合相位模糊度.卫星轨道固定为JPL发布的无基准事后精密星历，同时采用JPL提供的卫星钟差文件.解算中各类物理模型分别是固体潮和极潮改正的IERS2010模型，以地球质心为基准的TPXO7.0(CM)海潮模型及对流层折射改正的全球经验大气+斜向延迟映射函数(GMF/GPT)组合模型.最后利用JPL发布的单日坐标转换参数，将无基准的单日坐标统一转换至ITRF2014参考框架下，形成各个站点可供变形分析的时序(图2).本文气象局测站(也包括喜马拉雅西段的cGPS站)采用GAMIT软件处理原始数据，基于双差分定位模式解算出包含这些测站及周边测站的单日区域网解，然后用GLOBK软件将区域网与全球IGS网联网平差，求取这些测站在ITRF2008下坐标时序，相关信息详见有关文献(Zhao et al.，2017；Wang and Shen，2020).

图2 GPS单日解垂向分量的时间序列图中虚线为2004年和2015年发震时刻，绿色代表经历2015年尼泊尔地震的台站.浅蓝色代表震前台站，深蓝色为同一站点的GAMIT时序，红色为震后布设的台站.Fig.2 Time series of the vertical component of GPS daily solutionsThe dotted lines in the figure show the time of earthquakes in 2004 and 2015, and the green represent stations that cross the 2015 Nepal earthquake. The light blue represents the stations deployed before the earthquake, the dark blue represents the time series of the same stations obtained by GAMIT, and the red represents the stations deployed after the earthquake.

1.3 垂直坐标时序分析

本文假定GPS垂向时序中包含三种可分离的信号成分：1)与构造运动或与晚更新世冰后回弹有关的线性位移信号；2)与测站周边环境有关(如大气、水文负荷的四季轮替)，可简化为具有一整年、半年周期变动特征的三角函数信号；3)大地震引起的同震位移和持续震后位移，前者可用单边函数表示，后者以指数或对数函数表示.综合以上，本文垂向时序展开为以下方程式：

(1)

式中：h(t)为测站大地高观测值，t为观测历元，t0为首次观测时刻，ci为季节性波动幅度，φi为季节峰值相位，H(t)为单边函数，dj为同震位移，tj为发震时刻，ej为震后变形幅度，τj为震后变形的半衰期.本文以式(1)拟合实测数据，解算时依据每个数据的解算精度定权，建立相应的误差模型，以加权最小二乘法估算各运动参数，准确反映测站垂直速率的信噪比.本文只关注2004年苏门答腊MW9.2、2015年尼泊尔MW7.8两次大震对时序的扰动，忽略其他地震的影响.另外，本文季节性波动项只考虑周年和半年周期项.

1.4 误差模型

GIPSY给出的单日大地高程的解算精度(名义误差)约为3～4 mm，主要反映了观测当日载波相位中的随机噪声.其他因素，如星历、电离层、对流层折射、测站周边温度变化、多路径等因素引起的各种偏差不在其中，这些偏差对于单日坐标而言属于系统性的，相邻数日或数周内系统偏差变化不大，对站点坐标的扰动一致，在时序中表现为有色噪声(不同时间观测值的统计性质相关).如简单将实测数据误差视为白噪声，虽不影响对各个信号的提取，但大大高估其解算精度，故本文时序分析中的误差模型采用白噪声+有色噪声的组合方式(Bos et al.，2013)，并结合赤池信息准则和贝叶斯准则选择其中最合适的模型.

借鉴此前尼泊尔地区GPS数据分析现状(Ader et al.，2012)，本文的有色噪声模型为广义高斯-马尔科夫模型，白噪声、幂率噪声、闪烁噪声和随机游走噪声模型都是其特例(Bos et al.，2013)，因此理论上广义高斯-马尔科夫模型能更好反映有色噪声的实际状况.每个测站具体的有色噪声模型主要依据其拟合残差的功率谱指数而定，相关计算采用Hector软件，同时给出站点垂向速率、周期项参数及其不确定性(表2).

2 垂直运动速率及精度评估

表2和图2所展示各类参数及信息可简单归纳如下：1)经GIPSY/OASIS Ⅱ-5.0处理的各测站时序的标准差(RMS)在6～10 mm之间，时序拟合残差限于[-20 mm, +20 mm]条带内.内华达大学基于高版本的GIPSY X-1.0处理的同为ITRF2014框架下的模糊度固定解时序(http：∥www.geodesy.unr.edu)的标准差为3～6 mm，比本文模糊度浮点解略好，但对比本文与内华达大学计算的IISC、HYDE、LHAZ、CHLM的长期速率, 互差小于0.2 mm·a-1；2)经GAMIT处理的ITRF2008框架下时序标准差在3 mm左右，虽然双差精密定位较单点精密定位能更好消除一些共模误差，能更好分辨出微弱变形信号，但在提取长期速率上，GIPSY与GAMIT时序的一致性保持在0.05 mm·a-1以内；3)区域内升降速率约束在[-2 mm·a-1, +6 mm·a-1]区间(相对ITRF框架原点——地球质心)，估算精度多半优于0.5 mm·a-1；4)各测站的季节性波动信号与地震位移信号得到较好的分离，残差时序没有明显的扭曲和跳变；5)除少数几个测站外，各站的周年项波动幅度一般在6～10 mm，估值误差在2 mm以内，周年项峰值位于5～6月的印度季风期，喜马拉雅山前测站波动(最大15 mm)明显大于高原内部测站；6)半年项波动幅度不超过周年项的一半，但估值误差类似；7)尼泊尔地震导致震区内测站(最大1.3 m)永久性垂向位移，震后变形也导致部分测站短期升降10～60 mm，此外苏门答腊地震导致尼泊尔及印度次大陆测站隆升大致10～15 mm.

表2 cGPS时序分析结果Table 2 Time series analysis results of continuous GPS stations

续表2

2.1 垂直速率的区域特征

GPS测站的垂直运动特征概括为：恒河平原及喜马拉雅山前(海拔500 m以下)总体下降不大于3 mm·a-1，DNGD沉降幅度最大，速率为3.1±0.4 mm·a-1；低喜马拉雅地区(海拔500～1500 m)处于升降转换区，高喜马拉雅地区(海拔1500～4500 m)相对印度次大陆最大隆升接近2～8 mm·a-1，在30 km距离内，抬升速率从1 mm·a-1(NAGA)急剧增加到6 mm·a-1(GUMB)，与地形陡变完全一致. 藏南地区(海拔4500 m以上)的GPS站在大约50 km距离内抬升速率从4～6 mm·a-1(GHER、CHLM)减少到1～2 mm·a-1(LMJM、SYBC)，雅鲁藏布缝合线以北(平均海拔4500 m)GPS站逐渐从上升约1 mm·a-1转变至下降1 mm·a-1，XZNM(尼玛)沉降最大为1.5±0.4 mm·a-1.

喜马拉雅1975—1990年间两期精密水准数据显示高喜马拉雅与低喜马拉雅相对前陆盆地的抬升分别约为4±1 mm·a-1和2.0±0.3 mm·a-1(Herman et al.，2010).尼泊尔西部的InSAR观测更凸显这种升降变化(图3)，高喜马拉雅相对于恒河平原最大隆升达到8 mm·a-1(Grandin et al.，2012).尽管本文直接解算的尼泊尔cGPS垂向速率场及印度西北地区cGPS观测结果(Yadav et al.，2019)的空间分辨率较低，但结合InSAR、精密水准，现有大地测量综合结果完整反映了大陆俯冲带前陆盆地沉降、山前快速抬升、造山带后缘隆升延展、幅度急剧降低这一典型特征.

图3 喜马拉雅地区垂向速率剖线左图：InSAR图像、水准线路和GPS站分布图，紫色和蓝色方框代表A、B剖线，色棒标示InSAR和水准点速率，黑色箭头为相对欧亚板块的GPS水平运动速率(Wang and Shen，2020)；右图：GPS(红或紫色方框)，InSAR(棕色圆点)和水准(蓝色圆点)速率及地形合成图.GPS速率未做任何修正.图中近场和远场GPS采用两种距离尺度表示，水准和InSAR误差棒代表1 mm·a-1的假设误差，GPS误差为真实误差.绿色虚线为震间运动模型的预测速率(Grandin et al.，2012)，棕色圆环代表此前的GPS结果(Fu and Freymueller，2012).Fig.3 Section lines of vertical velocities in HimalayaLeft: InSAR image, leveling line and distribution of GPS stations. Blue and purple boxes represent section lines A and B. Color bars indicate InSAR and leveling benchmark velocity. Black arrows represent GPS horizontal motion rates relative to the Eurasian plate (Wang and Shen, 2020). Right: Composite image of GPS (red or purple boxes), InSAR (brown dots) and leveling (blue dots) rate with terrain. No correction has been made to the GPS rate. In the figure, near-field GPS and far-field GPS are represented by two distance scales. Leveling and InSAR error bars represent hypothetical errors with 1 mm·a-1, while the GPS error is real error. The green dotted line represents the predicted velocity of the interseismic motion model (Grandin et al., 2012), and the brown rings represent the previous GPS results (Fu and Freymueller, 2012).

2.2 内符合精度分析

此前尼泊尔境内站点及藏南部分GPS站垂直速率曾有文献记载(Fu and Freymueller，2012；Ader et al.，2012；Liang et al.，2013；Pan et al.，2018)，本文将处理结果与此前结果逐一对比(图4).加州理工最早公布测站速率误差一般为1～3 mm·a-1(Ader et al.，2012)，比本文大3倍.加州理工除顾及有色噪声外，还考虑了时序不连续可能引入的阶跃误差，因此总体上精度估计偏保守.其他三组与本文采用完全一致的误差评估算法(Fu and Freymueller，2012；Liang et al.，2013；Pan et al.，2018)，图4显示本文速率精度总体上最优，这得益于某些测站已有了更长时序，当时间跨度大于5年时(XZAR、XZNM、SNDL)，PPP给出的垂向速率精度在0.5～0.8 mm·a-1，大于10年时(KKN4、KAWA、CHLM)，误差进一步降低到0.3 mm·a-1，而超过20年观测(LHAS、HYDE、IISC)，误差甚至可低至0.1 mm·a-1.与此对照，陆态网络观测时长12年的区域站精度一般低于1 mm·a-1(Liang et al.，2013)，如此估计25年观测历史的流动站有望将误差降低到0.5 mm·a-1以内.

图4 4组GPS速率结果比对与此前公布的4组速率结果对比，西藏地区陆态站点用绿色点位表示.虚线阴影区代表互差在+/-1 mm·a-1以内区域.Fig.4 Comparison of GPS velocity results from four groupsComparison with the previously published four groups of velocity results, the CMONOC sites in Tibet are represented by green dots. The dotted shaded areas represent the areas with errors within +/-1 mm·a-1.

需要解释的是：本文中少数站点尽管具有更长时序，但有些站点的速率精度仍低于此前结果(Pan et al.，2018).由于利用了尼泊尔震后数据，为抵消同震和震后变形影响，本文引入了附加参数，相应降低了长期速率的估计精度.此外，GHER、LMJG震后位移量较小，从时序中分离震后、震间信号需要借助一定先验假定，比如，震后位移与同震位移方向一致，震后位移幅度不超过同震位移等，因此速率计算与先验信息有关，有一定的不确定性.即便如此，本文与此前结果没有系统性偏差，例如与最近一次的结果(Pan et al.，2018)比，速率互差在1 mm·a-1内，仅PYUT、DRCL、SNDL三站互差超过2 mm·a-1.

2.3 模型对比及外符合精度

喜马拉雅变形模型以往主要基于GPS水平位移、精密水准、InSAR三种资料(Bilham et al.，1997；Grandin et al.，2012)，GPS垂直位移数据因精度较差，对建模贡献甚小(Fu and Freymueller，2012).为客观评估本文结果，将GPS测站分东西二组(图3)，沿N17°E方向构建两条跨喜马拉雅、宽度近200 km的带状速率剖线(剖线A、B)，剖线A内GPS速率与C波段ENVISAT卫星InSAR视线方向速率对比(Grandin et al.，2012)，剖线B则是与尼泊尔一等精密水准复测结果对比(Jackson and Bilham，1994).

跨喜马拉雅的InSAR干涉图像覆盖83°～84°经度带，从低喜马拉雅延伸到雅鲁藏布，南北长250 km(图3)，是处理2003—2010年间29幅降轨原始图像，将不同时段、独立干涉图叠加而成(Grandin et al.，2012).变形起算点位于图像最北端，变形信号从北到南累积，南北两端像素点速率差接近10 mm·a-1，InSAR条带内各像素点变形精度估计为0.7～3.1 mm·a-1(Grandin et al.，2012)，变形分辨率为每公里0.2 mm·a-1.尽管理论上InSAR图像代表地表变形在卫星视线方向(LOS)投影，但实际上更多地反映了垂向变动的分布，就A剖线而言，LOS 速率与垂向速率的偏差不超过0.5 mm·a-1.本文将InSAR最北端像素点的速率设为JRGR速率值1.2 mm·a-1，推算整条剖线其他像素点在ITRF框架上的速率.

尼泊尔境内的精密水准线路南起低喜马拉雅斯瓦里克山前(Jackson and Bilham，1994)，向北经加德满都终于高喜马拉雅山前的西藏樟木口岸，总长350 km，沿线高程起伏不超过1500 m, 该线路1977—1991年间观测，观测误差为1.1 mm·km-1，按误差传播率估计，水准线路南、北两端的速率差相对精度为1.6～2.8 mm·a-1(Jackson and Bilham，1994).珠峰北坡一带的水准线路与尼泊尔独立，起算点位于定日(陈俊勇等，2001).本文将尼泊尔水准起始点速率等同于SIMP的速率值(-1.0±0.5 mm·a-1)，而将珠峰北水准线路的起算点设为最靠近XZZF(1.6±0.5 mm·a-1)的一个测点，相应推算各自路线上其他水准点的ITRF框架速率值.

对比GPS、精密水准、InSAR资料可见(图3)，三者在1 mm·a-1的精度范围内保持了一致，剖线B内只有DAMA、RMJT、SYBC偏差较大(1.5～2.0 mm·a-1)，如果将垂向变形的模型值(Grandin et al., 2012)与GPS实测值对比，则DAMA、TPLJ、SYBC偏差大于1.5 mm·a-1.偏差较大站点都远离水准线路，可能叠加了局部变形效应，与水准点揭示的变形特征有所不同.

剖线A中InSAR与GPS的一致性更佳，只有KLDN与InSAR偏差大于4 mm·a-1.原因可能是山前一带InSAR观测误差较大或KLDN叠加了前沿断裂变形的影响.与模型值(Grandin et al., 2012)对比，只有KLDN、YARE、XZZB有近2 mm·a-1的偏差，其中YARE、XZZB与所在区域的InSAR观测值都系统性偏离了模型值，且GPS、InSAR的变化趋势也相同，这种异常变形可能与2006—2008年仲巴序列地震有关.整体而言，本文比早期结果(Fu and Freymueller，2012)更接近模型值.

2.4 垂直速率修正

GPS测定的垂向变动主要包含三种运动成分：1)板块挤压引起的构造变形；2)地表水变化导致的升降，近年来喜马拉雅冰川加速融化和恒河平原地下水抽取加剧(Rodell et al.，2009；Yi and Sun，2014)，将导致区域地壳大幅调整；3)末次冰期以来北半球岩石圈松弛变形及青藏高原冰后回弹.依据近年来GPS和水准观测，印度—欧亚碰撞带内板块水平移动30～40 mm·a-1(Wang et al.，2001)，升降不超过6 mm·a-1(Jackson and Bilham，1994；王庆良等，2008；Ader et al.，2012；Liang et al.，2013)，其垂向运动受环境影响更大.

地表水变化引起的均衡调整一般根据时变重力场推算的区域地表载荷变化及弹性地球变形响应来估算(Fu and Freymueller，2012).本文利用最新的GRACE观测结果，假定时变重力场完全是地表水变化决定，忽略其他因素的贡献(Rao and Sun，2021)，估算喜马拉雅及藏南地区均衡隆升大致在0.2～0.4 mm·a-1范围，比此前估算0.4～0.8 mm·a-1低一半左右(Pan et al.，2018).两组结果皆有0.05 mm·a-1的标称精度，大大低于两组间的系统性偏差，说明用不同滤波算法平滑GRACE资料，给出的修正值存在至少0.2～0.3 mm·a-1的不确定性.

末次冰盛期导致的青藏高原隆升同样存在不确定问题.按极端区域性冰川模型推测，青藏高原中部地壳冰后回弹的升幅可达4～7 mm·a-1(Kuhle，1998；Kaufmann，2005)，但另有研究表明青藏高原冰盖规模不大、影响极小(Derbyshire，1991；Matsuo and Heki，2010；Liang et al.，2013)，本文也不加考虑.但北半球岩石圈松弛变形不容忽视，据现有大地测量约束下的全球冰川模型，亚洲大陆的冰后回弹为0.3～0.5 mm·a-1(Peltier，2004；Peltier et al.，2015)，为此本文统一用0.4 mm·a-1作为修正值.需要说明的是，现有修正值的误差难以估算，本文取修正值的一半.

与以往研究不同(Liang et al.，2013；Hao et al.，2016；Pan et al.，2018)，本文还考虑地心运动对GPS速率的影响.根据全球大地测量结果(Argus，2012；Métivier et al.，2020)，末次冰期以来北半球岩石圈松弛变形以及内部质量迁移导致地球形心相对其质心以0.5～1.0 mm·a-1速率向北极运移，由于ITRF框架的原点为地球质心，GPS实测速率相对于质心，而构造运动相对于地球形心，地心的视运动将导致中低纬度地区真实的垂向运动偏大0.4 mm·a-1左右(Ding et al.，2019)，必须加以修正.本文按0.8 mm·a-1地心速率乘测站纬度正弦值计算修正值，修正值的精度大致在0.2 mm·a-1左右.综合以上因素，现有GPS速率中非构造影响在1 mm·a-1左右，按保守估计，修正值的不确定度大致在0.3～0.4 mm·a-1.

3 讨论

本文研究涉及印度—欧亚板块碰撞三个不同的构造区(Molnar and Tapponnier，1975；Tapponnier et al.，2001)，分别为：刚性的印度板块，整体性向北挤压；位于板块边界的喜马拉雅，造山楔强烈挤压变形，向南仰冲、推覆在印度板块之上；青藏高原内部，处于挤压、拉张状态的拉萨、羌塘地块.印度板块的地壳厚度不到40 km，藏南(喜马拉雅北坡至雅鲁藏布)地区厚度接近75 km，增厚近一倍,高原中部地壳厚度减薄到65 km(Nábělek et al.，2009)，地壳从增厚到减薄转换对应了GPS垂向运动从下沉到隆升、再到下沉的空间展布.

3.1 印度次大陆

沿A剖线印度板块相对欧亚板块的水平速率从北到南为36～38 mm·a-1(Wang and Shen，2020)，如扣除非构造因素，现有GPS表明印度板块前缘部分已下沉2～3 mm·a-1，使其以3°～4°倾角向下插入欧亚板块.根据尼泊尔地震的同震破裂产状，位于低喜马拉雅底部的印度板块以7°倾角向北运移(Wang and Fialko，2015)，远震接收函数成像显示大约在高喜马拉雅到雅鲁藏布之间，印度板块的倾角陡变到11°～12°(Nábělek et al.，2009)，板块边界面的倾角每100 km就减小3°～4°(板块曲率(0.5～0.6)×10-6m)，为海洋俯冲带倾角的1/10～1/3(Bletery et al.，2016).考虑到恒河平原新生代沉积层最大深度为5～6 km(Decelles et al.，1998)，大陆俯冲的前陆伸展盆地至少100 km宽，从NEPJ延伸到BHUP.现有GPS结果还表明，德干高原北部为约500 km宽的前缘隆起区(从BHUP到JBPR)，德干高原南部(IISC、HYDE)几乎没有构造性垂向运动，GPS观测与造山理论的预测一致(Molnar and Lyon-Caen，1988)，展示了印度板块在青藏高原的重压下，挠曲由浅入深的几何结构(图5).

图5 青藏高原南缘垂向运动速度场GPS 实测速率经地心(坐标框架原点视)运动、冰川均衡调整模拟和地表水变化GRACE观测三种改正，改正后速率见表1.Fig.5 GPS-derived vertical velocity field in southern Tibetan PlateauThe GPS-derived rates have been corrected by geocenter motion (frame origin apparent motion), glacial isostatic adjustment simulation and GRACE-derived groundwater variation. The corrected rates are shown in Table 1.

3.2 喜马拉雅

如果GPS观测到的印度次大陆前缘部分沉降可视为永久性、不可恢复变形，但GPS观测的喜马拉雅变形大部分是弹性的、与板块边界错动及大震周期有关的变形.作为板块边界的喜马拉雅主冲断层最南缘的100 km段在两次大震间完全闭锁，低喜马拉雅地区的震间变形缓慢，高喜马拉雅地区快速隆升、变形累积(Bilham et al.，1997).这种隆升不是单调持续的，特大地震将导致大部分隆升回落，例如2015年尼泊尔地震时海拔3000 m以下的地带全面隆升，而海拔3000 m以上地区整体下降(Elliott et al.，2016)，凸显出地震在弹性变形从高海拔向低海拔区域的迁移过程中的关键作用.

GPS反映的弹性变形状态与瞬时地震观测的结果完全相反，加德满都KKN4沉降速率为0.2 mm·a-1，尼泊尔地震中却抬升1.3 m，中尼边界地带高海拔的CHLM隆升速率2.6 mm·a-1，地震时下降0.6 m，弹性变形集中在地形陡变带，永久性变形集中在低海拔的山前地带，喜马拉雅地区现今变形表现为双重性.喜马拉雅的大震周期平均500～1000年(Feldl and Bilham，2006)，2015年尼泊尔地震不能一次性将大震间积累的全部变形从高喜马拉雅转移至斯瓦里克，只有1934年比哈尔MW8.2地震乃至更大地震才有可能错断整个主冲断裂，至最南端主前缘断裂(Sapkota et al.，2013)，将震间绝大部分变形迁移至前缘地带，转换为永久性地形抬升.但高喜马拉雅地带GPS隆升速率变化与地形变化高度一致，表明这种转换并不彻底，仍有部分弹性变形永久存留在高海拔地带，维持其陡变地形(Meade，2010).

3.3 青藏高原

最早的跨青藏高原的精密水准显示，藏南(喜马拉雅北坡至雅鲁藏布)地区近10 mm·a-1隆升(张青松等，1991)，这种大幅度现代隆升得到早期流动GPS观测的佐证(Jackson and Bilham，1994).但本文结果，该地区GPS站仅有0～2 mm·a-1的上升幅度，显然与10～16 mm·a-1极快速隆升不符(Xu et al.，2000).实际上，初期(1960)跨青藏高原水准观测误差偏大，后两期(1980、1990)复测精度较好，其显示拉萨以南地区隆升约为2 mm·a-1(张全德等，2001)，与现有藏南一带的GPS速率十分接近.

藏南地区隆升固然来自喜马拉雅主冲断裂的活动，但仍有一部分与印度板块上地壳撕裂、增生到欧亚板块底部有关，GPS观测还难以区分这两种机制的具体贡献.此外沿雅鲁藏布缝合线XZAR、XZZB、YARE表现为较大幅度隆升(>2 mm·a-1)，难以完全依据板块边界的构造活动来解释，也许是雅鲁藏布缝合线附近的新构造活动所致(例如2006—2008年仲巴地震)，现有测站密度还远不足以证实.

雅鲁藏布缝合线以北的GPS速率扣除非构造因素后，基本表现为下沉1～3 mm·a-1，雅鲁藏布缝合线南北两侧由南升转为北降.现有结果一方面与青藏高原末次冰期巨厚冰盖消融导致的隆升分布特征不同，可以排除区域性大冰盖的可能性；另一方面，即使忽略北极冰盖消融的影响，雅鲁藏布缝合线以北地区也有最大2 mm·a-1下降，与拉萨地块(雅鲁藏布缝合线至班公—怒江缝合线)的南北向缩短幅度相当(Wang and Shen，2020)，依据印度板块插入青藏高原下方的倾角在雅鲁藏布缝合带逐渐变小不足以解释这种幅度的下降，可能是拉萨、羌塘地块东西向拉张变形相应导致地壳减薄，其幅度大大超过挤压，导致高原内部塌陷，相对于西部站点，88°E以东区域GPS测站普遍下沉.

4 结论

青藏高原及喜马拉雅地区连续GPS观测以优于1 mm·a-1的精度约束区域垂直位移场.GPS与精密水准和InSAR观测对比验证，比较完整地揭示了印度—欧亚碰撞带垂直形变特征：从印度次大陆到青藏高原中部，地表经历沉降-隆升-沉降，对应了地壳从增厚到减薄的转换过程.除雅鲁藏布缝合线以南的区域继续保持隆升状态，并向外扩展外，鉴于青藏高原内部GPS测站普遍下降1～2 mm·a-1，青藏高原的中南部可能已不再隆升.

致谢本文原始GPS数据来自中国大陆构造环境监测网络，中国气象局、加州理工学院，美国鲍灵-格林州立大学富宇宁提供了GRACE修正值.审稿人与编辑对本文提出了宝贵的修改意见和建议.在此一并表示感谢.