利用烈度數據點估算中強歷史地震震級與震中的方法研究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.