Combination of coseismic displacement fields : a geodetic perspective

This study provides the mathematical framework for the rigorous combination of coseismic offsets observed by a global positioning system (GPS) network and investigates the results obtained on the occasion of the recent Emilia earthquakes (Italy). This seismic sequence that affected northern Italy from May 20, 2012, allowed two offset fields to be computed, one with reference to the mainshock (M 5.9, followed by two other M 5.1 events on the same day), and a second with reference to the replicas that occurred on May 29, 2012 (M 5.8, M 5.3 and M 5.2; ISIDe data archive, http://iside.rm.ingv.it). The final displacement field is basically the result of a comparison and validation process with repeated feedback between the different analysis groups at the Istituto Nazionale di Geofisica e Vulcanologia (INGV; National Institute of Geophysics and Volcanology) that was established to obtain prompt coseismic displacement solutions, as precise as possible, and in the first days after an event. This is important for early seismic-source evaluation as it represents the most complete and validated dataset at the very early stage of a seismic crisis, and it is also extremely useful in reducing random and systematic errors in the estimated parameters. This study is the result of a cooperative effort that involved different research groups at INGV, with the sharing of all of the collected GPS data. The intention was to compare these results and thus reducing sources of error associated with individual processing strategies, to allow the final combination of the different displacement fields into a single consensus solution. The process assessed the robustness of each single GPS result, thus minimizing erroneous interpretations of individual solutions. […]


Introduction
This study provides the mathematical framework for the rigorous combination of coseismic offsets observed by a global positioning system (GPS) network and investigates the results obtained on the occasion of the recent Emilia earthquakes (Italy).This seismic sequence that affected northern Italy from May 20, 2012, allowed two offset fields to be computed, one with reference to the mainshock (M 5.9, followed by two other M 5.1 events on the same day), and a second with reference to the replicas that occurred on May 29, 2012 (M 5.8, M 5.3 and M 5.2; ISIDe data archive, http://iside.rm.ingv.it).The final displacement field is basically the result of a comparison and validation process with repeated feedback between the different analysis groups at the Istituto Nazionale di Geofisica e Vulcanologia (INGV; National Institute of Geophysics and Volcanology) that was established to obtain prompt coseismic displacement solutions, as precise as possible, and in the first days after an event.This is important for early seismic-source evaluation as it represents the most complete and validated dataset at the very early stage of a seismic crisis, and it is also extremely useful in reducing random and systematic errors in the estimated parameters.This study is the result of a cooperative effort that involved different research groups at INGV, with the sharing of all of the collected GPS data.The intention was to compare these results and thus reducing sources of error associated with individual processing strategies, to allow the final combination of the different displacement fields into a single consensus solution.The process assessed the robustness of each single GPS result, thus minimizing erroneous interpretations of individual solutions.For more information about the interpretation of the coseismic displacements, see Serpelloni et al. [2012, this volume].The combination procedure can be easily generalized to other geodetic parameters, and might be useful whenever a validated dataset is desirable or is intrinsically required (i.e., near real-time observations of unknown surface processes, tectonic motion of large geographic areas, angular momentum observations of the whole Earth, coseismic surface deformation patterns).
Other important advantages of the combination process include efficient outlier detection via cross-checking of independent solutions, and realistic estimates of the displacement uncertainties.It is important to emphasize here that in this kind of activity, the cooperative behavior between individual researchers and research groups should prevail over competitive and self-assertive behaviors, because the individuals (input solutions) become part of the same objective (validation).On the other hand, the validation process should not overlap or replace the research activity; i.e., there should always be space for individual initiatives.In this respect, the cooperative effort spent in constructing a prompt and validated product is a distinctive activity that involves both technological and scientific aspects, and it certainly stimulated the individual research activities.
Methods for comparison and combination of geodetic network solutions were introduced and developed in several studies in the early 1990's [Heflin et al. 1992, Boucher and Altamimi 1993, Davies and Blewitt 2000, and references therein].It is a particular least-squares task in which the input data types are the same as the output parameter types and the covariance matrices of the input datasets are fully populated (because each dataset is the output parameter set of an earlier least-squares analysis).The parameters that can be combined can be three-dimensional (3-D) station positions at a stated epoch, constant 3-D station velocities, Earth-orientation parameter time series, or coseismic station offsets.In principle, each of the aforementioned parameter types, or a subset of them, can be treated in the same mathematical framework and combined into a unique average solution with only slight modifications to the formulation given below.were computed in the International Terrestrial Reference Frame (ITRF)2008 reference framework [Altamimi et al. 2011], which provided three distinct time series of a regional GPS network located around the epicenter, in north-eastern Italy (see Figures 1, 2, for maps of the GPS network).To estimate the coseismic offsets we considered only a single day of observations after the event, to minimize the effects of early post-seismic displacements, which can be eventually caused by afterslip or poroelastic rebound [e.g., Peltzer et al. 1996, Mahsas et al. 2008].Fast exponential decays in the very early phase of the seismic cycle are commonly observed in a wide range of seismic events that cause deformations of up to several centimeters, depending on the magnitude and the driving mechanism [Segal et al. 2000, Yu et al. 2003, Devoti et al. 2012a].Therefore, the coseismic offsets were estimated using 15 days before the epoch of the earthquake (whenever possible), and one day after the event, with the following days being used only to check the consistency of the coseismic displacement at the given station.We checked for outliers in the time series using the detrended 3-D position residuals, using different criteria for before and after the event.Before the earthquake, a point was rejected if its residual exceeded the weighted-root-mean-squared (WRMS) of all of the residuals by three times; after the event instead, if the difference between consecutive points exceeded 3 times the WRMS, we computed the weighted average between them.This latter condition was applied to only a few sites in the near field: CONC, MSEL and SBPO.All of the offsets were estimated directly in the same least-squares inversion, to obtain the estimates of the vector and the associated network covariance matrix for each of the three processing schemes.Figures 3  and 4 show the residuals of the estimated station offsets (namely the WRMS of the residuals between pairs of solutions) in the near field (<50 km from the epicenters).This parameter measures the repeatability of the estimates across independent solutions, and it is extremely useful for the isolating of systematic errors of the individual solutions.For the May 20, 2012, event, the overall average WRMS was 2.1 mm in the horizontal plane and 4.2 mm in the vertical direction; for the May 29, 2012, event, the WRMS was 2.6 mm in both components.The worst case was observed for station MO05 (Finale Emilia), where the GIPSY solution differed on average by 15 mm from the other two, in the vertical direction.Nevertheless, given the high uncertainties of the vertical GPS estimates, this difference can still be considered acceptable.

Displacement field combinations
Consider the linearized observation equation: ( where y is the column vector of the observations, A is the design matrix, x is the column vector of the functional model   A y x = +o parameters, and o is a column vector of random errors.In our case, the observations and the model parameters are of the same nature, and therefore the linearized observation equations are simply: (2) where X s and X c are the offsets (observations) of the solution s and the combined parameters c (estimates), respectively, while the design matrix is a simple reordering matrix, as observed by Davies and Blewitt [2000], which reparameterizes the parameter list from the combined offsets to the solution offsets (i.e., the component a ij = 1 if element i of list c is the same as element j of list s).A notable property of the design matrix is that: (3) i.e., transposing the reordering matrix from list c to list s corresponds to reversing the ordering from list s to list c.To simplify the mathematical formulation, we adopt the convention that ; i.e., the design matrix of the i-th solution (the reordering matrix from solution s = i to the combined solution c) is denoted simply as A i .
The combination of the individual solutions is obtained by solving for the unknowns X c in a least-squares sense, and a useful recursive solution can be easily derived: (4).
From each input solution (usually in SINEX format), we can extract the estimate parameter vector X i and its covariance matrix C i , compute the reordering matrix A i , and finally cumulate the normal matrix and vector, as the first and second factors, respectively, of the right-hand side of Equation ( 4).
The relative weight of the input solutions is imposed by the condition that the contribution to the total | 2 of each individual solution is equally balanced [Devoti et al. 2012b]: (5).
From the first preliminary combination, a variance factor f i is estimated for each solution, and then applied to the covariance matrix of the corresponding solution in an iterative way, until the global | 2 reaches unity.This procedure assures that each solution contributes equally to the global | 2 , which avoids an individual solution from dominating or any suffering from the combination of the individual solutions (i.e., over-weighted or down-weighted solutions).
It might happen that the GPS stations are geographically clustered, especially in urban areas where there are different installations over inter-distances of a few kilometers.In this case, it is useful to constrain the clustered offsets, es-timating them only once.If the inter-distances of the GPS stations are small with respect to the coseismic deformation pattern, i.e., of the order of 1 km, then clustered observations will rather bias any inverse solution because repeated information at the same site might artificially over-weight that particular site, thus providing biased solutions.This can be handled during the combination estimation process by imposing the constraint: (6) whenever sites i and j have to be estimated together.The constraints were imposed using the classical Lagrange multipliers technique [e.g., Tarantola 2005] that provides the constrained least-squares solution with a fully populated covariance matrix of the combined parameters.

Results for the Emila earthquake
Figures 1 and 2 show the combined offsets (red arrows) of the two seismic events (the May 20 and 29, 2012, events, respectively) projected in the local horizontal and vertical directions.The maximum horizontal offset was observed at MO05, the closest site to the mainshock epicenter (horizontal, 30 mm; vertical, 73 mm).The offsets caused by the May 29, 2012, events were generally smaller than the offsets caused by the mainshock, although it has also to be noted that the available GPS network did not correctly sample the expected surface deformation, and therefore the measured offsets might be aliased.The individual solutions are also shown with different colors in Figures 1 and 2, to have a first glance of their repeatability.Apart from any apparently random differences, we will show that systematic bias might also characterize the individual solutions.A first validation of the individual solutions was carried out by plotting the offsets measured outside an area of 50-km radius from the epicenter.Ideally, the far-field offsets should be randomly distributed (if not affected by the earthquake), and they represent a good test case to determine any systematic errors.Figures 5 and 6 show the histograms of the far-field offsets for the two events, where the most notable features are non-zero medians for the individual solutions; e.g. for May 20, 2012, a median vertical component of -2.1 mm (GAMIT), 3.0 mm (GIPSY) and north component of 2.1 mm (GIPSY), and less pronounced bias for the May 29, 2012, event.The combined offsets show instead a general bias reduction of the far-field offsets, indicating that the combined solution is in fact the best unbiased representation of the displacement field.The combined solution also shows a significant reduction in the formal errors (Table 1): both events show (on average) 1.2 mm one-sigma error in the horizontal plane, and about 3.1 mm in the vertical direction, whereas the individual solutions show (on average) a 50% higher sigma.

Figure 1 .
Figure 1.Combined (red arrows) and individual solution offsets of the horizontal (top) and vertical (bottom) components associated with the May 20, 2012, seismic events.The individual solutions are represented in green (BERNESE), cyan (GIPSY) and blue (GAMIT).Stars, epicenters of the M >5 events; circle, delimitation of the 50-km range from the mainshock (M 5.9); error ellipses, one-sigma confidence regions.

Figure 2 .
Figure 2. Combined (red arrows) and individual solution offsets of the horizontal (top) and vertical (bottom) components associated with the May 29, 2012, seismic events.The individual solutions are represented in green (BERNESE), cyan (GIPSY) and blue (GAMIT).Stars, epicenters of the M >5 events; circle, delimitation of the 50-km range from the main event (M 5.8); error ellipses, one-sigma confidence regions.

Figure 3 .
Figure 3. Cumulative WRMS of the station offset differences between the individual May 20, 2012, solutions (BERNESE, GAMIT and GIPSY) for the horizontal (gray) and vertical (open bar) components.

Figure 4 .
Figure 4. Cumulative WRMS of the station offset differences between the individual May 29, 2012, solutions (BERNESE, GAMIT and GIPSY) for the horizontal (gray) and vertical (open bar) components.

Figure 5 .
Figure 5. Distribution of the far-field offsets for the individual (BERNESE, GAMIT and GIPSY) and COMBINED May 20, 2012, solutions.Red line, bestfit Gaussian curve.The median (m) and dispersion ( ) are computed for the three spatial directions (up, east, north).

Figure 6 .
Figure 6.Distribution of the far-field offsets for the individual (BERNESE, GAMIT and GIPSY) and COMBINED May 29, 2012, solutions.Red line, bestfit Gaussian curve.The median (m) and dispersion (v) are computed for the three spatial directions (up, east, north).

Table 1 .
Average formal errors (1-sigma) of the estimated coseismic offsets.The two values represent the horizontal and vertical components, respectively, of the uncertainties.