基于殘差修正的ARMA模型探測(cè)門(mén)源MS6.4地震前電離層TEC異常

李 磊 張 寧 尹淑慧 王軼卓

1 大連海事大學(xué)理學(xué)院，大連市凌海路1號(hào)，116026

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

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

(1)

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

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

(2)

式中，TECpre、TECIGS分別為T(mén)EC預(yù)報(bào)值和IGS發(fā)布值。

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

表1 3種方法預(yù)測(cè)TEC值的誤差統(tǒng)計(jì)

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

3 門(mén)源MS6.4地震電離層TEC異常探測(cè)

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

3.1 ARMA-RC模型計(jì)算TEC參考背景值

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

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

3.2 TEC異常限差判定策略

(3)

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

基于門(mén)源縣TEC數(shù)據(jù)，初步發(fā)現(xiàn)10 d共17個(gè)時(shí)刻出現(xiàn)TEC異常。為排除太陽(yáng)、地磁異常活動(dòng)等因素對(duì)TEC數(shù)據(jù)的影響，對(duì)地震前14 d各時(shí)段地磁活動(dòng)的Kp指數(shù)、Dst指數(shù)和太陽(yáng)輻射流量F10.7進(jìn)行分析。從圖2可以看出，在探測(cè)時(shí)段內(nèi)，F(xiàn)10.7遠(yuǎn)小于150 SFU，指數(shù)變化相對(duì)穩(wěn)定，可以排除太陽(yáng)活動(dòng)引起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地磁活動(dòng)Kp指數(shù)均超過(guò)3，表明存在一定的地磁干擾；除01-21外，其他日期Dst指數(shù)絕對(duì)值均在30 nT以內(nèi)。因此，剔除上述日期內(nèi)的TEC異常后，其他異常極可能由此次地震引起。

3.4 探測(cè)結(jié)果及分析

3.4.1 基于IGS數(shù)據(jù)作為T(mén)EC真值的結(jié)果分析

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

表2 門(mén)源地震震中上空TEC異常擾動(dòng)數(shù)據(jù)

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

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

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

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

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

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

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

圖4 ARMA-RC模型探測(cè)區(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ū)域，異常峰值點(diǎn)均位于震中偏南方向，同時(shí)震中位置并不是異常區(qū)域的中心點(diǎn)：UT 00:00發(fā)生負(fù)異常時(shí)，震中位于負(fù)異常區(qū)域的偏東方向；UT 06:00發(fā)生正異常時(shí)，震中位于正異常區(qū)域的偏南方向。上述分析特征與前人的研究結(jié)論基本一致[4-6, 11]，表明2016-01-16 UT 00:00和UT 05:00的電離層異?，F(xiàn)象與本次地震密切相關(guān)。

4 結(jié) 語(yǔ)

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

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

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