the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An improved pixelbased water vapor tomography model
Yibin Yao
Linyang Xin
Qingzhi Zhao
As an innovative use of Global Navigation Satellite System (GNSS), the GNSS water vapor tomography technique shows great potential in monitoring threedimensional water vapor variation. Most of the previous studies employ the pixelbased method, i.e., dividing the troposphere space into finite voxels and considering water vapor in each voxel as constant. However, this method cannot reflect the variations in voxels and breaks the continuity of the troposphere. Moreover, in the pixelbased method, each voxel needs a parameter to represent the water vapor density, which means that huge numbers of parameters are needed to represent the water vapor field when the interested area is large and/or the expected resolution is high. In order to overcome the abovementioned problems, in this study, we propose an improved pixelbased water vapor tomography model, which uses layered optimal polynomial functions obtained from the European Centre for MediumRange Weather Forecasts (ECMWF) by adaptive training for water vapor retrieval. Tomography experiments were carried out using the GNSS data collected from the Hong Kong Satellite Positioning Reference Station Network (SatRef) from 25 March to 25 April 2014 under different scenarios. The tomographic results are compared to the ECMWF data and validated by the radiosonde. Results show that the new model outperforms the traditional one by reducing the rootmeansquare error (RMSE), and this improvement is more pronounced, at 5.88 % in voxels without the penetration of GNSS rays. The improved model also has advantages in more convenient expression.
As the most active component in the troposphere, water vapor is one of the most difficult parameters to monitor and describe (Rocken et al., 1997). A good understanding of the spatiotemporal variation of water vapor is very helpful for improving weather forecasting and early warning of disastrous weather (Weckwerth et al., 2004).
The Global Navigation Satellite System (GNSS) technique can not only retrieve the precipitable water vapor (Bevis et al., 1994; Emardson et al., 1998; Baltink et al., 2002; Bock et al., 2005) but can also monitor the threedimensional water vapor distribution by using the GNSS tomography method (Flores et al., 2000; Seko et al., 2000; Macdonald et al., 2002).
Braun et al. (1999) first proposed the concept of reconstructing the tropospheric water vapor structure using 20 GPS stations in a regional observational network. Flores et al. (2000) first applied the tomography technique to obtain wet refractivity from the GNSS slant wet delay (SWD). In the same year, Hirahara (2000) used a different method to conduct GNSS tomography experiments, which also successfully obtained threedimensional water vapor fields. Since then, many scientists proposed new methods and applied them to GNSS water vapor tomography experiments (Rohm et al., 2014; Yao et al., 2016; Zhang et al., 2017; Ding et al., 2018; Zhao et al., 2018).
Hirahara (2000) conducted a fourdimensional tomography experiment and solved the tomography equations using the damped leastsquares method. Braun et al. (2003) and Braun (2004) overcame the sensitivity problem in GNSS tomography by using the extended sequential filtering method. Perler et al. (2011) presented a new parametric method for the water vapor retrieval. Nilsson and Gradinarsky (2006) obtained the wet refractivity directly from the GNSSphase observations using the Kalman filter method. Rohm and Bosy (2009) used the Moore–Penrose pseudoinverse of variance and covariance to solve the linear equations and emphasized the illposed tomography equation. Yao et al. (2016) obtained good tomographic results by using the optimal gridmaking method. Zhao and Yao (2017) proposed a method of using the sidepenetrating signals for tomography and improved the utilization rate of the GNSS rays. Aghajany and Amerian (2017) obtained the tomography results of water vapor profiles from ERAI numerical weather prediction data by applying a 3D ray tracing technique. Dong and Jin (2018) reconstructed the water vapor density using multiGNSS systems and showed that the accuracy of GNSS tomography results are improved by 5 % from the GPSonly system to the dual systems (GPS + GLONASS). Besides, the virtual reference station approach (Marel, 1998; Vollath et al., 2013), an effective method for attenuating the effects of atmospheric errors in longdistance dynamic positioning, was also used in GNSS tropospheric tomography.
In previous studies, since most of the GNSS tomography methods divided the troposphere of interest into finite voxels, and the water vapor density in each voxel is considered as constant, these methods with the above assumptions are defined as the pixelbased method. Apparently, this kind of method cannot retrieve the variations in voxels and breaks the continuous nature of the troposphere as well. Moreover, the pixelbased method requires each voxel to have a parameter to represent the water vapor density in it, which may lead to the situation that we have to use huge numbers of parameters when the research area is large and the expected resolution is high. Last, overparameterization may cause mathematical problems when we use limited observations to invert for the parameters that may be correlated. Therefore, this paper analyzes the limitations of the traditional pixelbased water vapor tomography method and proposes an improved model. The improved model uses the water vapor density obtained from the traditional model as the input value and outputs the fitting water vapor density by the layered optimal polynomial functions. This new model has the advantage of reflecting the variations in voxels and keeping the continuity of water vapor in the troposphere.
2.1 Establishment of the traditional pixelbased water vapor tomography model
2.1.1 Retrieval of SWV
For tropospheric tomography, the most important observation is the slant water vapor (SWV), which is related to the water vapor density and can be defined by
where s represents the path of the satellite signal ray, and ρ_{V} is the water vapor density (units: g m^{−3}). SWV can be obtained by the following method:
where ${k}_{\mathrm{2}}^{\prime}=\mathrm{16.48}$ K hPa^{−1}, ${k}_{\mathrm{3}}=\mathrm{3.776}\times {\mathrm{10}}^{\mathrm{5}}$ K^{2} hPa^{−1}, and R_{ω}=461 J kg^{−1} K^{−1}, which represent the specific gas constants for water vapor. T_{m} is the weighted mean tropospheric temperature, calculated from an empirical equation proposed by Liu et al. (2001) using the meteorological measurements. SWD is the slant wet delay, which may be given as
where elv is the satellite elevation, φ is the azimuth, m_{wet} is the wet mapping function, and ${G}_{\mathrm{NS}}^{\mathrm{w}}$ and ${G}_{\mathrm{EW}}^{\mathrm{w}}$ are the wet delay gradient parameters in the north–south and east–west directions, respectively. R refers to the nonmodeled zero difference residuals that may involve unmodeled influence on the threedimensional spatial water vapor distribution, which can make up for the lack of tropospheric anisotropy using only the gradient term (Bi et al., 2006). Since the GAMIT software only provides the double difference residuals, the zero difference residuals in this paper are obtained from the double difference residuals according to the method proposed by Alber et al. (2000). ZWD is the zenith wet delay, which is extracted from the zenith tropospheric delay (ZTD) by separating the zenith hydrostatic delay (ZHD) using equation ZWD = ZTD − ZHD. ZHD can be calculated precisely using surface pressure based on the Saastamoinen model (Saastamoinen, 1972):
where P_{s} is the surface pressure (unit: hPa), φ is the latitude of the station, and H is the geodetic height (unit: km). The unit of ZHD is meter.
Since the SWV is obtained, the tomographic area can be discretized into a number of voxels in which the water vapor density is a constant during a given period of time. Therefore, a linear equation relating the SWV and the water vapor density can be established as follows (Chen and Liu, 2014):
where SWV^{p} is the slant water vapor of ρth signal path (unit: mm). i, j, and k are the positions of discrete tomographic voxels in the longitudinal, latitudinal, and vertical directions, respectively. ${D}_{ijk}^{\mathrm{p}}$ is the distance of the ρth signal in voxel (i, j, k; unit – km). ρ_{ijk} is the water vapor density in a given voxel (i, j, k; unit – g m^{−3}). A matrix form of this observation equation can be rewritten as follows (Flores et al., 2000; Chen and Liu, 2014):
where m is the number of total SWVs, and n is the number of voxels in the tomographic area. y is the observed value here as the SWV, which penetrates the whole interest area, A is the coefficient matrix of the signal transit distances through the voxels, and ρ is the column vector of the unknown water vapor density.
2.1.2 Constraint equations of the tomography modeling
Solving for the unknown water vapor density in Eq. (6) is actually an inversion algorithm issue, as the design matrix A is a large sparse matrix whose normal equation is singular, leading to numerical problems when using a direct inversion method (Bender et al., 2011). To overcome this rank deficiency problem, constraint equations are often introduced into the tomography equation (Flores et al., 2000; Troller et al., 2002; Rohm and Bosy, 2009; Bender et al., 2011). In our study, the horizontal constraint equation is imposed by the Gaussian weighted functional method (Guo et al., 2016), and the vertical constraint equation is imposed by the functional relationship of the exponential distribution (Cao, 2012), respectively. The final tomography model is then obtained as
where H and V are the coefficient matrices of horizontal and vertical constrains, respectively. In order to obtain the inverse matrix shown in Eq. (7), singular value decomposition is used in this paper (Flores et al., 2000).
2.2 An improved pixelbased water vapor tomography model
The improved pixelbased water vapor tomography model proposed in this paper can take advantage of facilitating the continuity of the water vapor expression efficiently in spatiotemporal distribution and calculating the water vapor density conveniently. The improved tomography model firstly obtains the water vapor density from voxels penetrated by GNSS rays using the traditional pixelbased tomography model then obtains the optimal polynomial function of each layer through adaptive training. With known coefficients of the layered optimal polynomial functions, the water vapor density can finally be calculated given the latitude, longitude and the altitude. Specific steps are as follows.
First, use the traditional pixelbased water vapor tomography model to obtain the initial water vapor density from voxels penetrated by GNSS rays as the observation values for obtaining the optimal polynomial function coefficients of each layer.
Second, normalize the coordinates of each voxel center in the tomographic area, since the polynomial fitting of the water vapor at each tomographic layer is, in essence, for establishing the relationship between the latitude as well as the longitude of the tomographic region and the water vapor density. The general expression is:
where B is the latitude, L is the longitude, and V_{d} represents the water vapor density. Polynomial coefficients such as a_{i} are obtained via the leastsquares method. In the process of data solving, where the numerical values of the latitude and longitude may not be small, then the magnitude of multiple power may be larger than 10^{4}, which will lead to the illposed problem of the design matrix in the inversion process and eventually affect the reliability of the estimated coefficients. To ensure the design matrix will be relatively stable in the inversion process, the latitude and longitude coordinates B and L need to be normalized. The specific methods are as follows:
where B^{∗} and L^{∗} are the normalized latitude and longitude, respectively, and B and L are the latitude and longitude in the initial region range. μ is the average value of the latitude or longitude, and σ is the standard deviation of the latitude or longitude.
Third, determine the layered optimal polynomial functions of the improved model through adaptive training.

First, based on the size of the selected tomographic region, determine the highest polynomial fit order. In this paper, the highest polynomial fit order, chosen as 5, turns out to be generally sufficient.

Then set the water vapor density from voxels penetrated by GNSSsignal rays as the input value and keep trying out new polynomial functions, as the optimal polynomial function of each layer is obtained by adaptive training. What needs to be noted here is that the number of estimated coefficients needs to be less than that of the voxels penetrated by GNSS rays in each layer. Otherwise the overfitting problem will happen.

Finally, after comparing training results of multigroup polynomial functions at different levels, the polynomial function with the minimum rootmeansquare error (RMSE) value obtained from the water vapor density of the postfitting layer and that of the European Centre for MediumRange Weather Forecasts (ECMWF) is the best fitting equation for this layer. Each layer could have the individual optimal polynomial function in general.
Fourth, after finding the optimal polynomial function of each layer in different heights, using the latitude, longitude, and altitude information in the function could obtain the threedimensional water vapor distribution in the tomographic region. The continuous water vapor density can be easily described by broadcasting the estimated coefficients of the layered optimal polynomial functions.
2.3 The optimal polynomial selection based on adaptive training
Since the polynomial form can reflect the continuity of water vapor well and has the advantage of highefficiency computing as well as easy expression, this paper chooses the polynomial form as the layered fitting function. Based on adaptive training, the selection process of the layered optimal polynomial function is as follows.
First, construct a polynomial equation training library, which contains a wide variety of polynomial function forms of the latitude and longitude as independent variables while using the water vapor density in the voxels as the dependent variable. After many experiments, the maximum power of the latitude and longitude, found as 5, is sufficient in describing the water vapor changes. Therefore, the maximum power of the fitting function part is adopted as 5 in the training library.
Second, according to the water vapor density observations from the voxels penetrated by the GNSS signals at each level, the form of the candidate polynomial function of each layer is automatically determined from the polynomial function training library, ensuring that the number of observations at all levels is always greater than that of estimated coefficients of candidate polynomials.
Third, calculate the water vapor variation index (WVVI) of each layer in both east–west and north–south directions using the traditional water vapor tomography results as shown in Eq. (10):
where wv_{EW} and wv_{NS} are the water vapor density in east–west and north–south directions, respectively.
The WVVI is a changing rate indicator of the water vapor density in a given direction. According to the water vapor variation index of each layer in the east–west and north–south direction, whether the water vapor exists mainly in the east–west distribution or the north–south distribution can be determined. As an aid, WVVI can decide the main body of the alternative polynomial function with higher order, whether in longitude or latitude, and then efficiently find the layered optimal polynomial function. If the water vapor density of a layer indicates a horizontal gradient of east–west distribution, the polynomial function with higherorder term of the longitude should be given the priority. It suggests that when the water vapor shows an east–west gradient distribution, there is a better correlation between the longitude and the water vapor variation; furthermore the highorder term in longitude can better reflect the nuanced water vapor variation. A simple example of the polynomial function with a higherorder term in longitude is shown in Eq. (11):
Otherwise, when the water vapor density of a layer indicates a horizontal gradient of north–south distribution, the polynomial function with higherorder term of the latitude should be given priority. A simple example is shown in Eq. (12):
While the distribution regularities of the water vapor density gradient are not clear or obvious, the polynomial function with the same order of the latitude and longitude can be considered as the example shown in Eq. (13):
Fourth, the candidate polynomials of all levels screened by the WVVI gradient auxiliary information are used as the next comparative polynomials, and the required estimated coefficients of the comparative polynomial are solved according to the principle of least squares through Eq. (14) and automatically recorded into the coefficients data set. M is the matrix of the longitude and latitude, and the vector x comprises the unknown coefficients of the comparative polynomial functions as shown in Eq. (15):
Fifth, through the comparative polynomials with the estimated coefficients in each layer, the wholevoxel water vapor fitting of each layer is automatically fit with the information of the latitude and longitude. In order to obtain the RMSE, the fitting result would be compared with the ECMWF water vapor density of each layer in this period. The results are then saved to the accuracy data sets of each layer. The comparative polynomials with the estimated coefficients are constantly selected to train the fitting of the layered water vapor density and are then compared with the water vapor density of the ECMWF at each layer. Thus, large accuracy data sets of RMSE can be obtained where the smallest RMSE value of the comparative polynomial form can be chosen, and then the optimal polynomial of each layer could come into being. It is noteworthy that the optimal polynomial of each layer might be different. With the layered optimal polynomial, the continuous threedimensional water vapor density in the tomographic region can be expressed conveniently by transmitting the estimated coefficient's information.
3.1 Experimental description and dataprocessing strategy
To study whether the accuracy and stability of the improved water vapor tomography model are better than those of the traditional model, the following experiment is designed.
Tomographic data are obtained from the Hong Kong Satellite Positioning Reference Station Network (SatRef) from 25 March to 25 April 2014. Two epochs are taken each day (00:00 and 12:00 UTC). The corresponding meteorological data are also used to calculate the precipitable water vapor (PWV). The tomographic area ranges between latitude 22.24 to 22.54^{∘} N and longitude 113.87 to 114.29^{∘} E. Taking the mean sea level as the height of the base level, the vertical resolution is 0.8 km, and the total grid number is $\mathrm{5}\times \mathrm{7}\times \mathrm{13}$. In the selected area, a total of 11 GNSS stations and one radiosonde station (located at King's Park, Hong Kong) are selected, and the ECMWF grid data are extracted twice daily, at 00:00 and 12:00 UTC, from 25 March to 25 April 2014 (grid resolution of $\mathrm{0.125}{}^{\circ}\phantom{\rule{0.125em}{0ex}}\times \phantom{\rule{0.125em}{0ex}}\mathrm{0.125}{}^{\circ}$). See Fig. 1 for details.
According to the official website of the Hong Kong Observatory (http://www.weather.gov.hk/wxinfo/pastwx/metob201403c.htm, last access: 30 January 2019), in the weather review, Hong Kong had a total of 15 days of rainy weather from 25 March to 25 April 2014, as shown in Table 1.
In this paper, GAMIT (v10.50; Herring et al., 2010) software was used for processing the GPS observations based on the doubledifferenced model at a sampling interval of 30 s, and the global mapping function was adopted. The zenith total delay (ZTD) and wet horizontal gradient intervals were estimated at 0.5 and 2 h, respectively. Based on the surface pressure obtained from observed meteorological parameters, the ZHD could be obtained by the Saastamoinen model, and ZWD was isolated from ZHD. Via GMF projection, the SWD could be obtained by transforming the observed SWV.
3.2 Experimental introduction and comparison
The RMSE and bias of the improved tomography model residuals were calculated by subtracting the ECMWF water vapor density from the water vapor density of the improved pixelbased water vapor tomography model (hereinafter referred to as the improved model). In a similar way, the RMSE and bias of the traditional tomography model residuals can also be obtained from the difference between the ECMWF water vapor density and the water vapor density obtained by the traditional pixelbased water vapor tomography model (hereinafter referred to as the traditional model).
In order to evaluate the improved model, this paper investigates six scenarios, comprising the spatial distribution scenario, the everyday distribution scenario, the rainy scenario, and the nonrainy scenario. Moreover the scenarios of residuals of the water vapor density in voxels with and without penetrating GNSS rays are inspected. The definitions of six scenarios mentioned above are as follows.
The spatial distribution scenario is investigated by obtaining the RMSE and bias of the residuals from all ECMWF comparative points at all time intervals.
The everyday distribution scenario is found by obtaining the RMSE and bias of the residuals from all ECMWF comparative points in two epochs each day. In addition the overall accuracy of 32 days between 25 March and 25 April 2014 was calculated.
The rainy scenario is based on 15 rainy days between 25 March and 25 April 2014, referred to in Table 1. The RMSE and bias of the residuals are obtained from all ECMWF comparative points for all the epochs during rainy days. Similarly, the nonrainy scenario is found with the accuracy analysis of the nonrainy days.
The scenario of residuals of the water vapor density in voxels without GNSSray penetration is found by obtaining the RMSE and bias of the residuals from ECMWF comparative points without rays passing through for all the epochs each day. Conversely, the scenario with GNSSray penetration is found by obtaining the RMSE and bias of the residuals from ECMWF comparative points with ray penetration in all the epochs each day.
According to the above classifications, the accuracy of the improved model residuals and the traditional model residuals were calculated, and the accuracy of the improved model was compared with the traditional one to find out which one is better. Furthermore, the water vapor comparison with radiosonde data was designed to show if the improved model would be more efficient than the traditional one.
4.1 Accuracy information of the spatial distribution scenario
To verify whether the accuracy of the improved model is better than that of the traditional model, the layered RMSE and bias of the residuals are obtained from the tomography results (using both the optimal polynomial function of each layer and the traditional way) and the ECMWF data in all ECMWF comparative points (shown in Table 2). The calculation of the RMSE improvement percentage involved in the following tables is shown in Eq. (16):
where RMSE_{impr} is the RMSE value of the residuals calculated from the improved model, and RMSE_{trad} is the RMSE value of the residuals obtained from the traditional model.
Table 2 shows that RMSE and bias values obtained by the improved model are smaller than those of the traditional model, and the RMSE improvement percentage is positive, which indicates that the improved model has a higher accuracy than the traditional model in general. The reason of the appreciable RMSE improvement percentage in the upper region is that the value of the water vapor density in high altitudes is very small (see Fig. 2 for details), and even small water vapor changes could cause a large percentage fluctuation. In addition, the bias and RMSE in the bottom from Table 2 are not as good as those of the other higher layers, regardless of which model is used. These results could be mainly ascribed to a certain system deviation between the comparison data of ECMWF and the GNSS tomographic data. Besides, due to less voxels with GNSSray penetration in the lower layers, the observations are too insufficient to obtain good accuracy. Figure 2 also shows that the water vapor content in the bottom region is so abundant and changeable that tomography results from models could not reflect it accurately. The above reasons lead to large bias and RMSE values in the bottom tropospheric area.
4.2 The accuracy information of the everyday distribution scenario
To prove whether the accuracy of the improved model is better than that of the traditional model on the everyday time scale, the RMSE improvement percentage is obtained from all ECMWF comparative points (a total of 12) at two epochs each day using both the layered optimal polynomial functions and the traditional method. Figure 3 shows that the percentage of RMSE improvement per day is mostly positive, and the percentage of 11 April can even approach 12 %, indicating that the improvement seems to be appreciable. This improvement shows that the accuracy of the improved model is mostly superior to that of the traditional model in everyday distribution; however, on 7, 9, and 15 April, the RMSE improvement percentage is negative. This might be due to the heavy showers bringing rapid water vapor changes from 7 to 8 April and on 14 April, which is difficult for fitting the polynomial function well with the unstable water vapor. However, since negative percentages do not exceed −1 %, the accuracy of the improved model could be considered basically equivalent to that of the traditional model for these 4 days.
In addition, the overall RMSE and bias of the residuals are obtained from the ECMWF comparative points (a total of 12) in two epochs under the entire everyday distribution scenario. The statistical results are shown in Table 3 below.
Table 3 shows that the RMSE obtained by the improved model is 3.44 % smaller compared to that of the traditional one. The bias of the improved model is more close to zero, indicating that the improved model has better stability and less systematic deviation from the comparative data. The better accuracy shows the superiority of the improved model.
4.3 The accuracy information of rainy and nonrainy scenarios
To further analyze the reliability of the improved model compared with the traditional model in different weather conditions, according to the distribution of rainy days in Table 1, all the rainy days data and nonrainy days data are used separately for tomography to obtain the RMSE and bias of the residuals under corresponding weather conditions. The number of matching points is still 12 (see Fig. 1). The overall statistical results are shown in Table 4.
Table 4a shows that the RMSE and bias of the residuals calculated by the improved model are better than those of the traditional model using data from rainy days. The RMSE of the improved model is 3.68 % better than that of the traditional model, indicating that the accuracy of the new model is superior. The improved model bias is more close to zero than that of the traditional one, which means the improved model has an increase in stability and a reduction in the system error. Using data from nonrainy days, the RMSE and bias of the residuals calculated by the improved model are also better than those of the traditional model (see Table 4b). The RMSE improvement percentage is 3.21 %, also indicating the improved model has enhanced accuracy. Besides, the improved model bias is more close to zero, making the system error weak. Table 4 also shows that the RMSE improvement percentage of the rainy days is better than that of the nonrainy days. This finding shows that the improved model is more suitable for obtaining the tomographic results when heavy water vapor changes occur.
4.4 The accuracy information of voxels with and without GNSSray penetration scenarios
In the traditional pixelbased water vapor tomography model, the water vapor density in the voxels without GNSS rays passing through depends on the accuracy of the water vapor density in the adjacent voxels with GNSSray penetration. However, the improved model uses the layered optimal polynomial functions for overall fitting to obtain the water vapor density in voxels without penetrating GNSS rays. To determine if the layered optimal polynomial function of the improved model contributes better to the accuracy of the water vapor density, the scenarios of voxels with and without GNSSray penetration, as described in Sect. 3.2, were designed. After obtaining the RMSE and bias of the residuals using the improved and traditional tomography models separately under designed scenarios, the overall accuracy information of voxels with and without GNSSray penetration is shown in Table 5.
Table 5a shows that the RMSE and bias of the residuals calculated by the improved model are better than those of the traditional model in the scenario of voxels without GNSSray penetration. Moreover the RMSE of improved tomography model is 5.88 % better than that of the traditional model, and the bias decreased from 1.59 to 1.51 g m^{−3}. To a certain extent, results show that the improved model is more advantageous for obtaining the water vapor density from the voxels without GNSSray penetration, which is consistent with the initial hypothesis: the traditional model uses empirical constraint equations in Sect. 2.1.2, Eq. (7), which are unable to represent the actual distribution of the water vapor density from voxels without GNSSray penetration well. However, the improved model uses the relatively accurate water vapor density from voxels with GNSSray penetration as the observation values to further fit the water vapor density in voxels without GNSSray penetration. Therefore, the improved model can better reflect the actual layered situation of continuous water vapor changes, and the accuracy is naturally better. In addition, in the scenario of voxels with GNSSray penetration, the RMSE and bias obtained by the improved model are also superior to those of the traditional models (see Table 5b). The RMSE calculated by the improved model is 1 % higher than that of the traditional model, and the bias reduced from 1.7 to 1.65 g m^{−3}, which could prove the reliability of the improved model.
In order to doublecheck if the improved model in the scenario of voxels without GNSSray penetration shows a better result in the vertical water vapor distribution, the water vapor density profiles for different altitudes at specific times are given in Fig. 4. Two times (00:00 UTC on 11 April 2014 and 12:00 UTC on 11 April 2014) are chosen, since they correspond to the maximum percentage of RMSE improvement during the experiment period of 32 days. Figure 4 shows that in the scenario of voxels without GNSSray penetration, the water vapor profile of the improved model matches that of ECMWF data better than the traditional model, especially in the bottom layers, which again implies that the water vapor density derived from the improved model is superior to that of the traditional one in the scenario of voxels without GNSSray penetration.
Furthermore, to directly compare the vertical accuracy of the water vapor density derived from different altitudes in the scenario of voxels without penetrating GNSS rays, the tomographic results (25 March to 25 April 2014) from two different tomography models are analyzed. Figure 5 shows the percentage of RMSE improvement and the relative error of the water vapor density changing with altitude. The percentage of RMSE improvement in Fig. 5 is defined as the same as Eq. (16), and the relative error is defined by using Eq. (17):
where RE is the relative error, ρ represents the water vapor density derived from the traditional or improved tomography model, and ρ_{ECMWF} is the water vapor density derived from ECMWF grid data.
It can be observed in Fig. 5 that in the scenario of voxels without GNSSray penetration the percentage of RMSE improvement is positive in lower layers while negative in some middle and upper layers, which could prove that the improved model enhances the accuracy of tomography results in most layers when there are few voxels with GNSSray penetration, especially in the bottom layers. Due to the lack of GNSS observation data, the accuracy in the bottom layer of tomographic results is generally low. In addition, Fig. 5 shows in the scenario of voxels without GNSSray penetration, and the relative error begins to decrease with the altitude and then increases above 3 km. When the altitude is higher, the relative error becomes larger because of the small water vapor values in the upper layers, and a very tiny difference could cause a large relative error between the models and the ECMWF data.
As radiosonde data can provide fairly accurate vertical profiles of tropospheric water vapor (Niell et al., 2001), in this paper, the water vapor profiles derived from radiosonde data, as a reference, are used to validate the tomographic results from two models to show if the improved model would be more efficient than the traditional one. In Hong Kong, there is one radiosonde station located at King's Park (shown in Fig. 1) where radiosonde balloons are launched twice daily, at 00:00 and 12:00 UTC, respectively. Water vapor profiles derived from the improved model and the traditional model for the location of the radiosonde station are compared with that derived from radiosonde data at 00:00 and 12:00 UTC daily for the experimental period of 32 days. The overall statistical results are shown in Table 6. The RMSE and the bias of the improved model are 1.03 and −0.06 g m^{−3}, respectively, and the values of the traditional model are 0.82 and −0.17 g m^{−3}, respectively, which indicates that the RMSE of the improved model is not as good as the traditional model, while the bias of the improved model is a little better than that of the traditional one. The reason for poor accuracy of the improved model could be systematic differences between the training source ECMWF data and the radiosonde data. Besides, as shown in Fig. 1, the location of the radiosonde station is close to one GNSS station (HKSC), leading to the voxels for the location of the radiosonde station having GNSSray penetration. Since the improved model has the advantage of only obtaining water vapor density from voxels without GNSSray penetration, this situation cannot show the superiority of the improved model.
In addition, water vapor profiles obtained by two models and radiosonde data are compared for the specific two epochs, at 00:00 UTC on 25 March and 00:00 UTC on 7 April 2014, shown in Fig. 6. Those two times are selected because they correspond to the nonrainy day and heavy rainfall day, which could be more comprehensive and representative for the comparison results of water vapor profiles. It can be seen from Fig. 6 that two models in the nonrainy day match the radiosonde data a little better than those in the rainy day. The traditional model shows better comparison results in upper layers than those of the improved model, while the two models have almost the same comparison results in the middle and lower layers. The reason for poor performance in the lower layers might due to abundant water vapor in the bottom troposphere as well as the division of the vertical resolution. Compared to the radiosonde data, with almost the same accuracy and profile matching results as the traditional model, the improved model still has the advantage of the convenient and efficient expression.
In order to further verify the reliability of the improved model, 5 days are randomly selected from different weather conditions to make assessments on the reconstructed SWVs of two models (see Table 7). The comparison between measured SWVs and the ones derived from tomography results of two models is performed, and the average RMSE and bias are shown in Table 8. The RMSE and bias of SWVs obtained from tomography results of two models are almost the same under different weather conditions, which indicates that the reconstructed SWVs of the improved model have similar accuracy to that of the traditional one. Since the improved model has the advantage of expressing the water vapor distribution more expediently, the similar accuracy of two models in SWV comparison shows the reliability and superiority of the improved model.
In this paper, an improved pixelbased water vapor tomography model has been proposed that is much more concise and convenient in expression than the traditional one. Using only the layered optimal polynomial coefficients, the threedimensional water vapor distribution in the tomography region could be described. By using the SatRef GNSS network observation data in Hong Kong between 25 March and 25 April 2014, the RMSE and bias have been assessed in six scenarios. The scenarios include the spatial distribution scenario and the everyday distribution scenario, the rainy scenario and the nonrainy scenario, and the voxels with and without GNSSray penetration scenarios. The results demonstrate that in either case, the RMSE and bias of the improved model are better than those of the traditional model. Among these scenarios, when there are voxels without GNSSray penetration, the RMSE improvement percentage can be significantly increased up to 5.88 %, which shows that the improved model is more advantageous for obtaining the water vapor density from voxels without GNSSray penetration. Using radiosonde data for evaluation, it is proven that, with almost similar accuracies, the improved model is more efficient in expression than the traditional one. However, some shortcomings still remain in the improved model. For example, more or less, the water vapor accuracy of the improved model is still affected by the traditional model, and the layered optimal polynomial functions are limited by the size of the tomographic area and the situation of dividing voxels. In the future, the functionbased water vapor tomography model, which is free from the above limitations, should be studied. It is expected that the functionbased water vapor tomography model will be more conveniently used when the expression parameters of the function part could be obtained directly from SWVs.
All data can be accessed from the corresponding author upon request.
YY and QZ participated in the design of this study. LX carried out the literature search, data acquisition, data analysis, and paper preparation. LX and QZ carried out the study and collected important background information. All authors read and approved the final paper.
The authors declare that they have no conflict of interest.
This article is part of the special issue “Advanced Global Navigation Satellite Systems tropospheric products for monitoring severe weather events and climate (GNSS4SWEC) (AMT/ACP/ANGEO interjournal SI)”. It is not associated with a conference.
The authors would like to thank the ECMWF for providing access to the layered
meteorological data. The Lands Department of HKSAR is also acknowledged for
providing GPS data from the Hong Kong Satellite Positioning Reference
Station Network (SatRef) and corresponding meteorological data.
Edited by: Petr Pisoft
Reviewed by: three anonymous referees
Aghajany, S. H. and Amerian, Y.: Three dimensional ray tracing technique for tropospheric water vapor tomography using GPS measurements, J. Atmos. Sol.Terr. Phys., 164, 81–88, 2017.
Alber, C., Ware, R., Rocken, C., and Braun, J. J.: Obtaining single path phase delays from GPS double differences, Geophys. Res. Lett., 27, 2661–2664, 2000.
Baltink, H. K., Marel, H. V. D., and Der Hoeven, A. V.: Integrated atmospheric water vapor estimates from a regional GPS network, J. Geophys. Res.Atmos., 107, ACL 31–ACL 38, 2002.
Bender, M., Stosius, R., Zus, F., Dick, G., Wickert, J., and Raabe, A.: GNSS water vapour tomography – Expected improvements by combining GPS, GLONASS and Galileo observations, Adv. Space Res., 47, 886–897, 2011.
Bevis, M., Businger, S., Chiswell, S. R., Herring, T. A., Anthes, R. A., Rocken, C., and Ware, R.: GPS meteorology: mapping zenith wet delays onto precipitable water, J. Appl. Meteorol., 33, 379–386, 1994.
Bi, Y., Mao, J., and Li, C.: Preliminary results of 4D water vapour tomography in the troposphere using GPS, Adv. Atmos. Sci., 23, 551–560, 2006.
Bock, O., Keil, C., Richard, E., Flamant, C., and Bouin, M.: Validation of precipitable water from ECMWF model analyses with GPS and radiosonde data during the MAP SOP, Q. J. Roy. Meteor. Soc., 131, 3013–3036, 2005.
Braun, J. J.: Remote sensing of atmospheric water vapor with the global positioning system, Ph.D. thesis, University of Colorado, 2004.
Braun, J. J., Rocken, C., Meertens, C., and Ware, R.: Development of a Water Vapor Tomography System Using Low Cost L1 GPS Receivers, in: Ninth ARM Science Team Meeting Proceedings, San Antonio, Texas, 22–26 March, 1–6, 1999.
Braun, J. J., Rocken, C., and Liljegren, J. C.: Comparisons of lineofsight water vapor observations using the global positioning system and a pointing microwave radiometer, J. Atmos. Ocean. Tech., 20, 606–612, 2003.
Cao, Y.: GPS Tomographying ThreeDimensional Atmospheric Water Vapor and Its Meteorological Applications, Ph.D. Thesis, The Chinese Academy of Sciences, Beijing, China, 2012.
Chen, B. and Liu, Z.: Voxeloptimized regional water vapor tomography and comparison with radiosonde and numerical weather model, J. Geodesy, 88, 691–703, 2014.
Ding, N., Zhang, S. B., Wu, S.Q., Wang, X. M., and Zhang, K. F.: Adaptive node parameterization for dynamic determination of boundaries and nodes of GNSS tomographic models, J. Geophys. Res.Atmos., 123, 1990–2003, 2018.
Dong, Z. and Jin, S.: 3D Water Vapor Tomography in Wuhan from GPS, BDS and GLONASS Observations, Remote Sens., 10, 62, https://doi.org/10.3390/rs10010062, 2018.
Emardson, T. R., Elgered, G., and Johansson, J. M.: Three months of continuous monitoring of atmospheric water vapor with a network of Global Positioning System receivers, J. Geophys. Res., 103, 1807–1820, 1998.
Flores, A., Ruffini, G., and Rius, A.: 4D tropospheric tomography using GPS slant wet delays, Ann. Geophys. Ger., 18, 223–234, 2000.
Guo, J., Yang, F., Shi, J., and Xu, C.: An Optimal Weighting Method of Global Positioning System (GPS) Troposphere Tomography, IEEE J.STARS., 9, 5880–5887, 2016.
Herring, T. A., King, R. W., and McClusky, S. C.: Documentation of the GAMIT GPS Analysis Software release 10.4. Department of Earth and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, Massachusetts, 2010.
Hirahara K.: Local GPS tropospheric tomography, Earth Planet. Space, 52, 935–939, 2000.
Liu, Y., Chen, Y., and Liu, J.: Determination of weighed mean tropospheric temperature using ground meteorological measurements, Geospatial Inf. Sci., 4, 14–18, 2001.
Macdonald, A. E., Xie, Y., and Ware, R. H.: Diagnosis of ThreeDimensional Water Vapor Using a GPS Network, Mon. Weather Rev., 130, 386–397, 2002.
Marel, H. V. D.: Virtual GPS reference stations in the Netherlands, in: Proceedings of ION GPS98, Nashville, TN, 15–18 September, 49–58, 1998.
Niell, A. E., Coster, A. J., Solheim, F. S., Mendes, V. B., Toor, P. C., Langley, R. B., and Upham, C. A.: Comparison of measurements of atmospheric wet delay by radiosonde, water vapour radiometer, GPS, and VLBI, J. Atmos. Ocean. Tech., 18, 830–850, 2001.
Nilsson, T. and Gradinarsky, L.: Water vapor tomography using GPS phase observations: simulation results, IEEE T. Geosci. Remote, 44, 2927–2941, 2006.
Perler, D., Geiger, A., and Hurter, F.: 4D GPS water vapor tomography: new parameterized approaches, J. Geodesy, 85, 539–550, 2011.
Rocken, C., Van Hove, T., and Ware, R.: Near realGPS sensing of atmospheric water vapor, Geophys. Res. Lett., 24, 3221–3224, 1997.
Rohm, W. and Bosy, J.: Local tomography troposphere model over mountains area, Atmos. Res., 93, 777–783, 2009.
Rohm, W., Zhang, K., and Bosy, J.: Limited constraint, robust Kalman filtering for GNSS troposphere tomography, Atmos. Meas. Tech., 7, 1475–1486, https://doi.org/10.5194/amt714752014, 2014.
Saastamoinen, J.: Atmospheric Correction for the Troposphere and Stratosphere in Radio Ranging Satellites, The Use of Artificial Satellites for Geodesy, Geophysical Monograph Series, edited by: Henriksen S. W., Mancini A., and Chovitz B. H., America, Volume 15, https://doi.org/10.1029/GM015p0247, 1972.
Seko, H., Shimada, S., Nakamura, H., and Kato, T.: Threedimensional distribution of water vapor estimated from tropospheric delay of GPS data in a mesoscale precipitation system of the Baiu front, Earth Planet. Space, 52, 927–933, 2000.
Troller, M., Burki, B., Cocard, M., Geiger, A., and Kahle, H. G.: 3D refractivity field from GPS double difference tomography, Geophys. Res. Lett., 29, 21–24, 2002.
Vollach, U., Buecherl, A., Landau, H., Pagels, C., and Wagner, B.: MultiBase RTK Positioning Using Virtual Reference Stations, in: Proceedings of ION GPS2000, Salt Lake City, 19–22 September, 123–131, 2000.
Weckwerth, T. M., Parsons, D. B., Koch, S. E., Moore, J. A., Lemone, M. A., Demoz, B., Flamant, C., Geerts, B., Wang, J., and Feltz, W. F.: An overview of the international H_{2}O project (IHOP_2002) and some preliminary highlights, B. Am. Meteorol Soc., 85, 253–277, 2004.
Yao, Y. B., Zhao, Q. Z., and Zhang, B.: A method to improve the utilization of GNSS observation for water vapor tomography, Ann. Geophys., 34, 143–152, https://doi.org/10.5194/angeo341432016, 2016.
Zhang, B., Fan, Q., Yao, Y., Xu, C., and Li, X,: An Improved Tomography Approach Based on Adaptive Smoothing and Ground Meteorological Observations, Remote Sens., 9, 886, https://doi.org/10.3390/rs9090886, 2017.
Zhao, Q. and Yao, Y.: An improved troposphere tomographic approach considering the signals coming from the side face of the tomographic area, Ann. Geophys., 35, 87–95, https://doi.org/10.5194/angeo35872017, 2017.
Zhao, Q., Yao, Y., Yao, W., and Xia, P.: An optimal tropospheric tomography approach with the support of an auxiliary area, Ann. Geophys., 36, 1037–1046, https://doi.org/10.5194/angeo3610372018, 2018.
 Abstract
 Introduction
 An improved pixelbased water vapor tomography model
 Experiment
 Interpretation of six scenario results
 Water vapor comparison with radiosonde data
 SWV comparison
 Conclusions
 Data availability
 Author contributions
 Competing interests
 Special issue statement
 Acknowledgements
 References
 Abstract
 Introduction
 An improved pixelbased water vapor tomography model
 Experiment
 Interpretation of six scenario results
 Water vapor comparison with radiosonde data
 SWV comparison
 Conclusions
 Data availability
 Author contributions
 Competing interests
 Special issue statement
 Acknowledgements
 References