董 偉

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

1 概述

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

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

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

2 基本理論

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

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

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

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

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

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

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

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

3 實(shí)例分析

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

圖1 GPS時(shí)間序列小波分析流程圖

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

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

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

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

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

圖2 KKN4站原始時(shí)間序列和去趨勢(shì)項(xiàng)時(shí)間序列

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

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

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

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

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

硬閾值處理方法:

軟閾值處理方法:

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

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

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

式中，SNR為信噪比，Ps為原始信號(hào)功率，Pz為噪聲功率。

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

圖5 GPS時(shí)間序列與硬閾值法、軟閾值法去噪對(duì)比

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

4 結(jié)語(yǔ)

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

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

［1］李婧.“陸態(tài)網(wǎng)絡(luò)”基準(zhǔn)站坐標(biāo)時(shí)間序列變化特性分析［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］陳強(qiáng)，黃聲享，王韋.小波去噪效果評(píng)價(jià)的另一指標(biāo)［J］.測(cè)繪信息與工程，2008，(5):13－14.

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

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