Coseismic slip inversion based on InSAR arc measurements

We present a new method for inverting coseismic slip distribution based on arc measurements of InSAR interferograms. The method only solves the integer ambiguities on the selected arcs so that the challenging task from global unwrapping of low coherence interferograms can be avoided. The simulated experiment results show that the new method recovered the given slip distribution well at different coherence quality levels. However, the conventional method with global interferogram unwrapping fails when the interferogram has some isolated areas. In addition, the new method is capable of using surface rupture offset data gathered in the field. We apply the proposed method to study the 2010 Yushu, China Ms 7.1 earthquake. Inclusion of field data can help to enhance the results of fault slip inversion. It derives a maximum slip of∼ 3 m, larger than the published coseismic slip results on this event, but agreeing with the largest offset of 3.2 m from field investigation.


Introduction
Interferometric synthetic aperture radar (InSAR) coseismic deformation measurements have over the past two decades become one of the most important data sources for studying coseismic slip distributions.To apply InSAR measurement in coseismic slip inversion, most methods require unwrapping of the entire interferogram to get the displacement on each pixel.Current coseismic slip inversion studies with InSAR data mostly use the constraint from point displacements calculated from an interferogram (Simons et al., 2002;Wright et al., 2003;Wang et al., 2012a, b).
There are many InSAR phase unwrapping algorithms, among which the branch-cut region growing algorithm (Goldstein et al., 1988;Rosen et al., 1994) and the minimum cost flow (MCF) algorithm (Costantini, 1998) have been used most commonly.While most of the algorithms perform well for good-quality interferograms, it is difficult to unwrap an entire interferogram correctly with any of the algorithms when the coherence of the interferogram is low.One common problem is that some isolated regions would possibly be created and have no reliable estimates of phases from unwrapping processes (Goldstein et al., 1988;Chen and Zebker, 2000;Gens, 2003).
A solution to the problem is to use the wrapped phase values of interferograms for fault slip inversion directly, which has been proposed in recent years (e.g., Feigl and Thurber, 2009;Fornaro et al., 2011).A drawback of such algorithms is that the inversion problem becomes nonlinear.The computation is time-consuming and only a limited number of fault geometry parameters can be resolved (Feigl and Thurber, 2009).An improved version of the algorithms can resolve a slip distribution by adopting singular value decomposition (SVD) to cut down the number of slip parameters, where the decomposed parts with singular values lower than a fixed threshold are abandoned (Fornaro et al., 2011).The computation cost of the algorithm is however high and some information may be lost due to the use of the SVD.
We propose a method that is based on the use of phase measurements between point pairs (arc measurements in PSInSAR terminology) of an interferogram for coseismic slip inversion.The method unwraps the high-quality arcs but not the entire interferogram.It avoids the problems, e.g., signal loss or bias in isolated regions, from global unwrapping Published by Copernicus Publications on behalf of the European Geosciences Union.

C. Wang et al.: Coseismic slip inversion based on InSAR arc measurements
of interferograms while keeping the inverse problem linear.Additionally, the new method is capable of using surface rupture offset data gathered in the field because such data can be directly treated as arc measurements.Simulated experiments and the case of the 2010 Yushu Ms 7.1 earthquake will be used to validate and demonstrate the advantages of the method.

Method
Point displacement is the most direct way to describe the character of field deformation, but sometimes the absolute displacement on the point is hard to obtain, e.g., the displacement of points on the isolated fringes of the InSAR interferogram.Then we could also use the relative displacement between points to describe the deformation pattern.The conventional dislocation model links the point displacement with underground slip.Here, we make a modification to build a relationship between arc measurement and fault slip.In case that point displacement is not available, we could still use arc measurements to constrain the slip model.
The arc measurement between two pixels in a SAR interferogram can be written as where θ 1 and θ 2 are the wrapped phase observations of the two points, n is the differential phase integer ambiguity along the arc, and λ is the wavelength of the SAR sensor.With the locations of two points, we search the path between them on a flag map where branch cut and incoherent points have been marked out.As long as a path between the two points of an arc, which does not cross branch cuts or go through incoherent points, can be found, we consider it to be a high-quality arc and the integer ambiguity n of the arc can be easily resolved.This is particularly useful for an interferogram with many isolated regions where arc measurements can be unwrapped but the entire interferogram cannot.
The relationship between arc measurements d and point measurements d can be expressed by where A is a transformation matrix.Each row in A corresponds to a transformation from the displacement vector d to a point pair displacement d i .For example, considering the point pair displacement between points 3 and 5, d 35 , the corresponding row of A, takes the form [0 0 1 0 −1 0 0 . ..]where the third element is 1, the fifth element is −1 and all the other elements are 0. When assuming that Green's function linking the point displacements d and the slip vectors s is G, the relationship between d and s is d = AGs. (3) Considering that the point displacements are statistically independent and have the same variances σ 2 , the variance covariance matrix of d is then σ 2 I.According to the law of variance propagation, the variance covariance matrix of d is A Cholesky decomposition of AIA T is first carried out, A linear transformation of d can be followed, The variance covariance of d n is then A new set of observation equations can be formed: where G n = (L T ) −1 AG.The equations can be solved by adopting the widely used approaches for equal weight equation systems (e.g., Funning et al., 2005;Fialko, 2004).
It should be noted that the following condition should be satisfied for the m×n matrix A in order to be able to compute the Cholesky decomposition of AIA T , where m is the number of arcs.Cholesky decomposition requires AIA T to be a Hermitian, positive-definite matrix.Because the Gram matrix of a linearly independent vector is always positive-definite, AIA T is a positive-definite matrix when the row vector of A is linearly independent (the condition implied in Eq. 9).This means that in the network formed by the arcs, only one unique connection should exist between any two points.For example, Fig. 1a has some redundant arcs not satisfying the condition, while the cases in Fig. 1b and c do satisfy the condition.In an irrotational field with all the residues balanced by branch cuts, the unwrapped result for an arc is independent of the path; therefore, the observations on these redundant arcs can be derived from other arcs observations.This means that such redundant arcs do not contribute to the inversion but bring additional computational cost, so it is necessary to abandon them during the inversion, not only to fit the Cholesky decomposition condition but also to cut the computational burden.

Simulated experiments
To demonstrate the advantages of the proposed method, we attempt to recover a given slip distribution separately with our method and the conventional method.Here an ascendingpass Envisat interferogram due to a given fault slip model (Fig. 2a) was calculated by forward modeling (Fig. 2 and  c).We simulated four coherence maps with different quality levels and generated the corresponding Gaussian noises.The variance of noise is in inverse proportion to the cubic of coherence.Figure 2d shows the interrupted interferograms after adding noises (the noise level increases from 1 to 4).Branch cut and MCF algorithms are separately adopted to unwrap the interferograms (Fig. 2e and f).The regions with coherence lower than 0.6 are masked out in both algorithms.
To assess the unwrapping quality, we calculated the difference between unwrapped displacement and initialized displacement (Fig. 2h and 2i).A stable difference map would imply a successful unwrapping process, because the correct unwrapped map is equal to the initialized displacement map minus the reference point displacement.Figure 2h shows that the branch cut algorithm gives a satisfactory result at different noise levels, but the problem is that many isolated regions, which are blocked from the reference point by incoherent areas or branch cuts, are not unwrapped.Such regions become more when the noise increases.The MCF algorithm can unwrap all the isolated regions, but MCF is very likely to give a wrong unwrapped value in the isolated regions.Figure 2i shows that the unwrapped values differ a lot from the initialized displacement, which is even terrible at higher noise levels.Most of the error is due to the miscalculation of the offset between the isolated region and the reference point.Therefore, if we implement global unwrapping in the lowquality interferogram, the isolated regions are either blocked off (by the branch cut algorithm) or unwrapped mistakenly (by the MCF algorithm).
The above global unwrapping problems happened mainly because it is hard to estimate the offsets between isolated regions and reference points.When using the arc to constrain the slip model, we directly adopt the relative displacement inside the isolated regions but not the absolute displacement on the point.Therefore it is not required to know the offsets between isolated regions and reference points, and the global unwrapping problems described above do not exist.Figure 2g shows that arcs can be constructed in all isolated regions as long as there are satisfied fringes.The increas-ing noises only remove the arcs over the noisy region, while not affecting the arcs inside isolated regions.To apply the arc constraint to coseismic slip inversion, there are generally four steps as follows: (1) a quadtree algorithm is first used to down sample the interferogram with a given threshold (Jonsson et al., 2002).( 2) A network linking the sampled points is constructed by using a local Delaunay triangulation (Zhang et al., 2011).The arcs that pass through branch cuts and points with coherence lower than a given threshold are removed from the network; (3) the redundant arcs are removed by a minimum spanning tree algorithm (Kruskal, 1956) to satisfy the condition in Eq. ( 9). ( 4) The Laplacian smoothing constraint is added to the equation system, and a least squares method is used to solve the equations to obtain the slip distribution.
Figure 3 shows the derived slip model with constraints separate from arcs, branch cut unwrapped points (BP) and MCF unwrapped points (MP).In a high-quality interferogram (Level 1), all of them suggest very similar results close to the given slip.However, when the noise gradually increases, the three methods have different performances.The slip inversion with arc constraints performs the best, where two slip concentrations can still be recognized even in Level 4. The BP constraint derives a much lower resolution slip than the arc constraint.The derived north slip is hard to recognize in Level 2. In levels 3 and 4, the slip is over-smoothed into one concentration and the maximum slip decreases a lot, much less than the given slip value.MP still performs satisfactorily in Level 2; however, in Level 3 and Level 4, the slip distribution is biased.Its location and value deviate a lot from the given slip.
When interferogram coherence becomes worse, the arcs across the noisy regions are removed and the slip resolution is weakened partly.However, the arcs in the isolated regions, which contain many deformation signals, are still left (as seen in Fig. 2g).The general slip pattern could therefore be retrieved successfully by using an arc constraint.The number of BP drops significantly along with the decrease in coherence, mostly because the isolated regions are blocked off.Consequently the decreasing constraints would reduce the resolution of the result.The MCF algorithm has estimated the offset between isolated regions and the reference point.In  Level 2, the estimation did not deviate too much, so the result is generally correct.However, in Level 3 and Level 4, the estimation deviates dozens of centimeters, even reaching the maximum displacement, and therefore the derived slip would be destroyed.

Application to the Yushu earthquake
The proposed method is now applied to study the coseismic slip of the 13 April 2010 Yushu, China Ms 7.1 earthquake.This event occurred on the left lateral Ganzi-Yushu fault, the western part of the Yushu-Garze-Xianshuihe fault zone, causing the death of around 2700 people.A coseismic deformation interferogram is formed by using ascending ALOS PALSAR images acquired on 15 January 2010 and 17 April 2010, respectively (Fig. 4a).The interferogram can be satisfactorily unwrapped as the coherence of the interferometric pair is high over the entire interferogram.
5173 points are sampled by the quadtree algorithm with a threshold of 3 rad.As shown in Fig. 4c, 21 108 arcs have been built based on the 5173 points.The arcs that pass through branch cuts and points with a coherence lower than 0.6 are then removed from the network.The integer ambiguities along the 11 533 remaining arcs (Fig. 4d) are then obtained.5085 arcs are finally left after applying the minimum spanning tree algorithm to remove the redundant arcs (Fig. 4e).A least squares method is used to solve the equations to obtain the slip distribution as shown in Fig. 5a.The fault geometry used is adopted from Li et al. (2011).
The conventional approach to the fault slip geometry inversion can also be applied based on the surface deformation field determined from the interferogram as shown in Fig. 5b.A comparison of the results shows that the proposed method has produced almost the same slip distribution as that produced by the conventional method (Fig. 5a and b).The largest slip patch is in the eastern segment at a depth of about 4 km.A slight difference between the two sets of results is the maximum slip, 2.27 m from the proposed method and 2.33 m from the conventional method.The results are also very close in general to those given by Li et al. (2011), but differ somewhat from those given by Qu et al. (2012) and Tobita et al. (2011).The reason for the different results is considered to be mainly due to the fact that three fault segments are assumed in this study and Li et al. (2011), but two were used in the other two studies.
The interferogram used for the Yushu case study has good coherence.There are almost no isolated regions blocked off during phase unwrapping.Therefore, conventional methods with constraints from unwrapped points can already perform well.Using the new method, we derived almost the same slip distribution as a conventional one.Additionally, the result conforms well to another published result using the same fault geometry.Therefore, it is believed that the new method is a reliable way to derive slip models from InSAR interferograms.
After the occurrence of the earthquake a field investigation was carried out to measure the surface rupture displacements (Lin et al., 2011).The displacements were measured at 54 locations using a tape measure and an Advantage Laser Rangefinder, with an error of ±15 cm.The surface rupture distributed from the epicentral area to the eastern end of the Changgu Temple segment with a total length of 51 km (Fig. 5c and d).Surface rupture offset measurements are rarely used for constraining coseismic slips because they are not point observations, which are required for the conventional inversion method.In the Yushu case, these data are also not incorporated into the existing published coseismic slip studies.However, it is possible to incorporate directly such rupture offset measurements into the solution when the proposed method is applied by considering the rupture displacements as corresponding arc measurements.The computed slip distribution when incorporating the rupture displacements is shown in Fig. 3c.It can be seen from the results that the slip location has changed slightly but the maximum slip has increased to 3 m, very close to the largest offset of 3.2 m measured from the field investigation.Comparatively, Field investigation is often carried after a significant earthquake.Most studies only use the field rupture location to constrain the fault rupture trace, while rarely adopting the rupture displacement into the dislocation model because it is not a point observation.However, the surface rupture offset represents the displacement closest to the fault, where there is normally no InSAR record due to coherence loss; therefore it would have brought a strong constraint onto the slip model.The new method could easily handle such data and include it in the inversion scheme.In this case, the maximum slip has been more reasonable after including surface offset data.It implies that the new method could help to enhance the geodetic slip model by introducing surface rupture displacement constraints.

Conclusions
A new approach has been proposed for fault slip inversion based on arc measurements in SAR interferograms.A distinctive advantage of the method is that there is no need to apply global unwrapping of the interferogram.The method is especially useful for interferograms with isolated regions that are blocked from the reference point in global unwrapping.In addition, the method also allows field measurements of slip rupture displacements to be incorporated into the solution easily.
The proposed method has been applied to simulated experiments and the realistic case of the 2010 Yushu Ms 7.1 earthquake.The following main conclusions can be drawn from the results: 1.As indicated by both simulated tests and the case study, when the quality of an interferogram is high, the new method produces similar results as those from the conventional method that is based on point displacements from an unwrapped interferogram.
2. As the new method does not require global unwrapping of an interferogram, the possible errors brought by it can be avoided.The simulated test has shown that the new method performs much better than conventional ones in dealing with the low-quality interferograms where many isolated regions exist.
3. Surface rupture displacements gathered from field investigations can be very useful in enhancing the accuracy of slip inversion.Such data can be easily incorporated into the solution of the proposed method.After including surface offset data, the maximum slip of the Yushu earthquake reaches 3 m, much closer to the field investigation than the existing coseismic slip studies.

Fig. 1 .
Fig. 1.Networks of arc observations in a SAR interferogram.(a) Network with redundant arcs where there is more than one path to connect some of the point pairs.The rank of the m×n matrix A is smaller than m.(b, c) Networks without redundant arcs where there is only one path to connect any two points.The rank of A is equal to m.

Fig. 2 .
Fig. 2. (a) Assigned initial slip distribution.(b) Simulated displacement in ascending Envisat line of Sight.(c) Simulated interferogram.(d) The interferograms after adding different level of noises.The noises gradually increase from Level 1 to Level 4. (e) The unwrapped results from the branch cut algorithm.(f) The unwrapped results from the MCF algorithm.(g) The arcs used for inversion.(h) The difference between displacement unwrapped by the branch cut algorithm and the true displacement.(i) The difference between displacement unwrapped by the MCF algorithm and the true displacement.

Fig. 3 .
Fig. 3. Slip distribution results from different constraints.Left to right slip distributions correspond to the Level 1 to Level 4 interferograms.(a) Arc constraint.(b) BP constraint.(c) MP constraint.

Fig. 4 .
Fig. 4. (a) ALOS PALSAR coseismic deformation interferogram over Yushu.(b) 5173 points after down sampling by using the quadtree algorithm.(c) Local Delaunay triangulation network.(d) Arcs after solving for the phase integer ambiguities.(e) Network after applying the minimum spanning tree algorithm.

Fig. 5 .
Fig. 5. Results of coseismic slip distribution inversion.(a) Results from arc measurements.(b) Results from point displacements.(c) Results from arc measurements and surface rupture displacements.Red points indicate the 54 surface rupture locations recorded in Lin et al. (2011).(d) Coseismic surface rupture displacements given by Lin et al. (2011).The upper lines denote horizontal displacement.The lower lines denote vertical displacement.