Inverting geodetic time series with a principal componentanalysis-based inversion methodA. P. Kositsky1,2and J.-P. Avouac1Received 13 April 2009; revised 24 August 2009; accepted 16 September 2009; published 3 March 2010.[1] The Global Positioning System (GPS) system now makes it possible to monitordeformation of the Earth’s surface along plate boundaries with unprecedented accuracy. Intheory, the spatiotemporal evolution of slip on the plate boundary at depth, associated witheither seismic or aseismic slip, can be inferred from these measurements through someinversion procedure based on the theory of dislocations in an elastic half-space. Wedescribe and test a principal component analysis-based inversion method (PCAIM), aninversion strategy that relies on principal component analysis of the surface displacementtime series. We prove that the fault slip history can be recovered from the inversion of eachprincipal component. Because PCAIM does not require externally imposed temporalfiltering, it can deal with any kind of time variation of fault slip. We test the approach byapplying the technique to synthetic geodetic time series to show that a complicated sliphistory combining coseismic, postseismic, and nonstationary interseismic slip can beretrieved from this approach. PCAIM produces slip models comparable to those obtainedfrom standard inversion techniques with less computational complexity. We also comparean afterslip model derived from the PCAIM inversion of postseismic displacementsfollowing the 2005 8.6 Nias earthquake with another solution obtained from the extendednetwork inversion filter (ENIF). We introduce several extensions of the algorithm to allowstatistically rigorous integration of multiple data sources (e.g., both GPS andinterferometric synthetic aperture radar time series) over multiple timescales. PCAIM canbe generalized to any linear inversion algorithm.Citation: Kositsky, A. P., and J.-P. Avouac (2010), Inverting geodetic time series with a principal component analysis-basedinversion method, J. Geophys. Res., 115, B03401, doi:10.1029/2009JB006535.1. Introduction[2] Faults slip in a variety of ways, such as during suddenseismic events or as a result of aseismic creep. Fault sliprates can therefore vary over a wide range of timescales,from the typical 10– 100 s duration of large earthquakes, tothe weeks or years duration of slow earthquakes and post-seismic relaxation. Monitoring how fault slip varies with timeis thus key to improving our understanding of fault behavior.Fault slip at depth results in surface deformation that can beobserved with geodet ic techniques [e.g., Lisowski et al.,1991; Segall and Davis, 1997], paleogeodetic techniques[e.g., Taylor et al., 1987; Sieh et al., 1999], or remote sensingtechniques [e.g., Massonnet and Feigl, 1998]. How faults slipat depth can thus be derived indirectly through modeling ofsurface deformation.[3] Theoretical surface displacements expected fromsome fault slip at depth is generally computed based onthe theory of linear elasticity [e.g., Savage, 1983; Okada,1985; Cohen, 1999]. This formulation is linear and easilyinverted using standard algorithms. The distribution of faultslip is generally parameterized based on some discretizationof the fault geometry. The cumulative fault slip needed toexplain displacements that have occurred between twoepochs for which geodetic data are available can then beobtained from some least squares inversion. Because thenumber of parameters generally exceeds the number ofobservations, regularization constraints are generally added;for example, the roughness of the slip distribution can bepenalized or a positivity constraint can be added. One wayto invert geodetic time series for time-dependent slip dis-tribution thus consists in inverting the displacements mea-sured between each two successive epochs. This method iscomputationally very intensive when the number of epochsis large, especially when nonlinear regularization criteria areused. Furthermore, this method considers each epoch indi-vidually, so measurement errors at different time steps arenot properly balanced. In addition, the method also requiresgeodetic time series to be sampled at each site at the sameepochs, limiting the possibility of analyzing a mixed dataset which could include campaign data or interferometricsynthetic aperture radar (InSAR) data.[4] P. Segall and colleagues proposed a variation of theepoch-by-epoch inversion called the extended network inver-sion filter (ENIF) [Segall and Matthews, 1997; McGuire andSegall, 2003] specifically for GPS measurements. ENIF takesJOURNAL OF GEOPHYSICAL RESEARCH, VOL. 115, B03401, doi:10.1029/2009JB006535, 2010ClickHereforFullArticle1Division of Geological and Planetary Sciences, California Institute ofTechnology, Pasadena, California, USA.2Now at Ashima Research, Pasadena, California, USA.Copyright 2010 by the American Geophysical Union.0148-0227/10/2009JB006535$09.00B03401 1of19into account a stochastic description of local benchmarkmotion, a nonparametric description of slip rate as a functionof time, estimation and removal of reference frame errors,and furthermore makes use of an extended Kalman filter tosmooth out noise-related temporal variations. This approachhas been applied with great success in a number of studies[e.g., Miyazaki et al., 2004]. However, the method has somelimitations. One is that it involves a number of hyper-parameters necessary for the model, making it a cumbersometool requiring somewhat subjective choices which define thespace of possible solutions. Another is that the technique iscomputationally costly. The inversion of a fairly modest dataset can take hours or days to complete on a desktop machine.For example, the 400 epoch 10 continuous GPS (cGPS)station data set used to infer the afterslip distribution follow-ing the Nias earthquake [Hsu et al., 2006] took approximately2 h to run on a typical 2005 laptop. ENIF as formulated bySegall and Matthews [1997] is also restricted to the use ofGPS time series, though it could theoretically be extended toallow analysis of any single type of spatiotemporal data (e.g.,InSAR data). Another limitation is that the method is noteasily applicable to the analysis of complex time seriesthat would include interseism ic, coseismic and postseismicdeformation.[5] In this study we describe and then test a principalcomponent analysis-based inversion method (PCAIM)designed to overcome some of the aforementioned