Sentinel-1 P-SBAS data for the update of the state of activity of national landslide inventory maps

The redaction of landslide inventory is a fundamental task for risk management and territorial planning activities. The availability of synthetic aperture radar imagery, especially after the launch of Sentinel-1 mission, enables to systematically update landslide inventories covering wide areas in a reduced time frame and at different scales of analysis. In this work, SAR data processed from the fully automatic P-SBAS pipeline have been adopted to update the Italian national landslide database. Specifically, a matrix has been introduced by comparing past landslide state of activity obtained with Envisat data (2003–2010) and that computed with Sentinel-1 (2014–2018). The state of activity was defined by obtaining the projected velocity along the slope dip direction. The analysis involved about 56,000 landslides which showed at least one Sentinel-1 measurement point, of which 74% were classified as dormant, having annual average velocity < 7 mm/year (considering a value of two times the standard deviation) and 26% as active (mean velocity > 7 mm/year). Furthermore, a landslide reliability matrix was introduced on the landslide inventory updated with S1 data, using the measurement point (MP) density within each landslide and the standard deviation of the mean Vslope value of each landslide. In this case, the analysis revealed that more than 80% of landslides has values of reliability from average to very high. Finally, the 2D horizontal and vertical components were computed to characterize magnitude and direction of every type of landslides included in this work, showing that spreadings, deep-seated gravitation slope deformations, and slow flows showed a main horizontal movement, while complex and translational/rotational slides had more heterogeneity in terms of deformation direction. Hence, the work demonstrated that the application of fast and automatically nationwide Sentinel-1 MTInSAR (multi-temporal interferometry SAR) may provide a fundamental aid for landslide inventory update.


Introduction
Slope instabilities represent globally one of the most remarkable and widespread natural hazards, determining a considerable number of casualties and huge economic losses (Schuster 1996). Italy is one of the countries affected the most by landslides: in the last 5 years, 26 victims and 107 injuries were reported, while landslide events involved 19 regions (out of 20) and 365 municipalities (IRPI 2021). The vulnerability of Italy to landslides may be addressed by several factors. Its geological and geomorphological setting makes the territory very prone to such phenomena; in addition, an uncontrolled urban sprawl following the economic growth of the early 1960s influenced and still influences the vulnerability and exposure of physical and social elements with respect to different landslide hazard severity. In order to assess landslide-prone areas, trustworthy and upgraded landslide inventories constitute the chief element on which to base hazard and risk analysis, landslide localization, extent, distribution, and typology (Brabb 1991). Landslides leave discernible signs, most of which can be recognized, classified, and mapped (Guzzetti et al. 2012). Their detection and interpretation can be achieved in multiple ways, from in situ surveys, or through the implementation of remote sensing method such as visual interpretation of stereoscopic aerial photography (Nichol et al. 2006;Ardizzone et al. 2012;Holbling et al. 2016;Li et al. 2016;Del Soldato et al. 2018a, b), visual analysis of digital elevation models (DEMs) and derived products (Conforti et al. 2014;Lazzari et al. 2018;Pawluszek 2019), semi-automatic recognition, and object-based algorithm applied on high-resolution DEMs (Bunn et al. 2019;Pradhan et al. 2020) or on optical images (Martha et al. 2011). The use of satellite imagery has become a standard way of detecting and mapping landslides. Among the satellite technologies, the synthetic aperture radar (SAR) has been demonstrating to be an ideal solution for the detection and mapping of landslides, either using polarimetric techniques (Plank et al. 2016;Park and Lee 2019), coherence (Konishi and Suga 2019;Tzouvaras et al. 2020;Burrows et al. 2020), or amplitude of pre-and post-failure images (Raspini et al. 2015;Tessari et al. 2017;Mondini et al. 2019). More details on the exploitation of SAR imagery for landslide detection can be found in the review by Mondini et al. (2021). All the abovementioned techniques are capable of detecting the areas affected by landslides, however without providing (in many cases) accurate measurement of the ground surface displacements and without the possibility of reconstructing the displacement time series. In this sense, the adoption of multi-temporal differential interferometry SAR techniques (MTInSAR) can be advantageous not only for the static mapping of landslide phenomena but also to update geomorphological and multi-temporal landslide inventory maps (Guzzetti et al. 2012;Solari et al. 2020), determining the landslide state of activity (SoA) and its temporal evolution. Several MTInSAR applications for landslide mapping can be found in Italy, each of them with different approaches, data and technique implemented, and scale of analysis, from basin (Cascini et al. 2013 andDel Ventisette et al. 2014), regional (Rosi et al. 2018;Guerriero et al. 2019) and national scale (Not-Ordinary Plan of Remote Landslides (2023) 20:10831097-Sensing, Di Martire et al. 2017Costantini et al. 2017). Di Martire et al. (2016) developed an integrated system, based on mixing ground surveys and persistent scatterers (PS) InSAR data derived from COSMO-SkyMed imagery, to detect and update landslide inventory in Palermo province (Southern Italy). COSMO-SkyMed data were also implemented by Antonielli et al. (2019) to update the state of activity of landslide in 26 areas of Lombardia (Northern Italy), as well as by Bonì et al. (2018), in combination with ERS and Radarsat imagery, to update landslide inventories of an area of Piemonte Region (NW Italy). Being adopted for several years with standardized procedures, updating landslide inventory maps is now a consolidated practice. Despite this, the scale of analysis still ranges from basin to regional scale. In this work, the potential of MTInSAR approaches in wide-area mapping has been fully exploited for the first time at national scale.
Here, the C-Band Sentinel-1 (S1) data set acquired on the entire Italian territory up to December 2018 to update the national IFFI landslide inventory (Inventario Fenomeni Franosi in Italia, in Italian, Landslide Inventory in Italy, Trigila et al. 2010). In particular, the S1 data have been processed through the Parallel Small BAseline Subset (P-SBAS) InSAR algorithm (Berardino et al., 2002;Casu et al. 2014;Zinno et al. 2015;Manunta et al. 2019) in the framework of an Operative Agreement with the Italian Ministry of Economic Development (MiSE). To update the IFFI inventory, a matrix approach has been adopted by cross-comparing the landslide state of activity (SoA) that was defined through Envisat (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) data with the SoA computed with Sentinel-1 (2014-2018) data. To this, the mean velocity projected along the slope value (V slope ) was used. Finally, P-SBAS data were also decomposed in 2D to obtain the vertical and the horizontal E-W component, to statistically analyze the main components for each type of landslide considered.

IFFI landslide inventory
The IFFI project was launched in 1999 with the aim of identifying and mapping landslides throughout Italy on the basis of standardized criteria (Trigila et al. 2010). It is redacted by the ISPRA (Istituto Superiore per la Protezione Ambientale, in Italian, Superior Institute for the Environmental Protection) in collaboration with regional authorities and autonomous provinces. The IFFI project started in the aftermath of the Sarno and Quindici landslides of May 5, 1998, which killed 161 people. ISPRA has been collecting online data since 2005 with the objective of enhancing the dissemination and fruition of useful information to local administration, research groups, and technical staff in charge of design and urban planning (Trigila et al. 2010). All the data are publicly accessible through the IdroGEO platform (https:// idrog eo. ispra mbien te. it/ app/ iffi?@= 41. 55172 52589 4153,12. 57350 14838 1829,1). The classification of the landslides is based on schemes by Varnes (1978), Cruden and Varnes (1996), and recommendations by the International Association of Engineering Geology (1990), the International Geotechnical Societies UNESCO Working Party on Word Landslide Inventory (1990, 1991, 1993a, b, 1994). The IFFI database contains vector layers of landslides and main information such as landslide location and type of movement. For the latter aspect, 11 types of movement are distinguished: (i) fall/toppling, (ii) rotational/translational slide, (iii) spreading, (iv) slow flow, (v) rapid flow, (vi) sinkholes, (vii) complex, (viii) area with diffuse falls/topple, (ix) area with diffuse sinkholes, (x) area with diffuse shallow landslides, (xi) deep-seated gravitational slope deformation (DSGSD), (xii) types of landslides which cannot be determined (ND). The distribution of the type of movement divided between the North of Italy (Fig. 1a), Central Italy (Fig. 1b) and South of Italy (Fig. 1c)  It must be specified that for the analysis conducted and presented in this paper, only landslides characterized by a slow-and intermittent-moving behavior, namely rotational/translational slide, spreading, slow flow, complex, area with diffuse shallow landslides, DSGSDs, and ND, are considered due to the intrinsic limitations of interferometric products. Furthermore, not all the landslides are identified as polygons, being not always mappable at a 1:10,000 scale. Firstly falls/topple, rapid flow, sinkhole, area with diffuse falls/topple, and area with diffuse sinkholes were excluded from the analysis; for the remaining landslides, only polygons were used, for a total of 458,545 features ( Fig. 1). In Fig. 1, landslides selected according to their kinematic are defined as MTInSAR IFFI landslides, while non-MTInSAR landslides are all those landslides discarded at the beginning of the analysis.

Interferometric dataset used
Two different interferometric datasets have been used in the present work: Envisat data, covering the timespan 2003-2010, and Sentinel-1 data acquired from 2014 to the end of 2018.
The first one belongs to the Not-Ordinary Plan of Environmental Remote Sensing (PST-A in Italian), a national project issued by the Ministry of Environment and Protection of Land and Sea with the aim of constituting a national database of active or potential instability phenomena affecting the Italian territory, based on the exploitation of interferometric product (http:// www. pcn. minam biente. it/ mattm/ proge tto-piano-strao rdina rio-di-teler ileva mento/ and Di Martire et al. 2017). In Fig. 2, the distribution of the PSs in ascending and descending modes is visible. ENVISAT data were processed by means of PS-InSAR (Ferretti et al. 2001) and PSP-DIFfSAR (Constantini et al. 2009).
The second dataset includes the output of a multi-temporal analysis on images acquired by Sentinel-1 (S1), a constellation made of S1A satellite, launched on 3 April 2014, and its twin, S1B, launched on 25 April 2016, sharing the same orbital plane, promoted by the EU Copernicus program and managed by ESA (European Space Agency). The abovementioned multi-temporal analysis has been generated within the operative agreement between the IREA-CNR and Italian Ministry of the Economic Development (MiSE, https:// www. cnr. it/ it/ proge tti-di-ricer ca/ proge tto/ 17105/ accor do-opera tivo-mise-dgrme-e-cnr-irea-dit-ad012-028), which is aimed at generating the ground displacement time series of the entire Italian territory by exploiting the S1 constellation data. Differently from the Envisat data, S1 data were processed through the application of an efficient SBAS approach, namely P-SBAS, which allows retrieving ground displacement information also from distributed scattering (DS) areas (Zinno et al. 2018;Manunta et al. 2019;De Luca et al. 2017;2019). The P-SBAS processing has been carried out in a fully automated way, without any tuning tailored to the specific phenomena under investigation. S1 data, due to the very short orbital tube of the satellite platforms (mean diameter is below 400 m), are very suitable for the SBAS approach, relying on data characterized by minimized spatial and temporal baseline.
The P-SBAS technique is also capable of using distributed computing infrastructures, making this technique extremely suitable for fast processing of huge volumes of S1 data (for further details on P-SBAS, please refer to Manunta et al. 2019 andDe Luca et al. 2022). In green, types of landslides and relative numbers selected for the analysis, in red those discarded

Fig. 2
Envisat data over the Italian territory. On the left, descending dataset, on the right, ascending one. The purple lines indicate the regional borders. The coverage for ENVISAT imagery is for the total Italian peninsula (excluding thus Sardegna) in ascending mode, while in descending mode also part of Calabria, Puglia, Basilicata, and Campania (and Sardegna, as well) were excluded S1 data here used include six descending and five ascending tracks, divided into 19 and 17 frames, respectively, which correspond to approximately 17,000 S1 slices and more than 6000 acquisitions. The retrieved MTInSAR products have been obtained with a resulting pixel dimension of about 80 × 80 m, and ground pixels with a temporal coherence value greater than 0.9 were considered. Taking into account both ascending and descending MTInSAR nationwide datasets, a total of 15,587,069 MPs have been made available for this study (Fig. 3).

Landslide activity and reliability matrices
To update the IFFI landslide inventory over the national territory of Italy, mean velocities of deformation, derived from MTInSAR processing, were adopted. First, ENVISAT data were used to assess the landslide SoA in the timespan 2002-2010; then, S1 P-SBAS data were implemented to assign a SoA to each landslide, covering the time 2014-2018. In both cases, all the landslides with at least one Measure Point within their boundary were considered in the final computation. To discriminate between active and inactive landslides, the representative velocity of each landslide polygon had to be defined. For the IFFI update, the projected velocity along the slope (V slope ) was adopted. This assumes a greater significance when dealing with movements along the slopes since the use of V Los values enable the measurement of a percentage of the landslide movement, depending on the combination of slopes morphology (i.e., aspect and slope gradient) and SAR acquisition geometries (Plank et al. 2012). In detail, to obtain the V slope value, the computation of the so-called C coefficient (Notti et al. 2014) is carried out through the formula: where N, E, and H are LOS directional cosines and change according to the satellite data and to the geometry of the acquired scene, while slope and aspect angles derived from the 10 × 10 m DEM (digital elevation model) of Italy by Tarquini et al. (2012). The C coefficient has been limited to − 0.2 if − 0.2 < C < 0 and to 0.2 if 0 < C < 0.2 in order to not exaggerate the projection, as suggested by Herrera et al. (2013) and Bonì et al. (2018). Once the C coefficient for each MP is obtained, the V slope is computed through the following formula: To discriminate between stable (i.e., not showing any sign of displacement as retrievable from InSAR data) and unstable landslides, a mean velocity threshold has been chosen by using the double value of the MPs standard deviation (σ) of the mean velocity, which is equal to 3.5 mm/year. Therefore, landslides showing a mean velocity over 7 mm/year are considered unstable. For landslides with just one MP, the same value of deformation mean velocity was chosen.
(1) Fig. 3 Sentinel-1 data covering the Italian territory. On the left, descending dataset, on the right, ascending one. The purple lines indicate the regional borders To estimate the state of activity and to recognize the temporal evolution of each landslide, an "activity matrix," as implemented by Righini et al. (2011) and Del Ventisette et al. (2014), has been adopted (Fig. 4). The activity matrix consists of a grid in which nine different velocity combinations are identified, by crossing the values of mean velocity obtained with ENVISAT (2002ENVISAT ( -2010 and Sentinel-1 (2014Sentinel-1 ( -2018 data. Landslides are hence quantitively identified as follows: 1. Stable, in case both ENVISAT and S1 identify landslides with velocity values below the threshold (hereon defined as SoA1); 2. Dormant, in case ENVISAT data are showing active landslides which are instead stable with S1 data (SoA2); 3. Dormant, when landslides have no ENVISAT data but just S1 (below the threshold) (SoA3); 4. Reactivated, when landslides previously stable (with ENVISAT) are showing displacement over the threshold with S1 (SoA4); 5. Active continuous, when either ENVISAT and S1 data show values greater than the thresholds (SoA5); 6. Reactivated/active continuous, since no ENVISAT data are available, and only S1 data can define the active state (SoA6); 7. Stable/reactivated, since no S1 data are available, the state of activity can be defined up to 2010, with MPs under the threshold (SoA7); 8. Dormant/active continuous, since no S1 data are available, and landslides were active with ENVISAT data (SoA8).
A 9th class, related to all those landslides without any MP from both datasets, has been excluded from the final computation (ND, not determined).
Many limitations have to be considered when merging MTIn-SAR data and IFFI landslides. In detail, mean velocity variability can be observed among different MP within the landslide boundaries, which can be either due to different behavior of landslide sectors, to the different exposure of the slope or to errors of the phase unwrapping. Also, the MP distribution, in terms of number of targets within each landslide, can vary, depending on the size of the landslide or on the visibility of the slope, in terms of temporal coherence or geometric decorrelation.
A "reliability matrix" was generated to assess the trustworthiness of the landslides whose SoA has been defined using the S1 data. The MP density was considered, expressed as the number of MP per square kilometers, along with the standard deviation (DevSt) of the mean velocity of each landslide (Fig. 5). Nine classes were defined, by using half of the average and the average of both values to discriminate the classes. As regards the MP density, the average value is 151 MP/km 2 , and 150 and 75 were set to discriminate the classes, while regarding the DevSt of the mean velocity, with an average value of 1.14 cm/year; therefore, 1.2 and 0.6 cm/yr are the thresholds identified from the lowest to the highest.
The classes, from the lowest to the highest reliability, are defined as follows:

LOS velocity decomposition
By definition, MTInSAR measurements are generated in the LOS geometry only, thus providing a one-dimension measurement of the ground motion. Many works have shown the capability to project the LOS measurements along the vertical and the horizontal (East-West) component, by using both ascending and descending analyses (Casu and Manconi, 2016;Del Soldato et al. 2018a  Considering that the orbital path of the satellite is approximately parallel to the meridian, the satellite SAR technology can be considered almost blind to the motion in the N-S direction (Tamburini et al. 2010); therefore, the vertical velocity (V v ) and the horizontal one (along the E-W, V E ) are calculated according to Notti et al. (2014).
Values of V v and V E have been derived on a regular grid with a cell resolution of 100 × 100 m, where common targets, in ascending and descending mode, are selected, thus generating a synthetic MP with both values.

Results
Landslide state of activity and activity matrix Satellite SAR systems are not able to detect movements along the N-S direction, as well as geometric distortions can affect the measurements. Moreover, low coherent areas are discarded within each MTInSAR processing chain, through the implementation of a coherence threshold (Manunta et al. 2019). For these reasons, the final computation of the state of activity has regarded 46,217 and 56,133 landslides for the Envisat and the S1 datasets, respectively, showing at least 1 MP. To determine the SoA of each landslide, a stability threshold has been applied, as identified previously equal to 7 mm/year. The representative velocity of each landslide has been obtained by averaging the value of the velocity of all the MP within the landslide boundary. In this way, about 37,000 were classified as dormant and 8000 as active, according to the Envisat data, while for S1, about 41,000 landslides show a dormant state of activity, and 14,000 an active one (Fig. 7a). The database of landslides used in this work and correlated with the Sentinel-1 data to update their state of activity is accessible and downloadable in the Supplementary Information section. Slow flows and areas with shallow landslides show the highest percentage of active landslides, among the various types of landslides detected with S1.
As seen in Table 1, about 11% of the total landslide polygons have been detected with S1 data. On a regional basis (and it is resumed in Table 2), the region showing the highest number of landslides detected is Emilia Romagna with 7507 phenomena, while the lowest number is shown by Puglia Region (99 landslides). This is obvious considering the areal extension of the region and the Landslide Index, which sees Emilia Romagna as having one of the highest values. Normalizing the value of landslides detected according to the number of total landslides per region is it possible to highlight that in Calabria more than 50% of landslides were detected by S1 data, followed by Basilicata and Abruzzo (approximately 35% and 25%, respectively), while the lowest ratio between landslide detected and a total number of events is observed in Friuli Venezia Giulia and Toscana Regions, with a value of about 4%. Considering the sum of the areal extent of each landslide considered, the region with the largest area affected by landslides detected by S1 data is still Emilia Romagna (for the abovementioned reasons), while the smallest one is the Sardegna Region (890 vs. 39 km 2 ). Taking into account the ratio between the total landslide area covered by MP and the total landslide area in each region, about 31% of landslide areas in Italy contain S1 MPs. At a regional scale, the highest value is seen in Calabria Region (80%) while the lowest in Friuli Venezia Giulia Region (7.6%).
According to the landslide activity matrix shown in the "Estimation of the time-varying landslide stability" section, the SoA of each polygon has been determined and classified as follows: -18,345 landslides were classified as SoA1; -3026 landslides belong to the class SoA2; -for 20,236 landslides, class is SoA3.
The resume of the above-listed figures is mapped in Fig. 6 while a bar chart summarizes the above-listed figures (Fig. 7c).
Summarizing the main figures, 19,757 landslides did not show any change in the state of activity, keeping dormant or active from 2001 to 2018, 6635 landslides changed their state, passing from dormant to active or in the other way round during the years and for 49,566 of them, it was not possible to compare their activity due to the lack of Envisat or S1 data (Fig. 7b).

Landslide reliability matrix
Two parameters were considered to compute the landslide reliability matrix. The average value of MP density computed is 151 MP/km 2 while the average standard deviation of the velocity is 1.13 mm/year. Applying the reliability matrix as described in the "Estimation of the time-varying landslide stability" section, the 56,133 landslides determined with S1 data are classified as follows and resumed in Fig. 8: -Landslides with a very low level of reliability represent about 9% of the total, counting in total 5033 landslides; -A low level of reliability is recognized in 8% of the landslide database analyzed (for a total of 4594 landslides); -The classes with a medium value of reliability count 27,414 landslides, equal to 49% of the whole database considered; -A high level of reliability can be assigned to 9314 landslides (16.6%); -A very high level is attributed to 17.4% of the landslides (9778 polygons).

Analysis of vertical and horizontal (east-west) components
The availability of both ascending and descending datasets provided by the P-SBAS S1 nationwide elaboration has allowed projecting the velocity vectors along the vertical (Zenit-Nadir direction) and horizontal (east-west) directions. As already depicted in LOS velocity decomposition section, values of displacement mean velocity of each MP belonging to the ascending and descending datasets have been reported on a 100 × 100 m regular grid. Thus, a synthetic MP with both values of Vv and V E within each cell is generated. Finally, IFFI landslides used in this work have been associated with synthetic MPs and only landslides with at least one synthetic MP have been used to analyze the principal components. After this step, the final selection is of 22,816 polygons. The component analysis has been performed on each type of landslide considered in this work, as previously mentioned. A heatmap-like graph has been plotted ( Fig. 9) to report the values of Vv and V E for each type of landslide, considering also the density of the data. The higher densities, as inferred by the scatter plots, are generally visible in the areas with low values of both Vv and V E ; nonetheless, prevalent horizontal component landslides can be observed, with very different distribution and velocities according to each type considered, as well as a high dispersion can be observed in some data distributions.
To also assign a numerical value to the results obtained, the mean and the standard deviation of the ratio between the absolute values of Vv and V E have been calculated for each type of landslide considered. The higher the ratio, the higher the vertical velocity, while the higher the standard deviation, the more heterogeneous the distribution.

Discussion
The IFFI project, promoted by the ISPRA along with the Italian regions and autonomous provinces, is a project aimed at providing state-of-the-art on landslides on the whole Italian territory. It is an archive of more than 600,000 landslides; however, the date of the latest update varies according to the various regions, with most of the Italian territory updated in 2007. The P-SBAS technique, an evolution of the classical SBAS algorithm, developed to deal with BigData and process entire stacks of hundredth images, has repeatedly demonstrated to be an efficient algorithm capable of providing precise deformation maps over entire national or even continental territories, especially with S1 datasets (Lanari et al. 2020;Cigna and Tapete 2021). The peculiarity of the algorithm, along with the continuous availability of SAR data, acquired with a 6-day repeat pass in the case of S1 imagery (which enables a higher detection capability), may represent a key aspect for the recurrent and constant update of the landslide inventory in a fragile territory. In this sense, Italy is known as one of the most landslide-prone countries in Europe, with more than 60% of the European landslides (Herrera et al. 2018). This assumes major importance, especially considering that in Italy every year hundreds of main landslide events occur, causing casualties and extensive damage. The application of an activity matrix, as depicted in this work, is a worthwhile methodology for the rapid and cost-effective update of the state of activity of landslides characterized by slow and very slow kinematics, as also testified by the large number of works dealing with the same aspect and using similar approaches, however on smaller scales (basin-scale mostly) (Cigna et al. 2013;Del Ventisette et al. 2014;Novellino et al. 2017;Bouali et al. 2018;Pawluszek-Filipiak et al. 2021 and many others). Such tools are easy to use and provide fast and reliable updates on every scale of analysis. However, many limitations in the use of SAR interferometry have to be yet considered. First of all, the possibility of detecting movements, which is dependent on the geometry of SAR acquisition, the revisit time, and the characteristics of the slopes under investigation (both in terms of slope geometry and land cover). Indeed, MTInSAR suffers severe limitations in the capability to measure "fast" deformation phenomena due to the ambiguous nature of its observations, i.e., the wrapped interferometric phases (Crosetto et al. 2016). If the deformation phase retrieved by InSAR between two subsequent acquisitions is bigger than π (value of the phase), the actual deformation cannot be retrieved unambiguously. The limit of π on the differential phases corresponds to a maximum differential deformation of ∕ 4 (where λ is the wavelength) over the revisit interval, depending on the satellite. In this sense, the use of satellite InSAR for detecting faster landslides is still challenging and requires a very short revisit time to retrieve reliable information.
In this work, an analysis of the slow-moving landslides with at least an MP has been made region by region, as plotted in Fig. 10, to understand which is the capability of MTInSAR, and more specifically of the systematic P-SBAS technique, to monitor landslides in all the Italian regions. In this bubble chart (Fig. 10), the ratio between the total areal extent of landslides with at least an MP and the total landsliding area is reported on the x-axis, while on the y-axis the landslide index (total landsliding area divided for the entire areal extension of a territory) is expressed. The application of the systematic P-SBAS technique, carried out in a fully automated way without any specific tuning tailored to the phenomena under investigation and with an average pixel size of 80 × 80 m, is more able to capture landslides in the Apennine regions, such as Sicilia, Emilia Romagna, and Basilicata, which are characterized by gentler slopes and a more favorable land cover setting, while it demonstrates less capacity to detect landslides located over steeper slopes, higher altitudes or in very vegetated and snow covered areas, such as in Lombardia, Veneto, Trentino, or Toscana Regions. The bubble chart also evidences the morphology of each regional territory and the reliability of the official inventories, being the regions with lower landslide index either characterized by extensive flat territories (Puglia and in lower measure Veneto Regions) or with a low number of landslides reported (mostly Sicilia, Calabria, and Friuli-Venezia Giulia Region) due to incomplete or old landslide inventories. This underlines once more the need for a constant revision of the landslide situation in Italy. In the cases of low detectability, especially when landslides occur over non-coherent areas or along the north-south direction, the definition of the state of activity of landslides presents many constraints using only satellitebased data. As aforementioned, the near-polar orbiting of SAR missions leads to low sensitivity of the north component with respect to two other components (Rocca 2003;Wright et al. 2004. Mehrabi et al. 2019, making the update of landslides state of activity very limited when moving in the NS direction.
In any case, the computation of the reliability of the data used can provide a rapid assessment of the quality of the analysis.
In addition to the abovementioned and well-known limitations of InSAR, a further point is raised in this paper. During the timespan of S1 acquisitions, central Italy was struck by a sequence of earthquakes between August and October 2016, whose deformational pattern is evident in Fig. 3. A major assumption was taken to update the state of activity of landslides within the seismic area since it is not possible to distinguish the co-seismic and landslide deformation, especially when working at such a wide scale. The displacements ongoing within the mapped landslide polygons were all interpreted as slope movements, assuming that co-seismic motion within the landslide area induces a landslide movement triggered by the earthquakes. Time series of motion can be used to track changes in the rate of motion induced by earthquakes; however, a procedure to filter out such a component of the displacement Fig. 7 a Pie chart of the landslide SoA with S1 data; b pie chart of the comparison between Envisat and S1 landslides; c bar chart of the results of the landslide state of activity matrix would require time and a not-univocal interpretation, especially when dealing with very large areas.
MTInSAR-related parameters such as MP density or the standard deviation of the displacement velocity may provide significant support for the correct interpretation of MTInSAR results, which at times can be questionable and complicated due to, for instance, a not homogeneous distribution of targets within the landslide boundaries, the possible concurrence of different displacement phenomena, or possible errors in the measurements of some point-wise targets. However, an in-depth and more detailed analysis, in specific and very peculiar cases is always recommendable. The matrix approach proposed in this paper can be suitable also for other MTInSAR data and different scales of analysis and may rapidly support the correct analysis and understanding of the results of the landslide state of activity update. Nonetheless, a further refinement of such procedure can be achieved, when dealing with more specific assessments, with Heatscatter plot for each type of landslide considered. Each dot is a landslide; the color of the dots is dependent on the density of the distribution MPs with a higher spatial distribution or with a smaller scale of analysis, as also stated by Bonì et al. (2018), in which they propose a Homogeneity index on a basin-scale update of landslides.
Finally, the assessment of the 2D components of the deformation represents more precise and comprehensive support in the estimation of the magnitude and direction of each type of considered landslide. As observed in the "Mechanical parameters for the time-varying constitutive model" section, landslides with a higher horizontal component are typically spreading, DSGSDs, and slow flows, while slides and complex landslides have also a significant vertical component, especially in the source areas, as also observed by Meng et al. (2020) and Crippa et al. (2021). The standard deviation of the various distributions for each considered type, on the other hand, highlights the heterogeneity of rotational/translational slides and complex landslides. Indeed, the concurrence of both mostly vertical-moving (rotational) and horizontal-moving (translational) phenomena and the mixing of different types of landslides (as in the case of complex ones), respectively, reflects the statistical evidence but also indicates the complexity of the MTInSAR-based interpretation, especially in cases of different landslide patterns and evolutions within the same landslide body or category. Clear improvements in this sense may be achieved with more precise data from other sources, such as unmanned aerial vehicle (UAV) or Global Positioning Systems (GPSs) and many others, however, limiting this analysis to a very detailed scale. Indeed, considering the abovementioned limitations of InSAR, the integration of satellite information with fieldwork and ground-based measurements becomes essential to obtain a complete awareness of landslide scenario in a given area, detecting also phenomena which cannot be identified through spaceborne interferometry.
Despite this, the activity and the reliability matrixes, along with the 2D component analysis, may represent a fast tool and valuable support given the upcoming deliveries of continental-scale ground motion data by the EMGS (European Ground Motion Service, Crosetto et al. 2020) which will provide millions of MPs exploiting S1 data. The update of landslide inventories, in the light of constant and periodic MTInSAR products over such a large scale, is a paramount task, serving for Civil Protection purposes and policymakers devoted to territorial planning. All these activities devoted to the definition of the state of activity and to landslides characterization will provide fundamental support for territorial management and may be implemented and further improved in view of periodical releases of MTInSAR data over wide areas. These activities can also benefit from mapping procedures, to fill the gap of non-inventoried landslides and support field activities (Festa et al. 2022).

Conclusions
The update of landslide inventory maps is a necessary step for land management and for territorial planning activities, however requiring a periodical and constant revision to provide up-to-date information on the ground displacement affecting hilly and mountain slopes. Being Italy one of the hardest-hit European countries by Fig. 10 Bubble chart of the regional distribution of the total number of landslides and the landslides with at least an MP. On the x-axis, the ratio between the total area of landslides with at least an MP and the total landsliding area is reported; on the y-axis, the landslide index is reported. The radius of the bubbles is dependent on the sum of the areal extent of landslides in a given region landslides, IFFI landslide inventory (Inventario Fenomeni Franosi in Italia, in Italian, Landslide Inventory in Italy) represents the official repository of all the landslide phenomena, which, however, contains quite outdated information. In this sense, the adoption of MTInSAR and, in the specific, of P-SBAS (Parallel Small Baseline Subset) data may represent a valuable tool for providing information on the state of activity of landslides characterized by slow and very slow kinematic. Indeed, in this work, an activity matrix has been implemented, defining the baseline state of activity of slowmoving landslides with Envisat data and updating it up to December 2018 with Sentinel-1 data processed using P-SBAS algorithm.
The results showed that about 37,000 landslides were classified as dormant (having representative velocity along the slope, V slope , lower than 7 mm/year), while 14,000 as active (V slope > 7 mm/year). Moreover, about 20,000 landslides did not display any change in the state of activity, remaining dormant or active from 2001 to 2018, while about 7000 landslides changed their state of activity (SoA); finally, for about 49,000 landslides, it was not possible to evaluate their SoA due to lack of Envisat or S1 data. In order to assess the reliability of the data, a reliability matrix has been adopted, by using the MP density and the standard deviation of the landslide representative velocity. The application of the reliability matrix on the Sentinel-1-based landslides has shown that more than 80% of the landslides retain reliability values from average to very high. Finally, the main 2D components have been determined to characterize the magnitude and direction of every type of landslide included in this work. In this case, spreadings, DSGSDs, and slow flows showed a main horizontal movement, while complex and translational/rotational slides confirmed their heterogeneity in terms of deformation direction. Therefore, the implementation of nationwide Sentinel-1 MTInSAR data may constitute a significant instrument for the constant and continuous update of landslide inventories. Moreover, the approach can also be exploited for a fast and valuable landslide deformation characterization, also considering the regular acquisition plan of Sentinel-1 mission. Considering that landslide inventories are updated until 2007 in many Italian regions, a full revision of the IFFI database would constitute valuable support to stakeholders involved in territorial planning and civil protection.