劉 琦 張 晶

1 中國地震局地震預(yù)測(cè)研究所，北京市復(fù)興路63號(hào)，100036

近年來，隨著觀測(cè)儀器的大量布設(shè)和震例資料的不斷累積，地震前兆異常研究獲得越來越多的關(guān)注，其中破年變異常研究較為普遍[1-2]。該類異常過去主要以人工判別為主，為提升客觀性，相關(guān)學(xué)者開展定量提取算法研究，以正常年變背景擬合為基礎(chǔ)，提出包括距平法、最小二乘法、分段最小二乘法、滑動(dòng)傅里葉方法和奇異譜分析方法等典型算法[3-5]。其中，距平法、最小二乘法構(gòu)建的背景年變時(shí)序相對(duì)固定，無法適應(yīng)變化的背景年變；分段最小二乘法通過分段形式解決背景年變變化問題，但會(huì)在分段點(diǎn)處引入階躍；滑動(dòng)傅里葉方法對(duì)非正/余弦波及非平穩(wěn)時(shí)序周期成分的識(shí)別不太穩(wěn)定；奇異譜分析方法存在相移現(xiàn)象，且處理過程中部分參數(shù)需要根據(jù)實(shí)際情況進(jìn)行交互設(shè)置，不利于自動(dòng)化業(yè)務(wù)應(yīng)用。在預(yù)報(bào)效能定量評(píng)估方面，目前的主要依據(jù)包括有震報(bào)準(zhǔn)率、無震報(bào)準(zhǔn)率、漏報(bào)率、虛報(bào)率、時(shí)空占有率等參數(shù)，對(duì)上述參數(shù)進(jìn)行適當(dāng)組合，形成R值評(píng)分、Molchan圖表法、接收者操作特性(ROC)曲線等評(píng)估方法[6-8]。雖然上述方法的側(cè)重點(diǎn)不同，但相互之間具有一定的關(guān)聯(lián)，在實(shí)際預(yù)報(bào)業(yè)務(wù)中得到廣泛應(yīng)用。

本文提出一種基于時(shí)頻分析方法的破年變信息分離提取流程，并基于雙向非對(duì)稱閾值策略，結(jié)合R值評(píng)分及Molchan圖表法構(gòu)建預(yù)測(cè)指標(biāo)確定和效能定量評(píng)估方法。以新疆庫爾勒水平擺傾斜北南測(cè)項(xiàng)的處理分析為例，介紹該方法的實(shí)際應(yīng)用過程。

1 破年變信息提取及預(yù)測(cè)指標(biāo)效能評(píng)估方法

1.1 基于時(shí)頻分析方法的破年變信息提取流程

時(shí)頻分析方法在多組分信號(hào)非穩(wěn)態(tài)時(shí)頻演化過程研究中具備優(yōu)勢(shì)，比較適合破年變信息的分離和提取。本文將時(shí)頻分析類方法中應(yīng)用相對(duì)廣泛的S變換方法[9]作為流程構(gòu)建基礎(chǔ)，該方法具有良好的時(shí)頻聚集性和時(shí)頻分辨率。具體步驟如下：

1)S變換處理。在去臺(tái)階、去突跳點(diǎn)、缺數(shù)插值等預(yù)處理基礎(chǔ)上，針對(duì)日值采樣數(shù)據(jù)進(jìn)行數(shù)據(jù)篩選，選擇觀測(cè)時(shí)長較長且年變顯著的數(shù)據(jù)用于S變換計(jì)算(基于頻段幅值比進(jìn)行年變顯著程度的定量評(píng)估)。為降低邊界效應(yīng)，對(duì)數(shù)據(jù)兩側(cè)進(jìn)行數(shù)據(jù)總長度7.5%的擴(kuò)邊處理。

2)信息指標(biāo)建立?；谟?jì)算結(jié)果建立2個(gè)信息指標(biāo)，分別為年頻段歸一化平均幅度(ANA)和其他頻段歸一化平均幅度(ONA)。具體定義如下：

(1)

(2)

1.2 綜合R值和Molchan圖表的預(yù)測(cè)指標(biāo)確定和效能評(píng)估方法

預(yù)測(cè)指標(biāo)確定和效能評(píng)估最常用的方法為R值評(píng)分法。R值為扣除隨機(jī)概率后的預(yù)報(bào)成功率，R=0表示預(yù)報(bào)無效。當(dāng)同樣的R值基于的地震次數(shù)不同時(shí)，信度也不同，因此還需要考慮R0值。當(dāng)R值大于R0值時(shí)可認(rèn)為，該R值至少有97.5%的置信度[10]。Molchan圖表法同樣考慮了漏報(bào)率、預(yù)報(bào)占時(shí)率等要素，可通過概率增益、顯著性水平等要素對(duì)模型預(yù)測(cè)能力展開進(jìn)一步評(píng)價(jià)。

綜上所述，動(dòng)態(tài)設(shè)置一系列異常閾值，利用漏報(bào)率、預(yù)報(bào)占時(shí)率建立坐標(biāo)，可生成Molchan曲線?；谠撉€可確定特定規(guī)則下(特定震級(jí)范圍、預(yù)報(bào)范圍及預(yù)報(bào)時(shí)窗)的相對(duì)最優(yōu)閾值，即同時(shí)考慮距對(duì)角線最遠(yuǎn)點(diǎn)P1及距原點(diǎn)最近點(diǎn)P2，從中選擇概率增益更大點(diǎn)對(duì)應(yīng)的閾值參數(shù)(概率增益Gain=有震報(bào)準(zhǔn)率/預(yù)報(bào)占時(shí)率，可反映與隨機(jī)概率相比目前概率的提升倍數(shù))。P1實(shí)際對(duì)應(yīng)R值最大點(diǎn)，即漏報(bào)率、預(yù)報(bào)占時(shí)率L1范數(shù)最小點(diǎn)，P2則對(duì)應(yīng)漏報(bào)率、預(yù)報(bào)占時(shí)率L2范數(shù)最小點(diǎn)。當(dāng)Molchan曲線點(diǎn)位于顯著性水平α=2.5%參考線左側(cè)時(shí)，表示該點(diǎn)對(duì)應(yīng)的R值大于R0值(圖1(a))。采用雙向非對(duì)稱二維閾值策略(信息指標(biāo)小于負(fù)向閾值Th1或大于正向閾值Th2，均被判定為異常，Th2>Th1)，可充分考慮正向、負(fù)向異常差異性。同一漏報(bào)率可能會(huì)對(duì)應(yīng)多個(gè)預(yù)報(bào)占時(shí)率，需從中選擇預(yù)報(bào)占時(shí)率最小時(shí)所對(duì)應(yīng)的閾值參數(shù)組合。常用的雙向?qū)ΨQ閾值、單向閾值等一維參數(shù)均為該策略的特例，因此基于雙向非對(duì)稱閾值策略獲得的閾值參數(shù)預(yù)測(cè)效能最優(yōu)(圖1(b))。按照震級(jí)檔對(duì)空間范圍、預(yù)報(bào)時(shí)窗等參數(shù)進(jìn)行調(diào)整，重復(fù)上述相對(duì)最優(yōu)參數(shù)指標(biāo)計(jì)算過程，對(duì)比不同參數(shù)組合下的R值，選擇最大值對(duì)應(yīng)的空間范圍、預(yù)報(bào)時(shí)窗、閾值參數(shù)作為該震級(jí)檔的最優(yōu)預(yù)測(cè)指標(biāo)。若R值相同，則可進(jìn)一步依據(jù)高概率增益、低顯著性水平進(jìn)行篩選。

圖1 預(yù)測(cè)指標(biāo)確定和效能評(píng)估定量方法原理Fig.1 Quantitative methods for predictor determination and efficacy assessment

2 實(shí)際觀測(cè)數(shù)據(jù)應(yīng)用

新疆庫爾勒臺(tái)周邊中強(qiáng)地震相對(duì)較多，以該臺(tái)站水平擺傾斜北南測(cè)項(xiàng)為例，介紹上述方法的實(shí)際應(yīng)用過程。使用數(shù)據(jù)時(shí)段為2008-01-01～2022-02-22(圖2(a))，該數(shù)據(jù)頻段幅值比為2.02，具備相對(duì)清晰的年變背景(圖2(b))。

圖2 新疆庫爾勒臺(tái)水平擺傾斜北南測(cè)項(xiàng)觀測(cè)曲線及年變顯著程度檢測(cè)Fig.2 The NS component of the horizontal pendulum tiltmeter in the Korla station in Xinjiang and the detection of significant degree of annual variation

圖3(a)、圖4(a)分別為其他頻段歸一化平均幅度ONA(周期為14～500 d)以及年頻段歸一化平均幅度ANA(周期為250～500 d)曲線圖。ONA可顯示短周期信號(hào)的變化情況，由圖3(a)可見，2015年之前周期信號(hào)變化相對(duì)劇烈，2015年之后變化相對(duì)減緩；ANA可反映背景年變信號(hào)演化過程，由圖4(a)可見，ANA指標(biāo)幅度緩慢增大，并在2012-07前后達(dá)到峰值，之后逐漸減小，2016-08前后再次緩慢增大?；谠撔畔?，根據(jù)M≥5和M≥6兩個(gè)震級(jí)檔開展預(yù)測(cè)指標(biāo)確定和效能評(píng)估。圖3(c)和圖4(c)為不同預(yù)報(bào)時(shí)窗、不同空間范圍下相對(duì)最優(yōu)閾值組合對(duì)應(yīng)的R值分布情況，白色五角星為最優(yōu)空間范圍和預(yù)報(bào)時(shí)窗所處位置。

由圖3、4可見，5級(jí)以上地震的ANA和ONA信息最優(yōu)預(yù)報(bào)范圍均在200 km以內(nèi)，最優(yōu)預(yù)報(bào)時(shí)窗分別為930 d和840 d。最優(yōu)預(yù)測(cè)指標(biāo)均可正確預(yù)測(cè)觀測(cè)期間發(fā)生的5次5級(jí)以上地震(2012-01-08和碩5.0級(jí)、2012-06-15輪臺(tái)5.4級(jí)、2013-03-29昌吉5.6級(jí)、2015-06-25托克遜5.4級(jí)和2016-01-14輪臺(tái)5.3級(jí)地震)，2種信息最優(yōu)指標(biāo)的R值均大于R0值，統(tǒng)計(jì)意義上可信度較高。ONA信息最優(yōu)預(yù)測(cè)指標(biāo)對(duì)應(yīng)的預(yù)報(bào)占時(shí)率和顯著性水平更低，R值和概率增益更高，預(yù)測(cè)效能相對(duì)更優(yōu)。

6級(jí)以上地震的ANA和ONA信息最優(yōu)預(yù)報(bào)范圍均在250 km以內(nèi)，且最優(yōu)預(yù)報(bào)時(shí)窗均為7 d。最優(yōu)預(yù)測(cè)指標(biāo)均可正確預(yù)測(cè)觀測(cè)期間發(fā)生的2次6級(jí)以上地震(2012-06-30新源6.6級(jí)和2016-12-08呼圖壁6.2級(jí)地震)。由于該時(shí)段內(nèi)發(fā)震次數(shù)較少，對(duì)應(yīng)的R0值較高，因此只有當(dāng)ANA信息最優(yōu)指標(biāo)的R值大于R0值時(shí)，才會(huì)在統(tǒng)計(jì)意義上有較高的可信度。ONA信息最優(yōu)指標(biāo)對(duì)應(yīng)的預(yù)報(bào)占時(shí)率較高，導(dǎo)致R值偏低，未通過顯著性水平2.5%的檢驗(yàn)。

圖3 庫爾勒臺(tái)站周邊ONA信息的5級(jí)以上地震預(yù)測(cè)指標(biāo)及對(duì)應(yīng)效能Fig.3 Forecasting index and corresponding performance of ONA information for M≥5 earthquakes around the Korla station

圖4 庫爾勒臺(tái)站周邊ANA信息的6級(jí)以上地震預(yù)測(cè)指標(biāo)及對(duì)應(yīng)效能Fig.4 Forecasting index and corresponding performance of ANA information for M≥6 earthquakes around the Korla station

3 結(jié) 語

本文基于S變換時(shí)頻分析方法，構(gòu)建破年變信息分離提取新流程，該流程在提供短周期信號(hào)(ONA)的同時(shí)還可提供背景年變信號(hào)(ANA)的演化過程。采用雙向非對(duì)稱二維閾值策略，結(jié)合R值評(píng)分及Molchan圖表法重新構(gòu)建預(yù)測(cè)指標(biāo)確定及效能評(píng)估流程。閾值參數(shù)空間擴(kuò)展及評(píng)估因素的增加使得獲取的指標(biāo)相對(duì)更優(yōu)、評(píng)估更加合理。新疆庫爾勒水平擺傾斜北南測(cè)項(xiàng)的實(shí)際應(yīng)用效果較好，該測(cè)項(xiàng)ANA對(duì)臺(tái)站周邊250 km內(nèi)6級(jí)以上地震具有較好的預(yù)測(cè)效能，ONA則對(duì)200 km范圍內(nèi)5級(jí)以上地震具有相對(duì)更優(yōu)的預(yù)測(cè)效果。

需要說明的是，本文方法獲得的預(yù)測(cè)指標(biāo)均基于概率統(tǒng)計(jì)，該類指標(biāo)大多會(huì)面臨地震樣本數(shù)較少導(dǎo)致的穩(wěn)定性問題，且單個(gè)地震實(shí)況不一定與整體概率預(yù)測(cè)結(jié)果一致，因此在實(shí)際震情研判過程中還需要綜合考慮多方面因素。此外，為提升預(yù)測(cè)指標(biāo)的合理性，后續(xù)可圍繞干擾排除、震例相關(guān)性等進(jìn)行更深入的研究。