朱軍桃,趙苗興,龔朝飛,王 雷

(桂林理工大學(xué) 測繪地理信息學(xué)院,廣西 桂林 541006)

地震具有突發(fā)性強、破壞性大、防御難度大且社會影響深遠(yuǎn)等特點,嚴(yán)重危害人類的生命財產(chǎn)安全。據(jù)中國地震臺網(wǎng)(http://www.ceic.ac.cn)統(tǒng)計,2016年9月至今,在世界范圍內(nèi)超過6.0級的地震有198起,造成的人員傷亡和經(jīng)濟(jì)損失讓人們難以承受。由于地震產(chǎn)生原因復(fù)雜,影響因素眾多,地震預(yù)報特別是短臨預(yù)報目前仍然處于探索階段[1]。1964年美國阿拉斯加大地震時,Leonard等[2]首次發(fā)現(xiàn)電離層擾動與地震之間存在某些聯(lián)系,從此拉開了震前電離層異常擾動研究的序幕;Pulinets等[3-4]提出了利用衛(wèi)星數(shù)據(jù)進(jìn)行地震短臨預(yù)報,統(tǒng)計分析了10年的地震電離層數(shù)據(jù),確認(rèn)電離層在地震前后有顯著的異常現(xiàn)象;林劍等[5]利用GPS數(shù)據(jù)研究了汶川地震前后電離層的異常擾動,發(fā)現(xiàn)地震前后一周,孕震區(qū)上空連續(xù)出現(xiàn)電離層異常擾動,異常形態(tài)具有共軛結(jié)構(gòu),且呈現(xiàn)向赤道漂移趨勢;童楊津等[6]采用滑動四分位法發(fā)現(xiàn)墨西哥M7.2地震前電離層TEC出現(xiàn)正異?，F(xiàn)象;楊力等[7]采用滑動時窗法和相關(guān)性分析法研究了日本九州島地震震前電離層異常變化。隨著電離層探測技術(shù)的不斷發(fā)展與完善,孕震區(qū)電離層異常擾動成了當(dāng)前地震預(yù)報的熱點之一。除了傳統(tǒng)的電離層異常分析方法外,姚宜斌等[8]提出了一種基于奇異譜分析的電離層異常探測方法,結(jié)合三維層析發(fā)現(xiàn)2015年尼泊爾地震前出現(xiàn)大范圍電離層正異常,異常隨高度變化規(guī)律與電離層電子密度垂直分布規(guī)律相一致;張小紅等[9]提出一種利用時間序列法(ARIMA模型)進(jìn)行震前電離層異常探測的新方法,發(fā)現(xiàn)時間序列法預(yù)測背景值的精度要明顯高于傳統(tǒng)方法,且預(yù)報背景值的平均偏差要比傳統(tǒng)方法小2倍左右;黃海莎等[10]介紹了震前電離層TEC異常探測原理的研究進(jìn)展和主要的異常探測方法,總結(jié)了時間序列法、Kalman濾波和小波變換在地震電離層異常分析中的應(yīng)用, 為研究震前電離層異常擾動提供了新的思路; 楊可可等[11]利用IGS數(shù)據(jù)分析中心提供的全球電離層格網(wǎng)數(shù)據(jù),用雙線性插值法和滑動四分位法對日本九州島地震上空電離層變化特性進(jìn)行分析,探討地震對電離層變化情況的影響。

2017-08-08T21:19:46(年積日220 d, 協(xié)調(diào)世界時UT 13:19:46)九寨溝(33.20°N, 103.82°E)發(fā)生MS7.0級地震,震源深度20 km,造成重大生命財產(chǎn)損失。本文以此次地震為背景,利用GIM數(shù)據(jù)分析地震前后電離層異常擾動,充分考慮太陽活動、地磁變化等因素，分析電離層異常與地震的關(guān)聯(lián)性。

1 數(shù)據(jù)和方法

本文采用CODE(The Center for Orbit Determination in Europe)提供的GIM(global ionosphere maps)數(shù)據(jù)(ftp://cddis.gsfc.nasa.gov/)進(jìn)行電離層異常研究, 該數(shù)據(jù)時間分辨率為1 h, 空間分辨率為2.5°×5°, 單位為TECU(1 TECU=1016個電子/m2)。CODE位于瑞士伯尼爾大學(xué), 是IGS(International GNSS Service)的7個分析中心之一, 其發(fā)布的GIM數(shù)據(jù)經(jīng)過一定的插值及平滑處理,過濾掉了小范圍的電離層擾動,非常適合研究全球大尺度的震前電離層異常擾動。本次九寨溝地震震級高、強度大,根據(jù)以往的研究成果分析，其很可能引起了大尺度的電離層擾動。在CODE制作GIM的過程中,使用了武漢、拉薩、烏魯木齊、昆明、上海等GPS觀測站的數(shù)據(jù),因此可以采用CODE的GIM數(shù)據(jù)分析九寨溝地震期間的電離層異常擾動現(xiàn)象。在日地環(huán)境方面,本文采用中科院空間環(huán)境研究預(yù)報中心(space environment prediction center, SEPC)發(fā)布的赤道地區(qū)地磁活動指數(shù)(Dst)、 全球地磁活動指數(shù)(Kp)和10.7 cm射電流量(F10.7)( http://www.sepc.ac.cn/)反映太陽活動和地磁變化情況。

在異常檢測之前,采用“db4”正交小波和軟閾值函數(shù)對GIM數(shù)據(jù)進(jìn)行小波閾值去噪處理,去噪的原理及具體過程見文獻(xiàn)[12-13]。小波閾值去噪可以很好地保護(hù)有用的信號尖峰和突變信號,并且有效去除暫態(tài)信號和瞬態(tài)信號引起的噪聲。經(jīng)小波閾值去噪后,有效地去除了地震異常檢測中的噪聲干擾。

本文采用滑動四分位法[14-15]進(jìn)行震前電離層的異常檢測。在時窗長度的選取方面,考慮到TEC變化具有季節(jié)效應(yīng),時窗長度不宜過長,且本文主要研究震前的電離層TEC異?，F(xiàn)象,因此時窗長度選擇15天。選取年積日第205～223天共19 d的TEC數(shù)據(jù)作為待檢測數(shù)據(jù),數(shù)據(jù)包含震前15天和震后3天,滑動四分位法檢測TEC異常的具體過程如下:

將某一時刻待測數(shù)據(jù)及其前15天同一時刻數(shù)據(jù)按升序排列,得到I1、I2、 …、I16, 計算其四分位數(shù)Q1、Q2、Q3[11]

(1)

四分位(IQR)為

IQR=Q3-Q1,

(2)

該時刻TEC上限UP和下限LOW可表示為

(3)

其中k為倍常數(shù), 決定篩選異常的閾值范圍。本文k值取1.5,即取1.5倍的IQR作為異常篩選的上下限。 假設(shè)待檢測數(shù)據(jù)服從正態(tài)分布, UP和LOW用均值(μ)和標(biāo)準(zhǔn)差(δ)可表示為μ±2δ,即與平均值的偏差超過2倍標(biāo)準(zhǔn)差的值被認(rèn)定為異常值,置信度為95%[16]。若某一時刻待檢測數(shù)據(jù)高于上限或低于下限,則表明該時刻TEC出現(xiàn)正異?；蜇?fù)異常。

2 電離層異常擾動分析

2.1 日地環(huán)境分析

太陽活動、地磁變化等諸多因素都會引起電離層的異常擾動,因此,在研究地震與電離層異常的關(guān)系時首先要分析孕震區(qū)地震前后日地環(huán)境變化。本文采用滑動四分位法對年積日第205～223天的Dst、Kp和F10.7進(jìn)行了異常檢測(倍常數(shù)k=1.5),Dst、Kp、F10.7變化情況及異常檢測結(jié)果如圖1所示。

圖1 Dst、Kp、F10.7變化情況(a)及異常探測結(jié)果(b)Fig.1 Variation(a) and abnormal detection results(b) of Dst,Kp,F10.7 index

結(jié)果顯示,年積日第216天Dst和Kp均出現(xiàn)正異?，F(xiàn)象,說明該日地磁活動強烈且存在磁暴現(xiàn)象,若該日出現(xiàn)電離層TEC異常擾動現(xiàn)象需要特別關(guān)注。除年積日第216天以外,此次地震前后太陽及地磁活動平靜,日地環(huán)境平穩(wěn)。

2.2 TEC時間序列異常分析

本文選取震中(33.20°N, 10.3.82°E)作為檢測點, 震中TEC采用雙線性插值法內(nèi)插獲得。 同時, 選取同一緯度的點A(32.5°N, 100°E)和點B(32.5°N, 110°E)、 同一經(jīng)度的點C(35°N, 105°E)和點D(30°N, 105°E)作為參考點, 如圖2所示,

圖2 地震信息分布Fig.2 Seismic information distribution

孕震區(qū)半徑由公式R=100.43M計算得出(R為半徑,M為震級)。利用滑動四分位法對震中及各參考點TEC時間序列進(jìn)行檢測, 并繪制異常分布圖, 結(jié)果如圖3、 圖4所示, 其中黑色縱線為震發(fā)時刻。

圖3 震中TEC時間序列及異常分布異常Fig.3 Epicenter TEC time series and anomaly distribution

震中在年積日第216和218天檢測到TEC正異?，F(xiàn)象, 異常持續(xù)時間分別為7和5 h, 異常時刻集中在(UT)00:00—08:00, 即當(dāng)?shù)貢r間(LT)08:00—16:00, 其中, 第216天異常峰值為2 TECU, 第218天異常峰值較大, 達(dá)6 TECU。 在年積日第208天, 震中檢測到TEC負(fù)異?，F(xiàn)象, 異常值為-2 TECU, 異常時刻為(UT)23:00, 即當(dāng)?shù)貢r間(LT)次日的07:00。 整體來看, 震中在震前有2天出現(xiàn)TEC正異常, 有1天出現(xiàn)TEC負(fù)異常, 正異?？倳r長為12 h, 負(fù)異常總時長為1 h; 震后未發(fā)現(xiàn)TEC異?，F(xiàn)象。

分析圖4,參考點的檢測結(jié)果與震中存在高度的一致性。年積日第218天,參考點A和參考點D發(fā)現(xiàn)的TEC正異常規(guī)模大于其他檢測點; 年積日第216天, 震中和參考點C發(fā)現(xiàn)的TEC正異常峰值較大; 年積日第208天, 參考點C發(fā)現(xiàn)的TEC負(fù)異常達(dá)-4 TECU,參考點A和參考點D負(fù)異常只有-0.1 TECU左右。因此,震前出現(xiàn)的TEC異常區(qū)域不垂直對應(yīng)震中地區(qū)。

圖4 參考點TEC時間序列及異常分布Fig.4 TEC time series and anomaly distribution of reference point

綜上所述, 孕震區(qū)在當(dāng)?shù)貢r間2017年7月28日(年積日第208天)、 8月4日(年積日第216天)和8月6日(年積日第218天)檢測到TEC異?，F(xiàn)象,這些異?，F(xiàn)象是否由孕震引起還需結(jié)合日地環(huán)境和異?？臻g分布進(jìn)一步分析。

2.3 全球TEC異常分析

利用滑動四分位法只能進(jìn)行單點的TEC異常檢測,不能確定異常的形態(tài)與分布。為了進(jìn)一步判斷電離層異常與九寨溝地震的相關(guān)性,本文對GIM數(shù)據(jù)的每個格網(wǎng)點進(jìn)行異常檢測,并制得當(dāng)?shù)貢r間2017年7月28日、8月4日和8月6日的全球TEC異常分布圖,如圖5—7所示,其中“★”表示震中,圓圈代表孕震區(qū),橫線代表赤道。

圖5給出了2017年7月28日間隔2 h的全球TEC異常分布,可以看出,該日全球TEC異常主要為負(fù)異常且集中在赤道附近(15°S—15°N),并沿赤道自東向西移動,該現(xiàn)象屬于赤道異常,與孕震關(guān)系不大。孕震區(qū)周圍在LT06:00開始出現(xiàn)TEC負(fù)異常, 異常位于孕震區(qū)西南方向, 隨后異常逐漸擴(kuò)大、 向西移動并靠近赤道, 至LT10:00異常消失, 負(fù)異常最大約-2 TECU。 結(jié)合該日日地環(huán)境平穩(wěn),且其他區(qū)域未出現(xiàn)類似現(xiàn)象,因此,該日孕震區(qū)周圍出現(xiàn)的電離層TEC負(fù)異?，F(xiàn)象可能是孕育地震引起,異常區(qū)域沒有垂直對應(yīng)于震中,而是位于震中的西南方向并逐漸向赤道漂移。

圖5 2017年7月28日全球TEC異常分布圖Fig.5 Global TEC anomaly maps on July 28,2017

圖6給出了2017年8月4日間隔2 h的全球TEC異常分布圖, 該日全球出現(xiàn)大范圍的TEC正異?，F(xiàn)象,異常在全球范圍內(nèi)呈帶狀分布,南半球異常的規(guī)模、幅度及持續(xù)時間均大于北半球。結(jié)合前文對日地環(huán)境的分析,該日地磁活動強烈且存在磁暴現(xiàn)象,因此,該日電離層TEC正異?，F(xiàn)象可能是地磁擾動引起。雖然孕震區(qū)域也出現(xiàn)了長達(dá)8 h的TEC正異常,但是該異常現(xiàn)象是全球性的,在其他區(qū)域也有出現(xiàn),并不能作為地震前兆。

圖6 2017年8月4日全球TEC異常分布圖Fig.6 Global TEC anomaly maps on August 4,2017

圖7給出了2017年8月6日間隔2 h的全球TEC異常分布圖。除了赤道異常以外,北半球只有孕震區(qū)周圍在該日出現(xiàn)了長時間的TEC正異常現(xiàn)象,異常從LT04:00持續(xù)到LT14:00,最大異常值達(dá)10 TECU,最大異常出現(xiàn)在LT14:00,位于震中西南方,異常區(qū)域由震中的東北方逐漸向西南方移動,表現(xiàn)出向赤道漂移的趨勢。同時,孕震區(qū)對應(yīng)的赤道共軛區(qū)內(nèi)也發(fā)現(xiàn)明顯的TEC正異?，F(xiàn)象,赤道共軛區(qū)的異常范圍大于孕震區(qū)的異常范圍,異常區(qū)域也表現(xiàn)出向赤道漂移的趨勢?？紤]到該日日地環(huán)境平穩(wěn),所以孕震區(qū)周圍出現(xiàn)的電離層TEC正異?，F(xiàn)象可能是孕育地震引起,異常區(qū)域及異常峰值沒有垂直對應(yīng)于震中,異常區(qū)域自東北方向穿過震中逐漸向赤道漂移且具有共軛結(jié)構(gòu)。

圖7 2017年8月6日全球TEC異常分布圖Fig.7 Global TEC anomaly maps on August 6,2017

3 結(jié) 論

本文基于CODE提供的GIM數(shù)據(jù), 采用滑動四分位法作為異常檢測方法, 分析了2017年8月8日九寨溝MS7.0級地震前后電離層的異常情況。 通過對震中及周圍4個參考點TEC時間序列的異常檢測, 發(fā)現(xiàn)地震孕育區(qū)內(nèi), 2017年7月28日出現(xiàn)電離層TEC負(fù)異?，F(xiàn)象, 2017年8月4日和8月6日出現(xiàn)電離層TEC正異?，F(xiàn)象。 結(jié)合日地環(huán)境和全球TEC異常分布進(jìn)一步分析后得出: 2017年7月28日出現(xiàn)的電離層TEC負(fù)異?，F(xiàn)象和8月6日出現(xiàn)的電離層TEC正異常現(xiàn)象可能是本次地震的前兆,其中7月28日異常出現(xiàn)在LT06:00—10:00, 異常區(qū)域沒有垂直對應(yīng)于震中, 而是位于震中的西南方向并逐漸向赤道漂移, 8月6日異常出現(xiàn)在LT04:00—14:00,異常區(qū)域沒有垂直對應(yīng)于震中,異常區(qū)域自東北方向穿過震中逐漸向赤道漂移且具有共軛結(jié)構(gòu);而2017年8月4日出現(xiàn)的電離層TEC正異?，F(xiàn)象可能是由地磁擾動引起。

地震電離層異常的耦合機(jī)制是個非常復(fù)雜的物理和化學(xué)過程,目前,還沒有很好的物理模型對其進(jìn)行解釋?，F(xiàn)階段地震電離層異常擾動的研究仍處于探索階段,要將震前電離層異常擾動做為地震預(yù)報的一種手段還有許多問題亟待解決。本文的研究成果再次證明了強震前確實存在電離層異常擾動,為研究地震電離層異常的耦合機(jī)制提供了參考和依據(jù)。