小流域暴雨山洪水文模型与水动力学方法计算比较分析

刘昌军，文 磊，周 剑，赵悬涛，郭 良，魏永强

（1.中国水利水电科学研究院 减灾中心，北京 100038；2.河海大学 水文水资源学院，江苏 南京 210098；3.中国科学院 西北生态环境资源研究院，甘肃 兰州 730000；4.华北水利水电大学 水利学院，河南 郑州 450045；5.湖南省水利水电科学研究院，湖南 长沙 410007）

1 研究背景

我国山洪灾害点多面广、突发性强、破坏力大，且多发生在偏远山区、交通不便、通讯不畅。山洪灾害是严重威胁群众生命财产安全的主要灾种之一，日渐成为防洪减灾的短板，防御工作一直是我国防汛和防灾减灾工作中的难点和薄弱环节［1］。近年来无资料小流域暴雨山洪模拟研究得到快速发展［2］，目前常用的是采用水文学方法和水动力学方法［3］。基于分布式水文模型的小流域暴雨洪水模拟得到快速发展和广泛应用［4-5］。黄燕等［6］开展了水文模型在山洪模拟中的比较应用研究，研究结果发现现有水文模型在模拟短历时山洪效果一般。

随着计算机技术的快速发展，采用水动力学方法进行暴雨山洪过程的模拟是今后研究方向和热点［7-8］。平面二维数值模型已成为洪水演进模拟的主流方法，但间断水流、干湿边界处理和计算效率仍存在较多难点，是急需解决的关键问题［9］。潘佳佳等［10］认为基于运动波和扩散波等简化水动力学模型因忽略了流动的力学因子，不能较好的模拟山洪的形成与演化机理。张培文等［11］采用有限元法研究了降雨条件下坡面径流和入渗问题的耦合模拟；王娇进行了坡面径流—入渗耦合的解析解及无网格数值模拟研究［12］；基于Godunov格式的以求解黎曼近似解的浅水数值模拟技术已逐渐成熟，该方法能模拟光滑解，又能自动适应复杂流态变化，并得到广泛应用［13-14］。张大伟等［15］提出了新的基于Godunov格式的地表径流二维水动力模型，使干湿边界问题得到很好解决。张奕虹等［16］提出了新的水深重构方法，解决了复杂地形上暴雨山洪数值模拟常面临的计算水深小于零、非物理大流速等问题。阚光远等［17-18］在GPU并行算法方面作了大量研究，发现使用GPU并行算法可大幅提高计算效率。GPU、CPU等并行计算技术已经在水文和水动力学技术中得到初步应用，具有广阔的应用前景［19］。

针对暴雨山洪模拟计算问题，以湖南浏阳市宝盖寺小流域2012年降雨径流过程为例，本文采用笔者提出的时空变源混合产流模型和模块化小流域暴雨山洪分布式水文模拟模型以及基于Godunov格式的二维坡面径流数值模拟模型对小流域暴雨山洪进行了对比研究，从建模过程、计算效率和应用等方面，对比分析了两种方法计算精度和计算效率及计算结果合理性，系统分析了坡面产汇流机制和微地形对降雨径流过程影响，为无资料小流域暴雨洪水计算和预报预警提供技术借鉴。

2 水文模型及软件简介

在山洪灾害调查评价成果和实时水文气象大数据基础上，笔者提出了以实时卫星、雷达和台站等多源降水融合数据和预报降雨数据为驱动，以调查评价小流域拓扑关系及基础属性为基础，以7类水文单元（流域、河段、节点、分水、水源、洼地、水库）水力关系为核心，以历史暴雨洪水和率定参数为先验知识，以机器学习和并行计算为手段，以降雨蒸发、产汇流、河道演进、水库调蓄等水文过程为主线，基于人工智能和大数据技术，构建以小流域为计算单元，集成模型库、人工智能算法库、参数库和先验知识库为一体的模块化分布式水文模型（FFMS）［20］，模型计算流程和结构见图1。

图1 分布式水文模型流程和结构

2.1 时空变源混合产流模型围绕山丘区中小流域地形地貌多样、产流机制时空变化复杂等问题，划分了山丘区山坡地貌水文响应单元，采用基于一维入渗理论的包气带土壤非线性下渗计算方法，计算湿润锋下移过程及入渗量，提出了超渗、蓄满时段转换准则，建立了在平面、垂向、时段上时空变源混合产流模型［21］，以精确模拟山丘区中小流域产流时空分布过程。

（1）平面混合产流模型。提出了山坡地貌水文响应单元的划分标准，将山洪地面响应单元划分为快速响应单元（根据差异性又分为快、中、慢三种类型）、滞后响应单元（根据差异性又分为快、中、慢三种类型），贡献较小单元和其他单元等。研究不同山坡地貌响应单元下垫面参数的差异，建立各类山坡地貌水文响应单元与产流机制的对应关系，实现超渗/蓄满产流机制的平面上混合产流模型搭建。

（2）垂向混合产流模型结构。时空变源混合产流模型的垂向结构是综合考虑了截留、填洼、超渗产流、蓄满产流、优先流、壤中流、基流和渗漏出流等过程。根据上述产流组分形成过程，基于非饱和下渗理论和水量平衡，采用不同的概念水库来模拟土壤包气带和饱和带，从上至下依次是表层土壤的毛细水库、下层土壤的变动面积饱和产流水库和地下水库，各水库产流的相关变量及参数见表1，垂向模型结构如图2。

（3）非饱和土壤下渗。非饱和土壤下渗采用Talbot和Ogden［20］提出的一种新的一维入渗和水分再分配的计算方法。将土壤含水量范围离散成具有恒定含水量的区间，推导一维入渗方程以快速获取湿润锋垂向位移，并考虑不同含水量单元之间的互吸力，实现再分配湿润锋横向位移，实现高效稳定的包气带土壤下渗计算方法，获得包气带土壤的下渗过程；本文采用此方法进行小流域包气带土壤非线性下渗计算，该方法采用数值化方案，无需估算非线性梯度，使得最终的土壤下渗过程具有数值稳定性。

表1 时空混合产流模型主要变量及参数

图2 垂向混合产流模型结构

2.2 FFMS水文模型软件模块化小流域暴雨洪水分析软件FFMS由中国水利水电科学研究院刘昌军博士负责研发。该软件基于C++、FOTRAN和JAVA等语言开发完成，具有模块化、参数化、智能化、可视化和自动化特点，软件主界面见图3。软件集成了包括山洪灾害调查评价成果数据导入导出、小流域自动划分及参数提取算法、模块化建模环境、地貌水文响应单元自动划分、拖拽式手动和全自动建模、计算结果显示等分布式水文模型的前后处理模块，实现了自由组合建模和全自动模块化建模功能；集成了多种国内外分布式水文模型的气象、蒸发、产汇流、坡面侵蚀、河道侵蚀和洪水演进等50余个算法模块，具有参数自动率定和手动率定算法模块等多个参数全局自动优化算法。同时该软件还集成了各种模型工具库、机器学习参数库和专家知识库。详细功能模块见图4。

图3 FFMS软件主界面

图4 FFMS软件架构与功能模块

3 水动力学模型简介

3.1 控制方程本模型应用动力波方法进行小流域坡面山洪过程的模拟计算。采用Godunov格式的有限体积法将计算区域进行空间离散，通过HLLC近似黎曼求解器计算单元界面上的质量通量和动量通量［22］。本文应用具有守恒格式的平面二维浅水方程见式（1）—式（2），鉴于模型和算法较为成熟，本文不再详细介绍，模型详细求解方法见文献［22］。

式中：q为变量矢量，包括水深h，两个方向的单宽流量qx和qy；g为重力加速度；u和v分别为x、y方向的流速；F和G分别为x、y方向的通量矢量；S为源项矢量；zb为河床底面高程；为谢才系数，n为曼宁系数。

3.2 水动力学模型软件简介刘昌军等采用OPENGL和C++语言开发了三维可视化水动力学模拟软件FHMS，该软件具有点云数据处理、DEM处理，网格剖分、参数设置、模型计算、参数率定和计算成果可视化展示等功能。本软件可视化功能强，操作简单，具有广泛的推广应用前景，限于篇幅所限，不再详细介绍，软件主要功能界面见图5。

4 两种方法对比分析

4.1 流域基本情况以湖南宝盖寺小流域为例，开展小流域暴雨山洪水文水动力学计算对比研究。小流域内有宝盖洞清水水文站，上游建有寒婆坳雨量站，小流域集水面积22.1 km2，河长9.1 km。小流域下垫面、土壤类型见图6。选取2012年5月9日—5月10日的实测一场降雨径流过程分别进行水文和水动力学模型对比计算分析。

图5 三维水动学计算软件主要界面

图6 小流域土地利用和土壤质地类型

4.2 水文建模与参数基于无人机获取的1 m分辨率高精度DEM数据，利用FFMS软件全自动建立了宝盖寺小流域分布式水文计算模型，共划分15个模型计算单元，见图7，模型参数采用2012年6、7月份的两场实测暴雨洪水资料进行率定，得到此次模拟的模型参数取值，见表2。

表2 水文模型计算参数

图7 小流域高精度DEM、计算单元划分及坡面局部三维网格图

4.3 水动力学建模与参数基于高分辨率DEM数据，应用FHMS软件，建立了宝盖寺小流域水动力学计算模型，计算网格模型采用非结构化三角网，共划分网格单元1 764 215个，见图7。该流域面积较小，地貌均质性较好，以林地为主。糙率值是影响模拟结果的最重要参数，以模型确定性系数为衡量标准，采用人工试错法对糙率值进行率定，最终确定计算区域内的坡面糙率值为0.09，沟道糙率值为0.03。

图8 不同时间河道流量分布

图9 不同时间的坡面淹没水深

图10 水文模型与水动力学方法计算值和水文站实测值

图11 3个断面水文模型与水动力学模拟洪峰流量过程

4.4 计算结果分析采用水文学和水动力学方法分别对该场次洪水进行了计算，水文学方法计算得到不同时间各河段洪峰流量平面分布见图8，水动力学计算得到的河道淹没水深见图9。流域出口两种方法的洪峰流量计算值和实测值见图10。选取流域内三个典型断面对比分析两种方法的洪峰流量见图11，流域出口及各断面计算的洪水峰现时间、洪峰流量结果见表3。

表3 洪水要素对比

（1）从图8和图9计算结果可以看出，水文模型和水动力学模型计算得到各子流域河道洪峰过程和洪峰流量基本一致。在降雨2 h后，9时左右在水文站下游河道开始形成径流；10时小流域内各支流河道均局部形成径流，坡面局部地方也开始出现径流；11时至12时河道内洪水逐渐形成，并向下游演进。二维水动力学计算得到的坡面径流产生、汇流、演进过程符合实际情况，与水文学方法计算结果一致。

（2）从图10和表3可以看出，流域出口处两种方法的模拟结果和水文站实测结果基本一致。此次洪水过程最大实测洪峰流量为27.7 m3/s，水文模型计算最大洪峰流量为31.6 m3/s，误差约14%；水动力学模拟的最大洪峰流量为34.1 m3/s，误差为23%。水文模型计算得到的峰现时间与实测值基本一致，水动力模型计算得到的峰现时间提前约1 h。

（3）从图11中可以看出，三个典型断面的水文学和水动力学计算结果基本一致，水动力学在3个断面处均模拟出现了3次洪峰过程，而水文学模拟结果出现一次洪峰，过程曲线较为光滑，说明水动力学对洪峰过程的刻画更加精细。

（4）本次采用的分布式水文模型在产流计算方面采用时空变源混合产流模型，此模型根据不同山坡地貌响应单元的下垫面类型的差异，确定各类山坡地貌水文响应单元所对应的产流机制；建立并求解非饱和下渗曲线方程以快速获取湿润锋垂向位移，考虑不同含水量单元之间的互吸力，实现湿润锋横向位移再分配，得到包气带土壤水下渗过程；确定短历时暴雨洪水的土壤最大含水量参数，完成垂向和时段上超渗产流、浅层土壤蓄满产流、深层土壤蓄满产流机制时空转换；有效解决了山丘区中小流域洪水预报面临的水文数据精度低和产流机制时空变化问题，为无资料小流域的产流计算提供了新思路。

总的来看，水文学和水动力学方法对小流域暴雨山洪模拟结果都很好，说明自主研发的模块化小流域洪水分析软件和水动力计算方法均可用于暴雨山洪模拟。时空变源分布式水文模型模拟的结果与实际观测结果更加接近；地表径流二维水动力模型能够处理实际流域的复杂地形，且可以处理具有间断水流情况，具有较高的计算精度。

5 结论

本文基于激光高精度激光点云数据，应用自主研发的时空变源混合产流模型、模块化小流域分布式水文模型FFMS和二维水动力模型FHMS研究了宝盖寺小流域坡面暴雨山洪过程。采用小流域出口水文站实测暴雨洪水资料验证了水文学和水动力学方法计算结果，并对两种方法计算暴雨山洪的效果、效率和精度进行对比分析，两种方法均可用于小流域暴雨山洪计算分析，各有优缺点，在山洪预报预警方面水文学方法更加简单实用，在小流域暴雨洪水机理、演进过程和精细化分析方面水动力学更加精细和准确，本文提出的时空变源混合产流模型为无资料小流域暴雨洪水产流计算提供了新的研究思路，研发的水文学和水动力学计算软件操作简单、可视化功能强大、计算效率高，具有较大的推广应用价值。