An automatic hole filling method for InSAR terrain products

Abstract

The present invention provides an automated method for filling holes in InSAR terrain products, including: data preprocessing; excluding invalid data areas; adaptively dividing the data into low-difficulty areas and high-difficulty areas; using an adaptive geospatial interpolation method to fill holes in the low-difficulty areas; using automated external data to assist in filling holes in the high-difficulty areas; using an adaptive smoothing transition method to completely fill data holes in the entire terrain product; and automatically generating and saving a processing area indicator map. The present invention fully utilizes the interpolation effect advantages of the classic geostatistical method Kriging interpolation, reasonably draws on external terrain data trends to help establish some missing parts of its own terrain data product, uses an adaptive smoothing method to resolve data inconsistencies that may be caused by external data intervention, and organizes the above processes in an automated manner. It has the advantages of simple implementation, a large processing range, and a high degree of automation.

Description

Technical Field

[0001] The present invention relates to the field of DEM production technology based on InSAR technology, and in particular to an automatic hole filling method for InSAR terrain products. Background Art

[0002] Since the 1980s, with the advancement of earth science, global digitization, and informatization, exploring global change and sustainable development has become a key issue. Since the vast majority of information in our daily lives is spatially related, the production of global digital terrain products has become increasingly significant. Commonly used digital terrain products include DEMs (digital elevation models) and DSMs (digital surface models).

[0003] Common acquisition methods for digital terrain products include field surveying, photogrammetry, laser scanning, and Interferometric 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 producing a global DEM, InSAR technology is currently one of the best approaches, considering both time and cost.

[0004] However, due to the presence of surface obstructions, such as trees and buildings, as well as limitations of the SAR system itself, InSAR technology can suffer from data loss and noise interference during terrain data acquisition and processing. This missing data can lead to incomplete and inaccurate terrain data, thus affecting subsequent terrain analysis and applications.

[0005] Currently, there has been some research on the reconstruction of digital terrain products, primarily focusing on the use of interpolation methods, such as nearest neighbor interpolation and inverse distance weighted interpolation. However, simple interpolation alone is clearly insufficient for the complex and ever-changing terrains of various regions around the world, especially in complex areas where data gaps are prone to occur. Furthermore, the automated filling process for terrain data products has received insufficient research and attention.

[0006] By automatically filling data gaps, appropriately referencing existing external terrain data, and rationally combining interpolation and external data reference, the accuracy and completeness of terrain data can be improved, providing more reliable support for terrain analysis and applications. Therefore, researching and developing automated methods for filling data gaps in InSAR terrain products is of great significance and practical application value. Summary of the Invention

[0007] The purpose of the present invention is to address the deficiencies in the above-mentioned background technology, combine SAR imaging geometric parameters with externally disclosed terrain products, automatically fill the data gaps in the terrain products produced by InSAR, complete the range of terrain products, and improve the quality of terrain products.

[0008] In order to achieve the above object, the present invention provides an automatic hole filling method for InSAR terrain products, comprising the following steps:

[0009] S1, preprocessing terrain product data and auxiliary data;

[0010] S2, using the backscatter intensity characteristics of SAR intensity images to exclude invalid data areas;

[0011] S3, based on the different levels of terrain undulation and the size of the holes, comprehensively considers the difficulty of filling holes and adaptively divides them into low-difficulty areas and high-difficulty areas;

[0012] S4: Adaptive geospatial interpolation is used to fill gaps in low-difficulty areas; automated external data is used to assist in filling gaps in high-difficulty areas.

[0013] S5: For inconsistencies between the terrain data filled with external data and the original terrain data, an adaptive smooth transition method is used to completely fill the data holes in the entire terrain product.

[0014] S6, according to the automatic processing process of the terrain data holes, combined with the original terrain data, automatically gives and saves the processing area indication map.

[0015] Furthermore, in S1, based on the SAR imaging geometric information and publicly available terrain data, the shadow and overlap conditions caused by the SAR side-view imaging process are established, and the area where geometric distortion occurs is calculated; based on the calculated geometric distortion information, the corresponding regional data in the terrain product data is eliminated, and the remaining non-geometric distortion areas are retained as reliable data areas.

[0016] Furthermore, in S1,

[0017]

[0018] Among them, (i, j) represents the DEM grid, χ ij is the slope at each point, θ ij is the incident angle of each point during SAR imaging;

[0019] In the data, let the intensity of each point be P ij , the geometric distortion image value is D ij , under the following conditions, it is reserved as a valid data area:

[0020] P ij >0∩D ij =0(2).

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

[0022] S31, calculating a terrain-related factor that can represent the complexity of the terrain based on the terrain data, and calculating the elevation gradient:

[0023]

[0024] Among them, h i represents the elevation value, p and q are the horizontal and vertical gradients respectively, w is the ground resolution, r and s are the intermediate values of the second-order gradient; calculate the surface curvature C m , which indicates the terrain complexity of the DEM:

[0025]

[0026]

[0027] S32, calculate the interpolation effect A according to the size of the hole ij :

[0028]

[0029] Where α is the scale factor, (m, n) are the coordinates of the known point relative to the unknown point (i, j);

[0030] S33, combining S31 and S32, obtains the adaptive hole filling difficulty division formula:

[0031]

[0032] 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.

[0033] Furthermore, in S4, the low-difficulty areas are processed one by one, and the valid data points within the preset range around the low-difficulty areas are aggregated, and the elevations and positions of representative points with uniform distribution are screened out; the spatial semi-variable function is determined based on the elevations and positions of the screened points as known points, and the Kriging model is fitted; based on the fitted Kriging model, the low-difficulty areas to be filled are predicted, and the filling of one area is completed; it is applied to all low-difficulty filling areas to complete the filling of holes in the low-difficulty areas.

[0034] Furthermore, in S4, adaptive geospatial interpolation methods are used to fill holes in low-difficulty areas, including the following:

[0035] Select the elevation points around the cavity as preliminary candidate points;

[0036] From the candidate points, select reference points based on uniform distribution of directions for reference calculation:

[0037]

[0038] Kriging interpolation is used to fill the gaps, expressing the value of the unknown location as a weighted average of the values of the nearest known data points, with the weights varying with distance and direction:

[0039]

[0040] in, is the estimated value of the point (x0, y0), that is, h0 = (x0, y0); λ i is the weight coefficient, obtained by fitting the Kriging model.

[0041] Furthermore, in s4, for difficult areas, external reference terrain products are selected based on the geographic coordinate range. After unification of resolution and coordinate system and data registration, the external data are fitted to the original InSAR terrain data. The local benchmark offset H between the external data and the original data is calculated. offset , offset the external data benchmark to make it consistent with the benchmark of the original data; all the holes to be filled in the area are filled with the assistance of external data.

[0042] Furthermore, adaptive bilateral filtering is used in S5 to perform local smooth transition. During the filtering process, the distance weight is adjusted according to the elevation value difference between pixels to achieve adaptive smoothing.

[0043] Furthermore, in terrain data processing, the parameters include the smoothing radius win, the elevation difference threshold T0, and the elevation standard deviation STD; the smoothing radius represents the neighborhood size considered by the bilateral filter during the smoothing process, and β is a set constant; the elevation difference threshold is used to control the distance weight λ d When the elevation difference T is less than the threshold, the distance weight becomes smaller and the smoothing effect is stronger; when the elevation difference is greater than the threshold, the distance weight becomes larger and the smoothing effect is weaker.

[0044]

[0045]

[0046] The elevation standard deviation is used to control the Gaussian kernel size of the filter, which is set according to the resolution and noise characteristics of the terrain data, including:

[0047] Calculate the elevation standard deviation, h, of image or DEM data avgis the mean elevation, n is the total number of points;

[0048]

[0049] Determine a baseline Gaussian kernel size based on the resolution and noise characteristics of the image or DEM data; calculate the actual Gaussian kernel size based on the elevation standard deviation and the baseline Gaussian kernel G0 size;

[0050] G=G0+μ·STD (13)

[0051] Where μ is the adjustment coefficient, and G is the actual Gaussian kernel after adaptive adjustment.

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

[0053] The InSAR terrain product hole filling method provided by this invention is based on SAR imaging geometric parameters and externally available terrain products. Starting from InSAR terrain product data with holes, it fully utilizes the interpolation effect advantages of the classic geostatistical method Kriging interpolation, reasonably draws on the trends of external terrain data to help establish some missing parts of its own terrain data product, and timely proposes an adaptive smoothing method to resolve data inconsistencies that may be caused by external data intervention. The above process is organized in an automated manner, making the hole filling process fully automated. The entire process has a clear structure and has the advantages of simple implementation, a large processing range, and a high degree of automation.

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

[0055] Figure 1 It is a flowchart of the present invention;

[0056] Figure 2 This is a schematic diagram of the geospatial interpolation reference point screening of the present invention;

[0057] Figure 3 The data used in the specific embodiment of the present invention are shown in Figure 1, where (a) is the original DSM terrain product with holes obtained by InSAR processing of the case data, and (b) is the SAR intensity image of the case data;

[0058] Figure 4 Figure 1 shows the effect of automatic hole filling in panoramic data according to a specific embodiment of the present invention, where (a) shows the original InSAR terrain product data before hole filling, and (b) shows the complete InSAR terrain product data after hole filling.

[0059] Figure 5 for Figure 4A diagram showing the effect of automatic hole filling in a locally enlarged area, where (a) is the terrain product before local hole filling, and (b) is the terrain product after local hole filling;

[0060] Figure 6 1 is a diagram showing the area selected for local LiDAR accuracy verification in an embodiment of the present invention, where (a) is the terrain product of the LiDAR verification area before hole filling, and (b) is the terrain product of the LiDAR verification area after hole filling. DETAILED DESCRIPTION

[0061] 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.

[0062] 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.

[0063] 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.

[0064] like Figure 1As shown, an embodiment of the present invention provides an automatic hole filling method for InSAR terrain products, comprising the following steps:

[0065] S1, preprocessing terrain product data and auxiliary data.

[0066] In this embodiment, shadows and overlap caused by the SAR side-view imaging process are first established based on SAR imaging geometry and publicly available terrain data, and areas of geometric distortion are calculated. Based on this calculated geometric distortion information, the corresponding areas in the terrain product data are removed, retaining the remaining areas without geometric distortion as reliable data.

[0067] It should be noted that geometric distortion means that pixel information is unreliable, and therefore the calculated elevation values are also unreliable. SAR images are not regular rectangles in the geographic coordinate system, and data is stored in a rectangular grid, so it is necessary to mark invalid data areas.

[0068] The calculation of SAR geometric distortion requires obtaining SAR images and corresponding terrain data, calculating SAR system parameters, ground object height and pixel position, and finally comparing the pixel position with the actual ground object position to obtain the size and direction of SAR geometric distortion.

[0069] For slopes facing the SAR sensor, when the local terrain slope angle is less than the local angle of incidence, the slope appears shorter in the SAR image than flat terrain, resulting in poor range resolution on these slopes. Furthermore, when the local terrain slope angle exceeds the local angle of incidence, inverted images appear at the base and top of the slope, a phenomenon known as active stacking. For slopes facing away from the SAR sensor, when the local terrain slope angle is less than the complementary angle of the local angle of incidence, the range resolution of these slopes in SAR images is better than that of flat terrain. However, when the local terrain slope angle exceeds the complementary angle of the local angle of incidence, the SAR signal on such steep slopes is completely blocked from the mountain itself. This is the active shadowing effect that causes slopes to appear darker in SAR images.

[0070]

[0071] Among them, (i, j) represents the DEM grid, χ ij is the slope at each point, θ ij is the incident angle of each point during SAR imaging.

[0072] In the data, let the intensity of each point be P ij , the geometric distortion image value calculated above is D ij , can be reserved as a valid data area under the following conditions:

[0073] P ij>0∩D ij =0 (2)

[0074] S2, using the backscatter intensity characteristics of SAR intensity images to exclude invalid data areas.

[0075] In this embodiment, according to the backscatter intensity characteristics, the empty areas in the SAR intensity image are areas where no information is obtained. It is neither possible nor necessary to recover terrain information, and such invalid data areas need to be excluded.

[0076] It should be noted that the SAR intensity image indicates the intensity of each pixel. There is no intensity in the invalid data area. Therefore, the presence or absence of intensity can indicate the validity of the data. By combining the data quality damaged area and the data valid area, the data valid area can be retained.

[0077] S3, based on the different degrees of terrain undulation and the size of the holes, comprehensively considers the difficulty of hole filling and adaptively divides the area into low-difficulty area and high-difficulty area.

[0078] It's important to note that data holes in terrain products often appear due to low coherence or geometric distortion. Depending on imaging and terrain conditions, holes often exhibit distinct characteristics. First, when holes occur in rugged terrain, elevation variations are significant even within a small area, and using existing surrounding elevations often cannot restore the missing elevation values. Conversely, in relatively flat areas, elevation variations are consistent over a larger area, and using surrounding elevations farther away can restore the missing elevations. Second, when a contiguous hole is large, the influence of surrounding elevations decreases toward the center of the hole, or even becomes largely irrelevant. Conversely, when holes are small, surrounding elevations provide a significant reference for the elevation within the hole. When filling holes, these two influencing factors should not be considered in isolation, but rather jointly. Furthermore, with the support of the hole-filling algorithm, processing areas with varying degrees of difficulty should be adaptively divided to facilitate subsequent tiered processing.

[0079] Specifically in this embodiment, first, a terrain-related factor that can represent the complexity of the terrain is calculated based on the terrain data, and the elevation gradient is calculated:

[0080]

[0081] Table 1 Schematic diagram of adjacent elevations

[0082]

[0083]

[0084] Among them, h iThe 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:

[0085]

[0086]

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

[0088]

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

[0090] Combining the above two factors, we can get the adaptive hole filling difficulty division formula:

[0091]

[0092] 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.

[0093] S4: Adaptive geospatial interpolation is used to fill gaps in low-difficulty areas. Each low-difficulty area is processed one by one, and valid data points within a certain range around the low-difficulty area are aggregated. After screening, the elevation and location of representative points with uniform distribution are retained. The spatial semivariogram is determined using the elevation and location of the screened points as known points, and a Kriging model is fitted. Based on the fitted Kriging model, the low-difficulty area to be filled is predicted, completing the gap filling of one area. Based on the above process, this process is applied to all low-difficulty filling areas to complete the gap filling work in the low-difficulty area.

[0094] Low-difficulty areas typically have small holes or relatively consistent terrain trends. Therefore, integrating valid surrounding terrain data is sufficient to leverage the proximity correlation of terrain to recover missing data. Among geospatial interpolation methods, Kriging is known for its effectiveness, making it a reliable method for filling low-difficulty holes and ensuring more accurate results.

[0095] It should be noted that when using the geospatial interpolation method, the first problem faced is the selection of known points. Not all nearby points are suitable for reference. When selecting points, sampling should be done as evenly as possible around the hole to ensure the uniformity of the points and further ensure the reliability of the restored elevation. Figure 2 Kriging interpolation is a geostatistical method that predicts values at unknown locations based on known locations and their spatial relationships. Its basic principles include collecting data, selecting an interpolation model, calculating distances, choosing a variogram, estimating model parameters, predicting values at unknown locations, and evaluating accuracy. The accuracy of kriging interpolation predictions depends on factors such as the number and spatial distribution of known locations, the choice of interpolation model, and the parameters of the variogram.

[0096] Specifically in this embodiment, first, elevation points around the cavity are selected as preliminary candidate points.

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

[0098]

[0099] 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, the search range R should be expanded in this case.

[0100] 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:

[0101]

[0102] 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.

[0103] For high-difficulty areas, automated external data-assisted gap filling is used. Based on the terrain product data that has already been filled in the low-difficulty areas in the first part, external reference data for the high-difficulty filling areas is selected from publicly available external free terrain product data. Based on the surrounding data of the high-difficulty area and the local external reference data, the offset between the original data and the external data is determined, and the high-difficulty terrain gaps are filled by introducing terrain change trends.

[0104] It should be noted that, aside from low-difficulty areas, InSAR terrain products can also contain large-scale or extremely complex terrain, where geospatial interpolation methods will not be able to effectively recover the gaps. In such cases, it is possible to use globally available and freely available terrain product data, selecting strictly corresponding regional data as an external reference, and then filling in the gaps based on this.

[0105] It should be noted that in this embodiment, high-difficulty areas refer to areas where it is difficult to accurately restore the hole elevation based solely on existing surrounding elevation data. Therefore, under the premise that external public terrain data is freely available globally, local terrain reference information can be introduced to assist in filling the gaps. When filling the gaps, the external data cannot be simply "stitched" with the existing data. Instead, the local terrain change trends of the external data should be superimposed on the existing data based on the existing data, achieving true external data-assisted gap filling.

[0106] Specifically in this embodiment, the connected range of each block in the divided high-difficulty hole filling area is processed. First, the corresponding external reference terrain product is selected according to the geographic coordinate range. After unification of resolution and coordinate system and data registration, the external data is aligned with the original InSAR terrain data. Second, the local benchmark offset H between the external data and the original data is calculated. offset , make a corresponding size of benchmark offset on the external data to make it consistent with the benchmark of the original data. All the holes to be filled in the area are completed with the assistance of external data.

[0107] S5: For the phenomenon that the terrain data filled with the assistance of external data may be inconsistent with the original terrain data at the junction, an adaptive smooth transition method is adopted to solve it, and finally the data gaps in the entire terrain product are completely filled, thereby improving the integrity of the data.

[0108] It should be noted that after hole filling, the introduction of external data will inevitably cause some data inconsistencies in the edge areas. The most common approach is to use filters to achieve a smooth transition, but the terrain is complex and changeable, which means that data errors and terrain changes are common. Therefore, the filter must not destroy the original important terrain features while adapting to the complex and changing terrain environment. Adaptive smooth transition methods based on different data conditions are suitable solutions. Bilateral filtering is a filtering method based on Gaussian kernels and distance weights. It can smooth the image while preserving the edge information. Based on bilateral filtering and combined with the actual terrain data, adaptability can be introduced to adapt to the different terrain conditions and smoothing requirements of the edge area by adaptively adjusting the filter parameters.

[0109] Specifically in this embodiment, adaptive bilateral filtering is used to perform local smooth transition, and adaptive smoothing is achieved by adjusting the distance weight according to the elevation value difference between pixels during the filtering process. In terrain data processing, parameters usually include smoothing radius (win), elevation difference threshold (T0) and elevation standard deviation (STD). The smoothing radius represents the neighborhood size considered by the bilateral filter during the smoothing process, which can be set according to the characteristics of the terrain data and application requirements, and β is a set constant. The elevation difference threshold is used to control the distance weight (λ d ), when the elevation difference (T) is less than the threshold, the distance weight becomes smaller and the smoothing effect is stronger; when the elevation difference is greater than the threshold, the distance weight becomes larger and the smoothing effect is weaker.

[0110]

[0111]

[0112] The elevation standard deviation is used to control the Gaussian kernel size of the filter and is usually set according to the resolution and noise characteristics of the terrain data. Specifically, it includes:

[0113] Calculate the elevation standard deviation, h, of image or DEM data avg is the mean elevation, and n is the total number of points.

[0114]

[0115] Determine a baseline Gaussian kernel size based on the resolution and noise characteristics of the image or DEM data. Generally speaking, the Gaussian kernel size should be related to the data resolution, meaning that a smaller Gaussian kernel is used for higher-resolution data and a larger Gaussian kernel is used for lower-resolution data. Furthermore, the noise characteristics also affect the Gaussian kernel size setting; if the data contains a large amount of noise, a larger Gaussian kernel should be used to smooth the data.

[0116] The actual Gaussian kernel size is calculated based on the elevation standard deviation and the base Gaussian kernel G0 size.

[0117] G=G0+μ·STD (13)

[0118] Where μ is the adjustment coefficient, and G is the actual Gaussian kernel after adaptive adjustment.

[0119] Through adaptive parameter adjustment, it can be adapted to tasks under different terrain conditions, striking a balance between smooth transition and protection of important terrain features.

[0120] S6, based on the terrain data hole automatic processing process, combined with the original terrain data, automatically generates and saves a processing area indication map. Specifically, the processed area is automatically recorded using a mask file, and a processing area indication map is saved to indicate where the hole filling is performed.

[0121] In summary, the InSAR terrain product hole-filling method provided in this embodiment is based on SAR imaging geometric parameters and externally available terrain products. Starting from InSAR terrain product data with holes, it fully leverages the interpolation performance advantages of the classic geostatistical method, kriging interpolation. It rationally draws on external terrain data trends to help construct some missing parts of the native terrain data product. It also proposes an adaptive smoothing method to address data inconsistencies that may be caused by external data intervention. The above process is organized in an automated manner, making the hole-filling process fully automated. The entire process has a clear structure, simple implementation, a wide processing range, and a high degree of automation.

[0122] The following example further illustrates the effectiveness of this method. The Macao Special Administrative Region and Zhuhai City, Guangdong Province, are used as the experimental area (21°N-22.53°N, 113.25°E-113.64°E). The stripe-mode imagery of the experimental area, including a SAR intensity image and a raw DSM image, is used as the case study. The data details are shown in Table 2. The cropped experimental data are shown in Table 2. Figure 3 shown.

[0123] Table 2 Detailed information of experimental data

[0124]

[0125] The automatic filling effect of data holes on the scene data in the Macao Special Administrative Region was tested. Figure 3 The data used are shown, including terrain products and SAR intensity images. Figure 4 It shows that after the automatic hole filling, the terrain product with holes is completely filled into a hole-free terrain product. In order to show the filling effect details more clearly, the local enlarged filling result is shown in Figure 5 It can be clearly seen that after the hole filling, the original data holes are well filled, and there is no discontinuity of terrain such as cliffs and steps. From the visual effect, the accuracy is better.

[0126] 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.

[0127] Using case data, we calculated the changes in void rate before and after automated void filling:

[0128]

[0129] The RMSE is calculated as follows:

[0130]

[0131] Among them, h ij is the elevation of each point, h truth The ground truth elevation is used as a reference. In this case, some areas of this scene data (such as Figure 6 , Figure 6 (a) Before the cavity is filled, Figure 6 (b) after hole filling), with the support of high-precision LiDAR data as verification data, the RMSE change statistics were performed (excluding the water elevation), and the accuracy statistics results are shown in Table 3

[0132] Table 3 RMSE change statistics

[0133]

[0134] After the fully automated hole filling process, this method can achieve good data accuracy, so the method proposed in this invention shows good qualitative and quantitative accuracy in the automated hole filling terrain product.

[0135] Based on the same inventive concept, this embodiment further provides a computer-readable storage medium having a computer program stored thereon, which implements the aforementioned InSAR terrain product hole filling method when executed by a processor.

[0136] Computer-readable media include, but are not limited to, any type of disk (including floppy disks, hard disks, optical disks, CD-ROMs, and magneto-optical disks), ROMs, RAMs, EPROMs (Erasable Programmable Read-Only Memory), EEPROMs, flash memories, magnetic cards, or optical cards. In other words, computer-readable media include any medium that can store or transmit information in a form that can be read by a device (e.g., a computer).

[0137] The computer-readable storage medium provided in this embodiment has the same inventive concept and the same beneficial effects as the aforementioned method, and will not be described in detail here.

[0138] 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. An automatic hole filling method for InSAR terrain products, characterized by: The steps include: S1, preprocessing terrain product data and auxiliary data; S2, using the backscatter intensity characteristics of SAR intensity images to exclude invalid data areas; S3, based on the different levels of terrain undulation and the size of the holes, comprehensively considers the difficulty of filling holes and adaptively divides them into low-difficulty areas and high-difficulty areas; S3 specifically includes the following sub-steps: S31, calculating a terrain-related factor that can represent the complexity of the terrain based on the terrain data, and calculating the elevation gradient: Among them, h i represents the elevation value, p and q are the horizontal and vertical gradients respectively, w is the ground resolution, r and s are the intermediate values of the second-order gradient; calculate the surface curvature C m , which indicates the terrain complexity of the DEM: S32, calculate the interpolation effect A according to the size of the hole ij : Where α is the scale factor, (m,n) is the coordinate of the known point relative to the unknown point (i,j); S33, combining S31 and S32, obtains the adaptive hole filling difficulty division formula: 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; S4: Adaptive geospatial interpolation is used to fill gaps in low-difficulty areas; automated external data is used to assist in filling gaps in high-difficulty areas. S5: For inconsistencies between the terrain data filled with external data and the original terrain data, an adaptive smooth transition method is used to completely fill the data holes in the entire terrain product. S6, according to the automatic processing process of the terrain data holes, combined with the original terrain data, automatically gives and saves the processing area indication map.

2. The method for automatically filling holes in InSAR terrain products according to claim 1, characterized in that: In S1, based on the SAR imaging geometry information and publicly available terrain data, the shadow and overlap conditions caused by the SAR side-view imaging process are established, and the area where geometric distortion occurs is calculated; based on the calculated geometric distortion information, the corresponding regional data in the terrain product data is eliminated, and the remaining non-geometric distortion areas are retained as reliable data areas.

3. The method for automatically filling holes in InSAR terrain products according to claim 2, characterized in that: In S1, Among them, (i, j) represents the DEM grid, χ ij is the slope at each point, θ ij is the incident angle of each point during SAR imaging; In the data, let the intensity of each point be P ij , the geometric distortion image value is D ij , under the following conditions, it is reserved as a valid data area: P ij >0∩D ij =0。 4. The method for automatically filling holes in InSAR terrain products according to claim 1, characterized in that: In S4, the low-difficulty areas are processed one by one, and the valid data points within the preset range around the low-difficulty areas are aggregated, and the elevation and position of representative points with uniform distribution are screened out; the spatial semi-variable function is determined based on the elevation and position of the screened points as known points, and the Kriging model is fitted; based on the fitted Kriging model, the low-difficulty areas to be filled are predicted to complete the filling of an area; it is applied to all low-difficulty filling areas to complete the filling of holes in the low-difficulty areas.

5. The method for automatically filling holes in InSAR terrain products according to claim 4, characterized in that: In S4, adaptive geospatial interpolation methods are used to fill holes in low-difficulty areas, including the following: Select the elevation points around the cavity as preliminary candidate points; From the candidate points, select reference points based on uniform distribution of directions for reference calculation: Kriging interpolation is used to fill the gaps, expressing the value of the unknown location as a weighted average of the values of the nearest known data points, with the weights varying with distance and direction: in, is the estimated value of the point (x0, y0), that is, h0 = (x0, y0); λ i is the weight coefficient, obtained by fitting the Kriging model.

6. The method for automatically filling holes in InSAR terrain products according to claim 5, characterized in that: In S4, for difficult areas, external reference terrain products are selected based on the geographic coordinate range. After unification of resolution and coordinate system and data registration, the external data are fitted to the original InSAR terrain data. The local benchmark offset H between the external data and the original data is calculated. offset , offset the external data benchmark to make it consistent with the benchmark of the original data; all the holes to be filled in the area are filled with the assistance of external data.

7. The method for automatically filling holes in InSAR terrain products according to claim 6, characterized in that: S5 uses adaptive bilateral filtering for local smooth transition. During the filtering process, the distance weight is adjusted according to the elevation value difference between pixels to achieve adaptive smoothing.

8. The method for automatically filling holes in InSAR terrain products according to claim 7, characterized in that: In terrain data processing, the parameters include the smoothing radius win, the elevation difference threshold T0, and the elevation standard deviation STD; the smoothing radius represents the neighborhood size considered by the bilateral filter during the smoothing process, and β is a set constant; the elevation difference threshold is used to control the distance weight λ d When the elevation difference T is less than the threshold, the distance weight becomes smaller and the smoothing effect is stronger; when the elevation difference is greater than the threshold, the distance weight becomes larger and the smoothing effect is weaker. The elevation standard deviation is used to control the Gaussian kernel size of the filter, which is set according to the resolution and noise characteristics of the terrain data, including: Calculate the elevation standard deviation, h, of image or DEM data avg is the mean elevation, n is the total number of points; Determine a baseline Gaussian kernel size based on the resolution and noise characteristics of the image or DEM data; calculate the actual Gaussian kernel size based on the elevation standard deviation and the baseline Gaussian kernel C0 size; G=G0+μ·STD Where μ is the adjustment coefficient, and G is the actual Gaussian kernel after adaptive adjustment.

Images