A method to improve the utilization of GNSS observation for water vapor tomography

Existing water vapor tomographic methods use Global Navigation Satellite System (GNSS) signals penetrating the entire research area while they do not consider signals passing through its sides. This leads to the decreasing use of observed satellite signals and allows for no signals crossing from the bottom or edge areas especially for those voxels in research areas of interest. Consequently, the accuracy of the tomographic results for the bottom of a research area, and the overall reconstructed accuracy do not reach their full potential. To solve this issue, an approach which uses GPS data with both signals that pass the side and top of a research area is proposed. The advantages of proposed approach include improving the utilization of existing GNSS observations and increasing the number of voxels crossed by satellite signals. One point should be noted that the proposed approach needs the support of radiosonde data inside the tomographic region. A tomographic experiment was implemented using observed GPS data from the Continuously Operating Reference System (CORS) Network of Zhejiang Province, China. The comparison of tomographic results with data from a radiosonde shows that the root mean square error (RMS), bias, mean absolute error (MAE), and standard deviation (SD) of the proposed approach are superior to those of the traditional method.


Introduction
Atmospheric water vapor plays an important role in the atmospheric dynamics.This parameter is one of the most variable components with strong impact on weather prediction, water resource management, and reducing natural hazards.To describe the process of water vapor changes, the tempospatial distribution of water vapor needs to be known in advance.Many studies have been conducted (Flores et al., 2000;Troller et al., 2002;Nilsson et al., 2006;Champollion et al., 2009;Perler et al., 2011;Chen and Liu, 2014;Adavi and Mashhadi-Hossainali, 2014), since the tomography technique was first published as having used GPS data (Bevis et al., 1992).
Using slant water vapor (SWV) acquired from the Global Navigation Satellite System (GNSS) observations, a water vapor field can be reconstructed by tomography.This technique requires the research area to be spatially discretized into a number of voxels and assumes that the water vapor in each voxel is constant in a given period.Water vapor density can then be estimated using SWV signals which penetrate the whole study area.However, due to the influence of geometric distribution of satellite constellations and the geographic distribution of ground-based receivers in the research area, many voxels, especially those at the bottom of the research area, may not be crossed by any signal ray during a certain period.This may lead to a numerical problem in the inversion of the tomography equation (Chen and Liu, 2014).The most commonly used method to solve this problem is to introduce constraints into the tomography equation, including a horizontal constraint, a vertical constraint, and a boundary Published by Copernicus Publications on behalf of the European Geosciences Union.

Y. B. Yao et al.:
A method to improve the utilization of GNSS observation constraint (Flores et al., 2000).However, the quality of the tomographic results may be anamorphic if the imposed constraints are not reasonable (Bender and Raabe 2007;Rohm and Bosy, 2009;Rohm, 2013).
To improve the quality of a tomographic result, and reduce the dependency on various constraints, the number of voxels crossed by signal rays should be increased (Bender and Raabe 2007;Rohm, 2012) by combining the observation of a multi-sensor system, increasing the number of receiver stations, or setting a lower satellite cut-off elevation angle.However, a majority of receivers which belong to the CORS network are not able to receive multi-sensor signals at the same time.In addition, increasing the density of receiver stations would increase construction and maintenance costs, and lowering the satellite cut-off elevation angle would make observations more susceptible to multipath error (Duan et al., 1996;Ware et al., 1997;Braun et al., 2001).Therefore, the aforementioned methods of increasing the number of voxels crossed by satellite signal rays have not been widely used.
It is a prerequisite that signal rays cross the entire research area when establishing a tomographic observation equation.As a result, signal rays that do not penetrate from the top of the research area are excluded; this reduces the utilization of existing data and means that many voxels are not passed through by any signal rays, especially those voxels comprising the bottom layers.To solve this problem, an approach of using signals penetrating from both the side and top of the research area is proposed.A satisfactory result was obtained using the data from CORS Network of Zhejiang Province, China.

Building the observation equation
For the tomography of a water vapor field, the most important observation focuses on SWV, which is relevant to the water vapor density and can be defined by where s represents the path of the satellite signal ray, ρ v is the water vapor density (units: g m −3 ).The SWV can be obtained by following method: where R ω = 461 (J kg −1 K −1 ) is the water : gas ratio constant; k 2 = 16.48K hPa −1 , k 3 = (3.776± 0.014) × 10 5 K 2 hPa −1 , are constants; T m is the weighted mean temperature; and SWD is the slant wet delay.The SWD can be extracted from the slant total delay by excluding slant hydrostatic delay (Davis et al., 1985;Bevis et al., 1994;Bender et al., 2011a).In our study, GAMIT10.5 software is used to process GPS observations (King and Bock, 2001;King and Bock, 2006).The sampling rate of the SWV data was 30 min while the interval of ZTD and horizontal gradient parameters estimated were 0.5 h and 2 h, separately.So, based on the Saastamoinen model (Saastamoinen, 1972) by which the zenith hydrostatic delay (ZHD) is calculated precisely using surface pressure, the ZWD can be easily separated from ZTD using formula ZWD = ZTD − ZHD.The SWD can thus be derived as (Flores et al., 2000) SWD ele,ϕ = f (ele) where ele and ϕ are satellite elevation and azimuth, respectively.f is the wet mapping function (wet Niell mapping function adopted in GAMIT10.5).G w NS and G w WE are components of the wet delay gradient in the north-south and eastwest directions.The last term, R e , refers to the unmodeled zero difference residuals.However, only the double difference residuals can be obtained from the GAMIT software, the zero difference residuals were calculated according to the approach proposed by Alber et al. (2000) in our study.
A linear equation relating SWV and water vapor density can be established: where SWV p is the slant water vapor amount of the pth signal ray (unit: mm); A p ij k is the distance of the pth signals ray inside voxel (i, j, k); and X ij k is the water vapor density in voxel (i, j, k) (unit: g m −3 ).The above equation can then be written in matrix form: where y is a column vector with a set of SWV measurements; m is the total number of SWV measurements; A is a coefficient matrix with elements representing distance; n is the number of voxels in the atmospheric region of interest; and x is a column vector of unknown water vapor density.In this study, an elevation cut-off angle of 10 • for signals at each receiver is selected so as to ignore the influence of ray bending (Mendes, 1999).

Constraint conditions
In Eq. ( 5), solving for the unknown water vapor density is in fact an issue for the inversion algorithm.The design matrix A in Eq. ( 5), however, is a large sparse matrix.The normal equation of A is singular which leads to numerical problems when using a direct inversion method (Bender et al., 2011a).
To overcome this rank deficiency issue, constraint conditions are often introduced to the tomography equation (Flores et al., 2000;Troller et al., 2002;Rohm and Bosy 2009;Bender et al., 2011a).In this study, both horizontal and vertical constraints are imposed.In the horizontal direction, the Gaussweighted functional method is used, while in the vertical direction, the functional relationship of the exponential distribution is established.After these constraints are imposed, the tomography model (traditional method) is then obtained as where H is the coefficient matrix of horizontal constraints equation whose elements are computed by Gauss-weighted function (Song et al., 2006), while V is the coefficient matrices for vertical constraints whose elements are obtained based on the exponential law for the water vapor refractivity (Davis, et al., 1993;Elósegui et al., 1998).

An approach for water vapor tomography using signals crossing both the side and top of the research area
Conventional methods merely use signal rays crossing from the top of the research area for water vapor tomography, however, signals rays such as P1, P2, P3, and P4 penetrating from the side of the research area are excluded, as shown in Fig. 1.This not only decreases the utilization of GNSS observations, but also makes many voxels, especially those at the bottom or edge areas in the research area suffer from an absence of crossing signals.Consequently, the accuracy of tomographic results at the bottom is influenced; however, water vapor is mainly located in the bottom of the research area.
To solve this issue, an approach is proposed that allows all signals to be used for calculating the initial information constraint, including those penetrating from the side of the research area as well, and then the established initial constraint is imposed to tomographic model for water vapor tomography.The water vapor unit index is introduced based on the fact that many signal rays which penetrate from the sides of a research area partly cross it.Then, the contribution of the part of the signals which are located within the research area contributes to the final tomographic result.Both the radiosonde data obtained in the first 3 days and the signals that cross both from the side and top of the area for the tomographic period were utilized; the initial water vapor density value of the voxels was calculated and imposed to the tomographic equation as an initial information constraint.The advantage of the proposed approach lies in the effective use of observed information, like satellite signal rays P1 to P4 , which penetrate from the gray area, as shown in Fig. 2, can also be used.The use of satellite signals, therefore, is considerably improved and the number of voxels crossed by signals is also increased.Using the proposed method, the quality of the reconstructed water vapor field is found to be enhanced  at the bottom area and the quality of the entire tomographic area is also improved.

Calculation of initial water vapor density
The voxels in which the radiosonde is located, and those in the vertical direction, are called datum voxel, as shown by green voxels in Fig. 3a and b.The number of datum voxels is equal to the number of layers within the tomographic area.Assuming there is a water vapor unit index for every SWV signal in each voxel, by which the water vapor density can then be reflected, the water vapor unit index may be defined as the unit water vapor density of a slant path in a certain voxel and can be expressed as where is initial water vapor density value of datum voxel (i, j, k), r ij k is water vapor unit index in datum voxel (i, j, k), and SWV is the integrated slant water vapor.Equation (7) indicates that determining the water vapor unit index for SWV signals in each voxel is a prerequisite to calculating initial water vapor values.A water vapor unit index model can be built using water vapor unit indices of datum voxels assuming there is a strong spatial correlation among water vapor density values.Thereby, the water vapor unit index in non-datum voxels traversed by signals can then be replaced using the established model.1. Determining the water vapor unit index for each SWV signal in datum voxels using SWVs which only penetrate from the head of tomographic array and radiosonde data obtained in the first 3 days.
As shown in Fig. 3b, a signal ray p in station q penetrates the datum voxels.In Fig. 3c, the signal ray p crosses one of those datum voxels (i, j, k), and its distance in datum voxel (i, j, k) is d pq ij k .So, the water vapor unit index in datum voxel (i, j, k) for signal ray p in station q can be obtained from where γ pq ij k is a water vapor unit index of datum voxel (i, j, k) for signal ray p in station q; d pq ij k is the distance of signal ray p of station q in datum voxel (i, j, k); SWV pq is the total water vapor amount of signal ray p of station q; and SWV pq ij k is the amount of water vapor of signal ray p of station q in datum voxel (i, j, k), which can be expressed by where ρ RS ij k is the average water vapor density calculated based on radiosonde data in the first 3 days that in datum voxel (i, j, k).The water vapor unit index γ pq ij k can then be obtained by combining Eqs. ( 8) and (9): (10) 2. Establishing water vapor unit index model for every datum voxel.
In the analysis, the water vapor unit index is a function of elevation, so the water vapor unit index model for each layer can be established based on those water vapor unit indices of each datum voxel calculated in Step (1).The water vapor unit index model γ k ele can thus be expressed as where ele is the satellite elevation angle, and γ k ele is a water vapor unit index of the kth layer for the SWV signal with elevation ele.The water vapor unit index models are related to two different coefficients a1 and a2, so, two or more signal rays must cross the datum voxel if the model is to be established.

Calculating water vapor unit index for the non-datum
voxels in which signal rays crossed using the established water vapor unit index model in Step (2).Here, the water vapor unit index of the non-datum voxels is obtained based on the assumption that we regard the established model for every datum voxel in step (2) as the model for the whole layer which the datum voxel located.
4. Finally, the average initial water vapor value of the nondatum voxel (m, n, k) can be attained using the SWV signals which cross from both the side and top of the research area, and the water vapor unit index calculated in Step (3), the formula expresses as follows: where ρ initial mnk is the average initial water vapor value in non-datum voxel (m, n, k), γ pq mnk is the water vapor unit index of signal ray p of station q in non-datum voxel (m, n, k), q nmk is the number of stations used for water vapor tomography, pq mnk is the number of signal rays of station q, which penetrate the non-datum voxel (m, n, k), and n mnk is the total number of signal rays which cross the non-datum voxel (m, n, k).

Constructing constraint equation of initial water vapor density
The initial water vapor density values calculated in Sect.3.1 are regarded as an initial constraint equation and can be expressed as follows: where l is the total number of calculated initial water vapor values.Equation ( 13) can be written in matrix form: Consequently, the tomographic model of the proposed method can be realized by combining Eqs. ( 6) and ( 14) as follows: where the last row is an extension of Eq. ( 14).The initial elements which could not be estimated in Eq. ( 13) would be replaced by zeros in C m×n and ρ initial m×1 leading to rows with 0 • x i = 0.
In this study, the multiplicative algebraic reconstruction technique (MART) is used to solve the proposed tomographic model (Eqs.6 and 14) because of its rapid convergence (Bender et al., 2011b).In our experiment, the initial field of MART is calculated using the method proposed by Flores et al. (2000) to the tomographic model, and then the MART method is used to get the final tomographic result (Chen and Liu, 2014).The number of iteration of MART is obtained based on the root mean square (RMS) error between the observed SWV and reconstructed SWV and it will converge rapidly after several iterations.It should be pointed out that for the voxels with signals crossed, a more accurate tomographic result can be obtained while for the voxels without signals penetrated, a smoother final field would be calculated due to the fact that horizontal and vertical constraints distribute the information within those grids.

Accuracy analysis
Radiosonde data can provide accurate vertical water vapor information at various altitudes and often used as a reference to evaluate the quality of the water vapor field obtained from other methods (Niell et al., 2001;Adeyemi and Joerg, 2012;Liu et al., 2013), so the values derived from radiosondes are used to validate the proposed method.
For any water vapor tomography model, the accuracy of the tomographic result is key to evaluate its quality.Here, RMS, bias, MAE, and SD are used for this purpose (Rohm and Bosy, 2009;Shangguan et al., 2013).Water vapor density values derived from different tomographic models are compared with that from a radiosonde using the following equations (Guerova, 2003;Yao et al., 2013): Here, x i is the water vapor density value estimated by different tomographic methods, x 0 i is the water vapor density derived from a radiosonde, which is considered as a true value, and N is the number of sample elements.RMS is often used to verify the reliability of the proposed method while MAE can be used to evaluate the extent of any deviation between the estimated and true values.

Validation of the new method
Here, tomographic experiments are carried out using two different tomography models: -Method 1: using a tomography model based on the traditional method, as shown by Eq. ( 6); -Method 2: Using a tomography model based on the proposed method, as indicated by combining Eqs. ( 6) and ( 14).

Tomography strategy
The tomographic experiment is implemented using GPS data from 10 stations (as shown by in Fig. 4) from the CORS network of Zhejiang Province, China, for 31 days (1 to 31 May 2015).The period covered is 0.5 h for each step of tomographic solution and results are compared with the water vapor density derived from radiosonde data.As shown in Fig. 4, radiosonde station 58457 is located in the research area, where radiosonde balloons are launched twice daily at 00:00 and 12:00 UTC.
As mentioned before that radiosonde data is one of the most accurate means to obtain vertical water vapor profiles  which can reflect the water vapor distribution in atmosphere.Therefore, water vapor profile for different altitudes is calculated by interpolation method for every layer based on the exponential law proposed by Davis et al. (1993).Using the radiosonde data for 10 years from 2004 to 2014 at specific dates at 00:00 and 12:00 UTC, the average water vapor profile and SD are obtained for various altitudes (see Fig. 5). Figure 5 shows that the water vapor density and SD are both close to zero above 10 km.Thereby, the vertical boundary has been selected as the 10 km level surface, assuming no water vapor is above this altitude.The range of the research area is as follows: latitude 29.95 to 30.63 • N and longitude 119.95 to 120.70 • E, the vertical resolution is 0.8 km, and the total number of voxels is 4 × 5 × 13.
Initially, the average number of daily signals used and the average number of daily voxels in which the signal rays crossed for each 30 min from 1 to 31 May are analyzed for two methods.Figure 6a gives the average number of signals used and Fig. 6b shows the average number of voxels crossed by these signals.Figure 6a shows that the number of signals used increased by adding signals penetrating from the side of the research area, which improves the utilization of observations and allows a more accurate tomographic result.In addition, Fig. 6b shows that the number of voxels crossed by signals also increased by 7.38 %, from 56 to 63.38 %, according to the statistical analysis over a 31-day period.
In addition, the distribution of effective SWV signals used for water vapor tomography is analyzed for two methods over 31 days in May 2015. Figure 7 shows the average number of SWV signals used for various elevation angle ranges: as shown, whatever range of elevation angle, the number of effective SWVs for Method 2 is greater than that of Method 1.This outcome suggests that compared to the traditional method, the utilization of SWV signals has been improved significantly by using the proposed approach.

Accuracy analysis of the water vapor unit index model
To test the reliability of the water vapor unit index model, the data collected over 31 days are analyzed.The initial water vapor density value is calculated using the established model, and compared with water vapor density derived from radiosonde data at 00:00 and 12:00 UTC. Figure 8a shows the calculated RMS and Fig. 8b shows the number of initial water vapor density values calculated by the proposed method.
The average RMS (1.25 g m −3 ) calculated over the 31-day period proved that the established model was accurate.In addition, from Fig. 8b, it is concluded that more than half of the initial water vapor values (the average number is 165 of 260 for 31 days) can be obtained.

Validation and accuracy analysis
To validate the proposed method, tomographic results from different methods at 00:00 and 12:00 UTC for 31 days (1 to 31 May 2015) are analyzed.Water vapor densities at the radiosonde station are first calculated using voxel information obtained from different tomographic methods: RMS, bias, MAE, and SD are then obtained (Figs. 9 to 12), and the statistical results for these 31 days are summarized in Table 1.As seen from Figs. 9 to 12, the RMS, bias, MAE, and SD using Method 2 are smaller than those when using Method 1. Table 1 shows that, in terms of RMS, bias, MAE, and SD, the values from Method 2 (1.29, −0.23, 0.91, and 1.10 g m −3 ) are lower than those from Method 1 (1.61, −0.47, 1.10, and 1.38 g m −3 ).Compared to the traditional method, the proposed method is proved more consistent with radiosonde data: the proposed method is, therefore, both feasible and effective.
To compare directly the vertical accuracy of water vapor density derived from different altitudes, the tomographic results obtained using different tomographic methods (1 to  31 May 2015) are analyzed.Figure 13 shows the RMS and relative error change with altitudes.It can be observed in Fig. 13 that the RMS and relative error of method 2 is less than that of method 1 for different layers, especially at the lower layers, which is clear evidence that the proposed approach improves the accuracy of the final tomographic result.
In addition, it also can be seen from Fig. 13a that the RMS value of different methods is decreased with height.In contrast, the relative error in general decreases with height and then increases above 4.5 km.The maximal values of relative error are found at both the top and bottom layers of tomographic region, this is because in the upper layers the value of water vapor field is relatively low, while in the lower lay- ers the difference between the radiosonde and tomographic result is relatively large.Therefore, the conditions of even a relatively small discrepancy in upper layers and a relatively large water vapor density in each voxel in lower layers would result in a relatively large error.
In addition, the water vapor density profiles for different altitudes at individual dates are given in Fig. 14.Two dates (12:00 UTC 19 May 2015 and 00:00 UTC 25 May 2015) are selected for they correspond to the maximum and minimum RMS during the experiment period of 31 days.Figure 14 shows that the tomographic water vapor profiles of both method 1 and 2 have a good agreement with that from radiosonde data.It is clear that the water vapor density profile of Method 2 better matches that from radiosonde data than Method 1, especially in the layers of 1 to 5, which implies that the water vapor density derived from the proposed method is superior to that of the traditional method.
The samples (water vapor density value for different layers at different dates) are selected randomly from the tomographic results from the two methods over the trial period.The scatter plot of water vapor density and RMS using different methods is compared with that derived from radiosonde data (see Fig. 15).According to data from 403 samples (data sampled at date 12:00 UTC once per day for the experiment period from 1 to 31 May 2015), it is concluded that the tomographic quality of the water vapor field is improved by using the proposed method; the RMS has decreased from 1.60 to 1.28 g m −3 and the slope of the regression profile has improved from 0.78 to 0.84.In addition, Fig. 15 shows that the distribution of water vapor density in Method 2 is closer to the regression profile.

Conclusions
An approach has been proposed that uses signal rays penetrating both from the side and top of the research area.In this approach, satellite signals captured by ground-based receivers are fully utilized, which means that satellite signals that do not cross the entire research area are not wasted.This proposed approach improves the utilization of observed data and increases the number of voxels crossed by satellite signals.
The proposed approach is validated by tomographic experiments using observed data from the CORS network of Zhejiang Province, China for 31 days from 1 to 31 May 2015.The experimental results verified that the proposed approach is feasible and effective.Through comparison of the tomographic results with radiosonde data, results show that the RMS, SD, bias, and MAE of proposed method are 1.29, 1.10, −0.23 and 0.91 g m −3 , respectively, which are much smaller than that of traditional method (1.61, 1.38, −0.47 and 1.10 g m −3 , respectively).Comparing the water vapor density profile at different altitudes, we find that the profile of proposed method better matches with radiosonde, especially for the height of 1 to 4 km in the tomographic area.This implies the tomographic accuracy of a water vapor field had been improved, especially with regard to data in the bottom layers of the research area.

Figure 2 .
Figure 2. Plane distribution of satellite signals.

Figure 4 .
Figure 4. Distribution of receiver stations and radiosonde station.

Figure 5 .
Figure 5. Average water vapor density and SD derived from the 58457 radiosonde for the 10 years from 2004 to 2014.

Figure 6 .
Figure 6.The average number of daily signals used and the average number of daily voxels crossed by signals once per day for different methods for 31 days from 1 to 31 May 2015.

Figure 7 .
Figure 7. Distribution of effective SWV signals with elevation angle for different methods.

Figure 8 .
Figure 8. RMS for the water vapor unit index model and the number of initial values calculated using the established model for the 31 days from 1 to 31 May 2015.

Figure 9 .
Figure 9. RMS comparison for two methods from 1 to 31 May 2015.

Figure 10 .
Figure 10.Bias comparison of two methods from 1 to 31 May 2015.

Figure 11 .
Figure 11.MAE comparison of two methods from 1 to 31 May 2015.

Figure 13 .
Figure13.RMS and relative error change with heights (blue curve and red curve are derived from the differences between the profiles of method 1, method 2 and radiosonde, separately for 62 epochs spanning from 1 to 31 May 2015.

Figure 15 .
Figure 15.Scatter plot of water vapor density between different methods and radiosonde at date 12:00 UTC once per day from 1 to 31 May 2015.
Figure 12.SD comparison for two methods from 1 to 31 May 2015.