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

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

0 引 言

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

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

1 基本理論

1.1 傅里葉變換

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

(1)

其反變換為

(2)

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

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

1.2 小波變換

小波變換是一種時(shí)間-頻率分析方法,與傅里葉變換相比,在分析非平穩(wěn)時(shí)間序列方面優(yōu)勢明顯,因此成為當(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]．小波變換的這種多尺度分解能力可以將各種頻率疊加在一起的信號分解為不同頻率的子信號,使得一種分辨率下無法發(fā)現(xiàn)的特性在另一個(gè)分辨率下很容易被發(fā)現(xiàn),從而能夠有效地應(yīng)用于信號和噪聲分離、信號的分析與重構(gòu)、特征提取等問題[14-15]．

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

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

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

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

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

實(shí)驗(yàn)數(shù)據(jù)來源于國際地球參考框架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 測站遍布全球,各站地理環(huán)境不同導(dǎo)致坐標(biāo)時(shí)間序列數(shù)據(jù)質(zhì)量不一,時(shí)間序列所表現(xiàn)的特征信息如周期規(guī)律等也就存在差異,有的測站還存在時(shí)間間隔較長的間斷點(diǎn)和突變信息．為了使結(jié)論更具有普遍性,本文分別從低、中、高緯度各選取了3～4個(gè)測站進(jìn)行了實(shí)驗(yàn)分析, 所選的測站信息如表1所示．

表1 GPS測站位置

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

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

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

去除長期趨勢項(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．表明了小波變換對細(xì)節(jié)信息的把控能力比較強(qiáng)．

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

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

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

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

小波分析是一種時(shí)、頻MRA方法,它可以對函數(shù)和信號系列進(jìn)行多尺度細(xì)化分析,以分析不同尺度(周期)隨時(shí)間的演變情況．小波分析能將時(shí)間序列的頻率特征在時(shí)域上顯現(xiàn)出來,分析其主周期項(xiàng),而且能準(zhǔn)確地給出變化趨勢和突變點(diǎn)以及各時(shí)間尺度(周期)的長短和分布情況．將待分析時(shí)間序列通過合適的小波基函數(shù)和分解層次分解成低頻部分和高頻部分,并投射到相應(yīng)尺度上．對剩余細(xì)節(jié)信號采用傅里葉變換進(jìn)行周期分析,實(shí)現(xiàn)周期性信號的探測．如果有些層的周期性還不十分明朗,可對該層再進(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)勢,但分析發(fā)現(xiàn)很難確定主要周期項(xiàng)的振幅,而且對于一些測站存在長周期如10年、12年,以及短周期項(xiàng)如季節(jié)周期項(xiàng)、月周期項(xiàng)的情況,小波變換優(yōu)勢并不明顯．因此,在小波變換的基礎(chǔ)上,對分解得到的高頻或低頻信息進(jìn)行傅里葉變換,提取其周期項(xiàng)及相應(yīng)的振幅信息．如圖4在對BJFS站N方向上坐標(biāo)時(shí)間序列進(jìn)行小波分解時(shí),雖然能根據(jù)高頻部分d5看出周期為1年的周期項(xiàng),但是其他周期項(xiàng)則不明顯．鑒于快速傅里葉變換頻率特性較好,這里對低頻部分a5進(jìn)行傅里葉變換,可以得到周期為10年,振幅為1.44 mm的周期項(xiàng),以及周期為3年,振幅為1.39 mm的周期,如圖5所示．對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次分解得到的低頻和高頻分量頻譜圖

測站名稱低頻部分周期/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ù)處理流程和方法,本文對所選測站進(jìn)行了處理分析,鑒于篇幅,只對小波分解后的各層主要周期項(xiàng)及振幅變化進(jìn)行了統(tǒng)計(jì),結(jié)果如表3所示．

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

3 結(jié) 論

本文依據(jù)GPS測站殘差時(shí)間序列非線性變化的特點(diǎn)以及小波變換和傅里葉變換在坐標(biāo)時(shí)間序列分析中各自的優(yōu)勢,提出了一種小波分析和傅里葉變換相結(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)勢明顯．其中低頻分析可以直觀地得到“周年項(xiàng)”和 “兩年周期項(xiàng)”,但2年以上的長周期項(xiàng)和1年以下的短周期項(xiàng)則不容易發(fā)現(xiàn);而高頻分析能夠較準(zhǔn)確提取“半周年項(xiàng)”、“一季項(xiàng)”等短周期．但是通過小波變換得到的只是各次諧波的時(shí)域波形,還不能很精確地得到時(shí)間序列的頻率和振幅信息．

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