An automated large-scale DEM differential correction method based on adaptive data segmentation

Abstract

The invention discloses an automatic large-scale DEM differential correction method based on adaptive data segmentation, comprising the following steps: step 1: initial DEM differential processing; step 2: rough overall correction processing; step 3: adaptive data segmentation; step 4: accurate local correction processing and data splicing, so as to automatically complete the error correction in the large-scale DEM differential map; by taking into account the DEM error characteristics and the local heterogeneity characteristics of the data of the DEM differential map, a universal DEM differential post-processing correction method that does not require manual supervision and is suitable for large-scale multi-source error aliasing complex scenes is developed. The invention not only improves the data processing efficiency of the DEM differential method, but also improves the accuracy of land surface elevation change monitoring.

Description

Technical Field

[0001] The invention relates to the technical field of digital elevation model in land surface change monitoring, and in particular to an automated large-scale DEM differential correction method based on adaptive data segmentation. Background Art

[0002] Digital Elevation Model (DEM) is one of the most basic and important data in the field of earth and planetary science, and has been widely used in many aspects. In the past 20 years, with the rapid development of earth observation technology (i.e. satellite remote sensing), the ability to obtain satellite data and prepare DEM data products has been revolutionary. In 2000, the United States launched the Shuttle Radar Topography Mission (SRTM), which used radar interferometry technology to obtain digital elevation models of more than 80% of the global land surface in 11 days, opening the era of global DEM product preparation. Since then, the space agencies of Japan and the United States have used Terra-ASTER and ALOS-PRISM optical satellite sensors and prepared and released ASTERGDE M and AW3D DEM data products based on photogrammetry technology. In 2010, the German Aerospace Agency launched a five-year TanDEM high-resolution global mapping program, and for the first time used a single-station mode spaceborne radar interferometry system to complete the update and preparation of high-precision global DEM data products. In 2022, my country also launched the L-band Land Exploration No. 1 radar satellite system (consisting of two satellites), opening the era of global mapping in my country. In addition, with the rapid development of sub-meter ultra-high-resolution optical satellites (such as Pléiades, WorldView and GeoEye series, etc.) and near-ground unmanned aerial vehicle systems in recent years, high-resolution local area DEM data has been widely used in engineering construction, disaster prevention and mitigation, surface changes and other fields.

[0003] With the acquisition of multi-source remote sensing data from air and space and the preparation of a large number of DEM data, quantitative estimation of surface height changes can be achieved by differential calculation using two or more DEM data, and then various information such as the spatiotemporal pattern, volume or mass changes of different surface types can be obtained, providing scientific guidance for government decision-making, regional management, disaster prevention and control, national defense security, engineering construction, etc. In particular, under the current background of global climate change and increased human activities, the height of the surface has also undergone profound and significant changes in the past few decades. At present, the DEM difference method (also known as geodesy) has been widely used in land surface change monitoring, such as glacier changes, thermal melt landslides, river valley erosion, earthquakes, volcanic eruptions, forest changes and many other fields.

[0004] For DEM data products generated by satellite data, there are generally three types of errors: gross errors, strip errors, and spatial trend errors. For gross errors, they are generally related to incorrect phase unwrapping in the radar interferometry process (for radar sensors) and incorrect image matching in optical photogrammetry (for optical sensors). Strip errors are generally mainly affected by the disturbance of the flight attitude of the optical satellite when the original satellite image is acquired (i.e., pitch, roll, and tilt), especially for push-broom optical imaging systems (such as SPOT-5 / 6 / 7 and Terra-ASTER satellites). In addition, uncorrected splicing of multiple data sometimes leads to similar error forms (such as KH-9 satellites). For spatial trend errors, they are generally related to the orbital errors of radar sensors. Studies have shown that the original data products of SRTM DEM still have long-distance trend errors caused by radar antenna disturbances. Although the subsequent updated products have corrected and processed this error, from a large number of evaluation experiments, the latest high-quality DEM products still inevitably have multi-scale spatial trend errors. In addition, for optical satellite sensors, the non-uniform distribution of ground control points will also lead to the emergence of spatial error terms.

[0005] For the processing of these errors, in the current study, when faced with a small area and a small amount of DEM differential processing, gross errors can generally be removed by manually setting masks. Strip errors can be removed along a specific direction by fitting a high-order polynomial. Spatial trend errors can also be removed by fitting a polynomial of appropriate order. In addition, existing studies have shown that there is always a local offset between different DEM data products, which is mainly related to different data processing and coordinate system definitions. Therefore, before performing the above correction processing, DEM data registration is also very necessary.

[0006] For publicly released global DEM data products (such as AW3D30, SRTM DEM, Multi-Error-RemovedImproved-Terrain (MERIT) DEM, TanDEM-X 90 DEM, Copernicus DEM, etc.), due to the characteristics of their data fusion, multi-source errors have aliasing effects, which makes the above errors more complex in DEM data, bringing great challenges to the subsequent correction processing. Faced with this problem, on the one hand, it is obviously impossible to achieve high data processing accuracy with only one processing. On the other hand, the traditional method based on manual supervision correction is no longer suitable for processing large amounts of DEM data. Therefore, it is necessary to develop a more general and automated DEM differential data processing method. Summary of the invention

[0007] The present invention aims at the defects of the prior art and provides an automatic large-scale DEM differential correction method based on adaptive data segmentation.

[0008] In order to achieve the above invention object, the technical solution adopted by the present invention is as follows:

[0009] An automated large-scale DEM differential correction method based on adaptive data segmentation includes the following steps:

[0010] Step 1: Data preprocessing: Unify the coordinate system reference of two or more different DEM data, then perform direct difference processing on the DEM data to obtain a DEM difference map, and then combine the surface type boundary data to mask the potential change area to obtain the initial elevation change map of the stable area;

[0011] Step 2: Overall data rough correction processing: First, check and correct the large-scale spatial trend errors of the elevation change map in the stable area; then, on this basis, perform initial data rough registration on the elevation difference map;

[0012] Step 3: Adaptive data segmentation: First, the terrain slope information is used as a reference, and the elevation difference map after overall correction is initially segmented according to the quadtree segmentation technology; then, the slices with strong local heterogeneity are further segmented with the standardized standard deviation as a reference, so as to obtain the final data segmentation solution;

[0013] Step 4: Local data fine correction processing: According to the obtained data segmentation scheme, the elevation difference map after rough correction is segmented, and the data slices are first spatially detrended and corrected one by one, and then refined data registration is performed;

[0014] Step 5: Splice the block-corrected data to generate the final corrected DEM difference map.

[0015] Furthermore, the step 1 includes the following sub-steps:

[0016] Step 11: Convert the DEM ellipsoidal height with the earth ellipsoid as the vertical reference to the orthometric height data with the geoid as the reference;

[0017] Step 12: Use appropriate sampling methods to unify the resolution of DEM data;

[0018] Step 13: Perform differential processing on DEM data of the same size and in the same geographical range to obtain a DEM differential map;

[0019] Step 14: Mask the area with surface elevation changes to obtain the initial elevation difference map of the stable area;

[0020] Step 15: Remove outliers from the initial elevation difference map.

[0021] Furthermore, the step 2 specifically includes the following sub-steps:

[0022] Step 21: Use the elevation difference map of the stable area to perform initial overall spatial trend removal;

[0023] Step 22: Calculate terrain slope and aspect information using one of the DEMs in the differential processing as a reference;

[0024] Step 23: The elevation difference map after overall detrending is roughly aligned with the overall data using the slope and aspect information.

[0025] Furthermore, in step 2, formula (1) is used to simulate the spatially varying trend error.

[0026] dh(x,y)=a1+b1x+c1y+d1xy+e1x 2 +f1y 2 (1)

[0027] Among them, dh(x, y) represents the elevation difference of the point with coordinates (x, y) in the data. a1, b1, c1, d1, e1 and f1 are the parameters of the overall spatial trend to be estimated. The entire estimation process is calculated using the iterative least squares method.

[0028] Furthermore, the data rough registration described in step 23 is specifically as follows: an analytical three-dimensional registration method based on slope and aspect information is used to perform overall data registration on the obtained DEM difference map, and a relatively rough registration accuracy is set to achieve preliminary alignment of the data. This process can be specifically expressed as formula (2).

[0029]

[0030] Among them, slope and aspect represent the terrain slope and slope direction of a point in the data, respectively. A and B represent the magnitude and direction of the potential horizontal displacement between DEMs, respectively. C represents the overall vertical average deviation between DEMs, expressed as In addition, flag is a Boolean data, that is, yes or no, which is used to determine whether iterative processing is required during data processing. R is the spatial resolution, (A×cos(B)) / R and (A×sin(B)) / R represent the number of pixels moved in the north and east directions respectively. Core thre The threshold value of the number of pixels that need to be moved, i.e., the iteration index.

[0031] Furthermore, the step 3 specifically includes the following sub-steps:

[0032] Step 31: Taking the terrain slope as a reference, set an appropriate slope threshold, and use the quadtree segmentation technology to perform initial data segmentation on the elevation difference map;

[0033] Step 32: Calculate the standardized standard deviation of the local elevation difference map for the data slices whose size is larger than the minimum segmentation size;

[0034] Step 34: Using the standardized standard deviation as a reference, set a suitable threshold, and segment the data units above the threshold again, and segment them to the minimum segmentation size, so as to obtain the final data segmentation solution;

[0035] Furthermore, the step 31 is specifically as follows: using quadtree decomposition technology to perform multi-scale data segmentation on the DEM difference map obtained after rough registration. This process is specifically expressed as formula (3) and formula (4).

[0036]

[0037]

[0038] Where I represents the DEM difference map to be processed, Represents the quadtree segmentation operator. S(i, j, k) represents the quadtree structure of I, (i, j) is the coordinate of the upper left corner of the data slice after segmentation, and k represents the segmentation level (k∈N 0 The size of the segmented slice is M×M pixels (M=2 (m-k) ), of which 2 m and 2 n Respectively represent the maximum and minimum size of a slice during the cutting process. QTseg_or_not is a Boolean data. If the average slope of the terrain within the data slice is Greater than the slope threshold Slope thre , then further data segmentation is required; otherwise, no further segmentation is required.

[0039] Furthermore, the re-segmentation in step 34 is specifically as follows:

[0040] The standardized standard deviation SSTD is used to evaluate the spatial heterogeneity of each data slice. When SSTD is greater than a set data heterogeneity threshold Heter thre , then further data segmentation is required, and the data is directly segmented to the smallest slice size; otherwise, no segmentation is performed. This process is expressed as formula (5) and formula (6).

[0041]

[0042]

[0043] where dh p , dh p_max , dh p_min and N represent the elevation difference, maximum elevation difference, minimum elevation difference and total number of pixels in a data slice (S(i, j, k)) respectively. Reseg_or_not is a Boolean data used to determine whether the data needs to be segmented again.

[0044] Furthermore, the step 4 specifically includes the following sub-steps:

[0045] Step 41: According to the final data segmentation scheme, perform spatial trend correction again for each segmentation unit;

[0046] Step 42: According to the final data segmentation scheme, perform data precise registration again for each segmentation unit;

[0047] Step 43: The segmented and re-corrected data are spliced ​​according to the original positions, and the final corrected elevation difference map is output.

[0048] Further, step 41 is specifically as follows: segment and slice the data according to the final segmentation scheme, then simulate the linear plane trend error and the error along the height direction for each data slice one by one, and subtract the fitting error from the original data to complete the removal of the local trend term error. The simulated error is estimated using formula (7).

[0049] dh(x,y)=a2+b2x+c2y+d2n(x,y) (7)

[0050] Among them, a2, b2, c2 are the coefficients of the local spatial trend to be calculated, and d2 is the coefficient related to the height trend error, which can be estimated by iterative least squares method. η(x, y) represents the elevation information of a certain point.

[0051] Compared with the prior art, the advantages of the present invention are:

[0052] 1. The automated large-scale DEM differential correction method based on adaptive data segmentation of the present invention fully considers the error characteristics and local heterogeneity characteristics in the DEM differential map, can realize fully automated large-scale and large-batch DEM differential processing, and is suitable for correction of various error types, thereby improving the accuracy of land surface elevation change monitoring.

[0053] 2. The present invention not only effectively solves the current difficulty of large-scale DEM differential processing under complex error modes, but also improves the computational efficiency. In addition, the algorithm of the present invention is simple and applicable to all error types, and can also be extended to the processing of time series DEM differentials, thus helping to promote and promote the application of DEM differential technology in various fields. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 A schematic diagram of a flow chart of an embodiment of the present invention;

[0055] Figure 2 The experimental results of differential processing using SRTM-GL1 and SRTM-X DEM data obtained in the Tianshan area according to the embodiment of the present invention;

[0056] Figure 2 (a) is the topography and glacier boundary information map of the Tianshan Experimental Area in China;

[0057] Figure 2 (b) The initial elevation difference map is obtained based on SRTM GL1 and SRTM-X DEM data;

[0058] Figure 2 (c) is the elevation difference map after overall detrending and rough registration;

[0059] Figure 2 (d) is the statistical information of the stable area of ​​the elevation difference map before and after the overall processing;

[0060] Figure 2 (e) is the case of regular data segmentation with a size of 1024*1024 pixels (the base map is terrain slope information);

[0061] Figure 2 (f) shows the case of regular data segmentation with a size of 256*256 pixels (the base map is terrain slope information);

[0062] Figure 2 (g) is the case of data segmentation based on quadtree (the base map is terrain slope information);

[0063] Figure 2 (h) is the case based on adaptive data segmentation (the base map is terrain slope information);

[0064] Figure 2 (i) For the use of Figure 2 (e) Elevation difference map after local data correction by segmentation scheme;

[0065] Figure 2 (j) For the use of Figure 2 (f) Elevation difference map after local data correction by segmentation scheme;

[0066] Figure 2 (k) For use Figure 2 (g) Elevation difference map after local data correction by segmentation scheme;

[0067] Figure 2 (l) For the use of Figure 2 (h) Elevation difference map after local data correction by segmentation scheme;

[0068] Figure 2 (m) is Figure 2 (i) Statistical information of the corresponding stable region before and after local fine correction processing;

[0069] Figure 2 (n) Figure 2 (j) Statistical information of the corresponding stable region before and after local fine correction processing;

[0070] Figure 2 (o) Figure 2 (k) Statistical information of the corresponding stable region before and after local fine correction processing;

[0071] Figure 2 (p) is Figure 2 (l) Statistical information of the corresponding stable region before and after local fine correction processing;

[0072] Figure 3 The experimental results of differential processing using SRTM-GL1 DEM and CopernicusDEM data obtained in the central Pamir region by the embodiment of the present invention;

[0073] Figure 3 (a) is a topographic and glacier boundary information map of the experimental area of ​​the Pamir Plateau in Tajikistan;

[0074] Figure 3 (b) The initial elevation difference map is obtained based on SRTM GL1 and Copernicus DEM data;

[0075] Figure 3 (c) is the elevation difference map after overall detrending and rough registration;

[0076] Figure 3 (d) is the statistical information of the stable area of ​​the elevation difference map before and after the overall processing;

[0077] Figure 3 (e) is the case of regular data segmentation with a size of 1024*1024 pixels (the base map is terrain slope information);

[0078] Figure 3(f) shows the case of regular data segmentation with a size of 256*256 pixels (the base map is terrain slope information);

[0079] Figure 3 (g) is the case of data segmentation based on quadtree (the base map is terrain slope information);

[0080] Figure 3 (h) is the case based on adaptive data segmentation (the base map is terrain slope information);

[0081] Figure 3 (i) For the use of Figure 3 (e) Elevation difference map after local data correction by segmentation scheme;

[0082] Figure 3 (j) For the use of Figure 3 (f) Elevation difference map after local data correction by segmentation scheme;

[0083] Figure 3 (k) For use Figure 3 (g) Elevation difference map after local data correction by segmentation scheme;

[0084] Figure 3 (l) For the use of Figure 3 (h) Elevation difference map after local data correction by segmentation scheme;

[0085] Figure 3 (m) is Figure 3 (i) Statistical information of the corresponding stable region before and after local fine correction processing;

[0086] Figure 3 (n) Figure 3 (j) Statistical information of the corresponding stable region before and after local fine correction processing;

[0087] Figure 3 (o) Figure 3 (k) Statistical information of the corresponding stable region before and after local fine correction processing;

[0088] Figure 3 (p) is Figure 3 (l) Statistical information of the corresponding stable region before and after local fine correction processing. DETAILED DESCRIPTION

[0089] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and examples.

[0090] An automated large-scale DEM differential correction method based on adaptive data segmentation includes the following steps:

[0091] Step 1: Data preprocessing: unify the reference benchmark, data resolution, data geographic range and data size of two or more DEM data to be processed for DEM difference, and then perform difference operation and combine the boundary information of surface type to obtain the initial DEM difference map of the stable area;

[0092] Step 2: Perform overall spatial trend correction on the initial DEM difference map. According to the elevation difference information of the stable area (assuming that its theoretical change is zero), use formula (1) (quadratic polynomial) to simulate the trend term error of spatial change.

[0093] dh(x,y)=a1+b1x+c1y+d1xy+e1x 2 +f1y 2 (1)

[0094] Among them, dh(x, y) represents the elevation difference of the point with coordinates (x, y) in the data. a1, b1, c1, d1, e1 and f1 are the parameters of the overall spatial trend to be estimated. The entire estimation process is calculated using the iterative least squares method.

[0095] Step 3: Overall data coarse registration. The DEM difference map obtained in step (2) is registered with an analytical three-dimensional registration method based on slope and aspect information. At the same time, a relatively rough registration accuracy is set to achieve preliminary alignment of the data. This process can be specifically expressed as formula (2).

[0096]

[0097] Among them, slope and aspect represent the terrain slope and slope direction of a point in the data, respectively. A and B represent the magnitude and direction of the potential horizontal displacement between DEMs, respectively. C represents the overall vertical average deviation between DEMs, which can be expressed as In addition, flag is a Boolean data (yes or no) used to determine whether iterative processing is required during data processing. R is the spatial resolution, (A×cos(B)) / R and (A×sin(B)) / R represent the number of pixels moved in the north and east directions, respectively. Core thre It indicates the threshold of the number of pixels that need to be moved (i.e. iteration index).

[0098] Step 4: Based on the data slices of quadtree decomposition, the DEM difference map obtained after rough registration in step (3) is subjected to multi-scale data segmentation using quadtree decomposition technology. In a given data slice (also called a segmentation unit), the average slope information is used as a reference indicator to determine whether to perform data segmentation. This process can be specifically expressed as formula (3) and formula (4).

[0099]

[0100]

[0101] Where I represents the DEM difference map to be processed, Represents the quadtree segmentation operator. S(i, j, k) represents the quadtree structure of I, (i, j) is the coordinate of the upper left corner of the data slice after segmentation, and k represents the segmentation level (k∈N 0 The size of the segmented slice is M×M pixels (M=2 (m-k) ), where 2m and 2n represent the maximum and minimum sizes of a slice during the cutting process. QTseg_or_not is a Boolean data. If the average slope of the terrain within the data slice is If the slope is greater than the slope threshold (Slopethre), further data segmentation is required; otherwise, no further segmentation is required.

[0102] Step 5: Re-segmentation taking into account the local heterogeneity of the data. For the cells in the initial segmentation map obtained in step (4) whose data slice size is larger than the minimum segmentation size, a new indicator (i.e., standardized standard deviation (SSTD)) is further used to evaluate the spatial heterogeneity of the data slices one by one. When SSTD is larger than a set data heterogeneity threshold (Heter thre ), further data segmentation is required, and the data is directly segmented to the smallest slice size; otherwise, no segmentation is performed. This process can be expressed as formula (5) and formula (6).

[0103]

[0104]

[0105] where dh p , dh p_max , dh p_min and N represent the elevation difference, maximum elevation difference, minimum elevation difference and total number of pixels in a data slice (S(i, j, k)) respectively. Reseg_or_not is a Boolean data used to determine whether the data needs to be segmented again.

[0106] Step 6: Remove the local trend term error. Slice the data according to the final segmentation scheme obtained in step (5). Then simulate the linear plane trend error and the error along the height for each data slice one by one. Subtract the fitted error from the original data to remove the local trend term error. The simulated error can be estimated using formula (7).

[0107] dh(x,y)=a2+b2x+c2y+d2n(x,y) (7)

[0108] Among them, a2, b2, c2 are the coefficients of the local spatial trend to be calculated, and d2 is the coefficient related to the height trend error, which can be estimated by iterative least squares method. η(x, y) represents the elevation information of a certain point.

[0109] Step 7: Local data precision registration. For the data slices obtained after the processing of step (6), the method in step (3) is used to register each slice one by one, and a relatively high registration accuracy (i.e., a smaller registration threshold) is set.

[0110] Step 8: Data splicing and result output. The data slices obtained in step (7) are spliced ​​to obtain a complete corrected elevation difference map, and the results are output for subsequent quantitative and qualitative analysis of surface elevation changes.

[0111] The present invention aims to effectively solve the key problem in the current DEM difference technology - the correction of large-scale DEM difference images under complex error modes, while overcoming the problems of the existing correction methods based on manual supervision being time-consuming and inefficient in processing, and provides an automated large-scale DEM difference correction method based on adaptive data segmentation.

[0112] The automated large-scale DEM differential correction method based on adaptive data segmentation described in the above embodiment of the present invention can realize fully automated large-scale and large-batch DEM differential processing by fully considering the error characteristics and local data heterogeneity characteristics in the DEM differential map, and is suitable for correction of all error types. It has strong operability and scalability, and therefore will help to popularize and promote the application of DEM differential technology in various fields.

[0113] When implementing the present invention, three very representative global DEM product data are selected as test data, namely the global 1 arc second SRTM C-band DEM (abbreviated as SRTM-GL1), SRTM X-band DEM (abbreviated as SRTM-X) and Copernicus DEM (corrected version of TanDEM-X 90m DEM); at the same time, experiments are carried out in the central Asian glacier area where the surface elevation changes are more significant. One of the experimental areas is located near the Ili River Valley in the Tianshan Mountains of China (relatively flat terrain, few glaciers), and the other experimental area is located in the central area of ​​the Pamir Plateau in Tajikistan (rugged and complex terrain, many glaciers). The specific implementation steps of the present invention are as follows:

[0114] Step 1: In the Tianshan experimental area, SRTM-GL1 and SRTM-X DEM data are used to estimate the elevation difference in the glacier area (actually representing the difference in penetration depth of X-band and C-band radars in the glacier area). At the same time, SRTM-GL1 DEM and Copernicus DEM are used for differentiation to study the changes in glacier thickness in the Pamir experimental area. Both experimental areas use a 2°×2° geographic grid, so the three DEM products are first spliced ​​separately, and the SRTM-X DEM with the WGS84 ellipsoid as the vertical reference is benchmarked (subtracting the geoid anomaly), and the data resolution is resampled from 25m to 30m. Then, the DEM difference is performed on the two sets of data, and the initial DEM difference map of the non-glacier area (stable area) is obtained by combining the glacier boundary information, as shown in Figure 1. Figure 2 (b) and Figure 3 (b) as shown.

[0115] Step 2: Use formula (1) and formula (2) to perform overall rough correction (including spatial trend removal and rough registration) on the DEM difference maps of the two experimental areas obtained in step (1). In the data registration process, the threshold for stopping iteration (i.e., registration threshold) is set to 0.4 pixels. The results after overall correction are shown in Figure 2. Figure 2 (c) and Figure 3 (c) as shown.

[0116] Step 3: For the DEM differential map obtained in step 2, use formulas (3) and (4) to perform data segmentation based on quadtree decomposition, where the average terrain slope threshold in each data slice is set to 5°, and the maximum and minimum sizes of the slices are set to 1024*1024 and 256*256 pixels respectively. The segmentation of the DEM differential map based on quadtree decomposition is shown in the figure. Figure 2 (g) and Figure 3 (g) as shown.

[0117] Step 4: For the data slices whose size is larger than the minimum size obtained in step 3, further use formula (5) to calculate the standardized standard deviation of the elevation difference in each slice, and use formula (6) to determine whether to perform further data segmentation, where the threshold for the standardized standard deviation is set to 0.01. If cutting is required, it is directly cut to the minimum size. The final data cutting situation after the cutting plan is updated is as follows: Figure 2 (h) and Figure 3 (h) shown.

[0118] Step 5: Based on the final data cutting information obtained in step 4, use formula (7) to remove the local plane spatial trend and correct the error along the height direction for each data slice one by one.

[0119] Step 6: Use formula (2) again to perform local data registration on the data slices obtained in step 5 one by one. Here, the iterative threshold of registration is set to 0.1 pixel.

[0120] Step 7: Splice the data slices obtained in step 6 in the original way to obtain the final corrected DEM difference map, as shown in the figure. Figure 2 (l) and Figure 3 (l) as shown.

[0121] In order to more clearly illustrate the effect of the present invention, we also designed two experiments to perform data segmentation with regular sizes (i.e., 1024 pixels and 256 pixels), and the data segmentation effect diagrams are shown in FIG. Figure 2 (e) Figure 3 (e) and Figure 2 (f) Figure 3 As shown in (f), the final local correction effect images are as follows Figure 2 (i) Figure 3 (i) and Figure 2 (j) Figure 3 (j) as shown.

[0122] Figure 2 The experimental results of differential processing using SRTM-GL1 and SRTM-X DEM data in the Tianshan area are obtained. From the original DEM differential map (the value range is ±20 meters, Figure 2 (b)), the spatial trend error and high-frequency strip error are complex and obvious. After the overall data correction processing ( Figure 2 (c)), the offset in the entire vertical direction is corrected (the average height difference is reduced from -1.08m to 0.72m, and the standard deviation also shows a slight decrease). For the data segmentation process, when considering the heterogeneity of the data, the segmentation strategy is updated ( Figure 2 (g) and Figure 2 (h)). From the final result ( Figure 2 (i)-(p)) It can be seen that the segmentation strategy is effective. At the same time, among all the results, the method proposed by the present invention and the regular segmentation method of small size (256×256 pixels) showed the best results (standard deviations were 4.17 meters and 4.16 meters, respectively), followed by the method based on quadtree segmentation (standard deviation: 4.22 meters), and the regular segmentation of large size (1024×1024 pixels) performed the worst (standard deviation: 4.28 meters). Figure 2 (i) and Figure 2 As can be seen in (k), there are still obvious residual errors in the latter two results. However, it should be noted that, despite having the same performance, the calculation time of the method of the present invention is shorter than that of the small-size regular segmentation method. Finally, compared with the original elevation difference map, the standard deviation of the difference map corrected by the method of the present invention is reduced by 25% compared with the previous one.

[0123] Figure 3 The experimental results of differential processing using SRTM-GL1 DEM and Copernicus DEM data obtained in the central Pamir region are shown. Figure 3 As shown in (b), the original DEM difference map (value range ±60 meters) in this area contains serious gross errors. The overall data correction process only corrects the overall average vertical deviation between the two DEMs, and the standard deviation does not change significantly ( Figure 3 (c) and Figure 3 (d)). In addition, compared with the results based on quadtree segmentation, the results processed by the present invention (i.e., the adaptive data segmentation method) only show a slight update in the local flat area ( Figure 3 (g) and Figure 3 (h)). The final results show that the performance of the large-scale regular data segmentation method is relatively poor (standard deviation: 9.98 meters), and there is still a significant residual error in the corrected results ( Figure 3 (i)), which is similar to the situation in the Tianshan experimental area. However, for the other three methods, due to almost the same segmentation strategy, they show the same performance, that is, all standard deviations are 9.50 meters. This shows that in areas with rugged terrain, the adaptive segmentation method is almost equivalent to the small-size data segmentation method. Finally, relative to the original elevation difference map, the standard deviation of the difference map corrected by the method of the present invention is 10% lower than before. It can be seen from the above two cases that the method of the present invention is very effective for correcting DEM difference maps and can handle different types of errors.

[0124] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the implementation methods of the present invention, and should be understood that the protection scope of the present invention is not limited to such special statements and embodiments. Those skilled in the art can make various other specific variations and combinations that do not deviate from the essence of the present invention based on the technical revelations disclosed by the present invention, and these variations and combinations are still within the protection scope of the present invention.

Claims

1. An automated large-scale DEM differential correction method based on adaptive data segmentation, characterized in that: The following steps are involved: Step 1: Data preprocessing: Unify the coordinate system reference of two or more different DEM data, then perform direct difference processing on the DEM data to obtain a DEM difference map, and then combine the surface type boundary data to mask the potential change area to obtain the initial elevation change map of the stable area; Step 2: Overall data rough correction processing: First, check and correct the large-scale spatial trend errors of the elevation change map in the stable area; then, on this basis, perform initial data rough registration on the elevation difference map; Step 3: Adaptive data segmentation: First, the terrain slope information is used as a reference, and the elevation difference map after overall correction is initially segmented according to the quadtree segmentation technology; then, the slices with strong local heterogeneity are further segmented with the standardized standard deviation as a reference, so as to obtain the final data segmentation solution; Step 4: Local data fine correction processing: According to the obtained data segmentation scheme, the elevation difference map after rough correction is segmented, and the data slices are first spatially detrended and corrected one by one, and then refined data registration is performed; Step 5: Splice the block-corrected data to generate the final corrected DEM difference map.

2. The method for automatic large-scale DEM differential correction based on adaptive data segmentation according to claim 1, characterized in that: The step 1 comprises the following sub-steps: Step 11: Convert the DEM ellipsoidal height with the earth ellipsoid as the vertical reference to the orthometric height data with the geoid as the reference; Step 12: Use appropriate sampling methods to unify the resolution of DEM data; Step 13: Perform differential processing on DEM data of the same size and in the same geographical range to obtain a DEM differential map; Step 14: Mask the area with surface elevation changes to obtain the initial elevation difference map of the stable area; Step 15: Remove outliers from the initial elevation difference map.

3. The automated large-scale DEM differential correction method based on adaptive data segmentation according to claim 1, characterized in that: The step 2 specifically comprises the following sub-steps: Step 21: Use the elevation difference map of the stable area to perform initial overall spatial trend removal; Step 22: Calculate terrain slope and aspect information using one of the DEMs in the differential processing as a reference; Step 23: The elevation difference map after overall detrending is roughly aligned with the overall data using the slope and aspect information.

4. The method for automatic large-scale DEM differential correction based on adaptive data segmentation according to claim 3, characterized in that: In step 2, formula (1) is used to simulate the trend error of spatial variation; dh(x,y)=a1+b1x+c1y+d1xy+e1x 2 +f1y 2 (1) Among them, dh(x, y) represents the elevation difference of the point with coordinates (x, y) in the data; a1, b1, c1, d1, e1 and f1 are the parameters of the overall spatial trend to be estimated; the entire estimation process is calculated using the iterative least squares method.

5. The method for automatic large-scale DEM differential correction based on adaptive data segmentation according to claim 4, characterized in that: The data rough registration described in step 23 is specifically as follows: an analytical three-dimensional registration method based on slope and aspect information is used to perform overall data registration on the obtained DEM difference map, and a relatively rough registration accuracy is set to achieve preliminary alignment of the data. This process can be specifically expressed as formula (2); Among them, slope and aspect represent the terrain slope and slope direction of a point in the data respectively; A and B represent the magnitude and direction of the potential horizontal displacement between DEMs respectively; C represents the overall vertical average deviation between DEMs, which is expressed as In addition, flag is a Boolean data, that is, yes or no, which is used to determine whether iterative processing is required during data processing; R is the spatial resolution, (A×cos(B)) / R and (A×sin(B)) / R represent the number of pixels moving in the north and east directions respectively; Core thre The threshold value of the number of pixels that need to be moved, i.e., the iteration index.

6. The method of automatic large-scale DEM differential correction based on adaptive data segmentation according to claim 1, characterized in that: The step 3 specifically includes the following sub-steps: Step 31: Taking the terrain slope as a reference, set an appropriate slope threshold, and use the quadtree segmentation technology to perform initial data segmentation on the elevation difference map; Step 32: Calculate the standardized standard deviation of the local elevation difference map for the data slices whose size is larger than the minimum segmentation size; Step 34: Using the standardized standard deviation as a reference, set a suitable threshold, segment the data units above the threshold again, and segment them to the minimum segmentation size, thereby obtaining a final data segmentation solution.

7. The method for automatic large-scale DEM differential correction based on adaptive data segmentation according to claim 6, characterized in that: The step 31 specifically includes: using quadtree decomposition technology to perform multi-scale data segmentation on the DEM difference map obtained after rough registration. This process is specifically expressed as formula (3) and formula (4); Where I represents the DEM difference map to be processed, represents the quadtree segmentation operator; S(i, j, k) represents the quadtree structure of I, (i, j) is the coordinate of the upper left corner of the data slice after segmentation, and k represents the segmentation level (k∈N 0 ); the slice size after segmentation is M×M pixels (M=2 (m-k) ), of which 2 m and 2 n Respectively represent the maximum and minimum size of a slice during the cutting process; QTseg_or_not is a Boolean data. If the average slope of the terrain within the data slice is Greater than the slope threshold Slope thre , then further data segmentation is required; otherwise, no further segmentation is required.

8. The method of automatic large-scale DEM differential correction based on adaptive data segmentation according to claim 7, characterized in that: The re-segmentation described in step 34 is specifically as follows: The standardized standard deviation SSTD is used to evaluate the spatial heterogeneity of data slices one by one; when SSTD is greater than a set data heterogeneity threshold Heter thre , then further data segmentation is required, and the data is directly segmented to the smallest slice size; otherwise, no segmentation is performed; This process is expressed as formula (5) and formula (6); where dh p , dh p_max , dh p_min and N represent the elevation difference, maximum elevation difference, minimum elevation difference and total number of pixels within a data slice (S(i, j, k)) respectively; Reseg_or_not is a Boolean data used to determine whether the data needs to be segmented again.

9. The method of automatic large-scale DEM differential correction based on adaptive data segmentation according to claim 1, characterized in that: The step 4 specifically includes the following sub-steps: Step 41: According to the final data segmentation scheme, perform spatial trend correction again for each segmentation unit; Step 42: According to the final data segmentation scheme, perform data precise registration again for each segmentation unit; Step 43: The segmented and re-corrected data are spliced ​​according to the original positions, and the final corrected elevation difference map is output.

10. The method for automatic large-scale DEM differential correction based on adaptive data segmentation according to claim 9, characterized in that: Step 41 is specifically as follows: the data is segmented and sliced ​​according to the final segmentation scheme, and then the linear plane trend error and the error along the height direction are simulated for each data slice one by one, and the fitting error is subtracted from the original data to complete the removal of the local trend term error; the simulated error is estimated using formula (7); dh(x,y)=a2+b2x+c2y+d2n(x,y) (7) Among them, a2, b2, c2 are the coefficients of the local spatial trend to be calculated, d2 is the coefficient related to the height trend error, these coefficients can be estimated by iterative least squares method; η(x, y) represents the elevation information of a certain point.

Images