华南沿海湿润区洪水频率计算：非平稳性序列研究

早在1914年，就有学者通过一次或多次分离的处理来消除趋势项进行非平稳性洪水处理[70]，再用分布对消除趋势后的数据进行拟合。但该方法的缺点是在没有研究原始变量的情况下就对数据进行转换，转换后的数据可能不再是人们感兴趣的数据，因为这种转化可能已经改变了原始数据的结构。

国内很多学者对水文序列一致性处理做了大量有意义的工作[71-72]。目前较为普遍的做法是通过还原、还现的途径来处理非平稳性水文序列。还原是把非平稳后的序列修正到变化环境之前的状态；还现是指把非平稳后发生前的天然状态修正到变化环境后。还原、还现的方法一般包括变异点前后序列与某一因子相关法、水量平衡逐项还原法、时间序列的分解与合成法以及水文模型法[73]。水量平衡逐项还原法是根据水量平衡原理，把调查得到的农业、工业、生活、蒸发、蓄水、引水等方面的水量进行逐项还原，得到断面天然径流量[74]。谢平等[71]提出了非平稳性年径流序列频率计算的随机方法，并将此方法应用在潮白河流域，通过对确定性成分的预测值和随机性成分的设计值进行合成，得到不同时期合成的频率分布。由于水文极值的变化可能是若干因子的极端事件的复合作用结果，因此利用多元极值理论[75-76]，采用引入导致水文过程发生非平稳性变化的协变量的多元极值理论也是解决非平稳性条件下的水文频率分析的可行办法。成静清和宋松柏[77]对基于条件概率和混合分布的非平稳性洪水频率计算方法[78]进行了转换和改进，得到陕北及关中地区基于条件概率的非平稳性年径流频率计算方法；得到非平稳性径流序列变化前后的频率差异，提出了洪水对水文变异前后的敏感性分析结论。基于条件概率和混合分布的非平稳性年径流频率分析计算过程，建立在水文序列变异分析基础之上，包括成因和统计分析，使其比传统的频率分析方法更符合实际。

虽然国内对水文序列一致性处理做了大量有意义的工作，但非平稳性径流序列频率计算至今尚未形成被广泛接受的统一方法[72]。以上方法多少存在一些不足，如广泛采用的还原调查法无法还原由下垫面变化后蒸发量增大、产流量减少等问题所引起的变化水量。谢平等[15]的序列分解合成法的确定性成分预测一定程度上存在外延风险。通过流域水文模型考虑物理成因分离水文序列的确定性成分，方法灵活，但模型参数的率定均只反映历史某时期的流域物理情况，对变化环境后的物理条件考虑尚不周全。

国外不少学者则通常直接基于非平稳性极值序列本身进行水文频率分析，而不需要把观测序列转换成满足一致性条件的序列后再进行水文频率分析。针对变化环境下的水文统计特征的非平稳性，目前提出建立基于时变统计参数非平稳序列洪水频率分析[14]；也有根据洪水产生机制不同，直接采用混合分布模型和基于条件概率的频率分析方法[44，79]。另外，有些研究是建立分布参数与环境因子（如大气环流）之间的对应关系，分析低频气候现象与水文极值之间的关系，将环境因子与水文极值概率分布参数建立联系，由这些因子的变化规律推求分布参数的时变规律[80]。总结关于非平稳性洪水频率分析方法，主要包括时变矩（TVM）、区域非平稳频率分析、超定量（POT）、局部似然、分位数回归和基于条件概率的混合分布等方法[16]。更多非平稳性洪水频率分析方法详见文献[16，78]。

1.2.3.1　基于时变统计参数的非平稳性洪水分析

在变化环境情况下，洪水重现期往往不是描述一场洪水的一个固定不变的属性。这就意味着在描述重现期这个概念的时候有必要引入时间基准点，这有点类似于物理上描述物体运动需要参照物一样。Strupczewski W G等[14]基于时变矩的概念提出了一种针对单站的非平稳性洪水频率计算方法，这种方法考虑了洪水时间序列第一、第二阶矩（均值m和标准差σ）的趋势成分。假设某分布的概率密度函数表达形式为f=f（x；p），其中p为分布参数向量，与前两阶矩的关系为p=p（m，σ），则其密度函数可用前两阶矩表示为f=f（x；m，σ）。其中，m、σ考虑趋势成分，与时间有关，有m=m[t；θ（m）]和σ=σ[t；θ（σ）]，θ（m）和θ（σ）分别为m和σ的参数向量，则参数p随时间变化，即p=p（t；θ）。通过TVM法把模型需要估计的参数由原来的分布参数向量p转换为m和σ的参数向量θ（m）和θ（σ），即f=f（x，t；θ），θ为由θ（m）和θ（σ）组成的参数矩阵。模型的分布函数通过第一、第二阶矩进行描述，从而可得到设计值随时间的变化关系。同一设计标准P下，设计值x P是随时间t变化的。此后，TVM法被学者成功应用于其他地区，也得到一定标准P下，洪水设计值x P随时间存在显著的变化关系。

1.2.3.2　基于成因机制的非平稳性AM洪水序列频率计算

传统洪水频率分析要求应有洪水分布一致性的假定。然而洪水序列通常是由两种或者更多的总体形成的混合分布。各个总体之间有区别是由影响因素差异所致，包括流域下垫面条件多样化、季节性洪水形成机制的多样性、区域气候的改变等。传统洪水频率分析并没有考虑到这一物理改变过程，所以导致频率分布曲线不能很好地和实测点据拟合[72]。

Waylen和Woo[81]认为年最大洪水序列分布为FT：FT（x）=F1（x）F2（x），其中两个独立的过程分别服从F1和F2两个Gumbel分布。FT的获得应包括证明洪水序列的不同产生机制、分别对不同分布估参、分别求洪水的发生概率。美国水资源委员会（USWRC）和苏联国家水文研究所（USSRNHI）强调要基于水文学和统计学原理将年最大洪水按不同成因分类，提出将季节性洪水频率分布p1（x）和p2（x）（假设只分成两类）合成为年最大洪水频率分布p（x），其表达式为p（x）=p1（x）+p2（x）-p1（x）p2（x）[44，79]。

Singh等[78]提出了一种估计非平稳性洪水序列中特大洪水的超过值概率的方法，步骤为：①分析洪水产生机制，在确保每个季节洪水服从同一分布的前提下，对非平稳性洪水按照S个不同的洪水产生机制分组；②从已划分好的数据中选取最大值组成年最大洪水序列；③分别计算频率；④用条件概率推求年最大洪水频率分布并计算参数。Alila和Mtiraoui[82]采用美国亚利桑那州中部和东南部希拉河流域的中长期水文气象资料，对非平稳性洪水频率分析计算的范围和意义进行了探究，讨论了形成非平稳性年最大洪水频率分布的可能性原因，阐述了如何选择在水文上和统计上都适用的计算非平稳性洪水序列的频率模型。

1.2.3.3　洪峰-洪水总量-历时组合特征量联合分析

水文事件由多个方面的特征属性构建而成，必须通过多个方面的特征属性进行定义和描述。实际上，洪灾严重程度很大程度上取决于洪水量级和洪水历时。因此，洪水量级和洪水历时的特征量变化过程需同时顾及。但是，在大部分的洪水频率分析中，通常只选其中的特征变量洪峰或者洪量进行单变量频率分析，显然，采用此类单变量频率分析没有考虑到洪量、洪峰、洪水历时及洪峰出现时间等特征属性之间的相关关系。

有关洪水量级和历时的研究，有两种不同的方法可以应用：洪峰-洪量分析法和流量-历时-频率（QdF）分析法[83]。值得一提的是洪峰-洪量分析法在确定洪水时间的起始和结束方面存在主观因素[83]。QdF模型是标准洪水频率模型的扩展，该模型利用一个复合方程总结了不同洪水事件在发生概率和洪水历时的洪水频率。关于QdF方法过去已经有很多研究，本书仅列举部分进行讨论。QdF方法大概来源于40年以前[45]。在20世纪90年代，Sherwood[84]、Balocki和Burges[85]等奠定了当今QdF模型的基础。Javelle等[83]基于不同洪水分布在小重现期收敛的假设，在加拿大成功将复合的方法应用到QdF模型上。随后QdF模型被成功应用于世界上其他很多地区，如马提尼克岛、法国、布基纳法索和罗马尼亚等[16，86-88]。(www.xing528.com)

总结以上研究，对具有多种特征属性的水文事件而言，单变量频率分析法不能全面准确地描述以上特征变量的统计规律以及相互之间的关系，而且，变化环境下洪水历时也在发生改变。因此，需要进行考虑洪水历时变化的流量-历时-　频率（QdF）的分析计算研究。

1.2.3.4　超定量（POT）模型

变化环境下，洪水量级已经发生明显改变。同一量级洪水下，变化环境前定义为洪水，变化环境后可能不属于洪水范畴，即存在一个洪水门限值确定标准。类似于需要定义黄海基准面和南海基准面等不同基准面的情况。变化环境后超过特定门限值流量的洪水年均发生次数明显增多，洪涝极端事件主要反映在发生的量级和频率上[3]。采用传统年最大值（AMS）取样，会忽略变化环境后的年内第二、第三大洪峰，尽管它们普遍比变化环境前洪水量级要大得多。以POT系列为样本的洪水频率分析方法是选取大于指定门限值的洪水观测值为样本进行分析。POT模型同时考虑了超定量年发生次数分布模型和超定量分布模型，比AMS法能更完整、更灵活地描述洪水及其产生过程，所以具有更多的物理相关性[89]。其思想是先假设超定量年发生次数和洪水超定量序列的分布，再从中得出年最大超定量洪水分布[90]。通常，如果门槛值足够大，超定量序列分析能比年最大值频率分析有更精确的极值分位数估计值[45，91-92]。Wang[93]采用统计试验方法，假设POT/GP模型与AMS/GEV模型选择的样本数相同，比较了两种模型在概率权重矩法（PWM）估计参数时对指定频率设计值的估计，表明在阈值较高时，POT/GP模型与AMS/GEV模型有相当的精度。Khaliq等[16]借鉴IPCC研究成果，指出在变化环境下，传统频率分析的结果将受到质疑，并将POT模型作为一种适合变化环境的频率分析方法。以往POT洪水频率研究主要集中于洪峰样本独立性判别和POT系列门限值的确定方面[16，89]，门限值选取主要有年均超定量发生次数μ法、超定量样本均值法和分散指数法等[89]。Lang[89]等提出应通过超定量样本均值法和分散指数法确定门限值区间，并选择满足μ＞2或μ＞3的较大门限值。关于年均超定量发生次数也有不同看法，Cunnane[31]认为用指数分布时μ应大于1.65，Begueria指出μ常取1.2～5，中国学者多通过试算将μ控制在2～3。戴昌军[94]等总结POT系列频率分析的关键是门限值的选择和样本的筛选，对设计洪水计算精度影响很大。

1.2.3.5　局部似然法

为了在非平稳性背景下设计合适的趋势变化模式，经常采用趋势分析来诊断水文气象数据的变化情况。因此，如果有关于数据变化的先验信息则是非常有利的，而不应仅是假设基于预先设定的几种参数模式。这方面，Davison和Ramesh[95-96]关于探索性分析极值流量的半参数方法就非常有用。在半参数方法中，对年最大日流量和POT序列，在使用标准的趋势模型拟合时应采用局部的点据拟合，通过对数据进行合理的加权，使模型的参数在各时刻可以分别进行估计。这就产生了取决于时间和极值分位数局部估计的参数估计。这一方法的目的是提供一种评估参数变化/趋势的探索性工具，还可以用来优选合适的参数模型并检验拟合度。更多关于应用局部似然法来估计极值趋势的极值和点过程模型方面，可以具体参考文献[96]。

1.2.3.6　分位数回归法

分位数回归法是一种半参数方法，由Koenkar和Basset[97]发展而来。该方法是利用不对称加权绝对偏差总和最小来估计条件分位数，不对称加权绝对偏差总和最小的确定则可应用优化技术，给予正、负残差不同的权重来获得。目前有一个重要的综合分位数回归模型。例如，截尾分位数回归（Censored Quantile Regression）模型被Powell[98-99]在经济学领域所推荐。截尾分位数回归模型当因变量实测数据被审查时可以持续地估计条件分位数。这个模型拥有重要的发展潜力，特别在有限的水文气象实测数据情况下，可以利用历史信息来提高风险估计水平，如Benito等[100]对洪水风险的估计。更多有关于截尾分位数回归模型的介绍可以参考文献[98 -99，101]，现在这种方法仍然广泛地应用于水文气象数据分析领域。

1.2.3.7　非平稳性区域洪水频率分析方法

目前很多学者对非平稳性洪水频率分析和区域洪水频率分析都有不少研究，然而把这两个方面内容联系在一起考虑的研究还相对较少。从少数几个站点的有限长度水文序列中提取的信息并不能代表区域的整体情况，且还比较容易引起设计误差[5，102]。因此，必须足够重视在区域背景下建立非平稳性，使用合适的区域方法[21]。

关于区域洪水频率分析方面，Cunderlik和Burn[13]提出了二阶非平稳的方法来进行区域洪水频率分析，该方法假设单站洪水时间序列第一、第二阶矩存在非平稳性（如时间序列均值m和标准差σ）。单站水文序列非平稳性分位数方程包括一个单站的包含位置和尺度参数的时间依赖成分，还有一个区域成分，该成分被认为在二阶非平稳性假设情况下是独立的。因为存在单站的二阶非平稳性，单站的水文序列的信息被汇集起来发展成一个稳定的生长曲线，这在组内不同的水文站点中都适用。因此，该方法被命名为“非平稳性汇集洪水频率分析”。该方法是在TVM法的基础上进行深化而得。Markiewicz等[103]对Cunderlik和Burn的区域非一致性方法提出不同看法，存在分位数表达式中的量纲不和谐现象。Leclerc等[104]基于单站非平稳性频率分析和区域分析建立多元回归模型，与平稳性分析结果对比可知忽略水文序列变化趋势会导致设计值出现偏差。

总体而言，关于非平稳性水文序列频率分析方面尚处于起步阶段，急需开展变化环境下非平稳性水文序列计算的理论与方法研究。