時序數據年變信息提取及地震預測指標評估方法研究

劉 琦 張 晶

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

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

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

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

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

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

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

2)信息指標建立。基于計算結果建立2個信息指標，分別為年頻段歸一化平均幅度(ANA)和其他頻段歸一化平均幅度(ONA)。具體定義如下：

(1)

(2)

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

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

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

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

2 實際觀測數據應用

新疆庫爾勒臺周邊中強地震相對較多，以該臺站水平擺傾斜北南測項為例，介紹上述方法的實際應用過程。使用數據時段為2008-01-01～2022-02-22(圖2(a))，該數據頻段幅值比為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前后再次緩慢增大。基于該信息，根據M≥5和M≥6兩個震級檔開展預測指標確定和效能評估。圖3(c)和圖4(c)為不同預報時窗、不同空間范圍下相對最優閾值組合對應的R值分布情況，白色五角星為最優空間范圍和預報時窗所處位置。

由圖3、4可見，5級以上地震的ANA和ONA信息最優預報范圍均在200 km以內，最優預報時窗分別為930 d和840 d。最優預測指標均可正確預測觀測期間發生的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種信息最優指標的R值均大于R0值，統計意義上可信度較高。ONA信息最優預測指標對應的預報占時率和顯著性水平更低，R值和概率增益更高，預測效能相對更優。

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

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

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

3 結 語

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

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