引言

我国大陆地区自1976年唐山M_S7.8地震发生后，经历了12年无M_S≥7.0地震的平静期。从1988年起进入了M_S≥7.0地震的活跃期，到2017年底已发生14次M_S≥7.0地震。这些大地震均发生在青藏高原地区地壳构造块体的边界断层上，表明现阶段青藏高原地壳构造块体系统中，在块体边界断层的多个部位上，应变已积累到或接近地下岩石的破裂强度水平。特别是，从1996年至今，我国大陆地区所有M_S≥7.0地震，包括2次M_S≥8.0地震（2001年昆仑山口西M_S8.1地震、2008年汶川M_S8.0地震）和7次M_S≥7.0地震（1996年喀喇昆仑山M_S7.1地震、1997年玛尼M_S7.9地震、2008年于田M_S7.0地震、2010年玉树M_S7.3地震、2013年芦山 M_S7.0地震、2014年于田M_S7.0地震、2017年九寨沟M_S7.0地震），均发生在巴颜喀拉块体的边界断层上。

图1是我国大陆地区构造块体及1988年以来M_S≥6.9地震震中分布图，其中，块体的划分引自张培震等（2003），震中分布图引自本研究团队于2009年发表的论文（陈祖安等, 2009），但增补了自2009年以来发生的M_S≥7.0地震。由图1可见，青藏高原地区各构造块体的边界断层上大震活动从未有过的活跃，显示了该地区的大震活动与构造块体运动、变形的密切关联。

图  1  我国大陆地区构造块体及1988年以来M_S≥6.9地震震中分布图

Figure  1.  Distribution of tectonic blocks and M_S≥6.9 earthquakes since 1988

下载: 全尺寸图片 幻灯片

2008年汶川M_S8.0地震发生前，在巴颜喀拉地块的西、南、北边界断层上，分别发生了多次M_S≥7.0地震，唯独在其东边界断层龙门山断裂带上未发生M_S≥7.0地震（陈祖安等，2009）。2001年发生在巴颜喀拉地块北边界断层东昆仑断裂带上的昆仑山口西M_S8.1地震，其震源机制是左旋走滑错动，破裂过程成像显示该地震的破裂总体上自西向东扩展（许力生，陈运泰，2004）。近东—西走向的昆仑山口西地震发震断层南侧向东的错动，必然会推动巴颜喀拉地块进一步东移，从而加速龙门山断裂带的应变积累。从大震活动与构造块体运动变形相互影响的角度分析，青藏高原构造块体边界断层的某些地段，特别是龙门山断裂带，的中长期大震活动危险性应予以特别关注。自汶川大地震前两年起，作者以国家自然科学基金项目“玛尼、昆仑山大震发生对青藏高原构造块体稳定性影响的数值模拟研究”为题，开始对青藏高原构造块体系统的稳定性、各构造块体边界断层上失稳危险度的分布，以及1997年玛尼地震和2001年昆仑山口西地震的发生对构造块体系统稳定性的影响等方面的问题进行了研究，获得了一些有意义的结果（陈祖安等，2008）。

2008年汶川大地震发生近5年之后，2013年4月20日在龙门山断裂带西南段上汶川发震断层的西南端外侧发生了芦山地震。由地震波资料反演结果可知，这两次大震均为全球内陆地区罕见的低倾角逆冲型大地震（陈运泰等，2013）。有关这两次大震的孕震机理及震前中长期地震危险性，是值得深入研究的问题。近年来，本文第一作者及其研究团队成员运用三维流变非连续变形法（描述地壳块体边界）与有限元法（描述地壳块体内部）相结合的数值模拟方法，研究了青藏高原及其邻区（东侧四川盆地和鄂尔多斯地块）内共12个构造块体的运动、变形、相互作用及其与近几十年发生于该区的大地震之间的关系（陈祖安等，2008，2009，2011；Chen et al，2011 ；林邦慧等，2012，2014，2016）。本文将在系统总结上述研究工作的基础上，进一步研究汶川地震和芦山地震的孕震机理以及大震前的中长期地震危险性。

2.   方法与模型

印度板块向北挤压造成青藏高原的隆起并产生厚约70 km的地壳。由于地壳内温度上升，下地壳层速度低且高导层发育，介质呈较强黏塑性。在印度板块向北的推挤下，下地壳的黏塑性流动驱动着上地壳脆性层运动与变形，形成了该地区复杂的孕震环境。考虑到上述因素，宜采用描述块体边界相互作用十分有效的非连续变形法（discontinuous deformation analysis，简写为DDA）与描述块体内相互作用的有限元法（finite element method，简写为FEM）相结合的方法研究这一问题。非连续变形分析法是一个分析块体系统中块体相互作用的有效工具（石根华，1997），可以用来计算块体系统在动力及静力作用下块体的滑移和不连续变形。对于各块体，允许有滑移、压性或剪切变形；对于整个块体系统，允许块体边界断层有滑动。迄今，已有一些研究人员使用该方法进行了多项研究，并取得重要成果（Cai et al，2000 ；白武明等，2003），这些结果证明了该方法适合用于模拟大震的破裂过程和大震活动与构造块体相互作用的关系。

陈祖安等（2011）在石根华（1997）非连续变形法和二维情形程序的基础上，编制了三维流变DDA方法与FEM方法相结合的程序“三维DDA＋FEM程序”。该程序的关键之处在于将块体之间的相互作用转变为它们之间的接触问题（Perić，Owen，1992）；同时假设位移很小，块体之间的接触点已知，详细情况请参见陈祖安等（2011），这里从略。本文将这个程序用于汶川、芦山地震孕震机理等研究。

研究区包括青藏高原及其东侧的四川盆地和鄂尔多斯地区（图2a）。参考张培震等（2003）对我国大陆地区构造块体的划分结果，将研究区划分为塔里木、柴达木、阿拉善、祁连、昆仑、羌塘、拉萨、川滇、滇西南、鄂尔多斯、四川盆地及华南等12个地块。地块边界断层包括阿尔金、祁连山、西秦岭、东昆仑、西昆仑、六盘山、贺兰山、玛尼—玉树、喀喇昆仑、喜马拉雅、龙门山、鲜水河、三江、红河、安宁河—小江、澜沧江、秦岭大别山、汾渭和阴山等断裂带。

将上述区域垂直深度达120 km的空间划分为5层，在青藏高原东缘以西地区划分上地壳20 km，中地壳20 km，下地壳30 km及上地幔。受资料限制，为简化起见，除龙门山断裂带以外，其它边界断层均为直立。在青藏高原东缘及以东地区，考虑了龙门山断裂带东西两侧地势、地壳厚度和分层的明显变化（滕吉文等，2008）：地势从四川盆地海拔500 m向西到龙门山地区陡升至4 000 m；地壳厚度在四川盆地为45 km，向西逐渐增厚，在龙门山断裂带东侧地区为60 km，并逐渐向西增厚，其西侧的青藏高原地区地壳厚度为70 km。

我们参考了陈运泰研究团队有关汶川地震破裂过程的反演结果（陈运泰等，2008，2013；张勇等，2008），对本研究团队2009年论文中所使用的汶川地震的发震断层模型（陈祖安等，2009）进行了改进。考虑到龙门山断裂带上汶川大震的发震断层由3个相连接的子事件组成，倾角由西南向东北扩展且逐渐增大，南段倾角在中地壳为40°，在上地壳中为60°，北段的倾角接近直立。参考万永革等（2009）一文中的龙门山断裂带三维模型，改进并构制了龙门山断裂带上汶川发震断层的三维分布模型（图2b）。

图  2  研究区构造块体模型图（a）和龙门山断裂带上汶川M_S8.0地震发震断层的三维示意图（b）（引自陈祖安等，2011）

Figure  2.  Tectonic block system model used in this study （a） and 3D schematic diagram of the seismogenic faults of the 2008 Wenchuan M_S8.0 earthquake on the Longmenshan fault zone （b） （after Chen et al，2011 ）

下载: 全尺寸图片 幻灯片

3.   研究区块体边界的失稳危险度、应变率强度及应变能密度分布

This page contains the following errors:

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

This page contains the following errors:

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

我们模拟计算了研究区现今构造块体边界断层上的失稳危险度分布（林邦慧等，2012，2014）。计算结果表明，上、中地壳构造块体边界断层上的失稳危险度高的地段多数与近30年来发生M_S≥7.0地震的区域基本一致，包括1990年共和M_S6.9地震、1997年玛尼M_S7.9地震、2001年昆仑山口西M_S8.1地震、2008年汶川M_S8.0地震、2013年芦山M_S7.0地震、2014年于田M_S7.3地震和2017年九寨沟M_S7.0地震的发震断层。此外，龙门山断裂带（包括汶川地震和芦山地震的发震断层）除北西端外均为失稳危险度最高值（图4）。本文参考陈运泰等（2008，2013）和张勇等（2008）关于汶川地震破裂过程的反演结果提出的汶川发震断层的三维模型的模拟结果表明，边界断层特别是龙门山断裂带的失稳危险度分布明显地比2009年用较简单的汶川发震断层模型（陈祖安等，2009）结果更符合近30年M_S≥7.0地震的地震活动性。

图  3  由数值模拟计算得到的研究区地表速度场（a）和上地壳层震源机制（b）分布图（引自陈祖安等，2009）

Figure  3.  Distributions of the ground surface velocity calculated （a） and the focal mechanism of upper crust （b） in the study area （after Chen et al，2009）

下载: 全尺寸图片 幻灯片

图  4  研究区上地壳（a）和中地壳（b）构造块体边界断层上的失稳危险度分布

1. 1990年共和M_S6.9地震；2. 1997年玛尼M_S7.9地震；3. 2001年昆仑山口西M_S8.1地震；4. 2008年汶川M_S8.0地震；5. 2013年芦山M_S7.0地震；6. 2014年于田M_S7.3地震；7. 2017年九寨沟M_S7.0地震

Figure  4.  Distribution of instability risk factor of boundary faults between tectonic blocks of upper crust （a） and lower crust （b） in the studied area

1. M_S6.9 Gonghe earthquake in 1990；2. M_S7.9 Mani earthquake in 1997；3. M_S8.1 western Kunlun Mountain Pass earthquake in 2001；4. M_S8.0 Wenchuan earthquake in 2008；5. M_S7.0 Lushan earthquake in 2013；6. M_S7.3 Yutian earthquake in 2014；7. M_S7.0 Jiuzhaigou earthquake in 2017

下载: 全尺寸图片 幻灯片

为了进一步研究汶川大震的孕震机理，我们还计算了研究区上、中地壳的应变率强度分布，结果如图5所示。计算结果表明，青藏高原绝大部分地区均为高剪切应变率地区，其东缘基本上为应变率强度的急剧变化带，特别是龙门山断裂带东西两侧应变率强度的变化最为急剧，其西侧的应变率强度约是东侧的4倍。我们还计算了研究区应变能密度分布（陈祖安等，2009），结果如图6所示，可以看出：龙门山断裂带及邻近地区为宽度与龙门山断裂带宽度基本相同、走向一致的高应变能密度带，其东西两侧则是应变能密度较低的地区，且断裂带东侧应变能密度较西侧衰减得快。

图  5  研究区上地壳（上）和中地壳（下）应变率强度分布图（左）及其等值线图（右）（引自陈祖安等，2009）

Figure  5.  Distribution of strain-rate intensity （left） and contours of strain-rate （right） of the upper crust layer （upper） and the middle crust layer （lower） in the studied area （after Chen et al，2009）

下载: 全尺寸图片 幻灯片

图  6  研究区上地壳（上）、中地壳（下）应变能密度分布图（左）及其等值线图（右）（引自陈祖安等，2009）

Figure  6.  Distribution of strain energy-density （left） and contours of strain energy-density （right） of the upper crust layer （upper） and the middle crust layer （lower） in the studied area （afterChen et al，2009 ）

下载: 全尺寸图片 幻灯片

4.   2001年昆仑山口西地震和2008年汶川地震的影响

4.1   2001年昆仑山口西M_S8.1地震的影响

为了探讨2001年昆仑山口西M_S8.1地震对2008年汶川M_S8.0地震产生的影响，我们通过数值计算模拟了昆仑山口西地震的破裂过程，研究了该地震对研究区各块体边界断层库仑破裂应力分布的影响（陈祖安等，2011）。图7给出了昆仑山口西地震引起的研究区各边界断层上库仑破裂应力变化的分布，可见，该地震的发生除了引起其发震断层两端部库仑破裂应力增加外，对阿尔金断裂带、玛尼—玉树断裂带中段和西南地段，以及从龙门山断裂带的北端到南端，包括汶川地震和芦山地震的发震断层，库仑破裂应力均有较大的增加，增幅约为0.016 MPa。汶川地震和芦山地震的发震断层是同时具有高失稳危险度（林邦慧等，2012）且大震引起库仑破裂应力增加这两个特点的地带。计算得出，这两条断层上约0.016 MPa的库仑破裂应力增加相当于龙门山断裂带约两年的应力积累，也就是说，昆仑山口西地震的发生使汶川、芦山地震的发震断层失稳破裂提前了大约两年。

图  7  昆仑山口西M_S8.1地震对青藏高原及邻近地区上地壳（a）和中地壳（b）各块体边界断层上库仑破裂应力的影响（引自陈祖安等，2011）

Figure  7.  The influence of the 2001 western Kunlun Mountain Pass M_S8.1 earthquake on the Coulomb failure stress of the boundary faults of the tectonic blocks of upper crust （a） and middle crust （b） in the Qinghai-Tibet Plateau and its vicinity （after Chen et al，2011）

下载: 全尺寸图片 幻灯片

4.2   汶川地震的影响及2013年芦山M_S7.0地震的成功预测

汶川地震的发生对青藏高原及邻近地区各边界断层地震危险性影响是本文特别予以关注的问题。汶川地震发生后，我们根据上述汶川大震破裂过程由3个倾角不同的子事件组成的研究结果，提出了汶川地震的发震断层模型（图2b），数值模拟了汶川地震的破裂过程，研究了大震引起各边界断层应力状态变化的特征（陈祖安等，2011），并在其文中的图6加上2013年M_S7.0芦山地震和2017年M_S7.0九寨沟地震的位置，如图8所示。结果显示，在上地壳，汶川发震断层南端外侧到龙门山断裂带南端、汶川发震断层北端外侧、鲜水河断裂带南段、东昆仑断裂带东南段等地段的库仑破裂应力增加幅度最大，达到约0.2 MPa，而2013年芦山地震的发震断层正处于汶川发震断层南端外侧，即龙门山断裂带西南段的中间部位，龙门山断裂带西南段是满足失稳危险度大（林邦慧等，2012）且经历昆仑山口大震和汶川大震两次库仑破裂应力增幅较大这两个条件的地区，是大震发生的危险区，林邦慧等（2012）据此对这次地震成功地进行了中长期地震预测。

图  8  汶川地震的发生对青藏高原及邻近地区上地壳（a）和中地壳（b）各块体边界断层库仑破裂应力的影响（修改自陈祖安等，2011）

Figure  8.  The influence of the 2008 Wenchuan M_S8.0 earthquake on the Coulomb failure stress of the boundary faults of the tectonic blocks of upper crust （a） and middle crust （b） in the Qinghai-Tibet Plateau and its vicinity （modified from Chen et al，2011 ）

下载: 全尺寸图片 幻灯片

5.   讨论与结论

本文概要叙述了作者运用三维流变非连续变形分析法与有限元法相结合的方法研究青藏高原及其东侧邻区构造块体的运动、变形、相互作用及其与近30年来发生于该区的大地震活动关系以及汶川、芦山大震的孕震机理，得出：

1） 失稳危险度高的地段与近30年来发生的M_S≥7.0地震所在区域基本一致，其中龙门山断裂带是力学上失稳危险度最高的地段，也是应变率强度的急剧变化带，西侧应变率强度约是东侧应变率强度的4倍，而且断裂带东侧的应变率强度等值线衰减得比西侧快。

2） 龙门山断裂带在上、中地壳中均位于与龙门山断裂带宽度相同、走向一致的高应变能密度带中，在上地壳的这个带的东西两侧则是应变能密度较低的地区，而在中地壳其强度在断裂带东侧逐渐向东衰减；西侧应变能密度高，而东侧应变能密度较低。这表明，汶川、芦山大震前，其发震断层已积累了大量的应变能，失稳危险度高，处于力学上的不稳定状态。

3） 2001年昆仑山口西M_S8.1地震的发生使得低倾角的汶川发震断层和包含芦山发震断层在内的龙门山断裂带西南段的库仑破裂应力增大约0.016 MPa，相当于龙门山断裂带约两年的应力积累，也就是说，昆仑山口西M_S8.1地震的发生使汶川、芦山地震发震断层的失稳破裂大约提前了两年。

4） 2008年汶川M_S8.0地震的发生也使芦山大震的发震断层和龙门山断裂带西南段的库仑破裂应力在昆仑山口西大震时增大，之后再次增大，应变能积累增强，对已处于失稳危险度较高状态的芦山大震发震断层的提前失稳起到了促进作用，据此成功地作了中长期地震预测（陈祖安等，2011；林邦慧等，2012）。陈运泰等（2013）的研究结果表明，龙门山断裂带西南段上，位于芦山大震发震断层南端至龙门山断裂带南端之间，以及芦山发震断层北端至汶川发震断层南端之间存在两个有发生大震危险的地震矩释放“亏空区”，长度分别约为70 km和30 km。在上述研究结果中，这两个区正好是满足失稳危险度大和库仑破裂应力多次增大这两个条件的地区。从青藏高原地区大震活动与构造块体运动变形的相互影响角度来看，汶川、芦山大震后龙门山断裂带上这些相关地段仍存在明显的中长期地震危险性。

5） 2008年汶川地震的发生引起东昆仑断裂带东南段库仑破裂应力的增幅也最大，表明汶川地震对于2017年九寨沟大震发震断层的提前失稳起到了促进作用（陈祖安等，2011；林邦慧等， 2012，2016）。

作者自2001年昆仑山口西M_S8.1地震发生以来，便开始关注巴颜喀拉块体边界断层上发生的大震活动。2006年启动的国家自然科学基金项目“玛尼、昆仑山大震发生对青藏高原构造块体稳定性影响的数值模拟研究”，使得我们对各构造块体边界断层上失稳危险度分布以及玛尼、昆仑山口西地震的发生对构造块体系统稳定性的影响等问题获得了一些新的认识。从2008年汶川地震至今近十年来，我们对发生在巴颜喀拉块体边界断层上的一系列地震（包括1990年共和M_S6.9地震、1997年玛尼M_S7.9地震、2001年昆仑山口西M_S8.1地震、2008年汶川M_S8.0地震、2013年芦山M_S7.0地震、2014年于田M_S7.3地震以及2017年九寨沟M_S7.0地震）的研究结果表明了这些大地震的发生与高失稳危险度、高应变能密度分布以及高应变率强度分布的急剧变化密切关联，这些新的认识对于中长期地震危险性的研究和地震预测预报研究是有益的参考。

作者对陈运泰院士以及两位匿名审稿专家对本文的修改建议表示衷心的感谢。