Time- and Space-Varying Atmospheric Phase Correction in Discontinuous Ground-Based Synthetic Aperture Radar Deformation Monitoring

Ground-based synthetic aperture radar (GB-SAR) uses active microwave remote-sensing observation mode to achieve two-dimensional deformation measurement and deformation trend extraction, which shows great prospects in the field of deformation monitoring. However, in the process of GB-SAR deformation monitoring, the disturbances caused by atmospheric effect cannot be neglected, and the atmospheric phases will seriously affect the precision of deformation monitoring. In discontinuous GB-SAR deformation monitoring mode, the atmospheric phases are particularly affected by changes of time and space, so the traditional models of atmospheric phase correction are no longer applicable. In this paper, the interferometric phase signal model considering atmospheric phase is first established. Then, the time- and space-varying characteristics of the atmospheric phase are analyzed, and a novel time- and space-varying atmospheric phase correction algorithm, based on coherent scatterers analysis, is proposed. Finally, slope deformation monitoring experiments are carried out to verify the validity and robustness of the proposed algorithm.


Introduction
Natural disasters, such as landslides and mudslides, seriously threaten people's personal and property safety. Ground-based synthetic aperture radar (GB-SAR) can extract the deformation value and the deformation trend of the overall monitoring area, which has become an important technical means for the deformation monitoring of slopes, dams, bridges, etc. [1][2][3][4][5][6][7][8]. Compared with space-borne and airborne SAR, GB-SAR has the advantages of a short time baseline, little spatial baseline error, and flexible observation mode, which make it a powerful supplement to space-borne SAR and airborne SAR in the field of deformation measurement [9,10]. Although the observation distance of GB-SAR is closer than that of space-borne SAR and airborne SAR, the atmospheric disturbances can still bring interferometric phase delay, which will seriously affect the accuracy and reliability of the deformation production [10][11][12][13]. Therefore, the study of the estimation and correction algorithm of the GB-SAR atmospheric phase is important in the process of GB-SAR interferometry.
Atmospheric disturbances, caused by the difference in the atmospheric temperature, pressure, and humidity in the electromagnetic wave propagation, are shown as the delay phases in radar images. The atmospheric delay phase has always been an important factor affecting the development of the

Ground-Based Synthetic Aperture Radar Signal Model and Atmospheric Phase Model
The observation geometry of GB-SAR is shown in Figure 1. High-range resolution can be achieved by transmitting and receiving electromagnetic signals in the range direction, and the synthetic aperture is realized through moving along the horizontal rail. Therefore, GB-SAR can realize two-dimensional resolution and zero-baseline interferometry. The signal of GB-SAR after imaging processing can be expressed as where i is the SAR image series number, R i = (x i ) 2 + (y − y i ) 2 + H 2 is the distance from the radar to the target, |σ i | is the target scattering coefficient, f c is the radar center frequency, and c is the speed of electromagnetic wave. For ideal zero-baseline interferometric data, the interferometric processing and deformation retrieval model can be expressed as follows: where S* is the conjugate of S, ∆R is the displacement between two data acquisition, λ is the signal wavelength, and ∆ϕ in is the interferometric phase.
where i is the SAR image series number, is the distance from the radar to the target, i  is the target scattering coefficient, fc is the radar center frequency, and c is the speed of electromagnetic wave. For ideal zero-baseline interferometric data, the interferometric processing and deformation retrieval model can be expressed as follows:  (2) where S* is the conjugate of S, R  is the displacement between two data acquisition,  is the signal wavelength, and   in is the interferometric phase.
x y z Figure 1. Observation geometry of fround-based synthetic aperture radar (GB-SAR).
However, in GB-SAR field measurement data processing, the interferometric phase includes not only the deformation phase, but also phases such as the atmospheric phase, the noise phase, and the wrapped phase. The GB-SAR interferometric processing is to extract the deformation phase from these phases, and realize the process of deformation retrieval. The interferometric phase model can be expressed as follows: where  defo is the deformation phase,  atmo is the atmospheric phase,  noise is the noise phase,  rail is the rail error phase, and k is the wrapped phase cycle value.
The noise phase  noise can be removed by the phase filtering and high-quality pixel selection [20], which can improve the signal-to-noise ratio of the interferometric phases. The noise phase is caused by equipment re-installment in discontinuous GB-SAR; the correction method is presented in [21], and is not the focus of this paper. In order to simplify the analysis, we selected some pairs of images with a long time span to achieve the discontinuous GB-SAR data simulation in this paper, so the rail error phase  rail can be neglected. The wrapped phase cycle value k can be solved by the phase-unwrapping algorithm, and the related algorithms are sufficient [22][23][24][25][26]. This paper mainly However, in GB-SAR field measurement data processing, the interferometric phase includes not only the deformation phase, but also phases such as the atmospheric phase, the noise phase, and the wrapped phase. The GB-SAR interferometric processing is to extract the deformation phase from these phases, and realize the process of deformation retrieval. The interferometric phase model can be expressed as follows: where ϕ de f o is the deformation phase, ϕ atmo is the atmospheric phase, ϕ noise is the noise phase, ϕ rail is the rail error phase, and k is the wrapped phase cycle value. The noise phase ϕ noise can be removed by the phase filtering and high-quality pixel selection [20], which can improve the signal-to-noise ratio of the interferometric phases. The noise phase is caused by equipment re-installment in discontinuous GB-SAR; the correction method is presented in [21], and is not the focus of this paper. In order to simplify the analysis, we selected some pairs of images with a long time span to achieve the discontinuous GB-SAR data simulation in this paper, so the rail error phase ϕ rail can be neglected. The wrapped phase cycle value k can be solved by the phase-unwrapping algorithm, and the related algorithms are sufficient [22][23][24][25][26]. This paper mainly studies the influence of the atmospheric phase on the GB-SAR interferometric phase and its correction algorithm. Therefore, the interferometric phase model can be simplified as Atmospheric phases of GB-SAR are caused by the difference of atmospheric pressure, temperature, and humidity in different time and space [12]. Therefore, the atmospheric phase model can be deduced by using the atmospheric physical parameters in the deformation monitoring scene. The atmospheric refractive index N depends on the atmospheric pressure P, the temperature T, and the relative humidity H. N can be expressed as N = N dry + N wet (5) where N dry and N wet are the dry refractive index and the wet refractive index, respectively. The specific expressions are expressed as follows: Therefore, the refractive index N is a vector related to the distance vector (r) and time (t), which can be expressed as follows: Then, the atmospheric phase can be obtained by integrating on the path S of radar electromagnetic wave propagation [12].
According to the atmospheric physics model, Noferini built the atmospheric phase linear correction model based on the analysis of the permanent scatterers in [11], which assumes that the atmospheric refractive index was homogeneous in the azimuth direction and that the atmospheric phase had a linear relationship with the observation distance. The linear model is as follows: where K is an unknown constant and r is the distance from radar to observation point. The above model parameters can be calculated by using the phase value of the permanent scatterers in the scene, so the GB-SAR time series data is mandatory. The discontinuous GB-SAR has a long acquisition interval, so the dataset is small and the temporal decorrelation is serious. Therefore, the permanent scatterers analysis method is not suitable for atmospheric phase correction in discontinuous GB-SAR deformation monitoring.

Time-and Space-Varying Atmospheric Phase Correction Algorithm
In this section, a novel time-and space-varying atmospheric phase correction algorithm for discontinuous GB-SAR deformation monitoring is proposed. Discontinuous deformation monitoring is different from that of the continuous mode, and is suitable for the scene with slow deformation and no continuous real-time monitoring. The GB-SAR data in the discontinuous deformation monitoring has the characteristics of little data and a long time baseline, so the frequency interleaved sub-images entropy method is used to extract the CSs as stable control point candidates. Then, according to the local consistency analysis of the atmospheric phases, the candidate points satisfying the requirements are selected as SCPs, which can be used for discontinuous GB-SAR time-and space-varying atmospheric phase correction.

Coherent Scatterers Extraction Using a Single Ground-Based Synthetic Aperture Radar Pair
In the continuous deformation monitoring mode, GB-SAR can obtain time series data, so the SCPs can be extracted by using permanent scatterer (PS) methods [11]. For discontinuous deformation monitoring, little data and temporal decorrelation make the PS methods inapplicable. The proposed time-and space-varying atmospheric phase correction algorithm exploits frequency interleaved sub-images entropy method to extract the CSs from a single GB-SAR pair [27]. The GB-SAR system for deformation monitoring in this paper uses the stepped frequency continuous wave signal. When echo data acquisition is completed, the frequency spectrum of GB-SAR data can be achieved. Then the frequency spectrum can be divided into two sub-band spectrums through the frequency interleaved sampling method. The sub-band spectrums, after imaging processing, can be transformed into sub-images [28,29]. The generation process of frequency interleaved sub-images is shown in Figure 2.

Coherent Scatterers Extraction Using a Single Ground-Based Synthetic Aperture Radar Pair
In the continuous deformation monitoring mode, GB-SAR can obtain time series data, so the SCPs can be extracted by using permanent scatterer (PS) methods [11]. For discontinuous deformation monitoring, little data and temporal decorrelation make the PS methods inapplicable. The proposed time-and space-varying atmospheric phase correction algorithm exploits frequency interleaved sub-images entropy method to extract the CSs from a single GB-SAR pair [27]. The GB-SAR system for deformation monitoring in this paper uses the stepped frequency continuous wave signal. When echo data acquisition is completed, the frequency spectrum of GB-SAR data can be achieved. Then the frequency spectrum can be divided into two sub-band spectrums through the frequency interleaved sampling method. The sub-band spectrums, after imaging processing, can be transformed into sub-images [28,29]. The generation process of frequency interleaved sub-images is shown in Figure 2. Entropy is an alternative approach to estimate coherence in the spectral domain, which can make full use of the sub-image data information by calculating the covariance matrix. The entropy H between sub-images is defined by [18] where (i, j) is the pixel index, 2N is the total number of sub-images, and Pk(i, j) is where λk represents the eigenvalues of the covariance matrix C, and the covariance matrix is given by Entropy is an alternative approach to estimate coherence in the spectral domain, which can make full use of the sub-image data information by calculating the covariance matrix. The entropy H between sub-images is defined by [18] where (i, j) is the pixel index, 2N is the total number of sub-images, and P k (i, j) is where λ k represents the eigenvalues of the covariance matrix C, and the covariance matrix is given by where S represents the sub-image, and E{·} represents the expectation operation. The entropy value is 0 when the pixels of the sub-images are very similar, and 1 when they are completely unrelated. Therefore, in order to extract the CSs, an appropriate threshold must be set. An adaptive threshold can be obtained by averaging the entropy values of the corner reflectors in the monitoring scene. If the entropy of the pixel is lower than the threshold, the pixel will be selected as a CS. If not, the pixel will be excluded. Some research shows that atmospheric phases have correlation in certain spaces, so in order to prevent a false alarm and to further identify the SCPs, local spatial consistency analysis of the atmospheric phase is performed [15].

Local Spatial Consistency Analysis of the Atmospheric Phase
The local spatial consistency analysis of the atmospheric phase can use the method of StaMPS phase analysis. The StaMPS [30] method mainly uses the spatial correlation of the interferometric phase to evaluate the phase stability of selected pixels. Considering that the atmospheric phases are relevant in local space, the StaMPS analysis method can be used to analyze the local spatial consistency of the atmospheric phases. Firstly, the interferometric phase values of all CSs in a circular region with radius L are statistically averaged, which can be expressed as follows: where φ k is the interferometric phase value of the kth CS in the circular region, and M is the number of CSs in the circular region. If the number of sub-images is N, the time correlation coefficient of the kth CS in a region with radius L can be defined as where γ k ∈ [0, 1]. The better the local spatial consistency is, the greater its value. Otherwise, the smaller the value is, the worse the local spatial consistency of the CS is. The threshold value γ th can be obtained through iterative operation [30]. If γ k ≥ γ th , the CS can be considered as a stable control point-otherwise, the CS will be eliminated.

Time-and Space-Varying Atmospheric Phase Estimation and Correction.
After the detection of the SCPs, spatial interpolation is performed to the phase values of the detected SCPs for estimating the atmospheric phase screen (APS) of the whole monitoring scene. The spatial interpolation of a discrete point set is a classical problem in computational geometry. The Delaunay triangulation method is used for the spatial interpolation in this paper. Delaunay triangulation is an optimized spatial structure [31]. There will be no other points in the outer circle of any triangle in the generated triangle net, and the smallest angle of the generated triangle is as large as possible. Therefore, Delaunay triangulation can make the result automatically approach the regular triangle, which is beneficial to improve the interpolation precision. A diagrammatic sketch of Delaunay triangulation is shown in Figure 3.
The realization of Delaunay triangulation first needs the construction of a super triangle [32], which includes all the SCPs and puts their spatial position into the linked list of triangles. The pending interpolation points are inserted in order, and the triangles surrounding the insertion point are extracted from the triangle linked list. Then the common side of the extracted triangles is deleted. Finally, the insertion point and the extracted triangle vertices are connected, and the insertion of a point in the Delaunay triangle linked list is completed. The newly formed triangle can be optimized based on the optimization criterion, and the optimization results are put into the Delaunay triangle linked list. The above steps are performed in loops until all of the pending interpolation points are inserted. The inserted points are interpolated in space according to the cubic method and the phase values of the triangle vertexes where they are located. The realization of Delaunay triangulation first needs the construction of a super triangle [32], which includes all the SCPs and puts their spatial position into the linked list of triangles. The pending interpolation points are inserted in order, and the triangles surrounding the insertion point are extracted from the triangle linked list. Then the common side of the extracted triangles is deleted. Finally, the insertion point and the extracted triangle vertices are connected, and the insertion of a point in the Delaunay triangle linked list is completed. The newly formed triangle can be optimized based on the optimization criterion, and the optimization results are put into the Delaunay triangle linked list. The above steps are performed in loops until all of the pending interpolation points are inserted. The inserted points are interpolated in space according to the cubic method and the phase values of the triangle vertexes where they are located.
However, the SCPs in the monitoring scene often cannot cover the whole scene, so the Delaunay triangulation algorithm can only estimate the atmospheric phases of the regions surrounded by the SCPs, so the APS of the whole scene cannot be obtained. According to Equation (8) in Section 2, the atmospheric phase can be obtained by integrating the refractive index function on the propagation path, so the values of the atmospheric phases change smoothly in the range direction. The relationship between the atmospheric phase and distance can be expressed as follows: where T0 is the fitting order and At is the fitting coefficient. For the atmospheric phase values in the range direction, the polynomial fitting is performed after spatial interpolation. The objective function is to minimize the sum of distance squares from the atmospheric phase points 0 i P to the polynomial curve atmo φ . The objective function of the polynomial fitting is as follows: After polynomial fitting in the range direction, considering the space-varying characteristics of the atmospheric phases-and that polynomial fitting results also hardly completely cover the entire scene-the polynomial fitting in the azimuth direction is also performed. Then the APS of the monitoring scene is completely obtained, and the atmospheric phase can be corrected according to the interferometric phase model in Equation (4).
The flowchart of the proposed GB-SAR time-and space-varying atmospheric phase correction algorithm is shown in Figure 4. However, the SCPs in the monitoring scene often cannot cover the whole scene, so the Delaunay triangulation algorithm can only estimate the atmospheric phases of the regions surrounded by the SCPs, so the APS of the whole scene cannot be obtained. According to Equation (8) in Section 2, the atmospheric phase can be obtained by integrating the refractive index function on the propagation path, so the values of the atmospheric phases change smoothly in the range direction. The relationship between the atmospheric phase and distance can be expressed as follows: where T 0 is the fitting order and A t is the fitting coefficient. For the atmospheric phase values in the range direction, the polynomial fitting is performed after spatial interpolation. The objective function is to minimize the sum of distance squares from the atmospheric phase points P i 0 to the polynomial curve ϕ atmo . The objective function of the polynomial fitting is as follows: After polynomial fitting in the range direction, considering the space-varying characteristics of the atmospheric phases-and that polynomial fitting results also hardly completely cover the entire scene-the polynomial fitting in the azimuth direction is also performed. Then the APS of the monitoring scene is completely obtained, and the atmospheric phase can be corrected according to the interferometric phase model in Equation (4).
The flowchart of the proposed GB-SAR time-and space-varying atmospheric phase correction algorithm is shown in Figure 4.

Experimental Results
In order to verify the effectiveness and robustness of the time-and space-varying atmospheric phase correction algorithm proposed in this paper, two slope deformation monitoring experiments in different scenes and atmospheric conditions were carried out. The experimental sites have different climates and height differences. The first experimental site is located in north China, and is a mine park in the suburb of Beijing. The height difference of the observation scene is small and the climate is dry. The second experimental site is located in a mine in Fujian province, in the southern part of China. This scene has large height difference, and the climate in south China is humid. All of the experimental conditions are beneficial for the validation of the proposed time-and space-varying atmospheric phase correction algorithm.
The GB-SAR system used in the experiments is mainly composed of a horizontal rail, a radar transmitting/receiving (T/R) module, and a hardware and software control cabinet. The system uses a stepped frequency continuous wave signal. The maximum observation distance of the system is 5 km. The system parameters are shown in Table 1.

Experimental Results
In order to verify the effectiveness and robustness of the time-and space-varying atmospheric phase correction algorithm proposed in this paper, two slope deformation monitoring experiments in different scenes and atmospheric conditions were carried out. The experimental sites have different climates and height differences. The first experimental site is located in north China, and is a mine park in the suburb of Beijing. The height difference of the observation scene is small and the climate is dry. The second experimental site is located in a mine in Fujian province, in the southern part of China. This scene has large height difference, and the climate in south China is humid. All of the experimental conditions are beneficial for the validation of the proposed time-and space-varying atmospheric phase correction algorithm.
The GB-SAR system used in the experiments is mainly composed of a horizontal rail, a radar transmitting/receiving (T/R) module, and a hardware and software control cabinet. The system uses a stepped frequency continuous wave signal. The maximum observation distance of the system is 5 km. The system parameters are shown in Table 1.

First Case Study: Field Measurements in Northern China
The first experimental scene is shown in Figure 5. Figure 5a is the optical photo of the monitoring scene, and Figure 5b is the aerial view of the monitoring scene on Google Earth. The right area of the monitoring scene is prone to deformation, which is marked using a red color in Figure 5a. The left area of the monitoring scene is a stable area. The monitoring scene, combining the deformation area and the stable area, is very beneficial to GB-SAR deformation monitoring and data analysis. Two corner reflectors are arranged in the scene, which are used for the calibration and the verification of the later data processing. The arrangement of the corner reflectors is difficult in the deformation area, so all of the corner reflectors are arranged in the stable area.

First Case Study: Field Measurements in Northern China
The first experimental scene is shown in Figure 5. Figure 5a is the optical photo of the monitoring scene, and Figure 5b is the aerial view of the monitoring scene on Google Earth. The right area of the monitoring scene is prone to deformation, which is marked using a red color in Figure 5a. The left area of the monitoring scene is a stable area. The monitoring scene, combining the deformation area and the stable area, is very beneficial to GB-SAR deformation monitoring and data analysis. Two corner reflectors are arranged in the scene, which are used for the calibration and the verification of the later data processing. The arrangement of the corner reflectors is difficult in the deformation area, so all of the corner reflectors are arranged in the stable area. The deformation monitoring experiment lasted for one month and the slope was monitored in real time, which was the continuous deformation measurement mode. However, the monitoring scene is a slow deformation scene for the GB-SAR, and the high-frequency data acquisition will increase unnecessary costs. Therefore, the discontinuous deformation monitoring mode is a cost-effective solution. We focus on the discontinuous deformation monitoring data in this case. With respect to the slow deformation scene, the time baseline is increased to achieve the discontinuous deformation monitoring of the scene. A long time baseline will trigger the time-varying characteristics of the environmental parameters. Moreover, the atmospheric refractive indexes are heterogeneous in the azimuth, so the atmospheric phases are space-varying. Therefore, the time-and space-varying atmospheric phase correction algorithm proposed in Section 3 is used for the atmospheric phase correction of the discontinuous deformation monitoring data. Finally, the deformation phases can be obtained. The experimental results are demonstrated in Figure 6.  The deformation monitoring experiment lasted for one month and the slope was monitored in real time, which was the continuous deformation measurement mode. However, the monitoring scene is a slow deformation scene for the GB-SAR, and the high-frequency data acquisition will increase unnecessary costs. Therefore, the discontinuous deformation monitoring mode is a cost-effective solution. We focus on the discontinuous deformation monitoring data in this case. With respect to the slow deformation scene, the time baseline is increased to achieve the discontinuous deformation monitoring of the scene. A long time baseline will trigger the time-varying characteristics of the environmental parameters. Moreover, the atmospheric refractive indexes are heterogeneous in the azimuth, so the atmospheric phases are space-varying. Therefore, the time-and space-varying atmospheric phase correction algorithm proposed in Section 3 is used for the atmospheric phase correction of the discontinuous deformation monitoring data. Finally, the deformation phases can be obtained. The experimental results are demonstrated in Figure 6.

First Case Study: Field Measurements in Northern China
The first experimental scene is shown in Figure 5. Figure 5a is the optical photo of the monitoring scene, and Figure 5b is the aerial view of the monitoring scene on Google Earth. The right area of the monitoring scene is prone to deformation, which is marked using a red color in Figure 5a. The left area of the monitoring scene is a stable area. The monitoring scene, combining the deformation area and the stable area, is very beneficial to GB-SAR deformation monitoring and data analysis. Two corner reflectors are arranged in the scene, which are used for the calibration and the verification of the later data processing. The arrangement of the corner reflectors is difficult in the deformation area, so all of the corner reflectors are arranged in the stable area. The deformation monitoring experiment lasted for one month and the slope was monitored in real time, which was the continuous deformation measurement mode. However, the monitoring scene is a slow deformation scene for the GB-SAR, and the high-frequency data acquisition will increase unnecessary costs. Therefore, the discontinuous deformation monitoring mode is a cost-effective solution. We focus on the discontinuous deformation monitoring data in this case. With respect to the slow deformation scene, the time baseline is increased to achieve the discontinuous deformation monitoring of the scene. A long time baseline will trigger the time-varying characteristics of the environmental parameters. Moreover, the atmospheric refractive indexes are heterogeneous in the azimuth, so the atmospheric phases are space-varying. Therefore, the time-and space-varying atmospheric phase correction algorithm proposed in Section 3 is used for the atmospheric phase correction of the discontinuous deformation monitoring data. Finally, the deformation phases can be obtained. The experimental results are demonstrated in Figure 6.  In the continuous deformation monitoring mode, the data acquisition time of one GB-SAR image is 20 min. In order to obtain the discontinuous deformation data in this experiment, the interferometric data with three-day intervals are selected for processing. Figure 6a shows the reflectivity image of the monitoring scene. The slope area in the monitoring scene has achieved high imaging performance. The correlation coefficient of the GB-SAR data is shown in Figure 6b. The correlation coefficient is a parameter representing correlation of the interferometric pixels [33]. The larger the value of the correlation coefficient, the better the consistency and the phase quality is. Figure 6c is the sub-images entropy map of the GB-SAR interferometric data. The performance of the sub-images entropy is consistent with the correlation coefficient. The pixels in the slope area show good correlation and entropy, and the pixels in the other area are noise, which have poor correlation and entropy. The interferometric phases without atmospheric phase correction are shown in Figure 6d. It can be clearly seen that the atmospheric phases distributed in the scene need to be corrected. The unstable region on the right side of the slope has obvious deformation, which is also consistent with the marked area in Figure 5a.
In order to verify the robustness of the proposed algorithm, the atmospheric phase correction algorithm, based on a linear model, is performed for comparison. Figure 6e is the atmospheric phase correction result using the linear model presented in Section 2 [11]. Figure 6f is the result using the atmospheric phase correction algorithm proposed in this paper. It can be qualitatively seen that the correction method based on the linear model has achieved a certain correction effect. However, the atmospheric parameters of the scene in discontinuous GB-SAR deformation monitoring mode are time-varying, and the atmospheric phases have spatial heterogeneity. Therefore, the corrected result using the linear model shows that the atmospheric phases are not corrected completely and the deformation phases on the left side of the slope still have more residual atmospheric phases. The atmospheric phase correction algorithm proposed in this paper extracts CSs through the frequency interleaved sub-images entropy method. Then the local spatial consistency analysis is performed, and the APS can be obtained by the spatial interpolation, based on Delaunay triangulation and the polynomial fitting. Our proposed algorithm not only considers the heterogeneity of the atmospheric phase in azimuth, but also completes the atmospheric phase estimation and correction through detecting enough stable CSs. Therefore, it can achieve better atmospheric phase correction performance than the algorithm using the linear atmospheric phase correction model. A section of the stable area is used for the evaluation of the atmospheric phase residuals. The histograms of the atmospheric phase residual, using the two methods, are shown in Figure 7. The residual phases, using the linear model method as shown in Figure 7a, are mainly distributed in (0.15-0.45 rad). The result of our proposed method, presented in Figure 7b, has smaller residual phases that are mainly distributed in (0-0.1 rad), and the atmospheric phase correction effect is better, which is consistent with our previous conclusion. In the continuous deformation monitoring mode, the data acquisition time of one GB-SAR image is 20 min. In order to obtain the discontinuous deformation data in this experiment, the interferometric data with three-day intervals are selected for processing. Figure 6a shows the reflectivity image of the monitoring scene. The slope area in the monitoring scene has achieved high imaging performance. The correlation coefficient of the GB-SAR data is shown in Figure 6b. The correlation coefficient is a parameter representing correlation of the interferometric pixels [33]. The larger the value of the correlation coefficient, the better the consistency and the phase quality is. Figure 6c is the sub-images entropy map of the GB-SAR interferometric data. The performance of the sub-images entropy is consistent with the correlation coefficient. The pixels in the slope area show good correlation and entropy, and the pixels in the other area are noise, which have poor correlation and entropy. The interferometric phases without atmospheric phase correction are shown in Figure 6d. It can be clearly seen that the atmospheric phases distributed in the scene need to be corrected. The unstable region on the right side of the slope has obvious deformation, which is also consistent with the marked area in Figure 5a.
In order to verify the robustness of the proposed algorithm, the atmospheric phase correction algorithm, based on a linear model, is performed for comparison. Figure 6e is the atmospheric phase correction result using the linear model presented in Section 2 [11]. Figure 6f is the result using the atmospheric phase correction algorithm proposed in this paper. It can be qualitatively seen that the correction method based on the linear model has achieved a certain correction effect. However, the atmospheric parameters of the scene in discontinuous GB-SAR deformation monitoring mode are time-varying, and the atmospheric phases have spatial heterogeneity. Therefore, the corrected result using the linear model shows that the atmospheric phases are not corrected completely and the deformation phases on the left side of the slope still have more residual atmospheric phases. The atmospheric phase correction algorithm proposed in this paper extracts CSs through the frequency interleaved sub-images entropy method. Then the local spatial consistency analysis is performed, and the APS can be obtained by the spatial interpolation, based on Delaunay triangulation and the polynomial fitting. Our proposed algorithm not only considers the heterogeneity of the atmospheric phase in azimuth, but also completes the atmospheric phase estimation and correction through detecting enough stable CSs. Therefore, it can achieve better atmospheric phase correction performance than the algorithm using the linear atmospheric phase correction model. A section of the stable area is used for the evaluation of the atmospheric phase residuals. The histograms of the atmospheric phase residual, using the two methods, are shown in Figure 7. The residual phases, using the linear model method as shown in Figure 7a, are mainly distributed in (0.15-0.45 rad). The result of our proposed method, presented in Figure 7b, has smaller residual phases that are mainly distributed in (0-0.1 rad), and the atmospheric phase correction effect is better, which is consistent with our previous conclusion.

Second Case Study: Field Measurements in Southern China
The monitoring scene of the second case is shown in Figure 8a. The height difference between the monitoring slope and the position of the GB-SAR system is large. The large height difference and the humid climate in southern China make the second monitoring experiment have great significance on the study of atmospheric phase correction. In the first case study, we only conducted the qualitative analysis on the experimental results for the effectiveness validation of the proposed algorithm. In order to further quantitatively analyze the correction accuracy, corner reflectors were set in the monitoring scene, which is shown in Figure 8b. The corner reflectors used in the experiments are large and heavy, and we took strengthening measures to ensure enough stability. The total station instrument was used to measure the deformation data of the corner reflectors, which was essential for the comparison and verification of the GB-SAR results.
The second experiment lasted 20 days, which realized the real-time deformation monitoring of the slope. The interferometric data with two-day intervals were selected as the discontinuous deformation monitoring data. Due to the great distance from the slope to the location of the GB-SAR system, the closest area is not the region of interest. Therefore, the monitoring range is limited to (400-1400 m) by software settings. The data processing results of the experiment are shown in Figure 9.

Second Case Study: Field Measurements in Southern China
The monitoring scene of the second case is shown in Figure 8a. The height difference between the monitoring slope and the position of the GB-SAR system is large. The large height difference and the humid climate in southern China make the second monitoring experiment have great significance on the study of atmospheric phase correction. In the first case study, we only conducted the qualitative analysis on the experimental results for the effectiveness validation of the proposed algorithm. In order to further quantitatively analyze the correction accuracy, corner reflectors were set in the monitoring scene, which is shown in Figure 8b. The corner reflectors used in the experiments are large and heavy, and we took strengthening measures to ensure enough stability. The total station instrument was used to measure the deformation data of the corner reflectors, which was essential for the comparison and verification of the GB-SAR results.

Second Case Study: Field Measurements in Southern China
The monitoring scene of the second case is shown in Figure 8a. The height difference between the monitoring slope and the position of the GB-SAR system is large. The large height difference and the humid climate in southern China make the second monitoring experiment have great significance on the study of atmospheric phase correction. In the first case study, we only conducted the qualitative analysis on the experimental results for the effectiveness validation of the proposed algorithm. In order to further quantitatively analyze the correction accuracy, corner reflectors were set in the monitoring scene, which is shown in Figure 8b. The corner reflectors used in the experiments are large and heavy, and we took strengthening measures to ensure enough stability. The total station instrument was used to measure the deformation data of the corner reflectors, which was essential for the comparison and verification of the GB-SAR results.
The second experiment lasted 20 days, which realized the real-time deformation monitoring of the slope. The interferometric data with two-day intervals were selected as the discontinuous deformation monitoring data. Due to the great distance from the slope to the location of the GB-SAR system, the closest area is not the region of interest. Therefore, the monitoring range is limited to (400-1400 m) by software settings. The data processing results of the experiment are shown in Figure 9.  The second experiment lasted 20 days, which realized the real-time deformation monitoring of the slope. The interferometric data with two-day intervals were selected as the discontinuous deformation monitoring data. Due to the great distance from the slope to the location of the GB-SAR system, the closest area is not the region of interest. Therefore, the monitoring range is limited to (400-1400 m) by software settings. The data processing results of the experiment are shown in Figure 9.  Figure 9a shows the radar reflectivity image of the multi-stage slope in scene 2, and it can be seen that the slope stages have achieved good imaging performance. Figure 9b is the obtained frequency interleaved sub-images entropy map, which shows that the entropy values of the slope in the monitoring scene are small. Figure 9c is the interferometric phase before the atmospheric phase correction. It can be seen that serious atmospheric phases are shown in the uncorrected interferometric phase image, which is consistent with the humid climate in southern China and seriously affects the accuracy of deformation monitoring. Figure 9d shows the interferometric phase after the linear atmospheric phase correction, and the correction result is not satisfactory. Figure 9e shows the interferometric phase, using our proposed phase correction method, which achieved good atmospheric phase correction effect. The deformation map corresponding to our corrected phase is shown in Figure 9f. In order to verify the precision of the proposed atmospheric phase correction algorithm, quantitative analysis results are given in Table 2. The deformation results of the corner reflectors in the scene (1) with linear atmospheric phase correction; (2) with our proposed atmospheric phase correction; (3) without the atmospheric phase correction; and (4) the total station deformation measurement results are presented. The observation position of the total station is located at the center of the GB-SAR, and the two corner reflectors are arranged along the radar's line of sight. As the displacement direction is the radar's line of sight, the displacement is calculated by projecting the displacement of the total station into the direction of radar's line of sight. It can be seen that the deformation values of the corner reflectors after our atmospheric phase correction are closer to the measurement results of the total station than the results before correction and of the linear model correction. Actually, the corner reflectors are SCPs, so the theoretical deformation value during the period of the deformation monitoring should be 0, which is consistent with the deformation results after atmospheric phase correction. The above results quantitatively and qualitatively validate that the proposed algorithm can improve the accuracy and reliability of deformation monitoring.   Figure 9a shows the radar reflectivity image of the multi-stage slope in scene 2, and it can be seen that the slope stages have achieved good imaging performance. Figure 9b is the obtained frequency interleaved sub-images entropy map, which shows that the entropy values of the slope in the monitoring scene are small. Figure 9c is the interferometric phase before the atmospheric phase correction. It can be seen that serious atmospheric phases are shown in the uncorrected interferometric phase image, which is consistent with the humid climate in southern China and seriously affects the accuracy of deformation monitoring. Figure 9d shows the interferometric phase after the linear atmospheric phase correction, and the correction result is not satisfactory. Figure 9e shows the interferometric phase, using our proposed phase correction method, which achieved good atmospheric phase correction effect. The deformation map corresponding to our corrected phase is shown in Figure 9f. In order to verify the precision of the proposed atmospheric phase correction algorithm, quantitative analysis results are given in Table 2. The deformation results of the corner reflectors in the scene (1) with linear atmospheric phase correction; (2) with our proposed atmospheric phase correction; (3) without the atmospheric phase correction; and (4) the total station deformation measurement results are presented. The observation position of the total station is located at the center of the GB-SAR, and the two corner reflectors are arranged along the radar's line of sight. As the displacement direction is the radar's line of sight, the displacement is calculated by projecting the displacement of the total station into the direction of radar's line of sight. It can be seen that the deformation values of the corner reflectors after our atmospheric phase correction are closer to the measurement results of the total station than the results before correction and of the linear model correction. Actually, the corner reflectors are SCPs, so the theoretical deformation value during the period of the deformation monitoring should be 0, which is consistent with the deformation results after atmospheric phase correction. The above results quantitatively and qualitatively validate that the proposed algorithm can improve the accuracy and reliability of deformation monitoring.

Conclusions
Atmospheric phase correction plays an important part in GB-SAR deformation monitoring, which is directly related to the precision and reliability of deformation production. In the discontinuous GB-SAR deformation monitoring, the atmospheric phases have time-and space-varying characteristics, and the conventional linear atmospheric phase correction model is no longer applicable. Therefore, in this paper, a novel time-and space-varying atmospheric phase correction algorithm based on CSs analysis is proposed. First of all, the frequency interleaved sub-images entropy method is used to detect the CSs in the scene. Then the extraction of SCPs is completed through the spatial consistency analysis. Finally, the spatial interpolation based on Delaunay triangulation and the polynomial fitting are used to realize the correction of the atmospheric phases. The proposed algorithm can use only a single GB-SAR pair to correct the atmospheric phases in the discontinuous deformation monitoring mode. Two case study results are given for the qualitative and quantitative analysis. Through analyzing the experimental results, and comparing those results with those of the linear atmospheric phase correction model and total station, the validity and reliability of the proposed algorithm are verified.

Conflicts of Interest:
The authors declare no conflict of interest.