古登堡-里希特定律中的b值統(tǒng)計樣本量研究1

李世杰 呂悅軍 劉靜偉

?

古登堡-里希特定律中的值統(tǒng)計樣本量研究1

李世杰 呂悅軍 劉靜偉

（中國地震局地殼應(yīng)力研究所，北京 100085）

值是研究地震活動的重要指標，其廣泛應(yīng)用于地震危險性分析和地震預(yù)測研究之中，與實際資料的完整性、樣本量的大小、計算方法等因素有著重要的關(guān)系。常見的值計算方法有最小二乘法和最大似然法，樣本量的大小對這2種方法影響很大。本文利用蒙特卡羅模擬地震目錄和汾渭地震帶實際目錄作為樣本，從中抽取不同大小的樣本量進行計算，研究不同樣本量下這2種方法計算得到的值與設(shè)定值或真實值之間的差別。結(jié)果表明，最小二乘法需要的最低樣本量為1000，最大似然法為200；當樣本量達不到要求時，計算出的值是不可靠的；由于對樣本量的要求不同，前者適用于計算區(qū)域的整體值，而后者在研究某區(qū)域值在時間軸上的變化方面更有優(yōu)勢。本研究為確定2種值計算方法對樣本量的最低要求提供了參考依據(jù)。

值 最小二乘法 最大似然法 樣本量大小

引言

古登堡和里希特在1941年美國地質(zhì)學(xué)會會刊中提出全球的地震活動服從經(jīng)驗關(guān)系lg=-（式中表示震級，表示震級≥的地震次數(shù)，、是常數(shù)）。一般認為值代表地區(qū)的地震活動總體水平（Rundle，1989）；值代表地震活動的大小地震數(shù)量的比例，是地震活動研究中的重要參量。研究表明，值具有明確的物理意義，與地殼的介質(zhì)特性、應(yīng)力狀態(tài)和不均勻性有關(guān)，能反映所研究區(qū)域的地質(zhì)構(gòu)造特征（王熠熙等，2015；謝卓娟等，2015）以及地震的震源特征（Schorlemmer等，2005；Gulia等，2010；劉靜偉等，2016；張廣偉，2016）。因此，值廣泛應(yīng)用于地震危險性分析和地震預(yù)測研究之中。在地震危險性分析中，值和地震年發(fā)生率共同用于確定地震活動的水平（胡聿賢，1999），其取值對地震危險性分析結(jié)果的影響較大（鄢家全等，1996；黃瑋瓊等，1998；謝卓娟等，2013）；而在地震預(yù)測研究中，值作為基本的地震活動性參數(shù)，成為地震預(yù)測的常用指標參數(shù)（韓渭賓，2003；沈建文等，2007）。

通常根據(jù)實際地震資料統(tǒng)計得到值。目前，常用的值估算方法是最大似然估計法（Aki，1965）和最小二乘法等，這些方法不僅對實際資料的完整性和精度有一定要求，同時也需要足夠的地震資料。在實際工作中，對于歷史地震資料短缺或地震活動水平低的地區(qū)，在計算時常常將現(xiàn)代小震資料與歷史地震資料聯(lián)合使用，以彌補地震資料樣本量的不足（黃瑋瓊等，1989；鄢家全等，1996；胡聿賢，1999；潘華等，2006）。統(tǒng)計值對地震資料樣本量需求的定量研究，國內(nèi)尚無專門的研究報道。在國外，Nava等（2017）利用蒙特卡羅模擬地震目錄進行抽樣估計，研究最大似然法計算值時對樣本量的需求，得出計算值時樣本量和精度之間的相互關(guān)聯(lián)關(guān)系，但沒有應(yīng)用實際地震目錄進行分析。國內(nèi)僅有少部分學(xué)者的研究涉及相關(guān)內(nèi)容，如韓曉明等（2016）研究河套地震帶值時空變化特征的文章中，探討了地震前后值的變化規(guī)律，分別使用最小二乘法和最大似然法對值時間和空間進行了掃描計算，在最小二乘法計算中，設(shè)定每次計算的窗長內(nèi)包含的樣本數(shù)目不少于100，而最大似然法的掃描窗內(nèi)包含的樣本數(shù)目不少于20；劉方斌等（2017）在魯西南聊考斷裂帶地震危險性評價與活動性分布的研究中，利用最大似然法和最小二乘法估算值并進行了對比研究，但是沒有給出樣本量的具體數(shù)目，且沒有用分震級段的方法進行值估算。

本文采用Utsu（1965）提出的最大似然法和最小二乘法，利用模擬地震目錄和實際地震目錄，定量分析最小二乘法和最大似然法計算值時分別對地震資料樣本量的需求。

1 資料概況

1.1 模擬地震目錄

蒙特卡羅法是以抽樣和隨機數(shù)的產(chǎn)生為基礎(chǔ)的隨機性方法，也稱為隨機抽樣法、計算機隨機模擬法等。蒙特卡羅方法的基本原理是通過數(shù)字模擬試驗，得到所要求解的出現(xiàn)某種事件的概率作為問題的近似解。其基本思想是：為了求解數(shù)學(xué)、物理、工程技術(shù)以及管理等方面的問題，首先建立概率模型或隨機過程，使用相應(yīng)的參數(shù)，得到某些問題（如概率分布或數(shù)學(xué)期望等問題）的解；然后通過對模型或過程的觀察或抽樣試驗來計算所求參數(shù)的統(tǒng)計特征，并用算術(shù)平均值作為所求解的近似值。因此，只要盡可能完整地表示地震震源模型，即獲得某一區(qū)域內(nèi)地震事件的發(fā)生時間和空間上的分布規(guī)律，那么就可以直接利用蒙特卡羅方法產(chǎn)生合成地震序列。

蒙特卡羅方法的基本依據(jù)是隨機模擬次數(shù)足夠多的情況下，事件發(fā)生的頻率可以反映事件發(fā)生的概率。當模擬次數(shù)有限時，計算結(jié)果和真實值之間必然存在誤差，這種誤差隨著模擬次數(shù)的增加而減少。

蒙特卡羅模擬中隨機變量的簡單子樣（1，2，3，……，X）是獨立分布的，即每次模擬中事件發(fā)生的次數(shù)與其它任意1次中事件發(fā)生的次數(shù)無關(guān)，那么隨機變量就是服從泊松分布的，當隨機變量的期望值趨向于無窮時，泊松分布趨近于正態(tài)分布。

本文利用蒙特卡羅方法模擬地震目錄（張建中，1974a，1974b；任雪梅等，2011），研究不同樣本量下計算方法對值的影響。

首先，構(gòu)造概率分布模型：

其中，（）為震級的概率密度函數(shù)，表示為：

其次，將公式（2）帶入公式（1），得到關(guān)于模擬震級的函數(shù)：

設(shè)定=1、0=2.0、u=8.0，利用MATLAB隨機數(shù)程序模擬生成不同樣本量的地震目錄，樣本數(shù)分別為5000、4500、4000、3500、3000、2500、2000、1750、1500、1250、1000、750、500、250、100、80、60、40、20、10，每個目錄模擬10000組。由于實際大地震非常稀缺，所以震級上限的限制對計算結(jié)果沒有顯著影響。

1.2 汾渭地震帶實際地震目錄

本研究使用的實際地震資料來自中國地震臺網(wǎng)中心的地震數(shù)據(jù)庫。研究表明地震資料的完整性對地震活動分析非常重要（Wiemer等，2000；Rotondi等，2002）。經(jīng)統(tǒng)計分析（焦遠碧等，1990；黃偉瓊等，1994a，1994b；周公威等，2007；謝卓娟等，2012a；徐偉進等，2014），1970年以來汾渭地震帶2.0以上的地震資料基本完整，故本研究選取1970年1月—2010年12月，≥2.0的8184次地震資料進行研究分析，其中2.0—2.9地震7044次、3.0—3.9地震964次、4.0—4.9地震148次、5.0—5.9地震25次、6.0—6.9地震3次，所用的地震目錄均已刪除前震、余震。汾渭地震帶4.7級以上地震的震中分布見圖1。

圖1 汾渭地震帶范圍及震中分布

2 模擬地震目錄的結(jié)果分析

利用最小二乘法和最大似然法對模擬的地震目錄進行定量研究分析，比較2種方法在不同的樣本量下計算得到的平均值以及樣本量的大小對計算結(jié)果的影響。

文中使用的最大似然法為Utsu（1965）提出的公式：

采用上述最大似然法和最小二乘法計算，分別得到平均值隨樣本量的變化，如圖2、3所示。

圖2 不同樣本量下最大似然法計算得到的平均b值擬合圖

圖3 不同樣本量下最小二乘法計算得到的平均b值擬合圖

由圖2、3可知，當樣本量大于300時，最大似然法計算的平均值能取得符合預(yù)期的數(shù)值；樣本量大于1000時，最小二乘法計算的平均值趨于穩(wěn)定，向理論值收斂。

表1 最大似然法模擬結(jié)果

表2 最小二乘法模擬結(jié)果

當設(shè)定=1時，不同樣本下得到的值分布直方圖見圖4。由圖可知，當樣本量≥200時，最大似然法計算的值很好地向設(shè)定值收斂；樣本量≥500時，最小二乘法計算的值有較好的收斂效果。從估值概率來看，取得的值一般都偏低。

根據(jù)上述分析，2種方法雖然都能得到符合預(yù)期的數(shù)值，但最小二乘法計算的平均值通常比最大似然法的小，而標準差比最大似然法的大；樣本量小于1000時，利用最小二乘法得到的平均值精確度較低；樣本量大于1000時，2種計算方法得出的結(jié)果差別不大。進一步分析認為，雖然2種方法都需要一定的樣本量，且樣本量越大、得到的結(jié)果越準確，但最大似然法對樣本量的要求要比最小二乘法低，樣本量大于200時，計算得到的平均值與設(shè)定值一致性較好；樣本量大于1000時，基本等于設(shè)定值。

3 實際地震目錄的計算結(jié)果分析

根據(jù)汾渭地震帶地震資料的完整性和可靠性研究，自1970年以來，臺站記錄得到的≥2.0地震目錄基本完整，故本研究以汾渭地震帶1970年1月—2010年12月的地震記錄為計算樣本，研究實際地震目錄樣本量對值計算的影響。表3給出該地震帶不同震級檔的地震數(shù)目，并計算出不同震級檔的年平均發(fā)生率（黃瑋瓊等，1989；潘華等，2006；吳兆營等2005；謝卓娟等，2012b）。

圖4 不同樣本量的b值分布直方圖

表3 汾渭地震帶地震分檔統(tǒng)計和年平均發(fā)生率

續(xù)表

震級檔M地震個數(shù)年平均發(fā)生率 5.0—5.4150.36 5.5—5.9100.23 6.0—8.530.07

注：年平均發(fā)生率指震級≥的年均地震數(shù)，代表地震活動水平。

根據(jù)表3中統(tǒng)計的各震級檔的地震個數(shù)和年平均發(fā)生率，利用最小二乘法和最大似然法計算得到的值分別為0.75和0.75335。

為研究實際地震目錄下值計算對樣本量的要求，采用任雪梅（2011）的抽樣原則，對汾渭地震帶的地震目錄進行均勻抽樣分析，每次抽樣完將樣本放回進行下一次抽樣，得到樣本量分別為10、50、100、200、300、500、700、1000、5000的地震目錄，每個目錄重復(fù)抽樣10000次。同樣，采用2種方法分別計算不同樣本量下的平均值（表4），并分析不同樣本量下的值變化。

表4 汾渭地震帶b值擬合情況（1500—2010年）

由表4可以看出，2種計算方法得到值隨著樣本量的增加逐漸接近真實值；當樣本量大于200時，最大似然法能夠得出相對穩(wěn)定可靠的值；當樣本量大于500時，最小二乘法才能得到相對穩(wěn)定可靠的值。因此，對實際的地震目錄，最小二乘法對樣本量的要求也比較高，且受地震目錄完整性的影響較大。

4 結(jié)論

通過上述研究和對比分析，得到以下認識：

（1）利用最大似然法計算值時，樣本量至少要在200以上；對于最小二乘法，樣本量要求不少于1000；一般來說，計算的值都低于設(shè)定值。

（2）最小二乘法利用震級-頻度的線性關(guān)系來擬合計算值，樣本量不足會影響其線性關(guān)系。因此，最小二乘法受樣本量影響較大，樣本量小于1000時，值的計算值與設(shè)定值相差較大，數(shù)值也不穩(wěn)定。

（3）最大似然法方便快捷，受樣本量影響小，計算出的值相對穩(wěn)定，但誤差估計值偏大，其利用平均震級計算，只與地震個數(shù)有關(guān)，受數(shù)據(jù)質(zhì)量的影響較?。蛔钚《朔ㄊ軜颖玖坑绊戄^大，在樣本量充足的情況下計算的值比較準確，而在樣本量不足時值波動較大，誤差也隨著樣本量的減少而增加。在實際應(yīng)用中，以地震目錄充足為前提，研究不同區(qū)域的值優(yōu)先選用最小二乘法；而研究某區(qū)域值時間上的變化時，采用最大似然法估算的值相對穩(wěn)定，更能體現(xiàn)值長期變化的趨勢。

（4）從誤差和計算量來看，最大似然法比最小二乘法要小，但隨著樣本量的增加，2種方法計算結(jié)果的差異越來越小。但由于半對數(shù)坐標下不同震級檔數(shù)據(jù)權(quán)重的不對等性，國外已極少使用最小二乘法。

韓渭賓，2003．值在地震預(yù)測中的三類應(yīng)用及其物理基礎(chǔ)與須注意的問題．四川地震，（1）：1—5．

韓曉明，張文韜，王樹波等，2016．河套地震帶的值時空變化特征分析．中國地震，32（3）：522—532．

胡聿賢，1999．地震安全性評價技術(shù)教程．北京：地震出版社．

黃瑋瓊，時振梁，曹學(xué)鋒，1989．值統(tǒng)計中的影響因素及危險性分析中值的選取．地震學(xué)報，11（4）：351—361．

黃瑋瓊，李文香，曹學(xué)鋒，1994a．中國大陸地震資料完整性研究之一——以華北地區(qū)為例．地震學(xué)報，16（3）：273—280．

黃瑋瓊，李文香，曹學(xué)鋒，1994b．中國大陸地震資料完整性研究之二——分區(qū)地震資料基本完整的起始年分布圖象．地震學(xué)報，16（4）：423—432．

黃瑋瓊，李文香，1998．地震區(qū)劃中值統(tǒng)計時空范圍的確定．地震學(xué)報，20（5）：2—6．

焦遠碧，吳開統(tǒng)，楊滿棟，1990．我國地震臺網(wǎng)監(jiān)測能力及臺網(wǎng)觀測條件質(zhì)量評定．中國地震，6（4）：3—9．

劉方斌，曲均浩，田兆陽等，2017．魯西南聊考斷裂帶地震危險性評價與活動性分布．大地測量與地球動力學(xué)，37（8）：802—807．

劉靜偉，呂悅軍，2016．川滇地區(qū)值空間分布特征及其與震源類型關(guān)系的初步探討．震災(zāi)防御技術(shù)，11（3）：561—572．

潘華，李金臣，2006．地震統(tǒng)計區(qū)地震活動性參數(shù)值及4不確定性研究．震災(zāi)防御技術(shù)，1（3）：218—224．

任雪梅，2011．地震區(qū)劃中值統(tǒng)計的若干問題研究．北京：中國地震局地球物理研究所．

任雪梅，高孟潭，馮靜，2011．地震目錄的完整性對值計算的影響．震災(zāi)防御技術(shù)，6（3）：257—268．

沈建文，余湛，邱瑛，2007．地震安評中地震活動性的統(tǒng)計區(qū)域與值．國際地震動態(tài)，（3）：1—6．

王熠熙，張輝，劉雙慶等，2015．河北平原地震帶值時空變化特征．地震工程學(xué)報，37（01）：188—195.

吳兆營，薄景山，劉志平等，2005．東北地震區(qū)值和地震年平均發(fā)生率的統(tǒng)計分析．東北地震研究，21（3）：27—32，72．

謝卓娟，呂悅軍，張力方，等，2012a．基于現(xiàn)代地震資料確定汾渭地震帶分區(qū)及其地震活動性參數(shù)．地球物理學(xué)進展，27（3）：894—902．

謝卓娟，呂悅軍，彭艷菊等，2012b．東北地震區(qū)小震資料完整性分析及其對地震活動性參數(shù)的影響研究．中國地震，28（3）：256—265．

謝卓娟，呂悅軍，蘭景巖等，2013．值和4的統(tǒng)計分析及其不確定性對地震危險性分析結(jié)果的影響研究．地震研究，36（1）：86—92．

謝卓娟，李山有，呂悅軍，2015．滇西南地區(qū)主要活動斷裂的值空間分布特征．地球科學(xué)（中國地質(zhì)大學(xué)學(xué)報），40（10）：1755—1766.

徐偉進，高孟潭，2014．中國大陸及周緣地震目錄完整性統(tǒng)計分析．地球物理學(xué)報，57（9）：2802—2812.

鄢家全，韓煒，高孟潭，1996．地震活動性參數(shù)的不確定性及其對區(qū)劃結(jié)果的影響．中國地震，12（S1）：71—77．

張廣偉，2016．云南地區(qū)地震的重新定位及值研究．中國地震，32（1）：54—62．

張建中，1974a．蒙特卡洛方法(Ⅰ)．數(shù)學(xué)的實踐與認識，（01）：28—40.

張建中，1974b．蒙特卡洛方法(Ⅱ)．數(shù)學(xué)的實踐與認識，（02）：43—56.

周公威，張伯明，吳忠良等，2007．中國數(shù)字地震臺網(wǎng)（CDSN）的近期發(fā)展．地球物理學(xué)進展，22（4）：1130—1134．

Maximum likelihood estimate ofin the formula log=-and its confidence limits. Bulletin of the Earthquake Research Institute, 43(2): 237—239.

Gulia L., Wiemer S., 2010. The influence of tectonic regimes on the earthquake size distribution: A case study for Italy. Geophysical Research Letters, 37(10): L10305.

Nava F. A., Márquez-Ramírez V. H., Zú?iga F. R., et al., 2017. Gutenberg-Richter-value maximum likelihood estimation and sample size. Journal of Seismology, 21(1): 127—135.

Rundle J. B., 1989. Derivation of the complete Gutenberg-Richter magnitude-frequency relation using the principle of scale invariance. Journal of Geophysical Research: Solid Earth, 94(B9): 12337—12342.

Rotondi R., Garavaglia E., 2002. Statistical Analysis of the Completeness of a Seismic Catalogue. Natural Hazards, 25(3):245—258.

Schorlemmer D., Wiemer S., Wyss M., 2005. Variations in earthquake-size distribution across different stress regimes. Nature, 437(7058): 538—542.

Utsu T., 1965. A method for determining the value ofin a formula log=-showing the magnitude-frequency relation for earthquakes. Geophysical Bulletin of Hokkaido University, (13): 99—103.

Wiemer S., Wyss M., 2000. Minimum magnitude of completeness in earthquake catalogs: examples from Alaska, the Western United States, and Japan. Bulletin of the Seismological Society of America, 90(4): 859—869.

The Study of Sample Size on-value Statistics in the Gutenberg-Richter's Law

Li Shijie, Lü Yuejun and Liu Jingwei

(Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China)

The-value is an important indicator for evaluating the level of seismicity, and is widely used in seismic hazard analysis and earthquake prediction research. The-value is basically effected by some factors, such as the actual data integrity, the quantity of earthquake samples, calculation methods and so on. The most widely used two methods for calculatingvalues include the least square method (LSM) and the maximum likelihood method (MLM), and the results from both methods can remarkably vary with different the quantity of earthquake samples. In this paper, based on the real catalog of the Fen-Wei seismic zone and the Monte-Carlo simulated earthquake catalog, we use different number of samples to calculate thevalues and try to find the differences between the calculated results and the set values/true values. It turns out that the threshold value of samples is 1000 and 200 for LSM and MLM, respectively. When the sample size does not meet the requirements, the calculated-value is not reliable. It suggets that LSM is suitable for calculating the-value of the whole region, while MLM is more advantageous in studying the-value of a region on the time axis. This study provides a reference for the threshold values of the samples for two methods, which is of great significance to the seismic hazard analysis and earthquake prediction.

-value; The least square method; The maximum likelihood method; Sample size

李世杰，呂悅軍，劉靜偉，2018．古登堡-里希特定律中的值統(tǒng)計樣本量研究．震災(zāi)防御技術(shù)，13（3）：636—645．

10.11899/zzfy20180315

中國地震局地殼應(yīng)力研究所中央級公益性科研院所基本科研業(yè)務(wù)專項（ZDJ2016-03）

2018-02-28

李世杰，男，生于1991年。碩士研究生。主要從事工程地震研究。E-mail：729883299@qq.com

呂悅軍，男，生于1966年。研究員。研究領(lǐng)域：工程地震。E-mail：luyj1@263.net