第六章時間序列的小波分析時間序列(TimeSeries)是地學研究中經常遇到的問題。在時間序列研究中，時域和頻域是常用的兩種基本形式。其中，時域分析具有時間定位能力，但無法得到關于時間序列變化的更多信息；頻域分析(如Fourier變換)雖具有準確的頻率定位功能，但僅適合平穩時間序列分析。然而，地學中許多現象(如河川徑流、地震波、暴雨、洪水等)隨時間的變化往往受到多種因素的綜合影響，大都屬于非平穩序列，它們不但具有趨勢性、周期性等特征，還存在隨機性、突變性以及“多時間尺度”結構，具有多層次演變規律。對于這類非平穩時間序列的研究，通常需要某一頻段對應的時間信息，或某一時段的頻域信息。顯然，時域分析和頻域分析對此均無能為力。20世紀80年代初，由Morlet提出的一種具有時-頻多分辨功能的小波分析(WaveletAnalysis)為更好的研究時間序列問題提供了可能，它能清晰的揭示出隱藏在時間序列中的多種變化周期，充分反映系統在不同時間尺度中的變化趨勢，并能對系統未來發展趨勢進行定性估計。目前，小波分析理論已在信號處理、圖像壓縮、模式識別、數值分析和大氣科學等眾多的非線性科學領域內得到了廣泛的應。在時間序列研究中，小波分析主要用于時間序列的消噪和濾波，信息量系數和分形維數的計算，突變點的監測和周期成分的識別以及多時間尺度的分析等。—、小波分析基本原理小波函數小波分析的基本思想是用一簇小波函數系來表示或逼近某一信號或函數。因此，小波函數是小波分析的關鍵，它是指具有震蕩性、能夠迅速衰減到零的一類函數，即小波函數V(t)eL2(R)且滿足：j+如(t)dt=0 (1)-s式中，V(t)為基小波函數，它可通過尺度的伸縮和時間軸上的平移構成一簇函數系：Vab(t)=|a|-1/2V(tab)其中，a,beR,a。0 (2)式中，V(t)為子小波；a為尺度因子，反映小波的周期長度；b為平移因子，反應時間上的平移。a,b需要說明的是，選擇合適的基小波函數是進行小波分析的前提。在實際應用研究中，應針對具體情況選擇所需的基小波函數；同一信號或時間序列，若選擇不同的基小波函數，所得的結果往往會有所差異，有時甚至差異很大。目前，主要是通過對比不同小波分析處理信號時所得的結果與理論結果的誤差來判定基小波函數的好壞，并由此選定該類研究所需的基小波函數。小波變換若V(t)是由(2)式給出的子小波，對于給定的能量有限信號f(t)eL2(R)，其連續小波變換a,b(ContinueWaveletTransform簡寫為CWT)為：W(a,b)=|a-1/2jf(t)v(^)dt (3)f R ax—b、式中，Wf(a,b)為小波變換系數;f(t)為一個信號或平方可積函數;a為伸縮尺度;b平移參數;V(二)x—b、為V(——)的復共軛函數。地學中觀測到的時間序列數據大多是離散的，設函數f(kAt)，(k=1,2,?,N;Ata為取樣間隔），則式（3）的離散小波變換形式為：W（a,b）=|a|-i/2At丈f（kAt偵（kAt~b） （4）k=1 a由式（3）或（4）可知小波分析的基本原理，即通過增加或減小伸縮尺度a來得到信號的低頻或高頻信息，然后分析信號的概貌或細節，實現對信號不同時間尺度和空間局部特征的分析。實際研究中，最主要的就是要由小波變換方程得到小波系數，然后通過這些系數來分析時間序列的時頻變化特征。小波方差將小波系數的平方值在b域上積分，就可得到小波方差，即Var（a）=廣|w（a,b）|2db （5）—sf小波方差隨尺度a的變化過程，稱為小波方差圖。由式（5）可知，它能反映信號波動的能量隨尺度a的分布。因此，小波方差圖可用來確定信號中不同種尺度擾動的相對強度和存在的主要時間尺度，即主周期。二、小波分析實例-時間序列的多時間尺度分析（Multi-timescaleanalysis）例題河川徑流是地理水文學研究中的一個重要變量，而多時間尺度是徑流演化過程中存在的重要特征。所謂徑流時間序列的多時間尺度是指：河川徑流在演化過程中，并不存在真正意義上的變化周期，而是其變化周期隨著研究尺度的不同而發生相應的變化，這種變化一般表現為小時間尺度的變化周期往往嵌套在大尺度的變化周期之中。也就是說，徑流變化在時間域中存在多層次的時間尺度結構和局部變化特征。表1給出了某流域某水文觀測站1966-2004年的實測徑流數據。試運用小波分析理論，借助Matlab6.5、suffer8.0和相關軟件（Excel等），完成下述任務：⑴計算小波系數；⑵繪制小波系數圖（實部、模和模方）、小波方差圖和主周期變化趨勢圖，并分別說明各圖在分析徑流多時間尺度變化特征中的作用。表1某流域某水文觀測站1966-2004年實測徑流數據（X108m3）年份徑流量年份徑流量年份徑流量年份徑流量年份徑流量19661.43819742.23519820.77419901.80619981.70919671.15119754.37419830.36719910.44919990.00019680.53619764.21919840.56219920.12020000.00019691.47019772.59019853.04019930.62720012.10419703.47619783.35019860.30419941.65820020.00919714.06819792.54019870.72819951.02520033.17719722.14719800.80719880.49219960.95520040.92119733.93119810.57319890.00719971.341分析選擇合適的基小波函數是前提在運用小波分析理論解決實際問題時，選擇合適的基小波函數是前提。只有選擇了適合具體問題的基小波函數，才能得到較為理想的結果。目前，可選用的小波函數很多，如Mexicanhat小波、Haar小波、Morlet小波和Meyer小波等。在本例中，我們選用Morlet連續復小波變換來分析徑流時間序列的多時間尺度特征。原因如下：1.1徑流演變過程中包含“多時間尺度”變化特征且這種變化是連續的，所以應采用連續小波變換來進行此項分析。

1.2實小波變換只能給出時間序列變化的振幅和正負，而復小波變換可同時給出時間序列變化的位相和振幅兩方面的信息，有利于對問題的進一步分析。1.3復小波函數的實部和虛部位相差為n/，2能夠消除用實小波變換系數作為判據而產生的虛假振蕩，使分析結果更為準確。繪制小波系數圖、小波方差圖和主周期變化趨勢圖是關鍵當選擇好合適的基小波函數后，下一步的關鍵就是如何通過小波變換獲得小波系數，然后利用相關軟件繪制小波系數圖、小波方差圖和主周期變化趨勢圖，進而根據上述三種圖形的變化識別徑流時間序列中存在的多時間尺度。具體步驟數據格式的轉化邊界效應的消除或減小計算小波系數計算復小波系數的實部繪制小波系數實部等值線圖繪制小波系數模和模方等值線圖繪制小波方差圖繪制主周期趨勢圖下面，我們以上題為例，結合軟件Matlab6.5、Suffer8.0和Excel，詳細說明小波系數的計算和各圖形的繪制過程，并分別說明各圖在分析徑流多時間尺度變化特征中的作用。數據格式的轉化和保存將存放在Excel表格里的徑流數據（以時間為序排為一列）轉化為Matlab6.5識別的數據格式（.mat）并存盤。具體操作為：在Matlab6.5界面下，單擊“File-ImportData”，出現文件選擇對話框“Import”后，找到需要轉化的數據文件（本例的文件名為runoff.xls），單擊“打開”。等數據轉化完成后，單擊“Finish”，出現圖1顯示界面；然后雙擊圖1中的Runoff，彈出“ArrayEditor:runoff”對話框，選擇File文件夾下的“SaveWorkspaceAs”單擊，出現圖2所示的“SavetoMAT-File:”窗口，選擇存放路徑并填寫文件名（runoff.mat），單擊“保存”并關閉“SavetoMAT-File”窗口。iq-rri?iq-rrie*圖1數據格式的轉化 圖2iq-rri?iq-rrie*圖1數據格式的轉化 圖2數據的保存□a?■+折口figmj■■研□jrxyriS?詢Hifirarl!!ski別Ik耳OhW_u ±J?i邊界效應的消除或減小因為本例中的實測徑流數據為有限時間數據序列，在時間序列的兩端可能會產生“邊界效用”。為消除或減小序列開始點和結束點附近的邊界效應，須對其兩端數據進行延伸。在進行完小波變換后，去掉兩

端延伸數據的小變換系數，保留原數據序列時段內的小波系數。本例中，我們利用Matlab6.5小波工具箱中的信號延伸(SignalExtension)功能，對徑流數據兩端進行對稱性延伸。具體方法為：在Matlab6.5界面的“CommandWindow”中輸入小波工具箱調用命令“Wavemenu”，按Enter鍵彈“WaveletToolboxMainMenu”(小波工具箱主菜單)界面(圖3)；然后單擊“SignalExtension”，打開SignalExtension/Truncation窗口，單擊“File”菜單下的“LoadSignal”，選擇runoff.mat文件單擊“打開”，出現圖4信號延伸界面。Matlab6.5的ExtensionMode菜單下包含了6種基本的延伸方式(Symmetric>Periodic>ZeroPadding、Continuous>SmoothandForSWT)和Directiontoextend菜單下的3種延伸模式(Both、LeftandRight)，在這里我們選擇對稱性兩端延伸進行計算。數據延伸的具體操作過程是：在ExtensionMode下選擇“Symmetric”，Dircetiontoextend下選擇“Both”，單擊“Extend”按鈕進行對稱性兩端延伸計算，然后單擊“File”菜單下的“SaveTranformedSignal”，將延伸后的數據結果存為erunoff.mat文件。從erunoff文件可知，系統自動將原時間序列數據向前對稱延伸12個單位，向后延伸13個單位。圖5圖5小波變換菜單界面圖3小波工具箱主菜單 圖4徑流時間序列的延伸計算小波系數選擇Matlab6.5小波工具箱中的Morlet復小波函數對延伸后的徑流數據序列(erunoff.mat)進行小波變換，計算小波系數并存盤。小波工具箱主菜單界面見圖3,單擊“Wavelet1-D”下的子菜單"ComplexContinuousWavelet1-D”，打開一維復連續小波界面，單擊“File”菜單下的“LoadSignal”按鈕，載入徑流時間序列erunoff.mat(圖5)。圖5的左側為信號顯示區域，右側區域給出了信號序列和復小波變換的有關信息和參數，主要包括數據長度(DataSize)、小波函數類型(Wavelet：cgau、shan、fbsp和cmor)、取樣周期(SamplingPeriod)、周期設置(ScaleSetting)和運行按鈕(Analyze)，以及顯示區域的相關顯示設置按鈕。本例中，我們選擇cmor(1-1.5)、取樣周期為1、最大尺度為32，單擊“Analyze”運行按鈕，計算小波系數。然后單擊“File”菜單下的“SaveCoefficients”，保存小波系數為cerunoff.mat文件。注意：上面涉及到的數據保存，其格式均為.mat。4.計算Morlet復小波系數的實部將復小波系數轉存到Excel表格，去掉兩端延伸數據的小波系數，并計算小波系數實部。在Matlab6.5界面下的Workspace中將cerunoff.mat文件導入,E“Command?掉前12列和dcniblearray(complex)doublearrayuhararraydciublearray在“a”上右擊，選擇“Graph”，在下拉菜單中選擇“plot”，即出小波方差圖FileEditViewInsertToolsWindowHelpset(gca,'XTickLabel',{'1963','1968','1973','1978','1983','1988','199^E'CcfemandWindows”中直接輸入函數:3','1998','2003'｝）使橫坐標顯示為年份。FileEditViewInsertToolsWindowHelpt”中更改X,Y坐標及圖名。, 1968 1973 1978 1983 1988 1993 1998在海)FigureHoset(gca,'XTickLabel',{'0','2004-01','2004-10','2005-08','2006-06','2007-04','2008-02','2008-12','2009-10','2010-0