小波和傅里葉變換在坐標時間序列分析中的應用

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

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

0 引 言

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

本文綜合考慮小波變換時頻局部分析優勢和傅里葉變換頻域分析優勢,提出一種小波變換與傅里葉變換相結合的方法來對坐標時間序列進行分析．

1 基本理論

1.1 傅里葉變換

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

(1)

其反變換為

(2)

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

通過傅里葉變換,在時域中連續變化的信號可轉化為頻域中的信號,因此傅里葉變換反映的是整個信號在全部時間下的整體頻域特征,但不能反映信號的局部特征．

1.2 小波變換

小波變換是一種時間-頻率分析方法,與傅里葉變換相比,在分析非平穩時間序列方面優勢明顯,因此成為當前分析時間序列的有效工具．小波變換的定義如下[10]:

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

(3)

則有下式:

(4)

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

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

(5)

(6)

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

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

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

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

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

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

2.1 數據描述

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

表1 GPS測站位置

2.2 小波基函數的選取

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

2.3 隨機項的提取

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

在對坐標時間序列進行小波多尺度分解后,對于包含隨機項的各層,其在不同尺度上的成分一般相互獨立,而且經常是平穩的,因此可以分別進行處理．

圖2 URUM站垂向坐標時間序列

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

2.4 周期項的提取

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

如圖3所示,由低頻部分a3和a4可以看出URUM站垂向坐標時間序列中存在1年周期項,并且為其主周期項,由d4可以看出半年周期項,由a5第0～100周時間段可以看出時間序列中存在2年周期項．這與快速傅里葉變換的結果也是相吻合的．雖然小波變換存在諸多優勢,但分析發現很難確定主要周期項的振幅,而且對于一些測站存在長周期如10年、12年,以及短周期項如季節周期項、月周期項的情況,小波變換優勢并不明顯．因此,在小波變換的基礎上,對分解得到的高頻或低頻信息進行傅里葉變換,提取其周期項及相應的振幅信息．如圖4在對BJFS站N方向上坐標時間序列進行小波分解時,雖然能根據高頻部分d5看出周期為1年的周期項,但是其他周期項則不明顯．鑒于快速傅里葉變換頻率特性較好,這里對低頻部分a5進行傅里葉變換,可以得到周期為10年,振幅為1.44 mm的周期項,以及周期為3年,振幅為1.39 mm的周期,如圖5所示．對d5進行快速傅里葉變換可以確定其1年周期項的振幅為0.80 mm．以此類推,將BJFS站N方向小波分解得到的各層主要周期項及振幅進行統計結果如表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

采用上述數據處理流程和方法,本文對所選測站進行了處理分析,鑒于篇幅,只對小波分解后的各層主要周期項及振幅變化進行了統計,結果如表3所示．

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

3 結 論

本文依據GPS測站殘差時間序列非線性變化的特點以及小波變換和傅里葉變換在坐標時間序列分析中各自的優勢,提出了一種小波分析和傅里葉變換相結合的時頻分析方法,來提取坐標時間序列中的特征信息,得到了以下結論:

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

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

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