一种基于Copula函数的水文模型综合不确定性分析方法

摘要

本发明公开了一种基于Copula函数的水文模型综合不确定性分析方法，包括以下步骤：⑴收集流域的数据资料；⑵建立水文模型模拟流域出口断面流量过程；⑶确定实测流量和模拟流量的边缘概率分布函数；⑷利用Copula函数构建实测流量和模拟流量的联合概率分布函数；⑸求解给定模拟流量时实测流量条件概率分布函数；⑹获取实测流量的中位数和不确定区间。本发明可以更加全面地同时考虑模型参数和模型结构的不确定性，得到水文模型的综合不确定性。本发明独立于确定性水文模型，可以和任意复杂度的确定性水文模型协同集成，而不需对模型附加任何假定，为分析水文模型的综合不确定性提供了通用性的理论框架。

说明书

技术领域

本发明涉及流域水文模型，具体地说是一种基于Copula函数的水文模 型综合不确定性分析方法。

背景技术

流域水文模型是指用模拟方法将复杂的水文现象和过程经概化所建成 的数学模型。目前，水文模型广泛应用于流域防洪抗旱、水资源规划与管 理、水环境和生态系统保护、气候变化及人类活动对水资源影响分析等领 域。实际上，水文系统是复杂、非线性的过程，采用相对简单的数学公式 来描述真实复杂的流域水文过程往往会出现失真，导致水文模型不可避免 地存在不确定性。因此，深入分析水文模拟结果的不确定性，可以为决策 者提供更加充分的风险信息，具有重要的科学意义和应用价值。

水文模型的不确定性通常包括参数不确定性和结构不确定性两个方 面。针对模型参数的不确定性问题，英国水文学家Beven和Binley于1992 年提出了基于贝叶斯理论的普适似然不确定性估计(Generalized LikelihoodUncertaintyEstimation，GLUE)方法，在国内外得到了广泛 的应用(张利茹,管仪庆,王君等.GLUE法分析水文模型参数不确定性的 研究[J].水力发电,2010,36(5):14-16)。GLUE方法虽然原理简单，易 于操作，但存在并非经典的贝叶斯方法、主观判断参数可行域阈值、推求 的参数后验概率分布不具有明显的统计学意义等问题，使得该法的可靠性 和适用性遭到质疑。

为此，有学者提出了基于马尔科夫链蒙特卡罗(MarkovChainMonte Carlo，MCMC)随机抽样技术的贝叶斯方法来研究水文模型参数的不确定性 (梁忠民,李彬权,余钟波等.基于贝叶斯理论的TOPMODEL参数不确定 性分析[J].河海大学学报(自然科学版),2009,37(2):129-132)。MCMC 方法避免了直接求解积分困难，能够推导出具有显著统计学意义的模型各 参数后验分布，但将MCMC方法应用于实际水文问题时是否有足够的计算效 率、能否满足实际要求等问题都存在较大争议，有待进一步研究。

贝叶斯模型加权平均(BayesianModelAveraging，BMA)方法则是目 前应用最广泛的水文模型结构不确定性分析方法(董磊华,熊立华,万民. 基于贝叶斯模型加权平均方法的水文模型不确定性分析[J].水利学报, 2011,42(9):1065-1074)。然而，研究表明BMA方法实际上仍然是一种加 权平均方法，通常会对模拟流量结果进行平滑化处理，进而降低对流量峰 值的模拟效果。另外，用于求解BMA模型的期望最大化算法需要假设模型 预报流量均服从正态分布，而实际的水文过程往往是非正态的，必须通过 正态分位数转换到正态空间再进行处理，这势必会影响模型的精度和准确 性。此外，当前研究大多均只单独涉及水文模型参数或者模型结构的不确 定性问题，较少同时考虑耦合模型结构和模型参数的综合不确定性。因此， 亟待研究科学有效的水文模型综合不确定性分析方法。

水文模型参数和结构的综合不确定性最终集中表现为模型输出结果的 不确定性。因而，水文模型的综合不确定性分析实质上可以看作是求解给 定模型输出流量时实测流量的条件概率密度函数和分布函数。Copula函数 理论可以将多个随机变量的边缘分布连接起来构造联合分布，求解条件分 布的解析表达式，在水文水资源领域的得到了广泛的应用(郭生练,闫宝 伟,肖义等.Copula函数在多变量水文分析计算中的应用及研究进展[J]. 水文,2008,28(3):1-7)。它可以允许边缘分布为任意分布，能较好地模 拟水文过程的非线性和非正态特征。目前，没有文献将Copula函数引入水 文模型的综合不确定性分析研究中。

发明内容

本发明的目的是克服现有技术存在的不足，提供一种基于Copula函数 的水文模型综合不确定性分析方法。

本发明一种基于Copula函数的水文模型综合不确定性分析方法，包括 以下步骤：

步骤1，收集流域的实测降雨、蒸发和流量数据资料；

步骤2，选择符合流域产汇流特性的水文模型，根据步骤1中的降雨、 蒸发和流量数据资料，采用数学优化算法率定水文模型的参数，利用建立 好的水文模型模拟流域出口断面流量过程；

步骤3，根据步骤1中的实测流量和步骤2中得到的模拟流量数据资 料，选取适当的边缘概率分布函数线型，并估计边缘概率分布函数的参数；

步骤4，采用Copula函数构造实测流量和模拟流量的联合概率分布函 数，并估计Copula函数的参数；

步骤5，根据步骤3估计的边缘概率分布函数和步骤4构建的联合概 率分布函数推求给定模拟流量时实测流量条件概率分布函数的解析表达 式；

步骤6，依据步骤5所得的条件概率分布函数的解析表达式，根据数 理统计原理，计算得到实测流量的中位数作为确定性模拟结果，同时获取 给定置信水平下的不确定性模拟区间。

所述步骤2中，用户根据具体流域的实际情况，选择适当的流域水文 模型结构，为概念性水文模型或者分布式水文模型，包括但不限于新安江 模型、萨克拉门托模型、TOPMODEL模型、TANK模型、VIC模型或MIKESHE 模型。

所述步骤2中，率定水文模型参数时采用参数自动率定方法，选定的 目标函数是残差平方和最小准则，所采用的数学优化算法为SCE-UA算法。

所述步骤3中，将P-III型分布作为实测流量和模拟流量的边缘概率 分布函数线型。

所述步骤3中，采用线性矩法估计边缘概率分布函数的参数。

所述步骤4中，采用FrankCopula函数构造实测流量和模拟流量的联 合概率分布函数，采用Kendall秩相关性系数法估计FrankCopula函数的 参数。

本发明直接对水文模型的输出结果进行后处理，同时考虑了水文模型 结构和参数的不确定性。利用Copula函数构建实测流量和模拟流量的联合 概率分布函数，通过求解给定模拟流量时的实测流量条件概率分布函数， 从而获取实测流量的中位数和不确定性区间，据此分析水文模型的综合不 确定性。

与现有技术相比，本发明的有益效果在于：

1、与常规水文模型不确定性分析方法只能单独考虑模型参数或模型结 构的不确定性不同，本发明可以更加全面地同时考虑模型参数和模型结构 的不确定性，得到水文模型的综合不确定性。

2、本发明独立于确定性水文模型，可以和任意复杂度的确定性水文模 型协同集成，而不需对模型附加任何假定，为分析水文模型的综合不确定 性提供了通用性的理论框架。

3、本发明允许实测流量和模拟流量具有任何形式的边缘概率分布函 数，可以准确地捕捉实测流量和模拟流量之间的非线性和异方差相关性结 构。

4、本发明定量地以概率分布的形式描述水文模型的综合不确定性，并 给出指定置信水平下的不确定性区间，使用户在决策中能定量的考虑各种 不确定性，估计各种决策风险和后果，实现了模拟与决策过程的有机耦合。

附图说明

图1为本发明方法的流程图。

图2为实测流量和模拟流量散点图的示意图。

图3为实测流量理论边缘概率分布函数值(采用P-III型分布计算得 到)与经验边缘概率分布函数值对比情况的示意图。

图4为模拟流量理论边缘概率分布函数值(采用P-III型分布计算得 到)与经验边缘概率分布函数值对比情况的示意图。

图5为采用FrankCopula函数计算得到的理论联合概率分布函数值与 经验联合概率分布函数值对比情况的示意图。

图6为给定模拟流量时实测流量的条件概率分布函数曲线的示意图。

图7为实测流量、根据本发明方法计算得到的中位数模拟结果及90％ 不确定性模拟区间对比情况的示意图。

具体实施方式

下面通过实施例，并结合附图对本发明作进一步说明。

如图1-图7所示，一种基于Copula函数的水文模型综合不确定性分析 方法，收集流域的实测降雨、蒸发和流量资料，建立水文模型模拟流域出 口断面流量过程，在确定实测流量和模拟流量的边缘概率分布函数的基础 上，利用Copula函数构建实测流量和模拟流量的联合概率分布函数，通过 求解给定模拟流量时的实测流量条件概率分布函数，从而获取实测流量的 中位数和不确定性区间。图1是本实施例的计算流程图，按照以下步骤进 行：

1.收集流域的实测降雨、蒸发和流量数据资料。

本具体实施中实测降雨、蒸发和流量数据资料的时间尺度为日。降雨 资料指的是研究流域的面平均降雨量，通过流域上多个代表性降雨站点利 用泰森多边形法计算得到。流域蒸发资料可以从气象站的蒸发皿实测数据 获得。流量资料是指流域出口断面的代表性水文站的实测流量过程，从水 文站的水文年鉴获取。

2.建立水文模型模拟流域出口断面流量过程。

根据流域的气候、地质地貌实际情况，本具体实施中选用新安江模型 作为流域水文模型结构。

根据步骤1中的降雨、蒸发和流量数据资料，本具体实施中采用SCE-UA 优化算法自动率定所选水文模型的参数。

本具体实施中模型参数率定的目标函数是残差平方和最小准则，它可 以描述为实测流量与模拟流量之间的残差平方和最小，如下式所示：

$min\underset{j=1}{\overset{n}{\Sigma }}{(h_j-s_j)}^2---\left(1\right)$

其中，h_j和s_j分别为实测流量和模拟流量，n表示时段数。

如图2所示，给出了实测流量和模拟流量散点图，其中，模拟流量是 通过所率定好的新安江水文模型计算而得。

3.确定实测流量和模拟流量的边缘概率分布函数。

根据步骤1中的实测流量和步骤2中得到的模拟流量数据资料，选取 适当的边缘概率分布函数线型，并估计其参数，本步骤包括两个子步骤：

3.1选择边缘概率分布函数线型

由于实测流量和模拟流量的总体分布频率线型是未知的，通常选用能 较好拟合多数水文样本资料系列的线型。我国经过多年分析比较，发现 P-III型分布对于我国大部分河流的水文资料拟合较好，推荐在工程实践 中采用。

本具体实施中采用P-III型分布作为实测流量和模拟流量的边缘概率 分布函数线型。

3.2估计边缘概率分布函数线型的参数

当频率分布线型选定后，接下来需要进行估计频率分布的参数。目前 常用的方法主要有矩法、极大似然法、适线法、概率权重矩法、权函数法 和线性矩法(L-矩法)等。其中，L-矩法是目前国内外公认的有效参数估 计方法，最大特点是对序列的极大值和极小值没有常规矩那么敏感，求得 的参数估计值比较稳健。

本具体实施中采用L-矩法估计边缘概率分布函数线型的参数。

如图3和图4所示，分别给出了采用P-III型分布计算得到的实测流 量、模拟流量理论边缘概率分布函数值与经验边缘概率分布函数值对比图。 其中，经验边缘概率分布函数值采用一维数学期望公式计算得到。

4.利用Copula函数构建实测流量和模拟流量的联合概率分布函数。

根据步骤1中的实测流量、步骤2中得到的模拟流量数据资料以及步 骤3中估计的边缘概率分布函数，选取适当的Copula函数作为连接函数构 造实测流量和模拟流量的联合概率分布函数，并估计其参数，本步骤包括 两个子步骤：

4.1选择Copula函数

假设H、S分别表示实测流量和模拟流量，h、s分别为相应的实现值。 F_H(h)、F_S(s)是边缘概率分布函数，相应的概率密度函数为f_H(h)、f_S(s)。 由Sklar定理可知，H，S的联合概率分布函数可以用一个二维Copula函 数C表示：

F_H,S(h,s)＝C_θ(F_H(h),F_S(s))＝C_θ(u,v)(2)

其中，θ为Copula函数的参数；u＝F_H(h),v＝F_S(s)为边缘概率分布函 数。

本具体实施中，采用FrankCopula函数构造实测流量和模拟流量的联 合概率分布函数，其表达式如下：

$C_{\theta }(u,v)=-\frac{1}{\theta }ln[1+\frac{(e^{-\theta u}-1)(e^{-\theta v}-1)}{e^{-\theta }-1}],\theta \in R---\left(3\right)$

4.2估计Copula函数的参数

本具体实施中，采用Kendall秩相关性系数法估计FrankCopula函数 的参数。Kendall相关系数τ与参数θ的关系为：

$\tau =1-\frac{4}{\theta }[-\frac{1}{\theta }\underset{-\theta }{\overset{0}{\int }}\frac{t}{\exp \left(t\right)-1}\mathrm{dt}-1]---\left(4\right)$

令{(x₁,y₁),…,(x_n,y_n)}表示从连续随机变量(X,Y)中抽取的n个观测值 的随机样本，则在样本中有种不同的观测值组合(x_i,y_i)和(x_j,y_j)。样 本Kendall秩相关系数τ通过下式计算

$\tau =\frac{2}{n(n-1)}\underset{i=1}{\overset{n-1}{\Sigma }}\underset{j=i+1}{\overset{n}{\Sigma }}\mathrm{sign}\left[(x_i-x_j)(y_i-y_j)\right]---\left(5\right)$

其中，sign(·)是符号函数。

实际计算时，先采用式(5)计算两变量样本的Kendall秩相关系数τ， 再利用式(4)反算出参数θ。由于式(4)不存在解析解，本具体实施中 采用牛顿迭代法求解得到数值解。

如图5所示，给出了采用FrankCopula函数计算得到的理论联合概率 分布函数值与经验联合概率分布函数值的对比情况。其中，经验边缘概率 分布函数值采用二维数学期望公式计算得到。

5.求解给定模拟流量时的实测流量条件概率分布函数。

给定模拟流量S取值s时，所对应的实测流量H的取值并非唯一，而是可 大可小，只是出现不同取值的概率有所不同，存在着一个条件概率分布函 数

F_H|S(h)＝P(H≤h|S＝s)(6)

借助Copula函数，条件概率分布函数F_H|S(h)可以表示为：

$F_{H|S}\left(h\right)=P(U\le u|V=v)=\frac{\partial C_{\theta }(u,v)}{\partial v}=\frac{1}{\theta }\frac{e^{-\theta v}(e^{-\theta u}-1)}{(e^{-\theta }-1)+(e^{-\theta u}-1)(e^{-\theta v}-1)}---\left(7\right)$

如图6所示，给出了给定模拟流量时实测流量的条件概率分布函数曲 线。

6.获取实测流量的中位数和不确定性区间。

条件概率分布函数F_H|S(h)相应的概率密度函数为

f_H|S(h)＝c_θ(u,v)f_H(h)(8)

其中，c_θ(u,v)为Copula函数的密度函数，解析表达式为：

$c_{\theta }(u,v)=\frac{{\partial }^2{C}_{\theta }(u,v)}{\partial u\partial v}=\frac{-\theta (e^{-\theta }-1)e^{-\theta (u+v)}}{{[(e^{-\theta }-1)+(e^{-\theta u}-1)(e^{-\theta v}-1)]}^2}---\left(9\right)$

得到随机变量H的条件概率密度函数f_H|S(h)后，根据数理统计原理， 可以计算得到实测流量的中位数作为确定性模拟结果，同时获取给定置信 水平下的不确定性模拟区间。

实测流量的中位数h_m通过下式求解：

${\int }_0^{h_m}{f}_{H|S}\left(h\right)dh=0.5---\left(10\right)$

选择一定的置信水平(1-ξ)，令实测流量取值出现在分布两端的概率 为ξ，就可以定义实测流量的区间估计。随机变量H的置信下、上限分别 由以下两式给出：

${\int }_0^{h_l}{f}_{H|S}\left(h\right)dh={\xi }_1---\left(11\right)$

${\int }_0^{h_u}{f}_{H|S}\left(h\right)dh=1-{\xi }_2---\left(12\right)$

式中：ξ₁＋ξ₂＝ξ，为显著性水平；ξ₁和ξ₂可以根据实际问题任意选定，本具 体实施中取ξ₁＝ξ₂＝ξ/2。

因此有

P(h_l≤H≤h_u)＝1-ξ(13)

即[h_l,h_u]为随机变量H的置信水平(1-ξ)的区间估计，根据置信区间 可以对实测流量估计值的不确定性进行定量评价。

考虑到式(10)-(12)均无法获得解析解，本具体实施中采用二分法 试算求解得到数值解。

如图7所示，给出了实测流量、根据本发明方法计算得到的中位数模 拟结果及90％不确定性模拟区间对比情况。其中，中位数模拟结果为条件 概率分布函数50％的分位数；给定显著性水平ξ＝0.1，计算得到条件概率分 布函数5％和95％的分位数，它们分别给出了90％不确定性模拟区间的置信 下限和上限值。

综上，本发明基于流域实测降雨、蒸发和流量资料，建立水文模型模 拟流域出口断面流量过程，在确定实测流量和模拟流量的边缘概率分布函 数的基础上，利用Copula函数构建实测流量和模拟流量的联合概率分布函 数，通过求解给定模拟流量时的实测流量条件概率分布函数，从而获取实 测流量的中位数和不确定性区间。本发明以概率分布的形式描述水文模型 的综合不确定性，给出指定置信水平下的不确定性区间，使用户在决策中 能定量的考虑各种不确定性，实现了模拟与决策过程的有机耦合。