2008-05-12龙门山中央断裂带和山前断裂上发生汶川M_W7.9地震，研究汶川地震震后形变不仅对理解震后机制具有重要意义，而且对青藏高原东部地区的区域流变动力学研究具有重要价值^[1]。在现有针对汶川地震的研究中，部分研究为早期1~2 a短期变化，部分研究为长期持续跟踪分析^[1-3]，缺乏将短期数据与长期数据进行联合分析的研究。

本文首先采用位置分布均匀、广泛且时间范围为2010~2015年的GNSS形变观测数据，建立精细的三维粘弹性模型，基于有限元方法进行震后粘弹性松弛模拟，获得最佳的中下地壳粘滞系数，并研究单一粘弹性松弛机制对震后形变场的影响；然后采用2008~2009年的GNSS形变观测数据，利用具有上述最佳中下地壳粘滞系数的三维粘弹性模型计算震后1 a内的粘弹性松弛效应。将观测值减去粘弹性松弛模拟值的剩余形变量作为约束，建立粘弹性松弛和余滑的双机制联合模型，并与运动学反演获得的单一震后余滑模型进行比较，分析震后1 a内粘弹性松弛与震后余滑2种机制对汶川地震震后形变的相对贡献。

1 数据来源

本文收集的形变观测资料包含52个GNSS观测站数据，包括2008~2009年汶川地震震后形变监测网中21个观测站资料^[4]和2008~2015年陆态网中31个观测站资料^[2]。目前震后形变理论主要包含震后余滑和震后粘弹性松弛，震后余滑主要在震后1~2 a内存在，2 a后一般以震后粘弹性松弛为主。但汶川地震震后1 a内的形变影响因素既包含震后余滑也包含震后粘弹性松弛，为了能在震后1 a内的形变研究中将两者区分开，本文首先进行震后2~7 a的粘弹性松弛模拟研究，再进行震后1 a内的余滑研究。

2 震后粘弹性松弛

本文采用2010~2015年GNSS震后形变观测数据，考虑到该期间的震后形变主要受粘弹性松弛机制影响，因此可进行单一粘弹性松弛的有限元模拟计算，通过计算得到的震后形变场研究汶川地震震后2~7 a的形变特征。

2.1 模型构建方法

首先建立汶川地震的断层模型，可分为青川-南坝、南坝-北川、北川-映秀及彭灌区域4个部分，其中前3个部分区域均为铲状曲面，第4部分区域为光滑平面，且彭灌区域断层平面与北川-映秀区域铲状断层在深度约20 km处相交。断层模型的几何特征见表 1。

根据地质资料及断层破裂位置和几何特征建立长和宽均为3 000 km、深度为400 km的模型块体(图 1(a))，基于汶川断层模型对块体进行切割，在青藏高原东部区域设置弹性上地壳、粘弹性下地壳和上地幔，在四川盆地区域设置弹性上地壳和粘弹性上地幔。对网格进行划分，断层附近的网格比远场网格密集，既可满足精度需求又可提高计算效率，其中断层区域网格划分单元最小为3 km，边界网格单元最大为40 km，共包含334 651个单元和59 535个节点。保持模型上表面自由无应力，固定底部和四周的节点，模型距离断层边界位置足够远以减少固定条件引起的影响。

为了更准确地模拟震后形变场，本文根据文献[5]中的同震滑动分布模型，采用同震形变场数据，利用Okada弹性半空间位错模型计算格林函数G，然后将反演的同震滑动分布模型作为计算粘弹性松弛的源断层模型(图 1(b))。

2.2 弹性层厚度选取

在确定地震的破裂源后，地球表面的震后形变将取决于深度的流变结构。分析表明，四川盆地的岩石圈结构和力学性质与青藏高原东部存在很大差异。参考文献[6]中的建模分层方法，本文共建立4个三维有限元模型，设置青藏高原东部块体的弹性层分别为25 km、30 km、33 km、35 km，粘性块体介质主要采用麦克斯韦体，粘滞系数搜索范围为10¹⁷~10¹⁹ Pa ·s。根据介质的差异性将上地壳、中下地壳和上地幔赋予不同的介质参数(表 2)，其中各层参数主要参考文献[2]。利用有限元软件Pylith^[7]模拟计算因粘弹性松弛引发的汶川地震震后地表位移，通过网格搜索方法分别对4个模型的最佳粘滞系数进行搜索，并计算拟合误差RMS，拟合误差为最小值时对应的弹性层厚度即为本文青藏高原块体弹性层的最优估值(图 2)。

本文利用2010~2015年的震后形变数据约束青藏高原东部的岩石圈流变结构，认为该时段的震后形变主要受粘弹性松弛效应控制。计算最优估值时所用的GNSS站点中有2个测站位于断层附近，其余29个测站均距离断层区域较远。研究表明^[2]，增加100 km以内的近场测站进行反演得到的粘滞系数值与去除近场测站单独反演获取的结果差异较小，因此本文的测站分布不会影响反演结果。比较4个弹性层模型的拟合误差可知(图 2)，最优弹性层厚度为25 km，对应的最优粘滞系数为4×10¹⁸ Pa ·s，RMS值为6.01 mm。

根据采用有限元方法建立的弹性层厚度为35 km的三维模型^[1]计算单一粘弹性松弛，结果表明，青藏高原东部中下地壳最优粘滞系数为1×10¹⁸ Pa ·s，然后建立联合模型计算得到青藏高原东部中下地壳对应的最优粘滞系数为2×10¹⁸ Pa ·s。温扬茂等^[8]采用将同震和震后形变进行联合的方法，获得中下地壳粘滞系数下限为2.0×10¹⁸ Pa ·s；Huang等^[6]采用震后1 a的GNSS观测数据及震后1.5 a的InSAR观测资料，基于有限元方法构建龙门山断裂带两侧粘滞系数差异的岩石圈流变结构三维模型，得到四川盆地粘滞系数不低于1.0×10²⁰ Pa ·s，而川西高原下地壳瞬态和稳态粘滞系数分别为4.4×10¹⁷ Pa ·s和1.0×10¹⁸ Pa ·s。本文建立的4个三维有限元模型计算粘弹性松弛对应的最佳粘滞系数为2×10¹⁸~4×10¹⁸ Pa ·s，与前人的研究结果基本一致。

2.3 震后粘弹性松弛位移场

基于青藏高原东部区域最优弹性层厚度和最优中下地壳粘滞系数的比较，本文选择弹性层厚度为25 km的三维粘弹性模型，青藏高原东部区域中下地壳、上地幔和四川盆地粘弹介质均采用麦克斯韦体，青藏高原东部中下地壳和上地幔粘滞系数分别设置为4×10¹⁸ Pa ·s和1×10²⁰ Pa ·s，四川盆地上地幔粘滞系数设置为1×10²⁰ Pa ·s。将2010~2015年的GNSS震后形变观测值作为模拟值，计算震后2~7 a粘弹性松弛影响的震后形变位移(图 3)。

从图 3(a)可以看出，多数测站分布于龙门山断裂带上盘中远场、鲜水河断裂带附近及岷江断裂以东区域，少数测站分布于四川盆地区域，龙门山断裂带上盘区域站点的震后形变拟合值最大约为3.2 cm，最小约为2 mm，与GNSS站点震后形变观测值最大约为3.9 cm、最小约为2 mm一致。结合图 3(b)可知，龙门山断裂带远场区域的拟合效果比中场区域好，且距离断裂区域越近，拟合效果越差，考虑到距离断裂区域越近其站点位移量越大，从而导致拟合误差也偏大，也可能由该区域深部结构的流变性质分布不同所致。鲜水河断裂带附近站点的拟合剩余形变场偏大，可能与鲜水河断裂带附近复杂的构造环境有关，也可能受该区域其他地震及震后影响。站点H043和H044位于龙门山断裂带下盘区域，考虑到站点位置处于相对稳固的四川盆地块体，站点震后形变量约为4 mm，形变量偏小且拟合效果不明显，表明龙门山断裂带上下盘存在明显的不对称性，反映青藏高原东部区域和四川盆地区域地下深部的流变结构性质存在很大差异。在岷江断裂附近及其以东区域，站点形变拟合值比观测值小，可能受该区域介质粘滞系数偏大的影响；从站点形变拟合值的方向来看，存在沿断裂带向北逆时针旋转的运动趋势，符合汶川地震右旋走滑的特征。

3 震后余滑

研究表明^{[3, 6, 9]}，余滑模型能够很好地拟合近场观测数据，但在远场区域拟合效果不佳；粘弹性松弛模型则相反，能够较好地拟合中远场观测数据，却无法拟合近场数据。因此对于汶川地震，单一机制无法很好地解释观测到的震后形变场，本文将联合震后余滑与粘弹性松弛，建立合理的双机制联合模型进行反演。

3.1 单一余滑机制

震后余滑是地震后同震应力变化引起无震滑移的过程，可通过对GNSS测量的累积震后位移进行反演来分析破裂断层表面的余滑分布。假设可通过弹性半空间位错模型来描述余滑，反演的目标函数可表示为：

$ \boldsymbol{F}(\boldsymbol{s}, \beta)=\|\boldsymbol{W}(\boldsymbol{G} \boldsymbol{s}-\boldsymbol{d})\|+\beta\left\|\nabla^{2} \boldsymbol{s}\right\| $

(1)

式中，s为估计的滑动分量；W为根据观测不确定性构造的权重矩阵；d为大地测量数据矢量；G为Green函数矩阵；▽²为用于平滑滑动模型的Laplacian算子的有限差分，且β为调整平滑度的因子。通过使用有界可变最小二乘(BVLS)算法获得最佳模型解，β最佳值应在模型粗糙度和数据失配之间进行权衡^[10]。

本文单一余滑模型采用运动学反演方法，认为汶川地震震后1 a内的形变均由余滑作用引起，故将汶川地震震后1 a的GNSS形变观测数据作为输入值。利用反演函数计算震后1 a内破裂断层面的余滑分布(图 4(d))，得到地面站点1 a内的震后形变值，进而分析汶川地震在单一余滑机制下震后1 a内的形变特征(图 4(a))。

3.2 粘弹性松弛与余滑联合机制

在研究粘弹性松弛时，利用汶川地震震后2~7 a的形变数据测得青藏高原东部弹性层最优厚度为25 km，龙门山断裂带上盘中下地壳最佳粘滞系数为4.0×10¹⁸ Pa ·s，此时拟合误差最小，RMS值为6.01 mm。这些结果与前人的研究结论相近^{[1, 11]}，故可认为得到的中下地壳粘滞系数具有可靠性。假设汶川地震震后7 a内龙门山断裂带中下地壳的粘滞系数不变且为4.0×10¹⁸ Pa ·s，选择弹性层厚度为25 km的三维有限元模型对震后1 a内的粘弹性松弛进行正演计算，再将震后1 a的观测数据减去粘弹性松弛的正演结果，得到剩余形变场。在只考虑粘弹性松弛与余滑2种机制影响的条件下，去除粘弹性松弛形变后的剩余形变场可认为均由余滑作用引起，将剩余形变场作为输入值，采用上述反演程序得到破裂断层面的余滑分布(图 4(c))及由余滑影响的地面震后形变结果。将正演粘弹性松弛的地表形变结果与剩余形变场反演得到的余滑影响下的地表形变结果相加，得到2种机制的联合结果，然后与震后1 a内GNSS观测结果进行对比，研究双机制下汶川地震震后1 a内的形变场特征(图 4(a))。

根据图 4(a)所示，将所用的52个地表站点划分为3个区域，距离断层破裂区域最远的站点拟合效果最差，距离较近的站点拟合效果较好，位于汶川地震破裂面附近的站点拟合结果与震后形变观测值基本一致，2种拟合方法均表现出由近到远拟合误差越来越大的特点。H027、H028、H029等远场站点的拟合效果较差，这可能是因为断层面滑脱层的水平宽度与延伸深度较小，使龙门山断裂带上盘区域远场位移受到较大影响^[1]；而在鲜水河断裂带区域的站点，如H053、H054、H056、H062等站点的拟合值普遍偏小，这可能是鲜水河断裂带区域地壳结构复杂、介质属性与附近区域存在较大差异所致。从单一余滑机制来看，余滑存在于断层面上，其影响范围一般覆盖近场及中场区域，对于远场站点的形变影响较小，因此出现蓝色拟合值由近到远逐渐减小的情况。从粘弹性松弛与余滑联合机制来看，虽然也存在拟合值渐变的情况，但从中远场拟合效果及图 4(b)可以看出，联合机制的结果比单一余滑机制下的拟合值偏大，且更接近震后形变观测值，表明多机制联合比单一余滑机制能更好地解释汶川地震震后1 a内的形变特征。

由图 4(c)和4(d)可以看出，2种拟合方法的余滑结果基本一致，且主要分布在同震过程中未滑动或滑动量较小的区域，包括青川东北部断坡区域、南坝-北川段断坡区域、虹口-映秀深部滑脱层区域及映秀西南部断坡区域等，其中单一余滑机制反演的余滑模型最大滑动量约为0.62 m，双机制联合反演的余滑模型最大滑动量约为0.54 m。从图 4还可以看出，断层浅部和深部均存在余滑分布，说明断层的震后余滑不仅发生在断层浅部，同时也可能发生在深部。青川东北部断坡区域可能是由断层破裂后东北部区域向北走滑形成；北川断坡区域可能是由同震破裂后存在的余震作用引起；虹口-映秀深部滑脱层区域可能是由于地震发生后其逆冲趋势依然存在，导致深部继续出现余滑运动。上述分析表明，余滑形变继承了同震形变的特征，从西南部逆冲过渡到中部逆冲兼走滑，东北部为纯走滑运动。

4 讨论 4.1 震后粘弹性松弛与震后余滑时间划分分析

根据文献[2]计算震后粘弹性松弛效应的方法，本文采用震后2~7 a的形变数据研究该时间段的震后粘弹性松弛效应，同时采用双机制联合方法对震后2 a内的粘弹性松弛和余滑效应的时间序列进行模拟分析。由图 5可知，H035站点距离破裂面较近，震后粘弹性松弛形变极小，大部分为震后余滑成分；SEEG站点距离破裂面较远，震后粘弹性松弛所占分量较大；WARI站点距离破裂面最远，震后粘弹性松弛成分甚至比震后余滑多。从增长趋势可以看出，震后1 a内余滑比重较大且增长明显，震后1~ a内余滑增长量逐渐减少并趋于平缓；而震后粘弹性松弛一直呈递增趋势，尤其在中远场区域增长趋势明显。各站点粘弹性松弛残差可近似认为由余滑引起，根据震后形变时间序列估计，在震后2~7 a内中场区域(以SEEG站点为例)余滑效应占比约36%，远场区域(以TAGO站点为例)余滑效应占比约25%。综合图 3和5可知，2010~2015年中远场区域震后形变以粘弹性松弛效应为主，因此利用震后2~7 a的形变数据研究震后粘弹性松弛效应具有合理性。由于在研究震后粘弹性松弛效应时采用的测站均分布于中远场区域，故对近场区域不作考虑。

因此，本文采用2010~2015年中远场区域站点数据研究震后粘弹性松弛效应，采用2008~2009年数据研究震后余滑和震后粘弹性松弛的短期作用特征。

4.2 粘弹性松弛结果分析

由图 2可知，弹性层厚度为25 km时拟合误差最小，为6.01 mm。根据crust2.0地壳模型可知，北川-青川地区上地壳弹性层厚度为17 km，映秀地区为22 km。夏婷婷^[12]通过建模得出龙门山断裂带上地壳厚度为22 km，且自北向南厚度逐渐变薄，说明龙门山断裂带上盘弹性层厚度并非均匀分布，北川-青川地区较浅，映秀地区较深。考虑到改变龙门山断裂带上盘弹性层厚度会影响中下地壳最优粘滞系数估计和震后粘弹性松弛效应，从而影响震后形变特征，建立非均匀厚度的弹性层具有重要意义。

利用有限元方法对单一粘弹性松弛进行震后2~7 a的正演模拟，结果表明，最大形变量为32 mm，最小为2 mm，与震后形变观测值较为接近。拟合结果与观测值存在误差，可能是受计算误差的影响，也可能是由三维建模中横向、纵向的分层粗糙导致，还可能与站点位置区域的地壳构造影响有关。下一步可根据汶川地震的层析成像结果建立更精细的三维模型，在断层两侧施加不均匀分布的介质，进而完善模拟效果。

4.3 单一余滑与多机制联合结果分析

以往研究中，直接由地表震后位移通过线性反演得到震后余滑是解释观测结果最方便、最有效的方法，但该方法缺乏物理基础，需要设置大量的自由参数。此外，由于忽略了粘弹性松弛的贡献，获得的余滑在大小和深度上往往偏大^[1]。通过对比单一余滑机制与双机制联合的模拟结果及剩余形变量(图 4(a)、4(b))可知，断层面附近测站的模拟结果均与观测值高度拟合，说明在断层面附近的震后形变粘弹性松弛作用较小，余滑起主导作用；对于中远场区域，双机制联合的模拟结果均比单一余滑机制的模拟结果更接近震后形变观测值，说明在中远场区域粘弹性松弛起到一定作用。从整体上看，多机制联合的拟合效果比单一余滑机制的拟合效果好，单一的形变机制不能很好地解释汶川地震震后形变场，多机制联合的方法能更好地说明汶川地震震后的形变特征。

5 结语

本文采用位置分布均匀、广泛且时间范围分别为2010~2015年和2008~2009年的31个GNSS形变观测资料及21个GNSS形变观测资料进行震后形变分析研究，得到以下结论：

1) 龙门山断裂带两侧的站点震后形变具有明显的不对称性，下盘形变值小于上盘，表明断层两侧上下盘的深部流变性质存在巨大差异。

2) 震后2~7 a粘弹性松弛模拟结果表明，龙门山断裂带上盘最佳弹性层厚度为25 km，中下地壳最佳粘滞系数为4.0×10¹⁸ Pa ·s，此时拟合误差RMS为6.01 mm。

3) 震后2~7 a粘弹性松弛模型在中近场区域拟合误差偏大，在远场区域拟合效果较好。

4) 震后1 a内的形变影响中余滑占主导作用；双机制联合方法得到的震后形变结果比单一余滑机制更接近实际观测值，粘弹性松弛和余滑的联合机制能更有效地解释震后形变特征。