基于GPS數(shù)據(jù)的震前短臨異常檢測

李 南，林莉莉

（福建農(nóng)林大學(xué) 計(jì)算機(jī)與信息學(xué)院，福州 350000）

0 引 言

震前短臨異常檢測已成為防震減災(zāi)研究與應(yīng)用領(lǐng)域的研究熱點(diǎn)。當(dāng)前，國內(nèi)外學(xué)者大多采用射出長波輻射、電離層電子總含量等電、熱指標(biāo)，來進(jìn)行短臨異常檢測，以探測震前一段時間及一定區(qū)域內(nèi)可能發(fā)生的地球物理、化學(xué)變化。但這些方法普遍存在觀測數(shù)據(jù)不足、容易受到人為活動干擾等缺陷。

近年來，隨著全球?qū)Ш叫l(wèi)星系統(tǒng)技術(shù)（Global Navigation Satellite System，GNSS）的發(fā)展和GPS臺站的普及，利用GPS數(shù)據(jù)進(jìn)行短臨異常檢測已成為熱門研究方向。與傳統(tǒng)的電、熱指標(biāo)相比，利用GPS數(shù)據(jù)的方式，能更直接觀測到大地震發(fā)生前出現(xiàn)的地表中、低頻地形形變，具有較好的客觀性和穩(wěn)定性。但是，現(xiàn)有大多數(shù)研究方法普遍依賴地理學(xué)科領(lǐng)域?qū)＜业闹R和經(jīng)驗(yàn)，且僅用單個典型震例的GPS數(shù)據(jù)來驗(yàn)證異常檢測方法的有效性，使其存在主觀性較強(qiáng)、普適性較差等問題。

鞅理論作為現(xiàn)代概率和隨機(jī)過程的基礎(chǔ)，適用于時間序列數(shù)據(jù)分析場合，已被廣泛運(yùn)用于數(shù)據(jù)挖掘中的決策優(yōu)化、異常檢測等領(lǐng)域。因此，本文結(jié)合數(shù)據(jù)挖掘的相關(guān)知識，運(yùn)用鞅理論，提出一種基于GPS數(shù)據(jù)的震前短臨異常檢測算法（Anomaly Detection Algorithm based on GPS data，ADA）。

實(shí)驗(yàn)結(jié)果表明，ADA算法所識別的GPS數(shù)據(jù)中異常出現(xiàn)時間與地震發(fā)生時間存在顯著相關(guān)。相比于傳統(tǒng)的準(zhǔn)則分析方法、異常檢測模型ARIMA、單類別支持向量機(jī)OCSVM，以及基于兩階段聚類的異常檢測算法TSOD等，ADA算法能夠更直觀、準(zhǔn)確地反映震前GPS數(shù)據(jù)中出現(xiàn)的異常，可為地震預(yù)警減災(zāi)提供有效手段。

1 研究方法

本文短臨異常檢測算法包括：數(shù)據(jù)預(yù)處理、特征提取以及異常檢測3部分內(nèi)容。

1.1 數(shù)據(jù)預(yù)處理

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

首先，采用二階多項(xiàng)式擬合方法，依次對東西向和北南向上的GPS數(shù)據(jù)進(jìn)行缺失值填補(bǔ)。

1.2 特征提取

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

使用所有相關(guān)GPS臺站的數(shù)據(jù)，估算震中位置的每日坐標(biāo)（以東西向?yàn)槔?/p>

震中位置的經(jīng)度x，相關(guān)的個GPS臺站的經(jīng)度，，…，x，相鄰兩日各臺站東西向坐標(biāo)偏移量，，…，Δx。

震中東西向的每日坐標(biāo)。

針對各臺站每日坐標(biāo)位移數(shù)據(jù)(x，Δx)，1，2，…，，求解出使得二項(xiàng)式擬合函數(shù)的損失函數(shù)最小化時的權(quán)重。

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

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

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

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

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

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

使用線性回歸算法，計(jì)算第天前_天內(nèi)東西向、北南向坐標(biāo)：

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

計(jì)算滑動窗口內(nèi)，斜率S（）、S（）的變化范圍：

計(jì)算震中第天的綜合特征值：

1.3 異常檢測

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

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

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

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

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

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

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

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

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

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

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

使用某震中的綜合特征序列，計(jì)算該震中t天內(nèi)的冪鞅值M。

某震中的綜合特征{，，…，V}、停止參數(shù)、穩(wěn)定參數(shù)_。

某震中天內(nèi)的冪鞅值M。

設(shè)：_1。

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

根據(jù)C和D，采用公式（3）計(jì)算S。

根據(jù)S，采用公式（4）計(jì)算M。

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

1.4 ADA算法流程

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

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

2 實(shí)驗(yàn)設(shè)計(jì)與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)

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

綜合考慮GPS臺站位置和數(shù)據(jù)完整性，本文最終采用的震例數(shù)據(jù)見表1，相關(guān)信息來自美國地質(zhì)調(diào)查局網(wǎng)站（https：／／earthquake.usgs.gov／）。

表1 震例數(shù)據(jù)Tab.1 Earthquake data

2.2 對比算法與參數(shù)設(shè)置

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

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

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

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

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

2.3 性能對比分析

圖2～圖4、表2和圖5分別給出了準(zhǔn)則、ARIMA模型、OCSVM、TSOD算法以及ADA算法的運(yùn)行結(jié)果。

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

圖2 P500臺站各向kσ準(zhǔn)則分析結(jié)果Fig.2 Analysis result of kσmethod on P500 station

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

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

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

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

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

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

圖5 2010－04－04地震的ADA算法運(yùn)行結(jié)果Fig.5 Analysis result of ADA algorithm on 2010－04－04 earthquake

2.4 參數(shù)優(yōu)化分析

2.4.1 穩(wěn)定參數(shù)分析

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

圖6 基于不同穩(wěn)定參數(shù)的ADA算法結(jié)果Fig.6 Comparison of result with different stable_day

從圖6可看出，不同的_參數(shù)值，并不會對檢測結(jié)果造成太大影響。不同取值下，冪鞅值的變化趨勢均表現(xiàn)為在地震前較短的一段時間內(nèi)開始增加，并在地震后的一段時間內(nèi)達(dá)到波峰。這是由于當(dāng)?shù)貧み\(yùn)動相對穩(wěn)定時，前_天的GPS數(shù)據(jù)并不會發(fā)生太大變化，因而對結(jié)果的影響不大，不同取值下冪鞅值波峰出現(xiàn)的時間僅差距1～3天。但當(dāng)參數(shù)為5時，冪鞅值最早出現(xiàn)增加的趨勢。據(jù)此，實(shí)驗(yàn)中將穩(wěn)定參數(shù)_設(shè)置為5。

2.4.2 平滑窗口分析

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

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

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

2.4.3 停止參數(shù)分析

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

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

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

2.5 顯著性檢驗(yàn)

為了驗(yàn)證ADA算法所識別的GPS數(shù)據(jù)中存在的短臨異常與對應(yīng)地震之間存在關(guān)聯(lián)，本文使用Molchan圖表法，在表1所示的8個震例上進(jìn)行了統(tǒng)計(jì)顯著性檢驗(yàn)，結(jié)果如圖9所示。

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

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

3 結(jié)束語

震前短臨異常檢測是地震預(yù)警減災(zāi)的關(guān)鍵。本文提出的ADA算法能夠彌補(bǔ)現(xiàn)有方法存在的主觀性較強(qiáng)、普適性較差的問題。8個震例的實(shí)驗(yàn)結(jié)果證實(shí)了ADA算法檢測到的短臨異常與地震之間存在顯著相關(guān)。另外，本文也對算法的參數(shù)進(jìn)行了優(yōu)化分析。

然而，地震監(jiān)測預(yù)警是一項(xiàng)復(fù)雜的任務(wù)，會涉及與孕震相關(guān)的巖石圈－蓋層－大氣層－電離圈層等多個數(shù)據(jù)源。因此，如何結(jié)合這些異源數(shù)據(jù)來進(jìn)行異常檢測是下一步的研究方向。