1. Introduction
The Italian peninsula together with Sicily, Sardinia, and several minor islands, is one of the most representative European coastal countries with a coastline of more than 7000 km, which make it the sixth European region regarding coast extension. Also, more than 1000 rivers give rise to many sea mouths along the Adriatic, Ionian, and Ligurian–Tyrrhenian sides. Over the centuries, the coastal areas have represented essential outlets for trade and the economic development of the country. Even today, there are critical urban settlements located along the Italian coasts that are major centers for agriculture, fisheries, tourism, port trade, and industrial activities. However, from north to south, the Italian coastline is mostly affected by several natural and anthropogenic ground deformation phenomena [
1,
2,
3,
4], which can severely impact the environment, stability of the infrastructure, and safety of the people. Such phenomena, i.e., ground deformation and uplift, arises from both shallow (e.g., compaction of sediments and groundwater pumping) and deeper processes (e.g., gas/oil exploitation, volcanic and tectonic processes [
1,
4,
5,
6,
7]). On one side, the compaction of high-porosity near-surface sediments caused by groundwater overexploitation often dominates ground deformation, with rates up to tens of cm/year [
8,
9,
10,
11,
12]. On the other side, more in-depth activities contribute to coastal deformation with rates up to a few mm/year [
2,
11]. Moreover, several phenomena such as the inundation of lowlands, coastal erosions, the salinization of aquifers, storms, flooding, and sea level rise may contribute to increasing the ground deformation hazard along coastal areas [
5]. Therefore, in order to have a complete understanding of the coastal hazards, a proper knowledge of such phenomena, as well as investigating their potential relationship with the sea-level rise effect, is needed [
1,
2,
5].
In this context, direct observations from ground-based and remote sensing data are nowadays available for the analysis of most of the deformation phenomena occurring along the Italian coastlines, thus supporting risk prevention and mitigation [
13].
Over the past decades, synthetic aperture radar (SAR) has been proven as a powerful tool for the monitoring and the surveillance of coastal areas, with multi-temporal and multi-scale capabilities enabled by the concept of SAR constellations [
4,
14,
15,
16]. In this framework, advanced multi-temporal differential interferometric SAR (A-DInSAR) approaches were extensively and successfully exploited to investigate coastal ground deformation [
17,
18]. These techniques, such as Permanent Scatterers (PS) [
19], small baseline subsets (SBAS) [
20], and interferometric point target analysis (IPTA) [
21] use a medium-to-large dataset of SAR images acquired at different times over the observed area to follow the temporal evolution of a deformation phenomenon, retrieving the mean deformation rate and the time series for each point target. Moreover, they can take benefits from the massive availability of free-of-charge and commercial SAR data archives collected and provided by space agencies, e.g., the European Space Agency (ESA), Italian Space Agency (ASI), Japan Aerospace Exploration Agency (JAXA), Canadian Space Agency (CSA), and German Aerospace Center (DLR).
This paper is meant as a demonstrator to show the capability of InSAR data to provide a detailed understanding of the coastal hazards constraining ground deformation phenomena occurring along coastal areas. To these purposes, we present and discuss some meaningful InSAR-based case studies along the Italian coastline, where multi-temporal techniques were applied to X-band COSMO-SkyMed (CSK) and C-band ENVISAT SAR data to evaluate local natural and anthropogenic phenomena. We analyze the advantages and limits of the adopted sensors and techniques regarding spatial and temporal coverage, the density of targets, lack of data, etc.
Finally, we discuss the potentialities of Sentinel-1 (S1) missions to improve the coastal monitoring and sustainability. Indeed, S1 SAR data show some peculiarities that may be suitably exploited for performing a free-accessible near real-time surveillance service of coastal areas with SMART features (i.e., Sustainable, Measurable, Automatic, Realistic, and Time-oriented) to support decisions from public administration, government institutions, and civil protection entities regarding risk mitigation and prevention.
2. Methodology
In this study, we used IPTA Multi-baseline (IPTA-MB) and SBAS approaches implemented through the commercial GAMMA© and SARScape© software, respectively. The main differences between the two methods concern the rationale (step-by-step for IPTA or semi-automatic for SBAS) and the data format. Regarding products, the output of the two methods are equivalent, and there are no significative differences in the quality and accuracy of the derived deformation maps. According to the product provider expertise, one of the two approaches can be exploited.
In the IPTA approach, the temporal and spatial characteristics of SAR point targets are exploited to retrieve the mean surface deformation velocities and the related time series of displacement [
21]. The advantages and the benefits of the proposed technique are: (a) the step-by-step rationale that allows the user to accurately check and refine any processing step; (b) the use of point scatter targets to reduce geometrical decorrelation; (c) the potential to find scatterers in low coherence areas, thus filling the gaps in deformation maps and then resolving both high phase gradients and phase jumps; (d) the use of vector format data structures instead of the raster data format (typically used in conventional interferometry) that allows drastically reducing the disk space required for data storage; and (e) the possibility of obtaining deformation information, even though a reduced number of images is available.
The main processing steps of the GAMMA software package to implement the IPTA-MB technique are summarized below (
Figure 1): (a) assembling a stack of co-registered single-look complex (SLC) images with the related SLC and orbit parameters; (b) selecting a candidate list of point targets (plist) based on coherence and/or amplitude mean-to-sigma ratio thresholds; (c) selecting the interferometric maximum spatial and temporal baselines; (d) obtaining an initial estimate of interferometric parameters (e.g., deformation, topography, atmosphere, noise); (e) identifying the iterative improvement of different parameters based on the results of previous steps; and (f) finally, the provision of multi-temporal InSAR results, e.g., deformation maps and time series together with height and quality information by solving the system of deformations equations through the singular values decomposition (SVD) method. Further details about the IPTA technique can be found in Werner et al. [
21].
The SBAS approach is a multi-temporal InSAR approach that provides the mean velocity and time series of surface velocities at medium (about 80 × 80m) and full (about 10 × 10 m) spatial resolution [
20], depending on the adopted multi-look factors and native resolution of the SAR system. It preserves the full area coverage ensured by actual spaceborne SAR systems despite a coarse spatial resolution. The proposed technique provides a stack of InSAR pairs starting from coupled SAR images acquired along nearby orbits (small spatial baselines) with low revisit time (small temporal baselines), by setting the allowed maxima spatial/temporal baselines. This enables minimizing the spatial and temporal decorrelation issues affecting the InSAR pairs, and increasing the spatial density of the point measurements of deformation. Also, it aims at mitigating atmospheric artifacts and topographic or orbital errors by applying spatial/temporal filtering. The expected outputs are characterized by high measurement accuracy for both mean velocity (about 1–2 mm/year) and surface deformation (about 5 mm).
The main processing steps of the SARscape software to implement the SBAS technique are (
Figure 1): (a) the generation of the connection graph, which includes defining the reliable InSAR pair connections distributed in different subsets according to the user-defined temporal and spatial baselines’ maxima values; (b) the possibility of choosing an area of interest in the full SAR imagery; (c) the interferometric workflow, providing (i) the co-registration of the SAR stack list with respect to a “super master” SAR acquisition, (ii) the generation of interferograms, (iii) the phase noise filtering, and (iv) the phase unwrapping; (d) refinement and re-flattening to remove the residual phase ramps due to orbital uncertainties; (e) the first inversion step, providing (i) the velocity and residual topography estimation, (ii) the mean power image and precision component estimation, (iii) the second unwrapping phase; (f) the second inversion step, implementing (i) newly the refinement and re-flattening by means of a set of user-defined ground control points (GCPs), (ii) the displacement component estimation, (iii) the atmosphere pattern removal, and (iv) the coherence estimation; and (g) the geocoding of ground velocity and time series displacement. Further details of the SBAS technique are provided in Berardino et al. [
20].
3. Results
In this section, we show some local deformation phenomena detected by InSAR analysis in several Italian coastal areas effectively covered by SAR data. Six study areas are analyzed. They are the coastal zones that include and surround Ravenna, Fiumicino, Campi Flegrei, Sibari, Augusta, and Taranto (
Figure 2). We used C-band ENVISAT data provided by the European Space Agency (ESA) and X-band data acquired by the CSK missions of ASI. They are characterized by different resolutions, geometry of view, and spatial and temporal coverage, depending on the specific case study; thus, they provide an overview of the current and future scenarios of InSAR-based applications to the monitoring of coastal areas.
Figure 2 reports an overview map of the Italian peninsula with the case studies discussed in this paper, while
Table 1 summarizes the datasets, methods, and max deformation rates for each of the analyzed areas.
3.1. Ravenna
The coastal areas overlooking the Adriatic Sea in the proximity of Ravenna form a quaternary sedimentary basin bordered by the Northern Apennines and the Southern Alps [
22]. The complex seismotectonic framework of the observed area is interested in the compression and extension processes in the foreland and hinterland, respectively
, resulting in the actual seismicity distribution with many strong earthquakes [
23]. Moreover, the whole coastal area was interested in intensive groundwater pumping and hydrocarbon exploitation activities, which provided decreasing land settlement phenomena with subsidence rates on the order of magnitude of some cm/year [
4,
24].
The SAR dataset used for studying this area consisted of 97 images acquired in single-look complex (SLC) format by the C-band ENVISAT SAR mission along ascending and descending orbits in the time-span of 2003–2010. They provide a ground coverage of 100 × 100 km2 and a spatial resolution of 25 m × 5 m along the range and azimuth direction, respectively. The SBAS approach was used to retrieve the mean deformation velocity and time series. To ensure a reliable stack of interferometric pairs, temporal and spatial baselines were set to maximum values of 450 m and 700 days, respectively. A multi-look factor of 4 × 20 in the range and azimuth direction was chosen to reduce the speckle noise, adequately match the digital elevation model (DEM) (Shuttle Radar Topography Mission, SRTM 90m resolution) spatial resolution, and obtain a large density of points in each SAR image.
The results of multi-temporal InSAR processing are shown in
Figure 3. The Line-Of-Sight (LOS) mean velocity maps highlight the overall subsidence along both acquisition orbits in the whole time span 2003–2010. In detail, the mean negative deformation velocities ranged between 2.5
–5 mm/year, with maximum subsidence rates of about 12 mm/year. The latter can be observed both along the coastlines and the major urban settlements where anthropogenic activities are intensely concentrated.
3.2. Fiumicino
Fiumicino is the area located west of the city of Rome, in the proximity of the Tiber River delta plain. There, it is also hosted the main airport of the city, i.e., the international airport “Leonardo da Vinci”, which has more than 40 million passengers per year. Due to the soft alluvial sediments surrounding the Tiber River, many subsidence phenomena affect the area, making it potentially dangerous for civil infrastructures [
25,
26,
27].
To investigate the deformation phenomena occurring in this area, we used a dataset of 38 C-band ENVISAT SAR images collected from 2003 to 2010 along the descending orbit. Multi-look factors of two and 10 along the range and azimuth directions were applied, thus resulting in a pixel posting of about 40 m
× 40 m. The SRTM DEM was used for removing the topography. Goldstein filtering [
24] and unwrapping with minimum cost flow algorithm [
25] was also applied to each interferogram. We discarded the interferogram with unwrapping errors and large atmospheric artifacts, which can affect the results producing spikes and over and underestimation displacements in the InSAR time series. Then, by applying the IPTA-MB approach, the mean deformation velocity and time series were retrieved (
Figure 4). InSAR analysis reveals some of the signals of subsidence surrounding the Tiber Valley and in the proximity of the municipalities of Ostia and Fiumicino, which are close to the coastline (points P1 and P3).
The most significant subsidence was observed along the southeastern part of runway three of the Leonardo da Vinci airport (point P2), where a subsidence rate of more than 5 mm/year was recorded. The stratigraphy of subsoil showing different subsoil features explains the deformation gradient along the two sections of the runway [
26,
28].
3.3. Campi Flegrei
Campi Flegrei is a broad area affected by volcanic activity in the southern coast of Italy, close the city of Naples. It is a potentially dangerous site because of the complexity of the caldera and the density of the population. Over the centuries, it was affected by the phenomenon known as bradyseism, i.e., successions of great uplift and subsidence phases. In 2011–2015, it was characterized by a significant uplift phase with registered vertical uplift values up to 20 cm [
29]. In this study, we exploited the large SAR dataset acquired by the COSMO-SkyMed (CSK) constellation along the ascending and descending orbits in order to constrain the 2011–2015 Campi Flegrei uplift. The descending dataset consisted of 56 SAR single-look complex (SLC) images acquired between February 2011 and March 2015, whereas the ascending one was composed of 264 images spanning from January 2011 to September 2015. The datasets were processed with the IPTA-MB approach. The topographic phase contribution was removed by the 90-m Digital Elevation Model provided by the SRTM. In total, 681 and 283 interferograms were retrieved for the ascending and the descending orbit, respectively. The results regarding the deformation rate and time series are shown in
Figure 5. The ascending and descending data highlight a clear uplift signal surrounding the Solfatara area starting from 2011.
As shown by the analysis of time series, the uplift accelerated in 2013; then, it slowed down in 2014 up to the beginning of 2015, when it seemed to be going toward a new and fast uplift phase (points P1 and P3). Concerning the negative signal on the ascending data corresponding to the area of Bagnoli–Fuorigrotta (point P2), it is consistent with a significant eastward displacement component going away from the satellite due to the volcano inflation.
3.4. Sibari
We analyzed the surface movements occurring in the Sibari Plain (hereinafter SP). The SP is located in the northeastern part of the Calabria Region (Southern Italy) and extends over an area of ~500 km
2 [
30]. The SP has an important economic factor thanks to its agriculture, and there is also a cultural enticement due to the presence of the ancient Sybaris, which was a famous Greek colony during the Magna Grecia period. The studied site is also interesting due to its flooding hazard, even if it is reduced by a drainage channels network. Subsidence phenomena have been known in the Sybaris coast and its archaeological area since the 1950s. More than 20 cm were observed between 2003–2010, which was ascribed to groundwater exploitation and primary terrain consolidation.
We adopted the SBAS approach to a SAR image dataset acquired from the ESA ENVISAT satellite along the ascending track. This dataset (track 86, frame 789) contained 38 images in the interval 2003–2010, resulting in 135 pairs. In the processing chain, temporally decorrelated atmospheric artifacts were estimated using double-pass filtering in the temporal and spatial domain, respectively. The SRTM DEM (90 m) was selected to subtract the topography contribution. Precise orbit positions were considered to remove orbital residuals, and ground control points (GCPs) on stable areas were chosen mainly along the frame borders to remove residual ramps.
The surface deformation field was typified by common subsidence with a rate of up to ~15 mm/year (
Figure 6). Anthropic activities can be considered an additional factor for coastal erosion and subsidence phenomena. Indeed, higher velocities were observed in the more urbanized areas, where urban-induced loading magnified soil compaction. Some areas were characterized by extreme water withdrawals (e.g., Sybaris archaeological site); water table changes may increase the mean displacement velocity, which could be the reason for the set back consolidation process in soft materials.
3.5. Augusta
Here, we present the results of a multi-temporal InSAR analysis performed using SAR images from the ENVISAT sensor and covering the time interval 2003–2010. The study has the objective of identifying active ground movements of tectonic and anthropogenic origin in Augusta Bay in southeastern Sicily [
31,
32].
In this study, 33 and 55 SLC acquisitions were selected in the temporal spans of 27 August 2003—11 November 2009 and 26 February 2003—22 September 2010 for the descending and ascending orbit, respectively. We obtained evident subsidence in the Augusta graben (
Figure 7). The latter is a structural basin extended in the northwest (NW) direction and developed in the last two to three million years. The basin is enclosed to the northeast (NE) by the Mt. Tauro Horst, and detached from the River Marcellino graben to the southwest (SW) by the Petraro Horst. The retrieved pattern (
Figure 7) shows an elongated shape with the major axis along the NW–southeast (SE) direction, which is parallel to the maximum graben extent. The affected area shows an edged border along the eastern rim, which is characterized by the narrow escarpment of the Mt. Tauro fault, and a more gradual deformation gradient along its western limit. Generally, the subsidence phenomenon is ~8.5 km long (NW–SE) and ~4 km wide, with a mean velocity of –6 mm/yr. Moreover, large displacements occur along an NW–SE alignment, parallel to the Mt. Tauro fault, with the highest rates (up to –18 mm/year) close to the recent urbanization of Augusta city, built over redeemed saltworks. The displacement time series mainly shows a linear lowering trend reaching a total amount of ~70 mm and ~40 mm for the ascending and descending cases, respectively [
33]. The lowest value for the descending orbit is caused by the shorter considered period, while the minor linearity (especially in the last period) depends on the smaller number of SAR acquisitions in the InSAR processing.
3.6. Taranto
A dataset of ENVISAT images from the ascending orbit has been processed to retrieve the mean ground velocity of the Ionian coastal area of the Puglia region, around Taranto city (
Figure 8). It is composed of 38 acquisitions, spanning a period from May 2003 to September 2010, just before the ENVISAT deorbiting operation was done extending the lifetime of the ESA mission. The dataset was processed by applying the IPTA-MB method of GAMMA software. In particular, we calculated the interferograms by applying a multi-look factor of 1 × 5 in the range and azimuth directions, respectively, which leads to a ground pixel coverage of 20 m × 20 m. Then, only coherent points (i.e., persistent scatterers) characterized by the high signal-to-noise-ratio of the interferometric phase were selected. The selection of the interferometric pairs was made by setting the maxima values of the spatial and temporal baselines. Such values were 280 m and 1300 days, respectively, which returned a stack of 107 point-wise interferograms used for the SVD inversion.
The area did not show ongoing deformation processes (
Figure 8). There was no evidence of anthropogenic-induced displacement, such as subsidence due to water exploitation or building loading on soft soils, even though there were industrial activities in this region. Neither tectonic deformation nor significant sea-level rise changes seemed to be acting in this region [
34].
4. Discussion
InSAR LOS displacement measurements enable quantitatively analyzing the deformation processes along the Italian coastlines in the last decades, over small-to-large areas, and for short-to-long temporal spans. Moreover, these measurements can identify areas affected by intense ground deformation phenomena that can be strongly exposed to coastal erosion, inundation, and flooding events. The availability of multi-temporal InSAR data represents valid support to improve the knowledge of both slow and fast ground deformation phenomena, ensuring a synoptic view of the whole Italian coasts.
In this work, we showed several ground deformation phenomena detected by X and C-band InSAR data along the Italian coastline due to the complex interplay of natural and anthropogenic activities.
In detail, the Ravenna coastline and Po Plain were affected by land subsidence related to natural soil compaction and consolidation, as well as groundwater pumping activities and gas/oil extraction from on-shore and offshore wells.
The subsidence observed near Rome, in the area surrounding the Tiber River and Fiumicino town, are induced by soil compaction effects that are mostly present close the rivers due to soft alluvial sediments. Particularly compelling is the case of the “Leonardo da Vinci” airport, where a transition between the natural sediments, ranging from clay levels to medium-to-coarse sand, leads to a deformation gradient observed along runway three.
Moving to the south of Italy, we observed a different phenomenon in the Campi Flegrei area, west of Napoli city. This area is affected mainly by the bradyseism phenomenon, i.e., the succession of several inflation/deflation phases because of the presence of a nested volcanic district. Here, we showed, by InSAR analysis, the significant uplifting phase that occurred in the time interval 2011–2015 due to the magma recharging in the magma chamber of the central volcano of the area, e.g., the Solfatara of Pozzuoli.
Going further south, we observed non-negligible ground deformation along the Sibari Plain, along with the Ionic Coast of the Calabria region, due to groundwater overexploitation at nearby nautical centers, residential centers, and industrial zones. Looking at the southernmost Calabrian coastlines and the Sicily region, the soil compaction of sediments due to groundwater exploitation represents the primary cause for the observed coastal ground deformation, which was also observed in the case of Augusta Bay along the eastern coast of Sicily.
Concerning the Puglia coastline, we analyzed the Taranto gulf, detecting no evidence of significant deformation phenomena, despite the massive presence of anthropogenic industrial activities.
These case studies show how InSAR-based analysis enables observing ongoing deformation phenomena along the Italian coastline, and not only well-constrained, subsiding, and uplifting signals. InSAR retrievals, i.e., deformation rate maps, and time series can be appropriately coupled with sea level rise measurements, marine observations, and atmospheric modeling to forecast the future levels of sea rise, coastal erosion, and flooding, thus estimating the effects of the ongoing climate changes on coastal areas. To these purposes, the scientific community is moving toward designing an InSAR-based semi-automatic procedure for performing a near real-time coastline surveillance service.
It is worth noting that the test cases shown here were studied by X and C-band SAR data provided by the CSK (still active) and ENVISAT (inactive from 2010) missions. The main problem of such missions is the satellite revisit time of about 16/35 days, which prevents a near real-time updating of the results. In addition, they are characterized by different geometries of view (ENVISAT incidence angle from 15° to 45°; CSK incidence angle from 20° to 59°), area coverage (ENVISAT ~ 100 × 100 km; CSK ~ 40 × 40 km), and spatial resolution (ENVISAT ~ 20 × 20 m; CSK ~ 3 × 3 m in StripMap mode), which are indeed elements to account for.
In theory, such differences can be exploited as a sort of data cross-validation and improve the knowledge of a phenomenon according to its scale and orientation. However, they can create ambiguity and mislead non-expert users such as employees within the Civil Protection Department, Ministry of Environment, and so on. Therefore, in view of performing a near real-time coastline surveillance service, such differences in data and results have to be avoided, adopting the same protocol in all cases (including results updating, processing algorithms, and representation results). Moreover, other issues should be taken into account when considering SAR data: (a) the one-dimensionality of InSAR-based deformation measurements, i.e., along with the LOS direction. It does not enable providing a three-dimensional displacement field, especially when SAR acquisitions are available along one single orbit; (b) the limited spatial coverage (CSK ~ 40 km; ENVISAT ~ 100 km); (c) the processing resources, both in terms of hardware and software facilities, which require a high-performing workstation for the fast processing of long SAR image stacks, as well as resources for the storage of a great amount of SAR products (data and intermediate results).
Some of these constraints can be easily overcome through the joint use of S1 acquisitions and the cloud-computing resources. The schedule of S1 acquisitions goes in the direction of reducing the number of variables, thus ensuring the comparability and homogeneity of the results. Also, the S1 mission provides acquisitions with a spatial coverage of about 250 km in the default interferometric wide (IW) Swath TOPSAR acquisition mode (along both the range and azimuth directions) and a revisit time of six days, thanks to its two satellite constellations. This means ensuring the coverage of the whole Italian territory with only four S1 images, with less than one-week temporal sampling (updating) of SAR deformation maps that are able to follow a high rate of displacement phenomena. Also, they are free and can be easily accessed by anyone. The main critical issue concerns the limited simultaneous download from the ESA data hub (maximum two images at a time) and the storage of these data (each S1 image has a size of ~7 Gb each) and related derived products, which is becoming a predominant issue when dealing with the new ESA S1 missions. Then, it is strictly needed to have a high-power computing (HPC) service to automatically download the S1 data and upload semi-automatic processing tools to be run weekly or monthly according to the users’ requirements.
Once the infrastructure has been defined, the theoretical processing chain retrieving InSAR-based products, i.e., deformation rate maps and time series, could be as follows: (1) automatic S1 data download, every six days; (2) setting the maximum temporal baseline to be used for estimating the interferometric pairs and updating the results. For example, every new S1 SAR acquisition could be coupled with the six previous ones (i.e., maximum temporal baseline equal to 36 days) for the results updating. Restriction on the maximum perpendicular baseline should not be needed, since the S1 orbital tube is very stable; (3) choosing how often to update and provide the results (weekly, monthly, annually, etc.) and set the processing tools accordingly; (4) providing open access results through web-based geographical information system (GIS) infrastructure for monitoring and scientific purposes.
5. Conclusions
The combination of measurements, sensors, platforms, and processing techniques is fundamental for better understanding the processes ongoing along coastal areas such as land subsidence, sea-level rise, and volcanic activities. The spatial and temporal changes of such phenomena must be measured accurately, and the contribution from several driving factors must be estimated and separated. In such a context, we analyzed by InSAR data some case studies along the Italian coastlines affected by either small or extreme deformation rates, and strictly correlated with natural and human-based risks, such as volcanic and hydrogeological risks.
The experimental results suggest that InSAR data, together with geomorphological assessments and in situ measurements, contribute to unveil the complexity of geohazards along the Italian coastlines.
Then, in view of developing an appropriate infrastructure for the storage and semi-automatic processing of S1 SAR data, we believe that A-DInSAR methods have the capacity to perform a near real-time surveillance service if applied to SAR data with a short revisit time, i.e., six days.
Some steps toward this have already been done [
35,
36], but performing a robust and reliable service should be the aim of the entire SAR scientific community in order to support the responsible organizations in designing optimal risk mitigation strategies in coastal areas characterized by significant hazard levels.