Non-permanent GPS data for regional-scale kinematics : reliable deformation rate before the 6 April , 2009 , earthquake in the L ' Aquila area

A GPS-based geodetic study at a regional scale requires the availability of a dense network that is characterized by 10 km to 30 km spacing, typically followed in a few continuous GPS stations (CGPSs) and several nonpermanent GPS stations (NPSs). As short observation times do not allow adequate noise modeling, NPS data need specific processing where the main differences between NPSs and CGPSs are taken into account: primarily time-series length and antenna repositioning error. The GPS data collected in the 1999-2007 time-span from non-permanent measurement campaigns in the central Apennine area (Italy) that was recently hit by the Mw 6.3 L'Aquila earthquake (April 6, 2009) are here further analyzed to compute a reliable strain-rate field at a regional scale. Moreover, areas characterized by different kinematics are recognized, and a complete characterization of the regional-scale kinematics is attempted. These new data can be interpreted as indicators from the viewpoint of seismic risk assessment.


Introduction
Global positioning system (GPS) measurements obtained from continuous or campaign-style surveying allow the definition of surface deformation patterns at both regional and local scales.The availability of a continuous GPS station (CGPS) network provides an accurate estimation of the strain-rate field, although the actual internal network spatial resolution is not always adequate for the definition of a dense strain-rate pattern in highly complex tectonic settings [Gasparini et al. 1985, Galadini and Galli 2000, D'Agostino et al. 2001, Anzidei et al. 2005].A significant improvement in the mean CGPS network spacing is expensive in terms of human and financial resources, which are often not available.To increase the geodetic observation in the Italian territory, several research institutes and universities are working on the installation of new CGPSs and several non-permanent GPS stations (NPSs).A NPS is a stable structure on which a GPS antenna is mounted for short periods of time, which can range from a few days to a few months per year, depending on the campaign organization.
The observations arising from both NPS and CGPS networks are currently used to investigate geological and geophysical phenomena at a regional scale, e.g. to study deformation caused by tectonic movements [Gomez et al. 2007, Cenni et al. 2008, D'Agostino et al. 2008, Baldi et al. 2009], or at a local scale, e.g. to monitor gravitational macroscopic effects due to rock-mass collapse, landslide activation, or other instabilities [Pesci et al. 2005, Pesci and Teza 2007, Baldi et al. 2008].
The daily coordinates from NPSs and CGPSs processed using the currently available scientific packages that are optimized for CGPS data are generally characterized by an error of a few millimeters for the horizontal components, and about twice that for the vertical coordinates.The station velocity is generally computed using sequential least-squares adjustment procedures, and are also quantified taking into account white noise and time-correlated noise [Langbein 2008].The data processing procedures based on the hypothesis of equi-spaced coordinate time series can model the noise of CGPS, but are not adequate for NPS.Therefore, the velocity error estimates obtained from the standard processing of NPS data can be too optimistic.Also, the effects of irregularly repeated antenna installations, which often result in significant point-position instability [Leonard et al. 2007], are neglected in standard GPS data processing.Pesci el al. [2009] showed a method that was intended to provide reliable estimates of velocity uncertainties of a NPS using subsampling of data related to a nearby CGPS, taking into account the NPS occupation plan, i.e. the set of start and stop survey dates.In this way, the quantities that are

Article history
Received November 18, 2009;accepted May 13, 2010.Subject classification: GPS, Non-permanent station, Subsampling, Velocity, Strain field critically dependent on velocity uncertainties can be reliably computed, e.g. the strain.The algorithmic aspects of this method, and the package that implements it, are described in Teza et al. [in press].
The present study shows the results from a further analysis of data belonging to the Central Apennine Geodetic Network [CAGeoNet, Italy;Anzidei et al. 2005] using a CGPS subsampling-based method.Moreover, a new statistical approach is used to obtain more stable solutions, and therefore to allow stronger interpretation of the results.Finally, emphasis has been placed on interpretation of the data in light of the recent Mw 6.3 L'Aquila earthquake (April 6, 2009).

NPS simulation from CGPS data
The velocities of NPS and CGPS are usually estimated through interpolation of the daily or weekly position time series, computed within a suitable reference frame (e.g.ITRF 2005), and using methods mainly based on least-squares adjustments or similar procedures.In particular, the velocities of the studied sites are generally computed from daily position estimates using the forward Kalman filter ability of GLOBK, or of other GPS analysis packages, with the combining of the NPS and CGPS observations [McClusky et al. 2000, McClusky et al. 2003, Battaglia et al. 2004, Bendick et al. 2006, Casula et al. 2007, Aktug et al. 2009, D'Agostino et al. 2008, Alchalbi et al. 2010, Buble et al. 2010, Weber et al. 2010].
The GPS coordinate time series are affected by significant colored noise, i.e. time-correlated noise, which includes a large predominance of flicker noise.Ray et al. [2008] showed that a continuum flicker noise, i.e. noise with a spectrum of the form 1/f (where f is the frequency), describes the background spectrum down to periods of a few months.Moreover, if longer time-spans are considered, the signal-to-noise ratio is not improved.As indicated by previous studies [Langbein and Johnson 1997, Dixon et al. 2000, Williams et al. 2004, King and Williams 2009], flicker noise components of the colored noise are commonly observed in a wide variety of environmental phenomena, e.g.atmospheric loading, sunspot variations, wobble of the Earth axis, and undersea ocean currents.However, they also appear as bias in estimations of satellite atomic clocks, and in semidiurnal, diurnal, annual and semi-annual frequency bands.Moreover, random-walk noise, the spectrum of which has the form 1/f 2 , has also been seen in long GPS time series.This is probably related to the instability of the structure [Mao et al. 1999].These effects are observable in long time series, but cannot be detected if NPS data are used.
A numerical experiment was performed to evaluate the effects of white, flicker and random-walk noise in velocity estimations from NPS datasets.In particular, three synthetic 20-year-long time series were generated.Let n w (t), n f (t) and n r (t) be the normalized white, flicker and random-walk noise, respectively (the distribution of a normalized noise is 0-centered in the range [−1, 1] m).Three 20-year-long simulated daily time series (7,306 days) were generated: (i) a 12 mm/y linear trend s LT (t) with 3 mm amplitude white noise, i.e. s LT (t) + 0.003 n w (t), as 3 mm is a typical noise level for real observations [see, e.g.Cenni et al. 2008]; (ii) the signal s LT (t) + 0.003 [n w (t) + n f (t)] obtained by adding 3-mmamplitude flicker noise to the previous signal; and (iii) the time series s LT (t) + 0.003 [n w (t) + n f (t) + n r (t)].The power spectra of the added noise were deterministic, with phases uniformly distributed in the range [−r, r].It should be noted that some hypotheses about the phase are needed to obtain a time-domain signal from its power spectrum.Indeed, no information about the phase is provided by the power spectrum [Little et al. 2007].
The continuous time series were subsampled to simulate 2 weeks of yearly non-permanent observations, and therefore with 280 simulated acquisition days over 20 years.For each year, the data were extracted corresponding to the same time periods, according to the typical occupation plan of real campaigns, which are quasi-periodical.Moreover, the procedure was repeated 10 times, with the extraction of data on similar epochs, although with time-shifting to enlarge the velocity dataset, making the next statistical analysis more stable.For each simulation, the estimated velocity was compared to the reference velocity obtained using the whole dataset (the best velocity), which led to a set of velocity differences.The root mean square (rms) velocity differences were computed for each synthetic time series, with respect to the number of simulated campaigns; the results are shown in Figure 1.
It is clear that the colored noise affects the velocity estimates, and that more measurement campaigns are needed to make the rms on small values stable.The noise cannot be modeled if time series that are too short are considered.Therefore, to obtain reliable quantities of geodetic interest from the estimated velocities (e.g. the strain-rate field), the velocity uncertainties must be analyzed appropriately.
Using CGPS or NPS data changes the results significantly in terms of velocity estimates, as shown by applying the described subsampling approach to AQUI CGPS, a station in the central Apennines and in the middle of the CAGeoNet network (Table 1).This area is characterized by high seismic risk, and it was recently hit by the Mw 6.3 April 6, 2010, destructive earthquake that caused extensive damage to buildings and resulted in hundreds of human casualties [Anzidei et al. 2009].
Several simulations were carried out with different settings to study significant cases.The first simulation series was performed to explore the effects of campaign frequency.In particular, the cases of 1, 2 and 3 campaigns per year were considered.The campaign length was 14 days in all the cases.Also in this case, 10 simulations were generated, which led to the velocity difference rms shown in Figure 2.
The second simulation series was performed by varying the campaign duration.In particular, yearly campaign series of 3, 7, 14 and 30 days in length were simulated.As expected, the rms decreases with increasing campaign number (Figure 3).All these results show that as at least five campaigns are performed, the velocity rms lies below 0.5 mm/y, and it rapidly decreases with campaign number until negligible values are reached.In particular, two 7-day sessions per year provide 4 y results that are practically the same as those obtainable with one 30-day measurement session per year, since the difference is 0.5 mm/y (absolute value), even if the corresponding instrument days are 56 and 120, respectively.
At present, each NPS is generally equipped with a highaccuracy (sub-millimeter) structure for the GPS antenna.Nevertheless, even if the antenna installation/removal operations are carefully performed, some operator-induced errors can affect the results.These can depend on several factors, such as level calibration, a different antenna (even if  the antenna model is the same, there can still be small differences), and other factors.The third simulation series was performed by imposing the effects of antenna installation/removal operations for campaign-style measurement simulations.A synthetic operator-induced error, i.e. a random distribution of horizontal steps of ~1-2 millimeters was imposed on the AQUI station coordinates (positioning offsets), and the above analysis was repeated.Yearly campaigns were extracted with 7-day durations, and 10 simulations were performed that considered different epochs for data extraction (Figure 4).The rms of the two horizontal components goes to zero, increasing the number of occupations.A similar behavior appears if the antenna replacing errors are taking into account, even if the convergence to zero is slow.Also, the effects of the stepwise error for the vertical component are similar, although wider amplitudes are involved.
These results show that NPS data require careful analysis and that short coordinate time series that are well distributed in time are a lot better than relatively long time series with too long a returning time.Moreover, the millimeter uncertainty of NPS data estimated using standard GPS analysis software appears to be highly underestimated.Therefore, a valid method for uncertainty evaluation is needed.
A possible solution for the estimation of reliable velocity errors of a NPS was shown in Pesci et al. [2009], where a method based on subsampling of data provided by one or more CGPSs neighboring this NPS was described in detail, together with discussion and validation.This method was mainly based on the hypothesis that similar noise sources should affect the data provided by GPS stations inside a welllocalized area (e.g.10-50 km), depending on the size of the GPS network.If GPS stations are installed in areas characterized by high topographic heterogeneity, the atmospheric effects might not be completely modeled and removed, which would lead to degradation of the results in terms of coordinate time series obtained by standard operations of data processing [Wdowinski et al. 1997, Zhang et al. 1997, Dong et al. 2002].Actually, there is no unambiguous definition of the characteristic distances between the reference CGPSs and the NPSs that will allow modeling of the NPS behavior from the CGPS data.However, a reasonable solution appears to be the definition of the area containing the NPSs used and the choice of all of the available CGPSs placed inside and around the area as reference stations for the velocity uncertainty modeling.Clearly, the smaller the area is, the more representative the neighboring CGPSs are expected to be [King and Watson 2010].For a specified real NPS, the method follows three steps: 1) The search for one or more neighboring CGPSs where the time series can be modeled as linear trends plus noise, and the computation of the corresponding reference velocities and related standard deviations.It must also be noted that the CGPS time series used are already postprocessed data, that have been correctly analyzed, checked and corrected of evident anomalies (e.g.steps and outliers caused by earthquake or human damage, antenna changes, electronic problems).
2) For each CGPS, performing several simulations via CGPS data subsampling, taking into account the defined occupation plan and antenna replacement effects.If the occupation plan consists of n S surveys where the start and stop dates and corresponding duration are t Im , t Fm (t Fm = t Im + t Dm ) and t Dm , respectively (m = 1, 2,..., n S ), the simulated m-th survey of each simulation series has a start date randomly obtained from (or, alternatively, regularly shifted from) t Im and duration t Dm .In this way, for each number of simulated occupations, a set of NPS velocities is obtained.
3) The statistical analysis of the distribution of differences between the reference velocities and the velocities of the simulated NPSs.For each number of simulated occupations (on the basis of the specific occupation plan), the rms of these velocity differences can be considered as thresholds of the real NPS velocity uncertainty.If for a velocity component of the NPS, the standard GPS postprocessing provides an error smaller than this threshold, the threshold should be used as an estimate of the velocity error.Clearly, to have meaningful results from this subsamplingbased approach, the requirement that the standard deviation of the reference velocity is significantly smaller than the threshold obtained must be satisfied.
The results can be used to both analyze the available data (the aim of this study) and optimize the survey plans according to the accuracy needed and the human and financial resources available.

Velocity pattern and strain-rate estimation
The CAGeoNet was designed to measure both the subregional and near-field strain rates across the main seismogenetic structures and faults, which are believed to  respond to the crustal dynamics of the area.This network comprises NPSs measured for the first time during a session over a few days in the 1999-2001 time-span.The first estimation of the strain rate of the area was computed over a regular 10 km × 10 km grid using velocity data processed over three years.These preliminary results showed both extensional and compressive behavior at the same values along the chain axis, ranging from about −10 −7 y −1 to +10 −7 y −1 [Anzidei et al. 2005].
In this section, the main results of the processing and related analyses of the CAGeoNet data obtained for the timespan of 1999-2007 are presented.The GPS data were processed using the GAMIT 10.3 software [Herring et al. 2006a], and included 22 permanent stations located over the whole Italian territory.The registration in the IGS05 external reference frame was performed using about 26 International Global Navigation Satellite Systems (GNSS) Service (IGS) stations located both in Italy and Europe [Bruyninx 2004, Kenyeres andBruyninx 2004].
The IGS stations of the Global Permanent Network included in the regional computation were used, together with positions and velocities from the Scripps Orbit and Permanent Array Center (SOPAC) facility ftp site [SOPAC 2009].Robust combinations of our regional solutions with the global ones computed by SOPAC for the IGS GPS network were obtained using the GLOBK package [Herring et al. 2006b].A distributed session approach was applied that was similar to sequential least-squares minimization, using sequential Kalman filtering and considering the | 2 parameters to estimate the levels of data fitting [Dong et al. 1998, Casula et al. 2007].Finally, the total adjusted solution was obtained, in terms of the coordinate time series and velocities in the same reference frame.The inherently rigid rotation of the Eurasian plate was removed by applying the general definition of the Euler pole given by Altamimi in the ITRF2005 reference frame [Altamimi et al. 2007], and using the plate module of the GLOBK package.
The resulting intra-plate velocity vectors that ranged between zero and a few millimeters per year are listed in Table 1, and the residual velocities with respect to a fixed Eurasian frame [Altamimi et al. 2007] are shown in Figure 6.The velocity uncertainties were estimated using the statistical method proposed by Williams et al. [2004].Moreover, the time series were iteratively edited by means of the tsview MATLAB toolbox [Herring 2003], to estimate/remove outliers and offsets, and subsequently verified by means of the CATS package [Williams 2008].Taking into account the non-permanent nature of the measurements, the short session lengths (about 2-3 days), and the results provided by previous simulations, the statistical values obtained appeared to underestimate the errors.
The CGPS subsampling-based approach presented by Pesci et al. [2009] was applied to the time series of five CGPSs located inside and around the area investigated, i.e.AQUI, INGR, RSTO, UNPG and VVLO.The typical measurement  sessions of a CAGeoNet NPS characterized by 2-3 day duration were considered.Unfortunately, this was not the best choice for the stability of the velocity estimation, as shown in the previous section, but these were the only data available.Moreover, the campaigns were typically irregularly spaced in time.For this reason, for each real NPS, the uncertainty analyses were carried out on the basis of the occupation plan; in order to reduce the number of calculations, these NPSs were grouped with a criterion of similarity in the occupation plan.The thresholds of velocity errors were therefore obtained for each campaign number and for an assigned occupation plan.Moreover, Table 2 shows the results in the special case in which almost regularly spaced yearly campaigns are simulated in the whole network, to show how the errors change with respect to those provided by standard CGPS processing applied to the NPS time series.As shown by Pesci et al. [2009], the procedure is based on a regular time sequence of campaigns that provides the worst case for timemeasurement distribution.For this reason, the threshold cannot be underestimated.
The velocity uncertainties evaluated using the subsampling-based thresholding method are shown in Figure 5 (blue), where it can be seen that the statistical values obtained by standard data processing were not significantly underestimated if more than three campaigns were available.Clearly, if a subsampling-based threshold does not exceed the velocity error provided by the standard processing software, this threshold is not used in further computations.
The velocity pattern and the corresponding uncertainties estimated using both the standard data processing, where the NPSs are processed like CGPSs, and the subsampling-based thresholding method are shown in Figure 6, where it can be seen that in the latter case, the errors are increased.These velocity patterns were used to compute the strain rate on the nodes of a 10 km × 10 km regular grid distributed in the area.
A modified least-squares method was used, with rescaling of the covariance matrix of velocity data (observables), and considering a weighting function that took into account the distances between the grid node and the GPS stations [Shen et al. 1996].This method was implemented in the grid_strain package [Teza et al. 2008], and together with the estimation of the strain rate and the corresponding uncertainty, it also provides a geometrical significance to the results.The strain rate in a grid node is highly significant only if the stations are well distributed all around this node (requirement for angular distribution) and if their distances do not exceed about three times the mean inter-station distance (requirement for radial distribution).For each NPS station, the errors associated to the velocity solutions in the 1999-2007 time-span were checked on the basis of the thresholds obtained as a function of the number of end epochs of the campaigns performed.
The scale factor considered in modified least-squares computations was 25 km, according to the mean spacing between the GPS stations.The deformation appears to be in agreement with the estimations by Anzidei et al. [2005].The computed strain rates ranged from 1.2 × 10 −8 y −1 to 1.6 × 10 −8 y −1 and from −1.4 × 10 −8 y −1 to −0.3 × 10 −8 y −1 along the normal and parallel directions to the chain axis, respectively, with a 1.1 × 10 −8 y −1 mean error.
Figure 7 shows the effects of the thresholding-based method on computation of the principal strains.In particular, the maximum of the strain-rate tensor field computed from the standard and threshold-modified velocity fields are shown together with their differences, as also for the minimum eigenvalues.These differences are not zero centered, which indicates that the deformation tensor is globally changed as an effect of error thresholding, even if the values are small with respect to the mean computation error of the order of 10 −8 .

Regional velocity field
The velocity field with the subsampling-based error thresholding (Figure 6d) shows a high heterogeneity for both vector modules and directions.The CAGeoNet was designed to measure the regional and near-field strain rates across the main seismogenetic structures and faults.Thus, the aim was to highlight possible differences in the regional and subregional velocity fields.For this reason, a simple statistical method was applied to highlight areas with similar kinematics.The method was based on an iterative search in a confined area, to compare the direction and module of the GPS station velocity.Simple criteria were used to select stations with similar trends, in terms of direction and module.
For the agreement level between the velocity directions of two stations, a scalar parameter, D ij , is defined by the relationship: (1) where are the normalized East and North velocity components of the i and j stations, respectively.The discrepancy parameter D ij ranges from 0 (perfect agreement) to 1 (total disagreement).Figure 7 shows the D ij parameter as a function of the angular differences between the directions of two normalized vectors.In particular, the reference vector V i is fixed (unit vector of y-axis, i.e. where T indicates the vector transposition), while the vector V j is rotated, providing all of the possible directions in the [0,360˚] interval.Principal strain rates computed using the standard velocity patterns and the threshold-modified patterns.The differences between the modified and original eigenvalues (thick black lines) have means of 0.8 × 10 −8 y −1 (standard deviation 0.5 × 10 −8 y −1 ) for the maximum eigenvalues (left panel), and means of 0.4 × 10 −8 y −1 (standard deviation 0.3 × 10 −8 y −1 ) for the minimum strain (right panel).
The discrepancy between the velocity modules of the station are also estimated by comparing the single modules: (2) To select the stations characterized by similar trends, for each i-station, a circular area of radius r (search radius) is defined, and the set of the j-stations inside this is recognized.The velocity of the i-station is then compared to all of the jstation velocities, defining the D ij and D ij parameter dataset.Subsequently, to obtain a unique pair of values representative of the discrepancy level for each i-station, the mean and the rms of the D ij and D ij distributions are considered.The parameters D i and D i are defined as follows, in order to characterize the angular and modular distributions by means of a single value, while taking into account both the mean and the rms: The computation was performed using the total dataset, and varying the search radius from very local (10 km) to regional (90 km), as shown in Figure 9.
Figure 10 shows the D m and D m that are representative of the GPS in the whole network with respect to the radius (10-90 km).The lower value that corresponds to the 10 km radius indicates a more scattered data distribution.In particular, changes in the discrepancies parameters occur going from 10 km to 40 km.Longer ranges provide values of about 0.33 (D i distributions) and 0.97 (D i distributions).1.
This indicates that the sizes of the sub-areas affected by different kinematics are probably lower than 40 km.This statistical method characterizes velocity discrepancies, which depend on the chosen search-radius and threshold values for D i and D i , and it is aimed at extracting some kinematicsrelated information, and in particular to recognize possible different velocity patterns in the area, without the use of a tectonic or a geodynamic model.In this way, a-priori hypotheses that might force the results are avoided.
The discrepancy parameters shown above that characterize the network were used together with a search radius of 20 km, to select a reduced velocity dataset that was representative of a regional pattern (Figure 11, left panel).The results are superimposed on a map that shows the topographical information of the area, and they provide a detailed view of the 1999-2007 sub-regional and near-field velocity patterns of this region.Figure 11 right panel shows the complete velocity pattern and in both cases the strain rate tensors are shown.

Discussion
Although a geological interpretation of these data is not within the scope of the present study, which mainly concerns the methodological procedures for data treatment, some aspects can be observed and discussed (see Figure 6d and Figure 10).It has to be noted that a large portion of the velocity vectors are characterized by errors at the same level as the velocity magnitudes, which makes unambiguous data interpretation very difficult.The eastern part of the network, and in particular the Gran Sasso area, presents a NNE trend, which is in agreement with the extrusion pattern of this sector of the belt that has been suggested by some studies [e.g.Mantovani et al. 2006], and with the activity of the Assergi and Campo Imperatore faults [Anzidei et al. 2005].The sites located on the Umbria-Marche Apennines (UMA) and in the north-western sector of the Fucino plain might be compatible with the sinistral transtensional tectonics observed in the area [Galli et al. 2008], although more information cannot be obtained using the data processing proposed.A N-NE trend is seen in the north-eastern part of the Fucino Valley, although also in this case the results are relatively limited by the errors.The more significant station velocities in the Liri Valley (south-western part of the Fucino plain) are characterized by a N-NE-ward trend, and they are responsible for the compressive deformation pattern computed in the Fucino plain (Figure 6d and Figure 11).This result is clearly in disagreement with the extensional pattern expected in a basin area, but this is only a preliminary analysis, where the necessary geological and geophysical information about the local movement for each site is neglected.
Some stations that belong to CaGeoNet that were included in the data analysis are in zones affected by large gravitational effects that can lead to movements that are not related to tectonics.Indeed, a preliminary study that was also based on an analysis of photogrammetric images reported the stations affected by local effects, e.g.SECI, ACCI and CVSE [Galvani 2009].
Despite this partial disagreement related to local movements, the general results for the strain rate that were obtained using this thresholding-based approach agree substantially with the known regional kinematics, whereas the strain-rate field computed without the application of this method does not.This confirms that the standard GPS data processing optimized for CGPS data treatment cannot provide reliable velocity uncertainty estimates if NPS data are used.
Clearly, another important result of the velocity uncertainty analysis is that new surveys need to be planned and performed in the coming years, to improve the velocity, and therefore the strain-rate, patterns.
As well as the thresholding of velocity uncertainties, a simple but efficient method was developed and applied that was intended to recognize areas that are characterized by similar velocity vectors.In this way, the subareas that are kinematically homogeneous at the regional scale were easily recognized, and in each subarea the displacement vector is shown.On the other hand, the border between two subareas should be characterized by an abrupt spatial change in local strain rate, whereas possible changes in strain rate within an area should be relatively low.Figure 11 shows the strain-rate field from selected velocities and that from the whole errormodified velocity field, which allows a direct visual comparison.The first characterizes a regional trend, and therefore provides global information.This shows two main zones that are positively and negatively deformed.The second provides relatively local information because all of the velocity vectors are used.Both of these results need to be taken into account in the data interpretation.Figure 12 summarizes the main result using a contouring map of the computed strain rate from the threshold-corrected velocity field (using the whole velocity dataset according to Figure 6d).The rate of change in an area, i.e. the sum of the principal strain rates, is used to recognize possible lines that are characterized by high-strain gradients, which in extreme RELIABLE DEFORMATION RATE FROM GPS NPS DATA  cases might be related to the existence of faults or other discontinuities.
The data used in this study were derived from nonpermanent GPS campaigns that were performed from 1999 to 2007.Additional information relating to the displacements caused by the Mw 6.3 L'Aquila earthquake (April 6, 2009) was obtained by processing CGPS data from the weeks before and after these shocks, and this is superimposed on the map (Figure 12).The AQUI and INGP stations, which were affected by the largest displacements, are located in the area that is mainly characterized by a different strain-rate pattern, which coincides with the line of the Apennine chain.Moreover, the results from Anzidei et al. [2009] are included in Figure 12, providing further information.In Figure 12, the area in which the strain-rate computation is characterized by a high significance level is highlighted.

Conclusions
A CGPS subsampling-based method was applied to the CAGeoNet data to obtain reliable errors in NPS velocity estimates, and therefore to allow computation of a reliable strain-rate field at a regional scale.The results show that: (a) the strain rate computed using this method can account for the existence of seismogenetic features, and it appears to be more regular than the strain-rate field computed without the specific treatment; (b) only velocities obtained from more than three measurement campaigns that were characterized by session lengths of a few days led to errors lower than 1 mm/y; (c) the method provides operative suggestions for the design of a NPS network and related observation planning, and for the inclusion of NPS velocities in the tectonic study.
Moreover, a simple statistical method was developed that was intended to detect possible regular patterns in the velocity field related to a network with a large number of stations.Although this procedure is strongly user-dependent because of the thresholds imposed on the module for the direction variations and search radius, its application to the CaGeoNet network allowed the robust recognition of some areas that are characterized by different kinematic patterns.In this way, these local-scale and regional-scale behaviors can be studied and compared.
These data agree with those from previous studies, and they provide very useful information about kinematics at different scales.Moreover, the displacements due to the L'Aquila earthquake of April 2009 are in areas that are characterized by the main strain-rate variations, suggesting that reliably corrected NPS data can be effectively used for territory monitoring.Finally, these procedures indicate that the currently available information from CAGeoNet is still too poor to provide reliable constraints to geological and geodynamical interpretations.

Figure 1 .
Figure 1.Empirical relationship between the rms of the velocity differences and the campaign numbers using synthetic data from three signals: white noise; white + flicker noise; white + flicker + random-walk noise.

Figure 3 .
Figure 3. Mean rms of the velocity differences for the east and north components based on the campaign numbers.

Figure 2 .
Figure 2. Rms of the velocity differences obtained by simulating 14-day campaigns characterized by different frequencies (campaigns per year).Bottom right panel: summary of the data showing the means of east and north velocity rms.

Figure 4 .
Figure 4. Mean velocity rms from the simulation of 7-day yearly campaigns with antenna installation errors (antenna offsets).

Figure 5 .
Figure5.Errors of the east and north velocity components related to the CAGeoNet stations, including both the standard and the modified errors corrected with thresholds.Since the NPSs measured using three campaigns show standard errors lower than the computed thresholds, their errors are replaced by the corresponding thresholds.In contrast, NPSs measured more than three times are characterized by errors that overcome the thresholds, and their final errors (blue) coincide with the standard errors.The station numbers along the abscissa are as listed in Table1.The ordinate axes are labeled as «V error» to highlight the corrected errors dataset (blue line).

Figure 6 .
Figure 6.Velocity patterns and associated 95% confidence error ellipses obtained with standard (a) and threshold-modified (b) uncertainties.Strain-rate fields were computed using the same velocity patterns, but different associated errors, as standard (c) and threshold-modified (d).

Figure 9 .
Figure 9. Distributions of D i and D i computed on different search radii from 10 km to 90 km.Gray shading: number of j-stations around each i-stations based on the search radii.For station numbers, seeTable 1.

Figure 8 .
Figure 8.The D ij parameter shows the discrepancy level between vector V i and V j direction in the [−1, 1] interval (red squares).Abscissa: angle between the reference velocity V i = [0 1] T and V j .The vector components, i.e. cosine and sine of the angle, are shown by solid and dashed lines respectively.

Figure 10 .
Figure 10.The D m and D m values computed as the means of the D i and D i parameters shown in Figure 9, for each search radius.

Figure 11 .
Figure11.Left panel: velocity field obtained using the statistical method proposed, with a search radius of 20 km and the matching parameters D i = 0.33 mm/y and D i = 1.0 mm/y.Cyan blocks, sites with similar trends.The strain rate field estimated using a modified least-squares approach is also shown.Right panel: previous velocity fields (bold black arrows) and trends (blue arrows) of the stations ruled out of the analysis, together with modified errors.Umbria-Marche Apennines (UMA) region.

Figure 12 .
Figure 12.Change-in-area contour map from the 1999-2007 corrected velocity fields and the displacements caused by April 6, 2009, destructive earthquake.Computed velocities (cm/y; blue arrows) and data published by Anzidei et al. [2009] are also shown (gray arrows).

Table 1 .
Velocity pattern and velocity errors from standard (GLOBK-based) data processing, including station names, station coordinates (latitude and longitude), occupations, velocities and velocity errors.If CGPS stations are considered, the label «perm» is used instead of the integer number of the relevant campaign.

Table 2 .
East and north velocity error thresholds, defined from the subsampling procedure.The velocity thresholds for the relevant campaign numbers are given.