邵霄怡 王曉青 竇愛霞 袁小祥

中國地震局地震預(yù)測研究所，北京市復(fù)興路63號 100036

0 引言

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

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

基于力學(xué)與統(tǒng)計模型相結(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)地震活動模型（金欣等，2017；周仕勇，2008）。以這些力學(xué)模型所反映的應(yīng)力-應(yīng)變過程為基礎(chǔ)，結(jié)合細(xì)胞自動機(jī)模型（劉桂萍等，1995、2000）、耦合應(yīng)力釋放統(tǒng)計模型（AIC）（李紅等，2015）、布朗時間過程模型（BPT模型）、時間可預(yù)報模型或基于經(jīng)驗(yàn)的模型（Working Group on California Earthquake Probabilities，2003；Hainzl et al，2007）、EEPAS模型（Rhoades et al，2011；畢金孟等，2017）等，可實(shí)現(xiàn)對地震目錄的模擬。彈簧-滑塊模型或斷層系統(tǒng)模型具有比較明確的物理含義和較為嚴(yán)格的理論基礎(chǔ)，但因影響地震發(fā)生的因素較多，模型參數(shù)的設(shè)置與實(shí)際斷裂構(gòu)造間存在一定差異，模型越復(fù)雜，假設(shè)條件就越多，這造成模擬的地震目錄不確定性較大，模擬結(jié)果的檢驗(yàn)存在一定的困難（李紅等，2015）。

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

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

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

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

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

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

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

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

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

假定地震的震級分布滿足截斷的G-R關(guān)系（Cornell，1968），則地震震級分布函數(shù)F（m）（Youngs et al，1985）為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖8 不同時長地震目錄v0統(tǒng)計

3.2 b值檢驗(yàn)

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

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

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

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

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

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

圖10 模擬目錄不同時長b值統(tǒng)計

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

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

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

4 結(jié)論與討論

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

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

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

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