引言

地震研究的终极目标是对地震事件的发生时间、位置和震级进行准确的预测.目前,地震短期预测仍是具有高复杂度的世界性科学难题,现有的预测方法仍存在很大的可创新空间(周硕愚等,2019).对于地壳内部的构造地震而言,地震事件的发生是整体断层系统长期演化的结果.因此,地震预测方法的发展离不开对发震断层内部的摩擦过程及触发机制的研究.随着观测技术和岩石物理学实验的深入开展,人们对断层演化物理机制的研究取得了深入的认识,表明地震的演化过程就是断层在不同时空状态下克服摩擦的过程.基于近三十年实验结果所得到的摩擦本构关系,即速率-状态摩擦定律(Rate-andState-DependentFrictionLaw,以下简称RSF)为定量化描述断层内部的复杂属性提供了一个基本框架(Dieterich,1972,1978,1979a,1979b;Ruina,1980,1983;Rice,1983).前人应用RSF的数值分析成功地解释和重现了地震演化周期的各个阶段(Dieterich,1992,1994;RubinandAmpuero,2005;张龙等,2013),因此,RSF已成为了解地震成因机制和断层时空演化的重要依据.

前震泛指主震之前较短时间内,在主震附近发生的中小地震.而主震发生之前所出现的一系列前震活动则构成了整个的前震序列(周少辉和蒋海昆,2016).本文所讨论的前震是与主震发生时间、空间位置接近,和主震发生有着直接关系的“前震”(foreshock).在所有的地震所出现的短临前兆中,前震是公认最有效的强震预测指标之一(周少辉和蒋海昆,2016).例如1975年海城MS7.3地震的成功预报主要依靠于前震序列所提供的时间、空间和地震频次信息(Wangetal.,2006).因此,对前震成因机制的研究可为地震短临预测提供重要的参考.利用统计学方法,前人对前震活动进行了系统性的总结,其中释放加速的小震活动和古登堡-里克特定律中b值变小是前震序列最典型的两个特征(周少辉和蒋海昆,2016).从力学成因机制上考虑,统计规律中得到的这两个特征是否拥有共性仍需要进一步分析探讨,因此,了解其力学成因机制是尤为重要的.描述前震序列加速发生的经验模型可表示为N⋅=α t−ω,其中,N⋅为前震发生率,α为常数,ω为接近于1的常数,t为距主震发生的时间(Scholz,2002).经验模型明确地表达出了前震发生的加速特征,但未能解释前震的力学机制.Popov(2010)则认为前震是整个地震成核和快速蠕动过程中的一部分.基于RSF,Popov(2010)给出了前震发生率N⋅可近似于A/(Hlt),其中,t为距主震发生时间,H=B/Dc-k/σ,σ为断层上有效正应力,A和B为RSF中的摩擦参数,Dc为临界滑动位移,k=Gη/l,G为剪切模量,l为断层的特征尺度,η为描述裂纹几何状态的参数,一般近似于1.上述近似表达式虽然从地震成核的角度量化了时间与前震发生率的关系,表达出前震发生率随着主震的临近而逐渐升高,但仍不能完整地描述前震序列的基本活动性特征和包括“前震-主震-余震”在内的整个地震序列时空变化的特征.例如,无法反映出类似1975年海城MS7.3地震前震序列中大震前所出现暂时“平静”的现象(王琳等,2013).

基于RSF,本文将进一步探讨前震产生的力学成因机制,利用Dieterich(1992,1994)所提出的断层滑移函数研究前震序列中地震活动性增加,以及一些前震序列在大震前出现短暂“平静”现象的可能的原因,并探讨前震序列中古登堡-里克特定律中b值变小与前震序列加速发生所存在的可能的共性.

1、方法

1.1基于速率-状态摩擦定律分析前震发生的力学机制

基于RSF,对断层采用单自由度弹簧-滑块模型近似,当断层处于成核阶段并进入自加速过程时,断层上的滑移δ和滑移速率δ⋅随时间t的变化可表示为(Dieterich,1992,1994):

δ=−AHln{δ⋅0Hστ⋅[1−exp(τ⋅tAσ)]+1},(1)

δ⋅={[1δ⋅0+Hστ⋅]exp(−τ⋅tAσ)−Hστ⋅}−1,(2)

其中,H=B/Dc-k/σ,k=Gη/l,σ为作用于断层上的有效正应力,A和B为RSF中的摩擦参数,Dc为临界滑动位移,G为剪切模量,l为断层的特征尺度,η为描述裂纹几何状态的参数,远场剪应力加载速率τ⋅=kVpl,Vpl为远场对断层的加载速率,δ⋅0为t=t0时断层的滑移速率.在地震自加速过程中,当瞬态滑移速率接近厘米每秒或米每秒则表示断层出现失稳,相比于远场对断层的加载速率为厘米每年而言,此刻的瞬态滑移速率可视为无穷大(Dieterich,1994),即1/δ⋅≈0.从公式(2)中可知,断层失稳尚需的时间T可写为:

T=Aστ⋅ln(τ⋅Hσδ⋅0+1).         (3)

在描述前震成因机制的模型中,最具有代表性的是Dieterich基于RSF所提出的前震模型-Ⅰ与前震模型-Ⅱ(DieterichandKilgore,1996).前震模型-Ⅰ表示先前地震所产生的应力扰动改变了包括主震在内所有后续地震发生的可能性,即应力扰动影响断层的后续成核和失稳.在这种情况下每一个前震的发生都受先前地震产生的应力变化的影响,而这个地震发生后所产生的新应力扰动又继续影响着后续的地震发生.前震模型-Ⅱ则表示了主震成核区域的断层自加速滑移触发了前震,即主震成核过程中的主震断层自加速滑移驱动了前震加速成核(δ⋅≫Vpl).在这个模型中,主震成核所产生的引起前震发生的应力变化可以远大于远场加载的应力变化,而前震的发生并不影响到主震的发生.事实上,根据前震模型-Ⅰ得到的地震发生率就是应力阶跃后的地震发生率模型,相应的公式可以从Dieterich(1994)地震发生率模型给出.但该公式所呈现的地震发生率则随着时间而逐步衰减,无法描述出前震序列中地震发生率逐渐上升的过程.因此,若前震物理模型满足前述前震模型-Ⅱ,则前震成核区剪应力随着主震的临近而持续增加,其增加的幅度可远大于由远场加载速率Vpl所导致的剪应力增加.假设受主震和远场影响下的前震区剪应力的增加均是呈线性的,由远场引起的区域内剪应力加载速率为一常数τ⋅r,而由主震成核引起的剪应力加载速率为一常数τ⋅1.在前震孕育(主震成核)时,前震区剪应力加载速率由主震成核引起的剪应力τ⋅1主导.当断层失稳,即主震发生后,区域的剪应力由远场加载速率主导,即剪应力加载速率恢复至τ⋅r(如图1).我们可得到前震及余震发生率的解析解(详细的推导过程如附录中所示),由附录公式(A12)与(A17)可知:

图1(a)剪应力随时间变化示意图;(b)剪应力加载速率随时间变化示意图

R2R1=1(1−τ⋅rτ⋅1)exp(−tt′a)+τ⋅rτ⋅1,0

R3R1=11+Cexp[−(t−t1)ta],t1

C=(1−τ⋅rτ⋅1)exp(−t1t′a)+τ⋅rτ⋅1−1.(6)

在公式(4),(5)和(6)中,R1表示剪应力加载速率为τ⋅r时的地震发生率,即背景地震发生率.R2表示剪应力加载速率为τ⋅1时的地震发生率,即前震发生率.R3表示剪应力加载速率由τ⋅1恢复为τ⋅r后的地震发生率,即主震后余震的地震发生率.前震开始的时刻为t=0,主震发生的时刻为t1.ta=Aστ⋅r,表示地震成核的特征时间,也表示了主震后余震的持续时间.t′a=Aστ⋅1,表示剪应力加载速率为τ⋅1时前震成核的特征时间,t′a与ta的关系可写为:t′a=taτ⋅rτ⋅1.从公式(4)中可以看出,前震发生率R2主要取决于背景地震发生率R1,τ⋅rτ⋅1与ta的大小以及前震的持续时间t1.

若t1持续时间非常短,即t1

R2R1=exp(tt′a),0

主震的出现代表前震序列的结束,震源区的剪应力加载速率恢复到背景水平τ⋅r.若考虑到主震同震位移产生的静态应力扰动量为Δτ,根据RSF中的Dieterich状态演化方程,静态剪应力Δτ扰动前后滑移速率δ⋅a和δ⋅b满足δ⋅b=δ⋅aexp(Δτ/Aσ),由附录公式(A24)进而得受主震应力扰动Δτ后余震发生率为

R3R1=11+[(1−τ⋅rτ⋅1)exp(−ΔτAσ)exp(−t1t′a)+τ⋅rτ⋅1exp(−ΔτAσ)−1]exp[−(t−t1)ta],t1

若前震序列中小震发生时对震源区产生的应力扰动相对于主震成核对前震震源区造成的应力升高可以忽略,则利用前震模型-Ⅱ得到的公式(4)便可以描述前震序列发生率.但若前震序列中出现中强震时,比如1975年海城MS7.3地震前震序列中出现的MS5.0地震,此时MS5.0地震的发生对震源区造成的应力扰动不应该被忽略.假设前震序列中中强震的发生对震源区产生应力扰动Δτ1,如图1所示,则由附录公式(A28)可知受此应力扰动后前震发生率变为

R2bR1=1τ⋅rτ⋅1+{(1−τ⋅rτ⋅1)[exp(−Δτ1Aσ)exp(−t′1t′a)]−τ⋅rτ⋅1+τ⋅rτ⋅1exp(−Δτ1Aσ)}exp[−(t−t′1)t′a],t′1

其中R2b为应力扰动Δτ1产生后前震序列的地震发生率,t′1为Δτ1产生时刻.在这种情况下,主震发生并对震源区产生应力扰动Δτ后,由附录公式(A31)、(A32)和(A33)可知余震发生率可表示为

R3R1=11+[(1−τ⋅rτ⋅1)C1+τ⋅rτ⋅1exp(−ΔτAσ)C2−1]exp[−(t−t1)ta],t1

其中,

C1=exp(−ΔτAσ)exp(−Δτ1Aσ)exp(−t′1t′a)exp[−(t1−t′1)t′a],         (11b)

C2=1+exp[−(t1−t′1)t′a][exp(−Δτ1Aσ)−1].(11c)

从公式(10)中可以看出,受应力扰动后的前震发生率主要取决于τ⋅rτ⋅1,Δτ1Aσ和ta的大小,t1的持续时间以及t′1发生的时刻.从公式(9)与公式(11)中可以看出,影响余震发生率的物理量除了包括所有影响前震发生率的物理量外,还取决于ΔτAσ的大小.

1.2前震b值变化的力学机制

古登堡-里克特震级-频度关系(GutenbergandRichter,1944)可以表示为lgN(≥m)=a-bm,其中,N(≥m)表示震级大于等于m的地震总个数,a值代表所选研究区的地震活动总体水平,b值同区域介质特性、应力状态和不均匀性等有关,能反映所研究区域地震的震源特性(Rundle,1989;MarzocchiandSandri,2003).地震数量的增量表示成震级的函数为

n(m)=−dNdm=bln(10)⋅10a−bm,         (12)

其中n表示对于无穷小的震级增量所对应的地震个数增量,此时N的单位为个数/年.累积的地震矩率可以表示为

M⋅=∫mmax−∞M0(m)n(m)dm,         (13)

其中,M0(m)表示震级为m的地震的地震矩,mmax表示所研究序列的最大矩震级.地震矩与震级的关系可近似为M0(m)=10p+q·m(KanamoriandAnderson,1975),其中p与q为常数,q一般取值为1.5,m代表矩震级.根据Kostrov模型(Scholz,2002),形变率ε⋅与地震矩M0的关系可表示为

ε⋅ij=12G1ΔVΔT∑k=1NMk0ij,         (14)

其中ΔV为研究区体积,ΔT为时间,G为剪切模量.若区域内发震断层的类型基本近似,则剪应力加载率τ⋅可表示为

τ⋅=1ΔVM⋅,         (15)

其中,

M⋅=1ΔT∑k=1NMk0ij,         (16)

将公式(15)与(12)和(13)联立可得:

τ⋅=1ΔVb1.5−bM0(mmax)10a−b⋅mmax.         (17)

若在τ⋅1影响下研究区的最大震级与在τ⋅r影响下保持一致,即区域内最大震级保持不变,则进一步可得研究区内τ⋅1与τ⋅r之比为

τ⋅1τ⋅r=b1brq−brq−b110a1−ar⋅10(br−b1)mmax,         (18)

其中a1与b1表示在τ⋅1作用下a与b的取值,a2与br表示在τ⋅r作用下a与b的取值.令f=b1/br,且f·br≠1.5,当t1≪t′a时,假设a1近似等于a2,即地震产生率相同,则(18)式可写为

τ⋅1τ⋅r≈f(1.5−br)1.5−fbr⋅10br(1−f)mmax.         (19)

公式(19)建立了研究区内剪应力加载速率的变化与其同时域内b值的变化之间的关系.如图2所示,当系数f小于1,即b1

图2τ⋅1/τ⋅r随系数f的变化

2、结果分析

利用前震模型-Ⅱ得到的前震发生率公式(4)与余震发生率公式(5),本文得到如图3a所示的在不同持续时间t1下地震发生率特征.当t1≪t′a时,前震发生率R2随时间的增加而快速上升;当t1>t′a时,前震发生率R2便维持在一稳定值,此时R2≈R1τ⋅1τ⋅r;而当t>t1时,表示主震发生后的余震序列发生阶段,在此阶段内τ⋅1τ⋅r=1,地震发生率R3在此阶段初期迅速下降,在t～ta时逐渐回归到背景场R1的水平.若前震发生率R2还未升高到R1τ⋅1τ⋅r这一值时余震便开始发生,例如图3a中t1=0.005ta时,则整个地震序列表现出主震前有加速发生的前震,而主震后余震发生率逐渐降低,与统计模型所得到的结果相近(周少辉和蒋海昆,2016).

从近似公式(7)和(8)中可以看出,前震发生率的近似解同解析解(4)和(5)所展现的不同之处在于,近似解所表现出的前震发生率R2在主震发生之前是一直以e指数上升的趋势升高,并且主震发生后余震发生率R3最终会趋于0(图3b),这是因为近似解中忽略了τ⋅r的影响所致.但当t1≪t′a,t≪ta时,近似解与解析解所表现出的地震发生率是相近的.

结合公式(4)和(9)可以得到如图3c所示的前震-主震-余震型序列地震发生率随时间演化特征,当主震对震源区所产生应力扰动Δτ>0时,余震发生率如图3c中紫色线所示,即在t1时刻地震发生率陡增,随后余震发生率R3逐渐下降,这种地震发生率分布与Scholz(2002)书中所描述的前震-主震-余震型地震序列示意图(图4a)完全相似.当Δτ=0时,公式(4)与(9)所表现出的地震发生率特征同图3a完全相同,即在t1时刻地震发生率是连续的.而当应力加载Δτ<0时,t1时刻的地震发生率会出现较大的下降,但如图3c中所示,余震发生率R3仍高于背景地震发生率R1且回归到背景地震发生率所需要的时间尺度为ta,与应力扰动Δτ>0时余震发生率回归到背景场水平所需的时间尺度相同.与前人研究结果不同的是,Dieterich(1994)地震发生率模型给出,当应力加载Δτ<0时便会有地震影区的出现,即t1时刻地震发生率R3便小于背景地震发生率R1,随后R3逐渐上升并经历ta的时间尺度恢复到R1水平.而本文的结果表明应力加载Δτ<0时并非一定会出现地震影区,也有可能出现地震发生率逐渐下降的余震序列(图3c).只有当负向应力扰动Δτ足够大时,本文的结果才会出现地震影区.

在前人的研究中,前震模型-Ⅰ与前震模型-Ⅱ是相互独立的两个模型(DieterichandKilgore,1996),而本文将这两个模型结合并进行了定量化研究进而得出公式(4)、(10)与(11).当Δτ1=0时,这组公式所表现出的就是前震模型-Ⅱ,而当Δτ1≠0时,即在前震模型-Ⅱ所蕴含的力学机制基础上,考虑了前震模型-Ⅰ所蕴含的应力变化对前震序列演化特征的影响.如图3d所示,前震序列在t′1时刻前地震发生率特征同图3a中R1相同,当Δτ1>0时,t′1时刻地震发生率陡增到一定值后,后续的前震发生率R2b开始下降,在这种情况下整个前震序列演化特征类似于图4a所示的完整的前震-主震-余震型地震序列,但这段序列其实全部是由后面t1时刻发生的主震的成核所引起.当Δτ1<0时前震序列在t′1时刻地震发生率陡降,随后后续的前震序列地震发生率R2b再次逐渐上升.若|Δτ1|相对较大,如图3d中Δτ1/Aσ=-8这条曲线所示,此时t′1时刻的前震发生率R2与背景地震发生率R1之比小于0.1,即表明前震在此时“静止”,随后后续的前震序列R2b逐渐升高.若主震发生在“静止”期之间,则整个地震序列演化特征表现为:逐渐加速的前震→地震发生率平静→主震→余震,这与1975年海城MS7.3地震前震序列的演化特征十分相似(吴开统等,1976;王琳等,2013).若主震发生在前震“静止”期之后,则整个序列演化特征表现为:加速发生的前震→平静→加速发生的前震→主震→余震.图3d中Δτ1/Aσ=1,Δτ/Aσ=2这条曲线表示t′1时刻前震序列中产生的应力扰动Δτ1>0,且t1时刻主震发生并产生应力扰动Δτ>0时对余震发生率R3的影响,表明在这种情况下R3在t1时刻出现了陡增后,后续的余震发生率逐渐降低,并最终趋于背景地震发生率R1的水平.当Δτ<0时,本文例举Δτ1/Aσ=1,Δτ/Aσ=-2情况下余震的发生率,在这种条件下R3在t1时刻陡降后,后续的余震发生率逐渐降低并最终趋于R1的水平,若|Δτ|足够大,则R3可能在t1时刻下降到R1之下,并经过ta的时间尺度回归到R1的水平,即出现主震后余震的影区.

图3地震发生率随时间演化特征示意图

图4地震序列类型的示意图（引自Scholz,2002．图4.22)

3、讨论

慢地震是指以几个小时到几个月的时间尺度缓慢地发生了有限的断层位错,但没有高频地震波辐射的断层慢滑动事件(吴忠良和许忠淮,2013).前人已将RSF应用于慢地震的研究中(Segalletal.,2010;Ikarietal.,2013),其中Segall等(2006)在研究夏威夷Kilauea火山地区慢滑移触发地震时,证明了慢滑移所引起的应力变化的增加(τ⋅´1/τ⋅r≫1,其中τ⋅´1为由慢滑移导致的剪应力加载速率)将导致火山地区地震活动性的增加.慢滑移触发地震模型与前震模型-Ⅱ具有相似性,前者表明慢滑移导致断层处的剪应力变化持续增加,而后者则表明主震成核驱动震源区次级断层上的剪应力变化量增加.值得注意的是,Segall等(2006)所给出的由慢滑移影响所导致的地震发生率公式和本文所推导出的无Coulomb应力扰动的前震和余震发生率公式(4)和(5)具有相同的形式.这是因为慢滑移触发地震模型与前震模型-Ⅱ所蕴含的力学机制相似,并且在描述剪应力随时间增加函数时,Segall等(2006)也将其简化成τ=τ0+τ⋅t的形式,τ0为t=0时的剪应力大小.而Segall等(2006)的研究与本文的不同之处在于推导公式的途径与阐述的物理过程不尽相同,前者利用的是Dieterich(1994)地震发生率公式进行推导,地震发生率变化是由慢滑移过程所导致,而且这样的慢滑移一般来自于孕震层的下部,而本文则是依据断层滑移函数并针对前震的成因过程得到的这组公式.前震模型-Ⅱ与由慢滑移触发地震的机制都能表现出地震发生率增加这一现象,但由慢滑移所触发的地震序列常表现出无较大震级的群震,如图4c所示,而由主震成核驱动的前震序列往往会伴随较大震级的主震.主震成核所导致的剪应力加载速率与背景场剪应力加载速率之比τ⋅1/τ⋅r可能远大于由慢滑移所导致的τ⋅´1/τ⋅r,例如下文中拟合的中国台湾花莲地区2018年2月4日前震序列τ⋅1/τ⋅r≈3200,与由Segall等(2006)拟合的夏威夷Kilauea火山地区慢滑移引起的τ⋅´1/τ⋅r≈33相比约大两个数量级.当出现一段逐渐加速发生的地震序列时,由于慢滑移触发的地震与本文给出的未受应力扰动的前震的地震发生率均随时间逐渐上升,所以目前仍无法根据地震目录的特征迅速判断该序列是由慢滑移还是后续主震成核引起.若能够在加速的地震序列出现的初期判断此序列是群震还是前震,将会有效地减少资源浪费与人员伤亡,因此,未来还需对慢滑移与前震展开更细致的对比研究.

公式(4)是按照Dieterich和Kilgore(1996)前震模型-Ⅱ的力学机制进行的推导,但实际上如前震模型-Ⅰ所述,每一个地震的发生都会对后续的地震产生影响,只不过当这种影响相对于主震成核所造成的影响太小时可忽略,而当前震中出现相对较大震级地震时,就不应该忽略这次地震产生的应力扰动对后续地震的影响,也就表明此时需要将前震模型-Ⅰ的物理过程与前震模型-Ⅱ进行结合,由此本文推出前震公式(10).该公式表明若前震序列中出现的地震对整个震源区产生较大的应力扰动Δτ为负值,则有可能出现前震序列先密集后平静再主震的现象,这也是Scholz(2002)和Popov(2009)所给出N⋅∼t−1前震模型所无法表达的.实际上,前震序列中出现不可忽略的应力扰动Δτ也许不止一次,当考虑多次应力扰动对前震区的影响时,前震发生率在时间空间分布上呈现出的形式可能更为多样.因此,看似完整的一个地震序列(“前震-主震-余震”)也有可能是后续更大地震的一段前震序列,而如何识别不同触发机制的前震仍需要更为深入研究.若前震中出现的中强震产生的应力扰动Δτ足够大,那么其对整个的地震循环将会产生扰动并可影响主震发生的时间,当Δτ为正值时会缩短此次主震成核周期,导致主震提前发生,反之当Δτ为负值时则会导致主震发生时间推迟(解孟雨和史保平,2016).而Δτ的临界值与主震发生时间提前或推后量Δt的大小的确定及与地震预报融合,目前仍是一个难以解决的问题.另外,在给出地震发生所需时间T时(Dieterich,1992,1994),本文采用了Dieterich(1994)解中的近似假设,即地震发生时断层的滑移量和滑移速率趋于无穷大,但在实际地震中地震发生时断层滑移量和滑移速率应为一有限值,当滑移量和滑移速率为有限值的断层失稳时间T的解由Beeler和Lockner(2003)给出,但其对于现象的分析则更为复杂.

前人通过元分析(meta-analysis)方法对前震研究的结果表明,很多前震活动可能并未被检测到(Mignan,2014).最近,最新的高分辨率地震目录显示,在2008年至2017年间,南加利福尼亚州发生的大部分中强地震都伴有较高地震活动性的前震,而这些前震大多震级非常小,因此之前并未被地震网络探测到,这表明了自然界中的前震活动可能比以往认为的更普遍(TrugmanandRoss,2019).基于本文的方法,若区域内被观测到的背景地震发生率R1为0或者非常接近于0,则观测到的前震发生率R2同样为0或者接近于0,所以在实际应用中,本文的方法更适用于具有检测微小震级地震能力且背景地震活动丰富的地区.理论上,前震分布在主震成核所能够影响到的区域内,但这一区域的边界范围的计算方法暂时仍是一个未被解决的问题.

在Dieterich模型中,单自由度的弹簧滑块模型类比了远场对断层加载速率为Vpl的摩擦过程,即断层的黏滑机制.弹簧的有效刚度k表示为k=Gη/l,其中G为剪切模量,η为描述裂纹几何状态的参数,l为断层的特征尺度,剪切应力加载速率可以写为τ⋅=kVpl.如果设主震断层的刚度为km,特征尺度为lm,而次级断裂断层(前震断层)的刚度为kf,特征尺度为lf,将所发生地震的地震矩由圆盘断层模型表示,则地震矩可写为Mo=167Δσr3,其中Δσ为应力降,r为圆盘断层模型的半径(Brune,1970,1971),也可表示为断层的特征尺度.针对常数应力降模型(Δσ保持不变),由矩震级MW同地震矩的关系可得到(Kanamori,1977;HanksandKanamori,1979):rmrf=100.5(Mmw−Mfw),其中Mmw为主震矩震级,Mfw为前震矩震级,若主震与前震矩震级相差2,则主震断层尺度rm与前震断层尺度rf之比为10,若主震与前震矩震级相差4,则主震断层尺度rm与前震断层尺度rf之比为100.若主震断层与前震断层的剪切模量G与裂纹几何状态η相等时,有kf/km=rm/rf,进而当主震断层的特征尺度远大于前震断层的特征尺度(rm≫rf)时,前震断层与主震断层刚度的关系为:kf≫km.图5给出t/ta介于1.6至2.4之间时断层滑移速率δ⋅(δ⋅由公式(2)给出)的演化(Dieterich,1992,1994),其中,A和B的取值分别为0.01和0.02,有效正应力σ取值100MPa,Dc取值0.001m,Vpl取值10-9m·s-1,即Vpl≈3.2cm·a-1,δ⋅0取值10-10m·s-1,k/kc取值0.95,其中kc为临界刚度(kc=(B-A)σ/Dc).当t≈2.31ta时表示断层的失稳,而当t处于1.7ta～2.3ta之间时,δ⋅/Vpl的值≫1并且随着距主震发生时刻的临近而升高,结合前震断层与主震断层刚度关系(kf≫km),我们可以近似得到当断层处于成核过程中的自加速阶段时,前震区域内剪应力加载速率τ⋅1=kfδ⋅+kfVpl,从而有τ⋅1≫τ⋅r.所以,前震区域内剪应力加载速率升高同样可以从断层的特征尺度,主震断层滑移速率和Dieterich模型结合的角度来说明.

图5断层滑移速率示意图

古登堡-里克特定律中的a值为地震产生率,公式(19)更适用于背景地震和余震与前震a值相近的情况,即满足t1≪t′a.从公式(19)中可以明显看出b值的变化同剪应力加载速率变化的关系,与前震发生率公式相比较可以得出前震b值变小与发生率逐渐增加的共性可能是区域内剪应力加载速率的升高.Gulia和Wiemer(2019)的近期研究表明,震级大于MW6.0的主震的余震的b值通常比背景地震的b值大20%左右,而若MW6.0或更大震级的地震发生后伴随的“余震”序列的b值相对背景场b值变小或者不变,则可能发生更大震级的地震,并依此提出了FTLS(ForeshockTrafficLightSystem),以用于对地震风险的评估.Gulia和Wiemer(2019)认为前震序列b值减小也许是因为前震区域内应力的升高,而我们的成果可能会对前震序列中b值的减小的原因提供新的认识.

4、实例分析

2018年2月4日在我国台湾花莲东海岸米伦断层附近(23.83°N—24.39°N,121.38°E—122.1°E)发生了主震震级为ML5.5的前震-主震-余震型地震,ML5.5主震震中位置为24.18°N,121.73°E,震源深度为10km.在ML5.5主震发生2天后,即2018年2月6日该区域又发生主震震级为ML6.0的主震-余震型地震,ML6.0主震震中位置为24.14°N,121.69°E,震源深度为10km.自2月4日地震开始后至2月10日,该区域地震共造成17人死亡,295人受伤,4幢建筑坍塌和至少175幢建筑受损.台湾地处环太平洋地震带的西部,该地震带是菲律宾海板块与欧亚板块的边界构造带,菲律宾海板块以大约75mm·a-1的速度向西北方向俯冲到欧亚板块之下,产生强烈的逆冲和左旋剪切作用,形成海岸山脉和纵谷断裂(邓起东等,2002;USGS,2018),造成该地区地震活动极为频繁.而此次花莲地震序列是由菲律宾海板块和欧亚大陆板块之间的米伦断层左旋走滑所致(Kuoetal.,2018).本文的重点是研究2018年2月4日ML5.5地震序列,ML5.5地震前震、ML5.5地震及其余震的空间与时间分布如图6所示.从图6b与图6c中可以看出ML5.5地震前震发生率逐渐升高,ML5.5地震发生后地震频率开始降低.从图6a中可以看出ML5.5前震集中在ML5.5地震北方且处于ML5.5余震空间分布范围之内,ML5.5余震分布围绕在ML5.5地震10km内,其东北至西南走向展布约20km,西北至东南走向展布约10km.

图7表示了花莲ML5.5地震发生前2h的前震及其震后2h的余震地震发生率与时间的关系,单位为个数/小时.黑色圆圈代表其x轴对应的时间点前后半小时中研究区内发生地震的总个数,即实际地震发生率,单位为个数/小时.黑色实线代表利用本文地震发生率公式(4)与(5)拟合情况,黑色虚线表示利用发生率近似公式(7)与(8)的拟合情况.花莲地区背景地震平均每4天发生一个,即背景地震发生率R1约为0.0104个/小时.将t1设置为2h,则从公式(4)、(5)可知待拟合参数为τ⋅1与τ⋅r的比值与特征时间ta,而从近似公式可知待拟合参数只有t′a.

图6台湾地形图及花莲ML5.5地震序列时空分布图

图7ML5.5前震与余震发生率图

拟合结果表明,当ta取值730h,τ⋅1/τ⋅r取值3200时,利用公式(4)、(5)拟合的曲线与实际情况最接近,此时均方根误差为1.95.当t′a取值0.26h时,通过近似公式(7)、(8)拟合的曲线与实际情况误差最小,此时均方根误差为2.85.通过比较均方根误差的大小可以表明解析表达式的拟合情况要优于近似表达式,但若前震持续时间较短时,利用近似表达式拟合却更为简洁,其待拟合的参数只有t′a一项.

利用古登堡-里克特震级-频度关系对花莲地震分析,取完备震级Mc=3.5后利用最大似然法计算得到研究区内a值与b值的变化如表1和图8所示.从中可以看出,背景地震和ML5.5余震序列的b值皆接近于1.0,ML5.5前震序列的b值近似为0.5.利用本文公式(19),br取1.0,f取0.5.根据台湾气象局(2018)给出的“1901—2000的灾害性地震列表”及“2001迄今的灾害性地震列表”,自1901年以来台湾花莲地区地震的最大震级为1920年6月5日发生的震源深度为20km的ML8.3地震.所以研究区内最大震级取ML8.3,则可计算得到τ⋅1/τ⋅r≈3531,与采用公式(4)和(5)拟合得到的结果相近.

表1通过G-R定律拟合而得到的研究区内不同时间尺度下的a值与b值

5、结论

前震的发生可能由主震成核过程中断层自加速滑移所驱动,当前震成核区受到较大的应力扰动时,后续的前震序列受主震和这次扰动共同影响.本文由Dieterich(1992,1994)地震发生率模型给出了受主震成核影响下前震发生率的解析表达式与近似公式及受应力扰动影响时的前震发生率公式.解析表达式表明前震发生率由区域应力加载速率τ⋅1/τ⋅r,特征时间ta,前震的持续时间t1与背景地震发生率R1共同决定,前震发生率R2最初随着距主震时间的缩短逐渐升高,而后趋于稳定.本文给出的受应力扰动影响的前震发生率公式表明前震发生率的大小还受Δτ1/Aσ与应力扰动施加的时刻t′1决定,当应力扰动Δτ1为正值时,前震发生率会在t′1时刻出现陡增而后降低,但不会低于未受应力扰动时前震的发生率.当应力扰动Δτ1为负值时,前震发生率会在t′1时刻出现陡降而后升高,这一过程在实际地震目录中可表现为前震序列密集后的平静.b值降低是前震序列的典型特征(周少辉和蒋海昆,2016),本文给出剪应力加载速率的变化与古登堡-里克特定律中的b值变化的关系,表明前震序列中b值降低与前震区剪应力加载率升高有关,而前震发生率公式中剪应力的升高也是决定前震发生率大小的重要因素,所以,加速发生的前震与前震b值降低的共性可能是前震区剪应力加载速率的升高.

图8花莲地震序列发生后不同时间段内地震目录的震级-频度（M-F）图

尚园程,史保平.基于Dieterich地震发生率模型分析前震成因的力学机制[J].地球物理学报,2020,63(09):3355-3369.

基金:国家自然科学基金(41574040);国家国际科技合作专项(2015DFA21260)共同资助