时间序列小波分析的操作步骤及实例分析

　　小波分析实例

　　时间序列（Time Series）是地学研究中经常遇到的问题。在时间序列研究中，时域和频域是常用的两种基本形式。其中，时域分析具有时间定位能力，但无法得到关于时间序列变化的更多信息；频域分析（如Fourier变换）虽具有准确的频率定位功能，但仅适合平稳时间序列分析。

　　河川径流是地理水文学研究中的一个重要变量，而多时间尺度是径流演化过程中存在的重要特征。所谓径流时间序列的多时间尺度是指：河川径流在演化过程中，并不存在真正意义上的变化周期，而是其变化周期随着研究尺度的不同而发生相应的变化，这种变化一般表现为小时间尺度的变化周期往往嵌套在大尺度的变化周期之中。也就是说，径流变化在时间域中存在多层次的时间尺度结构和局部变化特征。

　　年降水量的变化趋势分析

　　该站点的年降水量变化情况（LI_plot函数）如图 0所示，其发展呈微微上升的趋势，降水量最高年份出现在2003年，全年累计达到1610.70mm，最低值出现在1965年，累计仅有748.90mm。

　　

　　具体步骤

　　1）数据格式的转化

　　2）边界效应的消除或减小

　　3）计算小波系数

　　4）计算复小波系数的实部

　　5）绘制小波系数实部等值线图

　　6）绘制小波系数模和模方等值线图

　　7）绘制小波方差图

　　8） 绘制主周期趋势图

　　1. 数据格式的转化和保存

　　将存放在Excel表格里的径流数据（以时间为序排为一列）转化为Matlab 6.5识别的数据格式（.mat）并存盘。

　　具体操作为：在Matlab 6.5 界面下，单击“File-Import Data”，出现文件选择对话框“Import”后，找到需要转化的数据文件（本例的文件名为runoff.xls），单击“打开”。等数据转化完成后，单击“Finish”，出现图1显示界面；然后双击图1中的Runoff，弹出“Array Editor： runoff”对话框，选择File文件夹下的“Save Workspace As”单击，出现图2所示的“Save to MAT-File：”窗口，选择存放路径并填写文件名（runoff.mat），单击“保存”并关闭“Save to MAT-File”窗口。

　　

　　2．边界效应的消除或减小

　　由于本例中的实测年降水量数据为有限时间数据序列，在时间序列的两端可能会产生“边界效用”。为消除或减小序列开始点和结束点附近的边界效应，须对其两端数据进行延伸。在进行完小波变换后，去掉两端延伸数据的小波变换系数，保留原数据序列时段内的小波系数。本例中，我们利用Matlab小波工具箱中的信号延伸（SignalExtension）功能，对数据两端进行对称性延伸。

　　具体步骤：在Matlab界面的“Command Window”中输入小波工具箱调用命令“wavemenu”，按Enter键弹“Wavelet Toolbox Main Menu”（小波工具箱主菜单）界面（图 1上）；然后单击“Signal Extension”，打开Signal Extension /Truncation窗口，单击“File”菜单下的“Load Signal”，选择Prec.mat文件单击“打开”，出现信号延伸界面。Matlab的Extension Mode菜单下包含多种延伸方式和Direction to extend菜单下的3种延伸模式（Both、Left and Right），在这里我们选择对称性两端延伸进行计算。具体操作过程是：在Extension Mode下选择“Symmetric（Whole-Point）”，Dircetion to extend下选择“Both”，单击“Extend”按钮进行对称性两端延伸计算（图 1下），然后单击“File”菜单下的“Save Tranformed Signal”，将延伸后的数据结果存为ePrec.mat文件。

　　从ePrec文件可知（entendvalue函数），系统自动将原时间序列数据向前对称延伸5个单位，向后延伸6个单位。

　　图 1

　　3．计算小波系数

　　选择Matlab小波工具箱中的Morlet复小波函数对延伸后的数据序列（ePrec.mat）进行小波变换，计算小波系数并保存。

　　小波工具箱主菜单界面见图1上，单击“One-Dimensional”下的子菜单“Complex Continuous Wavelet 1-D”，打开一维复连续小波界面，单击“File”菜单下的“Load Signal”按钮，载入时间序列数据ePrec.mat。图 2第一排的左侧为信号显示区域，右侧区域给出了信号序列和复小波变换的有关信息和参数，主要包括数据长度（Data Size）、小波函数类型（Wavelet：cgau、shan、fbsp和cmor）、取样周期（Sampling Period）、周期设置（ScaleSetting）和运行按钮（Analyze），以及显示区域的相关显示设置按钮。本例中，我们选择cmor （1-1.5）、取样周期为1、最大尺度为64（延伸后的时间序列的一半），单击“Analyze”运行按钮，计算小波系数。然后单击“File”菜单下的“Save Coefficients”，保存小波系数为cePrec.mat文件。

　　

　　图 2

　　4．计算Morlet复小波系数的实部

　　去除两端延伸数据的小波系数（entendvalue函数），并计算小波系数实部（real函数）。

　　5．绘制小波系数实部等值线图

　　这部分过程应用LI_contourf函数，年降水量小波系数实部等值线图见图 3。

　　如图3所示的小波系数实部等值线图。其中，横坐标为时间（年份），纵坐标为时间尺度，图中的等值曲线为小波系数实部值。小波系数实部等值线图能反映年降水量序列不同时间尺度的周期变化及其在时间域中的分布，进而能判断在不同时间尺度上，年降水量的未来变化趋势。

　　

　　图 3

　　由图 3可以看出降水量1894~2010年演化过程中存在着14~18a的主振荡周期（这一周期是多个振荡周期叠加的矢量和，随后将在小波方差图中对这一主周期进行分解，剥离出第一主周期等等）。在整个时间尺度上出现4个偏多中心和3个偏少中心，分别为1899、1938、1972、2007年和1919、1956、1990年。

　　6．绘制小波系数模和模方等值线图

　　首先，计算小波系数的模和模方，再利用LI_contourf函数绘制模和模方等值线图。

　　Morlet小波系数的模值是不同时间尺度变化周期所对应的能量密度在时间域中分布的反映，系数模值愈大，表明其所对应时段或尺度的周期性就愈强。从图 4可以看出，在降水量演化过程中，40~64年时间尺度模值最大（大于400），但1924~1994年之间模值小余200，说明在此时段内40~64年时间尺度的周期变化并不明显，1995年之后模值再次增大，此时期后40~64年时间尺度的周期变化趋于显著。

　　小波系数的模方相当于小波能量谱，它可以分析出不同周期的震荡能量。由图 5可知，40~64年时间尺度的能量最强、周期最显著，但它的周期变化具有局部性（1924年之前和1994年之后）；10~15年时间尺度能量虽然较弱，但周期分布比较明显，几乎占据整个研究时域（1894~2010年）。

　　

　　图 4

　　

　　图 5

　　7．绘制小波方差图

　　小波方差的计算过程参考：小波方差制作步骤（参考文献存在一个错误，小波方差未除N，本贴小波方差计算已有修正）。

　　小波方差图能反映降水量时间序列的波动能量随尺度（a）的分布情况。可用来确定降水量演化过程中存在的主周期。

　　

　　图 6

　　小波方差图中（图 6）存在4个较为明显的峰值，它们依次从小至大对应着8a、19a、30a和55a的时间尺度。其中，最大峰值对应着55a的时间尺度，说明55a左右（时间尺度）的周期震荡最强，为年降水量变化的第一主周期（注意：不是55a是第一主周期）；30a时间尺度对应着第二峰值，为第二主周期；第三、第四峰值分别对应着19a和8a的时间尺度，它们依次为降水量的第三和第四主周期。这说明上述4个周期的波动控制着降水量在整个时间域内的变化特征。

　　8．主周期趋势图的绘制及其在多时间尺度分析中的作用

　　根据小波方差检验的结果，绘制演变的第一和第二主周期小波系数图。

　　从主周期趋势图中可以分析出在不同的时间尺度下，降水量存在的平均周期及丰-枯变化特征。图 7显示，在55a特征时间尺度上，降水量变化的平均周期为35a左右，大约经历了3个丰—枯转换期；而在30a特征时间尺度上（图 8），平均变化周期为20a左右，大约6个周期的丰—枯变化。

　　

　　图 7

　　

　　图 8