基于輪廓似然估計(jì)的廣義極值分布在地震中長(zhǎng)期預(yù)測(cè)中的應(yīng)用＊

趙宜賓張艷芳王福昌任晴晴

（中國(guó)河北三河 065201 防災(zāi)科技學(xué)院）

引言

極端事件在隨機(jī)現(xiàn)象的研究過(guò)程中出現(xiàn)頻率很低，比如大地震、洪水、颶風(fēng)等，但這些事件一旦發(fā)生，將會(huì)給生產(chǎn)或生活帶來(lái)巨大的影響，因此對(duì)數(shù)據(jù)樣本中的極端變異性數(shù)據(jù)信息進(jìn)行深入研究成為眾多學(xué)者的共識(shí)，在此基礎(chǔ)上形成的極值統(tǒng)計(jì)分析方法成為災(zāi)害損失評(píng)估與風(fēng)險(xiǎn)評(píng)價(jià)的重要工具.

von Bortkiewicz （1922）是近代第一個(gè)明確提出極值問(wèn)題的學(xué)者，其在研究正態(tài)分布樣本的極差時(shí)，發(fā)現(xiàn)來(lái)自正態(tài)分布的樣本最大值具有新的分布；von Mises （1923）首次對(duì)正態(tài)樣本極值的漸近分布進(jìn)行了研究；Fréchet （1927）的研究表明，來(lái)自不同分布但有某種共同性質(zhì)的最大值可以有相同的漸近分布；Fisher和Tippett （1928）將次序統(tǒng)計(jì)量規(guī)范化，按中心極限定理的思想，提出了極值分析漸近原理的基礎(chǔ)性定理，完成了極值統(tǒng)計(jì)分析基礎(chǔ)理論的搭建；Gnedenko （1943）對(duì)極值分布的漸近原理給出了嚴(yán)格的證明.在此基礎(chǔ)上，眾多學(xué)者從理論分析（Jenkinson，1955；De Haan，1970，1971）與實(shí)際應(yīng)用（史道濟(jì)，2006；Bhunyaet al，2012）兩個(gè)方面開(kāi)展了對(duì)極值理論的大量研究.

最具破壞力的自然災(zāi)害之一的地震，尤其是超強(qiáng)震，對(duì)社會(huì)經(jīng)濟(jì)造成的損失是巨大的，因此對(duì)特定區(qū)域作地震風(fēng)險(xiǎn)估計(jì)對(duì)于宏觀決策和微觀管理都是必要的.Nordquist （1945）將極值理論引入到震級(jí)預(yù)測(cè)過(guò)程；多名研究人員（Epstein，Lomnitz，1966；陳培善，林邦慧，1973；Yegulalp，Kuo，1974）對(duì)最大震級(jí)預(yù)測(cè)模型作了系統(tǒng)的研究；高盂潭和賈素娟（1988）對(duì)極值理論在工程抗震中的應(yīng)用作了詳細(xì)論述；陳虹和黃忠賢（1995）及錢小仕等（2012，2013）給出了基于統(tǒng)計(jì)分析的地震危險(xiǎn)性評(píng)價(jià)指標(biāo)及其計(jì)算方法.

對(duì)于廣義極值模型，對(duì)模型分類和統(tǒng)計(jì)規(guī)律的表達(dá)起主要作用的是形狀參數(shù)，其數(shù)值變化是研究者比較關(guān)心的.極值分布主要用極大似然估計(jì)和矩估計(jì)法確定參數(shù)取值，對(duì)于參數(shù)的區(qū)間估計(jì)，特別是地震重現(xiàn)水平的區(qū)間估計(jì)，傳統(tǒng)的Delta法（Rao，1965）得到的對(duì)稱區(qū)間，隨著震級(jí)的提高，與預(yù)測(cè)結(jié)果在大于重現(xiàn)水平取值時(shí)有更大的不確定性的實(shí)際情況不符.為此，一些科研專家嘗試用輪廓似然估計(jì)代替極大似然估計(jì)，用以確定極值分布的參數(shù)值.Murphy和Van der Vaart （2000）證明了在通常情況下輪廓似然估計(jì)與極大似然估計(jì)是等價(jià)的.Tajvidi （2003）、Gilli和 K?llezi （2006）及魯帆等（2013）對(duì)各類實(shí)際問(wèn)題的風(fēng)險(xiǎn)評(píng)估，證明了輪廓似然區(qū)間估計(jì)的不對(duì)稱性恰是對(duì)極值變量的不確定性隨取值增加而變大的較好表達(dá).

綜上所述，在基于廣義極值分布對(duì)地震危險(xiǎn)性進(jìn)行評(píng)價(jià)的過(guò)程中，本文擬用輪廓似然函數(shù)進(jìn)行參數(shù)估計(jì)，以便更合理地表達(dá)強(qiáng)震預(yù)測(cè)的不確定性.

1 極值分布相關(guān)理論簡(jiǎn)述

作為風(fēng)險(xiǎn)管理和安全評(píng)價(jià)的重要分析工具，極值分析主要針對(duì)隨機(jī)變量中的極端變異性數(shù)據(jù)進(jìn)行系統(tǒng)建模，進(jìn)而對(duì)極端事件在未來(lái)發(fā)生的可能性給出預(yù)測(cè)評(píng)價(jià).為了更清晰地表述基于輪廓似然估計(jì)的廣義極值分布模型的構(gòu)建過(guò)程，對(duì)重要理論和概念概述如下：

1.1 廣義極值分布的基本理論

極值分析的基礎(chǔ)定理（Fisher，Tippett，1928）：若X1，X2，···，Xn，···是分布函數(shù)為F（x）的獨(dú)立同分布的隨機(jī)變量序列，對(duì)于任意n∈Z＋，令Yn＝max{X1，X2，···，Xn}，如果存在常數(shù)列{αn＞0}和{βn}，以及非退化函數(shù)Λ（x），使得

成立，則稱Λ（x）為極值分布，并稱F（x）屬于極值分布Λ（x）的最大值吸引場(chǎng)，記為F（x）∈MDA（Λ）.

吸引場(chǎng)的概念說(shuō)明相互獨(dú)立且服從相同分布函數(shù)F（x）的隨機(jī)變量序列，在滿足一定的條件下，序列的區(qū)組最大值近似服從廣義極值（generalized extreme value，縮寫為GEV）分布.對(duì)于驗(yàn)證條件αn和βn難于確定的問(wèn)題，F(xiàn)réchet （1927）證明了區(qū)組最大值分布規(guī)律是穩(wěn)定的：即使隨機(jī)變量分布不同，若區(qū)組最大值的近似分布存在，則其服從GEV分布，這是極值理論應(yīng)用的基礎(chǔ).

Jenkinson （1955）將Λ（x）的三種極值分布形式統(tǒng)一為一個(gè)表達(dá)式，稱之為廣義極值分布.具有位置參數(shù)和尺度參數(shù)的廣義極值分布的分布函數(shù)記為

式中，ξ為形狀參數(shù)，μ為位置參數(shù)，σ為尺度參數(shù).

GEV 分布的類型由形狀參數(shù)ξ的符號(hào)決定：當(dāng)ξ＞0時(shí)，Gev（x，ξ，μ，σ）表示 Fréchet分布；當(dāng)ξ＝0時(shí)，Gev（x，ξ，μ，σ）表示耿貝爾（Gumbel）分布；當(dāng)ξ＜0時(shí)，Gev（x，ξ，μ，σ）表示具有有限上端點(diǎn)的威布爾（Weibull）分布.推導(dǎo)過(guò)程史道濟(jì)（2006）一文中有詳細(xì)描述.

應(yīng)用極值理論進(jìn)行安全評(píng)價(jià)時(shí)，常用到分位數(shù)概念，定義如下：

設(shè)隨機(jī)變量X∽F（x），若xp滿足條件

則稱xp為F（x）的p分位數(shù).

從分位數(shù)的定義可知：xp＝F?1（p），因此，由 GEV 分布函數(shù)（式 2）可得

若ξ＜0，令p→1可得GEV分布理論的上端點(diǎn)

為了盡可能地滿足樣本之間相互獨(dú)立的條件，對(duì)時(shí)間序列進(jìn)行區(qū)組劃分時(shí)，應(yīng)設(shè)置足夠長(zhǎng)的時(shí)間間隔，以每個(gè)時(shí)間段中最大值序列作為GEV建模樣本，最大限度地滿足樣本的漸近獨(dú)立性.在此基礎(chǔ)上，進(jìn)行模型參數(shù)估計(jì)和適用性檢驗(yàn)，最后利用模型進(jìn)行安全性評(píng)價(jià).

1.2 GEV分布參數(shù)的輪廓似然估計(jì)

GEV分布的形狀參數(shù)ξ可以確定模型的分類及密度曲線形狀，而分位數(shù)是安全評(píng)價(jià)的重要指標(biāo)，本節(jié)將利用輪廓似然估計(jì)給出這兩個(gè)重要參數(shù)的估計(jì).

設(shè)隨機(jī)變量X∽f（x，θ），θ∈Θ，其中概率密度f(wàn)（x，θ）的形式已知，θ＝（θ1，θ2，···，θk）包含k個(gè)未知參數(shù)，Θ為θ的可能取值范圍.X1，X2，···，Xk是樣本，x1，x2，···，xk是樣本值，則（X1，X2，···，Xk）的聯(lián)合概率密度在（x1，x2，···，xk）的函數(shù)值L（θ）稱為樣本的似然函數(shù).

取 Θ0?Θ，令和分別為θ在Θ0和Θ中的極大似然估計(jì)值，稱λ為對(duì)數(shù)似然比值，相應(yīng)的λ（X）稱為對(duì)數(shù)似然比統(tǒng)計(jì)量.史寧中（2008）的研究表明，在一定條件下，當(dāng)Θ0＝{θ0} 時(shí) ，有＝θ0，此時(shí)有），即對(duì)數(shù)似然比近似服從自由度為k的χ2分布.這是輪廓似然估計(jì)法進(jìn)行參數(shù)區(qū)間估計(jì)的理論基礎(chǔ).

若將θ中的未知參數(shù)分成兩類， 則似然函數(shù)寫作L（），其中θi是研究者重點(diǎn)關(guān)注的分量，相應(yīng)的是θ的其它未知分量.則θi的輪廓似然函數(shù)定義為即取定θi時(shí)，函數(shù)值Lp（θi）是L）的最大值.

基于上述理論，在利用GEV分布進(jìn)行安全評(píng)價(jià)過(guò)程中主要處理ξ≠0的情況，所以下文僅對(duì)ξ≠0的情況進(jìn)行表述，ξ＝0的推導(dǎo)過(guò)程與之類似.

由GEV的概率密度公式

可得，GEV的對(duì)數(shù)似然函數(shù)為

因此GEV關(guān)于形狀參數(shù)ξ的輪廓似然函數(shù)可表示為即對(duì)于在可能取值范圍內(nèi)ξ的每個(gè)值，其函數(shù)值Lp（ξ）取在μ和σ定義范圍內(nèi)，使得lnL（ξ，μ，σ）取得最大值，即GEV輪廓似然函數(shù)對(duì)應(yīng)的點(diǎn)集為

由此可得ξ的輪廓似然估計(jì)值為

因?yàn)閷?duì)數(shù)似然比近似服從χ2分布，所以當(dāng)ξ＝ξi成立時(shí)，λ（ξ）＝2{Lp－Lp（ξ）}∽χ2（1），進(jìn)而可得ξ的置信水平為1－α的置信區(qū)間為

為了求分位數(shù)的置信區(qū)間，必須在似然函數(shù)中引入?yún)?shù)xp.由分位數(shù)的計(jì)算（式4）可得μ＝xp＋σ/ξ［1－（?lnp）?ξ］，將μ代入到對(duì)數(shù)似然表達(dá)式（式7）得lnL（ξ，xp，σ）.因此，GEV關(guān)于xp的輪廓似然函數(shù)即對(duì)于在可能取值范圍內(nèi)xp的每個(gè)值，其函數(shù)值Lp（xp）取在ξ，σ定義范圍內(nèi)，使得 lnL（ξ，xp，σ）取得最大值，即輪廓似然函數(shù)對(duì)應(yīng)的點(diǎn)集為

由此可得xp的輪廓似然估計(jì)值為

類似于ξ，xp的置信水平為1－α的置信區(qū)間為

1.3 地震發(fā)生的重現(xiàn)期與重現(xiàn)水平

假設(shè)X1，X2，···，Xn為某一特定區(qū)域以半年為間隔的地震震級(jí)最大值樣本，在進(jìn)行余震刪除工作后，鑒于時(shí)間間隔比較長(zhǎng)，可認(rèn)為樣本滿足漸近獨(dú)立性，且服從GEV分布.設(shè)第一次出現(xiàn)超閾值u的時(shí)間為τ1＝min{k：Xk＞u}，令P{X＞u}＝1－Gevξ（u）＝q，則有

另由幾何分布的性質(zhì)可得：相鄰兩次超閾值的時(shí)間與第一次出現(xiàn)超閾值的時(shí)間理論上是相等的，這是危險(xiǎn)事件重現(xiàn)期應(yīng)具備的基本屬性.依上述討論，重現(xiàn)期定義為給定時(shí)間序列X1，X2，···，Xn及閾值u，若序列中變量第一次出現(xiàn)超閾值u的平均時(shí)間為T（u），則稱u為重現(xiàn)水平，相應(yīng)的T（u）稱為重現(xiàn)期.由上述分析可得重現(xiàn)水平為u的重現(xiàn)期為

相應(yīng)地，給定重現(xiàn)期T，求重現(xiàn)期的反函數(shù)，可得相應(yīng)的重現(xiàn)最大震級(jí)（重現(xiàn)水平）為

實(shí)際上重現(xiàn)期為T的重現(xiàn)水平u（T）就是GEV分布的1 －1/T分位數(shù)，即u（T）＝x1－1/T.重現(xiàn)期與重現(xiàn)水平是進(jìn)行地震預(yù)報(bào)分析的兩個(gè)重要指標(biāo)，也是進(jìn)行地震應(yīng)急預(yù)案制定的重要參考因素.下文將以輪廓似然估計(jì)法為工具，對(duì)東昆侖地震帶的地震危險(xiǎn)性進(jìn)行綜合分析.

2 基于輪廓似然估計(jì)的東昆侖地震帶地震危險(xiǎn)性分析

2.1 地震數(shù)據(jù)樣本完整性分析

依據(jù)對(duì)大陸活動(dòng)地塊劃定的東昆侖地震帶區(qū)域邊界（張國(guó)民等，2005），以從國(guó)家地震科學(xué)數(shù)據(jù)共享中心獲取的正式地震目錄為源數(shù)據(jù)，提?。?2.8°E——104.2°E，33.5°N——37.1°N）自公元前193年至2019年12月的4 385條記錄作為研究樣本.為了最大程度地滿足地震記錄樣本的相互獨(dú)立性，采用震級(jí)相關(guān)時(shí)空窗法（陳凌等，1998）對(duì)研究樣本中的余震進(jìn)行剔除，得到東昆侖地震帶的震例數(shù)據(jù)3 616條.本文對(duì)地震大小的描述采用面波震級(jí)MS，如果面波震級(jí)缺失，通過(guò)公式MS＝1.13ML－1.08 （汪素云等，2010）將近震震級(jí)ML轉(zhuǎn)換為MS，進(jìn)行數(shù)據(jù)填充，進(jìn)而得到用于地震危險(xiǎn)性分析的完整數(shù)據(jù)樣本.

東昆侖地震帶的地震空間分布如圖1a所示，整個(gè)地震帶內(nèi)斷裂帶密集排布，在高海拔地區(qū)大體沿西北向東南展布.地震帶內(nèi)發(fā)震點(diǎn)的分布呈現(xiàn)東北部頻率高、震級(jí)低，中部和西南部頻率低、震級(jí)高的特點(diǎn).MS7.0以上強(qiáng)震基本上都發(fā)生在斷裂帶上，而在地殼隆起邊緣的斷裂帶上，地震發(fā)生頻率較高.

圖1 東昆侖地震帶的地震分布規(guī)律（a） 地震空間分布；（b） M-t圖Fig.1 Distribution law of earthquakes in East Kunlun seismic zone（a） Spatial distribution of earthquakes；（b） M-t diagram

東昆侖地震帶MS5.0以上地震從1937年開(kāi)始記錄（黃瑋瓊等，1994），MS6.0上地震從1926年開(kāi)始記錄，1970年中國(guó)地震臺(tái)網(wǎng)建立后，發(fā)震情況記錄則比較完整.1930——2019年的MS3.0以上震級(jí)與時(shí)間關(guān)系如圖1b所示.M-t圖表明東昆侖地震帶MS4.0和MS3.0以上地震分別從1950年和1965開(kāi)始才有相對(duì)完整的記錄.為了盡可能地保留數(shù)據(jù)信息的同時(shí)滿足模型分析對(duì)數(shù)據(jù)連續(xù)性的要求，本文以半年為區(qū)組間隔，提取了1950年后的震級(jí)最大值為地震危險(xiǎn)性分析樣本.

2.2 基于輪廓似然估計(jì)的東昆侖地震帶GEV分布模型構(gòu)建

按1.2節(jié)GEV分布參數(shù)的輪廓似然估計(jì)理論，以Matlab軟件為計(jì)算平臺(tái)，以遺傳算法作為數(shù)值計(jì)算的主要工具，按如下步驟估計(jì)GEV分布的形狀參數(shù)ξ：首先，設(shè)定ξ的可能取值范圍為［ξl，ξu］ ， 對(duì)于任意ξi∈ ［ξl，ξu］ ， 依據(jù)式（8），搜索μi，σi，使得并將滿足條件的點(diǎn)（ξi，μi，σi，Lp（ξi））加入輪廓似然估計(jì)點(diǎn)集PLξ；其次，依據(jù)式（9）確定i0，由ξ的輪廓似然估計(jì)值ξi0和相應(yīng)的μi0，σi0，確定GEV分布函數(shù)GEV（x；ξi0，μi0，σi0）；最后，依據(jù)式（11）確定ξ的置信水平為1－α的置信區(qū)間.

根據(jù)上述步驟得到的輪廓似然估計(jì)點(diǎn)集PLξ如圖2所示.圖2表明輪廓似然函數(shù)與形狀參數(shù)ξ存在類似拋物線的關(guān)系，在ξ＝?0.204 0時(shí)輪廓似然函數(shù)取得最大值，此時(shí)μ＝0.847 5，σ＝4.834 5.進(jìn)而可得，GEV 的概率密度函數(shù)中參數(shù)（ξ，μ，σ）的估計(jì)值為（?0.204 0，0.847 5，4.834 5）.

圖2 形狀參數(shù)與輪廓對(duì)數(shù)似然函數(shù)之間的關(guān)系Fig.2 Relationship between shape parameters and profile log likelihood function

基于輪廓似然估計(jì)構(gòu)建的GEV分布模型對(duì)樣本數(shù)據(jù)分析的適用性檢驗(yàn)如圖3所示.直方圖輪廓趨勢(shì)線與概率密度曲線基本吻合；P-P圖的散點(diǎn)在56°線附近小幅波動(dòng)，表明樣本與理論分布的分位數(shù)特征基本是一致的；且ξ?＜0表示震級(jí)的重現(xiàn)水平有理論上限.上述檢驗(yàn)表明，本文構(gòu)建的GEV分布模型適用于對(duì)昆侖山地震帶作地震危險(xiǎn)性分析.

圖3 GEV 分布模型適應(yīng)性檢驗(yàn)圖（a） 密度曲線與直方圖；（b） P-P 檢驗(yàn)Fig.3 Adaptability test of GEV distribution model（a） Density curves and histograms；（b） P-P test

由于ξ的輪廓置信上限小于0，則依據(jù)式（10）可求震級(jí)的理論上限.同時(shí)，在區(qū)組震級(jí)X∽GEV （x；ξi0，μi0，σi0）的條件下，利用E（X）＝μi0－σi［01－Γ（1－ξi0）］/ξi0，ξi0＜1，可求得最大震級(jí)理論平均值，其中Γ（x）是伽馬函數(shù).

為了比較參數(shù)估計(jì)方法的差異，表1分別列出了輪廓似然法和極大似然法（錢小仕等，2012）的估計(jì)結(jié)果，在進(jìn)行模型主要參數(shù)估計(jì)時(shí)，兩種估計(jì)方法的效果基本相同.

表1 輪廓似然估計(jì)與極大似然估計(jì)GEV的分布結(jié)果對(duì)比Table 1 Comparation of the profile likelihood estimation and maximum likelihood estimationof GEV distribution

通過(guò)上述分析可知，東昆侖地震帶每半年的最大震級(jí)平均約為MS5.2，理論上的最大震級(jí)約為MS9.0，屬于巨震級(jí)別，說(shuō)明這一區(qū)域的地質(zhì)條件很不穩(wěn)定，屬于大地震發(fā)生危險(xiǎn)性極高的地區(qū).

2.3 東昆侖地震帶重現(xiàn)水平及置信區(qū)間的輪廓似然估計(jì)

上文估計(jì)的半年最大震級(jí)理論均值和上限，是對(duì)地區(qū)發(fā)震情況的簡(jiǎn)要概述，對(duì)應(yīng)急預(yù)案制定參考意義不是很大.重現(xiàn)期和重現(xiàn)水平是評(píng)估地震危險(xiǎn)性的兩個(gè)核心指標(biāo)，在估計(jì)了GEV分布重要參數(shù)后，給定重現(xiàn)期T，按式（17）可以確定理論重現(xiàn)水平u（T）.為了滿足決策需要，有時(shí)需要求出重現(xiàn)水平的置信區(qū)間.利用信息矩陣獲得的置信區(qū)間（Rao，1965）是關(guān)于置信水平對(duì)稱的，但實(shí)際上隨著預(yù)測(cè)震級(jí)的提高，在高于置信水平的時(shí)候，震級(jí)選擇會(huì)更加離散，也就是說(shuō)置信區(qū)間關(guān)于置信水平對(duì)稱是不合理的.為此，下文將利用輪廓似然法來(lái)估計(jì)重現(xiàn)水平的置信區(qū)間.

在給定重現(xiàn)期的情況下，確定重現(xiàn)水平及置信區(qū)間可按如下步驟進(jìn)行：首先，對(duì)于給定重現(xiàn)期T，按式（17）可以確定理論重現(xiàn)水平u（T），根據(jù)u（T）值設(shè)定一個(gè)相對(duì)保守的重現(xiàn)水平取值范圍［XPl，XPu］.對(duì)于 任 意xpi∈［XPl，XPu］ ， 依據(jù) 式 （12），搜索ξi，σi， 使得L（pxpi）＝并將滿足條件的點(diǎn)（ξi，xpi，σi，Lp（xpi））加入輪廓估計(jì)點(diǎn)集PLxp；其次，依據(jù)式（13）確定i0，得到xp的輪廓似然估計(jì)值xpi0；最后，依據(jù)式（14）確定xp的置信水平為1－α的輪廓置信區(qū)間.

根據(jù)上述步驟得到的重現(xiàn)期分別為20年、50年、100年、500年的重現(xiàn)水平和置信區(qū)間如圖4所示.

圖4 重現(xiàn)水平及置信區(qū)間的輪廓似然估計(jì)（a） 20 年重現(xiàn)期；（b） 50 年重現(xiàn)期；（c） 100 年重現(xiàn)期；（d） 500 年重現(xiàn)期Fig.4 The reappearance level and confidence interval of the profile likelihood estimation（a） 20-year return period；（b） 50-year return period；（c）100-year return period；（d） 500-year return period

用Delta法和輪廓似然估計(jì)得到的估計(jì)結(jié)果對(duì)比分析詳見(jiàn)表2.估計(jì)結(jié)果從數(shù)值計(jì)算的角度佐證了對(duì)于重現(xiàn)水平的估計(jì)，輪廓似然法與極大似然法是等效的（Murphy，Vander Vaart，2000）.

表2的對(duì)比結(jié)果說(shuō)明，對(duì)于重現(xiàn)期10年以下的重現(xiàn)水平置信區(qū)間的估計(jì)，輪廓似然法與極大似然法的估計(jì)結(jié)果基本相同，即針對(duì)短期地震危險(xiǎn)性評(píng)估，兩種方法是等效的.但在進(jìn)行中長(zhǎng)期地震危險(xiǎn)性分析時(shí)，輪廓似然估計(jì)得到的估計(jì)區(qū)間整體向右偏移，即置信下限和上限比相應(yīng)的極大似然估計(jì)要高，而且重現(xiàn)水平右側(cè)的寬度隨著重現(xiàn)期的提高，與左側(cè)區(qū)間寬度之比越來(lái)越大，說(shuō)明隨著重現(xiàn)期的變長(zhǎng)，重現(xiàn)水平隨之提高，而且發(fā)生超重現(xiàn)水平的震級(jí)取值更分散，是對(duì)震級(jí)區(qū)間更為保守的預(yù)測(cè)，也是與實(shí)際情況相吻合的.兩種方法估計(jì)結(jié)果的直觀差異如圖5所示.

圖5 重現(xiàn)水平的輪廓似然估計(jì)與極大似然估計(jì)對(duì)比Fig.5 Comparation of the reproduction level between profile likelihood estimation and maximum likelihood estimation

表2 極大似然估計(jì)與輪廓似然估計(jì)的重現(xiàn)水平對(duì)比Table 2 Comparation of the recurrence level between profile likelihood estimation and maximum likelihood estimation

本文所用源數(shù)據(jù)是東昆侖地震帶截止到2019年12月的數(shù)據(jù)，為了對(duì)比兩種方法的分析效果，首先截取1920年后100年的震例數(shù)據(jù)，其中大于極大似然估計(jì)下限（MS7.20）和輪廓估計(jì)下限（MS7.29）的樣本數(shù)都是5個(gè)，分別為1924年和1973年的MS7.3、1937年和1997年的MS7.5以及2001年的MS8.1，其中MS8.1地震已經(jīng)超過(guò)極大似然估計(jì)上限MS7.95，但在輪廓似然估計(jì)上限之內(nèi).如果估計(jì)精度為MS0.1，則大于輪廓似然估計(jì)下限的樣本為3個(gè)，從數(shù)量來(lái)說(shuō)，優(yōu)于極大似然估計(jì).如果截取1420年后的500年的地震數(shù)據(jù)，大于極大似然估計(jì)下限（MS7.48）的樣本數(shù)為3個(gè)，分別是1937年和1997年的MS7.5，以及2001年的MS8.1，而超過(guò)輪廓似然估計(jì)下限（MS7.63）的樣本只有1個(gè)MS8.1震例，重現(xiàn)數(shù)量上優(yōu)于極大似然估計(jì).以上實(shí)例分析表明，輪廓似然估計(jì)作為震級(jí)重現(xiàn)水平區(qū)間估計(jì)的理論工具，比以信息矩陣為基礎(chǔ)的Delta法更加合理和適用.

3 討論與結(jié)論

本文對(duì)輪廓似然估計(jì)法在廣義極值（GEV）分布模型的參數(shù)估計(jì)中的應(yīng)用從理論分析到數(shù)值計(jì)算進(jìn)行了詳細(xì)地闡述，并將構(gòu)建的極值分布模型應(yīng)用于東昆侖地震帶的地震預(yù)報(bào).在對(duì)參數(shù)的點(diǎn)估計(jì)過(guò)程中，輪廓似然估計(jì)與極大似然估計(jì)效果相當(dāng)；在進(jìn)行重現(xiàn)水平置信區(qū)間估計(jì)過(guò)程中，如果重現(xiàn)期比較短，兩種方法估計(jì)效果幾乎無(wú)差異，但在進(jìn)行中長(zhǎng)期地震預(yù)報(bào)時(shí)，輪廓似然估計(jì)法得到的重現(xiàn)水平置信區(qū)間對(duì)于不確定性的表達(dá)更加充分，較基于信息矩陣的Delta法得到的對(duì)稱置信區(qū)間更為合理.

基于輪廓似然估計(jì)的GEV分布模型能夠?qū)Φ卣鹞ｋU(xiǎn)性作出相對(duì)比較客觀的評(píng)價(jià)，但該模型所用數(shù)據(jù)為區(qū)組最大值信息，在模型構(gòu)建過(guò)程中對(duì)觀測(cè)信息利用不夠充分，使得由模型得到的預(yù)測(cè)結(jié)果與實(shí)際情況存在偏差，作者后續(xù)將主要從歷史信息挖掘和模型調(diào)整與優(yōu)化兩個(gè)方向入手，以構(gòu)建更加有效的地震預(yù)報(bào)預(yù)測(cè)模型.