基于廣義極值理論的潛在地震海嘯源震級(jí)上限及強(qiáng)震重現(xiàn)水平的估計(jì)
——以琉球海溝俯沖帶為例

張錕 任魯川 田建偉 劉哲

中國地震局防災(zāi)科技學(xué)院，河北省三河市燕郊 065201

0 引言

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

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

1 廣義極值理論簡(jiǎn)介

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

根據(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)一為

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

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

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

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

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

2 基于廣義極值分布的震級(jí)上限估計(jì)

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

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

2.1 平均復(fù)發(fā)周期

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

2.2 T年內(nèi)的強(qiáng)震重現(xiàn)水平

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

2.3 震級(jí)上限的極大似然估計(jì)

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

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

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

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

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

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

前人研究成果表明，菲律賓板塊正在以55°左右的傾角向琉球群島和東海之下俯沖，這導(dǎo)致琉球島弧-海溝系之下存在一個(gè)明顯的貝尼奧夫地震帶，該地震帶以25°～75°的角度向西北傾斜，插入地殼下達(dá)280km，屬于環(huán)太平洋地震帶的一部分，是至今仍在活動(dòng)的強(qiáng)地震帶，地震活動(dòng)頻度高，其地震震中分布密度遠(yuǎn)遠(yuǎn)高于板塊邊界的大陸內(nèi)部，曾發(fā)生多次6級(jí)以上地震。（李乃勝，2000）（圖1）。因此該區(qū)域具備破壞性海嘯地震的誘發(fā)條件，且一旦發(fā)生海嘯有可能對(duì)中國東南沿海地區(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年的地震目錄，結(jié)合2張M-t圖可以看出，琉球海溝俯沖帶的地震活動(dòng)周期為10年左右，為保證所選的時(shí)間步長與地震活動(dòng)的平靜期和活躍期盡可能吻合以及所選時(shí)間步長下數(shù)據(jù)的完整性，確定時(shí)間步長t＝10年，并選取1910年為初始時(shí)間（表1），然后選取出每個(gè)時(shí)間步長內(nèi)的最大震級(jí)，考慮到不同的初始時(shí)間可能會(huì)對(duì)所選震級(jí)的序列產(chǎn)生影響，因此選擇1906和1908年作為初始時(shí)間，選出新的最大震級(jí)序列（表1），從表1可以看出，初始時(shí)間的變化導(dǎo)致震級(jí)序列發(fā)生了較小的變化。通過選取初始時(shí)間1908年的震級(jí)數(shù)據(jù)進(jìn)行計(jì)算并與初始時(shí)間為1910年的震級(jí)上限的計(jì)算結(jié)果進(jìn)行比較發(fā)現(xiàn)變化不大，說明計(jì)算結(jié)果比較穩(wěn)定可靠。

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

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

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

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

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

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

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

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

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

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

6 結(jié)語

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

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

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