2020 , Vol. 22 >Issue 3: 505 - 515

“数字地形分析”专栏

• 吴钱娇 ,
• 陈玉敏 , *

• 武汉大学资源与环境科学学院,武汉 430079
* 陈玉敏（1977— ）,女,湖北武汉人,教授,主要从事数字地形分析和空间统计等方面的研究。E-mail:
 吴钱娇（1990— ）,女,安徽铜陵人,博士生,主要从事数字地形分析方面的研究。E-mail:qianjiaowu@whu.edu.cn

要求修回日期: 2019-12-11

网络出版日期: 2020-05-18

A Parallel Method of Surface Flow Dynamics Simulation based on CUDA

• WU Qianjiao ,
• CHEN Yumin , *
Expand
• School of Resource and Environment Science, Wuhan University, Wuhan 430079, China
* CHEN Yumin, E-mail:

Request revised date: 2019-12-11

Online published: 2020-05-18

Supported by

National Natural Science Foundation of China(41671380)

Abstract

Simulating the surface flow dynamics is of great importance in disaster prevention and mitigation, which could benefit the land remediation, regional planning, environmental protection and water resource management. Therefore, in this paper the Compute-Unified-Device-Architecture (CUDA) is embedded into the TIN-based surface flow dynamic simulation algorithm to get a parallel method for simulating the surface flow discharge. The aim of this study is to rapidly and accurately perform the surface flow dynamic simulation at any position and any time to meet the actual application requirements. First, the proposed algorithm extracts the critical points and drainage network from a high-precision Digital Elevation Model (DEM) to generate a drainage-constrained Triangulated Irregular Network (TIN). Second, the flow direction of each triangular facet over the TIN is calculated by the coordinate data of the triangular facets. Third, the flow path network is traced by the flow direction and rainfall source points. Fourth, the flow velocity of rain dropsover the flow paths is obtained by the flow velocity calculation kernel function based on the manning formula. Finally, by combing the flow velocity with the pre-set time, the algorithm can rapidly simulate the flow discharge at any location by usingthe flow discharge count kernel function. The specific paralle lization process consists of the transmission mode of data, partition model of thread and implement ation of the flow velocity calculation kernel function and flow discharge count kernel function. Data transmission in the paralle lization process includes the input and output of data. It inputs the data of the DEM, rainfall source points and flow path network from the CPU to GPU and outputs the data of the flow discharge calculated by the above two kernel functions from the GPU to CPU. The two kernel functions are parallelized by the flow paths. Each thread handles a single flow path. As a result, flow paths are divided by the partitioning method to obtain numerous thread blocks and the number of the thread in each thread block is allocated by the computing power of the GPU. The Black Brook Watershed (BBW) located in the north-eastern of New Brun swick was selected as the study area. To validate its accuracy and efficiency, the proposed method was compared with TIN-based surface flow dynamic simulation method. The experimental results demonstrate that the proposed algorithm can greatly improve the simulation efficiency while maintaining the accuracy simultaneously and its acceleration ratio can reach up to 11.2. In addition, the parallel algorithm was compared with the Soil and Water Assessment Tool (SWAT) model to verify its precision. The experimental results indicate that the Nash coefficient of the parallel method is increased by 88% and the correlation coefficient is increased by 5%.

1 引言

SHE模型是由英国、法国、丹麦的科学家联合研制出的一种典型的三维分布式物理模型。该模型为了便于模型参数与降雨数据的结合[10],在平面上将流域划分为规则格网,在垂向上将每个格网划分为若干水平层,以完成产流和汇流等水文过程的模拟。该模型的物理过程描述明确,但模拟效率较低。TOPMODEL模型是1979年由Beven和Kirkby[11]提出的一种基于规则格网DEM获取地形因子,从而实现流域汇流动态模拟的水文模型。该模型的模拟精度较高,但因为未考虑到降雨和蒸发,只适用于非复杂过程的水文模拟,并且该模型的数据分析效率较低。SWAT模型是1994年由Arnold[12]提出的一种基于规则格网DEM提取水文响应单元（HRUs）的半分布式水文模型。该模型的模拟效率较高,但模拟精度却受到格网DEM分辨率的限制[13,14]。LISFLOOD模型是2000年由Bates等[15]提出的一种利用非线性分布函数来表达每个格网与汇流的对应关系的降雨—汇流模型。该模型比较适合大流域的地表水动态模拟。DTVGM模型是2002年由夏军等[16]提出的一种分布式水文模型,该模型从DEM上获取坡度、坡向、水流路径等地形因子,以水量平衡方程为核心来实现流域产流和汇流过程的模拟。该模型降低其模拟精度对模型参数的依赖性,适用于信息不全流域的模拟,但该模型的普遍性还未得到验证。DR2-2013模型是2014年由 Lopez等[9]提出的通过改进径流深度与流域水温之间的相关性,从而能够适用于不同季节的一种分布式水文模型。该模型能够完成不同空间模式的水文模拟,但模拟精度较低。

2 基于TIN的地表水动态模拟方法与CUDA程序设计原理

2.1 基于TIN的地表水动态模拟方法

$v = R 2 3 ∙ S 1 2 / n$

2.2 CUDA程序设计原理

图1 CUDA线程组织模型

Fig. 1 Model of the threads by the CUDA

3 基于CUDA的地表水动态模拟并行 方法

3.1 基于TIN的地表水动态模拟并行化

3.1.1 数据的传输

图2 基于CUDA的并行化过程中的数据传输方式

Fig. 2 Transmission mode of data during the parallelization process based on CUDA

（1）根据降雨源点数据、DEM数据和流水线网络数据,基于曼宁公式,利用流水线流速计算核函数得到每条流水线上雨滴的流速;
（2）利用汇流量统计核函数得到该时刻的地表汇流量;
（3）循环计算,即可得到任意时刻的地表汇 流量。
3.1.2 CUDA中线程的划分

图3 基于CUDA的并行化过程中的线程划分模型

Fig. 3 Partition model of thread during the parallelization process based on CUDA

3.1.3 核函数的实现

Input: number of the flow lines (e.g. flineNum), the location of the flow line (e.g. flineTrackSize[]), the break points of the flow line (e.g. flowTrack[]),range of the DEM (e.g. xMin, Ymax), the resolution of the DEM (e.g. dx, dy), the width of the DEM (e.g. imgWidth), the runoff of the rainfall source points (e.g. rVal)
1: idx = blockIdx.x * blockDim.x + threadIdx.x
2: idy = threadIdx.y
3: id= idy*(gridDim.x * blockDim.x) + idx
4: ifid>flineNum
5: terminate
6:ifid > 0
7: Calculate the index of the flow line over flow track network (e.g. flineStart = flineTrackSize[id-1])
8: Calculate the start point of the flow line over flow track network (e.g. sP = flowTrack[flineStart)
9: Calculate the column and row of the start point (e.g. col = (sP.x-xMin)/dx, row = (yMax-sP.y)/dy)
10: Calculate the runoff of the start point (e.g. pRVal = rVal[row*imgWidth+col])
11:ifpRVal> 0
12:Calculate the flow velocity by the Manning’s equation
13: Calculate the discharge of the start point’s location after time interval (e.g. fVal[]=pRVal)
Obtain: the discharge of all of the start points after everytime interval(e.g. fVal[])

Input: number of the flow lines (e.g. flineNum), the location of the flow line (e.g. flineTrackSize[]),the width of height of the DEM (e.g. imgWidth, imgHeight),range of the DEM (e.g. xMin, Ymax), the resolution of the DEM (e.g. dx, dy), the width of the DEM (e.g. imgWidth), the discharge of all of the start pointsafter every time interval(e.g. fVal[])
1: idx = blockIdx.x * blockDim.x + threadIdx.x
2: idy = threadIdx.y
3: id = idy*(gridDim.x * blockDim.x) + idx
4: if id >flineNum
5: terminate
6: Calculate the index of the flow line over flow track network (e.g. flineStart = flineTrackSize[id-1])
7: Calculating the index of all of the raindrops (e.g. runP=flowTracks[flineStart+rainIdx[]])
8: Calculate the row and height of all the raindrops(e.g. col= (runP.x-xMin)/dx, row= (yMax-runP.y)/dy)
9: Calculate the discharge of all grid cells (e.g. pOutDis[]+=fVal[])
Output: the discharge of all grid cells (e.g. pOutDis[])

3.2 并行方法评价

$N = 1 - ∑ t = 1 n ( Q o t - Q m t ) 2 / ∑ t = 1 n ( Q o t - Q o ¯ ) 2$
$R = ∑ t = 1 n ( Q m t - Q m ¯ ) ( Q o t - Q o ¯ ) / ∑ t = 1 n ( Q m t - Q m ¯ ) 2 ∑ t = 1 n ( Q o t - Q o ¯ ) 2$
$B = ∑ t = 1 n Q o t / ∑ t = 1 n Q m t$
$S = T c / T p$

4 实验与结果分析

4.1 研究区域与实验环境

图4 实验中BBW流域的DEM影像

Fig. 4 Test DEM of the BBW

图5 图4子区域内不同阈值下的TIN

Fig. 5 TINs on the sub-region under different thresholds

图6 图4子区域内不同降雨源点尺度下的流水线网络（阈值为0.5 m）

Fig. 6 Flow pathnetworks on the sub-region under different scales of the rainfall points （threshold is 0.5m）

4.2 基于并行方法的地表水动态模拟实验

表1 BBW流域的不同阈值下TIN的精度结果

Tab.1 The accuracy of TIN under different thresholds on the BBW

TIN表面三角面的数量/个 23 348 11 079 6854 5139 4285
RMSE/m 0.17 0.31 0.45 0.58 0.73

4.3 基于并行方法的地表水动态模拟结果

表2 BBW流域出口水汇流量模拟精度的评价对比结果

Tab. 2 Accuracy comparison of simulating the flow discharge at the outlet of the BBW

10 15 20 25 30
0.5 基于TIN的地表

N 0.77 0.77 0.77 0.77 0.77
R 0.91 0.91 0.91 0.91 0.91
B 1.30 1.31 1.31 1.30 1.30

R 0.91 0.91 0.91 0.91 0.91
B 1.30 1.31 1.31 1.30 1.30
1.0 基于TIN的地表

N 0.77 0.77 0.77 0.77 0.77
R 0.91 0.91 0.91 0.91 0.91
B 1.31 1.31 1.31 1.31 1.31

R 0.91 0.91 0.91 0.91 0.91
B 1.31 1.31 1.31 1.31 1.31
1.5 基于TIN的地表

N 0.77 0.77 0.77 0.77 0.77
R 0.91 0.91 0.91 0.91 0.91
B 1.31 1.31 1.32 1.32 1.31

R 0.91 0.91 0.91 0.91 0.91
B 1.31 1.31 1.32 1.32 1.31
2.0 基于TIN的地表

N 0.77 0.77 0.77 0.77 0.77
R 0.91 0.91 0.91 0.91 0.91
B 1.31 1.32 1.31 1.32 1.31

R 0.91 0.91 0.91 0.91 0.91
B 1.31 1.32 1.31 1.32 1.31
2.5 基于TIN的地表

N 0.77 0.77 0.77 0.77 0.77
R 0.91 0.91 0.91 0.91 0.91
B 1.31 1.31 1.32 1.32 1.32

R 0.91 0.91 0.91 0.91 0.91
B 1.31 1.31 1.32 1.32 1.32

图7 BBW流域出水口模拟汇流量的评价结果对比

Fig. 7 Comparison of the simulated results at the BBW's outlet

表3 SWAT模型的精度评价结果（DEM为30 m）

Tab. 3 Statistical factors utilized to assess the precision of SWAT model （scale of DEM is 30 m）

N R B
SWAT 0.41 0.87 0.96

表4 BBW流域出口水汇流量模拟的运行时间统计结果

Tab. 4 Computation performance of simulating the flow discharge at the outlet of the BBW (s)

10 15 20 25 30
0.5 基于TIN的地表水动态模拟算法 953.82 425.37 278.18 213.77 169.77

1.0 基于TIN的地表水动态模拟算法 912.16 407.40 278.21 203.25 159.87

1.5 基于TIN的地表水动态模拟算法 899.07 395.99 267.96 196.33 159.60

2.0 基于TIN的地表水动态模拟算法 893.92 396.22 270.80 198.07 159.21

2.5 基于TIN的地表水动态模拟算法 892.59 400.38 260.40 199.60 164.56

4.4 模拟结果对比分析

（1）模拟效率的提高

图8 不同降雨源点尺度下的加速比

Fig. 8 Speed ratio under different scales of the rainfall points

（2）模拟精度的保持

图9 模拟结果的Nash系数分布

Fig. 9 Nash coefficient distribution map of simulation results

图10 模拟结果的相关系数分布

Fig. 10 Correlation coefficient distribution map of simulation results

图11 模拟结果的水量平衡系数分布

Fig. 11 Balance coefficient distribution map of simulation results

5 结论与讨论

（1）通过对基于TIN的地表水动态模拟算法进行分析发现,该算法中模拟任意时刻的地表汇流量耗时较多且可并行,因此,利用CUDA将该汇流量模拟功能改进成核函数,得到地表水动态模拟并行方法;
（2）在加拿大新布伦瑞克西北部的BBW流域,通过改变特征点提取阈值和降雨源点尺度,对比分析该算法与基于TIN的地表水动态模拟方法的模拟精度和效率;
（3）在该实验区域,将该算法与经典的SWAT模型进行精度的对比分析。

 [1] Beven K J, Warren R, Zaoui J . SHE: Towards a methodology for physically-based distributed forecasting in hydrology[J]. Hydrological Forecasting (IAHS Publication), 1980,129:133-137.
 [2] Abbott M B, Bathurst J C, Cunge J A , et al. An introduction to the European hydrological system-systemehydrologiqueEuropeen, ‘SHE’, 1: History and philosophy of a physically-based, distributed modelling system[J]. Journal of Hydrology, 1986,87:45-59.
 [3] Bathurst J C, Wicks J M, Connell P E. The SHE/SHESED basin scale water flow and sediment transport modelling system. In: V.P. Singh, eds. Computer models of watershed hydrology[J]. Highlands Ranch, CO: Water Resources Publication , 1995: 563-594.
 [4] Refsgaard J C, Storm B, Mike S, Singh. Computer models of watershed hydrology[C]. Highlands Ranch, CO: Water Resources Publication, 1995: 809-846.
 [5] Beven K J, Kirkby M J . A physically based, variable contributing area model of basin hydrology[J]. Hydrological Sciences Bulletin, 1979,24:43-69.
 [6] Arnold J G, Srinivasan R, Muttiah R R , et al. Large area hydrologic modeling and assessment part I: model development[J]. Journal of American Water Resources Association, 1998,34(1):73-89.
 [7] Bates P D, DE ROO A P J. A simple raster-based model for flood inundation simulation[J]. Journal of Hydrology, 2000,236(1-2):54-77.
 [8] 王纲胜, 夏军, 谈戈 , 等. 潮河流域时变增益分布式水循环模型研究[J]. 地理科学进展, 2002,21(6):573-582. [ Wang G S, Xia J, Tan G , et al. A research on distributed time variant gain model: A case study on Chaoheriver basin[J]. Progress in Geography, 2002,21(6):573-582. ]
 [9] Lopez-Vicente M, Perez-Bielsa C, Lopez-Montero T , et al. Runoff simulation with eight different flow accumulation algorithms Recommendations using a spatially distributed and open-source model[J]. Environmental Modelling & Software, 2014,62:11-21.
 [10] Chen Y, Zhou Q, Li S , et al. The simulation of surface flow dynamics using a flow-path network model[J]. International Journal of Geographical Information Science, 2014,28(11):2242-2260.
 [11] Romanowicz R, Beven K, Moore R . GIS and distributed hydrological models[J] //Mather P. Geographical information handling-research and applications Chichester[C]. UK: John Wiley & Sons, 1993: 197-205.
 [12] Shoemaker L, Dai T, Koenig J . TMDL Model Evaluation and Research Needs[M]. Newyork: The US Environmental Protection Agency, 2005.
 [13] Kumar N, Singh SK, Srivastava PK , et al. SWAT Model calibration and uncertainty analysis for streamflow prediction of the Tons River Basin, India, using Sequential Uncertainty Fitting (SUFI-2) algorithm[J]. Modeling Earth Systems and Environment, 2017,3(1):30.
 [14] Das B, Jain S, Singh S , et al. Evaluation of multisite performance of SWAT model in the Gomti River Basin, India[J]. Applied Water Science, 2019,9(5):134.
 [15] Van Der Knijff J M, Younis J, De Roo A P J. Lisflood: A GIS-based distributed model for river basin scale water balance and flood simulation[J]. International Journal of Geographical Information Science, 2010,24(2):189-212.
 [16] 蔡明勇 . DTVGM模型参数的遥感驱动研究与应用[D]. 北京:北京师范大学, 2011: 14-16. [ Cai M Y . Remote sensing based DTVGM and its application[D]. Beijing: Beijing Normal University, 2011: 14-16.]
 [17] Tucker G E, Lancaster S T, Gasparini N M , et al. The channel-hillslope integrated landscape development (CHILD) model[M]. In: Harmon R S, W Doe. Landscape erosion and evolution modelling. Dordrecht: Kluwer Academic/Plenum Publishers, 2001: 349-388.
 [18] Vivoni E R, Teles V, Ivanov V Y , et al. Embedding landscape processes into triangulated terrain models[J]. International Journal of Geographical Information Science, 2005,19(4):429-457.
 [19] Zhou Q, Pilesjö P, Chen Y . Estimating surface flow paths on a digital elevation model using a triangular facet network[J]. Water Resources Research, 2011,47(7):40.
 [20] Zhou Q, Chen Y . Generalization of DEM for terrain analysis using a compound method[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2011,66(1):38-45.
 [21] 陈玉敏, 吴钱娇, 巴倩倩 , 等. 多尺度地表水动态模拟及应用[J]. 测绘学报, 2015,44(S0):36-41. [ Chen Y M, Wu Q J, Ba Q Q , et al. Application research of multi-scale simulation of surface flowdynamics[J]. Acta Geodaetica et Cartographica Sinica, 2015,44(S0):36-41.]
 [22] Schenk O, Christen M, Burkhart H . Algorithmic performance studies on graphics processing units[J]. Journal of Parallel and Distributed Computing, 2008,68(10):1360-1369.
 [23] Mielikainen J, Huang B, Wang J , et al. Compute unified device architecture (CUDA)-based parallelization of WRF Kessler cloud microphysics scheme[J]. Computers and Geosciences, 2013,52:292-299.
 [24] NVIDIA Corporation. NVIDIA CUDA C Best Practice Guide[M]. Santa Clara, California: NVIDIA Corporation, 2010.
 [25] NVIDIA Corporation. NVIDIA CUDA C Programming Guide[M]. Santa Clara, California: NVIDIA Corporation, 2010.
 [26] Abouali M, Timmermans J, Castillo J E , et al. A high performance GPU implementation of Surface Energy Balance System (SEBS) based on CUDA-C[J]. Environmental Modelling and Software, 2013,41:134-138.
 [27] Xia X, Liang Q . A GPU-accelerated smoothed particle hydrodynamics (SPH) model for the shallow water equations[J]. Environmental Modelling and Software, 2016,75:28-43.
 [28] Zhu H, He L, Leeke M , et al. WolfGraph: The edge-centric graph processing on GPU[J]. Future Generation Computer Systems, 2019.
 [29] 王浩, 王含宇, 杨名宇 , 等. Retinex图像增强在GPU平台上的实现[J]. 地球信息科学学报, 2019,21(4):623-629. [ Wang H, Wang H Y, Yang M Y , et al. Implementation of Retineximage enhancement algorithm on GPU platform[J]. Journal of Geo-information science, 2019,21(4):623-629. ]
 [30] Pascal G, Pascal O, Camille D , et al. Ray-tracing of GNSS signal through the atmosphere powered by CUDA, HMPP and GPUs technologies[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014,7(5):1592-1602.
 [31] Chang F, Zhao Y, Liu S . CUDA parallel algorithm for CVA change detection of remote sensing imagery[J]. Journal of Remote Sensing, 2016,20(1):114-128.
 [32] 刘扬, 付征叶, 郑逢斌 . 高分辨率遥感影像目标分类与识别研究进展[J]. 地球信息科学学报, 2015,17(9):1080-1091. [ Liu Y, Fu Z Y, Zheng F B . Review on high resolution remote sensing image classification and recognition[J]. Journal of Geo-information science, 2015,17(9):1080-1091. ]
 [33] Qin C, Zhan L . Parallelizing flow-accumulation calculations on graphics processing units from iterative DEM preprocessing algorithm to recursive multiple-flow-direction algorithm[J]. Computers & Geosciences, 2012,43:7-16.
 [34] Ortega L, Rueda A J . Parallel drainage network computation on CUDA[J]. Computers & Geosciences, 2010,36:171-178.
 [35] Rueda A J, Jose M N, Adrian L . A comparison of native GPU computing versus OpenACC for implementing flow-routing algorithms in hydrological applications[J]. Computers & Geosciences, 2016,87:91-100.
 [36] 赵向辉, 苗青, 付忠良 , 等. 基于CUDA的汇流分析并行算法的研究与方法[J]. 计算机应用研究, 2010,27(7):2445-2451. [ Zhao X H, Miao Q, Fu Z L , et al. Research and realization in parallel algorithm of confluence analysis based on CUDA[J]. Application Research of Computers, 2016,87:91-100. ]
 [37] Carlotto T, da Silva R V, Grzybowski J M V . A GPGPU-accelerated implementation of groundwater flow model in unconfined aquifers for heterogeneous and anisotropic media[J]. Environmental Modelling and Software, 2018,101:64-72.
 [38] Wu Q, Chen Y, Wilson J P , et al. An effective parallelization algorithm for DEM generalization based on CUDA[J]. Environmental Modelling & Software, 2019,114:64-74.
 [39] O’Callaghan J F, Mark D M . The extraction of drainage networks from digital elevation data[J]. Computer Vision, Graphics and Image Processing, 1984,28:328-344.
 [40] Chen Z, Guevara J. Systematic selection of very important points (VIP) fromdigital terrain model for constructing triangular irregular networks[C]. In: Proc.8th International Symposium on Computer-Assisted Cartography, AutoCarto 8, Baltimore, MD, 1987.
 [41] Chang K . Introduction to Geographic Information Systems, 4th ed[M]. New York: McGraw-Hill, 2007.
 [42] 国家测绘局测绘标准化研究所. CH/T1008—2001基础地理信息数字产品1:10 000, 1:50 000数字高程模型[S]. 北京:国家测绘局, 2001. [ State Bureau of Surveying and Mapping. 1:10000 and 1:50000 digital elevation model. The Basic Digital Geographical Information Product, Standard CH/T 1008-2001[S]. Beijing: The State Bureau of Surveying and Mapping of China, 2001. ]
 [43] 刘学军, 卢华兴, 卞璐 , 等. 基于DEM的河网提取算法的比较[J]. 水利学报, 2006,37(9):1134-1141. [ Liu X J, Lu H X, Bian L , et al. Comparison of algorithms for extracting drainage network from grid-based digital elevation model[J]. Journal of Hydraulic Engineering, 2006,37(9):1134-1141. ]
 [44] 金鑫, 金彦香, 杨登兴 . SWAT模型在土地利用/覆被变化剧烈地区的改进与应用[J]. 地球信息科学学报, 2018,20(8):1064-1073. [ Jin X, Jin Y X, Yang D X . Improved SWAT and its application to a region with severe land use/land cover changes[J]. Journal of Geo-information science, 2018,20(8):1064-1073. ]
 [45] Thavhana M P, Savage M J, Moeletsi M E . SWAT model uncertainty analysis, calibration and validation for runoff simulation in the Luvuvhu River catchment, South Africa[J]. Physics and Chemistry of the Earth, 2018,105:115-124.
Options

/

 〈 〉