楊建文 葉 泵 陳 佳 高 瓊 王 軍

1 中國地震科學(xué)實驗場大理中心，云南省大理市濱海大道， 671000 2 云南大理滇西北地殼構(gòu)造活動野外科學(xué)觀測研究站，云南省大理市濱海大道，671000

現(xiàn)階段，GNSS站點主要分為連續(xù)站和區(qū)域站，區(qū)域站網(wǎng)點分布密集，空間分辨較高，可通過定期監(jiān)測獲取板塊間的相對運動，來描述地殼趨勢運動，解釋區(qū)域地球動力學(xué)的背景等[1-2]；連續(xù)站觀測從區(qū)域站的定期變?yōu)闇蕦崟r，短周期信息較為豐富，具有較高的時間分辨率，能實時監(jiān)測地殼運動的動態(tài)變化特征。GNSS連續(xù)觀測的特點決定了其應(yīng)該具有捕獲地震短臨異常前兆的能力。

不少學(xué)者在GNSS用于短臨預(yù)測方面進行嘗試，取得了一些進展。張曉亮等[3]通過數(shù)學(xué)形變?yōu)V波的方式處理單條基線的變化，從中分離出與地震時間有關(guān)的形變分量，用于提取地震前的短期形變異常。吳云等[4]利用GNSS基準站的時間序列，采用小波變換分析部分震例認為，GNSS基準站連續(xù)觀測序列的低頻段中，地震前6個月或稍長時段有一定的前兆異常出現(xiàn)，可作為地震中短期預(yù)測的依據(jù)。張燕等[5]利用小波變換對GNSS時間序列進行分析，捕獲了大震前1～2 a各站點運動趨勢的改變，并認為該變化可能是一種地震前兆信號。邵德盛等[2]基于云南地區(qū)28個GNSS連續(xù)站2011-01～2016-04數(shù)據(jù)，以21°～29°N、97°～106°E范圍內(nèi)的20個5級以上地震為例，提出面應(yīng)變綜合預(yù)測指標法，獲取了適用于云南地區(qū)的地震短臨異常識別指標。

利用GNSS觀測資料，可進行應(yīng)變參數(shù)的求解，獲取區(qū)域應(yīng)變場。應(yīng)變參數(shù)主要包括面應(yīng)變和最大剪應(yīng)變等，各應(yīng)變參數(shù)實際上就是對不同特性形變信息的一種客觀分離[6]。面應(yīng)變參數(shù)能直接反映區(qū)域擠壓(收縮)或拉張(膨脹)強弱的特性，而最大剪應(yīng)變是反映剪切應(yīng)變強度的參數(shù)[6]。震源機制解研究表明，云南地區(qū)的M≥5.0地震多為走滑型地震，地震的破裂類型多含有剪切破裂，震前應(yīng)變積累包含一定的剪切應(yīng)變積累(走滑型地震以剪應(yīng)變?yōu)橹?，而非走滑型地震也帶有一定的走滑成?[6]，因而，地震發(fā)生前最大剪應(yīng)變會表現(xiàn)出一定的異常特征，利用最大剪應(yīng)變資料提取適用于云南地區(qū)的短臨異常指標并進行短臨預(yù)測應(yīng)該是可行的。

鑒于此，本文在前人研究的基礎(chǔ)上，基于云南地區(qū)的43個GNSS連續(xù)站2011-01～2018-12數(shù)據(jù)，分別提取適用于云南地區(qū)的面應(yīng)變及最大剪應(yīng)變短臨異常指標，并以其間發(fā)生的27個M≥5.0地震為樣本，采用R值評分法[7]對異常指標進行預(yù)測效能評價，以期將GNSS連續(xù)觀測資料更好地應(yīng)用于云南地區(qū)的震情跟蹤和預(yù)測預(yù)報工作。

1 數(shù)據(jù)收集整理

本文的研究范圍(21°～29°N，97°～107°E)內(nèi)共包含43個GNSS連續(xù)站。GNSS數(shù)據(jù)為中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(http:∥www.cgps.ac.cn)提供的站點位移時間序列原始數(shù)據(jù)，數(shù)據(jù)起止時間為2011-01～2018-12。該平臺采用美國麻省理工學(xué)院(MIT)和加州大學(xué)圣地亞哥分校Scripps海洋研究所(SIO)研制的GAMIT/GLOBK軟件(10.4版)進行數(shù)據(jù)處理。研究區(qū)2011-01～2018-12期間共發(fā)生27個M≥5.0地震。GNSS站點(點位信息來源于中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(http:∥www.cgps.ac.cn))及地震震中分布(地震信息來源于中國地震臺網(wǎng)(http:∥news.ceic.ac.cn))如圖1所示，相關(guān)地震參數(shù)見表1。

表1 云南地區(qū)27個M≥5.0地震事件列表 (2011-01～2018-12)

加號為GNSS連續(xù)站點；大圓點為M≥5.0地震震中；數(shù)字為格網(wǎng)編號；細線為活動斷裂分布圖1 GNSS連續(xù)站、地震震中分布及格網(wǎng)劃分結(jié)果Fig.1 The grid division and the distribution of GNSS continuous stations and earthquakes

2 GNSS應(yīng)變場計算與短臨異常指標提取

2.1 GNSS應(yīng)變場計算

使用云南地區(qū)43個GNSS連續(xù)站的站點位移時間序列原始數(shù)據(jù)進行應(yīng)變場的求解?？紤]到GNSS連續(xù)站分布的不均勻性，在獲取了各個站點位移時間序列的基礎(chǔ)上，采用克里金插值方法對點位數(shù)據(jù)進行1°×1°格網(wǎng)化插值(格網(wǎng)劃分結(jié)果及編號見圖1)，獲取均勻分布于80個格網(wǎng)的位移時間序列，再通過位移時間序列求取區(qū)域應(yīng)變場[2,8]，具體解算方法見文獻[8]。

通過計算可獲取面應(yīng)變、最大剪應(yīng)變、主應(yīng)變等應(yīng)變參數(shù)[6]。其中，面應(yīng)變θ和最大剪應(yīng)變γmax計算公式為[6,8]：

(1)

2.2 短臨異常指標提取

本文的短臨異常指標提取采用邵德盛等[2]提出的應(yīng)變異常識別方法。在§2.1應(yīng)變場計算的基礎(chǔ)上，分別提取80個格網(wǎng)的應(yīng)變(面應(yīng)變、最大剪應(yīng)變)時間序列，并進行線性去趨勢處理，以獲取長趨勢變化背景下的短期異常波動。對去除趨勢項后的時間序列，劃定±2倍標準差作為單條曲線的異常指標線[2]。

圖2為24號格網(wǎng)的面應(yīng)變和最大剪應(yīng)變時間序列圖。圖2(a)顯示，該格網(wǎng)存在一定的趨勢壓縮，年壓縮速率為-0.40×10-8，rms為±1.07。通過線性去趨勢處理，并設(shè)定±2倍標準差作為異常閾值，可以明顯看到(圖2(b))，該格網(wǎng)的面應(yīng)變時間序列在騰沖M5.2、寧蒗M5.7地震前存在較為明顯的超指標異常，即面應(yīng)變時間序列值超過±2倍標準差線。圖2(c)顯示，該格網(wǎng)的最大剪應(yīng)變年速率為3.41×10-8，rms為±0.94。通過去除趨勢項，并設(shè)定±2倍標準差作為異常閾值，可以明顯看到(圖2(d))，該格網(wǎng)的最大剪應(yīng)變時間序列在騰沖M5.2、彝良M5.7、洱源M5.0地震前存在較為明顯的超指標異常，即最大剪應(yīng)變時間序列值超過±2倍標準差線。

圖2 24號格網(wǎng)的面應(yīng)變和最大剪應(yīng)變時間序列Fig.2 Time series diagram of surface strain and maximum shear strain of grid No.24

在對應(yīng)變格網(wǎng)時序進行線性去趨勢處理后，需要對單個格網(wǎng)時間序列的映震能力進行評價。評價方式以±2倍標準差線作為閾值，超出該指標線的均為異常，并根據(jù)異常出現(xiàn)時間與地震事件的對應(yīng)關(guān)系，對面應(yīng)變、最大剪應(yīng)變曲線的映震能力進行評分[2]。

由于需獲取短臨異常信息，因此地震事件間對應(yīng)的時間間隔設(shè)定為不超過90 d(3個月)。假設(shè)在某時段內(nèi)第i個格網(wǎng)面應(yīng)變及最大剪應(yīng)變時序曲線報準次數(shù)為Ni，漏報次數(shù)為Li，虛報次數(shù)為Xi，那么該曲線的最終評價得分為[2]：

(2)

經(jīng)過以上處理，可分別獲取相應(yīng)時段80個格網(wǎng)的預(yù)測能力評分。若在該時段內(nèi)第i、j、l、m格網(wǎng)出現(xiàn)異常，則該時段的異常總評分為[2]：

S=Si+Sj+Sl+Sm

(3)

在此基礎(chǔ)上通過逐個時間點對所有格網(wǎng)進行時空掃描，可分別獲取云南地區(qū)的面應(yīng)變及最大剪應(yīng)變異常評分時間序列，稱之為面應(yīng)變及最大剪應(yīng)變短臨異常指標。當該異常評分值超出一定閾值后(此處設(shè)定2倍均值線作為異常閾值)，表示異常指標出現(xiàn)，認為云南地區(qū)存在發(fā)生M≥5.0地震的危險性。圖3為分別提取的2011-01～2018-12云南地區(qū)的面應(yīng)變及最大剪應(yīng)變短臨異常指標。

3 異常指標的預(yù)測效能評價

在分別提取了云南地區(qū)的面應(yīng)變及最大剪應(yīng)變短臨異常指標后(圖3)，需要結(jié)合相應(yīng)的地震事件對異常指標進行預(yù)測效能評價?，F(xiàn)階段應(yīng)用最為廣泛的預(yù)測效能評價方法是1973年由許紹燮院士提出的R值評分法[9]。R值評分法主要依據(jù)預(yù)測指標實際預(yù)測地震的有震報準率c和預(yù)測占時率(預(yù)測時間占有率)b[7]，反映的是預(yù)測方法或指標與地震的相關(guān)程度。R值的計算公式為：

圖3 云南地區(qū)的面應(yīng)變及最大剪應(yīng)變短臨異常指標(2011-01～2018-12)Fig.3 Short-term anomaly indicators of surface strain and maximum shear strain in Yunnan area (2011-01-2018-12)

R=c-b=

(4)

式中，R=1表示全報對，R=0表示預(yù)測沒有起作用。R值越大，預(yù)測效果越好。當預(yù)測效能R>R0(R0為保證97.5%置信水平的最低R值)時[10]，表明所評估的預(yù)測方法通過統(tǒng)計檢驗，高于隨機預(yù)測效能。

采用R值評分法對云南地區(qū)的面應(yīng)變及最大剪應(yīng)變短臨異常指標進行效能評價，兩種指標不同預(yù)測窗長的R值遍歷曲線如圖4所示。由圖可知，最大剪應(yīng)變短臨異常指標的R值分布較面應(yīng)變短臨異常指標更為集中，R值主要分布于預(yù)測窗長250 d以內(nèi)，最優(yōu)R值為0.360，所對應(yīng)的預(yù)測窗長為30 d。而面應(yīng)變短臨異常指標的R值主要分布于預(yù)測窗長450 d以內(nèi)，最優(yōu)R值為0.280，相對應(yīng)的預(yù)測窗長為40 d和70 d。

圖4 不同預(yù)測窗長的R值遍歷曲線Fig.4 R value traversal curves of different prediction window lengths

由于本文獲取的是云南地區(qū)的短臨異常指標，因此，在對異常指標進行預(yù)測效能評價時，預(yù)測窗長設(shè)定為不超過90 d(3個月)。在圖4的基礎(chǔ)上，我們分別對面應(yīng)變短臨異常指標進行預(yù)測窗長為30 d、40 d、60 d、70 d、90 d的預(yù)測效能評價，對最大剪應(yīng)變短臨異常指標進行預(yù)測窗長為30 d、60 d、90 d的預(yù)測效能評價，具體效能評價結(jié)果見表2和表3。

表2 面應(yīng)變短臨異常指標預(yù)測效能評價結(jié)果

表3 最大剪應(yīng)變短臨異常指標預(yù)測效能評價結(jié)果

由表2可知，面應(yīng)變短臨異常指標在預(yù)測窗長為30 d、40 d、60 d、70 d、90 d時都通過了效能檢驗。其中，當預(yù)測窗長設(shè)置為30 d時，預(yù)測效能評分R為0.23，R0為0.21，R>R0，對于實際發(fā)生的27個地震，該指標報對地震15次，漏報12次，準確率為55.56%，漏報率為44.44%；當預(yù)測窗長為60 d時，預(yù)測效能評分R為0.21，R0為0.20，R>R0，該指標報對地震19次，漏報8次，準確率為70.37%，漏報率為29.63%；當預(yù)測窗長為90 d時，預(yù)測效能評分R為0.24，R0為0.19，R>R0，該指標報對地震23次，漏報4次，準確率為85.19%，漏報率為14.81%；當根據(jù)最優(yōu)R值設(shè)置預(yù)測窗長為40 d和70 d時，該指標預(yù)測準確率分別為66.67%和81.48%，漏報率分別為33.33%和18.52%。通過進一步對比發(fā)現(xiàn)，從30 d到90 d，該指標預(yù)測準確率逐漸升高，到90 d時預(yù)測準確率達到85.19%。

由表3可知，當最大剪應(yīng)變短臨異常指標的預(yù)測窗長設(shè)置為30 d時(最優(yōu)R值對應(yīng)的預(yù)測窗長)，預(yù)測效能評分R為0.36，R0為0.20，R>R0，該預(yù)測指標通過了效能檢驗，對于實際發(fā)生的27個地震，該指標報對地震21次，漏報6次，準確率為77.78%，漏報率為22.22%；當預(yù)測窗長為60 d時，預(yù)測效能評分R為0.25，R0為0.19，R>R0，預(yù)測指標通過了效能檢驗，該指標報對地震23次，漏報4次，準確率為85.19%，漏報率為14.81%；當預(yù)測窗長為90 d時，預(yù)測效能評分R為0.12，R0為0.19，R

對比表2和表3可知，當預(yù)測窗長設(shè)置為30 d時，面應(yīng)變短臨異常指標的預(yù)測效能評分R為0.23，R/R0為1.12，預(yù)測準確率為55.56%；最大剪應(yīng)變短臨異常指標的預(yù)測效能評分R為0.36，R/R0為1.82，預(yù)測準確率為77.78%，前者的R值、R/R0值和準確率均低于后者。當預(yù)測窗長設(shè)置為60 d時，面應(yīng)變短臨異常指標的預(yù)測效能評分R為0.21，R/R0為1.03，預(yù)測準確率為70.37%；最大剪應(yīng)變短臨異常指標的預(yù)測效能評分R為0.25，R/R0為1.30，預(yù)測準確率為85.19%，前者的R值、R/R0值和準確率均低于后者。當預(yù)測窗長為90 d時，二者的預(yù)測準確率相同，但是最大剪應(yīng)變短臨異常指標未通過效能檢驗。

4 討 論

根據(jù)R值評分規(guī)則，在97.5%置信度下，當R>R0值時，預(yù)測效能檢驗通過，因此，當異常指標超過控制線時可認為，云南地區(qū)未來30～90 d內(nèi)發(fā)生M≥5.0地震的可能性較高(預(yù)測窗長為90 d時，最大剪應(yīng)變短臨異常指標除外)?？v觀不同預(yù)測窗長的R值大小發(fā)現(xiàn)，兩種異常指標的R值都不太高，主要原因是M≥5.0地震在云南地區(qū)的自然發(fā)生概率較高。

從兩種異常指標的預(yù)測效能評價結(jié)果來看，當預(yù)測窗長設(shè)置為30 d或60 d時，無論是R值、R/R0值還是預(yù)測準確率，最大剪應(yīng)變短臨異常指標均優(yōu)于面應(yīng)變短臨異常指標。究其原因，可能與本文所涉及到的27個M≥5.0地震多為走滑型地震，最大剪應(yīng)變對走滑型地震前的應(yīng)力變化過程更為敏感有關(guān)。

結(jié)合R值、R/R0值、預(yù)測準確率，對兩種異常指標進行不同預(yù)測窗長的對比發(fā)現(xiàn)，面應(yīng)變短臨異常指標更適合于預(yù)測窗長為90 d的短臨預(yù)測，而最大剪應(yīng)變短臨異常指標更適合預(yù)測窗長為60 d的短臨預(yù)測。

結(jié)合圖3中的短臨異常指標線來看，兩種異常指標與地震事件之間均存在較好的對應(yīng)關(guān)系。2014年云南地區(qū)相繼發(fā)生05-30盈江6.1級、08-03魯?shù)?.5級、10-07景谷6.6級地震。3次M≥6.0地震發(fā)生前，面應(yīng)變短臨異常指標分別于05-09、07-18、09-16出現(xiàn)明顯的超指標異常，最大剪應(yīng)變短臨異常指標分別于05-07、07-11、09-22出現(xiàn)明顯的超指標異常。兩種異常指標都對3次地震作出了準確的預(yù)測，應(yīng)該說3次M≥6.0地震的發(fā)生分別對兩種異常指標進行了3次很好的檢驗。

5 結(jié) 語

本文通過對云南地區(qū)43個GNSS連續(xù)站2011-01～2018-12站點位移時間序列原始數(shù)據(jù)進行深加工處理，在計算應(yīng)變參數(shù)格網(wǎng)時間序列的基礎(chǔ)上，分別提取適合云南地區(qū)M≥5.0地震短臨預(yù)測的面應(yīng)變及最大剪應(yīng)變短臨異常指標；并基于其間發(fā)生的27個M≥5.0地震，采用R值評分法對異常指標進行預(yù)測效能評價，形成以下認識：

1)面應(yīng)變、最大剪應(yīng)變都是重要的應(yīng)變參數(shù)，基于二者分別提取的面應(yīng)變及最大剪應(yīng)變短臨異常指標物理意義明確，兩種指標對云南地區(qū)的M≥5.0地震的發(fā)震時間均具有較好的指示意義。

2)通過對兩種指標進行不同預(yù)測窗長的效能檢驗發(fā)現(xiàn)，當預(yù)測窗長設(shè)置為30 d或60 d時，在相同預(yù)測窗長下，最大剪應(yīng)變短臨異常指標的預(yù)測準確率高于面應(yīng)變短臨異常指標；通過對同一異常指標進行不同預(yù)測窗長的對比發(fā)現(xiàn)，面應(yīng)變短臨異常指標的最佳預(yù)測窗長為90 d，最大剪應(yīng)變短臨異常指標的最佳預(yù)測窗長為60 d。

3)本文異常指標提取過程中設(shè)定±2倍標準差作為單個格網(wǎng)時間序列映震能力的異常閾值，設(shè)定2倍均值作為異常指標評分的異常閾值，主要是基于提高信噪比和異常識別能力而考慮的。對于選取更高或更低的異常閾值標準下，兩種指標的短臨預(yù)測和效能評價結(jié)果如何，有待進一步研究。

4)現(xiàn)階段，云南地區(qū)的GNSS連續(xù)站點密度較低，空間分辨率亦較低。隨著中國地震科學(xué)實驗場更多項目的開展，云南地區(qū)的GNSS連續(xù)站點數(shù)量將大幅提高，到時可基于更高空間分辨率的GNSS連續(xù)觀測數(shù)據(jù)提取具有更強短臨異常識別能力的綜合指標，這對于云南地區(qū)的震情跟蹤及預(yù)測預(yù)報工作具有重要的意義。

致謝：本文GNSS數(shù)據(jù)來源于中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(http:∥www.cgps.ac.cn)，數(shù)據(jù)分析軟件由云南省地震局信息中心提供，在此一并表示感謝。