Advances in three-dimensional deformation mapping from satellite radar observations: application to the 2003 Bam earthquake

ABSTRACT The advance of the surface deformation measurement from synthetic aperture radar interferometry (InSAR) provides an opportunity of reinterpretation for the past controversial event in the geoscience community. According to the development of the multiple-aperture interferometry (MAI) method, three-dimensional (3D) surface displacements can be estimated with few centimetres accuracy by integrating InSAR and MAI observations. In this study, we provided a renewed fault model of the 2003 Bam earthquake using the advanced method. The 3D deformation map showed the clear distribution pattern of the right-lateral strike-slip fault as well as the additional information of an asymmetry of the surface deformation. To determine the optimal model parameters, we employed a two-segment fault model considering the multiple segments. As a result, two sub-parallel fault segments showing N-S trend were obtained. The model parameters of the second segment have relatively large uncertainties though, the first segment which is presumed as the causative fault of the Bam event has been well-modelled with precise model parameters. The more constrained fault model based on the 3D deformation field enabled us to suggest a possibility of a new interpretation and the better understanding of the fault behaviour.


Introduction
It is well known that interferometric synthetic aperture radar (InSAR) has been widely used for mapping the Earth's surface deformation (Massonnet et al. 1993), but it is limited to onedimensional (1D) deformation mapping in the line-of-sight (LOS) direction. That means deformation mapping from the InSAR measurement cannot be available in the along-track direction, which in turn we cannot retrieve the precise 3D movements of the surface. On the other hand, for the better understanding of geological processes of natural hazards such as earthquakes and volcanic activities, three-dimensional (3D) surface deformation mapping has been becoming increasingly more essential.
The extension to two-dimensional (2D) deformation mapping in the LOS and azimuth directions has been accomplished by the integration of InSAR and multiple-aperture SAR interferometry (MAI) (Bechor and Zebker 2006). MAI is azimuth split-band interferometry, which is a technique where the Doppler spectrum is split into two (Bechor and Zebker 2006;Jung et al. 2009). From the two split spectra, forward-and backward-looking interferograms are generated and then the CONTACT Hyung-Sup Jung hsjung@uos.ac.kr azimuth surface deformation is estimated from the phase difference between the two sub-aperture interferograms. The measurement accuracy of the MAI method is lower than that of the InSAR method due to the smaller Doppler bandwidth (approximately 1600-6000 Hz). However, since the MAI measurement is not affected by the atmospheric phase delay due to water vapour (Bechor and Zebker 2006), the performance of the MAI method can be comparable with that of the InSAR method in the coherent areas. The 3D surface deformation field can be retrieved by integrating LOS and along-track deformation vectors measured from ascending and descending observations (Wright et al. 2004;Jung et al. 2011). Since most radar satellites have the near-polar orbit, east and up components are dominant in the LOS deformation vectors while along-track deformation is mostly induced by north component of the surface deformation. Thus, the along-track deformation from the MAI method plays a key role for retrieving 3D maps. Even though the offset tracking method has been used for measuring along-track displacements, the accuracy of the offset-derived azimuth deformation is much lower than that of the LOS deformation by InSAR, and hence the north component in 3D deformation is quite noisy. The MAI method, on the other hand, has improved the measurement accuracy almost doubled when comparing to the offset tracking method (Bechor and Zebker 2006;Jung et al. 2009;Jung et al. 2011;Jung et al. 2013;Jung et al. 2014;Jo, Jung, Won, Poland, Miklius, and Lu 2015;Jung, Yun, and Jo 2015). It was reported that the measurement accuracy by this approach was about 1.6, 3.6 and 2.1 cm for L-band ALOS PALSAR ) and about 0.86, 1.04 and 0.55 cm for X-band COSMO-SkyMed data in the east, north and up directions, respectively (Jo, Jung, Won, and Lundgren 2015). Thus, the advance in the 3D surface deformation measurement could be achieved by the MAI method.
The advanced 3D measurement has been applied to the seismic analysis (Hu et al. 2012;Wang et al. 2014;Boncori et al. 2015;Jo, Jung, Won, Poland, Miklius, and Lu 2015;Wang et al. 2015;Grandin et al. 2016), volcanic activity Jo, Jung, Won, and Lundgren 2015;Jo, Jung, Won, Poland, Miklius, and Lu 2015), etc. Retrieved 3D map will help us to better understand the geological processes of the lithosphere. In other words, the more constrained geological model could be obtained from the advanced 3D deformation field. Especially for the past geological events, we would be able to obtain the meaningful interpretations through the current technique.
In this study, we have retrieved the 3D surface displacement map by using Envisat ASAR ascending and descending images, and performed the 3D-based modelling about the Bam earthquake occurred on 26 December 2003. Although many previous studies have suggested the various scenarios of the Bam earthquake, we expect that the 3D measurements are able to provide a little different but meaningful interpretation about the fault behaviour.

Methods
While the InSAR method has been widely used for measuring surface displacements, the MAI method was not applied in many cases because it requires the high interferometric coherence. In recent studies, however, several researches showed the remarkable improvement of the MAI performance by correcting the phase distortions and residual errors (Jung et al. 2013;Jung et al. 2014;Liu, Jung and Lu 2014). This study provides the precise 3D deformation map for the 2003 Bam earthquake through the advanced observation method.

MAI basics
The MAI method has been developed based on the split-beam InSAR method. The conventional InSAR method utilizes a single-look complex (SLC) image generated from full-aperture signals, but the MAI method uses a SLC image from the sub-aperture processing. During the MAI processing, forward-looking and a backward-looking SLC images are generated by controlling Doppler centroids for each master and slave data. The MAI phase, which can be represented as the alongtrack deformation, is then calculated by the phase differences between forward-looking and backward-looking InSAR phases. The MAI phase (f MAI ) can be defined as follows: where x and l indicate the along-track displacement and the effective antenna length, respectively. A normalized squint represented as n is the faction factor of the full aperture bandwidth (Bechor and Zebker 2006). For ENVISAT ASAR with l = 10.0 m, 2p of the MAI phase represents 10 m of the along-track displacement in case of n = 0.5. Thus, a sensitivity of the MAI measurements is much lower than InSAR measurements because 2p of InSAR measurement represents about 2.8 cm deformation in LOS direction. The MAI measurement, however, is highly required for yielding a precise north component of the surface deformation. Eventually, this will help us to understand the surface motion three-dimensionally.

Retrieval of 3D deformation field
The 3D surface deformation can be reconstructed by integrating InSAR and MAI measurements which obtained at least two different geometries, commonly one ascending and one descending orbits. The imaging geometry for 3D deformation retrieval from LOS and along-track deformation vectors is shown in Figure 1. The observation vector (d) by using the InSAR and MAI method can be defined by (Wright et al. 2004): where d can be defined as where d asc;LOS and d dsc;LOS are the InSAR measurements obtained from the ascending and descending orbits, respectively. In the same way, d asc;AZ and d dsc;AZ are the MAI measurements from the ascending and descending orbits, respectively. u is the 3D components of the surface deformation as given by: where u e , u n and u u are the east, north and up components of the surface deformation, respectively. Design matrix (p) consists of the unit vectors for the InSAR and MAI observations as given by The length of unit vectors and observations are the same. 3D components of surface deformation u can be estimated by solving a least-square inversion problem as follows: where weighting matrix w is diagonal matrix, and can be normally estimated from the standard deviation or coherence of InSAR and MAI measurements. We also calculated phase standard deviations for both InSAR and MAI measurements from ascending and descending data (Jung et al. 2009). All the standard deviation maps were used in the 3D decomposition process as weighting factors. Moreover, we generated error propagation map based on the weighting factors, which in turn was used as the masking window in modelling process.

Study area
The city of Bam is located in the Kerman province of south-eastern Iran, about 900 km south-east of Teheran, Iranian capital. It is situated on the south-western margin of the Dasht-e Lut desert, which has two N-S strike-slip fault systems at the western and eastern sides that are Nayband-Gowk-Sabzevaran (NGS) and Neh strike-slip fault system, respectively ( Figure 2). These two active fault systems accommodate about 14 mm/year relative motion between central Iran and western Afghanistan (Vernant et al. 2004). Bam city is located near the south-eastern termination of the Gowk strike-slip fault which has been seismically active during the past couple of decades. The nearest active fault known by geological studies is the Bam fault (Berberian 1976;Hessami et al. 2003). This fault, located on 5 km south-east of the city, mainly has a right-lateral strike-slip motion and pass through between the Bam and Baravat city ( Figure 2). There were no records of historical earthquakes associated with the Bam fault within the last hundred years. The 2003 Bam earthquake (Mw 6.6), one of the most catastrophic earthquakes in the Iranian recorded history, occurred in the city of Bam, at 1:56 (UTC) on 26 December. The epicentre of this earthquake was located at 28.85 N and 58. 25 E in 10 km southwest of Bam, and its depth was 15 km (Bond ar et al. 2015). The focal mechanism solution indicated that this earthquake had a steep dipping, right lateral strike-slip faulting with an N-S trend (Eshghi and Zar e 2004). Numerous researchers have explored this seismic event in many different ways including the InSAR method  (Talebian et al. 2004;Wang et al. 2004;Fielding et al. 2005;Funning et al. 2005;Fielding et al. 2009). Many of the previous studies suggested that the Bam earthquake was caused by an unknown blind fault. This causative fault showing N-S trend lies at approximately 5 km west of the geological preexisting Bam fault ( Figure 2). Additionally, the deep part of the Bam fault also affected this seismic event. Since the mechanism of the Bam earthquake is quite complex, there were various interpretations and arguments about clarifying the causative fault and its behaviour.

Data processing
This study shows the 3D surface displacements induced by the 2003 Bam earthquake by using coseismic ascending and descending InSAR pairs obtained from the ENVISAT ASAR sensor. Table 1 shows the characteristics of a couple of interferometric pairs used in this study.
For reconstructing the 3D surface deformation, we measured the LOS and along-track displacements from ascending and descending interferometric pairs applying the InSAR and MAI methods. The MAI method first produced four sub-aperture SLC images from a couple of raw SAR signals by controlling the antenna look direction. The Doppler centroids of the forward-and backward-looking SLC images are determined by the effective Doppler bandwidth and a normalized squint (Jung et al. 2009). From the ascending pair, Doppler centroids of ¡556.9, ¡177.6 and ¡950.6 Hz were calculated for the centre-ward, forward and backward directions, and the sub-aperture processing bandwidth of 773.0 Hz was used. Similarly, the MAI processing of the descending pair was carried out by using the Doppler centroids of 28.0, 406.5 and ¡366.7 Hz. The sub-aperture processing bandwidth of 773.2 Hz was used for the ascending pair as well. The normalized quint of 0.5 and the effective antenna length of 10.0 m have been adopted for the MAI processing. After generating SLC images, forward-and backward-looking interferograms were generated, and the MAI interferogram was then produced by calculating phase differences between forward-and backward-looking interferograms. The initial MAI interferograms were enhanced by correcting the topographic phase distortions and removing the residual flat-Earth phases (Jung et al. 2009). The enhanced MAI method enabled us to decrease the phase standard deviation between 30% and 45% in coherent areas. The topographic phase distortions were not much because the perpendicular baselines were small enough, which were about 49.1 and 1.6 m for the ascending and descending pairs.
During the MAI and InSAR processing, two-step multi-looking operations were carried out. We performed 1 £ 5 multi-looking in the range and azimuth directions before the correction of the flat-Earth and topographic phase distortions. The second multi-looking of four looks respectively in the range and azimuth directions was operated. The pixel spacing of the final MAI and InSAR interferograms was approximately 80 m in both range and azimuth directions. For reducing the phase noises, the Goldstein adaptive filter with a window size of 32 was applied (Goldstein and Werner 1998). The non-local (NL) filter was especially employed because we were able to expect higher efficiency of edge detection (Buades et al. 2005). It was helpful to observe the discontinuous boundaries clearly including the fault ruptures.
The topographic phase of InSAR and MAI interferograms was removed by using 1-arc shuttle radar topography mission (SRTM) digital elevation model (DEM). This height map was interpolated with pixel spacing of 20 and 90 m for correcting the first and second multi-looking interferograms, respectively. Despite a small incidence angle about 20 of the ENVISAT SAR system, we were able to remove the topographic phases properly due to the very short perpendicular baselines and flat surface conditions of the study area.
After measuring the LOS and along-track displacements from two different tracks, the 3D components of surface deformation have been retrieved by integrating all the measurements. The 3D mapping was carried out only for the common area from both data-sets, especially around the vicinity of the epicentre (dashed-line box in Figure 2). The standard deviations of the MAI and InSAR phases were calculated (Jung et al. 2009), and the reverse value of square standard deviation was used as the weightings in the 3D decomposition. From the retrieved 3D maps, decorrelated areas below the criteria were masked during the modelling process due to their low reliability. Figure 3 shows the wrapped InSAR and MAI interferograms generated from ascending and descending data-sets. We can clearly identify the deformation centre by the rapid phase gradient of InSAR interferograms. Moreover, its opposite directional phase patterns mean that the fault slips were dominantly occurred in the horizontal direction. The perpendicular baseline of the forward-looking interferogram and the perpendicular baseline difference.

3D deformation fields
The MAI measurements from the ascending and descending paths also show the large horizontal movements in N-S direction (Figure 3(c,d)). We can see the clear rupture line of the Bam earthquake from the MAI measurements, and we are able to recognize the right-lateral shear sense of the fault by the moving directions in the along-track direction. In addition, we have found that the right-hand side of the fault shows the relatively large deformation compared to the left-hand side, especially from the descending path data (Figure 3(d)). This means that the fault plane possibly has been tilted to the east.
We have generated the 3D surface deformation maps by combining LOS and along-track displacements measured from two different geometries. Figure 4(a-c) displays the retrieved 3D components of surface deformation for the Bam earthquake. The east (u e ), north (u n ) and up (u u ) components showed the typical pattern of the right-lateral strike-slip motion. However, we can also see the asymmetric distributions of surface deformation in all the 3D component maps (Figure 4(ac)). Especially, the southern parts of the fault zone showed higher asymmetry than the northern parts. This asymmetry can be seen more clearly by the profiles of surface displacements across the rupture zone ( Figure 5). While the profile of A-A' on the 3D maps shows the similar level of surface displacements in both left and right sides of the rupture, the noticeable difference of surface displacements can be seen in the southern areas of the fault zone by the profile B-B'. The right side of the fault shows the surface displacements at least two times larger than those of left side. These results could be interpreted that the fault plane is a little tilted to the eastward direction.
While the surface deformation by the Bam earthquake showed quite typical pattern of a dextral strike-slip fault, we have also recognized that this event occurred by more than a single fault. Several previous studies also suggested a fault model consisting of two or three segments for the Bam earthquake (Wang et al. 2004;Xia 2005;Fielding et al. 2009). In this study, we confirmed three segments by the SAR observation. There were some ambiguities to determine the positions of the surface rupture from the InSAR interferograms, but the 3D maps enabled us to confirm the clear surface  (Figure 4). The southern and the northern segments have the similar strike, and relatively short segment connects the other two segments. This middle segment is placed on the collapsed Bam city, so it does not show a clear fault line. Based on the 3D observation and surface ruptures, it could be interpreted that two north-south direction faults were slipped and the city zone between active faults were destructed by a pull-apart extension.

3D-based fault modelling
There have been a number of studies to model the Bam earthquake based on the InSAR measurements or GPS observations. As was interpreted from the 3D measurements, it was strongly presumed that this seismic event was triggered by at least two adjacent and sub-parallel faults. A few studies performed modelling assuming a single fault model, but the modelled results showed misfits with large residuals (Stramondo et al. 2005;Peyret et al. 2007). Thus, many studies carried out twofault modelling by using InSAR measurements or complex observations of InSAR and azimuth offset measurements. Although these studies provided more probable fault models of the Bam earthquake, there were various interpretations about the causative fault location and its behaviour (Talebian et al. 2004;Wang et al. 2004;Funning et al. 2005;Fielding et al. 2009). In this study, we modelled the Bam earthquake using the 3D displacement map Jo, Jung, Won, and Lundgren 2015) assuming the two-segment fault model, and suggested the more constrained model of the earthquake. Basically, there is no intrinsic difference between using 3D measurement and using a combination of InSAR and MAI measurements for modelling ). However, it was more helpful to use the 3D map than LOS deformation for determining initial boundary conditions in model parameter estimation, because the 3D map enabled us to interpret surface movements easily. That means the 3D map could provide a better understanding for surface deformation.
We employed the Okada dislocation model under isotropic, homogeneous and elastic half-space conditions for determining the fault geometry and slip amounts induced by the Bam earthquake (Okada 1985(Okada , 1992. A Poisson's ratio of 0.25 was adopted in this modelling. Among the three dislocation components, only slip-related parameters were estimated assuming the opening was negligible. Since we confirmed a possibility of complex movements of multiple faults, we applied the dislocation model for two finite fault planes which lies on the north-south direction. The middle segment showing the NE-SW trend was not considered in the modelling because we could fairly perceive that the surface displacement patterns were formed by N-S trend dextral shear. To determine the best-fit parameters of the two-segment fault model, we performed Monte Carlo simulation with 10,000 iterations. Figure 6 displays the inversion results of the best-fitting dislocation model for the Bam earthquake, constrained by the 3D component maps. The east, north and up components of the surface deformation were jointly modelled to estimate the fault dimension, geometry and slip amounts. The east and north components were well-modelled with less residual, while the up component showed relatively large residuals. The total root-mean-square ratios between the residuals and measurements were 0.25, 0.23 and 0.60 for the east, north and up components, respectively. Residuals were mostly remained near the first segment from all 3D components. We were able to see the clear residual pattern by the underestimation from the modelled result of east component. This could be interpreted that the estimated model could not detect the slips on the shallow-depth zone. We also considered a feasibility that the undetectable deformation has been induced by the other source such as a small movement in shallow areas. Table 2 provides mean and standard deviation of the model parameters after removing the outliers from each parameter's histogram. Strikes of the optimal uniform-dislocation faults were  estimated as 3568 and 8.78, respectively. As can be seen from Figure 6, two faults showed an almost parallel N-S trend, but a little tilted to the opposite directions. We estimated a nearly vertical fault plane for the first segment which is lying in the South of Bam city. For the second segment, approximately 608 and 98.78 of dip and dip-direction were determined in this modelling. While the standard deviation of the first segment's model parameters was small enough, relatively large uncertainties of the model parameters were estimated for the second segment. This result coincided with the fact that it was difficult to find the clear fault trace of the second segment from the 3D surface displacement map. Based on the best-fit model parameters, the depths of the centre and the top part of the second fault segment was about 6 and 4 km depth, respectively, from the surface. We also presumed that two fault segments had a quite different slip amount in strike direction by the model, while dip-slip of the two segments was quite similar. Strike-slip of the main segment was approximately ¡2.8 m, and it was more than two timeslarger than strike-slip of the second one. Moreover, the slip centre of the first segment was about 4.1 km, which was almost 2 km shallower than the second one. From these modelled results, we could get an interpretation that the surface deformation was mostly induced by the first segment, and the second one made a distorted or asymmetrical pattern of the overall surface deformation. It is regarded that the second segment contributed to make complex surface deformation, but did not generate a dominant surface pattern of its own. The 3D map was not able to validate its accuracy due to the absence of in situ measurements, but it has a determinate meaning for providing a possibility of a new fault model. Especially, the 3D-based modelling is encouraging in the point that the precise model parameters were obtained for the first segment which is strongly presumed as the causative fault of the Bam earthquake. These findings are expected to allow the different interpretation by providing a more constrained fault model.

Conclusions
We have retrieved the 3D surface displacement maps for the 2003 Bam earthquake, by combining the LOS and along-track deformation from ascending and descending ENVISAT ASAR images. Compared to the previous studies, the 3D deformation map enabled us to identify the more realistic surface movements by the earthquake. The 3D deformation map has shown a clear surface deformation pattern induced by a right-lateral strike-slip fault. One notable point was that the southern part of the fault zone showed a high asymmetry from the 3D maps. It means the southern part of the fault could be affected by a low-dipping slip motion. Moreover, the 3D maps clearly demonstrated that this seismic event has been triggered by multiple fault segments.
Based on the 3D observation, we have determined the two-segment fault model of the Bam earthquake by employing the Okada dislocation model. The model parameters for the two sub-parallel segments showing N-S trend were obtained. As a result of iterative estimation by Monte Carlo simulation, the second fault segment has more than doubled uncertainties compared to the first segment. This result coincides with the fact that we were not able to observe the clear surface trace of the second segment. However, these findings enable us to estimate more constrained model parameters of the first fault segment which is strongly presumed as a causative fault. In conclusion, the 3D-based modelling allowed us to interpret this seismic event with more constrained model, and provided the new understanding for the mechanism and behaviour of the fault behaviour.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This study was funded by the Korea Meteorological Administration Research and Development Program [grant number KMI2017-9060].