基于块体模型的青藏高原及邻区地壳三维构造形变分析

党亚民，杨 强，王 伟，梁玉可

1. 中国测绘科学研究院，北京 100830； 2. 山东科技大学，山东 青岛 266590

青藏高原是世界屋脊，与世界上的其他高原不同，青藏高原地壳具有复杂、不均一、多地体拼贴、巨厚(60～70 km)等特点[1]，是地学研究的宝库。青藏高原地壳水平方向变形主要呈现东西向拉伸、南北向压缩趋势，自南向北运动速度逐渐衰减，高程方向表现为分阶段隆升特征[2-3]。

青藏高原的变形与隆升，与印度板块、欧亚板块的汇聚碰撞过程密切相关。喜马拉雅造山带位于印度板块与欧亚板块交接处，是青藏高于最南端，长约2500 km，宽300～500 km，珠穆朗玛峰(以下简称“珠峰”)是其中最高峰，也是我国和世界最高峰。喜马拉雅造山带是印度板块和欧亚板块汇聚的主要变形带，是研究大陆构造形变机制最具代表性的区域。

青藏高原地壳隆升和变形机制目前仍存在很多争论。文献[4]提出“双倍地壳模式”，随后众多学者先后提出多种理论和动力学假说，主要有“地壳增厚模式”[5-7]“大陆逃逸模式”[8-11]“双向俯冲作用模式”[12]，我国科学家也先后提出“多边大陆俯冲和地幔底辟模式”[13-15]“大陆注入模式”[16-17]等理论假说。上述多种假说均是基于不同地质或地球物理观测现象提出的青藏高原隆升演化变形模式，目前仍存在争议。

随着空间大地测量观测技术的发展，GNSS呈现出高精度、大范围、准实时的优点，自1988年我国首次开展GNSS地壳形变监测以来，GNSS就越来越广泛地应用于地壳形变监测中，使得地壳运动与形变监测逐步实现高精度、连续动态监测，实现了由定性向定量监测的飞跃。青藏高原地壳隆升与变形也逐步实现连续、定量监测。

我国于1988年在滇西地震试验场，首次将GNSS技术应用到地壳运动与形变监测中[18]。1991年，中美合作首先在青藏高原东部及邻区，开展了GPS观测。通过1991—1995年的4期观测，定量地监测了青藏高原东南缘现今构造变形特征[19]。同期，中美还合作开展了跨喜马拉雅GPS观测。美方在尼泊尔布设28个GPS站进行了观测，中方在西藏南部布设了5个GPS站，并联合格尔木地震台进行了同步观测[20]。

我国独立开展青藏高原GNSS测量始于1992年，国家攀登计划项目“现代地壳运动和地球动力学研究”于1992和1994年两期联测在青藏高原获取了拉萨-温泉等基线[21]。1993—1995年，国家地震局在“八五”期间，在青藏高原开展了两期GPS观测[22]。“九五”期间，我国开展了“中国地壳运动观测网络”(CMONOC-I)工程建设，2000年12月通过国家验收，正式投入运行[23-25]。“十一五”期间，在中国地壳运动观测网络的基础上，开展了“中国大陆构造环境监测网络”(CMONOC-Ⅱ)工程建设，该网络于2011年正式投入运行，包括260个基准站和近2000个区域站[24-26]。

1992—2020年，我国与意大利、美国合作，先后于1992、1998、1999年开展了珠峰高程GNSS测量，我国于2005、2020年自主开展了珠峰高程GNSS测量观测[27-29]。2020年珠峰高程测量中首次在峰顶使用了国产GNSS接收机，利用北斗卫星信号进行了峰顶高程测量。

GNSS应用于青藏高原地壳形变监测，从多层次、多方面研究了青藏高原地壳形变及其动力学机制[30-45]。文献[36—37]利用GNSS数据，发现印欧板块碰撞使青藏高原地壳缩短并增厚，但碰撞90%以上被喜马拉雅山及其边缘形变所吸收，从而为地壳增厚学说提供了支持。文献[38—39]利用青藏高原的553个GNSS站点，研究了青藏高原构造运动特征，发现青藏高原现今主要构造运动属于连续介质变形过程。文献[40]基于GNSS观测结果研究了青藏高原的整体运动状态，揭示了青藏高原物质运移机理。文献[41]利用10多年累积的GNSS观测资料，对青藏高原现今地壳运动进行研究，探讨了青藏高原现今隆升状况及其动力学含义。文献[42]通过GNSS应变率场，研究了青藏高原内部更加重要的地壳减薄过程。文献[43]利用GNSS数据结合GRACE重力卫星数据研究了青藏高原东缘地壳流速度。文献[44]利用GPS结合GRACE监测了喜马拉雅地区地壳水平形变。文献[45]利用GRACE研究了青藏高原地壳垂直形变，并与GPS结果进行了对比。

本文利用1992—2020年历次珠峰高程测量GNSS复测点监测数据，并结合自1999年以来，20多年的青藏高原GNSS观测数据，对青藏高原及周边地壳形变、地壳隆升进行了定量研究，分析了青藏高原及周边地壳三维形变特征及其动力学机制。

1 GNSS数据处理

本文收集了1992—2020年我国历次珠峰高程测量建立的区域地壳监测网GNSS数据6个复测站GNSS数据，1999—2020年陆态网络位于青藏高原及周边30个连续观测的CORS站、352个定期观测的区域站，2009—2020年尼泊尔14个CORS站GNSS观测数据。站点分布如图1所示。本文利用GAMIT软件对GNSS数据基线解算，解算方案主要参数设置见表1。

图1 青藏高原及周边GNSS站点分布

表1 基线解算设置

解算获取各站单日解时间序列，单天解精度水平方向优于3 mm，高程方向优于5 mm。利用GLOBK软件进行整体平差，获取ITRF2014框架下GNSS站点三维线性速率。

2 GNSS监测点垂直形变精化

为研究青藏高原地壳隆升，首先需要精确获取地壳垂直运动速率。GNSS时序不仅反映线性的构造运动，还包含了非线性变化特征，这种非线性变化一般是由大气、陆地水等地表质量负荷引起，属于非构造负荷形变，特别是垂直方向时序主要特征是非构造负荷导致的周期性变化。GNSS观测站单天解时序如图2所示。有学者研究发现，GRACE能够监测一定空间尺度上的陆地水变化，利用GRACE扣除陆地水负荷形变，可以消除约64%的GNSS共模误差[46-47]。本文利用GRACE时变重力场模型，分析非构造负荷对GNSS高程时序的影响，排除陆地水等非构造负荷影响，进一步精化站点垂直运动速率。

图2 GNSS时序

本文所采用的GRACE数据产品为美国德克萨斯大学空间研究中心(CSR)发布的GRACE Level-2(RL06)月重力场模型数据。主要计算过程为：利用卫星激光测距(SLR)获得的低阶重力场球谐系数，替换时变重力场模型(GSM)数据的C20项及一阶项；利用Fan组合滤波P3M15处理，去除噪声和泄露误差等影响；利用前后月份数据内插，补足时变重力场模型中缺失月份数据；以所有月份的数据均值为基准，从月重力场模型中扣除基准；根据扣除基准后的重力场球谐系数文件反演获得陆地水储量变化；根据负荷形变理论和相关计算公式，利用陆地水储量变化，计算获得其产生的垂直形变。

图3展示了GARCE获取地壳垂直运动时序与GNSS高垂直方向时序，可以看出两者在周期、振幅上基本一致，这表明GNSS年周期主要是由非构造负荷形变引起的。

图3 GNSS垂直时序GRACE非构造负荷形变改正

GRACE对于水平方向线性运动影响很小，在研究水平方向构造运动时可以忽略不计。但是对于垂直方向影响很大，特别是对于大部分定期观测的区域站，由于数据量较少，且每期观测时间并不完全固定，利用GRACE对GNSS垂直方向时序进行改正，可以有效提升线性速率提取的准确度和精度。

图4展示了利用GRACE对垂直方向时序进行改正前后，GNSS垂直方向线性速率及其精度，可以看出，改正后GNSS站点垂直方向线性速率精度明显提升，平均提升约40%，其中对于定期观测的区域站提升更为明显，平均达到约49%。为从不同角度研究青藏高原地壳运动特征，本文还将ITRF2014框架下速率转换为相对欧亚板块速率。站点三维运动速率如图5所示。

图4 改正前后GNSS站点垂直运动速率和精度

图5(a)包括了ITRF2014框架下各个站点水平方向运动速率(蓝色箭头)和相对欧亚板块的水平方向运动速率(红色箭头)。ITRF框架下站点水平速率显示出明显的分区块特征，自南向北运动速率逐渐递减，东向运动速率逐渐增加。

图5 GNSS观测站三维速率

而在扣除欧亚板块整体运动趋势后，水平方向运动量级明显减小，运动方向分化更为明显，东部右旋更加明显，西部则出现左旋趋势。

3 基于块体模型的地壳形变分析

3.1 块体模型构建

GNSS站点分布较为分散，且分布不均匀。为研究区域地壳整体运动特征，需要通过建立模型，实现格网化。青藏高原地壳运动最显著的特征之一是具有明显的分区块特征[48-49]，因此本文采用三维块体模型，通过构建地壳三维运动与形变模型，研究地壳三维形变特征[50-53]

Vs=Rt·Ω

(1)

(2)

(3)

式中，vn、ve、vu表示块体任一点(λ,φ)南北向、东西向和垂直方向速率；r为地球平均半径；(λo,φo)为块体几何中心坐标；Δx、Δy、Δz表示块体x,y,z方向的整体位移；ωx、ωy、ωz是欧拉旋转矢量；εe、εn、γen分别为块体相对于块体几何中心的东向、北向主应变和剪应变。

块体划分采用参考我国CPM-CGCS2000板块模型及地质研究成果[53-55]。本文主要研究范围主要包括拉萨块体、羌塘块体、巴颜喀拉块体和柴达木块体。

本文分别利用ITRF2014框架下速度场和相对欧亚板块速度场构建了块体模型，利用块体模型拟合站点三维速度场(红色箭头)，与实测GNSS速度场(蓝色箭头)进行对比，评估模型精度。结果如图5、图6所示，精度统计结果见表2、表3。

图6 实测GNSS站点水平速率与模型拟合速率

表2 ITRF 2014框架下块体模型精度

表3 相对欧亚板块运动框架下块体模型精度

图6显示，模型拟合水平方向速度场，与站点实测速率趋势、量级基本一致。表2和表3统计了拟合速度场与实测速度场之间差值的均值和均方差，均值反映模型无偏性，均方差则反映模型有效性。可以看出，各个块体差值N、E、U 3个方向均值接近零，说明本文构建的模型具有很好的无偏性；水平方向均方差最大约6 mm/a，最小仅有1 mm/a，显示模型在水平方向具有良好的有效性。模型在巴彦喀拉块体、羌塘块体在南北向均方差较大的主要原因是块体运动出现明显的右旋趋势，导致南北向速度场发生较大变化，反映到模型上，表现为站点南北向速率离散度较大。与其他学者采用欧拉矢量法[55]、有限元法[56]构建的速度场模型对比，本文构建的模型在无偏性更好，有效性方面在东西方向均方差优于明显优于其他结果，南北方向均方差则基本一致。这表明，青藏高原地壳采用非连续活动块体模型，考虑块体内部应变构建速度场模型，更能体现地壳实际运动与应变特征。

图7显示，模型拟合垂直方向速率与站点实测速率，部分站点差别较大，对模型在垂直方向有效性有一定影响，但表2、表3显示，模型垂直方向均值接近于0，表明模型无偏性较好，能够反映地壳整体垂直运动趋势。

图7 实测GNSS站点垂直速率与模型拟合速率

3.2 基于块体模型的地壳水平形变特征

块体模型可以反映地壳三维运动与应变特征。图8展示了模型构建的地壳水平形变场，包括格网速度场、块体间相对运动。图8蓝色箭头表示地壳水平速度场，红色箭头表示块体间相对运动，底图为水平运动角度。

图8 基于块体模型的地壳水平运动

图8(a)展示了ITRF框架下块体之间的相对运动。新生喜马拉雅块体和拉萨块体之间西段主要为东西向相对运动，中段则具有明显的挤压状态，东段又转为东西向滑动；拉萨块体与羌塘块体之间同样表现为明显的挤压，尤其是东段量级较大；羌塘块体和巴颜喀拉块体之间则主要表现为东西向相对滑动；巴颜喀拉块体和柴达木块体之间主要呈现挤压运动。

图8(b)显示了去除欧亚板块整体运动趋势后，青藏高原地壳水平运动空间分布特征，总体趋势基本一致，相对运动速率量级有所减小。

本文通过截取ITRF框架下东经90°、北纬33°剖面展示不同区域地壳运动量级、方向分布(图9)，喜马拉雅山脉南部运动量级约为50 mm/a，方向大概在NE40°，向北量级有所减弱，拉萨板块和羌塘板块运动为40～45 mm/a，运动方向则有明显变化，由NE30°逐渐旋转为大概在NE0°，至巴颜喀拉块体、柴达木块体则北向速率略有减小，东向速率基本稳定，量级在30～40 mm/a，方向在NE10～20°。

图10展示相对欧亚板块地壳水平运动空间变化特征，喜马拉雅山脉南部运动量级约为40 mm/a，方向达到NE70°，向北量级逐渐减弱，拉萨板块和羌塘板块运动为20～30 mm/a，运动方向则有明显变化，逐渐旋转为大概在NE30°，往北至巴颜喀拉块体、柴达木块体则北向速率略有减小，东向速率基本稳定，量级约在20 mm/a左右，方向则明显增大。

图9、图10中，北纬33°剖面图显示在90°～98°之间羌塘块体有一个明显速率下滑的过程，结合图8可以看出，该区域正是地壳沿玉树-鲜水河断裂逐渐右旋转向的过程，2010年玉树Ms7.1级地震就发生在这个区域。

图9 ITRF2014框架下地壳水平方向运动量级与角度

图10 相对欧亚板块框架下地壳水平方向运动量级与角度

不同框架下GNSS运动速率存在一定差异，而基于速度场的地壳应变张量则与框架选取无关，可以从不同方面反映地壳形变特征。图11展示了面膨胀率、最大剪应变率和最大最小主应变。

图11(a)面膨胀率显示，在印度板块与欧亚板块交界的喜马拉雅造山区域，存在明显的挤压变形，特别是在2015年尼泊尔地震所在地区到加德满都周边，挤压变形尤其明显。而在青藏高原中部的拉萨、羌塘、巴颜喀拉中部则面膨胀率不明显，平均约4.0±1.0nonstrain/a。

图11(b)显示，在剪应变率和最大最小主应变方面，加德满都周边形成一个明显的应变集中区域。在拉萨块体，主应变表现为南北向压缩与东西向拉伸，特别是在加德满都-日喀则方向是最大主应变方向，该方向面膨胀率为负值，表明该方向是拉萨块体主要挤压变形区域。

图11 地壳应变率场

3.3 基于块体模型的地壳垂直形变特征

图12展示了模型构建的地壳垂直形变场。不同框架下，垂直形变基本一致，相对欧亚板块框架下，模型垂直方向形变速率略低于ITRF2014框架下，而分布一致，包括块体之间垂直方向相对运动也基本一致。因此，在研究垂直形变时，基本可以不考虑框架影响。

图12显示，青藏高原总体隆升趋势明显。位于印度板块部分，垂直方向速率为负值，向北垂直方向逐步变为正值，地壳逐渐呈隆升特征，拉萨、羌塘块体整体隆升速率为2～3 mm/a，巴颜喀拉块体隆升不明显，平均速率低于1 mm/a，柴达木块体同样隆升不明显，而祁连块体则出现明显上升的趋势。

图12 基于块体模型地壳垂直形变

为研究青藏高原地壳整体形变特征，本文通过计算垂直速率与水平速率比值，研究了地壳运动倾斜度(图13)。倾斜度可以在一定程度上反映地壳在运动过程中出现的弯曲程度。由图13可以看出，在青藏高原中部拉萨块体和羌塘块体，特别是在两个块体的西部和中部，倾斜度明显高于周边地区，表明在该区域地壳运动中出现一定程度向上弯曲的变形特征。

图13 地壳运动倾斜度

本文截取了90°E、95°E，31°N、34°N等剖面，展示不同剖面地壳垂直运动速率分布，以及倾斜度等(图14)。

图14 不同经纬度剖面地壳运动垂直速率与倾斜度

从各个经度剖面可以看出，青藏高原在中部都有一个明显的隆起部分，可能的解释是青藏高原受南北压缩出现中部隆升大于两端的向上弯曲现象。从纬度剖面可以看出，青藏高原西部主要表现为整体隆升，而在东部则向东地壳隆升逐渐下降。

4 结论与讨论

目前，包括CORS站和定期观测站在内的GNSS监测网已覆盖了青藏及其周边地区，为精确定量研究青藏高原地壳形变、隆升及其动力学机制提供了丰富、可靠的数据资料。

本文基于GNSS长期监测数据，利用块体模型从多个方面分析了青藏高原地壳三维形变特征，探讨了地壳形变与隆升机制。

顾及块体内部应变的活动块体模型，模型拟合结果与实际GNSS观测结果差值均值接近于0，具有优良的无偏性，其有效性方法水平和垂直方向均方差均优于1 cm/a，可以在大区域内精确反映地壳实际运动特征。

GNSS水平速度场模型，特别是相对欧亚板块速度场表明，青藏高原受到印度板块南东方向挤压、塔里木块体北东方向挤压，从而出现西部边缘西向旋转挤出、中部和东部东向旋转的整体趋势。

利用GRACE精化GNSS垂直速度场，可以更精确地获取地壳垂直方向构造运动趋势，其中喜马拉雅山南北两侧临近区域，GNSS垂直速率改正后略有降低，而在羌塘块体、巴彦喀拉块体、柴达木块体等区域，改正后垂直速率有所上升，改正后垂直速率精度整体提升约40%，大部分GNSS站点垂直方向长期运动趋势均为隆升。

基于块体模型的地壳三维形变场显示，青藏高原南部拉萨块体，受到南北向挤压，南北向挤压收缩是块体主要应变特征，特别是中段具有明显的应变集中现象，2015年尼泊尔Ms8.1级地震就发生在该区域。同时南北向挤压也导致地壳隆升相对较大。

位于青藏高原中部的羌塘块体，受到拉萨块体和塔里木块体挤压，其水平方向自西向东沿黎嘉断裂右旋运动，东部则受到川滇块体阻挡。块体西部存在一定南北向挤压应变，中部和东部则没有明显的南北向挤压应变，地壳表面积略有膨胀；高程方向，羌塘块体地壳隆升总体大于周边块体，平均约2.0 mm/a，其垂直速率与水平速率比率大于周边地区，表现为青藏高原地壳运动过程中，在中部地区出现向上突出隆起的垂直形变特征。

巴颜喀拉块体整体东向旋转，其东部龙门山地区挤压作用强烈，为强烈挤压应变区；西部边界以伸展应变为主，北边界主要表现为东向滑移。水平运动整体具有顺时针旋转特征。地壳隆升，平均约为0.8 mm/a。

柴达木块体受塔里木块体东向挤压，整体以东向运动为主，同时又收到祁连块体北向挤压，导致块体整体略有南北向略有收缩，地壳隆升速率较巴彦喀拉块体略大，约1 mm/a，但低于祁连块体的约1.8 mm/a。

研究表明，青藏及周边地壳形变具有明显的分区特征。南北向收缩主要表现在拉萨块体，东西向伸展则主要表现在巴颜喀拉块体，中部羌塘块体没有明显南北向水平挤压应变，有长期东西向伸展趋势。垂直方向总体隆升，特别是在中部羌塘块体，地壳隆升速率最大、水平面上面积有扩大趋势，反映出青藏高原地壳有增厚的可能，这在一定程度上支持了地壳增厚学说。

致谢：感谢自然资源部第一大地测量队提供的珠峰高程测量区域监测网GNSS数据，感谢“中国地球运动观测网络”“中国大陆构造环境监测网络”提供了青藏高原GNSS观测数据，感谢UNAVCO提供了尼泊尔CORS站观测数据，感谢GAMIT/GLOBK软件和GMT软件团队。