王 剛 廖洪月 唐好叢

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

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

引言

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

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

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

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

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

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

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

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

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

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

1.2 分析方法

1.2.1 迭代法

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

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

圖2 汶川地震 TEC 活躍性指數(shù)變化圖（a）與 TEC 活躍指數(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 滑動四分位距法

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

四分位距是一種穩(wěn)健統(tǒng)計技術(shù)中用于表示數(shù)據(jù)離散度的一個量，常用來檢查數(shù)據(jù)的異常情況。四分位距法在數(shù)據(jù)分析中具有能較為準(zhǔn)確地獲得預(yù)測參考背景值，并可以比較準(zhǔn)確的計算出背景參考值上下限的特點。所謂的四分位數(shù)就是將數(shù)列分為4部分，一個數(shù)列可以有3個四分位數(shù)，即下四分位數(shù)、中位數(shù)和上四分位數(shù)。以一個16個數(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)計學(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倍，異常檢驗置信度約為95%。

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

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

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

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

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

圖4 北流地震震中周圍格網(wǎng)點TEC時間序列（黑色柱線為發(fā)震時間）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]可知，電離層異常變化易受到如太陽活動變化、地磁活動異常、天氣變化因素等多種因素的影響。因此，研究地震前后電離層異常變化需要考慮時間窗口內(nèi)影響因素的變化情況。在下面的異常分析中將結(jié)合F10.7指數(shù)、Ap指數(shù)、Kp指數(shù)、Dst指數(shù)等主要影響因素進(jìn)行分析。相關(guān)太陽活動和地磁活動數(shù)據(jù)均來源于國家科學(xué)院空間環(huán)境預(yù)報中心（圖5）。

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

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

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

通過剔除對電離層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)較短時間異常（圖6）。因此，可以認(rèn)為該數(shù)據(jù)分析異常不宜作為地震前兆異常。

圖6 10月5日數(shù)據(jù)異常高值時間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日異常高值時間給出10月2日2:00UT—12:00UT間隔兩小時震中周邊區(qū)域TEC空間分布圖?？梢悦黠@看出，TEC異常區(qū)域自2:00UT出現(xiàn)在震中東南方向，隨時間推移逐漸沿22°N由東向西移動，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間隔兩小時震中周邊區(qū)域TEC空間分布圖?？梢悦黠@看出，TEC異常區(qū)域變化同10月2日情況基本一致，除異常區(qū)域自東向西移動外，6:00UT異常區(qū)域出現(xiàn)在震中附近，直至8:00UT異常區(qū)域遠(yuǎn)離震中區(qū)域（圖7b）。

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

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

4 結(jié)語

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

致謝

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