1. 引言

随着我国各大流域联合防洪调度工作的逐渐深入，流域性洪水引发的多条河流的洪水遭遇分析工作越来越迫切。洪水遭遇风险研究对下游的防洪安全具有重要意义。传统的多变量水文分析计算方法主要包括多元正态分布 [1] [2] 、特定边缘分布构成的联合分布 [3] [4] 、非参数方法 [5] 、将多维联合分布转换成一维分布的方法 [6] 等。传统方法往往存在着一些缺陷，比如多元正态分布在数据转换时容易造成信息失真，非参数方法构造的联合分布的边缘分布类型未知等 [7]。目前，国内外针对上游有支流汇入的防洪控制点的风险估算，多是针对一条支流汇入的情况，对于多条河流汇入的情况研究较少；并且多是将干流洪水的随机性作为研究对象，区间洪水设为特定值进行计算，这样显然与实际有所偏差。因此，研究新的方法，充分考虑干支流洪水的随机性，对于上游有多条支流汇入的防洪控制点的防洪安全具有重要意义。

频率组合法是研究洪水遭遇的传统方法，目前多是针对二维水文随机变量进行实例研究，对于多维水文随机变量的研究相对较少，当涉及的维数较多时，用传统的解析方法来确定联合概率密度函数的解析式以及对其积分求解是非常困难的，因此，通常采用离散叠加法进行数值求解。但是，干支流洪水变量间往往存在相关性，在采用离散叠加法进行求解时，需要对变量进行独立性处理，此过程难免会出现数据失真。

Copula函数能够通过变量的边缘分布和相关性结构两部分来构造多维联合分布，形式灵活多样，是构建联合分布的一种非常有效的方法，近些年在水文分析计算领域中得到了广泛应用。熊立华等 [8] 采用Copula函数构造了两站年最大洪水联合分布函数，并且对经验概率和理论概率进行了拟合。李天元等 [9] 通过双参数Copula函数构造了两江洪峰之间的联合分布函数，估计了洪峰遭遇风险的条件概率。刘章君等 [10] 根据Copula函数对设计洪水地区组成进行区间估计，计算出地区组成的点估计值，并对估计的不确定性进行评价。

为此，本文结合Copula函数易于构造联合分布的特点，将Copula函数在构造联合分布函数的优势与频率组合法在计算多维联合分布函数的功能相结合，以雅砻江、鲜水河和庆大河洪水遭遇为例，进行了验证分析。相对于传统的频率组合法，此方法可以充分考虑干支流的随机性及相关性，且计算过程中不会出现信息失真；对变量边缘分布类型没有要求，适用性强；解决了频率组合法在多维问题上的求解问题，使其不再局限于求解二维问题。

2. Copula函数

Copula函数由Sklar在1959年提出，是定义在[0,1]上均匀分布的联合分布函数。其中Archimedean Copula函数结构简单，适应性强，在水文多变量联合分布中的应用较为广泛。n维Archimedean Copula函数表达式如下：

$C\left(u,\phi \right)=\phi \left\{{\phi }^{-1}\left(u_1\right)+{\phi }^{-1}\left(u_2\right)+\cdots +{\phi }^{-1}\left(u_n\right)\right\},u\in \left[0,1\right]$ (1)

其中，算子 $\phi$ 决定了Copula函数的类型。双参数三维Archimedean Copula函数，表达式如下：

$C\left(u_1,u_2,u_3\right)={\phi }_1\left({\phi }_1^{-1}\left(u_1\right)+{\phi }_1^{-1}{\phi }_2\left({\phi }_2^{-1}\left(u_2\right)+{\phi }_2^{-1}\left(u_3\right)\right)\right)$ (2)

双参数三维Archimedean Copula函数拥有两个发生元，因此可以较好地刻画出变量之间的相依结构。三维GH Copula、Clayton Copula、Frank Copula函数表达式分别如下：

$C\left(u_1,u_2,u_3;{\theta }_1,{\theta }_2\right)=\exp \left\{-{\left({\left[{\left(-\log u_3\right)}^{{\theta }_2}+{\left(-\log u_2\right)}^{{\theta }_2}\right]}^{{\theta }_1/{\theta }_2}+{\left(-\log u_1\right)}^{{\theta }_1}\right)}^{1/{\theta }_1}\right\}$ (3)

$C\left(u_1,u_2,u_3;{\theta }_1,{\theta }_2\right)={\left[{\left(u_3^{-{\theta }_2}+u_2^{-{\theta }_2}-1\right)}^{{\theta }_1/{\theta }_2}+u_1^{-{\theta }_1}-1\right]}^{-1/{\theta }_1}$ (4)

$C\left(u_1,u_2,u_3;{\theta }_1,{\theta }_2\right)=-{\theta }_1^{-1}\log \left\{1-c_1^{-1}\left(1-{\left[1-c_2^{-1}\left(1-{\text{e}}^{-{\theta }_2{u}_3}\right)\left(1-{\text{e}}^{-{\theta }_2{u}_2}\right)\right]}^{{\theta }_1/{\theta }_2}\right)\left(1-{\text{e}}^{-{\theta }_1{u}_1}\right)\right\}$ (5)

其中， $c_1=1-{\text{e}}^{-{\theta }_1}$， $c_2=1-{\text{e}}^{-{\theta }_2}$，参数 ${\theta }_1$ 用来刻画 $\left(u_1,u_2\right)$ 、 $\left(u_1,u_3\right)$ 之间的相依性，参数 ${\theta }_2$ 用来刻画 $\left(u_2,u_3\right)$ 之间的相依性。

3. 基于Copula函数的洪水遭遇风险估计模型

当防洪控制点上游有多条区间洪水时，由于区间洪水这一不确定性因素，必须考虑区间洪水与干流洪水组合后带来的风险，其组合形式有很多种，为下游防洪控制点的防洪安全带来了极大的威胁 [11] [12] [13]。

如图1所示，干流洪水为 $X_1$，区间洪水为 $X_2,X_3,\cdots ,X_n$，防洪控制点洪水为，则 $W=X_1+X_2+\cdots +X_n$。 $X_1,X_2,\cdots ,X_n$ 之间可以是相关关系或相互独立。 $n-1$ 条区间洪水在防洪控制点上游相继汇入，为下游防洪控制点的防洪安全带来了一定的威胁，为保证下游防洪控制点的防洪安全，需要求出下游防洪控制点组合洪水的分布函数，根据频率组合法，组合洪水洪峰流量超越下游安全泄量的概率分布函数可以表示为：

$P\left(w\ge w_0\right)=1-{\displaystyle \int \cdots }{\displaystyle {\int }_{D}h\left(x_1,x_2,\cdots ,x_n\right)\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_n}$ (6)

式中，D为 $x_1+x_2+\cdots +x_n\le w_0$， $h\left(x_1,x_2,\cdots ,x_n\right)$ 是随机变量 $x_1,x_2,\cdots ,x_n$ 的联合概率密度函数。当 $w_0$ 为下游防洪控制点的安全泄量时，防洪控制点发生洪水的概率 $P\left(w\ge w_0\right)$ 即为防洪风险率。

公式(6)中的联合概率密度函数 $h\left(x_1,x_2,\cdots ,x_n\right)$ 可以用Copula函数来表示，公式如下：

$H\left(x_1,x_2,\cdots ,x_n\right)=C\left[F\left(x_1\right),F\left(x_2\right),\cdots ,F\left(x_n\right)\right]$ (7)

$h\left(x_1,x_2,\cdots ,x_n\right)=\frac{{\partial }^{n}H\left(x_1,x_2,\cdots ,x_n\right)}{\partial x_1\partial x_2\cdots \partial x_n}=\frac{{\partial }^{n}C\left(u_1,u_2,\cdots ,u_n\right)}{\partial x_1\partial x_2\cdots \partial x_n}$ (8)

式中， $H\left(x_1,x_2,\cdots ,x_n\right)$ 为随机变量 $x_1,x_2,\cdots ,x_n$ 的联合分布函数， $C\left[F\left(x_1\right),F\left(x_2\right),\cdots ,F\left(x_n\right)\right]$ 为随机变量 $x_1,x_2,\cdots ,x_n$ 的Copula函数， $u_1=F\left(x_1\right),u_2=F\left(x_2\right),\cdots ,u_n=F\left(x_n\right)$ 分别为 $x_1,x_2,\cdots ,x_n$ 的边缘分布函数。其中，

$\frac{{\partial }^{n}C\left(u_1,u_2,\cdots ,u_n\right)}{\partial x_1\partial x_2\cdots \partial x_n}=\frac{{\partial }^{n}C\left(u_1,u_2,\cdots ,u_n\right)}{\partial u_1\partial u_2\cdots \partial u_n}\cdot \frac{\partial u_1}{\partial x_1}\cdot \frac{\partial u_2}{\partial x_2}\cdots \frac{\partial u_n}{\partial x_n}=c\left(u_1,u_2,\cdots ,u_n\right)\cdot f_{x_1}\left(x_1\right)\cdot f_{x_2}\left(x_2\right)\cdots f_{x_n}\left(x_n\right)$ (9)

式中， $c\left(u_1,u_2,\cdots ,u_n\right)$ 为 $C\left(u_1,u_2,\cdots ,u_n\right)$ 的密度函数， $f_{x_1}\left(x_1\right),f_{x_2}\left(x_2\right),\cdots ,f_{x_n}\left(x_n\right)$ 分别为 $x_1,x_2,\cdots ,x_n$ 的边缘密度函数。

由上述公式可知，下游防洪控制点组合洪水的分布函数可由Copula函数密度函数以及边缘分布密度函数表示，因此利用Copula函数易于将多个边缘分布构造成多维联合分布的特点，对公式(6)进行改进可以得到新的n维组合洪水的分布函数表达式，如式(10)所示：

$P\left(w\ge w_0\right)=1-{\displaystyle \int \cdots }{\displaystyle {\int }_{D}c\left(u_1,u_2,\cdots ,u_n\right)\cdot f_{x_1}\left(x_1\right)\cdot f_{x_2}\left(x_2\right)\cdots f_{x_n}\left(x_n\right)\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_n}$ (10)

式(10)具有以下优点：对于组合洪水问题推广至n维，适用性较强；随机变量的联合概率密度函数用Copula函数表示，可以充分考虑变量间的相关性；对于边缘分布类型没有限制，任意边缘分布都可以构造联合分布。

函数求解步骤如下：

①确定边缘分布，即洪水流量的分布函数。

②在分析边缘分布相关性的基础上，选择合适的Copula函数。

③对Copula函数密度函数以及边缘分布密度函数的乘积求n重积分。

式(10)具有以下优点：对于组合洪水问题推广至n维，适用性较强；随机变量的联合概率密度函数用边缘分布密度函数以及Copula函数密度函数表示，可以充分考虑变量间的相关性；对于边缘分布类型没有限制，任意边缘分布都可以构造联合分布。

4. 实例分析

雅砻江流域位于青藏高原东部，整个流域呈南北向条带状，流域平均长度约950 km，平均宽度约137 km，河系呈羽状发育。支流鲜水河与庆大河相继汇入雅砻江，其流域图如图2所示，对下游防洪控制点的防洪安全造成了一定的威胁，本文采用雅砻江新龙站、鲜水河道孚站、庆大河扎巴站历史洪水资料，选取同一时刻三个测站的历史流量组成三元序列对作为样本，充分考虑洪水遭遇的可能性，对下游组合洪水分布函数以及防洪风险率进行计算。

4.1. 边缘分布的确定

P-III分布曲线在水文变量频率分析中应用较为普遍，其对径流、暴雨的拟合结果较好。其概率密度函数如下：

$f\left(x\right)=\frac{{\beta }^{\gamma }}{\Gamma \left(\gamma \right)}{\left(x-\alpha \right)}^{\gamma -1}\exp \left[-\beta \left(x-\alpha \right)\right]$ (11)

其中， $\alpha$ 、 $\beta$ 、 $\gamma$ 分别为P-III型分布的位置、尺度、形状参数，计算公式如下：

$\alpha =\frac{4}{C_s^2},\beta =\frac{2}{\bar{x}{C}_v{C}_s},\gamma =\bar{x}\left(1-\frac{2C_v}{C_s}\right)$ (12)

根据上述公式，可以计算出雅砻江洪峰、鲜水河洪峰和庆大河洪峰的P-III分布参数，由表1可知，雅砻江、鲜水河和庆大河洪水流量均服从P-III分布，以P-III分布作为雅砻江、鲜水河和庆大河洪水流量的边缘分布是合理可行的。

4.2. Copula函数类型的选择

首先，对雅砻江洪峰、鲜水河洪峰、庆大河洪峰进行相关性分析，判断其相关性，可以初步确定适合本文研究所应用的Copula函数类型。本文应用SPSS软件中的相关性分析对雅砻江洪峰、鲜水河洪峰和庆大河洪峰进行相关性分析。经计算，雅砻江洪峰、鲜水河洪峰之间的相关系数为0.938，雅砻江洪峰、庆大河洪峰之间的相关系数为0.946，鲜水河洪峰、庆大河洪峰之间的相关系数为0.937，t检验的显著性概率均为0 < 0.01。由此可以判断，三者之间存在着较强的正相关性。由于GH Copula和Clayton Copula适用于变量之间存在正相关关系，Frank Copula函数适用于变量之间存在正或负相关关系。因此，本文选用GH Copula、Clayton Copula和Frank Copula函数进行组合洪水的研究。

4.3. Copula函数参数估计

参数估计是Copula函数的核心，本文采用极大似然法对Copula函数进行参数估计。双参数三维Copula函数的极大似然估计方法如下：

1) 构建似然函数

$L\left({\theta }_1,{\theta }_2\right)={\displaystyle \underset{i=1}{\overset{n}{\prod }}c\left(\left(u_{i1},u_{i2},u_{i3}\right);{\theta }_1,{\theta }_2\right)}$ (13)

$c\left(\left(u_1,u_2,u_3\right);{\theta }_1,{\theta }_2\right)=\frac{{\partial }^{3}C\left(\left(u_1,u_2,u_3\right);{\theta }_1,{\theta }_2\right)}{\partial u_1\partial u_2\partial u_3}$ (14)

$\ln \left(L\left({\theta }_1,{\theta }_2\right)\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\ln \left(c\left(\left(u_1,u_2,u_3\right);{\theta }_1,{\theta }_2\right)\right)}$ (15)

2) 似然函数最大值求解

$\frac{{\partial }^2\ln \left(L\left({\theta }_1,{\theta }_2\right)\right)}{\partial {\theta }_1\partial {\theta }_2}=0$ (16)

4.4. Copula函数拟合优度评价

拟合优度评价主要包括Genest-Rivest方法(图形评价法)、AIC信息准则法和OLS法(离差平方和最小准则)。

Genest-Rivest方法是将经验频率和理论频率进行绘图，根据点距的分布情况，来判断Copula函数对序列的拟合情况。

经验频率公式如下：

$F_{emp}\left(x_{i1},x_{i2},x_{i3}\right)=p\left(X_1\le x_{i1},X_2\le x_{i2},X_3\le x_{i3}\right)=\left({\displaystyle \underset{g=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{n}{\sum }}{\displaystyle \underset{l=1}{\overset{n}{\sum }}{n}_{gkl}-0.44}}}\right)/\left(n+0.12\right)$ (17)

式中， $F_{emp}$ 为经验频率， $n_{gkl}$ 为同时满足 $X_1\le x_{i1},X_2\le x_{i2},X_3\le x_{i3}$ 的观测值个数，n为系列长度。

AIC信息准则是衡量模型拟合优劣的一种标准，AIC值越小，说明相应的Copula函数拟合越优，其表现形式如下：

$\text{MSE}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(F_{emp}\left(x_{i1},x_{i2},x_{i3}\right)-C\left(u_{i1},u_{i2},u_{i3}\right)\right)}^2}$ (18)

$\text{AIC}=n\ln \left(\text{MSE}\right)+2k$ (19)

式中， $C\left(u_{i1},u_{i2},u_{i3}\right)$ 为理论频率值，k为模型参数的个数。

离差平方和最小准则的判断标准和AIC准则一样，OLS值越小，拟合越优。

$\text{OLS}=\sqrt{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(F_{emp}\left(x_{i1},x_{i2},x_{i3}\right)-C\left(u_{i1},u_{i2},u_{i3}\right)\right)}^2}}$ (20)

由图形评价法即可得到不同Copula函数对洪峰序列的拟合效果，如图3所示。

从图中可以看出，双参数GH Copula、Clayton Copula、Frank Copula函数所计算出的理论频率与经验频率均落在45˚线附近，且相差不大，说明三种函数对径流序列的拟合效果均较好，通过图形评价法不能直观的判断哪种函数对径流序列拟合最优，所以还需AIC信息准则和OLS法进行拟合优度评价。

计算得出不同Copula函数的AIC值和OLS值，根据AIC信息准则和OLS对拟合优度进行评价，AIC值、OLS值以及Copula函数参数θ₁、θ₂见表2。

由表2可知，GH Copula函数所计算的AIC和OLS值最小，Frank Copula函数和Clayton Copula函数所计算的AIC和OLS值均大于GH Copula函数所计算的AIC和OLS值。因此GH Copula函数对洪峰序列的拟合效果最优。所以根据AIC、OLS最小准则，选取双参数三维GH Copula函数进行组合洪水的计算。

4.5. 组合洪水的计算

根据GH Copula函数，可以得到洪峰联合分布(图4)，其中，x轴、y轴、z轴分别表示雅砻江洪峰、鲜水河洪峰、庆大河洪峰，组合洪水频率用颜色表示。由图4可以看出任意洪峰组合发生的频率及其变化趋势。根据公式(6)~(10)可以得到下游组合洪水洪峰w大于下游防洪安全泄量w₀的概率，设置不同的w₀，即可得到组合洪水分布(表3)及其分布曲线(图5)。

由组合洪水分布曲线可以看出，随着w₀的增加，组合洪水洪峰大于w₀的概率不断减小，在w₀小于3000 m³/s之内，组合洪水洪峰大于w₀的概率的减小与w₀的增加基本上呈线性关系，当w₀大于3000 m³/s后，组合洪水洪峰大于w₀的概率随着w₀的增加而减小的越来越缓慢，最终趋向于0。下游防洪控制点的安全泄量为4800 m³/s时，由表3可知，防洪风险率为1.43%。结果表明该地区的防洪安全不易受区间洪水的影响。

5. 结论

本文基于Copula函数和频率组合法，推导出n维组合洪水分布函数表达式，并构造出洪水遭遇风险估计模型。以雅砻江洪峰、鲜水河洪峰和庆大河洪峰为例，在分析其相关性基础上，采用双参数三维GH Copula、Clayton Copula、Frank Copula函数对其洪峰联合分布进行计算，通过对比分析，GH Copula函数对变量的拟合效果最优。建立了基于双参数三维GH Copula函数的三河流洪峰联合分布，然后计算得到组合洪水分布曲线，通过组合洪水分布曲线可以得到下游防洪控制点的防洪风险率以及各组合洪水发生的概率。通过实例分析证明本文所构造的模型具有可行性，且推导出的n维组合洪水分布函数对于易发生洪水遭遇地区的防洪安全提供了一定的参考。

基金项目

国家自然科学基金(51709105)；中央高校基本科研业务费专项资金资助(2019MS031)。

参考文献