1. 引言
精准的洪水预报可以支撑精细、科学的水库调度,利用洪水资源的同时规避调度中可能形成的防洪风险 [1] [2] 。为了加深对洪水规律认知、提升洪水预报水平,学者们致力于从水文统计、水文模型等多个方面开展系统性研究。洪水分类 [3] [4] [5] [6] 和相似性研究 [7] [8] 是其中一个重要方向,其学科仍处于一个探索与发展过程当中,张静 [9] 将分类理论的概念引入到洪水中,将洪水分类定义为根据洪水的共同点和差异点,采用某一标准和方法,依据分类原则,将洪水进行归类 [10] 。张艳平 [11] 则在其基础上进一步提出,洪水分类的原理需根据洪水本身特征及客观规律、分析洪水成因、时程分配和空间分布等相关性和差异性,以探寻场次洪水发生、发展规律及产生的灾害后果。洪水作为自然现象,不仅具有确定性的变化,同时还伴随有随机性变化规律,洪水特征具有模糊性。根据洪水过程的“三性”特点,洪水分类的方法大体可分为三大类:成因分析法、数理统计法和聚类分析法 [12] 。此外,在洪水分类学中,目前不同学者出于不同的研究目的 [13] [14] [15] [16] [17] ,选择的标准和指标也不尽相同,但原则是要易于识别与比较。
如今,随着大数据、机器学习等智能技术的不断进步 [18] ,在暴雨洪水分类、识别方面有了更多的方法选择,机器学习(Machine Learning, ML)作为人工智能的核心组件 [19] ,主要由机器自主“学习”的各类算法构成,计算机可以通过对大量样本数据进行分析并获取规律,然后利用规律对未知数据进行分类、预测。这与水文预报领域通过分析历史雨洪数据获取规律,并在其基础上通过模型概化或采用统计方式预测来水高度类似。目前机器学习技术已在气象灾害识别预测、天气指标分类、天气预报预测等方面取得了较好的应用并在不断发展进步当中 [20] [21] [22] 。传统的人工交互式预报存在操作效率低,主观性强等问题,机器学习则具备自主学习和分析的特质,并且效率高、稳定性强,大力发展以机器学习为核心的智能预报技术,是辅助预报工作开展、提升预报作业效率的有效途径。目前,水文预报领域所采用的机器学习算法多数都是有监督学习算法,即样本有标注,但是绝大多数可用数据都没有标签,而无监督学习则无需有标注的数据便能够从大量未标记的数据中自动学习。
因此,本文提出采用无监督学习方法,在洪水分类理论的基础上,开展相似洪水分析与智能推荐模型研究。一定程度上将丰富现有预报手段,补充现有预报方法,为洪水预报提供新的思路和方法,从而更好地支撑梯级水库群调度精细化管理,实现综合效益最大化。
2. 研究区概况及数据
嘉陵江是长江上游左岸的一级支流,流经陕西、甘肃、四川、重庆,干流全长1120 km,落差2300 m。全流域面积15.98万km2,自下而上的主要支流有西汉水、白龙江、东河、西河、渠江、涪江等。本文将嘉陵江流域划分为9个子流域作为研究区(图1)。
嘉陵江流域全年降雨主要集中在6~9月,其中7~9月发生暴雨的概率最大。由于流域内地形复杂,洪水特性受流域下垫面和支流洪水加入影响,使得干支流洪水的组成及遭遇情况各异。嘉陵江洪水过程多呈双峰或多峰形,北碚站单峰3~5 d,复峰可达7~12 d,峰现持续时间约4 h左右。
![](//html.hanspub.org/file/4-2411293x8_hanspub.png?20231218113043047)
Figure 1. Schematic diagram of the Jialing River basin and its main stations
图1. 嘉陵江流域及主要控制站点示意图
对于气象数据,本研究主要收集了东亚地区位势高度场、风场1951~2021年逐日再分析资料以及流域内49个雨量站点(图1) 1951~2021年逐日降雨资料;对于水文数据,选择嘉陵江涪江小河坝站、干流武胜站、渠江罗渡溪站、干流北碚站,收集了4站1951~2021年逐日流量资料。
3. 研究方法
本文以嘉陵江干流北碚站为例对相似洪水推荐模型进行分析验证。采用机器学习算法中两个常见的非监督学习算法(主成分分析和K-means法)分别对降雨洪水特征进行降维和聚类,从降雨洪水特征指标中提取关键信息;然后利用洪水模式进行场次洪水的相似类型判别和特征距离度量,并利用暴雨洪水两种模式建立暴雨洪水分类器,构建完整的相似洪水智能推荐模型。具体流程如图2所示。
![](//html.hanspub.org/file/4-2411293x9_hanspub.png?20231218113043047)
Figure 2. Similar flood recommendation technology roadmap
图2. 相似洪水推荐技术路线图
3.1. 暴雨及洪水特征指标相关性分析
本文在研究暴雨和洪水之间的关系时,为充分分析暴雨及洪水过程特征,选取的特征指标较多,例如对于一场洪水的降雨特征指标包括峰现时间前1 d、前2 d、前3 d、前5 d、前7 d及前10 d的各分区累计面雨量,其中流域分区包括子流域和全流域。如果基于全部指标进行分析,不仅难以顾全指标间的独立性,同时可能存在信息冗余等问题,影响分类结果的准确性和可靠性。而直接从众多的指标中直接剔除某些指标又会造成信息的丢失,因此,本研究在区分不同洪水类型之前,先对所有场次洪水的特征指标做一个相关性分析。
本文采用Pearson相关系数来衡量各个指标间的相关程度。Pearson相关系数表达式如下:
(1)
式中:ρ表示相关性系数,取值范围为−1 ≤ ρ ≤ 1。ρ的正负代表两个变量的变化方向是否相同,ρ的绝对值越大代表相关性越强。
3.2. 暴雨及洪水指标主成分分析
本文采用的降维算法是主成分分析方法(Principal Component Analysis, PCA),其主要原理是将原来的多维变量通过一个正交变换重新组合成一组新的、数量减少的、且互不相关的综合变量,并根据实际需要从这些综合变量中选取出几个贡献率最大的作为主要成分来反映原来变量信息。
假设数据集
。其中样本
。计算基于主成分的输入变量数据矩阵
(2)
由该矩阵计算主成分z1,z2,…,zm的数据矩阵。
3.3. 场次洪水聚类相似推荐模型
3.3.1. 基于K-Means算法的洪水聚类分析
K-means聚类算法是一种常见的非监督学习算法,常用于非线性降维和异常监测等,其原理就是随机选取k个对象作为初始的聚类中心,然后计算每个对象与各个种子聚类中心之间的距离,把每个对象分配给距离它最近的聚类中心,如图3所示。
![](//html.hanspub.org/file/4-2411293x14_hanspub.png?20231218113043047)
Figure 3. Diagram of K-means algorithm
图3. K-means算法示意图
计算场次洪水过程特征指标到聚类中心最短的距离,将其归为一类。再重新计算聚类中心进行聚类,直到没有对象被重新分配给不同簇或聚类中心不发生变化,即得到最终的聚类结果。根据输入洪水的特征指标到聚类中心的欧式距离最近,得到相似洪水搜索结果。做分类相似处理。具体步骤如下:
a) 根据上节得到的主成分分析结果。假定由n场洪水,每场洪水q个主成分,这n × q个数据构成一个洪水主成分矩阵,即:
(3)
b) 首先按经验给一个初始聚类类数k,然后将所有样本分成k个初始类,进而将这k个类的重心作为初始的类中心点。
c) 得到初始化聚类结果。计算每个对象到k个聚类中心得距离,把每个对象分配给离它最近的聚类中心所代表的类别中,全部分配完毕即得到初始化聚类结果,聚类中心连同分配给它的对象作为一类。
d) 重新计算聚类中心。得到初始化聚类结果后,重新计算每类的类中心点,得到新的聚类中心。
e) 迭代循环,得到最终聚类结果。重复第三步和第四步,直到满足迭代终止条件。
f) 统计出不同样本的归属类别。
3.3.2. 相似洪水推荐模型
在暴雨和场次洪水相似性度量的时候,选取场次洪水峰现时间前1 d流量与后1 d流量、前3 d流量与后3 d流量、洪峰流量M、洪峰总历时D、涨水时长Dz、涨水点到峰值点的斜率k为特征指标。
对于场次暴雨的每个特征而言,都是一个标量,单特征之间的距离只需要使用曼哈顿距离表示即可。对于完成降雨过程而言,可以用多特征向量
表示,其中i表示降雨过程的编号,Tji表示降雨过程i的第j个特征。那么对于两场洪峰的多个特征之间的距离采用2个向量之间的欧式距离表示。假设降雨过程1对应的多特征向量为
,降雨过程2对应的多特征向量为
,对数据进行z-score标准化后,那么两场洪峰之间的度量DH就可以表示为:
(4)
上述过程介绍了如何对洪水特征指标的选取以及如何进行相似度计算。按照场次洪水提取、雨洪特征提取、降维、聚类的思路,我们提出了一种雨洪聚类推荐模型(图4),这里只需要以降雨数据、水文数据为基础输入即可,可以针对雨洪模式实现聚类划分,和模式与模式之间的对应关系。这个分类器在划分雨型和洪型对应关系时(从时间上雨和洪有对应关系),可以从两者在成因上形成递进的经验概率,从而在推荐相似雨洪时,按照最有可能发生的一种进行推荐。
![](//html.hanspub.org/file/4-2411293x20_hanspub.png?20231218113043047)
Figure 4. Architecture diagram of the rain-flood recommendation model
图4. 雨洪推荐模型的架构图
4. 结果分析
4.1. 相关性分析结果
北碚站洪水特征指标选取起涨流量、洪峰流量、退水点流量、洪水总历时、涨水历时、退水历时、最大2 d、3 d、5 d、7 d洪量、场次洪水平均流量、涨水面平均流量、退水面平均流量、流量总涨幅、流量总退幅、流量起涨率、流量最大涨率、流量最大落率、高脉冲历时、变差系数、峰现时间比例、洪峰时间偏度、洪量集中度、涪江桥以下、风滩以上等分区7 d累计降雨等指标。对上述33个特征指标进行相关性分析,结果如图5所示。
![](//html.hanspub.org/file/4-2411293x21_hanspub.png?20231218113043047)
Figure 5. Correlation analysis results of the characteristics of heavy rain and flood indicators at the Beibei station
图5. 北碚站暴雨洪水特征指标相关性分析结果
由图5可知,在北碚站,最大2 d、3 d洪量与6个特征指标的相关系数绝对值大于0.9,洪峰流量与5个特征指标的相关系数绝对值大于0.9,最大5 d洪量、流量总涨幅、流量总退幅与4个特征指标的相关系数绝对值大于0.9。从上述结果可以发现,洪水特征指标内部之间的相关性较高,而降雨特征指标(分区面雨量)之间的相关性较低。对于洪水特征指标与降雨特征指标之间的相关性结果,在北碚站,亭子口以下7 d累计降雨与洪水特征指标间的相关性最高。可以看出站点的洪水的特征指标与所在流域的降雨特征相关性更强,这一结果也符合流域产汇流规律。
综上所述,在描述洪水特征的33项指标中,洪量特征指标之间的相关度高,说明描述洪量的特征存在信息冗余,因此在做后续分析前应删除部分特征变量。考虑到挑选出的部分场次洪水不足7 d甚至不足5 d,因此选择删除最大5 d洪量、最大7 d洪量两个特征指标,剩下的31个特征变量做后续研究分析。
4.2. 主成分分析结果
北碚站降雨特性指标主要包括峰现时间前1、3、5、7及10 d的流域分区累计面雨量,流域分区包括9个子流域以及全流域,因此每一场洪水有60个降雨特征指标。对北碚站的降雨特征指标和洪水特征指标分别进行主成分分析,如图6所示。
由图6可知,北碚站降雨特征指标前6个成分的方差累计贡献率接近0.85,说明6个成分可以保留原始的降雨特征信息,所以,本研究在北碚站均提取6个主成分因子;北碚站洪水特征指标前3个成分的方差累计贡献率都接近0.9,说明3个成分可以保留原来的洪水特征指标信息,所以,本研究在北碚站均对洪水特征指标提取3个主成分因子。
利用主成分分析方法,北碚站点部分场次洪水降雨指标降维后的6个主成分因子的结果,具体见表1。
![](//html.hanspub.org/file/4-2411293x22_hanspub.png?20231218113043047)
Figure 6. Cumulative contribution rate results of indicators for evolutionary characteristics (left) and flood characteristics (right) of Beibei Station
图6. 北碚站降雨特征(左)、洪水特征(右)指标累计贡献率结果
![](Images/Table_Tmp.jpg)
Table 1. Results of three principal component factors after dimensionality reduction of rainfall index
表1. 降雨指标降维后三个主成分因子结果
北碚站在对洪水特征指标选取3个主成分时,首先计算了主成分因子载荷矩阵,根据这个矩阵可以计算出各个站点的每场洪水的综合指标结果,见表2。
![](Images/Table_Tmp.jpg)
Table 2. Principal component factor loading matrix results of Beibei station
表2. 北碚站主成分因子载荷矩阵结果
由表2可知,对于第1个主成分,贡献最大的是洪峰流量、最大2 d、3 d洪量;对于第2个主成分,贡献最大的是退水历时、退水点流量、洪水总历时;对于第3个主成分,贡献最大的是涨水历时、流量起涨率、洪峰时间偏度。
根据北碚站的主成分因子载荷矩阵,计算出其所有场次洪水3个主成分因子的结果,见表3。
![](Images/Table_Tmp.jpg)
Table 3. Results of three principal component factors after dimensionality reduction for flood events
表3. 场次洪水降维后三个主成分因子结果
4.3. 聚类分析结果
根据相似洪水聚类分析方法,对北碚站挑选出的场次洪水做分类相似处理。分类结果如图7所示。
![](//html.hanspub.org/file/4-2411293x24_hanspub.png?20231218113043047)
Figure 7. Flood clustering results at Beibei station
图7. 北碚站场次洪水聚类结果
由图7可知,北碚站的所有洪水过程共分为15个类型,该站点洪峰排名靠前的洪水过程被划分到了类型4和类型9,而这两种类型的洪水主要差别体现在涨水历时这个特征上,类型9的洪水涨水历时更短。
4.4. 相似洪水推荐结果
以嘉陵江流域编号20110921的大洪水为例,需要匹配的时段为2011年9月15日至9月19日的暴雨天气过程,在历史上连续5 d的降雨过程匹配搜索中,与本次过程最为匹配的暴雨天气过程为1973年9月3日至9月8日(图8)。完成降雨匹配后,通过提取不同分区累计雨量的特征,并结合起涨点流量、洪峰流量、洪水总量等信息,通过对比多个洪水样本的相似度(欧式空间距离的度量值),模型推荐的19730909这场洪水相似度为1.197,优于其他选项(见表4),因此无论从降雨还是洪水过程,场次编号19730909洪水与20110921洪水所表现的特征最为接近,见图9。
![](Images/Table_Tmp.jpg)
Table 4. Comparison of similar flood characteristics between 20110921 and 19730909
表4. 20110921与19730909相似洪水特征对比
![](//html.hanspub.org/file/4-2411293x35_hanspub.png?20231218113043047)
Figure 9. Comparison of similar flood hydrographs between 20110921 and 19730910
图9. 20110921与19730910相似洪水过程线对比
针对嘉陵江流域编号20210907的大洪水同样进行了个例检验,需要匹配的时段为2021年9月2日至9月6日的暴雨天气过程,在历史上连续5 d的降雨过程匹配搜索中,与本次过程最为匹配的暴雨天气过程为2007年7月2日至6日(图10)。完成降雨匹配后,通过提取不同分区累计雨量的特征,并结合起涨点流量、洪峰流量、洪水总量等信息,通过对比多个洪水样本的相似度(欧式空间距离的度量值),模型推荐的为20070708这场洪水,见图11。
![](//html.hanspub.org/file/4-2411293x36_hanspub.png?20231218113043047)
Figure 10. Similar rainfall processes in September 2021 and July 2007
图10. 2021年9月与2007年7月相似性降雨过程
![](//html.hanspub.org/file/4-2411293x37_hanspub.png?20231218113043047)
Figure 11. Comparison of similar flood hydrographs between 20210907 and 20070708
图11. 20210907与20070708相似洪水过程线对比
5. 结论与展望
本文收集整理嘉陵江流域水文气象数据资料,依据雨洪特征刻画指标量化分析方法,结合无监督学习方法,实现了相似雨洪特征的识别、提取,构建相似雨洪特征工程库。采用降维和聚类算法实现了场次洪水的类型分类,北碚站的历史洪水过程可以通过前面保留的主成分因子聚类分成15个类型。完成分类后通过欧式空间距离计算的属性值来度量实体与同类型洪水的相似程度,实现了洪水相似推荐。最后,以北碚站的2场洪水作为相似洪水推荐演示,获得较好的相似洪水推荐结果,结果表明,该方法可以为实际业务中洪水预报应用提供参考。
随着水利信息技术的快速发展,业务需求、应用方向的不断深化,面向各类应用场景的人工智能技术研究也需不断提升和深化,特别是要进一步深挖已研究成果的落地应用,通过应用推动以知识和数据驱动的模型在预报调度领域的迭代升级,进而使现有的业务体系达到新的高度。基于水利智能化的发展将对防洪、预报提供新的方法和思路。
基金项目
国家重点研发计划项目(2022YFC3002705, 2022YFC3002701);中国长江三峡集团有限公司科研项目资助(0799251, 0704189);水利部重大科技项目(SKS-2022128)。