GPS tomography tests for DInSAR applications on Mt. Etna

Tropospheric artifacts of SAR images in a volcanic area like Mt. Etna cause ambiguity in the interpretation of deformations with such technique. It would be useful to measure the delay caused by tropospheric anomalies in synthetic aperture radar (SAR) satellite signals (phase of the back-scattered radar wave) that could be interpreted as deformation. From the delay estimated through the GPS data processing, the tropospheric tomography of electromagnetic waves refractivity, has been performed using the SIMULps12 software. The aim of this study was to perform software synthetic tests by using SIMULps12 applied to atmospheric tomography and to verify the influence of the different GPS geodetic network configurations on obtaining a reliable tomography. Three different anomalies of increasing complexity have been investigated in order to understand the representative parameters of a correct tomography, the best spatial resolution and the portions of space in which the tomography is reliable. The tests also focused on fixing/establishing the a-priori atmospheric model and the critical values of the main parameters involved in the tomographic inversion. To this end, we made a random choice of two days, necessary to define the tomographic problem. Three different network configurations with 15, 30 and 90 GPS receivers were studied. The results indicate that the well-resolved area of tomographic images increases with the number of GPS receivers not linearly, and that the actual GPS network of 42 receivers is capable of revealing/detecting the atmospheric anomalies.


Introduction
Differential SAR (Synthetic Aperture Radar) interferometry (DInSAR) is a satellite-based technique used frequently today to measure the shape of the Earth and its variation in time.
The measurement of ground deformations with DInSAR technique is based on the comparison between two SAR images acquired at different times (called interferograms).If ground deformation occurs between two satellite passes in the same area, the path differences between the sensor and the ground (through the two passes) lead to a difference in the phase of the back-scattered radar wave [Ulaby et al. 1986, Massonnet et al. 1993].
Among various perturbations affecting interferograms, atmospheric artifacts are one of the most significant and, probably, the most difficult to identify and reduce [Zebker et al. 1997, Delacourt et al. 1998, Williams et al. 1998].Indeed, electromagnetic (EM) waves propagate in the neutral atmosphere at velocities lower than the vacuum depending on the local characteristics of the medium.Variations of pressure, temperature and water vapor content cause changes in the effective refractive index along the ray path and therefore may corrupt any phase change measurement due to a deformation and degrade any geophysical parameter estimated from such measurements.The neutral atmosphere (stratosphere and troposphere) is the portion of the atmosphere where the variations that significantly affect the interferograms occur.In areas with prominent topography, such as the volcanoes, atmospheric heterogeneities are relatively frequent and yield the use of DInSAR technique for monitoring purposes in such areas critical.
Mt. Etna was one of the first volcanoes studied by using DInSAR techniques [Massonnet et al. 1995].Due to its prominent topography, several studies have been carried out on Mt.Etna since 1996 to better understand the influence of the neutral atmosphere on the measurement of ground deformations viewed from DInSAR.These studies investigate the possibility of evaluating and correcting the tropospheric effects of SAR images through different approaches: by using the integrated tropospheric delay obtained from the GPS analysis [Beauducel et al. 2000], starting from empirical data used to predict atmospheric conditions through numerical models [Delacourt et al. 1998, Wadge et al. 2002, Webley and Wadge 2004].All these studies attempt to estimate the atmospheric temporal variability and then compare the results with SAR interferograms.Given that the SAR satellites are side looking radar (SLR), the effects of atmospheric anomalies on DInSAR measure-

Subject classification:
Volcano monitoring, Satellite geodesy, GPS, Tomography, Etna, SAR, DInSAR, Refractivity, Delay.ments can be evaluated only through the knowledge of the 3D structure of the atmosphere.
Starting from the previous studies [Bonforte et al. 1999, Bonforte et al. 2001, Bruno et al. 2007], we propose some advances aimed at assessing a method to obtain a suitable tomography of EM waves propagation velocity of the lower atmosphere, from delays measured by the Mt.Etna GPS network.
In this paper, we present a methodological investigation of the minimum GPS network configuration required to obtain a 3D refractivity tomography over Mount Etna.

Method
Although GPS is commonly used to estimate the delays of EM waves through the atmosphere and some studies have recently exploited this characteristic to obtain tomographies [Hirahara 2000, Hoeven et al. 2002, Nilsson and Gradinarsky 2006, Nilsson et al. 2007], few studies have been conducted on Mount Etna.
Our research starts from previous achievements both in seismic tomography [Aloisi et al. 2002] and GPS atmospheric sounding [Bonforte et al. 1999[Bonforte et al. , 2001] ] to evaluate a method to produce a tomography of the lower atmosphere (neutral atmosphere) over the Mt.Etna area, by appropriately modifying a software originally implemented for seismic tomography, SIMULps12 [Bruno et al. 2007].
As in seismic tomography, where it is necessary to fix the seismometers and earthquake positions, as well as the a-priori wave velocity fields, here we have to set corresponding parameters to compute the atmospheric EM wave velocity tomography.We assume the satellites as source of signal and the GPS receivers as seismic stations [Bruno et al. 2007].In particular, GPS satellites are considered as points of shot in seismology, namely earthquakes with well-known position and time.To study the tomography problem, in this work we use an interpolative function defined by values specified at nodal points (knots) within a three-dimensional grid [Aki 1993, Lee andPereyra 1993]; a complete discussion of the problem is examined in Thurber and Ellsworth [1980], Thurber [1983], Thurber and Aki [1987], Thurber [1993] and Eberhart-Phillips [1993].Some processing parameters were analyzed and several spatial configurations were tested to determine the best knots position.A few tests were also carried out to assess whether benchmarks are suitable to obtain a "reliable tomography".

GPS satellites (position and selection)
As in local earthquake tomography (LET), in the case of atmospheric tomography we suppose to ideally enclose the satellites, the receivers and the portion of atmosphere we want to analyze, inside a parallelepiped.The base of the parallelepiped of the investigated area should be comparable with the Etna area.Although the typical GNSS applications use a cutoff angle of 10° or less, for our purpose we assume a cut-off angle of 15°, because in the GPS analysis the atmospheric anomalies are better modeled and also because the straight line propagation approximation can be applied.Since the GPS satellites orbit at altitudes between 18,000 and 22,000 km, assuming a cut-off angle of 15° for GPS measurements, the base of the parallelepiped results extremely large.For instance, if we consider a satellite altitude of 18,000 km, we obtain a base for the parallelepiped of 1.3*10 5 km 2 , which is oversized if compared to the Mt.Etna area, namely about 7.5×10 2 km 2 .
Because we are only interested in the neutral atmosphere, we reduced the dimension of our system by "shifting down" the satellites ideally along the directions connecting them with the GPS receivers.This shifting does not affect the goodness of the result since the first order ionospheric contribution on GPS delay, as explained in Aranzulla et al. [2013], is eliminated by an "iono-free" linear combination [Klobuchar 1996].The new satellite coordinates (here called synthetic coordinates) result from the intersection of the straight line (realistic approximation for a cutoff of 15°) between satellite and receiver with the ellipsoid 10 km away from the Earth ellipsoid (Figure 1).
Indicating with A(x A ,y A ) the actual satellite coordinates and with B(x B ,y B ) and C(x C ,y C ) the coordinates of two receivers in Figure 1, the points A´(x A´, y A´) or A˝(x A˝, y A˝) define the coordinates of the new satellite position for receiver B and C at the reference surface at 10 km of elevation.According to this model we have "reduced" the height of satellites to 10 km above the Earth's ellipsoid, inside a parallelepiped with a base of about 70×70 km 2 , which is much closer to the dimensions of the Mt.Etna area.The geometrical distance between the satellite and the GPS receiver will not therefore be the actual one, but depends on the new coordinates of the satellite.That mean the tropospheric delay is due only to the first 10 km of atmosphere (troposphere) and not to the entire neutral atmosphere (troposphere and stratosphere), which extends as far as about 50 km above sea level.This seems a reasonable assumption, because the troposphere contains about 90% of the air mass and 99% of water vapour of atmosphere, which is the main source of the atmospheric delays for the GPS EM waves [Essen and Froome 1951].

The atmospheric model
In SIMULps12, the wave velocity is computed on knots of a three-dimensional grid, in which velocity varies continuously in all directions, with linear interpolation among knots [Evans et al. 1994].Knots are defined by the intersections of a set of planes in three orthogonal directions, one horizontal and two vertical.The a-priori model we used is characterized by a grid in which velocity of EM waves in the atmosphere varies only with respect to the altitude.EM waves velocity depends on refraction index n (excluding the imaginary part of the refractive index that indicates the amount of absorption loss when the EM wave propagates through the material) according to the expression v=c/n, where c is the velocity of light in vacuum [Saastamoinen 1972[Saastamoinen , 1973]].Since the refraction index is related to refractivity N by: we can express EM waves velocity through the following expression: According to the Essen and Froome formula adopted by the International Association of Geodesy in 1963 for atmospheric refractivity applied to GPS satellites [Essen and Froome 1951], refractivity has the following expression: (3) where P is the partial pressure of dry air (expressed in mbar), T is the atmospheric temperature (expressed in Kelvin) and e is the partial pressure of water vapour (expressed in mbar), while the empirical constants k 1 , k 2 , k 3 , have the following numerical values [Rüeger 2002]: k 1 = 77.6890K/mba; k 2 =71.2952K/mba; k 2 =375463 K/mba.Equation ( 3) is valid assuming that air is an ideal gas and for our purposes this approximation is acceptable.By applying Equation (3) we may compute the refractivity at an arbitrary altitude.Unfortunately, we do not have regular atmospheric sounding on Mt.Etna and the nearest launch of weather balloons is performed from Trapani (about 200 km west of Mt.Etna).Thus we have to model P, T and e, to apply Equation (3) and to determine the apriori velocity model that we use for the tomography.In the present study, assuming a standard atmosphere for different altitudes, we used the "Normal Temperature and Pressure" (NTP) reference conditions (20 [°C], 293,15 [K] and 1013.25 [mbar]) at soil.
To evaluate the impact of the presence of water vapour on the EM wave velocity value, we calculate refractivity and consequently propagation velocity of radio waves at different heights, for the two extreme conditions: dry and saturated.Since the highest value of water vapour pressure at each height corresponds to the saturated water vapour pressure, by the wellknown Clausius-Clapeyron equation, we are able to compute the refractivity into two extreme conditions: the dry refractivity using the first term of Equation ( 3), the saturated refractivity using all the three terms of Equation ( 3).
At sea level, according to the atmospheric model build with the Essen and Froome Equation (3), the water vapour maximum impact on the EM wave velocity is about 3*10 4 m/s.The SIMULps input/output velocity value of the knots is expressed in km/s with 4 significant digits.Consequently, at this level of study, SIMULps reveals large actual anomalies due to tropospheric water vapor.Although this strong limitation will be resolved in the future, the main issue of determining the optimal network configuration for tomographic purposes, will be solved.Furthermore, despite this limitation, the assessment of the optimal network configuration keeps its general validity unchanged.

Knots fixing
As previously introduced, the tomographic system is made up of the satellites and receivers position together with the tomographic volume enclosing the atmosphere that we want to study.Ideally speaking, if an atmospheric anomaly is inside the tomographic volume during the investigated time, the GPS observations should be able to reveal it in the right place; in any place within the tomographic volume.This is not true in the actual case, because the properties of the tomographic volume (i.e., the knots velocity value) depend on the discretization of the investigated volume and (2) also on the amount and geometry of the recorded observables.Knots where velocity is computed are defined by the intersections of a set of planes in three orthogonal directions, one horizontal and two vertical.Furthermore, the position of the knots depends on the user's assumptions.Typically the planes are closely spaced in the centre of the target volume, where most rays should be, and more widely spaced in the periphery.
Every knot can either be held fixed or included in the inversion, except the outermost knots which are always fixed.The tomographic volume is set as a parallelepiped with the base of 60×60 km centered on the Mt.Etna central crater and the height of 10 km above ellipsoid level.SIMULps software performs a damped-leastsquares inversion, according to which the norm of the model perturbations (the complexity of the model) is weighted and combined with the squared data misfit through a parameter called "damping factor".
To fix the knots we use four parameters from the output SIMULps software: the hit count (KHIT), the derivative weight sum (DWS), the resolution (RS) and the most important spread function (SF).The KHIT represents the number of ray paths running near the considered knot.The DWS provides an average relative measure of the density of rays near a given velocity knot, describes the amount of data actually constraining the velocity at that knot.The magnitude of DWS depends on the selected step length of the arc connecting the start and end points, smaller step lengths yield larger DWS, therefore DWS provides only a relative measure of the ray distribution.To study how well each knot was illuminated, initially the KHIT and DWS values were investigated for each knot.The knots have been moved so that: a sufficient number of rays hit each knot, but also, to allow that a large number of knots were inserted in the tomographic computation.We fix a KHIT threshold value of 9, the knots where KHIT value is less than 9 are excluded from the inversion.Furthermore all knots with DWS < KHIT are excluded.Because DWS is affected by the distance of the rays from the knot, a large DWS indicates that the velocity at the knot is based on a large body of data [Toomey and Foulger 1989].The diagonal elements of resolution matrix represent the resolution (RS) of each model parameter.The RS provides the weights utilized in the averaging process of model parameters to estimate.Values close to 1 indicate well-resolved parameters, whereas lower values show us that the computed velocities are smeared over a larger model volume [Gökalp 2007].The analysis of the resolution matrix is given by the spread function [Menke 1989].The definition of a SF (Equation 4) provides a synoptic view of RS [Toomey and Foulger 1989]: (4) where s kj is an element of the RS matrix, N is the number of the parameters, D jk (Equation 5) is the distance from j-th to k-th knot. (5) The SF represents the spread of RS relative to the diagonal element for each parameter, and the averaging vectors, indicating the volumes over the velocities which are averaged.Several configurations were tested to determine the best knots position, analyzing for each of them the values of KHIT, DWS and SF.In a series of test inversions, the horizontal and vertical separation of knots was varied from 2 to 7 km.The results of simultaneous inversion for all the nodal configurations are better for the larger grid spacing, hence we chose a configuration grid having 245 knots.Concerning the horizontal layers, we fix them at 2, 3, 5, 7, 9 km above the ellipsoid level (layer A, B, C, D, E, respectively) (Figure 2b).The better configured layers for the inversion are the upper layers, specifically the layers at 5, 7 and 9 km from the ellipsoid level.The gap of the horizontal layers and the knots produced by the vertical planes intersection are reported in Figure 2a   where V 1 , V 2 , V 3 , V 4 , V 5 , V 6 , V 7 , V 8 are the EM wave velocity value of the knots surrounding the segment and X, Y, Z are: (7) The time spent to cover the distance between the new satellite position and receiver of the k-th ray is: where n is the number of segments Ds i in which the ray between the new satellite position and receiver, has been divided.

GPS network tests
We carried out a few tests to evaluate which GPS configuration (i.e. the minimum number and position of GPS receivers, satellites and knots) is suitable for the tomography and also to verify if the software can identify and correctly place some synthetic atmospheric anomalies.
Although today's GPS Etnean network is made up of about 42 permanent stations covering the volcanic edifice fairly uniformly, we tested three different GPS network configurations, with 15 and 30 and 90 GPS re-ceivers.The configuration with 15 stations corresponds to the GPS permanent network until 2005 (Figure 4a), while the one with 30 stations also includes stations until 2009 [Palano et al. 2010] (Figure 4b).The configuration with 90 stations is not shown in the figure because it is an ideal network that includes both permanent stations and all the benchmarks currently on Mt.Etna.Since an increase in the number of GPS stations improves the results, we focused on the minimum number of stations needed to reveal anomalies in troposphere.
In the test with 30 GPS receivers, two spatial configurations of GPS satellite were considered to investigate the variability of the well-resolved knots (and ; ; . GPS TOMOGRAPHY TESTS FOR DINSAR APPLICATIONS ON MT.ETNA  The position of satellites is obtained by reading the IGS precise ephemeris (http://ww.iers.org),which are provided with a sample rate of 15 minutes.A significant difference in shape of the areas is visible in Figure 5; this demonstrates that the shape and size of the investigated area depend on the positions of the satellites.The tests were carried out assuming a set of "structures" (i.e.volume in which the velocity of waves varies from the surrounding space) of well-defined shapes.The three structures considered are: points, cube and checkerboard.In particular, the anomaly called "point" represents an alteration of one knot on layer C, placed at SE of the central crater; the anomaly called "bubble" represents an alteration of eight knots (four for each layer) on layer C and D, placed from the central crater to SE direction; the anomaly called "checkerboard" represents the alteration of all the knots alternately.Although the last anomaly has no physical meaning, it is very important to test the entire tomographic system.Figures 6, 7 and 8 show the three anomalies (in the columns called "model"), specifying that the variations applied to the modified knots are equal to the 0.5% of the unperturbed model.We calculate the travel time (called "synthetic travel times") with the presence of one of the atmospheric anomalous structure above and then, starting from the network and satellites positions and unperturbed atmospheric model of Section 2, we compute the tomography.Because of the small perturbations introduced the initial model, the synthetic travel time for the three anomalies is approximately equal and varies from 2.5 to 9.8 sec.The reader remember that: the input velocity value of the knots in SIMULps is expressed in km/s with 4 significant digits.This means that, compared with the actual travel time, the travel time computed by SIMULps has to be multiplied of a 10 -5 factor.Furthermore the synthetic perturbations influence the wave velocity of the SIMULps model of about 15 m/s, about two orders of magnitude larger than the actual case.Although this is a drawback for our purposes, however, because we are testing the best network configuration, it does not represent an actual limitation on the assessment of the optimal network configuration.Unfortunately the number of significant digits used by SIMULps, does not allow to go further, so to overcome this drawback it is necessary to face with this problem with different SW tools, which it is what is currently ongoing, as further depicted.
The first analyzed configuration includes 15 stations and 8 satellites, it shows the pattern of the values of the velocity, for the three investigated structure of anomalies and the tomographic inversions, respectively.The results of the tests are shown as maps of gridded values for different layers, at different elevations, together with the contour plot of the SF using the satellite spatial configuration of June 22, 2005.It is worth noting that the GPS satellites geometry/configuration causes differences between the spread function shapes.We visualize data using a "natural neighbour" gridding method.As visible in Figure 6, the tomography provides information about the presence of atmospheric anomalies in the studied area.In the inversion of the "point" and "cube" structures, for instance, the tomography produces a sort of "reverberation" of the anomalies in the adjoining layers but the maximum magnitude of variations are revealed in the right places.The test on the "checkerboard" structure clearly shows that inside places where SF is less than 2 the anomalies were fairly well-resolved especially in the upper layers, while outside the anomalies become blurred and unrecognizable.Furthermore, this last test shows that the image quality increases with the altitude of the considered layer.
The second analyzed configuration includes 30 sta-tions.In accordance with the higher number of observables with respect to the case of 15 stations, Figure 7 suggests a better reconstruction of the anomalies.Indeed, although the reverberation of the anomalies in the adjoining layers is still present with the same magnitude of the previous case, the reconstructed values at different knots are generally closer to the model values than the previous corresponding tests.Furthermore, the contours of the SF of Figure 7, are wider than those reported in Figure 6.As concerns the "checkerboard inversion" in Figure 7, this test confirms that a value of 2 for SF is the most suitable for defining the volumes where the anomalous structures are better reconstructed.For the 30-stations case, we observe a substantial improvement in the reconstructed anomalies.In fact, comparing Figure 7 with Figure 6, we can see a better overlapping of the reconstructed knots value with respect to the model knots value.
As a last test, we checked the maximum input data that SIMULps12 may process.We processed information provided by a configuration that includes 90 receivers.Although this test may be considered an "unrealizable test case" it is useful to understand the limits of the adopted software.The reconstructed values from this simulation are shown in Figure 8.
The reverberation of the anomalies in the adjoining layers is still present with the same magnitude of the previous cases.The value of reconstructed knots is GPS TOMOGRAPHY TESTS FOR DINSAR APPLICATIONS ON MT.ETNA still different from the actual value and equal to the values found in the previous case.Even though the number of receivers is much higher than the previous simulations, the growth of the area covered by the con-tour plot of the SF is not comparable with the increasing of number of stations.
To assess and quantify the tomographic results, for each layer and network configuration, we compute a  sort of variance through the following equation: (9) where V tomo (x i , y i , z i ) represents the velocity of the i-th knot computed through the tomography, V model (x i , y i , z i ) is the expected velocity value of the i-th knot, m is the number of well-resolved knots in the considered layer X.Table 1 shows the results of all the tested network configuration and anomalies.For each layer we show the variance value, expressed in parts per million, related to each of the tested anomalies.The third column also shows how many knots have been resolved with the SF threshold set.As expected, the value of the variance decreases when the number of GPS receivers increases.
This test suggests that, also using a great number of receivers, the increase in RS is lower compared with the cost of the survey and the computational efforts.The tests above also led to conclude that 2 can be taken as the best value for the SF for the quality of the obtained results.This method will be exploited in future studies to obtain an estimation of the effects of anomalies in the distribution of water vapour in atmosphere on the SAR images.Indeed, knowing the spatial location of these anomalies, from simple geometrical considerations, it will be possible to estimate the location of the false deformation on SAR images and correct them.

Discussion and conclusion
The aim of this study was to evaluate the best GPS network configuration and the minimum number of GPS receivers needed to obtain the tropospheric "corrections" with respect to the starting velocity model and hence to carry out the atmospheric tomography.Although GPS is commonly used to estimate the delays of EM waves through the atmosphere, relatively few studies exploit this characteristic of GPS to obtain tomographies of the atmosphere in general, and on Mt.Etna in particular.
We have carried out many synthetic tests to modify the seismic tomography applications by considering the geometrical aspect of the problem and by optimizing the steps of data processing.The main modifications of the geometry of this problem with respect to the seismic one, concern the GPS satellite shifting down, which was necessary to reduce the base of the investigated area, and the definition of a 3D grid.We have tested three different network configurations with 15, 30 and 90 GPS receivers.Clearly, the results of the tests show that the well-resolved areas of tomographic images increase with the increasing number of GPS receivers but, in the last configuration with 90 GPS receivers, the wide well-resolved area is not comparable with the increase in number of the receivers.Results in The columns show: the name of layer, the network station configuration, the number of well resolved knots, the variance (in ppm) for the three anomalies respectively.
ers show that the gap between the 15 and 30 stations is wider than the gap between the 30 and 90 stations.Although the number of stations has a larger gap between 30 and 90 stations (Table 1, column 2), the number of well resolved knots does not have the same trend (Table 1, column 3).We can summarize this result by stating that: the well-resolved areas of the tomographic images do not increase with the number of GPS receivers linearly.Thus these tests suggest that current GPS network of 42 receivers, is able to reveal the atmospheric anomalies.The tests were also devoted to fix the a-priori atmospheric model and the critical values of main parameters involved in the tomographic inversion.The results obtained allowed us to understand what kind of anomalies can be revealed by the GPS network placed on Mt.Etna.Although the simulations show that the SIMULps software cannot give a measure (values) of the anomalies caused by water vapor, the results show that the SIMULps tomographic engine is able to reconstruct the shape and the location of complex anomalies.Based on these preliminary synthetic tests, we will apply this approach to the analysis of the atmosphere with geometric characteristics of the SAR satellites passing over Mt.Etna.To this end, we will estimate the path delays of electromagnetic waves by using the GPS data collected by the Mt.Etna network, processed by GAMIT software [Herring et al. 2010].The results of the tests we have performed suggest that: the modified version of SIMULps 12 is able to reconstruct the geometry of most relevant anomalies through the GPS data of Etna network, but it is an unuseful tool for estimating actual ones.Future studies, in fact, allow to fix the highlighted problems and overcame the actual drawbacks, by building a specific software that uses the tomographic approach of SIMULps but with better precision.The procedure we have presented in this work will be applied to actual cases by comparing the results with direct radiometric or satellites spectroscopy measurements.By using the obtained tomographic results, we will be able to compute the tropospheric delay, correcting the tropospheric effects on SAR images.

Figure 1 .
Figure 1.Scheme of satellite projection from the actual position A(x A ,y A ) to the positions A´(x A´, y A´) and A˝(x A˝, y A˝) with respect to stations B(x B ,y B ) and C(x C ,y C ).

Figure 2 .
Figure 2. (a) Plan view of the knots fixed for tomographic inversion.(b) A section of horizontal fixed layers.The values of DWS vary substantially from one layer to another with the biggest values in the upper layers.

Figure 3
Figure 3 shows an example of the value of DWS versus SF; the trend of this figure has a clear change around the value of 2.5.Because well-sampled knots have a large value of DWS versus smaller values of SF, we chose a spread function value of 2 as an upper limit of well-resolved knots [Toomey and Foulger 1989].

Figure 3 .
Figure 3. DWS versus SF; a clear change of trend around the abscissa value of 2.5 allows choosing the upper limit of well-resolved knots.

Figure 4 .
Figure 4. Plan view of Mt.Etna GPS network configuration (a) until 2005 with 15 benchmarks and (b) until 2009 with 30 benchmarks.

Figure 5 .
Figure 5. Plan view contour plots of the SF values obtained by using 30 GPS receivers.Dashed lines enclose areas with SF values less than 3, while solid lines denote area with SF values less than 2. (a) From the satellites position of June 22, 2005, at 09:00.(b) From the satellites position of July 27, 2005, at 09:00.

Figure 6 .
Figure 6.Tomographic images of June 22, 2005, at 09:00 with the 15 stations.Each reconstructed image shows also the spread function contour plot.Dashed lines enclose area as before with SF values less than 3, while continuous lines enclose area with SF values less than 2.

Figure 7 .
Figure 7. Tomographic images of June 22, 2005, at 09:00 with 30 stations.Each reconstructed image shows also the spread function contour plot.Dashed lines enclose area with SF values less than 3, while continuous lines enclose area with SF values less than 2.

Figure 8 .
Figure 8. Tomographic images of June 22, 2005, at 09:00 with 90 stations.Each reconstructed image shows also the spread function contour plot.Dashed lines enclose area with SF values less than 3, while continuous lines enclose area with SF values less than 2.

Table 1 .
Variance of all the three tested network configuration (15, 30 and 90 GPS stations) and anomalies (point, bubble and checkerboard).