引言

青藏高原是我国最主要的地震活动区，是我国板内6个Ⅰ级活动块体之一，也是我国地震活动最为强烈的块体． 该块体内自南向北可分为拉萨、 羌塘、 川滇、 巴颜喀拉、 东昆仑—柴达木和祁连山等6个Ⅱ级块体． 青藏块体地质构造变形复杂，断裂活动强烈，控制着一系列历史强震的发生. 1900年以来在青藏高原中部的巴颜喀拉块体内部及周缘发生了14次7级以上地震，尤其是近十几年来，接连发生了1997年11月8日西藏玛尼M_S7.5地震，2001年11月14日青海昆仑山口西M_S8.1地震，2008年5月12日四川汶川M_S8.0地震，以及2010年4月14日青海玉树M_S7.1地震，这一系列7级以上强震全部发生在巴颜喀拉块体边界断裂带上． 通过计算断层活动所引起的应力变化，研究人员发现，大地震之间是互相影响的(Stein，2003)． 而围绕巴颜喀拉块体边界发生的多次强震之间的相互影响和关联则是一个值得深入研究的科学问题，受到许多学者的关注． (邓起东等2010)指出，这一期地震活动与青藏高原中部巴颜喀拉块体的活动密切相关，青藏高原块体活动和地震活动与澳大利亚—印度板块对亚洲大陆的推挤作用相关，青藏高原和巴颜喀拉块体的强震活动与澳大利亚—印度板块边界苏门答腊强震活动相呼应，是未来一定时期内值得关注的主要地震活动区．

一些学者(王连捷等，2009，2010； 闻学泽等，2011； 戴黎明等，2011； 张晓亮等，2007)对近年来巴颜喀拉块体周缘发生的7级以上强震的发震机理、 发震过程及同震形变场进行了数值模拟分析研究． (郑勇等2006)利用有限元方法探讨了欧亚大陆的碰撞对中国大陆岩石圈形变和应力场的影响，并揭示了它们与强震活动性的关系，认为印度板块的碰撞对中国大陆的强震活动性有重要影响．

本文围绕巴颜喀拉块体，建立青藏高原目标区壳体有限元模型． 从区域动力学角度着手，使用有限元方法模拟澳大利亚—印度板块的推挤作用对巴颜喀拉块体内及边缘强震形势的影响，并对未来的地震危险区进行初步的探索．

1.   青藏高原线弹性壳体有限元模型

1.1   活动断裂

以青藏高原为目标建立了壳单元有限元模型，采用降刚法模拟强震的发生. 模型中的断裂主要是活动块体之间的活动断裂，包含有喜马拉雅俯冲带、 嘉黎断裂带、 可可西里断裂带、 甘孜玉树—鲜水河—安宁河—则木河—小江断裂带、 金沙江断裂带、 红河断裂带、 龙门山断裂带、 阿尔金断裂带、 东昆仑断裂带、 柴达木盆地北缘断裂带及祁连北缘断裂带等． 在模拟计算中，把活动断裂处理为一软弱带，通过物性参数的选取，降低其刚度，使其更加容易变形． 在二维有限元模型中为了更好地模拟活动断裂，以青藏高原块体分界线为目标，块体内的小断裂忽略处理．

1.2   本构关系

使用壳单元有限元模型，对复杂的地壳结构进行了简化处理，采用线弹性本构关系，应力-应变关系符合胡克定律，即

式中，σ_ij 为应力分量， D_ijkl 是应力应变矩阵，ε_kl 为应变分量． 对于各向同性材料的广义胡克定律可写为

This page contains the following errors:

Below is a rendering of the page up to the first error.

弹性模量、 泊松比等力学参数的关系为

式中，E为杨氏模量，ν为泊松比．

1.3   几何分区及其物性参数

根据张培震等(2003a)、 (邓起东等2010)对活动块体分区的研究结果，将青藏高原及邻区分为拉萨、 羌塘、 巴颜喀拉、 柴达木、 祁连、 塔里木、 川滇、 滇西、 华南和鄂尔多斯等构造区域，在本文中将所有大型断裂均处理为弱化带单元． 关于青藏高原地壳介质参数，本文主要依据前人的研究成果(李永华等，2006，2009； 刘宝峰等，2003； 苏伟等，2002； 王椿镛等，2008； 吴庆举，曾融生，1998； 胥颐等，2009； 张中杰等，2002)，文中所选的地震震源深度大多数在15 km左右，因此取15 km处的波速推算出杨氏模量及泊松比，按图1所示进行分区，弱化带杨氏模量为两侧分区杨氏模量平均值的1/3，弱化带内泊松比比周围区域稍高，介质物性参数见表1(弱化带介质参数略)．

图  1  模型分区示意图

Figure  1.  Sketch map of the area partition in the model

下载: 全尺寸图片 幻灯片

表  1  青藏高原及邻区地壳介质物性参数

Table  1.  Crustal physical properties in Qinghai Xizang plateau and surroundings

下载: 导出CSV 

| 显示表格

1.4   网格划分及边界条件

本文用ANSYS12.1软件进行数值模拟计算，整个有限元模型网格共划分为42 531个单元，21 544个节点，单元平均长度约20 km．

依据板块碰撞及活动块体理论，模拟印度板块推挤青藏高原，经过不同的尝试，在模型中施加了如下边界条件组合: 考虑到印度板块以N20°E的运动方向推挤中国大陆(张培震等，2002)，以及大西洋中脊向东的作用，对西边界进行法向约束，切向自由；北边界GPS位移较小，采用弹簧单元法向约束，但西段塔里木盆地较为稳定，西段边界施加的约束强于东段； 同样东边界四川盆地以北受鄂尔多斯的阻挡GPS位移很小，采用弹簧单元法向约束；四川盆地较为稳定，边界采用切向与法向全约束； 滇西块体西边界存在sagaing大型右旋走滑断层，采用法向约束，切向自由;作为印度板块推挤青藏高原绕喜马拉雅构造结物质东流的通道，东边界四川盆地以南和川滇块体南边界不施加特定约束;模型在垂直方向进行约束.经过上述组合约束后，模型的主体区域应力场和位移场与实际基本一致.

前人的研究结果认为印度板块以50mm/a的速度向北东方向推挤(Paul et al，2001;Banerjee，Bürgmann，2002)，本文中给定一个约5万年相对时间尺度，施加到实际模型中换算成力相当于在印度板块俯冲边界段施加100MPa的应力边界条件.以前的相关研究中，有学者也曾选取印度板块的挤压作用强度为80MPa的应力边界条件(张东宁，许忠淮，1999)和大于100MPa的应力边界条件(张东宁，高龙生，1989);Lu等(2011)与杨树新等(2012)也在印度边界施加了100MPa的应力.因此，100MPa只是一个相对的应力边界条件，在相对应力条件下模拟出的地震之间影响的强度也是相对的.

2.   数值模拟结果及分析

2.1   区域水平构造背景应力场模拟结果

在印度板块的推挤作用下，经过多年的变形过程，形成了以羌塘块体东部（玛尼—玉树—鲜水河断裂南侧)为中心地壳物质向东的流动带，反映了青藏高原的物质东移(张培震等，2003b； 陈连旺等，2010)． 本文在青藏高原巴颜喀拉块体及邻区即高原中部，从西向东横跨大型活动断裂带设计了5条路径线(图2)，分别描述水平构造背景场等效应力沿路径方向的变化特征．

图  2  路径示意图

Figure  2.  Sketch map of the path

下载: 全尺寸图片 幻灯片

等效应力可定义为

式中，σ₁，σ₂分别为最大及最小主应力．

等效应力的大小能够反映介质趋于屈服程度的强弱． 等效应力增大，有利于介质失稳而加速孕震进程； 反之，则等效应力减小，不利于介质失稳而减缓孕震进程．

1号路径位于巴颜喀拉块体北部，2号路径位于巴颜喀拉块体内部，3，4和5号路径纵跨巴颜喀拉块体，具体结果见图3a-e． 图3中分别给出了从西向东、 从南向北水平构造背景场等效应力沿路径节点的变化过程．

图  3  青藏高原中部等效应力变化曲线

(a) 1号路径 (b) 2号路径 (c) 3号路径 (d) 4号路径 (e) 5号路径

Figure  3.  Equivalent stress variation in the middle of Qinghai-Xizang plateau

(a) Path 1; (b) path 2; (c) path 3; (d) path 4; (e) path 5

下载: 全尺寸图片 幻灯片

1号路径线的模拟结果显示了巴颜喀拉块体的北部等效应力西部地区比东部大，大致从西向东逐渐递减； 2号路径线的模拟结果显示了巴颜喀拉块体内部等效应力西部大、 东部小，大致从西向东逐渐递减，但在块体的东边等效应力呈起伏变化； 3号与4号路径线的模拟结果显示了等效应力从巴颜喀拉块体南向北逐渐递减； 但4号比3号路径的等效应力低； 5号路径线的模拟结果显示巴颜喀拉块体的东部等效应力比较低． 图3中的单点等效应力变化较大是断层所在的位置．

2.2   强震序列相互作用的模拟结果

以青藏高原及邻区为背景建立壳单元有限元模型，从区域动力学角度着手，模拟澳大利亚—印度板块推挤对我国青藏块体内巴颜喀拉块体及边缘强震形势的影响． 杨立明(2009)对青藏高原中强地震动力学过程进行了研究，表明区域强震之间是相互影响相互作用的，一次强震的发生会引起区域应力场的调整，其发生、 发展及演变过程蕴含着大震孕育、 发展、 发生过程的重要信息，具有深刻的动力学内涵． 本文主要模拟巴颜喀拉块体内及周缘1900年以来发生的14次7级以上地震. 其震源参数见表2．

表  2  1900年以来巴颜喀拉块体内及周缘发生的M_S≥7.0地震目录

Table  2.  Earthquake catalogue(M_S≥ 7.0)in the Bayan Har block and surroundings since 1900

下载: 导出CSV 

| 显示表格

2.2.1   系列强震引起应力变化模拟结果

Lu等(2011)以及杨树新等(2012)在模拟过程中通过降低震源所在断层的杨氏模量来模拟震源区部分丧失承载能力引起的应力场的调整变化，本文模拟过程中也使用该方法.对于近年来发生地震的错动距离主要依据相关文献中给出的结果(陈兵等，2003;陈宇坤，陈杰，2004;陈运泰，2008;张桂芳等，2011)，发生较早的地震通过经验公式得到位错，在背景应力场下通过不断尝试使震源处由于降低刚度所产生的位移与地震水平错动距离基本一致.由于区域水平构造应力场大小､方向､各次地震震中区发震断层走向及地震震级大小不一致，上述14次地震震源处的弹性模量降低量也为3.5%—37%不等.从表2中可以看出，降低量最小的为亦基台错地震，最大的为汶川地震.在模拟过程中地震是按发生的时间顺序逐步累加的，用等效应力值的变化来分析地震之间的相互影响.

4次地震的模拟结果如图4a-n所示.图中的等效应力是模拟地震后各节点的等效应力减去背景场的等效应力，即所得到的等效应力为地震所引起的等效应力的变化.从图中可以清楚看出前面地震对后续地震的影响及未来相对危险的区域.

图  4  (a)-(h) 地震引起的等效应力变化图

(a)模拟第1次地震(道孚M_S7.0);(b)模拟第2次地震(炉霍M_S7.3);(c)模拟第3次地震(民丰M_S7.3);(d)模拟第4次地震(叠溪M_S7.5);(e)模拟第5次地震(阿兰湖M_S7.5);(f)模拟第6次地震(达日M_S7.7);(g)模拟第7次地震(康定M_S7.5);(h)模拟第8次地震(阿兰湖M_S7.0).图中虚线为0Pa线，图中黑色圆点表示下一次地震发生的位置

Figure  4.  (a)-(h) Equivalent stress variation caused by earthquakes

(a)The first simulated earthquake(Daofu M_S7 .0); (b)The second simulated earthquake (Luhuo M_S7 .3);(c)T he third sim ulated earthquake (Minfeng M_S7 .3); (d)T he fourth sim ulated earthquake (DiexiM_S7 .5); (e)T he fifth sim ulated earthquake (Alan Lake M_S7 .5); (f)T he sixth sim ulated earthquake(Dari M_S7 .7); (g)T he seventh sim ulated earthquake (Kangding M_S7 .5); (h)T he eighth sim ulatedearthquake(Alan Lake M_S7 .0). Dashed lines for 0 Pa and black point stand for the location of nextearthquake in the figure

下载: 全尺寸图片 幻灯片

图  4  (i)-(n) 地震引起的等效应力变化图

(i)模拟第9次地震(炉霍M_S7.6);(j)模拟第10次地震(亦基台错M_S7.3);(k)模拟第11次地震(玛尼M_S7.5);(l)模拟第12次地震(昆仑山口西M_S8.1);(m)模拟第13次地震(汶川M_S8.0);(n)模拟第14次地震(玉树M_S7.1).图中虚线为0Pa线，黑色圆点表示下一次地震发生的位置

Figure  4.  (i)-(n) Equivalent stress variation caused by earthquakes

(i)T he ninth sim ulated earthquake (Luhuo M_S7 .6); (j)T he tenth sim ulated earthquake (Yijitaicuo M_S7 .3);(k)T he eleventh sim ulated earthquake (M ani M_S7 .5); (l) T he T welfth sim ulated earthquake (K nlunM ountains Pass M_S8 .1); (m)T he thirteenth im ulated earthquake (W enchuan M_S8 .0); (n)T he fourteenthsim ulated earthquake (Y ushu M_S7 .1). Dashed lines for 0 Pa and black point stand for the location of nextearthquake in the figure

下载: 全尺寸图片 幻灯片

从模拟结果来看，有8次地震发生在等效应力增加的区域，分别为1923年炉霍M_S7.3，1924年民丰M_S7.3，1937年阿兰湖M_S7.5，1955年康定M_S7.5，1973年亦基台错M_S7.3，1997年玛尼M_S7.5，2001年昆仑山口西M_S8.1和2008年汶川M_S8.0地震;有3次发生在等效应力减小的区域，分别为1933年叠溪M_S7.5地震，1973年炉霍M_S7.6地震和2010年玉树M_S7.1地震;有2次发生在等效应力增减过渡的区域，分别为1947年达日M_S7.7地震和1963年阿兰湖M_S7.0.周仕勇(2008)通过定量计算断层破裂产生的库仑应力得出断层中发生的强震的相互影响，加载､卸载和影响很小的情况都存在.从本文的模拟结果来看这几种情况也都存在，总体来说后续强震主要发生在前面地震所产生的应力加载区.

2.2.2   历史地震对1970年后地震的影响

1970年至目前，巴颜喀拉块体周缘共发生了6次7级以上地震，分别为1973年炉霍M_S7.6，1973年亦基台错M_S7.3，1997年玛尼M_S7.5，2001年昆仑山口西M_S8.1，2008年汶川M_S8.0和2010年玉树M_S7.1地震.从图5中可以看出，炉霍地震前震源区由已发生的地震所引起的应力变化一直在下降，1955年康定M_S7.5地震使应力大幅增加，1963年阿兰湖M_S7.0地震对炉霍地震无影响(图5a);亦基台错、玛尼、昆仑山口西这3次地震前等效应力逐步增大，尤其是亦基台错离玛尼较近，亦基台错地震使玛尼的等效应力突增，也促使了这几次地震提前发生(图5b，c，d);叠溪、阿兰湖、达日、康定等地震使汶川震源处的等效应力大幅增加，1973年炉霍M_S7.6地震使其应力水平降低，汶川地震发生前应力水平不是在最大的时段(图5e);叠溪、阿兰湖、达日地震使玉树地震震源处的应力水平逐步降低，尤其是1947年达日M_S7.7地震引起的等效应力减小幅度最大(图5f)，昆仑山口西、汶川地震的发生也使玉树地震震源处的应力得到了释放，玉树地震是在应力水平降低过程中发生的，也就是说昆仑山口西、 汶川地震的发生可能延缓了玉树地震的发生．

图  5  历史地震(地震顺序见表2)对1970年以来地震的影响

(a) 1973年炉霍地震前震源处应力变化曲线; (b) 1973年亦基台错地震前震源处应力变化曲线;(c) 1997年玛尼地震前震源处应力变化曲线; (d) 2001年昆仑山口西地震前震源处应力变化曲线;(e) 2008年汶川地震前震源处应力变化曲线; (f) 2010年玉树地震前震源处应力变化曲线

Figure  5.  Influence of previous earthquakes on the earthquakes since 1970

(a) Stress change in the source region before 1973 Luhuo earthquake; (b) Stress change before 1973 Yijitaicuo earthquake; (c) Stress change before 1997 Mani earthquake; (d) Stress change before 2001 Kunlun Mountain Pass earthquake; (e) Stress change before 2008 Wenchuan earthquake; (f) Stress change before 2010 Yushu earthquake

下载: 全尺寸图片 幻灯片

3.   讨论与结论

1) 本文采用壳单元有限元模型，采用线弹性本构关系，在边界条件约束下，计算出的背景应力场特征大致为： 在青藏高原内部从西到东，从南到北等效应力逐渐减小．

2) 以青藏高原及邻区为背景建立有限元模型，模拟了青藏高原巴颜喀拉块体内及周缘强震之间的相互影响. 从模拟结果看，14次地震中有8次地震发生在等效应力增加的区域，3次发生在等效应力减小的区域，2次发生在等效应力增减过渡的区域． 说明强震之间的相互影响，加载、 卸载和影响很小的情况都存在，总体来说强震发生在应力增加区域的次数较多．

3) 模拟了历史地震对1970年后地震的影响. 结果表明，5次地震发生前震源区等效应力逐渐增大，其中4次震前震源区等效应力持续增加，分别为1973年炉霍M_S7.6，1973年亦基台错M_S7.3，1997年玛尼M_S7.5，2001年昆仑山口西M_S8.1地震，说明这几次地震前发生的强震对本次地震具有触发作用； 1次地震(2008年汶川M_S8.0地震)震前震源区等效应力略有增加，但变化幅度不大，说明汶川地震前的几次强震对汶川震源区影响较小； 1次地震(2010年玉树M_S7.1地震)震前震源区等效应力持续降低，也就是说昆仑山口西、 汶川等强震可能延缓了玉树地震的发生．

4) 本文中使用的是一种简化的模型 线弹性壳体模型，同时未考虑重力、 地壳介质的流变性等因素的影响． 等(程佳2011)计算了玛尼、 昆仑山口西、 汶川地震所形成的同震和震后形变场的变化过程与特征，本文结果与他们的同震结果是一致的，但是与黏弹性结果存在差别，显示了开展黏弹性数值模拟研究的重要性． 针对这些问题，作者正在开展三维黏弹性有限元的研究工作，以便把本研究推向深入．

5) 用等效应力来描述地震序列过程中的应力场变化，其物理意义是着重分析孕震区应力的积累与释放特征． 但是，由于不涉及地震类型，这种分析方法存在着一定的局限性．

本文是在中国地震局地壳应力研究所(地壳动力学重点实验室)完成，审稿专家提出了宝贵建议，在此一并表示感谢.