考虑非平稳性特征的雅砻江流域洪水量化分析

高 洁，杨 龙

(1.水电水利规划设计总院，北京 100120；2.南京大学地理与海洋科学学院，江苏 南京 210023)

0 引 言

设计洪水是确定水电工程规模和制定防洪标准的重要基础。基于频率分析的设计洪水计算在分布线型、参数估计等方面积累了丰富的研究成果[1]。气候变化和人类活动对水文时间序列的平稳性特征产生影响，传统水文设计方法面临变化环境带来的挑战[2]。郭生练等[1]将非一致性洪水频率分析、设计洪水不确定性研究、基于水文物理机制的洪水频率分析等内容列为当前设计洪水研究前沿与热点。鲁帆等[3]在水文气象领域非平稳时间序列极值模型研究综述中介绍，传统水文频率分析通常采用皮尔逊III型(P-III)，对数皮尔逊III型，广义极值分布(GEV)，广义帕累托分布(GPD)等理论线型。我国工程界设计洪水标准主要采用P-III分布。以年最大为样本的GEV分布也在水文极值统计中逐渐被采用。关于GEV分布参数，GEV-CDN模型[4]认为形状参数决定了分布的尾部形态，具有较大的不确定性，可将形状参数取值设定在[-0.5，0.5]。GEV模型应用于洞庭湖区、岷江流域年最大日降水系列[5- 6]，发现各地的形状参数差异较大，形状参数为0的极值I型Gumbel分布、形状参数<0的极值III型Weibull分布，形状参数>0的极值II型Frechet分布均有适用。通过对全国降雨资料统计，Yang等[7]发现极端降雨在长江流域呈现形状参数<0的薄尾分布；在环渤海地区形状参数>0，呈现厚尾分布。长江中下游控制站大通水文站[8]1950年至2012年历年最大和非汛期最大流量序列符合形状参数<0的Weibull分布，新疆天山开都河[9]1957年至2011年历年、春季和夏季最大日流量序列服从形状参数≥0的Gumbel分布和Frechet分布。通过对珠江流域北江水系武江上游坪石水文站和下游犁市(二)水文站实测年最大径流系列分别采用薄尾P-III分布、Weibull分布、厚尾广义逻辑分布(GLO)、三参数混合尾GEV、Gumbel、两参数对数正态分布(LN2)[10]拟合发现，厚尾GLO和三参数混合尾GEV分布在拟合时可使高水端上翘，提高了高水端拟合优度。

非平稳时间序列的概率分布均值、方差或协方差随时间变化，重现期和风险度有别于传统水文频率计算结果，需要对变化序列进行还原/还现，或针对时变序列进行频率计算。谢平等[11]按照“分离”“合成”的思路对水文序列非一致性进行诊断。考虑时间因素对传统频率分布的影响，刘德地等[12]推导了非平稳条件下洪水重现期的计算公式。顾西辉等[13]选取珠江流域28个水文站1951年至2010年历年最大流量数据，在突变点分析的基础上采用GEV、GPD、Gumbel、P-III、GLO、广义对数正态分布(GNO)6种分布进行拟合。研究发现，最优分布以GEV居多。珠江流域内西江、北江频率设计值增大，且随着重现期增加而增大，基于整体样本的水利水电工程规划设计，可能低估大洪水影响，导致防洪风险。顾西辉等[14]又将GEV-CDN模型[4]应用于珠江流域6个水文站近45年历年最大日流量序列。研究发现，50年一遇以下洪水，平稳和非平稳模型成果差异不大；但50年一遇以上洪水，不同模型频率设计值差异较大，存在潜在的防洪风险。梁忠民等[15]针对非平稳环境下工程水文设计值发生变化的问题，以等可靠度建立平稳序列与非平稳序列的联系，提出了待建工程、已建工程推求和协调频率设计值的思路。

考虑位置参数、尺度参数、形状参数的广义可变模型GAMLSS为变化环境水文频率计算提供了直接可靠的流程和方法。顾西辉等[16-17]将GAMLSS模型应用于全国839个气象站极端降水发生率及新疆塔里木河洪水发生率的非平稳性研究。高洁[18-19]将GAMLSS模型应用于雅砻江流域。研究发现，全流域降水特性无显著变化，但是流域中下游极端降水量增加，并伴随着最长连续干期变长、连续湿期变短。流域中上游极端降水平稳性较好，下游不同站点极端降水非平稳性特征差异较大。小流域及支流洪水的平稳性特征易受局地极端降水影响。虽然有建坝修库等人类活动，但随着支流入汇，流域干流水文站洪峰流量回复平稳状态。

关于水文时间序列，基于平稳假设和考虑非平稳性，均开展了大量的研究工作。根据传统规范要求，我国水电水利工程行业，水文设计选择P-III分布作为频率计算的基础。考虑极值理论和变化环境，传统频率计算成果与GEV分布、GAMLSS模型成果的差异尚缺乏定量分析。本文以雅砻江流域为研究对象，通过GEV统计参数和GAMLSS最优模拟，量化分析GEV分布、GAMLSS模型成果与传统P-III适线成果差异。雅砻江干流河道天然落差约3 870 m，水能资源丰富，是我国规划建设的十四大水电基地之一。根据河流水电规划，雅砻江流域水电技术可开发量约2 900万kW，主要集中在流域中下游河段，分布着众多水电水利工程。因此，充分认识流域洪水时变特征、洪水频率计算成果偏差，对于流域内众多涉水工程设计运行具有重要意义。

1 数据和方法

1.1 研究对象

雅砻江流域(96°52′E～102°48′E，26°32′N～33°58′N)是长江上游金沙江最大的一级支流，地处青藏高原东南部，位于金沙江和大渡河之间的狭长地带。流域面积约13万km2，谷岭高差大[20]，地形地貌以高山峡谷为主，基于90 m空间分辨率的DEM数据，海拔高程位于966～5 888 m[20]，地势西北高东南低。本研究收集了流域内8个水文站：甘孜、雅江、泸宁、道孚、朱巴、濯桑、庄房(二)、孙水关长系列年最大洪峰流量系列。

1.2 研究方法

1.2.1洪水系列趋势检验和突变分析

针对洪水系列的单调变化和突变特征，分别进行趋势检验和突变点分析。

MK检验作为广泛应用于水文气象时间序列趋势分析的非参数化方法，得到了国际气象组织的推荐[21]，被广泛应用。在突变点分析方面，分别采用Pettitt、AMOC和ECP三种检验方法。Pettitt法通过检验时间序列均值变化，确定序列跃变发生时间[22]。基于R语言平台，changepoint程序包[23]的AMOC检验法是基于似然函数框架的突变点检测，时间序列无需满足正态分布假设。ECP程序包[24]是基于凝聚和分裂算法的聚类分析法，可用于时间序列的多突变点检验。

1.2.2洪水系列非平稳性分析

Stasinopoulos和Rigby于2005年提出(半)参数回归模型GAMLSS(基于位置参数、尺度参数和形状参数的广义可加模型)，该模型[25]可模拟统计参数随时间的变化模式。本文为了简化模型结构，采用两参数分布函数，基于对数正态(LOGNO)、Gamma(GA)、Gumbel(GU)和Weibull(WEI)分布，拟定位置参数、尺度参数与时间的关系，通过AIC准则选择最优分布函数，计算水文序列的时变频率设计值。

1.2.3洪水系列适线

我国传统设计洪水频率适线通常采用P-III型曲线，极值统计近年来也得到了广泛关注[26]。GEV作为极值统计的一个重要分类，包含了Gumbel、Frechet和Weibull三种分布类型。

H(x;μ，σ，ξ)=exp{-[1+ξ(x-μ)/σ]-1/ξ}，[1+ξ(x-μ)]/σ>0

(1)

式中，μ为位置参数；σ为尺度参数；ξ为形状参数。μ和σ反映洪水的量级和变幅。ξ决定了分布函数的尾部形态[7]。ξ>0的极值II型(Frechet)呈现厚尾分布，尾部较长且无上边界。ξ<0极值的III型分布(Weibull)具有上端点μ-σ/ξ，呈现薄尾分布。当统计变量x趋近于无穷，ξ=0的极值I型分布(Gumbel)为无上下边界的薄尾分布。

表1 各水文站趋势检验和突变点分析

2 成果分析

2.1 时间序列特性

根据MK趋势检验，各水文站的年最大洪峰流量序列均无显著变化趋势(显著度水平5%)。由于收集的洪峰流量系列长度仅在37年至65年，本研究检测一个突变点，并根据变点位置将序列划分为前后两部分。经Pettitt、AMOC(mean)和ECP突变检验，甘孜、濯桑、庄房(二)、孙水关4站通过显著水平5%突变检验，雅江站满足显著水平10%突变检验。变化趋势和突变时间点如表1所示。

2.2 平稳性分析

对各水文站年最大洪峰流量系列进行GAMLSS非平稳性分析，在LOGNO、GA、GU和WEI四种分布中，以LOGNO分布占优。虽然各水文站年最大洪峰流量的MK趋势检验无显著变化，但是与甘孜、濯桑、庄房、孙水关、雅江5站通过突变检验相协调。GAMLSS模型显示，其中的甘孜、濯桑和孙水关3站的年最大洪峰流量呈现非平稳特征，如图1所示。

(1)甘孜站位置参数表征的年最大洪峰流量均值随时间单调减小，MK检验呈不显著减小趋势，1986年以后均值突变减小。

(2)濯桑站尺度参数表征的不同重现期年最大流量离散程度逐渐减小，体现为常遇频率设计值增加、稀遇频率设计值减小，并逐渐稳定。MK检验常遇频率设计值呈不显著增加。

(3)孙水关站位置参数表征的年最大洪峰流量均值在1980年～1990年出现拐点，先减小后增大。MK检验常遇频率设计值呈不显著增加，在1985年以后突变增加。

受限于观测资料多为常遇频率，MK趋势检验和突变分析主要反映的是常遇频率下洪水事件。GAMLSS成果可以体现不同重现期表征的常遇频率和稀遇频率设计值差异化的时变特征。年最大洪峰流量的非平稳性分析(见表2)，基本映射出渐变性和突变性检验中尚不显著的发展趋势，具有一定的前瞻性。

时间序列趋势检验无显著变化，采用传统的P-III适线和稳态GEV分布拟合仍然合适。基于GAMLSS非平稳性分析和考虑突变点的GEV分布，可反映变化环境下洪峰流量频率设计值随时间的变化。

表2 各水文站GAMLSS非平稳性分析成果

2.3 频率适线成果对比

P-III分布目估适线，通常微调Cv和Cs两个参数，控制适线的倾斜程度和弯曲程度，使高水端上翘，倾向于头部特大值，使稀遇频率设计偏安全，得到比较保守的频率设计值。GAMLSS两参数分布由于参数较少，灵活性低，不能人为提高高水端的拟合优度，可能出现低频稀遇设计洪水洪峰流量小于实测头部特大值，该方法调整裕度较少。GEV分布通过拟合位置参数、尺度参数和形状参数计算频率设计值。位置参数和尺度参数均与集水面积具有相关性，关于形状参数的影响因素，学术界尚无定论。但是，形状参数又对于GEV尾部形态具有关键作用，以正负为界，区别厚尾和薄尾特征，体现了分布极值的偏离程度和密集程度，对于研究极端事件具有重要意义。

2.3.1平稳序列适线成果对比

年最大洪峰流量序列相对平稳的雅江、泸宁、朱巴、道孚、庄房5站经三参数P-III目估适线、GEV分布拟合及两参数GAMLSS模型优选，除了道孚站在低频稀遇洪水时GEV分布推算频率设计值偏大，其他4站均为GEV分布和GAMLSS模型优选成果较为接近，普遍小于P-III适线成果，且部分站点(如泸宁、庄房)随着重现期增大，稀遇洪水设计值的偏差增加。三种方法成果对比见图2。

图2 平稳序列GAMLSS、P-III、GEV 3种方法计算成果比较

根据平稳序列5个水文站的拟合成果，三种方法对于20年一遇以下洪水洪峰流量的拟合偏差低于10%，50年一遇至200年一遇洪水洪峰流量的拟合偏差≤15%，1 000年一遇以下洪水洪峰流量的拟合误差基本不超过20%。对于5 000年、10 000年等稀遇特大洪水的拟合偏差在25%至30%之间(见图2f)。

2.3.2非平稳序列适线成果对比

通过GAMLSS模型推算年最大洪峰流量时间序列非平稳的甘孜、濯桑和孙水关3站频率设计值的时变特性。随着重现期增大，频率设计值变幅增大，体现了采用有限资料推算稀遇事件的外延性导致不确定性增加(见图3a、3b、3c)，可通过频率设计值时变范围与设计值的比例来反映。

甘孜站1956年～2016年洪峰流量频率设计值时变范围与均值比例，对于不同重现期基本在14.8%～15.5%。濯桑站1966年～2016年洪峰流量频率设计值随着重新期增大，时变范围比例从10年一遇洪水23%增至10 000年一遇洪水的57%。孙水关站时变范围比例从10年一遇洪水45%降至10 000年一遇洪水的39%。

濯桑站和孙水关站由于集水面积较小，实测年最大洪峰流量量级相对较小(均值分别为287 m3/s和610 m3/s)；因此，频率设计值时变范围比例较大(23%～57%)。GAMLSS频率设计值的时变特征是否具有尺度效应，小流域非平稳性普遍较大，有待更多站点验证。

随着重现期增大，甘孜站和孙水关站洪峰流量时变范围与均值的比例基本稳定。濯桑站洪峰流量时变范围比例明显增大，主要由于低频稀遇洪水时变减小趋势(见图1b、3b)，体现为10年一遇以上洪水随着重现期增大增幅减小，稀遇洪水致灾风险随时间降低。

P-III目估适线成果均在GAMSS测算区间内。但是，GEV分布在低频大值的估算上与P-III适线和GAMLSS模型成果(见图3a和3c)产生偏差。GEV分布对于非平稳系列的拟合较为敏感，需要考虑参数时变特性和尾部特点。

图3 非平稳序列GAMLSS、P-III、GEV三种方法计算成果比较

表3 各水文站突变前后及完整系列GEV参数

2.4 考虑突变的GEV变化分析

2.4.1GEV参数分析

Yang[7]发现中国极端降水GEV分布的位置参数、尺度参数与集水面积相关，长江流域极端降水的GEV形状参数以负值居多。与该结论一致，流域内8个水文站原始系列、突变前后系列的位置参数和尺度参数均与集水面积(见表3)呈现双对数相关(见图4)。本流域位于长江流域上游支流，年最大洪峰流量GEV分布也以形状参数<0的薄尾分布为主。仅道孚站和孙水关站的形状参数>0，呈现厚尾分布，由此造成两站GEV拟合频率设计值比P-III分布和GAMLSS上限值均偏大(见图2c、3c)。GEV对于极大值的过分高估，是否由于形状参数>0所致，值得研究。

图4 GEV统计参数与集水面积双对数线性相关

2.4.2突变前后GEV分析

选取形状参数<0的甘孜、雅江、濯桑、庄房4个站点，基于突变点分析，对比完整系列、突变前、突变后3种系列的GEV位置参数、尺度参数和形状参数，推算频率设计值。

(1)甘孜站和濯桑站分别在1986年突变后和1989年突变前，系列的GEV形状参数>0，由此导致GEV拟合值偏离较大，分析无效。

(2)雅江站在突变点1998年以后，各频率设计洪水均有增加(见图5a)。因此1998年以前设计的工程有必要复核现状防洪安全。

(3)庄房站在2006年以前，各频率设计洪水略大于现状水平，2006年以前的设计仍可满足现状相应的防洪标准(见图5b)。

图5 基于完整系列、突变前后系列推算设计洪水频率值

3 结 论

基于MK趋势检验，Pettitt、AMOC和ECP突变分析，研究雅砻江流域8个水文站年最大洪峰流量渐变和突变特征，基于GAMLSS模型分析流域洪水的平稳性特征。按照平稳序列和非平稳序列，对各水文站年最大洪峰流量系列进行P-III目估适线、GEV分布拟合和GAMLSS模型成果对比分析，探索GEV分布位置参数、尺度参数和形状参数特点，结合突变点分析，研究流域防洪安全。研究发现：

(1)GAMLSS拟合以对数正态最优分布居多。其中，甘孜站、濯桑站和孙水关站年最大洪峰流量时间序列非平稳。GAMLSS模型成果与实测系列渐变、突变分析统筹考虑，体现了不同重现期设计值时变特征，揭示了实测点难以反映的低频洪水特性。

(2)比较了3参数P-III目估适线、GEV分布和两参数GAMLSS模型成果，针对平稳序列5个站点，倾向于头部特大值的P-III目估适线成果普遍大于以穿过点据中央为拟合原则的GEV分布和GAMLSS模拟成果。3种途径推算20年一遇频率洪水偏差<10%，200年一遇频率洪水偏差<15%，1 000年一遇频率洪水偏差基本不超过20%。

(3)非平稳序列的3个站点，P-III适线成果基本在GAMLSS模型拟合时变范围内。GEV拟合非平稳序列建议考虑参数的时变特性和尾部特点。

(4)雅砻江流域年最大洪峰流量序列GEV分布形状参数以<0的薄尾分布为主。呈现厚尾分布的道孚站和孙水关站GEV频率设计值远大于其他模型成果，存在极大值被高估的现象。雅江站1998年以前设计的工程有必要复核现状防洪安全。