CN105205217B - 一种判别水文时间序列非平稳性的方法 - Google Patents

Abstract

本发明公开一种判别水文时间序列非平稳性的方法，包括：确定具体的集合经验模态分解方法；对序列进行分解；利用白噪声的能量扩散函数，识别出原序列中所有的周期成分；去除原序中的周期成分，将剩余成分作为新序列；选择合适的模型方程进行拟合；利用选择的模型方程对新序列进行单位根检验；利用自相关系数图和偏相关系数图，进一步分析新序列的统计特性；比较自（偏）相关系数结果与单位根检验结果的一致性；若所述的新序列表现出良好的相关性，且单位根检验特征统计量接受原假设，则表明序列具有非平稳性；若二者的结果存在不一致，则认为序列不具有非平稳性。本发明克服了序列周期项对单位根检验结果的影响，得到准确的序列非平稳性检验结果。

Description

一种判别水文时间序列非平稳性的方法

技术领域

本发明属于水文科学技术领域，尤其指代一种判别水文时间序列非平稳性的方法。

背景技术

揭示自然界水循环过程的变化特性和演变规律是解决复杂交叉水问题的重要基础，因此一直是水文学领域的核心研究课题。实际中，欲准确认识水循环过程的复杂变化特性是十分困难的科学问题。特别是近几十年来，受气候变化和人类活动的持续影响，许多流域的天然水循环过程受到破坏，各主要水循环要素的周期和趋势等特性发生变异，降雨——径流响应关系明显改变，洪涝灾害等极端水文事件频发，且未来旱涝灾害发生的不确定性和风险会愈来愈大。大量水文观测结果显示，许多流域实测水文资料呈现出明显的非平稳性，给水利工程设计、水文模拟计算、水资源规划等一系列涉水工作带来很大困难。

依据数理统计学理论，数据非平稳性主要是指大样本的统计推断基础——数据一致性要求遭到破坏。自2008年Milly等人在Science发表论文《Stationarity Is Dead:Whither Water Management？》以来，水文过程非平稳性正逐渐成为国际水文学界的研究热点，尤其是全球环境变化研究领域重点关注的难点问题之一。判别水文过程是否具有平稳性有十分重要的意义。水文过程平稳性可以简化水文分析计算工作的难度，提高序列统计特征量的估计精度。从随机水文学角度看，常用的水文时间序列模型是以通过揭示时间序列自身变化规律为主线发展起来的方法论。由于经典回归分析方法的一个重要假设是数据的平稳性，因此识别水文过程是否具有平稳性也是开展回归分析和水文时间序列模拟预报的前提。如果水文时间序列表现出明显的非平稳性，经典回归分析的基础和有效性已不适用，用样本均值和方差等统计指标推断随机变量的分布特征已无效，且会造成“伪回归”等严重后果。

数理统计方法是目前检验时间序列是否存在非平稳特性的常用方法。检验序列平稳性主要有两类方法，一类是利用时序图和相关性曲线图显示的特征进行判断的图检验方法；另一类是构造检验统计量进行假设检验的方法。图检验方法的优点是操作简单，但在判别时需要有足够经验。最常用的是通过自相关分析图和偏相关分析图判断序列的平稳性。如果一个序列在零均值化后自相关函数出现缓慢衰减或周期性衰减的情况，说明该序列可能存在某种趋势或周期波动特性(季节性变化和多年波动变化)，即表现出非平稳性。此外，也可以利用Barlett定理判断水文过程是否具有平稳性。根据Barlett定理，通过构造检验统计量并使用卡方检验识别序列的纯随机性。由于平稳序列具有短期相关性，因此具体检验过程中在短时间延迟范围内进行分析即可。单位根检验方法(Unit Root Test)是目前最常用的判断时间序列是否具有平稳性的假设检验方法，最常用的方法有Dickey-FullerTest(DF)、AugmentedDickey-Fuller Test(ADF)、Phillips-Perron test(PP)、KPSS test等检验方法，具体分析时需要查用统计特征值临界值表。

虽然平稳性和非平稳性在数学上有严格定义和界定，在其他领域也有较多的研究应用成果，但在水文学领域关于水文非平稳性问题缺乏系统性的研究，对若干基础和关键问题缺乏统一认识，且缺乏有效的水文时间序列非平稳性判别方法，导致大量的水文统计分析结果和工程水文计算结果出现很大偏差。目前关于水文非平稳性研究在认识上的最大误区是将水文过程平稳性与水文过程趋势等特性联系起来，即认为环境变化影响下某水文变量无趋势变化时属于平稳过程，当发生趋势变化时认为是非平稳过程。然而依据平稳性的严格定义可知，即使在无趋势情况下水文过程也可能会存在周期变化等确定性变化规律，或存在随机游走过程(即单位根)等假的趋势现象，并不完全属于平稳过程。因此，直接利用各种单位根检验方法判别水文时间序列的非平稳性时，由于水文序列含有不同时间尺度上的周期项，会直接影响单位根检验结果的准确性。尽管目前也有一些消除序列周期项影响的方法，但由于这些方法十分简单，无法准确刻画和描述实际水文时间序列中周期大小和周期波动幅度均复杂多变的周期项。

发明内容

针对于上述问题，本发明的目的在于提供一种判别水文时间序列非平稳性的方法，以解决现有技术中消除序列周期项影响的方法无法准确刻画和描述实际水文时间序列中周期大小和周期波动幅度均复杂多变的周期项。

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

1)选择合理的序列拟合方法、边界点处理方法和终止条件，确定具体的集合经验模态分解方法；

2)应用所确定的集合经验模态分解方法对序列x(t)进行分解，结果记为：

x(t)＝N(t)+P(t)+T(t)

其中，N(t)表示序列中含有的噪声成分，P(t)表示序列中含有的周期成分，T(t)表示序列中含有的趋势成分；

3)利用白噪声的能量扩散函数，识别出原序列中所有的周期成分P(t)；

4)去除原序列x(t)中的周期成分，将剩余成分作为新序列；

5)根据上述的新序列曲线，经判别后选择合适的模型方程进行拟合；

6)利用选择的模型方程对所述的新序列进行单位根检验，依据特征统计量判别其是否显示出非平稳特性；

7)利用自相关系数图和偏相关系数图，进一步分析所述的新序列的统计特性；

8)比较自相关系数和偏相关系数结果与单位根检验结果的一致性；若所述的新序列表现出良好的相关性，且单位根检验特征统计量接受原假设，则表明序列具有非平稳性；若二者的结果存在不一致，则认为序列不具有非平稳性。

进一步地，所述的步骤2)中集合经验模态分解方法对序列进行分解具体包括：

21)初始化：i＝1，并定义r₀＝x(t)；

22)对于r₀，识别所有的局部极值点(极大值和极小值)，然后利用三次样条曲线拟合方法分别拟合局部极大值点和极小值点，作为上包罗线和下包罗线；

23)对比并求解上下包罗线的均值曲线m_(j＝1)；

24)通过求解序列r₀与m_(j＝1)的差值，得到第一个序列h_(j＝1)；

25)将h_(j＝1)作为r₀然后重复步骤22)-24)，即不断j＝j+1直至上下包罗线关于横坐标对称，最后的h_(j)结果记为C_i；

26)重新定义r₀＝x(t)-C_i并i＝i+1，重复步骤21)-25)，当i＝N且剩余成分R_N变成一个单调函数，只包含内部极值点不能再识别任何IMF时，序列分解过程结束。

进一步地，所述的步骤5)中选择合适的模型方程进行拟合具体包括：

51)分析所述的新序列曲线随时间的变化特点，判断其是否存在明显异于零的截距和线性趋势项；

52)若所述的新序列不存在明显异于零的截距和线性趋势项，选择的模型方程为：

其中，为一系数，w(t)是一个白噪声序列；

53)若所述的新序列存在明显异于零的截距，选择的模型方程为：

其中，a表示截距项；

54)若所述的新序列存在明显异于零的截距且存在明显的线性趋势项，选择的模型方程为：

其中，b表示线性趋势项的斜率。

本发明的有益效果：

本发明的方法由于是首先分离出序列中的周期项，可以克服周期项对非平稳性判别的影响，因此较常规方法有更大的优势，可为准确进行水文统计分析和工程水文计算奠定良好基础。

附图说明

图1绘示本发明判别水文时间序列非平稳性的方法的流程图；

图2绘示算例分析中使用的具有不同周期幅值的第一类人工生成序列；

图3绘示算例分析中使用的具有不同周期大小的第二类人工生成序列。

具体实施方式

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

参照图1所示，本发明的判别水文时间序列非平稳性的方法，具体步骤包括：

1)选择合理的序列拟合方法、边界点处理方法和终止条件，确定具体的集合经验模态分解方法。

2)应用所确定的集合经验模态分解方法对序列x(t)进行分解，结果记为：

x(t)＝N(t)+P(t)+T(t)

其中，N(t)表示序列中含有的噪声成分，P(t)表示序列中含有的周期成分，T(t)表示序列中含有的趋势成分；

具体表现为：使用经验模态分解(Empirical Mode Decomposition，EMD)方法将水文序列的不同成分进行分离后，结果称为固有模态函数(IMF)。每个固有模态函数必须满足两个条件：(1)极值数目和过零点数目必须相等或至多差一个；(2)任何时间位置点，由局部极大值点决定的上包罗线和由局部极小值点决定的下包罗线必须关于横轴对称。根据IMF定义可知，每个IMF代表了序列包含的一个波动模式，而且每个IMF都是一阶平稳序列。通过将序列分解成一系列的IMFs，可揭示序列在不同时间尺度下的复杂变化特性。

假设水文时间序列由IMFs组成，可以利用EMD方法对序列进行分解。EMD方法是利用序列极值点进行分解的过程。对于水文时间序列x(t)，分解过程描述如下：

21)初始化：i＝1，并定义r₀＝x(t)；

22)对于r₀，识别所有的局部极值点(包括极大值和极小值)，然后利用适当的曲线拟合方法分别拟合局部极大值点和极小值点，并作为上包罗线和下包罗线；具体选用三次样条曲线拟合方法；

23)对比并求解上下包罗线的均值曲线m_(j＝₁₎；

24)通过求解序列r₀与m_(j＝₁₎的差值，得到第一个序列h_(j＝₁₎；

25)将h_(j＝₁₎作为r₀然后重复步骤22)-24)，即不断j＝j+1直至上下包罗线关于横坐标对称，最后的h_(j)结果记为C_i；

26)重新定义r₀＝x(t)-C_i并i＝i+1，重复步骤21)-25)。当i＝N且剩余成分R_N变成一个单调函数，只包含内部极值点不能再识别任何IMF时，序列分解过程结束。

由于EMD方法仅依靠序列局部极值点进行分解，因此具有很好的适用性，可用于分析任何具有复杂非线性和非平稳特性的时间序列。EMD方法序列分解结果中，最先分解出的IMF具有最高频率，对应着最小时间尺度。随着序列分解过程进行，分解出的IMF对应的频率逐渐降低。实际序列EMD分析过程中，会受到许多关键和难点因素的影响，例如选择合适的曲线拟合方法、边界效应、终止条件等。为提高EMD方法分析结果的可靠性，本发明中使用EMD的改进方法：集合经验模分解(Ensemble Empirical Mode Decomposition，EEMD)进行序列分解，其基本思路是不断分析叠加了噪声之后的序列，然后取平均值作为最后结果。由于EEMD方法的序列分解结果常包含一些缺乏实际物理意义的假IMFs。因此需要首先识别出真实的IMFs。此处，使用经验模态分解(EMD)方法来研究白噪声的特性，该方法依据白噪声在EMD分解之后表现出的特性，通过对比白噪声的能量扩散函数与时间序列的EMD结果，可以准确地识别出确定性成分，也可以对不确定性进行定量估计。

3)利用白噪声的能量扩散函数，识别出原序列中所有的周期成分P(t)。

4)去除原序列x(t)中的周期成分，将剩余成分作为新序列。

5)根据上述的新序列曲线，经判别后选择合适的模型方程进行拟合；

具体表现为：

由于水文时间序列常常表现出不同的复杂变化特性，因此需要根据待分析序列自身的特性选择相应合理的模型方程进行单位根检验，具体选择思路如下：

51)分析所述的新序列曲线随时间的变化特点，判断其是否存在明显异于零的截距和线性趋势项；

52)若所述的新序列不存在明显异于零的截距和线性趋势项，选择的模型方程为：

其中，为一系数，w(t)是一个白噪声序列；

53)若所述的新序列存在明显异于零的截距，选择的模型方程为：

其中，a表示截距项；

54)若所述的新序列存在明显异于零的截距且存在明显的线性趋势项，选择的模型方程为：

其中，b表示线性趋势项的斜率。

6)利用选择的模型方程对所述的新序列进行单位根检验，依据特征统计量判别其是否显示出非平稳特性。

7)利用自相关系数图和偏相关系数图，进一步分析所述的新序列的统计特性。

8)比较自(偏)相关系数结果与单位根检验结果的一致性；若所述的新序列表现出良好的相关性，且单位根检验特征统计量接受原假设，则表明序列具有非平稳性；若二者的结果存在不一致，则认为序列不具有非平稳性。

算例分析：

利用人工生成序列检验本发明方法的有效性，由于人工生成序列的周期等确定成分和是否存在单位根等情况已知，因此有利于不同方法之间进行对比分析，而实测水文时间序列是否含有单位根往往未知，无法准确判断各方法分析结果的准确性。为分析周期项对序列非平稳性检验结果的影响，在设计方案时生成两类序列。第一类序列具有相同的序列长度，且含有相同的周期值，但周期具有不同的幅值，分别记为S11、S12、S13和S14(图2)。第二类序列也具有相同的序列长度，但含有不同的周期，分别记为S21、S22、S23和S24(图3)。此外，所有序列均含有单位根过程(即随机游走过程)。选用PP test和KPSS test两种方法对各序列进行分析检验，结果见表1，其中符号“+”表示该序列中含有单位根，即表现出非平稳性；符号“-”表示该序列不含有单位根，是平稳时间序列。

表1不同方法对各人工生成序列的非平稳性检验结果

序列分析结果显示：周期值大小和周期幅值两个因素对序列非平稳性检验结果均具有很大影响。当序列中含有的周期项幅值较小时(例如S11)，常规的PP test和KPSS test方法均可以准确地识别出其含有的单位根；但是随着序列中周期项的幅值越来越大(例如S13和S14)，PP test和KPSS test方法均无法有效地识别出其含有的单位根，且PP test也无法有效地识别出S12序列中含有的单位根，因此会做出错误的判别结果。第二类序列中含有的周期项与S12序列周期项具有相同的幅值。当序列周期值较大时(例如S21和S22)，PPtest和KPSS test方法均可以有效地识别出其含有的单位根；随着序列周期值逐渐减小(例如S23和S24)，PP test方法无法识别出其含有的单位根；而KPSS test方法可以对第二类所有序列的单位根进行准确地识别。相比常规方法，本发明的方法由于首先利用EEMD方法准确识别并分离出序列中含有的周期项，因为可以克服序列周期项对单位根检验结果的影响，准确识别出序列中含有的随机游走过程，得到准确的序列非平稳性检验结果。

综合上述分析结果，可以得到关于时间序列非平稳性判别的以下几点：(1)序列中含有的周期项对序列非平稳性判别结果有很大影响，特别是周期项幅值较大或周期值较小时，往往会造成序列非平稳性判别的错误结果；(2)相比之下，KPSS test方法较PP test方法更不易受周期项的影响，因此前者效果更优，但实际中仍需要使用多种方法进行综合分析判断，以得到更准确的结果；(3)本发明方法由于是首先分离出序列中的周期项，可以克服周期项对非平稳性判别的影响，因此较常规方法有更大的优势，可为准确进行水文统计分析和工程水文计算奠定良好基础。

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

Claims (2)

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

1)选择合理的序列拟合方法、边界点处理方法和终止条件，确定具体的集合经验模态分解方法；

2)应用所确定的集合经验模态分解方法对序列x(t)进行分解，结果记为：

x(t)＝N(t)+P(t)+T(t)

其中，N(t)表示序列中含有的噪声成分，P(t)表示序列中含有的周期成分，T(t)表示序列中含有的趋势成分；

3)利用白噪声的能量扩散函数，识别出原序列中所有的周期成分P(t)；

4)去除原序列x(t)中的周期成分，将剩余成分作为新序列；

5)根据上述的新序列曲线，经判别后选择合适的模型方程进行拟合；

6)利用选择的模型方程对所述的新序列进行单位根检验，依据特征统计量判别其是否显示出非平稳特性；

7)利用自相关系数图和偏相关系数图，进一步分析所述的新序列的统计特性；

8)比较自相关系数和偏相关系数结果与单位根检验结果的一致性；若所述的新序列表现出良好的相关性，且单位根检验特征统计量接受原假设，则表明序列具有非平稳性；若二者的结果存在不一致，则认为序列不具有非平稳性；

其中，所述的步骤5)中选择合适的模型方程进行拟合具体包括：

51)分析所述的新序列曲线随时间的变化特点，判断其是否存在明显异于零的截距和线性趋势项；

52)若所述的新序列不存在明显异于零的截距和线性趋势项，选择的模型方程为：

其中，为一系数，w(t)是一个白噪声序列；

53)若所述的新序列存在明显异于零的截距，选择的模型方程为：

其中，a表示截距项；

54)若所述的新序列存在明显异于零的截距且存在明显的线性趋势项，选择的模型方程为：

其中，b表示线性趋势项的斜率。

2.根据权利要求1所述的判别水文时间序列非平稳性的方法，其特征在于，所述的步骤2)中集合经验模态分解方法对序列进行分解具体包括：

21)初始化：i＝1，并定义r₀＝x(t)；

22)对于r₀，识别所有的局部极大值点和极小值点，然后利用三次样条曲线拟合方法分别拟合局部极大值点和极小值点，作为上包罗线和下包罗线；

23)对比并求解上下包罗线的均值曲线m_(j＝1)；

24)通过求解序列r₀与m_(j＝1)的差值，得到第一个序列h_(j＝1)；

25)将h_(j＝1)作为r₀然后重复步骤22)-24)，即不断j＝j+1直至上下包罗线关于横坐标对称，最后的h_(j)结果记为C_i；

26)重新定义r₀＝x(t)-C_i并i＝i+1，重复步骤21)-25)，当i＝N且剩余成分R_N变成一个单调函数，只包含内部极值点不能再识别任何IMF时，序列分解过程结束。