CN105069309A - 一种识别水文时间序列非线性趋势的方法 - Google Patents

Abstract

本发明公开了一种识别水文时间序列非线性趋势的方法，包括：根据序列长度计算最大小波分解水平，确定具体的离散小波变换方法；得到不同分解水平上对应的子序列；计算各子序列的小波能量密度值，得到待分析水文时间序列的小波能量密度函数；利用离散小波变换方法对白噪声序列进行分解得到子序列；将各白噪声序列小波能量密度函数的均值作为标准小波能量密度函数，得到标准小波能量密度函数的置信区间；对比最大时间尺度上待分析水文序列子序列的小波能量密度值与标准小波能量密度函数置信区间的位置关系。本发明解决了小波分析方法在水文时间序列趋势识别方面缺乏可靠的水文物理基础，也无法有效估计水文序列非线性趋势的显著性和不确定性的问题。

Description

一种识别水文时间序列非线性趋势的方法

技术领域

本发明涉及水文科学技术领域，尤其指代一种识别水文时间序列非线性趋势的方法。

背景技术

水文时间序列分析是揭示和认识自然界水循环过程变化特性的重要手段和技术途径。实际水文时间序列分析过程中，趋势识别与提取是一项十分重要的内容，其主要目的是揭示水文变量在大时间尺度上的变化规律。此外，在水文时间序列相关性分析以及频谱分析过程中，也需要首先去除水文序列中的趋势项，防止序列相关性和频谱分析结果受到非零均值或趋势的影响。尽管目前关于水文时间序列趋势识别已有大量相关研究，但准确识别水文序列的趋势仍是一项较为困难的工作。

目前水文时间序列趋势识别方法大致可以分为四类：第一类是基于数据拟合的趋势识别方法，使用这类方法时一般需要根据个人经验事先给定先验函数去拟合趋势，由于受序列长度较短等因素制约，由最小二乘法或极大似然法得到的参数估计值往往具有不确定性，由此得到的趋势识别结果带有主观性且不可靠。第二类是基于时间域分析的趋势识别方法，这类方法中最常用的是Mann-Kendall(MK)趋势检验方法，其简单易操作且不会受序列缺值等不利因素的影响，但会受到序列相关性、序列长度、待识别趋势项强度等因素的影响；为克服MK趋势检验方法的缺陷，许多学者也提出了大量用于去除序列相关性的方法；Spearman秩次相关性检验方法类似于MK方法，但在实际中较少使用；线性回归是另外一类基于时间域的趋势识别方法，但由于水文时间序列常表现出不同时间尺度上的非线性和非平稳特性，因此线性趋势识别结果缺乏真实的物理依据且往往不合理。第三类是基于频率域分析的趋势识别方法，其中具有代表性的滑动平均方法也需要事先指定一个时间尺度，然而这个时间尺度是一个未知的先验信息，其它一些更加复杂的频率域趋势识别方法(例如傅里叶变换方法等)由于均基于平稳性和线性假设，因此也缺乏可靠的水文物理基础。相比较而言，基于时频域综合分析的第四类趋势识别方法的性能更优，其中具有代表性的是基于小波分析的趋势识别方法，因为该方法能够同时揭示序列在时域和频域内的非平稳变化特性，目前小波分析方法已广泛应用于识别水文时间序列的趋势变化。此外，在其他大量研究过程中，也常将小波分析方法与前面的三类方法联合使用识别水文时间序列的趋势项。

综合来看，大量实例分析验证了小波分析方法在水文时间序列趋势识别方面较常规方法具有很大的优势，然而实际小波分析结果会受到小波函数选择和分解水平选择等不利因素的影响，且许多基于小波分析的趋势识别方法尽管具有很好的数据理论，但缺乏可靠的水文物理基础，也无法有效的估计水文序列非线性趋势的显著性和不确定性。文献[SangY.F.,WangZ.G.,LiuC.M.,2013.Discretewavelet-basedtrendidentificationinhydrologictimeseries.HydrologicalProcesses,DOI:10.1002/hyp.9356](以下简称为“文献[1]”)中利用小波分析方法尝试判别水文序列趋势的显著性，但其能量曲线的计算公式存在较大缺陷，容易低估趋势结果的显著性。整体上，有效的水文时间序列趋势识别方法应该能够准确分离序列的非线性趋势，此外还应能够定量估计非线性趋势在统计意义上的显著性。

发明内容

针对于上述问题，本发明的目的在于提供一种识别水文时间序列非线性趋势的方法，以解决现有技术中小波分析方法在水文时间序列趋势识别方面缺乏可靠的水文物理基础，也无法有效的估计水文序列非线性趋势的显著性和不确定性的问题。

为达到上述目的，本发明的一种识别水文时间序列非线性趋势的方法，包括步骤如下：

1)检查待分析水文序列数据的一致性和可靠性，选择合理的小波函数与边界点处理方法，根据序列长度计算最大小波分解水平，确定具体的离散小波变换方法；

2)应用所确定的离散小波变换方法对水文时间序列进行分解，得到不同分解水平上对应的子序列，序列f(t)的分解结果记为：

$f\left(t\right)=\underset{i=1}{\overset{N}{\Σ}}{f}_i\left(t\right)+T_N$

其中，N表示最大小波分解水平，f_i(t)表示由高频小波系数重构得到的第i个子序列，T_N是最大分解水平上由低频小波系数重构得到的子序列，一般对应着序列趋势项；

3)计算各子序列的小波能量密度值，得到待分析水文时间序列的小波能量密度函数：

$S\left(i\right)=\frac{1}{n}\underset{t=1}{\overset{n}{\Σ}}{\left(f_i\right(t\left)\right)}^2$

其中，S(i)表示分解水平i上子序列f_i(t)的小波能量密度值，n表示序列长度；

4)利用Monte-Carlo方法生成与待分析水文序列相同长度的白噪声序列，利用离散小波变换方法对白噪声序列进行分解得到子序列，并计算对应的小波能量密度函数；

5)重复上述步骤4)，生成大量白噪声序列并分别计算其小波能量密度函数，直至白噪声序列小波能量密度函数的统计特性稳定；

6)将各白噪声序列小波能量密度函数的均值作为标准小波能量密度函数；通过计算各分解水平上白噪声序列小波能量密度值的95％置信区间，得到标准小波能量密度函数的置信区间；

7)对比最大时间尺度上待分析水文序列子序列T_N的小波能量密度值与标准小波能量密度函数置信区间的位置关系；若位于置信区间外，则表明该序列的非线性趋势在统计意义上显著，若位于置信区间内，则认为该序列的非线性趋势在统计意义上不显著。

进一步地，所述的步骤2)具体包括：

21)对于长度为n的水文时间序列，计算得到最大分解水平：

N＝[log₂(n)]；

22)利用二进制离散小波变换方法(DyadicDiscreteWaveletTransform)对该序列进行分析：

$W_f(i,k)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}f\left(t\right){{\mathrm{\ψ}}_{i,k}}^{*}\left(t\right)dtwith{\mathrm{\ψ}}_{i,k}\left(t\right)=2^{-i/2}\mathrm{\ψ}(2^{-i}t-k);$

其中，i表示分解水平，k表示时间位置因子；ψ^*(t)是小波函数ψ(t)的复共轭函数；W_f(i,k)是离散小波函数；

23)重构分解水平i上的子序列：

f_i(t)＝Σ_kW_f(i,k)ψ^*(2^-it-k)；

24)对不同水解水平上的子序列相加，得到原序列：

$f\left(t\right)=\underset{i=1}{\overset{N}{\Σ}}{f}_i\left(t\right)+T_N.$

本发明的有益效果：

(1)本发明可以自适应性的识别并提取出序列中的非线性趋势，而传统的MK检验等方法仅能得到线性趋势；

(2)本发明可以定量估计非线性趋势识别结果在统计意义上的显著性，而传统方法无法对非线性趋势的显著性进行定量判断；

(3)本发明不仅可以准确识别水文时间序列中的非线性趋势，还可以在考虑不确定性的基础上同时识别出水文序列中的周期等其他确定成分，进而为准确认识水文过程的确定性变化规律奠定良好基础。

附图说明

图1绘示本发明识别水文时间序列非线性趋势的方法的流程图。

图2绘示北半球89年月气温序列的小波能量密度函数。

图3绘示北半球89年月气温序列的非线性趋势识别结果。

具体实施方式

为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。

参照图1至图3所示，本发明的一种识别水文时间序列非线性趋势的方法，于实施例中，具体实施过程如下：

1.DWT离散小波变换方法

实测水文时间序列常是离散信号。令L²(R)表示定义在实轴上、可测的平方可积函数空间，信号f(t)∈L²(R)的离散小波变换(DiscreteWaveletTransform,DWT)可表示为：

$W_f(i,k)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}f\left(t\right){{\mathrm{\ψ}}_{i,k}}^{*}\left(t\right)dtwith{\mathrm{\ψ}}_{i,k}\left(t\right)={a_0}^{-i/2}\mathrm{\ψ}({a_0}^{-i}t-{\mathrm{kb}}_0)$

式中，a₀和b₀均为常数，i表示分解水平(DecompositionLevel,DL；也称时间尺度水平)，k为时间位置因子，可反映时间上的平移；W_f(i,k)是时间序列f(t)通过单位脉冲响应的滤波器输出，能同时反映时域参数b(或k)和频域参数a(或j)的特性。a较小时，信号分析在频域内的分辨率低，但在时域内的分辨率高；a增大时，信号分析在频域内的分辨率增高，但时域内的分辨率降低。因此，小波变换能满足窗口大小和形状可变的信号时频局部化分析的要求。

实际水文时间序列分析过程中，常选择使用二进制离散小波变换(DyadicDiscreteWaveletTransform)对其进行分析，即设定a₀＝2和b₀＝1，表达式如下：

$W_f(i,k)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}f\left(t\right){{\mathrm{\ψ}}_{i,k}}^{*}\left(t\right)dtwith{\mathrm{\ψ}}_{i,k}\left(t\right)=2^{-i/2}\mathrm{\ψ}(2^{-i}t-k)$

根据二进制离散小波变换基本原理，分解水平的理论最大值N可由如下式(1)求得：

N＝[log₂(n)](1)

其中，n表示序列f(t)的长度。

若使用的小波函数满足“规则性条件”式(2)，利用离散小波系数并通过小波逆变换可得到序列在不同分解水平上的子序列式(3)：

$\begin{array}{cc}{\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}{t}^k\mathrm{\ψ}\left(t\right)dt=0,& k=1,...,N-1\end{array}---\left(2\right)$

f_i(t)＝Σ_kW_f(i,k)ψ^*(2^-it-k)(3)

对不同分解水平上的子序列相加，即可重构得到原序列。

2.非线性趋势识别的方法

本发明的识别水文时间序列非线性趋势的方法，通过联合使用DWT方法和Monte-Carlo方法。其中DWT方法主要用于识别并分离出水文时间序列的趋势项，而Monte-Carlo方法主要用于定量判别非线性趋势在统计意义上的显著性；最后待分析水文时间序列的非线性趋势可被准确识别。具体参照图1，具体步骤包括如下：

1)检查待分析水文时间序列数据的一致性和可靠性，选择合理的小波函数与边界点处理方法，根据序列长度计算最大小波分解水平，确定具体的离散小波变换分析方法(DWT)；

2)应用所确定的离散小波变换方法对水文时间序列进行分解，得到不同分解水平上对应的子序列，序列f(t)的分解结果记为：

$f\left(t\right)=\underset{i=1}{\overset{N}{\Σ}}{f}_i\left(t\right)+T_N$

其中，N表示最大小波分解水平，f_i(t)表示由高频小波系数重构得到的第i个子序列，T_N是最大分解水平上由低频小波系数重构得到的子序列，一般对应着序列趋势项；

3)计算各子序列的小波能量密度值，得到待分析水文时间序列的小波能量密度函数：

$S\left(i\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{\left(f_i\right(t\left)\right)}^2$

其中，S(i)表示分解水平i上子序列f_i(t)的小波能量密度值，n表示序列长度；

4)利用Monte-Carlo方法生成与待分析水文序列相同长度的白噪声序列，利用离散小波变换方法对白噪声序列进行分解得到子序列，并计算对应的小波能量密度函数；

5)重复上述步骤4)，生成大量白噪声序列并分别计算其小波能量密度函数，直至白噪声序列小波能量密度函数的统计特性稳定；

6)将各白噪声序列小波能量密度函数的均值作为标准小波能量密度函数；通过计算各分解水平上白噪声序列小波能量密度值的95％置信区间，得到标准小波能量密度函数的置信区间；

7)对比最大时间尺度上待分析水文序列子序列T_N的小波能量密度值与标准小波能量密度函数置信区间的位置关系；若位于置信区间外，则表明该序列的非线性趋势在统计意义上显著，若位于置信区间内，则认为该序列的非线性趋势在统计意义上不显著。

3.算例分析

将上述本发明的方法应用于1901-1989年北半球月气温数据进行分析，参照图2-图3所示。利用MK方法对该序列进行趋势识别时发现，统计值高达17.4，表明该序列有明显的上升趋势。选择最大时间尺度为序列长度，计算得到最大分解水平为10。利用本发明的方法首先得到该序列的小波能量密度函数，由图2可以看出，该序列在分解水平2之后各子序列的小波能量密度明显高于95％置信区间，因此认为他们是该序列中的确定性成分，反映了该序列确定性的变化规律；其中分解水平9和10上子序列之和为该序列识别出的非线性趋势，且该趋势在统计意义上是显著的；由图3可以看出该结果与文献[1]中的结果存在一定差别，主要反映出小波能量密度函数计算方法不同导致的趋势识别结果之间存在差异。为进一步分析验证所提方法分析结果的可靠性，对原序列和趋势识别结果进行对比，可以看出相比于文献[1]中的结果，由该方法得到的趋势识别结果对应更大的互相关系数值(0.922>0.901)和更小的均方根误差值(0.010<0.013)，进一步验证了本发明方法的有效性。

本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims (2)

1.一种识别水文时间序列非线性趋势的方法，其特征在于，包括步骤如下：

1)检查待分析水文序列数据的一致性和可靠性，选择合理的小波函数与边界点处理方法，根据序列长度计算最大小波分解水平，确定具体的离散小波变换方法；

2)应用所确定的离散小波变换方法对水文时间序列进行分解，得到不同分解水平上对应的子序列，序列f(t)的分解结果记为：

$f\left(t\right)=\underset{i=1}{\overset{N}{\&Sigma;}}{f}_i\left(t\right)+T_N$

其中，N表示最大小波分解水平，f_i(t)表示由高频小波系数重构得到的第i个子序列，T_N是最大分解水平上由低频小波系数重构得到的子序列，一般对应着序列趋势项；

3)计算各子序列的小波能量密度值，得到待分析水文时间序列的小波能量密度函数：

$S\left(i\right)=\frac{1}{n}\underset{t=1}{\overset{n}{\&Sigma;}}{\left(f_i\right(t\left)\right)}^2$

其中，S(i)表示分解水平i上子序列f_i(t)的小波能量密度值，n表示序列长度；

4)利用Monte-Carlo方法生成与待分析水文序列相同长度的白噪声序列，利用离散小波变换方法对白噪声序列进行分解得到子序列，并计算对应的小波能量密度函数；

5)重复上述步骤4)，生成大量白噪声序列并分别计算其小波能量密度函数，直至白噪声序列小波能量密度函数的统计特性稳定；

6)将各白噪声序列小波能量密度函数的均值作为标准小波能量密度函数；通过计算各分解水平上白噪声序列小波能量密度值的95％置信区间，得到标准小波能量密度函数的置信区间；

7)对比最大时间尺度上待分析水文序列子序列T_N的小波能量密度值与标准小波能量密度函数置信区间的位置关系；若位于置信区间外，则表明该序列的非线性趋势在统计意义上显著，若位于置信区间内，则认为该序列的非线性趋势在统计意义上不显著。

2.根据权利要求1所述的识别水文时间序列非线性趋势的方法，其特征在于，所述的步骤2)具体包括：

21)对于长度为n的水文时间序列，计算得到最大分解水平：

N＝[log₂(n)]；

22)利用二进制离散小波变换方法对该序列进行分析：

$W_f(i,k)={\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}f\left(t\right){{\mathrm{\&psi;}}_{i,k}}^{*}\left(t\right)dtwith{\mathrm{\&psi;}}_{i,k}\left(t\right)=2^{-i/2}\mathrm{\&psi;}(2^{-i}t-k);$

其中，i表示分解水平，k表示时间位置因子；ψ^*(t)是小波函数ψ(t)的复共轭函数；W_f(i,k)是离散小波函数；

23)重构分解水平i上的子序列：

f_i(t)＝Σ_kW_f(i,k)ψ^*(2^-it-k)；

24)对不同水解水平上的子序列相加，得到原序列：

$f\left(t\right)=\underset{i=1}{\overset{N}{\&Sigma;}}{f}_i\left(t\right)+T_N.$