A DEM upscaling method considering multi-level expression of terrain structure lines

Abstract

The present invention discloses a DEM upscaling method that takes into account the multi-level expression of terrain structure lines, and relates to the field of digital elevation model construction and terrain analysis technology. The technical solution includes the following steps: S1, generating valley data under different runoff accumulation thresholds; S2, generating upscaled DEM initial values ​​from original DEM data; S3, constructing local morphological constraint equations in the upscaled DEM grid; S4, constructing the accuracy control equation of the upscaled DEM with different valley data; S5, combining the equations constructed in steps S3 and S4 to form an equation group, and solving the equation group to obtain the grid elevation value of the upscaled DEM. Different valley data participate in the calculation to obtain DEM upscaling results at different levels. The present invention is characterized in that the morphological accuracy of the DEM is ensured by constructing local morphological constraint equations, and the DEM upscaling results of the effective expression of multi-level terrain structure lines are achieved by means of valleys extracted under different runoff accumulation thresholds.

Description

Technical Field

[0001] The invention belongs to the technical field of digital elevation model construction and terrain analysis, and relates to a DEM upscaling method taking into account multi-level expression of terrain structure lines. Background Art

[0002] Digital Elevation Model (DEM) is a digital abstract expression of the spatial distribution of the earth's surface elevation. It is the basis of terrain analysis and the core content of national basic geographic data. With the rapid development of surface observation technology and computer hardware and software technology, the terrain data used to construct and produce DEM data is now multi-source, multi-scale and massive. However, it takes a lot of time and manpower to generate DEMs of different scales based on original terrain data. Therefore, it is a useful supplement to the construction of multi-scale DEM databases to obtain multi-scale DEM data based on a scale of DEM through scale transformation, especially DEM upscaling. Because of its relatively complete theoretical basis, it has become the main way to obtain large-scale DEM.

[0003] The so-called DEM upscaling is to convert a high-resolution DEM into a low-resolution DEM, or to change the size of the DEM grid cell from small to large. DEM upscaling and DEM terrain simplification are similar in concept, so many studies confuse the two. However, from the direct goal of the operation, the two concepts are different. DEM upscaling focuses on the change of grid cell size, while DEM terrain simplification simplifies the terrain expressed by the DEM. The DEM resolution can be changed or not. Therefore, the more accurate relationship between the two is that DEM upscaling is an important way to simplify DEM terrain. Because the size of the DEM grid cell changes, smaller terrain ups and downs cannot be effectively expressed, thus achieving the goal of terrain simplification.

[0004] Considering the close relationship between DEM upscaling and DEM terrain simplification, almost all current research on DEM upscaling is related to DEM terrain simplification, that is, DEM upscaling is regarded as an important way to achieve DEM terrain simplification. There are many studies related to terrain simplification based on DEM upscaling, and different scholars have proposed many upscaling methods. The more commonly used methods include interpolation based on data resampling, global filtering, and structured upscaling methods. Among the above methods, interpolation based on data resampling may destroy the high-order sequence relationship of DEM during the resampling process, and how to resample is also an issue that requires special attention. Simple regular sampling may destroy the expression of terrain skeleton feature lines, and some sampling strategies that take into account the expression of terrain structure lines will bring additional workload; the global filtering method is mainly through window analysis, with average filtering being the most commonly used, which can weaken the terrain mutation within the analysis window. The main problem of this method is that it performs filtering operations based on window analysis on the complete data area. The typical feature of this operation mode is the uniform filtering of terrains of different complexity. The result may be insufficient simplification of complex terrain areas and over-simplification of simple terrain areas, which may also cause problems such as valley line and ridge line offset. The structured DEM upscaling method incorporates typical terrain feature elements such as valleys and ridges, and considers the scale changes of terrain feature elements in the process of terrain simplification. The terrain simplification result obtained by this method is better in feature preservation, but its disadvantage is that it requires auxiliary data and the algorithm is relatively complex.

[0005] From the above analysis, it can be seen that the existing DEM upscaling methods have their own advantages and disadvantages, and most of the current DEM upscaling methods focus on the calculation of single grid cell values, but pay insufficient attention to the expression effect of the terrain after upscaling. In order to solve the above defects, the method of the present invention proposes a DEM upscaling method that takes into account the multi-level expression of terrain structure lines. The design idea of ​​this method is based on the following analysis: First, the terrain structure line is the control line of the terrain morphology, and the expression of the terrain structure line by the DEM upscaling result should be a gradual change. For example, a valley should gradually disappear in different DEM upscaling results; secondly, after DEM upscaling, the terrain undulations and mutations on the small scale are eliminated, so the terrain expressed by the DEM upscaling result should have better smooth characteristics; finally, no matter how much the DEM is upscaled, the elevation sequence relationship of different spatial positions recorded by it should remain unchanged, otherwise the DEM upscaling result will lose its basic application value.

[0006] The above-mentioned analyses are also the basic goals to be achieved by the method of the present invention in the process of DEM upscaling. To this end, the present invention designs a technical process for DEM upscaling: multiple morphological and elevation sequence relationship constraint equations are constructed on the local surface, and then the original DEM data of the valley lines at different scales are obtained as precision control points, and then the two are combined to calculate the DEM upscaling results. The DEM upscaling results of the technical process of the present invention can well express the changes in the terrain expression characteristics brought about by the DEM scale change, and the DEM has very good smoothness and correct elevation sequence relationship. The technology of the present invention is expected to provide important technical method support for the construction of a multi-scale DEM database. Summary of the invention

[0007] In view of the above problems, the present invention designs a DEM upscaling method which takes into account the multi-level expression of terrain structure lines. By constructing morphological constraint equations between local grid cells, the smoothness of the upscaled DEM and the orderly elevation relationship between different grid cells are guaranteed. In the upscaling process, terrain structure line information of different levels is added so that the upscaled DEM can well maintain the original basic undulating morphology. In addition, a sequence DEM upscaling data product is quickly constructed according to the different requirements of DEM data in different fields.

[0008] The technical solution adopted by the present invention is as follows:

[0009] A DEM upscaling method taking into account the multi-level expression of terrain structure lines includes the following steps: first, based on DEM data, the confluence process simulation related operations are performed to obtain the confluence accumulation raster data of the study area, and a threshold sequence is set to sequentially extract the raster values ​​above different threshold values ​​and convert them into vector point data; then, an upscaled result raster data set is created, and the elevation average values ​​of the original DEM data raster cells in the upscaled result raster cells are counted, and the value is assigned to the raster cells of the upscaled result raster data as the initial elevation value; then, for each raster cell of the upscaled result raster data, the surrounding raster cells are selected to participate in the construction of a local morphological constraint equation group, and the point data converted from the raster above a certain confluence accumulation threshold value are assigned the elevation value of the original DEM data at the corresponding position as the accuracy control condition, and a large linear equation group is formed by combining them; finally, the pre-processed conjugate gradient method is used to iteratively calculate and solve the large linear equation group to obtain the elevation value of each grid cell after upscaling, thereby obtaining the upscaled DEM data.

[0010] Furthermore, the method comprises the following steps:

[0011] S1, after the processes of depression filling and flow direction calculation, the runoff accumulation raster data of the study area is calculated, and an increasing runoff accumulation threshold sequence is set according to the upscaling data requirements, such as 500, 1000, 1500, etc.;

[0012] S2, taking each runoff accumulation threshold in turn, extracting the part of the runoff accumulation raster data that is greater than the threshold, converting the selected raster into vector point data, and then assigning the elevation value of the corresponding spatial position in the original DEM data to each vector point, and finally obtaining the vector point data set under different runoff accumulation thresholds;

[0013] S3, according to the target grid size of the upscaling result, the grid size, number of rows and columns, coordinates and other information of the upscaling result raster data are calculated according to the grid size, number of rows and columns, coordinates and other information of the original DEM data;

[0014] S4, for each grid cell of the upscaled result grid data, the average value of the grid cell elevation values ​​of the original DEM data falling within the cell is calculated as the initial value of the grid cell corresponding to the upscaled result raster data;

[0015] S5, for the upscaled grid data, traverse each grid cell in turn. If the current grid cell is not on the boundary, take the eight surrounding grid cells respectively, and combine the current grid cell to construct a local morphological constraint equation group;

[0016] S6, superimposing the vector point data under different runoff accumulation thresholds obtained in step S2 with the upscaled target grid data according to the position, taking the vector point data as the elevation information source, using the inverse distance weighted interpolation method to calculate the elevation value of the grid cell in the upscaled target grid data where the vector point falls, and constructing a set of precision control constraint equations;

[0017] S7, combining the equations constructed in step S5 and step S6 to form a large linear equation system, using the pre-conditioned conjugate gradient iteration method to solve the equation system, and outputting the solution vector of the equation system in a raster data format to obtain the final DEM upscaling result;

[0018] S8, reselect vector point data under different runoff accumulation thresholds, and repeat steps S6-S7 until all vector points are completed to obtain DEM upscaling results at different levels.

[0019] Furthermore, the specific process of step S5 is as follows:

[0020] S51, traverse each grid cell of the upscaled result grid data in sequence, if the current grid cell is in the first row, or the last row, or the first column, or the last column, then the grid is not processed, otherwise, execute steps S52-S54 in sequence;

[0021] S52, taking the current grid and its left grid and right grid, constructing a first morphological constraint equation, wherein the grid value is expressed as an initial value plus an unknown increment;

[0022] S53, taking the current grid and its upper and lower grids, constructing a second morphological constraint equation, wherein the grid value is expressed as an initial value plus an unknown increment;

[0023] S54, taking the current grid and its upper left grid, upper right grid, lower left grid and lower right grid, and constructing a third morphological constraint equation, wherein the grid value is expressed as the initial value plus an unknown increment;

[0024] S55, process other grid cells in turn, construct three morphological constraint equations for each non-boundary grid cell, and after the traversal is completed, combine the first morphological constraint equation established for each grid cell into the first large linear equation group, combine the second morphological constraint equation established for each grid cell into the second large linear equation group, and combine the third morphological constraint equation established for each grid cell into the third large linear equation group.

[0025] Furthermore, the specific process of step S6 is as follows:

[0026] S61, selecting a vector point data corresponding to a convergence accumulation threshold as precision control point data;

[0027] S62, performing a spatial overlay operation on the selected vector point data and the created upscaling result raster data;

[0028] S63, for the upscaled result grid data, searching for the grid cells where the vector points fall;

[0029] S64, for the grid cell where the vector point falls, the elevation value of the center point of the grid cell is calculated by using the inverse distance weighted interpolation method using the elevation value of the vector point in the grid;

[0030] S65, with the grid of calculated elevation values, establishes the accuracy constraint equations to form a large linear equation system.

[0031] Technical features and beneficial effects of the present invention:

[0032] (1) The DEM upscaling method proposed in this paper takes into account the multi-level expression of terrain structure lines. It ensures the morphological accuracy of the upscaled DEM by constructing a local morphological constraint equation, and uses the valley lines extracted under different runoff accumulation thresholds to ensure that the terrain structure lines can be correctly expressed after the DEM is upscaled.

[0033] (2) The DEM upscaling method proposed in the present invention takes into account the multi-level expression of terrain structure lines. By setting the threshold parameters of sequence changes for the runoff accumulation raster data, it is possible to achieve DEM upscaling products with multi-level expression of terrain structure lines, providing more data sources for the demand for DEM scale in different fields. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 It is a schematic diagram of the process of the present invention.

[0035] Figure 2 Schematic diagram of the positional relationship between the original DEM grid and the upscaled DEM grid in the method of the present invention.

[0036] Figure 3 Schematic diagram of the 3×3 analysis window in the method of the present invention.

[0037] Figure 4 It is a position relationship diagram of the original DEM grid center point and the upscaled DEM grid center point on the confluence path in the method of the present invention.

[0038] Figure 5 This is the DEM upscaling result corresponding to the runoff accumulation threshold parameter of 500 in the embodiment of the present invention.

[0039] Figure 6 This is the DEM upscaling result corresponding to the runoff accumulation threshold parameter of 5000 in the embodiment of the present invention.

[0040] Figure 7 This is the DEM upscaling result corresponding to the runoff accumulation threshold parameter of 15000 in the embodiment of the present invention.

[0041] Figure 8 This is the DEM upscaling result generated by the neighborhood sampling method in an embodiment of the present invention.

[0042] Fig. 9 It is the contour line extracted from the upscaling result corresponding to the confluence accumulation threshold parameter 500 in the embodiment of the present invention.

[0043] Fig.10 It is the contour line extracted from the upscaling result corresponding to the confluence accumulation threshold parameter 5000 in the embodiment of the present invention.

[0044] Fig.11 It is the contour line extracted from the upscaling result corresponding to the confluence accumulation threshold parameter 15000 in the embodiment of the present invention.

[0045] Fig.12 It is the contour line extracted by using the upscaling result of the neighborhood sampling method in the embodiment of the present invention. DETAILED DESCRIPTION

[0046] The principles and features of the present invention are described below in conjunction with the accompanying drawings. The examples given are only used to explain the present invention and are not used to limit the scope of the present invention.

[0047] like Figure 1 As shown, the specific steps of the DEM upscaling method taking into account the multi-level expression of terrain structure lines of the present invention are as follows:

[0048] Step 1: Read the DEM data of the study area, and implement the processes of depression filling, flow direction calculation, and runoff accumulation calculation in sequence to obtain the runoff accumulation value of each grid;

[0049] Step 2: Set a runoff accumulation threshold sequence, such as 1000, 2000, 3000, etc., and then extract grids with values ​​greater than the runoff accumulation threshold from the runoff accumulation grid data;

[0050] Step 3: Convert the extracted raster into vector point data, and read the elevation value from the same spatial position in the original DEM data and assign it to each point;

[0051] Step 4: Determine the upscaled DEM grid size based on the original DEM data grid size, and then create raster data to store the upscaled results;

[0052] Step 5: Overlay the original DEM data with the upscaled DEM data, and calculate the average elevation of the grid cells of the original DEM data that fall into each grid cell in the upscaled DEM ( Figure 2 ), assign values ​​to the upscaled DEM grid cells as initial values;

[0053] Step 6: Starting from the first non-boundary grid cell of the upscaled DEM grid data as the current processing unit, take the surrounding grid cells to form a 3×3 analysis window ( Figure 3 ), where the grid cell P 4 , P 5 , P 6 Construct the morphological constraint equation (1), the grid cell P 2 , P 5 , P 8 Construct the morphological constraint equation (2), the grid cell P 1 , P 3 , P 5 , P 7 , P 8 Construct the morphological constraint equation (3);

[0054]

[0055]

[0056]

[0057] In the above formula, f xx 、f yy 、f xy They represent the second-order partial derivative in the X direction, the second-order partial derivative in the Y direction, and the mixed second-order partial derivative in a certain grid cell, respectively. i Indicates the initial elevation value of the corresponding grid cell, Δfi It represents the deviation between the initial elevation value of the corresponding grid cell and the elevation value after upscaling, and h represents the grid cell size.

[0058] Step 7: Traverse all non-boundary grid cells and combine the three equations established for each grid cell into a system of equations:

[0059] AX=d (4)

[0060] BX=q (5)

[0061] CX=p (6)

[0062] In the above formula, A, B, and C represent the coefficient matrices of the three equations respectively, X represents the vector composed of unknown variables, and d, q, and p represent the vectors composed of constant terms in the three equations.

[0063] Step 8: Select the vector point data corresponding to a runoff accumulation threshold generated in step 2-3, and overlay it with the upscaled DEM grid data. For the grid cells where the vector points fall, use the inverse distance weighted method to calculate the elevation value of the grid cell with the elevation value of the vector point falling into a certain grid cell ( Figure 4 );

[0064] Step 9: For all grid cells whose elevation values ​​are recalculated, establish the precision control point equation in the form of (7), and combine them into the equation system (8);

[0065] fm i +Δf i =k i (7)

[0066] SX=K (8)

[0067] In the above formula, fm i Indicates the initial elevation value of a grid cell, Δf i Indicates the deviation between the initial elevation value of a grid cell and the elevation value of the upscaling result, k i It represents the elevation value of the upscaling result of a grid cell, S represents the coefficient matrix of the equation system, X represents the vector composed of unknown variables, and K represents the vector composed of constant terms in the equation system.

[0068] Step 10: Combine the precision control point equation group (8) established in step 9 and the equation groups (4), (5), and (6) established in step 7 to form a new linear equation group according to the least squares criterion:

[0069]

[0070] Convert the above formula into a system of equations:

[0071] WX=V (10)

[0072] in:

[0073] W=A T A+B T B+C T C+S T S

[0074] V=A T d+B T q+C T p+S T k

[0075] Step 11: Solve the equation group (10), and assign the solution vectors of the equation group to the upscaled grid cells in turn, so as to obtain a DEM upscaling result under a runoff accumulation threshold;

[0076] Step 12: Select the vector point data corresponding to the remaining runoff accumulation thresholds, repeat steps 8 to 11, and calculate the DEM upscaling results under other runoff accumulation thresholds.

[0077] Example

[0078] Existing DEM upscaling methods are either based on a resampling strategy at a certain distance interval, or combine terrain structure lines of different levels to generate DEMs of different scales. They cannot correctly express the ordered relationship of DEM elevation at different scales, nor can they express the gradual characteristics of terrain structure information caused by the change of DEM scale. In order to verify that the DEM upscaling results implemented based on the method of the present invention are more reasonable, the embodiment selected a small watershed in a certain area of ​​the Loess Plateau in northern Shaanxi as the experimental area. The gullies in this area are well developed and the terrain undulations are well layered. It is an ideal experimental area to verify the effect of DEM upscaling. The implementation steps are as follows:

[0079] Step 1: Read the DEM data of the experimental area, whose grid size is 5 meters, and perform the processes of filling depressions, flow direction calculation, and runoff accumulation calculation in sequence;

[0080] Step 2: Set the runoff accumulation threshold sequence parameters to 500, 5000, and 15000, and extract the grids whose runoff accumulation values ​​are greater than the threshold parameters in sequence;

[0081] Step 3: Convert the selected grid cells into points and assign the elevation values ​​of the original DEM;

[0082] Step 4: Generate DEM upscaling result raster data by calculating the average value. In this example, the grid size of the upscaling result is set to 25 meters.

[0083] Step 5: Construct the morphological constraint equation, and use the point data obtained in step 3 as the accuracy control condition to establish a large linear equation system. Solve the equation system to obtain the upscaled DEM result, such as Figure 5 , Figure 6 , Figure 7 As shown;

[0084] Step 6: Based on ArcGIS software, execute the upscaling method integrated in the software to obtain the upscaling result DEM as the control result, such as Figure 8 As shown;

[0085] Step 7: Using ArcGIS software, extract contour lines using different upscaled DEMs, with a contour line spacing of 10 meters. Compare and analyze the differences between the upscaling results of the method of the present invention and the upscaling results of ArcGIS software from the perspective of contour line fitting analysis. Fig. 9 , Fig.10 , Fig.11 , Fig.12 shown.

[0086] By comparing the upscaled DEM corresponding to different runoff accumulation thresholds generated by the method of the present invention and the upscaled DEM generated by the control method, as well as the contour lines extracted from different DEM results, it can be seen that the advantages of the method of the present invention are reflected in two aspects: first, the terrain continuity of the upscaled DEM result is good, especially the valleys reflecting the terrain skeleton characteristics are continuous; second, by setting different runoff accumulation threshold parameters, the method of the present invention can obtain DEM sequence results with different simplification degrees, which can meet the needs of different needs for different DEMs. In contrast, the control method can only generate one type of upscaled DEM data, and the terrain expressed by its results presents discontinuous characteristics, especially the phenomenon of discontinuous distribution of some valleys, and the surface morphology expression cannot meet the relevant application needs.

Claims

1. A DEM upscaling method that takes into account the multi-level expression of terrain structure lines, characterized in that: The steps include: S1, after the depression filling and flow direction calculation process, the runoff accumulation raster data of the study area is calculated, and an increasing runoff accumulation threshold sequence is set according to the upscaling data requirements; S2, taking each runoff accumulation threshold in turn, extracting the part of the runoff accumulation raster data that is greater than the threshold, converting the selected raster into vector point data, and then assigning the elevation value of the corresponding spatial position in the original DEM data to each vector point, and finally obtaining the vector point data set under different runoff accumulation thresholds; S3, according to the target grid size of the upscaling result, the number of rows and columns, and the coordinate information of the upscaling result raster data are calculated according to the grid size, number of rows and columns, and coordinate information of the original DEM data; S4, for each grid cell of the upscaled result grid data, the average value of the grid cell elevation values ​​of the original DEM data falling within the cell is calculated as the initial value of the grid cell corresponding to the upscaled result raster data; S5, for the upscaled grid data, traverse each grid cell in turn. If the current grid cell is not on the boundary, take the eight surrounding grid cells respectively, and combine the current grid cell to construct a local morphological constraint equation group; The specific process of S6 is as follows: S61, selecting a vector point data corresponding to a convergence accumulation threshold as precision control point data; S62, performing a spatial overlay operation on the selected vector point data and the created upscaling result raster data; S63, for the upscaled result grid data, searching for the grid cells where the vector points fall; S64, for the grid cell where the vector point falls, the elevation value of the center point of the grid cell is calculated by using the inverse distance weighted interpolation method using the elevation value of the vector point in the grid; S65, using the grid of calculated elevation values ​​to establish the accuracy constraint equation, forming a large linear equation system; S7, combining the equations constructed in step S5 and step S6 to form a large linear equation system, using the pre-conditioned conjugate gradient iteration method to solve the equation system, and outputting the solution vector of the equation system in a raster data format to obtain the final DEM upscaling result; S8, reselect vector point data under different runoff accumulation thresholds, and repeat steps S6-S7 until all vector points are completed to obtain DEM upscaling results at different levels.

2. A DEM upscaling method taking into account the multi-level expression of terrain structure lines according to claim 1, characterized in that: The specific process of step S5 is as follows: S51, traverse each grid cell of the upscaled result grid data in sequence, if the current grid cell is in the first row, or the last row, or the first column, or the last column, then the grid is not processed, otherwise, execute steps S52-S54 in sequence; S52, taking the current grid and its left grid and right grid, constructing a first morphological constraint equation, wherein the grid value is expressed as an initial value plus an unknown increment; S53, taking the current grid and its upper and lower grids, and constructing a second morphological constraint equation, wherein the grid value is expressed as an initial value plus an unknown increment; S54, taking the current grid and its upper left grid, upper right grid, lower left grid and lower right grid, and constructing a third morphological constraint equation, wherein the grid value is expressed as the initial value plus an unknown increment; S55, process other grid cells in turn, construct three morphological constraint equations for each non-boundary grid cell, and after the traversal is completed, combine the first morphological constraint equation established for each grid cell into the first large linear equation group, combine the second morphological constraint equation established for each grid cell into the second large linear equation group, and combine the third morphological constraint equation established for each grid cell into the third large linear equation group.

Images