1、時間序列的小波分析時間序列（Time Series）是地學研究中經常遇到的問題。在時間序列研究中，時域和頻域是常用的兩種基本形式。其中，時域分析具有時間定位能力，但無法得到關于時間序列變化的更多信息；頻域分析（如Fourier變換）雖具有準確的頻率定位功能，但僅適合平穩時間序列分析。然而，地學中許多現象（如河川徑流、地震波、暴雨、洪水等）隨時間的變化往往受到多種因素的綜合影響，大都屬于非平穩序列，它們不但具有趨勢性、周期性等特征，還存在隨機性、突變性以及“多時間尺度”結構，具有多層次演變規律。對于這類非平穩時間序列的研究，通常需要某一頻段對應的時間信息，或某一時段的頻域信息。顯然，時域分析和頻

2、域分析對此均無能為力。20世紀80年代初，由Morlet提出的一種具有時-頻多分辨功能的小波分析（Wavelet Analysis）為更好的研究時間序列問題提供了可能，它能清晰的揭示出隱藏在時間序列中的多種變化周期，充分反映系統在不同時間尺度中的變化趨勢，并能對系統未來發展趨勢進行定性估計。目前，小波分析理論已在信號處理、圖像壓縮、模式識別、數值分析和大氣科學等眾多的非線性科學領域內得到了廣泛的應。在時間序列研究中，小波分析主要用于時間序列的消噪和濾波，信息量系數和分形維數的計算，突變點的監測和周期成分的識別以及多時間尺度的分析等。一、小波分析基本原理1. 小波函數小波分析的基本思想是用一簇小

3、波函數系來表示或逼近某一信號或函數。因此，小波函數是小波分析的關鍵，它是指具有震蕩性、能夠迅速衰減到零的一類函數，即小波函數且滿足： （1）式中，為基小波函數，它可通過尺度的伸縮和時間軸上的平移構成一簇函數系： 其中， （2）式中，為子小波；a為尺度因子，反映小波的周期長度；b為平移因子，反應時間上的平移。需要說明的是，選擇合適的基小波函數是進行小波分析的前提。在實際應用研究中，應針對具體情況選擇所需的基小波函數；同一信號或時間序列，若選擇不同的基小波函數，所得的結果往往會有所差異，有時甚至差異很大。目前，主要是通過對比不同小波分析處理信號時所得的結果與理論結果的誤差來判定基小波函數的好壞，并

4、由此選定該類研究所需的基小波函數。2. 小波變換若是由（2）式給出的子小波，對于給定的能量有限信號，其連續小波變換（Continue Wavelet Transform，簡寫為CWT）為： （3）式中，為小波變換系數；f(t)為一個信號或平方可積函數；a為伸縮尺度；b平移參數；為的復共軛函數。地學中觀測到的時間序列數據大多是離散的，設函數，（k=1,2,N; 為取樣間隔），則式（3）的離散小波變換形式為： （4）由式（3）或（4）可知小波分析的基本原理，即通過增加或減小伸縮尺度a來得到信號的低頻或高頻信息，然后分析信號的概貌或細節，實現對信號不同時間尺度和空間局部特征的分析。實際研究中，最主要

5、的就是要由小波變換方程得到小波系數，然后通過這些系數來分析時間序列的時頻變化特征。3. 小波方差將小波系數的平方值在b域上積分，就可得到小波方差，即 (5)小波方差隨尺度a的變化過程，稱為小波方差圖。由式(5)可知，它能反映信號波動的能量隨尺度a的分布。因此，小波方差圖可用來確定信號中不同種尺度擾動的相對強度和存在的主要時間尺度，即主周期。二、小波分析實例-時間序列的多時間尺度分析(Multi-time scale analysis)例題河川徑流是地理水文學研究中的一個重要變量，而多時間尺度是徑流演化過程中存在的重要特征。所謂徑流時間序列的多時間尺度是指：河川徑流在演化過程中，并不存在真正意義

6、上的變化周期，而是其變化周期隨著研究尺度的不同而發生相應的變化，這種變化一般表現為小時間尺度的變化周期往往嵌套在大尺度的變化周期之中。也就是說，徑流變化在時間域中存在多層次的時間尺度結構和局部變化特征。表1給出了某流域某水文觀測站1966-2004年的實測徑流數據。試運用小波分析理論，借助Matlab R2012a、suffer 12.0和其他相關軟件（Excel、記事本等），完成下述任務：（1）計算小波系數；（2）繪制小波系數圖（實部、模和模方）、小波方差圖和主周期變化趨勢圖，并分別說明各圖在分析徑流多時間尺度變化特征中的作用。表1 某流域某水文觀測站1966-2004年實測徑流數據（

7、15;108m3）年份徑流量年份徑流量年份徑流量年份徑流量年份徑流量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

8、.17719722.14719800.80719880.49219960.95520040.92119733.93119810.57319890.00719971.341分析1. 選擇合適的基小波函數是前提在運用小波分析理論解決實際問題時，選擇合適的基小波函數是前提。只有選擇了適合具體問題的基小波函數，才能得到較為理想的結果。目前，可選用的小波函數很多，如Mexican hat小波、Haar小波、Morlet小波和Meyer小波等。在本例中，我們選用Morlet連續復小波變換來分析徑流時間序列的多時間尺度特征。原因如下：1.1 徑流演變過程中包含“多時間尺度”變化特征且這種變化是連續的，所以應

9、采用連續小波變換來進行此項分析。1.2實小波變換只能給出時間序列變化的振幅和正負，而復小波變換可同時給出時間序列變化的位相和振幅兩方面的信息，有利于對問題的進一步分析。1.3 復小波函數的實部和虛部位相差為/2，能夠消除用實小波變換系數作為判據而產生的虛假振蕩，使分析結果更為準確。2. 繪制小波系數圖、小波方差圖和主周期變化趨勢圖是關鍵當選擇好合適的基小波函數后，下一步的關鍵就是如何通過小波變換獲得小波系數，然后利用相關軟件繪制小波系數圖、小波方差圖和主周期變化趨勢圖，進而根據上述三種圖形的變化識別徑流時間序列中存在的多時間尺度。具體步驟1. 數據格式的轉化2. 邊界效應的消除或減小3. 計算

10、小波系數4. 計算復小波系數的實部、模、模方、方差5. 繪制小波系數實部、模、模方等值線圖6. 繪制小波方差圖7. 繪制主周期趨勢圖下面，我們以上題為例，結合軟件Matlab R2012a、suffer 12.0、Excel、記事本等，詳細說明小波系數的計算和各圖形的繪制過程，并分別說明各圖在分析徑流多時間尺度變化特征中的作用。1. 數據格式的轉化和保存將存放在Excel表格里的徑流數據（以時間為序排為一列）轉化為Matlab R2012a識別的數據格式（.mat）并存盤。具體操作為：在Matlab R2012a 界面下，單擊“File-Import Data”，出現文件選擇對話框“Impor

11、t”后，找到需要轉化的數據文件（本例的文件名為runoff.xls），單擊“打開”。等數據轉化完成后，單擊“Finish”，出現圖1顯示界面；然后雙擊圖1中的Runoff，彈出“Array Editor: runoff”對話框，選擇File文件夾下的“Save Workspace As”單擊，出現圖2所示的“Save to MAT-File:”窗口，選擇存放路徑并填寫文件名（runoff.mat），單擊“保存”并關閉“Save to MAT-File”窗口。圖1 數據格式的轉化圖2數據的保存2. 邊界效應的消除或減小因為本例中的實測徑流數據為有限時間數據序列，在時間序列的兩端可能會產生“邊界效

12、用”。為消除或減小序列開始點和結束點附近的邊界效應，須對其兩端數據進行延伸。在進行完小波變換后，去掉兩端延伸數據的小變換系數，保留原數據序列時段內的小波系數。本例中，我們利用Matlab R2012a小波工具箱中的信號延伸（Signal Extension）功能，對徑流數據兩端進行對稱性延伸。具體方法為：在Matlab R2012a界面的“Command Window”中輸入小波工具箱調用命令“Wavemenu”，按Enter鍵彈“Wavelet Toolbox Main Menu”（小波工具箱主菜單）界面（圖）；然后單擊“Signal Extension”，打開Signal Extensio

13、n / Truncation窗口，單擊“File”菜單下的“Load Signal”，選擇runoff.mat文件單擊“打開”，出現圖4信號延伸界面。Matlab R2012a的Extension Mode菜單下包含了6種基本的延伸方式（Symmetric、Periodic、Zero Padding、Continuous、Smooth and For SWT）和Direction to extend菜單下的3種延伸模式（Both、Left and Right），在這里我們選擇對稱性兩端延伸進行計算。數據延伸的具體操作過程是：Desired Length可以任意選，只要比原始信號長度大，建議在原

14、始信號的基礎上加20（這樣左右對稱地延伸10個數據），這里選擇默認的64；Dircetion to extend下選擇“Both”；Extension Mode下選擇“Symmetric”；單擊“Extend”按鈕進行對稱性兩端延伸計算，然后單擊“File”菜單下的“Save Tranformed Signal”，將延伸后的數據結果存為erunoff.mat文件。從erunoff文件可知，系統自動將原時間序列數據向前對稱延伸12個單位，向后延伸13個單位。圖3 小波工具箱主菜單圖4 徑流時間序列的延伸3. 計算小波系數圖5 小波變換菜單界面選擇Matlab R2012a小波工具箱中的Morle

15、t復小波函數對延伸后的徑流數據序列（erunoff.mat）進行小波變換，計算小波系數并存盤。小波工具箱主菜單界面見圖3，單擊“Wavelet 1-D”下的子菜單“Complex Continuous Wavelet 1-D”，打開一維復連續小波界面，單擊“File”菜單下的“Load Signal”按鈕，載入徑流時間序列erunoff.mat（圖5）。圖5的左側為信號顯示區域，右側區域給出了信號序列和復小波變換的有關信息和參數，主要包括數據長度（Data Size）、小波函數類型（Wavelet：cgau、shan、fbsp和cmor）、取樣周期（Sampling Period）、周期設置（

16、Scale Setting）和運行按鈕（Analyze），以及顯示區域的相關顯示設置按鈕。本例中，我們選擇cmor (1-1.5)、取樣周期為1、最大尺度為32，單擊“Analyze”運行按鈕，計算小波系數。然后單擊“File”菜單下的“Save Coefficients”，保存小波系數為cerunoff.mat文件。4. 計算Morlet復小波系數的實部、模、模方、方差在Matlab R2012a界面下的Workspace中將cerunoff.mat文件導入，見圖6。圖6 小波系數導入到Matlab然后雙擊“coefs”打開，刪掉掉延伸數據的小波變換系數（本例中去掉前12列和后13列），保存

17、。接下來開始計算Morlet復小波系數的實部、模、模方、方差，具體操作為：在“Command Windows”中直接輸入函數“shibu=real(coefs);”，點擊“回車”鍵，計算實部；輸入函數“mo=abs(coefs);”，點擊“回車”鍵，計算模；輸入函數“mofang=(mo).2;”，點擊“回車”鍵，計算模方；輸入函數“fangcha=sum(abs(coefs).2,2);”，點擊“回車”鍵，計算方差。見圖7。圖7計算出的實部、模、模方、方差成果注意：上面涉及到的數據保存，其格式均為.mat。5. 繪制小波系數實部、模、模方等值線圖實部、模、模方等值線圖的繪制方法一樣，這里僅以

18、實部等值線圖為例。5.1 小波系數實部等值線圖的繪制首先，將小波系數實部數據復制到Excel中按照圖8格式排列，其中列A為時間，列B為尺度，列C為不同時間和尺度下所對應的小波系數實部值。其次，將圖9數據轉化成Suffer 12.0識別的數據格式。具體操作為：在Surfer 12.0界面下，單擊“網格”菜單下的“數據”按鈕，在“打開”窗口選擇要打開的文件（小波系數實部.xls），單擊“打開”后彈出“網格化數據”對話框（圖10）。它給出了多種不同的網格化方法、文件輸出路徑及網格線索幾何學等信息。這里我們選擇“克里格“網格方法”，單擊“確定”，完成數據格式的轉化。圖9 Suffer 12.0可以識別

19、的數據格式列表圖8 小波系數實部數據格式圖10 小波系數實部數據格式轉化圖11 Suffer8.0中的小系數實部等值線圖最后，繪制小波系數實部等值線圖。在Surfer 12.0界面下，單擊“地圖”菜單下的“等值線圖-新建等值線圖”按鈕，彈出“打開網格”窗口后，選擇“小波系數實部.grd”文件，單擊“打開”，完成等值線圖的繪制并保存（圖11）。5.2 小波系數實部等值線圖在多時間尺度分析中的作用圖12 小系數實部等值線圖小波系數實部等值線圖能反映徑流序列不同時間尺度的周期變化及其在時間域中的分布，進而能判斷在不同時間尺度上，徑流的未來變化趨勢。為能比較清楚的說明小波系數實部等值線圖在徑流多時間尺

20、度分析中的作用，我們利用Surfer 12.0對其進一步處理和修飾，得到圖12顯示的小波系數實部等值線圖。其中，橫坐標為時間（年份），縱坐標為時間尺度，圖中的等值曲線為小波系數實部值。當小波系數實部值為正時，代表徑流豐水期，在圖中我們用實線繪出，“H”表示正值中心；為負時，表示徑流枯水期，用虛線繪出，“L”表示負值中心。由圖12可以清楚的看到徑流演化過程中存在的多時間尺度特征?？偟膩碚f，在流域徑流演變過程中存在著1832年，817年以及37年的3類尺度的周期變化規律。其中，在1832年尺度上出現了枯-豐交替的準兩次震蕩；在817年時間尺度上存在準5次震蕩。同時，還可以看出以上兩個尺度的周期變化

21、在整個分析時段表現的非常穩定，具有全域性；而310年尺度的周期變化，在1980s以后表現的較為穩定。5.3 小波系數模和模方等值線圖的繪制參考5.1，繪制小波系數模和模方等值線圖（圖13、14）。圖13 小波系數模等值線圖圖14 小波系數模方等值線圖5.4 小波系數模等值線圖在多時間尺度分析中的作用Morlet小波系數的模值是不同時間尺度變化周期所對應的能量密度在時間域中分布的反映，系數模值愈大，表明其所對應時段或尺度的周期性就愈強。從圖13可以看出，在流域徑流演化過程中，1832年時間尺度模值最大，說明該時間尺度周期變化最明顯，1822年時間尺度的周期變化次之，其他時間尺度的周期性變化較小；

22、 5.5 小波系數模方等值線圖在多時間尺度分析中的作用小波系數的模方相當于小波能量譜，它可以分析出不同周期的震蕩能量。由圖14知，2532年時間尺度的能量最強、周期最顯著，但它的周期變化具有局部性（1980s前）；1015年時間尺度能量雖然較弱，但周期分布比較明顯，幾乎占據整個研究時域(19742004年)。6. 繪制小波方差圖6.1小波方差圖的繪制在圖7的“fangcha”上右擊，選擇“Graph”，在下拉菜單中選擇“plot”，即出小波方差圖，見圖15，在Matlab中可繼續美化。也可雙擊“fangcha”，將數據復制到其他軟件（如Excel）中，以小波方差為縱坐標，時間尺度a為橫坐標，繪

23、制小波方差，如圖16。圖15 Matlab繪制的小波方差圖圖16 小波方差圖6.2小波方差圖在多時間尺度分析中的作用小波方差圖能反映徑流時間序列的波動能量隨尺度a的分布情況?？捎脕泶_定徑流演化過程中存在的主周期。流域徑流的小波方差圖中（圖15）存在4個較為明顯的峰值，它們依次對應著28年、14年、8年和4年的時間尺度。其中，最大峰值對應著28年的時間尺度，說明28年左右的周期震蕩最強，為流域年徑流變化的第一主周期；14年時間尺度對應著第二峰值，為徑流變化的第二主周期，第三、第三峰值分別對應著8年和4年的時間尺度，它們依次為流域徑流的第三和第四主周期。這說明上述4個周期的波動控制著流域徑流在整個

24、時間域內的變化特征。7. 主周期趨勢圖的繪制及其在多時間尺度分析中的作用根據小波方差檢驗的結果，我們繪制出了控制流域徑流演變的第一和第二主周期小波系數圖（圖17）。從主周期趨勢圖中我們可以分析出在不同的時間尺度下，流域徑流存在的平均周期及豐-枯變化特征。圖16a顯示，在14年特征時間尺度上，流域徑流變化的平均周期為9.5年左右，大約經歷了4個豐-枯轉換期；而在28年特征時間尺度上（圖16b），流域的平均變化周期為20年左右，大約2個周期的豐-枯變化。圖17 大沽夾河流域年徑流變化的13年和28年特征時間尺度小波實部過程線參考文獻王文圣，丁晶，李耀清. 2005. 水文小波分析M. 北京：化學工

25、業出版社曹素華等. 1998. 實用醫學多因素統計方法M. 上海：上海醫科大學出版社方開泰. 1989. 實用多元統計分析M. 上海：華東師范大學出版社何清波，蘇炳華，錢亢. 2002. 醫學統計學及其軟件包M. 上海：上海科學技術文獻出版社胡秉民. 1987. 微電腦在農業科學中的應用M. 北京：科學出版社孫尚拱. 1990. 實用多元變量統計方法與計算程序M. 北京：北京醫科大學、中國協和醫科大學聯合出版社唐守正. 1986. 多元統計分析方法M。北京：中國林業出版社王學仁. 1982. 地址數據的多變量統計分析. 北京：科學出版社徐振邦，金淳浩，婁元仁. 1986. 距離系數和距離系數尺度在聚類分析中的應用M. 趙旭東等主編，中國數學地質（1）. 北京：地質出版社於崇文. 1978. 數學地質的方法與應用M. 北