江小英,孔祥增,郭躬德,李 南,林 嶺

1(福建師范大學(xué) 數(shù)學(xué)與信息學(xué)院,福州 350117)

2(福建農(nóng)林大學(xué) 計算機與信息學(xué)院,福州 350002)

過去幾十年,通過研究大量震例的遙感數(shù)據(jù),研究者們發(fā)現(xiàn)大地震發(fā)生前常有各種異常現(xiàn)象出現(xiàn).如地表潛熱通量 (Surface Latent Heat Flux,SLHF)[1],熱紅外(ThermaL Infrared,TIR)異常[2],電離層總電荷量(TOtal Electron Content,TEC)[3]和射出長波輻射(Outgoing Longwave Radiation,OLR)[4]等,在地震發(fā)生之前都會出現(xiàn)明顯的變化.Qin 等[5]分析2010年新西蘭7.1 級地震發(fā)生前后3 個月內(nèi)的SLHF 數(shù)據(jù)的變化,發(fā)現(xiàn)在地震發(fā)生前一個月左右,震中的東北方向出現(xiàn)了一個孤立的高達(dá)160 W/m2的SLHF 正異常.結(jié)合GPS 位移觀測和構(gòu)造背景,確定SLHF 異?？赡芘c地殼運動時上升的地下熱物質(zhì)有關(guān),從而解釋了震中東北部,北島中心和南島西南部的局部溫度升高現(xiàn)象.張璇等[6]分析了尼泊爾8.1 級地震前熱紅外遙感信息,發(fā)現(xiàn)震前出現(xiàn)明顯的熱紅外亮溫異常變化.Aggarwal[7]分析了2010年4月13–14日青海地區(qū)發(fā)生的地震中電離層中總電荷數(shù)的變化情況,發(fā)現(xiàn)電離層總電荷數(shù)的異常變化或許可以作為這次地震相關(guān)前兆特征.Ouzounov 等[8]研究OLR 數(shù)據(jù),并發(fā)現(xiàn)震前會有連續(xù)的異常變化發(fā)生.

OLR 是地球-大氣系統(tǒng)透過大氣頂部向宇宙空間發(fā)射的主要波長在4～120 μm 的輻射.目前很多研究者研究OLR 數(shù)據(jù)與地震前兆的關(guān)系,并取得了一定的成果.Kong 等[4]采用基于Martingale 理論的方法分析了汶川和廬山兩次地震OLR 數(shù)據(jù)中的異常信息,發(fā)現(xiàn)在汶川地震和蘆山地震發(fā)生前后OLR 數(shù)據(jù)都會顯著增大,并在震后一段時間后恢復(fù)正常.Xiong 和Shen[9]提出了一種統(tǒng)計學(xué)方法,分析了從2007年9月1日至2015年5月23日,約3376 個震例的OLR 異常信息,統(tǒng)計結(jié)果表明震前靠近地震中心的地方有異?，F(xiàn)象出現(xiàn),并且這些異常與該地區(qū)即將發(fā)生的地震有關(guān).康春麗等[10]提出用渦度場的方法分析昆侖山口西8.1 級地震前后長波輻射強度OLR 的變化,結(jié)果顯示,地震活動過程中,不論是衛(wèi)星遙感長波輻射場強度,還是亮度溫度都出現(xiàn)了較明顯的異常變化.孫珂等[11]采用RST(Robust Satellite Technique)算法和異常指數(shù)算法,分別對2015年4月25日和5月12日尼泊爾發(fā)生的兩次大地震前后,衛(wèi)星長波輻射變化特征進(jìn)行分析,發(fā)現(xiàn)運用RST 算法,兩次地震前后都沒有檢測到明顯異常變化,而運用異常指數(shù)算法,兩次地震前都檢測到明顯的熱紅外異常.Xiong 等[12,13]使用小波變換和小波最大值的空間/時間連續(xù)性分析來分析連續(xù)OLR 數(shù)據(jù).他們從OLR 數(shù)據(jù)中發(fā)現(xiàn)了與地震前兆相對應(yīng)的奇點.這些研究成果都表明了在OLR 數(shù)據(jù)中隱藏著重要的前兆信息,地震前后,OLR 數(shù)據(jù)會出現(xiàn)明顯的異常變化.分析OLR 的變化與地震之間的關(guān)系有助于我們進(jìn)一步了解地震相關(guān)機理,給未來的防震減災(zāi)工作提供啟示.

然而目前基于OLR 數(shù)據(jù)的相關(guān)工作,大多只是基于地學(xué)領(lǐng)域的方法分析,算法上有一定的局限性.并且大多只是考慮地震中心所在區(qū)域產(chǎn)生的異常,而對周邊區(qū)域的影響和歷年背景環(huán)境產(chǎn)生的影響,以及季節(jié)性影響并未全面考慮,而這些因素對結(jié)果都會產(chǎn)生較大影響.基于以上兩點考慮,本文提出一種基于OLR數(shù)據(jù)的震前異常分析方法-時空分析法,結(jié)合統(tǒng)計學(xué)中的Martingale 理論[14–16],聚類分析,隨機漫步等計算機相關(guān)算法,并通過汶川地震,日本地震,以及蘆山地震3 個震例詳細(xì)闡述該算法的分析方法及實驗結(jié)果,探索OLR 數(shù)據(jù)中隱藏的前兆信息.

1 相關(guān)數(shù)據(jù)

由于2008年5月12日汶川8 級地震是10年來中國國內(nèi)發(fā)生的最大一次地震,而2013年4月20日蘆山7 級地震的震中區(qū)域與汶川地震在同個研究區(qū)域中,各種地理特性比較相似.2011年3月11日發(fā)生的日本9 級大地震是近年來全球范圍內(nèi)的一次特大地震.這3 個震例比較具有代表性,故本文采用時空分析法分析這3 個震例的OLR 數(shù)據(jù)中隱藏的前兆信息.OLR 數(shù)據(jù)由美國國家海洋和大氣管理局(National Oceanic and Atmosphere Administration,NOAA)衛(wèi)星監(jiān)測[17].具體震例信息如表1所示.

表1 震例相關(guān)信息

2 研究方法

本文提出的基于OLR 數(shù)據(jù)的時空分析法分為6 個部分:1)計算相對于背景數(shù)據(jù)的時間維度變化增量;2)計算空間維度相關(guān)變化增量;3)提取隨機漫步概率特征,4)獲得概率特征序列;5)計算各個點的異常度,統(tǒng)計異常概率;6)計算當(dāng)前點的異常值,與預(yù)設(shè)閾值進(jìn)行比較,判斷異常是否發(fā)生,若發(fā)生,則異常標(biāo)記為1,并重啟算法.OLR 數(shù)據(jù)集是將全球按照2.5?×2.5?分辨率劃分為 72×144個單元.由于OLR 數(shù)據(jù)分為日間數(shù)據(jù)和夜間數(shù)據(jù)兩類,故全球每一天的數(shù)據(jù)可用144×144 單元格表示.若用Data_All表示一天的數(shù)據(jù),則有數(shù)組:

其中,下標(biāo)x和y分別為每個數(shù)據(jù)點的對應(yīng)的單元格橫縱坐標(biāo).研究區(qū)域的橫縱坐標(biāo)一旦確定,從時間維度看,每個區(qū)域的數(shù)據(jù)Data_Datex,y,n也是一個數(shù)組,存儲每天該區(qū)域的OLR 值,則有Data_Datex,y,n=Z∪zi,其中Z={z1,···,zi?1}是同一個單元格中,隨時間變化的連續(xù)數(shù)據(jù)點,代表該區(qū)域的歷史數(shù)據(jù),而zi表示該區(qū)域的當(dāng)前數(shù)據(jù)點.數(shù)據(jù)子集Data_Datex,y,n=Z∪zi可以視為一個連續(xù)變化的時間序列.下面詳細(xì)描述具體的算法流程.

在計算相對于背景數(shù)據(jù)的時間維度變化增量時,先以地震中心所在單元格為中心,獲取 5×5個單元格的震前歷年數(shù)據(jù),根據(jù)式(1)計算地震當(dāng)年的OLR 序列相對于背景數(shù)據(jù)的增量值 ?i,t,獲取地震發(fā)生年的時間維度上的增量序列,稱為背景增量序列.

其中,num表示從可獲取數(shù)據(jù)的起始年份到地震發(fā)生年的總年數(shù).

由于地震發(fā)生時,其周邊環(huán)境也可能相互影響,出現(xiàn)異?，F(xiàn)象.故本文將空間因素考慮在內(nèi),研究地震中心區(qū)域相對于周圍區(qū)域的變化情況.根據(jù)式(2),處理背景增量序列,計算各個數(shù)據(jù)點的空間增量值? (?x,y,i).

得到具有時空相關(guān)性的新增量序列,記為Tn,則有Tn=T∪ti,其中T表示歷史數(shù)據(jù)點,ti表示當(dāng)前數(shù)據(jù)點.

在提取隨機漫步概率特征時,設(shè)置一個滑動窗口,窗口大小記為ws,判斷當(dāng)前數(shù)據(jù)點的值是否大于或等于前一個數(shù)據(jù)點的值,若是,則當(dāng)前監(jiān)測數(shù)據(jù)點的值定義為向右走,否則,當(dāng)前監(jiān)測數(shù)據(jù)點的值定義為向左走.并根據(jù)隨機漫步概率式(3)計算出隨機漫步概率特征序列.

此處,i∈[?ws,ws],q=0.5,而是數(shù)學(xué)中的組合公式.

計算異常度時,采用K-means 聚類方法,計算各個點到聚類中心的距離,記為si.再根據(jù)式(4)統(tǒng)計前i個點的異常概率,記為.

此處,#{}是一個用于統(tǒng)計滿足條件的數(shù)據(jù)點的下標(biāo)個數(shù).如統(tǒng)計滿足sj>si(j=1,2,···,i,i=1,2,···,n)的j的個數(shù),θ是[0,1]區(qū)間的一個隨機值.

最后結(jié)合Martingale 理論,并且為了減小滑動窗口大小設(shè)置不當(dāng)帶來的影響,設(shè)置滑動窗口值取值范圍,求多組滑動窗口計算得到的均值作為最終異常指標(biāo).即根據(jù)式(5)計算異常值,記為CD.

其中,ε是[0,1]區(qū)間的一個任意值.winbegin和winend是窗口的起止取值,本文設(shè)置滑動窗口取值范圍為[30,45].同時預(yù)設(shè)一個閾值h,將各個數(shù)據(jù)點的異常值CD與h進(jìn) 行比較.若從某點開始,CD值突然顯著增大并超過h,則可視為出現(xiàn)異常.此時,需標(biāo)記異常,并重啟算法,重新計算該點之后的數(shù)據(jù)點異常值.即算法重啟規(guī)則如表達(dá)式(6)所示.

綜上所述,震前異常時空分析法的算法流程圖如圖1所示.

3 結(jié)果與分析

實驗時,根據(jù)文獻(xiàn)[4],將閾值h設(shè)置為1000.參考文獻(xiàn)[18],取 ε=0.82.研究的數(shù)據(jù)時間周期為一年,研究數(shù)據(jù)從發(fā)生地震前一年的9月1日到地震發(fā)生當(dāng)年的8月31日.圖2中的豎線表示地震發(fā)生時間.

3.1 汶川地震

2008年5月12日在中國四川省的汶川縣發(fā)生了8 級大地震.地震中心所在位置的經(jīng)緯度是(31.0°N,103.3°E).故地震中心所在單元格的經(jīng)緯度為(28.75°N–31.25°N,102.5°E–105°E),實驗研究區(qū)域為(26.25°N–33.75°N,100°E–107.5°E),總共9 個單元格,地震中心位于第5 單元格中.圖2顯示了9 格單元格中的CD值變化情況.從中可以看出,在地震發(fā)生前,第2,3,4 單元格出現(xiàn)明顯的波峰.并且這些波峰都出現(xiàn)在汶川地震發(fā)生前3 個月內(nèi).這可能是由于地震前,大量能量積累引起的.第1,7 和9 單元格在震前未發(fā)現(xiàn)明顯變化,卻在震后出現(xiàn)了明顯的波峰.這可能是由于地震后大量能量釋放導(dǎo)致的.而第6 和8 單元格在地震前后都沒有明顯變化.

圖1 震前異常時空分析法流程圖

從圖3可以看到,第5 單元格,即地震中心所在單元格也在地震前一個月也有小波峰,但是顯然第5 單元格的異常趨勢沒有第1,2,3,4,7 和9 單元格的異常明顯.

圖2 汶川地震時空分析圖

圖3 汶川地震中心所在單元格的CD 值變化曲線

綜合圖2和圖3,可以看出,地震前后,震中區(qū)域OLR 值會發(fā)生顯著變化,并且地震中心周邊區(qū)域也會出現(xiàn)較大異?，F(xiàn)象.而且較大的異常現(xiàn)象可能并不是出現(xiàn)在震中所在區(qū)域,而是出現(xiàn)在其周邊區(qū)域.OLR 的變化在時間和空間上都有很大的關(guān)聯(lián)性.

3.2 日本地震

2011年3月11日在日本發(fā)生了9 級大地震,給日本東北部造成了毀滅性破壞.震中經(jīng)緯度是(38.1°N,1 4 2.6°E).故地震中心所在單元格的經(jīng)緯度為(36.25°N–38.75°N,142.5°E–145°E),實驗研究區(qū)域為(33.75°N–41.25°N,140°E–147.5°E),地震中心位于第5 單元格中.從圖4可以看出,9 個單元格除了第1 和3 單元格外,都出現(xiàn)了較明顯的波峰.其中第2,5,7,8 和9 都在地震發(fā)生前3 個月內(nèi)出現(xiàn)明顯異常.第4 單元格,雖然波峰出現(xiàn)在地震后,但在地震發(fā)生前兩個月內(nèi)CD值也有明顯波動的趨勢,并且3月11日,即地震當(dāng)天,至3月15日的CD值均大于100.第6 單元格,在地震后一個月內(nèi)出現(xiàn)了明顯的波峰.從圖5可以看出,震中所在區(qū)域,即第5 單元格,在地震發(fā)生前一個月內(nèi)急劇增大,并達(dá)到峰值.在2月10日出現(xiàn)最大值,為444.4.

圖4 日本地震時空分析

圖5 日本地震中心所在單元格的CD 值變化曲線

3.3 蘆山地震

為了探索異常大小與地震震級的關(guān)系,本文也分析了2013年4月20日的蘆山7 級地震.地震中心所在位置的經(jīng)緯度是(30.3°N,103.0°E).故地震中心所在單元格的經(jīng)緯度為(28.75°N–31.25°N,102.5°E–105°E),與汶川地震的震中位于同個單元格內(nèi).研究區(qū)域與汶川地震相同.從圖6可以看出,第1,2,3,5 和6 單元格都在蘆山地震發(fā)生前3 個月內(nèi)出現(xiàn)明顯的波峰,這與前面兩個震例發(fā)現(xiàn)的時間規(guī)律吻合.而第8 單元格的異常出現(xiàn)時間更早,在地震發(fā)生前5 個月已經(jīng)出現(xiàn)明顯的異常波峰.

圖6 蘆山時空分析圖

而從圖7可以看出,蘆山地震發(fā)生前4 個月已經(jīng)出現(xiàn)明顯,并在地震發(fā)生前一個月內(nèi)CD值再次急劇增大,并在震前半個月達(dá)到最大值.第4,7 和9 單元格震前均無明顯變化趨勢.

從圖2,圖4和圖6整體分析,可以發(fā)現(xiàn),蘆山地震的震級明顯小于汶川地震和日本3.11 地震,但震中所在區(qū)域的異常值卻并不比汶川地震和日本地震的異常小.這也可以看出,并非震級越大,異常值就越大.異常與震級并不是簡單的正比關(guān)系,還與其他因素有關(guān).

圖7 蘆山地震中心所在單元格的CD 值變化曲線

此外,結(jié)合圖1,圖3和圖5,從空間角度看,發(fā)現(xiàn)3 個震例的正北方向,即第2 單元格所在方向上,都在震前3 個月出現(xiàn)了明顯的波峰.其他震例的異?，F(xiàn)象是否也具有類似的方位規(guī)律,還有待驗證.

3.4 對比實驗

吳立新等[19]采用概率統(tǒng)計學(xué)中的μ +2σ作為地震異常的衡量指標(biāo),其中μ 表示當(dāng)?shù)卣鹎岸嗄甑腛LR 均值,σ表示方差.他們分析以震中區(qū)域的地震前后3 個月的OLR 序列,發(fā)現(xiàn)很多震例的地震前有些點的值會大于μ +2σ.為了驗證本文所提出的震前異常時空分析法的先進(jìn)性,本文采用與他們相同的方法及相似的實驗設(shè)置分析本文3 個震例的OLR 震前異常.受篇幅影響,僅以日本地震的結(jié)果為例作一說明.我們分析了第5 單元格2011年1月1日到4月1日共3 個月的OLR序列.實驗結(jié)果如圖8所示.震前僅有一個點的值大于μ+2σ,視為異常點.汶川和蘆山兩個震例的結(jié)果中震前也只能看到分散的幾個異常點出現(xiàn),但是這些零星分散的異常點不易被發(fā)現(xiàn)和注意.相對而言,本文所提出的時空分析法也同樣可以檢測到震前異常值,并且異常持續(xù)時間較長,可以很直觀地反映出了異常值的變化趨勢.因此,我們的方法可以更直觀有效的反映變化過程,呈現(xiàn)當(dāng)前數(shù)據(jù)點與它之前的數(shù)據(jù)之間的關(guān)系.

圖8 日本地震中心所在單元格的3 個月OLR 時間序列

3.5 參數(shù)討論

為了比較不同震例,不同年份,不同單元格的實驗結(jié)果,我們在實驗中設(shè)置了相同的參數(shù)初值.本文中設(shè)h=1000,ε=0.82,ws∈[30,45].為了說明本文所提出的方法對參數(shù)初值敏感性,以及參數(shù)初值對實驗結(jié)果的影響,本文設(shè)計了對比實驗,討論了3 個參數(shù)的初值對3 個震例的結(jié)果影響,受篇幅影響,以圖9–圖11的日本震例結(jié)果為例,做一說明.首先,固定ε=0.82,ws∈[30,45],分別令h=250,h=500,h=750,h=1000,h=1500 和h=2000.得到的結(jié)果如圖9所示,可以看出,h的初值大小僅改變波峰極大值出現(xiàn)的時間以及異常值的幅度,如果h設(shè)置太小,則會較早出現(xiàn)波峰極大值,如果h設(shè)置過大,則會導(dǎo)致極大值幅度很大,但是并不會改變異常變化的趨勢.本文設(shè)h=1000,大小較為適中,能較好地反映異常變化趨勢.然后,固定,h=1000,ws∈[30,45],分別令 ε=0.1,ε=0.2,ε=0.4,ε=0.6,ε=0.8和 ε=0.9,得到的結(jié)果如圖10所示,可以看出,ε的初值大小對異常值的變化幅度有影響,若 ε太小,可能檢測不到異常,而從ε=0.6開始,異常變化明顯;當(dāng) ε=0.9時,異常值有所減小,但是整體的變化趨勢并沒有太大影響,不同 ε得到的波峰極大值出現(xiàn)的時間很接近.本文根據(jù)適中原則以及實驗結(jié)果,設(shè)置 ε=0.82,可以檢測到較明顯的變化趨勢.最后固定h=1000,ε=0.82,分別討論ws=25,ws=30,ws=35,ws=40,ws=45 和ws=60.得到的結(jié)果如圖11所示.從中可以看出,不同的滑動窗口也會影響波峰出現(xiàn)的時間以及變化幅度.ws小,CD值的變化幅度也較小,變化趨勢不明顯,隨著ws的增大,CD值的變化幅度也會更大,趨勢更明顯.但是如果滑動窗口太大,會導(dǎo)致窗口起始值前面很多數(shù)據(jù)被浪費.一般來說,選取60 天以內(nèi)的窗口尺度比較合適.本文為了減少不同窗口初值對結(jié)果的影響,分別設(shè)置窗口值為30 到45,計算異常值,再取多窗口尺度的平均異常值作為最終的異常指標(biāo).

圖9 不同h 對日本地震結(jié)果的影響

圖10 不同ε 對日本地震結(jié)果的影響

圖11 不同ws 對日本地震結(jié)果的影響

綜合圖9–圖11可以發(fā)現(xiàn),不同的參數(shù)初值可能對實驗結(jié)果的變化幅度以及異常出現(xiàn)的時間先后有細(xì)微影響,但是不會改變異常的變化過程,整體變化趨勢是類似的.

4 結(jié)論

本文提出一種基于OLR 數(shù)據(jù)的震前異常分析方法-時空分析法,并結(jié)合汶川地震,日本3.11 地震和蘆山地震3 個震例闡述該算法的適用性.從實驗結(jié)果可以看出,時空分析法由于既考慮了不同年份,同一個地區(qū)背景數(shù)據(jù)的影響,又考慮了同一時間,周邊區(qū)域的影響,還考慮了同個季節(jié)內(nèi),歷史數(shù)據(jù)對當(dāng)天數(shù)據(jù)的影響.相對于只考慮其中一種因素影響的研究來說,本文提出的時空分析法考慮得更加全面,分析更合理.

通過分析以震中區(qū)域為中心的9 個單元格,可以發(fā)現(xiàn)幾個規(guī)律:1)地震發(fā)生前后,地震周邊區(qū)域的OLR 也會出現(xiàn)明顯的變化;2)震中區(qū)域出現(xiàn)的異常不一定會比周邊區(qū)域的異常明顯;3)并且一般在震前3 個月就已經(jīng)出現(xiàn)明顯的異常變化趨勢.4)震級大小與異常大小并不是簡單的正比關(guān)系,還與其他因素有關(guān);5)從方位角度分析,發(fā)現(xiàn)3 個震例的正北方向,即第2 單元格所在方向,都在震前3 個月出現(xiàn)了明顯的波峰.方位角度的這個規(guī)律是否具有普遍性,還需要后續(xù)進(jìn)一步驗證.受限于當(dāng)前研究的數(shù)據(jù)量和研究水平,目前只得到了一些初步的規(guī)律,要得到更精確的,更具有普遍性的結(jié)論還需要收集更多的震例數(shù)據(jù),進(jìn)行更完善的實驗,進(jìn)一步探索與震前異常相關(guān)的參數(shù),及各參數(shù)之間的關(guān)聯(lián)性.