龍門山地區(qū)的地震活動性廣義帕累托模型構(gòu)建

任夢依

摘要：利用1930-2016年龍門山地區(qū)M≥4.5歷史地震目錄，通過時空窗法減少余震對統(tǒng)計結(jié)果的影響，構(gòu)建了該地區(qū)廣義帕累托分布的超出量分布模型，估計了龍門山地區(qū)震級上限。結(jié)果表明：廣義帕累托分布較好地擬合了龍門山地區(qū)強震數(shù)據(jù)，形狀參數(shù)估計均為負(fù)值，最大震級分布存在有限上界，震級上限為8.3。

關(guān)鍵詞：地震活動性；廣義帕累托分布；震級上限；龍門山地區(qū)

中圖分類號：P315.5 文獻(xiàn)標(biāo)識碼：A 文章編號： 1000-0666（2018）02-0226-07

0 引言

地震活動性模型的構(gòu)建方法研究可為地震區(qū)劃和地震危險性分析提供重要的技術(shù)支撐。對于地震活動性模型的改進(jìn)方法研究，一直是地震學(xué)領(lǐng)域的研究熱點。概率方法是目前構(gòu)建地震活動性模型中最常使用的分析方法（胥廣銀，2003），如基于截斷G-R關(guān)系模型的方法（Cornell，1968）以及部分學(xué)者在此基礎(chǔ)上提出的改進(jìn)方法（徐偉進(jìn)，高孟潭，2012），美國現(xiàn)行的地震區(qū)劃圖以及第五代《中國地震動參數(shù)區(qū)劃圖》中均主要應(yīng)用該方法構(gòu)建地震活動性模型（潘華等，2009；潘華，李金臣，2016）。然而，使用基于截斷G-R關(guān)系的方法，需要先驗地選取起算震級和設(shè)定震級上限。大量的實際地震資料分析也表明，用G-R關(guān)系模型檢驗地震震級分布情況時，在許多情況下震級高端或低端會出現(xiàn)明顯的“掉頭”現(xiàn)象（胡聿賢，1999）。由于實際資料對G-R關(guān)系的偏離，特別是強震級段的偏離，一定程度上影響了利用該關(guān)系式推測未來強震發(fā)生危險性的可信度，因此，如何對上述關(guān)系進(jìn)行改善，特別是針對強震數(shù)據(jù)的擬合一直是地震學(xué)家研究關(guān)注的焦點。

Balkema和Haan（1974），Pickands（1975）研究發(fā)現(xiàn)，對于充分大的閾值，隨機變量超過閾值超出量的極值分布為廣義帕累托分布，這為研究隨機變量尾部分布特征提供了理論依據(jù)。近年來，由極值模型發(fā)展而來的廣義帕累托分布模型逐漸被應(yīng)用到潛源區(qū)地震活動性模型的構(gòu)建中，如Pisarenko和Sornette（2003）利用廣義帕累托分布分析了哈佛地震目錄中18個地震區(qū)的淺層地震地震矩分布；Huyse等（2010）利用太平洋地震工程研究中心的地面峰值加速度數(shù)據(jù)和殘差數(shù)據(jù)比較了對數(shù)正態(tài)分布與廣義帕累托分布的擬合優(yōu)度；錢小仕等（2013）利用廣義帕累托分布分析了各活動地塊邊界帶的強震震級分布特征；田建偉等（2017）利用廣義帕累托分布估計了馬尼拉海溝俯沖帶潛在地震海嘯源區(qū)的地震危險性。該方法利用一段時限內(nèi)潛源區(qū)的大震震級數(shù)據(jù)，不僅無需先驗地選定震級下限和震級上限，而且可充分利用大震級區(qū)段的信息，適合于歷史地震記錄時間長但低震級地震記錄缺失的地區(qū)，便于構(gòu)建潛源區(qū)強震活動模型。

在構(gòu)建廣義帕累托模型前，需充分了解、分析潛源區(qū)地震地質(zhì)信息，確定潛源區(qū)邊界范圍，對于地震樣本的選取也是十分重要的。另外，在利用廣義帕累托分布構(gòu)建模型時，需提前排除余震的影響，盡量滿足事件相互獨立的條件。本文選取地震活動較為頻繁的龍門山地區(qū)為研究區(qū)，在分析該地區(qū)地震地質(zhì)背景以及地震活動特征的基礎(chǔ)上，利用1930-2016年歷史地震數(shù)據(jù)的大震級區(qū)段，確定震級閾值，構(gòu)建了龍門山地區(qū)廣義帕累托分布的超出量分布模型，用于估計該區(qū)域的震級上限。

1 廣義帕累托模型構(gòu)建方法

1.1 廣義帕累托分布

設(shè)X1，X2，…，Xn為地震震級隨機變量列，假設(shè)它們相互獨立且服從同一分布，記其最大震級為Mn=max（X1，X2，…，Xn），若存在{an>0，bn∈R}和非退化分布·函數(shù)G（x），使得：則G稱為極值分布。Fisher和Tippett（1928）證明了極值分布的3大類型定理，指出G必屬于Gum-bel分布、Frechet分布、Weibull分布之一。廣義極值理論將這3種分布類型統(tǒng)一為1個分布類型，由位置參數(shù)μ、尺度參數(shù)σ、形狀參數(shù)ξ來表示（Coles，2001；史道濟，2006），如式（2）所示：選定震級閾值u，若Xi>u，稱Y=Xi-u為震級超出量。則超出量的分布函數(shù)可表示為：對于充分大的震級閾值u，震級超出量Y的分布函數(shù)近似服從廣義帕累托分布：式中：y>0，

廣義帕累托分布與廣義極值分布具有相同的形狀參數(shù)ξ。當(dāng)ξ≥0時，廣義帕累托分布均沒有上邊界；只有當(dāng)ξ<0時，廣義帕累托分布具有有限上界。

1.2 震級上限估計方法

震級上限是指潛在震源區(qū)內(nèi)可能發(fā)生的最大地震的震級，可以認(rèn)為在潛在震源區(qū)內(nèi)未來發(fā)生超過該震級地震的概率幾乎為零（胡聿賢，1999）。當(dāng)選定震級閾值u，且X>u時，則有：令ζu=Pr（X>u），表示為超過閾值的樣本個數(shù)與總樣本的比值，式（5）可寫為：考慮到潛在震源區(qū)的地震事件的震級必為有限值，因此只有當(dāng)形狀參數(shù)ξ<0時，廣義帕累托分布具有有限上界，才符合震級存在上限的條件，震級上限可表示為：利用最大似然估計方法可得到3個參數(shù)（ζu，σ，ξ）的估計值。根據(jù)Delta定理可知，3個參數(shù)的最大似然估計值近似符合正態(tài)分布，因此根據(jù)方差協(xié)方差矩陣可以獲得置信水平為1-α的震級上限的置信區(qū)間如式（8）所示：式中：Var（xN）≈▽xNTV▽xN；Y是（ζu，σ，ξ）參數(shù)估計的協(xié)方差矩陣；數(shù)（Coles，2001）。

2 龍門山地區(qū)地質(zhì)構(gòu)造

第五代《中國地震動參數(shù)區(qū)劃圖》中將龍門山地區(qū)劃分為地震構(gòu)造區(qū)，區(qū)內(nèi)的地震活動具有相似的地震構(gòu)造環(huán)境和發(fā)震構(gòu)造模型，因此，本文選擇龍門山地區(qū)為研究區(qū)。該區(qū)域塊體滑移方向為南南東，區(qū)域主壓應(yīng)力呈北西西一南西西向的水平擠壓，具有發(fā)生強震的構(gòu)造與動力學(xué)環(huán)境。該區(qū)主要包括岷山斷塊和龍門山斷裂帶（圖1），為巴顏喀拉地塊的東邊界，該邊界不僅是強烈的活動構(gòu)造帶，也是我國南北地震帶中段的一部分（Densmore，2007；聞學(xué)澤等，2011；易桂喜等，2011；孔軍，周榮軍，2014；高孟潭，2015）。

岷山斷塊東、西邊界分別為虎牙斷裂和岷江斷裂?；⒀罃嗔阎饕时北蔽飨蜓由?，自晚第四紀(jì)以來表現(xiàn)為逆一走滑運動性質(zhì)，且以水平走滑運動為主。岷江斷裂是劃分西部川青塊體和東部岷山斷塊的界限，自晚更新世以來顯示為逆一走滑運動性質(zhì)，且水平位錯量與垂直位錯量大致相當(dāng)。川青塊體向南東東方向的滑移過程現(xiàn)今仍在繼續(xù)，導(dǎo)致了岷山斷塊區(qū)強烈的近代地震活動（陳智梁等，1998；Chen et al，2000；周榮軍等，2006）。

龍門山斷裂帶是青藏高原東緣邊界山脈，不僅是中國東部和西部地形地貌的界限，還地處我國南北地震帶的中段，地震活動頻度高，強度大，是青藏高原北部地震亞區(qū)主要強震活動帶之一（鄧起東等，1994；李勇等，2006）。該斷裂帶主要由一系列大致平行的北東一南西向的逆沖斷裂組成，由西向東分別為：汶川—茂縣斷裂、映秀—北川斷裂和安縣—灌縣斷裂。汶川—茂縣和映秀一北川斷裂是控制龍門山快速隆升的重要斷裂。龍門山斷裂帶具有長期活動性，自晚三疊世至古近紀(jì)時期，以逆沖-左旋走滑作用為特征，而中新生代以來均顯示由北西向南東的逆沖運動，并具有顯著的右旋走滑分量（周榮軍等，2006；李勇等，2006）。

3 龍門山地區(qū)震級上限估計

3.1 地震活動

本文數(shù)據(jù)來源主要為國家地震科學(xué)數(shù)據(jù)共享中心（http：//data.earthquake.cn）的中國歷史地震目錄數(shù)據(jù)庫以及國家地震臺網(wǎng)地震目錄數(shù)據(jù)庫，黃瑋瓊等（1994）對中國大陸大部分地區(qū)的歷史地震資料作了較詳細(xì)的分析與研究，顯示龍門山區(qū)域MS≥4%地震資料在1931年之后基本完整。因此，筆者收集了1930-2016年龍門山地區(qū)（101°～106°E，29.5°～33.5°N）M≥4.5的歷史地震目錄。經(jīng)統(tǒng)計，1930年至今，龍門山地區(qū)發(fā)生過5次M--7.0強震，其中在龍門山斷裂帶上發(fā)生過2次，分別是2008年汶川8.0級地震和2013年蘆山7.0級地震；北部岷山地區(qū)發(fā)生過3次，分別是1933年疊溪7.5級地震和1976年松潘2次7.2級地震。

在利用廣義帕累托模型進(jìn)行地震危險性估計之前需要盡量刪除地震目錄中的余震。目前國內(nèi)外地震危險性分析研究中，考慮到計算的簡便和工程應(yīng)用上的可操作性，普遍采用Gardner和Knopoff（1974）提出的時空窗法來刪除余震（徐偉進(jìn)，2012）。在該方法中，空間窗口根據(jù)式（9）確定（Knopoff，Gardner，1972）：

lgR=aM+b（9）式中：a，b為固定參數(shù)值，在我國a，b的取值分別為0.5、-1.78（劉杰等，1996）。

時間窗口如表1所示。刪除余震后，共收集到龍門山地震構(gòu)造區(qū)內(nèi)地震事件118個，對收集到的地震資料統(tǒng)計發(fā)現(xiàn)，龍門山地震構(gòu)造區(qū)內(nèi)以淺源地震為主，其中M≥7地震5次，6≤M<7地震14次，5≤M<6地震34次。從圖1可看出，地震震中分布位置與斷裂帶的展布形態(tài)密切相關(guān)。

3.2 震級上限估計

利用廣義帕累托分布估計龍門山地區(qū)震級上限，首先需確定震級閾值u。閾值選取通過2種方法：（1）依據(jù)震級平均超出量分布函數(shù)圖，震級平均超出量與震級閾值應(yīng)滿足線性關(guān)系：（2）隨著閾值選取不斷增大，形狀參數(shù)ξ和修飾的尺度參數(shù)σ估計值應(yīng)基本保持穩(wěn)定不變（Coles，2001）。由圖2可以看出，當(dāng)震級閾值在（5.1，5.5）、（5.7，6.0）以及（6.3，6.5）范圍內(nèi)時，平均超出量近似呈線性分布。為了對閾值進(jìn)行最合理的選取，還需要綜合考慮樣本數(shù)據(jù)的實際情況，閾值選取過大，超出量樣本數(shù)過少，誤差增大。結(jié)合圖3分析，閾值在（5.1，5.5）范圍內(nèi)相應(yīng)的形狀參數(shù)和修飾的尺度參數(shù)的估計值基本保持不變，根據(jù)廣義帕累托分布應(yīng)用條件，選取盡可能大的閾值u=5.5。

在確定震級閾值u后，通過概率圖、分位數(shù)圖、重現(xiàn)水平圖以及概率密度圖可以對所擬合的廣義帕累托分布模型進(jìn)行檢驗。當(dāng)數(shù)據(jù)符合理論假設(shè)分布時，概率圖與分位數(shù)圖應(yīng)近似為直線。重現(xiàn)水平圖可作為模型檢驗的診斷工真，當(dāng)ξ=0時為直線，當(dāng)ξ<0時為凸曲線，當(dāng)ξ>0時為凹曲線。模型診斷如圖4所示。

圖4a是概率圖，橫坐標(biāo)表示實際數(shù)據(jù)的累積概率，縱坐標(biāo)表示選用分布模型的累積概率；圖4b是分位數(shù)圖，橫坐標(biāo)表示所選分布模型的分位數(shù)，縱坐標(biāo)表示實際數(shù)據(jù)利用經(jīng)驗分布函數(shù)計算得到的分位數(shù)。從圖4a，b可看出，大部分點均沿直線分布，因此不能否認(rèn)所擬合的模型。由于ξ的估計值為負(fù)，因此最大震級分布存在有限上界，故重現(xiàn)水平圖中的曲線為凸曲線，漸近地趨于有限值（圖4c）。概率密度圖中的密度曲線估計與直方圖也顯示一定的擬合（圖4d）。因此4個診斷圖均支持?jǐn)M合的廣義帕累托分布模型。

利用最大似然法估計廣義帕累托分布的參數(shù)，得到修飾的尺度參數(shù)和形狀參數(shù)的估計值為（σ，ξ）=（1.31，-0.47），其中形狀參數(shù)估計的標(biāo)準(zhǔn)差為0.18，其95%的置信區(qū)間為[-0.82，-0.12]，形狀參數(shù)的最大似然估計與置信區(qū)間都顯示為負(fù)值，這表明龍門山地區(qū)在研究時間尺度范圍內(nèi)震級超出量分布具有上限值，這也與潛在震源區(qū)一定具有震級上限的認(rèn)識相符。

將閾值、形狀參數(shù)與修飾的尺度參數(shù)估計值代入式（7）得到震級上限為8.3。對震級上限的誤差進(jìn)行分析，首先需計算參數(shù)估計的協(xié)方差矩陣：將式（10）協(xié)方差矩陣代入式（8），得到震級上限的標(biāo)準(zhǔn)差為0.50，由上限估計值和標(biāo)準(zhǔn)差可知震級上限置信度為95%的置信區(qū)間為[7.3，9.3]。

4 結(jié)論與討論

通過對龍門山地區(qū)歷史強震記錄進(jìn)行統(tǒng)計，從擬合結(jié)果可以看出震級超出量分布服從廣義帕累托分布。依據(jù)震級平均超出量分布函數(shù)圖以及參數(shù)估計值變化圖，綜合分析選取5.5為震級閾值，利用最大似然估計方法，獲得修飾的尺度參數(shù)以及形狀參數(shù)的估計值分別為1.31、-0.47。對龍門山地區(qū)強震危險性進(jìn)行分析，計算得到龍門山地區(qū)的震級上限為8.3，置信度為95%的置信區(qū)間為[7.3，9.3]。

在利用廣義帕累托分布計算龍門山地區(qū)震級上限時，用于統(tǒng)計的震級樣本具有相似的地震構(gòu)造環(huán)境和發(fā)震構(gòu)造模型，滿足構(gòu)建廣義帕累托模型時樣本屬于同分布的條件。另外，采用時空窗法相對地減少余震對統(tǒng)計結(jié)果的影響，也更加符合事件相互獨立的條件。雖然時空窗法的經(jīng)驗公式在刪除余震時難免會存在少量誤差，但本文在計算龍門山地區(qū)震級上限時，震級閾值設(shè)定為5.5，余震震級大于5.5的數(shù)量較少，因此，刪除余震的誤差對于計算結(jié)果的影響相對較小。關(guān)于龍門山地區(qū)地震目錄余震刪除方法在模型中的應(yīng)用還需進(jìn)一步的研究與探討。

本文利用廣義帕累托分布模型計算出的震級上限為8.3，高于第五代《中國地震動參數(shù)區(qū)劃圖》中龍門山構(gòu)造區(qū)的震級上限8.0。廣義帕累托分布可以充分利用歷史地震大震級區(qū)段的數(shù)據(jù)資料，在小震資料不充足的情況下也可構(gòu)建潛源區(qū)的地震活動模型，其統(tǒng)計結(jié)果可作為地震區(qū)劃以及地震危險性分析研究工作的參考。

本文在撰寫過程中得到蔣長勝研究員、邊銀菊研究員以及孔韓東博士的幫助，在此表示衷心感謝。

參考文獻(xiàn)：

陳智梁，劉宇平，趙濟湘，等，1998.青藏高原東部地殼運動的GPS測量[J].中國地質(zhì)，（5）：32-35.

鄧起東，陳社發(fā)，趙小麟.1994.龍門山及其鄰區(qū)的構(gòu)造和地震活動及動力學(xué)[J].地震地質(zhì)，16（4）：389-403.

鄧起東.2007.中國括動構(gòu)造圖（1：400萬）[M].北京：地震出版社.

高孟潭.2015.GB 18306-2015《中國地震動參數(shù)區(qū)劃圖》宣貫教材[M].北京：中國質(zhì)檢出版社.

胡聿賢.1999.地震安全性評價技術(shù)教程[M].北京：地震出版社，183-193.

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

孔軍，周榮軍.2014.龍門山和成都地震構(gòu)造區(qū)的劃分[J].震災(zāi)防御技術(shù)，9（1）：64-73.

李海兵，司家亮，付小方，等.2009.2008年汶川地震同震滑移特征、最大滑移量及構(gòu)造意義[J].第四紀(jì)研究，29（3）：387-402.

李勇，周榮軍，Denstaore A L，等.2006.青藏高原東緣龍門山晚新生代走滑—逆沖作用的地貌標(biāo)志[J].第四紀(jì)研究，26（1）：40-51.

劉杰，陳棋福，陳颙.1996.華北地區(qū)地震目錄完全性分析[J].地震，（1）：59-67.

潘華，高孟潭，李金臣.2009.新版美國地震區(qū)劃圖源及其參數(shù)模型的分析與評述[J].震災(zāi)防御技術(shù)，4（2）：131-140.

潘華，李金臣.2016.新-代地震區(qū)劃圖的地震活動性模型[J].城市與減災(zāi)，108（3）：28-33.

錢小仕，蔡曉光，任晴晴.2013.中國大陸活動地塊邊界帶強震震級分布特征研究[J].地震工程與工程振動，33（1）：212-220.

史道濟.2006.實用極值統(tǒng)計方法[M].天津：天津科學(xué)技術(shù)出版社.

田建偉，劉哲，任魯川.2017.基于廣義帕累托分布的馬尼拉海溝俯沖帶地震危險性估計[J].地震，37（1）：158-165.

聞學(xué)澤，杜方，張培震，等.2011.巴顏喀拉塊體北和東邊界大地震序列的關(guān)聯(lián)性與2008年汶川地震[J].地球物理學(xué)報，54（3）：706-716.

胥廣銀.2003.潛在震源三維空間模型及其在地震危險性概率分析中的應(yīng)用研究[D].北京：中國地震局地球物理研究所.

徐偉進(jìn)，高孟潭.2012.根據(jù)截斷的G-R模型計算東北地震區(qū)震級上限[J].地球物理學(xué)報，55（5）：1710-1717.

徐偉進(jìn).2012.地震危險性分析中地震時空統(tǒng)計分布模型研究[D].北京：中國地震局地球物理研究所.

易桂喜，聞學(xué)澤，辛華，等.2011.2008年汶川MS8.0地震前龍門山—岷山構(gòu)造帶的地震活動性參數(shù)與地震視應(yīng)力分布[J].地球物理學(xué)報，54（6）：1490-1500.

周榮軍，李勇，Alexander，等，2006.青藏高原東緣活動構(gòu)造[J].礦物巖石，26（2）：40-51.

Balkema A A，Haan L D.1974.Residual Life Time at Great Age[J].An-nals of Probability，2（5）：792-804.

Chen Z，Burchfiel B C，Liu Y，et al.2000.Global Positioning Systemmeasurements from eastern Tibet and their implications for India/Eurasia intercontinental deformation[J].Journal of Geophysical Be-search Solid Earth，105（B7）：16215-16227.

Coles S.2001.An Introduction to Statistical Modeling of Extreme Values[M].Berlin：Springer-Verlag.

Comell C A.196$.Engineering seismic risk analysis[J].Bulletin of theSeismological Society of America，58（5）：1583-1606：

Densmore A L，EllsMa，Li Y，et al.2007.Active tectonics of the Bei-chuan and Pengguan faults at the easternMargin of the Tibetan Plat-eau[J].Tectonics，26（4），TC4005.

Fisher R A，Tippett L H C.1928.Limiting forMS of the frequency distribu-tion of the largest or smallest member of a sample[J].MathematicalProceedings of the Cambridge Philosophical Society，24（2）：180-190.

Gardner J K，Knopolf L.1974.Is the sequence of earthquakes in southernCalifornia，with aftershocks removed，Poissonian？[J].Bulletin ofthe Seismological Society of America，64（5）：1363-1367.

Huyse L，Chen R，Stamatakos J A.2010.Application of Generalized Pare-to Distribution to Constrain Uncertainty in Peak Ground Accelera-tions [J].Bulletin of the Seismological Society of America，100（1）：87-101.

Knopoff L.Gardner J K.1972.Higher seismic activity during local nighton the raw worldwide earthquake catalogue[J].Geophysical Journalof the Royal Astronomical Society，28（3）：311-313.

Pickands J.1975.Statistical Inference Using Extreme Order Statistics[J].Annals of Statistics，3（1）：119-131.

Pisarenko V F，Sornette D.2003.Characterization of the Frequency of Ex-treme Earthquake Events by the Generalized Pareto Distribution[J].Pure & Applied Geophysics，160（12）：2343-2364.