INTEGRATION OF GIS AND ADVANCED REMOTE SENSING TECHNIQUES FOR LANDSLIDE HAZARD ASSESSMENT: A CASE STUDY OF NORTHWEST SYRIA

Landslide susceptibility and hazard mapping has been developed providing remarkable results through the integration of geographic information system (GIS) and remote sensing. In this regard, some approaches have considered the use of Sentinel-1 data and timeseries interferometric synthetic aperture radar (InSAR) techniques, such as differential InSAR (D-InSAR) and persistent scatterers interferometric (PSI), for providing precise information about total amount and velocity of ground-surface deformations and landslides within a specific area during a specific time period which is important for disaster management’s planning process. In this paper, artificial neural network (ANN) was used as a statistical analysis method for landslide susceptibility mapping in Northwest Syria using multi-layer perceptron (MLP) neural network on a training dataset of one dependent variable (landslide or non-landslide) and nine independent variables (slope, aspect, curvature, land cover, NDVI, lithology, distance from faults, distance from road, distance from stream networks). The resulting map of landslide susceptibility was validated using area under curve (AUC) analysis using a testing dataset which showed 90.28% of AUC value. Then, landslide susceptibility map was reclassified into high-moderate-low classes and integrated with intensity map of mean velocity of ground-surface deformations during the time period form (16 October 2018) until (21 March 2019) by using a landslide hazard matrix in a GIS environment in order to get landslide hazard map of the study area for that time period. The result shows that around 44.4%, 52.9% and 2.5% of total study area was classified as a high, moderate and low hazard zone of landslide, respectively. * Corresponding author


INTRODUCTION
Landslide hazard investigation, assessment, monitoring and mapping processes become increasingly important in most mountainous and hilly regions around the world where landslides are a major hazard (kervyn et al., 2015) and (Hammad et al., 2019). So, landslide studies in those places are important to evaluate the relative contributions of all elements and factors involved in landslide process within an area and to identify the possible locations that might be affected by landslides in order to reduce risk on people in those areas and help local authorities in the process of future regional planning.
Previous attempts of using Earth observation data for landslides investigation and mapping studies were typically based on field work and visual interpretation of stereoscopic aerial photography which are relatively time-consuming and not always easy to be carried out (Guzzetti et al., 2012). Later, with the start of the satellite data era, processing and analysing of optical satellite data has been applied to detect any relief changes between different dates. Furthermore, landslide investigation has been developed providing remarkable results, especially through the integration of optical remote sensing data and geographic information system (GIS) using different statistical analysis models in order to investigate and map susceptibility and hazard of landslides (Biswajeet, Saro, 2007). Recently, it has been shown that synthetic aperture radar (SAR) satellite data can be used as a complementary data source providing useful and precise information about ground-surface deformation and landslides (Hammad et al., 2018). In this direction, some approaches have already considered the use of time series interferometric SAR (InSAR) techniques for providing information about stability of areas suffering from ground-surface deformations and landslides by analysing velocity and amount of these ground-surface deformations (Hammad et al., 2018). The widely used approach applying time-series interferometry for ground-surface deformations detection depends mainly on persistent scatterers interferometric (PSI) technique (Hooper, 2006). Moreover, artificial neural network (ANN) analysis which can be applied to landslide susceptibility is a computational mechanism that can acquire, represent, and compute a mapping from one multivariate space of information to another, given a dataset representing that mapping (Garrett, 1994). For many years, artificial neural networks have been employed in a wide range of classification applications in Earth sciences (Van Leeuwen, 2012). Also, some studies have also applied artificial neural networks as a statistical analysis method for landslide susceptibility mapping (Ermini et al., 2005. New approaches using GIS and advanced remote sensing techniques are required for all researches in field of landslide investigation to take advantage of all freely available data and related techniques. Therefore, this paper attempts to integrate free interferometric data with help of GIS in order to investigate and determine locations of possible future landslides in northwest Syria and to assess landslide hazard during the study period.

Aim and objectives
Road networks in the study area are susceptible to be damaged annually due to unforeseen landslides after intense rainfall events. These unforeseen landslides threaten the lives and the properties of people who live there (Hammad et al., 2019).
The main direct losses are costs of removing debris from the roads and costs of roads construction, especially the international road between Latakia city, the main city in the Syrian coastal area, and the Turkish border near Kassab city in the north (Figure 1). While the main indirect losses are the disruption of the economic and social activities which can frustrate local and regional development. Some pervious geological and geotechnical studies were carried out after significant landslide events, but these studies are more descriptive geotechnical reports and are limited to specific sites and very small isolated spatial extents. Most of these geotechnical reports indicated that the ribbon radiolarite deposits were the main lithology at the sites where landslides occurred, and all landslides at those sites were rotational landslides. However, the study area lacks wide extent and continuous spatial studies of landslide susceptibility and landslide hazard which can be updated continuously in order to take the necessary precautions in possible landslides locations and minimize the risk of these landslides.
The overall aim of this research work is to get benefits from the integration of geographic information system (GIS) and remote sensing data, taking into consideration applying interferometric techniques in order to guide the landslides investigation process after defining the optimal landslide susceptibility map by using landslide statistical analysis in the study area, in order to achieve the landslide hazard map in northwest Syria as a final result. This will be achieved by meeting these objectives: -Assessing the suitability of using Sentinel-1 data and SAR interferometric techniques for landslide investigation.
-Evaluating the appropriateness of using SAR interferometric techniques results with landslide susceptibility map for the aim of landslide hazard mapping.

Study area
The study area extends over an area of 269 km 2 from 35.75 • in the south to 35.94 • at the Turkish border in the north, and from 36.00 • in the east to 35.80 • on the Mediterranean Sea in the west, which corresponds to the extension of Kassab topographic map 1:25,000 which was printed in 1990 after cartographic field works in 1984 based on aerial photos taken in 1983.
The elevation in the study area ranges between 0 and 1130 metres above sea level. It should be mentioned that adjacent to the study area and to the north is the highest mountain in the region, Jebel Al-Aqra mountain, rising from a narrow coastal plain to reach 1717 metres ( Figure 2). The highest monthly average of precipitation over study area is in December and January with around 125 mm, while maximum monthly precipitation reached 379 mm in January 2012. Also, it should be mentioned that many extreme precipitation events occur in the study area where the maximum daily precipitation reached around 70 mm.
In general, the study area consists of the Badrousieh basin, most of the Al-Bassit basin, a very small part of Alkabir Alshamali basin, and almost the half of Wadi Qandil basin within which the Balloran dam is located. The study area also has a good road network linking all the small towns and villages in the area with the main city of Latakia and with the Turkish border gate in Kassab. The main sub-districts in the study area are Qastal Maaf sub-district which consists of 19 localities with a collective population of 16,784 in 2004 official census, and Kassab sub-district which has a total population of 2,500 along with the surrounding villages.
According to many geotechnical reports carried out by the General Establishment of Geology and Mineral Resources in Syria for several landslide events in the study area, heavy rainstorms and extreme precipitation events are the main trigger for all of those landslides. In recent years, these events have increased, and different landslides have occurred at various locations in the study area and along its road network.
From geological aspect, the study area is characterized by spread of the Baer-Bassit ophiolitic complex deposits in addition to different Mesozoic and Cenozoic sedimentary rocks. The ophiolitic complex outcrop which represents fragments of oceanic crust and lithospheric mantle exposed on land by tectonic processes, is dominated by two massifs, Baer in the east and Bassit in the west near to the coast. The Baer massif rocks is relatively structurally intact, while the Bassit massif rocks are overthrusted by thin imbricate thrust sheets of pillow lavas (Al-Riyami et al., 2000).

DATA
The Kassab topographic map from the 1980s was used in order to represent the topographic characteristics of the study area similarly to that which were before the occurrence of all landslides participating in this study. This is very significant to know the real effects of the main topographic factors on the landslide occurrence and to get reliable results.
Optical multispectral data of Landsat-5 Thematic Mapper (TM) from 22 September 1984 with 0% cloud cover and 30 m resolution were used in order to produce land cover map and the normalized difference vegetation index (NDVI) map of the study area to reflect also the vegetation cover characteristics similarly to that which were before the occurrence of the landslides taken into consideration in this study. Furthermore, since lithology is constant in general during the study period, and in order to get benefits of high resolution optical multispectral data of the European Space Agency (ESA), Sentinel-2A image from 08 September 2017 with 0% cloud cover and 10 -20 m resolution in VNIR and SWIR range was used applying principal component analysis and band ratio techniques to prepare the lithology map with the help of the geological map which was also used for preparing the fault map.
Additionally, radar satellite data from 14 Sentinel-1B single look complex (SLC) images during the time period from 16 October 2018 until 21 March 2019 in interferometric wide (IW) swath mode and vertical-vertical (VV) polarization were used to produce the mean velocity map of the ground-surface deformation during this period of time using persistent scatterers interferometric (PSI) technique. The results were combined with the landslide susceptibility map to get the final landslide hazard map for the specified time period.

Landslide susceptibility mapping
Considering the local characteristics and different criteria such as type of landslides, mapping unit and availability of data taking into account similar landslide susceptibility studies (Guzzetti, 2005) and (Lee et al., 2012), nine main causative factors controlling the occurrence of landslides were considered. These factors are slope angle, slope aspect, terrain curvature, land cover, normalized differential vegetation index (NDVI), lithology, distance to lineaments, distance to roads and distance to stream networks. Note that within the relatively small study area, precipitation is likely consistent, therefore the use of the data for this factor within the analysis processes has been excluded in this research.
Each raster layer of the nine causative factors was reclassified into number of classes depending on its own criteria in order to be used in the ANN analysis method.
Basically, the first step of landslide susceptibility mapping is preparing a database of all previous landslides that happened in a given area during a specified period of time. Depending on this database, landslide inventory map of that area during that period of time will be prepared in order to be used with related causative factor maps through statistical analysis methods for defining the relationship between landslide occurrences and those related factors. Actually, landslide and non-landslide dataset was obtained using the information from the available technical reports and from local people regarding previous landslide events, and also from comparison between available older and recent optical satellite data and with the help of the ancillary information available through the historical imagery tool in Google Earth.
Then, the landslide and non-landslide dataset was randomly divided into two datasets; a 70% training dataset of landslides to be used in preparing the data from all causative factors during the training and testing processes of the artificial neural network (ANN) statistical analysis method, and a 30% testing dataset of landslides to be used for the validation process of the statistical analysis method (Figure 3). The pixel size for all maps was set to 12.5 m, same as the recommended grid resolution of the DEM extracted from the topographic map 1:25,000 used in this research (Hengl, 2006).
For the ANN analysis in this research, a feed forward multilayer perceptron (MLP) neural network architecture using a back-propagation (BP) learning algorithm was applied on the training dataset of landslides and non-landslides pixels and their related causative factors. The training dataset was prepared using the Raster-To-Point and Extract-Multi-Values-To-Points tools in ArcGIS in order to extract the data from all layers to a table in which each row stores one pixel's data of the dependent variable (landslide or non-landslide) and the nine independent variables (the causative factors). This table was imported to Matlab ® as a numeric matrix in order to set the inputs and the outputs (targets) and start the neural network training phase.
In fact, the number of hidden layers in the training phase depends mainly on the dependent and the independent variables and the complexity of the analysis. In this research, the structure of 9 input neurons, one hidden layer with 5 neurons, and 2 neurons in the output layer was selected for the architecture of the neural network ( Figure 4). The dataset was divided randomly into 70% for training, 15% for validation and 15% for testing. It should be mentioned that all the input training data used in the training phase were numeric, and these numeric values were not ordinal data but nominal data and they denote the different classes of the causative factors. The susceptibility result of ANN statistical analysis method was validated depending on the (30%) landslide testing dataset using the area under the curve (AUC) analysis which allows to describe the statistical analysis ability of correctly predicting occurrence of landslides. So, the landslide susceptibility map that resulted from the ANN analysis was reclassified into 100 class equally and sorted in descending order (it means classes having high susceptible to landslide should appear at top).
Then, the resulting map was used with the testing dataset to calculate the pixel-based landslide areas in each class. The percentage of landslides pixels in each landslide susceptibility index class of the total landslide pixels in all the test dataset were calculated. Then, the cumulative of landslide area was calculated in each class from this percentage. Finally, the AUC value in each class was calculated using the following equation (Vakhshoori, Zare, 2018): Where (X{i+1} -X{i}) = the difference value between two consecutive susceptibility classes; and (Y{i+1} + Y{i}), is the sum value of the cumulative percentage of landslides in these two relevant consecutive susceptibility classes.
The sum of all AUC values was calculated and divided by (100) to get the final value of AUC related to ANN analysis.

Landslide hazard mapping
Since a systematic record of landslides during the past decades in the study area is not available for landslide recurrence analysis and since the study area is small to be divided into different zones according to precipitation which is the main landslide trigger in the study area, landslide hazard mapping was carried out in this research depending on landslide hazard assessment using a landslide hazard matrix ( Figure 5) (Hürlimann et al., 2008) and (Lu et al., 2014). Figure 5. The high-moderate-low landslide hazard assessment using a landslide hazard matrix (after Hürlimann et al., 2008).
Basically, this matrix focuses on two parameters which are probability of landslide occurrence (landslide susceptibility), and intensity of ground-surface deformations. In this research, these two parameters were correlated using this matrix in order to represent the landslide hazard as low-moderate-high zones.
In general, the intensity of ground-surface deformations can be defined using various parameters including volume or velocity (Li et al., 2010). In this research, the intensity of ground-surface deformations was defined using the mean velocities of these deformations in the satellite line-of-sight (LOS) direction during the study time period which was chosen from October 2018 till March 2019 to cover one main precipitation season in the study area.
To produce the mean velocities map of ground-surface deformations, PSI technique was applied using radar satellite data of single look complex (SLC) imageries from the descending orbit, Sentinel-1B, in interferometry wide swath mode ( After preparing SAR data in SNAP as it is shown in figure (6), the resulted products were processed by using StaMPS in Matlab under Linux in order to perform the PS pixels' estimation process and calculate the mean velocity values of ground-surface deformations in each of these PS pixels averaged over the total study period in the satellite LOS direction by unwrapping all phases in the different produced interferograms. Figure 6. Flow chart of preparing SAR data in SNAP for PSI.
Eight steps were performed by StaMPS on the products resulted from the processing of SAR data. The first step was loading the data for initial candidates' selection process of PS pixels using an amplitude dispersion threshold value which was set as 0.4 (Hooper et al., 2013). In the second step which is an iterative step, the estimation process of the phase noise value for each candidate pixel in every interferogram was performed. Then, PS pixels and non-PS pixels were sorted on the basis of their noise characteristics. In the fourth step, the PS Pixels selected in the previous step were weeded, dropping those that are due to signal contribution from neighbouring ground resolution elements. After that, the wrapped phase of the selected pixels is corrected for spatially-uncorrelated look angle error. In the sixth step, filtering and unwrapping phase processes were performed using related algorithms. Then, spatially-correlated look angle error was calculated, and the master atmosphere and orbit error (AOE) phase was also estimated simultaneously. Finally, in the eighth step, the mean velocity values of groundsurface deformations were calculated in each of the PS pixels averaged over the total study period in the satellite LOS direction. As the result, a multipoint dataset of mean velocity values in millimetres was produced representing the deformations during the whole study period, 156 days, and the resulted data was exported to an ESRI shapefile in a point feature format in order to be processed after that in ArcGIS. In order to obtain a continuous spatial extent for the mean velocity values of ground-surface deformations and to create the highmoderate-low zonation areas map required with the total ground-surface deformation map for ranking the intensity of the ground-surface deformations in the study area during the study period, interpolation process of the velocity values was performed through inverse distance weighted (IDW) interpolation method since the multipoint dataset of the mean velocity values is interconnected dataset in the spatiotemporal domain (Liu et al., 2019).

RESULTS AND DISCUSSION
A landslide inventory map of all landslides that happened in the study area since 22 September 1984 was prepared in ArcGIS depending on all available data and information. Also, the nonlandslide dataset with the same size of the landslide dataset was also prepared depending on the related geotechnical reports and the knowledge gained through field works there. As a result, a dataset of 57 landslide polygons (Figure 7-a) and 57 nonlandslide polygons (Figure 7-b) were prepared and mapped.
Then, this dataset was randomly subdivided into two datasets, 70% training landslide dataset, to be used in preparing the data from all causative factors during the training process of the artificial neural network (ANN) method, and 30% testing landslide dataset, to be used later in the validation process.    All maps of the nine causative factors were prepared in ArcGIS and converted to a raster format, and clipped to the same spatial extent. Moreover, all the map layers were resampled using the same geometric projection WGS_1984_UTM_Zone_36N in order to be unified and have the same pixel size, 12.5 m. As a result, the study area has 1693340 pixels, and training dataset of landslides in the study area contains 1353 pixels.
After preparing the causative factors in ArcGIS in a (2706×10)  (11). Figure 11. The landslide susceptibility map of the study area using ANN analysis.
The validation process of the ANN analysis results showed an AUC value equal to 90.28% as it is shown in figure (12). March 2019 over the mountains in the north and around Balloran dam lake in the south due to the extreme precipitation events before these dates ( Figure 13) same as the velocity values of the persistent scatterers detected by PSI technique during the study period which show several extreme negative values in the northern mountains and the area around Balloran dam lake in the south (Figure 14).  The resulting data was exported to ArcGIS in order to be processed and interpolated using inverse distance weighted (IDW) interpolation method. As the result, a raster map representing the mean velocity map of the ground-surface deformation during the study period from 16 October 2018 until 21 March 2019 was produced as it is shown in figure (15). This map shows values between -47.70 mm/year and + 101.64 mm/year and also show more clearly the specific locations of extreme ground-surface deformation values which refer to the high landslide intensity areas. Figure 15. The interpolated mean velocity map of the groundsurface deformation produced using PSI technique.
Then, the mean velocities in the satellite LOS direction was reclassified into high-moderate-low classes in order to produce the intensity map of ground-surface deformations ( Figure 16) to be used with the landslide susceptibility map of the study area in the landslide hazard matrix to get the final map of the landslide hazard during that period of time.

CONCLUSION
This research demonstrates the potential and capability of radar satellite data and In-SAR time-series technique to investigate and monitor the ground-surface deformation, and also to measure millimetre velocity values in the line-of-sight direction over time using freely available data and software.
The integration of ANN analysis and PSI technique with the help of GIS can be considered as an attractive method for geological hazards tasks giving important results which can be updated regularly in order to be used in the disaster management's planning process.