Updating the subsidence map of Emilia-Romagna region (Italy) by integration of SAR interferometry and GNSS time series: the 2011–2016 period

The analysis of the vertical movements of the soil in the Po River plane of the Emilia-Romagna Region (Italy) was updated through an interferometric analysis referred to the 2011–2016 time-span. This activity is a continuation of previous studies on the state of knowledge of vertical soil movements in the same area, analyzed firstly by levelling and GNSS and more recently by SAR interferometry for the periods 1992–2000, 2002–2006, 2006–2011, on behalf of the Emilia-Romagna Region. The survey area analysed was approximately 13 300 km², which corresponds to the territory of the regional plain. The interferometric dataset was calibrated through the use of velocity time series of several permanent GNSS stations. Among the 36 stations analysed, 22 were included in the study area: 16 were used for the calibration and 6 as check points). The velocities required for the calibration of the SAR analysis were calculated in the period following the important seismic events that struck the territory of the Emilia Romagna Region in May 2012. The interferometric analysis was carried out by TRE ALTAMIRA using the SqueeSAR™ technology. In particular, in order to update the interferometric dataset to 2016, it was necessary to perform a joint processing of the available RADARSAT-1 data and of the data acquired by the RADARSAT-2 satellite using a specific operating mode of the SqueeSAR™ algorithm known as stitching; this approach allowed the joint processing of images acquired in the same geometry by these two satellites. The study of the time series of the GNSS permanent stations used to provide the velocity datum to the interferometric analysis, is described, and the results of the SqueeSAR™ interferometric processing are reported. Statistical analyses on the spatial distribution and the type of scatterers have been performed during the screening and validation procedures of the dataset, and for the identification and removal of the outliers. Finally, the resulting map is described in order to analyse the measured soil movements with respect to the results obtained in past analyses, and the possible geological and human-induced causes, which could have produced them.

1 Introduction

This work focuses on the alluvial plain of the Po River valley, Emilia Romagna region (Italy). This sedimentary basin has been affected by a widespread land subsidence of both natural and anthropogenic origin, studied since several decades (Caputo et al., 1970; Arca and Beretta, 1985). Ongoing natural factors (i.e. tectonic, consolidation, oxidation, organic soils shrinkage) still contribute a few mm yr⁻¹. Human-induced subsidence reaches higher values and it is mainly due to the groundwater exploitation. It was mostly intense after World War II, mainly because of groundwater pumping and, subordinately, because of the gas production from a number of deep onshore and offshore gas reservoirs (Teatini et al., 2006). During the second half of the 20th century, several cm per year were observed in the Po River delta and in the Bologna area, where a significant decrease of the subsidence is observed in the last period, as pointed out by this analysis.

2 Materials and methods

Since the 1950s, different agencies handled the subsidence monitoring in the Emilia Romagna region by geodetic levelling surveys. The campaigns have been initially performed in the localized areas where the phenomenon had become particularly evident, without the benefit of a consistent geodetic network design at regional scale. Over the years, different approaches were adopted, comprising three main techniques: geodetic levelling, GNSS and SAR interferometry.

2.1 Subsidence monitoring in Emilia-Romagna region: previous activities

In 1999, ARPA-Emilia-Romagna (Regional Agency for Environmental Prevention in Emilia-Romagna Region), in collaboration with University of Bologna, designed and instituted a network of 2300 levelling benchmarks, connected to 60 GNSS stations; it was the first integrated regional-scale monitoring geodetic network (Bitelli et al., 2000). Both the levelling and the GNSS networks were surveyed in 1999; the GNSS survey was repeated in 2002.

During the successive campaigns, different techniques have been integrated. In 2005, the high precision geodetic levelling of a subnet of the 1999 network (more than 1000 benchmarks), and the radar interferometric analysis (PSInSAR™ – Permanent Scatterers SAR Interferometry, by TRE – Tele-Rilevamento Europa) have been integrated to update the subsidence measurements (Bissoli et al., 2010). PSInSAR™ analysis was conducted using European Space Agency's (ESA) ERS1 and ERS2 data for the interval 1992–2000; data from ESA's Envisat and Canadian Space Agency's Radarsat (RSAT) satellites have been processed for the 2002–2006 period. The SqueeSAR™ technology (developed by TRE ALTAMIRA) was further applied on RSAT ascending data for the 2006–2011 campaign. No geodetic levelling surveys were performed during this period. Seventeen GNSS-derived positioning time series were processed both for inserting the whole survey in an international geodetic datum and for the calibration and validation of the SqueeSAR™ results. After sub-sampling on a 100×100 m grid based on the radar coherence, the final dataset was composed of about 320 000 PS (Permanent Scatterers) and DS (Distributed Scatterers) points, more than double the number of points derived from the 2002–2006 dataset. Geostatistical interpolation methods were applied to compute a dense regular grid (100 m × 100 m) of ground vertical movements covering the Emilia-Romagna Po River plain (Bitelli et al., 2014) and to subsequently generate a contour-based thematic map by isokinetic isolines with a contour interval of 2.5 mm yr⁻¹.

The workflow process is described in Fig. 1; the same approach has been adopted, with some fine-tuning, in the processing of the 2011–2016 campaign.

Figure 1 The workflow from data acquisition to final map production.

2.2 Data collection for 2011–2016 period

The interferometric analysis for the 2011–2016 period, covering the territory of the Emilia-Romagna regional Po River plain (red polygon in Fig. 2), was carried out by TRE ALTAMIRA using the SqueeSAR™ technology. The radar dataset includes imagery from RADARSAT-2 (RSAT) and COSMO-SkyMed (CSK) in ascending geometry with a resolution of 20 m × 5 m and 3 m × 3 m respectively (Table 1); the delimitation of the areas for the six tracks is shown in Fig. 2.

Table 1 The 6 tracks of processed radar with their area (km²) and point density (dens.: points km⁻²).

Figure 2 The area covered by the 6 tracks of processed radar data.

2.3 GNSS analysis and SAR data calibration

Radar data calibration was performed using measurements from GNSS stations, following the basic approach of the previous campaign. A network of 16 permanent GNSS stations within the study area was used to perform an accurate SAR calibration. An additional 6 permanent GNSS stations were used as control. The movements of the 22 permanent reference stations (inside the blue area of Fig. 3b) have been processed within a reference European network, depicted in Fig. 3a.

Figure 3 (a) GNSS permanent stations used as reference to frame the Emilia-Romagna geodetic network within the International Terrestrial Reference Frame (ITRF14). (b) Locations of the GNSS permanent stations processed for data SAR calibration. In blue the study area.

2.4 Outlier detection and geostatistical analysis

Prior to the generation of the subsidence map, a thorough statistical analysis was performed in order to identify outliers in the acquired dataset. Any potential outlier could degrade the subsidence estimation. In this scenario, an outlier is an anomaly of a single radar target, or a small group of radar targets, with a low spatial correlation to surrounding radar targets. Therefore, their vertical movement is deemed to be unrelated to subsidence phenomena, and ascribed to problems related to interferometric processing or local deformation such as thermal deformation or the construction of new structures, buildings or roads. Outlier detection is both a crucial and delicate process, carried out through the application of statistical procedures and an objective approach. The process needs the evaluation and control of an experienced operator.

The dataset was preliminary filtered by the selection of a subset of points that present coherence stability in the time series. The density of those points is obviously higher in urban areas than in the vegetated ones. In our case, the selected coherence threshold was set to 0.7 after several analysis on study area and from previous experiences: any scatterer, permanent or distributed, with a coherence (COH.) lower than 0.7 was discarded from further processing (Table 2).

Table 2 Comprehensive list of radar targets before (Origin) and after (COH.≥0.7) the application of coherence threshold.

The total number of points preserved for further processing was 1 936 132 after having discarded 1.93 % of points. Considering each track, the largest percentage of points were discarded from the Parma track, which is the track with the highest density realized by CSK. No points were discarded from the Mirandola track (Tables 1 and 2).

Figure 4 Spatial distribution of 338 isolated outliers within the study area.

Figure 5 Spatial distribution of outliers in cluster within the area of the case study.

An automatic procedure was used to isolate the outliers from the filtered dataset. The procedure estimated a predicted displacement speed at each point and compared its value against the measured displacement speed. The analysis was performed with the use of Kriging algorithm: an exact interpolator for the optimal estimation of a measurement distributed over an area. All the values are estimated using the available observations and providing a reliability information for each known value. Then, the results obtained by the interpolation derived from the Kriging processing are used as input to a Cross Validation (CV) phase. CV is the procedure that computes the value of each point and relies on spatially adjacent data; the difference between the predicted value and the measured value is defined as a residual.

The procedure was performed iteratively in two subsequent steps, with the results shown in Table 3. At each step, the points with the standard deviation of the residual value greater than the defined 5σ threshold were classified as outliers. The remaining points, below the threshold, were kept for the next step.

Table 3 Outliers detection in 2 steps.

This analysis identified a total of 21 638 outliers, equal to 1.1 % of the input dataset. Two procedures were applied for further validation before discarding the outliers. A total of 338 outliers were identified as isolated, with less than 5 valid points within a radius of 500 m (Fig. 4). After a careful validation by a direct inspection, 24 of them were reintroduced in the output dataset.

A total of 1784 outliers were identified as grouped in clusters, potentially representing significant local phenomena such as new construction. After an automatic approach using density maps in a GIS analysis and further validation, about 33 % of them were reintroduced in the output dataset for a total of 1 915 122 points.

Figure 6 Updated subsidence map of Emilia-Romagna region with isokinetic isolines.

In a final step, 2341 points were removed from the dataset after a supervised analysis that identified local phenomena unrelated to subsidence. The final dataset is then composed by 1 912 781 points.

3 Results and discussion

The 1 912 781 points of the final dataset were processed with geostatistical interpolation methods (Kriging) to produce a dense regular grid of ground vertical movements with a resolution of 100 m × 100 m, which is a useful resolution to represent the phenomenon. The grid was then clipped and only cells inside the Emilia-Romagna plain were preserved, as in the previous steps a narrow buffer around the in-land boundary was considered in order to reduce border effects. A thematic map showing isokinetic isolines was finally generated with a contour interval of 2.5 mm yr⁻¹ (Fig. 6).

Compared to the precedent 2006–2011 time interval investigated (Bitelli et al., 2014), subsidence has decreased in a large portion of the Emilia-Romagna plain, and in some places a small rate of uplift was indicated (green areas in Fig. 6). Quite impressive is the change in some parts close to Bologna area, where a strong subsidence has been recorded until the recent past (Bitelli et al., 2014). Uplift seems to be related to a large and fast rising in piezometric level, due to a reduction in pumping from water wells (ARPAE, 2018). Uplift recorded in the western part of the study area was also registered in all the previous surveys and is related to tectonic movements (Cenni et al., 2012).

4 Conclusions

The land subsidence results for 2011–2016 period, described in this paper, update the long history of measurements of vertical movements for the Po River plain that began in the 1950s. Different evolving approaches and technologies were adopted for these studies over this period. The use of modern SAR interferometry techniques, combined with accurate verification and validation procedures, demonstrate the strength of this approach over very large regions. In particular, the integration with GNSS was crucial for the calibration.

The results derived from this study indicate continued subsidence in some areas of the Emilia-Romagna region, whereas a reduction in the subsidence rate and even uplift was indicated in other areas.