An integrated method based on DInSAR, MAI and displacement gradient tensor for mapping the 3D coseismic deformation field related to the 2011 Tarlay earthquake (Myanmar)

Remote Sensing of Environment 170 (2015) 388–404

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

An integrated method based on DInSAR, MAI and displacement gradient tensor for mapping the 3D coseismic deformation field related to the 2011 Tarlay earthquake (Myanmar) Xiaowen Wang a, Guoxiang Liu a,b,⁎, Bing Yu a, Keren Dai a, Rui Zhang a, Deying Ma a, Zhilin Li a,b a b

Department of Remote Sensing and Geospatial Information Engineering, Southwest Jiaotong University, Chengdu, 610031, Sichuan, China State-Province Joint Engineering Laboratory of Spatial Information Technology for High-Speed Rail Safety, Chengdu 610031, Sichuan, China

a r t i c l e

i n f o

Article history: Received 23 March 2015 Received in revised form 14 September 2015 Accepted 25 September 2015 Available online 22 October 2015 Keywords: 3D coseismic deformation field DInSAR MAI Displacement gradient tensor

a b s t r a c t The satellite differential interferometric synthetic aperture radar (DInSAR) technology has been widely applied for mapping the ground deformation associated with the geophysical events such as earthquakes. However, the conventional DInSAR can measure only the one-dimensional (1D) ground displacement along the radar line of sight (LOS), and the crucial displacement measurements, e.g., near the epicenter and along the surface ruptures, are usually not available or degraded in quality due to the significant deformation gradients. With availability of satellite SAR images acquired along ascending and descending orbits, this paper proposes an integrated method to map the three-dimensional (3D) coseismic deformation field by combining DInSAR for detecting LOS displacements, the multiple aperture interferometry (MAI) for detecting along-track displacements, and the displacement gradient tensor (DGT) model for characterizing spatial correlation of ground displacements. The proposed method (termed as InSAR–DGT) was first tested through an experiment of recovering 3D displacements from a simulated coseismic deformation field. We then applied the proposed method to map the 3D coseismic deformation field related to the 2011 Tarlay earthquake (Myanmar) by using the ascending and descending ALOS PALSAR images. Both the results derived for the simulated experiment and the 2011 Tarlay earthquake show that the quality of the 3D coseismic displacements can be raised efficiently by the InSAR–DGT method, and the precisions in the east–west (E–W), north–south (N–S) and up–down (U–D) displacements for the 2011 Tarlay earthquake are increased 22%, 36% and 24%, respectively. The validation indicates that the improved coseismic deformation field for the 2011 Tarlay earthquake is in good agreement with the deformation field simulated from the existing optimized fault model. The number of the missing data points in the 3D coseismic deformation field can be reduced significantly by the InSAR–DGT method. It is revealed that the causative fault (i.e., the west part of the Nam Ma Fault) for the Tarlay earthquake generated a left-lateral slip that was accompanied by a minor normal dip-slip component. The careful interpretation demonstrates that the causative fault structures include the determined main fault segment and the two suspected branched segments at east of the Tarlay town. © 2015 Elsevier Inc. All rights reserved.

1. Introduction The satellite differential interferometric synthetic aperture radar (DInSAR) technology has proven to be very useful and potential for monitoring the regional ground deformation associated with the geophysical events such as earthquakes, volcanoes, landslides, ice flows and land subsidence (e.g., Massonnet & Feigl, 1998; Mohammed, Saeidi, Pradhan, & Yusuf, 2013; Prati, Ferretti, & Perissin, 2010; Rucci, Ferretti, Monti, & Rocca, 2012; Sansosti et al., 2014). However, one of the limiting factors in the conventional DInSAR is that it can measure ⁎ Corresponding author at: Department of Remote Sensing and Geospatial Information Engineering, Southwest Jiaotong University, Chengdu 610031, Sichuan, China. E-mail address: [email protected] (G. Liu).

http://dx.doi.org/10.1016/j.rse.2015.09.024 0034-4257/© 2015 Elsevier Inc. All rights reserved.

only one-dimensional (1D) ground displacement along the radar line of sight (LOS), thus resulting in biased source parameters for geophysical modeling (González, Fernández, & Camacho, 2009; Jung, Lee, & Zhang, 2014). In addition, the crucial displacement measurements, e.g., near the epicenter and along the surface ruptures, are usually not available or degraded in quality due to the significant deformation gradients. In recent years, the great efforts have been made to extend the capability of DInSAR to measure the three-dimensional (3D) ground displacements particularly for the coseismic deformation analysis (Hu et al., 2014). The most widely used method for mapping a 3D deformation field is based on a liner inversion model by combining the LOS and along-track (AT) displacements derived from the multi-track satellite SAR images (Erten, Reigber, & Hellwich, 2010; Fialko, Simons, &

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

Agnew, 2001; Hu et al., 2014; Jung, Lu, Won, & Miklius, 2011; Raucoules, de Michele, Malet, & Ulrich, 2013; Wang et al., 2014a). The existing studies show that the AT displacements can be obtained by either the pixel offset tracking (POT) method (Michel & Avouac, 1999) or the multiple aperture interferometry (MAI) method (Barbot, Hamiel, & Fialko, 2008; Bechor & Zebker, 2006; Jung, Won, & Kim, 2009). As the MAI and POT methods concentrate on analysis of interferometric phases and SAR amplitude information, respectively, the AT-displacement measurements derived by MAI are generally more precise than those derived by POT (Jung et al., 2009; Wang et al., 2014a). However, due to decreasing of the azimuth bandwidth, the quality of MAI interferograms may be degraded with decorrelation noise. It should be noted that the POT method can also be applied to the optical images to extract horizontal displacements, which can be combined with InSAR measurements for extracting the 3D deformation field (Barnhart, Briggs, Reitman, Gold, & Hayes, 2015; González, Chini, Stramondo, & Fernandez, 2010; Leprince, Barbot, Ayoub, & Avouac, 2007). Although the linear inversion method based on the LOS and AT displacements can be applied to reconstruct a 3D deformation field, the resultant 3D displacements may be degraded in quality by decorrelation noise or contain some missing data points due to the significant deformation gradients. In addition, it should be pointed out that the spatial correlation of ground displacements is not taken into account by the liner inversion method. According to the elastic deformation theory (Segall, 2010), the vectors of coseismic displacements at the adjacent points are spatially correlated in the statistical and elastic senses. The closer the two adjacent points, the more similar the ground displacements in magnitude and direction are. An alternative method for mapping a 3D deformation field is based on a joint use of DInSAR and GPS measurements (Gudmundsson, Sigmundsson, & Carstensen, 2002; Samsonov & Tiampo, 2006). However, the sparse GPS stations cannot match the high density of DInSAR measurements and the interpolation of GPS data has to be performed, thus resulting in unsatisfactory 3D displacements. In order to avoid the direct interpolation procedure, Guglielmino, Nunnari, Puglisi, and Spata (2011b) proposed a SISTEM technique to reconstruct the 3D deformation field. Taking into account the elastic sense of the deformation field (Segall, 2010), the SISTEM method estimates the unknown 3D displacement components at a point of interest by using the GPS or DInSAR measurements at the adjacent points and incorporating the displacement gradient tensor (DGT) parameters of the relevant points. The SISTEM method has been applied to investigate the ground deformation of the Mount Etna area and the 2009 L'Aquila earthquake (Italy) (Guglielmino et al., 2011a, 2011b). Nevertheless, the existing studies indicate that the dense and uniformly distributed GPS stations are required for guaranteeing the quality of SISTEM solution. In reality, such deployment of GPS stations is not available in many study areas. With availability of SAR images acquired along ascending and descending orbits, this paper proposes an integrated method to map the 3D coseismic deformation field by combining the DInSAR processing for detecting LOS displacements, the MAI processing for detecting AT displacements and the DGT model for characterizing spatial correlation of ground displacements, thus further overcoming the disadvantages in the aforementioned methods. The proposed method is termed as InSAR–DGT and also characterized by that the additional datasets such as GPS measurements are not needed, but only the SAR images are used for mapping the 3D coseismic deformation field. Based on the concept of the spatial correlation of ground displacements in the statistical and elastic senses by incorporating the DGT model, the InSAR–DGT method is therefore expected to improve the 3D coseismic deformation results. The InSAR–DGT method first generates the initial 3D coseismic deformation maps in a study area through the linear inversion method by combining the ascending/descending DInSAR and MAI measurements. The 3D displacements are then re-estimated through the weighted LS adjustment by incorporating the DGT model into the initial 3D displacements. For validation purpose, the InSAR–DGT method will

389

be first tested through an experiment of recovering 3D displacements from a simulated coseismic deformation field, and then applied to map the 3D coseismic deformation field related to the 2011 Tarlay earthquake (Myanmar) by using the ascending and descending ALOS PALSAR images. In this paper, the introduction is followed by proposing the InSAR–DGT method for mapping the 3D coseismic deformation field (Section 2). A simulated experiment is presented in Section 3 for testing the applicability of the InSAR–DGT method. Section 4 concentrates on the application of the InSAR–DGT method for mapping the 3D coseismic deformation field related to the Tarlay earthquake. The validation and interpretation of the coseismic displacements are performed in Section 5. The concluding remarks are given in Section 6. 2. The InSAR-DGT method for mapping the 3D coseismic deformation field In this Section, we will first briefly describe a linear inversion method based on the DInSAR and MAI processing, as well as the LS solution for extracting the initial 3D coseismic ground displacements. To improve the 3D coseismic deformation field, the InSAR– DGT method is then presented, which is based on the weighted LS solution by incorporating the DGT model into the initial 3D displacement measurements. 2.1. Extracting the initial 3D coseismic displacements by a linear inversion method Suppose that two interferometric pairs are available over an area affected by an earthquake. One pair of SAR images is acquired along the ascending orbits before and after the event, respectively, while another pair of SAR images is acquired along the descending orbits before and after the event, respectively. The two interferometric pairs can be used to extract the ground displacements. With use of ascending and descending SAR images, several research groups demonstrated that the 3D deformation field could be generated by a linear inversion method based on the DInSAR processing for detecting LOS displacements and the MAI processing for detecting AT displacements (Hu et al., 2014; Jung et al., 2011; Wang et al., 2014a). For better understanding, we will shortly brief the primary procedures of the linear inversion method. For extracting the LOS coseismic displacements from each interferometric pair, the DInSAR processing steps are conducted by using, for example, the ROI_PAC software provided by the Jet Propulsion Laboratory (JPL) (Rosen, Hensley, Peltzer, & Simons, 2004). The primary steps include the generation of raw interferogram, removal of flat-earth trend effects, removal of topographic effects, phase unwrapping and calculation of LOS displacements (e.g., Hanssen, 2001; Liu et al., 2010; Massonnet & Feigl, 1998), thus resulting in the two maps of LOS displacements of the study area corresponding to the ascending and descending cases, respectively. A digital elevation model (DEM) derived from the Shuttle Radar Topography Mission (SRTM) with 3 arcsec resolution (Farr et al., 2007) is generally used to remove the topographic effects from the interferograms. The relationship between the LOS coseismic displacement DLOS and the unwrapped interferometric phase φunw at a pixel can be expressed by DLOS ¼ 

λ ∙φ 4π unw

ð1Þ

where λ denotes the wavelength of the radar sensor (e.g., 23.6 cm for PALSAR). For extracting the AT coseismic displacements from each interferometric pair, we performed the MAI processing procedures (Bechor & Zebker, 2006). As the raw datasets recorded by a radar sensor are actually a collection of backscattering signals in the full synthetic aperture, the Doppler frequency changes regularly with the variation of squint

390

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

angle of the radar illumination (Hanssen, 2001). Depending on the Doppler information, one can split the full aperture into the forward and backward apertures to generate the sub-band single look complex (SLC) images (Barbot et al., 2008). The sub-band SLC images can be generated from the raw datasets by following the method by Jung et al. (2009), thus resulting in the forward- and backward-looking SLC images for each SAR image. For the ascending interferometric pair, we denote Mf and Mb as the forward- and backward-looking parts of the master image, respectively, and Sf and Sb as the corresponding parts of the slave image. Mf and Sf are used to generate the forward-looking interferogram by the conventional differential interferometry method, while Mb and Sb are used to generate the backward-looking interferogram. The ascending MAI interferogram is then constructed by a complex conjugate multiplication of the forward- and backward-looking interferograms. To derive the AT displacements, the flat-earth trend phases and topographic phases should be removed from the ascending MAI interferogram using a polynomial fitting method (Hu et al., 2014; Jung et al., 2009). The relationship between the AT displacement DAT and the final MAI phase φMAI at a pixel can be expressed by DAT ¼

L ∙φ 4π∙n MAI

ð2Þ

where L denotes the antenna length of the radar sensor (e.g., 8.9 m for PALSAR), and n is the fraction (generally 0.5) of the full-aperture width. The descending interferometric pair can be processed in the same way to extract the AT displacements. The previous studies indicate that the AT displacements derived by MAI can reach up to centimeter accuracy for the geophysical applications (Bechor & Zebker, 2006; Gourmelen et al., 2011; Liu et al., 2014). For example, the accuracy of 3–4 cm in the AT displacements can be achieved by the MAI processing with the L-band (long wavelength) ALOS PALSAR images in the case that the interferometric coherence γ is larger than 0.8 (Jung et al., 2014). After the ascending/descending LOS and AT measurements are geocoded and resampled into the same geographic grids (e.g., in the UTM coordinate system), one can generate the 3D deformation maps using the LS method. For an arbitrary point, the relationship between the 3D displacement components and the LOS and AT displacements can be represented by (Fialko et al., 2001; Wright, Parsons, & Lu, 2004) 3 2 A 3 2 A 3 DALOS Γ LOS δLOS 6 D 7 6 D 7 6 D 7 6 DLOS 7 6 Γ LOS 7 6 δLOS 7 6 A 7 ¼ 6 A 7∙d þ 6 A 7 4 DAT 5 4 Γ AT 5 4 δAT 5 DDAT Γ DAT δDAT

on a point-by-point basis, thus resulting in the initial 3D coseismic deformation maps of the study area. As a remark, if only three observaD tions (e.g., DALOS, DD LOS and DAT) are useful, one can calculate the unknown 3D displacement components by a direct solution, instead of the LS adjustment. It should be noted that the spatial correlation of ground displacements are not considered by the linear inversion method, and the missing data points generally appear in some parts of the coseismic deformation maps due to the coherence loss of interferometric phases or the remarkable deformation gradients.

2.2. Improving the 3D coseismic displacements by the InSAR–DGT method In order to better constrain the 3D displacement components derived interferometrically, the previous studies proposed to combine the in situ data (e.g., measured by GPS, leveling and total station) with the InSAR measurements for overcoming the weakness of the method as described in Section 2.1 (Guglielmino et al., 2011b; Muller et al., 2014). Without use of the sparse in situ data, we however re-estimate the 3D coseismic displacement components on a point-by-point basis through the weighted LS adjustment by incorporating the DGT model into the initial 3D displacements derived by the linear inversion method, thus interpolating the missing data points and improving the 3D deformation results. To re-estimate the 3D displacement components at an arbitrary point, the crucial procedures include the selection of valid points (VPs) with high quality of displacement measurements around P, the formation of observation equations based on the DGT concept and the weighted LS adjustment. Fig. 1 schematically illustrates how to select the VPs around an arbitrary point P. We select the VPs around P according to the three thresholds, i.e., correlation distance, interferometric coherence and displacement dispersion. To determine the correlation distance threshold, the experimental variogram should be first identified for the E–W, N–S and U–D displacements, respectively. As the areas far away from the epicenter can be viewed as the zero-deforming parts, we use the farfield “initial displacements” (actually noises) for estimation of the experimental variogram. Previous studies indicate that the characteristics of errors in InSAR measurements may be approximately represented by a power-law function (Hanssen, 2001; Sudhaus & Jonsson, 2009).

2

ð3Þ

where the superscripts A and D stand for the ascending and descending cases, respectively; d is the unknown 3D displacement vector at the point to be estimated and denoted by [de , d n, d u ] Τ (the superscripts e, n and u correspond to the E–W, N–S and U–D dimensions, respectively); δLOS and δAT are the errors in the LOS and AT displacements, respectively; ΓLOS and ΓAT are the transform vectors which can be expressed by       3π 3π ;  sinðθÞ∙ cos φ  ; cosðθÞ Γ LOS ¼  sinðθÞ∙ sin φ  2 2

ð4Þ

      3π 3π ; sin φ  ;0 Γ AT ¼  cos φ  2 2

ð5Þ

where θ is the radar off-nadir angle; φ is the heading angle of the satellite flight direction (measured clockwise from the north). Based on the observation Eq. (3), the unknown 3D displacement components at the point can be resolved by an LS solution (Fialko et al., 2001; Wang et al., 2014a). Such LS adjustment can be performed

Fig. 1. The schematic illustration for selecting the valid points around an arbitrary point P. The valid points are selected from the points within a circle centered at P with a radius of S (i.e., correlation distance threshold). The valid points, invalid points and missing data points are marked by the triangles, filled circles and squares, respectively.

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

Assuming that the errors in the initial displacements are stationary, their characteristics (i.e., variogram) in each of the E–W, N–S and U–D dimensions can be expressed by a 1-D exponential function (Funning, Parsons, Wright, Jackson, & Fielding, 2005; Jolivet et al., 2012; Sudhaus & Jonsson, 2009), i.e., C ðrÞ ¼ C 0  expðr=SÞ

ð6Þ

where r is the distance between the two sampled points; C0 and S are the model parameters representing the variance and the correlation distance, respectively. We estimate C0 and S through data fitting with the far-field initial displacements derived by the linear inversion method as described in Section 2.1. In this way, we can derive three values of correlation distance for the E–W, N–S and U–D dimensions. To balance between the computation cost and the estimation precision, we take the averaged correlation distance S as the final correlation distance threshold. If the distance between the point P and its neighboring point is not greater than S (see Fig. 1), and the interferometric coherence values for both the ascending and descending cases at the neighboring point are not smaller than 0.8, this neighboring point can be identified as a VP candidate. We then calculate the mean values and standard deviations of the E–W, N–S and U–D initial displacements, respectively, of all the VP candidates around P. If the initial displacement in each of the three dimensions of a VP candidate lies smaller than two times of standard deviations away from the mean value (i.e., with low displacement dispersion), the VP candidate is finally identified as the true VP. In this way, all the true VPs around P can be obtained. With availability of a DEM (i.e., the 3D location is known for each point) and a continuous 3D deformation field in a study area, previous studies indicate that the 3D strain tensor on the ground surface can be computed with the DGT model (Cardozo & Allmendinger, 2009; Guglielmino et al., 2011b; Pietrantonio & Riguzzi, 2004; Teza, Pesci, & Galgaro, 2008). To re-estimate the 3D coseismic displacements at P, we incorporate the DGT model into the initial 3D displacements at all the true VPs under the hypothesis of infinitesimal and homogeneous strain (Guglielmino et al., 2011b; Muller et al., 2014). Suppose that the K VPs have been selected around P. The UTM coordinates at P and the ith (i = 1 , 2 , … , K) VP are denoted by ξ0 =[x0, y0,z0]Τ and ξi = [xi,yi,zi]Τ, respectively. The 3D displacement vector to be estimated is denoted by d0 = [de0, dn0, du0]Τ, and the known initial 3D displacement vector at the ith VP is denoted by di = [dei , dni , dui ]Τ. By incorporating the DGT model, the relationship between the 3D displacement vectors at the ith VP and P can be represented by di ¼ H∙Δξi þ d0

ð7Þ

where Δξi is the coordinate-increment vector [Δxi, Δyi,Δzi]Τ and equal to [xi -x0, yi - y0, zi - z0]Τ, and H is the DGT matrix and can be expressed by 3 2 0 ε 11 ε12 ε13 H ¼ E þ Ω ¼ 4 ε12 ε22 ε23 5 þ 4 ω3 ε13 ε23 ε33 ω 2 32 ε 11 ε12  ω3 ε13 þ ω2 ¼ 4 ε12 þ ω3 ε22 ε23  ω1 5 ε13 ω2 ε23 þ ω1 ε33 2

ω3 0 ω1

3 ω2 ω1 5 0 ð8Þ

in which the symmetric (E) and anti-symmetric (Ω) matrix correspond to the infinitesimal strain tensor and rotation tensor, respectively. For all the K VPs, the overall observation equations can be formed according to Eq. (7), i.e., B

∙ X ¼ L þ V

3K12 121

3K1

3K1

ð9Þ

where L and V are the observation and residual vectors, respectively, of the initial 3D displacements at the K VPs; X is the unknown vector of the 3D displacements at P and the DGT parameters to be estimated; and B is

391

the coefficient matrix. The specific elements of the vectors or matrix can be expressed as follows:  e e e n n n u u u Τ L ¼ d1 ; d2 ; ⋯dK ; d1 ; d2 ; ⋯dK ; d1 ; d2 ; ⋯dK ;

ð10Þ

V ¼ ½v1 ; v2 ; ⋯v3K Τ ;

ð11Þ

 e n u Τ X ¼ d0 ; d0 ; d0 ; ε11 ; ε12 ; ε13 ; ε22 ; ε23 ; ε33 ; ω1 ; ω2 ; ω3 ;

ð12Þ

2

1 6⋮ 6 61 6 60 6 B¼6 6⋮ 60 6 60 6 4⋮ 0

0 ⋮ 0 1 ⋮ 1 0 ⋮ 0

0 Δx1 ⋮ ⋮ 0 ΔxK 0 0 ⋮ ⋮ 0 0 1 0 ⋮ ⋮ 1 0

Δy1 ⋮ ΔyK Δx1 ⋮ ΔxK 0 ⋮ 0

Δz1 ⋮ ΔzK 0 ⋮ 0 Δx1 ⋮ ΔxK

0 ⋮ 0 Δy1 ⋮ ΔyK 0 ⋮ 0

0 ⋮ 0 Δz1 ⋮ ΔzK Δy1 ⋮ ΔyK

0 ⋮ 0 0 ⋮ 0 Δz1 ⋮ ΔzK

0 ⋮ 0 Δz1 ⋮ ΔzK Δy1 ⋮ ΔyK

Δz1 ⋮ ΔzK 0 ⋮ 0 Δx1 ⋮ ΔxK

3 Δy1 ⋮ 7 7 ΔyK 7 7 Δx1 7 7 ⋮ 7 7 ΔxK 7 7 0 7 7 ⋮ 5 0

ð13Þ The unknown parameters can be resolved using the weighted LS method, i.e.,  1 X ¼ BΤ WB BT WL

ð14Þ

where W is the weighting matrix with a dimension of 3K × 3K, which reflects the variance–covariance information of the observations. It should be pointed out that we can calculate the elements of W using the empirical covariance function expressed by Eq. (6) by following the method by Jolivet et al. (2012). 2.3. Analysis of error propagation in the InSAR–DGT method It can be observed from Sections 2.1 and 2.2 that the errors in the 3D displacements (derived by the InSAR–DGT method) at an arbitrary point are originated from the InSAR observations and the LS solution by incorporating the DGT model. The InSAR–DGT results are affected not only by the errors in InSAR measurements (i.e., LOS displacements by DInSAR and AT displacements by MAI), but also by the number and distribution of the VPs selected for the weighted LS solution. According to the investigation by Liu et al. (2014), the number and distribution of the VPs directly govern the condition number of the coefficient matrix B in Eq. (9), which is a key factor indicating the robustness (illness) of a linear equation system. The numerical analysis on relation between the number of VPs and the precision of 3D displacements will be presented in Section 5.1 as a case study. Here we concentrate on discussion of how the errors in InSAR measurements might propagate through and bias the InSAR–DGT results. The errors in InSAR measurements are generally caused by several factors, including spatiotemporal decorrelation of SAR signals, uncertainty of orbital data and topographic data (i.e., DEM), atmospheric delay and phase unwrapping. These errors are propagated into the InSAR–DGT results first through the linear inversion procedure (i.e., regularization, see Eq. (3) for estimating the initial 3D displacements) and then trough the weighted LS solution procedure (see Eq. (14) for estimating the final 3D displacements). By following the work by Liu et al. (2014), one can group the errors in InSAR measurements into three categories, i.e., random error (δrand), systematic error (δsys), and gross error (δgross). The random error (e.g., due to decorrelation) always follows a Gaussian distribution, which adds variability to data but does not bias the final results. The systematic error (e.g., due to uncertainty of orbital data and topographic data, or atmospheric delay) can be modeled mathematically and therefore corrected, while the gross error (e.g., due to phase unwrapping) is generally difficult to detect. It should be noted that the gross error is equivalent to the systematic error in a localized area, and both the systematic and gross

392

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

errors should be carefully controlled and handled during the InSAR data processing steps. Otherwise, the InSAR–DGT results might be biased remarkably. According to the law of error propagation (Liu et al., 2014), the errors in InSAR measurements can be propagated into the InSAR–DGT results, i.e.,

δj

error

3K

¼ ∑ Nji ∙δi i¼1

2 rand

!2

3K

þ ∑ Nji ∙δi i¼1

sys

ð15Þ

where i and j are the number of the VPs and estimated parameters, respectively; δj_error is the synthetic variance of the jth parameter estimated by the InSAR–DGT method; δi_rand and δi_sys are the random error and the systematic error for the ith observation, respectively, and N = (BΤWB)-1BTW. As it is difficult to separate the random and systematic errors from the regularized errors, the variance of the 3D displacements derived by the InSAR–DGT method cannot be calculated quantitatively using Eq. (15). By considering that the regularized errors are spatially correlated, we therefore construct the variance–covariance matrix for the VPs, thus assessing the precision of the parameters estimated by InSAR–DGT using the following equation: DX ¼

1 V Τ WV  T B WB f

ð16Þ

where DX is a matrix whose diagonal elements denote the variances of the relevant unknown parameters, f is the redundancy number and equivalent to 3K - 12. 3. A simulated experiment for testing the InSAR–DGT method For testing the InSAR–DGT method proposed in Section 2, an experiment will be performed in this Section with a simulated coseismic deformation field. The experiment intends to make one convinced of both the necessity to implement the InSAR–DGT method and the accuracy of the derived 3D displacements themselves. The dataset will be first simulated with use of the Okada dislocation modelling (Okada, 1992). Both the linear inversion and InSAR–DGT methods are then applied to recover the 3D displacements for the simulated coseismic deformation field. The comparison and analysis will be made to demonstrate if and how the InSAR–DGT method can improve estimation of the 3D displacements and inversion of the fault slip model. 3.1. Dataset simulated from a synthetic fault slip model With use of the Okada dislocation modelling (Okada, 1992), we generated a 3D coseismic deformation field based on a synthetic rectangular fault plane in a homogenous half-space with a passion ratio of 0.25 and a shear model of 30 GPa. Assuming that the simple fault plane is perpendicular to the ground surface and strikes to the north direction, we set the geometric parameters of the fault plane to 10 km in length, 6 km in width and 0.5 km in depth. Assuming that a uniform slip occurs along the fault plane, we set the strike-slip component to 2 m and the dip-slip component to 0.5 m. Fig. 2a–c shows the maps of simulated coseismic displacements for the E–W, N–S and U–D dimensions, respectively, in an area of 40 × 40 km2 (corresponding to 200 × 200 pixels). The fault trace is marked by the white line in Fig. 2a–c. It should be noted that the 3D displacements shown in Fig. 2a–c will be used as the reference data (i.e., true values) for the subsequent comparison analysis. Suppose that the ALOS PALSAR is taken as the imaging radar system, we simulated the InSAR measurements, i.e., LOS and AT displacements by DInSAR and MAI, respectively, through using the nominal imaging parameters of ALOS PALSAR and projecting the simulated 3D displacements (Fig. 2a–c) onto the LOS and AT directions, respectively (see

Eqs. (3)–(5)). We then added the spatially correlated noises (with a correlation distance of 4 km) to the InSAR measurements simulated. The noises for either the DInSAR or MAI measurements were generated using a predefined variance–covariance model with a given variance value characterizing the precision in measurements (Jolivet et al., 2012; Sudhaus & Jonsson, 2009). As the DInSAR measurements are generally more accurate than the MAI measurements, we set the different variance values, i.e., 1 and 4 cm2, to generate noises for the LOS and AT displacements, respectively. For the convenience of subsequent analysis, 40% of all the pixels were randomly selected as the VPs without any added noises; the remaining pixels were the perturbed points (PPs, 58%) with the added noises and the null points (NPs, 2%) with the missing InSAR measurements. Fig. 2d–f shows the maps of noises generated for three cases of InSAR measurements, i.e., the descending LOS and AT displacements and the ascending AT displacements, respectively, while Fig. 2g–i shows the maps of InSAR measurements (with noises) for the same three cases, respectively. It should be noted that each of the NPs is indicated by a black dot in Fig. 2d–i. 3.2. Estimation and comparison of the 3D displacements By following the two methods (i.e., the linear inversion and InSAR– DGT methods) as presented in Section 2, we performed estimations to recover the 3D displacements using the simulated InSAR measurements (see Fig. 2g–i). Fig. 3a–c shows the initial E–W, N–S and U–D displacements, respectively, derived by the linear inversion method, while Fig. 3d–f shows the corresponding displacements derived by the InSAR–DGT method. As the InSAR–DGT method can conduct interpolation for the NPs, it can be observed that the InSAR–DGT results are cleaner and more continuous than the results derived by the linear inversion method. For further analysis, we calculated the residuals between the true displacements (see Fig. 2a–c) and the estimated displacements in the E–W, N–S and U–D dimensions, respectively. The maps of residuals obtained at the PPs or VPs are shown in Fig. 4. Fig. 4a–c shows the maps of residuals at the PPs for the case of the linear inversion method. Fig. 4d–f shows the maps of residuals at the PPs for the case of the InSAR–DGT method. Fig. 4g–i shows the maps of residuals at the VPs for the case of the InSAR–DGT method. A closer inspection to Fig. 4a–c versus Fig. 2d–f indicates that the residuals for the linear inversion case are similar in pattern to the spatially correlated noises. It is clear that the errors in the InSAR measurements are linearly propagated into the 3D displacements estimated by the linear inversion method. However, the InSAR–DGT residuals (see Fig. 4d–i) are smoothed due to the weighted LS solution. As the significant deformation gradients generally occur along the fault, some bigger residuals appear around the fault line for the InSAR–DGT results. The statistical calculation shows that the root mean square errors (RMSEs) at all the PPs for the linear inversion case are 1.52, 1.63 and 0.91 cm in the E–W, N–S and U–D dimensions, respectively, while the counterparts for the InSAR–DGT case are 0.56, 0.79 and 0.38 cm, respectively. This indicates that the accuracy in displacements derived by the InSAR–DGT method can be improved remarkably. In addition, the RMSEs at all the VPs are 0.12, 0.23 and 0.14 cm in the E–W, N–S and U–D dimensions, respectively, thus further demonstrating that the InSAR–DGT method can efficiently and accurately recover the 3D deformation field that is not smoothed significantly by the weighted LS solution. To check if the InSAR–DGT results are reasonable, we also calculated the displacements at the NPs (i.e., the black dots in Fig. 2d–i) by the Kriging-interpolation method (e.g., Kim et al., 2010; Yaseen, Hamm, Woldai, Tolpekin, & Stein, 2013). Fig. 5a–c shows the InSAR–DGT and Kriging-interpolation results at all the NPs versus the true displacements in the E–W, N–S and U–D dimensions, respectively. It can be observed that the Kriging-interpolation results divert much more from the true values than the InSAR–DGT results. The statistical calculation

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

393

Fig. 2. (a)–(c) show the maps of simulated coseismic displacements for the E–W, N–S and U–D dimensions, respectively, in an area of 40 × 40 km2 (corresponding to 200 × 200 pixels). The fault trace is marked by the white line. (d)–(f) show the maps of noises simulated for three cases of InSAR measurements, i.e., the descending LOS and AT displacements and the ascending AT displacements, respectively, while (g)–(i) show the maps of InSAR measurements with noises for the same three cases, respectively. Each of the null points (NPs) is indicated by a black dot in (d)–(i).

shows that the correlation coefficients between the InSAR–DGT results and the true displacements are 0.995, 0.997 and 0.991 in the E–W, N– S and U–D dimensions, respectively, while the counterparts between the Kriging-interpolation results and the true displacements are 0.967,

0.956 and 0.974, respectively. This indicates that the InSAR–DGT method is more advantageous than the Kriging-interpolation method for dealing with the NPs. This can be explained by that the former takes into account the spatial correlation in both the statistical and elastic

Fig. 3. (a)–(c) show the initial E–W, N–S and U–D displacements, respectively, derived by the linear inversion method, while (e)–(f) show the corresponding displacements derived by the InSAR–DGT method.

394

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

Fig. 4. The residuals between the true displacements (see Fig. 2a–c) and the estimated displacements in the E–W, N–S and U–D dimensions, respectively. (a)–(c) show the maps of residuals at the PPs for the case of the linear inversion method. (d)–(f) show the maps of residuals at the PPs for the case of the InSAR–DGT method. (g)–(i) show the maps of residuals at the VPs for the case of the InSAR–DGT method.

senses, while the latter considers the spatial correlation in only the statistical sense. 3.3. Estimation and comparison of the fault slips The coseismic 3D displacements can be used to estimate the fault slips based on the Okada inversion model (Okada, 1992). We performed the slip inversion using the three types of dataset (i.e., the true, linear inversion and InSAR–DGT 3D displacements), thus generating the three types of distribution of the fault slips. For comparison purpose, Fig. 6a

shows the map of distribution of fault slips (referred to as the true slips) derived from the true 3D displacements (see Fig. 2a–c), while Fig. 6b and c shows the maps of residuals between the true slips and the fault slips derived from the linear inversion and InSAR–DGT 3D displacements, respectively. The statistical calculation indicates that the RMSEs of the residuals for the linear inversion case and the InSAR– DGT case are 0.035 and 0.021 m, respectively, thus demonstrating that the accuracy in the fault slips derived from the InSAR–DGT displacements has been raised 40% if compared with that derived from the linear inversion displacements. Therefore, it can be concluded that the

Fig. 5. (a)–(c) show the InSAR–DGT and Kriging-interpolation displacements at all the NPs versus the true displacements in the E–W, N–S and U–D dimensions, respectively. The Kriginginterpolation results divert much more from the true values than the InSAR–DGT results.

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

395

1998), which runs about N70°E nearby the frontier between Myanmar and Laos. As shown in Fig. 7, the 2011 Tarlay earthquake was captured by the L-band PALSAR sensor onboard the Japanese ALOS satellite along the ascending and descending orbits shortly before and after the event. To map the 3D deformation field associated with the event, we will utilize the two ascending PALSAR images (track 486) acquired on February 16 and April 3 of 2011, respectively, and the two descending images (track 126) acquired on February 14 and April 1 of 2011, respectively, thus forming the ascending and descending interferometric pairs. The coverage for both the ascending and descending PALSAR images is indicated in Fig. 7. All the PALSAR images were collected at a radar off-nadir angle of 38.7° in HH polarization and are generated with a pixel spacing of 4.7 m in slant range and 3.2 m in azimuth (i.e., along-track). The heading angles of the ALOS flight direction along the ascending and descending orbits are 349.3° and 190.7°, respectively. The numerical analysis shows that the perpendicular baselines of the ascending and descending interferometric pairs are 37 and 438 m, respectively, while the two pairs have the same temporal baseline (i.e., time interval) of 46 days. The two interferometric pairs will be used as the basic datasets for estimating the coseismic deformation related to the Tarlay earthquake by following the procedures as described in Section 2. As a remark, the previous studies investigated the ground deformation pattern and the source mechanism of the Tarlay earthquake using the PALSAR data (Feng et al., 2013; Sun et al., 2013; Wang et al., 2014b). However, the 3D deformation field caused by the mainshock has not been documented yet. In Fig. 7, the study area is marked by the black rectangle for mapping the 3D displacements associated with the Tarlay earthquake. It can be seen that the topographic heights in the study area range between 174 and 2535 m. 4.2. The LOS and AT displacements caused by the Tarlay earthquake

Fig. 6. (a) shows the map of distribution of fault slips (referred to as the true slips) derived from the true 3D displacements, while (b) and (c) show the maps of residuals between the true slips and the fault slips derived from the linear inversion and InSAR–DGT 3D displacements, respectively.

InSAR–DGT method is obviously advantageous in terms of providing the high quality of 3D displacements for the coseismic modelling and analysis. 4. Mapping the coseismic deformation field of the 2011 Tarlay earthquake 4.1. The study area and data source On March 24 of 2011 (UTC Time 13:55:12), a large earthquake with a magnitude of Mw6.8 occurred in the Eastern Myanmar nearby the border between Thailand and Laos, causing 74 fatalities and injuring over 110 people (Sun, Shen, Burgmann, & Xu, 2013). The catalog of the Global Centroid Moment Tensor (GCMT) indicates that the epicenter (20.62°N, 100.02°E) of the earthquake is close to the Tarlay town, and this event was the largest earthquake for this region in the past 50 years (Feng et al., 2013). Fig. 7 shows the tectonic settings around the epicenter of the 2011 Tarlay earthquake. The previous studies demonstrated that a set of nearly parallel ENE left-lateral faults (marked by the black lines in Fig. 7) are distributed in this region due to both the tectonic effects of the India–Asia collision (i.e., north–south converging) and the eastward East Asia continental extrusion (Molnar & Tapponnier, 1975; Styron, Taylor, & Okoronkwo, 2010; Taylor & Yin, 2009; Wang, Lin, Simons, & Tun, 2014b). The Tarlay earthquake occurred in the west part of the Nam Ma Fault with a left-lateral geologic slip velocity of about 2.4 × 0.4 mm/yr (Lacassin, Replumaz, & Hervé,

We extracted the LOS and AT ground displacements caused by the Tarlay earthquake through the DInSAR and MAI processing, respectively, using the ascending and descending interferometric pairs as described in Section 4.1. Fig. 8 shows the resultant coseismic maps of LOS and AT displacements with spacing of 25 m by 25 m in the study area. The LOS displacements derived for the ascending and descending cases are depicted in Fig. 8a and b, respectively. It should be noted that the LOS displacements are re-wrapped with a color cycle of 11.8 cm. The AT displacements derived for the ascending and descending cases are depicted in Fig. 8c and d, respectively. The dashed black polygons marked in Fig. 8c and d corresponds to the areas deformed remarkably. It can be seen from Fig. 8a and b that both the ascending and descending interferometric pairs have captured the significant surface displacements along the radar LOS due to the mainshock. The ground moved in two reverse directions along a boundary line (corresponding to the west part of the Nam Ma Fault, marked by the dashed white line in Fig. 8) with a strike of about N70°E, which divides the fault structure into the northern wall and the southern wall. From the map of ascending LOS displacements (Fig. 8a), six fringes can be observed on both the northern and southern walls, reflecting that the relevant lands have moved about 70.8 cm toward (i.e., land uplifting along the radar LOS) and away from (i.e., land sinking along the radar LOS) the ALOS satellite, respectively. On the contrary, six fringes (about 70.8 cm of land sinking along the radar LOS) and seven fringes (about 82.6 cm of land uplifting along the radar LOS) can be observed on the northern and southern walls, respectively, from the map of descending LOS displacements (Fig. 8b). It can be seen from Fig. 8c and d that both the ascending and descending interferometric pairs have captured the surface displacements along the track due to the mainshock. As shown in Fig. 8c, the AT displacements derived from the ascending pair mainly appear on the northern wall, while the relatively smaller magnitude of displacements

396

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

Fig. 7. The tectonic settings around the epicenter of the 2011 Tarlay earthquake (the background is the shaded relief). The active and suspected fault traces (Wang et al., 2014b) are marked by the solid and dashed black lines, respectively. The location of the epicenter and the Tarlay town are indicated by the red star and the yellow square, respectively. The dashed blue lines represent the borders among the countries. The two dashed white rectangles mark the coverage of the ALOS PALSAR images acquired along the ascending and descending orbits, respectively. The black rectangle marks the study area for mapping the 3D coseismic deformation field.

can be observed on the southern wall. From the map of descending AT displacements (Fig. 8d), the obvious displacements can be observed on both the northern wall and the southern walls, which range between −82.6 and 91.2 cm. The most significant AT displacements with a maximum value of 91.2 cm appear in the two areas marked by the dashed black polygons in Fig. 8d nearby the predicted fault trace, and the displacement directions of the two areas are opposite. It should be pointed out that the ascending pair is less sensitive to the AT displacements caused by the mainshock, if compared with the descending pair. This is because the ascending flight direction of the ALOS satellite is almost perpendicular to the strike of the fault plane. In addition, the signal to noise ratio (SNR) in the LOS displacements derived by DInSAR is higher than that in the AT displacements derived by MAI. This is because the DInSAR interferograms are generated using the full-aperture SLC images with higher SNR, while the MAI interferograms are generated by the sub-aperture SLC images with lower SNR.

ascending DInSAR results (Fig. 8a) and both the descending DInSAR (Fig. 8b) and MAI (Fig. 8d) results were used to generate the maps of the initial 3D coseismic displacements. We further improved the 3D displacements in the study area on a point-by-point basis by using the InSAR–DGT method as described in Section 2.2. To select the VPs for an arbitrary point to be analyzed, we utilized the three thresholds of correlation distance, interferometric coherence and displacement dispersion (see Section 2.2). Fig. 9d–f shows the masked maps of the initial 3D displacements using the interferometric coherence threshold of 0.8. However, the correlation distance threshold should be determined through numerical analysis. Based on the initial far-field 3D displacements and Eq. (6), we determined the empirical covariance functions for the 3D components as follows: 8 e < C ðr Þ ¼ 1:02  expðr=5:26Þ C n ðr Þ ¼ 8:37  expðr=1:54Þ : : u C ðr Þ ¼ 3:11  expðr=2:13Þ

ð17Þ

4.3. The 3D deformation maps of the Tarlay earthquake By following the linear inversion method as described in Section 2.1, we estimated the initial 3D coseismic displacements through the DInSAR and MAI processing. Fig. 9a–c shows the maps of the initial 3D coseismic displacements in the study area caused by the Tarlay earthquake. As mentioned before, the ascending pair is not sensitive to the AT displacements, and the MAI results from this pair were not used for estimating the initial 3D displacements. It means that only the

From Eq. (17), one can see that the spatial correlation distance for the E–W, N–S and U–D dimensions are 5.26, 1.54 and 2.13 km, respectively. The average correlation distance of the three dimensions is 2.98 km, taking as the correlation distance threshold for the study area. For an arbitrary point, once the true VPs are determined, we computed the 3D displacements for the point by the LS adjustment. Fig. 9g–i shows the maps of the resultant 3D displacements improved by the InSAR–DGT method.

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

397

Fig. 8. The coseismic displacement maps of the Tarlay earthquakes. (a) and (b) are for the re-wrapped LOS displacements derived from the ascending and descending pairs, respectively, in which one fringe represents the LOS displacements of 11.8 cm. (c) and (d) are for the AT displacements derived from the ascending and descending pairs, respectively, in which the dashed polygons indicate the areas deformed remarkably. The dashed white line denotes the predicted fault trace associated with the Tarlay earthquake.

It can be seen from Fig. 9g–i that all the displacements in each dimension caused by the Tarlay earthquake are anti-symmetric with respect to the fault trace and similar to the prediction by a linear elastic deformation model related to a strike-slip fault (Fialko et al., 2001). The displacements in the E–W, N–S and U–D dimensions range between −116.3 and 124.5 cm, −56.2 and 60.3 cm, −23.6 and 26.2 cm, respectively. The strongest displacements occurred in the E–W dimension, while the relatively smaller displacements occurred in the N–S and U–D dimensions. The vertical displacements mainly occurred at both the ends of the fault trace, and this is very consistent with the field investigation by Tun et al. (2014). Moreover, the magnitude of vertical displacements on the southern wall is greater than that on the northern wall, indicating that the fault plane most likely dipped into the southeast and generated a minor normal slip component during the earthquake. The pattern of the resultant 3D deformation field reveals that the causative fault for the Tarlay earthquake had a left-lateral slip that was accompanied by a minor normal dip-slip component, and this coincides with the seismic inversion results documented by the USGS and GCMT (see Acknowledgments). The comparison between the maps of the initial (Fig. 9a–c) and improved (Fig. 9g–i) 3D displacements shows that the number of the missing data points near the fault trace are reduced apparently by the InSAR–DGT method. To check the precision in the improved 3D displacements, we randomly selected three far-field areas (each area with a size of 50 × 50 points) marked by the white rectangles Z1–Z3 in Fig. 9a–c and Fig. 9g–i, in which zero deformation can be assumed. It should be emphasized that the selected areas are far away from the fault trace and not within the significantly deformed zone marked by the dashed black ellipse in Fig. 9, thus guaranteeing the robust error analysis. The mean values and the standard deviations (SDs) of the E–W, N–S and U–D displacements in each of the selected areas were calculated for the two types of displacements derived by the linear inversion method and the InSAR–DGT method, respectively. Table 1 lists the comparison of the mean values and SDs calculated for the two types of displacements in each of the selected areas. Moreover, the results derived with all the far-field points show that the SDs of the E–W, N–S and U–D displacements derived by the linear inversion

method are 2.21, 8.83 and 1.78 cm, respectively, while those by the InSAR–DGT method are 1.73, 5.67 and 1.36 cm, respectively, thus indicating a precision improvement of 22%, 36% and 24%, respectively. Therefore, we can conclude that the InSAR–DGT method can significantly improve the precision in the final 3D displacements. 5. Validation and interpretation of the coseismic displacements 5.1. Relation between the number of VPs and the precision of 3D displacements As described in Section 2.2, selecting the VPs is the crucial step that has impact on the precision of the 3D deformation field derived by the InSAR–DGT method. Here we analyzed the relationship between the number of VPs and the SDs of the 3D displacements estimated by InSAR–DGT. It should be noted that the SDs can be calculated using Eq. (16) during the weighted LS adjustment. Fig. 10a shows the number of VPs obtained for each point in the study area. Fig. 10b shows the histogram of the VPs, representing the variation of frequencies of the points in the study area associated with the different number of VPs. Fig. 10c shows the variation of SDs in the E–W, N–S and U–D displacements, respectively, as a function of the number of VPs. The parts in red color in Fig. 10a corresponds to the large number of VPs at the estimated points, and the maximum number of VPs reaches up to 1845 in the study area. Comparing the initial (see Fig. 9a–c) and improved (see Fig. 9g–i) 3D displacement maps, it can be found that the noises in the parts of points with the large number of VPs have been reduced significantly by the InSAR–iDGT method. As shown in Fig. 10c, the SDs for each dimension decrease with an increase in the number of VPs, and such drop trend becomes insignificant when the number of VPs reaches to a certain value (e.g., 200 VPs for this study). It can be observed from Fig. 10c that the maximum SDs in the E–W, N–S and U–D dimensions are 1.62, 3.87 and 1.35 cm, respectively. If the number of VPs is more than 80, the SDs in the E–W and U–D displacements are less than 1 cm, and the cumulative frequency of SDs of less than 1 cm is 75.3% (see Fig. 10b and c). This means that 75.3% points in the study area have the SDs of less than 1 cm in the E–W and U–D

398

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

Fig. 9. The 3D coseismic deformation maps for the Tarlay earthquakes. (a)–(c) are for the initial 3D displacements derived by the linear inversion method. (d)–(f) are for the initial 3D displacements masked using the interferometric coherence threshold of 0.8. (g)–(i) are for the resultant 3D displacements improved by the InSAR–DGT method. For better comparison, the displacements in each dimension are scaled into the same color bar. The dashed black ellipse denotes the significantly deformed zone around the fault trace. Three far-field areas marked by the white rectangles Z1–Z3 are selected for precision comparison.

displacements. If the number of VPs is more than 166, the SDs in the N–S displacements are less than 2 cm, and the relevant cumulative frequency is 66.7%. However, it should be noted that the quality of displacements at a point estimated by the InSAR–DGT method is not only affected by the number of VPs, but also the reliability of the initial displacements at the VPs. In addition, it should be emphasized that the SDs in the E–W and U–D displacements are very close to each other, which are obviously smaller than the SDs in the N–S displacements. This can be ascribed to the E–W and U–D displacements are mainly constrained by the DInSAR measurements that are usually more precise than the MAI measurements (Bechor & Zebker, 2006; Wang et al., 2014a).

5.2. Validation on the fault slips with the existing studies In order to confirm the source of the Tarlay earthquake, a survey team conducted a field investigation about two weeks after the earthquake (Tun et al., 2014). The fieldwork showed that the fractures caused by the earthquake were widely distributed within a 2-km-width zone around the main fault trace with a strike of about N70°E. The interpretation with aid of the high spatial resolution (0.5 m) Worldview-2 images suggested that the ground rupture length and the maximum offset are about 30 km and 2 m, respectively (Tun et al., 2014).

For comparison between the results derived by SAR interferometry and Tun et al. (2014), we projected the E–W and N–S displacements (i.e., horizontal displacement) at each point onto the fault orientation (N70°E) to obtain the fault slip map of the Tarlay earthquake. Fig. 11a and b shows the fault slip maps in the study area calculated from the initial and improved horizontal displacements derived by the linear inversion method and the InSAR–DGT method, respectively. It can be seen that the fault slip map derived by InSAR–DGT is cleaner than that derived by the linear inversion method, although the fault trace can be easily identified from both of them. The slip offsets mainly distributed within a 10-km-width zone around the fault trace. The rupture length of about 29.18 km can be estimated from Fig. 11b, which is very consistent with the report by Tun et al. (2014). For further comparison of displacements, we selected two profiles marked in Fig. 11 by the dashed black lines AA′ and BB′ which are across the fault trace. The two points P1 and P2 (close to the fault trace) with the known slip offsets measured by the field survey (Tun et al., 2014) are included in the two reference lines selected. Fig. 12a and b shows the comparison between the profiles of slip offsets extracted from Fig. 11a (by linear inversion) and Fig. 11b (by InSAR–DGT) along AA′ and BB′, respectively. The shaded parts denote the SDs of slip offsets, and the shaded gray zone marks the fault location. It can be seen that the slip offsets derived from the improved displacements are smoother and more consistent near the fault trace than those derived from the

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404 Table 1 The comparison of the mean values and SDs calculated for the two types of 3D displacements in the far-field areas. Statistical item

Mean values (cm)

SDs (cm)

Dimension

E–W

E–W

For the area Z1 For the initial displacements For the improved displacements Improved percentage For the area Z2 For the initial displacements For the improved displacements Improved percentage For the area Z3 For the initial displacements For the improved displacements Improved percentage For all the far-field points For the initial displacements For the improved displacements Improved percentage

N–S

U–D

N–S

U–D

−2.67

5.45

0.78

1.56

9.81

1.67

−2.14

3.63

0.59

1.42

5.77

1.38

20%

33%

24%

9%

41%

17%

−2.99

−6.86

−1.31

1.02

7.25

1.03

−2.24

−2.59

−1.12

0.92

2.08

0.48

25%

61%

15%

10%

71%

53%

3.71

−4.66

−0.82

1.09

8.49

2.52

−3.21

−0.64

0.88

3.43

32%

31%

22%

−1.97

−2.86

−0.61

2.21

8.83

1.78

−1.48

1.67

−0.49

1.73

5.67

1.36

25%

42%

20%

19%

22%

60%

36%

1.87 1.34 28%

24%

initial displacements. We calculated the sinistral slips at P1 and P2 by differencing the mean slip offsets between the northern and southern walls. Based on the slip offsets derived by InSAR–DGT, it can be inferred that the sinistral slips at P1 and P2 are 76.8 ± 2.2 and 61.8 ± 3.9 cm, respectively. However, the field survey by Tun et al. (2014) indicates that the sinistral slips at P1 and P2 are 71.1 and 46.5 cm, respectively. It is clear that the two types of sinistral slips at P1 are generally consistent, and the mutual discrepancy is about 5.7 cm, while the discrepancy between the two types of sinistral slips at P2 is about 15.3 cm. The large discrepancy at P2 is most likely because P2 is located in the transtensional basin (i.e., northeast the Tarlay town) filled with lower brittle strength of the saturated fluvial sediments, where the liquefaction due to the mainshock might occur (Wang et al., 2014b). Moreover, the widely distributed off-fault displacements transecting the fault trace may be another factor for the larger sinistral slip at P2 by InSAR–DGT. It is also worth noting that the slip offsets measured by field survey might strongly be affected by the local lithology, topography and geomorphic features, probably resulting in biased analysis in many geophysical applications (Dolan & Haravitch, 2014). Nevertheless, the numerical comparison at P1 and P2 indicates that the slip offsets measured by

399

InSAR–DGT are comparable with the field surveys. The high density of displacements generated by InSAR–DGT is valuable for estimating the fault slip offsets without field trip. 5.3. Validation on the 3D deformation field with the existing optimized fault model The further validation was performed by comparing the results of 3D displacements derived by the InSAR–DGT method with those simulated from a fault model by Feng et al. (2013) for the Tarlay earthquake using the Okada elastic half-space dislocation theory (Okada, 1992). Based on DInSAR and MAI data processing, Feng et al. (2013) reported that one fault segment could explain the interferometric measurements of the Tarlay earthquake. Table 2 lists the optimized source parameters of the one-segment fault model derived by Feng et al. (2013) for the Tarlay earthquake. Fig. 13a–c shows the maps of 3D displacements (in the E–W, N–S and U–D dimensions, respectively) in the study area simulated using the source parameters as listed in Table 2. Fig. 13d–f shows the maps of residuals between the 3D displacements derived by simulation and the linear inversion method, while Fig. 13g–i shows the maps of residuals between the 3D displacements derived by simulation and by the InSAR–DGT method. It can be seen that the simulated 3D deformation maps (Fig. 13a–c) are in good agreement with those (see Fig. 9a–c and Fig. 9g–i) derived by the linear inversion and InSAR–DGT method. However, the comparison analysis shows that the improved 3D displacements have the smaller magnitude of residuals than the initial 3D displacements. Fig. 13j–l shows the histograms of residuals between the displacements derived by simulation and the InSAR–DGT method. The statistical calculation shows that the mean values of the discrepancies between the two types of displacements are 0.39, −1.76, and −0.16 cm in the E–W, N–S and U–D dimensions, respectively, and the root mean square errors (RMSEs) are 4.48, 5.78 and 2.61 cm in the E–W, N–S and U–D dimensions, respectively. The discrepancies between the two types of displacements in the 95% confidence interval are [−8.57, 9.35], [−13.32, 9.80] and [− 5.38, 5.06] cm in the E–W, N–S and U–D dimensions, respectively. Several factors are responsible for the discrepancies between the two types of displacements derived by simulation and InSAR–DGT, respectively. First, the parameters (e.g., length and strike) related to the fault geometry of the Tarlay earthquake was mapped using the AT displacements by Feng et al. (2013), which is slightly different from those mapped from the fault slip map derived by this study. Second, the causative fault system associated with the Tarlay earthquake may be more complicated than the one-segment fault taken by Feng et al. (2013), for example, the fault system may changes its orientation in

Fig. 10. (a) shows the number of VPs for each point in the study area using the form of counter lines. (b) shows the histogram of the VPs. (c) shows the variation of SDs in the E–W, N–S and U–D displacements, respectively, as a function of the VP numbers.

400

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

Fig. 11. The fault slip maps in the study area. (a) is obtained from the horizontal displacements derived by the linear inversion method. (b) is obtained from the horizontal displacements derived by the InSAR–DGT method. The white line denotes the fault trace and the white arrows denote the slip directions. For further analysis, the dashed black lines AA′ and BB′ are marked across the two points P1 and P2 with the known slip offsets measured by the field survey (Tun et al., 2014). The black rectangle marks the area that will be further analyzed in Section 5.4.

the eastern part of the Tarlay town (see Section 5.4). Finally, the errors (e.g., atmospheric effects, orbital errors and topographic errors) in the DInSAR and MAI measurements may be propagated into the 3D displacements derived by the InSAR–DGT method (see Section 2.3). 5.4. Interpretation of the causative fault system for the Tarlay earthquake The analysis in Section 5.3 indicates that the one-segment fault model can basically explain the 3D deformation field derived by our InSAR–DGT method. However, the discrepancies between the displacements in the E–W or N–S directions derived by simulation and the InSAR–DGT method are over about 10 cm (see Fig. 13j–l). In addition, the distribution of slips on the one-segment fault model inverted by Sun et al. (2013) shows that the fault patches on the east edge of the fault plane have slips as large as 2.0 m, even though the smoothing constraints were applied toward the fault plane edges (Sun et al., 2013).

Fig. 12. The comparison between the profiles of slip offsets extracted from the fault slip maps as shown in Fig. 11a (by linear inversion) and Fig. 11b (by InSAR–DGT) along AA′ and BB′ respectively. The shaded parts denote the SDs of slip offsets, and the shaded gray zone marks the fault location.

This implies that the fault plane probably extend further to east. With careful interpretation of the PALSAR interferograms, Sun et al. (2013) argued that a short second fault segment (with a strike of N58°E) connecting the east edge of the main segment (with a strike of N70°E) should be considered, which can be used to accommodate the large discrepancies as shown in Fig. 13j–l. In order to investigate whether the causative fault system for the Tarlay earthquake extends further to east and changes its orientation, a comprehensive interpretation is presented below by taking into account the fault slip map, the field survey and the tectonic settings. According to the fault slip map (Fig. 11b), the field survey by Tun et al. (2014) and the tectonic settings (Fig. 7), Fig. 14 exhibit the detailed fault slips and surface ruptures in the east part (marked by the black rectangle in Fig. 11b) of the main fault trace. Fig. 14a shows the enlarged fault slip map where two reference lines CD and CE and a feature point T are marked for further analysis. Fig. 14b and c shows the variation of slips along CD and CE, respectively. Fig. 14d shows the map of possible fault structures related to the Tarlay earthquake, where the red arrows denote the slip directions of the fault segments. Fig. 14e and f shows the photos taken around the point T by Tun et al. (2014), indicating the localized surface ruptures. One can observe from Fig. 14a that two separate linear features clearly form the branched pattern at about 5 km east of the Tarlay town. The northern branch strikes N40°E and extends further northeast, while the southern branch is almost along the E–W direction. It can be seen from Fig. 14b that the magnitude of fault slips gradually increases along CD, and there is an abrupt jump at the location of about 1 km away from point C. The slips along CE show a similar trend, and the abrupt jump appears at the location of about 1.5 km away from point C. Such features imply that the two linear branches may be of the extension of the main fault trace, and the two fault segments might also generate the left-lateral slip during the rupture process. Based on the interpretation, we propose a sketch (as shown in Fig. 14d) of the plausible fault structures related to the Tarlay earthquake. The field investigation by Tun et al. (2014) indicated that a group of ground failures (see Fig. 14e and f) occurred in the Tarlay basin with a strike of about 20–40° that is obviously different from the orientation (70°) of the surface rupture along the main fault trace. Tun et al. (2014) speculated that such discrepancy may be ascribed to the localized ground liquefaction in the transtensional Tarlay basin. A closer inspection with Fig. 14d indicates that a combined effect of the E–W extension and the N–S compression might take place around the northeast part of the Tarlay town, which may lead to the liquefaction of the saturated fluvial sediments that fill the basin. Due to the slip mechanism, a releasing band with a width of about 1.5 km appears on the north edge of the basin. We therefore conclude that the Tarlay

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

401

Table 2 The optimized source parameters of the one-segment fault model derived by Feng et al. (2013) for the Tarlay earthquake. Source parameters

Longitude (°)

Latitude (°)

Length (km)

Width (km)

Depth (km)

Strike (°)

Dip (°)

Rake (°)

Mw (−)

Values

99.995

20.674

22.4

8.4

4.6

69.7

87.3

1.8

6.8

earthquake occurred in the west part of the Nam Ma Fault and infer that the main fault most likely extends into the Tarlay basin and separates into two branches at about 5 km east of the Tarlay town. It should be pointed out that such speculation requires further investigation and verification, thus leaving a room for future work.

6. Conclusions To map the 3D deformation field related to a geophysical event, this paper proposes an integrated method termed InSAR–DGT by combining DInSAR for detecting LOS displacements, MAI for detecting AT

Fig. 13. (a)–(c) are for the maps of 3D displacements (in the E–W, N–S and U–D dimensions, respectively) in the study area simulated using the source parameters as listed in Table 2. (d)–(f) are for the maps of residuals between the 3D displacements derived by simulation and the linear inversion method. (g)–(i) are for the maps of residuals between the 3D displacements derived by simulation and by the InSAR–DGT method. (j)–(l) are for the histograms of residuals between the 3D displacements derived by simulation and the InSAR–DGT method.

402

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

Fig. 14. The detailed fault slips and surface ruptures in the east part (marked by the black rectangle in Fig. 11b) of the main fault trace. (a) is the enlarged fault slip map where two reference lines CD and CE and a feature point T are marked for further analysis. (b) and (c) are the variations of slips along CD and CE, respectively. (d) is the map of possible fault structures related to the Tarlay earthquake, where the red arrows denote the slip directions of the fault segments. (e) and (f) are the photos taken around the point T by Tun et al. (2014), indicating the localized surface ruptures.

displacements, and the DGT model for characterizing spatial correlation of ground displacements. The InSAR–DGT method takes into account the spatial correlation in both the statistical and elastic senses and establishes relationship between the high and low quality of InSAR observation points by incorporating the DGT model, thus improving estimation of the 3D displacements with a weighted LS solution. The proposed method was first tested through an experiment of recovering 3D displacements from a simulated coseismic deformation field. We then applied the proposed method to map the 3D coseismic deformation field related to the 2011 Tarlay earthquake (Myanmar) by using the ascending and descending ALOS PALSAR images. The testing results with the simulated experiment show that the InSAR–DGT method can efficiently and accurately recover the 3D deformation field that is not smoothed significantly by the weighted LS solution. The accuracy in the fault slips derived from the InSAR–DGT displacements can be raised 40% if compared with that from the linear inversion displacements. The testing results with the Tarlay earthquake show that the quality of the 3D coseismic displacements derived by InSAR–DGT can be raised efficiently, and the precisions in the E–W, N–S and U–D displacements are increased 22%, 36% and 24%, respectively. Moreover, the number of the missing data points (particularly near the fault trace) in the 3D coseismic deformation field can be reduced significantly by the InSAR–DGT method. The 3D coseismic deformation field improved by InSAR–DGT indicates that the peak-to-peak coseismic displacements caused by the

Tarlay earthquake are about 240.8, 116.5 and 49.8 cm in the E–W, N–S and U–D dimensions, respectively. It can be easily revealed from the 3D displacement maps that the causative fault (i.e., the west part of the Nam Ma Fault) for the Tarlay earthquake generated a left-lateral slip that was accompanied by a minor normal dip-slip component. In addition, it can be found from the fault slip map generated using the 3D displacements that the ground rupture length is about 29.18 km, which is very consistent with the result (about 30 km) derived by the field survey. The validation also shows that the improved deformation field is in good agreement with the deformation field simulated from the existing optimized fault model. The comprehensive interpretation demonstrates that the main fault most likely extends into the Tarlay basin and separates into two branches at about 5 km east of the Tarlay town. Such speculation requires further investigation and verification, thus leaving a room for future work. As a final remark, it should be emphasized that the InSAR–DGT method estimates the 3D displacements using the InSAR measurements at the adjacent points on the assumption of spatial correlation in the statistical and elastic senses. If such hypothesis cannot be completely satisfied for a study area with the special geological settings, the performance of the InSAR–DGT method might be degraded for reconstruction of 3D displacements. However, for the scenarios where the field trip is hard to carry out and only the satellite SAR images are available, the InSAR–DGT method is an effective approach for mapping the highresolution 3D deformation field. It can be applied to monitor and

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

analyze the complicated ground displacements associated with earthquakes and volcanoes, thus providing valuable data source for modeling the geographical events. Acknowledgments This work was jointly supported by the National Basic Research Program of China (973 Program) under Grant 2012CB719901, the National Natural Science Foundation of China under Grant 41474003, and the Program for Changjiang Scholars and Innovative Research Team in University under Grant IRT13092. We thank JAXA for providing the PALSAR data. We also thank JPL/Caltech and USGS for providing the software ROI_PAC and SRTM DEM, respectively, and thank the two anonymous reviewers for their valuable comments. Some parameters related to the Tarlay earthquake were obtained from the GCMT and USGS databases (available at www.globalcmt.org and http://earthquake.usgs. gov/earthquakes, respectively). Some figures in this paper were plotted using the Generic Mapping Tools (Wessel & Smith, 1998). References Barbot, S., Hamiel, Y., & Fialko, Y. (2008). Space geodetic investigation of the coseismic and postseismic deformation due to the 2003 Mw7.2 Altai earthquake: Implications for the local lithospheric rheology. Journal of Geophysical Research, 113, B03403. http://dx.doi.org/10.1029/2007jb005063. Barnhart, W. D., Briggs, R. W., Reitman, N. G., Gold, R. D., & Hayes, G. P. (2015). Evidence for slip partitioning and bimodal slip behavior on a single fault: Surface slip characteristics of the 2013 Mw7.7 Balochistan, Pakistan earthquake. Earth and Planetary Science Letters, 420, 1–11. http://dx.doi.org/10.1016/j.epsl.2015.03.027. Bechor, N. B. D., & Zebker, H. A. (2006). Measuring two-dimensional movements using a single InSAR pair. Geophysical Research Letters, 33http://dx.doi.org/10.1029/ 2006gl026883. Cardozo, N., & Allmendinger, R. W. (2009). SSPX: A program to compute strain from displacement/velocity data. Computers & Geosciences, 35, 1343–1357. http://dx.doi.org/ 10.1016/j.cageo.2008.05.008. Dolan, J. F., & Haravitch, B. D. (2014). How well do surface slip measurements track slip at depth in large strike-slip earthquakes? The importance of fault structural maturity in controlling on-fault slip versus off-fault surface deformation. Earth and Planetary Science Letters, 388, 38–47. http://dx.doi.org/10.1016/j.epsl.2013.11.043. Erten, E., Reigber, A., & Hellwich, O. (2010). Generation of three-dimensional deformation maps from InSAR data using spectral diversity techniques. ISPRS Journal of Photogrammetry and Remote Sensing, 65, 388–394. http://dx.doi.org/10.1016/ j.isprsjprs.2010.04.005. Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., ... Alsdorf, D. (2007). The shuttle radar topography mission. Reviews of Geophysics, 45, RG2004http://dx.doi. org/10.1029/2005RG000183. Feng, W. P., Li, Z. H., Elliott, J. R., Fukushima, Y., Hoey, T., Singleton, A., ... Xu, Z. (2013). The 2011 Mw6.8 Burma earthquake: Fault constraints provided by multiple SAR techniques. Geophysical Journal International, 195, 650–660. http://dx.doi.org/10. 1093/gji/ggt254. Fialko, Y., Simons, M., & Agnew, D. (2001). The complete (3-D) surface displacement field in the epicentral area of the 1999 Mw7.1 Hector Mine earthquake, California, from space geodetic observations. Geophysical Research Letters, 28, 3063–3066. http://dx. doi.org/10.1029/2001GL013174. Funning, G. J., Parsons, B., Wright, T. J., Jackson, J. A., & Fielding, E. J. (2005). Surface displacements and source parameters of the 2003 Bam (Iran) earthquake from Envisat advanced synthetic aperture radar imagery. Journal of Geophysical Research — Solid Earth, 110http://dx.doi.org/10.1029/2004jb003338. González, P., Chini, M., Stramondo, S., & Fernandez, J. (2010). Coseismic horizontal offsets and fault-trace mapping using phase correlation of IRS satellite images: The 1999 Izmit (Turkey) earthquake. IEEE Transactions on Geoscience and Remote Sensing, 48, 2242–2250. http://dx.doi.org/10.1109/tgrs.2009.2039352. González, P., Fernández, J., & Camacho, A. (2009). Coseismic three-dimensional displacements determined using SAR data: Theory and an application test. Pure and Applied Geophysics, 166, 1403–1424. http://dx.doi.org/10.1007/s00024-009-0500-7. Gourmelen, N., Kim, S. W., Shepherd, A., Park, J. W., Sundal, A. V., Björnsson, H., & Pálsson, F. (2011). Ice velocity determined using conventional and multiple-aperture InSAR. Earth and Planetary Science Letters, 307, 156–160. http://dx.doi.org/10.1016/j.epsl. 2011.04.026. Gudmundsson, S., Sigmundsson, F., & Carstensen, J. M. (2002). Three-dimensional surface displacement maps estimated from combined interferometric synthetic aperture radar and GPS data. Journal of Geophysical Research — Solid Earth, 107, 2250. http://dx.doi.org/10.1029/2001JB000283. Guglielmino, F., Bignami, C., Bonforte, A., Briole, P., Obrizzo, F., Puglisi, G., ... Wegmüller, U. (2011a). Analysis of satellite and in situ ground deformation data integrated by the SISTEM approach: The April 3, 2010 earthquake along the Pernicana fault (Mt. Etna — Italy) case study. Earth and Planetary Science Letters, 312, 327–336. http://dx.doi.org/10.1016/j.epsl.2011.10.028.

403

Guglielmino, F., Nunnari, G., Puglisi, G., & Spata, A. (2011b). Simultaneous and integrated strain tensor estimation from geodetic and satellite deformation measurements to obtain three-dimensional displacement maps. IEEE Transactions on Geoscience and Remote Sensing, 49, 1815–1826. http://dx.doi.org/10.1109/tgrs.2010.2103078. Hanssen, R. F. (2001). Radar interferometry: Data interpretation and error analysis. Springer. Hu, J., Li, Z. W., Ding, X. L., Zhu, J. J., Zhang, L., & Sun, Q. (2014). Resolving threedimensional surface displacements from InSAR measurements: A review. EarthScience Reviews, 133, 1–17. http://dx.doi.org/10.1016/j.earscirev.2014.02.005. Jolivet, R., Lasserre, C., Doin, M. P., Guillaso, S., Peltzer, G., Dailu, R., ... Xu, X. (2012). Shallow creep on the Haiyuan Fault (Gansu, China) revealed by SAR interferometry. Journal of Geophysical Research — Solid Earth, 117, B06401http://dx.doi.org/10.1029/ 2011JB008732. Jung, H. S., Lee, W. J., & Zhang, L. (2014). Theoretical accuracy of along-track displacement measurements from multiple-aperture interferometry (MAI). Sensors (Basel), 14, 17703–17724. http://dx.doi.org/10.3390/s140917703. Jung, H. S., Lu, Z., Won, J. S., Poland, M. P., & Miklius, A. (2011). Mapping three-dimensional surface deformation by combining multiple-aperture interferometry and conventional interferometry: Application to the June 2007 eruption of Kilauea Volcano, Hawaii. IEEE Geoscience and Remote Sensing Letters, 8, 34–38. http://dx.doi.org/10. 1109/Lgrs.2010.2051793. Jung, H. S., Won, J. S., & Kim, S. W. (2009). An improvement of the performance of multiple-aperture SAR interferometry (MAI). IEEE Transactions on Geoscience and Remote Sensing, 47, 2859–2869. http://dx.doi.org/10.1109/Tgrs.2009.2016554. Kim, S. W., Wdowinski, S., Dixon, T. H., Amelung, F., Kim, J. W., & Won, J. S. (2010). Measurements and predictions of subsidence induced by soil consolidation using persistent scatterer InSAR and a hyperbolic model. Geophysical Research Letters, 37, L05304http://dx.doi.org/10.1029/2009gl041644. Lacassin, R., Replumaz, A., & Hervé, L. P. (1998). Hairpin river loops and slip-sense inversion on Southeast Asian strike-slip faults. Geology, 26, 703–706. http://dx.doi.org/10. 1130/0091-7613(1998)026b0703:HRLASSN2.3.CO;2. Leprince, S., Barbot, S., Ayoub, F., & Avouac, J. P. (2007). Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements. IEEE Transactions on Geoscience and Remote Sensing, 45, 1529–1558. http://dx.doi.org/10.1109/tgrs.2006.888937. Liu, B., Zhang, J., Luo, Y., Jiang, W., Chen, X., & Li, Y. (2014). Error propagation analysis in three-dimensional coseismic displacement inversion. IEEE Geoscience and Remote Sensing Letters, 11, 1971–1975. http://dx.doi.org/10.1109/LGRS.2014.2315815. Liu, G. X., Li, J., Xu, Z., Wu, J. C., Chen, Q., Zhang, H. X., ... Luo, X. J. (2010). Surface deformation associated with the 2008 Ms 8.0 Wenchuan earthquake from ALOS L-band SAR interferometry. International Journal of Applied Earth Observation and Geoinformation, 12(6), 496–505. http://dx.doi.org/10.1016/j.jag.2010.05.005. Massonnet, D., & Feigl, K. L. (1998). Radar interferometry and its application to changes in the earth's surface. Reviews of Geophysics, 36, 441–500. http://dx.doi.org/10.1029/ 97rg03139. Michel, R., & Avouac, J. P. (1999). Measuring near field coseismic displacements from SAR images: Application to the Landers earthquake. Geophysical Research Letters, 26(19), 3017–3020. http://dx.doi.org/10.1029/1999GL900524. Mohammed, O. I., Saeidi, V., Pradhan, B., & Yusuf, Y. A. (2013). Advanced differential interferometry synthetic aperture radar techniques for deformation monitoring: A review on sensors and recent research development. Geocarto International, 1–18. http://dx. doi.org/10.1080/10106049.2013.807305. Molnar, P., & Tapponnier, P. (1975). Cenozoic tectonics of Asia: Effects of a continental collision: Features of recent continental tectonics in Asia can be interpreted as results of the India–Eurasia collision. Science, 189, 419–426. http://dx.doi.org/10.1126/science.189.4201.419. Muller, C., del Potro, R., Biggs, J., Gottsmann, J., Ebmeier, S. K., Guillaume, S., ... Van der Laat, R. (2014). Integrated velocity field from ground and satellite geodetic techniques: Application to Arenal volcano. Geophysical Journal International, 200, 863–879. http://dx.doi.org/10.1093/gji/ggu444. Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 82, 1018–1040. Pietrantonio, G., & Riguzzi, F. (2004). Three-dimensional strain tensor estimation by GPS observations: Methodological aspects and geophysical applications. Journal of Geodynamics, 38, 1–18. http://dx.doi.org/10.1016/j.jog.2004.02.021. Prati, C., Ferretti, A., & Perissin, D. (2010). Recent advances on surface ground deformation measurement by means of repeated space-borne SAR observations. Journal of Geodynamics, 49, 161–170. http://dx.doi.org/10.1016/j.jog.2009.10.011. Raucoules, D., de Michele, M., Malet, J. P., & Ulrich, P. (2013). Time-variable 3D ground displacements from high-resolution synthetic aperture radar (SAR). Application to La Valette landslide (South French Alps). Remote Sensing of Environment, 139, 198–204. http://dx.doi.org/10.1016/j.rse.2013.08.006. Rosen, P. A., Hensley, S., Peltzer, G., & Simons, M. (2004). Updated repeat orbit interferometry package released. Eos, Transactions American Geophysical Union, 85, 47-47. http://dx.doi.org/10.1029/2004EO050004. Rucci, A., Ferretti, A., Monti, G. A., & Rocca, F. (2012). Sentinel 1 SAR interferometry applications: The outlook for sub millimeter measurements. Remote Sensing of Environment, 120, 156–163. http://dx.doi.org/10.1016/j.rse.2011.09.030. Samsonov, S., & Tiampo, K. (2006). Analytical optimization of a DInSAR and GPS dataset for derivation of three-dimensional surface displacement. IEEE Geoscience and Remote Sensing Letters, 3, 107–111. http://dx.doi.org/10.1109/ LGRS.2005.858483. Sansosti, E., Berardino, P., Bonano, M., Calo, F., Castaldo, R., Casu, F., ... Lanari, R. (2014). How second generation SAR systems are impacting the analysis of ground deformation. International Journal of Applied Earth Observation and Geoinformation, 28, 1–11. http://dx.doi.org/10.1016/j.jag.2013.10.007.

404

X. Wang et al. / Remote Sensing of Environment 170 (2015) 388–404

Segall, P. (2010). Earthquake and volcano deformation. Princeton University Press. Styron, R., Taylor, M., & Okoronkwo, K. (2010). Database of active structures from the Indo-Asian collision. Eos, Transactions American Geophysical Union, 91, 181–182. http://dx.doi.org/10.1029/2010EO200001. Sudhaus, H., & Jonsson, S. (2009). Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: Application to the June 2000 Kleifarvatn earthquake, Iceland. Geophysical Journal International, 176, 389–404. http://dx.doi.org/10.1111/j.1365-246X.2008.03989.x. Sun, J., Shen, Z. K., Burgmann, R., & Xu, X. (2013). Coseismic slip distribution of the 24 March 2011 Tarlay (Myanmar) Mw6.8 earthquake from ALOS PALSAR interferometry. Bulletin of the Seismological Society of America, 103, 2928–2936. http://dx.doi. org/10.1785/0120120365. Taylor, M., & Yin, A. (2009). Active structures of the Himalayan–Tibetan orogen and their relationships to earthquake distribution, contemporary strain field, and Cenozoic volcanism. Geosphere, 5, 199–214. http://dx.doi.org/10.1130/ges00217.1. Teza, G., Pesci, A., & Galgaro, A. (2008). Grid_strain and grid_strain3: Software packages for strain field computation in 2D and 3D environments. Computers & Geosciences, 34, 1142–1153. http://dx.doi.org/10.1016/j.cageo.2007.07.006. Tun, S. T., Wang, Y., Khaing, S. N., Thant, M., Htay, N., Htwe, Y. M. M., ... Sieh, K. (2014). Surface ruptures of the Mw6.8 March 2011 Tarlay earthquake, Eastern Myanmar. Bulletin

of the Seismological Society of America, 104, 2915–2932. http://dx.doi.org/10.1785/ 0120130321. Wang, X. W., Liu, G. X., Yu, B., Dai, K. R., Zhang, R., Chen, Q., & Li, Z. L. (2014a). 3D coseismic deformations and source parameters of the 2010 Yushu earthquake (China) inferred from DInSAR and multiple-aperture InSAR measurements. Remote Sensing of Environment, 152, 174–189. http://dx.doi.org/10.1016/j.rse.2014.06.014. Wang, Y., Lin, Y. N. N., Simons, M., & Tun, S. T. (2014b). Shallow rupture of the 2011 Tarlay earthquake (Mw6.8), Eastern Myanmar. Bulletin of the Seismological Society of America, 104, 2904–2914. http://dx.doi.org/10.1785/0120120364. Wessel, P., & Smith, W. H. F. (1998). New, improved version of generic mapping tools released. Eos, Transactions American Geophysical Union, 79, 579-579. http://dx.doi.org/ 10.1029/98EO00426. Wright, T. J., Parsons, B. E., & Lu, Z. (2004). Toward mapping surface deformation in three dimensions using InSAR. Geophysical Research Letters, 31, L01607http://dx.doi.org/10. 1029/2003GL018827. Yaseen, M., Hamm, N. A. S., Woldai, T., Tolpekin, V. A., & Stein, A. (2013). Local interpolation of coseismic displacements measured by InSAR. International Journal of Applied Earth Observation and Geoinformation, 23, 1–17. http://dx.doi.org/10.1016/j.jag.2012. 12.002.