何 暢,路 茜,龍 鋒,宮 悅,李 麗
(1.四川省地震局,四川 成都 610041;2.山西省地震局,山西 太原 030021)
長期的地球物理場觀測與地震預報實踐表明,由于地球的不可入性以及地震孕育過程的復雜性(陳運泰,2008),單一手段的地震預測方法往往很少有較高的信度以支撐起對某次強震的準確預測,采用算法聯(lián)合多種監(jiān)測預報手段進行綜合概率預測是一條可行探索的途徑,已有眾多基于不同算法的實驗性工作得以開展并取得成效,如支持向量機回歸方法(蔣淳等,2006;衛(wèi)定軍等,2014)、概率增益建模(邱玉榮,2011)、遺傳神經(jīng)網(wǎng)絡(陳秀瓊等,2005)、模糊聚類法(狄莉莎等,2000)等。利用多手段綜合預報,可以“抹平”各手段的隨機預測噪聲,同時可在一定程度上增強前兆信號的強度,以提高預測水平。
在設計綜合概率預測方法的過程中,如何保證各單一指標以不同的權重(信度、預測效能等)參與進來,并產(chǎn)出符合先驗的、相對的強震發(fā)震期望(概率或某些可度量的數(shù)值)是這一方法的核心部分。本文在對四川地區(qū)開展地震預測實踐中,從貝葉斯概率出發(fā),推導出了單個指標的發(fā)生強震的后驗概率,采用聯(lián)合概率計算所有異常的綜合概率,并對2013~2019年四川區(qū)域的3次6級以上強震進行了回溯性研究,并嘗試提取不同震級的綜合概率值。
不同于古典頻率學派,貝葉斯概率公式側(cè)重于基于已有的先驗知識預測某件事情發(fā)生的概率(四川大學數(shù)學系高等數(shù)學教研室,1990),已經(jīng)在參數(shù)回歸(Hoff,2009)、風險決策(王學軍等,2003)、機器學習(茅偉強,2008)等方面得到廣泛應用。針對單個地震預測指標的貝葉斯概率的公式為:
式中:Y表示異常,D表示地震。P(Y)和P(D)分別是該指標異常出現(xiàn)以及地震發(fā)生的先驗概率,P(Y|D)是地震發(fā)生前必然存在異常的條件概率,是似然值。而所求的P(D|Y)是異常出現(xiàn)后必然出現(xiàn)地震的條件概率。由于地震預測是二元分類問題——即有異常、無異常、有地震、無地震,因此依據(jù)全概率公式,異常出現(xiàn)的總概率P(Y)共有兩項,可以寫為:
式中:N表示無地震。為了方便計算,同時考慮到異常和地震發(fā)生的低概率性,將它們離散到時間上(以天為單位),用時間的占有率來代替各自的概率。設總的研究時段為Tall,地震發(fā)生次數(shù)(或天數(shù))為Neq,異常出現(xiàn)次數(shù)(或天數(shù))為Na,異常出現(xiàn)后發(fā)生預測中的地震次數(shù)(即“報準”的次數(shù)或天數(shù))為Nb,則:
由(3)、(4)兩式可得:
另有:
P(Y|N)為未發(fā)生地震時出現(xiàn)異常的概率,可表示為未發(fā)生地震時未出現(xiàn)異常的概率的補集,即:
式中:X為未出現(xiàn)異常的時段,t為預報指標中單次預報時長。聯(lián)合(5)~(7)式,最終可得:
上式表明,某個指標異常出現(xiàn)后對應地震的后驗概率可從其具體的預報指標和歷史驗證效能中獲得。進一步將式(8)上下同時除以Na,則可得到:
式中:Rc為報準率,從上式可以看出高的報準率和短的預報時長會明顯提升預報概率,這與實際相符。另外也可看出式(9)存在分母為0的函數(shù)奇點。由于Rc∈[0,1],而t往往大于1,因此分母是否為0取決于t與Neq/Na的大小關系,也就是說,當預報時長與研究區(qū)總地震數(shù)和異常出現(xiàn)次數(shù)的比例接近時,概率值可能為無窮;且當這一比例明顯大于預報時長t時,概率值為負數(shù)。
針對四川地區(qū)目前仍在正常運作的各觀測項進行了6級以上歷史震例的回溯性研究,提取了能通過R值檢驗(許紹燮,1973①許紹燮.1973.地震技術資料匯編:震兆分析一例.北京:《地震戰(zhàn)線》.;許紹燮等,1989)的預測指標,最終獲得87個測項指標的貝葉斯后驗概率計算結果(表1)。在常規(guī)的地震預測實踐中,我們更期望獲得月這一時間尺度上的強震發(fā)生概率,因此在通過式⑻計算概率時,往往將指標的預測時間t歸算到月。在表1中,地震活動性的部分異常如空區(qū)、條帶、震群、b值、震源機制解一致性空間掃描、以及宏觀異?;蚱渌R時出現(xiàn)但具備一定預報意義等共計9項指標由于無法溯及過往總數(shù),無法計算P(D|Y)值,依相應學科專家依經(jīng)驗給定數(shù)值;另有9項指標由于P(D|Y)值計算為負數(shù)或無窮大,人為給定該值為0.010(表1中帶②的指標項)。
通過對P(D|Y)進行統(tǒng)計,所有指標(圖1a)與不考慮人為給定概率的指標(圖1b)在統(tǒng)計分布上并無明顯差異,后驗概率在0.1~0.2是優(yōu)勢概率區(qū)間,表明四川地區(qū)大部分預測指標在月尺度上對6級以上地震的預測并不具備強約束。具有較高預測概率的3個指標分別是巴塘川-58泉(305K)水溫預測概率為0.33、普格川-65泉水溫預測概率為0.6和康定二道橋川-55泉水溫預測概率為0.75,均屬于地下流體測項。從表1中震例回溯的統(tǒng)計數(shù)據(jù)來看,這3個測項的報準率并不高(占50%或更低),顯著低于部分跨斷層形變測量測項,如爾烏水準3-A測邊(85%以上)、龍燈壩基線A-D測邊(100%)等,且有著偏高的虛報率,但由于其較短的預報期限(3個月),使其在月尺度上有著較高預測概率值。
圖1 P(D|Y)值所有指標(a)和不考慮認為給定概率的指標(b)統(tǒng)計柱狀圖
表1 四川地區(qū)強震指標及相關震例
續(xù)表1
續(xù)表1
續(xù)表1
考慮到歷史強震震例存在部分地球物理場觀測項失效的問題,這里選擇針對2013年以來四川地區(qū)的3次6級以上震例進行了概率分析(表1)。統(tǒng)計結果顯示,2013年蘆山7.0級地震前存在異常40項,其中預報期限在6個月以內(nèi)的中短期異常12項。為了綜合考量所有異常指標對未來強震的預測概率,采用了聯(lián)合概率的方法,設由式(8)計算得到的第i個指標的預測概率為P ci,則n個異常指標在月尺度上6級以上強震的發(fā)震概率可表示為:
根據(jù)式(10),若假定P ci=0.01,需229項異??墒筆值達到0.9,而當P ci=0.1時,則需22項異常。
以2013年蘆山7.0級地震而言,震前40項異常計算出來的發(fā)震概率為0.992。繪制發(fā)震概率隨異常項數(shù)變化的曲線繪制出來(圖2中藍色實線),可以看出第7項異常(例如康定二道橋川-55泉水溫,P c=0.75;普格川-65泉水溫,P c=0.60等)對提升總體概率貢獻極大,從0.45提升至0.86。當異常項數(shù)達到15項而概率超過0.90后,概率曲線平緩,這是因為式(10)為極限為1的冪函數(shù),這意味著盡管異常項數(shù)增加,但最終發(fā)震概率值變化幅度較小,通常在小數(shù)點后兩位變化,這無益于異常信號的識別和預報指標的提取。為此,本文試圖從兩個方面解決這個問題:(a)僅使用中短期異常。蘆山7.0級地震前12項中短期異常的發(fā)震概率同樣達到了0.92(圖2中紅色實線),這樣既保證了概率值可達到較高的水平,同時相對較少的異常數(shù)量使得相鄰概率值之間的區(qū)分度較高。這種方法并不能保證未來隨著監(jiān)測能力的提升,更多的中短期異常項數(shù)依舊會使得概率曲線變得平坦。(b)考慮到概率值在0~1之間變化,通過反正弦函數(shù)將概率值轉(zhuǎn)換為角度,即:
圖2 蘆山地震前概率值和角度值隨異常項數(shù)變化曲線(實線為概率,虛線為角度)
可在一定程度上消弭概率曲線的平坦效應,增加概率的區(qū)分度。從圖2可以看出,對于蘆山7.0級地震前的40項異常,相較于概率值在第15~40項異常之間僅有0.07的變化幅度,相應區(qū)間的角度變化值有15°(藍色虛線),方便于指標提取時閾值的確立。蘆山7.0級地震前40項異常的角度值為82.79°(藍色虛線),如果只考慮12項中短期異常,則角度值為68.43°(紅色虛線)。從實際應用的效果來看,四川地區(qū)背景性異常長期維持在30項左右,這一數(shù)量使得發(fā)震概率P值始終在0.9左右,不利于區(qū)分當前強震的緊迫性。因此,僅采用6個月以內(nèi)的中短期異常來計算概率符合當前按月來追蹤震情形勢的需求,且概率變幅較大有助于強震緊迫性的區(qū)分。
本文收集整理了2014年康定6.3級和2019年長寧6.0級地震前的中短期異常,并計算了它們的發(fā)震概率值和概率角度值,其中康定6.3級地震前中短期異常5項,概率和概率角度值分別為0.52和31.10°,而長寧6.0級地震前中短期異常4項,概率和概率角度值分別為0.48和26.08°。
為了建立四川地區(qū)強震綜合概率預測指標,從貝葉斯概率出發(fā),推導出了單個指標出現(xiàn)異常后強震發(fā)生的后驗概率;基于聯(lián)合概率,引入多項指標出現(xiàn)異常后的綜合概率,并對蘆山7.0級、康定6.3級和長寧6.0級等強震開展了回溯性研究。獲得的結論如下:(1)單個指標的貝葉斯后驗概率的大小與報準率正相關、與預報時長負相關,這是符合經(jīng)驗直覺的。但該概率同時也被預報時長和總地震數(shù)及總預報數(shù)的比例的相互大小所約束,在計算過程中要規(guī)避極小可能出現(xiàn)的無窮或負值。(2)四川地區(qū)強震指標的后驗概率大多在0.2以內(nèi),說明大多數(shù)指標短期預測效能一般,對未來強震約束不強。(3)聯(lián)合概率是一個極限為1的冪函數(shù),當異常數(shù)較多時,少數(shù)異常的增減在概率上的區(qū)分度不高,不便于閾值的選取和指標的提取。為此,建議只采用6個月以內(nèi)的中短期異常進行聯(lián)合概率計算,又或者,通過反正弦函數(shù)將概率值轉(zhuǎn)換為角度值以獲得更好的區(qū)分度。(4)采用以上規(guī)則,震例回溯研究統(tǒng)計出:蘆山7.0級地震前有12項中短期異常,計算得到的聯(lián)合概率為0.92,概率角度值為68.43°;康定6.3級地震前中短期異常5項,概率和概率角度值分別為0.52和31.10°;而長寧6.0級地震前中短期異常4項,概率和概率角度值分別為0.48和26.08°。(5)通過這三次較少的震例,本文嘗試性提出了四川地區(qū)6級以上地震的綜合概率預報指標為:中短期異常在4項以上,聯(lián)合概率值大于0.45,概率角度值大于25°;而7級以上地震的綜合概率預報指標為:中短期異常在10項以上,聯(lián)合概率值大于0.90,概率角度值大于65°。(6)四川地區(qū)選取6級以上強震震例較少,且地震發(fā)生和地球物理場的監(jiān)測能力均存在時空不均勻性,因此需要開展更多的震例回溯總結和積累工作,以形成更穩(wěn)健的綜合概率預報指標體系。
致謝:本文在數(shù)據(jù)收集整理過程中,得到了四川地震臺張致偉、官致君、蘇琴、祁玉萍、馬伶俐、芮雪蓮、林圣杰、王迪等提供的大力幫助;中國地震局地球物理研究所蔣長勝研究員提供了部分思路,在此一并致謝。