基于輪廓似然估計的廣義極值分布在地震中長期預測中的應用＊

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

（中國河北三河 065201 防災科技學院）

引言

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

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

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

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

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

1 極值分布相關理論簡述

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

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

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

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

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

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

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

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

應用極值理論進行安全評價時，常用到分位數概念，定義如下：

設隨機變量X∽F（x），若xp滿足條件

則稱xp為F（x）的p分位數.

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

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

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

1.2 GEV分布參數的輪廓似然估計

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

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

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

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

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

由GEV的概率密度公式

可得，GEV的對數似然函數為

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

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

因為對數似然比近似服從χ2分布，所以當ξ＝ξi成立時，λ（ξ）＝2{Lp－Lp（ξ）}∽χ2（1），進而可得ξ的置信水平為1－α的置信區間為

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

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

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

1.3 地震發生的重現期與重現水平

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

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

相應地，給定重現期T，求重現期的反函數，可得相應的重現最大震級（重現水平）為

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

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

2.1 地震數據樣本完整性分析

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

東昆侖地震帶的地震空間分布如圖1a所示，整個地震帶內斷裂帶密集排布，在高海拔地區大體沿西北向東南展布.地震帶內發震點的分布呈現東北部頻率高、震級低，中部和西南部頻率低、震級高的特點.MS7.0以上強震基本上都發生在斷裂帶上，而在地殼隆起邊緣的斷裂帶上，地震發生頻率較高.

圖1 東昆侖地震帶的地震分布規律（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年中國地震臺網建立后，發震情況記錄則比較完整.1930——2019年的MS3.0以上震級與時間關系如圖1b所示.M-t圖表明東昆侖地震帶MS4.0和MS3.0以上地震分別從1950年和1965開始才有相對完整的記錄.為了盡可能地保留數據信息的同時滿足模型分析對數據連續性的要求，本文以半年為區組間隔，提取了1950年后的震級最大值為地震危險性分析樣本.

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

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

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

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

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

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

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

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

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

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

2.3 東昆侖地震帶重現水平及置信區間的輪廓似然估計

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

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

根據上述步驟得到的重現期分別為20年、50年、100年、500年的重現水平和置信區間如圖4所示.

圖4 重現水平及置信區間的輪廓似然估計（a） 20 年重現期；（b） 50 年重現期；（c） 100 年重現期；（d） 500 年重現期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法和輪廓似然估計得到的估計結果對比分析詳見表2.估計結果從數值計算的角度佐證了對于重現水平的估計，輪廓似然法與極大似然法是等效的（Murphy，Vander Vaart，2000）.

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

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

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

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

3 討論與結論

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

基于輪廓似然估計的GEV分布模型能夠對地震危險性作出相對比較客觀的評價，但該模型所用數據為區組最大值信息，在模型構建過程中對觀測信息利用不夠充分，使得由模型得到的預測結果與實際情況存在偏差，作者后續將主要從歷史信息挖掘和模型調整與優化兩個方向入手，以構建更加有效的地震預報預測模型.