(完整word版)时间序列的小波分析及等值线图、小波方差制作.docx

时间序列的小波分析时间序列（ Time Series）是地学讨论中常常遇到的问题；在时间序列讨论中，时域和频域是常用的两种基本形式； 其中，时域分析具有时间定位才能，但无法得到关于时间序列变化的更多信息；频域分析 （如Fourier 变换）虽具有精确的频率定位功能，但仅适合平稳时间序列分析；然而，地学中很多现象（如河川径流、地震波、暴雨、洪水等）随时间的变化往往受到多种因素的综合影响，大都属于非平稳序列，它们不但具有趋势性、周期性等特点， 仍存在随机性、 突变性以及 “多时间尺度” 结构， 具有多层次演化规律；对于这类非平稳时间序列的讨论，通常需要某一频段对应的时间信息，或某一时段的频域信息；明显，时域分析和频域分析对此均无能为力；20 世纪 80 岁月初，由 Morlet 提出的一种具有时 -频多辨论功能的小波分析（Wavelet Analysis ）为更好的讨论时间序列问题供应了可能，它能清晰的揭示出隐匿在时间序列中的多种变化周期，充分反映系统在不同时间尺度中的变化趋势，并能对系统将来进展趋势进行定性估量；目前，小波分析理论已在信号处理、图像压缩、模式识别、数值分析和大气科学等众多的非线性科学领域内得到了广泛的应；在时间序列讨论中，小波分析主要用于时间序列的消噪和滤波，信息量系数和分形维数的运算，突变点的监测和周期成分的识别以及多时间尺度的分析等；一、小波分析基本原理1.小波函数小波分析的基本思想是用一簇小波函数系来表示或靠近某一信号或函数；因此，小波函数是小波分析L2 ( R)( t)的关键，它是指具有震荡性、能够快速衰减到零的一类函数，即小波函数且满意：( t)dt0（ 1）式中， (t) 为基小波函数，它可通过尺度的伸缩和时间轴上的平移构成一簇函数系：( t b ) a1/ 2(t)a其中， a,bR, a0（ 2）a ,ba, b (t ) 为子小波； a 为尺度因子，反映小波的周期长度；式中，b 为平移因子，反应时间上的平移；需要说明的是，挑选合适的基小波函数是进行小波分析的前提；在实际应用讨论中，应针对具体情形挑选所需的基小波函数；同一信号或时间序列，如挑选不同的基小波函数，所得的结果往往会有所差异，有时甚至差异很大；目前，主要是通过对比不同小波分析处理信号时所得的结果与理论结果的误差来判定基小波函数的好坏，并由此选定该类讨论所需的基小波函数；2.小波变换L2 ( R) ，其连续小波变换如 a ,b ( t)f (t)是由（ 2）式给出的子小波，对于给定的能量有限信号（ Continue Wavelet Transform ，简写为 CWT ）为：tb-1/ 2Wf ( a, b )af(t)()dt（ 3）Ra( xb )式中，W (a, b) 为小波变换系数；f(t) 为一个信号或平方可积函数；a 为伸缩尺度； b 平移参数；fa( x b ) 的复共轭函数；a,N; t为设函数 f (kt ) ，（ k=1,2,地学中观测到的时间序列数据大多是离散的，为取样间隔） ，就式（ 3）的离散小波变换形式为：N( k t - b) a-1/ 2aW (a, b)tf(kt)（ 4）fk1由式（ 3）或（ 4）可知小波分析的基本原理，即通过增加或减小伸缩尺度信息，然后分析信号的概貌或细节，实现对信号不同时间尺度和空间局部特点的分析；a 来得到信号的低频或高频实际讨论中，最主要的就是要由小波变换方程得到小波系数，然后通过这些系数来分析时间序列的时频变化特点；3.小波方差将小波系数的平方值在 b 域上积分，就可得到小波方差，即2Wf (a,b) db(5) 可知，它能反映信号波动的能量随尺度Var(a)小波方差随尺度 a 的变化过程，称为小波方差图；由式(5)a 的分布；因此，小波方差图可用来确定信号中不同种尺度扰动的相对强度和存在的主要时间尺度，即主周期；二、小波分析实例 - 时间序列的多时间尺度分析例题(Multi-time scale analysis)河川径流是地理水文学讨论中的一个重要变量，而多时间尺度是径流演化过程中存在的重要特点；所谓径流时间序列的多时间尺度是指：河川径流在演化过程中，并不存在真正意义上的变化周期，而是其变化周期随着讨论尺度的不同而发生相应的变化，这种变化一般表现为小时间尺度的变化周期往往嵌套在大尺度的变化周期之中；也就是说，径流变化在时间域中存在多层次的时间尺度结构和局部变化特点；表 1 给出了某流域某水文观测站R2021a、suffer 12.0 和其他相关软件（1966-2004 年的实测径流数据；试运用小波分析理论，借助Matlab2）绘Excel、记事本等），完成下述任务：（1）运算小波系数；（制小波系数图（实部、模和模方）、小波方差图和主周期变化趋势图，并分别说明各图在分析径流多时间尺度变化特点中的作用；8 3表 1 某流域某水文观测站 1966-2004 年实测径流数据（ 10 m）年份19661967196819691970197119721973径流量1.4381.1510.5361.4703.4764.0682.1473.931年份19741975197619771978197919801981径流量2.2354.3744.2192.5903.3502.5400.8070.573年份19821983198419851986198719881989径流量0.7740.3670.5623.0400.3040.7280.4920.007年份19901991199219931994199519961997径流量1.8060.4490.1200.6271.6581.0250.9551.341年份1998199920002001200220032004径流量1.7090.0000.0002.1040.0093.1770.921分析1. 挑选合适的基小波函数是前提在运用小波分析理论解决实际问题时，挑选合适的基小波函数是前提；只有挑选了适合具体问题的基Mexican hat 小波、 Haar 小波、小波函数，才能得到较为抱负的结果；目前，可选用的小波函数很多，如Morlet 小 波和 Meyer 小波等； 在本例中， 我们选用 Morlet 连续复小波变换来分析径流时间序列的多时间尺度特点；缘由如下：1.1 径流演化过程中包含“多时间尺度”变化特点且这种变化是连续的，所以应采纳连续小波变换来进行此项分析；1.2 实小波变换只能给出时间序列变化的振幅和正负，而复小波变换可同时给出时间序列变化的位相和振幅两方面的信息，有利于对问题的进一步分析；1.3 复小波函数的实部和虚部位相差为使分析结果更为精确；π /，2 能够排除用实小波变换系数作为判据而产生的虚假振荡，2. 绘制小波系数图、小波方差图和主周期变化趋势图是关键当挑选好合适的基小波函数后，下一步的关键就是如何通过小波变换获得小波系数，然后利用相关软件绘制小波系数图、小波方差图和主周期变化趋势图，进而依据上述三种图形的变化识别径流时间序列中存在的多时间尺度；具体步骤1.2.3.4.5.6.7.数据格式的转化边界效应的排除或减小运算小波系数运算复小波系数的实部、模、模方、方差绘制小波系数实部、模、模方等值线图 绘制小波方差图绘制主周期趋势图Matlab R2021a 、suffer 12.0 、Excel、记事本等，具体说明小波系数下面，我们以上题为例，结合软件的运算和各图形的绘制过程，并分别说明各图在分析径流多时间尺度变化特点中的作用；1. 数据格式的转化和储存将存放在 Excel 表格里的径流数据 （以时间为序排为一列） 转化为 Matlab R2021a 识别的数据格式 （ .mat）并存盘；具体操作为： 在 Matlab R2021a 界面下， 单击“ File-Import Data”，显现文件挑选对话框 “ Import ”后，（本例的文件名为 runoff.xls ），单击“打开” ；等数据转化完成后， 单击“ Finish”，找到需要转化的数据文件显现图 1 显示界面；然后双击图1 中的 Runoff ，弹出“ Array Editor: runoff ”对话框，挑选 File 文件夹下的“ Save Workspace As ”单击，显现图（ runoff.mat ），单击“储存”并关闭“2 所示的“ Save to MAT-File:”窗口，挑选存放路径并填写文件名Save to MAT-File”窗口；图 1 数据格式的转化2. 边界效应的排除或减小图 2 数据的储存由于本例中的实测径流数据为有限时间数据序列，在时间序列的两端可能会产生“边界效用”；为消除或减小序列开头点和终止点邻近的边界效应，须对其两端数据进行延长；在进行完小波变换后，去掉两Matlab R2021a 小波工端延长数据的小变换系数，保留原数据序列时段内的小波系数；本例中，我们利用具箱中的信号延长（ Signal Extension ）功能，对径流数据两端进行对称性延长；具体方法为： 在 Matlab R2021a 界面的“ Command Window ”中输入小波工具箱调用命令“ Wavemenu”，按 Enter 键 弹“ Wavelet Toolbox Main Menu”（ 小波工具箱主菜单） 界面（图３）；然后单击“Signal Extension ”，打 开 Signal Extension / Truncation窗口，单击“ File ”菜单下的“ Load Signal ”，挑选 runoff.mat 文件单击“打开”，显现图 4 信号延长界面；Matlab R2021a 的 Extension Mode 菜单下包含了 6 种基本的延长方式（ Symmetric 、Periodic 、Zero Padding、Continuous 、Smooth and For SWT ）和 Direction to extend 菜单下的3 种延长模式（ Both 、Left and Right ），在这里我们挑选对称性两端延长进行运算；数据延长的具体操作过程是： Desired Length 可以任意选，只要比原始信号长度大，建议在原始信号的基础上加20（这样左右对称地延长 10 个数据），这里挑选默认的64； Dircetion to extend 下挑选“ Both ”； Extension Mode 下挑选“ Symmetric ”；单击“ Extend”按钮进行对称性两端延长运算，然后单击 “ File ”菜单下的 “ Save Tra。