2012, Vol. 55
2. 中国科学院研究生院, 北京 100049;
3. 中国科学技术大学地球与空间科学学院, 合肥 230026;
4. 中国地震局地球物理研究所, 北京 100081
, LI Jia-Hang3, HAN Li-Bo4, FANG Li-Hua4, YANG Song1,2, JIN Bi-Kai1,2, ZHENG Yong1, XIONG Xiong1
2. Graduate University, Chinese Academy of Sciences, Beijing 100049, China;
3. School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China;
4. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
2010年4 月14 日7 时49 分青海省玉树藏族自治州玉树县发生MS7.1(Mw6.8)级地震.根据中国地震台网中心测定结果，玉树地震震中位置为33.2°N，96.6°E，震源深度约14km，由于地震震中距离玉树州州府所在地的结古镇仅30km，地震造成了巨大的经济损失和人员伤亡[1-2].玉树地震发生在青藏高原巴颜格拉块体的中南部，发震断层为NWW-NW 走向的甘孜-玉树断裂.该断裂以左旋走滑运动为主，总长约500km，与鲜水河断裂在甘孜附近呈左阶错列，并与鲜水河断裂共同组成巴颜格拉块体的南边界[3].
和传统的板内地震相比，玉树MS7.1级地震的余震具有数量少､震级大的特点[4].截止到2010 年5月29日，仅记录到21次MS≥3.0级余震，但是记录到余震的最大震级达到MS6.3级，主震之后一个半月仍然发生了MS5.7级余震.研究玉树地震主震与余震之间的关系，对于我们了解余震的发震机理具有十分重要的参考价值.
近年来，越来越多的研究表明，地震的发生释放了断层上所积累的应力，但应力并不会完全凭空消失，部分应力转移到其他地区，导致应力聚集，进而诱发后续地震，这就是地震的应力触发理论[5].根据应力触发理论，库仑应力的增强相当于断层额外负荷的加载，促进地震的发生；反之，应力影区的库仑应力减弱相当于断层负荷部分卸载，其发生地震的可能性降低[6].应力触发理论在一定程度上揭示了余震的分布规律[7-8]，Deng和Sykes[9-10]的研究发现85%的余震发生在震后应力增强区.但是，另外的研究发现，某些地震的余震分布并不完全遵循地震触发规律，例如:1994 年Northridge地震之后的余震只有60%位于库仑应力增长区域[11]；台湾集集地震之后的两天内位于应力影区断层上的余震活动却十分活跃[12]；Marsan[13]发现LomaPrieta､Landers､Northridge等地震之后的100天内没有观测到平静的应力影区.玉树MS7.1级地震为我们研究主震对余震的触发作用提供了一次很好的震例.
此外，2011 年6 月17 日在囊谦县城以西发生了MS5.2级囊谦地震，中国地震台网中心给出的震中位置为32.4°N，95.9°E.由于地震发生的区域人口不密集，没有造成较大的人员伤亡.Shan 等[4]根据玉树地震所导致的同震应力变化场以及该地区的背景地震活动性认为囊谦地震受到了玉树地震的触发作用.但是在他们的研究中存在两个问题:第一，库仑应力变化的计算过程中始终选取玉树地震的震源机制解作为接收断层参数，没有考虑囊谦地震真实的震源信息，因此，计算得到的囊谦地震震中的库仑应力变化与实际情况可能存在误差；第二，基于弹性地壳模型，只计算了玉树地震所导致的静态库仑应力变化，没有考虑两次地震之间一年零两个月内的震后效应的影响.
本文利用弹性位错理论和分层岩石圈模型，计算玉树地震引起的同震及震后黏弹松弛应力场变化，讨论MS7.1级玉树地震对余震分布的影响以及与2011年囊谦MS5.2 级地震之间的触发关系，为深入研究地震的应力触发理论提供参考依据.
2 原理与模型 2.1 震后库仑应力变化的计算当断层上所积累的剪应力大到足以克服与正应力和摩擦系数有关的静摩擦力时，锁闭的断层就会发生破裂进而引发地震，这一地震发生的过程可以用库仑破裂准则来描述[14-15].虽然，地壳中构造应力的绝对值难以通过测量的手段直接获得，但是大地震之后地壳中的应力场变化却可以通过弹性位错理论来进行很好地估算[6].
利用弹性位错理论，可以计算出由静态地震滑移所引起的应力场变化[16]，投影到相邻的接收断层上就可以得到相应的正应力和剪应力变化，而根据库仑破裂定律
|
(1) |
|
图 1 玉树及周边地区构造背景 沙滩球为2010年玉树MS7.1级地震和2011年MS5.2级囊谦地震震源机制解，灰色圆圈为玉树地震之后2010年4月13日至6月17日间870次ML>1.0级余震，白色圆圈为过去40年内背景地震活动，白色五角星为该地区内主要城镇，黑色线段为主要活动断层. Fig. 1 Sketch map of tectonics in Yushu and its surrounding region Beachballs represent the focal mechanism of2010 MS=7.1 Yushu earthquake and 2011 MS=5.2 Nangqên earthquake .Gray circles represent the 870 after shocks withML>1.0 occurred between 13 April and 17 June，2010.White circles represen this torical earthquakes in past 40years.Whites stars represent maincities and countie sin this region.Blacklines represent faults. |
可计算出所研究断层面上库仑应力变化的大小.式中Δσf为计算断层面上库仑应力变化，Δτs 是断层面上剪应力的变化(与断层滑动方向相同时为正值)，μ′为有效摩擦系数，Δσn 是计算断层面上正应力变化(压应力方向为正值).Parsons等[17]认为，积累较大滑移量的断层等效摩擦系数较小(0.2~0.4)，而积累较小滑移量的断层等效摩擦系数则比较大(0.6~0.8).King等[7]认为有效摩擦系数对库仑应力变化的计算结果将产生中等大小的影响，因而在计算过程中可以选取0.4 的中间值，这一做法通常为大多数地球物理工作者所采用[18-19].在本文计算中，选取有效摩擦系数0.4，并在结果讨论中比较不同摩擦系数取值对计算结果的影响.
2.2 计算模型 2.2.1 Ms7.1级玉树地震有限断层模型随着大地测量技术的发展，越来越多的大地测量资料，例如GPS､InSAR 和SAR 等，具有了很高的空间覆盖率，为我们研究大地震引起的近场同震及震后形变提供了良好的约束条件.Shan等[4]利用玉树地震InSAR 资料反演了玉树地震的同震动态破裂过程.虽然，Shan 等[4]给出的玉树地震有限断层模型在反演过程中进行了一定简化，忽略了断层破裂面在西北端的弯曲，但是由于玉树地震的同震位错分布主要位于破裂面的东南段，因此，模型的简化对本研究中库仑应力计算结果的影响较小.有限断层反演结果显示，玉树地震的破裂过程以左旋走滑运动为主，最大位错量约为1.5m，位于玉树地震震中东南方向20~30km 处近地表区域，该结果与实际地表考察结果较为一致[3].
2.2.2 地壳分层模型Wang[20]提出了利用正交归一法计算地震应力场Green函数的方法，并在此基础上建立了黏弹松弛分层模型下的地震同震及震后形变模型，并发展了相应的数值方法.运用PSGRN/PSCMP 软件[21]，根据同震破裂模型[4]，我们计算了玉树MS7.1级地震导致同震及震后黏弹松弛应力场变化.分层地壳模型参数取自Crust2.0，参数如表 1 所示.根据公式(2)，分层模型各层的剪切模量可以通过密度和地震P波波速求得[22].
|
(2) |
分层地壳模型参数的不确定性会对模型计算结果产生一定影响.本研究根据Crust2.0 模型，给出的地壳厚度为70km，与研究区的实际结构可能会存在一定的差异.不同研究给出的地壳厚度也不完全一样，例如，王有学等[23]根据爆破测深资料的研究结果表明，金沙江断裂附近的玉树一带莫霍面深达78km；Chen等[24]用接收函数研究中国大陆地壳结构的计算结果表明这一地区的地壳厚度为65~70km；王未来等[25]根据玉树台接收函数h-k叠加计算的地壳厚度在70±2km 左右.因此，在表 1模型参数的基础上，我们改变地壳模型各层厚度±10%，研究了地壳分层模型参数不确定性对同震库仑应力变化计算结果的影响.结果表明，当模型中地壳厚度从63km 增加到77km 时，同震库仑应力变化场的空间分布基本没有发生改变，地壳模型参数变化的影响主要位于玉树地震的主震破裂面附近，并且远小于该区域同震库仑应力变化场的大小，所以本研究中采用Crust2.0 模型计算得到的结果是可以接受的.
|
表 1 分层地壳模型及参数 Table 1 Parameters of the multilayer crustal model |
我们假设玉树及周边地区震后黏弹松弛的过程主要发生在22km 的上地壳以下.震后黏弹松弛应力变化是指在一定深度下，受高温高压环境的影响，岩石层的力学性质逐渐由脆性向黏塑性转化，但由于地震发生时破裂的速度很快，在短时间内黏弹性的中下地壳以及上地幔的力学响应表现为弹性体，然而震后随着时间的推移，中下地壳以及上地幔所积累的应力､应变逐渐释放，并且向上传输到弹性的上地壳中，导致上地壳中应力状态发生改变[26].虽然现在很多研究认为标准线弹性体[27]､Burgers体[28]等复杂的岩石圈流变模型能够更好地解释震后黏弹松弛过程，但由于玉树及周边地区缺乏长期的形变观测资料和相关的震后形变研究成果，模型的复杂化将会带来更多的参数不确定性，因此在本研究中还是采用最简单和常用的Maxwell流变模型.由于玉树地震与囊谦地震均发生在羌塘块体上，根据邵志刚等[28]数值模拟的研究结果，计算取中､下地壳黏滞系数为5×1017Pa·s，上地幔黏滞系数取1×1020Pa·s，上地壳为完全弹性介质，黏滞系数取无穷大.文章的讨论部分将分析黏滞系数的取值对震后黏弹松弛应力变化计算结果的影响.
2.2.3 囊谦地震震源机制反演在Ms7.1级玉树地震同震及震后应力场变化的计算结果上，研究2010 年玉树地震与2011 年Ms5.2级囊谦地震之间的触发关系，需要了解囊谦地震相关的震源信息.本文采用国际上近年发展起来的剪切粘贴(CAP)方法对囊谦地震震源机制解进行了反演[29-30]，它利用近震地震数据，将宽频带数字地震记录分为体波部分(Pnl)和面波部分，分别计算它们的合成波形和真实记录的误差函数，在相关参数空间中搜索最佳解进行矩张量反演，目前在国内应用也较为广泛[31-34].本研究使用经震中距矫正后的绝对误差值作为误差函数[30]，定义为
|
(3) |
式中，r为震中距，r0 为选定的参考震中距，p为指数因子，一般而言，对体波p=1，面波p=0.5.u为理论合成地震波形，s为观测地震波形.
反演过程中使用了国家数字测震台网数据备份中心[35]提供的国家台网及青海､四川､西藏等区域台网的宽频带数字地震波形记录.按照方位角覆盖及信噪比的要求，挑选出符合要求的数据进行处理，对挑选出的宽频带数据首先去除仪器响应，并旋转至大圆路径，然后对波形中Pnl部分使用带宽为0.02~0.2Hz､面波部分使用带宽为0.01~0.1 Hz的带通滤波器进行滤波.相应地，采用相同的滤波参数，对计算得到的理论地震图进行滤波.计算理论地震图时，采用了目前广泛使用的频率-波数(F-K)法[36]，该方法适用于水平分层的简单地壳速度结构模型.
根据公式(3)给出的误差目标函数，在参数全空间范围内搜索最佳震源机制解､震源深度和矩震级.图 3a给出了观测波形与理论波形拟合误差随深度的分布，从图中可以看出反演收敛较好，误差函数在深度为9km 时达到极小，因此对应的震源深度9km即为反演所得最佳震源深度.而在最佳震源深度附近，震源机制解变化并不大，这也在一定程度上说明反演过程中震源机制解较为稳定.反演得到囊谦地震震源机制解的两个节面分别为:节面I，strike=259.8°，dip=84.4°，rake= -21.1°；节面Ⅱ，strike=352°，dip=69°，rake= -174°.图 3b 为波形拟合情况，大部分台站面波部分波形拟合较好，说明震源机制解较为可靠；而体波部分拟合效果并不理想，可能是所使用的地壳模型和实际速度结构之间的差异所致.
|
图 3 2011年6月26日囊谦MS5.2地震震源深度及震源机制解反演结果 (a)矩张量反演中波形拟合误差随深度变化；(b)囊谦地震矩张量反演理论地震波形(灰色)与实际观测地震波形(黑色)比较，波形图下方第一行数字为理论地震波形相对实际观测波形的移动时间，正值表示理论波形相对观测波形向左移动.第二行数字为理论波形与观测波形的相关系数(百分比). Fig. 3 Results of the focal depth and mechanism of 26 June 2011Ms=5.2 Nangqên earthquake (a)Misfit variation with focal depth；(b)Comparison between synthetic (gray) and observed (black) seismograms.The numbers on the lower left side of these ismograms are the time shifts (upper)and cross-correlation coefficient in percent(lower).Positive time shifts mean the synthetics have been delayed or shifted to the left. |
根据Ms7.1级玉树地震同震破裂模型以及岩石圈分层模型，运用PSGRN/PSCMP软件[21]，我们计算了玉树地震导致的同震库仑应力变化场.本研究中用到的玉树地震余震重定位的结果，时间段从2010年4月13日至6月17日，共870个地震，震级为ML>1.0 级[25].由于玉树地区地震台站较为稀疏，为了提高定位精度，王未来等[25]在采用固定台震相数据的基础上，还使用了玉树附近7 个流动台的震相数据.定位的过程分两步，首先使用单纯形法[37]进行绝对定位，然后再使用双差方法进行定位，定位时使用的速度模型参考了姜枚等[38]的研究结果.经度､纬度和深度方向的平均定位误差分别为0.15km､0.18km 和0.3km，平均走时残差为0.1s[25].由于2010年6月17日之后的余震数量较少，震级也较小，在本研究中就没包括6 月17 日以后的余震.从图 1可以看到，玉树地震之后2个月内的余震主要分布在沿破裂面附近.在此基础上，我们比较了计算过程中选取不同参数(有效摩擦系数取值､计算时的参考深度､是否考虑区域构造应力场)，库仑应力变化场与玉树地震余震分布之间的关系.
3.1.1 有效摩擦系数取值对结果的影响图 4为有效摩擦系数不同取值的情况下(0.0､0.2､0.4 和0.8)不考虑区域背景构造应力场时，10km深度上库仑应力变化图像，计算中接收断层参数的选取参考玉树地震震源机制解.图中，2010年MS7.1级玉树地震所导致的库仑应力变化场主要包括四个库仑应力增强的扇区，其中两处扇区位于主震破裂面东南和西北两端的相邻区域，另外两处扇区位于主震破裂面的东北和西南区域.图 4a为有效摩擦系数取0.0 时的库仑应力变化场，根据公式(1)有效摩擦系数为0的时候，正应力在库仑应力变化计算中所占的比重为0，此时的库仑应力变化场即为剪应力变化场.随着有效摩擦系数取值的增加，正应力变化在库仑应力计算中所占的比重逐渐升高.由于玉树地震的破裂过程以左旋走滑运动为主，正应力的变化主要位于破裂面两侧呈近似对称的分布[7]，因此，摩擦系数的变化主要影响位于主震破裂面两侧的东北､西南两个应力增强区，有效摩擦系数越大，两处扇区不仅面积逐渐增大，并且增长幅度显著增强.
|
图 4 不同有效摩擦系数取值下的10km 深度上的库仑应力变化场 各图中右上角为相应的有效摩擦系数取值.本图及以后图中，灰色圆点为2010年4月13日至6月17日间870次ML>1.0级余震，黑线为该地区主要构造断层，灰色五角星为该地区主要城镇(NQ:囊谦；JG:结古镇；CD:称多). Fig. 4 Coseismic Coulomb failure stress changes with different effective coefficient of friction on the depth of 10km |
本研究用到的玉树地震震后两个月内870 次ML>1.0级地震主要分布于沿主震破裂面附近以及位于破裂面东北端应力增强扇区的区域内.我们统计了870次地震的震中分布与库仑应力变化场之间的关系(表 2)，结果表明，随着有效摩擦系数的变化，位于库仑应力增强区内的余震数量占总余震数目的53.2%~55.7%.在模型计算中，有效摩擦系数的取值对同震库仑应力变化场与余震分布之间的关系影响并不大，这主要是由于库仑应力变化场中受有效摩擦系数影响较大的东北､西南扇区内几乎没有余震分布.
|
|
表 2 不同参数下库仑应力变化场与余震的关系 Table 2 The number on the upper right side in each figure represents the corresponding effective coefficient of friction in its numerical calculation.In current and following figures，black lines represent faults.Gray circles represent the 870 after shocks withML>1.0 occurred between 13 April and 17June，2010.Gray stars represent maincities and counties in this region (NQ，JG and CD represent Nangqên County，Jiegu Town and Chengduo County，respectively). |
Shan等[4]对过去40年间玉树及周边地区记录到的历史地震活动进行了时间和空间上平滑处理，得到该地区的地震背景活动率，单位为每10年100km2区域内地震发生的次数.从玉树及周边地区的地震活动性来看[4]，玉树Ms7.1级地震发生之前结古镇西北-西南扇形区域内的背景地震活动较其他区域要活跃得多.Dieterich[39]和Catalli等[40]根据岩石物理学实验以及实际观测结果认为大地震之后的地震活动速率由该地区的背景地震活动性以及应力的突然加载或卸载所共同作用.因此，我们认为玉树及周边地区地震活动背景场的横向非均匀性很可能是导致不同应力增强区内的余震发生率差异较大的主要原因.
3.1.2 是否考虑区域构造应力场对结果的影响不考虑区域构造应力场时，库仑应力计算中选取玉树地震的震源机制解作为接收断层参数.然而，许多地震应力触发的研究结果[7, 41]表明，考虑背景区域构造应力场情况下，位于库仑应力增强区域内的余震数量会有显著提升.考虑到大地震的发生通常就体现该区域的构造背景，因此在PSGRN/PSCMP软件的计算过程中，把玉树地震的震源机制当作背景构造应力场的最优化破裂方向，背景构造应力场所导致最优化破裂方向即为我们给定的接收断层参数.考虑背景构造应力场时，地壳中的绝对应力状态由背景构造应力场和地震所导致的应力场两部分组成.基于地壳中的绝对应力状态，根据库仑破裂准则计算出地震发生以后新的最优化破裂方向.由于地震导致的静态库仑应力场变化大小随着与破裂面之间的距离增加衰减较快，考虑背景区域构造应力场时，在远场区域，地震导致的库仑应力场变化远远小于背景构造应力场，库仑应力计算中作为接收断层参数的最优化破裂面方向接近于玉树地震主震的震源机制；而在近场区域，同震静态应力场变化较大，地震导致最优化破裂面方向发生了显著变化，对库仑应力场计算结果影响较大.
图 5a和图 5b分别为不考虑和考虑背景构造应力场情况下计算得到的同震库仑应力场变化，图 5c为两者之间的残差.和前面分析的一样，是否考虑背景构造应力场对库仑应力计算结果的影响主要位于近场区域，尤其是主震破裂面附近.考虑背景构造应力场之后，原来沿断层面附近的库仑应力影区变成了库仑应力显著增强的区域，而该区域内余震分布较为密集，位于库仑应力增强区域内的余震数目与不考虑背景应力场相比从479次增加到754次，占余震总数的86.7%(表 2).
|
图 5 比较是否考虑区域构造应力场影响下10km 深度上的库仑应力变化场 (a)不考虑构造应力场；(b)考虑构造应力场；(c)考虑区域构造应力时的库仑应力变化场与不考虑情况之间的残差. Fig. 5 Coseismic Coulomb failure stress changes whether considering the influence of regional stress or not (a) Ignoring the influence of regional stress；(b)Taking the influence of regional stress in to account；(c)Misfit between whether considering the influence of regional stress or not. |
在之前的计算中，我们假设玉树地震的震源机制为背景构造应力场所导致的最优化破裂方向，但是实际情况中断层的滑动方向并不一定是实际背景构造应力场的最优化破裂方向，两者之间可能存在一定差距，而接收断层方向的不确定性对库仑应力变化场的计算结果影响较大[42].玉树MS7.1 级地震震源机制解对应的背景构造应力场P 轴方向为北偏东74°，王辉等[43]应用地震学的方法研究了中国大陆活动地块的应力应变场，推断该地区平均应变状态P轴方位角北偏东38.1°､倾角8.7°，而根据GPS观测得到的P 轴方向为北偏东24.9°[43].我们计算并比较了考虑最优化破裂面情况下区域构造应力场方向不确定性对库仑应力变化场的影响.结果表明，随着区域构造应力场方向的不同，应力场分布图像发生较大改变.当接收断层参数所对应的P 轴方向从北偏东74°向北偏东24.9°逆时针旋转时，库仑应力变化场图像也相应地产生了逆时针的旋转，位于主震破裂面东北端及沿破裂面的库仑应力增强区面积增大，库仑应力变化场与余震分布之间的相关性增强，基于地震学和GPS研究给出的背景构造应力场的计算结果能够更好地解释余震的分布，但总体而言对之前的结论没有影响.
3.1.3 不同深度上的库仑应力变化场对结果的影响在库仑应力变化场与余震分布之间相关性的研究中，当近场台站覆盖较好时能够获得较为可靠的余震震源深度定位结果，考虑不同深度上的余震分布与相应深度上库仑应力变化场之间的关系，能够使得结果更为可靠[44].玉树及周边地区地处青藏高原，地震台站分布稀疏，余震震源深度约束较差.因此，取有效摩擦系数0.4，不考虑背景构造应力场，我们计算了不同深度截面上(5.0､10.0､15.0 和20.0km)的库仑应力场变化，讨论不同深度应力场变化与全部870次余震的关系，如图 6所示.
|
图 6 有效摩擦系数取0.4时，不同深度上的库仑应力变化场 各图中右上角为相应的计算深度. Fig. 6 Coseismic Coulomb failure stress changes on different depths with the effective friction coefficient of 0.4 The number on the upper right side in each figure represents the corresponding calculation depth. |
由于玉树地震的破裂过程以左旋走滑运动为主，因此计算得到的库仑应力变化场随着深度的增加变化并不明显.由表 2可见，5~15km 的深度，位于库仑应力增强区的余震随着深度增加数量有一定的提高，占余震总数的百分比从46.8% 上升为59.7%.但是，从20km 深度上的库仑应力变化场来看，沿断层破裂面的周边区域均为库仑应力的增强区，93.9%的余震均位于库仑应力变化的增强区.这主要是由于玉树地震破裂模型(图 2)中的位错分布集中位于15km 深度以上，从走滑型地震应力场变化的垂直剖面图来看[7]，断层破裂面正下方及周边区域为库仑应力增强区，因此玉树地震在20km深度上的断层面附近以库仑应力增强区为主.
3.2 玉树地震对2011 年MS5.2 级囊谦地震的触发作用
2011年6月17日囊谦MS5.2级地震发生在玉树MS7.1级地震的震后一年左右，地震位于囊谦县城的西北西方向，震中位置32.4°N，95.9°E.从图 4-图 6中的同震库仑应力变化场图像来看，囊谦地震位于玉树地震所导致的破裂面西南方向的库仑应力增强扇区.这是否说明囊谦地震受玉树地震的影响，是一次被触发的地震事件?我们根据反演得到的囊谦地震震源机制解以及震源深度，研究玉树地震对囊谦地震的同震及震后黏弹松弛应力作用.
3.2.1 接收断层参数取不同节面时对计算结果的影响从反演结果来看，囊谦地震震源深度为9km，地震以走滑运动为主，震源机制解的两个节面分别为:节面I，strike=259.8°，dip=84.4°，rake=-21.1°；节面Ⅱ，strike=352°，dip=69°，rake=-174°，与玉树地震的震源机制解存在一定差异.由于囊谦地震的震级相对较小，震中附近没有明显的活动断层，对地震孕育发生的构造背景缺乏清晰的认识，无法区分囊谦地震震源机制解的两个节面中哪一个是实际的地震发震面，因此我们分别选取囊谦地震震源机制解的两个节面作为库仑应力计算中的接收断层参数，计算了玉树地震导致的同震及震后一年零两个月黏弹松弛库仑应力变化场.因为囊谦地震距离玉树地震破裂面较远，根据前面的分析，是否考虑背景区域构造应力场在囊谦地震震中附近的影响很小，所以在后面的库仑应力计算中均不考虑背景构造应力场的影响.根据囊谦地震的震源深度，计算选取的深度截面为9km，有效摩擦系数取中间值0.4，结果如图 7所示.
|
图 7 有效摩擦系数取0.4时，9km 深度上，接收断层参数选取囊谦地震两个不同节面计算的同震及震后库仑应力变化场(节面I:strike=259.8°，dip=84.4°，rake=-21.1°；节面Ⅱ:strike=352°，dip=69°，rake=-174°) (a)取节面I时的同震应力变化场；(b)取节面I时14个月之后的应力场变化；(c)取节面Ⅱ时的同震应力变化场；(d)取节面Ⅱ时14个月之后的应力场变化.黑色沙滩球为囊谦地震的震源机制解. Fig. 7 Coseismic and post seismic Coulomb failure stress changes on two nodal planes with the effective friction coefficient of 0.4 on the depth of 9km (Nodal planeI:strike=259.8°，dip=84.4°，rake=-21.1°；nodal planeⅡ:strike=352°，dip=69°，rake=-174°) (a)Coseismic stress change on nodal planeI；(b)Post seismic stress change on nodal planeI after 14 months；(c)Coseismic stress change on nodal planeⅡ；(d)Post seismic stress change on nodal planeⅡ after14 months.Black beach ball represents the focal solution of Nangqên earthquake . |
与选取玉树地震震源机制解作为接收断层参数相比，囊谦地震震源机制解两个节面的走向与玉树地震破裂面的走向均呈近45°夹角.由于库仑应力计算中选取的接收断层参数发生了改变，图 7 中玉树地震导致的库仑应力变化场虽然仍分布着四处库仑应力增强的扇区，但与图 4- 图 6中的库仑应力增强扇区相比，均逆时针旋转了一定角度.选取囊谦地震震源机制解不同节面计算得到的库仑应力变化场差距并不明显，选取节面I时，囊谦地震位于玉树地震破裂面西南方向的应力影区边缘地区，在囊谦地震的震中，玉树地震导致的同震库仑应力变化为-0.002 MPa，震后一年零两个月的黏弹松弛应力变化为-0.003MPa；选取节面Ⅱ时，囊谦地震位于西南影区的中间区域，同震库仑应力变化为-0.003 MPa，到囊谦地震发生前降低为-0.009 MPa.选取不同的节面计算得到的同震及震后黏弹松弛库仑应力变化均表明囊谦地震位于玉树地震所导致的应力影区.
3.2.2 不同黏滞系数对计算结果的影响囊谦地震发生在玉树地震震后一年零两个月，考虑同震及震后黏弹松弛效应对应力场变化的影响时，中下地壳以及上地幔的黏滞系数的不同取值可能会改变囊谦地震震前的应力场分布.Shan等[4]在参数的敏感性分析中，讨论了黏滞系数取值对计算结果的影响，发现中下地壳黏滞系数取值的不同对计算结果的影响要远大于上地幔黏滞系数的取值.因此，我们选取不同的中上地壳黏滞系数，计算了囊谦地震震前黏弹松弛应力场变化，计算中不考虑背景构造应力场，深度截面为9km，有效摩擦系数取中间值0.4，以囊谦地震震源机制解的节面Ⅱ作为接收断层参数.
由图 8可见，中上地壳黏滞系数取值较小时，中下地壳以及上地幔所积累的应力､应变释放较快，对上地壳应力状态的影响较大.当中下地壳黏滞系数ηc 取5×1017 Pa·s时，囊谦地震震中库仑应力变化为- 0.009 MPa，其中的同震应力变化为-0.003 MPa(图 7)，震后黏弹松弛效应所导致的应力变化为-0.006 MPa；而当ηc 取1×1020 Pa·s时，同震库仑应力变化为-0.003 MPa，震后黏弹松弛效应几乎没有影响.考虑震后黏弹松弛效应时，虽然黏滞系数的取值改变囊谦地震震中位置库仑应力变化的大小，但对震后库仑应力变化场的空间分布几乎没有影响，囊谦地震始终位于玉树地震破裂面西南反向的应力影区中央.这是由于对倾角为90°的纯走滑断层而言，随着黏弹松弛的过程，上地壳中同震库仑应力变化为正的区域应力进一步增强，而库仑应力变化为负的区域应力得到进一步释放，应力场的空间分布变化较小；而对于正断层和逆冲断层，震后上地壳中应力的重新分配一般是相反的，同震应力增加的区域和应力影区在松弛作用下相互转化，应力场的空间分布改变较大[23, 45-46].
|
图 8 有效摩擦系数取0.4时，接收断层参数取节面Ⅱ，取不同的地壳黏滞系数，14个月以后9km 深度上的库仑应力变化场 各图中右上角为相应的中下地壳黏滞系数取值.黑色沙滩球为囊谦地震的震源机制解. Fig. 8 Considering the uncertainty of viscosity，post seismic Coulomb failure stress changes on nodal planeⅡ with the effective friction coefficient of 0.4 on the depth of 9km after 14 months The value on the upper right side in each figure represents the corresponding viscosity in numerical calculation.Black beach ball represents the focal solution of Nangqên earthquake . |
根据地震应力触发理论，本文利用弹性位错理论和分层岩石圈模型计算了玉树地震引起的同震及震后黏弹松弛应力场变化，并研究MS7.1级玉树地震对余震分布的影响以及与2011年囊谦MS5.2级地震之间的触发关系，得到的主要结论有:
(1) 根据同震库仑应力的计算结果，玉树地震导致了四处明显的库仑应力增强的扇区，2010年4月13日至6月17日的870次ML>1.0级余震主要分布于主震破裂面附近区域以及破裂面东北端的应力增强扇区.考虑同震库仑应力变化对余震分布的影响时，有效摩擦系数的取值对计算结果的影响较小；不同深度截面的计算结果显示，5､10､15km 深度上的库仑应力变化场随深度变化较小，20km 深度上在破裂面附近变化较大，这是由于玉树地震的位错分布主要集中在15km 深度以上；考虑区域构造应力场的影响主要改变了玉树地震破裂面附近库仑应力变化场，占总数的86.7%余震位于库仑应力增强区，地震应力触发理论较好地解释了余震的分布.在研究主震对余震的触发作用时，应该考虑区域构造应力场的影响.
(2) 选取囊谦地震震源机制解的两个节面作为库仑应力计算中的接收断层参数，并且考虑不同黏滞系数下的玉树地震同震及震后黏弹松弛效应，模型计算结果均表明囊谦地震位于玉树地震所导致应力影区，仅依靠地震的静态､震后黏弹松弛应力触发理论，无法解释囊谦地震的发生，说明该次地震很可能是一次独立的事件，与该地区的背景地震活动性有关.
本研究用到的870 次ML>1.0 级余震集中分布在主震破裂面附近以及破裂面东北端的应力增强扇区内，其他三处应力增强扇区内几乎没有记录到余震的分布，导致该现象的原因可能是因为同震静态应力变化虽然增加这些地区的应力积累，提高了余震发生的概率，但是地震是否发生､震级大小还主要取决于断层所积累的构造应力状态和大小，例如:位于玉树地震破裂面东南端的巴塘-洛须段记录到的上一次大地震的发生在1896年，陈立春等[3]根据断层的平均滑移速率推测下一次大地震可能将发生在公元2339±67年，当前断层上的应力积累还处于地震循环的初级阶段；而玉树MS7.1级地震发生在当江-玉树段，该段上一次的大地震发生在1738年，陈立春等[3]推断下一次大地震可能发生在公元2012±30年，虽然玉树地震使得当江-玉树段东南段发生破裂，释放了该段上的应力积累，但位于地震破裂面西北方向的当江-隆宝段仍然维持着较高的应力积累.虽然玉树地震同震库仑应力变化同时增加了破裂面东南端和西北端的应力积累，但由于此前西北端已经积累了较强的构造应力，因此附近区域内余震活跃，今后同样具有较高的地震危险性.
根据模型的计算结果，玉树地震释放了囊谦地震震中位置上的应力积累，此外，在位于玉树地震破裂面西北方向的应力影区内也分布有一些余震，从地震应力触发的理论无法解释这些地震的发生.之所以有部分余震仍然发生在库仑应力下降的应力影区，可能的解释主要有三种:第一，在库仑应力计算过程中所产生的误差，例如:经过简化的断层破裂模型､简单的地壳分层模型没有考虑地壳的各向异性､余震发生的微断层实际的震源机制与先前假设最优化方位角之间的差异；第二，库仑应力下降只是降低了应力影区中地震发生的概率，并没有完全抑制地震的发生，地震发生与否是该区域内背景地震活动性和库仑应力变化共同作用的结果[47-48]；第三，位于应力影区的地震可能是由其他地震触发机制所触发的，例如动态应力触发，目前通常认为动态应力触发与地震波的传播有关，与静态应力触发相比，地震导致的动态应力场变化不存在应力影区[49].
考虑区域背景构造应力场的库仑应力场变化计算结果显示，沿玉树地震主破裂带附近为库仑应力增强的扇区，统计得到玉树地震震后两个月内870次ML>1.0级余震中有86.7%位于库仑应力增强区域，同震库仑应力变化场很好地解释了余震的分布.但实际上，离断层破裂面越近，断层破裂的空间分布和几何形态对库仑应力计算的影响越大，计算的结果越不可靠，并且在主震破裂面附近区域内余震的发震机理也越复杂.考虑区域背景构造应力场时，地震主破裂带附近区域内的余震均位于库仑应力增强区，如果全部被认为是由玉树地震同震库仑应力所触发的话，可能比实际被触发的地震数目偏高.
此外，本研究还存在一些不足，从地震触发角度而言，本文仅从静态弹性以及震后黏弹松弛应力触发的角度出发，计算了MS7.1级玉树地震所导致的应力场变化，并且讨论玉树地震对周边地震活动的影响.其他效应，如基于孔隙回弹､震后余滑等效应的与时间有关的应力触发和动态应力触发[23]等的综合考虑和研究将使本文结果更加全面.
致谢感谢中国地震局地球物理研究所“国家数字测震台网数据备份中心"为本研究提供地震波形数据.本文图件采用GMT[50]绘制.
| [1] | 倪四道, 王伟涛, 李丽. 2010年4月14日玉树地震:一个有前震的破坏性地震. 中国科学D辑 , 2010, 53(6): 791–793. Ni S D, Wang W T, Li L. The April 14th, 2010 Yushu earthquake, a devastating earthquake with foreshocks. Sci. China Earth Sci. (in Chinese) , 2010, 53(6): 791-793. DOI:10.1007/s11430-010-0083-2 |
| [2] | 张勇, 许力生, 陈运泰. 2010年青海玉树地震震源过程. 中国科学D辑 , 2010, 53(9): 1249–1251. Zhang Y, Xu L S, Chen Y T. Source process of the 2010 Yushu, Qinghai, earthquake. Sci. China Earth Sci. (in Chinese) , 2010, 53(9): 1249-1251. DOI:10.1007/s11430-010-4045-5 |
| [3] | 陈立春, 王虎, 冉勇康, 等. 玉树Ms7. 1级地震地表破裂与历史大地震. 科学通报 , 2010, 55(31): 1200–1205. Chen L C, Wang H, Ran Y K, et al. The Ms7. 1 Yushu earthquake surface ruptures and historical earthquakes. Chinese Sci. Bull. (in Chinese) , 2010, 55(31): 1200-1205. DOI:10.1007/s11434-010-0293-1 |
| [4] | Shan B, Xiong X, Zheng Y, et al. The co-seismic Coulomb stress change and expected seismicity rate caused by 14 April 2010 Ms=7. 1 Yushu, China, earthquake. Tectonophysics , 2011, 510(3-4): 345-353. DOI:10.1016/j.tecto.2011.08.003 |
| [5] | Stein R S. Earthquake conversations. Sci. Am. , 2003, 288(1): 72-79. DOI:10.1038/scientificamerican0103-72 |
| [6] | 单斌, 熊熊, 郑勇, 等. 2008年5月12日Mw7.9汶川地震导致的周边断层应力变化. 中国科学D辑 , 2009, 52(5): 593–601. Shan B, Xiong X, Zheng Y, et al. Stress changes on major faults caused by Mw7.9 Wenchuan earthquake, May 12, 2008. Sci. China Earth Sci. (in Chinese) , 2009, 52(5): 593-601. DOI:10.1007/s11430-009-0060-9 |
| [7] | King G C P, Stein R S, Lin J. Static stress changes and the triggering of earthquakes. Bull. Seismol. Soc. Am. , 1994, 84(3): 935-953. |
| [8] | Stein R S. The role of stress transfer in earthquake occurrence. Nature , 1999, 402(6762): 605-609. DOI:10.1038/45144 |
| [9] | Deng J S, Sykes L R. Evolution of the stress field in southern California and triggering of moderate-size earthquakes: a 200-year perspective. J. Geophys. Res. , 1997, 102(B5): 9859-9886. DOI:10.1029/96JB03897 |
| [10] | Deng J, Sykes L R. Stress evolution in southern California and triggering of moderate-, small-, and micro-size earthquakes. J. Geophys. Res. , 1997, 102(B11): 24411-24435. DOI:10.1029/97JB02127 |
| [11] | Hardebeck J L, Nazareth J J, Hauksson E. The static stress change triggering model; constraints from two southern California aftershock sequences. J. Geophys. Res. , 1998, 103(B10): 24427-24437. DOI:10.1029/98JB00573 |
| [12] | Wang W H, Chen C H. Static stress transferred by the 1999 Chi-Chi, Taiwan, earthquake: effects on the stability of the surrounding fault systems and aftershock triggering with a 3D fault-slip model. Bull. Seismol. Soc. Am. , 2001, 91(5): 1041-1052. |
| [13] | Marsan D, Bean C J. Seismicity response to stress perturbations, analysed for a world-wide catalogue. Geophys. J. Int. , 2003, 154(1): 179-195. DOI:10.1046/j.1365-246X.2003.01963.x |
| [14] | Jaeger J C, Cook N G. Fundamentals of Rock Mechanics. 3rd ed. London: Chapman & Hall, 1979. |
| [15] | Scholz C H. The Mechanics of Earthquakes and Faulting. New York: Cambridge University Press, 1990 : 439 . |
| [16] | Okada Y. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. , 1992, 82(2): 1018-1040. |
| [17] | Parsons T, Stein R S, Simpson R W, et al. Stress sensitivity of fault seismicity: a comparison between limited-offset oblique and major strike-slip faults. J. Geophys. Res. , 1999, 104(B9): 20183-20202. DOI:10.1029/1999JB900056 |
| [18] | Heidbach O, Ben-Avraham Z. Stress evolution and seismic hazard of the Dead Sea Fault System. Earth and Planetary Science Letters , 2007, 257(1-2): 299-312. DOI:10.1016/j.epsl.2007.02.042 |
| [19] | Xiong X, Shan B, Zheng Y, et al. Stress transfer and its implication for earthquake hazard on the Kunlun Fault, Tibet. Tectonophysics , 2010, 482(1-4): 216-225. DOI:10.1016/j.tecto.2009.07.020 |
| [20] | Wang R J. A simple orthonormalization method for stable and efficient computation of Green's functions. Bull. Seismol. Soc. Am. , 1999, 89(3): 733-741. |
| [21] | Wang R J, Lorenzo-Matín F, Roth F. PSGRN/PSCMP-a new code for calculating co-and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory. Comput. Geosci. , 2006, 32(4): 527-541. DOI:10.1016/j.cageo.2005.08.006 |
| [22] | Aki K, Richards P G. Quantitative Seismology. 2nd ed. California: University Science Books, 2002: 73. http://www.oalib.com/references/18990157 |
| [23] | 王有学, 钱辉. 青海东部地壳速度结构特征研究. 地学前缘(中国地质大学, 北京) , 2000, 7(4): 568–579. Wang Y X, Qian H. Study of crustal velocity structure in east Qinghai. Earth Science Frontiers (China University of Geosciences, Beijing) (in Chinese) , 2000, 7(4): 568-579. |
| [24] | Chen Y L, Niu F L, Liu R F, et al. Crustal structure beneath China from receiver function analysis. J. Geophys. Res. , 2010, 115: B03307. DOI:10.1029/2009JB006386 |
| [25] | 王未来, 吴建平, 房立华. 2010年玉树Ms7.1地震及其余震的双差定位研究. 中国科学D辑 , 2012, 42(7): 1037–1046. Wang W L, Wu J P, Fang L H. Relocating of the mainshock and its aftershock of 2010 Ms7.1 Yushu earthquake by double difference method. Sci. China Earth Sci. (in Chinese) , 2012, 42(7): 1037-1046. |
| [26] | Freed A M. Earthquake triggering by static, dynamic, and postseismic stress transfer. Annu. Rev. Earth Planet. Sci. , 2005, 33(1): 335-367. DOI:10.1146/annurev.earth.33.092203.122505 |
| [27] | Pollitz F F, Bürgmann R, Segall P. Joint estimation of afterslip rate and postseismic relaxation following the 1989 Loma Prieta earthquake. J. Geophys. Res. , 1998, 103(B11): 26975-26992. DOI:10.1029/98JB01554 |
| [28] | 邵志刚, 傅容珊, 薛霆虓, 等. 昆仑山Ms8. 1级地震震后变形场数值模拟与成因机理探讨. 地球物理学报 , 2008, 51(3): 805–816. Shao Z G, Fu R S, Xue T X, et al. The numerical simulation and discussion on mechanism of postseismic deformation after Kunlun Ms8.1 earthquake. Chinese J. Geophys. (in Chinese) , 2008, 51(3): 805-816. |
| [29] | Zhao L S, Helmberger D V. Source estimation from broadband regional seismograms. Bull. Seismol. Soc. Am. , 1994, 84(1): 91-104. |
| [30] | Zhu L P, Helmberger D V. Advancement in source estimation techniques using broadband regional seismograms. Bull. Seismol. Soc. Am. , 1996, 86(5): 1634-1641. |
| [31] | 韦生吉, 倪四道, 崇加军, 等. 2003年8月16日赤峰地震:一个可能发生在下地壳的地震?. 地球物理学报 , 2009, 52(1): 111–119. Wei S J, Ni S D, Chong J J, et al. The 16 August 2003 Chifeng earthquake: Is it a lower crust earthquake?. Chinese J. Geophys. (in Chinese) , 2009, 52(1): 111-119. |
| [32] | 郑勇, 马宏生, 吕坚, 等. 汶川地震强余震(Ms≥5. 6)的震源机制解及其与发震构造的关系. 中国科学D辑 , 2009, 52(6): 739–753. |
| [33] | Zheng Y, Ni S D, Xie Z, et al. Strong aftershocks in the northern segment of the Wenchuan earthquake rupture zone and their seismotectonic implications. Earth Planets Space , 2010, 62(11): 881-886. DOI:10.5047/eps.2009.06.001 |
| [34] | 罗艳, 倪四道, 曾祥方, 等. 汶川地震余震区东北端一个余震序列的地震学研究. 中国科学D辑 , 2010, 53(11): 1655–1664. Luo Y, Ni S D, Zeng X F, et al. A shallow aftershock sequence in the north-eastern end of the Wenchuan earthquake aftershock zone. Sci. China Earth Sci. (in Chinese) , 2010, 53(11): 1655-1664. DOI:10.1007/s11430-010-4026-8 |
| [35] | 郑秀芬, 欧阳飚, 张东宁, 等. "国家数字测震台网数据备份中心"技术系统建设及其对汶川大地震研究的数据支撑. 地球物理学报 , 2009, 52(5): 1412–1417. Zheng X F, Ouyang B, Zhang D N, et al. Technical system construction of Data Backup Centre for China Seismograph Network and the data support to researches on the Wenchuan earthquake. Chinese J. Geophys. (in Chinese) , 2009, 52(5): 1412-1417. |
| [36] | Zhu L P, Rivera L A. A note on the dynamic and static displacements from a point source in multilayered media. Geophys. J. Int. , 2002, 148(3): 619-627. DOI:10.1046/j.1365-246X.2002.01610.x |
| [37] | 房立华, 吴建平, 张天中, 等. 2011年云南盈江Ms5. 8地震及其余震序列重定位. 地震学报 , 2011, 33(2): 262–267. Fang L H, Wu J P, Zhang T Z, et al. Relocation of mainshock and aftershocks of the 2011 Yingjiang Ms5.8 earthquake in Yunnan. Acta Seismologica Sinica (in Chinese) , 2011, 33(2): 262-267. |
| [38] | 姜枚, 王有学, 钱辉, 等. 造山的高原——青藏高原及其邻区的宽频地震探测与地壳上地幔结构. 北京: 地质出版社, 2009 : 69 . Jiang M, Wang Y X, Qian H, et al. Orogenic Plateau-Tibetan Plateau and Adjacent Areas of Broad Band Seismic Detection and Structure of the Crust and Upper Mantle (in Chinese). Beijing: Geological Publishing House, 2009 : 69 . |
| [39] | Dieterich J H. A constitutive law for rate of earthquake production and its application to earthquake clustering. J. Geophys. Res. , 1994, 99(B2): 2061-2618. |
| [40] | Catalli F, Cocco M, Console R, et al. Modeling seismicity rate changes during the 1997 Umbria-Marche sequence (central Italy) through a rate-and state-dependent model. J. Geophys. Res. , 2008, 113: B11301. DOI:10.1029/2007JB005356 |
| [41] | Xu C J, Wang J J, Li Z H, et al. Applying the Coulomb failure function with an optimally oriented plane to the 2008 Mw7.9 Wenchuan earthquake triggering. Tectonophysics , 2010, 491(1-4): 119-126. DOI:10.1016/j.tecto.2009.09.019 |
| [42] | Hainzl S, Zoller G, Wang R. Impact of the receiver fault distribution on aftershock activity. J. Geophys. Res. , 2010, 115: B05315. DOI:10.1029/2008JB006224 |
| [43] | 王辉, 张国民, 汪素云, 等. 应用地震学方法研究中国大陆活动地块应力应变场. 地球物理学报 , 2004, 47(6): 1035–1043. Wang H, Zhang G M, Wang S Y, et al. Stress and strain fields of active tectonic blocks in the China mainland deduced by seismological methods. Chinese J. Geophys. (in Chinese) , 2004, 47(6): 1035-1043. |
| [44] | Ma K F, Chan C H, Stein R S. Response of seismicity to Coulomb stress triggers and shadows of the 1999 Mw=7.6 Chi-Chi, Taiwan, earthquake. J. Geophys. Res. , 2005, 110: B05S19. DOI:10.1029/2004JB003389 |
| [45] | Freed A M, Lin J. Time-dependent changes in failure stress following thrust earthquakes. J. Geophys. Res. , 1998, 103(B10): 24393-24409. DOI:10.1029/98JB01764 |
| [46] | Nostro C, Piersanti A, Cocco M. Normal fault interaction caused by coseismic and postseismic stress changes. J. Geophys. Res. , 2001, 106(B9): 19391-19410. DOI:10.1029/2001JB000426 |
| [47] | Felzer K R, Brodsky E E. The absence of stress shadows. Seismol. Res. Lett. , 2003, 75: 285. |
| [48] | Mallman E P, Zoback M D. Testing Coulomb stress-transfer models with seismicity rates for the Landers and Kobe earthquakes. Seismol. Res. Lett. , 2003, 75: 285. |
| [49] | Harris R A. Introduction to special section: Stress triggers, stress shadows, and implications for seismic hazard. J. Geophys. Res. , 1998, 103(B10): 24347-24358. DOI:10.1029/98JB01576 |
| [50] | Wessel P, Smith W H F. New, improved version of Generic Mapping Tools released. Eos Trans AGU , 1998, 79(47): 579. DOI:10.1029/98EO00426 |