1、小波分析 原理与应用 Niu 二哥 CUIT 3S 集成 1 / 9 小波分析理论分析多时间尺度变化特征 时间序列 （ Time Series） 是地学研究中经常遇到的问题 。在时间序列 研究中， 时域和频域是常用的两种 基本 形式 。 其中， 时域分析具有时间定位能 力 ，但无法得到关于时间序列变化 的 更多信息 ； 频域分析 （ 如Fourier 变换） 虽具有准确的频率定位功能，但仅适合平稳时间序列分析。 然而 ， 地学中许多现象 （如河川径流、 地震波 、 暴雨 、 洪水 等） 随时间的变化 往往 受到多种因素的 综合 影响， 大都属于非平稳序列，它们不但具有 趋势性、周期性 等 特征
2、， 还 存在 随机性、突变性 以及 “ 多时间尺度 ” 结构，具有多层次演变规律 。对于这类非平 稳时间序列的研究，通常需要某一频段对应的时间信息，或某一时段的频域信息。 显然，时域分析和频域分析 对此 均无能为力。 20 世纪 80 年 代初，由 Morlet 提出的 一种 具有时 -频多分辨功能的 小波分析（ Wavelet Analysis） 为 更好的 研究 时间序列问题提供了可能， 它 能 清 晰 的揭示出隐藏在时间序列中的 多 种变化周期， 充分 反映系统在不同时间尺度中的变化趋势，并能对系统未来发展趋势进行定性估计。 目前， 小波分析理论已在信号处理、图像压缩、模式识别、数值分析
3、和大气科学等众多的非线性科学领域内 得到了广泛的应 。 在时间序列研究中 ， 小波分析 主要用于 时间序列的 消噪和滤波，信息量系数和分形维数的计算，突变点的监测和 周期 成分的识别以及多时间 尺度 的分析等。 一、小波分析基本原理 1. 小波函数 小波分析的基本思想是用一簇小波函数系来表示 或逼近 某一信号 或函数。因此，小波函数是小波分析的关键，它 是指具有震荡性、能够迅速衰减到零的一类函数， 即小波函数 )R(L)t( 2 且 满足 ： 0dt)t( （ 1） 式中 ， )t( 为 基 小波 函数 ， 它 可 通过尺度的伸缩和时间轴上的平移 构成一簇函数系： )a bt(a)t( 2/1
4、b,a 其中， 0aR,ba, （ 2） 式中， )t(b,a 为子小波； a 为尺度因子， 反映小波的周期长度； b 为平移因子， 反应时间上的平移 。 需要说明的是，选择合适的基小波函数是进行 小波分析的前提。 在实际应用研究中，应针对具体情况选择所需的基小波 函数 ；同一 信号或 时间序列，若选择 不同的基小波函数， 所得的 结果 往往会 有所差异 ，有时甚至差异很大 。 目前，主要是通过对比不同小波分析处理信号时所得的结果与理论结果的误差来判定基小波函数的好坏，并由此选定该类研究所需的基小波函数。 2. 小波变换 若 )t(b,a 是 由 （ 2）式给出的子小波，对 于 给定的能量有限
5、信号 )R(L)t(f 2 ，其连续小波变换（ Continue Wavelet Transform，简写为 CWT） 为： dt)a bt(f ( t )a)b,a(W R2/1-f （ 3） 式中， )b,a(Wf 为小波变换系数； f(t)为一个信号或平方可积函数 ； a 为伸缩尺度 ； b 平移参数 ； )abx( 为 )abx( 的复共轭 函数。 地学中 观测到的 时间序列 数据 大多是离散的， 设 函数 )tk(f ，（ k=1,2, ,N; t为取样间隔）， 则 式（ 3）的离散小波变换形式为： 小波分析 原理与应用 Niu 二哥 CUIT 3S 集成 2 / 9 )a b-tk
6、(t)f(kta)b,a(W N 1k2/1-f （ 4） 由式（ 3）或 （ 4）可知小波分析的基本原理，即 通过 增加或减小 伸缩尺度 a 来 得到信号的 低频或高 频信息 ，然后 分析信号的概貌 或 细节，实现对信号 不同时间尺度和空间局部特征的 分析。 实际研究中 ，最主要的就是 要 由小波变换方程得到小波系数 ，然后通过这些系数 来 分析时间序列的时频变化特征 。 3. 小波方差 将小波系数的平方值在 b 域上积分，就可得到小波方差，即 db)ba,(W)a(V a r 2f (5) 小波方差随尺 度 a 的变化过程，称为小波方差图。由式 (5)可知，它 能 反映 信号 波动的能量随
7、尺度 a 的分布。因此，小波方差图 可用来 确定 信号 中 不同 种尺度扰动的相对强度 和 存在的主要时间尺度，即主周期。 二、小波分析实例 -时间序列的多时间尺度分析 (Multi-time scale analysis) 河川径流是地理水文学研究 中 的一个重要变量 ，而多时间尺度是径流演化过程中 存在 的重要特征。所谓径流时间序列的多时间尺度是指：河川径流在演化过程中，并不存在 真正意义上的变化 周期， 而是其变化周期随着研究尺度的不同而 发生相应的 变化， 这种变化一 般 表现为小时间尺度的变化周期 往往 嵌套在大尺度的变化周期之中。 也就是说，径流变化在时间域中存在多层次的时间尺度结
8、构和局部变化特征。 表 1 给出了某流域某水文观测站 1966-2004 年的实测径流数据。试运用小波分析理论 ， 借助 Matlab6.5、suffer8.0 和相关软件 （ Excel 等 ） ， 完成下 述任务 ： 计算 小波系数； 绘制 小波系数图 （实部、模和模方） 、小波方差图和主周期变化趋势图 ，并分别说明各图在 分析 径流多时间尺度变化特征 中 的作用 。 表 1 某流域某水 文 观测 站 1966-2004 年 实测 径流 数据（ 108m3） 年份 径流量 年份 径流量 年份 径流量 年份 径流量 年份 径流量 1966 1.438 1974 2.235 1982 0.77
9、4 1990 1.806 1998 1.709 1967 1.151 1975 4.374 1983 0.367 1991 0.449 1999 0.000 1968 0.536 1976 4.219 1984 0.562 1992 0.120 2000 0.000 1969 1.470 1977 2.590 1985 3.040 1993 0.627 2001 2.104 1970 3.476 1978 3.350 1986 0.304 1994 1.658 2002 0.009 1971 4.068 1979 2.540 1987 0.728 1995 1.025 2003 3.177 1
10、972 2.147 1980 0.807 1988 0.492 1996 0.955 2004 0.921 1973 3.931 1981 0.573 1989 0.007 1997 1.341 分析 1. 选择合适的基小波函数 是 前提 在运用小波分析理论 解决 实际问题时， 选择合适的 基 小波函数 是前提。只有选择了适合具体问题的基小波函数，才能得到较为理想的结果。目前，可选用的小波函数很多，如 Mexican hat 小波、 Haar 小波、Morlet 小波和 Meyer 小波等。在本例中，我们选用 Morlet 连续复 小波变换来分析径流时间序列的多时间尺度特征 。原因如下： 1.
11、1 径流演变过程中包含“多时间尺度”变化特征且这种变化是连续的，所以应采用连续小波变换来进行此项分析。 1.2 实小波变换只能给出 时间序列变化的 振幅和正负 ，而 复 小波变换可同时给出时间序列变化的位相和振幅两方面的信息 ，有利于对问题的进一步分析。 小波分析 原理与应用 Niu 二哥 CUIT 3S 集成 3 / 9 1.3 复小波函数的实部和虚部位相差为 /2，能够消除用实小波变换系数作为判据而产生的虚假振荡，使分析结果更 为 准确。 2. 绘制 小波系数 图 、小波方差图 和主周期变化趋势图 是关键 当选择好合适的基小波函数后，下一步的关键就是如何通过小波变换获得小波系数，然后利用
12、相关软件绘制小波系数图、小波方差图和主周期变化趋势图， 进而 根据上述三种图形的变化识别径流 时间序列中存在的多时间尺度。 三、 具体 步骤 1. 数据格式 的 转化 2. 边界效应的消除 或 减小 3. 计算 小波系数 4. 计算复小波系数的实部 5. 绘制小波系数 实部等值线 图 6. 绘制小波系数模 和 模方等值线图 7. 绘制小波方差图 8. 绘制主周期趋势图 下面，我们 以上题为例， 结合软件 Matlab 6.5、 Suffer 8.0 和 Excel， 详细 说明小波系数的计算和 各 图形的绘制过程， 并分别说明各图在分析径流多时间尺度变化特征 中的 作用。 1. 数据格式 的
13、转化 和保存 将存放在 Excel 表格里的 径流 数据 （以时间为序排为一列）转化为 Matlab 6.5 识别的数据格式（ .mat）并存盘。 具体操作为： 在 Matlab 6.5 界面下， 单击 “ File-Import Data” ，出现 文件选择对话框 “ Import” 后，找到需要转化的数据文件 （ 本例的文件名为 runoff.xls） ， 单击 “ 打开 ” 。等数据转化完成后， 单击“ Finish” ，出现图 1 显示界面； 然后双击图 1 中的 Runoff， 弹出 “ Array Editor: runoff” 对话框，选择 File 文件夹下的“ Save Wo
14、rkspace As” 单击，出现 图 2 所示的“ Save to MAT-File:”窗口 ，选择存放路径 并填写 文件名（ runoff.mat），单击 “ 保存 ” 并关闭 “ Save to MAT-File” 窗口。 2. 边界效应的消除 或 减 小 因为本例中的实测径流数据为有限 时间 数据 序列 ，在 时间 序列 的两端可能会产生“边界效用”。为消除或减小 序列 开始点和结束点附近的边界效应，须对其两端数据进行延伸。在进行完小波变换后，去掉两端延伸数据的小变换系数，保留原数据序列时段内的小波系数。本例中，我们 利用 Matlab 6.5 小波 工具箱中的信号延伸 （ Signa
15、l Extension） 功能 ，对径流数据 两端 进行 对称性 延伸。 图 1 数据格式的转化 图 2 数据的保存 小波分析 原理与应用 Niu 二哥 CUIT 3S 集成 4 / 9 具体 方法 为 ：在 Matlab 6.5 界面的 “ Command Window” 中输入小波工具箱调用命令“ Wavemenu”，按 Enter 键弹 “ Wavelet Toolbox Main Menu”（ 小波工具箱主菜单 ） 界面（ 图 ） ；然后单击 “ Signal Extension” ，打开 Signal Extension / Truncation 窗口，单击 “ File”菜单下的“
16、 Load Signal” ，选择 runoff.mat 文件单击“ 打开 ” ，出现图 4 信号延伸界面 。 Matlab 6.5 的 Extension Mode 菜单下 包含 了 6 种 基本 的 延伸方式（ Symmetric、 Periodic、 Zero Padding、 Continuous、 Smooth and For SWT） 和 Direction to extend 菜单下 的3 种延伸模式（ Both、 Left and Right） ， 在这里我们选择 对称性 两端延伸 进行计算 。 数据延伸的具体操作过程 是 ： 在 Extension Mode 下选择 “ Sy
17、mmetric” ， Dircetion to extend 下选择 “ Both” ， 单击 “ Extend”按钮进行对称性两端延伸计算，然后单击 “ File”菜单下的“ Save Tranformed Signal”， 将延伸后的 数据 结果存为 erunoff.mat 文件。 从 erunoff 文件可知，系统自动将 原 时间序列 数据向 前 对称 延伸 12 个 单位 ，向后延伸 13 个 单位 。 3. 计算小波系数 选择 Matlab 6.5 小波 工具箱 中的 Morlet 复 小波 函数 对延伸后的径流数据 序列（ erunoff.mat） 进行小波变换，计算小波系数并存盘
18、。 小波工具箱主菜单界面见图 3，单击 “ Wavelet 1-D” 下的子菜单 “ Complex Continuous Wavelet 1-D” ，打开一维复连续小波界面，单击 “ File”菜单下的“ Load Signal”按钮， 载入径流时间序列 erunoff.mat（图 5） 。图 5 的 左侧为 信号 显示区域， 右 侧区域 给出了信号序列和复小波变换的有关信息和参数， 主要包括 数据长度 （ Data Size） 、小波函数类型（ Wavelet： cgau、 shan、fbsp 和 cmor）、取样周期 （ Sampling Period） 、周期设置 （ Scale Se
19、tting） 和 运行 按钮 （ Analyze） ， 以及显示区域的相 关显示设置按钮 。 本例中，我们选择cmor (1-1.5)、取样周期为 1、最大尺度为 32，单击“ Analyze” 运行按钮 ， 计算小波系数。 然后单击 “ File”菜单下的“ Save Coefficients” ，保存小波系数为cerunoff.mat 文件。 注意 ： 上面涉及到的 数据保存 ，其 格式均为 .mat。 4. 计算 Morlet 复小波系数的实部 将复小波系数转存到 Excel 表格，去掉两端延伸数据的小波系数，并计算小波系数实部。 在 Matlab 6.5 界面 下 的 Workspac
20、e 中将 cerunoff.mat 文件导入，然后 双击打开，全部复制到 Excel图 3 小波工具箱主菜单 图 4 径流时间序列的延伸 图 5 小波变换菜单界面 小波分析 原理与应用 Niu 二哥 CUIT 3S 集成 5 / 9 后去掉延伸数据的小波变换系数 （ 本例中 去掉 前 12 列和后 13 列） ， 或只复制原时间序列的小波变换系数到 Excel， 最后使 用 Excel 中的 IMREAL 函数计算原时间序列的小波 系数实部 （ 图 6） 。 图 6 复小波系数及实部计算 示意 图 Excel 中 IMREAL 函数的调用方法为：单击 “ 插入 ” 菜单下的“函数（ F）”按钮
21、，弹出图 7 所示的“插入函数”窗口， 在“搜索函数（ S）：”框中输入：“ IMREAL”后单击 “ 转到 ” ，再单击 “ 确定 ” ，出现函数参数窗口（图 8） 。 在 “ Inumber”一 栏 的 空白处 输入所要计算的数据 （图 6） ，单击 “ 确定 ” 即可 得到小波系数实部值 。 图 7 IMEAL 函数调用 图 8 IMEAL 函数 参数 小波分析 原理与应用 Niu 二哥 CUIT 3S 集成 6 / 9 需要说明的是， 从 cerunoff.mat 文件中转到 Excel 里的复小波系数， 在其实部和虚部 中间 包含许多“空格”，在计算之前需要先将其去掉。 5. 借助
22、Suffer 8.0， 绘制小波系数实部等值线图 5.1 小波系数实部等值线图的绘制 首先， 将小波系数实部数据 按照图 9 格式排列，其中列 A 为时间，列 B 为尺度，列 C 为不同时间和 尺度下所对应的小波系数实部值。 图 9 小波系数实部 数据格式 其次，将图 9 数据转化成 Suffer 8.0 识别的数据格式。具体操作为：在 Suffer 8.0 界面下，单击“网格”菜单下的“数据”按钮，在“打开”窗口选择要打开的文件（小波系数实部 .xls），单击“打开”后弹出“ 网格化数据 ” 对话框（图 10） 。 它给出了多种不同的网格化方法、文件输出路径及网格线索几何学 等信息 。 这里
23、 我们 选择 “ 克里格 “ 网格方法 ” ， 单 击 “ 确定 ” ，完成数据格式的转化 。 最后， 绘制 小波 系数实部 等值线 图。在 Suffer 8.0 界面下，单击 “ 地图 ” 菜单下的 “ 等值线图 -新建等图 10 小波系数实部数据格式转化 图 11 Suffer8.0 中的小 系数实部等值线图 小波分析 原理与应用 Niu 二哥 CUIT 3S 集成 7 / 9 值线图 ” 按钮 ， 弹出“打开网格”窗口 后 ，选择 “ 小波系数实部 .grd” 文件，单击 “ 打开 ” ，完成等值线图的绘制 并存盘 （图 11） 。 5.2 小波系数实部等值线图在多时间尺度 分析中的作用
24、 小波系数实部等值线图能 反映 径流序列不同 时间尺度的周期变化及其在时间域中的分布， 进而 能 判断在不同时间尺度上 ， 径流的 未来变化趋势 。 为 能 比较清楚的说明小波系数实部等值线图在 径流 多时间尺度分析中的作用，我们利用 Suffer 8.0 对其 进一步处理 和 修饰，得到图 12 显示的小波系数实部等值线图 。 其中，横坐标为时间 （年份） ，纵坐标为时间尺度，图中的等值曲线为小波系数实部值。 当 小波系数实部 值 为正时，代表 径流 丰水期， 在 图中 我们 用实线绘出， “H”表示正值中心；为负时，表示 径流 枯水期，用虚线绘出， “L”表示负值中心。 由图 12 可以
25、清楚的看到径流演化过程中存在的多时间尺度特征。 总的来说， 在 流域 径流 演变过程中 存在 着1832 年， 817 年以及 37 年的 3 类尺度的周期变化规律。其中， 在 1832 年 尺度上 出现 了 枯 -丰交替的准两次震荡 ；在 817 年时间尺度上 存在 准 5 次震荡 。 同时，还可以看出以上两个尺度的周期变化在整个分析时段表现的非常稳定，具有全域性 ； 而 310 年尺度的周期变化， 在 1980s 以后 表现的较为稳定 。 6. 绘制小波系数模和模方等值线图 6.1 小波系数模和模方等值线图的绘制 参考 4、 5 两步，绘制小波系数模 和 模方等值线图 （图 13、 14）
26、 。 说明：在 Excel 中，复数模的计算函数为 “ IMABS” 。 1970 1975 1980 1985 1990 1995 200051015202530(b)年份197 0 1975 19 80 198 5 1990 1995 200 051015202530(c)年份 图 13 小波系数模等值线图 图 14 小波系数模方等值线图 6.2 小波系数模 等值线图 在多时间尺度分析中的作用 Morlet 小波系数的模值是不同时间 尺度变化周期 所对应的 能量密度在时间域中分布的反映， 系数模值愈大，表明其所对应时段或尺度的周期性就愈强 。从图 13 可 以看出 ， 在流域径流演化过程中
27、， 1832 年时间尺度模值最大，说明该时间尺度周期变化最明显， 1822 年时间尺度的周期变化次之，其他时间尺度的周期性变化较小 ； 6.2 小波系数模方等值线图在多时间尺度分析中的作用 小波系数 的 模方相当于小波能量谱， 它可以 分析 出 不同周期的震荡能量。由图 14 知 ， 2532 年时间尺1970 1975 1980 1985 1990 1995 200051015202530(a)年份H HH HHH HHLLLLLLLL图 12 小系数实部等值线图 小波分析 原理与应用 Niu 二哥 CUIT 3S 集成 8 / 9 度的 能量最强 、 周期 最 显著 ，但 它的 周期变化具
28、有局部性 （ 1980s 前） ； 1015 年时间尺度能量 虽然较弱，但周期分布比较明显，几乎占据整个研究时域 (19742004 年 )。 7. 绘制小波方差图 7.1 小波方差图的绘制 将不同时间尺度下的小波系数代入 式 (5)可得径流变化的小波方差 ， 以小波方差为纵坐标，时间尺度 a为横坐标， 可 绘制小波方差图 （图 15） 。 7.2 小波 方差图 在多时间尺度分析中的作用 小波方差图能 反映 径流时间序列的 波动能量随尺度 a 的分布 情况。可用来 确定 径流演化过程 中 存在的主周期 。 流域径流的小波方差图中（ 图 15）存在 4 个 较为明显的 峰值， 它们依次 对应着
29、28 年、 14 年、 8 年和 4年的时间尺度 。其中， 最大峰值 对应着 28 年 的 时间尺度，说明 28 年左右的周期震荡最强，为流域年径流变化的第一主周期 ； 14 年时间尺度对应着第二峰值，为径流变化的第二主周期， 第三、第三峰值分别对应着8 年和 4 年 的 时间尺度 ，它们 依次为 流域径流的第三和第四主周期 。 这说明 上述 4 个周期的波动控制着流域径流在整个时间域内的变化特征。 气象家园提问截图，小波方差图： http:/ 8. 主周期趋势图 的绘制及 其 在多时间尺度分析中的作用 根据小波方差 检验 的 结果 ， 我们 绘制 出了控制流域径流演变的 第一和第二主周期小波
30、系数图（图 16） 。从主周期趋势图中我们可以分析出在不同的时间尺度下，流域径流存在的平均周期及丰 -枯变化特征。 图16a 显示 ，在 14 年特征时间尺度上， 流域径流变化的平均周期为 9.5 年左右， 大约经历了 4 个丰 -枯转换期 ；而在 28 年特征时间尺度上（图 16b）， 流域的平均变化周期为 20 年左右， 大约 2 个周期的丰 -枯变化 。 (d)0204060801001201400 5 10 15 20 25 30 35时间尺度/a小波方差图 15 小波方差图 小波分析 原理与应用 Niu 二哥 CUIT 3S 集成 9 / 9 a 14特征 时间尺度-1-0.8-0.
31、6-0.4-0.200.20.40.60.811965 1970 1975 1980 1985 1990 1995 2000 2005年份小波系数b 28特征 时间尺度-3-2-101231965 1970 1975 1980 1985 1990 1995 2000 2005年份小波系数图 16 大沽夹河流域年径流变化的 13 年和 28 年特征时间尺度小波实部过程线 参考文献 王文圣，丁晶，李耀清 . 2005. 水文小波分析 M. 北京：化学工业出版社 曹素华等 . 1998. 实用医学多因素统计方法 M. 上海：上海医科大学出版社 方开泰 . 1989. 实用多元统计分析 M. 上海：华
32、东师范大学出版社 何清波，苏炳华，钱亢 . 2002. 医学统计学及其软件包 M. 上海：上海科学技术文献出版社 胡秉民 . 1987. 微电脑在农业科学中的应用 M. 北京：科学出版社 孙尚拱 . 1990 实用多元变量统计方法与计算程序 M. 北京：北京医科大学、中国协和医科大学联合出版社 唐守正 . 1986. 多元统计分析方法 M。北京：中国林业出版社 王学仁 . 1982. 地址数据的多变量统计分析 . 北京：科学出版社 徐振邦，金淳浩，娄元仁 . 1986. 2 距离系数和 2 距离系 数尺度在聚类 分析中的应用 M. 赵旭东等主编，中国数学地质（ 1） . 北京：地质出版社 於崇
33、文 . 1978. 数学地质的方法与应用 M. 北京：冶金工业出版社 Anderson T. W. 1967. Introduction to multivariate statistical analysis, 2ndM. New York: Wiley Gauch H. G. J. 1982. Multivariate analysis in community ecologyM. Britain: Cambridge University Press Horel A. E. ,Wennard. R. W. and Baldwin K. F. 1975. regression: some
34、simulations. Communications in StatisticsJ, 4: 105123 练习 试运用小波分析理论， 分析某 市 年 平均 降水 过程中存在 的多时间尺度 变化 特征。 表 2 某市 1957-2004 年实测年均降水量（ mm） 年份 降水量 年份 降水量 年份 降水量 年份 降水量 1957 320.0 1969 324.8 1981 506.0 1993 384.4 1958 481.2 1970 412.3 1982 282.1 1994 503.9 1959 522.6 1971 366.5 1983 508.6 1995 406.7 1960 33
35、9.3 1972 262.4 1984 523.9 1996 465.7 1961 719.9 1973 521.9 1985 518.9 1997 345.3 1962 373.5 1974 351.7 1986 320.1 1998 454.9 1963 332.9 1975 398.4 1987 340.0 1999 327.9 1964 741.2 1976 320.2 1988 478.5 2000 406.2 1965 454.3 1977 445.4 1989 402.4 2001 404.7 1966 604.3 1978 534.8 1990 552.4 2002 401.9 1967 451.9 1979 509.7 1991 313.9 2003 605.1 1968 424.3 1980 395.5 1992 591.0 2004 385.4
