Digital terrain surface area downscaling calculation method

Abstract

The invention belongs to the technical field of surveying and mapping geographic information data processing, and particularly relates to a digital terrain surface area downscaling calculation method. The invention aims to solve the problems of strong dependency of regional surface area calculation work on a spatial data structure of a digital terrain, non-rigorous mathematical basis and high calculation complexity. On the basis of a geodetic coordinate system, a downscaling discrete grid is established in a sub-spatial scale and serves as a spatial basis for expressing topographic surface morphological characteristics; in the downscaling discrete grid, the subspace scale is used as the size length of a grid unit, a triangular face is constructed in the grid unit, and the area of the triangular face is calculated to serve as the area of the downscaling grid unit.

Description

Technical Field

[0001] The present invention belongs to the technical field of surveying and mapping geographic information data processing, and in particular relates to a method for downscaling calculation of digital terrain surface area. Background Art

[0002] The earth's surface is rich and varied, and the terrain of different geomorphic units is significantly different. Obtaining area values ​​on the complex and diverse earth's surface is a basic and important task. The ups and downs of the terrain surface determine that the area value of the terrain surface area is a dependent variable of the ground elevation value. In the field of surveying and mapping geographic information, the plane area or ellipsoid area is significantly different from the surface area. Therefore, strictly speaking, the plane area value or ellipsoid area value cannot replace the surface area value, and can only be approximated to a certain extent in limited application scenarios. The confusing use of area values, especially the use of plane or ellipsoid area values ​​to replace the surface area value in strictly quantitative scenarios, will lead to systematic errors and serious impacts on applications and research in geography, ecology, and natural resource management, such as surface cover area estimation / land use change analysis, carbon sink measurement, and ecosystem service value estimation.

[0003] The current surface area calculation methods can be categorized into three main types in principle: slope-based calculation, surface roughness-based calculation, terrain unit differential fitting calculation, terrain unit approximation calculation, and calculation based on the rectangular coordinates of the earth-fixed geocentric space. Calculating the surface area based on slope is based on the plane area, and is obtained by multiplying the plane area by the inverse of the cosine of the terrain slope value (see the literature: Liu Xingyu, Zhou Guangsheng, Lv Xiaomin, et al. Analysis of the difference between the actual surface area and the vertical projection area in the Hengduan Mountain area - a case study of Yajiang County [J]. Bulletin of Surveying and Mapping, 2021(8): 37-41.). In essence, this type of method cannot take into account the difference in area change between the terrain unit before and after projection. The error is not obvious in a small area, but the error accumulates to a large extent in a large area. The calculation method based on surface roughness starts from the concept of surface roughness and multiplies the plane area by the surface roughness conversion coefficient to obtain the area value of the rough terrain surface (see the literature: Zeng Zhen, Yang Benyong, Fan Jianrong, et al. Calculation of the true surface area based on the geoscientific significance of surface roughness [J]. Remote Sensing Technology and Applications, 2014, 29(5):846-852.). This type of method inherits the shortcomings of the calculation method based on slope, and also introduces errors in the calculation of surface roughness. The terrain unit differential fitting method is mostly to use mathematical functions to divide the local terrain undulation into micro-elements in a local area, replace the terrain surface micro-elements with geometric micro-elements, and obtain the area of ​​the local terrain by accumulating the areas of the geometric micro-elements (see the literature: Chen Jilong, Wu Wei, Liu Hongbin. Research on the application of DEM in forest surface area calculation [J]. Journal of Southwest Agriculture, 2008, (05): 1348-1352. Gong Jinqi, Chen Yonggang, Lu Jianqing, et al. Accurate calculation of forest surface area using vector extraction and dynamic segmentation [J]. Journal of Wuhan University·Information Science Edition, 2021.). This type of method has high surface area calculation accuracy, but its accuracy strictly depends on the accuracy of the terrain data. At the same time, the calculation complexity is high and the amount of calculation is large, which makes it not suitable for large-scale implementation. The terrain unit approximation calculation method uses terrain factors to encrypt the DEM and polygonal area boundaries, and then constructs a surface triangulation network to calculate the surface area (see the literature: Xue Shuqiang, Dang Yamin, Mi Jinzhong, et al. Surface area calculation considering nonlinear terrain factors [J]. Journal of Geodesy and Cartography, 2015, 44(03): 330-337. Han Maixia, Cheng Chuanlu, Wang Wei. Research on surface area calculation methods [J]. Bulletin of Surveying and Mapping, 2015, (09): 72-74+82.). This type of method requires first using the Taylor series approximation principle to perform least squares estimation on the micro-topographic factors, and then using them as input for surface area calculation. The steps are complicated and the calculation complexity is high. The calculation method based on the fixed geocentric space rectangular coordinate system converts the surface area into a spatial rectangular coordinate system, and then uses the terrain unit differential fitting to calculate the surface area.This type of calculation takes into account the actual undulating shape of the terrain unit, but is still limited by the complex and changeable geometric shapes and number of boundary points of the surface unit, resulting in the area calculation process being highly dependent on the distribution of surface unit boundary points and geometric surface division. The overall efficiency is low and the applicability is not strong (see the literature: Liu Jiping, Dong Chun, Kang Xiaochen, et al. Statistical Analysis of National Geographic Conditions in the Big Data Era [J]. Journal of Wuhan University (Information Science Edition), 2019, 44(01): 68-76+83.). Summary of the Invention

[0004] To address the problems of regional-level surface area calculations being highly dependent on the spatial data structure of digital terrain, having a loose mathematical foundation, and exhibiting high computational complexity, the present invention provides a method for downscaling digital terrain surface area. This method, based on a geodetic coordinate system, establishes a downscaled discrete grid at a subspatial scale as the spatial basis for expressing the surface morphological characteristics of the terrain. Within the downscaled discrete grid, the subspatial scale is used as the grid cell length, triangular faces are constructed within the grid cells, and the area of ​​the triangular faces is calculated as the area of ​​the downscaled grid cells.

[0005] In order to achieve the above-mentioned purpose, the present invention adopts the following technical solutions:

[0006] The present invention provides a method for downscaling surface area of ​​digital terrain, comprising the following steps:

[0007] Step 1: Based on the geodetic coordinate system, within the digital terrain range, refer to the spatial scale of the minimum digital terrain unit, calculate the sub-spatial scale R that is better than its spatial resolution, and use it to establish a downscaled discrete grid;

[0008] Step 2: Calculate the geodetic longitude and latitude coordinates of the grid points based on the established downscaled discrete grid;

[0009] Step 3: Perform elevation interpolation calculation on the grid points according to the elevation distribution within the neighborhood, and then form the geodetic coordinates of the grid points;

[0010] Step 4: converting the geodetic coordinates of the grid points into geocentric, earth-fixed space rectangular coordinates consistent with the geodetic coordinate system datum;

[0011] Step 5: In the downscaled discrete grid, the grid cell length is determined by the subspace scale, a triangular face is constructed inside the grid cell, and the area of ​​the triangular face is calculated;

[0012] Step 6: Accumulate the areas of all triangles in the grid unit as the downscaled grid unit area. The downscaled area of ​​the digital terrain is represented by the area of ​​each downscaled grid unit of all downscaled discrete grids.

[0013] Furthermore, in step 1, the sub-spatial scale R is obtained by multiplying the spatial resolution of the digital terrain by the percentage transformation ratio;

[0014] The downscaled discrete grid is a square global discrete grid with the sub-spatial scale R as the minimum grid cell size length.

[0015] Furthermore, the geodetic coordinate system in step 1 refers to the geocentric geodetic coordinate system in geodesy, which uses three components, geodetic longitude L, geodetic latitude B, and geodetic height H, to represent the position (B, L, H).

[0016] Digital terrain is a method of recording ground elevation values ​​in a regular or irregular form to represent the surface morphology, and records the geodetic elevation values ​​H of surface points at a certain spatial resolution.

[0017] The spatial scale of the smallest digital terrain unit is defined by the spatial resolution of the digital terrain.

[0018] Furthermore, the geodetic longitude and latitude coordinates of the grid points are calculated in step 2, specifically:

[0019] The grid points are the four corner points of each grid cell in the square global discrete grid. According to the row / column number (m, n) of the grid point in the downscaled discrete grid, the geographic longitude and latitude (B, L) coordinates are calculated by affine transformation by recording the geographic affine transformation parameters of the downscaled discrete grid range, as shown in formula (1):

[0020] (1)

[0021] Among them, a0 is the longitude coordinate value of the upper left starting point, a1 is the pixel block size in the east-west direction, a2 is the column rotation parameter, a3 is the latitude coordinate value of the upper left starting point, a4 is the row rotation parameter, and a5 is the pixel block size in the north-south direction.

[0022] Furthermore, the neighborhood range in step 3 is a circular area centered on the grid point and with a radius greater than or equal to the subspace scale R;

[0023] Elevation interpolation calculation uses the average of the digital terrain elevation values ​​within the neighborhood as the elevation value of the grid point.

[0024] Furthermore, the geodetic coordinates of the grid points in step 4 are the positions of the grid points in the geocentric geodetic coordinate system represented by coordinate values ​​(B, L, H);

[0025] The Earth-centered Earth-fixed space rectangular coordinates are expressed as coordinate values ​​(X, Y, Z). They are three-dimensional space rectangular coordinate values ​​with the same geodetic datum as the geodetic coordinate system. They are obtained by converting geodetic coordinates (B, L, H) into space rectangular coordinates (X, Y, Z), as shown in formula (2):

[0026] (2)

[0027] Where X, Y, and Z are the three components of the spatial rectangular coordinate system, N is the radius of curvature of the y-axis, and e is the first eccentricity. The calculations of N and e are as follows:

[0028] (3)

[0029] (4)

[0030] Among them, a is the major semi-axis of the ellipsoid of the geodetic coordinate system reference, and b is the minor semi-axis of the ellipsoid of the geodetic coordinate system reference.

[0031] Furthermore, the triangular face constructed inside the grid cell in step 5 is formed by connecting the four corner points of the square grid cell, connecting any diagonal of the square, and connecting the four sides of the square grid cell through one diagonal line.

[0032] The area of ​​a triangle is calculated from the spatial rectangular coordinates of the three corner points of the triangle.

[0033] Furthermore, in step 6, the accumulation of the areas of all triangular faces in the grid cell is the accumulation of the areas of the two triangular faces inside each grid cell (as the downscaled area of ​​the digital terrain within the corresponding grid cell), and this processing is performed on all grid cells of the downscaled discrete grid.

[0034] Compared with the prior art, the present invention has the following advantages:

[0035] 1) The vector data structure or raster data structure commonly used in the prior art is used as the spatial basis for expressing the morphological characteristics of the terrain surface. The present invention is based on the geodetic coordinate system and establishes a downscaled discrete grid at a sub-spatial scale as the spatial basis for expressing the morphological characteristics of the terrain surface. The prior art usually directly uses the polygon structure in the vector data structure or the pixel structure in the raster data structure to calculate the surface area. The present invention uses the sub-spatial scale as the grid unit size length within the downscaled discrete grid, constructs a triangular face inside the grid unit, and calculates the area of ​​the triangular face as the downscaled grid unit area.

[0036] 2) The present invention is not limited to vector data structures or raster data structures and has strong adaptability to basic spatial representations of terrain. The present invention has a strict mathematical foundation and clear algorithmic steps, which are beneficial to ensuring calculation accuracy and computer programming implementation. The present invention calculates surface area at sub-spatial scale grid cells, and can improve the spatial resolution of surface area cells and the expression accuracy of surface area through downscaling processing based on the original elevation data. The present invention uses a downscaled discrete grid as a spatial basis, which facilitates the parallelization and high-performance implementation of large-scale surface area calculations. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 Schematic diagram of the downscaled discrete grid, where a is the spatial resolution of the digital terrain data, b is the subspatial scale, and c is the downscaled discrete grid established with b as the subspatial scale.

[0038] Figure 2 The figure shows the spatial interpolation calculation of grid points on regularly distributed elevation data. Grid point A directly uses the value of elevation point (4, 2); grid point B uses the average of elevation points (4, 2), (4, 3), (3, 2), and (3, 3); grid point C uses the average of elevation points (3, 3) and (2, 3); and grid point D uses the average of elevation points (3, 3) and (3, 4).

[0039] Figure 3 The figure is a schematic diagram of spatial interpolation calculation of grid points on irregularly distributed elevation data. Grid point A takes the average of the elevation points within a circular neighborhood with A as the center and the subspatial scale as the radius, and grid point B takes the average of the elevation points within a circular neighborhood with a radius larger than the subspatial scale.

[0040] Figure 4 A schematic diagram of the internal triangulated surface constructed from the four corner points of a grid cell. DETAILED DESCRIPTION

[0041] In order to further illustrate the technical solution of the present invention, the present invention is further described below through examples. Example 1

[0042] A method for downscaling digital terrain surface area in this embodiment includes the following steps:

[0043] Step 1: Based on the geodetic coordinate system within the digital terrain range, refer to the spatial scale of the minimum digital terrain unit, calculate the sub-spatial scale R that is better than its spatial resolution, and use it to establish a downscaled discrete grid ( Figure 1 );

[0044] The subspatial scale R is obtained by multiplying the spatial resolution of the digital terrain and the percentage transformation ratio; the downscaled discrete grid is a square global discrete grid with the subspatial scale R as the minimum grid unit size length.

[0045] The geodetic coordinate system refers to the geocentric geodetic coordinate system in geodesy, in which the position is represented by three components: geodetic longitude L, geodetic latitude B, and geodetic height H (B, L, H). Digital terrain records ground elevation values ​​in a regular or irregular form to characterize the surface morphology, and records the geodetic height H of surface points at a certain spatial resolution. The spatial scale of the smallest digital terrain unit is defined by the spatial resolution of the digital terrain.

[0046] Step 2: Calculate the geodetic longitude and latitude coordinates of the grid points based on the established downscaled discrete grid;

[0047] The grid points are the four corner points of each grid cell in the square global discrete grid. According to the row / column number (m, n) of the grid point in the downscaled discrete grid, the geographic longitude and latitude (B, L) coordinates are calculated by affine transformation by recording the geographic affine transformation parameters of the downscaled discrete grid range, as shown in formula (1):

[0048] (1)

[0049] Among them, a0 is the longitude coordinate value of the upper left starting point, a1 is the pixel block size in the east-west direction, a2 is the column rotation parameter (usually 0), a3 is the latitude coordinate value of the upper left starting point, a4 is the row rotation parameter (usually 0), and a5 is the pixel block size in the north-south direction.

[0050] Step 3: Perform elevation interpolation calculation on the grid points according to the elevation distribution within the neighborhood ( Figure 2 and Figure 3 ), and then form the geodetic coordinates of the grid points;

[0051] The domain range is a circular area with the grid point as the center and a radius greater than or equal to the sub-spatial scale R; the elevation interpolation calculation uses the average of the digital terrain elevation values ​​within the neighborhood range as the elevation value of the grid point.

[0052] Step 4: converting the geodetic coordinates of the grid points into geocentric, earth-fixed space rectangular coordinates consistent with the geodetic coordinate system datum;

[0053] The geodetic coordinates of the grid points are the grid point positions in the geocentric geodetic coordinate system represented by the coordinate values ​​(B, L, H); the geocentric earth-fixed spatial rectangular coordinates are the three-dimensional spatial rectangular coordinate values ​​with the same geodetic datum as the geodetic coordinate system represented by the coordinate values ​​(X, Y, Z), which are obtained by converting the geodetic coordinates (B, L, H) into spatial rectangular coordinates (X, Y, Z), as shown in formula (2):

[0054] (2)

[0055] Where X, Y, and Z are the three components of the spatial rectangular coordinate system, N is the radius of curvature of the y-axis, and e is the first eccentricity. The calculations of N and e are as follows:

[0056] (3)

[0057] (4)

[0058] Among them, a is the major semi-axis of the ellipsoid of the geodetic coordinate system reference, and b is the minor semi-axis of the ellipsoid of the geodetic coordinate system reference.

[0059] Step 5: In the downscaled discrete grid, the grid cell length is determined by the subspace scale, a triangular face is constructed inside the grid cell, and the area of ​​the triangular face is calculated ( Figure 4 );

[0060] The triangular face constructed inside the grid cell is composed of two triangular faces connected by the four corner points of the square grid cell, connecting any diagonal of the square, and connecting through one diagonal and the four sides of the square grid cell. The area of ​​the triangular face is calculated by the spatial rectangular coordinate values ​​of the three corner points of the triangle.

[0061] Step 6: Accumulate the areas of all triangles in the grid cell as the downscaled grid cell area. The downscaled area of ​​the digital terrain is represented by the area of ​​each downscaled grid cell of all downscaled discrete grids.

[0062] The accumulation of the areas of all triangular faces in a grid cell is to accumulate the areas of the two triangular faces inside each grid cell, and this process is performed on all grid cells of the downscaled discrete grid.

[0063] More specifically, the digital terrain data used in this embodiment is digital elevation model data (ALOSDEM), which has a spatial resolution of 12.5 meters. Using this elevation data, a digital terrain surface area downscaling calculation is performed within the surface range represented by this data, including the following steps:

[0064] Step 1: Establish the WGS84 geodetic coordinate system. Based on the minimum spatial scale of the digital terrain unit, that is, the spatial resolution of the ALOS DEM data is 12.5 meters, determine the subspatial scale to be 6.25 meters, and establish a square downscaled discrete grid within the surface area with a grid size of 6.25 meters.

[0065] Step 2: Based on the row / column number (m, n) of the grid point in the downscaled grid, the geographic affine transformation parameters of the discretized grid range are recorded and the geodetic latitude and longitude (B, L) coordinate values ​​are obtained through affine transformation calculation;

[0066] Step 3: Perform elevation interpolation calculations on each grid point of the downscaled grid based on the elevation distribution within the neighborhood. In this embodiment, the surface elevation values ​​in the digital elevation model data are regularly distributed. Therefore, the elevation spatial interpolation calculation of the grid points is performed based on the distribution of the grid points and the elevation data points. The grid point elevation value H is obtained and the geodetic coordinates (B, L, H) of the grid point are formed.

[0067] Step 4: According to the coordinate conversion formula, traverse all grid points by row / number, and convert the geodetic coordinate values ​​(B, L, H) of the grid points into the Earth-centered Earth-fixed space rectangular coordinate values ​​(X, Y, Z) corresponding to the WGS84 ellipsoid;

[0068] Step 5: In the downscaled discrete grid, traverse all grid cells. Inside each grid cell, connect the two diagonal points of the grid cell to divide the grid cell into two triangles. Use the spatial rectangular coordinates (X, Y, Z) of the grid points and the formula for calculating the area of ​​a triangle using coordinates in spatial analytic geometry to calculate the area of ​​the triangle face S. t ;

[0069] Step 6: Accumulate the area S of the two triangles in the grid cell t1 and S t2 , as the downscaled grid unit area S, the surface area covering the surface range of the elevation data with a spatial resolution of 6.25m is obtained.

[0070] The foregoing shows and describes the principal features and advantages of the present invention. It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above and that the present invention can be embodied in other specific forms without departing from the spirit or essential characteristics of the present invention. Therefore, the embodiments should be considered in all respects as illustrative and non-restrictive, and the scope of the present invention is defined by the appended claims, not the foregoing description, and all variations that come within the meaning and range of equivalents of the claims are intended to be embraced therein.

[0071] In addition, it should be understood that although this specification is described in terms of implementation methods, not every implementation method contains only one independent technical solution. This narrative method of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment can also be appropriately combined to form other implementation methods that can be understood by those skilled in the art.

Claims

1. A method for downscaling surface area of ​​digital terrain, characterized in that: The following steps are involved: Step 1: Based on the geodetic coordinate system, within the digital terrain range, refer to the spatial scale of the minimum digital terrain unit, calculate the sub-spatial scale R that is better than its spatial resolution, and use it to establish a downscaled discrete grid; Step 2: Calculate the geodetic longitude and latitude coordinates of the grid points based on the established downscaled discrete grid; Step 3: Perform elevation interpolation calculation on the grid points according to the elevation distribution within the neighborhood, and then form the geodetic coordinates of the grid points; Step 4: converting the geodetic coordinates of the grid points into geocentric, earth-fixed space rectangular coordinates consistent with the geodetic coordinate system datum; Step 5: In the downscaled discrete grid, the grid cell length is determined by the subspace scale, a triangular face is constructed inside the grid cell, and the area of ​​the triangular face is calculated; Step 6: Accumulate the areas of all triangles in the grid unit as the downscaled grid unit area. The downscaled area of ​​the digital terrain is represented by the area of ​​each downscaled grid unit of all downscaled discrete grids.

2. The method for calculating surface area downscaling of digital terrain according to claim 1, characterized in that: In step 1, the sub-spatial scale R is obtained by multiplying the spatial resolution of the digital terrain by the percentage transformation ratio; The downscaled discrete grid is a square global discrete grid with the sub-spatial scale R as the minimum grid cell size length.

3. The method for calculating surface area downscaling of digital terrain according to claim 2, characterized in that: The geodetic coordinate system in step 1 refers to the geocentric geodetic coordinate system in geodesy, which uses three components, geodetic longitude L, geodetic latitude B, and geodetic height H, to represent the position (B, L, H). Digital terrain is a method of recording ground elevation values ​​in a regular or irregular form to represent the surface morphology, and records the geodetic elevation values ​​H of surface points at a certain spatial resolution. The spatial scale of the smallest digital terrain unit is defined by the spatial resolution of the digital terrain.

4. The method for calculating surface area downscaling of digital terrain according to claim 2, characterized in that: The geodetic longitude and latitude coordinates of the grid points calculated in step 2 are specifically: The grid points are the four corner points of each grid cell in the square global discrete grid. According to the row / column number (m, n) of the grid point in the downscaled discrete grid, the geographic longitude and latitude (B, L) coordinates are calculated by affine transformation by recording the geographic affine transformation parameters of the downscaled discrete grid range, as shown in formula (1): （1） Among them, a0 is the longitude coordinate value of the upper left starting point, a1 is the pixel block size in the east-west direction, a2 is the column rotation parameter, a3 is the latitude coordinate value of the upper left starting point, a4 is the row rotation parameter, and a5 is the pixel block size in the north-south direction.

5. The method for calculating surface area downscaling of digital terrain according to claim 1, characterized in that: The neighborhood range in step 3 is a circular area centered on the grid point and with a radius greater than or equal to the subspace scale R; Elevation interpolation calculation uses the average of the digital terrain elevation values ​​within the neighborhood as the elevation value of the grid point.

6. The method for calculating surface area downscaling of digital terrain according to claim 1, characterized in that: The geodetic coordinates of the grid points in step 4 are the grid point positions in the geocentric geodetic coordinate system represented by coordinate values ​​(B, L, H); The Earth-centered Earth-fixed space rectangular coordinates are expressed as coordinate values ​​(X, Y, Z). They are three-dimensional space rectangular coordinate values ​​with the same geodetic datum as the geodetic coordinate system. They are obtained by converting geodetic coordinates (B, L, H) into space rectangular coordinates (X, Y, Z), as shown in formula (2): （2） Where X, Y, and Z are the three components of the spatial rectangular coordinate system, N is the radius of curvature of the y-axis, and e is the first eccentricity. The calculations of N and e are as follows: （3） （4） Among them, a is the major semi-axis of the ellipsoid of the geodetic coordinate system reference, and b is the minor semi-axis of the ellipsoid of the geodetic coordinate system reference.

7. The method for calculating surface area downscaling of digital terrain according to claim 1, characterized in that: The triangular face constructed inside the grid cell in step 5 is formed by connecting the four corner points of the square grid cell, connecting any diagonal of the square, and connecting the four sides of the square grid cell through one diagonal. The area of ​​a triangle is calculated from the spatial rectangular coordinates of the three corner points of the triangle.

8. The method for calculating surface area downscaling of digital terrain according to claim 1, characterized in that: In step 6, the accumulation of the areas of all triangular faces in the grid cells is performed by accumulating the areas of the two triangular faces inside each grid cell. This process is performed on all grid cells of the downscaled discrete grid.

Images