Inversion of synthetic geodetic data for the 1997 Colfiorito events : clues on the effects of layering , assessment of model parameter PDFs , and model selection criteria

The 1997 September-October Umbria-Marche sequence has been extensively studied in the past by analyzing coseismic displacement data (GPS, leveling, SAR). Here we focus on synthetic data representative of the main event of the 1997 Umbria-Marche sequence and investigate the effects of a crustal layering proper to the Colfiorito area on surface displacements and inferred source features when inverting coseismic geodetic data without taking into account layering. We compare bootstrapping and NA-Bayes as tools for parameter uncertainty assessment and show how the Akaike Information Criterion can be used to select the model which is most likely to be correct. Since SAR images offer the most complete coverage of the study area, we use synthetic line-ofsight displacement data. Mailing address: Dr. Antonella Amoruso, Dipartimento di Fisica «E.R. Caianiello», Università degli Studi di Salerno, Via S. Allende, 84081 Baronissi (SA), Italy; e-mail: antonella.amoruso@sa.infn.it Vol51,2_3,2008 4-03-2009 10:28 Pagina 461


Introduction
The 1997 Umbria-Marche sequence includes two moderate magnitude earthquakes (September 26th, Mw = 5.8 and Mw = 6.0) that occurred near the village of Colfiorito followed by a sequence of aftershocks lasted several months.The aftershock activity included a Mw = 5.6 earthquake (that occurred near the village of Sellano on October 14th) and three events with magnitude larger than 5.2.The Colfiorito sequence has been extensively studied using geologic data, strong motion data, teleseismic and regional waveforms, GPS, levelling and SAR data (see e.g.Cinti et al., 1999;Capuano et al., 2000;Pino et al., 1999;Hunstad et al., 1999;Lundgren and Stramondo, 2002;De Martini et al., 2003;Crippa et al., 2006).Larger shocks share SW strikes and normal fault mechanisms.
After the September 26th event, 12 monuments from the IGMI95 GPS network within 30 km off the epicentre were reoccupied, including FOLI and COLF that were occupied continuously during the survey.Coseismic diplacements were obtained by comparing the 1997 survey coordinates with those collected in 1995.Two sites (CROC, -13.8 cm N and -1.9 cm E; PENN, 7.6 cm N and 8.0 cm E) showed significant horizontal displacements and one (CROC, -24.7 cm) showed significant vertical displacements (Hunstad et al., 1999).
The epicentral area is crossed by Line 21 of the Italian first-order levelling network (operated by IGM), which runs from Foligno, central Umbria, to Fiumesino, on the Marche coast; the first (southernmost) section of Line 21 trends almost N-S between Foligno and Fossato di Vico and includes about 40 bench marks.The route was measured in 1951, 1992 and was resurveyed in 1998.The data showed maximum subsidence of 7.8 cm, and maximum uplift of 1.6 cm (De Martini et al., 2003).
Differential interferograms related to the 1997 Colfiorito events have been obtained from ERS-1 and ERS-2.Because of the strong topographic relief of the Apennines and the significant vegetation cover of the epicentral area, even the best interferograms show large decorrelation areas (Salvi et al., 2000).For this reason, even though a large number of SAR data are available from July 1993 to October 1997, only few significant interferograms have been obtained, covering different periods of the sequence.The most studied interferogram is the 35-day ERS-2 one (e.g.Stramondo et al., 1999;Crippa et al., 2006) covering the period 7 September to 12 October 1997.
Geodetic data are usually affected by both correlated and uncorrelated noise.For example, levelling data are obtained by measuring the elevation difference between consecutive bench marks (section height differences).Bench mark elevations with respect to the reference bench mark can then be computed, and are often used instead of section height differences for modeling.This operation introduces nondiagonal terms in the covariance matrix of bench mark elevation differences, whose diagonal terms are related to uncorrelated noise (e.g.Amoruso and Crescentini, 2007).
As regards SAR data, the unwrapped differential phase is also affected by both uncorrelated and correlated noise.Correlated noise is caused by atmospheric disturbances, and can be expressed by an exponential autocorrelation function.This noise can be handled following two different schemes: using a non-diagonal covariance matrix in the data inversion proce-dure (Fukushima et al., 2005) or removing from the interferogram the influence of atmospheric heterogeneities after considering stable areas close to the epicentral area to analyse the spatial autocorrelation (e.g.Crosetto et al., 2002).We prefer the former scheme, since, as stressed in Crosetto et al. (2002) the removal procedure «suffers for its intrinsic limitation related to the non-stationarity of the analysed signal, which practically confines its use to the vicinity of the stable areas».To give an idea of overall noise in SAR data related to the Colfiorito events, Lundgren and Stramondo (2002) assigned a standard deviation of 1 cm to each contoured data point (i.e., points of line-of-sight, LOS, displacement contours at 28 mm intervals).
If both correlated and uncorrelated noise is present, measured data are not statistically independent.However, the data covariance matrix is symmetric and can be reduced to diagonal form by means of a rotation matrix; the same rotation transforms the data to independent form (i.e. a set of statistically independent data, which are linear combinations of the measured ones).The eigenvalues of the covariance matrix give the uncertainties of the transformed independent data.This procedure can also be used in case of non-normal distribution of errors and consequently for robust fitting (Amoruso and Crescentini, 2007).
Model parameters are usually obtained from experimental data through the appraisal of a cost function which measures the disagreement between data and model predictions.
Several different causes contribute to differences between predictions and observations (the residuals of the retrieved model): measurement errors, intrinsic misfits (e.g.misfits due to heterogeneities in slip distribution, non-planar faults, heterogeneities in the elastic properties of the medium, when the adopted model is too simple, like a planar uniform-slipping fault in a homogeneous half-space), local response effects (e.g.soil compaction, landslides, topographic effects).For example, Amoruso et al. (2004) generated synthetic geodetic coseismic data (GPS, SAR, and levellings) in a layered medium and performed inversions of the synthetics, under the assumption of homogeneous half-space.Horizontal displacements are more affected than vertical ones by the presence of a superficial soft layer, but differences with synthetics generated in a homogeneous half-space did not show any simple regular behaviour, even if a few features could be identified.Consequently, also retrieved parameters of the homogeneous equivalent fault obtained by unconstrained inversion of surface displacements did not show a simple regular behaviour.Amoruso et al. (2004) pointed out that the presence of a superficial layer may lead to misestimating several fault parameters both using joint and separate inversions of the three components of synthetic displacement.The effects of the presence of the superficial layer depend on whether all fault parameters are left free in the inversions or some of them are fixed a priori.In the inversion of any kind of coseismic geodetic data, fault size and slip could be largely misestimated, but the seismic potency (product of ruptured area and average slip) is well determined (within a few per cent) in most cases even neglecting layering.Of course, the amount of the effect of the superficial low-rigidity layer depend on the contrast of the elastic parameters.For example, Battaglia and Segall (2004) and Crescentini and Amoruso (2007) have shown that layering effects are negligible and relevant respectively when considering two different volcanic environments and layering properties (Long Valley and Campi Flegrei).Assessing the effects of layering is a case-to-case problem.
Intrinsic misfits and local response effects, as well as possible systematic errors, could be spatially coherent over small or large areas, and are not zero-mean random variables.Rigorous handling of data should take into account, and properly treat, all error sources (measurement errors, systematic errors, intrinsic misfits and so on), after a reliable a priori estimate of all the errors, and the determination of how much different error sources effectively affect data.This is quite a hard task, largely case-to-case dependent.In case of levelling data, if measurement errors can be reliably estimated, after inverting the transformed independent data, reduced chi-square of residuals can be computed and used to refine uncorrelated error variance when the number of degrees of freedom is easily computable, e.g. in case of one uniform-slip rectangular fault in a homogeneous half-space (Amoruso and Crescentini, 2007).
The cost function is usually obtained from maximum-likelihood arguments according to the assumed statistical distribution of the residuals.Misfit is written in terms of the transformed independent data set, and, usually, either the mean squared deviation of residuals (chi-square fitting, M 2, proper for normally distributed errors) or the mean absolute deviation of residuals M1, proper for two-sided-exponentially distributed errors and robust fitting) are used.If possible, it is preferable to compare results obtained using M 1 and M2.
Estimates of the model parameters can be obtained by the best-fitting approach, i.e. the minimization of the cost function.When inverting geodetic data, the cost function depends nonlinearly on parameters and is characterized by a rough landscape and several local minima.Monte Carlo techniques are the only efficient tool to search for the global minimum, but require computing a large number of forward models (data prediction given a set of model parameters).
Thus the use of fast codes for forward modelling is preferable if not essential.Among the Monte Carlo global optimization methods (e.g.uniform sampling, genetic algorithms, simulated annealing) we prefer to use Adaptive Simulating Annealing (ASA, Ingber 1993).Generally speaking, each step of any Simulated-Annealing algorithm replaces the current point in the parameter space by a random «nearby» point, chosen with a probability that depends on the difference between the corresponding costfunction values and on a global parameter (T) that is gradually decreased during the process according to some annealing schedule.ASA considers that different parameters might require different annealing schedules.The exponential annealing schedules used by ASA ensure ample global searching in the first phases of search and ample quick convergence in the final phases.Even if success in finding the global minimum is never guaranteed and inversions from different starting points and/or using different program options can lead to different endpoints, the ability of ASA to self optimize its program options recursively makes it a very powerful technique.
The best-fitting approach does not give parameter uncertainties and consequently is of limited significance.A commonly used technique (associated with global minimization) to assess acceptable parameter ranges consists in the use of the bootstrap percentile method.This method applies the best-fit technique to a large number of synthetic data sets, obtained by randomly resampling from the actual data set with replacement and is often used to estimate confidence intervals without making assumptions about the underlying statistics of errors.
Bootstrapping also estimates correlations between different parameters, by forming a scatter plot of all the estimates for each parameter pair and visualizing the correlations between the parameter pairs (e.g.Amoruso et al., 2005).The computed values for the model parameters from the synthetic data sets form an estimate of the sampling distribution of the parameter values, independently from the underlying statistics.This method gives reliable results only if the hypothesis of independent identically-distributed data is verified.As a consequence, in the presence of a non-diagonal data covariance matrix, bootstrapping can only be applied to the independent rotated data, not the measured data.
A different approach is used in the Neighbourhood Algorithm (NA-sampler, Sambridge, 1999a).NA generates ensembles of models which preferentially sample the good data-fitting regions of the parameter space, rather than seeking a single optimal model.The algorithm makes use of only the rank of a data fit criterion rather than the numerical value; in other words, points in the parameter space are listed according to their capability to fit the data (answering the question «does point A give a better fit to the data than point B?») without quantifying the difference in a precise way.With this rank-based approach, the weight of each of the previous sampled points in driving the search depends only on their position in the rank list.The companion code NA-Bayes (NAB, Sambridge 1999b) consists of an algorithm for using the entire ensemble of models produced by NA, and deriving information from them in the form of Bayesian measures of resolution, covariance and marginal probability density functions (PDF) etc.The input ensemble in princi-ple can follow any distribution and be generated by any search method (e.g. the NA-sampler algorithm, a genetic algorithm or simulated annealing algorithm).
In any case, we assume a priori a physical model to explain the data.In principle, any set of data can be almost perfectly fit by using a sufficiently complicated model, but it could be unrealistic and overmodelling (i.e., using more parameters than necessary) has to be avoided.Attempts to find the model that best explains the data with a minimum number of parameters accomplish the parsimony principle.
If nested models (the more complicated model includes the simpler one as a particular case) have to be compared, the F-test approach could be used.The F-test is based on traditional hypothesis testing.Under the normal distribution hyphothesis, F-test compares the differences in χ 2 between nested models with the difference expected by chance, i.e. it quantifies the relationship between the improvement in χ 2 and the relative increase in degrees of freedom g.The null hypothesis is that the simpler model is correct, i.e. that none of the additional parameters is significantly correlated with the model.Unfortunately, the test is unreliable if residuals depart from the normal distribution even slightly.Moreover, the F-test is not valid if the models are not nested or it is difficult to estimate the actual number of g (as in the case of distributed-slip faults, due to the presence of constraints like smoothness or slip positivity).
The Akaike Information Criterion (AIC, Akaike, 1974) is a measure of the goodness of fit of an estimated statistical model.It is based on information theory and does not use hypothesis testing, so there is no conclusion about statistical significance and rejection of a model.AIC attempts to find the model that best explains the data with a minimum number of parameters taking into account both the value of the likelihood function and the number of parameters in the model, i.e. trading off the complexity of an estimated model against how well the model fits the data.AIC includes a penalty that is an increasing function of the number of estimated parameters.The model with the smallest AIC is most likely to be correct.For a small number of observations, a second order correction is added to the general AIC; corrected AIC (AICc) converges to AIC as the number of observations increases (e.g.Hurvich and Tsai, 1989).In addition, the difficulty of estimating the actual number of degrees of freedom can be surmounted by using the Akaike Information Criterion.Even if the results are not automatically generalizable, Amoruso and Crescentini ( 2007) have shown, in case of levelling measurements and the one-fault model of the 1908 Messina earthquake, that the Akaike Information Criterion is not only an effective tool to discriminate if increasing the number of subfaults in a distributed-slip model really improves the model or does not, but also to estimate uncorrelated errors, which partially reflect model inadequacy.
Here we investigate the effects of a crustal layering proper to the Colfiorito area on surface displacements and inferred source features when inverting synthetic geodetic data representative of the main event of the 1997 Umbria-Marche sequence.We compare different methods for parameter uncertainty assessment and we show how the AICc can be used to select the distributed-slip model which is the most likely to be correct as a result of the inversion of available data.Since SAR images offer the most complete coverage of the study area, we use synthetic LOS displacement data.

Generation
To generate synthetic displacements we use the slip distribution in fig. 1, which resembles the recently published distributed-slip model of the main event in Crippa et al. (2006).The source fault (divided into 100-along-strike and 80-along-dip subfaults) is embedded both in a homogeneous half-space (HHS; Okada, 1985)  and the layered half-space (LHS; Wang et al., 2006) detailed in table I (Megna, 2007).Strike, dip and rake are 138°, 45°, and -75°respectively.We then add uncorrelated random noise.As previously mentioned (see Introduction) atmospheric disturbances generate correlated random noise that can be expressed as an exponential covariance function and would require transforming the data to independent form, but treatment of the effects of atmospheric disturbances is beyond the aims of this paper.Displacements (northward, eastward, vertical) were computed on a grid spanning 15000 m toward North, step 100 m, and 20000 m toward East, step 100 m.We simulated SAR-like unwrapped data through a proper combination of the displacement components taking into account the ERS-2 LOS.LOS displacement (positive toward the satellite) maps for HHS and LHS are in fig.2a and fig.2b respectively.Figure 2c shows the residuals of the linear least squares regression between LHS and HHS displacements.As expected (e.g.Amoruso et al., 2004) residuals are about 15% of the signal, and, being obtained using a linear best fit, are not merely due to a translation or a scale factor.LHS LOS displacement after the addition of noise (generated using random variates from a normal distribution with zero-mean and 1-cmstandard deviation) is shown in fig. 3. To reduce the number of data in the inversions, the displacement map was decimated (undersampling interval = 3 bins).
Even if noise added to synthetics was obtained from a normal distribution, model inadequacy is expected to enlarge tails of the residual distribution with respect to the normal distribution and we minimize the mean absolute devia-tion of residuals M 1, commonly used for robust fitting: (2.1) where n is the number of data and vi is the coseismic change estimated for element i of the data set.The data set must be composed of statistically independent data: if they are not (i.e. the covariance matrix is not diagonal) then they have to be transformed accordingly.Here experimental data are not rotated because of the absence of correlated noise.Model predictions v i (a) depend on model parameters a, and wi are the uncertainties of the independent data.

Inversion, uniform-slip model
At first, we invert noisy LHS LOS displacements for a uniform-slip rectangular fault embedded both in HHS and in LHS using ASA.This model involves a parameter set a which consists of the strike and dip angles of the fault, the rake angle of the slip vector, two geometrical fault dimensions (along-strike length and along-dip width), location of start point of the fault upper side along the strike direction (x, positive Northward; y, positive Eastward), depth of the upper side of the fault (z, positive downward), and magnitude of the slip.
Figure 4 shows the two best-fit uniform-slip faults over the distributed-slip source fault.Potency of the best-fit faults is the same inside 2% (Ӎ 2.65 × 10 7 m 3 ), about 9% lower than the potency of the source fault (2.89 × 10 7 m 3 ).The only noticeable difference between the two best-fit faults consists in the expected focusing effect in the HHS: while the LHS fault is centered on the slip patch of the source fault, the HHS fault is shifted by about 700 m upward along the dip direction, i.e. 500 m in depth.PDFs of model parameters were estimated using both bootstrapping and NAB.Strictly speaking, a full optimization procedure should be performed at each resample in the bootstrap procedure, but in our case computation time would be unacceptable.In this work we use a faster method (downhill simplex; Press et al., 1992), starting from the ASA best-fit model to explore the parameter region surrounding it.As regards NAB, the input ensemble has been generated by NA. Figure 5 shows estimated PDFs as well as best-fit parameter values.The true  source values are also shown when pertinent, i.e. for strike, dip, and rake.The best-fit values and PDFs are in agreement with the true values, but strike is somewhat overestimated, owing to the uniform-slip approximation, as it will be clarified below.PDFs obtained using the two techniques are generally in mutual agreement; apparently larger differences relate to depth of the fault upper side and along-dip width in the LHS case.However, fig.6 shows that there is a trade-off between them and the trade-off pattern is very consistent between bootstrapping and NAB.

Inversion, distributed-slip models
As a second step, we consider a small number of coplanar subfaults (2 along strike and 2 along dip, numbered along the strike direction, starting from the shallower part of the fault), which slip independently at the same rake angle, thus adding three additional parameters to the model.Figure 7 shows the two best-fit 2 × 2 faults over the distributed-slip source fault.Potency of the best-fit faults is the same inside 4% (Ӎ 2.6 × 10 7 m 3 ), about 10% lower than the potency of the source fault (2.89 × 10 7 m 3 ).As in the case of the uniform-slip model, the only noticeable difference between the two best-fit faults consists in the focusing effect.
Once again, PDFs of model parameters were estimated using both bootstrapping and NAB, as detailed in the case of the uniform-slip fault.Figure 8 shows estimated PDFs as well as best-fit parameter values and suggests similar comments with respect to the uniform-slip model.PDFs from NAB are more scattered than from bootstrapping, especially for subfault slips.This behaviour is probably related to the cost function roughness: on the one side NA could be in some trouble in sampling the parameter space efficiently, as suggested by its best-fit misfit (larger than that obtained by ASA, using the same number of forward models), on the other, bootstrapping uses simplex, which could be unable to escape from local minima.This problem deserves great attention and deep investigation while inverting real data.As a first simple application of model selection criteria (selecting the model that best ex-plain the data with a minimum number of parameters) we use AIC corrected for small sample sizes (Hurvich and Tsai, 1989) (2.2) to compare the uniform-slip fault model and the 2x2 one.Here n is the number of data, k is the number of parameters, and M is the chosen discrepancy (measure of lack of fit, e.g.Zucchini, 2000); in this case M = M 1 .Difference in AICc  is about -500 both for the HHS and the LHS, indicating that the uniform-slip model is too simple and requires complex models to be tested.
Increasing the number of subfaults would disclose the main features of the slip pattern and possible changes in fault geometry while relaxing the assumption of uniform slip (e.g.Amoruso et al., 2002), but also increases the difficulty of studying the cost function features.In this work, fault plane geometry and rake an-gle show very small changes when turning from the uniform-slip fault model to the 2 × 2 one, thus we use both the geometry and mechanism of the fault plane of the best fitting uniform-slip model to estimate the slip distribution, but we increase the along-strike length and the downdip width of the fault.To enlighten the effects of neglecting layering when inverting geodetic data, we consider the HHS case.The fault is divided into an even grid of p-along- strike x q-along-dip subfaults, with variable slip magnitude and constant rake.We minimize the cost function (2.3) using the algorithm for the least squares problem with linear inequality constraints of Lawson and Hanson (1995).Here ∇ 2 s is a finite difference approximation to the Laplacian of the slip distribution and A is the area of each subfault.It would also be interesting to use misfit M 1 of LOS displacements, but it leads to the prohibitive problem of nonlinear minimization in a huge space.
Finding the optimal model parameters thus requires the simultaneous optimization of two objective functions, namely the fit to the data and the roughness, whose relative importance is controlled by the adimensional smoothing parameter f.The choice of f is a basic problem in inverse theory.Traditionally, it is selected from a «trade-off» curve which plots the trade-off between the two objective functions when varying the value of f.The trade-off curve is monotonic and asymptotes to the minimum value of the roughness (which is zero) at one end and the minimum value of the fit to the data at the other.Choosing a suitable value of f is somewhat arbitrary, but the knee of the trade-off curve is a good compromise between the model complexity and the data fit.Cross validation is a more rigorous criterion (e.g., Matthews and Segall, 1993) but is computationally intensive.Tests performed by Arnadottir and Segall (1994) show that cross validation gives an optimal estimate of f, while the trade-off curve gives a slightly smoother solution.We use the trade-off  estimation since it is somewhat more conservative in smoothness and computationally less intensive than the cross-validation technique.The trade-off curve between the M 2 misfit of LOS displacements and the roughness of the slip distribution is shown in fig.9.It is noticeable that f ഠ 3 × 10 10 is in the knee region of the tradeoff curve independently of the number of subfaults (as far as the number of subfaults is large enough to introduce roughness) suggesting that f can be chosen using a single trade-off curve (i.e. a single p x q configuration) for each real case.Figure 10 shows the retrieved 40 × 32 slip distribution overlying the true one.Once again, the most relevant feature consists in the expected focusing effect in the HHS, since the maximum slip patch is shifted by about 1400 m upward along the dip direction, i.e. 1000 m in depth.
It is important to realize if so many subfaults (40 × 32) are really necessary to represent the slip distribution or, on the other end, if they are enough.We computed AICc for 6 × 4, 8 × 6, 10 × 8, 12 × 10, 15 × 12, 20 × 16, 40 × 32 subfaults.Minimum AICc is obtained for 12 × 10 subfaults.Because AICc is on a relative scale, here we give the AICc differences with respect to the 12 × 10 case, i.e. 739 for 6 × 4, 84 for 8 × 6, 10 for 10 × 8, 93 for 15 × 12, 384 for 20 × 16, 3778 for 40 × 32. Figure 11, upper plots, shows the capability of AICc to find the model that best explains the data with a minimum of parameters, by comparing contours of the slip distribution for 6 × 4, 12 × 10, and 40 × 32, subfaults.Lower plots in fig.11 show the importance of using the correct weighting factor from trade-off curves when combining different terms in the cost function: if f is too small (left plot) the retrieved slip distribution shows unreal small-scale features, if f is too large (right plot) the retrieved slip distribution is too smooth.The same approach should be followed in joint invertion of different kinds of data (see e.g.Amoruso et al., 2008, for EDM and levelling data).

Conclusions
Inversion of synthetic coseismic displace-ment data representative of the main event of the 1997 Umbria-Marche sequence has been used to enlight the effects of a crustal layering proper to the Colfiorito area on inferred source features when layering is not taken into account.As a test case we use a recently published distributed-slip model (Crippa et al., 2006) and we show that the focusing effect of the fault slip distribution can be as large as 1 km.
PDFs of model parameters obtained using bootstrapping and NAB are generally in mutual agreement; apparent larger differences relate to the presence of trade-off between selected model parameters (e.g.depth of the fault upper side and along-dip width).PDFs from NAB are usually more scattered than from bootstrapping, especially for subfault slips.This behaviour is probably related to the cost function roughness: on one side NA could be in some trouble in sampling the parameter space efficiently, on the other we use simplex, which could be unable to escape from local minima, for bootstrapping to reduce computing time.This problem deserves attention in the inversion of real data.
When inverting for the slip distribution, we use a roughness weight that appears independent of the number of subfaults, thus suggesting that its optimal value can be assessed by a single trade-off curve (i.e. a single p x q configuration) for each real case.The number of subfaults, that are necessary and sufficient (from the point of view of the information theory) to account for measurements can be obtained from AICc and in this work is as small as 12 × 10.

Fig. 2 .
Fig. 2. Map of synthetic LOS displacements generated in the HHS (left) and LHS (middle).Right, map of the residuals of the best-fit linear regression between LOS displacements in the HHS and LHS.

Fig. 3 .
Fig. 3. Map of noisy synthetic LOS displacements generated in the LHS.Horizontal projection of the source fault in black.

Fig. 4 .
Fig. 4. Best fit uniform-slip fault in the HHS (thin-line rectangle) and LHS (tick-line rectangle) over the distributed-slip source fault.Numbers give slip.

Fig. 9 .
Fig. 9.Trade-off curves between the M2 of LOS displacements and the roughness of the slip distribution, for different p × q-subfault models.Labels give the weight f of roughness in the cost function.