1、時間序列-小波分析時間序列是地學研究中經常遇到的問題。在時間序列研究中，時域和頻域是常用的兩種（如Fourier 變換）雖具有準確的頻率定位功能，但僅適合平穩時間序列分析。然而，地學中許多現象（如河川徑流、地震波、暴雨、洪水等）隨時間的變化往往受到多種因素的綜合影響，大都屬于非平穩序列，它 們20 世紀 80 年代初，由Morlet 為更在不同時間尺度中的變化趨勢，并能對系統未來發展趨勢進行定性估計。形維數的計算，突變點的監測和周期成分的識別以及多時間尺度的分析等。一、小波分析基本原理小波函數 (tL2(R且滿足： (t)dt 0（1）式中， (t) 為基小波函數，它可通過尺度的伸縮和時間軸上

2、的平移構成一簇函數系：a,b(t) a(tb)其中bR,a 0（2）a ab (t) 為平移因子，反應時間上的平移。基小波函數的好壞，并由此選定該類研究所需的基小波函數。小波變換若(t) 是由（2）式給出的子小波，對于給定的能量有限信號f(t)L2(R)，其連續小波變換a , b（Continue Wavelet Transform，簡寫為 CWT）為：W (a, b) a -1/ 2 t （3）fRa式中，Wf(a,b)ab(xb)a為(xb)f(t)（k=1,2; ta為取樣間隔，則式）的離散小波變換形式為：2Nkt -2W (a, b) at (fak1)（4）由式（3）或（4）可知小波

3、分析的基本原理，即通過增加或減小伸縮尺度a 來得到信號的低頻或高頻信息，然后分析信號的概貌或細節，實現對信號不同時間尺度和空間局部特征的分析。頻變化特征。小波方差2將小波系數的平方值在b 域上積分，就可得到小波方差，即2(5)Var(a) Wf(a, b) db小波方差隨尺度 a 的變化過程，稱為小波方差圖。由式(5)可知，它能反映信號波動的能量隨尺度 a 的分布。因此，小波方差圖可用來確定信號中不同種尺度擾動的相對強度和存在的主要時間尺度，即主周期。二、小波分析實例-時間序列的多時間尺度分析(Multi-time scale analysis)例題尺度的變化周期之中。也就是說，徑流變化在時間

4、域中存在多層次的時間尺度結構和局部變化特征。表1 給出了某流域某水文觀測站1966-2004 年的實測徑流數據Matlab6.5、suffer8.0 和相關軟件 等），完成下述任務：計算小波系數；繪制小波系數圖（實部、模和模方）、小波方差圖和主周期變化趨勢圖，并分別說明各圖在分析徑流多時間尺度變化特征中的作用。表 1 某流域某水文觀測站 1966-2004 年實測徑流數據（108m3）年份徑流量年份徑流量年份徑流量年份徑流量年份徑流量19661.43819742.23519820.77419901.80619981.70919671.15119754.37419830.36719910.449

5、19990.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分析選擇合適的基小波函數是前提在運用小波分析理論解決實際問題

6、時，選擇合適的基小波函數是前提。只有選擇了適合具體問題的基 小波函數，才能得到較為理想的結果。目前，可選用的小波函數很多，如Mexican hat 小波、Haar 小波、Morlet 小波和 Meyer 小波等。在本例中，我們選用 Morlet 連續復小波變換來分析徑流時間序列的多時間尺度特征。原因如下：行此項分析。和振幅兩方面的信息，有利于對問題的進一步分析。復小波函數的實部和虛部位相差為使分析結果更為準確。繪制小波系數圖、小波方差圖和主周期變化趨勢圖是關鍵存在的多時間尺度。具體步驟數據格式的轉化邊界效應的消除或減小計算小波系數計算復小波系數的實部繪制小波系數實部等值線圖繪制小波系數模和模方

7、等值線圖繪制小波方差圖繪制主周期趨勢圖下面，我們以上題為例，結合軟件Matlab Suffer 8.0 和 Excel，詳細說明小波系數的計算和各圖形的繪制過程，并分別說明各圖在分析徑流多時間尺度變化特征中的作用。數據格式的轉化和保存將存放在 Excel（以時間為序排為一列轉化為 Matlab 6.5（.mat） 并存盤。具體操作為：在 Matlab 6.5 界面下，單擊“File-Import Data需要轉化的數據文件（本例的文件名為runoff.xl），“ ”，出現圖 1 顯示界面；然后雙擊圖1 中的 RunoffEditor: ”對話框，選擇File 文件夾下的“Workspace 2

8、 所示的“Save to to MAT-File”窗口。圖1 數據格式的轉化圖2數據的保存邊界效應的消除或減小因為本例中的實測徑流數據為有限時間數據序列，在時間序列的兩端可能會產生“邊界效用”。為消除 Matlab 6.5 小波工具箱中的信號延伸 Extension）功能，對徑流數據兩端進行對稱性延伸。具體方法為：在 Matlab 6.5 界面的“Command 鍵彈“Wavelet Toolbox Main 小波工具箱主菜單（然后單擊“Signal 打開 Signal Extension / Truncationrunoff.mat 文件單擊“打4 Matlab 5的n Mode菜單下包含了

9、6 、 PeriodicZero ContinuousSmooth and For Direction to extend3 種延伸模式（Both、Left and 在這里我們選擇對稱性兩端延伸進行計算。數據延伸的具體操作過程是： 在 Extension Mode 下選擇“ Symmetric”，Dircetion to extend 下選擇“Both”，性ee d 文件。從 erunoff 文件可知，系統自動將原時間序列數據向前對稱延伸12 個單位，向后延伸 13 個單位。圖3 小波工具箱主菜單圖4 徑流時間序列的延伸計算小波系數選擇 Matlab 6.5 小波工具箱中的 Morlet 復小

10、波函數對延伸后的徑流數據序列（erunoff.mat）進行小波變換，計算小波系數并存盤。小波工具箱主菜單界面見圖 3，單擊“ Wavelet 1-D ” 下的子菜單“ Complex Continuous Wavelet File“Load Signal”按鈕，載入徑流時間序列 erunoff.mat（圖 5）。圖 5 的左側為信號顯示區域，右側區域給出了信號序列和復小波變換的有關信息和參數，主要圖 5 小波變換菜單界面（Data shanfbsp 和 （Sampling Period）、周期設置 和運行按鈕以及顯示區域的相關顯示設置按鈕。本例中，我們選擇cmor )、取樣周期為 、最大尺度為

11、 Coefficientscerunoff.mat注意：上面涉及到的數據保存，其格式均為.mat。Morlet將復小波系數轉存到 Excel 表格，去掉兩端延伸數據的小波系數，并計算小波系數實部。在 Matlab 6.5 界面下的 Workspace 中將 cerunoff.matExcel （本例中去掉前 12 列和后 13 列ExcelExcelIMREAL（圖6）。圖 6 復小波系數及實部計算示意圖Excel 中 IMREAL 函數的調用方法為：單擊“插入”菜單下的“函數按鈕，彈出圖 7 所示的“插入（圖8）。在“Inumber”一欄的空白處輸入所要計算的數據（圖6），單擊“確定”即可得

12、到小波系數實部值。圖7IMEAL 函數調用圖8IMEAL 函數參數需要說明的是，從 cerunoff.mat 文件中轉到 Excel 里的復小波系數，在其實部和虛部中間包含許多“空格”，在計算之前需要先將其去掉。Suffer 8.0，繪制小波系數實部等值線圖小波系數實部等值線圖的繪制首先，將小波系數實部數據按照圖9 格式排列，其中列A 為時間，列B 為尺度，列C 為不同時間和尺度下所對應的小波系數實部值。圖 9 小波系數實部數據格式9 數據轉化成 Suffer 8.0Suffer 8.0 界面下，單擊“網格” 菜單下的“數據”按鈕，在“打開”窗口選擇要打開的文件（.xls），單擊“打開”后彈出

13、“網格化數 據”對話框（圖 10）們選擇“克里格“網格方法”，單擊“確定”，完成數據格式的轉化。圖10 小波系數實部數據格式轉化圖11Suffer8.0中的小系數實部等值線圖Suffer 8.0.grd（圖 11）。小波系數實部等值線圖在多時間尺度分析中的作用小波系數實部等值線圖能反映徑流序列不同時間尺度的周期變化及其在時間域中的分布，進而能判斷在不同時間尺度上，徑流的未來變化趨勢。為能比較清楚30LHHHL的說明小波系數實部等值線圖在徑流多時間尺度分析中的25LL作用，我們利用 Suffer 8.0 對其進一步處理和修飾，得到2015圖 12 顯示的小波系數實部等值線圖。其中，橫坐標為時間10HLHLHLL（年份實部值。當小波系數實部值為正時，代表徑流豐水期，在HLHH5(a)圖中我們用實線繪出，“H”表示正值中心；為負時，表示徑流枯水期，用虛線繪出，“