張 帆, 韓曉明, 李 娟, 裴東洋

(內(nèi)蒙古自治區(qū)地震局, 內(nèi)蒙古 呼和浩特 010051)

0 引言

地震是群災(zāi)之首,地震災(zāi)害造成的死亡人數(shù)占自然災(zāi)害死亡人數(shù)總數(shù)的54%,地震在各種自然災(zāi)害中造成的傷亡和損失最大[1]。對(duì)未來(lái)可能發(fā)生強(qiáng)震的時(shí)間和位置判斷,是地震預(yù)測(cè)工作的主要目標(biāo)。經(jīng)過多年的科學(xué)研究和實(shí)踐探索,我國(guó)形成了“長(zhǎng)-中-短-臨”的地震預(yù)測(cè)工作思路。精準(zhǔn)的短臨預(yù)報(bào)仍然是科學(xué)難題,但基于現(xiàn)有技術(shù)的地震預(yù)測(cè)工作體系已經(jīng)在防震減災(zāi)中有一定成果。地震預(yù)報(bào)工作目前處于由經(jīng)驗(yàn)預(yù)報(bào)向物理預(yù)報(bào)發(fā)展的階段,但其實(shí)際工作仍然以經(jīng)驗(yàn)預(yù)報(bào)為主[2]。地震預(yù)測(cè)指標(biāo)體系建設(shè)是經(jīng)驗(yàn)預(yù)報(bào)階段最有效的途徑,在不同區(qū)域針對(duì)不同學(xué)科建立的地震預(yù)報(bào)指標(biāo)體系在地震預(yù)測(cè)工作中發(fā)揮了重要作用[3],可應(yīng)用于不同空間和時(shí)間尺度地震風(fēng)險(xiǎn)的估計(jì)。

地震預(yù)測(cè)指標(biāo)是根據(jù)觀測(cè)資料經(jīng)過數(shù)據(jù)處理后提取的具有物理意義的參數(shù),在地震發(fā)生前,預(yù)測(cè)指標(biāo)會(huì)發(fā)生某種程度的異常變化[5]。預(yù)測(cè)指標(biāo)的建立,包括指標(biāo)的選取、觸發(fā)條件的確定、預(yù)測(cè)規(guī)則的確定等。在經(jīng)驗(yàn)預(yù)報(bào)工作中,地震預(yù)測(cè)指標(biāo)主要是通過震例回溯和統(tǒng)計(jì)檢驗(yàn)來(lái)確定[5-8]。地震頻次反映了地震活動(dòng)水平和應(yīng)力水平的變化,被作為常用的一種地震預(yù)測(cè)指標(biāo)[7-9]。地震的孕育和發(fā)生是地殼應(yīng)力、應(yīng)變的積累和釋放過程,未來(lái)強(qiáng)震震源區(qū)域周邊一定范圍內(nèi)中等地震活動(dòng)的增強(qiáng)可能蘊(yùn)含著有關(guān)地震成核、孕育和發(fā)展的信息,是震源體在物理場(chǎng)上的體現(xiàn)。Mogi[10]認(rèn)為大地震前一段時(shí)間內(nèi),震源區(qū)周圍可能出現(xiàn)環(huán)形的增強(qiáng)活動(dòng)區(qū)域,其增強(qiáng)區(qū)的半徑比震源區(qū)的半徑大2～3倍。梅世蓉等[11]對(duì)中國(guó)6級(jí)以上地震的研究也表明構(gòu)造應(yīng)力的增加,并達(dá)到臨界狀態(tài)時(shí)會(huì)出現(xiàn)地震增強(qiáng)活躍現(xiàn)象。

晉冀蒙交界地區(qū)位于燕山地塊、華北地塊和鄂爾多斯地塊交界區(qū)域,在鄂爾多斯地塊逆時(shí)針旋轉(zhuǎn)運(yùn)動(dòng)作用下構(gòu)造活動(dòng)強(qiáng)烈,形成盆嶺構(gòu)造,歷史上發(fā)生過4次7級(jí)以上地震和多次6級(jí)以上地震,但部分?jǐn)嗔讶鄙?級(jí)以上地震,形成了破裂空段[12]。2019年以來(lái),晉冀蒙交界地區(qū)出現(xiàn)ML3.0以上地震顯著增強(qiáng)現(xiàn)象,這種平靜背景下的活躍增強(qiáng)對(duì)未來(lái)地震形勢(shì)的判定是重要的。本文擬將該增強(qiáng)區(qū)域作為研究對(duì)象,考察將地震的滑動(dòng)年頻次作為預(yù)測(cè)指標(biāo)的有效性,通過震例回溯和統(tǒng)計(jì)檢驗(yàn),對(duì)指標(biāo)的預(yù)測(cè)效能進(jìn)行評(píng)估,以實(shí)現(xiàn)預(yù)測(cè)規(guī)則參數(shù)的優(yōu)化。

1 地震活動(dòng)和預(yù)測(cè)指標(biāo)選取

本文數(shù)據(jù)使用了中國(guó)地震臺(tái)網(wǎng)地震目錄[13]和內(nèi)蒙古臺(tái)網(wǎng)觀測(cè)資料[14],地震目錄資料的完整性及可靠性分析是地震學(xué)研究的基礎(chǔ)工作。相關(guān)研究顯示,1975年以來(lái)中國(guó)大陸地震目錄完整性能夠達(dá)到ML3.0以下[15],而本文研究時(shí)空范圍內(nèi)地震目錄的完整性震級(jí)可以達(dá)到ML2.5以下[6,16],因此本文所采用的數(shù)據(jù)資料滿足了對(duì)地震目錄完整性研究的需求。對(duì)于地震活動(dòng)性分析中地震目錄是否需要?jiǎng)h除余震尚存爭(zhēng)議,一些觀點(diǎn)認(rèn)為余震掩蓋了地震作為獨(dú)立事件的統(tǒng)計(jì)特征,使得地震活動(dòng)的時(shí)空分布偏離了正?；顒?dòng)狀態(tài),增加了地震活動(dòng)的非平穩(wěn)性,因此在地震活動(dòng)分析中應(yīng)設(shè)法刪除余震活動(dòng)的影響[17]；另一些觀點(diǎn)認(rèn)為余震對(duì)概率地震危險(xiǎn)性計(jì)算結(jié)果不會(huì)產(chǎn)生顯著影響,因此在地震區(qū)劃或一般性地震危險(xiǎn)性分析中可考慮不用刪除余震[18-19]。本文采用未刪除余震的目錄。

2019年以來(lái),晉冀蒙交界地區(qū)出現(xiàn)了ML3.0以上地震顯著增強(qiáng)的現(xiàn)象,在考察范圍(圖1)(圓心:113.5°E,40.2°N,半徑:200 km)內(nèi),截至2020年5月1日已發(fā)生ML≥3.0地震共計(jì)20次。自1970年以來(lái),該區(qū)域年平均地震個(gè)數(shù)為11次,2020年后的地震活動(dòng)水平顯著高于背景水平。圖1右下角的子圖為2015年以來(lái)上述區(qū)域ML≥3.0以上地震的滑動(dòng)年頻次曲線可以看出,2019年上半年起滑動(dòng)年頻次逐漸升高,2020年3月達(dá)到20次。

圖1 2019年以來(lái)增強(qiáng)活動(dòng)區(qū)域(圓形內(nèi)部) ML≥3.0地震Fig.1 ML≥3.0 earthquakes in seismicity enhancement region (inside the circle) since 2019

根據(jù)地震活動(dòng)性和中強(qiáng)地震的時(shí)空相關(guān)性,選取晉冀蒙交區(qū)域(38°～42°N,111°～115°E)作為預(yù)測(cè)范圍,主要考察MS5.0以上地震的回溯情況。

圖2為預(yù)測(cè)區(qū)域內(nèi)目標(biāo)地震分布圖。圖3給出1970至今ML≥3.0地震滑動(dòng)年頻次時(shí)序曲線,圖中使用倒三角形標(biāo)注了目標(biāo)地震的時(shí)刻,可見滑動(dòng)年頻次的升高和中強(qiáng)地震的發(fā)生具有一定關(guān)聯(lián)性。

圖2 目標(biāo)地震分布圖Fig.2 Distribution of target earthquakes

圖3 增強(qiáng)區(qū)域ML≥3.0地震年滑動(dòng)頻次圖Fig.3 The annual slip frequency of ML≥3.0 earthquakes in the enhancement region

滑動(dòng)年頻次反映了中期地震活動(dòng)強(qiáng)弱變化,初步選取預(yù)測(cè)窗長(zhǎng)為2年,閾值14,將年滑動(dòng)頻次超過14次后的2年作為預(yù)測(cè)時(shí)段。使用R值評(píng)分法[20]對(duì)上述指標(biāo)對(duì)于初步選取預(yù)測(cè)指標(biāo)進(jìn)行檢驗(yàn)。R值的含義是扣除了隨機(jī)概率的預(yù)報(bào)成功率,當(dāng)R值大于零,表明預(yù)測(cè)方法有效,R值越高則預(yù)測(cè)效果越好。統(tǒng)計(jì)樣本的數(shù)量會(huì)影響R值的置信水平,許紹燮[20]編制了保證97.5%置信水平的最低R值表(R0)。如果實(shí)際計(jì)算R值大于表中對(duì)應(yīng)值,則認(rèn)為該R值有至少97.5%置信度。在上述預(yù)測(cè)規(guī)則下,采用的預(yù)測(cè)指標(biāo)報(bào)準(zhǔn)率為87%,時(shí)空占有率為42%,R值為0.45,對(duì)應(yīng)R0為0.37,可見其預(yù)測(cè)效能是顯著的。

1994年和2014年研究區(qū)域內(nèi)出現(xiàn)指標(biāo)達(dá)到閾值的情況,而在兩年內(nèi)未出現(xiàn)5級(jí)的中強(qiáng)地震,屬于虛報(bào)情況。地震預(yù)測(cè)指標(biāo)主要基于震前指標(biāo)異常和目標(biāo)地震發(fā)生的關(guān)聯(lián)性,并不是直接的因果關(guān)系,僅在現(xiàn)象上存在統(tǒng)計(jì)學(xué)關(guān)系,因此存在虛報(bào)現(xiàn)象。本文主要關(guān)注晉冀蒙交界地區(qū)中小地震活動(dòng)性和中強(qiáng)地震,而地震前兆存在“長(zhǎng)程關(guān)聯(lián)”和“遠(yuǎn)程觸發(fā)”的情況[21],特別是在相同的構(gòu)造體系中,晉冀蒙交界地區(qū)、河套地震帶、鄂爾多斯西北緣同屬鄂爾多斯地塊北緣構(gòu)造體系,因此其地震活動(dòng)存在相關(guān)性。1994年的地震活動(dòng)增強(qiáng)現(xiàn)象可能和發(fā)生在河套地震帶的1996年包頭6.4級(jí)地震有關(guān),2014年的地震活動(dòng)增強(qiáng)可能和鄂爾多斯地塊西北緣的2015年阿拉善左旗5.8級(jí)地震有關(guān)。

2 評(píng)估預(yù)測(cè)規(guī)則的優(yōu)化

雖然初步選取的預(yù)測(cè)規(guī)則可以通過統(tǒng)計(jì)檢驗(yàn),但其預(yù)測(cè)規(guī)則是否最優(yōu)還有待驗(yàn)證。預(yù)測(cè)規(guī)則確定后,根據(jù)指標(biāo)值和預(yù)測(cè)值進(jìn)行預(yù)測(cè),就可以確定報(bào)警范圍(時(shí)間段),如果一次目標(biāo)地震發(fā)生在報(bào)警范圍內(nèi),則該次地震預(yù)測(cè)成功,反之預(yù)測(cè)失敗,根據(jù)歷史地震的回溯情況可以評(píng)價(jià)預(yù)測(cè)規(guī)則的優(yōu)劣。本文選用的指標(biāo)對(duì)應(yīng)的預(yù)測(cè)規(guī)則主要包括閾值和預(yù)測(cè)窗長(zhǎng),閾值是預(yù)測(cè)指標(biāo)的觸發(fā)條件,當(dāng)頻次超過閾值進(jìn)入有效預(yù)測(cè)時(shí)段；預(yù)測(cè)窗長(zhǎng)為有效預(yù)報(bào)期的長(zhǎng)度。為了進(jìn)一步優(yōu)化預(yù)測(cè)規(guī)則,以R值為目標(biāo)函數(shù),改變預(yù)測(cè)規(guī)則的參數(shù),在參數(shù)空間內(nèi)搜索最優(yōu)的規(guī)則(對(duì)應(yīng)最高的R值)。同時(shí)使用ROC曲線和Molchan圖表法考察指標(biāo)的效果。

2.1 R值評(píng)分法

上文已經(jīng)對(duì)R值的含義進(jìn)行了介紹，其值定義為:

(1)

式中：等式右端第一項(xiàng)為報(bào)準(zhǔn)率,第二項(xiàng)為時(shí)空占有率。

首先固定預(yù)測(cè)窗長(zhǎng)為2年以考察閾值變化時(shí)的預(yù)測(cè)效果,計(jì)算其不同閾值的報(bào)準(zhǔn)率和時(shí)空占有率,并繪制了R值隨閾值變化的曲線，同時(shí)給出R0值曲線作為參考(圖4)。由圖4 看出,當(dāng)閾值從0增加至100,R值在-0.2至0.45范圍變化；當(dāng)閾值為14次時(shí),有最高的R值0.45,并且超過R0。因此,初步選取的閾值14次是合理的。

圖4 R值隨閾值變化Fig.4 Variation of R value with the threshold value

為了同時(shí)選取最優(yōu)的預(yù)測(cè)窗長(zhǎng)和閾值,繪制R值在閾值-預(yù)測(cè)窗長(zhǎng)構(gòu)成的二維參數(shù)空間的分布(圖5)。由圖5可以看出，較高的R值分布在閾值為10～14次,窗長(zhǎng)為100～400天范圍。結(jié)果表明初步選取的閾值14次是合理的,但預(yù)測(cè)窗長(zhǎng)2年不是最優(yōu)的,最佳預(yù)測(cè)窗長(zhǎng)在200～400天之間,因此預(yù)測(cè)窗長(zhǎng)設(shè)為1年更加合理。因此最后確定閾值為14次,預(yù)測(cè)窗長(zhǎng)為1年。

圖5 R值隨閾值和窗長(zhǎng)的變化Fig.5 Changes of R value with the threshold value and window length

2.2 ROC曲線

接受者操作特性曲線(ROC曲線),也稱感受性曲線,反映敏感性與特異性之間關(guān)系,經(jīng)常被應(yīng)用于預(yù)測(cè)模型的評(píng)估。ROC曲線中,橫坐標(biāo)為誤報(bào)率(FPR),也稱特異性,表示假陽(yáng)性的比例,橫坐標(biāo)接近零時(shí)準(zhǔn)確率較高;縱坐標(biāo)為敏感度(TPR),表示稱為真陽(yáng)性比例,縱坐標(biāo)較大時(shí)準(zhǔn)確率較好。ROC曲線將圖劃分成為兩部分,曲線下方部分的面積稱為AUC(Area Under Curve),AUC較大時(shí)準(zhǔn)確性較好。曲線越接近左上角(X越小,Y越大),預(yù)測(cè)準(zhǔn)確率越高。

將預(yù)測(cè)窗長(zhǎng)設(shè)置為1年和2年,在固定預(yù)測(cè)窗長(zhǎng)的情況下,改變閾值范圍(0～100),計(jì)算FPR和對(duì)應(yīng)的TPR,繪制出兩種窗長(zhǎng)的ROC曲線(圖6)。

圖6 兩種預(yù)測(cè)窗長(zhǎng)的ROC曲線Fig.6 ROC curves of the two forecast window lengths

圖6顯示,黑色實(shí)線(窗長(zhǎng)1年)在紅色虛線上方(窗長(zhǎng)2年),表明預(yù)測(cè)窗長(zhǎng)為1年的預(yù)測(cè)效果好于預(yù)測(cè)窗長(zhǎng)為2年。

在ROC曲線基礎(chǔ)上發(fā)展出來(lái)的Molchan圖表法[22]被應(yīng)用于地震前兆的預(yù)報(bào)效能的評(píng)估,將時(shí)空占有率τ和漏報(bào)率ν的關(guān)系在二維空間中表示出來(lái),當(dāng)異常閾值由最高改變到最低時(shí),時(shí)空占有率τ從0上升到1,漏報(bào)率ν從1下降到0；當(dāng)最大預(yù)測(cè)成功率和付出最小的時(shí)空代價(jià)時(shí)達(dá)到最優(yōu)的預(yù)測(cè)效能。將τ-ν曲線和顯著性水平α[式(2)]或概率增益Gain值[式(3)]對(duì)比,可以有效地評(píng)價(jià)預(yù)測(cè)方法的效果和顯著性。τ-ν曲線和坐標(biāo)軸包圍的面積可以反映不同預(yù)測(cè)方法的優(yōu)劣[23-24]。圖7為預(yù)測(cè)窗長(zhǎng)360天固定后,改變閾值繪制的Molchan曲線。其結(jié)果顯示，概率增益達(dá)到4,顯著性水平小于1%,但閾值下調(diào)至一定值時(shí),漏報(bào)率可以降至0,顯著性水平小于0.01,但概率增益略高于2,則最高的概率增益對(duì)應(yīng)的漏報(bào)率為0.6。

(2)

式中:N為地震數(shù)；h為擊中數(shù)。

(3)

圖7 Melchan圖表(預(yù)測(cè)窗長(zhǎng)360天)Fig.7 Melchan chart (The forecast window length is 360 days)

3 討論

R值的含義是扣除了隨機(jī)概率的預(yù)報(bào)成功率。本文以R值為評(píng)估標(biāo)準(zhǔn),在預(yù)測(cè)規(guī)則的參數(shù)空間搜索,以獲得對(duì)應(yīng)最高R值的最優(yōu)規(guī)則。概率增益是指標(biāo)觸發(fā)和未觸發(fā)情況下,地震發(fā)生概率的比值,也可以由報(bào)準(zhǔn)率和時(shí)空占有率計(jì)算得出。概率增益表示異常出現(xiàn)后目標(biāo)地震發(fā)生概率的放大倍數(shù),具有明確的統(tǒng)計(jì)學(xué)意義,也可以作為預(yù)測(cè)指標(biāo)的評(píng)估標(biāo)準(zhǔn)。如果選取概率增益作為評(píng)估標(biāo)準(zhǔn),也可獲得對(duì)應(yīng)于最大概率增益的最優(yōu)預(yù)測(cè)規(guī)則,結(jié)果與以R值為評(píng)估標(biāo)準(zhǔn)的結(jié)果會(huì)有差異。圖8給出概率增益在參數(shù)空間的分布。可以看出，概率增益最高值在預(yù)測(cè)窗長(zhǎng)1年和閾值45次附近。概率增益反映了收益成本比,在擊中率和時(shí)空占有率都很小時(shí)可能達(dá)到很高的概率增益,此時(shí)R值并未達(dá)到最高。如果選取概率增益最大時(shí)的預(yù)測(cè)規(guī)則,可能出現(xiàn)較高的漏報(bào)率。這種差異說(shuō)明,概率增益對(duì)成本敏感,而R值在成本和命中率這兩方面達(dá)到某種折中,因此在實(shí)際工作中選取R值高的指標(biāo)更合理。

圖8 概率增益值隨閾值和窗長(zhǎng)的變化Fig.8 Changes of probability gain value with the threshold value and window length

4 結(jié)論

本文對(duì)晉冀蒙地區(qū)中小地震活動(dòng)性和中強(qiáng)地震的對(duì)應(yīng)關(guān)系進(jìn)行分析,在比較歷史地震目錄和中強(qiáng)地震資料的基礎(chǔ)上建立了基于地震頻次的預(yù)測(cè)指標(biāo),并優(yōu)化預(yù)測(cè)規(guī)則。

在規(guī)則優(yōu)化中,以R值為目標(biāo)函數(shù),分別考察固定窗長(zhǎng)后R值隨閾值的變化和同時(shí)改變閾值/預(yù)測(cè)窗長(zhǎng)時(shí)R值在二維參數(shù)空間的變化。固定窗長(zhǎng)為2年后,最優(yōu)閾值為14次;同時(shí)改變窗長(zhǎng)和閾值,逐點(diǎn)計(jì)算預(yù)測(cè)規(guī)則參數(shù)空間，其結(jié)果顯示,較高的R值分布在閾值為10～14次,窗長(zhǎng)為100～400天的范圍。根據(jù)R值在預(yù)測(cè)參數(shù)空間的掃描結(jié)果,最終確定了以晉冀蒙交界增強(qiáng)區(qū)域ML≥3.0地震滑動(dòng)年頻為預(yù)測(cè)指標(biāo),閾值為14次,預(yù)測(cè)窗長(zhǎng)1年,對(duì)應(yīng)的R值達(dá)到0.5。

R值和概率增益也都可以作為指標(biāo)的評(píng)估函數(shù)。概率增益在擊中率和時(shí)空占有率都很小時(shí)可能達(dá)到高值,如果選取概率增益最大時(shí)的預(yù)測(cè)規(guī)則可能出現(xiàn)較高的漏報(bào)率。R值在成本和命中率這兩方面達(dá)到某種折中,在實(shí)際工作中選取R值高的指標(biāo)更合理。

使用ROC曲線和和Molchan圖表法對(duì)預(yù)測(cè)方法進(jìn)行評(píng)估。ROC曲線顯示預(yù)測(cè)窗長(zhǎng)為1年時(shí)的預(yù)測(cè)效果優(yōu)于窗長(zhǎng)2年的；Molchan圖表顯示預(yù)測(cè)指標(biāo)的最優(yōu)概率增益達(dá)到4,對(duì)應(yīng)的顯著性水平小于1%。

本文建立的基于單個(gè)預(yù)測(cè)指標(biāo)的預(yù)測(cè)方法可應(yīng)用于地震預(yù)測(cè)指標(biāo)體系中。本文所采用評(píng)估和優(yōu)化的方法亦可以應(yīng)用于各個(gè)學(xué)科的地震預(yù)測(cè)指標(biāo)評(píng)估中。