結合PCA和滑動時窗的電離層異常探測方法研究

謝 婷 陳必焰 匡翠林

1 中南大學地球科學與信息物理學院, 長沙市麓山南路932號, 410083

多年的觀測實踐與大量的震例研究[1-3]表明，地震發生前電離層會出現異常變化。目前，傳統電離層異常探測的方法主要包括中位數法、平均值法、滑動四分位距法和滑動時窗法等，而傳統電離層異常探測方法是基于參考背景值與標準差來進行探測的，這就決定了高精度的參考背景值和合理的異常判別限差是獲取精準異常探測結果的前提。基于此，研究者們陸續提出了許多異常探測新方法[4-5]，能更有效地提取信息、細化信號，并提高結果的可靠性，對地震異常探測具有積極的影響。隨著研究的深入，有學者也嘗試了一些無需構造背景場的異常探測方法，如主成分分析(PCA)方法[6-8]。

但是，利用PCA方法探測異常也存在一些不足之處：1)PCA方法無法提取出異常的空間特征，不利于震前異常研究；2)因為主特征值容易出現極端高值(如磁暴)，在進行長時間序列異常探測時，直接利用主特征值進行異常分析會削弱PCA方法的異常探測能力；3)目前的PCA方法直接采用經驗值(0.5)為閾值進行異常判斷，容易出現疏漏和失誤，且檢測得到的結果無法區分正負異常。針對這些問題，本文對PCA方法作出了改進，首先利用PCA方法初步過濾空間環境(地磁活動和太陽活動)等非震因素對電離層TEC擾動的影響，再采用滑動時窗方法進行電離層異常探測，同時與傳統探測方法(滑動四分位距法和滑動時窗法)進行對比分析，驗證結合PCA與滑動時窗方法的有效性，并利用本文方法初探了2020年瓦哈卡地震震前電離層時空變化特性，為后續地震研究工作提供可靠信息。

1 結合PCA與滑動時窗的探測方法

PCA是由Pearson[9]提出的，在數據挖掘和數據分析等領域得到了廣泛應用。PCA方法主要是通過特征分解得到數據的主成分(特征向量)與特征值，進而獲得數據的主要信息。

本文基于全球電離層網格(GIM)數據，研究PCA提取的TEC主特征值所代表的物理含義，以2014-01-12～02-12主特征值時空圖為例進行探討。從圖1可以看出，在緯度±15°處，主特征值與地磁線走勢吻合，且主特征值會隨時間和位置發生變化，其中02-02于田地震震中上方區域觀測到明顯的高值異常，該異常關于磁赤道線對稱，這與研究發現的震前電離層異常特征[10]是一致的。基于這些特征推測，PCA分離出的主特征值可能是基礎背景場與異常的疊加。

圖中，五角星為2014-02-12于田地震震中位置，圓圈為地震影響范圍，紅色曲線為緯度±15°處的地磁線

滑動時窗法是一種在給定的特定窗口大小的數組或字符串上執行要求操作的算法，其操作簡單、運行快速，是一種常用的電離層異常探測方法。但該方法探測結果的準確度取決于背景值的精度，易受太陽活動等非震因素的影響，因此其探測精度有限。而PCA提取的主信號能初步過濾掉太陽活動等非震因素對電離層TEC的影響，在PCA方法的基礎上采用滑動時窗法進行異常探測能更有效地檢測出地震引起的電離層異常，且滑動時窗主要是基于滑動窗口來進行異常探測的，可以很好地解決PCA方法時序探測有限的問題。

結合PCA和滑動時窗(PCA-滑動時窗)方法的具體數據處理流程為：1)基于單點TEC時間序列，利用PCA求解出特征值；2)對特征值進行排序，求得最大特征值(單天1個主特征值)，多天的數據處理可以得到TEC的主特征值序列；3)采用滑動時窗法對主特征值序列進行檢測，滑動窗口為16 d，以1.5倍均方差為探測依據，得到異常序列；4)由圖1可知，主特征值之間相差較大會導致異常大小無法被很好地對比分析，因此對異常數據進行標準化處理(該處理不影響異常探測結果)。其中，電離層異常的判定標準為：異常出現在震中位置附近獨立存在的區域，在位置或數值上與赤道電離層異常區[11](磁緯0°～±15°)有明顯分離[2]。

2 分析與討論

本文分別以2008年汶川地震、2010年玉樹地震及2014年于田地震為例，研究PCA-滑動時窗方法的異常探測能力，其中汶川地震和玉樹地震震前地磁活動頻繁，且玉樹地震震前地磁活動的劇烈程度較高，可用來探究PCA方法的抗地磁活動干擾能力；而于田地震震前的太陽活動可用來研究太陽活動對PCA方法探測異常的影響。最后，利用PCA-滑動時窗方法研究2020年瓦哈卡地震的震前異常現象，為后續地震研究工作提供可靠信息。4次地震的基本信息見表1。

表1 4次地震的基本信息

2.1 汶川、玉樹及于田地震分析

本文選取了2008-01-01～05-20、2010-01-01～05-20及2014-01-01～05-20的GIM數據，采用PCA-滑動時窗方法進行電離層時空異常探測。由于不同震級對應孕震區的范圍不同，因此采用經驗公式R=100.43M(R為孕震區半徑，M為震級)計算地震的影響范圍[12]。此外，本文還選取了太陽活動和地磁活動指數數據，以表征地磁及太陽的活動情況，其中Dst為中低緯地區地磁活動指數，Dst >-30 nT表明地磁處于平靜狀態，-50 nT < Dst <-30 nT表明可能發生弱磁暴，Dst <-100 nT表明可能發生大磁暴；KP為一種判斷地磁擾動總強度的指標，可將KP= 0～2、3、4、5、6、7～9的地磁活動分為平靜、不平靜、活動、小、大、強烈等不同等級；F10.7為表示太陽活動強弱的常用指標，其值通常在60～300之間，F10.7平穩微弱浮動表示太陽活動平靜。

2.1.1 3個地震的時間序列異常分析

滑動四分位距法[13]是目前被認可并廣泛使用的異常探測方法，同時滑動時窗法也是傳統方法中相對合理的一種方法，因此本文從已有研究中選取這2種方法的震前異常探測結果，并結合PCA方法的探測結果與PCA-滑動時窗方法進行驗證對比。由于目前PCA方法僅用于單點探測，因此本文基于GIM數據內插出震中TEC值的時間序列，再采用PCA方法進行異常探測，其中表2為4種方法探測得到的前3個地震的異常結果，滑動四分位距法和滑動時窗法的探測結果在不同研究中有所不同，但總體上是一致的。

表2 基于PCA-滑動時窗方法、PCA和前人研究資料匯總得到的3個地震的探測結果

表2中PCA-滑動時窗方法的探測結果沒有排除空間環境(太陽活動和地磁活動)的影響，但得到的異常值與前人研究中排除空間環境影響的結果基本一致，因此可初步推測PCA-滑動時窗方法探測異常不易受地磁活動和太陽活動的影響。分析探測結果可發現，PCA方法在汶川地震震前21 d和19 d探測到了異常，但經過對孕震區內多點的異常進行檢測發現該異常不具備區域性，即該異常只存在于相鄰兩點之間，這也表明僅利用震中上空單點進行異常探測存在局限性，需進一步開展時空異常探測。本文結果有別于湯俊等[6]利用PCA方法探測得到的結果，因為本文選取了2008-01-01～05-20共141 d的數據來進行異常探測。在進行長時間序列異常探測時，主特征值容易出現異常高值，直接利用其進行異常分析會削弱PCA方法的異常探測能力，這也是該方法的缺陷之一。從表2還可以看出，PCA-滑動時窗方法對這一缺陷有很好的改進，其結果與前人研究結果有很好的一致性。由于前人研究地震前兆一般只研究震前16 d左右異常，因此利用滑動四分位距法和滑動時窗法進行震前異常分析的資料匯總時間均不超過16 d。并且通常來說，如果某天有地磁或太陽活動，那么當天的電離層會出現明顯異于平靜期的擾動，且該擾動一般是全球尺度的，其中赤道區域電離層最易受太陽和地磁影響[17]，其擾動程度會更大。

2.1.2 汶川地震空間異常分析

圖2(a)～2(c)給出了2008-04-21～05-14期間各時段的KP指數、Dst指數及F10.7指數，由圖可知，汶川地震前后除了04-23、04-24、04-26、04-28、04-30及05-01和05-02時段的KP指數大于3，且04-23、04-24、04-25和04-26時段的Dst指數低于-30 nT外，其余時段的KP指數和Dst指數基本處于較低水平，而F10.7指數在整個時段內都非常穩定。

圖2(d)為PCA-滑動時窗方法探測得到的二維異常分布，可以看出，該方法探測到的異常時段為04-29(震前13 d)、05-02(震前10 d)、05-05(震前7 d)、05-06(震前6 d)、05-07(震前5 d)和05-09(震前3 d)，這與前人研究結果具有極高的一致性，其中與滑動四分位距法和滑動時窗法相比，PCA-滑動時窗方法檢測到的異常更全面。從圖2(e)可以看出，PCA方法檢測出的異常時段為04-21和04-23(震前21 d和19 d)。由圖2(a)～2(c)可知，04-23～26的Dst指數低于-30 nT，查詢空間環境預報中心的報告[18]發現，地磁活動在04-23處于中等擾動水平，對比圖2(d)的異常空間分布可知，04-23異常主要出現在磁赤道線附近。從異常情況來看，04-29和05-05～07異常減少，05-02南部異常減弱、東部異常增大，05-09異常增加，其中電離層異常出現在震中偏東南方向，呈共軛結構，這與前人采用其他方法探測得到的結果是一致的[1,15]，證實PCA-滑動時窗方法可準確區分正負異常，彌補了PCA方法的缺陷。

(a)Dst指數；(b)KP指數；(c)F10.7指數(虛線是各指數設定的閾值，Dst為-30 nT，KP為3)；(d)基于PCA-滑動時窗法探測得到的二維異常分布(五角星為2008-05-12汶川地震震中位置，紅色曲線為磁赤道線(磁緯0°)和磁緯±15°)；(e)利用PCA分析電離層TEC所得震中位置的主特征值序列(水平紅線為0.5閾值線，垂線紅線為地震發生日)

2.1.3 玉樹地震空間異常分析

通過分析圖3的太陽活動和地磁活動數據可以發現，04-05～08和04-11～12全球發生了中等地磁暴，但在圖3(d)對應時間段的異常分布中并沒有出現全球尺度或大范圍的異常，驗證了PCA方法探測異常不易受地磁活動干擾這一推測。根據空間環境預報中心發布的報告[18]可知，04-05全球發生了一次強烈的空間天氣事件，并從04-06開始，地球同步軌道連續6 d發生高能電子通量增強事件，命名為“清明節事件”。因此，該時段內的TEC異常很可能是地磁活動引起的，本文方法在04-05～10檢測到的異常位置主要位于赤道異常區和美洲扇區，與磁暴引起的電離層特征[17]吻合。相比于利用傳統方法探測到的全球尺度大范圍異常，PCA-滑動時窗方法會削弱地磁活動對電離層的影響。另外，整個時段內太陽活動水平較低，F10.7指數變化幅度小，整體維持在75～80 sfu附近，因此PCA-滑動時窗方法探測到的03-25(震前19 d)異常為正異常，03-30(震前14 d)異常為負異常，04-01(震前12 d)異常為正異常，這3次電離層異常都出現在震中偏東南方向，呈共軛結構。04-13地震當天震中以南部分區域出現TEC異常增加的現象，與前人研究獲取的地震電離層擾動特征非常相似[1-2]，很可能是地震引起的電離層異常。

(a)Dst指數；(b)KP指數；(c)F10.7指數(虛線是各指數設定的閾值，Dst為-30 nT，KP為3)；(d)基于PCA-滑動時窗法探測得到的二維異常分布(五角星為2010-04-13玉樹地震震中位置，紅色曲線為磁赤道線(磁緯0°)和磁緯±15°)；(e)利用PCA分析電離層TEC所得震中位置的主特征值序列(水平紅線為0.5閾值線，垂線紅線為地震發生日)

2.1.4 于田地震空間異常分析

于田地震前后太陽活動和地磁活動情況如圖4(a)～4(c)所示，從圖4(c)中可以看出，震前太陽活動并不穩定，在01-22、01-28、02-05和02-06出現了異常上升或下降，太陽活動頻繁，其中F10.7為200 sfu的峰值出現在震前6 d、7 d和8 d。

圖4(a)和4(b)分別顯示了Dst指數和KP指數的時間序列，綜合分析2個參數發現，在02-07之前，Dst指數大于-30 nT，KP指數小于3，但在02-07～10期間，Dst指數達到-40 nT，KP指數超過5，查詢空間環境預報中心[18]資料發現，該時段內地磁有3 h達到小磁暴水平，21 h達到活躍水平。

綜上所述，太陽活動頻繁期為01-22(震前21 d)、01-28(震前15 d)、02-04(震前8 d)、02-05(震前7 d)和02-06(震前6 d)，地磁活動頻繁期為02-07～10(震前2～5 d)。從圖4(d)的異常時空分布中可發現，在太陽活動頻繁期和地磁活動頻繁期內，異常主要集中在磁赤道附近和美洲扇區，這與電離層對地磁暴和太陽活動的主要響應特征一致，且本文方法在此期間并沒有探測到全球大范圍尺度的異常。結合前2個震例研究發現，PCA-滑動時窗方法探測異常不易受太陽活動和地磁活動的影響，因此本震例探測的02-02和02-12異常為正異常，且均能在磁共軛區被觀測到，該探測結果與前人研究異常分析的結果比較一致[3]。

(a)Dst指數；(b)KP指數；(c)F10.7指數(虛線是各指數設定的閾值，Dst為-30 nT，KP為3)；(d)基于PCA-滑動時窗法探測得到的二維異常分布(五角星為2014-02-12于田地震震中位置，紅色曲線為磁赤道線(磁緯0°)和磁緯±15°)；(e)利用PCA分析電離層TEC所得震中位置的主特征值序列(水平紅線為0.5閾值線，垂線紅線為地震發生日)

這3個震例驗證了PCA-滑動時窗方法探測電離層異常的有效性，并具有以下優勢：1)在數據處理過程中剔除掉了許多無關因素，計算速度快，探測到的異常情況比較直觀、簡潔；2)能進行長時間序列的時空異常探測，并能很好地區分正負異常；3)不易受地磁活動和太陽活動的影響。

2.2 2020年瓦哈卡地震研究

墨西哥位于三大構造板塊之上，地殼板塊的相對運動導致其頻繁發生地震和爆發火山。2020-06-23在墨西哥瓦哈卡附近發生了7.4級地震，目前還沒有關于該地震震前電離層異常的研究，因此本文利用PCA-滑動時窗方法對該地震的電離層TEC資料進行處理和分析，初探瓦哈卡地震震前電離層TEC的時空變化特性，為后續的地震研究工作提供可靠信息。

圖5為瓦哈卡地震震前21 d至地震后2 d共24 d的F10.7、Dst及KP等指數的變化情況，由圖可知，這24 d內空間天氣比較平靜。圖5(d)為異常探測結果，可以看出，06-03(震前20 d)TEC異常減少，06-10(震前13 d)TEC異常增加，06-11(震前12 d)、06-12(震前11 d)、06-14(震前9 d)及06-17～19(震前4～6 d)電離層TEC異常減少，但由于06-03在23.274°S、68.468°W發生了6.8級地震，當天異常可能與該地震有關。從空間上來看，地震異常位于震中偏南區域，主要表現為低值異常。

(a)Dst指數；(b)KP指數；(c)F10.7指數(虛線是各指數設定的閾值，Dst為-30 nT，KP為3)；(d)基于PCA-滑動時窗法探測得到的二維異常分布(五角星為2020-06-23瓦哈卡地震震中位置，紅色曲線為磁赤道線(磁緯0°)和磁緯±15°)；(e)利用PCA分析電離層TEC所得震中位置的主特征值序列(水平紅線為0.5閾值線，垂線紅線為地震發生日)

3 結 語

本文基于GIM數據，利用PCA-滑動時窗方法探測了4次地震的震前電離層異常，首先利用PCA-滑動時窗方法探測了3次強震事件(汶川地震、玉樹地震及于田地震)，并將結果與滑動時窗法和滑動四分位距法的結果進行了對比驗證，認為本文方法具有有效性及探測優勢。結果表明：1)PCA-滑動時窗方法能夠應用于地震電離層的時空異常探測，且探測結果直觀、簡潔；2)相較于PCA方法，PCA-滑動時窗方法能夠準確探測出正負異常，并能更好地應用于長時間序列的時空異常探測；3)相較于傳統方法，本文方法探測電離層異常時不易受太陽活動和地磁活動的影響，探測出的異常結果可信度高。本文在驗證了PCA-滑動時窗方法可行性的基礎上對2020年瓦哈卡地震進行了異常探測，結果表明，06-10(震前13 d)TEC異常增加，06-11(震前12 d)、06-12(震前11 d)、06-14(震前9 d)及06-17～19(震前4～6 d)電離層TEC異常減少，這些異常很可能是地震的前兆異常。從空間位置來看，異常位于震中偏南區域，TEC主要表現為低值異常。綜合發現，電離層TEC異常主要出現在震中偏南方向，且呈共軛結構。由于發震機制復雜，目前電離層與地震之間的耦合關系尚不清楚，因此將TEC異常作為一種前兆異常來進行地震預測非常重要。本文方法可很好地進行多震例異常特征(出現時間、空間方位、正負情況等)的提取及統計研究，從而更好地應用于電離層異常監測。