刁艳芳,王 蒙,2,王 昊,丛方杰,王 刚
(1.山东农业大学水利土木工程学院,山东泰安 271018;2.华北水利水电大学水利学院,郑州 450046;3.烟台大学网络与教育技术中心,山东烟台 264005)
近几年,极端气候所导致的洪水灾害频繁发生,造成了巨大的社会、经济损失。水库是保障防洪安全、提高洪水资源利用率和促进水生态文明建设的重要工程措施,科学合理的水库防洪调度更能有效促进其作用的发挥。目前,水库防洪调度的研究方法大体分为常规方法和系统分析方法两类。常规调度方法是借助水库防洪调度图、抗洪能力图等经验图表实施操作的一种半经验、半理论方法[1,2]。随着系统工程的发展和应用,许多学者尝试采用系统分析方法解决水库调度问题。早期主要是线性规划[3,4]、动态规划[5,6]、逐次优化算法[7]和大系统分解协调原理[8]等;然而面对由水库、湖泊、闸坝、河道及堤防、蓄滞洪区以及各种防洪工程组成的复杂防洪系统,基于生物学、物理学和人工智能的现代启发式智能优化算法迅速发展,其具有全局最优化、鲁棒性强、通用性强及高效性等优点,主要包括遗传算法[9]、蚁群算法[10]、粒子群算法[11]和人工神经网络[12]等。
然而无论是常规方法还是系统分析方法,均存在对历史调度数据信息利用不够的问题,而这些数据既包括了水库流域水文特征与规律,又包括了决策者多年调度和配置经验。在当前水利部门已积累了大量流域水文、地貌以及洪水调度数据的基础上,采用数据挖掘技术解决上述问题是主要手段。随机森林算法是数据挖掘技术的一种,在降雨预测[13]、径流预测[14,15]、水资源评价[16]以及洪灾预测[17]等诸多领域得到了广泛应用,但在水库调度中应用极少。大量研究成果表明,随机森林算法具有不易出现过拟合、泛化误差低及预测精度高等优点,被认为是当前最好的算法之一。因此,本文提出采用随机森林算法拟定水库防洪调度方案,并以山东省田庄水库为例验证了本方法的合理性。
田庄水库位于淮河流域沂河上游,流域面积424 km2,总库容为1.305 7 亿m3,兴利库容为0.684 亿m3,是一座以防洪、灌溉为主,结合水力发电、水面养殖、工业供水等综合利用的大(II)型水库。水库兴建于1958年,1960年建成蓄水,按百年一遇设计,五千年一遇校核,其设计洪水位、校核洪水位、汛限水位、防洪高水位(P=2%)及死水位分别为312.38、315.07、310.64、312.33及293.64 m。水库保护下游沂源县城及南麻、悦庄、鲁村等7 个乡镇30 万人口及10 万亩耕地。下游第一安全泄量为1 000 m3/s,田庄水库控制下泄流量为600 m3/s,相应标准为20年一遇;下游第二安全泄量为2 000 m3/s,田庄水库控制下泄流量为1 000 m3/s,相应标准为50年一遇。
田庄水库流域地处中纬度,受东亚季风和欧亚大陆的影响,属温带季风大陆性气候,四季温差大。流域多年平均降水量为730.6 mm,年际之间变化大,年内分配也不均匀,主要集中在6-9月份,约占多年平均年降水量的74.9%。
根据田庄水库1963-2019年的洪水调度数据,考虑到较大量级的洪水对防洪安全更具重要意义,选取了启用溢洪闸且最大泄流量超过50 m3/s 的18 场洪水过程,整编出洪水过程逐时段的入库流量、库水位、蓄水量及出库流量等数据资料,时段长为1小时,其洪号见表1所示。根据《水文情报预报规范》(GB/T 22482-2008),按照洪水总量W将18 场洪水划分为小(W<3 369 万m3)、中(3 369 万m3≤W<6 470 万m3)、大(6 470 万m3≤W<8 800 万m3)、特大(W≥8 800 万m3)4 个级别,列于表1 中。田庄水库流域内洪水主要由暴雨形成,由于天气系统原因导致暴雨持续时间一般较短;此外,田庄水库流域为纯山丘地区,山高坡陡,汇流历时较短,因此,造成大部分洪水过程为形状比较尖瘦的单峰过程,且洪水持续时间为3日左右。表1 所示的18场洪水过程均为符合上述特性的单峰洪水过程。田庄水库流域内共有6个雨量站,分别为包家庄站、草埠站、大张庄站、田庄水库站、徐家庄站及朱家庄站,整理出与18 场洪水过程相对应的各雨量站逐时段的降雨量资料,鉴于雨量站分布较均匀,故采用算术平均法计算流域平均降水量。
表1 场次洪水及级别表Tab.1 Floods and their grades
2.2.1 随机森林算法
随机森林(Random Forest,RF)算法是Breiman 在2001年提出的一种分类和预测的机器学习算法[18],其以Bagging 集成学习算法和随机空间算法为基础。RF 算法分为RF 分类和RF 回归两类,本文采用后者进行水库调度方案的拟定。与目前其他机器学习算法相比,RF 算法具有如下优点:①较强非线性模拟能力,可以高效处理多变量和大数据量的问题,且无需进行变量的删减和量纲统一处理,同时能够评估所有变量的重要程度;②模型参数少、运算效率高、数据挖掘能力强及预测精度高等;③不易出现过拟合现象,对异常值、缺失值、干扰值的容忍度高,对数据集特征的挖掘具有很好的鲁棒性。
(1)RF 回归算法。RF 回归算法的基本组成单元是回归树{t(x,θi),i= 1,2,…,k},其中x表示样本数据集,k表示回归树的个数,θi表示第i棵回归树的参数向量,算法实现过程如下,具体流程见图1所示。
图1 RF回归算法流程图Fig.1 Flow chart of RF regression algorithm
①通过自助法(Bootstrap)重抽样技术,从原始训练样本集N中有放回地重复随机抽取n个样本生成子训练集k个,每个子训练集构建1 棵回归树,每次未被抽到的样本称为袋外数据(Out-Of-Bag,OOB)。
②设原始训练样本集的变量个数为m,在每棵树的叶节点处等概率随机抽取mt个变量作为备选分枝变量,根据分枝优度准则计算最佳分枝,且在建树过程中mt保持不变。
③设叶节点的最小尺寸为nodesize,每棵树最大限度地生长,不进行干预和剪裁。
④k棵回归树组成随机森林回归模型,将验证样本集代入该模型,以每棵回归树输出结果的算数平均值作为预测结果。
(2)RF 回归算法的参数优化。由RF 回归算法的计算过程可以看出,回归树个数k、节点备选分枝变量mt以及叶节点最小尺寸nodesize是3 个主要参数,其取值影响该算法的效率和精度。已有研究表明,k太小会使训练不够充分,降低RF 的随机性和精度;太大会增加模型的计算量,导致计算效率低下、计算精度不高。mt太小会导致过拟合现象,预测分类的误差变大及预测精度降低;mt太大会使建模效率降低。nodesize对RF 回归算法性能的影响不显著。由此可以看出,RF回归算法的预测精度主要取决于k和mt,其取值采用OOB 均方误差进行评价,计算公式为:
式中:MSEOOB为OOB 均方误差;yi为OOB 中因变量的实际值;为OOB的预测值。
MSEOOB可以用来评估回归树的泛化误差,Breiman 通过大量实验数据证明OOB误差是一种无偏估计[19],不需要使用交叉检验进行估计。理论上,可以通过MSEOOB来选取最佳的k和mt数值组合;但是在具体实践中,往往是设定k值或mt值,之后对另一个参数进行适应性调整。
2.2.2 随机森林算法拟定调度方案的步骤
采用RF回归算法推求场次洪水调度方案的步骤如下:
(1)洪水优化调度数据集的生成。由于历史洪水实测调度未必是最优调度方案,故采用优化算法对历史洪水过程进行优化调度,生成水库洪水优化调度数据集。
(2)水库泄流量主要影响因素的选取及RF 回归算法最优参数的推求。通过查阅已有研究成果[20,21]可知,影响水库泄流量的因素包括降雨量、净雨量、库水位、入库流量及入库水量等,整理出各场次洪水优化调度过程逐时段的累计降雨量、累计净雨量、库水位、入库流量、累计入库水量及泄流量等数据,以泄流量作为决策属性,以影响因素作为条件属性,输入RF 回归算法中,由此计算出最优的k和mt值。同时,RF 回归算法可以计算出每个影响因素对每棵回归树的贡献,取均值后即为影响因素对RF 回归模型的贡献程度,由此可以判断出影响泄流量的主要因素。
(3)生成泄流量回归树。将所有场次洪水划分为训练样本和验证样本,将训练样本洪水优化调度过程逐时段的主要影响因素和泄流量输入RF回归算法生成泄流量回归树。
(4)由RF回归算法拟定泄流量。
①设i=1;
②将第i时段的主要影响因素代入泄流量回归树中,计算泄流量qi,并由水量平衡方程求得zi;
③为保护下游防洪安全,判断qi是否不大于该场洪水由水库调度规则调洪后的最大泄流量q*。如果是,转入第④步;否则,qi=q*;
④为确保水库泄流能力约束,判断qi是否不大于水位zi时的水库泄流能力q(zi)。如果是,转入第⑤步;否则,qi=q(zi);
⑤令i=i+1,如果i≤T,T为洪水总时段数,转入第②步;否则,调洪结束。
其流程见图2所示。
图2 RF回归算法拟定调度方案流程图Fig.2 Flow chart of operation scheme maked by RF regression algorithm
田庄水库为大(II)型水库,调节性能较高且承担下游防洪任务,故洪量在调洪过程中起主要控制作用。结合田庄水库特性和已有研究成果,选取累计净雨量、累计入库水量、起调水位、实时水位以及入库流量作为泄流量的影响因素。以田庄水库下泄流量最小为目标函数,采用粒子群算法对18场洪水进行优化调度,得到水库洪水优化调度数据集。以18场洪水优化调度过程逐时段的累计净雨量、累计入库水量、起调水位、实时水位、入库流量以及泄流量作为RF 回归算法的输入数据,回归树数量变化范围为100~1 000,计算MSEOOB绘制于图3 中。由图3看出,当k=700 时MSEOOB最小,为32.6;当k=100 时MSEOOB最大,为37.3,两者相差4.7,故确定最优回归树的数量k=700。同时,由RF 回归算法计算出的5 个影响因素对泄流量的影响程度分别为累计净雨量21.8%、累计入库水量16.1%、起调水位21.5%、实时水位16.0%、入库流量24.6%,影响程度最大的为实时水位,最小的为累计入库水量,两者相差8.6%,差距较小,故将上述5个因素均作为RF回归算法的输入变量,即mt=5。
图3 MSEOOB随回归树数量的变化曲线Fig.3 The variation curve of MSEOOB with the number of regression trees
由表1 看出,小、中、特大洪水的数量分别为11 场、6 场和1场,故发生不同级别洪水的概率由大到小依次为小洪水、中洪水、特大洪水;洪水级别越大,对水库及其上下游造成的威胁越大,因此,选择洪号为19640912 与20010804 的两场中洪水作为验证样本,其余16 场洪水作为训练样本,将训练样本逐时段的5个主要影响因素和泄流量作为输入,由RF 回归算法生成泄流量回归树,以此求解验证样本的泄流量。依据图2的计算流程,求得19640912 与20010804 两场洪水的泄流量及库水位,与实测入库流量、出库流量及库水位一并绘于图4、5 中。本文比较了两场洪水实测调度和RF 回归算法调度结果的最大泄流量和最高水位两个指标,计算了两者的绝度偏差和相对偏差,见表2。由表2 可知:①两场洪水最大泄流量的实测值均大于RF 回归算法的计算值,19640912 洪水的绝度偏差和相对偏差的绝对值均小于20010804 洪水,由此可见,从保证下游安全而言,RF回归算法的调度结果优于实测调度;②两场洪水最高水位的实测值均大于RF 回归算法的计算值,19640912 洪水的绝度偏差和相对偏差的绝对值均大于20010804 洪水,由此可见,从保证水库安全而言,RF回归算法的调度结果同样优于实测调度。
表2 RF回归算法计算值与实测值对比表Tab.2 Comparison between measured values and calculated values by RF regression algorithm
图4 19640912洪水过程线Fig.4 19640912 flood hydrographs
图5 20010804洪水过程线Fig.5 20010804 flood hydrographs
(1)本文首先对历史洪水过程进行优化调度生成优化调度数据集,然后采用RF 回归算法拟定水库泄流量回归树,最后由泄流量回归树拟定水库调度方案,指导实时防洪调度。该方法既能够提取影响水库泄流量的主要因素,又能够充分挖掘历史调度经验,开辟了水库防洪调度的新方法。
(2)田庄水库的实例可以看出,RF 回归算法制定的水库调度方案的最大泄流量、最高水位均小于实测值,表明RF 回归算法拟定的调度方案较实测调度更能保证水库本身及其下游防洪安全,故是一种合理的调度方案。鉴于田庄水库暴雨洪水特性,本文仅对单峰洪水开展了研究,对于多峰洪水过程有待于进一步的验证。 □