小波分析-經(jīng)典解讀

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

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

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

4、時(shí)甚至差異很大。目前，主要是通過(guò)對(duì)比不同小波分析處理信號(hào)時(shí)所得的結(jié)果與理論結(jié)果的誤差來(lái)判定基小波函數(shù)的好壞，并由此選定該類研究所需的基小波函數(shù)。2. 小波變換若屮(t)是由(2)式給出的子小波，對(duì)于給定的能量有限信號(hào)f(t)GL2(R)，其連續(xù)小波變換a,b(ContinueWaveletTransform，簡(jiǎn)寫為CWT)為：tbW(a,b)=|a|-1/2Jf(t(b)dtraxb式中，W(a,b)為小波變換系數(shù)；(t)為一個(gè)信號(hào)或平方可積函數(shù)a為伸縮尺度；平移參數(shù)W()faxb為屮()的復(fù)共軛函數(shù)。地學(xué)中觀測(cè)到的時(shí)間序列數(shù)據(jù)大多是離散的設(shè)函數(shù)f(kAt)(k=1,2,.,N;At為取樣間隔

5、)，則式(3)的離散小波變換形式為：W(a,b)二|a|-1/2Atff(kAt”(kAt-bfa由式(3)或(4)可知小波分析的基本原理，即通過(guò)增加或減小伸縮尺度a來(lái)得到信號(hào)的低頻或高頻信息，然后分析信號(hào)的概貌或細(xì)節(jié)，實(shí)現(xiàn)對(duì)信號(hào)不同時(shí)間尺度和空間局部特征的分析。實(shí)際研究中，最主要的就是要由小波變換方程得到小波系數(shù)，然后通過(guò)這些系數(shù)來(lái)分析時(shí)間序列的時(shí)頻變化特征。3. 小波方差將小波系數(shù)的平方值在b域上積分，就可得到小波方差，即(5)小波方差隨尺度a的變化過(guò)程，稱為小波方差圖。由式(5)可知，它能反映信號(hào)波動(dòng)的能量隨尺度a的分布。因此，小波方差圖可用來(lái)確定信號(hào)中不同種尺度擾動(dòng)的相對(duì)強(qiáng)度和存在的主

6、要時(shí)間尺度，即主周期。二、小波分析實(shí)例-時(shí)間序列的多時(shí)間尺度分析（Multi-timescaleanalysis）例題河川徑流是地理水文學(xué)研究中的一個(gè)重要變量，而多時(shí)間尺度是徑流演化過(guò)程中存在的重要特征。所謂徑流時(shí)間序列的多時(shí)間尺度是指：河川徑流在演化過(guò)程中，并不存在真正意義上的變化周期，而是其變化周期隨著研究尺度的不同而發(fā)生相應(yīng)的變化，這種變化一般表現(xiàn)為小時(shí)間尺度的變化周期往往嵌套在大尺度的變化周期之中。也就是說(shuō)，徑流變化在時(shí)間域中存在多層次的時(shí)間尺度結(jié)構(gòu)和局部變化特征。表1給出了某流域某水文觀測(cè)站1966-2004年的實(shí)測(cè)徑流數(shù)據(jù)。試運(yùn)用小波分析理論借助Matlab6.5、suffer8.

7、0和相關(guān)軟件（Excel等），完成下述任務(wù)：計(jì)算小波系數(shù)；繪制小波系數(shù)圖（實(shí)部、模和模方）、小波方差圖和主周期變化趨勢(shì)圖，并分別說(shuō)明各圖在分析徑流多時(shí)間尺度變化特征中的作用。表1某流域某水文觀測(cè)站1966-2004年實(shí)測(cè)徑流數(shù)據(jù)（Q08m3）年份徑流量年份徑流量年份徑流量年份徑流量年份徑流量19661.43819742.23519820.77419901.80619981.70919671.15119754.37419830.36719910.44919990.00019680.53619764.21919840.56219920.12020000.00019691.47019772.5901

8、9853.04019930.62720012.10419703.47619783.35019860.30419941.65820020.00919714.06819792.54019870.72819951.02520033.17719722.14719800.80719880.49219960.95520040.92119733.93119810.57319890.00719971.341分析1. 選擇合適的基小波函數(shù)是前提在運(yùn)用小波分析理論解決實(shí)際問(wèn)題時(shí)，選擇合適的基小波函數(shù)是前提。只有選擇了適合具體問(wèn)題的基小波函數(shù)，才能得到較為理想的結(jié)果。目前，可選用的小波函數(shù)很多，如Mexicanha

9、t小波、Haar小波、Morlet小波和Meyer小波等。在本例中，我們選用Morlet連續(xù)復(fù)小波變換來(lái)分析徑流時(shí)間序列的多時(shí)間尺度特征。原因如下：1.1 徑流演變過(guò)程中包含“多時(shí)間尺度”變化特征且這種變化是連續(xù)的，所以應(yīng)采用連續(xù)小波變換來(lái)進(jìn)行此項(xiàng)分析。1.2 實(shí)小波變換只能給出時(shí)間序列變化的振幅和正負(fù)，而復(fù)小波變換可同時(shí)給出時(shí)間序列變化的位相和振幅兩方面的信息，有利于對(duì)問(wèn)題的進(jìn)一步分析。1.3復(fù)小波函數(shù)的實(shí)部和虛部位相差為n/2，能夠消除用實(shí)小波變換系數(shù)作為判據(jù)而產(chǎn)生的虛假振蕩，使分析結(jié)果更為準(zhǔn)確。2. 繪制小波系數(shù)圖、小波方差圖和主周期變化趨勢(shì)圖是關(guān)鍵當(dāng)選擇好合適的基小波函數(shù)后，下一步的關(guān)

10、鍵就是如何通過(guò)小波變換獲得小波系數(shù)，然后利用相關(guān)軟件繪制小波系數(shù)圖、小波方差圖和主周期變化趨勢(shì)圖，進(jìn)而根據(jù)上述三種圖形的變化識(shí)別徑流時(shí)間序列中存在的多時(shí)間尺度。具體步驟1. 數(shù)據(jù)格式的轉(zhuǎn)化2. 邊界效應(yīng)的消除或減小3. 計(jì)算小波系數(shù)4. 計(jì)算復(fù)小波系數(shù)的實(shí)部5. 繪制小波系數(shù)實(shí)部等值線圖6. 繪制小波系數(shù)模和模方等值線圖7. 繪制小波方差圖8. 繪制主周期趨勢(shì)圖下面，我們以上題為例，結(jié)合軟件Matlab6.5、Suffer8.0和Excel，詳細(xì)說(shuō)明小波系數(shù)的計(jì)算和各圖形的繪制過(guò)程，并分別說(shuō)明各圖在分析徑流多時(shí)間尺度變化特征中的作用。1. 數(shù)據(jù)格式的轉(zhuǎn)化和保存將存放在Excel表格里的徑流數(shù)據(jù)

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

12、ile窗口。圖1數(shù)據(jù)格式的轉(zhuǎn)化圖2數(shù)據(jù)的保存2. 邊界效應(yīng)的消除或減小因?yàn)楸纠械膶?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ù)。本例中，我們利用Matlab6.5小波工具箱中的信號(hào)延伸(SignalExtension)功能，對(duì)徑流數(shù)據(jù)兩端進(jìn)行對(duì)稱性延伸。具體方法為：在Matlab6.5界面的“CommandWindow”中輸入小波工具箱調(diào)用命令Wavemenu,按Enter鍵彈WaveletToolboxMainMen

13、U(小波工具箱主菜單)界面(圖3)撚后單擊“SignalExtension,打開SignalExtension/Truncation窗口，單擊“File菜單下的“LoadSignal,選擇runoff.mat文件單擊“打開，出現(xiàn)圖4信號(hào)延伸界面oMatlab6.5的ExtensionMode菜單下包含了6種基本的延伸方式(Symmetric、Periodic、ZeroPadding、Continuous、SmoothandForSWT)和Directiontoextend菜單下的3種延伸模式(Both、LeftandRight)，在這里我們選擇對(duì)稱性兩端延伸進(jìn)行計(jì)算。數(shù)據(jù)延伸的具體操作過(guò)程是：

14、在ExtensionMode下選擇“Symmetric,Dircetiontoextend下選擇“Both，單擊“Extend按鈕進(jìn)行對(duì)稱性兩端延伸計(jì)算然后單擊“File菜單下的“SaveTranformedSigna將延伸后的數(shù)據(jù)結(jié)果存為erunoff.mat文件。從erunoff文件可知，系統(tǒng)自動(dòng)將原時(shí)間序列數(shù)據(jù)向前對(duì)稱延伸12個(gè)單位，向后延伸13個(gè)單位。圖3小波工具箱主菜單圖4徑流時(shí)間序列的延伸3. 計(jì)算小波系數(shù)圖5小波變換菜單界面選擇Matlab6.5小波工具箱中的Morlet復(fù)小波函數(shù)對(duì)延伸后的徑流數(shù)據(jù)序列(erunoff.mat)進(jìn)行小波變換，計(jì)算小波系數(shù)并存盤。小波工具箱主菜單界

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

16、取樣周期為1、最大尺度為32，單擊“Analyze運(yùn)行按鈕，計(jì)算小波系數(shù)。然后單擊“File菜單下的“SaveCoefficients，保存小波系數(shù)為cerunoff.mat文件。注意：上面涉及到的數(shù)據(jù)保存，其格式均為.mat。4. 計(jì)算Morlet復(fù)小波系數(shù)的實(shí)部將復(fù)小波系數(shù)轉(zhuǎn)存到Excel表格，去掉兩端延伸數(shù)據(jù)的小波系數(shù)，并計(jì)算小波系數(shù)實(shí)部。在Matlab6.5界面下的Workspace中將cerunoff.mat文件導(dǎo)入，然后雙擊打開，全部復(fù)制到Excel后去掉延伸數(shù)據(jù)的小波變換系數(shù)(本例中去掉前12列和后13列)，或只復(fù)制原時(shí)間序列的小波變換系數(shù)到Excel，最后使用Excel中的IM

17、REAL函數(shù)計(jì)算原時(shí)間序列的小波系數(shù)實(shí)部(圖6)。FlliCTOiEOfiernxaiiEf.tIg璉人需戛肖主件錦巻視團(tuán)迫插入格式加iMd）議據(jù)皿爾口on幫助on柑日日也召丨日丨眞丨曜丨刁丨皤鶴富SIW*辰-IIREAL(B2)CJE匚KI1避小頑系豪19莊L9&r19&91別L970197119721973Z1a0皿凸曠L3552e-oa.00305630.003S71-D.004780.0.011613-0.02TF50.a(ljO4EE-D.08733d1018272-0.-ClOOBi：iEl00041-a.2Z-1104鵡尸.3BZ6e-a-0.055T9B005100100643

18、6Z+CI11810.13554-Q.05363L-il.il3ZOL-nZ4L76-0Z3735+0.0.053540.DE15Z50.(430.OOCfeaaT-2.2F9&Q-0.010T：?3-D.O?Q2S&i-LL022&42+0.Q5$5啦dj.0M9WJ-Q.1025&-a.flSUtri-D.1957&0.0F0553-0.-0.040510.DDFt2?Q.1：5彳-0.1CES74：!.59G3E-OOQ.(14139-D.1159&iD.15415+0.-0.ise2+a.J94&L-a.1J3G-D.35a&2id：i.45719*0.1-0.54G07(i.2G99

19、7-a.:&5-0.2363-a.00J4L51a.iTTsrd.Eosr-jiD.13350.3OT3210.3S38-a.00604a.13rD.3435710.e閃？尹a2-0.2&7110.D63341-a.CL33E63-0.D0SZ3L31Iq.155TD.35Z161-030354+0.Z303310.41ZT0.2Z54104gzyp.OT|q.D40555D4&g47|-Q401+0.a|0.l&9714i-D.65193-0D.131T1-D.4&31*312024232&o.94id4+a.029Ttr1.1095-(1.i：i：!5242l0.6E13T-Q.OOfDS

20、口.D55SZ-0.3272740.2(iL59-0.27457+0.1T994-0.02483-0.0L.T3口.14E43-D.33ZBLo.uses-a.sa&psi0.027aa6i：i.32B2-Q.23T53L-D.7S3-fl.22222T.1導(dǎo)旳-Q.航-1.6193-0.030DFB-1.STffli-Q.0139EL-2.DS95-K).324-2LL-2.dlG7-K.2：?6B7l-巴44口(L:J2$B3i二瓷5406-i-Q.-z7S23+03524fiL-2.G29HQ.593D3l-2.37714a.50T4Ll-z.陽(yáng)4滬IL563Z5L-2SB440.5】0

21、2Fl91-2.317740.4T4DLl32-2.162M-H.513L61卅迪齊歎比鋼1|=IMFEA1（E2；IrinKEKLInimbar返回復(fù)數(shù)的實(shí)鄙Iznuihar求其實(shí)累軸壊匏.計(jì)算結(jié)果=0.CC213232.2.5$；:d：i53；86i2JSfiT-fl49SSai-2J24-4-12.G404-D.50&92i2.6T6-D.24B04i2.6T?-ri3Q31Si巳能陰T-E1155-1-2.2007-1-2.3191-1-EE535-12.525i-flami-22271-12.3217-D?L5937i2.226-D.LB529i71993-0-2.0B34-0X宜=

22、TO.LCei323-2.355CU2T433614&曲i-1.5141.91鈿2-404i.21323職消-1.Q265i-1.511-2.P61362fli-1.7201-1.5132i1&42J1499i-1.7074-1.5S2S-1.76?2.355395021i623331-1.Q-l.北I：打-1.71621.i1973-0.io.4ea-l-D.4aa&9-aPD.閃0.BD1胡FD14”-Q,5-d.a4ss4e-d.dfcj21-d7U*O.20.20075+-D.2O31&-0|079*O.j0.019694-0.349&-a42J1-0.3-Q.1885j-D.45B3

23、-q15+0.2-.3402-0.S3T99-n41120.1-0.29K7-D.270Ei-ai5556-a|-0.11436-D.D202304a&-0.00.206110.SBIStHo706-0.5Q.57594-40.E24&2-0n7：3-o.$m：汨is-1.dls：!-oaMl-1.21.0976-01.9902-01638-1.51.3161-1j.6E0&-01：344-l.q1.2755-1J.Brfi2-01202-2.111.434-1.J.9S91-020.64$i：iE45051d：i=;=;：B3-2.45L0.7fiflTfc-25抽“0099434-2?t0

24、.fl5!：3r-2.4519i-0.059-2.9G0.TTI663=2.21.2935-11.3441-1224021E；帕弓-12Cl刃2-2】?e9b-t29112-21.5G1G-L21.OOiJ-E.3OO4i7165752.d.32136-1.297E-11I.QB為-E.LlEiF35J1LE.43釘時(shí)】.汕:LdFl1.站51-2.L512i-D.299242.LI2341510.fl&45-1754149-2.0.2211甲0.-0J55754-2.d.0956-20.6915311.2G50-L.992511.2mi-L.B351i譏J3T319EEL568U6JiriM

25、TO圖6復(fù)小波系數(shù)及實(shí)部計(jì)算示意圖Excel中IMREAL函數(shù)的調(diào)用方法為：?jiǎn)螕簟安迦搿辈藛蜗碌摹昂瘮?shù)（F）”按鈕，彈出圖7所示的“插入函數(shù)”窗口，在搜索函數(shù)（S）:”框中輸入:“IMREAL”后單擊轉(zhuǎn)到”，再單擊“確定”，出現(xiàn)函數(shù)參數(shù)窗口（圖8）。在“Inumber”一欄的空白處輸入所要計(jì)算的數(shù)據(jù)（圖6）,單擊“確定”即可得到小波系數(shù)實(shí)部值。圖7IMEAL函數(shù)調(diào)用圖8IMEAL函數(shù)參數(shù)需要說(shuō)明的是，從cerunoff.mat文件中轉(zhuǎn)到ExceI里的復(fù)小波系數(shù)，在其實(shí)部和虛部中間包含許多“空格”，在計(jì)算之前需要先將其去掉。5. 借助Suffer8.0，繪制小波系數(shù)實(shí)部等值線圖5.1小波系數(shù)實(shí)部等值線圖的繪制首先，將小波系數(shù)實(shí)部數(shù)據(jù)按照?qǐng)D9格式排列，其中列A為時(shí)間，列B為尺度，列C為不同時(shí)間和尺度下所對(duì)應(yīng)的小波系數(shù)實(shí)部值。也文洋電I:嘶童）機(jī)曲世）rtA1）觀冠|0工