劉 琦 張 晶

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

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

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

1 破年變信息提取及預(yù)測指標效能評估方法

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

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

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

2)信息指標建立?；谟嬎憬Y(jié)果建立2個信息指標，分別為年頻段歸一化平均幅度(ANA)和其他頻段歸一化平均幅度(ONA)。具體定義如下：

(1)

(2)

1.2 綜合R值和Molchan圖表的預(yù)測指標確定和效能評估方法

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

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

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

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

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

圖2 新疆庫爾勒臺水平擺傾斜北南測項觀測曲線及年變顯著程度檢測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可顯示短周期信號的變化情況，由圖3(a)可見，2015年之前周期信號變化相對劇烈，2015年之后變化相對減緩；ANA可反映背景年變信號演化過程，由圖4(a)可見，ANA指標幅度緩慢增大，并在2012-07前后達到峰值，之后逐漸減小，2016-08前后再次緩慢增大?；谠撔畔ⅲ鶕?jù)M≥5和M≥6兩個震級檔開展預(yù)測指標確定和效能評估。圖3(c)和圖4(c)為不同預(yù)報時窗、不同空間范圍下相對最優(yōu)閾值組合對應(yīng)的R值分布情況，白色五角星為最優(yōu)空間范圍和預(yù)報時窗所處位置。

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

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

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

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

3 結(jié) 語

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

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