Inversion of Surface Deformation Data for Rapid Estimates of Source Parameters and Uncertainties: A Bayesian Approach

New satellite missions (e.g., the European Space Agency's Sentinel‐1 constellation), advances in data downlinking, and rapid product generation now provide us with the ability to access space‐geodetic data within hours of their acquisition. To truly take advantage of this opportunity, we need to be able to interpret geodetic data in a prompt and robust manner. Here we present a Bayesian approach for the inversion of multiple geodetic data sets that allows a rapid characterization of posterior probability density functions (PDFs) of source model parameters. The inversion algorithm efficiently samples posterior PDFs through a Markov chain Monte Carlo method, incorporating the Metropolis‐Hastings algorithm, with automatic step size selection. We apply our approach to synthetic geodetic data simulating deformation of magmatic origin and demonstrate its ability to retrieve known source parameters. We also apply the inversion algorithm to interferometric synthetic aperture radar data measuring co‐seismic displacements for a thrust‐faulting earthquake (2015 Mw 6.4 Pishan earthquake, China) and retrieve optimal source parameters and associated uncertainties. Given its robustness and rapidity in estimating deformation source parameters and uncertainties, our Bayesian framework is capable of taking advantage of real‐time geodetic measurements. Thus, our approach can be applied to geodetic data to study magmatic, tectonic, and other geophysical processes, especially in rapid‐response operational settings (e.g., volcano observatories). Our algorithm is fully implemented in a MATLAB®‐based software package (Geodetic Bayesian Inversion Software) that we make freely available to the scientific community.


Introduction
Geodetic observational data, most commonly global navigation satellite system (GNSS) and interferometric synthetic aperture radar (InSAR) measurements, are regularly used to infer information about sources of surface displacements and to understand the underlying processes. With these aims, inverse problem theory has been applied to geodetic data to study magmatic systems (Pinel et al., 2014, and references therein), the earthquake cycle (Elliott et al., 2016, and references therein), and many other geophysical phenomena that cause deformation of the Earth's interior and surface such as the response to ice load changes, changes in aquifer storage, and geothermal exploitation (e.g., Auriac et al., 2014;Juncu et al., 2017;Samsonov et al., 2014). However, many commonly employed inversion approaches aim at solely determining an optimal set of source parameters-for example, the weighted least-squares "best-fitting" model-by solving an optimization problem that minimizes the weighted misfit between measured and simulated surface displacements. Among these, the most commonly used methods are simulated annealing (e.g., Cervelli et al., 2001) and genetic algorithm (e.g., Currenti et al., 2005). These have shown to be successful in solving a variety of optimization problems to study different geophysical problems (Sambridge & Mosegaard, 2002), and a detailed analysis of these methodologies, applied to the inversion of InSAR data, is presented by Shirzaei and Walter (2009). Direct-search methods (e.g., simulated annealing and genetic algorithm) do not fully and appropriately characterize uncertainties associated with the source parameter estimates, with the risk that results can be inadequately interpreted. In fact, it is common that a wide range of model parameter values can adequately explain the observations, and it is therefore fundamental to know the credible interval of such values, especially if interpretations are used for the assessment and mitigation of natural hazards or in operational settings (e.g., volcano observatories). account uncertainties in the data (e.g., data errors and incompleteness) and any available prior information (in the form of a prior PDF). The Bayesian method provides the means to investigate a wealth of statistical inferences, such as point estimates (e.g., mean and median of posterior distributions), credible intervals (e.g., quantiles), and direct probability statements about parameters (e.g., the probability that a certain parameter is greater than a certain value). It also allows analyses of joint and conditional probabilities of pairs or sets of parameters and is particularly instructive in the case of non-Gaussian multimodal posterior PDFs. An optimal set of source parameters can also be extracted from the posterior PDF by finding the maximum a posteriori probability solution.
In this work we propose an approach, summarized in Figure 1, for inverting geodetic data-in particular those derived from InSAR measurements-using a Bayesian probabilistic inversion algorithm capable of including multiple independent data sets (e.g., González et al., 2015;Hooper et al., 2013;Sigmundsson et al., 2014). To efficiently sample the posterior PDFs, we implement a Markov chain Monte Carlo method (MCMC), incorporating the Metropolis-Hastings algorithm (e.g., Hastings, 1970;Mosegaard & Tarantola, 1995), with automatic step size selection. We then review and discuss existing methodologies to characterize errors in InSAR data and to subsample large data sets, which are both necessary steps to be performed prior to an inversion. The proposed method is applied to the inversion of synthetic InSAR and GNSS data to demonstrate the ability of the algorithm to retrieve known source parameters. Finally, as a test case, we invert InSAR data spanning the 2015 M w 6.4 Pishan (China) earthquake and determine the fault model parameters for this blind thrusting event and validate our results through the comparison with other independent studies (e.g., Ainscoe et al., 2017;He et al., 2016;Wen et al., 2016).
The proposed approach has been implemented in a software package (Geodetic Bayesian Inversion Software [GBIS], http://comet.nerc.ac.uk/gbis) that we make freely available to the scientific community. The software is written in MATLAB® (which is a commercial software and needs to be installed in order to run GBIS) and uses, among others, analytical forward models from the dMODELS software package (Battaglia et al., 2013). Simple mechanical models of crustal deformation that use closed-form analytical solutions for the characterization of magmatically or tectonically induced deformation processes (e.g., Lisowski, 2007;Segall, 2010) can Figure 1. Schematic representation of the proposed Bayesian inversion approach, including the Markov chain Monte Carlo method-Metropolis-Hastings iterative algorithm. Each step is fully described in the text. b is a random value from a uniform distribution within the range [0, 1]. Note that for i < 20,000, we perform a sensitivity test (see section 2.2) every N s iterations to tune the step size Δm j . GNSS = global navigation satellite system, InSAR = interferometric synthetic aperture radar.

10.1029/2018GC007585
Geochemistry, Geophysics, Geosystems compute surface displacements at 10 3 -10 5 observation points in 10 À3 -10 1 s on consumer-grade computers. These models can be used to place constraints on source location, geometry, orientation, and strength (e.g., volume changes in a magma reservoir and slip on a fault), and their rapidity in computing surface displacements makes them suitable for exploring large numbers (e.g., 10 6 ) of model parameter combinations. Numerical forward models (e.g., boundary elements method) can also be used to explore more complex source geometries or to take into account complexities in the Earth's crust (e.g., nonflat topographic surface; e.g., Bathke et al., 2015;Cayol & Cornet, 1997;Hooper et al., 2011). However, the increase in complexity in the forward model significantly increases the computational burden, making numerical models less efficient in operational or rapid-response investigations.
Finally, with this work we aim at proposing a detailed guideline for the use of geodetic measurements in inverse problems and present a potential standardized procedure that would allow appropriate comparisons of results obtained by different entities. Inversion results are often used to populate global data sets of deformation source models (e.g., Biggs & Pritchard, 2017;Ebmeier et al., 2018) and should therefore be of comparable quality and obtained with a congruent approach.

Bayesian Inversion
For a given discrete inverse problem, the data vector d = {d 1 , d 2 , …, d ND } is equal to a nonlinear model function, G, of the model parameters m = {m 1 , m 2 , …, m NM }, plus error ϵ: In a Bayesian framework, the posterior PDF, p (m|d), describes the probability associated with a given set of model parameters m that is based on how well such parameters can explain the data d given their uncertainties, while considering any prior information. The posterior PDF is calculated as follows: where p(d|m) is the likelihood function of m given d based on residuals between the data and the model prediction of the observations, p(m) expresses the prior information (in the form of a prior joint PDF) of the model parameters, and the denominator is a normalizing constant independent of m.
When the errors are multivariate Gaussian with zero mean and covariance matrix Σ d , ϵ~N(0, Σ d ), the likelihood function is calculated as follows: where N is the total number of data points and Σ À1 d is the inverse of the variance-covariance matrix for a given set of data. The data vector d can be formed from multiple data sets (e.g., multiple SAR interferograms or a combination of interferograms and GNSS data). The likelihood function for multiple data sets, assuming independence, can therefore be expressed as the product of the likelihoods of the data sets (e.g., Fukuda & Johnson, 2008): where N k is the number of data points and r k are the residual vectors (difference between observed and modeled data, r k = d k À G k m) associated with each kth data set and K is the total number of data sets. Similarly, the prior probability of the model vector, assuming independence for the model parameters, is the product of the prior probabilities of the different model parameters m j : where NM is the total number of model parameters.

10.1029/2018GC007585
Geochemistry, Geophysics, Geosystems When inverting geodetic data to infer model parameters for deformation sources, prior information is often not available, in which case a so-called uninformative Jeffreys prior is used (Jeffreys, 1939;Ulrych et al., 2001). This depends on the universe of possible outcomes and, in cases where all values for a given parameter are equally probable, for example, source location parameters, is simply uniform. If there is a finite range of possible values, upper and lower bounds can be added outside which p(m) = 0. However, some source parameters cannot assume negative values (for example, length and width of a rectangular dislocation or the aspect ratio describing a magmatic source) and must be treated accordingly. An appropriate uninformative prior for nonnegative parameters that could equally well be described by the reciprocal of the parameter is a logarithmic prior. In our generalized algorithm, we transform the prior for these cases into a uniform prior, by treating the logarithm of the parameter of interest as the actual model parameter.
Finally, there are other sources of error that are intrinsic to the chosen model (e.g., material properties of the crust) and its ability to accurately predict the data (known as model or prediction errors, e.g., Duputel et al., 2014). While the effect of model/prediction errors should still be considered in any ultimate interpretation (e.g., Jolivet et al., 2015), they are strongly case dependent and difficult to implement in a widely applicable approach. We therefore do not consider model/prediction errors further in this manuscript, as they are also not implemented in the GBIS software.

Sampling the Posterior Probability With MCMC
An efficient way to sample and therefore characterize a posterior probability distribution is through a MCMC sampling method (Mosegaard & Tarantola, 1995). Using the Metropolis-Hastings algorithm (Hastings, 1970;Metropolis et al., 1953), the sampling can be efficiently controlled so that after a sufficiently large number of iterations, the density of the samples approximates the posterior distribution.
The approach, summarized in Figure 1, starts with the selection of an initial set of model parameters m i = 0 from the prior distribution p(m)-either arbitrarily chosen or previously estimated using a direct-search method such as simulated annealing-and the estimation of the associated likelihood function p(d|m i ). A new trial set of model parameters is generated by taking a random step within p(m). When the prior for each model parameter is uniform and independent, this is achieved by perturbing each parameter in m i by an amount a n Δm j , where a n is a random value generated from a uniform distribution within the range [À1, 1] and Δm j is the maximum random walk step size for each parameter m j . If a model parameter of the new model trial falls outside the bounds of the uniform prior probability, we replace this value, m j TRIAL* , as follows: where UB j and LB j are the upper and lower bounds, respectively, of the uniform probability for a given model parameter m j . The likelihood for the new set of model parameters p(d|m i + 1 ) is then calculated, and if its value is greater than the previous one, p(d|m i + 1 ) > p(d|m i ), the step is taken and the trial model values are retained. If the new likelihood value is less than the previous one, p(d|m i + 1 ) < p(d|m i ), the step is taken with a probability equal to the ratio of the new likelihood and the previous one, that is, if p( where b is a random value from a uniform distribution within the range [0, 1]. Otherwise the previous set of model parameters is retained by setting m i + 1 to m i . i is then incremented, and the steps described above are iterated from the selection of a new trial set of model parameters. This approach guarantees that sets of model parameters that improve the probability of the current model are always retained, while those with lower probability are sometimes retained, allowing the algorithm to escape local minima. The process is repeated until a representative sampling of the posterior distribution is achieved (e.g., 10 5 -10 7 iterations).

Automatic Step Size Selection
The efficiency of the Metropolis-Hastings algorithm in sampling the posterior probability distribution depends on the maximum size of the random walk step for each model parameter, Δm j , that can be taken within p(m). If even one step is too small, the inversion will converge too slowly. Conversely, if one or more steps are too large, the algorithm will reject too many proposed sets of model parameters, also leading to slow convergence. To ensure an appropriate convergence and maximize the efficiency of the algorithm, we automatically tune the size of Δm j for each model parameter during the first 20,000 iterations by performing a sensitivity test every N s iterations (e.g., N s = 100 for i < 1,000; N s = 1,000 for 1,000 < i < 20,000). In a similar way to Amey et al. (2018), we aim at finding Δm j such that the mean chance of acceptance approaches the optimal acceptance rate of 23.4% (Roberts et al., 1997) and to ensure that all parameters approximately equally contribute to the changes in likelihood. The optimal acceptance/rejection ratios are maintained by monitoring the evolution of the perturbed probability ratio (PPR) for each parameter, defined as the ratio between the current posterior probability and the posterior probability calculated after perturbing the parameter value by half of the current step size. If the PPR is greater than one, its reciprocal is taken. A target PPR, PPR TARGET , is set as follows: where PPR TARGETi-Ns is the target PPR at the time of the previous sensitivity test, RR CURRENT is the current rejection rate, and RR OPTIMAL is the optimal rejection rate (76.6%). In this approach, if too many model trials are rejected, PPR TARGET will decrease, while if not enough model trials are kept, PPR TARGET will increase.
The step size for each parameter is then adjusted, with the aim that its PPR approaches PPR TARGET . Such changes are achieved by adjusting the maximum step size, Δm j , for each model parameter as follows: first, we calculate the PPR and subtract it from PPR TARGET . We define this difference as p DIFF . If p DIFF < 0, which means that the model parameter's contribution to the posterior probability change is too large, the step size must be reduced. If p DIFF > 0, the model parameter's contribution must be increased by increasing the step size. For each case, the step size is adapted as follows: After 20,000 iterations, the maximum step size, Δm j , is fixed to that tuned through the sensitivity tests for all remaining iterations. The aim is to achieve, at the end of the inversion, an acceptance rate between 15% and 50% (e.g., Roberts & Rosenthal, 2001).

Estimating Convergence and Burn-in Period
After all iterations have been completed (i = I) and before working on the statistics of posterior PDFs, it is important to check that the Markov chain has converged. If a Markov chain does not converge, it does not explore the parameter space sufficiently, and the sampled posterior distribution does not approximate the target distribution well. Note that the convergence for all parameters should be checked. While there are no exact ways to determine the convergence of a Markov chain, there are empirical tools to evaluate it. For example, the visual examination of the trace plots, which display the accepted values of a parameter as a function of iteration number, can be instructive ( Figure 2). Trace plots are also useful to determine the number of early samples that strongly depend on the choice of the starting value, known as the "burn-in" period, during which the MCMC algorithm walks from a highly improbable starting model to models with reasonable values of p(d|m), and starts adequately sampling the posterior PDF. Choosing a burn-in period that is too short may leave unrepresentative samples in the posterior PDF, and since these samples have parameter values that are typically far from those in high posterior probability areas, they may unrealistically increase the credible interval for a given parameter. Conversely, choosing a too large burn-in period will discard representative samples of the posterior PDF that could instead be used to improve the accuracy of the statistics (e.g., Sahlin, 2011). In Figure 2 we display examples of trace plots showing a chain that has converged after the burn-in period, and one that did not converge.
Finally, multiple, independent MCMC inversions can be run in sequence or in parallel, and the results can be combined in a final posterior PDF consisting of the combination of the independent Markov chains (e.g.,

10.1029/2018GC007585
Geochemistry, Geophysics, Geosystems Anderson & Poland, 2016). This step can also be used to check that PDFs for all chains converge to similar distributions, independently from the chosen initial set of model parameters.

Estimates of Data Errors
The Bayesian approach requires that we quantify errors in the data, which we assume are drawn from a multivariate Gaussian distribution. We achieve this by estimating the data variance and covariance, for each independent data set. In our approach, GNSS data are characterized by variances associated to each displacement component, namely, σ 2 N , σ 2 E , and σ 2 V , for the north, east, and vertical components at each site, respectively. Such values, obtained through standard GNSS data processing (e.g., Herring et al., 2015), populate the diagonal of the variance-covariance matrix Σ d , which has all remaining off-diagonal elements set to 0, assuming that no covariance exists between the three components of displacement or between individual GNSS measurement sites. However, when available, covariance between components can be placed in the variancecovariance matrix accordingly.
Conversely, we directly characterize the errors in InSAR data by experimentally estimating variance and covariance in each independent data set. Together with randomly distributed noise caused by phase decorrelation, it is well known that InSAR data are affected by spatially correlated phase delays that are mainly caused by spatiotemporal changes in water vapor content in the atmosphere, known as the "wet" tropospheric delay (e.g., Hanssen et al., 1999). Further spatially correlated errors can be caused by topographic residuals, which are phase change artifacts introduced during data processing. Although errors in InSAR data may sometimes have two-dimensional spatial structures (anisotropy; Knospe & Jónsson, 2010), for simplicity we assume that they are isotropic and stationary (e.g., errors estimated in non-deforming areas are the same as in adjacent deforming areas) and calculate theoretical semivariograms that have the same characteristics as the data (e.g., Lohman & Simons, 2005;Sudhaus & Jónsson, 2009).
A semivariogram measures the spatial variability of a regionalized variable by computing the dissimilarity between pairs of data values (e.g., Wackernagel, 1995). Usually, values at two places near to one another are similar, whereas those at more widely separated distances are less so. By expressing the lag distance h (Euclidean distance) between two observation points x i and x i + h and by averaging increments in equidistant bins, the experimental semivariogram b γ h ð Þ for a quantity Q is calculated as follows: where N is the number of data-point pairs in each distance bin. For second-order stationary processes, the semivariogram and the covariance function are related as follows: where C(0) is the covariance at lag distance h = 0-i.e., the variance-and C (h) is the covariance at any given distance h.
For InSAR data, we recommend calculating the experimental semivariogram over an area of the data set of at least the same size as the area subject to surface deformation under investigation and where deformation is also not expected and/or observed. An efficient way to achieve this is to mask deforming areas from the data set and to calculate the semivariogram using all remaining data points (e.g., Figure 3). It is also recommended to remove a linear ramp of the form Z(x,y) = ax + by, where a and b are linear coefficients of the x and y coordinates, respectively, to account for any residual orbital error or very long wavelength atmospheric delay across the entire InSAR data set (e.g., Sudhaus & Jónsson, 2009). The parameters for a linear ramp across an entire image can then be estimated directly from the data during the inversion.
In our approach, to keep the problem computationally efficient but still adequately representative of the entire data set, the semivariogram is calculated from 3,000 randomly drawn points over 30 equal distance Geochemistry, Geophysics, Geosystems bins (h = maximum distance between points/30). We successively estimate the least squares best fit to the experimental semivariogram using an unbounded exponential one-dimensional function with a nugget: where s 0 is the nugget variance-that is, the variance value as the lag distance approaches 0 and representing the variability at distances smaller than the sample spacing-s is the sill variance (including the nugget variance), and r is the the range (Figure 3). Note that for an unbounded exponential function, where the semivariogram increases asymptotically toward its sill value, the effective range r ε equals 3r and is the distance at which γ(h) has increased to 95% of the sill variance s. Based on equation (10), the covariance can be calculated for any distance between points as follows: A large number of theoretical functions can be used to fit the experimental semivariogram and to simulate the covariance in InSAR data (for example, see González & Fernández, 2011, and references therein). Among these models, we chose the unbounded exponential one-dimesional function with nugget because of its simple implementation and ability to well approximate the power law behavior exhibited by the spatial structure of atmospheric delay.

Subsampling of InSAR Data
InSAR data can provide surface displacement measurements over continuous large areas (for example, a Sentinel-1 wide-swath image product covers an area of~250 × 180 km 2 ) with over 10 8 measurement points when data are processed at full resolution. The inversion of such large data sets would be therefore extremely computationally expensive, making spatial subsampling a necessary step to achieve a tractable computational burden. Any subsampling approach must, on the other hand, maintain enough information as possible for the inversion to be successful.
Three main types of InSAR data downsampling methods can be applied: uniform sampling (e.g., Pritchard et al., 2002), resolution-based sampling (e.g., Lohman & Simons, 2005;Wang et al., 2014), and gradient-based sampling (e.g., Jonsson, 2002;Simons, 2002). Uniform sampling extracts data at regularly spaced locations across the image, and while it is the simplest approach, it is also the least suitable for an effective sampling of those areas where surface displacements are most diagnostic for estimates of deformation source parameters. To achieve the appropriate data point density in such areas, this approach often leads to a data vector that is still too large for efficient inversions. Conversely, resolution-based sampling uses the design of the inverse problem to determine the optimal data sampling density and is based on the data resolution matrix, with higher density of points in the near field and fewer points in the far field. This approach, despite being shown to be the best performing in the inversion of InSAR data (Wang et al., 2014), requires some knowledge of the source of deformation itself (e.g., the fault plane geometry) and is therefore not always applicable, especially when the estimation of source location and geometry is the main goal of the inversion. Finally, gradient-based sampling methods increase data point density in those areas characterized by higher displacement gradients and vice versa. While this approach may oversample areas where higher displacement gradients are not caused by deformation but by noise (e.g., atmospheric delay or topographic errors), it is the most generally suitable in cases where no prior knowledge on the source geometry and location is available. An appropriate characterization of the data errors (see section 3.1) can also mitigate the effect of oversampling noisy areas.
In our approach, we adopt an adaptive quadtree sampling (gradient-based method; Decriem et al., 2010) capable of dealing with points irrespective of whether they are regularly spaced (e.g., data in gridded format) or not (e.g., persistent scatterers). The algorithm recursively divides the data set into sets of four polygons (squares for regularly spaced data) until the phase variance of the points within a polygon is below a given threshold (Figures 4h and 4i). Areas with high variance (e.g., deforming areas, near field) will be subdivided more finely than areas with low variance (e.g., nondeforming areas, far field). Once the points within a polygon have variance lower than the selected threshold, the mean value of such points is assigned to a sample point with the coordinates of the centroid of the polygon. This strategy allows for nonsquare, irregularly aligned polygons, which reduce the effect of data gaps within and at the edges of InSAR images, a shortcoming of the "classic" quadtree approach. Small polygons with less than three points are eliminated to avoid sampling in areas with extreme deformation gradients, where these are likely to be inaccurately unwrapped during processing (e.g., very near field, surface rupture of faults). The variance Figure 4. Synthetic global navigation satellite system (GNSS) and interferometric synthetic aperture radar (InSAR) data sets used to assess the inversion approach. In panel (a) the horizontal components of the GNSS data are plotted with black arrows. Panels (a) and (d) show the InSAR simulated surface displacements and the linear ramp for the descending and ascending data sets, respectively. InSAR line-of-sight displacements are wrapped to 2.8-cm color cycles. Panels (b) and (e) show the simulated atmospheric noise and the coherence mask in white. Panels (c) and (f) are the final data sets. In panel (g) we show the experimental (diamonds) and theoretical (solid line) semivariograms for the two InSAR data sets. Panels (h) and (i) show the subsampled unwrapped InSAR data sets obtained using the quadtree algorithm.

10.1029/2018GC007585
Geochemistry, Geophysics, Geosystems threshold can be iteratively adjusted until the minimum possible data vector length is achieved while maintaining a sampling sufficiently high to characterize the deformation field.
After subsampling the InSAR data, the theoretical covariance function estimated using the experimental semivariogram (see section 3.1) is used to populate the variance-covariance matrix according to the distance between the polygon centroids.

Validation Using Synthetic GNSS and InSAR Data
We demonstrate and assess the validity of our Bayesian inversion approach using synthetic deformation data that simulate surface displacements caused by magmatic sources. In this example, we produce two InSAR line-of-sight (LOS) displacement maps, simulating images from a descending and from an ascending pass of the satellite (Figures 4a-4f), and GNSS three-dimensional displacements at 20 sites randomly distributed over the area of interest (black arrows in Figure 4a show the GNSS horizontal components). Surface displacements are generated using a set of arbitrarily chosen model parameters (Table 1) and forward models for a deflating finite spherical source (McTigue, 1987) and a rectangular dipping dike with uniform opening (Okada, 1985), assuming an isotropic elastic half-space with Poisson's ratio ν = 0.25.
To make the synthetic data more representative of real measurements, we randomly assign a standard deviation to each component of the GNSS data ranging between 1 and 9 mm (the standard deviation of the vertical component is on average five times larger than that of the horizontal components). In the case of InSAR data, we first add a linear ramp and a rigid offset in the LOS direction, we then simulate atmospheric noise using an isotropic two-dimensional fractal surface with a power law behavior (Hanssen, 2001), and finally we remove randomly distributed portions of the data sets to simulate low-coherence areas (Figures 4b  and 4e). Synthetic three-dimensional displacements are projected into the LOS direction using uniform incidence angles of 41°and 29°and heading angles of 191°and 349°for the descending and ascending data sets, respectively.
We calculate the experimental semivariogram of the added atmospheric noise in each InSAR data set and estimate best-fit exponential functions (Figure 4g) with sill variance s DESC = 4.85 mm 2 (s 0DESC = 4 × 10 À9 mm 2 ) and a range distance r DESC = 9,380 m for the descending data, and s ASC = 5.42 mm 2 (s 0ASC = 8 × 10 À9 mm 2 ) and r ASC = 8,350 m for the ascending data. We then apply the quadtree subsampling algorithm and retrieve the data vectors d DESC with 1,453 data points and d ASC with 1,292 data points (Figures 4h-4i).
The inversion algorithm is tested using five data combinations: GNSS data only, the descending InSAR data set, the ascending InSAR data set, both InSAR data sets, and finally all three data sets (2× InSAR + GNSS). We assume a uniform prior probability for all source parameters (logarithmic for nonnegative parameters) between reasonable bounds and perform 10 6 iterations to sample the posterior PDFs (the first 2 × 10 4 iterations are not retained as they represent the burn-in period/step size adjustment). All inversion results are reported in Table 1.
In all cases, the true input model parameters fall between the 2.5 and 97.5 percentiles of the posterior PDFs, which validates our inversion algorithm. As expected, the uncertainty in the model parameters varies as a function of the size of the data vector (e.g., looser when using only GNSS data and tighter for all combinations that use InSAR data) and the errors in the data. The duration of each inversion, when performed on a consumer-grade computer, varies between 3 min (GNSS only) and 3.5 hr (all data).

Application to the 2015 M w 6.4 Pishan Earthquake
As a test case, we apply the Bayesian inversion approach to three InSAR data sets spanning the 3 July 2015 M w 6.4 Pishan (China) earthquake ( Figure 5) Rectangular dike with uniform opening (Okada, 1985) Length ( Note. For each combination of data used in the inversion, we report maximum a posteriori probability (MAP) solutions, 2.5 and 97.5 percentiles of posterior PDFs of model parameters (in squared brackets). GNSS = global navigation satellite system, InSAR = interferometric synthetic aperture radar, LOS = line of sight.
a Duration: Time taken by inversion on an Apple MacBook Pro, 2.7 GHz Intel Core i5 processor, 8 GB 1,867 MHz DDR3 memory.
The three Sentinel-1 interferograms show coseismic deformation consistent with thrust faulting on an ENE-WSW striking structure. Ascending interferograms (Figures 5a and 5d) show two lobes with a characteristic "butterfly" shape, where LOS displacements are toward the satellite (a combination of uplift and predominantly northward horizontal motion) northwest of the epicenter and away from the satellite (a combination of subsidence and predominantly southward horizontal motion) west-southwest of the epicenter. The descending interferogram (Figure 5g), which only covers the northern portion of the coseismic displacement field, also shows LOS displacements toward the satellite north of the epicenter. Maximum LOS displacements of~0.12 m are recorded at the center of the uplifting lobe.
We first estimate the experimental semivariogram for deformation-free regions of each InSAR data set and fit the theoretical exponential function (see Table 2). We then subsample the data sets using a variance  threshold of 3 × 10 À5 m 2 and obtain a data vector composed of 1,146 data points (the contribution of each data set is reported in Table 2). The inversion is carried out using a kinematic forward model for a rectangular dislocation source (Okada, 1985) with nine source model parameters: length, width, depth of the lower edge, dip angle (positive upward from horizontal), strike (measured clockwise from north), X and Y coordinates of the midpoint of the lower edge, uniform slip in the strike direction, and uniform slip in the dip direction. We also estimate a rigid shift in the LOS direction since all InSAR measurements are relative to an arbitrary reference point, and the coefficients for a linear ramp, a (x) and b (y), for each interferogram. The model parameter vector m is therefore composed of 18 parameters, for which we assume uniform priors between reasonable bounds (logarithmic for nonnegative parameters; see Table 3). Priors for fault geometry are based on previous results (Ainscoe et al., 2017;He et al., 2016;Wen et al., 2016) and on the size of the measured deformation field, while no constraints are imposed onto fault strike (varying between 0°and 360°) and dip angles (varying between À90°and +90°from horizontal) and direction of slip, to explore all possible fault orientations and kinematics.
In Figure 6 we show the posterior PDFs for the nine fault source parameters obtained after 10 6 iterations (a burn-in period of 2 × 10 4 iterations is removed), with the bottom row showing histograms of marginal distributions for each parameter and the remaining rows showing the joint distributions between pairs of parameters. The maximum a posteriori probability solution and the 95% credible intervals are reported in Table 3.
The inversion reveals that coseismic surface displacements can be well explained by slip on a 19.5-22.2-kmlong and 7.6-11.7-km-wide fault, striking 113°-120°and shallowly dipping at 20°-29°. The slip direction is consistent with that of a thrust fault with a minor left-lateral component, with 0.4-0.6 m and 0.03-0.14 m of slip in the two directions, respectively. From the analysis of the marginal posterior probabilities (Figure 6), we can observe correlations between fault width, dip, and horizontal position. Conversely, we can well constrain the fault length, depth, strike direction, and the amount of slip in both the dip and strike directions. Assuming a shear modulus of 3.32 × 10 10 N/m 2 (Wen et al., 2016) and the marginal posterior probabilities for fault geometry and slip, the geodetic moment magnitude M w is equal to 6.3 for the entire probability distribution, similar to the U.S. Geological Survey-National Earthquake Information Center and global centroid moment tensor solutions, and in agreement with previous studies of this earthquake (Ainscoe et al., 2017;He et al., 2016;Wen et al., 2016).
In Table 4 we compare our inversion results with previous studies and show a consistent agreement in optimal source parameters. He et al. (2016) performed the source parameter optimization using a simplex algorithm but did not report how the associated uncertainties were estimated. The extremely narrow uncertainties they provide may represent an improper handling of errors in the data or an inappropriate approach in determining uncertainties. Wen et al. (2016) estimated optimal source parameters using a particle swarm optimization and Ainscoe et al. (2017) using a downhill Powell's algorithm. Both these last two studies estimated model parameter uncertainties by running 100 and 250 inversions, respectively, of the InSAR data sets perturbed with simulated noise (e.g., Wright et al., 2003). Uncertainties estimated by Wen et al. (2016) seem still too optimistic while those by Ainscoe et al. (2017) are similar to our estimates obtained using the proposed Bayesian approach.  (Okada, 1985) Length ( Note. We report maximum a posteriori probability solutions, 2.5 and 97.5 percentiles of posterior probability density functions of fault parameters. a Strike slip is positive if right lateral and negative if left lateral. b Dip slip is positive for thrust faulting and negative for normal faulting.

Discussion
Our Bayesian inversion framework is aimed at taking advantage of new opportunities offered by unprecedented improvements in the availability and resolution, both temporal and spatial, of geodetic data. We show that with our approach, the inversion of net surface displacements from GNSS data, which in many cases are telemetered in real-time to geophysical monitoring facilities (e.g., https://www.unavco.org/data/  Ainscoe et al. (2017). a Total slip is calculated from the two components. Uncertainties are calculated through propagation of errors associated with the two slip components. gps-gnss/real-time/real-time.html), can offer robust estimates of source parameters and uncertainties very rapidly (e.g., 3 min in the case of our synthetic GNSS data set). Similarly, we independently or jointly invert InSAR data to provide reliable estimates of deformation source parameters in 1.5-3.5 hr, depending on the size of the data sets. This becomes of unique value when using InSAR data from new satellite missions, such as Sentinel-1, which are acquired with much shorter revisit time (e.g., Sentinel-1, 6-12 days) and made available within a few hours of acquisition.
For example, in March 2017 we applied our approach to an episode of volcanic unrest at Cerro Azul volcano (Galápagos Islands, Ecuador). Using the GBIS software, which implements our Bayesian inversion algorithm, we inverted InSAR data spanning a period of sudden and frequent seismicity beneath the flanks of the volcano. Within 5 hr from being alerted by Ecuador's Instituto Geofisico-Escuela Politecnica National (Geophysical Institute-National Polytechnic School), we were able to access and process Sentinel-1 SAR data for the area of interest, which were acquired only a few hours earlier and which showed contemporary subsidence of the summit of the volcano and uplift of the lower flank (e.g., Figure 7). The unwrapped interferograms were then inverted to gain information on the sources of surface displacements that were accompanying this episode of volcanic unrest, and inversion results were reported to Ecuador's Instituto Geofisico-Escuela Politecnica National in less than 10 hr from receiving the first alert (http://www.igepn. edu.ec/servicios/noticias/1473-informe-especial-cerro-azul-no-2-2017). With this event, we demonstrated the potential of our approach for a quantitative interpretation of InSAR data in a rapid-response operational setting. In fact, our estimates of the credible interval of magma volumes involved in the volcanic unrest episode were used by Ecuador's Instituto Geofisico-Escuela Politecnica National to evaluate hazards associated with the event.
When inverting InSAR data, a frequently applied alternative approach for the estimate of uncertainties in model parameters, also used in previous studies of the 2015 Pishan earthquake (see section 5; Ainscoe et al., 2017;Wen et al., 2016), implies the use of a nonlinear inversion technique capable of finding the weighted least-squares best-fit solution given the data plus simulated noise (e.g., Wright et al., 2003). This process is repeated several times (e.g., 10 2 ) using different simulations of the noise in the data, with the resulting distribution of solutions representing the posterior probability distribution of model parameters.
While this approach and ours should provide similar results in cases were all priors are uniform, Hooper et al. (2010) showed that an MCMC approach, similar to the one presented here, is two orders of magnitude more efficient in terms of the number of forward models that need to be run, hence significantly faster. Therefore, the Bayesian inversion approach may be preferable in rapid-response or operational settings (e.g., volcano observatories) and in all those cases where fast and robust estimates of source parameters may be needed. If prior probabilities on model parameters are also available, then the Bayesian approach is the only one capable of including these in the inversion. The rapidity of the inversion approach in estimating source parameters, especially in the case of a limited number of GNSS sites (e.g., GNSS data only in the validation example), can also be of use in the planning and design of geodetic monitoring networks. Through the use of synthetic data sets, the effect of a given measurement site can be quantified in terms of its contribution to changes in source parameter uncertainties.
The proposed approach is aimed at the characterization of deformation source parameters through the inversion of static surface displacements spanning a given time interval. While this is of great value for early warning and in the rapid response to events of volcanic unrest and earthquakes, it is not optimized for the study of time-dependent dynamic processes. However, inversion results obtained through our approach are complementary and valuable in instructing other time-variable data assimilation algorithms (e.g., the ensemble Kalman Filter approach for volcano monitoring; Bato et al., 2017;Gregg & Pettijohn, 2016;Zhan & Gregg, 2017). Similarly, the Bayesian approach presented here can instruct or be extended to physical and dynamic models of geodetic and other geophysical measurements, as successfully demonstrated by Anderson and Segall (2013) in the study and forecasting of an episode of volcanic unrest.
Our Bayesian inversion framework aims at being applicable to different geophysical processes, in particular those related to tectonic and magmatic activity, and to both scientific research and natural hazard monitoring. To maintain the flexibility of the algorithm and its ability to efficiently invert different types of geodetic data for a multiplicity of deformation sources, we must rely on certain assumptions and reduce the level of complexity of certain steps. For example, a variety of one-and two-dimensional covariance functions could be applied to characterize errors in InSAR data (e.g., González & Fernández, 2011;Knospe & Jónsson, 2010), or different subsampling methods could be applied to subsample the data (see Section 3.2). While, at this stage, such complexities are not implemented in our approach and in the GBIS software, users can optionally adapt the algorithms to better fit their aims and the desired level of complexity.

Conclusions
We have presented a Bayesian approach for the simultaneous inversion of independent geodetic data sets, in particular those from GNSS and InSAR, which takes into account errors in the data and prior information on model parameters. The inversion algorithm, which we implemented in the freely available MATLAB®-based GBIS software, is designed to rapidly estimate optimal model parameters and associated uncertainties through an efficient sampling of the posterior PDFs. Such sampling is performed using a MCMC method incorporating the Metropolis-Hastings algorithm and with an automatic step-size selection. We have applied the inversion method to synthetic GNSS and InSAR data sets and demonstrated its ability to retrieve the true input model parameters. We have also applied the same approach to InSAR data spanning a thrust earthquake and retrieved source parameters for a rectangular fault with uniform slip. Our results are similar to those of previous studies that estimated uncertainties using the added simulated-noise method, but our methodology has been shown to be significantly faster in characterizing optimal model parameters and associated uncertainties, demonstrating its value in rapid-response and operational settings.