New GNSS tomography of the atmosphere method – proposal and testing

Paper is focused on GNSS meteorology which is generally used for the determination of water vapour distribution in the atmosphere from GNSS measurements. Water vapour in the atmosphere is an important parameter which influences the state and development of the weather. At first, the paper presents basics of the GNSS meteorology and tomography of the atmosphere and subsequently introduces a new GNSS tomography method which doesn't require an extensive network of GNSS receivers, but uses only a few receivers situated in a line. After a theoretical concept describing this method and used mathematical background, the results from a real experiment are shown and discussed. Unfortunately the results indicate that presented method is not able to provide credible outputs. Possibly the main problem lies in an insufficient number of available signals from current global navigation satellite systems (GPS and GLONASS) where the improvement could be expected after the start of Galileo and Compass. Potential ways how to improve the results without increasing the number of satellites are outlined in the last section.


Introduction
Signal from the GNSS satellite is during its path through the atmosphere influenced by troposphere and ionosphere.The influence of the ionosphere is possible to eliminate from the GNSS measurements using a suitable method of dual-frequency data processing.The troposphere cannot be simply eliminated and, due to its dependence on the variable amount of water vapour in the atmosphere, it cannot be modelled with sufficient accuracy based on external data.However, it is possible to determine the troposphere and thus also the water vapour from GNSS measurements.
Currently used methods for atmospheric water vapour determination (e.g.meteorological radiosondes, remote sensing satellites, water vapour radiometer) are connected with some limitations and disadvantages and these could be completed by processing GNSS data.On the other side, even this method has its own problems which cause a limitation when used in numerical weather prediction (NWP) models.The main disadvantage is that it can produce only an absolute value of atmospheric water vapour in a zenith direction above the GNSS receiver in any time (Zenith Total Delay, ZTD).Thus, the vertical distribution of water vapour above the receiver, which is the information very important for data assimilation in NWP models, is not known.
From this reason the method of GNSS tomography, which would allow a three-dimensional study of atmospheric water vapour, has been developed since 2001 (Flores et al., 2001, Noguchi et al., 2001).At the present day there exists a few approaches how to solve this task, but all of them require a dense network of GNSS reference stations to achieve a reasonable spatial and temporal resolution.Those networks had to often been built artificially in the past for the needs of tomographic research projects.However, nowadays a possibility of using national reference station networks for creating tomographic systems over area of a whole country becomes an interesting topic.Both local and national tomographic systems require accessibility to the data from a large number of GNSS reference stations and necessary technical background for its processing.
Because all contemporary tomographic projects are connected with the necessity of dense network of GNSS receivers, authors attempted to propose a new easier method.Its core concept is to use only one cross-section throw the atmosphere instead of solving a whole threedimensional network of voxels.Thanks to that, the necessity of large number of receivers in a network would be reduced to a few receivers situated in a single line in a terrain.The output would be then limited to the 2D tomographic grid instead of the three-dimensional water vapour reconstruction.

GNSS meteorology
It has been shown many times that the GPS system is useful for tropospheric parameters estimation -the approach called GPS meteorology (Bevis et al., 1992, Duan et al., 1996).In the past only the signals from system GPS were used for the determination of tropospheric parameters.But nowadays a combination of signals from more global navigation satellite systems (GPS + GLONASS) can be used for such a processing.This approach is analogically called GNSS meteorology.
The standard value of the GNSS signal delay due to the atmosphere in the zenith direction is at zero sea level 2.3 m (or 8 ns).Troposphere influences GNSS signals in two different ways.Firstly, signal bends in the response to gradients of the index of atmospheric refraction, travelling a curved path instead of a straight line.Secondly, the signal travels throw an environment with some density slower than would travel in a vacuum.The total delay is the sum of those two effects.For the purposes of GNSS meteorology, however, ZTD can be divided in a different way -to the part caused by the hydrostatic components of the atmosphere (Zenith Hydrostatic Delay, ZHD) and part caused by the non-hydrostatic components of the atmosphere (Zenith Wet Delay, ZWD).While ZHD depends mainly on the atmospheric pressure, ZWD depends mainly on the water vapour (Zhengdong, 2004).During the GNSS data processing both components are modelled separately and the relation is expressed in (1).From a general point of view the hydrostatic characteristics of the atmosphere are rather stable in a time and it is relatively easy to model them.On the other side, the atmospheric Kačmařík, M., Rapant, L.: New GNSS tomography of the atmosphere method water vapour amount is temporally highly instable and thus much more difficult for modelling.
On the basis of Zenith Total Delays (ZTD) and meteorological parameters observed in situ of GNSS measurements it is possible to precisely determinate the atmospheric water vapour (Precipitable Water Vapour, PWV).

Tomography of the Atmosphere
For the modern generation of NWP models with a horizontal resolution of a few kilometres the existing network of sites launching meteorological radiosondes is not able to provide corresponding information about the vertical distribution of atmospheric water vapour.Since the classical GNSS meteorology cannot provide vertical profile, but GNSS tomography of the atmosphere can what has been already proven by many projects (e.g.Gradinarsky a Jarlemark, 2004, Champollion et al., 2004, Troller, 2004, Bender et al., 2011, Notarpietro et al., 2011).
The tomography method is generally known and with a success it is used in many directions of human activity as medicine, rock and material studies, archaeology etc.The principle is a spatial discretization of the space above the network of GNSS receivers into a three-dimensional network of cells called voxels.Information about the water vapour is reconstructed for each voxel from all accessible slant wet tropospheric delays (SWD) observed between GNSS receivers and satellites.Usually, it is assumed that within one voxel the water vapour is constant.Slant total delay of the GNSS signal caused by atmosphere can be written (Bevis et al., 1992): where Geoinformatics FCE CTU 9, 2012 ST D slant total delay N atmospheric refraction s length of the signal path In GNSS tomography slant total delays are provided by independent observations to all visible GNSS satellites and the refraction parameters are thus reconstructed from a large number of such delays.This defines an inverse problem, which could results in more than one solution.
The primary problem is a limited number of signals given with the number of ground receivers and with a time variable number of distant satellites.Using a discretization of space above the ground network of GNSS receivers into voxels, the relation ( 2) can be written where m vector of i observations, i.e. slant path delays (SWD) x vector of j refractivities (N) in j voxels A kernel matrix, the mapping of the state x on the observations m.The matrix elements a ij are the subsections of the i-th slant path in the j-th grid cell.(Bender et al., 2011) For the tomography purposes the path of signal is assumed to be a straight line and an influence of signal bending is neglected therefore.Thanks to this assumption a whole problem of inversion becomes linear and leads to a system of linear equations.

Input data
As above mentioned, the slant tropospheric delays in various azimuthal and elevation angles are used for the purpose of GNSS tomography and not just ZTDs which are estimated in classical GNSS processing model.More concretely, the SWDs caused by the non-hydrostatic characteristics of the atmosphere are used instead of STD.The relationship between ZWD and SWD is given as follows: where Mapping function defines a relationship between tropospheric delay in a zenith angle and those observed in all elevation angles.The basis of any mapping function is 1/ sin ε, which corresponds to the simplest case.Currently, the most used mapping function is defined by Niell (1996).
Horizontal gradients are considered to bring in the model the information about an azimuthal anisotropy in the atmosphere and in most accurate application these are usually estimated too.Post-fit residual represents a difference between a measured observation and an observation adjusted during the data processing.It is expected to contain the anisotropic part of the • ZTD computation with optional determination of horizontal gradients or post-fit residuals • ZHD computation using an appropriate hydrostatic model and values of meteorological quantities measured at the GNSS receiver's location, derivation of ZWD from ZTD and ZHD • Computation of SWD for particular observations using appropriate mapping function and optionally horizontal gradients or post-fit residuals with or without multipath correction

General Quality of Tomographic Reconstructions
Quality evaluation of GNSS tomographic results is possible to perform on the basis of comparison with vertical profiles from meteorological radiosondes or NWP models.Regardless of used ground network of GNSS receivers or applied tomographic solution so far, results are confronted with two common trends: • during standard atmospheric conditions, when the water vapour descends approximately exponentially with the elevation, the GNSS tomography achieves very good results comparable with radiosondes measurements, • in cases of atmospheric inversion or with distinctively different vertical development of water vapour from the standard atmospheric conditions, the tomography can represent this trend only with strong difficulties.
This situation is partially given by the necessity of adding constraints in the tomographic system defining mutual relations between neighbouring voxels.If the values between voxels are more smoothed, it is impossible to catch such an inverse trend.On the other way, if the Kačmařík, M., Rapant, L.: New GNSS tomography of the atmosphere method weights in the tomographic system are set more freely, the system becomes less stable and is more prone to provide unreliable results.

Theoretical concept of new tomographic method
As already mentioned in the introduction, the main idea of new tomography method is to use only a limited number of GNSS receivers situated in a line and reconstruct vertical 2D water vapour distribution.The most signals are crossing throw the 2D voxels situated above the centre receiver and therefore the quality of tomographic reconstruction in these voxels should be the best.Therefore the primary output from the method would be a vertical profile for a receiver situated in the middle of the line.Vertical profile is an output produced also by meteorological radiosondes, which is perfectly suitable for NWP assimilation requirements.
The orientation of used line of receivers in a terrain is very important for the new tomographic method.It directly influences the number of visible satellites in a thin atmospheric strip demarcated parallely with the line of receivers.
It is possible to use standard applications for planning GNSS measurements to propose a suitable orientation of the line and demonstrate an affordable state.Fig. 1 shows a GPS satellites constellation for an interval 14.30 -16.00 UTC for selected day at the reference station VSBO.This graphical output was acquired from Topcon Occupation Planning tool.The orientation of notional line was chosen in the east-west direction and simulation of obstructions was provided.Thus, only parts of the sky with azimuthal angles between 88 • -92 • and 268 • -272 • remained visible.Six various GPS satellites were accessible in the narrow strip and selected interval.Similar, but mostly worse situations held for different parts of a day or different line orientations.Therefore proposed method would be generally able to produce results only during a few intervals in a day depending on the actual satellite constellation.Also radiosondes measurements are limited with such a phenomenon while being launched each six or twelve hours.
Besides the maximum number of visible satellites, very important is their elevation.The lowest atmospheric layer above the ground contains the most of water vapour and therefore the signals at low elevations are of the highest interest for GNSS tomography.While choosing an ideal line orientation and time interval for the reconstruction, it is important to try to achieve a state with the highest number of visible satellites distributed on various elevations including the lowest ones.

Terrain measurements, data selection
On August 2nd 2012 an 8-hour long testing measurement was done at six GNSS receivers situated almost in a line with a northwest-southeast direction.Five of those receivers were GNSS reference stations (TVID, BISK, TKRN, VSBO, CFRM) while the sixth receiver (Topcon Hiper GD) was stabilized for the measurement's time only.All receivers worked in a GPS and GLONASS mode.GLONASS measurements were included to raise a number of satellites and increase an accessible reconstruction quality.
Total length of the line was 113 km and distances between two adjacent stations ranged from 21 to 25 km.The line orientation wasn t ideal but came out from the accessibility of GNSS measurements from selected reference stations and overall feasibility to realize such For processing GNSS measurements Bernese GPS Software (Dach et al., 2007) was used.Accurate coordinates of temporarily stabilized receiver were determined using a Precise Point Positioning technique (Zumberge, 1997).Consequently a double differencing technique was used for ZTD processing using settings summarized in a Table 1 and a Europe-wide network containing approximately 30 reference stations.
For ZHD determination a Saastamoinen model and measured values of meteorological parameters were used.Slant wet were computed from ZWD estimated in 30-second intervals from GPS and GLONASS observations.Horizontal gradients or post-fit residuals weren t used for SWD computation.Values of atmospheric air pressure and temperature were measured only at the station VSBO and in a place of the temporarily stabilized receiver.For locations of rest reference stations the values of atmospheric air pressure were re-counted using the Babinet s formula given in (5) and values of pressure and temperature measured at the location of the nearest receiver or at a close synoptic meteorological station operated by the Czech Hydrometeorological Institute.After SWD computations while horizontal gradients or post-fit residuals had not been added the observations from a narrow strip in a vertical direction above the line of receivers were selected.It was done in a geometrical way.A plane was held through a centre of the Earth and two ending points of the line.Plane's equation was computed from the coordinates of these three points in a geocentric coordinate system ITRF.A test whether a distance between a particular satellite and the given plane is in a concrete time less than 1 160 km for GPS satellite and less than 1 113 km for the GLONASS satellite was held.Values defining those distances were obtained by delimiting an angle of 2.5 • on a distance of 26 560 km which represents a distance between the centre of the Earth and a GPS satellite on its orbit.For GLONASS satellites the distance of 25 510 km was used.List of satellites satisfying given condition and therefore located in a defined strip was saved into a text file.Satellites from this file were consequently confronted with the real data from the RINEX observation files to exclude satellites satisfying given condition but being located on the opposite site of the Earth during the testing terrain campaign.Positions of all satellites in ITRF coordinate system in 15-minute interval were acquired from files with precise ephemerides used previously for ZTD processing.
Finally, an interval between 18.45 and 19.45 UTC was chosen to test the tomographic reconstruction.During this selected period one GPS satellite and five various GLONASS satellites were observed at elevation angles between 7 • and 168 • .Generally it would be better to select a constellation with more GPS satellites because there are still some difficulties with processing GLONASS measurements (mainly it is not possible to fix the ambiguities of GLONASS measurements and the quality of some input products like satellites ephemerides and clocks is still lower for GLONASS than for GPS).SWD values for GPS satellites observations should be therefore on higher quality level.Unfortunately according to measured data it was not possible to select any hour interval where more than one GPS satellite had been observed.

Mathematical background of the reconstruction
Tomographic reconstruction is primarily a mathematical task.In our case certain numerical tools were used to compute 2D tomographic grid.Beside slant wet delays the ZWD values for particular receivers were used as input.
First of all the discretize of tomographic system was done by dividing atmosphere above the stations into rectangular network where vertical lines are placed between the stations and horizontal lines are placed at desired height levels.Coordinates of the receivers for this purpose were used and also the vertical restriction in the height of 10 km where amount of water vapour was considered to be null.Signals that left the tomographic system below this height of 10 km were excluded from the reconstruction.Generally such signals can be included in the reconstruction only if it is possible to determine impact of water vapour on the part of the signal beyond the tomographic system.This operation requires an external source of information about the vertical water vapour distribution.Because of relatively high distance between two adjacent stations in the presented solution there were only a few signals which had to be excluded due to this problem.All of them were measured on low elevation angles on receivers situated by the line ends.
From the discretization and slant SWD measurements, a set of linear equations defined by equation ( 3) can be created.This gives a matrix where k is number of slant SWD measurements, l is number of horizontal layers and s is number of measuring stations.Value of each element of the matrix is determined by length of trajectory of slant signal passing through corresponding rectangle.Rectangles are numbered from the bottom and the left (ak1 correspond to bottom leftmost rectangle, ak(s+1) corresponds to the leftmost rectangle in the second row and so on).Right side of this system is created from the SWD of the corresponding measurement This gives us a system of equations from (3) from which we can compute amount of water vapours in each rectangle.Graphical illustration about the tomographic discretization and input data presents Fig. 2.
Kačmařík, M., Rapant, L.: New GNSS tomography of the atmosphere method However, in the computations there are more slant measurements than the number of rectangles in the discretization.This creates a problem of over determined system of equations.Such a system has no exact solution.Usual method for solving over stated systems of linear equations is the method of least squares (Björck, 1996).
The method of least squares solves the problem where b is not in the field matrix values of A, i.e. there is no such a vector x that solves the equation.The method of least squares tries to find an alternative solution xLS, which minimizes It can be rewritten equivalently with vector e, which modifies the right side, as and where (b + e) is in the field matrix values of A. This transforms system of linear equation into task of finding elements of vector e that are minimal in the Euclidean norm and satisfy previous condition.It can be done by various means such as QR or SVD decompositions and by iterative methods.In our solutions we used LSQR iterative method.After using the method of least squares we obtain approximate solution of linear equations system (3).

Results
Using computed SWD values and specified system the tomographic reconstruction was done.However its results stated in Table 2 can t represent a real water vapour distribution.The Kačmařík, M., Rapant, L.: New GNSS tomography of the atmosphere method table shows ZWD values reconstructed for particular cells.Their horizontal centres were situated to locations of GNSS receivers and used tomographic system was divided only into three layers vertically (thickness of two lower vertical layers was 2.5 km and 5 km for the highest layer).Generally tomographic projects use about ten vertical layers to represent the vertical distribution of water vapour in the troposphere.The classification into only three vertical layers would be insufficient for real meteorological applications and was chosen just to test the potential of the new tomographic method.If the results were promising a grid with much higher vertical resolution would be tested consequently.
Because during standard atmospheric conditions approximately 50 % of water vapour is situated in the surface layer of the atmosphere below 1.5 km from the ground and above 5 km there is only 5 % of its amount (Seidel, 2002)

Future insights
Potential improvement in described problem with low number of GNSS signals can be expected after start of Galileo and Compass systems when the number of satellites in space will be practically doubled.Except the increase in an absolute number of GNSS signals it should also theoretically bring more time intervals during the day, when higher number of satellites will be situated in a narrow strip above the line of receivers.
More detailed work at optimization and broaden of the mathematical apparatus used for the first testing purposes would also help.The usage of iterative component trying to find the most ideal reconstruction result based on the input data and a given initial state of the atmosphere would possibly help and should be implemented in a near future.The initial state will be derived from the standard atmospheric conditions when the water vapour decreases exponentially with the increasing height.Some tomographic projects use external source of data about the water vapour (e.g.Numerical Weather Prediction models, radiosonde profiles) to support the GNSS tomographic reconstruction but we do not consider doing that in future.It would bring new problems with the necessity of assimilation such a data in the solution and also made the suggested method dependent on external sources.

Conclusion
GNSS tomography of the atmosphere is a modern method allowing spatial and temporal study of the water vapour distribution in the atmosphere.So far presented projects were focused on three-dimensional reconstructions above large networks of GNSS receivers.A new tomographic method based on measurements from a line of a few receivers and reconstructing a 2D vertical grid above the line was proposed and described in the paper.Unfortunately the results based on described terrain measurements and mathematic solution are not too much promising.Situation indicates that GPS and GLONASS systems can t provide sufficient number of observations for this method.Authors expect a possible improvement after an optimization of mathematic apparatus, considering further constraints or applying a priori data and using Galileo and Compass systems in future.If it is proven that new method can provide credible information about vertical distribution of water vapour it could became an effective source of data for NWP models and possibly reduce the necessity of meteorological radiosondes.

Figure 2 :
Figure 2: Basic graphical illustration about the tomographic system Kačmařík, M., Rapant, L.: New GNSS tomography of the atmosphere method wet delay, but also various other unmodeled errors which influence the GNSS measurements: multipath, errors of antenna phase centre variations or signal noise.Some tomography projects use the tropospheric horizontal gradients and/or post-fit residuals as additional information of water vapour when reconstructing the anisotropic part of the atmosphere(Flores et al., 2001, Noguchi etal., 2001, Gradinarsky a Jarlemark, 2004, Bender et al., 2011), but other authors (Gradinarsky and Jarlemark, 2004, Nilsson et al., 2005) are sceptical of such an approach.Kačmařík et al. (2012) used various approaches for deriving SWD values from GNSS reference station GOPE and confronted them with in situ water vapour radiometer measurements.The comparison shown that using the tropospheric horizontal gradients or post-fit residuals (with or without correcting for multipath) did not give improved estimated SWDs.Adding the post-fit residuals to the GPS SWD had only a very small positive impact on bias, while it negatively influenced the standard deviation.The situation slightly improved after using the stacking maps for removing the multipath effect from the residuals.The results show that the post-fit residuals, even after multipath elimination, do not represent only the anisotropic part of the water vapour distribution but represent further insufficiencies in the model.

Table 1 :
Basic Characteristics of ZTD Processing in Bernese GPS SW

Table 2 :
, reconstructed ZWD with the highest values in the top layer can't represent a real state of the atmosphere prevailing during the test measurements.Tomographic reconstruction results, ZWD values in [cm] for particular cells.The bottom layer is in the lowest row of the table.But the conditions of convergence for successful run of the least squares method were not fulfilled what denotes an insufficient amount of input SWD data.Whole situation indicates that two contemporary accessible Global navigation satellite systems GPS and GLONASS are not able to provide a sufficient number of signals for a possibly meaningful usage of proposed tomographic method.