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 ﬁrstly 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 2 , 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


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 −1 . 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.

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.

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 −1 .
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.

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.

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.

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).
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).
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.
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. 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.

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 −1 (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).

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.
Author contributions. Conceptualization of the study, GB, LV; interferometric analysis, SDC, FN; GNSS analysis and SAR data calibration, LV; outlier detection and geostatistical analysis, FF, AL; discussion and interpretation, GB, FB, PS, LV; writing-original draft, GB, FF, AL, LV.