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

Abstract

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.

1. 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.

2. 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.

3. The Proposed PSI Approach

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

• 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.
• 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.
• Generation of two redundant networks of interferograms: full-resolution (pixel footprint: 4 by 14 m) and 10 in range by 2 in azimuth (10 × 2) multi-look (pixel footprint: 40 by 28 m).
• Candidate PS selection using the Amplitude Dispersion index [18].
• The 2 + 1D phase unwrapping of the redundant 10 × 2 multi-look interferograms, see for details [16].
• Identification of stable areas in the surroundings of the area of interest.
• 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.
• Prediction and removal of the estimated atmospheric component from the original single-look interferograms.
• Using the atmospheric-free single-look interferograms, estimation of linear deformation velocity and RTE using the periodogram, see [23].
• Removal of the RTE from the atmospheric-free single-look interferograms.
• The 2 + 1D phase unwrapping of single-look (RTE- and atmospheric-free) interferograms.
• Generation of the deformation time series and estimation of the deformation velocity.
• Geocoding of the results: the deformation velocity and the deformation time series.

4. Data Description

The study area concerns the surroundings of the Glories Square in Barcelona, see Figure 1. This area is located in the Besòs river delta, that is in a sedimentary environment with (sub-)horizontal layers. Below a 6-m thick anthropogenic deposits layer, two stratigraphic levels can be identified: quaternary units and tertiary (Pliocene) units. The quaternary materials are mainly clays with sand and gravel interbeds. The tertiary units are composed of medium to coarse sands from NE to SE and from SW to NW around the Glories Square, and mainly of grey marls from NW to NE and from SE to SW. The processed area is bordered by the red line shown in Figure 1, while the main area of interest has a perimeter indicated by a yellow line. In order to define this area, based on aquifer hydraulic properties and piezometers located in the area of interest, it was assumed that water pumping can have a maximum influence area (i.e., the area where changes in the water table can have effect on the surface) with a radius of 1 km, centred in the middle of the Glories Square. The remaining areas, i.e., outside the influence area, were considered stable. In this specific case there were previous PSI processing results that indicate the stability of this area. However, as it is discussed later, the hypothesis of the stable area can be validated during data processing and analysis. The analysed dataset includes 78 descending IWS SLC Sentinel-1A images, which cover the period from March 2015 to November 2017, see

Table 1. Dates of the 78 descending IWS images used in this work. Relative orbit: 110.

#	Date	#	Date	#	Date	#	Date
1	20150306	21	20151219	41	20160827	61	20170424
2	20150318	22	20151231	42	20160908	62	20170506
3	20150330	23	20160112	43	20160920	63	20170518
4	20150411	24	20160124	44	20161002	64	20170530
5	20150505	25	20160205	45	20161014	65	20170611
6	20150517	26	20160217	46	20161026	66	20170623
7	20150529	27	20160229	47	20161107	67	20170705
8	20150610	28	20160312	48	20161119	68	20170717
9	20150704	29	20160324	49	20161201	69	20170729
10	20150716	30	20160405	50	20161213	70	20170810
11	20150728	31	20160417	51	20161225	71	20170822
12	20150809	32	20160429	52	20170106	72	20170903
13	20150821	33	20160511	53	20170118	73	20170915
14	20150914	34	20160523	54	20170130	74	20170927
15	20150926	35	20160604	55	20170211	75	20171009
16	20151008	36	20160628	56	20170223	76	20171021
17	20151101	37	20160710	57	20170307	77	20171102
18	20151113	38	20160722	58	20170319	78	20171114
19	20151125	39	20160803	59	20170331
20	20151207	40	20160815	60	20170412

. The dataset includes 1813 interferograms, which were generated using all possible image combinations, with a limit of one year for the temporal baselines.

5. 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\left(d_K\right)$, whose values are second order statistics [24]:

$$C\left(d_K\right)=\frac{1}{n_T}\cdot {\displaystyle \sum _{i=1}^{n_T}\left\{\left[M\left(P_i\right)-m_M\right]\cdot \frac{1}{n_J}\cdot {\displaystyle \sum _{j=1}^{n_J}\left[M\left(P_j\right)-m_M\right]}\right\}}$$

(1)

where $d_K=K\cdot \Delta$ is the distance from the origin; $K$ is an integer number; $\Delta$ is the function step; $n_T=l\cdot m$ is the total image pixel number, given by the number of line times the number of columns; $M\left(P_i\right)$ 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: $\left(K-1\right)\cdot \Delta <\Vert P_i-P_j\Vert <K\cdot \Delta$. 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. Main characteristics of the 27 phase images that are affected by the strongest atmospheric phase component. Values before atmospheric correction (orig) and after the correction (res).

  Image #	s_tot_orig	s_corr_orig	s_tot_res	s_corr_res	L_corr_orig	L_corr_res
  34	0.63	0.50	0.60	0.47	204	192
  68	0.71	0.51	0.69	0.49	24	0
  45	0.67	0.50	0.62	0.44	27	18
  62	0.64	0.51	0.64	0.51	255	183
  40	0.67	0.54	0.48	0.32	655	0
  39	0.69	0.54	0.63	0.47	26	18
  4	0.78	0.54	0.67	0.44	0	0
  72	0.67	0.54	0.49	0.33	417	0
  48	0.69	0.55	0.62	0.48	417	18
  16	0.72	0.56	0.54	0.35	281	0
  69	0.69	0.56	0.52	0.36	468	0
  37	0.74	0.59	0.55	0.36	24	0
  3	0.71	0.60	0.64	0.52	204	18
  24	0.71	0.61	0.40	0.22	1785	0
  42	0.80	0.63	0.62	0.41	28	0
  76	0.76	0.63	0.59	0.43	553	37
  35	0.77	0.67	0.67	0.55	765	329
  23	0.76	0.69	0.47	0.35	1029	18
  36	0.83	0.72	0.55	0.38	1029	0
  57	0.99	0.75	0.81	0.53	26	0
  41	0.94	0.80	0.59	0.40	1122	0
  67	0.97	0.82	0.64	0.45	842	0
  14	1.01	0.92	0.53	0.37	1496	0
  64	1.16	0.95	1.07	0.86	51	18
  70	1.17	1.08	0.60	0.46	1658	18
  12	1.33	1.21	1.10	0.97	468	384
  44	1.40	1.30	0.98	0.85	638	185

  . 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 that there are five images where the L_corr values remain quite high (i.e., above 180 m) after removing the atmospheric component (images 12, 34, 35, 44, and 62). In two cases (images 12 and 44), the corresponding σ_corr values are also high (0.85 and 0.97 rad). These two images represent a case where there is an important atmospheric component, which cannot be modelled by a linear atmospheric model. In this case there are two options: (i) discarding the images, especially if the dataset is big enough; (2) if the images cannot be eliminated, the deformation time series have to be interpreted with attention: the time series values in correspondence of the two images may contain spikes.
• The remaining images have σ_corr values ranging between 0.22 and 0.55 rad. These values indicate the dispersion of the residual atmospheric signal, which affect the corresponding deformation time series. In terms of displacements, the standard deviations of such a signal, range from 0.97 mm (best case) to 2.42 mm (worst case): these values seem to be acceptable for the purpose of the application at hand, as discussed later.

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.

6. 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 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.

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.

7. 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.