應(yīng)用統(tǒng)計(jì)物理方法分析川滇地區(qū)地震危險(xiǎn)性

王 凡 沈正康 王 敏 王閻昭 陶 瑋

1)中國地震局地質(zhì)研究所，地震動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京 100029

2)北京大學(xué)地球與空間科學(xué)學(xué)院，北京 100871

0 引言

川滇地區(qū)位于青藏高原的東緣，構(gòu)造運(yùn)動(dòng)強(qiáng)烈。受印度板塊NNE向推擠和青藏高原E向擠出的動(dòng)力學(xué)作用控制，川滇地區(qū)構(gòu)造形變場主要表現(xiàn)為繞東喜馬拉雅造山結(jié)的順時(shí)針大回轉(zhuǎn)及與之相伴的大小地塊間的斷層錯(cuò)動(dòng)，其主要活動(dòng)斷裂包括構(gòu)成青藏高原與四川盆地邊界的龍門山斷裂帶，高原內(nèi)部的龍日壩和岷江斷裂帶，構(gòu)成川滇“菱形地塊”東邊界與北邊界的鮮水河－小江斷裂帶及南邊界的紅河、南華－楚雄－建水、無量山斷裂帶，菱形地塊內(nèi)部的麗江－小金河斷裂帶，南部的龍陵－瑞麗斷裂帶，以及西南的怒江、瀾滄江斷裂帶等。伴隨著強(qiáng)烈的構(gòu)造運(yùn)動(dòng)，本地區(qū)地震活動(dòng)頻繁，據(jù)中國地震臺(tái)網(wǎng)中心記錄(圖1)，自1970年以來該地區(qū)發(fā)生5.5級(jí)以上地震達(dá)280次之多。

圖1 川滇地區(qū)活動(dòng)構(gòu)造及1500年以來5.0級(jí)以上地震Fig.1 Tectonic map of Sichuan-Yunnan region and M≥5.0 earthquakes of 1500—2010.

根據(jù)歷史記載，大地震通常會(huì)造成巨大的人員傷亡和財(cái)產(chǎn)損失。因此，川滇地區(qū)的地震危險(xiǎn)性分析一直是地震預(yù)測研究的熱點(diǎn)(聞學(xué)澤，2001;易桂喜等，2002，2004，2006，2008;蔣長勝等，2008，2010)。易桂喜等(2002)和徐錫偉等(2005)利用簡單的點(diǎn)過程統(tǒng)計(jì)方法研究地震復(fù)發(fā)行為，結(jié)果表明，這一區(qū)域活動(dòng)斷裂帶的強(qiáng)震復(fù)發(fā)表現(xiàn)出趨于隨機(jī)的叢集行為，其過程不存在良好的準(zhǔn)周期性;同時(shí)，基于地震復(fù)發(fā)間隔與震級(jí)的統(tǒng)計(jì)研究表明，該地區(qū)所發(fā)生的地震也不具有“強(qiáng)度可預(yù)報(bào)”或“時(shí)間可預(yù)報(bào)”的性質(zhì)。因此，基于斷裂帶整體強(qiáng)震復(fù)發(fā)性質(zhì)的中長期地震危險(xiǎn)性評(píng)估具有一定困難。

自20世紀(jì)80年代末，國際上逐漸將現(xiàn)代統(tǒng)計(jì)理論和方法應(yīng)用于地震預(yù)測問題的研究并取得多方面的進(jìn)展，已成為物理學(xué)和地震學(xué)之間一個(gè)十分活躍的交叉領(lǐng)域。Kagan等(1994)提出了基于地震空間相關(guān)性的方法，根據(jù)歷史地震目錄預(yù)測未來地震發(fā)生的概率，并將此方法應(yīng)用于日本(Jackson et al.，1999;Kagan et al.，2000)、中國(Rong et al.，2002)以及全球(Kagan et al.，2011)。Rundle等將圖像信息學(xué)(Pattern Informatics)算法應(yīng)用于地震活動(dòng)性圖像分析，并嘗試對(duì)加州的地震做出預(yù)測(Rundle et al.，2000，2003;Tiampo et al.，2002)。蔣長勝(2008)將此算法應(yīng)用于川滇地區(qū)，其回溯性檢驗(yàn)表明該方法對(duì)于川滇地區(qū)的中長期地震活動(dòng)具有一定的預(yù)測能力。Keilis-Borok等(1990)利用TIPs(Time of Increased Probability)方法對(duì)意大利中部地區(qū)進(jìn)行了長期強(qiáng)震預(yù)測。Huang等(2001)利用RTL(Regional-Time-Length)方法研究了M 7.2阪神地震造成的區(qū)域地震危險(xiǎn)性變化。除利用地震目錄進(jìn)行統(tǒng)計(jì)分析的方法外，Shen等(2007)提出了基于地殼形變場建立地震危險(xiǎn)性預(yù)測模型的方法，研究發(fā)現(xiàn)，過去50a以來南加州地區(qū)的強(qiáng)震分布與應(yīng)變率分布有很高的相關(guān)性。以上這些研究證實(shí)了基于統(tǒng)計(jì)理論的方法應(yīng)用于不同觀測數(shù)據(jù)進(jìn)行地震危險(xiǎn)性分析的有效性。

本文分別嘗試應(yīng)用基于地震空間相關(guān)性和地殼形變場的方法分析川滇地區(qū)的長期地震危險(xiǎn)性，并通過對(duì)預(yù)測結(jié)果的回溯性檢驗(yàn)分析該方法的有效性。

1 基于地震空間相關(guān)性的方法

1.1 基于地震空間相關(guān)性的地震危險(xiǎn)性模型

基于地震空間相關(guān)性的地震危險(xiǎn)性預(yù)測方法由Kagan等(1994)最先提出，其基本思想是根據(jù)歷史地震目錄建立統(tǒng)計(jì)學(xué)模型，從而得到未來發(fā)生地震的概率。假定地震發(fā)生概率與歷史地震發(fā)生頻度成正比，地震震級(jí)－數(shù)目分布遵從修正的古登堡－里克特關(guān)系(Bird et al.，2004)，地震發(fā)生為泊松過程，由此可以根據(jù)區(qū)域地震目錄給出強(qiáng)震發(fā)生的統(tǒng)計(jì)預(yù)測。地震概率密度函數(shù)可表示為

其中，Λ為單位時(shí)間、單位面積內(nèi)發(fā)生地震的概率密度，θ和φ分別表示經(jīng)度和緯度，m是震級(jí)，G為地震震級(jí)-數(shù)目分布函數(shù)，F(xiàn)為地震空間分布概率函數(shù)，是由根據(jù)歷史地震目錄得到的地震空間分布概率函數(shù)與背景場地震概率C的加權(quán)和。之所以采用加權(quán)和是為了避免某些沒有歷史地震的地區(qū)未來發(fā)生地震的概率為零。地震空間分布概率函數(shù)是歷史地震目錄的平滑函數(shù)，本文采用的平滑函數(shù)如下:

其中φi和θi分別是地震i震中的緯度和經(jīng)度，rs是尺度參數(shù)，它與很多因素有關(guān)，如地質(zhì)構(gòu)造、地殼圈層厚度以及地震定位的準(zhǔn)確度等，因此，不同區(qū)域的rs是不同的。

利用1500—2010年5.0級(jí)以上地震目錄提供約束，給出了川滇地區(qū)未來發(fā)生5.0級(jí)以上地震的概率(圖2)。在基于地震空間相關(guān)性的預(yù)測方法中，每一次地震都是獨(dú)立的，結(jié)果大震的前震和余震都作為獨(dú)立事件參與統(tǒng)計(jì)。因此，只選取5.0級(jí)以上的地震，既能保證有充分的數(shù)據(jù)約束預(yù)測模型，又可以在一定程度上去除部分前震以及余震以獲得盡可能合理的結(jié)果。此外，我們還在一定的范圍內(nèi)遍歷尺度參數(shù)rs，結(jié)果表明，當(dāng)rs=15時(shí)預(yù)測模型能更好地符合檢驗(yàn)階段的地震目錄。由圖2可以看出，地震危險(xiǎn)性較高的地區(qū)包括龍門山斷裂中2008年汶川地震破裂段及鄰域、虎牙斷裂、鮮水河斷裂、龍陵地區(qū)、馬邊斷裂及麗江－小金河斷裂中－南段。此外，地震危險(xiǎn)性較高地區(qū)還包括小江斷裂、楚雄－建水?dāng)嗔涯隙?、無量山斷裂中段及鄰域地區(qū)等。

圖2 由1500—2010年地震目錄約束得到的5級(jí)以上地震發(fā)生的概率，背景色表示地震概率Fig.2 Earthquake long-term potential based on earthquake catalog from 1500 to 2010.Colors show the long-term probability of earthquake occurrence.

1.2 模型結(jié)果的回溯檢驗(yàn)

基于地震空間相關(guān)性的預(yù)測方法能否準(zhǔn)確估計(jì)川滇地區(qū)長期地震危險(xiǎn)性?為驗(yàn)證其有效性，下面采用2種方法對(duì)模型的預(yù)測結(jié)果進(jìn)行了回溯性檢驗(yàn)。方法一為基于似然函數(shù)的檢驗(yàn)方法，即用似然函數(shù)比來檢驗(yàn)預(yù)測模型與實(shí)際發(fā)生地震的符合程度(Kagan et al.，1977)。其中，j是網(wǎng)格序號(hào)，N是網(wǎng)格數(shù)，vj是歸一化的網(wǎng)格內(nèi)發(fā)生地震的概率，τj為歸一化的網(wǎng)格面積，n為檢驗(yàn)階段地震目錄中的地震數(shù)，λi為預(yù)測模型得到的地震i的震中所在網(wǎng)格的單位面積單位時(shí)間內(nèi)發(fā)生地震的概率，ξ為單位面積、單位時(shí)間內(nèi)發(fā)生地震的統(tǒng)計(jì)值，I0為基于預(yù)測模型的地震目錄，I1為基于檢驗(yàn)階段的地震目錄，兩者之差反映了預(yù)測模型的準(zhǔn)確性(詳細(xì)介紹見Kagan，2009)。

檢驗(yàn)方法二為區(qū)域地震概率與地震發(fā)生數(shù)的累積函數(shù)檢驗(yàn)(Kagan，2009)。將研究區(qū)分成小的地域單元，按照地域內(nèi)地震發(fā)生概率由高到低排列，排列后的概率函數(shù)按面積積分，同時(shí)將單元內(nèi)發(fā)生的地震數(shù)累積。將2個(gè)函數(shù)在終點(diǎn)歸一化，得到概率積分函數(shù)與地震數(shù)累積函數(shù)。如果地震發(fā)生率符合地震概率預(yù)測模型，則2條函數(shù)曲線應(yīng)有很好的相似性;如果實(shí)際地震發(fā)生曲線在概率曲線之上，表明地震的發(fā)生概率高于概率模型預(yù)期，即在概率模型預(yù)期地震發(fā)生率高的地區(qū)實(shí)際地震發(fā)生數(shù)大于概率預(yù)期值，而在概率模型預(yù)期值低的地方則相反。

圖3 利用歷史目錄的地震預(yù)測模型Fig.3 Earthquake forecast test using historical catalog data.

將1500年以來5.0級(jí)以上地震分為1500—1999年和2000—2010年2個(gè)部分，第1部分用于約束地震概率估計(jì)，所得結(jié)果見圖3a;第2部分用于回溯檢驗(yàn)。第1種檢驗(yàn)方法的結(jié)果列于表1，結(jié)果中I0、I1差異較大，表明預(yù)測模型不夠準(zhǔn)確。第2種檢驗(yàn)方法的結(jié)果見圖3b，顯示實(shí)際地震發(fā)生曲線位于概率模型曲線之下，同樣表明概率模型(圖3a)的預(yù)測結(jié)果并不好?？疾爝^去10a地震發(fā)生的空間分布，發(fā)現(xiàn)確實(shí)有不少地震發(fā)生在預(yù)測概率相對(duì)高的地區(qū)，但仍有相當(dāng)多的地震發(fā)生在預(yù)測概率較低的地區(qū)，其中特別引人注目的是2008年的汶川地震。綜其原因可能是預(yù)測方法只是在大樣本空間統(tǒng)計(jì)意義上有效，當(dāng)統(tǒng)計(jì)樣本數(shù)較少時(shí)，如本例，小概率事件出現(xiàn)的可能性大大增加，因此對(duì)統(tǒng)計(jì)結(jié)果的影響極大。普遍認(rèn)為發(fā)生汶川地震的龍門山斷裂帶的特征地震復(fù)發(fā)周期在千a以上(Densmore et al.，2007;Burchfiel et al.，2008;Shen et al.，2009;Zhang et al.，2010)，有歷史記錄的500a以來極少有強(qiáng)震發(fā)生，故本方法得到的汶川地震發(fā)生地預(yù)測概率值很低，預(yù)測并不成功。

表1 基于似然函數(shù)的檢驗(yàn)結(jié)果Table 1 Forecast testing catalog and information scores

針對(duì)上述預(yù)測模型存在的樣本數(shù)較少的問題，進(jìn)行了另一個(gè)回溯檢驗(yàn)，將1500—1972年的地震目錄用于約束地震概率估計(jì)，得到的預(yù)測概率如圖4a所示，1973—2010年的地震目錄用于模型檢驗(yàn)。檢驗(yàn)方法一的結(jié)果見表一，I0和I1二者差異約為12%，表明模型預(yù)測效果較好。檢驗(yàn)方法二得到的結(jié)果如圖4b所示，從中可以看出預(yù)測函數(shù)曲線與實(shí)際地震發(fā)生函數(shù)曲線符合得較好，雖然二者在10%高概率地區(qū)的范圍仍有一定差距，但在20%高概率地區(qū)的范圍已相當(dāng)接近。結(jié)果進(jìn)一步顯示50%強(qiáng)震發(fā)生在模型給定20%高概率地區(qū)，而80%強(qiáng)震發(fā)生在50%高概率地區(qū)。2種檢驗(yàn)都表明在較為充分的樣本空間內(nèi)基于地震空間相關(guān)性的地震預(yù)測方法可以很好地反映地震發(fā)生率，雖然在中、短期地震預(yù)測中有一定局限性，但對(duì)于長期地震預(yù)測不失為一種有效方法。

圖4 a背景為利用1500—1972年地震目錄約束得到的5.0級(jí)以上地震發(fā)生的概率，黑色圓圈為1973—2010年發(fā)生的5.0級(jí)以上的地震;b區(qū)域地震概率與地震發(fā)生數(shù)的累積函數(shù)檢驗(yàn)，紅線為地震預(yù)測概率積分函數(shù)，藍(lán)線為實(shí)際發(fā)生地震的累計(jì)數(shù)Fig.4 The same as Fig.3 except that the earthquake catalog data are of 1500－1972 and the M ＞5.0 events are of 1973 －2010 in(a)，and the test curves in(b)are plotted accordingly.

2 基于地殼形變場的方法

2.1 基于地殼形變場的地震危險(xiǎn)性模型

基于地殼形變場的地震危險(xiǎn)性模型(Shen et al.，2007)假定區(qū)域長期地震危險(xiǎn)性與區(qū)域構(gòu)造應(yīng)變率成正比，地震發(fā)生數(shù)遵從修正的古登堡－里克特震級(jí)－地震數(shù)統(tǒng)計(jì)關(guān)系，則地震發(fā)生概率可表示為

其中A(x)為大地測量得到的地表應(yīng)變率場，G(m)為修正的古登堡－里克特震級(jí)－地震數(shù)統(tǒng)計(jì)關(guān)系。本文收集了研究區(qū)內(nèi)1999—2011年所有可利用的GPS觀測資料，采用非線性模型模擬臺(tái)站位置變化，盡可能地消除2001年MW7.8可可西里地震、2004年MW9.2蘇門答臘地震以及2008年MW7.9汶川地震的同震位移和震后形變的影響，得到了反映區(qū)域構(gòu)造運(yùn)動(dòng)的速度場(圖5)，具體方法見王敏等(2003)和Shen等(2011)。

由區(qū)域構(gòu)造運(yùn)動(dòng)速度場求解應(yīng)變率場，采用了一種連續(xù)滑移回歸方法(Shen et al.，1996;Jackson et al.，1997)。其主要原理是對(duì)于每個(gè)計(jì)算點(diǎn)，集合周邊臺(tái)站速度場數(shù)據(jù)，用最小二乘法估算平均應(yīng)變率場，所用速度場數(shù)據(jù)依據(jù)相鄰臺(tái)站的方位角大小和與計(jì)算點(diǎn)的距離加權(quán)，其中距離加權(quán)函數(shù)表示為exp(－Δ2/σ2)的形式，Δ為臺(tái)站至計(jì)算點(diǎn)距離，σ為平滑因子。不同σ所對(duì)應(yīng)的總權(quán)重和應(yīng)變率估值的形式誤差之間存在負(fù)相關(guān)關(guān)系，這一負(fù)相關(guān)關(guān)系中的最佳平衡點(diǎn)即為最優(yōu)的σ值。

圖5 由GPS觀測獲得的區(qū)域水平運(yùn)動(dòng)速度場(相對(duì)穩(wěn)定的歐亞板塊)Fig.5 GPS horizontal velocity field with respect to the stable Eurasian reference frame.

圖6 背景為應(yīng)變率第二不變量，白色菱形為1500—1972年間M＞6.0地震的震中，白色圓圈為1973年以來M＞5.5地震的震中Fig.6 Geodetic strain rate field and spatial distribution of strong earthquakes.The color background is the second invariant of strain rates，and white diamonds and white circles are epicenters of M ＞ 6.0 events occurring in 1500－1972 and M ＞5.5 events occurring in 1973－2010，respectively.

2.2 模型結(jié)果的回溯檢驗(yàn)

為了對(duì)地震的空間分布與應(yīng)變率場的相關(guān)性有更客觀的認(rèn)識(shí)，應(yīng)用1.2節(jié)中區(qū)域地震概率與地震發(fā)生數(shù)的累積函數(shù)檢驗(yàn)方法對(duì)上述預(yù)測模型進(jìn)行回溯性檢驗(yàn)。圖7a和圖7b分別為1973年以來5.5級(jí)以上地震和1500年以來6.0級(jí)以上地震與應(yīng)變率場的關(guān)系。在有地震儀器記錄之前(即1500—1973年)地震發(fā)生地點(diǎn)與震級(jí)量度都存在一定不確定性，并且可能還有未被記錄的地震?？紤]川滇地區(qū)歷史上居住人口較為稠密，過去500a中6.0級(jí)以上地震被遺漏的可能性不大，而導(dǎo)出應(yīng)變率場的平滑因子也在幾十km范圍，因此地點(diǎn)不確定性的影響不會(huì)很大。從圖7a可以看出，1973年以來發(fā)生在低應(yīng)變率地區(qū)的5.5級(jí)以上地震的數(shù)目明顯多于基于應(yīng)變率場的預(yù)期，說明應(yīng)變率場不能作為過去37年強(qiáng)震發(fā)生地點(diǎn)的指標(biāo)，而圖7b則顯示1500年以來6.0級(jí)以上地震的空間分布與應(yīng)變率場存在很好的相關(guān)性，即在數(shù)百a時(shí)間尺度上，應(yīng)變率的高低可以作為強(qiáng)震發(fā)生地點(diǎn)的指示因子。該地區(qū)10a尺度的地殼應(yīng)變率場實(shí)際上反映了本地區(qū)長期的構(gòu)造形變場，而這一形變場的主體為彈性形變，只能通過地震、特別是強(qiáng)震來釋放，但其釋放時(shí)間需要數(shù)百a。

圖7 應(yīng)變率與地震發(fā)生數(shù)的相關(guān)關(guān)系Fig.7 Geodetic strain rate vs.earthquake count.

3 討論與結(jié)論

Rong等(2002)應(yīng)用基于地震空間相關(guān)性的預(yù)測方法研究了整個(gè)中國大陸的地震危險(xiǎn)性，采用的地震目錄為1978—1999年，截止震級(jí)為M＞5.4，檢驗(yàn)時(shí)間區(qū)間僅為2a(2000—2002年)。相比之下，我們研究的川滇地區(qū)范圍小且地震頻度大，約束模型和檢驗(yàn)?zāi)Ｐ偷臄?shù)據(jù)區(qū)間都遠(yuǎn)遠(yuǎn)大于前者。因?yàn)槭褂玫膮?shù)不同，2個(gè)模型的預(yù)測結(jié)果不能嚴(yán)格比較，但從回溯檢驗(yàn)的結(jié)果不難看出(本文圖4b與Rong等(2002)圖3)，本文模型預(yù)測結(jié)果與實(shí)際發(fā)生地震情況符合程度明顯更好。這也再一次說明基于地震空間相關(guān)性的預(yù)測方法對(duì)歷史地震目錄有很強(qiáng)的依賴性，歷史地震目錄越長，記錄越完整，模型的預(yù)測結(jié)果越可靠。在較為充分的樣本空間內(nèi)，基于地震空間相關(guān)性的預(yù)測方法對(duì)于長期地震預(yù)測具有一定的有效性。

地震發(fā)生的過程很短，也就是幾十s，但是地震孕育的過程很長，通常要幾百a，甚至上千a。在幾十a(chǎn)時(shí)間尺度內(nèi)，地震的發(fā)生有叢集性的特征，而基于地震空間相關(guān)性的方法所根據(jù)的地震目錄只有近幾十a(chǎn)比較完整，所以對(duì)這一時(shí)段有較高的有效性;而數(shù)十a(chǎn)的GPS觀測獲得的區(qū)域應(yīng)變率場更能代表區(qū)域長期的應(yīng)變率場，在這個(gè)時(shí)間尺度內(nèi)，地震的發(fā)生主要表現(xiàn)為泊松過程。不難理解，相比過去30a發(fā)生的地震，過去500a間發(fā)生的地震與基于地殼形變場的地震預(yù)測結(jié)果有更好的一致性。因此，區(qū)域應(yīng)變率場可以作為長期地震危險(xiǎn)性預(yù)測的指標(biāo)，但需要特別關(guān)注2個(gè)方面:其一，由數(shù)十a(chǎn)尺度的大地測量資料獲得的區(qū)域應(yīng)變率場應(yīng)該能夠代表區(qū)域長期的應(yīng)變率場;其二，大地測量觀測得到的應(yīng)變通常包含了彈性應(yīng)變和非彈性應(yīng)變，要注意區(qū)分這2種應(yīng)變，因?yàn)橹挥袕椥詰?yīng)變能是與地震相關(guān)的。

本文的研究表明，基于地震空間相關(guān)性和地殼形變場的地震預(yù)測方法對(duì)于川滇地區(qū)的長期地震活動(dòng)具有一定的預(yù)測能力，這為地震危險(xiǎn)性分析和地震預(yù)測研究提供了新的線索。但無論如何，本文僅是概率性的地震預(yù)測方法在川滇地區(qū)的初步應(yīng)用探索，未來應(yīng)當(dāng)考慮開展更深入的研究。