(完整word版)時間序列的小波分析及等值線圖、小波方差制作

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

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

3、表示或逼近某一信號或函數(shù)。因此，小波函數(shù)是小波分析的關(guān)鍵，它是指具有震蕩性、能夠迅速衰減到零的一類函數(shù)，即小波函數(shù)y(t)eU(R)且滿足：TOC o 1-5 h zU(t)dt二0(1)g式中，屮(t)為基小波函數(shù)，它可通過尺度的伸縮和時間軸上的平移構(gòu)成一簇函數(shù)系：tb屮(t)=|a|t/2屮()其中，a,beR,a豐0(2)a,ba式中，屮(t)為子小波；a為尺度因子，反映小波的周期長度；b為平移因子，反應(yīng)時間上的平移。a,b需要說明的是，選擇合適的基小波函數(shù)是進行小波分析的前提。在實際應(yīng)用研究中，應(yīng)針對具體情況選擇所需的基小波函數(shù)；同一信號或時間序列，若選擇不同的基小波函數(shù)，所得的結(jié)果往

4、往會有所差異，有時甚至差異很大。目前，主要是通過對比不同小波分析處理信號時所得的結(jié)果與理論結(jié)果的誤差來判定基小波函數(shù)的好壞，并由此選定該類研究所需的基小波函數(shù)。小波變換若屮(t)是由(2)式給出的子小波，對于給定的能量有限信號f(t)eU(R)，其連續(xù)小波變換a,b(ContinueWaveletTransform簡寫為CWT)為：_tbTOC o 1-5 h zW(a,b)=la-2ff(t)7()dt(3)fraxb式中，Wf(a,b)為小波變換系數(shù);f(t)為一個信號或平方可積函數(shù);a為伸縮尺度;b平移參數(shù);屮()faxb為屮()的復(fù)共軛函數(shù)。地學(xué)中觀測到的時間序列數(shù)據(jù)大多是離散的,設(shè)函

5、數(shù)f(kAt)，(k=l,2,N;Ata為取樣間隔)，則式(3)的離散小波變換形式為：W(a,b)=|a|-i/2At丈f(kAt(kt-b)(4)fak=1由式(3)或(4)可知小波分析的基本原理，即通過增加或減小伸縮尺度a來得到信號的低頻或高頻信息，然后分析信號的概貌或細節(jié)，實現(xiàn)對信號不同時間尺度和空間局部特征的分析。實際研究中，最主要的就是要由小波變換方程得到小波系數(shù)，然后通過這些系數(shù)來分析時間序列的時頻變化特征。小波方差將小波系數(shù)的平方值在b域上積分，就可得到小波方差，即Var(a)=代|W(a,b)db(5)gf小波方差隨尺度a的變化過程，稱為小波方差圖。由式(5)可知，它能反映信號

6、波動的能量隨尺度a的分布。因此，小波方差圖可用來確定信號中不同種尺度擾動的相對強度和存在的主要時間尺度，即主周期。二、小波分析實例-時間序列的多時間尺度分析(Multi-timescaleanalysis)例題河川徑流是地理水文學(xué)研究中的一個重要變量，而多時間尺度是徑流演化過程中存在的重要特征。所謂徑流時間序列的多時間尺度是指：河川徑流在演化過程中，并不存在真正意義上的變化周期，而是其變化周期隨著研究尺度的不同而發(fā)生相應(yīng)的變化，這種變化一般表現(xiàn)為小時間尺度的變化周期往往嵌套在大尺度的變化周期之中。也就是說，徑流變化在時間域中存在多層次的時間尺度結(jié)構(gòu)和局部變化特征。表1給出了某流域某水文觀測站1

7、966-2004年的實測徑流數(shù)據(jù)。試運用小波分析理論，借助MatlabR2012a、suffer12.0和其他相關(guān)軟件(Excel、記事本等)，完成下述任務(wù)：(1)計算小波系數(shù)；(2)繪制小波系數(shù)圖(實部、模和模方)、小波方差圖和主周期變化趨勢圖，并分別說明各圖在分析徑流多時間尺度變化特征中的作用。表1某流域某水文觀測站1966-2004年實測徑流數(shù)據(jù)(X108m3)年份徑流量年份徑流量年份徑流量年份徑流量年份徑流量19661.43819742.23519820.77419901.80619981.70919671.15119754.37419830.36719910.44919990.000

8、19680.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分析選擇合適的基小波函數(shù)是前提在運用小波分析理論解決實際問題時，選擇合適的基小

9、波函數(shù)是前提。只有選擇了適合具體問題的基小波函數(shù)，才能得到較為理想的結(jié)果。目前，可選用的小波函數(shù)很多，如Mexicanhat小波、Haar小波、Morlet小波和Meyer小波等。在本例中，我們選用Morlet連續(xù)復(fù)小波變換來分析徑流時間序列的多時間尺度特征。原因如下：徑流演變過程中包含“多時間尺度”變化特征且這種變化是連續(xù)的，所以應(yīng)采用連續(xù)小波變換來進行此項分析。實小波變換只能給出時間序列變化的振幅和正負，而復(fù)小波變換可同時給出時間序列變化的位相和振幅兩方面的信息，有利于對問題的進一步分析。1.3復(fù)小波函數(shù)的實部和虛部位相差為n/2,能夠消除用實小波變換系數(shù)作為判據(jù)而產(chǎn)生的虛假振蕩，使分析結(jié)

10、果更為準(zhǔn)確。繪制小波系數(shù)圖、小波方差圖和主周期變化趨勢圖是關(guān)鍵當(dāng)選擇好合適的基小波函數(shù)后，下一步的關(guān)鍵就是如何通過小波變換獲得小波系數(shù)，然后利用相關(guān)軟件繪制小波系數(shù)圖、小波方差圖和主周期變化趨勢圖，進而根據(jù)上述三種圖形的變化識別徑流時間序列中存在的多時間尺度。具體步驟數(shù)據(jù)格式的轉(zhuǎn)化邊界效應(yīng)的消除或減小計算小波系數(shù)計算復(fù)小波系數(shù)的實部、模、模方、方差繪制小波系數(shù)實部、模、模方等值線圖繪制小波方差圖繪制主周期趨勢圖下面，我們以上題為例，結(jié)合軟件MatlabR2012a、suffer12.0、Excel、記事本等，詳細說明小波系數(shù)的計算和各圖形的繪制過程，并分別說明各圖在分析徑流多時間尺度變化特征中

11、的作用。數(shù)據(jù)格式的轉(zhuǎn)化和保存將存放在Excel表格里的徑流數(shù)據(jù)(以時間為序排為一列)轉(zhuǎn)化為MatlabR2012a識別的數(shù)據(jù)格式(.mat)并存盤。具體操作為：在MatlabR2012a界面下，單擊“File-ImportData”，出現(xiàn)文件選擇對話框“Import”后，找到需要轉(zhuǎn)化的數(shù)據(jù)文件(本例的文件名為runoff.xls)，單擊“打開”。等數(shù)據(jù)轉(zhuǎn)化完成后，單擊“Finish”，出現(xiàn)圖1顯示界面；然后雙擊圖1中的Runoff，彈出“ArrayEditor:runoff”對話框，選擇File文件夾下的SaveWorkspaceAs”單擊，出現(xiàn)圖2所示的“SavetoMAT-File:”窗口

12、，選擇存放路徑并填寫文件名(runoff.mat)，單擊“保存”并關(guān)閉“SavetoMAT-File窗口。a“IjJil-l-CffilLALlElKE-I.B-*nKiralaiq-mtKH-gSil屮i-b時Vi-CK1J13iEIKfd.ft*.lircmr!dl-rirtF,!S.ImH.圖1數(shù)據(jù)格式的轉(zhuǎn)化圖2數(shù)據(jù)的保存邊界效應(yīng)的消除或減小因為本例中的實測徑流數(shù)據(jù)為有限時間數(shù)據(jù)序列，在時間序列的兩端可能會產(chǎn)生“邊界效用”。為消除或減小序列開始點和結(jié)束點附近的邊界效應(yīng)，須對其兩端數(shù)據(jù)進行延伸。在進行完小波變換后，去掉兩端延伸數(shù)據(jù)的小變換系數(shù)，保留原數(shù)據(jù)序列時段內(nèi)的小波系數(shù)。本例中，我們利

13、用MatlabR2012a小波工具箱中的信號延伸(SignalExtension)功能，對徑流數(shù)據(jù)兩端進行對稱性延伸。具體方法為：在MatlabR2012a界面的“CommandWindow”中輸入小波工具箱調(diào)用命令“Wavemenu”，按Enter鍵彈“WaveletToolboxMainMenU(小波工具箱主菜單)界面(圖3)；然后單擊“SignalExtension”，打開SignalExtension/Truncation窗口，單擊“File”菜單下的“LoadSignal，選擇runoff.mat文件單擊“打開”，出現(xiàn)圖4信號延伸界面。MatlabR2012a的ExtensionMo

14、de菜單下包含了6種基本的延伸方式(Symmetric、Periodic、ZeroPadding、ContinuousSmoothandForSWT)和Directiontoextend菜單下的3種延伸模式(Both、LeftandRight)，在這里我們選擇對稱性兩端延伸進行計算。數(shù)據(jù)延伸的具體操作過程是：DesiredLength可以任意選，只要比原始信號長度大，建議在原始信號的基礎(chǔ)上加20(這樣左右對稱地延伸10個數(shù)據(jù))，這里選擇默認(rèn)的64；Dircetiontoextend下選擇“Both”；ExtensionMode下選擇“Symmetric”；單擊“Extend”按鈕進行對稱性兩端

15、延伸計算，然后單擊“File”菜單下的“SaveTranformedSignal”，將延伸后的數(shù)據(jù)結(jié)果存為erunoff.mat文件。從erunoff文件可知，系統(tǒng)自動將原時間序列數(shù)據(jù)向前對稱延伸12個單位，向后延伸13個單位。l3i|.k|aiui|uJl圖5小波變換菜單界面圖3小波工具箱主菜單圖4徑流時間序列的延伸計算小波系數(shù)選擇MatlabR2012a小波工具箱中的Morlet復(fù)小波函數(shù)對延伸后的徑流數(shù)據(jù)序列(erunoff.mat)進行小波變換，計算小波系數(shù)并存盤。小波工具箱主菜單界面見圖3,單擊“Wavelet1D”下的子菜單“ComplexContinuousWavelet1-D，

16、打開一維復(fù)連續(xù)小波界面，單擊“File”菜單下的“LoadSignal”按鈕，載入徑流時間序列erunoff.mat(圖5)。圖5的左側(cè)為信號顯示區(qū)域，右側(cè)區(qū)域給出了信號序列和復(fù)小波變換的有關(guān)信息和參數(shù)，主要包括數(shù)據(jù)長度(DataSize)、小波函數(shù)類型(Wavelet：cgau、shan、fbsp和cmor)、取樣周期(SamplingPeriod)、周期設(shè)置(ScaleSetting)和運行按鈕(Analyze)，以及顯示區(qū)域的相關(guān)顯示設(shè)置按鈕。本例中，我們選擇cmor(1-1.5)、取樣周期為1、最大尺度為32，單擊“Analyze”運行按鈕，計算小波系數(shù)。然后單擊“File”菜單下的“

17、SaveCoefficients”，保存小波系數(shù)為cerunoff.mat文件。計算Morlet復(fù)小波系數(shù)的實部、模、模方、方差見圖6。在MatlabR2012a界面下的Workspace中將cerunoff.mat文件導(dǎo)入,圖6小波系數(shù)導(dǎo)入到Matlab然后雙擊“coefs”打開，刪掉掉延伸數(shù)據(jù)的小波變換系數(shù)(本例中去掉前12列和后13列)，保存。接下來開始計算Morlet復(fù)小波系數(shù)的實部、模、模方、方差，具體操作為：在CommandWindows”中直接輸入函數(shù)“shibu=real(coefs);”，點擊“回車”鍵,計算實部；輸入函數(shù)“mo=abs(coefs);”，點擊“回車”鍵，計算

18、模；輸入函數(shù)“mofang=(mo).A2;”，點擊“回車”鍵，計算模方；輸入函數(shù)“fangcha=sum(abs(coefs).人2,2);”，點擊“回車”鍵，計算方差。見圖7。SelectdatatoplotNameValue-MinMaxwM2丸49doeplex已-OE31L3OS”,3.466-.5.464.32x1doubleu算49donbl日=u弓Toubl嚴(yán)fan日uh日0.83&54&3,8,EQfarigshitiu32x49doubleggggg2,341”,1S2-415.284.2.圖7計算出的實部、模、模方、方差成果注意：上面涉及到的數(shù)據(jù)保存，其格式均為.mat。

19、繪制小波系數(shù)實部、模、模方等值線圖實部、模、模方等值線圖的繪制方法一樣，這里僅以實部等值線圖為例。小波系數(shù)實部等值線圖的繪制首先，將小波系數(shù)實部數(shù)據(jù)復(fù)制到Excel中按照圖8格式排列，其中列A為時間，列B為尺度，列C為不同時間和尺度下所對應(yīng)的小波系數(shù)實部值。其次，將圖9數(shù)據(jù)轉(zhuǎn)化成Suffer12.0識別的數(shù)據(jù)格式。具體操作為：在Surfer12.0界面下，單擊“網(wǎng)格”菜單下的“數(shù)據(jù)”按鈕，在“打開”窗口選擇要打開的文件(小波系數(shù)實咅.Xls)，單擊“打開”后彈出“網(wǎng)格化數(shù)據(jù)”對話框(圖10)。它給出了多種不同的網(wǎng)格化方法、文件輸出路徑及網(wǎng)格線索幾何學(xué)等信息。這里我們選擇“克里格“網(wǎng)格方法”，單

20、擊“確定”，完成數(shù)據(jù)格式的轉(zhuǎn)化。理打0!5甘岸竺&!二申總魁直舉郡”埜1BCD11O-OO213222-O.04331030.CO3E34a-o.LM675-0-2063600.3336570.aB1.1i95g90.6219T100-071172A11-0.327271212-0.27-571313-0.燉4旳1-JC-.145-ta150.113B5坊-a.02?IT17-0.45T3913IB-G.6T33513-1-15052：-1.61392i-1-872722-2-039523-2-41672a-2_4dCl25-3.5-30735-2-793327-2-629129-3-3771

21、23-2-53-133：-2_MM3i-3-317732-3-132433lO-QOKiSa3-2-fl-D55a3530a01073355do.ouiss3750.1T730L/SkKTZ/ShHHt也丈洋狒很團叨fliAil：i稻式Q）工；L1因日直總歸I屁悒町圖境A1-右1日隨AILKecgnizedTyp礙(*.)ACCDBAccuseDitilise(*.iccdb)BLNGoldenSoftwareBlHiking(*.ilk)BNAAtlasEoundary(*.but)CSVCommaSeparatedVariables(*.csv)DAIDita(*.iatiDBFDatab

22、ase(*.ibfDXFAutoCADDrawingData你.dsf)LASLiDARBata(.las)MBBAccess97_2003Datalase(1=.mdb)SEG-P1DataEkchangeFormat保.seg)PlDataExchangeFormat傑.spl)SLKSylkSpreadsheet(*.slk)TXTTextD址玄(*.txt)XLSEclSprdzheet(*.kIzXLSXxcel2007Spreadsheet(.xIzk)XLSMxcel2007Spreadsheet(.xlzm)圖9Suffer12.0可以識別的數(shù)據(jù)格式列表jJa進狂列(1246.

23、過遢儡|此的1光生靈料人撲建汞舉霎鋼“整嫡十桶出網(wǎng)梧文件冋格議隸幾何學(xué)展小H方向|1.|-|.IH.律圖10小波系數(shù)實部數(shù)據(jù)格式轉(zhuǎn)化圖11Suffer8.0中的小系數(shù)實部等值線圖圖8小波系數(shù)實部數(shù)據(jù)格式最后，繪制小波系數(shù)實部等值線圖。在Surfer12.0界面下，單擊“地圖”菜單下的“等值線圖新建等值線圖”按鈕，彈出“打開網(wǎng)格”窗口后，選擇“小波系數(shù)實部grd”文件，單擊“打開”，完成等值線圖的繪制并保存（圖11）。小波系數(shù)實部等值線圖在多時間尺度分析中的作用30252015105LLLLHCVLL丿H1970197519801985199019952000年份圖12小系數(shù)實部等值線圖1970

24、197519801985199019952000年份u1970197519801985199019952000年份-U.5小波系數(shù)實部等值線圖能反映徑流序列不同時間尺度的周期變化及其在時間域中的分布，進而能判斷在不同時間尺度上，徑流的未來變化趨勢。為能比較清楚的說明小波系數(shù)實部等值線圖在徑流多時間尺度分析中的作用，我們利用Surfer12.0對其進一步處理和修飾，得到圖12顯示的小波系數(shù)實部等值線圖。其中，橫坐標(biāo)為時間（年份），縱坐標(biāo)為時間尺度，圖中的等值曲線為小波系數(shù)實部值。當(dāng)小波系數(shù)實部值為正時，代表徑流豐水期，在圖中我們用實線繪出，“H”表示正值中心；為負時，表示徑流枯水期，用虛線繪出，

25、“L”表示負值中心。由圖12可以清楚的看到徑流演化過程中存在的多時間尺度特征。總的來說，在流域徑流演變過程中存在著1832年，817年以及37年的3類尺度的周期變化規(guī)律。其中，在1832年尺度上出現(xiàn)了枯-豐交替的準(zhǔn)兩次震蕩；在817年時間尺度上存在準(zhǔn)5次震蕩。同時，還可以看出以上兩個尺度的周期變化在整個分析時段表現(xiàn)的非常穩(wěn)定，具有全域性；而310年尺度的周期變化，在1980s以后表現(xiàn)的較為穩(wěn)定。小波系數(shù)模和模方等值線圖的繪制參考5.1，繪制小波系數(shù)模和模方等值線圖（圖13、14）。圖13小波系數(shù)模等值線圖圖14小波系數(shù)模方等值線圖小波系數(shù)模等值線圖在多時間尺度分析中的作用Morlet小波系數(shù)的

26、模值是不同時間尺度變化周期所對應(yīng)的能量密度在時間域中分布的反映，系數(shù)模值愈大，表明其所對應(yīng)時段或尺度的周期性就愈強。從圖13可以看出，在流域徑流演化過程中，1832年時間尺度模值最大，說明該時間尺度周期變化最明顯，1822年時間尺度的周期變化次之，其他時間尺度的周期性變化較??；小波系數(shù)模方等值線圖在多時間尺度分析中的作用小波系數(shù)的模方相當(dāng)于小波能量譜，它可以分析出不同周期的震蕩能量。由圖14知，2532年時間尺度的能量最強、周期最顯著，但它的周期變化具有局部性（1980s前）；1015年時間尺度能量雖然較弱，但周期分布比較明顯，幾乎占據(jù)整個研究時域（19742004年）。繪制小波方差圖小波方差

27、圖的繪制在圖7的“fangcha”上右擊，選擇Graph”，在下拉菜單中選擇“plot”，即出小波方差圖，見圖15，在Matlab中可繼續(xù)美化。也可雙擊“fangcha”，將數(shù)據(jù)復(fù)制到其他軟件（如Excel）中，以小波方差為縱坐標(biāo)，時間尺度a為橫坐標(biāo)，繪制小波方差，如圖16。差方波小5101520253035時間尺度/a（d）4020002040oO86圖15Matlab繪制的小波方差圖圖16小波方差圖小波方差圖在多時間尺度分析中的作用小波方差圖能反映徑流時間序列的波動能量隨尺度a的分布情況?？捎脕泶_定徑流演化過程中存在的主周期。流域徑流的小波方差圖中（圖15）存在4個較為明顯的峰值，它們依次

28、對應(yīng)著28年、14年、8年和4年的時間尺度。其中，最大峰值對應(yīng)著28年的時間尺度，說明28年左右的周期震蕩最強，為流域年徑流變化的第一主周期；14年時間尺度對應(yīng)著第二峰值，為徑流變化的第二主周期，第三、第三峰值分別對應(yīng)著8年和4年的時間尺度，它們依次為流域徑流的第三和第四主周期。這說明上述4個周期的波動控制著流域徑流在整個時間域內(nèi)的變化特征。主周期趨勢圖的繪制及其在多時間尺度分析中的作用根據(jù)小波方差檢驗的結(jié)果，我們繪制出了控制流域徑流演變的第一和第二主周期小波系數(shù)圖（圖17）。從主周期趨勢圖中我們可以分析出在不同的時間尺度下，流域徑流存在的平均周期及豐-枯變化特征。圖16a顯示，在14年特征時

29、間尺度上，流域徑流變化的平均周期為9.5年左右，大約經(jīng)歷了4個豐-枯轉(zhuǎn)換期；而在28年特征時間尺度上（圖16b）,流域的平均變化周期為20年左右，大約2個周期的豐-枯變化。18642aaaa數(shù)系波小0-0.219-0.4-0.6-0.8rja14特征時間尺度o-98.IJ5891599121數(shù)系波小-1-3圖17大沽夾河流域年徑流變化的13年和28年特征時間尺度小波實部過程線參考文獻王文圣，丁晶，李耀清.2005.水文小波分析M.北京：化學(xué)工業(yè)出版社曹素華等.1998.實用醫(yī)學(xué)多因素統(tǒng)計方法M.上海：上海醫(yī)科大學(xué)出版社方開泰.1989.實用多元統(tǒng)計分析M.上海：華東師范大學(xué)出版社何清波，蘇炳華，錢亢.2002.醫(yī)學(xué)統(tǒng)計學(xué)及其軟件包M.上海：上?？茖W(xué)技術(shù)文獻出版社胡秉民.1987.微電腦在農(nóng)業(yè)科學(xué)中的應(yīng)用M.北京：科學(xué)出版社孫尚拱.1990.實用多元變量統(tǒng)