一种水文时间序列小波互相关分析方法

摘要

本发明公开了一种水文时间序列小波互相关分析方法，首先选择小波函数和时间尺度范围，然后对待分析的水文时间序列进行连续小波变换分析；计算水文时间序列在不同时间尺度上及不同时滞下的小波互协方差；求解两序列在不同时间尺度上及不同时滞下的小波互相关系数，计算小波互相关度以描述两序列在整体时间域上的互相关程度，绘制小波互相关系数等值线图，实现对时间序列之间互相关关系进行时频综合分析。实例分析结果显示了小波互相关分析方法的有效性和优越性，该方法能够分析和定量描述非平稳时间序列在特定时间尺度和指定时滞下的互相关关系，可克服传统互相关分析方法的局限，具有更好的灵活性和适用性。

说明书

技术领域

本发明涉及一种水文时间序列分析方法，具体是一种水文时间序列小波互相关分析方法。

背景技术

小波分析方法(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)之间的互相关关系时，实际中常用下式求解互相关系数。

$r_{\mathrm{xy}}\left(k\right)=C_{\mathrm{xy}}\left(k\right)/{\sigma }_x{\sigma }_y=\frac{\underset{t=1}{\overset{n-k}{\Sigma }}(x\left(t\right)-\bar{x})(y(t+k)-\bar{y})}{{\left[\underset{t=1}{\overset{n}{\Sigma }}{(x\left(t\right)-\bar{x})}^2\underset{t=1}{\overset{n}{\Sigma }}{(y\left(t\right)-\bar{y})}^2\right]}^{1/2}}$

式中，σ_x、σ_y分别表示序列x(t)和y(t)的均方差，分别表示x(t)和y(t)的均值，C_xy(k)表示两时间序列在时滞k下的互协方差，r_xy(k)为两时间序列在时滞(也称时移)k下的互相关系数，n表示时间序列的长度。

实际中，当需要定量描述各非平稳时间序列之间在某特定时间尺度上的互相关关系时，该式无能为力。

发明内容

本发明针对常规互相关分析方法存在的缺陷，并基于小波分析方法，提出了一种水文时间序列小波互相关分析方法，并绘制了小波互相关系数等值线图，通过该等值线图可达到对两时间序列互相关关系进行“时频综合分析”的目的。

本发明所述的一种水文时间序列小波互相关分析方法，其特征在于包括以下步骤：

(1)选择小波函数和时间尺度范围，然后对待分析的水文时间序列进行连续小波变换分析(continuous wavelet transform，CWT)；

(2)计算水文时间序列在不同时间尺度上及不同时滞下的小波互协方差；

(3)根据小波互协方差计算结果，求解两序列在不同时间尺度上及不同时滞下的小波互相关系数；

(4)求得不同时间尺度上及不同时滞下的小波互相关系数之后，计算小波互相关度；以描述两序列在整体时间域上的互相关程度；

(5)求得小波互相关系数和小波互相关度之后，绘制小波互相关系数等值线图，并通过详细分析小波互相关系数等值线图，掌握所研究序列之间由整体到局部的互相关关系；实现对时间序列之间互相关关系进行时频综合分析。

上述步骤(2)根据连续小波变换分析结果存在有实部和模两种不同的情况，分别定义了小波互协方差的求解方法：

${\mathrm{WC}}_{\mathrm{xy}}(a,k)=\sqrt{R{\left({\mathrm{Wcov}}_{\mathrm{xy}}(a,k)\right)}^2}---\left(1\right)$

或${\mathrm{WC}}_{\mathrm{xy}}(a,k)=\sqrt{R{\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}^2+I{\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}^2}---\left(2\right)$

式中   Wcov_xy(a，k)＝E[W_x(a，b)W_y(a，b+k)]                    (3)

上式中，W_x(a，b)和W_y(a，b)分别表示分析序列x和y时得到的连续小波变换系数，a表示时间尺度因子，b表示时间位置因子，k表示时滞。R()和I()分别表示求解括弧内复数的实部和虚部，E()表示求解均值。WC_xy(a，k)表示在时间尺度a和时滞k下求得的小波互协方差。

上述步骤(3)根据小波互协方差的两种不同的情况，分别定义了小波互相关系数的求解方法：

${\mathrm{WR}}_{\mathrm{xy}}(a,k)=\frac{R\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}{\sqrt{R\left(W{\mathrm{cov}}_{\mathrm{xx}}(a,0)\right)R\left(W{\mathrm{cov}}_{\mathrm{yy}}(a,0)\right)}}---\left(4\right)$

或${\mathrm{WR}}_{\mathrm{xy}}(a,k)=\frac{\sqrt{R{\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}^2+I{\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}^2}}{\sqrt{\left|W{\mathrm{cov}}_{\mathrm{xx}}(a,0)\right|\left|W{\mathrm{cov}}_{\mathrm{yy}}(a,0)\right|}}---\left(5\right)$

上式中，WR_xy(a，k)表示在时间尺度a和时滞k下求得的小波互相关系数，‖表示求解绝对值，其余公式符号同上。

上述步骤(4)根据小波互相关系数求解结果，定义了小波互相关度的求解方法，以描述两序列在整体时间域上的互相关程度，步骤如下：

(4.1)在求得两时间序列在尺度a和时滞k下小波互相关系数WR_xy(a，k)的基础上，通过积分求得两时间序列在时滞k下对应整体时间域上的小波互相关程度的总和：

WR_xy(k)＝∫WR_xy(a，k)²da                                     (6)

(4.2)然后，求解各时间尺度a下的小波互相关系数WR_xy(a，k)的权重系数：

f(WR_xy(a，k))＝WR_xy(a，k)²/WR_xy(k)                           (7)

(4.3)求解两时间序列在时滞k下的小波互相关度(wavelet cross-correlation degree，WCCD)为：

WCCD_xy(k)＝∫WR_xy(a，k)f(WR_xy(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利津站和花园口站年径流序列小波互相关度曲线和互相关系数曲线。

具体实施方式

以下具体介绍基于连续小波变换的小波互相关分析方法：

令L²(R)表示定义在实轴上、可测的平方可积函数空间，则对于时间序列x(t)∈L²(R)，其连续小波变换(continuous wavelet transform，CWT)公式为：

$W_x(a,b)={\int }_{-\infty }^{+\infty }x\left(t\right){{\psi }^{*}}_{a,b}\left(t\right)\mathrm{dt}$${\psi }_{a,b}\left(t\right)=|a{|}^{-1/2}\psi \left(\frac{t-b}{a}\right)$a，b∈R，a≠0                (9)

式中，ψ^*(t)为ψ(t)的复共轭函数，W_x(a，b)为x(t)经连续小波变换后的小波系数。a为时间尺度因子，b为时间位置因子。ψ_a，b(t)为小波函数，是由一满足“容许性”条件(式10)的小波母函数ψ(t)经尺度伸缩和平移后得到。式(10)中，ψ^*(ω)表示复共轭函数ψ^*(t)在频率ω处的Fourier变换。

$C_{\psi }={\int }_{-\infty }^{+\infty }\frac{{\left|{\psi }^{*}\left(\omega \right)\right|}^2}{\left|\omega \right|}\mathrm{d\omega }<\infty ---\left(10\right)$

小波变换的实质是采用一种大小可变，位置可动的变窗对时间序列进行谱分析，因此可满足序列时频局部化分析的要求。对水文序列进行小波变换后得到的小波系数是序列在不同时间尺度和不同时间位置上的投影，可用来刻画和描述水文序列的组成结构和多时间尺度变化特性。此外，通过研究两时间序列对应小波系数之间的关系，同样可以刻画和描述两序列之间的互相关关系。

1、基于CWT的小波互相关系数求解方法

依据传统时间序列互协方差的定义式，在此处定义小波互协方差(wavelet cross-covariance)。由于某些连续小波变换系数结果(例如用“Morlet”小波函数进行小波变换)有实部和模两个重要变量，实部表示不同特征时间尺度信号在不同时间位置上的分布和相位两方面的信息，模的大小主要表示特征时间尺度信号的强弱。因此，此处分别定义对应的小波互协方差WC_xy(a，k)：

${\mathrm{WC}}_{\mathrm{xy}}(a,k)=\sqrt{R{\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}^2}---\left(1\right)$

或${\mathrm{WC}}_{\mathrm{xy}}(a,k)=\sqrt{R{\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}^2+I{\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}^2}---\left(2\right)$

式中      Wcov_xy(a，k)＝E[W_x(a，b)W_y(a，b+k)]              (3)

其中，W_x(a，b)和W_y(a，b)分别表示时间序列x(t)和y(t)在尺度a下对应的连续小波变换系数，k表示时滞，R()和I()分别代表括号内变量的实部(real part)和虚部(imaginary part)，E[]表示方括弧内结果的均值，WC_xy(a，k)表示序列x(t)和y(t)在尺度a和时滞k下对应的小波互协方差。

在求得小波互协方差的基础之上，首先介绍小波局部互相关系数(wavelet local correlation coefficient，WLCC)，用于分析两序列在特定时间尺度和时间位置点上的互相关关系。

${\mathrm{WLCC}}_{\mathrm{xy}}(a,b)=\frac{R\left(W_{\mathrm{xy}}(a,b)\right)}{\left|W_x(a,b)\right|\left|W_y(a,b)\right|}=\frac{R\left({W_x}^{*}(a,b)W_y(a,b)\right)}{\left|W_x(a,b)\right|\left|W_y(a,b)\right|}---\left(11\right)$

然后对时间序列x(t)和y(t)进行小波互相关分析(wavelet cross-correlation，WCC)。此处同样针对连续小波变换系数结果的实部和模两个变量，分别定义对应的小波互相关系数WR_xy(a，k)：

${\mathrm{WR}}_{\mathrm{xy}}(a,k)=\frac{R\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}{\sqrt{R\left(W{\mathrm{cov}}_{\mathrm{xx}}(a,0)\right)R\left({\mathrm{Wcov}}_{\mathrm{yy}}(a,0)\right)}}---\left(4\right)$

或${\mathrm{WR}}_{\mathrm{xy}}(a,k)=\frac{\sqrt{R{\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}^2+I{\left(W{\mathrm{cov}}_{\mathrm{xy}}(a,k)\right)}^2}}{\sqrt{\left|W{\mathrm{cov}}_{\mathrm{xx}}(a,0)\right|\left|{\mathrm{Wcov}}_{\mathrm{yy}}(a,0)\right|}}---\left(5\right)$

WR_xy(a，k)定量描述了时间序列x(t)和y(t)在时间尺度a上和时滞k下相应的互相关程度。

2、时间序列小波互相关度

WCC方法主要用于分析两时间序列在特定时间尺度上和指定时滞下的互相关关系。为便于分析和描述时间序列之间在整体时间域上的互相关程度，笔者经分析探讨，在此处提出基于CWT的时间序列小波互相关度(wavelet cross-correlation degree，WCCD)的概念，用于刻画时间序列之间整体时间域上的互相关程度。具体如下：

在求得两时间序列在尺度a和时滞k下小波互相关系数WR_xy(a，k)的基础上，通过积分求得两时间序列在时滞k下对应整体时间域上的小波互相关程度的总和：

WR_xy(k)＝∫WR_xy(a，k)²da                        (6)

然后，各时间尺度a下的小波互相关系数WR_xy(a，k)的权重系数可定义为：

f(WR_xy(a，k))＝WR_xy(a，k)²/WR_xy(k)              (7)

定义两时间序列的小波互相关度(wavelet cross-correlation degree，WCCD)为：

WCCD_xy(k)＝∫WR_xy(a，k)f(WR_xy(a，k))da           (8)

WCCD_xy(k)表征了两时间序列在时滞k下整体时间域上互相关程度的大小，其实质是在相同时滞k下，对各时间尺度上的互相关程度求解其加权期望值，因此可综合反映两时间序列在各时间尺度上关于互相关性的大小和分布两方面信息。

此外，通过绘制WCCD_xy(k)随时滞k的变化曲线(WCCD_xy(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)求解小波互相关系数WR_xy(a，k)。结合文献[桑燕芳，王栋.连续小波变换在黄河河口地区特性分析中的应用研究[A].第五届中国水论坛论文集[C].北京：中国水利水电出版社，2007，766-770.]中的相关分析结果，此处同样选用“Morlet”小波函数对花园口54年年径流序列进行连续小波变换。然后依据CWT分析结果，分别求解各时间尺度和各时滞下的小波互相关系数，并绘制小波互相关系数等值线图(图2a)。其中，最大时间尺度取50a，最大时滞取20a。为便于分析说明，分别绘制时滞k＝0、2和5时的WR_xy(a，k)～a曲线(图2b)，和时间尺度a＝3、7、11和20时的WR_xy(a，k)～k曲线(图2c)。最后绘制小波互相关度WCCD_xy(k)～k曲线和互相关系数曲线(图3)。传统互相关分析方法计算结果显示(图3)：两年径流序列之间的最大互相关系数为r_xy(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)由实例分析可以看出，小波互相关分析方法可以揭示和描述非平稳时间序列在不同时间尺度上和不同时滞下的互相关关系，因此可克服传统互相关分析方法的局限，有助于对时间序列之间互相关关系进行全面细致的定量分析。