基于GPS數據的震前短臨異常檢測

李 南，林莉莉

（福建農林大學 計算機與信息學院，福州 350000）

0 引 言

震前短臨異常檢測已成為防震減災研究與應用領域的研究熱點。當前，國內外學者大多采用射出長波輻射、電離層電子總含量等電、熱指標，來進行短臨異常檢測，以探測震前一段時間及一定區域內可能發生的地球物理、化學變化。但這些方法普遍存在觀測數據不足、容易受到人為活動干擾等缺陷。

近年來，隨著全球導航衛星系統技術（Global Navigation Satellite System，GNSS）的發展和GPS臺站的普及，利用GPS數據進行短臨異常檢測已成為熱門研究方向。與傳統的電、熱指標相比，利用GPS數據的方式，能更直接觀測到大地震發生前出現的地表中、低頻地形形變，具有較好的客觀性和穩定性。但是，現有大多數研究方法普遍依賴地理學科領域專家的知識和經驗，且僅用單個典型震例的GPS數據來驗證異常檢測方法的有效性，使其存在主觀性較強、普適性較差等問題。

鞅理論作為現代概率和隨機過程的基礎，適用于時間序列數據分析場合，已被廣泛運用于數據挖掘中的決策優化、異常檢測等領域。因此，本文結合數據挖掘的相關知識，運用鞅理論，提出一種基于GPS數據的震前短臨異常檢測算法（Anomaly Detection Algorithm based on GPS data，ADA）。

實驗結果表明，ADA算法所識別的GPS數據中異常出現時間與地震發生時間存在顯著相關。相比于傳統的準則分析方法、異常檢測模型ARIMA、單類別支持向量機OCSVM，以及基于兩階段聚類的異常檢測算法TSOD等，ADA算法能夠更直觀、準確地反映震前GPS數據中出現的異常，可為地震預警減災提供有效手段。

1 研究方法

本文短臨異常檢測算法包括：數據預處理、特征提取以及異常檢測3部分內容。

1.1 數據預處理

由于數據采集設備、傳輸線路故障等原因，各GPS臺站的原始數據存在部分數據缺失的情況。另外，GPS臺站每日坐標包括東西、北南和垂直3個方向的數據，但垂直向數據通常誤差較大。因此，本文僅針對各GPS臺站東西向和北南向的坐標數據進行處理。

首先，采用二階多項式擬合方法，依次對東西向和北南向上的GPS數據進行缺失值填補。

1.2 特征提取

當同一個震例涉及多個GPS臺站，且不同臺站之間GPS數據出現異常的時間和強度存在較大差異時，會導致異常檢測結果出現較大偏差。為了彌補以上不足，本文基于同一震例的所有相關GPS臺站數據，采用二階多項式擬合方法估算相應震例震中位置的每日坐標。震中坐標估算過程如算法1所示。若某一震例只涉及一個臺站，則直接使用該臺站數據即可。

使用所有相關GPS臺站的數據，估算震中位置的每日坐標（以東西向為例）。

震中位置的經度x，相關的個GPS臺站的經度，，…，x，相鄰兩日各臺站東西向坐標偏移量，，…，Δx。

震中東西向的每日坐標。

針對各臺站每日坐標位移數據(x，Δx)，1，2，…，，求解出使得二項式擬合函數的損失函數最小化時的權重。

基于擬合函數(，)，輸入震中位置的經度x，獲得估算的偏移量( x，)。

根據x和(x，)，得到預估的震中東西向的坐標x＋y（x，）。

為了降低GPS數據中白噪聲、高斯噪聲等對檢測結果的影響，在運用算法1獲得震中位置各方向的時序坐標數據后，使用滑動窗口技術對數據進行降噪處理。固定大小的滑動窗口內樣本數據斜率的變化，不僅能有效刻畫數據在長趨勢變化下的短期特征，而且對噪聲具有一定的魯棒性。因此，本文使用GPS臺站東西向和北南向坐標數據的斜率變化范圍來提取震中每日的綜合特征。特征提取過程如算法2所示。

根據震中東西向、北南向的每日坐標（見算法1），提取震中的每日綜合特征。

震中東西向、北南向的每日坐標（）、（），滑動窗口大小_。

震中第天的綜合特征（）。

使用線性回歸算法，計算第天前_天內東西向、北南向坐標：

（）：_≤≤以及（）：_≤≤的斜率，記為S（）、S（）。

計算滑動窗口內，斜率S（）、S（）的變化范圍：

計算震中第天的綜合特征值：

1.3 異常檢測

地震的發生通常需要一定時間的能量累積，本文基于提取到的某震中綜合特征值，評估該震中第天的短臨異常程度，并在此基礎上利用鞅理論評估震中在某連續時間段內短臨異常程度。

設：C為前1天綜合特征{，，…，V}的均值（即中心值），即：

D為V相對于C的偏移程度，即：

其中，‖·‖表示歐式距離。

根據公式（2）得到的偏移程度，進一步計算V和{，，…，V}之間的相異度值S，即：

其中，是一個(0，1]之間的隨機數，（）是一個函數，返回滿足指定條件數據的數量。

如：( j｜D＞D)表示在｛，，…，V｝中D＜D，1，2，…，的數據數量。

從公式（3）可以看出，S∈(0，1]。根據同一分布中各樣本差異最小化原則，S越小V就越遠離前1天數據的中心值C，則V和{，，…，V}之間越不相似，表明震中第天的短臨異常程度越高。

鞅理論適合于刻畫時間序列數據的連續變化情況，使用統計量冪鞅值，可對持續一段時間內數據的異常程度進行量化。冪鞅值越高，越傾向于拒絕接受數據序列分布穩定的假設。本文采用鞅理論對數據{，，…，S}的分布情況進行量化分析，得到天內{，，…，S}的冪鞅值M。

其中，S為V和{，，…，V}之間的相異度值，根據文獻［1］的推論取值0.82。

從公式（4）可見，冪鞅值M值越大，說明天內頻繁出現值較小的情況，暗示t天內GPS數據頻繁出現異常的程度越高。為了避免公式（4）中冪鞅值M值無限增大，需引入一個停止參數作為M的閾值。此外，本文還引入一個穩定參數_，從第_1天開始計算冪鞅值，以避免過短的時間序列數據對分析結果造成誤差。異常檢測算法具體過程如算法3所示。

使用某震中的綜合特征序列，計算該震中t天內的冪鞅值M。

某震中的綜合特征{，，…，V}、停止參數、穩定參數_。

某震中天內的冪鞅值M。

設：_1。

根據{，，…，V}，采用公式（1）、（2）分別計算，得到C和D。

根據C和D，采用公式（3）計算S。

根據S，采用公式（4）計算M。

如果M≤，則1，重新執行步驟25，否則將第1天作為第1天，重新執行本算法。

1.4 ADA算法流程

基于算法1、算法2和算法3，則ADA算法具體流程如圖1所示。

圖1 ADA算法流程圖Fig.1 Flow chart of ADA algorithm

2 實驗設計與分析

2.1 實驗數據

本文研究對象為2001～2010年間，北美發生的震源深度小于60 km且震級大于6.0級的地震。GPS臺站時序坐標數據來自Nevada Geodetic Laboratory提供的數據共享服務網站（http：／／geodesy.unr.edu／）。選擇的GPS臺站，需處于受相應地震孕育影響的范圍（即震中半徑10之內，表示地震震級，的單位是km）。實驗選擇位于影響范圍內，最靠近震中的10個GPS臺站。由于地震孕育過程通常在地震前1～30天開始，因此所使用臺站的數據從地震發生前180天開始，到后30天結束。為了確保有足夠的臺站以供分析，單個臺站在這段時間內最多允許5%的數據缺失。從平臺獲得的數據是初步處理后的GPS臺站每日坐標分別是東西向、北南向和垂直向。文獻［3］中證實，GPS臺站的時序坐標數據在垂直方向的測量誤差遠大于水平方向。因此，實驗中只選用東西向、北南向的每日坐標作為研究數據，以保證分析結果的可靠性。

綜合考慮GPS臺站位置和數據完整性，本文最終采用的震例數據見表1，相關信息來自美國地質調查局網站（https：／／earthquake.usgs.gov／）。

表1 震例數據Tab.1 Earthquake data

2.2 對比算法與參數設置

為了驗證基于GPS數據的短臨異常檢測算法的有效性，將本文的ADA算法與傳統的準則分析方法、異常檢測模型ARIMA、單類別支持向量機OCSVM以及基于兩階段聚類的異常檢測算法TSOD進行性能對比。

（2）ARIMA模型：利用差分整合移動平均自回歸模型，得到一個預測值，通過預測值與實際值的誤差大小來判斷異常位置。本文ARIMA模型中的信息準則函數選用貝葉斯信息準則。

（3）OCSVM：單類別支持向量機，將異常檢測視為特殊的分類問題。在訓練過程中，只有一類數據，首先得到可以代表這部分數據的模型。在檢測過程中，判斷給定樣本是否屬于此類別。本文OCSVM算法中的核函數選用高斯核函數。

（4）TSOD算法：是基于兩階段聚類的多變量時間序列異常檢測算法。第一次聚類在各個變量上篩選初始異常時間，第二次聚類結合所有變量進行異常定位，以降低誤檢率。TSOD算法中第一次聚類使用基于混合高斯模型的EM算法，第二次聚類使用以全連接方式度量的層次聚類方法。

對比實驗和參數優化實驗中，使用的數據來自表1中震級最大的2010－04－04地震，準則、ARIMA模型、OCSVM以及TSOD算法僅使用距離地震震中最近的編號為P500的單個GPS臺站（Latitude：32.69N，Longitude：-115.30E）數據。本文ADA算法則涉及多個臺站的GPS數據，停止參數設置為2 000，窗口大小_設置為7，穩定參數_設置為5。

2.3 性能對比分析

圖2～圖4、表2和圖5分別給出了準則、ARIMA模型、OCSVM、TSOD算法以及ADA算法的運行結果。

表2 P500臺站使用TSOD算法分析結果Tab.2 Analysis result of TSOD algorithm on P500 station

圖2 P500臺站各向kσ準則分析結果Fig.2 Analysis result of kσmethod on P500 station

從圖3（a）～圖3（b）可見，震前45天左右，P500臺站東西向和北南向的坐標出現了ARIMA模型的預測值與真實值誤差較大（即異常）的情況，但隨著地震的臨近，誤差并沒有繼續保持在較高的水平。因此，不能完全確定此次異常是否與地震相關，也可能與噪聲有關。在垂直向上，ARIMA模型的預測值和真實值之間的誤差并沒有顯現出任何規律，這同準則的分析結果相一致。

圖3 P500臺站的ARIMA模型分析結果Fig.3 Analysis result of ARIMA model on P500 station

圖4（a）～圖4（c）中，橫坐標表示時間，縱坐標表示當天給定方向的坐標值，在OCSVM算法下的類別。1表示正常，－1表示異常。從圖4（a）可看出，P500臺站東西向異常最早出現在震前45天左右，在一周后斷斷續續出現并延續到震前。從圖4（b）～圖4（c）可看出，OCSVM算法在北南向和垂直向上，震前沒有發現明顯的異常顯現規律，并出現多次誤報。異常的不持續以及3個方向上異常出現時間的不統一都增加了結果分析的難度。

圖4 P500臺站的OCSVM算法分析結果Fig.4 Analysis result of OCSVM algorithm on P500 station

從表2可以看出，TSOD算法在P500臺站GPS數據上，最終檢測出3次異常。其中，距離地震發生最近的異常是在震前40天左右（2010－02－25），并出現了兩次明顯的誤報（在2009－12－22以及2010－01－28）。

圖5給出了2010－04－04地震ADA算法運行結果（即冪鞅值的變化趨勢）。從圖5可明顯看出，在震前絕大部分時間，冪鞅值始終保持在一個相對較小的區間內。由于大地震震前能量是一個累積的過程，冪鞅值從地震前較短的一段時間（約1個星期）開始緩慢增加，說明GPS數據開始出現異常，暗示震前局部應力場開始調整。地震后各個臺站的坐標發生了較大變化，因此冪鞅值的波峰是在地震后出現，且在地震后迅速超過預設的閾值。這說明ADA算法對2010－04－04地震的異常檢測是有效的，且比4種對比算法能更直觀地反映出震前短臨異常，不易出現誤報的情況。

圖5 2010－04－04地震的ADA算法運行結果Fig.5 Analysis result of ADA algorithm on 2010－04－04 earthquake

2.4 參數優化分析

2.4.1 穩定參數分析

為了分析穩定參數_對本文方法性能的影響，在2010－04－04地震上分別將_設置為5、7和9進行實驗，冪鞅值的變化趨勢如圖6所示。

圖6 基于不同穩定參數的ADA算法結果Fig.6 Comparison of result with different stable_day

從圖6可看出，不同的_參數值，并不會對檢測結果造成太大影響。不同取值下，冪鞅值的變化趨勢均表現為在地震前較短的一段時間內開始增加，并在地震后的一段時間內達到波峰。這是由于當地殼運動相對穩定時，前_天的GPS數據并不會發生太大變化，因而對結果的影響不大，不同取值下冪鞅值波峰出現的時間僅差距1～3天。但當參數為5時，冪鞅值最早出現增加的趨勢。據此，實驗中將穩定參數_設置為5。

2.4.2 平滑窗口分析

為了分析平滑窗口大小_對本文方法性能的影響，在2010－04－04地震上分別將_設置為5、7和10進行了實驗，冪鞅值的變化趨勢如圖7所示。

圖7 基于不同平滑窗口的ADA算法結果Fig.7 Comparison of result with different window_size

當平滑窗口_的取值較小時，在特征提取階段，計算第天的綜合特征值所需要的樣本數就越少，因此更容易受到單個樣本的影響，對異常的檢測也更敏感。從圖7可看出，對比_7和_10，當_5時，在地震發生前一個月就出現了冪鞅值緩慢增加的趨勢，相應冪鞅值的波峰也在地震發生前最早出現。而對于_7和_10，冪鞅值開始增加和波峰出現的時間并不存在顯著差別。因此，為了提高算法的魯棒性，實驗中將_設置為7更為合理。

2.4.3 停止參數分析

為了分析停止參數對本文方法性能的影響，在2009－07－02地震上分別將設置為500、1 000和2 000進行了實驗，冪鞅值的變化趨勢如圖8所示。

圖8 基于不同停止參數的ADA算法結果Fig.8 Comparison of result with different h

從圖8可看出，對于2009－07－02地震，3種停止參數設置下，冪鞅值均在地震發生前較短一段時間內顯著提高，這與文獻［1］中的結論相一致，即孕震活動最早在震前30天左右開始，并在震前幾天內表現最為活躍。值得注意的是，當參數設置為較小（500）時，冪鞅值的波峰多次出現，會導致異常的誤報。當檢測到GPS數據出現震前異常時，冪鞅值僅經過一天就從小于1 000增大到2 000以上。因此，參數1 000和2 000的冪鞅值曲線幾乎重合。由于較大的停止參數會減少預警時間，降低誤報可能，實驗中將設置為2 000。

2.5 顯著性檢驗

為了驗證ADA算法所識別的GPS數據中存在的短臨異常與對應地震之間存在關聯，本文使用Molchan圖表法，在表1所示的8個震例上進行了統計顯著性檢驗，結果如圖9所示。

圖9 8個震例的Molchan圖表分析結果Fig.9 Analysis result of ADA algorithm oneight earthquakes by Molchan error diagram

圖9中，橫坐標表示時間占有率，縱坐標表示相應的漏報率，使用的方法相比于隨機預測的優劣程度以曲線與圖表邊界線所包圍的面積來衡量，面積越小則說明預測效果越好。若測試的結果接近于圖9所示的對角線，則表示預測方法無統計顯著性。實驗中，若冪鞅值波峰出現在第天，那么將前天和后天作為變量，即以［，］作為預警時間范圍。若地震發生在此時間段內，則表示預警成功。通過調整和的取值以繪制圖表。從圖9中可以看出，ADA算法的時間占有率－漏報率曲線遠在對角線之下，說明所識別的短臨異常與對應地震之間存在顯著性關聯。

3 結束語

震前短臨異常檢測是地震預警減災的關鍵。本文提出的ADA算法能夠彌補現有方法存在的主觀性較強、普適性較差的問題。8個震例的實驗結果證實了ADA算法檢測到的短臨異常與地震之間存在顯著相關。另外，本文也對算法的參數進行了優化分析。

然而，地震監測預警是一項復雜的任務，會涉及與孕震相關的巖石圈－蓋層－大氣層－電離圈層等多個數據源。因此，如何結合這些異源數據來進行異常檢測是下一步的研究方向。