董 偉

(中鐵第四勘察設(shè)計院集團有限公司，湖北 武漢 430063)

1 概述

GPS已被廣泛運用于地殼形變的研究中。將GPS資料制成時間序列，通過分析GPS的連續(xù)觀測資料便可以得到測站點與時間及空間的關(guān)系，并推估出各種構(gòu)造，從而獲得與地殼形變之間的關(guān)聯(lián)性［1］。

目前常用的GPS時間序列分析方法對不平穩(wěn)的時間序列的分析無能為力，而小波分析能有效地分析地學(xué)現(xiàn)象中不平穩(wěn)的時間序列。根據(jù)小波的多尺度特點，可以將信號進行分解和重構(gòu)，對每一層進行處理，從而獲得所需要的部分。

本文在小波分析的基礎(chǔ)上對GPS時間序列進行了去除觀測噪聲的研究。

2 基本理論

小波分析的基礎(chǔ)是進行小波變換，而小波變換就是通過小波函數(shù)系去表示或逼近一個信號，可由基小波函數(shù)通過平移和伸縮構(gòu)成。

基小波函數(shù)可設(shè)為ψ(t)∈L2(R)，其中L2(R)表示平方可積的實數(shù)空間，小波變換基底則可采用下式進行定義:

式中，a為尺度因子，b為平移因子。

對于給定的能量有限信號或函數(shù) f(t)∈L2(R)，其連續(xù)小波變換可得到該信號或者函數(shù)與小波函數(shù)的內(nèi)積:

式中，Wf(a，b)為小波變換系數(shù)的復(fù)共軛函數(shù)。而對于離散的數(shù)據(jù)，只需要將上式的積分形式改成求和形式即可。該變換過程在實際運用中主要是為了獲得小波系數(shù)，從而通過利用這些系數(shù)來分析時間序列的時頻變化特征。

本文在進行數(shù)據(jù)處理時所采用的是morlet小波函數(shù)，在實際計算過程中可以通過選擇一套離散化的尺度變量來計算小波變換，選取尺度參數(shù)原則上可以通過下式確定［2］:

式中，a0為可分解的最小的尺度，δj為離散尺度之間的間距，J為決定最大尺度的因子，N為數(shù)據(jù)的個數(shù)，其中a0一般選擇兩倍的采樣間隔，即a0=2δt，選取較小的δj可以得到更好的尺度分辨率，但計算和畫圖會變得非常慢，對于morlet小波，本文將選擇δj=0.125。

由于GPS時間序列是一維數(shù)據(jù)，根據(jù)Mallat算法，如果已知雙尺度方程中的濾波器系數(shù)，就可以快速計算出各尺度的逼近和細(xì)節(jié)。

3 實例分析

本文使用的GPS時間序列為“陸態(tài)網(wǎng)絡(luò)”中KKN4站時間序列，利用小波分析GPS時間序列的主要工作流程見下圖1。

圖1 GPS時間序列小波分析流程圖

由于坐標(biāo)時間序列的小波分析要求原始坐標(biāo)時序具有零均值的特性，因此需要對部分原始的GPS時間序列進行消除線性趨勢項和常數(shù)項的處理，利用消除線性趨勢項和常數(shù)項得到小波方差圖和小波功率譜圖，從而可以判斷出數(shù)據(jù)中的周期性質(zhì)。

在對KKN4站原始時間序列數(shù)據(jù)進行分析時，從圖2可以看出KKN4站原始時間序列中豎直方向有著較為明顯的線性趨勢，在去除其線性趨勢項后的時間序列變化量在零的上下波動，且分布較為均勻，已經(jīng)具備零均值特性，可用于后續(xù)小波分析。

在選取KKN4站2010～2013年的數(shù)據(jù)進行小波分析時首先獲得小波方差圖和小波功率譜，由于小波函數(shù)通常是一個復(fù)數(shù)，所以一個實數(shù)序列經(jīng)過小波變換以后可以得到的小波系數(shù)也是復(fù)數(shù)，再利用實部、虛部可以得到小波功率譜，小波方差和功率譜的具體獲得過程在此不作詳細(xì)介紹。

通過小波方差圖分析可以確定信號中不同的尺度擾動相對強度及其存在的主要時間尺度，圖3為得到的小波方差圖和功率譜圖。從小波方差圖可以明顯看出，大致在1/4，1/2、1和2處出現(xiàn)波峰，這與GPS時間序列中存在的季節(jié)性、半周年、周年和兩周年周期相吻合，這也說明GPS時間序列中確實存在著多個不同時間尺度的周期，同時從小波功率譜圖也可以明顯看出其周期性。

根據(jù)圖3的分析結(jié)果，確定將GPS時間序列進行4層小波分解，并確定閾值對每一層的高頻部分進行去噪處理，之后再對時間序列進行重構(gòu)。

圖2 KKN4站原始時間序列和去趨勢項時間序列

圖3 KKN4站GPS時間序列小波方差和功率譜

根據(jù)小波的多尺度特點，可以將信號進行分解和重構(gòu)，即用小波將信號進行分解，然后對每一層進行相關(guān)處理，最后重構(gòu)信號。在對GPS的時間序列數(shù)據(jù)進行分析時，需要考慮長周期運動、同震變形以及震后變形。GPS單站、單分量坐標(biāo)序列一般用下列周期性模型來表示:

式中，ti為時間，以年為單位;a為地殼位置;b為線性變換率;c、d、e、f為周年和半周年運動振幅;gj為地震造成的同震偏移;H(ti)為階躍函數(shù);vi為殘差值，代表觀測值與預(yù)測值間的差異，也就是包含GPS時間序列的噪聲信號。圖4為KKN4站2010～2013年間原始序列及進行4層分解后的高頻、低頻序列。

圖4 GPS時間序列的原始序列和4層分解的高頻、低頻序列

本文分步采用軟閾值和硬閾值兩種方法進行去噪處理，然后根據(jù)信噪比大小來確定去噪的好壞。

硬閾值處理方法:

軟閾值處理方法:

式中，|x|為小波變換的系數(shù)，T為預(yù)先選定的閾值，考慮到T值的選取會有不同，本文利用GPS數(shù)據(jù)解算的3倍中誤差來限定，設(shè)定T=3σ，σ為中誤差。

由于影響小波去噪處理效果的因素非常多，選擇不一樣小波基函數(shù)、不一樣的閾值以及不一樣的分解尺度，最終的去噪效果會有一定的差異，所以在進行降噪處理時必須以一些詳細(xì)的指標(biāo)來衡量去噪效果。

本文中所采用的評價方法是信噪比。信噪比是衡量信號中噪聲量度的傳統(tǒng)方法，其定義式［3］為:

式中，SNR為信噪比，Ps為原始信號功率，Pz為噪聲功率。

本文利用軟閾值和硬閾值兩種方法進行去噪處理及信號重構(gòu)后信號機信噪比對比圖見圖5［4－5］。

圖5 GPS時間序列與硬閾值法、軟閾值法去噪對比

從圖5可以非常清楚地看到軟閾值法去噪結(jié)果比硬閾值法去噪結(jié)果更加平滑，采用前者比采用后者進行處理的信號的信噪比小，但采用前者進行降噪處理的同時也掩蓋了某些由于閃爍噪聲和含隨機漫步噪聲而導(dǎo)致的點位突變，濾掉了部分變化較大的坐標(biāo)值。雖然采用軟閾值進行降噪處理便于分析基準(zhǔn)站點的運動趨勢，但缺少對基準(zhǔn)站點整體穩(wěn)定性分析，因此在具體運用中需要根據(jù)實際情況選擇合適的軟閾值或硬閾值降噪，才能獲得更好的結(jié)果。

4 結(jié)語

在利用小波進行GPS時間序列分析中，小波分析可以對信號進行不同尺度的分析，而且還可以將不同特效的噪聲非常有效地進行分離。另外小波分析還可以通過對GPS時間序列采用高低頻分析相結(jié)合的方式進行數(shù)據(jù)處理，能有效地提取季節(jié)性、半周年、周年及兩周年的周期項，再通過對信號進行降噪處理，從而可以更好地根據(jù)實際需要來分析GPS時間序列。

本文將GPS時間序列進行多尺度分析，盡量保持了邊緣位置精確，同時在去噪方面做到了折中處理。分析處理結(jié)果表明，小波分析可以很好地用于GPS時間序列的分析中，但由于處理方式不同，去噪處理最終的結(jié)果也會不同，因此現(xiàn)場實際使用過程中需要根據(jù)具體情況選取合適的處理方式。至于如何深入地進行多尺度分解、去噪及重構(gòu)，還需要進一步完善。

［1］李婧.“陸態(tài)網(wǎng)絡(luò)”基準(zhǔn)站坐標(biāo)時間序列變化特性分析［D］.鄭州:解放軍信息工程大學(xué)，2013.

［2］Torrence C，Compo G P.A practical guide to wavelet analysis.［J］.Bulletin of the American Meteorological Society，1998，79(1):61－78.

［3］陳強，黃聲享，王韋.小波去噪效果評價的另一指標(biāo)［J］.測繪信息與工程，2008，(5):13－14.

［4］周亞，王立峰，張思慧，等.IGS連續(xù)運行參考站高程時間序列功率譜分析［J］.太赫茲科學(xué)與電子信息學(xué)報，2014，(1):103 －107.

［5］范朋飛.高精度GPS站點坐標(biāo)時間序列分析與應(yīng)用［D］.西安:長安大學(xué)，2013.