Preconditions to ground based GPS water vapour tomography

. The GPS water vapour tomography is a new technique which provides spatially resolved water vapour distributions in the atmosphere under all weather conditions. This work investigates the information contained in a given set of GPS signals as a precondition to an optimal tomographic reconstruction. The spatial distribution of the geometric intersection points between different ray paths is used to estimate the information density. Different distributions of intersection points obtained from hypothetical GPS networks with varying densities of GPS stations are compared with respect to the horizontal and vertical resolution of a subsequent tomographic reconstruction. As a result some minimum requirements for continuously operating extensive GPS networks for meteorological applications are given.


Introduction
Water vapour is one of the most important constituents of the troposphere which effects the weather in various ways.Temporally and spatially resolved humidity information are therefore essential in many fields of meteorological research and applications, such as weather forecast, especially precipitation forecast and nowcasting, hazard mitigation and water management.These applications would require much more detailed vertical profiles than the currently available radiosonde data with a temporal resolution between 6 and 12 h and a horizontal resolution of a few hundred kilometers.Such information can only be provided by remote sensing techniques, e.g.microwave radiometers, lidars or GPS.The rapidly developing field of GPS meteorology makes use Correspondence to: M. Bender (bender@uni-leipzig.de) of existing GPS satellites and ground stations and moisture data can be obtained with small efforts compared to the investments required by other techniques.The preconditions for the application of tomographic techniques have first been developed by Bevis et al. (1992) and have been utilized by a large number of experiments.Some GPS networks provide already moisture data (Dick et al., 2001;Troller et al., 2006) and plans to set-up GPS based water vapour retrieval systems exist in most European countries.
GPS signals carry information about the water vapour integrated along the transmitter-receiver path, i.e. the total water vapour between the GPS satellite and the GPS ground-station (Kursinski et al., 2000).The 3-D water vapour distribution can be reconstructed from a large number of such integrated values by means of tomographic techniques (Flores et al., 2000a;Subbarao et al., 1997;Kunitsyn and Tereshchenko, 2003).This requires the solution of an inverse ill-posed problem with incomplete data (Natterer, 2001).The term "incomplete data" refers to the fact that only a very limited number of line integrals (compared with more usual tomographic applications, such as CT or PET) can be obtained with a given constellation of GPS satellites and receivers.The reconstruction of the atmospheric water vapour from GPS data alone would therefore be limited to cases of very dense local networks with receiver distances ≤1 km.In more sparse regional or national networks a priori knowledge, constraints and meteorological observations become more important.The lack of GPS-data can to a certain extent be compensated by using e.g. a good initial distribution, by introducing inter-voxel constraints and by using additional non-GPS observations (Bevis et al., 1992).Due to their inverse character all reconstruction techniques are very sensitive to small variations in the input data.It is therefore important to verify the quality of the input data, i.e. the GPS data, before the results of the reconstruction can be analysed.Naturally, each single GPS observation must be processed and validated carefully to remove the usual GPS errors from the data, e.g.

Published by
Copernicus Publications on behalf of the European Geosciences Union.clock errors and multipath effects.This work concentrates on the quality of the GPS dataset as a whole which depends on the effective distribution of transmitter-receiver paths in the atmosphere.Different tomographic reconstruction techniques have successfully been used to obtain the water vapour distribution in the troposphere, each of them having its own benefits (Flores et al., 2000b;Bastin et al., 2005;Troller et al., 2006;Foelsche and Kirchengast, 2001).Most of them need a set of internal parameters which have to be optimised for atmospheric reconstructions.Furthermore, all these inverse methods depend strongly on the chosen initial conditions, e.g. the initial water vapour field, and require additional information which can be used during the reconstruction process.The finally reconstructed moisture distribution represents not only the quality of the GPS input data but also the skills of the reconstruction algorithm, the chosen parameters, the initial conditions, etc.
To validate the available GPS data independently from a specific reconstruction technique we developed a procedure which analyses the intersection points between the ray paths to find a first measure of the information provided by a certain GPS data set.Tomographic reconstruction algorithms require not only a large number of views through the region of interest but also views from different angles.A GPS station can usually receive a large number of GPS signals in short time but, unfortunately, most of the ray paths are almost parallel and cannot provide any information about the vertical structure of the atmosphere.The total number of ray paths in a given region is therefore no indication for a high amount of information about this region.Only if neighbored stations are close enough to add views from different angles through the same region the information would be sufficient to reconstruct vertical structures.The density of intersection points, together with the intersecting angles, can be used to quantify the information content.The spatial distribution of the intersection points gives a first approximation of the horizontal and vertical resolution and identifys regions with sparse data.In this work the results of a validation study will be presented and some basic requirements of a nationwide GPS water vapour retrieval system will be worked out.

Geometric constraints of ground-based GPS
The GPS tomography has to recover the humidity structure of the troposphere from a limited number of line integrals measured along a series of ray paths between the GPS satellites and the GPS ground stations.The main requirement is to get as many ray paths as possible which intersect any given region from a wide range of angles.Due to its integrated character the information provided by a single ray path cannot be located along its path.A large number of intersecting rays are required to locate the information.Therefore, it is advantageous to concentrate on the intersection points rather than on the ray paths itself.
Up to 50% of the atmospheric water vapour is located in the boundary layer.It is therefore important to estimate the performance of GPS techniques for altitudes below ∼2 km.One precondition is the existence of intersection points between different ray paths in this range.The minimum height of intersection points can be deduced from geometric considerations neglecting the deviations of Earth's shape from an ideal sphere and ray bending.
Assuming a situation as shown in Fig. 1 with a vertical ray through a station A and a ray with an elevation ε through a station B in a distance S AB leads to a height = 1 + 2 of the intersection point.The first term 1 describes the effect of Earth's curvature and depends on S AB only.The second term 2 is due to the ray elevation ε.Both quantities can be written as a function of the inter-station distance S AB : where R 0 =6371.23 km is Earth's radius.as derived from Fig. 1 gives the height of an intersection point between two rays of the elevation ε and π/2.The intersection point between two opposite coplanar rays both with the elevation ε is located below this hight and can be obtained from Eqs. ( 1) and ( 2 intersection points in Fig. 2 have been plotted using interstation distances S AB /2 and different values of ε. 2 is the major part which dominates even for very small elevations. As can be seen from Fig. 2 some first intersection points can be expected below altitudes of 1-2 km for inter-station distances below 50 km and/or ε min ≤10 • .However, the tomographic reconstruction would require an evenly distributed availability of intersection points in a given vertical range, not only some isolated points.A minimum elevation of ε min =7 • . . . 10 • as used by most GPS receivers would therefore provide insufficient information in the boundary layer below ∼2-3 km (Bevis et al., 1992).
The tomographic reconstruction of a given region requires not only a large number of intersecting ray paths but also views from a wide angular range.The angular range which can be observed by ground-based GPS is limited to 0 • . . .λ 2 , where λ 2 is given by (Fig. 1) The range not visible from Earth's surface depends mainly on ε min and varies between ∼5 • and ∼15 • (Bevis et al., 1992).The GPS tomography is therefore not substantially restricted by geometric constraints.In principle, any region in the troposphere can be investigated from a wide range of angles, even in the boundary layer.The only restriction is the very small number of satellites and the rather small number of ground stations.

Simulation of GPS networks
Various computer simulations have been performed to investigate the correlation between the receiver density and the spatial resolution of intersection points in the troposphere.The quality of the data provided by a certain GPS network depends on the extension of the GPS network with respect to the region of interest and the receiver density within the network.As ray paths which partially leave the reconstruction area must be rejected by the tomography, the number as well as the angular range of ray paths decreases rapidly at the boundaries and leads to artefacts in the reconstructed fields.Such artefacts could be minimized by extending the receiver network and the reconstruction area so far beyond the region of interest that low elevation ray paths pointing to the center of the reconstruction area can be processed.Assuming ε min =7 • additional 75 km would be ideal.Well resolved vertical water vapour profiles require a dense network of receivers, especially if a good representation of the boundary layer is expected.
A large network extended over 2000×2000 km has been simulated to minimise the boundary effects in relation to the complete area.Within this region the GPS stations have been placed on a regular grid with spacings between 100 km and 15 km (Table 1).The topography has not been considered and all stations are located on the surface of the WGS84 ellipsoid.All these networks are sparse considering the requirements given by Bevis et al. (1992) who claim for a receiver density which is comparable to the scale height of water vapour.
Such a hypothetical network could be related to the current situation in Germany where a network of about 280 GPS stations is operated (Reigber et al., 2004).The mean inter-station distance in this network is ∼40 km and the tomographic reconstruction on a horizontal 50 km grid might be realistic.
The ray paths have been simulated using realistic positions of GPS satellites obtained from ephemerides and modeled networks of GPS ground-stations.The positions of 29 GPS satellites at a certain date have been computed in hourly steps to illustrate the diurnal variation.Under these idealised conditions 6-10 satellites are visible above each GPS groundstation.A set of connecting lines between each station and all satellites visible in its local horizon system is computed in the first step.This set of ray paths is subsequently analysed.Due to the integral character of the GPS data the information is mainly located at the intersection points of two rays.All possible intersection points of any pair of ray paths are therefore computed using basic properties of the analytical geometry.As two lines seldom intersect in three dimensions an "intersection point" is assumed where the minimum distance between two lines is less than a given threshold value.
The horizontal and vertical thresholds have been chosen in analogy to a hypothetical reconstruction grid.
To estimate the preconditions to tomographic reconstructions the distribution of the intersection points within different spatial grids has been evaluated.Grids with horizontal spacings of 7.5 km, 15 km and 30 km have been used.The vertical spacing has been chosen to 500 m in all cases and the region up to the tropopause at 12 km has been considered.The total number of grid cells is 1 646 592, 443 688 and 102 912, respectively.
The number of navigation satellites will increase considerably in the next years as the European Galileo system will become operational and new GLONASS satellites will be started.To simulate the advantage of these future satellites a given GPS constellation has been extended by adding a second GPS constellation 4 h ahead.

Number and temporal variation of intersection points
The total number of intersection points in a given region gives a first measure of the achievable resolution.This number depends on the available satellites and ground stations and changes with the diurnal cycle of satellite orbits.The diurnal variation of the GPS data leads to a considerable variation in the total number of intersection points with about a factor of four between the minimum and maximum val- ues (Fig. 3).If only rays with a certain minimum elevation were considered, the diurnal variation became even more pronounced as ε min increased.A factor of 15 . . .20 between the "best" and the "worst" constellation is possible.
An important parameter is the achievable resolution of the reconstructed field.A low resolution guarantees a large number of intersection points in the majority of cells but provides little information about the spatial distribution while the rconstruction on a highly resolved grid might fail due to a large number of empty or sparsely filled cells.It is therefore difficult to find the best resolutions which guarantees reliable results.One parameter could be the number of "empty" cells.This number should obviously not be too high.According to Fig. 5 this can only be achieved if the horizontal resolution of the reconstruction is below the mean distance between ground stations.But even in this case most cells contain very few intersection points (Fig. 6).A larger number of satellites would improve the situation.A similar effect can be achieved by extending the observation period up to several hours.
One important parameter is the total number of available ray paths through the atmosphere.The number of intersection points is connected to the number of ray paths by a quadratic law (Fig. 4).Even small improvements in the GPS constellation either by adding new ground stations or by the start of further satellites produce much more detailed information.The total number of intersection points as obtained for networks with different densities of stations is shown in Fig. 3 (solid lines).The situation improved considerably if the number of satellites was doubled (Fig. 3, dashed lines).To estimate the resolution of the GPS data these numbers have to be related to a 3-D grid which might be used for the tomographic reconstruction.Unfortunately, the distribution Ann.Geophys is very inhomogeneous and average numbers of intersections per grid cell are not too meaningful.Much more significant is the portion of empty grid cells containing no intersections points (Fig. 5).Nearly "empty" regions are difficult to reconstruct and would require extra information to obtain reliable results.About 10% of empty cells might be tolerated as this is comparable to the boundary cells in large grids.Much larger portions indicate that major regions are not represented by the GPS data.As can be seen in Fig. 5 this would require a network with average inter-station distances not exceeding the horizontal dimensions of the grid.
A reliable tomographic reconstruction would require a large number of intersection points per grid cell.As shown in Fig. 6 this can only be achieved if the inter-station distances are well below the grid spacing.Otherwise, most cells contain only very few intersection points.With an increasing number of stations and/or satellites the maximum of the distributions moves to larger numbers of intersection points per cell.A similar effect can be achieved by extending the observation period and to use all data for the same reconstruction (Fig. 6).If the observation period fits into the assimilation interval of a numerical weather model, e.g. 3 h, no information is lost.The distribution broadens considerably if the observation period is extended, indicating a highly inhomogeneous spatial distribution of the GPS data.The inhomogeneity is another problem for tomographic techniques which iteratively approximate the best solution.Neighboring cells with a highly different number of intersection points often disequilibrate the reconstruction.The result represents in such cases only the inhomogeneity of the data but not the real distribution of the reconstructed quantity.

Vertical distribution of intersection points
The vertical distribution of intersection points is one of the most important parameters because it gives a first estimate of the vertical resolution which might be obtained from the data.Not only the total number of intersection points but also their distribution depends strongly on the elevation cutoff angle ε min (Fig. 7).The vertical distribution shows some very promising features if all geometrically possible intersection points are taken into account (ε min =0 • ): The distribution decreases slightly with increasing height, starting from a maximum at very low altitudes (<500 m).Unfortunately, most of these rays propagate almost horizontally through the atmosphere and can usually not be detected.With an increasing elevation ε min the number of available intersection points decreases rapidly and the maximum of the distribution moves up to several kilometers leaving little data in the boundary layer.About 60% of the data belong to ray paths with ε min <5 • , less than 15% are left if only elevations ≥10 • are considered.It can be seen from Fig. 7 that almost no data below 5 km (50 km inter-station distance) or 3 km (25 km) are available if ε min is set to 10 • .At the current time most receivers work with elevations between 7 • and 10 • .Advanced receiver and processing techniques should make elevations down to 2 • available (grayed area in Fig. 7) (Pany et al., 2001;Pany, 2002;Guo and Langley, 2003).
The number of intersection points and their vertical distribution changes significantly during the diurnal cycle of the satellite orbits (Fig. 3).The data presented in Fig. 7 belong therefore to a certain time (21:00 h in Fig. 3) and position but the effects discussed above show little variations with time.
The distribution in Fig. 7 shows somewhat more data at very low altitudes than could be expected from Fig. computed from any pair of rays with a vertical distance below 500 m.The "intersection point" has been placed in the middle of the connecting vector which results in a down shift of max.250 m.

Distribution of crossing angles
The crossing angle between intersecting ray paths gives a first measure of the angular range covered by the GPS data.Figure 8 shows a typical distribution of crossing angles with maxima near 90 • and very little data below ∼20 • (almost parallel rays) and above ∼160 • (almost antiparallel rays).The actual shape depends strongly on the satellite constellation but seems to be close to a Gaussian distribution.The elevation cut-off angle ε min and the density of ground stations have little effect on the angular distribution besides smoothing its shape.

Possible improvements
The quality and the amount of ground-based GPS data can substantially be improved in at least 4 ways.

Using more satellites
The small number of GPS satellites visible above a station (6-10) at a certain time restricts the possible views through the troposphere.Getting more satellites provides additional views and more information.The situation will therefore improve considerably until 2010: The European Galileo system will be operational at this time providing 30 satellites and the Russian GLONASS will be upgraded by new satellites and reference stations.About 90 satellites should be available in 2010, ∼30 of them visible to any station at any time.

Collecting data over a period of time
The ray paths to the fairly fast moving satellites traverse changing regions of the atmosphere.Collecting data over a limited period of time provides therefore little information about the temporal structures of the atmosphere but substantial spatial information.The resolution of the tomographic reconstruction can be improved by using data collected up to several hours.The reduced temporal resolution is of no consequence for several applications which require the data in certain intervals, e.g.data assimilation in numerical weather models.

Detecting satellites with low elevations
The amount and quality of data provided by GPS networks strongly depends on the elevation cut-off angle ε min used by the GPS receivers and the processing algorithms (Fig. 7).Reducing ε min would considerably increase not only the total amount of data but would also provide valuable information about the boundary layer.The latest development in the receiver hardware and processing techniques should lead to a decreasing ε min down to ∼2 • (Pany et al., 2001;Pany, 2002;Guo and Langley, 2003).In mountainous regions receivers can be placed on top of high mountains and detect even satellites with negative elevations.This leads to a well resolved boundary layer at least in certain directions (Zuffada et al., 1999;Aoyama et al., 2004).
Increasing the number of ground stations An increasing number of GPS ground stations would improve the data quality in several ways: The poor results at the outer boundary of the grids can be reduced by covering Ann.Geophys., 25, 1727Geophys., 25, -1734Geophys., 25, , 2007 www.ann-geophys.net/25/1727/2007/an extended area.A higher density of stations will provide fields with an enhanced resolution, especially in the vertical direction.The vertical resolution can be further improved by making use of the topographic conditions, e.g. at mountain sides.
Further improvements can be expected from low earth orbiters (LEOs) which provide horizontal views through the atmosphere.Several LEOs providing radio occultation data are currently in orbit, other missions are planned (Wickert et al., 2006;Gobiet and Kirchengast, 2004;Yunck, 2002).According to Foelsche (1999) up to 50 occultations above Europe could be expected within one hour if 20 LEOs where operating.Each occultation provides several nearly horizontal ray paths of about 100 km length through the atmosphere.These data carry valuable information about the vertical structure of the atmosphere but are not sufficient to guarantee a uniform vertical resolution all above Europe.

Conclusions
Atmosphere sounding by means of ground-based GPS can in principle provide the water vapour distribution in the troposphere with high temporal and spatial resolution.The spatial resolution is limited only by the density of the GPS receivers and the available navigation satellites.The horizontal resolution of the humidity fields would be restricted to the mean inter-station distance.The receiver density within a large GPS water vapour retrieval network will presumably be well below the scale hight of water vapour and it would therefore not be possible to reproduce the boundary layer very well using only GPS data.Additional meteorological information about the boundary layer would be required to obtain reliable vertical profiles.
The quantity as well as the quality of the GPS data depends strongly on the constellation of the GPS satellites and therefore on time.It has been shown that the number of intersection points and with it the available information changes significantly within rather short periods.It depends on the application how to deal with such poor constellations.If a high temporal resolution is required the network and the tomographic grid must be optimised.A high spatial resolution could be obtained by averaging over longer periods of time.
One of the most critical parameters is the elevation cutoff angle ε min used by the GPS processing.The number of available ray paths through the atmosphere increases heavily with decreasing ε min and provides more detailed information about the lower regions of the troposphere.
The results of this study cannot easily be applied to existing GPS networks.The exact position of the GPS receivers and their surrounding has a large impact on the information provided by a GPS data set.Furthermore, the absolute altitude of the GPS stations must be considered and should have a positive effect on the vertical resolution, especially in mountainous regions.A horizontal resolution comparable to the mean inter-station distance might be obtained for networks with rather homogeneously distributed stations.Randomly placed stations will reduce the overall resolution but some extra stations at optimized positions might enhance the situation considerably.The impact of the vertical resolution might be even larger.The realistic simulation of existing networks and their optimal extension is subject to future work.
The preconditions to GPS meteorology will improve considerably in the next few years as Galileo will start operating and GLONASS will be modernised.About 90 navigation satellites will be available multiplying the density of the GPS data.Future simulation studies will make use of the exact (computed) ephemerides of the new satellites and provide much more detailed information on the impact of additional satellites.

Fig. 1 .
Fig. 1.Idealised geometry used to estimate the altitude of intersection points.

Fig. 2 .
Fig. 2. Minimum altitude of intersecting ray paths as a function of the elevation ε min in the local horizon system and the inter-station distance S AB /2.

Fig. 4 .
Fig. 4. The number of intersection points increases with the number of ray paths according to ∼0.000137•N 2 ray .The ray paths have been computed for an increasing number of ground stations and two different satellite constellations with 29 (•) and 58 ( ) satellites.

Fig. 6 .
Fig. 6.Distribution of intersection points per gid cell within a 15×15×0.5 km grid.All distributions have been normalised to one.

Fig. 7 .
Fig. 7. Vertical distribution of intersection points for different elevation cut-off angles 0 • ≤ε min ≤10 • .N c is the number of intersection points in a 167 m vertical interval.

Fig. 8 .
Fig.8.Distribution of crossing angles between different intersecting ray paths for 15 km inter-station distances and ε min =0 • .N c gives the number of intersection points within a 4 • interval.A Gaussian fit has been applied to the data (dashed line).

Table 1 .
Total number of stations within a 2000×2000 km GPS network for inter-station distances between 100 km and 15 km.The number of rays is the mean value over a 12 h cycle considering all possible views (ε≥0 • ).
Diurnal variation of intersection points for different densities of GPS ground stations.A realistic GPS constellation with 29 satellites has been computed for a 12 h cycle (solid lines), a situation with 58 satellites has been simulated for 6 h (dashed lines).A 15×15×0.5 km grid has been used to calculate the number of intersection points.