Integration of an InSAR and ANN for Sinkhole Susceptibility Mapping: A Case Study from Kirikkale-Delice (Turkey)

: Suitable route determination for linear engineering structures is a fundamental problem in engineering geology. Rapid evaluation of alternative routes is essential, and novel approaches are indispensable. This study aims to integrate various InSAR (Interferometric Synthetic Aperture Radar) techniques for sinkhole susceptibility mapping in the Kirikkale-Delice Region of Turkey, in which sinkhole formations have been observed in evaporitic units and a high-speed train railway route has been planned. Nine months (2019–2020) of ground deformations were determined using data from the European Space Agency’s (ESA) Sentinel-1A/1B satellites. A sinkhole inventory was prepared manually using satellite optical imagery and employed in an ANN (Artiﬁcial Neural Network) model with topographic conditioning factors derived from InSAR digital elevation models (DEMs) and morphological lineaments. The results indicate that high deformation areas on the vertical displacement map and sinkhole-prone areas on the sinkhole susceptibility map (SSM) almost coincide. InSAR techniques are useful for long-term deformation monitoring and can be successfully associated in sinkhole susceptibility mapping using an ANN. Continuous monitoring is recommended for existing sinkholes and highly susceptible areas, and SSMs should be updated with new results. Up-to-date SSMs are crucial for the route selection, planning, and construction of important transportation elements, as well as settlement site selection, in such regions.


Introduction
Sinkholes are usually determined using traditional geomorphological methods and geospatial datasets such as aerial photographs, topographic maps, digital elevation models (DEMs), and field surveys. The effectiveness of these approaches has limitations in areas where the geomorphological features of sinkholes are suppressed by natural processes or anthropogenic activities [1]. In addition to these reasons, sinkhole mapping is often a difficult task when considering factors such as varying sizes, unpredictable nature, and a nonlinear rate. Although traditional geodetic techniques provide accurate results for sinkhole mapping and monitoring, they have significant drawbacks such as time-consuming fieldwork requirements, high costs in large areas, and the provision of deformation data via only a small number of ground-surveyed points [1]. For these reasons, it is generally not practical to apply traditional geodetic methods over large areas [2].
Satellite-based synthetic aperture radar (SAR) systems offer frequent revisit times, wide area coverage, and high precision and are frequently used in the determination and monitoring of deformations thanks to their advantages. In this context, interferometric SAR (InSAR) techniques and traditional geodetic techniques provide complementary ISPRS Int. J. Geo-Inf. 2021, 10, 119 2 of 16 information and have been used together in many studies [3][4][5]. Furthermore, InSAR techniques allow for the determination of deformation that has occurred in large areas for many years at centimeter level under suitable conditions (e.g., landslide and infrastructure monitoring, see [3][4][5][6][7][8][9][10][11]). On the other hand, InSAR techniques also have limitations such as data acquisition at fixed times, orbits, and geometric conditions, as well as the decrease in coherence values between different acquisitions depending on the temporal interval [12,13]. With the use of the InSAR and differential InSAR (DInSAR) techniques, highly precise deformation data can be obtained under suitable conditions.
Sinkhole research using SAR data has gained importance in the last two decades [1,2,[14][15][16]. One of the pioneering studies was conducted in 1996, aiming to determine the dissolution-induced subsidence on the Dead Sea Coast of Israel and Jordan [17,18]. Baer et al. [19] also investigated, in the same area, the rapid human-induced decrease in lake level, which caused continuous sinkhole formation, by using a total of 16 ERS-1 and ERS-2 SAR datasets acquired between 1992 and 1999. As a result, subsidence rates were determined in the range of 5 mm/year to 20 mm/year in large shallow depressions. Closson [20] determined subsidence rates as high as 24 mm/year in an area that collapsed as a result of salt evaporation in the year 2000 in the Southern Dead Sea Region of Jordan by using the ERS (European Remote-Sensing Satellite) images. For the same area, Closson et al. [21] argued that earthquake and fault activities played an important role in the occurrence of sinkhole formation in the Lynch Strait of the Dead Sea, based on InSAR-based deformation time series data acquired from the Japanese ALOS-2 (Advanced Land Observation Satellite-2) sensor. These studies have shown that InSAR-based techniques provide useful information for monitoring some structures such as dams and high-speed railways and estimating sinkhole subsidence that may occur by determining current displacements [1]. Paine et al. [22] used interferograms generated from ALOS-2 data to determine large sinkholes formed by salt dissolution in Texas, USA, and stated that the defined subsidence areas and deformation rate (up to 30 mm/year) were generally consistent with the data obtained from geodetic surveys. With the increasing temporal and spatial resolutions of SAR data and recent algorithmic developments in DInSAR processing and machine learning techniques such as PS-InSAR (Permanent Scatterers) and SBAS (Small BASeline Subset), their use in land subsidence and sinkhole susceptibility mapping has also increased (e.g., [2,13,[23][24][25][26][27][28]).
In this study, the possibility of integrating different InSAR techniques for sinkhole susceptibility mapping was investigated. Delice County in the Kirikkale Province of Turkey, where a high-speed train railway route was already planned, was selected as the study area. Sinkhole formations, which have been occasionally dissolving and collapsing, were observed in evaporitic units in this region. Ground deformations in the study area were determined by using data from the European Space Agency's (ESA) Sentinel-1A and 1B SAR satellites for a period of 9 months in 2019-2020. The morphological lineaments, the topographic conditioning factors obtained from the DEM, and a sinkhole inventory map were implemented in the susceptibility modelling performed by using the Artificial Neural Network (ANN) method. A total of 11 Sentinel-1 SAR datasets (two for DEM generation and nine for deformation monitoring) were obtained from the Copernicus Open Access Hub of the ESA and processed using the freely provided Sentinel Application Platform (SNAP) tool from the ESA [29]. The deformation map was produced from Sentinel-1 data. Finally, both results obtained from the ANN model and the deformation map were compared, and the integration potential of InSAR techniques with the sinkhole susceptibility mapping procedure was assessed. Additionally, based on the produced sinkhole susceptibility map (SSM), the suitability of the study area to the planned railway route was discussed.

Study Area and Sinkhole Inventory
The study area is located near Delice County in the Kirikkale Province of Turkey and in the north of Kirsehir Province. The location and extent of the study area, simplified geology, and sinkhole inventory map prepared by satellite images are presented as follows ( Figure 1). The area covers the Barakobasi, Halitli, Ocakbasi, Tavozu, and Tatlicak villages and is in the close vicinity of the Akcakent County of the Kirsehir Province and the Keskin County of the Kirikkale Province in the region.
The sinkhole formations observed in this study area develop within the Late Eocene-Early Miocene aged [30] Incik formation and the Sekili evaporite member of this formation. The lower parts of the Incik formation are medium-well coalesced, with fine-medium-thick horizontal bedded sandstones alternating with gypsum and anhydrite-rich mudstones, while the middle and upper layers are alternating with mudstones and increasingly crossbedded conglomerate and sandstone [30]. The red-brown, grey, light green, variegated white, medium-thick bedded evaporite, mudstone, and sandstone alternation, part of the Incik Formation, was named the Sekili evaporite member [31]. The Sekili evaporite member is generally composed of anhydrite and gypsum and includes sandstone and mudstone, while the top portion consists of mudstone and marl interlayers [30]. Donmez et al. [30] indicate that gypsum and anhydrides are white and yellowish-white, medium-thick bedded, generally pure and crystalline, and lumpy in places. Additionally, it was observed that gypsum and anhydrite get thinner laterally, and mudstones become thicker. The thickness of the Sekili member in the basin varies between 150 and 700 m [30]. Sandstones and mudstones are generally of thin-to-medium thickness and mottled in red-brown, light green, and green colors [30]. A general view of the Sekili evaporite member in the study area is also given ( Figure 2). The study area is located near Delice County in the Kirikkale Province of Turkey and in the north of Kirsehir Province. The location and extent of the study area, simplified geology, and sinkhole inventory map prepared by satellite images are presented as follows ( Figure 1). The area covers the Barakobasi, Halitli, Ocakbasi, Tavozu, and Tatlicak villages and is in the close vicinity of the Akcakent County of the Kirsehir Province and the Keskin County of the Kirikkale Province in the region.
The sinkhole formations observed in this study area develop within the Late Eocene-Early Miocene aged [30] Incik formation and the Sekili evaporite member of this formation. The lower parts of the Incik formation are medium-well coalesced, with fine-medium-thick horizontal bedded sandstones alternating with gypsum and anhydrite-rich mudstones, while the middle and upper layers are alternating with mudstones and increasingly cross-bedded conglomerate and sandstone [30]. The red-brown, grey, light green, variegated white, medium-thick bedded evaporite, mudstone, and sandstone alternation, part of the Incik Formation, was named the Sekili evaporite member [31]. The Sekili evaporite member is generally composed of anhydrite and gypsum and includes sandstone and mudstone, while the top portion consists of mudstone and marl interlayers [30]. Donmez et al. [30] indicate that gypsum and anhydrides are white and yellowishwhite, medium-thick bedded, generally pure and crystalline, and lumpy in places. Additionally, it was observed that gypsum and anhydrite get thinner laterally, and mudstones become thicker. The thickness of the Sekili member in the basin varies between 150 and 700 m [30]. Sandstones and mudstones are generally of thin-to-medium thickness and mottled in red-brown, light green, and green colors [30]. A general view of the Sekili evaporite member in the study area is also given ( Figure 2). A total of 25 sinkhole formations were mapped onto the evaporitic units in the study area (see Figure 1). The smallest sinkhole mapped is 169 m 2 , and the largest is 7265 m 2 ( Table 1). According to the classification of sinkholes given by Waltham et al. [32] and considering the dominant process behind the development of the sinkholes observed in A total of 25 sinkhole formations were mapped onto the evaporitic units in the study area (see Figure 1). The smallest sinkhole mapped is 169 m 2 , and the largest is 7265 m 2 ( Table 1). According to the classification of sinkholes given by Waltham et al. [32] and considering the dominant process behind the development of the sinkholes observed in the study area, sinkholes mapped during the inventory studies were classified as subsidence sinkholes (Figure 3a). Solution sinkholes were also observed in the region (Figure 3b). the study area, sinkholes mapped during the inventory studies were classified as subsidence sinkholes (Figure 3a). Solution sinkholes were also observed in the region ( Figure  3b).

Sentinel-1 Characteristics and the SAR Datasets
Sentinel-1 is a C-band radar instrument (with 5.4 GHz center frequency) on board the two-constellation satellites (Sentinel-1A&1B) operated by the European Commission and European Space Agency (ESA) in the Copernicus Global Earth Observation Project. They were launched on 3 April 2014 (1A) and 25 April 2016 (1B), respectively [33]. The satellites  the study area, sinkholes mapped during the inventory studies were classified as subsidence sinkholes (Figure 3a). Solution sinkholes were also observed in the region ( Figure  3b).

Sentinel-1 Characteristics and the SAR Datasets
Sentinel-1 is a C-band radar instrument (with 5.4 GHz center frequency) on board the two-constellation satellites (Sentinel-1A&1B) operated by the European Commission and European Space Agency (ESA) in the Copernicus Global Earth Observation Project. They were launched on 3 April 2014 (1A) and 25 April 2016 (1B), respectively [33]. The satellites

Sentinel-1 Characteristics and the SAR Datasets
Sentinel-1 is a C-band radar instrument (with 5.4 GHz center frequency) on board the two-constellation satellites (Sentinel-1A&1B) operated by the European Commission and European Space Agency (ESA) in the Copernicus Global Earth Observation Project. They were launched on 3 April 2014 (1A) and 25 April 2016 (1B), respectively [33]. The satellites supply four operational modes depending on the satellite position and scope of the observation modes, which are Stripmap (SM), Interferometric Wide (IW) Area, Extra Wide (EW), and Wave, and provide multi-polarization (VV, HH, VV + VH, and HH + HV) data that are combinations of vertical (V) and horizontal (H) polarization [34]. The temporal resolution is 12 days for each Sentinel-1 satellite, which results in 6 days for the combination [35]. Sentinel-1 data are free and publicly available through the Copernicus program. In this study, a total of 11 descending Sentinel images were used (i.e., six Sentinel-1A images and five Sentinel-1B images). The main characteristics (i.e., the satellite, acquisition date, data type and mode, pass direction, and usage purpose) of the SAR data employed in this study are shown in Table 2, together with the DS (Dataset ID) for all SAR data. For DEM production, a pair of SAR data (DS1 and DS2) with a temporal interval of six days, acquired in 2018, were used. A short temporal baseline was preferred to reduce atmospheric effects and temporal decomposition. However, it was also challenging to find a large perpendicular baseline between the image pairs since most of the perpendicular baselines between consecutive Sentinel-1 data are about 30 m, and they are primarily not designed for DEM generation [36]. The perpendicular baseline for the datasets used for DEM generation was 32 m, and the output DEM has a 10 m spatial resolution.
The ground deformation of the study area was analyzed using time series data with a 30-day temporal interval obtained between July 2019 and March 2020. In total, nine SAR datasets acquired at descending orbit (orbit number 167) in TOPSAR IW Swath Mode in Single Look Complex (SLC) format were used here. The VV (Vertical transmit and Vertical receive) polarization data were employed here due to the higher sensitivity of coherence information. The image pairs were selected to monitor monthly deformations on the surface. While selecting the data pairs, short perpendicular and temporal baselines were preferred. A short temporal baseline ensures high coherence between the pairs, since extensive changes on the surface due to land cover or human activities (e.g., crop growth) cause a decrease in consistency between the two acquisitions and may distort the results during the interferometric processing stage of the SAR data. The baseline information for the time series data used here were analyzed using ASF (Alaska Satellite Facility) Software [37] from the University of Alaska, and the results are shown graphically as follows ( Figure 4). The main methodological workflow consists of six main steps ( Figure 5). The DEM of the study area with a 10 m spatial resolution was produced using the SAR interferometry in SNAP software. One of the possible solutions for preparing a sinkhole inventory is the use of data extracted from publicly available incident reports [38]. Additionally, according to the methodology published by [39], identifying potential karst-related problems in remote sensing is also involved in the stages of site investigation for karst hazard assessment. In this study, the sinkhole inventory was produced by using satellite images provided on the Google Earth platform and checked through the short field studies performed in the study area. Sinkhole inventory data were used for training and test purposes in the ANN algorithm. Geological data were obtained from the WebGIS (Web Geographic Information System) portal of the General Directorate of Mineral Research and Exploration (GDMRE / MTA), Turkey, and the geological map of the region [30]. Topographic parameters, such as vertical distance to channel network (VDTCN), slope length, topographic wetness index (TWI), closed depressions, and convergence index, were produced from the DEM obtained from InSAR data using open-source SAGA (System for Automated Geoscientific Analyses) GIS software [40]. The theoretical background of the topographic parameters evaluated in SAGA GIS were provided by Olaya and Conrad [41] and Gruber and Pecham [42]. The ANN method was used as a predictor for sinkhole susceptibility using geological and topographic datasets together with sinkhole inventory data. The vertical deformation map produced from SAR time series data was used for verification and visual assessment of the susceptibility results. Further methodological details are presented in the following sub-sections.

Methodology
The main methodological workflow consists of six main steps ( Figure 5). The DEM of the study area with a 10 m spatial resolution was produced using the SAR interferometry in SNAP software. One of the possible solutions for preparing a sinkhole inventory is the use of data extracted from publicly available incident reports [38]. Additionally, according to the methodology published by [39], identifying potential karst-related problems in remote sensing is also involved in the stages of site investigation for karst hazard assessment. In this study, the sinkhole inventory was produced by using satellite images provided on the Google Earth platform and checked through the short field studies performed in the study area. Sinkhole inventory data were used for training and test purposes in the ANN algorithm. Geological data were obtained from the WebGIS (Web Geographic Information System) portal of the General Directorate of Mineral Research and Exploration (GDMRE / MTA), Turkey, and the geological map of the region [30]. Topographic parameters, such as vertical distance to channel network (VDTCN), slope length, topographic wetness index (TWI), closed depressions, and convergence index, were produced from the DEM obtained from InSAR data using open-source SAGA (System for Automated Geoscientific Analyses) GIS software [40]. The theoretical background of the topographic parameters evaluated in SAGA GIS were provided by Olaya and Conrad [41] and Gruber and Pecham [42]. The ANN method was used as a predictor for sinkhole susceptibility using geological and topographic datasets together with sinkhole inventory data. The vertical deformation map produced from SAR time series data was used for verification and visual assessment of the susceptibility results. Further methodological details are presented in the following sub-sections. lems in remote sensing is also involved in the stages of site investigation for karst hazard assessment. In this study, the sinkhole inventory was produced by using satellite images provided on the Google Earth platform and checked through the short field studies performed in the study area. Sinkhole inventory data were used for training and test purposes in the ANN algorithm. Geological data were obtained from the WebGIS (Web Geographic Information System) portal of the General Directorate of Mineral Research and Exploration (GDMRE / MTA), Turkey, and the geological map of the region [30]. Topographic parameters, such as vertical distance to channel network (VDTCN), slope length, topographic wetness index (TWI), closed depressions, and convergence index, were produced from the DEM obtained from InSAR data using open-source SAGA (System for Automated Geoscientific Analyses) GIS software [40]. The theoretical background of the topographic parameters evaluated in SAGA GIS were provided by Olaya and Conrad [41] and Gruber and Pecham [42]. The ANN method was used as a predictor for sinkhole susceptibility using geological and topographic datasets together with sinkhole inventory data. The vertical deformation map produced from SAR time series data was used for verification and visual assessment of the susceptibility results. Further methodological details are presented in the following sub-sections.

Processing of Sentinel-1 Datasets
The pre-processing steps applied to generate deformation and height information from Sentinel-1 SLC data are given in the workflow (Figure 6a). Data with suitable subswaths and bursts were selected from all S1 SLC data using the TOPSAR-Split in SNAP Tool, and precise orbit ephemeris was applied for accurate georeferencing before perform-

DEM Generation with the InSAR technique
Sinkhole inventory preparation using satellite images Geological and geomorphological data preparation ANN model construction Vertical deformation map, sinkhole susceptibility assessment

Processing of Sentinel-1 Datasets
The pre-processing steps applied to generate deformation and height information from Sentinel-1 SLC data are given in the workflow (Figure 6a). Data with suitable sub-swaths and bursts were selected from all S1 SLC data using the TOPSAR-Split in SNAP Tool, and precise orbit ephemeris was applied for accurate georeferencing before performing interferometric processing. Each dataset was stored as containing two split SLC images with Back Geocoding, according to the intended use shown in Table 2 (i.e., DEM generation or time series analysis of deformation). While DS-1 and DS-2 were used for DEM generation, DS3-4, DS4-5, DS5-6, DS6-7, DS7-8, DS8-9, DS9-10, and DS10-11 pairs were used to obtain monthly surface displacements. In order to produce height and displacement information from SAR images, ESD (Enhanced Spectral Diversity) was applied, and interferograms were generated for each dataset. With TOPSAR-Deburst, bursts in the selected sub-swaths were resampled in a common grid in range and azimuth directions. All bursts in a swath were merged into a single image with TOPSAR-Deburst, and topographic phase contributions were removed. Finally, for the pre-processing step, Goldstein phase filtering [43] and phase unwrapping processes were applied to correlate the interferometric phase with the topographic height.  (DS3-4, DS4-5, DS5-6, DS6-7, DS7-8,   DS8-9, DS9-10, DS10-11 DS1-2, DS3-4, DS5-6, DS6-7, DS7-8,   DS8-9, DS9-10, DS10-11  In the second step (Figure 6b), the unwrapped phase of each dataset was converted into elevation information (Phase to Elevation). DEM data were terrain-corrected and orthorectified using the Range Doppler Terrain Correction algorithm and 1 Arc-Second SRTM DEM. The DEM was later used for extracting the topographic derivative to utilize in the ANN method. The accuracy of the InSAR DEM is affected by imaging conditions such as atmospheric conditions, noise, and surface deformation, as stated by Braun [44]. In the coherence and interferograms for the interferometric pair, high coherence and low noise were observed. Even though the InSAR DEM accuracy is possibly equivalent to the accuracy of SRTM DEM or slightly better [45], the use of InSAR was favorable in the study since it represented the topographical state for the study date, and a higher spatial resolution (10 m) could be obtained.
In order to produce a vertical displacement map from the time series data (Figure 6c), unwrapped phases were transformed into displacement data and orthorectified with a process similar to DEM production. Capable of detecting only 1-D deformation in the LOS (line of sight) direction of the satellite, the DInSAR technique cannot be recommended for the determination of ground deformation with long time series based on two images [46,47]. Considering the relatively flat topography and the expectation of mainly vertical displacement in sinkhole-prone areas, further processing was carried out to obtain the vertical deformation component, as shown in Equation (1).
Division of LOS displacement to the cosine of the radar wave incident angle (cosϑ) was found as an easily applicable approach to determine vertical deformation in the literature [48,49]. The method ignores the contribution of horizontal deformation to LOS displacement [50]. Afterwards, the vertical displacement was georeferenced in the geographic coordinate system with WGS84 (World Geodetic System 1984) datum for further analysis. Vertical displacement maps were produced by this method for each data set. The steps of this process are shown in Figure 6c. Eight vertical displacement maps, each containing information on vertical displacement within a one-month period, were added cumulatively to determine the total vertical change over nine months (Figure 6c,d). A feature set with six topographic parameters, such as Convergence Index (CI), Closed Depressions (CD), TWI, Slope Length Factor (LS), and the VDTCN, were utilized in the ANN model together with the training data, including morphological lineaments and sinkhole inventory, to produce the final SSM (Figure 6d). The cumulative vertical displacement map was also employed for comparison analysis of the SSM results in sinkhole-prone areas.

Sinkhole Susceptibility Assessment with the ANN Model
Prediction performances in susceptibility assessments of natural hazards depend on the selected method, data quality, and field characteristics. Recent studies have shown that knowledge-based artificial intelligence methods provide a significant performance increase not only in susceptibility evaluations of different natural hazards but also in disaster management applications [51].
The ANN has become a popular method used by researchers in the last two decades, particularly in susceptibility mapping related to mass movements [52][53][54][55][56]. A neural network can be considered as a semantic extraction model based on the working of principles of the human brain [57]. Considering its structure, the human brain can be expressed as a highly complex, non-linear system with parallel-positioned information processors. In this case, an ANN can be described as a simplified reproduction of a mirror of the complex system. Hecht-Nielsen [58] and Schalkoff [59] describe ANNs as structures consisting of simple processing elements called neurons, which are densely connected to each other and have the capacity to perform parallel processing for data processing and information presentation.
The first step in the process of developing ANN models is to split the existing data into two separate groups as training and testing datasets. The training dataset used in the model construction is expected to involve the necessary number of parameters to solve the problem in question. In other words, in the ANN model construction, it is desired that the training data sets should have sufficient numbers of presence "1" and absence "0" data to successfully predict the susceptibility of mass movements. On the other hand, the testing dataset should be independent of the model construction process. The main purpose of using the testing dataset is to evaluate the performance and generalization capacity of the established ANN model by using the dataset that is not included in the model construction process. There is no exact mathematical approach for the determination of an appropriate structure of the model that should be implemented in training [60]. The second step in the ANN model construction is the assignment of initial values for the weights and thresholds that provide connections within the neural network structure. Negnevitsky [57] states that neural networks can be successfully developed even though the resulting model may consist of different weights and thresholds depending on different initial values.
To develop the ANN model, a program was written in Python by using the scikit-learn libraries. In the procedure, data are first normalized between the closed interval of [0, 1]. Afterwards, the data matrix is separated, depending on the presence and absence information of target. To achieve a ratio of 50% between the presence and absence information in the training data set, random sampling is performed from the absence ("0") by considering the number of rows in the presence ("1"). In this study, a total of six topographic parameters and one lineament parameter were used as independent features, and sinkhole formations were used as the dependent target in the construction of the ANN model. Bedrock type, proximity to geological boundaries and structures, topography, existing sinkholes, and drainage history can be considered to be the main conditioning factors in sinkhole formation [32]. In this study, unique, bedrock-type, and Late Eocene-Early Miocene-aged evaporites were evaluated. The topographic parameters, convergence index, closed depressions, topographic wetness index, and vertical distance to the channel network are completely related to the drainage history of the region. The topographic parameter of the slope length factor indicates the spatially sediment transport capacity in the region [61]. The areas with low slope length factor values, which mean sediment transport capacity is low, are considered to be more prone to sinkhole formation. Morphological lineaments, mapped using satellite imageries, were considered to be related with the geological boundaries and structures. Eventually, the data matrix consisted of topographic parameters, and the lineaments were produced. In the matrix, the dependent target was the field represented by existing sinkholes, and the independent features were the parameters that controlled sinkhole formation. Using these features and the target, the ANN model was developed. The training dataset is composed of a matrix including 258 grid cells with an equal portion of presence and absence information of the existing sinkholes. Ten percent of the training data set was also selected at random as the testing dataset. Additionally, in the construction of the model with the training dataset, k-fold cross validation was applied. In the analysis performed using 10 folds, the most successful model was obtained in the learning process, where the learning coefficient was 0.0001. Three hidden layers containing 100 neurons and the rectified linear unit (ReLU) transfer function were used.

Results
The DEM of the study area obtained from InSAR processing is shown together with the spatial distribution of the convergence index in Figure 7a,b. The convergence index is similar to the parameter plan curvature but provides smoother results. The positive and negative values reflect convergent and divergent flow conditions, respectively. When the convergence map is examined, it is revealed that the existing sinkholes are mainly observed in convergent flow conditions. Some other topographic parameters extracted from the DEM are also presented (Figure 7c-f). The closed depressions represent the closed contours in the form of a hollow (Figure 7c). The TWI defines the effect of topography on hydrological processes. The spatial distribution of the index is given in Figure 7d. When the TWI map was examined, it was observed that the existing sinkholes tended to be located on the average values of that parameter. The slope length factor describes the effect of topography on soil erosion ( Figure 7e). As mentioned previously, the parameter indicates spatially sediment transport capacity in the study area, and it is expected that the areas with low slope length factor values are more prone to sinkhole formation. According to the resultant slope length factor map, this expectation is verified: sinkholes are almost exclusively observed on the lower values of the parameter in the study area. The last topographic parameter derived from the DEM was the VDTCN (Figure 7f). This parameter indicates the relative position on the slope with respect to the channel basement. When the resultant VDTCN map was examined, it was revealed that existing sinkholes were observed at the lower values of those parameters. The result points out that surface waters have a critical effect on sinkhole formation in the region.   By considering the best parameters and structure explained in Section 3.2, the final ANN model was trained and tested. The Root Mean Square Error (RMSE) obtained from the model was 0.451. Considering the standard deviation value of the testing data, which is 0.476, and the calculated RMSE, the model was found to be successful. As a result, the constructed ANN model was run by using the data matrix of the study area, which involved 158,902 grid cells, and the resultant SSM was obtained ( Figure 9).   Sinkhole susceptibility was presented in the form of probability ranging between 0 and 1. The distribution of the recent sinkholes indicates the prediction performance of the ANN model. The blue regions on the SSM indicate a high susceptibility to sinkhole formation.
The negative and positive displacements on the ground were produced from InSAR time series data. Comparative evaluation of areas with high susceptibility to sinkhole formation and areas with surface deformations determined by the interferometry was carried out. When the sinkhole inventory, deformation, and SSMs were superimposed, the existing sinkholes, negative deformations in the interferometer, and areas with high sinkhole susceptibility values on the susceptibility map indicated the most critical conditions in terms of railway route selection ( Figure 10). Vertical deformations were observed in the existing sinkholes and in areas with high sinkhole susceptibility values. Based on the results, high deformation areas on the vertical displacement map and areas that are highly prone to sinkhole formation in the SSM almost coincide.

Discussion and Conclusions
In areas where evaporitic rock units are observed, sinkhole hazard assessment and prevention studies can be carried out in four stages: (i) inventory of existing sinkholes or developing depressions; (ii) field investigations such as drilling, geophysics, and preparation of high-resolution terrain models; (iii) prediction models; and (iv) mitigation [62]. In mitigation applications, the most effective approach is to avoid these subsidence areas or areas prone to sinkhole formation [63]. In this study, an approach based on sinkhole inventory, sinkhole formation susceptibility, and currently observed vertical deformation was proposed. The planned Kirsehir-Kayseri railway has to pass through evaporitic rock units. However, subsidence sinkholes were observed on the route due to these evaporitic units in the region. The sinkhole inventory was first delineated using recent optical satellite images. In this context, lithology was kept constant by focusing on a single formation Figure 10. Superimposition of the maps, sinkhole inventory, observed vertical deformation, and sinkhole susceptibility, indicating the most critical conditions in terms of railway route selection in the study area.

Discussion and Conclusions
In areas where evaporitic rock units are observed, sinkhole hazard assessment and prevention studies can be carried out in four stages: (i) inventory of existing sinkholes or developing depressions; (ii) field investigations such as drilling, geophysics, and preparation of high-resolution terrain models; (iii) prediction models; and (iv) mitigation [62]. In mitigation applications, the most effective approach is to avoid these subsidence areas or areas prone to sinkhole formation [63]. In this study, an approach based on sinkhole inventory, sinkhole formation susceptibility, and currently observed vertical deformation was proposed. The planned Kirsehir-Kayseri railway has to pass through evaporitic rock units. However, subsidence sinkholes were observed on the route due to these evaporitic units in the region. The sinkhole inventory was first delineated using recent optical satellite images. In this context, lithology was kept constant by focusing on a single formation (evaporitic rock units) throughout the study. Accordingly, geological and hydrological parameters that can be effective in the development of existing sinkholes were investigated. In order to evaluate the lineaments, a morphological lineament map was produced. Topographic parameters were determined from the DEM produced from Sentinel-1 SAR data. By using the DEM, the spatial distributions of the convergence index, the closed depressions, TWI, slope length factor, and VDTCN in the study area were determined. In addition, the distance of each grid cell to the morphological lineament was calculated.
As stated by Mohammadi et al. [64], Sentinel-1 data can be successfully used for DEM generation due to its high spatial resolution and C-band sensitivity to the surface structure. The accuracy of the Sentinel-1 DEM was found to be slightly worse than ALOS global DEM and SRTM [64], similar to the ASTER global DEM [65], and slightly better than the SRTM DEM [45]. It was also found that the baseline of the Sentinel-1 pair does not significantly affect the final accuracy of DEMs [65]. The main land cover types in the region are grasslands and agricultural fields, and the main products are wheat and maize. One study [66] has found that the C-band data of Sentinel-1 in VV polarization are able to penetrate the maize canopy, even when it is well developed. Thus, the vertical deformation obtained in the study can be associated with changes on the ground. Although this band is sensitive to well-developed grassland and wheat based on Normalized Difference Vegetation Index (NDVI) analysis results [66], considering the selected monitoring period (23 July 2019-19 March 2020), under continental climatic conditions, it can be said that the vegetation has low canopy development and thus that sufficient penetration to the ground is possible. However, it should be noted that land cover and land use changes may reflect other factors such as vegetation growth or other constructions in some areas, which may cause errors in the deformation map.
In this study, interferometric SAR data from Sentinel-1 satellites were utilized for the DEM generation and temporal monitoring of surface deformations, particularly in the vertical direction, in order to evaluate sinkhole susceptibility in part of the planned high-speed railway route in Turkey. The topographic features obtained from the DEM, the inventory of existing sinkholes, and the evaporitic units of the study area were employed in the ANN method for producing an SSM. The results were compared and validated with the vertical deformation map obtained in the study. The proposed approach was found to be successful for sinkhole susceptibility assessment by using InSAR data as input.
The results obtained within the scope of this study indicate that the proposed method can be used effectively in determining the routes of linear engineering structures. As a consequence, the railway route planned in critical conditions needs to be revised before construction. In future studies, deformation monitoring will continue on both existing sinkholes and areas with high sinkhole susceptibilities, and the sinkhole susceptibility model will be updated depending on new monitoring results.