A candidate depression comprehensive evaluation method based on water immersion simulation and double-layer decision model

Abstract

The present application relates to the technical field of topographic form analysis, and provides a candidate depression comprehensive evaluation method based on water immersion simulation and double-layer decision model, wherein a spatial form quantization index is constructed based on a watershed boundary, and quantitative expression of the form characteristics of the candidate depression is realized; further, a standard pole depression is taken as a reference to construct a form similarity and a squareness optimization double criterion, and a comprehensive ranking of multiple candidate depressions is realized through objective weighting of information entropy, so that rapid and quantitative selection of the candidate depression is realized, interference of a traditional analytic hierarchy process and other subjective weighting methods on the screening result is avoided, and the objectivity, reproducibility and engineering applicability of large-scale candidate depression screening are improved.

Description

Technical Field

[0001] This invention relates to the field of terrain morphology analysis technology, and in particular to a comprehensive evaluation method for candidate depressions based on water immersion simulation and a two-layer decision model. Background Technology

[0002] Karst depressions, with their natural electromagnetic shielding effect, relatively enclosed topographic environment, and smaller spherical fitting cut and fill volume, have become the preferred landform type for selecting candidate sites for large radio telescopes. However, karst depressions are widely distributed, numerous, and diverse in shape. To select the optimal candidate sites for radio telescopes, it is necessary to fully analyze the spatial morphology of the candidate depressions and formulate a scientific and reasonable selection mechanism.

[0003] In terms of evaluating and screening candidate depressions, existing technologies are mainly based on scoring mechanisms such as the analytic hierarchy process (AHP) and rely on technical personnel to verify the geological conditions of each candidate depression one by one, including water system distribution, transportation accessibility, statistics of the number of surrounding peaks, and interpretation of the location of mountain passes. The influence of human subjective judgment is strong, which can easily lead to unreasonable screening results. For large-scale depression screening scenarios, there are problems of large workload and difficulty in unifying evaluation standards.

[0004] Therefore, it is necessary to propose a comprehensive evaluation method for candidate depressions based on immersion simulation and a two-layer decision-making model to address or at least mitigate the aforementioned shortcomings. Summary of the Invention

[0005] The main objective of this invention is to provide a comprehensive evaluation method for candidate depressions based on immersion simulation and a two-layer decision-making model, in order to solve the technical problems of existing technologies having a strong influence of human subjective judgment, which can easily lead to unreasonable screening results; and having a large workload and difficulty in unifying evaluation standards for screening large-scale depressions.

[0006] To achieve the above objectives, this invention provides a comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision model, comprising the following steps: S1, Obtain the spatial morphological index dataset of the benchmark depression and multiple candidate depressions. The spatial morphological index dataset includes multiple spatial morphological quantification indices and volume indices generated based on the watershed boundaries and digital elevation model data of each candidate depression. S2, perform collinearity removal and principal component dimensionality reduction on the spatial morphology quantification index in the spatial morphology index dataset to obtain multiple principal components after dimensionality reduction, and obtain the principal component score vectors of the benchmark depression and each candidate depression based on the multiple principal components; S3, Construct a two-layer decision-making model, which includes a morphological distance layer and an economic attribution layer; S4, in the morphological distance layer, taking the benchmark depression as a reference, calculate the Mahalanobis distance of each candidate depression to the benchmark depression in the principal component space composed of multiple principal components, and determine the morphological similarity of each candidate depression based on the Mahalanobis distance, and use the morphological similarity as criterion A score. S5, in the economic attribution layer, the principal component scores of each candidate depression on the multiple principal components are used as independent variables, and the volume index of each candidate depression is used as dependent variable. A partial least squares regression mapping is established. Based on the partial least squares regression mapping, the predicted volume of each candidate depression is output, and the volume goodness of each candidate depression is determined based on the predicted volume. The volume goodness is used as the criterion B score. S6, normalize the scores of criterion A and criterion B, and use information entropy weighting to determine the criterion weights of criterion A and criterion B. S7. Based on the normalized score of criterion A, the normalized score of criterion B, the criterion weight of criterion A, and the criterion weight of criterion B, the comprehensive score of each candidate depression is obtained, and the multiple candidate depressions are sorted according to the comprehensive score, and the comprehensive evaluation result of the candidate depressions is output.

[0007] Preferably, the spatial morphology quantification index includes the closure coefficient, number of peaks, number of passes, maximum depth, net depth, length of the major axis of the depression, length of the minor axis of the depression, diameter of the great circle of the watershed, diameter of the great circle of the pass, roundness coefficient, slope, directionality index, and structural similarity index.

[0008] Preferably, step S2 involves collinearity removal and principal component dimensionality reduction of the spatial morphology quantification indices in the spatial morphology index dataset, including the following steps: S21, Calculate the variance inflation coefficient of the spatial morphology quantification index in the spatial morphology index dataset; S22, remove redundant indicators with a variance inflation coefficient greater than 10; S23, perform principal component analysis on the retained spatial morphology quantification index and extract principal components with a cumulative variance contribution rate greater than 90%.

[0009] Preferably, step S2 further includes the following step: The principal component loading matrix is ​​rotated to maximize variance to improve the interpretability of the principal components, and the dominant indices of each principal component are determined based on the rotated principal component loading matrix.

[0010] Preferably, in step S4, determining the morphological similarity of each candidate depression based on the Mahalanobis distance includes the following steps: S41, For each candidate depression, calculate the difference vector between the principal component score vector of the candidate depression and the principal component score vector of the benchmark depression; S42, calculate the Mahalanobis distance from the candidate depression to the benchmark depression in the principal component space based on the difference vector, the sample covariance matrix, and the inverse matrix of the sample covariance matrix; S43, the Mahalanobis distance is converted into morphological similarity through inverse transformation; Among them, the smaller the Mahalanobis distance between the candidate depressions, the higher the morphological similarity of the candidate depressions.

[0011] Preferably, the partial least squares regression mapping established in S5 includes the following steps: S51, the principal component scores of each candidate depression on the multiple principal components are used as the independent variables of the partial least squares regression mapping; S52, the volume index of each candidate depression is used as the dependent variable of the partial least squares regression mapping; S53, the optimal number of components in the partial least squares regression mapping is determined by leave-one-out cross-validation; S54. Based on the optimal number of components, establish a partial least squares regression mapping that predicts the volume index by the principal component scores of each candidate depression.

[0012] Preferably, the step S5, which determines the quantitative superiority of each candidate depression based on the predicted quantitative volume, includes the following steps: S55, perform inverse standardization on the predicted volume to obtain the volume quality of each candidate depression; The range of the quantity goodness is [0,1]. The closer the quantity goodness of a candidate depression is to 1, the smaller the predicted quantity of the candidate depression and the better its economy.

[0013] Preferably, step S5 further includes the following step: The regression residuals are output based on the partial least squares regression mapping, and the candidate depressions are marked as abnormal samples based on the regression residuals.

[0014] Preferably, step S5 further includes the following step: Calculate the variable projection importance index for each principal component, and determine the importance score of each spatial morphology quantification index for volume prediction based on the variable projection importance index for each principal component and the absolute value of the rotational load of each spatial morphology quantification index on each principal component.

[0015] Preferably, the determination of the criterion weights of criterion A and criterion B using information entropy weighting in step S6 includes the following steps: S61, normalize the scores of criterion A and criterion B for each candidate depression to the interval [0,1] respectively; S62, calculate the information entropy values ​​of criteria A and criteria B respectively; S63, Determine the discriminative power of criterion A based on the information entropy value of criterion A, and determine the discriminative power of criterion B based on the information entropy value of criterion B; S64. Determine the criterion weights of criterion A and criterion B based on the discriminative power of criterion A and criterion B. The sum of the criterion weights of criterion A and criterion weights of criterion B is 1. The overall score is obtained by weighted summation of the normalized scores of criterion A and criterion B.

[0016] Compared with the prior art, the present invention has the following beneficial effects: This invention constructs a spatial morphology quantification index based on watershed boundaries, realizing the quantitative expression of the morphological characteristics of candidate depressions. Furthermore, it proposes to construct a dual criterion of morphological similarity and volume superiority based on benchmark depressions, and to achieve comprehensive ranking of multiple candidate depressions through objective weighting by information entropy. This enables rapid and quantitative selection of candidate depressions, avoids interference from subjective weighting methods such as traditional analytic hierarchy process, and improves the objectivity, reproducibility, and engineering applicability of large-scale candidate depression selection. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the structures shown in these drawings without creative effort.

[0018] Figure 1 This is a schematic flowchart of one embodiment of the present invention.

[0019] Figure 2 This is a schematic diagram of the watershed boundary of a candidate depression obtained based on immersion simulation in one embodiment of the present invention.

[0020] The objectives, features, and advantages of this invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0021] It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0022] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0023] Furthermore, the technical solutions of the various embodiments can be combined with each other, but only if they are feasible for those skilled in the art. If the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such combination of technical solutions does not exist and is not within the scope of protection claimed by this invention.

[0024] It is worth noting that in the extraction of watershed boundaries in depressions, existing technologies mainly rely on traditional hydrological analysis methods to segment the entire digital elevation model (DEM) to obtain the watershed ranges of different catchment units. This method, based on the principle of "regional competition," has two main drawbacks: First, the computational load of the entire region segmentation is large, while the area of ​​the depressions that ultimately become candidates is often less than 1%, resulting in a large amount of computational resources being consumed in irrelevant areas and excessive computational costs. Second, the entire region segmentation is prone to oversegmentation, requiring manual operations such as merging regions and defining the ranges of watershed boundaries for candidate depressions in the later stages, increasing the workload of duplication.

[0025] Please refer to Figures 1 to 2 The present invention provides a method for extracting and morphologically quantifying candidate depression watersheds based on flooding simulation, comprising the following steps: S1. Obtain the digital elevation model (DEM) data corresponding to the candidate depressions and determine the center pixel in the DEM data corresponding to the candidate depressions. The candidate depressions are those initially delineated through field surveys, remote sensing interpretation, etc. The center pixel is the pixel in the DEM data corresponding to the center point of the candidate depression. The center of the candidate depression is usually located near the lowest point of the depression. Water injection simulation begins with the center pixel, with water spreading outwards from the deepest point of the depression to ensure the initial calculation area is located inside the depression and naturally connected to its catchment area. If starting from a slope or pass, the initial area may extend beyond the depression's boundaries, causing subsequent iterations of the inundation range to fail to accurately reflect the depression's catchment process.

[0026] S2, based on the elevation value of the central pixel and the elevation threshold used to define the initial calculation area of ​​the depression, extract the set of pixels connected to the central pixel to obtain the initial calculation area of ​​the depression, and use this initial calculation area as the currently determined depression area; the initial calculation area of ​​the depression is located inside the depression and does not touch the mountain pass boundary. Extracting the set of pixels connected to the central pixel is to ensure that the initial calculation area is a single continuous pool of water and is within the bottom range of the candidate depression. The connectivity condition excludes two types of irrelevant pixels: first, pixels that are separated from the bottom of the depression by elevation barriers and cannot be naturally reached by water flow; second, pixels that belong to adjacent independent depressions and are separated by ridges. Without the connectivity condition, the initial area may cross the ridge and include the bottom of adjacent depressions, and the flooding range of subsequent iterations will include the water bodies of two depressions, and the area growth rate and perimeter growth rate will not accurately reflect the overflow event of a single depression.

[0027] S3 updates the water level value according to the preset water level rise step size, extracts the set of pixels whose elevation is less than or equal to the updated water level value and is connected to the currently identified depression area, and obtains the candidate inundation area; in each iteration, pixels that meet the elevation condition after the water level rise may form multiple unconnected areas. Only the area connected to the currently identified depression area is the inundation range after the water level rises in this depression. The connectivity condition simulates the physical constraint that water flow is blocked by ridges and cannot cross watersheds to enter adjacent depressions.

[0028] S4. Based on the area growth rate and perimeter growth rate of the candidate inundation area relative to the currently identified depression area, determine whether the candidate inundation area has overflowed. During a normal water level rise, water spreads evenly along the slope, and the growth rates of area and perimeter are stable and small. When the water level reaches a pass and water flows through the pass into an adjacent depression, the inundation area suddenly incorporates a large area of ​​the adjacent depression, resulting in a step increase in area or perimeter. The area growth rate is sensitive to large-scale breaches, while the perimeter growth rate is sensitive to overflows through narrow channels. Overflow is determined if either indicator exceeds a preset threshold, avoiding the risk of missed detection that may occur with a single indicator.

[0029] S5, if overflow occurs in the candidate flooded area, extract the adjacent boundary pixels between the currently determined depression area and the candidate flooded area. Based on the overflow judgment results before and after excluding the adjacent boundary pixels from the candidate flooded area, determine the overflow pixels from the adjacent boundary pixels and exclude the overflow pixels in the water level iteration. The boundary pixel set is the boundary zone between the currently determined depression area and the candidate flooded area, that is, a ring of pixels that has just been submerged after the water level rises. The purpose of extraction is to narrow the search range of the overflow point. When overflow occurs, the water flows through this boundary zone into the adjacent depression. Narrowing the search range from the entire candidate flooded area to this ring of boundary pixels greatly reduces the computational load of subsequent verification. Each pixel on this ring of boundary pixels may be a normal slope (where the water flow is just normally flooding) or a pass (where the water flow is flowing out from the pass).

[0030] For each cell in the boundary cell set, it is temporarily released (set as unblocked), and the flooding simulation and overflow determination are re-executed to obtain the overflow determination result after release. Before this, when the cell was in a blocked state, the overflow determination result of the flooding simulation was either overflow or not overflow. Comparing the two overflow determination results before and after release, if the cell did not overflow before release due to other cells being blocked, but overflowed after release because the cell was released, then the cell is determined to be an overflow cell.

[0031] S6: If the candidate flooded area does not overflow, update the candidate flooded area to the currently determined depression area; S6 and S5 are parallel branches, and only one of them is executed in each iteration. When no overflow occurs, it indicates that the water flow is only expanding normally within the depression, and the new flooded area is confirmed as the reference area for the next iteration.

[0032] S7. Repeat steps S3 to S6 until the water level reaches the preset termination elevation. Obtain the watershed boundary of the candidate depression based on the overflow pixels. Then, based on the watershed boundary and the digital elevation model data, generate a quantitative result of the candidate depression's spatial morphology. Based on the watershed boundary and the digital elevation model data, quantitative morphological indicators such as the number of peaks, the number of passes, the lowest pass elevation, the maximum depth, the net depth, the major axis length, the minor axis length, the average slope, the roundness, and the cut / fill volume can be calculated. Further explanation will follow.

[0033] like Figure 2 As shown, Figure 2 The closed curve in the middle region marks the boundary line of the watershed.

[0034] This invention simulates the connectivity expansion process of the flooded area by raising the water level layer by layer, uses the area growth rate and perimeter growth rate as dual indicators to detect the overflow phenomenon, and combines the release of adjacent boundary pixels one by one to verify that it can distinguish between real mountain passes and pseudo-slope boundaries. This makes the extracted watershed boundaries distributed along the ridgeline and highly consistent with the actual water catchment units, avoiding the over-segmentation phenomenon in traditional depression watershed extraction methods.

[0035] This invention automates the entire process, from acquiring digital elevation model data of depressions, determining the initial calculation area of ​​depressions, water level iteration, overflow detection and overflow pixel location, to generating watershed boundaries and outputting spatial morphology quantification results. This reduces reliance on human interpretation experience and only analyzes elevation data around the target depression area. The algorithm's timeliness is exponentially improved compared to traditional methods, making it suitable for rapid processing of large-scale, multi-candidate depressions.

[0036] This invention adopts an elevation-layered iterative approach, where each water level iteration only requires connected domain marking calculations. Compared to the traditional multi-step iterative process involving filling depressions, calculating flow direction, and accumulating confluence, this reduces redundant calculation steps. At the same time, by setting an initial calculation area, it avoids initial overflow misjudgments and invalid iterations caused by starting from a single pixel, thereby improving the stability of overflow judgment and overall calculation efficiency.

[0037] In a preferred embodiment, S1 includes the following steps: Obtain the center coordinates of the candidate depression; wherein, the candidate depression is a depression that simultaneously meets the following conditions: the number of peaks on the enclosing ridge line surrounding the candidate depression is greater than 3, the height difference between the lowest point of the candidate depression and the enclosing ridge line is greater than 100 meters, and the ratio of the major axis length to the minor axis length of the area enclosed by the enclosing ridge line is less than a preset ratio threshold. Specifically, the candidate depressions are karst depressions initially screened based on field surveys, remote sensing satellite imagery, and topographic data, with their center coordinates accurate to the meter level. This step further filters the candidate depressions based on three criteria to eliminate those with lower morphological suitability.

[0038] The first condition is that the number of peaks on the enclosing ridgeline surrounding the candidate depression is greater than three. The enclosing ridgeline refers to the continuous ridge line that surrounds the depression, forming the boundary between the depression and external water flows. This ridgeline can be identified through remote sensing image interpretation, topographic map reading, or topographic data analysis. The peaks refer to pixels on the enclosing ridgeline whose elevation values ​​are local maxima. The number of peaks is obtained by counting the number of pixels on the enclosing ridgeline that satisfy the local maxima condition. In this embodiment, it is preferred to set the number to greater than three to screen candidate depressions with better enclosure. Generally, the more peaks there are, the better the enclosure and sealing of the depression, which is more conducive to forming the mountain shielding conditions required for radio telescopes. Limiting the number of peaks to greater than three ensures that the candidate depression has basic enclosure and excludes open terrain with sparse surrounding peaks and large openings.

[0039] The second condition is that the elevation difference between the lowest point of the candidate depression and the enclosing ridgeline is greater than 100 meters. Preferably, the elevation difference can be obtained by subtracting the elevation value of the lowest pixel within the depression from the average elevation value of each pixel on the enclosing ridgeline. Using the average elevation of the enclosing ridgeline as the representative elevation can comprehensively reflect the overall depth level of the depression and reduce the influence of individual extremely high or low points on the depth evaluation.

[0040] The third condition is that the ratio of the major axis length to the minor axis length of the area enclosed by the ridgeline is less than a preset ratio threshold. In this embodiment, the major axis length refers to the maximum length among all cross-sections passing through the lowest point of the depression within the polygonal area enclosed by the ridgeline; the minor axis length refers to the maximum length among cross-sections perpendicular to the major axis direction. The preset ratio threshold can be set according to actual site selection requirements. The closer the ratio is to 1, the closer the depression shape is to a circle. An approximately circular depression is beneficial for the geometric fitting of the radio telescope's spherical surface, reducing the amount of excavation and filling during construction and lowering project construction costs.

[0041] Based on the center coordinates, obtain digital elevation model data covering the candidate depressions; After determining the center coordinates of the candidate depression, a Digital Elevation Model (DEM) covering the candidate depression is obtained based on these coordinates. The DEM data is in raster format, with each cell containing an elevation value representing the altitude of the corresponding geographic location. The acquired DEM data should cover at least the entire area within a preset radius centered on the candidate depression to minimize the truncation of complete topographic information required for subsequent watershed extraction. Preferably, the preset radius is 1 kilometer. This value is based on the aperture of a large radio telescope (approximately 500 meters) and considers the topographic analysis buffer zone. With a 500-meter aperture, a 1-kilometer radius covers an area with a diameter of 2 kilometers, four times the telescope aperture, which can completely cover the watershed catchment unit of the candidate depression while avoiding the inclusion of too much irrelevant topographic data and unnecessary computational overhead. Preferably, the DEM data uses the ASTER GDEM V3 dataset, which has a spatial resolution of 30 meters, meeting the accuracy requirements for karst depression morphology analysis.

[0042] The digital elevation model data is subjected to coordinate transformation and spatial filtering. Coordinate transformation is used to unify DEM data to a single geographic coordinate system. Since candidate depression coordinates may originate from different data sources, to avoid discrepancies in subsequent calculations of geometric parameters such as area and perimeter due to inconsistent coordinate systems, all data needs to be transformed to the same coordinate system. Preferably, the coordinate transformation is based on the CGCS2000 geographic coordinate system.

[0043] Spatial filtering is used to reduce the impact of local noise in digital elevation model (DEM) data on subsequent connected component growth. During the acquisition and processing of DEM data, there may be random errors such as abnormally high or low elevation values ​​in individual pixels. These error pixels may form false barriers or false channels during subsequent connected component growth, interfering with the accuracy of watershed extraction. This step uses Gaussian filtering to convolve the elevation matrix of the DEM data. After Gaussian filtering, random spike noise in the DEM data is effectively suppressed, while the true depression morphology features with diameters of tens of meters or more are preserved, thereby improving the accuracy of subsequent watershed location.

[0044] The center coordinates are mapped onto the processed digital elevation model data to obtain the center pixel.

[0045] In a preferred embodiment, S2 includes the following steps: The elevation value of the center pixel is added to the elevation threshold to obtain the initial upper elevation limit; The elevation value of the center pixel can be obtained from the processed digital elevation model data, representing the altitude corresponding to the center point of the candidate depression. An elevation threshold is used to extend the elevation range upwards from the bottom of the depression to delineate an initial calculation area of ​​a certain size. Preferably, the elevation threshold is 10 meters. This value can be based on statistical observations of the morphological characteristics of the candidate depression: generally, the elevation at the lowest pass of the candidate depression is more than 20 meters higher than the lowest point. Taking 10 meters as the initial extension range helps to ensure that the initial calculation area of ​​the depression is located within the candidate depression and reduces the possibility that the initial calculation area of ​​the depression exceeds the lowest pass.

[0046] The initial elevation upper limit defines the elevation screening criteria for the initial calculation area. Pixels with elevations less than or equal to the initial elevation upper limit can be used as candidate pixels for the initial calculation area of ​​depressions, thus providing an initial calculation basis for subsequent iterations.

[0047] Eight-neighbor connectivity refers to checking eight adjacent pixels in the pixel grid of digital elevation model data, centered on the pixel being processed: top, bottom, left, right, upper left, upper right, lower left, and lower right. Connectivity determination requires two conditions to be met simultaneously: first, the pixel to be detected is adjacent to a pixel already included in the region in the eight-neighbor direction; second, the elevation value of the pixel to be detected is less than or equal to the initial upper elevation limit. This implementation uses the central pixel as the starting point for growth and employs a region growing algorithm (such as the flood fill method), which is a mature technology and will not be elaborated upon. The resulting set of pixels represents the largest connected region under the constraint of the initial upper elevation limit.

[0048] It is worth noting that the eight-neighborhood approach, by simultaneously detecting four diagonal directions, can capture continuous terrain extending along the diagonal, avoiding false connectivity breaks in gentle slope areas and more accurately simulating the physical process of water flowing naturally across the surface in all directions. If a cell is not connected, it means that although its elevation meets the initial upper limit of elevation, there is a cell with an elevation exceeding the upper limit blocking it from the center cell. Water cannot naturally flow from the bottom of the depression to that location, and therefore it cannot be included in the initial range of this depression; or the cell itself is located in another independent depression. Excluding unconnected cells ensures that the initial calculation area is limited to the bottom range of a single candidate depression, providing a correct initial benchmark for the area growth rate and perimeter growth rate to accurately reflect the occurrence of overflow events in subsequent iterations.

[0049] It's also worth noting that starting the calculation with the depression has two advantages: first, it reduces the number of iterations and algorithm time; second, if the water level increase starts directly from the center pixel of the depression, the initial area growth rate will be relatively large, easily leading to misjudgments of overflow. Using the core area of ​​the depression with a certain area as the initial value results in more robust area and perimeter growth rates, leading to higher accuracy in identifying overflow phenomena.

[0050] Calculate the area and perimeter of the initial calculation area of ​​the depression; Specifically, the area is calculated as follows: The total number of pixels within the initial calculation area is multiplied by the actual ground area represented by each pixel. The actual ground area of ​​a pixel is equal to the square of the spatial resolution of the digital elevation model data. The perimeter is calculated as follows: First, the boundary pixels of the initial calculation area are extracted. It should be noted that a boundary pixel is a pixel within the initial calculation area whose eight neighboring pixels do not belong to the initial calculation area. The total number of boundary pixels is counted and multiplied by the side length of a single pixel (i.e., spatial resolution, such as 30 meters) to obtain the perimeter of the initial calculation area.

[0051] The initial calculation area of ​​the depression is taken as the currently determined depression area during the first iteration.

[0052] The initial calculation region of the depression is assigned to the currently determined depression region as the initial state value for subsequent iterations from S3 to S6. During subsequent iterations, the currently determined depression region will be dynamically updated as the water level rises layer by layer. When a candidate flooding area does not overflow, it is updated to the currently determined depression region; when overflow occurs, the area redefined after processing through the overflow channel is used as the currently determined depression region. This initialization operation provides the basis for calculating the area growth rate and perimeter growth rate in each subsequent iteration.

[0053] In a preferred embodiment, S3 includes the following steps: Construct a blocking array of the same size as the digital elevation model data, and set the initial value of all elements in the blocking array to an unblocked state; Digital elevation model (DEM) data is stored in the computer as a two-dimensional raster matrix with a defined number of rows and columns. The clogging array is a two-dimensional array with the same number of rows and columns as the DEM data, with each element mapping one-to-one with a corresponding cell in the DEM data. The elements of the clogging array have two states: unclogged and clogged. All elements of the clogging array are initialized to the unclogged state because the flooding simulation has not yet begun; no water level rise or overflow detection has been performed, and no overflow points have been identified. Water can freely submerge any cell that meets the elevation conditions. As overflow points are identified in subsequent iterations, the corresponding elements in the clogging array are gradually updated to the clogging state. Cells in the clogging state cannot be submerged by water at any water level. Cells marked as clogging are no longer included in the flooding candidate cell set in subsequent iterations at higher water levels, thus simulating the effect of water being blocked by a dam and unable to overflow from the identified breach. When all mountain passes are blocked and the water level exceeds the highest elevation, all the blocked pixels form a closed ridgeline, which constitutes the watershed boundary.

[0054] In each iteration, the current water level value is updated according to the preset water level rise step size, and pixels that meet the following three conditions are selected from the digital elevation model data to obtain a set of flood candidate pixels: the elevation value of the pixel is less than or equal to the current water level value, the elevation value of the pixel is not empty, and the corresponding value of the pixel in the blocking array is in an unblocked state. The preset water level rise step size refers to the magnitude of the water level rise in each iteration. A smaller preset water level rise step size results in higher accuracy in extracting the watershed boundary, but also increases the number of iterations and computational cost. Preferably, the preset water level rise step size is 1 meter. This value is determined based on engineering experience with 30-meter resolution DEM data and typical morphological features of karst depressions. Those skilled in the art can also make adaptive settings according to actual needs. It is worth noting that in the first iteration, the water level value from the previous iteration is the upper limit of the initial elevation determined in step S2.

[0055] After determining the current water level, iterate through each cell in the DEM data and check whether each cell simultaneously meets the three screening criteria: Condition 1: The pixel's elevation value is less than or equal to the current water level. Condition 1 simulates the terrain area submerged by water after the water level rises; only locations with elevations lower than or equal to the water surface will be covered by water.

[0056] Condition 2: The elevation value of the pixel is not null. In DEM data, measurement data may be missing at certain locations (e.g., areas covered by water, data blind spots, etc.), and the pixels at these locations are marked as null. Null pixels do not contain valid elevation information and cannot participate in elevation comparison and connectivity analysis, so they are excluded. This ensures that each selected pixel has real and reliable elevation data, avoiding calculation anomalies caused by missing data.

[0057] Condition 3: The corresponding value of the cell in the blocking array is in an unblocked state. If a cell was identified as an overflow cell in a previous iteration and marked as blocked in the blocking array, then Condition 3 is no longer satisfied, and this cell will not be included in the flood candidate cell set in this iteration and all subsequent iterations at higher water levels. Condition 3 simulates the physical effect of damming a pass. Even if the water level continues to rise and exceeds the elevation of the cell, the water flow cannot pass through that location, and the blocked pass remains effective in subsequent water level iterations.

[0058] The set of candidate submerged pixels is labeled with eight neighboring connected regions, and the connected regions that are connected to the currently determined depression area are extracted as candidate submerged areas of the current water level layer.

[0059] Eight-neighbor connectivity labeling refers to grouping spatially connected cells in the flooded candidate cell set into a single connected region. Eight-neighbor connectivity labeling is a mature algorithm; this implementation preferably uses a connectivity labeling function based on the SciPy library. This function is compiled and implemented in C, offering high execution efficiency and enabling connectivity labeling of the entire DEM grid in a single function call, saving computation time for each iteration. After connectivity labeling, the flooded candidate cell set is divided into several independent, unconnected connected regions. In karst topography, the current water level may simultaneously flood multiple spatially unconnected low-lying terrains, such as pools within the current candidate depression, pools within adjacent depressions, and independent small pits, separated by ridges or highlands in the DEM discrete grid space, and thus unconnected.

[0060] It is worth noting that among the several independent connected regions obtained by the connectivity domain marking, not all connected regions belong to the catchment area of ​​this candidate depression. Although other connected regions are also below the current water level and are not blocked, they are separated from the currently identified depression area by ridges or highlands and are not spatially connected to each other. They belong to the catchment area of ​​adjacent depressions or other independent low-lying terrains and should not be included in the candidate inundation area of ​​this candidate depression.

[0061] In a preferred embodiment, S4 includes the following steps: Calculate the area and perimeter of the candidate flooded area, and obtain the area growth rate based on the ratio of the area of ​​the candidate flooded area to the area of ​​the currently determined depression area; The perimeter growth rate is obtained by the ratio of the perimeter of the candidate flooded area to the perimeter of the currently determined depression area; If the area growth rate is greater than or equal to a preset area growth rate threshold, or the perimeter growth rate is greater than or equal to a preset perimeter growth rate threshold, it is determined that the candidate flooded area has overflowed.

[0062] It is worth noting that when the water level rises normally along the slope inside the depression, both the area and perimeter of the inundated area increase at a relatively stable rate, remaining below their respective preset thresholds. However, when the water level reaches the pass height and overflow occurs, the water flows into adjacent low-lying areas through the narrow pass channel. At this point, the area growth may not be fully reflected due to the narrow channel, but the perimeter increases significantly due to the addition of a narrow overflow channel. If overflow determination is based solely on the area growth rate, overflow events with narrow channels and lagging area growth may be missed. By simultaneously monitoring both the area growth rate and the perimeter growth rate, the reliability and accuracy of overflow detection are significantly improved.

[0063] The preset area growth rate threshold and preset perimeter growth rate threshold can be determined empirically based on the DEM data resolution, preset water level rise step size, and typical morphological characteristics of karst depressions, with a value range generally between 1.1 and 1.5. Values ​​that are too low may misjudge normal slope expansion as overflow, while values ​​that are too high may miss overflow events at narrow passes. In this embodiment, both the preset area growth rate threshold and the preset perimeter growth rate threshold are preferably set to 1.2.

[0064] It is important to note that if the water level has not yet reached any pass height, the candidate flooded area mainly corresponds to the flooded range within the current candidate depression. The area growth rate and perimeter growth rate usually remain below the corresponding preset thresholds, and overflow determination is not triggered. If the current water level reaches the elevation near the pass, the connected region corresponding to the current candidate depression in step S3 may be connected to the connected regions corresponding to adjacent depressions through the pass pixels, and both may be extracted as candidate flooded areas. At this time, because the candidate flooded area contains some pixels of adjacent depressions, its area and boundary perimeter may be significantly larger than the currently determined depression area. Based on the comparison results of the area growth rate and perimeter growth rate with the corresponding preset thresholds, it is determined whether the candidate flooded area has overflowed. If overflow is determined, step S5 is entered to extract adjacent boundary pixels and determine overflow pixels; if overflow is determined, step S6 is entered to update the candidate flooded area to the currently determined depression area, and the next iteration begins.

[0065] In a preferred embodiment, S5 includes the following steps: Morphological dilation is performed on the currently identified depression area, and the intersection of the dilated currently identified depression area and the candidate flooding area is obtained to obtain the boundary cell set; The purpose of the dilation operation is to create a spatial intersection between the currently identified depression area and the candidate flooding area. In flooding simulation, the candidate flooding area is formed by expanding outward from the currently identified depression area. Therefore, the candidate flooding area contains all the pixels of the currently identified depression area. By performing morphological dilation on the currently identified depression area, its boundary is expanded outward by one pixel width. The expanded area and the candidate flooding area create a ring-shaped intersection zone with a width of approximately one pixel. This zone represents the leading edge of the water flow just as it overflows the old boundary and expands outward, and it is also the channel area through which the water flow passes when the overflow occurs.

[0066] The expanded, currently identified depression area is intersected with the candidate flooding area on a pixel-by-pixel basis. The pixels shared by both areas form the boundary pixel set, which significantly reduces the search range of the overflow point from the entire candidate flooding area to a narrow annular strip area. This significantly reduces the number of pixels that need to be processed in subsequent verification steps and improves the overall computational efficiency of the algorithm.

[0067] All boundary pixels in the boundary pixel set are marked as blocked in the blocking array; traverse each boundary pixel and perform the following operations on the currently detected boundary pixel: release only one pixel at a time for testing. If overflow occurs, the pixel can be uniquely determined as the overflow channel point because all other pixels are blocked and water can only flow out from this opening.

[0068] The currently detected boundary cell is temporarily marked as unblocked in the blockage array, while the blockage status of the remaining boundary cells remains unchanged. This step simulates opening a small hole at the currently detected boundary cell to allow water to flow through, while ensuring that water cannot flow out from other locations.

[0069] Based on the current water level, a new set of candidate flooding pixels is generated, and connected regions connected to the currently identified depression area are extracted. It is then determined whether overflow has occurred. Based on the current water level, step S3 (candidate flooding pixel screening, connected region marking, and candidate flooding area extraction) and step S4 (overflow determination) are re-executed. The currently detected pixel, because it is temporarily released, meets condition three of step S3 and can participate in the formation of the flooded area; other boundary pixels remain blocked and cannot participate. This design ensures that each verification only evaluates the overflow capacity of the current single pixel.

[0070] If an overflow is determined to have occurred, the currently detected boundary cell is identified as an overflow cell and marked as blocked in the blocking array; if no overflow is determined to have occurred, the currently detected boundary cell is marked as unblocked in the blocking array.

[0071] If overflow is detected, the currently detected boundary cell is identified as an overflow cell and marked as blocked in the blocking array, so that the overflow cell is excluded from the composition of the flood candidate cell set in subsequent water level iterations. If overflow is not detected, the currently detected boundary cell is marked as unblocked in the blocking array, so that the currently detected boundary cell can normally participate in the selection of the flood candidate cell set and the expansion of the connected region in subsequent water level iterations.

[0072] Through the above processing, the boundary cells identified as overflow cells remain blocked in the blocking array, which helps to distribute the finally extracted watershed boundaries along the ridgeline rather than shrinking towards the inward slope.

[0073] Further, S6 includes the following steps: If the candidate flooding area does not overflow, the candidate flooding area is determined as the currently determined depression area, and the area and perimeter of the candidate flooding area are calculated to update the geometric parameters of the currently determined depression area for the calculation of the area growth rate and perimeter growth rate in the next iteration.

[0074] Specifically, if the area growth rate in step S4 does not reach the preset area growth rate threshold and the perimeter growth rate does not reach the preset perimeter growth rate threshold, it indicates that after the current water level rises, the water flow only expands normally along the slope within the depression, and the flooded area does not enter adjacent depressions through any pass. At this time, there is no need to perform pass location and blocking operations, and the newly flooded area is identified as the reference area for the next iteration.

[0075] Calculate the area and perimeter of the candidate flooded areas to update the geometric parameters of the currently determined depression areas. Specifically, assign the area of ​​the candidate flooded areas to the area of ​​the currently determined depression areas, and assign the perimeter of the candidate flooded areas to the perimeter of the currently determined depression areas. The purpose of updating the geometric parameters is to provide a calculation basis for the next round of overflow determination. After the state update is complete, return to step S3 to proceed to the next iteration.

[0076] Further, S7 includes the following steps: Repeat steps S3 to S6, increasing the current water level value by the preset water level rise step size in each iteration until the current water level value reaches the preset termination elevation and the iteration terminates; wherein the preset termination elevation is greater than or equal to the highest elevation value in the digital elevation model data; Repeat steps S3 to S6, increasing the current water level value by a preset water level rise step size in each iteration until the current water level reaches a preset termination elevation, at which point the iteration terminates. The preset termination elevation is greater than or equal to the highest elevation value in the digital elevation model data to ensure that all potential passes are detected. Each iteration sequentially executes step S3 for generating the flooded area, step S4 for determining the overflow, and step S5 for processing the overflow channel or step S6 for updating the state, simulating the process of water being injected from the bottom of the depression until the water overflows the entire ridge.

[0077] The closed loop formed by all the cells marked as blocked in the blocking array is used as the watershed boundary of the candidate depression. Based on the watershed boundary and the digital elevation model data, the spatial morphology quantification results of the candidate depressions are generated.

[0078] As a preferred embodiment, generating the spatial morphology quantification result of the candidate depression based on the watershed boundary and the digital elevation model data includes the following steps: Step 1, Preparation of data for morphological index calculation: After obtaining the watershed boundary, organize the relevant data required for the spatial morphological quantification index calculation. The relevant data includes the watershed line elevation unfolding data, watershed boundary polygon, lowest pass contour polygon, spherical elevation matrix, and the center of the bottom of the sphere. The process of obtaining each data is as follows.

[0079] (1) Watershed line elevation unfolding data: Select any cell on the watershed boundary as the starting point, and count the elevation values ​​of each cell on the watershed boundary in a clockwise direction. Use the plane distance from the center point of each cell to the starting cell as the horizontal axis and the elevation value of each cell as the vertical axis to generate watershed line elevation unfolding data. (2) Watershed boundary polygon: Connect the center points of each pixel on the watershed boundary in sequence to form a closed watershed boundary polygon. The watershed boundary polygon is used to characterize the overall outline range of the depression. (3) Lowest pass contour polygon: The elevation value of the lowest pass is determined according to the elevation unfolded data of the watershed line. The area enclosed by the contour lines corresponding to the elevation value of the lowest pass in the depression is extracted to form the lowest pass contour polygon. The lowest pass contour polygon is used to characterize the core area range where the radio telescope's reflecting sphere can be deployed. (4) Spherical Elevation Matrix: Given a spherical radius and a cap radius, calculate the corresponding elevation value for each pixel location within the cap area based on the spherical equation, thus constructing the radio telescope's reflecting spherical elevation matrix. Pixel locations within the cap area are recorded as spherical elevation values, while pixel locations outside the cap area are marked as null values. Based on the reflecting spherical elevation matrix, a valid area mask matrix is ​​generated. Pixel locations corresponding to valid elevation values ​​are assigned preset valid marker values, while pixel locations corresponding to null values ​​are assigned preset invalid marker values. In the formula, Let be the spherical elevation value of the radio telescope's reflecting spherical elevation matrix at position (i,j); and Let (i,j) be the position coordinates of (i,j) in the matrix, NaN is a null value, R is the radius of the sphere, and r is the radius of the spherical cap.

[0080] (5) Center of the bottom of the sphere: 1. Calculate the difference between the highest elevation value in the spherical elevation matrix and the average elevation value of the spherical elevation matrix, and denote it as the spherical vertical average elevation difference; 2. Calculate the difference matrix: Using the spherical elevation matrix as a sliding window, traverse all pixels within the contour polygon of the lowest pass one by one, and calculate the difference between the maximum elevation and the mean of the adjacent elevation data of each pixel within the sliding window. This difference is denoted as the average vertical elevation difference of the terrain. If the average vertical elevation difference of the terrain is greater than the average vertical elevation difference of the spherical elevation matrix, it indicates that the terrain at the pixel location is relatively open, which can fit well with the radio telescope sphere and ensure that the height of the sphere does not exceed the height of the lowest pass. Subtract the averaged matrix of the digital elevation model data within the sliding window from the averaged matrix of the spherical elevation matrix to obtain the difference matrix. Otherwise, it indicates that the terrain at the pixel location is relatively close to the elevation of the lowest pass. Shift the elevation of the spherical elevation matrix until the elevation of the sphere equals the elevation of the lowest pass to obtain the adjusted spherical elevation matrix. Subtract the averaged matrix of the digital elevation model data within the sliding window from the adjusted spherical elevation matrix to obtain the difference matrix. 3. Calculate the center position of the sphere: The positive difference value in the difference matrix corresponds to the excavation volume, and the negative difference value corresponds to the filling volume. The sum of the absolute values ​​of the excavation volume and the filling volume is taken as the total volume of the current pixel position. After traversal, the total volume of all pixel positions within the polygon range of the lowest arch contour line is obtained, and the pixel position with the smallest total volume is selected as the center position of the sphere. 4. Calculate the elevation of the bottom of the sphere: If the average vertical elevation difference of the terrain corresponding to the center position of the sphere is greater than the average vertical elevation difference of the sphere, the average value of the adjacent elevation data of the pixel in the sliding window is used as the average elevation of the spherical elevation matrix to calculate the elevation of the bottom of the sphere; otherwise, the lowest pass elevation value is used as the spherical crown elevation to calculate the elevation of the bottom of the sphere. 5. Determine the coordinates of the center of the bottom of the sphere: Determine the center of the bottom of the sphere based on the position of the center of the sphere and the elevation of the bottom of the sphere.

[0081] In the search process for the center of the bottom of the sphere, the difference matrix between the spherical elevation matrix and the terrain is calculated using a mean-removing method. This allows us to obtain the optimal matching position between the sphere and the terrain in the vertical direction of the sliding window location, which is the spherical elevation position where the engineering excavation volume is minimized. Simultaneously, the elevation of the spherical crown is constrained by the lowest pass elevation value to ensure the enclosure and electromagnetic quietness of the depression after the radio telescope is built, thus completing the search for the center of the bottom of the sphere in the vertical direction. Furthermore, the total volume of all pixels within the contour polygon range of the lowest pass is calculated one by one, and the minimum value is selected to complete the search for the center of the bottom of the sphere in the horizontal direction. Thus, the optimal spatial location for the construction of the radio telescope is determined in three-dimensional space.

[0082] Step 2, Calculation of Spatial Morphology Quantification Index: Based on the data prepared in Step 1, calculate the following spatial morphology quantification index. The spatial morphology quantification index is divided into five categories according to the evaluation purpose: degree of enclosure, size characteristics, shape characteristics, roundness characteristics, and local matching. It is supplemented by an engineering economic index to form a quantitative evaluation system for the spatial morphology of depressions.

[0083] The first category is the enclosure index, which is used to assess the shielding ability of candidate depressions against external electromagnetic interference. The cosmic radio signals received by large radio telescopes are extremely weak, and the electromagnetic quietness of the site directly determines the observation performance. The mountains surrounding the depression form a natural shield against anthropogenic electromagnetic radiation. The higher the enclosure degree of the depression and the fewer the overflow channels (passes), the more difficult it is for external electromagnetic interference to penetrate. Therefore, the enclosure coefficient, the number of peaks, and the number of passes are used to characterize the enclosure degree of the depression.

[0084] Index 1. Closure Coefficient: Based on the watershed elevation unfolded data, the closure coefficient refers to the ratio of the area of ​​the closed region in the watershed elevation unfolded data to the total area. The closure coefficient is obtained by using the area of ​​the closed region calculated from the lowest elevation point of the candidate depression as the numerator and the total area of ​​the watershed elevation unfolded data as the denominator. The closure coefficient reflects, to a certain extent, the degree to which the candidate depression is protected from external environmental interference such as electromagnetic interference and wind loads.

[0085] Index 2. Number of Peaks: Based on the elevation data of the watershed line, calculate the difference between the elevation value at each location and the previous location, obtaining a difference sequence. Iterate through the difference sequence, and increment the number of peaks by one when the difference changes from a positive value to a negative value.

[0086] Index 3. Number of passes: Based on the difference sequence obtained from the elevation expansion data of the watershed line, the number of passes is incremented by one when the difference changes from negative to positive.

[0087] The second category is size characteristic indices, used to evaluate the range of radio telescope sizes that can be accommodated in candidate depressions. The depth of the depression determines the degree to which the reflecting sphere is embedded downwards, the axial length determines the maximum aperture that the reflecting sphere can be installed, and the net depth determines the effective clearance for sphere installation. Therefore, the maximum depth, net depth, major axis length, and minor axis length of the depression are used to characterize the size characteristics of the depression.

[0088] Index 4. Maximum Depth: The elevation difference between the highest elevation value within the boundary of the watershed and the lowest elevation value of the candidate depression.

[0089] Index 5. Net depth: The elevation difference between the lowest pass elevation and the lowest point elevation of the candidate depression within the boundary of the watershed.

[0090] Index 6. Major axis length of the depression: The length of the longest section among all vertical sections passing through the center of the bottom of the sphere within the polygonal boundary of the watershed.

[0091] Index 7. Minor axis length of the depression: The length of the section perpendicular to the major axis within the polygonal boundary of the watershed.

[0092] The third category, morphological characteristic indices, is used to characterize the regularity of the overall outline of candidate depressions and the usable space range. The diameter of the watershed great circle measures the size of the usable area from the overall outline scale of the depression, while the diameter of the pass great circle measures the range that can be deployed at the bottom of the sphere from the scale of the lowest overflow channel. Together, they reflect the regularity and expandability of the depression space.

[0093] Index 8. Great circle diameter of the watershed: The diameter of the largest inscribed circle of the boundary polygon of the watershed is the great circle diameter of the watershed.

[0094] Index 9. Great circle diameter of the pass: Calculate the largest inscribed circle of the polygon with the lowest contour line of the pass, and the diameter of the largest inscribed circle is the great circle diameter of the pass.

[0095] The fourth category is the roundness characteristic index, which is used to evaluate how close the shape of the candidate depression is to the ideal circle and the slope conditions. The higher the roundness coefficient, the better the fit between the depression outline and the circular reflective surface, which is more conducive to the concentric arrangement of the reflective sphere; the slope reflects the steepness of the slope inside the depression, which directly affects the construction difficulty and drainage organization.

[0096] Index 10. Roundness Coefficient: Calculate the ratio of the area of ​​the largest inscribed circle within the watershed boundary polygon to the area of ​​the watershed boundary polygon, and the ratio of the area of ​​the largest inscribed circle within the lowest pass contour polygon to the area of ​​the lowest pass contour polygon. Take the average of the two ratios as the roundness coefficient. The roundness coefficient is used to reflect the shape characteristics of the candidate depression.

[0097] Index 11. Slope: Calculate the slope value of all pixels within the watershed boundary and take the average of the slope values ​​of all pixels as the slope value of the candidate depression.

[0098] The fifth category, the local matching index, is used to quantitatively evaluate the degree of matching between the local topography of the candidate depression and the radio telescope's reflecting sphere from the perspective of spherical fit. The elevation value of the center of the bottom of the sphere is assigned to the spherical elevation matrix to obtain the translated spherical elevation matrix.

[0099] Index 12. Directional Index: Using the pixel at the center of the bottom of the sphere and the effective region mask matrix, a matching region matrix is ​​extracted between the candidate depression center region and the radio telescope reflecting sphere from the digital elevation model data. The row and column gradient components of the matching region matrix are calculated. Based on the squares of the row, column, and intersection gradient components, the directional index is calculated. The directional index analyzes the spatial morphology of the candidate depression from the perspective of anisotropy in the slope direction.

[0100] in, As a directional index, , These are the averages of the squares of the row-direction gradient components and the squares of the column-direction gradient components of the matching region matrix over all pixels. This is the average of the product of the row and column gradient components over all pixels. A larger directional index indicates stronger anisotropy in the slope direction of the candidate depression.

[0101] Index 13. Structural Similarity Index: The structural similarity index consists of three components: a brightness comparison function, a contrast comparison function, and a structure degree comparison function. The brightness comparison function is calculated based on the local mean of the digital elevation model (DEM) data and the local mean of the radio telescope's reflecting spherical elevation matrix. The contrast comparison function is calculated based on the local standard deviation of the DEM data and the local standard deviation of the reflecting spherical elevation matrix. The structure degree comparison function is calculated based on the local covariance of the DEM data and the reflecting spherical elevation matrix. The three comparison functions are combined with equal weights, and a constant term is introduced to prevent anomalies in the denominator, resulting in the structural similarity index. Based on the structural similarity index calculated using the matching region matrix and the translated spherical elevation matrix, it is used to describe the degree of matching between the local morphology of the candidate depression and the radio telescope's reflecting spherical surface.

[0102] The sixth category is engineering economic indexes, which are used to evaluate the amount of filling and excavation work and the economic efficiency of depression leveling construction.

[0103] Index 14. Volume Index: The volume index is calculated by successively subtracting the elevation values ​​of corresponding pixels in the matching region matrix of the digital elevation model and the translated spherical elevation matrix. All positive differences are summed and multiplied by the ground area corresponding to a single pixel to obtain the fill volume; the absolute values ​​of all negative differences are summed and multiplied by the ground area corresponding to a single pixel to obtain the cut volume. The sum of the fill and cut volumes is the volume index. The volume index is used to characterize the economic benefits of constructing a radio telescope in a candidate depression; the smaller the volume index, the lower the relative construction cost.

[0104] The local matching index and the volume index are calculated based on the spherical bottom center and spherical elevation matrix determined in the first step: taking the spherical bottom center as the spherical placement position, according to the effective region mask matrix, the matching region matrix corresponding to the radio telescope reflecting spherical surface is extracted from the digital elevation model data, and used for the calculation of the directional index, structural similarity index and volume index.

[0105] It is worth noting that in the evaluation and screening of candidate depressions, existing technologies are mainly based on scoring mechanisms such as the analytic hierarchy process (AHP) and rely on technical personnel to verify the geological conditions of each candidate depression one by one, including water system distribution, transportation accessibility, statistics of the number of surrounding peaks, and interpretation of the location of mountain passes. The influence of human subjective judgment is strong, which can easily lead to unreasonable screening results. For large-scale depression screening scenarios, there are problems of large workload and difficulty in unifying evaluation standards.

[0106] To address this, the present invention also provides a comprehensive evaluation method for candidate depressions based on immersion simulation and a two-layer decision-making model. This method involves acquiring spatial morphology quantification indices and volume indices of multiple candidate depressions to form a spatial morphology index dataset. The dual-criteria comprehensive ranking is implemented based on a two-layer decision-making model, which includes a morphological distance layer and an economic attribution layer. The comprehensive ranking of multiple candidate depressions is achieved through objective weighting using information entropy. Specific embodiments are described below.

[0107] Step 1: Taking the Dawodang Depression in Guizhou (the site of the 500-meter Aperture Spherical Radio Telescope (FAST)) as the benchmark depression, denoted as K0 (106.8560°E, 25.6520°N), and using 41 candidate depressions as analysis samples, numbered W1~W41, the above steps are used to calculate 13 spatial morphological quantification indices and 1 volume index for both the benchmark depression and each candidate depression, forming a spatial morphological index dataset. That is, the spatial morphological parameter set is generated based on the watershed boundaries and digital elevation model data of each candidate depression.

[0108] It should be noted that the spatial morphology quantification index system constructed in this invention covers multiple aspects such as the degree of depression closure, size characteristics, morphological characteristics, roundness characteristics, local matching, and engineering economics, and can comprehensively characterize the spatial morphology of candidate depressions. However, precisely because the index system is relatively comprehensive, information overlap and linear correlation inevitably exist between the various indices. For example, the number of mountain passes and the number of mountain peaks describe the troughs and peaks of the same watershed elevation development curve, and the length of the depression's major axis and the diameter of the watershed great circle both change in the same direction with the size of the depression. If the morphological distance is directly calculated based on the original indices or a volumetric regression model is established, the same morphological features reflected by the related indices will be repeatedly counted, leading to similarity distortion and unstable regression coefficients. Therefore, this invention first performs collinearity removal and principal component dimensionality reduction on the spatial morphology index dataset before sorting, converting the mutually correlated original indices into mutually independent principal components, and then constructing a dual criterion based on the principal components, thereby ensuring the objectivity and stability of the sorting results.

[0109] Step 2: Collinearity removal and principal component dimensionality reduction are performed on the spatial morphology parameter set to obtain multiple principal components after dimensionality reduction. The variance inflation coefficient (VIF) is calculated for the spatial morphology quantification index in the spatial morphology index dataset, and redundant indicators with VIF greater than 10 are removed. In this example, the number of mountain passes and the number of mountain peaks are completely collinear (VIF=∞), and the VIF of the long axis length of the depression is 23.77. After removing these two, 11 indicators are retained for subsequent analysis.

[0110] Principal component analysis was performed on the retained indicators to extract principal components with a cumulative variance contribution rate greater than 90%. In this example, a total of 6 principal components were extracted, with a cumulative variance contribution rate of 90.85% (PC1=32.14%, PC2=21.36%, PC3=11.29%, PC4=10.19%, PC5=8.10%, PC6=7.76%).

[0111] Furthermore, the principal component loading matrix is ​​rotated using Varimax to improve the interpretability of the principal components. After rotation, the dominant indices of each principal component are as follows: PC1 is the scale factor, consisting of the watershed great circle diameter (0.51), the depression minor axis length (0.50), and the number of peaks (0.43); PC2 is the depth morphology factor, consisting of the net depth (0.58), the slope (0.56), and the maximum depth (0.42); PC3 is the pass scale factor, consisting of the pass great circle diameter (0.94); PC4 is the directionality factor, consisting of the directionality index (0.89) and the structural similarity index (0.395); PC5 is the roundness factor, consisting of the roundness coefficient (0.96); and PC6 is the closure factor, consisting of the closure coefficient (-0.76) and the number of peaks (-0.57).

[0112] Step 2, Morphological Distance Layer (Criterion A): Using the benchmark depression K0 as a reference, calculate the Mahalanobis distance from each candidate depression to the benchmark depression in the principal component space composed of multiple principal components. Determine the Criterion A score for each candidate depression based on the Mahalanobis distance; the smaller the Mahalanobis distance of a candidate depression, the higher its Criterion A score. Using the benchmark depression K0 as a reference, calculate the Mahalanobis distance from each candidate depression to the benchmark depression K0 in the principal component space, and convert the Mahalanobis distance into morphological similarity as the first criterion (Criterion A). The essence of Mahalanobis distance is to eliminate the influence of dimensions and correlation under the premise of considering the covariance structure between variables, thereby obtaining a true statistical distance; morphological similarity is obtained by transforming the Mahalanobis distance by its reciprocal, with a value range of (0,1]. The closer to 1, the more similar the spatial morphology of the candidate depression is to the benchmark depression. Step 2 specifically includes the following sub-steps: Step 2.1, Distance Difference Vector Calculation: For each candidate depression, calculate the principal component score vector of the i-th candidate depression. Principal component score vector of benchmark depression difference vector In this example, each candidate depression and the benchmark depression are constructed into a 6-dimensional score vector using the scores of 6 principal components. Step 2.2, Mahalanobis distance calculation: Calculate the Mahalanobis distance from each candidate depression to the benchmark depression K0 in the principal component space using the following formula: In the formula, Let be the Mahalanobis distance from the i-th candidate depression to the benchmark depression K0 in the principal component space; the superscript T denotes the vector transpose; S is the sample covariance matrix. Let K0 be its inverse matrix. In this example, with K0 as the benchmark, the top 5 candidate depressions with the most similar shapes are: W4 (Mahathano distance 2.23, shape similarity 0.310), W37 (2.89, 0.257), W32 (3.10, 0.244), W1 (3.60, 0.217), and W10 (3.94, 0.202).

[0113] Step 2.3, Morphological Similarity Conversion: Convert the Mahalanobis distance into morphological similarity using the following formula via inverse transformation: In the formula, Let be the morphological similarity score of the i-th candidate depression, i.e., the score of criterion A, with a value range of (0,1]. The closer to 1, the more similar the spatial morphology is to the benchmark depression. In this example, the morphological similarity scores of the top 5 candidate depressions are W4 (0.310), W37 (0.257), W32 (0.244), W1 (0.217), and W10 (0.202), respectively. At the same time, the candidate depressions can be automatically grouped based on the KMeans clustering silhouette coefficient (preferably 5 groups in this example, with a silhouette coefficient of 0.2878) to assist in analyzing the type distribution of the candidate depressions.

[0114] Step 3, Economic Attribution Layer (Criterion B): Using the scores of each candidate depression on the multiple principal components as independent variables and the volume index of each candidate depression as dependent variables, a partial least squares regression mapping is established. Based on the partial least squares regression mapping, the predicted volume of each candidate depression is output. Based on the predicted volume, the Criterion B score of each candidate depression is determined, with a value range of [0,1]. The closer to 1, the smaller the predicted volume and the better the economy, which is used as the second criterion (Criterion B). At the same time, the regression residuals are output for the labeling of outlier samples.

[0115] It should be noted that partial least squares regression was used because the number of candidate depressions in this example is limited (only 41), resulting in a small sample size compared to the number of exponents. If ordinary least squares regression were used under small sample conditions, it would easily lead to unstable fitting coefficients and weak generalization ability. Partial least squares regression, while extracting components, also considers the explanatory power of the independent variables on the volume exponent of the dependent variable, enabling the establishment of a robust volume prediction mapping under small sample conditions. Combined with leave-one-out cross-validation, it can further suppress overfitting. With a sufficiently large sample size, other suitable regression methods can achieve the same results.

[0116] Step 3 specifically includes the following sub-steps: Step 3.1, Volume Prediction: Partial least squares regression predicts the volume index from the principal component scores of each candidate depression using the following formula: In the formula, Let be the predicted volume of the i-th candidate depression; The partial least squares component score for the i-th candidate depression; Let be the regression coefficient of the k-th component. is a constant term; m is the number of optimal components determined by leave-one-out cross-validation. In this example, the coefficient of determination is... =0.3616; Step 3.2, Goodness-of-Performance Transformation: The predicted quantities are inversely standardized using the following formula to obtain the goodness-of-performance: In the formula, Let be the goodness of measure for the i-th candidate depression, i.e., the score of criterion B; , These represent the maximum and minimum predicted volumes for all candidate depressions, respectively. The value range is [0,1], and the closer it is to 1, the smaller the prediction volume and the better the economy; Step 3.3, Key Pattern Index Importance Analysis: Calculate the variable projection importance index (VIP) for each principal component using the following formula: In the formula, is the variable projection importance index for the j-th principal component; p is the number of principal components; Let h be the variance explained by the partial least squares component for the volume index; Let h be the weight of the h-th partial least squares component on the j-th principal component; the left side of the equation is VIP. j The square root of the square gives the VIP. j In this example, PC4 (1.437) and PC2 (1.410) contributed the most to the VIP scores of each principal component, indicating that directionality and depth morphology are the most critical factors for volume prediction. Further analysis of the importance of quantitative indices for spatial morphology in predicting volume using the following formula: In the formula, The importance score for the l-th spatial morphology quantification index; Let l be the absolute value of the rotational load on the j-th principal component for the l-th spatial morphology quantification index. In this example, the top four most important morphology indices are, in order: directionality index (1.426), structural similarity index (1.371), net depth (1.328), and roundness coefficient (1.043).

[0117] Step 4, Decision Ranking Layer: After normalizing the scores of the two criteria, information entropy is used to determine the criterion weights. The weighted sum is then used to obtain the comprehensive score, and the data are ranked according to the comprehensive score. The reason for using information entropy weighting is that the criterion weights are entirely determined by the objective differences in the candidate low-lying data under each criterion, without relying on expert experience. This avoids the interference of human judgment on the ranking results in subjective weighting methods such as the analytic hierarchy process, ensuring consistent ranking conclusions for different personnel and different batches. Step 4 specifically includes the following sub-steps: Step 4.1, Score Normalization: Normalize the scores of criterion A and criterion B for each candidate depression to the [0,1] interval using the following formula: In the formula, Let be the normalized score of the i-th candidate depression under criterion j; The original score of the i-th candidate depression under criterion j; , These are the maximum and minimum scores of all candidate depressions under criterion j, respectively. Step 4.2, Information Entropy Calculation: Calculate the information entropy values ​​for the two criteria using the following formulas: In the formula, The information entropy value of criterion j; N is the number of candidate depressions; Let be the proportion of the normalized score of the i-th candidate depression under criterion j to the sum of the scores of all candidate depressions. When =0, let =0. The smaller the entropy value, the greater the data difference and the stronger the distinguishing ability of each candidate depression under the corresponding criterion; Step 4.3, Determine the weights of the criteria: Calculate the weights of each criterion using the following formula: In the formula, For the weight of criterion j, For the discriminative power of criterion j, the weights satisfy... + =1 and w j Greater than 0. In this example, the weight of criterion A (morphological similarity) is 0.5281, and the weight of criterion B (quantitative goodness) is 0.4719; Step 4.4, Overall Score and Ranking: The normalized scores of the two criteria are weighted and summed using the following formula to obtain the overall score for each candidate depression, and then ranked from highest to lowest overall score: In the formula, Let be the overall score of the i-th candidate depression. , These are the normalized scores for Criterion A and Criterion B, respectively. A higher overall score indicates that the candidate depression performs better in terms of both morphological similarity and volume economy. In this example, the top 5 candidate depressions are W20 (0.786), W6 (0.777), W38 (0.761), W15 (0.760), and W27 (0.754).

[0118] This invention constructs a spatial morphology quantification index based on watershed boundaries, realizing the quantitative expression of the morphological characteristics of candidate depressions. Furthermore, it proposes to construct a dual criterion of morphological similarity and volume superiority based on benchmark depressions, and to achieve comprehensive ranking of multiple candidate depressions through objective weighting by information entropy. This enables rapid and quantitative selection of candidate depressions, avoids interference from subjective weighting methods such as traditional analytic hierarchy process, and improves the objectivity, reproducibility, and engineering applicability of large-scale candidate depression selection.

[0119] The above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any equivalent structural or procedural transformations made based on the description and drawings of the present invention, or direct or indirect applications in other related technical fields, are similarly included within the scope of patent protection of the present invention.

Claims

1. A comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision-making model, characterized in that, Includes the following steps: S1, Obtain the spatial morphological index dataset of the benchmark depression and multiple candidate depressions. The spatial morphological index dataset includes multiple spatial morphological quantification indices and volume indices generated based on the watershed boundaries and digital elevation model data of each candidate depression. S2, perform collinearity removal and principal component dimensionality reduction on the spatial morphology quantification index in the spatial morphology index dataset to obtain multiple principal components after dimensionality reduction, and obtain the principal component score vectors of the benchmark depression and each candidate depression based on the multiple principal components; S3, Construct a two-layer decision-making model, which includes a morphological distance layer and an economic attribution layer; S4, in the morphological distance layer, taking the benchmark depression as a reference, calculate the Mahalanobis distance of each candidate depression to the benchmark depression in the principal component space composed of multiple principal components, and determine the morphological similarity of each candidate depression based on the Mahalanobis distance, and use the morphological similarity as criterion A score. S5, in the economic attribution layer, the principal component scores of each candidate depression on the multiple principal components are used as independent variables, and the volume index of each candidate depression is used as dependent variable. A partial least squares regression mapping is established. Based on the partial least squares regression mapping, the predicted volume of each candidate depression is output, and the volume goodness of each candidate depression is determined based on the predicted volume. The volume goodness is used as the criterion B score. S6, normalize the scores of criterion A and criterion B, and use information entropy weighting to determine the criterion weights of criterion A and criterion B. S7. Based on the normalized score of criterion A, the normalized score of criterion B, the criterion weight of criterion A, and the criterion weight of criterion B, the comprehensive score of each candidate depression is obtained, and the multiple candidate depressions are sorted according to the comprehensive score, and the comprehensive evaluation result of the candidate depressions is output.

2. The comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision-making model according to claim 1, characterized in that, The spatial morphology quantification index includes the closure coefficient, number of peaks, number of passes, maximum depth, net depth, length of the major axis of the depression, length of the minor axis of the depression, diameter of the great circle of the watershed, diameter of the great circle of the pass, roundness coefficient, slope, directionality index, and structural similarity index.

3. The comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision-making model according to claim 1, characterized in that, The S2 step of performing collinearity removal and principal component dimensionality reduction on the spatial morphology index dataset includes the following steps: S21, Calculate the variance inflation coefficient of the spatial morphology quantification index in the spatial morphology index dataset; S22, remove redundant indicators with a variance inflation coefficient greater than 10; S23, perform principal component analysis on the retained spatial morphology quantification index and extract principal components with a cumulative variance contribution rate greater than 90%.

4. The comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision-making model according to claim 3, characterized in that, S2 further includes the following steps: The principal component loading matrix is ​​rotated to maximize variance to improve the interpretability of the principal components, and the dominant indices of each principal component are determined based on the rotated principal component loading matrix.

5. The comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision-making model according to claim 1, characterized in that, In step S4, the morphological similarity of each candidate depression is determined based on the Mahalanobis distance, including the following steps: S41, For each candidate depression, calculate the difference vector between the principal component score vector of the candidate depression and the principal component score vector of the benchmark depression; S42, calculate the Mahalanobis distance from the candidate depression to the benchmark depression in the principal component space based on the difference vector, the sample covariance matrix, and the inverse matrix of the sample covariance matrix; S43, the Mahalanobis distance is converted into morphological similarity through inverse transformation; Among them, the smaller the Mahalanobis distance between the candidate depressions, the higher the morphological similarity of the candidate depressions.

6. The comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision-making model according to claim 1, characterized in that, The partial least squares regression mapping established in S5 includes the following steps: S51, the principal component scores of each candidate depression on the multiple principal components are used as the independent variables of the partial least squares regression mapping; S52, the volume index of each candidate depression is used as the dependent variable of the partial least squares regression mapping; S53, the optimal number of components in the partial least squares regression mapping is determined by leave-one-out cross-validation; S54. Based on the optimal number of components, establish a partial least squares regression mapping that predicts the volume index by the principal component scores of each candidate depression.

7. The comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision-making model according to claim 1, characterized in that, The step S5, which determines the quantitative superiority of each candidate depression based on the predicted quantitative volume, includes the following steps: S55, perform inverse standardization on the predicted volume to obtain the volume quality of each candidate depression; The range of the quantity goodness is [0,1]. The closer the quantity goodness of a candidate depression is to 1, the smaller the predicted quantity of the candidate depression and the better its economy.

8. The comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision-making model according to claim 1, characterized in that, S5 further includes the following steps: The regression residuals are output based on the partial least squares regression mapping, and the candidate depressions are marked as abnormal samples based on the regression residuals.

9. The comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision-making model according to claim 1, characterized in that, S5 further includes the following steps: Calculate the variable projection importance index for each principal component, and determine the importance score of each spatial morphology quantification index for volume prediction based on the variable projection importance index for each principal component and the absolute value of the rotational load of each spatial morphology quantification index on each principal component.

10. The comprehensive evaluation method for candidate depressions based on flooding simulation and a two-layer decision-making model according to claim 1, characterized in that, The S6 step, which uses information entropy weighting to determine the weights of criteria A and B, includes the following steps: S61, normalize the scores of criterion A and criterion B for each candidate depression to the interval [0,1] respectively; S62, calculate the information entropy values ​​of criteria A and criteria B respectively; S63, Determine the discriminative power of criterion A based on the information entropy value of criterion A, and determine the discriminative power of criterion B based on the information entropy value of criterion B; S64. Determine the criterion weights of criterion A and criterion B based on the discriminative power of criterion A and criterion B. The sum of the criterion weights of criterion A and criterion weights of criterion B is 1. The overall score is obtained by weighted summation of the normalized scores of criterion A and criterion B.

Images