A Procedure to Map Subsidence at the Regional Scale Using the Persistent Scatterer Interferometry ( PSI ) Technique

In this paper, we present a procedure to map subsidence at the regional scale by means of persistent scatterer interferometry (PSI). Subsidence analysis is usually restricted to plain areas and where the presence of this phenomenon is already known. The proposed procedure allows a fast identification of subsidences in large and hilly-mountainous areas. The test area is the Tuscany region, in Central Italy, where several areas are affected by natural and anthropogenic subsidence and where PSI data acquired by the Envisat satellite are available both in ascending and descending orbit. The procedure consists of the definition of the vertical and horizontal components of the deformation measured by satellite at first, then of the calculation of the “real” displacement direction, so that mainly vertical deformations can be individuated and mapped.


Introduction
Subsidence can be defined as the progressive lowering of the ground; it can be a rapid process, such as a sudden collapse (i.e., sinkholes), due to karst processes, or a slow process; in the latter case, this phenomenon can be due both to natural causes, such as active faults (e.g., [1]) or volcanic processes, and anthropic causes, such as settlements, induced by new buildings [2], mining [3] and fluid pumping from underground [4,5].
Subsidence can be successfully analyzed by means of the PSI (persistent scatterer interferometry) technique.
The PSI technique has shown its capability to provide information about ground deformations over wide areas with millimetric precision, making this approach suitable for both regional- [21] and local-scale [22] subsidence investigations.Through a statistical analysis of the signals backscattered from a network of individual, phase coherent targets, this approach retrieves estimates of the displacements occurring between different acquisitions by distinguishing the phase shift related to ground motions from the phase component, due to the atmosphere, topography and noise [23,24].
The application of the PSI technique to the analysis of subsidence is somehow facilitated by the fast data processing.Furthermore, subsidence often occurs in urban areas, where the high density of radar targets allows high accuracy measurements.Several applications of PSI to subsidence mapping can be found in the literature, such as [25][26][27][28].
The analysis of these phenomena using PSI is usually performed in plain areas and where their presence is known [6,29], leading to a lack of knowledge about the natural hazards affecting whole regions.
Data acquired by various kinds of instruments are suitable for subsidence mapping (e.g., GPS [30] or leveling [31]); however, most of them are useful only to study local phenomena, and they are usually not available at the regional scale.
This limitation has been overcome using radar satellite data, which allow measurements over large areas and, under suitable conditions, spatially dense information on slow ground surface deformations.
At the same time, the PSI technique provides high accuracy that ranges from 1 to 3 mm on single measurements in correspondence to each SAR acquisition and between 0.1 and 1 mm/y for the line of sight (LOS) deformation rate [24].
The proposed method has been tested in Tuscany, where several areas affected by subsidence have been recognized and analyzed by a number of authors (i.e., [6,29,32]), but nowadays, a unique archive of the subsidences present in the regional territory still does not exist.
A complete review of the existing works highlighted the lack of information for most part of the region, since almost all of the examined works were focused on the phenomena affecting the main alluvial plains of the region, where the major urban areas of Tuscany are located.
To perform a mapping of the subsidence phenomena of an entire region, the main problem is to discriminate areas where ground movements can be due to subsidence or to other causes (i.e., landslides), since, at this scale, the morphology of the territory may vary significantly and the assumption that subsidences occur only in plain areas cannot and should not be considered valid (e.g., [33,34]).Overcoming this problem requires data and an analysis procedure adequate to the scale of the work.
This work wants to provide a procedure to efficiently map subsidence at the regional scale using the PSI technique.The test area is the Tuscany region in Central Italy where several areas of different extensions have been affected by natural and anthropogenic subsidence.The procedure makes use of ascending and descending datasets of the Envisat satellite ranging in time from 2002 to 2010.This procedure can be reproducible also in other test sites and also with different datasets of PSI data, given the availability of ascending and descending data.

Study Area
The proposed procedure has been tested in the Tuscany region (Central Italy), which is about 23,000 km 2 wide.The territory is mainly hilly (66.5%), but it includes also mountainous areas (25.1%) and limited plains (8.4%) [35].
The geology and morphology of this region vary greatly from zone to zone: the main reliefs lie in the northern and eastern part of the region.
The northern part is characterized by high mountains (up to 2000 m a.s.l.) made up of metamorphic rocks (Apuan Alps) and by steep valleys filled up by alluvial deposits; the eastern part is characterized by mountains made up of sedimentary rocks and by intermontane basins filled up by alluvial deposits, where the main cities of the region are located.
The central and southern parts are characterized by hilly morphology with isolated volcanic reliefs (e.g., Monte Amiata) and spread alluvial plains (Figure 1).

Methodology
By the use of the PSI technique, a simple procedure to perform rapid mapping of subsidence phenomena has been developed.This methodology has been applied and tested in the Tuscany region (Italy), where PSI data acquired by the Envisat Satellite are available from 2002 to 2010, both in ascending (41 images, 1,142,707 PS) and descending orbit (41 images, 1,353,746 PS).
The PSI data have been processed following the PSInSAR™ approach.
The acquisition of the SAR satellite occurs along two different polar orbits, descending from north to south and ascending from south to north.Displacements are measured along a unit vector codirectional to the satellite, defined as the line of sight.Being that the orbit of SAR satellites is polar, it is impossible to estimate the component of displacement along the N-S direction on the azimuth plane.
Data acquired in both orbits are necessary to perform subsidence analysis by PSI data over areas with different terrain altitudes, because the deformations detected by this technique are measured along the line of sight of the satellite, so the real value of the deformation highly depends on the orientation of the terrain with respect to the satellite.
In fact, if deformations are detected in plain areas, landslide phenomena can be easily excluded, but this is not true if deformations are detected in hilly-mountainous areas, leading to possible misinterpretations of results.
If data acquired in both orbits are available, this problem can be overcome, since the velocity (or the rate) of vertical and horizontal (only in E-W direction) deformations can be calculated by solving the following equation system [26,36]: where Va and Vd represent the velocity of deformation in ascending and in descending orbit, VV and VE the velocity in the vertical and horizontal (E-W) planes and θa and θd the incident angles (both 23°) in ascending and descending orbit (Figure 2).The solution of this system requires that a single measurement point is recognized as a valid target from the satellite in both orbits, but this condition rarely occurs, so artificial persistent scatterers (called synthetic PS), have to be created.This operation can be performed dividing the space through a sampling grid with square cells and calculating the mean deformation velocity for each of the two orbits; the results of this operation are two spatially regular series of synthetic PS, each one corresponding to the centroid of every cell, which can be now combined to define the components of the deformation vector.
This procedure is shown in the flow chart reported in Figure 3: a sampling grid is applied to the same scene, for both orbits, and the mean deformation velocity is calculated for each cell.Now, these velocities (Va and Vd) can be used to calculate the vertical (VV) and horizontal (VE) velocities, by Equation (1).The definition of the cell size is not so obvious, since it should be sized appropriately.Having too few points within each cell should be avoided to use in the calculation, but also, points too far from each other, which could measure deformations due to different causes, should be avoided.Furthermore, using a mesh with too big of a dimension, it could be possible to merge the phenomena of the same type, but that should be considered as separate from a geological point of view (e.g., subsidences affecting two distinct river basins could be mapped as one).
The scale of the work [37] is another factor to be considered when dimensioning the grid: using small cells (e.g., 1 hm 2 ) for regional scale analyses could be excessively expensive, in terms of time and computational capability.
On the other side, using small cells can be useful for a detailed analysis of local phenomena, where high detail is necessary.
If in local analyses the cell sizes can be easily defined considering the range of autocorrelation given by the semivariogram of PS, in case of regional scale works, the experience of the operator plays a fundamental role, since a unique geometric criterion is not identifiable, but other factors, mainly the geology and geomorphology of the area, have to be considered.
Once the vertical and horizontal components of the deformation vector are known (positive velocities mean uplift or eastward deformations; negative velocities mean subsidence or westward deformations), it is possible to calculate the main direction of ground deformation (Vr), so as to discriminate areas affected by subsidence or by landslides.
The modulus of Vr can be easily calculated by Pythagoras' theorem, since the angle between Vv and VE is 90°.
The direction of Vr has been computed by means of trigonometric rules (Figure 4), at first solving the following equation (Equation 2) for both α and β to define the angle between VV (or VE) and Vr: where α is the angle between VE and Vr, β is the angle between VV and Vr, δ is the angle between VV and VE and it is 90° (sin δ = 1).Subsequently, the signs of VV and VE moduli have been used to define the direction of Vr.Solving Equation (1), upward or eastward movements will have positive moduli, and conversely westward or downward movements will have negative moduli.These values have been used to define in which quarter Vr will be located (Table 1).
Table 1.Definition of the Vr direction on the basis of VV and VE values.

V V and V E Values
Quarter both V V and V E positive Z-E quarter V V positive and V E negative Z-W quarter both V V and V E negative N-W quarter V V negative and V E positive N-E quarter In order to discriminate areas affected by subsidence or by landslides the areas with mainly vertical deformations have been selected, filtering all the Vr values which are within a range of ±45° to the vertical axis (Figure 5), since subsidences (or uplifts) can be considered as mainly vertical movements.The filtered values could be now used to individuate the main subsidence areas (Figure 6), but it must be borne in mind that the radar satellites, used to perform PS analysis, cannot evaluate N-S-oriented displacements [24].This is to mean that northward (or southward) landslides may look like subsidences, since satellites can detect only the vertical component of the displacement.Comparison between unfiltered and filtered PS.This image is an example of the result of a regional-scale analysis, performed by a grid with 1-km 2 cells.It is possible to notice the presence of a subsiding area (red box).

Figure 7.
Comparison between the aspect map and synthetic PS for the subsiding area highlighted in Figure 6.The analysis of deformation has been refined by the use of a grid with 1-hm 2 cells.Zone 1: two PS seem to show only vertical deformations, but they are in a northward slope, so horizontal deformations cannot be excluded.Zone 2: Almost all PS show only vertical deformations, regardless of the aspect of the slope; a subsidence can be recognized.A 100-m side sampling grid has been used.
To avoid misunderstandings of the results, Vr (or VV and VE) values can be compared with a map showing the orientation of the slopes (namely an aspect map), easily obtainable from a DTM of the study area (in this case, the DTM is 10 × 10 m).
This comparison allows excluding these values, which fall within the northern or southern slopes (Figure 7).
It is worth noticing that this comparison is especially useful for subsiding areas that are not too large (few km 2 ), in hilly-mountainous areas: in these areas, slope angles can vary from low (~15°) to moderate (~30°), and misunderstandings between landslide and subsidence can be frequent; consequently, the definition of the Vr direction can be useful to avoid misinterpretations.
In large subsiding areas, a comparison with the aspect map would not be necessary, since a large number of PS, consistent with each other, is expected, regardless of the morphology of the territory.

Discussion
The presented procedure allows for quickly identifying and mapping subsidence phenomena at all scales, from regional to local ones.The main advantage of this procedure is to perform subsidence mapping in areas characterized by high topographic variability, where misunderstandings of the results can be easy.
Since each synthetic PS is representative of a certain area, which may have different dimensions on the basis of the scale, this procedure can be more powerful when used in successive phases, i.e., using a sampling grid with broad cells (e.g., 1 km 2 ) to individuate the main subsiding areas, followed by a new analysis performed by a grid with smaller cells (e.g., 1 hm 2 ) for these areas.
In this way, it is possible to initially identify the areas affected by subsidence and then to perform an analysis and a precise mapping of the phenomena, taking advantage of the previous small-scale analysis.
Obviously, each one of these phases can be singly applied, on the basis of the goal of the work, i.e., if a subsidence phenomenon is well known, it would be a waste of time to perform a double analysis to locate and successively map the subsidence.
This procedure has some limitations and uncertainties related to the PSI technique.
The main limitation of this procedure is due to the impossibility of radar satellites to measure N-S-oriented deformation: in such a case, no analysis can be performed.Some uncertainties can be present when the analysis is performed in hilly-mountainous areas, where subsidences and landslides could be mixed up, but this problem can be overcome if a sufficient number of PS, with concordant velocities and directions of deformation, is available; the majority of these PS should be located on slopes not oriented in the N-S direction, otherwise, merely their high number would be useless for identifying the cause of the deformation.
A further limitation is the availability of data, since if the satellite data are not available, a long time is required for their acquisition.Similarly, if data acquired only in one orbit (ascending or descending) are available, the whole procedure cannot be performed.
Similarly, if few PS targets are available, the synthetic PS may not be representative of its cell, since it would be calculated only with few PS data (e.g., 2-3 targets for 1 km 2 ).The minimum number of required PS can be defined a priori, since it can vary on the basis of the scale of the work and the extension and morphology of the study area.
In this work, we used Envisat data, but this procedure can be applied also to data acquired by other SAR satellites, such as TerraSAR-X or COSMO-SkyMed.The temporal resolution of different satellites can influence the results of the procedure: for instance, data acquired by satellites with a shorter revisiting time (e.g., COSMO-SkyMed) can be used to analyze faster subsidence phenomena than those identified by the use of Envisat data.
The spatial resolution of the satellite (e.g., Envisat 5 × 20 m) must be considered when dimensioning the grid, to avoid useless analyses, such as, for instance, defining the cell size as smaller than the spatial resolution of the satellite.Furthermore, the ground resolution of the satellite (pixel size) influences the mapping precision when a detailed analysis is performed, while it has a lower influence in regional scale analyses.

Conclusions
Subsidence mapping and analysis by PSI are usually performed in areas where this phenomenon is already known.Furthermore, these analyses are rarely performed in hilly or mountainous areas.
This paper presents a procedure to identify and map subsidence using the PSInSAR technique at the regional scale and/or in areas characterized by high topographic variability, where, usually, the identification of subsidences could be difficult.
The discrimination of subsidence from landslides could be a useful tool for spatial planning and for risk management strategies, since good knowledge of the natural hazards involving a territory is essential to providing proper risk reduction measures.
The procedure is based on the comparison and combination of interferometric data, acquired in different satellite orbits, so as to discriminate subsiding areas from other areas.
The first step of the procedure is the calculation of the vertical and horizontal (in E-W direction) components of the deformation.This operation requires the definition of synthetic PSs (for both orbits), which can be made by dividing the space trough a sampling grid and assigning to each cell the mean deformation velocity of each PS.Once synthetics PS is defined, the two components can be calculated, and consequently, the main direction of ground deformation can be defined.
Until now, the traditional methods of combining ascending and descending orbits allow one to obtain two different maps related to VE and VV, which have to be compared to identify subsiding areas (e.g., [24]).The proposed approach goes further, since it allows one to calculate the displacement direction (Vr direction) from the values of VE and VV in order to have a unique parameter to distinguish subsidence from other ground deformations.
This discrimination is based on the principle that subsidences are characterized mainly (or only) by vertical deformations, so that, filtering out all of the deformations with the main horizontal components, the subsidences can be isolated and mapped.
The proposed technique is quite flexible, since it can be applied in various geological and geomorphological settings, and it can be used with different datasets of PSI.The procedure works at the regional scale, and it can be applied to the fast detection of subsidence phenomena for land planning and risk management purposes.

Figure 1 .
Figure 1.Location of the test area.

Figure 3 .
Figure 3. Flow chart of the procedure to calculate the vertical and E-W components of the deformation recorded by the satellite.PS, persistent scatterers.Vel, velocity.

Figure 4 .
Figure 4. Schema illustrating the angles between Vr, VV and VE.α is the angle between VE and Vr; β is the angle between VV and Vr; δ is the angle between VV and VE.

Figure 5 .
Figure 5. Filtering of the data, on the basis of the synthetic PS.The direction of Vr is calculated on the basis of the values of VV and VE (positive or negative).

Figure 6 .
Figure 6.Comparison between unfiltered and filtered PS.This image is an example of the result of a regional-scale analysis, performed by a grid with 1-km 2 cells.It is possible to notice the presence of a subsiding area (red box).