基于射出長波輻射數據的震前異常時空分析法①

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

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

2(福建農林大學 計算機與信息學院,福州 350002)

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

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

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

1 相關數據

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

表1 震例相關信息

2 研究方法

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

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

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

其中,num表示從可獲取數據的起始年份到地震發生年的總年數.

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

得到具有時空相關性的新增量序列,記為Tn,則有Tn=T∪ti,其中T表示歷史數據點,ti表示當前數據點.

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

此處,i∈[?ws,ws],q=0.5,而是數學中的組合公式.

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

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

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

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

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

3 結果與分析

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

3.1 汶川地震

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

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

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

圖2 汶川地震時空分析圖

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

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

3.2 日本地震

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

圖4 日本地震時空分析

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

3.3 蘆山地震

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

圖6 蘆山時空分析圖

而從圖7可以看出,蘆山地震發生前4 個月已經出現明顯,并在地震發生前一個月內CD值再次急劇增大,并在震前半個月達到最大值.第4,7 和9 單元格震前均無明顯變化趨勢.

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

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

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

3.4 對比實驗

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

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

3.5 參數討論

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

圖9 不同h 對日本地震結果的影響

圖10 不同ε 對日本地震結果的影響

圖11 不同ws 對日本地震結果的影響

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

4 結論

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

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