CN113469528A - 一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法 - Google Patents

Abstract

本发明公开了一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，包括S1、水库群防洪调度不确定性描述；S2、水库群多目标防洪调度模型建立及随机模拟；S3、水库群防洪调度风险分析与计算。本发明建立了考虑多不确定性、多目标的水库群防洪调度风险分析模型，耦合了洪水资源利用所导致的起调水位不确定性风险源，更符合我国汛期水库调度的实际情况；综合考虑了水库和下游防洪目标，可进行防洪风险在上下游的竞争转化关系分析，剖析水库群系统防洪风险的内在演变规律，为防洪调度策略制定和风险管理提供理论参考。

Description

一种考虑时空关联多维不确定性的水库群多目标防洪调度风 险分析方法

技术领域

本发明涉及水库群防洪调度风险分析方法，尤其是一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法。

背景技术

洪水作为主要自然灾害之一，严重影响着人类生命和财产安全。以防洪为主要目的的水库工程在过去多年间逐渐形成由点到面的全流域综合水库群系统，多水库联合调度能够充分利用各水库之间的防洪库容补偿调节作用，对于实现高效的洪水管理和水资源利用发挥着重要作用。但是水库群防洪调度面临着诸多不确定性，基于风险调控和管理的防洪调度对于耦合多种风险因子、确保洪水安全管理具有重要意义。

水库群实时防洪调度依据发布的洪水预报信息制定预报调度方案。由于降雨输入、模型结构、模型参数等不确定因素，洪水预报也不可避免地带有不确定性，预报低估会导致水库水位过高或下泄流量过大，从而造成风险事件。洪水预报误差是水库防洪调度的主要风险因子。实际上，在日益严峻的水资源供需矛盾下，水库工程汛期在以保证防洪安全为主要任务的前提下，也要兼顾洪水资源利用，充分发挥流域工程群的综合效益。目前的洪水资源利用主要以水库汛限水位动态控制为手段，即在预报无雨时，根据水库预报期内的预泄能力确定水位浮动范围，拦蓄洪水尾水进行超蓄。这种计算水位阈值的方法称为“预报预泄法”。由于预报误差的存在，基于预报信息所确定的汛限水位上限值也具有不确定性。若预报准确，在下一场洪水发生时水库水位可安全回落到汛限水位，而当降雨预报漏报或预报入流偏低时，原蓄存在防洪库容中的洪水无法及时下泄，将导致洪水发生时水库水位高于设计汛限水位，造成额外的防洪风险。目前的防洪风险研究大多以设计汛限水位作确定性起调水位，鲜少有考虑洪水资源利用所导致的起调水位不确定性。除此之外，由于工程测量、模型概化等误差导致的水位库容关系、水库泄流能力和洪水演进不确定性也是防洪调度过程中重要的风险因子。耦合各风险因子的调洪风险分析具有重要的学术意义和实际应用价值。

流域水库群防洪工程系统常以全流域、多时段信息进行流域级联合调度，涉及的多种不确定性因子受到水文、水力复杂关联关系影响也相互关联，比如不同时间(时)和不同水库(空)的预报误差，以及受预报误差影响的各水库的起调水位误差。忽视风险源内部的依赖关系可能导致防洪风险的低估，尤其是流域级大洪水所致的系统群发性风险事件。Copula函数是一种成熟的描述多变量相关性的有效手段，本发明采用Copula函数建立上述风险源的联合分布函数，用以表征风险因子的时空关联特征。

另外，水库群防洪调度本身也是协调上游与下游安全的多目标问题。将洪水尽量蓄存在水库则会增加大坝和库区的防洪损失，而将洪水过多地下泄出去又会造成下游河道的防洪压力增加。分析防洪风险在上下游之间的协调转化关系对于制定水库群防洪调度方案具有重要意义。

发明内容

发明目的：本发明的目的是提供一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，综合考虑了水库和下游防洪目标，可进行防洪风险在上下游的竞争转化关系分析，剖析水库群系统防洪风险的内在演变规律，为防洪调度策略制定和风险管理提供理论参考。

技术方案：本发明的一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，包括以下步骤：

S1、水库群防洪调度不确定性描述；

防洪调度不确定性包括水库起调水位不确定性，洪水预报的水文不确定性，水库泄洪的水力不确定性，水位库容关系不确定性和洪水传播演进不确定性；

S2、水库群多目标防洪调度模型建立及随机模拟；

调度模型考虑上下游两个目标，上游防洪目标以一次调洪中占用的防洪库容最小为原则，下游防洪目标以各水库最大泄流最小为原则；

调度模型约束条件包括：水量平衡约束、水库水位约束、出库流量约束、泄流能力约束、出库流量变幅约束和初始、边界条件；

采用NSGA-II智能算法求解水库群多目标防洪调度模型，将各水库出库流量作为决策变量，采用水库防洪调度中常用的分段试算法进行约束处理，种群数为非劣解数目；

S3、水库群防洪调度风险分析与计算；

水库群防洪调度风险分别从各水库的防洪风险、各下游防洪点的防洪风险和系统总风险进行定义，其中，水库和防洪点风险定义为超安全阈值的累积概率，系统风险定义为水库最大风险与下游最大风险之和。

进一步的，步骤S1中水库起调水位不确定性描述具体为：

S111、基于预报预泄法和水库聚合分解获取各水库起调水位的经验样本；预报预泄是指在预报无雨时，无雨期内以水库兴利用水做预泄流量，下一场洪水起涨的有效预见期内以下游安全泄量做预泄流量，以计算洪水退水段可超蓄洪水资源量，公式表示如下：

其中，ΔW为超蓄水量上限；q_m为水库兴利用水流量；T_d为预报无雨时长；T_f为洪水预报有效预见期；α为安全泄量折扣系数；q_s为下游河道安全泄量；Q_fl(t)为预报流量过程；当前时刻记为0；

记实际所需超蓄水量为ΔW_a，则在以预报信息计算得到的ΔW超蓄的前提下，以实际情况预泄ΔW_a，则起调水位Z₀为：

Z₀＝f(V_x+ΔW-ΔW_a) (2)；

其中，V_x为设计汛限水位对应的水库蓄水量；f(·)为水位库容关系；

对于包含m个水库的水库群系统，采用聚合分解法计算各水库的起调水位样本值：

其中S为样本数量；

S112、以Copula函数建立m个水库的起调水位Z_0,1,Z_0,2,...,Z_0,m的联合分布函数以进行随机模拟；具体为：

S1121、根据Z_0,1,Z_0,2,...,Z_0,m的经验样本值，采用非参数核密度估计法计算各水库起调水位的边缘密度函数f_i(Z_0,i)的核密度估计：

其中，

为水库i起调水位的第s个样本值，K(·)为核函数，h为窗宽；

采用R语言的kernel函数进行核密度及核分布估计值计算；

S1122、采用最大似然估计拟合备选的各Copula函数模型的参数；

S1123、计算拟合优度统计量赤池信息准则、均方根误差和KS检验统计量，以进行Copula模型优选；

赤池信息准则计算公式为：

AIC＝-2ln(L)+2v (4)；

均方根误差计算公式为：

KS检验统计量计算公式为：

D_max＝max|P_e,s-P_s| (6)；

其中，L为似然函数，v为分布参数个数，P_e和P分别为经验Copula函数和理论Copula函数，P_e,s和P_s分别为第s个样本的经验Copula函数和理论Copula函数值；

经验Copula函数定义如下：

其中，F_i(·)为各库起调水位的经验分布函数，I_[·]为示性函数，当

时，

否则

综合选取AIC、RMSE和D_max统计值最小的Copula模型为最优Copula模型；

S1124、采用R语言copula包中的rCopula函数，根据最优的Copula模型生成服从该联合分布的m个具有相关关系的、服从[0,1]均匀分布的随机向量θ₁,θ₂,…,θ_m，向量长度为L，即模拟组数；

S1125、根据所得的随机向量θ₁,θ₂,…,θ_m，对各变量边缘分布逆向计算求得各水库的起调水位模拟值

为水库i第l个起调水位模拟值，L为模拟组数。

进一步的，步骤S1中洪水预报的水文不确定性描述具体为：

S121、定义第i个水库、第j个防洪点区间在t时刻的预报误差e_i(t)、e_j(t)分别为预报值与实测值之差：

其中，e_i(t)为第i个水库在t时刻的预报误差，e_j(t)为第j个防洪点区间在t时刻的预报误差，QR_for,i(t)和QR_obs,i(t)分别为第i个水库所辖集水区间的预报来水和实际来水；QL_for,j(t)和QL_obs,j(t)分别为第j个防洪点区间的预报来水和实际来水；T为预报系统发布的预报总时长；

S122、对于包含m个水库、n个下游防洪点区间的水库群系统，预报误差表示为以下(m+n)*T的矩阵：

S123、采用t-Copula建立时空关联的预报误差联合分布，并基于Copula模型进行随机模拟，具体为：

S1231、根据历史预报和实测数据统计各集水区间的洪水预报误差，得到误差序列样本，同样采用非参数法建立其边缘分布；

S1232、采用最大似然法估计t-Copula函数的参数；

S1233、生成服从[0,1]均匀分布的各集水区间、各预见期的误差分布函数值向量θ₁(1),…,θ₁(T)；θ₂(1),…,θ₂(T)；…；θ_m+n(1),…,θ_m+n(T)，向量长度为L，即模拟组数；

S1234、根据上步的分布函数值向量和各变量的边缘分布，逆向计算得到各集水区间、各预见期的预报误差模拟值e₁(1)^l,…,e₁(T)^l；e₂(1)^l,…,e₂(T)^l；…；e_m+n(1)^l,…,e_m+n(T)^l；l＝1,2,…,L，e_i(t)^l，i＝1,2,…,m，为水库i第l个预报误差模拟值，e_j(t)^l，j＝m+1,m+2,…,m+n，为防洪点区间j第l个预报误差模拟值；

S1235、将步骤1234所得的各集水区间、各预见期的预报误差模拟值叠加至预报过程得到考虑预报不确定性关联的入库洪水系列和防洪区间来水系列分别为：

QR_i(t)^l＝QR_obs,i(t)+e_i(t)^l i＝1,...,m；t＝1,...,T；l＝1,...,L (10)；

QL_j(t)^l＝QL_obs,j(t)+e_j(t)^l j＝1,...,n；t＝1,...,T；l＝1,...,L (11)；

其中，QR_i(t)^l为第l组模拟的水库i在t时刻的入库流量；QL_j(t)^l为第l组模拟的防洪区间j在t时刻的来水流量。

进一步的，步骤S1中水库泄洪的水力不确定性，水位库容关系不确定性和洪水传播演进不确定性描述具体为：

S131、引入泄流能力随机误差因子ε_1,i，则实际泄流能力表示为：

q_max,i(t)＝f_Zq(Z_i(t))·(1+ε_1,i) (12)；

其中，q_max,i(t)为水库i在t时刻的泄流能力；Z_i(t)为水库i在t时刻的水位；f_Zq(·)为水库泄流能力曲线；正态分布是水力不确定性的最佳概率分布模型，即

σ_1,i为水库i的误差分布标准差；

S132、引入水位库容随机误差因子ε_2,i，则实际水位库容关系表示为：

V_i(t)＝f_ZV(Z_i(t))·(1+ε_2,i) (13)；

其中，V_i(t)为水库i在t时刻的库蓄水量；f_ZV(·)为水库水位库容曲线；ε_2,i同样假定为服从正态分布，

σ_2,i为水库i的误差分布标准差；

S133、洪水传播演进采用马斯京根洪水演进模型，引入洪水传播随机误差因子ε_3,i，则实测下游流量表示为：

qc_i(t)＝r(q_i(t))·(1+ε_3,i) (14)；

其中，qc_i(t)为t时刻水库i的出流在下游防洪点或下游水库处的流量响应；q_i(t)为水库i在t时刻的出库；r(·)为马斯京根法计算模型；ε_3,i也假定为服从正态分布，

σ_3,i为水库i的误差分布标准差；

S134、基于ε_1,i、ε_2,i和ε_3,i各正态分布，分别生成L组各库泄流能力、水位库容关系和洪水传播误差为：

进一步的，步骤S2中建立的水库群多目标防洪调度模型的第l个随机情景下的目标函数和约束条件为：

S21、目标函数：

上游防洪目标以一次调洪中占用的防洪库容最小为原则，表示为：

其中，F₁为一次调洪中占用的防洪库容，

为第l个随机情景下水库i在t时刻的水位；

下游防洪目标以各水库总最大泄流最小为原则，表示为：

其中，F₂为各水库最大泄流，

为第l个随机情景下水库i在t时刻的出库流量；

S22、约束条件：

1)水量平衡约束；

其中，

分别为第l个随机情景下t时段初和时段末水库i的库蓄量；

为与水库i有直接水力联系的上游水库集合；

为第k个上游水库出库流量在水库i处的响应流量；

和

分别为t时段水库i的入库和出库流量；△t为调度时段长；

2)水库水位约束；

其中，Z_i(t)、

分别为水库i在t时刻的水位上下限值；

3)出库流量约束；

其中，q_i(t)、

分别为水库i在t时刻的出库流量上下限值；

4)泄流能力约束；

5)出库流量变幅约束；

其中，Δq_i为水库i的允许最大泄流量变幅；

6)初始、边界条件；

其中，

为第l个随机情景下初始时刻的水库i起调水位，即t＝0时刻的水库i起调水位；

为水库i的第l个起调水位模拟值；

为第l个随机情景下调度期末的水库i期末水位，即t＝T时刻的水库i期末水位；Ze_i为水库i的目标期末水位，设为设计汛限水位，根据预报信息进行洪水资源利用，设为汛限水位动态控制上限值。

进一步的，步骤S2中水库群多目标防洪调度模型求解与随机模拟方法为：

采用NSGA-II智能算法求解，将各水库出库流量作为决策变量，采用水库防洪调度中常用的分段试算法进行约束处理，种群数设为N，即非劣解数目；将步骤S1生成的L组随机情景作为模型输入，分别得到相应于L组随机情景的Pareto前沿。

进一步的，步骤S3具体为：

水库群系统防洪风险包括系统内各水库的防洪风险、各下游防洪点的防洪风险和系统总风险，分别定义如下：

1)水库防洪风险：

其中，P_Ri为水库i的防洪风险；Z_max,i为水库i的最高水位；

为水库i在第l组模拟中的最高水位；Z_s,i为水库i的安全水位，根据实时防洪要求设定；

2)下游防洪风险：

其中，P_Qi为水库i的下游防洪风险；qc_max,i为水库i的下游最大流量；

为第l组模拟中水库i的下游最大流量；qc_s,i为水库i的下游安全流量；

3)系统总风险：

系统总风险P_S应兼顾水库和下游的防洪风险，定义为水库最大风险与下游最大风险之和：

P_S＝max(P_Ri)+max(P_Qi) i＝1,2,...,m (26)；

针对N个非劣解方案，可分别求得相应的N组风险值，代表不同目标偏好下的风险，用来分析防洪风险在上下游的变化规律，辅助决策。

有益效果：与现有技术相比，本发明建立的考虑多不确定性、多目标的水库群防洪调度风险分析模型，耦合了洪水资源利用所导致的起调水位不确定性风险源，更符合我国汛期水库调度的实际情况；综合考虑了水库和下游防洪目标，可进行防洪风险在上下游的竞争转化关系分析，剖析水库群系统防洪风险的内在演变规律，为防洪调度策略制定和风险管理提供理论参考。

附图说明

图1是本发明方法流程图；

图2是NSGA-II算法流程图；

图3是分段试算法具体实施步骤图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

如图1所示，本发明的考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，包括以下步骤：

S1、水库群防洪调度不确定性描述；

考虑水库起调水位不确定性，洪水预报的水文不确定性，水库泄洪的水力不确定性，水位库容关系不确定性和洪水传播演进不确定性。

其中，起调水位不确定性描述中，水库间起调水位误差考虑其空间相关性，洪水预报误差考虑其时程-空间相关性，根据历史预报数据和实测数据，基于预报预泄法和水库聚合分解得到起调水位经验样本，以Copula函数建立各水库起调水位的联合分布函数并进行随机模样本生成拟。

洪水预报不确定性描述中，采用t-Copula建立不同集水区间和不同预报时长的预报误差联合分布，并基于Copula模型进行随机模拟。

水库泄洪的水力不确定性，水位库容关系不确定性和洪水传播演进不确定性均以正态分布的随机误差描述。

以下分三部分描述上述不确定性随机模拟具体步骤：

S11、起调水位不确定性：

S111、基于预报预泄法和水库聚合分解获取各水库起调水位的经验样本。预报预泄是指在预报无雨时，无雨期内以水库兴利用水做预泄流量(兴利预泄)，下一场洪水起涨的有效预见期内以下游安全泄量做预泄流量(防洪预泄)，以计算洪水退水段可超蓄洪水资源量，公式表示如下：

其中，ΔW为超蓄水量上限；q_m为水库兴利用水流量；T_d为预报无雨时长；T_f为洪水预报有效预见期；α为安全泄量折扣系数；q_s为下游河道安全泄量；Q_fl(t)为预报流量过程。当前时刻记为0。

记实际所需超蓄水量为ΔW_a，则在以预报信息计算得到的ΔW超蓄的前提下，以实际情况预泄ΔW_a，则起调水位Z₀为：

Z₀＝f(V_x+ΔW-ΔW_a) (2)；

其中，V_x为设计汛限水位对应的水库蓄水量；f(·)为水位库容关系。

对于包含m个水库的水库群系统，采用聚合分解法计算各水库的起调水位样本值：

其中S为样本数量。

S112、以Copula函数建立m个水库的起调水位Z_0,1,Z_0,2,...,Z_0,m的联合分布函数以进行随机模拟。具体实施如下：

S1121、根据Z_0,1,Z_0,2,...,Z_0,m的经验样本值，采用非参数法计算其边缘分布函数值；

非参数法不需事先假设总体的分布形式，本发明采用非参数核密度估计法计算各水库起调水位的边缘密度函数f_i(Z_0,i)的核密度估计：

其中，

为水库i起调水位的第s个样本值，K()为核函数，h为窗宽。

采用R语言的kernel函数进行核密度及核分布估计值计算。

S1122、采用最大似然估计拟合备选的各Copula函数模型的参数，计算工具为R语言copula包中的fitCopula函数。备选Copula函数如下：

(1)Gaussian Copula函数

其中，Φ(·)为标准正态分布函数；ρ为N阶对称正定矩阵；W＝[w₁,w₂,…,w_N]为积分变量矢量。

(2)t-Copula函数

其中，T_υ(·)为自由度为υ的标准t分布函数；ρ为N阶对称正定矩阵；W＝[w₁,w₂,…,w_N]为积分变量矢量。

(3)Gumbel Copula函数

(4)Clayton Copula函数

(5)Frank Copula函数

S1123、计算拟合优度统计量赤池信息准则(Akaike Information Criterion，AIC)、均方根误差(The root mean square error，RMSE)和KS检验统计量(D_max)(见下式)，以进行Copula模型优选；

AIC＝-2ln(L)+2v (9)；

D_max＝max|P_e,s-P_s| (11)；

式中，L为似然函数，v为分布参数个数，P_e和P分别为经验Copula函数和理论Copula函数，P_e,s和P_s分别为第s个样本的经验Copula函数和理论Copula函数值。

经验Copula函数定义如下：

式中，F_i(·)为各库起调水位的经验分布函数，I_[·]为示性函数，当

时，

否则

综合选取AIC、RMSE和D_max统计值最小的Copula模型为最优Copula模型。

S1124、采用R语言copula包中的rCopula函数，根据最优的Copula模型生成服从该联合分布的m个具有相关关系的、服从[0,1]均匀分布的随机向量θ₁,θ₂,…,θ_m，向量长度为L(即模拟组数)；

S1125、根据所得的随机向量样本值，对各变量边缘分布逆向计算求得各水库的起调水位模拟值

为水库i第l个起调水位模拟值。

S12、洪水预报的水文不确定性：

S121、定义第i个水库、第j个防洪点区间在t时刻的预报误差e_i(t)、e_j(t)分别为预报值与实测值之差：

其中，e_i(t)为第i个水库在t时刻的预报误差，e_j(t)为第j个防洪点区间在t时刻的预报误差，QR_for,i(t)和QR_obs,i(t)分别为第i个水库所辖集水区间的预报来水和实际来水；QL_for,j(t)和QL_obs,j(t)分别为第j个防洪点区间的预报来水和实际来水；T为预报系统发布的预报总时长。

S122、对于包含m个水库、n个下游防洪点区间的水库群系统，预报误差表示为以下(m+n)*T的矩阵：

S123、椭圆型Copula可根据变量间的相关系数矩阵建立，在高维分布模型中具有一定的优势。因此，采用t-Copula建立时空关联的预报误差联合分布，并基于Copula模型进行随机模拟，具体步骤如下：

S1231、根据历史预报和实测数据统计各集水区间的洪水预报误差，得到误差序列样本，同样采用非参数法建立其边缘分布；

S1232、采用最大似然法估计t-Copula函数的参数；

S1233、生成服从[0,1]均匀分布的各集水区间、各预见期的误差分布函数值向量θ₁(1),…,θ₁(T)；θ₂(1),…,θ₂(T)；…；θ_m+n(1),…,θ_m+n(T)，向量长度为L(即模拟组数)；

S1234、根据步骤1233中的分布函数序列(即各集水区间、各预见期的误差分布函数值向量)和各变量的边缘分布，逆向计算得到各集水区间、各预见期的预报误差模拟值e₁(1)^l,…,e₁(T)^l；e₂(1)^l,…,e₂(T)^l；…；e_m+n(1)^l,…,e_m+n(T)^l；l＝1,2,…,L，e_i(t)^l(i＝1,2,…,m)为水库i第l个预报误差模拟值，e_j(t)^l(j＝m+1,m+2,…,m+n)为防洪点区间j第l个预报误差模拟值；

S1235、将步骤1234所得的误差序列(即各集水区间、各预见期的预报误差模拟值)叠加至预报过程得到考虑预报不确定性关联的入库洪水系列和防洪区间来水系列分别为：

QR_i(t)^l＝QR_obs,i(t)+e_i(t)^l i＝1,...,m；t＝1,...,T；l＝1,...,L (15)；

QL_j(t)^l＝QL_obs,j(t)+e_j(t)^l j＝1,...,n；t＝1,...,T；l＝1,...,L (16)；

其中，QR_i(t)^l为第l组模拟的水库i在t时刻的入库流量；QL_j(t)^l为第l组模拟的防洪区间j在t时刻的来水流量。

S13、水库泄洪的水力不确定性，水位库容关系不确定性和洪水传播演进不确定性：

S131、引入泄流能力随机误差因子ε_1,i，则实际泄流能力表示为：

q_max,i(t)＝f_Zq(Z_i(t))·(1+ε_1,i) (17)；

其中，q_max,i(t)为水库i在t时刻的泄流能力；Z_i(t)为水库i在t时刻的水位；f_Zq(·)为水库泄流能力曲线。研究表明，正态分布是水力不确定性的最佳概率分布模型，即

σ_1,i为水库i的误差分布标准差，可由水库实际泄流资料确定。

S132、引入水位库容随机误差因子ε_2,i，则实际水位库容关系表示为：

V_i(t)＝f_ZV(Z_i(t))·(1+ε_2,i) (18)；

其中，V_i(t)为水库i在t时刻的库蓄水量；f_ZV(·)为水库水位库容曲线。ε_2,i同样假定为服从正态分布，

σ_2,i为水库i的误差分布标准差，可由水库实际蓄水量观测资料确定。

S133、洪水传播演进采用马斯京根洪水演进模型，引入洪水传播随机误差因子ε_3,i，则实测下游流量表示为：

qc_i(t)＝r(q_i(t))·(1+ε_3,i) (19)；

其中，qc_i(t)为t时刻水库i的出流在下游防洪点或下游水库处的流量响应；q_i(t)为水库i在t时刻的出库；r(·)为马斯京根法计算模型。ε_3,i也假定为服从正态分布，

σ_3,i为水库i的误差分布标准差，可由下游流量的实际观测资料确定。

S134、基于上述ε_1,i、ε_2,i和ε_3,i各正态分布，分别生成L组各库泄流能力、水位库容关系和洪水传播误差为：

S2、水库群多目标防洪调度模型建立及随机模拟：调度模型考虑上下游两个目标，上游防洪目标以一次调洪中占用的防洪库容最小为原则，下游防洪目标以各水库最大泄流最小为原则；调度模型约束条件包括：水量平衡约束、水库水位约束、出库流量约束、泄流能力约束、出库流量变幅约束和初始、边界条件；

第l个随机情景下的目标函数和约束条件如下：

S21、目标函数：

上游防洪目标以一次调洪中占用的防洪库容最小为原则，表示为：

其中，F₁为一次调洪中占用的防洪库容，

为第l个随机情景下水库i在t时刻的水位。

下游防洪目标以各水库最大泄流最小为原则，表示为：

其中，F₂为各水库最大泄流，

为第l个随机情景下水库i在t时刻的出库流量。

S22、约束条件：

1)水量平衡约束；

其中，

分别为第l个随机情景下t时段初和时段末水库i的库蓄量；

为与水库i有直接水力联系的上游水库集合；

为第k个上游水库出库流量在水库i处的响应流量；

和

分别为t时段水库i的入库和出库流量；△t为调度时段长。

2)水库水位约束；

其中，Z_i(t)、

分别为水库i在t时刻的水位上下限值。

3)出库流量约束；

其中，q_i(t)、

分别为水库i在t时刻的出库流量上下限值。

4)泄流能力约束；

5)出库流量变幅约束；

其中，Δq_i为水库i的允许最大泄流量变幅。

6)初始、边界条件；

式中，

为第l个随机情景下初始时刻(t＝0)的水库i起调水位；

为水库i的第l个起调水位模拟值；

为第l个随机情景下调度期末(t＝T)的水库i期末水位，Ze_i为水库i的目标期末水位，一般设为设计汛限水位，可根据预报信息进行洪水资源利用，设为汛限水位动态控制上限值。

S23、模型求解与随机模拟：

本模型采用NSGA-II智能算法求解，即带有精英保留策略的快速非支配多目标优化算法。该算法是求解多目标优化问题的一个经典算法，其有效性已得到多次证明。将各水库出库流量作为决策变量，采用水库防洪调度中常用的分段试算法进行约束处理，种群数设为N(即非劣解数目)。将步骤S1生成的L组随机情景作为模型输入进行循环求解，分别得到相应于L组随机情景的Pareto前沿。

NSGA-II算法流程如图2所示，其具体计算流程如下：

(1)初始化种群，设置种群数为N，最大进化代数为maxGen；

(2)Gen＝1，采用分段试算法对初始种群进行约束处理和适应度值计算，并进行快速非支配排序和拥挤度计算，得到非支配集P(1)；

(3)对种群P(Gen)进行选择、交叉和变异操作，得到N个后代S(Gen)，并计算其适应度值；

(4)合并P(Gen)和S(Gen)，得到规模为2N的种群PS(Gen)，对PS(Gen)进行快速非支配排序和拥挤度计算，根据拥挤度比较算子选择较好的N个个体作为新的子代P(Gen+1)；

(5)若Gen<maxGen，令Gen＝Gen+1，返回步骤(3)；否则结束计算。

步骤(2)中分段试算法流程如图3所示，对于第l组模拟中水库i的具体计算步骤如下：

(21)读取随机入流情景

和水库起调水位

等初始条件，读取初始化的随机出库流量

(22)根据

和

按水量平衡调节计算，并逐时段检验泄流能力约束和出库允许变幅约束，若满足转(23)，否则：

①若

则令

②若

则令

③重新调节计算。

(23)检查最高水位约束，若满足转(24)，否则：

①计算满足最高水位约束所需调整的水量

t_m为计算库最高水位出现的时刻。按下式计算分段调整系数：

当t∈[0,t_m]时，

当t∈[t₀,t_m]时，

②调整下泄量

其中，T_m、T_e分别是[0,t_m]区间和[t_m,T]区间中未因泄流能力约束和出库变幅约束而调整的时段数。采用等比例缩放调整出库流量，目的是使调整后的出库流量过程尽可能保持与调整前的出库流量过程的基本形状相似。

③转(22)重新调节计算。

(24)检查期末水位约束，若满足则整理计算结果，结束计算。否则：

转(22)重新调节计算。

S3、水库群防洪调度风险分析与计算；水库群防洪调度风险分别从各水库的防洪风险、各下游防洪点的防洪风险和系统总风险进行定义，其中，水库和防洪点风险定义为超安全阈值的累积概率，系统风险定义为水库最大风险与下游最大风险之和。

水库群系统防洪风险包括系统内各水库的防洪风险、各下游防洪点的防洪风险和系统总风险，分别定义如下：

1)水库防洪风险：

其中，P_Ri为水库i的防洪风险；Z_max,i为水库i的最高水位；

为水库i在第l组模拟中的最高水位；Z_s,i为水库i的安全水位，可根据实时防洪要求设定。

2)下游防洪风险：

其中，P_Qi为水库i的下游防洪风险；qc_max,i为水库i的下游最大流量；

为第l组模拟中水库i的下游最大流量；qc_s,i为水库i的下游安全流量。

3)系统总风险：

系统总风险P_S应兼顾水库和下游的防洪风险，定义为水库最大风险与下游最大风险之和：

P_S＝max(P_Ri)+max(P_Qi) i＝1,2,...,m (31)；

针对N个非劣解方案，可分别求得相应的N组风险值，代表不同目标偏好下的风险，可用来分析防洪风险在上下游的变化规律，辅助决策。

本发明所提出的考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，耦合了起调水位、洪水预报误差、水位库容关系误差、泄流能力曲线误差和洪水传播误差多维不确定性，考虑水库间起调水位误差的相关关系和洪水预报误差的时程-空间相关性，建立相应的Copula模型以进行随机模拟的样本生成，建立水库群防洪调度多目标优化模型，分析在不同目标之间的风险变化关系，为水库群系统防洪调度风险调控提供了理论参考。

Claims (7)

1.一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，其特征在于，包括以下步骤：

S1、水库群防洪调度不确定性描述；

防洪调度不确定性包括水库起调水位不确定性，洪水预报的水文不确定性，水库泄洪的水力不确定性，水位库容关系不确定性和洪水传播演进不确定性；

S2、水库群多目标防洪调度模型建立及随机模拟；

调度模型考虑上下游两个目标，上游防洪目标以一次调洪中占用的防洪库容最小为原则，下游防洪目标以各水库最大泄流最小为原则；

调度模型约束条件包括：水量平衡约束、水库水位约束、出库流量约束、泄流能力约束、出库流量变幅约束和初始、边界条件；

采用NSGA-II智能算法求解水库群多目标防洪调度模型，将各水库出库流量作为决策变量，采用水库防洪调度中常用的分段试算法进行约束处理，种群数为非劣解数目；

S3、水库群防洪调度风险分析与计算；

水库群防洪调度风险分别从各水库的防洪风险、各下游防洪点的防洪风险和系统总风险进行定义，其中，水库和防洪点风险定义为超安全阈值的累积概率，系统风险定义为水库最大风险与下游最大风险之和。

2.根据权利要求1所述的一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，其特征在于，步骤S1中水库起调水位不确定性描述具体为：

S111、基于预报预泄法和水库聚合分解获取各水库起调水位的经验样本；预报预泄是指在预报无雨时，无雨期内以水库兴利用水做预泄流量，下一场洪水起涨的有效预见期内以下游安全泄量做预泄流量，以计算洪水退水段可超蓄洪水资源量，公式表示如下：

其中，ΔW为超蓄水量上限；q_m为水库兴利用水流量；T_d为预报无雨时长；T_f为洪水预报有效预见期；α为安全泄量折扣系数；q_s为下游河道安全泄量；Q_fl(t)为预报流量过程；当前时刻记为0；

记实际所需超蓄水量为ΔW_a，则在以预报信息计算得到的ΔW超蓄的前提下，以实际情况预泄ΔW_a，则起调水位Z₀为：

Z₀＝f(V_x+ΔW-ΔW_a) (2)；

其中，V_x为设计汛限水位对应的水库蓄水量；f(·)为水位库容关系；

对于包含m个水库的水库群系统，采用聚合分解法计算各水库的起调水位样本值：

其中S为样本数量；

S112、以Copula函数建立m个水库的起调水位Z_0,1,Z_0,2,...,Z_0,m的联合分布函数以进行随机模拟；具体为：

S1121、根据Z_0,1,Z_0,2,...,Z_0,m的经验样本值，采用非参数核密度估计法计算各水库起调水位的边缘密度函数f_i(Z_0,i)的核密度估计：

其中，

为水库i起调水位的第s个样本值，K(·)为核函数，h为窗宽；

采用R语言的kernel函数进行核密度及核分布估计值计算；

S1122、采用最大似然估计拟合备选的各Copula函数模型的参数；

S1123、计算拟合优度统计量赤池信息准则、均方根误差和KS检验统计量，以进行Copula模型优选；

赤池信息准则计算公式为：

AIC＝-2ln(L)+2v (4)；

均方根误差计算公式为：

KS检验统计量计算公式为：

D_max＝max|P_e,s-P_s| (6)；

其中，L为似然函数，v为分布参数个数，P_e和P分别为经验Copula函数和理论Copula函数，P_e,s和P_s分别为第s个样本的经验Copula函数和理论Copula函数值；

经验Copula函数定义如下：

其中，F_i(·)为各库起调水位的经验分布函数，I_[·]为示性函数，当

时，

否则

综合选取AIC、RMSE和D_max统计值最小的Copula模型为最优Copula模型；

S1124、采用R语言copula包中的rCopula函数，根据最优的Copula模型生成服从该联合分布的m个具有相关关系的、服从[0,1]均匀分布的随机向量θ₁,θ₂,…,θ_m，向量长度为L，即模拟组数；

S1125、根据所得的随机向量θ₁,θ₂,…,θ_m，对各变量边缘分布逆向计算求得各水库的起调水位模拟值

为水库i第l个起调水位模拟值，L为模拟组数。

3.根据权利要求1所述的一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，其特征在于，步骤S1中洪水预报的水文不确定性描述具体为：

S121、定义第i个水库、第j个防洪点区间在t时刻的预报误差e_i(t)、e_j(t)分别为预报值与实测值之差：

其中，e_i(t)为第i个水库在t时刻的预报误差，e_j(t)为第j个防洪点区间在t时刻的预报误差，QR_for,i(t)和QR_obs,i(t)分别为第i个水库所辖集水区间的预报来水和实际来水；QL_for,j(t)和QL_obs,j(t)分别为第j个防洪点区间的预报来水和实际来水；T为预报系统发布的预报总时长；

S122、对于包含m个水库、n个下游防洪点区间的水库群系统，预报误差表示为以下(m+n)*T的矩阵：

S123、采用t-Copula建立时空关联的预报误差联合分布，并基于Copula模型进行随机模拟，具体为：

S1231、根据历史预报和实测数据统计各集水区间的洪水预报误差，得到误差序列样本，同样采用非参数法建立其边缘分布；

S1232、采用最大似然法估计t-Copula函数的参数；

S1233、生成服从[0,1]均匀分布的各集水区间、各预见期的误差分布函数值向量θ₁(1),…,θ₁(T)；θ₂(1),…,θ₂(T)；…；θ_m+n(1),…,θ_m+n(T)，向量长度为L，即模拟组数；

S1234、根据上步的分布函数值向量和各变量的边缘分布，逆向计算得到各集水区间、各预见期的预报误差模拟值e₁(1)^l,…,e₁(T)^l；e₂(1)^l,…,e₂(T)^l；…；e_m+n(1)^l,…,e_m+n(T)^l；l＝1,2,…,L，e_i(t)^l，i＝1,2,…,m，为水库i第l个预报误差模拟值，e_j(t)^l，j＝m+1,m+2,…,m+n，为防洪点区间j第l个预报误差模拟值；

S1235、将步骤1234所得的各集水区间、各预见期的预报误差模拟值叠加至预报过程得到考虑预报不确定性关联的入库洪水系列和防洪区间来水系列分别为：

QR_i(t)^l＝QR_obs,i(t)+e_i(t)^l i＝1,...,m；t＝1,...,T；l＝1,...,L (10)；

QL_j(t)^l＝QL_obs,j(t)+e_j(t)^l j＝1,...,n；t＝1,...,T；l＝1,...,L (11)；

其中，QR_i(t)^l为第l组模拟的水库i在t时刻的入库流量；QL_j(t)^l为第l组模拟的防洪区间j在t时刻的来水流量。

4.根据权利要求1所述的一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，其特征在于，步骤S1中水库泄洪的水力不确定性，水位库容关系不确定性和洪水传播演进不确定性描述具体为：

S131、引入泄流能力随机误差因子ε_1,i，则实际泄流能力表示为：

q_max,i(t)＝f_Zq(Z_i(t))·(1+ε_1,i) (12)；

其中，q_max,i(t)为水库i在t时刻的泄流能力；Z_i(t)为水库i在t时刻的水位；f_Zq(·)为水库泄流能力曲线；正态分布是水力不确定性的最佳概率分布模型，即

σ_1,i为水库i的误差分布标准差；

S132、引入水位库容随机误差因子ε_2,i，则实际水位库容关系表示为：

V_i(t)＝f_ZV(Z_i(t))·(1+ε_2,i) (13)；

其中，V_i(t)为水库i在t时刻的库蓄水量；f_ZV(·)为水库水位库容曲线；ε_2,i同样假定为服从正态分布，

σ_2,i为水库i的误差分布标准差；

S133、洪水传播演进采用马斯京根洪水演进模型，引入洪水传播随机误差因子ε_3,i，则实测下游流量表示为：

qc_i(t)＝r(q_i(t))·(1+ε_3,i) (14)；

其中，qc_i(t)为t时刻水库i的出流在下游防洪点或下游水库处的流量响应；q_i(t)为水库i在t时刻的出库；r(·)为马斯京根法计算模型；ε_3,i也假定为服从正态分布，

σ_3,i为水库i的误差分布标准差；

S134、基于ε_1,i、ε_2,i和ε_3,i各正态分布，分别生成L组各库泄流能力、水位库容关系和洪水传播误差为：

l＝1,2,…,L。

5.根据权利要求1所述的一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，其特征在于，步骤S2中建立的水库群多目标防洪调度模型的第l个随机情景下的目标函数和约束条件为：

S21、目标函数：

上游防洪目标以一次调洪中占用的防洪库容最小为原则，表示为：

其中，F₁为一次调洪中占用的防洪库容，

为第l个随机情景下水库i在t时刻的水位；

下游防洪目标以各水库总最大泄流最小为原则，表示为：

其中，F₂为各水库最大泄流，

为第l个随机情景下水库i在t时刻的出库流量；

S22、约束条件：

1)水量平衡约束；

其中，

分别为第l个随机情景下t时段初和时段末水库i的库蓄量；

为与水库i有直接水力联系的上游水库集合；

为第k个上游水库出库流量在水库i处的响应流量；

和

分别为t时段水库i的入库和出库流量；△t为调度时段长；

2)水库水位约束；

其中，Z_i(t)、

分别为水库i在t时刻的水位上下限值；

3)出库流量约束；

其中，q_i(t)、

分别为水库i在t时刻的出库流量上下限值；

4)泄流能力约束；

5)出库流量变幅约束；

其中，Δq_i为水库i的允许最大泄流量变幅；

6)初始、边界条件；

其中，

为第l个随机情景下初始时刻的水库i起调水位，即t＝0时刻的水库i起调水位；

为水库i的第l个起调水位模拟值；

为第l个随机情景下调度期末的水库i期末水位，即t＝T时刻的水库i期末水位；Ze_i为水库i的目标期末水位，设为设计汛限水位，根据预报信息进行洪水资源利用，设为汛限水位动态控制上限值。

6.根据权利要求1所述的一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，其特征在于，步骤S2中水库群多目标防洪调度模型求解与随机模拟方法为：

采用NSGA-II智能算法求解，将各水库出库流量作为决策变量，采用水库防洪调度中常用的分段试算法进行约束处理，种群数设为N，即非劣解数目；将步骤S1生成的L组随机情景作为模型输入，分别得到相应于L组随机情景的Pareto前沿。

7.根据权利要求1所述的一种考虑时空关联多维不确定性的水库群多目标防洪调度风险分析方法，其特征在于，步骤S3具体为：

水库群系统防洪风险包括系统内各水库的防洪风险、各下游防洪点的防洪风险和系统总风险，分别定义如下：

1)水库防洪风险：

其中，P_Ri为水库i的防洪风险；Z_max,i为水库i的最高水位；

为水库i在第l组模拟中的最高水位；Z_s,i为水库i的安全水位，根据实时防洪要求设定；

2)下游防洪风险：

其中，P_Qi为水库i的下游防洪风险；qc_max,i为水库i的下游最大流量；

为第l组模拟中水库i的下游最大流量；qc_s,i为水库i的下游安全流量；

3)系统总风险：

系统总风险P_S应兼顾水库和下游的防洪风险，定义为水库最大风险与下游最大风险之和：

P_S＝max(P_Ri)+max(P_Qi) i＝1,2,...,m (26)；

针对N个非劣解方案，可分别求得相应的N组风险值，代表不同目标偏好下的风险，用来分析防洪风险在上下游的变化规律，辅助决策。

Images