较大流域下游洪水通常是由其上游不同来源洪水过程遭遇叠加而成，上下游洪水变量通常具有较强的因果性关联性^[1-2]，例如，三峡水库入库洪水往往由长江上游与水库区间流域的洪水叠加形成^[3]。一般来说，上游不同来源洪水变量，由于其产生的空间接近，气候条件与下垫面条件也较为相似，往往存在一定的相关性，能够反映洪水空间组成规律^[4-5]。研究流域不同来源洪水相关性及其下游洪水统计特征的影响有助于科学认识洪水空间组成规律对流域洪水风险的作用机制，对水利水电工程、防洪系统工程设计具有重要的参考价值。

目前关于流域上游不同来源洪水对下游洪水频率分布影响的研究主要集中于设计洪水地区组成^[6-9]，相关研究通常基于水量平衡原理，假设下游设计断面的洪水变量(通常为洪量)由上游不同来源的洪水变量直接线性叠加而成，因此可以通过不同来源洪水变量和的分布推求设计断面洪量的频率分布。宋松柏等^[10]运用Copula函数和数学积分变换原理，假设清江隔河岩水库洪量等于上游水布垭水库洪量加上水-隔区间洪量，严格地推导了二维相依随机变量和概率分布计算公式，并用于研究洪水地区组成。Merz等^[11]基于线性叠加的方法由河流干支流交汇点上游站点的洪峰流量频率分布参数推求了下游洪峰流量频率分布参数。上述研究中，下游断面洪水变量一般假设严格等于上游洪水变量的加和，然而实际中考虑到洪水的传播过程，上下游的洪水变量通常不能完全满足该条件，尤其是洪峰流量，往往不能用简单的线性叠加刻画上下游之间关系，现有基于水量平衡的方法适用性受到制约。

广义回归模型可以假设自变量与因变量服从任意概率分布，能够较为灵活地刻画不同随机变量之间的关系，在水文分析中得到广泛应用^[12-15]。Jiang等^[1]采用广义回归模型刻画上下游洪水变量的关系，基于层次模型由上游洪水频率分布推求了下游洪水频率分布，但该研究并未考虑上游干支洪水变量相关性对下游洪水频率分布的影响。干支流洪水相关性反映了流域不同来源洪水的遭遇规律，然而，现有的研究较少深入分析上游不同来源洪水变量空间相关性对下游洪水频率分布特征的影响。为此，本研究将基于广义回归模型构建上下游洪水变量之间的条件分布，然后采用Copula函数构建河流断面上游干支流洪水变量的联合概率分布，以此刻画流域不同来源洪水的相关性，最后根据全概率公式推导出下游洪水频率分布，在此基础上，将进一步探讨上游干支流洪水变量相关性强度与尾部相关性对下游洪水统计特征及其设计洪水的影响。

1 研究方法 1.1 基于层次模型的洪水频率分布估计

假设某水文站位于一条河流干支流交汇点的下游位置，该站点洪水变量定义为Y，交汇点上游两水文站的洪水变量分别定义为X₁和X₂，具体如图 1所示。

为建立上下游洪水频率分布之间关联，本研究基于层次模型由上游干支流洪水的联合概率分布推导下游洪水频率分布^[1]。洪水变量Y与上游洪水变量的关系可以由条件分布表示：

$ Y \mid X_1, X_2 \sim f_{Y \mid X_1, X_2}\left(y \mid x_1, x_2 ; \theta_{Y \mid X_1, X_2}\right) $

(1)

式中，x₁、x₂和y分别为X₁、X₂和Y的实际观测值，θ_{Y|X1, X2}为条件分布f_{Y|X1, X2}的参数向量，分别包含位置参数μ_{Y|X1, X2}、尺度参数σ_{Y|X1, X2}以及形状参数ν_{Y|X1, X2}(若分布包含该参数)。本研究采用广义回归模型构建条件分布f_{Y|X1, X2}参数与上游洪水变量的函数关系，以位置参数μ_{Y|X1, X2}为例，其表达如下：

$ \mu_{Y \mid X_1, X_2}=a_0+a_1 x_1+a_2 x_2 $

(2)

式中，a₀、a₁和a₂分别为广义回归模型的参数，采用极大似然法估计得到^[12]。采用水文频率分析中常用两参数正态分布、对数正态分布、伽马分布、韦伯分布以及三参数皮尔逊Ⅲ型(PⅢ)分布和广义极值(GEV)分布作为条件分布的备选理论分布^[1]，具体分布表达式如附表Ⅰ所示。采用赤池信息准则(AIC)^[12]优选出最优条件分布类型以及参数回归方程。

根据全概率公式，洪水变量Y的频率分布可由下式推导得到：

$ f_Y(y)=\iint f_{Y \mid X_1, X_2}\left(y \mid x_1, x_2 ; \theta_{Y \mid X_1, X_2}\right) h\left(x_1, x_2\right) \mathrm{d} x_1 \mathrm{~d} x_2 $

(3)

式中，h(x₁, x₂)表示洪水变量X₁和X₂的联合概率密度函数。根据Casella和Berger^[16]的证明，洪水变量Y的期望公式为：

$ \mathrm{E}(Y)=\mathrm{E}\left(\mathrm{E}\left(Y \mid X_1, X_2\right)\right) $

(4)

Y的方差表达式为：

$ \operatorname{var}(Y)=\mathrm{E}\left(\operatorname{var}\left(Y \mid X_1, X_2\right)\right)+\operatorname{var}\left(\mathrm{E}\left(Y \mid X_1, X_2\right)\right) $

(5)

1.2 基于Copula函数的两变量洪水联合概率分布构建

Copula函数可以比较灵活地描述多个水文变量的相关性，在水文领域具有广泛应用^[17-20]，因此本研究采用Copula函数构建X₁和X₂的联合概率分布H(x₁, x₂)，如下表示：

$ H\left(x_1, x_2\right)=C\left[F_{X_1}\left(x_1\right), F_{X_2}\left(x_2\right)\right]=C\left(u_1, u_2\right) $

(6)

式中，C(·)表示一个二维Copula函数，F_X1(·)和F_X2(·)分别表示X₁和X₂的边缘概率函数，u₁和u₂分别为x₁和x₂对应的边缘概率，即u₁=F_X1(x₁)和u₂=F_X2(x₂)，均服从区间为[0, 1]的均匀分布。

根据公式(6)，X₁和X₂联合概率密度函数如下表示：

$ h\left(x_1, x_2\right)=c\left(u_1, u_2\right) \cdot f_{X_1}\left(x_1\right) \cdot f_{X_2}\left(x_2\right) $

(7)

式中，c(·)表示Copula函数的概率密度函数，f_X1(·)和f_X2(·)分别表示X₁和X₂的概率密度函数。采用我国设计规范^[21]推荐的PⅢ型分布作为洪水变量X₁和X₂的频率分布类型(具体表达式见附表Ⅰ)，参数由极大似然法^[12]估计得到。

选取Archimedean Copula函数族中的Gumbel Copula、Frank Copula和Clayton Copula构建两变量洪水联合概率分布^[17]。为避免边缘分布计算对Copula函数参数θ估计结果的影响，本研究通过计算X₁和X₂的Kendall相关系数τ，由τ与θ的函数关系(附表Ⅱ)直接估计参数θ。基于理论频率值与实测经验频率值的均方根误差RMSE定量评估Copula函数的拟合效果，并选择拟合效果最佳的Copula函数构建两变量洪水联合概率分布。

1.3 模型求解与分布拟合优度检验

理论上，洪水变量Y的频率分布可以通过公式(3)求得，但该公式包含了PⅢ分布以及Copula函数的概率密度函数，难以求得解析解，因此本研究采用蒙特卡洛数值积分方法求解Y的推导频率分布^[1]。首先基于Copula函数生成N组边缘概率的随机样本(u_{1, i}^s, u_{2, i}^s)(i=1, 2, …, N)，并分别代入X₁和X₂分布的逆函数F_X1^-1(·)和F_X2^-1(·)求得N组两变量洪水随机样本(x_{1, i}^s, x_{2, i}^s)(i=1, 2, …, N)；然后基于广义回归模型结果，将(x_{1, i}^s, x_{2, i}^s)代入洪水变量Y条件分布参数的回归方程(公式(2))，最后根据每一组条件分布参数分别生成一个关于洪水变量Y的随机样本y_i^s(i=1, 2, …, N)。对于任意洪水变量观测值y，其累积概率可以通过经验概率公式由随机样本y_i^s(i=1, 2, …, N)推求得到：

$ F_Y(y)=\frac{1}{N+2} \sum\limits_{i=1}^N 1\left(y_i^s \leqslant y\right) $

(8)

洪水频率分布拟合优度采用Kolmogorov-Smirnov(KS)方法进行检验^[22]，如果检验统计量对应的P值大于某一阈值，则表明洪水频率分布可以通过拟合优度检验。本研究中KS检验P值的阈值取0.05。

2 研究区域与研究数据

本文选取抚河流域作为研究区域，抚河位于江西省东部，为全省五大河流之一，也是鄱阳湖主要支流，其主河道长348 km，流域面积16493 km²。抚河下游控制站为李家渡水文站，上游集水面积为15812 km²，占全流域面积的95%以上。李家渡上游干支流交汇点以上分别设有廖家湾水文站和娄家村水文站，控制流域面积分别为8723、4969 km²。抚河流域及其水文站点的空间分布如图 2所示。本研究收集了以上3座水文站1953—2013年间的年最大洪峰流量序列作为研究数据，其中李家渡站、娄家村站和廖家湾站的洪水变量分别定义为Y_d、X_c以及X_w。

3 结果分析与讨论 3.1 上下游洪水变量条件分布构建

基于广义回归模型构建李家渡洪水序列Y_d对上游娄家村站和廖家湾站洪水序列(X_w, X_c)条件分布，经过模型比选，伽马分布的AIC值最小，是拟合效果最好的理论分布(附表Ⅲ)。条件分布的位置参数(与期望值相关)与上游干支流洪水变量存在正相关关系，而条件分布的尺度参数为常数。如图 3所示，所选择的条件分布对洪水变量实测值具有较高的拟合程度，95%的分位数区间较窄，这表明李家渡站洪水主要受上游干支流的洪水叠加控制。

3.2 廖家湾站和娄家村站两变量洪水频率分布构建

采用PⅢ分布分别拟合廖家湾站和娄家村站洪水变量X_w与X_c的边缘分布，参数的极大似然估计结果如表 1所示。结果显示KS检验的P值都超过了0.05，表明PⅢ分布通过了显著性为0.05的拟合优度检验，对X_w与X_c都有很好的拟合效果。由于位于干流的廖家湾站集水面积(8723 km²)要远大于娄家村站的集水面积(4969 km²)，因此X_w的均值(2822 m³/s)明显高于X_c的均值(1905 m³/s)。图 4给出了PⅢ分布频率曲线对实测洪水点据的拟合结果，总体来说PⅢ分布对两组洪水序列的拟合较好。

在边缘分布估计的基础上，采用Copula函数构建两变量洪水(X_w, X_c)的联合概率分布，由Kendall相关系数计算Copula函数参数θ。Gumbel、Frank和Clayton三种Copula函数对应的RMSE值分别为0.044、0.038和0.049，很显然Frank Copula对实测序列的拟合效果最好，因此选用该Copula构建(X_w, X_c)的联合概率分布。图 5对比了该Copula函数计算的理论联合概率与经验联合概率，结果表明两者具有较好的吻合程度。

3.3 李家渡站洪水频率分布推导

根据洪水变量(X_w, X_c)的联合概率分布生成100000组随机样本，然后进一步由(X_w, X_c)的随机样本生成100000个李家渡洪水变量Y_d的随机样本，并以此来刻画Y_d的频率分布。结果显示，推导洪水频率分布的KS检验P值为0.521，能够通过显著性水平为0.05的拟合优度检验，因此由层次模型推求的洪水频率分布具有较好的拟合效果。

为与理论分布方法相对比，本研究还采用PⅢ分布拟合了李家渡洪水序列的频率分布，结果如表 2所示，总体来说2种方法得到的洪水频率分布统计参数基本一致。图 6对比了PⅢ分布推求的理论洪水频率分布曲线与层次模型推导得到的洪水频率分布曲线，总体来说两者差别不大。需要指出的是蒙特卡洛数值积分方法为一种近似的求解方法，存在一定的计算误差，这导致层次模型推导得到洪水频率分布曲线上尾部分不够平滑。

3.4 上游干支流洪水相关性对下游洪水频率分布特征影响

如公式(3)所示，下游洪水变量的频率分布取决于上下游洪水变量之间的条件分布与上游洪水变量的联合分布。根据前文结果，下游李家渡站洪水变量条件分布为伽马分布，其位置参数μ_{Yd|Xw, Xc}即为条件分布期望值E(Y_d|X_w, X_c)。李家渡洪水变量均值E(Y_d)有如下理论表达式：

$ \mathrm{E}\left(Y_{\mathrm{d}}\right)=\mathrm{E}\left(\mathrm{E}\left(Y_{\mathrm{d}} \mid X_{\mathrm{w}}, X_{\mathrm{c}}\right)\right)=\mathrm{E}\left(a_0+a_1 X_{\mathrm{w}}+a_2 X_{\mathrm{c}}\right)=a_0+a_1 \mathrm{E}\left(X_{\mathrm{w}}\right)+a_2 \mathrm{E}\left(X_{\mathrm{c}}\right) $

(9)

式中, E(X_w)和E(X_c)分别为廖家湾站与娄家村站洪水变量的均值。李家渡洪水变量的方差var(Y_d)有如下理论表达式：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\; \operatorname{var}\left(Y_{\mathrm{d}}\right)=\mathrm{E}\left(\operatorname{var}\left(Y_{\mathrm{d}} \mid X_{\mathrm{w}}, X_{\mathrm{c}}\right)\right)+\operatorname{var}\left(\mathrm{E}\left(Y_{\mathrm{d}} \mid X_{\mathrm{w}}, X_{\mathrm{c}}\right)\right)= \\ \mathrm{E}\left({\sigma_{Y_{\mathrm{d}} \mid X_{\mathrm{w}}, X_{\mathrm{c}}{ }}}^2 \cdot\left[a_0+a_1 X_{\mathrm{w}}+a_2 X_{\mathrm{c}}\right]^2\right)+\operatorname{var}\left(a_0+a_1 X_{\mathrm{w}}+a_2 X_{\mathrm{c}}\right)= \\ {\sigma _{{Y_{\rm{d}}}\mid {X_{\rm{w}}},{X_{\rm{c}}}}}^2 \cdot \left\{ \begin{array}{l} {\rm{E}}\left( {a_0^2} \right) + a_1^2{\rm{E}}\left( {X_{\rm{w}}^2} \right) + a_2^2{\rm{E}}\left( {X_{\rm{c}}^2} \right) + 2{a_0}{a_1}{\rm{E}}\left( {{X_{\rm{w}}}} \right) + 2{a_0}{a_2}{\rm{E}}\left( {{X_{\rm{c}}}} \right)\\ + 2{a_1}{a_2}\left[ {{\rm{E}}\left( {{X_{\rm{w}}}} \right) \cdot {\rm{E}}\left( {{X_{\rm{c}}}} \right) + {\mathop{\rm cov}} \left( {{X_{\rm{w}}},{X_{\rm{c}}}} \right)} \right] + \end{array} \right\} + \\ a_1^2{\mathop{\rm var}} \left( {{X_{\rm{w}}}} \right) + a_2^2{\mathop{\rm var}} \left( {{X_{\rm{c}}}} \right) + 2{a_1}{a_2}{\mathop{\rm cov}} \left( {{X_{\rm{w}}},{X_{\rm{c}}}} \right) \end{array} $

(10)

式中，var(X_w)和var(X_c)分别为廖家湾站与娄家村站洪水变量的方差，cov(X_w, X_c)为上述两个洪水变量的协方差。从公式(9)和(10)可以看出Y_d的均值仅取决于上游两变量洪水X_w和X_c的均值，Y_d的方差(或C_v)除与廖家湾站与娄家村站洪水变量各自统计特征值有关外，还受两者协方差cov(X_w, X_c)的控制，即X_w和X_c之间的相关性。

为进一步探讨上游干支流洪水变量相关性对下游李家渡站洪水变量统计特征及其洪水设计值的影响，本研究基于Frank Copula函数参数θ与Kendall相关系数τ的函数关系(附表Ⅱ)，生成τ不同取值条件下Y_d的频率分布。图 7给出Y_d统计特征值随(X_w, X_c)相关系数τ的变化，由于蒙特卡洛数值积分存在一定误差，Y_d统计特征值存在一定的随机波动，但可以明显看出Y_d分布的均值基本不随τ的变化而变化，而C_v值随着τ的增大呈现非常明显的增加趋势，上述发现与公式(9)和(10)的结果一致。如图 7b所示，当(X_w, X_c)的Kendall相关系数τ为0时，Y_d分布的C_v值约于0.385，而当τ为0.9时, C_v值增加至0.520。(X_w, X_c)相关性对C_s值的影响相对较小，当τ小于0.5时，C_s值基本没有变化，而当τ超过0.5时，C_s随相关性的增强而增大。

除了相关性强度，进一步探讨上游洪水变量(X_w, X_c)尾部相关性对下游洪水频率分布的影响。不同Copula函数对尾部特征的刻画差别较大^[17]，例如Gumbel Copula侧重上尾相关性，Clayton Copula侧重下尾相关性，而Frank Copula的上、下尾相关系数相等且都为零，因此将分别基于上述3种Copula函数生成廖家湾站与娄家村站洪水变量的随机样本，在此基础上分别推求李家渡站洪水变量Y_d对应的频率分布。不同Copula函数对应的洪水频率分布统计特征值对比如图 7所示，可见(X_w, X_c)尾部相关性对李家渡站洪水频率分布的均值没有影响，但对高阶参数C_v和C_s影响较为显著，尤其是相关系数τ介于0.1~0.8之间时。从图 7b可以看出，Gumbel Copula函数侧重上尾相关性，由该函数生成的李家渡站洪水频率分布C_v值最大，其次为无尾部相关性的Frank Copula，侧重下尾相关的Clayton Copula函数得到C_v值最小，偏态系数C_s也呈现类似的对比结果(图 7c)。同时还可以看出，尾部相关性对C_s的总体影响程度甚至要明显超过相关性强度。同时，还可以发现在上游洪水变量相关性极弱或极强的情况下，尾部相关性的差异对下游洪水频率分布的C_v和C_s的影响都不再明显。上述发现表明流域上游不同来源洪水发生的空间相关性强度与尾部相关性都是影响下游洪水频率分布规律的重要因素。

图 8给出了上游廖家湾和娄家村站洪水变量(X_w, X_c)相关性对不同重现期下李家渡站设计洪峰流量的影响，从中可以看出，上游洪水变量相关性强度变化对下游不同重现期设计洪水的影响有所不同。随着上游洪水变量(X_w, X_c)相关性增强，下游李家渡站2年一遇的设计洪水值呈现微弱的减小趋势，τ=0.9对应的设计值较τ=0的设计值减小了约240 m³/s，而对于5年及以上重现期，李家渡站水文设计值随(X_w, X_c)相关性增强呈现增加趋势，且靠近上尾部分(重现期较长)的极端水文设计值对(X_w, X_c)相关性变化更为敏感，随着(X_w, X_c)的相关系数τ由0增加至0.9，5、10、20、50、100年一遇的设计洪水增加幅度分别为400、900、1400、2000、2300 m³/s，相对增加幅度分别为7%、13%、18%、22%、24%。上游洪水变量(X_w, X_c)相关性越强意味着上游洪水的发生越趋向于同频率，上述结果进一步表明，流域内不同来源洪水的发生规律越趋向于同频率，流域的极端洪水风险就越高。

如图 8所示，随着重现期的变长，上游洪水变量(X_w, X_c)尾部相关性对李家渡站水文设计值的影响越来越显著。当设计值重现期超过50年后，侧重上尾相关的Gumbel Copula函数对应的水文设计值要显著大于无尾部相关性Frank Copula函数以及侧重下尾相关性Clayton Copula函数对应的水文设计值。上游洪水变量上尾相关性越强意味着同时发生极端洪水事件的可能性越高，相应的下游洪水极端洪水风险也就越高。因此，在极端水文事件分析中，尾部相关性对水文风险的刻画与计算具有不可忽略的作用。从图 8同样可以看出，在上游洪水变量相关性极弱或极强的情况下，尾部相关性的差异对下游水文设计值的影响不再显著。

4 结论与讨论

为研究上游不同来源洪水相关性对下游洪水频率分布特征的影响，选取鄱阳湖支流抚河流域作为研究区，基于层次模型由上游干支流洪水变量的联合概率分布推求了下游洪水频率分布，并探讨了上游洪水变量相关性强度以及尾部相关性对下游洪水统计特征的影响。主要研究结论为：由层次模型推导得到的抚河下游李家渡站洪水频率分布有较好的拟合效果，与传统的PⅢ分布拟合的频率曲线基本相当；上游不同来源洪水相关性是影响下游洪水频率分布的重要因素，理论公式推导与数值实验结果都表明，上游洪水变量相关性强度与尾部相关性对下游洪水频率分布均值没有影响，但均对高阶统计参数有明显的影响，下游洪水频率分布的C_v值会随上游洪水变量相关性的增强而显著增大，上尾相关性越强C_v和C_s越大，对应的流域洪水极端洪水风险也就越高；流域极端洪水设计值对流域不同来源洪水相关性强度及尾部相关性的变化更加敏感。

现有的设计洪水地区组成研究方法假定下游断面洪量为上游断面与区间洪量的加和，需要指出的是该假设仅适用于相应洪量情景，忽略了洪水的传播过程，这与实际情况并不完全相符，且无法用于洪峰的地区组成。本文采用的广义回归模型可以灵活描述上下游洪水变量的线性与非线性关系，为解决上述问题提供了一种可行的途径。此外，上游洪水变量作为下游洪水变量具有物理成因机制的解释变量，在上游的洪水序列的观测期长于下游的情况下，引入上游更长洪水样本信息可有效提高下游洪水频率分布估计成果的稳健性^[5]，本文采用的层次模型可以为上述思路的实现提供方法支撑^[23]。

本文提出的基于层次模型的洪水频率分布推求方法依赖于上下游洪水变量之间的条件分布、Copula函数、上游洪水频率分布的估计结果，另一方面，采用的蒙特卡洛数值积分为一种近似的求解方法，也存在一定的计算误差，上述因素对洪水频率分布估计成果的不确定性有待进一步深入研究。

5 附录

附表Ⅰ~Ⅲ见电子版(DOI: 10.18307/2024.0453)。