張錕 任魯川 田建偉 劉哲

中國地震局防災科技學院，河北省三河市燕郊 065201

0 引言

震級上限是指潛在震源區(qū)內可能發(fā)生的最大地震的震級，因此可以認為在潛在震源區(qū)內未來發(fā)生超過該震級地震的概率幾乎為零（胡聿賢，1999）。作為描述地震活動性的參數(shù)之一，震級上限在地震危險性分析中具有非常重要的作用（Ward，1997），因為震級上限的取值將直接影響地震危險性分析的結果。潛在震源區(qū)的震級上限主要是通過分析所選區(qū)域本身的地震活動性和地震活動構造特征來確定，目前主要有確定性估計和概率估計2種方法，前者是依據(jù)震級與斷層幾何尺度之間的經驗關系，將獲得的區(qū)域內的最大歷史地震作為震級上限來考慮（Kijko，2004）；后者是基于區(qū)域內歷史地震記錄建立震級和頻度的關系模型，將基于該模型估計所得的震級極限值作為震級上限，如基于G-R截斷模型的震級上限的估計方法（Cornell，1968），部分學者在此基礎上提出的改進方法（徐偉進等，2012）都屬于概率估計方法。這種基于G-R截斷模型的方法求得的震級上限不再是無限大，而是具有一定物理意義的極限值，但是該方法并不是利用強震目錄來研究震級上限，其結果的可靠性也稍差（任雪梅等，2012）。為解決這一問題，基于極值理論的震級上限估計方法開始得到應用，這類方法對數(shù)據(jù)精度要求不高，關注的是一段時間內的震級最大值，非常適合歷史地震記錄時間長但低震級地震記錄缺失的地區(qū)（胡聿賢，1999）。由于不同區(qū)域地震震級分布存在很大差異，因此極值類型分布在不同區(qū)域也存在不同的適用性（張衛(wèi)東等，2005）。近年，Pisarenko等提出基于廣義極值理論的震級上限的估計方法（Pisarenko et al，2008），該方法通過引入形狀參數(shù)ξ，將3類極值的漸近分布統(tǒng)一為一個分布模型，降低了模型選擇的風險，故使用極值方法分析地震危險性有著更大的適用性（陳虹，1996）。錢小仕等將該方法用于中國地震危險性分析的相關研究中（錢小仕等，2012），任晴晴等據(jù)此方法討論了中國大陸活動地塊邊界帶最大震級分布特征（任晴晴等，2013），但是沒有對區(qū)域地質背景和地震活動性特征作更深入的分析，本文補充了相關的研究并用于估計潛在地震海嘯源區(qū)震級上限。

潛在地震海嘯源區(qū)是因海域發(fā)生地震而可能觸發(fā)海嘯的潛在震源。本文采用廣義極值分布（Generalized extreme value distribution，簡稱GEV）來估計震級上限及強震重現(xiàn)水平。海嘯災害歷史記錄顯示，能觸發(fā)海嘯且導致災害的地震的震級要足夠大（MW≥6.5級）、震源足夠淺（一般為淺源地震），且震中區(qū)域海水足夠深（陳颙等，2005）。此外，海溝俯沖帶在潛在地震海嘯源位置界定中廣受關注（任魯川等，2014）。通過對琉球海溝俯沖帶地形地貌、水深、地震等相關資料的分析，認為該區(qū)域具備發(fā)生破壞性地震海嘯的條件，因此，本文選定該區(qū)域作為研究區(qū)。

1 廣義極值理論簡介

設x1，x2，…，xn是獨立同分布的隨機變量，具有共同的分布函數(shù)F（x），對自然數(shù)n，令Mn＝max｛X1，…，Xn｝，mn＝min｛X1，…，Xn｝，表示n個隨機變量的最大值和最小值，則

根據(jù) Fisher-Tippet的極值類型定理，如果存在常數(shù)列｛an＞0｝和｛bn｝，使得

成立，其中H（x）是非退化的分布函數(shù)，則H（x）必屬于下列3種函數(shù)類型，即Gumbel分布、Frechet分布、Weibull分布中的1類，這3類分布統(tǒng)稱為極值分布，引入位置參數(shù)μ、尺度參數(shù)σ及形狀參數(shù)ξ，將3種類型不同的極值分布統(tǒng)一為

式中，當ξ＝0，屬于 Gumbel極值分布；ξ＞0，屬于 Frechet極值分布；ξ＜0，屬于 Weibull極值分布（史道濟，2006；Hüsler，1984）

根據(jù)式（3）廣義極值的分布形式獲得的期望E（x）和方差Var（x）為

其中，Γ為 Gamma函數(shù)，其表達式為：

根據(jù)分位數(shù)的定義，廣義極值分布的p（0＜p＜1）處的分位數(shù)可表達為

當ξ＝0時，為Gumbel分布，其分位數(shù)為

2 基于廣義極值分布的震級上限估計

在估計T年強震重現(xiàn)水平時，依地震活動特性將時間域T年分為n段，取每段時間內的最大震級Mi，構成樣本Mmax＝（M1，M2……Mn），在進行時間域劃分和選擇震級時，要考慮 2點，一是地震的活躍周期和平靜周期，二是要剔除前震和余震。

構造示性函數(shù)I，當Mi≥X時，記為I＝1；Mi＜X時，記為I＝0。令則K表示樣本Mmax＝（M1，M2……Mn）中震級大于X的地震次數(shù)，則K符合二項分布，參數(shù)為n，p＝1－G（X），其中G（X）表示單位時間內發(fā)生1次M≤X地震的概率，這里的G（X）就是式（3）中的廣義極值分布公式。

2.1 平均復發(fā)周期

設最大震級超過X一次的地震平均復發(fā)周期為T（X），即T（X）年內最大震級超過X的平均次數(shù)為 1，即期望值E（k）＝np＝1，則T（X）［1－G（X）］＝1，故有

2.2 T年內的強震重現(xiàn)水平

式（9）中，X表示T年內的最大震級（錢小仕等，2012），當T→＋∞時，X即為所求潛在地震海嘯源區(qū)的震級上限。

2.3 震級上限的極大似然估計

由式（3）得到廣義極值分布的密度函數(shù)，求得其對數(shù)似然函數(shù)為

當 0＜p＜1時，根據(jù)式（9）可以得到參數(shù)σ，μ，ξ的最大對數(shù)似然估計值代入式（6）和（7），獲得分位數(shù)的極大似然估計

當ξ＜0時，根據(jù)式（3）可知，xp存在上限，從而獲得震級上限的極大似然估計值為

3 震級上限的置信區(qū)間估計

4 潛在地震海嘯源區(qū)——琉球海溝俯沖帶

琉球海溝俯沖帶屬于溝-弧-盆體系，其中的“溝”指的是琉球海溝，它是溝-盆-系與大洋盆地的天然分界，在地貌上表現(xiàn)為島坡坡麓的深溝，是1條向SE凸出，向WN方向傾沒的弧形海溝，呈環(huán)帶狀環(huán)繞琉球島弧延伸，北端以九州-帛琉海嶺為界，南端位于臺灣島中部外海，總長約1350km，平均寬度60km，平均水深大于6km，最大水深7.781km位于123°E附近。由于加瓜海脊向北俯沖到琉球島弧以下，琉球海溝在此處發(fā)生較大的變形，導致東西兩區(qū)呈不同的地球物理場特征，海底沉積物結構和物質組成也不同。在123°E以西，海溝寬度和深度逐漸減小，在地形上演變成海底峽谷及深海盆地；在123°E以東，海底地形平坦，呈倒“V”字型；在日本島的中部以南，帛琉-九州海脊的北段消失（張訓華，2008；劉宗臣等，2005）。琉球海溝的西坡是具有大陸性質琉球島?。ǜ呦榱?，2003），是由琉球諸島形成的島鏈，稱之為琉球島弧，北起九州島南端，南至臺灣島東部，全長1200km，屬于雙列島弧，內弧主要是古琉球火山帶，是一條水下火山脊，外弧是琉球群島的主體（張訓華，2008），其地貌分布格局主要由琉球群島隆褶帶、弧前盆地和八重山海脊帶控制（劉宗臣等，2005）。琉球島弧的西側是沖繩海槽，是正在發(fā)育的弧后裂谷盆地，目前仍然是大陸地殼（周祖翼等，2001；王述功等，1998）。

前人研究成果表明，菲律賓板塊正在以55°左右的傾角向琉球群島和東海之下俯沖，這導致琉球島弧-海溝系之下存在一個明顯的貝尼奧夫地震帶，該地震帶以25°～75°的角度向西北傾斜，插入地殼下達280km，屬于環(huán)太平洋地震帶的一部分，是至今仍在活動的強地震帶，地震活動頻度高，其地震震中分布密度遠遠高于板塊邊界的大陸內部，曾發(fā)生多次6級以上地震。（李乃勝，2000）（圖1）。因此該區(qū)域具備破壞性海嘯地震的誘發(fā)條件，且一旦發(fā)生海嘯有可能對中國東南沿海地區(qū)造成影響，所以筆者選定22°～32.5°N，120.5°～133°E作為研究區(qū)域。

圖1 琉球海溝俯沖帶地震分布（1900～2010年）（M W≥7.0）

根據(jù)USGS提供的地震目錄分別繪制了1900～2010年MW≥6.5（圖2）和1950～2010年的MW≥5.0地震的M-t圖（圖3）。通過分析1900～2010年的地震目錄，結合2張M-t圖可以看出，琉球海溝俯沖帶的地震活動周期為10年左右，為保證所選的時間步長與地震活動的平靜期和活躍期盡可能吻合以及所選時間步長下數(shù)據(jù)的完整性，確定時間步長t＝10年，并選取1910年為初始時間（表1），然后選取出每個時間步長內的最大震級，考慮到不同的初始時間可能會對所選震級的序列產生影響，因此選擇1906和1908年作為初始時間，選出新的最大震級序列（表1），從表1可以看出，初始時間的變化導致震級序列發(fā)生了較小的變化。通過選取初始時間1908年的震級數(shù)據(jù)進行計算并與初始時間為1910年的震級上限的計算結果進行比較發(fā)現(xiàn)變化不大，說明計算結果比較穩(wěn)定可靠。

圖2 1900～2010球海溝俯沖帶 M-t圖（M W≥6.5）

圖3 1950～2010球海溝俯沖帶M-t圖（M W≥5）

表1 琉球海溝沉降帶地震的最大震級數(shù)據(jù)

5 數(shù)據(jù)的分析和處理

根據(jù)式（9）進行最大似然估計，得到式（3）中的參數(shù)ξ、μ、σ，其值依次為－0.4162811，7.5876929，0.3366287。由于ξ＜0，屬于廣義極值第3類分布，因此具有上限值，可將獲得的參數(shù)代入式（12），獲得廣義極值模型估計值的上端點為8.4，即該區(qū)域的震級上限為8.4級，根據(jù)式（14）得置信水平為95%的置信區(qū)間為 8.4±0.3708。

對所得的估計結果進行擬合診斷，得到的結果如圖4。圖4（a）是概率P-P圖，橫坐標表示實際數(shù)據(jù)的累積概率，縱坐標表示選用極值模型的累計概率；圖4（b）是分位數(shù)Q-Q圖，橫坐標表示的是所選分布模型的分位數(shù)，縱坐標表示的是實際數(shù)據(jù)的分位數(shù)。從圖4（a）和圖4（b）可以看出，所有的點沿直線分布，未顯示一定的偏離性，因此不能否定所選用的分布模型。圖4（c）是重現(xiàn)水平圖，因為當ξ＝0時，分布模型為 Gumbel分布，其重現(xiàn)水平圖為1條直線，通過計算知，所得ξ趨近于0，因此其重現(xiàn)水平近似為1條直線，由于ξ＜0重現(xiàn)水平圖是1條凸曲線，位于中間的曲線表示的是強震重現(xiàn)水平，上、下2條曲線代表的是考慮置信水平為95%時的置信區(qū)間。圖4（d）是密度曲線和直方圖，從圖中可以看出概率密度曲線的估計和直方圖相吻合。圖4表明可以選用GEV模型進行擬合。

將所獲得的參數(shù)代入式（9），分別獲得自2010年起，研究區(qū)域在未來100年內的最大地震的重現(xiàn)水平為MW＝8.1，在未來50年內的最大地震的重現(xiàn)水平為MW＝8.0，在未來30年內最大地震的重現(xiàn)水平為MW＝7.8。通過廣義極值計算獲得式（13）中的協(xié)方差矩陣

根據(jù)公式（13）分別得出t＝30年、t＝50年、t＝100年的置信水平為95%的置信區(qū)間（表2）。通過對形狀參數(shù)ξ的置信區(qū)間進行計算，獲得置信水平為95%的置信區(qū)間為（－0.6508，－0.1738），根據(jù)廣義極值理論，ξ值小于0時為 Weibull分布，因此琉球海溝沉降帶震級上限符合Weibull分布。

圖4 廣義極值分布的擬合診斷

表2 t年強震重現(xiàn)水平的置信區(qū)間

6 結語

（1）采用廣義極值模型估計潛在地震海嘯源區(qū)震級上限，在選取時間步長時考慮了地震活動平靜期和活躍期，使所獲得的強震樣本與地震的活動特征相吻合，在選取潛在地震海嘯源區(qū)時充分考慮地質構造背景，因此本文的參數(shù)估計是在充分考慮所選區(qū)域的地震活動性特征和地質構造背景的前提下采用廣義極值模型進行的。

（2）采用廣義極值模型估計琉球海溝俯沖帶震級上限值為8.4級，對應的置信水平為95%的置信區(qū)間為 8.4±0.3708；估計琉球海溝俯沖帶在未來100年、50年、30年內的最大地震的重現(xiàn)水平分別為 8.1級、8.0級、7.8級，對應的置信水平為 95%的置信區(qū)間為 8.1±0.4208、8.0±0.4822、7.8±0.5417。

（3）所獲得的形狀參數(shù)ξ的估計值為－0.4162811，其對應的置信度為95%的置信區(qū)間為（－0.6508，-0.1738），根據(jù)廣義極值理論ξ＜0時，其分布屬于Weibull分布，因此可以推斷，琉球海溝俯沖帶強震震級符合Weibull分布，同時也說明所獲得的強震震級是有上限的，這與地震活動的特征相吻合。