Three-Dimensional Nepal Earthquake Displacement Using Hybrid Genetic Algorithm Phase Unwrapping from Sentinel-1A Satellite

Abstract

Introduction: Geophysicists had forewarned for decades that Nepal was exposed to a deadly earthquake, exceptionally despite its geology, urbanization and architecture. Gorkha earthquake is the most horrible natural disaster to crash into Nepal since the 1934 Nepal-Bihar earthquake. Gorkha earthquake occurred on April 25, 2015, at 11:56 NST and killed more than 10,000 people and injured more than 23,000 population. Objective: The main objective of this work is to utilize hybrid genetic algorithm for three-dimensional phase unwrapping of Nepal earthquake displacement using Sentinel-1A satellite. The three-dimensional best-path avoiding singularity loops (3DBPASL) algorithm was implemented to perform 3D Sentinel-1A satellite phase unwrapping. The hybrid genetic algorithm (HGA) was used to achieve 3DBPASL phase matching. Advancely, the errors in phase decorrelation were reduced by optimization of 3DBPASL using HGA. Results: The findings indicate a few cm of ground deformation and vertical northern of Kathmandu. Approximately, an area of 12,000 km2 has been drifted also the northern of Kathmandu. Further, each fringe of colour represents about 2.5 cm of deformation. The large amount of fringes indicates a large deformation pattern with ground motions of 3 m. Conclusion: In conclusion, HGA can be used to produce accurate 3D quake deformation using Sentinel-1A satellite.

1. Introduction

Satellite-based interferometric synthetic aperture radar (InSAR) is a potential tool for precise measurements of ground shifts triggered off by earthquakes or landslides [1–4]. Whilst a satellite sweeps over a territory, it can pick up the amplitude and phase of radar pulses backscatter off Earth’s surface [5–7]. Subtracting the phases logged in different passes over the same region acquiesces an interference pattern which is sensitive to changes in topography along the satellite’s line of sight within few millimetres [6, 8]. In other words, interferometric synthetic aperture radar (InSAR) techniques have excellent potentials to measure the millimetre scale of the Earth’s surface deformation. Nevertheless, the foremost difficult raised up in InSAR techniques is phase unwrapping [8, 9].

The keystone questions are how the earthquake can generate based on plate tectonic theory and how InSAR can be used to detect earthquake deformation. The following two sections are devoted to understand the mechanisms of earthquake and InSAR technology. In this paper, we address the question of utilization of inversion genetic algorithm (GA) as an optimization methodused for three-diensional modelling of rate deformation due to Nepal earthquake displacement. In this context, Lee [10] used three-dimensional sorting reliabilities algorithm (3D-SRA) phase unwrapping for modelling volume rate changes of shoreline. However, 3D-SRA was not able to remove the artefacts in digital elevation model (DEM) due to radar shadow, layover, multipath effects and image misregistration and finally the signal-to-noise ratio (SNR) [11–16]. In fact, 3D-SRA does not identify singularity loops at all. It be determined by completely a quality measure to unwrap the phase volume. Ignoring singularity loops may cause the unwrapping path to penetrate these loops, and errors may propagate in the unwrapped phase map [17].

The two-dimensional unwrapping approaches, however, could introduce discontinuous regions when the noise is high. The resulting inconsistent baselines within a slice would produce an incorrectly unwrapped baseline. Then, the one-dimensional baseline unwrapping could give incorrect results. Many of the methods apply to quality map to guide the unwrapping procedures. The quality map was defined with the quality of the edges that connect two neighbouring voxels and unwrap the most reliable voxels first [18, 19]. Therefore, three-dimensional phase unwrapping approach considers the temporal domain and the spatial domain restrictions simultaneously [17, 18]. The main contribution of this study is to combine hybrid genetic algorithm (HGA) with three-dimensional phase unwrapping algorithm of three-dimensional best-path avoiding singularity loops (3DBPASL) algorithm with InSAR technique. Two hypotheses examined are (i) the HGA algorithm that can be used as filtering technique to reduce noise in the 3D phase unwrapping and (ii) 3D Nepal earthquake displacement that can be reconstructed using satisfactory phase unwrapping of 3DBPASL by involving HGA optimization algorithm.

2. Mechanism of earthquake based on plate tectonic theory

The main question that can arise is how do earthquakes generate? The answer of the above question required a comprehensive understanding of the theory of plate tectonic. Truthfully, plate tectonic theory is a starting point to capture the mechanisms of tsunami generation. Consequently, the plate tectonic theory let us have inclusive thought about generation of earthquakes and volcano which are the foremost causes for tsunami. Consequently, there are three versions for the plate tectonic theory. The first version has discussed the lithosphere broken into strong, rigid moving plates that carry the continents on their backs ocean basins that come and go at mid-ocean ridges and subduction zones, respectively (Figure 1) [20].

Figure 1.

First version of plate tectonic theory.

Further, Marghany [21] has accepted that physical plate margins fluctuate significantly sideways belt (i.e. beginning location to location sideways their dimensions), combined plate margin kinds are usual. Therefore, plate margins are not completely unbending Successively and an entirely plate exchanges would correspondingly explicate, otherwise by slightest crack through, huge interior features similar mountain stripes and extensive faults intraplate distortions alike the rock-strewn mountains and washbasin. For instance, Hawaiian Islands represent the foremost mid-plate volcanic disruptions (Figure 2) [21]. This agrees with the study conducted by Marghany [20].

Figure 2.

Northern Pacific floor with Hawaii-Emperor and many other volcanic chains.

Marghany [22] stated that the major feature of plate tectonics is the lithosphere. With this regard, it is considered as unstable thermodynamic system which is staging to disintegrate the Earth’s crust and mantle thermal. Further, Nitecki et al. [23] reported that plates are dim in expansion and penetrable to liquefy from the original, hardly dense asthenosphere. The absolute dynamic influences ahead of plate movements are the vertical cooling fluctuations and gravity force. The lower mantle is heating the upper mantle through conductivity. Further, 60% of seafloors and subduction cool the lower mantle. Consequently, subduction which is correlated to joint rollback and overruling plate leeway which have disappeared, for instance, the essential profiles of the Earth’s surface. Extreme volcanism is considered as hot spots which concentrated by expansion of plate weakness. Besides, hot spots can be caused by limited and comparatively shallow thermal instabilities which are associated to thermal fluctuation of the upper mantle because of dynamic motions of adjacent plates [22].

Understanding the main features of plate tectonics can lead to extremely perceptive of earthquake occurrence mechanisms. According to Marghany [21], plate tectonics are associated by major features which involved:

• The surface of Earth is constructed of chains of huge plates.

• Per year, these plates are in steady movement drifting at a limited centimetres.

• The ocean bottoms are frequently shifting, scattering starting the midpoint and descending on the boundaries.

• The plates are moving in various paths due to the effects of convection current below the plates.

Consistent with above prospective, the radioactive decay is the main source of thermal fluctuations which is creating the convection current flows. This occurs deeply in the Earth which is deriving vertical density exchanges of fluid through the Earth. With this regard, less dense of hot rocks below the mantle rises above cooler rock which emerged down due to the impact of gravity on their high densities (Figure 3) [21].

Figure 3.

The convection currents of the Earth.

The comprehensive understanding of tectonic plate interactions is keystone of earthquakes, volcano and tsunami generations. With this regard, the principles of tectonic plate interactions are (i) divergent boundaries, (ii) convergent boundaries and (iii) transform boundaries (Figure 4) [20].

Figure 4

Three types of plate boundary.

Consequently, divergent boundaries represent zones where plates deviate from each other which are forming, for instance, (Figure 5) either mid-oceanic ridges or rift valleys. With these regards, divergent boundaries are also identified as constructive boundaries [20]. The Atlantic Ocean was created by this process. The Mid-Atlantic Ridge is an area where a new seafloor is being created [24–26].

Figure 5.

Divergent boundaries.

Convergent boundaries are well known as compressional or destructive margins. Along these boundaries, the plates approach each other and crash (Figure 6). These boundaries involved (i) subduction zones, (ii) obduction and (iii) orogenic belts. Subduction zones arise where an oceanic plate convenes a continental plate and is broke underneath it. Subduction zones are indicated by oceanic trenches (Figure 4). The descendent completion of the oceanic plate melts down and constructs compression in the mantle, triggering volcanoes to create.

Figure 6.

Convergent boundaries.

According to Marghany [24], the continental plate can make obduction arises up as soon as it is plugged underneath the oceanic crust. Nonetheless, this phenomena is considered unfamiliar, for instance, the plate tectonics have absolute densities which enable the oceanic crust’s subduction. This heritage, the oceanic crust to fasten and frequently upshot in a different fangled mid-ocean ridge, for instance, is evolving and spinning the seizure hooked on subduction [25–28].

Consequently, the collision boundaries are known as orogenic belts [24, 28]. These collision boundaries are produced because of collision of two plates which thrust upwards to create bulky mountain ranges, for instance, the Himalayas, the highest mountain range on Earth (Figure 7).

Figure 7.

Mountain ranges of the Himalayas.

Furthermore, Marghany [21, 24] and Ferraiuolo et al. [27] agreed that the most spectacular geological feature on the surface of the Earth is orogen restraint. This can be seen between crust boundaries of Indo-Australian crust and African crust and to the north of the Eurasian crust. This restraint runs from beginning of New Zealand, cutting edge of the East-South-East, from side to side of Indonesia, sideways of the Himalayas, from end to end of the Middle East and moves up to the Mediterranean in the West-Northwest. It is also termed the ‘Tethyan’ District, as it creates the zone lengthways which the earliest Tethys Ocean was malformed and vanished [25].

Figure 8 shows occurrence of transform boundaries when two plates suppress past each other with only partial convergent or divergent movement.

Figure 8.

Transform boundaries.

The San Andreas Fault in California, for instance, one of the well-known transform boundaries, is an active transform boundary. The Pacific Plate (carrying the city of Los Angeles) is moving northwards with respect to the North American Plate [26].

The critical question is how to prove the occurrence of plate tectonics? The landmasses appear to be appropriate and organized, similar a massive picture puzzle. In this respect, when anyone goes into a map, Africa appears to nestle agreeably addicted to the Caribbean Sea and the east coast of South America. In line with Ferraiuolo et al., [27], in 1912, Alfred Wegener suggested that these two landmasses were as soon as combined with each other formerly someway moved away from each other. He identified this procedure as Pangea. He assumed that Pangea was together up to 200 million years ago. The map of plate tectonics has been established by NASA and is shown in Figure 9 which has an excellent evidence of the existence of plate tectonics [27].

Figure 9.

Current map of plate tectonics.

3. Continental drift

In 1915, the theory of continental drift was proposed by the German geologist and meteorologist Alfred Wegener, which declares that parts of the Earth’s crust slowly drift atop a liquid core. The fossil record verifies and provides credibility to the theories of continental drift and plate tectonics. The theory of continental drift has been counted by the theory of plate tectonics, which clarifies how the continents shift (Figure 10). The impression that continents can drift about is termed, not surprisingly, continental drift. The old (and very strong!!) theory earlier was the ‘contraction theory’ which recommended that the Earth was once a molten ball and in the process of cooling the surface cracked and folded up on itself [28].

Figure 10.

Continental drift.

The critical issue with this clue was that all mountain ranges should be roughly of the same age, and this cannot be proper. It was suggested that as the continents relocated, the primary edge of the continent would confront resistance and thus constrict and bend upwards creating mountains close to the foremost edges of the drifting continents. Consistent with Ferretti et al., [28], there is contrariety that the key protest to continent drift was that there is no mechanism, and plate tectonics was recognized without a mechanism, to transfer the continents.

However, plate tectonics are the extensively conventional theory that Earth’s crust is broken into rigid, transferring plates. In the 1950s and 1960s, scientists revealed that the plate edges through magnetic surveys of the ocean floor and through the seismic heeding networks assembled to examine the nuclear trying. Discontinuous patterns of magnetic anomalies on the ocean floor specified seafloor spreading. With this regard, a new plate material is born. Magnetic minerals affiliated in ancient rocks on continents correspondingly exposed that the continents have lifted completely to one another [25–28]. For instance, India drifted northwards into the Asia forming the Himalayas and of course Mount Everest (Figure 11).

Figure 11.

Himalayas theory of drifting tectonic.

The East Asian Sea Plate was an unknown tectonic plate which has been swallowed up by the Earth. It can be found in the Philippine Sea. In fact, the Philippine Sea locates at the juncture of numerous foremost plate tectonics: (i) the Pacific, (ii) Indo-Australian and (iii) Eurasian plates. These set up numerous minor plates, containing the Philippine Sea Plate. Since approximately 55 million years ago, the Philippine Sea Plate has been drifting northwest [27].

4. Nepal earthquake

Gorkha earthquake is the most horrible natural disaster to crash into Nepal since the 1934 Nepal-Bihar earthquake. Gorkha earthquake occurred on April 25, 2015, at 11:56 NST, killed more than 10,000 people and injured more than 23,000 population. Its epicentre was the east of the District of Lamjung, and its hypocentre was with an approximately depth of 15 km with maximum magnitude of 8.1 M_w. Consequently, within 15–20 min, aftershock was struck across Nepal with a magnitude of 6.7 on the 26th of April at 12:54:08 NST. Thus, the epicentre of a foremost aftershock was close to the Chinese border between the capital of Kathmandu and Mount Everest (Figure 12) with a moment magnitude of 7.3 M_w [29].

Figure 12.

Location of Nepal’s earthquake.

Consistent with the USGS [30], the temblor was triggered through an abrupt push. Further, the accumulative energy underneath the surface of the Earth produced huge stress on the main fault which is associated with the Indian crust. This was contributed to drift the Indian crust further down the Eurasian crust. With this regard, Kathmandu sited on a slab of layer almost 120 km widespread and 60 km elongated. Within 30 s, the crust was drifted 300 cm southward.

Counter to Hussein et al. [31], Nepal situated southwards which is closed to the slab crust somewhere the Indian crust moved down of the Eurasian crust. This was the plate dominating the essential locality (Figure 13) within 800 km of the Himalayan arch. Organically, the Nepal Himalayas are segregated into five tectonic district morpho-geotectonic zones: (1) Terai valley, (2) Sub-Himalaya, (3) Lesser Himalaya (Mahabharat Range and mid-valleys), (4) Higher Himalaya and (5) Inner Himalaya (Tibetan-Tethys) [32, 59, 60, 61]. In central Nepal, the convergence rate between the plates is about 45 mm/year. According to USGS [30], the location, magnitude and focal mechanism of the earthquake indicate that it was triggered off by a slip along the Main Frontal Thrust (Figure 14).

Figure 13.

Topography map of Nepal.

Figure 14.

Nepal Himalayas are partitioned into five tectonic zones.

The earthquake’s effects were amplified in Kathmandu as it sits on the Kathmandu Basin, which contains up to 600 m of sedimentary rocks, representing the infilling of a lake. Wei and Cumming [33] stated that the latest quake follows the same pattern as a duo of big tremors that occurred over 700 years ago and results from a domino effect of strain transferring along the fault. The last time the fault ruptured at this location was back in 1344. It was preceded in 1255 by a big event to the east of Kathmandu. The last rupture there was in 1934, hinting strain might accumulate westwards. This means that 2015’s quake follows the pattern with a gap between events of 80 years or so. In other words, the Main Frontal Thrust, on average of a great earthquake, occurs every 750 ± 140 and 870 ± 350 years in the east Nepal region. Perhaps, the 700-years delay between earthquakes in the region. Finally, it can be suggested that the 2015 quake was due to the tectonic stress buildup from the 1934 quake [32].

5. Interferometric synthetic aperture radar for earthquake monitoring

The SAR interferometry is established on the SAR machinery. The abbreviation SAR tolerates for synthetic aperture radar. SAR mechanism is a function of signal transmission-reception of a percentage of energy that interrelates with the surface, which is denoted as backscattered, existing a determining strength and time delay of the received signals [34].

Consistent with Massonnet et al. [1], InSAR is a skeleton covering various approaches or procedures that expend phase knowledge obtained by monitoring phase variance or status of displacement of the microwave at the instantaneous which is obtained by the radar antennas between two SAR imageries. With this regard, it is recognized as master and slave images. Therefore, the interferometric phase is assimilated from diverse sensor locations [1, 3, 5, 35, 36]. In this concern, map variations of digital elevation model (DEM) can be created by employing the two single-look complex (SLC) synthetic aperture radar (SAR) data which are obtained by two or extra disconnect antennas. As stated by Qifeng et al. [37], the output phase scene’s result is generated by multiplying the coregistered complex conjugate pixels and complex radar data. By the way, the difference of the two SAR phase data is handled to obtain DEM and/or the earth’s surface distortion. Consistent with Qifeng et al. [37] and Nan et al. [36], the absolute value of the real phase difference between any two neighbouring pixels is less than π, and then the real phase can be acquired by integrating the wrapped values of the wrapped phase differences. Nevertheless, wrapped phase image is extremely challenged by several factors. Particularly, these critical factors are involved multiplicative speckle noises, shadow, foreshortening, layover, temporal, geometric and atmospheric decorrelations [1, 3, 37–39] which are negatively produced in the area of nonstandard phase, i.e. low-quality area [36, 40, 41]. In this respect, low-quality area can contribute to critical decorrelation issues in the phase unwrapping procedures.

Several phase unwrapping algorithms have been introduced to solve the critical issue of low-quality area and the decorrelation. These algorithms are categorized into (i) path-following algorithms and (ii) minimum-norm algorithms [42]. Subsequently, minimum-norm algorithms express the unwrapping issue in terms of minimization of the global function as compared to path-following algorithms. Conversely, the constraint of minimum-norm algorithms cannot be involved every individual pixel in the phase unwrapping procedures [63]. In contrast, path-following algorithms are extremely advanced compared to minimum-norm algorithms. The advances of path-following algorithms are (i) to identify the residues and use the quality map to guide the generation of branch cuts by implementing branch-cut [41] and mask-cut algorithm [43], (ii) to guide the path of integration which is function of quality-guided algorithm and (iii) to dismiss the total discontinuities of the unwrapped result by using minimum discontinuity algorithm [44, 45]. Besides, there are other phase unwrapping algorithms. Therefore, the conventional phases unwrapping region-growing and the least-square algorithms require accurate image coregistration, which should be about 1/10 to 1100 resolution cell size (i.e. SAR pixel) [46]. In this context, Hai and Renbiao [5] stated that the interferometric phase unwrapping method which is based on subspace projection can provide an accurate estimation of the terrain interferometric phase (interferogram) even if the coregistration error reaches one pixel. This can be achieved by phase unwrapping algorithms such as the branch-cut method [11, 46].

Commonly, the accuracy of the image coregistration is a serious issue for accurate interferometric phase unwrapping. It is well known that the performance of interferometric phase unwrapping suffers seriously from poor image coregistration. Therefore, TanDEM-X and TerraSAR-X have been implemented to maintain the issue of image coregistration. Further, ERS-1 and ERS-2 and TerraSAR-X in tandem mode are the excellent examples of short temporal resolution. In wide range of contexts, TanDEM-X and TerraSAR-X are imaging the terrain below them simultaneously, from different angles. These images are processed into precise elevation maps with a 12 m resolution and any vertical accuracy better than 2 m [1, 3, 15, 39, 46, 47].

Further, Gens [14] described that the period change of two SAR acquirements is an additional kind of decorrelation. Certainly, the period modifications despite the fact associate with SAR data groups through a comparable dimension of baseline which attained one and 35 days which causes temporal component decorrelation. Therefore, Dickinson [48] stated that the loss of coherence in the same repeat cycle in data acquisition is most likely because of baseline decorrelation and dense vegetation covers in such a tropical as Malaysia. According to Nizalapur et al. [7], uncertainties could arise in DEM [14] because of the limitation of the InSAR repeat passes [3, 9, 15, 23, 36, 49, 56, 58, 62]. A persistent phase alteration throughout the two SAR data which is produced by the straight regular atmosphere was above the dimension measure of an interferogram and vertically above the landscape. The atmosphere, nevertheless, is alongside dissimilar on dimension measures equally greater and lesser than usual distortion signals [8, 23].

Recently, Zhong [44] stated that Differential Interferometric Synthetic Aperture Radar (DInSAR) has recently been applied with success to investigate the temporal evolution of the detected deformation phenomena through the generation of displacement time series. With this regard, two foremost kinds of progressive DInSAR practices of distortion period-chain origination have been recommended in several scientific works, which are frequently mentioned as persistent scatterer (PS) and small baseline (SB) procedures, correspondingly. The PS algorithms choose completely the interferometric SAR data pairs through indication to a single SAR master scene, deprived of slightly restraint on the sequential and three-dimensional parting (baseline) amongst the orbits.

6. Phase unwrapping

Phase unwrapping is a key step in modelling surface deformation from interferometric synthetic aperture radar (InSAR) data. Thus, the measured phases are wrapped in the interval (2π, π], and phase unwrapping is a necessary step to recover the original full phase value [4]. In particular, phase unwrapping algorithms can be classified into two categories: (i) the path-following method and the (ii) minimum-norm approach. Goldstein’s method is a classic path-following phase unwrapping method which explores for excesses cutting edge of the phase scene and employs method of branch cuts amongst closed interferogram errors. In this manner, it eliminates the mixing route dependency once the phase scene is perfectly unwrapped. The phase alteration mixing departs from the beginning of the upper leftwards corner of the scene. It followed by unwrapping the phases through the hypothesis of that contiguous pixels which have phase variances contained by the array of [−π, π]. Similarly, the unwrapping route requests to assimilate branch cuts with the intention of eliminating the route dependence. Nevertheless, after the branch cuts are positioned in the incorrect site, the wrong unwrapping is produced [4].

The minimum-norm approaches usually require extra calculation period than the route-succeeding mixing procedures. In Flynn’s algorithms, for instance, Gao and Yin [46] attempt to diminish the cutoff in phase results. From the algorithm it was discovered that the pathway can use the loop of incoherence and can eliminate the gap by toting or eliminating 2π on the phase standards inside the ring section. The procedure is achieved by pending to no extra incoherence loops which are created. In this respect, the incoherence placements can be directed by using a quality map [40].

In proportion to ENVISAT and Karout [19, 40], the keystone of the quality map concept is firmly based on the approach which unwraps laterally a route initial through the utmost excellence pixel. With this regard, slight unwrapping is compulsory to the furthermost vague spaces everywhere the phases of contiguous recovered superiority pixels are implemented to achieve precisely unwrapping levels. The quality map can be a correlation map, phase derivative variance map, etc. When a correlation map is not available, the phase derivative variance map is generally used. In the word of Mitri [50], the phase estimation is a major challenge to determine more accurate DEM. This is because of the measured phase differences which are given as a wrapped phase field of the principal values of a range −π to π; thus, the existence is unspecified within multiples of 2 π [7]. This procedure produces phase jumps between neighbouring pixels. Smooth function is used to resolve phase jump by adding or subtracting multiples of 2 π. Consequently, multichannel MAP height estimator based on a Gaussian Markov random field (GMRF) has developed by Wikipedia [51] and Frankel [52] to solve the uncertainties of DEM reconstruction from InSAR technique. They found that the multichannel MAP height estimator has managed the phase discontinuities and improved the DEM profile [53]. Recently, Bollinger et al. [54] validated GMRF technique with a lower range of error (0.01±0.11 m) with 90% confidence intervals for DEM reconstruction using RADARSAT-1 SAR F1 mode data.

However, the two-dimensional unwrapping methods could introduce discontinuous regions when the noise is high. The resulting inconsistent baselines within a slice would produce an incorrectly unwrapped baseline. Then, the one-dimensional baseline unwrapping could give incorrect results. Many of the methods are applied to quality map to guide the unwrapping procedures. The quality map was defined with the quality of the edges that connect two neighbouring voxels and unwrap the most reliable voxels first [31, 32]. Therefore, three-dimensional phase unwrapping approach considers the temporal domain and the spatial domain restrictions simultaneously [10].

7. Deformation reconstruction using three-dimensional best-path avoiding singularity loops (3DBPASL) algorithm

Three methods are involved to perform InSAR from Sentinel-1A satellite data: (i) conventional InSAR procedures, (ii) three-dimensional phase unwrapping algorithm, i.e. three-dimensional best-path avoiding singularity loops (3DBPASL) algorithm [19, 40, 55], and [iii] hybrid genetic algorithm. Consistent with Zebker et al. [4], SAR interferometry (equivalently, the InSAR technique) is a technique to extract surface’s physical properties by using the complex correlation coefficient of two SAR signals. The complex correlation coefficient, γ, of the two SAR observations, s₁ and s₂, is defined as

$\gamma =\frac{E\left[s_1{s}_2^{*}\right]}{\sqrt{E\left[s_1{s}_1^{*}\right]E\left[s_2{s}_2^{*}\right]}}$E1

where E[] is the mathematical expectation (ensemble averaging) and * represents the complex conjugate [2]. Further, Hanssen [3] stated that the interferometric phase is defined as the phase of the complex correlation coefficient as

$\varphi =arg\left\{\gamma \right\}=arg\left\{E\left[s_1{s}_2^{*}\right]\right\}\text{,}$E2

and its two-dimensional map is called the interferogram. Moreover, Hanssen [3] and Zebker et al. [4] defined the coherence as the amplitude of the complex correlation coefficient which expressed as

$\rho =\left|\gamma \right|\text{,}$E3

and its two-dimensional map is presented in the coherence image [1, 3, 5]. According to Qifeng et al. [37], an interferogram contains the interferometric phase fringes from SAR geometry, together with those from topography and displacement of the surface. The level of coherence can give a measure of the quality of the interferogram. Initially, the InSAR techniques were mainly dedicated to topographic information retrieval from interferograms. Further development resulted in techniques to extract interferometric phase fringes from coherent block displacement of the surface [26, 51–53]. In Lee [10], the surface displacement can estimate using the acquisition times of two SAR data s₁ and s₂. The component of surface displacement, thus, in the radar-look direction, contributes to the further interferometric phase (φ) as [8, 12]

$\phi =\frac{4\pi \left(\Delta R\right)}{\lambda }=\frac{4\pi \left(B_{h}sin\theta -B_{v}cos\theta \right)}{\lambda }$E4

where ∆R is the slant range difference from satellite to target, respectively, at different times, θ is the look angle (20–46°) [19, 55], λ is Sentinel-1A satellite wavelength single-look complex (SLC) of C-band. Therefore, B_h and B_v are horizontal and vertical baseline components [10]. Then, the phase difference Δφ between the two Sentinel-1A data positions and the pixel of target of terrain point is given by

$\Delta \zeta =\frac{\lambda Rsin\theta }{4\pi B}\Delta \phi$E5

Eq. (5) is a role of standard baseline B and the range R. Furthermore, Eq. (5) can deliver evidence approximately of DEM and phase variance approximations. In truth, the predictable DEM of SAR images is a significant mission to construct the 3D distortion. However, this equation required short baseline and accurate image coregistration to acquire accurate quality of interferogram phase map. Phase unwrapping in Eq. (5) can be extended to three-dimensional to

$\sum _{i,j,k}{w}_{i,j,k}^x{\left|\Delta {\phi }_{i,j,k}^x-\Delta {\psi }_{i,j,k}^x\right|}^P+\sum _{i,j,k}{w}_{i,j,k}^y{\left|\Delta {\phi }_{i,j,k}^y-\Delta {\psi }_{i,j,k}^y\right|}^P+\sum _{i,j,k}{w}_{i,j,k}^z{\left|\Delta {\phi }_{i,j,k}^z-\Delta {\psi }_{i,j,k}^z\right|}^P$E6

where Δφ and Δψ are the unwrapped and wrapped phase differences in x, y and z, respectively, and w represents user-defined weights. The summations are carried out in between x, y and z directions over all i, j and k, respectively. A more advanced method developed by Haupt and Haupt [42] is L^p-norm which uses similar methods like the two previous least squares methods to solve the phase unwrapping problem. However, this method does not compute the minimum L²-norm but the general minimum L^p-norm. In essence, by computing the minimum L^p-norm where p≠2, this method can generate data-dependent weight unlike the weighted least squares method. The data-dependent weights can eliminate iteratively the presence of the residues in the unwrapped solution. This method is more robust than the previous mentioned least squares method, and it is more computationally intensive [45]. Then the phase unwrapping based on the quality map can be calculated using the following equation [10, 40]

$\begin{array}{l}{Q}_{m,n,l}=\frac{1}{m\times n\times l}\ast {\left({(\sum (\Delta {\phi }_{i,j,k}^x-\overline{\Delta {\phi }_{i,j,k}^x})}^2\right)}^{0.5}\\ +\left({(\sum (\Delta {\phi }_{i,j,k}^y-\overline{\Delta {\phi }_{i,j,k}^y})}^2\right)+{\left({(\sum (\Delta {\phi }_{i,j,k}^z-\overline{\Delta {\phi }_{i,j,k}^z})}^2\right)}^{0.5},\end{array}$E7

where Δφ^x, Δφ^y, and Δφ^z are the unwrapped phase gradients in the x, y and z directions, respectively. ${\overline{\Delta \phi }}^x,{\overline{\Delta \phi }}^y,\mathrm{and}{\overline{\Delta \phi }}^z$ are the mean of unwrapped phase gradient in m × n × l cube in Δφ^x, Δφ^y, and Δφ^z, respectively. i, j and k are neighbours’ indices of the voxel v_in m, n and l cube. Following ENVISAT [19], the maximum gradient of the voxel v_m, n, l_ can be obtained by estimating the maximum unwrapped phase gradient of in the x, y or z directions, as described in Eq. (8):

$X_{m.n,l}=max\left\{max\left\{\left|\Delta {\phi }_{i,j,k}^x\right|\right\},max\left\{\left|\Delta {\phi }_{i,j,k}^y\right|\right\},max\left\{\left|\Delta {\phi }_{i,j,k}^z\right|\right\}\right\}$E8

According to ENVISAT [19] the maximum gradient method indicates the badness rather than the goodness of the unwrapped phase data, so the quality is calculated using the reciprocal of the unwrapped phase gradient of Eq. (8). Further, ENVISAT and Karout [19, 40] and Lee [10] agreed that quality-guided phase unwrapping algorithms are the function of the quality of the voxels themselves to conduct the unwrapping path and to minimize error propagation during the unwrapping procedure. In this respect, the unwrapping path algorithm is the function of the quality of the edges as an intermediate stage, rather the quality of the voxels [19, 40]. Following Lee [10] and Karout [40], the quality map of voxels can be given by

$Q_{m,n,l}=\sqrt{{O_x}^2\left(m,n,l\right)+{E_y}^2\left(m.n,l\right)+L^2\left(m.n,l\right)}$E9

where O_x, E_y and L are the horizontal, vertical and normal second differences, respectively [10], where

$\begin{array}{l}{O}_x\left(m,n,l\right)=\Delta {{\phi }^x}_{m,n,l}{\left[{{\phi }^x}_{m+1,n,l}-{{\phi }^x}_{m,n,l}\right]}^{-1}\left[\partial {{\phi }^x}_{m-1,n,l}-\partial {{\phi }^x}_{m,n,l}\right]-d{{\phi }^x}_{m,n,l}{\left[{{\phi }^x}_{m+1,n,l}-{{\phi }^x}_{m,n,l}\right]}^{-1}\\ \left[\partial {{\phi }^x}_{m,n,l}-\partial {{\phi }^x}_{m+1,n,l}\right]\text{,}\end{array}$E10

$\begin{array}{l}{E}_y\left(m,n,l\right)=\Delta {{\phi }^y}_{m,n,l}{\left[{{\phi }^y}_{m,n+1,l}-{{\phi }^y}_{m+1,n,l}\right]}^{-1}\left[\partial {{\phi }^y}_{m,n-1,l}-\partial {{\phi }^y}_{m,n,l}\right]-d{{\phi }^y}_{m,n,l}{\left[{{\phi }^y}_{m,n+1,l}-{{\phi }^y}_{m,n,l}\right]}^{-1}\\ \left[\partial {{\phi }^y}_{m,n,l}-\partial {{\phi }^y}_{m,n+1,l}\right]\text{,}\end{array}$$$E11

$\begin{array}{l}L\left(m,n,l\right)=\Delta {{\phi }^z}_{m,n,l}{\left[{{\phi }^z}_{m,n,l+1}-{{\phi }^z}_{m,n,l}\right]}^{-1}\left[\partial {{\phi }^z}_{m,n,l-1}-\partial {{\phi }^z}_{m,n,l}\right]-d{{\phi }^z}_{m,n,l}{\left[{{\phi }^z}_{m,n,l+1}-{{\phi }^z}_{m,n,l}\right]}^{-1}\\ \left[\partial {{\phi }^z}_{m,n,l}-\partial {{\phi }^z}_{m,n,l+1}\right]\text{,}\end{array}$E12

where m, n and l are the neighbours’ indices of the voxel in 3 × 3 × 3 cube and dϕ^z_m,n,l[ϕ^z_m,n,l + 1 − ϕ^z_m,n,l]^− 1 defines an unwrapping operator that unwraps all values of its argument in the range [−π, π]. This can be done by adding or subtracting an integer number of 2π rad to its argument [10, 19, 51–53]. Eqs. (10)–(12) are signified 3D displays of the gradient in unwrapped phase variations of dϕ^x, dϕ^y, and dϕ^z. Furthermore, the determined phase incline estimates the scale of the biggest phase slope which is fractional derived or draped of the phase alteration in Vⁱ, V^j, and V^k dimensions [40].

In general, the quality of separately border of phase unwrapping is a function of the connection of two voxels in 3D Cartesian axis, e.g. x, y and z. Starting to optimize the unwrapping path from high-quality voxels to bad quality voxels [19] by using hybrid genetic algorithm (HGA). Indeed, the 3D best path which is avoiding singularity loops (3DBPASL) algorithm may not be the shortest when the residues are dense as the isolated districts enclosed by the branch cuts can easily appear. Searching for the shortest path in the path-following algorithm is in fact an optimization problem. Several optimization approaches, such as genetic algorithms, have been applied to 2D phase unwrapping. Therefore, hybrid genetic algorithm (HGA) has numerous advantages, for instance, global searching, robustness, parallelism. Further, it can be implemented with other algorithms.

8. Hybrid genetic algorithm

Following Ghiglia and Pritt [45], the HGA algorithm relies on estimating the parameters of an nth order of polynomial to approximate the unwrapped surface solution from the wrapped phase data. The coefficients of the polynomial that best unwrap the wrapped phase map are obtained by initial solution of GA algorithm to avoid a long time to converge to the global optimum solution. In this context, GA minimizes minimum 3DBPASL and Q_i,j,k errors between the gradient of the polynomial unwrapped surface solution and the gradient of the original wrapped phase map. In other words, more precision and lower minimum 3DBPASL and Q_i,j,k errors are achieved by increasing the order of the polynomial. This proposed algorithm is mainly applicable to adjoining phase distributions (albeit with gaps). Any optimization problem using a GA requires the problem to be coded into GA syntax form, which is the chromosome form. In this problem, the chromosome consists of a number of genes where every gene corresponds to a coefficient in the nth order of surface-fitting polynomial as described into Eq. (13):

$f:=n\to \sum _{k=0}^n\sum _{j=0}^n\sum _{i=0}^n{a}_{i,j,k}\Delta {\phi }_{i,j,k}^x\Delta {\phi }_{i,j,k}^y\Delta {\phi }_{i,j,k}^z$E13

where a[0.... n] is the parameter coefficients which are retrieved by the genetic algorithm to approximate the unwrapped phase that can achieve the minimum 3DBPASL and Q_i,j,k errors. Further, i, j and k are indices of the pixel location in the unwrapped phase, respectively, and n is the number of coefficients.

The initial population is generated by creating an initial solution using one of the quality-guided phase unwrapping algorithms (3DBPASL algorithm) [40]. Following Ghiglia and Pritt [45], the initial solution is approximated using a ‘polynomial surface-fitting-weighted least squares multiple regression’ method. The initial population is then generated based on the initial solution. In doing so, in every a_g in each chromosome in the population, a small number relying on the accuracy of the gene is added or subtracted to the value of the gene as given by [44]

$a_g=a_g+\left(\pm 1\right)\left\{{10}^{\left[log\left(a_g\right)+\Re \right]}\right\}$E14

where a_g is the coefficient parameter stored in gene g and ℜ is a random number generated between the values.

Calculate the objective values of chromosomes in the population, and record the Pareto-optimal solutions.

Let Q₀, Q₁, Q₂ ∈ Ρ and Ρ is an achievable locality. And, Q₀ is called the Pareto-optimal finding of Q₀ in the minimization problematic decision which is taken under the following satisfied circumstances.

1. If f(Q₁) is supposed to be moderately larger than f(Q₂), i.e. f_m(Q₁) ≥ f_m(Q₂), ∀ i = 1, 2, … , n and

$f_m\left(Q_1\right)>f_i\left(Q_2\right),\exists i=1,2,\dots ,n\text{,}$E15

then Q₁ is assumed to be subjugated by Q₂.

2. Suppose there is no Q ∈ P St. Q dominates Q₀, then Q₀ is the optimal solution of Pareto.

In this step, the Q quality solution is evaluated at every generation to determine the global optimum solution to the parameter estimation-phase unwrapping problem. Therefore, the genes of a chosen chromosome are substituted as coefficients in Eq. (15) to evaluate the approximated phase value at coordinate (i, j, k). Then, the obtained phase is subtracted from the contiguous pixel approximated phase value to retrieve the approximated unwrapped phase solution gradient. It is then subtracted from the gradient of the wrapped phase in the i, j and k directions [39, 44].

Following Hai and Renbiao [56], the two-point greedy continuous crossover is implemented in crossover operator. Therefore, crossover is less problem than the mutation operator. Thus, mutation operator concerns deliberate changes to a gene at random, to keep variation in genes and to increase the probability of not falling into a local minimum solution [41–44]. It involves exploring the search space for new better solutions. This proposed operator uses a greedy technique which ensures that only the best-fit chromosome is allowed to propagate to the next generation.

The accurate 3D phase unwrapping [39] is obtained by phase-matching algorithm proposed by ESA [55]. According to ESA [55], phase-matching algorithm is matched the phase of wrapped phase with unwrapped phase by the given equation:

${\psi }_{i,j,k}=\Delta {\phi }_{i,j,k}+2\pi \rho \left[\frac{1}{2\pi }\left({\widehat{\Delta \phi }}_{i,j,k}-\Delta {\phi }_{i,j,k}\right)\right]$E16

where ψ_i,j,k is the phase-matched unwrapped phase; i, j and k are the pixel positions in the quality phase map; Δϕ_i,j,k is the given wrapped phase; $\Delta {\widehat{\varphi }}_{i,j,k}$ is the approximated unwrapped phase ρ[.] is a rounding function which is defined by ρ[t] = ⌊t + 1/2⌋ for t ≥ 0 and ρ[t] = ⌊t − 1/2⌋ for t < 0 and i, j and k are the pixel positions in x, y and z directions, respectively [41–44].

9. Results and discussion

Figure 15 shows that the Sentinel-1A data were acquired pre-earthquake and post-earthquake on April 17 and 29, 2015, respectively. The urban zones and top of mountains are dominated with higher coherence of 0.8 and 1, respectively, as compared to vegetation and water. In contrast, the water has lower backscatter and coherence of −30 dB and 0.25, respectively (Figure 15b). In fact, Sentinel-1A beam mode of interferometric wide swath has spatial resolution of 5 × 20 m and swath width of 250 km with VV polarization.

Figure 15.

Sentinel-1A satellite data (a) pre- and (b) post-earthquake.

Sentinel-1 is the European Radar Observatory, representing the first new space component of the global monitoring for environment and security (GMES) satellite family, designed and developed by ESA and funded by the European Commission (EC). Sentinel-1 is composed of a constellation of two satellites, Sentinel-1A and Sentinel-1B, sharing the same orbital plane with a 180° orbital phasing difference [57]. In fact, the appropriate SAR polarization for earthquake research either VV or HH polarization modes in C-band SAR satellite data [57].

The specific needs of the four different measurement modes with respect to antenna agility require the implementation of an active phased array antenna. For each swath the antenna has to be configured to generate a beam with fixed azimuth and elevation pointing. Appropriate elevation beamforming has to be applied for range ambiguity suppression. In addition, the incident angle is ranged between 20° and 46°.

Figure 16 shows the interferogram produced by Marghany [57] by conjoining the two complex satellite Sentinel-1A SAR data which are acquired on April 17 and 29, 2015. Figure 16 illustrates huge ground deformation by close cycle of fringe pattern which stricken due to the April 25, 2015 earthquake that smashed into Kathmandu. Furthermore, a few centimeters of ground deformation was uplifted vertically and subsided in Northern Kathmandu. Approximately, an area of 12,000 km² has been drifted also in Northern Kathmandu. In addition, 3 cm of ground deformation is represented with every fringe colour (Figure 16).

Figure 16.

Sentinel-1A interferogram over Kathmandu produced by Marghany [57].

The interferogram fringes are produced by using the combination of three-dimensional best-path avoiding singularity loop (3DBPASL) algorithm and HGA algorithm. Clearly, the proposed algorithm for 3D phase unwrapping produced vibrant fringe cycles which indicate critical surface motion of 8.5 cm which is coincided with ground motion of 1.4 m north of Kathmandu (Figure 17). It is obvious that more than 1 m ground distortion is indicated by the highest concentration of fringe patterns. Moreover, the 3DBPASL shows a sharpest fringe patterns on the east-west of Kathmandu which confirms the dynamic uplift and ground subsiding across Kathmandu with an approximate value of 1.4 m.

Figure 17.

Interferometry produced by hybrid genetic algorithm.

This study confirms the work done by ENVISAT [19]. In fact, the 3DBPASL algorithm acquires an optimal unwrapping path, whereas the effect of singularity loops is also taken into account. In addition, zero-weighted edge is used in zero-weighted edges to adjust the optimal path and avoid these singularity loops [58].

In line with Karout [39], the quality map of phase unwrapping is precisely constructed by implementing algorithm of 3DBPASL. In fact, the proposed algorithm of 3DBPASL can determine the exact singularity loops and also remove the influences of distinctiveness loop uncertainties. Therefore, a combination of 3DBPASL for phase unwrapping with hybrid genetic algorithm produced more precisely fringe cycle. In this regard, hybrid genetic algorithm matches the phase of the wrapped phase with approximated unwrapped phase to establish the best representation of the unwrapped phase. This confirms the work done by Mughier et al. and Rajghatta [59, 60]. With this regard, a genetic algorithm is used to estimate the coefficient of an nth order of polynomial that best approximates the unwrapped phase map which minimizes the difference between the unwrapped phase gradient and the wrapped phase gradient. The genetic algorithm uses an initial solution to speed convergence. The initial solution is achieved by unwrapping using a simple unwrapping algorithm and estimating the parameters of the polynomial using weighted least squares multiple regression [44].

10. Conclusions

The deadly Nepal’s earthquake was caused by an abrupt push, or due to massive stress which was built-up sideways, the foremost fault line somewhere the Indian Plate, ringing India, is gradually plunging beneath the Eurasian Plate, loud abundant of Europe and Asia. Kathmandu, consequently, located on a block of plate about 120 km wide and 60 km stretched, allegedly moved 3 m to the south in just 30 sec. The study has demonstrated new approach for monitoring the Earth deformation using advanced synthetic aperture radar sensor of Sentinel-1 A. The conventional method of InSAR based on 2D phase unwrapping has been modified to 3D phase unwrapping. This study has implemented 3D phase unwrapping based on hybrid genetic algorithm for three-dimensional Nepal’s earthquake construction. In doing so, Sentinel-1A satellite data beam mode of interferometric wide swath during April 17 and 29, 2015 are used. Further, three-dimensional phase unwrapping is performed using three-dimensional best-path avoiding singularity loops (3DBPASL) algorithm. Then, phase matching is implemented with 3DBPASL using hybrid genetic algorithm (HGA). Further, polynomial surface-fitting-weighted least-square multiple regression method is implemented too as initial solution. The study shows that the top of mountain and urban areas have the highest value of backscatter of −5 dB than the surrounding environment. The study shows also that the interferogram produced by ESA has unclear pattern with 3 cm ground deformation and more than 1 m ground motion. However, 3DBPASL algorithm based on HGA provided the sharpest fringe pattern of the interferogram than ESA interferogram. In addition, the complete cycle of fringe patterns indicated 8.5 cm of the surface motion which was coincided with ground motion of 1.4 m north of Kathmandu. In conclusion, integration of the HGA with 3DBPASL phase unwrapping produces excellent 3D Nepal’s earthquake surface displacement using Sentinel-1 A satellite data beam mode of interferometric. The Sentinel-1A satellite data have a great promise for disaster monitoring specially earthquake occurrence.