DEM (Digital Elevation Model) rapid updating method based on minimum elevation difference mosaic

Abstract

The invention discloses a DEM rapid updating method and system based on minimum elevation difference mosaic, and belongs to the field of DEM dynamic updating, and the method comprises the steps: obtaining and processing natural resource monitoring data and regional historical DEM data; dEM change areas are extracted based on land type change rules, and pattern spots meeting area requirements are screened to serve as to-be-updated areas; generating a point cloud through dense matching of aerial triangulation achievements of the aerial image, and making a new DEM of a change area through processing; plane matching of the new DEM and the historical DEM is achieved through a SiftGPU, Hough transform and a UMEYAMA algorithm; setting a buffer area for a change range, calculating an elevation difference weight, generating an optimal mosaic line by using a minimum cost path algorithm, and implementing two-stage DEM seamless mosaic through feathering; and verifying the updated DEM precision based on field actual measurement points. The method and the system do not need a large amount of manual editing, and can improve the DEM updating efficiency.

Description

Technical Field

[0001] This invention relates to the field of dynamic DEM updating, and specifically to a method for rapid DEM updating based on minimum elevation difference mosaic. Background Technology

[0002] A Digital Elevation Model (DEM) simulates a continuous terrain surface using a finite number of data points with elevation attributes. It can recreate terrain features to a certain extent, and the degree of terrain reconstruction is a key indicator of DEM quality. Because DEMs describe ground elevation information, they are widely used in surveying, hydrology, meteorology, geomorphology, geology, soil science, engineering construction, communications, and military fields. These fields encompass national economic and defense construction, as well as humanities and natural sciences. Furthermore, DEMs have become one of the most important spatial information data in Geographic Information System (GIS) databases.

[0003] As one of the national basic geographic information data products, the production and updating cycle of DEMs is relatively long, and there are three main methods for their production and updating. The first method is production using topographic maps. This method is fast, low-cost, and relatively simple in process. However, due to limitations in the timeliness of topographic maps and the capacity of map elements, the accuracy of the final results is difficult to meet the accuracy requirements of large-scale datasets. The second method is production using ground surveying. This method yields high-accuracy results, but its production efficiency is low, and most of the work needs to be carried out in the field, resulting in high production costs. Therefore, it is not suitable for large-scale production and updating work. The third method is remote sensing surveying, specifically including aerial photogrammetry, satellite stereo imagery surveying, laser scanning point cloud surveying, and SAR image interferometry. This method can quickly acquire data over a large area, making it suitable for large-scale data production or updating work. Among these types of remote sensing measurements, satellite stereo image measurement and SAR image interferometry have lower data acquisition costs, but the accuracy of the results obtained by these two methods is low due to limitations in image resolution and positioning accuracy. On the other hand, laser scanning point cloud measurement and aerial photogrammetry have higher data acquisition costs, and these two methods also have higher requirements for data processing technology and a larger workload for data processing, so they are also not suitable for large-scale production work.

[0004] With the continuous development of the economy and society and the accelerating pace of urbanization, urban infrastructure construction is also developing rapidly. Numerous construction activities have altered the surface morphology of some areas, causing significant changes in local topography. This directly leads to inconsistencies between DEM data and the actual terrain. Furthermore, the update frequency of DEMs struggles to keep pace with urban development, failing to meet the needs of practical applications. Summary of the Invention

[0005] The purpose of this invention is to provide a method for rapid DEM updating based on minimum elevation difference mosaic, so as to achieve rapid DEM updating.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] A fast DEM update method based on minimum elevation difference mosaicking includes the following steps:

[0008] Acquire natural resource monitoring data and existing regional DEM data;

[0009] DEM change area extraction rules are established based on the natural resource monitoring data. Land use patches in the natural resource monitoring data are analyzed according to the extraction rules. Adjacent land use patches of different types are merged or integrated. Patches that meet the area requirements are selected to obtain the DEM area to be updated.

[0010] Point cloud data is generated using aerial triangulation results from aerial imagery. Ground points are extracted from the point cloud data by filtering. The point cloud data is then expanded, cropped, and ground points are edited according to the DEM to be updated area to obtain the changed area DEM.

[0011] Feature matching and gross error removal are performed on the existing DEM of the region and the DEM of the changed region. The transformation parameters of the DEM of the changed region and the existing DEM of the region are calculated based on the processed matching points. The planar position transformation of the DEM of the changed region is performed according to the transformation parameters to obtain the DEM after planar matching.

[0012] A buffer is set for the DEM to be updated area. The elevation difference between the existing DEM and the plane-matched DEM in the buffer is calculated as a reference weight for subsequent processing. Based on the reference weight map, the minimum cost path of the elevation difference between the two DEMs is found. The minimum cost path is spliced ​​to form a closed line. The closed line is smoothed to obtain the final mosaic line.

[0013] Based on the final mosaic line, the existing DEM and the plane-matched DEM in the region are replaced according to the weight calculation rules to obtain the updated DEM.

[0014] In one possible implementation, after obtaining the updated DEM, the three-dimensional coordinates of the verification points are obtained through field measurements, the elevations of the corresponding verification points in the updated DEM are extracted, the elevation errors of the verification points are calculated, and the accuracy of the updated DEM is verified based on the elevation errors.

[0015] In one possible implementation, when acquiring the natural resource monitoring data, annual land change survey data, routine monitoring patch data, and urban land space monitoring data are acquired. Attribute analysis is then performed on the annual land change survey data, the routine monitoring patch data, and the urban land space monitoring data to obtain the natural resource monitoring data containing land use change information.

[0016] When acquiring existing DEM data for the region, digital elevation model data containing different resolutions is acquired, and the digital elevation model data is merged to obtain existing DEM data for the region.

[0017] In one possible implementation, when establishing the DEM change area extraction rules, the land use types that are prone to causing changes in terrain elevation are analyzed, the mapping relationship between land use patches and land use types in the natural resource monitoring data is sorted out, and the DEM change area extraction rules are established based on the mapping relationship.

[0018] In one possible implementation, after generating the point cloud data, if ground points cannot be obtained or the ground point density is insufficient in a local area, auxiliary points are added according to the region using existing or surrounding ground points and vegetation height to fit the ground of the local area; if the point cloud has edge confluence or time difference issues with ground features, a small-scale aerial triangulation region is constructed by selecting the latest photo, and the point cloud data is re-generated by dense matching of the small-scale aerial triangulation region. Then, the point cloud data is expanded, cropped, and the ground points are edited to obtain the DEM of the changed region.

[0019] In one possible implementation, when performing feature matching between the existing DEM of the region and the changed DEM of the region, feature points are extracted using the SiftGPU algorithm, the matching range is determined using the initial positioning information of the existing DEM of the region and the changed DEM of the region, the feature points are described within the matching range, and the initial matching point is determined.

[0020] In one possible implementation, when removing gross errors from the initial matching points, the Hough transform is used to filter the initial matching points based on the overall geometric consistency constraint of the stereo image, removing matching points that do not conform to the overall geometric consistency, thus obtaining matching points after removing gross errors. When calculating the transformation parameters, the two sets of object-side 3D coordinates corresponding to the matching points after removing gross errors are arranged into matrices using the UMEYAMA algorithm, and singular value decomposition is performed on the matrices. The scale, rotation matrix, and translation are solved based on the decomposition results to obtain the transformation parameters.

[0021] In one possible implementation, when setting the buffer, buffer parameters are determined, and the buffer is generated around the DEM area to be updated. When searching for the minimum cost path based on the elevation difference, the elevation difference is used as a reference weight to generate a cost matrix, and a gradient difference matrix of the same size as the cost matrix is ​​generated. The cost matrix and the gradient difference matrix are normalized. The two normalized matrices are traversed, and the larger value is taken to update the cost matrix. The minimum cost path is then found using the A* algorithm in combination with the updated cost matrix.

[0022] In one possible implementation, when performing data replacement according to the weight calculation rules, the elevation of non-grid points involved in the replacement is obtained using bilinear linear interpolation; the weight is determined based on the distance from any point to the final mosaic line; when the distance is greater than or equal to the feathering parameter expansion distance, the weight is set to 1; when the distance is less than the feathering parameter expansion distance and greater than or equal to the negative value of the contraction distance, the weight is determined by a weight calculation function; when the distance is less than the negative value of the contraction distance, the weight is set to 0; the final elevation of the arbitrary point is calculated based on the weight, and the data replacement is completed.

[0023] In one possible implementation, when performing the accuracy verification, flat land, hills, mountains, and high mountains are selected as accuracy verification areas. Verification points are selected within the accuracy verification areas for field measurements to obtain the three-dimensional coordinates of the verification points. After calculating the elevation error, it is determined whether the elevation error meets the preset mean square error requirement and does not exceed the maximum limit. If it meets the requirement, the updated DEM accuracy is determined to be qualified. If it does not meet the requirement, the aforementioned steps are repeated until the updated DEM accuracy is qualified.

[0024] Compared with the prior art, the advantages of this invention are as follows:

[0025] This invention utilizes land use change rules in natural resource monitoring data (including annual land change survey data, routine monitoring patch data, and urban land space monitoring data) to automatically and quickly identify areas of terrain change without requiring a comprehensive, manual screening of the entire region, significantly reducing the workload of locating terrain change areas. Furthermore, during the two-phase DEM plane matching process, SiftGPU technology leverages the parallel computing capabilities of the graphics processor to significantly improve the processing efficiency of feature point extraction and matching. The Hough transform can automatically eliminate gross errors in feature point extraction, and the UMEYAMA algorithm can quickly solve for the optimal coordinate system transformation parameters. The combination of these three technologies enables rapid spatial alignment of the two-phase DEMs, avoiding a large amount of manual adjustment operations in the traditional matching process, further shortening the overall DEM update cycle and improving production efficiency. Attached Figure Description

[0026] To more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0027] Figure 1 This is a technical roadmap of the DEM fast update method based on minimum elevation difference mosaic in the embodiments of the present invention;

[0028] Figure 2 This is a logic diagram for extracting the terrain change range in this embodiment of the invention;

[0029] Figure 3 This is a flowchart of the model matching technology in an embodiment of the present invention;

[0030] Figure 4 This is a schematic diagram illustrating the extraction of the DEM variation range of the study area in an embodiment of the present invention;

[0031] Figure 5 This is a schematic diagram of manual editing and modification of DEM in an embodiment of the present invention;

[0032] Figure 6 This is a schematic diagram of the automatic generation and updating of DEM mosaic lines based on the minimum cost path algorithm of elevation difference weighted map in an embodiment of the present invention.

[0033] Figure 7 This is a schematic diagram of the inlay line difference in an embodiment of the present invention. Detailed Implementation

[0034] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of this application.

[0035] Example:

[0036] It should be noted that the terms "comprising" and "having" and any variations thereof in the embodiments of the present invention are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such processes, methods, products, or devices.

[0037] This invention provides a method for rapid DEM updating based on minimum elevation difference mosaicking, comprising the following steps:

[0038] Step 101: Obtain natural resource monitoring data and existing regional DEM data.

[0039] In this step, when acquiring the natural resource monitoring data, annual land change survey data, routine monitoring patch data, and urban land space monitoring data are acquired. Attribute analysis is performed on the annual land change survey data, the routine monitoring patch data, and the urban land space monitoring data to obtain the natural resource monitoring data containing land use change information.

[0040] Furthermore, when acquiring existing DEM data for the region, digital elevation model data containing different resolutions is acquired, and the digital elevation model data is merged to obtain existing DEM data for the region.

[0041] Step 102: Establish DEM change area extraction rules based on the natural resource monitoring data, analyze land use patches in the natural resource monitoring data according to the extraction rules, merge or integrate adjacent land use patches of different types, and filter patches that meet the area requirements to obtain the DEM area to be updated.

[0042] In this step, when establishing the DEM change area extraction rules, the land use types that are prone to causing changes in terrain elevation are analyzed, the mapping relationship between land use patches and land use types in the natural resource monitoring data is sorted out, and the DEM change area extraction rules are established based on the mapping relationship.

[0043] Step 103: Generate point cloud data using aerial triangulation results from aerial imagery, filter the point cloud data to extract ground points, and perform outward cropping and ground point editing on the point cloud data according to the DEM to be updated area to obtain the changed area DEM.

[0044] In this step, after generating the point cloud data, if ground points cannot be obtained in a local area or the ground point density is insufficient, auxiliary points are added according to the region using existing or surrounding ground points and vegetation height to fit the ground of the local area; if the point cloud has edge confluence or time difference issues with ground features, a small-scale aerial triangulation region is constructed by selecting the latest photo, and the point cloud data is re-generated by dense matching of the small-scale aerial triangulation region. Then, the point cloud data is expanded, cropped, and the ground points are edited to obtain the DEM of the changed region.

[0045] Step 104: Perform feature matching and gross error removal on the existing DEM of the region and the DEM of the changed region. Calculate the transformation parameters of the DEM of the changed region and the existing DEM of the region based on the processed matching points. Perform planar position transformation on the DEM of the changed region according to the transformation parameters to obtain the DEM after planar matching.

[0046] In this step, when performing feature matching on the existing DEM of the region and the changed DEM of the region, feature points are extracted using the SiftGPU algorithm. The matching range is determined using the initial positioning information of the existing DEM of the region and the changed DEM of the region. The feature points are described within the matching range to determine the initial matching point.

[0047] Furthermore, when removing gross errors from the initial matching points, based on the overall geometric consistency constraint of the stereo image, the Hough transform is used to filter the initial matching points, removing those that do not conform to the overall geometric consistency, thus obtaining the matching points after removing gross errors. When calculating the transformation parameters, the two sets of object-side 3D coordinates corresponding to the matching points after removing gross errors are arranged into matrices using the UMEYAMA algorithm. Singular value decomposition is performed on the matrices, and the scale, rotation matrix, and translation are solved based on the decomposition results to obtain the transformation parameters.

[0048] Step 105: Set a buffer for the DEM to be updated area, calculate the elevation difference between the existing DEM and the DEM after plane matching in the buffer, find the minimum cost path based on the elevation difference, stitch the minimum cost path to form a closed line, smooth the closed line to obtain the final mosaic line.

[0049] In this step, when setting the buffer, the buffer parameters are determined, and the buffer is generated around the DEM area to be updated. When finding the minimum cost path based on the elevation difference, the elevation difference is used as a reference weight to generate a cost matrix, and a gradient difference matrix of the same size as the cost matrix is ​​generated. The cost matrix and the gradient difference matrix are normalized. The two normalized matrices are traversed, and the larger value is taken to update the cost matrix. The minimum cost path is found by combining the updated cost matrix with the A* algorithm.

[0050] Step 106: Based on the final mosaic line, replace the data of the existing DEM and the plane-matched DEM in the region according to the weight calculation rules to obtain the updated DEM.

[0051] In this step, when replacing data according to the weight calculation rules, bilinear interpolation is used to obtain the elevation of non-grid points involved in the replacement; the weight is determined based on the distance from any point to the final mosaic line. When the distance is greater than or equal to the feathering parameter expansion distance, the weight is 1; when the distance is less than the feathering parameter expansion distance and greater than or equal to the negative value of the contraction distance, the weight is determined by the weight calculation function; when the distance is less than the negative value of the contraction distance, the weight is 0; the final elevation of the arbitrary point is calculated based on the weight, and the data replacement is completed.

[0052] Step 107: Obtain the three-dimensional coordinates of the verification points through field measurements, extract the elevation of the corresponding verification points in the updated DEM, calculate the elevation error of the verification points, and verify the accuracy of the updated DEM based on the elevation error.

[0053] In this step, when performing the accuracy verification, flat land, hills, mountains, and high mountains are selected as accuracy verification areas. Within the accuracy verification area, verification points are selected for field measurements to obtain the three-dimensional coordinates of the verification points. After calculating the elevation error, it is determined whether the elevation error meets the preset mean square error requirement and does not exceed the maximum limit. If it meets the requirement, the updated DEM accuracy is determined to be qualified. If it does not meet the requirement, the aforementioned steps are repeated until the updated DEM accuracy is qualified.

[0054] As a preferred example, see Figure 1 A method for fast DEM updating based on minimum elevation difference mosaicking includes the following steps:

[0055] Step 201: Data Acquisition and Processing

[0056] (1) Natural resource monitoring data

[0057] ①Annual land change survey data

[0058] The data is obtained from the National Land Use Change Survey issued by the Ministry of Natural Resources. This data is collected and implemented annually by the national natural resources system's survey and monitoring departments. The survey monitors changes in land use status across the country annually, extracting land category change information through comprehensive remote sensing monitoring. Utilizing existing data and combining it with routine changes and relevant natural resources management findings, a unified survey base map is created. Field investigations and evidence collection are conducted, and operational information is reported. The county-level land use status database is then updated. After review and approval by the county, city, province, and ministry levels, the annual national land use change survey results are formed. The data comprehensively captures the annual changes in land categories, areas, attributes, and related individual layer information, possessing advantages in accuracy, reliability, and authenticity.

[0059] ② Routine monitoring patch data

[0060] Using the previous year's annual land change survey images as the preceding time phase and the latest monitoring images acquired in real time as the following time phase, guided by the needs of natural resource business management, focusing on farmland protection, land law enforcement, reclamation of paddy fields, ecological restoration and other businesses, monthly change patches of different business types are extracted from the monitoring images covering the entire province, and timely pushed to the natural resource authorities at the city and county levels to carry out relevant work.

[0061] ③ Other data

[0062] Acquire urban territorial spatial monitoring data (formerly known as geographic national conditions monitoring data before 2022). This data is based on high-resolution remote sensing images acquired annually and the latest relevant thematic data, combined with field surveys and other work, to organize and carry out urban territorial spatial monitoring work in prefecture-level cities and above across the country.

[0063] (2) Existing DEM data of the region

[0064] The acquisition of Guangdong Province Digital Elevation Model (referred to as "historical DEM") provides topographic elevation data for Guangdong Province, including 30-meter, 12.5-meter, and higher resolution data.

[0065] Step 202: Rapid extraction of DEM variation areas

[0066] (1) Establish rules for extracting DEM change areas

[0067] Topographic elevation changes are influenced by both natural processes (such as crustal movement, erosion and deposition) and human activities. Among these, the choice of land use type and the intensity of development are key anthropogenic factors driving elevation changes. The main land use types that are prone to causing topographic elevation changes include: mining land (industrial land, mining land), construction land (residential, commercial, public facilities, etc.), transportation land (highways, railways, airports, etc.), agricultural land (facility agricultural land), and water conservancy project land (hydraulic engineering land, ditches, etc.). These land use types drive topographic elevation changes through direct filling and excavation (such as mining and construction) or indirect erosion-deposition (such as water conservancy). Based on the results of natural resource monitoring data, the land use patches were analyzed to identify various land use types that easily cause changes in topography and elevation. A mapping relationship was established between land use patches in the natural resource monitoring data and the aforementioned land use types for two periods. Land use patches that conform to the mapping relationship between the two periods were extracted. For example, if a land use patch in the previous period was not subject to construction activity, it may become a land use patch subject to construction activity in the current period. This preliminarily defines the initial extent of the topographic changes caused by the changes from the previous period to the current period. Figure 2 As shown.

[0068] (2) Extraction of DEM variation areas

[0069] Based on the established DEM change area extraction rules, the range of map patches that cause topographic changes over multiple years is obtained. Adjacent land patches of different types are merged or fused. Then, based on the change extraction results, the number of map patches extracted with different areas and the corresponding area proportions are statistically analyzed. The proportion of map patches with an area of ​​2000 square meters is more than 80%. Combined with actual production work, map patches with an area of ​​not less than 2000 square meters are selected to obtain the DEM to be updated area of ​​this invention.

[0070] Step 203: Creating a DEM of the changing region based on dense point cloud matching

[0071] Using existing high-resolution aerial imagery and aerial triangulation results, point cloud data is generated through dense matching technology. Ground points are extracted by filtering. The point cloud data of the DEM change area formed in step 202 is processed by outward cropping, ground point editing, etc., to produce the DEM result of the change range.

[0072] After classifying the point cloud, if there are local areas where ground points cannot be obtained due to shadows, dense vegetation, etc., or if the ground point density cannot represent the ground model, the ground in that area can be fitted by adding auxiliary points based on existing or surrounding ground points and vegetation height. In addition, since there are more than 3 overlapping aerial images, if the densely matched point cloud has problems such as edge matching and ground feature time difference due to overlapping areas of different models, time difference, etc., a small-scale aerial triangulation area can be constructed and densely matched point cloud again when densely matching point clouds.

[0073] Step 204: Two-phase DEM model matching based on three algorithms: SiftGPU, Hough Transform, and UMEYAMA

[0074] As a special type of image data, DEM possesses rich geometric and textural features. Based on these features, feature points can be quickly extracted from its texture, and these feature points can be used to complete the matching process of two-phase DEM models, thereby achieving the best fit in planar positions. The Scale-Invariant Feature Transform (SIFT) algorithm is a classic algorithm used for feature detection and description in the fields of image processing and computer vision. Its main principle is to search for feature points in multi-scale space and calculate the orientation information of the feature points. SIFTGPU technology significantly improves the processing efficiency of the SIFT algorithm by leveraging the parallel computing capabilities of graphics processing units (GPUs), enabling rapid feature point extraction and matching. Hough transform is one of the fundamental methods in image processing for identifying geometric shapes within an image. Its main function is to separate geometric shapes (such as lines) with similar characteristics from an image. In practical applications, it can be used to eliminate gross errors generated during feature point extraction. The UMEYAMA algorithm is used to calculate the optimal coordinate system transformation between two sets of point clouds or trajectories. It solves for the optimal transformation of the rotation matrix, translation vector, and scale factor using the least squares method, aligning two corresponding point sets to form a registration point pair. The specific implementation is as follows:

[0075] Using historical DEMs as a reference, the original data of each DEM is divided into small blocks according to a set size. Based on the location information of each small block, its corresponding DEM block to be matched, formed in step 203, is found. Then, feature points are extracted and described using the SIFTGPU algorithm to determine the initial matching points. Finally, the Hough transform method is used to remove gross errors. Then, based on the pixel coordinates of all corresponding points, the object-space 3D coordinates of each pixel are obtained from the DEM. Based on the two sets of 3D coordinates, the UMEYAMA algorithm is used to calculate the transformation parameters between the new DEM and the historical DEM. The DEM from step 203 is then transformed according to the transformation parameters to achieve spatial alignment between the two DEMs. Figure 3 As shown.

[0076] (1) Registration of DEM data based on SIFTGPU algorithm

[0077] DEMs are regular raster data, with grayscale values ​​corresponding to different terrain elevations. The classic image matching SIFTGPU algorithm can be used for matching. Furthermore, because DEMs themselves possess initial positioning information, this information can be fully utilized to determine the initial matching range, narrowing the search area for image matching and improving the matching success rate.

[0078] (2) Elimination of matching gross errors based on parameter space Hough transform method

[0079] Using the Hough transform method for candidate matching selection falls under the second type of problem addressed by the Hough transform: problems involving two sets of features that conform to a certain spatial correspondence, such as translation transformation. Due to differences in object elevation or depth of field, the scale, rotation angle, and disparity vary across stereo image pairs. However, from the perspective of the overall image, the features at each location roughly satisfy a planar similarity transformation; only the constraint needs to be appropriately relaxed. Therefore, stereo images exhibit "overall geometric consistency accompanied by local disparity differences." To address this characteristic, the Hough transform method is used to resolve overall geometric consistency, while a rigorous imaging relationship model is employed to resolve local disparity differences. Since overall geometric consistency is a global constraint between features, an error matching removal strategy based on global geometric consistency constraints is used to eliminate gross errors in the initial matching points, resulting in the final registration points.

[0080] (3) Calculate spatial similarity transformation parameters based on the UMEYAMA algorithm

[0081] In a 3-dimensional space, there are two sets of points, each containing There are 3 points, which are arranged into a matrix as follows:

[0082]

[0083]

[0084] Corresponding points satisfy a similarity transformation relationship, that is: To calculate the scale Rotation matrix Horizontal displacement This makes the total error Minimum:

[0085]

[0086] The error can be expressed in matrix form as follows:

[0087]

[0088]

[0089] The UMEYAMA algorithm uses the conclusion of the optimal rotation matrix between two matrices to solve the above problem: for The matrix, for The matrix, The singular value decomposition is . for The rotation matrix makes Minimal rotation matrix satisfy .in, for A matrix that satisfies the following two conditions:

[0090]

[0091] and:

[0092] (4) Planar matching of two-phase DEM models

[0093] Using the registration points generated in the above steps, and taking the historical DEM data as a reference, the newly generated DEM data in step 203 is subjected to planar matching to obtain the DEM data that best coincides with the historical DEM data in planar position after registration.

[0094] Step 205: Automatic generation and updating of DEM mosaic lines based on the minimum cost path algorithm of elevation difference weighted map

[0095] The A* algorithm is a commonly used, efficient shortest path algorithm in rasterized maps. It estimates the cost function and traverses the nodes to obtain the shortest path from the starting point to the ending point, and the path is not a closed line. The goal of this invention is to find the location with the smallest elevation difference between two DEM data periods, connect these locations to form a closed line, and use this as the mosaic line of the DEM update area.

[0096] The operation process is as follows: First, a buffer of a certain width is set for the range of changes extracted in step 202. Second, the elevation difference at each pixel location is calculated within the buffer range and used as a reference weight for subsequent processing. Based on this, and combined with the two-phase matching DEMs in step 204, a closed mosaic line is automatically generated on the weight map using the minimum cost path algorithm, according to the reference weight map. This mosaic line is formed by connecting the elevation points with the smallest elevation difference between the two phases of DEM data. Finally, the automatically generated mosaic line is manually adjusted and appropriately feathered to form the final mosaic line, thereby achieving a smooth transition of elevation between the two phases of DEM data.

[0097] (1) Generate a buffer

[0098] Set appropriate buffer parameters to generate a buffer around the original patch, which is used to determine the area of ​​the DEM and historical DEM data mosaic after matching in step 204.

[0099] (2) Calculate the cost matrix

[0100] Set up a cost matrix of the same size as the elevation matrix, and fill it with the absolute value of the elevation difference between the two DEM periods. Then set up a gradient difference matrix of the same size, normalize it and the cost matrix to the same range, iterate through the two matrices, and fill the cost matrix with the larger value.

[0101] (3) Find the path with minimum cost and automatically close to form a mosaic line.

[0102] In the cost matrix, the middle row is taken, dividing the matrix into upper and lower parts. Within the middle row, the two points with the lowest costs in the left and right halves of the matrix are selected. The A* algorithm, which combines cost calculation with the minimum cost calculation, is used to find the minimum-cost paths connecting these two points from the upper and lower halves of the matrix. The A* algorithm is a method for finding the shortest optimal path in a static grid; here, the cost matrix is ​​used to represent the path calculated through the cost of each grid cell. The two minimum-cost paths are then concatenated to form a closed line, which represents the mosaic line with the smallest difference in elevation between the two DEM phases.

[0103] (4) Smooth closed line

[0104] The closed lines obtained by splicing are smoothed to reduce unnecessary turns and fluctuations. The smoothed closed lines are the final inlay lines.

[0105] (5) Replace data

[0106] A portion of the DEM data registered in step 203 is selected and replaced with historical DEM data. The replacement principle can be calculated using a formula.

[0107]

[0108] Expression, where the function Representing data exist Elevation of the point Represents any point, for The final elevation used for the point This is old data. For new data, As weight, based on point The calculation covers the distances inside and outside the inlay line, as well as the distances to the inlay line, within the range of... Formulas can be used

[0109]

[0110] Expression, in which Point The distance to the tessellation line is defined as negative inside the tessellation line and positive outside the line. The outer distance of the feathering parameter. For the indentation distance, the function This is the weight calculation function. Without feathering, it is assumed... and If both are 0, then Outside the inlay line A value of 1 indicates that only the old data is used; a value of 0 indicates that only the new data is used within the mosaic lines. (Function) Retrieve data exist When determining the elevation of a point, if Point not in data When the grid points are on the grid, the elevation is obtained using bilinear linear interpolation.

[0111] Before updating, a smoothing scale is set for the mosaic lines so that the DEM data registered in step 204 can be continuously and smoothly updated into the historical DEM data, resulting in an updated DEM with a natural boundary transition.

[0112] By using the cost A* algorithm, the minimum cost path connecting two points in the upper and lower halves of the matrix is ​​found from the two DEMs respectively, generating the minimum elevation difference mosaic line. This allows the updated DEM to smoothly transition in areas where the elevation has not changed, achieving minimum elevation difference fitting and realizing seamless and efficient DEM mosaicking.

[0113] Step 206: Verification of the accuracy of DEM update results

[0114] Based on field-measured elevation points, the elevation accuracy of exposed surface areas in the updated DEM results is verified. The specific method is as follows:

[0115] (1) Select bare surface areas of different terrain types as accuracy verification areas, including plains, hills, mountains and high mountains.

[0116] (2) Within the accuracy verification area, select verification high points located on flat ground, bare ground, or other elevation changes, and conduct field measurements to obtain their precise three-dimensional coordinates (x1, y1, z1).

[0117] (3) Display the field measurement points on the map according to the three-dimensional coordinates, extract the elevation of the ground features at the same plane position in the updated DEM results, and obtain their three-dimensional coordinates (x1, y1, z2).

[0118] (4) Calculate the elevation error of the elevation verification points in the DEM update results, and perform elevation accuracy analysis according to different terrain regions. The mean error of the elevation error in the DEM update results should meet the corresponding accuracy requirements, and the elevation error should not exceed the maximum limit, i.e.: ±|z1-z2|<±2m0 (where m0 is the mean elevation error in the technical specifications).

[0119] The following explanation uses specific data.

[0120] Step 301: Data Acquisition and Preprocessing

[0121] (1) Natural resource monitoring data

[0122] ① Obtain historical land use change survey data for Guangdong Province. Analyze the land use category codes, names, areas, and other relevant attributes of land use changes in each year to determine the changes in map features in each year.

[0123] ② Collect routine monitoring patch data, analyze the relevant attributes such as patch type, code, and area in each period, and determine the changes in patches in each period.

[0124] ③ Obtain urban land space monitoring data, analyze the land type, area and other related attributes of the monitoring data for each year, and determine the changes in map patches.

[0125] (2) Existing DEM data of the region

[0126] The results of the Guangdong Provincial Digital Elevation Model (referred to as "historical DEM") were obtained. The coordinate system of this data is the CGCS2000 National Geodetic Coordinate System, the 1985 National Elevation Datum, and the grid size is 2 meters × 2 meters. The data is divided into 5km × 5km sections, and the divided data need to be merged.

[0127] Step 302: Rapid extraction of DEM variation areas

[0128] The land use patches obtained in step 301(1) are analyzed. If the update frequency is relatively high (once a year or several times a year), the non-construction land use patches that have been converted to construction land use in the previous year or period can be extracted to form the range of patches that caused topographic changes from the previous year or period to the current year or period. Table 1 shows the types of land use patches that caused topographic changes. If the update frequency is relatively low (once every five or ten years), the patches that have been converted from non-construction land use to construction land can be extracted year by year. The extracted patches are compared with the patches of the first year, and patches of the same land use type are removed to finally obtain the range of patches that caused topographic changes over multiple years. Finally, the patches of adjacent different land use types are merged or integrated, and patches with an area of ​​not less than 2000 square meters are selected to obtain the DEM area to be updated in this invention.

[0129] surface Topographical changes and land cover types

[0130]

[0131] Compared to extracting the range of terrain changes by using the elevation difference between two DEMs, this method can significantly improve the accuracy and update efficiency of terrain change detection, and is more accurate, reliable and realistic, effectively reducing production time and costs.

[0132] like Figure 4As shown, the location of changed patches obtained through the land use change extraction rules is relatively complete, and even small areas of patches can be extracted.

[0133] Step 303: Creating a DEM of the changing regions based on dense point cloud matching

[0134] Aerial triangulation results of aerial imagery are densely matched to obtain image point clouds. The DEM to be updated area formed in step 302 is expanded outward, and the point cloud data is cropped based on the expanded range. Then, the cropped point cloud is filtered, deredundantized, and classified. Finally, point cloud interpolation is performed to generate DEM data of the changed area of ​​this invention. For cases where ground points cannot be obtained in local areas due to shadows, the ground of that area can be fitted by adding auxiliary points according to the region using existing or surrounding ground points, such as... Figure 5 As shown.

[0135] This method offers advantages in point cloud creation, including guaranteed accuracy and low production costs. It fully utilizes existing aerial imagery to generate large-area point cloud data, significantly reducing production costs compared to aerial point cloud acquisition. Furthermore, the point cloud data obtained through dense matching of aerial imagery has a much higher accuracy than point clouds generated from stereo satellite imagery, ensuring data precision. The locally updated DEM accuracy meets the requirements for provincial-scale DEMs. While DEM creation based on densely matched point clouds suffers some loss in accuracy compared to acquiring elevation points from stereo image pairs, it still meets the accuracy requirements, greatly reducing workload, improving efficiency, and ensuring strong timeliness.

[0136] Step 304: Two-phase DEM model matching based on three algorithms: SiftGPU, Hough Transform, and UMEYAMA

[0137] This invention employs the classic SIFTGPU algorithm for two-stage DEM model matching. Based on the SIFTGPU algorithm, feature points are extracted and described from the data to determine initial matching points. Finally, the Hough transform method is used to remove gross errors. Then, based on the pixel coordinates of all corresponding points, the object-space 3D coordinates of each pixel are obtained from the DEM. Next, based on the two sets of 3D coordinates, the UMEYAMA algorithm is used to calculate the transformation parameters between the new and historical DEMs. The new DEM is then subjected to coordinate transformation based on these parameters to achieve spatial alignment between the new and old DEMs.

[0138] Due to differences in production technology and methods, the two DEM data sets exhibit varying accuracy. Direct and rapid updates can lead to deviations in planar position and elevation. This invention, based on the SIFTGPU algorithm, incorporates the UMEYAMA algorithm and Hough transform method to optimize the geometric matching model, resolving the geometric position deviation problem and improving the planar position accuracy of the updated DEM. Figure 6As shown, in some areas, the planar position deviation between the two DEM planar models was about 10 meters before matching, but after matching, the planar position deviation was almost zero.

[0139] Step 305: Automatic generation and updating of DEM mosaic lines based on the minimum cost path algorithm of elevation difference weighted map

[0140] A buffer zone of a certain width is set for the area of ​​the patches extracted in step 302. The elevation difference at each pixel location within the buffer zone is calculated and used as a reference weight for subsequent processing. A mosaic line is automatically generated on the weight map using the minimum cost path algorithm. This mosaic line is the line connecting the elevation points with the smallest elevation difference between the historical DEM and the DEM data matched in step 304. The automatically generated mosaic line is then manually adjusted and appropriately feathered to form the final mosaic line, thereby achieving a smooth transition of elevation between the two DEM data periods.

[0141] like Figure 7 As shown, the blue boundary line is a mosaic line drawn manually, with a maximum height difference of 2.5 meters between the two models. The red boundary line is a mosaic line with the minimum height difference generated by automatically finding the minimum cost path between the two DEMs using the cost A* algorithm, with a maximum height difference of 0.8 meters.

[0142] Step 306: Update accuracy verification

[0143] For the updated DEM data obtained in step 305, 325 points with elevation changes, mainly located in flat areas and exposed surfaces, were selected. Elevation accuracy was verified using field measurement points. The elevation accuracy meets the requirements of 1:2000 DEM in the "Technical Specifications for Production of Basic Geographic Information Digital Results 1:500, 1:1000, and 1:2000 Part 2: Digital Elevation Model" (CH / T 9020.2-2013). The statistical results are shown in Table 2 below.

[0144] Table 2 Statistical Table for Verification of Elevation Accuracy in Field Measurements

[0145]

[0146] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the present invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of different embodiments or examples.

[0147] The above embodiments are merely illustrative of the technical concept and features of the present invention, and are intended to enable those skilled in the art to understand the content of the present invention and implement it accordingly. They should not be construed as limiting the scope of protection of the present invention. All equivalent changes or modifications made based on the essence of the content of the present invention should be covered within the scope of protection of the present invention.

Claims

1. A DEM fast updating method based on minimum elevation difference tessellation, characterized in that, Includes the following steps: Acquire natural resource monitoring data and existing regional DEM data; DEM change area extraction rules are established based on the natural resource monitoring data. Land use patches in the natural resource monitoring data are analyzed according to the extraction rules. Adjacent land use patches of different types are merged or integrated. Patches that meet the area requirements are selected to obtain the DEM area to be updated. Point cloud data is generated using aerial triangulation results from aerial imagery. Ground points are extracted from the point cloud data by filtering. The point cloud data is then expanded, cropped, and ground points are edited according to the DEM to be updated area to obtain the changed area DEM. Feature matching and gross error removal are performed on the existing DEM of the region and the DEM of the changed region. The transformation parameters of the DEM of the changed region and the existing DEM of the region are calculated based on the processed matching points. The planar position transformation of the DEM of the changed region is performed according to the transformation parameters to obtain the DEM after planar matching. A buffer is set for the DEM to be updated area. The elevation difference between the existing DEM and the plane-matched DEM in the buffer is calculated as a reference weight for subsequent processing. Based on the reference weight map, the minimum cost path of the elevation difference between the two DEMs is found. The minimum cost path is spliced ​​to form a closed line. The closed line is smoothed to obtain the final mosaic line. Based on the final mosaic line, the existing DEM and the plane-matched DEM in the region are replaced according to the weight calculation rules to obtain the updated DEM.

2. The DEM fast updating method based on minimum elevation difference tessellation of claim 1, wherein, After obtaining the updated DEM, the three-dimensional coordinates of the verification points are obtained through field measurements. The elevations of the corresponding verification points in the updated DEM are extracted, the elevation errors of the verification points are calculated, and the accuracy of the updated DEM is verified based on the elevation errors.

3. The DEM fast updating method based on minimum elevation difference tessellation of claim 1, wherein, When acquiring the natural resource monitoring data, annual land change survey data, routine monitoring patch data, and urban land space monitoring data are acquired. Attribute analysis is performed on the annual land change survey data, the routine monitoring patch data, and the urban land space monitoring data respectively to obtain the natural resource monitoring data containing land category change information. When acquiring existing DEM data for the region, digital elevation model data containing different resolutions is acquired, and the digital elevation model data is merged to obtain existing DEM data for the region.

4. The DEM fast updating method based on minimum elevation difference tessellation of claim 1, wherein, When establishing the DEM change area extraction rules, the land use types that are prone to causing changes in terrain elevation are analyzed, the mapping relationship between land use patches and land use types in the natural resource monitoring data is sorted out, and the DEM change area extraction rules are established based on the mapping relationship.

5. The DEM fast updating method based on minimum elevation difference tessellation of claim 1, wherein, After generating the point cloud data, if ground points cannot be obtained or the ground point density is insufficient in a local area, auxiliary points are added according to the region using existing or surrounding ground points and vegetation height to fit the ground of the local area; if the point cloud has edge confluence or time difference issues with ground features, a small-scale aerial triangulation region is constructed by selecting the latest photo, and the point cloud data is re-generated by dense matching of the small-scale aerial triangulation region. Then, the point cloud data is expanded, cropped, and the ground points are edited to obtain the DEM of the changed region.

6. The DEM fast update method based on minimum elevation difference mosaicking according to claim 1, characterized in that, When performing feature matching on the existing DEM of the region and the changed DEM of the region, feature points are extracted using the SiftGPU algorithm. The matching range is determined using the initial positioning information of the existing DEM of the region and the changed DEM of the region. The feature points are described within the matching range to determine the initial matching point.

7. The DEM fast update method based on minimum elevation difference mosaicking according to claim 1, characterized in that, When performing gross error removal on the initial matching points, Hough transform is used to filter the initial matching points based on the overall geometric consistency constraint of the stereo image, and matching points that do not conform to the overall geometric consistency are removed to obtain matching points after gross error removal. When calculating the transformation parameters, the two sets of object-side 3D coordinates corresponding to the matching points after removing gross errors are arranged into matrices using the UMEYAMA algorithm. Singular value decomposition is performed on the matrices, and the scale, rotation matrix, and translation are solved based on the decomposition results to obtain the transformation parameters.

8. The DEM fast update method based on minimum elevation difference mosaicking according to claim 1, characterized in that, When setting the buffer, the buffer parameters are determined, and the buffer is generated around the DEM area to be updated. When finding the minimum cost path based on the elevation difference, the elevation difference is used as a reference weight to generate a cost matrix, and a gradient difference matrix of the same size as the cost matrix is ​​generated. The cost matrix and the gradient difference matrix are normalized. The two normalized matrices are traversed, and the larger value is taken to update the cost matrix. The minimum cost path is found by combining the updated cost matrix with the A* algorithm.

9. The DEM fast update method based on minimum elevation difference mosaicking according to claim 1, characterized in that, When performing data replacement according to the weight calculation rules, bilinear interpolation is used to obtain the elevation of non-grid points involved in the replacement; the weight is determined based on the distance from any point to the final mosaic line; when the distance is greater than or equal to the feathering parameter expansion distance, the weight is 1; when the distance is less than the feathering parameter expansion distance and greater than or equal to the negative value of the shrinkage distance, the weight is determined by the weight calculation function; when the distance is less than the negative value of the shrinkage distance, the weight is 0; the final elevation of the arbitrary point is calculated based on the weight to complete the data replacement.

10. The DEM fast update method based on minimum elevation difference mosaicking according to claim 2, characterized in that, When performing the accuracy verification, flat land, hills, mountains, and high mountains are selected as accuracy verification areas. Within the accuracy verification areas, verification points are selected for field measurements to obtain the three-dimensional coordinates of the verification points. After calculating the elevation error, it is determined whether the elevation error meets the preset mean square error requirements and does not exceed the maximum limit. If it meets the requirements, the updated DEM accuracy is determined to be qualified. If it does not meet the requirements, the aforementioned steps are repeated until the updated DEM accuracy is qualified.

Images