2010年4 月14 日青海玉树发生M_S7.1 地震(震中位置33.218°N，96.594°E).在中长期时间尺度上，国际上通用的一些算法给出的结果(例如Kossobokov等发表的CN 算法的全球预测http://www.phys.ualberta.ca/mirrors/mitp/[2010-06-25])均未能给出这次地震前“强震发生概率增长"(TIP)的信息，因而目前时间相依的地震危险性估计的方法，似乎不适用于这一地区.本文试图通过考察该地震前的加速矩释放(AMR)现象，对这些算法在这一地区的“失效"给出一个可能的解释.这种回溯性的系统研究，尽管在地震危险性分析和中长期地震预测的意义上仍价值有限，对于反思和改进目前的地震危险性估计方法，却有不可取代的参考价值.

物理上对AMR 现象的兴趣，至少部分地来源于Sornette和其他一些物理学家^[1]，他们将AMR现象与地震孕育过程中的“临界性"联系在一起.迄今国内外学者已将AMR 分析方法应用于诸多震例^[1~6]，但正面和反面的震例均有.定量上，Varnes^[7]、Bufe和Varnes^[2]和Bufe等^[3]将AMR 分析表示成“破裂时间分析"(Time-to-failureanalysis)方程:

(1)

其中Ω 是地震活动性的量度，例如地震矩(或地震能量)、地震数目或Benioff 应变等(本文中使用Benioff应变)，A和B为常数，m为标度常数，t_f 是大地震发生的时间.

由(1)式可知，标度常数m表示了地震矩释放曲线的类型，当m<1 时，累积Benioff应变释放(CBS)曲线呈现“加速"行为；m>1 时为“减速"；而当m~1时CBS释放曲线呈现线性行为.利用(1)式可研究强震的发生时间和震级^[1~6]，以及未来强震的发生地点^[8].m值有时也称为“地震矩释放程度"，是本文的讨论重点.

以AMR 为切入点考察玉树地震，有双重目的.一方面，目前关于AMR，不同的作者给出的结果很不相同，且争议很大^{[9， 10]}.解决这一争议的重要方法之一是针对更多的震例进行详细研究.另一方面，AMR 是目前应用较多的在中长期时间尺度上估计地震危险性的方法之一.从AMR 入手研究玉树地震前的地震活动异常，有助于阐明为什么一些国际上通用的算法(这些算法主要是以地震活动性为“输入"信息的)不适用于玉树地震的情况.

玉树M_S7.1地震发生在巴颜喀拉块体南边界的玉树-甘孜断裂上，该断裂同时也是中国大陆II级活动地块---羌塘地块与巴颜喀拉地块的分界线.受印度板块向北推挤和青藏高原东向逃逸运动影响，羌塘地块与巴颜喀拉地块水平运动方向均为SE 方向，但由于羌塘地块相对水平运动较快，玉树-甘孜断裂表现为走向NW 的左旋走滑断裂，平均运动速率7.3mm/y.

本文选择巴颜喀拉地块及附近(31°~38°N，90°~102°E)作为研究区，该区历史上曾发生多次M≥7以上强震，其中2008 年的汶川地震也属于一些国际通用的中长期预测算法失效的情况^[11].研究使用了中国地震台网中心提供的(1970-01-01~2010-04-14)《中国地震月报目录》和《全国编目系统速报目录》.图 1 给出了研究区的活动构造和1970年以来M_L2.5以上地震的分布.

图 1

Fig. 1

图 1 玉树Ms7.1地震周边地区活动构造和1970年以来ML2.5以上地震分布，震前20年的地震用绿色实心圆表示 Fig. 1 Distribution of earthquakes larger than ML2. 5 since 1970 and the main tectonics in the study region.Green solid circles show the events occurred within 20 years before the Yushu Ms 7. 1 earthquake

为得到研究区地震目录的最小完整性震级M_c, 分别采用了定性和定量分析方法，考察M_c 随时间的变化.定性分析中采用了震级-序号法^[12]，即按地震发生时间的先后顺序排序，考察不同震级事件的数密度分布，来定性分析M_c, 其中地震数密度较大的位置对应的震级即是最小完整性震级.使用地震序号，而不是地震发生时间，主要是为避免“丛集事件"的影响.图 2a给出了由震级-序号法描述的M_c时序变化，结果显示1988 年前后M_c 大致分别在M_L2.5~3.0和M_L2.0~2.5之间.

图 2

Fig. 2

图 2 玉树Ms7.1地震周边地区1970年以来地震目录的完整性分析 (a)震级-序号法给出的地震密度分布；（b)Mc-Best法给出的最小完整性震级Mｃ的时间变化. Fig. 2 Catalogue completeness from 1970 in the surrounding areas of the Yushu Ms7. 1 earthquake (a) Seismc rate distribution by using Magnitude-Rank analysis； (b) Minmum completeness magnitude Mc by using Mc-Best method, varying with tme.

定量评估采用了Mc-Best方法，即根据“最大曲率"MAXC 方法和拟合度分别为90% 和95% 的GFT方法^[13]分别计算并遴选出各时间段的最佳结果.其中MAXC 方法将震级-频度曲线的一阶导数最大值对应的震级作为M_c, 而GFT 方法通过搜索实际和理论震级-频度分布的拟合程度确定M_c.由于同时采用90%和95% 拟合度的GFT 方法，这里分别称为GFT-95 和GFT-90.计算中使用固定为800个地震事件的窗口选取数据，并进行滑动计算.这里按照[GFT-95优于GFT-90优于MAXC]的优先级顺序，选择每次滑动的计算结果，选择后的M_c如图中标注为Mc-Best的黑色曲线所示.Mc-Best的结果显示，1970年以来研究区M_L2.5以上的地震是完整的，如图 2b所示.

为考察玉树M_S7.1 地震前的AMR，分别考虑了不同时间尺度T、以震中为圆心的不同圆形研究区半径R、截止震级M_c 对计算结果的影响，表示为m(T-R-M_c)的分布.将(1)式中的t_f 固定为实际发震时刻；T以1y步长，震前1~20y变化；R以10km步长，10~250km 多尺度半径选取圆形区域；M_c 以0.1个震级单位步长，M_L2.5~4.0 多尺度变化.与Brehm 和Braile^[4]在研究新马德里地震带的AMR现象时的做法相同，在分析中删除与主震发生在同一时期、较近范围内、具有相似或更高震级的“干扰"事件，“干扰"事件的震级一般定义为与主震差1 个震级单位.

由此得到m(T，R，M_c)的分布如图 3a所示.作为示例，图 3(b, c)分别给出了对应T=18y、R=50km 和M_c=M_L2.5的情况下，去除“干扰"事件前后的Benioff应变曲线及其拟合的情况.模型计算中需考虑避免“过度拟合"的问题，而BIC 增益能较好地平衡拟合残差与模型自由度之间的关系.这里采用了Jiang 和Wu^[11]提出的AMR-BIC 判据，其中BIC 的计算采用了Seher和Main^[14]给出的简化的BIC 表达式.只有当非线性拟合(1)和线性拟合的ΔBIC≥0 时，才能说采用公式(1)而不是简单的线性拟合是合理的.与Bowman等^[15]早期建议的非线性拟合与线性拟合的均方根残差比(C值)和经验地给定的阈值(C小于0.7)相比，Jiang 和Wu^[11]的AMR-BIC 判据可以自然地避免拟合效果判定的主观性.图 3(b, c)中也同时标出了ΔBIC 的结果.为保证m值的可靠性，图上还删除了地震序列中数据点数少于5 个、或者ΔBIC<0、或者拟合误差Δm>0.5的计算结果，在图 3a中用空白色块表示.

图 3

Fig. 3

图 3 玉树Ms7.1地震前矩加速释放分析 (a)在（T-R-Mc)三维空间中m值的分布；（b)T =18y、R= 50km且未去除“干扰事件”的Benioff应变曲线；（C) T =18y、R= 50km且去除 “干扰事件”的Benioff应变曲线；（d)T =3y、R= 110 km且去除“干扰事件”的Benioff应变曲线.图中空白色块表示地震数目不足、ΔBIC＜0,或者m值拟合误差过大的情况，见正文. Fig. 3 AMR before the Yushu Ms7. 1 earthquake (a) The m value distribution in the (T-R-Mc) space； (b) Cumulative Moment Release (CBS) curve with T=18y,R=50 km. The crossshows the 4interfering event? with magnitude larger than Ml6. 1； (c) CBS curve with T=18y, R = 50 km, with 4 interfering events? removed； (c) CBS curve with T=3y, R=110 km, with interfering events? removed. Blanks pixels denote the cases that the number of earthquakes within the spatial-temporal range is less than 5, or ΔBIC＜0 or RMS residuals exceeding the limitation. See text for details.

由图 3a可见，截止震级M_c 更多地通过地震数目的多少影响m(T，R，M_c)的分布，而T和R的选取对m值影响较大.其中两个区域存在较稳定的AMR现象，分别对应T=10~20y、R=50~120km, 以及T=1~4y、R=120~250km.图 3d给出了与图 3c分属不同AMR 区域的Benioff应变曲线示例，图中T=3y、R=110km、M_c= M_L2.5.由图 3，回溯性地看，2010年玉树M_S7.1 地震前，AMR 的存在是客观的.

玉树M_S7.1主震前两小时，震中附近发生一次M_S4.6地震(33.113°N，96.597°E).为阐明这次地震对AMR 曲线的影响，将这一M_S4.6 地震去除后，重新计算m(T，R，M_c)的分布，如图 4所示.由图可见，T=1~4y、R=120~250km 对应的AMR区域几乎消失，这表明相应时空尺度内的AMR，基本上是这个M_S4.6 地震的作用.换句话说，如果这一前震不发生，则在T=1~4y、R=120~250km 的时空尺度上是看不到AMR 的.同时去掉这一地震后，另一时空尺度(T=10~20y、R=50~120km)的AMR 现象仍较为稳定.因此，较短时间尺度上看到的AMR，更多地反映了前震的作用.

图 4

Fig. 4

图 4 删除2010年4月14日Ms4.6地震后得到的三维空间m值分布.图注见图 3 Fig. 4 The m value distribution in the (T-R-Mc) space,with the 14 April,2010, Ms4. 6 earthquake removed. For captions see Figure 3

已有结果通过诸多震例研究，给出了“目标地震"的震级与临界时间尺度T和空间尺度R之间可能存在的定标率，以期用于向前预测.Bowman等^[15]通过对美国加州地区的AMR 震例分析，得到

(2)

式中R_C 为出现AMR 现象最明显的空间临界半径，M为相应的地震震级，c，d是常数，其中c=0.44；Brehm 和Braile^[4]对新马德里地震带进行了类似研究，得到c=0.75；Jaumé和Sykes^[16]综合前两者的研究，得出c=0.36.Zller等^[17]选用美国加州1952年以来M≥6.5 的所有地震进行系统分析，得到lgR_C ∝0.7M.Papazachos和Papazachos^[18]对希腊爱琴海地区的研究认为c=0.42，d= -0.68.对于AMR 的持续时间T与主震震级M，一般认为无明显的定标率存在^[4].

为考察玉树M_S7.1 地震震前出现AMR 现象的时、空尺度与前人结果的关系，这里根据蒋长胜和吴忠良^[19]收集整理的M6以上地震的AMR 震例及其参数^{[2~4， 15， 18， 20~23]}给出图 5 所示的定标率关系.应该指出，不同的作者选择“最佳"时空临界尺度的方法不尽相同.而与前人研究不同的是，本研究未刻意选取“最佳"的临界时、空尺度，而是将m(T，R，M_c)分布上出现稳定AMR现象的时空范围T=10~20y和R=50~120km 作为一个“范围"标注在图上.作为比较，图中还给出了采用类似的m(T-R-M_c)分布方法得到的2008 年于田M_S7.3 地震^[8]、2008年汶川M_S8.0地震^[11]的结果.其中，由于汶川M_S8.0地震的破裂尺度较大，且周边的强震活动较多，在考察其AMR 现象的时-空关系时，作为Brehm 和Braile^[4]去除“干扰事件"的方法的一个发展，使用了“蚀法"(eclipsemethod)删除周边主要活动断裂上的强地震事件.

图 5

Fig. 5

图 5 AMR研究中主震震级与时、空参数选取的关系 (a)主震震级与时间尺度T的关系；（b)主震震级与空间半径R的关系.图中橙色、紫色和红色线段分别对应2008年于田Ms7. 3地震、 2008年汶川Ms8. 0地震和2010年玉树Ms7. 1地震的AMR时、空范围，黑色直线为除上述三个震例外的线性拟合，灰色虚线为拟合的95%置信区间. Fig. 5 Relation between the magnitude of the 4 target earthquake7 (M) and the critical spatio-temporal scales (T and R) in the AMR study (a) M-T relationship; (b) M-R relationship. The orange, purple and red vertical lines represent the 2008 Yutian Ms7. 3 earthquake, the 2008 Wenchuan Ms8. 0 earthquake, and the 2010 Yushu Ms7. 1 earthquake, respectively. Black line shows the linear fit of all the cases except those three earthquakes, and the grey dashed lines show the 95% confidence interval.

由图 5a可见，与前人认识基本一致的是，采用m(T-R-M_c)分布方法得到的3次7级以上强震震级M与AMR 持续时间之间T不存在明显的定标关系.在图 5b中，3次强震对应的R则系统小于相同M下的其他震例的结果，与以往主要是在板块边缘地区得到的定标关系不同.这种定标率的不同，或许正是一些国际上普遍采用的算法未能给出玉树地震前“强震发生概率增长"信息的原因之一.另一个原因是，国际上普遍采用的算法，往往使用了比较大的截止震级^[11].也许仅3个地震尚不足以修订原来的定标率.但在今后的巴颜喀拉块体及其附近的地震危险性分析中，注意到这一定标率的差别，可能是重要的.

AMR 现象与强震的孕育发生是否物理上相关，一个需要考虑的方面是在假设空间各点均可能是(候选的)破裂成核点的情况下，考察不同时间尺度下m值的空间分布与“目标"地震的破裂成核点的对应关系^[8].为此，将研究区按照0.1°×0.1°的尺度进行网格化，假设空间各点均为可能的破裂成核点，并将式(1)中的tf 固定为玉树地震的发震时刻，空间尺度R固定为50km, 截止震级固定为M_c=M_L2.5，分别考察T=6y、10y、14y和18y情况下m值的空间分布，如图 6 所示.由图可见，在各时间尺度下，玉树M_S7.1 地震震中附近均出现稳定的AMR“热点".但值得注意的是，研究区其他地点同时也有大量的AMR“热点"存在，这就使得利用m值空间分布寻找可能的地震破裂成核点存在较大的困难.这个问题，对于其他的算法，例如PI算法^[24]，也是存在的.实际上，AMR 现象用于实际的地震预测预报的困难，除了计算所用时、空参数难以确定外，其与强震的发生不存在时、空上简单的一一对应关系^[10]，也是一个值得重视的方面.

从图 6，沿时间轴“反向"考察(即按照d-c-b-a的次序看AMR“热点"的变化)，可以看到随着时间的“逼近"，“热点"的范围趋向于向震中“聚集".而区分由于数据点数的变化而引起的“自然的"“聚集"和“真正的"“聚集"，就需要进一步的定量分析.这里，我们试图借鉴“迁移图像"的做法，考察是否在玉树地震前还能看到更多的东西.

图 6

Fig. 6

图 6 玉树Ms7. 1地震前多时间尺度的矩加速释放扫描 (a) T =6y；b) T =10y;(c) T =14y;(d) T =18y.图中空心六角星标出了玉树Ms7.1地震震中，十字标出计算中排除掉的“干扰事件”计算中R固定为50 km. Fig. 6 Seismic moment release mapping before the 2010 Yushu Ms7. 1 earthquake (a) T= 6y ； (b) T= 10y; (c) T=14y; (d) T= 18y. The epicenter of the italic>Ms7. 1 earthquake is marked by the hexagonal star； cross shows the removed ‘interfering events’;R is fixed as 50 km.

作为一种可能的地震前兆，Mogi^[25]认为强震前可出现地震活动异常向未来震中迁移的“迁移图像"(migrationpattern)，并将其用于识别大尺度的与地震孕育相关的应力积累和成核过程.近年来，这一概念重新引起关注.为定量描述强震前可能出现的“迁移"过程，Wu等^[26]通过定义“误差距离"ε 和“覆盖率"f来计算图像信息学(PI)算法给出的PI“热点"向未来震中迁移的程度，并对2006 年台湾屏东M_L6.4和M_L6.7双震进行了有趣的震例研究.其中ε为最接近震中的PI“热点"到震中的距离；覆盖率f=AH/A，AH为某一域值下的PI“热点"的面积，A为研究区的总面积.Wu等^[26]按照“两步法"来度量“迁移":(1)计算平均误差距离〈ε〉随f增加的衰减变化；(2)将〈ε〉-f的曲线积分后获得震前时间尺度T对应的ε_area(T)，如果ε_area(T)随T缩短并逐渐逼近发震时刻而衰减，即认为出现震前的“迁移".

本研究借鉴了“迁移图像"的做法，但所给出的并不是真正意义上的“迁移图像"---因为这里涉及的是累积AMR 的不同的时间尺度，而“迁移图像"给出的是线性时间轴上的空间分布的演化，因而，本文的方法所给出的结果可能在表观上更“清楚"一些，但并不是原来意义上的“迁移图像".为示区别，我们权将此图像称为“累计迁移图像"(MCP-MigratingCumulativePattern).在度量AMR“热点"的“迁移图像"时，在Wu等^[26]的做法基础上，在定义误差距离ε时，为降低个别“热点"对结果影响较大的可能，考虑了配以权重的所有AMR“热点"对空间参考点的平均作用.由于空间各“热点"到x_i点的距离r差异较大，这里在考虑空间全部“热点"对xi的平均距离时，采用高斯核函数对r进行平滑和归一，并直接将其作为“平滑距离":

(3)

其中σ 为用于平滑的参考距离，这里人为设定为50km.常数c被设定为0，即意味着对x_i处的“热点"不进行平滑.按照一般理解，“热点"到x_i的距离r越大，对该点的作用应该越小，所对应的误差距离ε则应该越大.这里定义空间第j个“热点"对xi的误差距离为ε_j=1-k_j.因此，在特定阈值的覆盖率f下的平均误差距离〈ε〉_f 可表示为:

(4)

上式中rj为第j个“热点"到xi点的距离，nf 为覆盖率f对应的“热点"数目.经过高斯核函数进行平滑和归一处理后，这里的〈ε〉成为无量纲变量.图 7给出了震中位置处的平均误差距离〈ε〉与f的关系曲线示例.

图 7

Fig. 7

图 7 玉树Ms7.1地震震中位置的“累积迁移”计算示例 (a)覆盖率f与平均距离误差(ε)曲线，图中不同颜色对应不同的震前时间尺度T;（b)积分后误差距离εarea 与震前时间尺度T的关系曲线，图中虚线为线性拟合，相应的斜率为一0.0066. Fig. 7 Migrating cumulative pattern priroi to the 2010 Yushu Ms 7. 1 earthquake (a) The coverage fraction f versus the average error distance (ε) , different colors indicate different time scales T；(b) The integrated error distance εarea versus T； dashed line shows the linear fit, corresponding to slope 一0. 0066.

与Wu等^[26]的做法类似，同样先对〈ε〉与f的关系曲线进行积分，得到时间尺度T下的“积分后的误差距离"ε_area(T).不断缩短T并使其逐渐逼近地震的实际发震时刻t_f, 对得到的一系列ε_area(T)进行线性拟合，其斜率表示了x_i点处的图像迁移程度.由于“迁移图像"描述的是异常点向参考点x_i的汇聚程度，这时ε_area(T)的斜率小于0.

因此在这里，我们采用空间网格化，假设任意格点/空间参考点x_i都是可能的破裂成核点，给出“迁移图像"程度的空间分布，进而考察与强震发生地点的关系.同时对ε_area(T)-T进行线性回归，将其斜率进行空间归一，用概率形式表示“累积图像迁移"的程度.对空间所有网格点分别计算ε_area(T)的斜率，对于斜率≥0的格点则统一设为0.将这种空间汇聚程度进行归一化处理，直接作为“迁移"意义上的地震危险性概率prob.对于玉树M_S7.1地震，分别计算T=6，7…18y对应下的ε_area(T)，并给出空间归一化的斜率或“累积图像迁移"意义上的地震危险性概率prob，结果如图 8所示.

图 8

Fig. 8

图 8 基于AMR图像迁移的地震危险性概率分布 图中空心六角星标出了 2010年4月14日玉树Ms7. 1地震的 震中，橙色曲线为研究区活动断裂.图中还用灰色、蓝色和紫色 曲线分别表示σ=20 km、σ=80 km和σ=110 km情况下的fprob= 0.7等值线. Fig. 8 Spatial distribution of seismic risk probability based on the AMR migrating cumulative pattern The hexagonal star marks the epicenter of the Ms7. 1 earthquake, orange curves show the active faults in the study region, and the gray, blue and purple curves indicate the contour lines of prob =0. 7, with、σ=20 km,、σ=80 km and、σ=110 km, respectively.

由图 8可见，玉树M_S7.1地震之前震中附近出现明显的prob高值分布.通过“迁移图像"的考察，研究区其他地方的AMR“热点"不具有随时间“聚集"于震中附近的性质.

每次地震的发生都会给地震危险性分析提出新的问题和新的研究线索，对地震孕育发生的认识的过程，也就是这种不断地积累经验、总结教训的过程.本文回溯性地考察了2010 年玉树M_S7.1 地震前的AMR 现象.在已知发震时刻和震中位置情况下，对(T-R-M_c)三维空间中矩释放程度m值的分布进行分析，得到在时间尺度T=10~20y 和空间尺度R=50~120km 范围内存在稳定的AMR；与前人研究AMR 获得的震级与AMR 时空尺度的定标率对比，利用m(T-R-M_c)分析获得的巴颜喀拉块体周边的3次地震(于田、汶川、玉树)给出的R值系统偏小；多时间尺度下的矩释放程度m值给出的AMR“热点"尽管在震中位置有分布，但空间上AMR“热点"分布较多，难以惟一地识别出未来地震的震中位置；借鉴“迁移图像"的方法、但在累积分布上考察AMR“热点"的演化，所给出的描述或可反映玉树地震的“逼近".因此，在通过一些地震得到AMR 的定标率，并同时考虑“迁移图像"后，巴颜喀拉块体周边十年尺度的时间相依的地震危险性，或可得到比较好的估计.

在“迁移图像"计算中将平滑参考距离σ 设定为50km 带有一定的主观性和任意性.事实上，统计预测模型中参数的合理与最优设定，通常是比较困难的事情，如果统计模型的预测效果对模型参数的设定过于敏感，则说明这种“前兆现象"不具稳定性进而无法用于实际的“向前"预测^[27].因此，我们还分别考察了σ=20km、σ=80km 和σ=110km 时地震危险性概率prob=0.7等值线的分布，结果如图 8所示.由图可见，平滑参考距离σ的选取对prob的分布可产生一定影响，但结果仍较为稳定.

必须提醒，地震现象很复杂，从玉树地震得到的经验教训，不一定适用于其他地震.但将由此得到的知识贡献于面向未来的时间相依的地震危险性研究，无论如何还是值得开展的探索.

中国地震台网中心提供了《中国地震月报目录》和《全国编目系统速报目录》，中央大学(台北)陈建志教授、数理统计研究所(东京)庄建仓博士、东京大学KazuyoshiNanjo博士与作者的讨论对本文启发很大，闻学泽研究员、张国民研究员和傅征祥研究员对本文提出有益的指导意见，两位评审专家的意见对稿件的修改帮助很大，相关工作得到中国地震局监测预报司“中国大陆未来10年M≥7.0地震危险区预测"工作专项的支持，作者谨表谢意.