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

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

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

0 引言

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

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

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

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

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

1 基于蒙特卡洛隨機獨立重復實驗的地震目錄模擬方法

蒙特卡洛方法（Monte Carlo method），也稱統計模擬方法或隨機抽樣技術，是以概率和統計理論為基礎的數值計算方法。該方法具有對小概率事件處理較好，受幾何條件限制小、收斂速度與問題的維數無關等特性，可用于地震、地震災害分析預測等領域（劉善琪等，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）和地震風險評估（Weatherill，2009；Yazdani et al，2012；Bourne et al，2015；劉甲美等，2016）等方面得到應用。

1.1 蒙特卡洛隨機獨立重復實驗方法

對于一個給定的物理（或數學）系統，如果系統中某一核心變量s與其概率分布函數F（s）的關系已知，即

函數F（s）為值域［0，1］的單調增函數。通過隨機獨立抽樣確定［0，1］區間均勻隨機數r，對式（1）進行反變換，得

即可得到核心變量s的取值。通過多次獨立重復抽樣，可以確定一系列核心變量值，完成該系統參數的隨機模擬實驗。上述方法稱為反抽樣法，是蒙特卡洛隨機實驗的基本方法之一（Huh et al，2016）。

1.2 指定年限內震級分布確定

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

其中，β＝2.3b，b為研究區G-R關系的系數；m0為震級下限，muz為震級上限。F（m）為［0，1］區間的單調增函數，令F（m）＝r，對式（3）進行變換，得

通過產生［0，1］區間的隨機數r，確定F（m）的值，再根據式（4）求得地震震級m。

1.3 指定年限內地震數量和相鄰兩次地震時間間隔確定

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

對式（5）進行求和，得到地震發生次數的概率分布函數為

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

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

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

因此，通過蒙特卡洛隨機獨立重復實驗，依據年發生率v0和b值，可確定區域地震震級分布以及地震事件數量k。在此基礎上，針對每次地震事件，確定震級m和發震時間t，得到滿足泊松分布假設的模擬地震目錄（圖1）。

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

圖1 地震目錄生成的技術流程

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

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

表1 汾渭地震統計區不同時段的地震年平均發生率

圖3 汾渭地震統計區b值及年發生率v0

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

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

3 模擬參數的檢驗及誤差估計

數學期望、方差、標準差是隨機變量的重要特性，反映了隨機變量取值的平均水平、穩定程度、集中與離散的程度。對序列xi，則有

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

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

對真值為的序列xi，標準差、殘差平方和以及決定性系數（判定擬合優度的參數）一般用來判別擬合效果

其中，yk為實測值，y＾為預測值。

3.1 G-R關系的年發生率檢驗

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

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

圖6 汾渭地震帶6組50年模擬地震目錄地震年發生次數泊松分布擬合

按上述方法計算2萬組地震目錄中每組目錄的v0，其均值為2.46，標準差為0.22，以0.2為間隔，統計其分布（圖 7）。由圖 7（a）可見，峰值對應的v0為 2.4～2.6，與設定值（2.5）基本接近。圖7（b）給出了 2萬組模擬地震目錄，不同v0和不同地震數目對應的地震目錄數量分布，當地震數目為 120～130、v0為2.4～2.6時，地震目錄組數最多，共2120組。

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

圖7 汾渭地震統計區樣本地震數目與v0分布關系（T＝50年）

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

圖8 不同時長地震目錄v0統計

3.2 b值檢驗

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

按相同方法計算2萬組地震目錄的b值，其均值為0.76，標準誤差為 0.13，以0.2為間隔，統計其分布（圖 9）。由圖 9（a）可見，峰值對應的b值為0.76，與設定值（0.78）接近。圖 9（b）給出了 2萬組模擬地震目錄在不同b值、地震數目時對應的地震目錄數量分布。由圖9（b）可見，當地震數目為 120～130、b值為0.8左右時，地震目錄組數最多，約540組。

表3 模擬目錄不同時長v0檢驗參數

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

表4 汾渭地震帶4組50年模擬地震目錄b值擬合檢驗參數

圖9 汾渭地震統計區年樣本地震數目與b值分布間的關系（T＝50年）

圖10 模擬目錄不同時長b值統計

模擬結果表明，無論是v0還是b值，模擬的地震目錄時間越長，標準誤差越小，均值越接近設定值，越能代表泊松模型描述的地震活動規律。

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

表5 模擬目錄不同時長b值檢驗參數表

4 結論與討論

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

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

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

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