利用GPS高精度监测数据开展青藏高原现今地壳运动与形变特征研究进展

瞿 伟,高 源,陈海禄,梁世川,韩亚茜,张 勤,王庆良,郝 明

(1. 长安大学 地质工程与测绘学院，陕西 西安 710054； 2. 中国地震局第二监测中心，陕西 西安 710054)

0 引 言

青藏高原是世界范围内地壳构造活动最为强烈的区域之一，受印度洋板块北推挤压动力学作用，青藏高原不断抬升并在其南缘形成了世界最高峰——珠穆朗玛峰(以下简称“珠峰”)，在大陆内部板块构造运动形变研究中具有重要地位，也是国际上开展地球动力学研究的热点区域[1-3]。在印度洋板块北推挤压与欧亚板块阻挡作用下，青藏高原内部及周缘发育有数条大型活动断裂带，如喜马拉雅碰撞带、雅鲁藏布江缝合带、班公湖—怒江缝合带、甘孜—玉树—鲜水河断裂带、东昆仑断裂带、阿尔金—海原断裂带等[4-5]。这些断裂构造活动导致青藏高原强震频发，曾发生过1920年海原Mw8.5地震、1950年墨脱Mw8.6地震、1976年龙陵Mw7.3地震、1988年澜沧—耿马Mw7.3地震、2001年昆仑山Mw8.1地震、2013年定西Mw6.0地震、2015年尼泊尔Mw7.9地震与2020年西藏尼玛Mw6.3地震等。因此，开展青藏高原构造活动研究对于认识该区域现今地壳形变与孕震机制具有重要的科学意义。

以GPS为代表的空间大地测量技术可以获取地表高精度、高时空分辨率的形变信息，在地壳形变监测研究中得到了广泛应用。自“九五”时期以来，中国地壳运动观测网络已在青藏高原及邻区布设了大量GPS观测站[6]，近20余年来持续高精度的监测为青藏高原地壳形变与动力学研究提供了丰富的基础数据。GPS速度场结果直观揭示了青藏高原SN向缩短与EW向伸展形变的现今构造形变模式[7-14]，地壳运动监测结果进一步结合运动学与动力学模型反演获得了区域内部应变与断裂运动特性及深部构造活动特征[15-27]。

本文系统总结了青藏高原及邻区GPS监测网络体系的构建与GPS监测高精度数据处理方案策略，详细阐述了当前研究青藏高原地壳运动形变的运动学和动力学模型，并结合研究成果深入剖析了青藏高原现今浅部与深部地壳构造运动与形变特征，从整个中国大陆构造动力学背景框架下进一步详细解译了青藏高原现今地壳构造活动动力学机制及区域动力学背景影响下的珠峰地区隆升机制，最后对利用GPS技术开展青藏高原地壳运动形变与动力学机制研究成果进行了总结，并对该技术与其他空间监测技术联合开展青藏高原构造活动研究的未来发展趋势进行了展望。

1 GPS监测网络与高精度数据处理

1.1 GPS监测网络体系

地壳形变过程中的精细运动图像可为研究地球动力学及地震孕育机制等相关地球科学问题提供十分重要的基础资料，得益于以GPS技术为代表的空间大地测量技术在大尺度、高精度观测地壳运动与形变状态的优势。中国地震局联合多部门于“九五”期间投资建设了国家重大科学工程“中国地壳运动观测网络”(CMONOC-Ⅰ)，其中GPS监测网络体系由27个连续观测站和1 056个非连续观测站构成，于1999年开始正常运行。

为进一步提升对中国大陆地壳形变监测能力，在中国地壳运动观测网络的基础上，中国地震局、中国人民解放军总参谋部测绘局、中国科学院、国家测绘局、中国气象局和教育部于“十一五”期间投资建设“中国大陆构造环境监测网络”(CMONOC-Ⅱ)。该网络于2009年开始正常运行，目前已有260个基准站与近2 000个区域站(图1)[27-28]，可为监测中国大陆地壳运动提供高精度、高时空分辨率的基础研究数据[29]。此外，鉴于中国大陆内部局部地区地震发生危险性较高，如华北/首都城市圈、四川盆地、山西断陷带、鄂尔多斯地块、青藏高原及其周边地区等，各省市地震局及相关科研院所在原有监测网络基础上进行了加密观测，特别是在青藏高原及其周边地区布设了大量加密监测站点，如中国地震局地质研究所联合四川省地震局在四川盆地和鲜水河断裂带布设了20个连续观测站，四川省地震局在四川盆地布设了17个连续观测站，中国科学院青藏高原研究所在藏南地区布设了6个连续观测站，中国气象局在青藏高原及其周边地区布设了23个连续观测站，美国加州理工学院在喜马拉雅地区布设了28个连续观测站，中国地震局第二监测中心在青藏高原东北缘布设了23个流动观测站，中国地震局地质研究所在云南地区布设了59个流动观测站等。

图1 青藏高原及其周边地区GPS连续观测站与流动观测站位置分布Fig.1 Distributions of GPS Continuous and Campaign Observation Stations in Qinghai-Tibet Plateau and Its Surrounding Areas

1.2 GPS监测资料高精度数据处理

目前国际上GPS数据处理主流软件主要有GAMIT/GLOBK、GIPSY/OASIS、Bernese与PANDA等。本文主要总结介绍利用GAMIT/GLOBK软件进行青藏高原GPS数据处理的流程与策略，数据处理大致分为4个主要步骤[16,30-31]：第一步，采用GAMIT软件中单日松弛解处理基线，采用sp3精密卫星星历削弱轨道误差，卫星轨道约束为10-8，光压模型采用BERNE模型，并进行海潮与固体潮改正，选用中国大陆及周边的国际IGS跟踪站(如BJFS、KUNM、WUHN、ULAB与CHAN等)进行联合解算获得基线单日松弛解h文件；第二步，利用GLOBK软件将区域单日松弛解文件与SOPAC(Scripps Orbit and Permanent Array Center)全球IGS跟踪站单日松弛解合并获取包含所有GPS观测站的全球分布松弛解；第三步，利用QOCA软件统一平差获取ITRF(International Terrestrial Reference Frame)参考框架下各站点GPS速度场(图2)，目前ITRF参考框架已更新到ITRF2014[32]；第四步，为更好展示出青藏高原内部的形变状况，进一步计算获得青藏高原相对于欧亚板块运动速度场，同时顾及观测时段内大震同震及震后对速度场结果影响，需进一步基于同震形变场模型进行数据处理，获取干净可靠的震间形变场[13,28,30]。

图3为中国地震局第二监测中心郝明副研究员基于青藏高原1999～2007年GPS监测数据，按照上述高精度数据处理步骤计算获得的相对于欧亚板块基准下的青藏高原及其周边地区GPS速度场[13]。该速度场清晰揭示出：印度洋板块北推挤压动力学作用下青藏高原具有NNE向运动趋势，在青藏高原西南缘环喜马拉雅造山带地壳挤压运动速率约为35 mm·年-1；由于受其北部稳定准格尔盆地与中亚造山带阻挡，青藏高原地壳运动速度向北逐渐减小，在青藏高原北缘地壳运动速率逐渐减小至约14 mm·年-1；又受其东北缘鄂尔多斯地块阻挡，地壳运动速度场由近EW向转为SEE向，同时受稳定四川盆地的阻挡，地壳运动速率减弱至约7 mm·年-1；在稳定华南地块阻挡下，青藏高原东南缘地壳运动速度场由北部近EW向转为SSE向，地壳运动速率约为9.4 mm·年-1，同时推挤华南地块朝SEE向逸出。

图件引自文献[13]图3 青藏高原及其周边地区GPS速度场(欧亚板块基准，1999～2007年)Fig.3 GPS Velocity Field in Qinghai-Tibet Plateau and Its Surrounding Areas (Eurasian Plate Benchmark, 1999-2007)

尽管相对于欧亚板块基准下的GPS速度场可较好地揭示青藏高原及其周边地区整体运动特征(图3)，但难以直观表征出区域内块体间运动形变差异，近年来相关研究也采用了区域参考基准获取区域内块体间形变差异特征。Ge等探讨了相对于印度洋板块基准下GPS速度场较欧亚板块基准可更加清晰地揭示出青藏高原与印度洋板块汇聚作用、青藏高原EW向伸展形变，特别是沿喜马拉雅东部的顺时针旋转作用[24]；邹镇宇等研究了南北地震带相对于稳定华南地块基准下GPS速度场在汶川地震前后的形变差异与主干断裂相对运动差异[33]。图4为青藏高原相对于华南地块基准下GPS速度场，青藏高原表现为显著SN向缩短与EW向伸展形变。青藏高原西南缘运动速率为近SN向约36 mm·年-1，西北缘运动速率减小为NE向约17 mm·年-1；受鄂尔多斯地块阻挡作用，其东北缘逆时针朝北运动，地壳运动速率约为5 mm·年-1，而东南缘相对华南地块则表现为顺时针旋转运动。

2 地壳运动与形变分析模型

2.1 运动学方法

区域地壳运动与形变分析的运动学方法可分为连续与非连续两种。连续分析方法主要将整个研究域视作一整体，进行应变率、滑动速率等反演计算，获取区域整体运动特征；而非连续分析方法则根据先验地质资料，将研究域进行块体单元划分，分别获取各单元运动与形变特征，然后再整合区域整体运动与形变特征。

图件引自文献[13]图4 青藏高原及其周边地区GPS速度场(华南地块基准，1999～2007年)Fig.4 GPS Velocity Field in Qinghai-Tibet Plateau and Its Surrounding Areas (South China Block Benchmark, 1999-2007)

2.1.1 连续分析方法

(1)多尺度球面小波模型。该模型基于小波理论从平面空间向球面空间延伸，已在地球重力场[34]、地壳运动速度场[35]等相关地球科学研究中得到了广泛应用。球面小波模型分析方法借助其对空间和频率局部化的优势，通过调整分解尺度控制信号频率信息，将观测数据所反映出的不同空间尺度地壳运动与形变信息进行刻画。大尺度对信号的整体特征更为敏感，主要体现出低频信息；而小尺度则对信号的局部特征较敏感，主要体现出高频信息。

以高斯函数差(Difference of Gaussians)构建的球面DOG小波[36]作为球面小波模型的基函数，其表达式为

(1)

相对于基函数，通过位置和尺度离散化得到的球面小波框架可更加灵活表示信号，一般可选择正十二面体球面三角剖分方式，将不同尺度因子取值所对应球面剖分的格网点作为中心极点，可得到不同尺度下球面小波框架(W)，其表达式为

W={ψo(q,i),a(o′),o(q,i)∈Gq,

q=0,1,2,…,qmax}

(2)

式中：o(q,i)为不同尺度对应的球面小波中心极；i为格网点序号；Gq为格网点集合。

考虑到球面格网点数会随着剖分层次而呈指数增长，按照“空间支撑”思想建立以站点分布为基础的框架选取方式。站点分布范围越大，尺度因子最小值(qmin)越小；站点分布密度越大，尺度因子最大值(qmax)越大。

GPS观测站点可用三维表示，以式(2)中得到的球面小波框架作为估计GPS速度场框架，可将速度场表示为

(3)

图5为不同尺度因子对应的球面小波函数图像。当q=3时，中心极为北极。

需要注意的是，上述模型是通过对球面小波函数的位置和尺度进行离散化提取小波框架，该策略可能存在冗余，进而引起速度场存在不唯一解的问题。此外，观测误差和数据分布异常也可能导致病态问题，但通过正则化方法可保证获得上述模型向量H的稳定、唯一解[37]，其表达式为

(4)

式中：L为已测站点的GPS观测值；Bnn为GPS观测向量的观测误差自协方差阵；ρ为正则化参数；G为模型系数阵；R为正则化矩阵；H为ak、bk和ck组成的模型向量。

(2)球面最小二乘配置模型。该模型将地壳形变分为刚体旋转和随机形变信号两部分，在观测数据信号空间连续性良好的前提下，可获取大区域应变场的低频分布，继而反映区域应变场整体特征。该模型解具有较高的抗差性和稳定性特征[20]，其函数模型为

L=AΩ+CZ+Δ

(5)

式中：L为GPS观测值；A为欧拉系数矩阵；Ω为欧拉参数；C由单位矩阵和零矩阵组成，分别对应已测点和未测点的形变信号；Z为形变信号；Δ为观测信号的噪声。

(6)

(7)

B=Boo+Bnn

(8)

式中：Boo为已测点信号协方差矩阵；Bou为未测点信号与已测点信号的协方差矩阵。

通常未测点与已测点信号的协方差是未知的。相关研究表明，在观测区域较大和观测数据较多的情况下，可利用信号统计特征提取形变信号的经验协方差分布函数。常用的经验协方差分布函数有高斯函数[40-41]、指数函数[42]和顾及距离尺度因子的映射函数[43]。其中，高斯函数为最常用的经验协方差分布函数。其模型形式为

f(d)=f(0)e-k2d2

(9)

式中：d为观测站点间距(以千米为单位)；f(0)和k为未知参数。

经验函数模型参数通过在给定距离范围内建立信号协方差与测点距离的统计关系，并基于最小二乘方法求解。通过上述流程可获取研究域连续的形变场数据，建立随机形变信号与位置函数关系，进一步微分可求解获得区域地壳形变特征的空间分布。

(3)反距离加权算法。该算法通过对不同测区速度给予不同权重，自动化实现观测数据分区处理[44]。反距离加权算法在二维空间以任意小的增量进行迭代，确保连续应变函数可解，在其插值位置

图5 不同尺度因子对应的球面小波函数图像Fig.5 Images of Spherical Wavelet Function in Different Scales

R处，GPS水平速度场被扩展为一阶导数，用刚性块体运动和均匀应变模型表示。线性函数模型为

L=Kh+Δ′

(10)

式中：h为未知参数向量，包括平移、旋转和应变参数；K为设计矩阵；Δ′为观测误差。

仅考虑水平形变时，h由平移分量(Ux和Uy)、应变分量(τxx、τxy和τyy)和旋转率(ω)组成。上述模型参数可通过式(11)求解。

(11)

进一步考虑到大区域范围内对加密观测数据插值的需求，基于“观测点距插值点越近，权重越大”的思想，对初始观测误差方差进行重构。其表达式为

(12)

(13)

(14)

式中：加权函数Gi=Li×Zi，Li和Zi的取值分别依赖于距离和观测数据的空间覆盖率[45]；D为平滑因子；Ri为观测点i和插值点的水平距离；ni为所选定的用于计算插值点位移的数据个数；θi为插值点附近选择的第i个和第i+1个数据之间的方位角。

(4)刀刃模型。Savage等提出了以地表位移为约束，断层滑动和闭锁深度为反演目标的断层剖面拟合弹性模型[46]。该模型假定断层滑动发生在闭锁深度以下。其表达式为

(15)

式中：V为平行于断层的速率；x为台站与断层的垂直距离；E为闭锁深度；V∞为断层滑动速率。

模型参数可通过格网搜索获得，其不确定性常采用蒙特卡洛随机试验结果统计分析做出评价[47]。弹性模型假设断层附近形变梯度大，远离断层形变梯度变小(图6)。

x代表平行于断层方向；y代表垂直于断层方向；m为最大有效距离；图件引自文献[46]图6 走滑断层形变模型示意图Fig.6 Schematic Diagram of Strike-slip FaultDeformation Model

考虑到平行断层间的相互作用，Savage等于1999年提出了对应的改进模型[48]。其表达式为

V(x)=b0-{b1arctan[(x-x1)/h1]+

b2arctan[(x-x2)/h2]+

b3arctan[(x-x3)/h3]}/π

(16)

式中：h1、h2和h3为各个断层闭锁深度；b1、b2和b3为各个断层在x1、x2和x3位置处的滑动速率；b0是与参考框架相关的常数；模型参数可用最小二乘拟合得到。

在上述模型基础上，Segall考虑到部分断层受浅层蠕滑影响[49]，对上述模型进行了进一步改进。其表达式为

(17)

式中：bicreep为断层蠕滑速率；hicreep为断层蠕滑深度；模型参数常用粒子群优化算法反演得到[50-51]。

2.1.2 非连续分析方法

(1)活动块体模型。区域活动块体一般可被视为刚体[52]，其函数模型[53]为

(18)

式中：ve和vn分别为板块任一位置(λ，φ)的东向与北向GPS观测值；r为地球半径；ωx、ωy和ωz为块体欧拉参数。

实际上，活动块体会受到周缘块体的挤压、碰撞，其整体运动不仅会发生旋转，其边缘和内部也会发生弹性形变。假设块体内部形变是均匀的，则地壳形变模型为整体旋转与均匀应变模型(REHSM)。其函数表达式[31,53]为

(19)

式中：εee为EW向线应变;εen为EW—SN向剪应变;εnn为SN向线应变；x0=rcosφ(λ-λ0)，y0=r(φ-φ0)，(λ0，φ0)为块体中心位置。

若进一步假设块体内部形变是非均匀的，可将其描述为位置与应变的线性函数，则地壳形变模型可表征为整体旋转与非线性应变模型(RELSM)。其函数表达式[31,53]为

(20)

式中：A0、B0、C0、ξ1、ξ2、ξ3、S1、S2、S3分别为块体应变参数。

上述3类模型均是在活动块体划分的基础上进行计算分析，考虑到活动块体内部形变特征具有的未知特性，通常需进一步利用统计理论对模型对应应变分量进行参数显著性假设检验[54]；线性假设为H0∶Y=0，当H0成立时，可构造F统计量为

(21)

对应的拒绝域为

F>F6,2n-12

(22)

(2)负位错模型。负位错模型假设地壳形变是由活动块体整体旋转、内部均匀应变和边界断层闭锁产生的负位错效应等3个部分组成的地表弹性形变之和(图7)[55-56]，考虑了活动块体的旋转效应、内部永久形变和Okada位错。其表达式为

Vsf=Vbr+Vis+Vfs

(23)

式中：Vsf表示地表形变；Vbr表示块体整体旋转引起的形变；Vis表示块体内部均匀应变引起的形变；Vfs为断层闭锁负位错效应引起的形变。

图件引自文献[57]图7 负位错模型示意图Fig.7 Schematic Views of Negative Dislocation Model

可以联合GPS速度场、InSAR数据、地震和地质资料作为约束，采用格网搜索法和模拟退火法反演块体旋转欧拉极、断层滑动速率和震间闭锁系数。

对于弹性形变部分，可以将断层进行离散化，形成一系列节点。每个节点是一个位错点源，每个点源的位错[u]是由节点处的闭锁程度(θ)和滑动矢量(Vθ)共同决定。其表达式为

[u]=-θ×Vθ

(24)

闭锁程度的取值范围为0≤θ≤1。当θ=1时，Vs=0,断层处于完全闭锁状态，应变积累能力达到最强；当θ=0时，Vs=Vb,块体边界处于自由滑动状态，无弹性形变积累。

相邻节点之间的断层分割为微小格网，选择双线性插值得到格网处的闭锁系数，从而得到断层的连续闭锁分布。负位错模型反演的条件方程为

Vk(X)=[ΩG×X]k+[ΩB×X]k+εkkDk+

(25)

式中：X为观测点坐标；Vk(X)为观测点X的GPS观测值；k为速度分量；ΩB为块体相对于参考框架的角速度；ΩG为GPS速度场相对于参考框架的角速度；εkk为应变参数；D为应变偏差；ΩF为断层下盘相对于上盘的欧拉极；N为断层走向的节点数；Qi为第i个节点的坐标；φi为第i个节点的闭锁系数；Gjk(X,Xi)为格林函数；j为断层面方向分量。

2.2 动力学方法

基于地表观测数据的运动学模型可较好地表征地壳浅部运动与形变特征，但难以深入揭示岩石圈深部动力过程及活动机理，因此，需进一步建立动力学分析模型探讨岩石圈深部构造活动影响下的地壳运动与应力应变特性。数值模拟方法顾及了岩石圈几何结构与物理结构的不均匀性以及主干构造断裂的空间展布特征，被广泛用于地壳形变与应力特征数值模拟中[58-62]。随着力学理论的成熟与计算机技术的迅速发展，模型几何结构也逐步由二维发展到三维(图8)，力学模型也由简单的弹性材料发展到考虑材料应力随时间松弛的黏弹性材料，再进一步发展到黏弹塑性材料。对于黏弹塑性材料而言，其应变增量[63]可表示为

{dε}={dεe}+{dεv}+{dεp}

(26)

{dεe}=[D]-1{dσ}

(27)

{dεv}=[Q]-1{σt}dt

(28)

(29)

式中：{}代表应变张量的6个分量；ε代表应变张量；上标e、v与p分别代表弹性、黏性与塑性应变；σt为t时刻的应力；dσ为应力增量；[D]为与弹性模量相关的弹性材料矩阵；[Q]为与黏滞系数相关的黏性材料矩阵；Gp为塑性势函数；dλ为塑性乘子。

断裂是地壳物质的重要组成部分，对地壳运动与岩石圈应变能累积有重要影响，因此，对断裂带的处理是数值模拟中的难点所在。目前大多数研究采用降低断裂带弹性强度的方法使其在模拟过程中更容易发生形变[64-66]，但这种物性参数的不连续性会导致断裂与地壳接触面出现人为的应力集中现象，可能未必与实际的力学情况相符[67]。另一种方式是采用摩擦接触分析，即认为断层上盘与下盘之间存在摩擦关系，采用库伦滑动定律表征接触面之间的黏着和滑移行为[68-70]，当滑动面上的剪应力超过极限应力，接触面发生相对位错；这种方法比采用弱化带更真实地模拟断层的相对运动特征，但在实际运算中易遇到模型收敛时间过长甚至不收敛的问题[63,71]。在黏弹塑性数值模拟中，给断裂带赋予较低的屈服强度而非弹性强度，采用断层区域的稳态塑性蠕变模拟断层长期滑动，当累积应力达到屈服强度时在断层外产生塑性应变。该方法已被广泛用于青藏高原的地壳形变与地震循环演化研究中[62,72-75]。

图件引自文献[75]图8 青藏高原三维有限元模型Fig.8 Three-dimensional Finite Element Model ofQinghai-Tibet Plateau

3 现今地壳构造运动与形变特征分析

3.1 浅部地壳运动与形变特征

考虑GPS速度场数据在研究区域分布的不均匀性，本文进一步反演计算了整个青藏高原及其周边地区地壳应变场分布特征[76]。

图9(a)显示了青藏高原及其周边地区的最大剪应变、主应变张量空间分布特征，其中环喜马拉雅和天山中部区域挤压形变最为显著。环喜马拉雅弧主应变方向近似垂直于板块边界走向，较好地反映了印度板块对欧亚板块的俯冲推挤作用；天山中部主应变近SN向分布，较好地反映了地壳近SN向缩短。最大主应变高值区的最显著区域位于喜马拉雅东构造结，平均值约为48 nonstrain·年-1，天山中部、环喜马拉雅弧和青藏高原东南缘也是最大剪应变高值区，平均值为25～45 nonstrain·年-1，由此反映了上述区域内现今地壳剧烈的构造形变特性。

图9(b)显示了青藏高原及其周边地区面应变分布特征。面应变正值表示地壳下陷运动，而面应变负值表示地壳隆升运动。从图9(b)可以看出：环喜马拉雅弧和天山中部面应变约为-80 nonstrain·年-1，预示着印度板块俯冲作用下青藏高原现今整体抬升的运动特性；青藏高原东北缘面应变约为-20 nonstain·年-1，表明该地区地壳也呈隆升状态，预示着青藏高原北向运动受阿拉善地块和鄂尔多斯地块阻碍而形成的地壳抬升作用；青藏高原中部部分地区和东南缘面应变约为20 nonstrain·年-1，表明这些区域地壳以下陷运动为主。

图9(c)显示了青藏高原及其周边地区地壳旋转率分布特征，有助于理解区域构造形变模式。青藏高原中部和南部地壳呈顺时针旋转运动，揭示了青藏高原物质的东向挤出作用，其北部、东南缘和天山中部地壳则呈逆时针旋转；喜马拉雅东构造结相对于欧亚板块顺时针旋转率最高可达60 nonstrain·年-1，青藏高原东南缘逆时针旋转率最高则可达30 nonastrain·年-1。

青藏高原在印度板块碰撞北推挤压作用下呈NNE向运动，位移由南向北逐渐减小[8，10]，其东南缘位移围绕东喜马拉雅顺时针旋转，反映了地壳SN向缩短和SE向逃逸的特征，区域运动与形变整体受控于印度板块北推挤压、印缅俯冲带深源俯冲以及缅甸微板块与巽他板块后撤的共同作用[77]；青藏高原西北缘构造形变主要通过天山地区地壳近SN向缩短而吸收调节[27]，而东北缘主要受周边阿拉善地块和鄂尔多斯地块阻挡呈NE向运动并进行顺时针旋转。值得注意的是，青藏高原内部大部分地区形变是连续的，只在垂直于欧亚板块碰撞方向的大型走滑断层上存在较大的形变梯度，但青藏高原东侧和西侧均发生了横向挤压，西侧挤出最高达6 mm·年-1，而东侧挤出最高达20 mm·年-1，东南部和南部应变为10～20 nanostrain·年-1，呈明

图9 青藏高原及其周边地区连续形变分布Fig.9 Distributions of Continuum Deformation of Qinghai-Tibet Plateau and Its Surrounding Areas

显扩张趋势，围绕喜马拉雅东部则主要做顺时针旋转运动[11，27]，但东侧挤出又受到华南地块的阻挡，地壳运动速度逐渐减小，在华南地块西部边界趋于0[33，78]；青藏高原喜马拉雅弧形区域在欧亚板块和印度板块的汇聚方向上受强烈挤压作用，其压缩应变高达30～60 nanostrain·年-1。

从应变特征来看，青藏高原主应变从西向东呈现出由NW向逐步转化为北向和NE向的变化特征；青藏高原东部断层应变积累速率明显高于西部，表现出显著的近SN向拉伸特征，但东南缘由于受华南地块阻挡，在多条大型断裂的控制作用下应变强烈，主应变方向变化复杂，EW向拉张扩散和SN向挤出压缩特征十分明显[77]。此外，鲜水河断裂带应变积累速率是青藏高原内最高区域[23]，值得进一步关注；青藏高原东北缘主应变主要表现为NE—SW向挤压和NW—SE向拉伸的特征，应变分布较均匀[79]，但剪应变在海原断裂带、东昆仑断裂带、鲜水河断裂带区域则非常集中，海原断裂带剪应变一直处于高值区，但剪应变在鄂尔多斯地块附近则迅速减小，原因可能是其运动形变被六盘山转换断层所吸收[26，79]，与之相反的是东昆仑断裂带剪应变由西向东逐渐减小，直至西秦岭附近消散。面膨胀率结果也同样表明青藏高原东向挤出在东昆仑断裂带附近被吸收，地壳明显增厚，而海原断裂带受到北部岩石圈阻挡，也呈现出高度收缩形态。青藏高原中西部存在显著正面膨胀区，而在喜马拉雅弧形区域和天山地区则存在较显著的负面膨胀区[27，40]；青藏高原北部和南部垂向应变分别为(8.9±0.8)和(7.4±1.2)nanostrain·年-1，进一步表明了青藏高原内部非常重要的减薄动力学过程[24，27]。

块体运动模型研究结果表明：青藏高原平均闭锁深度约为21 km[80]，东南缘各块体欧拉极处于93°E～96°E和22°N～25°N内，围绕东喜马拉雅顺时针旋转；此外，青藏高原内各断层走滑速率明显高于中国大陆其他区域，最大滑动速率位于青藏高原南边界[21，80-81]，反映了印度—亚洲板块挤压作用；喜马拉雅断裂带挤压速率由东至西从22～18.3 mm·年-1逐渐减小，东部断裂带主要呈走滑运动特征[22，80]，走滑速率最高可达40 mm·年-1，东部边界龙门山断裂带的挤压速率为3.2 mm·年-1，右旋走滑速率为1.2 mm·年-1，而川滇地区主要受鲜水河断裂带、安宁河和小江则木断裂的控制，以走滑剪切形变为主，表现出明显北强南弱、西强东弱的特点，鲜水河断裂带作为边界断裂，其左旋滑动速率为10～11 mm·年-1，对东喜马拉雅物质的北向扩张进行了有效限制[16]；青藏高原内部喀喇昆仑断裂带主要呈右旋走滑趋势，走滑速率为11.6 mm·年-1，挤压速率为4.3 mm·年-1，但其东段表现为相对一致的高速左旋走滑运动，表明东昆仑断裂带吸收了部分青藏高原的侧向挤出运动[82]；青藏高原北部阿尔金断裂带作为主控边界断裂，对中国大陆新构造运动影响显著，其挤压速率约为3 mm·年-1，平均走滑速率为6.9 mm·年-1，由西到东逐渐减小[11，25，80]；青藏高原东北缘的祁连—海原、六盘山断裂表现出中速左旋走滑特征，活动速率为3～5 mm·年-1，其中六盘山断裂具有一定的逆冲分量特征，而西秦岭活动速率较小[25，79]，青藏高原东北缘断层活动速率整体呈现自西向东减小趋势，也进一步表明了构造形变被其东端山脉隆升和逆冲断裂转换吸收。

3.2 深部地壳运动与形变特征

针对青藏高原深部地壳运动与形变机制，国内外学者采用大地测量监测数据结合地质、地球物理资料研究了印度板块北推挤压与欧亚板块阻挡作用下的青藏高原形变特性，并探讨了地壳有效黏滞系数与幂律指数的差异性对模拟结果的影响[2，58，83]。二维黏弹性有限元模型显示软弱下地壳流与壳幔运动的解耦作用导致了青藏高原东南缘的顺时针旋转特性与东部的伸展形变[84]；二维纵剖面黏弹性模型反演出的青藏高原深部黏滞系数特征表明，东部底层地壳内的流体以地表位移近8倍的速度围绕刚性的四川盆地和塔里木盆地流动[85]；基于GPS观测结果与震源机制解建立的三维黏弹性有限元模型显示出岩石圈横向于纵向流变特性，以及不同边界约束影响下的青藏高原伸展形变机制[9]；青藏高原近200 ka岩石圈运动形变过程也显示出，印度板块北推挤压作用对青藏高原隆升机制具有重要影响，但这不是导致青藏高原隆升的唯一机制，还可能包括深部地幔对流等因素[86]；青藏高原东南缘构造活动断裂分布于流变特性的三维有限元模型揭示出，当下地壳运动速率比地表GPS观测速度高10 mm·年-1时，数值模拟青藏高原东南缘地表位移与GPS观测结果会达到最佳拟合，青藏高原东南缘现今地壳形变模式主要受区域边界动力作用控制与主要活动断裂影响[87]；基于GPS速度场的柔性地壳动力学作用下，青藏高原整体抬升现象及绕喜马拉雅东构造结顺时针旋转特性显示当下地壳边界位移大于地表时，地表水平速度的均方根误差达到最小[88]；青藏高原下地壳不同流动方式与地表地貌关系的二维模型则揭示出地壳流变特性的横向非均匀性和地表构造特征的时序演化特征以及地壳黏性物质流动对青藏高原垂向形变的控制作用，表明下地壳管道流影响下的物质黏度变化会影响隆升速率[89]；热力学模型揭示出亚洲板块内部的软弱带在青藏高原形成过程中起作用，青藏高原尤其是其东南缘广泛存在的下地壳流也被数值模拟结果[90-93]所证实，而大地电磁成像结果也进一步表明青藏高原深部存在的下地壳流从青藏高原水平扩展800 km直至中国西南地区，并对青藏高原的隆升起到了促进作用[94]。

除下地壳流变作用影响青藏高原深部运动形变之外，基于震源机制解与地质调查资料的三维遗传有限元模型，研究青藏高原东南缘下地壳对上部地壳的拖拽作用对地表运动的影响，认为考虑下地壳对上地壳SSE向剪切作用的模型更为合理[95]；青藏高原东北缘岩石圈波速、流变结构与GPS观测结果的三维有限元模型模拟结果揭示了地幔对地壳对流拖拽力影响下的地块内部形变，更加明确了地幔拖拽作用对青藏高原运动形变的明显控制作用，并且不同地块壳幔耦合程度的差异性支持地震波各向异性这一特点[96]；地球物理探测结果进一步表明青藏高原内部温度高，中下地壳局部地区可能存在熔融，因此，在运动过程中下地壳更容易流动，从而对上地壳产生拖拽作用[97]。

4 现今区域构造动力学背景

4.1 构造动力学背景

中国大陆位于欧亚板块东南部，是一个晚第四纪和现代构造活动强烈的地区，其构造形变的动力源主要来自板块边界作用力：东部受西太平洋板块向西、菲律宾板块向西北的俯冲和削减联合作用，形成了一系列与弧后扩张有关的陆海缘伸展和断陷盆地，构造应力场的主体特征表现为NEE—SWW向挤压；西部和西南部主要受印度板块向北持续碰撞北推影响，形成了著名的“喜马拉雅—青藏高原”造山带，构造应力场的主体特征也表现为近SN—NNE向挤压(图10)[78，81，98]。

图件引自文献[81]图10 中国大陆及邻区主要活动断裂、活动地块分布Fig.10 Distribution of Main Active Faults and Blocks in Mainland China and Adjacent Areas

在中国大陆构造格局中，青藏高原处于十分特殊的构造位置，被誉为世界“第三极”，是现今大陆上最年轻、面积最大的地块，构造活动极为强烈且至今仍在继续。青藏高原块体内部可细划分为拉萨、羌塘、巴颜喀喇—松潘甘孜、昆仑—柴达木、祁连山以及川滇地块，每个地块均被狭长的活动构造所分割[99]。青藏高原周缘的活动构造主要有阿尔金活动构造带、祁连山—河西活动构造带、海原—六盘山活动构造带、龙门山活动构造带和安宁河—则木河—小江活动构造带。大多活动构造带表现出的逆冲或左旋逆冲几何形态很好地反映了印度板块与欧亚板块碰撞作用下青藏高原的受力状态。

印度板块与欧亚板块碰撞后，大陆板块之间的作用并没有停止，印度板块仍以44～50 mm·年-1的速率向北运动，青藏高原则受到由南向北运动挤压作用后向北扩展，而其北部又受到西伯利亚板块阻挡和塔里木地块向南的约束作用，致使青藏高原快速隆起。至少1 500 km的SN向缩短量被吸收，使青藏高原成为两倍于正常地壳的巨厚陆壳体，并在印度板块与西伯利亚板块之间形成了超大范围(SN向为2 000 km，EW向为3 000 km)的新生代陆内形变区域[100-101]。在青藏高原内部，早期由于沿几条EW向断裂带发生陆内俯冲，导致一系列推覆构造群和地壳隆升，后期则由于均衡调整，产生大幅度快速隆升[102]。Zhang等对青藏高原及其周边地区GPS观测资料的分析发现，印度板块与欧亚板块相对运移的70%～94%被青藏高原内部的构造形变所调整吸收[10]。在青藏高原北缘，阿尔金断裂带受到挤压，运动速率向北逐渐减小，被祁连山内部隆起和两侧新生代盆地形变引起的缩短所吸收，而在海原—祁连山断裂带则被断裂尾端的陇西盆地形变以及六盘山隆起所吸收[12，103-104]。

印度板块与欧亚板块碰撞的同时，受块体内部一系列大规模走滑断裂的影响，青藏高原应力场EW向呈现出明显的拉张作用[11，58]，主应力则为NE向[105]。块体内部物质在重力均衡作用下向东运移[15，106]，且不同地块间物质运动的方向和速率均存在一定差异：由南向北，北向地壳运动速率逐渐减小，运动方向从南部恒河平原约NE20°起逐渐向东偏转，东向运动速率逐渐增大，在青藏高原腹地玛尼—玉树区域达到最大，此时运动方向约为NE70°；再向北跨过昆仑和柴达木地块后，运动方向转为约NE60°，东向运动速率略有降低，总体上形成了以玛尼—玉树—鲜水河断裂带为中心的流动带[77，107]。

青藏高原东北缘由于受北部阿拉善地块和东部鄂尔多斯地块的阻挡，区域地壳在印度板块的挤压作用下不断缩短、抬升，印度板块与欧亚板块碰撞作用则由近SN向转换为NEE向，区域应力也由青藏高原西部的NNE向转变为NEE向，成为物质东流的重要汇聚带[108]。而较为稳定的扬子地块和鄂尔多斯地块则限制了青藏高原向东扩张，造成了青藏高原东缘抬升[85，109]。甘东南区域北部应力场呈NE—SW向挤压，而中部则转变为近NWW—SEE向挤压，这种应力场的变化直观地反映了鄂尔多斯地块对青藏高原向东扩张的阻挡作用；甘东南区域南部受扬子地块的阻挡，区域构造应力场则由近EW向挤压转变为NNW—SEE向挤压[110]。青藏高原向东运移的物质受到阻挡后，在东南缘发生顺时针旋转以SE向进入川滇地块南部后转为向南运动[80，111]，通过玉树—甘孜、鲜水河、安宁河—则木河—小江以及怒江、澜沧江等断裂活动吸收和调整[16，112-113]。同时，缅甸微板块向东俯冲导致川滇地块构造应力场形式较复杂，西北部川西地区压应力场主要表现为NNE—NE向，而川滇地块内部则主要表现为NW向，在块体南部边界附近则转为近NW向[114]。

4.2 区域动力学环境下珠峰地区构造隆升机制

地块隆升通常会伴随着构造特性的巨大改变。为进一步探讨青藏高原南部最为特殊的珠峰地区隆升机制和构造动力学特征，国内外学者在喜马拉雅地区开展了大量研究。1991～1995年，中美两国科学家首次将GPS技术应用到喜马拉雅地区板块形变监测，Bilham等分别对测量结果进行分析，得出喜马拉雅中段地壳现今汇聚速率为(17.5±2)mm·年-1[115]，西北段汇聚速率为14～17 mm·年-1[116-117]，东段汇聚速率为17～22 mm·年-1[117-118]。青藏高原是全球板块运动最为剧烈的地区之一，珠峰高程变化也一直是关注的焦点。继1975与2005年后，2020年12月8日中尼两国联合发布珠峰最新高程为8 848.86 m。相比于以往，此次珠峰高程测量的科学性更高、可靠性更强、创新性更高，此次获取的珠峰高程数据将为该地区的隆升机制研究提供重要的数据基础[119]。

国内学者基于青藏高原构造动力学背景，从多种不同角度分析了青藏高原区域构造动力学机制及其对喜马拉雅造山带的影响机制。傅容珊等以地壳短缩为基础,从地幔动力学角度论证了地块剥离隆升、挤压隆升到对流隆升达到目前的高度,且得出地块隆升的非稳性、多阶段、多机制驱动特征[120]；滕吉文等从力学角度提出了隆升、地壳短缩和增厚的动力学模式，论证了地块整体隆升的主要因素是印度板块与欧亚板块碰撞、挤压和长期楔入，并发现重力均衡作用和热作用虽在地块隆升和地壳短缩过程中发挥了重要作用，但不是主要动力源[121]；李廷栋认为地块隆升过程可总结为陆内汇聚、地壳分层加厚、重力均衡调整3种模式[102]；许志琴等认为地块隆升驱动力来自多种因素，印度板块俯冲、地块周缘克拉通作用及深部热驱动都是地块隆升不可忽视的地内因素，从而提出“周缘内向的陆内俯冲及腹地地幔底辟”隆升机制[122]；孔祥儒等通过分析西藏西部地球物理剖面与岩石圈结构特性，发现地块深部构造由具有不同特性的多个块体组合而成，认为青藏高原的形成具有多阶段、多形式、多机制与多块体的特性[123]；方小敏从青藏高原隆升的演化过程角度，将其划分为3个主要生长隆升期次，早期隆升阶段主要集中在地块中南部的拉萨地块和羌塘地块，中期在喜马拉雅和可可西里—昆仑山地区开始强烈隆升，晚期在喜马拉雅和地块东北部隆升显著加速[124]；葛伟鹏发现在地块内部和边缘地区的垂向运动存在很大差别，地块内部的地壳正在减薄，而边缘地区却存在不同程度的增厚与隆升[125]；Niu再次明确了俯冲板片的拖拽力是板块运动与喜马拉雅造山带隆升的主要驱动力[126]。

运用构造地质学、地震学对地壳运动进行研究仍是一种间接估算的方法，近年来在青藏高原开展的GPS高精度观测有力推动了青藏高原地壳运动研究的发展[127-130]。图2～4中得到的速度场能为青藏高原现今构造形变总体态势、长期构造演变历史等科学问题提供最直观并且有效的证据，也可在一定程度上为更好地解释珠峰地区隆升扩展机制，以及深入认知产生大陆碰撞带构造形变背后的大陆地球动力学背景提供必要且定量的参考信息。研究青藏高原地壳构造运动形变的核心是印度板块与欧亚板块碰撞动力学效应，因此，在广泛地质和地球物理观测的基础上，对该动力学过程进行数值模拟研究就显得尤为重要。一些学者利用数学物理方法对地壳形变分析和模拟进行了大量研究，但这也仅仅是研究的开始。对于现今青藏高原南缘特别是珠峰地区地壳形变来说，如何结合局部详细地质条件并且综合多技术手段深入分析断裂活动与块体运动形变特征，还需要科研工作者的不懈努力和探索。

5 结 语

5.1 结 论

本文系统总结了青藏高原GPS监测网络体系与高精度数据处理方法及策略，全面阐述了当前用于青藏高原地壳运动与形变分析的模型，结合3类模型计算结果从浅部与深部两方面深入分析了青藏高原现今地壳运动与形变特征，最后从板块运动构造动力学机制角度详细解译了青藏高原现今动力学环境及其影响下的珠峰隆升机制。GPS高精度地壳运动速度场定量表明了青藏高原现今地壳SN向挤压与EW向伸展运动的特征，主干断裂带走滑与挤压速率较高，深部地壳活动研究则进一步揭示了青藏高原下地壳流变特性以及地幔对地壳拖拽作用可能是导致青藏高原隆升的重要原因之一，而青藏高原与珠峰的隆升是大陆板块俯冲、区域重力加载与地壳深部地幔多尺度对流等多因素协同作用的结果。

5.2 展 望

(1)现阶段中国用于大陆地壳形变监测的“中国大陆构造环境监测网络”国家重大科技基础设施建设项目的建设完成与不断完善，可为青藏高原及其周边地区提供长期高精度形变监测数据。考虑到特殊需求，近年来国内外相关研究机构和高校科研院所也对青藏高原内部重点关注地区进行加密观测，丰富了区域地壳形变监测数据，为更精细获取青藏高原现今地壳运动信息提供了宝贵的基础资料。但受限于观测环境，青藏高原内部仍存在许多观测盲区，一些大型活动断层近场观测资料仍较缺乏，因此，今后需对青藏高原内部重点区域进一步加密台站。此外，融合多源观测资料获取更高精度、更高时空分辨率的监测信息，实现各种监测技术的优势互补，将成为未来青藏高原地壳形变监测研究的重要发展趋势。

(2)青藏高原下地壳与地表运动的解耦效应为数值模拟工作中深部位移约束带来了困难，因此，研究岩石圈深部运动特性仍是未来青藏高原深部地壳构造运动与形变研究的重点。此外，若采用近20年的GPS观测结果作为位移约束模拟青藏高原数十万年甚至百万年的构造形变特征，大地测量与地质学时间尺度的不对等问题也需进一步深入探讨。

(3)青藏高原处在全球构造运动最为活跃的地区之一，受多种驱动力源作用，板块边界构造机制复杂。虽然学者们对青藏高原形成演化机制与构造动力学特征进行了很多年研究，但目前仍存在许多争论。例如，青藏高原北边界活动断裂对地块隆升的动力学约束作用至今仍存在争议；青藏高原周缘不同地形特征的差异是受何种结构控制？此外，虽然相关研究也证实了青藏高原物质东流受阻挡后转向SE向，但物质的流动方式和流动深度尚未准确确定；青藏高原东北缘复杂的横向扩展模式详细机制等也欠缺探索。这一系列科学问题仍需今后深入研究。

谨以此文庆祝长安大学七十周年华诞。七十年风雨兼程，七十年青春如歌；年华流转，不变的是您永驻的学者心；岁月如流，永恒的是您铸造的师者魂！2001年8月，我从一个毗邻边境仅仅二十多公里的新疆边陲小镇吉木乃，拖着沉重行李箱，带着当时全县唯一一个高三理科班考取疆外重点大学的应届生“光环”，经历了短途汽车、长途大巴、绿皮火车近四天四夜的时间辗转来到了母校。至今还清晰记得2001年8月25日清晨天蒙蒙亮的时候,和父母一起站在雁塔校区东门口敲击传达室窗户请求开门的场景，当走进长安大学的一刹那，心中是莫名的感动，周围的一切是那么的新鲜，也是那么的让人憧憬！当偶然听到周围有人抱怨学校操场地面灰尘太大时，我淡然一笑，心中却已开始遐想自己驰骋沙场的英姿！当偶然听到周围有人抱怨教室设施不够现代化时，我会心一笑，心中却已开始憧憬自己坐在明亮教室里沐浴着师恩的教诲！当偶然听到周围有同学抱怨专业的辛苦而选择离开时，我微微一笑，心中却已开始期待若干年后自己的学业有成！对于母校，我内心只有感恩，感恩母校给予了我与我心中“女神”一场如此美妙的相遇！这一场相遇从大学到工作至今已有约二十年时间。时至今日，也难忘开学典礼时，1号食堂大礼堂墙壁上的红色醒目大字“今天我以长安大学为荣，明天长安大学以我为傲”，这时刻激励着我努力学习的热血。到现在一路在专业上的坚持以及与母校七千多个日夜的相守，“弘毅明德，笃学创新”的母校精神已深深融入我的生命，融入我的心灵，亦融入我的灵魂！