1. Introduction
Land surface dynamics can be the evidence of different processes (i.e., swelling/shrinkage of expansive soils, aquifer-system compaction, tectonic movements, and consolidation due to the load of new buildings) of natural and anthropic origins. In many cases, the movements are due to the complex interactions of multi-driving forces, which act at various spatial and temporal scales [
1,
2,
3]. The interference and/or overlap makes the effective decomposition of the components of motion difficult. For this reason, only a few authors have performed the decomposition of multi-driving force processes that trigger surface movements [
4,
5].
The areas affected by displacements need effective identification methods, semi-continuous spatio-temporal surface deformation monitoring and mechanism recognition analysis, in order to adopt the appropriate risk mitigation strategies, and to develop a sustainable management of natural resources.
Advanced Differential Interferometric SAR (A-DInSAR) is a powerful remote sensing tool, capable of mapping displacements over wide areas at very high spatial resolutions. The technique is based on the processing of multiple interferograms derived from a large set of SAR images in order to obtain a displacement time series, along the line of sight (LOS) of the satellite, of radar _targets on the Earth’s surface [
6,
7,
8]. To date, several studies have documented the application of A-DInSAR data to monitor and to improve knowledge regarding different types of land surface dynamics [
9,
10,
11].
Up to now, geohazard-mapping methodologies are mainly based on the visual inspection of average LOS velocity, such in the case of the ground stability areas in the framework of the PanGeo project [
12], or in hotspot analysis and hotspot and cluster analysis of average LOS velocity to detect slow-moving landslides [
13,
14].
Most of the studies aimed at identifying the areas affected by ground motion, at a regional and local scale, were focused on the average rates of displacement detected by the A-DInSAR technique [
15,
16,
17]. Meisina et al. (2008) [
16] implemented a systematic methodology to identify anomalous areas selected on the basis of the average LOS velocity, and Peduto et al. (2015) [
17] introduced the average LOS velocity by weighting the radar _target values on their coherence value. Both methodologies were mainly based on the average LOS velocity analysis, constrained by the fact that previous A-DInSAR data allowed the obtaining of noisy time series that are difficult to interpret, showing some limitations in the distinguishing of ground deformation due to different processes, or in detecting shallow deformations caused by seasonal processes.
Taking into account the recent improvements of A-DInSAR data acquired by the COSMO-SkyMed satellites and the current ESA Sentinel missions, which act at higher spatial and temporal resolution, it is necessary to develop an appropriate methodology to analyze extremely large datasets that consist of a large number of measuring points.
The goal of this study is to present a novel methodology to analyze multi-sensors and multi-temporal A-DInSAR data for the geological interpretation of areas affected by land subsidence/uplift and seasonal movements. The approach is aimed at (i) improving the detection of ground motion areas; (ii) understanding the spatio-temporal evolution of ground motion areas; and (iii) supporting the interpretation of driving-force mechanisms. The proposed methodology is addressed to overcome limitations, such as the analysis of large data sets, by allowing exploitation of the great potential contained in A-DInSAR time series. In this paper, for the first time, ground motion areas are identified on the basis of their peculiar ground deformation behaviors.
The procedure was applied in the Oltrepo Pavese (Po Plain, Italy) by use of ERS-1/2 and Radarsat data processed using the SqueeSAR™ algorithm, and covering the time spans from 1992 to 2000 and from 2003 to 2010, respectively. The test area is a representative site of the Po Plain, characterized by various geohazards with moderate rates of motion. The obtained results are helpful for scientists and authorities in charge of land use planning in order to improve public safety and to assess the management of groundwater resources.
2. Study Area
The study area is located in the central area of the Po Plain, within the Oltrepo Pavese (Italy), covering an area of about 440 km
2 (
Figure 1a). The plain of the Oltrepo Pavese is constituted of alluvial quaternary deposits, originating from the combined action of Apennine streams that form coalescent fans and of the Po River. These deposits overly the Miocene-Pliocene marine substratum, constituted by sandy-marls, sandstones, conglomerates, gypsy-marls, and calcareous-marls [
18,
19].
Three main geomorphological units composed of quaternary deposits were previously distinguished: post-Würmian alluvial deposits, Würmian–Holocene alluvial deposits, and pre-Würmian alluvial deposits [
20,
21].
The post-Würmian alluvial deposits were sedimented by the most recent depositional activity of the Po River, and by the recent and present deposition from flooding events of the main Apennine streams. These deposits are localized in proximity to the Po River, and they are composed of sand, sandy silt, and silt. The deposits from the Apennine streams are characterized by gravel and sand. In general, post-Würmian alluvial deposits contain a shallow phreatic aquifer with a water level near ground surface.
The Würmian–Holocene alluvial unit is the most prevalent geomorphologic unit of the study area. These deposits are made up of alternating sand and gravel, with interbedded clays or argillaceous silt. A shallow phreatic aquifer and deeper aquifers, of both phreatic and confined types, were distinguished in this unit. Deeply buried structures, consisting of a series of folds and fold-faults, shaped in the tertiary marine basement, have a direct control over aquifer geometry [
22]. These alluvial deposits show the constant presence of a tabular clayey-silty cover, which acts as a seal, and limits water infiltration [
21]. Groundwater exploitation of the Oltrepo Pavese is mainly derived from the aquifer contained in this unit.
Conversely, the pre-Würmian unit, located in the southern part of the study area, consists of older fluvial terraces, and is constituted of gravel and sand with a silty matrix, characterized by lower permeability [
19]. These deposits are strongly weathered and covered with clay soils. Locally, loess deposits overlap this unit. The depth of the water level in this unit reaches values up to 20 m in depth.
In the study area, quaternary sediment thickness decreases from west to east, with the minimum thickness and the outcrop of the marine substratum in correspondence with the Stradella thrust, evidence of the neotectonic activity that affects the area [
23]. The distensive regime of the neotectonic activity supported higher alluvial sediments accumulation in the northwestern sector of the plain with respect to the southeastern sector [
24]. The structural setting of the Oltrepo Pavese plain is conditioned by the presence of an important NE–SW tectonic discontinuity, known in the literature as the Vogherese Fault (
Figure 1b) [
25]. This tectonic lineament is a portion of the composite seismogenetic source called Rivanazzano-Stradella, belonging to the Northern Apennines outer thrust front [
26], which triggered an earthquake of Mw 3.94 ± 0.34 in 1971 [
27].
From the climatic point of view, the plain area of the Oltrepo Pavese is characterized by a humid continental climate, with average annual rainfall of around 700 mm. Two rainy seasons, with maxima in May and in October–November, have been detected by previous authors [
28]. In the last two decades, the region has experienced many drought periods, such as from March 1989 to August 1993, and from May 1998 to September 2000 [
28].
Here, the analyses of the drought periods, performed by Meisina et al. (2003) [
28] were extended from 1960 to 2010, using available measurements from 1960 and to match satellite observation time intervals. Therefore, climatic measurements acquired at the Voghera weather station were exploited, and the annual rainfall deficit was computed as the difference between the monthly water balance (difference between precipitation and potential evapotranspiration) and the average monthly water balance obtained for the period from 1960 to 2010. A significant drought period was identified from 2003 to 2007.
It is worth noting that the plain of the Oltrepo Pavese is characterized by discontinuous urban fabric, and 68% of the monitored area is covered by non-irrigated arable land [
29]. The study area is a representative site of similar geological contexts in the Po Plain, where geohazards, due to natural and anthropic factors, were previously recognized. Volume changes of clayey soils; shrinkage during drying periods, and swelling during wet ones, resulted in severe damage to overlying structures [
30].
Engineering Geological Units
Recent investigations have been performed to investigate the geotechnical proprieties of the quaternary deposits in the first 15 m of depth of the Oltrepo Pavese plain. On the basis of the analyses of penetrometer results (176 Dynamic Cone Penetration Test, 20 Standard Penetration Test and 266 Cone Penetration Test), calibrated against soil profile logs, six geotechnical classes of non-cohesive soils (from I1 to I6) and four geotechnical classes of cohesive soils (from C1 to C4), have been previously distinguished [
31].
Taking into account the geotechnical classification of these soils, eight engineering geological units were introduced in order to reduce the lithologic variability and to map the geotechnical behavior of representative and homogeneous geotechnical profiles (
Figure 1b,c). Unit 1 represents the post-Würmian alluvial deposit of the Po River, and is mainly constituted of non-cohesive soils (I3). Conversely, Units 2, 3, 4 and 7 are Würmian–Holocene alluvial deposits, characterized by cohesive soils. Unit 5 is composed of the alluvial fan of the Scuropasso River, and Units 6 and 8 of the pre-Würmian terraced deposits.
Figure 1b highlights the prevalent reported occurrence of damaged buildings [
31], located over Units 8 (34%) and 6 (29%), which correspond with pre-Würmian deposits, and over Units 5 (20%) and 4 (9%), which correspond, respectively, with the alluvial fan of the Scuropasso River and Würmian-Holocene deposits, which are characterized by a high content of clay deposits that are susceptible to swelling/shrinkage processes due to the moisture changes in the dry/wet periods [
31].
The majority of problems are experienced by single-story family residences, founded on conventional shallow strip, concrete footings, which generally extend to a depth of between 1 m and 2 m below ground level [
30]. In most of the cases, the remedial measures adopted were pile underpinnings, which reduces the value of the house by 20% [
30].
3. Advanced DInSAR Data
The Persistent Scatterer Interferometry (PSI) technique represents a class of the A-DInSAR techniques that uses radar _targets on the Earth’s surface to produce, starting from a set of images, the displacement time series along the line of sight of the satellite (LOS) of individual persistent scatterers (PS) [
6]. The technique allows the obtaining of data regarding land deformation over wide areas with millimeter precision [
6]. However, it has some limitations in rural areas due to low density of PS and phase ambiguities. A new algorithm, named “SqueeSAR™”, has been developed in order to overcome this problem and to extract the deformation of distributed scatterers (DS), increasing the density of measuring points; for this reason, it was exploited to process the SAR scenes used in this work [
32].
The input satellite data for this analysis consist of SAR images acquired by sensors operating in the C-band (wavelength, 5.6 cm; frequency, 5.3 GHz) onboard the ERS-1 and ERS-2 and Radarsat (RSAT) satellites. The scenes were acquired in the ascending mode, covering the time intervals from 9 July 1992 to 2 August 2000, and from 24 March 2003 to 5 May 2010. The dataset acquired in the descending mode covers the time intervals from 3 April 1992 to 7 January 2001, and from 28 April 2003 to 5 January 2009. The ERS-1/2 images were acquired with a nominal repeat cycle of 35 days, while the Radarsat images were acquired with a nominal repeat cycle of 24 days. These scenes were processed with the SqueeSAR™ technique by Tele-Rilevamento Europa (TRE). The algorithm allows the extraction of movement measurements, not only from traditional persistent scatterers (PS), such as anthropic structures or rocks, but also from distributed scatterers (DS), such as sparse vegetated areas [
32]. This permitted a high density of A-DInSAR data over non-urban areas. The algorithm improved the quality of the displacement time series; atmospheric effects can be better estimated and removed, and so the resulting time series are characterized by less noise. Available SAR images were processed with the SqueeSAR™ technique to identify acceleration and slowing of scatterers, in order to better observe the kinematic evolution of seasonal and non-linear processes.
The processing results have precision along the LOS, usually better than 1 mm/year, depending on the amount of available data, the local PS density (a key element in the estimation of spurious atmospheric phase components), and the distance from the reference point.
The geocoding accuracy of the PS locations is around +7 m in the eastern direction, and +2 m in the north–south direction, while the estimated accuracy of the elevation values was better than 1 m.
Additional details of the SAR datasets are summarized in
Table 1. It has been possible to detect more than 20,000 PS-DS over an area of around 430 km
2 (
Figure 2). It worth noting that, as explained in
Section 2, the area is mainly characterized by non-irrigated arable land, which was classified by Plank et al. (2009) [
33] as being not at all suitable for the A-DInSAR technique using the C-, X-, and L-band sensors. To provide a quantitative assessment of the PS-DS density for this area, given the land use, the Corine Land Cover (CLC), which was implemented by the Lombardia region using color and infrared orthophotos from IT2007 (made by Blom CGR; 50-cm pixels), and named “DUSAF 2007” [
29] was exploited. More precisely, such analysis was carried out by selecting, respectively, PS and DS in the non-irrigable arable land (CLC code 211). The results reveal that the ascending scenes show a PS density of 3.29 and 18.69 PS/km
2, respectively, for ERS-1/2 and Radarsat sensors. Conversely, the PS-DS densities for the same dataset are 20.74 and 39.42 PS-DS/km
2. In the case of the descending scenes, PS densities are 8.83 and 13.56 PS/km
2, respectively, for ERS-1/2 and Radarsat sensors, and the PS-DS densities for the same dataset are 37.63 and 29.10 PS-DS/km
2. By applying the SqueeSAR™ technique, an improvement of measuring-point density was obtained for the non-irrigable land. However, the PS-DS density for this land use coverage still shows low densities, which are comparable to values obtained using the PSI technique over the non-irrigable arable land in Britain [
34]. Higher PS-DS density was detected in the continuous and discontinuous urban fabric (CLC code 11), and the road and rail networks and associated land (CLC code 122), which reach an average of 334.27 ± 43.57 and 137.35 ± 31.96 PS-DS/km
2, respectively.
4. Methodology
Taking into account the quality improvement of the A-DInSAR time series, it is necessary to develop an appropriate methodology to exploit the great potential contained in the A-DInSAR data for the geological interpretation of ground-motion areas. A systematic methodology to analyze the temporal evolution of ground motion areas is proposed here (
Figure 3).
The procedure consists of three main phases. In the first phase, the vertical and E-W components of motion are disentangled, and a displacement time series (TS) accuracy assessment is performed. In the second phase, different statistical tests are applied in order to find the spatio-temporal pattern of the principal components of movement, and the kinematic model of the _targets. The result of this step is the identification of areas with significant movement, so-called “ground motion areas”. Ground motion areas correspond to a cluster of a minimum 3 of PS, with a maximum distance of 50 meters, characterized by the same trends (linear, non-linear, seasonal). Ground motion areas do not have a geological significance, but the systematic procedure represents a useful and fast approach to detect sectors where A-DInSAR analysis detects spatio-temporal ground deformation, and where the attention of scientists has to be focused. Finally, the third step consists of a data fusion of the A-DInSAR data and the geological data to determine the causes of ground motion processes. Multidisciplinary data, such as geotechnical, hydrogeological, and land use data, are fundamental to recognize ground motion mechanisms
Figure 3 illustrates the proposed methodological approach, which consists of the following steps:
Step 1. Decomposition of the vertical and horizontal components of motion (Phase I). A-DInSAR techniques allow the obtaining of the average LOS velocities of _targets. Smaller values of the satellite incidence angle (usually 23°–35°) make the measurements more susceptible to the vertical component of motion with respect to the horizontal; however, ignoring the horizontal components may implicate an over-/underestimation of the deformation, leading to misinterpretation of the results [
35]. By the use of the datasets acquired in the ascending and descending modes, the vertical and horizontal components of _target motion can be extracted using the following equations:
V
asc is the average LOS velocity acquired in the ascending mode, V
desc is the average LOS velocity acquired in the descending mode, V
VERT is the vertical component of motion, and V
EW is the E–W component of motion. The
easc,
hasc and
edesc,
hdesc are the LOS directional cosines, respectively, for ascending and descending orbits. First, the average LOS velocity of each dataset was interpolated by an inverse distance weighted (IDW) approach, and, then, the component of motion was decomposed using these interpolated maps. The method was previously used by several authors [
36,
37]. If the horizontal component of motion is absent, the vertical component of motion can be extracted using the incidence angle (θ) of the dataset acquired in one mode (ascending or descending) using the following equation:
Step 2. A-DInSAR displacement time series accuracy assessment (Phase I). The displacement time series detected by each _target may be affected by local error or regional trends, which can be observed in the whole dataset. Therefore, post-processing analyses of the A-DInSAR data are fundamental in order to avoid misinterpretation of unreal ground movements. The check of the displacement time series is performed by selecting the most coherent (>0.9) _targets with an average LOS velocity in the range ±0.5 mm/yr., where no significant movements are expected. More precisely, the approach proposed by Notti et al. (2015) [
38] was applied in order to correct problems due to phase unwrapping, regional unreal trends, and anomalous displacement detected on certain dates (i.e., unreal movements at the same time as meteorological events, such as snowfall).
Step 3. Principal component analysis (PCA) of the A-DInSAR displacement time series (Phase II). A-DInSAR measurements represent the cumulative ground movements (natural versus anthropic, superficial versus deep displacements, multi-year, and seasonal components); the overlapping of several ground motion components may make the accurate interpretation of ground deformation processes difficult. Recent studies have applied PCA to satellite-based time series analysis [
39,
40] to detect spatio-temporal deformation patterns.
Here, principal component analysis (PCA), implementing a matrix organization of location versus time (T-mode), was performed in order to decompose the long-term (multi-year) and the seasonal components of ground motion. A matrix, in which each column contains the LOS displacements for each SAR image, and where each row contains the displacement time series of the _targets (PS-DS), was formed. The main outcomes are the correlation and covariance matrices, the eigenvalues and eigenvectors, the percent variance that each eigenvalue captures, and the principal component (PC) score maps. In interpreting the principal components, PC scores related to each observation (PS-DS) are useful for knowing the correlations with the principal components in the whole dataset; the higher values correspond to higher correlations with the analyzed PC. To detect the displacement time series of the PC, the displacement time series can be multiplied for the PC eigenvectors of each SAR image date.
Step 4. Identification of the kinematic model of the _targets (Phase II). A-DInSAR data proves to be efficient for the estimation of linear and non-linear ground motion movements [
8]. Hence, automated classification of displacement time series (TS) in large datasets are needed for TS interpretation at a regional scale. Here, an automatic sequential procedure, based on statistical tests, the PS-time program implemented by Berti et al. (2013) [
41], was used to classify the TS into one of three predefined _target trends: uncorrelated, linear, and non-linear. The program permits the detection of the date (break) where abrupt changes in slope in the non-linear TS are recorded. This tool also assigns an additional parameter to the non-linear TS, index “BICW”, related to the bi-linearity of the TS (i.e., TS that records BICW higher than 1.2 indicates a strong bi-linearity [
38]).
Step 5. Identification of “ground motion areas” (Phase III). This task aims at detecting ground motion areas in order to focus the interpretation on significant sectors. Physical processes can be characterized as linear, non-linear, and seasonal ground movements. While the average velocity of _targets with a linear trend is useful to separate movement and non-movement areas, the same parameter is not efficient for recognizing ground motion areas affected by non-linear and seasonal movements. Even _targets with an average LOS velocity in the stable range (±1.5 mm/yr) can be affected by significant seasonal and/or non-linear movements. Taking into account that non-linear and seasonal movements may cause considerable damage to buildings and infrastructure [
42], PCA results were exploited to identify ground motion areas with linear, non-linear, and seasonal trends. The ground motion areas were identified using a spatial cluster of significant PC scores. First, we define the threshold of the PC scores by applying a statistical inspection of the PC score frequency distributions. If the histogram of the PC scores shows a normal distribution (skewness = 0), we consider the threshold equal to twice the standard deviation, while, if the PC scores shows asymmetrical distribution (skewness ≠ 0), the threshold is defined as being equal to the interquartile range (IQR). Then, PS-DS with PC scores higher than the threshold are selected, and a buffer zone of 50 m around the detected measuring points can be derived.
Therefore, ground motion areas consist of clusters of a minimum of 3 PS-DS with a maximum distance of 50 meters. The approach to extracting ground motion areas is based on the procedure proposed by Meisina et al. (2008) [
16], but, herein, the delineated areas are derived from the PCA analysis (
Figure 4). The ground motion areas are not found on the basis of unstable velocity; hence, even seasonal and non-linear processes can be identified and taken into account for the interpretation of displacement time series. Following this, the kinematic model of the _targets is combined with the ground motion areas in order to distinguish linear and non-linear processes. More precisely, the TS with BICW higher than 1.2 are selected to distinguish non-linear TS from the linear ground motion areas.
Step 6. Interpretation of “ground motion areas” (Phase III). Once the ground motion areas have been identified, the processes are interpreted by integrating them into a Geographical Information System (GIS) with ancillary data (such as geological, geotechnical hydrogeological, compressible thickness map, digital orthophoto, and damaged buildings distribution). Mechanism recognition is based on a cross-comparison of the representative subsoil geological profiles, and the relative displacement time series with multidisciplinary information. Analysis of the breaks and the detection of the deceleration and acceleration periods are crucial to correlate the processes with the triggering mechanisms. The final outcome is the attribution of driving-force systems and triggering mechanisms to ground motion areas.
5. Results
5.1. Post-Processing Elaborations (Phase I)
The decomposition of the horizontal and vertical components of motion highlights that the motion is mainly vertical, both in the period of 1992–2000 and 2003–2010 in the Oltrepo Pavese. The results are consistent with literature data over the study area [
43]. Therefore, the east–west component of motion was considered to be negligible. Following this, the procedure was implemented on the datasets that shows the higher PS-DS density and a higher number of scenes (
Table 1) in order to analyze the datasets with higher spatial and temporal resolutions. Hence, the ERS-1/2 descending and Radarsat ascending data were exploited to monitor the time intervals, from 1992 to 2000, and from 2003 to 2010.
A-DInSAR displacement time series accuracy assessment has provided the detection of anomalous displacements recorded on 9 March 1997 and 16 July 2000 by the ERS-1/2 descending datasets, and on 10 December 2008, for the Radarsat ascending dataset. The anomalous displacement identified in the Radarsat dataset might be related to the snowfall occurred on the day of SAR acquisition. Consequently, these scenes were not included in the following analyses.
5.2. Identification of the Ground Motion Areas (Phase II)
Once the post-processing elaborations have been performed on the A-DInSAR data, the improved displacement time series have been exploited in order to delineate ground motion areas.
In order to achieve the spatio-temporal extraction of the principal components of motion, a matrix of the ERS-1/2 desc. and Radarsat asc. dataset was formed.
The results show that the PC1 explains 79% of the variance, and PC2 and PC3 explain, respectively, 4.7% and 1.8% of the variance in the ERS-1/2 dataset. In the Radarsat, 69% of the variance is explained by PC1, while 5.7% and 1.7 % are explained by PC2 and PC3.
Figure 5a,b shows the spatial distribution of the principal components score units. The spatial pattern of the PC of the deformation highlights that the first and the third principal components of motion mainly affect the southwestern sector of the study area, while the second component is mainly localized in the northeastern zone. It is worth noting that the distribution of the components of motion is heterogeneous across the Oltrepo Pavese, and that the comparison between the ERS-1/2 and Radarsat data reveals that the spatial distribution of the components of motion detected during the period 1992–2000 is correlated with that of the period 2003–2010 (
Figure 5a,b).
A visual inspection of the principal components eigenvectors (
Figure 5c,d) highlights that in the ERS-1/2 and Radarsat datasets, the first component (PC1) corresponds to the long-term lowering of the Earth’s surface. The second corresponds to the long-term uplift ground motion, and the third PC highlights seasonal deformations.
Ground motion areas were delineated using the PC scores, applying the procedure described in
Section 4, in order to focus the interpretation of the processes on significant sectors. The results reveal that 30%–40% of the ground motion areas, delineated using the PC1 of the ERS-1/2 and Radarsat measurements, is located in proximity of the town of Voghera and along the railway of Voghera–Pavia (
Figure 6), and it records an average LOS velocity in the range from −2 to −3.7 mm/yr, in the period 1992–2000, and in the range of −2 to −4.8 mm/yr in the period 2003–2010. A distinct concentration of ground motion detected by PC1 is visible only in the time span of 1992–2000, along Emilia Road, from the town of Broni to the town of Redavalle, and along the railway from Broni to the town of Arena Po. The average LOS velocity ranges from −6.4 to 2 mm/yr.
Ground motion areas mapped via PC2, of both ERS-1/2 and Radarsat data, are centered in the area between the town of Stradella and the town of Monteacuto (
Figure 6). The average LOS velocity records an uplift maximum of 3.40 mm/yr, in the period 1992–2000, and 2.58 mm/yr, in the period 2003–2010.
Ground motion areas observed using PC3 are mainly located in the southern sector of Voghera and in proximity to Lungavilla, Codevilla, Casei Gerola, and Broni (
Figure 6). In the second period, ground motion areas via PC3 are not present in the southern sector of Voghera, and the movements are mainly localized along the Coppa River (see location in
Figure 1) and in proximity to Verrua Po. Seventy-four percent and 97% of the ground motion areas detected using PC3, respectively, for the ERS-1/2 and Radarsat data show an average LOS velocity in the range of ±1.5 mm/yr.
As was previously explained, an automated time series classification tool was used to distinguish uncorrelated, linear, and non-linear trends. The outcome of this step reveals that around 40% of the _targets are classified as uncorrelated, both in the periods 1992–2000 and 2003–2010. Thirty-five percent and 22% of the _targets are classified, respectively, as non-linear and linear during the time interval 1992–2000. In the following period (2003–2010), 28% and 27% of the _targets are classified as linear and non-linear, respectively.
Furthermore, the results of the kinematic model analysis were combined with the ground motion areas detected using the PCA approach.
Figure 6 shows the results of the overlapping of the ground motion areas, detected using PC1, PC2, and PC3, and the non-linear PS-DS characterized by a strong non-linearity (BICW > 1.2). Mainly, changes in trends in displacement time series were identified at the end of 1999 and 2008 in the areas of Broni, Voghera, Lungavilla, and Casei Gerola, and, in most of the cases, the non-linear _targets are superimposed on seasonal ground motion areas detected via PC3, which represent the seasonal components of motion.
Overall, the ground motion area detected at Voghera records a linear displacement time series that is active during both periods, while, in Broni, the ground motion occurs only in the period from 1992 to 2000.
The ground motion area centered in the north sectors of Voghera could be related to the presence of over 6–7 m of shallow thick clay deposits in the subsoil, and to groundwater exploitation for industrial supply. Unfortunately, the lack of historical data regarding the piezometric level variations did not allow us to analyze the effects of its variations on the land subsidence recognized in this area. The phenomenon is characterized by a linear trend for both 1992–2000 and 2003–2010, and no evident acceleration has been recognized (
Figure 6).
The ground motion areas along the railway Voghera–Pavia and along the railway from Broni to Arena Po show a linear trend and rates of average LOS velocities in the range of −5 to 1 mm/yr for both periods. Given the close association with the local transport network, these ground motion areas could be due to the consolidation of the railroad embankment, and anthropogenic activity could have an influence on the surface settlement of these areas.
Areas of moderate uplift were observed via PC2 of the ERS-1/2 and Radarsat data, in the sector from Stradella to Monteacuto, where geomorphic evidence of an active emergent thrust was previously observed [
43]. The terraced surfaces derived from topographic profiles, combined with inferences on the ages of the terraces of the exanimated area, were used to carry out estimations about the late quaternary uplift rates across the escarpment. The agreement between the distribution and trend of deformation with the uplift trend, found by previous authors in the same area, suggests that the deformation could be due to deep-ground motion related to tectonic movements.
5.3. Mechanism Recognition (Phase III)
A spatio-temporal analysis for mechanism recognition was performed, considering some factors and evidence, which might have relevance in explaining the patterns of ground motion, as observed in others ground motion areas by previous authors [
44,
45,
46]: (i) the geotechnical properties of the deposits; (ii) the thickness of the superficial clayey soils and the hydrogeological setting; (iii) the distribution of damaged buildings; and (iv) land use change effects.
Different natural and man-induced processes were recognized such as swelling/shrinkage of clayey soils located over engineering geological Unit 8, and seasonal ground motions due to seasonal groundwater level variations, which are located over engineering geological Units 3 and 4.
Land subsidence due to the load of new buildings was observed over engineering geological Units 1, 2 and 3, and moderate tectonic uplift was recognized in the sector from Stradella to Monteacuto, where geomorphic evidence of an active emergent thrust was previously observed [
43].
The main conditioning and triggering factors on ground motion are examined in detail in
Section 6.
6. Discussion
The approach of this work was to identify areas with significant movements where the attention of scientists has to be focused by the use of A-DInSAR time series analysis. Ground motion areas were mapped on the basis of their peculiar ground deformation behaviors using a systematic and reproducible procedure.
Although various methodologies were implemented to map the ground motion areas using the average velocity detected by A-DInSAR data [
12,
13,
14,
15,
16,
17], ground motion areas mapping on the basis of A-DInSAR time series is still not a common practice in the scientific community.
While the average velocity may be useful to investigate physical processes characterized by linear trends, the same parameter is not efficient to detect and interpret ground motion areas affected by non-linear and seasonal movements. In these cases, the use of the average velocity may result in misleading interpretations of the process.
Our intent is to propose a novel approach to take into account different deformational behaviors. Hence, the reliability of the results obviously depends on the quality of the dataset. However, the first phase of the procedure is focused on the time series check in order to correct problems due to phase unwrapping, regional unreal trends, and anomalous displacement detected on certain dates. In addition, if the whole dataset is affected by a relevant error, it can be easily detected by applying the PCA, being clearly visible as principal component of motion, and could be subtracted from the entire dataset.
This approach allowed us to detect ground motion areas due to different processes, even at low average LOS velocities in the range of ±1.5 mm/yr, which are affected by geohazards that may impact public safety.
In this section, the main conditioning and triggering factors on ground motion were discussed.
6.1. Comparison between Ground Motion and Engineering Geological Units
As regards geotechnical interpretation, the relationship between the engineering geological units and the ground motion was analyzed by computing the maximum and minimum of the displacement measured for each unit (
Figure 7). Furthermore, the PS-DS density for each unit was extracted in order to take into account the distribution of measures.
Figure 7a,b shows the results of the spatial analysis for the periods 1992–2000 and 2003–2010. It is worth noting that the extreme values of cumulated displacement observed in the first period are higher than those of the second one. A decrease of deformation was observed, and, in general, the cumulated displacement decreases from Unit 1 to 8. The PS-DS density reaches an average for the engineering geological units of 108 and 113 PS-DS/km
2, respectively, for the ERS-1/2 and Radarsat data.
Engineering geological Unit 1, which is composed of recent and actual alluvial deposits of the Po River, shows the lower PS-DS density (~40 PS-DS/km2), and, for this reason, the maximum cumulated land subsidence may be due to local movements.
Units 2, 3, and 4 are characterized by cohesive soils (50%–60% of the thickness) in the first five meters, and show an average PS-DS density of ~100 PS-DS/km2. These units record a maximum cumulated land subsidence that varies from 100 mm to 200 mm in the first period, and an average value of 150 mm in the second one.
Unit 5 is composed of the alluvial fan of the Scuropasso River, with heterogeneous soils in the first 15 m of depth. As shown in
Figure 7, this unit shows the highest PS-DS density (~300 PS-DS/km
2) and the maximum cumulated land subsidence reaches ~51 mm for both periods.
Units 6, 7, and 8 are characterized by cohesive soils, predominant only in the first five meters of depth, and gravel and sand in the sub-soil. The PS-DS detected over these units is around 100 PS-DS/km2 and the maximum cumulated land subsidence reaches values in the ranges of 100–180 mm and 50–100 mm, respectively, for the periods 1992–2000 and 2003–2010. Note that the maximum cumulated uplift of 159 mm was detected for Unit 6, in the period 1992–2000.
In addition, comparisons between the engineering geological units and the areas of ground motion detected using the principal component analyses were carried out for the monitored periods (
Figure 7c,d). The results of the analyses show that the area of ground motion reaches around 23 and 17 km
2, respectively, in the periods 1992–2000 and 2003–2010. Therefore, the ground motion area decreases from the first period to the second one. It is worth noting that the study area is mainly characterized by PC3, which corresponds to seasonal deformation. Seasonal deformations may be due to seasonal aquifer compaction and/or swelling/shrinkage of clayey soils.
The PC1 ground motion that corresponds to the long-term lowering trend is mainly localized on Units 2 and 3 for both periods, while PC2 ground motion that corresponds to the uplift long-term trend is mainly localized on Units 1 and 6 for the entire monitored period.
6.2. Comparison between Ground Motion and Hydrogeological Settings
Another important factor that controls ground motion is the presence, and the thickness, of clayey soils, which may be susceptible to the consolidation process. A thickness map of the upper clayey-silty deposits was elaborated by Pilla et al. (2007) [
21] by interpolating the thickness data available from around 490 water wells, performed in the Oltrepo Pavese (
Figure 8a). This map allows us to analyze the relationship between the average LOS velocity and the thickness of the upper compressible deposits. The Oltrepo Pavese is characterized by the presence of a tabular clayey-silty cover, which increases from the northeast to the southwest. The thickness of the clayey-silty deposits ranges from 0 to 40 m, and the maximum thickness is observed in the southeastern sector of Voghera, and in the proximity of the towns Broni and Stradella (
Figure 8b,c).
Figure 8d,e shows the comparison between the average LOS velocity, detected by ERS-1/2 and Radarsat data, and the thickness of the upper clayey-silty deposits along the N–S and E–W sections (
Figure 8a). The results reveal that the peaks of average LOS velocity are consistent with the maximum thickness of clayey-silty deposits. As a consequence, the thickness of the upper clayey-silty deposits seems to be a controlling factor of the rate of ground motion, as observed in others case histories, such as Murcia, Vega Media of the Segura River Basin, and Alto Guadalentín Basin in Spain [
46,
47,
48], and along the Venice coast and Sibari plain in Italy [
49,
50].
Furthermore, the constant presence of a tabular clayey-silty coverage, which acts as a seal, is a crucial factor in determining the hydraulic condition of aquifer state. Indeed, the clayey-silty cover limits infiltration, and influences aquifer recharge, which occurs in correspondence with the coalescent fans originating from the deposition of the Apennine streams [
21]. Analysis of the hydraulic condition role on ground motion was performed using the unsaturated zone thickness, measured in October 2005 (
Figure 9). The unsaturated zone thickness detected in October 2005, was established by Pilla et al. (2007) [
21] using a piezometric level obtained from around 60 water wells across the study area. For the computation of the unsaturated zone, the piezometric levels were subtracted to the base of the clayey-silty cover. When the values are negative, the aquifer is confined and the flow is mainly horizontal. Conversely, when the values are positive, the aquifer is phreatic.
Figure 9 shows the overlapping of the ground motion areas mapped using PC3 of the Radarsat data with the hydraulic condition of the aquifer in October 2005. The results reveal that the seasonal ground motions (PC3), which are located over engineering geological Units 3 and 4, occur where phreatic condition of the aquifer was observed. The role of the seasonal fluctuations of the groundwater level on the ground motion areas is clear. Seasonal ground motions related to fluctuations in the groundwater level were previously recognized by other authors in Santa Clara Valley and Los Angeles Basin in the United States [
51,
52]. Meanwhile, the ground motion areas located over engineering geological Unit 8 are due to other processes, which are explained in
Section 6.3.
6.3. Comparison between Ground Motion Areas and Damaged Buildings
Most locations in the Oltrepo Pavese experienced damage to buildings due to swelling and shrinkage of clayey soils, respectively, during wet and dry periods. In most cases, the damaged buildings are single-story family residences, built on conventional shallow strip concrete footings. The damage was mainly detected over engineering geological Unit 8, where clayey soils susceptible to swelling/shrinkage were previously identified [
28]. Meisina et al. (2006) [
30] implemented a damaged building database.
Here, comparisons between the damaged building distribution and ground motion areas were performed. The results of the analyses highlight that the occurrence of damaged buildings is mainly localized over the ground motion areas detected using PC3. Therefore, this is consistent with the hypothesis that the seasonal deformational behavior detected by the third component of ground motion in engineering geological Unit 8 may be explained by the swelling/shrinkage of clayey soils.
In Codevilla (see location in
Figure 6), located over engineering geological Unit 8, the higher density of damaged buildings (12 damaged buildings in an area of around 0.6 km
2) corresponds to ground motion areas delineated using the PC3 of ERS-1/2 and Radarsat data (
Figure 10a). In this area, the shallow alluvial deposits of the first 8 m of depth exhibit swelling potentials from medium to very high [
30].
The swelling potential appears to be in good agreement with the displacement time series. Average LOS displacement time series was computed in the buffer zone of 100 m from the damaged building with code 425 in order to consider a significant number of measurements (around 20 PS-DS). The results reveal that the seasonal component of motion detected in the TS is directly correlated with the effective rainfall detected at the nearest pluviometric station, located in Voghera (
Figure 10a,b). In particular, significant damage, which consists of cracks, appeared in this building in August 1999 (crack width 5–10 mm). The pre-damage period from 1998 to 1999 was recognized as a drought period [
30], and the average effective rainfall from May to November of 49.45 mm corresponds to around −6.14 mm of ground motion. In February and March 2001, partial pile underpinning was carried out as a remedial measure for this building. In the following drought period, 2003–2007, the variability of the effective rainfall corresponds to there being no clear cyclic displacement time series. However, the average effective rainfall from May to November, of 62.95 mm, corresponds to around −1.27 mm of ground motion.
The results confirm that the ground motion areas detected at Codevilla via PC3 are due to swelling/shrinkage phenomena triggered by the variation of water content in these soils.
Overall, Codevilla seems to be characterized by an average LOS velocity range that is commonly defined as stable (±1.5 mm/yr.) during the periods 1992–2000 and 2003–2010; however, in this case, the seasonal deformation is significant, and provides evidence of damage, such as cracks in buildings, dislocations, and the breaking of walls, which were recognized by Meisina et al. (2006) [
30,
53]. In particular, the damaged buildings required the installation of micropile construction up to 8 m in depth, in order to reinforce structures.
6.4. Comparison between Ground Motion and Land Use Change Effects
In order to define the influence of land use change on ground motion, a visual inspection was performed by comparing orthophoto 1988 and 2000 [
54] and orthophoto 2003 and 2007 provided by the Lombardia Region, with PC1 ground motion areas. In addition, Corine Land Cover changes (CLC) from 1990–2000 and 2000–2006, implemented by the Institute for Environmental Protection and Research (ISPRA) [
55,
56] that match our satellite data, were exploited. These layers were created by using a working scale of 100,000; therefore, the land cover cannot be mapped in all its complexity, and the definitions of the units were established by the combination of surfaces, which are homogeneous to a greater or lesser degree, whatever the scale used. From the total of the ground motion areas delineated by PC1 of the ERS-1/2 and Radarsat data, only two and three areas seem to be influenced by land use changes, respectively. More precisely, in the period 1992–2000, land use change from non-irrigated arable land to discontinuous urban fabric, which occurred over engineering geological Units 1 and 2, matches with PC1 ground motion areas detected by ERS-1/2 data. In the period 2000–2006, land use change from non-irrigated arable land to construction sites occurred over engineering geological Unit 3, which matches the PC1 ground motion area detected using Radarsat data.
It worth noting that these ground motion areas, detected during the period 1992–2000, are not present in the period 2003–2010. Therefore, the stress increment produced by a building’s load seems to have dissipated in the first period.
Figure 11 shows ground motion detected by the PC1 of ERS-1/2 data in Voghera. The average displacement time series of the PS-DS in the ground motion areas highlight the decrease of the average LOS velocity from −3.15 mm/yr during 1992–2000 to −0.80 for 2003–2010 (
Figure 11c).
7. Conclusions
A challenge to ground motion studies is the disentanglement of different phenomena and to assess their relative contributions. Indeed, specific causes of deformation phenomena in lowlands are presently unknown due to the interaction of either deep (“geological”) or superficial causes (withdrawal of water and gas, soil consolidation, load-induced compaction), which take place at different spatio-temporal scales and of which the mixed effects must be decorrelated. This paper demonstrated the ability of a novel methodology to detect and disentangle the contributions of different geohazards, by analyzing multi-sensors and multi-temporal A-DinSAR data. The principal novelty of the proposed approach is the exploitation of a full time series, which enables the detection of ground motion areas affected by linear, non-linear, and seasonal movements. The proposed methodology consists of three phases:
Post-processing elaborations; disentanglement of the vertical and horizontal components of motion and displacement time series accuracy assessments, performed by selecting the most coherent (>0.9) _targets with an average LOS velocity in the range ±0.5 mm/yr, where no significant movements are expected.
Identification of ground motion areas; application of two approaches for automatic A-DInSAR time series analyses. First, principal components analysis of the time series is performed, implementing a matrix organization of location versus time (T-mode). Then, PS-time software is used to find the kinematic behavior of the time series (i.e., uncorrelated, linear, and non-linear trend). Finally, ground-motion areas were identified using a spatial cluster with a buffer zone of 50 m around the _targets, characterized by a significant correlation with the principal components of motion. The results of both approaches are overlapped to distinguish linear and non-linear ground motion areas.
Mechanism recognition: Comparison between ground motion areas and ancillary data (i.e., geological, geotechnical, hydrogeological, and land use information) in order to interpret the driving forces of the subsidence mechanisms.
The methodology was tested in Oltrepo Pavese (Italy), which is a representative site of moderate subsidence of the Po Plain, using ERS-1/2 and Radarsat data covering a time span of around 20 years. Results confirm its suitability for the detection of the temporal evolution of deformation patterns due to different geohazards, not recognizable by conventional analyses of the average LOS velocity alone (i.e., seasonal components of motion). Three deformational behaviors were detected: linear, non-linear (accelerations and decelerations of the movements), and seasonal components of motion. The comparison between a kinematic model of the ground motion areas and hydrogeological and engineering geological data allow for an understanding of the causes of motions that are often superimposed. The main natural and man-induced processes are swelling/shrinkage of clayey soils, land subsidence due to the load of new buildings, moderate tectonic uplift, and seasonal ground motion due to seasonal groundwater level variations. Overall, the cumulated displacements observed in the 1992–2000 period are higher than those detected in the period 2003–2010, and a decrease of deformation was observed. The findings in this work confirmed that the combination of A-DInSAR time series analysis and hydrogeological data is capable of supporting knowledge of the relationship between groundwater level fluctuations and deformations, and of providing fundamental information about the aquifer state conditions.
The proposed methodology can be easily applied to other geological contexts, using different A-DInSAR data. The reliability of the results depends on the quality of the dataset. New opportunities will be provided by the improvement of the time series acquired by new sensors of the COSMO-SkyMed and Sentinel missions.
The procedure will also allow for specifying the priority area of where to address prevention activities, in order to optimize the costs and benefits and to draw a management plan of land use and groundwater resources, at national and regional scales. The procedure may also be used in land subsidence studies to discern the contributions of the seasonal components of motion from others processes, such as soil consolidation and solid and fluid extraction or injection. The approach will be tested, in the future, in geological contexts characterized by other geohazards, such as landslides.