【www.zhangdahai.com--其他范文】
何朝飞，王晓云，陈伏龙，唐 豪，龙爱华,3
(1.石河子大学水利建筑工程学院，新疆 石河子 832000；
2.中水北方勘测设计研究有限责任公司，天津 300222；
3.中国水利水电科学研究院流域水循环模拟与调控国家重点实验室，北京 100038)
融雪洪水作为西北干旱区一种典型的自然灾害事件，严重危害了人们的人身财产安全，如何控制洪水的发生并有效降低财产损失是灾害风险研究的重点。在随机性洪水事件中，频率分析是评估其风险的有效手段。近年来，随着全球气候变化及人类活动的影响，洪水序列的随机独立同分布假设面临极大的挑战[1]。如果直接利用实测洪水资料进行分布拟合会大大降低其结果的可靠性，因此必须进行序列变异诊断，并对已经发生变异的洪水序列进行一致性处理[2]。
随着水文频率分析研究的不断深入，非一致性频率分析方法及基于Copula函数的多变量组合分析方法得到了很好的应用研究。目前水文非一致性频率分析方法主要从两个角度进行，一是基于水文序列还原、还现的间接分析方法，即将水文序列修正到过去或现在的状态；
二是非平稳极值序列的直接分析[3-4]。Singh等[5]针对洪水序列的非一致性问题首次提出了混合分布模型，并有效应用于非一致性洪水序列频率分析中。该分布模型既能避免水文时间序列修正和分解等转化的烦琐步骤，又能克服人工适线法带来的较大误差，因此得到了较为广泛的应用[6-9]。郑锦涛等[10]在多种变异诊断方法的基础上，对玛纳斯河年径流量进行了混合分布模型和条件概率分布模型的频率分析，很好地描述了混合分布模型的拟合效果。
为解决单变量频率分析与风险评估方法不足以描述极端事件的失效概率和重现时间的问题[11-12]，陈晶等[13]将Copula函数引入到水文学中，并阐述了其在构建水文多变量联合分布方面的相关问题；
黄强等[14]在多变量分布模型的基础上，提出了二次重现期概念，并分析探讨了“且”“或”以及二次重现期在洪水风险分析中的可靠性；
高玉琴等[15]在G-H Copula函数的基础上对比分析了洪水单变量风险概率及二维、三维变量联合风险概率，发现它们之间存在明显差异。但以上联合分布研究均建立在水文序列一致的基础上，缺少对于非一致性条件的必要考虑。
现阶段对于干旱区融雪洪水峰量组合风险规律的研究仅以单一变量为主要参考值，其结果往往较为片面，且在联合分析中缺乏对变量间组合的风险量化，风险指标过于单一，不能较全面地反映其设计水平值。因此，本文选取玛纳斯河流域出山口肯斯瓦特水文站1957—2014年逐日径流量，并根据该站融雪型洪水一日一峰的典型特点，建立年最大洪峰和年最大1 d洪量的设计频率分析样本，结合序列变异诊断结果，选用不同子序列频率分布构建的混合分布为联合分布模型的边际函数。由联合函数特性分析及数据结构特征，将Gumbel Copula函数作为联合函数进行建模，并以流域500年一遇设计防洪标准为例，多角度分析组合风险类型，以期为玛纳斯河流域防洪工程设计提供参考。
玛纳斯河发源于新疆天山北麓的依连哈比尔尕山，东西长198.7 km，南北宽260.8 km，流域总面积为24 300 km2，平均气温为4.7～5.7℃。流域地处中纬度欧亚大陆中心带(85°01′E～86°32′E，43°27′N～45°21′N)，天山北麓，准噶尔盆地南部，具有明显的大陆性干旱气候特征。河流属于冰雪融水和雨水混合补给型，平均年径流总量约12.99亿m3，其中海拔3 600 m以上终年积雪覆盖，平均雪线高度为3 970 m，冰川面积多年维持在608.25 km2左右，占玛纳斯河径流补给量的35.3%[16]，是河流的主要补给源。本文研究对象为玛纳斯河肯斯瓦特水文站(85°57′E，43°58′N)控制流域(图1)，该水文站位于玛纳斯河出山口处，是干支流汇合后的主要控制点，控制面积为4 637 km2，它的建立为玛纳斯河流域融雪洪水观测及预报起到了重要的作用。
图1 肯斯瓦特水文站控制区域
2.1 累积距平法与滑动t检验
(1)
由式(1)可得n个时刻的累积距平值，其中序列突变年为绝对值最大指标所对应的年份。
滑动t检验是根据两组样本平均值的差异是否显著来判断其突变情况[18]。对于具有n个样本量的时间序列X，基于某一基准点前后两子序列X1和X2的样本分别为n1和n2，则统计量为
(2)
2.2 混合分布模型
混合分布模型可直接用于非一致性时间序列的频率分析。假设该水文时间序列是由若干个一致性子分布赋予一定的权重系数累加而成：
F(x)=α1F1(x)+α2F2(x)+…+αkFk(x) (3)
式中：F1(x)、F2(x)、…、Fk(x)为各子分布的累积分布函数；
α1、α2、…、αk为各子分布的权重，且各权重之和为1。
为保证模型参数估计的准确性，根据洪水形成机制，对洪水时间序列进行子序列的合理划分，使子序列的个数保持在最低限度，避免了模型的过度拟合。当非一致性融雪洪水特征时间序列X的样本容量为n，变异点为τ时，变异点之前的时间序列为F1，子序列服从概率密度函数为f1(x)的分布；
变异点之后的时间序列为F2，子序列服从概率密度函数为f2(x)的分布；
整体融雪洪水时间特征序列X服从概率密度函数为f(x)的混合分布，其表达式为
(4)
式中α′、1-α′为两个一致性子序列的权重系数。
2.3 Copula函数及尾部相关关系
近年来，随着水文研究工作的快速发展，多变量水文联合分析得到了深入的研究，其中，Copula函数作为连接或耦合两个或多个独立变量的数学函数[19]，得以普遍的应用。将一个k元联合分布函数用多个边缘分布和一个Copula函数连接起来构建多维分布，依次用来描述变量间的相关关系。本文利用二元Copula函数建立模型，其表达形式为
C(u,v)=φ-1[φ(u)+φ(v)]
(5)
式中：u、v分别为随机变量；
φ为连续、严格递减的凸函数，满足φ(0)=∞和φ(1)=0；
φ-1为φ的反函数，且满足φ-1(∞)=0和φ-1(0)=1。
常见的单参数Archimedean Copula函数有Clayton Copula、Frank Copula和Gumbel Copula。Clayton Copula函数对变量分布的下尾部变化十分敏感，能够快速捕捉到其下尾的变化，Gumbel Copula函数则对变量分布的上尾部变化十分敏感，能够快速捕捉到其上尾的变化。由于Frank Copula函数的密度函数具有对称性，无法捕捉随机变量间非对称的相关关系[20]，因此在洪水等极端事件中局限性较大。Archimedean族系的双参数Copula函数可以同时刻画上尾和下尾的相关性，因此本文采用Clayton Copula和Gumbel Copula两种单参数Copula函数以及BB1 Copula和BB5 Copula两种双参数Copula函数[21]进行分析。
2.4 组合分析计算
洪水的发生是由多种影响因素共同决定的，各变量间的相互作用也复杂多变。因此分析不同条件下的洪峰和洪量组合发生的概率，即峰量组合风险分析是极其有必要的。本文分别选取联合、同现、组合3种风险模型[22-24]，以多角度刻画不同量级峰量的风险概率来评价洪水设计值，其风险概率模型见表1。
表1 风险概率模型
表1中Q和W1分别为年最大洪峰和年最大1 d洪量； 1957—2014年玛纳斯河肯斯瓦特水文站控制流域年最大洪峰和最大1 d洪量的累积距平曲线如图2所示，可以很直观地看出峰量都具有先下降后上升的变化趋势，并在1993年出现了较大的转折。 图2 峰量累积距平变化 为了进一步确定洪水时间序列的突变年份，利用滑动t检验，取步长n1=n2=10，显著性水平α=0.01，结果如图3所示，其中年最大洪峰在1992年和1993年发生突变，年最大1 d洪量在1971—1975年、1992年和1993年发生突变。 图3 滑动t检验统计量变化 IPCC第五次评估报告显示，1980—2012年全球气温升高了0.85℃[25]。由于受冰川补给影响，从20世纪70年代至本世纪初，玛纳斯河流域冰川退缩面积为159.02 km2。在1990年前后，肯斯瓦特水文站控制流域内下垫面变化不大。因此，气温变化对玛纳斯河流域洪水时间序列变化产生了重要影响[26-27]。综合洪水特征分析以及滑动t检验，最终确定年最大洪峰序列和最大1 d洪量序列的突变年份发生在1993年，这与陈伏龙等[28]的研究结论一致，故将玛纳斯河肯斯瓦特控制流域洪水变化过程分为1957—1993年和1994—2014年两个阶段。 对于洪水变量的分布方式，本文分别就极端分布中的广义极值(generalized extreme value，GEV)分布，广义帕累托(generalized Pareto，GP)分布，对数逻辑斯谛克(log-logistic，LLG)分布和对数正态(log-normal，LOGNO)分布方式进行分析比较，其中前两个为三参数概率分布函数，后两个为两参数概率分布函数。分布函数参数采用最大似然法进行估计[29]，并以Kolmogorov-Smirnov[30](K-S)检验进行拟合优选。图4为各分布函数拟合优度D，经比较，最大洪峰1957—1993年和1994—2014年两个子序列Q1和Q2分别采用GEV分布和LLG分布； 图4 各分布函数拟合优度 对于混合分布模型，以频率离差绝对值和(ABS)最小为目标函数，采用模拟退火算法[6]对参数进行估计，其执行策略为：在矩法初步估计的基础上，给定一个最优初始点X0开始探测整个解空间，通过扰动该解产生一个新解XN，按照Metropolis准则判定是否接受新解，相应地降低控制温度，直到满足收敛准则，即没有新解产生或者控制参数小到一定程度[31]。表2给出了混合分布模型的参数估计结果，根据式(3)计算洪水时间序列混合分布理论频率。 表2 混合分布模型参数估计结果 混合分布和未考虑其变异特性的整体时间序列的边缘分布函数的拟合曲线如图5所示。就年最大洪峰而言，混合分布拟合曲线的上部与单一拟合函数GEV和LLG存在较大差异，但混合分布拟合曲线上部更贴近洪水实测值能较好地描述洪水时间序列； 图5 边缘分布概率拟合 考虑到变量间的相互依存关系，利用Pearson相关分析计算年最大洪峰和年最大1 d洪量的相关系数为0.846，由此可见，玛纳斯河流域肯斯瓦特水文站控制流域的年最大洪峰和年最大1 d洪量存在着较强的正相关性。利用贝叶斯框架中的马尔可夫链蒙特卡罗(Markov chain Monte Carlo，MCMC)[32-33]模拟来估计其参数值。同时采用AIC、贝叶斯信息准则(BIC)和Max-Likelihood(MAX-L)方法对联合概率分布拟合优度进行检验，其结果见表3。 表3 Copula函数参数估计及拟合优度 本文采用图形分析法直观地描述函数的拟合优劣程度。将理论联合概率值与经验联合概率值点绘到同一张图中，当点集越靠近45°线，表明联合分布模型拟合较好。其中两变量经验联合概率分布计算公式为 Femp(x,y)=P(X≤x,Y≤y)= (6) 式中：Femp为联合经验概率； 图6 经验联合分布与理论联合分布Q-Q图 Gumbel Copula函数与两参数函数BB1 Copula和BB5 Copula的理论与联合概率点均落在45°对角线附近，起到了较为满意的拟合效果； 分别取重现期10 a、20 a、50 a、100 a和500 a的年最大洪峰与年最大1 d洪量相互组合，并计算联合分布概率，根据表1风险概率模型计算其组合风险概率，以此确定任一重现期标准的设计洪峰和设计最大1 d洪量组合的风险情况。各量级洪峰和年最大1 d洪量组合的风险概率计算结果见表4。 表4 各量级峰量组合风险概率 表4中，当洪水峰量组合为10年一遇时，即洪峰(565.882 m3/s)与年最大1 d洪量(0.369亿m3)在相同低重现期的情况下联合风险(10.89%)与同现风险(8.45%)最高，说明这种组合设计标准最为危险，同时也反映出该组合下防洪设计标准设置过低。在洪峰防洪设计标准一定的情况下，只要稍微提高洪量设计标准，则3项风险指标均呈显著下降趋势。对于500年一遇峰量两变量防洪设计标准，联合风险与同现风险均达到最小值，反映这种组合的设计标准相对最为安全，因此在高防洪设计标准下，洪水峰量共同影响下的洪水发生概率较低。 本文组合风险概率是在“以洪峰为主”设计工况中洪峰变量Q没有超过设计值q，而年最大1 d洪量变量W1超过设计值w1的概率。在500年一遇洪峰组合10年一遇洪量中，由于其年最大1 d洪量的设计标准偏低，故组合风险增加，风险概率达到最大值9.60%； 以流域500年一遇设计防洪标准为例，洪水峰量联合重现期为322 a，联合风险概率为0.31%，即洪峰、洪量中至少有一个变量超过设定值时的概率为0.31%。其同现风险概率为0.21%，即同时遭遇超过设计标准洪峰及洪量的重现期为476 a。组合风险概率为0.01%，即洪峰量级未超过设计标准，但年最大1 d洪量超过设计标准时的概率为0.01%。这表明在利用水库对洪量进行调节来分担风险效果的同时，也不能忽略洪峰对水工建筑工程的不利影响。根据肯斯瓦特水利枢纽工程实际防洪需求，对多个洪水变量间的相关作用进行分析，能更准确地掌握洪水防范效果，为玛纳斯河流域风险评估提供参考。 a.玛纳斯河洪水峰量时间序列在1993年发生了显著性跳跃变化； b.对于玛纳斯河流域，以高重现期参照下的峰量两变量为防洪设计标准，联合风险与同现风险均能接近最小值，即在较高的防洪设计标准约束下，峰量共同影响的洪水发生概率较低。 c.根据实际工程需求，选择频率较为接近的洪水要素设计标准，可降低其组合风险概率，满足低成本承受较大风险的实际需求，考虑洪水多属性的联合特征可为防洪工程设计提供更为全面可靠的理论依据。
q和w1分别为年最大洪峰和年最大1 d洪量的设计值；
u=FQ(q)为年最大洪峰边缘分布函数；
v=FW1(w1)为年最大1 d洪量边缘分布函数；
C(u,v)=C(FQ(q),FW1(w1))为设计年最大洪峰和设计年最大1 d洪量的联合分布函数。P(W1>w1∪Q>q)为Q与W1中至少有一个变量超过设定值的概率；
P(W1>w1∩Q>q)为Q与W1两变量同时超过设计值的概率；
P(W1>w1|Q3.1 洪水序列特征
3.2 混合分布拟合
年最大1 d洪量1957—1993年和1994—2014年两个子序列W1-1和W1-2分别采用GP分布和LLG分布拟合。
年最大1 d洪量模拟中混合分布曲线的下部与单一拟合函数GP为一致，而其上部分布存在一定的差异。同时拟合曲线中部单一LLG分布拟合曲线有更大概率通过洪水实测值，但对于洪水极端事件，线性上尾分布是设计洪水拟合效果更佳的选择。由于混合分布模型考虑了洪水序列的非一致性，采用两条不同分布的子序列进行拟合分布比传统单一分布经验点据拟合效果更好，且更好地保证了洪水设计的安全需求。通过采用K-S拟合优度检验，确定统计数据是否服从于5%显著性水平的拟合分布，得出年最大洪峰与年最大1 d洪量混合分布对应的D值分别为0.054和0.062，均小于临界值D58,0.05=0.179，即混合分布函数均满足假设检验，因此考虑变异情况下的混合分布拟合曲线是可行的。
3.3 联合分布及组合风险
ngk为满足X≤x、Y≤y的观测值，其中g和k分别为两实测数据从小到大排列后的序号；
n为系列长度。其概率散点Q-Q图见图6。
而Clayton Copula函数高频率概率点偏差较大，拟合效果并不理想(图6)。由表3可以得出：Gumbel Copula函数拟合年最大洪峰与年最大1 d洪量效果最好，其中AIC值为-433.963，BIC值为-431.902，均为各函数最低，这与图形分析检验结果一致。因此，采用对变量分布上尾部变化较为敏感的Gumbel Copula函数进行联合是合理的。
若以10年一遇洪峰组合500年一遇洪量，500年一遇年最大1 d洪量远大于10年一遇洪峰设计标准，故此组合风险几乎可以忽略不计，从而反映出该组合的设计标准相对安全，但这并不符合防洪设计需求。因此，选择频率较为接近的洪峰和洪量设计标准时，其风险概率相对较小，且符合实际需求，如500年一遇的年最大洪峰和500年一遇的年最大1 d洪量的组合风险概率仅为0.01%。
混合分布模型能有效地对非一致性洪水时间序列进行模拟。
