Combined Use of Cand X-Band SAR Data for Subsidence Monitoring in an Urban Area

In this study, we present the detection and characterization of ground displacements in the urban area of Pisa (Central Italy) using Interferometric Synthetic Aperture Radar (InSAR) products. Thirty RADARSAT-2 and twenty-nine COSMO-SkyMed images have been analyzed with the Small BAseline Subset (SBAS) algorithm, in order to quantify the ground subsidence and its temporal evolution in the three-year time interval from 2011 to 2014. A borehole database was reclassified in stratigraphical and geotechnical homogeneous units, providing the geological background needed for the local scale analysis of the recorded displacements. Moreover, the interferometric outputs were compared with the last 30 years’ urban evolution of selected parts of the city. Two deformation patterns were recorded by the InSAR data: very slow vertical movements within the defined stability threshold (±2.5 mm/yr) and areas with subsidence rates down to −5 to −7 mm/yr, associated with high peak velocities (−15 to −20 mm/yr) registered by single buildings or small groups of buildings. Some of these structures are used to demonstrate that the high subsidence rates are related to the recent urbanization, which is the trigger for the accelerated consolidation process of highly compressible layers. Finally, this urban area was a valuable test site for demonstrating the different results of the Cand X-band data processing, in terms of the density of points and the quality of the time series of deformation.


Introduction
Subsidence is a broad term usually referring to the ground lowering induced by natural or anthropic factors.Among the anthropic causes, rapid urbanization is considered as one of the most important accelerating factors for the consolidation process, especially where unconsolidated alluvial deposits, peat-rich layers, and reclamation areas are present [1][2][3][4][5][6].
In this study, C-band RADARSAT-2 (RSAT) and X-band COSMO-SkyMed (CSK) data have been processed with the Small BAseline Subset (SBAS) approach, to provide the deformation maps of the urban area of Pisa (Central Italy) in the time interval 2011-2014.The area of interest has been characterized by exploiting the potential of the two satellites, at a regional scale for RADARSAT-2 and at a building scale for COSMO-SkyMed.New deformation-affected areas and trend changes in the time series of deformation of selected buildings have been examined, to demonstrate the existing relation between recorded displacements, local geological setting, and urban development.

Geological and Hydrogeological Setting
The Arno coastal plain, which is 450 km 2 wide, is a densely-populated area with several industrial activities and diffuse cultural heritages.The city of Pisa occupies a central position in this coastal plain, which represents the southern inshore portion of the Viareggio basin, a subsiding half-graben [36,37]; this tectonically formed basin was progressively infilled by marine and alluvial sediments after the Last Glacial Maximum (LGM) [36,37].The Quaternary sea level fluctuations produced an incised valley perpendicular to the present shoreline, with a maximum depth of 40 m [38][39][40].After the eustatic sea level rise post LGM, estuarine and coastal plain sedimentation took place, filling almost half of the valley [37].The last 20-25 m of the infilling sequence is composed of transgressive Holocene sediments, deposited after the phase of the decelerated sea level rise [41].
At its base, this succession has a level of clays and silty clays, locally known as "pancone", continuously widespread in the plain to the border of the coastal beach ridges system.This level, deposited in a lagoonal environment, is characterized by the presence of layers rich in organic matter [42].Above the pancone, a fluvio-deltaic sequence, 10-15 m thick, is found; it is composed of clay and silt with sporadic sand and silty sand bodies that reflect the coastal progradation of the modern Arno delta [37].
The stratigraphic structure of the Holocene succession leads to the presence of geotechnically-weak layers in the first 20 m of the underground beneath the city [43], producing a preferential predisposing factor for ground subsidence and the structural failure of buildings, as testified by the most famous cultural heritage of the city, the Leaning Tower [44].
From a hydrogeological point of view, a phreatic aquifer composed of sandy and silty-sand layers within the clayey sequence is found.This water reservoir is not usually used for water extraction because of its low transmissivity and quality, and its high vulnerability to contamination [45].The phreatic surface level is almost the same across the northern part of the city and in the city center, while the minimum depth (2 m below sea level) is registered in the southern part of the International Airport Galileo Galilei (Figure 1).The maximum level is found in the Ospedaletto district, reaching 1.5 m above sea level [45].The most important water resource for the city of Pisa is located at a higher depth, ranging from 50 to 100 m below the surface, and it is composed of a confined aquifer made up of Pleistocene gravels and sands, 10-20 m thick [46,47].This aquifer, generally characterized by water of a low quality, is mainly used for gardening, irrigation, and industry, while only a few pumping wells are exploited for drinking water [46].

Figure 1.
Geographical framework of the area of interest and the phreatic surface in the urban area of Pisa (from [45]).The black line represents the urban area.The contour map is overlaid on a 2010 digital orthophoto.

Interferometric Data
The two SAR datasets exploited for this study are analyzed with the advanced DInSAR technique, referred to as the SBAS approach [48][49][50][51][52][53], to retrieve the deformation time series in the temporal range between 2011 and 2014.The data processing of the SAR images has been performed by using the full Resolution SBAS-DInSAR processing chain, developed and owned by CNR-IREA (National Research Council-Institute for Electromagnetic Sensing of the Environment, Napoli, Italy).
The key feature of this algorithm's chain is the selection of the SAR images that will be used to generate the interferograms; this is done by choosing only those data pairs characterized by short temporal and perpendicular baselines [48].This procedure allows for a minimization of the noise effects, intended as decorrelation phenomena [52], thus maximizing the number of coherent pixels in the generated interferograms.Independent acquisition datasets, separated by large baselines, are properly linked by applying the Singular Value Decomposition (SVD) method, to retrieve deformation time series and velocity maps.The SBAS algorithm is also capable of producing results at two different scales [49]: a regional-scale, for investigating diffused deformation phenomena at a medium spatial resolution (the cell resolution ranges from 100 × 100 m to 30 × 30 m) [51]; and a local scale, which exploits the full resolution interferograms (benefiting from the high-resolution capability of the SAR satellite system) for generating products ideally suited for building analysis [53,54].The regional-scale product is used to estimate the low pass signal components, typical of the atmospheric artifacts and related to deformations over large areas [48].Conversely, the local-scale output is exploited to investigate the high frequencies, typical of localized displacements [49,51,53].
The COSMO-SkyMed X-band (λ = 3.1 cm) dataset is comprised of 29 ascending acquisitions collected between the 26/09/2012 and the 11/05/2014.The data have been acquired in Stripmap mode

Interferometric Data
The two SAR datasets exploited for this study are analyzed with the advanced DInSAR technique, referred to as the SBAS approach [48][49][50][51][52][53], to retrieve the deformation time series in the temporal range between 2011 and 2014.The data processing of the SAR images has been performed by using the full Resolution SBAS-DInSAR processing chain, developed and owned by CNR-IREA (National Research Council-Institute for Electromagnetic Sensing of the Environment, Napoli, Italy).The key feature of this algorithm's chain is the selection of the SAR images that will be used to generate the interferograms; this is done by choosing only those data pairs characterized by short temporal and perpendicular baselines [48].This procedure allows for a minimization of the noise effects, intended as decorrelation phenomena [52], thus maximizing the number of coherent pixels in the generated interferograms.Independent acquisition datasets, separated by large baselines, are properly linked by applying the Singular Value Decomposition (SVD) method, to retrieve deformation time series and velocity maps.The SBAS algorithm is also capable of producing results at two different scales [49]: a regional-scale, for investigating diffused deformation phenomena at a medium spatial resolution (the cell resolution ranges from 100 × 100 m to 30 × 30 m) [51]; and a local scale, which exploits the full resolution interferograms (benefiting from the high-resolution capability of the SAR satellite system) for generating products ideally suited for building analysis [53,54].The regional-scale product is used to estimate the low pass signal components, typical of the atmospheric artifacts and related to deformations over large areas [48].Conversely, the local-scale output is exploited to investigate the high frequencies, typical of localized displacements [49,51,53].
The COSMO-SkyMed X-band (λ = 3.1 cm) dataset is comprised of 29 ascending acquisitions collected between the 26 September 2012 and the 11 May 2014.The data have been acquired in Stripmap mode with HH polarization.The CSK scenes are characterized by an incidence angle at mid-range of around 30 • and by a spatial resolution of 3 m in both range and azimuth directions.The 29 CSK scenes have been firstly converted from the raw data L0 format into the Level1 SLC (Single Look Complex) format, and a quality assessment has also been applied to these images.Then, they have been co-registered with respect to the 11 July 2013 master image, chosen as reference geometry.The CSK dataset is characterized by a large spatial baseline distribution, since the orbital tube exceeds 2000 m and an average perpendicular baseline of 380 m, as well as by a quite dense temporal sequence, with an average sampling of 21 days (Figure 2a).The CSK images have been coupled, by applying the Delaunay triangulation method [55], into 88 interferometric pairs, characterized by maximum spatial and temporal baseline values of about 800m and 500 days, respectively.Following the interferogram generation step, the PhU (Phase Unwrapping) operation at the regional scale has been carried out by applying the method described in [55] and selecting the pixels that have a coherence, estimated in a box of 10 × 10 pixels, in the range and azimuth directions of greater than 0.35 in at least 30% of the interferograms [55].The final selection of the pixels correctly unwrapped has been carried out by considering pixels with a temporal coherence value [52] greater than 0.7.
The RADARSAT-2 C-band (λ = 5.6 cm) dataset is composed of 30 ascending orbit images acquired in Standard 3 (S3) beam mode (one of eight standard beam modes of the sensor, covering a nominal stripmap swath of 108 km and acquired with a range bandwidth of 11.56 MHz) and with HH polarization, covering the time interval 27 October 2011-26 April 2014.The RSAT scenes are characterized by a mid-range incidence angle of around 30 • and a spatial resolution of 12 × 5 m in the slant range and azimuth directions, respectively.The 30 RSAT images, in the Level1 SLC format, have been properly co-registered with respect to a common geometry (master image); the 14 November 2012 acquisition has been selected for this purpose.The RSAT dataset is characterized by a good spatial baseline distribution, with an orbital tube around 700 m and an average perpendicular baseline of 120 m, in addition to a temporal sequence with an average sampling of 31 days (Figure 2b).The RSAT images have been coupled into 87 interferometric pairs, characterized by maximum spatial and temporal baseline values of about 300 m and 900 days, respectively.Following the interferogram generation step, the PhU operation at the regional scale has been carried out by selecting the pixels that have a coherence, estimated in a box of 4 × 20 pixels, in the range and azimuth directions of greater than 0.35 in at least the 30% of the interferograms [55].The final selection of the pixels correctly unwrapped has been carried out by considering pixels with a temporal coherence value [52] greater than 0.7.
The results of the data processing of both satellite datasets show a very high coverage for the urban area of Pisa, with a lack of measurement points in only some areas, near the city limits, which are largely covered by crop fields (Figure 3).The two datasets show a marked difference in terms of the density of points; in fact, the density of the measurement points within the area of interest (33 km 2 ) is 910 points/km 2 for RSAT and 11,027 points/km 2 for CSK.This great increase in the number of points between the C-and X-band is consistent with the results obtained from other test sites using CSK data [14,15,51,56,57].The higher number of measurement points-i.e., coherent points-is achieved thanks to the lower revisit time and higher spatial resolution of the X-band sensors.Moreover, the shorter wavelength allows an increase in the sensitivity to LOS (Line of Sight) displacements, together with a higher capability of detecting very low displacement rates [14,15,51,56,57].The derived displacements are projected along the line ideally joining the sensor and ground target (LOS), mainly considering vertical (up-down) and East-West movements [58,59].Each coherent point is associated with a displacement value expressed in millimeters per year and a sign depending on the movement; positive if an uplift is registered (blue color in the deformation map), and negative if ground subsidence is affecting the point (red color in the deformation map).

Stratigraphical Information and Borehole Classification System
For better understanding the recorded displacements, detailed stratigraphical information about the subsurface setting is needed.This information has been extracted from the Pisa Province borehole database, that includes both the lithological descriptions and geotechnical parameters obtained from laboratory tests.Ninety-two boreholes in the urban area of Pisa were selected according to the unique restriction of being deeper than 10-15 m, in order to intercept at least the upper limit of the pancone deposit.For homogenizing the subsurface information from a geotechnical point of view, a classification system was applied.This classification, proposed by Sarti et al. [43] and recently

Stratigraphical Information and Borehole Classification System
For better understanding the recorded displacements, detailed stratigraphical information about the subsurface setting is needed.This information has been extracted from the Pisa Province borehole database, that includes both the lithological descriptions and geotechnical parameters obtained from laboratory tests.Ninety-two boreholes in the urban area of Pisa were selected according to the unique restriction of being deeper than 10-15 m, in order to intercept at least the upper limit of the pancone deposit.For homogenizing the subsurface information from a geotechnical point of view, a classification system was applied.This classification, proposed by Sarti et al. [43] and recently modified by Solari et al. [6,60], is composed of five facies associations.A detailed discussion on the characteristics of each class can be found in Solari et al. [6,60].Accordingly, we will summarize the class descriptions, focusing on the geotechnical implications.
Unit 1 consists of the pancone very soft clay and silt clay level.The thickness of this unit varies from 3 to 12 m and can be found at a depth of 8 to 20 m, depending on the position in the plain.From a geotechnical point of view, it is considered as a highly compressible level, with a high water content that is sometimes close to the liquid limit.These characteristics strongly depend on the brackish sedimentation environment in which the unit was formed, which also favored the development of thin organic layers that further decrease the general resistance of this unit.
Unit 2 consists of dark soft clay and silty clay with a high organic matter content and occasional peat layers.It is typically found on top of Unit 1, with a thickness ranging from 2 to 6 m.This level has similar characteristics to the previous one, with a low strength and high compressibility, sometimes improved by the presence of less organic layers with a higher consistency.
Unit 3 consists of a rhythmical alternation of clayey silt, sandy silt, and silty sand, with a discontinuous distribution in the plain; where found, it lays at 2-6 m below surface, with a thickness smaller than 4 m.The geotechnical parameters, such as the water content and plasticity index, show a high variability, reflecting the chaotic structure of this sequence in which the fine sand percentage fluctuates.This unit is characterized by better mechanical properties than the previous two; however, it must be considered a potentially compressible layer because of its degree of saturation.In fact, if saturated conditions occur, even these fine-grained sands and silts can produce a high loss in volume.
Unit 4 consists of stiff clay with a very low organic matter content and high consistency due to a subaerial depositional setting of the sediments that produced resistant pedogenetic surfaces.Unit 4 is the uppermost level in the sequence and has a thickness of 5-6 m.This unit has the highest consistency and shear strength among all of the clayey layers.
Unit 5 consists of fluvial channel fine to coarse sand, with a high variability in the thickness (1-10 m) and vertical position within the stratigraphical sequence (lenticular bodies can be found at the contact with the pancone level).In terms of its mechanical behavior, it is considered as the unit with the lowest settlement potential.
The potential settlements in the area are produced not only by the presence of Unit 1 alone, but are related to the simultaneous presence of more compressible layers.The worst-case scenario consists of this bottom-up stratigraphy: Units 1 and 2, overlapped by Unit 4, and partially replaced by a saturated Unit 3 level.This sequence, common in the area, is considered as the most susceptible to settlements [6].

Regional Scale Analysis
The analysis of the two A-DInSAR datasets shows that the urban area of Pisa is not characterized by the presence of regional scale deformation patterns, such as large subsidence bowls related to water overexploitation (Figure 3).This type of large scale deformation has been observed in other sectors of the Tuscany region, such as the urban area of Pistoia, where a water overexploitation-related subsidence was recorded [61], or the geothermal area of Larderello, where a large subsidence bowl, produced by geothermal activity, has been detected [62].The general trend of ground subsidence increases in the northern part of Pisa, where the city has developed outside of the historic center walls.In this portion of the city, a general lowering of 3.5-5 mm/yr is documented.On the other hand, the highest subsidence rates are only recorded in localized areas, corresponding to a single building or groups of buildings.The presence of these areas with high localized subsidence rates, coupled with phreatic surface reconstruction, which do not show any area of intense groundwater level lowering (Figure 1), suggests the exclusion of water overexploitation as a triggering factor for the measured displacements.Moreover, a possible role played by the second aquifer compaction was further excluded because of its high depth (>50 m below surface).
A stability threshold of ±2.5 mm/yr, which is an estimation of the accuracy of the measurements [54], has been defined on the basis of the standard deviation of the calculated velocity values of each point of the two datasets.This stability range is also coherent with various literature examples in the field of C-and X-band Interferometric Synthetic Aperture Radar (InSAR) data analysis, i.e., [12,14,16,19].Using this threshold, two areas have been identified as stable: the Galileo Galilei International Airport (Figure 3 point A) and the historic city center, where the greatest part of the Pisa cultural heritage is present (Figure 3, point D).The airport structures and the runaways are stable in the investigated time-period, showing mean subsidence rates of −0.3 mm/yr for RSAT and −0.8 mm/yr for CSK, respectively.Within the area delimited by the ancient city walls, the LOS velocities are very low, with average values of −1.4 mm/yr and −1.2 mm/yr for RSAT and CSK, respectively.As previously mentioned, the northeastern portion of the city shows a general subsidence of 3.5-5 mm/yr in both InSAR datasets, with high localized subsidence rates.For example, in the northern part of the La Fontina district (Figure 1), ground lowering of −13 mm/yr for RSAT and −12 mm/yr for CSK datasets have been recorded by commercial buildings.In the San Cataldo district (Figures 1  and 3-point C), the peak subsidence rates measured are in agreement with those of a university residence for students, with mean values of −13 mm/yr for both sensors.In the Cisanello district (Figures 1 and 3-point C), peak velocities of −12 mm/yr for both sensors have been recorded by a group of residential buildings.Another area where localized high subsidence rates have been recorded is the commercial district recently built outside of the western border of the International airport (Figures 1 and 3-point A); values of −28 mm/yr for CSK data have been detected.The RSAT data show a very low density of points in this area; for this reason, displacement information cannot be derived.In the Ospedaletto industrial and commercial district (Figures 1 and 3-point B), the greatest part of the buildings is within the defined stability threshold (±2.5 mm/yr), but the peak velocities can also be identified, with values of −12 mm/yr for the RSAT dataset and −14 mm/yr for the CSK dataset.
A stability threshold of ±2.5 mm/yr, which is an estimation of the accuracy of the measurements [54], has been defined on the basis of the standard deviation of the calculated velocity values of each point of the two datasets.This stability range is also coherent with various literature examples in the field of C-and X-band Interferometric Synthetic Aperture Radar (InSAR) data analysis, i.e., [12,14,16,19].Using this threshold, two areas have been identified as stable: the Galileo Galilei International Airport (Figure 3 point A) and the historic city center, where the greatest part of the Pisa cultural heritage is present (Figure 3, point D).The airport structures and the runaways are stable in the investigated time-period, showing mean subsidence rates of −0.3 mm/yr for RSAT and −0.8 mm/yr for CSK, respectively.Within the area delimited by the ancient city walls, the LOS velocities are very low, with average values of −1.4 mm/yr and −1.2 mm/yr for RSAT and CSK, respectively.As previously mentioned, the northeastern portion of the city shows a general subsidence of 3.5-5 mm/yr in both InSAR datasets, with high localized subsidence rates.For example, in the northern part of the La Fontina district (Figure 1), ground lowering of −13 mm/yr for RSAT and −12 mm/yr for CSK datasets have been recorded by commercial buildings.In the San Cataldo district (Figure 1 and Figure 3-point C), the peak subsidence rates measured are in agreement with those of a university residence for students, with mean values of −13 mm/yr for both sensors.In the Cisanello district (Figure 1 and Figure 3-point C), peak velocities of −12 mm/yr for both sensors have been recorded by a group of residential buildings.Another area where localized high subsidence rates have been recorded is the commercial district recently built outside of the western border of the International airport (Figure 1 and Figure 3-point A); values of −28 mm/yr for CSK data have been detected.The RSAT data show a very low density of points in this area; for this reason, displacement information cannot be derived.In the Ospedaletto industrial and commercial district (Figure 1 and Figure 3-point B), the greatest part of the buildings is within the defined stability threshold (±2.5 mm/yr), but the peak velocities can also be identified, with values of −12 mm/yr for the RSAT dataset and −14 mm/yr for the CSK dataset.

Movement Detection of Recently Built Buildings
The exploited datasets have allowed us to map active deformation areas, with subsidence rates higher than 10 mm/yr affecting the buildings which have been built in the last 10 years.Examples of these deformations can be found in the San Cataldo district and in the commercial area near the International Airport, where shopping malls and shipyards are present (Figure 3 points A, C).
In particular, the I Praticelli complex, a residential structure for students inaugurated in 2008, is one valuable example of a recently built building that displays high deformation rates (Figure 4a).Both SAR datasets registered similar subsidence rates; the RADARSAT-2 data ranged between −11 and −18 mm/yr with a mean value of −13.3 mm/yr for the entire complex, whereas the COSMO-SkyMed data ranged between −6 and −17 mm/yr with a mean value of −12.3 mm/yr (Figure 4b).It is interesting to note that the maximum subsidence rates for both datasets are registered in the southern part of the complex; this can be related to foundation problems or very localized variations of the geotechnical characteristics of the compressible levels that cannot be defined with the available subsurface data.This case study also shows the different performances of C-and X-band data in urban areas; in fact, the X-band data have a density of points for this single structure that is six times larger than the C-band data, passing from 75 to 500 points.Only three boreholes, localized at a mean distance of 200 m from the I Praticelli complex, are available for the geological reconstruction of the subsurface; this is reasonable considering that no significant stratigraphical variations can occur in this distance thanks to the almost homogeneous distribution of the main units.The upper limit of the most compressible succession (lagoonal clay and organic clay, Units 1 and 2) lies at a variable depth, from 4 to 5 m below the surface.This sequence, which has a total thickness of 20-25 m, is overlaid by a thin consolidated clay level (Unit 4).As mentioned by Solari et al. [6], this type of sequence is the most susceptible to settlement.The registered settlements are directly related to the natural evolution of the consolidation process, as defined by Terzaghi and Peck [63].The end of this time-dependent process is highly modified by the layers' thickness, hydraulic conductivity, and drainage conditions, as well as by the time and magnitude of the load imposition.Accordingly, the high deformation rates registered by this building are therefore related to a triggering anthropogenic factor (the recently applied load) and to a geological predisposing factor, i.e., the presence of shallow highly compressible levels (Units 1 The registered settlements are directly related to the natural evolution of the consolidation process, as defined by Terzaghi and Peck [63].The end of this time-dependent process is highly modified by the layers' thickness, hydraulic conductivity, and drainage conditions, as well as by the time and magnitude of the load imposition.Accordingly, the high deformation rates registered by this building are therefore related to a triggering anthropogenic factor (the recently applied load) and to a geological predisposing factor, i.e., the presence of shallow highly compressible levels (Units 1 and 2).Moreover, foundation design problems cannot be excluded.

Buildings with Trend Changes
The RADARSAT-2 and COSMO-SkyMed datasets have been exploited to detect possible trend changes in the temporal deformation patterns registered by buildings, in the San Cataldo and Ospedaletto districts (Figure 1), built in the eighties-early nineties.A complete time series from 1992 is available for some of these buildings, thanks to the exploitation of data from the ERS 1/2 and Envisat archives, already analyzed by Solari et al. [6].
An example of integration between historical SAR datasets and the new C-and X-band data is the analysis of the displacement rates recorded in correspondence of the CNR building complex.This three-floor construction, built at the end of the eighties, occupies a total area of 30,000 m 2 .A high subsidence rate was recorded for the ERS 1/2 and Envisat C-band datasets, with a maximum magnitude of −13.4 mm/yr in the period covered by the ERS 1/2 data (1992-2000) [6].A decrease in magnitude of the mean LOS velocity was then registered, up to the value of −9 mm/yr for the period 2003-2010, covered by the Envisat data [6].The RADARSAT-2 data confirms this decreasing subsidence rate, recording a mean magnitude of −6 mm/yr with a maximum of −8.7 mm/yr in the central part of the complex (Figure 5a).A similar trend, with a comparable magnitude (mean velocity of about −6 mm/yr) and spatial distribution to the highest velocity values, is recorded by the COSMO-SkyMed data (Figure 5a).
The stratigraphical sequence below the analyzed building, with the presence of the Unit 1 and 2 compressible layers at a mean depth of 5m and with the only resistant Unit 4 level partially substituted by a saturated silty level (Unit 3), is again highly prone to ground settlements (Figure 5b).According to this setting, the deformation rates registered in the first part of the time series are directly related to compaction, due to the applied load of the highly compressible layers present in the sub-surface of the city (Figure 5c).The decrease of LOS velocity recorded by the CNR building represents the natural evolution of the consolidation process that slows down over time [63].As observed by Peduto et al. [64], a decrease in the registered velocities can be related to external anthropic factors, such as restoration works to the foundation of buildings.Moreover, a deeper analysis can be performed to characterize the building settlements in further detail, through the implementation of specific indexes of structural health for each building (see for example [65,66]).
In the Ospedaletto commercial district, a different behavior of the structures has been observed in the monitored period.Several deformations, that reach velocities higher than −10 mm/yr, have been detected by ERS 1/2 and Envisat InSAR data, covering the period 1992-2010 [6].All of the eastern part of this area has been subjected to recent urban expansion, with the construction of numerous commercial buildings since the end of the 1990s.Many of the structures built in the first part of this period, which recorded the above mentioned high deformation rates, appear now as stable or near to the defined stability threshold (±2.5 mm/yr).One example of this behavior is building A (Figure 6a); as shown by the proposed time series (Figure 6c), the average subsidence rates drastically decrease from 12.3 mm/yr for the period covered by the Envisat data (2003-2010), to 1.4 mm/yr in the RADARSAT-2 analysis and 1.0 mm/yr for the COSMO-SkyMed dataset.A stratigraphical explanation for this behavior must be considered: the most compressible sequence has its upper boundary deeper than in the previous example, passing from 5 to 10 m below surface.This ground evidence implies that the compressible sequence is less influenced by the stress increment produced by the commercial building loads, generating a lower compaction of the compressible strata.On the other hand, if we consider structures which have been built very recently, the recorded subsidence rates reach 15 mm/yr, as shown in Figure 6a by building B. One example of this behavior is building A (Figure 6a); as shown by the proposed time series (Figure 6c), the average subsidence rates drastically decrease from 12.3 mm/yr for the period covered by the Envisat data (2003-2010), to 1.4 mm/yr in the RADARSAT-2 analysis and 1.0 mm/yr for the COSMO-SkyMed dataset.A stratigraphical explanation for this behavior must be considered: the most compressible sequence has its upper boundary deeper than in the previous example, passing from 5 to 10 m below surface.This ground evidence implies that the compressible sequence is less influenced by the stress increment produced by the commercial building loads, generating a lower compaction of the compressible strata.On the other hand, if we consider structures which have been built very recently, the recorded subsidence rates reach 15 mm/yr, as shown in Figure 6a

Relation between InSAR Displacements and Urbanization
For a smaller scale use of the displacement information retrieved for the single buildings, the LOS velocities have been averaged in the area occupied by 500 different structures of selected districts, in which the highest peak velocities and subsidence rates have been recorded (San Cataldo, Cisanello, La Fontina and Ospedaletto districts, see Figure 1).
In Figure 7A, the C-and X-band InSAR-derived LOS velocity values are plotted versus the age of construction of the selected buildings (Pre 1978, 1978-1988, 1988-2000, Post 2000).Four temporal classes have been defined on the basis of multi-temporal ortophotographical information.The most recent class does not present any ERS 1/2 data because the end of the acquisition period of the sensor (December 2000), that from 2001, even if still operational, does not provide reliable acquisitions because it has been operating without gyroscopes.As expected, the data show an increase in the mean LOS velocity values starting from the oldest temporal class (Figure 7B).The RSAT and CSK

Relation between InSAR Displacements and Urbanization
For a smaller scale use of the displacement information retrieved for the single buildings, the LOS velocities have been averaged in the area occupied by 500 different structures of selected districts, in which the highest peak velocities and subsidence rates have been recorded (San Cataldo, Cisanello, La Fontina and Ospedaletto districts, see Figure 1).
In Figure 7a, the C-and X-band InSAR-derived LOS velocity values are plotted versus the age of construction of the selected buildings (Pre 1978, 1978-1988, 1988-2000, Post 2000).Four temporal classes have been defined on the basis of multi-temporal ortophotographical information.The most recent class does not present any ERS 1/2 data because the end of the acquisition period of the sensor (December 2000), that from 2001, even if still operational, does not provide reliable acquisitions because it has been operating without gyroscopes.As expected, the data show an increase in the mean LOS velocity values starting from the oldest temporal class (Figure 7b).The RSAT and CSK data show a general lower dispersion of the points for the three oldest classes (see variance values in Figure 7b); this is related not only to the general decrease in the LOS velocities of the two central classes (1978-1988 and 1988-2000), but also to the higher density of measurement points for each building.In fact, the higher redundancy of points in these two datasets allows for the minimization of the effects of anomalous velocity values within the building perimeter.The high variance (Figure 7b) of the points constituting the Post 2000 class is mainly related to the presence of outliers characterized by high subsidence rates; these points correspond to the peak velocities recorded in the analyzed districts for the studied time-periods.These evidences are coherent with the InSAR-derived deformations obtained in the city of Florence (located 80km eastern than Pisa) by Pratesi et al. [66,67], where a similar stratigraphical context and urban evolution in the last 30 years is present.As demonstrated by the I Praticelli complex example, these high subsidence rates are directly correlated to compaction, due to the recently applied load, of the highly compressible levels present in the stratigraphic setting of the city at a shallow depth.Moreover, the general decrease in the mean velocities and variance values in the two central classes demonstrate that the subsidence in the urban area is slowing down.This is testified, for example, by the CNR complex in the San Cataldo district and by many structures in the commercial area of Ospedaletto.
Considering the statistical parameters, the dispersion of the oldest class (Pre 1978) points suggests the presence of a natural subsidence, with values ranging from −2.5 to 5 mm/yr; this is more evident if only the RSAT and CSK are considered (Figure 7a).In this class, points with a displacement magnitude above 5 mm/yr are outliers related to buildings affected by the ongoing consolidation process due to geological (clayey levels with high organic matter content) or constructive (bad foundation design) characteristics.
Geosciences 2017, 7, 21 13 of 17 data show a general lower dispersion of the points for the three oldest classes (see variance values in Figure 7B); this is related not only to the general decrease in the LOS velocities of the two central classes (1978-1988 and 1988-2000), but also to the higher density of measurement points for each building.In fact, the higher redundancy of points in these two datasets allows for the minimization of the effects of anomalous velocity values within the building perimeter.The high variance (Figure 7B) of the points constituting the Post 2000 class is mainly related to the presence of outliers characterized by high subsidence rates; these points correspond to the peak velocities recorded in the analyzed districts for the studied time-periods.These evidences are coherent with the InSAR-derived deformations obtained in the city of Florence (located 80km eastern than Pisa) by Pratesi et al. [66,67], where a similar stratigraphical context and urban evolution in the last 30 years is present.As demonstrated by the I Praticelli complex example, these high subsidence rates are directly correlated to compaction, due to the recently applied load, of the highly compressible levels present in the stratigraphic setting of the city at a shallow depth.Moreover, the general decrease in the mean velocities and variance values in the two central classes demonstrate that the subsidence in the urban area is slowing down.This is testified, for example, by the CNR complex in the San Cataldo district and by many structures in the commercial area of Ospedaletto.
Considering the statistical parameters, the dispersion of the oldest class (Pre 1978) points suggests the presence of a natural subsidence, with values ranging from −2.5 to 5 mm/yr; this is more evident if only the RSAT and CSK are considered (Figure 7A).In this class, points with a displacement magnitude above 5 mm/yr are outliers related to buildings affected by the ongoing consolidation process due to geological (clayey levels with high organic matter content) or constructive (bad foundation design) characteristics.

Conclusions
In this paper, A-DInSAR data from different sensors operating in C-and X-band have been exploited, to describe the temporal and spatial evolution of the land subsidence in the urban area of Pisa.The ground movements were related to a predisposing factor (presence of highly compressible layers of clays and organic clays at shallow depth), and to a triggering factor (recent urbanization of part of the city).
The C-band RADARSAT-2 and the X-band COSMO-SkyMed data, analyzed with the SBAS algorithm, show that the highest deformation registered in the investigated time-period (2011-2014) affects single buildings in the recent urban expansion areas of the city.These buildings recorded subsidence rates higher than −15 mm/yr.Considering this localized subsidence, a building-scale analysis was carried out.In the San Cataldo and Ospedaletto districts, three examples have been provided to display the different behaviors of the selected buildings.In these areas, a widespread compressible sequence was detected, mainly composed of organic rich clay and silty clay classified in two geotechnical units (Unit 1 and 2).This sequence supports the development of an accelerated consolidation process in the case of a load application; this is the case of the I Praticelli complex, a recent structure that registered subsidence rates of up to 20 mm/yr.The CNR building and the commercial structures of the Ospedaletto district show, when also exploiting older C-band data (ERS 1/2 and Envisat), the slow-down in the subsidence phenomena during the analyzed period.Finally, using the averaged LOS velocities extracted from several buildings in four selected districts, the correlation between registered subsidence rates and the age of construction of buildings was confirmed and strengthened, highlighting the role of urbanization as a triggering factor for the registered rapid subsidence.
This work has also shown the potential of DInSAR information extracted from sensors of different bands in urban areas, providing data at both a regional and building scale, with a high density of points.

Figure 1 .
Figure 1.Geographical framework of the area of interest and the phreatic surface in the urban area of Pisa (from [45]).The black line represents the urban area.The contour map is overlaid on a 2010 digital orthophoto.

Figure 2 .
Figure 2. Synthetic Aperture Radar (SAR) data representation in the temporal/perpendicular baseline plane for the COSMO-SkyMed (a) and RADARSAT-2 (b) datasets.Each black triangle represents a SAR acquisition, with the code-name of the image, and the arcs that are the interferometric pairs selected for the Small BAseline Subset (SBAS) processing.

Figure 2 .
Figure 2. Synthetic Aperture Radar (SAR) data representation in the temporal/perpendicular baseline plane for the COSMO-SkyMed (a) and RADARSAT-2 (b) datasets.Each black triangle represents a SAR acquisition, with the code-name of the image, and the arcs that are the interferometric pairs selected for the Small BAseline Subset (SBAS) processing.

Figure 4 .
Figure 4. (a) LOS displacements affecting the analyzed building.The displacement map is overlaid on a 2013 digital orthophoto; (b) Geological reconstruction of the stratigraphical asset, the classification is the one proposed by Sarti et al. [43]; (c) RADARSAT-2 and COSMO-SkyMed average time series for the "I Praticelli" complex.

Figure 4 .
Figure 4. (a) LOS displacements affecting the analyzed building.The displacement map is overlaid on a 2013 digital orthophoto; (b) Geological reconstruction of the stratigraphical asset, the classification is the one proposed by Sarti et al. [43]; (c) RADARSAT-2 and COSMO-SkyMed average time series for the "I Praticelli" complex.

Figure 5 .
Figure 5. (a) LOS displacements affecting the National Research Council (CNR) complex.The displacement map is overlaid on a 2013 digital orthophoto; (b) Geological reconstruction of the stratigraphical asset, the classification is the one proposed by Sarti et al. [43]; (c) ERS 1/2, Envisat, RADARSAT-2 and COSMO-SkyMed average time series for the analyzed building.The ERS 1/2 and Envisat time series have been modified after Solari et al. [6].The RADARSAT-2 and COSMO-SkyMed time series are shifted by a constant value for a better visualization; (d) RADARSAT-2 and COSMO-SkyMed average time series for the investigated period.The y axis units refer to the displacement in millimeters.

Figure 5 .
Figure 5. (a) LOS displacements affecting the National Research Council (CNR) complex.The displacement map is overlaid on a 2013 digital orthophoto; (b) Geological reconstruction of the stratigraphical asset, the classification is the one proposed by Sarti et al. [43]; (c) ERS 1/2, Envisat, RADARSAT-2 and COSMO-SkyMed average time series for the analyzed building.The ERS 1/2 and Envisat time series have been modified after Solari et al. [6].The RADARSAT-2 and COSMO-SkyMed time series are shifted by a constant value for a better visualization; (d) RADARSAT-2 and COSMO-SkyMed average time series for the investigated period.The y axis units refer to the displacement in millimeters.
by building B.

Figure 6 .
Figure 6.(a) LOS displacements affecting the Ospedaletto district.The displacement map is overlaid on a 2013 digital orthophoto; (b) Geological reconstruction of the stratigraphical asset, the classification is the one proposed by Sarti et al. [43]; (c) Envisat, RADARSAT-2 and COSMO-SkyMed average time series for the analyzed building.The Envisat time series have been modified after Solari et al. [6].The RADARSAT-2 and COSMO-SkyMed time series are shifted by a constant value for a better visualization.

Figure 6 .
Figure 6.(a) LOS displacements affecting the Ospedaletto district.The displacement map is overlaid on a 2013 digital orthophoto; (b) Geological reconstruction of the stratigraphical asset, the classification is the one proposed by Sarti et al. [43]; (c) Envisat, RADARSAT-2 and COSMO-SkyMed average time series for the analyzed building.The Envisat time series have been modified after Solari et al. [6].The RADARSAT-2 and COSMO-SkyMed time series are shifted by a constant value for a better visualization.

Figure 7 .
Figure 7. (a) Age of construction of the selected buildings compared to the mean LOS velocity of the DInSAR data; (b) Mean, variance and standard deviation values for each class and each sensor.The ERS 1/2 and Envisat series are modified from Solari et al. [6].The vertical black line represents the defined stability threshold (ST).

Figure 7 .
Figure 7. (a) Age of construction of the selected buildings compared to the mean LOS velocity of the DInSAR data; (b) Mean, variance and standard deviation values for each class and each sensor.The ERS 1/2 and Envisat series are modified from Solari et al. [6].The vertical black line represents the defined stability threshold (ST).