Matlab: 小波分析—时间序列的多时间尺度分析（一）

一则小波分析帖子在我的科学网博客上获得了非常大的关注，至今这一帖子保持着留言询问最多的地位。不过，由于2015年11月的一场浩劫，这个帖子附带的示例代码就都烟消云散了，其后，还是有很多同学留言询问，我在这里把代码再回顾一遍，文字基本不变。

本部分对应——MATLAB：小波分析—时间序列的多时间尺度分析。

Introduction

时间序列（Time Series）是地学研究中经常遇到的问题。在时间序列研究中，时域和频域是常用的两种基本形式。其中，时域分析具有时间定位能力，但无法得到关于时间序列变化的更多信息；频域分析（如Fourier变换）虽具有准确的频率定位功能，但仅适合平稳时间序列分析。然而，地学中许多现象（如河川径流、地震波、暴雨、洪水等）随时间的变化往往受到多种因素的综合影响，大都属于非平稳序列，它们不但具有趋势性、周期性等特征，还存在随机性、突变性以及“多时间尺度”结构，具有多层次演变规律。对于这类非平稳时间序列的研究，通常需要某一频段对应的时间信息，或某一时段的频域信息。显然，时域分析和频域分析对此均无能为力。

20世纪80年代初，由Morlet提出的一种具有时—频多分辨功能的小波分析（Wavelet Analysis）为更好的研究时间序列问题提供了可能，它能清晰的揭示出隐藏在时间序列中的多种变化周期，充分反映系统在不同时间尺度中的变化趋势，并能对系统未来发展趋势进行定性估计。

本帖以美国某气象站1894~2010年连续的年降水量为例，试用小波分析，完成如下任务：①小波变换系数；②绘制小波系数实部等值线图；③绘制小波系数模和模方等值线图；④绘制小波方差图；以及⑤绘制不同时间尺度的小波实部过程线。所谓年降水量时间序列的多时间尺度是指：年降水量在演化过程中，并不存在真正意义上的变化周期，而是其变化周期随着研究尺度的不同而发生相应的变化，这种变化一般表现为小时间尺度的变化周期往往嵌套在大尺度的变化周期之中。也就是说，年降水量变化在时间域中存在多层次的时间尺度结构和局部变化特征。

小波分析的计算过程请参考：小波分析—经典、小波方差制作步骤等。

小波分析的基本理论在此不多叙述，请参考其他文献。本帖主要介绍小波分析的一般过程，数据及部分代码将附在文末。

0．年降水量的变化趋势分析

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

Fig. 0

1．数据的加载

加载多年降水量至MATLAB使用Load命令。

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

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

具体步骤：在Matlab界面的“Command Window”中输入小波工具箱调用命令“wavemenu”，按Enter键弹“Wavelet Toolbox Main Menu”（小波工具箱主菜单）界面（Fig. 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”按钮进行对称性两端延伸计算（Fig. 1下），然后单击“File”菜单下的“Save Tranformed Signal”，将延伸后的数据结果存为ePrec.mat文件。

Fig. 1

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

3．计算小波系数

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

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

Fig. 2

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

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

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

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

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

Fig. 3

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

红框处出现4个偏多中心和3个偏少中心，一对临近的偏多和偏少波动就可以被识别为一个周期，所以是28~36a主震荡周期。

周期判断可以从正弦曲线上类比，如下图，两个最邻近波峰（或波谷）之间是一个完整的周期，或者两个最邻近的与x轴向下（或向上）交点也是一个完整的周期。

那么按照正弦曲线的周期检测思路，两个最邻近的偏多中心（或偏少中心）就是一个完整的周期长度。

File.rar.

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

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

Fig. 4

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

Fig. 5

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

7．绘制小波方差图

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

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

Fig. 6

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

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

Fig. 7

Fig. 8

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

File.rar.

More

[1] A Practical Guide to Wavelet Analysis.

[2] 初学应用举例（小波分析，实部画图，小波方差画图及小波模画图等）.