Reservoir Risk Operation for Flood Control Under Design Flood Uncertainty
-
摘要: 水库防洪调度中，设计洪水的推求对于水库防洪安全至关重要，而在其实际的推求过程中存在着诸多不确定性因素影响。本文针对设计洪水洪峰、洪量以及洪水过程线多重不确定性，提出了一种水库防洪风险调度模型并进行求解，探求设计洪水不确定性下的水库防洪调度过程。应用Copula函数建立洪峰与最大3日洪量的联合分布模型，采用蒙特卡洛重抽样方法，获取一系列峰量联合设计值；针对传统设计洪水过程线选择的单一性问题，对水库实测洪水过程进行了分类，并对不同类型的洪水过程线进行随机模拟，生成多维不确定性下的设计洪水过程；引入经济学指标条件风险值（conditional value at risk，CVaR）来衡量水库防洪调度过程中超过调洪最高洪水位的风险，建立了考虑CVaR的水库防洪调度模型，通过不同的风险系数取值获得不同的水库防洪调度过程。以安康水库洪水过程为例，二次重现期标准下百年一遇设计洪水联合设计值存在着较大的不确定性；采用K-means聚类法把洪水过程分成了3类，并应用蒙特卡洛法对不同类型洪水过程进行随机模拟，获取了考虑不确定性的大量水库设计洪水过程，推求了安康水库不同风险系数下水库防洪调度规则。本文所建立的考虑设计洪水不确定性的水库防洪风险调度模型可为不确定性条件下水库调度规则的制定提供指导。Abstract: Design flood determines reservoir safety in the context of reservoir flood control operation. However, many uncertainties come with design flood calculation process. This paper proposes a reservoir risk operation model for flood control under design flood uncertainty and explores reservoir flood control operations. Specifically, the Copula function is employed to establish the joint distribution model of flood peak and the maximum three-day flood volume. The Monte Carlo method is used to resample the bivariate design flood under specific return period. Moreover, in order to lower the impact from selection singleness in the traditional design flood hydrograph, we classify the measured flood processes into different types and perform random simulations to generate large amount of design flood processes. Finally, the economic index conditional value at risk (CVaR) is introduced to measure the risk of the maximum reservoir level exceeding the flood control water level in the process of reservoir flood control. A reservoir risk operation model for flood control employing CVaR is established. Different reservoir flood control operations can be obtained under different risk coefficient values. The Ankang reservoir is selected as the case study. It can be seen that design flood uncertainty is considerably large with secondary type return period of one hundred years. The K-means clustering method divides the measured flood processes into three categories. Then, the Monte Carlo method is used to simulate different types of flood processes and obtain a large amount of design flood processes. Facing the uncertain design flood processes, we derive reservoir flood control operations under different risk coefficients. This paper establishes the reservoir risk operation model for flood control which can effectively copes with the uncertain design flood. The results can provide support for making reservoir flood control operations under design flood uncertainty.
-
Keywords:
- Copula function /
- Monte Carlo /
- random simulation /
- CVaR /
- flood control operation
-
设计洪水是水利工程规划设计与运行管理的重要依据，水库工程规模的确定、防洪调度规则的制定等均以确定性的设计洪水为前提条件推求得到，但设计洪水的计算受多重不确定性因素影响，如样本不确定性、典型洪水过程不确定性等，而这些不确定性因素对于水库防洪调度也存在着一定的影响，导致水库防洪调度过程中有潜在的风险。因此如何量化不确定性因素带来的风险，并据此推求在设计洪水不确定性因素下的水库防洪调度规则是确保防洪调度规则在实际应用中稳健性的必然途径。
迄今为止，国内外关于设计洪水不确定性的研究已形成了许多成果，主要是采用贝叶斯理论或蒙特卡洛方法揭示概率分布函数参数不确定性、概率分布函数类型不确定性及样本不确定性等因素对设计洪水峰值或洪量的影响，部分研究在采用Copula函数模拟设计洪水峰值和洪量相关性的基础上，开展多维设计洪水不确定性研究。杜涛等[1]在非一致性条件下设计洪水的计算中引入了水文风险的概念，并结合大气环流模型推求了特定水文风险下的设计洪水值。Reis等[2]采用马尔科夫—蒙特卡洛（markov chain monte carlo，MCMC）抽样方法对分布参数进行抽样，定量评估了参数的不确定性并证明了贝叶斯方法的有效性。Yin等[3]采用Copula函数构建洪峰流量与洪量的联合分布模型，采用参数Bootstrap抽样法评估了样本不确定性对两变量联合洪水设计值计算的影响。Guo等[4]通过引入最大熵理论，分析了频率计算中边缘分布类型不确定性与样本不确定性对洪峰与降雨2维联合设计值的影响。鲁帆等[5]以丹江口水库为例，采用基于Metropolis-Hastings抽样算法的贝叶斯MCMC方法估计设计洪水的后验分布，并据此进行频率分析。梁忠民等[6]针对水文设计值计算中概率分布函数的参数和函数类型的不确定性进行了研究，提出了基于贝叶斯理论的水文频率分析方法。胡义明等[7]利用Bootstrap方法，研究了样本不确定性对水文设计值的影响，并且对其进行了定量评价。郑永恒等[8]提出了考虑不确定性的防洪工程设计洪水风险评估方法，以混合概率分布拟合非一致性洪水序列，对混合分布重抽样，采用偏差校正与加速算法评估样本不确定性对设计洪水风险计算的影响。除却设计洪水峰值与量级的不确定性外，设计洪水过程线的不确定性（如洪峰发生时刻、洪水过程线形状等）也是影响水库防洪调度规则的重要因素，以往研究在同时考虑设计洪水峰值、量级及其过程线的不确定性方面较少涉及。
如上所述，设计洪水的不确定性必然影响水利水电工程设计规模及其运行管理[9]。李大鸣等[10]构建了洪峰和洪量的联合分布模型，并进行随机抽样，基于不同的洪水类别生成了洪水过程，采用蒙特卡洛法计算了不同风险因子组合所对应的水库调度风险，为合理利用洪水资源提供了一定的参考。尹家波[11]等建立了描述双变量设计洪水不确定性的C-PBU（copula-based parametric bootstrap uncertainty）模型，分析了联合设计值的不确定性对水库最高调洪水位的影响，并对比了不同典型洪水过程下的水位不确定性。阎晓冉等[12]提出了一种考虑峰型和频率的洪水随机模拟方法，引入了峰型系数和峰现时间，进行洪水特征量与不同类型洪水过程线的随机模拟，将其应用于防洪调度计算，对于水库防洪调度规则的制定及风险分析具有重要意义。Requena等[13]采用Copula函数构建了洪峰洪量的联合分布，生成了多场基于实测洪水过程的设计洪水，评估了大坝的漫顶风险。综上，国内外学者在评估设计洪水的不确定性对水库防洪安全的影响方面开展了大量研究，但设计洪水不确定性下，关于水库防洪调度规则应该如何制定的研究较少。
因此，本文以安康水库为例，针对设计洪水计算过程中的峰量设计值和设计洪水过程线选择的不确定性，引入CVaR指标量化设计洪水不确定性下水库防洪失事风险，建立基于CVaR的水库防洪优化调度模型并进行求解，并根据不同的风险偏好程度制定水库防洪调度决策。
1. 多维不确定性下的设计洪水计算
1.1 Copula函数
传统的单变量设计洪水将洪峰、洪量等特征参数进行单独的控制，忽略了峰量之间的随机组合的关系。Copula函数能够联立多个变量之间复杂的关系，边缘分布能够根据应用进行灵活的选择，在水文研究中有着广泛的应用。
Copula函数能够将服从指定边缘分布的若干个变量连接起来，可以用来构造洪峰洪量的联合分布。Copula函数是定义域为[0,1]均匀分布的多维联合分布函数，2维Copula函数为[14]：
$$ F\left( {x,y} \right) = {C_\theta }\left( {{F_X}\left( x \right),{F_Y}\left( y \right)} \right) = {C_\theta }\left( {u,v} \right) $$ (1) 式中，C为Copula函数，θ为Copula的参数，u=FX(x)、v=FY(y)分别为随机变量X和Y的边缘分布。它们的联合概率密度函数为：
$$ f\left( {x,y} \right) = c\left( {u,v} \right){f_X}\left( x \right){f_Y}\left( y \right) $$ (2) 式中，f (x,y)为联合概率密度函数，c (u,v)为Copula概率密度函数，fX (x)、fY (y)分别为边缘概率密度函数。
Copula函数总体来说有3类，分别是椭圆形、二次型、Archimedean型。其中，Archimedean copula函数经常运用于水文频率分析，其常用的函数形式见表1。
表 1 常用的Archimedean Copula函数Table 1 Commonly used Archimedean Copula functionsCopula函数类型 Cθ(u,v) 参数θ
范围Gumbel-Hougaard
(GH)${\rm{exp} }\left\{ { - { {\left[ { { {\left( { - {\rm{ln} }\;u} \right)}^\theta } + { {\left( { - {\rm{ln} }\;v} \right)}^\theta } } \right]}^{\frac{1}{\theta } } } } \right\}$ (1, ∞) Frank $- \dfrac{1}{\theta }{\rm{ln} }\left[ {1 + \dfrac{ {\left( { { {\rm{e} }^{ - \theta u} } - 1} \right)\left( { { {\rm{e} }^{ - \theta v} } - 1} \right)} }{ { { {\rm{e} }^{ - \theta } } - 1} } } \right]$ R Clayton $ {\left( {{u^{ - \theta }} + {v^{ - \theta }} - 1} \right)^{ - \frac{1}{\theta }}} $ (0, ∞) Copula参数
$ \theta $ 的估计采用极大似然法。似然函数的表达式为：$$\begin{aligned}[b] l\left( {\boldsymbol{\varTheta }} \right) =& \sum\limits_{i = 1}^N {\ln c\left( {{F_1}\left( {{x_1};{\theta _1}} \right),{F_2}\left( {{x_2};{\theta _2}} \right), \cdots ,{F_i}\left( {{x_i};{\theta _i}} \right);{\theta _0}} \right)} +\\& \sum\limits_{j = 1}^N {\sum\limits_{i = 1}^d {\ln \left( {{f_i}\left( {{x_i};{\theta _i}} \right)} \right)} } \\[-18pt]\end{aligned}$$ (3) 式中：N为样本容量；d为Copula函数的维数；θi为边缘分布的参数，i=1,2,···,d；θ0为联合分布的参数；Fi为第i个边缘分布的概率分布函数；fi为第i个边缘分布的概率密度函数；c为Copula函数的概率密度函数。
将似然函数关于参数
${{{{\boldsymbol{\varTheta}}}}}$ 最大化，得到参数向量${{{{\boldsymbol{\varTheta}}}}}$ 的评估值为[15]：$$ \hat {{{\boldsymbol{\varTheta}}}} = \arg \max l\left( {{{\boldsymbol{\varTheta}}}} \right) $$ (4) 1.2 二次重现期
在水文计算中，通常采用“重现期”来表示某一特定大小洪水出现的平均时间间隔，它是衡量洪水量级大小的一种方法。在计算二次重现期时，引入了kendall测度Kc，Kc只与联合分布函数Copula有关，表示在给定概率水平t∈(0,1)下，特征量联合分布概率W小于或等于t的概率[16]，表达式为：
$$ {K_{\rm{c}}}\left( t \right) = P\left( {W \leqslant t} \right) = P\left( {C\left( {u,v} \right) \leqslant t} \right) $$ (5) 式中，C(u,v)为Copula函数，Kc为kendall测度。
则以Kc表示的kendall重现期Tk为：
$$ {T_k} = \frac{1}{{P\left( {C\left( {u,v} \right) \geqslant t} \right)}} = \frac{1}{{1 - {K_{\rm{c}}}\left( t \right)}} $$ (6) 采用最可能组合法计算洪水联合设计值。最可能联合设计值是指满足指定的防洪标准T下，洪峰洪量联合概率密度最大对应的联合设计值[17]，计算公式为：
$$ \left( {{x_m},{y_m}} \right) = \mathop {\arg \max }\limits_{\left( {x,y \in S\hat p} \right)} f\left( {x,y} \right) $$ (7) $$ f\left( {x,y} \right) = c\left( {u,v} \right)f\left( x \right)f\left( y \right) $$ (8) 式中，(xm，ym)为推求的洪峰洪量最可能设计值，f (x)、f (y)分别为边缘概率密度函数，
$ S\hat p $ 为具有相同重现期的峰量组合构成的等值线，c(u，v)为Copula的概率密度函数，f (x，y)为洪峰、洪量联合概率密度函数。2. 设计洪水的随机模拟
蒙特卡洛随机抽样以概率统计理论为基础，是一种对随机变量进行数理统计实验及分布概率模拟，从而近似求解得到预测值的方法[18]。因此本文采用蒙特卡洛随机抽样并结合高超等[19]论文中提及的不同类型洪水过程线的随机模拟方法，进行多维不确定性下设计洪水的计算。步骤如下：
1）用最大似然法估计联合分布函数H(x,y;
${\boldsymbol{\varTheta }}$ )=C(F(x;θx),F(y;θy);θ0)的参数向量${{{{\boldsymbol{\varTheta}}}} }$ =($ \theta $ x,$ \theta $ y;$ \theta $ 0)，其中，变量x与y为年最大洪水峰值与量级。2）基于上述联合分布函数，应用蒙特卡洛方法，从中抽取与实测数据长度相等的2维样本，随机抽取N次。
3）采用最大似然法估计每组模拟样本的参数，得到N组联合分布函数参数向量。
4）在每一种参数向量下，计算不同重现期类型下的2维联合设计值。
5）为了移除洪峰、洪量等特征量的影响，将流域内收集的洪水过程线进行无量纲化处理。
$$ \tau = \frac{t}{T},{F_\tau } = \frac{{{W_t}}}{{{W_T}}} $$ (9) 式中：T为洪水历时；τ为时刻t的无量纲时间，τ∈(0, 1]；WT为一场洪水的总洪量；Wt为时刻t的累积洪量；Fτ为无量纲累积洪量，即洪量随时间的累积百分比，Fτ∈(0, 1]。
6）将无量纲后的洪量累积曲线均分为K个时段，此时，无量纲时间τi（τi = i/k，i=1, 2, ···, K）对应的无量纲累积洪量为Fi，则每个时段的无量纲洪量为Pi = Fi–Fi-1。将每场洪水过程每个时段的Fi或Pi作为输入值，采用基于欧式距离的K-means聚类法[20]对其进行分类，这种方法不仅简单、容易操作，并且分类效果较好。
7）对于分类后的洪水过程，由于其无量纲化后的洪水过程线通常是相关的非正态变量，直接模拟比较困难。因此，先将其进行对数转换将受约束的相关非正态多变量转换为不受约束的相关非正态多变量，采用Johnson系统函数进行正态转换将其转换为相关标准正态多变量，再通过Cholesky分解进行正交转换将其转换为独立的标准正态多变量，最后采用蒙特卡洛模拟生成独立的标准正态多变量，进行逆转换。
8）将生成的无量纲洪水过程线与模拟得到的洪峰、洪量采用变倍比放大的方式融合，从而生成完整的洪水过程。
经过以上步骤，则可以获得若干条设计洪水序列，且与实测洪水过程接近。
3. 水库防洪风险调度
3.1 CVaR理论
CVaR是Rockafellar和Uryasev提出来的一种风险度量技术，常用于刻画尾部风险，是在风险值（value at risk，VaR）的基础上得到的，表示的是超过设定VaR时可能遭受的平均潜在损失[21]，常被用于金融和市场方面的风险分析。CVaR克服了VaR不满足凸性与一致性风险测度所要求的次可加性，并且缺乏对尾部风险的控制这些缺陷，是一种更为有效的风险管理方法。CVaR理论同样可以推广到水库调度的风险分析中，用来描述水库损失风险。
在给定的置信度
$ \;\beta \in \left(\mathrm{0,1}\right) $ 下，由定义可得条件风险值CVaR的公式为：$$ \begin{aligned}[b] {\rm{CVa}}{{\rm{R}}_\beta }\left( x \right) =& {\rm{Va}}{{\rm{R}}_\beta }\left( x \right) + E[ f\left( {x,y} \right) - {\rm{Va}}{{\rm{R}}_\beta }\left( x \right)| f\left( {x,y} \right) \geqslant \\& {\rm{Va}}{{\rm{R}}_\beta }\left( x \right) ] = E[ {f\left( {x,y} \right)| {f\left( {x,y} \right) \geqslant {\rm{Va}}{{\rm{R}}_\beta }\left( x \right)}} ] \end{aligned}$$ (10) 采用概率密度函数的表达形式为：
$$ {\rm{CVa}}{{\rm{R}}_\beta }\left( x \right) = \frac{1}{{1 - \beta }}\int_{f\left( {x,y} \right) \geqslant {\rm{Va}}{{\rm{R}}_\beta }\left( x \right)} {f\left( {x,y} \right)} p\left( y \right){\rm{d}}y $$ (11) 广义的风险是指某一事件在特定的条件下所发生的不利事件及其概率并由此产生的损失程度。水库防洪调度过程中的风险是指由于各种不确定因素的影响，在水库防洪调度的过程中发生库水位或者下泄量高于某一规定值等非期望事件的可能性[22]。水库防洪风险是一个复杂的系统，会受到各种不确定性因素的影响，比如入库洪水过程、调度滞时、预报误差、水位泄流能力等，各种因素互相组合作用于防洪目标，使其产生风险，本文仅考虑在设计洪水不确定的条件下所造成的风险。
在防洪调度过程中首先要考虑的问题是确保大坝自身的安全，即遭遇洪水时，水库最高库水位不得超过校核洪水位。常规的水库防洪调度中风险事件通常指水库水位超过校核或设计洪水位的可能性。因此在不确定因素的作用下，可以选择调度过程中最高调洪水位的CVaR作为刻画水库失事风险的评价指标。
则在给定风险厌恶程度β∈(0,1)下,最高调洪水位的条件风险价值定义为：
$$ {\rm{Va}}{{\rm{R}}_\beta }\left( {{Z_{{\rm{max}}}}} \right) = E\left[ {{Z_{{\rm{max}}}}\left| {{Z_{{\rm{max}}}} \geqslant {\rm{VaR}}\left( {{Z_{{\rm{max}}}}} \right)} \right.} \right] $$ (12) $ \;\beta $ 越小，决策者风险厌恶程度越高。3.2 基于CVaR的防洪优化调度模型
水库防洪调度在多重不确定因素的作用下，不仅要追求整体的期望风险最小，还要尽可能地规避潜在的风险。风险度量在不确定性条件的优化问题中起着至关重要的作用[23]。因此，本文在常规水库防洪优化调度目标的基础上，额外增加了风险衡量方法CVaR进行风险度量，即同时考虑水库调洪过程中由于水位过高所造成的水库期望损失以及最高损失风险。
目标函数：
$${\qquad {\rm{min}}G\left( x \right) = \left( {1 - \lambda } \right)E\left( x \right) + \lambda {\rm{CVaR}}\left( x \right) }$$ (13) 约束条件：
1）水量平衡约束：
$$ {V_{t + 1}} = {V_t} + \left( {{Q_t} - {q_t}} \right) \times \Delta t $$ (14) 2）水库库容约束：
$$ {V_{\min }} \leqslant {V_t} \leqslant {V_{\max }} $$ (15) 3）水库水位约束：
$$ {Z_{\min }} \leqslant {Z_t} \leqslant {Z_{\max }} $$ (16) 4）泄流能力约束：
$$ {q_t} \leqslant f\left( {{q_t}} \right) $$ (17) 5）泄量变幅约束：
$$ \left| {{q_t} - {q_{t - 1}}} \right| \leqslant \Delta q $$ (18) 式（13）～（18）中：x为每场洪水的最高调洪水位；G(x)为目标函数；E(x)为最高调洪水位的期望；CVaR(x)为CVaR度量准则下的最高调洪水位；λ定义为风险系数，取值范围为λ∈(0,1)，用来权衡防洪调度最高水位的均值和水位变动的风险之间的关系；Vt、Zt为当前时段水库水位和蓄水量；Vt+1为下一时段水库蓄水量；Qt、qt分别为当前时段平均入库流量和出库流量；Vmin、Vmax和Zmin、Zmax分别为水库蓄水量和水位的最大值与最小值；f(qt)为t时刻水库的泄流能力；∆q水库允许的最大下泄流量变化量，m3/s。
当λ取值为0时，此时目标函数为最高调洪水位的期望，调度人员为风险中立者，即调度人员不关心来水不确定性因素造成的水位波动风险，既不规避风险，也不主动追求风险，只考虑总体期望最小。λ的加入使得在调洪过程中进一步得考虑到了潜在的风险，当
$ \lambda $ 取值较小时，调度人员开始尝试规避风险，此时CVaR值较大，运行人员追求较小的防洪库容但面临较大的波动风险。随着$ \lambda $ 的取值增加，调洪潜在风险所占的比重在目标函数中也随之增加，当$ \lambda $ 取值较大时，调度人员极力规避风险，此时CVaR值较小，运行人员追求较低的波动风险但面临使用较高的防洪库容规模。故而，防洪库容的运用和$ \lambda $ 的取值有关，调度人员可以根据风险偏好程度制定合理的防洪调度决策，当调度人员追求使用较小的防洪库容而不考虑来水不确定性所造成的水库库容使用的波动风险时，可选择偏小的$ \lambda $ 值；当调度人员考虑到由于入库洪水不确定性所造成的防洪库容使用情况时，可选择较大的$ \lambda $ 值进行风险规避，此时防洪库容占用规模会较大。4. 实例应用
4.1 研究区域概况
本研究以汉江上游流域安康水库为研究对象，是一座具有发电、航运、防洪、养殖、旅游等综合效益的大型水电枢纽工程，地理位置如图1所示。安康水库坝址以上流域面积35 700 km2，多年平均流量608 m3/s。水库正常蓄水位330 m，相应库容25.85×108 m3，死水位300 m，调节库容16.7×108 m3。汉江上游流域洪水成因多为暴雨，呈山溪性陡涨陡落特点，故夏季入库洪水历时短，峰型尖瘦，秋季洪水历时长，峰型胖。历史上发生过几次大洪水，如：1983年7月上游流域发生了200年一遇的大洪水，安康老城遭受淹没之灾，造成重大损失。同时，汉江上游作为国家南水北调中线工程和陕西省引汉济渭工程重要的水源地，不确定洪水事件的发生将直接影响到工程的运行和效益以及安康等城市的防洪安全。
图 1 汉江上游流域图Fig. 1 Map of the upper reaches of the Han River basin
下载:
全尺寸图片
4.2 双变量设计洪水的计算
本文以安康水库1956年到2013年共58年的最大洪峰和最大3日洪量资料作为研究材料。从P3（Pearson type 3 distribution）分布、Gamma分布、Logn（Lognormal distribution）分布、Gp（Generalized Pareto distribution）分布、Gev（Generalized Extreme Value distribution）分布、Norm（Normal distribution）分布中分别对于洪峰、最大3日洪量数据选择拟合程度较高的边缘分布函数，采用均方根误差准则RMSE（Root Mean Squared Error）、AICc准则（corrected Akaike Information Criterion）和BIC准则（Bayesian Information Criterion）进行拟合优度检验。
拟合优度检验及K-S（Kolmogorov-Smirnov）检验值见表2。可以看出，Gp分布未通过显著性检验，所以不能作为洪峰的备选边缘分布类型。选择Gamma分布作为洪峰与最大3日洪量的边缘分布类型。图2为洪峰和最大3日洪量的边际分布曲线。
表 2 边缘分布拟合优度检验Table 2 Fit test results of marginal distribution指标 边缘分布
函数类型K-S检验 RMSE AICc BIC 洪峰 Gev 0 0.0323 1163.6489 1168.9414 Logn 0 0.0399 1162.1070 1165.7915 Gp 1 0.1033 1186.6071 1190.2916 Gamma 0 0.0305 1160.0285 1163.7130 Norm 0 0.0309 1166.6515 1170.3361 P3 0 0.0305 1162.4810 1167.7735 最大3
日洪量Gev 0 0.0408 954.1038 959.3963 Logn 0 0.0430 954.1739 957.8584 Gp 0 0.0574 956.8394 960.5239 Gamma 0 0.0373 950.4106 954.0951 Norm 0 0.0429 954.6240 958.3085 P3 0 0.0373 952.8400 958.1325 
图 2 边际分布曲线Fig. 2 Marginal distribution curves下载:
全尺寸图片
边缘分布函数类型选定后，下一步就是从Clayton、Frank、GH 3种类型中选择合适的联结函数，采用极大似然法估计联合分布的参数。之后，采用AICc准则和BIC准则（Bayesian 信息准则）选择拟合程度较高的Copula函数类型，计算结果见表3。从中选择AICc和BIC值较小的分布函数组合，因此，选择GH copula函数作为洪峰和最大3日洪量的联合分布函数。
表 3 联合分布拟合优度检验Table 3 Fit test results of joint distribution fit accuracy testCopula函数类型 AICc BIC Gumbel-Gougaard(GH) –84.89 –84.92 Frank –82.63 –82.66 Clayton –67.16 –67.20 根据得到的洪水变量的分布函数类型，结合二次重现期的计算公式，采用最可能组合法计算得到二次重现期标准下的联合设计值，计算结果见表4。采用二次重现期作为计算标准所推求的联合设计值略小于单变量设计值，与安康水库原设计值进行对比，洪峰流量在高重现期时较小，低重现期时较大，3日洪量设计值均大于原设计值。其考虑了峰量之间的相关性，能够更加合理地描述多变量洪水事件，为防洪工程设计提供了新的理论基础。
表 4 不同重现期下设计值对比Table 4 Design flood values under different return periods重现期/a 安康水库原设计值 单变量设计值 二次重现期 洪峰/
(m3·s–1)最大3
日洪量/
(108 m3)洪峰/
(m3·s–1)最大3
日洪量/
(108 m3)洪峰/
(m3·s–1)最大3
日洪量/
(108 m3)10000 45000 70.39 43936 73.50 43443 72.63 1000 36700 57.66 36427 60.21 35910 59.30 100 28100 44.16 28430 46.17 27874 45.19 20 21500 33.88 22303 35.52 21714 34.50 5 15100 23.97 16246 25.14 15665 24.14 联合设计值是否准确关系到流域水利工程是否经济合理，运行是否安全可靠。采用设计标准过高，会造成工程规模过大，导致投资增加；如果采用设计标准过低，则会造成水利工程的规模较小，不利于水利工程发挥自身功能，严重的话可能导致工程事故，以及造成巨大的生命财产的损失[7]。由于安康水库的防洪标准为百年一遇，因此本文以重现期为百年一遇的设计洪水为例，采用蒙特卡洛方法，对原始样本进行多次重抽样，获得二次重现期下的峰量联合设计值，结果见图3。抽样值分散在采用最可能组合法计算的联合设计值周围，可以看出，重现期为100年时，联合设计值约跨越了100至500 a的重现期，存在较大的不确定性。
图 3 重现期为100 a一遇的抽样值Fig. 3 Uncertain design flood with 100-year return period下载:
全尺寸图片
4.3 洪水过程的随机模拟
选取安康水库1954年至2009年36场完整的入库洪水过程，为了移除洪水特征量对洪水形状的影响，首先对其进行无量纲化处理，将洪水过程线划分为12段。之后采用K-means聚类法对其进行分类，将安康站的洪水过程类型分为3类，结果如图4所示。其中Ⅰ类和Ⅱ类洪水均表示洪峰出现时间靠前的类型，不同的是Ⅰ类洪水的峰现时间比Ⅱ类更靠前，Ⅲ类洪水表示的是洪峰出现时间中间偏后的类型。通过对不同类型洪水的分析，可以看出，安康站的洪峰基本没有出现在比较靠后的时间段，符合该地区的洪水特征。
图 4 洪水过程线的分类Fig. 4 Classification of flood process下载:
全尺寸图片
采用上述所描述的方法对不同类型的洪水过程进行随机模拟。首先，采用对数转换将无量纲洪量累积曲线转换为不受约束的相关非正态多变量，并采用Johnson分布系统进行正态转换。其次，采用蒙特卡洛随机模拟正态随机变量进行逆转换。最终，求得不同类型的洪水过程。为了评价模拟精度的好坏，将模拟的无量纲累积洪量的均值与实测洪水的无量纲累积洪量的均值进行比较，结果见表5。以Ⅱ类洪水过程为例进行模拟，模拟结果如图5所示，从图5可以看出，洪水模拟的结果较好，能够运用到实际中。
表 5 洪水过程线模拟结果评价Table 5 Evaluation of flood process line simulation results洪水类型 RMSE Ⅰ类 0.0256 Ⅱ类 0.0582 Ⅲ类 0.0553 
图 5 Ⅱ类洪水过程模拟结果Fig. 5 Simulation results of Category Ⅱ flood process下载:
全尺寸图片
将模拟的无量纲洪水过程线与洪水特征量采用变倍比的方法融合，从而形成完整的洪水过程。该方法既能满足洪峰洪量的要求，又能保持洪水过程线的形状不变。分别选取不同类别的一场洪水进行融合，结果如图6所示。其中：图6（a）洪峰流量为28 409.90 m3，洪量为44.04×108 m3，峰现时间靠前，洪水陡涨陡落，属于Ⅰ类洪水；图6（b）洪峰流量为28 588.29 m3，洪量为47.37×108 m3，洪峰出现时间为中间靠前，属于Ⅱ类洪水；图6（c）洪峰流量为28 488.39 m3，洪量为43.11×108 m3，峰现时间为中间偏后，属于Ⅲ类洪水过程。经过上述步骤，任意场不同类型的洪水过程均可随机模拟，克服了传统设计洪水过程中典型洪水过程线选取的唯一性，可以为防洪规划设计提供依据。
图 6 不同类别的洪水模拟过程Fig. 6 Different types of flood simulation process下载:
全尺寸图片
4.4 水库防洪风险调度
为了验证基于CVaR准则的防洪优化调度模型的有效性，选取重现期为100 a一遇的每种类型的100场洪水作为输入。分析不同风险系数对水库调度决策及相应的均值–CVaR目标的影响，令参数λ={0，0.2，0.4，0.6，0，8，1.0}，风险厌恶程度
$ \;\beta $ 取固定值0.95，采用粒子群算法求解模型，从而统计得到调洪高水位的期望、CVaR与相应的目标值，结果见表6。表 6 均值–CVaR模型与常规调度模型下水库水位结果Table 6 Reservoir level results obtained by mean-CVaR model and general scheduling modelm 风险系数 (1-λ)E+λCVaR 均值 CVaR 0 328.23 328.23 333.87 0.2 329.21 328.30 332.86 0.4 330.24 328.52 332.82 0.6 331.15 328.77 332.74 0.8 331.80 328.87 332.53 1.0 332.28 329.17 332.28 常规调度 — 328.56 332.60 图7为安康水库坝前水位为均值、CVaR与风险系数的关系。
图 7 安康水库坝前水位均值、CVaR与风险系数的关系Fig. 7 Relationship between the mean, CVaR of Ankang reservoir level and risk factor下载:
全尺寸图片
由图7和表6的结果可知，随着风险系数
$ \lambda $ 的增加，最高调洪水位的CVaR随之减少，期望值随之增大，说明条件风险值占的比例越大，对均值–CVaR的影响就越高，而最高调洪水位的条件风险值总是大于其期望风险值。当λ=0时，等价于期望风险模型，此时，CVaR为328.23 m，最高调洪水位期望值E(Z)为333.87 m，决策者只考虑期望风险而忽视条件风险值，导致此时CVaR值达到最大；当λ=1时，等价于条件风险模型，此时，CVaR为329.17 m，E(Z)为332.28 m，此时期望风险达到最大，但CVaR值达到最小。在不确定性随机优化过程中，$ \lambda $ 的取值决定了调度过程中调度策略的选取，$ \lambda $ 取值越大，调度过程中会倾向于追求规避风险的调度策略。相反，$ \lambda $ 取值越小，调度目标越偏向于期望风险最小，然而却很难抵御不确定性洪水过程对调度目标的影响。当遭遇100 a一遇洪水时，安康水库坝前最高洪水位不得超过防洪高水位330 m，这是安康水库洪水调度的主要原则之一，而不论$ \lambda $ 取值如何，其最高调洪水位的期望总是小于330 m，说明在不确定性设计洪水的影响下，安康水库仍具有一定抵御100 a一遇洪水的能力。与水库原有防洪调度规则调洪结果相比，当风险系数减小时，水位的期望值也随之减少，但CVaR值却随之增大，并且当其取0.4以下时，期望值均小于常规调度的期望值，在风险系数等于0.8时，水位的CVaR值小于常规调度的结果。从图7可以看出，伴随着风险系数的变化，均值–CVaR模型能够有效降低水库调度过程中由于水位过高所造成的风险，一是水库和大坝本身的安全风险，二是迎接后续洪水的防洪调度风险。 并且风险系数对水库的调洪过程有一定的影响，不同的风险系数对应不同的风险喜好程度，同时也会对应不同的水库防洪调度过程。因此，水库在实际中面对来临洪水的不确定时，应同时考虑调洪过程中水位的期望风险与条件风险值，才能减少水库调度过程中因为来水不确定性所造成的水库损失。风险系数的选择决定了调度人员风险态度的考虑，可进一步为水库操作带来的损失风险提供有效的决策支持。
典型洪水过程选择的原则为对水库防洪不利，主峰靠后、峰高量大，Ⅲ类洪水为峰现时间偏后的洪水过程，因此以Ⅲ类洪水过程为例，选择100场100 a一遇设计洪水过程带入调度模型，统计得到不同时刻水库水位的变化过程，如图8所示。从图8可以看出，在不同的风险系数下，水库水位在涨洪初期的变化并不大，而在洪峰出现时期水库水位有着明显的波动，并且维持较高水位的时间较长，在调度期未很难回到汛限水位，存在着较大的风险。随着不同风险系数的增加，水库最高调洪水位也随之变化，洪峰附近调洪水位的异常值随之减少，说明水库调洪风险也随之减小。在调洪过程中，若是考虑期望风险较多，则风险系数可取小值，在考虑不确定性洪水过程的调度中可降低整体风险；反之，若是注重考虑条件风险，则风险系数取大值，可降低不确定性洪水对应的调洪过程所造成的极端水位，从而保障水库自身安全，更好地应对后续可能发生的洪水过程。
图 8 不同风险系数下水库水位变化结果Fig. 8 Reservoir level under different risk factors下载:
全尺寸图片
从水库规划的角度来看，设计洪水的不确定性对于防洪安全有着极大的影响。综上考虑，在如此高的不确定性下，水库若需维持原设计标准，可以通过提前降低汛限水位以增加防洪库容，或者通过加强洪水预报的精度使水库进行提前下泄而腾出库容，方能保证防洪安全。
5. 结 论
本文以汉江上游安康水库为例，在基于设计洪水的特征量与洪水过程线不确定性的条件下，引入了条件风险值CVaR表示防洪调度中所产生的风险，建立了考虑风险厌恶程度的水库防洪调度模型，并求解了不同风险厌恶系数下水库相应的调度规则，得到以下结论：
1）采用Copula函数建立了边缘分布均为Gamma分布的峰量联合分布模型，采用最可能组合法计算二次重现期水平下的洪水设计值，并采用蒙特卡洛方法对指定二次重现期标准下的设计洪水进行了重抽样，结果表明100 a一遇的二次重现期标准下的联合设计值存在着较大的不确定性，对水利工程的建设提出了巨大挑战。
2）对流域内的洪水过程进行了分类，针对传统设计洪水过程线选择的单一性，采用蒙特卡洛法进行不同类型洪水过程的随机模拟，该方法模拟的洪水过程更接近实际，对防洪规划具有十分重要的意义。
3）针对传统的水库防洪调度中风险的定义，引入了经济学指标条件风险值CVaR的概念来衡量风险，建立了水库防洪风险调度模型，得到了不同风险系数下水库相应的调度过程，能够有效降低不确定条件下水库调度过程中大坝失事的风险。并以其中一种类型的100 a一遇设计洪水过程为例，统计得到其调度过程，可以看出在洪峰附近水库水位的变化较大，水库可采取提前预泄的方式保证防洪安全。
-
图 1 汉江上游流域图
Fig. 1 Map of the upper reaches of the Han River basin
下载:
全尺寸图片
图 2 边际分布曲线
Fig. 2 Marginal distribution curves
下载:
全尺寸图片
图 3 重现期为100 a一遇的抽样值
Fig. 3 Uncertain design flood with 100-year return period
下载:
全尺寸图片
图 4 洪水过程线的分类
Fig. 4 Classification of flood process
下载:
全尺寸图片
图 5 Ⅱ类洪水过程模拟结果
Fig. 5 Simulation results of Category Ⅱ flood process
下载:
全尺寸图片
图 6 不同类别的洪水模拟过程
Fig. 6 Different types of flood simulation process
下载:
全尺寸图片
图 7 安康水库坝前水位均值、CVaR与风险系数的关系
Fig. 7 Relationship between the mean, CVaR of Ankang reservoir level and risk factor
下载:
全尺寸图片
图 8 不同风险系数下水库水位变化结果
Fig. 8 Reservoir level under different risk factors
下载:
全尺寸图片
表 1 常用的Archimedean Copula函数
Table 1 Commonly used Archimedean Copula functions
Copula函数类型 Cθ(u,v) 参数θ
范围Gumbel-Hougaard
(GH)${\rm{exp} }\left\{ { - { {\left[ { { {\left( { - {\rm{ln} }\;u} \right)}^\theta } + { {\left( { - {\rm{ln} }\;v} \right)}^\theta } } \right]}^{\frac{1}{\theta } } } } \right\}$ (1, ∞) Frank $- \dfrac{1}{\theta }{\rm{ln} }\left[ {1 + \dfrac{ {\left( { { {\rm{e} }^{ - \theta u} } - 1} \right)\left( { { {\rm{e} }^{ - \theta v} } - 1} \right)} }{ { { {\rm{e} }^{ - \theta } } - 1} } } \right]$ R Clayton $ {\left( {{u^{ - \theta }} + {v^{ - \theta }} - 1} \right)^{ - \frac{1}{\theta }}} $ (0, ∞) 表 2 边缘分布拟合优度检验
Table 2 Fit test results of marginal distribution
指标 边缘分布
函数类型K-S检验 RMSE AICc BIC 洪峰 Gev 0 0.0323 1163.6489 1168.9414 Logn 0 0.0399 1162.1070 1165.7915 Gp 1 0.1033 1186.6071 1190.2916 Gamma 0 0.0305 1160.0285 1163.7130 Norm 0 0.0309 1166.6515 1170.3361 P3 0 0.0305 1162.4810 1167.7735 最大3
日洪量Gev 0 0.0408 954.1038 959.3963 Logn 0 0.0430 954.1739 957.8584 Gp 0 0.0574 956.8394 960.5239 Gamma 0 0.0373 950.4106 954.0951 Norm 0 0.0429 954.6240 958.3085 P3 0 0.0373 952.8400 958.1325 表 3 联合分布拟合优度检验
Table 3 Fit test results of joint distribution fit accuracy test
Copula函数类型 AICc BIC Gumbel-Gougaard(GH) –84.89 –84.92 Frank –82.63 –82.66 Clayton –67.16 –67.20 表 4 不同重现期下设计值对比
Table 4 Design flood values under different return periods
重现期/a 安康水库原设计值 单变量设计值 二次重现期 洪峰/
(m3·s–1)最大3
日洪量/
(108 m3)洪峰/
(m3·s–1)最大3
日洪量/
(108 m3)洪峰/
(m3·s–1)最大3
日洪量/
(108 m3)10000 45000 70.39 43936 73.50 43443 72.63 1000 36700 57.66 36427 60.21 35910 59.30 100 28100 44.16 28430 46.17 27874 45.19 20 21500 33.88 22303 35.52 21714 34.50 5 15100 23.97 16246 25.14 15665 24.14 表 5 洪水过程线模拟结果评价
Table 5 Evaluation of flood process line simulation results
洪水类型 RMSE Ⅰ类 0.0256 Ⅱ类 0.0582 Ⅲ类 0.0553 表 6 均值–CVaR模型与常规调度模型下水库水位结果
Table 6 Reservoir level results obtained by mean-CVaR model and general scheduling model
m 风险系数 (1-λ)E+λCVaR 均值 CVaR 0 328.23 328.23 333.87 0.2 329.21 328.30 332.86 0.4 330.24 328.52 332.82 0.6 331.15 328.77 332.74 0.8 331.80 328.87 332.53 1.0 332.28 329.17 332.28 常规调度 — 328.56 332.60 -
[1] Du Tao,Xiong Lihua,Li Shuai,et al.Risk-based nonstationary design flood and uncertainty analysis[J].Journal of Hydraulic Engineering,2018,49(2):241–253. doi: 10.13243/j.cnki.slxb.20160862 [2] Reis Jr D S,Stedinger J R.Bayesian MCMC flood frequency analysis with historical information[J].Journal of Hydrology,2005,313(1/2):97–116. [3] Yin J,Guo S,Liu Z,et al.Uncertainty analysis of bivariate design flood estimation and its impacts on reservoir routing[J].Water Resources Management,2018,32(5):1795–1809. doi: 10.1007/s11269-018-1904-x [4] Guo A,Chang J,Wang Y,et al.Maximum entropy-copula method for hydrological risk analysis under uncertainty:A case study on the loess plateau,China[J].Entropy,2017,19(11):609. doi: 10.3390/e19110609 [5] 鲁帆,严登华.基于广义极值分布和Metropolis-Hastings抽样算法的贝叶斯 MCMC洪水频率分析方法[J].水利学报,2013,44(8):942–949. doi: 10.13243/j.cnki.slxb.2013.08.004 Lu Fan,Yan Denghua.Bayesian MCMC flood frequency analysis based on generalized extreme value distribution and Metropolis-Hastings algorithm[J].Journal of Hydraulic Engineering,2013,44(8):942–949 doi: 10.13243/j.cnki.slxb.2013.08.004 [6] 梁忠民,李磊,王军,等.考虑参数和线型不确定性的水文设计值估计的贝叶斯方法[J].天津大学学报,2010,43(5):379–384. doi: 10.3969/j.issn.0493-2137.2010.05.001 Liang Zhongmin,Li Le,Wang Jun,et al.Bayesian method for hydrological frequency analysis considering uncertainties of parameter and model[J].Journal of Tianjin University,2010,43(5):379–384 doi: 10.3969/j.issn.0493-2137.2010.05.001 [7] 胡义明,梁忠民,王军,等.考虑抽样不确定性的水文设计值估计[J].水科学进展,2013,24(5):667–674. doi: 10.14042/j.cnki.32.1309.2013.05.018 Hu Yiming,Liang Zhongmin,Wang Jun,et al.Determination of design hydrologic characteristics withsampling uncertainty considerations[J].Advances in Water Science,2013,24(5):667–674 doi: 10.14042/j.cnki.32.1309.2013.05.018 [8] 郑永恒,郭爱军.考虑不确定性的变化环境下防洪工程设计洪水风险评估[J].电网与清洁能源,2018,34(5):69–73. Zheng Yongheng,Guo Aijun.Uncertainty evaluation of risk analysis of design flood for flood control projects under changing environment[J].Power System and Clean Energy,2018,34(5):69–73 [9] 康有,张军良.基于混合抽样的洪水设计值抽样误差计算方法[J].水电能源科学,2019,37(11):87–91. Kang You,Zhang Junliang.Sampling error calculation of design flood based on the hybrid sampling[J].Water Resources and Power,2019,37(11):87–91 [10] 李大鸣,顾利军,高正廉,等.基于Copula函数的桃林口水库防洪风险分析[J].水利水电技术,2017,48(3):158–164. doi: 10.13928/j.cnki.wrahe.2017.03.028 Li Daming,Gu Lijun,Gao Zhenglian,et al.Copula function-based risk analysis for flood control of Taolinkou Reservoir[J].Water Resources and Hydropower Engineering,2017,48(3):158–164 doi: 10.13928/j.cnki.wrahe.2017.03.028 [11] 尹家波,郭生练,吴旭树,等.两变量设计洪水估计的不确定性及其对水库防洪安全的影响[J].水利学报,2018,49(6):715–724. doi: 10.13243/j.cnki.slxb.20180038 Yin Jiabo,Guo Shenglian,Wu Xushu,et al.Uncertainty of bivariate design flood estimation and its impact on reservoir flood prevention[J].Journal of Hydraulic Engineering,2018,49(6):715–724 doi: 10.13243/j.cnki.slxb.20180038 [12] 阎晓冉,王丽萍,张验科,等.考虑峰型及其频率的洪水随机模拟方法研究[J].水力发电学报,2019,38(12):61–72. doi: 10.11660/slfdxb.20191207 Yan Xiaoran,Wang Liping,Zhang Yanke,et al.Study on stochastic flood simulation method considering peak shape and its frequency[J].Journal of Hydroelectric Engineering,2019,38(12):61–72 doi: 10.11660/slfdxb.20191207 [13] Requena A I,Mediero L,Garrote L.A bivariate return period based on copulas for hydrologic dam design:Accounting for reservoir routing in risk estimation[J].Hydrology and Earth System Sciences,2013,17(8):3023–3038. doi: 10.5194/hess-17-3023-2013 [14] 刘和昌,梁忠民,姚轶,等.基于Copula函数的水文变量条件组合分析[J].水力发电,2014,40(5):13–16. doi: 10.3969/j.issn.0559-9342.2014.05.004 Liu Hechang,Liang Zhongmin,Yao Yi,et al.Analysis on conditional compositions of hydrological variants based on Copula function[J].Water Power,2014,40(5):13–16 doi: 10.3969/j.issn.0559-9342.2014.05.004 [15] 李天元.基于Copula函数的设计洪水计算方法研究[D].武汉:武汉大学,2014. Li Tianyuan.Design flood estimation based on Copulas[D].Wuhan:Wuhan University,2014. [16] 范嘉炜,黄锦林.基于Kendall重现期的降雨潮位风险分析[J].水电能源科学,2017,35(5):21–24. Fan Jiayi,Huang Jinlin.Risk analysis of rainstorm and tidal level based on kendall return period[J].Water Resources and Power,2017,35(5):21–24 [17] 刘章君,郭生练,许新发,等.两变量洪水结构荷载重现期与联合设计值研究[J].水利学报,2018,49(8):956–965. doi: 10.13243/j.cnki.slxb.20180466 Liu Zhangjun,Guo Shenglian,Xu Xinfa,et al.Bivariate structure load return period and joint flood quantile estimation[J].Journal of Hydraulic Engineering,2018,49(8):956–965 doi: 10.13243/j.cnki.slxb.20180466 [18] García–Alonso C R,Arenas–Arroyo E,Pérez–Alcalá G M.A macro-economic model to forecast remittances based on Monte-Carlo simulation and artificial intelligence[J].Expert Systems with Applications,2012,39(9):7929–7937. doi: 10.1016/j.eswa.2012.01.108 [19] 高超,朱聪,泮苏莉,等.不同类型洪水过程线的随机模拟[J].应用基础与工程科学学报,2018,26(4):767–779. doi: 10.16058/j.issn.1005-0930.2018.04.008 Gao Chao,Zhu Cong,Pan Suli,et al.Stochastic simulation of different flood hydrographs[J].Journal of Basic Science and Engineering,2018,26(4):767–779 doi: 10.16058/j.issn.1005-0930.2018.04.008 [20] Zahmatkesh Z,Karamouz M,Nazif S.Uncertainty based modeling of rainfall-runoff:Combined differential evolution adaptive metropolis (DREAM) and K-means clustering[J].Advances in Water Resources,2015,83(Sep):405–420. doi: 10.1016/j.advwatres.2015.06.012 [21] Rockafellar R T,Uryasev S.Optimization of conditional value-at-risk[J].Journal of Risk,2000,2(3):21–42. doi: 10.21314/JOR.2000.038 [22] 杨礼祥,何子干.水库防洪调度风险评估[J].人民长江,2007,38(8):69–71. doi: 10.16232/j.cnki.1001-4179.2007.08.021 Yang Lixiang,He Zigan.Risk assessment of reservoir flood control operation Yangtze River[J].Yangtze River,2007,38(8):69–71 doi: 10.16232/j.cnki.1001-4179.2007.08.021 [23] 张萌萌,董军.基于CVaR的灵活综合能源系统随机调度优化模型[J].全球能源互联网,2020,3(3):301–309. Zhang Mengmeng,Dong Jun.Stochastic scheduling model for flexible integrated energy system based on CVaR[J].Journal of Global Energy Interconnection,2020,3(3):301–309
