0 引言

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

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

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

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

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

(1)

(2)

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

(3)

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

(4)

(5)

(6)

在公式(4)，(5)和(6)中，R₁表示剪应力加载速率为时的地震发生率，即背景地震发生率.R₂表示剪应力加载速率为时的地震发生率，即前震发生率.R₃表示剪应力加载速率由恢复为后的地震发生率，即主震后余震的地震发生率.前震开始的时刻为t=0，主震发生的时刻为t₁.，表示地震成核的特征时间，也表示了主震后余震的持续时间.，表示剪应力加载速率为时前震成核的特征时间，t′_a与t_a的关系可写为：.从公式(4)中可以看出，前震发生率R₂主要取决于背景地震发生率R₁，与t_a的大小以及前震的持续时间t₁.

若t₁持续时间非常短，即t₁ < t′_a，且由主震成核引起的前震成核区剪应力上升率远大于由远场加载引起的剪应力上升率，即≫，在此近似条件下，，其中为t₁时刻断层滑移速率，T₁/t_a→0，exp(T₁/t_a)≈1+T₁/t_a.由附录公式(A19)与(A21)可知前震及余震发生率可近似为

(7)

(8)

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

(9)

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

(10)

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

(11a)

(11b)

(11c)

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

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

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

(12)

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

(13)

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

(14)

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

(15)

其中,

(16)

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

(17)

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

(18)

其中a₁与b₁表示在作用下a与b的取值，a₂与b_r表示在作用下a与b的取值.令f=b₁/b_r，且f·b_r≠1.5，当t₁≪t′_a时，假设a₁近似等于a₂，即地震产生率相同，则(18)式可写为

(19)

公式(19)建立了研究区内剪应力加载速率的变化与其同时域内b值的变化之间的关系.如图 2所示，当系数f小于1，即b₁ < b_r时，f的取值越小的比值越高，表明了前震期间前震区内剪应力加载率的升高与b值变小相关.前人的研究结果确定了应力对b值的影响(Schorlemmer et al., 2005)，低b值往往代表高应力，并表明了b值可作为推断应力的方法来使用(Narteau et al., 2009).而本文的结果表示剪应力加载速率的变化同样会影响b值的变化，若已知研究区域内发生过的最大震级m_max，背景地震的b值以及前震期间b值，则可利用本文的方法估算前震期间剪应力加载速率相对于背景场剪应力加载速率的升高量.

2 结果分析

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

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

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

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

3 讨论

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

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

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

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

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

4 实例分析

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

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

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

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

5 结论

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

附录  前震发生率公式推导

将RSF定律与单自由度弹簧-滑块模型相结合(Dieterich, 1992, 1994)可得到如公式(2)所示的断层滑移速率函数.当单个断层处于成核-自加速过程时，不同时域内断层的滑移速率为

(A1)

(A2)

(A3)

假定t=0为前震开始时刻，t₁为前震结束主震发生时刻，为t=0时断层滑移速率，为t₁时刻断层滑移速率，为背景场下断层滑移速率，为前震断层滑移速率，为前震结束后余震断层滑移速率.不同时域内地震发生所尚需的时间为(Dieterich, 1994)：

(A4)

(A5)

(A6)

在公式(A4)，(A5)和(A6)中，T₁表示以初始速度在作用下背景地震发生所需要的时间，T₂表示以初始速度在剪应力加载速率由升高到后，地震发生所需要的时间，T₃表示以初始速度在剪应力加载速率由又恢复到后，地震发生所需的时间.联立公式(A4)与(A5)有：

(A7)

根据Dieterich模型所述的地震群体化模型(Dieterich, 1994; Kaneko and Lapusta, 2008)，如附图A1所示，设定任一时间段Δt_i=ΔT′_i+1－ΔT′_i内地震发生的个数只同Δt_i的长短相关，即地震发生频率始终为一常数值，则可将这种状态定义为背景场，由R₁来表示.当研究区域内未受到应力扰动时，该区域的背景地震发生率R₁保持不变，于是有：

(A8)

其中n_i+1为T′_i到T′_i+1时间段内地震发生的个数，T′_i为背景场下断层失稳的时间(附图A1a).而当研究区域内不同尺度的断层受到应力扰动影响后，每个断层的失稳时间均发生变化，这一现象类似于Coulomb应力加载模型中的时钟提前或推后(ΔCFF/，其中ΔCFF为Coulomb应力变化).设定受应力扰动后断层失稳的时间变为F(T′_i)，则新的地震发生率R可以表示为(附图A1b)：

(A9)

将公式(A9)代入公式(A8)中，可得到与背景地震发生率相关的受应力扰动影响后的地震发生率为：

(A10)

而当T'_i+1－T'_i趋近于无穷小时，其微分解可表示为：

(A11)

则：

(A12)

公式(A3)可写为

(A13)

当地震发生时，→∞，则：

(A14)

由公式(A4)得：

(A15)

进而有：

(A16)

则：

(A17)

其中R₂表示前震发生率，R₃表示主震后的余震发生率.若t₁持续时间非常短，即t₁≈Δt，且主震成核远大于由板块构造运动对前震成核区的影响，即≫.在此近似条件下，，则公式(A7)可写为

(A18)

进而可得：

(A19)

同理，T₃/Aσ≈1/Hσ，则：

(A20)

进而可得：

(A21)

若t₁时刻主震发生后对研究区域产生的应力扰动为Δτ，并且研究区域内剪应力加载速率恢复到，则t₁时刻断层滑移速率可表示为(Dieterich, 1992, 1994)：

(A22)

将公式(A22)代入公式(A3)中，可得断层失稳时有：

(A23)

将公式(A4)代入公式(A23)中并整理可得施加应力扰动Δτ后余震发生率R₃为

(A24)

若前震序列中产生中等强度地震，比如1975年海城M_S7.3地震前震序列中的M_S5.0地震，并对前震成核区造成的应力扰动为Δτ₁，而前震成核仍受主震成核影响时，断层滑移速率可表示为

(A25)

(A26)

其中t′₁为前震序列中应力扰动Δτ₁产生时刻，为t′₁时断层滑移速率，为主震之前，应力扰动Δτ₁产生后断层滑移速率.断层失稳时有：

(A27)

将公式(A4)代入公式(A27)中并整理可得：

(A28)

其中T_2b为应力扰动Δτ₁产生后前震序列中地震发生尚需时间，R_2b为应力扰动Δτ₁产生后后续前震序列的地震发生率.进一步考虑主震发生后产生的应力扰动Δτ对后续余震序列的影响，则断层滑移速率可表示为

(A29)

(A30)

其中′_1b为主震发生时断层的滑移速率.同理，将公式(A4)代入(A30)式并整理可得：

(A31)

(A32)

(A33)