以极端水灾害发生规律与机理为重大科学问题的全球变化研究已成为当今重大科学前沿^[1-3]。Milly等指出：流域径流序列形成的环境背景“一致性”已不复存在,采用现有的工程水文分析方法设计的流域开发利用工程、防洪和抗旱工程及其运行调度等,将面临由变化环境带来的设计频率失真风险^[4]。水文频率分析计算需满足独立一致同分布假设,识别、检验并分离非一致性水文序列的暂态成分是正确评估极值水灾害的基础。很多学者对非一致性径流序列频率计算做了大量有意义的工作,常用方法包括：变异点前后系列与某一参数的关系分析法、时间系列的分解与合成方法以及水文模型方法3种^[5]。但非一致性径流序列的频率计算至今尚未形成被广泛接受的统一方法^[5]。变化环境导致非一致性条件下的水文频率研究近年来在国外受到格外关注^[6,7]。针对水文序列非一致性,可建立基于时变统计参数的非一致性洪水频率方法^[6]。Waylen和Woo等提出直接基于非同分布的极值样本系列进行频率分析的混合分布模型和基于条件概率的分析方法^[8,9]。更多非一致洪水频率分析方法详见^[9,10]等文章。由于非一致性洪水序列的统计参数（如平均值和标准差）和分布线型自身时刻发生改变,相应设计流量的不确定性也将发生变化^[6]。传统极值流量分析的“极值理论”根据历史序列率定一套固定分布参数,推求得的分布曲线将局限于序列所处历史时期的流域物理条件,不能反映变化环境下洪水频率分布曲线的变化。东江流域在气候变化和人类活动尤其是水库调蓄下,洪水序列样本的产生条件和统计特性已呈现非一致性^[11]。本文根据洪水序列变化综合诊断结果,试图从分离变化环境后水文极值演变规律来探索非一致性序列的频率计算问题：揭示变化环境下非一致性洪水序列频率分布参数及线型响应规律;并探讨非一致性洪水序列洪水设计值演变特征及影响机制。研究对于科学理解在当前气候变化与人类活动双重影响下流域的防洪安全、生态环境演变具有重要科学与现实意义。

1 研究区与数据来源

东江干流长度562 km,流域总面积27 040 km²,属南亚热带季风气候区。流域水资源已经被高度开发利用,水库控制全流域面积达51%^[12]。

径流资料为流域上游龙川和中游河源2个水文站的逐日流量数据,源自广东省水文局（表1）。流域28个逐日降雨站点资料源于国家气象中心,相关主要水库信息见表2。

2 研究方法

2.1 洪水序列非一致性识别

基于秩的非参数Mann-Kendall（M-K）趋势突变检验法不要求分析数据服从某一概率分布且具有和参数检测方法相同的检测能力,受到国际水文组织强烈推荐^[13,14]。在1%显著水平下,若 Z>2.56,则序列没有趋势的假设被拒绝,序列产生非一致性。

2.2 洪水频率线型与参数估计

选用常用的皮尔逊Ⅲ分布（PT3）、广义逻辑分布（GLO）、广义极值分布（GEV）、威布尔分布（Weibull）、耿贝尔分布（Gumbel）、二参数/三参数对数正态分布（LN2/LN3）和对数皮尔逊Ⅲ分布（LP3）8种分布函数（表3）拟合东江年最大流量序列^[15,16]。

极大似然法估参（MLE）表现出来的有偏性比其他方法要弱^[17],且能提供更加一致的参数估计,受到美国联邦应急管理署（FEMA）推荐。采用MLE法估算参数。

2.3 最优模型选择

优选分布线型方法主要包括基于假设检验和一定准则的拟合优选^[18]。假设检验法存在一定主观性和显著性水平的不确定性^[18]。选用3种性能较好的信息准则方法优选分布线型（表4）,包括AIC(Akaike information criterion),BIC（Bayesian information criterion）和AICc（Corrected Akaike Information）^[19],量符号分别表示为A_AIC、B_BIC和A_AICc。

3 结果分析

3.1 年最大流量变化过程

M-K检验对年最大流量序列进行非一致性识别（图1）。龙川、河源站均通过显著性水平α=0.01检验,呈现非一致性。龙川突变点出现在1974年,1974~2004年流量有显著下降的趋势;河源站具有更显著的下降变化趋势。

显示原图| 下载原图ZIP| 生成PPT

图1   东江流域年最大流量M-K检验

Fig.1   M-K test of annual maximum daily flow series in Dongjiang River Basin

东江年最大流量变化时间与水库建设时间比较吻合（表2）。根据M-K检验并结合水利工程建设时间,确定龙川站以1974年,河源站以1964年为变化环境的时间分割点。

3.2 概率分布函数选择

分别以年最大流量全部序列龙川（1954~2009）、河源（1951~2010）和变化环境后年最大流量序列龙川（1975~2009）、河源（1964~2010）拟合配线。AIC、BIC和AICc拟合优度检验（表5）,MLE法估参见（表6）。

年最大流量全部序列和变化环境后序列拟合得到的最优分布线型无差异（表5）。龙川站LN2分布拟合最好,GEV分布拟合最差;河源站LN2分布拟合最好,PT3分布拟合最差,LN2分布较适合该地区。2个站最优分布对水文极值的描述更符合水文极值实际分布情况（图2）。本文运用LN2分布来研究年最大流量变化特征。

显示原图| 下载原图ZIP| 生成PPT

图2   东江年最大流量频率分布曲线对比

Fig.2   Comparison the flood frequency distribution curve of annual maximum flow in Dongjiang River

洪水频率曲线选择及其高水尾端特征与序列的偏态系数C_s、变差系数C_v有较大关系。通常,正偏情况,当C_s愈大,频率曲线中部愈向左偏,且高水上段越陡,下段越平缓。随C_s和C_v增大,频率曲线高水端会越来越陡^[21]。采用30 a滑动平均法,即每30 a作为一个数据样本,逐年依次向后推算,往后逐年滑动计算龙川站序列起始年为（1954~1980年）和河源站（1951~1981）年最大流量序列C_v、C_s滑动值。C_v和Cs变化过程相似,本文只列C_s变化过程。龙川和河源站洪水序列C_s值陡涨陡落（图3）,最大C_s值龙川和河源站分别达到3.07和2.35,曲线高水尾端应该不太平缓。因此,拟合结果显示薄尾PT3和Weibull分布均不适合,而混尾LN2分布更适合东江流域（表5）。

显示原图| 下载原图ZIP| 生成PPT

图3   东江年最大流量序列参数30 a滑动C_s变化过程

Fig.3   Change process of 30 a moving C_s value of annual maximum flow series in Dongjiang River

从年最大流量的30 a滑动C_s值变化过程探讨频率分布线型特征演变规律。年最大流量变化环境后序列拟合线型比全部序列拟合线型高水端存在由“陡峭”变为“平缓”的趋势（图2）。龙川站1966年（指代序列1966~1995年,下同）年最大流量30 a滑动C_s值显著下降、河源站自1964年C_s值显著下降（图3）;龙川站C_s值变化要比河源站更为剧烈,这与河源站有更大的水库调蓄库容有关。参数变化趋势表明选用不同流量时间序列估计的线型参数有较大差异。C_s减小使变化环境后序列拟合的线型比全部洪峰序列拟合得到线型尾部更低（图2）,导致设计洪水流量也不一致。

3.3 变化环境前后各重现期对应年最大流量(m³/s)的设计值

龙川和河源站变化环境（下简称“变化”）后年最大流量序列与全部序列设计洪水流量有显著差异（图4,表7）。全部序列龙川（1954~2009）、河源（1951~2010）拟合线型指定重现期的设计流量均显著大于变化后序列龙川（1975~2009）、河源（1964~2010）设计流量(表7)。

显示原图| 下载原图ZIP| 生成PPT

图4   重现期100 a的设计洪峰流量变化趋势

Fig.4   Change process of design flood peak flow (Return period=100 years)

3.4 指定标准下设计流量变化过程

推求30 a滑动序列频率曲线,对比30 a滑动序列设计洪水流量值变化过程。对同一设计标准P（百年一遇）,东江非一致性洪水序列洪水设计值随时间由大变小（见图4）。

龙川站设计流量值先是平稳,逐渐减小至趋于稳定,再略有回升。设计流量的最大值为9 189 m³/s（序列起始时间1957年）;1972年达到最小值2 305 m³/s。与龙川站类似,河源站设计流量最大值为11 125 m³/s（1953年）;1972年达到最小值4 072 m³/s。

3.5 频率分布特征与设计流量变化成因分析

变化环境后东江特大洪水比历史上任何时段都少。以龙川站全部序列拟合分析,发生6次大于十年一遇的洪水中在建库后（1975年）仅有1次。1954~1974年年最大日流量为2 374 m³/s,远大于变化后的1 211.4 m³/s。水文环境变异原因可归结为气候变化和人类活动影响。

1） 降雨变化影响。洪水频率分布参数和线型随年最大流量序列变化而改变,而年最大流量变化与区域降雨变化密切相关。将雨量站点资料用算术平均法获得面降雨资料。结合流域汇流历时,取年最大流量前7 d内的降雨总量（前7 d降雨）进行趋势分析（图5）。

显示原图| 下载原图ZIP| 生成PPT

图5   东江流域年最大流量及前7 d降雨量的变化过程

Fig.5   Change process of annual maximum flow and sum of previous seven-days rainfall value in Dongjiang River Basin

龙川和河源站年最大流量序列与前7 d降雨序列变化趋势一致（图5）。但年最大日流量变化趋势均达到0.01显著性水平,kendall趋势检验龙川和河源分别为-3.18和-4.20;而前7 d降雨下降趋势不显著,趋势检验龙川和河源仅为-1.51和-1.66。说明东江变化环境下洪水频率分布线型改变除降雨影响外,还有其他因素。

同样计算前7 d降雨自1956年以来C_s值30 a的滑动值变化,得到洪水和降雨序列的C_s滑动值序列（图6）。龙川站前7 d降雨序列与年最大流量序列的30 a滑动C_s值变化趋势基本一致,但存在不同。龙川洪水30 a滑动C_s值在1966年前（序列1966~1995年）为缓慢上升,而龙川降雨30 a滑动C_s值则在1966年前为下降,说明洪水只是一定程度上受降雨影响;河源站年最大流量序列比前7 d降雨序列影响变化趋势更为剧烈,1964年前上升明显,受降雨影响大;1964年后显著下降,随后平稳,受降雨影响程度减弱（图6）。说明东江分布线型特征改变影响因子中降雨只是部分原因。

显示原图| 下载原图ZIP| 生成PPT

图6   东江年最大流量和前7 d降雨量序列参数30 a滑动C_s变化过程

Fig.6   Change process of 30 a moving Cs value of annual maximum daily flow series and previous seven-days rainfall series in Dongjiang River

2） 土地利用变化影响。归一化植被指数（NDVI）为美国国家航天航空局的全球植被指数变化研究数据GIMMS,可反映流域植被覆盖变化情况。经过大气校正、MVC法等前期处理后,得到东江流域1982~2006年的逐年NDVI数字影像集。东江1982~2006年植被覆盖率存在不显著下降趋势,NDVI龙川和河源站Kendall趋势检验值分别为-0.841和-1.588。流域内1980年代来经济发展迅速,农田面积减少,地面不透水面积增加^[12]。森林覆盖率降低,不透水面积增加缩短了洪水汇流时间,洪水强度理应增强,土地利用变化对东江年最大流量变化趋势应该有上升影响。然而,东江年最大流量却存在显著下降趋势。因此,土地利用变化不是影响洪峰流量变化的主要原因。

3） 水利工程调蓄影响。前面分析发现年最大流量和统计参数主要集中在序列起始时间为1964年左右的流量序列中发生突变。考虑到新丰江水库于1963年建成开始下闸蓄水,可见流量序列和统计参数发生变化的时间与水库建设的时间非常吻合。据统计新丰江水库建成前后（至1979年）的6次大洪水分析,河源站削减洪峰流量1 570~6 560 m³/s,降低水位1.47~4.96 m^[12]。大中型水库对洪峰有较大的削减作用,因此河源年最大流量在1964年之后、龙川站在1974年后波动幅度减小,洪水频率最优线型由陡峭变为平缓。

4 结 论

运用8种概率分布函数对东江流域龙川和河源站极值流量进行系统分析。

1） M-K检验显示东江年最大流量显著下降,序列呈现非一致性。降雨和流域水库蓄水是造成年最大流量显著下降的主要原因。

2） 龙川和河源站年最大流量LN2混尾分布拟合最好,变化环境后洪水频率最优分布线型基本一致,但流量变小造成分布参数改变已导致分布线型高水尾部特征由“陡峭”变为“平缓”。

3） 随序列参数及分布线型特征演变,指定标准P（百年一遇）下,其相应设计流量值也时刻发生改变,量级由大逐渐减小。

4） 在影响洪水频率的因素中,大型水利蓄水工程对洪水重现期的影响不可忽视。流域径流序列形成的环境背景“一致性”已不复存在,采用现有工程水文分析方法将面临由变化环境带来的设计频率失真风险。

5） 变化环境背景下,洪水重现期往往不是描述一场洪水的一个固定不变的属性。用水文情势变化前估计的洪水重现期往往不能很好地描述变化环境后洪水频率的特征,需进一步采用非一致性水文序列频率计算。

The authors have declared that no competing interests exist.

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子