A method for automatic selection of ground control points based on permanent scatterers

Abstract

The present invention relates to a method for automatically selecting ground control points based on permanent scatterers, comprising: obtaining multi-time series SAR images; extracting pixels that meet the permanent scatterer selection strategy from the multi-time series SAR images as ground control points; the permanent scatterer selection strategy includes: calculating the average coherence of each pixel, selecting pixel points whose average coherence is greater than a preset low coherence threshold; selecting pixel points whose amplitude is greater than the minimum amplitude mean; selecting pixel points whose amplitude deviation index is less than a preset amplitude deviation index threshold; selecting pixel points whose average coherence is less than a preset high coherence threshold; selecting pixel points whose slope is less than a preset terrain slope threshold; and further performing thinning processing according to the density, distance and terrain conditions of the points. Compared with the prior art, the present invention comprehensively considers the influencing factors such as phase stability, scattering intensity, terrain slope and GCP distribution, and the accuracy of DEM obtained in different regions is improved by about 20% to 30%, and the performance is best in plain areas.

Description

Technical Field

[0001] The invention relates to the technical field of ground control point selection, and in particular to an automatic ground control point selection method based on permanent scatterers. Background Art

[0002] Interferometric synthetic aperture radar (InSAR) is a high-precision technology that obtains phase information from single-look complex (SLC) images to extract ground elevation and monitor deformation information. The essence of InSAR generation of digital elevation model (DEM) is phase measurement and phase-to-height conversion. Since the working range of SAR is the microwave band, SAR can be used under any meteorological conditions and has the ability to work all day and all weather. Although InSAR has the advantages of large area and spatial continuous coverage, the process of InSAR generation of DEM is easily affected by various factors, such as the selection of interference pairs, baseline and atmospheric environment. Introducing high-quality ground control points (GCPs) in the process of track refinement and re-leveling can correct baseline parameters and remove the flat ground effect, effectively improving the accuracy of DEM. However, commonly used GCP selection methods such as manual selection may reduce the accuracy of DEM. It is necessary to study new methods to select appropriate GCPs for track refinement and re-leveling to obtain high-precision DEM.

[0003] A variety of methods, including manual selection, existing survey data, artificial corner reflectors, and permanent scatterers (PS), can be used to select GCPs for track refinement and re-leveling. Manual selection of GCPs for track refinement and re-leveling is a common method in sedimentation monitoring and DEM generation. Such GCPs are usually extracted from SAR images by operators. Since it is difficult to identify points of the same name in SAR images, the quality of manually selected GCPs is highly related to the operator's learning level of SAR images, familiarity with the study area, and SAR processing experience. In addition, manual selection is very time-consuming when processing a large number of SAR images.

[0004] Existing measurement data can also be used to extract high-quality GCPs to improve baseline parameters for orbit refinement and remove interferometric phase errors to generate high-precision InSAR DEM. Although this is a reliable method, not every SAR image has measurement data, and field measurements at each research site are time-consuming, labor-intensive and costly. Artificial corner reflectors can provide the highest quality GCPs for SAR geometric calibration, radiation calibration and InSAR processing, but artificial corner reflectors need to be professionally designed and manufactured for different satellites and placed on site to obtain their accurate coordinates and attitudes. Summary of the invention

[0005] The purpose of the present invention is to overcome the defect that in the process of generating a digital elevation model (DEM) using interferometric synthetic aperture radar (InSAR), ground control points (GCPs) are generally determined by manual selection or field measurement for track refinement and re-leveling, but this cannot ensure that sufficient and suitable GCPs are provided for subsequent interferometric processing, and to provide a method for automatic selection of ground control points based on permanent scatterers.

[0006] The purpose of the present invention can be achieved by the following technical solutions:

[0007] A method for automatically selecting ground control points based on permanent scatterers comprises the following steps:

[0008] Acquire multi-time series SAR images;

[0009] Extract pixels satisfying the permanent scatterer selection strategy from the multi-time series SAR image as ground control points;

[0010] The permanent scatterer selection strategy includes:

[0011] The first selection sub-strategy: calculate the average coherence of each pixel and select the pixel points whose average coherence is greater than a preset low coherence threshold;

[0012] The second selection sub-strategy: select pixels whose amplitude is greater than the minimum amplitude mean;

[0013] The third selection sub-strategy: selecting pixels whose amplitude deviation index is less than a preset amplitude deviation index threshold;

[0014] The fourth selection sub-strategy: selecting pixels whose average coherence is less than a preset high coherence threshold;

[0015] The fifth selection sub-strategy: selecting pixel points whose slope is less than a preset terrain slope threshold.

[0016] Furthermore, the value of the low coherence threshold is within the range of 0.37-0.43.

[0017] Furthermore, the amplitude deviation index threshold value is within the range of 0.15 to 0.25.

[0018] Furthermore, the value of the high coherence threshold is in the range of 0.90 to 0.98.

[0019] Furthermore, the terrain slope threshold is within the range of 13° to 17°.

[0020] Furthermore, the calculation expression of the average coherence is:

[0021]

[0022]

[0023] Where γ is the coherence, M(i,j) is the complex number of the main image, S(i,j) is the complex number of the auxiliary image, the multi-time series SAR image includes the main image and the auxiliary image, i and j are the rows and columns in the square window respectively, the asterisk is the conjugate operator of the complex number, a and b represent the window size, K is the number of SAR images in the multi-time series SAR image, k is the kth image, and γ pix is the average coherence, γ pix-thd is the low coherence threshold.

[0024] Furthermore, the calculation expression of the minimum amplitude mean is:

[0025]

[0026] Where (m,n) represents the total number of rows and columns of SAR images covering the same area, A k(e,f) is the amplitude value of the pixel in row e and column f in the kth SAR image, A img-thd is the minimum value of the amplitude mean, A img is the amplitude value.

[0027] Furthermore, the calculation expression of the amplitude deviation index is:

[0028]

[0029] In the formula, A pix-disp is the amplitude deviation index, A pix-σ is the standard deviation of the amplitude, A pix is the average amplitude of each pixel at the same position in all images, A pix-disp-thd is the amplitude deviation index threshold.

[0030] Furthermore, the method further comprises:

[0031] According to the acquired ground control points, the distance and number between each ground control point are calculated, so as to further extract evenly distributed ground control points, and then the slope of each extracted ground control point is calculated, and the slope is compared with the adjacent ground control points in the preset first area, and the ground control points with lower slopes reaching the preset first number are retained. If the slopes of the ground control points in the first area are the same, they are randomly selected to obtain the final ground control points.

[0032] Furthermore, the selected ground control points are used for track refinement and re-leveling in the process of generating a digital elevation model.

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

[0034] In the process of generating digital elevation models (DEMs) using interferometric synthetic aperture radar (InSAR), ground control points (GCPs) are generally determined by manual selection or field measurement for track refinement and re-leveling. However, this cannot ensure that sufficient and appropriate GCPs are provided for subsequent interferometric processing.

[0035] Therefore, the present invention proposes an automatic selection method of GCPs based on permanent scatterers (PS) for InSARDEM generation. The auto-PS-GCP method selects PS as GCP by automatically defining thresholds of coherence, amplitude, and amplitude deviation index. The auto-PS-GCP-Thin method further performs thinning based on auto-PS-GCP according to point density, distance, and terrain conditions.

[0036] The present invention uses a three-stage evaluation method in the embodiment to evaluate the results, including: 1) evaluating the phase stability and strength of the GCP, 2) the RMSE of the elevation difference between the GCP and the same-name points in the reference DEM, and 3) the accuracy of the generated InSAR DEM. The present invention selects three representative areas of plains, hills and mountains to verify the proposed method. The results show that compared with manual selection, the DEM accuracy obtained by auto-PS-GCP-Thin in different regions is improved by about 20% to 30%, and it performs best in the plains with a DEM accuracy of 4.7051m. In addition, compared with auto-PS-GCP, the DEM accuracy obtained by auto-PS-GCP-Thin is improved by about 20%-40%. It can be concluded that the proposed method is effective and reliable in areas with different terrain conditions. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 A schematic diagram of the geometric relationship of a space-borne InSAR provided in an embodiment of the present invention;

[0038] Figure 2 A flow chart of an automatic selection method of ground control points based on permanent scatterers provided in an embodiment of the present invention and used for track refinement and re-leveling of InSAR DEM generation;

[0039] Figure 3 A GCP thinning flow chart provided in an embodiment of the present invention;

[0040] Figure 4 Schematic diagram of three research areas provided in the embodiments of the present invention. Figure 4 In the figure, (a) DEM map of China, (b) Beijing city center in the plain area, (c) Pingliang city in the hilly area, and (d) Zhangjiakou city in the mountainous area; the base map of (bd) is a composite map of Landsat-8OLI image;

[0041] Figure 5 Schematic diagram of the spatial baseline (ac) and time baseline (df) of the three study areas provided in the embodiment of the present invention; the main image acquisition date of the plain area is January 29, 2020, the hilly area is October 20, 2020, and the mountainous area is September 1, 2020;

[0042] Figure 6 The interference pattern (ac) and the coherence pattern (df) provided in the embodiment of the present invention; the interference pair in the plain area is the main image on January 29, 2020 and the auxiliary image on February 22, 2020, the interference pair in the hilly area is the main image on October 20, 2020 and the auxiliary image on November 13, 2020, and the interference pair in the mountainous area is the main image on September 1, 2020 and the auxiliary image on September 25, 2020;

[0043] Figure 7 The GCP distribution diagram of the manual, auto-PS-GCP and auto-PS-GCP-Thin methods provided in the embodiments of the present invention; the selected GCP is overlaid on the corresponding coherence coefficient diagram, Figure 7 For the three study areas, (ac) indicates manual selection, (df) indicates auto-PS-GCP selection, and (gi) indicates auto-PS-GCP-Thin selection;

[0044] Figure 8 Phase histograms of GCPs obtained manually (ac), auto-PS-GCP (df) and auto-PS-GCP-Thin (gi) provided in the embodiments of the present invention;

[0045] Fig. 9 The DEM generated by the three selected GCPs provided in the embodiment of the present invention is: Fig. 9 Among them, (ac) is the DEM generated by manual point selection, (df) is the DEM generated by auto-PS-GCP point selection, and (gi) is the DEM generated by auto-PS-GCP-Thin point selection;

[0046] Fig.10 This is a schematic diagram of random points provided in an embodiment of the present invention superimposed on an external SRTM DEM. DETAILED DESCRIPTION

[0047] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings here can be arranged and designed in various different configurations.

[0048] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0049] It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, further definition and explanation thereof is not required in subsequent drawings.

[0050] In the description of the present invention, it should be noted that the terms "first" and "second" are used for descriptive purposes only and cannot be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include one or more of the features. In the description of this application, the meaning of "plurality" is two or more, unless otherwise clearly and specifically defined.

[0051] Example 1

[0052] GCPs can also be extracted from PS. PS refers to various ground targets with strong scattering properties, high coherence and stable phase. In SAR images, typical PS usually include the top corners of bridges, buildings and exposed rocks. Such PS have been applied to GB-SAR to generate high-precision DEM. Due to the spatiotemporal stability of PS, they can be used for SBAS-InSAR for sedimentation monitoring. PS technology is also combined with interference processing to extract DEM from TanDEM-X. Since PS has strong scattering characteristics and can maintain phase stability in long time series, this method can produce enough GCPs for track refinement and re-flattening to reduce interferometric phase errors and thus improve DEM accuracy. PS is usually identified and extracted based on coherence, amplitude mean and amplitude deviation index. Some scholars define the coherence threshold for extracting PS, while others use the appropriate amplitude deviation index as the threshold. In fact, a single criterion such as coherence or amplitude mean only considers one aspect of PS. For example, the phase dispersion index threshold only considers phase stability, so the extracted PS may be sensitive to decorrelation and system noise. Coherence considers a high signal-to-noise ratio (SNR). The amplitude dispersion index mainly considers phase stability, which may result in some PS that perform poorly in orbit refinement and re-flattening.

[0053] The selection of PS needs to consider phase stability and scattering intensity, so appropriate thresholds should be defined in advance to extract suitable points and discard unreliable points. And the generation of interferograms is related to the influence of terrain. Valid GCPs are usually distributed in flat areas. Appropriate terrain slope can be used as a metric to help select high-quality GCPs. The distribution and density of the selected GCPs will also significantly affect track refinement and re-flattening in InSAR. Calibrating GCPs by removing points may improve the effect of the application. Therefore, the selection of GCPs should fully consider the influence of terrain conditions and GCP distribution. In addition, there are several stages related to GCP selection in InSAR DEM generation, and each stage should be fully evaluated to ensure its verification. This study aims to solve the key issues of GCP selection: 1) Can this embodiment combine terrain conditions with SAR images to select suitable GCPs for track refinement and re-flattening? 2) How to calibrate the distribution of selected GCPs? 3) How to evaluate GCPs at each stage of InSAR DEM generation?

[0054] In order to solve the above problems, this embodiment provides a new method: an automatic ground control point selection method based on permanent scatterers (auto-PS-GCP), which takes into account the phase, scattering characteristics, coherence, terrain slope and point distribution to extract more reliable and useful GCPs, which are applied to the InSAR process to generate high-precision DEM. This method is applied to the InSAR DEM generation in three test areas: plains (central Beijing), hilly areas (Pingliang) and mountainous areas (Zhangjiakou). Then a three-stage evaluation method is used to evaluate the phase and intensity of the GCP, the accuracy of the GCP itself and the accuracy of the generated DEM.

[0055] This embodiment provides a method for automatically selecting ground control points based on permanent scatterers, comprising the following steps:

[0056] Acquire multi-time series SAR images;

[0057] Extract pixels satisfying the permanent scatterer selection strategy from the multi-time series SAR image as ground control points;

[0058] Permanent scatterer selection strategies include:

[0059] The first selection sub-strategy: calculate the average coherence of each pixel and select the pixel points whose average coherence is greater than a preset low coherence threshold;

[0060] The second selection sub-strategy: select pixels whose amplitude is greater than the minimum amplitude mean;

[0061] The third selection sub-strategy: selecting pixels whose amplitude deviation index is less than a preset amplitude deviation index threshold;

[0062] The fourth selection sub-strategy: selecting pixels whose average coherence is less than a preset high coherence threshold;

[0063] The fifth selection sub-strategy: selecting pixel points whose slope is less than a preset terrain slope threshold.

[0064] The execution order of the first selection sub-strategy, the second selection sub-strategy, the third selection sub-strategy, the fourth selection sub-strategy and the fifth selection sub-strategy is not limited, and each subsequent sub-strategy further extracts pixel points from the selection result of the previous sub-strategy when executed.

[0065] Preferably, the first selection sub-strategy, the second selection sub-strategy, the third selection sub-strategy, the fourth selection sub-strategy and the fifth selection sub-strategy are executed in sequence.

[0066] As a preferred implementation, the low coherence threshold is set within the range of 0.37-0.43 to remove low coherence points.

[0067] As a preferred implementation, the amplitude deviation index threshold is in the range of 0.15 to 0.25 to extract stable GCPs.

[0068] As a preferred implementation, the value of the high coherence threshold is in the range of 0.90 to 0.98 to further select GCPs.

[0069] As a preferred embodiment, the terrain slope threshold is within the range of 13° to 17° to further identify GCPs in flat areas. Incorporating the slope into GCP selection removes points affected by overlapping and shadow effects.

[0070] As a preferred embodiment, the method further comprises:

[0071] According to the acquired ground control points, the distance and number between each ground control point are calculated, so as to further extract evenly distributed ground control points, and then the slope of each extracted ground control point is calculated, and the slope is compared with the adjacent ground control points in the preset first area, and the ground control points with lower slopes reaching the preset first number are retained. If the slopes of the ground control points in the first area are the same, they are randomly selected to obtain the final ground control points.

[0072] As a preferred implementation, the selected ground control points are used for track refinement and re-leveling in the process of generating a digital elevation model.

[0073] Specifically, the implementation process of this embodiment includes the following steps:

[0074] 1. Method execution process

[0075] 1.1 InSAR DEM generation method

[0076] The basic principle of using InSAR technology to extract DEM is to use two SAR antennas with interferometric imaging capabilities (or one antenna for repeated observation) to obtain two SLC images with a certain viewing angle difference in the same area, and invert the elevation information of the surface based on their interferometric phase information to generate the DEM of the corresponding area. There is a phase difference between the main and auxiliary SAR images, which can be combined with the orbital parameters of the observation satellite to generate high-precision and high-resolution ground elevation information. Figure 1 As shown in the figure, due to the different distances between the satellite and the target during the two observations, the phase ( and ) are different. Interference phase It can be given by the following formula:

[0077]

[0078] Where φ1 and φ2 are the phases of the echo signals from the target, λ is the radar wavelength, R1 and R2 are the slant ranges from the two antennas to the observed target P, and ΔR is the slant range difference.

[0079] The interference phase (φint) consists of the terrain phase (φtop), terrain deformation phase (φdef), flat ground phase (φflat), atmospheric phase (φatm) and noise phase (φnoi):

[0080]

[0081] After removing the flat ground effect, atmosphere and noise phase, the real terrain phase of the target can be obtained. based on Figure 1 The geometric relationship shown, the elevation (h) of the target (P) can be given by:

[0082]

[0083] Where B is the SAR antenna A 1 and A 2 The spatial baseline between H and θ is A 1 is the height and viewing angle, and α is the baseline inclination.

[0084] In the process of InSAR DEM generation, stable and high-quality GCPs can effectively improve the performance of track refinement and re-flattening, thereby improving the accuracy of the generated DEM. This example proposes an auto-PS-GCP method for GCP selection. The workflow is as follows: Figure 2 As shown: a) InSAR DEM generation; b) GCP selection using auto-PS-GCP and auto-PS-GCP-thin methods; c) three-stage evaluation. Finally, the auto-PS-GCP and auto-PS-GCP-Thin methods were applied to multi-temporal SAR images of three test areas (located in plains, hills and mountains, respectively) to verify the proposed method. This example evaluates the phase stability and scattering intensity of the selected GCPs, their elevation accuracy, and the accuracy of the generated InSAR DEM according to the three-stage evaluation method.

[0085] 1.2. Auto-PS-GCP selection method

[0086] The auto-PS-GCP selection method selects the appropriate number of PSs by considering both temporal stability and spatial continuity. Figure 2 As shown, GCPs are selected by automatically determining the coherence, amplitude mean and amplitude deviation index thresholds and thinning. To ensure that the amplitudes of all images can truly reflect the scattering of ground targets, this embodiment performs radiation calibration before selecting GCPs.

[0087] 1.2.1 Coherence coefficient threshold

[0088] Coherence (γ) is the first criterion for extracting PS. γ is calculated by the complex number of adjacent pixels in a square window (e.g. 3×3). The calculation method can be given by the following formula:

[0089]

[0090] Where M(i,j) is the complex number of the primary image, S(i,j) is the complex number of the secondary image, these two images form an interference image pair, i and j are the rows and columns in the square window, the asterisk is the conjugate operator of the complex number, and a and b represent the window size (here a=b=3). The SNR of a pixel can be estimated by its coherence:

[0091]

[0092] A higher coherence coefficient indicates a higher SNR, which can produce better interferogram quality. Considering the coherence (γ) of the same pixel in all K-scene images in the time series, each pixel has an average coherence (γ pix ):

[0093]

[0094] Where K represents the number of SAR images and k is the kth image. Select the image with a coherence coefficient greater than the predefined threshold (γ pix-thd ) point as PS.

[0095] A low threshold can only remove severely decoherent points, while a high threshold may miss useful PS. Therefore, this embodiment uses a low coherent threshold (γ pix-thd1 ) Remove areas such as water bodies and vegetation covered areas ( Figure 2 Step 1 in , and then use a high coherence threshold (γ pix-thd2 )Further select PS( Figure 2 4 in the previous step).

[0096] 1.2.2 Amplitude Mean Threshold

[0097] The amplitude mean of each SAR image can represent its scattering intensity. For long-time series of SAR images, there is a minimum value (A img-thd ). This embodiment defines the minimum amplitude mean ( Figure 2 Step 2) in is used as the threshold for extracting PS, and pixels above the threshold are retained as PS:

[0098]

[0099] Where (m,n) represents the total number of rows and columns of SAR images covering the same area, A k(e,f) is the amplitude value of the pixel in row e and column f in the k-th SAR image. img-thd ) can remove the points with low amplitude but artificially high coherence in Step 1.

[0100] 1.2.3. Amplitude Deviation Index Threshold

[0101] Amplitude deviation index is the third criterion for selecting PS ( Figure 2 Step 3 in ). This index can be used to reflect phase stability. Amplitude deviation index (A pix-disp ) can be given by the following formula:

[0102]

[0103] Among them A pix-σ is the standard deviation of the amplitude, A pix is the average amplitude of each pixel at the same position in all images. This embodiment defines a suitable threshold (A pix-disp-thd ) is used to extract points whose amplitude deviation index is less than the threshold as PS.

[0104] This method is not affected by neighboring pixels because it is calculated from pixels of different images at the same location, and this embodiment takes into account the phase stability of each pixel and removes points with severe decoherence from the PS candidate points. This method in formula (8) is suitable for pixels with high SNR, which can be extracted according to formula (4).

[0105] 1.2.4. Terrain slope threshold

[0106] This embodiment uses the terrain slope as one of the criteria for PS selection ( Figure 2 Step 5 in the figure is used to reduce the impact of steep terrain on interferometry accuracy. Based on the PS with good phase stability and strong scattering intensity generated in Step 1-4, points located in the flat area are further extracted. Therefore, a suitable threshold (SLOPE thd ) to extract the slope less than SLOPE thd PS. The slope can be calculated as follows:

[0107]

[0108] Where p and q are elevation gradients, representing the rate of change of elevation in the east-west and north-south directions respectively, and 180 / π is the conversion factor from radians to degrees.

[0109] 1.2.5 GCP Sparseness

[0110] Initial GCPs (i.e., PS) are selected based on coherence, amplitude mean, amplitude deviation index, and slope. When applied to the InSAR process to generate DEM, these GCPs are used to correct baseline parameters and remove flat ground effects. In addition, the number and distribution of GCPs will also affect the removal of baseline parameters and flat ground effects.

[0111] For areas with high-density GCPs obtained using the above method, this embodiment considers the density of points and takes the slope as a priority to optimize the points by eliminating some of them ( Figure 3 ). First, considering the point distance and the number of points, appropriate GCPs are selected from the candidate points to ensure that the points are as evenly distributed as possible. Then, the slope of each GCP is calculated and compared with its adjacent pixels, and the points with lower slopes are retained. When these GCPs have similar slopes, this embodiment randomly selects from the candidate points.

[0112] The average nearest neighbor (ANN) metric is applied to analyze the distribution of GCPs before and after thinning, which is calculated based on the average distance of each GCP to its nearest neighbor GCP. The observed average distance (NNObserved) is calculated from the distance between each GCP and its nearest neighbor GCP. The z-score is a measure of statistical significance, indicating whether the GCPs are randomly distributed. At a 99% confidence level, if the z-score is greater than 2.58, the GCPs are dispersed, and if the z-score is less than -2.58, the GCPs are clustered.

[0113] 1.3. Three-stage evaluation

[0114] This example proposes a three-stage evaluation method to evaluate the selection of GCPs and the generated DEM ( Figure 2 c). In the first stage, the phase and intensity of the selected GCP are evaluated. If the interference phase of the selected GCP is concentrated around 0 rad and the intensity of the GCP is more than 5 times higher than the intensity mean of the image, the selected GCP is considered to be valid. In the second stage, the elevations of the selected GCP and the points of the same name in the reference DEM are compared. This embodiment considers that GCPs with an RMSE of less than 10 m in flat areas and less than 20 m in hilly and mountainous areas are valid in InSAR processing. In the third stage, the external high-precision DEM is compared with the InSAR DEM generated by the auto-PS-GCP and auto-PS-GCP-Thin methods.

[0115] 3. Study area and dataset

[0116] 3.1 Study area and dataset

[0117] This example selects three research areas: plains, hills and mountains ( Figure 4) to test the auto-PS-GCP and auto-PS-GCP-Thin methods. The plain area is the center of Beijing located in the North China Plain ( Figure 4 b), the hilly area is Pingliang City, Gansu Province, located in the Loess Hills ( Figure 4 c), the mountainous area is Zhangjiakou City, Hebei Province, located in Yanshan Mountain ( Figure 4 d). The urban area of ​​Beijing is about 500 square kilometers, and the Pingliang and Zhangjiakou areas are about 600 square kilometers. The vegetation in these three areas is sparse, which reduces the impact of vegetation cover on interferometry.

[0118] This embodiment uses Sentinel-1A SLC data from the dual-polarization C band. Table 1 shows the orbit-raising data of 24 scenes in IW mode in three test areas, with a resolution of 5m×20m and a polarization mode of VV (scihub.copernicus.eu). The accuracy of orbital parameters has a great influence on the results of InSAR processing. The orbital error will exist in the interference pattern in the form of residual fringes. Therefore, this embodiment uses POD Precise Orbit Ephemerides data (slqc.asf.alaska.edu) to correct the orbital parameters.

[0119] This example uses the high-resolution DEM obtained by the PALSAR sensor of the Advanced Land Observing Satellite 1 Project (ALOS) of the Japan Aerospace Exploration Agency (search.asf.alaska.edu) as the reference DEM. The reference DEM has a terrain resolution of 12.5m and a vertical resolution of 12m in Fine Beam Single polarization mode. This example selects SRTM DEM as the external DEM for accuracy analysis. SRTM is equipped with two antennas during the 11-day mission and provides high-resolution DEMs on a nearly global scale for free.

[0120] Table 1 Sentinel-1A SAR dataset

[0121]

[0122]

[0123] 2.2 Data Preprocessing

[0124] The workflow of InSAR DEM generation is as follows Figure 2 This example considers both temporal and spatial baselines and selects the main image ( Figure 5), the other 23 SAR images are auxiliary images, and then 23 interferograms are generated for each study area. For the plain area and mountainous area, the average spatial baseline is 48.5435m and 47.5130m respectively, and the time baseline of the two study areas is suitable. However, for the hilly area, the average spatial baseline is 48.2428m, which is more suitable for generating interferograms. However, the maximum time baseline is 479 days, which may reduce the coherence of the interferogram.

[0125] After the master and slave images are registered, interference processing is performed to generate an interference map, and the interference map is interference filtered using the Goldstein filter method. The filtered interference map is then phase unwrapped using the minimum cost flow method to obtain the true phase corresponding to the terrain information. Then, the selected GCP is imported for track refinement and re-leveling to improve the interference phase. Since the selection of GCP is carried out in the geographic coordinate system, this embodiment converts the GCP into a slant range coordinate system to generate an InSAR DEM. Finally, the DEMs of the three study areas are obtained through phase elevation conversion and geocoding.

[0126] Figure 6 The following are typical interference patterns and coherence patterns obtained from 23 interferometric pairs. Together, they reflect the undulating terrain information. The higher correlation coefficient and clearer interference fringes indicate that the interferometric results are more reliable. Figure 6 a&d) and mountainous areas ( Figure 6 c&f), the interference fringes are clear and most areas have high coherence. The overall coherence of the hilly area is lower than that of the other two areas, but most areas still show clear interference fringes ( Figure 6 b) and high coherence ( Figure 6 e). Therefore, InSAR DEM can be efficiently generated using these interferograms and coherence maps.

[0127] 4. Test and results

[0128] 4.1. GCP selection based on auto-PS-GCP method

[0129] Before geocoding the unwrapped interferometric phase, the selected GCPs are imported into the InSAR process for track refinement and re-flattening to improve the phase. According to the proposed auto-PS-GCP method, this embodiment selects GCPs for each study area by combining the coherence coefficient, amplitude mean, amplitude deviation index, terrain slope and distribution of GCPs, as shown in Table 2. pix-thd1 ) is defined as 0.4, and the amplitude mean threshold (A img-thd ) is the minimum value of the amplitude mean A img , the amplitude deviation index thresholds of the three regions (A pix-disp-thd) are 0.25, 0.15 and 0.20 respectively, and then define the high threshold of the coherence coefficient (γ pix-thd2 ) to select more useful GCPs from the PS candidate points obtained above. Among them, the plain study area is an urban area with a large number of strong scattering objects, so the high threshold of the coherence coefficient (γ pix-thd2 ) is set to 0.98 to get a suitable number of GCPs. Finally, GCPs with slopes below 15° are selected to ensure that these points are located in flat areas. The GCPs derived using auto-PS-GCP show relatively strong clustering, which can be weakened using the rarefaction method.

[0130] Table 2 Evaluation of thresholds defined by the auto-PS-GCP method and the distribution of GCPs

[0131]

[0132] Figure 7 The distribution of GCPs by manual, auto-PS-GCP and auto-PS-GCP-Thin methods. This example uses the manual and auto-PS-GCP-Thin methods to obtain the same number of GCPs, which is convenient for comparing their distribution and the accuracy of the subsequently generated DEM. There are many unchanging scatterers in the center of Beijing, resulting in a large number of GCPs in the plain area ( Figure 7 d). After thinning the GCP of auto-PS-GCP ( Figure 7 df), the points of auto-PS-GCP-Thin show a more reasonable distribution ( Figure 7 gi); specifically, the number of GCPs in plains, hills, and mountains is reduced by about 70%. Table 2 shows that the z-scores of manual and auto-PS-GCP-Thin are greater than 2.58, indicating that their GCPs are dispersed; while the z-score of auto-PS-GCP is less than -2.58, indicating that their GCPs are clustered. The thinning process transforms the automatically selected GCPs from clustered to scattered, thereby improving the GCP distribution.

[0133] 3.2. Phase I Evaluation of GCP

[0134] This example performs interference processing on all SAR image pairs and evaluates the intensity and phase (first stage evaluation). Due to page limitations, this example only presents the results of one image pair for each study area (Table 3). Since the time interval between each image pair is very short, the influence of terrain deformation, weather conditions and other external factors is very small. Table 3 compares the average intensity between the entire SAR image and the GCPs selected by the three methods. The average intensity of the GCPs selected by auto-PS-GCP and auto-PS-GCP-Thin is higher than that of the entire SAR image, and the average intensity of the manually selected GCPs is not much different from the average intensity of the entire image. Taking the plain area as an example, Table 3 shows that the average intensity of the manually selected GCPs is twice that of the entire image, while the intensity of the GCPs of auto-PS-GCP and auto-PS-GCP-Thin is much greater, more than 20 times that of the entire image. Similar situations are found in the average intensity of hilly and mountainous areas, indicating that the GCPs selected by these two methods have strong scattering properties.

[0135] Table 3 Evaluation of average strength of GCPs

[0136]

[0137]

[0138] Figure 8 The phase histogram of the GCPs in the interferometric results is shown. The phases of the manually selected GCPs do not show a consistent distribution in the three study areas, while the phases of the GCPs selected by auto-PS-GCP and auto-PS-GCP-Thin are all around 0 rad. This shows that the GCPs selected by auto-PS-GCP and auto-PS-GCP-Thin have good phase stability, demonstrating the effectiveness and feasibility of the two methods. For the auto-PS-GCP-Thin method, the selected GCPs not only have high scattering and phase stability, but also show a nearly dispersed distribution in the flat area, indicating that these GCPs can be used in the subsequent InSAR DEM generation.

[0139] 3.3. GCP Phase II Evaluation

[0140] The GCPs selected by the three methods were imported into the InSAR process for track refinement and re-flattening to generate DEM. Fig. 9The DEM generated using InSAR and the selected GCPs is shown. This example compares the elevations of the GCPs in the generated DEM with the points of the same name in the reference DEM and calculates their RMSE (Table 4). The GCPs selected by auto-PS-GCP-Thin have the smallest RMSE, showing the best performance; however, the RMSE of the GCPs selected by auto-PS-GCP is greater than that of manual selection, indicating that the aggregation of GCPs may lead to reduced DEM accuracy. Although the number of GCPs selected manually and auto-PS-GCP-Thin is the same, the InSAR DEM generated using GCPs with strong scattering characteristics and dispersed distribution works better. This shows that the proposed auto-PS-GCP-Thin method has good performance.

[0141] Table 4 RMSE of selected GCPs and the same-name points in the reference DEM

[0142]

[0143] 3.4. Third-stage evaluation of generated DEM

[0144] In order to evaluate the accuracy of DEM, this example uses SRTM DEM with a resolution of 30m as the external DEM and generates 100 random points in each study area ( Fig.10 ). The elevations of these random points were extracted from both types of DEMs and compared to analyze the accuracy (Table 5). The DEMs associated with the GCPs generated by auto-PS-GCP-Thin had the highest accuracy in all three study areas. Due to the influence of GCP clustering, the DEMs associated with the GCPs generated by auto-PS-GCP had the lowest accuracy in all three study areas. The comparison of the accuracy before and after thinning showed that the thinned GCPs produced a DEM with higher accuracy. Among the three landform types, the accuracy was the best in the plain area and the accuracy was poor in the hilly area; however, as confirmed by the second stage evaluation, auto-PS-GCP-Thin performed best in all three study areas.

[0145] Table 5 Analysis of the accuracy of elevation differences between generated DEMs and external DEMs

[0146]

[0147] 4. Discussion

[0148] Pixels with stable phase and strong scattering properties are crucial for track refinement and re-flattening in InSAR. This embodiment proposes a PS-based automatic selection method (auto-PS-GCP-Thin) to extract GCPs for track refinement and re-flattening. GCP selection is performed considering coherence, amplitude, amplitude deviation index and terrain slope. In particular, the clustered GCPs are thinned considering point density and slope priority. Three terrain areas, plains, hills and mountains, are selected as the study area, and these selected GCPs are imported into InSAR processing to generate DEM. A three-stage evaluation is used to verify the proposed method and the results obtained. The first stage shows that the GCPs obtained by auto-PS-GCP-Thin have good phase stability and strong intensity, the second stage shows that the RMSE of the GCPs obtained by auto-PS-GCP-Thin for different terrain areas is about 6 to 12m, and the third stage shows that the RMSE of the DEM generated by auto-PS-GCP-Thin is about 4 to 15m in different areas.

[0149] The definition of thresholds is crucial for the selection of GCPs. This embodiment considers five aspects, namely, coherence, amplitude, amplitude deviation index, terrain slope, and GCP distribution. The low threshold of coherence is set to 0.4 to remove low coherence points, and the amplitude mean threshold is set to the minimum amplitude mean to remove points that show false high coherence. The threshold of the amplitude deviation index is usually set to less than 0.25 to ensure that points that can maintain good stability in long time series are obtained. In this study, this embodiment defines the amplitude deviation index threshold as 0.15-0.25 to extract stable GCPs. Then this embodiment defines the high coherence threshold as 0.90-0.98 to further select GCPs. These selected thresholds are suitable for different fields. In addition, the selection of GCPs is also affected by the overlap and shadow effects of SAR imaging, which is greatly related to the terrain slope. Therefore, this embodiment refers to the early terrain classification research and defines the slope threshold as 15° to further identify GCPs in flat areas. Incorporating slope into GCP selection removes points affected by overlap and shadow effects.

[0150] In summary, the selection of GCPs is jointly determined by phase stability, scattering intensity, and GCP distribution. The GCPs selected by auto-PS-GCP have stable phases and high scattering intensity, while the manually selected GCPs show a dispersed (i.e., more appropriate) distribution. According to the second-stage evaluation, for example, in the plain area, the RMSE of the GCPs selected by auto-PS-GCP is 9.2138m, the RMSE of the GCPs selected by manual selection is 8.5438m, and the RMSE of the GCPs selected by auto-PS-GCP-Thin is only 6.5562m. The results show that the GCP distribution is a more important factor affecting the subsequent InSAR DEM generation. The average GCP intensity of auto-PS-GCP-Thin is lower than that of auto-PS-GCP (see Table 3), which may be due to the deletion of some points with high scattering intensity that are concentrated in a very small area. The average scattering intensity of the GCPs obtained by these two methods is still greater than that of the manually selected GCPs.

[0151] It is well known that terrain conditions have a great impact on the accuracy of InSAR DEM, which has been confirmed in this study. The results show that the three methods have the highest accuracy in the plain area, among which auto-PS-GCP-Thin performs the best with an elevation accuracy of 4.7051m. However, the DEM obtained by this method in the hilly area has an RMSE of 14.9749m, which is the lowest accuracy reported among all three regions (see Table 5). In addition, the coherence in the hilly area is the weakest among all three regions (see Figure 6 ), and accordingly, the selected GCPs are the least. Compared with other methods, the auto-PS-GCP-Thin method can produce DEMs with higher accuracy. For example, the RMSE obtained in the case of generating the DEM of Bushehr city using Sentinel-1A data is 18.53m.

[0152] In order to better compare the three methods, this embodiment adopts a three-stage evaluation method to analyze the phase and intensity of GCP, RMSE, and DEM accuracy. For manual selection, although its GCP is available, the manual selection is variable, so the RMSE is different each time. auto-PS-GCP-Thin takes into account phase stability, scattering intensity, terrain slope, and GCP distribution, so the results are consistent and can be applied to a variety of situations. In this study, this embodiment applies auto-PS-GCP-Thin to three areas with different terrains: plains, hills, and mountains. These areas include areas with dense buildings and vegetation coverage, verifying the proposed GCP selection and three-stage evaluation method.

[0153] 5. Conclusion

[0154] High-quality GCPs for track refinement and re-flattening are essential for InSAR DEM generation. GCPs can be extracted from PS candidate points, which are usually ground targets that maintain strong phase stability and scattering intensity in long-time series images. In this study, two GCP selection methods (auto-PS-GCP and auto-PS-GCP-Thin) are proposed and tested in three regions: plains, hills, and mountains. First, low coherence coefficient and amplitude mean are used to obtain PS candidate points. Then, amplitude deviation index and high coherence coefficient are used to extract higher-quality PS from the above PS candidate points. The influence of terrain slope is further considered to remove high-slope points. Finally, a point thinning method is proposed to correct the distribution of GCPs by removing some points that are clustered in a small area. The results show that compared with manual, the DEM accuracy obtained by auto-PS-GCP-Thin is improved by about 20% to 30% in different regions, among which the DEM accuracy in the plain area is 4.7051m, which is the best effect. In addition, compared with auto-PS-GCP, auto-PS-GCP-Thin improves the DEM accuracy by about 20%-40%.

[0155] The GCP selection method of this embodiment comprehensively considers the influencing factors such as phase stability, scattering intensity, terrain slope and GCP distribution. Three case studies show that auto-PS-GCP-Thin is a reliable GCP selection method that can be applied to InSAR DEM generation. Therefore, the proposed method is effective in various areas with different terrain conditions.

[0156] The preferred specific embodiments of the present invention are described in detail above. It should be understood that a person skilled in the art can make many modifications and changes based on the concept of the present invention without creative work. Therefore, any technical solution that can be obtained by a person skilled in the art through logical analysis, reasoning or limited experiments based on the concept of the present invention on the basis of the prior art should be within the scope of protection determined by the claims.

Claims

1. A method for automatically selecting ground control points based on permanent scatterers, characterized in that: The following steps are involved: Acquire multi-time series SAR images; Extract pixels satisfying the permanent scatterer selection strategy from the multi-time series SAR image as ground control points; The permanent scatterer selection strategy includes: The first selection sub-strategy: calculate the average coherence of each pixel and select the pixel points whose average coherence is greater than a preset low coherence threshold; The second selection sub-strategy: select pixels whose amplitude is greater than the minimum amplitude mean; The third selection sub-strategy: selecting pixels whose amplitude deviation index is less than a preset amplitude deviation index threshold; The fourth selection sub-strategy: selecting pixels whose average coherence is less than a preset high coherence threshold; The fifth selection sub-strategy: selecting pixel points whose slope is less than a preset terrain slope threshold; When each subsequent sub-strategy is executed, further pixel extraction is performed on the selection result of the previous sub-strategy; The method further comprises: According to the acquired ground control points, the distance and number between each ground control point are calculated, so as to further extract evenly distributed ground control points, and then the slope of each extracted ground control point is calculated, and the slope is compared with the adjacent ground control points in the preset first area, and the ground control points with lower slopes reaching the preset first number are retained. If the slopes of the ground control points in the first area are the same, they are randomly selected to obtain the final ground control points.

2. The method for automatically selecting ground control points based on permanent scatterers according to claim 1, characterized in that: The value of the low coherence threshold is within the range of 0.37-0.

43.

3. The method for automatically selecting ground control points based on permanent scatterers according to claim 1, characterized in that: The amplitude deviation index threshold value is within the range of 0.15 to 0.

25.

4. The method for automatically selecting ground control points based on permanent scatterers according to claim 1, characterized in that: The value of the high coherence threshold is in the range of 0.90 to 0.

98.

5. The method for automatically selecting ground control points based on permanent scatterers according to claim 1, characterized in that: The terrain slope threshold is within the range of 13° to 17°.

6. The method for automatically selecting ground control points based on permanent scatterers according to claim 1, characterized in that: The calculation expression of the average coherence is: Where γ is the coherence, M(i,j) is the complex number of the main image, S(i,j) is the complex number of the auxiliary image, the multi-time series SAR image includes the main image and the auxiliary image, i and j are the rows and columns in the square window respectively, the asterisk is the conjugate operator of the complex number, a and b represent the window size, K is the number of SAR images in the multi-time series SAR image, k is the kth image, and γ pix is the average coherence, γ pix-thd is the low coherence threshold.

7. The method for automatically selecting ground control points based on permanent scatterers according to claim 1, characterized in that: The calculation expression of the minimum amplitude mean is: Where m and n represent the total number of rows and columns of SAR images covering the same area, and A k(e,f) is the amplitude value of the pixel in row e and column f in the kth SAR image, A img-thd is the minimum value of the amplitude mean, A img is the amplitude value.

8. The method for automatically selecting ground control points based on permanent scatterers according to claim 1, characterized in that: The calculation expression of the amplitude deviation index is: In the formula, A pix-disp is the amplitude deviation index, A pix-σ is the standard deviation of the amplitude, A pix is the average amplitude of each pixel at the same position in all images, A pix-disp-thd is the amplitude deviation index threshold.

9. The method for automatically selecting ground control points based on permanent scatterers according to claim 1, characterized in that: The selected ground control points are used for track refinement and re-leveling in the process of generating a digital elevation model.

Images