Evolution of the techniques for subsidence monitoring at regional scale: the case of Emilia-Romagna region (Italy)

Abstract. The recent decades have seen a significant evolution of the methodologies and techniques for the monitoring of subsidence on a regional scale: from the traditional levelling technique to GNSS and finally to SAR interferometry. The case study of Emilia-Romagna, Italy, is a prime example of this evolution. As known, the Emilia-Romagna plain is subject to a phenomenon of subsidence with a natural and an anthropogenic component, both of varying amounts depending on the area. The first contributes a few mm/year; the second, particularly evident in the last 60 years, is mainly correlated to excessive withdrawal of fluids from underground and reaches higher values (in the past, subsidence rates of several cm per year were observed in the Po delta and near Bologna). The geodetic monitoring of subsidence started in the 1950s by different entities, establishing and measuring levelling networks of varying size and with various characteristics, mainly located where the phenomenon was most clearly manifest. These local initiatives were not able to provide a consistent understanding of the phenomenon throughout the entire Emilia-Romagna plain. The first regional-scale monitoring of the Emilia-Romagna plain was initiated in 1999, with a large levelling network (about 3000 km) and a coupled network of 60 GNSS points. In subsequent years, the monitoring approach has mainly focused on the most modern remote sensing techniques integrated with each other, with the adoption of the method DInSAR calibrated to a GNSS Continuous Operating Reference Stations (CORS) database. The application of DInSAR methods resulted in subsidence maps with a greater level of detail. The paper analyzes the methodology choices made during 1999–2012, through three successive campaigns that adopted and integrated the different techniques.


Introduction
Monitoring of subsidence at regional level was done in the last decades with different approaches and methodologies, following the development of the technology in the land surveying disciplines.
From the scientific literature regarding subsidence measurement, and more specifically from the papers presented at the latest editions of the International Symposia on Land Subsidence, the evolution is evident from the three main techniques adopted in this process: geodetic levelling, GNSS and SAR interferometry.
Whereas geodetical infrastructures are available, the integration of the three techniques can constitute the best approach to the problem, when datum definition and transformations permit all the datasets can be aligned into a single reference frame.Large levelling networks are still very widespread and utilize a well-established method, especially useful when historical data for a large number of benchmarks are available and can be updated; the accuracy level is lower than 1 mm per km.GNSS networks (Hofmann-Wellenhof et al., 2008), mainly implemented today by permanent stations, can provide point information with a lower spatial density and good accuracy, depending on different surveying strategies and the time span of the acquisition.For example, using current technologies of measurements and of data analysis, GNSS permanent stations are able to provide measurements of velocity components with an accuracy of about 1 mm year −1 in the vertical component, after at least 3 years of continuous acquisition (Kontny and Bogusz, 2012).
SAR interferometry, using the Persistent Scatterers approach (Ferretti et al., 2001;Hanssen, 2001;Ketelaar, 2009), can provide high density point information about the rate of the phenomenon with precision on the order of 1-2 mm per year.
The Emilia-Romagna region, Italy, is a prime example of this evolution.The Emilia-Romagna plain is subject to the phenomenon of subsidence with a natural and an anthropogenic component, both of varying amounts depending on the areas.The first contributes a few mm/year; the second, particularly evident in the last 60 years, is mainly correlated to excessive withdrawal of fluids from underground and reaches higher values.
Since the 1950s different agencies have designed and managed subsidence monitoring networks, measured by geodetic levelling, in the areas in which the phenomenon had become particularly evident (Caputo et al., 1970;Pieri and Russo, 1984).These local monitoring networks highlighted maximum subsidence rates of more than 25 cm year −1 during 1951-1962 in the Po delta, more than 11 cm year −1 during 1974-1981 near Bologna and more than 4 cm year −1 during 1972-1977 in the Ravenna area (Bissoli et al., 2010).More recently, rates of subsidence have significantly decreased but still remain important in some places.
This paper describes the experience carried out using geodetic techniques to monitor the phenomenon and the results of the three regional campaigns completed since 1999 by ARPA-Emilia-Romagna (Regional Agency for Environmental Prevention in Emilia-Romagna Region) on behalf of the Emilia-Romagna Region and in collaboration with the DICAM Department of Bologna University.

The first campaign (1999 and 2002)
The first campaign designed and instituted the geodetic infrastructure composed by two networks (Fig. 1): a high precision levelling network composed of over 2300 benchmarks, and a GNSS network consisting of about 60 points to be surveyed by long static sessions (Benedetti et al., 2000;Bitelli et al., 2000).
The levelling network was conceived with the aim to connect, where possible, a number of existing networks, dealing with a very heterogeneous situation: each authority or agency used different specifications and, therefore, different precision for their networks, different observation periods and different reference benchmarks.The 1999 levelling survey, accomplished in 75 days, was one of the first uses of large networks entirely surveyed by digital levels; the a-posteriori mean square error was less than 1 mm per km.
The scope of GNSS surveys was two-fold: to improve the gravimetric geoid model in the area and to establish quick and inexpensive control of some points in order to decide, for example, the measurement repetition frequencies of some levelling lines.
Both the levelling and the GNSS networks were measured for the first time in 1999; the measurement of the GNSS network was repeated in 2002 (Bonsignore et al., 2003).
A first attempt to evaluate the subsidence from this first survey was done by comparing some benchmarks heights with those from previous surveys, after some homogenization procedures in the vertical datum.The results largely showed the main trends of the phenomenon.
A specialized relational database was developed as an archiving and querying system, coupled with a WebGIS implementation, to improve the management of the huge amount of data coming from historical and 1999 and 2002 records.

The second campaign (2005)
In 2005, the knowledge of the subsidence phenomenon was updated through the integration of two techniques: high precision geodetic levelling of a subnet of the regional network (more than 1000 benchmarks, about 50 % of the 1999 network; Fig. 1), and radar interferometric analysis, adopting the PSInSAR ™ technique (Ferretti et al., 2001) in the standard PS analysis low resolution method.The study covered the whole plain territory of the region, about 11 000 km 2 .The interferometric analysis focused, in particular, on two distinct periods: the first time interval 1992-2000 refers to the processing of data from the European Space Agency's ERS1 and ERS2 satellites, while the second one, for 2002-2006, refers Proc.IAHS, 372, 315-321, 2015 proc-iahs.net/372/315/2015/to the processing of data from Envisat (ESA) and Radarsat (Canadian Space Agency) satellites.
The levelling results and the analysis of the discrepancies between forward and backward runs of the benchmark height differences showed the good quality of the campaign (Fig. 2), with the indication that it was possible to reduce the threshold level on the discrepancy to the value of 2 × √ L. A critical aspect for the comparison of results from the levelling campaign and from the interferometry was the choice of a reference point to be considered fixed in subsequent measurements or processing: the origin was fixed at a benchmark situated near Sasso Marconi (Bologna), in the Apennines foothills, located in a position almost barycentric (in a longitudinal sense) with respect to the network and considered stable in previous large surveys.
During the analysis of the results (Bissoli et al., 2010), some tests were performed -within a GIS -to have an objective evaluation of the intrinsic reliability of the radar results and of the comparison with results of the levelling campaigns for similar periods (levelling 1999-2005 vs. PSInSAR ™  2002-2006).The interferometric technique provided measurements of vertical relative motion characterized by a precision of the same order of those obtainable using levelling.The distribution of measurement points, in the case of urban environments, however was considerably higher than the number of benchmarks usually measured per square kilometer using levelling techniques.Standard deviation and coherence values derived from the interferometric processing for the PS points were determined for the three types of radar sensors used.As expected, the estimated velocities for the period 2002-2006 showed higher standard deviation values and lower coherence than during 1992-2000, because of the reduced number of SAR images available.
The results associated with PS points were converted to a grid model by geostatistical interpolation, in order to provide a more expressive representation of the subsidence velocity for the entire Emilia-Romagna plain (Fig. 3).

The third campaign (2011)
The third and last campaign, which will be described in more detail, was accomplished by integrating SAR interferometry and GNSS data collected at permanent stations (Bitelli et al., 2014).SAR analysis was performed by TRE -Tele-Rilevamento Europa through the SqueeSAR ™ technology, patented algorithm PSInSAR ™ second generation (Ferretti et al., 2011), using a RADARSAT-1 ascending dataset, acquired between 2006 and 2011.The analysis led to the identification of more than 2 000 000 measurement points, with an average point density of 200 MP km −2 .The full resolution data was sub-sampled on a 100 × 100 m grid using the radar coherence as selection parameter in order to obtain around 320 000 PS (Permanent Scatterers) and DS (Distributed Scatterers) points, more than double the previous campaign.Based on previous knowledge of the area, supported by GNSS measurements, it was assumed that horizontal motion was negligible, and the velocity of movement derived by SAR had been reprojected from the Line Of Sight to the vertical.

SqueeSAR TM dataset assessment and outliers rejection
In order to filter the outliers present in the SqueeSAR dataset which degrade the results, some procedures were implemented to identify points characterized by "abnormal" velocities compared to their surroundings, indicative of different phenomena not related to the subsidence of a regional nature, object of this study.
These "abnormal" points may arise due to extremely localized movements related to sagging of individual structures, settling of recent constructions, or in some cases, particularly in agricultural areas, they could be attributable to changes in soil moisture that involve phase shifts in the SAR signal, erroneously identified as movements.
Identifying "abnormal" points is an important process, which must be conducted with objective methods based on statistics.Given the large amount of data to be processed, the main part of the analysis must be automated, but ultimate identification of outliers should be supported by specific assessments guided by a visual inspection.

Definition and calibration of velocity datum
In planning this campaign, a first attempt was to overcome the previous situation for the definition and calibration of the velocity datum.The study of subsidence, in a similar way as any other geodetic phenomenon that is described with respect to a geographically defined reference frame, requires the definition of a geometrical and velocity datum to allow the comparison between observations collected at different times (epochs).This datum can be defined according to geometrical criteria or external measurements (i.e.GNSS time proc-iahs.net/372/315/2015/Proc.IAHS, 372, 315-321, 2015 series), or according to more qualitative criteria related to geological setting.
For the levelling networks, the vertical datum was defined, for example, through the a priori knowledge of the elevation of at least one reference point, taking into account its vertical velocity if it is not 0 mm year −1 .As previously mentioned, a levelling network was designed to monitor the regional subsidence, which was referenced to a benchmark located in the Apennines foothills close to Bologna.The stability of this benchamrk was verified through the comparison of vertical velocity obtained through the geodetic observation collected at the co-located GNSS-VLBI Observatory of Medicina (Bo), with respect to the estimation of the vertical velocity obtained from the repeated geodetic levelling between the reference benchmark and the GNSS-VLBI Observatory.
In the third campaign, no geodetic levelling measurement were done, and the datum was defined, as much as possible, in an "absolute" way for the PSInSAR ™ analysis through a network of 17 GNSS permanent stations (Fig. located in the study area and belonging to three different geodetic networks: -three stations belonging to global networks IGS-EUREF (Global and European Service); -four stations belonging to the "Rete Nazionale Integrata Gps" (RING) of the INGV Institute, materialized mainly for studies of Geodynamics; -ten stations belonging to the network of permanent GNSS stations established by Foger -Foundation of professional land surveyors of Emilia-Romagna, for technical purposes and, in particular, for real time positioning.
GNSS-derived positioning time series were processed from acquisitions made every 30 s during 2007-2011.The daily solutions of each station were computed with respect to ITRF2008.For each component (North, East, Up) of the time series, an average linear velocity (mm year −1 ) was estimated.
Six stations served initially to establish the datum of the interferometric analysis, and the other 11 were used for validation through the comparison between the rate of vertical movement provided by GNSS and rates derived from the interferometric measurement points (MP) distributed close to each GNSS station.These points were individuated by a GIS procedure, because the locations were not generally coincident.
The selection of the MP for the comparison was made for a small area analyzed at full resolution using the SqueeSAR ™ method in the neighbourhood of the GNSS stations (Fig. 5).The velocities of the selected points (subjected to a preliminary visual screening) were then averaged, giving each point a weight based on its coherence.The resulting values of the differences with GNSS measurements for the six calibration stations were then calculated.Since InSAR results can be affected by low frequency artefacts due to uncompensated atmospheric component or orbital inaccuracy, a first order surface was estimated and removed from the differences between InSAR and GNSS velocities.The cal-  ibrated vertical velocity was then validated with the remaining 11 GNSS stations, adopting the same procedure for the identification of GNSS and MP match.The validation highlighted the very good agreement (Table 1) between the calibrated Squee-SAR ™ data and the GNSS network data (average residual = −0.36 ± 1.39 mm year −1 ).Some exception could be related to different causes: GNSS problems, low number or radar quality of the selected MP, relative position between MP and GNSS antenna, estimation inaccuracy, pres-ence of horizontal motion component, spatial and/or temporal motion variations in the monitored period.Considering all the factors, an overall uncertainty of the whole analysis is expected, as a precaution, to be within the value of about  From the final data set of 315 371 points, geostatistical interpolation methods were applied to compute a dense regular grid (100 m × 100 m) of ground vertical movements for the Emilia-Romagna plain for 2006-2011 and to subsequently generate a contour-based thematic map by isokinetic isolines with a contour interval of 2.5 mm year −1 (Fig. 6).The results showed in different areas a reduction of the subsidence phenomenon.

Conclusions
A short description of the subsidence monitoring for the plain area of the Emilia-Romagna Region (Italy) was presented, describing the methodologies used in the 1999-2012 period.
During three successive studies, different approaches were adopted integrating the main available geodetic techniques.Some aspects of the whole experience are described, highlighting the main findings in terms of better spatial density, accuracy, definition and calibration of a velocity datum.The next survey is scheduled in 2016, following primarily the method adopted in the 2011 campaign.

Figure 1 .
Figure 1.Levelling network and GNSS points for the 1999 campaign.

Figure 2 .
Figure 2. Levelling campaign 2005: discrepancies vs. distance for the forward and backward runs of the section height differences.

Figure 4 .
Figure 4. Positions of the GNSS permanent stations used.The symbols with a circle indicate the stations used as calibration to establish the datum of the interferometric analysis.

Figure 5 .
Figure 5. Distribution of measurement points identified by the analysis SqueeSAR ™ at full resolution in the neighborhood of two GNSS stations; selected points are in red.The positions of GNSS antennas are marked with red triangles.

Table 1 .
Comparison of the residual in the velocities (mm year −1 ) between the values obtained from the GNSS permanent stations and those derived from the analysis of the closest radar measurement points.