顧及分段相似性的同震形變時間序列恒星日濾波優(yōu)化算法

沈 飛，李建成，張小紅，束蟬方

1.武漢大學(xué) 測繪學(xué)院，湖北 武漢 430079；2.江蘇北斗衛(wèi)星應(yīng)用產(chǎn)業(yè)研究院，江蘇 南京 210032；3.南京工業(yè)大學(xué) 測繪學(xué)院，江蘇 南京 210009

1 引 言

隨著GPS接收機技術(shù)以及數(shù)據(jù)存儲技術(shù)的不斷發(fā)展，高頻（1Hz）和超高頻（20～50Hz）GPS接收機相繼出現(xiàn)，同時單歷元GPS處理技術(shù)亦逐漸成熟。GPS接收機被許多國內(nèi)外學(xué)者作為一種“地震儀”來實時監(jiān)測地殼運動，因此，一門新的前沿交叉學(xué)科——GPS地震學(xué)成為當前地震監(jiān)測研究的熱點［1－6］。單歷元高頻GPS定位精度受跟蹤GPS衛(wèi)星數(shù)及分布情況、軌道誤差、GPS和接收機鐘差、大氣延遲、天線相位偏差、多路徑效應(yīng)等眾多因素影響［7－8］，其中很多影響因素都與GPS衛(wèi)星—接收機的幾何關(guān)系有關(guān)，因而可利用精確的GPS衛(wèi)星運動周期來有效消除和減弱這些誤差的影響。因GPS衛(wèi)星重復(fù)周期約為一個恒星日，此濾波方法也稱為恒星日濾波（sidereal filtering）。國內(nèi)外學(xué)者對恒星日濾波及其應(yīng)用進行大量研究，其中主要出發(fā)點是確定準確的GPS衛(wèi)星重復(fù)周期，進而來削弱單歷元定位噪聲及分離多路徑誤差［9－14］。

恒星日濾波的基本思想就是假設(shè)地震前后兩天經(jīng)過GPS衛(wèi)星周期平移后的坐標時間序列的相關(guān)誤差影響是相同的，理想情況下，震前、震中、震后3天數(shù)據(jù)的解算過程中，各時刻應(yīng)采用相同的GPS衛(wèi)星，相同的解算策略，解算相同的參數(shù)個數(shù)。但現(xiàn)實情況往往難以滿足此條件，即使任何一顆GPS衛(wèi)星的模糊度解算出現(xiàn)問題都可能導(dǎo)致最終的濾波出現(xiàn)較大的偏差；更極端的情況是，發(fā)生強震后由于地表建筑物等發(fā)生的巨大變化，地震前后兩天的坐標時間序列大相徑庭。這時仍然采用常規(guī)的恒星日濾波，非但無法提高精度，反而還可能會帶來新的誤差。文獻［13］曾提出通過預(yù)處理觀測數(shù)據(jù)等方法解決上述問題，但其計算工作量將十分巨大，而且極端的強震影響通過預(yù)處理也無法解決。

針對這一問題，本文提出一種顧及時間序列分段相似性的恒星日濾波優(yōu)化算法，其基本思想是在構(gòu)建濾波殘差序列之前，首先確定用于構(gòu)建濾波的坐標殘差時間序列與地震當天坐標殘差時間序列的相似性，根據(jù)相似性度量指標的大小確定權(quán)重后再進行濾波，濾波時設(shè)定一定的閾值，當相似度量指標低于該閾值時，該天的坐標殘差時間序列將不參與構(gòu)建濾波時間序列。

2 同震坐標時間序列的相似性分析

2.1 時間序列相似性定義

時間序列相似性定義為［15］：給定兩個時間序列Q＝｛q1，q2，…，qn｝和C＝ ｛c1，c2，…，cn｝，如果dist（Q，C）≤α，就認為Q和C是相似的，其中dist（Q，C）是距離度量函數(shù)，比如歐氏距離；α是相似值閾值，用來調(diào)節(jié)相似程度。

時間序列的相似性分析是時間序列數(shù)據(jù)挖掘的基礎(chǔ)問題，用于分析不同數(shù)據(jù)段之間的相似性和差異性，提取出有價值的信息。幾乎所有的時間序列數(shù)據(jù)挖掘方法都和相似性分析有關(guān)，比如相似性搜索、聚類、分類、異常檢測、分段和主題發(fā)現(xiàn)等。

2.2 時間序列相似性度量指標

衡量兩個時間序列之間的相似性一般通過距離或者相似性度量函數(shù)來衡量，相似性與距離成反比關(guān)系，距離越近越相似，否則越不相似［16］。眾多學(xué)者已從不同的角度提出了多種距離或者相似性度量函數(shù)，距離函數(shù)主要有歐氏距離、動態(tài)彎曲距離、編輯距離等，相似性度量函數(shù)主要有余弦函數(shù)相似度、皮爾遜相關(guān)系數(shù)、通用影像質(zhì)量評價指數(shù)。本文介紹其中比較常用的3個相似性度量指標。

2.2.1 歐氏距離（Euclidean distance，ED）

歐氏距離［17］是時間序列相似性研究中最為廣泛采用的相似性度量，其定義為

歐氏距離越短，表示序列之間越相似。但其最大的不足之處就是對時間軸上的偏移變化非常敏感，在進行相似性時間序列搜索時，時間軸上微小的偏差往往會導(dǎo)致搜索結(jié)果的錯誤。

2.2.2 皮爾遜相關(guān)系數(shù)（Pearson correlation coefficient，PCC）

亦稱簡單相關(guān)系數(shù)，是統(tǒng)計學(xué)中常說的相關(guān)系數(shù)，它描述的是兩個定距變量間聯(lián)系的緊密程度。其定義為

2.2.3 通用影像質(zhì)量評價指數(shù)（universal image quality index，UIQI）

UIQI主要應(yīng)用于影像質(zhì)量評價，可認為是對皮爾遜相關(guān)系數(shù)的改進，其定義為［18］

此式可分解為

從3個相似性度量的定義可以看出，歐氏距離越小表明相似性越高，而皮爾遜相關(guān)系數(shù)和UIQI絕對值越大相似性越高。

2.3 時間序列分段相似性的確定

在研究地震同震形變時，其一般的研究對象是GPS站點坐標殘差的時間序列，即站點坐標時間序列減掉坐標參考值后的殘差序列［19－20］。文中沒有特殊說明的時間序列均指坐標殘差時間序列。考慮到地震引起的形變，地震前后的站點坐標參考值分別通過地震發(fā)生前后一段時間的靜態(tài)GPS觀測得到。當?shù)卣鸢l(fā)生時，地震波會導(dǎo)致站點的坐標時間序列出現(xiàn)較大的擺動波動，若直接求定地震前后的時間序列和地震當天序列的相似性，必然會造成相似性度量的不準確、不客觀，同時也無法分辨強震導(dǎo)致震后相似性發(fā)生變化的情況。為解決該問題，本文提出采用分段相似性度量確定的方法，具體計算過程如下。

假定地震當天、前一天及后一天的坐標時間序分別為Xe（ti）、Xb（ti）和Xp（ti），（i＝1，2，…，n），地震波到達的時刻為tk，持續(xù)歷元為s，即tk～tk＋s為受地震波影響的時間區(qū)間。根據(jù)地震波持續(xù)時間，將Xe（ti）、Xb（ti）和Xp（ti）分別劃分為兩段（根據(jù)實際情況可分為若干段）

一般情況下，地震前后兩段時間序列的相似性基本相同，同時由于單歷元定位誤差的周期性，地震當天的時間序列與前后兩天序列的相似性也基本相同，且會高于某一閾值ε，即

①ε前的符號是對皮爾遜相關(guān)系數(shù)或UIQI等函數(shù)而言，當S指的是歐氏距離時，符號需反號，即“≥”將變換為“≤”

式中

3 顧及分段相似性的恒星日濾波優(yōu)化算法

3.1 基本思路

顧及分段相似性的恒星日濾波優(yōu)化算法的基本思路，是在利用地震前后兩天的時間序列構(gòu)建濾波算子前，進行分段相似性度量的求定和判斷，根據(jù)相似性度量S的大小進行定權(quán)，然后進行加權(quán)平均構(gòu)建濾波時間序列。給定相似性度量的閾值ε，當時間序列的相似性度量低于ε時即認為不相關(guān)，不參與恒星日濾波的構(gòu)建。具體相似性度量的情況及其處理方式總結(jié)見表1所示。第4類情況即為發(fā)生強震的情況，由于無法分辨出環(huán)境發(fā)生變化的時刻，因而無法進行恒星日濾波。

表 1分段相似性度量及濾波處理方案Tab.1 Different segments similarity measure and its filtering solution

根據(jù)大量的試驗統(tǒng)計和相關(guān)文獻的研究結(jié)果［21］，本文中對不同的相似性度量選取的閾值均為ε＝0.5。由于不同的相似性度量的含義不同，對歐氏距離而言通過的是S≤0.5，而對皮爾遜相關(guān)系數(shù)和UIQI而言通過的則是｜S｜≥0.5。

3.2 計算步驟

根據(jù)上述基本思路，顧及分段相似性的恒星日濾波優(yōu)化算法的具體計算步驟見圖1。

4 算例分析

為評估相似性度量指標對恒星日濾波的影響，采用本文提出的方法對江蘇省CORS網(wǎng)某GPS基準站的坐標殘差時間序列進行濾波分析。試驗數(shù)據(jù)對應(yīng)的GPS年積日為2011年069、070和071 3d，亦即2011年3月11日日本地震前后3d的數(shù)據(jù)，069和071的時間序列已經(jīng)根據(jù)計算得到的GPS衛(wèi)星重復(fù)周期進行了相應(yīng)平移。對站點071天N、E、U方向坐標殘差時間序列分別加入了－4.0～4.0cm區(qū)間的隨機噪聲，模擬GPS解算衛(wèi)星、解算策略及解算參數(shù)個數(shù)不同等原因?qū)е碌亩ㄎ徽`差或粗差。為了增強模擬算例的可靠性，模擬了兩組不同的數(shù)據(jù)。站點3天各方向的坐標殘差原始時間序列及加入噪聲的時間序列見圖2所示。圖中橫軸坐標是相對于2011年3月11日UTC05∶30∶05的歷元數(shù)。3個子圖分別代表了N、E、U方向的時間序列，每個圖中從上至下分別為年積日069、070、071、模擬數(shù)據(jù)一和模擬數(shù)據(jù)二的坐標殘差時間序列，同時給出了坐標殘差的中誤差統(tǒng)計值。為了合理評定解算的中誤差，在統(tǒng)計前先將發(fā)生地震時1450～1850歷元的數(shù)據(jù)剔除。從圖中可以明顯看出，069、070和071 3d的坐標殘差時間序列的變化趨勢非常相近，這也驗證了前述周期性誤差的存在，071天加入噪聲后相似度明顯降低。

圖1 顧及分段相似性的恒星日濾波優(yōu)化算法流程Fig.1 Process steps of improved sidereal filtering accounting for segments'similarity

對這些坐標殘差時間序列利用常規(guī)恒星日濾波以及分別采用歐氏距離、相關(guān)系數(shù)和UIQI等相似性度量的恒星日濾波優(yōu)化算法進行處理。表2為各種濾波處理后的坐標殘差時間序列的統(tǒng)計。從表2可以看出：

（1）常規(guī)的恒星日濾波方法應(yīng)對特殊情況的能力最差，濾波經(jīng)常沒有任何作用，如模擬算例1的N方向濾波前后精度保持不變；甚至有時會出現(xiàn)原始序列濾波后精度變差的情況，如模擬算例2的U方向精度從原始序列的±2.23cm降低為濾波后的±2.56cm。

圖2 站點3d的坐標殘差時間序列及第3天加入噪聲的模擬序列Fig.2 Coordinate residual time series of 3days and the 3th day's series with the simulation of the noise

（2）顧及分段相似性的優(yōu)化算法則大大改善了恒星日濾波的效果，提高了濾波的可靠性，特別是采用皮爾遜相關(guān)系數(shù)和UIQI作為相似性度量的優(yōu)化算法，對于兩個模擬算例都能得到滿意的結(jié)果，且二者的精度相當。

（3）選擇歐氏距離作為相似性度量，雖然有時能夠得到比較滿意的結(jié)果，但有時精度較差，表現(xiàn)不穩(wěn)定。

（4）由于UIQI方法不僅顧及了時間序列間的相關(guān)系數(shù)，而且考慮了序列自身的均值和方差，因而本文建議采用UIQI作為恒星日濾波優(yōu)化算法的相似性度量。

5 結(jié) 論

本文將數(shù)據(jù)挖掘理論中的時間序列相似性分析引入恒星日濾波，提出了基于分段相似性度量的同震形變恒星日濾波優(yōu)化算法，并研究了不同相似性度量指標對恒星日濾波優(yōu)化算法的影響。研究結(jié)果表明，當由于GPS解算相關(guān)參數(shù)不同等原因引起解算誤差較大或存在粗差時，或者強震后地表建筑物等發(fā)生了巨大變化而引起多路徑效應(yīng)等周期性誤差發(fā)生根本性改變時，常規(guī)恒星日濾波非但不能起到提高精度的作用，有時還會適得其反。利用本文提出的顧及分段相似性度量的恒星日濾波優(yōu)化算法則能有效克服傳統(tǒng)恒星日濾波的缺點，大大提高了恒星日濾波的精度和可靠性。同時也發(fā)現(xiàn)，選擇歐氏距離作為相似性度量的恒星日濾波算法精度不穩(wěn)定；采用皮爾遜相關(guān)系數(shù)和UIQI作為相似性度量的優(yōu)化算法，二者的精度較高且精度相當；由于UIQI方法不但顧及了時間序列間的相關(guān)系數(shù)，而且還考慮了序列自身的均值和方差，因而建議采用UIQI作為恒星日濾波優(yōu)化算法的相似性度量。

表2 恒星日濾波優(yōu)化算法的精度統(tǒng)計Tab.2 Accuracy statistics of improved sidereal filtering cm

［1］ BILICH A，CASSIDY J F，LARSON K M.GPS Seismology：Application to the 2002Mw7.9Denali Fault Earthquake［J］.Bulletin of the Seismological Society of America，2008，98（2）：593－606.

［2］ MENG Guojie，REN Jinwei，JIN Honglin，et al.Data Processing Methods of High Rate GPS and Its Application to Seismology［J］.Recent Developments in World Seismology，2007（7）：613－620.（孟國杰，任金衛(wèi)，金紅林，等.GPS高頻數(shù)據(jù)處理方法及其在地震學(xué)中的應(yīng)用研究進展［J］.國際地震動態(tài)，2007（7）：26－30.）

［3］ XU Caijun，LIU Yang，WEN Yangmao.Mw7.9Wenchuan Earthquake Slip Distribution Inversion from GPS Measurements［J］.Acta Geodaetica et Cartographica Sinica，2009，38（3）：195－201.（許才軍，劉洋，溫揚茂.利用GPS資料反演汶川Mw7.9級地震滑動分布［J］.測繪學(xué)報，2009，38（3）：195－201.）

［4］ ZHONG Ping，DING Xiaoli，ZHENG Dawei，et al.Filter－based GPS Structural Vibration Monitoring Methods and Comparison of Their Performances［J］.Acta Geodaetica et Cartographica Sinica，2007，36（1）：31－36.（鐘萍，丁曉利，鄭大偉，等.GPS結(jié)構(gòu)振動監(jiān)測數(shù)據(jù)濾波方法及其性能實驗研究［J］.測繪學(xué)報，2007，36（1）：31－36.）

［5］ BOCK Y.Continuous Monitoring of Crustal Deformation［J］.GPS World，1991，2（6）：40－47.

［6］ GENRICH J F，BOCK Y.Rapid Resolution of Crustal Motion at Short Ranges with the Global Positioning System［J］.Journal of Geophysical Research，1992，97 （B3）：3261－3269.

［7］ HUANG Dingfa，DING Xiaoli，CHEN Yongqi，et al.Wavelet Filters Based Separation of GPS Multipath Effects and Engineering Structural Vibrations［J］.Acta Geodaetica et Cartographica Sinica，2001，30（1）：36－41.（黃丁發(fā)，丁曉利，陳永奇，等.GPS多路徑效應(yīng)影響與結(jié)構(gòu)振動的小波濾波篩分研究［J］.測繪學(xué)報，2001，30（1）：36－41.）

［8］ XIONG Yongliang，DING Xiaoli，HUANG Dingfa，et al.Integrated Single Epoch Algorithm Based on Wavelet Transform and Its Application to Structural Vibration Monitoring［J］.Acta Geodaetica et Cartographica Sinica，2005，34（3）：202－207.（熊永良，丁曉利，黃丁發(fā)，等.基于小波變換的單歷元算法及其在結(jié)構(gòu)振動監(jiān)測中的應(yīng)用研究［J］.測繪學(xué)報，2005，34（3）：202－207.

［9］ SEEBER G，MENGE F，VOLKSEN C，et al.Precise GPS Positioning Improvements by Reducing Antenna and Site Dependent Effects［C］∥Advances in Positioning and Reference Frames.Berlin：Springer，1997：237－244.

［10］ PARK K D，NEREM R S，SCHENEWERK M S，et al.Site－specific Multipath Characteristics of Global IGS and CORS GPS Sites［J］.Journal of Geodesy，2004，77（12）：799－803.

［11］ CHOI K，BILICH A，LARSON M，et al.Modified Sidereal Filtering：Implications for High－rate GPS Positioning［J］.Geophysical Research Letters，2004，31（22）：1－4.

［12］ AGNEW D C，LARSON K M.Finding the Repeat Times of the GPS Constellation［J］.GPS Solutions，2006，11（1）：71－76.

［13］ LARSON K M，BILICH A，AXELRAD P.Improving the Precision of High－rate GPS［J］.Journal of Geophysical Research：Solid Earth，2007，112（B05）：1－11.

［14］ YIN Haitao，GAN Weijun，XIAO Genru.Modified Sidereal Filter and Its Effect on High－rate GPS Positioning［J］.Geomatics and Informaton Science of Wuhan University，2011，36（5）：609－611.（殷海濤，甘衛(wèi)軍，肖根如.恒星日濾波的修正以及對高頻GPS定位的影響研究［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2011，36（5）：609－611.）

［15］ FANG Jiaguo.Time Series Data Mining Based on Similarity Analysis［D］.Hangzhou：Zhejiang University，2011.（方加果.基于相似性分析的時間序列數(shù)據(jù)挖掘算法研究［D］.杭州：浙江大學(xué)，2011.）

［16］ YANG Kang，LI Manchun，LIU Yongxue，et al.Accumulated Similarity Surface for Spatial Weights Matrix Construction［J］.Acta Geodaetica et Cartographica Sinica，2012，41（2）：259－265.（楊康，李滿春，劉永學(xué)，等.基于累積相似度表面的空間權(quán)重矩陣構(gòu)建方法［J］.測繪學(xué)報，2012，41（2）：259－265.）

［17］ ZHANG Jun.Research on Data Mining Technology of Pattern－based Similarity Search in Time Series Database［D］.Nanjing：Southeast University，2006.（張軍.基于時間序列相似性的數(shù)據(jù)挖掘方法研究［D］.南京：東南大學(xué)，2006.）

［18］ WANG Z，BOVIK A C.A Universal Image Quality Index［J］.IEEE Signal Processing Letters，2002，9（3）：81－84.

［19］ FANG Rongxin，SHI Chuang，GU Shengfeng.Precise Point Positioning with High－rate GPS Data Applied to Seismic Displacements Analysis［J］.Geomatics and Informaton Science of Wuhan University，2009，34（11）：1340－1343.（方榮新，施闖，辜聲峰.基于PPP動態(tài)定位技術(shù)的同震地表形變分析［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2009，34（11）：1340－1343.）

［20］ FANG Rongxin.High－rate GPS Data Non－difference Precise Processing and Its Application on Seismology ［D］.Wuhan：Wuhan University，2010.（方榮新.高采樣率GPS數(shù)據(jù)非差精密處理方法及其在地震學(xué)中的應(yīng)用研究［D］.武漢：武漢大學(xué)，2010.）

［21］ ZHANG Qian，WANG Zhihe，ZHANG Guozhi.An Algorithm for Discovering Positive and Negative Association Rules Based on Correlation Coefficient［J］.Journal of Shanxi University of Technology，2005，21（4）：35－38.（張倩，王治和，張國治.基于相關(guān)系數(shù)的正、負關(guān)聯(lián)規(guī)則挖掘算法［J］.陜西理工學(xué)院學(xué)報，2005，21（4）：35－38.）