A Persistent Scatterer Interferometry Procedure Based on Stable Areas to Filter the Atmospheric Component

This paper describes a Persistent Scatterer Interferometry (PSI) procedure to monitor the land deformation in an urban area induced by aquifer dewatering and the consequent drawdown of the water table. The procedure, based on Sentinel-1 data, is illustrated considering the construction works of Glories Square, Barcelona (Spain). The study covers a period from March 2015 to November 2017, which includes a dewatering event in spring 2017. This paper describes the proposed procedure, whose most original part includes the estimation of the atmospheric phase component using stable areas located in the vicinity of the monitoring area. The performances of the procedure are analysed, characterising the original atmospheric phase component and the residual one that remains after modelling the atmospheric contribution. This procedure can work with any type of deformation phenomena, provided that its spatial extension is sufficiently small. The quality of the obtained time series is illustrated discussing different deformation results, including a validation result using piezometric data and a thermal expansion case.


Introduction
Deformation monitoring of urban areas is an important tool for city management and asset maintenance. An important application is the monitoring of the deformation caused by construction works that involve aquifer dewatering, which can affect buildings and infrastructures. In this paper we describe one of this type of application: the monitoring of the deformation associated with the construction works related to the transformation process of the Glories Square, located in the centre of Barcelona (Spain). These works involve the construction of underground tunnels, which requires aquifer dewatering. They are monitored with a set of in situ measurements, e.g., inclinometers, topographic surveys, levelling, which are mainly located in the area of the Glories Square. Such measurements are complemented with Persistent Scatterer Interferometry (PSI) observations, which aim at achieving a global view of the deformation phenomena occurring in the square and, especially, its surroundings, where in situ measurements are not available.
The PSI monitoring is based on C-band data acquired by the Sentinel-1A sensor. The Sentinel-1 mission offers significant improvements, with respect to previous European Space Agency missions, in terms of revisiting time, spatial coverage and quality of the Synthetic Aperture Radar (SAR) imagery. The Interferometric Wide Swath (IWS) acquisition mode of Sentinel-1 images provides a 250-km swath with a repeat cycle of 6 days, considering both 1A and 1B. The Sentinel-1 mission has been especially designed for massive wide-area monitoring. An important advantage of the Sentinel-1 mission is that the data are freely available for both scientific and commercial applications. Several studies based on Sentinel-1 data have been devoted to urban deformation monitoring. Relevant examples include the monitoring of Mexico City [1], Madrid [2], Wuhan [3], Shanghai [4], Ravenna [5], the San Francisco Bay Area [6], Beijing [7], the Lanzhou New District [8], Florence [9], and Istanbul [10]. The monitoring of a slow-moving urban landslide area is described in [11], while the monitoring of sink-holes in urban areas is described in [12]. In [13] a high-speed railway bridge is studied, while in [14] the monitoring of a network of roads and railways is focused on.
This paper is organised as follows. Section 2 examines some important characteristics of PSI monitoring in urban areas. Section 3 describes the PSI approach proposed by the authors to monitor the area of interest. Section 4 discusses the results achieved over the area of interest. Section 5 includes the conclusions of this work.

Urban PSI Monitoring
Several PSI deformation monitoring approaches have been proposed in the last two decades, see [15] for a review. The PSI approach usually used by the authors is described in [16]. However, in some cases, it is possible to tailor the PSI data processing and analysis to the specific characteristics of the application at hand. This was the case in the analysis of the Glories Square area, which is described in this work. In the following, we briefly describe the key characteristics of the proposed processing.
The monitoring of the Glories Square and its neighbouring areas primarily concerns buildings and infrastructures. In order to obtain measurement points (hereafter called Persistent Scatterers, PSs) over such elements, it is important to properly model the Residual Topographic Error (RTE) component of the PSI observations for two reasons: for PSI modelling purposes and to properly geocode the PSI results. A second aspect regards the thermal expansion displacements. Some approaches have been proposed to explicitly model and estimate such displacements, e.g., see [17]. It is worth mentioning that these displacements can be neglected if one could filter out the PSs characterised by high RTE values (e.g., those above 10 m). However, in this work, this approach was discarded because the PSs located on the top of the buildings are of major interest, and we decided to keep the thermal displacements together with the ground deformation. This has the advantage of performing a standard PSI analysis, i.e., avoiding the computational costs of explicitly modelling the thermal expansion component. The analysis of the thermal expansion displacements and their separation from the deformation was carried out during data interpretation.
The third and most important aspect concerns the atmospheric phase component. Most of the PSI approaches use sets of spatial and temporal filters to estimate the atmospheric phase component [18][19][20]. The critical issue is to separate two low-pass spatial contributions: the deformation and the atmospheric components. The assumption usually employed is that the former one is temporally correlated, while the latter component is not correlated in time. In the targeted application this approach has two limitations. The first one is the difficulty to properly separate subtle deformation from the atmospheric contribution, with the risk to underestimate the former one. The second limitation is that no prior information is available about the dewatering plans and hence on the expected ground deformation. For instance, water pumping can be intermittent, implying a series of ups and downs of the ground, which cannot be smooth in time. For these reasons, we decided to adopt an alternative approach, which does not make use of spatio-temporal filters and, more importantly, does not make any assumption associated with such filters. We use a different philosophy based on stable areas, which allows us to estimate the atmospheric component without making any assumption concerning the deformation at hand. This has two advantages: we can correctly estimate sudden deformation, without the disadvantage of filtering out the high-frequency temporal components of the deformation; in addition, we can avoid the subjective decisions often associated with filtering. We use a similar procedure to process interferometric ground-based SAR data [21].
This approach requires the availability of known stable areas in the vicinity of the area of interest. This needs external information, which however can be validated during the data analysis. The stable areas are used to estimate the atmospheric component, which is then predicted and removed from the PSI observations over the area of interest. Several approaches can be used for this purpose, e.g., kriging, least squares collocation (e.g., see [22]), polynomials, etc.

The Proposed PSI Approach
In this section we detail, step by step, the proposed PSI approach.

1.
Acquisition of a set of N interferometric Sentinel-1 SAR images that cover the area of interest.
In this work, a minimum of 25 IWS images were used.

2.
Precise co-registration of the entire burst stack that covers the area of interest. This is based on the information provided by the precise orbits associated with the images.
Identification of stable areas in the surroundings of the area of interest. 7.
Estimation of the atmospheric phase component over the stable areas. In the current implementation of the monitoring, this step is performed assuming a linear phase model. The residuals of such models are used to validate the hypothesis regarding the stable areas. 8.
Prediction and removal of the estimated atmospheric component from the original single-look interferograms. 9.
Using the atmospheric-free single-look interferograms, estimation of linear deformation velocity and RTE using the periodogram, see [23]. 10. Removal of the RTE from the atmospheric-free single-look interferograms. 11. The 2 + 1D phase unwrapping of single-look (RTE-and atmospheric-free) interferograms. 12. Generation of the deformation time series and estimation of the deformation velocity. 13. Geocoding of the results: the deformation velocity and the deformation time series.

PSI Results over the Test Area
First of all, we analysed the characteristics of the output of step 5, described in the previous section. This output includes a set of 78 of 10 x 2 multi-look unwrapped phase images, where the first image values are set to zero because they represent the temporal reference. The objective of the analysis focused on the stable areas was the inspection of the characteristics of the atmospheric phase component. We assumed that, in the stable areas, this component is the only one spatially correlated (i.e., the RTE phase component, the thermal phase component and the phase noise are assumed to be spatially uncorrelated).
The data were analysed using the empirical autocorrelation function, EAF, C(d K ), whose values are second order statistics [24]: where d K = K · ∆ is the distance from the origin; K is an integer number; ∆ is the function step; n T = l · m is the total image pixel number, given by the number of line times the number of columns; M(P i ) is the image value in the pixel P i ; m M is the mean image value; and n J is the total number of pixels P j that, for a given P i , satisfy the condition: Two examples of EAF are shown in Figure 2. From this function, the following information can be derived: (1) σ tot , the total standard deviation of the phase image; (2) σ corr , the standard deviation of the spatially correlated part of the phase image; (3) σ noise , the standard deviation of the spatially uncorrelated part of the phase image; (4) L corr , the correlation length, i.e., the distance from the origin where the EAF has a correlation which is half of that in the origin.    The main results of the analysis are summarised below: • A number of 24 images out of 77 (i.e., 31.1% of the total) have σ corr <0.4 rad. This indicates a rather weak atmospheric component. In terms of displacement this corresponds to a standard deviation below 1.76 mm. A number of 14 images out of 24 were acquired during winter or late autumn: this confirms that this is the period of the year when less atmospheric turbulence occurs. With the exception of one image, with L corr = 680 m, this group is characterised by zero or negligible L corr . • A number of 26 images out of 77 (i.e., the 33.8% of the total) have σ corr between 0.4 and 0.5 rad.
With the exception of one image, with L corr = 833 m, this group is characterised by moderated L corr values, in the range between 20 and 50 m. • Finally, 27 out of 77 images (i.e., the 35% of the total) have σ corr above 0.5 rad, see their main characteristics in Table 2. A number of 18 images out of 27 were acquired in summer or late spring: this is the period when there is maximum atmospheric turbulence. With the exception of 7 images, which have negligible L corr values, all the remaining images have correlation lengths above 50 m, with a maximum value of 1785 m. As mentioned in the previous section, the atmospheric component was estimated over the stable areas using a linear phase model, which represents an easy-to-implement and robust modelling approach over small areas. The original 78 phase images, the atmospheric models (the linear planes) and the residual phase images after removing the models are shown in Figure 3. The hypothesis of the stable area can be validated by analysing the latter ones. In fact, a deformation signal would typically appear as a persistent pattern in consecutive images. This is not the case in the analysed dataset. The main results of the EAF analysis of the residual phase images displayed in Table 2 are:

•
Compared to the original data, the average reduction of σ corr is 30.5%. The most relevant result is a drop of 90% of the average L corr values. This is an important indicator to judge the goodness of the proposed method. An example of correlation drop is shown in Figure 2: in this case the L corr values drop from 1496 m to zero. • It is worth noting from Table 2   The last part of the analysis concerned the 1-km radius area of interest (influence area, Figure 1). The objective was to assess the effects of the residual atmospheric component in the deformation time series (computed in the step 12 of the proposed PSI procedure). The analysis concerned the stable PSs of the main area of interest: they were selected by only considering those points that have an absolute The last part of the analysis concerned the 1-km radius area of interest (influence area, Figure 1). The objective was to assess the effects of the residual atmospheric component in the deformation time series (computed in the step 12 of the proposed PSI procedure). The analysis concerned the stable PSs of the main area of interest: they were selected by only considering those points that have an absolute deformation velocity below 1 mm/yr. The time series were analysed using the EAF: all the 78 images were included in the analysis (no images were discarded). These are the main results:

•
The average σ tot of all the time series is 1.90 mm. Only 127 out of 3862 time series have σ tot above 3 mm: this represents 3.3% of the PSs. It is worth observing that the σ tot includes, among others, the residual (non-modelled) atmospheric effects and the noise of the observations, which depends on the PS quality. • As expected, the L corr is close to zero for the great majority of time series. This is key to detect subtle (temporally correlated) deformation using the time series.

Deformation Results and Discussion
We discuss in the following some examples of time series, to illustrate the goodness of the proposed procedure. Figure 4 shows the line-of-sight (LOS) deformation time series of a point located close to the Glories Square. In this figure, during the period from March to June 2017, there is a terrain subsidence (up to about −10 mm), followed by an uplift to roughly recover to the original height. It is worth underlining that this result was achieved using exclusively SAR data, without any additional information. In the interpretation of the data, it was identified that this behaviour was due to aquifer dewatering, and the subsequent recovery of the water level due to the stop of the pumping activity. This is evident from the piezometer data plotted in Figure 4, which were acquired in a location close to the point, and which match very well with the estimated deformation time series. This represents an example of validation of the PSI results.

Deformation Results and Discussion
We discuss in the following some examples of time series, to illustrate the goodness of the proposed procedure. Figure 4 shows the line-of-sight (LOS) deformation time series of a point located close to the Glories Square. In this figure, during the period from March to June 2017, there is a terrain subsidence (up to about −10 mm), followed by an uplift to roughly recover to the original height. It is worth underlining that this result was achieved using exclusively SAR data, without any additional information. In the interpretation of the data, it was identified that this behaviour was due to aquifer dewatering, and the subsequent recovery of the water level due to the stop of the pumping activity. This is evident from the piezometer data plotted in Figure 4, which were acquired in a location close to the point, and which match very well with the estimated deformation time series. This represents an example of validation of the PSI results. To illustrate the performance of the proposed procedure, in Figure 4 the above time series is compared with a solution based on spatio-temporal filters, see the procedure described in [16]. In this case we used a 96-day moving average. In this example, it is evident that the first deformation time series matches the sudden drop of the piezometric data, while the second solution is basically missing such a drop, providing a biased temporal low-pass solution. One may observe that in the time series there are several peaks that are not linked to a particular period of dewatering. They have an amplitude of ± 2-3 mm with respect to the global trend of the time series. They are not due to thermal expansion: they are examples of the residual atmospheric effects, discussed earlier in this paper. Figure 5 shows the LOS accumulated deformation maps corresponding to the maximum of the displacement (12 April 2017) and to the recovery of the displacements (3 September 2017). From the first map it is possible to assess the actual dimension of the area affected by ground deformation In the above image, the three rectangles (grey-zone 1, orange-zone 2 and pink-zone 3) show the location of the three zones shown in Figure 6. The green circle A shows the location of the point considered in Figure 4. The white circle shows the location of the piezometer.  In the above image, the three rectangles (grey-zone 1, orange-zone 2 and pink-zone 3) show the location of the three zones shown in Figure 6. The green circle A shows the location of the point considered in Figure 4. The white circle shows the location of the piezometer.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 12 Figure 5. Examples of LOS accumulated deformation maps corresponding to the maximum of the displacement (12 April 2017, above) and to the recovery of the displacements (3 September 2017, below). In the above image, the three rectangles (grey-zone 1, orange-zone 2 and pink-zone 3) show the location of the three zones shown in Figure 6. The green circle A shows the location of the point considered in Figure 4. The white circle shows the location of the piezometer.  To illustrate the performance of the proposed procedure, in Figure 4 the above time series is compared with a solution based on spatio-temporal filters, see the procedure described in [16]. In this case we used a 96-day moving average. In this example, it is evident that the first deformation time series matches the sudden drop of the piezometric data, while the second solution is basically missing such a drop, providing a biased temporal low-pass solution. One may observe that in the time series there are several peaks that are not linked to a particular period of dewatering. They have an amplitude of ± 2-3 mm with respect to the global trend of the time series. They are not due to thermal expansion: they are examples of the residual atmospheric effects, discussed earlier in this paper. Figure 5 shows the LOS accumulated deformation maps corresponding to the maximum of the displacement (12 April 2017) and to the recovery of the displacements (3 September 2017). From the first map it is possible to assess the actual dimension of the area affected by ground deformation induced by dewatering. Figure 6 shows some additional examples of LOS deformation time series of three zones located in the deformation area shown in Figure 5. The time series display the median values of the points contained in each zone. One may appreciate that there is a good agreement between the measured time series and the piezometric data. Finally, Figure 7 shows the displacement time series of PSs located outside the maximum potential influence area. The location of such points is shown in Figure 1. The time series correspond to basically stable areas, with most of the time series values between ±2 mm. The time series include some spikes, with absolute values up to 6 mm, which are mainly due to residual atmospheric effects. Another example of deformation time series is shown in Figure 8. In the same figure, the temperatures of the scenes are plotted in correspondence to the days of acquisition of the SAR images. One may appreciate a strong correlation between deformation and temperature. This is clearly a displacement behaviour induced by thermal expansion. This type of result is possible only if an appropriate estimation of the atmospheric phase contribution is carried out: this confirms the goodness of the procedure proposed in this work.  Another example of deformation time series is shown in Figure 8. In the same figure, the temperatures of the scenes are plotted in correspondence to the days of acquisition of the SAR images. One may appreciate a strong correlation between deformation and temperature. This is clearly a displacement behaviour induced by thermal expansion. This type of result is possible only if an appropriate estimation of the atmospheric phase contribution is carried out: this confirms the goodness of the procedure proposed in this work.
Another example of deformation time series is shown in Figure 8. In the same figure, the temperatures of the scenes are plotted in correspondence to the days of acquisition of the SAR images. One may appreciate a strong correlation between deformation and temperature. This is clearly a displacement behaviour induced by thermal expansion. This type of result is possible only if an appropriate estimation of the atmospheric phase contribution is carried out: this confirms the goodness of the procedure proposed in this work. This paper has focused on a PSI procedure to monitor the land deformation associated with construction works that involve water pumping and hence lowering of the water table. The case study focuses on the construction works related to the transformation process of the Glories Square, located in the centre of Barcelona (Spain). The PSI monitoring was based on C-band data acquired by Sentinel-1 sensors. The used PSI procedure was tailored to the specific characteristics of the application at hand. The performances of the proposed procedure were analysed. This involved the characterization of the original atmospheric phase component of 78 Sentinel-1 SAR images using the EAF and, in particular, the standard deviation of the spatially correlated part of the phase images, σcorr, and their correlation length, Lcorr. In the analysed case study, the atmospheric component was estimated and removed from the original PSI observations using a linear phase model. The same EAF analysis was carried out on the residual atmospheric phase component. The most important aspect is a drop of 90% of the average Lcorr values: this indicates that the proposed atmospheric estimation procedure works properly. However, in the analysed case, at least two images (2.5% of the images) show relatively high σcorr and Lcorr values: they represent two cases where a strong atmospheric component cannot be properly modelled by the linear atmospheric model used. However, this 2.5% does not compromise the quality of the estimated deformation time series. An EAF analysis of such This paper has focused on a PSI procedure to monitor the land deformation associated with construction works that involve water pumping and hence lowering of the water table. The case study focuses on the construction works related to the transformation process of the Glories Square, located in the centre of Barcelona (Spain). The PSI monitoring was based on C-band data acquired by Sentinel-1 sensors. The used PSI procedure was tailored to the specific characteristics of the application at hand. The performances of the proposed procedure were analysed. This involved the characterization of the original atmospheric phase component of 78 Sentinel-1 SAR images using the EAF and, in particular, the standard deviation of the spatially correlated part of the phase images, σ corr , and their correlation length, L corr . In the analysed case study, the atmospheric component was estimated and removed from the original PSI observations using a linear phase model. The same EAF analysis was carried out on the residual atmospheric phase component. The most important aspect is a drop of 90% of the average L corr values: this indicates that the proposed atmospheric estimation procedure works properly. However, in the analysed case, at least two images (2.5% of the images) show relatively high σ corr and L corr values: they represent two cases where a strong atmospheric component cannot be properly modelled by the linear atmospheric model used. However, this 2.5% does not compromise the quality of the estimated deformation time series. An EAF analysis of such time series was carried out: the average σ tot of the time series of stable points is 1.90 mm. Only 3.3% of the time series have σ tot above 3 mm. The quality of the time series was further illustrated considering different deformation result examples. The most important one concerns a validation result, where the PSI estimated deformation matches well with external piezometric data. This result is confirmed by other time series coming from different locations in the main deformation area. A final example regards a thermal expansion case. The strong correlation with the temperature is only possible with an appropriate estimation of the atmospheric phase contribution: this confirms the goodness of the approach proposed in this work.

Conclusions
The most original aspect of this procedure is the estimation of the atmospheric phase component using stable areas located in the vicinity of the monitoring area. This approach overcomes some important limitations of PSI techniques that use sets of spatial and temporal filters. In particular, it avoids the assumption that the low-pass spatial deformation is temporally correlated. This implies that, with the proposed procedure, sudden deformation can properly be estimated. In addition to the PSI approach described in this paper, the authors routinely apply this procedure to process interferometric ground-based SAR data.
The proposed PSI procedure was developed and tested for a specific application, i.e., the monitoring of urban land deformation due to water extraction. However, the same procedure can work with any type of deformation phenomena, provided that its spatial extension is sufficiently small. This is an important characteristic of the proposed procedure: it works only over relatively small areas, where the atmospheric component, estimated over the stable areas, can be interpolated in the area of interest. The larger the area of interest, the bigger will be the error in the estimation of the atmospheric component over this area (i.e., the larger will be the residual atmospheric component). It is worth noting that the interpolation of the atmospheric component is only possible if stable areas surround the area of interest. In some applications this cannot be possible: in this case the atmospheric component needs to be extrapolated. This implies larger errors in the estimation of the atmospheric component.