Statistical and Independent Component Analysis of Sentinel-1 InSAR Time Series to Assess Land Subsidence Trends

Abstract

Advanced statistics can enable the detailed characterization of ground deformation time series, which is a fundamental step for thoroughly understanding the phenomena of land subsidence and their main drivers. This study presents a novel methodological approach based on pre-existing open-access statistical tools to exploit satellite differential interferometric synthetic aperture radar (DInSAR) data to investigate land subsidence processes, using European Ground Motion Service (EGMS) Sentinel-1 DInSAR 2018−2022 datasets. The workflow involves the implementation of Persistent Scatterers (PS) time series classification through the PS-Time tool, deformation signal decomposition via independent component analysis (ICA), and drivers’ investigation through spatio-temporal correlation with geospatial and monitoring data. Subsidence time series at the three demonstration sites of Bologna, Ravenna and Carpi (Po Plain, Italy) were classified into linear and nonlinear (quadratic, discontinuous, uncorrelated) categories, and the mixed deformation signal of each PS was decomposed into independent components, allowing the identification of new spatial clusters with linear, accelerating/decelerating, and seasonal trends. The relationship between the different independent components and DInSAR-derived displacement velocity, acceleration, and seasonality was also analyzed via regression analysis. Correlation with geological and groundwater monitoring data supported the investigation of the relationship between the observed deformation and subsidence drivers, such as aquifer resource exploitation, local geological setting, and gas extraction/reinjection.

1. Introduction

Land subsidence is a gradually evolving, but often very impactful, geological phenomenon influenced by various natural and anthropogenic factors [1,2]. Recent advancements in remote sensing technologies, especially satellite interferometric synthetic aperture radar (InSAR) [3], have significantly enhanced our capability to monitor and analyze subsidence magnitude, patterns and temporal trends. Particularly, differential InSAR (DInSAR) is a technique that exploits the information contained in the radar phase of at least two complex SAR images acquired at different times over the same area to identify and estimate any terrain displacement potentially occurring between the satellite passes [4]. To discriminate the displacement from other phase components (such as topographic errors, atmospheric effects, decorrelation noise), as well as to generate displacement time series, different advanced DInSAR data processing and analysis workflows have been developed [5,6,7,8,9,10]. For instance, [5] effectively separates displacement signals by analyzing pixels with low phase noise, considering only targets coming from a strong reflecting object constant over time (called Permanent or Persistent Scatterers, PS; PSInSAR technique), or including those targets where scattering is constant over time, but comes from different small scattering objects (Distributed Scatterers, DS; SqueeSAR technique) [6]. The accuracy of advanced DInSAR datasets against independent geodetic data, such as leveling and global navigation satellite systems (GNSS), can reach a few mm for displacement records and less than 1 cm/year for displacement velocities [11,12,13,14,15].

The increasing SAR data availability, the development of new advanced DInSAR data processing techniques, and the increase of computational capability has enhanced the capability of DInSAR analysis to cover wide areas, up to the continental scale [16,17,18]. It is in this context that the European Ground Motion Service (EGMS) was developed, a wide-area deformation monitoring system based on advanced DInSAR techniques, part of the Copernicus Land Monitoring Service [16,19]. It builds upon big data stacks of C-band Sentinel-1A/B Interferometric Wide (IW) swath SAR images, processed by the combination of advanced PS and DS techniques (see [5,6,19] and the authors cited therein), and its time series datasets have recently started to be exploited by the scientific community for urban instability and geohazard investigations using semi-automated tools [20].

On the other hand, different methodologies have been developed in recent years to analyze the resulting multi-temporal DInSAR time series data, aimed to the geological interpretation of areas affected by land subsidence/uplift and seasonal movements [21,22,23,24,25]. By employing advanced statistical and machine learning tools, these methodologies succeed in identifying distinct spatio-temporal trends that characterize land deformation. Distinguishing between different components of motion, such as long-term and seasonal behaviors, is fundamental to analyzing the different natural and anthropogenic processes that may be contributing to the subsidence process. For example, [26] employed a principal component analysis (PCA) to extract the spatial patterns of deformation embedded in an InSAR time series and related them to land elevation changes caused by groundwater withdrawal and aquifer system recharge. On the other hand, [27] characterized areas of accelerating and decelerating displacement behavior, identified zones of sagging and hogging using horizontal displacement data, and marked areas with potentially yet unmapped ground discontinuities by analyzing deformation trends on structurally-controlled subsidence.

Some of the methodologies are based on well-known statistical algorithms, that were not originally designed for InSAR data evaluation but proved to be useful for deformation time series characterization. This is the case with the independent component analysis (ICA), an advanced statistical method that decomposes a mixed signal into a linear combination of many independent components. ICA has been used with InSAR data for the evaluation of volcanic deformation signals [28,29], to correct atmospheric phases [30], and can be used to retain the different components of motion that could be embedded in a deformation time series [31]. Conversely, some authors focused on developing algorithms specifically for DInSAR time series evaluation and subsidence characterization. In this framework, [21] developed a semi-automated classification tool that performs a series of statistical tests to classify different deformation time series into predefined categories of motion. This methodology has proved effective for detailed landslide and subsidence analysis [27,32] and was even included in combined methodological approaches proposed by other authors [23,24]. However, the full potential of these statistical methodologies for a complete evaluation of subsidence processes is yet to be explored.

With the aim to contribute to this research line, the scientific goal of this study is to bring novelty to the field by establishing and evaluating for the first time the benefits of the sequential application of the following: (1) a semi-automated statistical test of displacement time series for trend classification and mapping; and (2) an ICA for trend and pattern retrieval to assess land subsidence and identify spatially the different ground motion trends affecting a given study area. In addition to the statistical analysis of deformation time series, a spatio-temporal correlation with geological and groundwater monitoring records is performed for drivers’ identification. The proposed methodological approach has been applied on three representative areas prone to land subsidence in the Po Plain, a sedimentary basin in northern Italy: Ravenna, Bologna, and Carpi, which are well known in the scientific literature for the subsidence phenomena affecting them over the last decades [33]. As such, they offer the ideal testing ground for the approach, given the associated variety of temporal trends, including accelerating, decelerating, seasonal, quadratic and linear behaviors, and a combination of natural and human-related drivers. Although each method has its own limitations, a combination of them that is experimented in this work results in a powerful and novel analytical methodology that enhances a comprehensive analysis of the subsidence phenomena necessary for land and urban planning, and the implementation of risk mitigation strategies. The application of the method for the characterization of the three land subsidence study sites enabled a better understanding of the phenomena, as well as the identification of the potential and limitations of the proposed workflow, and key aspects and directions for future research in the field of land subsidence characterization with advanced statistics and computational tools, using free and open data such as the EGMS products.

2. Study Areas

Ravenna, Bologna, and Carpi are located in the Po Plain, the most extensive (~46,000 km²) and densely populated (~450 inhabitants/km²) plain in Italy (Figure 1). Located in the northern part of the country, the plain is crossed by main roads and highways, and many urban centers and industrial settlements that are widely developed across a vast portion of its territory, while the rest of it is characterized by intensive farming activities [34,35]. The plain spreads west to east from the western Alps to the Adriatic Sea, representing the foreland sedimentary basin system of two fold-and-thrust belts (Northern Apennines and Southern Alps) [36,37]. The upper part of its stratigraphic sequence is characterized by the presence of Middle Pleistocene–Holocene alluvial fine and unconsolidated sediments that were deposited by the Po River and its tributaries, and consist of an alternation of sands, silts, and clays [38]. This sedimentary sequence hosts a complex multi-aquifer system, which has been pumped through the years [39].

A natural subsidence of 2.5–3 mm/year was recognized in the Po Plain due to tectonics and sediment load/compaction [40,41]. Other studies have identified higher rates of up to 5 mm/year through the geological history of the sedimentary basin, observing the additional effect of short-term glacial cycles to the long-term processes (tectonics, geodynamics, and sediment load/compaction) [42]. On the other hand, significant anthropogenically-driven subsidence, with rates over an order of magnitude higher than that associated with natural causes, has been recognized in the plain for years, and has been principally attributed to intense groundwater pumping [39,43]. For example, the last report published by the Regional Agency for Prevention, Environment and Energy of Emilia-Romagna (ARPAE)’s subsidence monitoring project [33], which encompasses deformation phenomena over the 2016−2021 period for the Emilia-Romagna region, reveals that over 6% and 3% of the regional territory experience subsidence rates from −5 to −10 mm/year and from −10 to −20 mm/year, respectively, attributed to anthropogenic causes [33]. The selected study sites of Bologna, Ravenna, and Carpi have been historically recognized among the main areas substantially affected with anthropogenically-driven deformation over the Po Plain [43] and were also identified as key subsidence hotspots by this report [33].

Figure 1. The study areas of (1) Ravenna, (2) Bologna, and (3) Carpi–Correggio–Soliera: (a) geographical location in Italy; (b) extent of the European Ground Motion Service (EGMS) Level-3 (L3) and Level-2b (L2b) dataset footprints used for the statistical analysis, overlapped onto the Copernicus Global Digital Elevation Model [44]; and (c) detail of the mean vertical deformation velocity from EGMS L3 datasets, overlapped onto a Google satellite imagery basemap.

Figure 1. The study areas of (1) Ravenna, (2) Bologna, and (3) Carpi–Correggio–Soliera: (a) geographical location in Italy; (b) extent of the European Ground Motion Service (EGMS) Level-3 (L3) and Level-2b (L2b) dataset footprints used for the statistical analysis, overlapped onto the Copernicus Global Digital Elevation Model [44]; and (c) detail of the mean vertical deformation velocity from EGMS L3 datasets, overlapped onto a Google satellite imagery basemap.

3. Materials and Methods

3.1. DInSAR Datasets

The input InSAR datasets used in this study include the most updated EGMS version currently released by Copernicus as of mid-2024, which provides ground displacement data for a time span from January 2018 to December 2022.

Level 3 (L3) “Ortho” products containing vertical and east-west displacement information, resampled to a 100 m × 100 m grid, were used (Figure 1b). These are obtained by combining ascending and descending orbit InSAR datasets, and assuming that no north-south displacement had occurred. A total of 6 L3 100 km × 100 km tiles were downloaded from the Copernicus EGMS Explorer portal, mosaicked and clipped to the extent of the three demonstration sites of Bologna, Ravenna, and Carpi (Figure 1c).

To detail the analysis at site scale, point-wise Level 2b (L2b) “Calibrated” products containing PS datasets with line-of-sight (LOS) annual deformation velocity, time series, and acceleration (among other values) were assessed. A total of five Sentinel-1 IW bursts were examined covering the three study sites (Figure 1b): three descending datasets to cover Bologna and Carpi, and two ascending datasets for Ravenna. The datasets involve, respectively, Sentinel-1 track 168 with IW sub-swath 1 (incidence angle varying from ~30° to ~36°, from the near to the far range) and track 117 with sub-swath 2 (incidence angle from ~36° to ~41°). These datasets are calibrated with a model derived from high-quality GNSS data and they are not relative to a local reference point (i.e., they display absolute displacement values). Mean annual velocity and acceleration values are calculated by the EGMS production team on the time series residuals after applying a regression model of a first and second order polynomial, respectively, plus a seasonal (sinusoidal) component [45]. Positive LOS displacement rates indicate movements towards the satellite sensor (or uplift, in case of vertical rates), while negative LOS rates indicate movements away from the sensor (or subsidence, in case of vertical rates). Likewise, positive acceleration values at subsiding sites indicate the deceleration of the subsidence process (i.e., an increase in the displacement rate towards positive values), whereas negative acceleration at subsiding sites indicates the transition towards lower negative values, hence an acceleration of the subsidence process.

To enable a comparison of L2b time series with groundwater monitoring records, LOS deformation estimates (D_LOS) contained in the EGMS L2b products were reprojected into vertical displacements (D_U), by using the simplified assumption of absence of horizontal displacement, and therefore by dividing D_LOS by the local value of the directional cosine along the vertical direction:

$$D_U={\displaystyle \frac{D_{LOS}}{cos\theta }}$$

(1)

where θ is the incidence angle of the LOS. The assumption was deemed suitable for the vast majority of the study areas, given that east-west displacements were generally absent according to the EGMS east-west ortho product or, when present, they are characterized by a lower order of magnitude with respect to vertical displacement records.

3.2. Groundwater and Geological Datasets

Groundwater monitoring data were provided by ARPAE, which oversees surveys at a network of groundwater monitoring wells over the whole Emilia-Romagna region, dividing them according to their use (e.g., irrigation, industrial, domestic). The data include histories of piezometric levels (height of the groundwater level above the datum) and groundwater depths to the topographic level at each well, which are manually recorded by ARPAE’s team during regular surveys at the regional network (using pressure transducers, electric probes or phreatimeters) or, at some piezometers, registered hourly with advanced equipment installed on site. The depths of the different piezometers and their filters vary (hence they allow recording data for different layers of the aquifer system), and some of them also feature open screens across multiple layers.

On the other hand, geological datasets used for the statistical correlation of ground deformation with lithology were downloaded from the MinERva Portal, also managed by the Emilia-Romagna Region service [46]. These include surface geology maps at the 1:25,000 scale, depicting depositional environments and a lithological description. Sedimentary deposits are categorized according to their texture in gravel, sand, silt, and clay, including mixtures of the different granulometries (e.g., silty sand).

3.3. Time Series Semi-Automatic Classification (PS-Time)

“PS-Time” is a MATLAB-based semi-automatic classification tool developed by [21] that recognizes different motion trends in the input InSAR time series through the application of a conditional sequence of statistical tests, such as multiple regression and Analysis Of Variance (ANOVA). The program has been used for the evaluation of DInSAR time series in many different geohazard applications, proving to be significantly useful for land movement characterization [24,27,47,48]. It is included in the methodology for detection of ground motion areas proposed by [23], who succeeded in identifying different components of motion in the Oltrepò Pavese (Po Plain, Italy) related to natural and man-induced processes (such as swelling/shrinkage of clayey soils, land subsidence due to load of new buildings, and seasonal behaviors in response to groundwater level variations).

The tool conducts a series of statistical tests to categorize the time series into six predefined categories: uncorrelated, linear, quadratic, bilinear, discontinuous with constant velocity, and discontinuous with variable velocity. The latter four categories are collectively referred to as “non-linear” motion trends. Additionally, it evaluates the presence of any cyclic fluctuations in the displacement time series, characterized by a periodicity of around one year, by computing its amplitude.

As explained in detail in [21], PS-Time first runs an ANOVA F test to assess the significance of the linear regression. If there is no relationship between displacement and time, the slope of the linear regression β1 equals 0 (i.e., null hypothesis). If the probability of such condition (p1) exceeds the level of significance of the linear regression (α1) that is set by the user (p1 > α1), the point is classified as “uncorrelated”. Otherwise, if p1 ≤ α1, the null hypothesis is rejected, and the existence of a linear relationship is confirmed. That being the case, a segmented linear regression is run to assess the probability of a bilinear motion trend. The Bayesian information criterion (BIC) [49] is exploited to verify the existence of a breakpoint, and the evidence ratio B_w [50] is calculated and compared with the selected level of significance of the segmented regression (B_th). If the segmented regression is not significant (B_w < B_th), the time series is further analyzed using quadratic regression and another ANOVA F test. If the probability of not having a quadratic term (p12) does not exceed the predefined threshold (α12) that was set by the user (p12 ≤ α12), the deformation behavior of the point is classified as “quadratic”. Otherwise, when p12 > α12 the “linear” motion trend is confirmed. On the other hand, if B_w ≥ B_th, a breakpoint is identified, and the series is classified as “bilinear”.

At that stage, the presence of any vertical jump in the series is also evaluated by a discontinuity test with confidence bands at 95%, and an F test for the equality of slopes [51] is applied to identify any regression slope change before vs. after the breakpoint (which would classify the series as “discontinuous”, with either constant or variable velocity).

The outcomes of the automatic classification tests are highly influenced by the chosen significance level for each statistical analysis, which can be determined before executing the program. In this study, we set such thresholds to α1 = 0.05, B_th = 1.2 and α12 = 0.01, by building upon both the recommendations of the developers [21] and the results of the extensive trials presented in [32]. The latter tested four different combinations of the three thresholds (α1, B_th, α12) to classify a ~5500 PS time series extracted from the EGMS InSAR dataset for the city of Ravenna as follows: (i) α1 = 0.01, B_th = 1 and α12 = 0.01; (ii) α1 = 0.05, B_th = 1.2 and α12 = 0.01; (iii) α1 = 0.05, B_th = 1.5 and α12 = 0.05; and (iv) α1 = 0.05, B_th = 1.2 and α12 = 0.05. By training the tool with an assorted selection of linear, quadratic, bilinear, and uncorrelated time series, set (ii) provided the best classification performance.

3.4. Independent Component Analysis (ICA)

ICA is a blind signal separation method that consists of decomposing a mixed signal into a linear combination of statistically independent sources, under the assumption that each component has a non-Gaussian probability distribution [52,53]. The central limit theorem underpins this assumption, suggesting that a mixture of non-Gaussian signals tends towards a Gaussian distribution, making it feasible to identify the original non-Gaussian sources by maximizing their statistical independence [52,53]. Therefore, the problem is expressed mathematically as:

$$X_{(t\ast p)}=A_{(t\ast n)}.S_{(n\ast p)}$$

(2)

where X is the input mixed signal matrix, with the number of dates (t) considered in the rows and the number of pixels/points (p) as columns, S is the matrix of independent components (IC), and A is the mixing matrix, which describes the relative contribution of each source to the input mixed signal. The number of IC retrieved is n.

The authors of [52] introduced the FastICA algorithm, which performs a fast fixed-point iteration to retrieve independent sources in the mixed signal. The algorithm starts by centering and whitening the observations, achieved by subtracting the mean from the mixed signal matrix, making the variables zero-mean. This converts the mixing matrix into an orthogonal matrix, reducing the number of free parameters. FastICA preconditions the centered observations using principal component analysis (PCA) [54], which is also helpful for reducing the dimensionality and noise of the data. Following this, the mixed signals are linearly transformed into uncorrelated variables with a variance of 1 through whitening, and the problem becomes U = D∙S, where U is the centered and whitened mixed matrix and D is the orthogonally adjusted mixing matrix. Consequently, the source matrix can be estimated using the equation S = D⁻¹∙U, where D⁻¹, also known as the unmixing matrix W.

The number of principal components to be retained for the ICA, which is consequently equal to the number of resulting ICs, can be defined by the user while running the FastICA algorithm. In this work, a trial-and-error approach was applied to select the minimum number of components necessary to accurately describe the different deformation patterns embedded in the total deformation signal, trying to avoid overfitting (e.g., [28]). Starting from two, the number of components was iteratively increased in each trial until no more significant new trends were recovered. A temporal ICA was applied to the deformation time series through the freely available FastICA MATLAB code, distributed by Aalto University [55]. The output number (unique identifier) assigned to each component was chosen randomly by the program, without following a specific order [30] (e.g., the linear IC could be assigned the identifier IC1, IC2, or even a different one, not necessarily placed at the beginning of the IC list).

4. Results

4.1. Subsidence Time Series Trend Classification

The vertical velocity provided by L3 EGMS datasets at 100 m resolution (Figure 1c), and the mean LOS velocity and acceleration values provided by L2b EGMS point-wise products enable a first-level assessment of the minimum and maximum values estimated at the three study sites as follows: (1) in the coastal area of Ravenna, LOS subsidence rates are lower than –5 mm/year for ~30% of the PS targets, with a peak at −37.1 mm/year and vertical rates as low as −50.0 mm/year, and accelerating trends reach up to −7.1 mm/year²; (2) within the wide area comprising Bologna, LOS velocities are lower than −5 mm/year for ~40% of the PS targets and reach up to −35.8 mm/year corresponding to −42.5 mm/year vertical rates, and deceleration and acceleration are up to +7.6 and −8.2 mm/year², respectively; and (3) in the industrial area of Carpi–Correggio–Soliera, another significant mean deformation velocity cluster is identified, with LOS values lower than −5 mm/year for ~60% of the PS targets, with a peak at −25.2 mm/year and vertical rates of −25.4 mm/year, and acceleration/deceleration values reaching −5.7 and +5.6 mm/year², respectively.

To fully understand and evaluate the spatial variability of subsidence during the five year-long period, LOS deformation velocity, acceleration, seasonality, and trend classification maps were created and analyzed in detail for each of the three study sites (Figure 2, Figure 3 and Figure 4). The first two parameters were extracted from the L2b EGMS products, while the seasonality (i.e., amplitude of the annual periodicity component) and the time series trend classification were obtained from the PS-Time analysis.

In the area of Ravenna (Figure 2a), the LOS deformation velocity map displays values higher than −5.0 mm/year towards the city center, aligning with the natural subsidence rates revised in the literature [43] and indicating a general stability in that area. On the other hand, lower velocities are observed in the western agricultural fields and over the coastal zone, with values varying between −5.0 and −10.0 mm/year. The industrial district and a small area identified by Copernicus’ Urban Atlas Land Cover inventory [56] as a dump site (see location in Figure 2a) exhibit the lowest deformation velocities, with peaks of −28.0 and −37.1 mm/year, respectively.

The acceleration map (Figure 2b) highlights a cluster of decelerating points in the southern coastal region, with values between +3 and +6 mm/year². Conversely, the area of the agricultural fields depicts an accelerating tendency, with values up to −7.1 mm/year², as well as some annual periodic components, with amplitudes varying between 2 and 8 mm (Figure 2c). The seasonality map (Figure 2c) also shows a small cluster of points with amplitude exceeding 6 mm within the industrial area of Ravenna, particularly noticeable over a large fabric metal rooftop. This seasonal pattern may be attributable to the possible thermal expansion of the material in response to seasonal temperature changes.

Finally, the PS-Time trend analysis classified the majority of the points as quadratic (~67%) or linear (~32%), with less than 1% falling into other categories (uncorrelated, bilinear, or discontinuous). A cluster of quadratic trends was noted exclusively in the agricultural fields zone (Figure 2d), aligning with the accelerating cluster previously mentioned (Figure 2b).

The maps from the area of Bologna show stability only in the southern area of the city center that is towards the Apennine border, with LOS deformation velocity values between −3 and +3 mm/year, among the natural subsidence values registered for the region. On the other hand, a notorious concentric subsidence cluster is observed in the center and northern part of the study area (Figure 3a), with the highest deformation rates (−35.8 mm/year) recorded along Reno river plain. The general mean deformation velocity cluster broadly coincides with the decelerating cluster shown in Figure 3b, with values generally in the range of +1 to +3 mm/year², and a localized peak of +7.6 mm/year² in the town of San Giovanni in Persiceto. On the other hand, two accelerating clusters were identified, one located in the western sector of the study area, around the industrial area of Ponte Samoggia, and another towards the eastern part, between Ozzano dell’Emilia and Minerbio towns, with acceleration values between −3.0 and −8.2 mm/year². The sector of Minerbio is also characterized by a strong seasonal component, as can be seen in Figure 3c, representing the only seasonally deforming cluster recognizable in the study area.

Similar to what was observed in the area of Ravenna, PS-Time mainly identified linear (~50%) and quadratic (~44%) trends, with clusters of quadratic points aligning with the most accelerating and decelerating areas (Figure 3d). In addition, points located over the Apennines boundary are mostly classified as uncorrelated (~6%).

The industrial area of Carpi, Correggio, and Soliera displays a general instability, also concentrically distributed, with subsidence velocities ranging between −3.0 and −20.0 mm/year (Figure 4a). The lowest deformation rate (−25.2 mm/year) is located along the A22 highway, in the north-western sector of Carpi. However, this area also experiences the strongest decelerating tendency, of around +5.6 mm/year², according to the results observed in the acceleration map (Figure 4b). On the other hand, a strongly focalized accelerating cluster is shown over the town of Soliera, with a peak value of −5.7 mm/year² (Figure 4b). No seasonally deforming clusters were recognized in the area (Figure 4c).

The trend classification map shows a general linear (~59%) and quadratic (~40%) time series behavior, with a cluster of quadratic points coinciding with the area of Carpi where the highest subsidence velocities and decelerating values were observed (Figure 4d). Only less than 1% of the time series exhibits different behavior. Among them, a cluster of bilinear points (0.14%) was identified over Soliera, aligning spatially with the highest accelerating area identified in Figure 4b. For these points, PS-Time software v.1.1 (64-bit) identified breakpoints in the time series between May 2021 and December 2021 (see Figure A1 in Appendix A).

4.2. Retrieval of Different Deformation Patterns Through ICA

Taking into consideration the results obtained using PS-Time, local scale analysis has been performed with the FastICA algorithm to evaluate the different deformation patterns and trends embedded in the deformation time series, as well as to detail the understanding of their spatial distribution, over the three study sites. Three small subsets of ~22 km² were selected, one within each study area (the industrial area of Ravenna, and the areas around Minerbio and Soliera; Figure 2a, Figure 3c and Figure 4a), plus a wider 52 km² testing sector north-west of Bologna (“Area 2”, between San Giovanni in Persiceto and Ponte Samoggia; Figure 3b). These sub-areas were selected according to the places where the most variable time series behaviors were observed, in order to fully assess their deformation pattern (see Section 4.1).

Figure 5, Figure 6 and Figure 7 depict the various independent signals identified for each area (i.e., namely, the rows of the mixing matrices, A in Equation (2)), along with the spatial distributions of the components (i.e., S in Equation (2)). The latter is represented with a color scale in the form of “scores”, which describe how similar the time series of a certain PS point is to the one retrieved in a certain component: higher scores indicate a higher match between the time series, while scores approaching 0 indicate no correlation, and negative scores highlight that the PS follows a trend that is opposite to the one depicted by that component. To simplify the following discussion, they are codified with the following structure: initials of the area (Ra, Ravenna; Bo, Bologna; So, Soliera) and number of IC (1, 2, 3, 4). The different independent signals mixed in the total ground motion time series generally showed the presence of a linear component, a quadratic one, plus one or two sinusoidal signals. For Ravenna, San Giovanni in Persiceto and Ponte Samoggia (Bologna, Area 2), and Soliera four components were retained (Figure 5, Figure 6a and Figure 7), while at Minerbio (Bologna, Area 1; Figure 6b) the deformation patterns were described by using two components only. This methodological choice enabled the identification of the most significant components, but included some noisy signals (e.g., at Soliera, So-IC4; Figure 7), and less relevant or mixed components, which are discussed in detail below.

By analyzing the time series of the various ICs within the four testing areas, linear deformation trends were identified in Ra-IC2, Bo1-IC1, Bo2-IC1, and So-IC3, corresponding to the areas of Ravenna, Bologna (1 and 2), and Soliera, respectively. On the other hand, Ra-IC4, Bo2-IC3, and So-IC1 (and, less pronouncedly, Bo1-IC1) depict quadratic deformation signals. Two sinusoidal signals of different periodicity were retained for Ravenna in Ra-IC1 and Ra-IC3 (with a yearly, and a multiyear frequency, respectively), whilst just one periodic signal was identified for the area of Soliera in So-IC2. Sinusoidal components were also retained in Bo1-IC2, Bo2-IC2, and Bo2-IC4 for Areas 1 and 2 of Bologna, although the latter component mentioned is significantly affected by noise and shows a slight shift in the temporal location of the peak values of the sinusoid with respect to Bo2-IC2.

When displaying the spatial distribution of the ICs (maps in Figure 5, Figure 6 and Figure 7), it is interesting to notice that there is a great spatial correspondence in most of the cases with the clusters observed in the mean deformation velocity, acceleration, and seasonality maps (Figure 2, Figure 3 and Figure 4). The areas that display high scores for linear components align with the ones where lower deformation velocities were recorded (i.e., Figure 3a with Bo1-IC1 and Bo2-IC1, Figure 4a with So-IC3, Figure 2a with Ra-IC2). Something similar happens with quadratic components, whose scores generally align with the spatial distribution observed in acceleration maps (i.e., Figure 3b with Bo2-IC3, Figure 2b with Ra-IC4, Figure 4b with So-IC1). This is the case for Ra-IC4 and So-IC1, which present a quadratic trend that scores positively in the areas where accelerating clusters were previously identified (see Figure 2b and Figure 4b). On the other hand, although the trend retrieved by the mixing matrix for Bo2-IC3 is flipped, positive IC3 scores are observed for the accelerating area of Ponte Samoggia, and strong negative scores are found for the decelerating areas of San Giovanni in Persiceto and Anzola dell’Emilia (Figure 3b). Therefore, it is reasonable to think that negative scores indicate time series trends that behave exactly oppositely to the signal retained in the component.

Sinusoidal components are more variable within the observed time series. Although Ra-IC1 aligns spatially with the seasonality amplitude cluster identified in Figure 2c, Ra-IC3 shows no clustering or specific correlation with any parameter determined before. Consequently, it can be assumed as a not relevant, noisy component. A similar situation occurs in Area 2 of Bologna, where Bo2-IC4 represents a noisy component, with no apparent spatial clustering on the map. However, it is interesting to notice that Bo2-IC2 depicts a cluster of points with seasonal time series that was not recognized in the area using the EGMS dataset information or PS-Time (Figure 3b,c), therefore suggesting the importance of implementing ICA to assess similarly complex trends or any hidden patterns. As an example, the time series of one of these points is shown in Figure A2 in Appendix A, where some very noisy seasonality can be observed onto a seeming quadratic behavior. In the area of Soliera, So-IC2 represents the seasonal component, but no spatial clustering is recognized, coinciding with the results observed in Figure 4c.

To confirm the relationship between the different components retained for each area and the mean deformation velocity, acceleration, and seasonality parameters previously analyzed, a regression analysis was pursued and R² values were calculated (Figure 8, Figure 9, Figure 10 and Figure 11). Linear fittings were performed for mean deformation velocity and acceleration correlations, while a quadratic fitting was applied for the amplitude of the annual periodicity component (APC) parameter. In all cases, the components that displayed linear trends (Ra-IC2, Bo1-IC1, Bo2-IC1, and So-IC3) are accurately correlated with the mean deformation velocity measured in each area, showing a negative linear correlation, with R² values ranging between 0.517 and 0.857. Similarly, quadratic ICs (Ra-IC4, Bo2-IC3, and So-IC1), which spatially seemed to correlate with the accelerating and decelerating trends observed on the maps, show a clear negative linear correlation with the accelerating values calculated from EGMS products, with slightly lower R² values, ranging between 0.518 and 0.709.

When comparing the seasonal ICs with the annual seasonality calculated through PS-Time, points are generally displayed following a parabolic distribution with axes of symmetry at IC scores equal to 0 (e.g., Bo1-IC2 and Bo2-IC2 in Figure 9 and Figure 10, and So-IC2 in Figure 11). Similarly, evidence of a bilinear correlation is displayed for Ra-IC1 in Figure 8, with breakpoint at a score of 0. The parabolic and bilinear correlations in this case are believed to happen because some points score positively to a certain periodic pattern retained by ICA, as they follow exactly the same trend; but others, although being very seasonal, show the opposite periodicity (e.g., positive peak in the deformation time series, corresponding to a negative peak in the IC mixing matrix series).

4.3. Subsidence Drivers Analysis

4.3.1. Ravenna

Historically, subsidence in Ravenna has been linked to the combined effects of groundwater extraction and methane withdrawal from both onshore and offshore gas reservoirs, which leads to pressure loss and subsequent compaction [40,41,57]. For the 2018−2022 period, Porto Corsini Mare and Angela Angelina offshore fields are the only gas extraction fields with ongoing operative concessions in the study area, according to the data published by the Italian Ministry of Environment and Energy Security [46]. In order to explore a possible connection between the Angela Angelina reinjection well and the deceleration pattern observed in Figure 2b, changes in acceleration rates within a 3 km buffer around the well were analyzed (Figure A3, Appendix A). The area closest to the well exhibited the most significant deceleration tendency (i.e., +1.64 mm/year² for that point, equivalent to declining subsidence velocity over time), and lower deceleration rates were observed as the distance from the well increased. These results suggest a possible correlation between gas reinjection activities and the localized decrease in deformation velocities near the well.

Comparison of piezometric records at ARPAE’s monitoring wells RA49-00 and RA29-00 with InSAR-derived ground deformation in contiguous PS–DS points provides insights into the correlation between groundwater withdrawal and accelerating subsidence at the agricultural fields and over the coastal area (Figure 12b). RA49-00 is drilled in the alluvial plain, and reaches the lower confined aquifer, with filter at 205−229 m; whereas RA29-00 is drilled in the coastal alluvial plain and supplies water for irrigation, from the confined aquifer, with filter at 233−251 m. In the first case, some correlation between the deformation time series of a representative PS point in the vicinity of the well and the measured piezometric change was observed, despite the limited piezometric and displacement changes recorded (i.e., an 80 cm decrease in underground water levels, accompanied by ~20 mm cumulative deformation during 2018−2022). Additionally, seasonal variations are observed in the displacement time series (with ~1.0 mm of amplitude according to PS-Time, overlapping onto a noisy quadratic behavior), though they are less distinguishable in the piezometric records, given the less frequent temporal sampling), that may suggest the existence of an elastic deformation behavior overlapping onto the main process of compaction. Peak values occur in winter and spring (potentially due to rainfall-driven aquifer recharge and short-term deceleration of the subsidence process, e.g., in January–May 2022) and are followed by acceleration in summer and autumn (e.g., in July–September 2022). A different situation is observed in RA29-00, where groundwater levels showcase an overall increase (though very limited, only in the range of less than 1 m), contrary to what is observed in the deformation time series, which depicts a negative cumulative deformation of 31 mm.

Descriptive statistical analysis enabled the identification of any possible correlation between the observed deformation velocity and the geological framework (Figure 12c). Key statistics (mean, median, mode, standard deviation, maximum, and minimum) for the velocity values of the points located within four different textural categories reported in the geological mapping of Ravenna were calculated. Areas encompassed by the lithified dune and the silty sands display similar statistics, with a mean of −3.8 mm/year that could be associated with the natural subsidence rates experimented by the region. The silty clays and sands, on the other hand, highlight a lower mean, median and mode, and a higher standard deviation (mean of −5.3 and −6.2 mm/year, respectively; standard deviation of 2.3 and 1.8 mm/year). When analyzing these statistics, it is worth noting that the density of PS points across the four lithologies varied, due to the different distribution of urban structures and reflective radar targets across the landscape of the study site. The silty clay lithology recorded sparser densities (855.7 PS/km²) with respect to the silty sand, sand, and lithified dune (1532.8, 1374.4, and 1319.1 PS/km², respectively); therefore, this might influence the confidence levels of the extracted statistics for the four lithologies. Sands are distributed mainly along the coast and present the lowest average, with a minimum of −24 mm/year, aligning with the coastal subsidence rates observed in Figure 2a. Localized peaks in the observed displacements could link with anthropogenic drivers, such as land use change. For instance, the statistics reveal the highest standard deviation and the most extreme subsidence values for the silty-clayey terrain. Considering that the industrial area of Ravenna and the dump site are located onto this geological unit (see Figure 12a), the high velocity could result from compaction by the loading of the buildings in the clayey sediments layer. Moreover, Urban Atlas Land Cover maps show an expansion of the dump site during the period 2012−2018 [56], modifying a previous forestry land cover and, therefore, probably triggering the observed subsidence process.

4.3.2. Bologna

When analyzing the change in piezometric levels (∆h_i) recorded in 2018−2022 across the area of Bologna (Figure 13a), it is possible to notice a deepening of groundwater levels over the central-western part of the city and in the area of San Giovanni in Persiceto. Most of the wells are located over the Reno River alluvial fan, known for its historical pumping record [59,60]. On the other hand, a recovery in the piezometric levels is observed in the northern, eastern, and western sectors of the study area. Even accounting for the fact that the different piezometers operate at different depths and, as such, at different layers of the multi-aquifer system of Bologna, the spatial pattern of changes in piezometric levels resembles the observed subsidence velocity distribution across the area (Figure 3a), suggesting the correlation between groundwater extraction activities and subsidence.

Comparison of piezometric level variations with deformation time series at neighboring PS–DS points in two industrial and one monitoring well (BO05-00, BOF6-00, and BOF9-00, respectively), showcase the local relationship between groundwater level changes and ground deformation (Figure 13b), initially investigated in [61]. During 2018−2022, a deepening of −10.0 m in piezometric levels was recorded in well BOF6-00 (drilled within the Reno-Lavino cone, reaching the unconfined aquifer, with filter at 90−131 m), and a cumulative vertical ground deformation of −17.8 mm was estimated in the area. The seasonal oscillation observed in both time series, with peak values during the winter–spring season and lower values during summer–autumn, suggests a close relationship between piezometric level variations and the observed deformation. On the other hand, although well BOF9-00 (located in the Apennine alluvial plain, reaching the upper confined aquifer, with filter at 100−130 m) presents a sudden recovery of the piezometric levels in the last measurement (recorded on 21 November 2022), leading to a positive Δh_i of +4.4 m, no remarkable change is observed in the deformation time series, which experiences a cumulative vertical displacement of −63.6 mm. Finally, well BO05-00 (within the alluvial plain, reaching the lower confined aquifer, with filter at 168−172 m) experiences a Δh_i = −14.7 m, which corresponds to −16.2 mm vertical displacement.

The study area’s morphology and geology are primarily influenced by the alluvial fans of the Reno and Savena rivers (Figure 13c). These fans consist of layered coarse-grained deposits that transition to finer, less amalgamated sediments in their distal areas [62]. Following the workflow applied for the analysis of Ravenna, a descriptive statistical analysis was performed to assess the correlation between deformation and different textural classes (Figure 13d). The geological units were retained from Bologna geological mapping at scale 1:25,000 [46], and the density of PS for each lithology was: 0.3 PS/km² for clays, 901.6 PS/km² for sands, 1573.2 PS/km² for gravels, and 1057 PS/km² for silt. Similar to what was observed for Ravenna, clay and sand lithologies exhibit lower subsidence velocities, with a mean of −7.1 and −6.5 mm/year, respectively. Silt and gravel lithologies have their mean among the natural subsidence values of the Emilia-Romagna region (up to −3 mm/year).

4.3.3. Carpi, Correggio, and Soliera

The industrial area of Carpi, Correggio, and Soliera encompasses the north-eastern part of Reggio Emilia province and the north-western sector of Modena province (Figure 14a). In order to evaluate the relationship between groundwater extraction activities and the ground deformation observed during the studied time period (2018−2022), a comparison between change in the piezometric levels, recorded in MO10-01 monitoring well (drilled within the Apennine alluvial plain, reaching the upper confined aquifer at 120 m depth), which is used for irrigation purposes, and cumulative ground deformation in a contiguous PS–DS point was pursued (Figure 14b). A decrease of −2.64 m in the piezometric levels is observed, with an associated cumulative ground deformation of −82 mm. Deformation time series display a slightly periodic component, which aligns with seasonality observed in piezometry variations from the second half of 2020 but displays an apparently opposite tendency in the previous period (2018−2019 to mid-2020).

Similar to the other sites analyzed, the area is characterized by the presence of highly compressible materials, with clays and silts being the only textural classes present among the lithologies (Figure 14a). The statistical correlation analysis between mean deformation velocities and the different geological layers showed negative values in all the cases, below the ones considered as natural subsidence rates (<−3 mm/year) (Figure 14c). While the sandy silt lithology shows the lowest mean subsidence velocity (−7.7 mm/year), clay and sandy silt display higher deformation rates, with mean velocities of −4.4 and −7.0 mm/year, respectively. The PS density is higher for the sandy silt (1524.3 PS/km²) and presents sparser values for the silt and clay lithologies (941.8 and 480.6 PS/km², respectively).

5. Discussion

5.1. PS-Time and ICA for Deformation Time Series Analysis

EGMS datasets offer an unprecedented temporal sampling, with calibrated displacement records available up to the weekly scale, using both acquisition geometries, and derived vertical and east-west datasets. This characteristic clearly opens new perspectives for the detailed assessment and characterization of land subsidence phenomena, as well as many other geohazard processes, the evolution and dynamics of which can be characterized using DInSAR data. Yet, the notable potential of EGMS data, on the other hand, calls for semi-automated and statistical approaches capable of supporting their advanced analysis, to ensure a detailed investigation of displacement trends and patterns of interest.

The integration of advanced statistical methods such as PS-Time and ICA into DInSAR data analysis for the three study sites in this work has significantly enhanced the identification of different patterns of motion in the deformation time series. The application of PS-Time proved to be effective for classifying the various deformation trends, especially for distinguishing clusters of points that experiment quadratic behaviors. In all cases, quadratic points were clustered in the areas where EGMS products displayed accelerating or decelerating behaviors, indicating a strong correlation between this category of the classification and a change in deformation velocities during the time series (Figure 2b–d, Figure 3b–d, and Figure 4b–d). Additionally, the tool has succeeded in identifying the existence of a breakpoint in the time series for a specific period of time, displaying a cluster of bilinear points over a small area of Soliera (Figure 4d). This functionality is very interesting, as not only enhances the identification of strongly focalized deformation hotspots over the study area, but also provides the specific time frame of when the change occurred. However, the classification provided by the program is highly dependent on the specific parameters set for the statistical analysis [32]. Therefore, it is highly recommended to run numerous trials with subsets of the input datasets, to correctly define the most appropriate thresholds for the analysis. Furthermore, as the PS-Time platform analyzes each point individually, executing all the statistical analysis in the series, the computational demands when processing large datasets are significant, leading to extended processing times. Given the point-by-point analysis approach, however, the computational load could easily be distributed by splitting the input dataset into subsets and using multiple processing machines, or via parallelization of the code and the exploitation of high-performance computing infrastructure.

On the other hand, ICA proved to be effective for separating the different motion trends embedded in the deformation mixed signal, successfully isolating linear, seasonal, accelerating, and decelerating patterns with a high degree of precision (Figure 5, Figure 6 and Figure 7). As can be appreciated from the correlation analysis shown in Figure 8, Figure 9, Figure 10 and Figure 11, a great correspondence between the retained components and the different motion trends identified was observed. Independent components showing linear or seasonal trends were spatially distributed in concordance with mean deformation velocity and seasonality maps, respectively, while the components showing a quadratic behavior spatially aligned with accelerating and quadratic PS-Time clusters. However, the analysis performed by the ICA algorithm also presented some limitations. Even though ICA is considered a blind signal separation method, when applying temporal ICA for evaluation of time series, some a priori knowledge is required, in order to properly define the number of components to be retained from the mixed signal. In this study, the number of PCs selected for each study area (which consequently correspond to the number of ICs retained by the program) was evaluated with a trial-and-error approach, taking into consideration the results obtained from PS-Time and EGMS products. Although the components managed to identify the different trends of motion, a trade-off had to be made in order not to retain too many noisy components, which would have made the interpretation of the results more challenging. Consequently, some components still showed mixed signals (e.g., slightly seasonal signals embedded within quadratic components in Bo1-IC1 and in So-IC1), although one of them is remarkably more significant than the other. Additionally, ICA successfully detected the accelerating trends (e.g., So-IC1), but did not identify any clear breakpoints that were, on the other hand, undoubtedly highlighted by PS-Time. Moreover, as the components could display mirrored deformation behaviors (e.g., Bo2-IC3), a priori knowledge is advisable to avoid any possible misinterpretation of the results.

Some of the limitations of PS-Time were therefore overcome by adding ICA to the analysis workflow. In the areas where quadratic clusters were identified by PS-Time, ICA complemented the analysis by differentiating between accelerating and decelerating behaviors, scoring the points positively or negatively in comparison to the retained quadratic signal. This is the case, for example, for Area 2 in Bologna, where the decelerating area of San Giovanni in Persiceto scored negatively for an accelerating quadratic trend of Bo2-IC3. Additionally, ICA has successfully identified another seasonal, but noisy, cluster of points for the area of Anzola dell’Emilia, which was not noticeable in the previous seasonality analysis. Hence, the combination of these two methodologies has the potential to fully identify and interpret the different components of motion that compose the mixed deformation signal.

5.2. Subsidence in Ravenna, Bologna and Carpi–Correggio–Soliera

Apart from the −3 to −5 mm/year natural subsidence rate affecting the whole Po Plain [40,42], significant anthropogenically-driven subsidence has been recorded through the years by many authors for the three land subsidence prone areas assessed in this study.

Previous investigations reported deformation velocities between −10 and −20 mm/year in Ravenna over the coastal area and related them with the combined effects of groundwater extraction and methane withdrawal from both onshore and offshore gas reservoirs, which leads to pressure loss and subsequent compaction [38,39,40,41,48,57]. The most recent report published by the ARPAE subsidence monitoring project [33], which investigated the evolution of subsidence over the Emilia-Romagna region between 2016 and 2021, recognized coastal deformation velocities in Ravenna between −7.5 and −10 mm/year, aligning with what was observed in Figure 2a. Even though subsidence resulting from groundwater exploitation has been recorded since the early 1950s [39] and different mitigation plans have been implemented through time, involving the decrease of pumping rates, recent studies prove that the recorded subsidence values along the coastline show no significant correlation with water withdrawal activities, consequently recognizing methane extraction as the main driver [48]. In accordance with what was observed by these authors, piezometric level variations in well RA29-00 seemed to have no correlation with the contiguous subsidence values recorded (Figure 12b). However, the opposite situation was observed for well RA49-00, where ground deformation showed some correlation with piezometric level variations. This could explain subsidence in the accelerating agricultural fields area, which has already been recognized in the past by [35] as the ground motion area of Savarna and Mezzano, groundwater withdrawal being the most statistically significant driver identified.

Recognized as one of the most important subsidence hotspots over the Emilia-Romagna region, the city of Bologna experienced peak subsiding velocities of up to −100/−110 mm/year in 1970−1982, and a lessening in ground deformation rates has been observed during the last 30 years, mainly due to an increase in regional groundwater pumping regulations and subsidence monitoring activities [59]. Although a significant recovery of deformation rates was observed in 2011−2016, with areas that even showed uplifting behaviors [63], a decrease of up to −20 mm/year in subsidence velocities has been observed in 2016−2021, according to the most recent ground velocity change map published by ARPAE [33]. Subsidence of up to −27.7 mm/year was recognized in the area north of Bologna for the 2016−2021 period [33], aligning with the subsiding values depicted in Figure 3a, that range between −20 and −30 mm/year for the same place, with peak velocities of −35.8 mm/year. The correlation between the subsidence phenomena and groundwater withdrawal activities in Bologna has been observed by different authors in the literature [38,59,60], even focusing on the same study period [64]. The comparison between the change in piezometric levels and the equivalent subsidence velocity implemented in this study showed a good correlation in the majority of cases (except for well BOF9-00, Figure 13b), suggesting that pumping activities continue to be one of the main drivers for the subsidence phenomena affecting the broad area. However, the geological setting also conditions the occurrence of ground deformation, as historically the highest deformation velocities in the area occur mainly associated with the Reno river alluvial fan [38,59,60,62], which is also the place where most of the pumping activities of the area are concentrated. The results of the correlation between the recorded subsidence velocities and the different lithological classes (Figure 13d) suggest that sediment compaction is more predominant in clayey materials, an observation that coincides with what was observed by other authors [59]. Nevertheless, this observation is based on the sole analysis of statistics across the geological layers, with no account for the spatially variable thickness of geological deposits or their geotechnical properties; therefore, more investigation should be done on the detailed correlation with these data to improve this initial assessment. Moreover, a detailed contextualization of the results with respect to the hydrogeology of the area would also be beneficial to complement the investigation.

Finally, subsidence over the industrial area of Carpi, Correggio, and Soliera was also recognized by [65], who recorded deformation velocities of up to −15.0 mm/year over the area for the time span between 2012 and 2022. Moreover, ARPAE’s last report depicts subsidence values of −27.3 mm/year, −15.0 mm/year, and −14.9 mm/year for each town, respectively, for 2016−2021. The values identified by the literature are in accordance with the ones identified in this study (Figure 4a), that present subsidence values for the area that vary between −10.0 and −20.0 mm/year, with peak velocities of up to −25.2 mm/year in the north-western sector of Carpi. Although subsidence has diminished for Correggio and Soliera in the last 5 years, Carpi continues to experience lower deformation velocities [33,63] The presence of this phenomenon has been attributed to groundwater withdrawal activities together with overload due to industrial expansion, in an area characterized by the presence of highly compressible materials [65]. Correlation analysis carried out in Figure 14b,c support this hypothesis, as a deepening in piezometric levels is accompanied by ground deformation in the surrounding area, and all the mapped lithologies in the area, apart from being highly compressible (clays and silts), show values below the natural subsidence rate.

6. Conclusions

Thoroughly characterizing the evolution of subsidence deformation in a certain period is important to recognize its main drivers, their identification being crucial for the implementation of appropriate mitigation measures. By employing different statistical tools, it is possible to evaluate acceleration or deceleration in the deformation time series, as well as cyclic or seasonal trends. In this study, the evolution of the subsidence phenomena between 2018 and 2022 in three land subsidence prone areas of the Po Plain was accurately characterized through the application of a three-step workflow, which includes the following: (1) the application of the semi-automated algorithm “PS-Time” [11], (2) an independent component analysis (ICA), and (3) a correlation with groundwater and geological data. The input data included EGMS L3 and L2b products, with the latter being particularly important because, in addition to providing the complete deformation time series of each PS–DS point, they calculate the acceleration values at each single target, serving as validation for the following statistical analysis.

The PS-Time tool [11] was applied in the three selected study areas (the coastal area of Ravenna, a broad area comprising Bologna, and the industrial area of Carpi, Correggio, and Soliera) in order to classify the time series of the points according to their characteristics and to identify any seasonal trends. Mean deformation velocity, acceleration, seasonality, and trend classification maps were generated for each study site. These maps not only allowed for the identification of the places affected by subsidence, including the recognition of local hotspots with the lowest deformation velocities in each area (e.g., the dump site at Ravenna with deformation rates of −37.1 mm/year, or an area in the northern part of Bologna with up to −35.8 mm/year) but also enhanced the recognition of ground motion patterns that differ from the classical (and often blindly assumed) linear trend. For example, a cluster of deformation trends with significant seasonality was identified in localized areas of Ravenna and Bologna (with peak amplitudes of 6 mm and 9 mm, respectively), and accelerating deformation time series were identified for an area with agricultural fields in Ravenna, for the industrial area of Ponte Samoggia and around Ozzano dell’Emilia in the area of Bologna, and over a strongly localized area of Soliera. Conversely, a small coastal area to the south of Ravenna, the city of Bologna, San Giovanni in Persiceto, and the industrial area of Carpi depicted decelerating time series. Moreover, trend classification maps showed clusters of quadratic points in the areas where accelerating/decelerating time series were identified, with even a bilinear cluster of points in a strongly accelerating area in Carpi–Correggio–Soliera, emphasizing the potential of this classification tool for identifying nonlinear time series.

Areas showing interesting deformation patterns were retained to perform an ICA. While a single window area of ~22 km² was carefully selected for each study site, an additional ~52 km² testing sector was retained for Bologna, encompassing the areas of San Giovanni in Persiceto, Ponte Samoggia, and Anzola dell’Emilia (named “Area 2”). In all cases, the different deformation trends embedded in the time series were thoroughly retained by using four components, except from “Area 1” of Bologna, which was completely described by using only two components. A linear, a quadratic, and one or two seasonal components were retained for each area. Moreover, a noisy component was recovered for all cases, either clearly represented as noise (e.g., So-IC3) or as a noisy sinusoidal trend (e.g., Bo2-IC4). By spatially comparing the distribution of the components, it was noticeable that components displaying linear trends correlate with the deformation velocity clusters previously identified, quadratic components with accelerating/decelerating trends and sinusoidal components with seasonal time series, validating the previously obtained results and discovering new clusters (e.g., a new cluster of seasonal time series in Bologna). This appreciation was confirmed by correlation analysis, which showed R² between 0.4 and 0.9.

Finally, the results were compared with piezometric levels recorded by ARPAE and with geological layers downloaded from MinERva Portal [46]. Generally, in places where high subsidence values were identified, a decrease in the piezometric levels was also recorded, underscoring the role of aquifer recharge and withdrawal on land subsidence. On the other hand, the correlation with geology revealed that lower deformation velocities are mainly associated with the most compressible materials (e.g., silts and clays).

The conjunct use of PS-Time and ICA for time series evaluation that has been experimented in this work, as well as the implementation of spatio-temporal correlation analysis for drivers’ evaluation at the three subsiding sites, has proved to be a powerful methodology to assess the evolution of the subsidence phenomena over a specific time period, with a particular strength on identifying deformation patterns that vary from linear trends. Despite the discussed limitations of each methodology itself, their conjunct use provides truthworthy identification and full characterization of linear, accelerating, decelerating, and seasonal deformation trends, as well as their spatial distribution. The evolution of ground motion in an area over time is closely related to the drivers that derive the subsidence phenomena, and future research may therefore focus on the more extensive implementation of such methods across wider areas, for instance across the whole territory of the Emilia-Romagna region to enable regional authorities and end-users to easily identify anthropogenic activities affecting the area and, therefore, inform land and groundwater resource management.