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

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由GEV的概率密度公式

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

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

由此可得ξ的輪廓似然估計值為

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

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

由此可得xp的輪廓似然估計值為

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

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

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

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

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

2 基于輪廓似然估計的東昆侖地震帶地震危險性分析

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

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

東昆侖地震帶的地震空間分布如圖1a所示，整個地震帶內(nèi)斷裂帶密集排布，在高海拔地區(qū)大體沿西北向東南展布.地震帶內(nèi)發(fā)震點的分布呈現(xià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年開始記錄（黃瑋瓊等，1994），MS6.0上地震從1926年開始記錄，1970年中國地震臺網(wǎng)建立后，發(fā)震情況記錄則比較完整.1930——2019年的MS3.0以上震級與時間關(guān)系如圖1b所示.M-t圖表明東昆侖地震帶MS4.0和MS3.0以上地震分別從1950年和1965開始才有相對完整的記錄.為了盡可能地保留數(shù)據(jù)信息的同時滿足模型分析對數(shù)據(jù)連續(xù)性的要求，本文以半年為區(qū)組間隔，提取了1950年后的震級最大值為地震危險性分析樣本.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 重現(xiàn)水平及置信區(qū)間的輪廓似然估計（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法和輪廓似然估計得到的估計結(jié)果對比分析詳見表2.估計結(jié)果從數(shù)值計算的角度佐證了對于重現(xiàn)水平的估計，輪廓似然法與極大似然法是等效的（Murphy，Vander Vaart，2000）.

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

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

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

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

3 討論與結(jié)論

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

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