Karst large round depression identification method based on DEM data

Abstract

The invention belongs to the technical field of geographic information, and discloses a Karst large round depression identification method based on DEM data. The identification method comprises the following specific steps of S1, DEM data preprocessing; S2, standard template construction; S3, sliding window mechanism design; S4, round depression form template construction; through a standardized data processing flow and an automatic coordinate conversion technology, full-flow automation from data input to identification output is realized, the identification accuracy is high, benefit from an advanced discrimination method integrating a plurality of characteristic indexes and accurate space positioning is achieved, a unified quantitative evaluation standard is established for ensuring the reliability of a result, and the method is suitable for popularization and application. A complete coordinate conversion and data processing mechanism is built in, in addition, wide applicability is shown, and the method can flexibly adapt to DEM data of different geographic reference systems and resolutions.

Description

Technical Field

[0001] The present invention belongs to the field of geographic information technology, and specifically relates to a method for identifying large karst dolines based on DEM data. Background Art

[0002] Traditional methods relying on manual visual interpretation not only involve a large amount of work, but also lack objectivity and repeatability in the interpretation results, making it difficult to meet the requirements of rapid assessment in engineering site selection; the parameter settings of existing automatic recognition methods are too empirical and cannot adapt to the complex and diverse terrain features in karst areas, affecting the scientific nature of site selection decisions; at the same time, existing technologies lack a reliable evaluation system and a complete visualization analysis function, especially in engineering applications such as multi-site comparison and dynamic monitoring, there are obvious deficiencies, and it is difficult to provide effective support for engineering site selection. The traditional method of manually visually interpreting large karst dolines is laborious and the results are subjective and non-repeatable, unable to meet the rapid assessment requirements of engineering site selection. Although existing automatic recognition technologies have made some progress, their parameter settings are too dependent on experience and it is difficult to cope with the complex terrain in karst areas, affecting the scientific nature of site selection decisions. In addition, existing technologies also lack a reliable evaluation system and a comprehensive visualization analysis function, especially there are shortcomings in multi-site comparison and dynamic monitoring, and it is difficult to provide solid technical support for engineering site selection. Summary of the Invention

[0003] The purpose of the present invention is to provide a method for identifying large karst dolines based on DEM data to solve the problems raised in the above background art.

[0004] To achieve the above purpose, the present invention provides the following technical solution: A method for identifying large karst dolines based on DEM data, and the specific steps of the identification method are as follows:

[0005] S1: DEM data preprocessing: Read DEM data, clean and process it, and convert the coordinate system to ensure that the data accurately corresponds to the ground surface;

[0006] S2: Standard template construction: Define hemisphere parameters, convert coordinates, perform template matching, and determine the optimal model parameters;

[0007] S3: Design a sliding window mechanism: Design windows and overlapping areas to ensure the complete identification of dolines and balance efficiency and accuracy;

[0008] S4: Construct a doline morphology template: Establish a coordinate grid, calculate the distance field, and construct an optimized doline morphology model;

[0009] S5: Implement a multi-feature matching algorithm: Implement a multi-feature matching algorithm, comprehensively analyze morphology, local minimum, and area-depth, and obtain a matching score;

[0010] S6: Result verification and visualization enhancement: 3D reconstruction and fitting analysis, 2D contour analysis, grayscale image analysis, comprehensively evaluate candidate circular depressions;

[0011] S7: Parameter setting and threshold selection: The parameters and thresholds need to be set precisely, including a radius of 350 and a flatness ratio of 0.6, to ensure the morphology and accuracy match.

[0012] Preferably, the specific steps of DEM data preprocessing in S1 are as follows:

[0013] Step 1: Data reading and verification

[0014] Read DEM data in.tif format

[0015] Use the readgeoraster function to read DEM data in.tif format;

[0016] Verify data format and georeference information

[0017] Verify whether the read data is in a two-dimensional matrix format and check whether the georeference information is correct;

[0018] Step 2: Data cleaning

[0019] Identify and process invalid values

[0020] Use functions in the NumPy library to identify and process NaN values and outliers, replace NaN values with the average of neighboring pixels or perform other reasonable interpolation processing;

[0021] Convert data type and smooth processing

[0022] Ensure the data type is double precision and perform data smoothing processing;

[0023] Step 3: Coordinate system conversion

[0024] Based on the average radius of the Earth being approximately 6371000 meters and combined with the resolution of the digital elevation model, that is, each pixel represents 30 meters of actual surface distance, convert the geographic coordinates to actual distance coordinates on the surface. The conversion process takes into account the curvature of the Earth and the spatial resolution of the DEM data.

[0025] Preferably, the specific steps of standard template construction in S2:

[0026] Step 1: Basic parameter definition and model establishment

[0027] Define the basic parameters of the hemisphere: Through data analysis, the optimal hemisphere parameters are determined, radius: 350 meters; height: 210 meters; flatness ratio: Number of grid subdivisions: 100;

[0028] Define the parameters of the coordinate system; azimuth angle (θ): θ ∈ [0, 2π]; polar angle (φ):

[0029] Step 2: 3D coordinate transformation and template matching

[0030] Perform 3D coordinate transformation:

[0031] The standard template describes the surface of the hemisphere using parametric equations, and the flatness parameter is introduced to achieve an accurate fit to the actual terrain:

[0032]

[0033] where X, Y, Z: positions in the space rectangular coordinate system; r: radius of the hemisphere; f: flatness; H: height of the reference plane;

[0034] Apply the 3D coordinate transformation results for template matching:

[0035] After completing the 3D coordinate transformation, use the obtained 3D coordinates for template matching; mainly compare the actual terrain data with the hemisphere model precisely. By calculating errors, performing correlation analysis, or using machine learning algorithm methods, to find the parameter configuration of the hemisphere model that best matches the actual terrain. Finally, the goal is to determine the hemisphere model parameters that can best reflect the characteristics of the actual terrain, thus laying a solid foundation for subsequent terrain analysis and modeling.

[0036] Preferably, the specific steps for designing the sliding window mechanism in S3 are as follows:

[0037] Step 1: Window size design

[0038] Determine the expected radius of the circular depression and the resolution of the DEM data:

[0039] According to the analysis requirements, clarify the radius of the expected circular depression;

[0040] Obtain the resolution of the DEM data, that is, the actual distance represented by each pixel;

[0041] Calculate the window size:

[0042] Use the window size calculation formula, divide the expected diameter of the circular depression by the resolution of the DEM data, and round up to ensure that the window can completely cover the expected circular depression shape;

[0043] The specific formula of the window size calculation formula is as follows:

[0044]

[0045] Step 2: Overlap area design

[0046] Set the overlapping ratio:

[0047] To avoid recognition deviation caused when the circular depression is located at the window boundary, an overlapping scanning mechanism is introduced; the formula is as follows:

[0048]

[0049] Determine the overlapping ratio between adjacent windows, and the overlapping ratio is 87.5%;

[0050] Calculate the step size and set the scanning mechanism:

[0051] Calculate the step size according to the overlapping ratio. The window size is 24x24 pixels, and the overlapping ratio is 87.5%. Then the step size is set to one-eighth of the window size to achieve 87.5% overlap.

[0052] Preferably, the specific steps for constructing the circular depression shape template in S4 are as follows:

[0053] Step 1: Establish a spatial coordinate system and generate a grid

[0054] Define the coordinate range:

[0055] According to the theoretically maximum circular depression range, set the boundary of the coordinate system; on this basis, establish a local and accurate coordinate system to describe the shape of the circular depression:

[0056] Coordinate range:

[0057]

[0058] Generate a high-precision grid: Use a linear spacing function to generate a high-precision grid within the coordinate system; [template_x,template_y]=meshgrid(linspace(-r,r,template_size)). The resolution of the grid should be determined according to the analysis requirements and data accuracy to ensure that the subtle shape changes of the circular depression can be captured;

[0059] Step 2: Calculate the distance field

[0060] Calculate the Euclidean distance:

[0061] On the established grid, calculate the Euclidean distance from each grid point to the center of the circular depression;

[0062] Store the calculated distance values on the corresponding grid points to form a distance field;

[0063] The calculation formula of the Euclidean distance is as follows:

[0064]

[0065] Step 3: Model Construction and Optimization

[0066] The model is constructed using a piecewise function method:

[0067] According to the distance field, a piecewise function method is used to construct the morphological template of the circular depression; the piecewise function defines different morphological features according to different intervals of distance, including depth and curvature;

[0068]

[0069] The oblateness parameter is introduced to adjust the depth distribution, enabling the model to more flexibly adapt to circular depressions of different morphologies; at the same time, a smooth transition zone is used to handle the boundary, avoiding sudden changes at the boundary and achieving adaptive curvature adjustment;

[0070] Optimize boundary handling:

[0071] During the model construction process, special attention is paid to the handling of the boundary. The smooth transition of the boundary is achieved by using a smooth transition characteristic function, avoiding unnatural sudden changes at the boundary;

[0072] Conduct necessary verification and adjustment on the constructed model:

[0073] The formula of the characteristic function is as follows:

[0074]

[0075] Preferably, the specific steps of implementing the multi-feature matching algorithm in S5 are as follows:

[0076] Step 1: Calculate the morphological similarity

[0077] The Pearson correlation coefficient is used to calculate the matching degree between the template and the actual terrain:

[0078]

[0079] Among them, the complete calculation formula of the Pearson correlation coefficient is:

[0080]

[0081] Among them, xi, yi: the observed values of the two variables; The average value of the two variables; n: the sample size;

[0082] Step 2: Verify the local minimum

[0083] In the actual karst landform, the center point of the circular depression usually shows the lowest point in the local area; first, the δ neighborhood of the point x0 is defined theoretically:

[0084]

[0085] There exists a neighborhood \(N_{\delta}(x^{*})\) of the point \(x^{*}\) such that:

[0086]

[0087] Then the point \(x^{*}\) is called a local minimum point of the function \(f(x)\);

[0088] Based on this theoretical basis, an improved local minimum determination algorithm is adopted. This algorithm not only considers the minimum characteristics but also introduces mean comparison to improve the reliability of the recognition result:

[0089]

[0090] Step 3: Area-depth analysis

[0091] Depth calculation:

[0092] depth = max(window) - min(window)

[0093] Contour threshold setting: Take the contour line 20% above the lowest point as the cut surface

[0094] contour_level = min(window) + depth·0.2

[0095] Area calculation: Use the contour line method to calculate the area of the depression area, which can more accurately reflect the shape of large circular depressions and is insensitive to noise and local undulations

[0096]

[0097] Area-depth ratio:

[0098]

[0099] where: Setting the parameter to 0.1 is to avoid division by zero

[0100] Step 4: Comprehensive scoring mechanism

[0101] The final matching score is obtained through the weighted combination of multiple feature indicators:

[0102]

[0103] where:

[0104]

[0105] Within the experimental sample area, the area-depth ratio limits the selection range of circular depression candidate points.

[0106] Preferably, the specific steps for result verification and visualization enhancement in S6 are as follows:

[0107] Step 1: 3D reconstruction and fitting analysis of the ideal model

[0108] Data extraction and 3D reconstruction

[0109] Extract data of the extended area from the actual terrain data;

[0110] Use digital elevation model technology for 3D reconstruction to generate a 3D representation of the actual terrain;

[0111] Construct an ideal depression fitting model

[0112] Construct an ideal depression model based on the recognition parameters; among them, the depression model formula is as follows:

[0113]

[0114] Where: base_elevation: elevation of the local reference plane; local_depth: recognized depression depth; local_radius: equivalent radius calculated based on the area; offset: model offset, default is 280 meters;

[0115] Perform fitting comparison analysis between the ideal model and the actual 3D terrain to evaluate the difference between the actual terrain and the ideal model;

[0116] Visualization enhancement and candidate point selection

[0117] Perform visualization enhancement on the 3D model, including applying smooth shading, adding lighting effects, setting transparency effects, and adding coordinate axes and title annotations,

[0118] Randomly select 4 candidate points as the analysis objects;

[0119] Step 2: 2D contour analysis

[0120] Local area extraction and data smoothing

[0121] Extract the local area from the overall terrain data to determine the analysis scope;

[0122] Calculation of the expansion coefficient: extension_factor = 2 · window_size

[0123] Determination of the extraction range:

[0124]

[0125] Use the Gaussian kernel function to smooth the extracted data to reduce noise interference. The Gaussian kernel function calculation formula:

[0126]

[0127] where σ = 1 and the size of the kernel matrix is 5×5;

[0128] Calculation of the smoothed elevation value:

[0129]

[0130] Calculation of the contour interval:

[0131]

[0132] Contour elevation value sequence:

[0133] contour_levels = {h i |h i = min(DEM local ) + i·interval, i ∈ [0, 25]}

[0134] Calculation of the boundary circle coordinates:

[0135]

[0136] Radius calculation and model verification

[0137] Calculate the radius of the circular depression model from the area information in the contour map;

[0138] Compare the calculated radius with the equivalent radius in the ideal model to verify the accuracy of the model;

[0139] Step 3: Grayscale image analysis and high-resolution interpolation

[0140] High-resolution interpolation processing

[0141] Perform high-resolution interpolation processing on the original terrain data using a bicubic interpolation polynomial to improve the accuracy and resolution of the data;

[0142] Bicubic interpolation polynomial:

[0143]

[0144] Determine the size of the interpolation grid to ensure the accuracy and reliability of the interpolation result;

[0145] Interpolation grid size formula:

[0146] grid_size new = 0.5·grid_size original

[0147] Dense Contour Generation and Smoothing Optimization

[0148] The optimization formula for contour interval is as follows:

[0149]

[0150] Smoothing coefficient:

[0151]

[0152] Where d is the distance between adjacent points

[0153] Calculate and generate dense contours based on the interpolated elevation values;

[0154] Optimize the contour interval and smoothing coefficient to reduce the fluctuations and overlaps of the contours;

[0155] Gray-scale Image Visualization Enhancement and Annotation

[0156] Convert the contour map into a gray-scale image for further visualization analysis;

[0157] Normalize the gray-scale image to improve the contrast and clarity of the image;

[0158] Set the annotation interval of the contours to accurately annotate the elevation values of the contours on the gray-scale image;

[0159] Gray-scale value normalization:

[0160]

[0161] Contour annotation interval:

[0162]

[0163] Preferably, the parameter setting and threshold selection in S7 refer to determining key parameters during the parameter setting and threshold selection process: in terms of template parameters, the radius is set to 350, the flatness ratio is set to 0.6 to match the morphological features, and the grid subdivision is selected as 100 to achieve the balance of accuracy. As for the algorithm parameters, the sliding step length needs to be set according to specific situations, the morphological similarity threshold is set to not less than 0.7 to ensure a high degree of morphological matching, and the area-depth ratio constraint is used as another important algorithm constraint condition, and its specific value needs to be carefully adjusted according to actual situations;

[0164] Area-depth ratio constraint formula:

[0165]

[0166] The beneficial effects of the present invention are as follows:

[0167] Through a standardized data processing process and automatic coordinate conversion technology, the present invention realizes the full-process automation from data input to recognition output. It has high recognition accuracy, thanks to an advanced discrimination method that combines multiple feature indicators and precise spatial positioning. To ensure the reliability of the results, a unified quantitative evaluation standard is established, and a complete coordinate conversion and data processing mechanism is built-in. In addition, it shows wide applicability and can flexibly adapt to DEM data with different geographical reference systems and resolutions. According to actual needs, the standard model can be designed to be multi-scale, effectively supporting the precise recognition of various terrain features such as mountains and plains, demonstrating strong application potential. Brief Description of the Drawings

[0168] Figure 1 It is a flow chart of the method for identifying large karst dolines based on DEM data of the present invention;

[0169] Figure 2 It is a schematic diagram of the three-dimensional hemisphere model of the present invention;

[0170] Figure 3 It is an analysis chart of the area-depth ratio of karst feature regions of the present invention;

[0171] Figure 4 It is a distribution map of karst feature regions with a matching score ≥ 70% of the present invention;

[0172] Figure 5 It is a three-dimensional comparison chart of karst feature regions of the present invention;

[0173] Figure 6 It is a two-dimensional map for the analysis of doline terrain of the present invention. Detailed Embodiment

[0174] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0175] As Figures 1 to 6 shown, the embodiment of the present invention provides a method for identifying large karst dolines based on DEM data. The specific steps of the identification method are as follows:

[0176] S1: DEM data preprocessing: Read the DEM data, clean and process it, and convert the coordinate system to ensure that the data accurately corresponds to the ground surface;

[0177] S2: Standard template construction: Define the hemisphere parameters, convert the coordinates, perform template matching, and determine the optimal model parameters;

[0178] S3: Design a sliding window mechanism: Design the window and overlapping areas to ensure the complete recognition of circular depressions and balance efficiency and accuracy;

[0179] S4: Construct a morphological template for circular depressions: Establish a coordinate grid, calculate the distance field, and construct an optimized morphological model of circular depressions;

[0180] S5: Implement a multi-feature matching algorithm: Implement a multi-feature matching algorithm, comprehensively analyze morphology, local minimum, and area-depth, and obtain a matching score;

[0181] S6: Result verification and visualization enhancement: Conduct 3D reconstruction and fitting analysis, 2D contour analysis, and grayscale image analysis to comprehensively evaluate candidate circular depressions;

[0182] S7: Parameter setting and threshold selection: The parameters and thresholds need to be set precisely, including a radius of 350 and a flatness ratio of 0.6, to ensure the matching of morphology and accuracy.

[0183] Among them, the specific steps of DEM data preprocessing in S1 are as follows:

[0184] Step 1: Data reading and verification

[0185] Read DEM data in.tif format

[0186] Use the readgeoraster function (assuming it is a function in a certain geographic information system or Python library, such as GDAL, rasterio, or a custom function) to read DEM data in.tif format. Ensure that the required libraries are installed and correctly configured;

[0187] Verify the data format and georeference information

[0188] Verify whether the read data is in a two-dimensional matrix format and check whether the georeference information (such as coordinate system and projection information) is correct;

[0189] Step 2: Data cleaning

[0190] Identify and process invalid values

[0191] Use functions in the NumPy library to identify and process NaN values and outliers, and replace NaN values with the average of adjacent pixels or perform other reasonable interpolation processing;

[0192] Convert the data type and smooth the data

[0193] Ensure that the data type is double precision (double), and perform data smoothing (such as Gaussian filtering);

[0194] Step 3: Coordinate system conversion

[0195] Based on the average radius of the Earth being approximately 6,371,000 meters and combined with the resolution of the Digital Elevation Model (DEM), that is, each pixel represents an actual surface distance of 30 meters, the geographical coordinates (longitude and latitude) are converted into actual distance coordinates on the surface. The conversion process takes into account the curvature of the Earth and the spatial resolution of the DEM data, thereby ensuring that each pixel can accurately correspond to its actual surface distance in the geographical space.

[0196] Among them, the specific steps for constructing the standard template in S2 are as follows:

[0197] Step 1: Definition of basic parameters and model establishment

[0198] Define the basic parameters of the hemisphere: Through data analysis, the optimal parameters of the hemisphere are determined. Radius (r): 350 meters (optimal value determined based on statistical analysis); Height (H): 210 meters (considering actual terrain features); Flattening ratio (f): Number of grid subdivisions (numDiv): 100 (balancing calculation efficiency and accuracy); The three-dimensional hemisphere model is as Figure 2 shown;

[0199] Define the coordinate system parameters; Azimuth angle (θ): θ ∈ [0, 2π]; Polar angle (φ):

[0200] Step 2: Three-dimensional coordinate transformation and template matching

[0201] Perform three-dimensional coordinate transformation:

[0202] The standard template describes the hemisphere surface using parametric equations, and the flattening ratio parameter is introduced to achieve precise fitting to the actual terrain:

[0203]

[0204] Among them, X, Y, Z: Positions in the spatial rectangular coordinate system; r: Hemisphere radius; f: Flattening ratio; H: Datum plane height;

[0205] Apply the three-dimensional coordinate transformation results for template matching:

[0206] After completing the three-dimensional coordinate transformation, the obtained three-dimensional coordinates are used for template matching; mainly, the actual terrain data is accurately compared with the hemisphere model. By calculating errors, performing correlation analysis, or using machine learning algorithm methods, the parameter configuration of the hemisphere model that best matches the actual terrain is found. Finally, the goal is to determine the hemisphere model parameters that can best reflect the characteristics of the actual terrain, thereby laying a solid foundation for subsequent terrain analysis and modeling.

[0207] Among them, the specific steps for designing the sliding window mechanism in S3 are as follows:

[0208] Step 1: Window Size Design

[0209] Determine the expected radius of the circular depression and the resolution of the DEM data:

[0210] According to the analysis requirements, clarify the radius of the expected circular depression (such as 350 meters);

[0211] Obtain the resolution of the DEM (Digital Elevation Model) data, that is, the actual distance represented by each pixel (such as 30 meters / pixel);

[0212] Calculate the window size:

[0213] Use the window size calculation formula, divide the diameter of the expected circular depression by the DEM data resolution, and round up to ensure that the window can completely cover the expected circular depression shape; for example, 2 * 350 meters / 30 meters / pixel ≈ 23.33, rounded up to 24, so the window size should be designed as 24x24 pixels.

[0214] The specific formula of the window size calculation formula is as follows:

[0215]

[0216] Step 2: Overlap Region Design

[0217] Set the overlap ratio:

[0218] To avoid recognition deviation caused by the circular depression being located at the window boundary, introduce an overlapping scanning mechanism; its formula is as follows:

[0219]

[0220] Determine the overlap ratio between adjacent windows, and the overlap ratio is 87.5%; this means that each window will have a large overlap with adjacent windows to ensure the continuity of boundary features.

[0221] Calculate the step size and set the scanning mechanism:

[0222] Calculate the step size according to the overlap ratio. When the window size is 24x24 pixels and the overlap ratio is 87.5%, the step size can be set to one-eighth of the window size (i.e., 3 pixels) to achieve 87.5% overlap;

[0223] Set the scanning mechanism, move the window according to the calculated step size, and scan the same area multiple times to improve the accuracy and reliability of recognition. At the same time, pay attention to balancing the coverage integrity and computational efficiency to ensure the complete recognition of the target features without excessive increasing the computational burden.

[0224] Among them, the specific steps for constructing the circular depression form template in S4 are as follows:

[0225] Step 1: Establishment of spatial coordinate system and generation of grid

[0226] Define the coordinate range:

[0227] According to the theoretically maximum circular depression range, set the boundaries of the coordinate system; ensure that this range can completely cover the expected circular depression form and reserve a certain margin for subsequent normalization processing and distance calculation. On this basis, establish a local and accurate coordinate system to describe the form of the circular depression:

[0228] Coordinate range:

[0229]

[0230] Generate a high-precision grid: Use the linear spacing function to generate a high-precision grid within the coordinate system; [template_x, template_y]=meshgrid(linspace(-r, r, template_size)). The resolution of the grid should be determined according to the analysis requirements and data accuracy to ensure that the subtle form changes of the circular depression can be captured.

[0231] Step 2: Calculation of distance field

[0232] Calculate the Euclidean distance:

[0233] On the established grid, calculate the Euclidean distance from each grid point to the center of the circular depression; implement it through programming, and use loops or matrix operations to traverse the grid points and calculate their distances to the center.

[0234] Store the calculated distance values on the corresponding grid points to form a distance field;

[0235] The calculation formula of the Euclidean distance is as follows:

[0236]

[0237] Step 3: Model construction and optimization

[0238] Adopt the piecewise function method to construct the model:

[0239] According to the distance field, adopt the piecewise function method to construct the form template of the circular depression; the piecewise function defines different form characteristics according to different distance intervals, including depth and curvature;

[0240]

[0241] The oblateness parameter is introduced to adjust the depth distribution, enabling the model to more flexibly adapt to circular depressions of different shapes; meanwhile, a smooth transition zone is used to handle the boundary, avoiding sudden changes at the boundary and achieving adaptive curvature adjustment;

[0242] Optimize boundary handling:

[0243] During the model construction process, special attention is paid to the handling of the boundary. The smooth transition of the boundary is achieved by using a smooth transition characteristic function, avoiding unnatural sudden changes at the boundary;

[0244] Necessary verification and adjustment are carried out on the constructed model to ensure that it can accurately reflect the gradual change characteristics of the actual terrain and provide a reliable basis for subsequent analysis and modeling;

[0245] The formula of the characteristic function is as follows:

[0246]

[0247] Among them, the specific steps of implementing the multi-feature matching algorithm in S5 are as follows:

[0248] Step 1: Calculate the morphological similarity

[0249] The Pearson correlation coefficient is used to calculate the matching degree between the template and the actual terrain:

[0250]

[0251] Among them, the complete calculation formula of the Pearson correlation coefficient is:

[0252]

[0253] Among them, xi, yi: the observed values of two variables; the average values of the two variables; n: the number of samples;

[0254] The strength determination criteria of the correlation coefficient:

[0255] Absolute value of correlation coefficient Degree of correlation Practical significance 0.8-1.0 Strong correlation There is a strong linear relationship between variables 0.5-0.8 Moderate correlation There is an obvious linear relationship between variables 0.3-0.5 Weak correlation There is a weak linear relationship between variables 0.0-0.3 Very weak correlation There is hardly any linear relationship between variables

[0256] Step 2: Verify the local minimum

[0257] In the actual karst landform, the center point of the circular depression usually shows the lowest point in the local area; first, the δ neighborhood of point x0 is defined theoretically:

[0258]

[0259] There exists a certain neighborhood N_δ(x*) of point x*, such that:

[0260]

[0261] Then the point x* is called the local minimum point of the function f(x);

[0262] Based on this theoretical basis, an improved local minimum determination algorithm is adopted. This algorithm not only considers the minimum characteristics but also introduces mean comparison to improve the reliability of the recognition result:

[0263]

[0264] Step 3: Area-depth analysis

[0265] Depth calculation:

[0266] depth = max(window) - min(window)

[0267] Contour threshold setting: Take the contour line 20% above the lowest point as the intercept plane

[0268] contour_level = min(window) + depth · 0.2

[0269] Area calculation: Use the contour line method to calculate the area of the depression area, which can more accurately reflect the shape of large circular depressions and is not sensitive to noise and local undulations

[0270]

[0271] Area-depth ratio:

[0272]

[0273] where: Setting the parameter to 0.1 is to avoid division by zero

[0274] Step 4: Comprehensive scoring mechanism

[0275] The final matching score is obtained through the weighted combination of multiple feature indicators:

[0276]

[0277] where:

[0278]

[0279] Within the experimental sample area, the area-depth ratio limits the selection range of circular depression candidate points, as Figure 3 shown; A total of 47 circular depression candidate points with a matching score ≥ 70% appeared, as Figure 4 shown.

[0280] Among them, the specific steps of result verification and visualization enhancement in S6 are as follows:

[0281] Step 1: 3D Reconstruction and Fitting Analysis with Ideal Model

[0282] Data Extraction and 3D Reconstruction

[0283] Extract data of the extended area from the actual terrain data.

[0284] Use digital elevation model (DEM) technology for 3D reconstruction to generate a 3D representation of the actual terrain.

[0285] Construct an ideal depression fitting model

[0286] Based on the identification parameters (including the local datum elevation base_elevation, the identified depression depth local_depth, the equivalent radius local_radius calculated based on the area, and the default model offset offset of 280 meters), construct an ideal depression model; where the depression model formula is as follows:

[0287]

[0288] Where: base_elevation: local datum elevation; local_depth: identified depression depth; local_radius: equivalent radius calculated based on the area; offset: model offset, default is 280 meters.

[0289] Perform fitting and comparative analysis between the ideal model and the actual 3D terrain to evaluate the differences between the actual terrain and the ideal model.

[0290] Visualization Enhancement and Candidate Point Selection

[0291] Perform visualization enhancement on the 3D model, including applying smooth shading, adding lighting effects, setting transparency effects, and adding coordinate axes and title annotations.

[0292] Randomly select 4 candidate points as the analysis objects, as Figure 5 shown, for further detailed analysis.

[0293] Step 2: 2D Contour Analysis

[0294] Local Area Extraction and Data Smoothing Processing

[0295] Extract the local area from the overall terrain data to determine the analysis scope;

[0296] Calculation of the extension coefficient: extension_factor = 2 · window_size

[0297] Determination of the extraction range:

[0298]

[0299] The extracted data is smoothed using a Gaussian kernel function (kernel matrix size is 5×5) to reduce noise interference. The calculation formula of the Gaussian kernel function is:

[0300]

[0301] where σ = 1 and the kernel matrix size is 5×5;

[0302] Calculation of the elevation value after smoothing:

[0303]

[0304] Calculation of the contour interval:

[0305]

[0306] Contour elevation value sequence:

[0307] contour_levels = {h i |h i = min(DEM local ) + i·interval, i ∈ [0, 25]}

[0308] Calculation of the boundary circle coordinates:

[0309]

[0310] Radius calculation and model verification

[0311] Calculate the radius of the circular depression model from the area information in the contour map;

[0312] Compare the calculated radius with the equivalent radius in the ideal model to verify the accuracy of the model;

[0313] Step 3: Grayscale image analysis and high - resolution interpolation

[0314] High - resolution interpolation processing

[0315] The original terrain data is processed by high - resolution interpolation using a bicubic interpolation polynomial to improve the accuracy and resolution of the data;

[0316] Bicubic interpolation polynomial:

[0317]

[0318] Determine the size of the interpolation grid to ensure the accuracy and reliability of the interpolation result;

[0319] Interpolation grid size formula:

[0320] grid_size new = 0.5·grid_size original

[0321] Dense Contour Generation and Smoothing Optimization

[0322] The contour interval optimization formula is as follows:

[0323]

[0324] Smoothing coefficient:

[0325]

[0326] where d is the distance between adjacent points

[0327] Calculate and generate dense contour lines based on the interpolated elevation values;

[0328] Optimize the contour interval and smoothing coefficient to reduce the fluctuations and overlaps of the contour lines; where the smoothing coefficient d is a function of the distance between adjacent points.

[0329] Gray-scale Map Visualization Enhancement and Annotation

[0330] Convert the contour map into a gray-scale map for further visual analysis.

[0331] Normalize the gray-scale map to improve the contrast and clarity of the image.

[0332] Set the annotation interval of the contour lines to accurately annotate the elevation values of the contour lines on the gray-scale map.

[0333] Gray-scale value normalization:

[0334]

[0335] Contour line annotation interval:

[0336]

[0337] where the parameter setting and threshold selection in S7 refer to determining the key parameters during the parameter setting and threshold selection process: in terms of template parameters, the radius is set to 350 (this value is obtained based on statistical analysis), the flatness ratio is set to 0.6 to match the morphological features, and the grid subdivision is selected as 100 to achieve the balance of accuracy. As for the algorithm parameters, the sliding step length needs to be set according to specific situations, the morphological similarity threshold is set to not less than 0.7 to ensure a high morphological match, and the area-depth ratio constraint is used as another important algorithm constraint condition, and its specific value needs to be carefully adjusted according to the actual situation;

[0338] Area-depth ratio constraint formula:

[0339]

[0340] It should be noted that in this article, relational terms such as first and second are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the term "comprising", "including" or any other variant thereof is intended to cover non-exclusive inclusion, so that a process, method, article or device comprising a series of elements not only includes those elements, but also includes other elements not expressly listed, or further includes elements inherent to such process, method, article or device.

[0341] Although the embodiments of the present invention have been shown and described, it will be understood by those of ordinary skill in the art that various changes, modifications, substitutions and variations can be made to these embodiments without departing from the principles and spirit of the present invention, and the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for identifying large circular depressions in karst based on DEM data, characterized by: The specific steps of this identification method are as follows: S1: DEM data preprocessing: read DEM data, clean it, convert the coordinate system, and ensure that the data accurately corresponds to the surface; S2: Standard template construction: define hemisphere parameters, transform coordinates, template matching, and determine optimal model parameters; S3: Design sliding window mechanism: Design windows and overlapping areas to ensure complete identification of circular depressions and balance efficiency and accuracy; S4: Construct circular depression morphology template: establish coordinate grid, calculate distance field, and construct optimized circular depression morphology model; S5: Implement multi-feature matching algorithm: Implement multi-feature matching algorithm, integrate morphology, local minimum, area-depth analysis, and obtain matching scores; S6: Result verification and visualization enhancement: 3D reconstruction and fitting analysis, 2D contour analysis, grayscale image analysis, and comprehensive evaluation of candidate circular depressions; S7: Parameter setting and threshold selection: Parameter and threshold settings must be accurate, including a radius of 350 and a flattening rate of 0.6, to ensure that the shape and accuracy match.

2. The method for identifying large circular depressions in karst based on DEM data according to claim 1, characterized in that: The specific steps of DEM data preprocessing in S1 are as follows: Step 1: Data reading and verification Read DEM data in .tif format Use the readgeoraster function to read DEM data in .tif format; Validate data format and georeferencing information Verify that the read data is in 2D matrix format and check if the geo-reference information is correct; Step 2: Data cleaning Identifying and handling invalid values Use functions in the NumPy library to identify and handle NaN values ​​and outliers, replacing NaN values ​​with the average value of neighboring pixels or performing other reasonable interpolation processing; Convert data types and smooth them Ensure that the data type is double precision and perform data smoothing; Step 3: Coordinate system conversion Based on the average radius of the earth of approximately 6,371,000 meters and the resolution of the digital elevation model, that is, each pixel represents an actual surface distance of 30 meters, the geographic coordinates are converted into actual distance coordinates on the surface. The conversion process takes into account the curvature of the earth and the spatial resolution of the DEM data.

3. The method for identifying large circular depressions in karst based on DEM data according to claim 1, characterized in that: The specific steps of constructing the standard template in S2 are: Step 1: Basic parameter definition and model establishment Define the basic parameters of the hemisphere: Through data analysis, the optimal hemisphere parameters were determined, with a radius of 350 meters; Height: 210 meters; Flatness: Grid subdivision number: 100; Define the coordinate system parameters; azimuth (θ): θ∈[0,2π]; polar angle (φ): Step 2: 3D coordinate transformation and template matching Perform three-dimensional coordinate transformation: The standard template uses a parametric equation to describe the surface of a hemisphere and introduces a flattening parameter to achieve an accurate fit to the actual terrain: Where, X, Y, Z: position in the spatial rectangular coordinate system; r: radius of the hemisphere; f: flattening rate; H: height of reference plane; Apply the three-dimensional coordinate transformation results for template matching: After completing the three-dimensional coordinate conversion, the obtained three-dimensional coordinates are used for template matching; mainly, the actual terrain data is accurately compared with the hemispherical model, and the hemispherical model parameter configuration that best matches the actual terrain is found by calculating errors, implementing correlation analysis, or using machine learning algorithms. The ultimate goal is to determine the hemispherical model parameters that best reflect the actual terrain characteristics, thereby laying a solid foundation for subsequent terrain analysis and modeling.

4. The method for identifying large circular depressions in karst based on DEM data according to claim 1, characterized in that: The specific steps of designing the sliding window mechanism in S3 are as follows: Step 1: Window size design Determine the expected circular depression radius and DEM data resolution: According to the analysis requirements, the radius of the expected circular depression is determined; Get the resolution of DEM data, that is, the actual distance represented by each pixel; Calculate the window size: Use the window size calculation formula to divide the expected circular depression radius by the DEM data resolution and round up to ensure that the window can completely cover the expected circular depression shape; The specific formula for calculating the window size is as follows: Step 2: Overlapping area design Set the overlap ratio: In order to avoid the recognition deviation caused by the circular depression at the window boundary, an overlapping scanning mechanism is introduced; its formula is as follows: Determine the overlap ratio between adjacent windows, which is 87.5%; Calculate the step size and set up the scanning mechanism: The step size is calculated based on the overlap ratio. If the window size is 24x24 pixels and the overlap ratio is 87.5%, the step size is set to one eighth of the window size to achieve 87.5% overlap.

5. The method for identifying large circular depressions in karst based on DEM data according to claim 1, characterized in that: The specific steps of constructing the circular depression morphology template in S4 are as follows: Step 1: Establishing the spatial coordinate system and generating the mesh Define the coordinate range: According to the theoretical maximum circular depression range, the boundary of the coordinate system is set; on this basis, a local and accurate coordinate system is established to describe the shape of the circular depression: Coordinate range: Generate a high-precision grid: Use the linear interval function to generate a high-precision grid in the coordinate system; [template_x, template_y] = meshgrid (linspace (-r, r, template_size)) The resolution of the grid should be determined according to the analysis requirements and data accuracy to ensure that the subtle morphological changes of the circular depression can be captured; Step 2: Distance Field Calculation Calculate the Euclidean distance: On the established grid, calculate the Euclidean distance from each grid point to the center of the circular depression; The calculated distance values ​​are stored at corresponding grid points to form a distance field; The calculation formula of Euclidean distance is as follows: Step 3: Model building and optimization The model is constructed using the piecewise function method: According to the distance field, a piecewise function method is used to construct the morphological template of the circular depression; the piecewise function defines different morphological features according to different distance intervals, including depth and curvature; The flattening parameter is introduced to adjust the depth distribution, so that the model can more flexibly adapt to circular depressions of different shapes; at the same time, a smooth transition zone is used to process the boundary to avoid sudden changes at the boundary and realize adaptive curvature adjustment; Optimized border processing: In the process of model construction, special attention is paid to the treatment of boundaries. By adopting the smooth transition indicator function, the smooth transition of boundaries is achieved to avoid unnatural mutations at the boundaries. Perform necessary verification and adjustments on the constructed model: The schematic function formula is as follows:

6. The method for identifying large circular depressions in karst based on DEM data according to claim 1, characterized in that: The specific steps of implementing the multi-feature matching algorithm in S5 are as follows: Step 1: Morphological similarity calculation The Pearson correlation coefficient is used to calculate the matching degree between the template and the actual terrain: Among them, the complete calculation formula of Pearson correlation coefficient is: Among them, xi, yi: observed values ​​of two variables; The average of two variables; n: sample size; Step 2: Local Minimum Verification In actual karst landforms, the center point of a circular depression is usually the lowest point in the local area. First, the δ neighborhood of point x0 is theoretically defined: There exists a neighborhood N_δ(x*) of a point x* such that: Then the point x* is called the local minimum point of the function f(x); Based on this theoretical foundation, an improved local minimum determination algorithm is adopted, which not only considers the minimum value feature, but also introduces mean comparison to improve the reliability of the recognition result: Step 3: Area-Depth Analysis Depth calculation: depth=max(window)-min(window) Contour threshold setting: Take the contour line 20% above the lowest point as the intercept surface contour_level=min(window)+depth·0.2 Area calculation: The area of ​​the depression is calculated using the contour method, which more accurately reflects the morphology of large circular depressions and is insensitive to noise and local fluctuations. Area-depth ratio: The parameter is set to 0.1 to avoid division by 0. Step 4: Comprehensive scoring mechanism The final matching score is obtained through a weighted combination of multiple feature indicators: in: In the experimental sample area, the area-depth ratio limits the range of candidate points for circular depressions.

7. The method for identifying large circular depressions in karst based on DEM data according to claim 1, characterized in that: The specific steps of result verification and visualization enhancement in S6 are as follows: Step 1: 3D reconstruction and ideal model fitting analysis Data extraction and 3D reconstruction Extracting data of the extended area from actual terrain data; 3D reconstruction using digital elevation model technology to produce a 3D representation of the actual terrain; Constructing an ideal concave fitting model Construct an ideal concave model based on the identification parameters; The concave model formula is as follows: Among them: base_elevation: local base elevation; local_depth: identified depression depth; local_radius: equivalent radius based on area calculation; offset: model offset, the default is 280 meters; The ideal model is fitted and compared with the actual three-dimensional terrain to evaluate the difference between the actual terrain and the ideal model; Visualization enhancement and candidate point selection Enhance the visualization of 3D models by applying smooth shading, adding lighting effects, setting transparency, and adding axis and title labels. Four candidate points were randomly selected as analysis objects; Step 2: 2D Contour Analysis Local area extraction and data smoothing Extract local areas from the overall terrain data to determine the scope of analysis; Extension factor calculation: extension_factor = 2·window_size Extraction range determination: Use the Gaussian kernel function to smooth the extracted data to reduce noise interference. The Gaussian kernel function calculation formula is: Where σ = 1, the kernel matrix size is 5 × 5; Calculation of elevation value after smoothing: Contour interval calculation: Sequence of contour elevation values: contour_levels={h i |h i =min(DEM local )+i·interval,i∈[0,25]} Bounding circle coordinate calculation: Radius calculation and model verification Calculate the radius of the circular depression model from the area information in the contour map; The calculated radius was compared with the equivalent radius in the ideal model to verify the accuracy of the model; Step 3: Grayscale image analysis and high-resolution interpolation High-resolution interpolation The original terrain data is interpolated with high resolution using bicubic interpolation polynomials to improve the accuracy and resolution of the data. Bicubic interpolation polynomial: Determine the size of the interpolation grid to ensure the accuracy and reliability of the interpolation results; Interpolation grid size formula: grid_size new =0.5·grid_size original Dense contour generation and smoothing optimization The contour interval optimization formula is as follows: Smoothing factor: Where d is the distance between adjacent points Calculate and generate dense contour lines based on interpolated elevation values; Optimize the interval and smoothing coefficient of contour lines to reduce the fluctuation and overlap of contour lines; Grayscale image visualization enhancement and annotation Convert contour maps to grayscale for further visualization analysis; Normalize the grayscale image to improve the contrast and clarity of the image; Set the labeling interval of the contour lines so that the elevation values ​​of the contour lines can be accurately marked on the grayscale map; Grayscale value normalization: Contour labeling interval:

8. The method for identifying large circular depressions in karst based on DEM data according to claim 1, characterized in that: The parameter setting and threshold selection in S7 refer to determining key parameters in the parameter setting and threshold selection process: in terms of template parameters, the radius is set to 350, the flattening rate is set to 0.6 to match the morphological features, and the grid subdivision is selected to 100 to achieve a balance in accuracy. As for the algorithm parameters, the sliding step size needs to be set according to the specific situation, the morphological similarity threshold is set to no less than 0.7 to ensure a high degree of morphological matching, and the area-depth ratio constraint is another important algorithm constraint condition, and its specific value needs to be carefully adjusted according to the actual situation; Area-to-depth ratio constraint formula:

Images