印度板块以每年约5 cm的速度向北漂移，持续碰撞并挤压中国大陆所在的欧亚板块 ，致使青藏高原内部物质向东南方向侧流^[1-7]，侧流过程中遇到其东边稳固的四川盆地阻挡，从而在龙门山断裂带上积累了大量应力。一旦断层面累积应力超过断层可承受最大阈值，会触发地震事件，释放出积蓄的巨大能量，引发各种自然灾害。2008年5月12日，位于青藏高原东北缘的龙门山断裂带上，以汶川县映秀镇为震中发生了一次具有右旋走滑兼逆冲运动性质的Mw7.9级特大地震事件^[1]。

大震发生后数年至数十年甚至更长时间尺度内，在震区及周边会观测到随时间呈指数或对数衰减的地表位移^[8-12]，即震后位移。引起震后位移效应的机制主要有孔隙回弹、余滑和黏弹性松弛。已有文献[13—20]研究发现，由地震导致流体流动引发的震后孔隙回弹通常在含有大量地下水的区域对垂直向形变有较大影响；震后余滑通常只局限在震后几个月到1年左右的短时间尺度内，且主要作用于同震破裂下倾的近场区域；震后黏弹性松弛对中、远场区域有持续影响，其作用时间可以延续几年、几十年甚至更久。

震后形变因包含有地球介质的流变信息，蕴含着地壳运动和深部构造信息，因此是研究地球介质流变性质及岩石圈动力学关系的一种重要手段^[18-21]。以震后位移为地表约束，结合不同的震后形变机制，可以对断层的动力学性质、地壳深部的流变结构等地球动力学问题有进一步认知。由于地球内部的未知性，近年来，随着大地测量观测技术不断发展进步，有关科研人员开始利用GNSS观测资料研究地球内部介质的差异性。发生在青藏高原东边界龙门山断裂带的汶川Mw7.9级地震为研究青藏高原东缘的岩石圈流变结构及龙门山断裂带的摩擦性质提供了难得的机会。文献[16]研究汶川震后8个月47个测站震后GNSS形变序列。文献[14]计算汶川震后1年16个测站三维位移，并以此进行反演约束。文献[21]利用GPS和InSAR数据联合反演同震滑动分布和黏弹性松弛，基于模型反演给出汶川震后2年时间尺度的InSAR形变场分布。文献[20]采用震后1.5年的InSAR资料和震后一年少量GNSS数据对汶川地震震后形变的一级特征进行了分析，发现川西高原和四川盆地的黏滞系数存在显著差异，相差两个数量级。尽管汶川地震已经过去十年了，但现有成果主要集中在震后早期1~2年短时间尺度内，缺少对汶川震后形变长时间尺度持续跟踪分析，并且所使用的GNSS数据主要集中在近场，缺乏中远场观测资料，对黏滞系数的约束存在不确定性。本文根据2008年汶川地震震后龙门山断裂带上、下盘区域收集到的109个测站观测资料，通过高精度GNSS数据处理并对震前、震后时间序列进行分析，获取各测站震后位移。根据位移场分布，进一步研究龙门山断裂上、下盘青藏高原东部和四川盆地深部介质差异性，推算汶川地震在未来几十年内对周边区域地壳形变的影响。

1 GNSS数据处理与时间序列分析 1.1 站点选取与数据处理

龙门山断裂带是南北地震带的重要组成部分，位于巴颜喀拉块体和四川盆地交界处，有历史地震记载以来，地震多发生于断裂带中南段及其北边岷江断裂间^[1]。本文使用109个GNSS测站观测资料，时间跨度为1999—2015年，GNSS测站观测资料统计见表 1。龙门山断裂带区域构造背景与GNSS测站分布情况如图 1所示，其中，连续观测站点共计28个，流动观测站点共计81个。采用高精度数据处理软件GAMIT/GLOBK10.4对收集的109个GNSS站点进行高精密数据处理，获取单日松弛解，并引进IGS站作为公共测站联合解算，通过坐标转换将解算结果归算得到ITRF08参考框架下的站点坐标^[23]。

1.2 时间序列分析

GNSS坐标时间序列包含站点长期构造运动、周期性季节变化、已知和未知阶跃及震后形变等。长期观测的GNSS单站、单分量坐标时间序列可用式(1)表示

式中，x(t)为历元t时刻测站某一分量坐标位置；C₁为坐标初始位置；C₂为测站长期线性速率；C₃和θ为周年变化的振幅与相位；C₄和ϕ为半周年变化的振幅与相位；D_i是t_i时刻因地震、仪器变更等原因产生的阶跃；H为阶梯函数；f为震后形变。

震后GNSS观测到的坐标时间序列包含震后形变f，从GNSS坐标时间序列中获取震后形变，需扣除测站长期构造运动、测站受周边地震与仪器更换等产生的阶跃、周期性季节性变化等影响，其中季节性变化主要影响垂直方向，本文研究水平方向坐标时间序列，故季节性变化忽略不计。在汶川地震震后形变研究过程中，根据1997—2007年中国大陆连续跟踪站和区域站的平均运动速率作为本文地壳运动长期运动速率构造背景值^[23]，采用小区域均匀应变率方法对震后新建GNSS观测站点进行内插处理，得到全部测站在震前的长期构造运动速率值。

本文研究2010—2015年期间震后形变，该时段内于2013年4月20日在龙门山前缘构造带南段雅安地区发生Mw6.6级芦山地震。通过计算发现，本次所选取站点中有9个测站受芦山地震同震影响明显。图 2(a)列举出了其中2个连续站点SCXJ和SCMX扣除长期线性速率后的坐标时间序列，图 2(b)列出了其中10个区域站点扣除长期线性速率后的坐标时间序列，图 2中可以看出所示测站具有明显的震后形变趋势。

2 震后形变计算

描述震后形变过程的函数模型主要有3类，分别为幂函数、指数函数和对数函数，3种函数模型各有优势。分析余震活动时常用幂指数函数；在多数岩体力学试验中，通常用指数函数模型来描述黏弹性松弛特性；对于理想的黏弹性Maxwell体，可用对数函数模型来描述其震后松弛变化过程。假设川西地区地壳结为Maxwell体，本文震后形变衰减过程可采用如下对数函数描述^[20]

式中，t为震后某一观测时刻；t_EQ为地震发生时刻；D为测站震后形变幅度值；τ为震后松弛时间；f(t, τ)表示扣除长期线性速率、已知或未知阶跃等影响后反映震后形变的坐标时间序列。不同测站震后衰减呈现不一致性，本文综合考虑，选取时间序列质量比较好且数据覆盖时段相对较长的站点来获取可靠的震后松弛时间τ。通过计算发现，测站取各自最优的τ值，与取最佳的公共τ=34d时，计算的2010—2015年震后累计形变值间差异性很小，大部分在0.5 mm以内，最大差异不超过1 mm。根据拟合残差最小原则求解出的最优τ值为34天，部分测站时间序列拟合如图 2所示，计算得到2010—2015年间各测站累计震后水平位移，如图 3所示。

根据图 3所示，在2010—2015年期间震后形变量较大的站点有H045、J417、H050等，其水平方向累计形变达到5~7 cm，主要分布在龙门山断裂带中段的上盘近场区域。龙门山断裂上盘，鲜水河断裂以东及岷江断裂以西区域中、远场震后形变普遍较大，达到2~5 cm，并随距离的增加呈衰减趋势。龙门山断裂上盘、岷江断裂以东，除青川附近靠近地震破裂带的少数站点形变稍大，如H010、H032等，其余站点的震后形变都普遍较小，基本在2 cm范围内，且有沿断裂带向北方向呈逆时针旋转的运动趋势，与同震形变运动趋势比较一致。龙门山断裂带南段与鲜水河断裂、安宁河断裂交界的三岔口区域，大部分站点震后形变量不足1 cm，且没有固定的运动趋势，说明该区域内震后形变比较微弱。相对上盘，龙门山断裂下盘计算的测站震后形变则非常微小，除个别测站观测到了比较明显的震后形变，其余测站震后形变基本都在5 mm以内，这与下盘测站位于相对稳定的四川盆地上比较一致。对图 3所示结果进行分析发现，以龙门山断裂带为界，断裂带上盘震后形变普遍比下盘大，表明龙门山断裂带上、下盘具有很强的不对称性，反映青藏高原东部和四川盆地深部流变性质的显著差异性。

3 地壳结构反演

大地震发生后引发的震后黏弹性松弛位移，是由地震破裂产生的同震应力作用于下地壳/上地幔的软流物质上，造成深部物质的流变过程，从而引起弹性层持续、缓慢的地表位移^[16]。黏弹性是地球介质的一种重要性质，黏滞系数是描述介质黏弹性质的重要参数，介质的黏弹性质对地球动力学过程及其演化具有重要作用^[17]。黏弹性质是随时间的变化而不断变化的，某一时间内计算出的黏滞系数只能反映特定时间段内该区域的流变结构。

本文拟研究汶川地震震后2~7年，即2010—2015年间震后形变，认为影响震后形变的主要机制为震后黏弹性松弛。已有相关研究表明，汶川地震震后松弛主要发生在上盘地下30~60 km深处^[18]。根据全球地壳模型crust1.0，并参考前人研究成果^[14-23]，建立龙门山上盘地壳分层结构模型(表 2)，h和H分别为中上地壳、下地壳/上地幔的厚度。根据介质的差异性将弹性的中上地壳划分为多层，赋予不同的介质参数。弹性层厚度h的搜索范围为30~60 km，黏滞系数搜索范围为10¹⁷~10²¹ pa·s。采用Wang等的PSGRN/PSCMP程序^[24]模拟计算汶川震后，龙门山上盘区域因黏弹性应力松弛所引发的地表位移，通过对参数区间(h，η)进行搜索计算，并计算拟合误差RMS，拟合误差取最小值所对应的弹性层厚度和黏滞系数即为本文龙门山断裂上盘区域最优估值。

文中采用2010—2015年间累计震后形变，约束川西高原的岩石圈流变一级结构，认为该时段的震后形变主要受震后黏弹性松弛效应控制。此外，为了避免靠近断层的GPS测站受震后余滑影响，还剔除距同震破裂区域约100 km以内的近场测站，最终选取龙门山断裂上盘中、远场40个测站进行反演约束并估计最优弹性层厚度和黏滞系数。震后形变的大小和分布特征除了受岩石圈流变结构影响外，还与同震破裂模型息息相关。为此，笔者考查了3个不同同震破裂模型对模拟结果的影响。模型一为文献[7]公布的断层模型参数、模型二为文献[25]公布的破裂模型、模型三为USGS公布的单断层模型(http://earthquake.usgs.gov/eqcenter/eqintheneus/2008/us2008ryan/timite_fault.php, 2008.)。其中模型一、二是基于GNSS、InSAR等大地测量资料约束反演得到的断层滑动模型，模型三是基于地震波数据模拟的断层破裂及静态滑动分布。分别采用3种破裂模型进行模拟计算，得到3种不同同震模型对应的最优弹性层厚度和黏滞系数估值(见表 3)，以及震后形变观测值与模拟值对比情况，如图 4所示。

图 4显示3种不同破裂模型计算的震后形变拟合值和模拟值均具有较好的一致性，除局部个别站点外，整体差异性不大。从残差分布图可以看出，龙门山断裂上盘、鲜水河断裂以东及岷江断裂以西区域测站模拟结果一致性较好，其中模型三在该区域内的拟合一致性最佳。岷江断裂及其周边区域地质构造环境极其复杂，断裂以东区域测站除H024、H007等个别站点模拟值相对观测值偏小，其余测站均偏大，推测可能是岷江断裂东、西两侧地下介质流变性质的差异性导致，岷江以东区域黏滞系数比岷江以西可能要稍大；H007、H004等测站拟合结果较差，可能是受地壳分层结构的影响，也可能是受区域站观测精度影响。鲜水河断裂周边测站距离汶川震中较远，不同破裂模型引起的震后形变变化不大；测站震后位移模拟值普遍小于观测值，由于观测值已经扣除2013年芦山Mw6.6地震同震形变影响，因而可能是该区域深部流变性质差异，以及鲜水河断裂长期构造运动导致。不同模型远场震后形变基本无差异，模型的差异主要影响上盘近场及岷江断裂以东区域，模型一考虑滑脱层的存在，故J412、J413、H037等测站黏弹性模拟值被低估；另更加准确的滑动模型在岷江断裂附近模拟效果改善有限，除H034拟合较好外，其余站点残差依然较大。尽管本文反演前已剔除同震破裂区域100 km以内的测站，但通过计算发现，模型一、二中少数靠近破裂带的测站其黏弹性模拟值与观测值仍然不一致，模拟值被明显低估，尤其是J413、J412；距同震破裂区域较近测站的残差近似认为是由余滑引起，估算模型一中测站最大余滑量占比约30%，模型二最大余滑量占比约25%。综合图 4得出，本文所研究的汶川地震破裂区100 km以外区域在2010—2015年的震后地壳形变以黏弹性松弛效应为主，震后余滑只在同震破裂有限区域内对少数测站有较明显影响。

4 分析与讨论 4.1 流变结构

图 3揭示川西高原和四川盆地震后形变呈现强烈的非对称性，反映青藏高原东部和四川盆地深部流变性质的显著差异，震后形变主要集中在龙门山断裂上盘区域，下盘四川盆地几乎没有观测到较明显的震后形变。四川盆地内的GNSS测站震后形变远小于川西高原观测值，说明川西高原的黏弹性系数相对较高，这与文献[20]分析震后1年GPS观测值得到的认识是一致的。本文选取的3种同震破裂模型计算出的震后黏弹性松弛效应可以较好地解释2010—2015年期间龙门山断裂上盘中远场震后形变，但不能很好地解释该时间段内靠近岷江断裂附近的测站，如SCSP、H030等。岷江断裂附近测站震后形变模拟值分布与文献[20]的比较相似，测站中出现的观测值与模拟值方位角偏差可能是受青藏高原东部横向黏度不均匀影响。同时，龙门山断裂以北区域发育有NWW向的东昆仑断裂、近南北向的岷江断裂、NNW向虎牙断裂等，其中东昆仑断裂东端向东南方向呈马尾状散开，该区域构造比较复杂，3种不同模型计算的拟合值与模拟值对比趋势比较一致，可能与该区域自身长期构造运动相关。

文献[28]采用三维有限元方法模拟了汶川震后初期14天的黏弹性松弛效应，给出成都平原与川西高原地区中下地壳黏滞系数分别为9×10¹⁸ pa·s和4×10¹⁷ pa·s，得出龙门山断裂两侧黏滞系数差异至少有1个量级。文献[21]采用简单的二元分层、横向均匀的地球模型，通过InSAR数据约束，采用同震和震后形变联合模型反演，初步获得中下地壳黏滞系数下限为为2.0×10¹⁸ pa·s。该黏滞系数结果与文献[28]川西高原的结果比较一致，但与成都平原地区相差一个量级。文献[20]采用震后1年的GPS观测数据以及震后一年半的InSAR观测资料系统研究了汶川地震震后形变机理，通过三维有限元构建龙门山断裂两侧黏滞系数差异的岩石圈流变结构，大地测量资料约束的四川盆地的黏滞系数不低于10²⁰ pa·s，而川西高原下地壳瞬态和稳态黏滞系数分别为4.4×10¹⁷ pa·s和1.0×10¹⁸ pa·s。文献[14]采用两层的地壳结构分层模型，以震后1年16个测站的GNSS形变为约束，得到弹性层最佳黏滞系数为1.8×10¹⁹ pa·s。本文选取的3种不同破裂模型计算给出的龙门山断裂上盘川西高原的黏滞系数量级均为10¹⁹ pa·s，与文献[28]和文献[14]的结果比较吻合。

本文筛选上盘距汶川地震破裂区100 km以外的40个测站进行反演约束，为了验证测站选取对计算本文黏弹系数结果的影响，同时增加了同震破裂区域10个近场测站的反演结果。以模型二为例，当RMS=16.8 mm时，h=31 km，η=1.6×10¹⁹ pa·s，弹性层厚度和黏滞系数取最优估值。通过计算发现，增加100 km以内的近场测站一起反演计算得到的黏滞系数值和扣除近场测站单独反演获取的结果变化不大。通过对格网搜索值结果进行分析发现，黏滞系数比弹性层厚度对反演结果要更加敏感，本文研究给出川西高原地区下地壳/上地幔黏弹性层黏滞系数下限为10¹⁹ pa·s。根据增加10个近场测站反演约束计算的结果，对比图 4分析发现，龙门山上盘余滑区域主要分布在汶川地震破裂带100 km以内，且有向东北方向的余滑趋势，受影响的测站有H010、H033、J417、SCMX、H045、H050，主要分布在四川茂县、青川、映秀等区域。

4.2 震后黏弹性松弛效应预测

选取上述模型二反演获取的最优估值，基于震后黏弹性松弛模型，分别计算2015—2018年3年尺度、2018—2028年10年尺度、2028—2038年10年尺度、2038—2058年20年尺度水平向震后形变，如图 5所示。从图中可以看出中场震后形变影响比较明显，远场震后形变量级相对较小。不同时间尺度下测站累计最大/最小震后形变及年平均速率见表 4。结果显示，随时间推移，震后形变年平均速率呈逐年减小趋势，直至趋于稳定。尽管由震后黏弹性松弛效应引起的震后位移年变化量较小，但长时间尺度累计下的震后形变绝对量依然较大，不可忽视。通过计算震后2018—2058年40年时间尺度累计位移，发现J413、H046、H037等测站震后位移最大能达到19 cm，远场测站H027、H028、H041等累计震后位移有5~8 cm。

5 结论

本文通过对收集的1999—2015年间GNSS观测资料进行高精度数据处理和时间序列分析，并以2010—2015年震后形变为约束，基于震后黏弹性模型松弛进行反演，得出如下结论：

(1) 通过时间序列分析，获取最佳震后松弛时间为34 d，各测站具有明显的震后形变趋势(图 2)。给出研究区域内109个测站在震后2010—2015年间累计震后形变(图 3)，其中水平向震后形变较大测站位于断裂带中部的上盘近、中场区域，形变量最大达7 cm。

(2) 不同区域测站震后形变的差异反映地下深部流变性质的差异，分析发现龙门山断裂上、下盘震后形变不在同一量级(图 3)，具有明显的强不对称性，表明在下地壳/上地幔这一层，下盘黏滞系数要明显高于上盘。

(3) 基于震后黏弹性松弛模型，以龙门山断裂上盘40个测站2010—2015年累计GNSS震后形变为约束，根据3种不同破裂模型，分别反演上盘中上地壳弹性层厚度和下地壳/上地幔黏滞系数最优估值，并计算3种不同模型下2010—2015年震后形变模拟值和观测值对比情况及残差分布(图 4)。

(4) 3种不同模型计算的2010—2015年震后形变拟合值和模拟值均具有较好的一致性，对不同破裂模型计算结果进行比较分析发现，不同模型对近场震后形变有一定影响，对远场形变几乎无影响，除个别测站有明显变化外，整体差异性不大。龙门山断裂上盘岷江断裂以西、鲜水河断裂以东区域，测站拟合情况基本一致。岷江断裂以东区域模拟值整体有偏大趋势，尤其在岷江断裂附近，表明岷江断裂东西两侧地下流变性质可能有一定差异。模型一考虑滑脱层的存在，J412、J413、H037等测站黏弹性模拟值被低估，可能含有少量震后余滑效应，其中最大余滑量级约30%；另更加准确的滑动模型，在岷江断裂附近模拟效果改善有限，除H034拟合较好外，其余站点残差依然较大。

(5) 选取上述模型二反演的最优估值进行震后黏弹性松弛效应进行分析预测(图 5)，其中h=33 km，η=1.6×10¹⁹ pa·s。汶川震后形变随时间呈递减趋势，尽管年平均速率较小，但随时间推移累计震后形变量依然较大，不容忽视，仅2018—2058年40年尺度累计震后位移最大能达到19 cm。

龙门山断裂带覆盖地域较广，汶川地震破裂情况比较复杂，尽管本文使用的GNSS观测资料在数量和时间尺度上有限，且受地壳分层结构局限，但根据反演得到的中上地壳弹性层厚度和上地壳/下地幔黏弹性层黏滞系数与已有研究结果基本一致，表明本文计算、模拟具有一定可靠性。通过计算分析发现，汶川地震在未来几十年内主要对龙门山断裂带上盘区域地壳形变影响较大，震后余滑局限在同震破裂100 km范围内分布。对汶川地震震后形变的持续深入研究，可以为研究龙门山断裂带运动学和区域地球动力学过程及其演化提供重要依据。

谨以此文献给汶川地震10周年纪念！

致谢: 感谢中国大陆构造环境监测网络(CMONOC)提供的GNSS观测数据。