A Compensation Method for a Time–Space Variant Atmospheric Phase Applied to Time-Series GB-SAR Images

: An atmospheric e ﬀ ect is a main error source that a ﬀ ects interferometric measurements. When a ground-based multiple-input multiple-output (GB-MIMO) radar, i.e., a speciﬁc type of GB-synthetic aperture radar (GB-SAR), was utilized to continuously monitor an open-pit mine, the interferometric phases of some interferograms were complexly space-variant due to time-variant weather conditions. The conventional method of atmospheric phase (AP) compensation was no longer applicable. This paper proposes an improved compensation method of a time-space variant AP applied to time-series GB-SAR images. The permanent scatterers (PSs) were classiﬁed into three types based on their di ﬀ erent spatial properties: The noise-dominant PS (NPS), the deformation-dominant PS (DPS), and the atmospheric e ﬀ ect-dominant PS (APS). The NPSs were ﬁrstly rejected based on the di ﬀ erential phase analysis of neighboring PSs. The DPSs were then rejected based on the cluster partition and selection. With the APSs, the space-variant AP was estimated with a spatial interpolation. To validate the feasibility of the proposed method, short-term and long-term experimental datasets were processed. Comparisons with a conventional method proved that the proposed method can well reduce AP errors and avoid the misunderstanding of motional areas.


Introduction
As a developed type of high-accuracy deformation monitoring system, ground-based synthetic aperture radar (GB-SAR) has been widely utilized to monitor the surface deformations of natural or geological slopes [1,2]. With two radar images acquired at different times and from the same position, their differential phases can be explored to measure 1D deformations based on the differential interferometry technique [3,4]. Since atmospheric conditions (temperature, humidity, pressure etc.) change over time and space, changes which can affect the propagation velocity of a radar wave through the troposphere, a transmission delay is caused due to atmospheric effects. Atmospheric disturbance is one major error source that affects measurement accuracy [5].
Different methods have been proposed to compensate for the atmospheric phase (AP). A common method is to build a meteorological station to measure some important weather parameters. Based on an empirical model of the refractivity index, the AP can be mathematically calculated [6]. The second method is to manually set or specially select some control points with high phase stability. Then, a spatial interpolation can be utilized to compensate for the AP of the whole image [7]. The third method is based on the permanent scatterer (PS) technique, where a parametric model is built according to the spatial characteristics of the atmosphere, and this technique used to estimate the atmospheric parameters [8].
When the PS technique is utilized to compensate for the AP, a proper parametric model should be firstly adopted to simulate the AP [9,10]. In general, the assumption of spatial homogeneity for the AP could work perfectly when the AP model is built as a linear phase ramp over the range axis. The PSs located in the motionless area could be selected to estimate the model parameters. If the motionless area is unknown, those PSs with a high coherence could be employed for a rough estimation. Then, those PSs with a large deviation from the model can be filtered out. A better estimation can be achieved by carrying out a second estimation. This method utilizes a large number of PSs to estimate the atmospheric parameters, and its estimation accuracy can be well guaranteed.
However, when a GB-MIMO (ground-based multiple-input multiple-output) radar was utilized to monitor an open pit for about 12 days, experimental datasets showed that the interferometric phases of some interferograms were complexly space-variant, which was mainly caused by two aspects. Firstly, the AP might have inhomogeneously changed due to bad weather conditions such as rainfall, snowfall, and wind. Secondly, a large motional area could have occurred and caused a nonlinear phase variation along the range axis. Moreover, for some interferograms acquired with a same master image and different slave images within a short time-interval, their phases also changed significantly. The AP is time-space variant for time-series interferograms. The applicability of the conventional compensation method based on the PS technique is doubtful.
This paper introduces an improved method to compensate for the time-space variant AP applied to time-series GB-SAR images. Taking advantage of the different spatial properties between the AP, the DP (deformation phase), and the noise phase (NP), the PSs were classified into three types: Noise-dominant PS (NPS), deformation-dominant PS (DPS) and atmospheric effect-dominant PS (APS). Firstly, a triangulation network was constructed to connect the neighboring PSs. The NPSs were rejected by applying a phase deviation analysis to each pair of neighboring PSs. Then, the remaining PSs were divided into a certain number of clusters, and those clusters located in the motional areas were selected to reject the DPSs. With the residual APSs, a number of control points was determined and utilized to estimate the space-variant AP of all the PSs with a spatial interpolation. With the improved method, short-term and long-term experimental datasets acquired with a GB-MIMO system were processed. Comparisons with the conventional method validated the feasibility and effectiveness of the proposed method.  Figure 1a shows a photo of the experimental scene and Figure 1b shows its top view taken from the Google satellite map. Apparently, its upper side has an approximate elliptical shape, and the lengths of its long and short axes are, respectively, about 1050 and 680 m. The maximum mining depth of the open pit is about 400 m, and the slope angles are 38 • -47 • for most of the regions. A yellow rectangular box shows the radar position, which is located right above the pit and directly opposite the monitoring area. A dotted ellipse shows a construction area where a large deformation has occurred during the experimental period.  Figure 2 shows the photo of the GB-MIMO radar, which utilizes 16 transmitting antennas to constitute two dense transmitting sub-arrays and 16 receiving antennas to constitute one sparse receiving array. An equivalent large aperture with 256 azimuth sampling points can be synthesized. The transmitted signal is a frequency-modulated continuous wave, and the range resolution was obtained using the de-ramping compression technique. The system uses a Ku band with a wavelength of 1.86 cm. Its synthetic aperture is 1.138 m, and the azimuth angular resolution is 0.466 • . The accuracy of the deformation measurement can reach the sub-millimeter level. The radar system was independently implemented by the Beijing Institute of Technology. Generic GB-SAR systems have a large aperture based on the mechanical movement of antennas along a rail track, whereas the GB-MIMO system utilizes the waveform diversity technique to synthesize a large aperture with a multiple antenna structure [11,12]. Their largest difference is the way they synthesize a large aperture. The GB-MIMO radar has the advantages of fast image acquisition and better repeat-pass ability. Other than those, it is an equivalent measuring tool with the traditional GB-SAR and can be regarded as a specific type of GB-SAR. They both face the same compensation problem of a time-space variant AP. Therefore, the proposed method was developed based on GB-MIMO radar and is universal for its applicability of GB-SAR. Figure 3a shows a power image of the monitoring area at the polar coordinate, which was dB processed with respect to the average noise level and could be regarded as a signal-to-noise ratio (SNR) map [13]. Since the beam width in the elevation direction is limited, the bottom of the pit could not be illuminated, which caused those pixels located in the bottom left side of the image to be of a power smaller than −10 dB. For those pixels located in the upper right side whose power was between −10 and 10 dB, they were covered by sparse vegetation or shadowed by terrain. For most pixels located on the pit slope, their power was larger than 10 dB. Based on the ADI (amplitude dispersion index) method, the PS selection was taken with the first 30 radar images [14]. As shown in Figure 3b, the PSs with an ADI lower than 0.1 made up about 82.6%, and those with an ADI lower than 0.15 about 95.9%. Moreover, most PSs were located right on the pit slope. For pixels located on the road, they were not selected as PSs due to the observation angle. The experiment lasted for about 11 days and 3 h, and 5070 radar images were acquired and processed on site. To satisfy the demand of real-time deformation measuring, a specific processing framework for time-series images was adopted. More detailed information can be found in a published article by Hu et al. [15]. The average processing period for every image was about 3.16 min, which included the image acquisition time and real-time processing time. However, on-site phase maps sometimes presented complex variation and could not be well-compensated for, as discussed in Section 3. As can be seen in the long-term deformation map, it seemed that almost the whole pit slope had motioned up to several millimeters. This measurement result was inconsistent with the on site investigation of the mine management.

Experimental Information
The 5070 radar images were then post-processed with a simplified framework, as shown in Figure 4. All N images were divided into an interferometric group (InG). For each InG, a PS set was selected with all the N images. Taking the k th InG as an example, its image sequences were from (k − 1)N + 1 to kN. With the last image (k − 1)N of the (k − 1) th InG being a master image and N images of the k th InG being slave images, N complex interferograms could be constructed. Though the PS selection was taken with just N slave images and without their common master image, its effect could be negligible due to the highly temporal stability of the PS set. In this paper, the proposed method of atmospheric phase compensation (APC) was firstly validated with the framework in Figure 4. As for the real-time processing, every time a new image was acquired together with (N -1) preceding images, they constituted a new InG. Then, the proposed method could be utilized to compensate for the AP for the newly acquired image. Based on the framework shown in Figure 4, every group of N = 30 images was designated s one InG. The time duration of one InG, which was also the temporal baseline of its 30th interferogram, was about 1 h and 40 min; 169 InGs could be acquired with the 5070 images. For each InG, the same ADI threshold of 0.2 and the same dB threshold of 10 dB were utilized for PS selection.

Method Principle
Before taking APC to a phase interferogram, it should be phase unwrapped, where a method based on network flow programming for sparse data could be adopted. The unwrapped phase of a PS could be modeled as: where ϕ defo is the phase component related to the deformation, ϕ atm is the phase component due to the atmospheric effects during image acquisitions, and ϕ noi is the noise term [16,17]. For every point in 3D space and for every moment in time, the electromagnetic properties concerning radio wave propagation are governed by the dimensionless reflectivity N → r , t . Therefore, the AP term in Equation (1) could be described as: where λ is the wavelength. ϕ atm is obtained from the integration of the refractivity variation ∆N along the path L, where ∆N changes with time t and space → r . The atmospheric effect can be distinguished as two types: Turbulent mixing and vertical stratification. Turbulent mixing is a result of different tropospheric processes such as solar heating, differences in wind direction or velocity, fractional drag, and large-scale weather systems. It causes 3D heterogeneity in refractivity during two SAR acquisitions and affects both flat and mountainous terrain. Vertical stratification is a result of different vertical refractivity profiles during two SAR acquisitions, assuming there are no heterogeneities within the horizontal layers. It affects mountainous terrain only, and it is correlated with topography [18,19].
The predominant part of the atmospheric effect is caused by water vapor distribution in the low troposphere. Moreover, its statistics are correlated with the wind vortices [20]. Under the condition of good weather, ∆N could be assumed to be spatially homogeneous and temporally stochastic. Then, ϕ atm can be modeled as a linear phase ramp along the range R: where β o (t) is a constant and β 1 (t) is a linear coefficient, and both are variant for different phase interferograms. When the conventional method is utilized, a system of linear equations could be built: ∆Φ is a n × 1 vector containing the unwrapped phase of n PSs, X is a n × 2 matrix with the observation variables, and ε is a n × 1 vector containing random errors. ∆ϕ i and R i are, respectively, the phase and slant range of the i th PS. β is a 2 × 1 vector to be estimated, where β 0 is a constant and β 1 is the linear coefficient. The unknown vector β can be estimated with a least square regression: where T denotes the matrix transposition. The estimated components ∆Φ AP of the AP are given by: The difference between ∆Φ and ∆Φ AP is the compensated phase. Firstly, all the n PSs are utilized to estimate the atmospheric parameters. Then, to achieve a better estimation, some unreliable PSs with a large deviation from the model should be filtered out, where the judging condition is shown in Equation (7), and the value range of ∆T Atm should be set to 0.1-0.2 rad. With the remaining PSs, a second estimation can be made:

Conventional APC
On each of experimental days of 2 April 2018 and 4 April 2018, a typical InG of 30 radar images was chosen for APC, and these InGs were defined as Group A and Group B. Figure 5a The unwrapped phases of all the PSs in the four interferograms is between -1 and 2 rad, where their phases varied linearly along the range axis. Moreover, the linear slope increased with the temporal baseline. Therefore, during the acquisition period of Group A, the weather condition was stable and the assumption of spatial homogeneity for the AP could work perfectly.  B . Apparently, the unwrapped phases of the PSs in the four interferograms were complexly space-variant. Due to the fact that the AP and the DP were both highly spatially correlated, the differences of their time-varying characteristics could be utilized to determine the reason for space-variant phases in these four interferograms. The AP was stochastic in time domain, while the DP had a low frequency temporal structure, which was often built as a constant velocity model with a sinusoidal component for periodic fluctuations. B . Large variations occurred within the time duration of one InG and showed no variation trends over time. Therefore, the space-variant phase was caused by a rapidly changing weather condition during the acquisition period of Group B, which was probably caused by a sudden rainfall. The assumption of spatial homogeneity for the AP could no longer be applicable. A special area is marked with a red ellipse in Figure 6c,e,g, where the unwrapped phases of those PSs located inside this elliptical area were obviously different from that of other PSs. The reason for this is that large deformations occurred on this elliptical area during the period of Group B.   B . It can be noted that the compensated phases of most PSs were between −1 and 1 rad. Except for those PSs located inside the elliptical area, some PSs located in partial areas also had compensated phases that obviously deviated away from 0 rad. This illuminates that the AP components of IM were not effectively compensated for.

Problem Analysis
There are two main reasons that the conventional method could not be utilized to effectively compensate for the AP components of the map IM (30) B . Firstly, those PSs located in the elliptical area with large deformations would affect the accurate estimation of the AP parameters. Secondly, the assumption of spatial homogeneity for the AP was not applicable, which meant that the liner phase model could no longer hold. A general solution to discriminate those PSs with large deformations is coherence selection. The coherence γ of a PS between two complex images I M and I S could be calculated with a window function method: where * denotes the complex conjugate. L 1 and L 2 are, respectively, the length and width of the rectangular window. With the PS to be estimated as the window center, two windows with the size of L 1 × L 2 are selected from I M and I S , and all the L 1 × L 2 pixels inside one window are used to calculate its coherence. Figure 9a shows the coherence map calculated based on the interferogram IM where the window size was set as 5 × 5. The PSs with coherences larger than 0.9 made up 96.15%. Figure 9b shows the PSs whose coherence was lower than 0.9. Only some PSs located in the margin of the elliptical area were selected, while those PSs located in the middle of the motional area had large coherence due to the high spatial correlation of the DP. Therefore, the spatial coherence was not effective enough to discriminate those PSs located inside the motional area. Since the liner phase model could not effectively simulate the AP components in the interferogram IM (30) B , other parametric models were utilized. Equation (9) describes a second order model of range, and Equation (10) describes a first order model of the range-angle [21,22]. These two models were utilized to compensate for the AP.
where R and θ are, respectively, the slant range and azimuth angle of a PS at the polar coordinate. β o , β 1 and β 2 are the atmospheric parameters to be estimated. Figure 10a,b, respectively, show the compensated maps with the two models. Apparently, the AP was still not well compensated for. This illuminates that the space-variant AP due to bad weather conditions could not be effectively simulated with these common parametric models.

Improved Method
When processing time-series GB-SAR images, since the weather conditions have strong time-varying characteristics and may inhomogeneously vary in space, the conventional method which builds a parametric model to simulate the AP cannot always be effective and may cause large compensation errors. This paper refers to the method of spatial interpolation based on some stable control points to estimate the AP and then proposes an improved method [23]. According to Equation (1), PSs can be classified into three types: NPS, DPS and APS.
• NPS: A PS with a large NP ϕ noi , i.e., an unstable noise-dominant PS. • DPS: A PS with a small NP ϕ noi but a large DP ϕ defo , i.e., a deformation-dominant PS. • APS: A PS with a small NP ϕ noi and a small DP ϕ defo , i.e., an atmospheric effect-dominant PS.
Taking advantage of the different spatial properties between the AP, the DP and the NP, those APSs were selected, and a number of control points were determined. Then, a spatial interpolation method was utilized to estimate and compensate for the space-variant AP. The improved method mainly includes three steps: NPS rejection, DPS rejection, and APS compensation.

NPS Rejection
According to Equation (1), the unwrapped phase of a PS could be modelled as the sum of the DP, the AP and the NP. Moreover, these components are spatially and temporally independent of each other. Therefore, the standard deviation σ InP of phase sequence of a PS could be described as: where σ defo , σ atm and σ noi , respectively, denote the deviation of the DP, the AP and the NP. Since the AP is space-variant, σ atm also changes with the PS position.

For the 30 interferograms IM
(1) A of Group A, the standard deviation of the phase sequence of every PS was calculated. Figure 11a,c, respectively, show the deviation maps before and after APC with the conventional method. It is worth noting that the conventional method models the AP as a linear phase ramp over the range axis, which is the same explanation in the following text. Before APC, the deviations of most PSs were between 0.1 and 0.3 rad, which was mainly caused by the AP and the NP. Moreover, their deviations increased with the PS range. This illustrates that σ atm is space-variant. After APC, the deviations of most PSs were lower than 0.1 rad. Therefore, phase fluctuation caused by the NP was typically smaller than 0.1 rad. However, the deviations of some PSs located in the upper side of the open pit were large, which could be rejected with Equation (7) when taking a second estimation of the atmospheric parameters.   Figure 11b,d, respectively, show the deviation maps before and after APC. Before APC, for those PSs located inside the elliptical motional area, their deviations were within the range of 0.6-1 rad. For other PSs, their deviations were basically larger than 0.2 rad, which was much larger than those of Group A and was mainly caused by the space-variant AP. After APC, the PSs with deviations larger than 0.2 rad still made up about 32.1%. Since that the conventional method was not effective enough, some AP components, i.e., the APC errors, still remained in the phase sequences of the PSs, and the NPSs could not be rejected just by setting a deviation threshold.
Taking two neighboring PSs, PS (m) and PS (n) (which are rather close in space), their differential phase ∆ϕ (m,n) could be described as: The atmospheric terms ϕ Inp of the differential phase sequence of two neighboring PSs could be modeled as: A Delaunay triangulation network was constructed to connect all the scattered PSs. The Delaunay triangulation has the property that the circumcircle of any triangle in the network contains no other points in its interior [24]. We analyzed the side length according to the PS positions and deleted those long sides whose lengths were larger than 3 m. We defined any two PSs which were connected by one of the remaining short sides as a pair of neighboring PSs. If a PS had no neighboring PS, i.e., isolated PS, it was regarded as a NPS and rejected from the PS set. We calculate the standard deviation SD NBPS of the differential phase sequence of a pair of neighboring PSs. If a PS had more than one neighboring PS, the deviation mean was calculated. The deviation SD NBPS could be utilized to evaluate the phase otherness between a PS and its neighboring PSs. Figure 12a,b show the deviation maps calculated with the 30 interferograms of Group A and B, respectively. Figure 12c,d show their scatter diagrams, where both of their deviation SD NBPS increased with the PS range. The reason for this is that the SNR decreased with the PS range and caused the increase of the noise deviation. Therefore, a deviation threshold which was related to the PS range could be utilized to select the NPS. To simplify this problem, a threshold linear with the PS range was set which varied from 0.1 rad at the nearest range of 400 m to 0.2 rad at the farthest range of 850 m, as shown with the red dotted lines. Figure 12e,f, respectively, show the selection results, where the red pixels denote the NPSs and green pixels denote other PSs. For Group A, the total PS number was 51,235 and 731 NPSs, which were mainly located in the margin of the open pit, were selected as shown with a red ellipse in Figure 12e. For Group B, the total PS number was 69,579 and 2525 NPSs, which were stochastically distributed on the whole image, were selected. Due to the sudden rainfall during the acquisition period of Group B, the surface water might have affected the scattering properties of some PSs and caused a large noise deviation. It is worth noting that for the elliptical area in Figure 12f, only some NPSs were selected due to the fact that the deviation σ (m,n) defo in Equation (13) could not be qualitatively analyzed if a pair of neighboring PSs was located inside the motional area.

DPS Rejection
To reject the DPS, the differences between the spatial properties of the AP and the DP were utilized. Though the AP variation might have been complicated on one image, its variation sequences of different PSs even far apart from each other in a group of images were still highly correlated [25]. By contrast, the DP sequences were only locally correlated for those DPSs.

Cluster Partition
After the NPS rejection, the remaining PSs could be divided into a certain number of clusters based on the K-means clustering algorithm, which is a method commonly utilized to automatically partition a point set into k clusters [26]. It proceeds by selecting k initial cluster centers and then iteratively refining them until there is no further change in the assignment of points into clusters.
Taking the green PSs in Figure 12f as an example, a total number of 67,054 PSs remained after the NPS rejection. According to the coordinates (R sin θ, R cos θ) of every PS, where R and θ are its slant range and azimuth angle, these PSs could constitute a point set. By using the K-means algorithm and setting the mean point number belonging to each cluster as 50, they were divided into 1341 clusters. Figure 13a shows the cluster partition map, where pixels with the same color belong to a same cluster and the black square points denote the cluster centers. These clusters could be classified into three types:

•
Type I: A cluster with PSs located entirely outside the motional area.

•
Type II: A cluster with PSs located entirely inside the motional area.

•
Type III: A cluster with PSs located across the motional and motionless areas.  Figure 13b shows an enlarged map of the elliptical motional area in Figure 11b. Each convex polygon is the minimum enclosing polygon of a cluster. The cluster labeled with blue polygon belongs to Type I, that labelled with the red one belongs to Type II, and that labelled with black one belongs to Type III.

Cluster Selection
A Delaunay triangulation network was constructed to connect all the cluster centers. Figure 14a,b, respectively, show the triangulation networks constructed with the cluster centers of Group A and Group B. We calculated the side length that connected any two cluster centers and deleted those long sides whose lengths were larger than 30 m. We defined any two clusters connected by one triangle's side as a pair of neighboring clusters. If a cluster had no neighboring cluster, it was reconnected with a nearest one. For a cluster, its mean phase sequence was acquired by averaging the phase sequences of all the PSs belonging to it. Apparently, the NP could be well filtered out with a spatial phase average. We then calculated the deviation of the difference sequence of two neighboring clusters and denoted it as SD NCS . Three cases exists regarding the value of SD NCS , as follows: 1.
Case 1: If two neighboring clusters both belong to Type I, the AP and the NP are almost eliminated. Therefore, SD NCS is usually smaller than 0.1 rad.

2.
Case 2: If two neighboring clusters belong to Type I and Type III, the DP in the difference sequence cannot be neglected and causes SD NCS to be large.

3.
Case 3: In other cases, the DP in the difference sequence cannot be quantitatively analyzed, and the value range of SD NCS is unknown. Figure 14a shows the deviation map of Group A, where the color of every triangle side represents the value of SD NCS . For most pairs of neighboring clusters, SD NCS was smaller than 0.1 rad. We set a simplified phase threshold linear with the range axis, which varied from 0.1 rad at the nearest range of 400 m to 0.2 rad at the farthest range of 850 m. No side could be selected. As for the deviation map of Group B shown in Figure 14b, some sides located around the elliptical motional area were selected. Figure 14c shows an enlarged map of the red dotted box in Figure 14b. Twenty-six sides which connected 26 clusters were selected, as shown in Figure 14d. These black square points denote the cluster centers, of which the red and green points denote the selected clusters. Moreover, the phase deviation of the green cluster was larger than that of its connected red cluster. Apparently, these green clusters were all located inside the elliptical motional area. However, at least one cluster that is denoted with a blue circle, which was located inside the motional area and belongs to Type II, was not selected. Such a case corresponds to Case 3, and a second selection was necessary to accurately locate the DPS.
The connectivity of the selected sides was utilized to judge the number of motional areas in the map. Then, every motional area was separately selected with those sides belonging to it. As shown in Figure 14d, any two of these 26 sides could be connected together with other sides. Therefore, there was only one motional area, and its shape could be approximated by connecting those marginal clusters, as shown with a black dotted polygon in Figure 15. Three other clusters denoted with the blue square points are inside it. They are also part of the motional area and should be selected secondly.

DPS Selection
Those PSs belonging to the selected 29 clusters of Group B were exacted, as shown in Figure 16a. The color of each PS denotes the deviation σ InP in Equation (11) of its phase sequence. The black square points denote the 29 cluster centers, and their minimum enclosing polygon is shown with the red lines. For the 12 marginal clusters which are crossed by the polygon, they belong to Type III, and some PSs with no DP were rejected. For a marginal cluster, after calculating the deviation SD MCS of its mean phase sequence, those PSs belonging to it with the deviation σ InP smaller than SD MCS were regarded as APSs and were rejected from it. Taking the same process to every marginal cluster, the selection result of the DPSs is shown in Figure 16b, where the red lines denote the minimum enclosing polygon of these DPSs, and this polygon area could be regarded as the selected motional area.

APS Compensation
After the rejection of the NPSs and DPSs, the residual PSs were APSs, and they were reclassified into a certain number of clusters. These cluster centers were regarded as control points. For one interferogram, the mean phase ϕ (m) CP of the m th control point was acquired by averaging the phases of all the APSs belonging to the m th cluster. With a spatial phase average, the NP was effectively filtered out. Based on these control points, the AP of all the PSs was estimated with an inverse distance weighting interpolation (IDWI). The estimated APφ (n) AP of the n th PS could be calculated with the IDWI as: where λ m is the weighting factor, |d m | denotes the distance between the n th PS and the m th control point, M is the total number of the control points, and u is the power exponent whose common value is 2.
To simplify the IDWI, just three control points were used for every PS to estimate the AP. A Delaunay triangulation network was reconstructed to connect the M control points. Every vertex of a triangle corresponded to a control point. If a PS was located inside one triangle, three vertexes of the triangle were selected. If a PS was located outside of any triangle, three closest vertexes were selected. With the selected three vertexes, the AP of a PS was estimated with Equation (14).  Figure 17b shows the triangulation network that connected these control points. Figure 17c shows the estimated AP map. Figure 17d shows the compensated map with the improved method after removing the estimated AP from the map IM  With the improved method, the interferogram IM

(30)
A was also processed. Figure 18a,b, respectively, show the estimated AP map and the compensated map. It seems that the improved map was much similar than the conventional one shown in Figure 7a.

Short-Term APC
With the conventional and improved methods, the 30 interferograms IM of Group B were both processed. To compare the compensation results of the short-term AP with both methods, some PSs with a low ADI were firstly selected to compare their phase curves. Then, the deviation of every PS's phase sequence in these compensated interferograms was calculated. The deviation maps were utilized to make a second comparison.
Seven PSs A 1 to A 7 located at different ranges with an ADI lower than 0.05 were selected from Group A, as marked with the red square points in Figure 18b. Figure 19a,c,e, respectively, show their phase variation curves before the APC and after the APC with the conventional and improved methods. Their phase variations were remarkably similar before the APC and within ±0.1 rad after APC with both methods. Another seven PSs B 1 to B 7 were selected from Group B, as marked in Figure 17d. Figure 19b,d,f, respectively, show their phase variation curves before the APC and after the APC with both methods. Their phase variations were also remarkably similar before the APC, but they were different after the APC with both methods. As for the conventional result, the phase variations could reach ±0.7 rad. Considering a phase error of 0.7 rad, a deformation error of about 0.93-1.39 mm could be caused for a GB-SAR system which works with the Ku band (12)(13)(14)(15)(16)(17)(18). What is particularly noteworthy is that the residual APC error was still spatially correlated. This spatial correlation would make it easy for users to mistake some areas for motional areas. As for the improved result, except for the seventh PS B 7 whose phase fluctuation was large, the other six PSs B 1 to B 6 had compensated phases within ±0.15 rad.  To better compare the compensated curves with both methods, other equipment should be used. However, since the mine slopes were rather steep and the slope angles could reach about 50 • , equipment which must be installed on the slope surface, such as a Global Positioning System, levelling, and total stations could not be used. Remote monitoring techniques such as laser scanning and photogrammetry, which only operate at good weather conditions and have no continuous monitoring ability, were not installed during the monitoring period. Therefore, quantitative comparisons were made by analyzing the standard deviations of the compensated curves of these selected PSs A 1 to A 7 and B 1 to B 7 . Since the time-duration of one InG was just about 1.5 h, it was impossible for these selected PSs that were separated far from each other to simultaneously move within a short time-interval. The deviation comparisons could be utilized to validate the superiority of the improved method to some extent. Table 1 shows the deviation comparisons with both methods. As for Group A, the phase deviations of these seven selected PSs A 1 to A 7 with both methods were small-about 0.010 and 0.006 rad for the conventional and improved methods, respectively. Therefore, for the interferograms of Group A whose phases could be well modelled as a linear ramp over the range axis, both methods could achieve good APC results. As for Group B, the deviations of the other seven PSs B 1 to B 7 with the conventional method were larger than 0.100 rad and much larger than those with the improved one. Only B 7 seemed to move slightly during the period of Group B. Therefore, for the interferograms of Group B whose phases presented complicated spatial variations, the improved method could achieve much better APC results than the conventional one.  Figure 20a,b, respectively, show the deviation maps calculated with every PS's phase sequence of Group A and Group B after the APC with the improved method. Figure 11c,d, respectively, show the deviation maps of Group A and Group B with the conventional method. As for Group A, the deviation maps with both methods were almost the same, and only a small number of PSs located at the far range had deviations larger than 0.1 rad. As for Group B with the conventional method, the PSs with deviations smaller than 0.1 rad made up just about 7.61%, and those with deviations smaller than 0.2 rad made up about 70.41%. With the improved method, the PSs with deviations smaller than 0.1 rad made up about 59.98%, and those with deviations smaller than 0.2 rad made up about 92.88%. The PS deviations were significantly reduced with the improved method.

Long-Term APC
The 5070 GB-MIMO images acquired during the experimental period of 11 days were processed with both methods. Figure 21a shows the cumulative phase map before the APC, where the elliptical motion area mentioned above is masked for a better display. The phase map is rather noisy and shows complicated spatial variation due to the long-term atmospheric effects. Figure 21c shows the compensated phase map with the conventional method. Even after APC, the compensated phase of most PSs still deviated from 3 rad, whose corresponding deformations could reach 3.98-5.97 mm for a Ku-band GB-SAR system. As can be seen from analyzing the map intuitively, a large portion of the pit slope above the middle road is light yellow, and deformations of several millimeters have occurred on its surface. Moreover, the phase map seems to be uniformly variant along the range axis, which is important for determining whether it was caused by the DP or the AP. Firstly, for the pit slope above the middle road, its unstable state on such a large scale was questioned by the mine management. Secondly, although the residual error of APC was complexly space-variant over time, as shown in Figure 8a, its local variation could be reduced with the temporal cumulation to some extent. Thirdly, it was impossible for the DP to occur on nearly the whole image. Therefore, the uniformly variant phase was mainly caused by the residual error of APC. Figure 21e shows the compensated phase map with the improved method. Its overall color is green, which stands for 0 rad, and only a small number of PSs with other colors randomly distributed on the phase map.
Six PSs P 1 to P 6 located at different ranges with high coherence during the experimental period were selected, as marked with the red square points in Figure 21a. Figure 21b shows their phase variation curves before the APC, whose variation trends over time were remarkably similar. Since these selected six PSs were far apart from each other, the similar trends were mainly caused by the spatially correlated AP. Apparently the AP could not be described as a zero-mean random variable over time. The phase error caused by the time-variant weather conditions could reach several or dozens of radians over a long time period, which fully demonstrates the necessity of APC. Figure 21d shows the compensated curves of P 1 to P 6 with the conventional method. Before 4 April 2018, their phase randomly varied around 0 rad but later gradually cumulated and positively increased with range, which was caused by the residual error of APC. Consequently, the residual compensating error accumulated with time. The reason for this is that the elliptical area started to slide after 4 April 2018, and the weather conditions in those days were bad, which caused the conventional method to not be applicable. By contrast, Figure 21f shows the compensated curves of P 1 to P 6 with the improved method. For the six PSs P 1 to P 6 , their phase curves were close to random variation around 0 rad during the whole period of 11 days. Figure 21. Cumulative phase maps: (a) Before the APC, after the APC (c) with the conventional method, and after the APC (e) with the improved method. Phase variation curves of six selected PSs P 1 to P 6 : (b) Before the APC, after the APC (d) with the conventional method, and after the APC (f) with the improved method. Figure 22a shows the cumulative deformation map of the monitoring scene with the proposed method. The maximum deformation for the elliptical area mentioned above could reach about −430 mm, where the minus sign represented the direction approaching the radar. Then, the deformation map was back projected onto a 3D terrain cloud that was acquired with a laser scanner to better locate the motional area, as shown in Figure 22b. The mine workers started to excavate the pit bottom after 4 April 2018. Then, a continuous sliding was caused along the pit slope and formed an elliptical motion area. Figure 22c

Discussion
The atmospheric condition is time-space variant, and its variation could cause a severe phase error. In general, the AP is spatially homogeneous, and it could be modelled as a linear phase ramp over the range and be well compensated for with the conventional method based on the PS technique, as shown in Figure 7. However, due to occurrences of rainfall, snowfall, etc., the AP might be complexly space-variant and cannot be described with a proper parametric model, as shown in Figure 6.
An improved method to compensate for the time-space variant AP applied to time-series GB-SAR images has been proposed. Firstly, based on the fact that the NP in Equation (1) was spatially uncorrelated, the standard deviation SD NPS of the differential phase sequence of a pair of neighboring PSs was defined. By setting a proper linear threshold with the PS range, the NPSs were rejected and are shown in Figure 12. Then, by taking advantage of the differences between the spatial properties of the AP and the DP, the remaining PSs were partitioned into a large number of clusters according to their spatial positions. With the cluster selection and analysis, those DPSs were selected and are shown in Figure 16. With the APSs, the space-variant AP was compensated for with a spatial interpolation method, as shown in Figure 17.
To validate the effectiveness and feasibility of the proposed method, comparisons were made with the conventional method from two aspects. As for the short-term APC, the interferograms of Group A and Group B were processed with both methods. Seven reference PSs were selected to intuitively compare their compensated curves, as shown in Figure 19; a comparison of their deviations is shown in Table 1. As for the long-term APC, the 5070 GB-MIMO images were processed with both methods. For the phase map compensated for with the conventional method, as shown in Figure 21c, it seems that deformations of several mm occurred on a large portion of the open pit, and these were caused by the residual APC error. Six reference PSs were also selected to make a further comparison. The compensated curves shown in Figure 21d,f prove that the time-variant AP could be better compensated for with the improved method.

Conclusions
This paper emphatically analyzed an improved compensation method of a time-space variant AP applied to time-series GB-SAR images. With a GB-MIMO radar system, 5070 images of an open-pit mine were acquired during an experimental period of 11 days. The analysis results of a time-series of phase maps show that in a few cases, the AP might have been complexly time-space variant and could not be simulated with a proper parametric model where the conventional method based on the PS technique was no longer applicable. The proposed method takes full advantage of the different spatial properties between the AP, the DP and the NP to classify the NPS, the DPS and the APS. With the APSs, a certain number of control points were determined and utilized to compensate for the space-variant AP of all the PSs with a spatial interpolation. Short-term and long-term experimental datasets were both processed with the conventional and improved methods. Moreover, quantitative and qualitative comparisons were made to validate the effectiveness of the proposed method. For a short-term AP, the phase deviations of some reference PSs with a low ADI could be effectively decreased. For a long-term AP, the space-variant AP could be well compensated for. The improved method could effectively improve the measurement accuracy and avoid the misunderstanding of motional areas.