余 娜 張曉清 袁伏全 趙燕杰
(中國西寧 810000 青海省地震局)
近年來,基于統計物理學與地震學的圖像信息學(Pattern Informatics,簡稱PI)算法在地震較活躍地區的中長期預測研究中得到了廣泛應用。Rundle 等(2002)將PI 算法應用到美國南加州地區地震的預測研究中,認為該方法對加州地區地震活動有較好的預測效果;國內一些研究者也將PI 算法應用到中國臺灣(Chen et al,2005)、中國大陸西部(蔣長勝等,2008;Zhang et al,2010,2013;Jiang et al,2010,2011,2013;張小濤等,2014;袁伏全等,2016)等地區地震的研究中,認為PI 算法的預測效果顯著優于隨機概率法。Zhang 等(2013)和張小濤等(2014)利用PI 算法分別對2008 年汶川8.0 級、2014 年于田7.3 級、2013 年蘆山7.0 級地震進行地震危險性的回溯性研究,結果表明模型計算參數的選取對PI 算法的預測結果有一定影響,對于7 級以上大地震,選擇較大的網格尺度和較長的預測時間窗可能會取得更好的預測效果;袁伏全等(2016)利用PI 算法對青海地區中強地震進行地震危險性的回溯性研究,通過采用不同的網格尺度進行回溯性預測檢驗,認為對于青海地區采用0.2°×0.2°網格尺度時預測效果更好。
印度板塊對歐亞板塊的擠壓作用使青藏塊體區域構造運動強烈,地震頻發,中國大陸7 級以上強震主要發生在青藏塊體及鄰近區域(吳哲,2018)。本文研究區(30.0°—41.0°N,88.0°—105.0°E)包括羌塘活動地塊、巴顏喀拉活動地塊、柴達木活動地塊和祁連活動地塊等4 個Ⅱ級活動地塊。活動塊體邊界帶是我國大陸強震的主體帶和集中區,因此,研究活動地塊邊界帶的強震趨勢顯得尤為重要(張浪平等,2010)。本文在前人工作的基礎上,利用PI 算法對研究區MS≥6.0 地震進行回溯性研究,對不同時間尺度下的預測結果進行ROC 檢驗,并給出適合該地區的模型計算參數。
PI 算法的實質是通過對研究區地震活動增強和平靜進行分析,計算并預測在中長期時間尺度上可能發生中強地震的概率。主要步驟:首先,將研究區按一定尺度劃分空間網格,對落入網格的不小于截止震級的地震事件構建時間序列;然后,對地震活動強度變化進行歸一化處理,計算每個網格發震概率;最后,發震概率減去背景概率得到發震概率較高的區域,即“PI 熱點分布圖”(Tiampo et al,2002;Nanjo et al,2006;蔣長勝等,2008)。具體計算步驟如下。
(1)對研究區進行空間網格劃分,每一網格為xi。
(2)對落入網格且不小于截止震級的地震事件的網格構建1 個時間序列Ni(t),其中,Ni(t)為中心坐標和相鄰的8 個格點(Moor 近鄰)(Moore,1962)在時間t 內單位時間的地震數目;t0為研究區地震目錄的起始時刻;tb為滑動變化各時間序列的起始時刻;t1為地震活動中異常學習的起始時刻;t2和t3分別為預測時段的起始時刻和終止時刻(圖1)。

圖1 空間網格劃分及時間序列構建Fig.1 Division of spatial grid and construction of time series
(3)地震活動強度函數Ii(tb,t)為從時刻tb到t 的單位時間內發生在網格i 中的不小于截止震級的平均地震數

(4)對不同時間段的地震活動強度求平均,再除以標準偏差,進行標準化處理

(5)計算地震活動強度函數平均變化量,以較少隨機擾動的影響

(6)未來強震發生概率Pi(t0,t1,t2)為地震強度函數平均變化量的平方。第i 網格的概率值減去所有網格概率的平均值得到強震發生在第i 網格的概率

將ΔPi(t0,t1,t2) >0的格點視為研究區的地震熱點。在PI 圖像中將地震活動性較強或地震頻次較高的30%格點畫出(陳建志等,2013),由其反映預測時間窗內空間上發生目標地震概率較高的區域,用lg(ΔP/ΔPmax)的值表示其發震概率。
ROC 檢驗方法是在國際合作項目(CSEP)中進行科學預測試驗的方法之一,其目的是檢驗預測結果的優劣程度。ROC 檢驗既要考慮“命中率”,同時也要考慮“虛報率”。有意義的預測必須是ROC 值大于0.5。ROC 值為ROC 曲線與隨機預測曲線所圍成的有效面積,ROC 值越大,預測效能越好。
選取青藏塊體為(30.0°—41.0°N,88.0°—105.0°E)研究對象(圖2)。歷史上該區域曾發生過多次7 級以上強震。本文使用中國地震臺網中心提供的1970—2019 年8 月全國ML≥2.0 地震目錄,共計20 896 次地震。其中,ML2.0—2.9 地震3 286 次;ML3.0—3.9 地震14 611 次;ML4.0—4.9 地震2 307 次;ML5.0—5.9 地震605 次;ML6.0—6.9地震74 次;ML≥7.0 地震13 次。Jiang 等(2011)基于川滇地區強震的回溯性研究討論了強余震對PI 算法的影響,認為余震對PI 算法計算結果的影響時間不超過1 年;Tiampo 等(2002)認為震后余震反映了區域高應力的釋放。因此,本文沒有對地震目錄進行剔除余震處理。

圖2 研究區空間位置矩形框為研究區域Fig.2 The spatial location of the study area
截止震級Mc是PI 算法中重要的參數之一,而截止震級Mc的選取與該地區的最小完備震級有關(Woessner et al,2005)。在計算中,若截止震級Mc選取過高,則參與計算的地震數目會減少;若截止震級Mc選取過低,則有些區域因監測能力較弱而記錄不到地震。因此,本文將截止震級Mc設定為ML3.0,這既能保證充足的地震數據量,又能滿足截止震級應至少小于目標震級2 個震級單位的要求(Holliday et al,2005,2006)。
一 些 研 究 者(Rundle et al,2002;Chen et al,2005;Nanjo et al,2006;Jiang et al,2010;袁伏全等,2016)利用PI 算法對中國大陸西部地震進行回溯性檢驗,認為網格參數選取為0.2°×0.2°,預測時間窗選為8 a 時,PI 算法的預測效果較好。本文參照前人的研究結果(Rundle et al,2002;Chen et al,2005;Nanjo et al,2006;Jiang et al,2010;袁伏全等,2016),選取0.2°×0.2°空間網格尺度,以30 d 為步長連續向前滑動,目標震級≥MS6.0。研究區域中強地震頻發,2009 年8 月至2019 年8 月研究區共發生9 次MS≥6.0 地震(圖3,表1),為了考察研究區MS≥6.0 地震在不同模型參數下的回溯性預測效果,在研究區范圍、網格尺度、截止震級以及滑動步長不變的情況下,選取預測時間窗為10 a、5 a、3 a 進行回溯性研究,模型參數設置結果見表2。

圖3 2009 年8 月至2019 年8 月研究區MS≥6.0 地震震中分布Fig.3 Distribution of MS≥6.0 earthquakes since August 2009

表1 2009 年8 月至2019 年8 月研究區MS≥6.0 地震參數Table 1 MS≥6.0 earthquakes parameters in the study area since 2009

表2 模型計算參數設置Table 2 Selection of calculation parameters
根據表2 的模型計算參數,得到了3 個回溯性預測檢驗時段的PI 地震熱點分布圖(圖4)。

圖4 研究區回溯性PI 預測圖像及ROC 檢驗結果(a)、(b)2009 年8 月1 日—2019 年8 月1 日;(c)、(d)2014 年8 月1 日—2019 年8 月1 日;(e)、(f)2016 年8 月1 日—2019 年8 月1 日Fig.4 PI retrospective forecast and ROC test for the research area
選定研究區地震目錄的起始時刻 t0為1970 年1 月1 日,分別選取地震活動中異常學習的起始時刻t1為1999 年8 月1 日、2009 年8 月1 日和2013 年8 月1 日,預測時段的起始時刻t2分別為2009 年8 月1 日、2014 年8 月1 日 和2016 年8 月1 日,預測時段的終止時刻t3為2019 年8 月1 日,因此共3 個回溯性預測時間段。圖4(a)為10 a 預測時間窗(2009 年8 月1 日至2019 年8 月1 日)內的地震熱點分布。由圖4(a)可見,異常出現在大柴旦—宗務隆山斷裂、鄂拉山斷裂、達布遜湖—霍布遜湖斷裂、祁連山北緣斷裂、岷江斷裂和汶川—北川斷裂附近。根據預測規則可知,這些區域是預測時間窗內空間上相對危險的地區,發生目標地震的概率較高。實際上,該預測時間窗內研究區共發生了9 次MS≥6.0地震,其中有4 次地震震中落在由PI 算法計算得出的地震熱點區域,分別為2009 年海西6.4級、2013 年蘆山7.0 級、2016 年門源6.4 級、2017 年九寨溝7.0 級地震。
從5 a 預測時間窗(2014 年8 月1 日至2019 年8 月1 日)內的PI 預測圖像[圖4(c)]可見,2016 年門源6.4 級、2017 年九寨溝7.0 級地震震中均位于PI 預測圖像的地震熱點區域;從3 a 預測時間窗(2016 年8 月1 日至2019 年8 月1 日)內的PI 預測圖像[圖4(e)]可見,2017 年九寨溝7.0 級地震震中位于PI 預測圖像的地震熱點區域。
采用ROC 方法回溯性檢驗不同預測時間窗的預測效果,發現10 a 時間尺度的ROC 值為0.610,5 a 時間尺度的為0.725,3 年時間尺度的為0.793。可以看出,預測時間窗內的強震震中基本都在PI 預測圖像的地震熱點區域,而且PI 算法預測結果明顯優于隨機概率法[圖4(b),4(d),4(f)]。雖然可供參考的強震事件數較少,但在對PI 算法的回溯性檢驗中發現,2009 年海西6.4 級、2016 年門源6.4 級、2013 年蘆山7.0 級、2017年九寨溝7.0 級地震震中都落在地震熱點的叢集區,這表明PI 算法對青藏塊體的強震具有一定的預測效能,并且將預測時間窗設定為3 a 時間尺度時其預測效能較好。
2009 年8 月至2019 年8 月研究區共發生9 次MS≥6.0 地震,其中,6 次發生在活動塊體邊界帶上(圖3),分別為2009 年海西6.4 級、2010 年玉樹7.1 級、2013 年蘆山7.0級、2014 年康定6.3 級、2016 年門源6.4 級、2017 年九寨溝7.0 級地震。考慮到研究區強震大部分發生在活動塊體邊界帶,且活動塊體邊界帶也是我國大陸強震的主體帶和集中區,因此定性分析了活動塊體邊界帶“目標”地震與地震熱點叢集間的關系,發現青藏塊體MS≥6.0 地震與活動塊體邊界帶上地震熱點間對應關系較好,活動塊體邊界帶的“目標”地震都發生在地震熱點附近區域,而在活動塊體內部的MS≥6.0 地震震前無地震熱點情況。僅就PI 算法而言,很難對該現象給出具體原因。但在地震中長期預測中應用PI 算法時,將其計算結果與活動塊體邊界帶結合起來考慮強震相對危險區域,在地震預測研究中可能有一定參考價值。
本文利用中國地震臺網中心提供的1970 年以來全國地震目錄,在分析研究區最小完備震級的基礎上,應用PI 算法,選取空間網格尺度0.2°×0.2°,截止震級為ML3.0,目標震級≥ML6.0,將預測時間窗設定為10 a、5 a、3 a 不同時間尺度進行回溯性研究,并對預測結果進行ROC 檢驗,發現PI 算法對青藏塊體強震有一定預測效能,預測效果明顯優于隨機概率法,并且將預測時間窗設定為3 a 時間尺度時預測效能較好。