基于蒙特卡洛法的水库汛限水位动态控制风险分析

建剑波 曾智珍 卢金阁 孙志伟 何芳婵

（1河南省河口村水库管理局；2河海大学水利水电工程学院；3南京水利科学研究院；4河南省水利科学研究院）

1 引言

中国现阶段绝大部分水库多以静态的方式对汛限水位进行控制，该方式虽然能够保证水库的防洪安全，但在汛前或汛期常常要进行大量弃水以降低至汛限水位，从而造成洪水资源大量浪费，不利于水库的兴利调度。中国是一个水资源缺乏和水资源时空分布不均的国家，近年来随着经济的快速发展和城镇化建设，水资源的供需矛盾越来越突出。为充分利用水资源，满足可持续发展的需求，尽可能利用现有的降雨预报和洪水预报技术，优化水库调度运行，实施水库汛限水位的动态控制，对提高水库水资源的利用率，发挥水库综合效益，以及水资源可持续发展战略都具有重要的意义。

文章以基于预报预泄的水库洪水动态水位控制为研究对象，结合河口村水库实例分析，对水文气象预报等不确定性信息进行风险估计，为决策者提供理论依据和决策信息。

2 汛限水位动态控制风险定义及风险源分析

2.1 汛限水位动态控制风险定义

水库汛限水位动态控制风险指的是由于水库调度过程中自然、人为等不确定因素的影响，水库系统在实施汛限水位动态控制前后的风险变化，主要从风险率变化和后果变化两方面得以反映。

已建成的水库，工程规模和功能是确定的。假设面临调度期内调洪最高水位Z为荷载，水库运行阶段允许最高水位Zm为承载能力，定义“Z＞Zm”为水库洪水动态水位控制导致的水库安全风险事件，则风险率可表达为：

式中字母含义同下。不考虑调度规则的变化，当起调水位Z0超过原设计汛限水位Z汛时，水库最高水位Zmax可由式（2）～（4）确定：

2.2 汛限水位动态控制风险源分析

基于预报预泄的水库汛期水位动态控制的实质，是利用洪水尾水实现超原设计汛限水位蓄水，关键问题是分析实时增蓄水量的大小。增蓄水量在上一场洪水退水段形成，其消落有兴利预泄与防洪预泄两种方式。如果增蓄后无雨期很长，则增蓄的水量可以通过兴利预泄消落，在兴利预泄来不及时，可通过防洪预泄消落。根据基于预报预泄的水库洪水动态水位风险控制方法的原理可知，导致水库水位不能回落到原设计汛限水位的原因主要有以下方面：

一是入库流量过程预报结果偏小，实际来水比预报来水大；二是下一场洪水发生的时间提前，水库丧失了部分或全部兴利预泄时间；三是下一场洪水入库时间比预计的快，防洪预泄水量比预计的减少。

可见预报入库水量的误差、预泄采用的连续无雨天数，以及预泄采用的洪水预报预见期是导致水库水位可能高于的三个不确定性（风险）因子，这三个因子既可能单独，也可能组合发生作用。

3 基于蒙特卡洛模拟的风险分析方法

汛限水位动态控制的各风险变量之间存在着比较复杂的影响机制，而使用蒙特卡洛模拟方法（MC法）来解决复杂系统的风险分析是一种行之有效的方法。MC法在水库洪水动态水位控制风险分析中按照以下步骤进行：①根据实验数据或者经验判断确定水库洪水动态水位控制中主要不确定性因素的概率分布。②随机生成符合相应分布的各主要不确定性因素的数值。③将不确定性因素叠加到水库洪水动态水位控制的过程中，调洪演算得到防洪目标（调洪最高水位）值。④重复步骤②～③，可得到目标样本，运用统计学方法分析得到在指定条件下防洪目标破坏的频率，作为风险计算的基本数据。

3.1 风险因子的随机模拟

3.1.1 预报入库水量的随机数生成

式中：εt为服从水量预报相对误差系列分布的正态随机数；ηt为标准化正态随机数；为水量预报相对误差系列均值；σ为水量预报相对误差系列均方差。

3.1.2 洪水预报预见期的随机数生成

因洪水预报预见期亦服从正态分布，故其随机数生成方法同3.1.1。

3.1.3 连续无雨天数的随机数生成

连续无雨天数服从负指数分布Exp（λ，μ），即

作如下变换

当y 为［0，1］均匀分布随机数时，d 即为服从式（7）的负指数分布。

因此，要生成符合负指数分布的随机数，首先需要生成符合U［0，1］均匀分布的随机数系列{yt}，将其代入式（7）中即可得到服从负指数分布Exp（λ，μ）的随机数系列{dt}。

3.2 基于蒙特卡洛模拟的风险率计算

第一，在确定性预报洪水过程Q（t）的基础上，考虑预报误差构建不同的入库流量过程。

计算洪水退水阶段及下一场洪水起涨阶段的洪量：

式中：t1、tm分别为预报洪水过程的起点与终点。

考虑洪量预报相对误差，计算修正的入库洪量W'：

利用同倍比法放大，计算考虑预报误差后修正的流量过程：

第二，用蒙特卡洛模拟法生成服从正态分布的洪水预报预见期系列{εt}，以及服从负指数分布的连续无雨天数系列{dt}。

第三，将考虑预报误差的流量过程系列、洪水预报预见期系列以及连续无雨天数系列随机组合，得到组合随机数系列。

第四，将组合随机数中的流量过程、洪水预报预见期以及连续无雨天数代入式（12）～（13）中。

计算得到从t0到t2水库的蓄量变化△W'。已知t0时刻水库完成洪水动态水位控制的预蓄过程且水位达到Zms，根据式（14）～（15）计算得到下一场洪水起涨时刻即t2时刻的水位Z0。然后，根据式（2）～（4），得到最高水位Zmax。

式中：V0为t2时刻的水库蓄量。

第五，考虑不同的组合随机数，重复步骤（4），得到最高水位Zmax系列。

第六，假设决策者所能接受的安全水位指标为Zd，则相应的风险率为：

式中：m 为随机试验总次数；mf为调洪最高水位大于安全水位（防洪控制最高水位）的次数。

4 实例分析

单独利用预蓄预泄确定河口村水库汛期水位动态控制的上限值为246.00 m。单独利用实时预报调度技术确定常遇洪水（20 a一遇以下）的汛期水位控制上限值为250.59 m，综合考虑后将河口村水库的主汛期水位动态控制的上限值取为246.00 m。

表1 洪水预报样本选取表

选取河口村水库1980-2003年共8场洪水资料作为分析对象，选择预见期T=4 h，分析预见期内实测洪量与预报洪量的误差。

8场洪水共获得有效样本434个，经统计4 h预见期洪量相对误差分布见表2。

表2 预见期4 h预报洪量相对误差分布统计表

假定洪量预报误差符合正态分布，求得4 h预见期洪量预报误差服从正态分布N（65.71，75.232）。同时，参考河口村水文资料可知，洪水预报预见期系列服从正态分布N（9.80，3.212），连续无雨天数系列（P≤3 mm）服从负指数分布Exp（0.19，1）。

考虑洪量预报误差、连续无雨天数以及洪水预报预见期三因子随机组合，采用蒙特卡洛方法生成10 000组随机数，由于19830907 号洪水滚动预报长度不够，故舍弃该场洪水，选取19800728、19920811、20010726 三场洪水进行风险分析，根据基于蒙特卡洛模拟的风险分析原理，分别得到下一场洪水起涨时刻的水位Z0系列，偏安全考虑，选取每种随机模拟情景下三场洪水中的最高水位作为最终的结果，结果如图1所示。

图1 Z0系列水位分布直方图与累积频率曲线图

在10 000次随机模拟中，下一场洪水起涨时刻的水位Z0超过汛限水位238.00 m的次数为237次，概率为2.37%，图1为Z0超过238.00 m的分布情况。

按照起调水位为238.00 m，遭遇不同频率设计洪水时的调洪最高水位作为特征水位，特征水位取值见表3，求得对应不同频率设计洪水的最高水位Zmax系列，安全水位取防洪高水位285.43 m，统计最高水位超过安全水位的概率，即为风险率，结果如表4所示。

表3 前汛期起调水位238.00 m对应的水库各频率洪水特征值（初步设计）表

表4 基于蒙特卡洛模拟法的风险率计算表

从表4可以看出，在主汛期水位动态控制值246.00 m的方案下，若考虑预报入库水量的误差ε1、预泄采用的连续无雨天数d，以及预泄采用的洪水预报预见期 这三个不确定性因子，该方案未增加500 a一遇及以下量级的洪水的风险，稍增加2 000 a一遇洪水的风险，但增加的概率也仅为1.19×10-3%。

5 总结

阐述了汛限水位动态控制风险分析的基本流程，包括对洪量预报误差、洪水预报预见期及连续无雨天数等风险因子的随机数模拟和基于蒙特卡洛模拟的汛限水位动态控制风险率计算。

以河口村水库为例，采用蒙特卡洛模拟方法分析了水库汛期水位动态控制的风险，得到在汛限水位抬高至246.00 m情况下，水库在下一场洪水起涨时刻超原汛限水位的概率仅为2.37%，该方案不增加500 a一遇及以下量级的洪水风险，稍增加2 000 a 一遇洪水的风险，汛期水位动态控制产生的风险极低。