吴清 高孟潭

中国地震局地球物理研究所，北京市海淀区民族大学南路5号 100081

0 引言

地震监测预报和工程地震是地震工作者为了减轻地震灾害面临的两大艰巨任务（谢毓寿，1991），地震活动性研究是上述工作的基础。由于地震仪器记录的历史较短，不足以在短时间、小尺度范围内弄清楚地震活动规律，因而对历史地震的研究非常重要。尤其是对于6.5级以上的历史强震，其参数的不确定性，会影响人们对某个地点或地区地震活动状况的认识，影响某项重大工程的地震危险性评估结果，更可能成为地震科学研究或防震减灾工作中一些重大问题的关键。

6.5级以上强震的烈度一般在破裂方向上衰减较慢，在垂直破裂方向上迅速衰减，地震影响场呈狭长的带状分布，普通的圆烈度模型和椭圆模型都不再适用，更适用的是断层破裂模型，但是断层破裂模型的带状跑道型分布较为复杂，在数学上不易处理。根据6.5级以上强震烈度分布的图像特征，本文提出了基于烈度数据点的考虑断层破裂长度的烈度椭圆分布模型。即在烈度椭圆分布模型的基础上，加入断层破裂长度的约束，要求最内圈烈度估计等值线长轴长度大于断层破裂长度，使得烈度椭圆分布模型更接近于强震的真实烈度分布。基于既有仪器测定记录又有宏观考察数据的现代大震建立适用于强震的烈度分布模型，进而构建强震参数估算方程，直接利用全部烈度点的空间分布信息来估算历史强震的震级和震中，为历史强震参数的确定提供新的思路与方法，进一步提升地震史料的可用性，提高历史地震参数的可靠性，对于6.5级以上历史大震参数的校核具有重要意义。

1 基于离散烈度点考虑断层破裂长度的强震烈度分布模型

1.1 地表破裂长度与震级的经验关系式

Tocher（1958）是第一个把地震震级与相应的断层地表破裂长度联系起来的。他根据美国加利福利亚州及内华达州10次地震的数据，利用最小二乘法建立了形如M＝a+b log L的关系式，式中M为震级，L为地震引起的地表破裂长度（km）。后来，不少人相继提出震级与地震引起的地表破裂长度的类似关系式。陈达生（1984）对这类关系式做过系统介绍，并指出M＝a+b log L与log L＝c+dM这2个回归方程之间不是反函数关系。根据回归分析的定义，前者适用于估算已知断层地表破裂长度对应的地震震级；而后者适宜于估算可能与某地震震级伴生的地表破裂长度。

学者们曾建立了不少破裂长度与震级的经验关系，但或多或少都用到了历史地震数据或者跨国跨区域的地震资料。本文旨在基于既有仪器测定记录又有宏观考察数据的现代大震建立适用于我国的强震烈度分布模型，因此收集新的现代大震资料重新进行经验关系的统计拟合。自20世纪60年代以来，中国相继发生了好几次伴有断层的地表破裂的大地震，我们收集了这些大地震数据并按照log L＝c+dM进行了统计回归，以估算可能与某地震震级伴生的地表破裂的长度（表1）。

表1 用于统计破裂长度与震级经验关系的地震列表

将表1中所列地震的破裂长度与震级进行最小二乘统计回归，得到

据此，在已知震级时可估算伴随该地震的地表破裂长度。

1.2 强震烈度分布模型及其回归方法

地震烈度在地面的空间分布会呈现出一定的几何图像，依据空间统计分析的理论，可以建立一系列的地震烈度分布模型。

地震烈度分布模型描述烈度随震级（或震中烈度）和距离而发生的变化，即假若已知某次地震的震级，应能给出各地表点的烈度。目前主要的地震烈度分布模型有圆模型、椭圆模型和断层破裂模型。对于震级较小的地震，常常看不出等震线的方向，因此常常不考虑烈度分布的方向性而寻求一种平均关系，即圆模型。当震级较大时，等震线形状常呈现一种椭圆形，可考虑在长轴和短轴方向有不同的分布，即椭圆模型。对于发震断层出露地表的地震，高烈度区通常沿发震断层分布为带状，沿断层垂直方向衰减很快，为此研究者们提出了断层破裂模型。

对于6.5级以上的强震，比较适用于断层破裂模型，但是带状分布的模型较为复杂，在数学上不易处理，这里我们提出基于烈度数据点的考虑断层破裂长度的烈度椭圆分布模型。构建该模型的步骤为：

①由1.1节拟合得到的破裂长度-震级关系式（1）估算地震可能伴随的地表破裂长度；

②直接从原始烈度点的空间展布出发，用最小二乘法对地震各烈度区原始烈度点空间分布进行椭圆拟合，得到各烈度区烈度点空间展布的椭圆拟合线，赋予这条椭圆拟合线该烈度区的烈度值，本文称之为烈度估计等值线。要求最内圈烈度估计等值线椭圆长轴要大于由公式（1）估计的地表破裂长度。

③以烈度估计等值线的长短半轴长度为基础数据进行统计回归得到烈度椭圆分布模型。

此模型是对原始烈度点空间展布的综合估计，充分利用了所有烈度点的原始空间分布信息，更接近真实地震的烈度分布。

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

1.2.1 强震烈度估计等值线拟合

中国目前公开发表的烈度调查资料大多采用等震线的形式发布，原始烈度调查点通常难以获得。本文选取已经出版的《中国震例（1966—2002）》（张肇诚等，1988、1990a、1990b、1999、2000；陈棋福等，2002a、2002b、2003、2008）中 2008年以后的强震烈度资料从中国地震局官方网站的地震专题上获取（http：／／www.cea.gov.cn／publish／dizhenj／468／553／index.htm l），收集了20世纪60年代以来全国所发生的M≥6.5的地震等震线图，提取了其中既有仪器测量数据又有宏观烈度调查数据的现代地震资料，并将烈度分布图数字化。这些地震的等震线图上部分明确标示了宏观烈度调查点位置；对于未标示烈度调查点的，本文将能在等震线图上明确获取烈度信息的城、县、镇、村所在地作为烈度点。经数字化配准后获取各地震烈度点的空间分布（表2）。

考虑中心点和方向性的椭圆参数方程为

式中，（xi，yi）为椭圆上点的坐标，（x0，y0）为椭圆中心坐标；θ（0°≤θ≤180°）是椭圆长轴逆时针与x轴正方向的夹角；Ra、Rb分别为椭圆长半轴和短半轴长度，ti（0≤θ≤2π）为椭圆上各点所对应的参数。

表2 地震数据及各烈度区椭圆烈度估计等值线长、短半轴长度和椭圆中心

由式（1），可以估计相应震级的地震可能伴随的地表破裂长度L，以此约束椭圆强震烈度分布模型，要求最内圈烈度估计等值线椭圆长轴要大于由式（1）估计的地表破裂长度，即

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

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

图1是 2014年10月7日云南普洱景谷6.6级地震烈度估计等值线拟合结果，以宏观震中为平面坐标原点。图中小三角是烈度数据点，椭圆长轴方向为极震区走向，烈度估计等值线是对原始烈度点空间分布进行最小二乘拟合所得，可以看到烈度估计等值线趋向于原始烈度点空间分布的平均估计。

图1 2014年10月7日云南普洱景谷6.6级地震烈度估计等值线分布图

1.2.2 强震烈度分布模型

本文据陈达生等（1989）关于地震烈度椭圆衰减关系长短轴统一回归的思想，采用烈度分布模型

式中，I为地震烈度，M为震级；系数 a、b、C1、C2均为回归常数；R0a、R0b分别为椭圆长、短轴两方向烈度近场饱和因子；ε为回归分析中表示不确定性的随机变量，通常假定为对数正态分布，其均值为0，标准差为σ。该模型在震中处的烈度值随震级的变化率为常数b，而中间距离仍保持长短轴烈度的差别，同时在远场也使烈度分布成圆形。需要强调指出的是，这里Ra、Rb分别是烈度I的烈度估计等值线的长、短半轴长度，不再是等震线的长、短半轴长度。

将表2中各地震各烈度区的烈度估计等值线，按照式（4）进行最小二乘统计回归。由于拟合的烈度估计等值线已经是对烈度点分布的综合估计，所以没有采取近场补点。但为了体现远场区发震构造影响消失，烈度分布趋于圆形的特点，汪素云等（1993）提出取有感范围的半径作为远场控制点，有感烈度值通常为Ⅲ～Ⅳ度，在计算中取为Ⅲ度半（3.5度），称为远场补点。这里提出的有感半径有外包线的性质，而本文烈度模型主要是对烈度分布的综合估计，因此计算时远场控制烈度取为Ⅲ度。有感半径与震级的关系见表3。

表3 有感半径与震级的关系（汪素云等，1993）

由原始烈度点拟合得到烈度区的估计等值线，进而统计回归得到地震烈度分布模型为

2 地震参数估计方法及其不确定性分析

2.1 地震参数估计方法

本文在上节中建立了如下烈度分布椭圆模型：

将（2）式改写成椭圆标准数学方程式（即（9）式）

用以估计地震震级与宏观震中。式中，（xi，yi）为椭圆上点的坐标，（x0，y0）为椭圆中心坐标；θ（0°≤ θ≤180°）是椭圆长轴逆时针与x轴正方向的夹角。

已知烈度数据点平面坐标（xi，yi）和烈度值Ii，联立烈度分布模型及考虑中心点和方向性的地震参数估计方程（7）、（8）、（9），即可计算椭圆中心坐标即地震震中平面坐标（x0，y0）和震级M，以及椭圆长轴与x轴夹角θ。

由式（7）得

由式（8）得

将式（10）和（11）带入到式（9）中得

在同一个烈度点应有

将式（13）带入式（12）得

2.2 求解参数的基本条件

2.2.1 烈度点数量

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

2.2.2 烈度点空间分布

参数估计方程是以椭圆数学方程为基础，因此用于计算的烈度点不能分布于一条或接近于一条的直线上，否则在极端情况下将造成方程无解，或所得数值解不满足地震参数的物理意义。

2.2.3 适用范围

由于在建模过程中，8.0级以上特大地震的资料太少，因此本方法只适用于6.5～8.0级的历史强震，对于8.0级以上的特大地震需要另行讨论。

2.3 震例验算

2.3.1 内符检验

带入由1.2.2节得到的强震烈度分布模型系数，式（14）即变为

为了验证此算法的可行性，以参与强震烈度分布模型拟合的其中7个地震为例，将各自的烈度数据点（xi，yi，Ii）直接带入到式（15）中进行试算，计算时所有烈度点以各自地震的宏观震中为原点（0，0）转换到平面坐标下，计算结果见表4。

表4 内符验算计算结果

由表4可以看到，参与验算的7个地震，计算震级与仪器震级的差值不超过0.5级，计算震中到宏观震中的距离不超过25km，由此可见此方法对参与模型拟合的地震集是可行的。

2.3.2 外推检验

为了验证地震参数估计方程的外推有效性，我们提取了3个没有参与模型拟合的强震烈度数据点，将所有烈度数据点以各自地震的宏观震中为原点（0，0）转换到平面坐标下（xi，yi，Ii），然后直接带入到式（15）中进行验算，计算结果见表5。

表5 外推验算计算结果

由表5同样可以看到，参与验算的3个地震，计算震级与仪器震级的差值不超过0.7级，计算震中到宏观震中的距离不超过15km，由此可见此方法外推也是可行的。

2.3 不确定性分析

在震中位置和震级的估算过程中，烈度点的数量、分布和精度好坏将直接影响估算结果。由于历史地震烈度点相对较少，为了对历史地震参数进行较好的估计，本文运用蒙特卡洛方法，对信息丰富的现代地震烈度点抽样，模拟历史地震烈度点分布情况，进而分析本文所建方法的不确定性。

由前文知，至少需要3个烈度点才能确定一次地震的震级、震中和方向。采用蒙特卡洛方法，先从同一地震的烈度数据点中随机抽取3个点，并且有放回的重复抽样1000次。增加控制信息，逐点添加烈度点数，随机重复抽取4，5，6……个点各1000次，这样模拟出大量的烈度点组样本。每个烈度点有可能被多次选择或者不被选择，每次抽样代表不同的烈度点空间分布，而不同的烈度点数代表了烈度信息的多寡，真实历史地震资料记载只相当于随机抽样中的一种情况。通过蒙特卡洛方法模拟出大量烈度点组样本带入参数估计方程，完成计算后分析此方法计算地震震级与震中的不确定性。可以推测，烈度点数越多，分布越均匀，得出的震中与震级越精确，不确定性越小。

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

本文以2014年10月7日云南普洱景谷6.6级地震为例，给出抽样计算结果。首先获取地震的烈度数据点，然后以其宏观震中为原点转换到平面坐标下，再从地震烈度平面坐标点里随机抽取3点、4点、5点……20点各1000次，分别带入式（15）进行解算，图2是计算所得震中的分布图。由计算所得震级与仪器震级之差（ΔM）绘成图3。可以看到，随着烈度点的增多，计算震中逐渐向宏观震中靠拢，计算震级与仪器震级的差值也逐渐收敛。

将参与内符检验的7个地震和参与外推检验的3个地震均进行抽样计算，统计不同抽样点数情况下计算震中到宏观震中距离的均值与标准差，以及计算震级与仪器震级差值的均值与标准差（图4）。

图4的（a）和（b）是内符检验地震的计算结果。由图4（a）可以看到，随着烈度点的增多，计算震中到宏观震中距离的均值和标准差都逐渐减小。取3个点时，距离均值最大接近26km，标准差最大接近15km；均值基本都在25km以下，取20个点时，标准差都在10km以下。

图2 2014年10月7日云南普洱景谷6.6级地震抽样计算震中分布

图3 2014年10月7日云南普洱景谷6.6级地震计算震级与仪器震级差（ΔM）的直方图

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

图4的（c）和（d）是外推检验地震的计算结果。由图4（c）同样可以看到，随着烈度点的增多，计算震中到宏观震中距离的均值和标准差都逐渐减小。取3个点时，距离均值最大不足25km，标准差最大接近11km；取6点个时，标准差在10km以下。

图4（d）与图4（b）的趋势类似，可以看到随着烈度点的增加，震级差值的均值有个增大过程，达到一定点数之后趋于平稳。随着烈度点的增加，震级差值的标准差整体上呈减小趋势。对于参与试算的3个地震，取3个点时，震级差值的标准差最大接近0.45级；取7个点时标准差都小于0.4级。

度量不确定性的指标就是不确定度。综合不确定度是随机误差和系统误差的综合度量，通常用均方误差来表示，其表达式为

图4（a） 内符检验地震不同烈度点数情况下计算震中到宏观震中的距离（ΔR）的均值及其标准差分布

图4（b） 内符检验地震不同烈度点数情况下计算震级与仪器震级之差（ΔM）的均值及其标准差分布

图4（c） 外推检验地震不同烈度点数情况下计算震中到宏观震中的距离（ΔR）的均值及其标准差分布

图4（d） 外推检验地震不同烈度点数情况下计算震级与仪器震级之差（ΔM）的均值及其标准差分布

式中，σ为精度，代表误差分布的密集或离散程度，即观测结果与其数学期望的接近程度；ε为准确度，代表观测量的真值与数学期望之差，也就是系统误差；综合不确定度D代表的是精确度，是精度和准确度的合成，即观测值与真值的接近程度，包括观测结果与其数学期望的接近程度和数学期望与其真值的偏差。本文以综合不确定度来度量不确定性。

对于既有仪器记录又有宏观考察的现代地震，本文将仪器震级和考察所得的宏观震中作为真值，以均方误差来估计蒙特卡洛方法所得计算结果的不确定度。由于本文运用蒙特卡洛方法时采取的是随机抽样，因此烈度点的组合难免会抽到极端情况，使得所得结果不满足地震参数的物理意义，这里按照粗差处理的方式予以剔除。然后对处理后的结果估计其综合不确定度，结果见图5和6。

图5 内符检验地震参数估计的综合不确定度

图5和6显示的趋势大致相同。对于地震震中综合不确定度，除个别情况外，所有点数情况下震中不确定度都在25km以内。历史地震震中精度分为1类≤10km，2类≤25km，3类≤50km，4类≤100km，5类＞100km。对应到历史地震震中精度，按照本文算法得到的地震震中可以直接给定为2类精度。由于震中位置主要受烈度点分布空间位置的限制，约束性很高，因此所得估算结果很好。

图6 外推检验地震参数估计的综合不确定度

对于震级综合不确定度，随着烈度点数的增多，一部分地震的综合不确定度越来越小，一部分地震的综合不确定度先减小后增大，另一部分地震的综合不确定度则越来越大。综合不确定度包含了随机误差和系统误差。随着烈度点数的增多随机误差越来越小，因此影响综合不确定度大小的主要是系统误差。由于所建模型是从很多地震里提取的平均规律，具有离散性，与任何一个真实地震之间都存在着系统偏差。由图4（b）、4（d）知道地震震级的估算随着烈度点数的增多有一个先增大后平稳的过程。当用最少的烈度点计算出来的震级值都比仪器震级值大时，随着烈度点数的增多，系统偏差会越来越大，导致综合不确定度随着烈度点数的增加而增大；当用最多的烈度点计算出来的震级值都比仪器震级值小时，随着烈度点数的增多，计算震级会越来越靠近仪器震级，显示为综合不确定度随着烈度点数的增多而减小；而用最少烈度点计算出来的震级比仪器震级小，最多烈度点计算出来的震级比仪器震级大时，系统偏差有一个先减小后增大的过程，就导致了综合不确定度随着烈度点数的增多先减小后增大。

造成震级不确定度无法控制的主要原因有：①所建模型时根据的是一个地区众多地震中的平均规律，模型回归需要有足够的地震样本，因此选定的研究区域比较大；而研究区域越大，则场地条件等因素就越偏离某个具体地震的实际情况。为了使模型更能适应具体地震情况，就需要缩小研究区域，但代价是研究区域内发生的地震太少，造成模型回归的离散性变大。两相钳制，较难取舍。②由烈度推算震级，震源深度等参数对结果影响较大。本文所建模型里暂时没有考虑震源深度等其他因素的影响，但每个具体地震的深度等参数仍然是有差别的，造成具体地震对所建模型而言各有偏差。

对地震震级的估算历来都是一个难题，本文算法能同时估算出地震的震级和震中。由综合不确定度来看，震中的精度可以定为2类，但是震级的精度很难直接给出。对于参与验算的地震，随着烈度点数的增多，震级值趋于平稳时，不确定度都不超过0.8级，多数在0.5级以下

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

3 历史大震参数估计

通过前面的分析，本文建立了一套基于烈度数据点考虑断层破裂长度的烈度分布椭圆模型以估算地震震中与震级的方法，我们以此来估算华北地区几个6.5～8.0级历史大震，结果见表6。历史大震的烈度数据点均从《中国历史地震图集》（明、清时期）（国家地震局地球物理研究所等，1986、1990）中提取。由表6可以看到，计算震级与图集震级相差均在1级以内，计算震中与图集震中相差不超过25km，对应于历史地震2类震中精度。对于历史大震，这样的计算结果比较可靠。

表6 华北地区历史大震参数估算结果

本文用5个华北地区历史强震证明了所建算法的可行性。由以上震例可见，此方法对历史强震参数的估定颇为有效。处理过程直接明了，减少了等震线勾画过程中的主观不确定性。

4 结论与讨论

本文收集现代大震资料，统计回归了断层破裂长度-震级经验关系；针对6.5级以上强震烈度分布特征，提出了基于烈度数据点的考虑断层破裂长度的烈度分布椭圆模型，在此基础上联立椭圆数学方程建立了强震参数估计方法，通过內符检验和外推检验证明了算法的可行性。采用蒙特卡洛方法，通过大量的抽样模拟，定量分析了所得参数的不确定性，结果表明：计算震级精度在1级以内，计算震中精度为历史地震的2类精度。进而选取了华北地区5个历史强震进行了地震参数估算。本方法直接通过烈度数据点来估算历史强震的震级与震中，回避了对等震线的勾画，处理过程直接明了，一定程度上减少了主观不确定性，提高了科学性。由于在建模过程中，8.0级以上特大地震的资料太少，因此本方法只适用于6.5～8.0级的历史强震，对于8.0级以上的特大地震需要另行讨论。

不论是破裂长度与震级的经验关系，还是强震烈度分布模型都应该具有区域性，应该根据研究区域的地震相关资料统计得到。但在一个研究区内很难找到足够多的现代大震统计资料，所以不得不采用全国的强震数据建立经验公式，那么应用到某个研究区时必然会产生某些偏差。地震烈度分布模型是据多个地震的平均规律而建立的，与地震资料之间存在着离散性，与任何一个真实地震之间都有系统偏差。模型是原型的抽象和简化表征，建模本身是一个持续的永无止境的活动集合。随着技术的发展、科学的进步以及基础数据的积累，需要对模型不断的改进以期更接近真实情况。

笔者认为，今后应多引入现代地震学的方法来处理历史地震，例如早期国外台站的波形数据、震源机制等，以提高历史地震参数的可靠性

致谢：衷心感谢金严研究员、汪素云研究员、俞言祥研究员和潘华研究员在本文研究过程中给出的有益指导与建议！