Measurement method for isoseismal major and minor axes based on ellipse fitting algorithm
-
摘要:
提出了一种改进的拟合椭圆算法用于测量等震线长短轴，即将等震线测量的模式和直接测量算法的概念以约束的形式加入拟合椭圆算法，并使用最小二乘法求解椭圆系数，在与传统等震线测量方法保持一致的基础上，减少了直接测量方法的主观不确定性。在此基础上讨论了拟合椭圆算法的适用性和鲁棒性，将等震线长短轴的测量结果与前人的结果进行了对比。结果表明：拟合椭圆算法的适用性较好，可用于等震线的拟合；拟合椭圆算法的鲁棒性较好，计算结果较为稳定；对于未闭合且形状复杂的等震线，计算结果相对离散，要单独验证结果是否合理；利用拟合椭圆算法得到的等震线长短轴结果与前人的结果较为一致，可用于建立烈度衰减关系。
-
关键词:
- 等震线长短轴测量方法 /
- 拟合椭圆算法 /
- 最小二乘拟合 /
- 算法测试
Abstract:As the fundamental data for establishing the seismic intensity attenuation relationship, the reliability of isoseismal data directly affects the accuracy and reliability of the intensity attenuation relationship. The acquisition of traditional isoseismal parameters is based on direct measurement, which inevitably introduces a degree of subjectivity. A more objective approach to measuring the parameters of isoseismals has been proposed. However, it is inevitable that conceptual differences will exist due to algorithmic differences and the lack of links between the novel algorithm and commonly used direct measurement methods. In view of this, this paper introduces a novel approach to measuring the major and minor axes of isoseismals. The mode of isoseismal measurement and the concepts from the direct measurement method are incorporated into the ellipse fitting algorithm as mathematical constraints, thereby providing a systematic method for measuring the radii of isoseismal major and minor axes. In order to reduce the subjective influence of isoseismal major and minor-axis measurements, the radii of isoseismal major and minor axes are measured using the best-fitting ellipse of the isoseismal. This approach is based on the consistency and inheritance of the methods and concepts of the direct measurements. The resulting data provide the basic data for the establishment of intensity attenuation relationship.
The measurement modes of isoseismals are classified into two types according to whether the major-axis strike is fixed or not. In order to accommodate this distinction, the strike of major axis is also incorporated into the algorithm as a constraint. Additionally, the area of isoseismal is incorporated into the algorithm as a constraint. This leads to the proposal of four ellipse fitting algorithms in total: ① An unconstrained fitting algorithm, which implies that no additional constraints are imposed during the model fitting process. ② A fitting algorithm with major-axis strike constrained, which involves the addition of a priori major-axis strike information to the model fitting process. ③ A fitting algorithm with area and major-axis strike constrained, which incorporates both a priori major-axis strike and area into the model fitting process. ④ A fitting algorithm with area constrained, which entails the addition of a priori model area to the model fitting process.
It should be noted that the results of the algorithm will degrade from an ellipse to a circle when the isoseismal shape is circular （i.e. a special ellipse）. Therefore, it is necessary to construct circular isoseismal and simulate different levels of sampling noises so as to test the applicability of the algorithm. The results of the algorithm test demonstrate that all four algorithms yield satisfactory outcomes when the noise is minimal. However, the algorithm with area constraints exhibits a significantly superior performance compared with the algorithm without area constraints when the noise is elevated. Nevertheless, the applicability of the algorithms for degradation to a circle is deemed acceptable, given that extreme noise cases are unlikely to occur in practical scenarios.
The isoseismal data used by the ellipse fitting algorithm is subject to sampling bias, which may affect the output of algorithm. Therefore, it is necessary to analyze the impact of data sampling differences on the results of algorithm. The robustness of the algorithm was evaluated by randomly sampling points of isoseismal to simulate different sampling scenarios. The results indicate that the parameters show normally distributed, with the majority falling within one standard deviation of the mean value. Furthermore, the algorithm is robust. For the unclosed isoseismal with complex shapes, the calculation results are relatively discrete, and thus require separate verification. Furthermore, the measurements obtained by the ellipse fitting algorithm are consistent with previous results and can be used to establish the intensity attenuation relationship.
-
引言
烈度衰减关系在地震科学研究、地震灾害评估、工程抗震等方面发挥着重要作用，而等震线资料作为建立烈度衰减关系的基础数据，其可靠性和精度直接影响着烈度衰减关系的准确性和可靠性（胡聿贤，2007）。
国外学者一般采用圆模型建立烈度衰减关系（Howell，Schultz，1975；Chandra，1979；Musson，2005；Sørensen et al，2010；Prieto et al，2011；Allen et al，2012；Ahmadzadeh et al，2020），等震线图或烈度点是其使用的基础数据，其中等震线图测量算法大致可以分为三类：① 测量不同走向的等震线半径。Anderson （1978）从最内圈等震线中心绘制16条射线，测量每条射线与每条等震线相交的距离，推导出美国修正麦卡利烈度（modified Mercalli intensity，缩写为MMI）烈度衰减关系。Ghosh和Mahajan （2013）根据Anderson （1978）的方法，搜集了计算喜马拉雅西北部10次地震（MS≥4.9）的等震线资料，结合不同方向不同震级地震的衰减模式建立了两种烈度衰减关系。② 基于等震线面积计算等效圆的半径，即根据每条等震线包围的面积，使用等效圆的方法得到每条等震线的距离参数，并将其用于推导区域烈度衰减关系（Howell，Schultz，1975；Chandra，1979；Casado et al，2000；Ahmadzadeh et al，2020）。③ 基于等效圆半径和不同走向平均半径的混合测量方法。Ambraseys和Douglas （2004）重新评估了印度北部地区的烈度和等震线分布，完善了等震线等效半径的测量方法，得到开放等震线的等效半径，并导出利用等震线半径估计面波震级MS的公式。Prieto等（2011）使用Ambraseys和Douglas （2004）的方法，收集并计算了1766—2004年哥伦比亚及其邻区68次地震事件的等震线图数据，给出了该区域的烈度衰减概率分布。
等震线的长短轴半径是椭圆烈度衰减模型的基础数据，基于该模型的等震线长短轴测量方法大致可以分为三类：① 基于等震线图资料的直接测量方法。陈达生和刘汉兴（1989）基于等震线图的比例信息给出了长短轴测量的常用算法：直观法和等面积法，并由此建立了华北地区的烈度衰减关系。直观法作为传统方法被用来测量等震线长短轴半径及建立中国分区烈度衰减关系（汪素云等，2000；雷建成等，2007；肖亮，俞言祥，2011；俞言祥等，2013；Xu et al，2021）。② 基于GIS平台的数字化等震线测量方法。李永强等（2006）借助GIS平台制作云南及其邻区175次强震的等震线图，并测量等震线的长短轴和面积；孙继浩（2011）借助ArcGIS平台的椭圆功能近似等震线，之后测量近似椭圆的长短轴半径，并对川滇及其邻区中强地震烈度衰减关系的适用性进行了研究。③ 基于拟合椭圆测量等震线长短轴的方法。郁曙君（1993）提出使用椭圆投影法求取等震线的长短轴半径，利用正交回归建立川滇地区震级分档烈度衰减关系；沙海军等（2010）提出利用拟合椭圆算法测量等震线长短轴，该算法利用椭圆内一点到椭圆上每点的距离与椭圆长轴和短轴半径的关系建立方程，使用Matlab的nlinfit函数进行求解，长短轴半径的计算结果与前人测量结果一致；吴清（2013）使用拟合椭圆算法得到基于烈度点的椭圆烈度估计线，并基于此建立椭圆烈度分布模型。
等震线的直接测量方法操作简单易于实现，但无法避免主观因素的影响。汪素云等（2000）指出其使用的等震线长短轴半径与第三代区划图的资料相比有所变动，研究者们也意识到这一问题，逐渐采取更加客观的方式测量等震线长短轴半径。国家地震局（1996）分析我国地震烈度等震线特征时曾指出，烈度等震线在近场通常呈椭圆形，随着距离的增大，在远场逐渐变为圆形，因此一般选择椭圆对等震线进行建模（郁曙君，1993；沙海军等，2010；吴清，2013）。然而，由于所使用的拟合椭圆算法并不相同，且未与等震线测量的模式和常用的直接测量方法建立联系，因此势必存在数据定义不一致的问题。鉴于此，本文尝试将等震线测量的模式以及直接测量方法中的概念作为数学约束加入到拟合椭圆算法，系统地给出基于拟合椭圆测量等震线长短轴半径的方法，旨在与传统等震线测量的模式和概念保持一致的基础上，利用等震线的最优拟合椭圆测量等震线长短轴半径，以期减少等震线长短轴测量时的主观影响，进而为建立分区烈度衰减关系提供基础数据。
1. 等震线测量理论基础
通过椭圆参数${\boldsymbol{\beta}} $＝$ ［ $a， b， θ， x0， y0$ ］ $可以确定空间位置中的唯一椭圆，其中a为长轴半径，b为短轴半径，θ为长轴走向（图1），（x0，y0）为中心坐标。
![]()
图 1 椭圆长轴走向示意图（a） 走向范围为0°—90°；（b） 走向范围为−90°—0°Figure 1. Diagram of the major-axis strike of the ellipse（a） Strike ranges 0°−90°；（b） Strike ranges −90°−0°$$ \frac{{{x^2}}}{{{a^2}}} ＋ \frac{{{y^2}}}{{{b^2}}} ＝ 1 $$ (1) 为标准椭圆方程，该方程经过旋转、平移得到一般椭圆方程如下：
$$ \begin{split} f （ {{\boldsymbol{\tau}} }\text{，} {{\boldsymbol{x}}} ） ＝& {\boldsymbol{\tau}} \cdot {\boldsymbol{x}} ＝ A{x^2} ＋ B{y^2} ＋ Cxy ＋\\& Dx ＋ Ey ＋ F ＝ 0 ， \end{split} $$ (2) 式中τ＝$ ［ $A， B， C， D， E， F $ ］ $T，x＝$ ［ $x2， y2， xy， x， y， 1$ ］ $T，
$$ \begin{split} \left\{ \begin{array}{l} A ＝ {a^2}{\sin ^2}\theta ＋ {b^2}{\cos ^2}\theta ，\\ B ＝ {b^2}{\sin ^2}\theta ＋ {a^2}{\cos ^2}\theta ，\\ C ＝ 2 （ {a^2} － {b^2} ） \sin \theta \cos \theta ，\\ D ＝ - 2A{x_0} － C{y_0}，\\ E ＝ - C{x_0} － 2B{y_0}，\\ F ＝ Ax_0^2 ＋ C{x_0}{y_0} ＋ By_0^2 － {a^2}{b^2}. \end{array} \right. \end{split}$$ (3) 通过式（3）由椭圆参数可得到一般椭圆方程的系数。在此基础上，通过
$$ \left\{ \begin{gathered} {a^2} ＝ \frac{{2 （ Ax_0^2 ＋ By_0^2 ＋ C{x_0}{y_0} － F ） }}{{A ＋ B ＋ \sqrt {{{ （ A － B ） }^2} ＋ {C^2}} }} \\ {b^2} ＝ \frac{{2 （ Ax_0^2 ＋ By_0^2 ＋ C{x_0}{y_0} － F ） }}{{A ＋ B － \sqrt {{{ （ A － B ） }^2} ＋ {C^2}} }} \\ \frac{{2\tan \theta }}{{{{\tan }^2}\theta － 1}} ＝ \frac{C}{{A － B}} \\ {x_0} ＝ \frac{{CE － 2BD}}{{4AB － {C^2}}} \\ {y_0} ＝ \frac{{CD － 2AE}}{{4AB － {C^2}}} \\ \end{gathered} \right. $$ (4) 即可得到椭圆参数${\boldsymbol{\beta }} $。因此使用拟合椭圆测量等震线的整体思路如下：首先使用椭圆拟合算法得到等震线的最佳拟合椭圆，之后利用式（4）计算得到等震线的长短轴半径。
2. 拟合椭圆算法
在模式识别和计算机视觉中，经常需要拟合或检测椭圆形物体或特征（Rosin，1993；Fitzgibbon et al，1999），其重点在于拟合散乱数据中的椭圆并避免平凡解，这与本文的需求和数据情况并不相同。因此，本文仅借鉴其拟合思想针对等震线椭圆拟合的具体需求对算法进行调整。
本文拟合椭圆算法使用的数据为等震线在投影坐标系下的点坐标，拟合时进行去均值处理，这样做的优点是：① 直接使用平面坐标系，避免地球曲率变换计算带来的误差；② 公里网坐标和面积可直接代入运算；③ 去均值处理提高计算的数值稳定性，且不影响计算结果。
$f （ {{{\boldsymbol{x}}}}_{i} ） $为点（xi， yi）到椭圆方程${{f}} （ {\boldsymbol{\tau}} \text{，} {\boldsymbol{x}} ） $的代数距离，借用残差平方和的概念定义代数距离平方和$D （ {\boldsymbol{\tau}} ） $如下：
$$ D （ {\boldsymbol{\tau}} ） ＝ \sum\limits_{i ＝ 1}^N {f^2{{ （ {{{x}_i}} ） }}} . $$ (5) 式（5）表示椭圆对数据点的拟合程度，其值越小，对数据的拟合能力越好。使用最小二乘方法求解即得到拟合椭圆系数${\boldsymbol{\tau }} $，进而利用式（4）计算椭圆参数${\boldsymbol{\beta}} $。
定义
$$ \rho ＝ \frac{{D （ {{{{\boldsymbol{\tau}} }_{{\mathrm{ALGO1}}}}} ） }}{{D （ {{{{\boldsymbol{\tau}} }_{{\mathrm{ALGO2}}}}} ） }} $$ (6) 比较不同算法拟合结果对等震线的近似程度，ρ＞1表明算法1 （ALGO1）结果的近似程度较差，ρ＝1表明两种算法结果的近似程度相当，ρ＜1表明算法2 （ALGO2）结果的近似程度较差。
等震线的测量模式分为长轴转向和长轴不转向两种（陈达生，刘汉生，1989），分别定义为：① 长轴转向，即以每条等震线最长方向的半径为长轴半径，与其相垂直的半径为短轴半径，长轴方向可能改变； ② 长轴不转向，即以最内圈等震线的长轴方向的半径为长轴半径，以与其相垂直的半径为短轴半径，长轴方向固定不变。两种测量模式结果的示意图如图2所示， 《中国地震烈度区划图（1990）概论》（国家地震局，1996）表明，两种测量模式的数据在地震烈度衰减关系的统计中，拟合结果相差不大，因此多采用长轴转向的模式。直观法和等面积法是传统的直接测量方法（国家地震局，1996），其中等面积法能够在一定程度上减少等震线测量的主观不确定性。因此，将等震线测量模式的定义和传统测量方法的思想以约束的形式添加到椭圆拟合方法中，得到四种椭圆拟合算法：① 无约束拟合算法；② 约束长轴走向的拟合算法；③ 约束面积与长轴走向的拟合算法；④ 约束面积的拟合算法。通过上述四种算法得到两种等震线测量模式的具体流程如下：
![]()
图 2 长轴转向（左）和不转向（右）模式的地震等震线图拟合结果（a，b） 1976年7月28日唐山MS7.8地震；（c，d） 2008年5月12日汶川MS8.0地震Figure 2. Fitting results of the isoseismals of the earthquakes with major-axis steering （a） and non-steering （b） modes（a，b） Tangshan MS7.8 earthquake on July 28，1976； （c，d） Wenchuan MS8.0 earthquake on May 12，20081） 长轴转向模式。如果是闭合等震线，使用约束面积的拟合算法计算，否则使用无约束拟合算法。
2） 长轴不转向模式。首先计算最内圈等震线，获得最内圈等震线长轴走向（如果是闭合等震线，使用约束面积的拟合算法计算，否则使用无约束拟合算法计算）；根据上一步计算得到的最内圈等震线长轴走向，如果等震线闭合，使用约束面积和长轴走向的拟合算法计算，否则使用约束长轴走向的拟合算法。
图2以1976年7月28日唐山MS7.8地震等震线图和2008年5月12日汶川MS8.0地震等震线图为例展示不同等震线测量模式下的结果。
我们选择了四类具有代表性的等震线：简单闭合等震线（以2013年4月20日芦山MS7.0地震Ⅶ度等震线为例，图3a）、复杂闭合等震线（以1945年9月23日滦县MS6.3地震Ⅶ度等震线为例，图3b）、简单未闭合等震线（以1986年1月28阳江MS4.8地震Ⅳ度等震线为例，图3c）、复杂未闭合等震线（以1976年7月28唐山MS7.8地震Ⅵ度等震线为例，图3d）。
![]()
图 3 代表性的等震线（a） 简单闭合等震线：2013年芦山MS7.0地震Ⅶ度等震线；（b） 复杂闭合等震线：1945年滦县MS6.3地震Ⅶ度等震线；（c） 简单未闭合等震线：1986年阳江MS4.8地震Ⅳ度等震线；（d） 复杂未闭合等震线：1976年唐山MS7.8地震Ⅵ度等震线Figure 3. Representative isoseismals（a） Closed isoseismal with simple shape：Isoseismal Ⅶ of 2013 Lushan MS7.0 earthquake；（b） Closed isoseismal with complex shape：Isoseismal Ⅶ of 1945 Luanxian MS6.3 earthquake；（c） Unclosed isoseismal with simple shape：Isoseismal Ⅳ of 1986 Yangjiang MS4.8 earthquake；（d） Unclosed isoseismal with complex shape：Isoseismal Ⅵ of 1976 Tangshan MS7.8 earthquake2.1 无约束拟合算法
无约束拟合算法是在拟合时不施加约束得到等震线的最优拟合椭圆，可视为等震线长轴转向模式下直接测量方法中的直观法在拟合椭圆算法中的延伸。Rosin （1993）提出按F＝1的方式将椭圆方程式（2）进行归一化，基于最小二乘法的思想求解式（5）得到拟合椭圆，经过验证得到结论：数据的平移和旋转对拟合结果几乎无影响。Fitzgibbon等（1999）提出椭圆拟合的直接最小二乘法，将拟合结果限定在椭圆而非一般二次曲线，通过广义本征系统求解，其结果表明该算法对噪声具有高鲁棒性，且计算效率较高。
以Rosin （1993）算法为算法1、Fitzgibbon等（1999）算法为算法2 $ ［ $式（6）$ ］ $，比较四类等震线在两种算法下的拟合结果，如图4所示。可见：对于简单的等震线（图4a，c），两种算法拟合结果的差异较小，所得椭圆参数基本相同；对于复杂的等震线（图4b，d），两种算法拟合结果差别较大，其中未闭合等震线的拟合结果差异最大；两种算法下四类等震线的拟合结果显示，Fitzgibbon等（1999）算法更适合等震线的椭圆拟合。因此，选择Fitzgibbon等（1999）算法作为无约束拟合算法。
![]()
图 4 四类等震线在不同无约束算法下的拟合结果的对比（a） 简单闭合等震线：2013年芦山MS7.0地震Ⅶ度等震线；（b） 复杂闭合等震线：1945年滦县MS6.3地震Ⅶ度等震线；（c） 简单未闭合等震线：1986年阳江MS4.8地震Ⅳ度等震线；（d） 复杂未闭合等震线：1976年唐山MS7.8地震Ⅵ度等震线Figure 4. Comparison of the fitting results of four isoseismals by the different unconstrained algorithms（a） Closed isoseismal with simple shape：Isoseismal Ⅶ of 2013 Lushan MS7.0 earthquake；（b） Closed isoseismal with complex shape：Isoseismal Ⅶ of 1945 Luanxian MS6.3 earthquake；（c） Unclosed isoseismal with simple shape：Isoseismal Ⅳ of 1986 Yangjiang MS4.8 earthquake；（d） Unclosed isoseismal with complex shape：Isoseismal Ⅵ of 1976 Tangshan MS7.8 earthquake2.2 约束长轴走向的拟合算法
等震线长轴不转向模式要求测量每条等震线在最内圈等震线长轴走向上的长轴半径，因此将最内圈等震线长轴走向作为约束条件，得到等震线在约束长轴走向下的最优拟合椭圆。利用椭圆系数-参数方程$ ［ $式（4）$ ］ $对椭圆方程$ ［ $式（2）$ ］ $的系数进行约束，回归得到约束长轴走向的拟合椭圆。约束长轴走向的拟合算法如下：
给定长轴走向θ*，由式（4）定义
$$ \lambda ＝ \frac{{2\tan {\theta _*}}}{{{{\tan }^2}{\theta _*} － 1}} ＝ \frac{C}{{A － B}} ，$$ (7) 从而得到椭圆方程的系数关系$C＝\lambda （ A－B ） $，将其带入椭圆方程$ ［ $式（2）$ ］ $，有
$$ A （ {{x^2} ＋ \lambda xy} ） ＋ B （ {{y^2} － \lambda xy} ） ＋ Dx ＋ Ey ＋ F ＝ 0. $$ (8) 根据Rosin （1993）算法思想，对式（8）进行F＝1归一化处理，回归即得到约束长轴走向的椭圆方程。
根据长轴不转向模式的定义，选取最内圈等震线的长轴走向为约束。以约束长轴走向的拟合算法为算法1、无约束拟合算法为算法2 $ ［ $式（6）$ ］ $，比较四类等震线长轴走向约束和无约束的结果，如图5所示。可见：约束长轴走向的拟合算法能够较好地按照设置的长轴走向进行拟合，得到约束长轴走向的最佳拟合椭圆；添加长轴走向作为约束条件后，椭圆结果的近似程度较差，这是将最内圈等震线长轴走向设置为目标走向所致。
![]()
图 5 四类等震线在无约束拟合算法和约束长轴走向拟合算法下拟合结果的对比（a） 简单闭合等震线：2013年芦山MS7.0地震Ⅶ度等震线；（b） 复杂闭合等震线：1945年滦县MS6.3地震Ⅶ度等震线；（c） 简单未闭合等震线：1986年阳江MS4.8地震Ⅳ度等震线；（d） 复杂未闭合等震线：1976年唐山MS7.8地震Ⅵ度等震线Figure 5. Comparison of the fitting results of the isoseismals by the unconstrained algorithm with those by the algorithm with major-axis strike constrained（a） Closed isoseismal with simple shape：Isoseismal Ⅶ of 2013 Lushan MS7.0 earthquake；（b） Closed isoseismal with complex shape：Isoseismal Ⅶ of 1945 Luanxian MS6.3 earthquake；（c） Unclosed isoseismal with simple shape：Isoseismal Ⅳ of 1986 Yangjiang MS4.8 earthquake；（d） Unclosed isoseismal with complex shape：Isoseismal Ⅵ of 1976 Tangshan MS7.8 earthquake2.3 约束面积和长轴走向的拟合算法
直接测量方法中的等面积算法使用等震线包围的面积S，测量长轴半径a后，利用椭圆面积公式
$$ S ＝ \text{π} ab $$ (9) 计算短轴半径b，等面积法一定程度上减少了测量时的主观因素。基于此种思想，本文尝试将长轴走向约束和面积约束带入椭圆方程$ ［ $式（2）$ ］ $，由此得到约束面积与长轴走向的等震线最佳逼近椭圆。将面积约束带入椭圆系数方程
$$ \left\{\begin{array}{l} A＝a^2\sin^2\theta＋\left(\dfrac{S}{\text{π}a}\right)^2\cos^2\theta\text{，} \\ B＝\left(\dfrac{S}{\text{π}a}\right)^2\sin^2\theta＋a^2\cos^2\theta\text{，} \\ C＝2\left[a^2－\left(\dfrac{S}{\text{π}a}\right)^2\right]\sin\theta\cos\theta\text{，} \\ D＝－2Ax_0－Cy_0\text{，} \\ E＝－Cx_0－2By_0\text{，} \\ F＝Ax_0^2＋Cx_0y_0＋By_0^2－\left(\dfrac{S}{\text{π}a}\right)^2，\end{array}\right. $$ (10) 由于椭圆系数${\boldsymbol{\tau}} $的耦合，参数a和b难以通过简单的替换实现2.2节的效果。确定长轴半径a ＋，等震线面积S和长轴走向θ＋后，可以确定系数A＋，B＋和C＋$ ［ $式（10）$ ］ $，未知系数即D，E和F，对应未知参数即x0和y0。D，E和F可按
$$ Dx＋Ey＋F＝－ （ A_＋x^2＋B_＋y^2＋C_＋xy ） $$ (11) 回归得到，进而计算参数x0和y0。因此，问题可以转换为：约束长轴走向θ ＋ 后，通过一系列长轴半径预设值a ＋ ，得到椭圆方程$ ［ $式（11）$ ］ $，利用代数距离平方和$D （ {\boldsymbol{\tau}} ） $判断拟合椭圆对等震线的近似程度，$D （ {\boldsymbol{\tau}} ） $越小，近似程度越好。
综上可知，通过面积S可得到一系列长轴半径的预设值a+ （根据长短轴定义a≥b），选择使$D （ {\boldsymbol{\tau}} ） $最小的情况作为最优估计值，这种遍历算法记为遍历算法1。根据李建文等（2010）提出的自适应极值搜索算法对遍历算法进行优化，优化后的算法记为遍历算法2。
比较遍历算法1与遍历算法2，结果如图6所示，其中遍历算法1对长轴半径的遍历范围为1—5倍的等效圆半径（步长设置为0.1 km），遍历算法2设置四个迭代起始点（图6红色空心三角形）。可见：两种算法得到的参数最优解相同（图6红色方框内红色空心圆），遍历算法2的精度更高；遍历算法2不依赖于迭代初始值，结果的稳定性较好（图6红色方框内红色空心圆）；遍历算法2计算量较小，计算速度较快。因此选择遍历算法2作为约束面积与长轴走向的拟合算法。
![]()
图 6 遍历算法1和遍历算法2计算过程及结果对比（a） 2013年4月20日芦山MS7.0地震Ⅶ度等震线；（b） 1945年9月23日滦县MS6.3地震Ⅶ度等震线Figure 6. Comparison of calculation process and result of two traversal algorithms（a） Isoseismal Ⅶ of Lushan MS7.0 earthquake on April 20，2013；（b） Isoseismal Ⅶ of Luanxian MS6.3 earthquake on September 23，1945以约束面积与长轴走向的拟合算法为算法1、约束长轴走向的拟合算法为算法2$ ［ $式（6）$ ］ $，比较2013年4月20日芦山MS7.0地震和1945年9月23日滦县MS6.3地震的Ⅶ度等震线在面积与长轴走向约束和长轴走向约束下的结果，如图7所示。可见：约束面积与长轴走向的拟合算法能够较好地完成约束条件下的拟合；对于简单闭合等震线，两种算法结果对等震线的逼近程度相近，约束面积与长轴走向的拟合算法对复杂闭合等震线的近似程度相对较好。总体而言，增加面积约束可提高算法拟合结果对等震线的近似程度。
![]()
图 7 约束长轴走向的拟合算法和约束面积与长轴走向的拟合算法所得拟合结果的对比（a） 2013年4月20日芦山MS7.0地震Ⅶ度等震线；（b） 1945年9月23日滦县MS6.3地震Ⅶ度等震线Figure 7. Comparison of fitting results obtained by the algorithm with major-axis strike constrained with those by the algorithm with area and major-axis strike constrained（a） Isoseismal Ⅶ of Lushan MS7.0 earthquake on April 20，2013；（b） Isoseismal Ⅶ of Luanxian MS6.3 earthquake on September 23，19452.4 约束面积的拟合算法
约束面积与长轴走向的拟合算法是等面积法在长轴不转向模式下的扩展算法，是对等面积法在长轴转向模式下的扩充，约束面积的拟合算法不约束长轴走向。本文首先确定等震线的最佳长轴走向，之后使用约束面积与长轴走向的拟合算法得到面积约束下的等震线最佳逼近椭圆。
以约束面积的拟合算法为算法1、无约束拟合算法为算法2 $ ［ $式（6）$ ］ $，比较简单闭合等震线和复杂闭合等震线在面积约束和无约束下的结果，如图8所示。可见：约束面积的拟合算法能够较好地完成约束条件下的拟合；对于简单闭合等震线，不同算法对等震线的近似程度相当，而对于复杂闭合等震线，约束面积的拟合算法对等震线的近似程度相对较好。总体而言，增加面积约束可提高算法拟合结果对等震线的近似程度。
![]()
图 8 无约束拟合算法和约束面积的拟合算法所得拟合结果的对比（a） 2013年4月20日芦山MS7.0地震Ⅶ度等震线；（b） 1945年9月23日滦县MS6.3地震Ⅶ度等震线Figure 8. Comparison of fitting results obtained by the unconstrained algorithm with those by the algorithm with area constrained（a） Isoseismal Ⅶ of Lushan MS7.0 earthquake on April 20，2013；（b） Isoseismal Ⅶ of Luanxian MS6.3 earthquake on September 23，19453. 讨论
3.1 适用性测试
前文结果已经说明本文算法能够藉由拟合等震线得到最佳椭圆。从算法适用性角度，当等震线形状趋向于圆（特殊的椭圆）时，算法结果应从椭圆退化为圆。因此，对于构造半径R为200 km的圆等震线，用高斯噪声$ ［ $d～N （0， σ2），其中d为采样噪声，N为正态分布，σ为采样误差水平$ ］ $模拟等震线的取样误差$ ［ $式（12）$ ］ $，分别选取σ＝20和σ＝100用于检验算法适用性（图9，图10）。
![]()
图 9 误差水平σ＝20时四种算法的适用性检验（a） 无约束拟合算法；（b） 约束长轴走向的拟合算法；（c） 约束面积与长轴走向的拟合算法；（d） 约束面积的拟合算法Figure 9. Applicability test of four algorithms when error σ is taken as 20（a） Unconstrained fitting algorithm；（b） Fitting algorithm with major-axis strike constrained；（c） Fitting algorithm with area and major-axis strike constrained；（d） Fitting algorithm with area constrained![]()
图 10 误差水平σ＝100时四种算法的适用性检验（a） 无约束拟合算法；（b） 约束长轴走向的拟合算法；（c） 约束面积与长轴走向的拟合算法；（d） 约束面积的拟合算法Figure 10. Applicability test of four algorithms when error σ is taken as 100（a） Unconstrained fitting algorithm；（b） Fitting algorithm with major-axis strike constrained；（c） Fitting algorithm with area and major-axis strike constrained；（d） Fitting algorithm with area constrained$$ {x_i} ＝ R\cos {\theta _i} ＋ {d_i} ， {y_i} ＝ R\sin {\theta _i} ＋ {d_i}， $$ (12) 式中：xi和yi为模拟得到的边界取样点的横纵坐标，R为构造的圆等震线半径，参数θ的取值范围为0°—360°，步长为1°，di为噪声模拟取样误差。
σ反映了取样的误差水平，σ＝20相对R＝200 km较低，仅能代表等震线数字化时的一般情况。图9的结果表明，当取样误差水平较低（σ＝20）时，四种算法都能得到较为理想的结果。误差水平σ＝100相对R＝200 km较高，是一种十分极端的情况（如图10所示，取样点相对等震线较为离散）。由图10可见，四种算法都能根据数据情况得到拟合结果，其中：约束面积与长轴走向的拟合算法和约束面积的拟合算法得到的结果较为理想，分析是源于施加的面积约束效果较强；无约束拟合算法和约束长轴走向的拟合算法虽然也给出了相应的拟合结果，但由于数据点过于离散，拟合效果较差，这种极端情况在实际等震线数字化中几乎不会出现，因此认为该算法对于退化为圆的适用性可以接受。
3.2 鲁棒性测试
本文通过椭圆拟合算法得到最佳拟合椭圆并求解椭圆参数，但算法使用的等震线数据在取样上具有一定的主观因素，因此我们需要对算法的鲁棒性进行检验，以此分析数据取样的主观性对算法的影响。图9和图10的对比已经说明，对于理想等震线（圆形），算法的鲁棒性较好，因此我们还需要利用实际的等震线进行鲁棒性测试。对已有等震线数据进行抽样模拟不同的等震线边界取样情景（抽样的示意图见图11），抽样率在0.80—0.99之间随机选取，模拟500次，将取样结果作为数据进行计算，比较每次计算的椭圆长短轴半径a和b，以检验算法的鲁棒性。由前文分析可知，简单闭合等震线和简单未闭合等震线形状简单，不同算法所得结果的一致性较好。因此我们选择形状更为复杂的闭合等震线和复杂未闭合等震线对鲁棒性进行检验。考虑算法的适用场景，使用未闭合复杂等震线检验无约束拟合算法（图12a，图13a）和约束长轴走向的拟合算法（图12b，图13b），使用复杂闭合等震线检验约束面积与长轴走向的拟合算法（图12c，图13c）和约束面积的拟合算法（图12d，图13d）。
![]()
![]()
12 算法鲁棒性检验（参数结果）灰色三角形为每次抽样计算的结果，红色虚线为抽样计算参数的一倍标准差，红色实线为抽样计算的平均结果，黑色虚线为全部数据的计算结果，直方图给出参数计算结果的统计情况（a） 无约束拟合算法；（b） 约束长轴走向的拟合算法12. Robustness test of the four algorithms （parameter results）The grey triangle is the result of each sampling calculation，the dashed red line is one standard deviation of the sampling calculation，the solid red line is the average result of sampling calculation，the black dashed line is the calculation result of the overall data，and the histogram gives the statistics of the parameter calculations（a） Unconstrained fitting algorithm；（b） Fitting algorithm with major-axis strike constrained![]()
12 算法鲁棒性检验（参数结果）灰色三角形为每次抽样计算的结果，红色虚线为抽样计算参数的一倍标准差，红色实线为抽样计算的平均结果，黑色虚线为全部数据的计算结果，直方图给出参数计算结果的统计情况（c） 约束面积与长轴走向的拟合算法；（d） 约束面积的拟合算法12. Robustness test of the four algorithms （parameter results）The grey triangle is the result of each sampling calculation，the dashed red line is one standard deviation of the sampling calculation，the solid red line is the average result of sampling calculation，the black dashed line is the calculation result of the overall data，and the histogram gives the statistics of the parameter calculations（c） Fitting algorithm with area and major-axis strike constrained；（d） Fitting algorithm with area constrained![]()
图 13 算法鲁棒性检验的结果示意图（a） 无约束拟合算法；（b） 约束长轴走向的拟合算法；（c） 约束面积与长轴走向的拟合算法；（d） 约束面积的拟合算法Figure 13. Robustness test of the four algorithms （diagram of the results）（a） Unconstrained fitting algorithm；（b） Fitting algorithm with major-axis strike constrained；（c） Fitting algorithm with area and major-axis strike constrained；（d） Fitting algorithm with area constrained图12中各子图依次是每次抽样试验的抽样率、长轴半径a、短轴半径b以及长轴走向θ的计算结果，图13展示的拟合椭圆结果。
由图12可见：长短轴半径计算结果大致稳定，参数计算结果大都落在平均值正负一倍标准差范围内，算法鲁棒性较好（图12a）；约束长轴走向的拟合椭圆算法能够较好地按照设定长轴走向得到拟合椭圆，计算得到的椭圆长短轴半径结果较为稳定（图12b）；闭合等震线长轴走向与目标走向相差较大时，附加面积约束可能会导致拟合结果偏向于圆（图12c）；面积约束算法的鲁棒性较好，计算结果较为稳定（图12d）。
鲁棒性检验结果说明：本文的四种算法鲁棒性较好，不同抽样率随机抽样计算结果大致相同，结果离散性较小，较为稳定；复杂未闭合等震线计算结果相对离散（图12b），在直接测量方法中这种问题依然存在，这是等震线自身情况造成的，对于是否使用这类等震线要单独进行判断，使用时要保证有足够的取样点；由于算法的约束性，当约束条件与实际等震线相差较大时，结果可能偏向于圆。
3.3 与已有方法测量结果的对比
将本文椭圆拟合计算的等震线长短轴数据与前人测量得到的结果进行对比。以1976年7月28日唐山MS7.8地震为例（图2a，b），选择长轴转向模式的计算结果（图2a）与陈达生和刘汉生（1989）直接测量法得到的结果以及沙海军等（2010）拟合椭圆法测量的结果进行对比，结果列于表1。以2008年5月12日汶川MS8.0地震为例（图2c，d），选择长轴转向模式的计算结果（图2c）与肖亮和俞言祥（2011）直接测量的结果对比，结果列于表2。
表 1 1976年7月28日唐山MS7.8地震等震线测量结果对比Table 1. Comparison of isoseismal measurement results of MS7.8 Tangshan earthquake on July 28，1976表 2 2008年5月12日汶川MS8.0地震等震线测量结果对比Table 2. Comparison of isoseismal measured by MS8.0 Wenchuan earthquake on May 12，2008等震线
烈度本文 肖亮和俞言祥（2 011） 长轴/km 短轴/km 长轴/km 短轴/km XI 33.61 10.15 33 10 XI 44.51 7.34 44 7.5 X 121.36 14.15 112.1 14 IX 161.47 25.24 159 22.5 VIII 212.09 60.41 206.5 57.5 VII 280.26 140.76 283 133.5 VI 471.23 295.57 468 298 由表1和表2所示结果可见：本文方法得到的椭圆长短轴参数与直接测量方法的结果大致相同，本文算法与直接测量方法中的直观法和等面积法（陈达生，刘汉生，1989）保持了继承性；对于类似椭圆分布的等震线（图2c，表2），拟合椭圆算法与直接测量方法得到的结果相差相对较小，只有唐山Ⅵ度等震线长短轴半径与直接测量方法的结果差别较大，在3.2节（图13a，b）对该等震线已进行了详细的计算分析。虽然唐山Ⅵ度等震线的整体形状较为复杂，但仍有趋向椭圆的趋势，因此选择椭圆进行建模是可以接受的。由于该等震线缺失部分较多且不完整，本文拟合椭圆方法效果较差。此外，将等震线拟合为椭圆也有其不足，如图2a和图2b中的烈度异常区，其形状与椭圆差异较大，虽然使用椭圆拟合略显不足，但本文的重点在于使用椭圆拟合的结果计算等震线的长轴和短轴数据，从而减少人为测量的主观性，而选择椭圆拟合等震线的优势在于更易于保持对直接测量方法的继承性。
4. 讨论与结论
本文考虑到与直接测量方法的一致性和继承性，以及国内常用的椭圆烈度衰减关系所需的数据形式，选择椭圆拟合等震线获取长短轴数据，将等震线测量的模式和直接测量方法中直观法和等面积法的概念以约束的形式加入到椭圆拟合算法，系统地给出了基于拟合椭圆测量等震线长短轴的方法。选择不同形状的等震线对椭圆拟合算法进行适用性测试和鲁棒性测试，并将等震线的测量结果与已有测量结果（陈达生，刘汉兴，1989；沙海军等，2010；肖亮，俞言祥，2011）进行对比，结果表明：① 本文的椭圆拟合算法能够根据等震线实际情况进行拟合，适用性较好，能够适应一般的噪声情况，鲁棒性较好，其中带有面积约束的算法稳定性更好；② 长轴不转向模式下需要约束模型的长轴走向，当约束条件与实际等震线之间差异较大时，拟合结果趋向于圆；③ 等震线长短轴测量结果与已有测量结果大致相同，相比直接测量方法，本文通过最小二乘法求解，减少了主观不确定性，而与郁曙君（1993）、沙海军等（2010）和吴清（2013）的椭圆拟合算法相比，本文提出的拟合算法更加全面，最大程度与常用的直接测量方法保持了衔接。
本文的研究能够完成大多数等震线长短轴数据的测量，而对于未闭合且形状较为复杂的等震线，当取样较少时，计算结果稳定性相对较差，需单独判断测量结果是否合理，直接测量方法中也存在这种问题（肖亮，俞言祥，2011）。利用更合理的模型和稳定的算法获取等震线长短轴数据是下一步工作的重点。此外，本研究结果可以为新一代烈度衰减关系的更新提供数据支持。
-
![]()
图 1 椭圆长轴走向示意图
（a） 走向范围为0°—90°；（b） 走向范围为−90°—0°
Figure 1. Diagram of the major-axis strike of the ellipse
（a） Strike ranges 0°−90°；（b） Strike ranges −90°−0°
![]()
图 2 长轴转向（左）和不转向（右）模式的地震等震线图拟合结果
（a，b） 1976年7月28日唐山MS7.8地震；（c，d） 2008年5月12日汶川MS8.0地震
Figure 2. Fitting results of the isoseismals of the earthquakes with major-axis steering （a） and non-steering （b） modes
（a，b） Tangshan MS7.8 earthquake on July 28，1976； （c，d） Wenchuan MS8.0 earthquake on May 12，2008
![]()
图 3 代表性的等震线
（a） 简单闭合等震线：2013年芦山MS7.0地震Ⅶ度等震线；（b） 复杂闭合等震线：1945年滦县MS6.3地震Ⅶ度等震线；（c） 简单未闭合等震线：1986年阳江MS4.8地震Ⅳ度等震线；（d） 复杂未闭合等震线：1976年唐山MS7.8地震Ⅵ度等震线
Figure 3. Representative isoseismals
（a） Closed isoseismal with simple shape：Isoseismal Ⅶ of 2013 Lushan MS7.0 earthquake；（b） Closed isoseismal with complex shape：Isoseismal Ⅶ of 1945 Luanxian MS6.3 earthquake；（c） Unclosed isoseismal with simple shape：Isoseismal Ⅳ of 1986 Yangjiang MS4.8 earthquake；（d） Unclosed isoseismal with complex shape：Isoseismal Ⅵ of 1976 Tangshan MS7.8 earthquake
![]()
图 4 四类等震线在不同无约束算法下的拟合结果的对比
（a） 简单闭合等震线：2013年芦山MS7.0地震Ⅶ度等震线；（b） 复杂闭合等震线：1945年滦县MS6.3地震Ⅶ度等震线；（c） 简单未闭合等震线：1986年阳江MS4.8地震Ⅳ度等震线；（d） 复杂未闭合等震线：1976年唐山MS7.8地震Ⅵ度等震线
Figure 4. Comparison of the fitting results of four isoseismals by the different unconstrained algorithms
（a） Closed isoseismal with simple shape：Isoseismal Ⅶ of 2013 Lushan MS7.0 earthquake；（b） Closed isoseismal with complex shape：Isoseismal Ⅶ of 1945 Luanxian MS6.3 earthquake；（c） Unclosed isoseismal with simple shape：Isoseismal Ⅳ of 1986 Yangjiang MS4.8 earthquake；（d） Unclosed isoseismal with complex shape：Isoseismal Ⅵ of 1976 Tangshan MS7.8 earthquake
![]()
图 5 四类等震线在无约束拟合算法和约束长轴走向拟合算法下拟合结果的对比
（a） 简单闭合等震线：2013年芦山MS7.0地震Ⅶ度等震线；（b） 复杂闭合等震线：1945年滦县MS6.3地震Ⅶ度等震线；（c） 简单未闭合等震线：1986年阳江MS4.8地震Ⅳ度等震线；（d） 复杂未闭合等震线：1976年唐山MS7.8地震Ⅵ度等震线
Figure 5. Comparison of the fitting results of the isoseismals by the unconstrained algorithm with those by the algorithm with major-axis strike constrained
（a） Closed isoseismal with simple shape：Isoseismal Ⅶ of 2013 Lushan MS7.0 earthquake；（b） Closed isoseismal with complex shape：Isoseismal Ⅶ of 1945 Luanxian MS6.3 earthquake；（c） Unclosed isoseismal with simple shape：Isoseismal Ⅳ of 1986 Yangjiang MS4.8 earthquake；（d） Unclosed isoseismal with complex shape：Isoseismal Ⅵ of 1976 Tangshan MS7.8 earthquake
![]()
图 6 遍历算法1和遍历算法2计算过程及结果对比
（a） 2013年4月20日芦山MS7.0地震Ⅶ度等震线；（b） 1945年9月23日滦县MS6.3地震Ⅶ度等震线
Figure 6. Comparison of calculation process and result of two traversal algorithms
（a） Isoseismal Ⅶ of Lushan MS7.0 earthquake on April 20，2013；（b） Isoseismal Ⅶ of Luanxian MS6.3 earthquake on September 23，1945
![]()
图 7 约束长轴走向的拟合算法和约束面积与长轴走向的拟合算法所得拟合结果的对比
（a） 2013年4月20日芦山MS7.0地震Ⅶ度等震线；（b） 1945年9月23日滦县MS6.3地震Ⅶ度等震线
Figure 7. Comparison of fitting results obtained by the algorithm with major-axis strike constrained with those by the algorithm with area and major-axis strike constrained
（a） Isoseismal Ⅶ of Lushan MS7.0 earthquake on April 20，2013；（b） Isoseismal Ⅶ of Luanxian MS6.3 earthquake on September 23，1945
![]()
图 8 无约束拟合算法和约束面积的拟合算法所得拟合结果的对比
（a） 2013年4月20日芦山MS7.0地震Ⅶ度等震线；（b） 1945年9月23日滦县MS6.3地震Ⅶ度等震线
Figure 8. Comparison of fitting results obtained by the unconstrained algorithm with those by the algorithm with area constrained
（a） Isoseismal Ⅶ of Lushan MS7.0 earthquake on April 20，2013；（b） Isoseismal Ⅶ of Luanxian MS6.3 earthquake on September 23，1945
![]()
图 9 误差水平σ＝20时四种算法的适用性检验
（a） 无约束拟合算法；（b） 约束长轴走向的拟合算法；（c） 约束面积与长轴走向的拟合算法；（d） 约束面积的拟合算法
Figure 9. Applicability test of four algorithms when error σ is taken as 20
（a） Unconstrained fitting algorithm；（b） Fitting algorithm with major-axis strike constrained；（c） Fitting algorithm with area and major-axis strike constrained；（d） Fitting algorithm with area constrained
![]()
图 10 误差水平σ＝100时四种算法的适用性检验
（a） 无约束拟合算法；（b） 约束长轴走向的拟合算法；（c） 约束面积与长轴走向的拟合算法；（d） 约束面积的拟合算法
Figure 10. Applicability test of four algorithms when error σ is taken as 100
（a） Unconstrained fitting algorithm；（b） Fitting algorithm with major-axis strike constrained；（c） Fitting algorithm with area and major-axis strike constrained；（d） Fitting algorithm with area constrained
![]()
![]()
12 算法鲁棒性检验（参数结果）
灰色三角形为每次抽样计算的结果，红色虚线为抽样计算参数的一倍标准差，红色实线为抽样计算的平均结果，黑色虚线为全部数据的计算结果，直方图给出参数计算结果的统计情况（a） 无约束拟合算法；（b） 约束长轴走向的拟合算法
12. Robustness test of the four algorithms （parameter results）
The grey triangle is the result of each sampling calculation，the dashed red line is one standard deviation of the sampling calculation，the solid red line is the average result of sampling calculation，the black dashed line is the calculation result of the overall data，and the histogram gives the statistics of the parameter calculations（a） Unconstrained fitting algorithm；（b） Fitting algorithm with major-axis strike constrained
![]()
12 算法鲁棒性检验（参数结果）
灰色三角形为每次抽样计算的结果，红色虚线为抽样计算参数的一倍标准差，红色实线为抽样计算的平均结果，黑色虚线为全部数据的计算结果，直方图给出参数计算结果的统计情况（c） 约束面积与长轴走向的拟合算法；（d） 约束面积的拟合算法
12. Robustness test of the four algorithms （parameter results）
The grey triangle is the result of each sampling calculation，the dashed red line is one standard deviation of the sampling calculation，the solid red line is the average result of sampling calculation，the black dashed line is the calculation result of the overall data，and the histogram gives the statistics of the parameter calculations（c） Fitting algorithm with area and major-axis strike constrained；（d） Fitting algorithm with area constrained
![]()
图 13 算法鲁棒性检验的结果示意图
（a） 无约束拟合算法；（b） 约束长轴走向的拟合算法；（c） 约束面积与长轴走向的拟合算法；（d） 约束面积的拟合算法
Figure 13. Robustness test of the four algorithms （diagram of the results）
（a） Unconstrained fitting algorithm；（b） Fitting algorithm with major-axis strike constrained；（c） Fitting algorithm with area and major-axis strike constrained；（d） Fitting algorithm with area constrained
表 1 1976年7月28日唐山MS7.8地震等震线测量结果对比
Table 1 Comparison of isoseismal measurement results of MS7.8 Tangshan earthquake on July 28，1976
下载: 导出CSV
表 2 2008年5月12日汶川MS8.0地震等震线测量结果对比
Table 2 Comparison of isoseismal measured by MS8.0 Wenchuan earthquake on May 12，2008
等震线
烈度本文 肖亮和俞言祥（2 011） 长轴/km 短轴/km 长轴/km 短轴/km XI 33.61 10.15 33 10 XI 44.51 7.34 44 7.5 X 121.36 14.15 112.1 14 IX 161.47 25.24 159 22.5 VIII 212.09 60.41 206.5 57.5 VII 280.26 140.76 283 133.5 VI 471.23 295.57 468 298
下载: 导出CSV
-
陈达生,刘汉兴. 1989. 地震烈度椭圆衰减关系[J]. 华北地震科学,7(3):31–42. Chen D S,Liu H X. 1989. Elliptical attenuation relationship of earthquake intensity[J]. North China Earthquake Science,7(3):31–42 (in Chinese).
国家地震局. 1996. 中国地震烈度区划图(1990)概论[M]. 北京:地震出版社:137−144. State Seismological Bureau. 1996. Introduction to China Seismic Intensity Zoning Maps (1990)[M]. Beijing:Seismological Press:137−144 (in Chinese).
胡聿贤. 2007. 地震安全性评价技术教程[M]. 北京:地震出版社:282−320. Hu Y X. 2007. A Technical Course on Seismic Safety Evaluation[M]. Beijing:Seismological Press:282−320 (in Chinese).
雷建成,高孟潭,俞言祥. 2007. 四川及邻区地震动衰减关系[J]. 地震学报,29(5):500–511. doi: 10.3321/j.issn:0253-3782.2007.05.007 Lei J C,Gao M T,Yu Y X. 2007. Seismic motion attenuation relations in Sichuan and adjacent areas[J]. Acta Seismologica Sinica,29(5):500–511 (in Chinese).
李建文,张成现,李颀,马学宗. 2010. 极小值搜索与平方和转化求解多元非线性方程组[J]. 陕西科技大学学报(自然科学版),28(3):143–147. doi: 10.3969/j.issn.1000-5811.2010.03.034 Li J W,Zhang C X,Li Q,Ma X Z. 2010. An algorithm of multivariant nonlinear equations by searching smallest number and quadratic sum converting[J]. Journal of Shaanxi University of Science &Technology,28(3):143–147 (in Chinese).
李永强,龚强,王景来. 2006. 基于GIS的数字等震线模型[J]. 地震研究,29(4):401–406. doi: 10.3969/j.issn.1000-0666.2006.04.015 Li Y Q,Gong Q,Wang J L. 2006. Digital isoseisms model based on GIS[J]. Journal of Seismological Research,29(4):401–406 (in Chinese).
沙海军,吕悦军,彭艳菊,唐荣余. 2010. 测量等震线长短轴的拟合椭圆法及其Matlab语言实现[G]//地壳构造与地壳应力文集(22). 北京:中国地震局地壳应力研究所:19–23. Sha H J,Lü Y J,Peng Y J,Tang R Y. 2010. A method of ellipse fitting and its apllication in measurement of long and short axes of isoseismal[G]//Bulletin of the Institute of Crustal Dynamics (22). Beijing:Institute of Crustal Dynamics,China Earthquake Administration:19–23 (in Chinese).
孙继浩. 2011. 川滇及邻区中强地震烈度衰减关系的适用性研究[D]. 北京:中国地震局地震预测研究所:14−41. Sun J H. 2011. Study of Moderate-Strong Seismic Intensity Attenuation Relations in Sichuan-Yunnan and Its Adjacent Areas[D]. Beijing:Institute of Earthquake Science,China Earthquake Administration:14−41 (in Chinese).
汪素云,俞言祥,高阿甲,阎秀杰. 2000. 中国分区地震动衰减关系的确定[J]. 中国地震,16(2):99–106. doi: 10.3969/j.issn.1001-4683.2000.02.001 Wang S Y,Yu Y X,Gao E J,Yan X J. 2009. Development of attenuation relations for ground motion in China[J]. Earthquake Research in China,16(2):99–106 (in Chinese).
吴清. 2013. 基于空间离散烈度点椭圆模型的历史地震参数估计方法研究[D]. 北京:中国地震局地球物理研究所:62−77. Wu Q. 2013. Study on Parameters Estimation Method of Historical Earthquakes Based on Discrete Intensity Points and Elliptical Model[D]. Beijing:Institute of Geophysics,China Earthquake Administration:62−77 (in Chinese).
肖亮,俞言祥. 2011. 中国西部地区地震烈度衰减关系[J]. 震灾防御技术,6(4):358–371. doi: 10.3969/j.issn.1673-5722.2011.04.002 Xiao L,Yu Y X. 2011. Earthquake intensity attenuation relationship in western China[J]. Technology for Earthquake Disaster Prevention,6(4):358–371 (in Chinese).
郁曙君. 1993. 确定烈度衰减关系的椭圆投影两步拟合法[J]. 地震学报,15(1):109–114. Yu S J. 1993. Two-step curve fitting with elliptic projection for determining intensity attenuation law[J]. Acta Seismologica Sinica, 15 (1):109−114 (in Chinese).
俞言祥,李山有,肖亮. 2013. 为新区划图编制所建立的地震动衰减关系[J]. 震灾防御技术,8(1):24–33. doi: 10.3969/j.issn.1673-5722.2013.01.003 Yu Y X,Li S Y,Xiao L. 2013. Development of ground motion attenuation relations for the new seismic hazard map of China[J]. Technology for Earthquake Disaster Prevention,8(1):24–33 (in Chinese).
Ahmadzadeh S,Doloei G J,Zafarani H. 2020. New intensity prediction equation for Iran[J]. J Seismol,24(1):23–35. doi: 10.1007/s10950-019-09882-7
Allen T I,Wald D J,Worden C B. 2012. Intensity attenuation for active crustal regions[J]. J Seismol,16(3):409–433. doi: 10.1007/s10950-012-9278-7
Ambraseys N N,Douglas J. 2004. Magnitude calibration of north Indian earthquakes[J]. Geophys J Int,159(1):165–206. doi: 10.1111/j.1365-246X.2004.02323.x
Anderson J G. 1978. On the attenuation of modified Mercalli intensity with distance in the United States[J]. Bull Seismol Soc Am,68(4):1147–1179.
Casado C L,Palacios S M,Delgado J,Peláez J A. 2000. Attenuation of intensity with epicentral distance in the Iberian Peninsula[J]. Bull Seismol Soc Am,90(1):34–47. doi: 10.1785/0119980116
Chandra U. 1979. Attenuation of intensities in the United States[J]. Bull Seismol Soc Am,69(6):2003–2024. doi: 10.1785/BSSA0690062003
Fitzgibbon A,Pilu M,Fisher R B. 1999. Direct least square fitting of ellipses[J]. IEEE Trans Pattern Anal Mach Intell, 21 (5):476–480. doi: 10.1109/34.765658
Ghosh G K,Mahajan A K. 2013. Intensity attenuation relation at Chamba–Garhwal area in north-west Himalaya with epicentral distance and magnitude[J]. J Earth Syst Sci,122(1):107–122. doi: 10.1007/s12040-012-0261-z
Howell B F Jr,Schultz T R. 1975. Attenuation of modified Mercalli intensity with distance from the epicenter[J]. Bull Seismol Soc Am,65(3):651–665. doi: 10.1785/BSSA0650030651
Musson R M W. 2005. Intensity attenuation in the U.K.[J]. J Seismol,9(1):73–86. doi: 10.1007/s10950-005-2979-4
Prieto J A,Foschi R O,Ventura C E,Finn W D L,Ramos A M,Prada F. 2011. Probability distribution of intensity attenuations for Colombia and western Venezuela[J]. Bull Seismol Soc Am,101(2):495–505. doi: 10.1785/0120090269
Rosin P L. 1993. A note on the least squares fitting of ellipses[J]. Pattern Recogn Lett,14(10):799–808. doi: 10.1016/0167-8655(93)90062-I
Sørensen M B,Stromeyer D,Grünthal G. 2010. Intensity attenuation in the Campania region,Southern Italy[J]. J Seismol,14(2):209–223. doi: 10.1007/s10950-009-9162-2
Xu W X,Wang C L,Yang W S,Yu D H,Zhang J G. 2021. Macroseismic intensity attenuation in western China[J]. J Seismol,25(2):711–731. doi: 10.1007/s10950-020-09973-w
计量
- 文章访问数: 163
- HTML全文浏览量: 40
- PDF下载量: 51
