基于烈度數(shù)據(jù)點(diǎn)的歷史強(qiáng)震參數(shù)估計(jì)研究

吳清 高孟潭

中國(guó)地震局地球物理研究所，北京市海淀區(qū)民族大學(xué)南路5號(hào) 100081

0 引言

地震監(jiān)測(cè)預(yù)報(bào)和工程地震是地震工作者為了減輕地震災(zāi)害面臨的兩大艱巨任務(wù)（謝毓壽，1991），地震活動(dòng)性研究是上述工作的基礎(chǔ)。由于地震儀器記錄的歷史較短，不足以在短時(shí)間、小尺度范圍內(nèi)弄清楚地震活動(dòng)規(guī)律，因而對(duì)歷史地震的研究非常重要。尤其是對(duì)于6.5級(jí)以上的歷史強(qiáng)震，其參數(shù)的不確定性，會(huì)影響人們對(duì)某個(gè)地點(diǎn)或地區(qū)地震活動(dòng)狀況的認(rèn)識(shí)，影響某項(xiàng)重大工程的地震危險(xiǎn)性評(píng)估結(jié)果，更可能成為地震科學(xué)研究或防震減災(zāi)工作中一些重大問(wèn)題的關(guān)鍵。

6.5級(jí)以上強(qiáng)震的烈度一般在破裂方向上衰減較慢，在垂直破裂方向上迅速衰減，地震影響場(chǎng)呈狹長(zhǎng)的帶狀分布，普通的圓烈度模型和橢圓模型都不再適用，更適用的是斷層破裂模型，但是斷層破裂模型的帶狀跑道型分布較為復(fù)雜，在數(shù)學(xué)上不易處理。根據(jù)6.5級(jí)以上強(qiáng)震烈度分布的圖像特征，本文提出了基于烈度數(shù)據(jù)點(diǎn)的考慮斷層破裂長(zhǎng)度的烈度橢圓分布模型。即在烈度橢圓分布模型的基礎(chǔ)上，加入斷層破裂長(zhǎng)度的約束，要求最內(nèi)圈烈度估計(jì)等值線長(zhǎng)軸長(zhǎng)度大于斷層破裂長(zhǎng)度，使得烈度橢圓分布模型更接近于強(qiáng)震的真實(shí)烈度分布?；诩扔袃x器測(cè)定記錄又有宏觀考察數(shù)據(jù)的現(xiàn)代大震建立適用于強(qiáng)震的烈度分布模型，進(jìn)而構(gòu)建強(qiáng)震參數(shù)估算方程，直接利用全部烈度點(diǎn)的空間分布信息來(lái)估算歷史強(qiáng)震的震級(jí)和震中，為歷史強(qiáng)震參數(shù)的確定提供新的思路與方法，進(jìn)一步提升地震史料的可用性，提高歷史地震參數(shù)的可靠性，對(duì)于6.5級(jí)以上歷史大震參數(shù)的校核具有重要意義。

1 基于離散烈度點(diǎn)考慮斷層破裂長(zhǎng)度的強(qiáng)震烈度分布模型

1.1 地表破裂長(zhǎng)度與震級(jí)的經(jīng)驗(yàn)關(guān)系式

Tocher（1958）是第一個(gè)把地震震級(jí)與相應(yīng)的斷層地表破裂長(zhǎng)度聯(lián)系起來(lái)的。他根據(jù)美國(guó)加利福利亞州及內(nèi)華達(dá)州10次地震的數(shù)據(jù)，利用最小二乘法建立了形如M＝a+b log L的關(guān)系式，式中M為震級(jí)，L為地震引起的地表破裂長(zhǎng)度（km）。后來(lái)，不少人相繼提出震級(jí)與地震引起的地表破裂長(zhǎng)度的類似關(guān)系式。陳達(dá)生（1984）對(duì)這類關(guān)系式做過(guò)系統(tǒng)介紹，并指出M＝a+b log L與log L＝c+dM這2個(gè)回歸方程之間不是反函數(shù)關(guān)系。根據(jù)回歸分析的定義，前者適用于估算已知斷層地表破裂長(zhǎng)度對(duì)應(yīng)的地震震級(jí)；而后者適宜于估算可能與某地震震級(jí)伴生的地表破裂長(zhǎng)度。

學(xué)者們?cè)⒘瞬簧倨屏验L(zhǎng)度與震級(jí)的經(jīng)驗(yàn)關(guān)系，但或多或少都用到了歷史地震數(shù)據(jù)或者跨國(guó)跨區(qū)域的地震資料。本文旨在基于既有儀器測(cè)定記錄又有宏觀考察數(shù)據(jù)的現(xiàn)代大震建立適用于我國(guó)的強(qiáng)震烈度分布模型，因此收集新的現(xiàn)代大震資料重新進(jìn)行經(jīng)驗(yàn)關(guān)系的統(tǒng)計(jì)擬合。自20世紀(jì)60年代以來(lái)，中國(guó)相繼發(fā)生了好幾次伴有斷層的地表破裂的大地震，我們收集了這些大地震數(shù)據(jù)并按照l(shuí)og L＝c+dM進(jìn)行了統(tǒng)計(jì)回歸，以估算可能與某地震震級(jí)伴生的地表破裂的長(zhǎng)度（表1）。

表1 用于統(tǒng)計(jì)破裂長(zhǎng)度與震級(jí)經(jīng)驗(yàn)關(guān)系的地震列表

將表1中所列地震的破裂長(zhǎng)度與震級(jí)進(jìn)行最小二乘統(tǒng)計(jì)回歸，得到

據(jù)此，在已知震級(jí)時(shí)可估算伴隨該地震的地表破裂長(zhǎng)度。

1.2 強(qiáng)震烈度分布模型及其回歸方法

地震烈度在地面的空間分布會(huì)呈現(xiàn)出一定的幾何圖像，依據(jù)空間統(tǒng)計(jì)分析的理論，可以建立一系列的地震烈度分布模型。

地震烈度分布模型描述烈度隨震級(jí)（或震中烈度）和距離而發(fā)生的變化，即假若已知某次地震的震級(jí)，應(yīng)能給出各地表點(diǎn)的烈度。目前主要的地震烈度分布模型有圓模型、橢圓模型和斷層破裂模型。對(duì)于震級(jí)較小的地震，常?？床怀龅日鹁€的方向，因此常常不考慮烈度分布的方向性而尋求一種平均關(guān)系，即圓模型。當(dāng)震級(jí)較大時(shí)，等震線形狀常呈現(xiàn)一種橢圓形，可考慮在長(zhǎng)軸和短軸方向有不同的分布，即橢圓模型。對(duì)于發(fā)震斷層出露地表的地震，高烈度區(qū)通常沿發(fā)震斷層分布為帶狀，沿?cái)鄬哟怪狈较蛩p很快，為此研究者們提出了斷層破裂模型。

對(duì)于6.5級(jí)以上的強(qiáng)震，比較適用于斷層破裂模型，但是帶狀分布的模型較為復(fù)雜，在數(shù)學(xué)上不易處理，這里我們提出基于烈度數(shù)據(jù)點(diǎn)的考慮斷層破裂長(zhǎng)度的烈度橢圓分布模型。構(gòu)建該模型的步驟為：

①由1.1節(jié)擬合得到的破裂長(zhǎng)度-震級(jí)關(guān)系式（1）估算地震可能伴隨的地表破裂長(zhǎng)度；

②直接從原始烈度點(diǎn)的空間展布出發(fā)，用最小二乘法對(duì)地震各烈度區(qū)原始烈度點(diǎn)空間分布進(jìn)行橢圓擬合，得到各烈度區(qū)烈度點(diǎn)空間展布的橢圓擬合線，賦予這條橢圓擬合線該烈度區(qū)的烈度值，本文稱之為烈度估計(jì)等值線。要求最內(nèi)圈烈度估計(jì)等值線橢圓長(zhǎng)軸要大于由公式（1）估計(jì)的地表破裂長(zhǎng)度。

③以烈度估計(jì)等值線的長(zhǎng)短半軸長(zhǎng)度為基礎(chǔ)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)回歸得到烈度橢圓分布模型。

此模型是對(duì)原始烈度點(diǎn)空間展布的綜合估計(jì)，充分利用了所有烈度點(diǎn)的原始空間分布信息，更接近真實(shí)地震的烈度分布。

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

1.2.1 強(qiáng)震烈度估計(jì)等值線擬合

中國(guó)目前公開(kāi)發(fā)表的烈度調(diào)查資料大多采用等震線的形式發(fā)布，原始烈度調(diào)查點(diǎn)通常難以獲得。本文選取已經(jīng)出版的《中國(guó)震例（1966—2002）》（張肇誠(chéng)等，1988、1990a、1990b、1999、2000；陳棋福等，2002a、2002b、2003、2008）中 2008年以后的強(qiáng)震烈度資料從中國(guó)地震局官方網(wǎng)站的地震專題上獲?。╤ttp：／／www.cea.gov.cn／publish／dizhenj／468／553／index.htm l），收集了20世紀(jì)60年代以來(lái)全國(guó)所發(fā)生的M≥6.5的地震等震線圖，提取了其中既有儀器測(cè)量數(shù)據(jù)又有宏觀烈度調(diào)查數(shù)據(jù)的現(xiàn)代地震資料，并將烈度分布圖數(shù)字化。這些地震的等震線圖上部分明確標(biāo)示了宏觀烈度調(diào)查點(diǎn)位置；對(duì)于未標(biāo)示烈度調(diào)查點(diǎn)的，本文將能在等震線圖上明確獲取烈度信息的城、縣、鎮(zhèn)、村所在地作為烈度點(diǎn)。經(jīng)數(shù)字化配準(zhǔn)后獲取各地震烈度點(diǎn)的空間分布（表2）。

考慮中心點(diǎn)和方向性的橢圓參數(shù)方程為

式中，（xi，yi）為橢圓上點(diǎn)的坐標(biāo)，（x0，y0）為橢圓中心坐標(biāo)；θ（0°≤θ≤180°）是橢圓長(zhǎng)軸逆時(shí)針與x軸正方向的夾角；Ra、Rb分別為橢圓長(zhǎng)半軸和短半軸長(zhǎng)度，ti（0≤θ≤2π）為橢圓上各點(diǎn)所對(duì)應(yīng)的參數(shù)。

表2 地震數(shù)據(jù)及各烈度區(qū)橢圓烈度估計(jì)等值線長(zhǎng)、短半軸長(zhǎng)度和橢圓中心

由式（1），可以估計(jì)相應(yīng)震級(jí)的地震可能伴隨的地表破裂長(zhǎng)度L，以此約束橢圓強(qiáng)震烈度分布模型，要求最內(nèi)圈烈度估計(jì)等值線橢圓長(zhǎng)軸要大于由式（1）估計(jì)的地表破裂長(zhǎng)度，即

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

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

圖1是 2014年10月7日云南普洱景谷6.6級(jí)地震烈度估計(jì)等值線擬合結(jié)果，以宏觀震中為平面坐標(biāo)原點(diǎn)。圖中小三角是烈度數(shù)據(jù)點(diǎn)，橢圓長(zhǎng)軸方向?yàn)闃O震區(qū)走向，烈度估計(jì)等值線是對(duì)原始烈度點(diǎn)空間分布進(jìn)行最小二乘擬合所得，可以看到烈度估計(jì)等值線趨向于原始烈度點(diǎn)空間分布的平均估計(jì)。

圖1 2014年10月7日云南普洱景谷6.6級(jí)地震烈度估計(jì)等值線分布圖

1.2.2 強(qiáng)震烈度分布模型

本文據(jù)陳達(dá)生等（1989）關(guān)于地震烈度橢圓衰減關(guān)系長(zhǎng)短軸統(tǒng)一回歸的思想，采用烈度分布模型

式中，I為地震烈度，M為震級(jí)；系數(shù) a、b、C1、C2均為回歸常數(shù)；R0a、R0b分別為橢圓長(zhǎng)、短軸兩方向烈度近場(chǎng)飽和因子；ε為回歸分析中表示不確定性的隨機(jī)變量，通常假定為對(duì)數(shù)正態(tài)分布，其均值為0，標(biāo)準(zhǔn)差為σ。該模型在震中處的烈度值隨震級(jí)的變化率為常數(shù)b，而中間距離仍保持長(zhǎng)短軸烈度的差別，同時(shí)在遠(yuǎn)場(chǎng)也使烈度分布成圓形。需要強(qiáng)調(diào)指出的是，這里Ra、Rb分別是烈度I的烈度估計(jì)等值線的長(zhǎng)、短半軸長(zhǎng)度，不再是等震線的長(zhǎng)、短半軸長(zhǎng)度。

將表2中各地震各烈度區(qū)的烈度估計(jì)等值線，按照式（4）進(jìn)行最小二乘統(tǒng)計(jì)回歸。由于擬合的烈度估計(jì)等值線已經(jīng)是對(duì)烈度點(diǎn)分布的綜合估計(jì)，所以沒(méi)有采取近場(chǎng)補(bǔ)點(diǎn)。但為了體現(xiàn)遠(yuǎn)場(chǎng)區(qū)發(fā)震構(gòu)造影響消失，烈度分布趨于圓形的特點(diǎn)，汪素云等（1993）提出取有感范圍的半徑作為遠(yuǎn)場(chǎng)控制點(diǎn)，有感烈度值通常為Ⅲ～Ⅳ度，在計(jì)算中取為Ⅲ度半（3.5度），稱為遠(yuǎn)場(chǎng)補(bǔ)點(diǎn)。這里提出的有感半徑有外包線的性質(zhì)，而本文烈度模型主要是對(duì)烈度分布的綜合估計(jì)，因此計(jì)算時(shí)遠(yuǎn)場(chǎng)控制烈度取為Ⅲ度。有感半徑與震級(jí)的關(guān)系見(jiàn)表3。

表3 有感半徑與震級(jí)的關(guān)系（汪素云等，1993）

由原始烈度點(diǎn)擬合得到烈度區(qū)的估計(jì)等值線，進(jìn)而統(tǒng)計(jì)回歸得到地震烈度分布模型為

2 地震參數(shù)估計(jì)方法及其不確定性分析

2.1 地震參數(shù)估計(jì)方法

本文在上節(jié)中建立了如下烈度分布橢圓模型：

將（2）式改寫(xiě)成橢圓標(biāo)準(zhǔn)數(shù)學(xué)方程式（即（9）式）

用以估計(jì)地震震級(jí)與宏觀震中。式中，（xi，yi）為橢圓上點(diǎn)的坐標(biāo)，（x0，y0）為橢圓中心坐標(biāo)；θ（0°≤ θ≤180°）是橢圓長(zhǎng)軸逆時(shí)針與x軸正方向的夾角。

已知烈度數(shù)據(jù)點(diǎn)平面坐標(biāo)（xi，yi）和烈度值Ii，聯(lián)立烈度分布模型及考慮中心點(diǎn)和方向性的地震參數(shù)估計(jì)方程（7）、（8）、（9），即可計(jì)算橢圓中心坐標(biāo)即地震震中平面坐標(biāo)（x0，y0）和震級(jí)M，以及橢圓長(zhǎng)軸與x軸夾角θ。

由式（7）得

由式（8）得

將式（10）和（11）帶入到式（9）中得

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

將式（13）帶入式（12）得

2.2 求解參數(shù)的基本條件

2.2.1 烈度點(diǎn)數(shù)量

在式（14）中，C1a、C1b、C2、C3a、C3b、C0a、C0b由烈度分布模型統(tǒng)計(jì)回歸可得，見(jiàn)式（5）和（6）。若已知一個(gè)地震的多個(gè)烈度信息點(diǎn)（xi，yi，Ii），帶入到（3.14）式，即可聯(lián)立方程組求得震中（x0，y0）、震級(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.2.2 烈度點(diǎn)空間分布

參數(shù)估計(jì)方程是以橢圓數(shù)學(xué)方程為基礎(chǔ)，因此用于計(jì)算的烈度點(diǎn)不能分布于一條或接近于一條的直線上，否則在極端情況下將造成方程無(wú)解，或所得數(shù)值解不滿足地震參數(shù)的物理意義。

2.2.3 適用范圍

由于在建模過(guò)程中，8.0級(jí)以上特大地震的資料太少，因此本方法只適用于6.5～8.0級(jí)的歷史強(qiáng)震，對(duì)于8.0級(jí)以上的特大地震需要另行討論。

2.3 震例驗(yàn)算

2.3.1 內(nèi)符檢驗(yàn)

帶入由1.2.2節(jié)得到的強(qiáng)震烈度分布模型系數(shù)，式（14）即變?yōu)?/p>

為了驗(yàn)證此算法的可行性，以參與強(qiáng)震烈度分布模型擬合的其中7個(gè)地震為例，將各自的烈度數(shù)據(jù)點(diǎn)（xi，yi，Ii）直接帶入到式（15）中進(jìn)行試算，計(jì)算時(shí)所有烈度點(diǎn)以各自地震的宏觀震中為原點(diǎn)（0，0）轉(zhuǎn)換到平面坐標(biāo)下，計(jì)算結(jié)果見(jiàn)表4。

表4 內(nèi)符驗(yàn)算計(jì)算結(jié)果

由表4可以看到，參與驗(yàn)算的7個(gè)地震，計(jì)算震級(jí)與儀器震級(jí)的差值不超過(guò)0.5級(jí)，計(jì)算震中到宏觀震中的距離不超過(guò)25km，由此可見(jiàn)此方法對(duì)參與模型擬合的地震集是可行的。

2.3.2 外推檢驗(yàn)

為了驗(yàn)證地震參數(shù)估計(jì)方程的外推有效性，我們提取了3個(gè)沒(méi)有參與模型擬合的強(qiáng)震烈度數(shù)據(jù)點(diǎn)，將所有烈度數(shù)據(jù)點(diǎn)以各自地震的宏觀震中為原點(diǎn)（0，0）轉(zhuǎn)換到平面坐標(biāo)下（xi，yi，Ii），然后直接帶入到式（15）中進(jìn)行驗(yàn)算，計(jì)算結(jié)果見(jiàn)表5。

表5 外推驗(yàn)算計(jì)算結(jié)果

由表5同樣可以看到，參與驗(yàn)算的3個(gè)地震，計(jì)算震級(jí)與儀器震級(jí)的差值不超過(guò)0.7級(jí)，計(jì)算震中到宏觀震中的距離不超過(guò)15km，由此可見(jiàn)此方法外推也是可行的。

2.3 不確定性分析

在震中位置和震級(jí)的估算過(guò)程中，烈度點(diǎn)的數(shù)量、分布和精度好壞將直接影響估算結(jié)果。由于歷史地震烈度點(diǎn)相對(duì)較少，為了對(duì)歷史地震參數(shù)進(jìn)行較好的估計(jì)，本文運(yùn)用蒙特卡洛方法，對(duì)信息豐富的現(xiàn)代地震烈度點(diǎn)抽樣，模擬歷史地震烈度點(diǎn)分布情況，進(jìn)而分析本文所建方法的不確定性。

由前文知，至少需要3個(gè)烈度點(diǎn)才能確定一次地震的震級(jí)、震中和方向。采用蒙特卡洛方法，先從同一地震的烈度數(shù)據(jù)點(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ī)抽樣中的一種情況。通過(guò)蒙特卡洛方法模擬出大量烈度點(diǎn)組樣本帶入?yún)?shù)估計(jì)方程，完成計(jì)算后分析此方法計(jì)算地震震級(jí)與震中的不確定性?？梢酝茰y(cè)，烈度點(diǎn)數(shù)越多，分布越均勻，得出的震中與震級(jí)越精確，不確定性越小。

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

本文以2014年10月7日云南普洱景谷6.6級(jí)地震為例，給出抽樣計(jì)算結(jié)果。首先獲取地震的烈度數(shù)據(jù)點(diǎn)，然后以其宏觀震中為原點(diǎn)轉(zhuǎn)換到平面坐標(biāo)下，再?gòu)牡卣鹆叶绕矫孀鴺?biāo)點(diǎn)里隨機(jī)抽取3點(diǎn)、4點(diǎn)、5點(diǎn)……20點(diǎn)各1000次，分別帶入式（15）進(jìn)行解算，圖2是計(jì)算所得震中的分布圖。由計(jì)算所得震級(jí)與儀器震級(jí)之差（ΔM）繪成圖3?？梢钥吹?，隨著烈度點(diǎn)的增多，計(jì)算震中逐漸向宏觀震中靠攏，計(jì)算震級(jí)與儀器震級(jí)的差值也逐漸收斂。

將參與內(nèi)符檢驗(yàn)的7個(gè)地震和參與外推檢驗(yàn)的3個(gè)地震均進(jìn)行抽樣計(jì)算，統(tǒng)計(jì)不同抽樣點(diǎn)數(shù)情況下計(jì)算震中到宏觀震中距離的均值與標(biāo)準(zhǔn)差，以及計(jì)算震級(jí)與儀器震級(jí)差值的均值與標(biāo)準(zhǔn)差（圖4）。

圖4的（a）和（b）是內(nèi)符檢驗(yàn)地震的計(jì)算結(jié)果。由圖4（a）可以看到，隨著烈度點(diǎn)的增多，計(jì)算震中到宏觀震中距離的均值和標(biāo)準(zhǔn)差都逐漸減小。取3個(gè)點(diǎn)時(shí)，距離均值最大接近26km，標(biāo)準(zhǔn)差最大接近15km；均值基本都在25km以下，取20個(gè)點(diǎn)時(shí)，標(biāo)準(zhǔn)差都在10km以下。

圖2 2014年10月7日云南普洱景谷6.6級(jí)地震抽樣計(jì)算震中分布

圖3 2014年10月7日云南普洱景谷6.6級(jí)地震計(jì)算震級(jí)與儀器震級(jí)差（ΔM）的直方圖

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

圖4的（c）和（d）是外推檢驗(yàn)地震的計(jì)算結(jié)果。由圖4（c）同樣可以看到，隨著烈度點(diǎn)的增多，計(jì)算震中到宏觀震中距離的均值和標(biāo)準(zhǔn)差都逐漸減小。取3個(gè)點(diǎn)時(shí)，距離均值最大不足25km，標(biāo)準(zhǔn)差最大接近11km；取6點(diǎn)個(gè)時(shí)，標(biāo)準(zhǔn)差在10km以下。

圖4（d）與圖4（b）的趨勢(shì)類似，可以看到隨著烈度點(diǎn)的增加，震級(jí)差值的均值有個(gè)增大過(guò)程，達(dá)到一定點(diǎn)數(shù)之后趨于平穩(wěn)。隨著烈度點(diǎn)的增加，震級(jí)差值的標(biāo)準(zhǔn)差整體上呈減小趨勢(shì)。對(duì)于參與試算的3個(gè)地震，取3個(gè)點(diǎn)時(shí)，震級(jí)差值的標(biāo)準(zhǔn)差最大接近0.45級(jí)；取7個(gè)點(diǎn)時(shí)標(biāo)準(zhǔn)差都小于0.4級(jí)。

度量不確定性的指標(biāo)就是不確定度。綜合不確定度是隨機(jī)誤差和系統(tǒng)誤差的綜合度量，通常用均方誤差來(lái)表示，其表達(dá)式為

圖4（a） 內(nèi)符檢驗(yàn)地震不同烈度點(diǎn)數(shù)情況下計(jì)算震中到宏觀震中的距離（ΔR）的均值及其標(biāo)準(zhǔn)差分布

圖4（b） 內(nèi)符檢驗(yàn)地震不同烈度點(diǎn)數(shù)情況下計(jì)算震級(jí)與儀器震級(jí)之差（ΔM）的均值及其標(biāo)準(zhǔn)差分布

圖4（c） 外推檢驗(yàn)地震不同烈度點(diǎn)數(shù)情況下計(jì)算震中到宏觀震中的距離（ΔR）的均值及其標(biāo)準(zhǔn)差分布

圖4（d） 外推檢驗(yàn)地震不同烈度點(diǎn)數(shù)情況下計(jì)算震級(jí)與儀器震級(jí)之差（ΔM）的均值及其標(biāo)準(zhǔn)差分布

式中，σ為精度，代表誤差分布的密集或離散程度，即觀測(cè)結(jié)果與其數(shù)學(xué)期望的接近程度；ε為準(zhǔn)確度，代表觀測(cè)量的真值與數(shù)學(xué)期望之差，也就是系統(tǒng)誤差；綜合不確定度D代表的是精確度，是精度和準(zhǔn)確度的合成，即觀測(cè)值與真值的接近程度，包括觀測(cè)結(jié)果與其數(shù)學(xué)期望的接近程度和數(shù)學(xué)期望與其真值的偏差。本文以綜合不確定度來(lái)度量不確定性。

對(duì)于既有儀器記錄又有宏觀考察的現(xiàn)代地震，本文將儀器震級(jí)和考察所得的宏觀震中作為真值，以均方誤差來(lái)估計(jì)蒙特卡洛方法所得計(jì)算結(jié)果的不確定度。由于本文運(yùn)用蒙特卡洛方法時(shí)采取的是隨機(jī)抽樣，因此烈度點(diǎn)的組合難免會(huì)抽到極端情況，使得所得結(jié)果不滿足地震參數(shù)的物理意義，這里按照粗差處理的方式予以剔除。然后對(duì)處理后的結(jié)果估計(jì)其綜合不確定度，結(jié)果見(jiàn)圖5和6。

圖5 內(nèi)符檢驗(yàn)地震參數(shù)估計(jì)的綜合不確定度

圖5和6顯示的趨勢(shì)大致相同。對(duì)于地震震中綜合不確定度，除個(gè)別情況外，所有點(diǎn)數(shù)情況下震中不確定度都在25km以內(nèi)。歷史地震震中精度分為1類≤10km，2類≤25km，3類≤50km，4類≤100km，5類＞100km。對(duì)應(yīng)到歷史地震震中精度，按照本文算法得到的地震震中可以直接給定為2類精度。由于震中位置主要受烈度點(diǎn)分布空間位置的限制，約束性很高，因此所得估算結(jié)果很好。

圖6 外推檢驗(yàn)地震參數(shù)估計(jì)的綜合不確定度

對(duì)于震級(jí)綜合不確定度，隨著烈度點(diǎn)數(shù)的增多，一部分地震的綜合不確定度越來(lái)越小，一部分地震的綜合不確定度先減小后增大，另一部分地震的綜合不確定度則越來(lái)越大。綜合不確定度包含了隨機(jī)誤差和系統(tǒng)誤差。隨著烈度點(diǎn)數(shù)的增多隨機(jī)誤差越來(lái)越小，因此影響綜合不確定度大小的主要是系統(tǒng)誤差。由于所建模型是從很多地震里提取的平均規(guī)律，具有離散性，與任何一個(gè)真實(shí)地震之間都存在著系統(tǒng)偏差。由圖4（b）、4（d）知道地震震級(jí)的估算隨著烈度點(diǎn)數(shù)的增多有一個(gè)先增大后平穩(wěn)的過(guò)程。當(dāng)用最少的烈度點(diǎn)計(jì)算出來(lái)的震級(jí)值都比儀器震級(jí)值大時(shí)，隨著烈度點(diǎn)數(shù)的增多，系統(tǒng)偏差會(huì)越來(lái)越大，導(dǎo)致綜合不確定度隨著烈度點(diǎn)數(shù)的增加而增大；當(dāng)用最多的烈度點(diǎn)計(jì)算出來(lái)的震級(jí)值都比儀器震級(jí)值小時(shí)，隨著烈度點(diǎn)數(shù)的增多，計(jì)算震級(jí)會(huì)越來(lái)越靠近儀器震級(jí)，顯示為綜合不確定度隨著烈度點(diǎn)數(shù)的增多而減小；而用最少烈度點(diǎn)計(jì)算出來(lái)的震級(jí)比儀器震級(jí)小，最多烈度點(diǎn)計(jì)算出來(lái)的震級(jí)比儀器震級(jí)大時(shí)，系統(tǒng)偏差有一個(gè)先減小后增大的過(guò)程，就導(dǎo)致了綜合不確定度隨著烈度點(diǎn)數(shù)的增多先減小后增大。

造成震級(jí)不確定度無(wú)法控制的主要原因有：①所建模型時(shí)根據(jù)的是一個(gè)地區(qū)眾多地震中的平均規(guī)律，模型回歸需要有足夠的地震樣本，因此選定的研究區(qū)域比較大；而研究區(qū)域越大，則場(chǎng)地條件等因素就越偏離某個(gè)具體地震的實(shí)際情況。為了使模型更能適應(yīng)具體地震情況，就需要縮小研究區(qū)域，但代價(jià)是研究區(qū)域內(nèi)發(fā)生的地震太少，造成模型回歸的離散性變大。兩相鉗制，較難取舍。②由烈度推算震級(jí)，震源深度等參數(shù)對(duì)結(jié)果影響較大。本文所建模型里暫時(shí)沒(méi)有考慮震源深度等其他因素的影響，但每個(gè)具體地震的深度等參數(shù)仍然是有差別的，造成具體地震對(duì)所建模型而言各有偏差。

對(duì)地震震級(jí)的估算歷來(lái)都是一個(gè)難題，本文算法能同時(shí)估算出地震的震級(jí)和震中。由綜合不確定度來(lái)看，震中的精度可以定為2類，但是震級(jí)的精度很難直接給出。對(duì)于參與驗(yàn)算的地震，隨著烈度點(diǎn)數(shù)的增多，震級(jí)值趨于平穩(wěn)時(shí)，不確定度都不超過(guò)0.8級(jí)，多數(shù)在0.5級(jí)以下

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

3 歷史大震參數(shù)估計(jì)

通過(guò)前面的分析，本文建立了一套基于烈度數(shù)據(jù)點(diǎn)考慮斷層破裂長(zhǎng)度的烈度分布橢圓模型以估算地震震中與震級(jí)的方法，我們以此來(lái)估算華北地區(qū)幾個(gè)6.5～8.0級(jí)歷史大震，結(jié)果見(jiàn)表6。歷史大震的烈度數(shù)據(jù)點(diǎn)均從《中國(guó)歷史地震圖集》（明、清時(shí)期）（國(guó)家地震局地球物理研究所等，1986、1990）中提取。由表6可以看到，計(jì)算震級(jí)與圖集震級(jí)相差均在1級(jí)以內(nèi)，計(jì)算震中與圖集震中相差不超過(guò)25km，對(duì)應(yīng)于歷史地震2類震中精度。對(duì)于歷史大震，這樣的計(jì)算結(jié)果比較可靠。

表6 華北地區(qū)歷史大震參數(shù)估算結(jié)果

本文用5個(gè)華北地區(qū)歷史強(qiáng)震證明了所建算法的可行性。由以上震例可見(jiàn)，此方法對(duì)歷史強(qiáng)震參數(shù)的估定頗為有效。處理過(guò)程直接明了，減少了等震線勾畫(huà)過(guò)程中的主觀不確定性。

4 結(jié)論與討論

本文收集現(xiàn)代大震資料，統(tǒng)計(jì)回歸了斷層破裂長(zhǎng)度-震級(jí)經(jīng)驗(yàn)關(guān)系；針對(duì)6.5級(jí)以上強(qiáng)震烈度分布特征，提出了基于烈度數(shù)據(jù)點(diǎn)的考慮斷層破裂長(zhǎng)度的烈度分布橢圓模型，在此基礎(chǔ)上聯(lián)立橢圓數(shù)學(xué)方程建立了強(qiáng)震參數(shù)估計(jì)方法，通過(guò)內(nèi)符檢驗(yàn)和外推檢驗(yàn)證明了算法的可行性。采用蒙特卡洛方法，通過(guò)大量的抽樣模擬，定量分析了所得參數(shù)的不確定性，結(jié)果表明：計(jì)算震級(jí)精度在1級(jí)以內(nèi)，計(jì)算震中精度為歷史地震的2類精度。進(jìn)而選取了華北地區(qū)5個(gè)歷史強(qiáng)震進(jìn)行了地震參數(shù)估算。本方法直接通過(guò)烈度數(shù)據(jù)點(diǎn)來(lái)估算歷史強(qiáng)震的震級(jí)與震中，回避了對(duì)等震線的勾畫(huà)，處理過(guò)程直接明了，一定程度上減少了主觀不確定性，提高了科學(xué)性。由于在建模過(guò)程中，8.0級(jí)以上特大地震的資料太少，因此本方法只適用于6.5～8.0級(jí)的歷史強(qiáng)震，對(duì)于8.0級(jí)以上的特大地震需要另行討論。

不論是破裂長(zhǎng)度與震級(jí)的經(jīng)驗(yàn)關(guān)系，還是強(qiáng)震烈度分布模型都應(yīng)該具有區(qū)域性，應(yīng)該根據(jù)研究區(qū)域的地震相關(guān)資料統(tǒng)計(jì)得到。但在一個(gè)研究區(qū)內(nèi)很難找到足夠多的現(xiàn)代大震統(tǒng)計(jì)資料，所以不得不采用全國(guó)的強(qiáng)震數(shù)據(jù)建立經(jīng)驗(yàn)公式，那么應(yīng)用到某個(gè)研究區(qū)時(shí)必然會(huì)產(chǎn)生某些偏差。地震烈度分布模型是據(jù)多個(gè)地震的平均規(guī)律而建立的，與地震資料之間存在著離散性，與任何一個(gè)真實(shí)地震之間都有系統(tǒng)偏差。模型是原型的抽象和簡(jiǎn)化表征，建模本身是一個(gè)持續(xù)的永無(wú)止境的活動(dòng)集合。隨著技術(shù)的發(fā)展、科學(xué)的進(jìn)步以及基礎(chǔ)數(shù)據(jù)的積累，需要對(duì)模型不斷的改進(jìn)以期更接近真實(shí)情況。

筆者認(rèn)為，今后應(yīng)多引入現(xiàn)代地震學(xué)的方法來(lái)處理歷史地震，例如早期國(guó)外臺(tái)站的波形數(shù)據(jù)、震源機(jī)制等，以提高歷史地震參數(shù)的可靠性

致謝：衷心感謝金嚴(yán)研究員、汪素云研究員、俞言祥研究員和潘華研究員在本文研究過(guò)程中給出的有益指導(dǎo)與建議！