A parallel analysis method for water flow direction and cumulative flow based on CUDA architecture D8 algorithm

Abstract

The present invention provides a parallel analysis method for watershed flow direction and cumulative flow using the D8 algorithm in a CUDA architecture. The method comprises: constructing a read and storage model for watershed DEM raster data based on GPU texture memory; setting parallel grids and thread blocks based on GPU global memory, mapping threads to DEM raster cells, and obtaining the current thread identification number; each thread continuously picking the elevation values of the center grid and adjacent grids of the DEM raster data according to the read and storage model; calculating the maximum distance-weighted drop based on the elevation values of the center grid and adjacent grids, and determining the flow direction of the center grid based on the maximum distance-weighted drop; constructing a flow direction grid cell based on the flow direction of the center grid, mapping threads to flow direction grid cells, and reading thread identification numbers and flow directions; and each thread performing serial flow calculations according to the flow direction to obtain a cumulative flow matrix. The present invention proposes a parallel method for fast computation of the D8 algorithm in a CUDA environment, solving the problems of raster pixel access and computation conflicts in parallel computation of the D8 parallel algorithm.

Description

Technical Field

[0001] The present invention relates to the technical field of confluence analysis, and in particular to a parallel analysis method for water area flow direction and cumulative confluence volume based on a CUDA architecture D8 algorithm. Background Art

[0002] Slope hydrological and hydrodynamic processes are the carriers of associated processes such as floods, soil erosion, and non-point source pollution. The process mechanisms and simulations are widely used in the early warning, forecasting, and response to water disasters. The flow direction and flow algorithm is a basic algorithm for slope hydrological and hydrodynamic simulation. It is used to quantify flow direction and flow and water systems, hydrological response unit division and runoff simulation, flood disaster-pregnancy processes and their temporal and spatial distribution patterns, sediment transport ratio and its LUCC effect, and other processes. The flow direction and flow algorithm mainly includes single-flow method and multi-flow method. The single-flow algorithm means that the grid only flows to the surrounding lowest elevation grid, mainly including the D8 algorithm, Rho4 / Rho8 algorithm, and Lea algorithm; the multi-flow algorithm means that the grid distributes flow to the downstream grid according to the diversion ratio, mainly including the DEMON algorithm, Dinf algorithm, and MFD algorithm. Among them, the D8 algorithm is more widely used.

[0003] The velocity of water flow on slopes is low and varies little, and the primary driving force behind the flow direction is potential energy. Single-direction and multi-direction algorithms use different flow direction selection methods, but water flows in the direction of decreasing potential energy, that is, toward surrounding grids with lower elevations. Calculating cumulative flow involves inter-grid operations, structurally presenting a serial operation from the water source to the boundary. Calculating flow direction and cumulative flow, grid-by-grid from the DEM raster data boundary, is time-consuming and inefficient. Some existing studies have attempted to address data reading, transmission, and storage issues within the CUDA architecture. However, when the D8 algorithm is run in parallel, grid calculations can conflict, and an effective solution has yet to be found.

[0004] To this end, the present invention proposes a new parallel analysis method for water area flow direction and cumulative flow based on CUDA architecture D8 algorithm. Summary of the Invention

[0005] To solve the above problems, the present invention proposes a parallel method for fast operation of the D8 algorithm in the CUDA environment, which solves the problem of raster operation conflicts in the parallel operation of the D8 parallel algorithm.

[0006] To achieve the above objectives, the present invention provides the following technical solutions.

[0007] A parallel analysis method for water flow direction and cumulative flow using a CUDA architecture D8 algorithm includes the following steps:

[0008] Construct a reading and storage model for watershed DEM raster data based on GPU texture memory;

[0009] Based on the GPU global memory, set the parallel grid and thread block, map the thread to the DEM grid cell, and obtain the current thread identification number;

[0010] Each thread continuously picks up the central grid and adjacent grid elevation values of the collected watershed DEM raster data according to the read storage model;

[0011] The maximum distance weighted drop is calculated based on the elevation values of the central grid and the adjacent grids, and the flow direction of the central grid is determined based on the maximum distance weighted drop;

[0012] Construct flow direction grid cells according to the center grid flow direction, map threads to flow direction grid cells, and read thread identification numbers and flow directions;

[0013] Each thread performs serial confluence calculations according to the flow direction. Access and operation conflicts on the same grid are resolved through the CUDA atomic add function to obtain the cumulative confluence matrix.

[0014] Preferably, the construction of the read storage model includes the following steps:

[0015] Declare a two-dimensional CUDA array and copy the watershed DEM raster data from the GPU global memory to the array;

[0016] Set the return data type to a two-dimensional floating-point global variable;

[0017] Set the non-normalized texture coordinate reading mode and boundary and picking method properties to determine whether the texture coordinates are out of bounds;

[0018] Bind the CUDA array to a texture reference, and fetch the elevation values of the center grid and its adjacent grids from the texture.

[0019] Preferably, the method further comprises: preprocessing the watershed DEM grid data, and preprocessing the DEM grid data of the filled flat land respectively.

[0020] Preferably, determining the central grid flow direction according to the maximum distance weighted drop specifically includes the following steps:

[0021] Get the current thread identification number and map the thread to the DEM grid cell;

[0022] If the current thread identification number is less than the total number of DEM grids, parallel calculation will begin;

[0023] Calculate the two-dimensional coordinates of the central grid unit and continuously pick the elevation values of the central grid and adjacent grids;

[0024] The maximum grid distance weight drop is calculated based on the elevation values of the central grid and adjacent grids to determine the flow direction of the central grid;

[0025] If the current thread identification number is greater than the total number of DEM grids, the parallel calculation of this thread ends.

[0026] Preferably, each thread continuously traverses and calculates through the CUDA atomic add function to obtain the cumulative sink flow, specifically including:

[0027] Get the current thread identification number and map the thread to the flow grid unit;

[0028] If the current thread identification number is less than the total number of traffic grids, parallel traffic calculation begins;

[0029] Calculate the two-dimensional coordinates of the central grid unit and continuously pick the central grid flow direction;

[0030] According to the flow direction of the central grid, the CUDA atomic add function is used to traverse and calculate the current grid cumulative flow;

[0031] If the picked flow direction grid is out of bounds, the loop calculation ends;

[0032] If the current thread identification number is greater than the total number of traffic grids, the parallel calculation of this thread ends.

[0033] Preferably, the method further comprises: counting the river lengths below the flow threshold according to the cumulative flow matrix, and calculating the river network density and the drainage area.

[0034] Beneficial effects of the present invention: The present invention proposes a parallel method for fast calculation of the D8 algorithm under the CUDA environment, which effectively improves the parallel calculation rate of flow direction and cumulative flow, and solves the problems of parallel calculation of raster pixel access and calculation conflicts in the D8 parallel algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 is a flow chart of the overall method of an embodiment of the present invention;

[0036] Figure 2 This is a parallel flow diagram of the D8 algorithm according to an embodiment of the present invention;

[0037] Figure 3 This is a parallel graph of D8 algorithm traffic according to an embodiment of the present invention;

[0038] Figure 4 The Ganjiang River Basin water system and boundaries under the D8 classic and parallel algorithms of the embodiment of the present invention;

[0039] Figure 5 It is the CUDA parallel speedup response of the Ganjiang River Basin in the embodiment of the present invention;

[0040] Figure 6 This is a curve showing the marginal change of the acceleration ratio with the number of grids under the Ganjiang River Basin and other thread blocks in an embodiment of the present invention. DETAILED DESCRIPTION

[0041] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0042] Example 1

[0043] The present invention proposes a parallel analysis method of water flow direction and cumulative flow based on CUDA architecture D8 algorithm. Figure 1 As shown, the specific steps include:

[0044] S1: Build a reading and storage model for watershed DEM raster data based on GPU texture memory.

[0045] S2: Based on the GPU global memory, set the parallel grid and thread block, map the thread to the DEM grid cell, and obtain the current thread identification number.

[0046] S3: Each thread continuously picks up the elevation values of the center grid and adjacent grids of the collected watershed DEM grid data according to the read storage model; calculates the maximum distance weight drop based on the elevation values of the center grid and adjacent grids, and determines the flow direction of the center grid based on the maximum distance weight drop. The process is as follows Figure 2 shown.

[0047] Specifically:

[0048] S3.1: Get the current thread identification number and map the thread to the DEM grid cell;

[0049] S3.2: If the current thread identification number is less than the total number of DEM grids, start parallel computing;

[0050] S3.3: Calculate the two-dimensional coordinates of the central grid cell and continuously pick the elevation values of the central grid and adjacent grids;

[0051] S3.4: Calculate the maximum grid distance weight drop based on the elevation values of the central grid and adjacent grids to determine the flow direction of the central grid;

[0052] S3.5: If the current thread identification number is greater than the total number of DEM grids, the parallel calculation of this thread ends.

[0053] S4: Construct flow grid cells based on the center grid flow direction, map threads to flow grid cells, and read thread identification numbers and flow directions. Each thread performs serial flow calculations according to the flow direction. Conflicts in accessing and calculating the same grid are calculated using the CUDA atomic add function to obtain the cumulative flow matrix. Figure 3 shown.

[0054] Specifically:

[0055] S4.1: Get the current thread identification number and map the thread to the flow grid unit;

[0056] S4.2 If the current thread identification number is less than the total number of traffic grids, start parallel traffic calculation;

[0057] S4.3 calculates the two-dimensional coordinates of the central grid unit and continuously picks the central grid flow direction;

[0058] S4.4 calculates the current grid cumulative flow according to the flow direction of the central grid and uses the CUDA atomic add function to traverse and calculate the current grid cumulative flow;

[0059] S4.5 If the picked flow direction grid is out of bounds, the loop calculation ends;

[0060] S4.6 If the current thread identification number is greater than the total number of traffic grids, the parallel calculation of this thread ends.

[0061] In this embodiment:

[0062] (1) Evaluation of water system extraction effect of D8 parallel algorithm under CUDA environment:

[0063] Based on the DEM data of the Ganjiang River Basin, the Wan'an Hydropower Station, Ji'an and Waizhou Hydrological Stations were selected as outlet sections, and the cumulative runoff threshold was set to 300,000. The ArcGIS hydrological analysis tool (AG), Matlab software D8 classic algorithm (ML) and CUDA-supported D8 parallel algorithm were used to extract the main tributary river systems and boundaries of the basin ( Figure 4 ). By counting the river length, river network density and basin area under the convergence threshold, it can be seen that the length of the upper Ganjiang River basin under ArcGIS hydrological analysis, Matlab software D8 classic algorithm and CUDA parallel algorithm are 1991km, 1992km and 1994km respectively, and the maximum value relative to the minimum value change rate is 0.2%; the length of the middle and lower Ganjiang River basin is 3288km, 3277km and 3278km respectively, and the maximum value relative to the minimum value change rate is 0.3%; the length of the Ganjiang River basin is 3693km, 3682km and 3684km respectively, and the maximum value relative to the minimum value change rate is 0.3%. The areas of the upper, middle and whole Ganjiang River basins are equal under the three methods, which are 39142km respectively. 2 、61819km 2 、79137km 2The river network densities of the upper reaches, the upper and middle reaches, and the entire basin of the Ganjiang River are equal, being 5.1 km / km², 61819 km / km², and 79137 km / km² respectively. The results show that the lengths of the river systems, the basin areas, and the river network densities extracted by the ArcGIS hydrological analysis tool (AG), the D8 classical algorithm of Matlab software (ML), and the D8 parallel algorithm supported by CUDA are close to each other.

[0064] (2) Parallel acceleration effect and scale effect of the D8 algorithm under the CUDA environment

[0065] Set the scenarios of the number of thread blocks and computational grids, and draw the response surface of the acceleration ratio of the D8 parallel algorithm under the CUDA environment relative to the D8 classical algorithms of ArcGIS and Matlab ( Figure 5 ) The results show that: in the upper reaches of the Ganjiang River, the acceleration ratio of the CUDA parallel operation varies between 6 and 53 relative to ArcGIS; and between 53 and 458 relative to Matlab. In the upper and middle reaches of the Ganjiang River, the acceleration ratio of the CUDA parallel operation varies between 5 and 41 relative to ArcGIS; and between 55 and 472 relative to Matlab. In the Ganjiang River basin, the acceleration ratio of the CUDA parallel operation varies between 5 and 45 relative to ArcGIS; and between 55 and 500 relative to Matlab. Comparing the time consumption of water system extraction in the upper reaches, the upper and middle reaches, and the entire basin of the Ganjiang River, the D8 algorithm under the CUDA environment < the ArcGIS hydrological analysis tool < the D8 classical algorithm of Matlab software in terms of time consumption; the acceleration ratio increases with the increase in the number of thread blocks and grids, and is more sensitive to the number of thread blocks.

[0066] Draw the marginal change curve of the acceleration ratio with the number of grids under the same number of thread blocks ( Figure 6) and calculated the optimal CUDA parallel speedup ratios for different thread blocks and grids (Table 1). In the upper Ganjiang River, relative to ArcGIS, the optimal speedup ratios for CUDA thread blocks of 4 to 128 and 256 to 1024 were 17 to 36 (grid number ≤ 1024) and 39 to 53 (grid number ≥ 65536), respectively. Relative to Matlab, the optimal speedup ratios for thread blocks of 4 to 128 and 256 to 1024 were 148 to 314 (grid number ≤ 1024) and 334 to 458 (grid number ≥ 65536), respectively. In the upper and middle reaches of the Ganjiang River, relative to ArcGIS, the optimal speedup ratios under the CUDA environment are 12-28 (grid number ≤ 1024) and 29-41 (grid number ≥ 65536) when the thread block size is 4-128 and 256-1024, respectively. Compared with Matlab, the optimal speedup ratios under the CUDA environment are 143-321 (grid number ≤ 1024) and 334-472 (grid number ≥ 65536) when the thread block size is 4-128 and 256-1024, respectively. In the Ganjiang River Basin, relative to ArcGIS, the optimal speedup ratios for CUDA thread block sizes of 4-128 and 256-1024 were 13-30 (grid number ≤ 1024) and 31-45 (grid number ≥ 65536), respectively. Relative to Matlab, the optimal speedup ratios for thread block sizes of 4-128 and 256-1024 were 147-328 (grid number ≤ 1024) and 343-500 (grid number ≥ 65536), respectively. The results show that in the upper, middle, and entire Ganjiang River basins, the optimal speedup ratio occurs when the number of threads is 128 or less, with a grid number below 1024; and when the number of threads is greater than 128, the optimal speedup ratio occurs when the number of grids is 65536 or greater.

[0067] Using the upper reaches of the Ganjiang River as a benchmark, the speedup ratios of the upper, middle, and entire Ganjiang River basins were compared under different threading conditions. The results showed that the speedup relative to ArcGIS for the upper and middle reaches of the Ganjiang River decreased by 23–30%, with an average decrease of 26%. The speedup relative to Matlab for the upper Ganjiang River varied by -6–3%, with an average of -0.3%. Within the Ganjiang River basin, the speedup relative to ArcGIS for the upper Ganjiang River decreased by 15–24%, with an average decrease of 20%. The speedup relative to Matlab for the upper Ganjiang River varied by -2–9%, with an average of 3%. The results showed that the speedup ratio of the D8 parallel algorithm relative to ArcGIS decreased significantly under the CUDA environment with increasing spatial scale, with the speedup ratio decreasing by over 20% in the upper, middle, and entire Ganjiang River basins. The speedup ratio of the D8 parallel algorithm relative to Matlab for the upper Ganjiang River basin did not change significantly under the CUDA environment.

[0068] Table 1 CUDA parallel optimal speedup ratio and grid of different thread blocks in Ganjiang River Basin

[0069]

[0070] The above are only preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A parallel analysis method for water flow direction and cumulative flow based on CUDA architecture D8 algorithm, characterized by: The following steps are involved: Construct a reading and storage model for watershed DEM raster data based on GPU texture memory; Based on the GPU global memory, set the parallel grid and thread block, map the thread to the DEM grid cell, and obtain the current thread identification number; Each thread continuously picks up the central grid and adjacent grid elevation values of the collected watershed DEM raster data according to the read storage model; The maximum distance weighted drop is calculated based on the elevation values of the central grid and adjacent grids, and the flow direction of the central grid is determined based on the maximum distance weighted drop; Construct flow direction grid cells according to the center grid flow direction, map threads to flow direction grid cells, and read thread identification numbers and flow directions; Each thread performs serial confluence calculations according to the flow direction. Conflicts in accessing and calculating the same grid are calculated through the CUDA atomic add function to obtain the cumulative confluence matrix. The construction of the read storage model includes the following steps: Declare a two-dimensional CUDA array and copy the water area DEM raster data from the GPU global memory to the array; Set the return data type to a two-dimensional floating-point global variable; Set the non-normalized texture coordinate reading mode and boundary and picking method properties to determine whether the texture coordinates are out of bounds; Bind the CUDA array to a texture reference, and fetch the elevation values of the center grid and its adjacent grids from the texture.

2. The parallel analysis method for water area flow direction and cumulative sink volume based on the CUDA architecture D8 algorithm according to claim 1 is characterized in that: Also includes: Preprocessing of water area DEM raster data, preprocessing of DEM raster data of filled flat land respectively.

3. The parallel analysis method for water area flow direction and cumulative sink volume based on the CUDA architecture D8 algorithm according to claim 1 is characterized in that: Determining the center grid flow direction according to the maximum distance weight drop specifically includes the following steps: Get the current thread identification number and map the thread to the DEM grid cell; If the current thread identification number is less than the total number of DEM grids, the parallel calculation starts; Calculate the two-dimensional coordinates of the central grid unit and continuously pick the elevation values of the central grid and adjacent grids; The maximum grid distance weight drop is calculated based on the elevation values of the central grid and adjacent grids to determine the flow direction of the central grid; If the current thread identification number is greater than the total number of DEM grids, the parallel calculation of this thread ends.

4. The parallel analysis method for water area flow direction and cumulative sink volume based on the CUDA architecture D8 algorithm according to claim 1 is characterized in that: Each thread continuously traverses and calculates through the CUDA atomic add function to obtain the cumulative sink flow, specifically including: Get the current thread identification number and map the thread to the flow grid unit; If the current thread identification number is less than the total number of traffic grids, parallel traffic calculation begins; Calculate the two-dimensional coordinates of the central grid unit and continuously pick the central grid flow direction; According to the flow direction of the central grid, the CUDA atomic add function is used to traverse and calculate the current grid cumulative flow; If the picked flow direction grid is out of bounds, the loop calculation ends; If the current thread identification number is greater than the total number of traffic grids, the parallel calculation of this thread ends.

5. The parallel analysis method for water area flow direction and cumulative sink volume based on the CUDA architecture D8 algorithm according to claim 1 is characterized in that: Also includes: According to the cumulative flow matrix, the length of rivers below the flow threshold is counted, and the river network density and the drainage area are calculated.

Images