An automatic editing method for converting InSAR-DSM to DEM

Abstract

The present invention provides an automatic editing and acquisition method for InSAR-DSM to DEM. The method utilizes SAR intensity images, SAR coherence information, satellite imaging information and reliable public external data, and aims to fully automatically complete the editing process of InSAR-DSM to DEM. The method can effectively and quickly complete the data anomaly identification and repair of InSAR-DSM, and on this basis automatically identify water systems and built-up areas, and perform adaptive elevation editing processing on water areas, built-up areas and vegetation areas respectively. Under the premise of ensuring the integrity of InSAR-DSM data, the method realizes fully automatic elevation editing of water bodies and built-up areas respectively. The whole process structure is clear, and has the advantages of clear process, simple implementation, large processing range, high degree of automation, etc.

Description

Technical Field

[0001] The present invention relates to the technical field of producing DEM based on InSAR technology, and in particular to an automatic editing and acquisition method for converting InSAR-DSM to DEM. Background Art

[0002] Since the 1980s, with the advancement of earth science and global digitization and informatization, exploring global change and sustainable development has become a key issue. "Digital Earth" can be understood as the informatization of the Earth, using a set of numbers to describe the spatial characteristics of every corner of the planet. Since the vast majority of information humans receive in their daily lives is related to spatial information, the production of global digital terrain products is becoming increasingly significant. Commonly used digital terrain products are digital elevation models and digital surface models.

[0003] A Digital Elevation Model (DEM) can be defined as a regular two-dimensional array that describes the varying elevations of exposed ground, representing the true relief of the terrain. It is typically stored as raster data and serves as one of the most important foundational terrain data, serving many scientific research and engineering applications, including ecology, hydrology, geomorphology, and other disciplines. It plays a fundamental role in understanding the Earth's surface, expressing landform morphological characteristics, and revealing geological mechanisms. It is an important geospatial data model widely used in Earth science analysis and modeling. A Digital Surface Model (DSM) is also a regular two-dimensional array used to describe surface elevation values. However, unlike a DEM, a DSM describes the elevation of all features on the Earth's visible surface.

[0004] DSM and DEM share the same production methods, data storage methods, and presentation formats. However, they differ in the objects represented, manually edited content, and application directions. DSM contains many elevations of artificial buildings, vegetation, and other objects. Common acquisition methods include field surveying, photogrammetry, laser scanning, and Synthetic Aperture Radar (InSAR). InSAR technology, as an active microwave remote sensing tool, can observe the Earth's surface around the clock, regardless of weather conditions. It offers advantages such as high precision and a wide range. Compared to traditional geodetic techniques, it can measure large areas at a fraction of the cost. For the production of a global DEM, InSAR technology is currently one of the best methods, considering both time and cost.

[0005] It's worth noting that what's obtained using InSAR technology is not a DEM but a DSM. However, DEM is more widely used than DSM in the "Digital Earth" strategy, and exploring how to obtain DEM from DSM has become a hot topic in this strategy. First, InSAR-DSM has many DSM anomalies and data missing areas caused by poor coherence due to geometric distortion. There's also a lot of noise in water areas, making it not a true DEM. In particular, the DSM in urban and forested areas differs significantly from the DEM, as the DSM includes many artificial buildings, artificial and natural vegetation, and so on. Therefore, how to repair the InSAR-DSM and edit the elevation of water areas, built-up areas, and vegetation areas requires special attention. Summary of the Invention

[0006] The purpose of the present invention is to provide an automatic editing and acquisition method for InSAR-DSM to DEM in response to the deficiencies in the above-mentioned background technology, so as to repair the elevation anomalies of InSAR-DSM and complete the elevation editing of the area covered by the ground object.

[0007] In order to achieve the above object, the present invention provides an automatic editing and acquisition method for InSAR-DSM to DEM, comprising the following steps:

[0008] S1, combines prior DEM information, SAR image intensity information, InSAR interferometric coherence coefficient information and satellite imaging parameters to accurately identify geometric distortion areas; identify InSAR-DSM holes and low-quality areas, and together with the geometric distortion areas, form a DSM data anomaly distribution mask;

[0009] S2, based on the anomaly distribution mask of the DSM data, the InSAR-DSM hole area is obtained, and the hole is filled by combining the public terrain product and the geospatial interpolation method to obtain the complete InSAR-DSM;

[0010] S3, based on the SAR coherence map and SAR intensity image, combined with external DEM information, the single threshold segmentation method is used to extract the SAR feature representation of the water body and obtain the water body mask;

[0011] S4, based on the SAR coherence map and SAR intensity image, the SAR feature representation of the building is extracted based on the image texture features and machine learning to obtain the building area mask;

[0012] S5, locates the water body position according to the water body mask, and classifies the water body into three categories according to its size and shape: rivers, lakes, and oceans. Based on the different characteristics of different water bodies in terms of elevation changes, automatic elevation editing is performed on each of them;

[0013] S6, locates the building position according to the built-up area mask, divides the ground points and non-ground points, and performs automatic elevation editing through interpolation based on the distribution characteristics of ground points in built-up areas of different terrains and sizes;

[0014] S7, locate the vegetation boundary position according to the vegetation mask, and screen the control points within the image range, calculate the vegetation height difference correction surface based on the inverse distance weighted method of the control points, and perform automatic height editing processing;

[0015] S8, based on the potential edge differences, an adaptive smooth transition method is adopted to obtain a DEM with reasonable transition information.

[0016] Furthermore, S1 specifically includes:

[0017] Obtain satellite incident angle and azimuth information from SAR satellite imaging parameter files;

[0018] Reconstruct the resolution of the public terrain products obtained and make the SAR image resolution within the study area the same by interpolation;

[0019] Calculate the slope, aspect, and elevation of the InSAR-DSM pixel by pixel, and determine whether it is a shadow or overlap area based on its relationship with the satellite's incidence angle and azimuth parameters;

[0020] Using SAR imaging parameters, we can derive a terrain product data quality indicator map to identify low-quality areas.

[0021] The missing areas of InSAR-DSM data are identified, and together with the geometric distortion areas and low-quality areas, they constitute the data anomaly areas, forming the DSM data anomaly distribution mask.

[0022] Furthermore, S2 specifically includes:

[0023] Holes formed in areas indicated by abnormal distribution masks of DSM data are eliminated. Based on the different correlation strengths of surrounding terrain data under different terrain complexity and the different terrain correlation strengths caused by holes of different sizes, the holes to be filled are adaptively divided into low-difficulty areas and high-difficulty areas.

[0024] The difficulty of filling holes is determined based on the complexity of the local terrain. For low-difficulty areas, geospatial interpolation is used. For high-difficulty areas, public external terrain products are introduced and the corresponding regional data is selected as an external reference to assist in the automatic filling of holes.

[0025] Furthermore, S3 specifically includes:

[0026] SAR intensity images and SAR coherence maps were used to identify water bodies in the study area;

[0027] Based on the SAR intensity image, the grayscale histogram distribution information of the intensity image is statistically analyzed, and the water system is identified according to the bimodal method. The lowest point of the image valley between the water body and the land is found as the critical value, and the image is segmented based on this threshold.

[0028] The SAR intensity image, SAR coherence map and slope map generated based on external DEM information are combined by band combination to obtain a pseudo-color image of the study area. This image contains the above three types of information, and water bodies are extracted from the pseudo-color image.

[0029] Furthermore, S4 specifically includes:

[0030] Perform filtering preprocessing on the SAR intensity image, extract the building area by counting the image texture features and determining the extraction threshold;

[0031] The texture features of SAR intensity images are extracted using gray-level co-occurrence matrix. The appropriate feature expression operators are screened through principal component analysis and Bhattacharyya distance evaluation, and the feature expression operators with strong feature expression capabilities are integrated.

[0032] Through machine learning classification, the fused features are classified, and building extraction is completed with the assistance of external DEM information.

[0033] Furthermore, S5 specifically includes:

[0034] Using the extracted water mask files and global high-resolution coastline vector data, the oceans are classified according to their geographical locations; the remaining water bodies are classified into lakes and rivers according to their shapes and areas.

[0035] Based on the water body classification results, elevation editing is performed on the three types of water bodies: oceans, lakes and rivers, in view of the different characteristics of different types of water bodies and the characteristics of actual data.

[0036] Furthermore, S6 specifically includes:

[0037] The extracted built-up area mask files are morphologically preprocessed, and the built-up areas are divided into large-area built-up areas and small-area built-up areas based on the degree of terrain relief and area. The building heights of the large-area built-up areas and the small-area built-up areas are adaptively filtered out respectively.

[0038] Based on the echo signal and geometric characteristic factors, a joint constraint condition is constructed to distinguish ground points from non-ground points. InSAR ground points are extracted by combining error statistics and point cloud cleaning.

[0039] Based on the selected InSAR ground points, the elevation recovery process is achieved by using an adaptive interpolation method based on the different distribution characteristics of ground points in urban areas of different areas and the original terrain characteristics.

[0040] Based on the texture features of the original DSM, texture restoration is performed after elevation restoration.

[0041] Furthermore, S7 specifically includes:

[0042] The obtained vegetation mask product was reconstructed and interpolated to the same resolution as the SAR image in the study area, and the scattered small forests were removed and the bare land in the large forests was filled.

[0043] The generated vegetation mask was used as the forest boundary control, and the ICESat-2ATL08 control points within the image range were screened;

[0044] The vegetation elevation correction surface is calculated based on the inverse distance weighted method of ATL08 control points, and the vegetation elevation is edited.

[0045] The above solution of the present invention has the following beneficial effects:

[0046] The automatic editing and acquisition method from InSAR-DSM to DEM provided by the present invention fully utilizes SAR intensity images, SAR coherence information, satellite imaging information, and reliable public external data, and aims to fully automatically complete the editing process from InSAR-DSM to DEM. It can effectively and quickly complete the identification and repair of InSAR-DSM data anomalies, and on this basis, fully automatically identify the distribution of water bodies and built-up areas. While ensuring the integrity of InSAR-DSM data, it can realize fully automatic elevation editing of water bodies and built-up areas respectively. The entire process structure is clear, and has the advantages of clear process, simple implementation, large processing range, and high degree of automation.

[0047] Other beneficial effects of the present invention will be described in detail in the subsequent specific implementation section. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 is a flow chart of the steps of the present invention;

[0049] Figure 2 Schematic diagram of geometric distortion of the present invention, wherein (a) is shadow and (b) is masking;

[0050] Figure 3 The data used in the specific embodiments of the present invention include (a) the original DSM data of the experimental area, (b) the intensity image of the experimental area, and (c) the coherence map of the experimental area.

[0051] Figure 4This is a demonstration of the geometric distortion area recognition effect in a specific embodiment of the present invention;

[0052] Figure 5 This is a demonstration of the effect of automatic panoramic hole filling in a specific embodiment of the present invention, where (a) is before hole filling and (b) is after hole filling;

[0053] Figure 6 The water extraction effect of a specific embodiment of the present invention is demonstrated, where (a) is the standard water mask result published by DLR, and (b) is the water extraction result of the present invention;

[0054] Figure 7 This figure shows the building extraction results in a specific embodiment of the present invention, where (a) is the standard building mask result published by DLR, and (b) is the building extraction result of the present invention;

[0055] Figure 8 This is a display of the water body classification results in a specific embodiment of the present invention;

[0056] Figure 9 This is a demonstration of the water body elevation editing effect in a specific embodiment of the present invention, where (a) is before the water body elevation editing, and (b) is after the water body elevation editing;

[0057] Figure 10 This is a display of the building area elevation editing effect in a specific embodiment of the present invention, where (a) is the location of the ROI, (b) is the building area elevation after editing, and (c) is the building area elevation before editing;

[0058] Figure 11 This is a display of vegetation editing and filtering elevation in a specific embodiment of the present invention;

[0059] Figure 12 This is a demonstration of the overall effect of DEM editing in a specific embodiment of the present invention, where (a) is before automatic editing and (b) is after automatic editing. DETAILED DESCRIPTION

[0060] The following describes the embodiments of the present disclosure through specific examples, and those skilled in the art can easily understand other advantages and effects of the present disclosure from the contents disclosed in this specification. Obviously, the described embodiments are only a part of the embodiments of the present disclosure, rather than all of the embodiments. The present disclosure can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present disclosure. It should be noted that, in the absence of conflict, the following embodiments and features in the embodiments can be combined with each other. Based on the embodiments in the present disclosure, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present disclosure.

[0061] It should be noted that various aspects of the embodiments within the scope of the appended claims are described below. It should be apparent that the aspects described herein can be embodied in a wide variety of forms, and any specific structure and / or function described herein is merely illustrative. Based on this disclosure, it should be understood by those skilled in the art that an aspect described herein can be implemented independently of any other aspect, and two or more of these aspects can be combined in various ways. For example, any number of aspects described herein can be used to implement an apparatus and / or practice a method. In addition, other structures and / or functionalities other than one or more of the aspects described herein can be used to implement this apparatus and / or practice this method.

[0062] It should also be noted that the diagrams provided in the following embodiments are merely schematic illustrations of the basic concepts of the present disclosure. The diagrams only show components relevant to the present disclosure and are not drawn according to the number, shape, and size of components in actual implementation. In actual implementation, the configuration, quantity, and proportion of each component may be varied at will, and the component layout may be more complex. Furthermore, in the following description, specific details are provided to facilitate a thorough understanding of the examples. However, those skilled in the art will appreciate that the described aspects may be practiced without these specific details.

[0063] like Figure 1 As shown, an embodiment of the present invention provides an automatic editing and acquisition method for converting InSAR-DSM to DEM, comprising:

[0064] S1 combines prior DEM information, SAR image intensity information, InSAR interferometry coherence coefficient information, and satellite imaging parameters to accurately identify geometric distortion areas; identify InSAR-DSM holes and low-quality areas, and together with the geometric distortion areas, form a DSM data anomaly distribution mask.

[0065] In this embodiment, this step obtains satellite incident angle and azimuth information from the SAR satellite imaging parameter file; reconstructs the resolution of the obtained public terrain product and selects an appropriate interpolation method to make it the same as the SAR image resolution in the study area; calculates the slope, aspect, and elevation of the InSAR-DSM pixel by pixel, and determines whether it is a shadow or overlap area based on the relationship between it and the satellite incident angle and azimuth parameters; uses SAR imaging parameters to derive a terrain product data quality indicator map to determine low-quality areas; identifies InSAR-DSM data missing areas, and together with geometric distortion areas and low-quality areas, forms data anomaly areas to form a DSM data anomaly distribution mask.

[0066] It is understandable that geometric distortion means that the pixel information is unreliable, and the elevation values calculated from it are also unreliable. In actual processing, a threshold is often set before phase unwrapping to mask areas with too low coherence to prevent areas with too low coherence from causing distortion in the unwrapping of the entire image. This will cause the DSM generated by InSAR to contain many "holes" with missing data information. In order to accurately identify low-quality areas, this work can be carried out from coarse to fine, first using statistical analysis methods for coarse identification, and then using spatial analysis methods on this basis for fine determination. This two-level processing ensures reliability. The geometric distortion area, hole area, and low-quality area together constitute the data anomaly distribution mask.

[0067] Figure 2 The figure shows the reflection of SAR echo signal in the face of very steep terrain. In this case, Figure 2 In (a), the terrain EFG segment cannot receive the signal sent by the radar, resulting in point EG not being imaged. This is the reason for the existence of the shadow area. The slope EF is defined as the “active shadow”, while the part of FG that is located on the slope EF and cannot reach the radar beam because it is blocked by the hillside EF is called the “passive shadow”. Figure 2 As can be seen in (b), the imaging order of slope BC is reversed and superimposed on the images of slices AB and CD. This is the overlap region, where slope BC is also called "active overlap", while the AB segment and CD segment at the two ends of slope BC are "near overlap" and "far overlap", respectively.

[0068] The judgment of shadow and overlapping areas also needs to be determined based on the satellite flight mode:

[0069] For satellites in orbit raising mode:

[0070] Judgment 1:

[0071]

[0072] Where aspect(i,j) is the aspect angle at pixel (i,j), measured clockwise from true north; slope(i,j) is the slope angle at pixel (i,j); azi_angle(i,j) is the SAR satellite's flight direction at the pixel, measured clockwise from true north; and inc_angle(i,j) is the SAR satellite's angle of incidence at the pixel. It should be noted that the angle of incidence and azimuth are identical for each pixel in an image, as specified in the satellite imaging parameter file.

[0073] For this judgment, if the mutual relationships satisfy the above three equations at the same time, it can be considered that the pixel point belongs to the distorted area.

[0074] Judgment 2:

[0075]

[0076] slope(i,j)>inc_angle(i,j) (3)

[0077] For the second judgment, if the relationship between its parameters satisfies any one or two of the equations (6) and (7), the pixel point can be considered to belong to the distorted area.

[0078] For satellites in descending orbit mode:

[0079] Judgment 1:

[0080]

[0081] For this judgment, if the above three equations are satisfied at the same time, the pixel point can be considered to belong to the distorted area.

[0082] Judgment 2:

[0083]

[0084] slope(i,j)>90-inc_angle(i,j) (6)

[0085] For the second judgment, if the relationship between its parameters satisfies any one or two of the equations (5) and (6), the pixel point can be considered to belong to the distorted area.

[0086] S2, based on the anomaly distribution mask of the DSM data, the InSAR-DSM hole area is obtained, and the hole is filled by combining the public terrain product and the geospatial interpolation method to obtain a complete InSAR-DSM with high reliability.

[0087] In this embodiment, this step first eliminates holes formed in areas indicated by abnormal distribution masks of DSM data. Based on the different correlation strengths of surrounding terrain data under different terrain complexity and the different terrain correlation strengths caused by holes of different sizes, the holes to be filled are adaptively divided into low-difficulty areas and high-difficulty areas. The difficulty of hole filling is determined according to the complexity of local terrain, and geospatial interpolation processing is adopted for low-difficulty areas. For high-difficulty areas, public external terrain products are introduced, and strictly corresponding regional data are selected as external references to assist in completing reliable automatic hole filling.

[0088] It's important to note that data holes in terrain products often arise from low coherence or geometric distortion. Holes often exhibit distinct characteristics depending on imaging and terrain conditions. First, when holes occur in rugged terrain, elevation variations are significant even within a small area, making it difficult to recover missing elevation values using existing surrounding elevation data. Conversely, in relatively flat areas, elevation variations are consistent over a larger area, and using surrounding elevation data from farther away can recover missing elevations. Second, when a contiguous hole is large, the influence of surrounding elevation data decreases toward the center of the hole, often becoming virtually irrelevant. However, when a hole is small, surrounding elevations provide a significant reference for the elevation within the hole. These two influencing factors should not be considered in isolation when filling holes, but rather in combination. Furthermore, the hole-filling method should adaptively partition processing areas into different difficulty levels to facilitate subsequent hierarchical processing. For areas with low difficulty, consider using geospatial interpolation to recover the hole elevation based on existing surrounding elevation data. Faced with the highly difficult areas that have been demarcated, it is difficult to accurately restore the void elevation by relying solely on the existing surrounding elevation data. Therefore, under the premise that external open terrain data is available free of charge worldwide, local terrain reference information can be introduced to assist in filling the gaps.

[0089] In this embodiment, a terrain-related factor that can represent the complexity of the terrain is calculated based on the terrain data. Here, the elevation gradient is calculated:

[0090]

[0091] Table 1 Schematic diagram of adjacent elevations

[0092] <![CDATA[h1]]> <![CDATA[h2]]> <![CDATA[h3]]> <![CDATA[h4]]> <![CDATA[h5]]> <![CDATA[h6]]> <![CDATA[h7]]> <![CDATA[h8]]> <![CDATA[h9]]>

[0093] Among them, h i The relationship is shown in Table 1, where p and q are the horizontal and vertical gradients, w is the ground resolution, and r and s are the intermediate values of the second-order gradient. The terrain gradient can be further used to calculate the surface curvature C. m , which indicates the terrain complexity of the DEM:

[0094]

[0095]

[0096] Secondly, the effectiveness of the interpolation algorithm will change depending on the size of the hole. Specifically, the effect A ij is inversely proportional to the distance:

[0097]

[0098] Where α is the scaling factor and (m,n) are the coordinates of the known point relative to the unknown point (i,j).

[0099] Combining the above two factors, we can get an adaptive filling difficulty division method:

[0100]

[0101] Where λ is the weighting factor, tre is the determined threshold, below which the area is divided into the low-difficulty area C1, and above which the area is divided into the high-difficulty area C2.

[0102] When performing automatic filling of geographic spatial interpolation in low-difficulty areas, first, the elevation points around the holes are selected as preliminary candidate points.

[0103] Secondly, from the candidate points, reference points are evenly distributed according to the direction (as shown in Table 1) and used as reference calculations:

[0104]

[0105] And ensure the number of selected reference points P r Greater than the lowest reference point limit P min , because too few reference points will seriously affect the stability of the results, in this case the search range R should be expanded.

[0106] Then, under the premise of determining the reference point, Kriging interpolation is used to fill the gaps. The value of the unknown location is expressed as the weighted average of the values of the nearest known data points. The weight changes with distance and direction. The specific formula is as follows:

[0107]

[0108] in, is the estimated value of the point (x0, y0), that is, h0 = (x0, y0); λ i is a weight coefficient that needs to be obtained by fitting a kriging model, rather than simply a function related to distance. The elevation of holes in the area can be calculated using this method.

[0109] For challenging areas, automated external data-assisted hole filling is used. Specifically, within the identified difficult hole filling areas, each connected area is processed block by block. First, an external reference terrain product is selected based on the geographic coordinate range. After unifying the resolution and coordinate system and performing data registration, the external data is aligned with the original InSAR terrain data. Second, the local baseline offset between the external data and the original data is calculated, and the corresponding baseline offset is applied to the external data to bring it into line with the baseline of the original data.

[0110] S3, based on the SAR coherence map and SAR intensity image, combined with external DEM information, uses the single threshold segmentation method to extract the SAR feature representation of the water body and obtain the water body mask, providing a data basis for the subsequent use of different editing methods for the water body.

[0111] In this embodiment, this step first uses SAR intensity images and SAR coherence maps to identify water bodies in the study area, providing data support for the subsequent special water area editing module; based on the SAR intensity image, the grayscale histogram distribution information of the intensity image is statistically analyzed, and then the water system is identified according to the bimodal method, and the lowest point of the image valley existing in the middle is found as the critical value between the water body and the land, and the image is segmented based on this threshold; the SAR intensity image, the SAR coherence map and the slope map generated according to the external DEM information are band combined to obtain a pseudo-color image of the study area, which contains the above three types of information, and the pseudo-color image is used for water extraction.

[0112] Water bodies, primarily due to their low scattering rate for satellite microwaves, produce weak echo signals, appearing as dark areas on SAR intensity images, making them easily identifiable. The single threshold segmentation method primarily uses the bimodal method for threshold identification and water extraction. Histogram analysis of SAR intensity images reveals the presence of a minimum point between two maximum points in the histogram curve. The bimodal approach to single threshold segmentation involves finding the lowest point in the image valley between two maxima, which serves as the critical value separating water from land. This threshold is then used to segment the image, enabling water extraction. The combined water extraction scheme, which combines SAR intensity images, SAR coherence maps, and external DEM information, addresses the problem of weak echo signals from water systems in SAR satellite imaging due to specular reflection. This results in water systems appearing darker in SAR intensity images, making them easier to distinguish from other land features. Furthermore, water systems exhibit significantly lower coherence compared to other land features on SAR coherence maps, providing the basis for combined water extraction using SAR coherence maps and SAR intensity images. The SAR intensity image, SAR coherence map and slope map generated based on the external reference DEM are combined by band combination to obtain a pseudo-color image of the study area. The image contains the above three types of information, and water bodies are extracted from the pseudo-color image.

[0113] S4, based on the SAR coherence map and SAR intensity image, image texture features and machine learning, extracts the SAR feature representation of the building and obtains the building area mask, providing a data basis for the subsequent use of different editing methods for the building.

[0114] In this embodiment, this step first performs filtering preprocessing on the SAR intensity image, and extracts the building area by counting the image texture features and determining the extraction threshold; uses the gray level co-occurrence matrix to extract the texture features of the SAR intensity image, and screens out suitable feature expression operators through principal component analysis and Bhattacharyya distance evaluation, and fuses feature expression operators with strong feature expression capabilities; uses machine learning classification to classify the fused features, and uses external DEM information to assist in building extraction.

[0115] It should be noted that statistical analysis shows that buildings exhibit distinct texture characteristics in SAR intensity images compared to other landforms. Because built-up areas are typically interspersed with man-made structures and roads, one of the prominent texture features of these areas in SAR intensity images is alternating light and dark patterns. Man-made structures, due to their regular structure and building materials, appear as bright lines in SAR intensity images, while roads, due to their two-dimensional shape, structure, and building materials, appear as dark lines. The combined distribution of these two features results in alternating light and dark lines in SAR intensity images. In this embodiment, based on this characteristic, a gray-level co-occurrence matrix is used to calculate the grayscale distribution within a certain window size, and an appropriate step size is selected for traversal to achieve grayscale information statistics for the entire image. For different feature expression operators, the Bhattacharyya distance is used as a measure of feature representation capability, and principal component analysis is applied to remove inter-feature correlations, effectively fusing the features. Based on this, K-means clustering machine learning is used to classify the fused image. Using unsupervised learning, the final classification result is obtained through continuous iterative learning until the inter-class difference is maximized. Building areas are extracted by exploiting their high brightness on SAR intensity images. Incorporating SAR coherence maps and external DEM information can effectively improve building extraction accuracy. Given the high coherence of buildings, misidentified areas can be eliminated by removing coherence below a certain threshold. Furthermore, buildings are rarely located in areas with steep slopes. Slope calculation and slope filtering based on external DEM information can effectively improve building extraction accuracy.

[0116] It should be noted that the Gray Level Co-occurrence Matrix (GLCM) counts the joint probability distribution of the positions of two pixels. Let S be the set of pixel pairs with specific spatial connections in the target area R, then the co-occurrence matrix P can be defined as:

[0117] P(i,j)=[p(i,j,d,θ)] (14)

[0118] This formula can be described as: the frequency of simultaneous occurrence of pixel pairs with grayscale values i and j (i, j = 0, 1, 2, 3, ..., N) in the image, satisfying the conditions of step size d, direction θ. θ is the angle of the line connecting the two pixels, and texture features are usually obtained along four angles: 0°, 45°, 90°, and 135°.

[0119] Texture features are further extracted from it, and variance, difference, and uniformity are used as features. The calculation formulas are as follows:

[0120]

[0121]

[0122]

[0123] S5 locates the water body according to the water body mask, and further divides the water body into three categories: rivers, lakes, and oceans according to the size and shape of the water body. Considering the different characteristics of different water bodies in terms of elevation changes, automatic elevation editing is performed on each of them.

[0124] In this embodiment, this step first uses the extracted water body mask file and global high-resolution coastline vector data to classify the oceans in the water bodies according to their geographical locations, and then classifies the remaining water bodies into lakes and rivers according to their shapes and areas; based on the water body classification results, the three types of water bodies, oceans, lakes and rivers, are subjected to elevation editing processing according to the different characteristics exhibited by different types of water bodies and the characteristics of actual data.

[0125] Understandably, oceans and lakes are similar in size and area, making them difficult to distinguish. Therefore, external data is needed to assist in classification. While currently available water products include water masks for lakes and rivers, these masks are often not categorized into multiple water bodies, have a long data update cycle, or only exist for specific regions. Therefore, they are not suitable for global DEM production. However, global coastline data products do exist. While coastlines often change over time, these changes are negligible relative to the overall ocean area. Therefore, when classifying oceans, high-resolution global coastline vector data is used to first classify oceans. The remaining water masks, excluding oceans, need to be classified into lakes and rivers. Lakes are stationary water bodies, while rivers are flowing water bodies. In theory, lake and river classification can be done based on the magnitude of elevation changes. However, due to low backscatter caused by specular reflection and temporal decoherence, water bodies appear as random noise or elevation anomalies in the raw InSAR DSM, failing to reflect the true water elevation. This makes it difficult to distinguish lakes from rivers based on elevation. Therefore, in actual classification, lakes and rivers are primarily classified based on their morphology and area. Accordingly, this solution sets appropriate classification indicators. Because lakes are generally circular and rivers are generally elongated, the ratio of water area to maximum extension can be used to classify lakes and rivers. Based on a large amount of actual data, an appropriate threshold is then determined to classify the remaining water bodies, excluding the ocean, into lakes and rivers.

[0126] Specifically, in this embodiment, a rectangular vector is first determined based on the DSM map range. This rectangular vector is then intersected with the global coastline vector data to obtain a planar vector of the coastline within the map range. This planar vector is then converted into a raster file consistent with the DSM map resolution. This results in an ocean mask raster file consistent with the DSM map range and resolution. The oceans within the water body are then classified based on the ocean mask file. Based on an index t, which represents the ratio of the classified water area to the maximum extension of the water body, connected domains with an index greater than a threshold α (i.e., t>α) are classified as lakes; otherwise, they are classified as rivers.

[0127] In this embodiment, the EGM2008 (Earth Gravitational Model 2008) earth gravity field model is used to calculate the geoid height and edit the ocean elevation. The geoid height is also called the geoid difference, which refers to the distance between a point on the geoid and the corresponding position on the reference ellipsoid projected along the normal. The geoid refers to the gravity equipotential surface in the earth's gravity field that coincides with or is closest to the mean sea level in a free and static state. The geoid calculated using gravity data is called the gravity geoid. The gravity geoid height is often calculated using the earth's gravity field model.

[0128] The height of the geoid is calculated according to the following formula:

[0129]

[0130] Among them, N GM is the height of the model geoid, r, θ, and λ are the distance from the center of the Earth, colatitude, and longitude, respectively. GM is the product of the gravitational constant and the mass of the Earth, and a is the semi-major axis of the reference ellipsoid. and is the fully normalized coefficient, To fully normalize the associated Legendre function, N max The highest order of expansion of the Earth's gravity field model.

[0131] The fully normalized associated Legendre function in the above formula The cross-order recursive method is used for calculation, and the principle is as follows:

[0132] Where n is the order and m is the degree, and its initial values are as follows:

[0133]

[0134] When m=0,1, the basic recursive formula is used for calculation, namely the following formula:

[0135]

[0136] In the above formula,

[0137]

[0138]

[0139] When m≥2, the cross-order recursive formula is used, which is the following formula:

[0140]

[0141] In the above formula:

[0142]

[0143]

[0144]

[0145]

[0146] When m≥2, The value of can be obtained by linear combination, and the recurrence coefficient α nm , β nm , γ nm The values of are all less than 1, so the cross-order recursion is stable and reliable.

[0147] In this embodiment, the elevation values of river bank pixels are obtained and denoised based on the principle of connected domain expansion. Kriging interpolation is then used to interpolate the river pixel elevations based on the denoised river bank pixel elevation values. After kriging interpolation of the river pixel elevation values, the interpolated elevation values are subjected to multiple local smoothing filters until the river elevation changes smoothly.

[0148] In this embodiment, according to the principle of connected domain expansion, the elevation values of the pixels on the lake shore are obtained, and their distribution histogram is drawn. The elevation value at one-quarter of the peak drop (towards a lower elevation value) in the histogram is selected as the unified elevation of the lake.

[0149] S6 locates the building position according to the built-up area mask, divides the ground points and non-ground points according to rules, assumptions, error theory and constraints, considers the distribution characteristics of ground points in built-up areas of different terrains and sizes, and uses multi-class interpolation method for automatic elevation editing.

[0150] In this embodiment, this step first performs morphological preprocessing on the extracted building area mask file, and divides the building area into large-area building area and small-area building area based on the terrain undulation and area, and adaptively filters out the building heights of the large-area building area and the small-area building area respectively.

[0151] It should be noted that different ground objects have different echo characteristics due to factors such as material, shape, and surface roughness. This is an important basis for radar remote sensing technology to distinguish ground object types. Considering the simple types of ground objects in built-up areas, they can be distinguished based on clustering. First, a nonlinear gain amplitude equalization method is used to enhance the strength of the target signal and reduce the impact of background clutter and noise, thereby improving target detection performance. Specifically, this embodiment uses nonlinear gain equalization to perform a nonlinear transformation on the signal, applying different gains to the target signal and noise signal, thereby improving the signal-to-noise ratio and target detection performance. Then, the K-means clustering method is used to classify the signal strength in the built-up area. The classification results will serve as the first constraint condition for distinguishing ground points from non-ground points. The geometric feature descriptor is constructed using the curve, surface, and neighborhood features of the DSM, and geometric feature factors for points and domains are proposed. The second constraint condition is formed by constraining the geometric factors. By combining these two constraints and setting an appropriate constraint strength, ground points can be selected. While the ground points selected based on the constraint equations satisfy the interpolation requirements, efficiency is an unavoidable consideration in actual production and application. Therefore, this embodiment proposes a uniform cleaning solution that downsamples the ground points to improve interpolation efficiency without affecting the interpolation effect. Considering that the same interpolation method may produce inconsistent results in scenes of varying complexity, after obtaining evenly distributed ground points, an adaptive interpolation method is selected based on their distribution and number to complete elevation recovery.

[0152] As described above, in this embodiment, the extracted building area intensity image is classified using the K-means algorithm.

[0153] The criterion for clustering is to minimize the sum of squared distances of each cluster through continuous iteration, that is:

[0154]

[0155] A certain range (circular range or rectangular range) is selected as the neighborhood range of point P, and the feature information of the geometric descriptor within the neighborhood is calculated.

[0156] For a surface, the normal vector of a point is the vector perpendicular to the tangent plane at that point. The specific calculation is as follows:

[0157] Assume that the surface equation is F(x,y,z)=0, and P is any point on the surface that satisfies F(x p ,y p ,z p )=0, if we want to find the normal vector of point P, we only need to find its partial derivatives in three directions:

[0158]

[0159] Then convert the normal vector n into a unit normal vector:

[0160]

[0161] Count the Z-axis component e of the unit normal vector within the neighborhood of point P z The histogram distribution of H Pe .

[0162] For any tangent line at point P, let the curve equation be y = f(x), and have a second-order derivative, we can get the curvature calculation formula:

[0163]

[0164] For a curved surface, the calculation of curvature is relatively complicated. In actual calculation, in order to simplify the calculation and improve efficiency, the directions of the principal curvature radius are defined as the x-axis direction and the y-axis direction, which are called the improved principal curvature radius directions. The formula for calculating the average curvature is improved to the form of formula (32), and the improved average curvature at point P is obtained as:

[0165]

[0166] The improved mean curvature histogram distribution within the neighborhood of point P is obtained by counting the H Pκ .

[0167] The geometric feature factors of point P are constructed based on the descriptor as follows:

[0168] M(P)={[e p κ p ]H Pe H Pκ H Pz} (33)

[0169] The domain geometric characteristic factors of the entire neighborhood are established as follows:

[0170] M all ={H pe H Pκ H Pz} (34)

[0171] According to the empirical prior information of urban ground points and non-ground points, the geometric characteristic factors M(G) of ground points and M(N) of non-ground points can be obtained; as well as the geometric characteristic factors M(N) of the ground area. all (G) and the geometric characteristic factor M of the non-ground area all (N).

[0172] Combining the above two constraints, the constraint equation is established as follows:

[0173]

[0174] dis(M all ,M all (G)) <tre all (36)

[0175] Where dis(·) represents the difference between geometric feature factors; tre represents the difference threshold of point geometric feature factors; tre all represents the difference threshold of domain geometric characteristic factors; c g Indicates the category to which the ground point belongs. Only points that satisfy both equations (35) and (36) can be selected as ground points.

[0176] It should be noted that commonly used methods for interpolating discrete ground points to generate a DEM include: inverse distance weighted method, Kriging interpolation method, minimum curvature method, linear interpolation triangulation method, moving average method, local polynomial method, etc. Different interpolation methods can be implemented based on different interpolation weighting methods. The same interpolation method will produce inconsistent results when used in scenes of different complexity. Therefore, different weighting methods are used for different ground point distributions and numbers.

[0177] S7, locates the vegetation boundary position according to the vegetation mask, and selects the ICESat-2 ATL08 control points within the image range. The vegetation height difference correction surface is calculated based on the inverse distance weighted method of the ATL08 control points, and automatic elevation editing is performed.

[0178] In this embodiment, this step first reconstructs the resolution of the obtained vegetation mask product, selects an appropriate interpolation method to make it the same as the SAR image resolution in the study area, removes scattered small woodlands, and fills the bare land in the large woodlands; uses the generated vegetation mask as the forest boundary control, and screens the ICESat-2ATL08 control points within the image range, calculates the vegetation height difference correction surface based on the inverse distance weighted calculation of the ATL08 control points, and performs elevation editing on the vegetation.

[0179] Specifically, in this embodiment, the resolution of the TanDEM-X Forest / Non-ForestMap vegetation mask product (abbreviated as TanDEM-X vegetation mask) provided by the German Aerospace Center (DLR) is reconstructed, and a suitable interpolation method is selected to make it the same as the SAR image resolution in the study area. There are woodlands of different sizes in the TanDEM-X vegetation mask. In order to facilitate subsequent processing, a certain threshold is set in this embodiment to remove scattered small woodlands and fill the bare land in large woodlands, thereby generating a mask map for vegetation editing. The generated vegetation mask is used as the forest boundary control, and the ICESat-2ATL08 control points of the vegetation area in the vegetation mask are selected. The vegetation height difference correction surface is calculated based on the ATL08 control points, and the vegetation in the original DSM is subjected to elevation editing processing. Finally, the DEM after filtering out the vegetation height is output.

[0180] The inverse distance weighted process is:

[0181] Calculate the distance from an unknown point to all points:

[0182]

[0183] Calculate weights:

[0184]

[0185] Compute the values at the interpolation points:

[0186]

[0187] S8, considering the potential edge differences before and after processing, adopts the adaptive smooth transition method to obtain a DEM with reasonable transition information.

[0188] In this embodiment, considering that different algorithms are used to edit the elevations of different land features, elevation jumps may occur at the edges between different types of land features, the edge positions of the processing areas are indicated by the above-mentioned multiple mask files, and the smooth edge method is used to complete the edge transition processing to eliminate the elevation jump phenomenon.

[0189] In summary, the method provided in this embodiment makes full use of SAR intensity images, SAR coherence information, satellite imaging information, and reliable public external data, aiming to fully automatically complete the editing process from InSAR-DSM to DEM. It can effectively and quickly complete the identification and repair of InSAR-DSM data anomalies, and on this basis, fully automatically identify the distribution of water bodies and built-up areas. While ensuring the integrity of InSAR-DSM data, it can realize fully automatic elevation editing of water bodies and built-up areas respectively. The entire process structure is clear and has the advantages of clear process, simple implementation, large processing range, and high degree of automation.

[0190] The following is a case study to further illustrate the effectiveness of this method. The experimental area is the Macao Special Administrative Region of the People's Republic of China and Zhuhai City, Guangdong Province (21°N~22.53°N, 113.25°E~113.64°E). The stripe-mode image of the experimental area is used to obtain a SAR intensity image and a raw DSM data. The data details are shown in Table 2. The experimental data are shown in Table 2. Figure 3 As shown, Figure 3 (a) is the original DSM data of the experimental area, Figure 3 (b) is the intensity image of the experimental area. Figure 3 (c) is the coherence diagram of the experimental area.

[0191] Table 2 Detailed information of experimental data

[0192]

[0193] The geometric distortion results of the data in this area are as follows: Figure 4 As shown, the white area is the result of geometric distortion.

[0194] The automatic filling effect of data holes in the data of this area was tested. Figure 5 It shows that after automatic hole filling, the terrain product with holes is completely filled into a hole-free terrain product.

[0195] The water bodies in the data of this area were extracted. Figure 6 Shows a schematic diagram of the water extraction effect. Figure 6 (a) is the standard water mask result published by DLR. Figure 6 (b) is the water extraction result of this embodiment.

[0196] The building area on the data in this area was extracted. Figure 7 The figure shows the effect of building extraction. Figure 7 (a) is the standard building mask result released by DLR. Figure 7 (b) is the building extraction result of the present invention.

[0197] The water body classification results are as follows Figure 8 As shown in the figure, it can be seen that water bodies are divided into three categories: oceans, lakes and rivers. The water body editing results are as follows Figure 9 As shown, it can be seen that there is a lot of noise in the water body in the original DSM, or the elevation value is a null value (outlier). After the water body elevation editing process, the oceans, lakes and rivers are given elevation values close to their shore pixels, and the effect is very smooth.

[0198] Select a building area in the region as ROI for effect display, such as Figure 10 As shown, it can be seen that Figure 10 (b) There are also some terrain undulations caused by buildings. Figure 10 (c) is much smoother, and it can be clearly seen that the elevation information of the building has been filtered out.

[0199] The elevation filtered out by vegetation editing is as follows Figure 11 As shown, it can be seen that after vegetation editing, the elevation of the vegetation area is reduced, indicating that the elevation information of the vegetation has been filtered out.

[0200] The overall effect of automatic editing from InSAR-DSM to DEM is as follows Figure 12 As shown, Figure 12 (a) is the effect before automatic editing. Figure 12 (b) is the result after automatic editing. It can be clearly seen that the DEM after automatic editing is not only more complete than the DSM before automatic editing, but also has good elevation editing effects in the areas where various land features belong.

[0201] The quantitative accuracy of hole filling in digital terrain products can be evaluated from two aspects from the mainstream point of view. The first aspect is the change in the hole rate, and the second aspect is the change in the root mean square error (RMSE) of the data.

[0202] Using case data, the changes in void rate before and after automated void filling are statistically analyzed, as shown in Table 3.

[0203] Table 3 Statistics of void ratio changes

[0204]

[0205] The RMSE is calculated as follows:

[0206]

[0207] Among them, h ij is the elevation of each point, h truth is the ground truth elevation for reference.

[0208] With the support of high-precision LiDAR data as verification data, the RMSE change statistics were performed (excluding the water elevation). The accuracy statistical results are shown in Table 4.

[0209] Table 4 Statistics of RMSE changes

[0210]

[0211] After a fully automated hole filling process, this method can achieve better data accuracy.

[0212] In this quality assessment scheme, the accuracy of feature extraction is evaluated using internal and external accuracy. Internal accuracy is calculated using the statistical standard deviation within a class, which can be used to determine the difference within a class. The smaller the difference within a class, the higher the classification accuracy. The specific calculation formula is as follows:

[0213]

[0214]

[0215] Among them, x i is the gray value of the i-th pixel belonging to a certain category, is the average grayscale value of pixels belonging to this category, and σ is the internal accuracy.

[0216] External accuracy is measured using overall accuracy and the Kappa coefficient. For the building category, the 2019 global built-up area mask file, produced by the World Settlement Footprint (WSF) released by DLR in 2022 based on Sentinel data, was used as the actual result for comparison. For the water category, the global water mask file (WAM) released by DLR was used as the actual result for comparison. For a given classification result, assume that the classification result is as shown in the following table:

[0217] Table 5 Classification results

[0218]

[0219] Strictly speaking, the overall accuracy OA is defined as the ratio of the number of correctly extracted objects to the total number of objects. The specific calculation formula is as follows:

[0220]

[0221] Among them, TP, TN, FP, and FN are the numbers of pixels that meet the classification requirements. The Kappa coefficient K is a ratio that represents the reduction in errors compared to completely random classification. Its value usually falls between 0 and 1. Its numerical range is as follows: 0.0-0.20 is very low consistency; 0.21-0.40 is general consistency; 0.41-0.60 is moderate consistency; 0.61-0.80 is high consistency; 0.81-1.0 is almost perfect consistency. The specific calculation formula is as follows:

[0222]

[0223]

[0224] The accuracy evaluation of the feature extraction part is mainly divided into two aspects: intra-class standard deviation evaluation and overall accuracy and Kappa coefficient evaluation.

[0225] The main steps of intra-class standard deviation evaluation are: 1) obtaining the object classification results; 2) counting the grayscale values of the intensity map belonging to a certain category; 3) performing numerical analysis on the obtained grayscale values belonging to the category, calculating the standard deviation of this batch of data, and determining the standard deviation as the result as the internal accuracy of the extraction of this type of object.

[0226] The main steps for evaluating the overall accuracy and Kappa coefficient are as follows: 1) obtaining a standard mask result based on the scope of the study area and sampling it to the same accuracy as the extraction result; 2) calculating the ratio of the number of correctly classified pixels (including pixels belonging to built-up areas and non-built-up areas) to the total number of pixels in the image to obtain the overall accuracy of the method; and calculating the number of pixels in the standard result and classification result of each category to obtain the Kappa coefficient of the method.

[0227] An accuracy assessment of the ground feature extraction using the above scheme revealed an internal accuracy of 63.576 for the building category and external accuracy of 0.92 and 0.72, respectively. This indicates a high degree of consistency between the extracted building category and the true results. The internal accuracy for the water category was 13.855, while the external accuracy was 0.96 and 0.90, respectively. This indicates near-perfect consistency between the extracted water category and the true results.

[0228] This method aims to eliminate the noise and outliers (including null values) in the original DSM while ensuring the authenticity of the edited water body elevation, so that the edited water body elevation presents a smooth change effect, which will lead to a decrease in the amplitude of the water body elevation value change. Therefore, the standard deviation change before and after the water body elevation is edited can be used to quantitatively evaluate the effect of water body elevation editing.

[0229] The table below compares the standard deviations of the oceans, lakes, rivers, and the entire water body in the experimental area before and after elevation editing. It can be seen that after water body editing, the standard deviations of the elevations of various water bodies and the entire water body decreased. The standard deviations of oceans and rivers increased significantly, indicating that water body elevation editing has achieved very positive results.

[0230] Table 6 Comparison of elevation standard deviation before and after water body editing

[0231]

[0232] LiDAR point cloud data covering the Macao Special Administrative Region of the People's Republic of China and parts of Zhuhai City, Guangdong Province, collected in 2010, with a spatial resolution of 5m and a vertical accuracy of 1m, was compared with ground points identified using the proposed method from InSAR data to verify the accuracy of ground point extraction. Considering the significant differences between LiDAR and InSAR data, the ground points identified using the algorithm from the LiDAR data were compared with the ground point extraction results from the InSAR data for the same area.

[0233] A 2×2 confusion matrix is used to display the ground point extraction results. The elements in the first row and first column of the confusion matrix are true positives (TP); the elements in the first row and second column are false negatives (FN); the elements in the second row and first column are false positives (FP); and the elements in the second row and second column are true negatives (TN). Based on the confusion matrix, the precision (PPV) and false positive rate (FPR) are used as evaluation indicators to quantitatively analyze the results. The calculation formula is as follows:

[0234]

[0235]

[0236] The calculation results are shown in Table 7:

[0237] Table 7 Evaluation indicators of ground point extraction accuracy

[0238]

[0239] As can be seen from Table 7, the proportion of ground points extracted by this method that are within the range of real ground points is 90.96%. The proportion of ground points extracted by this method that are within the range of real non-ground points is only 11.02%. Overall, the method of extracting ground points of the present invention is considered reliable.

[0240] The above is a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as within the scope of protection of the present invention.

Claims

1. A method for automatically editing and acquiring InSAR-DSM to DEM, characterized in that: The steps include: S1, combines prior DEM information, SAR image intensity information, InSAR interferometric coherence coefficient information and satellite imaging parameters to accurately identify geometric distortion areas; identify InSAR-DSM holes and low-quality areas, and together with the geometric distortion areas, form a DSM data anomaly distribution mask; S2, based on the anomaly distribution mask of the DSM data, the InSAR-DSM hole area is obtained, and the hole is filled by combining the public terrain product and the geospatial interpolation method to obtain the complete InSAR-DSM; S3, based on the SAR coherence map and SAR intensity image, combined with external DEM information, the single threshold segmentation method is used to extract the SAR feature representation of the water body and obtain the water body mask; S4, based on the SAR coherence map and SAR intensity image, the SAR feature representation of the building is extracted based on the image texture features and machine learning to obtain the building area mask; S5, locates the water body position according to the water body mask, and classifies the water body into three categories according to its size and shape: rivers, lakes, and oceans. Based on the different characteristics of different water bodies in terms of elevation changes, automatic elevation editing is performed on each of them; S6, locates the building position according to the built-up area mask, divides the ground points and non-ground points, and performs automatic elevation editing through interpolation based on the distribution characteristics of ground points in built-up areas of different terrains and sizes; S7, locate the vegetation boundary position according to the vegetation mask, and screen the control points within the image range, calculate the vegetation height difference correction surface based on the inverse distance weighted control points, and perform automatic height editing; S8, based on the potential edge differences, an adaptive smooth transition method is adopted to obtain a DEM with reasonable transition information.

2. The automatic editing and acquisition method of InSAR-DSM to DEM according to claim 1 is characterized in that: S1 specifically includes: Obtain satellite incident angle and azimuth information from SAR satellite imaging parameter files; Reconstruct the resolution of the public terrain products obtained and make the SAR image resolution within the study area the same by interpolation; Calculate the slope, aspect, and elevation of the InSAR-DSM pixel by pixel, and determine whether it is a shadow or overlap area based on its relationship with the satellite's incidence angle and azimuth parameters; Using SAR imaging parameters, we can derive a terrain product data quality indicator map to identify low-quality areas. The missing areas of InSAR-DSM data are identified, and together with the geometric distortion areas and low-quality areas, they constitute the data anomaly areas, forming the DSM data anomaly distribution mask.

3. The automatic editing and acquisition method of InSAR-DSM to DEM according to claim 2 is characterized in that: S2 specifically includes: Holes formed in areas indicated by abnormal distribution masks of DSM data are eliminated. Based on the different correlation strengths of surrounding terrain data under different terrain complexity and the different terrain correlation strengths caused by holes of different sizes, the holes to be filled are adaptively divided into low-difficulty areas and high-difficulty areas. The difficulty of filling holes is determined based on the complexity of the local terrain. For low-difficulty areas, geospatial interpolation is used. For high-difficulty areas, public external terrain products are introduced and the corresponding regional data is selected as an external reference to assist in the automatic filling of holes.

4. The automatic editing and acquisition method of InSAR-DSM to DEM according to claim 3 is characterized in that: S3 specifically includes: SAR intensity images and SAR coherence maps were used to identify water bodies in the study area; Based on the SAR intensity image, the grayscale histogram distribution information of the intensity image is statistically analyzed, and the water system is identified according to the bimodal method. The lowest point of the image valley between the water body and the land is found as the critical value, and the image is segmented based on this threshold. The SAR intensity image, SAR coherence map and slope map generated based on external DEM information are combined by band combination to obtain a pseudo-color image of the study area. This image contains the above three types of information, and water bodies are extracted from the pseudo-color image.

5. The automatic editing and acquisition method of InSAR-DSM to DEM according to claim 4 is characterized in that: S4 specifically includes: Perform filtering preprocessing on the SAR intensity image, extract the building area by counting the image texture features and determining the extraction threshold; The texture features of SAR intensity images are extracted using gray-level co-occurrence matrix. The appropriate feature expression operators are screened through principal component analysis and Bhattacharyya distance evaluation, and the feature expression operators with strong feature expression capabilities are integrated. The fused features are classified through machine learning classification, and building extraction is completed with the assistance of external DEM information.

6. The automatic editing and acquisition method of InSAR-DSM to DEM according to claim 5 is characterized in that: S5 specifically includes: Using the extracted water mask files and global high-resolution coastline vector data, the oceans are classified according to their geographical locations; the remaining water bodies are classified into lakes and rivers according to their shapes and areas. Based on the water body classification results, elevation editing is performed on the three types of water bodies: oceans, lakes and rivers, according to the different characteristics of different types of water bodies and the characteristics of actual data.

7. The automatic editing and acquisition method of InSAR-DSM to DEM according to claim 6 is characterized in that: S6 specifically includes: The extracted built-up area mask files are morphologically preprocessed, and the built-up areas are divided into large-area built-up areas and small-area built-up areas based on the degree of terrain relief and area. The building heights of the large-area built-up areas and the small-area built-up areas are adaptively filtered out respectively. Based on the echo signal and geometric characteristic factors, a joint constraint condition is constructed to distinguish ground points from non-ground points. InSAR ground points are extracted by combining error statistics and point cloud cleaning. Based on the selected InSAR ground points, the elevation recovery process is achieved by using an adaptive interpolation method based on the different distribution characteristics of ground points in urban areas of different areas and the original terrain characteristics. Based on the texture features of the original DSM, texture restoration is performed after elevation restoration.

8. The automatic editing and acquisition method of InSAR-DSM to DEM according to claim 7 is characterized in that: S7 specifically includes: The obtained vegetation mask product was reconstructed and interpolated to the same resolution as the SAR image in the study area, and the scattered small forests were removed and the bare land in the large forests was filled. The generated vegetation mask was used as the forest boundary control, and the ICESat-2ATL08 control points within the image range were screened; The vegetation elevation correction surface is calculated based on the inverse distance weighted method of ATL08 control points, and the vegetation elevation is edited.

Images