Layout 6 Non-permanent GPS data for regional-scale kinematics: reliable deformation rate before the 6 April, 2009, earthquake in the L'Aquila area Arianna Pesci1,*, Giordano Teza2, Giuseppe Casula1, Nicola Cenni3, Fabiana Loddo1 1 Istituto Nazionale di Geofisica e Vulcanologia, sezione di Bologna, Italy 2 Dipartimento di Geoscienze, Università degli Studi di Padova, Italy 3 Dipartimento di Scienze della Terra, Università degli Studi di Siena, Italy ANNALS OF GEOPHYSICS, 53, 2, APRIL 2010; doi: 10.4401/ag-4740 ABSTRACT 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 non- permanent 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 55 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 nw (t), nf (t) and nr (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 sLT (t) with 3 mm amplitude white noise, i.e. sLT (t)+ 0.003 nw (t), as 3 mm is a typical noise level for real observations [see, e.g. Cenni et al. 2008]; (ii) the signal sLT (t) + 0.003 [nw (t) + nf (t)] obtained by adding 3-mm- amplitude flicker noise to the previous signal; and (iii) the time series sLT (t) + 0.003 [nw (t) + nf (t) + nr (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 RELIABLE DEFORMATION RATE FROM GPS NPS DATA 56 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. 57 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 high- accuracy (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 RELIABLE DEFORMATION RATE FROM GPS NPS DATA Figure 3. Mean rms of the velocity differences for the east and north components based on the campaign numbers. 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. 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 well- localized 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 nounambiguous 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 post- processed 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 nS surveys where the start and stop dates and corresponding duration are tIm, tFm (tFm= tIm + tDm) and tDm, respectively (m = 1, 2,..., nS), the simulated m-th survey of each simulation series has a start date randomly obtained from (or, alternatively, regularly shifted from) tIm and duration tDm. 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 post- processing 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 subsampling- based 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 sub- regional and near-field strain rates across the main seismogenetic structures and faults, which are believed to RELIABLE DEFORMATION RATE FROM GPS NPS DATA 58 Figure 4.Mean velocity rms from the simulation of 7-day yearly campaigns with antenna installation errors (antenna offsets). 59 RELIABLE DEFORMATION RATE FROM GPS NPS DATA T ab le 1 .V el oc it y pa tt er n an d ve lo ci ty e rr or s fr om s ta nd ar d (G L O B K -b as ed ) d at a pr oc es si ng , i nc lu di ng s ta ti on n am es , s ta ti on c oo rd in at es (l at it u de a nd lo ng it u de ), o cc u pa ti on s, v el oc it ie s an d ve lo ci ty e rr or s. If C G P S st at io ns a re c on si de re d, t he la be l « pe rm » is u se d in st ea d of t he in te ge r nu m be r of t he r el ev an t ca m pa ig n. N um N am e La t ( °) Lo n (° ) C am p V e (m m /y ) V n (m m /y ) !! V e (m m /y ) !! V n (m m /y ) 1 A C C I 13 .7 44 42 .1 78 3 -1 .2 3. 4 1. 0 1. 3 2 A N A T 13 .5 49 42 .1 82 3 0. 7 0. 9 0. 6 0. 7 3 B A N O 13 .5 82 42 .3 37 3 1. 0 1. 7 0. 7 0. 6 4 B S P I 13 .6 50 42 .3 06 3 1. 3 1. 4 0. 8 0. 7 5 B S S I 13 .8 15 42 .2 14 3 0. 8 3. 2 0. 9 0. 8 6 C D A Q 12 .9 91 42 .5 90 3 2. 7 2. 7 1. 1 1. 1 7 C ID A 12 .9 07 42 .5 16 3 1. 1 3. 4 0. 8 0. 6 8 C IN C 13 .4 05 42 .0 08 3 -0 .1 -0 .1 0. 6 0. 6 9 C N T A 12 .9 16 42 .4 76 3 1. 3 0. 6 1. 2 1. 3 10 C O C A 13 .3 75 42 .3 02 3 0. 9 2. 1 0. 8 0. 8 11 C P P D 13 .2 69 41 .9 98 3 -0 .3 1. 7 0. 5 0. 4 12 C P S E 12 .9 07 42 .6 04 3 0. 1 2. 8 0. 7 1. 0 13 C U M U 13 .0 28 42 .5 69 3 -0 .9 1. 5 0. 7 1. 0 14 C V A L 13 .8 11 41 .9 84 3 -0 .2 3. 3 0. 5 0. 5 15 C V S E 13 .7 45 42 .1 31 3 -0 .6 3. 8 0. 5 0. 6 16 F C LM 13 .4 59 42 .1 11 3 -0 .4 2. 5 0. 6 0. 7 17 F R C A 13 .6 78 42 .0 59 3 1. 3 1. 1 0. 7 0. 5 18 F R LC 12 .9 75 42 .6 96 3 -1 .6 3. 2 1. 5 1. 7 19 F T T I 13 .3 34 41 .8 92 3 -0 .2 0. 9 0. 7 0. 6 20 G R E C 12 .7 31 42 .4 45 3 -0 .8 -0 .1 0. 8 0. 6 21 M A M A 13 .2 53 42 .2 09 3 -1 .1 -0 .6 0. 9 0. 8 22 M E R 2 13 .0 18 42 .6 87 3 -0 .9 3. 8 0. 9 1. 0 23 M LN N 13 .7 05 41 .8 22 3 -1 .1 1. 8 0. 9 0. 8 24 M P E T 12 .8 33 42 .6 00 3 -0 .1 1. 3 0. 8 0. 6 25 M R P N 13 .6 85 41 .8 86 3 -0 .7 1. 7 0. 6 0. 8 26 M S LV 13 .7 42 42 .0 51 3 1. 9 0. 3 0. 6 0. 5 27 M V IP 13 .3 77 41 .9 29 3 1. 1 1. 9 0. 6 0. 6 28 N U R I 13 .0 61 42 .3 75 3 5. 0 -2 .4 1. 0 0. 8 N um N am e La t ( °) Lo n (° ) C am p V e (m m /y ) V n (m m /y ) !! V e (m m /y ) !! V n (m m /y ) 29 O C R E 12 .9 71 42 .5 97 3 -2 .3 -0 .4 1. 6 1. 5 30 O R T O 13 .7 25 41 .9 99 3 -1 .2 3. 4 0. 7 0. 6 31 P O G B 12 .8 73 42 .5 15 3 0. 5 -0 .1 1. 0 1. 1 32 P P E Z 13 .4 26 42 .1 83 3 0. 8 2. 0 0. 8 0. 7 33 P R C A 13 .4 94 41 .7 82 3 1. 1 2. 7 0. 6 0. 6 34 P S C A 13 .1 25 42 .1 28 3 -1 .2 1. 8 0. 8 0. 6 35 P S IR 13 .5 93 42 .1 80 3 -1 .1 1. 4 0. 7 0. 7 36 P S M A 13 .5 81 42 .1 27 3 1. 1 1. 6 0. 7 0. 6 37 Q R T O 13 .3 99 42 .2 84 3 -1 .3 0. 4 0. 8 0. 6 38 R D C A 13 .4 89 42 .2 28 3 2. 3 0. 5 0. 6 0. 9 39 R E N D 13 .4 65 41 .8 34 3 1. 1 3. 6 0. 8 0. 6 40 R ID O 13 .6 05 41 .7 99 3 0. 8 2. 5 1. 0 1. 0 41 R IF P 13 .1 76 42 .7 63 3 0. 5 3. 2 1. 3 1. 4 42 R O C A 13 .6 97 42 .3 28 3 0. 9 2. 2 0. 5 0. 5 43 S E C I 13 .6 70 42 .1 48 3 0. 5 1. 5 0. 8 0. 7 44 S E T C 12 .7 96 42 .4 68 3 0. 5 -0 .2 1. 2 1. 4 45 S IE R 13 .6 68 41 .9 25 3 0. 0 1. 3 0. 5 0. 5 46 S LL I 13 .1 21 42 .7 27 3 0. 3 2. 7 1. 4 1. 0 47 S O R B 13 .3 17 42 .0 82 3 -0 .9 1. 5 0. 6 0. 6 48 S S 83 13 .7 80 41 .8 42 3 -0 .2 1. 0 0. 7 0. 8 49 S S E B 13 .7 52 41 .9 34 3 -1 .3 1. 0 0. 7 0. 5 50 S S T S 13 .6 51 42 .3 60 3 1. 3 1. 0 0. 5 0. 8 51 T E R R 12 .7 79 42 .4 26 3 -1 .1 4. 8 0. 9 1. 2 52 T N E R 13 .5 32 42 .2 37 3 -1 .1 -0 .6 0. 8 0. 7 53 T R A S 13 .5 43 41 .9 54 3 0. 1 2. 3 0. 5 0. 7 54 T R M T 13 .2 01 42 .0 96 3 -0 .5 2. 3 0. 6 0. 6 55 T V IO 12 .9 88 42 .6 56 3 1. 8 1. 5 1. 6 1. 0 56 V C A R 12 .9 96 42 .5 66 3 -3 .2 0. 2 1. 8 2. 4 N um N am e La t ( °) Lo n (° ) C am p V e (m m /y ) V n (m m /y ) !! V e (m m /y ) !! V n (m m /y ) 57 V LA G 13 .8 29 41 .9 38 3 0. 9 -0 .5 1. 0 0. 9 58 V R C E 13 .2 40 42 .0 39 3 0. 1 1. 6 0. 8 0. 4 59 V T F U 12 .9 18 42 .5 64 3 2. 0 2. 4 0. 7 0. 8 60 C LA C 12 .9 76 42 .5 46 4 -0 .9 0. 1 0. 6 0. 9 61 P R E T 13 .3 16 42 .3 82 4 -1 .0 -1 .6 1. 5 1. 3 62 B O R B 13 .1 62 42 .5 11 5 0. 1 0. 7 1. 1 1. 1 63 C A S B 12 .8 49 42 .3 90 5 -0 .3 0. 7 1. 3 1. 0 64 C R O G 13 .4 85 42 .5 86 5 0. 2 2. 8 1. 0 1. 2 65 LA C U 13 .1 12 42 .4 88 5 1. 5 0. 0 1. 1 1. 2 66 M O S P 12 .9 49 42 .6 47 5 0. 1 0. 4 0. 8 0. 6 67 M S A N 13 .1 54 42 .7 61 5 -0 .1 2. 7 1. 5 1. 3 68 P O C A 13 .3 26 42 .5 71 5 0. 5 4. 0 1. 2 1. 3 69 R O F A 13 .5 41 42 .3 97 5 -0 .3 1. 0 0. 7 0. 6 70 R O IO 13 .3 86 42 .3 27 5 -0 .1 1. 5 1. 5 0. 9 71 S C U O 13 .3 59 42 .6 29 5 -0 .9 1. 5 1. 3 0. 9 72 S F R A 13 .4 08 42 .4 60 5 1. 0 3. 1 2. 0 1. 5 73 S R O T 13 .1 42 42 .6 28 5 -2 .7 1. 9 1. 2 1. 5 74 V P E Z 13 .4 85 42 .1 82 5 -1 .3 0. 5 1. 0 2. 1 75 C E P P 12 .8 55 42 .5 30 6 -0 .8 1. 5 0. 6 0. 9 76 F V IO 13 .0 87 42 .5 35 6 -0 .1 0. 3 0. 9 1. 6 77 T R LO 13 .0 10 42 .4 69 7 -1 .6 0. 8 1. 1 1. 1 78 A Q U I 13 .3 50 42 .3 68 pe rm 0. 6 2. 0 0. 6 0. 5 79 IN G P 13 .3 16 42 .3 82 pe rm 1. 9 3. 8 0. 6 0. 4 80 IN G R 12 .5 15 41 .8 28 pe rm -0 .2 1. 7 0. 5 0. 5 81 IP R A 12 .7 05 42 .4 84 pe rm 1. 6 0. 2 0. 8 1. 3 82 R S T O 14 .0 01 42 .6 58 pe rm 1. 5 2. 5 0. 5 0. 5 83 U N P G 12 .3 56 43 .1 19 pe rm -1 .0 1. 0 0. 5 0. 6 84 V V LO 13 .6 23 41 .8 70 pe rm -1 .8 2. 5 0. 6 0. 6 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−7y−1 to +10−7y−1 [Anzidei et al. 2005]. In this section, the main results of the processing and related analyses of the CAGeoNet data obtained for the time- span 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 and Bruyninx 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 RELIABLE DEFORMATION RATE FROM GPS NPS DATA 60 Campaigns Ve rms (mm/y) Vn rms (mm/y) 2 7.6 11.7 3 1.8 2.1 4 1.4 1.2 5 0.7 0.7 6 0.5 0.5 7 0.4 0.3 8 0.3 0.4 9 0.2 0.4 Table 2. East and north velocity error thresholds, defined from the subsampling procedure. The velocity thresholds for the relevant campaign numbers are given. Figure 5.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 Table 1. The ordinate axes are labeled as «V error» to highlight the corrected errors dataset (blue line). 61 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 time- measurement 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 RELIABLE DEFORMATION RATE FROM GPS NPS DATA 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). 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−8y−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 sub- regional 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, Dij, 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 Dij ranges from 0 (perfect agreement) to 1 (total disagreement). Figure 7 shows the Dij parameter as a function of the angular differences between the directions of two normalized vectors. In particular, the reference vector Vi is fixed (unit vector of y-axis, i.e. where T indicates the vector transposition), while the vector Vj is rotated, providing all of the possible directions in the [0,360˚] interval. RELIABLE DEFORMATION RATE FROM GPS NPS DATA 62 Figure 7. 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−8y−1 (standard deviation 0.5 × 10−8y−1) for the maximum eigenvalues (left panel), and means of 0.4 × 10−8y−1 (standard deviation 0.3 × 10−8y−1) for the minimum strain (right panel). 2 1 ,D V V V V2 2ij E= +- -iN ijN jE^ ^h h , , andV V V VEiN i jN jE / ;V V V V V( N N E= + =; ;N E /V V V )E N E= +; ; 1 0V ji T= =t 6 @ 63 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 j- station velocities, defining the Dij and Dij 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 Dij and Dij distributions are considered. The parameters Di and Di 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: Di = (Di0+ Di1)/2, Di = (Di0+ Di1)/2. 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 Dm and Dm that are representative of the GPS velocity 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 (Di distributions) and 0.97 (Di distributions). RELIABLE DEFORMATION RATE FROM GPS NPS DATA Figure 9.Distributions of Di and Di 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, see Table 1. Figure 8.The Dij parameter shows the discrepancy level between vector Vi and Vj direction in the [−1, 1] interval (red squares). Abscissa: angle between the reference velocity Vi = [0 1] T and Vj. The vector components, i.e. cosine and sine of the angle, are shown by solid and dashed lines respectively. V V V V2 2ij iE iN E Nj j= + +D -^ ^h h 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 Di and Di, and it is aimed at extracting some kinematics- related 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 RELIABLE DEFORMATION RATE FROM GPS NPS DATA 64 Figure 10. The Dm and Dm values computed as the means of the Di and Di parameters shown in Figure 9, for each search radius. 65 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 error- modified 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 Figure 11. Left panel: velocity field obtained using the statistical method proposed, with a search radius of 20 km and the matching parameters Di = 0.33 mm/y and Di = 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. 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). cases might be related to the existence of faults or other discontinuities. The data used in this study were derived from non- permanent 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. References Aktug, B., J. M. Nocquet, A. Congöz, B. Parsons, Y. Erkan, P. England, O. Lenk, M.A. Gürdal, A. Kilicoglu, H. Akdeniz and A. Tekgül (2009). Deformation of western Turkey from a combination of permanent and campaign GPS data: Limits to block-like behavior, J. Geophys. Res., 114, B10404, doi:10.1029/2008JB006000. Alchalbi A., M. Daoud, F. Gomez, S. McClusky, R. Reilinger, M. Abu Romeyeh, A. Alsouod, R. Yassminh, B. Ballani, R. Darawcheh, R. Sbeinati, Y. Radwan, R. Al Masri, M. Bayerly, R. Al Ghazzi and M. Barazangi (2010). Crustal deformation in northwestern Arabia from GPS measurements in Syria: Slow slip rate along the northern Dead Sea Fault, Geophysical Journal International, 180 (1),125-135. Altamimi, Z., X. Collilieux, J. Legrand, B. Garayt and C. Boucher (2007). ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters, Journal of Geophysical Research, 112, B09401, doi:10.1029/2007JB004949. Anzidei, M., P. Baldi, A. Pesci, A. Esposito, A. Galvani, F. Loddo, P. Cristofoletti, A. Massucci and S. Del Mese (2005). Geodetic deformation across the Central Apennines from GPS data in the time-span 1999-2003, Annals of Geophysics, 48 (2), 259-271. Anzidei, M., E. Boschi, V. Cannelli, R. Devoti, A. Esposito, A. Galvani, D. Melini, G. Pietrantonio, F. Riguzzi, V. Sepe and E. Serpelloni (2009). Coseismic deformation of the destructive April 6, 2009 L'Aquila earthquake (central Italy) from GPS data, Geophysics Research Letters, 36, L17307, doi:10.1029/2009GL039145. Baldi, P., N. Cenni, M. Fabris and A. Zanutta (2008). Kinematics of a landslide derived from archival photogrammetry and GPS data, Geomorphology, 102 (3-4), 435-444. Baldi, P., G. Casula, N. Cenni, F. Loddo and A. Pesci (2009). GPS-based monitoring of land subsidence in the Po Plain (Northern Italy), Earth and Planetary Science Letters, 288 (1-2), 204-212, doi: 10.1016/j.epsl.2009.09.023. Battaglia, M., M. H. Murray, E. Serpelloni and R. Bürgmann (2004). The Adriatic region: An independent microplate within the Africa-Eurasia collision zone, Geophysical Research Letters, 31, L09605, doi:10.1029/2004GL019723. Bendick, R., S. McClusky, R. Bilham, L. Asfaw and S. Klemperer (2006). Distributed Nubia-Somalia relative motion and dike intrusion in the Main Ethiopian Rift, Geophysical Journal International, 165 (1), 303-310. Buble, G., R. A. Bennett and S. Hreinsdóttir (2010). Tide gauge and GPS measurements of crustal motion and sea level rise along the eastern margin of Adria, Journal of Geophysical Research, 115, B02404, doi:10.1029/2008JB006155, 2010. Bruyninx, C. (2004). The EUREF Permanent Network: a multi- disciplinary network serving surveyors as well as scientists, GeoInformatics, 7, 32-35. RELIABLE DEFORMATION RATE FROM GPS NPS DATA 66 67 Casula, G., M. Dubbini and A. Galeandro (2007). Modeling environmental bias and computing velocity field from data of Terra Nova Bay network in Antartica by means of a quasi-observation processing approach, U.S. Geological Survey and the National Academies, Short research paper, USGS OF-2007-1041, doi:10.3133/of2007-1047.srp054. Cenni, N., P. Baldi, E. Mantovani, M. Ferrini, M. Viti, V. D'Intinosante, D. Babbucci and D. Albarello (2008). Short- term (geodetic) and long-term (geological) deformation pattern in the Northern Apennines, Bollettino della Società Geologica Italiana, 127 (1), 93-104. D'Agostino, N., R. Giuliani, M. Mattone and L. Bolci (2001). Active crustal extension in the Central Apennines (Italy) inferred form GPS measurements in the interval 1994- 1999, Geophysical Research Letters, 28 (19), 2121-2124. D'Agostino, N., A. Avallone, D. Cheloni, E. D'Anastasio, S. Mantenuto and G. Selvaggi (2008). Active tectonics of the Adriatic region from GPS and earthquake slip vectors, Journal of Geophysical Research, 113, B12, doi:10.1029/2008JB005860. Dixon T.H., M. Miller, F. Farina, H. Wiang, and D. Johnson (2000). Present-day motion of the Sierra Nevada block and some tectonic implications for the Basin and Range province, North American Cordillera, Tectonics, 19(1), 1-24. Dong, D., T.A. Herring and R.W. King (1998). Estimating regional deformation from a combination of space and terrestrial geodetic data, Journal of Geodesy, 72 (4), 200-214. Dong, D., P. Fang, Y. Bock, M.K. Cheng and S. Miyazaki (2002). Anatomy of apparent seasonal variation from GPS– derived site position, Journal of Geophysical Research, 107 (B4), 2075, doi:10.1029/2001JB000573. Galadini, F. and P. Galli (2000). Active Tectonics in the Central Apennines (Italy): Input Data for Seismic Hazard Assessment, Natural Hazards 22, 225–270. Galli, P., F. Galadini and D. Pantosti (2008). Twenty years of paleoseismology in Italy, Earth-Sciences Reviews, 88 (1- 2), 89-117. Galvani, A. (2009). Indagini geodetiche e geomorfologiche in Appennino Centrale per la caratterizzazione della tettonica attiva, Tesi di Dottorato di Ricerca in Geodinamica, XX Ciclo, Università degli Studi Roma Tre. Gasparini C., G. Iannaccone and R. Scarpa (1985). Fault-plane solutions and seismicity of the Italian Peninsula, Tectonophysics, 117, 59-78. Gomez, F., G. Karam, M. Khawlie, S. McClusky, P. Vernant, R. Reilinger, R. Jaafar, C. Tabet, K. Khair and M. Barazangi (2007). Global Positioning System measurements of strain accumulation and slip transfer through the restraining bend along the Dead Sea fault system in Lebanon, Geophysical Journal International, 168 (3), 1021-1028. Herring, T.A. (2003). MATLAB Tools for viewing GPS velocities and time series, GPS Solutions, 7 (3), 194-199. Herring, T.A., R.W. King and S.C. McClusky (2006a). GPS Analysis at MIT, GAMIT Reference Manual, Release 10.3 (Department of Earth, Atmospheric, and Planetary Sciences Massachusetts Institute of Technology, Cambridge MA, USA). Available online at: http://chandler.mit.edu/~simon/gtgk/GAMIT_Ref_1 0.3.pdf, accessed October 12, 2009. Herring, T.A., R.W. King and S.C. McClusky (2006b). Global Kalman filter VLBI and GPS analysis program, GLOBK Reference Manual, Release 10.3. (Department of Earth, Atmospheric, and Planetary Sciences Massachusetts Institute of Technology, Cambridge MA, USA). Available online at: http://chandler.mit.edu/~simon/gtgk/GLOBK_Ref_10.3. pdf, accessed October 12, 2009. Kenyeres, A. and C. Bruyninx (2004). Monitoring of the EPN Coordinate Time Series for Improved Reference Frame Maintenance, GPS solutions, 8 (4), 200-209. King, M. A. and S.D.P. Williams (2009). Apparent stability of GPS monumentation from short-baseline time series, Journal of Geophysical Research, 114, B10403, doi:10.1029/2009JB006319. King, M.A. and C.S. Watson (2010). Long GPS coordinate time series: multipath and geometry effects, Journal of Geophysical Research, 115, B04403, 1–23, doi:10.1029/2009JB006543. Langbein, J. and H. Johnson (1997). Correlated errors in geodetic time series: implications for time-dependent deformation, Journal of Geophysical Research, 102, 591-603. Langbein, J (2008). Noise in GPS displacement measurements from Southern California and Southern Nevada, Journal of Geophysical Research, 113, B05405, doi:10.1029/2007JB005247. Leonard, L.J., R.D. Hyndman, S. Mazzotti, L. Nykolaishen, M. Schmidt and S. Hippchen (2007). Current deformation in the northern Canadian Cordillera inferred from GPS measurements, Journal of Geophysical Research, 112 (B11401), doi: 10.1029/2007/2007JB005061. Little, M.A., P.E. McSharry, S.J. Roberts, D.A. Costello and I.M. Moroz (2007). Exploiting nonlinear recurrence and fractal scaling properties for voice disorder detection, BioMedical Engineering OnLine, 6 (23), doi:10.1186/1475- 925X-6-23. Mantovani, E., D. Babbucci, M. Viti, D. Albarello, E. Mugnaioli, N. Cenni and G. Casula (2006). Post-late miocene kinematics of the adria microplate: inferences from geological, geophysical and geodetic data, Proceedings of the NATO Advanced Research Workshop on The Adria Microplate: GPS Geodesy, Active Tectonics and Hazards (Veszprem, Hungary, 4-7 April 2004), Berlin, 51-69. Mao, A., C.G.A. Harrison and T.H. Dixon (1999). Noise in GPS coordinate time series, Journal of Geophysical Research, 104 (B2), 2797-2816. McClusky, S., S. Balassanian, A. Barka, C. Demir, S. Ergintav, I. Georgiev, O. Gurkan, M. Hamburger, K. Hurst, H. Kahle, K. Kastens, G. Kekelidze, R. King, V. Kotzev, O. RELIABLE DEFORMATION RATE FROM GPS NPS DATA Lenk, S. Mahmoud, A. Mishin, M. Nadariya, A. Ouzounis, D. Paradissis, Y. Peter, M. Prilepin, R. Reilinger, I. Sanli, H. Seeger, A. Tealeb, M.N. Toksöz and G. Veis (2000). Global Positioning system constraints on plate kinematics and dynamics in the eastern Mediterranean and Caucasus, Journal of Geophysical Research, 105, 5695-5719. McClusky, S., R. Reilinger, S. Mahmoud, D. Ben Sari and A. Tealeb (2003). GPS constraints on Africa (Nubia) and Arabia plate motions, Geophysical Journal International, 155, 126-138. Pesci, A., P. Baldi, A. Bedin, G. Casula, N. Cenni, M. Fabris, F. Loddo, P. Mora and M. Bacchetti (2005). Digital elevation models for landslide evolution monitoring: application on two areas located in the Reno river valley (Italy), Annals of Geophysics, 47 (4), 1339-1353. Pesci, A. and G. Teza (2007). Strain rate analysis over the central Apennines from GPS velocities: the development of a new free software, Bollettino di Geodesia e Scienze Affini, 56 (2), 69-88. Pesci, A., G. Teza and G. Casula (2009). Improving strain rate estimation from velocity data of non-permanent GPS stations: the Central Apennine study case (Italy), GPS Solutions, 13 (4), 249-261. Ray, J., Z. Altamimi, X. Collilieux and T. van Dam (2008). Anomalous harmonics in the spectra of GPS position estimates, GPS Solutions, 12 (1), 55-64. Shen, Z-K., D.D. Jackson and B.X. Ge (1996). Crustal deformation across and beyond the Los Angeles basin from geodetic measurements, Journal of Geophysical Research, 101 (B12), 27957-27980. SOPAC (2009). SOPAC facility ftp site, available online at: http://sopac.ucsd.edu/other/facilities.html, accessed October 14, 2009. Teza, G., A. Pesci and A. Galgaro (2008). Grid_strain and grid_strain3: software packages for strain field computation in 2D and 3D environment, Computers & Geosciences, 34 (9), 1142-1153. Teza, G., A. Pesci and G. Casula (in press). SURMODERR: A MATLAB toolbox for estimation of velocity uncertainties of a non-permanent GPS station, Computers & Geosciences, accepted. Wdowinski, S., Y. Bock, J. Zhang, P. Fang and J. Genrich (1997). Southern California permanent GPS geodetic array: spatial filtering of daily positions for estimating coseismic and postseismic displacements induced by the 1992 Landers earthquake, Journal of Geophysical Research, 102(B8), 18057-18070, doi: 10.1029/97JB01378. Weber, J., P. -Prešeren, T. Dixon, Y. Jiang and B. Stopar (2010). GPS-derived motion of the Adriatic microplate from Istria Peninsula and Po Plain sites, and geodynamic implications, Tectonophysics, 483(3-4), 214-222. Williams, S.D.P., Y. Bock, P. Fang, P. Jamason, R.M. Nikolaidis, L. Prawirodirdjo, M. Miller, and D.J. Johnson (2004). Error analysis of continuous GPS position time series, Journal of Geophysical Research, 109, B03412, doi:10.1029/2003JB002741. Williams, S.D.P. (2008). CATS: GPS coordinates time series analysis software, GPS Solutions, 12 (2), 147-153. Zhang, J., Y. Boch, H. Johnson, P. Fang, J. Genrich, S. Williams, S. Wdowinski and J. Beh (1997). Southern California permanent GPS geodetic array: error analysis of daily position estimates and site velocities, Journal of Geophysical Research, 102(B8), 18035-18055. doi: 10.1029/97JB01380. *Corresponding author: Dr. Arianna Pesci, Istituto Nazionale di Geofisica e Vulcanologia, sezione di Bologna, Italy; e-mail: pesci@bo.ingv.it © 2010 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved. RELIABLE DEFORMATION RATE FROM GPS NPS DATA 68 Pavlov ic cs s