距今约5500万年以来，青藏高原在印度板块的挤压碰撞下发生形变、隆升，直至形成今天全球平均海拔最高的高原(Tapponnier et al., 2001).这一形变过程的物理机制是地学界长期争论的重要科学问题.其中一种颇具代表性的模型是黏性下地壳流模型(如：England and Molnar, 1997; Clark et al., 2005; Royden et al., 2008).黏性下地壳流模型认为青藏高原岩石圈下地壳为广泛分布的黏性流体，由于流体的不可压缩性，当下地壳流受到西南印度板块的推挤和内陆刚性块体的阻挡时，介质便发生挤压方向的缩短和垂向增厚，并驱动上地壳形变和隆升.由于下地壳流的广泛分布，高原地表形变也呈连续分布状态.这一模型成功地解释了青藏高原高海拔的形成，但无法很好地解释青藏高原沿其北缘大型走滑断裂带的东向运动.

与之相反的另一种模型为块体运动模型(如：Tapponnier et al., 1982; Avouac and Tapponnier, 1993).这一模型认为青藏高原上、下地壳耦合构成近刚性的块体，相邻块体的相互作用驱动地表产生弹性形变，且形变主要集中在块体边界.当印度板块向青藏高原推挤时，其形变主要由青藏高原地壳的东向逃逸和高原内部跨逆冲断裂的全地壳挤压缩短来吸收.因此，该模型可以很好地解释青藏高原沿其北缘大型走滑断裂带的东向运动，但难以解释高原内部广泛分布的形变场和小规模断裂.

上述两种模型可以认为是描述青藏高原地壳形变模式的两个极端模型，真实的地壳形变场应当介于两者之间(王敏等, 2003; Zhang et al., 2004).

决定地壳形变机制的关键在于地壳与岩石圈介质的力学结构，特别是下地壳与上地幔的流变学结构.由于缺乏真实的野外观测，人们对下地壳与上地幔的流变学结构所知甚少，这也是有关青藏高原地壳形变模式长期争论未决的主要原因.迄今为止，人们对岩石圈流变学结构的认识大多来源于地震学与实验室岩石力学实验的间接推测(如：Kirby and Kronenberg, 1987)，结果并不可靠.近些年来卫星大地测量在地壳形变监测领域的广泛应用为相关研究提供了另一种可能：地震造成区域应力场的改变，进而触发下地壳和上地幔介质的黏弹性松弛引起地表形变; 运用地球物理模型模拟震后形变场，便可反演该地区岩石圈的流变学结构.例如，Rousset等(2012)模拟GPS观测到的1999年台湾集集地震的震后形变场，推断中央山脉下方介质存在一个低黏度(黏滞系数为5×10¹⁷Pa·s)的中地壳层；Freed等(2006)模拟2002年美国阿拉斯加Denali地震的GPS震后形变观测，得到该地区60 km以下可能存在低黏度(黏滞系数约为3×10¹⁸~4×10¹⁸Pa·s)上地幔层的结论.

2001年11月14日青藏高原北部东昆仑断裂带上发生的M_W7.8(www.seismology.harvard.edu/ CMTsearch.html)/M_S8.1 (中国地震台网)可可西里地震是该断裂带过去100年来记录到的最大地震，造成藏北地区岩石圈应力场的很大改变，为研究青藏高原北部流变学结构和形变机制提供了难得的机会.本文通过模拟震后形变研究青藏高原北部地区岩石圈流变学结构，并据此判断下地壳流模型在该地区的适用性.

前人关于青藏高原北部岩石圈介质结构和物理性质的许多研究都表明跨昆仑山地区岩石圈介质结构和物理性质存在显著的横向不均一性.例如，Karplus等(2011)采用人工地震源数据反演得到昆仑山以南的羌塘-巴颜喀拉块体地壳厚度大于其以北的柴达木盆地地壳厚度，前者厚达近70 km，后者则约50 km；Sun等(2010)的剪切波速度成像结果显示巴颜喀拉-羌塘块体下地壳S波速度显著低于柴达木盆地下地壳S波速度；Unsworth等(2004)的大地电磁反演结果也表明巴颜喀拉-羌塘块体下地壳的电阻抗明显小于柴达木盆地下地壳电阻抗.因此，我们的模型考虑跨昆仑山两侧介质的横向差异，试图通过大地测量的手段认识和验证青藏高原北部地区地壳介质结构的复杂性及其对青藏高原形变场的影响，增进我们对于青藏高原地壳与岩石圈形变模式的认识.

我们知道，震后形变的主要机制包括下地壳与上地幔介质的黏弹性松弛和断层面上的震后余滑.在断层近场和震后短时间周期，震后形变场主要来源于断层面的震后余滑；而在远场及震后长时间周期，介质的黏弹性松弛对震后形变场的贡献更为显著.但这两种形变机制所引起的地表形变都表现为类似的随时间呈对数或指数函数衰减的特征，仅在时间域内难以从观测中将两者区分(Savage, 1990).因此，为简单起见，许多震后形变的模拟研究往往只考虑一种震后形变机制，比如只考虑介质黏弹性松弛(Ryder et al., 2011; Zhang et al., 2007)，或震后余滑(Ryder et al., 2011).本文则同时考虑黏弹性松弛和震后余滑两种机制的贡献，建立联合反演模型.反演结果可对东昆仑断裂两侧下地壳-上地幔黏弹性介质的流变性质提供迄今为止最为精确的约束.

2001年可可西里地震发生后，中国地震局展开了地震综合科学考察，其中包括沿青藏公路布设与观测的由4个临时固定站和12个流动点构成的跨断层GPS剖面和在震区较大范围内开展的GPS流动点重复观测(任金卫和王敏，2005).地震发生后的第二年我们沿断层两侧又布设了30个GPS流动观测点，使震后形变观测点的总数增加至60个.基于流动点复测期次的考虑，本研究共选用其中45个点；且考虑到本地区震后形变场可能受到2008年3月21日M_W7.3新疆于田地震和2008年5月12日M_W7.9四川汶川地震的影响，仅使用于田地震之前的观测数据分析和表征可可西里地震的震后形变场，数据的时间跨度最大为6.4年.

如图 7所示，GPS观测点主要集中在地震破裂带的近场，但中远场也有一定的分布.其中，约35%的点从震后数天至3个月内开始观测，到于田地震前有7~10期观测，多位于中场和远场；约65%的点观测始于震后3个月，在4.5~6.4年的时间内有4~5期观测，多数为近场点.此外，震后布设的临时固定站(BDGD、KLGD、JB30)维持了近1年的时间，有200~300天连续观测，此后JB30改为流动观测.DLHA是中国地壳运动观测网络的基准站，数据有很好的连续性.GPS观测资料在时间与空间的分布为区分介质黏弹性松弛和震后余滑两种震后形变机制提供了可能.

图 7

Fig. 7

图 7 GPS观测点位分布与震后1年地表位移拟合 橘黄色实线表示同震破裂断层地表迹线，橘黄色虚线为盆地南缘位置；蓝色三角表示GPS观测点位置，蓝色和红色箭头分别表示震后1年地表位移的观测值和模型预测值. Fig. 7 Distribution of GPS sites and data fitting of first year observed and predicted postseismic displacements The orange solid curve indicates earthquake surface rupture trace and the orange dash line denotes the southern boundary of the Qaidam Basin. Blue triangles represent GPS sites, and blue and red arrows are first year observed and predicted postseismic displacements, respectively.

由于GPS观测点绝大部分是震后布设的，为了区分地壳长期构造形变和可可西里地震造成的震后形变，我们基于本地区已有的长期构造运动速度场(图 1，王敏等，2003)插值得到GPS点的震前构造运动速度(插值方法见Shen et al., 2015)，进而获得这些点的震后位移时间序列.进一步，可可西里地震以走滑破裂为主，造成的地表垂直形变很微弱，而GPS垂向观测精度、特别是流动点的垂向观测精度不足以监测到这样微弱的形变.因此，本研究只采用水平向观测量，即东向和北向两个分量的数据为模型提供约束.

图 1

Fig. 1

图 1 青藏高原北部区域构造、地震活动和GPS震间水平运动速度场 蓝色箭头为相对于稳定欧亚板块的速度场(王敏等，2003)，其误差椭圆的置信水平为90%；过去40年MS≥6.0地震的震源机制来自http://www.globalcmt.org/CMTsearch.html. Fig. 1 Regional tectonics, earthquakes and GPS derived interseismic horizontal velocity field in northern Tibetan Plateau Blue vectors are GPS velocities relative to the stable Eurasia plate (Wang et al., 2003) with confidence level of the error ellipses 90%. Focal mechanism solutions of MS≥6.0 events in the past 40 years are from http://www.globalcmt.org/CMTsearch.html.

大地震的发生造成岩石圈应力场的改变，从而引起下地壳与上地幔介质的黏弹性松弛和断层面的震后余滑.计算下地壳与上地幔介质的黏弹性松弛，需要选定一个同震破裂模型作为驱动源，同时，反演震后余滑也需要先验的断层几何形态.

2001年可可西里地震破裂过程相当复杂，因此，前人基于不同方法获得的认识也不尽相同.例如，Lin等(2002)的野外地质考察结果显示断层最大同震破裂为~17 m，位于库赛湖附近.Xu等(2002)的独立考察结果与Lin等(2002)的大致相似，但最大破裂确定为~6 m，并发现太阳湖段有拉张分量.Klinger等(2005)基于高分辨率卫星光学影像数据得到最大走滑位错达9 m，且发现库赛湖段有逆冲分量.Lasserre等(2005)根据野外地质考察和高分辨率卫星光学影像测量结果假定断层为纯直立断层，并根据InSAR观测数据反演得到同震破裂分布；反演结果显示有3处破裂极大值，且最大走滑位错为~8 m.万永革等(2008)采用GPS和InSAR数据联合反演位错分布，得到位错最大值为~7 m，同时还反演得到库赛湖段和昆仑山口段的倾角分别为~81°和~80°，与库赛湖段存在逆冲位错(Klinger et al., 2005)的研究结果相符；基于Xu等(2002)的野外考察结果将太阳湖段的倾角设定为90°~45°的渐变，反演得到明显的拉张分量(图 2).

图 2

Fig. 2

图 2 可可西里地震的同震破裂模型(万永革等，2008) Fig. 2 Coseismic rupture model of the Kokoxili earthquake (Wan et al., 2008)

对比Lasserre等(2005)和万永革等(2008)两个基于大地测量数据反演得到的同震破裂模型，由于采用GPS数据和InSAR数据的联合约束，并同时反演断层倾角，万永革等(2008)的模型具有更好的解析度和可靠性，与野外地质观测符合更好.因此，本研究采取万永革等(2008)的同震破裂模型作为下地壳和上地幔介质黏弹性松弛的驱动源和震后余滑的先验断层面.其中，断层分为主断层段和西端太阳湖段，主断层段倾角80°~81°，是同震破裂的主要段落；西端太阳湖段是支断层，显著南倾.主断层段和太阳湖段只在地表贯通，分隔点在91.1562°E处.

震后GPS观测点主要位于主断层段两侧，太阳湖段几乎没有观测点(图 7)，因此太阳湖段的震后余滑分布难以得到有效约束；且绝大多数GPS点至太阳湖段的距离远大于太阳湖段自身尺度，太阳湖段的震后余滑对观测点震后形变的贡献由于距离遥远可忽略不计.所以本研究只反演同震破裂主断层段的震后余滑.由于震后余滑可能扩展到同震破裂区域之外，又将主断层段继续向东延伸40 km，并向深部延伸到50 km(图 4).

图 4

Fig. 4

图 4 (a) 同震库仑应力变化；(b)同震库仑应力变化与震间累积应力之和 注意二图中库仑应力色标范围的不同. Fig. 4 Coseismic Coulomb stress change (a) and sum of Coulomb stress change and interseismic accumulated stress (b) Note the color bars of the two figures are different in denoting the level of Coulomb stress.

岩石圈由上、下地壳和上地幔组成.一般认为上地壳主要呈现弹性特征，用弹性体模拟.由于高温高压，下地壳和上地幔的介质行为在长时间尺度内表现为黏性而非弹性.当地壳受到同震应力扰动产生弹性形变时，黏弹性的下地壳和上地幔在同震应力扰动作用下持续形变，即黏弹性松弛.

在青藏高原北部地区，岩石圈的弹性介质结构和物理性质存在显著的横向不均一性(Karplus et al., 2011; Sun et al., 2010; Unsworth et al., 2004)，预示着其流变性质也可能存在横向不均一性.进一步，研究区域的介质分界面与断层位置也不重合，即介质横向差异的分界面位于昆仑山北缘、柴达木盆地南缘(Karplus et al., 2011)，而破裂断层则位于该分界线以南约80 km处(图 3和图 7).位于断层北侧的25个观测点中有多达12个点位于昆仑山区域，只有10个点位于盆地内，另有3个位于交界地带(图 7).这意味着断层以北多数观测点的下方介质更接近断层南侧介质的结构和性质.

图 3

Fig. 3

图 3 研究区域介质结构示意图 Fig. 3 Sketch of medium structure in the study area

本文采用PSGRN/PSCMP程序(Wang et al., 2006)计算黏弹性松弛引起的震后形变和反演震后余滑时所需要的弹性位错格林函数矩阵.PSGRN/PSCMP要求一个平面成层的介质结构，当断层两侧介质横向不均一时，通常采取近似方法，即对于坐落在某介质之上的台站，采用相应介质参数代入模型进行计算.因此，本研究将断层近似视为介质分界面，假设断层两侧的介质具有独立的黏滞系数，但断层北侧介质的弹性结构和性质都等同于断层南侧.基于这样的假设建立的流变学结构模型称为模型A.

为了分析介质结构和性质对地表数据模拟的影响，我们又建立了两个对照模型，分别称作模型B和模型C.其中，模型B完全不考虑介质的横向差异，断层两侧均采用南侧介质的结构和性质；模型C则同时考虑介质弹性和黏弹性性质的横向差异，但忽略断层面与介质间断面的不一致性.

模型中，地壳介质厚度和P波、S波速度结构均参考Karplus等(2011)和Mechie等(2012)根据InDepth Ⅳ项目人工反射地震波数据分析得到的结果.其中，上地壳为弹性，用弹性体模拟，根据余震定位的深度结果(Wei et al., 2010)和同震破裂的深度范围(Lasserre et al., 2005; 万永革等, 2008)，设置其厚度为20 km.下地壳和上地幔为黏弹性，由于GPS数据覆盖震后的短期与长期形变场，因此采用Burgers体模拟，并假定下地壳和上地幔黏滞系数相同.Burgers体由Maxwell体(刚性模量k₁的弹簧和黏滞系数η₁的黏壶串联)和Kelvin体(刚性模量k₂的弹簧和黏滞系数η₂的黏壶并联)串联组成(Bürgmann and Dresen, 2008)，有效刚性模量k₁由介质的S波速度和密度求得，未完全松弛的刚性模量k₂=k₁α/(1-α)，根据Ryder等(2011)的研究设α=0.69；η₁和η₂分别表示介质长期和短期黏滞系数.因此，模型A和C均有4个未知参数待定，即断层南北两侧的短期黏滞系数η_2s、η_2n和长期黏滞系数η_1s、η_1n；模型B则有2个未知参数待定，即研究区域介质平均的短期黏滞系数η₂和长期黏滞系数η₁.

大地震产生的应力扰动不仅作用于区域介质，而且作用于断层面，因而产生的形变也与断层面的力学性质相关.一般认为，断层在脆性破裂层以下呈现速率强化的特征，即摩擦系数与错动速率有关，断层错动速率越大，摩擦系数越大，滑动的阻力越大(Scholz, 1998).在同震应力扰动的驱动下，断层面震后会呈现出快速的驰豫滑移，且一般主要集中在断层的脆韧性转换带上.

震后余滑在地表产生的位移场随时间呈对数函数衰减的特征，可以采用(1)式拟合：

(1)

其中，地震发生时刻t=0，C_ij和D_ij(t)分别表示观测点i震后初始观测时刻和t时刻的震后位移j分量；A_ij和τ分别为震后位移随时间衰减的对数函数振幅和时间特征常数，其物理意义为：在t=τ(e-1)时，观测点i的震后位移j分量为A_ij.C_ij和A_ij由最小二乘法通过数据反演得到，而最佳τ值则由网格搜索法获得.

由地表在t时刻的位移D_ij(t)反演相应时刻断层面上的滑移分布.根据弹性位错理论，观测方程可写为：

(2)

其中，G为弹性位错格林函数矩阵，S为待反演的震后余滑参数向量，由于地震以左旋破裂为主，故假定只存在走滑分量；d₁为由震后t时刻的位移D_ij(t)组成的准观测向量，本文选定t=1年，ε₁为其误差向量，E₁为ε₁的协方差矩阵.

为了弥补观测点空间分布的不足，更好地区分两种震后形变机制造成的地表形变场，利用先验条件约束震后滑移是个行之有效的办法.如前所述，断层面的震后余滑与断层面的力学性质相关，震间构造运动在断层面上积累了应力，这部分应力一部分随同震释放，所余部分驱动断层面发生余滑.因此，除了采用运动学约束(式(2))外，本研究还采用了库仑应力的先验约束，即假定断层面上各单元的滑动量近似正比于该单元的库仑应力，而库仑应力为震间累积应力与同震库仑应力变化之和.类似的应力约束模型也被前人用于其他地震的震后形变研究中，例如Johnson等(2009)假定震后滑移与库仑应力成正比，其比例为速率-状态相依本构关系的经验常数a-b，并由此反演了2002年美国阿拉斯加Denali地震的震后余滑分布.

假定震间断层面上积累的构造应力呈均匀分布，其中脆性层内(0~20 km深度)为一常数；转换带20 km至35 km随深度呈线性衰减至0；35 km深度以下为0.进一步假定震间断层上积累的构造应力在最大同震破裂处完全释放，等同于同震库仑应力释放的最大值即11 MPa(由万永革等(2008)同震破裂模型计算得到)，但方向相反；断层面上总应力分布为同震库仑应力变化和震间累积应力之和(图 4).

由图 4可以看到同震库仑应力变化的分布存在两条较为明显的分界线，一条是在东经93.5°—94°之间，另一条在20 km深度线上.这是由同震滑移的分布和断层面的延展造成的.对于一个断层单元而言，同震库仑应力变化源于本身滑移造成的卸载和周边单元滑移贡献的加载，因此，同震库仑应力变化的量值与同震滑移的空间梯度相关，即在同震滑移不均一处，高滑移量一侧产生相对卸载，而低滑移量一侧为相对加载；滑移的空间梯度越高库仑应力变化的程度就越大.

因此，应力约束的观测方程为：

(3)

其中，F为总库仑应力向量，k为一待定常数，d₂为零向量，表示震后余滑与总库仑应力近似成正比，ε₂为d₂的误差向量，其协方差矩阵为E₂，其主对角元大小可由解的分辨率与数据拟合残差间的折衷曲线获得.

(2)、(3)两式合并后写为：

(4)

令，由d₂= 0，则式(4)的最小二乘解为：

(5)

解的分辨率矩阵R可表示为(Jackson and Matsu′ura, 1985)：

(6)

其中，各断层单元滑动量的分辨率为其所对应的R的对角元素，值在0~1之间.若对应的对角元素为零，则表示观测数据对该参数的分辨率为零，该参数的确定完全来源于库仑应力的先验约束；反之，若对角元素为1，则该参数完全由观测数据确定.

本文尝试利用黏弹性松弛和震后余滑两种机制联合解释震后的地表位移.相对于地表位移，断层面的余滑是线性参量，而介质的黏滞系数是非线性参量，因此，我们采用正演和反演相结合的方法.具体的步骤如下：

首先给定一组介质黏滞系数(例如本文流变学模型A有2个短期和2个长期黏滞系数)，采用PSGRN/PSCMP程序(Wang et al., 2006)正演计算黏弹性松弛引起的观测点的位移时间序列.

其次，从各观测点的观测时间序列中减去黏弹性松弛引起的位移时间序列，得到剩余位移时间序列，并用对数函数(公式(1))拟合得到时间特征常数和振幅, 由此推算震后1年的位移；然后采用库仑应力先验约束模型(公式(4))反演震后余滑的分布.

最后，计算模型时间序列(黏弹性松弛导致的位移时间序列+震后余滑导致的位移时间序列)对观测时间序列的拟合残差RMS值.各黏滞系数采用网格式搜索，其中，短期黏滞系数搜索范围设为10¹⁷~ 10²⁰Pa·s，长期黏滞系数搜索范围为10¹⁸~10²¹Pa·s，取拟合残差RMS值最小所对应的黏滞系数组为最优解.

利用模型A反演得到的区域介质黏滞系数最优解如图 5所示：东昆仑断层北侧柴达木盆地地区下地壳-上地幔短期和长期黏滞系数分别为5×10¹⁸ Pa·s和1.5×10²⁰ Pa·s，东昆仑断层南侧巴颜喀拉-羌塘地区下地壳-上地幔短期和长期黏滞系数分别为1.5 ×10¹⁸ Pa·s和1.5×10¹⁹ Pa·s.我们通过F检验方法(Shen et al., 1996)确定各黏滞系数的置信区间，得到柴达木盆地地区下地壳-上地幔短期和长期黏滞系数的90%置信区间分别为4×10¹⁸~9×10¹⁸Pa·s和＞6×10¹⁹Pa·s；而巴颜喀拉-羌塘地区下地壳-上地幔短期和长期黏滞性系数的90%置信区间分别为8×10¹⁷~2×10¹⁸Pa·s和1×10¹⁹~2×10¹⁹Pa·s.可见，观测数据对断层南侧介质的黏滞系数有较好的约束，但对断层北侧介质的长期黏滞系数只能分辨其90%置信区间的下限，而难以确定其上限.

图 5

Fig. 5

图 5 数据残差随黏滞系数变化的分布图 (a)断层南侧; (b)断层北侧.红色圆点表示最佳解位置. Fig. 5 Data residual chi-square distribution with respect to viscosity parameter realizations Panels (a, b) show the results of regions south and north of the fault respectively. The red dot denotes the location of best solution.

上述结果表明柴达木盆地下方介质的流变性质显著弱于巴颜喀拉-羌塘块体.并且，无论是柴达木盆地还是巴颜喀拉-羌塘块体，其下地壳-上地幔的短期黏滞系数都大大小于长期黏滞系数，这表明震后黏弹性松弛使得区域应变率提高，介质的有效黏滞系数大大减小.

与介质黏滞系数最优解相对应的震后余滑分布如图 6所示.从图 6a可以看到震后余滑在20~25 km深度即断层的脆韧性转换带上有集中的分布，这和Scholz(1998)提出的震后余滑机制相吻合；而在断层20 km深度内，震后余滑分布(图 6a)和同震位错分布(图 2)有一定的互补性，这与在反演震后余滑时采取了库仑应力的先验约束有关，也符合前人对震后余滑分布特征的普遍认识.此外，在93°E—93.5°E之间0~5 km处出现了反向滑移.对比图 2所示的同震破裂模型，可以看到此处正是同震位错的极大值所在位置，可能意味着同震位错出现了过冲.断层深部绝大多数区域反演得到的滑移量非常小，仅为几个毫米，完全在数值误差范围内，并不能说明相应区域存在余滑.余滑参数的分辨率显示，数据的约束作用主要集中在断层面东部的近地表区域(图 6b)，这是由GPS观测点多数位于近场及断层偏东一侧决定的；断层面大多数区域的震后余滑分布，尤其是断层深部的震后余滑分布，主要依赖于库仑应力的先验约束.

图 6

Fig. 6

图 6 震后余滑分布(a)和分辨率(b) Fig. 6 Inverted afterslip distribution (a) and resolution (b)

由于震后观测位移也不可避免地包括了余震的贡献，因此，反演得到的部分震后余滑可能反映了该区域余震的贡献.如果相应区域有余滑但并没有发生余震，则意味着该区域为速率强化区，具有一定的流变性质，体现了断层面上力学性质的复杂性，也意味着该区域不是孕震的危险区，而会通过无震蠕滑释放积累的能量.

图 7显示了两种震后形变机制联合反演模型对震后1年地表位移观测数据的拟合情况.可以看到，在断层以北、盆地以南的区域，特别是NK02、NK03、NK04和NK05几个观测点，观测值普遍大于模型预测值.这可能是由于这一区域下方介质的流变学性质更接近断层南侧，即黏滞性更低，但在模型中设定这一区域下方介质的流变学性质与盆地内一致，因而是高估了这一地区介质的黏滞系数.在断层南侧远场，也可看到观测值显著大于模型预测值，且越往南，差异越明显.这一方面可能由于震前速度场在此处存在较大误差，另一方面也可能意味着南侧羌塘地区介质的流变学性质比北侧巴颜喀拉地区更强，进一步反映了藏北地区介质流变学性质的横向不均一性.

图 8给出了模型对6个观测点(位置参见图 7，分别来自断层南北两侧近场和远场)观测数据的拟合情况，从中可以看到两种形变机制对不同点震后位移的贡献是不同的，即近断层点的位移主要来源于断层面震后余滑(EK06, BDGD)，而远场点的位移主要来源于下地壳-上地幔介质的黏弹性松弛(DLHA，JB49)，符合不同机制在时间域和空间域的作用模式，表明联合反演的结果具有很好的合理性.

图 8

Fig. 8

图 8 GPS观测点震后位移时间序列(相对于第一个观测时刻) 观测(蓝)、介质的黏弹性松弛(黑)和模型预测(红), 观测点位置见图 7. Fig. 8 Postseismic displacement time series of GPS sites(relative to the first observation time) Observed (blue) and predicted time series based on the viscoelastic relaxation model (black) and the joint model of two mechanisms (red). The locations of GPS sites are in Fig.7.

如前所述，为了分析介质结构和性质对地表数据模拟的影响，除模型A外，本研究还建立了两个对照模型B和C.其中，模型B不考虑介质性质的横向差异，采用均一的分层结构；模型C则同时考虑介质弹性和黏弹性性质的横向差异，并采取不同的分层结构.

采用与模型A同样的反演方法和步骤，模型B的黏滞系数最优解为：下地壳-上地幔短期和长期黏滞系数分别为1×10¹⁸ Pa·s和4×10¹⁹ Pa·s；模型C的黏滞系数最优解为：东昆仑断层北侧柴达木盆地地区下地壳-上地幔短期和长期黏滞系数分别为5×10¹⁸ Pa·s和1.5×10²⁰ Pa·s，东昆仑断层南侧巴颜喀拉-羌塘地区下地壳-上地幔短期和长期黏滞系数分别为1.5×10¹⁸ Pa·s和1.5×10¹⁹ Pa·s.模型A、B和C数据拟合的加权残差平方和分别为614、770和615.

比较三种不同的流变学结构模型所对应的最优解和数据的拟合残差，可以发现模型A和C的拟合残差明显小于模型B的拟合残差，这意味着柴达木盆地和巴颜喀拉-羌塘地区之间流变学性质的横向差异不可忽略.模型B的结果反映的是东昆仑断层两侧介质流变学性质的平均状态，因此其黏滞系数介于模型A和C所得到的断层两侧介质黏滞系数之间，即小于断层北侧介质的黏滞系数，是对断层北侧介质流变学性质的高估，同时大于断层南侧介质的黏滞系数，是对断层南侧流变学性质的低估.这种系统偏差导致介质黏弹性松弛引起的地表位移部分残存于震后余滑反演的准观测值中，但并不能被震后余滑模型所解释，因此，震后余滑模型对准观测值的拟合在断层两侧也呈现出可辨识的系统偏差，即在断层以南准观测值均大于模型预测值，而在断层以北准观测值则小于模型预测值.

模型A和C所得到最优解相同，数据拟合残差也基本一致，这可以从数据分布和模型设置得到解释.模型A和模型C的差异在于断层以北地壳厚度的不同(模型A断层以北地壳厚度与断层南侧一致均为68 km，模型C则为50 km)和介质弹性参数的不同.由于本研究假定下地壳和上地幔具有相同的黏滞系数，因此，地壳厚度的差异归根结底也表现为介质弹性参数的差异.弹性介质的横向差异会造成弹性位错的地表形变场存在横向差异，如陶玮等(2007)和万永革等(2008)发现GPS和InSAR观测得到的可可西里地震的同震形变场，断层南侧比断层北侧大10%~30%，但并不会对介质的黏弹性松弛效应产生影响.理论上讲，震后余滑造成的地表形变场也应反映出弹性介质的横向差异，但由于震后余滑的量级很小，现有的数据分布(断层以北近半数观测点位于真实的介质分界面以南)和精度尚无法对其有效识别.所以，我们认为本研究的数据可以区分断层南北二侧介质黏性性质的差异，但无法区分弹性性质的差异，从数据拟合的角度看，模型A和C并无本质上的差别.但考虑到所用数据的分布，模型A的弹性分层和参数设定是最为合理的.

2001年可可西里地震造成青藏高原北部地区大范围的震后形变，为模拟研究该地区的流变学性质提供了机遇，一些研究成果陆续发表(如：Zhang et al., 2007; Diao et al., 2011; Ryder et al., 2011; Wen et al., 2012).但是，由于数据在时间和空间上的分布不同，考虑的形变机制也不同，这些研究(包括本研究)所获得的区域介质的流变学参数并不具有直接的可比性.

例如Ryder等(2011)在假定断层两侧介质横向性质均一的条件下，以区域介质的黏弹性松弛模拟解释震后第1年的GPS观测数据和震后2~5年的InSAR数据，得到青藏高原北部下地壳的长期黏滞系数为1×10¹⁹Pa·s，小于本研究中同样没有考虑介质横向差异所得到的结果4×10¹⁹Pa·s(模型B)，这一差别应该与Ryder等(2011)忽略震后余滑效应有关.

Diao等(2011)同样假定介质横向性质均一，虽然同时考虑了黏弹性松弛和震后余滑两种机制，但由于只使用了震后4个月的数据，反演得到的黏滞系数可能只反映介质震后短期的流变性质，其下地壳黏滞系数~10¹⁸ Pa·s的结果与本研究所获得的断层南侧介质短期黏滞系数1.5×10¹⁸Pa·s (模型A和C)基本一致也说明了这一点.

Wen等(2012)也采用了震后余滑和黏弹性松弛两种机制解释震后形变场，其所主张的下地壳和上地幔黏滞系数2×10¹⁹Pa·s的结果与本文同样不考虑介质横向差异所得到的结果4×10¹⁹Pa·s(模型B)量级相当；同时，源于InSAR数据比GPS数据具有更高的空间密度，其在二阶平滑约束下反演得到的震后余滑分布也具有一定的合理性.但是，由于Wen等(2012)是先单独采用介质的黏弹性松弛解释形变数据，这一机制不能解释的部分即数据残差再用断层面的余滑解释，因此造成部分余滑导致的地表形变可能被介质的黏弹性松弛所解释，在理论上会低估介质的黏滞系数，这也可能就是其结果略小于本文结果(模型B)的原因之一.

综合分析前人及本研究的结果(参见表 1)可以发现：忽略震后余滑对震后形变场的贡献会低估介质的黏滞系数，而且，数据覆盖的震后时间越短低估的程度越严重；如果不考虑断层两侧介质的横向性质差异就不能获得对青藏高原北部地区流变学性质的精确认知，该地区的形变模式与特征也很难得到解释；震后形变观测数据在时间上和空间上的分布是区分黏弹性松弛和震后余滑两种机制的关键所在，因此也成为模拟研究青藏高原北部地区流变学性质的必要条件.

表 1 Table 1

表 1 基于可可西里地震震后形变的区域流变学参数估计 Table 1 Viscosity estimates of northern Tibet based on postseismic deformation of the Kokoxili earthquake 模型 震后机制 黏滞系数(Pa·s) 备注(模型设置与数据) 断层北侧 断层南侧 Zhang et al., 2007 黏弹性松弛 η1=2×1017 黏弹性层20 km以下；4个GPS连续站震后1年数据 Diao et al., 2011 黏弹性松弛 η1=(5~6)×1017 黏弹性层32 km以下；沿青藏公路GPS流动站震后4个月数据 震后余滑+黏弹性松弛 η1=1×1018 Ryder et al., 2011 黏弹性松弛 η2= 9×1017 η1=1×1019 黏弹性层16 km以下；沿青藏公路GPS流动站震后1年数据和震后2~5年InSAR数据 Wen et al., 2012 黏弹性松弛 η1=(2~5)×1019 黏弹性层深度15 km以下；震后2~6年InSAR数据 本研究模型A 震后余滑+黏弹性松弛 η2=5×1018 η1=1.5×1020 η2=1.5×1018 η1=1.5×1019 忽略介质分界面，断层北侧的介质分层和弹性参数都与南侧一致，黏弹性层20 km以下 本研究模型B 震后余滑+黏弹性松弛 η2=1×1018 η1=4×1019 同上 本研究模型C 震后余滑+黏弹性松弛 η2= 5×1018 η1=1.5×1020 η2=1.5×1018 η1=1.5×1019 同时顾及断层两侧弹性介质的横向差异，但将断层视为介质分界，黏弹性层20 km以下

表 1 基于可可西里地震震后形变的区域流变学参数估计 Table 1 Viscosity estimates of northern Tibet based on postseismic deformation of the Kokoxili earthquake

相比于前人的研究，我们的研究所采用的数据具有很好的时空分布(观测数据更多，在近场和中远场都有一定分布，且观测持续达6年多之久)，而且通过断层面库仑应力的先验约束有效地区分黏弹性松弛和震后余滑两种机制对震后形变的贡献；考虑断层南北两侧介质流变性质的差异，更符合区域介质结构和性质的客观存在.因此，我们的研究结果应当是更精确和可靠的.

本研究最终得到柴达木盆地下地壳的长期黏滞系数为1.5×10²⁰Pa·s，巴颜喀拉-羌塘块体下地壳的长期黏滞系数为1.5×10¹⁹Pa·s.这一结果比下地壳流模型所主张的该地区下地壳黏滞系数(10¹⁷Pa·s量级，如Shen et al., 2001)高2~3个数量级，意味着黏性下地壳流在这一地区难以普遍存在，也不是这一地区地壳形变的主要机制.当然，由于本研究假定下地壳与上地幔黏滞系数相等，因而这一组黏滞系数的结果是下地壳和上地幔的平均，进一步的研究需要构建区分下地壳和上地幔的流变学模型.但由于观测数据在空间上的覆盖范围决定了其对下地壳黏滞系数比上地幔更为敏感，因此结果对于下地壳黏滞系数的约束更为可靠.

同时，昆仑山以南巴颜喀拉-羌塘块体的下地壳黏滞系数比昆仑山以北柴达木盆地的下地壳黏滞系数小大约1个量级，意味着巴颜喀拉-羌塘块体下地壳的流变性质更强.这与其增厚的地壳对应的高温高压导致的部分熔融相关，可以解释该地区下地壳剪切波速度和电阻率都低于柴达木盆地下地壳相关数值的现象(Unsworth et al., 2004; Sun et al., 2010)；也可以解释介质各向异性跨昆仑山的横向不均一性：在巴颜喀拉-羌塘块体下地壳S波分裂相当强烈，但在柴达木盆地下地壳就大为减弱(McNamara et al., 1994; Li et al., 2011)，这应当是由于巴颜喀拉-羌塘块体下地壳的相对高温导致在构造应力作用下晶体更倾向于定向排列的结果.

巴颜喀拉-羌塘块体和柴达木盆地的下地壳流变性质差异较大，因而其形变模式也有所差别.Ge等(2015)根据覆盖青藏高原的GPS观测计算得到青藏高原地表的应变率分布，该结果显示巴颜喀拉-羌塘块体的北西西-南东东向拉张率和北北东-南南西向缩短率相当，同时柴达木盆地的北北东-南南西向缩短率则显著大于北西西-南东东向拉张率.这一结果与沈正康等(2003)得到的结论一致.这可以解释为，由于相对于柴达木盆地，巴颜喀拉-羌塘块体下地壳介质流变性质更强，其形变更趋近于连续形变模式，在介质体积不变的情况下，主应变方向上的缩短率和其正交方向上的拉张率几乎相等；柴达木盆地介质强度高，其形变模式更接近块体运动模式，因而北西西-南东东向的拉张形变率主要集中在断层附近，在盆地内则迅速衰减.因而，由于介质的高强度，柴达木盆地的北西西-南东东向拉张率小于北北东-南南西向缩短率，区域内北北东向的构造挤压可能主要由跨活动构造(包括盆地边缘与内部断层)的全地壳缩短(图 1)而不是南东东向的伸展变形来吸收.

本研究通过采用GPS震后观测和库仑应力先验约束的方法模拟2001年可可西里地震震后形变，得到以下主要结论：

(1) 青藏高原北部下地壳、上地幔介质流变性质存在显著的横向不均一性.其中，东昆仑断层北侧柴达木盆地地区下地壳-上地幔短期和长期黏滞系数分别为5×10¹⁸ Pa·s和1.5×10²⁰ Pa·s，90%置信区间分别为4×10¹⁸~9×10¹⁸Pa·s和＞6×10¹⁹Pa·s；东昆仑断层南侧巴颜喀拉-羌塘块体下地壳-地幔短期和长期黏滞系数分别为1.5×10¹⁸ Pa·s和1.5×10¹⁹ Pa·s，90%置信区间分别为8×10¹⁷~2×10¹⁸Pa·s和1×10¹⁹~2×10¹⁹Pa·s.这一结果表明柴达木盆地下方介质的流变性质显著弱于巴颜喀拉-羌塘块体；并且，无论是柴达木盆地还是巴颜喀拉-羌塘块体，其下地壳的短期黏滞系数都大大小于长期黏滞系数，反映了由于震后黏弹性松弛使得应变率提高，介质的有效黏滞系数大大减小.

(2) 本研究得到的下地壳介质黏滞系数比下地壳流模型所主张的10¹⁷Pa·s高2~3个量级，意味着黏性下地壳流在巴颜喀拉-羌塘块体和柴达木盆地难以普遍存在，也不是该地区地壳形变的主要机制.相对于柴达木盆地，巴颜喀拉-羌塘块体下地壳介质流变性质更强，其形变更趋近于连续形变模式，区域内北北东向的构造挤压主要由南东东向的伸展变形来吸收.而柴达木盆地介质强度高，其形变模式更趋近块体运动模式，区域内北北东向的构造挤压可能主要由跨活动构造(包括盆地边缘与内部断层)的全地壳缩短来吸收.

感谢中国地震局第一监测中心、中国地震局第二监测中心和中国地震局武汉地震研究所的同仁为采集GPS数据所付出的艰辛和努力.