利用烈度数据点估算中强历史地震震级与震中的方法研究1

吴 清 高孟潭

（中国地震局地球物理研究所，北京 100081）

引言

历史地震是分析区域地震活动特征，预测地震活动趋势不可缺少的基础资料。工程抗震设计工作中提出的地震危险性分析方法，其历史重演原则和地质构造类比原则也在很大程度上依赖于历史地震。近年来，历史地震研究取得了丰硕的成果，但始终存在着史料不均衡、不完备，地震参数误差较大等问题。即使相同的史料，不同专家的理解也可能产生分歧，不同版本地震目录所确定的同一地震参数可以出现较大差异。而新一代地震区划图的编制，对历史地震参数精度、可信度都有更高的要求。因此，有必要探讨估算历史地震参数的新思路和新方法。

地震烈度是指某一地区地面及各人工建筑物遭受一次地震影响的强烈程度（胡聿贤，2006）。公元1900年以前，在没有仪器可以记录地震动各种物理参数的时候，地震烈度资料就成为记录历史地震的唯一指标。历史文献中关于地震破坏的记录，以及震后烈度调查资料，是某个具体地点的破坏程度或烈度记录，可以据其评定该地点的烈度值，称为烈度数据点。烈度数据点从某种意义上说与地震台站有一定的相似处，仿照台站数据处理方法，用烈度数据点对一些主要地震参数进行定量估计，理论上切实可行（谢毓寿，1991），更是现今研究历史地震的重要手段。

基于Peruzza（1992）提出的残差最小化方法，Bakun等（1997）以烈度数据点为基础，利用衰减关系通过网格搜索来试算历史地震震级和震中区域，这种方法曾应用于美国加州和美国西北太平洋地区（Bakun等，1997；Bakun，1999；Bakun等，2003；Bakun，2006a；Diane，2009）、日本（Bakun，2005）、德国（Hinzen等，2001）、法国（Bakun等，2006b）及我国华北地区（张扬等，2009）。但该方法是以圆烈度衰减关系为基础，没有考虑烈度衰减的方向性问题，而且宏观震中与仪器震中有一定的混淆。Gasperini等（1999）用一套基于烈度数据点的Boxer方法来估算大震级历史地震的震源参数。但该方法首先将最高烈度点的几何中心定为震中，再分步得到其它参数，没有考虑其它烈度点分布对震中位置的影响；而且震级也是用圆烈度衰减关系进行估算。Musson等（2008）综合了以上几种方法的优点，采用Kövesligethy（1906）模型网格搜索最小均方根来估计震中，然后用 Frankel（1994）的有感面积方法确定震级大小，这个过程仍然没有考虑地震衰减的方向性。

一般来说，烈度随震级的加大而加大，随着到震中距离的加大而减小，而且震害分布总是呈现出明显的椭圆形状。我国研究者注意到这个现象并用椭圆模型来建立烈度衰减关系（沈建文等，1989；陈达生等，1989；霍俊荣，1989；胡聿贤，1999；汪素云等，2000）。但这些衰减关系是以等震线为基础数据回归得来，而等震线是原始烈度点的外包线，不是对烈度值的平均估计。霍俊荣（1989）曾提出，相同距离处另赋等震线低一度值的方法，近似等效烈度数据点来回归椭圆衰减关系，但这种等效法只是在等震线基础上进行数据处理，不是对真实烈度点的应用。

因此，建立一种直接基于烈度数据点的方法，充分考虑地震烈度的椭圆分布，尽可能利用全部烈度数据点来估算历史地震宏观震中和震级，对于提高历史地震参数的精度具有重要意义。

1 方法原理

本文首先以烈度数据点为基础，建立一种椭圆地震烈度分布模型来描述烈度随震级和距离而变化的关系，即假若已知一次地震的震级，此模型应能给出各地表点的烈度。然后基于所建烈度分布模型，联立考虑中心点和方向性的地震参数估计方程，代入烈度数据点信息，估算地震震级和震中。主要的技术路线如图1所示。

图1 本文技术路线示意图Fig. 1 The schematic diagram of the technology in this paper

一般来说，对于某一次地震，震级越大，距离震中越近，其地震烈度越大。我国地质条件复杂，烈度分布随地区变化差异较大。一般看来，地震烈度呈椭圆分布，即沿发震断层方向上烈度分布范围较广，而与断层垂直方向上烈度分布较窄，这种方向上的差别到了远场逐渐消失，最终趋于圆形。为了反映出这种烈度分布，本文借鉴陈达生等（1989）提出的地震烈度椭圆衰减关系长短轴统一回归的思想，采用如下椭圆烈度分布模型：

写成长、短轴方向的形式为：

沿长轴方向

沿短轴方向

式中，I为地震烈度；M为震级；系数1aC 、1bC 、2C、3aC 、3bC 均为回归常数；oaR 、obR 分别为椭圆长、短轴两个方向上烈度近场饱和因子；ε为回归分析中表示不确定性的随机变量，通常假定为对数正态分布，其均值为0，标准差为σ。

这里需要强调指出的是，这个模型以烈度数据点为基础，先用最小二乘法对地震各烈度区原始烈度点空间分布进行椭圆拟合，得到各个烈度区烈度点空间展布的椭圆拟合线，赋予这条椭圆拟合线该烈度区的烈度值，本文称之为烈度估计线。此时的aR、bR分别是烈度I的椭圆估计线长、短半轴长度。

该模型在震中处的烈度值随震级的变化率保持为常数C2，所以在任意震级长、短轴方向上震中处的烈度值均相等，从而保证回归后不需要通过参数调整来解决近场烈度差异的问题（肖亮，2011），而中间距离仍保持长、短轴烈度的差别，同时在远场也使烈度分布成圆形。

考虑中心点和方向性的椭圆数学方程见式（4），本文用此方程来估计地震震级与宏观震中：

式中，(xi, yi)为椭圆上的点坐标； ( x0, y0)为椭圆中心坐标； θ (0° ≤ θ ≤ 1 80°)是椭圆长轴逆时针与x轴正方向的夹角； Ra、 Rb分别为椭圆长半轴和短半轴长度。

由式（2）得：

由式（3）得：

将式（5）和式（6）代入到式（4）中得：

在同一个烈度点应有：

将式（8）代入式（7）得：

在式（9）中， C1a、 C1b、 C2、 C3a、 C3b、 Roa、 Rob由椭圆烈度分布模型统计回归可得。若已知一个地震的多个烈度信息点，带入到式（9），即可联立方程组求得震中震级M和方向θ。有4个未知数，则至少需要4个方程确定一组解，但 x0与 y0不是相互独立的，因此在已知椭圆烈度分布模型下，3个烈度信息点就能确定一组震中、震级和方向。烈度信息点增多，变成求解超定非线性方程组，采用通用全局优化法的数值计算方法寻求最优解，可非常直观地求得地震宏观震中位置、震级和方向。

2 地震数据及椭圆烈度分布模型

我国目前公开发表的烈度调查资料大多采用等震线的形式发布，原始烈度调查点通常难以获得。本文选取最新出版的4册《中国震例（1992—2002）》（陈棋福，2002a；2002b；2003；2008），提取其中既有仪器测量数据又有宏观烈度调查数据的现代地震，将烈度分布图数字化。这些地震的等震线在图上部分明确标示了宏观烈度调查点位置；对于未明示烈度调查点位置的，本文将等震线在图上能明确获取烈度信息的城、县、镇、村所在地算作烈度点。经数字化配准后，获取各地震烈度点的空间分布。受到提取烈度点有效数字信息的目的所限，本文获取的地震基本分布在青、甘、川、滇4省，地震震级集中在MS4.9—6.6级，如表1所示。

将式（4）改写为考虑中心点和方向性的椭圆参数方程：

对于既有仪器测量数据又有宏观调查烈度数据的现代地震，可获得地震的宏观震中位置(x0, y0)、极震区走向θ和烈度分布点坐标根据各烈度区调查点的空间分布基于式（10）利用最小二乘原理采用通用全局优化算法拟合出各烈度区椭圆估计线长半轴和短半轴长度 Ra和 Rb，由此获得各烈度区的椭圆烈度估计线。由于地震的宏观震中并不一定都是该地震各烈度区的几何中心，因此对某些地震需要通过式（10）同时拟合出各烈度区的几何中心。这里的烈度估计线长轴方向均取极震区走向θ。

以各自宏观震中为原点（0，0），将各地震对应烈度点都转换到了平面坐标下。表1列出了拟合所得的各地震、各烈度区椭圆烈度估计线的长半轴和短半轴长度，拟合烈度估计线几何中心不同于宏观震中的也同样列在了表1中。表1中拟合烈度估计线的几何中心与宏观震中重合的，坐标即为（0，0）；几何中心与宏观震中不重合的，再具体列出。

表1 地震数据及各烈度区椭圆烈度估计线长、短半轴长度和中心点Table 1 List of parameters of the earthquakes in the study area

续表

续表

将表1中各地震、各烈度区的椭圆烈度估计线的长、短半轴长度及对应烈度值和震级值，代入到式（2）、（3）进行最小二乘统计回归，得到地震烈度分布模型：

不同地区的震源特性、传播介质与场地条件都可能不同，地震烈度分布规律自然也可能不同。因此，烈度分布模型以按大地区分别使用为宜。受获取地震原始烈度调查点所限，本文所得的烈度分布模型仅适用于青、甘、川、滇4省的中强地震。

3 震例验算

由上节得到椭圆烈度分布模型的回归系数后，式（9）即变为：

为了验证地震参数估计方程的有效性，本文搜集了《2001—2005年中国大陆地震灾害损失评估汇编》（中国地震局震灾应急救援司，2010）里的地震烈度资料，提取了5个发生于模型适用区域的地震的烈度数据点，将所有烈度数据点以各自地震的宏观震中为原点（0，0）转换到平面坐标下（xi，yi，Ii），然后直接代入到式（13）中进行验算，计算结果见表2。这5个地震空间位置在研究区内分布合理，震级也具有一定的代表性。

表2 现代震例计算结果Table 2 Estimated results of the modern earthquakes

由表2可以看到，参与验算的5个地震，计算震级与仪器震级的差值不超过0.6级，计算震中到宏观震中的距离不超过6km，由此可见此方法是可行的。

4 不确定性分析

在震中位置和震级的估算过程中，烈度点的数量、分布和精度好坏将直接影响估算结果。由于历史地震烈度点相对较少，为了对历史地震参数进行较好的估计，本文运用蒙特卡洛方法，对信息丰富的现代地震烈度点抽样，模拟历史地震烈度点分布情况。由于至少需要3个烈度点才能确定一组解，先从烈度点中随机抽取3个点，并且有放回的重复抽样1000次。增加控制信息，逐点添加烈度点数，随机重复抽取4，5，6……个点各1000次，这样模拟出大量的烈度点组样本。每个烈度点有可能被多次选择或者不被选择，每次抽样代表不同的烈度点空间分布，而不同的烈度点数代表了烈度信息的多寡，真实历史地震资料记载只相当于随机抽样中的一种情况。通过蒙特卡洛方法模拟出大量烈度点组样本代入式（13）中计算，由此来分析此方法计算地震震级与震中的不确定性。可以推测，烈度点数越多，分布越均匀，得出的震中与震级越精确，不确定性越小。

式（13）是超越方程，当参与计算的烈度点数多于3个点时，更成为超定的超越方程组。这类矛盾方程组在严格的数学意义下无解，但是可以在最小二乘意义下，最大限度拟合数据，寻求最优解。随机抽样得到的烈度点组会存在不同的空间分布，利用数学软件采用全局最优化的数值计算方法，可以得到每组抽样点满足式（13）的最优数学解，但得到的数学解有可能并不满足所求参数的物理意义。因此，所得的结果需要剔除与地震参数物理意义不符的数学解，然后对最终结果进行统计分析。

本文以2001年5月24日云南宁蒗—四川盐源5.8级地震为例，给出抽样计算结果。宁蒗—盐源地震有50个烈度点，以宏观震中为原点转换到平面坐标下，然后从这50个烈度点里随机抽取3点、4点、5点……各1000次，分别带入到式（13）进行解算，计算震中分布见图2。将计算震级与仪器震级做差值，差值分布见图3。由图2可以看到，随着烈度点的增多，计算震中逐渐向宏观震中靠拢。由图3可以看到，随着烈度点的增多，计算震级与仪器震级的差值也逐渐收敛。

图2 2001年5月24日云南宁蒗—四川盐源5.8级地震抽样计算震中分布（R为计算震中到宏观震中平面距离）Fig. 2 Distribution of calculated epicenters of Ninglang-Yanyuan MS 5.8 earthquake occured in Yunnan-Sichuan on 24 May, 2001（R is the distance between calculated epicenter and macro-epicenter）

图3 2001年5月24日云南宁蒗—四川盐源5.8级地震计算震级与仪器震级差值分布Fig. 3 The distribution of the differences between calculated magnitudes and instrument magnitude of Ninglang-Yanyuan MS 5.8 earthquake occured in Yunnan-Sichuan in 24 May，2001

将参与试算的5个地震均进行抽样计算，统计不同抽样点数情况下，计算震中到宏观震中距离的均值与标准差，以及计算震级与仪器震级差值的均值与标准差，见图4（a）、（b）。

由图4（a）可以看到，随着烈度点的增多，计算震中到宏观震中距离的均值和标准差都逐渐减小。3点时，距离均值最大有18km，标准差最大接近12km；6点时，标准差到10km以下；8点时，均值也到10km以下。历史地震震中精度分布为1类≤10km，2类≤25km，3类≤50km，4类≤100km，5类>100km。所以按照本文所提算法，一般达到6点时，标准差就能达到10km以下，即达到历史地震的1类精度。

由图4（b）可以看到，随着烈度点的增加，震级差值的均值有个增大过程，达到一定点数之后趋于平稳。这是因为烈度点数较少时，代表的地震影响范围相对较小，由此估算所得的震级自然偏小；而达到一定点数之后，烈度点分布所能反映的影响范围趋于稳定，计算出来的震级也就趋于稳定。模型本身是对一个地区烈度分布的综合估计，与任何一个真实地震之间都存在着系统偏差。不同地震之间的震源特性与场地条件不同，系统偏差也各不相同。而随着烈度点的增加，震级差值的标准差整体上呈减小趋势。对于参与试算的 5个地震，3点时，震级差值的标准差最大接近0.45级；8点时，标准差都集中到0.3级以下。

通过对《中国历史强震目录（公元前23世纪—公元1911年）》（国家地震局震害防御司，1995）中历史强震烈度点数的统计，发现历史地震烈度记载点普遍较少，18个及以上烈度点的地震，仅占历史强震总数的8%。因此，本文暂且分析20个及以下烈度点，计算地震参数不确定性的情况。

图4 （a） 不同烈度点数情况下计算震中到宏观震中距离的均值与标准差分布图Fig. 4（a） Distribution of the mean and standard errors of distances between calculated epicenters and macro-epicenter under different numbers of intensity data points

图4 （b） 不同烈度点数情况下计算震级与仪器震级之差值的均值与标准差分布图Fig. 4（b） Distribution of the mean and standard errors of the magnitude difference between calculated and instrument recorded under different number of intensity data points

不同地震、不同烈度点数情况下，不确定性各有不同，但整体分布趋势基本相同。因此，可以按照烈度点数情况，对估定历史地震参数的不确定性做粗略估计。

5 历史地震参数的估算

通过前面的分析，本文建立了一套基于烈度数据点，用椭圆烈度分布模型估算地震宏观震中与震级的方法，旨在用此方法来估算只有烈度记载且烈度数据点相对较少的历史地震参数。以云南地区为例，将此方法用于历史地震进行试算。从《中国历史地震图集》（明、清时期）（下文简称《图集》）（国家地震局地球物理研究所等，1986；1990）中提取4个烈度点较少的中强地震，获取有历史记载的烈度数据点，然后用本文所提算法估算历史地震参数。

5.1 1599年10月16日云南嵩明5.0级地震

万历二十七年八月二十八日（公元1599年10月16日）云南嵩明地震，《图集》给出地震参数为震中（103.0°E，25.3°N），震级5.0级，有5个烈度记载点，图5是其烈度记载点分布图。用本文的算法计算，计算震中为（102.99°E，25.28°N），与《图集》震中相距2.42km；计算震级为5.2级，比《图集》震级大0.2级。

图5 云南嵩明5.0级地震烈度点分布图Fig. 5 Distribution of intensity data points of Songming MS 5.0 earthquake in Yunnan

5.2 1696年7月7日云南昆明5.5级地震

康熙三十五年六月初九（公元1696年7月7日）云南昆明地震，《图集》给出地震参数为震中（102.7°E，25.1°N），震级5.5级，有12个烈度记载点，图6是其烈度记载点分布图。用本文的算法计算，震中为（102.68°E，25.07°N），与《图集》震中相距 3.8km；计算震级为5.2级，比《图集》震级小0.3级。

图6 云南昆明5.5级地震烈度点分布图Fig. 6 Distribution of intensity data points of Kunming MS 5.5 earthquake in Yunnan

5.3 1876年8月5日云南永平6.0级地震

光绪二年六月十六日（公元1876年8月5日）云南永平地震，《图集》给出地震参数为震中（99.5°E，25.5°N），震级6.0级，有8个烈度记载点。图7是其烈度记载点分布图。用本文的算法计算，计算震中为（99.45°E，25.55°N），与《图集》震中相距 7.04km；计算震级为5.9级，比《图集》震级小0.1级。

图7 云南永平6.0级地震烈度点分布图Fig. 7 Distribution of intensity data points of Yongping MS 6.0 earthquake in Yunnan

5.4 1909年5月11日云南米勒西南6.5级地震

宣统元年三月二十二日夜（公元1909年5月11日）云南米勒、宁州地震，《图集》给出地震参数为震中（103.1°E，24.3°N），震级6.5级，有16个烈度记载点，图8是其烈度记载点分布图。用本文的算法计算，计算震中为（103.11°E，24.30°N），与《图集》震中相距0.63km；计算震级为6.8级，比《图集》震级大0.3级。

图8 云南米勒西南6.5级地震烈度点分布图Fig. 8 Distribution of intensity data points of Mile MS 6.5 earthquake in Yunnan

由以上震例可见，本文的方法对烈度数据点较少、等震线不易勾画的历史地震颇为有效。处理过程直接明了，减少了等震线勾画过程中的主观不确定性，可重复和再现。而且对于1599年10月16日云南嵩明5.0级地震和1876年8月5日云南永平6.0级地震，这类烈度点分布极不均匀、基本集中在极震区一侧的地震，用本文的算法也能得到较为可信的地震参数。由此可见，本文的方法对烈度信息空间分布散乱的稳定性。

6 结论与讨论

本文首先基于地震原始烈度调查点，用椭圆数学方程对烈度点空间分布进行最小二乘拟合，得到各烈度区椭圆烈度估计线。进而对烈度估计线进行统计回归，得到一组适用于青、甘、川、滇4省中强地震的椭圆烈度分布模型该模型的建立充分利用了大量空间分布的原始烈度信息，通过数学手段得到烈度估计线，避免了等震线圈定过程中的主观经验因素，可还原和重现。相较于烈度外包线的地震等震线，基于烈度点的椭圆估计线更能代表烈度分布的真实情况。

根据所建立的椭圆烈度分布模型，联立考虑中心点和方向性的椭圆数学方程，代入烈度数据点信息，直接求解历史地震震中和震级，过程直接明了。通过蒙特卡洛方法，定量分析所得参数的不确定性，可以看到随着烈度点的增多，计算震中到宏观震中距离的均值和标准差都逐渐减小，计算震级与仪器震级差值的均值有个增大过程，达到一定点数之后趋于平稳。

本文的算法符合椭圆模型的数学、物理概念，满足烈度分布的地震学规律，充分利用了大量空间分布的原始烈度信息，过程简单直接。算例表明，计算结果与实际地震记录参数在一定范围内吻合。将此方法应用到仅有烈度资料的历史地震，可以减少人为判定震中位置和震中烈度，进而估算震级过程中人为的主观不确定性，提高历史地震震中和震级的精度。

由于中国大陆的地震绝大部分发生在地壳以内，其震源深度差别不大（刘百篪等，2002），因此本文暂且不考虑震源深度对烈度分布的影响。

地震震害一般呈椭圆分布，但不一定是数学上的椭圆。这里将地震烈度分布近似为椭圆系列时，已经包含了数学抽象的不确定性。但即要建立模型，简化及随之而来的不确定性或误差就不可避免（沈建文等，1989）。地震烈度分布模型本身是从多个地震中提取的平均规律，与地震资料之间存在着离散性，与任何一个真实地震之间都有系统偏差，在以后的工作中需要进一步改进模型，以期减小这种系统偏差。

本文所建模型至少需要3个烈度点才能求解，而实际历史地震资料中，仅有1个或2个烈度点的地震占了资料的半数以上，这类地震震中参数该如何获取，需要进一步讨论。另外，式（13）中同时解算出了椭圆长轴与x轴正方向的夹角θ，这个夹角θ应该与发震构造或者破裂参数有关，但相关程度与其不确定性还有待进一步研究。

最后，本文所提算法的基础是地震烈度数据点。但是现代地震的烈度评定与历史地震的烈度评定之间也存在着差异，如何将现代地震烈度与历史地震烈度扩展衔接，也是值得深入思考的问题。

陈达生，刘汉兴，1989. 地震烈度椭圆衰减关系. 华北地震科学，7（8）：31—42.

陈棋福，郑大林，车时，2002a. 中国震例（1992—1994）. 北京：地震出版社.

陈棋福，郑大林，刘桂萍等，2002b. 中国震例（1995—1996）. 北京：地震出版社.

陈棋福，郑大林，高荣胜，2003. 中国震例（1997—1999）. 北京：地震出版社.

陈棋福，郑大林，车时等，2008. 中国震例（2000—2002）. 北京：地震出版社.

国家地震局地球物理研究所，复旦大学中国历史地理研究所，1986. 中国历史地震图集（明时期）. 北京：地图出版社.

国家地震局地球物理研究所，复旦大学中国历史地理研究所，1990. 中国历史地震图集（清时期）. 北京：中国地图出版社.

国家地震局震害防御司，1995. 中国历史强震目录（公元前23世纪-公元1911年）. 北京：地震出版社.

胡聿贤，1999. 地震安全性评价技术教程. 北京：地震出版社.

胡聿贤，2006. 地震工程学. 北京：地震出版社.

霍俊荣，1989. 近场强地面运动衰减规律的研究 [博士论文]. 黑龙江：国家地震局工程力学研究所.

刘百篪，郑文俊，郭华等，2002. 活断层工作方法在中早期历史地震研究中的作用. 中国地震，18（3）：283—288.

沈建文，华宜平，1989. 关于地震烈度衰减模型的系统偏差. 地震学报，11（1）：38—45.

汪素云，俞言祥，高阿甲等，2000. 中国分区地震动衰减关系的确定. 中国地震，16（2）：99—106.

谢毓寿，1991. 历史地震参数的估定. 东北地震研究，7（5）：1—6.

肖亮，2011. 水平向基岩强地面运动参数衰减关系研究 [博士论文]. 北京：中国地震局地球物理研究所.

张扬，马干，史保平等，2009. 华北地区烈度衰减模型建立及其用于震中区域和震级的定量估算. 地震学报，31（3）：290—306．

中国地震局震灾应急救援司，2010．2001—2005年中国大陆地震灾害损失评估汇编．北京：地震出版社，1—560．

Bakun W.H.，Wentworth C.M.，1997. Estimating Earthquake Location and Magnitude from Seismic Intensity Data.Bull. Seism. Soc. Amer.，87（6）：l502—1521．

Bakun W.H.，1999. Seismic Activity of the San Francisco Bay Region. Bull. Seism. Soc. Amer.，89（3）：764—784.

Bakun W.H.，Johnston A.C.，Hopper M.G..，2003. Estimating Locaitons and Magnitudes of Earthquakes in Eastern North America from Modified Mercalli Intensities. Bull. Seism. Soc. Amer.，93（1）：190—202.

Bakun W.H.，2005. Magnitude and Location of Historical Earthquakes in Japan and Implications for the 1855 Ansei Edo Earthquake. J. Geophys. Res.，110：B02304.

Bakun W.H.，2006a. Estimating Location and Magmotudes of Earthquakes in the Southern California from Modified Mercalli Intensities. Bull. Seism. Soc. Amer.，96（4）：1278—1295.

Bakun W.H.，Oona Scotti，2006b. Regional intensity attenuation models for France and the estimation of magnitude and location of historical earthquakes. Geophysics J. Int.，164：596—610.

Diane I. Doser，2009. Estimating Magnitude and Location of Alaskan Earthquakes Using Intensity Data. Bull.Seism. Soc. Amer.，99（6）：3430—3453.

Frankel A.，1994. Implications of Felt Area-magnitude Relations for Earthquake Scaling and the Average Frequency of Perceptible Ground Motion. Bulletin of the Seismological Society of America，84：462—465.

Gasperini P.，Bernardini F.，Valensise G. et al.，1999 Defining Seismogenic Sources from Historical Earthquake Felt Reports. Bull. Seism. Soc. Amer.，89（1）：94—110．

Hinzen K.G.，Oemisch M.，2001. Location and Magnitude from Seismic Intensity Data of Recent and Historic Earthquake in the Northern Rhine Area，Central Europe. Bull. Seism. Soc. Amer.，91（1）：40—56．

Kövesligethy R.D.，1906. A makroszeizmikus rengések feldolgozása. Mathematikai és Természettudományi Értesítõ，24：349—368.

Musson R.M.W. and Jeménes M.J.，2008．Macroseismic estimation of earthquake parameters．Sixth Framework Programmer report．

Peruzza L.，1992. Procedure of macroseismic epicentre evaluation for seismic hazard purposes. See：Proceedings of the XXIII General Assembly of the ESC，Prague. 2：434—437.