基于GPU并行加速的水库溃坝洪水演进快速模拟方法

基于gpu并行加速的水库溃坝洪水演进快速模拟方法
技术领域
1.本发明属于数据处理及洪涝预测领域，具体公开了一种基于gpu并行加速的水库溃坝洪水演进快速模拟方法。


背景技术：

2.水库大坝具有防洪、灌溉、发电等功能，对推动社会经济发展起着重要的作用。随着城镇化进程的不断加快与扩张，一些原在郊区的水库逐渐或已进入城区，这些位于生命财产集中度更高的城区内的水库大坝，更是重点警惕的对象。溃坝虽然是一种低频率、高风险的灾害，但是大坝一旦发生溃决，将对下游地区带来灾难性影响，因此必须高度重视并采取相应的预防和处理措施。大坝因设计或施工不当，存在危险隐患，在隐患发展到一定程度后就可能产生溃坝。大坝除了质量问题导致溃坝外，地震、战争、恐怖袭击等不可抗力因素以及人为破坏也有可能导致大坝溃决。在这样的背景下，对溃坝后洪水下泄可能造成的破坏性影响和淹没情况进行模拟和预测，对下游可能被淹的地区制订适当的防护措施和紧急撤离预案是非常重要的。
3.自上世纪50年代以来，随着计算机运算能力不断提高以及数值计算技术进一步发展，数值模拟逐渐成为研究溃坝洪水演进的主要手段。数值模拟包括一维、二维以及三维三种方式，现今，一二维模拟的方式比较成熟，而二维模拟克服了一维模拟只能模拟几个点的缺陷，具有可以模拟较为真实地形、能计算平面参数特征、边界处理方便、结合gis等技术可视化效果较好、方便直观显示模拟结果、计算结果精确等优点而被大量应用。水库溃坝洪水的演进可以采取数值模拟方法进行预测，相对一维模拟，二维模拟能提供更丰富的淹没效果分析；相对三维模拟，二维模拟所需的计算时间更少，研究范围更广。因此二维模拟是当前溃坝洪水演进的一种常用的数值模拟方法。实施二维模拟的工具包括mike21、infoworks icm、lisflood-fp模型、swm模型、floodmap、二阶 muscl 模型等，它们一般通过求解完整或者简化的浅水方程模拟洪水的演进过程，但直接计算需要消耗大量时间，难以满足溃坝洪水预测和模拟的时效性,且商业化性质限制了像是mike21和infoworks icm的推广以及二次开发潜力。与传统二维水动力模型不同，wca2d（weighted cellular automata 2d inundation model）模型是由michele guidolin与albert s. chen等人开发，是一种基于权重转换规则的元胞自动机模型，尽管没有考虑惯性项和动量守恒，但具备模拟二维水动力过程的能力，并且模拟效率很高；该模型能够用来处理多种网格( 例如矩形网格、六边形网格或三角形网格) 以及不同的元胞邻域类型( 例如四元胞的冯-诺依曼邻域或八元胞的摩尔邻域) ，而模型的gpu并行加速技术目前也已经取得突破。在水库溃决这类突发性的灾难面前，时间就是生命，数值模拟的计算效率是影响灾害防御决策的关键因素。
4.鉴于此，采取全开源的wca2d模型对水库溃坝洪水的演进开展数值模拟理论上是可行的，是能满足其时效性与推广要求的。然而，目前应用wca2d模型对水库溃坝洪水演进实施快速模拟的相关技术尚未见报到，因此在该模型基础上发明一种基于gpu并行加速的水库溃坝洪水演进快速模拟方法显得十分迫切。


技术实现要素：

5.本发明的目的在于提供一种可实现水库溃坝洪水演进二维快速模拟方法的方法，为城区溃坝洪水演进模拟提供一种操作简单、计算精度高、计算效率高、可视化效果好的工具。
6.本发明至少通过如下技术方案之一实现。
7.基于gpu并行加速的水库溃坝洪水演进快速模拟方法，包括以下步骤：s1、收集水库及其下游的基础数据，s2、对地形与土地利用数据进行预处理，建立数字高程模型；s3、选定工况，并根据该工况对应的水库库容、溃决水深、溃口宽度及溃决方式推求溃口流量过程；s4、在使用图形管理器（gpu）并行计算的wca2d模型中输入步骤s2以及步骤s3中处理好的地形和流量过程数据文件，并设置wca2d模型的模拟参数后进行模拟；s5、将模拟结果进行水量平衡、流态分析，并验证、率定、优化数值模型，根据水量平衡、流态分析的结果与经验值进行对比分析计算结果可靠性和合理性；s6、输出wca2d模型模拟结果，并基于地理信息系统（gis）的空间分析工具和r语言对结果进行处理，实现溃坝洪水演进过程的可视化。
8.优选的，所述基础数据不限于水库库容曲线、水库大坝相关数据、地形数据、土地利用数据；优选的，输出的模拟结果包括溃坝后下游分时的洪水流速、淹没水深及淹没范围以及整个洪水演进过程中下游淹没最大范围的最大水深。
9.优选的，所述所述预处理不限于对地形数据进行修正并叠加建筑和河道数据。
10.优选的，步骤s2包括基于地理信息系统，并结合土地利用数据、卫星图对建筑物、河道、街道、公路的地形，建立包括建筑物的形状、实际高程的数字高程模型，以模拟建筑物对溃坝洪水的延缓作用以及因建筑物之间的空间与间隙的涌水作用。
11.优选的，水库溃决方式包括大坝瞬间溃决和逐渐溃决两种溃坝形式。
12.优选的，所述瞬间溃决根据溃口宽度不同，分为瞬间全部溃决和瞬间部分溃决，以模拟包括因恐怖袭击、地震的突发事件造成的水库溃决洪水演进过程；逐渐溃决模拟溃口由小到大的发育发展过程。
13.优选的，步骤s3包括以下步骤：s31、根据要模拟工况，选定对应的参数，包括水库溃决方式、溃口宽度、溃决水深及对应水库库容；s32、计算最大溃决流量，其中，瞬间溃决的最大溃决流量为：式中q
max
、g、b 、bm、h0分别表示大坝溃决最大流量、重力加速度、坝长、最终溃口宽度、溃决水深；逐渐溃决的大坝先渗透变形破坏再瞬间溃塌，渗透变形破坏流量为：
式中：h为水库水位，即水库计算水位高程；a为水流所流过的管道的断面面积；h
p
为管道中心线高程；f 为达西摩擦系数，取决于d50粒径，可以由moody 曲线计算，l 为管道沿水流方向的长度；d为管道的直径或宽度；s33、溃口洪水流量过程假设大坝最大溃决流量为q
max
，溃决洪量为已知的水库库容w，洪水历时tn待求，则t 时刻流量q为：。
14.优选的，所述wca2d模型用于模拟地表水流，主要过程包括：（1）获取相邻元胞水流交换的权重：（1）获取相邻元胞水流交换的权重：（1）获取相邻元胞水流交换的权重：式中：m表示与中心元胞相邻的元胞总数；表示编号为0的中心元胞与第i个相邻元胞的水位差值；表示相应的体积差值；表示水位差值大于临界值时才计算下游元胞的获得水量；表示第i个相邻元胞的权重；表示中心元胞的权重；表示与中心元胞相邻的体积差值总和；表示与中心元胞相邻的元胞最小的体积差值；（2）获取元胞间交互水量：式中：表示时刻下游第i个元胞获得的交换水量；表示在时刻从中心元胞转移的水量；表示中心元胞面积；表示中心元胞水深；表示权重；表示中心元胞与该元胞的水位差值；表示中心元胞中点到该元胞中点的距离；表示中心元胞水流转移到该元胞的最大允许流速；表示时间步长，s；表示该元胞与中心元胞的交界长度；表示t时刻从中心元胞转移的水量；n为曼宁系数；g为重力加速度；下标为m的变量都是针对拥有最大权重的元胞；所有不含上标的变量都假定在t时刻；（3）计算水深：
式中：表示时刻中心元胞的水深；表示t时刻中心元胞水深；表示中心元胞获得的流入水量；表示中心元胞流出的水量；（4）流速计算公式为：式中：表示时刻中心元胞水流转移到第i个元胞的流速；表示时刻中心元胞和第i个元胞水深的算术平均值；表示第i个元胞与中心元胞的交界长度。
15.优选的，采用opencl异构算法优化wca2d模型，获得gpu加速版本的wca2d模型，以提高对水库溃坝洪水演进的模拟效率。
16.优选的，所述基础数据不限于水库库容曲线、水库大坝相关数据、地形数据、土地利用数据。
17.优选的，所述模拟参数包括初始条件、边界条件和糙率、步长、迭代次数、模拟时间。
18.基于arggis 平台，实现洪水演进的可视化，将大量的计算结果数据以溃坝洪水的动态演示形式表现，利于从视角角度上直观了解水流动态特性。叠加卫星图，直观清晰地呈现分时下游地区的淹没范围及淹没深度。
19.与现有的技术相比，本发明的有益效果为：wca2d模型是一种基于权重转换的元胞自动机模型，模拟效率很高，弥补了传统二维水动力模型运行效率低下的不足。此外，wca2d模型可以实现并行gpu计算的功能，可以进一步提高模拟的效率。模型输出的下游分时淹没水深、淹没范围以及水流速度结果，可在gis平台很好的完成可视化，方便对数据进行后处理，对于制定应急预案减少损失具有重大的意义。此模型全开源，利于推广应用以及二次开发。可基于arcgis 平台，实现洪水演进的可视化，方便对数据进行后处理，对于制定应急预案减少损失具有重大的意义。
附图说明
20.图1为本发明实施例基于gpu并行加速的水库溃坝洪水演进快速模拟方法流程图；图2为本发明实施例一洪水流量过程线图；图3为本发明实施例一工况1下lisflood、wca2d模型、wca2d并行cpu加速模型、wca2d并行gpu加速模型最大水深示意图；图4为本发明实施例一工况2下lisflood、wca2d模型、wca2d并行cpu加速模型、wca2d并行gpu加速模型最大水深示意图；图5为本发明实施例一工况3下lisflood、wca2d模型、wca2d并行cpu加速模型、wca2d并行gpu加速模型最大水深示意图；
图6为本发明实施例二洪水流量过程线图；图7为本发明实施例二wca2d并行gpu加速模型最大水深示意图；图8为本发明实施例三洪水流量过程线图；图9为本发明实施例三wca2d并行gpu加速模型最大水深示意图；图10为本发明实施例四洪水流量过程线图；图11为本发明实施例四wca2d并行gpu加速模型最大水深示意图。
具体实施方式
21.为了使本技术领域的人员更好地理解本发明方案，下面结合附图和具体实施方式对本发明作进一步的详细说明。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
22.实施例一，以广东省廉江市江头水库为例。
23.如图1所示，本发明公开了一种基于gpu并行加速的水库溃坝洪水演进快速模拟方法，包括以下的步骤：s1、收集水库及其下游的基础数据，包括但不限于水库库容曲线、水库大坝相关数据、地形数据、土地利用数据。
24.本具体实施例以广东省廉江市江头水库为例，其水库坝长175m，坝底高程21m，正常蓄水位42.65m，对应库容825万m3，校核洪水位41.92m，对应库容为1172万m3，漫顶水位42.65m，对应库容为1215万m3。
25.s2、对地形与土地利用数据进行预处理，包括但不限于对地形数据进行修正并叠加建筑和河道数据，转化为模型能用的格式。基于arcgis，并结合土地利用数据、卫星图等对建筑物、河道、街道、公路等地形对原始数字高程模型(dem)进行精确处理，处理后的dem可以精确把建筑物的形状、实际高程表现出来，以精确模拟建筑物对溃坝洪水的延缓作用以及因建筑物之间的空间与间隙的涌水作用。
26.s3、选定工况，并根据该工况对应的水库库容、溃决水深、溃口宽度及溃决方式，通过经验公式推求溃口流量过程。
27.s31、根据要模拟工况，选定对应的参数，包括水库溃决方式（包含大坝瞬间溃决和逐渐溃决两种方式）、溃口宽度、溃决水深及对应水库库容。
28.水库溃决方式同时考虑了水库瞬间溃决和逐渐溃决两种溃坝形式。瞬间溃决可以根据溃口宽度不同，分为瞬间全部溃决或是瞬间部分溃决，可以模拟因恐怖袭击、地震等突发事件造成的水库溃决洪水演进过程，是目前大多数溃坝计算的模拟基础。而大坝溃口逐渐溃决模型可较真实模拟溃口由小到大的发育发展过程，与大部分溃坝的实际情况相符。逐渐溃决情景下，大坝先渗透变形破坏，当渗透变形破坏断面发展到一定大小再瞬间溃塌，即先从小孔逐渐发展成为大孔，然后再瞬间溃塌。
29.在本具体实施例，选择三个工况如表1所示，溃口宽度都设置为175m，即整个大坝全溃，溃坝方式选择为瞬间全溃。三种工况分别采用正常蓄水位、校核水位以及漫顶水位及其对应的库容。
30.s32、计算最大溃决流量。其中，瞬间溃决的最大溃决流量计算方式为：式中q
max
、g、b 、bm、h0分别表示大坝溃决最大流量(m3/s)、重力加速度(m/s
2 )、坝长(m)、最终溃口宽度(m)、溃决水深(m)。
31.s33、溃口洪水流量过程假设洪峰为qmax(m3/s)，溃决洪量为已知的水库库容w(m3)，洪水历时tn(s) 待求，则t 时刻流量 q 的公式：最后得到三种工况下的溃口流量过程线如图2所示。
32.s4、在使用并行gpu计算的wca2d模型版本中输入处理好的地形和流量过程等数据文件，并设置模型的初始条件、边界条件和其他模拟参数后进行模拟。
33.wca2d（weighted cellular automata 2d inundation model）模型是一种基于权重转换规则的元胞自动机模型。与传统二维模拟模型不同，该模型没有考虑惯性项和动量守恒，故而模拟二维水动力过程的效率很高。模型能够用来处理多种网格以及不同的元胞邻域类型。 wca2d模拟地表水流的主要计算过程包括：
①
相邻元胞水流交换的权重计算；
②
元胞间交互水量计算；
③
水深计算；
④
流速计算。 计算公式如下，所有不含上标的变量都假定在t时刻：（1）权重计算公式为：（1）权重计算公式为：式中：m表示与中心元胞相邻的元胞总数；表示中心元胞（编号为0）与第i个相邻元胞的水位差值，m；表示相应的体积差值，m
³
；表示水位差值大于临界值时才计算下游元胞的获得水量；表示第i个相邻元胞的权重；表示中心元胞的权重。
34.（2）元胞间交互水量：
式中：表示时刻下游第i个元胞获得的交换水量，m
³
，表示在时刻从中心元胞转移的水量，m
³
；表示中心元胞面积，m2；表示中心元胞水深，m；表示权重；表示中心元胞与该元胞的水位差值；表示中心元胞中点到该元胞中点的距离，m；表示中心元胞水流转移到该元胞的最大允许流速，m/s；表示时间步长，s；表示该元胞与中心元胞的交界长度，m；表示t时刻从中心元胞转移的水量，m
³
；n为曼宁系数；g为重力加速度。所有下标为m的变量都是针对拥有最大权重的元胞。
35.（3）水深计算公式为
ꢀꢀꢀꢀꢀ
（5）式中：表示时刻中心元胞的水深；表示t时刻中心元胞水深；表示中心元胞获得的流入水量，比如降雨，m
³
；表示中心元胞流出的水量，比如下渗水量，m
³
。
36.（4）流速计算公式为式中：表示时刻中心元胞水流转移到第i个元胞的流速，m/s；表示时刻中心元胞和第i个元胞水深的算术平均值，m；表示第i个元胞与中心元胞的交界长度，m。
37.在此案例中，采用opencl异构算法优化wca2d模型，获得gpu加速版本的wca2d模型，以提高对水库溃坝洪水演进的模拟效率。opencl可以实现异构系统（一般指cpu与gpu（显卡））之间的通讯与数据传送，能够开辟、释放和初始化设备中的存储空间，能够实现主机和设备间的数据的传输。
38.为深入分析基于gpu加速的wca2d模型在水库溃坝洪水演进的模拟效率，将基于openmp算法（cpu加速）的wca2d模型、未进行加速的wca2d模型以及二维水动力模型lisflood的模拟效率作为对比。
39.s5、将数值模拟的结果进行水量平衡、流态分析，和多个模型进行对比，验证、率定、优化数值模型，分析计算结果可靠性和合理性。
40.在本具体实例中，通过对比wca2d的gpu并行加速版与cpu并行加速版、普通版本以
及二维水动力模型lisflood版本（自带cpu加速）的结果进行相互验证，验证结果如图3a-d、图4a-d、图5a-d所示。结果表明，四种模型的结果在淹没范围以及淹没深度的分布上（图3a-d、图4a-d、图5a-d、表2）近似，可相互印证。其中，四种模型的运行时间如表3所示（所用设备配置为：intel(r) core(tm) i7-8700 cpu），在进行gpu并行加速前，wca2d的模拟效率比lisflood差，但在进行gpu加速后，效率得到了极大的提升，运行速度快于lisflood的模拟速度。
41.表2 lisflood模型、ca模型以及ca加速模型下三种工况水深面积分布(km
²
)表3 三种模型下各工况模拟时长（秒）s6、输出并处理wca2d模型模拟的结果，包含溃坝后下游分时的洪水流速、淹没水深及淹没范围以及整个洪水演进过程中下游淹没最大范围的最大水深。
42.基于gis平台，实现洪水演进的可视化，将大量的计算结果数据以溃坝洪水的动态演示形式表现，利于从视角角度上直观了解水流动态特性。叠加卫星图，直观清晰地呈现分时下游地区的淹没范围及淹没深度。
43.实施例二，以广东省廉江市长青水库为例。本实施例中的步骤与实施例一基本相同，不同之处在于：s1、收集水库及其下游的基础数据，包括但不限于水库库容曲线、水库大坝相关数据、地形数据、土地利用数据。
44.本具体实施例以广东省廉江市长青水库为例，其主要水库岭北下水库坝顶高程50m，校核水位48.92m，校核库容12450万方，汛限水位为45m，汛限库容为7100万方，坝长525m。
45.s3、选定工况，并根据该工况对应的水库库容、溃决水深、溃口宽度及溃决方式，通过经验公式推求溃口流量过程。
46.s31、根据要模拟工况，选定对应的参数，包括水库溃决方式（包含大坝瞬间溃决和逐渐溃决两种方式）、溃口宽度、溃决水深及对应水库库容。
47.在本具体实施例，选择一个工况，溃决方式为瞬间全溃。其中漫顶水位50.5m，溃口溃决宽度为414m，溃到底，溃口为梯形。溃口流量过程如图6所示。
48.本实施例二的受灾区受淹没区域的最大淹没水深及范围图如图7所示。
49.实施例三，以广东省廉江市武陵水库为例。本实施例中的步骤与实施例一基本相同，不同之处在于：s1、收集水库及其下游的基础数据，包括但不限于水库库容曲线、水库大坝相关数据、地形数据、土地利用数据。
50.本具体实施例以广东省廉江市武陵水库为例，水库坝顶高程42m，校核水位41.26m，校核库容9930万方，汛限水位为38m，汛限库容为6380万方，坝长150m。
51.s3、选定工况，并根据该工况对应的水库库容、溃决水深、溃口宽度及溃决方式，通过经验公式推求溃口流量过程。
52.s31、根据要模拟工况，选定对应的参数，包括水库溃决方式（包含大坝瞬间溃决和逐渐溃决两种方式）、溃口宽度、溃决水深及对应水库库容。
53.在本具体实施例，选择一个工况，溃决方式为逐渐溃坝再瞬间溃坝。其中逐渐溃坝，渗透变形破坏为矩形，初始渗透变形破坏长度为0.1m，校核洪水位42.60m，渗透变形破坏高程30.95m，溃口线性发展，溃决宽度为75m。溃口流量过程如图8所示。
54.本实施例三的受灾区受淹没区域的最大淹没水深及范围图如图9所示。
55.实施例四，以广东省广州市龙洞水库为例。本实施例中的步骤与实施例一基本相同，不同之处在于：s1、收集水库及其下游的基础数据，包括但不限于水库库容曲线、水库大坝相关数据、地形数据、土地利用数据。
56.本具体实施例以广东省广州市龙洞水库为例。水库集雨面积6.36km2，正常蓄水位为64.86m(珠基百年一遇设计洪水位为64.90m，千年一遇校核洪水位为66.10m），正常库容237万m3,总库容250万m3。水库坝长164m,坝顶高程67.00m，防浪墙高0.5m,最大坝高为22.00m。
57.s3、选定工况，并根据该工况对应的水库库容、溃决水深、溃口宽度及溃决方式，通过经验公式推求溃口流量过程。
58.s31、根据要模拟工况，选定对应的参数，包括水库溃决方式（包含大坝瞬间溃决和
逐渐溃决两种方式）、溃口宽度、溃决水深及对应水库库容。
59.在本具体实施例，选取龙洞水库在校核水位为66.1m时由地震、滑坡等突发状况引起主坝瞬间溃决，溃口宽度为大坝长度164 m，水库库容250万m3。溃口流量过程如图10所示。
60.本实施例四的受灾区受淹没区域的最大淹没水深及范围图如图11所示。
61.以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为所述的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。