鄒斌,郭金運,常曉濤,朱廣彬,李武東

(1.山東科技大學 測繪科學與工程學院,青島 266590;2.國家測繪地理信息局 衛(wèi)星測繪應(yīng)用中心,北京 100048)

?

基于主成分分析與滑動四分位法的震前TEC異常探測對比分析

鄒斌1,2,郭金運1,常曉濤2,朱廣彬2,李武東1,2

(1.山東科技大學 測繪科學與工程學院,青島 266590;2.國家測繪地理信息局 衛(wèi)星測繪應(yīng)用中心,北京 100048)

通過利用IGS提供的全球電離層格網(wǎng)數(shù)據(jù)對四次震級Mw>7.0的地震(汶川大地震,玉樹地震,東日本大地震,于田地震)分別利用主成分分析(PCA)和滑動四分位法進行震前的TEC異常探測。結(jié)果顯示,兩種方法都能有效探測到震前TEC異常,且異常有穩(wěn)定的時間尺度。兩種方法相比,PCA方法更優(yōu),其探測到的異常比較直觀、簡明,而滑動四分位法則探測到的異常天數(shù)較多,但兩者總體趨勢一致,且兩種方法得到的結(jié)果與其他學者分析得到的結(jié)果也一致。從探測的效果上來看PCA不易受太陽地磁活動影響,這也表明PCA的確可用于震前TEC異常分析,分析時將兩種方法有效結(jié)合,相互補充,能對震前的電離層異常分析起到更好的效果。

地震;TEC;主成分分析;滑動四分位法;電離層異常

0　引　言

地震臨震前后電離層擾動以及地震電離層耦合機制一直是地震前兆研究的熱點問題[1],許多學者通過大量的地震——電離層異?，F(xiàn)象的實驗研究,發(fā)現(xiàn)地震孕育過程中,在巖石圈、大氣圈以及電離層確實存在電磁異?，F(xiàn)象,且電離層前兆的出現(xiàn)有穩(wěn)定的時間尺度,這一突出特點將使它在短臨地震預(yù)報中更具實用價值[2]。近年來,由于GPS技術(shù)的快速發(fā)展,使得電離層探測變得更為簡單,很多的電離層探測方法都是基于GPS觀測數(shù)據(jù)進行的。而典型的傳統(tǒng)分析方法是滑動時窗法,它是對TEC時間序列做統(tǒng)計分析進行異常檢查,進而判斷異常,許多學者都采用此方法來探測震前TEC異常情況[3-6]。隨著研究工作的不斷深入,一些學者也在不斷嘗試新的異常探測方法,張小紅等[7]提出一種利用時間序列法(ARIMA模型)進行震前電離層異常探測的新方法,取得很好的效果;Lin[8]利用主成分分析(PCA)對臺灣2002-2003年發(fā)生的12次M≥5.0地震進行異常分析,發(fā)現(xiàn)PCA可用于異常探測且異常主要集中在震前5天;隨后Lin[9]又利用PCA結(jié)合圖像處理技術(shù)分析了2008年汶川地震,對5月4日、8日、9日的全球電離層格網(wǎng)數(shù)據(jù)(GIM)進行主成分分析,結(jié)果顯示PCA能夠有效探測到8日和9日的TEC異常信號并能確定異常的大致位置。在本文研究中,將采用滑動四分位法和PCA對2008-2015年的四次Mw>7.0地震進行對比分析,以研究震前TEC擾動與地震之間的關(guān)系。

1　數(shù)據(jù)及異常檢測方法

1.1數(shù)據(jù)

本文分析的四次地震的基本信息如表1所示,數(shù)據(jù)信息來源于中國地震局官方網(wǎng)站(http://www.cea.gov.cn/),表1中所采用的時間系統(tǒng)為世界協(xié)調(diào)時(UTC)。

表1　2008-2015年4次地震基本信息

分析電離層異常時,采用的電離層數(shù)據(jù)來源于IGS提供的全球電離層格網(wǎng)數(shù)據(jù)(ftp://cddis.gsfc.nasa.gov/gnss/products/ionex/),該數(shù)據(jù)的時間分辨率為2 h,空間分辨率為5°(lon)×2.5°(lat),使用時將其線性內(nèi)插到1 h間隔,通過對震中周圍格網(wǎng)點TEC值進行內(nèi)插得到震中的TEC時間序列[10],本文分析時提取了震中地區(qū)震前15天及震后3天的TEC時間序列。

此外,由于太陽活動和地磁場的變化都會導(dǎo)致電離層的異常變化,因此分析電離層異常時需結(jié)合該時期的空間環(huán)境狀況綜合分析該段時期的電離層異常變化原因。為此,提取了震前15天及震后3天的太陽風速度Vsw、太陽射電量F10.7指數(shù)、地磁Dst指數(shù)、地磁Kp指數(shù)的時間序列。太陽地磁數(shù)據(jù)來源于美國航空航天局(NASA)的OMNIWeb公開數(shù)據(jù)(http://omniweb.gsfc.nasa.gov/),采用數(shù)據(jù)的時間分辨率為1 h.各指數(shù)的評判標準如下:太陽風速度在地球附近的平均值為450 km/s左右;F10.7根據(jù)其大小將太陽活動劃分為三個等級,即F10.7>150SFU,100SFU

1.2主成分分析法

主成分分析法(PCA)又叫KLT被廣泛應(yīng)用于數(shù)據(jù)處理與分析。其主要目的是通過分離不相關(guān)成分提取主要信息并進行數(shù)據(jù)降維[15-16]。

設(shè)數(shù)據(jù)矩陣為

(1)

J(u)=pTp=uTAATu.

(2)

通過拉格朗日乘數(shù)法引入一個拉格朗日乘子,得:

(3)

令

(4)

則J(u)的特征值為λ1≥λ2≥…≥λm,相應(yīng)的特征向量為u1,u2,…,um,最大特征值λ1也即是主特征值,代表主特征信號。

利用PCA分析TEC異常時,每天的數(shù)據(jù)維數(shù)為m=1,n=24,即將經(jīng)時間插值后每天1 h間隔的TEC數(shù)據(jù)代入矩陣A,計算時先將其中心化,然后代入式(4)則可計算出一個主特征值,代表當天的電離層TEC主信號,最后將每天的主特征值組成一個時間序列進行TEC異常分析。

1.3滑動四分位法

四分位距 (IQR)是一種穩(wěn)健統(tǒng)計技術(shù)中用于表示數(shù)據(jù)離散度的一個量,常用來檢查數(shù)據(jù)的異常情況[5-6]。本文考慮到TEC變化具有季節(jié)效應(yīng),選擇窗口時間長度不宜過長,因此在本文中選擇27d[17]作為時間窗口,在分析每天的數(shù)據(jù)時取其前27天同一時刻的TEC組成一組時間序列,并按照從小到大順序?qū)ζ溥M行排列為x1,x2,…,x27則:

(5)

于是四分位距IQR=Q3-Q1,按照標準化正態(tài)分布進行統(tǒng)計計算四分位距與標準差的比值為IQR=1.34σ.本文定義Q2±1.5IQR作為TEC異常檢驗上下邊界的判斷閥值,如果TEC值沒有超出上下界則ΔTEC=0,超出上界ΔTEC>0,表示正異常,低于下界值ΔTEC<0,表示負異常。

2　分析與討論

采用上述兩種異常探測方法對表1中列舉的四次大地震進行TEC異常探測,所得的結(jié)果如圖2、3、5、7所示。從圖中可以看出兩種方法均能探測到震前TEC有明顯的異常,對于滑動四分位法來說其探測出的異常比較容易看出,異常有正有負,只是異常幅度有所差異,而對于PCA來說,其異常表現(xiàn)為在短時間內(nèi)出現(xiàn)顯著的特征值極大值,且該極大值一般就是該時間序列的最大值。

2.1汶川地震

圖1示出了對汶川地震利用滑動四分位法和PCA探測后的結(jié)果。從圖中可以看出在整個時間序列內(nèi)TEC異常正負交替,出現(xiàn)異常的天數(shù)共有4天,分別是震前2天、3天、8天、13天,其中震前3天出現(xiàn)的正異常幅度最大,接近3TECU; 從圖1中PCA得出的結(jié)果,可以看出在整個時間序列內(nèi)主特征值變化比較平穩(wěn),只有在震前3天、震前9天和震前14天有過明顯變化,其中兩個明顯的峰值點分別是在震前3天和震前9天,說明這兩天TEC出現(xiàn)異常。

圖1　2008/05/12汶川地震TEC時間序列及滑動四分位法和PCA結(jié)果(a)震中TEC值及其上下界; (b)滑動四分位結(jié)果; (c)主特征值(圖中豎直線表示地震發(fā)生當天，圖中的特征值序列是歸1化后結(jié)果)

相對于滑動四分位法來說,PCA探測到的異常天數(shù)比之較少,而且兩者并沒有完全對應(yīng),在圖1中滑動四分位法檢測到震前8、13天均出現(xiàn)小的異常,而從PCA得到的結(jié)果來看這兩天的主特征值都比較小,表明這幾天TEC沒有異常;但兩種方法共同探測到震前3天TEC出現(xiàn)異常,其中震前3天的異常在滑動四分位法探測下異常幅度最大,在PCA探測下主特征值大小也僅次于震前9天,因此該天異常最特殊,通過查詢這兩天的太陽地磁活動指數(shù)(Vsw,F10.7,Dst,Kp),結(jié)果顯示這兩天太陽地磁活動都很平靜,將兩種方法綜合對比分析后可進一步認為震前3天的異常與隨后的地震有關(guān),這與余濤、毛田等[18]通過分析汶川震前的主要參量(TEC,foF2)變化得出震前3天的異常與時間太陽地磁活動時間情況一致。

2.2玉樹地震

圖2示出了玉樹地震的探測結(jié)果,從圖2中可看出玉樹地震在滑動四分位法探測下TEC出現(xiàn)異常的天數(shù)是3天,分別是震前6、8、12天,其中異常較顯著的是震前7天和震前8天。圖2也示出了PCA探測結(jié)果,在整個時間序列內(nèi)除了地震當天出現(xiàn)了一個小峰值以及震前12天出現(xiàn)的一個異常顯著的峰值外,其余時間主特征值都很小很平穩(wěn),結(jié)合滑動四分位的探測結(jié)果,兩者共同探測到震前異常只有震前12天,因此震前12天的TEC異常最為特殊。

圖2　2010/04/13玉樹地震TEC時間序列及滑動四分位法和PCA結(jié)果(a)震中TEC值及其上下界; (b)滑動四分位結(jié)果; (c)主特征值;(圖中豎直線表示地震發(fā)生當天，圖中的特征值序列是歸1化后結(jié)果)

圖3　2010/04/13玉樹地震太陽地磁活動時間序列圖(a)太陽網(wǎng)速度Vsw; (b)太陽射點亮F10.7； (c)Dst指數(shù)；(d)Kp指數(shù)(堅線表示地震發(fā)生時刻，橫線是各指數(shù)設(shè)定的閾值)

為了進一步分析震前TEC異常原因,圖3示出了震前15天及震后3天太陽地磁活動時間序列,從圖中可以看出震前7、8天除了F10.7指數(shù)在設(shè)定閥值以內(nèi)外,其余各項指數(shù)均超標:震前8天的太陽風速度最大值接近900km/s, Kp指數(shù)最大值也接近8,震前7天的Dst指數(shù)最低時更是接近-90nT,說明這三天的太陽地磁活動非常強烈,因此圖2用滑動四分位法探測出這幾天的TEC異常與太陽地磁活動干擾密切相關(guān)。

對比圖2中PCA的探測結(jié)果,可以發(fā)現(xiàn)這兩天的主特征值卻很小很穩(wěn)定,因此不難發(fā)現(xiàn)利用PCA探測TEC異常是不受太陽地磁活動干擾的,而滑動四分位法卻受其影響。

再看震前12天太陽風速度在450km/s以下,F10.7低于100SFU,Dst≥-30nT,Kp指數(shù)也低于4,可見該天太陽地磁活動很平靜,但兩種方法都檢測到該天有異常,且PCA檢測時該天的異常最顯著,通過兩種方法對比分析后可認為震前12天的TEC異常是此次地震的前兆,這與劉軍[20]和湯俊[21]用的滑動標準差探測結(jié)果一致。

2.3東日本大地震

在滑動四分位法探測下日本311地震震前電離層異常主要集中出現(xiàn)在震前1～8天,如圖4所示。從圖中可以看出震前3至8天異常最為密集,且震前3天異常幅度最大,接近8TECU,但這幾天的異常過于頻繁,前面也提到滑動四分位法易受太陽地磁活動影響,為此特意查詢了這段時間的太陽地磁活動的情況。

由查詢結(jié)果可知，震前4到9天太陽風速度Vsw和F10.7指數(shù)均超過閾值,可見這幾天太陽活動劇烈,Dst指數(shù)除了在震前第9天的時候低于-30nT外,其余時段還是處于比較穩(wěn)定的狀態(tài),不過Kp指數(shù)一直很接近4,說明這幾天還是有弱的地磁擾動,因此由滑動四分位法探測到的震前4～9天的TEC異常與太陽活動密切。

而震前3天除F10.7指數(shù)在150SFU左右,其余指數(shù)也均在正常范圍內(nèi),該天的異常比較特殊。

圖4中也示出了PCA探測結(jié)果,從圖中可以清晰地看到在震前3天有一個顯著的主特征值極值,說明該天檢測到異常,用PCA探測到震前4、7、9、10天這幾天的主特征值比較小,說明并沒有異常,但用滑動四分位法卻探出這幾天TEC存在異常,而這些天正是太陽地磁活動比較劇烈的時候,這也進一步說明利用PCA分析震前TEC異常不易受太陽地磁活動影響,通過兩種方法對比分析后,可將兩者共同探測到的震前3天的異常作為此次地震的前兆,這與姚宜斌[2]和馬一方[22]等將震前3天的異常作為此次地震前兆的結(jié)果也是一致的。

圖4　2011/03/11日本311大地震TEC時間序列及滑動四分位法和PCA結(jié)果 (a)震中TEC值及其上下界； (b)滑動四分位結(jié)果; (c)主特征值(圖中豎直線表示地震發(fā)生當天，圖中的特征值序列是歸1化后結(jié)果)

2.4于田地震

圖5示出了于田地震的TEC探測結(jié)果,由圖5中可以看出,在滑動四分位法探測下TEC異常主要集中出現(xiàn)在震前8、10天,其中震前10天的異常最為特殊,其最大異常幅度接近16TECU.圖5中也示出了PCA探測結(jié)果,從圖中可清晰地看到震前10天出現(xiàn)一個顯著的特征值極值,另外震前8天與地震當天的主特征值也比較明顯,說明這幾天的TEC有異常。通過查詢這幾天的太陽地磁活動情況,發(fā)現(xiàn)震前2天和4天的太陽地磁活動比較活躍,因此由滑動四分位法探測到的這兩天的小異?？烧J為是太陽地磁活動造成的,再看這兩天的主特征值都比較平穩(wěn),說明PCA并沒有探測到這兩天有異常,這再次驗證了PCA分析震前TEC異常時不受太陽地磁活動影響。而震前8天、10天這兩天除F10.7超出100SFU外,其余各指數(shù)(包括太陽風速度Vsw,Dst,Kp)均在設(shè)定閥值內(nèi),說明太陽地磁活動并不活躍,并且兩種方法都探測到這兩天有異常,但震前10天在兩種方法中都是異常最顯著的,因此可認為與隨后的地震有關(guān),這與李旺[23]等分析得出的結(jié)果也一致。

圖5　2014/02/12新疆于田地震TEC時間序列及滑動四分位法和PCA結(jié)果(a)震中TEC值及其上下界; (b)滑動四分位結(jié)果; (c)主特征值(圖中豎直線表示地震發(fā)生當天，圖中的特征值序列是歸1化后結(jié)果)

3　結(jié)束語

通過對2008-2015年發(fā)生的四次大地震分別利用滑動四分位法和PCA進行震前TEC異常探測,發(fā)現(xiàn)兩種方法均能有效探測到震前的TEC異常信息;從探測結(jié)果來看,滑動四分位法探測到的異常天數(shù)比PCA多,這可能是因為PCA分析異常時不易受太陽地磁活動影響,而滑動四分位法易受其影響;而從異常幅度上來看,PCA探測到的異常更為明顯,一般就是出現(xiàn)峰值時刻且該峰值一般是整個時間序列內(nèi)的最大值,而滑動四分位法往往因其它各種因素的影響使得其探測到的異常幅度不一定是最大的;兩種方法相比,PCA檢測的結(jié)果更清晰明了,而用滑動四分位法得到的結(jié)果較多,不容易剔除無關(guān)因素,無法快速確定具體哪天的TEC異常與地震相關(guān);但將兩種方法有效結(jié)合,通過對比分析并結(jié)合太陽地磁活動排除干擾因素,將這樣得到的結(jié)果作為地震前兆信息會更加可靠更具說服力。

[1] 蔡軍濤,趙國澤,詹艷, 等.地震期間電離層擾動現(xiàn)象研究[J]. 地球物理學進展, 2007, 22(3):695-701.

[2] 姚宜斌,陳鵬,吳寒, 等. 2011年3月11日日本地震震前電離層異常變化分析[J]. 科學通報, 2012, 57(5): 355-365.

[3] 姚璐,申旭輝,張學民. 玉樹MS7. 1 地震前電離層異常擾動分析[J]. 地震, 2014, 34(3):74-85.

[4] 張學民,劉靜,趙必強, 等. 玉樹地震前的電離層異?，F(xiàn)象分析[J]. 空間科學學報, 2014, 34(6): 822-829.

[5]GUOJ,LIW,LIUX, et al.OnTECanomaliesasprecursorbeforeMW8.6SumatraearthquakeandMW6.7MexicoearthquakeonApril11, 2012 [J].GeosciencesJournal, 2015, 19(4):721-730.

[6]GUOJ,LIW,YUH, et al.ImpendingionosphericanomalyprecedingtheIquiqueMw8. 2earthquakeinChileon2014April1 [J].GeophysicalJournalInternational, 2015, 203(3):1461-1470.

[7] 張小紅,任曉東,吳風波,等. 震前電離層TEC異常探測新方法[J]. 地球物理學報, 2013 (2): 441-449.

[8]LINJW.Ionospherictotalelectroncontent(TEC)anomaliesassociatedwithearthquakesthroughKarhunen-LoéveTransform(KLT) [J].Terrestrial,AtmosphericandOceanicSciences, 2010, 21(2):253-265.

[9]LINJW.Two-dimensionalionospherictotalelectroncontentmap(TEC)seismo-ionosphericanomaliesthroughimageprocessingusingprincipalcomponentanalysis[J].AdvancesinSpaceResearch, 2010, 45(11): 1301-1310.

[10]李征航,張小紅. 衛(wèi)星導(dǎo)航定位新技術(shù)及高精度數(shù)據(jù)處理方法[M]. 武漢大學出版社, 2009.

[11]FEJERBG,DEPAULAER,GONZLEZSA,etal.AverageverticalandzonalF-regionplasmadriftsoverJicamarca[J].JournalofGeophysicalResearch, 1991, 96(A8):13901-13906.

[12]FEJERBG,PAULAERD,BATISTAIS, et al.EquatorialFregionverticalplasmadriftsduringsolarmaxima[J].JournalofGeophysicalResearch, 1989, 94(A9):12049-12054.

[13]沈曉飛,倪彬彬, 顧旭東, 等. 第23太陽活動周期太陽風參數(shù)及地磁指數(shù)的統(tǒng)計分析[J]. 地球物理學報, 2015, 58(2):362-370.

[14]LINJW.PotentialreasonsforionosphericanomaliesimmediatelypriortoChina’sWenchuanearthquakeon12May2008detectedbynonlinearprincipalcomponentanalysis[J].InternationalJournalofAppliedEarthObservationandGeoinformation, 2012, 14(1):178-191.

[15]LINJW.Principalcomponentanalysismethodinthedetectionoftotalelectroncontentanomaliesinthe24hrspriortolargeearthquakes[J].ArabianJournalofGeosciences, 2011, 4(3-4):575-593.

[16]GMEZLONDOE,CASTILLOLPEZL,THASDESOUZAK.UsingtheKarhunen-Loóvetransformtosuppressgroundrollinseismicdata[J].EarthSciencesResearchJournal, 2005, 9(2):139-147.

[17]PANCHEVAD,LASTOVICKAJ.Solarormeteorologicalcontroloflowerionosphericfluctuations(2-15and27days)inmiddlelatitudes[J].HandbookforMAP, 1989(29):210-214.

[18]余濤,毛田,王云岡, 等. 汶川特大地震前電離層主要參量變化[J]. 科學通報, 2009 (4):493-499.

[19]LINJW.Ionosphericprecursorof2008ChinaWenchuanearthquakeusingtwo-dimensionalprincipalcomponentanalysis[J].HKIETransactions, 2014, 21(3):192-194.

[20]劉軍,柴洪洲,劉長建,等. 基于CODEGIM的震前電離層TEC異常分析[J]. 大地測量與地球動力學, 2011, 31(6):39-42.

[21]湯俊,姚宜斌,陳鵬, 等. 利用PCA分析震前電離層TEC異常[J]. 大地測量與地球動力學, 2013, 33(4):22-25.

[22]馬一方,匡翠林,周曉慧. 基于GPS數(shù)據(jù)的震前電離層異常分析[J]. 測繪通報, 2015, 4:001.

[23]LIW,GUOJY,YUXM, et al.AnalysisofionosphericanomalyprecedingtheMw7.3Yutianearthquake[J].GeodesyandGeodynamics, 2014, 5(2):54-60.

Analysis of TEC Anomalies before Earthquake Based on Principal Component Analysis and the Sliding Inter Quartile Range Method

ZOU Bin1,2,GUO Jinyun1,CHANG Xiaotao2,ZHU Guangbin2,LI Wudong1,2

(1.CollegeofGeodesyandGeomatics,ShandongUniversityofScienceandTechnology,Qingdao266590,China;2.SatelliteSurveyingandMappingApplicationCenter,NationalAdministrationofSurveying,MappingandGeo-information,Beijing100048,China)

This paper using the global ionospheric map data from IGS and then using principal component analysis method and the sliding inter quartile range method to investigate precursors for 4 earthquakes of Richter magnitude scale Mw≥7.0(Wenchuan earthquake, Yushu earthquake, Japan 311 earthquake, Yutian earthquake). Results shows that both methods can effectively detect pre-earthquake TEC anomaly and the anomaly has stable time scale. Compared two methods, PCA is better, it detects the anomaly more intuitive and concise, but the sliding inter quartile range method detects more abnormal days than PCA. However, the results of two methods has a common trend, and the results of two methods is consistent with the results from other scholars. In addition, PCA is not affected by solar geomagnetic activity, this finding verify that PCA can be used in pre-earthquake ionospheric anomaly analysis. So it obvious that combine two methods effectively when analysis the pre-earthquake ionospheric anomaly can have a better effect.

Earthquake; Total Electron Content; principal component analysis; sliding inter quartile range method; ionospheric anomaly

10.13442/j.gnss.1008-9268.2016.04.014

2016-03-09

國家自然科學基金(批準號:41374009,41204007); 山東省自然科學基金(批準號:ZR2013DM009); 國家科技基礎(chǔ)性工作專項(編號:2015FY310200); 國家重點基礎(chǔ)研究發(fā)展計劃(編號:2013CB733302); 測繪地理信息公益性行業(yè)科研專項項目(編號:201512012)

P228.4

A

1008-9268(2016)04-0063-08

鄒斌(1992-),男,江西撫州人,碩士生,主要從事空間大地測量研究。

郭金運(1969-),男,山東巨野人,博士、教授、博導(dǎo),主要從事空間大地測量、海洋大地測量和物理大地測量等研究。

常曉濤(1972-),男,河南人,研究員,主要從事衛(wèi)星大地測量方面的研究。

聯(lián)系人： 鄒斌 E-mail: binchow77tt@163.com