River channel confluence stabilization and efficient calculation method of distributed hydrological model CASC2D

Abstract

The invention relates to a river channel confluence stabilization and efficient calculation method of a distributed hydrological model CASC2D, and the method comprises the steps: employing a semi-implicit local inertial wave method in a river channel unit to carry out the stable solving of a water flow evolution process, effectively inhibiting the numerical oscillation, and improving the model calculation efficiency; meanwhile, parallel computing is carried out on the double-layer nested grid structure of the model through OpenMP, a thread safety control and dynamic scheduling mechanism is matched, synchronous waiting and data competition are remarkably reduced, and the parallel utilization rate is increased. The method is particularly suitable for flood simulation of small and medium-sized watersheds in mountainous and hilly areas, the model calculation efficiency can be improved by 15 times while high-precision simulation is kept, and the risk of numerical instability under the short-time high-intensity rainfall condition is remarkably reduced. The method does not need to depend on a GPU or a distributed cluster platform, has high universality, good portability and engineering realizability, and is suitable for various application scenes such as flood forecasting, risk analysis and real-time scheduling.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of hydrological simulation and high-performance computing, and particularly relates to a river confluence stabilization and efficient calculation method of a distributed hydrological model CASC2D. BACKGROUND

[0002] The distributed hydrological model is widely used because it can better describe the spatial and temporal variation process of water cycle. The distributed hydrological model CASC2D using the mechanism of excess infiltration runoff generation takes a grid as a spatial discrete unit, and performs confluence calculation along the flow direction by grid, so that the spatial and temporal description of the hydrological process of a small river can be more accurate. The existing CASC2D model has the following defects: the river confluence is unstable: the one-dimensional explicit finite difference diffusion wave method is used for river confluence calculation, the numerical stability is poor, a shorter river confluence time step is needed to maintain the stability of the calculation, which increases the number of time cycles and the running time cost; the serial calculation mode is difficult to support the real-time flood forecasting demand of high resolution: the application of the CASC2D model at the present stage is mainly limited by the calculation efficiency, and the calculation efficiency mainly depends on the precision interval of the grid unit, the smaller the scale division, the more accurate the geographical features represented, and the more accurate the simulation results, so that the calculation amount and the calculation time of the CASC2D model increase in geometric multiples, resulting in very low calculation efficiency and stability of the simulation. SUMMARY

[0003] In view of the defects of the prior art, the present application provides a river confluence stabilization and efficient calculation method of a distributed hydrological model CASC2D, which can improve the stability and efficiency of the model calculation.

[0004] To achieve the above object, the technical scheme adopted by the present application is as follows: A river confluence stabilization and efficient calculation method of a distributed hydrological model CASC2D, comprising the following steps: Step S1, constructing a watershed grid structure based on a digital elevation model DEM, soil type and land use type, using Arcgis software to extract slope and river unit information, and generating a model input file; Step S2, using a semi-implicit local inertial wave equation to solve the river flow, wherein the semi-implicit local inertial wave equation is calculated by using a semi-implicit difference format of water depth and flow as the main variable; Step S3, using OpenMP to realize multi-thread parallel processing of the model calculation process, and optimizing the traversal and processing of rainfall interpolation, slope runoff, slope confluence and river confluence; Step S4, estimating and rating the model parameters, and outputting the spatial and temporal distribution runoff results in the simulation process.

[0005] Further, in the step S1, the basin grid structure is constructed based on a digital elevation model DEM, a soil type and a land use type, slope surface and river channel unit information is extracted by using Arcgis software, and a model input file is generated, specifically as follows: S1.1, depression filling processing, in the ArcGIS software, opening the spatial analysis tool in the toolbox, selecting the depression filling function in the hydrological analysis module; taking the basin digital elevation model DEM data as input data, performing depression filling processing, generating the elevation data file Fill.tif after depression filling; converting the Fill.tif file into ASCII format to obtain the Fill.asc file; S1.2, flow direction analysis, in the spatial analysis tool of the ArcGIS software, selecting the flow direction function in the hydrological analysis module; taking the Fill.tif file obtained in the step S1.1 as input data, performing flow direction analysis by using the D8 algorithm, generating the flow direction data file Dir.tif; converting the Dir.tif file into ASCII format to obtain the Dir.asc file; S1.3, flow accumulation calculation, in the spatial analysis tool of the ArcGIS software, selecting the flow accumulation function in the hydrological analysis module; taking the Dir.tif file obtained in the step S1.2 as input data, performing flow accumulation calculation, generating the flow accumulation data file Acc.tif; converting the Acc.tif file into ASCII format to obtain the Acc.asc file; S1.4, river channel extraction, in the spatial analysis tool of the ArcGIS software, selecting the raster calculator function in the map algebra module; taking the Acc.tif file obtained in the step S1.3 as input data, setting a predetermined flow accumulation threshold T, extracting the river channel network by raster calculation, generating the river channel data file Str.tif; converting the Str.tif file into ASCII format to obtain the Str.asc file; S1.5, river link generation, in the spatial analysis tool of the ArcGIS software, selecting the river link function in the hydrological analysis module; taking the Str.tif file obtained in the step S1.4 and the Dir.tif file obtained in the step S1.2 as input data, generating the river link data file Link.tif; converting the Link.tif file into ASCII format to obtain the Link.asc file; S1.6, node coding, based on the Acc.asc file obtained in the step S1.3, for the nodes with the same Link code, the flow accumulation is sorted according to the size, the Node code is generated for each river link, and the node coding data file Node.asc is output; S1.7, river network classification, in the spatial analysis tool of ArcGIS software, select the river network classification function of hydrological analysis module; take the Str.tif file obtained in step S1.4 and the Dir.tif file obtained in step S1.2 as input data, select the Strahler river classification method to classify the river, and generate a classified data file riverorder.tif; Step S1.8, according to the soil tif file, using the spatial analysis tool in the toolbox, execute the mask extraction command to extract the soil type of the basin, and then convert the file to soil.asc file; according to the land use tif file, using the spatial analysis tool in the toolbox, execute the mask extraction command to extract the land use type of the basin, and then convert the file to landuse.asc file; Step S1.9, the river information file required by the input file includes the river width, river depth and river roughness information of each river node; the river depth is set to the maximum value, and the river width of each river grid is calculated.

[0006] Specifically, in step S1.9, the river roughness is the Manning roughness coefficient of the river, which is determined in combination with the Strahler river classification obtained in step S1.7, and the Strahler river classification includes I, II, III, IV, V, VI and 0, a total of six levels, each level corresponds to a respective roughness.

[0007] Specifically, in step S1.9, the river width of each river grid is calculated according to the following formula: ; In the formula, represents the maximum river width, m; represents the minimum river width of the critical area, m; represents the critical area, m 2 ; that is, the minimum catchment area required to form a river, m 2 ; represents the upstream catchment area of the river grid numbered, m 2 ; represents the total area of the whole basin, m 2 ; in Arcgis software, the upstream catchment area of the grid numbered is obtained by multiplying the flow accumulation value of each river grid by the area of each grid , m 2 .

[0008] Further, in step S2, the semi-implicit local inertia wave equation is calculated by using a semi-implicit difference format with water depth and flow As the primary calculation variable, the calculation expression for the semi-implicit difference scheme is: In the formula: express time The calculated flow rate at point m 3 / s; express time The calculated flow rate at point m 3 / s; express time The water depth at point i is in meters; represents the water surface gradient at time t+1 / 2. It is the acceleration caused by gravity, m 3 / s; S The water surface gradient; The water depth at the center of the river segment, in meters; subscript. They are Pixel spatial index in the direction; It is time; It is the time step, in seconds; It is the spacetime step, in meters.

[0009] Furthermore, in step S3, OpenMP is used to implement multi-threaded parallel processing of the model calculation process, optimizing the traversal and processing of rainfall interpolation, slope runoff generation, slope confluence, and river confluence. Specifically, this includes the following steps: Step S3.1, Rainfall Data Processing Stage: Before the parallel region begins, the rain station data reading operation is completed, and the rainfall intensity data is stored in the global array rrg, i.e., the rain station data array; when the rainfall calculation flag icall is equal to 1, the rainfall calculation function rain is called; the rainfall calculation rain function performs independent calculations for each grid cell, using the rain station horizontal coordinate array xrg and the rain station vertical coordinate array yrg to calculate the grid position, and performing a distance-weighted average based on the rain station data array rrg, and storing the calculation result in the grid rainfall array rainfall; Step S3.2, Slope runoff calculation stage: The update of the grid water depth array h adopts grid-level parallel computing, wherein: each grid cell independently calculates the single grid runoff dqq, and the unit is converted by the grid width parameter w and the time step parameter dt; the update of the minimum rainfall statistics aminrain and the maximum rainfall statistics amaxrain adopts reduction operation, and the reduction operation clause of OpenMP is used to ensure thread safety, and the reduction type is set to minimum value min and maximum value max; Step S3.3, slope confluence calculation stage, in the slope confluence function OvrlRout: parallelize the outer loop j with OpenMP; set the scheduling strategy schedule as dynamic scheduling dynamic, in the single grid confluence calculation function ovrl, the update of the inter-grid runoff flux array dqov is protected by atomic operation: #pragma omp atomic dqov[j][k]-= dqq; / / j,k for the current grid coordinates #pragma omp atomic dqov[jj][kk]+= dqq; / / jj,kk for adjacent grid coordinates Step S3.4, river confluence calculation stage, in the river confluence main function ChannRout, parallelize the river link number loop ic, each river grid belongs to only a single link, and ensure no data race; in the single river calculation function chnchn, the update of the river runoff flux array dqch is protected by atomic operation: #pragma omp atomic dqch[j][k]-= dq; / / j,k for the current river grid coordinates #pragma omp atomic dqch[jj][kk]+= dq; / / jj,kk for downstream river grid coordinates.

[0010] Further, in step S4, the model parameters are estimated and calibrated, and the spatial and temporal distribution of runoff results in the simulation process is output, and the specific steps are as follows: Step S4.1, estimation of vegetation interception parameters, based on the land use type corresponding to the grid, the vegetation interception depth is set, the land use type includes forest, water surface, farmland, grassland; different land use types correspond to different vegetation interception depths; the land use type data is derived from the landuse.asc file, and the vegetation interception depth is used to describe the initial rainfall loss; Step S4.2, estimation of infiltration parameters, the Green-Ampt model is used to simulate the infiltration process, which needs to input the following three types of physical parameters: saturated hydraulic conductivity, capillary water, and soil initial water deficit, the parameters are assigned according to the soil type, the soil type refers to the Rawls experimental data, including Calloway, Fallaya, Grenada, Loring, Collins, Memphis, Gullied Land seven types, which correspond to different physical parameter values and soil particle size composition, the initial water deficit of all grids in the watershed is uniformly set to an empirical value; Step S4.3, Slope runoff parameter estimation: The Manning roughness coefficient of the land surface is determined based on the land use classification in the landuse.asc file. The specific correspondence is as follows: The Manning roughness coefficient corresponding to forest land is... The Manning roughness coefficient corresponding to the water surface is: The Manning roughness coefficient corresponding to cultivated land is: The Manning roughness coefficient corresponding to grassland is: ; All parameters can be assigned values ​​in ArcGIS software through raster reclassification and spatial linking, and then used in model calibration. The calibration process can combine observed runoff processes or peak flow rates, and use automated algorithms or manual empirical adjustment methods to optimize key parameters.

[0011] Compared with the prior art, the technical solution of the present invention has the following beneficial technical effects: significantly improves the simulation accuracy of rainstorm floods in small watersheds in mountainous areas, especially under short-term and local heavy rainfall conditions; effectively avoids problems such as negative water depth and flow oscillation in traditional methods; supports larger time step calculations and improves overall calculation efficiency; the model is robust and suitable for terrain with elevation differences and complex river network structures. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the process of the present invention; Figure 2 This is a schematic diagram of the DEM and rain gauge stations in the Jiangwan watershed of this invention; Figure 3 These are the flood process curves simulated by the original model and the improved model. Detailed Implementation

[0013] The present invention will be further described below with reference to specific embodiments: like Figure 1 As shown, a method for stabilizing and efficiently calculating river confluence in the distributed hydrological model CASC2D includes the following steps: Step S1: Construct a watershed grid structure based on the Digital Elevation Model (DEM), soil type, and land use type. Use ArcGIS software to extract slope and channel unit information and generate the model input file. This includes the following steps: Step S1.1, Depression Filling Process: In ArcGIS software, open the Spatial Analysis tool in the toolbox, select the Depression Filling function in the Hydrological Analysis module; use the watershed digital elevation model (DEM) data as input data, perform depression filling processing, and generate the depression-filled elevation data file Fill.tif; convert the Fill.tif file to ASCII format to obtain the Fill.asc file; Step S1.2, flow direction analysis, in the spatial analysis tool of ArcGIS software, select the flow direction function of the hydrological analysis module; take the Fill.tif file obtained in step S1.1 as input data, and use the D8 algorithm to perform flow direction analysis to generate a flow direction data file Dir.tif; convert the Dir.tif file into ASCII format to obtain a Dir.asc file; Step S1.3, accumulation calculation, in the spatial analysis tool of ArcGIS software, select the accumulation function of the hydrological analysis module; take the Dir.tif file obtained in step S1.2 as input data to perform accumulation calculation to generate an accumulation data file Acc.tif; convert the Acc.tif file into ASCII format to obtain an Acc.asc file; Step S1.4, river extraction, in the spatial analysis tool of ArcGIS software, select the raster calculator function of the map algebra module; take the Acc.tif file obtained in step S1.3 as input data, set a predetermined accumulation threshold T, and extract the river network by raster calculation to generate a river data file Str.tif; convert the Str.tif file into ASCII format to obtain a Str.asc file; Step S1.5, river link generation, in the spatial analysis tool of ArcGIS software, select the river link function of the hydrological analysis module; take the Str.tif file obtained in step S1.4 and the Dir.tif file obtained in step S1.2 as input data to generate a river link data file Link.tif; convert the Link.tif file into ASCII format to obtain a Link.asc file; Step S1.6, node coding, based on the Acc.asc file obtained in step S1.3, sort the nodes with the same Link code according to the accumulation size, generate Node code for each river chain, and output a node coding data file Node.asc; Step S1.7, river network grading, in the spatial analysis tool of ArcGIS software, select the river network grading function of the hydrological analysis module; take the Str.tif file obtained in step S1.4 and the Dir.tif file obtained in step S1.2 as input data, select the Strahler river grading method to grade the river, and generate a grading data file riverorder.tif; Step S1.8, according to the soil tif file, using the spatial analysis tool in the toolbox, the command of extraction by mask is executed to extract the soil type of the watershed, and then the file is converted into a soil.asc file; according to the land use tif file, using the spatial analysis tool in the toolbox, the command of extraction by mask is executed to extract the land use type of the watershed, and then the file is converted into a landuse.asc file; Step S1.9, the river channel information file required by the input file mainly includes the river width, river depth and river roughness information of each river node. In mountainous areas, the floodplain is not considered for the time being, and the maximum value of 1000m is set for the river depth; the river width of each river grid is calculated according to formula one: ; In the formula, represents the maximum river width, m; represents the minimum river width of the critical area, m; represents the critical area, m 2 ; that is, the minimum catchment area required to form a river channel, m 2 ; represents the upstream catchment area of the river grid numbered, m 2 ; represents the total area of the entire watershed, m 2 ; in Arcgis software, the upstream catchment area of the grid numbered is obtained by multiplying the convergence accumulation value of each river grid by the area of each grid , m 2 .

[0014] River roughness is the Manning roughness coefficient of the river channel, which is determined in combination with the Strahler river classification obtained in step 1.7, as shown in Table 1 for the river roughness coefficient; Table 1 ;

[0015] Step S2, the semi-implicit local inertia wave equation is used to solve the river flow, which specifically includes: the semi-implicit difference format is used to solve the local inertia wave equation, with water depth and flow as the main calculation variables, and the calculation expression of the semi-implicit difference format is: ; In the formula, represents the calculation flow at time , m 3 / s; represents the calculation flow at time , m 3 / s; denotes time water depth at time t, m; denotes water surface slope at time t+1 / 2; is the acceleration due to gravity, m 3 / s; S is the water surface slope; is the water depth at the center of the reach element, m; the subscript are pixel space indices in the direction; is time; is the time step, s; is the space-time step, m.

[0016] Step S3, using OpenMP to realize multi-thread parallel processing of the main calculation process of the model, optimizing the traversal and processing of the calculation modules such as rainfall interpolation, slope runoff, slope confluence, and river confluence, specifically including the following steps: Step S3.1, rainfall data processing stage, complete rainfall station data reading operation before starting the parallel region, store the rainfall intensity data to the global array rrg (rainfall station data array); when the rainfall calculation flag icall (rainfall call flag) is equal to 1, call the rain (rainfall calculation function); the rain function internally performs independent calculation for each grid element, wherein the grid position is calculated using xrg (rainfall station horizontal coordinate array) and yrg (rainfall station vertical coordinate array), distance weighted average is performed based on rrg (rainfall station data array), and the calculation result is stored in rainfall (grid rainfall array); Step S3.2, slope runoff calculation stage, the update of h (grid water depth array) adopts grid-level parallel calculation, wherein: each grid element independently calculates dqq (single grid runoff), and the unit is converted through w (grid width parameter) and dt (time step parameter); the update of aminrain (minimum rainfall statistical value) and amaxrain (maximum rainfall statistical value) adopts reduction operation, and the reduction type is set to min (minimum value) and max (maximum value) using the reduction (reduction operation) clause of OpenMP to ensure thread safety.

[0017] Step S3.3, slope confluence calculation stage, in OvrlRout (slope confluence function): parallelize the outer loop j (grid row number loop) with OpenMP; set the schedule (scheduling strategy) to dynamic (dynamic scheduling), and in ovrl (single grid confluence calculation function), the update of dqov (inter-grid runoff flux array) adopts atomic (atomic operation) protection: #pragma omp atomic dqov[j][k] -= dqq; / / j,k for the current grid coordinates #pragma omp atomic dqov[jj][kk] += dqq; / / jj,kk for the adjacent grid coordinates Step S3.4, river confluence calculation stage, in ChannRout (river confluence main function), ic (river link number loop) is parallelized, each river grid belongs to only a single link, ensuring no data competition; in chnchn (single river calculation function), the update of dqch (river runoff flux array) adopts atomic operation: #pragma omp atomic dqch[j][k] -= dq; / / j,k for the current river grid coordinates #pragma omp atomic dqch[jj][kk] += dq; / / jj,kk for the downstream river grid coordinates Step S4, estimation and calibration of model parameters, output the spatial and temporal distribution of runoff results in the simulation process, the specific steps are as follows: Step S4.1, estimation of vegetation interception parameters According to the land cover type (landcover.asc) corresponding to the grid, the vegetation interception depth I is set, which is used to describe the initial rainfall loss. Land cover types are divided into forest, water surface, farmland and grassland, which correspond to different measured interception values, as shown in the following Table 2 for plant interception depth: Table 2 ;

[0018] Estimation of infiltration parameters The Green-Ampt model is used to simulate the infiltration process. The model needs to input the following three types of physical parameters: saturated hydraulic conductivity, capillary water, and initial soil water deficit. The parameters are assigned according to the soil type (soil.asc). The soil type is based on the experimental data of Rawls, including Calloway, Fallaya, Grenada, Loring, etc. Seven types correspond to different physical parameter values and soil particle size composition. The initial water deficit of all grids in the basin can be uniformly set to an empirical value (e.g. 0.2); as shown in Table 3 for soil experimental values and composition percentage: Table 3 ;

[0019] Step S4.3, estimation of slope confluence parameters The Manning roughness coefficient of the ground surface is determined in combination with land use according to the value proposed by Woolhiser in 1975, as shown in Table 4, the ground roughness coefficient; Table 4 ;

[0020] All parameters can be assigned by means of grid reclassification and spatial connection in ArcGIS, and then participate in model calibration; the calibration process can combine the observed runoff process or flood peak flow, use automatic algorithm (such as genetic algorithm, SCE-UA) or manual experience adjustment method, and optimize the key parameters.

[0021] Embodiment 1 Taking the Jiangwan basin in Zhejiang as an example, the basin is located in the Dongqianxi basin, which contains 10 rainfall stations. The DEM and rainfall station distribution of Jiangwan are shown in Figure 2 The annual average precipitation is 1580mm, and the Meiyu in May and June, typhoon rain in August and September, and thunderstorm in summer are the main precipitation in the rainy season. There is a small amount of snow in winter. The precipitation distribution is still uniform, but the rainfall decreases with the height due to the influence of topography. The digital elevation data (DEM) in the study area of this embodiment uses the 90m resolution STRIM (Shuttle Radar Topography Mission) data provided by the United States Space Agency (NASA) and the National Mapping Bureau (NIMA) of the Ministry of Defense.

[0022] Step S1, based on the digital elevation model DEM, soil type and land use type, the grid structure of the basin is constructed, the slope and river channel unit information is extracted by using Arcgis software, and the model input file is generated, as follows: Step S1.1, fill-in processing, in the ArcGIS software, open the spatial analysis tool in the toolbox, select the fill-in function in the hydrological analysis module; take the digital elevation model DEM data of the basin as the input data, execute the fill-in processing, and generate the fill-in elevation data file Fill.tif; convert Fill.tif file to ASCII format to obtain Fill.asc file; Step S1.2, flow direction analysis, in the spatial analysis tool of ArcGIS software, select the flow direction function of the hydrological analysis module; take the Fill.tif file obtained in step S1.1 as the input data, and use D8 algorithm to perform flow direction analysis, to generate flow direction data file Dir.tif; convert Dir.tif file to ASCII format to obtain Dir.asc file; Step S1.3, accumulation calculation. In the spatial analysis tool of ArcGIS software, the accumulation function of the hydrological analysis module is selected. The Dir.tif file obtained in step S1.2 is taken as input data, and the accumulation calculation is performed to generate an accumulation data file Acc.tif. The Acc.tif file is converted into ASCII format to obtain an Acc.asc file. Step S1.4, river extraction. In the spatial analysis tool of ArcGIS software, the raster calculator function of the map algebra module is selected. The Acc.tif file obtained in step S1.3 is taken as input data, and a predetermined accumulation threshold T is set. The river network is extracted by raster calculation to generate a river data file Str.tif. The Str.tif file is converted into ASCII format to obtain a Str.asc file. Step S1.5, river link generation. In the spatial analysis tool of ArcGIS software, the river link function of the hydrological analysis module is selected. The Str.tif file obtained in step S1.4 and the Dir.tif file obtained in step S1.2 are taken as input data to generate a river link data file Link.tif. The Link.tif file is converted into ASCII format to obtain a Link.asc file. Step S1.6, node coding. Based on the Acc.asc file obtained in step S1.3, the nodes with the same Link code are sorted according to the accumulation size, and the Node code is generated for each river chain. The node coding data file Node.asc is output. Step S1.7, river network grading. In the spatial analysis tool of ArcGIS software, the river network grading function of the hydrological analysis module is selected. The Str.tif file obtained in step S1.4 and the Dir.tif file obtained in step S1.2 are taken as input data, and the Strahler river grading method is selected for river grading to generate a grading data file riverorder.tif. Step S1.8, according to the soil.tif file, the spatial analysis tool in the toolbox is used to execute the mask extraction command to extract the soil type of the basin, and then the file is converted into a soil.asc file. According to the land use.tif file, the spatial analysis tool in the toolbox is used to execute the mask extraction command to extract the land use type of the basin, and then the file is converted into a landuse.asc file. Step S1.9, the river information file required by the input file mainly includes the river width, river depth, and river roughness information of each river node. In mountainous areas, the floodplain is not considered, and the river depth is set to a maximum value of 1000m. The river width of each river grid is calculated according to formula one: ; wherein, Represents the maximum width of the river channel, in meters (m). The minimum channel width representing the critical area, in meters; Represents the critical area, m 2 That is, the minimum catchment area required to form a river channel, in meters. 2 ; The upstream catchment area of ​​the river grid with the number is represented in m. 2 ; Represents the total area of ​​the entire watershed, m 2 In ArcGIS software, the upstream catchment area of ​​the grid numbered cell is calculated by multiplying the cumulative runoff value of each river grid by the area of ​​each grid cell. m 2 .

[0023] Step S2: Solve the river flow equation using a semi-implicit local inertial wave equation, specifically including: solving the local inertial wave equation using a semi-implicit difference scheme, with water depth as the key factor. and traffic As the primary calculation variable, the calculation expression for the semi-implicit difference scheme is: In the formula: express time The calculated flow rate at point m 3 / s; express time The calculated flow rate at point m 3 / s; express time The water depth at point i is in meters; represents the water surface gradient at time t+1 / 2. It is the acceleration caused by gravity, m 3 / s; S The water surface gradient; The water depth at the center of the river segment, in meters; subscript. They are Pixel spatial index in the direction; It is time; It is the time step, in seconds; It is the spacetime step, in meters.

[0024] Step S3: Utilize OpenMP to implement multi-threaded parallel processing of the main calculation processes of the model, optimizing the traversal and processing of calculation modules such as rainfall interpolation, slope runoff generation, slope confluence, and river confluence. Specifically, this includes the following steps: Step S3.1, rainfall data processing stage, before parallel region starts, complete rainfall station data reading operation, store rainfall intensity data to global array rrg (rainfall station data array); when rainfall calculation flag icall (rainfall call flag) is equal to 1, call rain (rainfall calculation function); the rain function internally calculates each grid cell independently, wherein the grid position is calculated using xrg (rainfall station horizontal coordinate array) and yrg (rainfall station vertical coordinate array), distance weighted average is performed based on rrg (rainfall station data array), and the calculation result is stored in rainfall (grid rainfall array); Step S3.2, slope runoff calculation stage, update of h (grid water depth array) uses grid-level parallel calculation, wherein: each grid cell independently calculates dqq (single grid runoff), and units are converted through w (grid width parameter) and dt (time step parameter); the update of amainrain (minimum rainfall statistical value) and amaxrain (maximum rainfall statistical value) uses a reduction operation, and the reduction (reduction operation) clause of OpenMP is used to ensure thread safety, and the reduction type is set to min (minimum value) and max (maximum value).

[0025] Step S3.3, slope confluence calculation stage, in OvrlRout (slope confluence function): parallelize the outer loop j (grid row number loop) with OpenMP; set schedule (scheduling strategy) to dynamic (dynamic scheduling), and in ovrl (single grid confluence calculation function), the update of dqov (inter-grid runoff flux array) uses atomic (atomic operation) protection: #pragma omp atomic dqov[j][k]-= dqq; / / j,k for current grid coordinates #pragma omp atomic dqov[jj][kk]+= dqq; / / jj,kk for adjacent grid coordinates Step S3.4, river confluence calculation stage, in ChannRout (river confluence main function), parallelize ic (river link number loop), each river grid only belongs to a single link, and ensure no data competition; in chnchn (single river calculation function), the update of dqch (river runoff flux array) uses atomic operation: #pragma omp atomic dqch[j][k]-= dq; / / j,k for current river grid coordinates #pragma omp atomic dqch[jj][kk]+= dq; / / jj,kk for the downstream river grid coordinates Step S4, estimate and calibrate the model parameters, and output the spatial and temporal distribution of runoff results in the simulation process. The parameter estimation and calibration results are shown in Table 5 (land use parameter calibration results), Table 6 (soil type parameter calibration results); Table 5 ; Table 6 ; Table 7 ; Note: The CPU used is 11th Gen Intel(R) Core(TM) i9-11950H.

[0026] By comparing the simulation results of the original CASC2D model and the improved model with the introduction of local inertial wave method combined with OpenMP parallel optimization in two typical flood processes (1983100400 / 121h and 1985073008 / 168h), the improvement effect of the model in calculation efficiency and accuracy is evaluated. As shown in Table 7 (statistical table of acceleration results of CASC2D model and improved model), the maximum calculation time step of the original model is 2s, while the improved model can stably run at 18s step, showing higher numerical stability. In terms of simulation accuracy, the determination coefficients (Nash efficiency coefficients) of both models are above 0.92, and the determination coefficients of the improved model and the original model in the second flood are both 0.97, which shows that it does not cause significant error increase while ensuring simulation accuracy. Through Figure 3 The flood hydrograph simulated by the improved model and the original model can also show that the simulation accuracy of the improved model is similar to that of the original model. In terms of calculation efficiency, the calculation time of the improved model in the two flood processes is shortened to 2.1 minutes and 2.9 minutes, respectively, compared with 34.7 minutes and 44.83 minutes of the original model, and the calculation acceleration ratio is 16.5 times and 15.5 times, respectively, significantly improving the simulation efficiency. In summary, the improved model has strong stability and practicality.

[0027] In addition, it should be noted that each specific technical feature described in the above specific embodiments can be combined in any appropriate manner without contradiction. In order to avoid unnecessary repetition, the present application will not describe various possible combinations again.

[0028] The above shows and describes the basic principles, main features and advantages of the present application. Those skilled in the art should understand that the present application is not limited to the above-mentioned embodiments, and the above-mentioned embodiments and descriptions in the specification are only to illustrate the principles of the present application. Without departing from the spirit and scope of the present application, various changes and improvements can be made to the present application, and these changes and improvements all fall within the scope of the claimed present application. The scope of protection of the present application is defined by the appended claims and their equivalents.

Claims

1. A river channel confluence stabilization and efficient calculation method of a distributed hydrological model CASC2D, characterized in that, Comprise the following steps: Step S1, based on digital elevation model DEM, soil type, land use type constructs the basin grid structure, adopts Arcgis software to extract the slope surface and river channel unit information, generates the model input file; Step S2, the semi-implicit local inertia wave equation is used to solve the river flow, and the semi-implicit local inertia wave equation is calculated by using the semi-implicit difference format of water depth and flow as the main variable; Step S3, the OpenMP is used to realize the multi-thread parallel processing of the model calculation process, and the traversal and processing of rainfall interpolation, slope runoff, slope confluence and river confluence are optimized; Step S4, the model parameters are estimated and rated, and the time and space distribution runoff results in the simulation process are output.

2. The river channel confluence stabilization and efficient calculation method of the distributed hydrological model CASC2D according to claim 1, characterized in that, In the step S1, based on digital elevation model DEM, soil type, land use type constructs the basin grid structure, adopts Arcgis software to extract the slope surface and river channel unit information, generates the model input file, as follows: S1.1, filling treatment, in ArcGIS software, open the spatial analysis tool in the toolbox, select the filling function in the hydrological analysis module;With the digital elevation model DEM data of the basin as input data, execute the filling treatment, generate the elevation data file Fill.tif;Fill.tif file is converted into ASCII format, and Fill.asc file is obtained; S1.2, flow direction analysis, in the spatial analysis tool of ArcGIS software, select the flow direction function of the hydrological analysis module;With the Fill.tif file obtained in step S1.1 as input data, the flow direction analysis is carried out by using D8 algorithm, and the flow direction data file Dir.tif is generated;Dir.tif file is converted into ASCII format, and Dir.asc file is obtained; S1.3, flow accumulation calculation, in the spatial analysis tool of ArcGIS software, select the flow accumulation function of the hydrological analysis module;With the Dir.tif file obtained in step S1.2 as input data, the flow accumulation calculation is carried out, and the flow accumulation data file Acc.tif is generated;Acc.tif file is converted into ASCII format, and Acc.asc file is obtained; S1.4, river extraction, in the spatial analysis tool of ArcGIS software, select the raster calculator function of the map algebra module;With the Acc.tif file obtained in step S1.3 as input data, set the predetermined flow accumulation threshold T, extract the river network by raster calculation, generate the river data file Str.tif;Str.tif file is converted into ASCII format, and Str.asc file is obtained; S1.5, river link generation, in the spatial analysis tool of ArcGIS software, select the river link function of the hydrological analysis module; take the Str.tif file obtained in step S1.4 and the Dir.tif file obtained in step S1.2 as input data, and generate a river link data file Link.tif; convert the Link.tif file into ASCII format to obtain a Link.asc file; S1.6, node coding, based on the Acc.asc file obtained in step S1.3, sort the nodes with the same Link code according to the size of the confluence accumulation, generate Node code for each river link, and output a node coding data file Node.asc; S1.7, river network grading, in the spatial analysis tool of ArcGIS software, select the river network grading function of the hydrological analysis module; take the Str.tif file obtained in step S1.4 and the Dir.tif file obtained in step S1.2 as input data, select the Strahler river grading method to grade the river, and generate a grading data file riverorder.tif; Step S1.8, according to the soil.tif file, use the spatial analysis tool in the toolbox to execute the mask extraction command to extract the soil type of the basin, and then convert the file to soil.asc file; According to the land use.tif file, use the spatial analysis tool in the toolbox to execute the mask extraction command to extract the land use type of the basin, and then convert the file to landuse.asc file; Step S1.9, the input file required river information file includes the river width, river depth and river roughness information of each river node; set the river depth to the maximum value, and calculate the river width of each river grid.

3. The method for river channel confluence stabilization and efficient calculation of the distributed hydrological model CASC2D according to claim 2, characterized in that, In step S1.9, the river roughness is the Manning roughness coefficient of the river, which is determined in combination with the Strahler river grading obtained in step S1.7, and the Strahler river grading includes I, II, III, IV, V, VI and 0, a total of six levels, each level corresponds to its own roughness.

4. The river channel confluence stabilization and efficient calculation method of the distributed hydrological model CASC2D according to claim 2, characterized in that, In step S1.9, the river width of each river grid is calculated according to the following formula: ; where, represents the maximum river width, m; represents the minimum river width of the critical area, m; represents the critical area, m 2 ; that is, the minimum catchment area required to form a river, m 2 ; represents the upstream catchment area of the river grid numbered, m 2 ; represents the total area of the entire watershed, m 2 ; in Arcgis software, the upstream catchment area of the grid numbered is calculated by multiplying the convergence accumulation value of each river grid by the area of each grid , m 2 .

5. The method for river confluence runoff stabilization and efficient calculation of the distributed hydrological model CASC2D according to claim 1, characterized in that, In step S2, the semi-implicit local inertia wave equation is calculated by using a semi-implicit difference scheme with water depth and flow rate as main variables. Specifically, the semi-implicit difference scheme is calculated by using water depth and flow rate as main variables, and the calculation expression of the semi-implicit difference scheme is as follows: In the formula: express time The calculated flow rate at point m 3 / s; express time The calculated flow rate at point m 3 / s; express time The water depth at point i is in meters; represents the water surface gradient at time t+1 / 2. It is the acceleration caused by gravity, m 3 / s; S The water surface gradient; The water depth at the center of the river segment, in meters; subscript. They are Pixel spatial index in the direction; It is time; It is the time step, in seconds; It is the spacetime step, in meters.

6. The method for river confluence runoff stabilization and efficient calculation of the distributed hydrological model CASC2D according to claim 1, characterized in that, In step S3, OpenMP is used to realize multi-thread parallel processing of the model calculation process, and the traversal and processing of rainfall interpolation, slope runoff, slope confluence and river confluence are optimized, which includes the following steps: Step S3.1, rainfall data processing stage, complete rainfall station data reading operation before starting parallel region, store rainfall intensity data to global array rrg, that is, rainfall station data array; when the rainfall calculation flag icall is equal to 1, call the rainfall calculation function rain; the rainfall calculation rain function calculates each grid unit independently, wherein the grid position is calculated using the rainfall station horizontal coordinate array xrg and the rainfall station vertical coordinate array yrg, the distance weighted average is calculated based on the rainfall station data array rrg, and the calculation result is stored in the grid rainfall array rainfall; Step S3.2, the slope runoff calculation stage, the update of the grid water depth array h adopts grid level parallel calculation, wherein: each grid unit independently calculates single grid runoff dqq, which is converted by grid width parameter w and time step parameter dt; the update of the minimum rainfall statistical value aminrain and the maximum rainfall statistical value amaxrain adopts reduction operation, the OpenMP reduction operation clause is used to ensure thread safety, and the reduction type is set to minimum value min and maximum value max; Step S3.3, the slope confluence calculation stage, in the slope confluence function OvrlRout: the outer loop j is parallelized by OpenMP; the scheduling strategy schedule is set to dynamic scheduling dynamic, and in the single grid confluence calculation function ovrl, the update of the inter-grid runoff flux array dqov adopts atomic operation atomic protection: #pragma omp atomic dqov[j][k] -= dqq; / / j,k for the current grid coordinates #pragma omp atomic dqov[jj][kk] += dqq; / / jj,kk for adjacent grid coordinates Step S3.4, the river confluence calculation stage, in the river confluence main function ChannRout, the river link number loop ic is parallelized, each river grid belongs to only a single link, and no data competition is ensured; in the single river calculation function chnchn, the update of the river runoff flux array dqch adopts atomic operation: #pragma omp atomic dqch[j][k] -= dq; / / j,k for the current river grid coordinates #pragma omp atomic dqch[jj][kk] += dq; / / jj,kk for the downstream river grid coordinates.

7. The method for river confluence runoff stabilization and efficient calculation of the distributed hydrological model CASC2D according to claim 1, characterized in that, In step S4, the model parameters are estimated and calibrated, and the spatial and temporal distribution of runoff results in the simulation process is output, and the specific steps are as follows: Step S4.1, vegetation interception parameter estimation, the vegetation interception depth is set based on the land use type corresponding to the grid, and the land use type includes forest, water surface, farmland and grassland; different land use types correspond to different vegetation interception depths; the land use type data is derived from the landuse.asc file, and the vegetation interception depth is used to describe the initial rainfall loss; Step S4.2, the estimation of the infiltration parameter, the Green-Ampt model is used to simulate the infiltration process, which needs to input the following three types of physical parameters: saturated hydraulic conductivity, capillary water, and soil initial water deficit, the parameters are assigned according to the soil type, the soil type refers to the Rawls experimental data, including Calloway, Fallaya, Grenada, Loring, Collins, Memphis and GulliedLand seven types, which correspond to different physical parameter values and soil particle size composition, and the initial water deficit of all grids in the basin is uniformly set to an empirical value; Step S4.3, slope confluence parameter estimation, the surface Manning roughness coefficient is determined according to the land use classification in the landuse.asc file, and the specific corresponding relationship is as follows: the surface Manning roughness coefficient corresponding to forest land is ; the surface Manning roughness coefficient corresponding to water surface is ; the surface Manning roughness coefficient corresponding to farmland is ; and the surface Manning roughness coefficient corresponding to grassland is ; All parameters can be assigned by reclassifying the raster and spatial joining in ArcGIS software, and then participate in the model calibration. The calibration process can combine the observed runoff process or peak flow, using an automated algorithm or manual experience adjustment method to optimize the key parameters.

Images