Integration of GPS and InSAR Data for Resolving 3-Dimensional Crustal Deformation

,


Introduction
The Global Positioning System (GPS) and the Interferometric Synthetic Aperture Radar (InSAR) are two satellite geodesy methods that have been widely used in recent years to measure crustal deformation.The GPS method can be used to precisely measure 3-dimensional positions and displacements at discrete locations, with up to one millimeter accuracy in horizontal directions and several millimeters accuracy in vertical direction (Bock and Melgar, 2016).The InSAR techniques can be used to measure areal displacements in the direction of radar line-of-sight (LOS) up to several millimeters to centimeter accuracy (Gens and Van Genderen, 1996).These two methods are therefore complementary to each other for crustal deformation monitoring, and efforts have been made to combine these two kinds of observations with common spatial and temporal span, for better spatial and temporal resolution than using either one of them.Such kinds of efforts include: 1) Construct a 3-dimential velocity field using a GPS derived velocity model to control the long-wavelength deformation and InSAR data to constrain the shortwavelength deformation (e.g.Tong et al., 2013).2) Interpolate 3-dimensional GPS velocity and combine that with the InSAR LOS rate data by point-by-point least-squares regression (e.g.Samsonov et al., 2007Samsonov et al., , 2008)).3) Integrate 3-dimensional GPS time series at discrete locations with 1-dimensional InSAR LOS time series data for 3-dimensional continuous time series (e.g.Gudmundsson et al., 2002).In this study we focus on the approach 2), and develop an algorithm to optimally integrate GPS and InSAR data sets for the production of 3-dimensional crustal velocity solution.We will also demonstrate the usefulness of the algorithm with a case study at a selected region in southern California.This method can be extended further to the combination of GPS and InSAR time series data.The code to perform the combination is released to interested users as a supporting information dataset to this paper.

Methodology
GPS station velocities can be derived from position time series of either campaign or continuous GPS observations.In this study as an example, we use velocity solutions of continuous GPS (CGPS) sites produced by the MEaSUREs project (ftp://sopacftp.ucsd.edu/pub/timeseries/measures/ats/),and of campaign GPS sites from the SCEC Crustal Motion Map version 4 (CMM4) solution (Shen et al., 2011) (Figs. 1 and 2).The CMM4 velocities are rotated to align with the CGPS solution which is referenced to the stable North America reference frame (SNARF) (Herring et al., 2008).We divide the GPS data into two groups.The first group utilizes the 3-dimential velocity components for solution, which includes most of the CGPS sites from the MEaSUREs project.The second group utilizes only the horizontal velocity components, which includes the CMM4 sites and a small portion of the CGPS sites whose vertical velocities show anomalously large and possibly non-tectonic signals.
Both data sets are screened to remove outliers, and 1052 horizontal and 542 vertical site velocities are employed.Separate interpolations are performed for the horizontal and vertical velocity fields, to account for different data populations.
In our algorithm of GPS and InSAR data integration, point-based discrete GPS velocities are first interpolated to produce continuous 3-D vector map at chosen grids.The interpolation is based on an algorithm of Shen et al. (2015), which takes into account GPS network density and configuration for data weighting.A Gaussian distance weighting function (w d ) and a Voronoi At each grid point σ k is adjusted to meet W, which is a predetermined constant.With this adjustment, less smoothing is performed and better resolution is achieved for grids with denser data coverage, and vice versa.This approach can also effectively smooth out the incoherencies in discretized GPS velocity data and make robust joint inversion result.Selection of the parameter W allows an overall control of the degree of smoothing for the solution.Greater W would result in more sites included for interpolation and more smoothed solution with less resolution, and smaller W would result in less sites included and less smoothed solution with better resolution.
An optimal balance can be achieved by assessing the overall data strength of the project.
To combine the interpolated GPS data with InSAR data, we need adequate estimates of GPS velocity uncertainties, to be used as data weighting in the combination.Formal GPS velocity uncertainties deduced in the interpolation process, however, are not fit for the job because they are largely determined by the amount of a priori information (i.e. the degree of smoothing) imposed during interpolation, which varies from grid to grid.It usually leads to apparently unreasonable results, that regions with sparser data points would have smaller uncertainty than regions with denser data points, and vice versa.To overcome the problem, we propose to propagate errors from GPS data input to interpolation output using the same interpolation functional form and least-square procedure as before, but not to alter the smoothing distance parameter σ k .Instead, σ k will be kept as a constant σ 0 for all the region, usually chosen as an average of previous σ k in relatively denser network region.In this way the same kind of a priori assignment algorithm will be applied for all the grids, and the only difference reflected in the output uncertainty estimates will be the in situ data strength; the denser the local observation network is, the smaller the uncertainty will be, and vice versa.The overall weighting for GPS data is also rescaled proportionally by a constant factor applied to the solution uncertainties to , where α i is the input horizontal or vertical velocity uncertainty for the i-th GPS site in the region, and β i is the rescaled corresponding uncertainty for the interpolated velocity component at the site, respectively.The sum is over all the GPS sites in the study region.The only free parameter in the derivation is therefore σ 0 , which is set by the user and should be on the order of average spacing of the GPS station network.

InSAR data processing, LOS rate and uncertainty estimation
Here we briefly describe InSAR processing and analysis steps for the InSAR data used in the case study for southern California.We processed the raw SAR data of ERS-1,2 and Envisat satellites from 1992 to 2010 for interferograms using a modified version of JPL/Caltech ROI_PAC software package.Major processing steps include interferometric phase flattening using precise orbit, topography phase correction based on 2-arc SRTM digital elevation model (DEM), baseline re-estimation for orbital error correction when needed, phase unwrapping, filtering and geocoding.For the ERS-2 data after 2001 that have Doppler issue due to gyroscope failure, we employ a maximum entropy approach to resolve Doppler ambiguity and identify all usable ERS-2 interferometric pairs.For Envisat ASAR sensors, we correct temporally correlated range ramp error due to long-term local oscillator frequency drift by adopting an empirical approach (Marinkovic and Larsen, 2013).Comparison with in-situ GPS shows that such a correction works well and reduces the RMS error between InSAR and GPS velocities to less than 2 mm/yr (Liu et al., 2014).
We use a variant of the Small Baseline Subset InSAR time series approach to solve for InSAR LOS time series and mean velocity (e.g., Sansosti et al., 2010).We incorporate topography dependent troposphere delay correction, residual DEM error and earthquake offset estimate, and employ spatiotemporal filtering to remove high frequency turbulent troposphere noise (Samsonov, 2010;Liu et al., 2014).Since orbital ramp error for data from the same track is typically limited to a few acquisitions (e.g., Fattahi & Amelung, 2014) and small, we correct only affected interferograms through baseline re-estimation with the constraint of a priori GPS based deformation model.The number of pairs with such correction is much less than the total number of interferograms that went into the analysis.This ensures that the influence of a priori model is negligible.Hundreds of interferograms that meet spatial and temporal baseline criteria are formed and used in the time series inversion.The InSAR data are weighted by their LOS uncertainties.To characterize the uncertainties associated with InSAR deformation map, we adopt a Jackknife variance estimation approach [Efron and Stein, 1981], which provides a reasonable way to account for uncertainties arisen from lacking or missing dates, uncorrected residuals or other noises, and/or the influence of reference pixel and date.

GPS and InSAR velocity data combination
We combine GPS interpolated velocities and InSAR LOS rate data to produce a spatially continuous 3-dimensional velocity field.We first divide the region into rectangle grid cells.At each grid cell, all of the available InSAR LOS rate data from different tracks (with different viewing geometries) are used.For each of the LOS rate images all the pixel data within the grid cell are averaged to produce a mean rate, weighted by the uncertainties.The binned averages are also made for azimuth angle, look direction, and LOS uncertainty (which is averaged the same way as the other observables) associated with the LOS measurements, and used as data inputs of subsequent analysis.
Because of relative measurements and selections of different reference regions, the InSAR LOS velocities usually have offsets between different tracks.The residual orbital error and/or remaining atmospheric phase noise that are not fully corrected may also introduce some subtle ramp difference between tracks.The first step in GPS/InSAR combination is to solve for the offsets/ramps of InSAR images.Since InSAR data provide only LOS measurements from ascending and/or descending viewing geometry, the offset/ramp parameters have to be solved together with the 3-D deformation components, and some GPS data and their interpolated values are therefore needed in the estimate to stabilize the inversion.Because these offset/ramp parameters are correlated with all the deformation parameters in the study area, an optimal estimate of the offsets/ramps means a global solution for all the parameters involved.However, the number of parameters for the 3-D velocity field can be huge, up to millions or even billions depending on the scope of the study area and the size of the grid cells provided, thus it may not be practical and/or even necessary to solve for all the parameters in a single least-squares solution.We therefore include GPS data at only a limited number of grid points in the solution in this step.Two groups of grid points are accounted: the first group includes all the grid points containing direct GPS velocity observations, and the second group involves decimated grid points with multiple InSAR data entries.Incorporation of the data in the second category helps reinforce the solution for the offsets/ramps, but only at decimated grid points (e.g. by a factor of 10 in each dimension in the overlapped regions) would be sufficient for the purpose.
In the second step the components of offsets/ramps are removed from the InSAR LOS data, and the 3D velocity is solved for each grid cell through least-squares regression, with GPS interpolated velocity and LOS data for the cell incorporated.The adaptive GPS and InSAR data uncertainties are used to weight the data input.The GPS vertical data may or may not be used to constrain the final solution, depending on the quality and reliability of the data.

GPS-InSAR combination in Southern California
We apply the GPS-InSAR combination method to a region in southern California covered by 4 ground tracks of ERS and Envisat satellites (Fig. 1).The GPS velocity dataset used in the study is shown in Fig. 2, along with the interpolated velocities and their uncertainties.As described in the previous section, we adopt an algorithm to determine uncertainties of the interpolated velocities, which employs the same degree of smoothing for all the grid cells, and considers GPS network distribution and site-specific uncertainties to determine uncertainties of the velocity solution.Here we set the data weighting threshold W for southern California to be 3, and the smoothing constant σ 0 for uncertainty evaluation to be 20 km.These parameters are determined after some trial-and-error, in accordance with the network spatial density in southern California.We also assign the lower thresholds of uncertainties for horizontal and vertical GPS velocity data input as 0.5 mm/yr and 1.0 mm/yr respectively, considering epistemic errors of the CGPS site velocity estimates.The final solution does not seem to be sensitive to the changes of these values within the same order of magnitude.
Four tracks of InSAR data sets are used in the study (Fig. 3).The data are the LOS velocities from our previous InSAR time series analysis (Liu et al., 2014)

Result discussion
In this study we explored six models of GPS-InSAR data combination with different options of model constraints, such as (a) whether to estimate the ramps of satellite orbits, (b) whether to use default of estimated uncertainties to condition the InSAR data, and (c) whether to use GPS vertical velocities to constrain the final solution of vertical velocity field.Comparing all the solutions, we find that the GPS-InSAR combined horizontal velocity fields of the six models are very similar to the GPS interpolated horizontal velocity field, and the differences are at the submillimeter per year level for all the data points.The results suggest that the horizontal velocity solution is mostly resolved by the GPS data, and contribution from the InSAR data is relatively minor.Consistency of all the model results also suggests that InSAR and GPS observations are in good agreement in documenting the horizontal deformation field, with both velocity solutions deduced using data of overlapped time span of 6-20 years.The difference in model constraints and/or parameterization, however, can have significant impact on vertical velocity solution and its error assessment.One of the factors involved in the combination is whether to use the GPS vertical data to constrain the pixel solution.Figs. 4 and 7 show the velocity solutions of Models A and C, for which all the parameterizations are the same except that Model C incorporated GPS vertical data to constrain the model and Model A did not.
Comparison of the two solutions reveals that, although inclusion of the GPS vertical data has provided additional constraints to the solution, its lack of detailed spatial resolution smeared and missed some regional deformation signals.For example, up to 8 mm/a subsidence is shown in the Lancaster, Los Angeles basin, and San Gabriel basin regions in the solutions without using GPS vertical data as constraints (Model A, Fig. 4), which are however absent or significantly suppressed in the solutions using GPS vertical data constraints (Model C, Fig. 7).These signals, detected by InSAR observations are caused by ground water withdrawal and shallow aquifer compaction (e.g., Galloway et al., 1998;Hoffmann et al., 2003) and are real, but cannot be picked up by GPS due to limited network spatial coverage (or missed time window).
The GPS data, on the other hand, provide effective constraints for mid to long range vertical deformation (>100 km in scale), associated with earthquake cycle and tectonic deformation.This is evidenced by the vertical deformation pattern shown in Fig. 4, which is similar to that reported t by Howell et al. (2016).We therefore use GPS vertical data to correct for the offsets/ramps of the InSAR data and to stabilize the long range deformation, but not to use that to constrain the local deformation.
Two sets of InSAR LOS data errors are adopted to constrain the solution in this study.One set of solutions assumes a uniform data error of 2 mm/a (Models B, D, and F), which is a common practice when no detailed error analysis is available.Another set of solutions takes the estimated uncertainties derived using the Jackknife variance estimation approach (Models A, C, and E).
Using the estimated uncertainties to weight the InSAR data, the result shows no noticeable difference from the one assuming uniform InSAR data uncertainty (e.g.Fig. 4 vs.Fig.6).
However, solution uncertainties hence derived for the two kinds of models are quite different.
The models assuming uniform LOS rate error deduce the uncertainty estimates with a spatial pattern dictated mostly by InSAR data coverage, i.e. the redundancy distribution of the observations (e.g.Fig. S1f).The models using the estimated LOS rate error yield the uncertainty estimates taking into account of the InSAR data quality and observation history, reflecting better the true data strength and weakness.For example, the solution uncertainty estimates shown in Fig. 5f illuminates not only the impact of spatial pattern of InSAR data redundancy, but also the relative strength of the data input, such as the largest uncertainties (up to 4 mm/yr) in the northwest corner of the studied region, resulted from weak data entry of track 120.
We test different ways to remove the orbital effect from the InSAR data, and examine how that affect the data fitting of the model.Two model parameterizations are tested, one is to solve for an offset (i.e.Models E and F), and another is to solve for a ramp and an offset (i.e.Models A and B) for each of the InSAR data images respectively.The results show that by adding two free parameters for each track of the InSAR data, the orbital ramp model is able to reduce the data postfit residual chisquares by half with respect to the orbital offset model (see statistics in Table 1), attesting the necessity of ramp correction in a joint inversion involving multiple InSAR data entries.Significant jumps can also be seen for vertical solutions across some of the InSAR data boundaries for the models adopted offset correction only (Figs.S5 and S7), which however are much reduced for the vertical solutions of the models adopted ramp corrections (Figs. 4 and 6).
And for the data postfit residual plots, the ones for the offset removal model show significant jumps at the edges of image overlaps (Figs.S6 and S8), which however are much suppressed for the ones in the ramp removal model (Figs. 5 and S1).

Result interpretations
Based on the above discussion, we think that Model A takes the most optimal approach, and its result is therefore the basis for our following interpretation (Fig. 4).
For the region in southern California under investigation, the horizontal velocity solution is mostly determined by GPS data, with the formal uncertainties below 0.7 mm/yr for most of the area (Fig. 5).The highest velocity gradient appears across the San Andreas and San Jacinto fault system, consistent with previous findings about the deformation pattern in southern California (e.g.Feigl et al., 1993;McCaffrey, 2005;Zeng and Shen, 2017).The formal uncertainties for the vertical component are mostly determined by InSAR data, with < 1.5 mm/yr uncertainties for most of the region with more than one LOS data entry, and < 2.5 mm/yr for most of the regions with only one LOS data entry (Fig. 5).One significant exception is for the region near southeast Central Valley at the northwest corner of the study region, where the formal uncertainties are up to 4.5 mm/yr.This is due to relatively short duration and fewer observations for the A120 track of InSAR data.Close to zero residuals usually appear at edges of the study region, where only data from a single InSAR LOS image are available, and the solution uncertainties are relatively larger.
Local subsidence is found at several locations in southern California, such as the Los Angeles basin, Lancaster area in western Mojave, Coso geotherm site, Searles Lake, San Gabriel basin, Death Valley, Palm Springs, and area spanning the southern sections of the San Jacinto and Elsinore faults (Fig. 4).The subsidence ranges 3-8 mm/yr, and possibly caused by the loss of ground water or contraction of geothermal/volcanic activity.Most of these subsidence features are recorded by more than one SAR images, and reliable.About 1-2 mm/yr subsidence appears across the southern plate boundary fault system including the San Andreas, San Jacinto, and Elsinore faults, which is slightly higher than most of the GPS observed vertical velocities.The result is mainly derived from the southeast edge of the image of track A349, and suffers from relatively larger uncertainties (~1.8 mm/yr, Fig. 5f).More SAR data coverage in the region is needed to further confirm the feature of deformation.
Scattered uplift of about 1-3 mm/yr appears in southern Great Valley, and may be due to hydrologic effect associated with drought and crust rebound of the region (Fig. 4).The solution uncertainties however are ~2-3 mm/yr and impede further interpretations.Uplift of about 1-3 mm/yr is also found from northern San Jacinto Mountains across the Banning and Northern San Andreas faults to southern Mojave Desert.The result is in general consistent with the GPS vertical measurements and seems to be credible, with the solution uncertainties less than 1 mm/yr (Fig. 2).The area around the east end of the Garlock fault shows 2-4 mm/yr uplift, which however is not consistent with the GPS vertical velocities in the region.This deformation pattern is solved with InSAR data from the descending track 399 only, with the solution uncertainties of ~2 mm/yr.Input of more InSAR data from this area will help resolve deformation pattern of the region.

Conclusions
We devise an algorithm to optimally combine GPS and InSAR data and produce 3-dimensional velocity field at Earth's surface.At the locations where both InSAR and interpolated GPS data are available, optimal 3-dimensional velocity components are derived using a weighted leastsquare method.Both GPS and InSAR data uncertainties are used to weight the observables in joint inversion.A GPS-InSAR combination code is provided for public use.This algorithm is applied to modeling deformation field at a selected region in southern California.Conclusions are the following.
1. Using optimally estimated GPS and InSAR uncertainties to weight the data provides proper accounting of the solution uncertainties, and helps adequately assess the solution quality and reliability.
2. Including InSAR data from both ascending and descending viewing geometry, if available, provides improved constraint on the 3-D deformation when integrating with GPS data.
3. The approach of using GPS vertical data to constrain deformation field should be subject to evaluation of data quality and deformation pattern.In southern California, the current GPS network is still too sparse to adequately detect localized vertical deformation, particularly in regions affected by hydrologic processes.Existence of certain outliers in the dataset makes identification of localized deformation even more challenging.The optimal approach is therefore to use the GPS vertical data to constrain the satellite orbital ramps only, and leave the localized vertical deformation solved by InSAR, aided by GPS horizontal constraints.

Figure 1 .
Figure 1.Study area in southern California.Black curves are active faults, red and blue squares are GPS sites whose 3D and 2D (horizontal only) data are used in this study respectively.The green frames denote the imprints of 4 InSAR tracks whose data are used in this study.

Figure 2 .
Figure 2. GPS velocities and interpretation result.(a) White vectors are GPS horizontal velocities in SNARF reference frame that are used in the combination with InSAR data.The background colors denote the amplitudes of interpolated horizontal velocity field.(b) Filled circles are GPS vertical observations, and the background colors denote the interpolated vertical velocity field.(c) and (d) are uncertainties of east and up components of interpolated GPS velocities, respectively.Magenta triangles denote the locations of GPS sites.

Figure 3 .
Figure 3. InSAR LOS data from 4 selected tracks of D170, D399, A349, and A120 that are used in the combination.The upper panel shows the LOS velocities, and the lower panel shows the corresponding uncertainties, respectively.
, including the following: (a) descending track 170 derived from ERS/Envisat data over the period of 1992-2010; (b) descending track 399 from Envisat over the period of 2003-2010; (c) ascending track 349 from Envisat over the period of 2003-2010; and (d) ascending track 120 from Envisat over the period of 2003-2010.The lower panel of Fig.3shows the estimated uncertainties for the LOS velocity data, which shows that although uncertainties are relatively uniform for track 170, they vary considerably for tracks 399, 349, and 120.For tracks 399, 349, and 120, we only use Envisat data for interseismic velocity estimates as these tracks spanning the East California Shear Zone (ECSZ).The ERS data from these tracks are not used because they are likely affected by postseismic deformation following the 1992 Landers and 1999 Hector Mines earthquakes in the ECSZ area.This is resulted in fewer SAR images for tracks 399, 349, and 120 than for track 170, and among which a few images were affected significantly by atmospheric disturbance with strong spatial variations.This is particularly true for track 120, with the residual atmosphere noise resulting in the largest errors for the northern part of the track where the LOS rate uncertainties are up to 4 mm/yr.A suite of combination models are tested with various selection of model parameters, including the choices of InSAR data uncertainties, the use of GPS vertical data for model constraints, and the InSAR offset/ramp estimation, etc.

Figure 4 .
Figure 4. Combined GPS and InSAR 3-D velocities for Model A, with estimated InSAR data uncertainties to weight the data and SAR satellite orbital ramps estimated.(a) shows amplitudes of the horizontal components, and (b) the vertical components, respectively.Round dots in (b) are GPS vertical velocities, which are used in the orbital ramps estimation but not the 3-D velocity solution.Name abbreviations: CG, Coso Geotherm site; DV, Death Valley; LAB, Los Angeles Basin; LC, Lancaster; PS, Palm Springs; SB, San Gabriel Basin; SL, Searles Lake.

Fig. 4
Fig.4shows the result of model A, which has 3 common parameters solved for each InSAR

Figure 5 .
Figure 5. Data postfit residuals and solution uncertainties for Model A. (a), (b), (c), and (d) are InSAR LOS postfit residuals for tracks D170, A349, A120, and D399, respectively.(e) and (f) are solution uncertainties for the east and vertical components respectively.

Figure 6 .
Figure 6.Combined GPS and InSAR 3-D velocities for Model B, with default InSAR data uncertainties (2 mm/yr) to weight the data and SAR satellite orbital ramps estimated.(a) shows amplitudes of the horizontal components, and (b) the vertical components, respectively.Round dots in (b) are GPS vertical velocities, which are used in orbital ramps estimation but not the 3-D velocity solution.

Figure 7 .
Figure 7. Combined GPS and InSAR 3-D velocities for Model C, with estimated InSAR data uncertainties to weight the data and SAR satellite orbital ramps estimated.(a) shows amplitudes of the horizontal components, and (b) the vertical components, respectively.Round dots in (b) are GPS vertical velocities, which are used in both estimation of orbital ramps and the final 3-D velocity solution.
4. The GPS and InSAR data are generally consistent for the horizontal velocities at submillimeter per year level.The vertical velocity field is determined much better of the combined solution than using GPS data only, especially for regions experiencing localized deformation.These regions include the Los Angeles basin, San Gabriel basin, Lancaster, Palm Springs, Searles Lake, and Death Valley where hydrologic processes caused induced subsidence of up to 3-8 mm/yr.They also include the southern Great Valley region which underwent drought related uplift of 2-3 mm/yr.Uplift of 1-3 mm/yr is detected across a transect from the northern San Jacinto Mountains to southern Mojave Desert.
The InSAR data delineate an area of approximately 32.5°-36°N, 116°-118.5°W,covering the southern part of the Eastern California Shear Zone, the Mojave Desert, the central Transverse Ranges, and the coastal area from Los Angeles to San Diego.Active faults in the region include the Mojave and San Bernardino segments of the San Andreas, Garlock, Mojave Shear Zone, Owens Valley, San Jacinto, and

Table 1 .
Combination model results.