小波分析實(shí)例

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

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

　　年降水量的變化趨勢(shì)分析

　　該站點(diǎn)的年降水量變化情況（LI_plot函數(shù)）如圖 0所示，其發(fā)展呈微微上升的趨勢(shì)，降水量最高年份出現(xiàn)在2003年，全年累計(jì)達(dá)到1610.70mm，最低值出現(xiàn)在1965年，累計(jì)僅有748.90mm。

　　具體步驟

　　1）數(shù)據(jù)格式的轉(zhuǎn)化

　　2）邊界效應(yīng)的消除或減小

　　3）計(jì)算小波系數(shù)

　　4）計(jì)算復(fù)小波系數(shù)的實(shí)部

　　5）繪制小波系數(shù)實(shí)部等值線圖

　　6）繪制小波系數(shù)模和模方等值線圖

　　7）繪制小波方差圖

　　8） 繪制主周期趨勢(shì)圖

　　1. 數(shù)據(jù)格式的轉(zhuǎn)化和保存

　　將存放在Excel表格里的徑流數(shù)據(jù)（以時(shí)間為序排為一列）轉(zhuǎn)化為Matlab 6.5識(shí)別的數(shù)據(jù)格式（.mat）并存盤。

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

　　2．邊界效應(yīng)的消除或減小

　　由于本例中的實(shí)測(cè)年降水量數(shù)據(jù)為有限時(shí)間數(shù)據(jù)序列，在時(shí)間序列的兩端可能會(huì)產(chǎn)生“邊界效用”。為消除或減小序列開始點(diǎn)和結(jié)束點(diǎn)附近的邊界效應(yīng)，須對(duì)其兩端數(shù)據(jù)進(jìn)行延伸。在進(jìn)行完小波變換后，去掉兩端延伸數(shù)據(jù)的小波變換系數(shù)，保留原數(shù)據(jù)序列時(shí)段內(nèi)的小波系數(shù)。本例中，我們利用Matlab小波工具箱中的信號(hào)延伸（SignalExtension）功能，對(duì)數(shù)據(jù)兩端進(jìn)行對(duì)稱性延伸。

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

　　從ePrec文件可知（entendvalue函數(shù)），系統(tǒng)自動(dòng)將原時(shí)間序列數(shù)據(jù)向前對(duì)稱延伸5個(gè)單位，向后延伸6個(gè)單位。

　　圖 1

　　3．計(jì)算小波系數(shù)

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

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

　　圖 2

　　4．計(jì)算Morlet復(fù)小波系數(shù)的實(shí)部

　　去除兩端延伸數(shù)據(jù)的小波系數(shù)（entendvalue函數(shù)），并計(jì)算小波系數(shù)實(shí)部（real函數(shù)）。

　　5．繪制小波系數(shù)實(shí)部等值線圖

　　這部分過程應(yīng)用LI_contourf函數(shù)，年降水量小波系數(shù)實(shí)部等值線圖見圖 3。

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

　　圖 3

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

　　6．繪制小波系數(shù)模和模方等值線圖

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

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

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

　　圖 4

　　圖 5

　　7．繪制小波方差圖

　　小波方差的計(jì)算過程參考：小波方差制作步驟（參考文獻(xiàn)存在一個(gè)錯(cuò)誤，小波方差未除N，本貼小波方差計(jì)算已有修正）。

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

　　圖 6

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

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

　　根據(jù)小波方差檢驗(yàn)的結(jié)果，繪制演變的第一和第二主周期小波系數(shù)圖。

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

　　圖 7

　　圖 8