基于主成分分析与滑动四分位法的震前TEC异常探测对比分析

邹斌,郭金运,常晓涛,朱广彬,李武东

(1.山东科技大学 测绘科学与工程学院,青岛 266590;2.国家测绘地理信息局 卫星测绘应用中心,北京 100048)



基于主成分分析与滑动四分位法的震前TEC异常探测对比分析

邹斌1,2,郭金运1,常晓涛2,朱广彬2,李武东1,2

(1.山东科技大学 测绘科学与工程学院,青岛 266590;2.国家测绘地理信息局 卫星测绘应用中心,北京 100048)

通过利用IGS提供的全球电离层格网数据对四次震级Mw>7.0的地震(汶川大地震,玉树地震,东日本大地震,于田地震)分别利用主成分分析(PCA)和滑动四分位法进行震前的TEC异常探测。结果显示,两种方法都能有效探测到震前TEC异常,且异常有稳定的时间尺度。两种方法相比,PCA方法更优,其探测到的异常比较直观、简明,而滑动四分位法则探测到的异常天数较多,但两者总体趋势一致,且两种方法得到的结果与其他学者分析得到的结果也一致。从探测的效果上来看PCA不易受太阳地磁活动影响,这也表明PCA的确可用于震前TEC异常分析,分析时将两种方法有效结合,相互补充,能对震前的电离层异常分析起到更好的效果。

地震;TEC;主成分分析;滑动四分位法;电离层异常

0　引　言

地震临震前后电离层扰动以及地震电离层耦合机制一直是地震前兆研究的热点问题[1],许多学者通过大量的地震——电离层异常现象的实验研究,发现地震孕育过程中,在岩石圈、大气圈以及电离层确实存在电磁异常现象,且电离层前兆的出现有稳定的时间尺度,这一突出特点将使它在短临地震预报中更具实用价值[2]。近年来,由于GPS技术的快速发展,使得电离层探测变得更为简单,很多的电离层探测方法都是基于GPS观测数据进行的。而典型的传统分析方法是滑动时窗法,它是对TEC时间序列做统计分析进行异常检查,进而判断异常,许多学者都采用此方法来探测震前TEC异常情况[3-6]。随着研究工作的不断深入,一些学者也在不断尝试新的异常探测方法,张小红等[7]提出一种利用时间序列法(ARIMA模型)进行震前电离层异常探测的新方法,取得很好的效果;Lin[8]利用主成分分析(PCA)对台湾2002-2003年发生的12次M≥5.0地震进行异常分析,发现PCA可用于异常探测且异常主要集中在震前5天;随后Lin[9]又利用PCA结合图像处理技术分析了2008年汶川地震,对5月4日、8日、9日的全球电离层格网数据(GIM)进行主成分分析,结果显示PCA能够有效探测到8日和9日的TEC异常信号并能确定异常的大致位置。在本文研究中,将采用滑动四分位法和PCA对2008-2015年的四次Mw>7.0地震进行对比分析,以研究震前TEC扰动与地震之间的关系。

1　数据及异常检测方法

1.1数据

本文分析的四次地震的基本信息如表1所示,数据信息来源于中国地震局官方网站(http://www.cea.gov.cn/),表1中所采用的时间系统为世界协调时(UTC)。

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

分析电离层异常时,采用的电离层数据来源于IGS提供的全球电离层格网数据(ftp://cddis.gsfc.nasa.gov/gnss/products/ionex/),该数据的时间分辨率为2 h,空间分辨率为5°(lon)×2.5°(lat),使用时将其线性内插到1 h间隔,通过对震中周围格网点TEC值进行内插得到震中的TEC时间序列[10],本文分析时提取了震中地区震前15天及震后3天的TEC时间序列。

此外,由于太阳活动和地磁场的变化都会导致电离层的异常变化,因此分析电离层异常时需结合该时期的空间环境状况综合分析该段时期的电离层异常变化原因。为此,提取了震前15天及震后3天的太阳风速度Vsw、太阳射电量F10.7指数、地磁Dst指数、地磁Kp指数的时间序列。太阳地磁数据来源于美国航空航天局(NASA)的OMNIWeb公开数据(http://omniweb.gsfc.nasa.gov/),采用数据的时间分辨率为1 h.各指数的评判标准如下:太阳风速度在地球附近的平均值为450 km/s左右;F10.7根据其大小将太阳活动划分为三个等级,即F10.7>150SFU,100SFU

1.2主成分分析法

主成分分析法(PCA)又叫KLT被广泛应用于数据处理与分析。其主要目的是通过分离不相关成分提取主要信息并进行数据降维[15-16]。

设数据矩阵为

(1)

J(u)=pTp=uTAATu.

(2)

通过拉格朗日乘数法引入一个拉格朗日乘子,得:

(3)

令

(4)

则J(u)的特征值为λ1≥λ2≥…≥λm,相应的特征向量为u1,u2,…,um,最大特征值λ1也即是主特征值,代表主特征信号。

利用PCA分析TEC异常时,每天的数据维数为m=1,n=24,即将经时间插值后每天1 h间隔的TEC数据代入矩阵A,计算时先将其中心化,然后代入式(4)则可计算出一个主特征值,代表当天的电离层TEC主信号,最后将每天的主特征值组成一个时间序列进行TEC异常分析。

1.3滑动四分位法

四分位距 (IQR)是一种稳健统计技术中用于表示数据离散度的一个量,常用来检查数据的异常情况[5-6]。本文考虑到TEC变化具有季节效应,选择窗口时间长度不宜过长,因此在本文中选择27d[17]作为时间窗口,在分析每天的数据时取其前27天同一时刻的TEC组成一组时间序列,并按照从小到大顺序对其进行排列为x1,x2,…,x27则:

(5)

于是四分位距IQR=Q3-Q1,按照标准化正态分布进行统计计算四分位距与标准差的比值为IQR=1.34σ.本文定义Q2±1.5IQR作为TEC异常检验上下边界的判断阀值,如果TEC值没有超出上下界则ΔTEC=0,超出上界ΔTEC>0,表示正异常,低于下界值ΔTEC<0,表示负异常。

2　分析与讨论

采用上述两种异常探测方法对表1中列举的四次大地震进行TEC异常探测,所得的结果如图2、3、5、7所示。从图中可以看出两种方法均能探测到震前TEC有明显的异常,对于滑动四分位法来说其探测出的异常比较容易看出,异常有正有负,只是异常幅度有所差异,而对于PCA来说,其异常表现为在短时间内出现显著的特征值极大值,且该极大值一般就是该时间序列的最大值。

2.1汶川地震

图1示出了对汶川地震利用滑动四分位法和PCA探测后的结果。从图中可以看出在整个时间序列内TEC异常正负交替,出现异常的天数共有4天,分别是震前2天、3天、8天、13天,其中震前3天出现的正异常幅度最大,接近3TECU; 从图1中PCA得出的结果,可以看出在整个时间序列内主特征值变化比较平稳,只有在震前3天、震前9天和震前14天有过明显变化,其中两个明显的峰值点分别是在震前3天和震前9天,说明这两天TEC出现异常。

图1　2008/05/12汶川地震TEC时间序列及滑动四分位法和PCA结果(a)震中TEC值及其上下界; (b)滑动四分位结果; (c)主特征值(图中竖直线表示地震发生当天，图中的特征值序列是归1化后结果)

相对于滑动四分位法来说,PCA探测到的异常天数比之较少,而且两者并没有完全对应,在图1中滑动四分位法检测到震前8、13天均出现小的异常,而从PCA得到的结果来看这两天的主特征值都比较小,表明这几天TEC没有异常;但两种方法共同探测到震前3天TEC出现异常,其中震前3天的异常在滑动四分位法探测下异常幅度最大,在PCA探测下主特征值大小也仅次于震前9天,因此该天异常最特殊,通过查询这两天的太阳地磁活动指数(Vsw,F10.7,Dst,Kp),结果显示这两天太阳地磁活动都很平静,将两种方法综合对比分析后可进一步认为震前3天的异常与随后的地震有关,这与余涛、毛田等[18]通过分析汶川震前的主要参量(TEC,foF2)变化得出震前3天的异常与时间太阳地磁活动时间情况一致。

2.2玉树地震

图2示出了玉树地震的探测结果,从图2中可看出玉树地震在滑动四分位法探测下TEC出现异常的天数是3天,分别是震前6、8、12天,其中异常较显著的是震前7天和震前8天。图2也示出了PCA探测结果,在整个时间序列内除了地震当天出现了一个小峰值以及震前12天出现的一个异常显著的峰值外,其余时间主特征值都很小很平稳,结合滑动四分位的探测结果,两者共同探测到震前异常只有震前12天,因此震前12天的TEC异常最为特殊。

图2　2010/04/13玉树地震TEC时间序列及滑动四分位法和PCA结果(a)震中TEC值及其上下界; (b)滑动四分位结果; (c)主特征值;(图中竖直线表示地震发生当天，图中的特征值序列是归1化后结果)

图3　2010/04/13玉树地震太阳地磁活动时间序列图(a)太阳网速度Vsw; (b)太阳射点亮F10.7； (c)Dst指数；(d)Kp指数(坚线表示地震发生时刻，横线是各指数设定的阈值)

为了进一步分析震前TEC异常原因,图3示出了震前15天及震后3天太阳地磁活动时间序列,从图中可以看出震前7、8天除了F10.7指数在设定阀值以内外,其余各项指数均超标:震前8天的太阳风速度最大值接近900km/s, Kp指数最大值也接近8,震前7天的Dst指数最低时更是接近-90nT,说明这三天的太阳地磁活动非常强烈,因此图2用滑动四分位法探测出这几天的TEC异常与太阳地磁活动干扰密切相关。

对比图2中PCA的探测结果,可以发现这两天的主特征值却很小很稳定,因此不难发现利用PCA探测TEC异常是不受太阳地磁活动干扰的,而滑动四分位法却受其影响。

再看震前12天太阳风速度在450km/s以下,F10.7低于100SFU,Dst≥-30nT,Kp指数也低于4,可见该天太阳地磁活动很平静,但两种方法都检测到该天有异常,且PCA检测时该天的异常最显著,通过两种方法对比分析后可认为震前12天的TEC异常是此次地震的前兆,这与刘军[20]和汤俊[21]用的滑动标准差探测结果一致。

2.3东日本大地震

在滑动四分位法探测下日本311地震震前电离层异常主要集中出现在震前1～8天,如图4所示。从图中可以看出震前3至8天异常最为密集,且震前3天异常幅度最大,接近8TECU,但这几天的异常过于频繁,前面也提到滑动四分位法易受太阳地磁活动影响,为此特意查询了这段时间的太阳地磁活动的情况。

由查询结果可知，震前4到9天太阳风速度Vsw和F10.7指数均超过阈值,可见这几天太阳活动剧烈,Dst指数除了在震前第9天的时候低于-30nT外,其余时段还是处于比较稳定的状态,不过Kp指数一直很接近4,说明这几天还是有弱的地磁扰动,因此由滑动四分位法探测到的震前4～9天的TEC异常与太阳活动密切。

而震前3天除F10.7指数在150SFU左右,其余指数也均在正常范围内,该天的异常比较特殊。

图4中也示出了PCA探测结果,从图中可以清晰地看到在震前3天有一个显著的主特征值极值,说明该天检测到异常,用PCA探测到震前4、7、9、10天这几天的主特征值比较小,说明并没有异常,但用滑动四分位法却探出这几天TEC存在异常,而这些天正是太阳地磁活动比较剧烈的时候,这也进一步说明利用PCA分析震前TEC异常不易受太阳地磁活动影响,通过两种方法对比分析后,可将两者共同探测到的震前3天的异常作为此次地震的前兆,这与姚宜斌[2]和马一方[22]等将震前3天的异常作为此次地震前兆的结果也是一致的。

图4　2011/03/11日本311大地震TEC时间序列及滑动四分位法和PCA结果 (a)震中TEC值及其上下界； (b)滑动四分位结果; (c)主特征值(图中竖直线表示地震发生当天，图中的特征值序列是归1化后结果)

2.4于田地震

图5示出了于田地震的TEC探测结果,由图5中可以看出,在滑动四分位法探测下TEC异常主要集中出现在震前8、10天,其中震前10天的异常最为特殊,其最大异常幅度接近16TECU.图5中也示出了PCA探测结果,从图中可清晰地看到震前10天出现一个显著的特征值极值,另外震前8天与地震当天的主特征值也比较明显,说明这几天的TEC有异常。通过查询这几天的太阳地磁活动情况,发现震前2天和4天的太阳地磁活动比较活跃,因此由滑动四分位法探测到的这两天的小异常可认为是太阳地磁活动造成的,再看这两天的主特征值都比较平稳,说明PCA并没有探测到这两天有异常,这再次验证了PCA分析震前TEC异常时不受太阳地磁活动影响。而震前8天、10天这两天除F10.7超出100SFU外,其余各指数(包括太阳风速度Vsw,Dst,Kp)均在设定阀值内,说明太阳地磁活动并不活跃,并且两种方法都探测到这两天有异常,但震前10天在两种方法中都是异常最显著的,因此可认为与随后的地震有关,这与李旺[23]等分析得出的结果也一致。

图5　2014/02/12新疆于田地震TEC时间序列及滑动四分位法和PCA结果(a)震中TEC值及其上下界; (b)滑动四分位结果; (c)主特征值(图中竖直线表示地震发生当天，图中的特征值序列是归1化后结果)

3　结束语

通过对2008-2015年发生的四次大地震分别利用滑动四分位法和PCA进行震前TEC异常探测,发现两种方法均能有效探测到震前的TEC异常信息;从探测结果来看,滑动四分位法探测到的异常天数比PCA多,这可能是因为PCA分析异常时不易受太阳地磁活动影响,而滑动四分位法易受其影响;而从异常幅度上来看,PCA探测到的异常更为明显,一般就是出现峰值时刻且该峰值一般是整个时间序列内的最大值,而滑动四分位法往往因其它各种因素的影响使得其探测到的异常幅度不一定是最大的;两种方法相比,PCA检测的结果更清晰明了,而用滑动四分位法得到的结果较多,不容易剔除无关因素,无法快速确定具体哪天的TEC异常与地震相关;但将两种方法有效结合,通过对比分析并结合太阳地磁活动排除干扰因素,将这样得到的结果作为地震前兆信息会更加可靠更具说服力。

[1] 蔡军涛,赵国泽,詹艳, 等.地震期间电离层扰动现象研究[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] 张学民,刘静,赵必强, 等. 玉树地震前的电离层异常现象分析[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]李征航,张小红. 卫星导航定位新技术及高精度数据处理方法[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太阳活动周期太阳风参数及地磁指数的统计分析[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数据的震前电离层异常分析[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); 国家科技基础性工作专项(编号:2015FY310200); 国家重点基础研究发展计划(编号:2013CB733302); 测绘地理信息公益性行业科研专项项目(编号:201512012)

P228.4

A

1008-9268(2016)04-0063-08

邹斌(1992-),男,江西抚州人,硕士生,主要从事空间大地测量研究。

郭金运(1969-),男,山东巨野人,博士、教授、博导,主要从事空间大地测量、海洋大地测量和物理大地测量等研究。

常晓涛(1972-),男,河南人,研究员,主要从事卫星大地测量方面的研究。

联系人： 邹斌 E-mail: binchow77tt@163.com