我国西南山区地质构造活跃、地形变化剧烈、雨热空间分异明显，人口和经济高密度区域常同山洪灾害风险区重叠．在全球气候变暖使得气候变化更加剧烈、极端降雨事件的发生更加频繁[1-2]的背景下，极端降雨伴随着的山洪灾害越来越成为防灾减灾所关注的热点问题[3]．据不同自然灾害损失数据统计，山洪灾害造成的损失在全球105个国家的自然灾害损失中排前两位，年均经济损失超过460亿美元，占全球自然灾害类型的50%[4]．自20世纪50年代以来，山洪灾害每年都会给我国造成数以百计的人员伤亡和巨大的财产损失，是我国自然灾害中的主要灾种之一[5]．1950—1990年，因山洪灾害直接或间接死亡人数达15.2万；90年代之后，全国防洪体系逐渐建立，因山洪灾害死亡人数呈现逐年下降趋势．然而近年来，伴随着极端气候的频发，山洪灾害导致的人员伤亡和经济损失依然居高不下[2,6]．山区小流域环境与平原、丘陵等地貌环境具有显著区别，受山地下垫面复杂环境条件强烈影响，当前暴雨天气条件下山洪灾害的演进和预警预报能力欠缺，灾害预报主要以基于临界降雨的宏观等级(高、中、低)预报为主，无法给出流域的具体灾情(发生时间、危害范围和风险程度)．通常，由于山区小流域地形坡度大，降雨-洪水的演变快速，且山区基础数据较为缺乏，因此山区小流域洪水比平原地区或大江大河流域更加难以预测[7-8]．伴随着计算能力的提升和小流域监测数据的增多，洪水动力学数值模拟理论与方法都得到了快速的发展[9-11]．其中，二维浅水波方程在洪水模拟中得到了很好的发展[12-15]并出现了多款商业软件，比如Massflow[16]，MIKE URBAN[17]，Infoworks CS[18]，还有开源软件HEC-RAS[19]．然而，山区小流域洪水的动力学研究依然面临着很多挑战．一方面，传统山洪计算模型主要聚焦于径流的演进过程而不考虑植被截留和沟道基流对径流的影响．但是在植被茂密区，植被对降雨的截留可超过20 mm[20-21]，所占比例为一次特大暴雨过程累积雨量的15%以上，其在雨洪过程中的作用不可忽略[22]．在雨季，山区小流域沟道内的基流较大，基流能显著改善汇流-径流过程的稳定性，对径流的计算时须考虑其影响．另一方面，一个州或者县域范围内包含成百上千条小流域山洪沟，同时山洪沟道狭窄，需要精细化的网格才能有效描述，因此对计算机硬件和计算效率要求极高．为此，近年来基于GPU(图像处理器)并行计算的数值计算方法已经取得了非常显著的效果并成为一种趋势．例如，Sanders等[23]采用高性能水动力模型模拟了大流域洪水动力学过程，Liang等[24-26]实现了城市洪水事件的实时模拟，Kalyanapu等[27]基于GPU实现了15倍真实事件的仿真速度．随着超算的高速发展，高效利用超算平台上成百上千的GPU资源是解决大区域众多小流域山洪实时快速计算的保障；因此，本研究基于北京超算中心GPU计算平台，结合当前的研究不足，实现从降雨到沟道径流的全物理过程山洪动力学模型的高效求解，相比传统的降雨入渗径流过程增加并耦合了植被截留和沟道基流，提升了小流域山洪模型的计算精度．相关研究能够为洪水研究过程的机理认识提供参考，同时也能为未来基于流域下垫面动力学全过程的山洪预警预报提供技术手段．1 模型构建山区小流域山洪灾害是暴雨驱动的快速演变过程，形成和演进涉及小流域下垫面内复杂多物理过程．开展小流域山洪演进过程的模拟须综合考虑从降雨开始到沟道基流、植被截留、坡面入渗，再到径流的产生和运动演进整个物理过程(如图1所示)．在扩散波模型基础上引入基流、Aston植被截留模型和Green-Ampt坡面入渗模型，模拟从降雨到山洪演进整个物理过程．10.13245/j.hust.239410.F001图1考虑下垫面过程演进的小流域山洪模型1.1　物理过程和模型a．坡面汇流-洪水径流过程模型坡面汇流-洪水径流过程采用基于纳维-斯托克斯方程简化后的扩散波模型，基底阻力采用曼宁模型计算，其控制方程为∂h∂t+∂hu∂x+∂hv∂y=R-V-I；(1)∂zs∂x=-n2uu2+v2h4/3；(2)∂zs∂y=-n2vu2+v2h4/3，(3)式中：t为时间；x和y为笛卡尔坐标系坐标；h，u和v分别为流体高度和在x，y方向的速度分量；R，V和I为降雨、植被截留和坡面入渗引入的物源项；n为曼宁粗糙系数；Zs=Zb+n为流体自由面的高程，Zb为地面高程．联立方程(2)和(3)求解可以把流体速度表示成当前坡度、流体深度和曼宁系数的函数，然后通过方程(1)描述的质量守恒方程求解下一时刻的流体深度．扩散波模型通过忽略动量守恒方程中的惯性项，解耦质量和动量守恒方程，极大地简化了模型求解，成为山洪演进的主要模拟方法之一．b．植被截留基于叶面积指数的植被截留模型[28]为V=Smax(1-e-kPc/Smax)，(4)式中：Smax为植被最大拦截量，根据不同的植被类型和叶面积指数确定[29]；Pc为当前累计降雨量；k为与植被郁闭度Co和叶面积指数(L)相关的参数，k=1-e-CoL．(5)叶面积指数是指单位水平地面面积上单面叶面积的总和，在植被截流模型中决定着最大的截流总量．可采用地表实测叶面积指数与高分辨率遥感影像提取的植被指数进行反演确定，L=13.746 9N3.646 4，(6)式中N为植被指数．c．地表入渗降雨条件下水的入渗过程决定了土体中水的含量和地表径流，是影响汇流和山洪径流的关键物理过程．对于入渗模型，Green-Ampt模型有明确的物理基础，改进后的模型在不均匀降雨入渗计算中有很好的应用，表达式为I=dfdt=ks(ψf+h)θs-θif+1，(7)式中：f为累计入渗深度；ks为饱和导水率；ψf为土壤湿润锋压力水头；θs为土壤湿润锋后方的饱和含水量；θi为土壤湿润锋前方的土壤初始含水量．d．沟道基流根据监测数据统计分析可知，在无强降雨情况下，日均基流量基本恒定，采用沟道平均基流系数qbase定量描述，表达式为qbase=Q/l，(8)式中：Q为雨季无降雨时沟口稳定流量；l为流域内沟道长度．1.2　数值求解与高效并行实现扩散波方程通过一阶迎风有限差分格式求解，动量方程中流体深度与速度是解耦的．可先通过连续方程得到下一时刻的流体深度，再求解动量方程获得流速．采用迎风格式求解方程时，凸地形节点本身的通量为零，无法求出特征波的方向，该点的流体高度只会随着降雨持续而增加，并不会向周围流动．因此，须对凸地形点进行流深修正，以x方向为例，流出到左结点(i-1，j)的修正高度为hcorr(i-1,j)=Δtdx[-qx(i-1,j)]，(9)向右节点(i+1，j)为hcorr(i+1,j)=Δtdx[qx(i+1,j)]，(10)点(i，j)自身为hcorr(i,j)=Δtdx[-q(-1,j)i]+Δtdx[qx(i+1,j)]，(11)时间步长须满足Courant-Freidrichs-Levy条件Cr=SΔtΔx，(12)式中：Cr为无量纲库朗数，为保证求解格式的稳定性，小于1；∆t为时间步长；∆x为网格大小；S为特征速度，对于忽略惯性项的浅水波方程，可表示为S=gh，(13)g为重力加速度，最大时间步长为Δtmax=αΔx/ght，(14)其中，α为系数，为保证在大多数流域内保持稳定，取值0.2~0.7；ht为t时刻流深．GPU是一块PCIe板卡，由于设计功能具有差异性，因此GPU显存与CPU存储空间互相独立并通过PCIe插槽与CPU连接．对于CPU+GPU异构并行计算平台，目前主要有NVIDIA公司搭建的CUDA(computer unified device architecture)和AMD公司搭建的ROCM(radeon open computing)两个集成开发生态．基于北京超算中心GPU显卡，采用CUDA编写数值计算程序，程序流程图如图2所示．任务从读取输入文件开始，在读入数字高程模型(DEM)文件的基础上，对网格和边界条件进行初始化．植被截留和土壤入渗参数由土地利用类型确定．参数初始化后，CPU自动分配主机和设备内存，然后调用内核函数．在计算过程中和计算完成后，由输出参数控制输出模型将流量和流速、水深、影响范围等写入输出文件中．10.13245/j.hust.239410.F002图2程序流程图在程序并行过程中主要涉及以下几点关键技术．a．研究区任务自动划分：单个设备显存大小固定且有限，当计算区域数据存储开销过多，会造成数据溢出而无法计算．同时，计算负载过大也会降低设备计算速度．因此，面对任务处理数据量大的情形，须对计算区域合理划分，之后再由不同的设备处理不同细分区域，这一过程须综合设备能力和数据大小，在程序算法设计时予以考虑．b．主机、设备异步多流并行：一个流仅是设备上一系列有顺序的操作，不同流的操作可以交错执行，某些情况下还可以重叠执行，合理利用此性质可隐藏主机与设备间的数据传输消耗．一般情况下，GPU上的所有操作默认在零号流上．结合主机端CPU和设备端GPU的计算特性，程序计算框架考虑为异构多流并行．2 模型验证2.1　一维构造地形实验一维构造地形实验是一个典型的计算分析案例[30-31]．通过比较当前扩散波方程与浅水波方程[32]进行验证．本案例中降雨范围设置在上游(0 m≤x≤20 m)，降雨强度为100 mm/h，持续整个计算过程(50 h)．地形上下游均设为开边界．曼宁摩擦系数设置为0.015 s/m1/3，x和y方向的网格长度均为1 m．一维地形基底表达式为z=42-0.5sinπ10x-0.005x．(15)本案例中两个模型的降雨漫流演进过程对比如图3和4所示．扩散波方程能够提供一个近似于复杂浅水方程模型的演进过程，惯性项的忽略简化了计算复杂性提高了求解速度．10.13245/j.hust.239410.F003图3一维构造地形流深随时间变化（扩散波方程）10.13245/j.hust.239410.F004图4一维构造地形流深随时间变化（浅水方程）2.2　V形槽数值实验V形槽数值实验[33-35]流量过程曲线可以用于验证程序数值求解精度和稳定性．具体模型设计如下：坡面地形坡度为0.05，河道坡度为0.02，河道宽20 m，模型设计示意如图5．坡面和沟道曼宁系数分别取0.015和0.15 s/m3．网格单元为1 m的均匀矩形网格．全流域降雨强度为10.8 mm/h，持续时间为1.5 h，总计算时间为2 h，坡脚和沟道出水口流量(Q)过程曲线如图6所示．10.13245/j.hust.239410.F005图5V型槽设置示意图10.13245/j.hust.239410.F006图6V型槽数值计算结果和解析解对比由图6可以看出：扩散波模型数值计算结果稳定，能够很好地反映坡脚和沟道出水口流量变化，且流量过程曲线在降雨持续阶段和降雨停止后均能够很好地匹配解析解结果，仅出水口位置在流量稳定前期同解析解相比略有误差，且最大相对误差小于5%，充分验证了当前计算模型和算法的稳定性和精度．3 都江堰龙溪河流域山洪过程分析3.1　都江堰龙溪河流域位置为进一步验证当前模型对于长历时降雨山洪事件的模拟适用性和计算性能，以都江堰龙溪河2018年6月26日发生的山洪为例进行数值模拟．龙溪河流域位于四川成都市都江堰市区西北部，属于长江水系的岷江一级支流，地理位置见图7，流域总面积96 km2．10.13245/j.hust.239410.F007图7龙溪河流域地理位置3.2　基础数据及结果分析根据调查，龙溪河土地利用类型主要包括林地、草地、裸地、人造地表、水体，详细地土地利用类型分布见图8．10.13245/j.hust.239410.F008图8龙溪河流域土地利用类型分布结合流域内多个雨量监测点R1~R9(见图7)降雨数据记录，对龙溪河流域2018年6月26日山洪事件进行模拟．龙溪河上游建设有小型水库，为方便计算，库口边界处理时将其设为墙边界．裸地、水体和人造地表不考虑入渗，曼宁粗糙度系数根据土地利用类型确定[36]，在事件发生初期，雨被冠层截取(设L均值为31)．本次模拟不同土地利用类型使用到的水文参数见表1．曼宁系数取值参考文献[37-38]，其他参数根据经验和校准取值．10.13245/j.hust.239410.T001表1龙溪河流域模拟水文参数参数林地草地裸地人造地表水体饱和导水系数/(m∙h-1)2.21.5土壤饱和含水率0.450.40初始土壤含水率0.200.15曼宁系数/(s∙m-1/3)0.1800.1000.0900.0800.015鉴于本次研究模型中新增了基流和植被截留两个因素，须对基流和植被截留对径流的影响进行定量评估．在评估分析中分别对仅考虑基流、仅考虑植被截留及全要素模拟3种工况的降雨径流曲线进行对比．本次研究获得图7中R7点支沟和R9点主沟附近流量监测点数据．据图9主沟R9点流量过程分析(H为降雨量)，对比全要素(同时考虑植被截留和基流)和监测径流曲线，耦合植被截留和沟道基流的山洪动力学模型通过径流过程、洪峰流量及峰现时间3个指标均能较好地反映真实山洪演进过程；不考虑植被截留影响的径流偏大且洪峰出现时间更靠前；不考虑基流影响的模型径流偏小且峰现时间有延迟．据图10支沟R7点流量过程分析可知，第一次降雨过程(6月25日15∶00开始)，仅考虑基流因素时沟内模拟流量出现明显上涨，但监测流量涨落不高，这可能同入渗参数的取值有关，入渗量偏小，但两次明显降雨过程中间时段监测流量同模拟流量匹配性好．第二次降雨过程开始后模拟流量上升过程快于监测流量，造成后期径流消散过快，流量峰值偏低．图11为龙溪河6∙26山洪事件中不同时刻的模拟的流深分布．对本次山洪过程模拟的不同时刻流深分析，6月25日15点(图11(a))降雨开始之前，根据监测数据显示流域内仅存在基流，第一次降雨过程后，6月25日20点(图11(b))流域内径流已经明显增加，6月26日凌晨2点(图11(c))正处于第二次强降雨过程，流域内支沟也形成了较大径流，6月26日5点(图11(d))第二次强降雨过程临近结束径流持续增加，6月26日7点(图11(e))降雨结束后，支沟内径流逐渐汇聚到主沟而减小，主沟径流增大．图11(f)为本次事件模拟记录的流域内历史最大流深分布．10.13245/j.hust.239410.F009图9龙溪河6.26山洪事件主沟R9点流量过程曲线10.13245/j.hust.239410.F010图10龙溪河6.26山洪事件支沟R7点流量过程曲线图11龙溪河6⋅26山洪事件中不同时刻的模拟流深分布10.13245/j.hust.239410.F11a110.13245/j.hust.239410.F11a2综上所述，本次模拟的流量过程曲线和流深分布总体上能够反映真实山洪形成和演进过程．植被截留和沟道基流对支沟的径流过程的影响相对主沟较小，但支沟径流模拟值同监测值的匹配度相比总体优于主沟．本模型引入植被截留和基流提高了模拟流量过程曲线、峰值流量和峰现时间与监测值的符合度，可以明确该模型相对传统模型能够在无资料长历时复杂地形条件下提供更为稳定和准确的结果．3.3　计算加速比分析为分析基于异构的并行技术对水动力学模型计算带来的加速效果，将本文案例基于异构并行计算的时间同CPU计算的时间进行对比，加速效果分析结果见表2．本次研究采用设备为：GPU NVIDIA TESLA-V100(GPU包含5 120个CUDA核)；CPU Intel®Xeon®E5-2650v4@2.20 GHz．10.13245/j.hust.239410.T002表2异构并行加速效果分析案例网格精度/m网格量/103模拟时间/103sCPU时间/103sGPU时间/s加速比一维地形1.010.251180.0004.711 215.1312V型槽1.01 620.00012.60026.478 973.3361龙溪河4.57 978.00064.80010 080.0对于一维构造地形和V型槽案例，对比CPU计算用时，基于异构并行技术为模型计算带来了超过300倍的加速效果．在龙溪河山洪模拟案例中，由于在CPU端计算时间过长，未能统计CPU计算时间，仅在NVIDIA V100进行了计算，执行显示NVIDIA V100仅需2.8 h完成了18 h的事件模拟，验证了长历时计算时异构并行加速方案能显著提升计算效率，实现模拟的超前计算，在进一步实现基于动力学过程的山洪预报目标方面展现了巨大潜力．4 结论为减轻山洪灾害对人民群众的生命财产安全的威胁，本研究建立了从降雨到沟道径流的全物理过程的山洪动力学模型，采用二维扩散波模型对小流域山洪演进过程进行模拟，通过一阶迎风差分格式求解并基于异构(CPU+GPU)并行加速技术实现高性能计算，得出以下主要结论．a．建立了从降雨到沟道径流的山洪动力学模型，相比传统的降雨入渗径流过程物理模型增加并耦合了植被截留和沟道基流．对比模型实验和真实山洪案例表明：耦合植被截留和沟道基流的山洪动力学模型通过径流过程、洪峰流量及峰现时间3个指标均能更好地反映真实山洪事件演进过程；不考虑植被截留影响的径流偏大且洪峰出现时间更靠前；不考虑基流影响的模型径流偏小且峰现时间有延迟．新模型对真实山洪演进过程的揭示更全面和准确．b．通过一维构造地形、V型槽案例和龙溪河山洪过程案例分别验证模型求解和应用，显示模型求解的精度和应用的稳定性较高．建立了CPU+GPU异构并行技术的高性能求解计算方案，同时结合基于任务划分的异步多流并行方法和主机设备之间高效通信优化方法，使新构建的物理模型相对传统串行计算方案在计算效率上相对CPU单核提升300倍左右，表明构建的高性能动力学数值模拟系统在山洪预报领域具有潜在应用价值．