Geoinformatics-based assessment of land deformation and damage zonation for Gorkha earthquake, 2015, using SAR interferometry and ANN approach

The Gorkha earthquake 2015 was one of the largest disastrous events that occurred in the main Himalayan thrust (MHT) region with epicenter at Gorkha region and magnitude of 7.8 which caused severe causality to life as well as property. The spatial statistics on vertical displacement and extent of damage zone is still too scarce to provide strong evidences of its hazard potential. In the present study, quantitative assessment on surface deformation has been carried out to compare land displacement of two different regions of Nepal (the Central Nepal and the Eastern Nepal), which are located at different distances from the epicenter, using InSAR technique on post- and pre-earthquake images from Sentinel 1A SLC product. The Central Nepal experiences an upliftment of 1.1 m and land subsidence of − 0.61 m, whereas for the Eastern Nepal the estimated upliftment and subsidence were 1.0 m and − 0.33 m, respectively. Further a regional earthquake-prone zone map was generated using the historical earthquake epicenter data and geographic information system (GIS) to understand the major vulnerability zones in the area. A total of 564 earthquake events were reported by USGS in Nepal region during 2000–2019, of which 476 (84.39%) were of magnitude greater than 4 on Richter scale and 376 events (66%) occurred at depth greater than 15 km. The damage assessment was done using machine learning (artificial neural network) back-propagation model in which the satellite imagery retrieved from the optical satellite Landsat 8 OLI sensor and digital elevation model was used to map slope, aspect, relief, drainage and lineament to be used as input layers to generate damage proxy map. The result obtained from ANN illustrated that despite being located comparatively at more distance from epicenter, the Eastern Nepal exhibited more damage-prone area (587 sq. km) in comparison with Central Nepal with 457sq. km damage prone in similar zone. Central Nepal evidences more damage-prone areas over compact build-up in contrast to Eastern Nepal, making greater risk potential in urban areas of Central Nepal during earthquake activity.


Introduction
Earthquake is defined as a tectonic or volcanic phenomenon that represents the movement of rock and generates shaking or trembling of the Earth's surface [4]. It is considered as one of the most severe natural hazards which can neither be predicted nor controlled. Earthquakes occur on faults, representing a clear link between seismology and structural geology [8]. The major earthquakes occur mainly in belts coinciding with the margins of tectonic plates. An inter-plate earthquake is caused mainly due to three types of faults, namely Normal Fault, Reverse (Thrust) Fault and Strike slip Fault. An earthquake's point of initial fissure is called its focus or hypocenter, whereas the point at the ground level beeline above the hypocenter is called the epicenter of the earthquake which generally experiences the more damage.
Many studies have been carried out to estimate the surface deformation, assess the damage caused due to earthquake and to identify the earthquake vulnerable areas using different remote sensing and GIS techniques. The most important and well-known technology for surface deformation estimation is SAR Interferometry (InSAR). Furuya [10] defined InSAR as a radar technique to image topography and ground displacements using phase difference of two or more temporal SAR images. In principle, different distributions of the receiving/ transmitting antennas define distinctive interferometric configurations [16]. Massonnet et al. [13] demonstrated the potential of InSAR for seismological application by mapping the co-deformation caused by 1992 Landers earthquake, where the co-seismic horizontal displacement was 3 m measured geodetically. They confirmed that InSAR could be used as a tool to obtain location, magnitude and type of an earthquake which otherwise could be only possible with seismological observation. Sreejith et al. [22] utilized InSAR technique for the coseismic and post-seismic investigations of the Gorkha earthquake. Satyabala and Bilham [21] successfully used this approach to study co-seismic deformation associated with the 1995 Chamoli earthquake. Chandrasekhar et al. [5], Saraf et al. [20], Rastogi et al. [17], Zia et al. [26] applied this technique to map post-seismic deformation correlated with 2001 Bhuj earthquake. Bhattacharya et al. [3] approximated interseismic deformation along parts of Himalayan Frontal Thrust using InSAR. [24] used radar observations and produced damage proxy maps derived from temporal changes in interferometric SAR coherence. Further Chatterjee et al. [6] also utilized this procedure to map surface deformation linked with fluid extraction. The accuracy of the InSAR measurements is affected by phase decorrelation, orbital errors, topographic residuals, phase unwrapping errors, imperfect simulation of the imaging geometry, and extra path delay due to the propagation of the microwave signal through the atmosphere [9]. In 2014 Fattahi & Amelung derived formulas for the uncertainty of InSAR velocity fields as a function of the orbital uncertainties. In 2015 Fattahi & Amelung concluded that the systematic component of the tropospheric delay arises because of seasonal variations in the moisture content of the atmosphere which biases the InSAR displacements and velocities. The magnitude of the bias depends on the amplitude of the seasonal variations and the sampling of the seasonal cycle by SAR data acquisitions.
Literature suggested that magnitude and intensity of the earthquake help to assess the vulnerability to earthquake. The information pertaining to earthquake epicenter, magnitude, depth to hypocenter, elevation, demography, etc., is analyzed in geospatial environment to deduce its relation with geotectonic settings and possible risk in the area [12]. Zonation is important as it provides information regarding decision making for urban regional planning. Kumar and Pandey, [12] analyzed that the earthquake events in Nepal were discovered at location deeper than 40 Kms during 1900-2014 as correlated to the event of 2015, where majority of earthquake events were recorded at below 10 km depth.
An artificial neural network (ANN) is a massively parallel-distributed-information-processing system that has certain performance characteristics resembling biological neural network of human brain [11]. It makes use of different layers of mathematical processing and activation functions to make generalization of the information it's fed [2]. Though being an important area of research in all fields in twenty-first century, its utilization to earthquake damage detection is comparably new area of study. Narayanakumar and Raja [15] worked on the earthquake magnitude prediction in Himalayas using back-propagation artificial neural network model and concluded that this model provides higher prediction efficiency for earthquake of magnitude ranging between 3 and 5, magnitude ranging between 3 and 5. Cooner et al. [7] used the feed forward ANN to detect the urban damage due to 2010 Haiti Earthquake and found that its error rate was very low, i.e. < 40%.
The most recent and devastating earthquake that occurred on 25th April 2015, also known as Gorkha Earthquake, left remarkable disfigurement to the people of Nepal as well as a challenge to the modern science and technology too since it came without semaphoring or detecting where and when it striked. Earthquake occurred at 11:56 a.m. local time with a Moment magnitude (Mw) of 7.8 on Richter scale having epicenter near Barpak village of Gorkha district which is about 181 km northwest of Kathmandu, with a focal depth of 15 km [1]. Out of 75 districts, 31 were deemed 'most affected' , among which 14 districts namely Gorkha, Dhading, Rasuwa, Nuwakot, Kathmandu, Lalitpur, Bhaktapur, Kavrepalanchowk, Sindhupalchowk, Dolakha, Sindhuli, Makawanpur, Ramechhap, Okhaldunga were 'severely affected' . It also affected some parts of India, Bangladesh and Tibet autonomous region of China. Tremors were also felt in Bhutan and Pakistan. This event is regarded as the largest event to have occurred in the zone of Main Himalayan Thrust (MHT) since the great 1934 Bihar-Nepal earthquake of the magnitude (Mw) 8.2. This earthquake, called the 2015 Gorkha earthquake, was supervened by large aftershocks to the west and east of the main shock rupture [14]. Over thirty-eight aftershocks of magnitude 4.5 or greater occurred in the day following the initial earthquake, including the one of magnitude 6.8 [23]. The earthquake and its aftershocks were the result of thrust faulting in the Indus-Yarlung Suture Zone, a thin east-west region spanning roughly the length of the Himalayan ranges. According to NPC 2015 (https:// www. world bank. org/ conte nt/ dam/ World bank/ docum ent/ SAR/ nepal/ PDNA% 20Vol ume% 20A% 20Fin al. pdf ) report about half of a million houses were destroyed and more than 250,000 houses were partially damaged leaving 8 million people homeless while the economic loss was estimated at $7 billion [18]. Compared with the fatality rates of the 1934 earthquake, the fatalities of the 2015 earthquake are significantly smaller for a given estimated intensity, except in some hill districts [19].
The present study is carried out to assess the displacement and damage caused due to the Gorkha Earthquake, 2015, by performing SAR interferometry technique over SAR SLC images of pre-and post-earthquake and ANN approach, respectively. This study will help in assessment of land deformation and the infrastructures displacement. It can also be helpful in determining the severity of damage as well as to analyze the extent of future damage if such events take place.
The objectives of the present study were: (i.) to prepare a regional earthquake prone zone map using epicenter data, (ii.) to assess the deformation and relative change in elevation by generating interferogram, and (iii.) to prepare a damage proxy map using ANN approach.

Study area
The present study was carried out in eastern and central regions of Nepal which are spatially located between 86°40′ E and 88°0′ E longitude and 26.22° N-28°0′ N latitude and 84°27′ E to 86°10′ E longitude and 27°0′ N-29°0′ N latitude, respectively, covering area of 28,456 sq.km and 27,410 sq.km individually.
Nepal occupies nearly one-third of 2400-km-long Himalayas with a small part of southern most Nepal stretching into the Indo-Gangetic Plain and two districts in the northwest stretching up to the Tibetan plateau. Since it lies completely within the collision zone, this makes it more prone to earthquake.
As per the census 2011, the population of Eastern Nepal was 5,811,555 (which is 21.9% of the total population of Nepal) with population density of 204 person per sq.km. In contrast the Central Nepal population was estimated as 9,656,985 (36.44% of total population of the country) with population density of 352 person per sq.km. About 18.6% of the total population of Nepal was living in urban areas and the rest 81.4% comprised of rural population.
During summer season, the temperature in both regions varies between 28 °C and 40 °C, whereas average rainfall was recorded to be 200-375 mm. Kosi, a tributary of Ganga river, and Mechi, a tributary of Mahananda river, are the two most important rivers of Eastern Nepal, whereas Bagmati river, which is a tributary of Kosi river and drain a part of Ganga basin, comprises the most important river in Central Nepal.

Materials
In this study for the preparation of regional earthquake prone zone map, earthquake epicenter data obtained from NASA website (www. earthquake.usgs.gov) has been used for the period from 2000 to 2019 and the map was prepared in ARCGIS 10.3 environment. For the interferometric process Sentinel-1 SLC images were used. For the Central Nepal (Study Area A), the VV polarization images of 17th April 2015 and 23rd May 2015 were used, respectively, as pre-and post-earthquake data. Similarly, for monitoring the land subsidence in Eastern Nepal (Study Area B) VV polarization SLC images of 24th April 2015 and 6th May 2015 were used, respectively, as pre-and post-earthquake data sets. Since the magnitude of both the earthquakes of 25th April and 12th May was almost same so taking datasets of different time interval and dates will not cause any error in the final result. The entire interferometric process was carried out in SNAP software. Snaphu export has been executed in Snaphu-v1.4.2_win64. For subtraction of the topography, SRTM 1 s HGT DEM was used. SRTM DEM images were also used to map slope, relief, aspect as well as to extract drainage, whereas LANDSAT-8 images were considered to map lineaments. Optical satellite images were processed using ERDAS IMAGINE (ver 2014). To prepare damage proxy map, ANN was employed in Rstudio (ver 3.5.1).

Sentinel 1A SAR satellite data
Sentinel-1 is a space mission from ESA of the Copernicus Programme, consisting of a constellation of two satellites. The payload of Sentinel-1 is a Synthetic Aperture Radar in C-band that provides continuous imagery (day, night and all weather) with repeat cycle of 12 days. Single Look Complex (SLC) product with VV polarization has been used to acquire displacement. SLC product was used as it includes single look in each dimension using the full transmit signal bandwidth and consists of complex samples preserving the phase information. The SAR data was acquired from Alaska Satellite Facility (ASF).

SRTM 1 arc second global DEM data
SRTM 1 Arc Second Global elevation data offer worldwide coverage of void filled data at a resolution of 1 arc-second (30 m) and provide open distribution of this high-resolution global data set. SRTM 30 m DEM has been used to obtain slope, aspect, relief and to delineate the drainage network using ARC spatial analyst tools.

Landsat 8 optical satellite imagery
The Landsat 8 series optical satellite data were obtained from USGS-Earth Explorer (earthexplorer.usgs.gov). Data was used to delineate lineament of both the study areas. The two main sensors for Landsat 8 are the Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS). The Operational Land Imager (OLI) produces 9 spectral bands (Band 1-9) at 15, 30 and 60 m resolution. Then, the Thermal Infrared Sensor (TIRS) consists of 2 thermal bands with a spatial resolution of 100 m.

Methods
In this study to map regional earthquake-prone zone, the USGS earthquake point data were interpolated using IDW interpolation in ARC map environment. To acquire the displacement, the Sentinel 1 SLC pre-and post-images were processed in the SNAP software. Further for assessing the damage caused ANN approach has been applied using Rstudio. The assessment was made based on 6 parameters which included slope, relief, aspect, drainage network, lineament and the displacement obtained after processing SAR data of each study area.

Processing of SAR data
Since, SAR interferometry requires pixel-to-pixel match between common features in the SAR image pairs, so coregistration of two SAR images is an essential step for the accurate determination of phase difference. It assures that each ground target subsidizes to the same pixel in both master and slave image. For co-registering the images of Central Nepal, the image acquired on 17th April 2015 was used as master image, whereas 23rd May 2015 has been considered as slave image. Similarly for coregistering Eastern Nepal SLC images the images obtained on 24th April 2015 and 6th May 2015 were considered as master and slave images, respectively.
After co-registration of the SLC images, multi-look complex interferogram is generated by multiplying the master image with its adjoint image. The interferogram shows the fringes which is related with the information of the topography on the Earth.
Since SLC products are composed of several burst overlapping in azimuth time for each sub swath, separated by black lines debursting is important. The deburst operation consists of generating a continuous image by removing black separation lines. The generated interferogram has been flattened by removing the topographic phase; this process is known as Differentia SAR interferometry (DInSAR) technique. SRTM 1 s HGT DEM was used for this purpose.
Debursting process is followed by coherencing and filtering of images. Coherence is defined as the correlation between two SLC SAR images. The correlation coefficient (q) can be used to quantify the coherence and is defined using the given equation: where E[.] is the mathematical expectation operator and is approximated by the mean of a finite number of elements, Z1 and Z2 are the complex SAR images, |q| is the absolute coherence measuring interferometric phase quality, and u is the phase of the complex correlation coefficient q.
Further the SAR image quality has been improved using multi-looking process. It is done to produce a product with nominal pixel image and also improves the radiometric resolution.
After getting the smoothened interferogram, phase unwrapping of the image has been carried out. Phase unwrapping is the consequence of the fringe counting problem in fringe pattern processing. When a large discontinuity occurs in the reconstruction, a 2π or multiples of 2π is added to the adjoining data to remove the discontinuity.
So, the key to phase unwrapping is the reliable detection of the 2π phase jumps. It can be done using SNAPHU tool which was used in windows command prompt by exporting the filtered image and then imported after processing. The absolute calibrated and unwrapped phase values are converted to a scalar unit representing displacement in meters. Range-Doppler terrain correction is used to correct such geometric distortion which is used to get an orthorectified end product.

Preparation of damage proxy map
Damage proxy map of both the regions was prepared using artificial neural network (ANN) approach in Rstudio. ANN was performed by importing the required libraries. The factors used are Slope, Drainage, Relief, Aspect, Displacement, Lineament buffer. Instable slope is considered as one of the major causes of earthquake. Many studies suggested that fracturing of the rock mass within active slope instabilities can amplify shaking up to 8 (Moore 2012). As the angle of slope (i.e. steeper slope) increases, it increases the chances of severe earthquake the region. Aspect values indicate the directions the physical slopes face. The aspect directions can be classified based on slope angle with a descriptive direction. Relief represents the elevation of a region. High elevation area is more prone to earthquakes in comparison with the low elevated area (Chuang 2017). The drainage system is among the controlling factor of earthquakes. The appropriate drainage system promotes infiltration of water into slope; higher number of drainage system running on or along the base of the system will bring about more desirable system for water infiltration. The earthquakes generally occur more near the first order stream. Surface displacement due to earthquake tells about the severity of damage caused (Tables 1, 2 and 3). 1 and 0 have been assigned using array format in r-Studio statistical software for each classes of all the factors. The result weight of all classes along with their generalized weights was obtained followed by error and steps of backpropagation neural network (BPNN) (Figs. 1, 2, 3).

Earthquake-prone zone map
Earthquake-prone zone map represents the potential for earthquake occurrence in a particular area. The map shows the earthquakes of varying magnitude in different areas. In this study earthquake-prone zone maps of Nepal region using magnitude and depth data for the period from 2000 to 2019 were prepared (Fig. 4a, b).
In the map generated using epicenter depth data, the area has been divided into 5 zones on the basis of depth (< 5, 5-15, 15-25, 25-35, > 35). The result exhibits that the majority of earthquake events took place at 5-15 km  (Table 4).
For the earthquake magnitude zone map, the area has been divided into 5 regions on the basis of their magnitudes. The 5 zones are < 2, 2-3, 3-4, 4-5, > 5. In the map it can be seen that the Central Nepal and far Western Nepal are prone to earthquake, in comparison with other regions. It was found that the earthquake of magnitude between 4 and 5(422 numbers) struck the region the most and the earthquake of magnitude greater than 5 hit the region the least number of times, i.e., 54. However the number of earthquake of magnitude less than 4 was only 88 during this period of time (Table 3).
To detect the relationship between the magnitude and depth of the earthquake, based on the data for the period 2000-2019, it was found that no significant relationship exists between the two factors (Fig. 5).

Generation of interferogram and assessing ground deformation
An interferogram is fundamentally a pattern formed by wave interference; it represents the variation in the phase of the radar waves before and after earthquake. In the figures shown, it can be seen that the fringes are prominent in the northern part of Central Nepal and southern part of Eastern Nepal (Fig. 6a, b). Here the fringes show the different phase cycle. One complete cycle of phase corresponds to one revolution of a color wheel. The repetition of these colored represents wrapping of the phase values. Each cycle represents a rise in elevation of 1/2ƛ. In figure region with wider space between the fringes indicates gentle change. Those with closer fringes show steep changes. And areas with no fringes are the places where no deformation took place.
Moving from centre to outward direction if the color of the cycle of the fringe is in order of cyan, yellow and magenta, then it is subsidence of the surface, and if the order is magenta, yellow, cyan, then it demonstrates the region where there is upliftment of the land surface. Ground deformation is defined as the Earth's surface movement (gradual settling, sudden sinking, or uplift) due to the subsurface movement of Earth materials, which has become a common geological hazard that has caused potential threats to public safety, bringing about a series of harmful impacts.
In the maps shown below, the terrain corrected displacement images have been interpreted. The slant range in the Central Nepal is between − 0.61 and 1.1 m. The south eastern region shows subsidence while the north eastern region shows the upliftment (Fig. 7), whereas in case of Eastern Nepal the slant range varies between − 0.33 and 1 m. The Eastern region shows downliftment and the southern region indicates the upliftment (Fig. 7). This variation in slant range might be due to difference in the slope, rainfall rate or presence of different rocks in the regions. The displacement value may include the effect of tropospheric water vapor

Damage proxy map
The damage proxy map of both the regions was prepared using slope, drainage, vertical displacement, relief, aspect, lineament buffer as the factors influencing damage during an earthquake (Fig. 8a, b). The bivalued function was assigned to the following classes of various factors, i.e. 1 for those who have experienced earthquake and 0 for those who do not experience any earthquake, using point shapefile obtained from www. earthquake.usgs.gov. The weights of various factors have been processed into first hidden layer along with the inputs, and then it has been processed into second hidden layer and then gets back propagated if error values are detected. It propagates 53 times and provides error of 0.018943 in case of Central Nepal, while for Eastern Nepal propagation time and error provided is 46 and 0.01295, respectively (Fig. 9a, b).
The ratio of the training set to test set was 75:25 percentage. And in the validation of the damage map, the artificial neural network showed prediction accuracy of 88.01%. The overall weight obtained is in optimized and in normalized form (0-1) influenced by biases and removal of errors in many steps. The weights to different classes of different factors have been assigned on the basis of the intermediate weights shown in the ANN architecture. And the overall weight of each factor is obtained by calculating the mean of weights of different classes of each parameter (Table 5). Table 6 represents the percentage influence of various factors obtained on the basis of overall weight of an individual factor. The factor with more overall weight has been assigned high percentage (%) influence and vice versa for each study area. The sum of the percentage influence of individual study area should always be equal to 100.
The final damage map was plotted by performing weighted overlay analysis of different parameters in Arc Environment. Taking the normalized weight obtained using ANN and percentage influence of all the factors in account.
The damage proxy map of both the regions was divided into five zones, viz. low zone, moderate zone, high zone, very high zone (Fig. 10) using natural break method. The value between 1 and 2 was categorized as low damage zone, 3-5 as moderate zone and 6-7 and 8-9 as high and very high zone, respectively, for both central and eastern Nepal region. In the Central Nepal damage proxy map moderate zone covers the maximum area (505.66 sq.km). The central region shows high to very high damage zones covering 364.28 sq. km and 93.14 sq.km (Table 7) area individually. It was found that the area covering compact build-up sighted very high damage. These area included Kathmandu, Kirtipur, Lalitpur, Deurali in southern central part and in northern central part it comprises of Laqu.Qu, Paiku Co., Chanong-Maierle, Se Marancuo. Along with the compact coverage of build-up, the area experiencing high damage shows high slope, high relief and high vertical displacement (upliftment) of surface, while the area revealing low vertical displacement (subsidence) experiences moderate damage. These areas also exhibit high number of first order drainage stream. In Eastern Nepal damage proxy map, the north western region (comprising of Bhojpur, Wana, Sisuwakhola, Chaukidada, Saphapokhari) and north eastern region (Phawakhola, Dummrise, Tiringe,  Mamangkhe, Dokhu) witnessed high (317.29 sq. km) to very high (270.33 sq. km), which has high slope, high relief (Table 7). Area experiencing high as well as low vertical displacement of surface shows very high damage. The area revealing low to moderate damage zone is the southern region which has dense build-up but is topographically plane (low slope, low relief ).

Conclusion
The present study analyzed the past century earthquake events of Nepal at various levels to understand the patterns in spatio-temporal dimension. The result indicated that majority of earthquakes event was of magnitude > 4 and the maximum depth of the earthquake was found to be between 5 and 15 km. A comparative study, first of its kind, demonstrates the methodology to estimate the surface deformation and damage zonation due to Gorkha Earthquake, 2015, in Central and Eastern Nepal using interferometric Synthetic Radar (InSAR) on Sentinel 1 SLC data sets. The Eastern Nepal showed an upliftment of 1 m and 0.33 m subsidence whereas in Central Nepal an upliftment of 1.1 m and subsidence of 0.61 m clearly reveal more surface displacement in area near to the epicenter of the earthquake event.  Damage proxy map prepared using artificial intelligence (ANN) back-propagation model shows that the damage caused due Gorkha earthquake 2015 was more intense in Eastern Nepal when compared to Central Nepal as the area under high damage was 270.33 sq.km in Eastern Nepal while in Central Nepal it was 93.14 sq.km.
A relation between the buildup cover and damage caused has also been divulged with compact buildup area suffering from more damage compared to the sparse buildup area. The study also affirmed that the damage depends on many factors including slope, drainage, relief, aspect, presence of lineament, surface displacement, magnitude and depth of the past earthquakes in the area.
Limitations related to SAR data's sensitivity to changes in acquisition geometry and surface variations were challenging. Further, in view of the contemporary seismicity and conspicuous displacements, study of long-term observations of this surface movement    is recommended to reveal the dynamic process in this region. The study could be carried out on other highresolution SAR sensors such as ALOS PALSAR, TER-RASAR-X and COSMO-SkyMed to obtain better accuracy in results.
Data Availability Data can be made available through user request.
Code availability Not applicable.

Conflict of interest
The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.