一种水文时间序列小波互相关分析方法 【技术领域】
本发明涉及一种水文时间序列分析方法, 具体是一种水文时间序列小波互相关分析方法。 背景技术 小波分析方法 (wavelet analysis, WA) 具有对非平稳时间序列进行时频综合分 析的能力 ( 崔锦泰 . 小波分析导论 . 西安 : 西安交通大学出版社, 1995), 因此适合于研 究具有多时间尺度变化特性的复杂水文水资源系统 ( 王文圣, 丁晶, 李跃清 . 水文小波分 析 . 北京 : 化学工业出版社, 2005 ; Labat D.Recent advances in wavelet analyses : Part 1.A review of concepts.Journal of Hydrology, 2005, 314 : 275-288)。随着理论研究 的深入和解决实际问题的需要, WA 在水文水资源学中的应用日益增多 ( 王文圣, 丁晶, 向 红莲 . 小波分析在水文学中的应用研究及展望 . 水科学进展, 2002, 13(4) : 515-520)。综 合分析可以看出, 目前 WA 主要用于揭示和描述水文系统自身的内部结构和变化特性。而 揭示各水文要素之间相互关系 ( 如降雨与径流序列、 气象因子与降雨和径流序列、 水位与 流量等 ) 是认识水文循环过程和揭示水文演变机制的另一重要途径。传统互相关分析方 法 ( 包括互谱分析方法 ) 由于存在以下主要缺陷 ( 丁晶, 邓育仁 . 随机水文学 . 成都 : 成 都科技大学出版社, 1988) : (1) 仅适用于平稳各态历经序列 ; (2) 无法揭示序列在不同时间 尺度范围内的互相关关系, 使其在实际应用中有较大局限。小波互相关分析方法 (wavelet cross-correlation, WCC) 能够实现对两非平稳时间序列在特定时间尺度和指定时滞下互 相关关系的定量描述, 具有更大的适用性和优越性, 因此可以很好地用于研究和揭示各水 文要素之间的相互关系。然而, 由于针对 WCC 的系统研究偏少, 目前缺乏统一的求解公式和
分析方法, 且本身理论体系还不甚完善, 仅在经济学、 信号处理、 临床医学等方面得到一定 的应用, 而在水文水资源学中的研究和应用非常少, 且国内目前未见有相关报道。
为此, 本发明旨在探讨适合于研究水文水资源学问题的小波互相关分析理论 和方法体系。首先经分析和整理, 系统地介绍用于水文序列分析的基于连续小波变换 (continuous wavelet transform, CWT) 的小波互相关分析方法 ; 同时定义了基于 CWT 的小 波互协方差和小波互相关度两个定量指标, 用于描述两时间序列在整体时间域上的互相关 关系 ; 然后提出绘制小波互相关系数等值线图的方法, 通过该等值线图可达到对两时间序 列互相关关系进行 “时频综合分析” 的目的。最后结合具体实例加以简要分析, 以显示小波 互相关分析方法用于研究水文时间序列的适用性和优越性。
分析实测时间序列 x(t) 和 y(t) 之间的互相关关系时, 实际中常用下式求解互相 关系数。
式中, σx、 σy 分别表示序列 x(t) 和 y(t) 的均方差,4分别表示 x(t) 和 y(t)102033851 A CN 102033854说明书2/7 页的均值, Cxy(k) 表示两时间序列在时滞 k 下的互协方差, rxy(k) 为两时间序列在时滞 ( 也称 时移 )k 下的互相关系数, n 表示时间序列的长度。
实际中, 当需要定量描述各非平稳时间序列之间在某特定时间尺度上的互相关关 系时, 该式无能为力。 发明内容 本发明针对常规互相关分析方法存在的缺陷, 并基于小波分析方法, 提出了一种 水文时间序列小波互相关分析方法, 并绘制了小波互相关系数等值线图, 通过该等值线图 可达到对两时间序列互相关关系进行 “时频综合分析” 的目的。
本发明所述的一种水文时间序列小波互相关分析方法, 其特征在于包括以下步 骤:
(1) 选择小波函数和时间尺度范围, 然后对待分析的水文时间序列进行连续小波 变换分析 (continuous wavelet transform, CWT) ;
(2) 计算水文时间序列在不同时间尺度上及不同时滞下的小波互协方差 ;
(3) 根据小波互协方差计算结果, 求解两序列在不同时间尺度上及不同时滞下的 小波互相关系数 ;
(4) 求得不同时间尺度上及不同时滞下的小波互相关系数之后, 计算小波互相关 度; 以描述两序列在整体时间域上的互相关程度 ;
(5) 求得小波互相关系数和小波互相关度之后, 绘制小波互相关系数等值线图, 并通过详细分析小波互相关系数等值线图, 掌握所研究序列之间由整体到局部的互相关关 系; 实现对时间序列之间互相关关系进行时频综合分析。
上述步骤 (2) 根据连续小波变换分析结果存在有实部和模两种不同的情况, 分别 定义了小波互协方差的求解方法 :
或式中 Wcovxy(a, k) = E[Wx(a, b)Wy(a, b+k)] (3)
上式中, Wx(a, b) 和 Wy(a, b) 分别表示分析序列 x 和 y 时得到的连续小波变换系 数, a 表示时间尺度因子, b 表示时间位置因子, k 表示时滞。R() 和 I() 分别表示求解括弧 内复数的实部和虚部, E() 表示求解均值。WCxy(a, k) 表示在时间尺度 a 和时滞 k 下求得的 小波互协方差。
上述步骤 (3) 根据小波互协方差的两种不同的情况, 分别定义了小波互相关系数 的求解方法 :
或 上式中, WRxy(a, k) 表示在时间尺度 a 和时滞 k 下求得的小波互相关系数, ‖ 表示求解绝对值, 其余公式符号同上。
上述步骤 (4) 根据小波互相关系数求解结果, 定义了小波互相关度的求解方法, 以描述两序列在整体时间域上的互相关程度, 步骤如下 :
(4.1) 在求得两时间序列在尺度 a 和时滞 k 下小波互相关系数 WRxy(a, k) 的基础 上, 通过积分求得两时间序列在时滞 k 下对应整体时间域上的小波互相关程度的总和 :
WRxy(k) =∫ WRxy(a, k)2da (6)
(4.2) 然后, 求解各时间尺度 a 下的小波互相关系数 WRxy(a, k) 的权重系数 : 2
f(WRxy(a, k)) = WRxy(a, k) /WRxy(k) (7)
(4.3) 求解两时间序列在时滞 k 下的小波互相关度 (wavelet cross-correlation degree, WCCD) 为 :
WCCDxy(k) =∫ WRxy(a, k)f(WRxy(a, k))da (8)
步骤 (5) 绘制小波互相关系数等值线图时, 以横轴表示时滞 k 的取值, 纵轴表示时 间尺度 a 的取值, 图中的某点数值表征了对应尺度 a 和时滞 k 下两序列的互相关系数的大 小。
步骤 (5) 绘制小波互相关系数等值线图之后, 对两时间序列之间的互相关关系进 行详细分析, 主要步骤如下 : (5.1) 通过对小波互相关系数等值线图进行垂向截取, 分析在固定时滞下, 两序列 在各时间尺度上互相关程度大小的变化情况 ;
(5.2) 通过对等值线图进行横向截取, 分析在固定时间尺度上, 两序列在各时滞下 互相关程度的变化情况 ;
(5.3) 通过分析各时间尺度上小波互相关系数值的正负性, 掌握两序列在各时间 尺度上互相关性的正负变化情况 ;
(5.4) 通过对比分析各时间尺度上小波互相关系数绝对值的大小, 识别并提取出 对应某个或若干个互相关性明显的时间尺度范围 ;
(5.5) 通过对比分析在各时滞下小波互相关系数值的大小, 识别出两时间序列之 间最显著的时间延迟关系。
本发明提供了一种适合于研究水文水资源学问题的小波互相关分析理论和方法 体系。首先经分析和整理, 系统地提出了用于水文序列分析的基于连续小波变换 (CWT) 的 小波互相关分析方法 ; 同时定义了基于 CWT 的小波互协方差和小波互相关度两个定量指 标, 用于描述两时间序列在整体时间域上的互相关关系 ; 并提出绘制小波互相关系数等值 线图的方法, 通过该等值线图可达到对两时间序列互相关关系进行 “时频综合分析” 的目 的。实例分析结果显示了小波互相关分析方法的有效性和优越性, 应用其可实现由整体到 局部认识时间序列互相关关系的目的。 小波互相关分析方法能够分析和定量描述非平稳时 间序列在特定时间尺度和指定时滞下的互相关关系, 可克服传统互相关分析方法的局限, 具有更好的灵活性和适用性。因此, 通过应用小波互相关分析方法对各水文要素进行分析 和描述, 更有助于深入认识复杂的水文循环过程和揭示水文演变机制, 同时也可促进小波 互相关分析自身理论和方法体系的不断完善。
附图说明 图 1 小波互相关系数等值线图绘制流程。
图 2 利津站和花园口站年径流序列小波互相关系数等值线图 ( 图 2a)、 不同时滞下 ( 图 2b) 和不同时间尺度上 ( 图 2c) 小波互相关系数曲线。
图 3 利津站和花园口站年径流序列小波互相关度曲线和互相关系数曲线。
具体实施方式
以下具体介绍基于连续小波变换的小波互相关分析方法 :
令 L2(R) 表 示 定 义 在 实 轴 上、 可 测 的 平 方 可 积 函 数 空 间, 则对于时间序列 2 x(t) ∈ L (R), 其连续小波变换 (continuous wavelet transform, CWT) 公式为 :
a, b ∈ R, a≠0(9)式中, ψ*(t) 为 ψ(t) 的复共轭函数, Wx(a, b) 为 x(t) 经连续小波变换后的小波 系数。a 为时间尺度因子, b 为时间位置因子。ψa,b(t) 为小波函数, 是由一满足 “容许性” * 条件 ( 式 10) 的小波母函数 ψ(t) 经尺度伸缩和平移后得到。式 (10) 中, ψ (ω) 表示复 * 共轭函数 ψ (t) 在频率 ω 处的 Fourier 变换。
小波变换的实质是采用一种大小可变, 位置可动的变窗对时间序列进行谱分析, 因此可满足序列时频局部化分析的要求。 对水文序列进行小波变换后得到的小波系数是序 列在不同时间尺度和不同时间位置上的投影, 可用来刻画和描述水文序列的组成结构和多 时间尺度变化特性。 此外, 通过研究两时间序列对应小波系数之间的关系, 同样可以刻画和 描述两序列之间的互相关关系。
1、 基于 CWT 的小波互相关系数求解方法
依 据 传 统 时 间 序 列 互 协 方 差 的 定 义 式, 在 此 处 定 义 小 波 互 协 方 差 (wavelet cross-covariance)。由于某些连续小波变换系数结果 ( 例如用 “Morlet” 小波函数进行小 波变换 ) 有实部和模两个重要变量, 实部表示不同特征时间尺度信号在不同时间位置上的 分布和相位两方面的信息, 模的大小主要表示特征时间尺度信号的强弱。 因此, 此处分别定 义对应的小波互协方差 WCxy(a, k) :
或式中 Wcovxy(a, k) = E[Wx(a, b)Wy(a, b+k)] (3)
其中, Wx(a, b) 和 Wy(a, b) 分别表示时间序列 x(t) 和 y(t) 在尺度 a 下对应的连 续小波变换系数, k 表示时滞, R() 和 I() 分别代表括号内变量的实部 (real part) 和虚部 (imaginary part), E[] 表示方括弧内结果的均值, WCxy(a, k) 表示序列 x(t) 和 y(t) 在尺 度 a 和时滞 k 下对应的小波互协方差。
在求得小波互协方差的基础之上, 首先介绍小波局部互相关系数 (wavelet local correlation coefficient, WLCC), 用于分析两序列在特定时间尺度和时间位置点上的互相关关系。
然 后 对 时 间 序 列 x(t) 和 y(t) 进 行 小 波 互 相 关 分 析 (wavelet cross-correlation, WCC)。此处同样针对连续小波变换系数结果的实部和模两个变量, 分 别定义对应的小波互相关系数 WRxy(a, k) :
或WRxy(a, k) 定量描述了时间序列 x(t) 和 y(t) 在时间尺度 a 上和时滞 k 下相应的 互相关程度。
2、 时间序列小波互相关度
WCC 方法主要用于分析两时间序列在特定时间尺度上和指定时滞下的互相关关 系。 为便于分析和描述时间序列之间在整体时间域上的互相关程度, 笔者经分析探讨, 在此 处提出基于 CWT 的时间序列小波互相关度 (wavelet cross-correlation degree, WCCD) 的 概念, 用于刻画时间序列之间整体时间域上的互相关程度。具体如下 : 在求得两时间序列在尺度 a 和时滞 k 下小波互相关系数 WRxy(a, k) 的基础上, 通过 积分求得两时间序列在时滞 k 下对应整体时间域上的小波互相关程度的总和 :
WRxy(k) =∫ WRxy(a, k)2da (6)
然后, 各时间尺度 a 下的小波互相关系数 WRxy(a, k) 的权重系数可定义为 : 2
f(WRxy(a, k)) = WRxy(a, k) /WRxy(k) (7)
定义两时间序列的小波互相关度 (wavelet cross-correlation degree, WCCD) 为:
WCCDxy(k) =∫ WRxy(a, k)f(WRxy(a, k))da (8)
WCCDxy(k) 表征了两时间序列在时滞 k 下整体时间域上互相关程度的大小, 其实质 是在相同时滞 k 下, 对各时间尺度上的互相关程度求解其加权期望值, 因此可综合反映两 时间序列在各时间尺度上关于互相关性的大小和分布两方面信息。
此外, 通过绘制 WCCDxy(k) 随时滞 k 的变化曲线 (WCCDxy(k) ~ k 曲线 ), 可以分析 两序列在整体时间域上随时滞变化时互相关性的变化特点和规律。
3、 小波互相关系数等值线图
类似于通过绘制小波系数等值线图可以分析时间序列自身的结构和多时间尺度 变化特性, 在此处提出绘制小波互相关系数等值线图的方法, 并利用小波互相关系数等值 线图由整体到局部定量分析时间序列之间在各时间尺度和各时滞下的互相关关系, 即达到 对时间序列之间互相关关系进行 “时频综合分析” 的目的。
小波互相关系数等值线图的绘制方法 ( 图 1) 如下 : (1) 首先选择合理的小波函数 和时间尺度范围, 对所研究各时间序列进行连续小波变换, 得到各序列对应的小波系数结 果; (2) 利用式 (4)( 或式 (5) 求解两时间序列在某时间尺度 a 和时滞 k 下对应的小波互相
关系数 ; (3) 依次取不同的 a 和 k 值, 分别对两序列进行小波互相关分析, 最终得到两时间 序列在各尺度和各时滞上对应的小波互相关系数 ; (4) 绘制小波互相关系数等值线图, 其 中横轴表示时滞 k 的取值, 纵轴表示时间尺度 a 的取值, 图中的某点数值表征了对应尺度 a 和时滞 k 下两序列的互相关系数的大小 ; (5) 通过定量分析小波互相关系数等值线图, 掌握 所研究序列之间由整体到局部的互相关关系 ; (6) 通过求解 WCCD, 掌握两序列不同时滞下 在整体时间域上的互相关程度。
小波互相关系数等值线图对描述和刻画两时间序列之间的相互关系具有重要的 应用意义。主要如下 : (1) 通过对小波互相关系数等值线图进行垂向截取, 可以分析在固定 时滞下, 两序列在各时间尺度上互相关程度大小的变化情况 ; (2) 通过对等值线图进行横 向截取, 可以分析在固定时间尺度上, 两序列在各时滞下互相关程度的变化情况 ; (3) 通过 分析各时间尺度上小波互相关系数值的正负性, 可以掌握两序列在各时间尺度上互相关性 的正负变化情况 ; (4) 通过对比分析各时间尺度上小波互相关系数绝对值的大小, 可以识 别并提取出对应某个 ( 或若干个 ) 互相关性明显的时间尺度范围 ; (5) 通过对比分析在各 时滞下小波互相关系数值的大小, 可以识别出两时间序列之间最显著的时间延迟关系等。
实例分析
选择黄河花园口站和利津站两水文站点实测的 54 年 (1950-2003) 年径流序列为 例进行分析。研究两站年径流序列之间的相互关系, 对认识黄河下游的径流过程和水文变 化特性具有重要的意义。
应用前述的基于 CWT 的小波互相关分析方法, 深入分析在不同时间尺度上以及不 同时滞下两年径流序列之间的互相关关系。此处选用式 (4) 求解小波互相关系数 WRxy(a, k)。结合文献 [ 桑燕芳, 王栋 . 连续小波变换在黄河河口地区特性分析中的应用研究 [A]. 第五届中国水论坛论文集 [C]. 北京 : 中国水利水电出版社, 2007, 766-770.] 中的相关分析 结果, 此处同样选用 “Morlet” 小波函数对花园口 54 年年径流序列进行连续小波变换。然 后依据 CWT 分析结果, 分别求解各时间尺度和各时滞下的小波互相关系数, 并绘制小波互 相关系数等值线图 ( 图 2a)。其中, 最大时间尺度取 50a, 最大时滞取 20a。为便于分析说 明, 分别绘制时滞 k = 0、 2 和 5 时的 WRxy(a, k) ~ a 曲线 ( 图 2b), 和时间尺度 a = 3、 7、 11 和 20 时的 WRxy(a, k) ~ k 曲线 ( 图 2c)。最后绘制小波互相关度 WCCDxy(k) ~ k 曲线和互 相关系数曲线 ( 图 3)。传统互相关分析方法计算结果显示 ( 图 3) : 两年径流序列之间的最 大互相关系数为 rxy(0) = 0.9655。
综合上述两径流序列互相关关系的分析结果, 可得到如下主要结论 :
(1) 对小波互相关系数等值线图进行分析, 可以识别出两径流序列主要在四个时 间尺度上存在较为明显的互相关关系 : a = 3、 7、 11 和 20。此四个时间尺度对应的也正好是 各径流序列自身四个明显的周期值。
(2) 对小波互相关系数等值线图进行垂向截取, 可以了解两序列在固定时滞下各 时间尺度上的互相关关系。图 2b 显示, ①在同一时滞下, 不同时间尺度上两序列互相关性 程度的大小不同 ; ②在不同时间尺度范围内两序列互相关的正负特性也不相同。以 k = 5 为例, 在时间尺度 4a、 8a、 和 50a 上的小波互相关系数分别为 0.58、 -0.95 和 0.88 ; ③随着 时滞 k 由 0 增大到 5, 在时间尺度 5-50a 范围内两序列互相关关系有逐渐减弱的趋势, 而在 1-5a 时间尺度范围内时滞 5 时的互相关关系要较时滞 2 的互相关关系更明显。(3) 对小波互相关系数等值线图进行横向截取, 可以了解两序列在固定时间尺度 上各时滞下的互相关关系。图 2c 显示, ①在相同时间尺度上, 两径流序列互相关性的大小 随时滞变化时存在波动, 且存在相关性逐渐减弱的趋势。以 a = 7 为例, 两序列在时滞 k = 0、 4、 9、 14 和 19 时, 互相关系数的波动变化为 0.94、 -0.88、 0.78、 -0.59、 0.58 ; ②随着时间 尺度的逐渐增大, 两径流序列互相关系数随时滞变化的波动周期也逐渐增大。
(4) 小波互相关度分析结果图 3 显示, 两序列整体时间域上在不同时滞下的互相 关程度不同。小波互相关度在时滞 k = 0 时取得最大值 0.96, 与传统互相关分析方法的结 果一致。之后随时滞增大时小波互相关度逐渐减小, 在 k = 13 时小波互相关度又达到另一 极大值 0.80。 与互相关系数曲线对比可以看出, 小波互相关度曲线更加光滑、 波动规律更加 明显, 这主要是由于小波互相关度曲线是上述各时间尺度上互相关关系的加权平均值, 因 此可以更有效地描述两序列在整体时间域上随时滞推移时的互相关关系变化情况。 而传统 互相关分析方法无法考虑和定量描述不同时间尺度上序列互相关关系之间的差异, 且曲线 中得到的多个峰值点缺乏实际的物理意义。
(5) 综合图 2 和图 3 的分析结果可以看出 : ①利津站和花园口站两个实测年径流 序列在时滞为 0 时的互相关关系最明显, 这与径流由花园口站至利津站的实际时间相符 合; ②两年径流序列在时间尺度为 3a、 7a、 11a 和 20a 时的互相关关系最明显, 即表明在序列 自身的四个明显周期变化时间尺度范围内, 两序列的互相关关系也最明显。 (6) 由实例分析 可以看出, 小波互相关分析方法可以揭示和描述非平稳时间序列在不同时间尺度上和不同 时滞下的互相关关系, 因此可克服传统互相关分析方法的局限, 有助于对时间序列之间互 相关关系进行全面细致的定量分析。