A.Liu M.Yamada

引言

在主震之后的地震活跃期间，可在不同的位置几乎同时发生多个地震。在这种情况下，由地面传感器测到的地震波包含不止一个震源的混合信号。如果探测算法假设仅为一个地震，那么估计出的地震参数（比如位置和震级）就不会精确。而这些不精确的估计也会导致虚报，这种情况在大地震后经常出现。

2011年发生在太平洋沿岸近海的东北地震（此后称东北地震），摧毁了本州岛东北部的大部分区域。在P波达到后的8s，即发震时刻后的31s，对东北地区的公众就发布了地震预警（EEW。Hoshiba and Iwakiri，2011；Hoshiba et al，2011；Sagiyaet al，2011）。在这个区域中没有盲区，也就是说，在S波达到之前，所有位置都接收到了预警，这是因为该地震离海岸相当远。

该主震之后有大量的余震发生，使得在主震后的第一个月内发布了70次预警 ［Japan Meteorological Agency（JMA），2011］。这些预警63%都有重大的误差，其中估计的地震烈度要比观测烈度至少大了两度。相比而言，在东北地震之前，仅有29%的预警含有这种误差。事后分析表明，这些误差有73%可归因于即不能根据由短时间分开的同一震源也不能根据空间遥远的成因来区分多个同时发生的地震（JMA，2011）。这些虚报的主要原因之一是当前的方法主要使用P波到时来估计震源。

本文中，我们提出一种新方法来检测和区分当前日本气象厅（JMA）系统框架内的多个同时发生的地震。我们引入近似的贝叶斯方法来估计同时发生余震的位置、震级和发震时刻。不同于当前的日本气象厅系统，这种方法对发生时间靠近的地震产生多组估计结果。根据几个典型例子研究的实验结果说明，这种方法可以成功地检测并估计多个同时发生的地震的参数。

表1 本文研究的地震信息汇总

1 数据和处理

本文使用的数据是由日本气象厅地震台站在东北地震期间及其后观测的强震数据。我们把这些记录按如下简述分成三段来对新分类法进行评估。对每一段，都将含有日本气象厅地震预警的数值与表1所示日本气象厅地震目录中出现的数值对比。

1.1 数据集

例子1：2011年3月15日，1∶36∶00～1∶38∶00（两个小地震）。根据检测到P波之后21s估计出日本气象厅震级为5.9而向公众发布了预警（见数据与来源一节）。然而观测的最大地震烈度只为日本气象厅地震烈度表的2度。如表1所示，在15s以内至少发生了2个相距200km左右、震级分别为2.5和3.3的地震事件。由于第二个地震开始的时间与第一个地震波至的时间非常接近，地震预警系统将这两个独立的事件当成了一个地震事件，从而导致高估了震级。

例子2：2011年3月20日，14∶19∶00～14∶21∶00（两个小地震）。日本气象厅地震预警系统在检测到初至P波之后6.6s估计出震级为7.6，并向公众发布了预警（见数据与来源一节）。然而，观测到的最大地震烈度仅为3度。同样，如表1所示，这一高估或许可归因于将在5s以内发生的相距约150km的2个小地震（日本气象厅震级分别为3.0和4.7）错误地当成了一个大地震，因为第二个地震发生的时间与第一个地震波至的时间非常接近。

例子3：2011年3月11日，14∶46∶00～14∶49∶00（东北地震）。为了说明这个方法也能处理单个地震事件的类别，我们还对东北地震（MW9.0）进行了分析。在P波初至后约8s，即发震时刻后的31s，对东北地区的公众发布了预警。

1.2 处理过程

本文使用的是来自200个地震台站、采样率为100Hz的三分量加速度数据。首先，将加速度数据转换为地震解析码格式，并用100为因子进行抽取，将抽样频率减少到1Hz。抽取的过程并不是必须的，但可用来减少计算时间。然后，将抽取的加速度的每个分量k（t）转换为位移A（t）。这个转换通过使用递归数字滤波器和机械地震仪的频率响应对k（t）进行两次积分完成。

式中，函数增益因子g0以及滤波常数h0，h1，h2和h3与抽样频率﹑阻尼常数和地震仪的固有周期有关。对于100Hz抽样率﹑0.55阻尼常数和6s固有周期的日本气象厅地震仪，这些值为：

下面的区分方法使用了三分量位移A（t）的矢量和以及加速度k（t）的垂直分量。采集用k（t）的短期平均／长期平均完成，短期窗为1s，长期窗为10s。该方法也计算预期的P波和S波到时（tP和tS）以确定台站是否观测到了P波、S波，或是都未观测到。这些到时用日本气象厅的一维分层速度结构计算（Ueno et al，2002）。

2 贝叶斯方法

多重地震事件的连续参数估计问题可用公式表示为贝叶斯推理问题。令θ为表征事件的参数矢量，Θ为用θ参数化的一组事件，从而Θ＝ ｛φ， ｛θ1｝，… ｛θ1，θ2，…｝｝。设z1：t为所有台站直到当前时间t的整个观测历史，后验概率P（Θt｜z1：t）则揭示在给定证据和先验信息的情况下在时间t当前发生事件的信息分布。

式中，P（zt｜Θt）代表似然函数，通常表示为L，L（zt｜Θt）＝P（zt｜Θt）。而 P（Θt｜z1：t－1）则是在时间t的预先更新的概率，

而P（Θ0｜z0）≡P（Θ0）是Θ的先验分布。

2.1 质点滤波

通常情况下，公式（3）并没有闭型解，而有几个近似后验分布的次优解（Arulampalamet al，2002），其中之一是网格搜索。网格搜索，虽然实现起来简单，但还是遭遇一些问题。首先，当参数是连续的、且没有严格限制时，这个方法由于只能搜索到有限的网格，因此并不能覆盖所有的参数空间。其次，网格大小是预先设定的，因此，为了获得理想分辨率的良好覆盖，这个方法在实现过程中需要大量的网格。

另一个解法是质点滤波（PF），即序贯蒙特卡罗方法，用一组加权质点近似后验分布（Doucet et al，2001）。当质点的数目趋于无穷时，由质点滤波得到的解也趋于最优解。关于质点滤波及其变种的文献很多（Liu and Chen，1998；Doucet et al，2001；Arulampalamet al，2002）。这个方法的基本过程简述如下做参考。

2.1.1 抽样

在每次迭代开始时，每个质点的值都根据重要的密度函数得到，这个密度函数为：

式中的～表示抽样Θit根据分布q（·）得到。

2.1.2 权重更新

质点滤波以一群加权质点近似后验：

式中wit是质点i在时间t的权重。所有权重的和归一化为1。

当新证据zt进入并且在每次更新后都重新归一化时，所有质点的权重都更新：

式中q（·）表示在抽样步骤中出现的同一重要密度。为了简化计算过程，q（·）常被当作P（｜）的替换。由于右边的项抵消，因此新的权重与似然函数L（zt｜Θit）呈正比关系。

2.1.3 重抽样

因为后验是以离散的质点近似的，在权重集中在一小部分质点时经过几次迭代更新后，该方程组就会遭遇抽样简并。权重方差的减少可确定简并度，可用近似（Arulampalamet al，2002）：

每一次迭代过程都包含一次采样和一次权重更新。仅当N︵eff减小到某一阈值时，才出现重采样。

2.2 模型

接下来，我们讨论基于质点滤波实时参数估计方程组在多重地震事件实际应用的细节。我们需要估计的参数为θ＝ ［x，y，D，M，t0］，其中，x指经度，y指纬度，D指深度，单位为km，M为日本气象厅震级，t0为发震时刻。完整的伪代码（图A1中的算法1）见附录。

2.2.1 先验分布

先验P（θ）决定质点的初始化状态。好的先验包含有地理信息，例如附近断层线相对于第一个被触发台站的位置，以及在断层线上产生的地震最常见震级。这些信息都可以根据每个台站的历史地震目录得到，并在对质点滤波进行初始化时实时使用。如果先验信息缺失，则可以使用平坦先验来进行替代。先验分布的选择会影响估计的质量和收敛速度。大覆盖的先验分布可能会使初始估计不稳定，这是因为几乎没有证据存在。小覆盖的先验分布则可能导致收敛速度慢或收敛错误（即在错误值上收敛）。这就需要根据经验来进行权衡。在本文中，我们使用统一的平坦先验，即位置为±100km，深度为±10km，震级为±1级，事件发震时刻为±10s。

2.2.2 似然函数

用质点滤波方法进行参数估计时，算法的性能在很大程度上取决于似然函数的设计。除了从当前日本气象厅方法使用的触发台站的到时和测量的振幅以外，我们的似然函数还使用了非触发台站上的相关信息，这是因为它们同样含有地震事件的重要信息。

在本文中，我们使用由日本气象厅研究出的用于震级估计的衰减关系。该衰减关系表述如下（Hoshiba and Ozaki，2013）。设Amax为地震事件发生后用地震仪测得的最大位移。地震的P波和S波震级MP、MS可以表示为从台站到震源的线性距离（R）、震源深度（D）、P波的最大位移（）或整个持续时间的最大位移）的函数：

图1 模型一节所用参数之间的关系说明。（a）震源和地震台站；（b）振幅和到时。tP和tS表示自地震在t0时刻发生后P波和S波的到时，tP≤tS

参数之间的关系在图1说明。这些公式是专门针对日本的地质情况定制的（见数据与来源一节）。由于位移的散射较小，因此P波和S波震级以最大位移Amax的形式表示，而不是最大加速度或最大速度。

假设如公式（10）和（11）且位移遵从对数正态分布A～lnN（μ，σ2），我们可以得到单一台站的如下似然函数：

此式中的Aexp是Amax的预期值，σ是位移测量值的标准差。依赖于观测P波和S波的地震台站，预期的最大位移及其标准差是不同的。为了方便起见，通过重新组织公式（10）和（11），我们可以计算下面三种情况的Aexp和σ。

注意公式（12）是基于振幅的方法，不同于基于标准到时的方法。采用这一方法的主要原因是发现，无震动信息是区分在时空上发生相近的多个地震事件的关键。这将在后面的讨论一节中进一步阐述。

·没观测到任何地震波：

·观测到P波：

·观测到S波：

Anoise和σnoise是在位移测量中由近期环境噪声带来的噪声，可以通过移动窗对每个台站单独计算出来。σP和σS可以根据历史地震数据先计算出来。对于单个台站，决定用哪个Aexp公式来计算取决于P波、S波是否到达或均未到达该台站。给定台站到震源的相对位置（x，y，D），P波和S波的预期走时（tP和tS）就可用射线理论计算得到。比较到时tP、tS，绝对当前时间t和绝对事件发震时刻t0，就可直接估计出一个台站使用哪一种Aexp计算。图2给出了这些设计想法的概要说明。

图2 单台站似然函数的设计说明概览。给定震源估计，由某一台站进行的预期观测取决于其是否可以观测到P波、S波或两者都观测不到

似然函数的这种设计依据的是P波或S波震动期间由地震仪测得的最大位移Amax。然而，地震仪或许不能在波至之后迅速观测到最大位移。在这种情况下，利用这种似然函数得到的初始估计值很有可能是不正确的。这时，我们可以用一个简单的延迟函数g（·）来近似最大位移被观测到之前的瞬时位移，

式中t和t0分别是绝对当前时间和绝对事件发震时刻。tP是预期P波走时。g（t）的形式是左移的S形函数。

文中，音译并非先生所青睐，先生认为主张音译的人,其实主张的是以音译为主,义译为辅。他还逐一批驳了玄奘的五不翻原则，大概随着时代的变革发展，这些翻译原则也不再那么绝对了。当提及章士钊先生所说的音译的两大好处不滥和持久时，先生坦然承认章先生说得足够透彻，但随即指出这两大好处也敌不过其坏处，而且义译也有自己的好处。此外音译也绝非不能用,以下两种情况可用音译。一为“所重在音的”，如人名和大部分地名。二为“意义暧昧的”，意义不够明确的。看惯了先生情感细腻的散文，第一次读到先生如此犀利、环环相扣的议论，难免讶异，却也欣喜，先生的文风谦逊、不卑不亢，说理非但不枯燥还让人忍俊不禁。

在每个时间步长应用似然L（·｜·）来更新每个质点的权重。假设每个台站都独立观测，所有台站的观测集族为z，则整个似然函数为：

式中，n指的是台站数。独立假设是我们简化的，因为邻近台站可能会有相关的观测资料。

2.3 广义质点滤波

质点依据参数的先验分布初始化。因为我们是用有界和离散的空间近似无界连续的五维空间，因此我们必须要确保质点有足够的覆盖范围和界内有所需要的质点数目。这在地震应用中尤其重要，因为可取的参数数目和值的范围都很大。以有限质点数目确保质点多样性的一个方法是采用广义质点滤波（RPF）方法（Arulampalamet al，2002）。

广义质点滤波方法与质点滤波方法的不同仅仅在于重采样阶段。常规质点滤波方法是根据后验密度P（·｜z）（公式6）的离散近似进行采样（Musso et al，2001），而广义质点滤波方法是从连续近似中进行采样，

图3 广义质点滤波中一些常用的平滑核。每个核函数的积分均为1，以确保所得到的密度仍然是概率密度函数

式中，nx为参数空间的维数，Cnx为单位球体的体积。图3给出了文献中一些常用的核。

带宽矢量h可以通过质点群的方差依比例进行设定，这主要是通过对协方差矩阵进行Cholesky分解计算得到（Bickel and Levina，2008）。

2.4 适用于多个同时发生地震的近似方法

当贝叶斯推断问题的准确结果很难得到时，我们可以使用质点滤波方法来进行解决；然而，为了使估计结果接近最优解，所需质点的数目必须与事件数目呈指数式增长关系。

所幸根据历史记录我们可知，当n较大时（n＞3），在60s的时间窗内n个地震事件同时发生的概率极小。因此将这一信息纳入到先验分布中，可以显著地减小状态空间的大小。然而即便这样，对于需要高效实时的计算而言，其状态空间还是太大了。例如，假设地震事件可由5个参数组成的矢量θ来表述，即θ＝ ［xyDMt0］T，其中， ［xyD］T是指经度、纬度和深度，M是指地震的震级，t0是指地震的发震时刻。如果仅有3个地震事件同时发生，那么需要搜索的就会是5×3＝15维的空间。

图4 2011年3月15日1时36分7秒之后2 000个质点分布的可视化图：（a）1s，（b）2s，（c）14s，（d）17s。时间对应于第一次检测到P波后的秒数。官方给出的这两次地震的震中位置用星号表示并在（d）中标出以供参考。本图的彩色版只适用于电子版本

作为首次近似，这个试探法是将所有可能的地震事件初始化为离散的质点滤波，即pf1（θ1），pf2（θ2），…，而不是将其初始化为一 整 个 质 点 滤 波：pf（Θ ＝ ｛θ1，θ2，…｝）。在每次更新结束后，每个质点滤波就会将其当前的参数估计传递给所有其他质点滤波。在t时刻，每个质点滤波（pfi）的后验概率计算公式如下：

这种近似分解了5n状态空间，其中n为同时发生的地震的数目。这就有效地减少了计算量。然而，由于所有质点包含的仅仅是整个参数空间的一小部分，因此它只是较优的方法。

这个试探法通过较大的阈值来过滤掉噪声使每个台站的P波记录初始化成一个新的质点滤波。由于可以根据被另一质点滤波追踪的现有事件进行局部监测，因此有必要对单个度量进行新的初始化。这个度量的自然选择是P ［z｜］，即触发台站的测量可由现有事件解释的概率。我们可以根据单个台站的似然函数（公式12）直接计算出这个度量。然而在这种情况下确定Aexp是非常重要的，因为它涉及计算多个波前干扰的加性效应。在这里我们提出了另一种可替代的度量，可以加快计算的速度。这种度量是由现有任何事件和最大概率时的阈值引起的震动的概率，公式如下：

通过调整阈值τ，我们可以使得这个系统在描述新的地震事件时更加完善。完整的算法在附图A1列出以供参考。

3 结果

我们根据数据与处理一节表述的数据，采用第一个被触发台站周围的平坦先验和每个质点滤波定义1 000个质点，进行了质点滤波参数估计近似。这个方法以1s的时间间隔进行更新，而且所有的实验都在模拟的实时中进行。

例子1：2011年3月15日，1∶36∶00～1∶38∶00（两个小地震）。

在这个时段我们进行了20次实验。其中一次实验的质点分布快照如图4所示。图5给出了这20次实验的参数估计结果（实线表示）和日本气象厅统一地震目录中公布的官方结果（虚线表示），其中，标准差用误差棒表示，x轴表示时间（秒）。根据实验结果，第一个质点滤波开始于P波初至，15s以后形成另一个质点滤波。这个方法成功地监测到了两个独立地震事件。另外，所有的估计参数在初始化后10s内收敛。平均来说，这个方法能够确定20km内的震中并可产生相对于日本气象厅统一地震目录误差为±1的震级估计值。

图5 2011年3月15日1时36分07秒到1时36分37秒期间20次独立运行实验的结果。时间历程为：（a）定位误差，（b）震级，（c）震源深度，（d）地震事件的发震时刻。这两个地震事件对应于图4（d）中的两个地震事件。这20次实验的平均运行时间历程标为实线，出现在日本气象厅地震目录中的官方数值标为虚线。所有实验的标准偏差用误差棒表示。在x轴显示的时间与最早事件的初至波对应

例子2：2011年3月20日，14∶19∶00～14∶21∶00（两个小地震）。

我们对例子2的数据进行了同样的分析，其中发生了相隔时间为5s的两次地震。相应的质点分布快照图和参数估计结果见图6和图7。在这个例子中，由于第一次地震发生在近海且近源记录很少，因此震源定位和参数估计较例子1更具挑战性。实验结果显示，参数估计的收敛速度很慢（事件A的收敛时间约为30s），定位误差也相当大（事件A的误差约为80km）。尽管如此，利用这个算法仍然可以有效地区分这两个独立地震事件，且估计的震级误差为±0.5之内。

例子3：2011年3月11日，14∶46∶00～14∶49∶00（东北地震）。

我们用东北地震的数据集说明，这个方法对于单个地震事件也是有效的。相应的质点分布快照图和估计参数的时间序列见图8和图9。由于这次地震为近海地震，因此在参数的初步估计中，有相当大的定位误差。但是在初至P波到达后40s，平均误差开始下降。随着地震破裂扩展，震级估计值由6.0增大到了8.4，这与地震破裂的机制相吻合。而估计的5个参数值均与日本气象厅地震目录中的值接近。

图6 2011年3月20日14时19分56秒之后2 000个质点分布的可视化图：（a）2s，（b）7s，（c）17s，（d）37s。图中的符号与图4相同。本图的彩色版只适用于电子版本

图7 2011年3月20日14时19分56秒到14时20分36秒期间20次独立运行实验的结果。各子图及其内的符号与图5定义的相同

4 讨论

当多个地震事件的发生在时间和空间上相距较远时，目前检测并将多个地震相关联的日本气象厅方法效果很好。然而当它们之间发生的时间和空间非常接近时，这些方法就会显示出产生很多虚报（Sagiyaet al，2011）。经验研究表明，基于质点的方法可以成功地区分多个同时发生的地震，并能提供合理的参数估计结果。同时，如果在似然中加入P波到时，即实测与预测P波到时之间的残差，收敛速度也可以得到有效的提高。这些结果说明，对于内陆地震，估计的参数收敛时间小于10s。而对于近海地震，估计收敛时间则在20～30s之间。至于定位误差，我们发现内陆地震小于20km，近海地震在20～80km之间。

为了区分多个同时发生的地震，使用非触发台站的信息非常重要。当前的日本气象厅地震预警系统在计算震源中仅使用触发台站波的到时。因此当多个地震近乎同时发生且后面的地震发生时间与前面地震的波的到时非常接近时，地震预警系统就将这些事件当成一个单一地震事件处理。如果出现这种情况且后面事件周围的台站观测到不可忽视的震动，则当前的预警系统就可能会高估其震级，因为这些台站距估计的震源（即前面事件的位置）很远。在我们提出的方法中，似然函数使用的信息不仅来自于触发台站，同时也来自非触发台站。这种设计加上Anoise的自适应测度就可以使该算法识别地震事件之间未受影响的区域，因而对于区分多个同时发生的地震是非常关键的。

图8 2011年3月11日14时46分46秒之后1 000个质点分布的可视化图：（a）2s，（b）7s，（c）13s，（d）22s，（e）32s，（f）62s。图中符号的定义与图4相同。本图的彩色版只适用于电子版本

我们方法的另一个优势是通过使用归一化质点滤波，可以有效地规避传统网格搜索中所需要的复杂计算。尽管如模型一节所述仍然需要先验分布，但这种分布可以根据历史记录得到。另外，也可以用初步的测量结果来选取合适的先验分布，从而得到更好的性能（Liu et al，2011）。但是这个方法也存在一些不足的地方，例如算法对先验分布的选择、质点的数目、以及Anoise、σnoise、σP和σS等的取值都非常敏感。尽管这些值可以进行实时调整，但是却需要大量的经验研究和对历史记录的分析。结果一节所述的一些慢收敛速度和高方差结果也可能是这些参数的次优选择造成的。

图9 2011年3月11日14时46分46秒到14时48分46秒期间15次独立运行实验的结果。各子图及其内的符号与图5定义的相同

在本文中，我们仅使用了3个例子来测试所提出的算法，因此我们还需要利用过去几年的多个同时发生的地震实例来进行实验研究，以评估这个方法。

另外，当有多个波前重叠发生时，如何模拟地震动会影响到多个地震事件的参数估计效果。在本文模型一节提出的算法（见图A1）中，并没有考虑这种模型。尽管在地震事件发生空间相距很远（大于100km）的例子研究中，这个省略只会带来很小的不同，但如果我们想用同一方法来区分发生时间非常接近的主震和余震，则应考虑这种模型。

5 结论

在地震活跃期间，震源类似或相距遥远的多个地震可以几乎在同一时间发生。如果不能将它们准确地区分开来，就会导致相应参数估计错误，进而还会导致虚报。本文中我们研究了发生时间接近的多个地震的鉴别和区分问题。基于考虑在任意给定时间有一个以上事件发生可能性的贝叶斯公式，我们提出了一种适用于区分多个同时发生地震的似然函数，同时还提出了一种序贯蒙特卡罗探试法，其复杂性与地震事件的数目呈线性增加关系。该探试法的性能用2011年东北地震后的3组日本气象厅地震记录数据进行了实验测试。初步测试表明，该方法能够成功区分在时间和空间上非常接近的多个地震，同时还可以实时估计相应的参数，且这些参数值与日本气象厅公布的官方数据在精度上相当。尽管该方法在用于实时检测之前还需要进行全面的验证，但初步研究表明，我们的方法可以减少地震震级高估的机会，因而有助于设计更好的地震预警系统。

6 数据与来源

本研究所用的波形数据选自日本气象厅（JMA）强震台网之内台站的连续记录。日本气象厅地震预警三个例子的性能来自http：／／www.data.jma.go.jp／svd／eew／data／nc／pub ＿ hist／2011／03／20110315013605／content／content＿out.html，http：／／www.data.jma.go.jp／svd／eew／data／nc／pub ＿hist／2011／03／20110320141959／content／content＿out.html，和 http：／／www.data.jma.go.jp／svd／eew／data／nc／pub ＿hist／2011／03／20110311144640／content／content＿out.htm（最后访问时间2013年7月）。数据处理过程我们使用了地震分析程序软件（http：／／www.iris.washington.edu／software／sac／manual／fileformat.html， 最 后 访 问 时间2013年7月）。日本气象厅的衰减关系来自第二届日本气象厅地震预警评估委员会的报 告 （http：／／www.data.jma.go.jp／svd／eqev／data／study－panel／eew－hyoka／t02／shiryou.pdf，最后访问时间2013年7月）。

7 附录

图A1中的算法描述了多个地震事件检测的归一化质点滤波的概要。 “CONVERGED”（收敛）的判据可由理想条件来取代（例如在估计‖＾θt－10－＾θt－1‖＜δ 时发生变化）。

附图A1 适用于多个地震事件的贝叶斯参数估计系统的质点滤波伪代码