1、时间序列-小波分析时间序列(TimeSeries)是地学研究中经常遇到的问题。在时间序列研究中，时域和频域是常用的两种基本形式。其中，时域分析具有时间定位能力，但无法得到关于时间序列变化的更多信息；频域分析(如Fourier变换)虽具有准确的频率定位功能，但仅适合平稳时间序列分析。然而，地学中许多现象(如河川径流、地震波、暴雨、洪水等)随时间的变化往往受到多种因素的综合影响，大都属于非平稳序列，它们不但具有趋势性、周期性等特征，还存在随机性、突变性以及“多时间尺度”结构，具有多层次演变规律。对于这类非平稳时间序列的研究，通常需要某一频段对应的时间信息，或某一时段的频域信息。显然，时域分析和频域
2、分析对此均无能为力。20世纪80年代初，由Morlet提出的一种具有时-频多分辨功能的小波分析(WaveletAnalysis)为更好的研究时间序列问题提供了可能，它能清晰的揭示出隐藏在时间序列中的多种变化周期，充分反映系统在不同时间尺度中的变化趋势，并能对系统未来发展趋势进行定性估计。目前，小波分析理论已在信号处理、图像压缩、模式识别、数值分析和大气科学等众多的非线性科学领域内得到了广泛的应。在时间序列研究中，小波分析主要用于时间序列的消噪和滤波，信息量系数和分形维数的计算，突变点的监测和周期成分的识别以及多时间尺度的分析等。一、小波分析基本原理1.小波函数小波分析的基本思想是用一簇小波函数
3、系来表示或逼近某一信号或函数。因此，小波函数是小波分析的关键，它是指具有震荡性、能够迅速衰减到零的一类函数，即小波函数中(t)wL2(R)且满足：承(t)dt=0(1)式中，中(t)为基小波函数，它可通过尺度的伸缩和时间轴上的平移构成一簇函数系：/21t-b中a,b(t)=a中()其中，a,bWR,a=0(2)a式中，+a,b(t)为子小波；a为尺度因子，反映小波的周期长度；b为平移因子，反应时间上的平移。需要说明的是，选择合适的基小波函数是进行小波分析的前提。在实际应用研究中，应针对具体情况选择所需的基小波函数；同一信号或时间序列，若选择不同的基小波函数，所得的结果往往会有所差异，有时甚至差
4、异很大。目前，主要是通过对比不同小波分析处理信号时所得的结果与理论结果的误差来判定基小波函数的好坏，并由此选定该类研究所需的基小波函数。2.小波变换若Wa,b(t)是由(2)式给出的子小波，对于给定的能量有限信号f(t)WL2(R),其连续小波变换(ContinueWaveletTransform,简写为CWT为：-1/2/t-bWf(a,b)=|a1M产()dt(3)Ra式中，Wf(a,b)为小波变换系数；f(t)为一个信号或平方可积函数；a为伸缩尺度；b平移参数；xbxb中()为“)的复共轲函数。地学中观测到的时间序列数据大多是离散的，设函数f(kN),（k=1,2,N;At为取样间隔），
5、则式（3）的离散小波变换形式为:-1/2NkAt-b(4)Wf(a,b)=a.：tf(k.：t),-()kja由式（3）或（4）可知小波分析的基本原理，即通过增加或减小伸缩尺度a来得到信号的低频或高频信息，然后分析信号的概貌或细节，实现对信号不同时间尺度和空间局部特征的分析。实际研究中，最主要的就是要由小波变换方程得到小波系数，然后通过这些系数来分析时间序列的时频变化特征。3.小波方差将小波系数的平方值在b域上积分，就可得到小波方差，即二2Var（a）=JWf（a,b）db小波方差随尺度a的变化过程，称为小波方差图。由式（5）可知，它能反映信号波动的能量随尺度a的分布。因此，小波方差图可用来确
6、定信号中不同种尺度扰动的相对强度和存在的主要时间尺度，即主周期。二、小波分析实例-时间序列的多时间尺度分析（Multi-timescaleanalysis）例题河川径流是地理水文学研究中的一个重要变量，而多时间尺度是径流演化过程中存在的重要特征。所谓径流时间序列的多时间尺度是指：河川径流在演化过程中，并不存在真正意义上的变化周期，而是其变化周期随着研究尺度的不同而发生相应的变化，这种变化一般表现为小时间尺度的变化周期往往嵌套在大尺度的变化周期之中。也就是说，径流变化在时间域中存在多层次的时间尺度结构和局部变化特征。表1给出了某流域某水文观测站1966-2004年的实测径流数据。试运用小波分析理
7、论，借助Matlab6.5、suffer8.0和相关软件（Excel等），完成下述任务：计算小波系数；绘制小波系数图（实部、模和模方）、小波方差图和主周期变化趋势图，并分别说明各图在分析径流多时间尺度变化特征中的作用。表1某流域某水文观测站1966-2004年实测径流数据（X10稿）年份径流量年份径流量年份径流量年份径流量年份径流量19661.43819742.23519820.77419901.80619981.70919671.15119754.37419830.36719910.44919990.00019680.53619764.21919840.56219920.12020000.0
8、0019691.47019772.59019853.04019930.62720012.10419703.47619783.35019860.30419941.65820020.00919714.06819792.54019870.72819951.02520033.17719722.14719800.80719880.49219960.95520040.92119733.93119810.57319890.00719971.341分析1 .选择合适的基小波函数是前提在运用小波分析理论解决实际问题时，选择合适的基小波函数是前提。只有选择了适合具体问题的基小波函数，才能得到较为理想的结果。目前，
9、可选用的小波函数很多，如Mexicanhat小波、Haar小波、Morlet小波和Meyer小波等。在本例中，我们选用Morlet连续复小波变换来分析径流时间序列的多时间尺度特征。原因如下：1.1 径流演变过程中包含“多时间尺度”变化特征且这种变化是连续的，所以应采用连续小波变换来进行此项分析。1.2 实小波变换只能给出时间序列变化的振幅和正负，而复小波变换可同时给出时间序列变化的位相和振幅两方面的信息，有利于对问题的进一步分析。1.3 复小波函数的实部和虚部位相差为兀/2,能够消除用实小波变换系数作为判据而产生的虚假振荡，使分析结果更为准确。2 .绘制小波系数图、小波方差图和主周期变化趋势图
10、是关键当选择好合适的基小波函数后，下一步的关键就是如何通过小波变换获得小波系数，然后利用相关软件绘制小波系数图、小波方差图和主周期变化趋势图，进而根据上述三种图形的变化识别径流时间序列中存在的多时间尺度。具体步骤.数据格式的转化.边界效应的消除或减小.计算小波系数.计算复小波系数的实部.绘制小波系数实部等值线图.绘制小波系数模和模方等值线图.绘制小波方差图.绘制主周期趋势图下面，我们以上题为例，结合软件Matlab6.5、Suffer8.0和Excel,详细说明小波系数的计算和各图形的绘制过程，并分别说明各图在分析径流多时间尺度变化特征中的作用。.数据格式的转化和保存将存放在Excel表格里的
11、径流数据（以时间为序排为一列）转化为Matlab6.5识别的数据格式（.mat）并存盘。具体操作为：在Matlab6.5界面下，单击File-ImportData”,出现文件选择对话框Import”后，找到需要转化的数据文件（本例的文件名为runoff.xls）,单击“打开”。等数据转化完成后，单击“Finish:出现图1显示界面；然后双击图1中的Runoff,弹出ArrayEditor:runoff”对话框，选择File文件夹下的SaveWorkspaceAs”单击，出现图2所示白SSavetoMAT-File:窗口，选择存放路径并填写文件名（runoff.mat）,单击保存并关闭Savet
12、oMAT-File”窗口。图1数据格式的转化图2数据的保存边界效应的消除或减小因为本例中的实测径流数据为有限时间数据序列，在时间序列的两端可能会产生“边界效用”。为消除或减小序列开始点和结束点附近的边界效应，须对其两端数据进行延伸。在进行完小波变换后，去掉两端延伸数据的小变换系数，保留原数据序列时段内的小波系数。本例中，我们利用Matlab6.5小波工具箱中的信号延伸(SignalExtension)功能，对径流数据两端进行对称性延伸。具体方法为：在Matlab6.5界面的CommandWindoW中输入小波工具箱调用命令Wavemenu,按Enter键弹WaveletToolboxMainM
13、enu”(小波工具箱主菜单)界面(图3);然后单击“SignalExtension”,打开SignalExtension/Truncation窗口，单击File菜单下的LoadSignal,选择runoff.mat文件单击打开，出现图4信号延伸界面。Matlab6.5的ExtensionMode菜单下包含了6种基本的延伸方式(Symmetric、Periodic、ZeroPadding、Continuous、SmoothandForSWT和Directiontoextend菜单下的3种延伸模式(Both、LeftandRight),在这里我们选择对称性两端延伸进行计算。数据延伸的具体操作过程是
14、在ExtensionMode下选择Symmetric”,Dircetiontoextend下选择Both,单击Extend”按钮进行对称性两端延伸计算，然后单击File菜单下的SaveTranformedSignal,将延伸后的数据结果存为erunoff.mat文件。从erunoff文件可知，系统自动将原时间序列数据向前对称延伸12个单位，向后延伸13个单位。图3小波工具箱主菜单图4径流时间序列的延伸图5小波变换菜单界面计算小波系数选才MMatlab6.5小波工具箱中的Morlet复小波函数对延伸后的径流数据序列(erunoff.mat)进行小波变换，计算小波系数并存盘。小波工具箱主菜单界面
15、见图3,单击“Wavelet1-D”下的子菜单“ComplexContinuousWavelet1-D”,打开一维复连续小波界面，单击“File”菜单下的LoadSignal”按钮，载入径流时间序歹Uerunoff.mat(图5)。图5的左侧为信号显示区域，右侧区域给出了信号序列和复小波变换的有关信息和参数，主要包括数据长度(DataSize)、小波函数类型(Wavelet:cgau、shan、fbsp和cmor)、取样周期(SamplingPeriod)、周期设置(ScaleSetting)和运行按钮(Analyze),以及显示区域的相关显示设置按钮。本例中，我们选择cmor(1-1.5)、
16、取样周期为1、最大尺度为32,单击“Analyze”运行按钮，计算小波系数。然后单击File菜单下的SaveCoefficients”,保存小波系数为cerunoff.mat文件。.mat。注意：上面涉及到的数据保存，其格式均为计算Morlet复小波系数的实部将复小波系数转存到Excel表格，去掉两端延伸数据的小波系数，并计算小波系数实部。在Matlab6.5界面下的Workspace中将cerunoff.mat文件导入，然后双击打开，全部复制到Excel后去掉延伸数据的小波变换系数（本例中去掉前12列和后13歹U）,或只复制原时间序列的小波变换系数到Excel,最后使用Excel中的IMRE
17、AL数计算原时间序列的小波系数实部（图6）。当主快（11蝌疆工）吼图0电人口格式缶工具（D期居工）窗匚如晶耨M34XJ病=IIRBAL(B2)1H00上Qirt.$5建它一口。Qd$gBDQQasri2一口E款3中.另KERw一1-I(155155005105-0-M4i511I,ITTSVDi.05731-0g4T993a咽才词F,。14日W萼Q目榜D耻4田方口.11到10.1型斯口。3谶机F.nmLFZ4LF5D115*0.CdTlfl;D1553*0.z071121-.l2l.194&l-fl.11皤叩瘴如器O.SE：3S-(.046)4C.5453711371018272-0.-D.O
18、CEO-3IJD041C,工口Z53oJ1JIB44D213Z5LI0LOF0SS9-O.-0.M051V.0“打JJ4571中口175+6071第既T7.；1m*”.2-0.miaDGlMi-fl,14-tZVO.OTa0*0555*04b34r-(JJ.LLE111.。绮Lq0.52157-0.IEU.CifiItJUU51SJF3/如4千口.？打熊：-V.41r-02Z144bUl70.MM1中rlHKEL2d10140_001S-971+-).C54W-I旧TAHJQ归用7433&0-A13270.55601T-D)14i4-011-0.55727+0.2015?二,002123-2.
19、359：LTrajnberyaL1j4rs4t+上Lu700.20M07*-3之m_LIERT5产Q,13024345-0.OLfS14u.i.4fe*4y*o.j.13eln*o.aloOieeM-9mb*74K+a.3-D.1835-1：145B33-CLf_L11E0.11t9)Ll&ko,027011.EIE1?-d4UWF.tjTSSl求其雷我颤的曳泡411S*0,-0.M66?-，打叫-L!5K62*0-0.1H36-L12023I.J425-a.aZ0611-G3ii8r-a.sIL-U鼻湿5U22JZ2ar1-1.1WW.加卜1国瓢F,髓21-1.gff湖11双123-2.05
20、0.IHflL嚣-2-41S7+C.工褥i4a44tti4.ai-23-254C6*a.3$kLL汨-2SMKLS7FMD抗弱-21smQ.5(17411.於-z5M3+A.W3J5L九-2定U/L31之81T7+H,47411l能-.1。烈TQ,51稣Gilb计算使累3.3ZC13E3碎京I度清|iklffiXALiBt-.5储-QEDA口英-2.K畲IF518S1FTJftiT-l4MaA1-2.64-5口8。21-EM，T耳74&Q4i-2.匚R2-IJ3H3L51T69sami-2.321j-D7iY.猊安FL65谢一”55-LEW-t116-1.45361-m-13r“2307-L
21、23W1-1,2li一IW】-/-ii-i7IMM.96021i-2Q&科一&MSi-l.酢3h.Q滤葡-1.511-5.0863-1si4：-1!：-1.6GS-S.imi-I.T6kl.用&*3i7门Jd-1匕-1,QCd-l.MOJs-1,n除L才调i17-0.fl541-1-2631.5珈7M2&S-2.15K-2-2-d.f43,t3-S+50&I333力？437彳钝行451比7耨酎望南-.C655-ELUZ31一口声MLN.卜LM51ISL鼬|n.99而-j?-J,9851-0旬M3-E.731m-l吃亍Ij：-Q.河丁4E13TI3oiM-oa.1I302-01-513-161
22、11.6SDB-DUEL4-1.mi-02(I3-l15441-12.JIMOrEIB*59-l221M79-21r蚪Al2llis-ei5frlC-l2.1Q呢1洸7川再7J.:147斤】30l3nzXJa上学life-，取扪*i5yEI”，4)图6复小波系数及实部计算示意图Excel中IMREAL函数的调用方法为：单击“插入”菜单下的“函数（F）”按钮，弹出图7所示的“插入函数”窗口，在“搜索函数（S）:框中输入：“IMREAL后单击“转到”，再单击“确定”，出现函数参数窗口（图8）。在“Inumber”一栏的空白处输入所要计算的数据（图6）,单击“确定”即可得到小波系数实部值。图7IM
23、EAL函数调用图8IMEAL函数参数需要说明的是，从cerunoff.mat文件中转到Excel里的复小波系数，在其实部和虚部中间包含许多“空格”，在计算之前需要先将其去掉。借助Suffer8.0,绘制小波系数实部等值线图波系数实部等值线图的绘制首先，将小波系数实部数据按照图9格式排列，其中列A为时间，列B为尺度，列C为不同时间和尺度下所对应的小波系数实部值。图9小波系数实部数据格式其次，将图9数据转化成Suffer8.0识别的数据格式。具体操作为：在Suffer8.0界面下，单击“网格”菜单下的“数据”按钮，在“打开”窗口选择要打开的文件（小波系数实部.xls）,单击“打开”后弹出“网格化数
24、据”对话框（图10）。它给出了多种不同的网格化方法、文件输出路径及网格线索几何学等信息。这里我们选择“克里格“网格方法”，单击“确定”，完成数据格式的转化。图10小波系数实部数据格式转化图11Suffer8.0中的小系数实部等值线图一1.二二3图12小系数实部等值线图13、14）。最后，绘制小波系数实部等值线图。在Suffer8.0界面下，单击“地图”菜单下的“等值线图-新建等值线图”按钮，弹出“打开网格”窗口后，选择“小波系数实部.grd”文件，单击“打开”，完成等值线图的绘制并存盘（图11）。波系数实部等值线图在多时间尺度分析中的作用小波系数实部等值线图能反映径流序列不同时间尺度的周期变化
25、及其在时间域中的分布，进而能判断在不同时间尺度上，径流的未来变化趋势。为能比较清楚的说明小波系数实部等值线图在径流多时间尺度分析中的作用，我们利用Suffer8.0对其进一步处理和修饰，得到图12显示的小波系数实部等值线图。其中，横坐标为时间（年份），纵坐标为时间尺度，图中的等值曲线为小波系数实部值。当小波系数实部值为正时，代表径流丰水期，在图中我们用实线绘出，H表示正值中心；为负时，表示径流枯水期，用虚线绘出，“L”表示负值中心。由图12可以清楚的看到径流演化过程中存在的多时间尺度特征。总的来说，在流域径流演变过程中存在着1832年，817年以及37年的3类尺度的周期变化规律。其中，在183
26、2年尺度上出现了枯-丰交替的准两次震荡；在817年时间尺度上存在准5次震荡。同时，还可以看出以上两个尺度的周期变化在整个分析时段表现的非常稳定，具有全域性；而310年尺度的周期变化，在1980s以后表现的较为稳定。绘制小波系数模和模方等值线图波系数模和模方等值线图的绘制参考4、5两步，绘制小波系数模和模方等值线图（图说明：在Excel中，复数模的计算函数为IMAB三。3025201510图13小波系数模等值线图302520151051970197519801985199019952000年伤图14小波系数模方等值线图图15小波方差图波系数模等值线图在多时间尺度分析中的作用Morlet小波系数的
27、模值是不同时间尺度变化周期所对应的能量密度在时间域中分布的反映，系数模值愈大，表明其所对应时段或尺度的周期性就愈强。从图13可以看出，在流域径流演化过程中，1832年时间尺度模值最大，说明该时间尺度周期变化最明显，1822年时间尺度的周期变化次之，其他时间尺度的周期性变化较小；波系数模方等值线图在多时间尺度分析中的作用小波系数的模方相当于小波能量谱，它可以分析出不同周期的震荡能量。由图14知，2532年时间尺度的能量最强、周期最显著，但它的周期变化具有局部性（1980s前）；1015年时间尺度能量虽然较弱，但周期分布比较明显，几乎占据整个研究时域（19742004年）。绘制小波方差图波方差图的
28、绘制将不同时间尺度下的小波系数代入式（5）可得径流变化的小波方差，以小波方差为纵坐标，时间尺度a为横坐标，可绘制小波方差图（图15）。波方差图在多时间尺度分析中的作用小波方差图能反映径流时间序列的波动能量随尺度a的分布情况。可用来确定径流演化过程中存在的主周期。流域径流的小波方差图中（图15）存在4个较为明显的峰值，它们依次对应着28年、14年、8年和4年的时间尺度。其中，最大峰值对应着28年的时间尺度，说明28年左右的周期震荡最强，为流域年径流变化的第一主周期；14年时间尺度对应着第二峰值，为径流变化的第二主周期，第三、第三峰值分别对应着8年和4年的时间尺度，它们依次为流域径流的第三和第四主
29、周期。这说明上述4个周期的波动控制着流域径流在整个时间域内的变化特征。主周期趋势图的绘制及其在多时间尺度分析中的作用根据小波方差检验的结果，我们绘制出了控制流域径流演变的第一和第二主周期小波系数图（图16）。从主周期趋势图中我们可以分析出在不同的时间尺度下，流域径流存在的平均周期及丰-枯变化特征。图16a显示，在14年特征时间尺度上，流域径流变化的平均周期为9.5年左右，大约经历了4个丰-枯转换20年左右，大约2个周期的丰-枯变化。期；而在28年特征时间尺度上（图16b）,流域的平均变化周期为a14特征时间尺度864220000-0数系波小4680.0-3b28特征时间尺度图16大沽夹河流域年
30、径流变化的13年和28年特征时间尺度小波实部过程线参考文献王文圣，丁晶，李耀清.2005.水文小波分析M.北京：化学工业出版社曹素华等.1998.实用医学多因素统计方法M.上海：上海医科大学出版社方开泰.1989.实用多元统计分析M.上海：华东师范大学出版社何清波，苏炳华，钱亢.2002.医学统计学及其软件包M.上海：上海科学技术文献出版社胡秉民.1987.微电脑在农业科学中的应用M.北京：科学出版社孙尚拱.1990.实用多元变量统计方法与计算程序M.北京：北京医科大学、中国协和医科大学联合出版社唐守正.1986.多元统计分析方法M。北京：中国林业出版社王学仁.1982.地址数据的多变量统计分
31、析.北京：科学出版社徐振邦，金淳浩，娄元仁.1986.22距离系数和邛2距离系数尺度在聚类分析中的应用M.赵旭东等主编，中国数学地质（1）.北京：地质出版社於崇文.1978.数学地质的方法与应用M.北京：冶金工业出版社AndersonT.W.1967.Introductiontomultivariatestatisticalanalysis,2ndM.NewYork:WileyGauchH.G.J.1982.MultivariateanalysisincommunityecologyM.Britain:CambridgeUniversityPressHorelA.E.,Wennard.R.W.
32、andBaldwinK.F.1975.regression:somesimulations.CommunicationsinStatisticsJ,4:105123练习试运用小波分析理论，分析某市年平均降水过程中存在的多时间尺度变化特征。表2某巾1957-2004年实测年均降水量（mm年份降水量年份降水量年份降水量年份降水量1957320.01969324.81981506.01993384.41958481.21970412.31982282.11994503.91959522.61971366.51983508.61995406.71960339.31972262.41984523.91996465.71961719.91973521.91985518.91997345.31962373.51974351.71986320.11998454.91963332.91975398.41987340.01999327.91964741.21976320.21988478.52000406.21965454.31977445.41989402.42001404.71966604.31978534.81990552.42002401.91967451.91979509.71991313.92003605.11968424.31980395.51992591.02004385.4
免费区 
