Method for automatically extracting elevation points of DEM (Digital Elevation Model)

Abstract

The invention belongs to the technical field of elevation data analysis, and provides an automatic elevation point extraction method for a DEM (Digital Elevation Model), which comprises the following steps of: clustering elevation data of sampling points by setting different clustering radiuses, judging terrain integrity of each clustering range, and dividing subspaces; analyzing reference values of elevation interpolation of different subspaces on the projection position point to obtain a reference elevation value of the projection position point; then analyzing the distance interval condition of each elevation data in the same region in different subspaces, and adjusting the coordinate position of the projection position point to obtain a corrected DEM space model; and finally, analyzing the difference between the water flow trend and the actual terrain in the corrected DEM space model, and adjusting the reference elevation values of partial projection position points to obtain a secondary corrected DEM space model. By means of the method, the accuracy of obtaining interpolation elevation points at different positions is improved, and then the accuracy of constructing the DEM space model is improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of elevation data analysis, and particularly relates to a DEM spatial model elevation point automatic extraction method. BACKGROUND

[0002] In a digital elevation model (DEM), it is crucial to ensure the accuracy of elevation points. High-precision elevation data is the basis for terrain analysis (such as slope, slope direction, and watershed division), and errors will be significantly transmitted to all derived analysis results. Accurate elevation points can more realistically reflect the details of the topography, ensuring the scientificity and reliability of applications such as hydrological modeling, flood simulation, engineering design, earthwork calculation, and disaster assessment, and avoiding decision-making errors or resource waste caused by deviations in basic data. Therefore, pursuing higher precision of elevation points is a core requirement for improving the application value of DEM.

[0003] The traditional method obtains an interpolated elevation point by using a global spatial data as a reference through a Kriging algorithm, and the influence of the global data on a single terrain condition is not clear enough, and the obtained interpolated elevation point does not conform to the actual situation of the terrain. SUMMARY

[0004] To solve the above technical problems, the present application provides a DEM spatial model elevation point automatic extraction method.

[0005] According to the DEM spatial model elevation point automatic extraction method provided by the present application, the method comprises:

[0006] obtaining elevation data of a sampling point;

[0007] setting different clustering radii, clustering the elevation data to obtain a clustering range corresponding to each clustering radius;

[0008] analyzing the occurrence of the boundary of each clustering range and the difference of the elevation data on both sides of the boundary of each clustering range under different clustering radii, judging the terrain integrity of each clustering range, and obtaining a plurality of subspaces;

[0009] performing an interpolation operation on the subspaces to obtain a projection position point and an initial elevation value thereof;

[0010] combining the projection area of the subspace, analyzing the reference value of different subspace for calculating the elevation interpolation of the projection position point according to the terrain integrity, and combining the initial elevation value to obtain a reference elevation value of the projection position point;

[0011] analyzing the distance interval of each elevation data in the same region in different subspaces, adjusting the coordinate position of the projection position point in combination with the reference value, and obtaining a corrected DEM spatial model;

[0012] The difference between the water flow trend in the modified DEM spatial model and the actual terrain is analyzed, and the reference elevation value of some projection position points is adjusted to obtain a second modified DEM spatial model.

[0013] In some embodiments of the present application, the terrain integrity of each cluster range is determined by analyzing the occurrence of each cluster range boundary under different cluster radii and the difference of elevation data on both sides of each cluster range boundary, including:

[0014] The total number of times each cluster range boundary appears under all cluster radii is counted to quantify the occurrence of each cluster range boundary under different cluster radii.

[0015] Under the cluster radius containing the same cluster range boundary, the absolute value of the difference between the average elevation data in the cluster range surrounded by the cluster range boundary and the average elevation data in the adjacent cluster range is calculated to obtain the difference degree of the elevation data on both sides of each cluster range boundary.

[0016] The terrain integrity of the cluster range surrounded by each cluster range boundary is obtained in combination with the occurrence and the difference degree.

[0017] In some embodiments of the present application, a plurality of subspaces are obtained, including:

[0018] An integrity threshold is set.

[0019] It is determined whether the terrain integrity is greater than the integrity threshold.

[0020] If yes, the cluster range surrounded by the cluster range boundary is taken as a subspace for elevation interpolation calculation by the Kriging algorithm to obtain a plurality of subspaces.

[0021] In some embodiments of the present application, an interpolation operation is performed on the subspace to obtain a projection position point and its initial elevation value, including:

[0022] An interpolation elevation point is obtained by performing an interpolation operation on the subspace by the Kriging algorithm.

[0023] For the subspace obtained when the cluster radius value is the smallest, an interpolation elevation point corresponding to a projection position point is obtained by the Kriging interpolation algorithm.

[0024] For all subspaces containing a single projection position point, the elevation value of the nearest interpolation elevation point obtained by each subspace is obtained to obtain the initial elevation value corresponding to the projection position point in each subspace.

[0025] In some embodiments of the present application, according to the topographic integrity, in combination with the projection area of the subspace, the reference value of the interpolation height of the projection position point is analyzed in different subspaces, and in combination with the initial height value, the reference height value of the projection position point is obtained, including:

[0026] According to the topographic integrity, in combination with the projection area of the subspace, the reference value of the interpolation height of the projection position point is analyzed in different subspaces;

[0027] According to the reference value and the initial height value, the corresponding weighted height value of the projection position point in each subspace is obtained, and the reference height value of the projection position point is obtained by summing all subspaces containing the projection position point.

[0028] In some embodiments of the present application, the distance interval of each height data in the same region in different subspaces is analyzed, the coordinate position of the projection position point is adjusted in combination with the reference value, and the corrected DEM space model is obtained, including:

[0029] The second Euclidean distance between the projection position point and the corresponding projection position point in different subspaces is calculated, and the movement vector of the projection position point to the corresponding projection position point in different subspaces is obtained in combination with the reference value;

[0030] The movement vectors corresponding to all the subspaces are summed to obtain the movement reference total vector of the projection position point;

[0031] According to the movement reference total vector, the coordinate position of the projection position point is adjusted to obtain the corrected DEM space model.

[0032] In some embodiments of the present application, the difference between the water flow trend in the corrected DEM space model and the actual topography is analyzed, including:

[0033] The area difference degree between the water area in the corrected DEM space model and the actual water area is analyzed, and the flow velocity difference degree between the water flow velocity in the corrected DEM space model and the actual water flow velocity is analyzed.

[0034] In some embodiments of the present application, the area difference degree between the water area in the corrected DEM space model and the actual water area is analyzed, including:

[0035] The contour line is generated by using the corrected DEM space model, the model water area is converted into a polygon model water area according to the height step, the area of the model water area is calculated, and the model water area in the corrected DEM space model is obtained;

[0036] The actual water area corresponding to the model water area is obtained.

[0037] Calculate the difference between the model water area and the actual water area, and obtain the area difference degree of the water area in the corrected DEM spatial model.

[0038] In some embodiments of the present application, the flow velocity difference degree of the water flow in the corrected DEM spatial model and the actual water flow is analyzed, comprising:

[0039] Based on the corrected DEM spatial model, the flow direction matrix is calculated using the flow direction algorithm, and the flow accumulation matrix is generated in combination with the slope. The model water flow velocity of the model water flow in the corrected DEM spatial model is obtained by the accumulation amount of the flow accumulation matrix and the slope.

[0040] Obtain the actual water flow velocity corresponding to the model water flow;

[0041] Obtain the first Euclidean distance between the interpolated elevation point on the terrain of the model water flow and the center point of the model water area, wherein the model water flow is connected with the model water area;

[0042] Calculate the difference between the model water flow velocity and the actual water flow velocity, and obtain the flow velocity difference degree of the water flow in the corrected DEM spatial model and the actual water flow in combination with the first Euclidean distance.

[0043] In some embodiments of the present application, the arrangement method of the sampling points is:

[0044] The sampling points are arranged by the feature point priority method, and the sampling points are arranged every 5-10 meters on the ridge and valley line, and the sampling points are arranged every 1-3 meters at the slope mutation, and the sampling points are arranged every 20-50 meters in the flat area.

[0045] From the above embodiments, the DEM spatial model elevation point automatic extraction method provided by the embodiments of the present application has the following beneficial effects:

[0046] The application judges the terrain integrity of the cluster range surrounded by each cluster range boundary by setting different cluster radii, clustering the elevation data of all sampling points, and analyzing the occurrence of each cluster range boundary and the difference of the elevation data on both sides of each cluster range boundary under different cluster radii, and divides the subspace when the Kriging algorithm is used for elevation interpolation calculation; the subspace is subjected to interpolation operation to obtain the projection position point and the initial elevation value thereof, and the reference elevation value of the projection position point is obtained by combining the terrain integrity and the projection area of the subspace, analyzing the reference value of the different subspace for the elevation interpolation of the projection position point, and combining the initial elevation value; the coordinate position of the projection position point is adjusted to obtain the corrected DEM spatial model by analyzing the distance interval of the elevation data in the same region in different subspaces and combining the reference value; finally, the reference elevation value of part of the projection position points is adjusted to obtain the secondary corrected DEM spatial model by analyzing the difference between the water flow trend and the actual terrain in the corrected DEM spatial model. Through the division of the subspace and the analysis of the influence of each subspace on the projection position point, the coordinate position and the elevation value of the projection position point corresponding to the interpolation elevation point are adjusted, the accuracy of the interpolation elevation point obtained at different positions is improved, and the accuracy of the DEM spatial model is improved.

[0047] It should be understood that the foregoing general description and the following detailed description are only exemplary and explanatory and are not restrictive of the application. BRIEF DESCRIPTION OF DRAWINGS

[0048] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, and the advantages thereof, a brief introduction will be given to the drawings needed in the embodiments or the prior art description. Obviously, the drawings in the following description are only some embodiments of the present application, and for those skilled in the art, other drawings can be obtained without creative labor based on these drawings.

[0049] Figure 1 A basic flowchart of a DEM spatial model elevation point automatic extraction method provided by the embodiment of the present application is shown in the figure.

[0050] Figure 2 A moving reference total vector diagram of a projection position point provided by the embodiment of the present application is shown in the figure. DETAILED DESCRIPTION

[0051] In order to further clarify the technical means and effects taken by the present application to achieve the predetermined inventive objectives, the following describes in detail the specific implementation, structure, features and effects of a DEM spatial model elevation point automatic extraction method according to the present application, in combination with the accompanying drawings and preferred embodiments. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. In addition, the specific features, structures or characteristics in one or more embodiments can be combined in any suitable form.

[0052] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. The use of the terms "including", "containing" or any other variation thereof is intended to cover a non-exclusive inclusion, such that a circuit structure, item or apparatus that includes a list of elements does not only include those elements but can also include other elements not expressly listed or inherent to such item or apparatus. Without more limitations, the element defined by the phrase "including a" does not exclude the presence of additional identical elements in the item or apparatus including the element. The relational terms "first" and "second" and the like are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between the entities or operations.

[0053] The following will describe in detail a DEM spatial model elevation point automatic extraction method provided by the present embodiment in combination with the accompanying drawings.

[0054] Please refer to Figure 1 , which shows the basic flow of a DEM spatial model elevation point automatic extraction method provided by one embodiment of the present application.

[0055] As shown in Figure 1 , the DEM spatial model elevation point automatic extraction method provided by one embodiment of the present application specifically includes the following steps:

[0056] S100: Obtain the elevation data of the sampling points.

[0057] The sampling points are arranged by the feature point priority method, which requires sampling points every 5-10 meters on ridge lines and valley lines, and sampling points are encrypted to 1-3 meters at slope mutation points, and sampling points are supplemented every 20-50 meters in flat areas.

[0058] The DEM is measured by using a total station, and high-precision instruments (such as Leica TS60, angle measurement 0.5'') are selected. Control points need to be arranged 3-5, uniformly covering the measurement area and being in visual communication with the measurement station. When measuring, the instrument is stably erected (centering error <1mm), and a single precision measurement mode (2-3 times to take the average value) is adopted, and complex terrains such as cliffs are measured by multi-station intersection. Finally, the elevation data of multiple sampling points are obtained.

[0059] S200: different clustering radii are set, the elevation data is clustered, and the clustering ranges corresponding to different clustering radii are obtained.

[0060] The Kriging algorithm for calculating the elevation points in the DEM spatial model may cause the smoothing of ridges and valleys due to the modeling deviation of the variogram function. Therefore, the DEM spatial model can be divided into multiple subspaces with consistent terrain for calculating the elevation data of each position point.

[0061] Therefore, in the embodiments of the present application, first, different clustering radii are set, the elevation data is clustered, and the clustering ranges corresponding to different clustering radii are obtained. Specifically, the DBSCAN density clustering method is used to cluster the elevation data of different sampling points, and the clustering range is changed by modifying the radius parameter ε value. Among them, the value of the radius parameter ε is set to 30m, 25m, 20m, 15m and 10m, respectively. The elevation data of all measured sampling points is clustered with different radius parameter ε values, and the clustering range corresponding to each radius parameter ε value is obtained. It should be noted that the clustering ranges corresponding to different radius parameter ε values may or may not have the same clustering range.

[0062] S300: the occurrence of each clustering range boundary under different clustering radii and the difference of the elevation data on both sides of each clustering range boundary are analyzed, the terrain integrity of each clustering range is judged, and multiple subspaces are obtained.

[0063] After obtaining the range by clustering, the occurrence of each clustering range boundary under different clustering radii and the difference of the elevation data on both sides of each clustering range boundary are analyzed, wherein the clustering range boundary refers to the separation line between different clusters (different clustering ranges) in the clustering result, the terrain integrity of the clustering range surrounded by each clustering range boundary is judged, and the subspace for Kriging algorithm for elevation interpolation calculation is obtained. Further including:

[0064] First, the total number of times each clustering range boundary appears under all clustering radii is counted, and the occurrence of each clustering range boundary under different clustering radii is quantified. Specifically, for the boundaries of multiple clustering ranges obtained under different radius parameters ε values, the total number of times n jThe occurrence of each cluster range boundary under different cluster radii is quantified.

[0065] Then, under the cluster radius containing the same cluster boundary, the absolute value of the difference between the mean elevation data within the cluster bounded by the cluster boundary and the mean elevation data within its adjacent clusters is calculated to obtain the degree of difference in elevation data on both sides of each cluster boundary. Specifically, for clustering operations with different radius parameter ε values ​​for cluster boundary j, the mean elevation data of all sampling points contained in the cluster bounded by cluster boundary j and its multiple adjacent clusters are calculated separately. Then, the sum of the absolute values ​​of the differences between the mean elevation data of the cluster bounded by cluster boundary j and the mean elevation data of all other adjacent clusters is calculated, denoted as H. j,ε Furthermore, iterate through all radius parameter ε values ​​to obtain the H value corresponding to each radius parameter ε value. j,ε Then calculate H corresponding to all radius parameter ε values. j,ε and (z represents the number of values ​​for the radius parameter ε), thus obtaining the degree of difference in elevation data on both sides of the boundary j of each cluster range.

[0066] Furthermore, by combining the occurrence and degree of difference, the topographic integrity of the cluster range enclosed by the boundary of each cluster range is obtained. Specifically, the more times (nx) different cluster range boundaries j appear under all radius parameter ε values, and the greater the difference between the elevation performance of the cluster range corresponding to the cluster boundary j and the average elevation of other adjacent cluster ranges under different radius parameter ε values, the higher the topographic integrity of the cluster range. The larger the value, the greater the difference between the terrain shown by the cluster range enclosed by the cluster range boundary j and the surrounding terrain, and the more complete the representation of the terrain of a single region by the cluster range. Therefore, the completeness of the terrain of the cluster range enclosed by the cluster range boundary j can be obtained as follows:

[0067]

[0068] In the formula, Q j H represents the degree of topographic integrity of the cluster range enclosed by the cluster range boundary j; j,ε This represents the sum of the absolute values ​​of the differences between the mean elevation of the cluster bounded by cluster boundary j and the mean elevation of all other adjacent clusters under clustering operations with radius parameter ε; z represents the number of values ​​for radius parameter ε; n j n represents the total number of times the cluster boundary j appears under all radius parameter ε values. j ;norm represents the maximum-minimum linear normalization function.

[0069] Finally, according to the terrain integrity of the cluster range surrounded by each cluster range boundary, a subspace for the Kriging algorithm to perform the elevation interpolation calculation is obtained. Specifically, a completeness threshold is set, which can be 0.5; it is judged whether the terrain completeness is greater than the completeness threshold; if yes, i.e. Q j > 0.5, the cluster range surrounded by the cluster range boundary j is taken as the subspace for the Kriging algorithm to perform the elevation interpolation calculation. It should be noted that the corresponding subspace can be obtained under different radius parameters ε values.

[0070] S400: performing an interpolation operation on the subspace to obtain the projection position point and the initial elevation value thereof.

[0071] After obtaining the subspace for the Kriging algorithm to perform the elevation interpolation calculation, an interpolation operation is performed on the subspace to obtain the projection position point and the initial elevation value thereof. Further, it includes:

[0072] Firstly, an interpolation operation is performed on the subspace by the Kriging algorithm to obtain an interpolation elevation point, and this step is the prior art and will not be described here.

[0073] Then, since the interpolation result of the interpolation operation on the subspace by the Kriging algorithm is highly sensitive to the selection of the subspace, the difference between the subspaces will change the spatial structure expression of the sample points, affect the weight distribution and the model parameters, and finally lead to the difference between the projection position points of the same interpolation elevation point in different subspaces. Therefore, in some embodiments of the present application, the projection position point corresponding to the interpolation elevation point obtained by the Kriging interpolation algorithm for the subspace obtained when the cluster radius ε value is the smallest is taken as the priority analysis object.

[0074] Finally, for all subspaces m (subspaces obtained when the cluster radius ε value is the smallest) containing a single projection position point p, the elevation value of the projection position point p m corresponding to the interpolation elevation point obtained by each subspace m and closest to the projection position point p is obtained, which is taken as the initial elevation value of the projection position point p in each subspace m, denoted as h p,m .

[0075] S500: according to the terrain integrity, combining the projection area of the subspace, analyzing the reference value of different subspaces for calculating the elevation interpolation of the projection position point, and combining the initial elevation value, obtaining the reference elevation value of the projection position point.

[0076] According to the terrain integrity, combining the projection area of the subspace, analyzing the reference value of different subspaces for calculating the elevation interpolation of the projection position point, and combining the initial elevation value, obtaining the reference elevation value of the projection position point. Further, it includes:

[0077] First, according to the terrain integrity, combined with the projection area of the subspace, the reference value of different subspace to the interpolation height of the projection position point is obtained.

[0078] Specifically, the terrain integrity Q m of all subspace m is obtained m (The subspace is three-dimensional, the projection is the projection on the x-y plane, and the projection area can be directly obtained by the system). When the terrain integrity Q m of the subspace m is greater, and the projection area s m is greater, the existing sampling points in the subspace m are taken as the reference of Kriging algorithm, the interference of other terrain is smaller, and when the projection area s m is greater, the existing sampling points are more, and the reference to the height of the interpolation point is more meaningful. Thus, the initial height value h p,m of the projection position point p obtained by the subspace m is obtained. p,m The reference value of the projection position point p to the interpolation height is:

[0079] W m =norm(Q m ×s p,m )

[0080] In the formula, W p,m represents the reference value of the subspace m to the interpolation height of the projection position point p, that is, the initial height value h m of the projection position point p obtained by the subspace m to the interpolation height of the projection position point p; Q m represents the terrain integrity of the subspace m; s p,m represents the projection area of the subspace m; and norm represents the maximum and minimum linear normalization function.

[0081] When W p,m is greater, the reference value of the initial height value h m of the projection position point p obtained by the subspace m to the interpolation height of the projection position point p is greater.

[0082] Then, according to the reference value and the initial height value, the corresponding weighted height value of the projection position point in each subspace is obtained, the sum of all subspace containing the projection position point is summed to obtain the reference height value of the projection position point. Specifically, for a single projection position point p, the initial height value h p,m of the projection position point p obtained by all subspace m (terrain integrity Q p,m > 0.5) containing the projection position point p is obtained.The multiplication is carried out to obtain the weighted elevation value corresponding to the projection position point p in each subspace m; all the subspaces m containing the projection position point p are traversed, and the weighted elevation values are summed to obtain the reference elevation value of the single projection position point p, denoted as h p .

[0083] S600: The distance interval conditions of the elevation data in the same region in different subspaces are analyzed, the coordinate positions of the projection position points are adjusted in combination with the reference values, and a corrected DEM spatial model is obtained.

[0084] The projection position point of the interpolation elevation point obtained by the clustering range with the minimum clustering radius ε value is adjusted, and the corresponding relationship between the obtained reference elevation value and the projection position point may not be accurate enough. The reference value of the calculated projection position point interpolation elevation and the distance relationship between the projection position points can be used to reasonably correct the projection position of the adjusted interpolation elevation point (the reference elevation value of the projection position point).

[0085] Based on the above analysis, in some embodiments of the present application, the distance interval conditions of the elevation data in the same region in different subspaces are analyzed, the coordinate positions of the projection position points are adjusted in combination with the reference values, and a corrected DEM spatial model is obtained. Further comprising:

[0086] First, for the subspace m containing the projection position point p, the second Euclidean distance between the projection position point p and the corresponding projection position point p m (the projection position point p m corresponding to the interpolation elevation point closest to the projection position point p in the subspace m) is calculated, denoted as l p,m ; then the reference value W p,m of the calculated projection position point p interpolation elevation in different subspaces m containing the projection position point p is obtained; p,m When the reference value W p,m is larger, the second Euclidean distance l m is larger, the contribution value of the corresponding projection position point p m of the projection position point p in the subspace m to the projection position point p is larger, and the projection position point p should move more distance in the direction of the corresponding projection position point p p,m in the subspace m. Therefore, the second Euclidean distance l p,m corresponding to the subspace m is combined with the reference value W m to obtain the movement vector of the projection position point p to the corresponding projection position point p m in different subspaces m, wherein the direction of the movement vector is the direction of the corresponding projection position point p p,nThe product of the reference value W p,m .

[0087] Then, the moving vectors corresponding to all the subspaces m containing the projection position point p are added to obtain the moving reference total vector of the projection position point p, as shown in the following formula. Figure 2 Further, the coordinate position of the projection position point (the coordinate position of the interpolation elevation point corresponding to the projection position point p in the x-y plane) is adjusted according to the moving reference total vector. And the corrected DEM spatial model is obtained based on the adjusted coordinate position.

[0088] S700: Analyzing the difference between the water flow trend and the actual terrain in the corrected DEM spatial model, adjusting the reference elevation value of part of the projection position points, and obtaining a second corrected DEM spatial model.

[0089] The difference between the terrain in the corrected DEM spatial model obtained in step S600 and the actual terrain may cause the water flow path and the like obtained based on the DEM spatial model to be different from the actual situation. In order to make the terrain in the DEM spatial model more consistent with the actual terrain, the reference elevation value of different elevation points should be adjusted through the simulated hydrological information.

[0090] Based on the above analysis, in the embodiments of the present application, the reference elevation value of part of the projection position points is adjusted by analyzing the difference between the water flow trend and the actual terrain in the corrected DEM spatial model, and a second corrected DEM spatial model is obtained. Further, it includes: analyzing the area difference degree of the water area in the corrected DEM spatial model and the actual water area, and analyzing the flow rate difference degree of the water flow rate in the corrected DEM spatial model and the actual water flow rate, adjusting the reference elevation value of part of the projection position points, and obtaining a second corrected DEM spatial model. Wherein:

[0091] The area difference degree of the water area in the corrected DEM spatial model and the actual water area is analyzed, and the specific implementation is as follows: first, the contour lines are generated by using the corrected DEM spatial model (such as the Contour List tool), and the water area is converted into a polygon model according to the elevation step, the area of the model water area is calculated, and the model water area in the corrected DEM spatial model is obtained; at the same time, the area of the actual water area corresponding to the model water area is obtained (obtained through existing data records); finally, the difference between the model water area and the actual water area (model water area-actual water area) is calculated, and the area difference degree of the water area in the corrected DEM spatial model and the actual water area is obtained.

[0092] The flow velocity difference degree of the water flow in the corrected DEM spatial model and the actual water flow is analyzed, and the specific implementation manner is as follows: first, based on the corrected DEM spatial model, the water flow direction matrix is calculated by using the flow direction algorithm (such as D8 or single flow direction algorithm), and then the water flow accumulation matrix is generated in combination with the slope, and the model water flow velocity of the model water flow in the corrected DEM spatial model is obtained through the accumulation amount of the water flow accumulation matrix and the slope; at the same time, the velocity of the actual water flow corresponding to the model water flow is obtained (obtained through the existing data record); and the first Euclidean distance between the projection position point corresponding to the interpolated elevation point on the terrain of the model water flow and the center point of the model water area is obtained, wherein the model water flow is connected with the model water area. Finally, the difference between the model water flow velocity and the actual water flow velocity is calculated, and the first Euclidean distance is combined to obtain the flow velocity difference degree of the water flow in the corrected DEM spatial model and the actual water flow.

[0093] The area difference degree of the water area in the corrected DEM spatial model and the actual water area, and the flow velocity difference degree of the water flow in the corrected DEM spatial model and the actual water flow are obtained, and then the adjustment degree of the projection position point corresponding to the interpolated elevation point is calculated by combining the area difference degree and the flow velocity difference degree.

[0094] U k =tanh{exp(s d -S d )×exp[(V b -v b )×l k,b,d ])

[0095] In the formula, U k represents the elevation value adjustment degree of the projection position point k; s d represents the model water area area of the model water area d in the corrected DEM spatial model; S d represents the actual water area area of the actual water area corresponding to the model water area d in the corrected DEM spatial model; v b represents the model water flow velocity of the model water flow b in the corrected DEM spatial model; V b represents the actual water flow velocity of the actual water flow corresponding to the model water flow b in the corrected DEM spatial model; l k,b,d represents the first Euclidean distance between the projection position point k on the terrain of the model water flow b and the center point of the model water area d (the model water flow b is connected with the model water area d); tanh represents the hyperbolic tangent function, and exp represents the exponential function with the natural constant e as the base number.

[0096] s d -S drepresents the area difference degree of the water area in the modified DEM spatial model and the actual water area, when the difference is positive and larger, the elevation value of the interpolation elevation point near the water flow inlet in the DEM spatial model should be larger, so as to reduce the area of the model water area.(V b -v b )×l k,b,d represents the flow velocity difference degree of the water flow in the modified DEM spatial model and the actual water flow, when the flow velocity difference (V b -v b ) is positive and larger, in order to improve the flow velocity of the model water flow in the DEM spatial model, the elevation value of the interpolation elevation point far from the model water area position should be increased more, and further, when the elevation value adjustment degree U k is larger, the elevation value of the interpolation elevation point should be increased more.

[0097] Therefore, the elevation adjustment value of the projection position point k is obtained as follows:

[0098] h k '=h k ×(1+U k )

[0099] In the formula, h k ' represents the elevation adjustment value of the projection position point k; h k represents the reference elevation value of the projection position point k; and U k represents the elevation value adjustment degree of the projection position point k.

[0100] Similarly, the interpolation elevation points in the area where the model water area does not match the actual water area are corrected, and based on the elevation adjustment value of the projection position point, a secondary modified DEM spatial model is finally obtained.

[0101] By the above method, the interpolation elevation points at different positions are accurately obtained, and an accurate DEM spatial model is constructed. The obtained DEM spatial model and the data of a plurality of interpolation elevation points are stored in a database in correspondence.

[0102] The positions of different interpolation elevation points in the DEM spatial model are visually displayed in the form of a table by using the SQL query statement, as shown in the following table.

[0103] Interpolated elevation point Position / m 001 (25,45,36) 002 (25,67,36) … …

[0104] It should be noted that the above-mentioned embodiment order of the application is only for description, and does not represent the advantages and disadvantages of the embodiments. The processes depicted in the drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multi-task processing and parallel processing are also possible or may be advantageous.

[0105] The various embodiments described in this specification are presented by way of example, and each embodiment is not necessarily composed of all features described with respect to other embodiments.

Claims

1. A method for automatically extracting elevation points of a DEM spatial model, characterized in that, The method comprises: acquiring elevation data of sampling points; setting different clustering radii, clustering the elevation data to obtain clustering ranges corresponding to different clustering radii; analyzing occurrence of each clustering range boundary under different clustering radii and difference of elevation data on both sides of each clustering range boundary, judging terrain integrity of each clustering range, and obtaining multiple subspaces; performing interpolation operation on the subspaces to obtain projection position points and initial elevation values of the projection position points; according to the terrain integrity, combining projection areas of the subspaces, analyzing reference values of different subspaces for calculating interpolation elevation of the projection position points, and combining the initial elevation values, obtaining reference elevation values of the projection position points; analyzing distance interval conditions of elevation data in the same region in different subspaces, combining the reference values, adjusting coordinate positions of the projection position points, and obtaining a modified DEM spatial model; analyzing difference conditions of water flow trends and actual terrain in the modified DEM spatial model, adjusting the reference elevation values of part of the projection position points, and obtaining a second modified DEM spatial model.

2. The method of claim 1, wherein, The analyzing of occurrence of each clustering range boundary under different clustering radii and difference of elevation data on both sides of each clustering range boundary, and the judging of terrain integrity of each clustering range, comprises: counting total times of occurrence of each clustering range boundary under all clustering radii, and quantifying occurrence of each clustering range boundary under different clustering radii; under clustering radii containing the same clustering range boundary, calculating absolute value of difference between mean value of elevation data in a clustering range surrounded by the clustering range boundary and mean value of elevation data in an adjacent clustering range, and obtaining difference degree of elevation data on both sides of each clustering range boundary; combining the occurrence and the difference degree, obtaining terrain integrity degree of the clustering range surrounded by each clustering range boundary.

3. The method of claim 2, wherein, The obtaining of multiple subspaces comprises: setting an integrity threshold; judging whether the terrain integrity degree is greater than the integrity threshold; if yes, taking the clustering range surrounded by the clustering range boundary as a subspace for elevation interpolation calculation by a Kriging algorithm, and obtaining multiple subspaces.

4. The method of claim 1, wherein, The performing of interpolation operation on the subspaces to obtain projection position points and initial elevation values of the projection position points comprises: performing interpolation operation on the subspaces by a Kriging algorithm to obtain interpolation elevation points; for the subspace obtained when the clustering radius value is the smallest, obtaining a projection position point corresponding to the interpolation elevation point by a Kriging interpolation algorithm; for all subspaces containing a single projection position point, obtaining elevation values of the interpolation elevation points closest to the projection position point in each subspace, and obtaining initial elevation values of the projection position point in each subspace.

5. The method of claim 1, wherein, The obtaining of reference elevation values of the projection position points according to the terrain integrity, combining projection areas of the subspaces, analyzing reference values of different subspaces for calculating interpolation elevation of the projection position points, and combining the initial elevation values, comprises: according to the terrain integrity, combining projection areas of the subspaces, obtaining reference values of different subspaces for calculating interpolation elevation of the projection position points; According to the reference value and the initial elevation value, a corresponding weighted elevation value of the projection position point in each subspace is obtained, summation of all the subspace containing the projection position point is obtained, and a reference elevation value of the projection position point is obtained.

6. The method of claim 1, wherein, The distance interval of each elevation data in the same region in different subspaces is analyzed, and the coordinate position of the projection position point is adjusted in combination with the reference value to obtain a corrected DEM spatial model, including: A second Euclidean distance between the projection position point and the corresponding projection position point in different subspaces is calculated, and a moving vector of the projection position point to the corresponding projection position point in different subspaces is obtained in combination with the reference value; The moving vectors corresponding to all the subspaces are summed to obtain a moving reference total vector of the projection position point; According to the moving reference total vector, the coordinate position of the projection position point is adjusted to obtain a corrected DEM spatial model.

7. The method of claim 1, wherein, The difference between the water flow trend in the corrected DEM spatial model and the actual terrain is analyzed, including: The area difference degree of the water area in the corrected DEM spatial model and the actual water area, and the flow rate difference degree of the water flow rate in the corrected DEM spatial model and the actual water flow rate are analyzed.

8. The method of claim 7, wherein, The area difference degree of the water area in the corrected DEM spatial model and the actual water area, including: An isoline is generated using the corrected DEM spatial model, and a polygon model water area is converted according to an elevation step, the area of the model water area is calculated, and the model water area in the corrected DEM spatial model is obtained. The actual water area corresponding to the model water area is obtained. The difference between the model water area and the actual water area is calculated to obtain the area difference degree of the water area in the corrected DEM spatial model and the actual water area.

9. The method of claim 8, wherein, The flow rate difference degree of the water flow rate in the corrected DEM spatial model and the actual water flow rate, including: Based on the corrected DEM spatial model, a flow direction matrix is calculated using a flow direction algorithm, and a flow accumulation matrix is generated in combination with a slope, and the model water flow rate of the model water flow in the corrected DEM spatial model is obtained through the accumulation amount of the flow accumulation matrix and the slope. The actual water flow rate corresponding to the model water flow is obtained. The first Euclidean distance between the interpolation elevation point on the terrain of the model water flow and the center point of the model water area is obtained, wherein the model water flow and the model water area are connected. The difference between the model water flow rate and the actual water flow rate is calculated, and the flow rate difference degree of the water flow rate in the corrected DEM spatial model and the actual water flow rate is obtained in combination with the first Euclidean distance.

10. The method of claim 1, wherein, The arrangement method of the sampling points is: The sampling points are arranged by a feature point priority method, which requires that the sampling points are arranged every 5-10 meters on ridge and valley lines, and are arranged every 1-3 meters on slope mutation places, and are arranged every 20-50 meters on flat areas.

Images