CN110377989A - 基于水热耦合平衡的非一致性两变量设计洪水推求方法 - Google Patents

Abstract

本发明涉及水库防洪安全设计技术领域，公开了基于水热耦合平衡的非一致性两变量设计洪水推求方法，通过收集研究流域的日降水、气温、风速、相对湿度、日照时长、实测蒸散发及多年长序列观测资料，建立流域水热耦合平衡方程，分析流域年均下垫面特征参数；依据水热耦合平衡方程推求非一致性条件下洪水单变量分布函数；基于洪水峰量相关结构的非一致性特征，采用Gumbel‑Hougaard Copula(G‑H Copula)方法建立非一致性条件下的洪峰‑洪量联合概率分布函数，求取非一致性条件下设计洪水峰量设计值。本发明能够充分考虑气候变化和人类活动对洪水特征的影响，同时采用洪峰‑洪量最可能组合模型，为变化环境下非一致性水库两变量设计洪水的推求提供重要且可行的参考依据。

Description

基于水热耦合平衡的非一致性两变量设计洪水推求方法

技术领域

本发明涉及水库防洪安全设计技术领域，具体涉及基于水热耦合平衡的非一致性两变量设计洪水推求方法。

背景技术

设计洪水是确定水利水电工程建设规模及制定运行管理策略的重要依据，对指导水库运行调度和水资源综合规划利用具有重要意义。我国传统的设计洪水计算一般假定水文系列服从一致性，且采用峰量同频率组合法来确定设计值。同频率组合法未能充分挖掘水文事件的内在规律，难以科学评估洪水风险。此外，由于气候变化及人类活动的影响，使得下垫面条件发生改变，影响了水文循环过程，改变了不同时期水文极值(如洪水)的时空分布规律，如何推求变化环境下的自适应设计洪水成为当前的研究难点。

为了解决以上问题，授权号为ZL201510106099.8的发明专利公开了基于多变量条件最可能组合推求设计洪水的方法，该专利基于Copula函数和单变量条件重现期提出了多变量条件最可能组合法，克服了传统同频率假设的局限，但该发明仍然假定水文系列服从一致性条件。授权号为ZL201610111489.9的发明专利公开了一种自适应气候变化的水库多变量设计洪水推求方法，该方法采用全球气候模式、统计降尺度模型和水文模型，考虑了未来气候变化情景下的设计洪水最可能组合模式，也未能考虑变化环境对洪水边缘分布以及相关性结构的影响。

水热平衡方程Budyko公式充分考虑区域的水量平衡和能量平衡，能够较好地反映区域下垫面条件和产汇流特征，近年来广泛应用于研究流域水循环过程。但是，目前应用Budyko公式开展变化环境下非一致性两变量设计洪水计算的研究未见文献报道。

发明内容

基于以上问题，本发明提供了基于水热耦合平衡的非一致性两变量设计洪水推求方法，能够充分考虑气候变化和人类活动对洪水特征的影响，同时采用洪峰-洪量最可能组合模型，为变化环境下非一致性水库两变量设计洪水的推求提供重要且可行的参考依据。

为了实现上述目的，本发明提供了基于水热耦合平衡的非一致性两变量设计洪水推求方法，包括如下步骤：

步骤1、基础气象水文数据采集与取样：收集流域的日降水、气温、风速、相对湿度、日照时长、实测蒸散发及控制水文站流量的多年长序列观测资料，从中统计年最大洪峰、年最大7日径流量、饱和水气压、实际水气压、地表净辐射、土壤热通量、干湿表常数；

步骤2、基础数据处理：采用Budyko公式的特征参数作为解释变量，建立流域水热耦合平衡方程，获取流域时变下垫面特征参数年均序列；

步骤3、建立边缘概率分布：基于水热耦合平衡方程分别推求非一致性条件下洪峰和洪量的单变量分布函数；

步骤4、建立联合概率分布：针对洪水峰量相关结构的非一致性，以水热耦合平衡方程中下垫面特征参数为解释变量，建立非一致性条件下基于G-H Copula的联合概率分布函数，求取非一致性条件下考虑洪峰-洪量相关关系的水库洪水设计值；

步骤5、输出洪水最可能联合设计值：考虑洪水峰量最可能组合模式，基于G-Hcopula函数，推求气候变化与人类活动影响下不同重现期的水库两变量洪水设计值。

进一步地，步骤2包括以下子步骤：

2a：通过彭曼公式计算潜在蒸散量PET，再通过水量平衡方程获取流域实际蒸散发ET；所述彭曼公式为：

所述水量平衡方程为ET＝P-R，P为实际降水量，R为实测径流量；

式中：PET为潜在蒸散量(mm/d)；Δ代表饱和水气压曲线斜率(kPa/℃)；R_n为地表净辐射(MJ/(m·d))；G为土壤热通量(MJ/(m2·d))；γ为干湿表常数(kPa/℃)；T_mean代表日平均温度(℃)；u₂指距离地表2米高处风速(m/s)；e_s指饱和水气压(kPa)；e_a指实际水气压(kPa)。

2b：选取一定的时间窗口，通过最小二乘法，率定水热耦合平衡方程的下垫面特征参数年均值序列，记为ω；分析下垫面特征参数ω与洪水系列的相关性，用以验证ω作为洪水序列解释变量的可靠性；所述年均水热耦合平衡方程为：

进一步地，步骤3包括以下内容：

3a：基于步骤1中取样所得年最大洪峰和洪量长系列观测值，以常用于洪水频率分析的P-III型概率分布函数为洪水单变量理论分布，通过时变矩法分别建立洪峰Q、洪量W单变量边缘分布参数的时变模型；其中，洪峰Q采用年最大洪峰，洪量W采用年最大7日径流量；一致性条件下，三参数P-III型分布的概率密度函数f(x)为：

式中，α为尺度参数，β为形状参数，μ为位置参数。

对于非一致性条件下，P-III型分布参数随时间发生变化，记t时刻对应的尺度参数为α_t、形状参数为β_t、位置参数为μ_t，假设三个时变参数的时变矩可通过解释变量ω^t的单调函数g(θx)来表示，即：

式中：g(θx)表示单调连续函数，具体形式通过统计参数θ_X的定义域来确定；当θ_x≤0时，g(θ_x)＝θ_x，当θ_X＞0时,则g(θ_X)＝ln(θ_X)；ω^t表示解释变量(budyko方程中下垫面特征参数，ω在时间t状态下的值；α_i(i＝1,2,3；10,20,30)表示模型的参数，其取值根据GAMLSS模型率定；

3b：将步骤3a的率定结果代入非一致性条件下P-III型分布的概率密度函数f(x)，分别获得洪峰和洪量的理论边缘分布；所述非一致性条件下P-III型分布形式为

进一步地，步骤4包括以下子步骤：

4a：基于步骤3b构建的洪峰Q与洪量W理论边缘分布函数，通过Gumbel-HougaardCopula(G-H Copula)方法，构造洪峰Q与洪量W的联合概率分布函数；所述一致性条件下G-HCopula方法表达式如下：

式中：为一致性条件下Copula联合分布函数，θ_c范围Copula联合分布函数参数，定义域为(1,∞)；u^t,v^t分别为洪峰Q、洪量W时变边缘分布函数

将G-H Copula联合分布函数的参数θ_c替换为时变参数联合步骤3中的洪峰Q、洪量W时变边缘分布，基于Copula函数的定义，非一致性时变两变量Copula函数可表示为：

式中：F^t(q^t,w^t)表示时变情形下洪峰Q和洪量W的联合分布函数；和分别表示Q和W变量的时变边缘分布函数和时变参数；代入G-H Copula表达式，可获得如下表达式：

其中，为时变情形下Copula联合分布函数F^t(q^t,w^t)，为时变参数，其范围为(1,∞)；u^t,v^t意义同上；

4b：使用G-H Copula函数建立其时变参数与解释变量之间的线性关系；所述线性关系可表示为：

式中：表示G-H Copula函数的联结函数；当时,b₀,b₁∈R分别为模型的参数，在无特别说明情况下，参数范围均为实数范围，取值根据GAMLSS模型率定。

进一步地，步骤5包括以下子步骤：

5a：记和分别表示洪峰Q和洪量W的时变边缘分布函数，基于步骤4b所得洪峰Q和洪量W的时变联合分布函数F^t(q^t,w^t)，则有对应于本发明采用的设计洪水两变量统计值(洪峰和洪量)，本发明采用OR(或)重现期作为水库防洪标准的度量指标；其定义如下：

T_or ^t(q^t,w^t)＝1/[1-F^t(q^t,w^t)]；

式中：T_or ^t(q^t,w^t)为时变OR联合重现期，以年为单位；

5b：基于给定的联合重现期，本发明采用最可能组合模式从时变联合概率分布中计算设计洪峰-洪量组合；所述最可能组合模式是指在重现期等值线上联合概率密度函数最大的组合(q^*(t),w^*(t))，为下述方程的解：

式中：f^t(q^t,w^t)表示Q和W时变联合分布的密度函数；表示时变Copula函数的密度函数；和分别表示和的密度函数；

5c：采用拉格朗日乘数法求解步骤5b中建立的最可能组合模式方程，构造以下求解方程：

式中：λ_t表示时间状态t对应的拉格朗日乘子。

上述求解方程可转化为如下方程组，其解即为联合重现期对应的水库洪水的洪峰、洪量设计值：

式中： 分别为洪峰和洪量的概率密度函数的导函数。该非线性方程组为基于Copula函数推求的联合重现期水平下峰量最可能组合法应满足的计算通式。本发明采用多变量同频率组合计算得到的峰量设计结果作为初始解，采用调和平均数牛顿法进行迭代求解，得到水库设计洪水峰量最可能组合(q*,w*)。

与现有技术相比，本发明有如下优点和效果：

1)科学合理、贴近工程实际：本发明充分考虑水文系列的非一致性特征，采用Budyko公式的特征参数作为解释变量，考虑水文系列的非一致性构建时变Copula模型，能够充分考虑气候变化和人类活动对洪水特征的影响。

2)考虑最可能组合法，可为变化环境下非一致性水库两变量设计洪水的推求提供重要且可行的参考依据。

附图说明

图1为基于水热耦合平衡的非一致性两变量设计洪水推求方法的流程图；

图2为Budyko水热平衡方程的示意图；

图3为对应的重现期等值线及峰量组合模式的示意图；

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本发明，并不作为对本发明的限定。

参见图1，基于水热耦合平衡的非一致性两变量设计洪水推求方法，包括以下步骤：

步骤1，基础气象水文数据采集与取样

收集流域的日降水、气温、风速、相对湿度、日照时长、实测蒸散发及控制水文站流量的多年长序列观测资料，从中统计年最大洪峰、年最大7日径流量、饱和水气压、实际水气压、地表净辐射、土壤热通量、干湿表常数。

步骤2，基础数据处理：采用Budyko公式的特征参数作为解释变量，建立流域水热耦合平衡方程，获取流域时变下垫面特征参数年均序列。

如图2所示，给出了Budyko水热平衡方程的示意图；该步骤包括的子步骤有：

2a：通过彭曼公式计算潜在蒸散量PET，再通过水量平衡方程获取流域实际蒸散发ET；所述彭曼公式为：

所述水量平衡方程为ET＝P-R，P为实际降水量，R为实测径流量；

式中：PET为潜在蒸散量(mm/d)；Δ代表饱和水气压曲线斜率(kPa/℃)；R_n为地表净辐射(MJ/(m·d))；G为土壤热通量(MJ/(m²·d))；γ为干湿表常数(kPa/℃)；T_mean代表日平均温度(℃)；u₂指距离地表2米高处风速(m/s)；e_s指饱和水气压(kPa)；e_a指实际水气压(kPa)。

2b：选取一定的时间窗口，通过最小二乘法，率定水热耦合平衡方程的下垫面特征参数年均值序列，记为ω；分析下垫面特征参数ω与洪水系列的相关性，用以验证ω作为洪水序列解释变量的可靠性；所述年均水热耦合平衡方程为：

特别地，本实施例选取的时间窗口为11年。

步骤3，基于水热耦合平衡方程分别推求非一致性条件下洪峰和洪量的单变量分布函数：

基于步骤1中取样所得年最大洪峰和洪量长系列观测值，以常用于洪水频率分析的P-III型概率分布函数为洪水单变量理论分布，通过时变矩法分别建立洪峰Q、洪量W单变量边缘分布参数的时变模型；其中，洪峰Q采用年最大洪峰，洪量W采用年最大7日径流量；一致性条件下，三参数P-III型分布的概率密度函数f(x)为：

式中，α为尺度参数，β为形状参数，μ为位置参数。

对于非一致性条件下，P-III型分布参数随时间发生变化，记t时刻对应的尺度参数为α_t、形状参数为β_t、位置参数为μ_t，假设三个时变参数的时变矩可通过解释变量ω^t的单调函数g(θx)来表示，即：

式中：g(θx)表示单调连续函数，具体形式通过统计参数θ_X的定义域来确定；当θ_x≤0时，g(θ_x)＝θ_x，当θ_X＞0时,则g(θ_X)＝ln(θ_X)；ω^t表示解释变量(budyko方程中下垫面特征参数，ω在时间t状态下的值；α_i(i＝1,2,3；10,20,30)表示模型的参数，其取值根据GAMLSS模型率定；

3b：将步骤3a的率定结果代入非一致性条件下P-III型分布的概率密度函数f(x)，分别获得洪峰和洪量的理论边缘分布；所述非一致性条件下P-III型分布形式为

步骤4，基于水热耦合平衡建立非一致性条件下两变量洪水联合概率分布函数：

选取在一致性条件下对洪水系列相关结构拟合较好的Gumbel-Hougaard Copula(简称G-H Copula)函数为联合概率分布函数，针对洪水峰量相关结构的非一致性，仍采用热耦合平衡方程中参数ω为解释变量，建立非一致性条件下基于Copula的联合概率分布函数，求取非一致性条件下设计洪水峰量设计值。

子步骤主要包括：

4a：基于步骤3b构建的洪峰Q与洪量W理论边缘分布函数，通过Gumbel-HougaardCopula(G-H Copula)方法，构造洪峰Q与洪量W的联合概率分布函数；所述G-H Copula方法表达式如下：

式中：为一致性条件下Copula联合分布函数，θ_c范围Copula联合分布函数参数，定义域为(1,∞)；u^t,v^t分别为洪峰Q、洪量W时变边缘分布函数

将G-H Copula联合分布函数的参数θ_c替换为时变参数联合步骤3中的洪峰Q、洪量W时变边缘分布，基于Copula函数的定义，非一致性时变两变量Copula函数可表示为：

式中：F^t(q^t,w^t)表示时变情形下洪峰Q和洪量W的联合分布函数；和分别表示Q和W变量的时变边缘分布函数和时变参数；代入G-H Copula表达式，可获得如下表达式：

其中，为时变情形下Copula联合分布函数F^t(q^t,w^t)，为时变参数，其范围为(1,∞)；u^t,v^t意义同上；

4b：使用G-H Copula函数建立其时变参数与解释变量之间的线性关系；所述线性关系可表示为：

式中：表示G-H Copula函数的联结函数；当时,b₀,b₁∈R分别为模型的参数，在无特别说明情况下，参数范围均为实数范围，取值根据GAMLSS模型率定。

步骤5、输出洪水最可能联合设计值计算结果：考虑洪水峰量最可能组合模式，基于G-H copula函数，推求不同重现期水平下，气候变化与人类活动影响下的水库两变量洪水设计值。

基于Copula函数，建立考虑下垫面人类活动的水库设计洪水峰量最可能组合法计算模型。

5a、记和分别表示洪峰Q和洪量W的时变边缘分布函数，基于步骤4b所得洪峰Q和洪量W的时变联合分布函数F^t(q^t,w^t)，则有对应于本发明采用的设计洪水两变量统计值(洪峰和洪量)，本发明采用OR(或)重现期作为水库防洪标准的度量指标；其定义如下：

T_or ^t(q^t,w^t)＝1/[1-F^t(q^t,w^t)] (10)

式中：T_or ^t(q^t,w^t)为时变OR联合重现期，以年为单位。

如图3所示，给出了重现期等值线及峰量组合模式的示意图；图中给出了常用的同频率组合联合设计值，以及本实施例采用的最可能组合模式。

5b、基于给定的联合重现期，本发明采用最可能组合模式从时变联合概率分布中计算设计洪峰-洪量组合；所述最可能组合模式是指在重现期等值线上联合概率密度函数最大的组合(q^*(t),w^*(t))，为下述方程的解：

式中：f^t(q^t,w^t)表示Q和W时变联合分布函数的密度函数；

表示时变copula函数的密度函数；和分别表示和的密度函数。

5c、本实施例提出采用拉格朗日乘数法求解最可能组合问题的方法，给定联合重现期T_or，构造以下求解方程：

式中：λ_t表示时间状态t对应的拉格朗日乘子。

上述求解方程可转化为如下方程组，其解即为联合重现期对应的水库洪水的洪峰、洪量设计值：

式中：； 分别为相应概率密度函数的导函数。

非线性方程组(13)为基于Copula函数推求的联合重现期水平下峰量最可能组合法应满足的计算通式。本发明采用多变量同频率组合计算得到的峰量设计结果作为初始解，采用调和平均数牛顿法进行迭代求解，得到水库设计洪水峰量最可能组合(q*,w*)。

本实施例的计算思路为：通过流域实测气象变量(气温、风速、相对湿度、日照时长)长系列，基于彭曼公式获取流域潜在蒸散发数据；通过水量平衡方程推求流域实际蒸散发长系列；基于budyko公式建立考虑流域水循环的水热平衡模型。利用该水热平衡模型中代表流域产汇流特征的w参数，基于时变矩法，建立变化环境下考虑水文系列非一致性的时变Copula模型。考虑洪水峰量最可能组合模式，基于G-H copula函数，推求不同重现期水平下，气候变化与人类活动影响下的水库两变量洪水设计值，为变化环境下非一致性水库两变量设计洪水的推求提供重要且可行的参考依据。

如上即为本发明的实施例。上述实施例以及实施例中的具体参数仅是为了清楚表述发明验证过程，并非用以限制本发明的专利保护范围，本发明的专利保护范围仍然以其权利要求书为准，凡是运用本发明的说明书及附图内容所作的等同结构变化，同理均应包含在本发明的保护范围内。

Claims (5)

1.基于水热耦合平衡的非一致性两变量设计洪水推求方法，其特征在于包括以下步骤：

步骤1、基础气象水文数据采集与取样：收集流域的日降水、气温、风速、相对湿度、日照时长、实测蒸散发及控制水文站流量的多年长序列观测资料，从中统计年最大洪峰、年最大7日径流量、饱和水气压、实际水气压、地表净辐射、土壤热通量、干湿表常数；

步骤2、基础数据处理：采用Budyko公式的特征参数作为解释变量，建立流域水热耦合平衡方程，获取流域时变下垫面特征参数年均序列；

步骤3、建立边缘概率分布：基于水热耦合平衡方程分别推求非一致性条件下洪峰和洪量的单变量分布函数；

步骤4、建立联合概率分布：针对洪水峰量相关结构的非一致性，以水热耦合平衡方程中下垫面特征参数为解释变量，建立非一致性条件下基于G-H Copula的联合概率分布函数，求取非一致性条件下考虑洪峰-洪量相关关系的水库洪水设计值；

步骤5、输出洪水最可能联合设计值：考虑洪水峰量最可能组合模式，基于G-H copula函数，推求气候变化与人类活动影响下不同重现期的水库两变量洪水设计值。

2.如权利要求1所述的方法，其特征在于步骤2包括以下子步骤：

2a：通过彭曼公式计算潜在蒸散量PET，再通过水量平衡方程获取流域实际蒸散发ET；所述彭曼公式为：

所述水量平衡方程为ET＝P-R，P为实际降水量，R为实测径流量；

式中：PET为潜在蒸散量(mm/d)；Δ代表饱和水气压曲线斜率(kPa/℃)；R_n为地表净辐射(MJ/(m·d))；G为土壤热通量(MJ/(m²·d))；γ为干湿表常数(kPa/℃)；T_mean代表日平均温度(℃)；u₂指距离地表2米高处风速(m/s)；e_s指饱和水气压(kPa)；e_a指实际水气压(kPa)。

2b：选取一定的时间窗口，通过最小二乘法，率定水热耦合平衡方程的下垫面特征参数年均值序列，记为ω；分析下垫面特征参数ω与洪水系列的相关性，用以验证ω作为洪水序列解释变量的可靠性；所述年均水热耦合平衡方程为：

3.如权利要求1所述的方法，其特征在于步骤3包括以下子步骤：

3a：基于步骤1中取样所得年最大洪峰和洪量长系列观测值，以常用于洪水频率分析的P-III型概率分布函数为洪水单变量理论分布，通过时变矩法分别建立洪峰Q、洪量W单变量边缘分布参数的时变模型；其中，洪峰Q采用年最大洪峰，洪量W采用年最大7日径流量；一致性条件下，三参数P-III型分布的概率密度函数f(x)为：

式中，α为尺度参数，β为形状参数，μ为位置参数。

对于非一致性条件下，P-III型分布参数随时间发生变化，记t时刻对应的尺度参数为α_t、形状参数为β_t、位置参数为μ_t，假设三个时变参数的时变矩可通过解释变量ω^t的单调函数g(θx)来表示，即：

式中：g(θx)表示单调连续函数，具体形式通过统计参数θ_X的定义域来确定；当θ_x≤0时，g(θ_x)＝θ_x，当θ_X＞0时,则g(θ_X)＝ln(θ_X)；ω^t表示解释变量(budyko方程中下垫面特征参数，ω在时间t状态下的值；α_i(i＝1,2,3；10,20,30)表示模型的参数，其取值根据GAMLSS模型率定；

3b：将步骤3a的率定结果代入非一致性条件下P-III型分布的概率密度函数f(x)，分别获得洪峰和洪量的理论边缘分布；所述非一致性条件下P-III型分布形式为

4.如权利要求1所述的方法，其特征在于步骤4包括以下子步骤：

4a：基于步骤3b构建的洪峰Q与洪量W理论边缘分布函数，通过Gumbel-HougaardCopula(G-H Copula)方法，构造洪峰Q与洪量W的联合概率分布函数；所述G-H Copula方法表达式如下：

式中：为一致性条件下Copula联合分布函数，θ_c范围Copula联合分布函数参数，定义域为(1,∞)；u^t,v^t分别为洪峰Q、洪量W时变边缘分布函数

将G-H Copula联合分布函数的参数θ_c替换为时变参数联合步骤3中的洪峰Q、洪量W时变边缘分布，基于Copula函数的定义，非一致性时变两变量Copula函数可表示为：

式中：F^t(q^t,w^t)表示时变情形下洪峰Q和洪量W的联合分布函数；和分别表示Q和W变量的时变边缘分布函数和时变参数；代入G-H Copula表达式，可获得如下表达式：

其中，为时变情形下Copula联合分布函数F^t(q^t,w^t)，为时变参数，其范围为(1,∞)；u^t,v^t意义同上；

4b：使用G-H Copula函数建立其时变参数与解释变量之间的线性关系；所述线性关系可表示为：

式中：表示G-H Copula函数的联结函数；当时,b₀,b₁∈R分别为模型的参数，在无特别说明情况下，参数范围均为实数范围，取值根据GAMLSS模型率定。

5.如权利要求1所述的方法，其特征在于步骤5包括以下子步骤：

5a：记和分别表示洪峰Q和洪量W的时变边缘分布函数，基于步骤4b所得洪峰Q和洪量W的时变联合分布函数F^t(q^t,w^t)，则有对应于本发明采用的设计洪水两变量统计值(洪峰和洪量)，本发明采用OR(或)重现期作为水库防洪标准的度量指标；其定义如下：

T_or ^t(q^t,w^t)＝1/[1-F^t(q^t,w^t)]；

式中：T_or ^t(q^t,w^t)为时变OR联合重现期，以年为单位；

5b：基于给定的联合重现期，本发明采用最可能组合模式从时变联合概率分布中计算设计洪峰-洪量组合；所述最可能组合模式是指在重现期等值线上联合概率密度函数最大的组合(q^*(t),w^*(t))，为下述方程的解：

式中：f^t(q^t,w^t)表示Q和W时变联合分布的密度函数；表示时变Copula函数的密度函数；和分别表示和的密度函数；

5c：采用拉格朗日乘数法求解步骤5b中建立的最可能组合模式方程，构造以下求解方程：

式中：λ_t表示时间状态t对应的拉格朗日乘子。

上述求解方程可转化为如下方程组，其解即为联合重现期对应的水库洪水的洪峰、洪量设计值：

式中： 分别为洪峰和洪量的概率密度函数的导函数。