孫 哲，劉智敏，杜自豪，孔曉宇

(山東科技大學 測繪科學與工程學院，山東 青島 266590)

電離層是一種由大量電子組成的殼層，容易受到各種自然現象的干擾，如火山爆發(fā)、太陽耀斑、地震等。研究表明地震孕育過程中，在巖石圈、大氣圈以及電離層確實存在電磁異常現象，且電離層前兆的出現有穩(wěn)定的時間尺度，這一突出特點將使它在短臨地震預報中更具實用價值[1]。隨著全球定位系統(tǒng)(Global Positioning System，GPS)的發(fā)展，基于GPS觀測數據反演獲取電子總含量(Total Electron Content，TEC)的成本逐漸降低，討論震前TEC異常變化的規(guī)律性成為熱點[2-6]。

探測數據異常性的方法有很多，傳統(tǒng)靜態(tài)探測法有中位數法、平均值法等，這些方法不能剔除數據中的異常數據，使得計算的背景參考值存在較大偏差，導致探測結果不準確，為了彌補傳統(tǒng)靜態(tài)探測法背景參考值偏差較大的問題，學者提出動態(tài)探測法—滑動時窗法[7]、滑動四分位法[8]。這兩種滑動時間窗口方法作為現在主流的動態(tài)探測法，雖然較傳統(tǒng)靜態(tài)探測法提高了背景參考值的精度，但對主要異常時段不夠突出且探測時窗較長。為提高對異常時段的敏感度和縮短時間窗口長度，本文提出利用平滑Z分法進行震前分析并對2017年4月29日4時23分發(fā)生的菲律賓棉蘭老島7.2級地震(震中位于5.51°N，125.08°E，震源深度為57 km)進行震前異常探測，探測結果與滑動時窗法和滑動四分位距法的探測結果基本一致，結合太陽地磁場活動指數、TEC異?？臻g變化后推測3種方法探測到的異常與此次地震有關，進一步證實了平滑Z分法的可行性、有效性。

1 數據獲取及異常檢測方法

1.1 數據獲取

本文在分析電離層異常時，采用的電離層數據來源于IGS提供的全球電離層格網數據，時間分辨率為 2 h，空間分辨率為 5°×2.5°(經度×緯度)，該數據是通過全球分布的IGS基準站觀測資料計算得到的。根據文獻[9]提出的孕震區(qū)半徑公式R=100.43M(M代表震級，R代表孕震區(qū)半徑(km))，可以得出菲律賓棉蘭老島地震孕震區(qū)半徑約為1 247 km，全球電離層格網數據的空間分辨率可以滿足此次地震電離層異常擾動分析的需求。地震前電離層異常區(qū)域往往偏離震中區(qū)域一定距離[10]，利用格網數據選取震中附近4個格網點A(125°E，5°N)，B(130°E，5°N)，C(125°E，7.5°N)，D(130°E，7.5°N)震前15 d和震后5 d的TEC時間序列進行分析效果會更好。震中以及4個格網點分布情況如圖1所示，其中紅色邊框五角星表示震中點(紅色實心五角星表示震中點在亞洲板塊的位置)，紅色圓表示格網點。

圖1 震中及異常觀測點示意圖

由于太陽活動和地磁場的變化都會導致電離層的異常變化，在分析電離層異常時需結合該時期的空間環(huán)境狀況綜合分析該段時期的電離層異常變化原因。本文提取了震前15 d及震后5 d的太陽射電量F10.7指數、地磁Dst指數、地磁Kp指數的時間序列，各指數的評判標準已具備[11-12]。

1.2 異常探測方法

1.2.1 滑動時窗法

滑動時窗法是以滑動均值為背景參考值，滑動標準差為異常檢測依據的探測方法，對于待探測的時間序列取l天作為時窗窗長，本文選擇15 d作為時間窗口來計算均值和均方差[13]。

均值：

(1)

均方差：

(2)

在進行異常探測處理時，通常以倍均方差為探測依據，這里取2倍均方差，置信度超過95%，探測上下限：

(3)

待探測天數同一時刻若超出探測上下限值，即可視該天該時刻的TEC值為異常。若探測第i+1天某一時刻的異常，需計算前一天同一時刻的平均數及均方差，作為背景參考值，按照探測上下限對第i+1天進行判斷。

1.2.2 滑動四分位法

結合以往的經驗，綜合考慮均值法和中位數法的不足，臺灣國立中央大學的劉正彥教授提出四分位距法[2]。四分位距( Inter Quartile Rang，IQR) 是一種穩(wěn)健統(tǒng)計技術中用于表示數據離散度的一個量，常用來檢查數據的異常情況。在選擇時間窗口時長度不宜過長，因此本文中選擇15 d作為時間窗口，取其前15 d同一時刻的TEC組成一組時間序列，并按照從小到大順序對其進行排列，為x1，x2，…，x15，則：

(4)

(5)

式中：Q1為上四分位，Q2為中位數，Q3為下四分位，k為常倍數，IQR為四分位距，LB為下邊界，UB為上邊界。本文定義Q2±1.5IQR作為TEC異常檢驗上下邊界的判斷閥值，約為2倍均方差，置信度在95%左右。如果TEC值沒有超出上下界則ΔTEC=0，超出上界ΔTEC>0，表示正異常，低于下界值ΔTEC<0，表示負異常。

1.2.3 平滑Z分法

針對實時數據中的峰值信號檢測，外國學者Jean-Paul構建了一個基于分散原則的穩(wěn)健算法，稱為平滑Z分法(Smoothed Z-score Method)。該算法使用探測時窗掃描一組數據點并計算移動平均值和標準差，高于設定閾值的Z分會被檢測出來，Z分的方程：

(6)

平滑Z分法需要設定3個參數：探測時窗(l)，閾值(th)和影響因子(in)。影響因子設置在0～1之間，表示以前觀測點對新觀測點的影響比。平滑Z分法的組成部分：

(7)

(8)

(9)

式中：xi為檢測到的異常值；si為平滑后的值?？紤]本文數據特點，參考Mehmet等學者的研究，本文選取影響因子in=0.5，閾值th=2，時間窗口l=4 d，提取震前15 d和震后5 d每個時間節(jié)點(全球格網數據2 h為一個時間節(jié)點)的TEC值與對應前15 d同一時間節(jié)點TEC序列的中值做差得到趨勢序列M，將趨勢序列M作為觀測值進行異常檢測。異常判定從觀測值第一個數據點開始，將前4 d的趨勢序列作為背景參考值，檢測到異常值進行平滑后繼續(xù)向前判定，直至觀測值判定完畢。

2 分析與討論

2.1 時間序列異常分析

本文提取震前15 d到震后5 d的TEC時間序列作為分析對象，采用上述3種方法對A，B，C，D4個點進行TEC異常探測，結果如圖2—圖5所示?；瑒铀姆治环ńY果顯示4個格網點大部分時間TEC時間序列都在上下邊界之內，從震前10 d開始TEC異常出現并逐漸增大，在震前9 d TEC異常達到最大，震前8 d逐漸縮小，震前7 d達到第二個峰值，震前6 d逐漸縮小，震前5 d到震后4 d之間縮小速度增快，震后5 d之后異常消失，4個格網點表現出一致性，整體呈現雙峰結構，震前9 d在A點出現最大正異常約為9TECu，震后3 d在C點出現最大負異常約為-4TECu，A，C點異常值比B，D點高；滑動時窗法結果與滑動四分位法結果大致相同，4個格網點在震前10～6 d期間呈現出雙峰結構，震前9 d在A點達到最大正異常約5TECu，震后3 d在A點達到最大負異常約-1TECu，相比于滑動四分位法，檢測到的TEC異常節(jié)點相對較少，最大正異常和最大負異常比滑動四分位法都低，A，C點異常值高出B，D點程度和對異常時段的敏感度比滑動四分位法都低；平滑Z分法結果與前兩種方法大致相同，震前9 d和7 d出現兩個峰值，其中震前9 d在A點達到最大正異常約11TECu，震后3 d在C點達到最大負異常約-3TECu，相對于滑動四分位法和滑動時窗法異常程度更高，自身A，C點異常值較B，D異常值高，與前兩種方法一致，探測到的異常節(jié)點較前兩種方法相對較少，但震前9 d和震前7 d的異常節(jié)點占比提高，更加突出異常時段。觀察圖1可知A，C點比B，D距震中更近，參考3種方法的異常結果推測格網點異常值與距震中距離相關；圖2—圖5顯示震前9 d和7 d兩個異常峰值最為特殊，是否與此次地震相關還需進一步結合地磁活動指數和空間異常進一步分析。

圖2 2017-04-29菲律賓棉蘭老島地震A點(125°E，5°N)TEC時間序列及滑動四分位法、滑動時窗法和平滑Z分法法結果(圖中0點表示地震當天)

圖4 2017-04-29菲律賓棉蘭老島地震C點(125°E，7°N)TEC時間序列及滑動四分位法、滑動時窗法和平滑Z分法法結果(圖中0點表示地震當天)

圖5 2017-04-29菲律賓棉蘭老島地震D點(130°E，7°N)TEC時間序列及滑動四分位法、滑動時窗法和平滑Z分法結果(圖中0點表示地震當天)

結合太陽活動方面因素選取正異常最明顯的震前9 d和7 d來進一步分析震前TEC異常原因，本文提取了震前15 d及震后5 d太陽地磁活動時間序列數據如圖6所示。從圖中可以看出，震前9 d和7 d中只有F10.7指數在設定閾值以內，其余各項指數均有所超標：震前7 d大部分時間的Dst指數超標，最大值達到6，Kp指數也超標，最低值接近-50 nT，說明這幾天的太陽地磁活動較強烈，推測圖2—圖5中用滑動四分位法、滑動時窗法和平滑Z分法探測出震前7 d的TEC異常與太陽地磁活動干擾密切相關。

圖6 2017-04-29菲律賓棉蘭老島地震地磁活動時間序列圖(圖中0點表示地震當天)

震前9 d出現很短時間的小幅度太陽活動異常，在UTC 08:00—09：00時Dst≈-40 nT，Kp=4，其余時間Dst≥-30 nT，Kp<4，F10.7全部低于100SFU，可見這天大部分時間太陽活動很平靜，但上述3種方法均檢測到這天UTC 00：00—14：00出現大量正異常，是否由孕育區(qū)地震引起還需從空間角度進一步研究。

2.2 空間分布異常分析

區(qū)分地震擾動與磁暴等活動的一個很重要的評判標準為電離層TEC異常是否具有局地性特征[14]。為了進一步分析震前9 d和震前7 dTEC異常是否與此次地震有關，需要進行空間上大范圍的異常檢測，為盡可能多的探測異常時段，突出空間表現效果，本文分別給出滑動四分位法在震前9 d和震前7 d，90°E～160°E，10°S～30°N 范圍內的電離層異常圖。

震前9 d UTC 08：00—18：00間隔2 h的TEC異常情況如圖7所示，其中紅色五角星為震中位置，由圖7可知在UTC 08：00時，震中東南部出現了明顯的正異常，異常幅度在6～8TECu；隨后在UTC 10：00正異常向震中移動并且異常值逐漸增大，達到峰值，最大異常值約10TECu；UTC 12：00—14：00時正異常開始削減并向西北方移動；UTC 16：00—18：00正異常逐漸削弱直至消失。這一天TEC異?，F象在UTC 08：00持續(xù)到UTC 16：00，時間較長，異常大部分出現在震中區(qū)域附近，而在其他地方沒有出現明顯的異常，推測震前9 d的電離層異常由此次地震引起。

圖7 震前9 d(2017-04-20)TEC異常分布圖(紅色五角星為震中位置)

震前7 d UTC 02：00—12：00間隔2 h的TEC異常情況如圖8所示，其中紅色五角星為震中位置，由圖8可知在UTC 04：00之前，整個區(qū)域中絕大部分的電離層很平靜，只有在震中東南方向出現了微弱的電離層異常，異常值為5TECu左右。UTC 04：00在震中東西部開始出現電離層異常，此后，異常區(qū)域瞬間向外擴散且異常值也顯著增加，UTC 08：00電離層異常達到峰值為12TECu左右，此時異常區(qū)域覆蓋了整個分析區(qū)域的絕大部分。UTC 10：00之后，電離層異常迅速減弱消失，UTC 12：00整個分析區(qū)域只有很小一部分微弱異常。此次電離層異常較震前9 d強，覆蓋范圍廣，異常時間縮短，結合圖6，可知震前7 d地磁活動劇烈，推測震前7 d電離層異常與地磁活動有關。

圖8 震前7 d(2017-04-22)TEC異常分布圖(紅色五角星為震中位置)

3 結束語

利用IGS提供的全球電離層格網數據對2017年菲律賓棉蘭老島地震分別用滑動四分位法、滑動時窗法和平滑Z分法進行震前TEC異常探測，主要結果如下：

1)3種方法均能探測到震前的TEC值異常信息。3種方法在震前5～10 d 4個點均呈現出明顯的正異常，震前9 d和7 d均出現峰值，震前9 d在A點均出現最大正異常，探測結果顯示3種方法在4個點A，B，C，D中表現一致，格網點異常強度與距震中距離有關，震前9 d的TEC值異常與地震相關；平滑Z分法較傳統(tǒng)的動態(tài)測量方法探測到的異常天數較少，但與地震相關異常天數占比更高，對異常數據進行平滑處理，提高對異常時段的敏感度，探測時窗上更短為4 d，有利于短臨地震的實時預報。

2)在空間異常分布的圖像上來看，震前9 d正異常在震中東南方出現逐漸向震中靠攏并且異常值逐漸增大，隨后逐漸向震中西北方散去，異常值也逐漸削減直至消失，期間異常值絕大部分時間都集中在震中附近，具有區(qū)域性；震前7 d異常在震中附近突然出現，增加速度極快，隨后迅速消失，移動方向無規(guī)律，分布范圍覆蓋整個研究區(qū)，具有全域性。

通過對比3種方法的異常探測結果并結合空間異常分布和太陽地磁活動排除干擾因素，推測震前9 d的正異常與此次地震有關，震前7 d的正異常與太陽活動有關。此次探測結果證實了平滑Z分法的可行性和較傳統(tǒng)動態(tài)測量方法的優(yōu)越性，將幾種方法互為補充可突出主要異常時段，對于短臨地震預報有一定價值。