王 剛 廖洪月 唐好叢

1)西安市地震監(jiān)測(cè)中心，西安 710007

2)西安市地震局，西安 710007

引言

中強(qiáng)地震具有較強(qiáng)的破壞性，給人們生命財(cái)產(chǎn)帶來(lái)巨大的損失。然而，地震產(chǎn)生的機(jī)理相當(dāng)復(fù)雜，世界各國(guó)對(duì)地震的預(yù)測(cè)預(yù)報(bào)研究一直處于探索階段。如何提升地震監(jiān)測(cè)預(yù)報(bào)的有效性一直是各國(guó)研究人員面對(duì)的難題，眾多的地震工作者不斷努力的利用各種分析方法在不同領(lǐng)域進(jìn)行探索。

自1964年美國(guó)阿拉斯加大地震時(shí)，Leonard等[1]發(fā)現(xiàn)地震期間電離層有異常擾動(dòng)現(xiàn)象后，各國(guó)學(xué)者便開(kāi)始從事電離層的變化與地震之間關(guān)系的研究工作，并逐漸成為當(dāng)前研究的熱點(diǎn)之一。Antsilevich[2]發(fā)現(xiàn)1966年塔什干地震期間電子含量有明顯增加的現(xiàn)象；Weaver等[3]發(fā)現(xiàn)1969年千島群島地震期間電離層也出現(xiàn)了異常擾動(dòng)。大量研究發(fā)現(xiàn)，大于5級(jí)的地震發(fā)生前幾天到幾個(gè)小時(shí)，震中附近上空電離層都會(huì)出現(xiàn)異?，F(xiàn)象[4-6]。

近年來(lái)，我國(guó)有眾多相關(guān)領(lǐng)域的研究人員在關(guān)于電離層TEC變化與地震的相關(guān)性方面做了大量的研究工作。馬一方等[7]研究發(fā)現(xiàn)蘆山地震前在震中附近上空東向和北向梯度VTEC均出現(xiàn)了負(fù)異常。姚宜斌等[8]利用震中附近GNSS觀測(cè)數(shù)據(jù)，采用滑動(dòng)窗口法發(fā)現(xiàn)2011年3月11日本地震前電離層有擾動(dòng)現(xiàn)象。張小紅等[9]利用IGS提供的震中附近4個(gè)格網(wǎng)點(diǎn)TEC數(shù)據(jù)，采用時(shí)間序列法得到2012年1月10日蘇門答臘島7.2級(jí)地震前13天電離層TEC時(shí)間序列，并分析了時(shí)間序列法、傳統(tǒng)滑動(dòng)時(shí)窗法和四分位法在預(yù)測(cè)電離層TEC參考背景值和精度方面的差異。姚璐等[10]和張明敏等[11]先后研究了玉樹地震和九寨溝地震前后電離層TEC的變化異常。

本文采用迭代法和滑動(dòng)四分位距法對(duì)同一震例進(jìn)行電離層TEC數(shù)據(jù)分析。說(shuō)明不同方法在具體震例中的應(yīng)用及各自的特點(diǎn)，為通過(guò)分析電離層TEC異常預(yù)測(cè)地震發(fā)生積累經(jīng)驗(yàn)。

1 數(shù)據(jù)震例及分析方法

1.1 數(shù)據(jù)及震例

本文采用經(jīng)NASA處理后的IGS全球電離層格網(wǎng)數(shù)據(jù)進(jìn)行電離層異常分析，其空間分辨率為5°（經(jīng)度）×2.5°（緯度），時(shí)間分辨率為 2 小時(shí)。

選取震例為北京時(shí)間2019年10月12日22時(shí)55分發(fā)生在廣西壯族自治區(qū)玉林市北流市（22.18°N，110.51°E）的5.2級(jí)地震（以下簡(jiǎn)稱北流地震），震源深度10 km，圖1為震中及格網(wǎng)點(diǎn)示意圖。

圖1 北流地震震中及格網(wǎng)點(diǎn)示意圖Fig.1 Schematic diagram of epicentral and grid points of Beiliu earthquake

選取該震例主要是為了便于在迭代法研究預(yù)測(cè)發(fā)震區(qū)域的基礎(chǔ)上，采用不同方法對(duì)于其預(yù)測(cè)區(qū)域發(fā)生地震震中周邊電離層TEC數(shù)據(jù)進(jìn)行數(shù)據(jù)分析，在較短時(shí)間序列內(nèi)發(fā)現(xiàn)電離層TEC變化異常，從而可以組成中長(zhǎng)期預(yù)測(cè)指標(biāo)和短臨預(yù)測(cè)指標(biāo)相結(jié)合的前兆異常分析方法。

1.2 分析方法

1.2.1 迭代法

迭代法是一種不斷用變量的舊值遞推新值的過(guò)程。使用迭代法進(jìn)行電離層TEC數(shù)據(jù)分析過(guò)程中，空間分辨率為 5°（經(jīng)度）×2.5°（緯度），時(shí)間分辨率為 2 h，連續(xù)對(duì) 20 年來(lái)全球經(jīng)度（?180°—180°），緯度（?87.5°—87.5°）范圍內(nèi)的 5183（73×71）個(gè)格點(diǎn) TEC 數(shù)據(jù)進(jìn)行系統(tǒng)處理分析。作者研究中為了更好的在較大區(qū)域范圍內(nèi)找到顯著的電離層TEC異常，經(jīng)過(guò)多次嘗試后最終選擇7×7組成的大格網(wǎng)（經(jīng)度跨度35°，緯度跨度17.5°）區(qū)域進(jìn)行分析應(yīng)用。在數(shù)據(jù)分析中，將各格點(diǎn)TEC值數(shù)據(jù)轉(zhuǎn)化成數(shù)據(jù)變化率，在大格網(wǎng)區(qū)域內(nèi)，以本格網(wǎng)為中心，將本格網(wǎng)變化率數(shù)據(jù)與周邊48個(gè)網(wǎng)格單元對(duì)比，若本格網(wǎng)數(shù)據(jù)全部大于或小于周邊網(wǎng)格，則判斷當(dāng)前格網(wǎng)出現(xiàn)了一次異常，然后以365日為滑動(dòng)窗口計(jì)算格網(wǎng)的異常頻次年均線，最后將年均值與歷史最低年均值比較，當(dāng)年均線值比歷史最低年均線值增大1倍時(shí)視為均線異常。

由文獻(xiàn)[12]可知，在TEC數(shù)據(jù)分析中，以年均值線的突變現(xiàn)象作為判斷標(biāo)準(zhǔn)，突升是TEC活躍性增強(qiáng)異常，反之則是TEC活躍性減弱異常，突變指數(shù)越大異常越明顯。經(jīng)對(duì)比研究發(fā)現(xiàn)：日本3?11大地震、汶川地震、玉樹地震、雅安地震、九寨溝地震等地震震前震中區(qū)或周鄰地區(qū)均出現(xiàn)了明顯的TEC活躍性變化特征。以2008年5月12日發(fā)生的汶川8.0地震（31°N，103°E）為例，地震前后震中所在區(qū)域TEC活躍性發(fā)生明顯變化，從2005年開(kāi)始緩慢上升，2007年上升加劇，短時(shí)間內(nèi)就達(dá)到分析方法設(shè)定的異常判斷標(biāo)準(zhǔn)，即增大1倍，2007年底達(dá)到峰值，隨后呈下降趨勢(shì)。之所以出現(xiàn)這種情況，根據(jù)地震發(fā)生的機(jī)理判斷可能是因?yàn)?007年地下板塊運(yùn)動(dòng)達(dá)到頂峰，可運(yùn)動(dòng)空間變小，導(dǎo)致活動(dòng)性減弱所致。TEC活躍性至地震發(fā)生時(shí)仍處于異常判斷標(biāo)準(zhǔn)之上，地震發(fā)生后其活躍性快速下降至異常判斷標(biāo)準(zhǔn)之下（圖2a）。分析以汶川為中心的大區(qū)域TEC活躍性，發(fā)現(xiàn)汶川周圍地區(qū)震前出現(xiàn)高值異常，且這種現(xiàn)象在地震前持續(xù)時(shí)間超過(guò)1年（圖2b）。

圖2 汶川地震 TEC 活躍性指數(shù)變化圖（a）與 TEC 活躍指數(shù)時(shí)段平面圖（2008 年 5 月） （b）Fig.2 Changes of TEC activity index in Wenchuan earthquake （a） and plan diagram of TEC activity index period（May 2008） （b）

1.2.2 滑動(dòng)四分位距法

作者在進(jìn)行迭代法研究電離層變化異常工作的同時(shí)，提出采用滑動(dòng)四分位距法對(duì)預(yù)測(cè)區(qū)域發(fā)生地震前后電離層數(shù)據(jù)進(jìn)行分析的思路，力求尋找電離層TEC變化短臨異常與中強(qiáng)地震的相關(guān)性。

四分位距是一種穩(wěn)健統(tǒng)計(jì)技術(shù)中用于表示數(shù)據(jù)離散度的一個(gè)量，常用來(lái)檢查數(shù)據(jù)的異常情況。四分位距法在數(shù)據(jù)分析中具有能較為準(zhǔn)確地獲得預(yù)測(cè)參考背景值，并可以比較準(zhǔn)確的計(jì)算出背景參考值上下限的特點(diǎn)。所謂的四分位數(shù)就是將數(shù)列分為4部分，一個(gè)數(shù)列可以有3個(gè)四分位數(shù)，即下四分位數(shù)、中位數(shù)和上四分位數(shù)。以一個(gè)16個(gè)數(shù)的數(shù)列為例，將其從小到大排列為x1，x2，……，x16，則

其中，Q1表示在該數(shù)值以下的數(shù)據(jù)占總數(shù)的25%；Q2表示在該數(shù)值以下的數(shù)據(jù)占總數(shù)的50%；Q3表示在該數(shù)值以下的數(shù)據(jù)占總數(shù)的75%。

在統(tǒng)計(jì)學(xué)上IQR=1.34σ，即四分位距的期望值是標(biāo)準(zhǔn)差的1.34倍。在本文數(shù)據(jù)分析中將Q2±1.5IQR作為判定TEC是否為異常值的標(biāo)準(zhǔn)控制值，即上限=Q2+1.5IQR，下限=Q2?1.5IQR。其閾值約為標(biāo)準(zhǔn)差的2倍，異常檢驗(yàn)置信度約為95%。

滑動(dòng)四分位距法可以看做滑動(dòng)窗口法與四分位距法的結(jié)合，就是對(duì)一段時(shí)間序列內(nèi)的數(shù)據(jù)采用四分位距法進(jìn)行數(shù)據(jù)處理，從而得到較為合理的控制限值用來(lái)判讀變化異常。由于電離層變化受季節(jié)影響較大，選擇時(shí)間窗口不宜過(guò)長(zhǎng)。因此，在具體數(shù)據(jù)分析中以15天作為時(shí)間窗口，即分析每天數(shù)據(jù)時(shí)取當(dāng)天及前14天數(shù)據(jù)組成數(shù)列，進(jìn)行數(shù)據(jù)分析，實(shí)際觀測(cè)值超出上下限值即為異常值。

2 TEC 數(shù)據(jù)分析處理

北京時(shí)間2019年10月12日22時(shí)55分在廣西玉林市北流市（22.18°N，110.51°E）發(fā)生 5.2 級(jí)地震，震源深度 10 km。

迭代法分析結(jié)果顯示，自2018年中廣西北部灣出現(xiàn)TEC異常且持續(xù)，直至2019年7月異常達(dá)到近10年來(lái)最大值，隨后異常出現(xiàn)降低，但仍處于分析方法設(shè)定值以上，直至地震發(fā)生，地震發(fā)生后TEC異常呈快速降低形態(tài)，2020年初基本回歸原活動(dòng)范圍（圖3a）。分析以震中區(qū)域?yàn)橹行牡拇髤^(qū)域TEC活躍性，發(fā)現(xiàn)震中周圍地區(qū)震前出現(xiàn)高值異常，且這種現(xiàn)象在地震前持續(xù)時(shí)間超過(guò)1年（圖3b）。

圖3 北流地震 TEC 活躍性指數(shù)變化圖（a）與 TEC 活躍指數(shù)時(shí)段平面圖（2019 年 7 月） （b）Fig.3 Variation chart of TEC activity index of Beiliu earthquake （a） and time interval plan of TEC activity index（July 2019） （b）

選取2019年9月1日—10月16日期間，地震發(fā)生時(shí)間前15天，發(fā)生后3天組成的時(shí)間序列，以15日為時(shí)間窗口應(yīng)用滑動(dòng)四分位法進(jìn)行數(shù)據(jù)分析處理。下面圖中依次顯示了震中周圍各格網(wǎng)點(diǎn)的TEC異常變化情況。通過(guò)數(shù)據(jù)分析可以得出，各格網(wǎng)點(diǎn)TEC異常變化具有較好的一致性，震前10—14天、7天、3天、2天、當(dāng)日均出現(xiàn)正值異常，震前1天出現(xiàn)負(fù)值異常（圖4）。

圖4 北流地震震中周圍格網(wǎng)點(diǎn)TEC時(shí)間序列（黑色柱線為發(fā)震時(shí)間）Fig.4 TEC time series of grid points around the epicenter of the Beiliu earthquake （the black column is the time of seismogenesis）

根據(jù)相關(guān)研究文獻(xiàn)[13]可知，電離層異常變化易受到如太陽(yáng)活動(dòng)變化、地磁活動(dòng)異常、天氣變化因素等多種因素的影響。因此，研究地震前后電離層異常變化需要考慮時(shí)間窗口內(nèi)影響因素的變化情況。在下面的異常分析中將結(jié)合F10.7指數(shù)、Ap指數(shù)、Kp指數(shù)、Dst指數(shù)等主要影響因素進(jìn)行分析。相關(guān)太陽(yáng)活動(dòng)和地磁活動(dòng)數(shù)據(jù)均來(lái)源于國(guó)家科學(xué)院空間環(huán)境預(yù)報(bào)中心（圖5）。

圖5 9月27日—10月15日北流地震前后太陽(yáng)活動(dòng)和地磁變化情況Fig.5 Changes of solar activity and geomagnetism before and after the Beiliu earthquake from September 27 to October 15

根據(jù)國(guó)家科學(xué)院空間環(huán)境預(yù)報(bào)中心對(duì)2019年9—10月份的太陽(yáng)活動(dòng)和地磁觀測(cè)數(shù)據(jù)報(bào)告可知，9月份太陽(yáng)活動(dòng)水平較低，無(wú)C級(jí)以上級(jí)別耀斑產(chǎn)生；地磁活動(dòng)9月27—30日出現(xiàn)小磁暴現(xiàn)象，對(duì)TEC變化有所影響。10月份太陽(yáng)活動(dòng)水平依然較低，無(wú)C級(jí)以上級(jí)別耀斑產(chǎn)生；地磁活動(dòng)10月9日Kp指數(shù)、Ap指數(shù)處于高值，10月10—11日，太陽(yáng)風(fēng)速上升至510 m/s左右，地磁活動(dòng)水平較高，會(huì)對(duì)TEC變化有所影響。

因此，可以判斷震前11—14天，震前1—3天TEC變化異常受到地磁影響因素較大，不宜作為地震前兆異常。

通過(guò)剔除對(duì)電離層TEC異常影響因素后，獲得電離層TEC異常后還有震前10天、7天及地震當(dāng)日3組異常，即10月2日、10月5日、10月12日。

3 TEC 空間分布異常分析

為了進(jìn)一步分析震前3組TEC異常是否與此次地震有關(guān)，本文分別根據(jù)數(shù)據(jù)分析結(jié)果給出了3組異常高值前后（5°—40°N，75°—125°E）范圍內(nèi)電離層TEC異常分布圖。

根據(jù)數(shù)據(jù)分析結(jié)果給出的10月5日異常高值出現(xiàn)在14:00UT，據(jù)此給出10月5日10:00UT—20:00UT之間的TEC空間分布圖，發(fā)現(xiàn)震中附近區(qū)域僅在18:00UT出現(xiàn)較短時(shí)間異常（圖6）。因此，可以認(rèn)為該數(shù)據(jù)分析異常不宜作為地震前兆異常。

圖6 10月5日數(shù)據(jù)異常高值時(shí)間TEC空間分布圖（黃色標(biāo)記為震中位置）Fig.6 Spatial distribution of TEC with abnormally high data on October 5 （epicenter location is marked in yellow）

根據(jù)數(shù)據(jù)分析10月2日異常高值時(shí)間給出10月2日2:00UT—12:00UT間隔兩小時(shí)震中周邊區(qū)域TEC空間分布圖。可以明顯看出，TEC異常區(qū)域自2:00UT出現(xiàn)在震中東南方向，隨時(shí)間推移逐漸沿22°N由東向西移動(dòng)，8:00UT異常區(qū)域出現(xiàn)在震中附近，直至12:00UT遠(yuǎn)離震中區(qū)域（圖7a）。

圖7 10月2日、12日TEC空間分布圖（黃色標(biāo)記為震中位置）Fig.7 The spatial distribution of TEC on October 2 and 12 （the epicenter is marked in yellow）

根據(jù)數(shù)據(jù)分析10月12日異常高值給出10月12日2:00UT—12:00UT間隔兩小時(shí)震中周邊區(qū)域TEC空間分布圖。可以明顯看出，TEC異常區(qū)域變化同10月2日情況基本一致，除異常區(qū)域自東向西移動(dòng)外，6:00UT異常區(qū)域出現(xiàn)在震中附近，直至8:00UT異常區(qū)域遠(yuǎn)離震中區(qū)域（圖7b）。

綜合考慮可以認(rèn)為，10月2日與10月12日的電離層TEC異常與廣西北流地震有關(guān)。

因此，對(duì)中強(qiáng)地震發(fā)生前后數(shù)天內(nèi)電離層TEC數(shù)據(jù)采用滑動(dòng)四分位距法進(jìn)行分析，盡管易受到太陽(yáng)活動(dòng)和地磁變化的影響，但將去除干擾因素后的異常作為地震前兆短臨判斷指標(biāo)，仍具有一定的有效性。這一結(jié)論與其他研究人員分析地震前電離層TEC變化異常的研究結(jié)果基本一致[14-17]。

4 結(jié)語(yǔ)

本文采用迭代法和滑動(dòng)四分位距法對(duì)同一震例中電離層TEC異常進(jìn)行了分析。迭代法具有分析較長(zhǎng)時(shí)間內(nèi)大區(qū)域電離層TEC異常的特點(diǎn)，滑動(dòng)四分位距法可以對(duì)單個(gè)格網(wǎng)較短時(shí)間序列內(nèi)電離層TEC異常進(jìn)行分析。在地震監(jiān)測(cè)預(yù)報(bào)日常工作中可以嘗試將迭代法和滑動(dòng)四分位距法分析電離層TEC異常的方法綜合使用來(lái)進(jìn)行地震前兆判斷研究，以達(dá)到中長(zhǎng)期預(yù)報(bào)和短臨預(yù)報(bào)相結(jié)合的目的。但由于地震發(fā)生的機(jī)理復(fù)雜多變，利用電離層擾動(dòng)對(duì)地震進(jìn)行預(yù)測(cè)存在明顯的局限性，還有待于參考多種不同分析方法進(jìn)行研究和震例印證，來(lái)提高地震預(yù)測(cè)的有效性。

致謝

感謝陜西省地震局方偉高工給予的指導(dǎo)和支持，以及西安市地震監(jiān)測(cè)中心同事在研究過(guò)程中提供的幫助。