吳 清 高孟潭

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

引言

歷史地震是分析區(qū)域地震活動(dòng)特征，預(yù)測(cè)地震活動(dòng)趨勢(shì)不可缺少的基礎(chǔ)資料。工程抗震設(shè)計(jì)工作中提出的地震危險(xiǎn)性分析方法，其歷史重演原則和地質(zhì)構(gòu)造類比原則也在很大程度上依賴于歷史地震。近年來，歷史地震研究取得了豐碩的成果，但始終存在著史料不均衡、不完備，地震參數(shù)誤差較大等問題。即使相同的史料，不同專家的理解也可能產(chǎn)生分歧，不同版本地震目錄所確定的同一地震參數(shù)可以出現(xiàn)較大差異。而新一代地震區(qū)劃圖的編制，對(duì)歷史地震參數(shù)精度、可信度都有更高的要求。因此，有必要探討估算歷史地震參數(shù)的新思路和新方法。

地震烈度是指某一地區(qū)地面及各人工建筑物遭受一次地震影響的強(qiáng)烈程度（胡聿賢，2006）。公元1900年以前，在沒有儀器可以記錄地震動(dòng)各種物理參數(shù)的時(shí)候，地震烈度資料就成為記錄歷史地震的唯一指標(biāo)。歷史文獻(xiàn)中關(guān)于地震破壞的記錄，以及震后烈度調(diào)查資料，是某個(gè)具體地點(diǎn)的破壞程度或烈度記錄，可以據(jù)其評(píng)定該地點(diǎn)的烈度值，稱為烈度數(shù)據(jù)點(diǎn)。烈度數(shù)據(jù)點(diǎn)從某種意義上說與地震臺(tái)站有一定的相似處，仿照臺(tái)站數(shù)據(jù)處理方法，用烈度數(shù)據(jù)點(diǎn)對(duì)一些主要地震參數(shù)進(jìn)行定量估計(jì)，理論上切實(shí)可行（謝毓壽，1991），更是現(xiàn)今研究歷史地震的重要手段。

基于Peruzza（1992）提出的殘差最小化方法，Bakun等（1997）以烈度數(shù)據(jù)點(diǎn)為基礎(chǔ)，利用衰減關(guān)系通過網(wǎng)格搜索來試算歷史地震震級(jí)和震中區(qū)域，這種方法曾應(yīng)用于美國加州和美國西北太平洋地區(qū)（Bakun等，1997；Bakun，1999；Bakun等，2003；Bakun，2006a；Diane，2009）、日本（Bakun，2005）、德國（Hinzen等，2001）、法國（Bakun等，2006b）及我國華北地區(qū)（張揚(yáng)等，2009）。但該方法是以圓烈度衰減關(guān)系為基礎(chǔ)，沒有考慮烈度衰減的方向性問題，而且宏觀震中與儀器震中有一定的混淆。Gasperini等（1999）用一套基于烈度數(shù)據(jù)點(diǎn)的Boxer方法來估算大震級(jí)歷史地震的震源參數(shù)。但該方法首先將最高烈度點(diǎn)的幾何中心定為震中，再分步得到其它參數(shù)，沒有考慮其它烈度點(diǎn)分布對(duì)震中位置的影響；而且震級(jí)也是用圓烈度衰減關(guān)系進(jìn)行估算。Musson等（2008）綜合了以上幾種方法的優(yōu)點(diǎn)，采用K?vesligethy（1906）模型網(wǎng)格搜索最小均方根來估計(jì)震中，然后用 Frankel（1994）的有感面積方法確定震級(jí)大小，這個(gè)過程仍然沒有考慮地震衰減的方向性。

一般來說，烈度隨震級(jí)的加大而加大，隨著到震中距離的加大而減小，而且震害分布總是呈現(xiàn)出明顯的橢圓形狀。我國研究者注意到這個(gè)現(xiàn)象并用橢圓模型來建立烈度衰減關(guān)系（沈建文等，1989；陳達(dá)生等，1989；霍俊榮，1989；胡聿賢，1999；汪素云等，2000）。但這些衰減關(guān)系是以等震線為基礎(chǔ)數(shù)據(jù)回歸得來，而等震線是原始烈度點(diǎn)的外包線，不是對(duì)烈度值的平均估計(jì)?；艨s（1989）曾提出，相同距離處另賦等震線低一度值的方法，近似等效烈度數(shù)據(jù)點(diǎn)來回歸橢圓衰減關(guān)系，但這種等效法只是在等震線基礎(chǔ)上進(jìn)行數(shù)據(jù)處理，不是對(duì)真實(shí)烈度點(diǎn)的應(yīng)用。

因此，建立一種直接基于烈度數(shù)據(jù)點(diǎn)的方法，充分考慮地震烈度的橢圓分布，盡可能利用全部烈度數(shù)據(jù)點(diǎn)來估算歷史地震宏觀震中和震級(jí)，對(duì)于提高歷史地震參數(shù)的精度具有重要意義。

1 方法原理

本文首先以烈度數(shù)據(jù)點(diǎn)為基礎(chǔ)，建立一種橢圓地震烈度分布模型來描述烈度隨震級(jí)和距離而變化的關(guān)系，即假若已知一次地震的震級(jí)，此模型應(yīng)能給出各地表點(diǎn)的烈度。然后基于所建烈度分布模型，聯(lián)立考慮中心點(diǎn)和方向性的地震參數(shù)估計(jì)方程，代入烈度數(shù)據(jù)點(diǎn)信息，估算地震震級(jí)和震中。主要的技術(shù)路線如圖1所示。

圖1 本文技術(shù)路線示意圖Fig. 1 The schematic diagram of the technology in this paper

一般來說，對(duì)于某一次地震，震級(jí)越大，距離震中越近，其地震烈度越大。我國地質(zhì)條件復(fù)雜，烈度分布隨地區(qū)變化差異較大。一般看來，地震烈度呈橢圓分布，即沿發(fā)震斷層方向上烈度分布范圍較廣，而與斷層垂直方向上烈度分布較窄，這種方向上的差別到了遠(yuǎn)場(chǎng)逐漸消失，最終趨于圓形。為了反映出這種烈度分布，本文借鑒陳達(dá)生等（1989）提出的地震烈度橢圓衰減關(guān)系長短軸統(tǒng)一回歸的思想，采用如下橢圓烈度分布模型：

寫成長、短軸方向的形式為：

沿長軸方向

沿短軸方向

式中，I為地震烈度；M為震級(jí)；系數(shù)1aC 、1bC 、2C、3aC 、3bC 均為回歸常數(shù)；oaR 、obR 分別為橢圓長、短軸兩個(gè)方向上烈度近場(chǎng)飽和因子；ε為回歸分析中表示不確定性的隨機(jī)變量，通常假定為對(duì)數(shù)正態(tài)分布，其均值為0，標(biāo)準(zhǔn)差為σ。

這里需要強(qiáng)調(diào)指出的是，這個(gè)模型以烈度數(shù)據(jù)點(diǎn)為基礎(chǔ)，先用最小二乘法對(duì)地震各烈度區(qū)原始烈度點(diǎn)空間分布進(jìn)行橢圓擬合，得到各個(gè)烈度區(qū)烈度點(diǎn)空間展布的橢圓擬合線，賦予這條橢圓擬合線該烈度區(qū)的烈度值，本文稱之為烈度估計(jì)線。此時(shí)的aR、bR分別是烈度I的橢圓估計(jì)線長、短半軸長度。

該模型在震中處的烈度值隨震級(jí)的變化率保持為常數(shù)C2，所以在任意震級(jí)長、短軸方向上震中處的烈度值均相等，從而保證回歸后不需要通過參數(shù)調(diào)整來解決近場(chǎng)烈度差異的問題（肖亮，2011），而中間距離仍保持長、短軸烈度的差別，同時(shí)在遠(yuǎn)場(chǎng)也使烈度分布成圓形。

考慮中心點(diǎn)和方向性的橢圓數(shù)學(xué)方程見式（4），本文用此方程來估計(jì)地震震級(jí)與宏觀震中：

式中，(xi, yi)為橢圓上的點(diǎn)坐標(biāo)； ( x0, y0)為橢圓中心坐標(biāo)； θ (0° ≤ θ ≤ 1 80°)是橢圓長軸逆時(shí)針與x軸正方向的夾角； Ra、 Rb分別為橢圓長半軸和短半軸長度。

由式（2）得：

由式（3）得：

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

在同一個(gè)烈度點(diǎn)應(yīng)有：

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

在式（9）中， C1a、 C1b、 C2、 C3a、 C3b、 Roa、 Rob由橢圓烈度分布模型統(tǒng)計(jì)回歸可得。若已知一個(gè)地震的多個(gè)烈度信息點(diǎn)，帶入到式（9），即可聯(lián)立方程組求得震中震級(jí)M和方向θ。有4個(gè)未知數(shù)，則至少需要4個(gè)方程確定一組解，但 x0與 y0不是相互獨(dú)立的，因此在已知橢圓烈度分布模型下，3個(gè)烈度信息點(diǎn)就能確定一組震中、震級(jí)和方向。烈度信息點(diǎn)增多，變成求解超定非線性方程組，采用通用全局優(yōu)化法的數(shù)值計(jì)算方法尋求最優(yōu)解，可非常直觀地求得地震宏觀震中位置、震級(jí)和方向。

2 地震數(shù)據(jù)及橢圓烈度分布模型

我國目前公開發(fā)表的烈度調(diào)查資料大多采用等震線的形式發(fā)布，原始烈度調(diào)查點(diǎn)通常難以獲得。本文選取最新出版的4冊(cè)《中國震例（1992—2002）》（陳棋福，2002a；2002b；2003；2008），提取其中既有儀器測(cè)量數(shù)據(jù)又有宏觀烈度調(diào)查數(shù)據(jù)的現(xiàn)代地震，將烈度分布圖數(shù)字化。這些地震的等震線在圖上部分明確標(biāo)示了宏觀烈度調(diào)查點(diǎn)位置；對(duì)于未明示烈度調(diào)查點(diǎn)位置的，本文將等震線在圖上能明確獲取烈度信息的城、縣、鎮(zhèn)、村所在地算作烈度點(diǎn)。經(jīng)數(shù)字化配準(zhǔn)后，獲取各地震烈度點(diǎn)的空間分布。受到提取烈度點(diǎn)有效數(shù)字信息的目的所限，本文獲取的地震基本分布在青、甘、川、滇4省，地震震級(jí)集中在MS4.9—6.6級(jí)，如表1所示。

將式（4）改寫為考慮中心點(diǎn)和方向性的橢圓參數(shù)方程：

對(duì)于既有儀器測(cè)量數(shù)據(jù)又有宏觀調(diào)查烈度數(shù)據(jù)的現(xiàn)代地震，可獲得地震的宏觀震中位置(x0, y0)、極震區(qū)走向θ和烈度分布點(diǎn)坐標(biāo)根據(jù)各烈度區(qū)調(diào)查點(diǎn)的空間分布基于式（10）利用最小二乘原理采用通用全局優(yōu)化算法擬合出各烈度區(qū)橢圓估計(jì)線長半軸和短半軸長度 Ra和 Rb，由此獲得各烈度區(qū)的橢圓烈度估計(jì)線。由于地震的宏觀震中并不一定都是該地震各烈度區(qū)的幾何中心，因此對(duì)某些地震需要通過式（10）同時(shí)擬合出各烈度區(qū)的幾何中心。這里的烈度估計(jì)線長軸方向均取極震區(qū)走向θ。

以各自宏觀震中為原點(diǎn)（0，0），將各地震對(duì)應(yīng)烈度點(diǎn)都轉(zhuǎn)換到了平面坐標(biāo)下。表1列出了擬合所得的各地震、各烈度區(qū)橢圓烈度估計(jì)線的長半軸和短半軸長度，擬合烈度估計(jì)線幾何中心不同于宏觀震中的也同樣列在了表1中。表1中擬合烈度估計(jì)線的幾何中心與宏觀震中重合的，坐標(biāo)即為（0，0）；幾何中心與宏觀震中不重合的，再具體列出。

表1 地震數(shù)據(jù)及各烈度區(qū)橢圓烈度估計(jì)線長、短半軸長度和中心點(diǎn)Table 1 List of parameters of the earthquakes in the study area

續(xù)表

續(xù)表

將表1中各地震、各烈度區(qū)的橢圓烈度估計(jì)線的長、短半軸長度及對(duì)應(yīng)烈度值和震級(jí)值，代入到式（2）、（3）進(jìn)行最小二乘統(tǒng)計(jì)回歸，得到地震烈度分布模型：

不同地區(qū)的震源特性、傳播介質(zhì)與場(chǎng)地條件都可能不同，地震烈度分布規(guī)律自然也可能不同。因此，烈度分布模型以按大地區(qū)分別使用為宜。受獲取地震原始烈度調(diào)查點(diǎn)所限，本文所得的烈度分布模型僅適用于青、甘、川、滇4省的中強(qiáng)地震。

3 震例驗(yàn)算

由上節(jié)得到橢圓烈度分布模型的回歸系數(shù)后，式（9）即變?yōu)椋?/p>

為了驗(yàn)證地震參數(shù)估計(jì)方程的有效性，本文搜集了《2001—2005年中國大陸地震災(zāi)害損失評(píng)估匯編》（中國地震局震災(zāi)應(yīng)急救援司，2010）里的地震烈度資料，提取了5個(gè)發(fā)生于模型適用區(qū)域的地震的烈度數(shù)據(jù)點(diǎn)，將所有烈度數(shù)據(jù)點(diǎn)以各自地震的宏觀震中為原點(diǎn)（0，0）轉(zhuǎn)換到平面坐標(biāo)下（xi，yi，Ii），然后直接代入到式（13）中進(jìn)行驗(yàn)算，計(jì)算結(jié)果見表2。這5個(gè)地震空間位置在研究區(qū)內(nèi)分布合理，震級(jí)也具有一定的代表性。

表2 現(xiàn)代震例計(jì)算結(jié)果Table 2 Estimated results of the modern earthquakes

由表2可以看到，參與驗(yàn)算的5個(gè)地震，計(jì)算震級(jí)與儀器震級(jí)的差值不超過0.6級(jí)，計(jì)算震中到宏觀震中的距離不超過6km，由此可見此方法是可行的。

4 不確定性分析

在震中位置和震級(jí)的估算過程中，烈度點(diǎn)的數(shù)量、分布和精度好壞將直接影響估算結(jié)果。由于歷史地震烈度點(diǎn)相對(duì)較少，為了對(duì)歷史地震參數(shù)進(jìn)行較好的估計(jì)，本文運(yùn)用蒙特卡洛方法，對(duì)信息豐富的現(xiàn)代地震烈度點(diǎn)抽樣，模擬歷史地震烈度點(diǎn)分布情況。由于至少需要3個(gè)烈度點(diǎn)才能確定一組解，先從烈度點(diǎn)中隨機(jī)抽取3個(gè)點(diǎn)，并且有放回的重復(fù)抽樣1000次。增加控制信息，逐點(diǎn)添加烈度點(diǎn)數(shù)，隨機(jī)重復(fù)抽取4，5，6……個(gè)點(diǎn)各1000次，這樣模擬出大量的烈度點(diǎn)組樣本。每個(gè)烈度點(diǎn)有可能被多次選擇或者不被選擇，每次抽樣代表不同的烈度點(diǎn)空間分布，而不同的烈度點(diǎn)數(shù)代表了烈度信息的多寡，真實(shí)歷史地震資料記載只相當(dāng)于隨機(jī)抽樣中的一種情況。通過蒙特卡洛方法模擬出大量烈度點(diǎn)組樣本代入式（13）中計(jì)算，由此來分析此方法計(jì)算地震震級(jí)與震中的不確定性。可以推測(cè)，烈度點(diǎn)數(shù)越多，分布越均勻，得出的震中與震級(jí)越精確，不確定性越小。

式（13）是超越方程，當(dāng)參與計(jì)算的烈度點(diǎn)數(shù)多于3個(gè)點(diǎn)時(shí)，更成為超定的超越方程組。這類矛盾方程組在嚴(yán)格的數(shù)學(xué)意義下無解，但是可以在最小二乘意義下，最大限度擬合數(shù)據(jù)，尋求最優(yōu)解。隨機(jī)抽樣得到的烈度點(diǎn)組會(huì)存在不同的空間分布，利用數(shù)學(xué)軟件采用全局最優(yōu)化的數(shù)值計(jì)算方法，可以得到每組抽樣點(diǎn)滿足式（13）的最優(yōu)數(shù)學(xué)解，但得到的數(shù)學(xué)解有可能并不滿足所求參數(shù)的物理意義。因此，所得的結(jié)果需要剔除與地震參數(shù)物理意義不符的數(shù)學(xué)解，然后對(duì)最終結(jié)果進(jìn)行統(tǒng)計(jì)分析。

本文以2001年5月24日云南寧蒗—四川鹽源5.8級(jí)地震為例，給出抽樣計(jì)算結(jié)果。寧蒗—鹽源地震有50個(gè)烈度點(diǎn)，以宏觀震中為原點(diǎn)轉(zhuǎn)換到平面坐標(biāo)下，然后從這50個(gè)烈度點(diǎn)里隨機(jī)抽取3點(diǎn)、4點(diǎn)、5點(diǎn)……各1000次，分別帶入到式（13）進(jìn)行解算，計(jì)算震中分布見圖2。將計(jì)算震級(jí)與儀器震級(jí)做差值，差值分布見圖3。由圖2可以看到，隨著烈度點(diǎn)的增多，計(jì)算震中逐漸向宏觀震中靠攏。由圖3可以看到，隨著烈度點(diǎn)的增多，計(jì)算震級(jí)與儀器震級(jí)的差值也逐漸收斂。

圖2 2001年5月24日云南寧蒗—四川鹽源5.8級(jí)地震抽樣計(jì)算震中分布（R為計(jì)算震中到宏觀震中平面距離）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級(jí)地震計(jì)算震級(jí)與儀器震級(jí)差值分布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個(gè)地震均進(jìn)行抽樣計(jì)算，統(tǒng)計(jì)不同抽樣點(diǎn)數(shù)情況下，計(jì)算震中到宏觀震中距離的均值與標(biāo)準(zhǔn)差，以及計(jì)算震級(jí)與儀器震級(jí)差值的均值與標(biāo)準(zhǔn)差，見圖4（a）、（b）。

由圖4（a）可以看到，隨著烈度點(diǎn)的增多，計(jì)算震中到宏觀震中距離的均值和標(biāo)準(zhǔn)差都逐漸減小。3點(diǎn)時(shí)，距離均值最大有18km，標(biāo)準(zhǔn)差最大接近12km；6點(diǎn)時(shí)，標(biāo)準(zhǔn)差到10km以下；8點(diǎn)時(shí)，均值也到10km以下。歷史地震震中精度分布為1類≤10km，2類≤25km，3類≤50km，4類≤100km，5類>100km。所以按照本文所提算法，一般達(dá)到6點(diǎn)時(shí)，標(biāo)準(zhǔn)差就能達(dá)到10km以下，即達(dá)到歷史地震的1類精度。

由圖4（b）可以看到，隨著烈度點(diǎn)的增加，震級(jí)差值的均值有個(gè)增大過程，達(dá)到一定點(diǎn)數(shù)之后趨于平穩(wěn)。這是因?yàn)榱叶赛c(diǎn)數(shù)較少時(shí)，代表的地震影響范圍相對(duì)較小，由此估算所得的震級(jí)自然偏??；而達(dá)到一定點(diǎn)數(shù)之后，烈度點(diǎn)分布所能反映的影響范圍趨于穩(wěn)定，計(jì)算出來的震級(jí)也就趨于穩(wěn)定。模型本身是對(duì)一個(gè)地區(qū)烈度分布的綜合估計(jì)，與任何一個(gè)真實(shí)地震之間都存在著系統(tǒng)偏差。不同地震之間的震源特性與場(chǎng)地條件不同，系統(tǒng)偏差也各不相同。而隨著烈度點(diǎn)的增加，震級(jí)差值的標(biāo)準(zhǔn)差整體上呈減小趨勢(shì)。對(duì)于參與試算的 5個(gè)地震，3點(diǎn)時(shí)，震級(jí)差值的標(biāo)準(zhǔn)差最大接近0.45級(jí)；8點(diǎn)時(shí)，標(biāo)準(zhǔn)差都集中到0.3級(jí)以下。

通過對(duì)《中國歷史強(qiáng)震目錄（公元前23世紀(jì)—公元1911年）》（國家地震局震害防御司，1995）中歷史強(qiáng)震烈度點(diǎn)數(shù)的統(tǒng)計(jì)，發(fā)現(xiàn)歷史地震烈度記載點(diǎn)普遍較少，18個(gè)及以上烈度點(diǎn)的地震，僅占?xì)v史強(qiáng)震總數(shù)的8%。因此，本文暫且分析20個(gè)及以下烈度點(diǎn)，計(jì)算地震參數(shù)不確定性的情況。

圖4 （a） 不同烈度點(diǎn)數(shù)情況下計(jì)算震中到宏觀震中距離的均值與標(biāo)準(zhǔn)差分布圖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） 不同烈度點(diǎn)數(shù)情況下計(jì)算震級(jí)與儀器震級(jí)之差值的均值與標(biāo)準(zhǔn)差分布圖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

不同地震、不同烈度點(diǎn)數(shù)情況下，不確定性各有不同，但整體分布趨勢(shì)基本相同。因此，可以按照烈度點(diǎn)數(shù)情況，對(duì)估定歷史地震參數(shù)的不確定性做粗略估計(jì)。

5 歷史地震參數(shù)的估算

通過前面的分析，本文建立了一套基于烈度數(shù)據(jù)點(diǎn)，用橢圓烈度分布模型估算地震宏觀震中與震級(jí)的方法，旨在用此方法來估算只有烈度記載且烈度數(shù)據(jù)點(diǎn)相對(duì)較少的歷史地震參數(shù)。以云南地區(qū)為例，將此方法用于歷史地震進(jìn)行試算。從《中國歷史地震圖集》（明、清時(shí)期）（下文簡稱《圖集》）（國家地震局地球物理研究所等，1986；1990）中提取4個(gè)烈度點(diǎn)較少的中強(qiáng)地震，獲取有歷史記載的烈度數(shù)據(jù)點(diǎn)，然后用本文所提算法估算歷史地震參數(shù)。

5.1 1599年10月16日云南嵩明5.0級(jí)地震

萬歷二十七年八月二十八日（公元1599年10月16日）云南嵩明地震，《圖集》給出地震參數(shù)為震中（103.0°E，25.3°N），震級(jí)5.0級(jí)，有5個(gè)烈度記載點(diǎn)，圖5是其烈度記載點(diǎn)分布圖。用本文的算法計(jì)算，計(jì)算震中為（102.99°E，25.28°N），與《圖集》震中相距2.42km；計(jì)算震級(jí)為5.2級(jí)，比《圖集》震級(jí)大0.2級(jí)。

圖5 云南嵩明5.0級(jí)地震烈度點(diǎn)分布圖Fig. 5 Distribution of intensity data points of Songming MS 5.0 earthquake in Yunnan

5.2 1696年7月7日云南昆明5.5級(jí)地震

康熙三十五年六月初九（公元1696年7月7日）云南昆明地震，《圖集》給出地震參數(shù)為震中（102.7°E，25.1°N），震級(jí)5.5級(jí)，有12個(gè)烈度記載點(diǎn)，圖6是其烈度記載點(diǎn)分布圖。用本文的算法計(jì)算，震中為（102.68°E，25.07°N），與《圖集》震中相距 3.8km；計(jì)算震級(jí)為5.2級(jí)，比《圖集》震級(jí)小0.3級(jí)。

圖6 云南昆明5.5級(jí)地震烈度點(diǎn)分布圖Fig. 6 Distribution of intensity data points of Kunming MS 5.5 earthquake in Yunnan

5.3 1876年8月5日云南永平6.0級(jí)地震

光緒二年六月十六日（公元1876年8月5日）云南永平地震，《圖集》給出地震參數(shù)為震中（99.5°E，25.5°N），震級(jí)6.0級(jí)，有8個(gè)烈度記載點(diǎn)。圖7是其烈度記載點(diǎn)分布圖。用本文的算法計(jì)算，計(jì)算震中為（99.45°E，25.55°N），與《圖集》震中相距 7.04km；計(jì)算震級(jí)為5.9級(jí)，比《圖集》震級(jí)小0.1級(jí)。

圖7 云南永平6.0級(jí)地震烈度點(diǎn)分布圖Fig. 7 Distribution of intensity data points of Yongping MS 6.0 earthquake in Yunnan

5.4 1909年5月11日云南米勒西南6.5級(jí)地震

宣統(tǒng)元年三月二十二日夜（公元1909年5月11日）云南米勒、寧州地震，《圖集》給出地震參數(shù)為震中（103.1°E，24.3°N），震級(jí)6.5級(jí)，有16個(gè)烈度記載點(diǎn)，圖8是其烈度記載點(diǎn)分布圖。用本文的算法計(jì)算，計(jì)算震中為（103.11°E，24.30°N），與《圖集》震中相距0.63km；計(jì)算震級(jí)為6.8級(jí)，比《圖集》震級(jí)大0.3級(jí)。

圖8 云南米勒西南6.5級(jí)地震烈度點(diǎn)分布圖Fig. 8 Distribution of intensity data points of Mile MS 6.5 earthquake in Yunnan

由以上震例可見，本文的方法對(duì)烈度數(shù)據(jù)點(diǎn)較少、等震線不易勾畫的歷史地震頗為有效。處理過程直接明了，減少了等震線勾畫過程中的主觀不確定性，可重復(fù)和再現(xiàn)。而且對(duì)于1599年10月16日云南嵩明5.0級(jí)地震和1876年8月5日云南永平6.0級(jí)地震，這類烈度點(diǎn)分布極不均勻、基本集中在極震區(qū)一側(cè)的地震，用本文的算法也能得到較為可信的地震參數(shù)。由此可見，本文的方法對(duì)烈度信息空間分布散亂的穩(wěn)定性。

6 結(jié)論與討論

本文首先基于地震原始烈度調(diào)查點(diǎn)，用橢圓數(shù)學(xué)方程對(duì)烈度點(diǎn)空間分布進(jìn)行最小二乘擬合，得到各烈度區(qū)橢圓烈度估計(jì)線。進(jìn)而對(duì)烈度估計(jì)線進(jìn)行統(tǒng)計(jì)回歸，得到一組適用于青、甘、川、滇4省中強(qiáng)地震的橢圓烈度分布模型該模型的建立充分利用了大量空間分布的原始烈度信息，通過數(shù)學(xué)手段得到烈度估計(jì)線，避免了等震線圈定過程中的主觀經(jīng)驗(yàn)因素，可還原和重現(xiàn)。相較于烈度外包線的地震等震線，基于烈度點(diǎn)的橢圓估計(jì)線更能代表烈度分布的真實(shí)情況。

根據(jù)所建立的橢圓烈度分布模型，聯(lián)立考慮中心點(diǎn)和方向性的橢圓數(shù)學(xué)方程，代入烈度數(shù)據(jù)點(diǎn)信息，直接求解歷史地震震中和震級(jí)，過程直接明了。通過蒙特卡洛方法，定量分析所得參數(shù)的不確定性，可以看到隨著烈度點(diǎn)的增多，計(jì)算震中到宏觀震中距離的均值和標(biāo)準(zhǔn)差都逐漸減小，計(jì)算震級(jí)與儀器震級(jí)差值的均值有個(gè)增大過程，達(dá)到一定點(diǎn)數(shù)之后趨于平穩(wěn)。

本文的算法符合橢圓模型的數(shù)學(xué)、物理概念，滿足烈度分布的地震學(xué)規(guī)律，充分利用了大量空間分布的原始烈度信息，過程簡單直接。算例表明，計(jì)算結(jié)果與實(shí)際地震記錄參數(shù)在一定范圍內(nèi)吻合。將此方法應(yīng)用到僅有烈度資料的歷史地震，可以減少人為判定震中位置和震中烈度，進(jìn)而估算震級(jí)過程中人為的主觀不確定性，提高歷史地震震中和震級(jí)的精度。

由于中國大陸的地震絕大部分發(fā)生在地殼以內(nèi)，其震源深度差別不大（劉百篪等，2002），因此本文暫且不考慮震源深度對(duì)烈度分布的影響。

地震震害一般呈橢圓分布，但不一定是數(shù)學(xué)上的橢圓。這里將地震烈度分布近似為橢圓系列時(shí)，已經(jīng)包含了數(shù)學(xué)抽象的不確定性。但即要建立模型，簡化及隨之而來的不確定性或誤差就不可避免（沈建文等，1989）。地震烈度分布模型本身是從多個(gè)地震中提取的平均規(guī)律，與地震資料之間存在著離散性，與任何一個(gè)真實(shí)地震之間都有系統(tǒng)偏差，在以后的工作中需要進(jìn)一步改進(jìn)模型，以期減小這種系統(tǒng)偏差。

本文所建模型至少需要3個(gè)烈度點(diǎn)才能求解，而實(shí)際歷史地震資料中，僅有1個(gè)或2個(gè)烈度點(diǎn)的地震占了資料的半數(shù)以上，這類地震震中參數(shù)該如何獲取，需要進(jìn)一步討論。另外，式（13）中同時(shí)解算出了橢圓長軸與x軸正方向的夾角θ，這個(gè)夾角θ應(yīng)該與發(fā)震構(gòu)造或者破裂參數(shù)有關(guān)，但相關(guān)程度與其不確定性還有待進(jìn)一步研究。

最后，本文所提算法的基礎(chǔ)是地震烈度數(shù)據(jù)點(diǎn)。但是現(xiàn)代地震的烈度評(píng)定與歷史地震的烈度評(píng)定之間也存在著差異，如何將現(xiàn)代地震烈度與歷史地震烈度擴(kuò)展銜接，也是值得深入思考的問題。

陳達(dá)生，劉漢興，1989. 地震烈度橢圓衰減關(guān)系. 華北地震科學(xué)，7（8）：31—42.

陳棋福，鄭大林，車時(shí)，2002a. 中國震例（1992—1994）. 北京：地震出版社.

陳棋福，鄭大林，劉桂萍等，2002b. 中國震例（1995—1996）. 北京：地震出版社.

陳棋福，鄭大林，高榮勝，2003. 中國震例（1997—1999）. 北京：地震出版社.

陳棋福，鄭大林，車時(shí)等，2008. 中國震例（2000—2002）. 北京：地震出版社.

國家地震局地球物理研究所，復(fù)旦大學(xué)中國歷史地理研究所，1986. 中國歷史地震圖集（明時(shí)期）. 北京：地圖出版社.

國家地震局地球物理研究所，復(fù)旦大學(xué)中國歷史地理研究所，1990. 中國歷史地震圖集（清時(shí)期）. 北京：中國地圖出版社.

國家地震局震害防御司，1995. 中國歷史強(qiáng)震目錄（公元前23世紀(jì)-公元1911年）. 北京：地震出版社.

胡聿賢，1999. 地震安全性評(píng)價(jià)技術(shù)教程. 北京：地震出版社.

胡聿賢，2006. 地震工程學(xué). 北京：地震出版社.

霍俊榮，1989. 近場(chǎng)強(qiáng)地面運(yùn)動(dòng)衰減規(guī)律的研究 [博士論文]. 黑龍江：國家地震局工程力學(xué)研究所.

劉百篪，鄭文俊，郭華等，2002. 活斷層工作方法在中早期歷史地震研究中的作用. 中國地震，18（3）：283—288.

沈建文，華宜平，1989. 關(guān)于地震烈度衰減模型的系統(tǒng)偏差. 地震學(xué)報(bào)，11（1）：38—45.

汪素云，俞言祥，高阿甲等，2000. 中國分區(qū)地震動(dòng)衰減關(guān)系的確定. 中國地震，16（2）：99—106.

謝毓壽，1991. 歷史地震參數(shù)的估定. 東北地震研究，7（5）：1—6.

肖亮，2011. 水平向基巖強(qiáng)地面運(yùn)動(dòng)參數(shù)衰減關(guān)系研究 [博士論文]. 北京：中國地震局地球物理研究所.

張揚(yáng)，馬干，史保平等，2009. 華北地區(qū)烈度衰減模型建立及其用于震中區(qū)域和震級(jí)的定量估算. 地震學(xué)報(bào)，31（3）：290—306．

中國地震局震災(zāi)應(yīng)急救援司，2010．2001—2005年中國大陸地震災(zāi)害損失評(píng)估匯編．北京：地震出版社，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.