小波分析實例

　　時間序列（Time Series）是地學研究中經常遇到的問題。在時間序列研究中，時域和頻域是常用的兩種基本形式。其中，時域分析具有時間定位能力，但無法得到關于時間序列變化的更多信息；頻域分析（如Fourier變換）雖具有準確的頻率定位功能，但僅適合平穩時間序列分析。

　　河川徑流是地理水文學研究中的一個重要變量，而多時間尺度是徑流演化過程中存在的重要特征。所謂徑流時間序列的多時間尺度是指：河川徑流在演化過程中，并不存在真正意義上的變化周期，而是其變化周期隨著研究尺度的不同而發生相應的變化，這種變化一般表現為小時間尺度的變化周期往往嵌套在大尺度的變化周期之中。也就是說，徑流變化在時間域中存在多層次的時間尺度結構和局部變化特征。

　　年降水量的變化趨勢分析

　　該站點的年降水量變化情況（LI_plot函數）如圖 0所示，其發展呈微微上升的趨勢，降水量最高年份出現在2003年，全年累計達到1610.70mm，最低值出現在1965年，累計僅有748.90mm。

　　具體步驟

　　1）數據格式的轉化

　　2）邊界效應的消除或減小

　　3）計算小波系數

　　4）計算復小波系數的實部

　　5）繪制小波系數實部等值線圖

　　6）繪制小波系數模和模方等值線圖

　　7）繪制小波方差圖

　　8） 繪制主周期趨勢圖

　　1. 數據格式的轉化和保存

　　將存放在Excel表格里的徑流數據（以時間為序排為一列）轉化為Matlab 6.5識別的數據格式（.mat）并存盤。

　　具體操作為：在Matlab 6.5 界面下，單擊“File-Import Data”，出現文件選擇對話框“Import”后，找到需要轉化的數據文件（本例的文件名為runoff.xls），單擊“打開”。等數據轉化完成后，單擊“Finish”，出現圖1顯示界面；然后雙擊圖1中的Runoff，彈出“Array Editor： runoff”對話框，選擇File文件夾下的“Save Workspace As”單擊，出現圖2所示的“Save to MAT-File：”窗口，選擇存放路徑并填寫文件名（runoff.mat），單擊“保存”并關閉“Save to MAT-File”窗口。

　　2．邊界效應的消除或減小

　　由于本例中的實測年降水量數據為有限時間數據序列，在時間序列的兩端可能會產生“邊界效用”。為消除或減小序列開始點和結束點附近的邊界效應，須對其兩端數據進行延伸。在進行完小波變換后，去掉兩端延伸數據的小波變換系數，保留原數據序列時段內的小波系數。本例中，我們利用Matlab小波工具箱中的信號延伸（SignalExtension）功能，對數據兩端進行對稱性延伸。

　　具體步驟：在Matlab界面的“Command Window”中輸入小波工具箱調用命令“wavemenu”，按Enter鍵彈“Wavelet Toolbox Main Menu”（小波工具箱主菜單）界面（圖 1上）；然后單擊“Signal Extension”，打開Signal Extension /Truncation窗口，單擊“File”菜單下的“Load Signal”，選擇Prec.mat文件單擊“打開”，出現信號延伸界面。Matlab的Extension Mode菜單下包含多種延伸方式和Direction to extend菜單下的3種延伸模式（Both、Left and Right），在這里我們選擇對稱性兩端延伸進行計算。具體操作過程是：在Extension Mode下選擇“Symmetric（Whole-Point）”，Dircetion to extend下選擇“Both”，單擊“Extend”按鈕進行對稱性兩端延伸計算（圖 1下），然后單擊“File”菜單下的“Save Tranformed Signal”，將延伸后的數據結果存為ePrec.mat文件。

　　從ePrec文件可知（entendvalue函數），系統自動將原時間序列數據向前對稱延伸5個單位，向后延伸6個單位。

　　圖 1

　　3．計算小波系數

　　選擇Matlab小波工具箱中的Morlet復小波函數對延伸后的數據序列（ePrec.mat）進行小波變換，計算小波系數并保存。

　　小波工具箱主菜單界面見圖1上，單擊“One-Dimensional”下的子菜單“Complex Continuous Wavelet 1-D”，打開一維復連續小波界面，單擊“File”菜單下的“Load Signal”按鈕，載入時間序列數據ePrec.mat。圖 2第一排的左側為信號顯示區域，右側區域給出了信號序列和復小波變換的有關信息和參數，主要包括數據長度（Data Size）、小波函數類型（Wavelet：cgau、shan、fbsp和cmor）、取樣周期（Sampling Period）、周期設置（ScaleSetting）和運行按鈕（Analyze），以及顯示區域的相關顯示設置按鈕。本例中，我們選擇cmor （1-1.5）、取樣周期為1、最大尺度為64（延伸后的時間序列的一半），單擊“Analyze”運行按鈕，計算小波系數。然后單擊“File”菜單下的“Save Coefficients”，保存小波系數為cePrec.mat文件。

　　圖 2

　　4．計算Morlet復小波系數的實部

　　去除兩端延伸數據的小波系數（entendvalue函數），并計算小波系數實部（real函數）。

　　5．繪制小波系數實部等值線圖

　　這部分過程應用LI_contourf函數，年降水量小波系數實部等值線圖見圖 3。

　　如圖3所示的小波系數實部等值線圖。其中，橫坐標為時間（年份），縱坐標為時間尺度，圖中的等值曲線為小波系數實部值。小波系數實部等值線圖能反映年降水量序列不同時間尺度的周期變化及其在時間域中的分布，進而能判斷在不同時間尺度上，年降水量的未來變化趨勢。

　　圖 3

　　由圖 3可以看出降水量1894~2010年演化過程中存在著14~18a的主振蕩周期（這一周期是多個振蕩周期疊加的矢量和，隨后將在小波方差圖中對這一主周期進行分解，剝離出第一主周期等等）。在整個時間尺度上出現4個偏多中心和3個偏少中心，分別為1899、1938、1972、2007年和1919、1956、1990年。

　　6．繪制小波系數模和模方等值線圖

　　首先，計算小波系數的模和模方，再利用LI_contourf函數繪制模和模方等值線圖。

　　Morlet小波系數的模值是不同時間尺度變化周期所對應的能量密度在時間域中分布的反映，系數模值愈大，表明其所對應時段或尺度的周期性就愈強。從圖 4可以看出，在降水量演化過程中，40~64年時間尺度模值最大（大于400），但1924~1994年之間模值小余200，說明在此時段內40~64年時間尺度的周期變化并不明顯，1995年之后模值再次增大，此時期后40~64年時間尺度的周期變化趨于顯著。

　　小波系數的模方相當于小波能量譜，它可以分析出不同周期的震蕩能量。由圖 5可知，40~64年時間尺度的能量最強、周期最顯著，但它的周期變化具有局部性（1924年之前和1994年之后）；10~15年時間尺度能量雖然較弱，但周期分布比較明顯，幾乎占據整個研究時域（1894~2010年）。

　　圖 4

　　圖 5

　　7．繪制小波方差圖

　　小波方差的計算過程參考：小波方差制作步驟（參考文獻存在一個錯誤，小波方差未除N，本貼小波方差計算已有修正）。

　　小波方差圖能反映降水量時間序列的波動能量隨尺度（a）的分布情況。可用來確定降水量演化過程中存在的主周期。

　　圖 6

　　小波方差圖中（圖 6）存在4個較為明顯的峰值，它們依次從小至大對應著8a、19a、30a和55a的時間尺度。其中，最大峰值對應著55a的時間尺度，說明55a左右（時間尺度）的周期震蕩最強，為年降水量變化的第一主周期（注意：不是55a是第一主周期）；30a時間尺度對應著第二峰值，為第二主周期；第三、第四峰值分別對應著19a和8a的時間尺度，它們依次為降水量的第三和第四主周期。這說明上述4個周期的波動控制著降水量在整個時間域內的變化特征。

　　8．主周期趨勢圖的繪制及其在多時間尺度分析中的作用

　　根據小波方差檢驗的結果，繪制演變的第一和第二主周期小波系數圖。

　　從主周期趨勢圖中可以分析出在不同的時間尺度下，降水量存在的平均周期及豐-枯變化特征。圖 7顯示，在55a特征時間尺度上，降水量變化的平均周期為35a左右，大約經歷了3個豐—枯轉換期；而在30a特征時間尺度上（圖 8），平均變化周期為20a左右，大約6個周期的豐—枯變化。

　　圖 7

　　圖 8