Quick calculation method for storage capacity curve

Abstract

The invention relates to the field of emergency dispatching calculation of a barrier lake, and discloses a rapid calculation method for a reservoir capacity curve of a barrier lake, which comprises the following steps of: acquiring watercourse satellite DEM or DSM topographic data of a barrier lake reservoir area in a dry season, and performing multi-source heterogeneous data fusion to obtain topographic data of a submerged area of the barrier lake reservoir area; acquiring water level elevation and calculation precision according to landform data of the submerged area of the reservoir area of the barrier lake, and tracking and calculating upper and lower boundaries of the submerged area of the barrier lake at the current water level by adopting an iterative calculation method; precision parameters are calculated according to the reservoir capacity curve, iteration is conducted on all water level elevations below the water surface, the grid volume sum between different water levels is calculated according to the number of regional strips where each layer of grids intersect with the river channel satellite DEM, and the reservoir capacity water level curve of the barrier lake is obtained. A space subdivision and iteration intersection method is adopted, the reservoir capacity water level line is obtained, submerging range polygon calculation is conducted, the calculation efficiency and accuracy are high, and the method is suitable for being applied to the emergency rescue scene of the barrier lake.

Description

Technical Field

[0001] The invention relates to the field of reservoir dispatching calculation, and in particular to a method for quickly calculating a reservoir capacity curve. Background Art

[0002] After a landslide occurs, it is necessary to quickly calculate the reservoir capacity and reservoir curve in the shortest time to calculate key indicators such as reservoir capacity, breach losses, and upstream flooding range. However, since landslide lakes are often located in remote, sparsely populated, and high mountain valleys, there is usually a lack of comprehensive perception data, especially data such as underwater topography and hydrological sections. This will seriously affect the calculation speed of reservoir capacity and water level, and is a key factor restricting the accuracy of emergency rescue. The existing traditional methods for calculating reservoir capacity and water level of landslide lakes, such as hydrological extrapolation and volume calculation, have the defects of complex data acquisition, low calculation efficiency, and low accuracy. Summary of the invention

[0003] In view of the above-mentioned deficiencies in the prior art, the present invention provides a method for quickly calculating the reservoir capacity curve of a barrier lake, thereby realizing high-precision calculation of the barrier reservoir capacity water level curve, so as to solve the defects of the prior art such as difficulty in data acquisition, low calculation efficiency and accuracy.

[0004] In order to achieve the above-mentioned object of the invention, the technical solution adopted by the present invention is:

[0005] A method for quickly calculating a reservoir capacity curve of a barrier lake comprises the following steps:

[0006] S1. Obtain satellite DEM or DSM terrain data of the river channel in the dry season of the barrier lake reservoir area, add bathymetric points and hydrological section data, perform multi-source heterogeneous data fusion, and obtain the terrain data of the submerged area of ​​the barrier lake reservoir area;

[0007] S2. Obtain water level elevation and calculation accuracy based on topographic data of the submerged area of ​​the barrier lake reservoir, divide the space below the dam crest of the barrier lake reservoir into uniform hexahedrons using a three-dimensional spatial division method, and use an iterative calculation method to track and calculate the upper and lower boundaries of the submerged area of ​​the barrier lake under the current water level in combination with the ridge line and terrain elevation;

[0008] S3. Calculate the accuracy parameters according to the reservoir capacity curve, iterate all water level elevations below the water surface, mark the position of each grid at the current water level based on the water level elevation below the dam crest, and calculate the volume and grid volume between different water levels based on the number of regional strips where each layer of grid intersects with the river satellite DEM to obtain the reservoir capacity water level curve of the barrier lake.

[0009] Furthermore, step S1 specifically includes:

[0010] S11. Obtain satellite DEM or DSM terrain data of the river channel during the dry season at the location of the barrier lake, with an accuracy of more than 10m;

[0011] S12, converting the DEM or DSM terrain data into a spatial three-dimensional point set, and merging the hydrological section and temporary bathymetric point data into the spatial three-dimensional point set;

[0012] S13, constructing a B-spline surface using all the merged three-dimensional point sets as control points, and removing points outside the B-spline surface;

[0013] S14. Afterwards, the B-spline surface is converted back into preliminary DEM data of the bottom topography of the barrier lake through spatial interpolation and data sampling. The preliminary DEM data is used as the topographic data of the submerged area of ​​the barrier lake reservoir.

[0014] Furthermore, step S2 specifically includes:

[0015] S21, constructing a spatial three-dimensional coordinate system, obtaining the area where the water level elevation of the barrier lake is located according to the topographic data of the submerged area of ​​the barrier lake reservoir, and calculating the number of strips by using the strip method; dividing the area where the water level elevation of the barrier lake is located into a number of strips, and calculating the area M contained in each strip according to the total area S of the area where the water level elevation of the barrier lake is located and the number N of strips;

[0016]

[0017] Among them, a1 and b1 represent the abscissa and ordinate of the first vertex of the irregular area in the area where the water level of the barrier lake is located, respectively; a2 and b2 represent the abscissa and ordinate of the second vertex of the irregular area in the area where the water level of the barrier lake is located, respectively. n-1 、b n-1 They represent the horizontal and vertical coordinates of the n-1th vertex of the irregular area in the area where the water level of the barrier lake is located, respectively. n 、b n They represent the horizontal and vertical coordinates of the nth vertex of the irregular area in the area where the water level of the barrier lake is located, respectively. n+1 、b n+1 They respectively represent the abscissa and ordinate of the n+1th vertex of the irregular area in the area where the water level of the barrier lake is located;

[0018] S22, determining the maximum and minimum values ​​of the coordinates of the irregular area in the y direction;

[0019] S23, taking the area between the maximum value and the minimum value of the irregular area in the y direction as the scanning area, and using the numerical differentiation method to take the y-axis direction as the main stepping direction of the scanning, and calculating the abscissa of each intersection of the ordinate and the scanning area;

[0020] S24, according to the calculated ordinate on the irregular area and the abscissa of the intersection of the scanned area, filter out the grids between the two intersections corresponding to the same ordinate in the cut-and-fill area, and execute step S25;

[0021] S25. Obtain the coordinates of the two intersection points and establish a binary index of the coordinates of the two intersection points, namely:

[0022] {y i ,(x n ,x n+1 )}

[0023] Among them, y i represents the i-th ordinate, x n Indicates the i-th ordinate y i The horizontal coordinate of the nth intersection point, x n+1 Indicates the i-th ordinate y i The horizontal coordinate of the coordinates of the n+1th intersection point, n and n+1 represent two intersection points respectively;

[0024] S26, according to the binary index of the two intersection coordinates established in step S25, obtain the fill-in and cut-out area, wherein each binary index {y i ,(x n ,x n+1 )} represents the horizontal coordinate x corresponding to the i-th vertical coordinate n to x n+1 The points in between are flooded areas;

[0025] S27. Define y m Indicates the upper boundary of the mth strip, calculates the ordinate y m The sum of the flooded areas of the grid numbers contained in the corresponding binary index S m , if S m If it is smaller than the flooded area contained in each strip, S m The inundated area corresponding to the ordinate of the next row is accumulated, otherwise, the lower boundary coordinates of each strip m and the upper boundary coordinates of strip m+1 are obtained;

[0026] S28. Repeat steps S21-S27, and obtain the lower boundary of each grid layer from the water surface to the water bottom in the direction perpendicular to the water surface through iterative calculation, so as to obtain the upper and lower boundaries of the submerged area of ​​the barrier lake under the current water level.

[0027] Furthermore, step S3 specifically includes:

[0028] S31, establishing a flooding polygon according to the upper and lower boundaries of the flooded area of ​​the barrier lake under the current water level obtained in step S2, and using an iterative calculation method to calculate the current flooding range of the reservoir water level to obtain a flooding range polygon; specifically including:

[0029] S311. According to the reservoir water level line in the topographic data of the submerged area of ​​the barrier lake reservoir, the number of layers is divided in the three directions of x, y and z within the two cross-sectional spaces of the river channel, wherein the number of layers in the x direction within the two cross-sectional spaces of the river channel is a, the number of layers in the y direction within the two cross-sectional spaces of the river channel is b, and the number of layers in the z direction within the two cross-sectional spaces of the river channel is c;

[0030] S312, re-mark all grids of layer c, and iterate all grids of layer c according to rows and columns, with the iteration direction being in the vertical direction of the river channel, and find two or more boundary points by looping each row to obtain the coordinates of the boundary points;

[0031] S313, connecting the found boundary points using a boundary tracing method, and obtaining a flooded range polygon by distinguishing the area of ​​the river surface;

[0032] S32, using the X-scan line method to trace from a boundary point of the flooded range polygon to eight grid directions of the surrounding grids, and after tracing in one direction is closed, a separate cross-section polygon is formed, wherein the cross-section polygon is the internal cavity of the water surface of the river channel, and a new flooded range polygon is obtained by iterating all cross-section polygons;

[0033] S33, calculating the number of grids below the new flood range polygon and above the underwater terrain, and obtaining the reservoir capacity value of the current water level elevation;

[0034] S34. According to the calculation accuracy, iterate all water level values ​​to the bottom elevation position of the barrier dam to obtain the reservoir capacity values ​​of all water level elevations, and use a smooth curve to connect the water level elevation combination to obtain the reservoir capacity water level curve.

[0035] Furthermore, the criteria for determining the boundary points of the submerged area of ​​the barrier lake under the water level in step S3 are:

[0036] When the previous grid and the next grid of the cycle are at different positions above water, underground and in the water body of the river, the boundary points are determined and marked as 1, 2, and 3, respectively, where 1 represents the underwater boundary point, 2 represents the underground boundary point, and 3 represents the water boundary point.

[0037] Furthermore, the two or more boundary points found in step S312 are divided into three cases, namely, 2 intersection points, 2*n intersection points, and 3+2*n intersection points, specifically:

[0038] When two intersection points are found, it means there is no island between the two points on both sides of the river, and the river has no forks or big bends. At this time, the marks are 1, 2, and 1 respectively;

[0039] When 2*n intersections are found, it indicates that there is an island in the middle of the polygonal area of ​​the flooded range of the river. According to the distribution of boundary points, the boundary of the river contains the boundary of the island. At this time, the marks are 1, 2, 1, ..., 2, 1, and they are divided into inner boundary and outer boundary. The two intersections of the inner boundary form the boundary of an island, and the grid between the boundaries uses the current medium. At the same time, the content of the current medium is obtained in the x and y directions.

[0040] When 3+2*n intersection points are found, it indicates that the river has a large turn or branch, which is divided into upward and downward situations, marked as 2, 1, 2, 1.

[0041] The present invention has the following beneficial effects:

[0042] A method for quickly calculating a reservoir capacity curve proposed in the present invention uses satellite river digital elevation model data in the dry season as basic terrain data, updates the river terrain by combining temporary bathymetry and hydrological section data, obtains the reservoir capacity water level line, and performs polygonal calculation of the flooded range. The method has simple data acquisition and high calculation efficiency and accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 A schematic diagram of a flow chart of a method for quickly calculating a storage capacity curve proposed by the present invention;

[0044] Figure 2 This is a data diagram of the cut-and-fill area within the strip;

[0045] Figure 3 Schematic diagram of the process of scanning boundary points using the boundary tracking method.

[0046] Figure 4 This is a grid diagram of the river cross section.

[0047] Figure 5 A grid map of the river boundary.

[0048] Figure 6 Grid map of the river flooding area.

[0049] Figure 7 This is a river ridgeline grid map.

[0050] Figure 8 It is a terrain elevation grid map.

[0051] Fig. 9 Schematic diagram of underwater, underground and water boundary points.

[0052] Fig.10 This is the calculated reservoir water level curve. DETAILED DESCRIPTION

[0053] The specific implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0054] like Figure 1 As shown, a method for quickly calculating a storage capacity curve includes the following steps S1-S3:

[0055] S1. Obtaining satellite DEM or DSM terrain data of the river channel in the dry season of the barrier lake reservoir area, adding sounding points and hydrological section data, and fusing multi-source heterogeneous data to obtain terrain data of the submerged area of ​​the barrier lake reservoir area; Step S1 specifically includes:

[0056] S11. Obtain satellite DEM or DSM terrain data of the river channel during the dry season at the location of the barrier lake, with an accuracy of more than 10m;

[0057] S12, converting the DEM or DSM terrain data into a spatial three-dimensional point set, and merging the hydrological section and temporary bathymetric point data into the spatial three-dimensional point set;

[0058] S13, constructing a B-spline surface using all the merged three-dimensional point sets as control points, and removing points outside the B-spline surface;

[0059] S14. Afterwards, the B-spline surface is converted back into preliminary DEM data of the bottom topography of the barrier lake through spatial interpolation and data sampling. The preliminary DEM data is used as the topographic data of the submerged area of ​​the barrier lake reservoir.

[0060] S2. According to the terrain data of the submerged area of ​​the barrier lake reservoir, the water level elevation and calculation accuracy are obtained. The space below the dam top of the barrier lake reservoir is divided into uniform hexahedrons using the three-dimensional spatial division method, and combined with the ridge line and terrain undulation elevation, such as Figure 7 and Figure 8 As shown, an iterative calculation method is used to track and calculate the upper and lower boundaries of the submerged area of ​​the barrier lake under the current water level; step S2 specifically includes:

[0061] S21, constructing a spatial three-dimensional coordinate system, obtaining the area where the water level elevation of the barrier lake is located according to the topographic data of the submerged area of ​​the barrier lake reservoir, and calculating the number of strips by using the strip method; dividing the area where the water level elevation of the barrier lake is located into a number of strips, and calculating the area M contained in each strip according to the total area S of the area where the water level elevation of the barrier lake is located and the number N of strips;

[0062]

[0063] Among them, a1 and b1 represent the abscissa and ordinate of the first vertex of the irregular area in the area where the water level of the barrier lake is located, respectively; a2 and b2 represent the abscissa and ordinate of the second vertex of the irregular area in the area where the water level of the barrier lake is located, respectively. n-1 , b n-1 They represent the horizontal and vertical coordinates of the n-1th vertex of the irregular area in the area where the water level of the barrier lake is located, respectively. n , b n They represent the horizontal and vertical coordinates of the nth vertex of the irregular area in the area where the water level of the barrier lake is located, respectively. n+1 , b n+1 They respectively represent the abscissa and ordinate of the n+1th vertex of the irregular area in the area where the water level of the barrier lake is located;

[0064] like Figure 2 As shown, Figure 2 It is a data diagram of the fill and cut area within the strip; since the data reading and writing amount and the calculation amount are basically equal for the river areas with equal area, load balancing can be achieved by dividing the areas with equal area. Therefore, this embodiment divides the raster data into strips with equal fill and cut areas, and only analyzes the data of the fill and cut areas in the strips to obtain the reservoir water level line.

[0065] S22, determining the maximum and minimum values ​​of the coordinates of the irregular area in the y direction;

[0066] S23, taking the area between the maximum value and the minimum value of the irregular area in the y direction as the scanning area, and using the numerical differentiation method to take the y-axis direction as the main stepping direction of the scanning, calculate the abscissa of each intersection of the ordinate and the scanning area; Figure 3 As shown, a schematic diagram of the process of scanning boundary points using the boundary tracking method is shown;

[0067] S24, according to the calculated ordinate on the irregular area and the abscissa of the intersection of the scanned area, filter out the grids between the two intersections corresponding to the same ordinate in the cut-and-fill area, and execute step S25;

[0068] S25. Obtain the coordinates of the two intersection points and establish a binary index of the coordinates of the two intersection points, namely:

[0069] {y i ,(x n ,x n+1 )}

[0070] Among them, y i represents the i-th ordinate, x n Indicates the i-th ordinate y i The horizontal coordinate of the nth intersection point, x n+1Indicates the i-th ordinate y i The horizontal coordinate of the coordinates of the n+1th intersection point, n and n+1 represent two intersection points respectively;

[0071] S26, according to the binary index of the two intersection coordinates established in step S25, obtain the fill-in and cut-out area, wherein each binary index {y i ,(x n ,x n+1 )} represents the horizontal coordinate x corresponding to the i-th vertical coordinate n to x n+1 The points in between are flooded areas;

[0072] S27. Define y m Indicates the upper boundary of the mth strip, calculates the ordinate y m The sum of the flooded areas of the grid numbers contained in the corresponding binary index S m , if S m If it is smaller than the flooded area contained in each strip, S m The inundated area corresponding to the ordinate of the next row is accumulated, otherwise, the lower boundary coordinates of each strip m and the upper boundary coordinates of strip m+1 are obtained;

[0073] S28. Repeat steps S21-S27, and obtain the lower boundary of each grid layer from the water surface to the water bottom in the direction perpendicular to the water surface through iterative calculation, so as to obtain the upper and lower boundaries of the submerged area of ​​the barrier lake under the current water level.

[0074] S3. According to the calculation accuracy parameters of the reservoir capacity curve, all water level elevations below the water surface are iterated. Based on the water level elevation below the dam crest, the position of each grid at the current water level is marked. According to the number of regional strips where each layer of grid intersects with the river channel satellite DEM, the volume and grid volume between different water levels are calculated to obtain the reservoir capacity water level curve of the barrier lake, such as Fig.10 Step S3 specifically includes:

[0075] S31, establishing a flooding polygon according to the upper and lower boundaries of the flooded area of ​​the barrier lake under the current water level obtained in step S2, and using an iterative calculation method to calculate the current flooding range of the reservoir water level to obtain a flooding range polygon; specifically including:

[0076] S311. According to the reservoir water level line in the topographic data of the submerged area of ​​the barrier lake reservoir, the number of layers is divided in the three directions of x, y and z within the two cross-sectional spaces of the river channel, wherein the number of layers in the x direction within the two cross-sectional spaces of the river channel is a, the number of layers in the y direction within the two cross-sectional spaces of the river channel is b, and the number of layers in the z direction within the two cross-sectional spaces of the river channel is c;

[0077] S312, re-mark all grids of layer c, and iterate all grids of layer c according to rows and columns, with the iteration direction being in the vertical direction of the river channel, and find two or more boundary points by looping each row to obtain the coordinates of the boundary points;

[0078] S313, the boundary points found are connected by using the boundary tracing method, and the flooding range polygon is obtained by distinguishing the area of ​​the river surface, such as Figure 5 and Figure 6 shown.

[0079] S32, using the X-scan line method to trace from a boundary point of the flooded range polygon to eight grid directions of the surrounding grids, and after tracing in one direction is closed, a separate cross-section polygon is formed, wherein the cross-section polygon is the internal cavity of the water surface of the river channel, and a new flooded range polygon is obtained by iterating all cross-section polygons;

[0080] S33, calculating the number of grids below the new flood range polygon and above the underwater terrain, and obtaining the reservoir capacity value of the current water level elevation;

[0081] S34. According to the calculation accuracy, iterate all water level values ​​to the bottom elevation position of the barrier dam to obtain the reservoir capacity values ​​of all water level elevations, and use a smooth curve to connect the water level elevation combination to obtain the reservoir capacity water level curve.

[0082] The criteria for determining the boundary points of the submerged area of ​​the barrier lake under the water level in step S3 are:

[0083] like Figure 4 As shown, when the previous grid of the cycle and the next grid are at different positions above water, underground and in the water body of the river, the boundary points are determined and marked as 1, 2, and 3, respectively, where 1 represents an underwater boundary point, 2 represents an underground boundary point, and 3 represents a water boundary point.

[0084] The two or more boundary points found in step S312 are divided into three cases, namely, 2 intersection points, 2*n intersection points, and 3+2*n intersection points. Fig. 9 As shown, specifically:

[0085] When two intersection points are found, it means there is no island between the two points on both sides of the river, and the river has no forks or big bends. At this time, the marks are 1, 2, and 1 respectively;

[0086] When 2*n intersections are found, it indicates that there is an island in the middle of the polygonal area of ​​the flooded range of the river. According to the distribution of boundary points, the boundary of the river contains the boundary of the island. At this time, the marks are 1, 2, 1, ..., 2, 1, and they are divided into inner boundary and outer boundary. The two intersections of the inner boundary form the boundary of an island, and the grid between the boundaries uses the current medium. At the same time, the content of the current medium is obtained in the x and y directions.

[0087] When 3+2*n intersections are found, it indicates that the river channel has a large bend or branch, which is divided into two situations: upward and downward, marked as 2, 1, 2, 1. Each boundary point is connected to form multiple flooding ranges.

[0088] A method for quickly calculating a reservoir capacity curve proposed in the present invention uses satellite river digital elevation model data in the dry season as basic terrain data, updates the river terrain by combining temporary bathymetry and hydrological section data, obtains the reservoir capacity water level line, and performs polygonal calculation of the flooded range. The method has simple data acquisition and high calculation efficiency and accuracy.

[0089] The present invention uses specific embodiments to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core idea. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as a limitation on the present invention.

[0090] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the principles of the present invention, and should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific variations and combinations that do not deviate from the essence of the present invention based on the technical revelations disclosed by the present invention, and these variations and combinations are still within the protection scope of the present invention.

Claims

1. A rapid calculation method for reservoir capacity curve of barrier lake, characterized in that: The following steps are involved: S1. Obtain satellite DEM or DSM terrain data of the river channel in the dry season of the barrier lake reservoir area, add bathymetric points and hydrological section data, perform multi-source heterogeneous data fusion, and obtain the terrain data of the submerged area of ​​the barrier lake reservoir area; S2. Obtain water level elevation and calculation accuracy based on topographic data of the submerged area of ​​the barrier lake reservoir, divide the space below the dam crest of the barrier lake reservoir into uniform hexahedrons using a three-dimensional spatial division method, and use an iterative calculation method to track and calculate the upper and lower boundaries of the submerged area of ​​the barrier lake under the current water level in combination with the ridge line and terrain elevation; S3. Calculate the accuracy parameters according to the reservoir capacity curve, iterate all water level elevations below the water surface, mark the position of each grid at the current water level based on the water level elevation below the dam crest, and calculate the volume and grid volume between different water levels based on the number of regional strips where each layer of grid intersects with the river satellite DEM to obtain the reservoir capacity water level curve of the barrier lake.

2. A method for rapid calculation of storage capacity curve according to claim 1, characterized in that: The step S1 specifically includes: S11. Obtain satellite DEM or DSM terrain data of the river channel during the dry season at the location of the barrier lake, with an accuracy of more than 10m; S12, converting the DEM or DSM terrain data into a spatial three-dimensional point set, and merging the hydrological section and temporary bathymetric point data into the spatial three-dimensional point set; S13, constructing a B-spline surface using all the merged three-dimensional point sets as control points, and removing points outside the B-spline surface; S14. Afterwards, the B-spline surface is converted back into preliminary DEM data of the bottom topography of the barrier lake through spatial interpolation and data sampling. The preliminary DEM data is used as the topographic data of the submerged area of ​​the barrier lake reservoir.

3. A method for rapid calculation of storage capacity curve according to claim 1, characterized in that: The step S2 specifically includes: S21, constructing a spatial three-dimensional coordinate system, obtaining the area where the water level elevation of the barrier lake is located according to the topographic data of the submerged area of ​​the barrier lake reservoir, and calculating the number of strips by using the strip method; dividing the area where the water level elevation of the barrier lake is located into a number of strips, and calculating the area M contained in each strip according to the total area S of the area where the water level elevation of the barrier lake is located and the number N of strips; Among them, a1 and b1 represent the abscissa and ordinate of the first vertex of the irregular area in the area where the water level of the barrier lake is located, respectively; a2 and b2 represent the abscissa and ordinate of the second vertex of the irregular area in the area where the water level of the barrier lake is located, respectively. n-1 、b n-1 They represent the horizontal and vertical coordinates of the n-1th vertex of the irregular area in the area where the water level of the barrier lake is located, respectively. n 、b n They represent the horizontal and vertical coordinates of the nth vertex of the irregular area in the area where the water level of the barrier lake is located, respectively. n+1 、b n+1 They respectively represent the abscissa and ordinate of the n+1th vertex of the irregular area in the area where the water level of the barrier lake is located; S22, determining the maximum and minimum values ​​of the coordinates of the irregular area in the y direction; S23, taking the area between the maximum value and the minimum value of the irregular area in the y direction as the scanning area, and using the numerical differentiation method to take the y-axis direction as the main stepping direction of the scanning, and calculating the abscissa of each intersection of the ordinate and the scanning area; S24, according to the calculated ordinate on the irregular area and the abscissa of the intersection of the scanned area, filter out the grids between the two intersections corresponding to the same ordinate in the cut-and-fill area, and execute step S25; S25. Obtain the coordinates of the two intersection points and establish a binary index of the coordinates of the two intersection points, namely: {y i ,(x n ,x n+1 )} Among them, y i represents the i-th ordinate, x n Indicates the i-th ordinate y i The horizontal coordinate of the nth intersection point, x n+1 Indicates the i-th ordinate y i The horizontal coordinate of the coordinates of the n+1th intersection point, n and n+1 represent two intersection points respectively; S26, according to the binary index of the two intersection coordinates established in step S25, obtain the fill-in and cut-out area, wherein each binary index {y i ,(x n , x n+1 )} represents the horizontal coordinate x corresponding to the i-th vertical coordinate n to x n+1 The points in between are flooded areas; S27. Define y m Indicates the upper boundary of the mth strip, calculates the ordinate y m The sum of the flooded areas of the grid numbers contained in the corresponding binary index S m , if S m If it is smaller than the flooded area contained in each strip, S m The inundated area corresponding to the ordinate of the next row is accumulated, otherwise, the lower boundary coordinates of each strip m and the upper boundary coordinates of strip m+1 are obtained; S28. Repeat steps S21-S27, and obtain the lower boundary of each grid layer from the water surface to the water bottom in the direction perpendicular to the water surface through iterative calculation, so as to obtain the upper and lower boundaries of the submerged area of ​​the barrier lake under the current water level.

4. The method for rapid calculation of subsurface spatial grid of barrier lake according to claim 2, characterized in that: Step S3 specifically includes: S31, establishing a flooding polygon according to the upper and lower boundaries of the flooded area of ​​the barrier lake under the current water level obtained in step S2, and using an iterative calculation method to calculate the current flooding range of the reservoir water level to obtain a flooding range polygon; specifically including: S311. According to the reservoir water level line in the topographic data of the submerged area of ​​the barrier lake reservoir, the number of layers is divided in the three directions of x, y and z within the two cross-sectional spaces of the river channel, wherein the number of layers in the x direction within the two cross-sectional spaces of the river channel is a, the number of layers in the y direction within the two cross-sectional spaces of the river channel is b, and the number of layers in the z direction within the two cross-sectional spaces of the river channel is c; S312, re-mark all grids of layer c, and iterate all grids of layer c according to rows and columns, with the iteration direction being in the vertical direction of the river channel, and find two or more boundary points by looping each row to obtain the coordinates of the boundary points; S313, connecting the found boundary points using a boundary tracing method, and obtaining a flooded range polygon by distinguishing the area of ​​the river surface; S32, using the X-scan line method to trace from a boundary point of the flooded range polygon to eight grid directions of the surrounding grids, and after tracing in one direction is closed, a separate cross-section polygon is formed, wherein the cross-section polygon is the internal cavity of the water surface of the river channel, and a new flooded range polygon is obtained by iterating all cross-section polygons; S33, calculating the number of grids below the new flood range polygon and above the underwater terrain, and obtaining the reservoir capacity value of the current water level elevation; S34. According to the calculation accuracy, iterate all water level values ​​to the bottom elevation position of the barrier dam to obtain the reservoir capacity values ​​of all water level elevations, and use a smooth curve to connect the water level elevation combination to obtain the reservoir capacity water level curve.

5. A method for rapid calculation of storage capacity curve according to claim 3, characterized in that: The criteria for determining the boundary points of the submerged area of ​​the barrier lake under the water level in step S3 are: When the previous grid and the next grid of the cycle are at different positions above water, underground and in the water body of the river, the boundary points are determined and marked as 1, 2, and 3, respectively, where 1 represents the underwater boundary point, 2 represents the underground boundary point, and 3 represents the water boundary point.

6. A method for rapid calculation of storage capacity curve according to claim 4, characterized in that: The two or more boundary points found in step S312 are divided into three cases, namely, 2 intersection points, 2*n intersection points, and 3+2*n intersection points, specifically: When two intersection points are found, it means there is no island between the two points on both sides of the river, and the river has no forks or big bends. At this time, the marks are 1, 2, and 1 respectively; When 2*n intersections are found, it indicates that there is an island in the middle of the polygonal area of ​​the flooded range of the river. According to the distribution of boundary points, the boundary of the river contains the boundary of the island. At this time, the marks are 1, 2, 1, ..., 2, 1, and they are divided into inner boundary and outer boundary. The two intersections of the inner boundary form the boundary of an island, and the grid between the boundaries uses the current medium. At the same time, the content of the current medium is obtained in the x and y directions. When 3+2*n intersections are found, it indicates that the river has a large bend or branch, which is divided into two situations: upward and downward, marked as 2, 1, 2, 1. Finally, the boundary points marked are connected to obtain the water surface boundary range.

Images