工程科学与技术   2017, Vol. 49 Issue (2): 69-76
设计洪水峰量最可能组合法的计算通式
尹家波, 郭生练, 刘章君, 李丹, 陈柯兵     
武汉大学 水资源与水电工程科学国家重点实验室, 湖北 武汉 430072
基金项目: 国家自然科学基金重大资助项目(51539009);十三五国家重点研发计划资助项目(2016YFC0402206)
摘要: 对于传统的多变量分析方法,给定某一联合重现期,满足防洪标准的峰量组合形式有无数多种,但不同联合设计值得到的洪水过程线及水库最高调洪水位均存在较大的差异,如何在联合重现期水平下选取合理的设计值尤为关键。本文利用Copula函数建立洪峰与不同时段洪量的多变量联合分布,采用OR重现期作为防洪标准,联合概率密度函数值作为洪水发生相对可能性的度量指标;以满足防洪标准为约束条件,构建峰量最可能组合的通用模型,并通过拉格朗日乘数法进行求解。采用最可能组合方法分别计算了汉江流域丹江口水库的3维和4维最可能联合设计值,并与单变量同频率法、多变量同频率法设计值对比。结果表明:单变量同频率设计值达不到防洪标准,其假设洪峰与洪量完全相关,各个洪水特征量均采用单变量概率分布来描述,未能充分考虑各个特征量的内在相关性,并不符合汛期洪水发生的内在规律;最可能组合法与多变量同频率法相比,洪峰偏小,长历时洪量偏大,丹江口水库具有多年调节能力,洪量起主要作用,所以最可能组合法推求的设计洪水结果偏安全。峰量最可能设计值更符合丹江口水库的应用要求,因此,本文所提方法具有较强的统计基础,更加符合水文现象的内在规律,可用于水库设计洪水计算。
关键词: 设计洪水    联合重现期    最可能组合法    计算通式    Copula函数    
General Formula of the Most Likely Composition Method for Flood Peaks and Volumes Estimation
YIN Jiabo, GUO Shenglian, LIU Zhangjun, LI Dan, CHEN Kebing     
State Key Lab. of Water Resources and Hydropower Eng. Sci., Wuhan Univ., Wuhan 430072, China
Abstract: For conventional multivariate frequency analysis methods, there are an infinite number of combinations of flood peaks and volumes satisfying the flood prevention standard under the given joint return period (JRP); however, different joint design values would lead to differences of flood hydrograph and maximal water level which occurs in reservoir routing.Hence, how to select reasonable design values under given JRP is very important. In this paper, the Copula functions were used to construct the multivariate joint distributions of flood peaks and volumes, and the OR JRP was adopted as flood prevention standard.The normalized joint probability density function value was adopted to measure the relative occurrence likelihood of flood events, and the general model of Most Likely Composition (MLC) of flood peaks and volumes were established subject to the flood prevention standard.The Lagrange multiplier method was adopted to solve the model.The trivariate and four-dimension MLC design floods of Danjiangkou reservoir in Hanjiang basin were estimated using the proposed method, and the quantiles estimated by univariate identical frequency (UIF) and multivariate identical frequency (MIF) methods were compared with those calculated by MLC.The results indicated that the UIF estimations could not satisfy the flood prevention standard; the UIF method assumed that the flood peak is fully dependent upon the flood volume and all the variables are characterized based on univariate distributions.The UIF method cannot take the correlations of hydrologic variables into consideration, and it does not satisfy the inherent rule of flood events. It is also indicated that the MLC method had a smaller flood peak and larger long-duration volume than those by the MIF method.The Danjiangkou reservoir has the capacity of multi-year regulation and the flood volumes are mainly responsible for flood prevention safety, and thus the MLC estimation values are safer.The joint design values of MLC method are more rational according to the application requirements in Danjiangkou reservoir, and consequently the proposed method had a strong statistical basis and could reflect the inherent rules of hydrological events, and it could be used to estimate the design floods of reservoir.
Key words: design flood    joint return period    most likely composition method    general formula    Copula function    

中国《水利水电工程设计洪水计算规范》(以下简称《规范》) 规定,在推求设计洪水过程线时,采用同倍比或同频率放大方法[1]。显然,认为防洪安全标准等于设计洪水过程线某一个特征量的频率 (同倍比法) 或者多个特征量的频率 (同频率法),这仅是一种假设,它能否使得设计洪水过程线整体达到指定的防洪标准仍有待商榷[2];并且,同频率方法仅控制各个特征量的重现期分别等于洪水的整体设计标准,而忽略了洪水过程中各个特征量之间的相关性,在水文机理上也存在一定的局限性。近年来不少学者将Copula函数应用到多变量洪水频率分析中,推求了满足防洪标准的联合重现期等值面[3-4]

在联合重现期等值面上,洪峰和不同时段洪量的组合形式有无数种,存在较大的任意性和不确定性。Volpi等[5]认为并非所有的的洪水组合都符合水文事件的内在规律,如何合理地在重现期等值曲面上选取联合设计值成为研究难点。冯平等[6]基于Copula函数构造了王快水库洪峰与6 d洪量的两变量联合分布,并推求了两变量重现期下的洪水设计值;李天元等[7]构造了洪峰与时段洪量之间的多变量联合分布,确定了三变量重现期下的峰量设计值。然而,他们为了在联合重现期等值线或等值曲面上选取确定的设计标准,仍然采用了同频率假设来确定联合设计值。针对该问题,Xu等[8]推导了条件最可能组合和条件期望组合法的计算通式,但该方法仅考虑单一变量的重现期,不能适用于多变量联合重现期。近年来,多变量水文事件的最可能组合模式受到关注,刘章君等[9]研究了地区洪水的最可能组合,但是该方法将梯级水库的水力联系假定为线性方式,不能解决多变量洪水事件在联合重现期约束下的非线性难题。

为了解决在联合重现期水平下选取设计值存在的任意性和不确定性,本文构建了满足防洪标准约束的峰量最可能组合法通用模型,并通过拉格朗日乘数法进行求解,以汉江流域丹江口水库为例进行验证,为水库的设计洪水计算提供一种新思路。

1 水库洪水的联合重现期

假设Q表示洪峰流量, Wi(i=1, 2, …, n) 表示不同时段的洪量,n表示水库推求设计洪水需要的时段洪量数。按照中国《规范》的要求,本文采用P-Ⅲ线型构建洪峰、洪量的边缘分布[1],分别记为FQ(q)、FWi(wi),相应的密度函数分别记为fQ(q)、fWi(wi)。

在传统的单变量洪水频率分析框架中,洪峰或洪量的重现期一般定义为其发生频率的倒数,即

$ \left\{ \begin{align} &{{T}_{Q}}\left( q \right)=\frac{1}{1-{{F}_{Q}}\left( q \right)}, \\ &{{T}_{{{W}_{i}}}}\left( {{w}_{i}} \right)=\frac{1}{1-{{F}_{{{W}_{i}}}}\left( {{w}_{i}} \right)} \\ \end{align} \right. $ (1)

而对于多变量洪水频率分析而言,存在多种联合重现期的定义方法,目前在洪水频率分析中,主要采用AND、OR或Kendall重现期[7-8]。对于水库的防洪安全设计问题,国内外学者一般关注其中的OR重现期[5-7],其定义如下:

$ {{T}^{\cup }}\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right)=1/\left[1-F\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right) \right] $ (2)

式中:T(q, w1, w2, …, wn) 为联合重现期,a;F(q, w1, w2, …, wn) 为联合分布函数。

2 设计洪水峰量最可能组合法的数学描述

理论上,对于一定的联合重现期,满足式 (2) 约束的峰量组合 (q, w1, w2, …, wn) 有无数多个联合设计值,但不同联合设计值得到的洪水过程线及水库最高调洪洪水位存在较大的差异[5]。同频率组合法是目前使用最为广泛的方法,但该法只是一种人为假设,并不符合流域设计洪水的客观规律,难以满足水库设计洪水计算的实际需要[8]

在工程设计中,通过分析实际发生洪水的时空特性规律,通常关心最可能发生的洪水事件。因此,为了减小联合设计值选取的盲目性和任意性,本文从洪水发生可能性最大的角度,建立水库设计洪水联合重现期水平下峰量最可能组合法的通用模型。

不同峰量组合发生的相对可能性大小,可以用QWi(i=1, 2, …, n) 的联合概率密度函数值f(q, w1, w2, …, wn) 来度量。若设计值的联合概率密度函数值越大,则表明该组合发生的可能性越大。对于峰量最可能组合的设计值,需要求解f(q, w1, w2, …, wn) 在满足式 (2) 条件下的最大值,即:

$ \left\{ \begin{align} &\max \ f\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right), \\ &\text{s}\text{.t}\ \ \ \ \ {{T}^{\cup }}\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right)=1/\left[1-F\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right) \right] \\ \end{align} \right. $ (3)

本文将方程组 (3) 称为峰量最可能组合法的通用模型。QWi(i=1, 2, …, n) 的联合分布函数用F(q, w1, w2, …, wn) 表示:

$ F\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right)=P $ (4)

式中,qQ, w1W1, w2W2, …, wnWn

相应的概率密度函数为:

$ \begin{align} &f\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right)= \\ &{{\partial }^{n+1}}F\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right)/\partial q, \partial {{w}_{1}}, \partial {{w}_{2}}, \cdots, \partial {{w}_{n}} \\ \end{align} $ (5)
3 基于Copula函数推导设计洪水峰量最可能组合法的计算通式 3.1 Copula函数

Copula函数可以将若干个随机变量的边缘分布连接起来,从而构造联合分布。令H(x1, x2, …, xn) 为一个n维分布函数,其边缘分布分别为F1(x1), F2(x2), …, Fn(xn), 那么可寻求一个n维Copula函数C,使得对任意xRn(xn维向量,Rnn维实数空间)[10]

$ H\left( {{x}_{1}}, {{x}_{2}}, \cdots, {{x}_{n}} \right)={{C}_{\theta }}\left( {{F}_{1}}\left( {{x}_{1}} \right), {{F}_{2}}\left( {{x}_{2}} \right), \cdots, {{F}_{n}}\left( {{x}_{n}} \right) \right) $ (6)

式中,θ为Copula函数的相关性参数。

采用AIC信息准则和均方根误差RMSE准则评价Copula函数的拟合情况,计算公式如下[11]

$ RMSE=\sqrt{\frac{1}{L}\sum\limits_{i=1}^{L}{{{\left( {{P}_{ei}}-{{P}_{i}} \right)}^{2}}}} $ (7)
$ AIC-L\ln \left( \frac{1}{L}\sum\limits_{i=1}^{L}{{{\left( {{P}_{ei}}-{{P}_{i}} \right)}^{2}}} \right)+2m $ (8)

式中, PeiPi分别为经验频率与理论频率, L为资料系列长度, m为Copula函数的参数个数。

3.2 推导峰量最可能组合法的计算通式

借助Copula函数,可以将联合分布函数F(q, w1, w2, …, wn) 表示为[8]

$ F\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right)=C\left( u, {{v}_{1}}, {{v}_{2}}, \cdots, {{v}_{n}} \right) $ (9)

式中:C(u, v1, v2, …, vn) 为Copula函数;u=FQ(q)、vi=FWi(wi) 分别为水库的洪峰、洪量的边缘分布函数。

相应的联合概率密度函数可以表示为[8]

$ f\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right)=c\left( u, {{v}_{1}}, {{v}_{2}}, \cdots, {{v}_{n}} \right){{f}_{Q}}\left( q \right)\prod\limits_{i=1}^{n}{{{f}_{{{W}_{i}}}}}\left( {{w}_{i}} \right) $ (10)

式中, c(u, v1, v2, …, vn) 为Copula函数的密度函数, fQ(q) 与fWi(wi) 分别为QWi (i=1, 2, …, n) 的概率密度函数。

本文的峰量最可能值组合为 (Q, W1, W2, …, Wn) 在满足防洪标准的条件下,f (q, w1, w2, …, wn) 取最大值时的洪水组合模式。为此,本文提出采用拉格朗日乘数法求解最可能组合通用模型。首先构造以下求解方程:

$ \left\{ \begin{align} &\max \ f\left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right)= \\ &\ \ c\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right), {{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right), \cdots, {{F}_{{{W}_{n}}}}\left( {{w}_{n}} \right) \right){{F}_{Q}}\left( q \right)\prod\limits_{i=1}^{n}{{{f}_{{{W}_{i}}}}}\left( {{w}_{i}} \right), \\ &C\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right), {{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right), \cdots, {{F}_{{{W}_{n}}}}\left( {{w}_{n}} \right) \right)=1-1/{{T}^{\cup }} \\ \end{align} \right. $ (11)

函数f(q, w1, w2, …, wn) 在n维空间的定义域内是连续的,故必然存在最小值和最大值。给定联合重现期T,构造拉格朗日函数如下:

$ \begin{align} &\varphi \left( q, {{w}_{1}}, {{w}_{2}}, \cdots, {{w}_{n}} \right)=c\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right), {{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right), \cdots, \right. \\ &\ \ \ \left. {{F}_{{{W}_{n}}}}\left( {{w}_{n}} \right) \right)\cdot {{f}_{Q}}\left( q \right)\prod\limits_{i=1}^{n}{{{f}_{{{W}_{i}}}}}\left( {{w}_{i}} \right)+\lambda \left[C\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right), \right. \right. \\ &\ \ \ \ \ \ \ \ \left. {{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right), \cdots, \left. {{F}_{{{W}_{n}}}}\left( {{w}_{n}} \right) \right)-1+1/{{T}^{\cup }} \right] \\ \end{align} $ (12)

分别对q, w1, w2, …, wnλ求偏导,并令其为0,即可求得拉格朗日函数φ(q, w1, w2, …, wn) 的所有极值点,由于式 (12) 的求解比较复杂,遵循由低维到高维的求解思路,先研究n=1~3的情形,再通过数学归纳法求解n个洪量时峰量最可能组合法的计算通式。

1)n=1。由式 (12) 求偏导,化简后得到如下方程组:

$ \left\{ \begin{align} &{{\varphi }_{q}}=\left( {{c}_{q}}f_{Q}^{2}\left( q \right)+cf_{Q}^{2}\left( q \right)+c{{f}_{Q}}^{'}\left( q \right) \right){{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right)+ \\ &\ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right) \right)}{\partial q}{{f}_{Q}}\left( q \right)=0, \\ &{{\varphi }_{{{w}_{1}}}}={{f}_{Q}}\left( q \right)\cdot \left( {{c}_{{{w}_{1}}}}f_{{{W}_{1}}}^{2}\left( {{w}_{1}} \right)+c{{f}_{{{W}_{1}}}}^{'}\left( {{w}_{1}} \right) \right)+ \\ &\ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right) \right)}{\partial {{w}_{1}}}{{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right)=0, \\ &{{\varphi }_{\lambda }}=C\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right) \right)-{{p}_{T}}=0 \\ \end{align} \right. $ (13)

式中:λ为拉格朗日乘子;$c=c\left( u, {{v}_{1}} \right), {{c}_{q}}=\partial c/\partial u, {{c}_{{{w}_{1}}}}=\partial c/\partial {{v}_{1}};{{p}_{T}}=1-1/{{T}^{\cup }};{{f}_{Q}}^{'}、{{f}_{{{W}_{1}}}}^{'}$分别为相应概率密度函数的导函数。

2)n=2。化简后得到如下方程组:

$ \left\{ \begin{align} & {{\varphi }_{q}}=\left( {{c}_{q}}f_{Q}^{2}\left( q \right)+c{{f}_{Q}}^{'}\left( q \right) \right){{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right){{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right)+ \\ & \ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right),{{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right),{{F}_{{{W}_{2}}}}\left( {{w}_{2}} \right) \right)}{\partial q}{{f}_{Q}}\left( q \right)=0, \\ & {{\varphi }_{{{w}_{1}}}}={{f}_{Q}}\left( q \right){{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right)\cdot \left( {{c}_{{{w}_{1}}}}f_{{{W}_{1}}}^{2}\left( {{w}_{1}} \right)+c{{f}_{{{W}_{1}}}}^{'}\left( {{w}_{1}} \right) \right)+ \\ & \ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right),{{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right) \right)}{\partial {{w}_{1}}}{{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right)=0, \\ & {{\varphi }_{{{w}_{2}}}}={{f}_{Q}}\left( q \right){{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right)\cdot \left( {{c}_{{{w}_{2}}}}f_{{{W}_{2}}}^{2}\left( {{w}_{2}} \right)+c{{f}_{{{W}_{2}}}}^{'}\left( {{w}_{2}} \right) \right)+ \\ & \ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right),{{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right),{{F}_{{{W}_{2}}}}\left( {{w}_{2}} \right) \right)}{\partial {{w}_{2}}}{{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right)=0, \\ & {{\varphi }_{\lambda }}=C\left( {{F}_{Q}}\left( q \right),{{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right),{{F}_{{{W}_{2}}}}\left( {{w}_{2}} \right) \right)-{{p}_{T}}=0 \\ \end{align} \right. $ (14)

式中:c=c(u, v1, v2)= 3C(u, v1, v2)/ u v1 v2cw2= c/ v2, fW2′为fW2的导函数,其他与前文相同。

3)n=3。化简后得到如下方程组:

$ \left\{ \begin{align} &{{\varphi }_{q}}=\left( {{c}_{q}}f_{Q}^{2}\left( q \right)+c{{f}_{Q}}^{'}\left( q \right) \right){{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right){{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right){{f}_{{{W}_{3}}}}\left( {{w}_{3}} \right)+ \\ &\ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right), {{F}_{{{W}_{2}}}}\left( {{w}_{2}} \right), {{F}_{{{W}_{3}}}}\left( {{w}_{3}} \right) \right)}{\partial q}. \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{f}_{Q}}\left( q \right)=0, \ \ \\ &{{\varphi }_{{{w}_{1}}}}={{f}_{Q}}\left( q \right){{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right){{f}_{{{W}_{3}}}}\left( {{w}_{3}} \right)\cdot \left( {{c}_{{{w}_{1}}}}f_{{{W}_{1}}}^{2}\left( {{w}_{1}} \right)+c{{f}_{{{W}_{1}}}}^{'}\left( {{w}_{1}} \right) \right)+ \\ &\ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right), {{F}_{{{W}_{2}}}}\left( {{w}_{2}} \right), {{F}_{{{W}_{3}}}}\left( {{w}_{3}} \right) \right)}{\partial {{w}_{1}}}. \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right)=0, \\ &{{\varphi }_{{{w}_{2}}}}={{f}_{Q}}\left( q \right){{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right){{f}_{{{W}_{3}}}}\left( {{w}_{3}} \right)\cdot \left( {{c}_{{{w}_{2}}}}f_{{{W}_{2}}}^{2}\left( {{w}_{2}} \right)+c{{f}_{{{W}_{2}}}}^{'}\left( {{w}_{2}} \right) \right)+ \\ &\ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right), {{F}_{{{W}_{2}}}}\left( {{w}_{2}} \right), {{F}_{{{W}_{3}}}}\left( {{w}_{3}} \right) \right)}{\partial {{w}_{2}}}. \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ {{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right)=0, \\ &{{\varphi }_{{{w}_{3}}}}={{f}_{Q}}\left( q \right){{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right){{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right)\cdot \left( {{c}_{{{w}_{3}}}}f_{{{W}_{3}}}^{2}\left( {{w}_{3}} \right)+c{{f}_{{{W}_{3}}}}^{'}\left( {{w}_{3}} \right) \right)+ \\ &\ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right), {{F}_{{{W}_{2}}}}\left( {{w}_{2}} \right), {{F}_{{{W}_{3}}}}\left( {{w}_{3}} \right) \right)}{\partial {{w}_{3}}}. \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ {{f}_{{{W}_{3}}}}\left( {{w}_{3}} \right)=0, \\ &{{\varphi }_{\lambda }}=C\left( {{F}_{Q}}\left( q \right), {{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right), {{F}_{{{W}_{2}}}}\left( {{w}_{2}} \right), {{F}_{{{W}_{3}}}}\left( {{w}_{3}} \right) \right)-{{p}_{T}}=0 \\ \end{align} \right. $ (15)

式中,c=c(u, v1, v2, v3)= 4C(u, v1, v2, v3)/ u v1 v2 v3cw3= c/ v3, fW3′为fW3的导函数,其他与前文相同。

由以上结果看出,n=1~3时峰量最可能组合满足的方程组有明显规律,通过数学归纳法,可以得到n为任意正整数时的方程组通式:

$ \left\{ \begin{align} & {{\varphi }_{q}}=\left( {{c}_{q}}f_{Q}^{2}\left( q \right)+c{{f}_{Q}}^{'}\left( q \right) \right)\prod\limits_{i=1}^{n}{{{f}_{{{W}_{i}}}}\left( {{w}_{i}} \right)}+ \\ & \ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right),{{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right),{{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right),\cdots ,{{F}_{{{W}_{n}}}}\left( {{w}_{n}} \right) \right)}{\partial q}{{f}_{Q}}\left( q \right)=0, \\ & {{\varphi }_{{{w}_{1}}}}={{f}_{Q}}\left( q \right)\prod\limits_{i=2}^{n}{{{f}_{{{W}_{i}}}}\left( {{w}_{i}} \right)}\cdot \left( {{c}_{{{w}_{1}}}}f_{{{W}_{1}}}^{2}\left( {{w}_{1}} \right)+c{{f}_{{{W}_{1}}}}^{'}\left( {{w}_{1}} \right) \right)+ \\ & \ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right),{{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right),{{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right),\cdots ,{{F}_{{{W}_{n}}}}\left( {{w}_{n}} \right) \right)}{\partial {{w}_{1}}}{{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right)=0, \\ & {{\varphi }_{{{w}_{2}}}}={{f}_{Q}}\left( q \right){{f}_{{{W}_{1}}}}\left( {{w}_{1}} \right)\prod\limits_{i=3}^{n}{{{f}_{{{W}_{i}}}}\left( {{w}_{i}} \right)}\cdot \left( {{c}_{{{w}_{2}}}}f_{{{W}_{2}}}^{2}\left( {{w}_{2}} \right)+c{{f}_{{{W}_{2}}}}^{'}\left( {{w}_{2}} \right) \right)+ \\ & \ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right),{{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right),{{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right),\cdots ,{{F}_{{{W}_{3}}}}\left( {{w}_{3}} \right) \right)}{\partial {{w}_{2}}}{{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right)=0, \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdots \\ & {{\varphi }_{{{w}_{n-1}}}}={{f}_{Q}}\left( q \right){{f}_{{{W}_{n}}}}\left( {{w}_{n}} \right)\prod\limits_{i=1}^{n-2}{{{f}_{{{W}_{i}}}}\left( {{w}_{i}} \right)}\cdot \left( {{c}_{{{w}_{n-1}}}}f_{{{W}_{n-1}}}^{2}\left( {{w}_{n-1}} \right)+ \right. \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. c{{f}_{{{W}_{n-1}}}}^{'}\left( {{w}_{n-1}} \right) \right)+ \\ & \ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right),{{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right),{{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right),\cdots ,{{F}_{{{W}_{n}}}}\left( {{w}_{n}} \right) \right)}{\partial {{w}_{n-1}}}{{f}_{{{W}_{n-1}}}}\left( {{w}_{n-1}} \right)=0 \\ & {{\varphi }_{{{w}_{n}}}}={{f}_{Q}}\left( q \right)\prod\limits_{i=1}^{n-1}{{{f}_{{{W}_{i}}}}\left( {{w}_{i}} \right)}\cdot \left( {{c}_{{{w}_{n}}}}f_{{{W}_{n}}}^{2}\left( {{w}_{n}} \right)+c{{f}_{{{W}_{n}}}}^{'}\left( {{w}_{n}} \right) \right)+ \\ & \ \ \ \ \ \ \ \ \lambda \frac{\partial C\left( {{F}_{Q}}\left( q \right),{{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right),{{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right),\cdots ,{{F}_{{{W}_{n}}}}\left( {{w}_{n}} \right) \right)}{\partial {{w}_{n}}}{{f}_{{{W}_{n}}}}\left( {{w}_{n}} \right)=0, \\ & {{\varphi }_{\lambda }}=C\left( {{F}_{Q}}\left( q \right),{{F}_{{{W}_{1}}}}\left( {{w}_{1}} \right),{{f}_{{{W}_{2}}}}\left( {{w}_{2}} \right),\cdots ,{{F}_{{{W}_{n}}}}\left( {{w}_{n}} \right) \right)-{{p}_{T}}=0 \\ \end{align} \right. $ (16)

非线性方程组 (16) 为联合重现期水平下最可能组合法的计算通式。式 (16) 有 (n+2) 个未知数λ, q, w1, w2, …, wn和 (n+2) 个方程,根据问题的实际意义,f(q, w1, w2, …, wn) 函数的最大值存在且唯一,因此该方程组必定有唯一解,因此由拉格朗日乘数法得到的极值点即为最大值点[9]

针对式 (16),可采用调和平均数牛顿法[12]进行迭代求解,从而得到水库设计洪水峰量最可能组合 (q*, w1*, w2*, …, wn*)。

4 实例研究

汉江流域位于东经106°15′~114°20′、北纬30°10′~34°20′。多年平均降水量897.2 mm,由上游向下游增大,暴雨集中发生在7至9月,尤以7、9月居多,5至10月径流量占全年的80%。汉江洪水由暴雨形成,集中在夏、秋两季。丹江口水库是南水北调中线工程的水源水库,也是治理开发汉江的关键性工程。大坝以上控制流域面积9.52×104 km2,具有防洪、供水、发电和航运等综合效益。

本文选用丹江口水库1954—2014年61 a的实测入库流量系列,研究其设计洪水峰量联合设计值,并分别采用单变量同频率、多变量同频率和最可能组合等三种方法估计丹江口水库的设计洪水成果。

4.1 确定峰量的边缘分布

采用年最大值取样方法,得到丹江口水库年最大洪峰Qmax、年最大3 d洪量W3d、年最大7 d洪量W7d和年最大15 d洪量W15d。采用《规范》推荐的方法分别估计P-Ⅲ分布函数的参数,结果见表 1。采用χ2检验法对其进行假设检验,在5%的显著性水平下,自由度为k-r-1(kχ2检验的分组数,r为参数个数) 的χ2检验接受域为小于等于临界值,表 1χ0.052为临界值,从表 1可以看出4个随机变量的P-Ⅲ型分布均通过了检验。图 1给出了丹江口水库QmaxW3dW7dW15d的频率曲线,从图 1可以看出实测值与理论曲线的拟合效果较好。

表1 丹江口水库洪峰洪量系列统计特征值和P-Ⅲ型分布参数估计结果 Tab. 1 Estimated sample statistic values and P-Ⅲ distribution parameters in the Danjiangkou reservoir

图1 丹江口水库年最大洪峰、3 d洪量、7 d洪量和15 d洪量的P-Ⅲ分布频率曲线图 Fig. 1 Plots of P-Ⅲ probability curves with observed flood peak, 3 d-flood, 7 d-flood and 15 d-flood volume

4.2 建立多维Copula联合分布函数

Archimedean族的Copula函数适用于变量间存在正相关或负相关的情形,并且能描述水文变量间的相关性,其中的χ0.052Gumbel-Hougaard (G-H) Copula函数在描述洪水峰量关系中应用广泛[11, 13-14]。采用Kendall秩相关性系数[9-10]描述丹江口水库QmaxW3dW7dW15d的相关性,分别为0.80(Qmax-W3d),0.76(Qmax-W7d),0.73(Qmax-W15d),0.86(W3d-W7d),0.79(W3d-W15d) 和0.85(W7d-W15d)。可以看出:随着洪水历时的增加,洪峰与洪量的相关性逐渐减弱;而不同洪量的相关性随着历时差的增大而减弱。本文采用表 2中6种不同的Copula函数建立QmaxW3dW7d的3维联合分布;分别采用4维非对称和对称型G-H函数建立QmaxW3dW7dW15d的联合分布函数; 采用极大似然法[13]估计参数,各类Copula函数的函数表达式详见文献[14]。参数估计结果见表 2

表2 Copula函数参数的估计和检验结果 Tab. 2 Estimated parameter and evaluation results for Copula function

为了检验所选取的Copula函数是否恰当,通过K-S检验法进行拟合检验,统计量D值越小,则说明所选Copula函数的拟合效果更优;利用AIC信息准则与均方根误差RMSE准则进行拟合优度评价,AIC值和RMSE值越小,则说明相应的Copula函数越优。表 2也给出了3维、4维对称和非对称Copula函数的K-S检验统计量DAIC值和RMSE值。取显著性水平为α=0.05,与n=61对应的分位数为0.212 4,若D值小于0.212 4则Copula函数通过K-S检验。

表 2可知,本文考虑的各类Copula函数均能通过检验,3维非对称型G-H Copula函数的D值、RMSE值和AIC值均小于其他5类3维函数;4维非对称型G-H Copula函数的统计值小于4维对称型G-H Copula函数。图 23分别给出了3维、4维非对称型G-H Copula函数计算得出的理论联合分布与经验联合分布值;图 45以另一种形式展示了3维、4维非对称型G-H Copula联合函数的理论联合分布值与经验联合分布值。从图 2~5可以看出,3维、4维非对称型Copula联合函数所建立的联合分布均拟合较好,都是合理可行的。再考虑到一般情况下构造多变量峰量关系较多采用非对称Copula函数[3-4, 15],因此选择3维、4维非对称G-H Copula函数构造丹江口水库峰量的多变量联合分布。

图2 3维非对称型G-H Copula函数的理论频率和联合观测值的经验频率相关关系 Fig. 2 Correlation between the empirical and theoretical values of trivariate asymmetrical G-H Copul

图3 4维非对称型G-H Copula函数的理论频率和联合观测值的经验频率相关关系 Fig. 3 Correlation between the empirical and theoretical values of four-dimensional asymmetric G-H Copula

图4 3维非对称型G-H Copula函数的理论频率和联合观测值的经验频率比较 Fig. 4 Comparison of the empirical plots and theoretical curve of four-dimensional asymmetrical G-H Copula

图5 4维非对称型G-H Copula函数的理论频率和联合观测值的经验频率比较 Fig. 5 Comparison of the empirical plots and theoretical curve of four-dimensional asymmetric G-H Copula

4.3 推求联合设计值 4.3.1 3维联合设计值

采用单变量同频率组合法和多变量同频率组合法[7]推求了丹江口水库Qmax-W3d-W7d的联合设计值,计算结果见表 3。令非线性方程组式 (16) 中n=2, 即为洪峰与两个控制时段洪量的情形, 为了比较不同初值对本文模型及求解方法的影响,分别设置单变量同频率和多变量同频率设计值为初值,得到峰量最可能组合法的联合设计值,计算结果也列于表 3。可以看出,两种初值方案下的最可能联合设计值没有多大区别,说明本文求解方法的稳健性强,对初值的要求不高。但是,采用多变量同频率法为初值,求解的收敛速度较快。

表3 丹江口水库三变量联合设计洪水值 Tab. 3 Results of trivariate joint design flood values in Danjiangkou reservoir

表 3中还可以看出,单变量同频率法估计结果与多变量同频率法设计值相比明显偏小,这说明现阶段采用单变量同频率方法估计设计洪水,降低了防洪标准。实际上,传统的单变量同频率法假设洪峰与洪量完全相关,各个洪水特征量均采用单变量概率分布来描述,该方法未能充分考虑各个特征量的内在相关性,并不符合汛期洪水发生的内在规律。

最可能组合法与多变量同频率法的联合设计值相比,洪峰偏小,洪量偏大。对于围堰或小水库,其调洪结果主要取决于洪峰流量;而对于调节能力大的水库,洪量起主要作用[9, 16]。丹江口水库具有多年调节能力,所以最可能组合法推求的设计洪水结果偏安全,更符合丹江口水库的应用实践。

4.3.2 4维联合设计值

为了分析本文递推通式的可靠性,令非线性方程式 (16) 中n=3,以多变量同频率设计值为初始值,求解了丹江口水库Qmax-W3d-W7d-W15d的联合设计值,结果如表 4所示。

表4 丹江口水库4维联合设计洪水值 Tab. 4 Results of four-dimensional joint design flood values in Danjiangkou reservoir

表 4中可以看出,最可能组合法的设计值与多变量同频率法的设计值相比,洪峰和3 d洪量偏小,7、15 d洪量均偏大,符合联合重现期的估计规律。与3维情形类似,长历时洪量偏大,丹江口水库的调洪结果偏于安全[16]

5 结论

1) 单变量同频率法设计值达不到防洪标准,该法假设洪水的峰量完全相关,未能充分考虑水文特征量之间的内在相关性,不符合洪水发生的内在规律。

2) 峰量最可能组合法适用于洪水峰量具有任意相关的情形,且组合方案唯一。利用Copula函数建立洪水峰量的联合分布,采用联合概率密度函数作为洪水发生相对可能性的度量指标,通过拉格朗日乘数法推导了联合重现期水平下的峰量最可能组合法的计算通式。

3) 应用汉江流域丹江口水库实例进行验证,结果显示,最可能组合法与多变量同频率法相比,推求的洪峰设计值偏小,长历时洪量偏大。该法具有较强的统计基础,更加符合水文现象的内在规律和工程实际的要求,为水库的设计洪水计算提供了一条新思路。

参考文献
[1]
Guo Shenglian, Liu Zhangjun, Xiong Lihua. Advances and assessment on design flood estimation methods[J]. Journal of Hydraulic Engineering, 2016, 47(3): 302-314. [郭生练, 刘章君, 熊立华. 设计洪水计算方法研究进展与评价[J]. 水利学报, 2016, 47(3): 302-314.]
[2]
Qin Guanghua, Song Kechao, Zhou Zejiang, et al. Research on annual runoff prediction based on WA-GRNN model[J]. Journal of Sichuan University (Engineering Science Edition), 2013, 45(6): 39-46. [覃光华, 宋克超, 周泽江, 等. 基于WA-GRNN模型的年径流预测[J]. 四川大学学报 (工程科学版), 2013, 45(6): 39-46.]
[3]
Sraj M, Bezak N, Brilly M. Bivariate flood frequency analysis using the copula function:A case study of the Litija station on the Sava River[J]. Hydrology Process, 2015, 29(2): 225-238. DOI:10.1002/hyp.v29.2
[4]
Ganguli P, Reddy M J. Probabilistic assessment of flood risks using trivariate copulas[J]. Theoretical and Applied Climatology, 2013, 111(1/2): 341-360.
[5]
Volpi E, Fiori A. Design event selection in bivariate hydrological frequency analysis[J]. Hydrological Sciences Journal, 2012, 57(8): 1506-1515. DOI:10.1080/02626667.2012.726357
[6]
Feng Ping, Li Xin. Bivariate frequency analysis of non-stationary flood time series based on Copula methods[J]. Journal of Hydraulic Engineering, 2013, 44(10): 1137-1147. [冯平, 李新. 基于Copula函数的非一致性洪水峰量联合分析[J]. 水利学报, 2013, 44(10): 1137-1147.]
[7]
Li Tianyuan, Guo Shenglian, Yan Baowei, et al. Derivative design flood hydrograph based trivariate joint distribution[J]. Journal of Hydroelectric Engineering, 2013, 32(3): 10-14. [李天元, 郭生练, 闫宝伟, 等. 基于多变量联合分布推求设计洪水过程线的新方法[J]. 水力发电学报, 2013, 32(3): 10-14.]
[8]
Xu Changjiang, Yin Jiabo, Guo Shenglian, et al.Deriving design flood hydrograph based on conditional distribution:A case study of danjiangkou reservoir in Hanjiang Basin[J/OL].Mathematical Problems in Engineering, 2016.http://dx.doi.org/10.1155/2016/4319646.
[9]
Liu Zhangjun, Guo Shenglian, Li Tianyuan, et al. General formula derivation of most likely regional composition method for design flood estimation of cascade reservoirs system[J]. Advances in Water Science, 2014, 25(4): 575-584. [刘章君, 郭生练, 李天元, 等. 梯级水库设计洪水最可能地区组成法计算通式[J]. 水科学进展, 2014, 25(4): 575-584.]
[10]
Xue Jinping, Hu Zhigen, Liu Quan. Risk analysis of the river diversion under construction of cascade hydropower stations[J]. Journal of Sichuan University (Engineering Science Edition), 2014, 46(1): 75-80. [薛进平, 胡志根, 刘全. 梯级水电站建设施工导流风险分析[J]. 四川大学学报 (工程科学版), 2014, 46(1): 75-80.]
[11]
Zhang Q, Xiao M Z, Singh V P. Uncertainty evaluation of copula analysis of hydrological droughts in the East River basin, China[J]. Global Planet Change, 2015, 129: 1-9. DOI:10.1016/j.gloplacha.2015.03.001
[12]
Özban A Y. Some new variants of Newton's method[J]. Applied Mathematical Letter, 2004, 17: 677-682. DOI:10.1016/S0893-9659(04)90104-8
[13]
Zhang Zhenghao, Zhang Qiang, Xiao Mingzhong, et al. Reservoir-based ecological water operation in the Liaohe river basin characterized by synchronous occurrence of wet/dry events[J]. Acta Ecologica Sinica, 2016, 36(7): 1-10. [张正浩, 张强, 肖名忠, 等. 辽河流域丰枯遭遇下水库调度[J]. 生态学报, 2016, 36(7): 1-10.]
[14]
Song Songbai, Nie Rong. Asymmetric archimedean Copulas for multivariate hydrological drought frequency analysis[J]. Journal of Hydroelectric Engineering, 2011, 30(4): 20-29. [宋松柏, 聂荣. 基于非对称阿基米德Copula的多变量水文干旱联合概率研究[J]. 水力发电学报, 2011, 30(4): 20-29.]
[15]
Wu Tao, Chen Xi, Yan Yusong. Study of the binary correlation quantum-behaved PSO algorithm[J]. Journal of Sichuan University (Engineering Science Edition), 2014, 46(4): 103-110. [吴涛, 陈曦, 严余松. 2元相关性量子行为粒子群优化算法研究[J]. 四川大学学报 (工程科学版), 2014, 46(4): 103-110.]
[16]
Zhou Ying, Song Dedun. A model for the optimal design of flood control reservoir[J]. Advances in Water Science, 1994, 5(3): 167-173. [邹鹰, 宋德敦. 水库防洪优化设计模型[J]. 水科学进展, 1994, 5(3): 167-173.]