1. 引言
在变化环境下，由于气候变化和人类活动干扰，许多流域的水文系统已经发生了明显变化[1] ，导致水文时间序列发生了趋势或者跳跃的变化，即所谓的非一致性[2] [3] 。传统水文频率分析方法的一个基本前提假设是水文时间序列必须满足一致性的要求，即序列的分布类型和统计参数不随时间变化[4] 。由于非一致性的存在，传统基于一致性假设的水文频率分析方法的适用性受到挑战。近些年来关于非一致性水文频率分析的研究已经成为水文学科的重点研究领域。根据水文序列的维数，非一致性水文频率分析可分为单变量非一致性水文频率分析和多变量非一致性水文频率分析，目前研究的重点主要集中于单变量水文序列。本文综合国内外最新研究成果，对该领域几个重点、难点和热点方面的研究进行归纳梳理，并对其中的问题进行探讨、评述。
2. 单变量水文序列非一致性诊断
Salas将单变量一致性水文序列定义为不存在趋势、突变和周期性的序列[5] 。非一致性是相对于一致性的概念，因此非一致的水文序列可以理解为存在趋势、突变或周期性的序列。在实际中，水文序列的非一致性一般被描述为序列统计参数随时间的变化[6] 。目前，围绕水文序列中的趋势、突变成分，国内外水文学者提出了许多检验方法，其中应用最广泛的是非参数方法，包括Mann-Kendall、Spearman、贝叶斯方法等[7] -[12] 。各种检验方法的结果可能存在差异，为此，谢平等[12] 提出了水文变异诊断系统，对各种方法的结果进行综合，以期得到结果更加可靠。
然而，随着对水文序列非一致性问题研究的深入，仅从一些统计检验结果判断序列是否存在非一致性的做法开始受到质疑。非一致性检验结果很容易受到水文序列时间尺度的影响[13] ，实际中实测水文序列的长度一般从几十年到一两百年不等，远远短于流域系统历史时期完整的水文过程，因而水文序列的代表性可能并不充分。如图1中所示，随机生成长度为500且满足一致性的径流序列，截取一段长度为50序列，发现其存在显著下降的趋势，由此可见通过整个序列中较短的片段来判断水文序列是否存在非一致性可能存在偏差。正如Montanari
Figure 1. Impact of time scale on the results of nonstationarity detection. Symbol “*” indicates that the examined series has a significant linear trend at the significance level of 0.05
图1. 时间尺度对水文序列非一致性检验结果的影响，图中“*”表示序列的线性趋势能够通过显著性水平为0.05的检验
和Koutsoyiannis [14] 所指出的，变化并不意味着非一致性，而一致性序列也不是一个一成不变的过程。因此，水文序列的非一致性并不能简单地根据统计检验结果得出，还需要一个明确的水文过程变化来进行验证。
3. 单变量水文序列非一致性数学描述与归因分析
3.1. 基于时间变量的非一致性数学描述
对单变量非一致性水文序列进行频率分析最常用的两种方法为还原/还现方法[7] 和时变矩法[15] -[17] 。然而，无论哪种方法其前提都是对水文序列的非一致性进行数学描述，这直接关系到频率分析的结果。谢平等[7] 对还原/还现方法做了较为详细的理论探讨及应用研究。在还原/还现方法中，认为非一致性水文序列由确定性成分和随机性成分构成，确定性成分通常被定义为序列的非一致性成分，而随机性成分为序列的一致性成分。在对非一致性水文序列进行频率分析之前要先确定序列中的确定性成分，然后通过还原或还现的方式剔除确定性成分，进而将原序列转换为满足一致性的序列，最后应用基于一致性假设的方法对序列进行频率分析。Strupczewski等[15] -[17] 最先提出时变矩法，并将线性趋势和二次三项式趋势直接嵌入到极值洪水序列分布的一阶和二阶矩中来描述水文序列的非一致性，从而能够直接对原序列进行频率分析。在Strupczewski等[15] -[17] 工作基础上，Katz等[18] 进一步提出直接将趋势性嵌入到极值洪水序列分布的统计参数当中，原理更加清晰直观。
在描述水文序列趋势性变化中，最常用的方法是将水文序列的分布参数
写成时间
的某种函数，如线性、指数、多项式、对数、样条插值函数等，表达式如下：
(1)
使用时间作为水文序列分布参数的解释变量可以比较直观地描述实测期内水文序列的趋势变化，在实际中已经被广泛应用。Kharin和Zwiers [19] 选取时间作为解释变量，在考虑GEV分布位置和尺度参数随时间变化的情况下，研究了非一致性条件下全球气温和降水极值在未来时期的变化情况。Cunderlik和Ouarda [20] 将时变统计参数的思路引入到区域洪水-历时-频率模型中，结果表明，相比于不考虑非一致性情况，该方法较好地克服了某一重现期下设计洪水存在的偏差问题。杜涛等[21] 选取时间为协变量，研究了渭河流域暴雨时间序列统计分布随时间的变化情况，最终发现，未来时期的设计暴雨量级有显著增大的趋势。
2005年，Rigby和Stasinopoulos [22] 提出了位置、尺度和形状的广义加法模型(GAMLSS，Generalized Additive Models for Location Scale and Shape)，该模型是一种(半)参数回归模型可以描述随机变量序列任何统计参数与解释变量之间的线性或非线性关系。同时，GAMLSS模型可以描述的随机变量分布类型范围比较广泛，包括一系列高偏度和高峰度的离散和连续分布，尤其适合拟合具有超峰度和平顶峰度、高度正偏/负偏而不服从传统指数分布的随机变量序列。近些年，以R软件为平台的GAMLSS工具包[23] 的出现，为时变矩法在非一致性水文频率分析中的应用提供了强有力的工具。Villarini等[24] [25] 应用基于耿贝尔、伽马、对数正态等两参数分布的GAMLSS模型分别对美国20世纪的年最大洪水序列以及罗马地区长期的降水和气温序列进行了非一致性研究。江聪和熊立华[26] 应用GAMLSS模型对我国长江宜昌水文站1882~2009间的年平均流量序列和年最小月流量序列的一阶、二阶和三阶矩分别进行趋势分析。
单纯地建立水文序列统计参数与时间变量的函数关系只能描述序列的趋势性变化，但无法描述序列的跳跃式变化(突变)。另一方面，由于时间变量与水文序列之间并不存在物理成因上的相关关系，因此无法解释导致水文序列非一致性的原因。此外，将序列的非一致性描述为时间的函数很可能带来比一致性模型更大的不确定性，尤其是用于对未来的趋势变化预测方面[27] 。
3.2. 基于物理因子的非一致性数学描述与归因分析
近些年来，国内外水文学者开始将一些与水文序列存在相关关系的物理因子作为水文序列统计参数的解释变量
(2)
式中
为m个解释变量，其中包括气候指数或人类活动指数。气候指数一般包括降水、气温以及大气环流指数等，人类活动指数主要考虑人口、土地利用变化、水库影响以及农业灌溉引水等。Villarini等[28] 以遥相关大气环流因子作为解释变量对洪水序列的非一致性进行了分析。López和Francés [29] 对比了分别以时间和物理因子(包括气候因子及水库系数)作为解释变量的模型在模拟和预测水文序列非一致性中的效果，结果表明物理因子要明显优于时间变量。Xiong等[30] 基于时变矩模型，建立了渭河年径流序列统计参数与降水、气温、灌溉面积等因子的函数关系，分析了年径流序列非一致性的成因。Liu等[31] 选取了太平洋年代际涛动指数(PDO)、海平面气温差异指数(NINO3)、太阳黑子数、面平均气温以及面总降水等作为解释变量，对长江流域宜昌站枯水流量进行非一致性频率分析。熊立华和江聪[32] 以前期影响雨量和水土保持面积作为解释变量，分析了导致渭河年最大洪峰流量序列非一致性的原因。这些研究都表明，相比于简单的时间变量，物理因子能够更好地描述以及预测水文序列的非一致性，而且还能够从物理成因的角度对序列中的非一致性进行归因分析。
在表述水文序列分布参数与解释变量关系过程中，最常用的方法是统计方法，即建立分布参数对解释变量的回归方程[28] -[32] 。这种方法固然简单，但是往往忽略了水文过程中的内在机理。近些年，一些基于水文物理机制的统计模型开始被应用于水文序列分布的推导中。Gottschalk等[33] 以及Yu等[34] 基于退水模型，建立了干旱天数和退水参数这两个随机变量的联合分布，推求得到低定量取样枯水流量样本的频率分布。Xiong等[35] 基于Budyko假设，由降水、潜在蒸散发等随机变量推导了年径流序列的频率分布，该方法证明可以用于无资料地区径流序列的频率分布推求。在对水文序列的非一致性进行数学描述及其归因分析中，如何考虑水文过程中的内在机制将是今后重要的研究内容。
4. 非一致性条件下单变量随机事件重现期定义与估计
4.1. 一致性条件下重现期定义与估计
如今，有两种较为常用的单变量水文极值事件重现期的定义。其一是Wigley [36] [37] 提出的期望等待时间概念：即从初始年起，直到下一次出现超过一定量级的极值事件的期望等待时间为
年。另一种是Parey等[38] [39] 提出的期望超过次数概念：即从初始年起，在未来
年内出现超过一定量级的极值事件的次数为1。在一致性条件下，通常认为水文极值序列的分布类型以及统计参数都不随时间变化(图2)。以量级为
的洪水事件为例，两种重现期的估计方法如下：
1) 期望等待时间概念。用
表示从初始年
起，首次出现超过洪水量级
的年份。容易得知
服从几何分布[40] ：
(3)
于是有
在未来时期的重现期为[41] [42] ：
(4)
2) 期望超过次数概念。用
表示从初始年
起，在未来
年内出现超过洪水量级
的次数。于是有
，
为指示函数。可知
服从二项分布[43] ：
(5)
于是有：
(6)
由图2可知，一定量级洪水事件
所对应的超过概率
不随时间变化，因此两种概念下相应的重现期均为
，这与传统估计方法得到的结果相同。
4.2. 非一致性条件下重现期估计方法
在非一致性条件下，由于水文序列的概率分布随时间发生变化，导致一定量级洪水事件的超过概率随时间发生变化(图3)，根据传统的重现期估计方法
，给定量级的洪水事件
在未来每一年都有一个重现期值与之相对应 [44] ，很难应用于工程实际。在此情况下，Salas和Obeysekera [40] 将期望等待时间概念扩展到非一致性情况，并分别选取具有上升、下降趋势或者跳跃突变的水文气象序列进行实例验证，结果表明考虑非一致性情况下所得重现期计算结果更为合理可信。Cooley [43] 同时将期望等待时间和期望超过次数两种概念下的重现期应用于非一致性水文极值序列，并进一步推求出非一致性条件下与给定重现期相对应的设计洪水量级。根据文献Salas和Obeysekera [40] ，两种重现期在非一致性条件下的估计方法如下：
1) 期望等待时间概念。
定义同式(3)，在非一致性条件下，由图3可明显发现，量级为
的洪水事件在未来各年内超过概率不再为常数，于是
的概率密度函数为：
(7)
Figure 2. Probability distribution of flood series under stationary condition
图2. 一致性条件下洪水序列的概率分布
Figure 3. Probability distribution of flood series under nonstationary condition
图3. 非一致性条件下洪水序列的概率分布
则第一次出现量级超过洪水
的期望等待时间为[39] ：
(8-1)
(8-2)
2) 期望超过次数概念。
定义同式(5)，在非一致性条件下，
将不再服从二项分布，此时有：
(9)
解上式即可得到非一致性条件下量级为
的洪水事件的重现期。
根据公式(8)和(9)，在非一致性条件下，任意给定的洪水事件的重现期依然是一个确切的值，更加符合工程实际需要。Du等[45] 在此基础上，将气象协变量引入到非一致性枯水流量重现期及风险分析中，结果表明考虑气象协变量的非一致性设计结果相比于一致性情况和仅以时间为协变量的非一致性情况更加合理可信。
5. 多变量非一致水文序列频率分析
水文事件(过程)一般具有多个方面的特征属性，是一个包含频域、时域和空间域的复杂过程，单一水文变量很难描述出整个水文事件的真实特征，例如在一个洪水事件通常考虑洪峰、洪水总量和洪水历时三个要素，一个干旱事件一般也要考虑干旱烈度干旱历时等要素，因此单变量频率分析也就不能完整地评估水文事件的发生概率。近几年，一些研究开始考虑多变量水文序列的非一致性问题。Chebana等[46] 应用类似于非参数的Mann-Kendall和Spearnman检验方法检验了多变量水文序列的趋势性。Ben Aissia等[47] 基于多元线性回归中多变点检验的贝叶斯方法分析了具有相关性的两变量洪水序列的变点。
多变量水文序列可以被认为是几个具有相关关系的单变量水文序列的组合，因此多变量水文序列的概率分布包括单个变量对应的边缘分布和单个变量间的相关性结构两方面内容。很显然多变量水文序列的非一致性也应包括两方面的内容，即边缘分布的非一致性和相关结构的非一致性。边缘分布的非一致性问题本质上就是单变量水文序列的非一致性问题，因此多变量水文序列的非一致性的核心是相关结构的非一致性。由于Copula函数可以描述水文变量间的相关性结构，能够构建任意边缘分布的联合分布，已经被广泛用于多变量水文频率分析计算[48] -[59] 。一些基于Copula函数的方法已经提出来应对多变量非一致性水文频率分析问题。冯平和李新[60] 在对单变量洪水序列进行变点分析的基础上用混合分布分别拟合了洪峰-洪量的边缘分布，并用Copula函数构建了峰量的联合分布。但是以上研究仅考虑了洪峰-洪量序列边缘分布的非一致性，并未考虑相关性结构的非一致性。Bender等[61] 通过50年滑动窗描述了莱茵河1820~2011间洪峰-洪量序列边缘分布参数和Copula函数参数的时变性，发现两者均存在非一致性。Jiang等[62] 分别考虑以时间和水库系数作为边缘分布参数和Copula函数参数的解释变量来描述汉江干流两个相邻水文站的年枯水序列的非一致性，结果发现水库系数较时间能够更好地描述两变量枯水序列的非一致性。根据文献[62] 的思路，两变量非一致性水文序列(
)的联合分布由如下表示
(10)
式中，
和
分别为
和
时变的边缘分布，其中
和
分别为
和
分布的时变参数向量，
是以
为时变参数的Copula函数。由公式(10)可见，多变量水文序列的非一致性可以通过边缘分布参数和Copula函数参数的时变性来进行描述。
需要指出的是，现有基于Copula函数的多变量非一致水文频率分析方法主要针对两变量水文序列，然而许多水文事件需要三个甚至三个以上的随机变量进行描述，因此三变量乃至更高维数水文序列的非一致性需要进一步研究。多变量水文序列重现期的计算相比与单变量重现期更加复杂，目前最常用的多变量重现期计算方法有“AND”、“OR”以及“KEN”等三种重现期[54] 。但是，这些重现期在非一致条件下的计算问题还尚未被研究。
6. 结论与展望
本文在总结国内外非一致性水文频率分析最新研究进展的基础上，归纳梳理了该领域几个重点、难点和热点方面，并对其中的问题进行了初步探讨，主要结论与展望如下：
1) 时间尺度是影响水文序列非一致性检验的重要因素，水文序列非一致性诊断并不是简单的统计检验问题，检验的结果必须与水文过程的具体变化相结合，从机理方面支撑检验结果。
2) 时变矩模型是描述单变量水文序列非一致性的有力数学工具，以物理因子作为水文序列统计参数的解释变量较时间变量能够更好地描述与预测水文序列中的非一致性，并且可以对非一致性进行归因分析。但现有的研究往往忽略了水文过程中的物理机制。引入基于水文物理机制的统计模型来描述及其解释水文序列的非一致性需要重点考虑。
3) 在非一致性条件下，传统计算重现期的估计方法可能会失效，如何将Salas和Obeysekera [40] 提出的非一致性条件下单变量水文事件重现期估计方法应用到工程实践中需要进一步研究探讨。
4) 目前国内外对多变量非一致性水文频率计算方法的研究还处于起步阶段。针对多变量水文序列非一致性的诊断、数学描述与归因分析及其非一致条件下多变量水文事件联合重现期的推求都需要更为广泛和深入的研究。
基金项目
国家自然科学基金重大项目(51190094)；国家自然科学基金项目(51479139)。