全球气候变化加大了降水径流等气象水文要素的空间分布差异，影响着多江流域水资源及气象水文极值的组合特征^[1].丰枯遭遇实质是气象水文变量的多区域组合问题，涉及到降水的分区组合^[2-4]、流域干流和其支流的径流组合^[5]、多条河流的径流组合^[6]或极值流量组合^[7]、调水工程水源区和受水区的降水或径流组合^[8-9]、多个水库的来水量组合^[10]等.

通常认为丰枯遭遇的各区域气象水文变量是个随机过程，但它们之间又存在相依性.多维随机变量构造函数Copula在反映单个随机变量分布特征的同时又能呈现变量之间的相依性，现已被广泛应用于洪水特性分析、干旱特性分析、水文过程模拟、气象水文变量多维联合分布模拟等^[11-21].其中，基于Copula函数的径流丰枯遭遇研究^[5]，是通过构建多条河流的径流Copula模型，分析径流丰枯组合的联合概率特征.径流丰枯遭遇联合分布模拟多为二维模型和三维模型^[6]，冯平^[22]等运用四元Copula函数研究了南水北调西线工程水源区河流与受水区黄河上游的丰枯遭遇问题，同一流域三变量以上的丰枯遭遇研究较少.同一流域涉及多条河流来水时，变量的维数变高，问题变复杂，Copula函数的应用使丰枯遭遇的问题简单化.

鄱阳湖为中国最大的淡水湖，接纳饶河、抚河、信江、赣江和修水的五江来水，经过湖盆调蓄后汇入长江^[23].作为国际重要湿地和长江干流重要调蓄性湖泊，鄱阳湖水量的丰枯变化影响着湖区生物多样性保护及长江中下游的洪水调蓄等^[24].前人对鄱阳湖水量的研究多集中在流域单一变量的水文特征分析方面，如径流变化^[25-27]、枯水特征^[28]、水文干旱特征^[23]和流域旱涝规律^[28]等.近些年，有学者也开始利用Copula函数研究鄱阳湖极值流量的遭遇特征和枯水多特征的联合分布^[29]等，但均为两变量联合分布模拟.受全球气候变化的影响，鄱阳湖流域水资源量发生变化，流域内出现极端干旱事件，旱涝灾害频繁^[30].研究表明，5大水系来水变化与湖区干旱等水文极值事件有密切关系^[31].因此，研究5大水系来水丰枯规律，能加强流域内的水量调控，对区域防洪抗旱有重要意义.为加强流域内水量调控，需分析流域各条河流来水量之间的相互关系，传统方法是建立变量之间的相关关系或选取若干典型年^[32].为保证调水方案的合理性，需要分析各水系水量的丰枯遭遇概率，传统方法无法回答遭遇概率大小的问题，Copula函数的引入能很好地解决此问题.对于多条河流来水的鄱阳湖水系来说，多维分布可对鄱阳湖流域局部至流域整体的径流规律进行研究.

因此，本文基于鄱阳湖5大水系径流，运用Copula函数构建二维至五维的联合分布模型，探讨多条河流来水丰枯遭遇的多维模拟方法，从多个维度上分析鄱阳湖水系径流丰枯遭遇特征变化，为流域水资源管理及水旱灾害预防提供依据.

1 方法和数据 1.1 径流多维联合分布

某一河流在一定时间尺度上的径流，如年径流、汛期径流和非汛期径流，其逐年变化一般认为是个随机过程.一个流域面积较大的河湖水系(如鄱阳湖水系)，其水量变化同时受到多条河流的影响.受地理位置及气候条件影响，同属一个水系的多条河流之间，其径流过程具有一定的相依性.因此，同一流域多条河流径流丰枯遭遇分析的前提是构建多维联合分布模型^[6].现已广泛应用于多维联合分布的Copula函数，在水文多变量联合分布中模拟较好^[33]，可用于径流丰枯遭遇的多维联合分布.

有d条河流的径流系列(X₁, …, X_d)，其边缘分布函数分别为F(x₁), …, F(x_d).根据Sklar定理，d维联合分布函数H(x₁), …, F(x_d)可定义为：

$ H\left( {{x_1}, \cdots , {x_d}} \right) = C\left[ {F\left( {{x_1}} \right), \cdots , F\left( {{x_d}} \right)} \right] = C\left( {{u_1}, \cdots , {u_d}} \right) $

(1)

式中，C(u₁), …, F(u_d)为边缘分布在[0, 1]区间的Copula联结函数；u₁=F(x₁), …, u_d=F(x_d)为各单变量的边缘分布，即累积频率；d为联合分布的维数，即变量个数.

Copula函数构造多样，其中Meta-Gaussian Copula在降水径流的多变量联合分布模拟中具有良好效果^[4]. d维Meta-Gaussian Copula的分布函数C(u₁), …, F(u_d)和密度函数c(u₁), …, F(u_d)分别为：

$ C\left( {{u_1}, \cdots , {u_d}} \right) = \int_{ - \infty }^{{b_1}} { \cdots \int_{ - \infty }^{{b_d}} {g\left( {{\omega _1}, \cdots , {\omega _d}} \right){\rm{d}}{\omega _1} \cdots {\rm{d}}{\omega _d}} } $

(2)

$ c\left( {{u_1}, \cdots , {u_d}} \right) = {\left| \Sigma \right|^{ - 1/2}}\exp \left( { - {\zeta ^T}{\Sigma ^{ - 1}}\zeta /2 + {\zeta ^T}\zeta /2} \right) $

(3)

其中：

$ g\left( {{\omega _1}, \cdots , {\omega _d}} \right) = {\left( {2{\rm{ \mathsf{ π} }}} \right)^{ - d/2}}{\left| \Sigma \right|^{ - 1/2}}\exp \left( { - {\omega ^T}{\Sigma ^{ - 1}}\omega /2} \right) $

(4)

式中，b₁=Φ^-1(u₁), …, b_d=Φ^-1(u_d)；ω为被积函数变量矩阵ω=[ω₁, …, ω_d]^T；ζ=[b₁, …, b_d]^T; Σ为相关系数矩阵，即：

$ \Sigma = \left[ {\begin{array}{*{20}{c}} 1& \cdots &{{\rho _{ij}}}\\ \vdots&\ddots&\vdots \\ {{\rho _{ij}}}& \cdots &1 \end{array}} \right], {\rho _{ij}} = \left\{ {\begin{array}{*{20}{c}} {1, i = j}\\ {{\rho _{ij}}, i \ne j} \end{array}} \right. $

(5)

式中，ρ_ij∈[-1, 1];i=1, …, d; j=1, …, d, 为任意两个变量之间的相关系数.

各径流系列的边缘分布可采用水文统计分析较为常用模型进行拟合，即伽马分布(GAMMA)、对数正态分布(LOGN)、韦布尔分布(WBL)、广义极值分布(GEV)、广义Logistic(GLO)、广义正态(GNO)、皮尔逊Ⅲ型(P-Ⅲ).采用Pearson-γ计算河流之间的径流相关性.多个维度的联合分布及其边缘分布模型拟合，分别通过基于Rosenblatt变换的Cramervon Mises检验之一的S_n^(B)[34]和Anderson-Darling(A-D)法的A_n²进行检验.边缘分布模型综合考虑均方根误差(RMSE)、赤池信息准则(AIC)和统计量的P值，进行优选^[35].联合分布及边缘分布模型参数均采用最大似然法估计.经验分布频率采用次序统计累积频率并进行适当的变换^[36-37].

1.2 径流丰枯分类

在水资源分析中，基于降水或径流的丰枯分类并没有统一的划分方法和标准.一种方法是基于距平百分率，划分丰、平、枯的标准^[38].更为常用的方法是根据径流频率，如以频率37.5 %和62.5 %划分丰、平、枯3种水平年^[2]，也有的采用25 %和75 %作为分界线划分丰、平、枯^[5].

在《水文基本术语和符号标准》(GB/T 50095—2014)设计年径流中，提出了特枯水年(年降水量或年河川径流量为历年最小值或接近历年最小值的年份)、枯水年(年降水量或年河川径流量显著小于多年均值的年份)、平水年(年降水量或年河川径流量接近多年均值的年份)、丰水年(年降水量或年河川径流量显著大于多年均值的年份).这说明在年径流设计实践中，有时需要更为精细分类，以关注特枯水年或特丰水年的水资源状况.

综上分析，本文拟采用从小到大的累积频率12.5 %、37.5 %、62.5 %和87.5 %为分界值，将径流系列分为特枯水年P≤12.5 %、偏枯水年12.5 % < P≤37.5 %、平水年37.5 % < P≤62.5 %、偏丰水年62.5 % < P≤87.5 %和特丰水年P>87.5 %.

1.3 径流丰枯遭遇

单一河流的径流丰枯划分，呈现的是径流的概率区间.丰枯遭遇是指多条河流的径流在不同概率区间的组合事件.若d条河流各自的径流概率区间边界为[a_i, b_i], i=1, 2…, d，a_i∈[0, 1]、b_i∈[0, 1]且a_i≤b_i，d维丰枯遭遇概率则为：

$ P\left( {{a_1} \le {u_1} \le {b_1}, \cdots , {a_i} \le {u_i} \le {b_i}, \cdots , {a_d} \le {u_d} \le {b_d}} \right) = \\ \sum\limits_{{i_1}}^2 { \cdots \sum\limits_{{i_j}}^2 { \cdots \sum\limits_{{i_d}}^2 {{{\left( { - 1} \right)}^{{i_1} + \cdots {i_j} + \cdots {i_d}}}C\left( {{u_{1, {i_1}}}, \cdots {u_{i, {i_j}}}, \cdots {u_{d, {i_d}}}} \right)} } } $

(6)

式中，i_j=1, 2；j=1, 2, …, d；u_{i, 1}=a_i；u_{i, 2}=b_i.

根据概率容斥原理^[39]，式(6)可解释为：

$ \begin{array}{l} P\left( {{a_1} \le {u_1} \le {b_1}, \cdots , {a_i} \le {u_i} \le {b_i}, \cdots , {a_d} \le {u_d} \le {b_d}} \right)\\ = C\left( {{b_1}, {b_2}, \cdots , {b_d}} \right)\\ \;\;\; - \left[ {C\left( {{a_1}, {b_2}, {b_3}, \cdots , {b_d}} \right) + C\left( {{b_1}, {b_2}, \cdots , {b_{d - 1}}, {a_d}} \right)} \right]\\ \;\;\; + \left[ {C\left( {{a_1}, {b_2}, {b_3}, \cdots , {b_d}} \right) + C\left( {{b_1}, {b_2}, \cdots , {b_{d - 2}}, {a_{d - 1}}, {a_d}} \right)} \right]\\ \;\;\;\; \vdots \\ \;\;\; + {\left( { - 1} \right)^k}\left[ {C\left( {{a_1}, {a_2}, \cdots , {a_k}, {b_{k - 1}}, \cdots , {b_d}} \right) + C\left( {{b_1}, {b_2}, \cdots , {b_{d - k}}, {a_{d - k + 1}}, \cdots , {a_d}} \right)} \right]\\ \;\;\;\; \vdots \\ \;\;\; + {\left( { - 1} \right)^d}C\left( {{a_1}, {a_2}, \cdots , {a_d}} \right) \end{array} $

(7)

式中，k=0, 1, 2, …, d，第k+1项中有$个联合概率值C相加.

依据上述的特枯、偏枯、平水、偏丰和特丰，共5个概率区间，d条河流的丰枯遭遇情况共有5^d个组合，即2、3、4和5条河流的丰枯遭遇组合数分别为25、125、625和3125个.为分析不同组合的概率变化，对多条河流不同概率区间组合的编号顺序按照每条河流从特枯至特丰依次变化.以d=2为例，丰枯组合及编号如图 1所示.因此，多条河流概率区间的组合即同特枯、同偏枯、同平、同偏丰和同特丰的编号分别是：2条河流为1、7、13、19、25，3条河流为1、32、63、94、125，4条河流为1、157、313、469、625，5条河流为1、782、1563、2344、3125.

1.4 流域概况及数据

鄱阳湖位于长江中游江西省北部，属于长江流域(图 2).鄱阳湖水系为辐射状，流域面积约为16.22万km².由于地处中国亚热带湿润季风气候区，雨量充沛，但存在明显的年内变化，通常每年4—9月为汛期，10月至次年3月为非汛期.

鄱阳湖水系五大河流的水文控制站分别为饶河的渡峰坑和虎山、信江的梅港、抚河的李家渡、赣江的外洲、修水的万家埠和虬津.选用各站1955—2001年共47年(水文年)的逐月径流资料，分析年径流、汛期径流和非汛期径流3个系列的丰枯遭遇特征.部分缺失数据通过选取邻近的水文站资料采用线性插值进行插补，饶河采用渡峰坑和虎山两站径流叠加数据.由于虬津站资料缺失，修水选用其支流潦河万家埠站的径流数据，集水面积约占修水的26.4 %.

2 结果与讨论 2.1 各条河流的径流常规统计特征

鄱阳湖五大水系的径流基本特征见表 1，将饶河、信江、抚河、赣江和修水5条河流分别编号.从径流量均值上看，由于赣江集水面积较大，径流量远大于其他河流；修水支流潦河集水面积小，其径流量远小于其他河流.信江和抚河的集水面积相差很小，但前者的径流量明显大于后者.从径流模数上看，信江最大，其次为饶河和修水，抚河和赣江的径流模数较小.从径流量的年际变差系数来看，赣江明显小于其他河流；非汛期径流变差系数大于年径流系列和汛期径流系列.以上结果表明：在单位面积上单位时间所产生的径流量和径流年际变化程度上，北部河流大于南部河流；在时间变化上，非汛期径流的年际变化程度明显大于年径流和汛期径流.

各径流系列的箱形统计如图 3.模比系数等于1的虚线为样本的均值线，箱内实线为中位值. 5条河流的年径流均值与中位值非常接近，非汛期径流均值均大于中位值.抚河的汛期径流均值略小于中位值，其他河流的汛期径流均值大于中位值. 图 3中上下框边界分别对应上四分位数Q₁和下四分位数Q₃，两端沿虚线延伸出的边界分别对应着Q₁-1.5(Q₃-Q₁)和Q₃+1.5(Q₃-Q₁)，在此区间内的值通常认为是正常值，区间外的则被认为是异常值(如图中十字符号所示).除了非汛期径流出现了异常最大值，大部分径流数据均在正常值范围内.

2.2 径流统计分布模型拟合

每条江各径流系列所有模型拟合检验P值均大于显著性检验水平0.05(表 2)，表明常用的7种单变量分布模型均可作为备选的边缘分布.基于RMSE和AIC最小原则，饶河、抚河、信江、赣江和修水的径流系列最优分布：年径流分别为GLO、GLO、GLO、GLO、GNO，汛期径流分别为GNO、P-Ⅲ、GEV、GEV、GNO，非汛期径流分别为GLO、LOGN、GAMMA、GNO、LOGN.为便于对比各条河流径流分布特征差异，同一种径流系列选用同一种模型模拟，即5条河流的年径流、汛期径流、非汛期径流的边缘分布分别选用GLO、GNO、LOGN模型.选用模型和最优模型的RMSE和AIC值差别均较小，对径流系列概率计算影响甚微.

5条河流的年径流形态参数均为负值，其绝对值较小.抚河的汛期径流形态参数为很小的正值，其他河流均为负值，其中饶河、信江和修水的汛期径流形态参数的绝对值较大(表 3).结合选用模型的概率密度曲线(图 4)，年径流趋向于对称的分布，汛期径流、非汛期径流均呈现正偏态分布(右偏)，水文现象中大多数都属于正偏.

两条河流之间的径流相关性及联合分布参数如表 4所示，其中饶河与抚河、赣江的相关系数相对较小，其次为修水与抚河、赣江.但从总体上来说，鄱阳湖水系5条河流之间呈现较高的相关性，其中北部与南部河流之间的相关性相对要低些.另外，河流之间的非汛期径流相关性高于年径流和汛期径流. Gaussian Copula模型参数反映了各变量之间的相关性，与相关系数接近.

多维径流联合分布模拟检验P值，均大于显著性检验水平0.05(表 5)，表明Gaussian Copula能够很好地模拟鄱阳湖水系多条河流之间的径流联合分布.

2.3 多维径流丰枯遭遇分析

基于各径流系列联合分布，2~5条河流多个维度的丰枯遭遇概率如图 5~8.中空圆圈标注的是同一概率区间遭遇的联合概率，即丰枯同步概率，沿横坐标方向依次表示为同特枯、同偏枯、同平、同偏丰和同特丰.由于1条河流的特枯和特丰的区间概率相等、偏枯和偏丰的区间概率相等，分别为0.125和0.25，选用的Gaussian Copula函数具有对称性，即以同平组合为中心，丰枯两侧对称.因此，多条河流组合的同特枯和同特丰的联合概率相等、同偏枯和同偏丰的联合概率相等.以此类推，如2条河流的特枯与偏枯组合的联合概率等于偏枯与特枯组合的联合概率，5条河流的特枯、偏枯、平水、偏丰与特丰组合的联合概率等于特丰、偏丰、平水、偏枯与特枯组合的联合概率.

对比不同概率区间的同步联合概率，2条河流组合的联合概率：同偏枯或同偏丰>同平>同特枯或同特丰；3条河流组合的联合概率：同偏枯或同偏丰>同特枯或同特丰>同平；4条或5条河流组合的联合概率：同特枯或同特丰>同偏枯或同偏丰>同平.也就是说，尽管对每1条河流概率区间设置时，偏枯、平水和偏丰的区间概率相对较大，均为0.25；特枯和特丰的区间概率相对较小，为0.125；但是2条或3条河流的同偏枯或同偏丰的联合概率明显大于同平的联合概率，4条或5条河流同特枯或同特丰、同偏枯或同偏丰的联合概率均明显大于同平的联合概率，且5条河流的同特枯或同特丰的联合概率也明显大于同偏枯或同偏丰的联合概率.以上结果表明，多条河流丰枯遭遇随着维数的增加，丰枯同步的最大联合概率值由同偏枯或同偏丰趋向于丰枯两端的同特枯或同特丰.另外，对于相同的概率区间，非汛期径流的丰枯同步联合概率明显大于年径流和汛期径流，而年径流和汛期径流之间的丰枯同步联合概率差别较小.

2条河流组合丰枯完全异步即1条河特丰1条河特枯概率区间对应的联合概率接近于0.1条河流特丰(特枯)另1条河流偏枯(偏丰)组合的联合概率相比丰枯同步概率相差较大.相邻2条河流组合的丰枯异步概率小于非相邻2条河流组合丰枯异步的概率.对比图 5~8，多条河流丰枯遭遇随着维数的增加，丰枯异步的组合变复杂，但遭遇概率变化趋势是一致的.

5条河流多维的丰枯同步联合概率计算结果见表 6.

从2条河流丰枯遭遇空间组合上来看，年径流丰枯同步最大联合概率为饶河与修水同偏枯或同偏丰的0.121，其次为信江与抚河的同偏枯或同偏丰的0.119；汛期径流丰枯同步最大联合概率为信江与抚河的同偏枯或同偏丰的0.126，其次为饶河与修水的0.123；非汛期径流丰枯同步最大联合概率为信江与抚河同偏枯或同偏丰的0.142，其次为抚河与赣江的0.140.年径流、汛期径流和非汛期径流的丰枯同步最小联合概率，均为饶河与赣江的同特枯或同特丰，分别为0.037、0.043和0.050；其次均为饶河与抚河的同特枯或同特丰，分别为0.041、0.045和0.058.总体上来说，饶河与信江、饶河与修水、信江与抚河、抚河与赣江之间的丰枯同步概率相对大些，而河流之间其他组合情况下的丰枯同步概率相对小些.

对于3条河流的丰枯遭遇空间组合，年径流丰枯同步最大联合概率为饶河、信江与修水的同偏枯或同偏丰的0.058，其次为信江、抚河与赣江的同偏枯或同偏丰的0.056；汛期径流丰枯同步最大联合概率为信江、抚河与赣江的同偏枯或同偏丰的0.061，其次为饶河、信江与修水的0.058；非汛期径流丰枯同步最大联合概率为信江、抚河与赣江的同偏枯或同偏丰的0.081，其次为饶河、信江与抚河的0.073.年径流和汛期径流的丰枯同步最小联合概率均为饶河、抚河与赣江的同平，分别为0.026和0.028；其次均为抚河、赣江与修水的同平，分别为0.028和0.029；非汛期径流的丰枯同步最小联合概率为饶河、抚河与修水同平的0.031，其次为饶河、赣江与修水同平的0.032.

对于4条河流的丰枯遭遇空间组合，年径流、汛期径流和非汛期径流的丰枯同步最大联合概率均为信江、抚河、赣江与修水的同特枯或同特丰，分别为0.030、0.031和0.042；丰枯同步最小联合概率均为饶河、抚河、赣江与修水的同平，分别为0.011、0.012和0.016.

对于5条河流即全流域的丰枯遭遇，年径流、汛期径流和非汛期径流丰枯同步最大联合概率分别为0.020、0.023和0.032，均为同特丰或同特枯；丰枯同步最小联合概率分别为0.006、0.007和0.011，均为同平.

以上结果表明，同处于流域北部河流(如饶河、信江与修水)或南部河流(如抚河与赣江)或相邻河流(如信江与抚河)，它们之间组合的同步联合概率相对要大些，而南北河流组合(如饶河与赣江、修水与赣江)的同步联合概率相对要小.

另外，多条河流丰枯遭遇随着维数的增加，丰枯组合数也增加，丰枯同步的联合概率明显下降.

3 结论

运用Gaussian Copula函数，构建了鄱阳湖流域5条河流径流二维至五维的联合分布.基于特枯、偏枯、平水、偏丰和特丰的径流丰枯分类，计算了年径流、汛期径流和非汛期径流的多条河流的多维丰枯遭遇联合概率，分析探讨了多维丰枯同步联合概率的变化特征.主要结论如下：

鄱阳湖水系河流之间具有较高的相关性，Gaussian Copula函数能较好地模拟二维至五维的径流联合分布.由于选用的Gaussian Copula函数具有对称性，即以同平组合为中心，丰枯两侧对称.

多条河流的丰枯遭遇随着维数的增加，丰枯同步的最大联合概率趋向于丰枯两端的概率区间.如5条河流的同特枯或同特丰的联合概率明显大于同偏枯或同偏丰的联合概率，同偏枯或同偏丰的联合概率明显大于同平的联合概率.

对于相同的概率区间，非汛期径流的丰枯同步联合概率明显大于年径流和汛期径流的同步联合概率，而年径流和汛期径流之间的丰枯同步联合概率差别较小.

同处于流域北部河流(如饶河、信江与修水)或南部河流(如抚河与赣江)或相邻河流(如信江与抚河)之间组合，其同步联合概率相较于其他组合大，而南部与北部河流组合(如饶河与赣江、修水与赣江)的同步联合概率相对较小.另外，多条河流丰枯遭遇随着维数的增加，丰枯组合数也增加，丰枯同步的联合概率明显下降.

总之，本文通过构建多维的径流联合分布模型，从多个维度分析了水系径流丰枯遭遇，主要是同步联合概率的情况及空间特征，对区域防洪抗旱以及流域内水资源协调利用有重要意义.但由于Gaussian Copula函数的对称性限制，研究结论未能反映丰枯两侧的差异性.因此，在多维Copula函数的选择方面，需要进一步研究和探讨.