小波和傅里葉變換在坐標(biāo)時(shí)間序列分析中的應(yīng)用

蔡富，孫付平，戴海亮，朱新慧，張龍龍

(信息工程大學(xué) 地理空間信息學(xué)院, 河南 鄭州 450001)

0 引 言

由于各種地球物理效應(yīng),大部分GPS坐標(biāo)時(shí)間序列存在明顯的非線性變化,尤以垂直方向最為明顯[1]．如何有效地提取非線性變化中的有用信息是在進(jìn)行坐標(biāo)時(shí)間序列分析時(shí)經(jīng)常遇到的問題．傳統(tǒng)的諧波分析方法是建立在傅里葉變換的基礎(chǔ)上的,傅里葉變換具有較好的頻域分辨率,能夠精確地確定時(shí)間序列中各次諧波的頻率和振幅．但是傅里葉變換反映的是信號(hào)總體在整個(gè)時(shí)間段的頻域特征,當(dāng)將時(shí)間序列從時(shí)域變換到頻域時(shí),其在時(shí)間域上的信息就丟失了,因此不能很好地反映時(shí)間序列的細(xì)節(jié)特征．而且基于傅里葉變換的時(shí)間序列分析方法一般是建立在信號(hào)平穩(wěn)的假設(shè)基礎(chǔ)之上,而絕大部分的時(shí)間序列其整體趨勢(shì)是不平穩(wěn)的[2]．故而傅里葉變換具有一定的局限性．而小波變換則可以改變這種局限性,研究表明,小波變換可以對(duì)許多含有周期性規(guī)律的時(shí)間序列進(jìn)行時(shí)頻分析,在地球物理等領(lǐng)域應(yīng)用廣泛[3-5]．利用多尺度分析的原理,小波變換可以將時(shí)間序列分解到不同尺度上進(jìn)行研究,并能進(jìn)行不失真重構(gòu),能夠得到時(shí)間序列的細(xì)節(jié)信息,因此對(duì)處理非平穩(wěn)時(shí)間序列優(yōu)勢(shì)明顯[6-7]．但是通過小波變換得到的只是各次諧波的時(shí)域波形,還不能很精確地得到時(shí)間序列的頻率和振幅信息．而且實(shí)驗(yàn)表明,當(dāng)所處理的時(shí)間序列中含有較長(zhǎng)周期項(xiàng)時(shí),僅基于小波變換這一單一技術(shù)很難對(duì)其進(jìn)行準(zhǔn)確確定．目前，綜合傅里葉變換和小波分析的研究現(xiàn)狀來看,單一的研究途徑和方法不再適合于復(fù)雜的時(shí)間序列分析與預(yù)測(cè),而多種理論和方法的有機(jī)結(jié)合能為正確分析和解決問題提供有效途徑．

本文綜合考慮小波變換時(shí)頻局部分析優(yōu)勢(shì)和傅里葉變換頻域分析優(yōu)勢(shì),提出一種小波變換與傅里葉變換相結(jié)合的方法來對(duì)坐標(biāo)時(shí)間序列進(jìn)行分析．

1 基本理論

1.1 傅里葉變換

傅里葉變換能夠?qū)⑿盘?hào)在時(shí)域和頻域中相互轉(zhuǎn)化, 其實(shí)質(zhì)是把時(shí)域中的原始信號(hào)f(t)分離為多種不同周期(頻率)的正、余弦波之和,進(jìn)而對(duì)時(shí)間序列進(jìn)行分析．這種變換是可逆的且能量保持不變[8]．定義如下式:

(1)

其反變換為

(2)

式中：f(t)是時(shí)間域的信號(hào)；F(jω)是快速傅里葉變換的結(jié)果,為函數(shù)f(t)的頻譜密度函數(shù),f(t)為F(jω)的原函數(shù)[9]．根據(jù)定義可知,原始信號(hào)f(t)在經(jīng)過傅里葉變換后其時(shí)間域的信息就消失了,即F(jω)只有頻率特性,且其值由f(t)在整個(gè)時(shí)段的特性所決定,利用傅里葉變換這一特性可獲取信號(hào)所有頻率．

通過傅里葉變換,在時(shí)域中連續(xù)變化的信號(hào)可轉(zhuǎn)化為頻域中的信號(hào),因此傅里葉變換反映的是整個(gè)信號(hào)在全部時(shí)間下的整體頻域特征,但不能反映信號(hào)的局部特征．

1.2 小波變換

小波變換是一種時(shí)間-頻率分析方法,與傅里葉變換相比,在分析非平穩(wěn)時(shí)間序列方面優(yōu)勢(shì)明顯,因此成為當(dāng)前分析時(shí)間序列的有效工具．小波變換的定義如下[10]:

設(shè)f(t)是一平方可積函數(shù),記作f(t)∈L2(R),ψ(t)為母小波,如果ψ(t)滿足容許條件:

(3)

則有下式:

(4)

式(4)為函數(shù)f(t)的小波變換,式中：a為尺度因子；b為位移因子．

另外,小波變換也可以寫成式(5)所示的卷積形式:

(5)

(6)

1.3 小波變換的多分辨率分析

1988 年Mallat提出的多分辨率分析(MRA)是小波分析的重要概念,又稱為多尺度分析．小波變換的MRA能將時(shí)間序列在不同分辨率上進(jìn)行分解,分別得到低頻部分和高頻部分,其中低頻部分能夠反映時(shí)間序列的概貌,而高頻部分又叫細(xì)節(jié)部分,能夠刻畫時(shí)間序列的細(xì)節(jié)信息,因而將小波變換譽(yù)為“數(shù)學(xué)顯微鏡”[11-13]．小波變換的這種多尺度分解能力可以將各種頻率疊加在一起的信號(hào)分解為不同頻率的子信號(hào),使得一種分辨率下無法發(fā)現(xiàn)的特性在另一個(gè)分辨率下很容易被發(fā)現(xiàn),從而能夠有效地應(yīng)用于信號(hào)和噪聲分離、信號(hào)的分析與重構(gòu)、特征提取等問題[14-15]．

小波多分辨率分析主要是通過空間分解,將信號(hào)分解在子空間列中,方便進(jìn)行信號(hào)分析[16]．小波多分辨率分解示意圖如圖1所示,其中S表示原始時(shí)間序列,A1、A2、A3表示小波分解后各層時(shí)間序列的低頻部分,D1、D2、D3則表示小波分解后各層時(shí)間序列的高頻部分．小波變換能將信號(hào)不斷分解成低頻和高頻兩部分,進(jìn)而由粗到細(xì)地逐步觀察時(shí)間序列,具有多尺度分析的特點(diǎn)．此外,各部分與原始信號(hào)具有以下關(guān)系:S=A3+D3+D2+D1.以此類推,若要進(jìn)行進(jìn)一步的分解,則可以把低頻部分A3作為高頻部分D4和低頻部分A4的疊加．可以將小波分析視為近似無損的信號(hào)處理方法．

圖1 小波多分辨率分解示意圖

大多數(shù)情況下,由于受到不同地球物理因素的影響,測(cè)站時(shí)間序列表現(xiàn)出來的隨機(jī)性較強(qiáng),異常變化在部分時(shí)段內(nèi)時(shí)有發(fā)生．頻譜分析雖然能夠分析出時(shí)間序列中有色噪聲的頻率特征,并從頻域發(fā)現(xiàn)時(shí)間序列的周期性趨勢(shì)項(xiàng),但是卻不能有效探測(cè)異常區(qū)間．通過小波分析,就能夠根據(jù)小波多尺度分析的原理對(duì)分解后的各層時(shí)間序列進(jìn)行分析,確定異常區(qū)間,從而達(dá)到提取測(cè)站時(shí)間序列特征信息的目的[16]．

2 小波和傅里葉變換相結(jié)合的方法及分析

2.1 數(shù)據(jù)描述

實(shí)驗(yàn)數(shù)據(jù)來源于國(guó)際地球參考框架ITRF2008的GPS坐標(biāo)殘差時(shí)間序列數(shù)據(jù),可在官網(wǎng)(http://itrf.ensg.ign.fr/ITRF-solutions/2008/)上下載,而且殘差時(shí)間序列已經(jīng)剔除了固體潮、極潮等部分環(huán)境負(fù)荷的影響[17-18]．由于GPS 測(cè)站遍布全球,各站地理環(huán)境不同導(dǎo)致坐標(biāo)時(shí)間序列數(shù)據(jù)質(zhì)量不一,時(shí)間序列所表現(xiàn)的特征信息如周期規(guī)律等也就存在差異,有的測(cè)站還存在時(shí)間間隔較長(zhǎng)的間斷點(diǎn)和突變信息．為了使結(jié)論更具有普遍性,本文分別從低、中、高緯度各選取了3～4個(gè)測(cè)站進(jìn)行了實(shí)驗(yàn)分析, 所選的測(cè)站信息如表1所示．

表1 GPS測(cè)站位置

2.2 小波基函數(shù)的選取

選取不同的小波基函數(shù),小波變換的結(jié)果也不同．目前廣泛使用的小波基函數(shù)有haar小波、dbN小波、symN小波、coif小波、bior小波、rbio小波等,不同的小波基函數(shù)有不同的對(duì)稱性、消失矩和正則性等性能特點(diǎn),在波形上也表現(xiàn)出了很大的差別[19]．本文在對(duì)坐標(biāo)時(shí)間序列數(shù)據(jù)進(jìn)行大量的實(shí)驗(yàn)分析表明選用dbN小波效果較好．dbN小波作為稀疏基所引入的光滑誤差很難被探測(cè),因此信號(hào)重構(gòu)過程比較光滑．伴隨著階次的增大,消失矩階數(shù)亦隨之增大,光滑性變得越好,頻域的局部化能力就越強(qiáng),劃分頻帶效果越好,但時(shí)域緊支撐性減弱,同時(shí)計(jì)算量劇增,實(shí)效性變差[20]．另外,dbN小波沒有明確的表達(dá)式,除N=1外(N=1時(shí)即為Haar小波)．

2.3 隨機(jī)項(xiàng)的提取

去除長(zhǎng)期趨勢(shì)項(xiàng)后,坐標(biāo)時(shí)間序列中的其它層就是隨機(jī)項(xiàng)和周期項(xiàng),而且在地球物理現(xiàn)象和經(jīng)濟(jì)領(lǐng)域等實(shí)際問題中,隨機(jī)因素往往是小尺度成分,因此,在大多數(shù)情況下,第一次小波分解的高頻部分中包含的隨機(jī)項(xiàng)的成分最多,其它層中包含的周期項(xiàng)成分較多[21]．如圖2所示的URUM站垂向坐標(biāo)時(shí)間序列中,既有周期性規(guī)律，還存在著一些隨機(jī)變動(dòng),使用具有正交性和高度緊支撐性的db4小波,分解5層分別得到a1～a5的低頻部分和d1～d5的高頻部分,隨機(jī)變動(dòng)信息通常包含在每一層的高頻部分中．如圖3所示,可以發(fā)現(xiàn),高頻部分d1、d2中存在著一些隨機(jī)變動(dòng)信息,即隨機(jī)項(xiàng),而且由d1可以看出在URUM站垂向坐標(biāo)時(shí)間序列在第475周附近存在跳變,量值達(dá)到8 mm．表明了小波變換對(duì)細(xì)節(jié)信息的把控能力比較強(qiáng)．

在對(duì)坐標(biāo)時(shí)間序列進(jìn)行小波多尺度分解后,對(duì)于包含隨機(jī)項(xiàng)的各層,其在不同尺度上的成分一般相互獨(dú)立,而且經(jīng)常是平穩(wěn)的,因此可以分別進(jìn)行處理．

圖2 URUM站垂向坐標(biāo)時(shí)間序列

圖3 小波分解后的各層低頻和高頻部分

2.4 周期項(xiàng)的提取

小波分析是一種時(shí)、頻MRA方法,它可以對(duì)函數(shù)和信號(hào)系列進(jìn)行多尺度細(xì)化分析,以分析不同尺度(周期)隨時(shí)間的演變情況．小波分析能將時(shí)間序列的頻率特征在時(shí)域上顯現(xiàn)出來,分析其主周期項(xiàng),而且能準(zhǔn)確地給出變化趨勢(shì)和突變點(diǎn)以及各時(shí)間尺度(周期)的長(zhǎng)短和分布情況．將待分析時(shí)間序列通過合適的小波基函數(shù)和分解層次分解成低頻部分和高頻部分,并投射到相應(yīng)尺度上．對(duì)剩余細(xì)節(jié)信號(hào)采用傅里葉變換進(jìn)行周期分析,實(shí)現(xiàn)周期性信號(hào)的探測(cè)．如果有些層的周期性還不十分明朗,可對(duì)該層再進(jìn)行小波分解,直到能夠方便地分析該層的波形為止．

如圖3所示,由低頻部分a3和a4可以看出URUM站垂向坐標(biāo)時(shí)間序列中存在1年周期項(xiàng),并且為其主周期項(xiàng),由d4可以看出半年周期項(xiàng),由a5第0～100周時(shí)間段可以看出時(shí)間序列中存在2年周期項(xiàng)．這與快速傅里葉變換的結(jié)果也是相吻合的．雖然小波變換存在諸多優(yōu)勢(shì),但分析發(fā)現(xiàn)很難確定主要周期項(xiàng)的振幅,而且對(duì)于一些測(cè)站存在長(zhǎng)周期如10年、12年,以及短周期項(xiàng)如季節(jié)周期項(xiàng)、月周期項(xiàng)的情況,小波變換優(yōu)勢(shì)并不明顯．因此,在小波變換的基礎(chǔ)上,對(duì)分解得到的高頻或低頻信息進(jìn)行傅里葉變換,提取其周期項(xiàng)及相應(yīng)的振幅信息．如圖4在對(duì)BJFS站N方向上坐標(biāo)時(shí)間序列進(jìn)行小波分解時(shí),雖然能根據(jù)高頻部分d5看出周期為1年的周期項(xiàng),但是其他周期項(xiàng)則不明顯．鑒于快速傅里葉變換頻率特性較好,這里對(duì)低頻部分a5進(jìn)行傅里葉變換,可以得到周期為10年,振幅為1.44 mm的周期項(xiàng),以及周期為3年,振幅為1.39 mm的周期,如圖5所示．對(duì)d5進(jìn)行快速傅里葉變換可以確定其1年周期項(xiàng)的振幅為0.80 mm．以此類推,將BJFS站N方向小波分解得到的各層主要周期項(xiàng)及振幅進(jìn)行統(tǒng)計(jì)結(jié)果如表2所示,“周期/a”中a指周年．

表2 小波分解得到的各層主要周期及振幅

圖4 BJFS站N方向小波分解后的各層低頻和高頻部分

圖5 小波第5次分解得到的低頻和高頻分量頻譜圖

測(cè)站名稱低頻部分周期/a振幅/mm周期/a高頻部分振幅/mm周期/a振幅/mm TUVA-N10.90.50.50.20.3 TUVA-E101.50.50.60.30.4 TUVA-H12.40.520.251 BOGT-N11.20.50.40.20.23 BOGT-E10.90.50.60.250.3 BOGT-H14.40.52.30.21.2 LHAZ-N130.510.30.2 LHAZ-E101.20.50.480.250.3 LHAZ-H160.51.70.30.5 ARTU-N10.80.50.30.250.2 ARTU-E11.70.50.30.250.13 ARTU-H14.40.510.250.5 MAC1-N100.90.50.60.250.3 MAC1-E210.50.40.250.3 MAC1-H12.40.50.60.250.8 BJFS-N101.40.50.40.250.2 BJFS-E11.50.50.30.250.2 BJFS-H14.60.52.00.250.6 URUM-N21.40.50.430.250.5 URUM-E510.50.40.250.3 URUM-H14.20.510.20.8 MTJO-N110.50.90.20.2 MTJO-E11.50.50.20.250.3 MTJO-H130.510.250.7 STR1-N10.60.50.50.240.3 STR1-E10.80.50.60.230.3 STR1-H14.20.451.30.30.5 BAN2-N11.50.50.430.230.5 BAN2-E1.51.50.50.80.30.2 BAN2-H150.530.250.8 KELY-N10.40.50.250.250.3 KELY-E10.50.50.360.240.3 KELY-H12.40.510.250.5

采用上述數(shù)據(jù)處理流程和方法,本文對(duì)所選測(cè)站進(jìn)行了處理分析,鑒于篇幅,只對(duì)小波分解后的各層主要周期項(xiàng)及振幅變化進(jìn)行了統(tǒng)計(jì),結(jié)果如表3所示．

綜合上述實(shí)驗(yàn)分析,高低頻小波MRA對(duì)坐標(biāo)時(shí)間序列周期項(xiàng)的提取各有優(yōu)缺點(diǎn)．低頻分析可以較容易地得到“周年項(xiàng)”和 “兩年周期項(xiàng)”,但較長(zhǎng)周期項(xiàng)(2年以上周期項(xiàng))和短周期項(xiàng)(1年以下周期項(xiàng))則不易發(fā)現(xiàn),可以結(jié)合傅里葉變換進(jìn)一步求取;高頻分析能夠較準(zhǔn)確提取“半周年項(xiàng)”、“一季項(xiàng)”等短周期．基于小波變換和傅里葉變換相結(jié)合的綜合分析方法，可以對(duì)小波分解的每層信號(hào)進(jìn)行分析,能夠結(jié)合兩者的優(yōu)勢(shì),達(dá)到提取坐標(biāo)時(shí)間序列中特征信息的目的．

3 結(jié) 論

本文依據(jù)GPS測(cè)站殘差時(shí)間序列非線性變化的特點(diǎn)以及小波變換和傅里葉變換在坐標(biāo)時(shí)間序列分析中各自的優(yōu)勢(shì),提出了一種小波分析和傅里葉變換相結(jié)合的時(shí)頻分析方法,來提取坐標(biāo)時(shí)間序列中的特征信息,得到了以下結(jié)論:

1)利用傅里葉變換可以得到時(shí)間序列中各次諧波的頻率和振幅,但當(dāng)將時(shí)間序列從時(shí)域變換到頻域時(shí),容易丟失其在時(shí)間域上的信息,不能很好地反映時(shí)間序列的一些細(xì)節(jié)特征,如間斷點(diǎn)和突變信息等,因此具有一定的局限性．

2)小波MRA可以將坐標(biāo)時(shí)間序列分解到不同尺度上,得到時(shí)間序列的細(xì)節(jié)特征和一些周期信息,在處理非平穩(wěn)時(shí)間序列時(shí)優(yōu)勢(shì)明顯．其中低頻分析可以直觀地得到“周年項(xiàng)”和 “兩年周期項(xiàng)”,但2年以上的長(zhǎng)周期項(xiàng)和1年以下的短周期項(xiàng)則不容易發(fā)現(xiàn);而高頻分析能夠較準(zhǔn)確提取“半周年項(xiàng)”、“一季項(xiàng)”等短周期．但是通過小波變換得到的只是各次諧波的時(shí)域波形,還不能很精確地得到時(shí)間序列的頻率和振幅信息．

3)與單一采用小波變換或傅里葉變換相比,小波和傅里葉變換相結(jié)合的方法既能保留測(cè)站坐標(biāo)時(shí)間序列的時(shí)間信息,又能在不同尺度上對(duì)時(shí)間序列進(jìn)行分析,解決了小波變換部分周期項(xiàng)難提取和傅里葉變換時(shí)間信息容易丟失的問題,非常適用于對(duì)復(fù)雜的時(shí)間序列進(jìn)行分析和研究．而且,這種方法在一定程度上也克服了小波基函數(shù)的選取問題,因此具有較高的研究?jī)r(jià)值.