The 2015–2016 Ground Displacements of the Shanghai Coastal Area Inferred from a Combined COSMO-SkyMed/Sentinel-1 DInSAR Analysis

Abstract

In this work, ground deformation of the Shanghai coastal area is inferred by using the multiple-satellite Differential Synthetic Aperture Radar interferometry (DInSAR) approach, also known as the minimum acceleration (MinA) combination algorithm. The MinA technique allows discrimination and time-evolution monitoring of the inherent two-dimensional components (i.e., with respect to east-west and up-down directions) of the ongoing deformation processes. It represents an effective post-processing tool that allows an easy combination of preliminarily-retrieved multiple-satellite Line-Of-Sight-projected displacement time-series, obtained by using one (or more) of the currently available multi-pass DInSAR toolboxes. Specifically, in our work, the well-known small baseline subset (SBAS) algorithm has been exploited to recover LOS deformation time-series from two sets of Synthetic Aperture Radar (SAR) data relevant to the coast of Shanghai, collected from 2014 to 2017 by the COSMO-SkyMed (CSK) and the Sentinel-1A (S1-A) sensors. The achieved results evidence that the Shanghai ocean-reclaimed areas were still subject to residual deformations in 2016, with maximum subsidence rates of about 30 mm/year. Moreover, the investigation has revealed that the detected deformations are predominantly vertical, whereas the east-west deformations are less significant.

1. Introduction

Differential interferometric synthetic aperture radar (DInSAR) [1,2,3] plays a significant role for the monitoring of Earth’s surface ground deformation, being able to retrieve accurate and timely information on the sensor’s line-of-sight-projected components of the terrain displacements. Over recent years, the availability of large archives of Synthetic Aperture Radar (SAR) images collected from complementary viewing angles and the exploitation of multi-pass DInSAR techniques [4,5,6,7] have made possible the discrimination of the three-dimensional (3D) components of the deformation which has occurred [8,9,10,11,12,13,14,15,16,17,18,19,20]. Theoretically, 3D deformation components can be recovered based on the combination of data produced by (at least) three complementary observation angles. Nevertheless, commonly-used radar sensors are mounted onboard platforms that fly near-polar orbits (i.e., satellite azimuth direction is approximately parallel to north-south direction), and this, unfortunately, makes the detection of the north-south surface displacement components unreliable. Furthermore, SAR data are not typically acquired in the same dates. Thus, an assumption as to the temporal evolution model for the observed deformation has to be made in order to effectively combine Line-of-Sight (LOS) deformation measurements collected from different orbital positions and acquisition times [17,18]. Some technical solutions for combining multiple-satellite LOS measurements, as well as for merging DInSAR data products with other external information (such as that derived from GPS stations), have been found [15,16,17,18,19,20]. Among these solutions, the approach known as Minimum Acceleration (MinA) combination technique, originally presented in [19], implements a straightforward combination of independently-produced multiple-satellite LOS displacement time-series, so as to retrieve the time-series of the 3D components of the terrain displacement.

In this work, a simplified version of the MinA technique is applied to two sets of SAR data collected respectively by the COSMO-SkyMed (X-band) and the Sentinel-1A (C-band) sensors, over the coastal area of the Shanghai metropolis. The megacity is characterized by high urban density and ground subsidence phenomena due to groundwater extraction and consolidation of high-tall buildings, as revealed by several analyses [21,22,23,24,25,26,27,28]. Most of these investigations, however, have in the last decade focused on the highly-urbanized downtown area, using both persistent scatters [4,7] and small baseline (SB)-oriented DInSAR methods [6]. In particular, the mentioned analyses have been conducted using, at most, ERS-1/2 [21], ENVISAT [22,23], ALOS-1 [22,24], TerraSAR-X [25,26,28], and COSMO-SkyMed [23,27,28] synthetic aperture radar (SAR) images. In particular, a multiple-sensor least squares (LS) combination approach, based on the implementation of the Temporarily Coherent Point SAR interferometry processing, has been proposed in [28]. Conversely, the study of the displacement signals affecting the alongshore ocean-reclaimed platforms through DInSAR-based approaches has been revealed to be more challenging. This is because the ocean-reclaimed areas are subject to broad temporal modifications associated with the consolidation phase of soils used for land-reclamation procedures [29], thus leading to severe temporal decorrelation noise effects [30,31] in the generated interferograms. Nonetheless, a few DInSAR investigations focused on the coastal district of the Shanghai megacity have been recently performed. For instance, an analysis of spatial and temporal behavior of deformation at Shanghai ocean-reclaimed platforms, with a focus on the Pudong International Airport area, has been proposed [32] by applying a modified persistent scatterers (PS) method to a high-resolution TerraSAR-X satellite image stack, collected for the period from 2011 to 2012. In order to partially circumvent the limitations due to the scarcity of high-stable PSs in ocean-reclaimed lands, small baseline approaches, which are able to investigate the deformation of distributed scatterers (DS) on the ground, have also been applied [29,33]. These latter investigations relied on the joint use of DInSAR-derived deformation measurements, derived by processing long-lasting sequences of ENVISAT and COSMO-SkyMed differential SAR interferograms, and empirical geotechnical models [29,33,34,35,36]. As a result, some insights on the future evolution of the deformations in the coastal area of Shanghai have been determined [29,33]. It is worth emphasizing that the analyses provided in [29,33] rely on the principal assumption that ocean-reclaimed areas of Shanghai are mostly affected by vertical deformations (i.e., along the up-down direction) that are ascribable to soil compaction mechanisms [34]. On the one hand, the lack of a wide network of Global Positioning Systems (GPSs) over the entire zone of reclamation has prevented, until now, the direct validation/refutation of some deformation models [29,33]. In particular, the hypotheses made in [29,33] are about the direction of surface deformations, which were assumed to be mostly vertical as testified by the results of some geological surveys [37]. On the other hand, a unique opportunity for quantifying the extent, shape, and direction of ground deformation in the ocean-reclaimed areas of Shanghai is offered nowadays by the exploitation of new-generation SAR data collected through the Sentinel-1 mission [38]. Sentinel-1 provides SAR data for land and ocean monitoring applications, and it is composed of a constellation of two satellites, Sentinel-1A and Sentinel-1B, which were launched on 3 April 2014 and 25 April 2016, respectively. The family of two satellites has been developed within the European environmental monitoring program Copernicus. It is worth noting that the repeat cycle of the Sentinel1A-B twin-sensors constellation is reduced to six days, thus permitting suitable and continuously-updated DInSAR analyses.

In this work, the LOS deformation time-series of the Shanghai coastal area, recovered by applying the small baseline subset (SBAS) algorithm to the two sets of COSMO-SkyMed (descending passages) and Sentinel-1A (ascending passages) data, have been combined to retrieve the up-down and the east-west components of deformation from 26 February 2015 to 4 March 2016. To this aim, the MinA technique [19] is applied. The main objectives of our investigation are: (i) to confirm, with experimental data, that the deformation of ocean-reclaimed platforms is mostly vertical; (ii) to corroborate the validity of some of the forecasting deformation models in [29,33], which are useful for the evaluation of the risk due to possible inundations of the highly-urbanized coastal area; (iii) to investigate the displacement signals that affect the public facilities of the Pudong International Airport. The presented results demonstrate the ocean-reclaimed land is actually subject to prevalent up-down movements (with a maximum subsidence rate of 30 mm/year), with the exception of some isolated areas, just in the proximity of the coastline, characterized by maximum lateral east-west deformation rates of 10 mm/year during the 2015–2016 observation period.

2. Study Area

Shanghai is located at the midpoint of the north-south coastline of China on the alluvial plain of the Yangtze River Delta. The megacity is bound in the north by the Yangtze River estuary, to the east by the East China Sea, and to the south by Hangzhou Bay, see Figure 1a. The land area covers about 6000 km² and is mostly flat, with altitudes ranging between 2 and 6 m above sea level. Under the combined influences of the Yangtze River runoff and tidal currents, large amounts of sediments have been deposited in this region, forming a large intertidal area. In order to alleviate the shortage of land resources caused by urban development, a large area has been reclaimed since 1950. Up to now, the new land formed by the reclamation project has reached up to 175 km², mainly including Lingang New City and Pudong International Airport, see Figure 1a,b. Furthermore, in the near future, more new land will be reclaimed in the eastern coastal area to make room for new buildings and public infrastructures.

Of particular interest for our study is the area of the Pudong International Airport, which is one of the most important airports in China. It is located in the eastern sector of the metropolitan area of Shanghai. More than 60 million passengers traveled through the airport in 2016 and it is expected the airport will handle 80 million passengers in 2025. The first phase of the airport construction began in 1997. It was opened to the public in 1999, replacing Shanghai Hongqiao International Airport as Shanghai's international airport. The first runway, located in the western half of the airport, was part of the first construction phase and has been in use since 1999. In 2005, the second runway was constructed and opened. The second round of expansion, including a second terminal and the third runway, also began in 2005 and started operation in 2008 for the Bejing 2008 Summer Olympics. A new round of expansion—including land reclamation, the fourth runway, and the fifth runway—was approved in 2011, and the fourth and fifth runways were completed in 2015 and 2017, respectively. The location of the five runways, which are represented with different colors representing their different construction times, are shown in Figure 1b. In terms of traffic, the new expansion will make the airport become one of the world's top three, after its completion in 2019. The airport is built on the tidal flats of the south bank of the Yangtze River estuary. The western region of the airport is composed of natural sediments, distributed sandy silt, and silty clay [39,40]. The eastern part of the airport is an artificial foundation, which is filled with loose sandy soil (see Figure 1c). Affected by land reclamation and the filling soil properties, obvious ground deformation phenomenon will inevitably occur. Previous studies [38,39] have shown that there are paleo-rivers within the airport region, and the foundation of the area around these paleo-rivers is more unstable than the other regions, which could intensify ground deformations.

3. Data and Methods

3.1. SAR Data

Our investigation relies on the use of two independent SAR datasets collected over the coast of Shanghai by the COSMO-SkyMed (CSK) and Sentinel-1A (S1-A) sensors constellations, respectively. The CSK dataset consists of 61 SAR scenes (descending passes, HH polarization, with a side-looking angle of about 29° and a satellite heading angle of about 198°) collected from 7 December 2013 to 18 March 2016. The second dataset, acquired through the Sentinel-1A mission from 26 February 2015 to 4 April 2017, consists of 33 SAR scenes (ascending passes, VV polarization, with a side-looking angle of about 39° and a satellite heading angle of about 348°). Sentinel-1A data were available in the single-look-complex (SLC) format and acquired through the interferometric wide swath (IW) mode by employing the terrain observation by progressive scans (TOPS) acquisition mode [41], which provides large swath widths of 250 km at ground resolutions of 5 m × 20 m. The full lists of CSK and S1-A SAR data used are included in Table S1 and Table S2 of the Supplementary Material.

3.2. SAR Data Processing

The two SAR datasets were independently processed by employing the multi-pass small baseline subset algorithm (SBAS) [6]. More precisely, CSK data were elaborated using the in-house SBAS processing chain developed at the Italian National Council Research laboratory [42], and the Sentinel-1A data were processed by the SBAS toolbox implemented within the commercial ENVI’s SARScape modules from EXELIS VIS Information Solutions [43]. It is worth noting that the potential use of different DInSAR tools for the retrieval of LOS deformation measurements represents a key issue of the applied MinA technique. SBAS is a well-established technique that allows the detection of the temporal evolution of Earth’s surface deformation by generating mean LOS velocity maps as well as LOS displacement time-series. In particular, SBAS is based on the use of multiple-master multilook interferograms [44] generated after a proper selection of Small Baseline (SB) SAR data pairs. LOS displacement time-series are computed by solving a least-squares (LS) minimization problem, based on the application of the singular value decomposition (SVD) method [45], to the sequence of unwrapped multilook interferograms. The SBAS approach also includes a filtering operation for the atmospheric phase components [4,6]. Interested readers can find additional details on SBAS and its inherent implementation within InSAR processing tools in [6,46]. Starting with available CSK and S1-A SAR data, two groups of small baseline (SB) interferometric SAR data pairs were first selected, and the relevant sequences of SB interferograms were subsequently generated. The distribution of available SAR acquisitions in the temporal/perpendicular baseline plane as well as the selected SB interferometric SAR data pairs are depicted in Figure 2.

Therefore, for each detected coherent pixel, the sensor line-of-sight (LOS) displacement time-series were recovered in correspondence to detected coherent pixels (in radar coordinates) for the CSK and S1-A datasets, separately. Obtained LOS displacement time-series were then combined as described in Section 3.3 to recover the east-west and up-down time-series of deformation relevant to coherent distributed scatterers (DS) on the ground.

3.3. Multiple-Satellite Combination Technique for the Retrieval of 2D Displacement Time-Series

In this subsection, we summarize the basic rationale of the MinA multiple-satellite DInSAR combination technique [19]. As described earlier, after the application of the SBAS technique, the LOS-projected deformation time-series are recovered. LOS measurement products are originally provided in the native radar coordinate grids of both datasets. Accordingly, in order to properly combine the LOS displacement products, so as to discriminate the 2D components of deformation (with respect to the up-down and east-west directions, respectively), the measurements have to be preliminarily referred to a common geocoded grid by applying geocoding operations [47]. Since the sensitivity of LOS measurement to north-south projection is very weak, hereinafter, we present a modified version of the MinA algorithm, which focuses exclusively on the discrimination of the 2D (i.e., up-down and east-west) components of deformation, by completely neglecting the presence of a deformation component along the north-south direction.

To introduce the right mathematical framework, let us refer to the simple case that only two independent sets of multi-orbits SAR data (i.e., in our case the CSK and the S1-A one) are available (this is the case in the Shanghai case study). Let us assume that SAR data are collected at the ordered times ${\mathit{t}}^{\left(\mathbf{1}\right)}\equiv {\left[t_0^{\left(1\right)},t_1^{\left(1\right)},\dots ,t_{Q_1-1}^{\left(1\right)}\right]}^T$ and ${\mathit{t}}^{\left(\mathbf{2}\right)}\equiv {\left[t_0^{\left(2\right)},t_1^{\left(2\right)},\dots ,t_{Q_2-1}^{\left(2\right)}\right]}^T$, respectively. Let $T\equiv {\bigcup }_{j=1,2}{\mathit{t}}^{\left(\mathit{j}\right)}={\left[T_0,T_1,\dots ,T_{Q-1}\right]}^T$ be the whole set of $\mathrm{Q}=Q_1+Q_2$ ordered time acquisitions, and assume that a suitable time overlap between the two SAR datasets is present. Notice that time-overlap between the two datasets is a fundamental requirement of any multiple-orbit combination algorithm. To accomplish this issue, a preliminary operation step that identifies the overlapping time-period and extracts the LOS time-series over the overlapped time is required (see the experimental Section 4). To describe how the algorithm works, let us suppose for the sake of simplicity, that only two different SAR sensors are available and ${\mathit{t}}^{\left(\mathbf{1}\right)}$ and ${\mathit{t}}^{\left(\mathbf{2}\right)}$ are properly time-overlapped. Moreover, for the sake of simplicity, let us also assume that the first dataset is acquired through ascending passes and the second through descending passes. Henceforth, let ${\mathit{d}}^{\left(\mathbf{1}\right)}\left(P\right)\equiv {\left[d_0^{\left(1\right)}\left(P\right),d_1^{\left(1\right)}\left(P\right),\dots ,d_{Q_1-1}^{\left(1\right)}\left(P\right)\right]}^T$ and ${\mathit{d}}^{\left(\mathbf{2}\right)}\left(P\right)\equiv {\left[d_0^{\left(2\right)}\left(P\right),d_1^{\left(2\right)}\left(P\right),\dots ,d_{Q_2-1}^{\left(2\right)}\left(P\right)\right]}^T$ be the geocoded LOS time-series computed in correspondence to the generic homologous pixel P.

Let us start by relating the LOS measurements collected over ascending and descending orbits, namely $d_{LOS,asc}$ and $d_{LOS,desc}$, to the inherent 2D displacement components with respect to up-down and east-west directions. Taking into account the acquisition geometry [19,28] of the ascending and descending passes depicted in Figure 3, it is easy to prove the validity of the following equations:

$$\begin{array}{c}\begin{array}{c}{d}_{LOS,asc}=\mathit{d}·{\widehat{\mathit{i}}}_{\mathit{a}\mathit{s}\mathit{c}}=d_{East}·sin{\theta }_{asc}+d_{Up}·\cos {\theta }_{asc}\end{array}\\ d_{LOS,desc}=\mathit{d}·{\widehat{\mathit{i}}}_{\mathit{d}\mathit{e}\mathit{s}\mathit{c}}=-d_{East}·sin{\theta }_{desc}+d_{Up}·\cos {\theta }_{desc}\end{array}$$

(1)

with ${\theta }_{asc}$ and ${\theta }_{desc}$ being the incidence angles of the electromagnetic waves on the ground for the ascending and descending data tracks, respectively. Note that ${\widehat{\mathit{i}}}_{\mathit{a}\mathit{s}\mathit{c}}$ and ${\widehat{\mathit{i}}}_{\mathit{d}\mathit{e}\mathit{s}\mathit{c}}$ are the versors for the ascending and descending line-of-sight directions. In order to retrieve the time-series of the east-west and up-down components, the available LOS time-series through ascending and descending passes are hence jointly exploited. In particular, the relative LOS displacement at the generic time $t_k^{\left(j\right)},j=1,2$ calculated with respect to the first acquisition time of the relevant SAR dataset $t_0^{\left(j\right)},j=1,2$ , that is $d_k^{\left(j\right)}-d_0^{\left(j\right)},j=1,2$ , is first extracted from the available ascending and descending LOS time-series. Subsequently, it is related to the unknown east-west and up-down velocity components, namely ${\mathit{V}}_{\mathit{E}}={\left[V_{E_1},V_{E_2},\dots ,V_{E_{Q-1}}\right]}^T$ and ${\mathit{V}}_{\mathit{U}}={\left[V_{U_1},V_{U_2},\dots ,V_{U_{Q-1}}\right]}^T$, between the available $\mathrm{Q}$ adjacent acquisition times (related to the whole set $\mathrm{T}$), following the lines described in [19]. As a result, a system of $\mathrm{Q}-2$ linear equations with respect to $\mathrm{M}=2\left(\mathrm{Q}-1\right)$ unknowns is obtained, which can be expressed using matrix formalism as:

$$B·\left[\begin{array}{c}{\mathit{V}}_{\mathit{E}}\\ {\mathit{V}}_{\mathit{U}}\end{array}\right]=\left[\begin{array}{c}{\mathit{d}}^{\left(\mathbf{1}\right)}\\ {\mathit{d}}^{\left(\mathbf{2}\right)}\end{array}\right]$$

(2)

where B is the incidence-like matrix of the adopted linear transformation, which is defined taking into account Equation (1) as follows:

$$B=\left[\begin{array}{cc}{\mathit{B}}^{\left(\mathbf{1}\right)}\sin {\theta }_{asc}& {\mathit{B}}^{\left(\mathbf{1}\right)}\cos {\theta }_{asc}\\ -{\mathit{B}}^{\left(\mathbf{2}\right)}\sin {\theta }_{desc}& {\mathit{B}}^{\left(\mathbf{2}\right)}\cos {\theta }_{desc}\end{array}\right]$$

(3)

with $B^{\left(\mathbf{j}\right)}\mathit{j}=\mathbf{1},\mathbf{2}$ being the matrixes related to the linear transformation from LOS displacements to velocity deformation rates (see [6,7,8,9,10,11,12,13,14,15,16,17,18,19] for additional details).

The system of Equation (2) is undetermined; thus it does not admit a unique solution. The strategy adopted by the MinA technique is to introduce additional equations imposing that the (combined) East-West and Up-Down deformation components are with minimum acceleration. That is to say that the difference of velocity between consecutive time intervals of the whole set of ordered times T is minimal.

This can be mathematically expressed as follows:

$$\{\begin{array}{c}\alpha \left(V_{E_{i+1}}-V_{E_i}\right)=0i=1,\dots ,Q-2\hfill \\ \alpha \left(V_{U_{i+1}}-V_{U_i}\right)=0i=1,\dots ,Q-2\hfill \end{array}$$

(4)

where $\alpha$ is a regularization factor used to balance the relative weight of the minimum acceleration constraints with respect to the minimization of the residual norm of the system in Equation (3).

Accordingly, the regularized system of equations becomes:

$$\left[\begin{array}{c}\mathit{B}\\ \mathit{C}\end{array}\right]·\left[\begin{array}{c}{\mathit{V}}_{\mathit{E}}\\ {\mathit{V}}_{\mathit{U}}\end{array}\right]=\left[\begin{array}{c}{\mathit{d}}^{\left(\mathbf{1}\right)}\\ {\mathit{d}}^{\left(\mathbf{2}\right)}\\ \mathbf{0}\end{array}\right]$$

(5)

where C is the incidence-like matrix related to the linear transformation in Equation (4). The solution of the regularized problem in Equation (5) is, thus, eventually obtained in the least-squares (LS) sense. Once obtained, the Up-Down and East-West velocity vectors are integrated with respect to the time T by assuming deformation signals related to the overall first SAR acquisition at time $T_0$ are zero. As a quality index of the reconstruction, which takes mostly into account the congruence between obtained 2D displacement time-series and the original LOS measurements, the root mean square error (RMSE) of the residuals of the subsystem in Equation (3), namely $\sum ={\Vert B·V-d\Vert }_2$, is calculated (see the experiments shown in Section 4), with $d={\left[{\mathit{d}}^{\left(\mathbf{1}\right)},{\mathit{d}}^{\left(\mathbf{2}\right)}\right]}^T$ and $V={\left[{\mathit{V}}_{\mathit{E}},{\mathit{V}}_{\mathit{U}}\right]}^T$. Once the velocity estimates are obtained, the relevant 2D displacement time-series are simply recovered by time-integration. A comprehensive description of the MinA algorithm, also including its extension for the retrieval of the north-south displacement components, as well as for the integration of additional measurements (e.g., pixel-offset-driven measurements recovered from sequences of amplitude SAR images [48]), can be found in [19].

4. Experimental Results

In this section, we present the experiments performed for detecting and characterizing the deformation phenomena that occur in the coastal area of Shanghai by combining CSK and S1-A displacement time-series.

4.1. SBAS Analysis

We applied the SBAS algorithm to the available CSK and S1-A SAR datasets shown in Section 3.1. For the CSK datasets, a group of 155 InSAR data pairs were selected by imposing constraints on the maximum spatial and temporal baseline separation between orbits, which have been set, in this case, to 800 m and 1000 days, respectively. Concerning the S1-A dataset, a group of 368 InSAR data pairs were selected. As evidenced by the distribution of S1-A SAR images in the temporal/perpendicular baseline plane, which is shown in Figure 2, the S1-A system is characterized by a very short orbital tube, with a nominal diameter of about 200 m, i.e., much smaller than the relevant critical baseline value [44]. Therefore, S1-A InSAR data pairs were selected by simply imposing a temporal baseline less than one year. Differential interferograms were generated by computing the phase difference between the co-registered SAR image pairs, and by subtracting the relevant topographic phase contributions, as synthesized using the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) (with a spatial sampling of 90 × 90 m) and the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEM (30 × 30 m) of the area, for the CSK and the S1-A SAR datasets, respectively. To mitigate the effects of decorrelation noise into the generated interferograms, we also conducted a complex multi-look operation (with 30 looks in the azimuth and 30 looks in the range direction for the CSK dataset, and five looks in the azimuth and 20 looks in the range direction for the S1-A dataset) [41], so as to end up with a radar pixel spacing for both SAR datasets of about 100 × 100 m. The differential interferograms were noise-filtered [45,49], and the relevant spatial coherence maps were generated.

By inspecting the generated maps of spatial coherence, the group of coherent pixels common to all interferograms was identified, and the subsequent phase unwrapping steps [50] were eventually performed on this sparse grid of coherent pixels. By applying the SBAS technique, the relevant CSK and S1-A LOS-projected displacement time-series were recovered.

Since the LOS measurements related to each SAR dataset correspond to different radar pixel coordinates, the CSK and S1-A data products were resampled to a common geocoded grid, so that the subsequent InSAR combination technique could be applied to the CSK and S1-A SBAS displacement time-series relevant to homologous ground targets. Maps showing the LOS-projected mean displacement velocities for both datasets, superimposed on a Landsat image of the investigated area, are shown in Figure 4. LOS displacement values are calibrated with respect to a common reference point (its location is identified by the red star shown in Figure 4), whose deformation is assumed zero. To ensure reliable results, we only displayed the deformation values relevant to highly-coherent targets on the ground, which have been selected for the CSK dataset by identifying the points characterized by a temporal coherence value [51] greater than 0.65, and for the S1-A dataset by showing only the pixels with an average spatial coherence greater than 0.55.

4.2. Two-Dimensional (2D) Deformation Time-Series Extraction

To retrieve the time-series of the two-dimensional (2D) surface displacement components, we firstly extracted from the achieved SBAS time-series the CSK and S1-A LOS displacement time-series related to the overlapped time periods (i.e., from 28 February 2015 to 18 March 2016 for the CSK dataset, and from 26 February 2015 to 4 March 2016 for the S1-A dataset). Then, we re-calculated the LOS mean deformation velocity and the corresponding time-series over the 2015–2016 overlapped time-periods (i.e., from 26 February 2015 to 4 March 2016). Maps of the 2015–2016 SBAS LOS mean deformation velocities as detected by processing the CSK and S1-A datasets are shown in Figure 5a,b, respectively. Extracted LOS displacement time-series were then combined by employing the MinA algorithm [19]. In particular, a regularization factor $\mathsf{\alpha }=0.1$ was used.

As a result, the east-west and up-down displacement time-series were recovered. Figure 6 portrays the maps of the 2D displacement velocity components. As evidenced, most of the Shanghai area was almost stable, with the only exception of the ocean-reclaimed platforms down the coast. This finding represents one the main achievements of our study, and, to the best of our knowledge, it represents the first direct evidence of the actual vertical movements of the ocean-reclaimed platforms documented through DInSAR-based approaches. Furthermore, in order to quantify the correctness of the achieved 2D InSAR products, we also produced, for each point of the common CSK/S1-A geocoded grid, the Root Mean Square Error (RMSE) of the residuals between the original LOS-projected displacement time-series and the ones reconstructed from the retrieved east-west and up-down time-series. The map of the RMSE is portrayed in Figure 7. A zoomed view of the east-west and up-down mean displacement maps over the reclaimed areas is shown in Figure 8.

5. Discussion

Let us now focus on the analysis of the achieved 2D displacement time-series. Figure 9 shows the plots of the east-west and up-down time-series of deformation relevant to eight selected pixels located along the coast. Pixels were selected among those with a low value of RMSE of the residuals between the original and the reconstructed LOS time-series (see Figure 7). The analysis of Figure 9 reveals that reconstructed east-west deformation signals over the ocean-reclaimed land are (on average) less prevalent with respect to the up-down deformation signals recorded in homologous locations. The detected vertical movements are the effects on the ground of the consolidation mechanisms of the soils used for the reclamation of land from the ocean, and have also been investigated in previous studies [28,32]. This outcome is in general agreement with the results of studies [27,31], based upon considerations of the geology of the area and on the mechanisms of soil consolidation in reclaiming land. A discussion on the deformation signals detected over the Lingang New City and the Pudong International Airport is provided in the following subsections.

5.1. Lingang New City

Lingang New City, near the Dishui Lake, is characterized by a quite unique deformation scenario. Its western sector, as already anticipated by previous studies [28,32], based on the use of ASAR/ENVISAT SAR data (from 2007 to 2010) and COSMO-SkyMed data (from 2014 to 2016), was almost stable from 2015 to 2016 with maximum up-down and east-west deformation rates of a few millimeters per year (see for instance the plots relevant to points (e) and (f) shown in Figure 9).

Conversely, the eastern sector of the city was still affected during the 2015–2016 period by significant up-down residual deformations, reaching along the coast values up to 25–30 mm/year (see the plots relevant to pixels (g) and (h) shown in Figure 9). This finding is in general agreement with the investigations performed in [28,32].

5.2. Pudong International Airport

The achieved results, which are portrayed in Figure 10, reveal that actually most of the deformation affecting the airport is vertical (i.e., it is a subsiding area, as expected). In Figure 10, six points labeled (a–f) have also been selected as representative of the deformations of the zones where they are located. Relevant two-dimensional east-west and up-down deformation time-series are plotted in Figure 11. Our results complement and update the analyses presented in [31], based on processing—with a modified persistent scatters interferometry (PSI) technique—a set of high-resolution TerraSAR-X data satellite acquired from September 2011 to October 2012. With respect to [31] our updated results reveal that the soil consolidation process in correspondence to runway lines R1–R3 were still active in 2016, but with rates less than 10 mm/year, while the runway line R4, built over an area that was later reclaimed (starting from 2006), was severely affected by vertical residual subsidence deformations with rates up to 30 mm/year. In general agreement with [31], the present investigation shows that the deformation of the first runway foundations was almost completed by March 2016, and only a very tiny deformation subsidence (with a maximum rate of 5 mm/year) was still visible (with respect to the maximum rate of 15 mm/year measured in [31]) in 2016. See, for instance, the plots of east-west and up-down time-series related to point (a) shown in Figure 11. These results are in accordance with the geological setting of the zone reported in [31], considering that the first runway is built on a typical soft soil foundation with complete soil layers from top to bottom, and is the most stable of the five runways. The deformation related to the foundations of the first and third runways is quite similar. The foundations of the third runway were completed nine years later than the first runway. However, point (b) is located in the southern area of the third runway, where the foundation is mainly composed of silty clay and sandy silt, and point (a) is located in the middle of the first runway, where the foundation mainly consists of silty clay. Previous research shows that the water content, void ratio, and coefficient of compressibility of the silty clay in this area are 29.3%, 0.84, and 0.18 MPa⁻¹ respectively, while the parameters of the sandy silt are 24.2%, 0.72, and 0.25 MPa⁻¹, implying that the silty clay is softer and weaker than sandy silt [52]. Hence, the third runway is characterized by moderate-to-low deformation signals, which, from 2015 to 2016, were mostly vertical with a maximum subsidence rate of about 10 mm/year. See the plots related to point (b) shown in Figure 11. Conversely, the second runway, which has been in operation since 2005 and whose foundations reclamation project started later, only in 2001 [38], was still subject from 2015 to 2016 to residual deformations, but with maximum subsidence deformation rates less than those observed in [31] on the order of 10 mm/year. See the time-series of east-west and up-down deformation component related to point (c) shown in Figure 11. This outcome was rather expected due to the more reduced consolidation time of the second runway, as well as to its softer geology [31], formed by blown sand and backfill, with the addition of partial backfilled fine-silty sand [39]. The fourth runway, which was only completed in 2015, was markedly affected by vertical deformation signals with rates on the order of −25 to −35 mm/year and no appreciable east-west displacements, as documented by the plots of east-west and up-down displacement time-series related to point (d) that are shown in Figure 11. Concerning the fifth runway, which is the most proximal to the coast and the last runway built over reclaimed lands, its deformations were unfortunately not fully detectable. This is probably due to the temporal decorrelation noise artifacts that corrupt the interferograms, as certified by the reduced average spatial coherence of its nearby area, and it represents indirect proof of the circumstance that on-going deformation phenomena are very active.

However, in the near proximity of R5 detected pixels evidenced high deformation rates greater than –25 mm/year (see the up-down deformation time-series related to point (e) of Figure 11). The availability of additional data in the next years will help in better quantifying the entity of deformation near the fifth runway. Nonetheless, other structures (see the area marked with red points in Figure 10) were coherent, and the relevant deformations have been well detected. In this case, the measured deformation is predominantly vertical (i.e., a subsidence signal) with rates also larger than 40 mm/year, with the exception of the points located just along the coastline, which exhibit a more complex behavior and are also characterized by a residual movement towards the west in the northern sector and towards east in the southern sector. This behavior can be explained taking into account that the land-reclamation procedures for expanding the Pudong Airport facilities were in progress at the time of our investigation, and the soil was not compacted. For example, see the east-west and up-down deformation time-series relevant to the point labeled as (f) in Figure 11.

5.3. Additional Remarks

To better emphasize the relationships between the measured displacement rates and the construction time of the airport runways, we have performed this additional analysis, the results of which are shown in Figure 12. In particular, by identifying the well-detected, coherent points around each runway, we have calculated and displayed the average values of the east-west and up-down mean displacement velocity over the different runways. Note that although the third runway was put into use nine years after the first runway, the first and third runways both have smaller values of average east-west displacement and subsidence rate, of about 2 mm/year and 3.85 mm/year, respectively as shown in Figure 12a,b. This is because the foundations of the first and third runways are composed of similar natural sediments. Moreover, the foundation of the two runways (R1 and R3) was completed for the runway operation in 1999. After more than a decade of sediment consolidation, the residual deformations from 2015 to 2016 are small, and the foundation is almost consolidated and stable. Conversely, residual ground deformation is more significant in correspondence to facilities built on ocean-reclaimed land. This is, in particular, the case of the second (R2) and fourth (R4) runways, which were built on reclaimed-land platforms and completed at a later time (in 2005 and 2015, respectively). Indeed, they are characterized by larger measured displacement velocity values. This finding was somewhat expected as the soil consolidation phase was responsible for large deformation signals during the first stages after the end of the reclamation procedures [28]. With regard to the fifth runway (R5), the reclaiming process and the construction time of this runway occurred at the latest point, and (as expected) the measured deformation velocity values are the largest ones, reaching average subsidence rates and average east-west displacement rates of about 35 mm/year and 4 mm/year as shown in Figure 12a,b, respectively.

6. Conclusions and Future Perspectives

In this work, we have performed an up-to-date investigation of deformations of Shanghai coastal areas, and in particular of Lingang New City and of the infrastructures of the Pudong International Airport. The multi-angle/multi-sensor DInSAR approach, known as Minimum Acceleration (MinA) combination and based on the joint exploitation of LOS-projected deformation time-series obtained by applying the SBAS technique [6], has been used. The presented results have permitted the retrieval of the time-evolution of the on-going deformation phenomena by discriminating the east-west and the up-down displacement time-series, with a particular focus on the area of Pudong International Airport. The main outcome of our investigation is that ocean-reclaimed land is subject mostly to subsidence, while lateral (east-west) movements are less significant. In general, over ocean-reclaimed land east-west deformation rates are one-tenth of the up-down deformation rates. To the best of our knowledge, this represents the first direct evidence of the negligibility of lateral movements in reclaimed land documented through DInSAR-based approaches. As a future perspective, we plan to extend the presented analyses in the next few years to corroborate the validity of the foreseen deformation models discussed in [28,32]. An important point of criticism, however, needs to be raised: the European Space Agency acquisition plan of Sentinel-1A SAR data over the Shanghai area through descending orbits (which combined with the ascending data could allow a continuous mapping of the 2D displacements of the imaged areas) is not regular as in the ascending S1-A dataset. In fact, at the time this investigation was performed (i.e., in the second half of the year 2017), only a few isolated Interferometric Wide Swath (IW) descending SAR data (six images) were available in the ESA Hub catalogue. This circumstance evidently represents a critical issue to be taken seriously into account in the perspective of future investigations. This problem can be partially overcome by using SAR data acquired through other multiple-orbits SAR platforms (e.g., COSMO-SkyMed, TerraSAR-X, ALOS-2), but at the expense of a worse temporal resolution of achievable results and higher costs.