Recent advances in SAR interferometry time series analysis for measuring crustal deformation

Tectonophysics 514-517 (2012) 1–13

Contents lists available at SciVerse ScienceDirect

Tectonophysics journal homepage: www.elsevier.com/locate/tecto

Review Article

Recent advances in SAR interferometry time series analysis for measuring crustal deformation Andrew Hooper ⁎, David Bekaert, Karsten Spaans, Mahmut Arıkan Delft University of Technology, The Netherlands

a r t i c l e

i n f o

Article history: Received 19 March 2011 Received in revised form 11 October 2011 Accepted 14 October 2011 Available online 25 October 2011 Keywords: Time series InSAR PS-InSAR Small baseline InSAR Crustal deformation Eyjafjallajökull Guerrero slow slip

a b s t r a c t Synthetic aperture radar (SAR) interferometry is a technique that permits remote detection of deformation at the Earth's surface, and has been used extensively to measure displacements associated with earthquakes, volcanic activity and many other crustal deformation phenomena. Analysis of a time series of SAR images extends the area where interferometry can be successfully applied, and also allows detection of smaller displacements, through the reduction of error sources. Here, we review recent advances in time series SAR interferometry methods that further improve accuracy. This is particularly important when constraining displacements due to processes with low strain rates, such as interseismic deformation. We include examples of improved algorithms applied to image deformation associated with the 2010 eruption of Eyjafjallajökull volcano in Iceland, slow slip on the Guerrero subduction zone in Mexico, and tectonic deformation in western Anatolia, Turkey. © 2011 Elsevier B.V. All rights reserved.

Contents 1. 2.

Introduction . . . . . . . . . . . . . . . . . . InSAR and time series InSAR . . . . . . . . . . 2.1. Persistent scatterer InSAR . . . . . . . . 2.2. Small baseline InSAR . . . . . . . . . . 3. Recent advances in time series InSAR . . . . . . 3.1. Combined time series InSAR . . . . . . . 3.2. Estimation of tropospheric signal . . . . 3.3. Processing of large data sets . . . . . . . 3.4. Improvement in phase-unwrapping . . . 3.5. Pixels with variable coherence in time . . 3.6. Time series analysis with wide-swath data 4. Conclusions and future directions . . . . . . . . Acknowledgements . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

1. Introduction Since the first interferogram showing differential motion of agricultural fields a little over two decades ago (Gabriel et al., 1989), synthetic aperture radar (SAR) interferometry, or InSAR, has revolutionised the field of crustal deformation studies. Currently there are several satellites ⁎ Corresponding author. Tel.: + 31 15 278 2574. E-mail address: [email protected] (A. Hooper). 0040-1951/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2011.10.013

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. 1 . 2 . 4 . 4 . 5 . 5 . 7 . 8 . 9 . 10 . 11 . 11 . 11 . . 11

acquiring SAR data suitable for InSAR, notably the European Space Agency's Envisat, the Canadian Radarsat-1 and 2, the German TerraSAR-X and the Italian Cosmo-SkyMed satellites, and the number of planned future SAR missions is greater still. However, despite the many successes of InSAR, problems due to changes in the scattering properties of the Earth's surface with time and incidence angle limit the applicability of the technique. Where measurement is possible, the signal due to displacement of the ground is overprinted by signal due to variation in atmospheric properties and errors in both satellite orbit and surface elevation

2

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

determination. This makes detection of slow deformation processes in particular, challenging to recover by standard InSAR techniques. By the end of the 1990s various groups were looking at ways of addressing the limitations of conventional InSAR by processing multiple acquisitions in time. One approach involves identifying “persistent scatterer” pixels, whose scattering characteristics remain stable in time and when viewed from different angles (e.g., Ferretti et al., 2001; Hooper et al., 2004; Kampes, 2005). Another approach involves taking many interferograms formed in the conventional way and inverting them to derive incremental displacements with time (e.g., Lundgren et al., 2001; Berardino et al., 2002; Schmidt and Bürgmann, 2003; Usai, 2003). This is known as the “small baseline” approach. These two broad categories of methods rely on the exploitation of resolution elements on the ground with two different endmembers of scattering properties. The persistent scatterer approach is optimised for resolution cells containing a single point scatterer, whereas small baseline methods are optimised for resolution cells containing a distribution of scatterers. More recently, methods have been developed that take advantage of both types of scattering (Ferretti et al., 2011; Hooper, 2008). These different techniques for InSAR time series analysis have been applied to measure mean velocities and incremental displacements in diverse applications that include tectonic motion associated with faults, volcano deformation associated with movement of magma, landsliding, and subsidence due to the extraction of hydrocarbons and water from the ground.

a

The purpose of this paper is to give an overview of recent advances in InSAR time series analysis for measuring deformation associated with processes that deform the crust. InSAR has applications beyond crustal deformation studies, but we restrict our discussion here to the advances that benefit this field of research. 2. InSAR and time series InSAR Synthetic aperture radar (SAR) is a technique that allows highresolution radar images to be formed from data acquired by sidelooking radar instruments carried by aircraft or spacecraft (Curlander and McDonough, 1991). The amplitude of a SAR image can be interpreted in terms of the scattering properties of the Earth's surface. The phase, on the other hand, is essentially random, as it represents a weighted average of the phase delay between transmitting and receiving for all scatterers on the ground within a resolution element. However, the difference in phase between two images can be interpreted in terms of the change in range from the instrument to the ground, as long as the scattering characteristics of the ground remain approximately the same. Interferometry is the process of multiplying one SAR image by the complex conjugate of a second SAR image resulting in an “interferogram”, the phase of which is the phase difference between the images (Hanssen, 2001; Rosen et al., 2000). An interferogram can be formed between two images acquired at the same time from different positions, but if the aim is to measure

b Line-of-sight Displacement (m) -2 -1 0

39

38

0

50 142

0

50 14

Fig. 1. An example interferogram displaying cosiesmic deformation for the Mw = 9.0 Tohoku-oki earthquake, Japan, which occurred on 11 March 2011. The SAR data were acquired by the L-band ALOS satellite on 28 October 2010 and 15 March 2011. In a the interferometric phase is displayed with each colour cycle representing 11.8 cm of displacement away from the satellite, which was moving in the direction of the white arrow and looking in the direction of the black arrow. In b the phase is integrated and converted to line-of-sight displacement, indicating that the southernmost point moved more than 2.5 m away from the satellite, relative to the northern end. The orbit data provided for 15 March were preliminary and so some orbital signal remains. Signal is partially lost in the mountainous regions due to decorrelation, caused chiefly by snow.

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

deformation of the ground, the images must be acquired at different times. Although InSAR is possible with airborne SARs (Hensley et al., 2009), extensive series of data have yet to be acquired, so time series analysis is only currently possible with spaceborne InSAR. A change in the position of the satellite between the two acquisitions leads to a geometric contribution to the phase change, which can be approximately corrected for knowing the positions of the satellite and the surface topography. The process of correcting for topography is sometimes referred to as differential InSAR (DInSAR). What remains in the interferogram is a contribution due to displacement of the ground between acquisitions, plus some other nuisance terms (Hooper et al., 2007) n o φ ¼ W φdef þ φatm þ Δφorb þ Δφθ þ φN ;

ð1Þ

where φdef is the phase change due to movement of the pixel in the satellite line-of-sight (LOS) direction, φatm is the difference in atmospheric phase delay between passes, Δφorb is the residual phase due to orbit errors, Δφθ is the residual phase due to look angle error (commonly referred to as DEM error, although there is also a contribution from the phase centre sub-pixel position), φN is the phase noise due to both variability in scattering and thermal noise, and W{⋅} is the wrapping operator that drops whole phase cycles, because phase can only be measured in terms of the fractional part of a cycle. What one sees in an interferogram are phase cycles of 2π radians, usually represented by colour “fringes” (Fig. 1a). Phase delay is translated into distance by multiplying by − λ/4π where λ is the wavelength of the SAR system, typically 3.1 cm for X-band, 5.7 cm for Cband and 23.6 cm for L-band systems. In other words, each additional λ/2 of range change results in the identical interferometric phase. As the change in phase delay due to changes in atmospheric properties usually exceeds one cycle, it is not possible to interpret interferometric phase in terms of absolute range change. However, it is possible to estimate the relative range change for two points within an interferogram, by integrating the number of fringes between them (Fig. 1b). This process of estimating the integrated phase is known as phaseunwrapping. The primary limitation of InSAR for most sensors is the phase noise due to changes in the scattering properties. In a SAR image, the amplitude and phase for each pixel comes from the coherent

sum of contributions from all scatterers within the associated ground resolution element (Fig. 2a). Relative movement of these scatterers, or a change in the look or squint angle, causes the scatterer contributions to sum differently, an effect known as decorrelation (Zebker and Villasenor, 1992). The degree of relative movement of the scatterers depends on their size, as larger scatterers tend to be more stable. As more energy is returned from scatterers of about the same size as the wavelength of the radar system, longer wavelength systems, such as L-band, display the least decorrelation. If the decorrelation term is a significant fraction of a phase cycle, the integration of phase difference between points becomes unreliable. This effect can be mitigated somewhat in two ways, each at the cost of resolution. Firstly, by bandpass filtering each image prior to interferogram formation (Gatelli et al., 1994). A change in look or squint angle can be interpreted in terms of a shift in bandwidth, and filtering ensures that only the overlapping frequencies are retained. Secondly, by filtering after interferogram formation (Goldstein and Werner, 1998). The latter can also be achieved by summing the interferometric values of many neighbouring resolution elements, which is known as “multilooking”. As long as the signal does not vary significantly across the area in the multilooked element, the signal in each element reinforces, whereas the decorrelation noise does not. Filtering prior to interferogram formation only reduces decorrelation due to the change in squint angle whereas the second approach also reduces decorrelation due to the relative movement of scatterers. In the case of extreme decorrelation however, such as when the scatterers are inherently non-stationary objects such as leaves, even this fails. The second limitation of InSAR, particularly in the case of small strain, due to processes such as interseismic plate motion and glacio-isostatic adjustment, is that the non-deformation signal can swamp the deformation signal. One approach for extracting small displacements is summing or “stacking” of many conventionally formed interferograms (Simons and Rosen, 2007; Wright et al., 2001; Zebker et al., 1994). The deformation signal reinforces, whereas other signals typically do not. However this approach is only appropriate when the deformation is purely steady-state, with no seasonal deformation, and even then is not optimal, as the non-deformation signals are reduced only by averaging rather than by explicit estimation. Algorithms for time series analysis of SAR data have been developed to better address these two limitations of conventional InSAR. The first limitation is tackled by using phase behaviour in time to

b

Phase

a

Phase

3

0

0 Acquisition

100

0

0 Acquisition

100

Fig. 2. Phase simulations for (a) a distributed scatterer pixel and (b) a persistent scatterer pixel. The cartoons above represent the scatterers contributing to the phase of one pixel in an image and the plots below show simulations of the phase for 100 acquisitions, with the smaller scatterers moving randomly between each iteration. The brighter scatterer in b is three times brighter than the sum of the smaller scatterers.

4

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

Perpendicular Baseline (m)

a

600 400 200 0 200 400 600 800 1000 2003 2004 2005 2006 2007 2008 2009 2010

Perpendicular Baseline (m)

b

600 400 200 0 200 400 600 800 1000 2003 2004 2005 2006 2007 2008 2009 2010

Fig. 3. An example baseline plot for (a) the persistent scatterer method and (b) the small baseline approach. Red circles represent SAR images and blue lines indicate the interferograms that are formed. Perpendicular baseline refers to the component of the satellite separation distance that is perpendicular to the look direction, and is proportional to the difference in look angle.

select pixels where decorrelation noise is minimised. The second limitation is addressed by estimating the non-deformation signal by a combination of modelling and filtering of the time series. The time series algorithms fall into two broad categories, the first being persistent scatterer InSAR, which targets pixels with consistent scattering properties in time and viewing geometry, and the second being the more general small baseline approach.

2.1. Persistent scatterer InSAR Decorrelation is caused by contributions from all scatterers within a resolution cell summing differently, due to relative movement of the scatterers and/or a change in the looking direction of the radar platform. If, however, one scatterer returns significantly more energy than other scatterers within the cell, the decorrelation phase is much reduced (Fig. 2b). This is the principle behind a “persistent scatterer” (PS) pixel, also referred to as a “permanent scatterer”. In urban environments, the dominant scatterers are commonly roofs oriented such that they reflect energy directly backwards, like a mirror, or the result of a double-bounce, where energy is reflected once from the ground, and once from a perpendicular structure (Perissin and Ferretti, 2007). Dominant scatterers can also occur in areas without man-made structures, e.g., appropriately oriented rocks, but there are fewer of them, and they tend to be less dominant. No filtering or multilooking is applied in PS processing as these techniques degrade resolution, thereby adding more scatterers to each resolution element. As nondominant scatterers are considered as noise sources for PS pixels, increasing their number can lead to an increase in decorrelation noise.

PS algorithms operate on a time series of interferograms all formed with respect to a single “master” SAR image (Fig. 3a). It is ultimately the level of decorrelation noise that defines whether pixels are PS pixels or not, but an initial selection of candidate PS pixels can be made using various proxies, the most common of which is amplitude dispersion (Ferretti et al., 2001). There are then essentially two approaches for determining the level of decorrelation noise for each of the candidate pixels. The first relies on modelling the deformation in time (Adam et al., 2003; Crosetto et al., 2003; Ferretti et al., 2001; Kampes, 2005; Lyons and Sandwell, 2003; Werner et al., 2003). Atmospheric and orbit error signal is reduced, initially, by taking the phase difference between nearby candidate pixels. Contributions to this “double-difference” phase from deformation and DEM error are then modelled for the whole time series, with the residuals between the model and the doubledifference phases providing an estimate of the noise level. Network adjustment is applied to determine the noise level for each candidate pixel. The second approach for estimating decorrelation noise relies on the spatial correlation of most of the phase terms (Hooper et al., 2004; van der Kooij et al., 2006). Spatial filtering is applied to estimate the spatially-correlated terms, including the deformation, atmospheric, and orbit error phase, for each PS candidate. The spatiallycorrelated phase is subtracted and the residual contribution from DEM error in the remaining phase is modelled for the whole time series, with the residual between the phase and the model providing an estimate of the noise for the pixel. In the first approach, the phase is unwrapped during the selection process, by fitting a temporal model of evolution to the doubledifference phase, whereas in the second approach a phaseunwrapping algorithm is applied to the selected pixels without assuming a particular model for the temporal evolution (Hooper, 2010). In both approaches, deformation phase is then separated from atmospheric phase, and noise, by filtering in time and space; the assumption is that deformation is correlated in time, atmosphere is correlated in space but not in time, and noise is uncorrelated in space and time. An example of mean velocities estimated from persistent scatterers is shown in Fig. 4. In comparative studies between the two approaches, estimates for the deformation estimates tend to agree quite well, but the second approach tends to result in better coverage, particularly in rural areas (Doin et al., 2010; Sousa et al., 2011). It is also more suited for deformation that is strongly non-linear in time, as no specific temporal model is required (e.g., Fig. 5). The persistent scatterer InSAR technique has the advantage of being able to associate the deformation with a specific scatterer, rather than a resolution cell of dimensions dictated by the radar system, typically on the order of several metres. This allows for very high resolution monitoring of infrastructure. From the point of view of crustal deformation studies, however, this level of detail is usually not required, although it can be useful in separating crustal deformation from the local deformation of specific structures. The accuracy of the technique depends on the sensor, the number of images, the time over which the images are acquired, the distance from the reference point, and the coherence of the PS. In the case where deformation is linear in time, accuracy can be better than 1 mm/year, however (Adam et al., 2009). The precision of individual measurements was assessed in an experiment using corner reflectors to represent ideal PS (Marinkovic et al., 2008), and the relative precision between reflectors 200 m apart was found to be 1.6 mm, for data acquired by the Envisat satellite. 2.2. Small baseline InSAR For resolution elements containing no dominant scatterer, phase variation due to decorrelation is often large enough to obscure the underlying signal. However, by forming interferograms only between

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

5

Fig. 4. Line-of-sight velocities in the San Francisco Bay Area determined from PS analysis. The colour of each point indicates its measured velocity toward or away from the ERS satellite flying toward 193° and looking down from the east at a look angle of ~ 23°. SAF, HF and CF denote the locations of the San Andreas, Hayward and Calaveras faults, respectively. CB, SL, AL, TI and BH show locations of the Cupertino and San Leandro basins, Alameda, Treasure Island and the Berkeley Hills, respectively. Range change rates gradually vary across the region due to elastic strain accumulation about the major plate-bounding faults, but step abruptly across the Hayward fault, which slips aseismically along its surface trace. Large subsidence rates due to settling are observed alongside San Francisco Bay such as on Treasure Island and in Alameda. From Ferretti et al., 2004.

181.1 181

Leveling GPS EDM PS master

Height Difference (m)

images separated by a short time interval and with a small difference in look and squint angle, decorrelation is minimised, and for some resolution elements can be small enough that the underlying signal is still detectable. Decorrelation is further reduced by spectral filtering in range (Gatelli et al., 1994) and discarding of the non-overlapping Doppler frequencies in azimuth. Pixels for which the filtered phase decorrelates little over short time intervals are the targets of small baseline methods. Interferograms are formed between SAR images that are likely to result in low decorrelation noise, in other words, those that minimise the difference in time, look angle and squint angle (Fig. 3b). Obviously it is not possible to minimise each of these simultaneously, so assumptions have to be made about the relative importance, based on the scattering characteristics of the area of interest. In many small baseline algorithms, the interferograms are then multilooked to further decrease decorrelation noise (e.g., Berardino et al., 2002; Schmidt and Bürgmann, 2003). However, there may be isolated coherent resolution ground elements that are surrounded by incoherent elements, such as a clearing in a forest, for which multilooking will increase decorrelation noise. Therefore, other algorithms have been developed that operate at full resolution (Hooper, 2008; Lanari et al., 2004). Pixels are selected based on their estimated spatial coherence

180.9 180.8

1986

1988

1990

1992

1994

1996

1998

2000

in each of the interferograms, using either standard coherence estimation (Li and Goldstein, 1990) or enhanced techniques, in the case of full-resolution algorithms. The phase is then unwrapped either spatially in two dimensions (e.g. Chen and Zebker, 2000), or using the additional dimension of time in 3-D approaches (Hooper, 2010; Pepe and Lanari, 2006). At this point the phase can be inverted to give the phase at each acquisition time with respect to a single image, using least-squares (Schmidt and Bürgmann, 2003), singular value decomposition (Berardino et al., 2002), or minimisation of the L 1-norm (Lauknes et al., 2011). Separation of deformation and atmospheric signals can be achieved by filtering the resulting time series in time and space, as in the PS approach. Alternatively, if an appropriate model for the evolution of deformation in time is known, the different components can be directly estimated from the small baseline interferograms (Biggs et al., 2007). The small baseline method can achieve accuracies similar to the PS technique, on the order of ~ 1 mm/year (Lanari et al., 2007), although this again depends on the number of images, the time over which the images are acquired and the distance from the reference point or area. 3. Recent advances in time series InSAR There have been many advances in time series analysis since the first algorithms were developed. Some of them, however, have limited applicability when the signal of interest is deformation of the crust, rather than deformation of man-made structures, and we do not discuss these further here. Examples in this category are the separation of scatterers within a resolution cell through tomography (Zhu and Bamler, 2010) and the use of polarimetric InSAR for separating different classes of scatterers (Schneider et al., 2006).

Year 3.1. Combined time series InSAR Fig. 5. Relative vertical motion between two benchmarks (23EG and G916) at Long Valley volcanic caldera, as measured by PS analysis in red, and compared to in situ measurements. The error bars represent 68% confidence bounds. The EDM measurements are scaled line length change between benchmarks CASA and KRAK, which is used as a proxy for the vertical deformation. The motion represents inflation and deflation of the resurgent dome in the caldera.

As discussed above, persistent scatterer and small baseline approaches are optimised for different models of surface scattering. The former technique targets resolution cells dominated by a single scatterer and the latter targets cells with many scatterers, none of

6

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

Fig. 6. Comparison of pixels selected by PS and small baseline methods from data acquired by the C-band ERS satellites, on and around Eyjafjallajökull volcano, Iceland. Left, pixels selected by aPS method and middle, pixels selected by a full-resolution small baseline method. The pixels are plotted on topography in shaded relief, with white representing the approximate area of permanent ice cover. The location of the area analysed is shown left inset. 27 images were used in the analysis although only one interferogram is shown here, which spans 27 June 1997 to 10 October 1999, and shows deformation due to the intrusion of a sill at 5.7 ± 0.5 km. Each colour fringe represents 2.8 cm of displacement in the lineof-sight. Right is a comparison of estimated coherence magnitude (γx) for all pixels selected by either, or both, methods. These values are estimated from the residual phase after subtraction of the spatially-correlated phase and correction for look angle (DEM) error (Hooper et al., 2007). A higher coherence magnitude indicates less phase noise. From Hooper, 2008.

which dominate. One might suppose that a full-resolution small baseline approach is equally good as a persistent scatterer approach also for cells dominated by a single scatterer, but in fact there are two advantages that the PS approach has in this case. The first is that all interferograms can be created with respect to a single master. This allows for a reduction of the noise contribution of the master image

a

prior to phase unwrapping, as it is present in all interferograms. The second is that no spatial filtering is applied, which avoids increasing the noise contribution of non-dominant scatterers by coarsening of the resolution. Therefore the two approaches can be considered as complementary, in the usual case where a data set contains pixels with a range of scattering characteristics. By combining both

b –10–5 0 5 1015 displacement (cm)

c Mexico City

0

50

Fig. 7. Topographic phase delay correction for a PS interferogram spanning the 2006 slow slip event on the Guerrero subduction zone. The interferogram was formed from data acquired by the C-band Envisat satellite and includes deformation between 6 December 2005 and 16 March 2007. a shows the unwrapped interferometric phase converted to line-of-sight displacement, with positive values indicating movement towards the satellite, which was moving in the direction of the white arrow and looking in the direction of the black arrow. b shows the differential tropospheric phase delay estimated using the method described in Section 3.2. c shows the result of subtracting the correction in b from the interferogram in a. The Mexico City area is undergoing rapid deformation not connected with the slow slip event and is therefore masked out. The differential tropospheric phase delay is correlated locally with topography, as expected, but a longer wavelength component uncorrelated with long wavelength topography is also captured.

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

approaches signal can be extracted from more pixels overall, improving the spatial sampling, and the signal-to-noise ratio can be improved for pixels that are selected by both approaches. Improvement of the spatial sampling is important not only because the resolution of any deformation signal is increased, but also because it allows for more reliable phase-unwrapping. An algorithm that combines both PS and small baseline approaches was developed by Hooper(2008). PS pixels are selected based on the method of Hooper et al.(2007) and a full resolution small baseline method is applied to extract coherent distributed scattering targets (Fig. 6). 3-D phase unwrapping is then applied to the combined data set (Hooper, 2010) and further processing for isolation of the deformation signal proceeds as for small baseline processing. Although the coherence for the single image shown looks superficially similar for both approaches, the right panel of Fig. 6 shows that the coherence estimated from the entire time series of images actually differs between approaches. The primary reason for these differences is due to the spectral filtering applied in the small baseline processing. In the case of distributed scatterers, the coherence increases (i.e., the phase noise decreases) for each interferogram, due to removal of non-overlapping bandwidth between the two constituent SAR images. In the case of pixels with a strongly dominant scatterer, the coherence decreases due to the noise contributed by the additional scatterers included in the larger resolution element. An alternate algorithm for combining PS and small baseline approaches is presented in Ferretti et al.(2011). In Wegmuller et al. (2010), a hybrid method is developed that is essentially a PS approach but interferograms are formed with small temporal baselines to aid in phase-unwrapping in the case of high deformation rates (>50 cm/year).

a

5

05

b

10

7

3.2. Estimation of tropospheric signal The propagation of the signal through the atmosphere influences the phase delay, principally due to interaction with the ionosphere and the troposphere. This effect varies over an image and is commonly referred to as an “atmospheric phase screen” (APS). Thus interferograms also contain a phase contribution due to the difference in the APS between acquisitions. At most latitudes, significant contribution to the APS from the ionosphere is only present at wavelengths of 100s of km, although shorter wavelength variation can sometimes be observed, particularly when using L-band data. In this section we consider the tropospheric contribution to the APS, which is typically most significant on length scales of a few km to 10s of km (Fig. 7). Standard time series algorithms rely on separating the APS from deformation by filtering in time and space. However, disentangling non-steady deformation from APS is challenging, and even in the case of a steady deformation rate, significant improvements in deformation accuracy can be achieved by reducing the APS before filtering. In the case of flat areas, the tropospheric APS is due only to lateral variation of the line-of-sight integrated water vapour distribution. The strength of this variation is typically limited to a few cm of delay over 10s of km (Hanssen, 2001). However, in areas of relief, there is an additional contribution due to the fact that the higher the ground is, the shorter the fraction of the troposphere that is traversed by the signal. Contribution to the phase delay from this effect can be 10s of cm (Liu et al., 2009), even when considering only the difference between the APS of two acquisitions, present in an interferogram. The APS contribution that correlates with topography can be estimated using weather models (Pinel et al., 2011; Wadge et al., 2002), continuous GPS stations (Onn and Zebker, 2006; Williams et al., 1998), spectrometer measurements, optionally combined with

c

d Mexico City

displacement (cm)

A

B

C

0

50 –10

0 cm

10

–10

0 cm

10

Fig. 8. Estimated deformation due to the 2006 slow slip event on the Guerrero subduction zone. a and c show line-of-sight displacements estimated from unwrapped PS interferograms, before and after correction for topographic phase delay using the method described in Section 3.2. b and d show a comparison between estimated PS displacements (grey dots) along the profile marked in a and c and displacements from nearby GPS stations (black triangles) mapped into the line-of-sight (red error bars). The slow slip displacement and the interseismic displacement rate were simultaneously estimated using weighted least-squares from 18 PS interferograms formed from data acquired by Envisat, five of which span the event. The satellite was travelling in the direction of the white arrow and looking in the direction of the black arrow. The signal in the uncorrected estimate is dominated by spurious long wavelength signals, whereas the deformation in the south dominates after correction.

8

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

five of the interferograms. Fig. 8 shows a comparison between displacement estimates from standard PS processing and when tropospheric APS is first estimated and subtracted. Errors with wavelengths of 10s of km are reduced by the additional processing, and the deformation at the southern end is preserved, such that the agreement with displacements measured at nearby GPS stations is improved. Fig. 9 shows time series for 3 locations before and after the correction. In general, scatter is reduced, leading to lower uncertainties for the estimated secular motion and slow slip displacement. The estimates themselves are also more realistic when compared to the GPS data. When using these data to solve for slip on the subduction interface, the resulting maximum likelihood distribution is significantly different after correction, demonstrating the importance of reducing the tropospheric signal.

weather models or GPS (Li et al., 2009; Puysségur et al., 2007), or the interferometric phase of non-deforming areas (Cavalie et al., 2007; Wicks et al., 2002). All of these methods have had some success, but also have limitations. Global weather models have a coarse resolution (>50 km). Although they can be run in a nested fashion to deliver finer resolution, the assimilated data that constrain them are still coarse, resulting in large uncertainties in their predictions. Continuously acquired GPS data only provide an estimate of tropospheric delay above the station itself, and the density of GPS stations is again coarse for most of the world. Coincident spectrometer measurements are only available for SAR data acquired by Envisat, and fail to provide estimates of water vapour at night or where clouds are present. Finally, when using the interferometric phase itself, the extent of any non-deforming areas may be limited, and only an average function relating interferometric APS to elevation is estimated for the whole image. We have developed a new method that uses interferometric phase to estimate a laterally-varying function relating differential tropospheric phase delay, and can be applied in areas that are deforming (Bekaert, 2011; Bekaert et al., 2010). A power law function is assumed to relate height to interferometric phase, and is estimated for many small windows over each interferogram. The height at which water vapour becomes negligible is estimated from sounding balloon data. To avoid deformation signal biasing the estimate, the interferogram and elevation model are first filtered spatially into several bands, and the function is estimated for every window, for each spatial frequency band. For each window, we retain estimates only for frequency bands that show consistency. These estimates are then combined and extended over the whole interferogram by interpolation. Lin et al.(2010) have also developed a method based on filtering the data spatially, but estimate only a single linear function relating height and interferometric phase over the whole interferogram. We have applied our method to PS time series, but the method is equally applicable in a small baseline approach, in which case network adjustment should be carried out to ensure consistency in the APS contribution of each image to each interferogram (Lin et al., 2010). The result for a single PS interferogram covering the Guerrero region in Mexico is shown in Fig. 7. Because the effects of orbit errors on interferometric phase also partially correlate with topography, it is possible that some of the estimated topographic phase delay is an orbital effect but, in any case, this is still a nuisance signal that we wish to subtract. We applied the same method to all 18 interferograms formed during PS processing, and estimated the displacement due to slow slip that occurred on the subduction interface in 2006, which was present in

3.3. Processing of large data sets The processing of large spatial data sets presents a number of challenges. From a computational perspective it can be exhausting both in terms of computational power, memory and storage. These issues can be addressed in terms of computational strategy. A more significant challenge relates to increasing nuisance contributions in the phase delay. Orbit errors, the influence of the ionosphere and the influence of the troposphere all increase in significance with distance. The first two are typically ignored or estimated in terms of a simple bilinear function for small areas, but this is no longer appropriate for areas spanning 100s of km. The reduction of tropospheric errors is treated separately in Section 3.2. For large data sets it is not typically possible to hold all the data in memory at full resolution during processing. A single full-resolution 100 × 1000 km interferogram acquired by ERS or Envisat is ~9 GB, and the analysis may include more than a hundred such images. One approach is to divide the data into patches spatially, which can then be processed individually. This has the added advantage that each patch can be processed on a separate computer node, allowing for easy parallelisation. As discussed in Section 2, it is optimal to process data at high resolution, but this only applies during pixel selection. Once pixels with useable signal have been selected, it usually makes sense to downsample the signal spatially to a rate that is suitable for sampling the deformation signal of interest. We have implemented an algorithm to downsample the data to a fixed sampling distance, but in principle a variable sampling strategy such as quadtree processing (Jónsson et al., 2002) could be applied. The phase at each sample point is calculated as the weighted sum of nearby pixels, with the weight based on the signal-to-noise ratio of each pixel as estimated from the time series processing. This is akin to smart multilooking; because pixels with zero

Before Correction

After Correction

Line of sight displacement (cm)

20

A

A

B

B

C

C

10

0

-10

-20 2005

2006

2007

2008

2009

2005

2006

2007

2008

2009

Fig. 9. Time series of line-of-sight displacements for positions A, B and C marked in Fig. 8, before and after correction for topographic phase delay. The red lines show the weighted least-squares model for secular motion and displacement during the slow slip event. Before correction a significant slow-slip displacement is estimated at position A, where none is expected, and the displacement at position C is underestimated compared to a nearby GPS site (Fig. 8). The residuals between the data and the model are also reduced for all positions after correction.

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

9

Fig. 10. Line-of-sight displacement rates for a region in western Anatolia, Turkey. (a) The results from PS processing of 25 SAR images acquired by the ERS satellite between 25 January 1993 and 28 December 2000. (b) The modelled long wavelength error, estimated by applying the method described in Section 3.3, using horizontal GPS velocities from Aktug et al. (2009). (c) The remaining signal after subtraction of the long wavelength signal. Arrows indicate GPS velocities and InSAR rates are relative to the most northwesterly measurements. The resulting displacement rates show small variations across some of the graben-bounding faults in the south and large variations across the North Anatolian Fault in the north, due mainly to the Izmit earthquake which occurred in 1999 (see Fig. 11).

signal-to-noise ratio are not included, this results in better results than can be achieved by starting the analysis with low-resolution (multilooked) data. Once the data have been downsampled, individual patches can then be recombined before phase unwrapping. To account for the influence of orbit errors and the ionosphere in large data sets, one approach is to estimate and subtract quadratic functions (Pritchard et al., 2002). Fournier et al.(2011) show that quadratic functions are sufficient to model orbit error effects even over distances greater than 1000 km. The disadvantage with this approach is that some of the long wavelength deformation associated with processes such as interseismic strain and postseismic relaxation may also be subtracted. When available, other deformation data such as GPS measurements can be used to constrain the orbit errors and the ionosphere (e.g., Argus et al., 2005;Brooks et al., 2007; Lundgren et al., 2009; Pritchard et al., 2006). In Arıkan et al.(2010a), a method for correcting the mean displacement rates over hundreds of km using GPS is presented. The InSAR displacement rates are modelled with a 2-D third-order polynomial, the results of which are then subtracted from the displacement rates. GPS horizontal displacement rates are projected onto the local line-of-sight vector, and also modelled with a third-order polynomial, the results of which are then added to the residual InSAR displacement rates. Thus the long

wavelength deformation is constrained by the GPS, and the shorter wavelength deformation is constrained by the InSAR (Figs. 10 and 11). In cases where significant vertical displacement is also expected in the long wavelength deformation signal, e.g., at subduction zones, the GPS vertical rates should also be included in the line-of-sight displacement rate calculation. An alternative approach is to use both InSAR and GPS displacement rates to constrain an Earth model. Two approaches have been developed to achieve this. The first divides the area into triangular elements and solves for the horizontal velocity of each in a weighted least-squares fashion, using Laplacian smoothing as an extra constraint (Wright and Wang, 2010). The second approach solves for the strain tensor at any point using weighted least-squares and assuming an exponential covariance function for the strain (Guglielmino et al., 2011). All of these techniques are not limited to the processing of SAR data from a single track, but can also include multiple tracks in the inversion scheme. When deformation is not steady-state, the GPS data can instead be used to constrain each acquisition. Gourmelen et al.(2010) present a method to do this by solving for the orbit error at each GPS point in each image, and modelling the orbit error as a one-dimensional second-order polynomial.

40

Izmit Earthquake

Line-of-sight displacement (mm)

3.4. Improvement in phase-unwrapping

20 0 –20

A

–40 –60

B 1993

1994

1995

1996

1997

1998

1999

2000

Fig. 11. Time series of line-of-sight displacements across the North Anatolian Fault after correction of the long wavelength signal using GPS. Red crosses and blue circles represent the phase of PS pixels within 1000 m of the positions marked A and B, respectively, in Fig. 10c. The red and blue lines show the best-fitting linear displacement rates prior to, and after, the Izmit earthquake, showing a clear postseismic response.

The original persistent scatterer algorithms unwrap the phase in the time dimension. Because changes in atmospheric delay effectively decorrelate the phase in time, algorithms usually rely on taking the phase difference between nearby PS pixels (double-difference phase) and unwrapping that instead. This is achieved by minimising the phase difference between the double-difference phase and a parametric function, such as a linear phase evolution in time (Ferretti et al., 2001). DEM error is also solved for simultaneously. The unwrapped double-difference phases are then tested for network consistency and integrated spatially. In the case where doubledifference deformation is not easily modelled in parametric form, or the evolution of the deformation is not known a priori, this approach may not be successful. Caro Cuenca et al.(2011) improve the chance

10

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

of success by combining, in a Bayesian fashion, the probability density function (pdf) for each parameter derived from the double-difference phases, with a pdf derived from the estimated parameter values for other pixels. The parameter values that maximise the posterior probability are then selected. Other approaches have been developed that unwrap the phase without the assumption of a particular phase evolution in time. Instead, the algorithms incorporate both the temporal and spatial dimensions in the phase-unwrapping problem. Hooper(2010) uses the phase evolution in the temporal dimension to guide unwrapping in the spatial dimension. Double-difference phases are filtered temporally to give an estimate of unwrapped displacement phase for each acquisition plus an overall estimate of the phase noise. This is used to build a pdf for each unwrapped double-difference phase in each interferogram, and an efficient algorithm is then applied to find the unwrapped phase that maximises the joint probability for each interferogram (Chen and Zebker, 2000). In the case where an area of interest is small enough that the atmospheric signal does not vary significantly, fully three-dimensional algorithms may be applied. These treat the phase difference between points in the temporal dimension in the same way as the doubledifference phases, although with different weighting. A method based on minimisation of the L ∞-norm of the residuals between unwrapped and wrapped double-difference phases was proposed by Hooper and Zebker(2007). More recently Costantini et al.(2010) and Shanker and Zebker(2010) have developed approaches that minimise the L 1-norm, using efficient linear programming algorithms. In most cases minimisation of the L 1-norm is expected to give more reliable results. Their algorithms also have the benefit of being able to incorporate other data as extra constraints, such as GPS data. In the original small baseline algorithms phase-unwrapping is achieved for each small baseline interferogram individually using a

2-D algorithm applicable to sparse data (e.g., Costantini and Rosen, 1999). Because a redundant network of interferograms is formed, small baseline algorithms have the advantage that phase-unwrapping of individual interferograms can be checked for network consistency (Pepe and Lanari, 2006). An alternative approach uses the temporal dimension to guide phase-unwrapping in the spatial dimensions, in the same way as for persistent scatterer time series (Hooper, 2010). 3.5. Pixels with variable coherence in time The time series methods outlined above rely on identifying pixels for which the phase can be tracked in all interferograms. However, there are often pixels which may be coherent in some interferograms but not in others (Colesanti et al., 2004). This can happen when the scattering characteristics change suddenly in time, for example when a field is ploughed or when snow covers the ground. In the former case, interferometric phase for a pixel between two acquisitions prior to ploughing or two acquisitions post ploughing may be coherent, but the phase between images spanning the ploughing event will be incoherent. In the latter case, the interferometric phase formed with one of the snow covered acquisitions will be incoherent, whereas any pair of snow-free images may result in coherent phase. Depending on how the selection parameters are set, these pixels will either not be selected by conventional algorithms, or will be selected, but contribute only noise in certain interferograms. In some cases it may be enough to identify and drop interferograms in which these “partial-PS” (otherwise known as temporary PS and semi-PS) are not coherent (Arıkan et al., 2010b). In other cases, retaining all interferograms may be preferable in order to better constrain time-variable deformation. One approach to identifying partial-PS pixels is to estimate the coherence of every pixel in every possible interferogram that can be

Fig. 12. Partial-PS processing results for persistent scatterer interferograms covering the Eyjafjallajökull eruption, formed from 29 images acquired by TerraSAR-X. The master image was acquired on 3 September, 2009. For a and c the second image was acquired 31 March, 2010 during the initial effusive eruption and the deformation is due mainly to preeruptive emplacement of sills and dikes. For b and d the second image was acquired on 5 June, 2010 and the deformation includes additional deformation due to the reduction of pressure of a magma body beneath the caldera, which occurred during the explosive eruption. Each colour cycle represents 1.55 cm of displacement in the line-of-sight. a and b show the standard PS selection results, whereas c and d show the partial-PS selection results, which are less noisy. A different subset of pixels is selected in each case due to different reasons for the decorrelation. On 31 March, there was snow cover causing decorrelation at higher elevations. On 5 June the snow cover was largely gone (except on the ice cap), but erupted tephra was present over much of the volcano. As well as decorrelated pixels being dropped, partial-PS processing actually leads to an increase in coverage for some areas, as pixels that were dropped by standard PS processing due to their decorrelation in other images are now being selected.

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

Pependicular Baseline [m]

1000

11

with stripmap images (Bertran Ortiz and Zebker, 2007), leading to an enhanced network (Fig. 13). The rest of the processing follows a standard small baseline approach.

500 4. Conclusions and future directions

0

-500

-1000 2003

2004

2005

2006 2007 Time [year]

2008

2009

Fig. 13. An example baseline plot demonstrating the improvement in the small baseline network when utilising both standard image mode and wide-swath data. Red triangles represent standard image mode acquisitions and the magenta crosses represent the wide-swath acquisitions. The black lines show the small baseline interferograms that can be made with only the standard mode images and the dashed blue lines show the extra interferograms that can be made when including the wide-swath images. Modified after Pepe et al., 2011.

formed from the set of images (Ferretti et al., 2011). An alternative method is to identify the partial-PS pixels only in the subset of interferograms that you have chosen to form, which has the advantage of not significantly increasing computing time. We have developed an algorithm to do the latter (Spaans, 2011). This is achieved by first identifying the set of interferograms in which all partial-PS pixels are coherent, using a standard PS processing algorithm (Hooper et al., 2007). For every selected PS pixel, the noise contribution is then estimated in each of the remaining interferograms. This is achieved by calculating the difference between the pixel phase and the spatially correlated phase, after correction for DEM error. If the estimated noise is below a certain threshold, which is chosen such that the noise level of all interferograms is equalized, the pixel is retained for that interferogram. The improvement in the number of pixels selected and the reduction in the decorrelation noise for interferograms spanning two different periods of the 2010 eruption of Eyjafjallajökull in Iceland is demonstrated in Fig. 12. 29 images were used in the analysis, although only two of the PS interferograms are shown in the figure. 3.6. Time series analysis with wide-swath data Some satellites currently in operation have a ScanSAR mode, which enables data acquisition over a wider swath on the ground by illuminating several parallel tracks (Tomiyasu, 1981). Beam elevation steering is used to cycle through the tracks, with each track illuminated by a burst of pulses on each cycle. Envisat in particular has operated in this mode extensively in recent years. Standard time series methods use data acquired in the same mode, but improved temporal sampling can be achieved by combining data acquired in standard image (stripmap) mode with data acquired in ScanSAR mode. Because each track is only illuminated in bursts, with data gaps in between, the resolution of ScanSAR data is degraded in azimuth and combining the two data sets in a PS algorithm does not therefore make sense. However, it is perfectly possible to combine the data sets within a small baseline method (Pepe et al., 2011). A network of small baseline interferograms between stripmap acquisitions from a single track is formed as usual. The same track is extracted from the ScanSAR data. In order to avoid problems due to bursts not exactly overlapping between multiple ScanSAR images, ScanSAR images are only used to form small baseline interferograms

We have presented here some of the recent developments in time series algorithms that enable measurement of crustal deformation to be achieved with improved accuracy. The accuracies achieved in the best case are still on the order of 1 mm/year, but the improvements come in the more common non-best case scenarios. This is helpful for the quantification of all deformation processes, but particularly so for processes that cause small strains, as InSAR is only able to measure relative displacement. Until now, InSAR has principally been used as a tool for analysing deformation processes some time after the fact, but with the upcoming launch of the Sentinel-1 satellites by ESA (planned for 2012 and 2013), SAR data will be acquired for almost every point on Earth at least once every six days, at a similar resolution to that of Envisat. When taken together with data from other SAR satellites, this opens up the possibility of using InSAR for near-realtime monitoring. To achieve this, time series algorithms will need to be adapted to incorporate each new image in an efficient and optimal way, without having to begin the processing from scratch. The noise and other error sources present in interferograms will also continue to be addressed from a technological point of view in future missions. Decorrelation noise is reduced by acquiring images more frequently, increasing bandwidth and using a longer wavelength. Ionospheric phase delay can be estimated using a splitbandwidth system (Brcic et al., 2010; Rosen et al., 2010), as it is frequency dependent, and tropospheric phase delay could be estimated using a system that simultaneously acquires data in a forward and backward looking direction. The accuracy of precise orbits continues to improve, as do elevation models, reducing residual geometric errors. Nevertheless, none of these error sources will be eliminated completely, and time series analysis will therefore always provide an improvement in accuracy, allowing us to detect ever more subtle deformation processes in the future. Acknowledgements We thank Timothy Horscroft for soliciting this review. ERS and Envisat data were provided by the European Space Agency (ESA). TerraSAR-X data were provided by the German Aerospace Centre (DLR). ALOS data were provided by the Japanese Aerospace Exploration Agency (JAXA), and made available through the Geohazards Supersites initiative. Some figures were prepared using the publicdomain GMT software (Wessel and Smith, 1998). We thank Miguel Caro Cuenca and Samiei Esfahany for helpful comments and Matt Pritchard plus an anonymous reviewer for their critical reviews. References Adam, N., Kampes, B., Eineder, M., Worawattanamateekul, J., Kircher, M., 2003. The Development of a Scientific Permanent Scatterer System. ISPRS Hannover Workshop Proceedings. Adam, N., Parizzi, A., Eineder, M., Crosetto, M., 2009. Practical persistent scatterer processing validation in the course of the Terrafirma project. Journal of Applied Geophysics 69 (1), 59–65. Aktug, B., Nocquet, J.M., Cingöz, A., Parsons, B., Erkan, Y., England, P., Lenk, O., Gürdal, M.A., Kilicoglu, A., Akdeniz, H., Tekgül, A., 2009. Deformation of western Turkey from a combination of permanent and cam- paign GPS data: Limits to block-like behavior. J. Geophys. Res. 114 (B13), 10404. Argus, D., Heflin, M., Peltzer, G., Crampé, F., Webb, F., 2005. Interseismic strain accumulation and anthropogenic motion in metropolitan Los Angeles. Journal of Geophysical Research 110, 1–26. Arıkan, M., Hooper, A., Hanssen, R., 2010a. Large-scale tectonics in western Anatolia from time series InSAR. 7th EGU General Assembly 2010, Vienna 12, EGU2010-13909. Arıkan, M., Hooper, A., Hanssen, R., 2010b. Radar time series analysis over West Anatolia. European Space Agency (Special Publication) ESA SP-677.

12

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13

Bekaert, D., January 2011. InSAR Time Series Analysis of the 2006 Slow Slip Event on the Guerrero Subduction Zone, Mexico. Master's Thesis, Delft University of Technology, The Netherlands. Bekaert, D., Hooper, A., Pathier, E., Yun, S., 2010. InSAR Time Series Analysis of the 2006 Slow Slip Event on the Guerrero Subduction Zone, Mexico. Eos Trans. AGU, Fall Meet. Suppl., 91 (52). G12A-03. Berardino, P., Fornaro, G., Lanari, R., Sansosti, E., 2002. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Transactions on Geoscience and Remote Sensing 40 (11), 2375–2383. Bertran Ortiz, A., Zebker, H., 2007. ScanSAR-to-stripmap mode interferometry processing using ENVISAT/ASAR data. IEEE Transactions on Geoscience and Remote Sensing 45 (11), 3468–3480. Biggs, J., Wright, T., Lu, Z., Parsons, B., 2007. Multi-interferogram method for measuring interseismic deformation: Denali fault, Alaska. Geophysical Journal International 170, 1165–1179. Brcic, R., Parizzi, A., Eineder, M., Bamler, R., Meyer, F., 2010. Estimation and compensation of ionospheric delay for SAR interferometry. 2010 IEEE Int. Geosci. and Remote Sensing Symposium (IGARSS), pp. 2908–2911. Brooks, B., Merrifield, M., Foster, J., Werner, C., Gomez, F., Bevis, M., Gill, S., 2007. Space geodetic determination of spatial variability in relative sea level change, Los Angeles basin. Geophysical Research Letters 34 (1), L01611. Caro Cuenca, M., Hooper, A.J., Hanssen, R.F., 2011. A new method for temporal phase unwrapping of persistent scatterer InSAR time series. IEEE Transactions on Geoscience and Remote Sensing 49, 4606–4615. Cavalie, O., Doin, M., Lasserre, C., Briole, P., 2007. Ground motion measurement in the Lake Mead area, Nevada, by differential synthetic aperture radar interferometry time series analysis: probing the lithosphere rheological structure. Journal of Geophysical Research 112, B03403. Chen, C.W., Zebker, H.A., 2000. Network approaches to two-dimensional phase unwrapping: intractability and two new algorithms March Journal of the Optical Society of America. A 17 (3), 401–414. Colesanti, C., Ferretti, A., Perissin, D., Prati, C., Rocca, F., 2004. Evaluating the effect of the observation time on the distribution of SAR permanent scatterers. European Space Agency (Special Publication) ESA SP-550. Costantini, M., Rosen, P.A., 1999. Generalized phase unwrapping approach for sparse data. 1999 IEEE Int. Geosci. and Remote Sensing Symposium (IGARSS), 1, pp. 267–269. Costantini, M., Malvarosa, F., Minati, F., 2010. A general formulation for robust integration of finite differences and phase unwrapping on sparse multidimensional domains. European Space Agency (Special Publication) ESA SP-677. Crosetto, M., Arnaud, A., Duro, J., Biescas, E., Agudo, M., 2003. Deformation monitoring using remotely sensed radar interferometric data. Proc., 11th FIG Symposium on Deformation Measurements, Santorini. Curlander, J.C., McDonough, R.N., 1991. Synthetic Aperture Radar — Systems and Signal Processing. John Wiley & Sons, New York. Doin, M., Lopez-Quiroz, P., Yan, Y., Bascou, P., Pinel, V., 2010. Time series analysis of Mexico City subsidence constrained by radar interferometry. 7th EGU General Assembly 2010, Vienna 12, EGU2010-12031. Ferretti, A., Prati, C., Rocca, F., 2001. Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing 39 (1), 8–20. Ferretti, A., Novali, F., Bürgmann, R., Hilley, G., Prati, C., 2004. InSAR permanent scatterer analysis reveals ups and downs in the San Francisco Bay Area. Eos 85 (34), 317–324. Ferretti, A., Fumagalli, A., Novali, F., Prati, C., Rocca, F., Rucci, A., 2011. A new algorithm for processing interferometric data-stacks: SqueeSAR. IEEE Transactions on Geoscience and Remote Sensing 49 (9), 3460–3470. Fournier, T., Pritchard, M.E., Finnegan, N., 2011. Accounting for atmospheric delays in InSAR data in a search for long-wavelength deformation in South America. IEEE Transactions on Geoscience and Remote Sensing 49 (10), 3856–3867. Gabriel, A.K., Goldstein, R.M., Zebker, H.A., 1989. Mapping small elevation changes over large areas: differential radar interferometry. Journal of Geophysical Research 94 (B7), 9183–9191. Gatelli, F., Guamieri, A.M., Parizzi, F., Pasquali, P., Prati, C., Rocca, F., 1994. The wavenumber shift in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing 32 (4), 855–865. Goldstein, R.M., Werner, C.L., 1998. Radar interferogram filtering for geophysical applications Nov Geophysical Research Letters 25 (21), 4035–4038. Gourmelen, N., Amelung, F., Lanari, R., 2010. Interferometric synthetic aperture radar– GPS integration: interseismic strain accumulation across the Hunter Mountain fault in the eastern California shear zone. Journal of Geophysical Research 115, B09408. Guglielmino, F.G.N., Puglisi, G., Spata, A., 2011. Simultaneous and integrated strain tensor estimation from geodetic and satellite deformation measurements (SISTEM) to obtain three-dimensional displacements maps. IEEE Transactions on Geoscience and Remote Sensing 49 (6), 1815–1826. Hanssen, R.F., 2001. Radar Interferometry Data Interpretation and Error Analysis. Springer. Hensley, S., Zebker, H., Jones, C., Michel, T., Muellerschoen, R., Chapman, B., 2009. First deformation results using the NASA/JPL UAVSAR instrument. 2nd Asian-Pacific Conference on Synthetic Aperture Radar, Oct 2009, pp. 1051–1055. Hooper, A., 2008. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches. Geophysical Research Letters 35, L16302. Hooper, A., 2010. A statistical-cost approach to unwrapping the phase of insar time series. European Space Agency, (Special Publication) ESA SP-677. Hooper, A., Zebker, H., 2007. Phase unwrapping in three dimensions with application to InSAR time series. Journal of the Optical Society of America A 24, 2737–2747.

Hooper, A., Zebker, H., Segall, P., Kampes, B., 2004. A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers. Geophysical Research Letters 31, 23. Hooper, A., Segall, P., Zebker, H., 2007. Persistent scatterer InSAR for crustal deformation analysis, with application to Volcán Alcedo, Galápagos. Journal of Geophysical Research 112, B07407. Jónsson, S., Zebker, H., Segall, P., Amelung, F., 2002. Fault slip distribution of the 1999 Mw 7.1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements. Bulletin of the Seismological Society of America 92 (4), 1377–1389. Kampes, B. M., 2005. Displacement Parameter Estimation Using Permanent Scatterer Interferometry. Ph.D. Thesis, Delft University of Technology. Lanari, R., Mora, O., Manunta, M., Mallorqui, J.J., Berardino, P., Sansosti, E., 2004. A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms. IEEE Transactions on Geoscience and Remote Sensing 42 (7), 1377–1386. Lanari, R., Casu, F., Manzo, M., Zeni, G., Berardino, P., Manunta, M., Pepe, A., 2007. An overview of the small baseline subset algorithm: a DInSAR technique for surface deformation analysis. Pure and Applied Geophysics 164 (4), 637–661. Lauknes, T., Zebker, H., Larsen, Y., 2011. InSAR deformation time series using an L1Norm small-baseline approach Jan IEEE Transactions on Geoscience and Remote Sensing 49 (1), 536–546. Li, F.K., Goldstein, R.M., 1990. Studies of multibaseline spaceborne interferometric synthetic aperture radars Jan IEEE Transactions on Geoscience and Remote Sensing 28 (1), 88–97. Li, Z., Fielding, E.J., Cross, P., Preusker, R., 2009. Advanced InSAR atmospheric correction: MERIS/MODIS combination and stacked water vapour models. International Journal of Remote Sensing 30 (13), 3343–3363. Lin, Y., Simons, M., Hetland, E., Muse, P., 2010. A multiscale approach to estimating topographically correlated propagation delays in radar interferograms. Geochemistry, Geophysics, Geosystems 11 (9), Q09002. Liu, S., Hanssen, R., Mika, A., 2009. On the valuse of high-resolution weather models for atmospheric mitigation in SAR interferometry. 2009 IEEE Int. Geosci and Remote Sensing Symposium (IGARSS), 45 (12), pp. 4043–4058. Lundgren, P., Usai, S., Sansosti, E., Lanari, R., Tesauro, M., Fornaro, G., Berardino, P., 2001. Modeling surface deformation observed with synthetic aperture radar interferometry at Campi Flegrei caldera. Journal of Geophysical Research 106 (B9), 19,355–19,366. Lundgren, P., Hetland, E., Liu, Z., Fielding, E., 2009. Southern San Andreas-San Jacinto fault system slip rates estimated from earthquake cycle models constrained by GPS and interferometric synthetic aperture radar observations. Journal of Geophysical Research 114 (B2), B02403. Lyons, S., Sandwell, D., 2003. Fault creep along the southern San Andreas from interferometric synthetic aperture radar, permanent scatterers, and stacking. Journal of Geophysical Research 108 (B1), 2047–2070. Marinkovic, P., Ketelaar, G., van Leijen, F., Hanssen, R., 2008. InSAR quality control: analysis of five years of corner reflector time series. European Space Agency, (Special Publication) ESA SP-649. Onn, F., Zebker, H.A., 2006. Correction for interferometric synthetic aperture radar atmospheric phase artifacts using time series of zenith wet delay observations from a GPS network. Journal of Geophysical Research 111. Pepe, A., Lanari, R., 2006. On the extension of the minimum cost flow algorithm for phase unwrapping of multitemporal differential SAR interferograms. IEEE Transactions on Geoscience and Remote Sensing 44 (9), 2374–2383. Pepe, A., Bertran Ortiz, A., Lanari, R., Lundgren, P., Rosen, P., Lanari, R., 2011. The Stripmap-ScanSAR SBAS approach to fill gaps in stripmap deformation time series with ScanSAR data. IEEE Trans. Geosci. Remote Sens 49 (12), 4788–4804. Perissin, D., Ferretti, A., 2007. Urban-target recognition by means of repeated spaceborne SAR images. IEEE Transactions on Geoscience and Remote Sensing 45 (12), 4043–4058. Pinel, V., Hooper, A., De la Cruz-Reyna, S., Reyes-Davila, G., Doin, M.P., 2011. The challenging retrieval of displacement field from InSAR data for andesitic stratovolvanoes: case study of Popocatepetl and Colima Volcano, Mexico. Journal of Volcanology and Geothermal Research 200 (1–2), 49–61. Pritchard, M., Simons, M., Rosen, P., Hensley, S., Webb, F., 2002. Co-seismic slip from the 1995 July 30 Mw = 8.1 Antofagasta, Chile, earthquake as constrained by InSAR and GPS observations. Geophysical Journal International 150 (2), 362–376. Pritchard, M.E., Ji, C., Simons, M., 2006. Distribution of slip from 11 Mw > 6 earthquakes in the northern Chile subduction zone. Journal of Geophysical Research 111 (B10), B10302. Puysségur, B., Michel, R., Avouac, J., 2007. Tropospheric phase delay in interferometric synthetic aperture radar estimated from meteorological model and multispectral imagery. Journal of Geophysical Research 112, B05419. Rosen, P., Hensley, S., Joughin, I., Li, F., Madsen, S., Rodriguez, E., Goldstein, R., 2000. Synthetic aperture radar interferometry. Proceedings of the IEEE 88 (3), 333–382. Rosen, P., Hensley, S., Chen, C., 2010. Measurement and mitigation of the ionosphere in L-band interferometric SAR data. 2010 IEEE Radar Conference, pp. 1459–1463. Schmidt, D.A., Bürgmann, R., 2003. Time-dependent land uplift and subsidence in the Santa Clara valley, California, from a large interferometric synthetic aperture radar data set. Journal of Geophysical Research 108 (B9), 2416–2428. Schneider, R., Papathanassiou, K., Hajnsek, I., Moreira, A., 2006. Polarimetric and interferometric characterization of coherent scatterers in urban areas. IEEE Transactions on Geoscience and Remote Sensing 44 (4), 971–984. Shanker, P., Zebker, H., 2010. Edgelist phase unwrapping algorithm for time series InSAR analysis. Journal of the Optical Society of America 27, 605–612. Simons, M., Rosen, P.A., 2007. Interferometric synthetic aperture radar geodesy. In: Schubert, G. (Ed.), Treatise on Geophysics, Vol. 3. Elsevier Press, pp. 391–446.

A. Hooper et al. / Tectonophysics 514-517 (2012) 1–13 Sousa, J.J., Hooper, A.J., Hanssen, R.F., Bastos, L.C., Ruiz, A.M., 2011. Persistent scatterer InSAR: a comparison of methodologies based on a model of temporal deformation vs. spatial correlation selection criteria. Remote Sensing of Environment 115 (10), 2652–2663. Spaans, K., March 2011. Deformation at Katla Volcano, Iceland, 2003–2009: Disentangling Surface Displacements Due to Ice Load Reduction and Magma Movement Using InSAR Time Series Analysis. Master's Thesis, Delft University of Technology, The Netherlands. Tomiyasu, K., 1981. Conceptual performance of a satellite borne, wide swath synthetic aperture radar. IEEE Transactions on Geoscience and Remote Sensing GE-19 (2), 108–116. Usai, S., 2003. A least squares database approach for SAR interferometric data. IEEE Transactions on Geoscience and Remote Sensing 41 (4), 753–760. van der Kooij, M., Hughes, W., Sato, S., Poncos, V., 2006. Coherent target monitoring at high spatial density: examples of validation results. European Space Agency, (Special Publication) ESA SP-610. Wadge, G., Webley, P.W., James, I.N., Bingley, R., Dodson, A., Waugh, S., Veneboer, T., Puglisi, G., Mattia, M., Baker, D., Edwards, S.C., Edwards, S.J., Clarke, P.J., 2002. Atmospheric models, GPS and InSAR measurements of the tropospheric water vapour field over Mount Etna. Geophysical Research Letters 29 (19). Wegmuller, U., Walter, D., Spreckels, V., Werner, C., 2010. Nonuniform ground motion monitoring with TerraSAR-X persistent scatterer interferometry. IEEE Transactions on Geoscience and Remote Sensing 48 (2), 895–904. Werner, C., Wegmüller, U., Strozzi, T., Wiesmann, A., 2003. Interferometric point target analysis for deformation mapping. 2003 IEEE Int. Geosci. and Remote Sensing Symposium (IGARSS).

13

Wessel, P., Smith, W.H.F., 1998. New, improved version of the generic mapping tools released. EOS, Transactions American Geophysical Union 79, 579. Wicks, C., Dzurisin, D., Ingebritsen, S., Thatcher, W., Lu, Z., Iverson, J., 2002. Magmatic activity beneath the quiescent Three Sisters volcanic center, central Oregon Cascade Range, USA. Geophysical Research Letters 29 (7), 1122. http://www.agu. org/pubs/crossref/2002.../2001GL014205.shtml. Williams, S., Bock, Y., Fang, P., 1998. Integrated satellite interferometry: tropospheric noise, GPS estimates and implications for interferometric synthetic aperture radar products. Journal of Geophysical Research 103. Wright, T., Wang, H., 2010. Large-scale crustal velocity field of western Tibet from InSAR and GPS reveals internal deformation of the Tibetan plateau. 7th EGU General Assembly 2010, Vienna 12, EGU2010-7092. Wright, T., Fielding, E., Parsons, B., 2001. Triggered slip: observations of the 17 August 1999 Izmit (Turkey) earthquake using radar interferometry. Geophysical Research Letters 28 (6), 1079–1082. Zebker, H.A., Villasenor, J., 1992. Decorrelation in interferometric radar echoes. IEEE Transactions on Geoscience and Remote Sensing 30 (5), 950–959. Zebker, H.A., Rosen, P.A., Goldstein, R.M., Gabriel, A., Werner, C.L., 1994. On the derivation of coseismic displacement-fields using differential radar interferometry — the Landers earthquake Oct Journal of Geophysical Research 99 (B10), 19,617–19,634. Zhu, X., Bamler, R., 2010. Very high resolution spaceborne SAR Tomography in urban environment. IEEE Transactions on Geoscience and Remote Sensing 48 (12), 4296–4308.