Measuring Coseismic Deformation With Spaceborne Synthetic Aperture Radar: A Review

25 space-borne Synthetic Aperture Radar imagery has become an increasingly available data source for the study of crustal deformation associated with moderate to large earthquakes (M > 4.0). Coseismic surface deformation can be measured with several well-established techniques, the applicability of which depends on the ground displacement pattern, on several radar parameters, and on the surface properties at the time of the radar acquisitions. The state-of-the-art concerning the measurement techniques is reviewed, and their application to over 100 case-studies since the launch of the Sentinel-1a satellite is discussed, including the performance of the different methods and the data processing aspects, which still constitute topics of ongoing research.


INTRODUCTION
Coseismic deformation measurements based on spaceborne SAR sensors have been carried out since 25 years (Massonnet et al., 1993). Furthermore, high-quality SAR data suitable for crustal deformation studies, in particular with low spatial and temporal baselines and precise orbital and timing control, has never been so abundant since the launch of the Sentinel-1a and Sentinel-1b satellites in Apr. 2014 and 2016 respectively, and of ALOS-2 in May 2014.
The maturity of the processing techniques together with the quality and availability of SAR data have contributed to the standardization of several processing steps and to the development of automated processing chains for Sentinel-1 data Li et al., 2016b). However, a given application and/or dataset can still present some challenges for the current state-of-the-art processing techniques, or require ad-hoc decisions on behalf of the image analyst. Furthermore, some specific properties as well as the wide-area coverage of the main Sentinel-1 and ALOS-2 acquisition modes, and the free and open data policy of the Sentinel-1 missions, has stimulated new algorithmic developments, which are partially still ongoing.
The measurement principles and state of the art of current SAR-based measurement techniques used for coseismic deformation applications are discussed in section SAR Deformation Measurements. A review of coseismic deformation studies concerning earthquakes, which occurred after the launch of Sentinel-1, is then carried out, discussing the performance of the measurement techniques (section Technique Application) as well as the data processing aspects, which are still challenging and not fully standardized (section Current Challenges). Finally, future developments are discussed in section Future Developments.

Sensitivity to Motion
SAR-based ground deformation measurements rely on repeat-pass acquisitions carried out by C-, L-, and Xband microwave SAR sensors operating in Stripmap and ScanSAR acquisition modes (Cumming and Wong, 2005) and, since the launch of Sentinel-1a, also on imagery acquired in Terrain Observation with Progressive Scans (TOPS) mode (De Zan and Monti Guarnieri, 2006). An overview of spaceborne SAR sensors and modes, which have been used for coseiemic deformation measurements since the launch of Sentinel-1a is provided in Table 1. The properties of other sensors and acquisition modes are reviewed in Sansosti et al. (2014).
A schematic of the above-mentioned acquisition modes is shown in Figure 1. In Stripmap mode, the SAR maintains a fixed side-looking beam pointing throughout a data take, yielding spatial resolutions in the order of a few meters and image sizes of tens of km in both ground range and azimuth ( Table 1). ScanSAR and TOPS modes trade-off azimuth resolution, which is degraded to tens of meters, for an increased range coverage, which can reach several hundreds of km (Table 1). This is achieved by illuminating the ground with a fixed antenna beam orientation only for the duration of a so called data burst, after which the beam is electronically steered in range to cover an adjacent sub-swath. In the TOPS case the antenna is also steered in the azimuth direction within each burst, mainly to reduce an undesired intensity modulation known as image scalloping (Meta et al., 2008). In the ScanSAR case, and to a lesser extent also in the TOPS case, an azimuth overlap exists between consecutive bursts of the same swath, which can be exploited for deformation measurements as discussed in section Split-Bandwidth Interferometry.
For all acquisition modes sensitivity to motion arises from the positioning information contained in the amplitude and phase of focused SAR imagery. The data of current spaceborne sensors is distributed in zero-Doppler geometry, implying that a point on ground is positioned about its zero-Doppler coordinates (η 0 , τ 0 ) in the SAR image. Coordinate η 0 represents the time at which the derivative of the azimuth phase modulation between consecutive pulses backscattered from the same object on ground is zero, whereas τ 0 represents the two-way range-delay at time η 0 . In an ideal scenario, and in particular neglecting the atmospheric effects discussed in sections Tropospheric Propagation and Ionospheric Propagation, the Doppler coordinates have a geometric interpretation shown in Figure 1A, in which η 0 represents the azimuth time of closest approach between a point on ground and the SAR antenna phase center (IEEE Standard Definitions of Terms for Antennas, 1993) and τ 0 = 2r 0 /c , where r 0 is the one-way range of closest approach and c is the speed of light.
Due to the finite resolution of the SAR in both the range and azimuth dimensions, a point on ground (target) actually covers a two-dimensional locus h(τ ,η) in the SAR image, which for lowsquint systems can be written as Holzner and Bamler (2002); Cumming andWong (2005), andDe Zan et al. (2014): where τ and η represent the slant-range delay and azimuth time variables, A and ψ are amplitude and phase terms related to the electromagnetic properties of the target, p r () and p a () represent sharply peaked (sinc-like) functions, f c is the carrier frequency of the radar, τ phase 0 is the range phase delay which differs slightly from the range group delay τ 0 due to ionospheric propagation (section Ionospheric Propagation), and f dc is the instantaneous Doppler centroid. The latter is a slowly varying function of (η 0 , τ 0 ) in the Stripmap case, whose values are typically <100 Hz in magnitude for current systems. For a single ScanSAR burst (Holzner and Bamler, 2002) and for TOPS acquisitions (De Zan et al., 2014) the Doppler centroid also depends on the azimuth time distance from the burst center (T c ). In the ScanSAR case values in the order of 1 KHz in magnitude can be reached, whereas for Sentinel-1 TOPS the maximum instantaneous Doppler centroid magnitude is around 2.6 KHz (Yague-Martinez et al., 2016).
Considering now two SAR images acquired respectively before and after an earthquake, based on Figure 1A, a 3D coseismic deformation, represented by its components ( e, n, u) in a local East-North-Up Cartesian system, will cause the following zero-Doppler deformations in the range and azimuth directions ( r defo , a defo ), counted positive respectively away from the radar and along the flight path: where θ is the local incidence angle counted positive from the vertical (U), and α is the angle between North (N) and the ground-projection of the flight direction, counted positive counter-clockwise. Typical incidence angle values are between 20 deg and 50 deg (Table 1), whereas α is in the order of 10 deg to 15 deg at mid-latitudes. r defo is more commonly referred to in literature as the lineof-sight (LoS) deformation. For a given point, such deformation components will contribute to the following zero-Doppler time coordinate variations between the pre-and post-event image: where V r is an effective rectilinear velocity (Cumming and Wong, 2005), which is in the order of 7100 m/s for remote sensing satellites in near-polar orbit.

Measurement Techniques
Three SAR-based techniques are currently used for coseismic deformation measurement, namely Differential SAR Interferometry (DInSAR), which measures the LoS deformation component, offset-tracking, which provides both the LoS and the azimuth components, and Split-bandwidth Interferometry (SBI), which can be separately applied to the range and azimuth dimensions to yield respectively the LoS and the azimuth deformation. These methods, which are detailed in the following sub-sections, measure the displacement occurred between two SAR acquisitions, and the permanent coseismic deformation is typically obtained by applying these techniques to a single image pair spanning the earthquake. Methods based on a redundant network of acquisition pairs, namely Multi-Temporal DInSAR (MT-DInSAR) techniques, have been used in a limited amount of cases to reduce measurement noise (e.g. Fattahi et al., 2015;Fielding et al., 2017;Grandin et al., 2017;Lee et al., 2017). In section Recovering the 3D Deformation Field methods to obtain the 3D deformation field from suitable LoS and azimuth measurements are discussed.

Differential SAR Interferometry
DInSAR exploits the range-dependent phase term in (1) to measure the LoS deformation component. Its application to Stripmap SAR imagery has been the topic of several review papers (Bamler and Hartl, 1998;Massonnet and Feigl, 1998;Bürgmann et al., 2000;Rosen et al., 2000;Zhou et al., 2009;Simons and Rosen, 2015). The main processing steps, based on developments of the original two-pass approach proposed by Massonnet et al. (1993), are shared by virtually all processing chains (cfr. Figure 2 in Bürgmann et al., 2000). These include image coregistration, interferogram formation and filtering, removal of Frontiers in Earth Science | www.frontiersin.org 3 February 2019 | Volume 7 | Article 16 the topographic phase contribution using an external Digital Elevation Model (DEM), 2D phase unwrapping and conversion to displacement. The latter step may optionally preceded by a calibration procedure to estimate and remove image-wide error trends and refer the measurements to a reference point or set of points of known motion (e.g., GPS stations or areas which are assumed to be stationary). Considering a pixel located at zero-Doppler coordinates τ 1 0 , η 1 0 in image 1, and whose corresponding position in image 2 is τ 2 0 , η 2 0 , orbital and timing uncertainties, as well as ionospheric propagation effects discussed in section Ionospheric Propagation, will allow only estimates of these quantities to be determined, namely τ 1 0 ,η 1 0 and τ 2 0 ,η 2 0 . If a deformation in the LoS and azimuth directions τ defo , η defo has caused the actual location in image 2 to be τ 2 0 + τ defo , η 2 0 + η defo , the DInSAR interferogram pixel value I DinSAR will be: in which: conj() represents complex conjugation; τ phase is the range phase delay difference (7), due to deformation τ defo , tropospheric propagation τ tropo , and ionospheric propagation τ phase iono f c ; η is the azimuth position difference (8), due to the underlying deformation η defo , and to the orbital uncertainty ( η orb ) and ionospheric propagation contributions (η iono ) to the overall azimuth co-registration error η coreg , which is given by (9). It is further assumed that in the differential interferogram formation process the flat earth and topographic contributions were perfectly removed. Any residual errors would appear as additional error terms in (7).
The prerequisite for DInSAR is interferometric coherence, which occurs when the target-dependent amplitude and phase terms, namely A and ψ in (1), are statistically similar in the two acquisitions. In this case the corresponding phase difference ψ can be considered a random variable with a distribution, which is sharply peaked around a null mean value (Bamler and Hartl, 1998), so that in the error-free case, the LoS deformation r defo can be related to the phase of a DInSAR interferogram as follows: where λ is the radar wavelength, arg() is a function returning the argument of a complex number in the (-π,π] interval and unwrap() is a function performing phase unwrapping. The last equality in (10) holds only if the Doppler-centroid dependent phase term in (6) is negligible, which is always the case for Stripmap imagery, whereas as discussed below, this assumption may not hold for pixels at the edge of TOPS bursts (large f dc ) in the presence of significant horizontal deformation, particularly in the north-south direction [large η based on (4) and (5)]. DInSAR application to burst-mode and full-aperture ScanSAR image pairs (Holzner and Bamler, 2002) as well as to ScanSAR-Stripmap imagery (Bertan-Ortiz and Zebker, 2007), was successfully demonstrated for several past missions, such as Radarsat-1, ENVISAT, and ALOS-1 (Bamler and Holzner, 2004;Guccione, 2006;Motagh et al., 2008;Tong et al., 2010;Xu et al., 2011;Fielding et al., 2013), even though these had not been originally designed to support ScanSAR interferometry. Currently ScanSAR DInSAR is highly relevant for coseismic applications, in particular due to the routine observations carried out by the ALOS-2 sensor, which maintains a high level of burst alignment since Feb. 2015 (Natsuaki et al., 2017). Specific DInSAR processing issues related to the full-aperture product of this sensor are discussed in Liang and Fielding (2017a).
In recent years further algorithmic developments were required to apply DInSAR to the Interferometric Wideswath (IW) TOPS imaging mode of the Sentinel-1a and Sentinel-1b sensors, for which the azimuth beam steering introduces a significant coupling between azimuth image registration and interferometric phase (Prats-Iraola et al., 2012) through the azimuth-dependent phase term in (6). Methods have been developed to meet the much more stringent azimuth coregistration requirements compared to the Stripmap case (Yague-Martinez et al., 2016) and adaptations have been proposed to handle spatially varying azimuth co-registration errors related to orbital/timing uncertainties over long data takes  and (presumably) ionospheric propagation effects (Wang et al., 2017c). As previously mentioned, the TOPS DInSAR phase is also sensitive to azimuth motion, through the η defo and η terms in (8) and (6) respectively. While on one side it has rightly been observed that this sensitivity is just a property of the TOPS acquisition mode (De Zan et al., 2014), in the presence of significant azimuth motion it may cause phase discontinuities at the burst boundaries due to the change of sign of and f dc in (6), which in turn complicates the phase unwrapping step (González et al., 2015;Scheiber et al., 2015). However, as reviewed in section 3, so far this did not prevent application of TOPS imagery also to very large earthquakes.
For each imaging mode, namely Stripmap, ScanSAR and TOPS, DInSAR is the measurement technique, which typically provides the highest accuracy, which can reach a fraction of the centimetric SAR wavelength, provided interferometric coherence between the two acquisitions is retained and the error sources discussed in section Current Challenges are not dominant compared to the deformation. Concerning coherence, for current sensors the main limitations for coseismic applications are often surface changes due to land cover (e.g., vegetation and water bodies) and weather (e.g., snow in mountainous regions), geometric distortions (radar Frontiers in Earth Science | www.frontiersin.org 4 February 2019 | Volume 7 | Article 16 layover and shadow) due to topographic relief, as well as high strain rates, such as those typically observed close to fault surface ruptures (Massonnet and Feigl, 1998;Bürgmann et al., 2000). In particular if the strain rate exceeds a deformation of λ/4 within a SAR resolution element, either in the range or in the azimuth dimension, coherence will be completely lost.
The achievable spatial resolution of DInSAR LoS measurements in each image dimension is in the order of several radar resolution cells, which in turn are between 2 and 80 m in each image dimension for current SAR acquisition modes relevant for coseismic studies (Table 1). Finally, an important property of DInSAR measurements, regardless of the imaging mode, is that they are spatially relative. Depending on the processing approach, the output deformation map is either referred to a reference pixel with respect to which the phase gradients were unwrapped, or to the average motion of the tie-points used to calibrate the measurements after the phase unwrapping step.
A first DInSAR application example is shown in Figure 2, where two Stripmap InSAR pairs are considered, namely a descending one acquired by COSMO-SkyMed-1/2 on 02 Feb. 2014 and 10. Feb. 2014 respectively, and an ascending image pair acquired by the TanDEM-X satellite on 28 Jan. 2014 and 08. Feb. 2014. Both pairs span the Mw 5.9 Feb 3rd, 2014 Cephalonia, Greece, earthquake (Merryman Boncori et al., 2014). Considering the descending dataset, Figure 2A shows a coherence map prior to phase filtering, in which low coherence levels correspond to vegetated and cultivated areas, to surface ruptures in the northern part of the Paliki peninsula, and of course to water bodies. Figure 2B shows the wrapped and filtered DInSAR phase, in which each phase cycle represents a LoS deformation of half a radar wavelength based on (10), which amounts to 1.55 cm for the X-band COSMO-SkyMed radars ( Table 1). Figure 2E shows the unwrapped LoS DInSAR deformation, in which data gaps are present for areas whose phase noise levels are too high even after phase filtering, and which can therefore not be reliably unwrapped. In this particular case, the coseismic deformation field causes a DInSAR phase discontinuity of several cycles in the northern part of the Paliki peninsula, so that phase unwrapping had to be performed separately for the western and eastern sectors of the image, using respectively R1 and R2 as spatial references. As detailed further in section 4.3, offset-tracking results were then used to recover the deformation gradients between these reference points and generate the LoS deformation map shown in Figure 2E.
Further examples of DInSAR interferograms at are shown in Figures 3A,E, 4A for the C-band Sentinel-1a sensor, operating in TOPS mode, and in Figures 5B,D for L-band ALOS-1 PALSAR sensor, acquiring in Stripmap mode. These examples will be referred to further in the following sections. At this stage it is just worth noting that the phase cycles (fringes) represent half a wavelength of LoS deformation (or more in general phase delay), which amounts to 1.55 cm for Figures 3B,D (Xband), to 2.8 cm for Figure 4A (C-band) and to 11.8 cm for Figures 5B,D.

Offset Tracking
Offset-tracking methods, also referred to as pixel-offset or pixel-tracking techniques, exploit the range-and azimuthdependent amplitude terms in (1) to measure both the LoS and the azimuth displacement components between two SAR acquisitions. This is done by estimating the 2D misregistration (i.e., group delay change) of corresponding pixels by applying several image matching techniques (Brown, 1992) to track amplitude (or intensity) features and/or coherent SAR speckle. Gray et al. (1998) first demonstrated the use of amplitude cross-correlation for ice velocity measurements, whereas Michel et al. (1999) and Peltzer et al. (1999) first carried out coseismic deformation measurements with phase correlation and normalized cross-correlation respectively. These latter techniques are still currently the most widely used, although recently some algorithmic variations have also been proposed (e.g., Wang and Jónsson, 2015).
An offset tracking measurement can be considered the result of a search of the registration shifts, which maximize a similarity parameter, e.g., Normalized Cross Correlation (NCC): Compared to the differential phase delay in (7), the group delay has an inverted ionospheric contribution (section Ionospheric Propagation), and an additional term related to LoS position uncertainties (τ orb ), due to orbit inaccuracies: The azimuth mis-registration η is the same as in (8). Based on (5), in the error-free case the LoS deformation r defo and azimuth deformation a defo can be related to the respective range and azimuth mis-registration as follows: The achievable accuracy of offset-tracking techniques is typically between 1/10th and 1/100th of the spatial resolution in each SAR image dimension, and thus at least an order of magnitude worse than for DInSAR. Accuracy depends on interferometric coherence, if coherent speckle is being tracked, or on the prominence of common features if these are present (Bamler and Eineder, 2005). The spatial resolution of the measurements is in the worst case equal to twice the matching window size, which is typically tens of radar resolution cells in each image dimension, and in the best case equal to the distance (posting) between neighboring measurements, which is typically chosen to be a fraction of the matching window size (Pritchard et al., 2005). Thus, compared to DInSAR, also the spatial resolution of the measurements is at least an order of magnitude worse.  Despite the above mentioned drawbacks, offset tracking also offers some advantages and complementarities with respect to DInSAR: the measurement of the azimuth deformation component; the possibility of succeeding also in the absence of interferometric coherence, provided trackable amplitude features are present; a greater algorithmic robustness, due to the fact of not requiring coregistration nor phase unwrapping, which may be delicate processing steps for specific sites and datasets (section Phase Unwrapping).
An example of the complementarities of offset-tracking with respect to DInSAR is shown in Figure 2 for the Mw 5.9 Cephalonia, Greece, event. Concerning the LoS measurements, the coarse spatial resolution and higher noise level of offset tracking compared to DInSAR is apparent comparing Figure 2I with Figures 2E,K with Figure 2G. For these Stripmap datasets ( Table 1) the spatial resolutions of DInSAR measurements after phase averaging are between 10 and 30 m for the descending and ascending datasets respectively, and between 100 and 800 m for the offset-tracking measurements. However, offset-tracking also provides LoS measurements where the DInSAR phase cannot be unwrapped, since amplitude features provide cross-correlation peaks in (11), even in the absence of interferometric coherence. Furthermore, offset tracking provides also the azimuth motion components (Figures 2J,L), which in this case are essential to reveal the complexity of the fault geometry (Merryman Boncori et al., 2014). Finally, it should be noted that although co-registration is not required by offset-tracking, any image-wide coregistration errors, e.g., due to orbital and timing uncertainties and/or atmospheric propagation, will translate directly into a bias in the measurements. As for DInSAR, a tie-point calibration procedure based on points of known motion or areas assumed to be stationary can be used to compensate such effects.

Split-bandwidth
Interferometry techniques measure displacement in each dimension through a double difference interferometric procedure. Interferograms between corresponding upper and lower spectral subbands in each image are formed (I u and I l ), followed by a second difference interferogram, the unwrapped phase of which is proportional to the sought displacement. For the TOPS mode case, the images should be deramped to remove the phase modulation due to the azimuth antenna scanning before forming the sub-band interferograms (I u and I l ).
In the azimuth dimension, assuming an upper and a lower Doppler frequency sub-bands centered about f u dc and f l dc respectively, the following relations hold:  from which the azimuth deformation is obtained in the error free case as: In the range dimension, assuming upper and lower rangefrequency sub-bands centered around f u c and f l c respectively, the corresponding split-bandwidth interferograms are given by: For TOPS imagery, unlike the azimuth-SBI case, it is not necessary to deramp the input images before forming the subband interferograms in (18) and (19). Neglecting ionospheric propagation effects (section Ionospheric Propagation), the LoS deformation is given by: SBI was originally proposed in the range dimension for absolute phase estimation (Madsen and Zebker, 1992) and later both in the range and in the azimuth dimensions for Stripmap image co-registration (Scheiber and Moreira, 2000), where it was named Spectral Diversity. Its potential for geophysical applications, was recognized following the work of Bechor and Zebker (2006), who reformulated azimuth Spectral Diversity in terms of squint angle rather than of Doppler frequency, dubbing their method Multi Aperture Interferometry (MAI), and applying it to the Mw 7.1 1999 Hector Mine event.
In the latter study, the azimuth sub-bands were obtained by dedicated raw data focusing steps, whereas the same method was applied to conventionally focused Stripmap imagery by Barbot et al. (2008)  use of an adaptive phase filtering algorithm widely used for DInSAR (Goldstein and Werner, 1998). Bamler and Eineder (2005) derived theoretical upper bound and empirical curves for the accuracy of the technique as a function of coherence and phase estimation window size (multi-looking factors), as well as the optimum frequency separation between sub-bands. De Zan et al. (2015) discussed the effects of multi-looking the subband interferograms prior to formation of the final difference interferogram (referred to as early multi-looking) as opposed to forming the sub-band interferograms at a full resolution and averaging only the final difference interferogram (referred to as late multi-looking). Recent developments of azimuth split-bandwidth interferometry addressed its application to Sentinel-1 TOPS and ALOS-2 ScanSAR imagery. Concerning Sentinel-1, spatially discontinuous measurements were obtained by applying the method to the burst-overlap regions Jiang et al., 2017a,b), which have a "naturally" large Doppler frequency separation due to the azimuth antenna steering, and thus a high motion sensitivity (10). This method will be referred to as Burst-Overlap MAI (BO-MAI) in the following. Jiang et al. (2017a) also reported the first application, to the author's best knowledge, of range-SBI, which was demonstrated for the 2016 Kumamoto sequence on Sentinel-1 TOPS data, the large range bandwidth of which allows a high sub-band frequency separation and thus good motion sensitivity (20). Concerning ALOS-2 ScanSAR, Liang and Fielding (2017b) applied azimuth SBI to data acquired from different bursts [with Doppler frequency separation given by (3)], after extracting these from the full-aperture product, to obtain a continuous azimuth deformation map for the 2015 Mw 7.8 Gorkha event. This technique can be considered the ScanSAR equivalent of BO-MAI, and will be referred to as such in the following.
From the performance point of view, SBI shares the same interferometric coherence requirements as DInSAR, being itself an interferometric method, including the limitations due to the maximum observable strain rates. However, when applied to Stripmap images or single ScanSAR or TOPS bursts, phase unwrapping is almost never required due to the limited frequency separations in (19) and (20), which in turn can lead to improved spatial coverages compared to DInSAR. Compared to offsettracking, SBI typically provides greater accuracies and spatial resolutions in coherent areas, whereas being a phase-based method it cannot carry out measurements in completely decorrelated areas.
For the azimuth-SBI case (MAI), an example of these considerations is shown in Figure 2. The spatial coverage of the MAI measurements (Figures 2F,H) has less data voids compared to their DInSAR counterparts (Figures 2E,G), since some areas had to be masked out to avoid phase unwrapping errors. Compared to offset-tracking (Figures 2J,L), an increase in spatial resolution can be noticed. The latter is about 60 m for MAI compared to 100-800 m for offset-tracking as mentioned in section Offset tracking. The lower noise level compared to offset-tracking is also apparent.
Concerning the range-SBI case, an application example is shown in Figure 3 for a Sentinel-1a TOPS image pair used to study the Mw 7.1 2016 Kumamoto, Japan, event. The coverage difference with respect to DInSAR is very significant in this case (cfr. Figures 2F,E), given the fringe density close to the fault (Figure 2A), which causes the phase unwrapping to be problematic. In contrast, where DInSAR is applicable, its significantly higher spatial resolution and reduced noise level compared to range SBI is obvious. Compared to LoS offset-tracking, Figure 2G, range-SBI provides an increased spatial resolution as well as a reduced noise level. However, offset-tracking provides an even better coverage close to the fault, since range-SBI fails where coherence is lost, e.g., due to the above-mentioned maximum strain-rate constraints.
Finally, an example of TOPS BO-MAI is shown in Figures 3D,H. The spatial coverage is discontinuous, although compared to conventional MAI applied to each TOPS burst ( Figure 2C), higher accuracies and resolutions can be achieved due to the much higher Doppler frequency separation (17), which is between 4.4 and 5.2 KHz in the BO-MAI case and at most 150 Hz in the conventional MAI.

Recovering the 3D Deformation Field
Based on (4), each LoS and azimuth deformation measurement provides one equation in 3 unknowns, namely the Cartesian components of the underlying 3D deformation. Provided at least 3 linearly independent LoS and/or azimuth components are available, a Weighted Least Squares (WLS) problem may be formulated and solved (Wright et al., 2004), in which the weights are typically based on the variance of the measurements in the deformation far-field (Funning et al., 2005). Improvements to the WLS approach to account for the spatial correlation among measurements have been recently proposed . Regardless of the inversion approach, observables are provided either by LoS and azimuth observations on ascending and descending passes of right-looking acquisitions (e.g., Jo et al., 2017;Morishita et al., 2018;Wang et al., 2018b;Xu et al., 2018b), or by LoS observations alone, taken from left-and right-looking acquisitions .
An application example of the WLS method is shown in Figure 2, where LoS and azimuth measurements available from DInSAR, MAI and offset-tracking (Figures 2E-L) are jointly inverted to solve for the Cartesian deformation components (Figures 2M-O).
For cases in which only DInSAR observations of the LoS deformation r defo are carried out from ascending and descending passes, which is often the case (e.g., Solaro et al., 2016;Wen et al., 2016;Cheloni et al., 2017), the so-called 2.5D deformation field can be computed, which consists in the east and vertical deformation components, respectively e and u, computed by assuming no north contribution ( n = 0). This assumption is motivated by the lower sensitivity of LoS measurements to the north motion component due to the nearpolar satellite orbit, which can be quantified for a given SAR image pixel using (4). Assuming a mid-latitude value of α = 10 • , and θ = 38 • as in the center of a Sentinel-1 IW2 swath, the From (21) it is seen that LoS measurements are 6 to 7 times more sensitive to east and vertical motion respectively, compared to the north direction. Weston et al. (2012) and Weston et al. (2011) Table 1), mainly owing to its higher accuracy and spatial resolution, and thus broader application range, compared to offsettracking and split-bandwidth techniques. Such differences are enhanced for TOPS and ScanSAR-based measurements, compared to the Stripmap case, since these acquisition modes are designed to trade-off in particular azimuth spatial resolution in favor of swath (range) coverage. Of course DInSAR requires a sufficient level of interferometric coherence to be applied, and if this requirement is not met, e.g., due to dense vegetation, earthquakes may be undetected (Funning and Garcia, 2019).

TECHNIQUE APPLICATION
The possibility of complementing DInSAR LoS measurements with azimuth deformation components derived with offsettracking and/or MAI depends on the fault mechanism of the earthquake as well as on the coherence and spatial resolution of the SAR images. Concerning Sentinel-1 IW and ALOS-2 WB modes, azimuth offset-tracking and direct application of MAI by splitting the azimuth bandwidth available within a single burst, restricts the application of these techniques to events associated with large azimuth deformations (>1 m). Successful examples include the 2015 Mw 7.2 Murghab, Tajikistan (Metzger et al., 2017), the 2015 Mw 7.8 Gorkha, Nepal (Grandin et al., 2015;Elliott et al., 2016;Wei et al., 2018) and the 2016 Mw 7.8 Kaikōura, New Zealand events (Hamling et al., 2017;Wang et al., 2018b. BO-MAI provides an increased sensitivity to azimuth motion, and has been applied to the Illapel , Gorkha (Jiang et al., 2017b) and Kumamoto events  concerning Sentinel-1 TOPS, and to the Gorkha  and Kaikōura events (Hamling et al., 2017), concerning ALOS-2 ScanSAR. A limitation in the Sentinel-1 case, is the spatial discontinuity of the azimuth measurements, which can be derived in 1.5 km bands spaced 20 km apart .
For higher resolution Stripmap acquisitions, such as those provided by the ALOS-2 Ultrafine 3 m resolution modes (UBS/UBD), azimuth deformation gradients of 20 to 30 cm can be easily measured in coherent imagery, as demonstrated for the 2014 Mw 6.2 Nagano, Japan event  and for the 2016 Mw 6.0 and Mw 6.2 foreshocks of the Kumamoto, Japan sequence (Kobayashi, 2017). A similar performance can be expected also from Sentinel-1 Stripmap mode data, which, although not routinely acquired, has been successfully used to derive azimuth measurements for the 2014 Mw 6.0 South Napa, USA event (Jo et al., 2017).
Another complementarity among measurement techniques, which has long been exploited (Peltzer et al., 1999), concerns the potential of offset-tracking measurements to provide LoS measurements closer to fault ruptures, where high strain rates often make DInSAR phase unwrapping unfeasible or cause loss of coherence as mentioned in section Differential SAR Interferometry and shown in Figure 3 for the 2016 Mw 7.1 Kumamoto, Japan earthquake. For Sentinel-1 IW and ALOS-2 WB acquisitions offset-tracking was carried out for events associated with large (>1 m) surface displacements, Examples include the 2015 Mw 7.2 Murghab, Tajikistan event (Metzger et al., 2017), the 2015 Mw 7.8 Gorkha, Nepal event (Elliott et al., 2016) and the 2016 Mw 7.8 Kaikōura, New Zealand event (Hamling et al., 2017).
For large-deformation applications/areas where sufficient coherence levels are retained, an even better complement to DInSAR may be represented by split-bandwidth methods, which can provide higher spatial resolutions and accuracy compared to offset tracking, with greatly reduced phase unwrapping requirements (if any) compared to DInSAR. An example concerning the LoS deformation component is given by Jiang et al. (2017a), in which DInSAR, offset-tracking, and range-SBI are applied to Sentinel-1 TOPS data covering the 2016 Mw 7.0 Kumamoto, Japan mainshock and foreshocks. In the azimuth dimension, the effects of ALOS-2 ScanSAR BO-MAI for the aforementioned Kaikōura event (Hamling et al., 2017), based on the approach of Liang and Fielding (2017b), can be compared to the ScanSAR azimuth offset-tracking results of Xu et al. (2018b). It should be noted however, that in cases/areas where interferometric coherence is lost, split-bandwidth methods, like DInSAR, will not yield useful results, whereas offset-tracking of intensity features might still succeed. An example of this is provided by the above-mentioned Gorkha event, comparing the azimuth offset-tracking results of Elliott et al. (2016) with the MAI results of Jiang et al. (2017b), the coverage of which is limited by loss of coherence due to snow in the Himalayas.

CURRENT CHALLENGES
The processing methods described in section Measurement Techniques are implemented in several commercial and freely available software packages, and, as mentioned in section Introduction, there are several existing systems, which carry out DInSAR processing of Sentinel-1 data in an automated fashion. On the other hand only few of the coseismic deformation studies in recent literature (Supplementary Table 1) are based on automated processing chains. For a given earthquake, area of interest and available SAR data, an image analyst must typically take some adhoc decisions. These include what techniques to apply, and assessing whether the results elicit additional processing steps, in particular concerning the mitigation of tropospheric and ionospheric propagation effects and, in the case of, also phase unwrapping errors.

Tropospheric Propagation
Tropospheric propagation affects all LoS measurement techniques, namely DInSAR, offset-tracking (range offsets only) and range-SBI, by introducing a spatially varying differential range delay: which depends on the vertical profiles of partial pressure of dry air, P d (z), Temperature, T(z), and partial pressure of water vapor, e(z), as well as on topographic height, z 0 , and local incidence angle, θ (Hanssen, 2001). In particular (13) can be considered the sum of turbulent component, uncorrelated with topography and characterized by significant horizontal spatial variability, and of a stratified delay component, strongly correlated with topography and slowly varying horizontally (Doin et al., 2009). Tropospheric propagation can be comparable in magnitude to the coseismic deformation and make it difficult or impossible to detect the deformation pattern (Funning and Garcia, 2019). An example is shown in Figure 4A for the 2017 Mw 6.4 Nyingchi, China, earthquake. The effects of tropospheric propagation on azimuth measurements are instead negligible due to the low height of the troposphere compared to the satellite orbit, which implies that different radar pulses backscattered from a same resolution cell traverse virtually the same troposphere (Meyer et al., 2006).
An extensive literature exists concerning both error characterization and mitigation strategies for DInSAR (e.g., Scott and Lohman, 2016;Yu et al., 2018, and references therein), whereas for offset-tracking (and the seldom applied range-SBI) tropospheric propagation effects often lie below the measurement noise floor and are therefore ignored. In recent coseismic deformation studies (Supplementary Table 1) three different approaches are followed to deal with this error source in DInSAR: no action is taken; an empirical model correction is applied; a correction based on a numerical weather model is applied.
Empirical model corrections address topography correlated tropospheric delays. The simplest models are based on the estimation of a linear phase-elevation relation in the farfield, where no deformation is assumed (Wicks et al., 2002), possibly coupled with the estimation of the parameters of a conic to model image-wide error trends, due to atmospheric propagation and/or orbit uncertainties (Cavalié et al., 2007). More general approaches to improve the robustness of stratified delay estimation in the presence of deformation, turbulence, and/or image-wide errors due to orbital uncertainties (Lin et al., 2010;Shirzaei and Bürgmann, 2012;Bekaert et al., 2015). In recent coseismic literature (Supplementary Table 1), linear phase-elevation correction models are applied by Feng et al. (2018), Xu et al. (2018a), Yi et al. (2018), Wang et al. (2017a), Metzger et al. (2017), Xu et al. (2017), Wen et al. (2016), Motagh et al. (2015), Copley et al. (2015). Ganas et al. (2016) analyze the effects of different linear corrections and conclude that they have no impact for that specific deformation study (2015 Mw 6.5 Lefkada, Greece). Non-linear phase-elevation models are applied by Yang et al. (2018a) and Nocquet et al. (2016).
Numerical weather prediction models provide estimates of atmospheric parameters, including those in (13), at a set of pressure levels, several times per day. Correction procedures have been proposed based on global reanalysis data, such as ERA-Interim and MERRA (Doin et al., 2009;Jolivet et al., 2011Jolivet et al., , 2014, on operational analysis models such as HRES-ECMWF , and on meso-scale models (Foster et al., 2006;Puyssegur et al., 2007;Wadge et al., 2010;Kinoshita et al., 2013;Jung et al., 2014). In recent coseismic studies, Jiang et al. (2018) tested the application of ERA-Interim corrections based on the interpolation approach of Jolivet et al. (2011) and of the HRES-ECMWF global model, based on the Iterative Tropospheric Decomposition model of Yu et al. (2017), provided through the Generic Atmospheric Correction Online Service for InSAR (GACOS). The latter was also used by Albano et al. (2018). Feng et al. (2017a) tested the application of ERA-Interim and MERRA global models (Jolivet et al., 2011) as well as empirical corrections (Bekaert et al., 2015). Finally, meso-scale models were used by Kobayashi et al. (2018) and Ozawa et al. (2016).
Although there is clearly not a consensus on the best tropospheric correction approach for a given area of interest, nor on whether corrections should be applied systematically, there is a good agreement in literature concerning the pros and cons of the above-mentioned methods. Empirical approaches may outperform weather models in situations where stratified delays prevail (Kinoshita et al., 2013;Bekaert et al., 2015), but are more limited in dealing with spatial variability of both the turbulent and stratified components, and cannot resolve a potential correlation of deformation with height. Global atmospheric models are readily available, also in near real time if operational models are used (Li et al., 2016b), and have been proved successful in many cases (Jolivet et al., 2014). However, there may be areas and/or conditions under which they do not perform satisfactorily (Walters et al., 2013). Meso-scale models provide the highest spatial resolutions and have been proved very effective (Wadge et al., 2010), but also require expertise in numerical weather modeling to be exploited successfully.
An example of weather-model and empirical corrections is shown in Figure 4. Comparing the correction based on a simple exponential delay vs. elevation model, Figure 4C2, to that based on Iterative Tropospheric Decomposition interpolation of HRES-ECMWF profiles, Figure 4B2, it is clear that the horizontal variability plays an important role for this dataset and in this case it is well-captured by weather model data, provided this more complex interpolation scheme is used to separate the turbulent and stratified delay contributions, rather than performing a more simple bilinear interpolation ( Figure 4D2).

Ionospheric Propagation
Ionospheric propagation affects both LoS and azimuth deformation measurements, regardless of which technique is applied. In the range dimension, a differential phase and group delay with opposite signs are introduced (Meyer et al., 2006): where TEC is the Total Electron Content (TEC) difference between the two acquisitions, f c is the carrier frequency of the radar, and θ the local incidence angle. Azimuth variations of TEC in (22) on spatial scales smaller than the synthetic aperture length (∼1 to 20 km depending on acquisition mode and carrier frequency) also cause a local shift in the temporal azimuth zero-Doppler coordinates given by Gray et al. (2000) and Liang and Fielding (2017b): Such shifts are spatially correlated in the range direction, and appear as "azimuth-streaks" of variable orientation in the azimuth measurements carried out with offset-tracking (13) and MAI (16), (e.g., Gray et al., 2000;Raucoules and De Michele, 2010;Jung et al., 2013). An example of the phase delay (22) and azimuth streak effects (24) is shown in Figure 5, in which L-band ALOS-1 Stripmap data pairs spanning the Mw 6.9 Apr. 13th 2010 Yushu, China, earthquake are analyzed. Figures 5B,C show the DInSAR phase and azimuth offset-tracking results for an image pair in which the coseismic deformation patterns are visible, alongside severe ionospheric effects, which amount to several tens of cm in terms of DInSAR LoS deformation ( Figure 5B) and several meters in terms of apparent azimuth deformation ( Figure 5C). In this case the azimuth streak orientation is actually parallel to the fault. Figures 5D,E show the DInSAR phase and azimuth offset-tracking results for an image pair with less ionospheric propagation effects, in which the coseismic deformation patterns are predominant.
Ionospheric propagation effects received an increased interest in recent years, due to the abundance of data acquired with wide-area acquisitions modes, in particular Sentinel-1 TOPS and ALOS-2 ScanSAR, in which the effects of ionospheric phase delay (22) can be better appreciated, even at C-band where they were previously considered negligible . In recent coseismic deformation studies (Supplementary Table 1) three different approaches are followed to deal with this error source, both concerning the LoS and the azimuth measurements: no action is taken; empirical model corrections are applied; corrections based on the split-spectrum approach are applied (Brcic et al., 2010;Rosen et al., 2010;Gomba et al., 2016. Although not applied in recent literature, other correction methods for the ionospheric phase delay (22) have been proposed, and are reviewed in .
Empirical corrections for the range phase delay (22) are considered only for DInSAR, since such effects can seldom be appreciated by offset-tracking (or range-SBI) measurements, due to the higher measurement noise floors. The parameters of a low-order polynomial are typically estimated in the far field of the coseismic deformation, possibly at the locations of GPS stations, and the derived error model is then removed from the entire interferogram. In recent studies, first or second-order polynomial calibrations are applied to DInSAR measurements by Wei et al. (2018) Galetzka et al. (2015), Kobayashi et al. (2015). DInSAR measurements are instead fitted to GPS data by Morishita et al. (2018) and Kobayashi (2017).
Empirical azimuth offset-tracking corrections to reduce the effects of azimuth streaks (24) consist in directional spatial filtering (Raucoules and De Michele, 2010;Chae et al., 2017) or band-pass filtering (Kobayashi et al., 2009). These approaches are equally applicable to MAI.
The split-spectrum approach conceptually consists in solving Equations (17,18) for τ defo and TEC, using the known frequency dependence of ionospheric delay (22) (Brcic et al., 2010;Rosen et al., 2010). Stripmap and TOPS split-spectrum algorithms are detailed in Gomba et al. (2016, whereas ALOS-2 ScanSAR application is discussed in Liang and Fielding (2017b). In recent literature split-spectrum corrections were applied by Hamling et al. (2017), Yue et al. (2017b), and Liang and Fielding (2017b). In the latter study, BO-MAI ionospheric corrections are also derived based on (24), using the ionospheric phase isolated with the split-spectrum technique.

Phase Unwrapping
Two-dimensional phase unwrapping is a processing step required by DInSAR and, in general, by range and azimuth split-bandwidth methods, as discussed in section Measurement Techniques. Errors occur locally due to phase variations greater than π in magnitude between adjacent pixels. If such variations occur for most of the paths connecting two regions, due to loss of coherence (e.g., due to water bodies, permanent snow cover, vegetation, radar layover/shadow) or discontinuity of the underlying deformation, then a phase jump of one or several multiples of 2π between all pixels in these regions is likely to arise.
In recent coseismic literature (Supplementary error-prone integration paths, e.g., across fault ruptures (Yue et al., 2017b) or high fringe-rate areas (Wang et al., 2018b); manual phase-jump corrections,.between mainland and islands (Moreno et al., 2018) or across coherence gaps Nocquet et al., 2016;Bie et al., 2018;Wei et al., 2018); use of offset-tracking range offsets, either to support manual correction of phase jumps (Xu, 2017) or to reduce ("flatten") the phase-gradients prior to DInSAR phase unwrapping (Baek et al., 2018). An example in which an ad-hoc phase unwrapping approach was used is shown in Figure 3E, in which the presence of water and of an east-west phase discontinuity due to the coseismic deformation pattern would made it impossible to unwrap the phase without violating the assumption that phase variations greater than π radians do not occur between adjacent pixels. Thus, unwrapping was performed separately for the western and eastern sectors of the wrapped interferogram in Figure 2B, with respect to reference points R1 and R2 in Figure 2E. The deformation of the reference points was inferred from the local median of the LoS offset-tracking measurements in Figure 2J.

FUTURE DEVELOPMENTS
The long-term commitment of the Sentinel-1 missions up to at least 2030, together with its free and open data policy, implies that such data will continue to be among the main source for coseismic deformation studies for many years to come. Full exploitation of the properties of its main acquisition mode over land, namely the TOPS IW mode, is unlikely to have been reached by the current state of the art. Concerning the measurement techniques, application of range-SBI and BO-MAI are still rightly considered innovative , and their use can expected to become more systematic in the future. It has also been pointed out in section Differential SAR Interferometry, that the effects of azimuth motion on the DInSAR phase have not been a limitation for recent coseismic studies, although theoretically they could be relevant (González et al., 2015). It could therefore not be excluded that applications will arise, which elicit a better procedure to account for this property of TOPS acquisition, either in the generation or in the subsequent interpretation/modeling of Sentinel-based deformation maps (De Zan et al., 2014).
Concerning the error sources, as reviewed in section Current Challenges, it is not seldom that an image analyst processing a given site and dataset is faced with measurement challenges, which require ad-hoc solutions. Concerning tropospheric propagation effects, new strategies are still being proposed  and still do not appear solve all the observed problems, as discussed in section Tropospheric Propagation. Regarding ionospheric propagation effects, the split-spectrum approach, the application of which to Sentinel-1 and ALOS-2 ScanSAR modes is also recent , could be used systematically to avoid ad-hoc calibration of so called "long-wavelength" errors, and could also be systematically exploited for azimuth-streak corrections of offset-tracking results in coherent areas and MAI results . Concerning phase unwrapping errors, a more systematic comparison with offset-tracking and range-SBI results could provide an additional tool to make decisions less subjective, in particular for Sentinel-1 TOPS IW data, which provides a large range bandwidth and is less sensitive to ionospheric propagation effects compared to L-band systems.
Finally, in a context of standardization of the processing algorithms and increased data availability, in particular in connection with the systematic wide-area products which will be provided openly by the NASA-ISRO SAR (NISAR) mission, the appeal of higher-level products (e.g., interferometric) generated by automated state-of-the-art cloud-based processing systems such as those described in Feng et al. (2016) and Li et al. (2016b) can be expected to increase significantly.

AUTHOR CONTRIBUTIONS
JM carried out the review of processing techniques and applications described in this paper and wrote all parts of the manuscript.