ASSESSING RESILIENCE OF INFRASTRUCTURES TOWARDS EXOGENOUS EVENTS BY USING PS-INSAR-BASED SURFACE MOTION ESTIMATES AND MACHINE LEARNING REGRESSION TECHNIQUES

Technologically advanced strategies in infrastructural maintenance are increasingly required in countries such as Italy, where recovery and rehabilitation interventions are preferred to new works. For this purpose, Interferometric Synthetic Aperture Radar (InSAR) techniques have been employed in recent years, achieving reliable outcomes in the identification of infrastructural instabilities. Nevertheless, using the InSAR survey exclusively, it is not feasible to recognize the reasons for such vulnerabilities, and further in-depth investigations are essential. The primary purpose of this paper is to predict infrastructural displacements connected to surface motion and the related causes by combining InSAR techniques and Machine Learning algorithms. The development and application of a Regression Tree-based algorithm have been carried out for estimating the displacement of road pavement structures detected by the Persistent Scatterer InSAR technique. The study area is located in the province of Pistoia, Tuscany, Italy. Sentinel-1 images from 2014 to 2019 were used for the interferometric process, and a set of 29 environmental parameters was collected in a GIS platform. The database is randomly split into a Training (70%) and Test sets (30%). With the Training set, through a 10-Fold Cross-Validation, the model is trained, validated, and the Goodness-of-Fit is evaluated. Also, with the Test set, the Predictive Performance of the model is assessed. Lastly, we applied the model onto a stretch of a two-lane rural road that crosses the area. Results show that the suggested procedure can be used for supporting decision-making processes on planning road maintenance by National Road Authorities. * Corresponding author


INTRODUCTION
One of the most prominent challenges in infrastructure Pavement Management Systems (PMSs) is timely detection for applying preventive actions and early recovery. Indeed, accurate planning of infrastructure maintenance enhances the service life and reduce the total maintenance cost. Non-Destructive High-Performance Techniques (NDT), such as network-scale monitoring devices, are required for allowing the identification of deficiencies and reducing the number of inspections directly to the sites. Specifically, SAR-based systems enable detection of surface motions on and around infrastructures. Therefore, such systems can increase the effectiveness of the maintenance planning (Fagrhi and Ozden, 2015). The scope of this research is to evaluate the possibility of using SAR-based systems for planning the strategy of infrastructure maintenance. In combination with SAR, we want to define a procedure that exploits the capability of advanced statistical modeling, such as Machine Learning (ML) techniques, that can associate conditioning factors with the target variable examined. In this research, the target variable is the surface motion detected by a SAR sensor, while factors are those related to exogenous events of the infrastructure. Exogenous events can cause modification in infrastructure conditions and are connected to major extreme natural events, such as earthquakes, landslides, subsidence, sinkholes, and floods.
These events can be related to Topological, geomorphological, geomorphometric, hydrological, and social systems of the surrounding environment of the infrastructure. Being able to correlate factors and surface motion enhance the possibility of defining a maintenance strategy accurately, identifying the most appropriate interventions that infrastructure requires. Therefore, two main research questions arise:
What are the exogenous factors most affecting the instabilities of infrastructures?
This paper attempts to answer to both research questions under the assumption that surface motion estimations are intended as a proxy of infrastructure instabilities. As regards the first research question, the Predictive Performance of the Regression Tree-based (RT) model demonstrates the correlation between surface motion and conditioning factors.
For what concerns the second research question, the RT model also provides an estimation of the predictor importance, identifying the features that most influence infrastructure instabilities.
Being able to answer both questions, we can validate the assumption on which this paper is based.

InSAR Techniques: D-InSAR and PS-InSAR
Differential InSAR (D-InSAR) has been promoted to measure gradual surface motion between a couple of SAR image acquisitions. It has been used broadly for subsidence monitoring (Benattou et al., 2018;Del Soldato et al., 2018;Rosi et al., 2016;Solari et al., 2018), landslides (Bianchini et al., 2015;Wasowski and Bovenga, 2014), and sinkholes (Hoppe et al., 2016). D-InSAR allows monitoring large areas (Costantini et al., 2017), at relatively low cost with accuracy within centimeters to millimeters. Nevertheless, DInSAR is influenced by two primary sources of errors: temporal and geometric decorrelation and phase distortion due to atmospheric conditions, which might decrease the monitoring accuracy.
In order to overcome the weaknesses of this technique, Persistent Scatterer InSAR (PS-InSAR) has been developed to map surface motion over time, exploiting a long stack of coregistered SAR image acquisitions (Ferretti et al., 2001(Ferretti et al., , 2000. The PS-InSAR technique relies on the use of the socalled Persistent (or Permanent) Scatterers (PS), that are onground items for which spectral response does not change substantially during various SAR image acquisitions. Phase information backscattered from these coherent PS is used to determine the magnitude and the temporal evolution of the surface motion along the Line of Sight (LOS) of the SAR sensor. The displacement values for each PS are expressed as a velocity in millimeters for a year. However, the technique has some usage limitations; indeed, the primary weakness of this technique is the almost-absence of PS points in agricultural and forested areas, where variations in geometry between SAR acquisitions can cause phase decorrelation, and make the selection of PS difficult (Solari et al., 2016).
Conversely, PS-InSAR shows its best performance in recognizing PS on urbanized areas and infrastructures (both roads and railways).

PS-InSAR for Infrastructure Maintenance Planning
The last decade has seen the increased use of InSAR techniques in the field of maintenance of infrastructures. Studies made use of InSAR techniques for identifying areas in which to build new infrastructures (planning) (Balz and Düring, 2017), while other studies use InSAR techniques to be able to prevent any critical infrastructural deficiency (prevention), identifying sections in which there are substantial movements (Bakon et al., 2014;Wasowski and Bovenga, 2014). Generally, studies are focused on application of InSAR techniques both for roads (D'Aranno et al., 2019;Murdzek et al., 2018;Wasowski et al., 2017;Xing et al., 2019), railways (Peduto et al., 2017), and tunnels (Perissin et al., 2012). In (Xing et al., 2019), the authors used InSAR techniques for understanding the seasonal movements of a highway in China, providing different types of maintenance interventions based on the magnitude of the detected surface motion. Also, InSAR has been studied under the point of view of economic sustainability. Indeed, Ozden et al. (Ozden et al., 2016) proposed a benefit/cost analysis of InSAR analyses for infrastructure, stating that SAR-based monitoring is useful as a complementary tool to improve the effectiveness of overall monitoring system and reduce the total cost.

Environmental Modelling with ML
ML algorithms have been widely used for modeling of environmental phenomena and attempting to understand the reason of why they happen. It is possible to find applications of ML for regression or classification tasks in an extensive set of topics. Classification tasks are those related to the prediction of discrete values, i.e., to assign a label (or class) to the response variable depending on a set of independent variables (also called features of factors). Conversely, regression tasks, are those related to the prediction of continuous values for the response variable.
As regards the purpose of this paper, bibliographic research related to ML has been carried out on the fields of extreme natural events modeling. The purpose of this phase was to identify the potential set of features that can be related to surface motions, the ML models mainly used, the procedure for developing the models, and the performance metrics used for evaluating the Goodness-of-Fit and the Predictive Performance. So far, there are no studies that involve the prediction of PS-InSAR measurement based on environmental features; then, bibliographic research has been focused on recent studies that account for surface motion modeling and occurring of natural phenomena. Several studies exploit ML techniques to model extreme natural events such as occurring of landslides (Al-Najjar et al., 2019;Dou et al., 2019;Hong et al., 2017;Tien Bui et al., 2016;Xie et al., 2017), gully erosion by stream power (Arabameri et al., 2019(Arabameri et al., , 2018Gayen et al., 2019), occurring of floods (Cian et al., 2019;Khosravi et al., 2019;Rahmati et al., 2019b), and land subsidence (Rahmati et al., 2019a). In all these studies, considering the geospatial nature of the information used to develop the models, the authors proposed susceptibility maps for providing the predictions. The task of each of the mentioned paper was to predict a class. For example, landslides modeling usually attempts to predict if each pixel of the susceptibility map has to be a "landslidepixel" or "non-landslide-pixel" (binary classification). Conversely, regression ML models have been used for the prediction of Safety Factor of slope stability (Bui et al., 2019;Samui, 2008). In such studies, authors attempted to predict a numeric value. Results are proposed as a Scatterplot, in which observation and prediction are compared onto a Cartesian Plane. Also, a set of performance parameters are generally used to assess the Goodness-of-Fit and Predictive Performance of the ML models, such as Correlation Coefficient (R 2  There is not a rule of thumb in determining the most appropriate model. Indeed, it depends on the framework of each different study. Accordingly, this paper proposes the calibration and application of a RT-based model. It was chosen considering its reasonable simplicity of implementation, its capacity of handling numerical and categorical input factors, its computational efficiency (time to train the model), and the possibility of quantifying the importance of the input factors that can potentially influence the surface motions. The outcomes have to be considered as a first result of the suggested procedure, that can be undoubtedly improved by developing and comparing more type of ML models.

STUDY AREA DESCRIPTION
The study area extends for 11 km x 12 km = 132 km 2 (red rectangle in Fig. 1 -b). It is located in the Province of Pistoia (green areas of Figure 1), in the North of the Tuscany Region ( Fig. 1 -a), central Italy. Other authors investigated the same area (Rosi et al., 2016) because of the strong subsidence effects that affect the city of Pistoia located at the center of the study area. The reason for such effects has been evaluated in , showing that they are related to the overexploitation of groundwater by several greenhouses and nurseries opened from the beginning of the 21 st century. The combination of the high request of water and soft layers compaction could have caused a decrease in the groundwater level, resulting in subsidence. Figure 1 also shows three regional two-lane rural roads that cross the study area. Twolane rural roads are identified as single carriageway roads with one lane for each direction of travel. In Italy, the width of the lanes should be 3.75 meters, the paved shoulder 1.50 meters, the radius of the circular curves at least 118 meters, and the maximum slope of 7%.

RT-based ML Algorithms
A RT-based ML algorithm is a hierarchical supervised learning model composed of decision rules that recursively split independent variables into homogeneous zones. The decision rules are learned by inferring directly from the available data. Once a Tree-based model is trained, the decision rules can be used to make predictions by using a new set of data. Tree-based algorithms are used both for classification and regression tasks, if the target outcomes are discrete or continuous, respectively. Historically, RT algorithms are introduced and defined in (Breiman et al., 1984) under the name of Classification and Regression Trees (CARTs). In this paper, the same modeling procedure is used. In (Torgo, 2017), advantages and limits in using RT have been described, and below they are briefly reported.
Features of RT make it an exciting approach in modeling multiple regression problems, also finding an easy interpretable solution for the users. Indeed, RT provides automatic variable selection making them insensitive to irrelevant variables, computational efficiency even in large problems, handling of unknown variable values, handling of both numerical and categorical predictors, insensitivity to the scales of predictors, and interpretable models by using treebased graph visualization (Torgo, 2017). As regards for its limits, RT could have poor prediction accuracy considering the piecewise constant approximation. Also, RT could be unstable: small changes in the training data can lead to significantly different models.
The RT used in this study is binary. Therefore, we refer to a Tree-based model in which the root node and each branch node is split into two different branches according to a specific decision rule. In order to define the decision rules, the Recursive Partitioning (RP) algorithm (Breiman et al., 1984;Loh and Shih, 1997) has been used. Starting from the root node, RP splits the node into two branches by evaluating all the possible cut point values for each predictor, identifying the split that minimizes the so-called splitting criterion. In this case, the splitting criterion is the Mean Square Error. Therefore, according to the MSE criterion, the error in a given node is given by: Where: Dt = sample of cases in node t nt = cardinality of this set y-bar = average target variable of the cases in Dt xi = sample of predictors in Dt yi = sample of target variables in Dt A logical test s split the cases in Dt into two partitions, DtL and DtR. The resulting error of such split s is given by: ( Where: ntL/nt and ntR/nt = proportions of cases after the split s that move to the left and right child nodes of t, respectively Therefore, it is possible to assess the goodness of the split s by the relative error reduction: Finding the best split test s for a node involves evaluating all possible tests using Equation 3, that is to evaluate all possible cut point for each predictor and compute the relative error reduction. As said, for each new child node, the RP has been applied again. This procedure moves on until a termination criterion is satisfied. In this study, the termination criterion is a Mean Square Error threshold equal to 10 -6 . Once the termination criterion is satisfied, the node becomes a leaf node. The estimation of the target variable is the average value of all the yi included in the leaf node.

ML Modeling Procedure
The ML approach involves three main modeling phases: training, validation, and testing phase. Therefore, the database has been randomly split into two parts: the training set (70%) and the test set (30%). In order to develop a model able of better generalizing and avoiding overfitting, a 10-Fold Cross-Validation was carried out; the training set is then divided into 10 folds. The model is trained (i.e., the set of hyperparameters is estimated) using 9 folds and is validated on the remaining one. The process is repeated 10 times, allowing each observation in training set to be used for both training and validation. The set of hyperparameters that provide the least MSE in predicting the target variable constitutes the set of the best hyperparameters of the model. The following are the initial Hyperparameters, the modeling process setup, and its evaluation metrics.  The maximum number of splits is the maximum number of decision splits (or branch nodes). The minimum of leaf size is the minimum number of leaf node observations. The minimum of parent size is the minimum number of branch node observations. Also, RT merges leaves that originate from the same parent node and yield a sum of risk values greater than or equal to the risk associated with the parent node. The Goodness-of-Fit evaluation is necessary for assessing how well the model fits the training dataset. However, such performance cannot be used for evaluating the prediction and generalization abilities of the model. Therefore, the Predictive Performance of the model has been evaluated using R 2 , RMSE, MAE, using the testing set instead of the training one.

RT input factors: Environmental Features
A set of environmental features were collected and preprocessed in a GIS platform in order to obtain the RT input factors related to Topography, Geomorphology, Geomorphometry, and Hydrology. The Tuscany Region Administration provided them through its Geoscope (http://www502.regione.toscana.it/geoscopio/cartoteca.html). Below the factors are reported together with their nomenclature (in round brackets), their processing (Gathered, G, or Computed through GIS, C), type of input factor (Numerical, Num, Categorical, Cat, Ordinal, Ord, or Binary, Bin), unity of measure [in square brackets] and, if possible, a reference for their definition and use in ML modeling. It is worth mentioning that all the gathered data have been provided by map with a scale equal to 1:10.000.
Elevation ( In order to solve the different semantic relationships between the data, we made use of the cell-based structure of the Raster format, which allows homogenizing the various sources. Indeed, the collected features were both Shapefiles (Point, Polyline, and Polygon information) and Rasters. All the resulting input were sampled with the same resolution (10 meters x 10 meters), thus being able to stack the input information and proceed with the definition of the database. The size of the Rasters is 1100 cells x 1200 cells.
Point-based Shapefiles (Rain gauges and Earthquakes localization): they have been transformed into Raster through the use of geospatial tools: occurred earthquakes (their magnitude), and rainfall (cumulative yearly rainfall) caught by rain gauges have been interpolated by spherical ordinary kriging (Xie et al., 2017) to generate the Seismic susceptibility and Cumulative Average Yearly Rainfall maps, respectively. Polyline-based Shapefiles (river graph): river graph has been transformed into a Raster by measuring the Euclidean distance from the rivers. It has also been computed the River Density, converting the river graph into pixels (obtaining river-and non-river-pixels), and then calculating the density of riverpixels/km^2.
Polygon-based Shapefiles (landslide localization, drain, floodSusc, erosionSusc, lsSusc, LU, UA, sand, silt, clay, org): landslide map has been transformed into a Raster by measuring the Euclidean distance from the landslides. Flood susceptibility, erosion susceptibility, drainage capacity of the soil, landslide susceptibility, and Land Use map were converted into Raster preserving their ordinal or categorical values. Urban areas map was converted into a binary Raster: each cell has to be 0 (absence of urban area) or 1 (presence of urban area). Shapefiles related to the composition of the subsoil were converted into Raster, preserving their numerical values.
Raster-based files (Elev, aspect, slope, curv, conI, SL, TPI, VTR, TRI, TWI, SPI, difIns, dirIns, WE): Elevation has been collected as a Raster with a resolution of 10m x 10m. The other input factors that derive directly from the Elevation are all Rasters with the same cell resolution.

RT output: Average Velocity Estimation
As regards the interferometric process, we used an outcome provided by TRE Altamira (https://site.tre-altamira.com) and available for free on https://geoportale.lamma.rete.toscana.it. They used a stack of co-registered Sentinel-1 ascending orbit acquisitions that cover a period from 12/12/2014 to 24/08/2019. The stack is composed of 210 images. An amount of 52257 PS (about 400 PS/km 2 ) with a coherence greater than 0.9 were selected for developing the model. In order to develop the mode, an interpolated surface of PS average velocity is required since PS-InSAR provides a discrete point-sampling outcome, and PS are not distributed onto all cells of the input factors. Therefore, an Inverse Distance Weighting (IDW) procedure has been applied (Bianchini et al., 2015;Peduto et al., 2015) for computing this interpolated map. The IDW technique is based on the principle that neighboring elements have a more significant correlation than distant ones. Consequently, IDW predicts a measurement for each cell in which there are no PS by exploiting those present, assuming that each measured point has a radial influence that decreases with distance.

Regression Tree Hyperparameters
Once the RT has been trained, the final Hyperparameters can be reported. They are described below: The number of nodes of the RT is 12865, and the hierarchical levels are 5675. The computational time required for training the model has been 5760 seconds.

Surface Motion Prediction
In order to realize the predicted map of the surface motion, we applied the trained RT onto the whole study area. Figure 2 shows both the map of the observed surface motions and the predicted one, which is the outcome of the RT. It is worth to remind that the observed map has been obtained after an IDW procedure.
(a) (b) Figure 2. Surface Motion, (a) Observation, (b) Prediction Qualitatively speaking, the two maps seem to be very similar in the whole study area, both in the shapes and in the magnitude of subsidence effects. The training and testing metrics show similar performance in both phases. This fact should ensure that the modeling has been carried out correctly. Indeed, a strong decrease in performance from the training phase to the testing one could have revealed a potential overfitting issue.

Residual Map
By taking advantage of the spatial nature of the input and output factors, it is possible to analyze the residual errors from a spatial point of view. The residual error is defined as the difference between the prediction of the model (Fig. 2 -b) and the observation (Fig 2 -a). The residual spatial distribution (Figure 3) can be useful to observe if, where, and how much the model overestimates or underestimates an observation.

Figure 3. Residual map [mm/year]
Qualitatively, observing Figure 3, most of the area is affected by residues whose value is around zero. The error peaks are due to the more urbanized areas, where there are a considerable number of PS. In these areas, the model is unable to predict the slightest difference between PS measurements. Moreover, from a quantitative point of view, the residual frequency distributions have been computed, both for the training and testing phase. The resulting distributions were Gaussian, with a mean equal to 2.3108e-5 (Training) and 2.4410e-4 (Testing), respectively. The standard deviations are equal to the RMSE parameters showed before.

Two-lane Rural Road Application
A stretch of a regional two-lane rural road that crosses the study area was chosen as a test site for the application of the model. Once the predicted map has been defined ( Figure 2 b), and under the assumption that road damages are related to surface motion estimations, road instabilities can be detected overlaying the predicted map and the road graph. Figure 4 shows the road condition, where "condition" has to be intended as road average yearly displacement.

Predictor Importance
As said, RT allows evaluating the importance of the predictors used in the modeling process. The importance related to a split of each node due to a specific factor is computed as the difference between MSE for the parent node and the total MSE for the two child nodes. Therefore, Predictor Importance (PI) is evaluated by computing all the increases in MSE due to splits for each predictor and dividing the sum by the number of branch nodes. The higher is the increase in MSE, and the higher is the importance of the factor considered. Figure 6 shows the standardized PI of the RT. Figure 6. Predictor Importance Figure 6 demonstrates how much subsidence is related to the various input factors. Nonetheless, it should be noted that the term "Importance" is related to the splitting process of the RT. If a factor is used early to split a node compared to another one, it gains more importance, even if it is not strongly related to the dependent variable. Therefore, Figure 6 has to be analyzed carefully and qualitatively. However, it seems that subsidence is more affected by hydrological factors than the other ones, considering that Cumulative Average Yearly Rainfall, distance from rivers, river density, and land use (that is associated with the overexploitation of water by nurseries) fall in the top zone of the ranking. The most important factor is the distance from landslides. Also, elevation and wind effects have significant importance. Concerning the composition of the subsoil, it seems that the percentage of organic substances occupies the dominant position.

Use of the Procedure by Road Authorities
Once the model has been tested and the reliability of the output is verified, it can be employed by Road Authorities on new data in different ways. Principally, it is possible to predict the subsidence rate in areas where there is not enough PS data. Being able to produce a comprehensive mapping of the subsidence rate is essential for drawing up a list of on-site inspections and planning priority road maintenance interventions. Furthermore, once the factors affecting the subsidence phenomenon have been determined, it is possible to design road maintenance interventions aimed at solving specific issues related to those conditioning factors.

CONCLUSION
Through the combination of SAR-based monitoring devices and ML techniques, we suggested a procedure for providing a road maintenance strategy that accounting for exogenous events of the infrastructure, such as subsidence effects. Qualitatively and quantitatively speaking, the research has shown that Regression Tree-based models can be considered a satisfactory tool for the prediction of PS-InSAR-based surface motion estimations. The outcomes of the model enable practitioners to assess infrastructural deficiencies by knowing the hydrological, geomorphometric, geomorphological, and social characteristics of the surrounding environment. Moreover, Tree-based models provide an estimation of the importance of predictor input factors. In the framework of this paper, the hydrological factors seem to be the most significant ones connected with subsidence effects. Further studies could be proposed for developing other types of Machine Learning models, making a comparison between them, and determining the most suitable one for this purpose. Also, the field of feature engineering and feature selection could be investigated to define a set of input factors in which the features are not correlated between each other or redundant. We can wrap up this study by asserting that the use of advanced network-scale monitoring devices in combination with a Machine Learning approach seems to be a promising work methodology for the improvement of processes on planning road maintenance by National Road Authorities. Dou, J., Yunus, A.P., Tien Bui, D., Merghadi, A., Sahana, M., Zhu, Z., Chen, C.-W., Khosravi, K., Yang, Y., Pham, B.T., 2019. Assessment of advanced random forest and decision tree algorithms for modeling rainfall-induced landslide susceptibility in the Izu-Oshima Volcanic Island, Japan. Sci. Total Environ. 662, 332-346.