统计地震学的两个基本定律对应力的依赖关系＊

C lém entNarteau Svetlana Byrdina Peter Shebalin Danijel Schorlemm er

余震衰减速率Λ由大森-宇津定律描绘,即Λ(t)=k(t+c)-p,这里 t是距离主震的时间,k是余震序列的发生率,p是幂律指数,c是幂律余震衰减速率起始前的延时。对于 c值的争议很大[4-10],研究人员[11-12]断定 c接近于零,大于零值就是 (至少部分是)由于主震后 (特别是主震尾波期间)数据不完整造成的人为结果。2004年加州帕克菲尔德M6.0地震[13]和 2004年日本新泻县中越M6.6地震[8]后的研究表明,主震后立即对微震进行仔细识别,估算出的 c值较小,但仍非零。我们对南加州和日本的研究结果表明,将分析范围限定在较小主震震级和较大余震震级会得出代表早期余震系统性特征的非零 c值,而不是一个完全的人为结果。

对于南加州,我们使用 1984—2003年有震源机制解[14]并经重新定位的地震目录和南加州地震台网 (SCSN)自 2003年以来的地震目录。我们只选用那些解的误差 ＜0.2、台站分布系数≥0.5的高质量地震事件。这一子目录包含 12000个事件。为识别主震,我们采用基于幂律关系的时间和空间窗算法[15-16],此关系反映了地震现象的分级特性。没有震源机制解的相应余震按相同规则取自全部 SCSN目录 (1984—2007年, ＞385000次事件)。

和南加州一样,我们对于日本也使用两个不同的地震目录。我们将 F-Net地震台网的目录与 JMA(日本气象厅)的目录相结合,分离出自 1997年以来主震 (F-Net)及其各自余震 (JMA)的震源机制。为剔除日本不同火山区的群发性地震活动,我们还去除了距主震时间的几何平均值超过 4小时的所有余震序列。因此也就去除了那些余震衰减速率不明显的群发性空间地震活动。

为了避免由重叠记录引起的人为结果,我们重点考虑由中等震级主震产生的余震序列,而不考虑大地震。同样,我们只考虑较大余震,并根据主震时间对其进行叠加处理,以弥补每一序列中地震数量的不足。尤其是必须分别确定主震和余震的两个震级范围:我们将主震的最小震级阈值设定为因为这是重新定位的目录的总体完整性水平[3]。然后,我们确定为重叠记录开始对南加州目录中M≥1.8余震的地震频度-震级分布产生系统性影响的主震震级 (见附 2:方法)。余震震级范围是一个可以改变从而产生不同 c值的自由参数。在所有情况下,我们都只选择完整性地方震级小于(余震的最小震级阈值[17])的余震序列。对于每一次主震,利用拟合优度检验法[18]对过去 3年间半径为 15 km的圆内这一完整性震级进行计算。由此,值不仅可以确定我们在单个序列中考虑的事件数量,而且可以确定被考虑的序列数量 (例如,南加州取则会在最初的 12小时期间产生近 2800次主震和 4000多次余震)。使用 A ki规则中的滑动角λ对主震的断裂模式进行分级,即走滑型事件λ≈0°或λ≈180°,逆冲型事件λ≈90°,正断层事件λ≈ -90°[3]。在此,我们既用滑动角也用滑动窗,其中滑动角的宽度为 60°,滑动窗的步长为 5°。这就产生了 73个主震等级。对于每一个等级和不同的震级区间,我们都根据余震距离各自主震的时间对余震进行叠加处理,并利用由模拟退火[4]实现的连续极小化和最大似然法对 c值进行估算。此外,我们还利用蒙特卡罗 (Monte Carlo)方法评估最大似然法估算 c值的不确定性,证实对第一个和第二个滑动角分别进行处理得出的观测结果都是一致的。我们也用优化程序的不同起始时间进行了一系列系统性的计算,以测试c值的稳定性,因此去除了与目录不完整相关的潜在的人为因素 (见附 2:方法)。

对于南加州,图 1(a、b)显示了 c值相对于余震震级分别为 1.8～2.8和 2.4～3.4的滑动角的演化。对于这些震级区间,图 1(c～h)显示了每一滑动角所选余震的数量以及所选主震和余震的平均震级。c值的不确定性对叠加中所用事件的数量有很强的依赖关系 (图 1c、d),但即使是在震级阈值很高的情况下也会出现信号。走滑型地震的 c值比逆冲断层型的高,但比正断层型的低。这意味着逆冲断层型地震的幂律余震衰减速率起始前的延时较短,正断层型的较长,走滑型的居中。两个震级范围的特性也是如此,只不过会出现一般性的偏移值,这可由 c值对余震平均震级的依赖关系来解释 (见附 2:方法)。逆冲型、走滑型和正断层型地震之间 c值的差别不能用余震平均震级的变化(即相应的 b值)来解释,因为当按照滑动角绘制时,它呈现的不是一个系统性模式 (图1e、f)。然而,两个震级范围的主震的平均震级都表现出类似特性 (图 1g、h)。与走滑型地震相比,逆冲断层型地震的主震较大,而正断层型地震的主震较小。这是南加州的普遍特点,可由地震频度-震级分布的差异解释。

图 1 南加州地区断裂模式对幂律余震衰减速率起始前延时的影响。(a、b)分别与余震震级范围(1.8,2.8)、(2.4,3.4)的滑动角相对应的 c值的对数。红、绿、蓝分别表示走滑型、正断层型和逆冲断层型事件的滑动角区间。误差棒由Mon te Ca r lo方法估算,相当于每个点上 500个合成余震序列的最大似然法估算 c值的 16%～84%的分位数。(c～h)主震和两个余震范围内余震的统计资料。(c、d)所选余震的数量(N)。(e、f)所选余震的平均震级。(g、h)所选主震的平均震级

对于日本,为了证明缺失的余震不会影响 c值的估算,我们也对涉及震级为 3.0～4.5的 77次主震的 JMA目录的完整性与Peng目录[7]的完整性进行比较。然后,正如在南加州一样,我们也发现相同的余震和主震震级范围内 c值对滑动角有着相同的依赖关系 (图 2a、b)。因为 F-Net目录只包含少数几次M＜3.5事件[3],所以所选主震的平均震级比南加州的要大。因此,估算的 c值总体较小。

图 2 在日本地区断裂模式对幂律余震衰减速率起始前延时的影响。(a、b)分别与余震震级范围(1.8,2.8)、(2.4,3.4)的滑动角相对应的 c值的对数。对 Peng[7]和 JMA目录进行比较可以发现,对于 1.8≤MA≤2.8的余震,优化程序在主震发生后 5×10-4天开始。误差棒由Mon te Car lo方法估算,相当于每个点上 500个合成余震序列的最大似然法估算 c值的 16%～84%的分位数

几乎在所有情况下,地震频度-震级分布都遵循一个以 b值 (Gutenberg-Richter频度-震级关系指数)为特性的幂律关系。对于用于识别日本和南加州主震的相同目录,b值会因断裂模式的不同而发生系统性变化。正如文献 [3]中所描述,走滑型事件的 b值(～0.9±0.05)确实比逆冲型事件的 (～0.7±0.05)高,而比正断层事件的 (～1.1±0.05)低。将两个观测结果结合起来看,c值和 b值都依赖于可由差异剪应力解释的震源机制,具体结果如下。

根据摩尔-库伦 (Mohr-Cou lom b)理论,均匀介质中受均匀三轴应力 (σ1,σ2,σ3)作用的剪切破坏应发生在最优取向平面上。这些平面与中等剪应力σ2平行,且与最大压应力σ1形成角Θ。因此,断裂运动状态受控于阿蒙顿 (Amonton)定律:

这里μ是静摩擦系数,σn和τ是在断层面上确定的正应力和剪应力。假设其中一个主应力是垂直的,另外两个位于一个水平面内,那么就有可能表示出相对于μ而言断层作用出现时的角Θ和差异剪应力 (σ1-σ3)。在此情况下,主要剪应力的取向不仅决定着断裂机制,而且也决定着差异剪应力的大小。例如,垂直应力对应的是正断层作用的最大压应力和逆冲断层作用的最小压应力。因此,对于相同的静摩擦系数,逆冲断层作用的差异剪应力要比正断层作用的高得多。在这些挤压和伸展状态之间,走滑断层作用会取 (σ1-σ3)的中间值。因此,c值与滑动角之间的关系表明余震衰减速率起始前的延时与差异剪应力呈显著的正相关。

考虑到余震源自主震附近应力的阶梯状扰动,有关余震产生的各种模型预示这种扰动的幅度控制着短时期内余震衰减速率的非幂律状态的持续时间。例如,我们可以看到一个增加的应力步长在速率-状态摩擦模型[20]、静疲劳模型[4,21]和破坏机理模型[22-23]中是如何加速幂律衰减速率的起始的。此外,经典裂纹扩展实验表明应力扰动的幅度是远场外加应力的递增函数[24]。图1a、b和图 2a、b中所表现出的 c值的特性印证了差异应力、应力重分布强度以及余震序列初始阶段之短暂性之间的这种关系。然而,个别余震序列却依然有高度分散,这极有可能是由应力不均匀性效应和测量不确定性引起。因此,c值对差异剪应力的依赖关系只能通过余震的叠加来描绘。

b值决定着小地震和大地震数量之间的比率,其变化可与成核状态下地震破裂能(低 b值 )否 (高 b值)传播相关。由此,不同应力状态下 c值与 b值演变之间的正相关在地震传播和余震产生期间都表现出相同的对应力的依赖关系。在需要较高的差异剪应力才能引起滑动的区域,孔隙度会降低,断裂带不同区段之间的相互作用增强,裂隙在不同方向上的大范围内蔓延速度更快[25-26]。于是,地震传播期间,在较高围限应力下的过程区内会呈现较短的破裂、滑移或破坏的特征时间尺度。因此,逆冲断层型地震有可能比正断层型地震传播得更远,这一结果与挤压状态下的 b值可能比伸展状态下的 b值高的观测结果一致。如果我们假设由于结构和组分的不均匀性而导致破裂终止,那么同样的推理依然有效,在高应力环境下使余震成核的时间也就更短。由此可更快地实现稳态松弛,逆冲断层型地震的 c值也就比正断层型的 c值更小。

地震频度-震级分布和余震衰减速率对应力的依赖关系为地震危险性分析和地震物理学提供了新的视点,因为它表明地震活动性的两种独立情况可用来约束脆性孕震地壳内的应力状态。例如,地震活动背景和余震发生的准系统性可能不仅有助于量化沿主要断层的应力积累[27],而且还可以量化一次潜在地震的特征震级[28]。

附 1:方法概述

为识别主震,我们没有选取那些震级小于M的地震,即在一次震级为M的地震事件发生后的 0.125×100.55M天、一个半径为0.02×100.50Mkm的圆内发生的地震[15-16]。利用同样的空间尺度标准,在较大事件前不到 12小时内发生的地震被视为可能的前震,并被从主震目录中删除。这样一来,在最初的 12小时内就不会出现重叠的余震序列。与所选主震相对应的没有滑动角的余震取自全部 SCSN目录 (1984—2007年, ＞385000次事件)和 JMA目录 (1997—2007年, ＞1100000次事件)。特别值得注意的是,我们选取了震级为M的主震后 12小时内在一个半径为0.02×100.50Mkm的圆内发生的震级较小的地震。最后,只有 2%的所选余震被记录在有震源机制解的地震活动目录中。因此,这里所分析的数据集不同于确定 b值和滑动角之间的关系所使用的地震活动目录[3]。

我们研究的只是来自特定震级范围的地震,即使是在余震序列的早期,这些地震的资料也非常齐全。有如下两个主要约束条件:(1)主震应该足够小,小到可产生较短尾波,这样才能对随后的余震进行可靠的探测;(2)余震应该足够大,以确保从余震序列的早期开始该特定震级水平上资料的完整性。这些震级范围由地震学和统计学约束确定 (见附 2:方法)。

附 2:方法

不用区分断裂模式,我们分别求出余震不同震级范围的 c值。由此,我们观测到幂率衰减之前的特征时间随着平均余震震级的增大而持续增加。在加州和日本的单个地震事件中也有类似相关性报道[5,9]。这些观测结果可以解释图 1a、b和图 2a、b中不同余震震级范围的平均 c值的变化。那么,平均 c值的绝对变化可以完全脱离 c值对滑动角的依赖关系。

在用最大似然法估算 c值时,优化程序取决于两个值,一个是 Tcomp,即目录完整的时间 (所有余震都有记录),另一个是 Tstart,即开始拟合的时间。很显然,如果 c＜Tcomp,估算的 c值就不会被认为是可靠的。如果c＜Tstart,c值误差就会很大,而且还会随着Tstart以指数方式增大。在 Tstart＜Tcomp＜c的情况下,c值可能会因为目录的不完整而被高估。因此,估算 c值的首选条件是Tcomp＜Tstart＜c。

对于日本,我们利用一个由人工挑选的目录[7]来确定不受余震缺失影响的 c值和Tcomp值。对于不同余震震级范围,我们证实c＞Tcomp。然后,采用 JMA目录中列出的地震事件,我们用增大的 Tstart值向 c值 (即由人工挑选出来的 c值)的合理估算靠拢。这种靠拢取决于记录事件的比例,也取决于它们的时间分布。对于 2.4≤MA≤3.4的余震,我们发现 JMA目录可被视为完整,我们所取的 Tstart值可以非常小。对于 1.8≤MA≤2.8的余震,我们推荐取 Tstart=5×10-4天,因为此时记录事件的比例从零迅速转变为 1。那么,还是这些地震事件,正断层型事件的 c值较大,逆冲型事件的较小,走滑型事件的居中(图2a)。

最重要的是,对于任何其完整性水平尚不确定的余震目录,可以使用优化程序的起始时间来估算 c值的稳定性,并证明缺失的余震不会对该参数的估算有任何干扰。本研究中,我们对南加州的情况进行了类似分析。我们发现在 Tstart=10-4天 (9 s)、Tstart=5 ×10-4天 (45 s)、Tstart=10 ×10-4天 (90 s)和 Tstart=14×10-4天 (2 m in)情况下,所有等级的滑动角的 c值估算都没有明显变化。由于MA≥1.8的余震目录在一次MM＜4.5的主震后 2m in时是完整的[10],所以 c值的这种稳定性清楚表明缺失的余震在这一特例中 (即 Tcomp＜Tstart＜c,包括不同滑动角)不会对该参数的估算产生影响。

[1]GutenbergB,R ichterC F.Frequency of earthquakes in California.Bu ll.Seismo l.Soc.Am.,1944,34:185-188

[2]Utsu T.A ftershocks and earthquake statistics.J.Fac.Sci.Hokkaido Univ.,1965,Ser.V II3:379-441

[3]Schorlemm erD, Wiem er S, WyssM.Variations in earthquake-size distribution across different stress regim es.Natu re,2005,437:539-542

[4]Narteau C,Shebalin P,HolschneiderM.Temporal lim itsof the power law aftershock decay rate.J.Geophys.Res.,2002,107,B2359,doi:10.1029/2002JB001868

[5]Shcherbakov R,TurcotteD L,Rund le JB.A generalized Omori’s law for earthquake aftershock decay.Geophys.Res.Lett.,2004,31,L11613,doi:10.1029/2004GL019808

[6]V idale J E,Peng Z,IshiiM.Anom alous aftershock decay rates in the firsthundred seconds revealed from the H i-netborehole data.Eos Trans.AGU 85(FallMeet.Supp l.),2004,abstr.S23C-07

[7]Peng Z G,V idale J E,IshiiM,et al.Seism icity rate imm ediately befo re and afterm ain shock rup tu re from high frequencywaveform s in Japan.J.Geophys.Res.,2007,112,B03306,doi:10.1029/2006JB004386

[8]Enescu B,Mori J,MiyasawaM.Quantifying early aftershock activity of the2004m id-N iigata Prefecture earthquake(MW6.6).J.Geophys.Res.112,B04310,do i:10.1029/2006JB004629(2007).

[9]Nanjo K Z,Enescu B,Shcherbakov R,etal.Decay of aftershock activity for Japanese earthquakes.J.Geophys.Res.,2007,112,B08309,doi:10.1029/2006JB004754

[10]Kilb D,Martynov V,Vernon F.A ftershock detection asa function of tim e:results from theANZA seismic network following the 31October2001ML5.1 Anza,Califo rnia,earthquake.Bull.Seismo l.Soc.Am.,2007,97:780-792

[11]Kagan Y Y.Short-term p roperties of earthquake catalogs andmodels of earthquake source.Bull.Seismo l.Soc.Am.,2004,94:1207-1228

[12]Lo lliB,GasperiniP.Com paring differen tmodelsof aftershock rate decay:the role of catalog incompleteness in the first tim es afterm ain shock.Tectonophysics,2006,423:43-59

[13]Peng ZG,V idale JE,Houston H.Anomalousearly aftershock decay rate of the 2004MW6.0 Parkfield,California,earthquake.Geophys.Res.Lett.,2006,33,doi:10.1029/2006GL026744

[14]Hauksson E.Crustal structure and seismicity distribution ad jacent to the Pacific and North America plate boundary in sou thern Califo rnia.J.Geophys.Res.,2000,105:13875-13903

[15]Gardner J,Knop off L.Is the sequence of earthquakes in southern Californiaw ith aftershocks removed Poissonian?Bull.Seismo l.Soc.Am.,1974,64:1363-1367

[16]Reasenberg P.Second-order moment of central California seismicity,1969—1982.J.Geophys.Res.,1985,90:5479-5495

[17]Schorlemmer D, Woessner J.Probability of detecting an earthquake.Bull.Seismo l.Soc.Am.,2008,98,2103-2117

[18] Wiemer S, WyssM.Minim um magnitude of completeness in earthquake catalogs:examples from Alaska,the western United States,and Japan.Bu ll.Seismol.Soc.Am.,2000,90:859-869

[19]Sibson R H.Frictional constraintson thrusts,wrench and no rm al fau lts.Natu re,1974,249:542-544

[20]Dieterich J.A constitutive law for rate of earthquake production and itsapp lication to earthquake clustering.J.Geophys.Res.,1994,99:2601-2618

[21]Scholz C.Microfractures,aftershocks,and seism icity.Bu ll.Seismo l.Soc.Am.,1968,58:1117-1130

[22]Shcherbakov R,TurcotteD L.A dam agemechanicsmodel foraftershocks.Pure App l.Geophys.,2004,161:2379-2391

[23]Ben-Zion Y,Lyakhovsky V.Analysis of aftershocks in a litho spheric model with seismogenic zone governed by dam age rheology.Geophys.J.Int.,2006,165:197-210

[24]Atkinson B K.Subcritical crack grow th in geological materials.J.Geophys.Res.,1984,89:4077-4114

[25]Amitrano D.Brittle-ductile transition and associated seismicity:experimental and numerical studies and relationship with the b-value.J.Geophys.Res.,2003,B108,2044,doi:10.1029/2001JB000680

[26]O’ConnellD R H,Ma S,A rchuleta R J.Influence of dip and velocity heterogeneity on reverse-and no rm alfau lting rup ture dynam ics and near-fault ground motions.Bu ll.Seismo l.Soc.Am.,2007,97:1970-1989

[27]Narteau C,Shebalin P,Ho lschneiderM.Loading rates in California inferred from aftershocks.Non lin.Process.Geophys.,2008,15:245-263

[28]Schorlemm erD, Wiem er S.Microseism icity data forecast rup ture area.Nature,2005,434:1086

[29] WellsD L,Coppersm ith K J.New empirical relationships among magnitude,rupture length,rupture width,rupture area and surface displacement.Bull.Seismol.Soc.Am.,1994,84:974-1002