一種基于蒙特卡羅模擬的發(fā)震概率計算方法*

郭　星　潘　華

1) 中國北京100082環(huán)境保護部核與輻射安全中心2) 中國北京100081中國地震局地球物理研究所

?

一種基于蒙特卡羅模擬的發(fā)震概率計算方法*

郭星1)潘華2),*

1) 中國北京100082環(huán)境保護部核與輻射安全中心2) 中國北京100081中國地震局地球物理研究所

針對大震發(fā)生概率計算過程中的不確定性，本文分別對不確定性及其處理方法進行了探討． 考慮到不確定構成的復雜性，提出了一種基于蒙特卡羅模擬的大震發(fā)生概率計算方法, 并以東昆侖斷裂帶塔藏段為計算實例，利用蒙特卡羅法處理發(fā)震概率計算過程中的各種不確定性. 結果表明, 古地震數(shù)據(jù)的不完整性對計算結果的影響很大. 本文采用邏輯樹法考慮古地震數(shù)據(jù)的不完整性，得到塔藏段未來100年的大震發(fā)生概率為0.12．

大震發(fā)生概率蒙特卡羅不確定性

引言

考慮到大地震的記憶性，Utsu(1972)，Rikitake(1974)和Hagiwara(1974)基于Reid(1910)提出的彈性回跳理論，提出了一種更新模型，該模型假定斷層上大地震的復發(fā)符合更新過程． 國外研究人員先后提出了多種概率分布來描述這種更新模型，具體的概率分布模型包括高斯分布(Rikitake，1974)、韋布爾分布(Hagiwara，1974)、對數(shù)正態(tài)分布(Nishenko，Buland，1987)、布朗過程時間(Brownian passage time，簡寫為BPT)分布(Ellsworthetal，1999； Matthewsetal，2002)等模型．

對于任何一種概率分布模型，若已知大震復發(fā)的概率密度函數(shù)f(T)和最近一次大震的離逝時間Te，則可利用下式計算斷層源在未來一段時間ΔT內(nèi)發(fā)生大震的條件概率(Wesnousky，1986)：

(1)

不確定性一般可分為認知不確定性和隨機不確定性，本文在此基礎上又對不確定性的定量化過程和處理過程分別進行了深入分析． 按照概率分布類型不同，將不確定性分布分為連續(xù)型分布和離散型分布; 按照認知主體的不同，將不確定性的處理分為對個體自身認知不確定性的處理和對不同個體之間認知差異性的處理． 最后進一步探討了具體的不確定性處理方法和原則．

此外，大震發(fā)生概率計算過程中的不確定性構成往往很復雜，例如，不確定性因素較多，而不同不確定性因素之間又存在相關性，這使得利用數(shù)值方法實現(xiàn)起來比較困難． 針對該問題，本文首先提出利用蒙特卡羅模擬法對大震發(fā)生概率模型中的各種不確定性分布直接進行反復隨機抽樣； 然后，通過對大量的隨機模擬結果進行統(tǒng)計分析，得到最終的大震發(fā)生概率計算結果； 最后，以東昆侖斷裂帶塔藏段的大震發(fā)生概率為計算實例，在充分考慮各種不確定性的基礎上，建立了塔藏段的大震發(fā)生概率模型，并利用蒙特卡羅法計算塔藏段未來100年的大震發(fā)生概率．

1　大震發(fā)生概率模型的不確定性分析

1.1不確定性分類

在各種地震危險性模型中，對于不確定性的劃分，一般可分為隨機不確定性和認知不確定性(Bommer，2003)． 隨機不確定性反映的是自然本身所固有的各種隨機性； 認知不確定性反映的是基于不完備的信息和知識，人的主觀認知和判斷的不準確性．

大震危險性分析中的不確定性分析是對不確定性進行定量化處理，進而在大震危險性分析中考慮各種不確定性因素． 隨機不確定性和認知不確定性有時可以分離，有時則很難分離． 因此，為了探討不確定性的具體處理方法，本文在不確定性處理過程中，分別按照概率分布類型的不同和認知主體的不同進行劃分，而未將隨機不確定性和認知不確定性分離開進行討論．

首先，按照概率分布類型不同，本文將不確定性分布分為連續(xù)型分布和離散型分布． 對于不確定性的定量化描述，根據(jù)具體情況的不同，可以采用連續(xù)型分布，也可以采用離散型分布． 其中，連續(xù)型分布更為客觀，但建立相對可靠的分布模型需要足夠的樣本支持，同時連續(xù)型的概率分布不方便計算； 離散型分布的確定相對主觀，但卻便于計算． 其次，按照認知主體的不同，本文將不確定性的處理分為對個體自身認知不確定性的處理和對不同個體之間認知差異性的處理． 在不確定性的實際處理中，不僅要考慮研究個體，同時還需考慮不同研究人員之間所得出的不確定性定量化結果的差異性(SeniorSeismicHazardAnalysisCommittee，1997)． 這種差異性反映的是人與人之間存在的認知分歧，為了避免計算結果過于主觀和片面，需要在大震危險性分析中考慮這種分歧．

圖1　不確定性的處理過程Fig.1　The process of uncertainty

1.2不確定性的處理方法

為便于討論不確定性的具體處理方法，根據(jù)本文所采用的分類方式，并參考前人在概率地震危險性分析中不確定性處理方面的研究成果(Senior Seismic Hazard Analysis Committee，1997)，本文對不確定性的處理過程參見圖1.

1.2.1采用連續(xù)型分布來定量化不確定性

連續(xù)型分布較離散型分布更為精確，但其計算量較大，因此如果計算條件允許，可以采用連續(xù)型分布來定量化不確定性． 連續(xù)型分布的具體分布形態(tài)又分為均勻分布和非均勻分布．

1) 均勻分布. 客觀信息量越少，不確定性越大，這是個不變的原則，必須體現(xiàn)在大震危險性分析中的每個具體環(huán)節(jié)中． 對于基本不可知因素或完全不可知因素的不確定性評價，一般采用均勻加權法，其本質(zhì)上屬于拉普拉斯決策法，即等可能決策法． 均勻分布在實際地震危險性分析中應用廣泛，例如： 在潛在震源區(qū)內(nèi)任何地方發(fā)生地震的可能性是相同的(Cornell，1968)； 古地震發(fā)生在其年代測定不確定性的上下限內(nèi)任意一年的可能性也是相同的．

2) 非均勻分布. 若要用非均勻分布來表示不確定性，首先必須假設一個非均勻分布模型，然后通過參數(shù)估計法確定該分布模型的參數(shù)． 具體的參數(shù)估計法包括矩估計、最大似然估計、貝葉斯估計等． 其中矩估計和最大似然估計多用于計算某些自然隨機分布的參數(shù)，而貝葉斯估計則多用于計算某些模型參數(shù)的認知不確定性分布． 如果有連續(xù)性的證據(jù)或數(shù)據(jù)，則可以采用非均勻分布來表示不確定性，但必須先定義一個分布模型，甚至是出于個體的主觀認識．

1.2.2采用離散型分布來定量化不確定性

邏輯樹是一種常用的離散化不確定性處理方法． 該方法一般分為兩種： 一種是加權邏輯樹法，一種是不加權邏輯樹法．

1) 不加權邏輯樹． 對于個體的認知不確定性處理，類似于連續(xù)型分布中的均勻分布，不加權邏輯樹也是一種等可能加權． 對于不同個體之間的認知不確定性的差異性處理，可以采用不加權邏輯樹，這意味著在決策過程中不同個體的認知被賦予相同的權重，即不同個體之間是平等的．

2) 加權邏輯樹． 利用加權邏輯樹對個體認知不確定性的處理，是單獨地研究個體基于有限的地震地質(zhì)資料和信息，對大震危險性評估過程中某個環(huán)節(jié)的認知不確定性分布所作的離散處理，其中每個邏輯樹分支均被主觀地賦予一定的權重． 根據(jù)有限的地震地質(zhì)資料和信息，對兩種或兩種以上可能(或離散型的參數(shù)值)的權重進行主觀評定． 最終評定結果的確定須參考信息的完備程度，信息量越少，未知因素越多，不確定性也就越大． 例如，某個不確定性環(huán)節(jié)包括A和B兩種可能，即使所有的已知信息均支持A，但如果信息量非常有限，也不可以給出確定的結論，而應參考信息量的完備程度，在充分考慮不確定性的基礎上進行權重評估．

對于個體之間的認知不確定性的差異性處理，也可以采用加權邏輯樹法，其綜合了不同個體認知的權重不平等性和單獨個體的認知不確定性． 此外，在不確定性處理過程中還需注意以下兩點：

1) 必須保證權重賦予和不確定性分布范圍確定過程的一致性原則． 一個獨立的人或者獨立的小組，在大震危險性評估過程中每個具體的不確定性環(huán)節(jié)，必須遵循相同的不確定分級原則或權重分配原則．

2) 在評估不確定性分布范圍的過程中，需要保守考慮，但又不可以無限地保守． 一些概率極低的極端事件對計算結果的貢獻很小，也很難進行概率評估和計算，所以，在實際大震危險性評估過程中可以對不確定性分布作適當?shù)慕財嗵幚恚?/p>

2　計算實例

本文以東昆侖斷裂帶東段塔藏段的大震發(fā)生概率為計算實例，具體研究在大震發(fā)生概率計算過程中如何考慮不確定性．

2.1東昆侖斷裂帶塔藏段的古地震

東昆侖斷裂帶東段位于巴顏喀拉地塊與西秦嶺地塊邊界斷裂的東段，全長約330 km，自西向東可劃分為3段： 瑪沁段、瑪曲段和塔藏段(青海省地震局，中國地震局地殼應力研究所，1999)． 國內(nèi)外研究人員在東昆侖斷裂帶東段上獲得了許多古地震數(shù)據(jù)，為評估該段的大地震發(fā)生概率提供了重要依據(jù)． 其中，瑪沁段和瑪曲段上的古地震數(shù)據(jù)比較豐富，而塔藏段上的古地震數(shù)據(jù)量則比較少，完整性較差，最近一次大震的離逝時間也存在很大的不確定性．

根據(jù)李正芳等(2012)對前人古地震研究結果的分析，本文給出了塔藏段古地震的發(fā)生時間分別為距今(4693±151)年，(7304±500)年和(9136±131)年，其中存在4700年左右的大地震空白期，很可能存在歷史地震漏記的現(xiàn)象． 另根據(jù)李正芳等(2012)對四川松潘(位于塔藏段東南方向約50 km)歷史文獻記載的分析，推斷塔藏段最近一次地震的發(fā)生時間應在公元638年之前，即最近一次地震的離逝時間應介于1377—(4693±151)年之間．

2.2不確定性分析

由于塔藏段古地震數(shù)據(jù)量比較少，完整性較差，不管是大震復發(fā)間隔還是大震離逝時間，均存在著非常大的不確定性． 具體的不確定性因素包括以下3個方面：

1) 在離逝時間Te的確定過程中，需要考慮古地震記錄和歷史地震記載的不完整性所造成的最近一次大地震發(fā)生年代的不確定性.

首先，古地震的完整性研究一直是古地震研究的難題，不同斷裂上的古地震研究深度不同，地層的完整性情況也不同，造成不同斷裂和不同時段的古地震記錄的完整性情況相差較大(冉勇康，鄧起東，1999)； 其次，我國歷史地震資料保存得相對較多，但就各個歷史時期和不同區(qū)域來說，歷史地震記載的完整性程度相差較大(陳春梅，任雪梅，2014)． 因此，無法通過客觀的統(tǒng)計分析方法得到最近一次大地震發(fā)生時間的不確定性分布，也無法直接判定某時間段內(nèi)的大震記錄是否完整或缺失的數(shù)量，只能主觀評估已知信息量能否足夠(或充分)判定給定時間段的大地震記載是完整的．

本文提出采用邏輯樹法來分析由古地震記錄和歷史地震記載的不完整性所造成離逝時間Te的認知不確定性． 邏輯樹在此處共有兩個分支： ① 根據(jù)已知的區(qū)域古地震研究成果和區(qū)域歷史地震記載情況，無法確定給定時間段的大地震記載是完整的，即信息量有限，無法得出確定性的結論； ② 根據(jù)已知的區(qū)域古地震研究成果和區(qū)域歷史地震記載情況，可以確定給定時間段的大地震記載是完整的，即信息量是足夠充分的，可以得出確定性的結論． 根據(jù)已知的數(shù)據(jù)和資料，并征求相關方面專家的意見，對兩個邏輯樹分支的權重進行主觀賦值．

2.3考慮不確定性的塔藏段大震發(fā)生概率計算模型

在Reid(1910)彈性回跳理論的基礎上，國內(nèi)外研究人員先后提出了多種符合更新過程的概率模型，其中比較常用的概率分布模型有對數(shù)正態(tài)分布模型(Nishenko，Buland，1987)、BPT(Ellsworthetal，1999； Matthewsetal，2002)模型等．

本文未考慮模型本身的不確定性，直接采用BPT模型計算大震發(fā)生概率. BPT分布也被稱作逆高斯分布(Seshadri，1983)，其概率密度函數(shù)為

(2)

參考對不確定性處理方法的分類討論，本文采用加權邏輯樹法來考慮古地震的不完整性． 考慮到每個邏輯樹分支所對應的計算過程均存在一定差異，本文分成A 和B兩個相互獨立的蒙特卡羅模擬過程分別進行計算． 其中分支A假定根據(jù)有限的信息，無法確定給定時間段的大地震記載是完整的； 分支B假定信息量足夠充分，可以確定給定時間段的大地震記載是完整的． 根據(jù)已知的數(shù)據(jù)和資料，并征求相關方面專家的意見，對兩個邏輯樹分支的權重進行主觀賦值．

圖2　塔藏段大震發(fā)生概率計算的分支A流程圖Fig.2　The process chart of branch A for calculating the occurrence probability of large earthquakes on the Tazang fault

分支A的具體實現(xiàn)步驟如圖2所示：

(3)

而復發(fā)間隔為T22的概率密度函數(shù)為

(4)

(5)

(6)

從保守和計算效率的角度考慮，判定過程中并未考慮記載空段內(nèi)漏記兩次或兩次以上古地震的小概率特殊事件的發(fā)生．

圖3　塔藏段大震發(fā)生概率計算的分支B流程圖Fig.3　The process chart of branch B for calculating the occurrence probability of large earthquakes on the Tazang fault

6) 反復模擬10萬次，計算出10萬次模擬結果的平均值P．

分支B的具體實現(xiàn)步驟與分支A類似(圖3)，只是第4步有所差異，考慮到分支B假定不可能出現(xiàn)古地震漏記，則直接取大震離逝時間Te=Te1．

2.4大震發(fā)生概率計算結果

分支A的大震發(fā)生概率計算結果為0.11； 分支B的大震發(fā)生概率計算結果為0.16． 而對于A和B兩個邏輯分支的權重賦值則需要綜合考慮以下幾個方面： ① 研究區(qū)域的人類活動歷史和歷史地震記載的完整性情況； ② 研究斷裂上的地層完整性情況； ③ 研究斷裂上的古地震研究深度．

基于上述3個方面的信息，并征求相關方面專家的意見，本文給出的A和B兩個邏輯分支的權重賦值WA和WB分別為0.8和0.2，最后計算出的塔藏段未來100年的大震發(fā)生概率為

P=PAWA+PBWB=0.11×0.8+0.16×0.2=0.12.

(8)

3　討論與結論

本文對大震發(fā)生概率評估過程中的各種不確定性進行了系統(tǒng)分析，在認知不確定性和隨機不確定性的基礎上對不確定性的定量化處理過程進行了分類研究，并對不確定性的具體定量化方法進行了梳理．

考慮到大震復發(fā)概率計算過程中存在大量的不確定性，且不確定性的構成又較為復雜，數(shù)值方法實現(xiàn)起來比較困難，故此類不確定性往往被忽略掉． 本文首先明確了大震復發(fā)概率模型中的各種不確定性，特別是常被忽略的不確定性，然后利用蒙特卡羅模擬法對大震發(fā)生概率模型中的各種定量化的不確定性直接進行隨機抽樣，最后通過對大量隨機模擬結果的統(tǒng)計分析得到了最終的大震發(fā)生概率計算結果． 雖然計算過程相對復雜，但可以盡可能地減少一些主觀因素，使得計算結果更為客觀．

本研究的重點在于探討如何利用蒙特卡羅法考慮發(fā)震概率計算過程中的各種不確定性，對于各個不確定性分支的具體賦值過程則未進行更深入的研究． 而不確定性分支的權重將直接影響最終的計算結果，需要作大量的地震地質(zhì)調(diào)查和統(tǒng)計分析工作，也需要采取科學有效的方法對多位專家的意見進行綜合． 在實際大震危險性評估中，怎樣考慮和處理這種人與人之間的認知分歧，也需進一步研究．

陳春梅，任雪梅. 2014. 我國大陸5個地區(qū)歷史地震資料記載的完整性分析與比較[J]. 防災減災學報，30(2): 66--70.

Chen C M，Ren X M. 2014. Integrity analysis and comparison of five regional historical earthquake records in mainland China[J].JournalofDisasterPreventionandReduction，30(2): 66--70 (in Chinese).

郭星，潘華. 2015. 強震復發(fā)間隔變異系數(shù)的一種計算方法[J]. 地震學報，37(3): 411--419.

Guo X，Pan H. 2015. A method for computing the aperiodicity parameter of the strong earthquake recurrence interval[J].ActaSeismologicaSinica，37(3): 411--419 (in Chinese).

郭星，潘華. 2016. 強震復發(fā)概率模型中的參數(shù)不確定性研究[J]. 地震學報，38(2): 298--306.

Guo X，Pan H. 2016. Parameter uncertainty analysis on probability model of strong earthquake recurrence[J].ActaSeismologicaSinica，38(2): 298--306 (in Chinese).

李正芳，周本剛，冉洪流. 2012. 運用古地震數(shù)據(jù)評價東昆侖斷裂帶東段未來百年的強震危險性[J]. 地球物理學報，55(9): 3051--3065.

Li Z F，Zhou B G，Ran H L. 2012. Strong earthquake risk assessment of eastern segment on the East Kunlun fault in the next 100 years based on paleo-earthquake data[J].ChineseJournalofGeophysics，55(9): 3051--3065 (in Chinese).

毛鳳英，張培震. 1995. 古地震研究中的逐次限定法與新疆北部主要斷裂帶的古地震研究[G]∥活動斷裂研究(4). 北京: 地震出版社： 153--164.

Mao F Y，Zhang P Z. 1995. Progressive constraining method in paleo-earthquake study and the research on paleo-earthquakes along the major active faults in northern Xinjiang[G]∥ResearchonActiveFault(4). Beijing: Seismological Press: 153--164 (in Chinese).

青海省地震局，中國地震局地殼應力研究所. 1999. 東昆侖活動斷裂帶[M]. 北京: 地震出版社： 1--186.

Earthquake Administration of Qinghai Province，Institute of Crustal Dynamics，China Earthquake Administration. 1999.EastKunlunActiveFaultZone[M]. Beijing: Seismological Press： 1--186 (in Chinese).

冉勇康，鄧起東. 1999. 古地震學研究的歷史、現(xiàn)狀和發(fā)展趨勢[J]. 科學通報，44(1)： 12--20.

Ran Y K，Deng Q D. 1999. History，status and trend about the research of paleoseismlogy[J].ChineseScienceBulletin，44(10)： 880--889.

Bommer J J. 2003. Uncertainty about the uncertainty in seismic hazard analysis[J].EngGeol，70(1/2): 165--168.

Cornell C A. 1968. Engineering seismic risk analysis[J].BullSeismolSocAm，58(5): 1583--1606.

Ellsworth W L，Matthews M V，Nadeau R M，Nishenko S P，Reasenberg P A，Simpson R W. 1999.APhysically-BasedEarthquakeRecurrenceModelforEstimationofLong-TermEarthquakeProbabilities[R]. Reston，Virginia: US Geological Survey Open-File Report： 99--522.

Hagiwara Y. 1974. Probability of earthquake occurrence as obtained from a Weibull distribution analysis of crustal strain[J].Tectonophysics，23(3): 313--318.

Matthews M V，Ellsworth W L，Reasenberg P A. 2002. A Brownian model for recurrent earthquakes[J].BullSeismolSocAm，92(6): 2233--2250.

Nishenko S P，Buland R A. 1987. A generic recurrence interval distribution for earthquake forecasting[J].BullSeismolSocAm，77(4): 1382--1399.

Reid H F. 1910.TheMechanicsoftheEarthquake，TheCaliforniaEarthquakeofApril18，1906[R]. Washington，DC: State Earthquake Investigation Commission: 43--47.

Rikitake T. 1974. Probability of earthquake occurrence as estimated from crustal strain[J].Tectonophysics，23(3): 299--312.

Senior Seismic Hazard Analysis Committee. 1997.RecommendationsforProbabilisticSeismicHazardAnalysis:GuidanceonUncertaintyandUseofExperts[R]. Washington，DC: US Nuclear Regulartory Commission: 1--256.

Seshadri V. 1983. The inverse Gaussian distribution: Some properties and characterizations[J].CanadianJStatist，11(2): 131--136.

Utsu T. 1972.LargeEarthquakesNearHokkaidoandtheExpectancyoftheOccurrenceofALargeEarthquakeOffNemuro[R]. Hokkaido: Coordinating Committee for Earthquake Prediction: 7--13.

Wesnousky S G. 1986. Earthquakes，Quaternary faults，and seismic hazard in California[J].JGeophysRes，91(B12): 12587--12631.

A method for calculating occurrence probability of large earthquakes based on Monte Carlo simulation

Guo Xing1)Pan Hua2),*

1)NuclearandRadiationSafetyCenter，MinistryofEnvironmentProtection，Beijing100082，China2)InstituteofGeophysics，ChinaEarthquakeAdministration，Beijing100081，China

According to the uncertainty in the process of calculating the occurrence probability of large earthquakes，a study is made on uncertainty and its dealing methods. Considering the complexity of uncertainty，this paper presents a method for calculating the occurrence probability of large earthquakes based on Monte Carlo simulation. With the Tazang segment of eastern Kunlun fault zone as an example，we deal with different kinds of uncertainties in calculating the occurrence probability of large earthquakes using Monte Carlo method. The result shows that the incompleteness of paleo-earthquakes data has great effect on the calculation result. With the logical tree to deal with the incompleteness of paleo-earthquakes data，the occurrence probability of large earthquakes is 0.12 in the next 100 years on the Tazang segment.

large earthquake； occurrence probability； Monte Carlo； uncertainty

國家科技支撐項目(2012BAK15B01-08)資助.

2016-03-01收到初稿，2016-07-18決定采用修改稿.

e-mail: panhua.mail@163.com

10.11939/jass.2016.05.012

P315.5

A

郭星，潘華. 2016. 一種基于蒙特卡羅模擬的發(fā)震概率計算方法. 地震學報, 38(5): 785--793. doi:10.11939/jass.2016.05.012.

Guo X, Pan H. 2016. A method for calculating occurrence probability of large earthquakes based on Monte Carlo simulation.ActaSeismologicaSinica, 38(5): 785--793. doi:10.11939/jass.2016.05.012.