基于蒙特卡洛方法的地震目錄模擬及相符性檢驗(yàn)
——以汾渭地震帶為例

邵霄怡 王曉青 竇愛(ài)霞 袁小祥

中國(guó)地震局地震預(yù)測(cè)研究所，北京市復(fù)興路63號(hào) 100036

0 引言

地震目錄是地震活動(dòng)性研究和地震危險(xiǎn)性分析的重要資料?，F(xiàn)代儀器記錄地震目錄的歷史只有幾十年，雖然基于歷史文獻(xiàn)記載和古地震研究可以將地震目錄時(shí)間延伸，但仍不足以研究復(fù)發(fā)周期長(zhǎng)、概率小的大地震事件，無(wú)論是基于統(tǒng)計(jì)學(xué)方法來(lái)重現(xiàn)周期長(zhǎng)達(dá)百年的地震預(yù)測(cè)或地震區(qū)劃研究，還是對(duì)基于一定物理模型假設(shè)的地震前兆或預(yù)測(cè)模型進(jìn)行檢驗(yàn)和修正，都是不夠的（周仕勇，2008）。因此，開展地震目錄模擬研究具有重要意義。

人工地震目錄模擬已有數(shù)十年的發(fā)展歷史，可采用基于力學(xué)實(shí)驗(yàn)或數(shù)值模擬、統(tǒng)計(jì)學(xué)，以及力學(xué)與統(tǒng)計(jì)模型相結(jié)合等方法進(jìn)行。Mogi（1962）通過(guò)巖石破裂力學(xué)實(shí)驗(yàn)，記錄實(shí)驗(yàn)過(guò)程中微觀破裂事件產(chǎn)生的聲發(fā)射，依據(jù)微破裂事件序列方法模擬地震目錄，并根據(jù)破裂事件（可包含多個(gè)“前震”或“余震”）序列，研究了主震型、前震-主震-余震型和群震型地震的破裂過(guò)程（唐春安，1997；唐春安等，1997）。這些模擬結(jié)果對(duì)于研究地震孕育、發(fā)生的成因機(jī)理具有重要啟示，但無(wú)法模擬由一系列獨(dú)立地震事件組成的長(zhǎng)地震序列目錄。

基于力學(xué)與統(tǒng)計(jì)模型相結(jié)合的地震目錄模擬，建立在地震孕育發(fā)生機(jī)制理論模型（Sólnes et al，1997；石耀霖等，1998；Console et al，2016；Zhou et al，2016）假設(shè)的基礎(chǔ)上，如基本的彈簧-滑塊模型（Burridge et al，1967）或改進(jìn)的彈簧-滑塊-阻尼器模型（朱元清等，1991），以及在此基礎(chǔ)上建立的斷層（相互作用）系統(tǒng)地震活動(dòng)模型（金欣等，2017；周仕勇，2008）。以這些力學(xué)模型所反映的應(yīng)力-應(yīng)變過(guò)程為基礎(chǔ)，結(jié)合細(xì)胞自動(dòng)機(jī)模型（劉桂萍等，1995、2000）、耦合應(yīng)力釋放統(tǒng)計(jì)模型（AIC）（李紅等，2015）、布朗時(shí)間過(guò)程模型（BPT模型）、時(shí)間可預(yù)報(bào)模型或基于經(jīng)驗(yàn)的模型（Working Group on California Earthquake Probabilities，2003；Hainzl et al，2007）、EEPAS模型（Rhoades et al，2011；畢金孟等，2017）等，可實(shí)現(xiàn)對(duì)地震目錄的模擬。彈簧-滑塊模型或斷層系統(tǒng)模型具有比較明確的物理含義和較為嚴(yán)格的理論基礎(chǔ)，但因影響地震發(fā)生的因素較多，模型參數(shù)的設(shè)置與實(shí)際斷裂構(gòu)造間存在一定差異，模型越復(fù)雜，假設(shè)條件就越多，這造成模擬的地震目錄不確定性較大，模擬結(jié)果的檢驗(yàn)存在一定的困難（李紅等，2015）。

基于統(tǒng)計(jì)學(xué)的地震目錄模擬方法主要依據(jù)地震活動(dòng)的概率分布特點(diǎn)，采用隨機(jī)實(shí)驗(yàn)方法模擬地震目錄，如泊松模型（趙宏等，2015）、對(duì)數(shù)正態(tài)模型（Nishenko et al，1987）、韋布爾分布模型（Hagiwara，1974）、伯努利模型（傅征祥等，1995）、ETAS模型（蔣長(zhǎng)勝等，2013；徐偉進(jìn)等，2017）等。這些方法沒(méi)有直接與地震的成因過(guò)程相結(jié)合，而是通過(guò)結(jié)合地震活動(dòng)期、幕等時(shí)間活動(dòng)特征或成組活動(dòng)規(guī)律、空間叢集以及活動(dòng)構(gòu)造塊體及其邊緣活動(dòng)等特征，在一定程度上貼近實(shí)際地震時(shí)空活動(dòng)規(guī)律。

蒙特卡洛方法使用隨機(jī)抽樣技術(shù)來(lái)實(shí)現(xiàn)對(duì)物理問(wèn)題的解決（Huh et al，2016），該方法較逼真地描述了具有隨機(jī)性質(zhì)事物的特點(diǎn)及物理實(shí)驗(yàn)過(guò)程。本文將采用泊松模型，基于蒙特卡洛隨機(jī)獨(dú)立重復(fù)實(shí)驗(yàn)，結(jié)合地震震級(jí)分布和時(shí)間分布規(guī)律，模擬區(qū)域地震目錄，并對(duì)生成的目錄進(jìn)行符合性檢驗(yàn)，進(jìn)而以此目錄為基礎(chǔ)，對(duì)研究區(qū)未來(lái)地震活動(dòng)危險(xiǎn)性進(jìn)行估計(jì)。

1 基于蒙特卡洛隨機(jī)獨(dú)立重復(fù)實(shí)驗(yàn)的地震目錄模擬方法

蒙特卡洛方法（Monte Carlo method），也稱統(tǒng)計(jì)模擬方法或隨機(jī)抽樣技術(shù)，是以概率和統(tǒng)計(jì)理論為基礎(chǔ)的數(shù)值計(jì)算方法。該方法具有對(duì)小概率事件處理較好，受幾何條件限制小、收斂速度與問(wèn)題的維數(shù)無(wú)關(guān)等特性，可用于地震、地震災(zāi)害分析預(yù)測(cè)等領(lǐng)域（劉善琪等，2013；廖景高等，2014；Iervolino et al，2016）。近年來(lái)，在地震活動(dòng)性（Ebel et al，1999；Musson，2000；Shaw et al，2007；Yazdani et al，2012；Assatourians et al，2013；Mohammed et al，2014；郭星，2015；郭星等，2016）、地震危險(xiǎn)性（Wang et al，2014；Pavel et al，2017；Setiawan，2017）和地震風(fēng)險(xiǎn)評(píng)估（Weatherill，2009；Yazdani et al，2012；Bourne et al，2015；劉甲美等，2016）等方面得到應(yīng)用。

1.1 蒙特卡洛隨機(jī)獨(dú)立重復(fù)實(shí)驗(yàn)方法

對(duì)于一個(gè)給定的物理（或數(shù)學(xué)）系統(tǒng)，如果系統(tǒng)中某一核心變量s與其概率分布函數(shù)F（s）的關(guān)系已知，即

函數(shù)F（s）為值域［0，1］的單調(diào)增函數(shù)。通過(guò)隨機(jī)獨(dú)立抽樣確定［0，1］區(qū)間均勻隨機(jī)數(shù)r，對(duì)式（1）進(jìn)行反變換，得

即可得到核心變量s的取值。通過(guò)多次獨(dú)立重復(fù)抽樣，可以確定一系列核心變量值，完成該系統(tǒng)參數(shù)的隨機(jī)模擬實(shí)驗(yàn)。上述方法稱為反抽樣法，是蒙特卡洛隨機(jī)實(shí)驗(yàn)的基本方法之一（Huh et al，2016）。

1.2 指定年限內(nèi)震級(jí)分布確定

假定地震的震級(jí)分布滿足截?cái)嗟腉-R關(guān)系（Cornell，1968），則地震震級(jí)分布函數(shù)F（m）（Youngs et al，1985）為

其中，β＝2.3b，b為研究區(qū)G-R關(guān)系的系數(shù)；m0為震級(jí)下限，muz為震級(jí)上限。F（m）為［0，1］區(qū)間的單調(diào)增函數(shù)，令F（m）＝r，對(duì)式（3）進(jìn)行變換，得

通過(guò)產(chǎn)生［0，1］區(qū)間的隨機(jī)數(shù)r，確定F（m）的值，再根據(jù)式（4）求得地震震級(jí)m。

1.3 指定年限內(nèi)地震數(shù)量和相鄰兩次地震時(shí)間間隔確定

假定地震活動(dòng)遵從泊松分布，則對(duì)于年平均發(fā)生率v0的區(qū)域，在［0，t］時(shí)間段內(nèi)發(fā)生k次m0級(jí)以上地震的概率（潘華等，2013）為

對(duì)式（5）進(jìn)行求和，得到地震發(fā)生次數(shù)的概率分布函數(shù)為

通過(guò)隨機(jī)抽樣，確定一個(gè)［0，1］的隨機(jī)數(shù)作為概率分布函數(shù)P（n≤k）的值，根據(jù)給定的年發(fā)生率v0和時(shí)間t，由式（6）運(yùn)用反抽法就可以確定該時(shí)段地震發(fā)生的次數(shù)k。

由于地震活動(dòng)服從泊松分布，任意2次地震的時(shí)間間隔τj＝tj-tj-1則遵循指數(shù)分布（胡聿賢，2006），其分布函數(shù)為

仿照反抽法，可以得到地震事件的時(shí)間間隔序列，以第1個(gè)地震為起點(diǎn)，根據(jù)時(shí)間間隔tj，即可得到第j次地震事件的發(fā)震時(shí)間τj。

因此，通過(guò)蒙特卡洛隨機(jī)獨(dú)立重復(fù)實(shí)驗(yàn)，依據(jù)年發(fā)生率v0和b值，可確定區(qū)域地震震級(jí)分布以及地震事件數(shù)量k。在此基礎(chǔ)上，針對(duì)每次地震事件，確定震級(jí)m和發(fā)震時(shí)間t，得到滿足泊松分布假設(shè)的模擬地震目錄（圖1）。

2 汾渭地震帶地震目錄模擬

圖1 地震目錄生成的技術(shù)流程

圖2 汾渭地震帶歷史破壞性地震震中分布（公元780～2014年，m≥4.6）

汾渭地震帶（統(tǒng)計(jì)區(qū)）位于鄂爾多斯塊體東、南緣，南起渭河盆地，北止于懷來(lái)-延慶盆地（王繼等，2002）。截至2014年，共發(fā)生5級(jí)以上地震130次，其中，6.0～6.9級(jí)地震23次，7.0～7.9級(jí)地震7次，8級(jí)地震2次（1556年2月2日華縣級(jí)地震和1303年9月25日洪洞8級(jí)地震）。強(qiáng)震主要發(fā)生于山西斷陷帶的忻定盆地、臨汾盆地以及渭河斷陷帶的東部地區(qū)（圖2）。第5代地震動(dòng)參數(shù)區(qū)劃圖在收集汾渭地震帶地震資料、分析地震活動(dòng)性特征的基礎(chǔ)上，充分考慮了地震資料的不完備性和認(rèn)識(shí)的不確定性，并結(jié)合歷史地震活動(dòng)水平與特征（高孟潭，2015；汪素云，1999；陜西省地震局，2005），確定了汾渭地震帶地震活動(dòng)性參數(shù)為b＝0.78，v0＝2.5，m0＝4，muz＝8.5（表 1、圖 3）。

表1 汾渭地震統(tǒng)計(jì)區(qū)不同時(shí)段的地震年平均發(fā)生率

圖3 汾渭地震統(tǒng)計(jì)區(qū)b值及年發(fā)生率v0

按照前述地震目錄模擬方法，模擬汾渭地震帶未來(lái)30、50、70、100年等不同時(shí)間尺度的地震目錄。每一時(shí)間尺度分別模擬2萬(wàn)組。圖4給出了隨機(jī)抽取其中6組時(shí)間尺度為50年的模擬地震目錄的m-t關(guān)系。

統(tǒng)計(jì)2萬(wàn)組50年模擬地震目錄的4級(jí)以上地震數(shù)目，其均值約為125，與理論值125（v0×T））非常接近。以10為間隔，統(tǒng)計(jì)其頻率分布，結(jié)果如圖5所示。由圖5可見，峰值對(duì)應(yīng)的地震數(shù)目區(qū)間在120～130，最大頻率為0.17。

3 模擬參數(shù)的檢驗(yàn)及誤差估計(jì)

數(shù)學(xué)期望、方差、標(biāo)準(zhǔn)差是隨機(jī)變量的重要特性，反映了隨機(jī)變量取值的平均水平、穩(wěn)定程度、集中與離散的程度。對(duì)序列xi，則有

圖4 汾渭地震帶6組50年模擬地震目錄震級(jí)-時(shí)間關(guān)系圖

圖5 汾渭地震帶50年2萬(wàn)組模擬地震目錄地震數(shù)目頻率分布

對(duì)真值為的序列xi，標(biāo)準(zhǔn)差、殘差平方和以及決定性系數(shù)（判定擬合優(yōu)度的參數(shù)）一般用來(lái)判別擬合效果

其中，yk為實(shí)測(cè)值，y＾為預(yù)測(cè)值。

3.1 G-R關(guān)系的年發(fā)生率檢驗(yàn)

為檢驗(yàn)生成的地震目錄是否符合設(shè)定的地震活動(dòng)性參數(shù)，以及是否滿足地震活動(dòng)泊松分布假設(shè)，分別統(tǒng)計(jì)每組目錄每年4級(jí)以上地震發(fā)生次數(shù)，得到50個(gè)樣本（T＝50年），其中，發(fā)生地震i次（取值 0，1，2，…）時(shí)，則地震年發(fā)生次數(shù)的概率密度ρ（i）為

圖6給出了由6組隨機(jī)選擇的50年模擬地震目錄用非線性擬合方法得到的4級(jí)以上年地震發(fā)生次數(shù)泊松分布擬合，表2為每組擬合得到的地震年發(fā)生率v0及其誤差參數(shù)。模擬地震年發(fā)生概率密度與擬合泊松概率密度值的擬合優(yōu)度（決定性系數(shù)）均在大于0.8，殘差平方和均在10-2數(shù)量級(jí)，由此可見，隨機(jī)實(shí)驗(yàn)所生成的目錄能較好符合泊松分布。

圖6 汾渭地震帶6組50年模擬地震目錄地震年發(fā)生次數(shù)泊松分布擬合

按上述方法計(jì)算2萬(wàn)組地震目錄中每組目錄的v0，其均值為2.46，標(biāo)準(zhǔn)差為0.22，以0.2為間隔，統(tǒng)計(jì)其分布（圖 7）。由圖 7（a）可見，峰值對(duì)應(yīng)的v0為 2.4～2.6，與設(shè)定值（2.5）基本接近。圖7（b）給出了 2萬(wàn)組模擬地震目錄，不同v0和不同地震數(shù)目對(duì)應(yīng)的地震目錄數(shù)量分布，當(dāng)?shù)卣饠?shù)目為 120～130、v0為2.4～2.6時(shí)，地震目錄組數(shù)最多，共2120組。

表2 汾渭地震帶6組50年模擬地震目錄泊松分布檢驗(yàn)與平均年發(fā)生率估計(jì)

圖7 汾渭地震統(tǒng)計(jì)區(qū)樣本地震數(shù)目與v0分布關(guān)系（T＝50年）

為了分析不同目錄時(shí)長(zhǎng)對(duì)于結(jié)果的影響，對(duì)30～150年不同時(shí)長(zhǎng)各2萬(wàn)組的地震目錄計(jì)算v0，并計(jì)算不同時(shí)長(zhǎng)的均值、中位值、4分位區(qū)間、極值以及標(biāo)準(zhǔn)誤差（圖8、表3）。由圖8、表3可見，目錄時(shí)間越長(zhǎng)，結(jié)果越收斂，均值越接近理論值，標(biāo)準(zhǔn)誤差越小。

圖8 不同時(shí)長(zhǎng)地震目錄v0統(tǒng)計(jì)

3.2 b值檢驗(yàn)

表4給出了模擬生成的2萬(wàn)組50年地震目錄中，隨機(jī)選擇4組地震目錄采用最小二乘法擬合計(jì)算得到的b值及參數(shù)。由表4可見，擬合優(yōu)度均大于0.9，殘差平方和均小于0.4。

按相同方法計(jì)算2萬(wàn)組地震目錄的b值，其均值為0.76，標(biāo)準(zhǔn)誤差為 0.13，以0.2為間隔，統(tǒng)計(jì)其分布（圖 9）。由圖 9（a）可見，峰值對(duì)應(yīng)的b值為0.76，與設(shè)定值（0.78）接近。圖 9（b）給出了 2萬(wàn)組模擬地震目錄在不同b值、地震數(shù)目時(shí)對(duì)應(yīng)的地震目錄數(shù)量分布。由圖9（b）可見，當(dāng)?shù)卣饠?shù)目為 120～130、b值為0.8左右時(shí)，地震目錄組數(shù)最多，約540組。

表3 模擬目錄不同時(shí)長(zhǎng)v0檢驗(yàn)參數(shù)

對(duì)30～150年不同時(shí)長(zhǎng)各2萬(wàn)組的地震目錄計(jì)算b值，求不同年限的均值、中位值、4分位區(qū)間、極值以及標(biāo)準(zhǔn)誤差（圖 10、表5）。由圖10、表5可見，目錄時(shí)間越長(zhǎng)，b值越收斂，標(biāo)準(zhǔn)誤差和方差越小，越接近理論值。

表4 汾渭地震帶4組50年模擬地震目錄b值擬合檢驗(yàn)參數(shù)

圖9 汾渭地震統(tǒng)計(jì)區(qū)年樣本地震數(shù)目與b值分布間的關(guān)系（T＝50年）

圖10 模擬目錄不同時(shí)長(zhǎng)b值統(tǒng)計(jì)

模擬結(jié)果表明，無(wú)論是v0還是b值，模擬的地震目錄時(shí)間越長(zhǎng)，標(biāo)準(zhǔn)誤差越小，均值越接近設(shè)定值，越能代表泊松模型描述的地震活動(dòng)規(guī)律。

地震中長(zhǎng)期預(yù)測(cè)是對(duì)某一地區(qū)今后數(shù)年到數(shù)十年地震形勢(shì)和危險(xiǎn)性的估計(jì)與預(yù)測(cè)（王曉青等，1996、2006），其結(jié)果既可直接服務(wù)于社會(huì)，又是短臨預(yù)測(cè)的基礎(chǔ)依據(jù)。本文采用上述模擬生成的2萬(wàn)組100年地震目錄，統(tǒng)計(jì)指定震級(jí)的目錄數(shù)占總目錄的比例，得到汾渭地震帶未來(lái)100年7級(jí)以上、8級(jí)以上地震的發(fā)震概率分別為0.66、0.19，與泊松模型理論計(jì)算的結(jié)果0.67、0.18基本一致。

表5 模擬目錄不同時(shí)長(zhǎng)b值檢驗(yàn)參數(shù)表

4 結(jié)論與討論

蒙特卡洛模擬通過(guò)對(duì)研究對(duì)象進(jìn)行隨機(jī)獨(dú)立抽樣，使其分布滿足已知或設(shè)定概率分布來(lái)實(shí)現(xiàn)模擬。該方法通過(guò)大量重復(fù)實(shí)驗(yàn)，能夠真實(shí)地接近實(shí)際物理過(guò)程，受幾何條件限制小，對(duì)于小概率事件處理較好（Huh et al，2016）。通過(guò)改變地震事件的震級(jí)-頻度關(guān)系形式、時(shí)間分布模型，采用蒙特卡洛模擬可得到各種復(fù)雜模型條件下的模擬地震目錄。模擬目錄在一定程度上填補(bǔ)了由于監(jiān)測(cè)手段不足或歷史記錄不完整（無(wú)時(shí)間）、年代記載空缺所造成的地震目錄缺失，使得那些長(zhǎng)期無(wú)地震活動(dòng)地區(qū)的地震序列更加完整。

基于地震活動(dòng)的泊松分布模型、古登堡-里克特震級(jí)-頻度（G-R）關(guān)系假設(shè)，運(yùn)用反抽樣法和直接抽樣法生成隨機(jī)數(shù)模擬地震目錄，依據(jù)汾渭地震帶的地震活動(dòng)性參數(shù)，模擬了該地震帶30、50、100年等不同時(shí)間長(zhǎng)度的地震目錄，檢驗(yàn)結(jié)果表明，用該方法得到的地震目錄體現(xiàn)了地震活動(dòng)統(tǒng)計(jì)特征。在此基礎(chǔ)上，應(yīng)用模擬地震目錄計(jì)算得到了汾謂地震帶未來(lái)100年發(fā)生7級(jí)以上、8級(jí)以上地震的概率分別為 0.66、0.19。以1815年平陸 6? 級(jí)地震以來(lái)202年無(wú)6.5級(jí)以上地震和1998年張北6.2級(jí)地震以來(lái)19年無(wú)6級(jí)以上地震為條件，未來(lái)10年發(fā)生6.5級(jí)以上、6.2級(jí)以上地震的條件概率分別為0.22、0.31。

由于文中使用的時(shí)間上的泊松模型與真實(shí)地震活動(dòng)常表現(xiàn)的叢集模型不一定吻合，考慮到地震之間的相互激發(fā)問(wèn)題（泊松+叢集），今后可改變抽樣函數(shù)，引入 ETAS模型（Marzocchi et al，2009），以期在古地震記錄較少地區(qū)以及地質(zhì)結(jié)構(gòu)復(fù)雜條件下能提供更精確的結(jié)果（Parsons，2008）。

致謝：感謝中國(guó)地震局地球物理研究所提供的第5代地震動(dòng)參數(shù)區(qū)劃圖資料。感謝陜西省地震局對(duì)地震目錄資料收集給予的支持和幫助，感謝審稿專家的意見和建議。