GPS约束下九寨沟地区断裂带现今运动速率的非连续接触模拟研究
孟庆筱 党学会
摘要：以中国地壳运动观测网络2009—2013年GPS观测数据为边界条件，使用非连续接触有限元技术构建九寨沟地区二维有限元模型，在不确定性分析的基础之上，计算区内主要断裂带现今运动速率。研究结果表明：在巴颜喀拉块体整体近于NE向的推挤过程中，九寨沟地区的塔藏断裂、虎牙断裂、树正断裂均呈现为较高的左旋走滑兼具挤压的现今运动特征；岷江断裂、龙日坝断裂和龙门山断裂则呈现为右旋走滑兼有挤压的运动特性。结合区域主应变率计算结果，发现九寨沟地区仍然具有较高的应变积累背景。树正断裂作为2017年8月8日九寨沟M7.0地震的发震断层，其现今左旋滑动速率为3.0 mm/a，与东昆仑断裂带玛沁—玛曲段附近的左旋走滑速率4.1mm/a基本匹配，说明该断裂可能是东昆仑断裂带东端分支断裂之一，而东昆仑断裂与虎牙断裂之间的历史地震空区可能已被九寨沟地震事件贯通。
关键词：九寨沟地震；GPS；有限元；接触单元；树正断裂
中图分类号：P315.725 文献标识码：A 文章编号：1000-0666（2018）03-0390-08
0 引言
2017年8月8日四川省阿坝州九寨沟县发生M7.0地震，震中位于（33.22°N、103.83°E），震源深度约20 km。该地震是继2008年汶川M8.0地震和2013年芦山M7.0地震之后发生在巴彦喀拉块体东部的又一强震事件。巴彦喀拉块体在印度—欧亚板块碰撞及高原侧向扩展作用下呈现为现今整体NE向的运动特征，其东部受到华南地块的阻挡，地壳呈近EW向挤压缩短变形，造成块体东缘密集分布的活动构造和频繁的强震活动（徐锡伟等，2017；易桂喜等，2017；屈勇，朱航，2017）。随着区域构造应力的长期积累、集中和加强，应变能最终可能在特定构造部位处以中强地震事件的形式突然释放出来。
中国地壳运动观测网络GPS观测结果可以提供高精度、大范围和准实时的地壳运动定量数据，利用复测获得的地壳运动速度场，通过研究断裂带运动学特征，可以为认识块体内部动力学驱动机制提供重要依据，对于未来区域地震危险性判断具有重要意义（陈长云等，2012）。因此，本文尝试以中国地壳运动观测网络2009—2013年GPS观测数据为基础数据，在忽略深度对断层活动差异性影响的基本假设下，使用非连续接触技术构建九寨沟地区二维有限元模型。在对断裂带摩擦系数等重要参数的不确定性分析基础之上，将有限元模拟结果作为GPS速度场强约束下复杂断裂带滑动速率的有效近似，分析九寨沟地区的强震危险性。
1 模型与方法
1.1 九寨沟地区二维非连续接触有限元模型
考虑到研究区内断裂构造十分发育，为了提高模型的可分析性，在保证模型可以客观反映基本问题的前提下对次级断裂带进行简化，纳入模型的断层共11条，如图1所示。为了清晰地展示九寨沟地区有限元模型，图2给出了格网空间步长为5 km的示意图，而实际后续研究中笔者使用更加精细的有限元模型，空间步长取0.5 km。在ANSYS有限元环境下，本文模型中构造块体和断裂过程带采用弹性连续本构关系，使用平面应变单元对研究区进行单元剖分。
为了获得九寨沟地区断裂带的现今形变特征，使用接触单元构建断裂带。在对非连续接触的摩擦方面进行控制时，为了处理断裂带两盘的接触和相对滑动问题，采用基于接触力的对称罚函数法，并将接触摩擦力考虑为库仑摩擦，使用ANSYS中Conta175和Targe169接触单元进行建模（安关峰等，2001）。
由于断裂带内介质常常发育有糜棱岩、断层泥等软弱介质，其物理力学性质与周围地质体具有较大差异，为了更加真实地反映断裂带附近破碎岩体较围岩强度更低的实际情况，本文通过介质物性参数弱化的方法在模型中构建了断裂过程带，如图3所示（Vermilye，Scholz，1998）。
1.2 介质物理力学参数
为了尽可能降低有限元模型的参數不确定性，笔者使用中国地震局工程力学研究所国家强震动台网中心提供的2017年8月8日九寨沟M7.0地震观测波形计算了研究区域内的P波与S波波速（Kenichi et al，2010），用于有限元弹性介质的物理力学参数设置。图4给出了P波与S波波速计算结果及其所用强震观测时程。
根据图4，将研究区域基岩块体的P波和S波波速取为6.1 km/s和3.5 km/s。参考徐晶等（2017）给出的巴彦喀拉块体东部分层介质模型参数，将介质密度取作2.8×103 kg/m3。根据Kenichi等（2010）对2000年日本鸟取MW6.6地震事件的研究结果，破裂过程带中存在的大量裂隙将会导致S波波速显著降低，可以将破裂过程带的S波波速设置为基岩块体区域的构建破裂过程带。因此，本文将九寨沟地区断裂带的破裂过程区的S波波速设置为基岩块体的0.5倍，即1.75 km/s。
断裂过程带宽度设置方面，Vermilye和Scholz（1998）将断裂带最大裂隙密度作为断裂过程带宽度的阈值标记，给出了破裂过程带宽度与断层长度之间的统计分布，如图5所示。据此，本文将九寨沟地区断层破裂过程带宽度设置为断裂带长度的1%。进一步参考雷显权等（2011）的研究，将断裂带接触单元法向接触刚度取为75 GPa/m，切向接触刚度取为3.5 GPa/m。根据刘晓红等（1986）使用双剪法摩擦实验得到的断层泥摩擦系数平均值，将区内各断裂的摩擦系数统一为0.7。
1.3 边界条件
本文采用中国地壳运动观测网络2009—2013年GPS观测数据，利用GAMIT/GLOBK和QOCA软件进行处理，得到了相应的速度场结果，其精度约为1 mm/a。为了尽可能剔除2008年汶川M8.0地震同震形变的影响，解算过程中通过插值获得各个GPS流动站处汶川M8.0地震的同震位移场，进而对同震位移大于3 mm/a的测站加以改正（国家重大科学工程“中国地壳运动观测网络”项目组，2008；江在森等，2009）。通过坐标旋转，得到区域无旋转坐标框架下GPS速度场（陈长云等，2012），如图6所示。将靠近模型矩形区域边界处的GPS测点的实际观测速度作为距离最近处节点的位移约束，通过双线性插值方法得到关键控制点之间其他节点的边界条件。
本文研究区域距离汶川地震震中较近，虽然在解算过程中对汶川地震同震形变进行了剔除，但是考虑到震后影响的复杂性，未能找到有效模型对汶川地震震后效应进行改正。因此本文所使用的GPS观测数据既包含了汶川地震后的变形，也包含了一部分汶川地震事件中尚未扣除干净的同震位移，这可能会导致最终断裂带滑动速率计算结果偏大（武艳强等，2012）。
2 计算结果与讨论
2.1 不确定性分析
为了分析参数不确定性对有限元模拟结果的影响，图7给出了研究区内21个GPS测站有限元模拟结果和GPS观测结果在EW分量和SN分量上的误差对比情况。由图可见，边界条件附近两者一致性较好，而区域中心位置附近两者的差异性则相对有所增大，但是模拟值总体上处于GPS观测结果解算误差范围之内，说明当前参数条件下的模拟结果具有较好的可靠度，可以用于进一步的分析。
本文尝试计算了断裂带接触单元摩擦系数取0.5，0.6，0.7和0.8这4种情况下的地壳运动速度场的模拟结果，并检测与实际观测结果的一致性。结果表明，摩擦系数取0.7时，模拟结果与观测结果的一致性最好，但是摩擦系数取值对计算结果的影响作用较小，其原因可能在于本文模拟过程未能进行断层闭锁到滑动的动态模拟。
2.2 区内断裂运动速度分析
本文利用九寨沟地区二维有限元模型中断裂带两盘接触单元的计算得到断裂现今运动特征与运动速率，主要包括平行于断层走向的走滑分量（右旋/左旋）Δvτ以及垂直于断层走向的压张分量（挤压/拉张）Δvn。
图8给出了九寨沟地区断裂带现今走滑速率及压张速率，其中线宽代表速率大小，断裂线条越宽则表示断裂此处运动速率越大。按照图8的现今运动性质计算结果，研究区内断裂带可以分为左旋走滑兼有挤压、左旋走滑兼有轻微拉张、右旋走滑兼有挤压3组。其中，东昆仑断裂、白龙江断裂、光盖山—迭山断裂、哈南—青山湾—稻畦子断裂、塔藏断裂、虎牙断裂、树正断裂为左旋走滑运动兼有挤压的特征，左旋走滑速率分别为4.1，0.8，1.5，2.5，4.0，2.5，3.0 mm/a，垂直于断裂走向的挤压速率分别为0.5，2.0，0.5，2.4，2.3，2.8，0.2 mm/a；雪山梁子断裂为左旋走滑兼有轻微拉张的现今运动特性，其左旋走滑速率为0.8 mm/a，拉张运动速率为0.2 mm/a；岷江断裂、龙日坝断裂和龙门山断裂则呈现为右旋走滑兼有挤压的运动特性，右旋走滑速率分别为2.1，4.8，1.5 mm/a，垂直于断裂走向的挤压速率分别为0.5，1.5，3.1 mm/a（表2）。
本文关于东昆仑断裂、白龙江断裂、光盖山—迭山断裂、哈南—青山湾—稻畦子断裂、塔藏断裂、虎牙断裂、龙日坝断裂、龙门山断裂无论是现今运动性质还是运动速率的分析，均与前人使用最邻近断层两侧的GPS测站位移速率计算方法以及断裂位错模型得到研究结果具有较好的一致性（程佳等，2009；陈长云等，2012；易桂喜等，2017），符合前人研究中巴彥喀拉块体向东方向运动过程中受到华南块体阻挡而强烈变形的相关认识。岷江断裂的运动特性存在争议，本文结果与杜方等（2009）和陈长云等（2012）根据GPS观测数据得到的岷江断裂右旋走滑速率分别为2.0，2.5 mm/a一致，而与赵小麟等（1994）利用地质、地貌资料得到的左旋走滑运动特征相左。考虑到地震地质资料反映的是断裂活动的长期运动特征，而GPS观测数据获得的是断裂带现今运动速率，该差异性说明岷江断裂除具有挤压特性以外，自形成以来可能经历了左旋与右旋走滑兼有的多期运动（陈长云等，2012；徐锡伟等，2017）。
2.3 树正断裂与九寨沟地震
树正断裂作为2017年8月8日九寨沟M7.0地震的发震断层，是位于塔藏断裂与岷江断裂之间的一条NW—SE向无名断裂，九寨沟地震之前该断裂无相关研究，震后被正式命名为树正断裂（易桂喜等，2017）。本文计算得到的该断裂左旋走滑速率为3.0 mm/a，同时伴有挤压作用。图9给出了九寨沟地区区域水平主应变场的有限元模拟结果。
由图9可知，2009—2013年九寨沟地区的水平主压应变率约为（2.0～9.0）×10-8/a，水平主张应变率约为（2.5～6.5）×10-8/a，该结果与陈长云等（2012）利用GPS观测获得的结果基本一致，说明巴彦喀拉块体整体处于近于NE向的推挤过程中，九寨沟地区仍然具有较高的应变积累背景。该区域构造背景下，2017年九寨沟M7.0地震的震源机制解为纯走滑型，且与余震带长轴方向一致的节面走向均值约为152.3°（易桂喜等，2017），与树正断裂较高的左旋走滑速率以及微弱的挤压速率特征一致。
徐锡伟等（2017）认为其成因机制主要在于树正断裂北段走向偏西而南段偏北，在巴彦喀拉块体整体近于向东推挤过程中，NNW向的北段左旋运动转化为近EW向的地壳缩短和岷山垂直隆升作用。同时考虑到塔藏断裂、虎牙断裂以及本文得到的树正断裂左旋滑动速率3.0 mm/a与东昆仑断裂带玛沁—玛曲段附近左旋走滑速率4.1 mm/a基本匹配，说明塔藏断裂、虎牙断裂以及树正断裂可能是东昆仑断裂带东端分支断层之一。与此同时，2017年九寨沟M7.0地震及其余震密集带恰好填补了虎牙断裂带西北段历史地震空区，故树正断裂可能是虎牙断裂的NW延伸的隐伏部分，而九寨沟地震的发生则可能贯通了虎牙断裂与东昆仑断裂（季灵运等，2017；徐锡伟等，2017），其区域地震危险性值得进一步关注。
3 结论
本文利用中国地壳运动观测网络2009—2013年GPS观测数据，使用二维非连续接触数值模拟方法计算了九寨沟地区断裂构造现今运动形变特征，得到以下主要结论：
（1）东昆仑断裂、白龙江断裂、光盖山—迭山断裂、哈南—青山湾—稻畦子断裂、塔藏断裂、虎牙断裂、树正断裂为左旋走滑兼有挤压的运动特征，左旋走滑速率分别为4.1，0.8，1.5，2.5，4.0，2.5，3.0 mm/a，垂直于断裂走向的挤压速率分别为0.5，2.0，0.5，2.4，2.3，2.8，0.2 mm/a；岷江断裂、龙日坝断裂和龙门山断裂则呈现为右旋走滑兼有挤压的运动特性，右旋走滑速率分别为2.1，4.8，1.5 mm/a，垂直于断裂走向的挤压速率分别为0.5，1.5，3.1 mm/a。结合区域主应变率计算结果，说明九寨沟地区仍然具有较高的应变积累背景。
（2）樹正断裂作为2017年8月8日九寨沟M7.0地震的发震断层，其现今左旋滑动速率约为3.0 mm/a，与东昆仑断裂带玛沁—玛曲段附近的左旋走滑速率4.1 mm/a基本匹配，说明该断裂可能是东昆仑断裂带东端分支断裂之一，可能是虎牙断裂北西段的隐伏部分，九寨沟地震事件可能已将东昆仑断裂与虎牙断裂贯通，其区域地震危险性值得进一步关注。
本文在撰写过程中得到中国地震局第一监测中心陈长云博士和占伟博士的帮助，在此向他们表示衷心感谢。
参考文献：
安关峰，高大钊，赵炯洋.2001.接触面弹粘塑性本构关系研究[J].土木工程学报，34（1）：88-91.
陈长云，任金卫，孟国杰，等.2012.巴颜喀拉块体北东缘主要断裂现今活动性分析[J].大地测量与地球动力学，32（3）：27-30.
陈顒，黄庭芳.2009.岩石物理学[M].合肥：中国科学技术大学出版社.
程佳，甘卫军，王泽河，等.2009.2001年昆仑山口西MS8.1地震前背景形变场的模拟研究[J].地震地质，31（1）：97-111.
杜方，闻学泽，张培震，等.2009.2008年汶川8.0级地震前横跨龙门山断裂带的震间形变[J].地球物理学报，52（11）：2729-2738.
国家重大科学工程“中国地壳运动观测网络”项目组.2008.GPS测定的2008年汶川MS8.0级地震的同震位移场[J].中国科学：地球科学，38（10）：1195-1206.
季灵运，刘传金，徐晶，等.2017.九寨沟MS7.0地震的Insar观测及发震构造分析[J].地球物理学报，60（10）：4069-4082.
江在森，方颖，武艳强，等.2009.汶川8.0级地震前区域地壳运动与变形动态过程[J].地球物理学报，52（2）：505-518.
雷显权，陈运平，赵炯洋.2011.天山现今地壳变形的非连续接触模型模拟[J].中南大学学报（自然科学版），42（3）：2754-2762.
刘晓红，方亚如，蔡戴恩，等.1986.海原断裂带断层泥摩擦特性的研究[J].西北地震学报，8（3）：50-54.
屈勇，朱航.2017.巴颜喀拉块体东—南边界强震序列库仑应力触发过程[J].地震研究，40（2）：216-225.
武艳强，江在森，杨国华，等.2012.南北地震带北段近期地壳变形特征研究[J].武汉大学学报（信息科学版），37（9）：1046-1048.
徐晶，邵志刚，刘静，等.2017.基于库仑应力变化分析巴颜喀拉地块东端的强震相互关系[J].地球物理学报，60（10）：4056-4068.
徐锡伟，陈桂华，王启欣，等.2017.九寨沟地震发震断层属性及青藏高原东南缘现今应变状态讨论[J].地球物理学报，60（10）：4018-4026.
易桂喜，龙锋，梁明剑，等.2017.2017年8月8日九寨沟M7.0地震及余震震源机制解与发震构造分析[J].地球物理学报，60（10）：4083-4097.
赵小麟，邓起东，陈社发，等.1994.岷山隆起的构造地貌学研究[J].地震地质，16（4）：429-439.
Kenichi F，Yuji K，Fusanori M.2010.Influence of fault process zone on ground shaking of inland earthquakes：Verification of Mj=7.3 Western Tottori Prefecture and Mj=7.0 West Off Fukuoka Prefecture earthquakes，southwest Japan[J].Engineering Geology，116（8）：157-165.
Vermilye J M，Scholz C H.1998.The Process Zone：a Microstructural View of Fault Growth[J].Journal of Geophysical Research，103（B6）：12223-12237.
Abstract The interpolated GPS observation data during the 2009-2013 period of the China crustal movement observation network was applied to gain the boundary conditions of 2-D element discontinuous model which was used to simulate the slip rates of the fault system in the Jiuzhaigou area and undertake uncertainty analysis over it.The study turned out that the NE direction movement of Bayan Har block had been hindered by the South China block.The eastern Kunlun fault，Tazang fault，Huya fault and Shuzheng fault showed a left-lateral strike slip and compressional characteristic，when the Minjiang fault，Longriba fault and Longmenshan fault showed a right-lateral strike slip and compressional characteristic.Combined with the calculation results of regional principal strain rate，it shows that the Jiuzhaigou area still has a high background of strain accumulation.As the seismogenic fault of the August 8，2017 Jiuzhaigou M7.0 earthquake，the left-lateral slip rates of Shuzheng fault is about 3 mm/year，which basically matched the slip rate of the eastern Kunlun fault.It is suggested that the fault may be one of the branch of the eastern Kunlun fault zone and the gap between the east Kunlun fault and the Huya fault may have been split by the Jiuzhaigou earthquake.
Keywords：Jiuzhaigou earthquake；GPS finite element；contact node；Shuzheng fault
摘要：以中国地壳运动观测网络2009—2013年GPS观测数据为边界条件，使用非连续接触有限元技术构建九寨沟地区二维有限元模型，在不确定性分析的基础之上，计算区内主要断裂带现今运动速率。研究结果表明：在巴颜喀拉块体整体近于NE向的推挤过程中，九寨沟地区的塔藏断裂、虎牙断裂、树正断裂均呈现为较高的左旋走滑兼具挤压的现今运动特征；岷江断裂、龙日坝断裂和龙门山断裂则呈现为右旋走滑兼有挤压的运动特性。结合区域主应变率计算结果，发现九寨沟地区仍然具有较高的应变积累背景。树正断裂作为2017年8月8日九寨沟M7.0地震的发震断层，其现今左旋滑动速率为3.0 mm/a，与东昆仑断裂带玛沁—玛曲段附近的左旋走滑速率4.1mm/a基本匹配，说明该断裂可能是东昆仑断裂带东端分支断裂之一，而东昆仑断裂与虎牙断裂之间的历史地震空区可能已被九寨沟地震事件贯通。
关键词：九寨沟地震；GPS；有限元；接触单元；树正断裂
中图分类号：P315.725 文献标识码：A 文章编号：1000-0666（2018）03-0390-08
0 引言
2017年8月8日四川省阿坝州九寨沟县发生M7.0地震，震中位于（33.22°N、103.83°E），震源深度约20 km。该地震是继2008年汶川M8.0地震和2013年芦山M7.0地震之后发生在巴彦喀拉块体东部的又一强震事件。巴彦喀拉块体在印度—欧亚板块碰撞及高原侧向扩展作用下呈现为现今整体NE向的运动特征，其东部受到华南地块的阻挡，地壳呈近EW向挤压缩短变形，造成块体东缘密集分布的活动构造和频繁的强震活动（徐锡伟等，2017；易桂喜等，2017；屈勇，朱航，2017）。随着区域构造应力的长期积累、集中和加强，应变能最终可能在特定构造部位处以中强地震事件的形式突然释放出来。
中国地壳运动观测网络GPS观测结果可以提供高精度、大范围和准实时的地壳运动定量数据，利用复测获得的地壳运动速度场，通过研究断裂带运动学特征，可以为认识块体内部动力学驱动机制提供重要依据，对于未来区域地震危险性判断具有重要意义（陈长云等，2012）。因此，本文尝试以中国地壳运动观测网络2009—2013年GPS观测数据为基础数据，在忽略深度对断层活动差异性影响的基本假设下，使用非连续接触技术构建九寨沟地区二维有限元模型。在对断裂带摩擦系数等重要参数的不确定性分析基础之上，将有限元模拟结果作为GPS速度场强约束下复杂断裂带滑动速率的有效近似，分析九寨沟地区的强震危险性。
1 模型与方法
1.1 九寨沟地区二维非连续接触有限元模型
考虑到研究区内断裂构造十分发育，为了提高模型的可分析性，在保证模型可以客观反映基本问题的前提下对次级断裂带进行简化，纳入模型的断层共11条，如图1所示。为了清晰地展示九寨沟地区有限元模型，图2给出了格网空间步长为5 km的示意图，而实际后续研究中笔者使用更加精细的有限元模型，空间步长取0.5 km。在ANSYS有限元环境下，本文模型中构造块体和断裂过程带采用弹性连续本构关系，使用平面应变单元对研究区进行单元剖分。
为了获得九寨沟地区断裂带的现今形变特征，使用接触单元构建断裂带。在对非连续接触的摩擦方面进行控制时，为了处理断裂带两盘的接触和相对滑动问题，采用基于接触力的对称罚函数法，并将接触摩擦力考虑为库仑摩擦，使用ANSYS中Conta175和Targe169接触单元进行建模（安关峰等，2001）。
由于断裂带内介质常常发育有糜棱岩、断层泥等软弱介质，其物理力学性质与周围地质体具有较大差异，为了更加真实地反映断裂带附近破碎岩体较围岩强度更低的实际情况，本文通过介质物性参数弱化的方法在模型中构建了断裂过程带，如图3所示（Vermilye，Scholz，1998）。
1.2 介质物理力学参数
为了尽可能降低有限元模型的参數不确定性，笔者使用中国地震局工程力学研究所国家强震动台网中心提供的2017年8月8日九寨沟M7.0地震观测波形计算了研究区域内的P波与S波波速（Kenichi et al，2010），用于有限元弹性介质的物理力学参数设置。图4给出了P波与S波波速计算结果及其所用强震观测时程。
根据图4，将研究区域基岩块体的P波和S波波速取为6.1 km/s和3.5 km/s。参考徐晶等（2017）给出的巴彦喀拉块体东部分层介质模型参数，将介质密度取作2.8×103 kg/m3。根据Kenichi等（2010）对2000年日本鸟取MW6.6地震事件的研究结果，破裂过程带中存在的大量裂隙将会导致S波波速显著降低，可以将破裂过程带的S波波速设置为基岩块体区域的构建破裂过程带。因此，本文将九寨沟地区断裂带的破裂过程区的S波波速设置为基岩块体的0.5倍，即1.75 km/s。
断裂过程带宽度设置方面，Vermilye和Scholz（1998）将断裂带最大裂隙密度作为断裂过程带宽度的阈值标记，给出了破裂过程带宽度与断层长度之间的统计分布，如图5所示。据此，本文将九寨沟地区断层破裂过程带宽度设置为断裂带长度的1%。进一步参考雷显权等（2011）的研究，将断裂带接触单元法向接触刚度取为75 GPa/m，切向接触刚度取为3.5 GPa/m。根据刘晓红等（1986）使用双剪法摩擦实验得到的断层泥摩擦系数平均值，将区内各断裂的摩擦系数统一为0.7。
1.3 边界条件
本文采用中国地壳运动观测网络2009—2013年GPS观测数据，利用GAMIT/GLOBK和QOCA软件进行处理，得到了相应的速度场结果，其精度约为1 mm/a。为了尽可能剔除2008年汶川M8.0地震同震形变的影响，解算过程中通过插值获得各个GPS流动站处汶川M8.0地震的同震位移场，进而对同震位移大于3 mm/a的测站加以改正（国家重大科学工程“中国地壳运动观测网络”项目组，2008；江在森等，2009）。通过坐标旋转，得到区域无旋转坐标框架下GPS速度场（陈长云等，2012），如图6所示。将靠近模型矩形区域边界处的GPS测点的实际观测速度作为距离最近处节点的位移约束，通过双线性插值方法得到关键控制点之间其他节点的边界条件。
本文研究区域距离汶川地震震中较近，虽然在解算过程中对汶川地震同震形变进行了剔除，但是考虑到震后影响的复杂性，未能找到有效模型对汶川地震震后效应进行改正。因此本文所使用的GPS观测数据既包含了汶川地震后的变形，也包含了一部分汶川地震事件中尚未扣除干净的同震位移，这可能会导致最终断裂带滑动速率计算结果偏大（武艳强等，2012）。
2 计算结果与讨论
2.1 不确定性分析
为了分析参数不确定性对有限元模拟结果的影响，图7给出了研究区内21个GPS测站有限元模拟结果和GPS观测结果在EW分量和SN分量上的误差对比情况。由图可见，边界条件附近两者一致性较好，而区域中心位置附近两者的差异性则相对有所增大，但是模拟值总体上处于GPS观测结果解算误差范围之内，说明当前参数条件下的模拟结果具有较好的可靠度，可以用于进一步的分析。
本文尝试计算了断裂带接触单元摩擦系数取0.5，0.6，0.7和0.8这4种情况下的地壳运动速度场的模拟结果，并检测与实际观测结果的一致性。结果表明，摩擦系数取0.7时，模拟结果与观测结果的一致性最好，但是摩擦系数取值对计算结果的影响作用较小，其原因可能在于本文模拟过程未能进行断层闭锁到滑动的动态模拟。
2.2 区内断裂运动速度分析
本文利用九寨沟地区二维有限元模型中断裂带两盘接触单元的计算得到断裂现今运动特征与运动速率，主要包括平行于断层走向的走滑分量（右旋/左旋）Δvτ以及垂直于断层走向的压张分量（挤压/拉张）Δvn。
图8给出了九寨沟地区断裂带现今走滑速率及压张速率，其中线宽代表速率大小，断裂线条越宽则表示断裂此处运动速率越大。按照图8的现今运动性质计算结果，研究区内断裂带可以分为左旋走滑兼有挤压、左旋走滑兼有轻微拉张、右旋走滑兼有挤压3组。其中，东昆仑断裂、白龙江断裂、光盖山—迭山断裂、哈南—青山湾—稻畦子断裂、塔藏断裂、虎牙断裂、树正断裂为左旋走滑运动兼有挤压的特征，左旋走滑速率分别为4.1，0.8，1.5，2.5，4.0，2.5，3.0 mm/a，垂直于断裂走向的挤压速率分别为0.5，2.0，0.5，2.4，2.3，2.8，0.2 mm/a；雪山梁子断裂为左旋走滑兼有轻微拉张的现今运动特性，其左旋走滑速率为0.8 mm/a，拉张运动速率为0.2 mm/a；岷江断裂、龙日坝断裂和龙门山断裂则呈现为右旋走滑兼有挤压的运动特性，右旋走滑速率分别为2.1，4.8，1.5 mm/a，垂直于断裂走向的挤压速率分别为0.5，1.5，3.1 mm/a（表2）。
本文关于东昆仑断裂、白龙江断裂、光盖山—迭山断裂、哈南—青山湾—稻畦子断裂、塔藏断裂、虎牙断裂、龙日坝断裂、龙门山断裂无论是现今运动性质还是运动速率的分析，均与前人使用最邻近断层两侧的GPS测站位移速率计算方法以及断裂位错模型得到研究结果具有较好的一致性（程佳等，2009；陈长云等，2012；易桂喜等，2017），符合前人研究中巴彥喀拉块体向东方向运动过程中受到华南块体阻挡而强烈变形的相关认识。岷江断裂的运动特性存在争议，本文结果与杜方等（2009）和陈长云等（2012）根据GPS观测数据得到的岷江断裂右旋走滑速率分别为2.0，2.5 mm/a一致，而与赵小麟等（1994）利用地质、地貌资料得到的左旋走滑运动特征相左。考虑到地震地质资料反映的是断裂活动的长期运动特征，而GPS观测数据获得的是断裂带现今运动速率，该差异性说明岷江断裂除具有挤压特性以外，自形成以来可能经历了左旋与右旋走滑兼有的多期运动（陈长云等，2012；徐锡伟等，2017）。
2.3 树正断裂与九寨沟地震
树正断裂作为2017年8月8日九寨沟M7.0地震的发震断层，是位于塔藏断裂与岷江断裂之间的一条NW—SE向无名断裂，九寨沟地震之前该断裂无相关研究，震后被正式命名为树正断裂（易桂喜等，2017）。本文计算得到的该断裂左旋走滑速率为3.0 mm/a，同时伴有挤压作用。图9给出了九寨沟地区区域水平主应变场的有限元模拟结果。
由图9可知，2009—2013年九寨沟地区的水平主压应变率约为（2.0～9.0）×10-8/a，水平主张应变率约为（2.5～6.5）×10-8/a，该结果与陈长云等（2012）利用GPS观测获得的结果基本一致，说明巴彦喀拉块体整体处于近于NE向的推挤过程中，九寨沟地区仍然具有较高的应变积累背景。该区域构造背景下，2017年九寨沟M7.0地震的震源机制解为纯走滑型，且与余震带长轴方向一致的节面走向均值约为152.3°（易桂喜等，2017），与树正断裂较高的左旋走滑速率以及微弱的挤压速率特征一致。
徐锡伟等（2017）认为其成因机制主要在于树正断裂北段走向偏西而南段偏北，在巴彦喀拉块体整体近于向东推挤过程中，NNW向的北段左旋运动转化为近EW向的地壳缩短和岷山垂直隆升作用。同时考虑到塔藏断裂、虎牙断裂以及本文得到的树正断裂左旋滑动速率3.0 mm/a与东昆仑断裂带玛沁—玛曲段附近左旋走滑速率4.1 mm/a基本匹配，说明塔藏断裂、虎牙断裂以及树正断裂可能是东昆仑断裂带东端分支断层之一。与此同时，2017年九寨沟M7.0地震及其余震密集带恰好填补了虎牙断裂带西北段历史地震空区，故树正断裂可能是虎牙断裂的NW延伸的隐伏部分，而九寨沟地震的发生则可能贯通了虎牙断裂与东昆仑断裂（季灵运等，2017；徐锡伟等，2017），其区域地震危险性值得进一步关注。
3 结论
本文利用中国地壳运动观测网络2009—2013年GPS观测数据，使用二维非连续接触数值模拟方法计算了九寨沟地区断裂构造现今运动形变特征，得到以下主要结论：
（1）东昆仑断裂、白龙江断裂、光盖山—迭山断裂、哈南—青山湾—稻畦子断裂、塔藏断裂、虎牙断裂、树正断裂为左旋走滑兼有挤压的运动特征，左旋走滑速率分别为4.1，0.8，1.5，2.5，4.0，2.5，3.0 mm/a，垂直于断裂走向的挤压速率分别为0.5，2.0，0.5，2.4，2.3，2.8，0.2 mm/a；岷江断裂、龙日坝断裂和龙门山断裂则呈现为右旋走滑兼有挤压的运动特性，右旋走滑速率分别为2.1，4.8，1.5 mm/a，垂直于断裂走向的挤压速率分别为0.5，1.5，3.1 mm/a。结合区域主应变率计算结果，说明九寨沟地区仍然具有较高的应变积累背景。
（2）樹正断裂作为2017年8月8日九寨沟M7.0地震的发震断层，其现今左旋滑动速率约为3.0 mm/a，与东昆仑断裂带玛沁—玛曲段附近的左旋走滑速率4.1 mm/a基本匹配，说明该断裂可能是东昆仑断裂带东端分支断裂之一，可能是虎牙断裂北西段的隐伏部分，九寨沟地震事件可能已将东昆仑断裂与虎牙断裂贯通，其区域地震危险性值得进一步关注。
本文在撰写过程中得到中国地震局第一监测中心陈长云博士和占伟博士的帮助，在此向他们表示衷心感谢。
参考文献：
安关峰，高大钊，赵炯洋.2001.接触面弹粘塑性本构关系研究[J].土木工程学报，34（1）：88-91.
陈长云，任金卫，孟国杰，等.2012.巴颜喀拉块体北东缘主要断裂现今活动性分析[J].大地测量与地球动力学，32（3）：27-30.
陈顒，黄庭芳.2009.岩石物理学[M].合肥：中国科学技术大学出版社.
程佳，甘卫军，王泽河，等.2009.2001年昆仑山口西MS8.1地震前背景形变场的模拟研究[J].地震地质，31（1）：97-111.
杜方，闻学泽，张培震，等.2009.2008年汶川8.0级地震前横跨龙门山断裂带的震间形变[J].地球物理学报，52（11）：2729-2738.
国家重大科学工程“中国地壳运动观测网络”项目组.2008.GPS测定的2008年汶川MS8.0级地震的同震位移场[J].中国科学：地球科学，38（10）：1195-1206.
季灵运，刘传金，徐晶，等.2017.九寨沟MS7.0地震的Insar观测及发震构造分析[J].地球物理学报，60（10）：4069-4082.
江在森，方颖，武艳强，等.2009.汶川8.0级地震前区域地壳运动与变形动态过程[J].地球物理学报，52（2）：505-518.
雷显权，陈运平，赵炯洋.2011.天山现今地壳变形的非连续接触模型模拟[J].中南大学学报（自然科学版），42（3）：2754-2762.
刘晓红，方亚如，蔡戴恩，等.1986.海原断裂带断层泥摩擦特性的研究[J].西北地震学报，8（3）：50-54.
屈勇，朱航.2017.巴颜喀拉块体东—南边界强震序列库仑应力触发过程[J].地震研究，40（2）：216-225.
武艳强，江在森，杨国华，等.2012.南北地震带北段近期地壳变形特征研究[J].武汉大学学报（信息科学版），37（9）：1046-1048.
徐晶，邵志刚，刘静，等.2017.基于库仑应力变化分析巴颜喀拉地块东端的强震相互关系[J].地球物理学报，60（10）：4056-4068.
徐锡伟，陈桂华，王启欣，等.2017.九寨沟地震发震断层属性及青藏高原东南缘现今应变状态讨论[J].地球物理学报，60（10）：4018-4026.
易桂喜，龙锋，梁明剑，等.2017.2017年8月8日九寨沟M7.0地震及余震震源机制解与发震构造分析[J].地球物理学报，60（10）：4083-4097.
赵小麟，邓起东，陈社发，等.1994.岷山隆起的构造地貌学研究[J].地震地质，16（4）：429-439.
Kenichi F，Yuji K，Fusanori M.2010.Influence of fault process zone on ground shaking of inland earthquakes：Verification of Mj=7.3 Western Tottori Prefecture and Mj=7.0 West Off Fukuoka Prefecture earthquakes，southwest Japan[J].Engineering Geology，116（8）：157-165.
Vermilye J M，Scholz C H.1998.The Process Zone：a Microstructural View of Fault Growth[J].Journal of Geophysical Research，103（B6）：12223-12237.
Abstract The interpolated GPS observation data during the 2009-2013 period of the China crustal movement observation network was applied to gain the boundary conditions of 2-D element discontinuous model which was used to simulate the slip rates of the fault system in the Jiuzhaigou area and undertake uncertainty analysis over it.The study turned out that the NE direction movement of Bayan Har block had been hindered by the South China block.The eastern Kunlun fault，Tazang fault，Huya fault and Shuzheng fault showed a left-lateral strike slip and compressional characteristic，when the Minjiang fault，Longriba fault and Longmenshan fault showed a right-lateral strike slip and compressional characteristic.Combined with the calculation results of regional principal strain rate，it shows that the Jiuzhaigou area still has a high background of strain accumulation.As the seismogenic fault of the August 8，2017 Jiuzhaigou M7.0 earthquake，the left-lateral slip rates of Shuzheng fault is about 3 mm/year，which basically matched the slip rate of the eastern Kunlun fault.It is suggested that the fault may be one of the branch of the eastern Kunlun fault zone and the gap between the east Kunlun fault and the Huya fault may have been split by the Jiuzhaigou earthquake.
Keywords：Jiuzhaigou earthquake；GPS finite element；contact node；Shuzheng fault
