李 南，林莉莉

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

0 引 言

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

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

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

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

1 研究方法

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

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

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

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

1.2 特征提取

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 異常檢測(cè)

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

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

D為V相對(duì)于C的偏移程度，即：

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

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

其中，是一個(gè)(0，1]之間的隨機(jī)數(shù)，（）是一個(gè)函數(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í)間序列數(shù)據(jù)的連續(xù)變化情況，使用統(tǒng)計(jì)量?jī)琪敝?，可?duì)持續(xù)一段時(shí)間內(nèi)數(shù)據(jù)的異常程度進(jìn)行量化。冪鞅值越高，越傾向于拒絕接受數(shù)據(jù)序列分布穩(wěn)定的假設(shè)。本文采用鞅理論對(duì)數(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值無限增大，需引入一個(gè)停止參數(shù)作為M的閾值。此外，本文還引入一個(gè)穩(wěn)定參數(shù)_，從第_1天開始計(jì)算冪鞅值，以避免過短的時(shí)間序列數(shù)據(jù)對(duì)分析結(jié)果造成誤差。異常檢測(cè)算法具體過程如算法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ù)

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

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

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

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

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

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

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

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

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

2.3 性能對(duì)比分析

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

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

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

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

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

圖4（a）～圖4（c）中，橫坐標(biāo)表示時(shí)間，縱坐標(biāo)表示當(dāng)天給定方向的坐標(biāo)值，在OCSVM算法下的類別。1表示正常，－1表示異常。從圖4（a）可看出，P500臺(tái)站東西向異常最早出現(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個(gè)方向上異常出現(xiàn)時(shí)間的不統(tǒng)一都增加了結(jié)果分析的難度。

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

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

圖5給出了2010－04－04地震ADA算法運(yùn)行結(jié)果（即冪鞅值的變化趨勢(shì)）。從圖5可明顯看出，在震前絕大部分時(shí)間，冪鞅值始終保持在一個(gè)相對(duì)較小的區(qū)間內(nèi)。由于大地震震前能量是一個(gè)累積的過程，冪鞅值從地震前較短的一段時(shí)間（約1個(gè)星期）開始緩慢增加，說明GPS數(shù)據(jù)開始出現(xiàn)異常，暗示震前局部應(yīng)力場(chǎng)開始調(diào)整。地震后各個(gè)臺(tái)站的坐標(biāo)發(fā)生了較大變化，因此冪鞅值的波峰是在地震后出現(xiàn)，且在地震后迅速超過預(yù)設(shè)的閾值。這說明ADA算法對(duì)2010－04－04地震的異常檢測(cè)是有效的，且比4種對(duì)比算法能更直觀地反映出震前短臨異常，不易出現(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ù)_對(duì)本文方法性能的影響，在2010－04－04地震上分別將_設(shè)置為5、7和9進(jìn)行實(shí)驗(yàn)，冪鞅值的變化趨勢(shì)如圖6所示。

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

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

2.4.2 平滑窗口分析

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

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

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

2.4.3 停止參數(shù)分析

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

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

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

2.5 顯著性檢驗(yàn)

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

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

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

3 結(jié)束語

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

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