Analysis and application of DEM

Fitting of topographic surface

The most basic application of DEM is to find the elevation of any point within the DEM range, and then perform terrain attribute analysis on this basis. Since the elevation of a finite number of grid points is known, it is possible to fit a terrain surface using these grid point elevations to derive the elevation of any point within the area. The surface fitting method can be regarded as a special case of spatial interpolation of a known regular grid point data, the distance reciprocal weighted average method, Kriging interpolation method, spline function and other interpolation methods can be used.

Stereoscopic perspective

Drawing perspective stereogram from digital elevation model is an extremely important application of DEM. Perspective stereogram can better reflect the three-dimensional shape of terrain and is very intuitive. Compared with using contour lines to represent topographic morphology, which has its own unique advantages and is closer to people’s visual perception. Especially with the enhancement of computer graphics processing and the development of screen display system, the production of stereo graphics has greater flexibility, people can make different stereo displays for the same terrain shape according to different needs. For example, local enlargement, change the magnification of elevation Z to exaggerate the three-dimensional shape; change the position of the viewpoint to observe from different angles, and even make the three-dimensional figure rotate, so that people can better study the spatial shape of the terrain.

From a spatial three-dimensional digital elevation model to a planar two-dimensional perspective, the essence is a perspective transformation. Considering “Viewpoint” as the “photography center”, you can directly apply the collinear equation to calculate the “pixel” coordinates (X,Y)from the object point (X,Y,Z). Another problem in the perspective is the problem of “blanking”, which deals with the problem of the foreground.

By adjusting the values of viewpoints, perspectives and other parameters, animation can be made from different perspectives with different shapes drawn from different directions and distances. When the computer speed is high enough, animated DTM perspective can be generated in real time.

Visibility analysis

Insight analysis has a wide range of applications. A typical example is the observation of the position of the post, the position of the observation post should be located to monitor an area of interest and the line of sight cannot be blocked by the terrain obviously. This is the typical point-to-area viewing problem in the analysis. Similar problems are the setting of fire monitoring points in the forest, the setting of the wireless tower, and so on. Sometimes it is also possible to analyze the invisible area, for example, when the low-altitude reconnaissance aircraft is flying, it should avoid the capture of the enemy radar as much as possible, the flight obviously chooses the radar blind spot flight. The problem of viewing can be divided into five categories [Lee, J. (1991)]:

1. Find out the visible area of a terrain by knowing one or a group of observation points.

2. To observe all the topographic surfaces of an area, the minimum number of observation points is calculated.

3. On the premise of a certain number of observation points, the maximum observation area can be calculated.

4. Constructing observation towers at minimum cost requires all areas to be visible.

5. On the premise of a given construction cost, the maximum visible area is obtained.

According to the different output dimension of the problem, the pass-through can be divided into point view, line view and plane view. Point perspective refers to the visibility problem between calculating viewpoints and points to be determined; line perspective refers to the problem of calculating the visual field of a given viewpoint; area perspective refers to the problem of calculating the set of terrain surface regions with known viewpoints and visible viewpoints. The methods based on grid DEM model and DEM based on TIN model are quite different.

1) Point-to-point communication

Based on the grid DEM problem, in order to simplify the problem, the grid point can be used as the calculation unit. Such a point-to-point view problem is simplified as the intersection of a discrete space line and a terrain profile line. (Figure 9-13)

The coordinates of the known viewpoint V are (x_0, y_0, z_0) and the coordinates of point P (x_1, y_1, z_1). DEM is a two-dimensional array Z [M] [N], then V is (m_0, n_0, Z [m_0, n_0]), P is (m_1, n_1, Z [m_1, n_1]. The calculation process is as follows:

(1.1) Using the Bresenham line algorithm, generate a set of projection line points {x , y} of V to P, K=||{x , y}||, and get a set of straight points {x , y } Corresponding elevation data {Z[k], (k=1,…K-1)}, thus forming a DEM profile curve from V to P.

(1.2) With the projection line from V to P as the X-axis and the projection point of V as the origin, the linear equation of the line of sight in the X-Z coordinate system is obtained.

(0<k<K)

`` K` is the number of discrete points on the projection line from V to `P’.

(1.3) Compare the values of the corresponding elements in the array H[k] and the array Z[k], if

If Z[k]>H[k] exists, then V and P are not visible, otherwise they are visible.

2) Point-to-line view

Point-to-line vision is actually the field of view for finding points. It should be noted that the points on the surface of any terrain outside the horizon are invisible, but the points in the horizon may or may not be visible. The algorithm based on grid DEM point-to-line is as follows:

(2.1) Let P point be a point moving clockwise along the edge of DEM data, which is similar to the calculation point-to-point perspective, the point set {x, y} on the projection line from viewpoint to P point is obtained, and the corresponding topographic profile {x, y, Z(x, y)} is obtained.

(2.2) Calculate the angle between the viewpoint and each Z-axis:

(2.3) Find. The corresponding point is a point of view horizon.

(2.4) Move P-point and repeat the above process until P-point returns to its original position and the algorithm ends.

3) Point-to-area visibility

Point-to-area algorithm is an extension of point-to-point algorithm. As with point-to-line vision, point P moves clockwise along the data edge. Check point by point whether the point on the line from viewpoint to point P is visible. An improved algorithm is that the occlusion point from viewpoint to point P is most likely the point with the highest elevation on the topographic section line. Therefore, the points on the section line can be sorted according to the elevation value, and each point after sorting can be checked in descending order whether it is visible. As long as one point does not satisfy the condition of visibility, the rest points will not be checked. The essence of peer-to-area vision is still peer-to-peer vision, only adding the sorting process.

Watershed feature geomorphology extraction and terrain automatic segmentation

Terrain factors are important factors affecting the process of watershed geomorphology, hydrology, biology, etc., the spatial distribution characteristics of terrain attributes have always been an important indicator used by people to describe these spatial process changes. High-precision DEM data and high-resolution, hyperspectral, multi-period remote sensing images provide an increasingly rich source of data for quantitative description of watershed spatial change processes, and people understand the spatial variation mechanisms of process geomorphology, hydrology and biology, constantly deepening, it can be said that human beings have entered an era of “space simulation.” Automatic extraction of watershed geomorphology based on DEM data and automatic segmentation of watershed terrain are the basic techniques for basin space simulation.

The automatic extraction of watershed feature geomorphology and automatic terrain segmentation based on grid DEM mainly includes two aspects: 1) The definition of watershed geomorphological structure, defining the characteristic landform that can reflect the structure of the watershed, and establishing the micro-geomorphology corresponding to the DEM of the grid. 2) Automatic extraction of feature geomorphology and automatic terrain segmentation algorithm, grid DEM data is some discrete elevation point data, and each data itself does not reflect the actual surface complexity, in order to obtain the watershed geomorphological structure from the grid DEM data, a clear watershed geomorphic structure model must be adopted, and then an automatic extraction algorithm is designed for the structural model.

1) Watershed structure definition

You can use a tree diagram with roots to describe the watershed structure [Shreve], which is currently used by most algorithms. In this structure, there are mainly three parts, namely, a node set, a boundary set, and a confluence zone set. As shown in Figure 9-14.

Its specific content includes several concepts:

1. Valley line segment: a line segment with confluence zones on both sides;

2. Water dividing section: a section with water dividing areas on both sides;

3. Valley junction: the intersection of two or more Valley lines;

4. The junction of the dividing line: the intersection of two or more dividing lines;

5. Valley source: the upstream starting point of the valley;

6. Source point of watershed: intersection point between watershed and watershed boundary;

7. Internal confluence area: the confluence area whose boundary does not include part of the boundary of the basin;

8. External confluence area: the boundary of confluence area includes the confluence area of some basin boundary.

The valley nodes and the valley source points together form a valley node set, and all the valley segments form a valley group to form a valley network; all the water lines form a set of watershed segments to form a watershed network. The valley segment set and the watershed segment set together divide the watershed into a sink zone set.

The gully section is the smallest gully unit, and the gully section can be divided into an inner gully section and an outer gully section. The inner gully section connects the two gully nodes, and the outer gully section connects a gully node and a valley source point. Similarly, the watershed segment is the smallest watershed unit and is also divided into an inner partial waterline segment and an outer partial waterline segment. The inner partial water line segment connects the two water line nodes, and the outer partial water line segment connects one water line node and one water line source point.

Each valley segment in the confluence network has a confluence area that is controlled by a watershed set. The outer gully section has an external confluence zone, and the inner gully section has two internal confluence zones distributed on either side of the inner gully section, the entire watershed is divided into sub-basins, each sub-basin is like a “leaf” on the tree.

2) Automatic extraction of watershed feature geomorphology and automatic terrain segmentation

Definition and extraction of feature geomorphology: According to the relationship between the elevation of the grid point and the surrounding elevation value, the grid points are divided into slopes, depressions, watersheds, valleys, terraces and saddles. First calculate the elevation difference between the center point and the eight neighbor points, then sort the elevation difference, and then assign a feature code to the center point grid according to the characteristics of the elevation difference sequence. Then, through the combination of a series of feature codes, the grid points are divided into known feature landform categories by pattern recognition.

Ridge line and valley line extraction: Automatic detection of ridge lines and valley lines is actually an automatic search for pits and bumps. The simpler operator is the local operator of 2*2. The operator is slid in the DEM data to compare the elevation of each grid point with the adjacent grid points on the rows and columns, and to identify the grid points where the elevation is the smallest (probing the valley line) or the elevation is the largest (detecting the ridge line). After calculating the entire DEM data, the remaining unmarked grid points are grid points on the ridgeline or valley line.

Automatic segmentation of watershed terrain: The goal of automatic segmentation of watershed terrain is to divide the entire watershed into sub-sink areas. Most algorithms use 3*3 windows to calculate flow direction and based on the “overflow tracking” algorithm to determine the sink network, the algorithm process is as follows:

(2.1) Definition of grid point flow direction

The flow direction of each grid point in the maximum slope direction is calculated by searching in 8 directions with 3*3 window, different codes are assigned in 8 directions, as shown in the right figure, each grid has a value from 1 to 9 representing its direction to adjacent pixels, if the pixel is a concave point, the value is 5 (Fig. 9-15).

Handling of several exceptions:

1. If the maximum gradient grid point of a grid point has the same elevation value as that of the grid point, and no other grid point has flowed to the adjacent grid before, it is forced to flow to it. If another grid point flows to the adjacent grid, the current grid point is a concave point.

2. When the maximum gradient of two or more adjacent grid points is equal, first compare the gradient of each adjacent grid point, if it is still unresolved, continue to compare the relative gradient of grid points and decide to assign a flow direction.

3. For regions with the same elevation value, the search window radius is enlarged with 7*7 window, also can use a larger window if you want.

4. Add a circle of grid points with an elevation value of 0 on the periphery of the DEM data, forcing its maximum slope to flow outside the study area.

When all grid points are processed, generated a flow chart encoding 1-9.

(2.2) Concave processing algorithm

Due to the existence of the pits, some of the flow paths do not flow to the outlet of the watershed, but end to the pits, so the pits are processed before the automatic segmentation of the watershed. The pits in the basin can be either real pits or due to interpolation errors, so you can’t use a simple filtering or smoothing function to remove all the pits, the purpose is to connect the breaks caused by the pits to the main valley network.

Search for the lowest neighbors of all pits (sometimes there may be multiple elevations with equal elevations), as the overflow point of the pit, continue searching for neighbors that are lower or equal to its elevation from the point of overflow (already searched) Do not ignore), determine whether there is a grid point lower than the original pit, if not, repeat the above search process with the pit point of the pit as the starting point; if a grid point lower than the original pit is found, the direction of the pit and the lowest neighbor are reversed. As shown in Figure 9-15: the point with the elevation of 48 is a pit, and the lowest neighbor with the highest elevation is 49, continue searching with it as the starting point, find the elevation point 49, still higher than the original pit height, continue search, find another elevation point 49, and then find the elevation 47 point, lower than the original pit elevation, end the search, modify the flow direction according to the search direction, as shown by the direction of the solid arrow in Figure 9-16.

(2.3) Extraction of confluence network

According to the modified flow chart, given a point, the sum of all grid points flowing to it is the confluence area of the point. The calculation method is to search eight neighboring points for a given point, record all the positions of the grid points flowing to it, and then continue to search for the grid points flowing to it based on the found grid points until there is no new confluence point, all recorded grid points constitute the confluence area of the point.

Generally, the area of the confluence area of the valley is larger than the area of the confluence area of other grid points, and a grid value of the area of the confluence area larger than the threshold can be identified as a valley point by setting a threshold. Obviously, the complexity of the valley network obtained by different thresholds is different, although this method provides flexibility for determining the complexity of the valley network, it also makes the determination of the valley network too large and random.

After obtaining the valley network, the valley network can be coded. Firstly, the valley node is coded. The whole confluence network is searched and traversed from the outlet of the river basin, and the upstream and downstream nodes of each valley segment are coded and identified, the coded values of the valley segment are the coded values of the valley segment, and the locations of these nodes are recorded. Secondly, each grid point in the valley segment is identified as the coding value of the valley segment, thirdly, according to the type of the upstream node of the valley section, the valley section is determined as internal or external Valley section.

(2.4) Extraction of water distribution network

The confluence area of each grid point in the valley section is searched recursively, and the grid points in the confluence area are assigned to the identification value of the valley section to form the sub-confluence area of each valley section. Then, boundary tracking is carried out to extract the boundary line of sub-confluence area as watershed, and the watershed network is obtained. Finally, the valley network, watershed network and sub-confluence area are topologically coded to achieve automatic watershed terrain segmentation.

DEM calculates terrain properties

Terrain attribute data derived from DEM can be divided into single element attribute and compound attribute, the former can be directly calculated from the elevation data, such as the slope factor and the aspect. The latter is a composite index composed of several single-element attributes combined in a certain relationship to describe the spatial variation of a process, this combination relationship is usually an empirical relationship, and a simplified natural process mechanism model can also be used.

Single-factor terrain attributes can be easily computed by computer programs, including:

1) Slope, aspect

Slope is defined as the tangent between the horizontal plane and the local surface. It consists of two components: the slope - the maximum ratio of height changes (often called slope); slope - the direction of the maximum value of the change ratio. Geomorphology analysis may also use second-order differential pits and convexities. A more general metric is: the slope is measured as a percentage, the slope is measured at an angle from the north, and the convexity is measured as the degree of inclination within the unit distance.

The calculation of slope and aspect usually uses 3*3 window, after the window moves continuously in the DEM elevation matrix, the calculation of the whole image is completed, the slope is calculated as follows:

Slope direction is calculated as follows:

In order to improve the calculation speed and accuracy, GIS usually uses the second-order difference to calculate the slope and aspect, the simplest finite second-order difference method is to calculate the slope of the point i, j in the x direction:

In the formula, the grid spacing (along the diagonal line should be multiplied by). This method calculates the slope of eight directions, and the calculation speed is much faster. However, the local errors of ground elevation will cause serious errors in slope calculation. The better results can be obtained with digital analysis method, the formula for calculating slope in east-west direction with digital analysis method, which is as follows:

Similarly, formulas for calculating slope in other directions can be written out.

2)Area, volume

(2.1) Section area

According to the engineering design route, the intersection point P_i (X_i, Y_i, Z_i) with DEM grid edges can be calculated, then the line section area is

Where n is the number of intersections; D_i,i+1is the distance between P_iand P _i+1. The same can be used to calculate any cross section and its area.

(2.2) Volume

The DEM volume is obtained by accumulating the quadrangular prism (uncharacteristic grid) and the triangular prism volume. The upper surface of the quadrilateral cylinder is fitted with a parabolic hyperboloid, the upper surface of the triangular prism is fitted with an oblique plane, and the lower surface is a horizontal plane or reference plane, the calculation formulas are as follows:

Where S₃and S₄are the bottom areas of the triangular prism and the quadrangular prism, respectively.

According to the two DEMs, the amount of excavation, filling and soil loss can be calculated.

3) Surface area

For a grid containing features, it is decomposed into triangles, for a featureless grid, the elevation of the four corner points can be averaged, that is, the center point elevation, and then the grid is divided into four triangles. The slope passing through the three vertices is calculated from the coordinates of the three corner points of each triangle (x_i, y_i, z_i) The area of the inner triangle is finally accumulated to obtain the surface area of the ground.

Triangular network interpolation

After the TIN is created, the elevation of any point in the area can be solved by the TIN, the interpolation of TIN has different characteristics from the interpolation of rectangular grid, and the retrieval of points for interpolation is more complicated than the retrieval of grid. In general, only linear interpolation, that is, the inclined plane determined by the three points of the triangle is used as the ground surface, so that only the ground is continuous and the smoothness cannot be ensured. To perform triangulation interpolation, generally go through the following steps:

1) Searching for grid points

Given a plane coordinate P(x, y) of one point, based on the elevation Z of the point at which the TIN is interpolated, it is first determined which triangle of the TIN the point P falls in. The general approach is to calculate the distance and get the closest point according to point P, set to Q₁. Then you need to determine the triangle where P is located. Take the triangle with Q₁as the vertex in turn, and judge whether P is inside the triangle. It can be determined whether P is equal to each side of the triangle on the same side of the opposite side of the vertex (the coordinates of the point are the same as the value sign obtained by the straight line equation of the side).If P is not in any triangle with Qsub:`1`vertex, then the nearest grid point from P is taken and the above processing is repeated until the triangle where P is located is retrieved, that is, three grid points for interpolating the elevation of P point are retrieved.

2) Elevation interpolation

If the triangle of P (x, y) is_Q_1 Q_2 Q_3, the coordinates of three vertices are (x_1, y_1, z_1), (x_2, y_2, z_2) and (x_3, y_3, z_3), the plane equation determined by Q_1, Q_2 and Q_3 is

Or

Order

Then the elevation of point P is

Contour tracing

Based on TIN, the contours are drawn directly from the original observation data, avoiding the loss of accuracy of DTM interpolation, so the accuracy of the contours is higher; the shorter closed contours near the elevation annotation points can also be drawn; the contours drawn are distributed in the sampling area without requiring regular quadrilateral boundaries in the sampling area. The contour of the same elevation passes through a triangle at most once, so the program design is simpler. However, because of the different storage structure of TIN, the tracking algorithms of contour lines are different.

The contour drawing algorithm based on triangle search is as follows:

For TIN that records triangle tables, it searches in triangle order. Its basic process is as follows:

1. For a given contour elevation h, z_i (i=1,2,… For comparison, if z_i = h, z_i is added (or subtracted) to a small positive number e > 0 (e.gε=10_-4’), so that the programming is simple without affecting the accuracy of contours.

2. Set up an array of triangle markers, whose initial value is zero, and each element corresponds to a triangle. If the triangle is processed, the markers will be set to 1, and will not be processed until the contour elevation changes.

3. Determine in sequence whether two sides of each triangle have contours passing through. If the two ends of a triangle are P_1 (x_1, y_1, z_1), P_2 (x_2, y_2, z_2), then (z_1-h) (z:_2-h) < 0 indicates that there is a contour point on the side; (z_1-h) (z_2-h) > 0 indicates that there is no contour point on the side.

Until the first intersection point between the contour line and the mesh edge is found, the point is called the starting point of the search. It is also the plane coordinate (x, y) of the current triangle’s contour entering the edge and interpolating the point linearly.

1. Search for the contour line on the off side of the triangle, that is, the incoming side of the adjacent triangle, and interpolate its plane coordinates. The search and interpolation methods are the same as the above search starting points, but only the other two sides of the triangle are processed.

2. Enter the adjacent triangle and repeat step (4) until there is no adjacent triangle on the leaving edge (the contour is an open curve at this time) or the adjacent triangle is the triangle where the search starting point is located (the contour is a closed curve at this time).

3. For open curves, the searched contour points are reversed and searched in the other direction until the boundary is reached (that is, there is no adjacent triangle on the left side).

4. When all contours are tracked, they are output smoothly. The method is the same as that of drawing rectangular grid contours. Then continue to search for triangles until all triangles are processed, then change the contour elevation, repeat the above process until all contours are drawn. The process of generating contours with a value of 50 using triangular mesh is described in the figure.