2.1. 强震烈度估计等值线拟合

我国目前公开发表的烈度调查资料大多采用等震线的形式发布，原始烈度调查点通常难以获得。本文选取已经出版的《中国震例（1966-2002）》[2–10]，2008年后的强震烈度资料从中国地震局官方网站的地震专题上获取(http://www.cea.gov.cn/publish/dizhenj/468/553/index.html),收集了20世纪60年代以来中国西部地区所发生的震级大于等于6.5级地震的等震线图，提取了其中既有仪器测量数据又有宏观烈度调查数据的现代地震资料(表1)，并将烈度分布图数字化。这些地震的等震线图上部分明确标示了宏观烈度调查点位置；对于未明示烈度调查点的，本文将能在等震线图上明确获取烈度信息的城、县、镇、村所在地作为烈度点。经数字化配准后获取各地震烈度点的空间分布。

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

式中，（x_i,y_i）为椭圆上点的坐标，(x₀, y₀)为椭圆中心坐标；θ(0° ≤ θ ≤ 180°)是椭圆长轴逆时针与x轴正方向的夹角；R_a、R_b分别为椭圆长半轴和短半轴长度，t_i(0 ≤ θ ≤ 2π)为椭圆上各点所对应的参数。

吴清等[1]根据我国现代大震资料统计拟合了地震破裂长度-震级关系式

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

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

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

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

图1

2013年7月22日甘肃岷县6.6级地震烈度估计等值线分布图.

2.2. 中国西部地区强震烈度分布模型

本文据陈达生等[11]关于地震烈度椭圆衰减关系长短轴统一回归的思想，采用如下烈度分布模型：

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

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

由原始烈度点拟合得到烈度区的估计等值线，进而统计回归得到地震烈度分布模型，其中I_a为长轴方向，I_b为短轴方向：

3.1. 强震参数估计方法

吴清等[1]给出了强震参数估计方法：

式中，C_1a、C_1b、C₂、C_3a、C_3b、R_0a、R_0b由烈度分布模型统计回归可得，见式（5）（6）。(x_i, y_i, I_i)为烈度数据点，(x₀, y₀)为震中，M为震级，θ为方向。若已知一个地震的多个烈度信息点(x_i, y_i, I_i)，带入到（7）式，即可联立方程组求得震中(x₀, y₀)、震级M和方向θ。有4个未知数，则至少需要4个方程确定一组解，但x₀与y₀不是相互独立的，因此在已知烈度分布模型下，3个烈度信息点就能确定一组震中、震级和方向。烈度信息点增多，变成求解超定非线性方程组，采用通用全局优化法的数值计算方法寻求最优解，可非常直观的求得地震宏观震中位置、震级和方向。

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

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

中国大陆的地震绝大部分发生在地壳以内，其震源深度差别不大[13]，因此本文暂不考虑震源深度对烈度分布的影响。

3.2.1. 內符检验

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

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

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

3.2.2. 外推检验

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

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

与吴清等[1]得到的适用于全国的烈度分布模型相比，由适用于西部地区的烈度分布模型估算得来的地震参数较为准确。