PSI Clustering for the Assessment of Underground Infrastructure Deterioration

Abstract

Remote sensing images find application in several different domains, such as land cover or land usage observation, environmental monitoring, and urbanization. This latter field has recently witnessed an interesting development with the use of remote sensing for infrastructural monitoring. In this work, we present an analysis of Sentinel-1 images, which were used to monitor the Italian provinces of Bologna and Modena located at the Emilia Region Apennines foothill. The goal of this study was the development of a machine learning-based detection system to monitor the deterioration of public aqueduct infrastructures based on Persistent Scatterer Interferometry (PSI). We evaluated the land deformation over a temporal range of five years; these series feed a k-means clustering algorithm to separate the pixels of the region according to different deformation patterns. Furthermore, we defined the critical areas as those areas where different patterns collided or overlapped. The proposed approach provides an informative tool for the structural health monitoring of underground infrastructures.

1. Introduction

The monitoring of infrastructures plays an important role for private companies and public investments, a timely detection of critic infrastructural deterioration can (in most dramatic cases) save lives and countless economic resources; an estimate of the gross domestic product fraction allocated for structural health monitoring (SHM) can be found in the several works [1,2]. The attention is generally focused on major infrastructures, such as bridges and large buildings; however, more recently, growing attention has been given to the development of monitoring systems for underground facilities. Several approaches have been proposed from different domains. The combination of blockchain technologies and IoT networks, for example, has been proposed for remote and real time monitoring [3]; other studies have pointed out the importance to monitor traffic vibrations, of course, in highly populated areas [4], which can cause severe damages, even at low frequencies. Although these advances offer novel instruments for SHM, some fundamental aspects are still left to good practices in policy, administration, and law [5].

The use of remote sensing data for SHM has a long and heterogeneous tradition in providing applications to urban areas, bridges, tunnels, railways, and dams [6,7,8,9,10,11]. The basic idea behind these applications is that Synthetic Aperture Radar (SAR) observations can be used to estimate small movements or vibrations which can eventually lead to critical deformations and collapses; a particularly dramatic example was given by the collapse of the Italian Morandi bridge [12]. Following an analogous approach, we use Persistent Scatterer Interferometry (PSI) in order to estimate the terrain displacements in the Italian provinces of Bologna and Modena. As this area is particularly affected by subsidence phenomena [13,14], it is not uncommon to observe large parts of the territory moving downwards, while the neighbor regions remain almost unaffected.

Here, we hypothesize that the frictions that are caused by inhomogeneous movements of the terrain can stress underground infrastructures and cause severe damages. In particular, we present here the challenging problem of monitoring the structural health of an underground aqueduct network, specifically with the attempt to detect early signs of deterioration. This application involves several interests; the use of satellite data, together with an automatic workflow, can improve standard maintenance activities. Besides, typical inspection activities are based on both periodical programmed checks or alerts provided by the citizen. In the first case, the inspection can be useless (if there are no damages); the second case can present an even worst scenario, as the damage found can be very expensive to fix. For this reason, the early detection of deterioration signs can help the experts to construct a priority map and yield huge economic advantages.

Using a well-established clustering technique, we segment the region of interest in different clusters according to the terrain movement; we used observations covering a time range of 5 years in order to detect robust patterns of terrain displacements. As there were no a priori hypotheses about the number of clustering classes to consider, we also performed a stability analysis, exploring different configurations, in order to determine the optimal number of classes. Finally, we designed two distinct procedures to determine which areas could be probably interested by damages: (i) on one hand, we focused on the edge regions between two or more classes; and, (ii) on the other hand, we introduced a purity measure to evaluate the presence of different classes in a particular region and selected the regions showing high inhomogeneity. Finally, we evaluated our predictions with a set of ground truth labels provided by human experts who directly observed and assessed the presence or not of structural damage.

2. Materials and Methods

The aim of this work is presenting a methodology for the structural health monitoring of underground infrastructures. Our approach is based on three distinct steps: (a) collection of remote sensing data for Persistent Scatterer Interferometry; (b) study of temporal series terrain deformations; and, (c) unsupervised clustering for detection of similarity patterns. An outline of the proposed approach is presented in Figure 1.

In particular, we examined the case of the underground aqueduct network of the Italian provinces of Modena and Bologna. This network consists of 58,392 georeferenced segments that were described by a set of 27 variables accounting for technical properties, like the material or the diameter, or geographical information, like the related municipality; although some of these features could be relevant for this study, we did not consider any of them to keep our approach based as much as possible on the sole PSI. For 47 aqueduct stretches it is also known if any maintenance was required after the observation time, which ended in 2019. Thus, while using this ground truth, we were also able to evaluate the accuracy of the proposed approach.

2.1. PSI Analysis

In this work, we used Sentinel-1 C-band images (central frequency $5.4$ GHz and wavelength $5.6$ cm), covering the Northern Italy and, in particular, the Bologna and Modena provinces, see Figure 2.

The Sentinel-1 constellation is currently composed of two twin satellites (Sentinel-1A and Sentinel-1B, respectively); the first one is active from October 2014 and the second one from September 2017, they observe the Earth from an altitude of about 693 km, at a nominal ground resolution of about $5\times 20$ m² (range × azimuth) and with a revisit time of 6 days at the equator. The study area is covered along two satellite tracks, in both ascending and descending geometries; in order to estimate the terrain movements, we used ascending observations from 30 March 2015 to 24 January 2019 and descending ones from 12 October 2014 to 25 January 2019 corresponding to 180 and 185 acquisitions, respectively. The ascending orbit consisted of ∼3.0 $\times {10}^6$ observations while for the descending orbit we had ∼3.9 $\times {10}^6$ observations; for both geometries the additional information about the Digital Elevation Model, latitude and longitude, coherence, head angle and incident angle were also available. We separately used the two geometries for two reasons: (i) the combination of ascending and descending orbit cause loss of spatial resolution; and, (ii) the combination of different geometries may hide inhomogeneous movements.

We used the SPINUA (Stable Point INterferometry even over Un-urbanised Areas) processing chain [15] to evaluate terrain displacements. PSI are well-established methods that are based on the processing of temporal series of co-registered SAR images acquired over the same target area. The main idea of PSI is that terrain movements can be derived by phase-shift differences, thus identifying a grid of movements for point-like targets (corresponding to single pixels) usually referred as Persistent Scatterers (PS). An important aspect of these measures is their phase coherence over the entire observation period; in fact, it is possible to assess the coherence with dedicated statistical analyses and neglect those measures, which do not reach a satisfactory value. For the present analyses, we selected the time-series whose phase coherence exceeded the $0.7$ threshold value [16], which ensures, in this case, an RMSE below 4 mm, for each displacement measurement.

There is another class of measurement points that can be extracted by using the SPINUA processing chain, they are named Distributed Scatterers (DS) and are coherent measurements obtained by spatial averaging statistical homogeneous pixels. The use of groups of pixels in the case of DS is dictated by the need to operate spatial averages between adjacent pixels of the same nature in order to reduce the presence of noise. This technique is particularly useful in order to extract measurement points in non-urban areas, i.e., arid areas with low vegetation and debris. On the contrary, in the case of the PS, the target inside the single pixel is sufficiently dominant with regard to the noise and, therefore, it is not necessary to further improve the quality of the phase response; this also allows for preserving the maximum spatial resolution of the target itself. In the case of DS, instead, spatial average operations introduce an unavoidable spatial resolution loss compared to the native resolution of the radar image. For that reason, only measurements derived from point-like targets have been taken into account in the present work, see Appendix A for further details.

The time series obtained for each PS within the examined grid were then clusterized with the k-means algorithm; the backbone of our approach consists in detecting similar patterns of PS and then grouping them together.

2.2. K-Means Clustering

Clustering techniques are unsupervised learning methods that basically exploit the informative content provided by the descriptive features of a dataset to outline the presence of similarities between different observations. Clustering techniques are divided into hierarchical and partitioning. The main advantage of hierarchical algorithms is that the number of the clusters directly emerges from data. However, these algorithms have a quadratic computational complexity [17,18], which restricts their application to small sets of data. This is why we used a partitioning algorithm with a linear computational complexity, like k-means.

In the k-means algorithm [19,20], the initial step consists in randomly setting N points, called centroids, whose number is equal to number of clusters that the data will be partitioned in. Of course, this is a parameter that needs to be tuned. For each cluster $C_k$, the centroid ${\overline{\mathbf{z}}}_{\mathbf{k}}$ is defined as the mean value of all observations ${\mathbf{z}}_i$ assigned to the cluster:

$${\overline{\mathbf{z}}}_{\mathbf{k}}=\frac{1}{|C_k|}\sum _{i\in C_k}{\mathbf{z}}_i$$

(1)

The choice to assign an observation to a specific cluster is taken by the minimization of a cost function $W\left(C\right)$:

$$W\left(C\right)=\sum _{k=1}^K|C_k|\sum _{i\in C_k}\left|\right|{\mathbf{z}}_i-{\overline{\mathbf{z}}}_k{\left|\right|}^2$$

(2)

which is the sum of the squared Euclidean distances between the objects ${\mathbf{z}}_i$ end their centers ${\overline{\mathbf{z}}}_k$.

Table 1. Iterative minimization procedure of the function of cost W.

Step 1.	Initialize the set of clusters $C=\{C_1,C_2,\cdots C_K\}$
Step 2.	Compute the new centers in the feature space:
	${\overline{\mathbf{z}}}_k={\sum }_{i\in C_k}{\mathbf{z}}_i/\left|C_k\right|$
Step 3.	Search for the set of clusters $C^{\ast }$ such that $W\left(C\right)$ is minimum
	i.e.,
	$C^{\ast }{\ni }^{\prime }{C}^{\ast }\left(i\right)=arg{min}_k\left(\left|\right|{\mathbf{z}}_i-{\overline{\mathbf{z}}}_k\left|\right|\right)$ per $i=1,2,\dots N$
Step 4.	If $C\left(\right)\ne C^{\ast }\left(\right)\wedge |W\left(C\right)-W\left(C^{\ast }\right)|>ϵ$,   repeat steps 2 and 3.
	Otherwise stop iterating.

outlines the iterative procedure adopted to minimize the function of cost W.

The main drawback of k-means is that its partitioning is strongly affected by the random initialization of centroids. For this reason, we performed 10,000 repetitions and averaged their results in order to obtain our final partition. Finally, as the examined time series consist of terrain displacements over time, the patterns detected by k-means can physically be interpreted as group of PS moving with the same velocity of the cluster centroid.

We explored a number of k clusters varying from 2 to 10. To evaluate the clusters’ separation, we used three distinct indicators:

• the between-sum-of-squares/total-sum-of-squares ratio ($BSS/TSS$);
• the Adjusted Rand Index ($ARI$); and,
• the Silhouette statistics S.

The $BSS/TSS$ ratio basically compares the distance of points within a class with all the distances in the data; the $BSS/TSS$ ratio ranges from 0 (perfect overlap) to 1 (disjoint clusters). The ARI indicator is a measure of agreement between two partitions of the same dataset. It is used here in order to measure the agreement between different k-means initializations. The ARI indicator ranges from $-1$ to 1. Perfect agreement corresponds to 1, when $ARI=0$ the agreement can be considered to be a random occurrence, while $-1$ denotes a perfect disagreement of the clusters. Finally, the silhouette S measures the clustering consistency, in general a good experimental value should lie around $S=0.5$, further details and references regarding these metrics are presented in Appendix B.

Once each PS within the observation grid was assigned to a cluster, we were able to detect critical areas.

2.3. Detection of Critical Areas

The underlying hypothesis of our approach consists in assuming that terrain deformations stress underground infrastructures, eventually causing severe damages or disruptions that require maintenance. Accordingly, we use the clustering analysis to detect which are the most probable regions where these phenomena can occur. In this work, we examined two distinct situations. On one hand, it is reasonable to assume that aqueduct stretches whose ends lie in different clusters; therefore, moving with different velocities will be subject to tension forces and frictions which account for the measured velocity difference and could yield infrastructural damages. This approach basically consists in determining a spatial gradient. On the other hand, even if the ends of a stretch belong to the same cluster, it is possible that the stretch experiences tensions and frictions, for example, because its central part lies in a region moving with different velocity; a figurative representation is shown in Figure 3.

The latter case presents infinite possibilities to examine; nevertheless, we can model the problem with a slight change of perspective. Instead of focusing on the aqueduct stretch, we consider a region surrounding the stretch, which we name buffer, and compute this region purity in terms of the number of PS falling into this region and the number of different clusters they belong to.

Let us examine in detail the gradient approach. Once the PS have been assigned to a cluster we have a grid of labeled points over the region of interest. As previously mentioned, physically we can interpret this map as a discrete velocity field. In order to obtain a continuous velocity field, the most adopted choice would be the kriging method [21], but the computational requirements suggested us an alternative strategy. Here, we used another consolidated interpolation procedure, the inverse distance weighting (IDW) interpolation [22]. Given a set of discrete observations ${\mathbf{x}}_{\mathbf{i}}$ with $i=1,2,\dots ,N$, it is possible to define the interpolating function $f\left(\mathbf{x}\right)$, such that:

$$f\left(\mathbf{x}\right)=\left\{\begin{array}{cc}{f}_i,\hfill & if\exists i{\ni }^{\prime }\mathbf{x}={\mathbf{x}}_{\mathbf{i}}\hfill \\ \frac{{\sum }_{i=1}^N{w}_i\left(\mathbf{x}\right)f_i}{{\sum }_{i=1}^N{w}_i\left(\mathbf{x}\right)},\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(3)

where $w_i$ is a similarity function usually dependent on the reciprocal of a distance, as:

$$w_i\left(\mathbf{x}\right)={\left(\frac{1}{d({\mathbf{x}}_{\mathbf{i}},\mathbf{x})}\right)}^p$$

(4)

where ${\mathbf{x}}_{\mathbf{i}}$ is the i-th observation and $d(\mathbf{x},{\mathbf{x}}_{\mathbf{i}})$ is its distance from $\mathbf{x}$. In this study we adopted $p=3$ to shorten the spatial influence of the measured grid. Finally, the gradients ${\nabla }_{x}f$ and ${\nabla }_{y}f$ of the interpolated field f along the x and y directions were estimated with a $5\times 5$ Sobel filter. In fact, we observed that, with a $3\times 3$ filter, the edges appeared affected by a salt-and-pepper-like noise, while greater windows did not yield any improvement in both the image quality or in computational times.

Finally, the magnitude of the gradient $|\nabla f|$ is computed, as follows:

$$|\nabla f\left(\mathbf{x}\right)|=\sqrt{{\left({\nabla }_{x}f\left(\mathbf{x}\right)\right)}^2+{\left({\nabla }_{y}f\left(\mathbf{x}\right)\right)}^2}=\sqrt{{(F_x\ast f\left(\mathbf{x}\right))}^2+{(F_y\ast f\left(\mathbf{x}\right))}^2}$$

(5)

The status of the segments of the water-supply network is then considered to require an inspection for maintenance by comparing the gradient value with a threshold $|{\nabla }_{TH}|$.

For what concerns the buffer approach, the initial step, analogously to what we did with the gradient approach, consists in considering the PS clustering. In this case, there is no need for spatial interpolation, as we are only interested in evaluating the PS near an aqueduct stretch. To gain statistical robustness, we observed that a consistent number of PS was obtained for each stretch if we considered a buffer area (centered on the aqueduct stretch) having a 30 m radius. We considered the distribution of PS velocities within the buffer area and determined their variance, as a measure of homogeneity, the less the variance the greater the homogeneity; therefore, a critical aqueduct stretch was pointed out when the variance exceeded a threshold value ${\Delta }_{TH}$.

It is worth noting that both the gradient and buffer approaches yield a score that can easily be interpreted as a classification score; in fact, the gradient itself and the variance within the buffer region can be thresholded, using $|{\nabla }_{TH}|$ and ${\Delta }_{TH}$, respectively, resulting in a binary response about the aqueduct stretch, damaged or not. It is important to note, to underline the clustering robustness, that these thresholds affect the classification of the aqueduct stretch, not the k-means partitioning. These thresholds could be tuned according to several different considerations, in particular, as maintenance inspections are usually expensive both in terms of time and cost, we needed to consider the trade-off between the possibility to deliver warnings when maintenance was unnecessary (false positives) and neglect critical situations (false negatives); in other situations, for example, when large databases are available both for validation and testing accordingly, which is not the case here, these thresholds could be tuned according to classification performance. Thus, as performance metrics, we used in this work the Area Under the receiver operating characteristic Curve (AUC), which does not depend on the choice of a specific threshold value.

3. Results

3.1. How Many PS Patterns?

Our first analysis aimed at detecting how many distinct PS patterns could be detected in the region of interests. Figure 4 shows the comparison of the BSS/TSS ratio and the ARI for ascending orbit.

We found that the BSS/TSS ratio reached the higher values from a number of clusters $k\ge 5$; on the contrary, the ARI showed the best separation of clusters when $k\le 5$. Therefore, these two metrics would suggest using $k=5$. Besides, we investigated the silhouette S, for $k=5$, we found $S=0.48$. We also observed that, except the case $k=2$ showing $S=0.63$, according to the silhouette metrics all the examined k values show a good consistency ($S\simeq 0.45$). Figure 5 presents the results.

Therefore, the subsequent analyses were performed when considering 5 different clusters. It should be noted that, when repeating the analysis with the descending geometry, we obtained the same result.

3.2. The Physical Interpretation of Clusters

An indirect validation of the previous result about the number of detected PS patterns was obtained by inspecting the distribution of the average line of sight (LOS) velocities within each cluster, see Figure 6.

We estimated the velocity of each PS in a cluster and therefore the cluster velocity distribution. We observed a nice separation for these distributions, so that it was possible to consider each cluster physically associated to a well-defined terrain movement in terms of average cluster velocity: “very high subsidence” ($-9.5$ mm/year), “high subsidence” ($-5.3$ mm/year), “natural subsidence” ($-2.6$ mm/year), “stable” ($-0.5$ mm/year), and “stable or uplifting” ($+1.2$ mm/year). Figure 6 also shows the percentage of PS points belonging to each cluster.

A visual inspection of the region of interest shows that the clusters have typical dimensions of about 2∼10 km, see Figure 7.

All of the detected patterns have in general a spatial coherence so that it is difficult, for example, to find an uplifting region whose edges border on a very high subsidence region. Interestingly, at small regions (<1 km) a substantial inhomogeneity emerges, in these small regions is not unusual to find PS belonging even to 3 or 4 different classes.

3.3. Damage Detection on Underground Aqueduct Stretches

In order to validate our pipeline, we collected and examined 47 on field assessments provided by human experts. These evaluations showed in 44 cases the presence of a structural damage while in 3 cases there was no manifest need of intervention. Figure 8 presents the predictions obtained while using the ascending orbit data.

Both the gradient and buffer methods performed better than chance with AUC $=0.72\pm 0.13$ and AUC $=0.78\pm 0.11$, respectively. The reported standard error was calculated with the Hanley formula [23]. At the operating points, the gradient method returned sensitivity $65.9\%$ and specificity $100\%$ while we obtained sensitivity $81.4\%$ and specificity $66.7\%$ for the buffer method.

According to a z-test, no significant difference can be detected between the two methods. A magnified view of a subsidence region with a detected damage is presented in Figure 9.

3.4. Replication Study on the Descending Orbit

We repeated our classification analyses with the descending orbit data in order to confirm the robustness of the performed analyses and the performance accuracy, especially for the damage detection, see Figure 10.

Even in this case, we observed no significant difference between the two approaches with a z-test. The same conclusions hold for a comparison between the descending orbit and ascending orbit results.

4. Discussion

The analyses determining the optimal number of clusters were entirely based on the concepts of stability and reproducibility; accordingly, we explored three different metrics: BSS/TSS ratio, ARI and silhouette S. While BSS/TSS and ARI suggested that the best choice consisted in 5 clusters the S best value was obtained for $k=2$. However, it is worth noting that S values for all partitions were ≥0.45, which means that the clusters are consistent, even if there is not a huge separation. These results suggest that a separation based only on PSI data can be affected by a considerable uncertainty; time series covering a wider temporal range or dedicated methodologies could, in principle, improve the clustering results [24,25,26]. Nevertheless, it is worth noting that, despite 10,000 different initializations, the clusters remained stable and their classification results reproducible while using two distinct datasets, the ascending and the descending geometries.

The clusters that were detected from our analysis have a nice physical validation. By inspecting the distribution velocities within each class, we found that each cluster could be characterized by a different LOS velocity. This is somehow expected, in that a reasonable clustering result should detect different terrain movements. In particular, we found that using the two available geometries, the optimal number of classes, 5, remained unchanged; this should not be taken for granted considering that the presence of horizontal (East-West) displacements can significantly change the LOS measurements obtained with the two geometries. In this regard, it must be considered that the region of interest is affected by tectonic motions in the NE direction (with velocity ranging from 3 mm/year to 6 mm/year) [27,28] which can yield significant changes to estimated velocities [29]. This could explain the different average within-cluster LOS velocities obtained in the two geometries.

Finally, we have designed two different approaches, the gradient and the buffer methods, to evaluate whether terrain movements estimated by PSI can predict structural damages. We evaluated the classification accuracy of these methods while using a labeled dataset of maintenance interventions on the aqueduct stretches of the Italian provinces of Bologna and Modena, Northern Italy. We had 47 maintenance interventions by human experts that evaluated the presence or the absence of a structural damage on an underground aqueduct stretch. The overall accuracy in terms of AUC was higher than $0.70$ for both methods and for both ascending and descending orbit data. Our findings were statistically indistinguishable.

Despite the good level of accuracy, there are some critical aspects that should be taken into account: first of all, the limited size of labeled data; this prevents us from obtaining more precise estimates of the accuracy, in fact the standard errors obtained using the Hanley formula are $0.11$–$0.13$. Another aspect to be considered concerns the classification errors that can affect this kind of analysis. In general, classification errors are of two type, false positives and false negative. A performance analysis in terms of sensitivity and specificity showed how, despite an analogous classification performance, the two proposed methods have completely different behaviors, while the gradient is more accurate for the detection of not damaged stretches, the buffer on the contrary is more accurate for the damage detection. Future studies could address the possibility to combine these approaches and even improve the classification performance. In this analysis, we only examined the terrain movements as a possible cause of structural damage, but, of course, this is an oversimplification. Neglecting other possible damage causes affects the methodology preventing the detection of damages which therefore are false negatives. At the same time, the frictions and forces developed by in-homogeneous terrain movements can yield micro-damages which are not detected by human experts at the present but can accumulate over time eventually ending with an overt structural damage; these situation represent the false positives of this analysis. Based on these considerations, we are aware that the estimated AUC should be considered cum grano salis; nevertheless, our findings suggest that a structural health monitoring of underground infrastructures could be founded on PSI analyses.

As previously mentioned, the presented results would benefit from a larger labeled dataset. Another possible improvement concerns data resolution. Sentinel-1 resolution is adequate for the detection of the pipelines that are subject to ground motion, as the PS resolution ($5\times 20$ m²) matches the typical dimensions of aqueduct stretches and reduces the regions to monitor, a fundamental aspect for operative requirements. In principle, high-resolution data (like Cosmo-SkyMed and TerraSAR-X) could improve the methodology. However, in Italy, these missions cannot ensure the systematic global coverage and the same revisit time provided by Sentinel-1; despite these limitations, future studies should evaluate to which extent high-resolution data can improve the proposed approach.

Finally, it could be of fundamental importance to design a framework of analysis that is able to manage heterogeneous data, not only from PSI, to take into account damages from other sources. With this regard, we aim at collecting more data and develop alternative strategies exploiting supervised learning [30,31,32,33]. In particular, novel graph-based approaches [34,35] and deep learning strategies [36,37,38] seem to be able to bring significant advances to the field. These approaches are naturally suited to manage large and heterogeneous datasets, therefore they could represent the best options to improve the findings that are presented in this work and pave the way for future applications in underground SHM.

5. Conclusions

In this work, we proposed a PSI processing pipeline for the monitoring and automated detection of underground infrastructural damages, specifically we investigated the case of aqueduct stretches in the Italian provinces of Bologna and Modena, Northern Italy. Based on the well-established analysis techniques of PSI and unsupervised clustering k-means, the proposed approach detects coherent PS movements and clusters them in groups whose physical interpretation is given in terms of similar LOS velocity. The border regions between different clusters and the small inhomogeneous regions that include PS of from different clusters are considered potentially dangerous for infrastructures; in these regions, the frictions and forces that arise from the different terrain velocities can weaken the structural health of the aqueduct stretches. The proposed approach exploits the Sentinel-1 data. The use of high-resolution data or the combination of PSI with in situ measures could represent a great opportunity to further improve the reliability of the methodology; however, the method attested an accuracy in terms of AUC $\sim 75\%$. Besides, it should not be underestimated the fact that, the current high-resolution missions have two major drawbacks related to (i) the economic cost of the data, that in some cases can be difficult to afford and (ii) the reduced spatial coverage capability, which prevents a systematic acquisition of data with global coverage and with a short revisit time, as instead happens with the Sentinel-1 mission. The reported performance can be largely affected by other factors, such as an increment of the available data. An intrinsic and unavoidable limitation is caused by the presence of damages that are caused by other sources and the damages whose signals cannot be detected from the proposed approach. Nevertheless, we demonstrated that the proposed approach is able to detect the relation between terrain movements and underground structural damages.