基于时变矩模型的湄公河上游非一致性洪水频率分析

龙云舒， 叶长青， 李运刚， 熊合勇， 赵承远

(1.云南大学 云南省国际河流与跨境生态安全重点实验室， 云南 昆明 650504； 2.海南大学 生态与环境学院， 海南 海口 570228； 3.华能澜沧江水电股份有限公司， 云南 昆明 650214)

1 研究背景

在水利水电工程设计中，水文频率分析是设计洪水计算及确定重现期的关键[1-2]。传统水文频率计算建立在水文时间序列一致的基础上，近些年来，在气候变化和人类活动的影响下许多河流的水文序列一致性遭到破坏，传统水文频率分析方法的适用性受到质疑[3-5]。变化环境下非一致性水文频率分析已成为水工设计规划、洪水分析和风险管理中的关键技术环节，也是水文科学研究中的热点问题[6-7]。

非一致性洪水频率研究主要集中在水文序列非一致性诊断和非一致性洪水频率分析两个方面。国内外代表性的非一致性检验方法有Mann-Kendall、Pettitt、Spearman、Brown-Forsythe检验和贝叶斯方法等[7-8]。在非一致性水文频率分析方面，常用水文极值序列还原/还现的方法，如降雨径流关系法、时间序列的分解与合成法及水文模型法等[9-10]。这类方法对水文序列进行还原/还现处理后，采用一致性假设的方法对水文序列进行频率计算，然而还原/还现方法的一致性修正成果往往存在不确定性[6-7]。因此，基于非平稳极值序列的直接水文频率分析逐渐受到重视，如混合分布法[11]、条件概率分布法[12]、时变矩模型[13]和GAMLSS(generalized additive models for location, scale and shape)模型[14]等。其中，时变矩模型考虑均值和方差的趋势性，将趋势性成分嵌入到了概率分布的一、二阶矩中(时变矩)，可得到设计值随时间的变化关系。时变矩模型是描述单变量水文序列非一致性的有效工具，已经被应用到许多地区[15-17]。

澜沧江-湄公河流经中国、缅甸、老挝、泰国、柬埔寨和越南6个国家，是东南亚最重要的国际河流。在气候变化和人类活动影响尤其是大规模梯级水电开发背景下，澜沧江-湄公河干流径流特征已经改变，水文序列一致性遭到破坏[8,18]。本研究以干流泰国清盛水文站为例，利用Mann-Kendall趋势突变检测方法对年最大日流量序列进行非一致性分析，基于时变矩模型选择6种概率分布和7种趋势模型进行组合，产生42种竞争模型进行比较择优，分析变化环境下洪水设计值的响应规律，以期为流域水工规划、防洪安全和跨境水资源利用管理等提供科学依据。

2 数据来源与研究方法

2.1 研究区概况及数据来源

澜沧江-湄公河河长4 880 km, 流域面积为81×104km2，多年平均径流量为4 750×108m3(图1)。流域主要受东南季风和西南季风的影响，降水量丰富但季节分配不均，有着较为明显的干湿季。2008年全流域水库有效库容(86×108m3)仅占多年平均径流量的2%，而预计到2025年水库有效库容将增加到868×108m3，达到多年平均径流量的19%[19]。截止目前中国境内澜沧江干流已建成11个梯级大坝，总库容约444.88×108m3(表1)。位于泰国境内的清盛水文站是湄公河干流上离中国最近的水文站，控制面积为18.9×104km2，多年平均流量为2 653 m3/s。

表1 澜沧江干流已建梯级电站情况

图1 澜沧江-湄公河上游流域水系及已建水电站分布

本研究中采用的清盛站1960-2019年逐日流量数据来源于湄公河委员会(https://www.mrcmekong.org/)。

2.2 研究方法

首先对水文序列采用Mann-Kendall法[20]进行趋势分析和突变检验。若序列的一致性遭到破坏，则选择时变矩方法进行分析。时变矩方法主要分析水文频率曲线特征参数随时间变化的影响，认为水文序列的第一阶矩均值(m)和第二阶矩标准差(σ)随时间具有线性或抛物线性趋势变化特征，并引入时间t对水文频率曲线进行函数描述[12]。时变矩模型可由不同的概率分布线型和水文序列的第一、二阶矩的趋势模型相互组合得到。本研究选择P-Ⅲ型分布(P-Ⅲ)、广义极值分布(GEV)、广义逻辑斯谛分布(GLO)、正态分布(NORM)、二参数对数正态分布(LN2)和耿贝尔分布(GMB)共6种概率分布函数。为分析前两阶矩随时间的变化关系，选取合适的趋势模型嵌入到分布模型中。考虑到曲线的外延性，采用线性(L)和抛物线性(P)趋势。将趋势嵌入频率曲线的第一、二阶矩中可得到7类趋势模型，具体描述参考文献[13]、[16]。采用极大似然法进行参数估计，并使用AIC(Akaike information criterion)准则选择最优模型作为时变矩模型进行相应的频率计算。时变矩方法具体的计算流程框架如图2所示。

图2 时变矩方法计算流程框架

3 结果与分析

3.1 洪水系列非一致性诊断

图3为1960-2019年清盛水文站年最大日流量序列变化过程及Mann-Kendall突变检验。由图3(a)可以看出，1960-2019年年最大日流量呈“显著下降-平稳-下降”的变化过程。其中，1966和2006年分别发生了洪量为23 500和29 300 m3/s的特大洪水。自2008年后，年最大日流量虽有波动，但均小于10 000 m3/s。图3(b)中趋势和突变分析结果表明，清盛水文站的年最大日流量呈显著的下降趋势(p<0.01)，突变点出现在2008年。

图3 1960-2019年清盛水文站年最大日流量变化及Mann-Kendall突变检验

3.2 时变矩模型优选

基于时变矩方法，选择6种概率分布线型和7种趋势模型相互组合，共得到42种备选模型。采用极大似然法给出最优线型的参数估计值，并计算各种模型组合的AIC准则值，结果见表2。根据AIC准则，选择AIC值最小的模型组合作为最优模型组合，则清盛水文站年最大日流量序列最优拟合分布为GLO分布，最优拟合趋势模型为AL趋势模型，时变矩最优模型为GLO-AL模型，即均值具有线性趋势而标准差不考虑趋势变化。

表2 清盛水文站时变矩模型AIC拟合检验值

3.3 水库修建不同阶段洪水线型响应规律及重现期变化

3.3.1 不同时间基点下洪水线型响应规律 结合Mann-Kendall法突变点检验并考虑澜沧江梯级水电工程的建设时间(表1)，选择时间基准点分析水电站工程建设前后洪水线型的变化情况。澜沧江干流第1个梯级水电站漫湾水电站于1986年开工，1992年一期工程结束；第2个梯级水电站大朝山水电站于2003年完工投产，其他电站于2008年以后陆续投产运行。因此选择水库修建前的1970、1985年以及1994、2005和2019年作为时变矩模型的基准时间点。

图4为清盛水文站不同基准时间点的洪水拟合线型变化情况。在尚未进行梯级水电开发时期即1960-1970年， GLO-AL线型高水尾端位于所有线型的最上方，同一量级洪水出现的概率偏大。1971-1985年GLO-AL曲线下移，但下移的幅度不大。在梯级水电站建成影响期，1986-1994年(漫湾水电站运行)GLO-AL曲线有所下降；1995-2005年(漫湾和大朝山水电站联合运行)GLO-AL曲线相较于1994年又有所下降；2006-2019年(11个水电站陆续运行)GLO-AL曲线位于所有时间基点线型的下方，同量级洪水发生概率最低。澜沧江中下游以小湾、糯扎渡水电站为核心的梯级水电开发对洪水线型的影响较小，但水库联合运行以后清盛站同量级洪水发生的概率降低，设计洪水量级也明显减小。这说明澜沧江水电站的建设对清盛水文站洪水具有明显的调控作用。

图4 清盛水文站不同基准时间点的洪水拟合线型变化

3.3.2 设计洪峰流量与重现期的变化过程 用传统频率分析方法计算得到清盛水文站100年一遇洪水的设计值为22 881 m3/s，用时变矩模型分析该设计值下重现期的变化，如图5所示。由图5可见，与传统频率分析下某一指定流量标准值重现期不变相比，采用时变矩模型得到的重现期随时间呈上升趋势。如1960-1965年为小于100年一遇，1966-2011年为大于100年一遇且小于200年一遇，2012年后为大于200年一遇。选取指定重现期T=100 a，采用时变矩模型分析指定重现期下设计流量的变化，如图5所示。

图5 清盛水文站设计流量下重现期及指定重现期下设计流量变化

由图5可见，与传统频率分析下指定重现期标准下的设计流量值不变相比，采用时变矩方法计算得到的设计流量随时间呈减小趋势。清盛水文站100年一遇洪水量级从1960年的23 274 m3/s减小至2019年的19 085 m3/s。如果不考虑序列非一致性处理，得到的100年一遇设计值为22 881 m3/s，该值比时变矩模型求得的2014-2019年糯扎渡水电站运行后的设计值高出18%～20%，即传统频率分析方法会高估洪水量级。从时变矩方法计算得到的设计流量来看，澜沧江-湄公河上游梯级水库尤其是糯扎渡水库运行后明显减小了清盛水文站100年一遇的洪峰流量，对下游防洪具有积极作用。

4 讨 论

气候和人类活动是河川径流变化的两大主要驱动因素。基于站点观测数据分析表明，近50年来澜沧江-湄公河上游地区降水量减少趋势不明显，气温升高和潜在蒸散发增加是径流量显著减少的重要因素[21-22]。相对于气候变化，人类活动对径流的影响在特定时期内则更为明显。研究表明，澜沧江允景洪水文站1987-2014年多年年均径流量较1980-1986年减少了6%，人类活动对径流量变化的贡献率从1987-2007年的43%增加到2008-2014年的95%，人类活动对径流量变化的影响逐渐加剧[23]。20世纪80年代以来，澜沧江梯级水库建设导致径流年内分配呈现均匀化，下游清盛站年最大日流量呈显著的减少趋势[24-25]。总的来说，变化环境尤其是梯级水库的建设对允景洪-万象区间径流年内分配的影响显著[26]。

本研究结果表明，1960-2019年来清盛水文站年最大日流量序列呈现出显著减小趋势，水文极值序列突变年份(2008年)与景洪水电站蓄水、小湾水电站投产时间吻合，说明变化环境下水文极值序列发生了变异。未来湄公河上游地区极端降水量将趋于增多且降水集中程度将增大，流域洪水风险势必增加[27]。流域内水库联合调度可以在一定程度上缓解洪水的风险[28-29]，如全流域水库联合调度可以将湄公河干流洪水从200年一遇减至20～50年一遇[30]。在未来气候变化和人类活动加剧的背景下，使用传统频率分析方法得到的重现期和洪水设计值将会“失真”。本研究采用时变矩方法分析指定设计流量下重现期和指定重现期下设计流量，发现重现期和设计流量均随时间变化。在非一致背景下，如果采用传统水文频率分析方法而不对序列的非一致性进行考虑，将会错估洪水重现期及设计洪水量级，因此传统洪水频率计算方法的使用应该慎重。

5 结 论

(1)1960-2019年清盛水文站年最大日流量序列呈显著减小的趋势并在2008年发生突变，水文极值序列呈现非一致性。

(2)采用时变矩方法选择6种概率分布和7种趋势模型进行组合，得到42种竞争模型进行比较，结果表明清盛水文站广义逻辑斯谛分布搭配均值具有线性趋势的模型拟合最优。

(3)基于传统频率分析方法求得的100年一遇洪水设计值，在时变矩模型下重现期由水电开发前的小于100年一遇变化到2012年后的大于200年一遇；100年一遇洪水设计值也随时间呈减小趋势。

(4)湄公河上游梯级水库尤其是糯扎渡水库运行后明显减小了清盛站100年一遇的洪峰流量，对下游防洪具有积极作用。