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

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

1 引 言

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

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

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

2 同震坐標(biāo)時(shí)間序列的相似性分析

2.1 時(shí)間序列相似性定義

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

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

2.2 時(shí)間序列相似性度量指標(biāo)

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

2.2.1 歐氏距離（Euclidean distance，ED）

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

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

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

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

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

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

此式可分解為

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

2.3 時(shí)間序列分段相似性的確定

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

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

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

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

式中

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

3.1 基本思路

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

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

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

3.2 計(jì)算步驟

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

4 算例分析

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

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

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

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

圖2 站點(diǎn)3d的坐標(biāo)殘差時(shí)間序列及第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)化算法，對(duì)于兩個(gè)模擬算例都能得到滿意的結(jié)果，且二者的精度相當(dāng)。

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

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

5 結(jié) 論

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

表2 恒星日濾波優(yōu)化算法的精度統(tǒng)計(jì)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.（孟國(guó)杰，任金衛(wèi)，金紅林，等.GPS高頻數(shù)據(jù)處理方法及其在地震學(xué)中的應(yīng)用研究進(jìn)展［J］.國(guó)際地震動(dòng)態(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.（許才軍，劉洋，溫?fù)P茂.利用GPS資料反演汶川Mw7.9級(jí)地震滑動(dòng)分布［J］.測(cè)繪學(xué)報(bào)，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)振動(dòng)監(jiān)測(cè)數(shù)據(jù)濾波方法及其性能實(shí)驗(yàn)研究［J］.測(cè)繪學(xué)報(bào)，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)振動(dòng)的小波濾波篩分研究［J］.測(cè)繪學(xué)報(bào)，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)振動(dòng)監(jiān)測(cè)中的應(yīng)用研究［J］.測(cè)繪學(xué)報(bào)，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)軍，肖根如.恒星日濾波的修正以及對(duì)高頻GPS定位的影響研究［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2011，36（5）：609－611.）

［15］ FANG Jiaguo.Time Series Data Mining Based on Similarity Analysis［D］.Hangzhou：Zhejiang University，2011.（方加果.基于相似性分析的時(shí)間序列數(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］.測(cè)繪學(xué)報(bào)，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í)間序列相似性的數(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動(dòng)態(tài)定位技術(shù)的同震地表形變分析［J］.武漢大學(xué)學(xué)報(bào)：信息科學(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ó)治.基于相關(guān)系數(shù)的正、負(fù)關(guān)聯(lián)規(guī)則挖掘算法［J］.陜西理工學(xué)院學(xué)報(bào)，2005，21（4）：35－38.）