一种判别水文时间序列非平稳性的方法
【技术领域】
[0001] 本发明属于水文科学技术领域，尤其指代一种判别水文时间序列非平稳性的方 法。
【背景技术】
[0002] 揭示自然界水循环过程的变化特性和演变规律是解决复杂交叉水问题的重要基 础，因此一直是水文学领域的核心研究课题。实际中，欲准确认识水循环过程的复杂变化 特性是十分困难的科学问题。特别是近几十年来，受气候变化和人类活动的持续影响，许 多流域的天然水循环过程受到破坏，各主要水循环要素的周期和趋势等特性发生变异，降 雨一一径流响应关系明显改变，洪涝灾害等极端水文事件频发，且未来旱涝灾害发生的不 确定性和风险会愈来愈大。大量水文观测结果显示，许多流域实测水文资料呈现出明显的 非平稳性，给水利工程设计、水文模拟计算、水资源规划等一系列涉水工作带来很大困难。
[0003] 依据数理统计学理论，数据非平稳性主要是指大样本的统计推断基础一一数据 一致性要求遭到破坏。自2008年Milly等人在Science发表论文《StationarityIs Dead:WhitherWaterManagement?》以来，水文过程非平稳性正逐渐成为国际水文学界的 研究热点，尤其是全球环境变化研究领域重点关注的难点问题之一。判别水文过程是否具 有平稳性有十分重要的意义。水文过程平稳性可以简化水文分析计算工作的难度，提高序 列统计特征量的估计精度。从随机水文学角度看，常用的水文时间序列模型是以通过揭示 时间序列自身变化规律为主线发展起来的方法论。由于经典回归分析方法的一个重要假设 是数据的平稳性，因此识别水文过程是否具有平稳性也是开展回归分析和水文时间序列模 拟预报的前提。如果水文时间序列表现出明显的非平稳性，经典回归分析的基础和有效性 已不适用，用样本均值和方差等统计指标推断随机变量的分布特征已无效，且会造成"伪回 归"等严重后果。
[0004] 数理统计方法是目前检验时间序列是否存在非平稳特性的常用方法。检验序列 平稳性主要有两类方法，一类是利用时序图和相关性曲线图显示的特征进行判断的图检验 方法；另一类是构造检验统计量进行假设检验的方法。图检验方法的优点是操作简单，但 在判别时需要有足够经验。最常用的是通过自相关分析图和偏相关分析图判断序列的平 稳性。如果一个序列在零均值化后自相关函数出现缓慢衰减或周期性衰减的情况，说明 该序列可能存在某种趋势或周期波动特性（季节性变化和多年波动变化），即表现出非平 稳性。此外，也可以利用Barlett定理判断水文过程是否具有平稳性。根据Barlett定 理，通过构造检验统计量并使用卡方检验识别序列的纯随机性。由于平稳序列具有短期 相关性，因此具体检验过程中在短时间延迟范围内进行分析即可。单位根检验方法（Unit RootTest)是目前最常用的判断时间序列是否具有平稳性的假设检验方法，最常用的方 法有Dickey-FullerTest(DF)、AugmentedDickey-FullerTest(ADF)、Phillips-Perron test(PP)、KPSStest等检验方法，具体分析时需要查用统计特征值临界值表。
[0005] 虽然平稳性和非平稳性在数学上有严格定义和界定，在其他领域也有较多的研究 应用成果，但在水文学领域关于水文非平稳性问题缺乏系统性的研究，对若干基础和关键 问题缺乏统一认识，且缺乏有效的水文时间序列非平稳性判别方法，导致大量的水文统计 分析结果和工程水文计算结果出现很大偏差。目前关于水文非平稳性研究在认识上的最 大误区是将水文过程平稳性与水文过程趋势等特性联系起来，即认为环境变化影响下某水 文变量无趋势变化时属于平稳过程，当发生趋势变化时认为是非平稳过程。然而依据平稳 性的严格定义可知，即使在无趋势情况下水文过程也可能会存在周期变化等确定性变化规 律，或存在随机游走过程（即单位根）等假的趋势现象，并不完全属于平稳过程。因此，直 接利用各种单位根检验方法判别水文时间序列的非平稳性时，由于水文序列含有不同时间 尺度上的周期项，会直接影响单位根检验结果的准确性。尽管目前也有一些消除序列周期 项影响的方法，但由于这些方法十分简单，无法准确刻画和描述实际水文时间序列中周期 大小和周期波动幅度均复杂多变的周期项。
【发明内容】
[0006] 针对于上述问题，本发明的目的在于提供一种判别水文时间序列非平稳性的方 法，以解决现有技术中消除序列周期项影响的方法无法准确刻画和描述实际水文时间序列 中周期大小和周期波动幅度均复杂多变的周期项。
[0007] 为达到上述目的，本发明的一种判别水文时间序列非平稳性的方法，包括步骤如 下：
[0008] 1)选择合理的序列拟合方法、边界点处理方法和终止条件，确定具体的集合经验 模态分解方法；
[0009] 2)应用所确定的集合经验模态分解方法对序列x(t)进行分解，结果记为：
[0010]x(t)= N(t)+P(t)+T(t)
[0011] 其中，N⑴表不序列中含有的噪声成分，P⑴表不序列中含有的周期成分，T(t) 表示序列中含有的趋势成分；
[0012] 3)利用白噪声的能量扩散函数，识别出原序列中所有的周期成分P(t);
[0013] 4)去除原序列x(t)中的周期成分，将剩余成分作为新序列；
[0014] 5)根据上述的新序列曲线，经判别后选择合适的模型方程进行拟合；
[0015]6)利用选择的模型方程对所述的新序列进行单位根检验，依据特征统计量判别其 是否显示出非平稳特性；
[0016]7)利用自相关系数图和偏相关系数图，进一步分析所述的新序列的统计特性；
[0017] 8)比较自（偏）相关系数结果与单位根检验结果的一致性；若所述的新序列表现 出良好的相关性，且单位根检验特征统计量接受原假设，则表明序列具有非平稳性；若二者 的结果存在不一致，则认为序列不具有非平稳性。
[0018] 进一步地，所述的步骤2)中集合经验模态分解方法对序列进行分解具体包括：
[0019] 21)初始化：i=1，并定义！〇= x(t);
[0020] 22)对于r。，识别所有的局部极值点（极大值和极小值），然后利用三次样条曲线 拟合方法分别拟合局部极大值点和极小值点，作为上包罗线和下包罗线；
[0021] 23)对比并求解上下包罗线的均值曲线mu = 1);
[0022] 24)通过求解序列r。与的差值，得到第一个序列hQ ;
[0023] 25)将1^ = 1)作为。然后重复步骤22)-24)，即不断」=」+1直至上下包罗线关 于横坐标对称，最后的结果记为Ci;
[0024]26)重新定义rQ=x(t)_C;并i=i+1，重复步骤21)_25)，当i=N且剩余成分 RN变成一个单调函数，只包含内部极值点不能再识别任何MF时，序列分解过程结束。
[0025] 进一步地，所述的步骤5)中选择合适的模型方程进行拟合具体包括：
[0026] 51)分析所述的新序列曲线随时间的变化特点，判断其是否存在明显异于零的截 距和线性趋势项；
[0027] 52)若所述的新序列不存在明显异于零的截距和线性趋势项，选择的模型方程 为：
[0028]
[0029] 其中，-为一系数，w(t)是一个白噪声序列；
[0030] 53)若所述的新序列存在明显异于零的截距，选择的模型方程为：
[0031]
[0032] 其中，a表示截距项；
[0033] 54)若所述的新序列存在明显异于零的截距且存在明显的线性趋势项，选择的模 型方程为：
[0034]
[0035] 其中，b表示线性趋势项的斜率。
[0036] 本发明的有益效果：
[0037] 本发明的方法由于是首先分离出序列中的周期项，可以克服周期项对非平稳性判 别的影响，因此较常规方法有更大的优势，可为准确进行水文统计分析和工程水文计算奠 定良好基础。
【附图说明】
[0038] 图1绘示本发明判别水文时间序列非平稳性的方法的流程图；
[0039]图2绘示算例分析中使用的具有不同周期幅值的第一类人工生成序列；
[0040]图3绘示算例分析中使用的具有不同周期大小的第二类人工生成序列。
【具体实施方式】
[0041] 为了便于本领域技术人员的理解，下面结合实例与附图对本发明作进一步的说 明，实施方式提及的内容并非对本发明的限定。
[0042]参照图1所示，本发明的判别水文时间序列非平稳性的方法，具体步骤包括：[0043] 1)选择合理的序列拟合方法、边界点处理方法和终止条件，确定具体的集合经验 模态分解方法。
[0044] 2)应用所确定的集合经验模态分解方法对序列x(t)进行分解，结果记为：
[0045]x(t) =N(t)+P(t)+T(t)
[0046] 其中，N⑴表不序列中含有的噪声成分，P⑴表不序列中含有的周期成分，T(t) 表示序列中含有的趋势成分；
[0047] 具体表现为：使用经验模态分解（EmpiricalModeDecomposition，EMD)方法将 水文序列的不同成分进行分离后，结果称为固有模态函数（mf)。每个固有模态函数必须 满足两个条件：（1)极值数目和过零点数目必须相等或至多差一个；（2)任何时间位置点， 由局部极大值点决定的上包罗线和由局部极小值点决定的下包罗线必须关于横轴对称。根 据頂F定义可知，每个頂F代表了序列包含的一个波动模式，而且每个頂F都是一阶平稳序 列。通过将序列分解成一系列的IMFs，可揭示序列在不同时间尺度下的复杂变化特性。
[0048] 假设水文时间序列由頂Fs组成，可以利用EMD方法对序列进行分解。EMD方法是 利用序列极值点进行分解的过程。对于水文时间序列x(t)，分解过程描述如下：
[0049] 21)初始化：i=1，并定义1~。=x(t);
[0050] 22)对于r。，识别所有的局部极值点（包括极大值和极小值），然后利用适当的曲 线拟合方法分别拟合局部极大值点和极小值点，并作为上包罗线和下包罗线；具体选用三 次样条曲线拟合方法；
[0051] 23)对比并求解上下包罗线的均值曲线mu = 1);
[0052] 24)通过求解序列r。