akita moき乃はち:地震的余震是怎样形成的？

一个地震序列中最强的地震称为主震；主震前在同一震区发生的较小地震称为前震；主震后在同一震区陆续发生的较小地震称为余震。

余震其实是能量物质燃烧爆炸的继续，那是因为还有部分能量物质在接近地壳的部位上继续被引爆
其实余震就是地震发完后还有的震动叫做余震
只不过它的 能量是慢慢释放罢了

地震能量释放未完全，多次释放就是余震

就象你把一棵石头丢在水里一样,荡起的波纹和地震是一个道理

Omori (1894)研究了1891年发生在日本中部的Nobi地震后的余震月频次和半月频次随时间的衰减后发现，余震在单位时间内的发震频次可以表示为

, (1)

其中t是以主震发震时刻为起点的时间，K和h为常数。现在多用字母c来代替h。

Omori曾尝试用指数衰减的形式来拟合余震频次的衰减，结果不是太满意。

Utsu (1961)研究了多个地震的余震活动后发现，余震的频次衰减要比Omori公式所描述的快些，用

(2)

来拟合余震频次更为合理。上式称为修正Omori公式。

Utsu (1995)总结了根据全球200个余震序列的p值发现，p值的范围为0.6~2.5，中位值大约为1.1，奇怪的是似乎找不到p值和主震震级的相关性。

二、高阶余震序列和Omori公式的扩充

Omori公式发表后，也有一些文章指出Omori公式和修正Omori公式并不适合某些地震的余震活动，因为这些地震的余震活动并不随着时间而正常衰减。产生这些情况的原因多是余震序列中存在有强余震，而这些强余震又伴有次一级的余震活动。Utsu (1970)和Ogata (1983)用修正Omori公式的和式来表述这种含有次阶余震活动的序列的发生频次，即

， (3)

其中H(t)为单位阶跃函数，T为主震发震时间，T1、T2、…、Tn分别为第1、2、…n个强余震的发震时间，K、c、p、Ki、ci、pi为常数。可以看出此公式表示的余震序列是由主震引发的余震活动和强余震带来的高阶余震活动合成的。

三、ETAS模型的简单形式

Ogata (1988, 1989, 1994) 将自相似的思想引人了Omori公式，认为不仅仅是余震序列中较强的余震能产生高阶余震，而是余震序列中任何一个地震均能产生自己的高阶余震。他采用灾害函数，更确切点说是条件强度函数来表示这一模型、由条件强度的定义

(4)

给出模型的条件强度函数

5)

式中的μ为常数，代表背景发生率；i取遍所有的地震。上式中右边和式的每一项代表每一事件对地震发生率所做的贡献。很自然地，常数Ki应该依赖于第i个地震的震级。Ogata在模型中取Ki只和第i个地震震级有关，并且为指数形式

(6)

其中，M0为参考震级，一般也取为震级下限。这种取法的根据是Utsu和Seki (1957)的经验公式

(7)

即余震区面积与主震震级M的关系。由上式可以大致得出余震的总数大致和主震震级的指数函数也正比，即

(8)

于是(3.5)式最终为

(9)

Ogata (1989)称满足(9)式的地震活动为“标准的”(Standard)的地震活动，而模型的名称为“传染型余震序列模型”。式中的参数p表示余震的衰减率，a表示由不同震级在产生余震能力方面的差异。他通过研究发现，震群型或前震型的地震往往有较小的a，甚至于会得到p<1 (Ogata, 1988)。

这种模型的思想和Vere-Jones和Davis的触发模型(Vere-Jones and Davis, 1966)相似。触发模型的思想也是认为序列中任一地震均能触发新的地震。不同的是，触发模型没有将震级和产生余震的能力联系起来。有意思的是，Vere-Jones和Davis利用触发模型拟合新西兰的地震活动性后发现，反幂律衰减(即修正Omori公式)优于指数衰减。这种并非巧合的事实，即余震频次在时间上的反幂律衰减，可能是能量扩散的结果(Kagan and Knopoff, 1980)。

四、ETAS模型的空间形式

在Ogata将时间上的ETAS模型推广至空间形式之前， Musmeci和Vere-Jones (1992)就用一种时空的成丛模型来分析意大利的地震活动。这种模型的条件强度函数为

(10)

或

， (11)

其中A、a、C、sx、sy、Cx和Cy为常数。这两种函数分别对应为扩散型的丛和Cauchy型的丛。对于固定的点(x, y)，当t®¥时，以上两式中余震频次随时间的衰减规律分别为t-1e-ct和t-2e-ct。

Kagan (1991)和Rathbun (1993)也给出了不同的时空上的条件强度模型。

Ogata (1998)将空间ETAS模型写成以下的形式：

, (12)

式中 (13)

这里的 为震级为M的地震产生余震的期望数目， 为余震时间分布的概率密度函数， 为余震空间分布的概率密度函数， j(M)为震级分布的概率密度函数，一般采用Gutenberg-Richter关系。Ogata (1998)给出了三个例子：

1.

2.

3.

其中H(u)为单位阶跃函数，d、a和q为常数。Ogata (1998)也将此模型推广到各向异性的情形，即用 来代替(13)式中的 ，其中

，

g、s1和s2为常数表示余震区为椭圆。

五、 参数估计和模型选择

在实际应用中，给定一组的观测数据，如何确定模型的系数呢？在统计学中常常用最大似然法进行参数估计。对于具有条件强度函数的点过程而言，其似然函数的形式为

(14)

其中[0, T]为观测的事件段，N(T)表示[0, T]内发生事件的数目，[0, ¥) 为震级范围，S为研究区域。为了简便起见，常常将上式写成对数形式

(16)

对于描述同一自然现象或者拟合同一组数据，均可以采用多个统计模型。Akaike (1974)提出的AIC方法(Akaike信息准则)就是如何在这些模型中选出最佳模型，方法是：假定logL为模型的对数似然值，则

(3.60)

其中k为模型中采用的参数的个数，取具有最小AIC值的模型作为最佳模型。可以看出，AIC方法既考虑模型对现象的拟合程度的好坏，也对通过无限制增加参数个数以促进拟合程度的行为进行了惩罚。

六、 参数估计和模型选择

利用ETAS模型可以直接进行的实际应用包括：拟合地震活动性、模拟地震目录、预报地震等。下文将对这些应用逐一做简单介绍。

6.1 拟合地震活动性

Utsu等(1995)曾将日本地震区划分为16区域，用简单形式的ETAS模型对每个区域的浅震活动进行了拟合，并用AIC准则说明了ETAS模型要比严格触发模型(restricted trigger model) 更适合描述这些地区的地震活动性。Ogata (1998) 用空间ETAS模型拟合了Tohoku和本州岛中南部地区地震活动性发现，虽然余震区面积在尺度上遵从Ustu-Seki公式，但余震频次在空间上并不是呈指数衰减或衰减更快，而是呈反幂律衰减。他猜想余震区应当分为两个部分：近场和远场。其中近场就是传统意义上的余震区，即断层区；而远场则对应与“广义上的余震”，如余震活动的迁移性，或者由地震断层破裂后大地应力场调整引起的地震。

ETAS模型也可以用来研究局部地震活动性。Guo (1997) 用简单ETAS模型研究了34个地震及其余震序列后发现，ETAS模型中的p值(用pE表示)比Omori公式中的p值(用pO表示)要大或者两者相等。板内地震的a值要远小于板间地震的。他还研究了各个参数之间的相关性，如a和pE是正相关的；Gutenberg的b值和pE对于板间地震是正相关的，而对于板内地震是负相关的。

Zhuang (1998)用简单ETAS模型拟合了新西兰Cape Palliser发生在1990年的ML5.3地震及其余震前后各个时期的地震活动性，发现了此地震前后存在着不同的活动特征，称之为活动相，并假想把地震活动在时间上分为中间相 (interseismic)、转换相 (preseismic)、活跃相 (coseismic)和调整相 (postseismic)四种。

6.2 模拟地震目录

利用点过程模型可以进行人工合成地震目录，这一点在Zheng和Vere-Jones (1994)，Musmeci和Vere-Jones (1992)的文中都有阐述。Ogata (1998)也给出了用空间ETAS模型进行数值模拟的算法。这一方法是基于“瘦化”法 (thinning method)的算法，需要进行高维数值积分而耗掉大量的计算时间。这里给出的是Davis等人(1997)的算法。这一算法是利用ETAS模型具有良好的分支过程(Athrey and Ney, 1972)的结构特征，并不需要进行数值积分。根据ETAS模型的条件强度公式，可以归结出以下特征方程

(16)

其中由(13)式知k(M)为为震级为M的地震产生余震的期望数目，j(M)为震级的概率密度函数，一般采用Gutenberg-Richter关系，即

， (17)

而b=bln10，b为Gutenberg-Richter关系中的b值，v(M)为特征函数。可以得出特征值

, (18)

r即临界性参数：若r<1，这个过程是亚临界的，事件总数目随着时间的增长而线性增长；若r>1，这个过程是超临界的，事件总数目随着时间的增长呈指数性或爆炸性增长。r的另一个含义是所有事件中的任一个事件产生余震的期望数目。模拟算法如下：

设进行模拟的目标区域为A，B是一个包含A并且比A大的一个区域。

(1) 根据背景发生率，为简单起见，这里仅设其为一常数m，模拟出区域B内发生的地震总数目N。N是服从期望为mST的Poisson分布的随机变量，其中S为区域B的面积，[0, T]为所要模拟的时间段。

(2) 模拟这N个地震的发震时间，即N个服从[0, T]上的均匀分布的随机变量。

(3) 模拟这N个地震的震中位置，方法同步骤1可简单设其为B上的均匀分布。

(4) 利用(17)式模拟出N个地震的震级。

(5) 设这N个地震为第0代。可按如下步骤模拟出下一代地震的时间、地点和震级。

(6) 对于第0代中的每个地震，模拟出各个事件的余震数目，即服从期望为k(M)的Poisson分布的随机变量。

(7) 针对第0代中每一地震的每一个余震，首先模拟出发震时间，即服从反幂律分布或Pareto分布；然后根据余震位置分布的p.d.f模拟出其震中位置，最后根据余震的震级分布震级。称这样产生的地震为第1代地震。

(8) 针对第1代地震，重复与6和7类以的步骤，如此产生第2代、第3代，…，直至再没有后代产生或所有地震的发震时间均落在[0, T]的外面。

(9) 将上面得到的各代地震合并在一起，按时间排序后，挑选出A中的即为所求的地震目录。

这一方法和传统的点过程的模拟方法，如反演法和瘦化法是等价的。Ogata(1998)利用瘦化法和空间ETAS模型进行模拟，成功地再现了日本本州岛西南部的地震活动特征。庄建仓、Vere-Jones和Savage (1997)利用空间ETAS模型拟合了新西兰北部1966年至1990年的浅震活动后，用所得得参数模拟出了一个新的合成地震目录。对比合成目录和真实目录后发现，合成目录中出现了较多的前震，产生这一结果的原因可能是因为真实地震目录中的余震震级可能和主震震级有关造成的。尽管如此，仍说明ETAS模型对于研究前震仍是极其有用的。

6.3 地震预报

Vere-Jones (1998)曾提出了基于点过程模型的概率预报方案。这一方案的要求是这一点过程模型必须有显式的条件强度函数。预报方法如下：

假定t为当前时间，(t, t+Dt)为预测的时间区间。

(1) 根据时间t以前的观测结果，拟合观测结果得出模型中的参数。

(2) 根据求出的模型参数，模拟(t, t+Dt)的地震活动。

(3) 重复步骤(2)多次，如至少1000次。计算(t, t+Dt)内有地震发生模拟次数在全部模拟次数中所占的比率。这一比率可以看作是所求的发震概率。

这一方案有这样的优点，即随着时间的推移，可以将新发生的地震添加到观测结果中，重新估计模型参数 (称为学习过程)，然后用新的时间参数去预测下一个时间区间的发震概率。

Vere-Jones还给出了这一方法的评分准则，这里不再重复。

将这一方案应用到ETAS模型上的效果并不理想。Vere-Jones (1998)用拟合新西兰惠灵顿地区5年内的区域地震活动性得到的简单ETAS模型的参数模拟了一个人工合成地震目录，把上面的预测方案应用到预测时间区间为2天的情形上。结果表明，ETAS模型对大震的预报效果不太理想，但是对大震之后的余震预测效果很好。对大震预测效果不理想的原因可能是因为各个预测区间多开始于平静时刻，而非余震序列内。由于ETAS模型在无震区间预测较Poisson模型好，其预测的总体效果仍好于Poisson 模型。尽管如此，这一模型和这一方案在大震现场用来进行预测余震仍应当是非常有效的。

6.4 相对平静

Ogata (1986)的相对平静来源于这样的思想，主震发生之后，其衰减规律应按照Omori公式或ETAS模型的模式进行衰减，实际地震活动有时会偏离这一规律，即期望在某一时刻发生的地震没有发生，那么会在不久的将来发生一个强余震或一个由数目较多的余震组成的地震丛来补足这一“缺陷”。这一思想对预报震后强余震往往很有效，Ogata 曾成功地用根据这一思想发展起来的理论预报了兵库地震的一个强余震。

检验这种相对的平静是有一定困难的。因为各个地震的余震的发生时间并非等间隔分布。随着时间的推移，相邻余震之间的时间间隔越来越长，将缺震引起的相对平静和这种越来越长的时间间隔区分开来会有一定的难度。在这里Ogata考虑了残差变换，这种残差变换可以将一个具有条件强度函数的非Poisson点过程的事件序列变换成一个遵从单位发生率的Poisson点过程的事件序列。残差变换的定义为：

。 (19)

经过残差变换后，得到的新序列{ti}称为转换时间序列。任何一个用于检验Poisson序列的方法均可以用来进行检验相对平静是否存在。

上述方法可以应用到(9)式中的ETAS模型所定义的标准地震活动的平静的判别上。首先给出两个模型假设：

1. 地震活动在所研究的整个时间段[0, T]内遵从单一的ETAS模型，即各个参数在此期间无显著变化；

2. 地震活动在时间点T0Î[0, T]发生变化，表现在模型拟合上，ETAS模型的各个参数从时间段[0, T0]到[T0, T]时各个参数发生了显著变化；

以上两个模型对应的AIC分别为

其中，p(N)是当T0作为可调整的参数时的惩罚项，和事件的总数N有关，如果T0的选择是科学根据而非从数据本身得到的，则p(N)=0；

；

，

即lnL(q; S, T)为参数q对应的似然函数值。

AICII<AICI说明地震活动性在T0前后截然不同。

这种截然不同可能说明地震发生率有所变化。Ogata就是在满足AICII<AICI的条件下，检验经过残差变换后得到的转换时间序列的发生率是否有所改变得到的。具体做法是：用ETAS模型分别拟合两个时间段和总的时间段的地震活动性，如果AICII<AICI，则用拟合前一时间段数据所得的参数针对整个时间段计算转换时间序列，然后将转换时间序列的频次累积曲线画出来，如果在T0点后的积累速度明显下降，则可以认为有相对平静发生了。Ogata (1992)曾用此方法验证了多个大地震之前相对平静的存在，如1960年的智利MS7.0地震，1923年的日本东京MS8.3地震。

地壳释放能量的过程