基于殘差修正的ARMA模型探測門源MS6.4地震前電離層TEC異常

李 磊 張 寧 尹淑慧 王軼卓

1 大連海事大學理學院，大連市凌海路1號，116026

震前電離層擾動與電離層耦合機制是地震前兆研究的熱點之一[1-3]。研究表明，在5級以上地震發(fā)生前，震中及其附近地區(qū)的電離層TEC會存在異常擾動的現(xiàn)象[2-6]。傳統(tǒng)的探測震前TEC異常的方法主要包括中位數(shù)法、四分位距法和滑動時窗法等[6]，雖然這些方法已成功分析了許多震前TEC異常現(xiàn)象，但其計算TEC參考背景值的精度較低，同時又缺乏統(tǒng)一的限差判決門限，因此會影響電離層異常擾動事件的準確判斷[5]。基于上述原因，許多新的震前TEC異常探測方法[6-8]被提出，其中基于時間序列分析理論的方法受到大量關注。但該方法主要用于7級以上的地震震例分析，對5、6級地震震例的研究很少。本文利用ARMA模型對TEC預測數(shù)據(jù)進行殘差修正，獲得TEC參考背景值，提出相對誤差和絕對誤差相互制約的限差判決門限，并以2016-01-21青海門源MS6.4地震為例，研究震前14 d和地震當天震中及周圍區(qū)域的TEC異常。

1 殘差修正的ARMA模型(ARMA-RC)

(1)

2 ARMA-RC模型與傳統(tǒng)方法預測精度比較

基于IGS提供的2014-11-01～12-29網(wǎng)格點(0°N, 90°E)TEC數(shù)據(jù)，分別采用ARMA-RC模型、滑動時窗法(時窗長度為10 d)和四分位距法預測2014-12-25～29共120 h的TEC小時數(shù)據(jù)。引入相對誤差(Ur)和絕對誤差(Δ)對結果進行分析：

(2)

式中，TECpre、TECIGS分別為TEC預報值和IGS發(fā)布值。

由表1可知，在120 h預報時間內，無論是從|Δ|<2 TECu的數(shù)據(jù)量占全部數(shù)據(jù)的比例，還是|Δ|<1 TECu所表示的較高精度來看，ARMA-RC模型的預測精度都最高。圖1橫軸為從預報開始時刻起算的小時數(shù)，可以看出，3種方法中，ARMA-RC模型預測TEC值的絕對誤差和相對誤差波動范圍都最小，說明其預測TEC的精度高于傳統(tǒng)方法。

表1 3種方法預測TEC值的誤差統(tǒng)計

圖1 ARMA-RC模型與傳統(tǒng)方法預測TEC的誤差對比Fig.1 Comparison of errors of predicted TEC between ARMA-RC model and traditional methods

3 門源MS6.4地震電離層TEC異常探測

門源MS6.4地震發(fā)生于2016-01-21 UT 01:13，震中位置為37.7°N、101.6°E。基于IGS發(fā)布的時間分辨率為1 h的TEC數(shù)據(jù)，通過雙線性插值求出震中位置TEC值，數(shù)據(jù)時段為2015-11-14～2016-01-21。

3.1 ARMA-RC模型計算TEC參考背景值

將震前數(shù)據(jù)分成3組：2015-11-14～2016-01-06、2015-11-19～2016-01-11、2015-11-24～2016-01-16，分別作為訓練數(shù)據(jù)用于參數(shù)估計及模型建立。將震前14 d和地震當天分為3個時段：2016-01-07～11、2015-11-12～2016-01-16、2015-11-17～2016-01-21，其相應數(shù)據(jù)作為觀測值。

基于訓練數(shù)據(jù)，參照文獻[7]的計算流程，利用ARMA-RC模型，分別對未來5 d的TEC小時數(shù)據(jù)進行滑動遞推預測，得到2016-01-07～11、2015-11-12～2016-01-16、2015-11-17～2016-01-21時段TEC參考背景值。對新加入的訓練數(shù)據(jù)進行異常剔除，處理方法為利用相鄰內插法計算正常值并將其替換[6]。

3.2 TEC異常限差判定策略

(3)

3.3 震前太陽和地磁環(huán)境

基于門源縣TEC數(shù)據(jù)，初步發(fā)現(xiàn)10 d共17個時刻出現(xiàn)TEC異常。為排除太陽、地磁異常活動等因素對TEC數(shù)據(jù)的影響，對地震前14 d各時段地磁活動的Kp指數(shù)、Dst指數(shù)和太陽輻射流量F10.7進行分析。從圖2可以看出，在探測時段內，F(xiàn)10.7遠小于150 SFU，指數(shù)變化相對穩(wěn)定，可以排除太陽活動引起TEC異常的可能性；01-07、

圖2 震前(01-07～01-21)Dst、Kp和F10.7指數(shù)變化Fig.2 Variation of Dst, Kp and F10.7 before Menyuan earthquake (from January 7 to 21)

01-11～13、01-19～21地磁活動Kp指數(shù)均超過3，表明存在一定的地磁干擾；除01-21外，其他日期Dst指數(shù)絕對值均在30 nT以內。因此，剔除上述日期內的TEC異常后，其他異常極可能由此次地震引起。

3.4 探測結果及分析

3.4.1 基于IGS數(shù)據(jù)作為TEC真值的結果分析

表2為2016-01-07～21時段TEC異常的相關數(shù)據(jù)，其中，ΔTEC=TECIGS-TECARMA，ΔTEC>0表示正異常，ΔTEC<0表示負異常；正負異常持續(xù)時間(ΔT)為ΔTEC連續(xù)同號的小時數(shù)。圖3(a)、3(b)分別為利用滑動時窗法和ARMA-RC模型探測的TEC異常結果，數(shù)字1～7表示異常數(shù)據(jù)點，EQ為地震發(fā)生時刻位置，綠色區(qū)域表示TEC參考背景值上下限值區(qū)間，黑色曲線為實際探測的TEC值。

表2 門源地震震中上空TEC異常擾動數(shù)據(jù)

圖3 門源地震前(2016-01-07～21)TEC異常擾動結果Fig.3 The results of TEC anomalies before Menyuan earthquake (from January 7 to 21,2016)

由圖3可知，2種方法探測的結果存在明顯不同。滑動時窗法由于預報的參考背景值精度較低，其發(fā)生異常的天數(shù)明顯多于ARMA-RC模型，且異常一旦出現(xiàn)，隨后幾天均會出現(xiàn)異常，同時正負異常無明顯分布規(guī)律，該結論與前人的討論結果相符[1,6]。由于滑動時窗法探測結果的準確性存在一定問題，因此本文只對ARMA-RC模型的探測結果作進一步分析。

由表2及圖3(b)可知，震前14 d內共發(fā)生4次正異常和3次負異常，正負異常次數(shù)基本相同；正異常最長持續(xù)時間為3 h，在震前第5天UT 05:00達到此次正異常峰值(位置7)，該點也為震前正異常的最大值時刻；負異常最長持續(xù)時間為3 h，在震前第11天UT 09:00達到此次負異常峰值(位置4)，但負異常最大值時刻并未與該點重合，而是發(fā)生在震前第5天UT 00:00(位置6)。對于此次地震而言，正負異常的最大值發(fā)生處均為最臨近地震的2次異常，并且負異常發(fā)生時間早于正異常。

3.4.2 考慮IGS數(shù)據(jù)解算精度的結果分析

IGS發(fā)布的電離層TEC數(shù)據(jù)的解算精度為2～8 TECu，為進一步驗證表2中TEC異常結果的可靠性，對上述TEC異常作進一步篩選。綜合考慮IGS數(shù)據(jù)的解算精度，制定如下異常篩選方案：1)將表2中|ΔTEC|與IGS發(fā)布的RMS值進行對比，剔除|ΔTEC|遠小于RMS的數(shù)據(jù)點；2)通過區(qū)域二維電離層地圖進行異常分析。根據(jù)方案1，表2中序號6、7的數(shù)據(jù)點可作為進一步待分析的異常候選點，其他數(shù)據(jù)點受IGS數(shù)據(jù)解算精度所限而無法有效判斷，暫不進行分析。

為進一步分析震前TEC異常的二維空間分布情況，利用ARMA-RC模型對0°～60°N、50°～150°E范圍內的GIM電離層地圖進行異常探測。探測數(shù)據(jù)時段與震中相同，時間分辨率為1 h，空間分辨率為經(jīng)度5°、緯度2.5°，共計525個探測點。圖4為2016-01-15 UT 23:00～2016-01-16 UT 07:00的TEC異常分布圖，圖中五角星為震中位置，圖4(b)、4(g)分別對應表2中序號6、7的震中異常時刻。

從圖4可以看出，UT 23:00震中區(qū)域電離層非常平靜，震中西南方向(15°N，70°E)附近出現(xiàn)明顯負異常，最大異常峰值為-2.9 TECu；隨著時間的推移，異常區(qū)域向震中方向移動并逐漸增大，且呈現(xiàn)集中現(xiàn)象， 在UT 00:00形成覆蓋震中的負異常區(qū)域，震中位于該區(qū)域的偏東方向，最大異常峰值為-2.6 TECu，位于震中以南；隨后異常區(qū)域遠離震中并向東南方向擴散，在UT 01:00形成較大范圍、分散的負異常區(qū)域，最大異常峰值為-3.7 TECu，并在10°N、150°E附近出現(xiàn)明顯正異常區(qū)域，最大異常峰值為4.4 TECu；UT 02:00震中東南方向出現(xiàn)異常峰值為3.8 TECu和-4.6 TECu的2個異常區(qū)域，查詢2016年全球地震動態(tài)可知，2016-01-27加羅林群島西部(8.3°N，137.8°E)、2016-01-19中國臺灣臺東縣海域(22.9°N，121.3°E)分別發(fā)生5.2級和5.8級地震，2次地震震中恰好位于圖4(d)中正、負異常區(qū)域，2次異常可能與這2次地震有關，具體關聯(lián)還需進一步研究；UT 03:00震中附近電離層又趨于平靜；UT 04:00震中西部出現(xiàn)峰值為-2.7 TECu的負異常；UT 05:00形成覆蓋震中的正異常區(qū)域，震中位于該區(qū)域的偏南方向，正異常峰值同樣位于震中以南，最大異常峰值為2.8 TECu；隨后異常區(qū)域向東北方向移動并逐漸減小，震中附近區(qū)域異常再次消失。

圖4 ARMA-RC模型探測區(qū)域電離層異常(2016-01-15 UT 23:00～2016-01-16 UT 07:00)Fig.4 TEC anomalies based on ARMA-RC model(from UT 23:00, January 15,2016 to UT 07:00, January 16,2016)

綜上所述，2016-01-16 UT 00:00和UT 05:00在震中及其附近形成范圍非常集中的電離層異常區(qū)域，異常峰值點均位于震中偏南方向，同時震中位置并不是異常區(qū)域的中心點：UT 00:00發(fā)生負異常時，震中位于負異常區(qū)域的偏東方向；UT 06:00發(fā)生正異常時，震中位于正異常區(qū)域的偏南方向。上述分析特征與前人的研究結論基本一致[4-6, 11]，表明2016-01-16 UT 00:00和UT 05:00的電離層異常現(xiàn)象與本次地震密切相關。

4 結 語

本文利用ARMA模型對TEC預測殘差進行修正，得到TEC參考背景值，并設定相對誤差和絕對誤差相互制約的限差判定策略，探測門源MS6.4地震前TEC異常。結果表明：1)震前第5天存在明顯的TEC異常現(xiàn)象，且異常與太陽活動、地磁活動等影響無關，從而驗證了本方法對震前TEC異常探測的有效性；2)震前TEC正負異常發(fā)生次數(shù)基本相同，并且負異常發(fā)生在正異常之前；3)震中及附近更大區(qū)域內形成范圍非常集中的電離層異常區(qū)域，異常峰值點位于震中以南，震中位置與異常區(qū)域的中心點并不重合，位于負異常區(qū)域偏東方向和正異常區(qū)域偏南方向。

本方法可為震前TEC異常探測中參考背景值的確定和限差判決門限的設定提供重要參考。為了更好地驗證本文方法的有效性，下一步將選取更多、更經(jīng)典的震例進行地震前后異常情況探測。

致謝:感謝IGS、日本京都地磁中心、中國科學院空間中心為本文提供電離層、Kp、Dst和F10.7數(shù)據(jù)。