地球物理学报  2015, Vol. 58 Issue (6): 1919-1930   PDF    
基于经验模型和物理模型研究2013 MS7.0芦山地震余震序列
米琦1, 申文豪1,2, 史保平1    
1. 中国科学院大学地球科学学院, 北京 100049;
2. 中国地震局地壳应力研究所, 北京 100085
摘要:基于中国地震台网中心2013 MS7.0芦山地震余震数据我们首先确定了余震空间分布范围并根据G-R关系计算了主震后半小时内的完备震级Mc=3.5, 并且得到了ML≥3.5和ML≥3.0的地震在2001年至芦山地震前的背景场地震发生率.通过Omori-Ustu经验定律和两种Dieterich模型对芦山地震余震发生率的拟合, 我们发现阶梯型Dieterich模型只能模拟p=1的情况, 从而造成了模拟曲线与观测数据的差别;前人研究表明震后滑移同样是产生余震的原因, 如果假设余震序列由主震静态剪应力Δτ和震后滑移共同作用所产生, 我们数值模拟得到的对数型Dieterich模型能够较好地推断余震发生率R随时间t增加而衰减的趋势, 能够从物理机制上解释MS7.0芦山地震余震序列衰减指数大于1这一现象.通过对数型Dieterich模型的拟合并结合Andrews的方法, 我们还得到MS7.0芦山地震约为0.155 MPa, ta约为8.4年, 这一值与前人研究结果十分接近.
关键词余震衰减     芦山地震     Dieterich模型     震后滑移    
Aftershock decay of the 2013 Lushan MS7.0 earthquake derived from the empirical and physical models
MI Qi1, SHEN Wen-Hao1,2, SHI Bao-Ping1    
1. College of Earth Science, University of Chinese Academy of Sciences, Beijing 100049, China;
2. Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China
Abstract: Accurately predicting aftershocks decay after the main shock over time is conducive to learn the aftershocks triggering physical mechanism. On the other hand, the predicted duration of aftershocks is also an important aspect of seismic hazard analysis. In this study, by fitting the seismicity decay with both the empirical and physical models, we try to study Lushan earthquake triggering mechanism, and figure out how coseismic stress, afterslip and other parameters affect aftershocks decay.
    First, based on the seismicity data from the China Earthquake Network Center, the magnitude of completeness Mc for the early part was estimated using a maximum likelihood procedure. We then compute the observed seismicity decay and fit it with both the empirical and physical models. By using a grid search method, parameters of predication models could be determinated.
    The magnitude of completeness Mc for the early part (~0.02 day) was estimated at 3.5 using a maximum likelihood procedure. We found that, with the stress step formed Dieterich model, a p value larger than 1 could not be explained, which was in contradiction with empirical observations of p>1 aftershock decays. We have shown that by modeling a coseismic step and a postseismic stress change with a constant background stress rate, the Omori exponent larger than one for the MS7.0 Lushan Earthquake aftershock sequence could be well explained and fitted. With Andrews' method, we calculated the distributed stress changes on the fault plane, and we inferred ~0.155 MPa, also, the logarithmic stressing Dieterich model predicts a more reasonable aftershocks duration ta~8.4 years, which is close to former researchers' result.
    By showing the merits and limitations of the modified Omori law and the step function formed Dieterich model, we find that the stress step formed Dieterich model cannot explain an Omori law decay with p>1, thus, this physical model may overestimate the aftershock duration. We show that modeling by a coseismic step, a postseismic stress change with a constant background stress rate, the Omori exponent larger than one for MS7.0 Lushan earthquake could be explained.
Key words: Aftershock decay     Lushan earthquake     Dieterich model     Afterslip    
1 引言

余震是指主震发生后,围绕主震断层一定尺度区域范围内发生的一系列地震,主要表现为在该区域范围内地震活动性的显著增加,且随着时间增加而衰减.余震震级通常比主震要小,其影响范围在空间上超过了主震断层的破裂尺度.余震序列的持续时间通常为数年乃至数十年不等,在此期间区域地震活动性明显高于背景场水平,并随时间逐渐衰减.Omori(1894)最先发现余震发生率R正比于1/t,其中t为时间,自此这一用于拟合余震发生率随时间衰减的经验模型及Utsu(1961)提出的改进模型(Modified Omori Law,MOL)

已经被成功用于全球范围内多次余震序列的拟合.地震学首先是一门观测学科,在观测现象未获得明确物理解释之前,像Omori-Ustu定律这种地震统计学观测定律,或者说经验模型对于认知地震活动性、了解地震活动规律发挥着举足轻重的作用(Ogata,1988).在公式(1)中,R(t)表示余震发生率,即单位时间里余震发生个数;公式右端t为主震后时间,K,cp为待确定的系数.在这三个参数中,p值的大小反映了余震发生率随时间衰减的快慢,称之为衰减指数,其取值一般在0.8~1.5之间,并且该参数与起算震级的选取没有必然联系(Utsu et al.,1995).也就是说,对于同一个主震引发的余震序列,不同起算震级得到的p值应该是相同的.通过分析美国加利福尼亚地区地震数据,Ouillon和Sornette(2005)发现p值随主震震级M的增加而增加.近期Hainzl和Marson(2008)利用非均匀的断层面应力变化模型和Dieterich 速率和状态模型给出了对于上述现象可信的物理解释.K值一般被视为余震序列产出量的大小(productivity),与主震震级M存在正相关关系:K~10αM,其中α为大于零的常数(Hainzl and Marsan,2008).在三个参数中,c值一直备受争议,通常认为c值表示主震后高水平余震发生率持续的时间(饱和期),其值在0.01~1.0 day之间.Enescu等(2007)Narteau等(2002)的研究认为主震后初始阶段出现的高水平余震发生率持续期(即c值)是余震序列衰减的真实特征,而Kagan(2004)则认为该现象并非是真实的衰减特征,而是由于主震后较短时间内地震目录(尤其是中小地震目录)缺失所导致的饱和现象.

随着震源物理研究的深入,余震触发的物理学机制在近些年来得到广泛关注并取得一系列新的成果,许多物理模型开始被用于解释余震发生率衰减 现象.由于在统计地震学中的成功应用,Omori-Ustu定律已经成为这些物理模型的“试金石”,因此几乎所有物理模型都试图得到同Omori-Ustu定律相似样式的衰减形式,并赋予K,c,p以明确的物理内涵.在这些物理模型中,本文采用的是速率和状态依赖性摩擦本构关系并结合一维弹簧-滑块模型所得到的理论解,通常称为Dieterich 模型,选用该模型的原 因是近些年来该模型在地震学中广泛而成功的应用.

汶川地震发生后,中国地震局及四川省地震局投入了大量的人力物力用于龙门山断裂带的地震活动性监测,相较于汶川地震之前该区域的监测能力有了明显提高.2013年4月20日芦山地震发生后,得益于监测能力的改善,中国地震台网中心为我们提供了一套相对完备的余震数据.在本文中,借助这一余震数 据体,我们将系统分析该余震序列的时空分布特征.首先,确定余震区域的范围和主震后半小时(~0.02 day)内截止震级的大小,以确保数据体的完备和有效.其次,通过Omori-Ustu定律和两种Dieterich模型对MS7.0芦山地震余震序列衰减的拟合,我们将分别讨论经验模型和阶梯型Dieterich模型存在的优缺点,进而发现结合同震应力变化和震后剪应力变化率衰减的Dieterich模型可以很好地解释和模拟MS7.0芦山地震余震序列衰减p值大于1这一现象.另一方面,我们还利用这一模型讨论了芦山地震的余震持续时间及Dieterich摩擦定律中参数的确定.

2 2013 MS7.0 芦山地震余震序列时空分布特征

准确地描述余震序列的时空特征对于我们正确认知断层内部演化状态,了解构造环境与地震活动性的关系都起着重要作用.2013 MS7.0芦山地震的发震断层一度被认为是龙门山断裂带南端的大川—双石断裂,而李传友等(2013)徐锡伟等(2013)的研究结果则认为芦山地震的发震断层为一条现今尚未出露地表、其上断点仍埋藏在地下9 km以下地壳中的一条盲逆断层,该发震构造是芦山之下、大川—双石断裂和新开店断裂之间的龙门山前缘滑脱带,此滑脱带在该段的运动导致了这次地震的发生,并可能带动了它上面的大川—双石和新开店等断裂的活动.选取29.8°N—30.5°N,102.5°E—103.3°E作为研究区域,该区域从2001年到汶川地震发生之前7年间ML≥3.5 地震有4次,而在汶川地震发生后到芦山地震发生之前5年间ML≥3.5 地震有3次,基本维持在背景场水平,也就是说汶川地震的发生并未造成该区域内地震活动性的增加,芦山地震后该区域内的地震数目激增是由芦山主震触发的.图 1给出了芦山地震后到2014年5月25日ML≥3.5 余震分布,可以看到余震呈密集条带状分布,其东北至西南方向走向展布约40 km,西北至东南向展布约25 km.图 1右下角的子图给出的是MS7.0芦山地震和MS8.0汶川地震的余震(ML≥4.0,2008-05-12—2014-05-25)分布图,可以看到两次地震的余震密集区相距约50 km(杜方等,2013).

图 1MS7.0芦山地震震中位置及其余震序列分布. 统计时间为2013年4月20日到2014年5月25日, 震级ML≥3.5
蓝色五角星为主震震中位置,红色圆点为余震位置.
Fig. 1Map of the location of MS 7.0 Lushan earthquake and its aftershocks distribution from the mainshock to May 25, 2014 for ML≥3.5
Blue star indicated epicenter location of the mainshock, and red dots indicated aftershock locations.

此外,由于在主震发生后的数小时内中低震级余震的数目缺失,我们必须求出在主震后短时间内最小完备震级(minimum magnitude of complete recording),Mc.过高的Mc有可能导致一些可用数据的丢失,而如果低估了Mc则会得到错误的地震活动性参数值,从而导致分析中出现偏差(Mignan and Woessner,2012).Mc通常可由古登堡-里克特(G-R)定律拟合观测到的震级-频度分布(FMD)来估计,并且取两者相符部分的起始震级作为Mc的估值(Zúniga and Wyss,1995).表 1给出了主震后半小时内(~0.02 day)余震序列目录(中国地震台网中心,2014).

表 1 MS 7.0芦山地震后半小时内(~0.02 day)主震及余震序列目录 Table 1 Catalogue of MS 7.0 Lushan earthquake and aftershock sequence in first half hour

利用表 1的余震序列数据,我们估算了在主震后半小时内(~0.02 day)余震序列的完备震级Mc 的取值为3.5,并得到了G-R关系中a和b的取值(图 2b).同样,根据2001年至芦山地震前的背景场地震数据我们得到了这段时间的完备震级为1.4(图 2a),明显小于主震后半小时内余震序列的完备震级Mc=3.5,说明在背景场地震数据中ML≥3.5的地震目录是完整的,根据G-R关系我们估算得到ML≥3.5的地震在2001年至芦山地震前的背景场发生率为0.0016/day.

图 2(a) 2001年至主震发生前研究区域震级-频度分布图; (b) MS7.0芦山地震序列主震后半小时内震级-频度分布图 Fig. 2(a) Magnitude-frequency distribution in the study area from 2001 to the mainshock; (b) Magnitude -frequency distribution of aftershocks for MS 7.0 Lushan earthquake in the first half hour after the main shock
3 利用经验模型拟合芦山地震余震序列

选取研究区域内2013年4月20日主震发生后至2014年5月25日ML≥3.5的余震数据,并采用Narteau等(2002)给出的计算方法,得到了该段时间余震发生率随时间的变化(图 3).从图 3中可以看到实测余震发生率随时间快速衰减,但是在震后约1.5天和震后20天左右出现了两次异常突增,前一次地震活动性的突增主要是由于在2013年4月21日17时5分左右发生的一次ML=5.4的余震引起的,此次地震也是芦山地震迄今为止最大的余震;第二次则是由于在2013年5月11日一天内发生了三次4.0级以上的余震,引起了该时间段内地震活动性的突增.通过最小二乘法的拟合,我们得到公式(1)与观测数据拟合最优时K,c,p的取值分别为28.12,1.21和0.05(图 3).从图 3中可以看到利用Omori-Ustu经验公式可以很好地拟合实测数据,但是其参数c的取值在学术界存在较大争议.Vidale等(2003)发现通过对在主震后数分钟的波形数据进行滤波得到的余震事件数目要数倍于地震目录记录到的余震事件,因此c值是震后地震目录缺失产生的.Kagan(2004)通过计算地震矩释放率发现c的取值接近0,并且应该是一个负值.这些发现都更新了研究者对Omori-Ustu定律的认识,而正如前文中提到的那样,Omori-Ustu定律最大的问题在于其虽然对物理现象描述很好但是终究缺乏明确的物理机制解释.

图 3MS7.0芦山地震主震后ML≥3.5余震发生率随时间变化示意图 Fig. 3Observations and modeling of seismicity rate during Lushan MS7.0 earthquake aftershocks sequence for earthquake ML≥3.5
Circles indicated observed data, the green solid line indicated the fitting curve of Omori-Ustu law, the red and purple solid line indicated fitting curve of Eqs. (5) and (7) respectively, and the blue solid line indicated the curve fitting by the logarithmic stressing Dieterich model.
4 利用物理模型拟合芦山地震余震序列

近些年来随着地震科学的发展,越来越多的研究者致力于从震源物理机制上解释主震后余震活动性的时空衰减规律,这些物理机制包括震后蠕变,流体扩散,速率和状态依赖性摩擦定律,应力腐蚀,损伤力学和亚临界裂纹扩展机制等.在本研究中针对MS7.0芦山地震我们使用的是速率和状态依赖性摩擦定律即Dieterich模型,该模型为我们提供了一个对断层属性进行复杂定量试验观测的基本框架(刘博研等,2013).如果给定不同的选取参数,Dieterich既能模拟断层摩擦滑动的稳定状态又能模拟非稳定状态.在失稳状态下,Dieterich模型能够提供不同应力加载形式与地震活动性之间的对应关系,而这种关系可以用于预测不同应力变化所触发的余震发生率.这些应力变化的来源是多方面的,包括主震造成的静态应力或动态应力,震后滑移,岩浆入侵或火山喷发造成的瞬态形变,慢地震或者潮汐.

Dieterich(1994)引入一个应力演化参数γ用来 描述断层上加载的剪应力随时间变化的过程(假设 正应力不变),

其中A是断层本构参数,为实验室常数,σ为正应力,为剪应力变化率.根据公式(2)再结合断层失稳时间的表达式就可以得出该模型下主震后余震发生率R的解析表达式(Dieterich,1994):

式中r为背景场地震发生率,r为背景场剪应力变化率.在公式(3)等号右边,对一个特定研究区域来说rr均为常数,可见余震发生率R的变化取决于应力演化参数.同样在公式(2)中,一般取为常数,所以γ的变化取决于加载的剪应力变化率γ的初值γ0.因此,不同的剪应力变化率加载形式和不同的应力演化参数初值γ0决定着余震发生率R的变化.在本文中我们将以两种Dieterich模型来拟合MS7.0芦山地震余震序列的观测数据,并讨论两种模型的差别.

4.1 阶梯型Dietrich模型拟合

假设余震序列仅由主震造成的静态剪应力Δτ产生,则主震后的应力演化参数γ初值可以写为:

另外再假设剪应力变化率在主震前后不变,由公式(2),(3)和(4)可得该情况下余震发生率R的解析表达式,

其中,ta=Aσ/r,表示余震发生率R衰减到背景场水平r的特征时间尺度,即余震持续时间.由于在这两种假设下剪应力随时间变化为阶梯型函数,因此称该模型为阶梯型Dieterich模型.由公式(5),在t远小于ta时,可以得到如下近似:

将公式(6)重新整理可以得到:

观察公式(7),我们发现等号右端具有Omori-Ustu定律相同的样式,对比公式(1)可以得到p=1,c=ta/[exp(Δτ/)-1],K=tar/[1-exp(-Δτ/)].因此阶梯型Dieterich模型为Omori-Ustu定律提供了基于速率和状态依赖性摩擦定律下的物理解释,并为其参数赋予物理意义.值得注意的是,在这种情况下c值是明确存在的,并且不为0或者负值,也就是说在阶梯型Dieterich模型下,余震发生率在开始衰减之前会出现一定时间的高水平余震发生率持续期(饱和期),其值决定于Δτr.对于一个特定研究区域来说r均为常数,因此c值的大小取决于Δτ的大小.相同的分析同样适用于K值,我们会发现K值的大小同样也与Δτ密切相关.

针对MS7.0芦山地震余震序列我们分别用公式(5)和公式(7)来拟合,在这两个公式中需要调整两个参数值,一个是ta,另一个是Δτ/的值.通过网格搜索法,我们得到拟合优度最高时ta=2.3×104 day,约63年;Δτ/=13.5(图 3).从图 3中可以看到在余震序列的前一部分公式(5)和公式(7)与观测数据符合较好,而在后一部分明显与观测数据产生分离.由Omori-Ustu定律的拟合我们知道MS7.0芦山地震余震发生率的衰减指数p=1.21,而公式(5)和公式(7)却只能模拟p=1的情况,因此造成了模拟曲线与观测数据的分离,事实上这一点正是阶梯型Dieterich模型在应用中所存在的最大问题.大量的地震观测数据表明衰减指数p值介于0.8~1.5(Utsu et al.,1995),阶梯型Dieterich模型显然不能解释p≠1的情况.为解决这一问题 Hainzl和Marson(2008)Helmstetter和Shaw(2006)等研究者引入了运动学有限断层模型(如k-2模型),将主震后的静态应力变化视为不均匀的(有正有负)而非单一值,这样就可以得到p<1的阶梯型Dieterich模型,但是这种模型依然不能解释p>1的情况,因此在该模型下得到的余震持续时间 往往会被高估,进而导致余震地震危险性分析的误差.

4.2 对数型Dietrich模型拟合

除了余震现象之外,在主震后还会发生明显的震后滑移,而大量研究则表明余震的衰减与震后位移的衰减遵从相同的时间衰减规律(Perfettini and Avouac,2007),这一现象使得许多研究者开始关注震后滑移与余震之间的关系,并且越来越多的证据 表明震后滑移能够触发余震(Perfettini and Avouac,2007; Hsu et al.,2006; Bourouis and Bernard,2007).Dieterich(1994)Helmstetter和Shaw(2009)以及Hainzl和Marsan(2008)等研究给出了震后滑移引起的剪应力变化率解析形式:

式中0为常数,t*为特征时间,q为衰减指数.公式(8)既能表示剪应力加载(>0)也能表示剪应力卸载(<0).在本研究中我们假设q=1,对公式(8)积分即可得到剪应力随时间变化的关系:

假设余震序列由主震静态剪应力Δτ和震后滑移共同作用产生,即主震后剪应力变化率按公式(8)的形式变化,那么将公式(8)代入公式(2)得到在该假设下对应的应力演化参数γ为:

其中m=0t*/(Aσ)为待确定的参数.考虑到主震静态剪应力Δτ的作用,公式(10)中的γ0与公式(4)相同.将公式(10)和公式(4)代入公式(3)我们就得到了另外一种形式的Dieterich余震发生率衰减模型,称之为对数型Dieterich模型:

在对数型Dieterich模型的应用中,我们在公式(8)中还增加了背景场剪应力变化率的影响,即将公式(8)写为:

,并利用Runge- Kutta法数值求解余震发生率R的结果.图 4给出了利用数值计算得到的余震发生率随时间变化示意图(=0.1 MPa,Δτ=0.92 MPa,t*=1.0 day,r=10-6 MPa/day),很显然m的取值对衰减曲线有重要影响:当m为正值时,0>0,表示应力加载,在这种情况下主震后的余震发生率R先增加到一个峰值然后再衰减,并且衰减指数p=1;当m=0时,=0,表示仅有主震静态应力和r作用而没有震后滑移的影响,这种情况等同于阶梯型Dieterich模型,因此可以说阶梯型Dieterich模型实际上是对数型Dieterich模型的一个特例;当m为负值时,0<0,表示应力卸载,余震发生率R会以大于1的衰减指数衰减,并且m越小衰减越快.

图 4基于对数型Dieterich模型得到的地震发生率示意图.=0.1 MPa,Δτ=0.92 MPa,t*=1.0 day,r=10-6 MPa/day,m分别取-1.5,-0.5,0,0.5和1.5.m= 0时表示仅由同震静态应力诱发余震 Fig. 4Seismicity rate triggered by a continuous stress change given by =0.1 MPa,Δτ=0.92 MPa,t*=1.0 day,r =10-6 MPa/day and for different amplitudes of the postseismic stress change, from m= -1.5 (bottom) to m=1.5 (top). The case m=0 represents seismicity rate triggered only by the coseismic stress step

为了用对数型Dieterich模型拟合MS7.0芦山地震余震序列,有五个参数需要确定:Δτmt*r,显然如果同时调整五个参数来拟合实际观测数据会有很大难度,因此我们先确定其中几个参数的值.利用快速傅里叶变换方法(Ripperger and Mai,2004),我们根据张勇等(2013)王卫民等(2013)给出的MS7.0芦山地震滑移模型计算了主震后断层面上的静态应力变化,如图 5所示.Helmstetter和Shaw(2006)指出主震后短时间内余震发生率的激增主要是受最大应力变化影响而非平均应力变化,这里所说的应力变化均为应力增加,即应力变化为正,因为只有在应力变化为正的区域才有可能产生余震.无论是张勇模型还是王卫民模型我们得到最大应力变化Δτ≈2.0 MPa.另外,背景 场剪应力变化率r可以近似等于Δτs/tr(Dieterich,1994),Δτs为主震应力降,本文取4.5 MPa(Shen et al.,2014);tr为地震复发周期,可以通过用地震同震位移除以断裂滑动速率近似求取得到:地震地质结果表明龙门山断裂带中单条断裂的滑动速率为~1.0 mm·a-1(张培震等,2008),而张勇模型和王卫 民模型给出的最大同震位移分别为1.30 m和1.59 m,因此我们得到MS7.0芦山地震复发周期为1300~1590年.由此r可以求得为7.75×10-6 ~9.48×10-6 MPa/day,本文取其平均值8.61×10-6 MPa/day作为输入参数.最后通过网格搜索法,得到该模型与观测数据拟合优度最高时=0.155 MPa,m=-1.0,t*=5.0 day(图 3).

图 5MS7.0芦山地震主震后断层面上的滑动量及静态应力变化
(a) 张勇等(2013)给出的滑移模型;(b) 王卫民(2013)等给出的滑移模型;(c) 基于张勇滑移模型得到的断层应力变化示意图;(d) 基于王卫民滑移模型得到的断层应力变化示意图.
Fig. 5Slip models and stress changes on the fault of MS7.0 Lushan earthquake
(a) Slip model of Zhang et al. (2013); (b) Slip model of Wang et al. (2013); (c) Stress changes calculated from Zhang model; (d) Stress changes calculated from Wang model.

除了ML≥3.5的地震之外,我们还给出了ML≥3.0 余震序列实测数据与经验模型及物理模型的对比(图 6),根据该区域G-R关系我们计算得到ML≥3.0余震序列背景场地震发生率为0.0037/day.在使用阶梯型Dieterich模型和对数型Dieterich模型进行拟合时,所用到的参数除背景场地震发生率r外其他参数与ML≥3.5余震序列拟合中采用参数完全一致.另外,用Omori-Ustu经验定律对ML≥3.0余震序列拟合时K值调整为58.12.对比图 3图 6,可以看到两种截止震级得到的余震衰减指数一致,均为p=1.21,并未随截止震级改变而改变,不同的是ML≥3.0余震序列在主震后初期由于数据不完整,发生率明显低于模型预测结果.从图 3图 6中可以看到对数型Dieterich模型预测的余震发生率衰减趋势与观测数据符合较好,并且 衰减指数与经验模型一致,因此我们可以推测MS7.0 芦山地震余震序列是由主震静态应力和震后滑移共同作用产生的.

图 6MS7.0芦山地震主震后ML≥3.0余震发生率随时间变化示意图 Fig. 6Observations and modeling of seismicity rate during Lushan Ms7.0 earthquake aftershocks sequence for earthquake ML≥3.0
Circles indicated observed data, the green solid line indicated the fitting curve of Omori-Ustu law, the red and purple solid line indicated fitting curve of Eqs. (5) and (7) respectively, and the blue solid line indicated the curve fitting by the logarithmic stressing Dieterich model.
5 结果与讨论

在对数型Dieterich模型中的5个物理参数包括背景场剪应力变化率r、特征时间t*、系数m、摩擦阻力参数和静态应力变化Δτ,将这些参数代入数值计算后即会得到余震发生率模拟曲线,而描述模拟曲线衰减特征主要有三个指标:初始时刻地震发生率大小、高地震发生率饱和期时间长短、余震发生率随时间衰减快慢,这三个指标即对应着Omori定律中的参数K,c,p.因此,我们可以通过调试这5个参数来研究它们对模拟曲线衰减特征的影响.首先给定一组标准取值表 2.

表 2 对数型Dieterich模型模拟余震衰减使用的标准参数 Table 2 Standard parameters used in the afterslip triggering Dieterich Model

为分离各个参数对模拟曲线的影响,我们设定在每个方案中只变化一个参数的量,其他参数采用表 2给出的标准值.每个方案所调整的参数名和值如表 3所示,每一方案的最终结果由图 7表示.

表 3 参数划分方案 Table 3 Parameters scenarios

图 7的子图中,不同的颜色分别代表同一种参数的不同取值模拟生成的余震发生率衰减曲线.通过图 7我们发现每一种参数对模拟曲线衰减特征的影响是不同的:背景场剪应力变化率r既影响衰减指数p又影响c的长短,并且rc成反比而与p成正比,对初始地震发生率K则没有影响;特征时间t*只影响c;m只影响衰减指数p;而摩擦阻力参数和静态应力变化Δτ只对K有影响.事实上,Dieterich(1994)指出对于震后滑移引起的余震模型,即公式(11)中右边括号里第一项在前期起到主导作用,于是R/r ≈exp(Δτ/)(1+t/t*)m,很显然从该近似表达式中可以看到t*m分别控制着c值大小和衰减指数p,而和Δτ只对K有影响.

图 7不同参数对对数型Dieterich模型余震衰减特征的影响 Fig. 7Comparison of the decay for afterslip triggering Dieterich Model based on different parameters

无论是阶梯型Dieterich模型还是对数型Dieterich模型,又或是其他基于速率和状态依赖性摩擦定律的物理模型,断层摩擦本构参数A和正应力σ的取值都起着非常关键的作用.现有研究表明,A值主要与温度和断层岩性有关,而正应力σ则与深度、区域应力场、断层方位和孔隙压力等因素有关系.在实际应用中一般被视为常数,其中A的参考值0.005~0.012是在实验室得到的,而在天然地震应用中如何取值一直没有明确限制.在许多基于速率和状态依赖性摩擦定律的研究中,的估算值为0.01~0.1 MPa(Hainzl et al.,2010),Perfettini和Avouac(20042007)针对1992年L and ers MW7.3地震和1999年Chi-Chi MW7.6地震分别得到约为0.47~0.53 MPa和0.34~1.5 MPa.Cochran等(2004)在对潮汐诱发地震的研究中给出的区间为0.048~0.11 MPa.在本文中针对MS7.0芦山地 震,基于阶梯型Dieterich模型,得到Δτ/()=13.5,根据上述已经求出Δτ约为2.0 MPa,所以=0.148 MPa,基于对数型Dieterich模型,求得=0.155 MPa.如果取A约为0.005~0.012,可以得到断层面上的正应力σ约为12.3~31.0 MPa.

Dieterich(1994)发现浅源地震的余震持续时间ta在10.2年左右,Parsons(2002)的研究发现全球范围内主震诱发的余震持续时间ta约为7~11年.在图 3图 6中,对数型Dieterich模型和Omori-Ustu定律预测的ta值为8.4年,这一值与上述结果十分接近,而阶梯型Dieterich模型的结果明显过大.另外从图中还可以看到,在衰减到背景场地震发生率r后,对数型Dieterich模型预测R会继续衰减一定时间,然后又逐渐增加恢复到背景场水平,也就是说在主震发生8.4年后研究区域有可能出现地震发生率低于背景场水平的现象,形成地震空区,而阶梯型Dieterich模型则不能有效解释地震空区现象.

尽管Dieterich 模型在实际应用中取得巨大的成功,但是由于其本身来源于室内的岩石物理学实验,在该模型中一些重要本构参数的取值(如A)是否能直接应用于天然地震仍存在较大争议,另外这些本构参数是否为常数又或者与构造环境有关系也存在不确定性.在运用该模型进行模拟时我们还忽略了流体压力、动态应力变化等因素的影响,这也有可能给结果带来不确定性.

6 结论

通过经验性的Omori-Ustu定律和两种Dieterich 模型对芦山地震余震发生率的拟合我们发现阶梯型Dieterich模型只能模拟p=1的情况,因此造成了模拟曲线与观测数据的分离;而如果假设余震序列由主震静态剪应力Δτ和震后滑移共同作用所产生,这样得到的对数型Dieterich模型能够较好预测芦山余震发生率衰减趋势,与观测数据符合较好,能够从物理机制上解释芦山地震余震序列衰减指数大于1这一现象.除此之外对数型Dieterich模型还可以合理解释地震空区现象,有效地弥补了阶梯型Dieterich模型的缺陷.综上所述,通过经验模型和物理模型的对比,我们有理由相信主震静态应力和震后滑移共同作用产生了芦山地震余震序列,我们得到的模型参数结果对于研究芦山地区地震危险性 以及分析该区域未来地震活动性具有重要参考意义.

参考文献
[1] Bourouis S, Bernard P. 2007. Evidence for coupled seismic and aseismic fault slip during water injection in the geothermal site of Soultz (France), and implications for seismogenic transients.Geophysical Journal International, 169(2): 723-732, doi: 10.1111/j.1365-246X.2006.03325.x.
[2] Cochran E S, Vidale J E, Tanaka S. 2004. Earth tides can trigger shallow thrust fault earthquakes. Science, 306(5699): 1164-1166.
[3] Dieterich J H. 1994. A constitutive law for rate of earthquake production and its application to earthquake clustering. Journal of Geophysical Research, 99(B2): 2601-2618.
[4] Du F, Long F, Ruan X, et al. 2013. The M7.0 Lushan earthquake and the relationship with the M8.0 Wenchuan earthquake in Sichuan, China. Chinese Journal of Geophysics (in Chinese), 2013, 56(5): 1772-1783, doi: 10.6038/cjg20130535.
[5] Enescu B, Mori J, Miyazawa M. 2007. Quantifying early aftershock activity of the 2004 mid-Niigata Prefecture earthquake (Mw6.6). Journal of Geophysical Research: Solid Earth, 112(B4), doi: 10.1029/2006JB004629.
[6] Hainzl S, Steacy S, Marson D. 2010. Seismicity models based on Coulomb stress calculations. Community Online Resource for Statistical Seismicity Analysis, doi: 10.5078/corssa-32035809. Available at http://www.corssa.org.
[7] Hainzl S, Marsan D. 2008. Dependence of the Omori-Utsu law parameters on main shock magnitude: Observations and modeling. Journal of Geophysical Research: Solid Earth, 113: B10309, doi: 10.1029/2007JB005492.
[8] Helmstetter A, Shaw B E. 2006. Relation between stress heterogeneity and aftershock rate in the rate-and-state model. Journal of Geophysical Research, 111: B07304, doi: 10.1029/2005JB004077.
[9] Helmstetter A, Shaw B E. 2009. Afterslip and aftershocks in the rate-and-state friction law. Journal of Geophysical Research, 114: B01308, doi: 10.1029/2007JB005077.
[10] Hsu Y J, Simons M, Avouac J P, et al. 2006. Frictional afterslip following the 2005 Nias-Simeulue earthquake, Sumatra. Science, 312(5782): 1921-1926, doi: 10.1126/science.1126960.
[11] Kagan Y Y. 2004. Short-term properties of earthquake catalogs and models of earthquake source. Bulletin of the Seismological Society of American, 94(4): 1207-1228.
[12] Li C Y, Xu X W, Gan W J, et al. 2013. Seismogenic structures associated with the 20 April 2013 Ms7.0 Lushan earthquake, Sichuan province. Seismology and Geology (in Chinese), 35(3): 671-683.
[13] Liu B Y, Shi B P, Lei J S. 2013. Effect of Wenchuan earthquake on probabilities of earthquake occurrence of Lushan and surrounding faults. Acta Seismologica Sinica (in Chinese), 35(5): 642-651.
[14] Mignan A, Woessner J. 2012. Estimating the magnitude of completeness for earthquake catalogs. Community Online Resource for Statistical Seismicity Analysis, doi:10.5078/corssa-00180805. Available at http://www.corssa.org.
[15] Narteau C, Shebalin P, Holschneider M. 2002. Temporal limits of the power law aftershock decay rate. Journal of Geophysical Research: Solid Earth, 107(B12): ESE 12-1-ESE 12-14, doi: 10.1029/2002JB001868.
[16] Ogata Y. 1988. Statistical models for earthquake occurrences and residual analysis for point processes. Journal of the American Statistical Association, 83(401): 9-27.
[17] Omori F. 1894. On after-shocks of earthquakes. Journal of the College of Science, Imperial University of Tokyo, 7: 111-200.
[18] Ouillon G, Sornette D. 2005. Magnitude-dependent Omori law: empirical study and theory. Journal of Geophysical Research: Solid Earth, 110: B04306, doi: 10.1029/2004JB003311.
[19] Parsons T. 2002. Global Omori law decay of triggered earthquakes: Large aftershocks outside the classical aftershock zone. Journal of Geophysical Research, 107(B9): ESE 9-1-ESE 9-20, doi: 10.1029/2001JB000646.
[20] Perfettini H, Avouac J P. 2004. Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan. Journal of Geophysical Research, 109: B02304, doi: 10.1029/2003JB002488.
[21] Perfettini H, Avouac J P. 2007. Modeling afterslip and aftershocks following the 1992 Landers earthquake. Journal of Geophysical Research, 112: B07409, doi: 10.1029/2006JB0041399.
[22] Ripperger J, Mai P M. 2004. Fast computation of static stress changes on 2D faults from final slip distributions. Geophysical Research Letter, 31(18): L18610, doi: 10.1029/2004GL020594.
[23] Shen W H, Qiu Z, Shi B P. 2014. Strong ground motion simulation for the 2013 Lushan Mw6.6 earthquake, Sichuan, China, based on the inverted and synthetic slip models. Earthquake Science, 27(4): 377-389, doi: 10.1007/s11589-013-0057-5.
[24] Utsu T. 1961. A statistical study on the occurrence of aftershocks. Geophysical Magazine, 30: 521-605.
[25] Utsu T, Ogata Y, Ritsuko S, et al. 1995. The centenary of the Omori formula for a decay law of aftershock activity. Journal of Physics of the Earth, 43(1): 1-33.
[26] Vidale J E, Cochran E S, Kanamori H, et al. 2003. After the lightning and before the thunder; non-Omori behavior of early aftershocks? Eos Trans AGU Fall Meet Suppl Abstract S31A-08, 84(46).
[27] Wang W M, Hao J L, Yao Z X. 2013. Preliminary result for rupture process of Apr. 20, 2013, Lushan earthquake, Sichuan, China. Chinese Journal of Geophysics (in Chinese), 56(4): 1412-1417, doi: 10.6038/cjg20130436.
[28] Xu X W, Chen G H, Yu G H, et al. 2013. Seismogenic structure of Lushan earthquake and its relationship with Wenchuan earthquake. Earth Science Frontiers (in Chinese), 20(3): 11-20.
[29] Zhang P Z, Xu X W, Wen X Z, et al. 2008. Slip rates and recurrence intervals of the Longmen Shan active fault zone and tectonic implications for the mechanism of the May 12 Wenchuan earthquake, 2008, Sichuan, China. Chinese Journal of Geophysics (in Chinese), 51(4): 1066-1073, doi: 10.3321/j.issn:0001-5733.2008.04.015.
[30] Zhang Y, Xu L S, Chen Y T. 2013. Rupture process of the Lushan 4.20 earthquake and preliminary analysis on the disaster-causing mechanism. Chinese Journal of Geophysics (in Chinese), 56(4): 1408-1411, doi: 10.6038/cjg20130435.
[31] Zúniga F R, Wyss M. 1995. Inadvertent changes in magnitude reported in earthquake catalogs: Their evaluation through b-value estimates. Bulletin of the Seismological Society of American, 85(6): 1858-1866.
[32] 杜方, 龙锋, 阮祥等. 2013. 四川芦山7.0级地震及其与汶川8.0地震的关系. 地球物理学报, 56(5): 1772-1783, doi: 10.6038/cjg20130535.
[33] 李传友, 徐锡伟, 甘卫军等. 2013. 四川省芦山Ms7.0地震发震构造分析. 地震地质, 35(3): 671-683.
[34] 刘博研, 史保平, 雷建设. 2013. 汶川地震对芦山地震及周边断层发震概率的影响. 地震学报, 35(5): 642-651.
[35] 王卫民, 郝金来, 姚振兴. 2013. 2013年4月20日四川芦山地震震源破裂过程反演初步结果. 地球物理学报, 56(4): 1412-1417, doi: 10.6038/cjg20130436.
[36] 徐锡伟, 陈桂华, 于贵华等. 2013. 芦山地震发震构造及其与汶川地震关系讨论. 地学前缘, 20(3): 11-20.
[37] 张培震, 徐锡伟, 闻学泽等. 2008. 2008年汶川8.0级地震发震断裂的滑动速率、复发周期和构造成因. 地球物理学报, 51(4): 1066-1073, doi: 10.3321/j.issn:0001-5733.2008.04.015.
[38] 张勇, 许力生, 陈运泰. 2013. 芦山4. 20地震破裂过程及其致灾特征初步分析. 地球物理学报, 56(4): 1408-1411, doi: 10.6038/cjg20130435.
[39] 中国地震台网中心. 2014. http://www.csndmc.ac.cn/newweb/catalog_direct_link.htm.