A method for estimating sand dune density based on DEM

Abstract

The present invention discloses a method for estimating dune density based on DEM, which relates to the field of information technology. The method includes: preprocessing the DEM data; extracting dune vertices by using an eight-neighborhood Laplacian operator based on the preprocessed DEM data; conducting connectivity analysis of the dune area and preliminary estimation of dune density based on the dune vertices; further calculating the preliminary estimation result of dune density by combining the four-direction topographic spectrum decomposition method; after estimating the dune density by the four-direction topographic spectrum decomposition method, adopting a machine learning-based method for correction and optimization; estimating and fusing the dune density at multiple scales by performing multi-scale decomposition on the dune density; comprehensively analyzing the multi-scale dune density results and making a visual display. Through higher-precision vertex extraction, more reasonable connectivity analysis, and the introduction of machine learning, the present invention can provide more accurate and comprehensive results in the estimation of dune density at multiple scales and in multiple directions.

Description

Technical Field

[0001] The invention relates to the technical field of information technology, and in particular to a dune density estimation method based on DEM. Background Art

[0002] Northwest my country is one of the two main sources of sandstorms in East Asia. The region has a dry climate and little rainfall. In spring and summer, strong winds often cause dust storms, which are sudden and local. At present, reliable sandstorm forecasts can guide various industries to take countermeasures in advance and reduce the harm caused by sandstorms. The forecast of sandstorms mainly relies on numerical models, but the physical scheme of the model does not consider the impact of small-scale terrain fluctuations in the desert on wind speed, resulting in a large error between the simulated wind speed in the desert area and the actual observed value. Therefore, it is very important to finely characterize the characteristics of the desert underlying surface in the numerical model. Traditional dune density estimation methods mainly rely on manual statistics, which is to manually count the dune density by establishing sample plots in specific areas and continuously observing them, or to obtain dune images in fixed areas by drones and other aerial cameras, and then rely on manual visual interpretation to count the dune density. However, the environment in desert areas is harsh. The method of establishing sample plots to count the density of sand dunes is time-consuming and labor-intensive, and can only count the density of sand dunes in a certain area. It is impossible to obtain the density of sand dunes in uninhabited areas such as the deep desert. Therefore, it is impossible to count the density of sand dunes in a large area. Estimating the density of sand dunes using aerial photos can only rely on human visual interpretation, which requires a high level of professionalism and prior knowledge.

[0003] With the advancement of high-resolution satellite remote sensing technology, the dune density of large dunes can be obtained by using satellite remote sensing images combined with geographic information software. The use of satellite images for dune morphology research is limited by image resolution, and can only obtain the density of large fixed dunes while ignoring the density of small and medium-sized dunes. With the improvement of DEM (digital elevation model) data resolution, it is possible to calculate the morphological parameters of small and medium-sized dunes using high-resolution DEM, but the dune density has not been estimated. Traditional dune density statistical methods are time-consuming and inefficient, and it is difficult to quickly estimate the dune density of large areas. Summary of the invention

[0004] In view of the above existing problems, the present invention is proposed.

[0005] Therefore, the present invention provides a dune density estimation method based on DEM to solve the problems of insufficient accuracy of dune density estimation and imperfect multi-scale analysis in the prior art.

[0006] In order to solve the above technical problems, the present invention provides the following technical solutions:

[0007] In a first aspect, an embodiment of the present invention provides a method for estimating sand dune density based on DEM, which comprises:

[0008] Based on the preprocessed DEM data, the dune vertices are extracted using the eight-domain Laplacian operator;

[0009] Connectivity analysis of dune areas and preliminary estimation of dune density based on dune vertices;

[0010] The preliminary estimation results of dune density were further calculated by combining the four-directional terrain spectrum decomposition method;

[0011] After estimating the dune density through the four-directional terrain spectrum decomposition method, a machine learning-based method was used for correction and optimization;

[0012] By decomposing the dune density at multiple scales, the dune density at multiple scales can be estimated and integrated;

[0013] Comprehensively analyze the multi-scale dune density results and make a visual display.

[0014] As a preferred solution of the DEM-based dune density estimation method described in the present invention, the preprocessing of DEM data includes reading and preliminary inspection of DEM data, outlier detection and correction, resolution standardization of DEM data, fusion of multi-source DEM data, terrain filtering and noise elimination.

[0015] As a preferred solution of the dune density estimation method based on DEM of the present invention, the dune vertices are extracted using the eight-domain Laplacian operator based on the preprocessed DEM data, including the following steps:

[0016] Apply the eight-domain Laplace operator to the preprocessed DEM data to calculate the second-order derivative value of each pixel. The expression is:

[0017]

[0018] Among them, h(i,j) represents the elevation value at position (i,j) in the DEM data;

[0019] This calculation will generate a numerical value for each pixel, representing its local elevation change characteristics;

[0020] Use Gaussian filter to smooth the candidate vertices, and filter the smoothed vertices to remove those vertices with insignificant height changes;

[0021] After smoothing and filtering, in order to avoid too dense vertex distribution, a minimum spacing d can be set min , if the distance between two vertices is less than d min , then keep the vertex with higher height and delete the other one. The distance calculation formula is:

[0022]

[0023] Among them, d(P i ,P j ) is point P i and point P j The Euclidean distance between two points in a two-dimensional plane, point P i The coordinates of (x i ,y i ), click P j The coordinates of (x j ,y j ), x i and i P i The coordinates, x j and j P j The coordinates of

[0024] In order to avoid too sparse vertex distribution, a maximum spacing d can be set max If the distance between vertices in a region is too large, then additional vertices are inserted in the region. The positions of the inserted vertices can be determined by the interpolation method. The interpolation formula is:

[0025]

[0026] Among them, P new is the position where the vertex is inserted;

[0027] Mark all dune vertices to generate a binary image, where the vertex position is marked as 1 and other areas are marked as 0;

[0028] The extracted dune vertex data are saved in a standard geographic information format to generate a vertex distribution map.

[0029] As a preferred solution of the dune density estimation method based on DEM described in the present invention, the connectivity analysis of the dune area and the preliminary estimation of the dune density based on the dune vertices include the following steps:

[0030] After extracting the dune vertices, a distance-based clustering algorithm is used to cluster the dune vertices;

[0031] By setting a radius ∈ and the minimum number of neighbors MinPts, the points containing at least MinPts vertices within the given radius are regarded as core points. The direct or indirect connections between core points form a cluster, which is expressed as,

[0032] Cluster(P)={P′|dist(P,P′)≤∈and|N ∈(P)|≥MinPts}

[0033] Among them, P and P' are two vertices, dist(P,P') represents the distance between the two vertices, and N ∈ (P) represents the set of vertices in the neighborhood with P as the center and radius ∈;

[0034] The parameters ∈ and MinPts need to be set according to the actual scale and density of the dunes. After clustering, each cluster represents a dune area;

[0035] For each cluster, the convex hull algorithm is used to identify the boundaries of the dunes;

[0036] After identifying each dune area, the connected component labeling algorithm was used to perform topological connectivity analysis between different dunes;

[0037] Based on the identified dune areas, the area of ​​each dune area is calculated by integration method, and the expression is:

[0038] Area(A)=∫∫ A dxdy

[0039] Where Area(A) is the area of ​​region A, A represents the two-dimensional projection of a dune region, dx and dy represent the small length increments in the x and y directions, respectively, and the integration range is the boundary of the region;

[0040] The dune density is calculated using area density to reflect the density of dunes in the study area;

[0041] Combining the connectivity analysis results and density calculation results of the dune area, a connectivity index was introduced to quantify the degree of connectivity of the dune group. The higher the connectivity index, the higher the degree of connectivity of the dunes.

[0042] As a preferred solution of the dune density estimation method based on DEM described in the present invention, the further calculation of the preliminary estimation result of the dune density by combining the four-directional terrain spectrum decomposition method includes the following steps:

[0043] Perform a two-dimensional fast Fourier transform on the DEM data to convert the terrain data in the spatial domain into the frequency domain. The calculation formula is:

[0044]

[0045] Among them, h(x,y) is the elevation value in the spatial domain, H(u,v) is the Fourier coefficient in the frequency domain, M and N are the number of rows and columns of DEM data respectively;

[0046] Using Gaussian bandpass filter, four directional filters are constructed, the expression is,

[0047]

[0048] Among them, G(u,v) is the value of the two-dimensional Gaussian function, u and v are the coordinates in the two-dimensional space, u 0 and v 0 is the center coordinate of the Gaussian function, (u 0 ,v 0 ) is the frequency center, θ is the direction angle of the filter, σ u and σ v is the standard deviation of the Gaussian function in the u direction and the v direction, which controls the bandwidth of the filter in the main direction and the vertical direction respectively;

[0049] Four directional filters are applied to the terrain data in the frequency domain to obtain the spectrum components in four directions. The result in each direction represents the main high-frequency components of the terrain in that direction. These components mainly correspond to the morphological characteristics of the sand dunes. Suppose the result after directional filtering is:

[0050] H dir (u,v)=H(u,v)·G dir (u,v)

[0051] Among them, H dir (u,v) is the Fourier coefficient in the two-dimensional frequency domain, G dir (u,v) represents a directional filter in a specific direction;

[0052] Perform inverse two-dimensional Fourier transform on the spectrum components in four directions respectively, and transform the filtering results in the frequency domain back to the spatial domain. The formula is:

[0053]

[0054] Among them, h dir (x, y) is the function value in the two-dimensional space domain;

[0055] Through inverse transformation, the spatial distribution characteristics of sand dunes in all directions were obtained;

[0056] For each direction, the high frequency component h dir (x,y) performs local extreme value detection, the formula is,

[0057] h max (x,y)={h dir (x,y)|h dir (x,y)>h N(x,y)}

[0058] Among them, h N(x,y) Represents the set of elevation values ​​within the neighborhood;

[0059] Calculate the local dune density in each direction using the formula,

[0060]

[0061] Among them, A is the neighborhood area, δ(h max (i,j)) is the indicator function, if h max If (i,j) has a local extreme value, the value is 1, otherwise it is 0;

[0062] The weighted average and maximum method are used to merge the dune densities in four directions to obtain the final dune density estimation result.

[0063] As a preferred solution of the dune density estimation method based on DEM described in the present invention, after the dune density is estimated by the four-directional terrain spectrum decomposition method, a method based on machine learning is used for correction and optimization, including the following steps:

[0064] Collect sample data for training machine learning models, and normalize, clean and partition the collected sample data;

[0065] Extracting topographic features, texture features and directional dune density from DEM data and its derivative data, which are helpful for dune density correction;

[0066] Select the extracted features, calculate the correlation coefficient between each feature and the target output, select the features with high correlation, use a tree-based model to evaluate the importance of each feature, and delete the features with low importance;

[0067] Select a set of optimal features as input to the machine learning model;

[0068] According to the characteristics of the data and the requirements of the correction task, a random forest regression model is selected, and the selected random forest regression model is trained using the training set;

[0069] The prediction formula of random forest is:

[0070]

[0071] in, is the final prediction value, T is the number of decision trees, f t (x) is the predicted value of the tth decision tree;

[0072] During the training process, the validation set is used to evaluate the model and adjust the model's hyperparameters;

[0073] After the model training is completed, the trained random forest model is used to predict the data of the test set and the actual dune area to obtain the corrected dune density value;

[0074] The neighborhood average-based smoothing method is used to optimize and smooth the dune density results after preliminary correction. The expression is:

[0075]

[0076] in, is the density of sand dunes after smoothing, N(x,y) is the neighborhood of point (x,y), and |N| is the number of pixels in the neighborhood.

[0077] As a preferred solution of the dune density estimation method based on DEM described in the present invention, the method of estimating and fusing the dune density at multiple scales by multi-scale decomposition of the dune density comprises the following steps:

[0078] The DEM data is decomposed into multiple scales using discrete wavelet transform to obtain detail coefficients and approximation coefficients at different scales.

[0079] The scale coefficients obtained by wavelet decomposition are subjected to inverse wavelet transform to reconstruct the corresponding terrain features at different scales;

[0080] In the reconstructed terrain at each scale, the local extremum detection method is used to identify the dune vertices. Based on the local extremum detection at each scale, the dune density at each scale is calculated, which is defined as the number of dunes per unit area.

[0081] The weighted average method is used to merge the dune densities at different scales to obtain a comprehensive dune density estimate.

[0082] As a preferred solution of the dune density estimation method based on DEM described in the present invention, the comprehensive analysis of multi-scale dune density results and the visual display thereof include the following steps:

[0083] Comprehensively analyze the fused dune density results to identify the distribution patterns and characteristics of dunes at different scales;

[0084] The dune density results after multi-scale fusion are compared with the density results obtained by actual observation data and other methods to verify the effectiveness of the method;

[0085] In order to intuitively display the overall distribution characteristics of sand dunes, a comprehensive sand dune density distribution map is generated;

[0086] Based on the results of density fluctuation analysis, a new image is generated to show the scale effect of dune density;

[0087] The comprehensive dune density distribution map, density fluctuation distribution map and actual geographic information are superimposed to form a comprehensive visual display.

[0088] In a second aspect, an embodiment of the present invention provides a computer device, comprising a memory and a processor, wherein the memory stores a computer program, wherein: when the computer program is executed by the processor, any step of the DEM-based dune density estimation method as described in the first aspect of the present invention is implemented.

[0089] In a third aspect, an embodiment of the present invention provides a computer-readable storage medium having a computer program stored thereon, wherein: when the computer program is executed by a processor, any step of the DEM-based dune density estimation method as described in the first aspect of the present invention is implemented.

[0090] The beneficial effects of the present invention are as follows: by preprocessing DEM data, the quality improvement and standardization of the original data are achieved, and the accuracy of subsequent dune vertex extraction and density estimation is effectively improved; by using the eight-domain Laplace operator to perform second-order derivative calculations, the vertex features of the dunes are extracted, and the accurate description of the microscopic morphology of the dunes is achieved, ensuring the basic data quality of subsequent density estimation, avoiding misjudgment and loss in the vertex extraction process, and improving the accuracy of dune density estimation; by performing distance-based clustering analysis on dune vertices, the connectivity analysis and preliminary density estimation of the dune area are achieved, and the density of the dune area is accurately estimated; by using the four-directional terrain spectrum decomposition method, the local extreme value detection of high-frequency components in different directions is performed to identify the dune structural characteristics in each direction, and by fusing the results in different directions, a more comprehensive dune density distribution map is obtained, which greatly improves the accuracy of dune density estimation, avoids the limitations of single-scale or single-directional analysis, and ensures the multi-dimensional accuracy of dune density estimation results. BRIEF DESCRIPTION OF THE DRAWINGS

[0091] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings required for use in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other accompanying drawings can be obtained based on these accompanying drawings without paying creative work.

[0092] Figure 1 This is a flow chart of the dune density estimation method based on DEM in Example 1.

[0093] Figure 2 This is a flow chart for calculating the final dune density estimation result in Example 1. DETAILED DESCRIPTION

[0094] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the specific implementation methods of the present invention are described in detail below in conjunction with the accompanying drawings.

[0095] In the following description, many specific details are set forth to facilitate a full understanding of the present invention, but the present invention may also be implemented in other ways different from those described herein, and those skilled in the art may make similar generalizations without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0096] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The term "in one embodiment" that appears in different places in this specification does not necessarily refer to the same embodiment, nor does it refer to a separate or selective embodiment that is mutually exclusive with other embodiments.

[0097] Example 1, reference Figure 1 , which is the first embodiment of the present invention, provides a method for estimating sand dune density based on DEM, comprising the following steps:

[0098] S1. Preprocess DEM data;

[0099] Read DEM data from storage media or databases. Common DEM data formats include GeoTIFF, ASCIIGrid, etc. When reading data, ensure the integrity of the data and record its metadata information, such as resolution, projection coordinate system, coverage, etc.;

[0100] After reading the data, a preliminary check is performed to ensure that the data is not missing or damaged. The check includes: invalid value detection: identifying and recording invalid values ​​(such as NaN or -9999), which usually represent areas with no data or sensor failures; range consistency check: checking whether the data is within a reasonable elevation range to avoid extreme values ​​caused by erroneous data affecting subsequent processing;

[0101] Use statistical methods to detect global outliers, which are usually significant deviations caused by sensor errors or data conversion errors; use neighborhood averaging to detect local outliers, which are points where the elevation value deviates significantly from that of its neighbors;

[0102] When processing multi-source DEM data, the resolution of different data may be different. In order to unify the resolution, the interpolation method can be used to increase the low-resolution data to the scale of high-resolution data. Common interpolation methods include bilinear interpolation and bicubic interpolation.

[0103] When using DEM data from multiple sources, data fusion is required, taking into account the accuracy differences and coverage areas of different data sources. The weighted average method is used for fusion, and the weight can be determined according to the credibility or resolution of the data source.

[0104] Since DEM data may contain noise, filtering is required. Morphological filtering is an effective method that can retain the main features of the terrain while removing noise. After filtering, the slope and curvature of the DEM data are calculated as auxiliary information for subsequent terrain analysis.

[0105] S2, based on the pre-processed DEM data, the dune vertices are extracted using the eight-domain Laplacian operator, including the following steps:

[0106] Apply the eight-domain Laplace operator to the preprocessed DEM data to calculate the second-order derivative value of each pixel. The expression is:

[0107]

[0108] Among them, h(i,j) represents the elevation value at position (i,j) in the DEM data;

[0109] according to The calculation result is That is, the height of the middle point h(i,j) is greater than that of the other surrounding points, indicating that point h(i,j) is the top of the dune. At this time, the point is assigned a logical value of 1. Similarly, when When , it means that point h(i,j) is the bottom of the dune, and the point is assigned a logical value of 0;

[0110] Through the extracted logical values ​​of 0 and 1, the dune vertices are discretized into two dimensions, that is, the three-dimensional DEM data representing the height is mapped into two-dimensional space;

[0111] This calculation will generate a numerical value for each pixel, representing its local elevation change characteristics;

[0112] The candidate vertices are smoothed using a Gaussian filter. The Gaussian filter can smooth out high-frequency noise but retain the main terrain features. The smoothed vertices are filtered to remove those vertices with insignificant height changes.

[0113] Find the eight-way connected area for the extracted dune vertices, that is, find the Euclidean distance between the extracted vertices, and determine the dune vertex with the shortest distance that can form the circumscribed rectangle of the eight-way connected area as a dune;

[0114] After smoothing and filtering, in order to avoid too dense vertex distribution, a minimum spacing d can be set min , if the distance between two vertices is less than d min , then keep the vertex with higher height and delete the other one. The distance calculation formula is:

[0115]

[0116] Among them, d(Pi ,P j ) is point P i and point P j The Euclidean distance between two points in a two-dimensional plane, point P i The coordinates of (x i ,y i ), click P j The coordinates of (x j ,y j ), x i and i P i The coordinates of x j and j P j The coordinates of

[0117] In order to avoid too sparse vertex distribution, a maximum spacing d can be set max If the distance between vertices in a certain area is too large, additional vertices are inserted in the area to ensure the uniformity of vertex distribution. The position of the inserted vertex can be determined by the interpolation method. The interpolation formula is,

[0118]

[0119] Among them, P new is the position where the vertex is inserted;

[0120] Mark all dune vertices to generate a binary image, where the vertex position is marked as 1 and other areas are marked as 0;

[0121] The extracted dune vertex data can be saved in a standard geographic information format for integrated analysis with other geographic data. At the same time, a vertex distribution map can be generated to visualize the spatial distribution of dune vertices.

[0122] S3, based on the dune vertices, the connectivity analysis of the dune area and the preliminary estimation of the dune density include the following steps:

[0123] After extracting the dune vertices, a distance-based clustering algorithm is used to cluster the dune vertices to identify the vertex groups belonging to the same dune;

[0124] By setting a radius ∈ and the minimum number of neighbors MinPts, the points containing at least MinPts vertices within the given radius are regarded as core points. The direct or indirect connections between core points form a cluster, which is expressed as,

[0125] Cluster(P)={P′∣dist(P,P′)≤∈and|N ∈ (P)|≥MinPts}

[0126] Among them, P and P' are two vertices, dist(P,P') represents the distance between the two vertices, and N ∈ (P) represents the set of vertices in the neighborhood with P as the center and radius ∈;

[0127] The parameters ∈ and MinPts need to be set according to the actual scale and density of the dunes. Usually, ∈ can be selected according to the average scale of the dunes, and MinPts can be set according to the expected minimum scale of the dunes. After clustering, each cluster represents a dune area.

[0128] For each cluster, the boundary of the dune is identified by the convex hull algorithm. The convex hull is the smallest convex polygon containing all vertices. The convex hull algorithm can effectively identify the outer boundary of the dune area, but due to the complex morphology of the dunes, this may cause some areas to not be correctly identified. Therefore, it is necessary to further correct the boundary in combination with terrain features (such as slope and curvature);

[0129] After identifying each dune area, the connected component labeling algorithm is used to perform topological connection analysis between different dunes. The connected component labeling algorithm assigns connected pixels to the same label. For each dune area, check whether its boundary is connected or overlapped with the boundaries of other dune areas. If so, they are considered as part of the same connected component.

[0130] Based on the identified dune areas, the area of ​​each dune area is calculated by integration method, and the expression is:

[0131] Area(A)=∫∫Adxdy

[0132] Where Area(A) is the area of ​​region A, A represents the two-dimensional projection of a dune region, dx and dy represent the small length increments in the x and y directions, respectively, and the integration range is the boundary of the region;

[0133] The dune area can be calculated based on the resolution of the DEM data by counting the number of grid cells in the area and multiplying the area of ​​a single cell;

[0134] The dune density is calculated using area density to reflect the density of dunes in the study area;

[0135] Combining the connectivity analysis results and density calculation results of the dune area, a connectivity index is introduced to quantify the connectivity of the dune group. The higher the connectivity index, the higher the connectivity of the dunes. By analyzing the connectivity, large-scale dune groups can be identified, while the density analysis helps to understand the concentration of dunes.

[0136] Due to the extremely irregular shapes of dunes, some secondary dunes have repeated vertices, which leads to an overestimation of dune density. By determining the connected areas of dunes, the overestimation of dune density is greatly reduced and the accuracy of the estimation is improved.

[0137] S4. Further calculation of the preliminary estimation results of sand dune density is performed by combining the four-directional terrain spectrum decomposition method, including the following steps:

[0138] The four-directional terrain spectrum decomposition method requires the selection of four main directions for spectrum decomposition. The directions usually selected are east-west (EW), north-south (NS), northeast-southwest (NE-SW) and northwest-southeast (NW-SE). These four directions can cover the main terrain change trends and help to fully capture the structural characteristics of sand dunes.

[0139] Perform a two-dimensional fast Fourier transform on the DEM data to decompose the terrain data into components of different frequencies for subsequent directional analysis and filtering, and convert the terrain data in the spatial domain to the frequency domain. The calculation formula is:

[0140]

[0141] Among them, h(x,y) is the elevation value in the spatial domain, H(u,v) is the Fourier coefficient in the frequency domain, M and N are the number of rows and columns of DEM data respectively;

[0142] Through two-dimensional fast Fourier transform, the terrain data representation in the frequency domain is obtained, and these data will be used for the next step of directional filtering processing;

[0143] Using Gaussian bandpass filter, four directional filters are constructed for spectrum decomposition in EW, NS, NE-SW and NW-SE directions respectively. The expression is:

[0144]

[0145] Among them, G(u,v) is the value of the two-dimensional Gaussian function, u and v are the coordinates in the two-dimensional space, u 0 and v 0 is the center coordinate of the Gaussian function, (u 0 ,v 0 ) is the frequency center, θ is the direction angle of the filter, σ u and σ v is the standard deviation of the Gaussian function in the u direction and the v direction, which controls the bandwidth of the filter in the main direction and the vertical direction respectively;

[0146] Four directional filters are applied to the terrain data in the frequency domain to obtain the spectrum components in four directions. The result in each direction represents the main high-frequency components of the terrain in that direction. These components mainly correspond to the morphological characteristics of the sand dunes. Suppose the result after directional filtering is:

[0147] H dir (u,v)=H(u,v)·G dir (u,v)

[0148] Among them, H dir (u,v) is the Fourier coefficient in the two-dimensional frequency domain, G dir (u,v) represents a directional filter in a specific direction;

[0149] Perform inverse two-dimensional Fourier transform on the spectrum components in the four directions, transform the filtering results in the frequency domain back to the spatial domain, and thus obtain the main high-frequency components of the terrain in each direction. The formula is:

[0150]

[0151] Among them, h dir (x, y) is the function value in the two-dimensional space domain;

[0152] Through inverse transformation, the spatial distribution characteristics of sand dunes in all directions were obtained;

[0153] For each direction, the high frequency component h dir (x, y) performs local extreme value detection to extract possible dune vertices. These local extreme value points correspond to the peaks of the dunes and are usually manifested as local maxima of high-frequency components. The formula is:

[0154] h max (x,y)={h dir (x,y)|h dir (x,y)>h N(x,y)}

[0155] Among them, h N(x,y) Represents the set of elevation values ​​within the neighborhood;

[0156] Calculate the local dune density in each direction. The local dune density can be defined as the number of dunes per unit area, or the distribution density of dune vertices. The calculation formula is:

[0157]

[0158] Among them, A is the neighborhood area, δ(h max (i,j)) is the indicator function, if h max If (i,j) has a local extreme value, the value is 1, otherwise it is 0;

[0159] The weighted average and maximum methods are used to merge the dune densities in four directions to obtain the final dune density estimation result. The weighted average method can balance the influence of various directions, while the maximum method can highlight the direction with the highest local dune density.

[0160] S5. After estimating the dune density by the four-directional terrain spectrum decomposition method, a machine learning-based method is used for correction and optimization, including the following steps:

[0161] Collect sample data for training machine learning models. The sample data should include the following parts: input features, including local terrain features of DEM data, directional dune density estimation results, terrain slope, curvature and other information; target output, actual observed or accurately annotated dune density data, which can be obtained through field measurements or high-resolution remote sensing image analysis;

[0162] The collection of sample data should cover different types of dune areas to ensure the generalization ability of the model;

[0163] Normalize input features to the same magnitude to prevent large numerical differences from causing unstable model training. Remove outliers from the data, fill in missing values, ensure data integrity and consistency, and divide the data into training set, validation set, and test set. The training set is used for model training, the validation set is used for parameter tuning, and the test set is used for final model evaluation.

[0164] Extract terrain features that are helpful for dune density correction from DEM data and its derivative data, such as slope, curvature, elevation, eight-domain Laplacian operator results, etc.; texture features, such as texture features (contrast, energy, entropy, etc.) based on gray-level co-occurrence matrix (GLCM), which can reflect the surface details of the terrain; directional dune density, which comes from the dune density estimation value of the four-directional terrain spectrum decomposition method;

[0165] Select the extracted features, calculate the correlation coefficient between each feature and the target output, select the features with high correlation, use a tree-based model to evaluate the importance of each feature, and delete the features with low importance;

[0166] Select a set of optimal features as input to the machine learning model;

[0167] According to the characteristics of the data and the requirements of the correction task, a random forest regression model is selected and trained using the training set. The training process involves the construction of multiple decision trees and reduces the overfitting risk of a single decision tree through ensemble learning.

[0168] The prediction formula of random forest is:

[0169]

[0170] in, is the final prediction value, T is the number of decision trees, f t (x) is the predicted value of the tth decision tree;

[0171] During the training process, the model is evaluated using the validation set and the model's hyperparameters (such as the number of trees, maximum depth, etc.) are adjusted to optimize model performance;

[0172] After the model training is completed, the trained random forest model is used to predict the data of the test set and the actual dune area to obtain the corrected dune density value;

[0173] Since the prediction results of the random forest model may have local discontinuities or fluctuations, a smoothing method based on neighborhood average can be used to optimize and smooth the dune density results after preliminary correction. The expression is:

[0174]

[0175] in, is the density of sand dunes after smoothing, N(x,y) is the neighborhood of point (x,y), and |N| is the number of pixels in the neighborhood;

[0176] Through smoothing, local anomalies in dune density can be eliminated and the continuity and rationality of density distribution can be improved.

[0177] S6. Estimating and integrating the dune density at multiple scales by decomposing the dune density at multiple scales includes the following steps:

[0178] The DEM data is decomposed into multiple scales by using discrete wavelet transform to obtain detail coefficients and approximation coefficients at different scales, which reflect the dune structures at different scales.

[0179] The scale coefficients obtained by wavelet decomposition are subjected to inverse wavelet transform to reconstruct the corresponding terrain features at different scales. Through inverse transform, the terrain features at each scale can be restored, thus constructing a multi-scale terrain image.

[0180] In the reconstructed terrain at each scale, the local extremum detection method is used to identify the dune vertices, which represent the main positions of the dunes at that scale. Based on the local extremum detection at each scale, the dune density at each scale is calculated, which is defined as the number of dunes per unit area.

[0181] The weighted average method is used to merge the dune densities at different scales to obtain a comprehensive dune density estimate. Through multi-scale density fusion, the dune characteristics of each scale can be comprehensively considered to obtain a more accurate and comprehensive dune density distribution.

[0182] S7. Comprehensively analyze the multi-scale dune density results and make a visual display, including the following steps:

[0183] Comprehensively analyze the fused dune density results to identify the distribution patterns and characteristics of dunes at different scales. For example, the spatial distribution of dunes: by analyzing the density distribution in different areas, identify the concentrated and sparse areas of dunes; scale effect of dunes: by comparing the density estimation results at different scales, analyze the scale effect of dunes and identify the main scale range of dunes;

[0184] The mean square error and determination coefficient were used as evaluation indicators to compare the dune density results after multi-scale fusion with the density results obtained by actual observation data and other methods to verify the effectiveness of the method.

[0185] In order to intuitively display the overall distribution characteristics of sand dunes and generate a comprehensive sand dune density distribution map, the image can use color gradients to represent the density, and the darker the color, the higher the density;

[0186] Based on the results of the density fluctuation analysis, a new image is generated to show the scale effect of dune density. Use a contour map or color map to show density fluctuations, highlighting those areas where density changes dramatically, and choose appropriate colors or contour intervals so that the observer can quickly identify areas of density fluctuations. This image can be used to show the standard deviation of density through contour maps or color maps;

[0187] The comprehensive dune density distribution map and density fluctuation distribution map are superimposed with the actual geographic information to form a comprehensive visualization. Select appropriate transparency, superimpose the geographic information layer with the dune density distribution layer, and adjust the superimposed map so that the geographic information and dune density information can coexist harmoniously and complement each other to generate the final comprehensive visualization image, which can be displayed to users or used for further analysis. This image can intuitively show the relationship between dune density and geographical environment.

[0188] This embodiment also provides a computer device, which is suitable for the case of a DEM-based dune density estimation method, including: a memory and a processor; the memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions to implement the DEM-based dune density estimation method proposed in the above embodiment.

[0189] The computer device may be a terminal, and the computer device includes a processor, a memory, a communication interface, a display screen and an input device connected via a system bus. The processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system and a computer program. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The communication interface of the computer device is used to communicate with an external terminal in a wired or wireless manner, and the wireless manner can be achieved through WIFI, an operator network, NFC (near field communication) or other technologies. The display screen of the computer device may be a liquid crystal display screen or an electronic ink display screen, and the input device of the computer device may be a touch layer covering the display screen, or a key, trackball or touchpad provided on the housing of the computer device, or an external keyboard, touchpad or mouse, etc.

[0190] This embodiment also provides a storage medium on which a computer program is stored. When the program is executed by a processor, the method for estimating the dune density based on the DEM as proposed in the above embodiment is implemented; the storage medium can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as static random access memory (Static Random Access Memory, referred to as SRAM), electrically erasable programmable read-only memory (Electrically Erasable Programmable Read-Only Memory, referred to as EEPROM), erasable programmable read-only memory (Erasable Programmable Read Only Memory, referred to as EPROM), programmable read-only memory (Programmable Red-Only Memory, referred to as PROM), read-only memory (Read-Only Memory, referred to as ROM), magnetic storage, flash memory, magnetic disk or optical disk.

[0191] In summary, the present invention achieves quality improvement and standardization of raw data by: preprocessing DEM data, effectively improving the accuracy of subsequent dune vertex extraction and density estimation; extracting the vertex features of dunes by using eight-domain Laplace operator for second-order derivative calculation, achieving accurate description of dune micromorphology, ensuring the basic data quality of subsequent density estimation, avoiding misjudgment and loss in the vertex extraction process, and improving the accuracy of dune density estimation; realizing connectivity analysis and preliminary density estimation of dune areas by performing distance-based clustering analysis on dune vertices, accurately estimating the density of dune areas; detecting local extreme values ​​of high-frequency components in different directions by four-directional terrain spectral decomposition method, identifying dune structural features in various directions, and obtaining a more comprehensive dune density distribution map by fusing the results in different directions, which greatly improves the accuracy of dune density estimation, avoids the limitations of single-scale or single-directional analysis, and ensures the multi-dimensional accuracy of dune density estimation results.

[0192] Example 2, referring to Table 1, is the second example of the present invention. To further verify the technical solution of the present invention, experimental simulation data of a dune density estimation method based on DEM is provided.

[0193] This embodiment aims to verify the effectiveness of a dune density estimation method based on DEM, and compares it with the prior art to show the superiority of this method in dune density estimation. A typical desert area was selected as the test area, with an area of ​​100 square kilometers and a resolution of 30 meters for DEM data. First, the acquired DEM data was preprocessed, and the preprocessing process included: data reading, outlier detection and correction, resolution standardization (all data were standardized to a resolution of 30 meters), data fusion and noise elimination. Through these steps, the accuracy and consistency of the DEM data were ensured, laying the foundation for subsequent analysis.

[0194] Next, based on the preprocessed DEM data, the eight-domain Laplacian operator is used to extract the dune vertices. The specific process includes calculating the second-order derivative of each pixel to locate the area with significant local elevation changes, and combining it with a Gaussian filter for smoothing to further remove insignificant vertices. In order to avoid too dense or too sparse vertex distribution, the minimum spacing is set to 50 meters and the maximum spacing is set to 200 meters to ensure a reasonable vertex distribution.

[0195] After vertex extraction, a distance-based clustering algorithm was used to cluster the vertices. The radius parameter was set to 100 meters and the minimum number of neighbors was 3 vertices to define the core point. After clustering, the convex hull algorithm was used to determine the boundaries of the dune area. Based on the area of ​​each dune area, the preliminary dune density was calculated, and the interconnection between the dunes was quantified through connectivity analysis.

[0196] Subsequently, the preliminary estimate of dune density was further calculated by combining the four-directional terrain spectrum decomposition method. This step uses the two-dimensional fast Fourier transform to convert the DEM data into the frequency domain, applies four directional filters to capture the terrain features in different directions, and returns the results to the spatial domain through the inverse Fourier transform. By fusing the dune densities in four directions, a more accurate dune density estimate was obtained.

[0197] Finally, the accuracy of density estimation was further improved by using machine learning methods for correction and optimization, and using random forest regression models for training and testing based on data from actual dune areas. Finally, the dune density was decomposed and fused at multiple scales to generate dune density estimation results at multiple scales, and visualization technology was used to intuitively display the dune distribution and density changes.

[0198] Table 1 Comparison of sand dune density estimates

[0199]

[0200] It can be seen from the experimental results that the method of the present invention is superior to the prior art in multiple key indicators, and the specific analysis is as follows:

[0201] The two methods are consistent in data accuracy, both maintaining a resolution of 30 meters. This is because the resolution standardization operation in the preprocessing step makes the spatial resolution of the DEM data uniform, ensuring the comparability of subsequent analysis.

[0202] The vertex extraction accuracy of the prior art is 85%, while the method of the present invention significantly improves the accuracy of dune vertex extraction to 95% by combining the eight-domain Laplacian operator and the Gaussian filter. The reason for this difference is that the method of the present invention effectively avoids the phenomenon of overly dense or sparse vertex distribution by reasonably controlling the vertex density (setting the minimum and maximum spacing), thereby improving the overall extraction accuracy. More accurate vertex extraction provides a reliable data basis for subsequent dune density estimation.

[0203] The connectivity index of the prior art is 0.65, while the method of the present invention achieves a connectivity index of 0.85 through a distance-based clustering algorithm and a convex hull algorithm. The improvement in the connectivity index shows that the method of the present invention can more accurately capture the connection relationship between sand dunes, thereby better reflecting the overall structure of the dune group. This improvement is of great significance in practical applications, because dune connectivity is a key indicator for evaluating dune migration and desertification progress.

[0204] The preliminary estimate of sand dune density in the prior art is 12 sand dunes / square kilometer, while the method of the present invention is 15 sand dunes / square kilometer. This difference is mainly due to the improvement of the method of the present invention in vertex extraction and regional connectivity analysis. Through more accurate vertex positioning and cluster analysis, the sand dune area can be more completely identified, avoiding the situation where the sand dunes are underestimated.

[0205] After correction and optimization by machine learning, the final dune density of the prior art is estimated to be 14 dunes / square kilometer, while the density after correction by the method of the present invention is 18 dunes / square kilometer. The advantage of the method of the present invention is that it combines the machine learning algorithm and further improves the accuracy of density estimation by multi-scale decomposition and fusion of dune density. Compared with the prior art, the correction result of the method of the present invention is closer to the actual observation data, reflecting its superiority in complex terrain.

[0206] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention, which should all be included in the scope of the claims of the present invention.

Claims

1. A method for estimating sand dune density based on DEM, characterized by: include, Preprocess DEM data; Based on the preprocessed DEM data, the dune vertices are extracted using the eight-domain Laplacian operator; Based on the dune vertices, the connectivity analysis of the dune area and the preliminary estimation of the dune density were carried out; specifically, after the dune vertices were extracted, the distance-based clustering algorithm was used to cluster the dune vertices, and each cluster represented a dune area; after the various dune areas were identified, the connected component labeling algorithm was used to perform topological connection analysis between different dunes; based on the identified dune areas, the area of ​​each dune area was calculated by the integral method, and the dune density was calculated by the area density; The preliminary estimation results of sand dune density are further calculated in combination with the four-directional terrain spectrum decomposition method. Specifically, the DEM data is subjected to a two-dimensional fast Fourier transform to decompose the terrain data into components of different frequencies, and the terrain data in the spatial domain is converted to the frequency domain. A Gaussian bandpass filter is used to construct four directional filters. Four directional filters are applied to the terrain data in the frequency domain to obtain the spectrum components in four directions. The result in each direction represents the main high-frequency component of the terrain in that direction, corresponding to the morphological characteristics of the sand dunes. The spectrum components in the four directions are subjected to an inverse two-dimensional Fourier transform to transform the filtering results in the frequency domain back to the spatial domain to obtain the spatial distribution characteristics of the sand dunes in each direction. The local sand dune density in each direction is calculated. The local sand dune density can be defined as the number of sand dunes per unit area, or the distribution density of the dune vertices. The weighted average and maximum method is used to fuse the sand dune densities in the four directions to obtain the final sand dune density estimation result. After estimating the dune density through the four-directional terrain spectrum decomposition method, a machine learning-based method was used for correction and optimization; By decomposing the dune density at multiple scales, the dune density at multiple scales can be estimated and integrated; Comprehensively analyze the multi-scale dune density results and make a visual display.

2. The method for estimating sand dune density based on DEM as claimed in claim 1, characterized in that: The DEM data preprocessing includes reading and preliminary checking of DEM data, detection and correction of outliers, resolution standardization of DEM data, fusion of multi-source DEM data, terrain filtering and noise elimination.

3. The method for estimating sand dune density based on DEM as claimed in claim 2, characterized in that: The method of extracting the dune vertices based on the pre-processed DEM data using the eight-domain Laplace operator includes the following steps: Apply the eight-domain Laplace operator to the preprocessed DEM data to calculate the second-order derivative value of each pixel. The expression is: Among them, h(i,j) represents the elevation value at position (i,j) in the DEM data; This calculation will generate a numerical value for each pixel, representing its local elevation change characteristics; Use Gaussian filter to smooth the candidate vertices, and filter the smoothed vertices to remove those vertices with insignificant height changes; After smoothing and filtering, in order to avoid too dense vertex distribution, a minimum spacing d can be set min , if the distance between two vertices is less than d min , then keep the vertex with higher height and delete the other one. The distance calculation formula is: Among them, d(P i ,P j ) is point P i and point P j The Euclidean distance between two points in a two-dimensional plane, point P i The coordinates of (x i ,y i ), click P j The coordinates of (x j ,y j ), x i and i P i The coordinates of x j and j P j The coordinates of In order to avoid too sparse vertex distribution, a maximum spacing d can be set max If the distance between vertices in a region is too large, then additional vertices are inserted in the region. The positions of the inserted vertices can be determined by the interpolation method. The interpolation formula is: Among them, P new is the position where the vertex is inserted; Mark all dune vertices to generate a binary image, where the vertex position is marked as 1 and other areas are marked as 0; The extracted dune vertex data are saved in a standard geographic information format to generate a vertex distribution map.

4. The method for estimating sand dune density based on DEM as claimed in claim 3, characterized in that: The connectivity analysis of the dune area and the preliminary estimation of the dune density based on the dune vertices include the following steps: After extracting the dune vertices, a distance-based clustering algorithm is used to cluster the dune vertices; By setting a radius ∈ and a minimum number of neighbors MinPts, points containing at least MinPts vertices within a given radius are considered core points. The direct or indirect connections between core points form a cluster, expressed as, Cluster(P)={P′∣dist(P,P′)≤∈and|N ∈ (P)|≥MinPts} Among them, P and P' are two vertices, dist(P,P') represents the distance between the two vertices, and N ∈ (P) represents the set of vertices in the neighborhood with P as the center and radius ∈; The parameters ∈ and MinPts need to be set according to the actual scale and density of the dunes. After clustering, each cluster represents a dune area; For each cluster, the convex hull algorithm is used to identify the boundaries of the dunes; After identifying each dune area, the connected component labeling algorithm was used to perform topological connectivity analysis between different dunes; Based on the identified dune areas, the area of ​​each dune area is calculated by integration method, and the expression is: Area(A)=∫∫ A dxdy Where Area(A) is the area of ​​region A, A represents the two-dimensional projection of a dune region, dx and dy represent the small length increments in the x and y directions, respectively, and the integration range is the boundary of the region; The dune density is calculated using area density to reflect the density of dunes in the study area; Combining the connectivity analysis results and density calculation results of the dune area, a connectivity index was introduced to quantify the degree of connectivity of the dune group. The higher the connectivity index, the higher the degree of connectivity of the dunes.

5. The method for estimating sand dune density based on DEM as claimed in claim 4, characterized in that: The further calculation of the preliminary estimation result of the sand dune density by combining the four-directional terrain spectrum decomposition method includes the following steps: Perform a two-dimensional fast Fourier transform on the DEM data to convert the terrain data in the spatial domain into the frequency domain. The calculation formula is: Among them, h(x,y) is the elevation value in the spatial domain, H(u,v) is the Fourier coefficient in the frequency domain, M and N are the number of rows and columns of DEM data respectively; Using Gaussian bandpass filter, four directional filters are constructed, the expression is, Among them, G(u,v) is the value of the two-dimensional Gaussian function, u and v are the coordinates in the two-dimensional space, u0 and v0 are the center position coordinates of the Gaussian function, (u0,v0) is the frequency center, θ is the direction angle of the filter, σ u and σ v is the standard deviation of the Gaussian function in the u direction and the v direction, which controls the bandwidth of the filter in the main direction and the vertical direction respectively; Four directional filters are applied to the terrain data in the frequency domain to obtain the spectrum components in four directions. The result in each direction represents the main high-frequency components of the terrain in that direction. These components mainly correspond to the morphological characteristics of the sand dunes. Suppose the result after directional filtering is: H dir (u,v)=H(u,v)·G dir (u,v) Among them, H dir (u,v) is the Fourier coefficient in the two-dimensional frequency domain, G dir (u,v) represents a directional filter in a specific direction; Perform inverse two-dimensional Fourier transform on the spectrum components in four directions respectively, and transform the filtering results in the frequency domain back to the spatial domain. The formula is: Among them, h dir (x, y) is the function value in the two-dimensional space domain; Through inverse transformation, the spatial distribution characteristics of sand dunes in all directions were obtained; For each direction, the high frequency component h dir (x,y) performs local extreme value detection, the formula is, h max (x,y)={h dir (x,y)∣h dir (x,y)>h N(x,y) } Among them, h N(x,y) Represents the set of elevation values ​​within the neighborhood; Calculate the local dune density in each direction using the formula, Among them, A is the neighborhood area, δ(h max (i,j)) is the indicator function, if h max If (i,j) has a local extreme value, the value is 1, otherwise it is 0; The weighted average and maximum method are used to merge the dune densities in four directions to obtain the final dune density estimation result.

6. The method for estimating sand dune density based on DEM as claimed in claim 5, characterized in that: After the dune density is estimated by the four-directional terrain spectrum decomposition method, a machine learning-based method is used for correction and optimization, including the following steps: Collect sample data for training machine learning models, and normalize, clean and partition the collected sample data; Extracting topographic features, texture features and directional dune density from DEM data and its derivative data, which are helpful for dune density correction; Select the extracted features, calculate the correlation coefficient between each feature and the target output, select the features with high correlation, use a tree-based model to evaluate the importance of each feature, and delete the features with low importance; Select a set of optimal features as input to the machine learning model; According to the characteristics of the data and the requirements of the correction task, a random forest regression model is selected, and the selected random forest regression model is trained using the training set; The prediction formula of random forest is: in, is the final prediction value, T is the number of decision trees, f t (x) is the predicted value of the tth decision tree; During the training process, the validation set is used to evaluate the model and adjust the model's hyperparameters; After the model training is completed, the trained random forest model is used to predict the data of the test set and the actual dune area to obtain the corrected dune density value; The neighborhood average-based smoothing method is used to optimize and smooth the dune density results after preliminary correction. The expression is: in, is the density of sand dunes after smoothing, N(x,y) is the neighborhood of point (x,y), and |N| is the number of pixels in the neighborhood.

7. The method for estimating sand dune density based on DEM according to claim 6, characterized in that: The method of decomposing the dune density at multiple scales, estimating and integrating the dune density at multiple scales includes the following steps: The DEM data is decomposed into multiple scales using discrete wavelet transform to obtain detail coefficients and approximation coefficients at different scales. The scale coefficients obtained by wavelet decomposition are subjected to inverse wavelet transform to reconstruct the corresponding terrain features at different scales; In the reconstructed terrain at each scale, the local extremum detection method is used to identify the dune vertices. Based on the local extremum detection at each scale, the dune density at each scale is calculated, which is defined as the number of dunes per unit area. The weighted average method is used to merge the dune densities at different scales to obtain a comprehensive dune density estimate.

8. The method for estimating sand dune density based on DEM according to claim 6, characterized in that: The comprehensive analysis of multi-scale dune density results and the visual display thereof include the following steps: Comprehensively analyze the fused dune density results to identify the distribution patterns and characteristics of dunes at different scales; The dune density results after multi-scale fusion are compared with the density results obtained by actual observation data and other methods to verify the effectiveness of the method; In order to intuitively display the overall distribution characteristics of sand dunes, a comprehensive sand dune density distribution map is generated; Based on the results of density fluctuation analysis, a new image is generated to show the scale effect of dune density; The comprehensive dune density distribution map, density fluctuation distribution map and actual geographic information are superimposed to form a comprehensive visual display.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the steps of the method for estimating sand dune density based on DEM according to any one of claims 1 to 8 are implemented.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the method for estimating sand dune density based on DEM according to any one of claims 1 to 8 are implemented.

Images