Method for optimizing rain gauge sites of mountain torrent disaster based on grid coverage

Abstract

The application discloses a mountain torrent disaster rainfall station optimization method based on grid coverage, belongs to the technical field of mountain torrent disaster research, and comprises the following steps: S1. extracting a grid related to an upstream of a mountain torrent disaster prevention object; S2. calculating a dynamic monitoring radius and a coverage area; S3. spatial intersection analysis; S4. calculating a coverage of the mountain torrent disaster prevention object and calculating a regional comprehensive coverage; and S5. rainfall station network layout optimization. The application solves the problems of insufficient coverage of an existing rainfall station network and failure to combine the characteristics of a basin, improves the scientificity of coverage evaluation, optimization efficiency and scheme economy, and is suitable for rainfall station network planning of a small-basin mountain torrent disaster prevention area.

Description

Technical Field

[0001] This invention belongs to the field of flash flood disaster research technology, specifically involving an optimization method for flash flood disaster rainfall stations based on grid coverage. Background Technology

[0002] Flash floods specifically refer to sudden floods that occur in hilly areas due to heavy rainfall, characterized by their suddenness and destructive power. Flash floods in small watersheds directly endanger the lives and property of residents, and their frequent occurrence has become a prominent issue that urgently needs to be addressed in current disaster prevention and mitigation efforts. A well-planned rain gauge network is crucial for obtaining regional characteristics of rainstorms and floods and improving the accuracy of flash flood monitoring and early warning systems.

[0003] Although numerous methods exist for rain gauge deployment, such as spatial interpolation, statistical analysis, and entropy theory, and domestic scholars have proposed methods like the cumulative difference method, the storm center control zone method, the cone method, the watershed hydrological model method, and the sampling method, methods for assessing rain gauge network coverage based on the characteristics of flash flood disaster prevention targets in small watersheds are still scarce. Therefore, a comprehensive solution is needed for evaluating small watershed rain gauge networks specifically for flash flood disaster prevention targets, aiming to accurately and efficiently optimize rain gauge network layout and enabling its widespread application in practical engineering projects. Summary of the Invention

[0004] To address the aforementioned issues, embodiments of the present invention propose an optimization method for rainfall stations in flash flood disasters based on grid coverage, which solves the problems of insufficient coverage of rainfall station networks and inability to incorporate watershed characteristics of the study area in existing technologies.

[0005] The present invention provides a method for optimizing rainfall stations for flash flood disasters based on grid coverage, comprising the following steps:

[0006] S1. Extraction of upstream associated rasters for flash flood disaster prevention objects: Create a direction mapping dictionary, use the reverse index of the direction mapping dictionary, and use a breadth-first search algorithm to extract all upstream rasters for each flash flood disaster prevention object;

[0007] S2. Calculation of dynamic monitoring radius and coverage area: Divide each mountain flood disaster prevention target into 8 sector sub-regions, calculate the slope standard deviation of each sector sub-region and the monitoring radius of each sector direction, and combine the monitoring radii of the 8 sector sub-regions to obtain the total coverage area of ​​each rain gauge station;

[0008] S3. Spatial Intersection Analysis: Based on the sector dynamic monitoring radius obtained in S2, construct its dynamic monitoring range polygon; create an R-tree spatial index to efficiently calculate the spatial intersection relationship between the dynamic monitoring range of each rain gauge station and the upstream area of ​​the flash flood disaster prevention object extracted in S1, and query and record the intersecting raster data and the corresponding type information;

[0009] S4. Calculation of coverage of mountain torrent disaster prevention targets and calculation of regional comprehensive coverage: Based on the spatial intersection results obtained in S3, the coverage calculation formula considering the attenuation coefficient is used to calculate the coverage of the current rain gauge network for mountain torrent disaster prevention targets and the regional comprehensive coverage.

[0010] S5. Rain gauge network layout optimization: Optimize the network round by round according to the optimization model. Under the premise of meeting the preset comprehensive coverage target, minimize the number of rain gauges to be retained. At the same time, find the solution that maximizes the watershed coverage in the optimal round and output the optimal rain gauge layout scheme under the target coverage.

[0011] The specific steps for extracting upstream associated raster data for flood disaster prevention targets in S1 are as follows: A direction mapping dictionary is pre-created for the flow direction values ​​of the D8 algorithm, converting the D8 algorithm flow direction values ​​into corresponding row and column coordinate offsets; flow direction raster data is imported, the flow direction value of the central cell of the flood disaster prevention target is read, a reverse index is used with the direction mapping dictionary, and then a breadth-first search algorithm is used to extract all upstream raster data for each flood disaster prevention target. The row and column coordinates, geographic coordinates, and flood disaster prevention target type of the raster data are stored in the database; this achieves the extraction of upstream catchment area raster data for any point source.

[0012] The construction steps of the direction mapping dictionary reverse index are as follows: First, map the coordinates of the flood disaster prevention target points to the raster row and column numbers of the DEM to determine their positions in the raster matrix; create a direction mapping dictionary in advance for the flow direction values ​​of the D8 algorithm, converting the flow direction values ​​of the D8 algorithm into corresponding row and column coordinate offsets; when extracting upstream rasters, scan the flow direction rasters, obtain the flow direction value of their center cells, and then convert the direction through the direction mapping dictionary, directly performing calculations on the row and column coordinates to calculate the position of the downstream cells; record the current cell in the reverse index of the downstream cells, repeat this operation until a mapping relationship of "downstream cells → all upstream cells" is established, thereby constructing a reverse index dictionary of "target raster to all source rasters".

[0013] The calculation of the dynamic monitoring radius and coverage area of ​​S2 specifically refers to dividing the monitoring area around each mountain flood disaster prevention target into 8 fan-shaped sub-regions with a radius of 5 kilometers at 45°, taking each target as the center; based on the slope grid data, reading the slope values ​​of all grids in each fan-shaped sub-region, and further calculating the slope standard deviation of the region to quantify the terrain complexity in each direction; obtaining the monitoring radius of each fan-shaped direction based on the obtained slope standard deviation of each fan-shaped region; and then combining the monitoring radii of the 8 fan-shaped regions to obtain the total coverage area of ​​each rain gauge station.

[0014] The spatial intersection analysis in S3 includes two stages: creating a spatial index and spatial querying. In the spatial index creation stage, an independent R-tree spatial index is established for the upstream raster set of each disaster prevention object, using the hierarchical structure of the R-tree: "root node → non-leaf node → leaf node". In the spatial querying stage, candidate rasters are first quickly queried in the R-tree using the bounding rectangle of the sector buffer. Then, the ray method is used to make a precise geometric relationship judgment on the candidate rasters to determine which rasters are actually located within the sector. Finally, all rasters within the sector are retained to form an "overlapping raster set".

[0015] The classification method for mountain flood disaster prevention objects in S4 is as follows: an automatic identification and classification method for disaster prevention object types is adopted. The surface vector data and DEM raster data of the built-up area of ​​the mountain flood disaster prevention object are intersected to obtain the comprehensive raster data of the built-up area of ​​the mountain flood disaster prevention object. The catchment area and spatial location are calculated based on this data, and the mountain flood disaster prevention objects are divided into river confluence type disaster prevention objects and slope overflow type disaster prevention objects.

[0016] The formula for calculating the coverage of the river confluence-type disaster prevention targets is as follows:

[0017]

[0018] In the formula, To ensure coverage of disaster prevention targets in river confluence areas, G represents the number of grid cells where the catchment area intersects with the coverage area of ​​the rain gauge, and G represents the number of grid cells upstream of the catchment area at or above the cross section of the flash flood disaster prevention target.

[0019] The formula for calculating the coverage of the slope-flow type disaster prevention object is as follows:

[0020]

[0021] In the formula, denoted as the coverage of the slope flooding disaster prevention object, d is the straight-line distance from the slope flooding disaster prevention object to the nearest rain gauge station, r is the monitoring radius of the nearest rain gauge station of the slope flooding disaster prevention object, D is the maximum influence distance of the target rain gauge station on the slope flooding disaster prevention object, and k is the attenuation coefficient.

[0022] Overall coverage The calculation formula is:

[0023]

[0024] In the formula, , The weights are respectively for river confluence-type disaster prevention targets and slope overflow-type disaster prevention targets. This represents the coverage of the i-th river confluence-type disaster prevention object. N represents the coverage of the i-th slope runoff type disaster prevention object. A N represents a subset of disaster prevention targets related to river confluence. B This represents a subset of disaster prevention targets for slope runoff, m A m is the number of river confluence-type disaster prevention targets. B This refers to the number of disaster prevention targets for slope runoff.

[0025] The objective function of the optimization model in S5 is:

[0026]

[0027] In the formula, Z represents the comprehensive target value for rain gauge network optimization, which is a comprehensive balancing indicator between the construction cost of the rain gauge network and the watershed coverage. Let the unit cost of the j-th rain gauge be... , indicating the retention status of the j-th rain gauge, where 1 indicates retention and 0 indicates removal. The economic costs of retaining rain gauges throughout their entire lifecycle were quantified. Indicates the overall coverage of the watershed. The coverage gap is represented by λ, which is a weighting coefficient used to balance the priority of rain gauge cost and coverage.

[0028] The constraints of the objective function include:

[0029] Overall coverage constraints:

[0030]

[0031] Constraints for dynamic monitoring radius calculation:

[0032]

[0033] Deployment state binary constraints:

[0034]

[0035] In the formula, To achieve the preset coverage target, , The weights are respectively for river confluence-type disaster prevention targets and slope overflow-type disaster prevention targets. This represents the coverage of the i-th river confluence-type disaster prevention object. N represents the coverage of the i-th slope runoff type disaster prevention object. A N represents a subset of disaster prevention targets related to river confluence. B This represents a subset of disaster prevention targets for slope runoff, m A m is the number of river confluence-type disaster prevention targets. B The number of flood-prone disaster prevention targets on slopes. This represents the dynamic monitoring radius of the k-th sector area at the j-th rain gauge station, in km. This represents the standard deviation of the slope in the k-th monitoring area at the j-th rain gauge station. Let $\frac{j}{k}$ be the standard deviation of the slope in the $k$-th monitoring area of ​​the $j$-th existing rain gauge station, where $k = 1, 2, ..., 8$, and the unit is °. This represents the retention status of the j-th existing rain gauge, where 1 = retained and 0 = removed. S represents the set of candidate rain gauge locations. R represents the dynamic monitoring radius R of the k-th sector of any rain gauge j. j,k Slope standard deviation T j,k All must satisfy the dynamic monitoring radius calculation constraint; for each rain gauge j in the candidate rain gauge set S, its deployment state It can only be 0 or 1.

[0036] The optimal solution of the optimization model is calculated as follows:

[0037] Starting with the initial set of all rain gauges, the search is performed iteratively round by round, gradually approaching the optimal solution of the objective function. In each round, one rain gauge is removed from the set of all rain gauges, and the watershed coverage under the current condition is calculated. If a feasible solution that meets the target coverage is found in the current round, the evaluation of the remaining combinations in the current round is terminated, and the calculation is carried out in the next round. If no feasible solution is found, the current round continues. When no feasible solution can be found in a certain round, the process backtracks to the previous round and selects the scheme with the highest coverage as the final optimization result, thereby calculating the optimal solution of the objective function of the optimization model.

[0038] The beneficial effects of this invention are that it constructs a complete, efficient, and accurate technical process for optimizing rain gauge networks for flash flood disasters based on coverage. This solves the problems of insufficient coverage and failure to incorporate watershed characteristics in existing rain gauge networks, significantly improving the scientific nature of coverage assessment, optimization efficiency, and economic efficiency of the scheme. It is applicable to rain gauge network planning in flash flood disaster prevention areas of small watersheds, providing reliable technical support for the scientific planning of rain gauge networks in flash flood disaster prevention areas. Attached Figure Description

[0039] Figure 1 This is a flowchart of the method for optimizing rainfall stations for flash flood disasters based on grid coverage, as proposed in this invention.

[0040] Figure 2 This is a schematic diagram of the elements of the Shiwen River Basin in this invention.

[0041] Figure 3 This is a schematic diagram illustrating the principle of a grid extraction method for upstream catchment areas of an arbitrary point source.

[0042] Figure 4A schematic diagram showing the raster extraction results of the upstream catchment areas of various mountain flood disaster prevention targets.

[0043] Figure 5 This is a schematic representation of the coverage of various mountain flood disaster prevention targets after optimization by the present invention.

[0044] Figure 6 This is a schematic diagram of the optimized rain gauge layout according to the present invention. Detailed Implementation

[0045] Embodiments of the present invention are described in detail below, examples of which are illustrated in the accompanying drawings. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0046] like Figure 1 As shown, the method for optimizing rainfall stations for flash flood disasters based on grid coverage according to the present invention includes the following steps:

[0047] S1. Upstream associated raster extraction of flash flood disaster prevention targets: First, map the coordinates (x, y) of the flash flood disaster prevention target points to the raster row and column numbers (row, col) of the DEM to determine their positions in the raster matrix; create a direction mapping dictionary in advance for the flow direction values ​​of the D8 algorithm, and convert the flow direction values ​​of the D8 algorithm into the corresponding row and column coordinate offsets; the direction mapping dictionary is shown in Table 1.

[0048] Table 1 Direction Mapping Dictionary

[0049]

[0050] When extracting upstream rasters, the flow direction rasters are scanned to obtain the flow direction value of their center cells. Then, the direction is transformed through the direction mapping dictionary, and the calculation is performed directly on the row and column coordinates to calculate the position of the downstream cells. The current cell is recorded in the inverted index of the downstream cells. This operation is repeated until a mapping relationship of "downstream cells → all upstream cells" is established, thereby constructing an inverted index dictionary from the target raster to all source rasters.

[0051] When processing raster-associated disaster prevention targets along river confluences, since the flood risk mainly originates from precipitation runoff in the upstream confluence area and the river's convergence process, it is necessary to statistically analyze the complete upstream raster region to improve optimization efficiency. Starting from this raster, a breadth-first search (BFS) algorithm is used to retrieve the set of upstream raster lines confluencing the flash flood disaster prevention target and store them in the database. Database fields are used to bind the upstream raster set to the raster associated with the flash flood disaster prevention target, forming a one-to-many mapping relationship. Starting from this raster, BFS is used to quickly retrieve all upstream raster sets confluencing this point from the inverted index and store them in the database. When processing raster-associated disaster prevention targets along slope runoffs, their flood risk is mainly affected by local slope runoff; therefore, only the local precipitation at the target's location needs to be considered, and only the raster containing the disaster prevention target needs to be extracted.

[0052] S2. Calculation of dynamic monitoring radius and coverage area:

[0053] Centered on each flash flood disaster prevention target, the surrounding monitoring area is divided into 8 fan-shaped sub-regions with a radius of 5 kilometers at 45°. Based on slope raster data, the slope values ​​of all grids in each fan-shaped sub-region are read, and the slope standard deviation of the region is further calculated to quantify the terrain complexity in each direction. The monitoring radius of each fan-shaped direction is obtained based on the obtained slope standard deviation of each fan-shaped region. On this basis, the total coverage area of ​​each rain gauge is obtained by combining the monitoring radii of the 8 fan-shaped regions.

[0054] S3. Spatial Intersection Analysis: Based on the sector dynamic monitoring radius obtained in S2, construct its dynamic monitoring range polygon; create an R-tree spatial index to efficiently calculate the spatial intersection relationship between the dynamic monitoring range of each rain gauge station and the upstream area of ​​the flash flood disaster prevention object extracted in S1, and query and record the intersecting raster data and the corresponding type information.

[0055] Spatial intersection analysis includes two stages: creating a spatial index and spatial querying.

[0056] During the spatial index creation phase, an independent R-tree spatial index is established for the upstream raster set of each flash flood disaster prevention object. The storage content of the root node, leaf nodes, and non-leaf nodes follows the hierarchical characteristics of the spatial index. Leaf nodes store the minimum boundary rectangle of the raster cell (used to map the actual raster data) and the identification information of the raster cell (used to associate with the upstream raster set of the flash flood disaster prevention object point and read the specific attribute data of the raster cell). Non-leaf nodes store the covering rectangle (i.e., the overall boundary range of all features (or child nodes) contained in the lower-level node (calculated by merging the MBR of the child nodes)) and the reference pointer of the lower-level node (a reference pointing to the next level node). The root node stores the minimum boundary matrix of all child nodes and the reference pointer of the child nodes. By adopting the hierarchical structure of R-tree "root node → non-leaf node → leaf node", the spatial information of massive raster points is aggregated and stored according to the boundary range, realizing fast filtering of the query range, efficient location of raster candidate sets that may overlap with the fan-shaped buffer, and significantly reducing the amount of computation for subsequent accurate judgment.

[0057] The spatial query phase consists of two steps: coarse filtering (using the bounding rectangle of the sector buffer to quickly query candidate rasters in the R-tree) and precise positioning (performing precise geometric relationship judgments on candidate rasters to determine which rasters are truly located within the sector). In the coarse filtering phase, since range queries using rectangles are extremely efficient in the R-tree (only requiring comparison of boundary coordinates), rasters that do not intersect the sector at all can be quickly filtered out using the bounding rectangle, significantly reducing the computational load for subsequent precise judgments. Therefore, the bounding rectangle of the sector buffer (the smallest rectangle that can completely enclose the sector polygon) is first created. The boundary coordinates are calculated based on parameters such as the center, radius, and angle range of the sector, and represented by a quadruple (min_x, min_y, max_x, max_y), representing the minimum and maximum x / y coordinates among all vertices of the sector, respectively. Then, the bounding rectangle quadruple is passed into the R-tree as the query condition, and starting from the root node, the aggregation boundary of the child nodes is checked layer by layer to see if it intersects with the bounding rectangle. If the nodes do not intersect, the subtree is directly excluded; if they intersect, the process continues to delve into its lower-level nodes until the leaf nodes are reached, ultimately returning all raster indices (candidate set) that meet the criteria. Because the circumscribed rectangle contains sectors, a raster intersecting the rectangle may be inside the sector or inside the rectangle but not within the sector, requiring precise judgment in the next stage. In the precise positioning stage, a ray casting method is used (radiating a ray from the center point of the raster in any direction, counting the number of intersections between the ray and the sector boundary; if the number is odd, the point is inside the polygon; if it is even, the point is outside the polygon) to verify whether the raster is inside the sector polygon. Finally, all rasteres inside the sector are retained, forming an "overlapping raster set".

[0058] S4. Calculation of coverage of mountain torrent disaster prevention objects and calculation of regional comprehensive coverage: Classify mountain torrent disaster prevention objects, substitute them into the formula according to different types of mountain torrent disaster prevention objects to calculate the coverage of mountain torrent disaster prevention objects, and calculate the regional comprehensive coverage based on the coverage of mountain torrent disaster prevention objects.

[0059] An automatic identification and classification method for disaster prevention object types was used to classify objects for flash flood disaster prevention. First, the built-up area vector file of the disaster prevention object and the DEM raster data of the study area were imported. Through spatial overlay analysis, the vector files were converted into a raster mask corresponding to the disaster prevention object, and all raster cells covering the disaster prevention object were extracted. The row and column numbers, geographic coordinates, and elevation values ​​of each raster were recorded. The total number of raster cells N for the disaster prevention object was counted, and the actual built-up area Area = N × raster resolution was calculated. 2 Using pre-extracted watershed river network vectors (or rasterized river networks), for each raster within the disaster prevention target area, calculate its Euclidean distance to the nearest linear feature of the river network, and take the minimum value d among all raster values. min The location of the nearest grid cell is recorded. Based on the constructed flow-direction grid and reverse index dictionary, for all grid cells within the built-up area of ​​the disaster prevention target, the upstream grid set of each grid cell is extracted by tracing upstream. Then, the union of these sets is taken to obtain the "all upstream grid set of the built-up area of ​​the disaster prevention target". The total number of grid cells M in this set multiplied by the grid area is the runoff area A. up (Unit: km²)

[0060] For d min And A up Make a judgment; if the spatial distribution is located within the 500m buffer zone of the river channel (d min (≤500m), catchment area generally ranges from 5 to 200km². 2 Between (5≤A) up ≤200km 2 If the location is outside the 500m buffer zone of the river channel, it is considered a target for flood prevention in mountainous areas prone to flash floods. These targets are primarily threatened by direct floods from streams and rivers. min >500m), the catchment area is generally between 0 and 5km². 2 Between (0≤A) up ≤5km 2 If the landslide is detected, it is classified as a surface runoff type. These landslide targets are primarily affected by surface runoff induced by short-duration heavy rainfall. If other conditions occur (such as d...), the landslide is classified as a surface runoff type. min ≤500m but A up <5; or d min Greater than 500 m but A up If ≥5), the name of the defending object and d will be output. min A upInformation is provided to help individuals determine the category of objects to be protected against flash floods.

[0061] The methods for determining the coverage of different categories of flash flood prevention targets are not entirely the same: For flash flood prevention targets with large catchment areas, the coverage needs to consider upstream and downstream relationships, and the overall coverage of the catchment area above the target section must be taken into account when calculating the coverage; for flash flood prevention targets with small catchment areas, the coverage only needs to consider the coverage of the catchment area where the target is located. The formula for calculating the coverage area ratio and coverage of each flash flood prevention target is as follows:

[0062] The formula for calculating the coverage of disaster prevention targets in river confluence areas is:

[0063]

[0064] In the formula, The coverage of river confluence-type disaster prevention targets, with a value ranging from 0 to 1. G represents the number of grid cells where the runoff area intersects with the coverage area of ​​the rain gauge, and G represents the number of grid cells upstream of the runoff area at or above the cross section of the flash flood disaster prevention target.

[0065] The formula for calculating the coverage of the slope-flow type disaster prevention object is as follows:

[0066]

[0067] In the formula, denoted as , where d is the straight-line distance from the slope flooding disaster prevention object to the nearest rain gauge station (km), r is the monitoring radius of the nearest rain gauge station for the slope flooding disaster prevention object (km), D is the maximum influence distance of the target rain gauge station on the slope flooding disaster prevention object (km), and k is the attenuation coefficient, representing the rate at which the influence coverage decreases with distance d.

[0068] when That is, when the target of flash flood disaster prevention is located within the monitoring range of the nearest rain gauge, the coverage is 1; when When the target of flash flood disaster prevention is outside the monitoring coverage of the rain gauge but not beyond the maximum impact distance, the coverage is calculated using the natural exponential decay function. This results in rapid decay near the boundary and slow decay far from the boundary. The monitoring impact weakens rapidly outside the core coverage area, and the decay rate gradually slows down as the distance increases. This means that if the target of flash flood disaster prevention is located outside the monitoring influence limit of the rain gauge, the coverage value is 0, reflecting that there is no monitoring correlation at a long distance.

[0069] Overall coverage The calculation formula is:

[0070]

[0071] In the formula, , The weights are respectively for river confluence-type disaster prevention targets and slope overflow-type disaster prevention targets. This represents the coverage of the i-th river confluence-type disaster prevention object. N represents the coverage of the i-th slope runoff type disaster prevention object. A N represents a subset of disaster prevention targets related to river confluence. B This represents a subset of disaster prevention targets for slope runoff, m A m is the number of river confluence-type disaster prevention targets. B This refers to the number of disaster prevention targets for slope runoff.

[0072] S5. Rain gauge network layout optimization: Optimize the network round by round according to the optimization model. Under the premise of meeting the preset comprehensive coverage target, minimize the number of rain gauges to be retained. At the same time, find the solution that maximizes the watershed coverage in the optimal round and output the optimal rain gauge layout scheme under the target coverage.

[0073] The objective function for optimizing the model is:

[0074]

[0075] In the formula, Z represents the comprehensive target value for the optimization of the rain gauge network, which is a comprehensive balancing indicator of the construction cost of the rain gauge network and the watershed coverage. The unit cost of the j-th rain gauge is expressed in ten thousand yuan, covering the entire process of construction, equipment, operation and maintenance, and data transmission, reflecting the cost heterogeneity of different stations. , indicating the retention status of the j-th rain gauge, where 1 indicates retention and 0 indicates removal; The economic costs of retaining rain gauges throughout their entire lifecycle were quantified. This indicates the overall coverage of the watershed; the higher the value, the more comprehensive the coverage. The value represents the coverage gap; the larger the gap, the larger this value. λ is a weighting coefficient, where λ > 0. It is used to balance the priority between rain gauge costs and coverage. It can be dynamically adjusted according to the watershed risk level. For example, a larger value is taken for high-risk areas to prioritize watershed coverage, while a smaller value is taken for low-risk areas to prioritize controlling the cost of station network construction.

[0076] The objective function ultimately aims to minimize the number of retained rain gauges while meeting preset coverage requirements, and simultaneously finds the solution that maximizes watershed coverage in the optimal round. Optimal resource allocation is achieved by balancing the two objectives through weighted coefficients.

[0077] The constraints of the objective function include:

[0078] Comprehensive coverage constraints (core constraints):

[0079]

[0080] Constraints for dynamic monitoring radius calculation:

[0081]

[0082] Deployment state binary constraints:

[0083]

[0084] In the formula, To achieve the preset coverage target, , The weights are assigned to river confluence-type disaster prevention targets and slope runoff-type disaster prevention targets, respectively, reflecting the importance of mountain flood disaster prevention targets at this node. This represents the coverage of the i-th river confluence-type disaster prevention object. N represents the coverage of the i-th slope runoff type disaster prevention object. A N represents a subset of disaster prevention targets related to river confluence. B This represents a subset of disaster prevention targets for slope runoff, m A m is the number of river confluence-type disaster prevention targets. B The number of flood-prone disaster prevention targets on slopes. This represents the dynamic monitoring radius of the k-th sector area at the j-th rain gauge station, in km. This represents the standard deviation of the slope in the k-th monitoring area at the j-th rain gauge station. Let $\frac{j}{k}$ be the standard deviation of the slope in the $k$-th monitoring area of ​​the $j$-th existing rain gauge station, where $k = 1, 2, ..., 8$, and the unit is °. This represents the retention status of the j-th existing rain gauge, where 1 = retained and 0 = removed. S represents the set of candidate rain gauge locations. R represents the dynamic monitoring radius R of the k-th sector of any rain gauge j. j,k Slope standard deviation T j,k All must satisfy the dynamic monitoring radius calculation constraint; for each rain gauge j in the candidate rain gauge set S, its deployment state It can only be 0 or 1.

[0085] The overall coverage constraint means that the overall coverage must meet the preset coverage target. The dynamic monitoring radius calculation constraint indicates that the dynamic radius is determined by the standard deviation of the slope of the sector region. The decision ensures that the radius calculation is correlated with the slope characteristics; the binary constraint of the deployment status indicates that each candidate location can only "deploy" or "not deploy" a rain gauge.

[0086] The optimization logic of the optimization model is as follows: starting from the initial set of all rain gauge stations, the model searches in cycles, gradually approaching the optimal solution of the objective function. In each cycle, one rain gauge station is removed from the set of all rain gauge stations, and the watershed coverage under the current condition is calculated. If a feasible solution that meets the target coverage is found in the current cycle, the evaluation of the remaining combinations in the current cycle is terminated, and the next cycle is performed. If no feasible solution is found, the current cycle continues. When a feasible solution cannot be found in a certain cycle, the model backtracks to the previous cycle and selects the scheme with the highest coverage as the final optimization result, thereby calculating the optimal solution of the objective function of the optimization model.

[0087] Example

[0088] like Figures 2-6 As shown, a method for optimizing rainfall stations for flash flood disasters based on grid coverage is described, with the following specific steps:

[0089] S1. Classification of Mountain Flood Disaster Prevention Targets and Upstream Raster Extraction: Taking the Shiwen River Basin in Shandong Province as an example, data on 197 mountain flood disaster prevention targets within the basin were collected. Based on the spatial relationship between the targets and the river channel, they were classified into 167 river confluence-type targets and 30 slope runoff-type targets. To optimize the overly dispersed distribution of mountain flood disaster prevention targets in the Shiwen River Basin, a K-means clustering algorithm was used. Based on the elbow rule, the targets were simplified into 7 river confluence-type targets requiring calculation of the entire upstream area and 5 slope runoff-type targets considering only the point location itself. The river confluence-type targets were coded as A1-A7, and the slope runoff-type targets as B1-B5. The DEM data of the Shiwen River Basin was imported into ArcGIS, and the flow direction raster was obtained through flow direction analysis in ArcGIS. Import flow raster data, pre-scan the entire flow raster, and construct a reverse index dictionary from the target raster to all source rasters. When processing any flash flood disaster prevention object, use a reverse index and breadth-first search algorithm to efficiently extract all upstream rasters for each flash flood disaster prevention object, and store the raster's row and column coordinates, geographic coordinates, and flash flood disaster prevention object type in the database;

[0090] S2. Dynamic Monitoring Radius and Coverage Area Calculation: Centered on each rain gauge, the surrounding monitoring area is divided into 8 fan-shaped sub-regions with a radius of 5 kilometers at 45°. Slope grids are obtained through ArcGIS surface analysis. Based on the slope grid data, the slope values ​​of all grids in each fan-shaped sub-region are read, and the slope standard deviation of the region is further calculated. After the calculation is completed, the slope standard deviation is stored in the database to quantify the terrain complexity in each direction. The slope standard deviation of each fan-shaped region is read, and the monitoring radius of each fan-shaped direction is obtained according to the formula.

[0091]

[0092]

[0093] In the formula, T k This represents the standard deviation of the slope of the k-th monitoring area at the rain gauge station, in degrees (°); t k R represents the standard deviation of the slope in the k-th monitoring area of ​​the rain gauge; k A represents the coverage radius of the k-th monitoring area of ​​the rain gauge, in km; R This indicates the coverage area of ​​the rain gauge monitoring area, in km². 2 Based on this, the total coverage area of ​​each rain gauge station was obtained by combining the monitoring radii of the eight sectors.

[0094] S3. Spatial Intersection Analysis: Based on the sector-shaped dynamic monitoring radius of each rain gauge obtained in S2, construct its dynamic monitoring range polygon; establish an R-tree spatial index for the upstream raster set of each flash flood disaster prevention object: group and store the upstream rasters of the watershed according to the flash flood disaster prevention object, establish a flash flood disaster prevention object-raster dictionary structure, and construct an independent R-tree spatial index for the upstream raster dataset of each flash flood disaster prevention object. Then, utilize the spatial pruning characteristics of the R-tree to quickly eliminate obviously irrelevant rasters; then perform precise geometric relationship calculations on the candidate rasters to quickly filter out candidate rasters located within the sector-shaped buffer bounding box; next, perform spatial inclusion judgment on the candidate set to determine whether the rasters are within the monitoring radius of the rain gauge; finally, store the data in the database.

[0095] S4. Calculation of coverage of mountain torrent disaster prevention objects and calculation of regional comprehensive coverage: Substitute each mountain torrent disaster prevention object into the coverage calculation formula to calculate the coverage of each object; to simplify the calculation of regional comprehensive coverage, substitute the cluster points of the mountain torrent disaster prevention objects that have undergone clustering into the coverage calculation formula to calculate the regional comprehensive coverage.

[0096] S5. Rain gauge network layout optimization: Taking optimization model parameters This reduces the weight of coverage in the objective function, allowing it to maximize the reduction of the number of rain gauges and maximize watershed coverage in the optimal cycle. , All values ​​are set to 1, ensuring consistency in importance across all flood prevention targets; a preset coverage target is established. The value is 0.8. The calculation is performed based on the optimization model. Starting from the initial set of all rain gauges, the search is performed iteratively round by round. In each round, one rain gauge is removed from the set of all rain gauges, and the watershed coverage under the current condition is calculated according to the coverage calculation formula.

[0097] As the number of iterations increases, the number of candidate points gradually decreases, and the watershed coverage shows a decreasing trend, with the objective function also decreasing accordingly. When the optimization runs for the 8th time, the watershed coverage falls below the preset target value for the first time, indicating that there is no feasible solution in the current round. Although the objective function value in this round reaches the minimum, the watershed coverage is lower than the preset coverage, so this value is discarded. The 9th round of optimization backtracks to the 7th round and iterates through all possible subsets in this round to find the optimal watershed coverage in this round, while the objective function value reaches the minimum. The final optimized scheme retains 5 rain gauges. Its watershed coverage... =0.82, which meets the requirements. Furthermore, the objective function value is minimized at this point, indicating that the optimal balance between cost and coverage performance has been achieved while satisfying the coverage constraint.

[0098] Although the above embodiments have been shown and described, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Any changes, modifications, substitutions and variations made to the above embodiments by those skilled in the art are within the protection scope of the present invention.

Claims

1. A method for optimizing rainfall stations for flash flood disasters based on grid coverage, characterized in that, Includes the following steps: S1. Extraction of upstream associated grid cells for flash flood disaster prevention objects: Create a direction mapping dictionary, use the reverse index of the direction mapping dictionary, and use a breadth-first search algorithm to extract the upstream grid cells of flash flood disaster prevention objects; S2. Calculation of dynamic monitoring radius and coverage area: Divide each mountain flood disaster prevention target into multiple fan-shaped sub-regions, calculate the slope standard deviation of each fan-shaped sub-region and the monitoring radius of each fan direction, and combine the monitoring radii of all fan-shaped sub-regions to obtain the total coverage area of ​​each rain gauge station; S3. Spatial Intersection Analysis: Based on the sector dynamic monitoring radius obtained in S2, construct its dynamic monitoring range polygon; create an R-tree spatial index to efficiently calculate the spatial intersection relationship between the dynamic monitoring range of each rain gauge station and the upstream area of ​​the flash flood disaster prevention object extracted in S1, and query and record the intersecting raster data and the corresponding type information; S4. Calculation of coverage of mountain torrent disaster prevention targets and calculation of regional comprehensive coverage: Based on the spatial intersection results obtained in S3, the coverage calculation formula considering the attenuation coefficient is used to calculate the coverage of the current rain gauge network for mountain torrent disaster prevention targets and the regional comprehensive coverage. S5. Rain gauge network layout optimization: Optimize the network round by round according to the optimization model. Under the premise of meeting the preset comprehensive coverage target, minimize the number of rain gauges to be retained. At the same time, find the solution that maximizes the watershed coverage in the optimal round and output the optimal rain gauge layout scheme under the target coverage.

2. The method for optimizing rainfall stations for flash flood disasters based on grid coverage according to claim 1, characterized in that, The specific steps for extracting upstream associated raster data for flood disaster prevention targets in S1 are as follows: A direction mapping dictionary is pre-created for the flow direction values ​​of the D8 algorithm, converting the flow direction values ​​of the D8 algorithm into corresponding row and column coordinate offsets; flow direction raster data is imported, the flow direction value of the central cell of the flood disaster prevention target is read, a reverse index is used with the direction mapping dictionary, and finally, a breadth-first search algorithm is used to extract all upstream raster data for each flood disaster prevention target, storing the row and column coordinates, geographic coordinates, and flood disaster prevention target type of the raster data in the database; this achieves the extraction of upstream catchment area raster data for any point source.

3. The method for optimizing rainfall stations for flash flood disasters based on grid coverage according to claim 2, characterized in that, The construction steps of the direction mapping dictionary reverse index are as follows: First, map the coordinates of the flood disaster prevention target points to the raster row and column numbers of the DEM to determine their positions in the raster matrix; create a direction mapping dictionary in advance for the flow direction values ​​of the D8 algorithm, converting the flow direction values ​​of the D8 algorithm into corresponding row and column coordinate offsets; when extracting upstream rasters, scan the flow direction rasters, obtain their flow direction values, and then convert the direction through the direction mapping dictionary, directly performing calculations on the row and column coordinates to calculate the position of the downstream cells; record the current cell in the reverse index of the downstream cells, repeat this operation until a mapping relationship of "downstream cells → all upstream cells" is established, thereby constructing a reverse index dictionary of "target raster to all source rasters".

4. The method for optimizing rainfall stations for flash flood disasters based on grid coverage according to claim 1, characterized in that, The calculation of the dynamic monitoring radius and coverage area of ​​S2 specifically refers to dividing the monitoring area around each mountain flood disaster prevention target into 8 fan-shaped sub-regions with a radius of 5 kilometers at 45°, taking each target as the center; based on the slope grid data, reading the slope values ​​of all grids in each fan-shaped sub-region, and further calculating the slope standard deviation of the region to quantify the terrain complexity in each direction; obtaining the monitoring radius of each fan-shaped direction based on the obtained slope standard deviation of each fan-shaped region; and then combining the monitoring radii of the 8 fan-shaped regions to obtain the total coverage area of ​​each rain gauge station.

5. The method for optimizing rainfall stations for flash flood disasters based on grid coverage according to claim 1, characterized in that, The spatial intersection analysis in S3 includes two stages: creating a spatial index and spatial querying. In the spatial index creation stage, an independent R-tree spatial index is established for the upstream raster set of each disaster prevention object, using the hierarchical structure of the R-tree: "root node → non-leaf node → leaf node". In the spatial querying stage, candidate rasters are first quickly queried in the R-tree using the bounding rectangle of the sector buffer. Then, the ray method is used to make a precise geometric relationship judgment on the candidate rasters to determine which rasters are actually located within the sector. Finally, all rasters within the sector are retained to form an "overlapping raster set".

6. The method for optimizing rainfall stations for flash flood disasters based on grid coverage according to claim 1, characterized in that, The classification method for mountain flood disaster prevention objects in S4 is as follows: an automatic identification and classification method for disaster prevention object types is adopted. The surface vector data and DEM raster data of the built-up area of ​​the mountain flood disaster prevention object are intersected to obtain the comprehensive raster data of the built-up area of ​​the mountain flood disaster prevention object. The catchment area and spatial location are calculated based on this data, and the mountain flood disaster prevention objects are divided into river confluence type disaster prevention objects and slope overflow type disaster prevention objects.

7. The method for optimizing rainfall stations for flash flood disasters based on grid coverage according to claim 6, characterized in that, The formula for calculating the coverage of the river confluence-type disaster prevention targets is as follows: In the formula, To ensure coverage of disaster prevention targets in river confluence areas, G represents the number of grid cells where the catchment area intersects with the coverage area of ​​the rain gauge, and G represents the number of grid cells upstream of the catchment area at or above the cross section of the flash flood disaster prevention target. The formula for calculating the coverage of the slope-flow type disaster prevention object is as follows: In the formula, denoted as the coverage of the slope flooding disaster prevention object, d is the straight-line distance from the slope flooding disaster prevention object to the nearest rain gauge station, r is the monitoring radius of the nearest rain gauge station of the slope flooding disaster prevention object, D is the maximum influence distance of the target rain gauge station on the slope flooding disaster prevention object, and k is the attenuation coefficient. Overall coverage The calculation formula is: In the formula, , The weights are respectively for river confluence-type disaster prevention targets and slope overflow-type disaster prevention targets. This represents the coverage of the i-th river confluence-type disaster prevention object. N represents the coverage of the i-th slope runoff type disaster prevention object. A N represents a subset of disaster prevention targets related to river confluence. B This represents a subset of disaster prevention targets for slope runoff, m A m is the number of river confluence-type disaster prevention targets. B This refers to the number of disaster prevention targets for slope runoff.

8. The method for optimizing rainfall stations for flash flood disasters based on grid coverage according to claim 1, characterized in that, The objective function of the optimization model in S5 is: In the formula, Z represents the comprehensive target value for rain gauge network optimization, which is a comprehensive balancing indicator between the construction cost of the rain gauge network and the watershed coverage. Let the unit cost of the j-th rain gauge be... , indicating the retention status of the j-th rain gauge, where 1 indicates retention and 0 indicates removal. The economic costs of retaining rain gauges throughout their entire lifecycle were quantified. Indicates the overall coverage of the watershed. The coverage gap is represented by λ, which is a weighting coefficient used to balance the priority of rain gauge cost and coverage.

9. The method for optimizing rainfall stations for flash flood disasters based on grid coverage according to claim 8, characterized in that, The constraints of the objective function include: Overall coverage constraints: Constraints for dynamic monitoring radius calculation: Deployment state binary constraints: In the formula, To achieve the preset coverage target, , The weights are respectively for river confluence-type disaster prevention targets and slope overflow-type disaster prevention targets. This represents the coverage of the i-th river confluence-type disaster prevention object. N represents the coverage of the i-th slope runoff type disaster prevention object. A N represents a subset of disaster prevention targets related to river confluence. B This represents a subset of disaster prevention targets for slope runoff, m A m is the number of river confluence-type disaster prevention targets. B The number of flood-prone disaster prevention targets on slopes. This represents the dynamic monitoring radius of the k-th sector area at the j-th rain gauge station. This represents the standard deviation of the slope in the k-th monitoring area at the j-th rain gauge station. Let be the standard deviation of the slope of the k-th monitoring area at the j-th existing rain gauge station, where k = 1, 2, ..., 8. This indicates the retention status of the j-th existing rain gauge, where 1 = retained and 0 = removed, and S represents the set of candidate rain gauge locations.

10. The method for optimizing rainfall stations for flash flood disasters based on grid coverage according to claim 9, characterized in that, The optimal solution of the optimization model is calculated as follows: starting from the initial set of all rain gauge stations, the search is performed cyclically round by round to gradually approach the optimal solution of the objective function. In each round, one rain gauge station is removed from the set of all rain gauge stations, and the watershed coverage under the current condition is calculated. Once a feasible solution that meets the target coverage is found in the current round, the evaluation of the remaining combinations in the current round is terminated, and the next round of calculation is performed. If no feasible solution is found, the current round continues. When a feasible solution cannot be found in a certain round, the previous round is backtracked and the scheme with the highest coverage is selected as the final optimization result, thereby calculating the optimal solution of the objective function of the optimization model.

Images