Optimization of GPS water vapor tomography technique with radiosonde and COSMIC historical data

The near-real-time high spatial resolution of atmospheric water vapor distribution is vital in numerical weather prediction. GPS tomography technique has been proved effectively for three-dimensional water vapor reconstruction. In this study, the tomography processing is optimized in a few aspects by the aid of radiosonde and COSMIC historical data. Firstly, regional tropospheric zenith hydrostatic delay (ZHD) models are improved and thus the zenith wet delay (ZWD) can be obtained at a higher accuracy. Secondly, the regional conversion factor of converting the ZWD to the precipitable water vapor (PWV) is refined. Next, we develop a new method for dividing the tomography grid with an uneven voxel height and a varied water vapor layer top. Finally, we propose a Gaussian exponential vertical interpolation method which can better reflect the vertical variation characteristic of water vapor. GPS datasets collected in Hong Kong in February 2014 are employed to evaluate the optimized tomographic method by contrast with the conventional method. The radiosonde-derived and COSMICderived water vapor densities are utilized as references to evaluate the tomographic results. Using radiosonde products as references, the test results obtained from our optimized method indicate that the water vapor density accuracy is improved by 15 and 12 % compared to those derived from the conventional method below the height of 3.75 km and above the height of 3.75 km, respectively. Using the COSMIC products as references, the results indicate that the water vapor density accuracy is improved by 15 and 19 % below 3.75 km and above 3.75 km, respectively.


Introduction
The temporal and spatial variation of water vapor is quite important in numerical weather forecasting especially for several hours of small-scale disastrous-weather monitoring and forecasting (Hamill and Church, 2000;Hart and Forbes, 1999;Posada et al., 2012).Accurate and reliable weather forecasting requires high-accuracy water vapor distribution information in both horizontal and vertical directions.Traditionally, the radiosonde sounding, microwave radiometer, meteorological satellite and laser radar technology are employed to obtain water vapor vertical distribution (Brettle and Galvin, 2003).These observation means have some drawbacks such as low temporal or spatial resolutions, high cost and weather dependence.In recent years, GPS tomography technique has been demonstrated as an effective means to acquire the three-dimensional (3-D) distribution of water vapor and can compensate for these disadvantages (Troller et al., 2006;Bender and Raabe, 2007;Bender et al., 2011;Champollion et al., 2009;Lutz et al., 2010;Perler et al., 2011;Rohm, 2013;Jiang et al., 2014;Benevides et al., 2016).
In the GPS tomography technology, the input observations are the slant water vapor (SWV).To obtain the SWV, the zenith troposphere delay (ZTD) is firstly determined using GPS double-difference methods.The ZTD is composed of two parts, namely zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD) (Davis et al., 1985).The ZHD can be obtained by empirical models and then the ZWD is obtained by removing the ZHD from the ZTD.Afterwards, the ZWD can be projected to the slant wet delay (SWD) along the line of sight using a Niell mapping function (Niell, 1996).Eventually, the SWD is converted to the SWV with a humidity transforming factor.In order to reconstruct the 3-D distribution of the water vapor, the tomography area is discretized Published by Copernicus Publications on behalf of the European Geosciences Union.
into voxels in both horizontal and vertical directions.Then, the SWV can be expressed by an integral of the water vapor density along the ray path from the receiver to the top boundary of the grid (Flores et al., 2000;Troller et al., 2006;Bender and Raabe, 2007;Bender et al., 2011;Perler et al., 2011).
In the tomography processing, the wet tropospheric delay, humidity transforming factor, grid division and tomography model are four crucial aspects affecting the tomography results.To obtain the ZWD, the ZHD is essential to remove from the ZTD.The ZHD is usually computed using the classical models such as Saastamoinen model, Hopfield model and Black model (Saastamoinen, 1972;Hopfield, 1971;Black, 1978).However, the classical models are inaccurate because they require measurements of pressure with a precision of 0.1 hPa, which can be obtained from local meteorological stations coupled to GPS stations (Tregoning and Herring, 2006).However, these measurements are usually obtained from climatological models such as ERA (Each Re-Analysis)-Interim and NCEP (National Center for Environmental Prediction) with lower precision.As a result, the accuracy of estimated ZWD will be degraded.
The humidity conversion factor is dependent on the atmospheric weighted mean temperature (Mateus et al., 2014), which can usually be estimated utilizing the surface temperature measurements.In the actual calculation of the SWV, the humidity conversion factor is usually determined using empirical models or even simply regarded as a constant (Shi and Gao, 2009;Bosy et al., 2012;Shangguan et al., 2013;Adavi and Mashhadi-Hossainali, 2014;Chen and Liu, 2014).However, the conversion factor is varied in different areas and seasons (Jiang et al., 2014).
In the tomography model, the grid is usually divided in an even vertical height with a constant top height of 10 km and the water vapor density is considered to be same inside each voxel (Flores et al., 2000;Hirahara, 2000;Troller et al., 2006;Nilsson and Gradinarsky, 2006;Bastin et al., 2007).However, such a grid model does not take into account the significant difference of actual water vapor distribution inside or outside each voxel for the vertical direction.
In this study, an optimized method for near real-time GPS 3-D troposphere tomography is developed using the external radiosonde and COSMIC historical data.Specifically, radiosonde and COSMIC products are utilized to improve the estimation accuracy of tropospheric ZHD as well as humidity conversion factor.Furthermore, the water vapor layer top (WVLT) and the dense layer of vapor top (DLVT) are dynamically determined for the purpose of grid division with an uneven voxel height and an unfixed grid top.In addition, a Gaussian exponential interpolation method is proposed to consider the water vapor variations within each vertical voxel.The ground-based GPS observation data from the Hong Kong SatRef network in February 2014 are used to verify the feasibility and superiority of the optimized tomography method.
The rest of the paper is organized as follows.Section 2 presents the improved methods in the aspects of ZHD estimation, humidity conversion factor, grid division and interpolation function.Section 3 analyzes the result improvement with the proposed tomography methods.Section 4 presents the validation of the tomographic results.The conclusions are included in Sect. 5.

Tomography model
The traditional troposphere tomography model based on discrete voxels was first proposed by Flores et al. (2000).It can be described by the formula below: where SWV is the slant water vapor; ρ w is the density of liquid water; s is the slant path that GNSS signals pass through troposphere; ρ(s) is the water vapor density.After being discretized by the fourth-order Newton-Cotes formula, the equation can be expressed as (Schwarz, 1997;Stoer and Bulirsch, 1980) where C (4) ij is the fourth-order coefficients of the j th segment point within the ith grid; ρ(s i j ) is the corresponding water vapor density values; m denotes the number of the grids which the signals pass through; S represents the distance that the GPS signals span the voxels; SWV is the noise.
ZWD is obtained using GPS double-difference method and then is projected to the SWD through the Niell mapping function (Niell, 1996).Further, the SWD can be transformed to the SWV using a conversion factor (Song, 2004): where is the humidity conversion factor, which is geographically dependent and varied from location to location.
Due to the nearly cone geometry of the GPS observations, the GPS signals cannot pass through all voxels (Bender and Raabe, 2007;Benevides et al., 2016).As a result, the tomographic system cannot be inverted due to too many zeros in the design matrix.Therefore, proper constraints have to be imposed to solve this issue.In horizontal direction, the Gaussian weighted method (Song, 2004) can be used for horizontal constraint.Given that the vertical water vapor is usually decreased with the increase of the height, the Gaussian exponential model is used to model the vertical distribution as where ρ(h) denotes the water vapor density at the height of h; h z represents the height index of water vapor; ρ C is a constant value of water vapor density; h 0 is a constant.ρ C and h 0 may be determined using radiosonde or COSMIC products.Based on Eq. (4), Eq. ( 5) is used as the vertical constraint to establish the relationship between water vapor density in adjacent vertical layers: where the subscripts "i", "j " and "k" denote the index of voxel in the east-west, north-south and vertical directions, respectively; h k is the height of the kth voxel.
The a priori water vapor information is necessary to be used as the initial values of tomography solutions.In order to obtain more accurate a priori information of water vapor density, it will be helpful to improve the accuracy of atmospheric elements by assimilating surface meteorological observations into tomographic system (Perler et al., 2011;Jiang et al., 2014).Since the synoptic observation products contain the pressure, temperature and relative humidity of the station, we interpolate the temperature and relative humidity into all of the grids.In the beginning, we interpolate horizontal relative humidity using Gaussian weighted constraint method and vertical relative humidity using Eq. ( 5).It is similar to the horizontal temperature interpolation, but the vertical temperature is interpolated by computing the temperature decay rate using radiosonde products.After the temperature and relative humidity are obtained, the water vapor density of all the voxels can be computed using the temperature and relative humidity and utilized as the prior values in the tomography (Hardy, 1998;Song, 2004).

ZHD computation
As mentioned in Sect. 1, the ZHD is required in order to achieve the ZWD which is associated with the tomography input observations of SWV.The computation of ZHD is usually made using empirical models with input of atmospheric pressure and temperature.However, the empirical models are not accurate enough in some areas where there exists significant height difference (Liu et al., 2000).In these cases, the radiosonde and radio occultation data products can be utilized to compensate the shortcomings of the empirical models.Equation ( 6) shows how to calculate the ZHD using the radiosonde products and GPS radio occultation products (Liu et al., 2000): where I denotes the height of the atmosphere whose pressure is less than 0.02 hPa; I 0 denotes the initial height of computing the ZHD; h is the geopotential height; N d is the dry air refractivity; t and T are the temperature in degrees Celsius and kelvin, respectively; P d denotes the dry pressure.
Via the stepwise regression analysis, we find that the offsets of Saastamoinen-derived and radiosonde-derived ZHD mainly correlate with the temperature and pressure.In order to improve the accuracy of the ZHD, we compensate the Saastamoinen-derived ZHD with a model at a similar form to the Black formula (Black, 1978).
where ZHD is the ZHD corrections; P is pressure; c and d are two parameters of the calibration model to be determined by radiosonde and COSMIC profile historical data.

Humidity conversion factor
The humidity conversion factor varies in different areas and seasons, and it can be expressed as a function of the tropospheric weighted mean temperature T m (Bevis et al., 1994): where ρ w is the density of liquid water; k 1 , k 2 and k 3 are constants as k 1 = 77.6K hPa −1 , k 2 = 70.4K hPa −1 and k 3 = 3.739 × 10 5 K hPa −1 (Bevis et al., 1994); m d and m w are the molar masses of dry atmosphere and water vapor, respectively; R indicates universal gas constant; P w means water vapor pressure whose unit is hPa.
The parameter conversion factor mainly depends on T m .Traditionally, T m is determined by a single means such as ECMWF (European Centre for Medium-Range Weather Forecasts), NCEP, radiosonde sounding and so on.Compared to the radiosonde technique, the accuracies of the ECMWF and NCEP are relatively lower.However, the temporal resolution of the radiosonde is low and the distribution of radiosonde stations is sparse.In this case, we jointly use the radiosonde and COSMIC products to estimate the T m for the purpose of improving its temporal resolution.The conversion factor is then determined using the T m .

Tomography grid division
In general, the lower limit and upper limit of the tomography grid is the height from the ground to tropopause.But in S. Ye et al.: GPS water vapor tomography technique fact the water vapor is mainly concentrated at a height significantly lower than the tropopause.If the tropopause is chosen as the grid top, the solutions of the tomography inversion are possibly negative as the water vapor is very sparse near the height of the tropopause (Flores et al., 2000).In this study, we define a varied water vapor layer top (WVLT) as the upper limit of the tomography grid based on the precipitable water vapor (PWV) variations acquired from the radiosonde and COSMIC data.Thus, the height of the grid top is decreased, conversely increasing the effective number of the satellite rays.In the tomography, only the rays that penetrate into the grid from the top boundary are used for the tomography processing.
The purpose of GPS tomography is to reconstruct the vertical distribution of water vapor density.The vertical grid division is vital to affect the tomography solutions.Traditionally, there are two ways to divide the grid.One is uniform division (Flores et al., 2000;Xia et al., 2013) and the other is non-uniform division (Perler et al., 2011;Rohm, 2012Rohm, , 2013;;Chen andLiu, 2014, andJiang et al., 2014).Considering the actual distribution characteristics of water vapor density that is sparse in high layers and dense in low layers, we use nonuniform division in this study.

Interpolation function inside each voxel
Perler et al. ( 2011) proposed a numerical integration model parameterization (NIMP) tomography method.The basic idea is that the water vapor within each voxel is considered to be unevenly distributed and can be calculated by the horizontal and vertical interpolation methods.
In this study, the difference of water vapor density within each voxel is also taken into account.The horizontal and vertical interpolations are made to obtain the water vapor density values of interior voxels.In the horizontal direction, the Barnes interpolation algorithm (Jiang et al., 2014) is normally used.In the vertical direction, there are many interpolation methods such as linear interpolation and cubic spline interpolation (Perler et al., 2011).The linear interpolation algorithm is simple but its accuracy is relatively poor.The cubic spline interpolation is complex in implementation.Besides, these traditional interpolation means belong to pure mathematics methods, without taking into account the inherent characteristics of the water vapor density variations in the vertical direction.In this study, a Gaussian exponential-based interpolation algorithm derived from Eq. ( 5) is proposed for the vertical interpolation, as shown in Eq. ( 11).
where ρ i k (h) denotes the water vapor density in the "kth" voxel; i represents the "ith" GPS ray; h k and h k+1 denote the corresponding height of the bottom and top of the kth voxel, respectively; ρ i k (h k ) and ρ i k+1 (h k ) represent the water vapor density on the bottom and top of the "kth" voxel, respectively.
3 Processing results and analysis

Data description
The COSMIC (Constellation Observation System for Meteorology Ionosphere and Climate) occultation is a joint Taiwan-US science mission for weather, climate, space weather and geodetic research (Schreiner et al, 2007;Anthes et al., 2008).The radiosonde sounding technique is an important means for meteorology study from the ground to the lower stratosphere.Both COSMIC and radiosonde observations have been the important data source for weather research and climate analysis (Kuo et al., 2005;Kishore et al., 2011).
The COSMIC Data Analysis and Archive Center (CDAAC) provides both the real-time and post-processed data products (www.cosmic.ucar.edu).The post-processed profiles that are available in 6-week latency are used for this study.As one type of the post-processed data products, "wetPrf" profile offers water vapor pressure, temperature and refractivity, which are almost evenly distributed around the globe.Radiosonde measurement has a high accuracy, high vertical resolution, long-term stability and all-weather capability for obtaining water vapor density, pressure and temperature in the troposphere and lower stratosphere.However, these measurements are only distributed in the limited areas where there are the radiosonde stations.Particularly, they are rare in the Southern Hemisphere and cannot cover the oceans.
Radiosonde data was collected at the 45004th radiosonde station equipped with a Vaisala RS92 sensor, which offers the world's highest level of performance in measuring the meteorological parameters, including pressure, temperature and humidity.The atmospheric pressure and temperature measured by the sensor have accuracies of better than 1 hPa and 0.5 • C, respectively.For the humidity sensor, the total uncertainties in sounding and repeatability are 5 and 2 %, respectively (www.vaisala.com).The COSMIC RO uses a compact, low-power and low-cost sensor, which provides the high-accuracy meteorological elements with averaged profiles temperature of less than 0.1 K. (http://www.cosmic.ucar.edu/ro.html) .
The ground-based GPS observation data are collected from the Hong Kong SatRef network, which consists of a total number of 12 continuous operational stations with inter-station distance of 10-15 km, as seen in Fig. 1.LEICA GRX1200 + GNSS receivers are equipped at all 12 stations with a data sampling rate of 5s.One-month GPS datasets on 1-28 February 2014 were collected in Hong Kong.In addition, the "wetPrf" profiles whose address is the same as Hong Kong and radiosonde products at the "45004th" station in February of 2007-2013 were used as historical data to optimize the tomography solutions.There are 81 "wetPrf" profiles during this period of time.The quality flags of all "wet-Prf" profiles are examined and the flags indicate "bad = 0", suggesting that these profiles have passed quality control successfully.
In order to simulate a near-real-time GPS tomographic experiment, it is good to use a sliding time window strategy (Foster et al., 2005).We use a 6 h interval time window and step forward 1 h at a time (Benevides et al., 2015).The ZTD estimates are obtained using GAMIT software in this study (Herring et al., 2010;Benevides et al., 2015).

ZHD calibration
The historical radiosonde data collected at the "45004th" radiosonde station in February from 2007 to 2013 are used to calibrate the ZHD model.The radiosonde product only provides the atmosphere element below the height of 30 km.Therefore, we use the ERA-Interim product (http://apps.ecmwf.int/datasets/data/interim-full-daily) to obtain atmosphere element above the height of 30km.The radiosondederived ZHD is used as references to obtain the ZHD offsets for the computation of c and d in Eq. ( 8).The established ZHD model can be expressed as where ZHD s is derived from the Saastamoinen model.In order to assess the feasibility of Eq. ( 8), the correlation between the ZHD offsets ( ZHD) and the pressure and temperature are displayed in Table 1, respectively.
As seen in Table 1, the ZHD are highly correlated to T , P and P /T , which explains why we express the ZHD as a function of the T , P and P /T in Eq. ( 8).Using the established model of Eq. ( 12), the calibration effect is assessed using datasets in February 2014.
Figure 2 shows the ZHD offsets before and after calibrations with respect to radiosonde-derived ZHD.The largest ZHD offsets using the calibrated model are only 5.2 mm, and the rms of the ZHD offsets is 3.2 mm, which is significantly smaller than the ones before calibration at 14.1 mm.From the statistical results, it is concluded that the established ZHD calibration model by using the radiosonde historical data improves the ZHD accuracy.

Humidity conversion factor determination
According to Eqs. ( 3) and ( 9), if the conversion accuracy is better than 1 mm, the precision of the weighted mean temperature T m must be less than 3 K.Since T m varies with seasons and region, we combine radiosonde product and COS-MIC "wetPrf" profile to calculate T m based on Eq. ( 10) using measurements in February from 2007 to 2013.Then, we fit the relationship between T m and the surface temperature T as follows: After obtaining T m using Eq. ( 13), the conversion factor can be determined using Eq. ( 9).In order to examine the improvement effect of conversion factor, we compare the PWV in Hong Kong in February 2014 using the different conversion factors obtained from the traditional method and our improved method.In the traditional method, we use a constant value of 0.1538 as the conversion factor (Shi and Gao, 2009).The PWV derived from the "45004th" radiosonde sta-  As an example, Table 2 displays the derived humidity conversion factors from 1 to 7 February 2014 using our improved method.The humidity conversion factors vary around 0.161, which is different from the constant value of 0.1538.
Table 3 shows the PWV comparison based on datasets in entire February 2014 in Hong Kong using different conversion factors.From the comparison results, it is found that PWV rms is significantly smaller using the improved method for humidity conversion factor computation.

Tomography grid division
It is well known that WVLT varies with regions and seasons.Since the radiosonde product and COSMIC occultation product provide the vertical distribution of PWV, we use radiosonde product and COSMIC product to estimate the changes of the PWV for determining the WVLT.Specific steps are as follows: (1) use the CPT (cold point tropopause) method to determine tropopause; (2) estimate the PWV at different heights; (3) calculate the proportion of PWV at different height over the sum at all heights; and (4) take the height as WVLT if the proportion is less than or equal to 0.001.Figure 3 shows the obtained WVLT using radiosonde and COSMIC products, respectively.
From the blue curves in Fig. 3, it can be seen that the PWV approaches 0 from the height of WVLT to tropopause.The estimated PWV using radiosonde product and COSMIC product exhibits an exponential distribution as seen in the blue curve.We define DLVT (dense layers of vapor top) as the corresponding height of the blue points whose distance is the nearest from the coordinate origin, as seen in Fig. 3. From the ground to the DLVT, the decreasing rate of the PWV becomes faster with the increase in height, whereas the decreasing rate of the PWV is slower from the DLVT to WVLT.Since the GPS network in Hong Kong area is dense, we divide grids into 10 km × 10 km squares in horizon direction, as seen in Fig. 1.We utilize the radiosonde and COSMIC historical data in February from 2010 to 2013 to determine the heights of DLVT and WVLT, and they are 3.88 and 3.61 km for the DLVT and 7.5 and 7.47 km for the WVLT, respectively.Then we take their average values, which are 3.75 and 7.5 km as the final DLVT and WVLT.Thus, we divide the vertical layers into two stages.The first stage is from the ground to the DLVT (3.75 km).We divide it into 10 layers and each layer height is 375 m.The second stage is from the DLVT (3.75 km) to the WVLT (7.5 km).We divide it into five layers and each layer's height is 750 m.

Result validations
The same datasets as described in Sect.3.1 are used for the result validations.The radiosonde balloon is carried aloft once every 12 h at the 45004th station, and the configured sensors on the radiosonde measure profiles of temperature, pressure, relative humidity and so on.In addition, the "wet-Prf" profiles provide the water vapor pressure, temperature, pressure and so on.Thus, the "wetPrf" profiles collocated in Hong Kong and the 45004th radiosonde products in February 2014 are used to evaluate the accuracy of the tomography results.We utilize GAMIT software to obtain the ZTD based on GPS data from the Hong Kong SatRef network with IGS (International GNSS Service) ultra-rapid product orbit file.The Saastamoinen model is calibrated using observations from the "45004th" radiosonde station and the humidity conversion factor is determined using the COSMIC products and radiosonde products in Hong Kong.Then, the ZHD is obtained using Eq. ( 12), and the ZWD is obtained by deducting the ZHD from the ZTD.The SWV is computed using Eq. ( 3) with optimized conversion factors.Finally, we estimate the 3-D distribution of atmospheric water vapor using parameterized approaches in which Eq. ( 11) is used for interpolation in the vertical direction and Eq. ( 5) is applied as the vertical constraints.Kalman filtering algorithm is used for tomography solutions.In detail, the unknowns are the water vapor densities at the grid points.The state transition matrix is a unit matrix.The process noise matrix is acquired from the statistics of radiosonde-derived water vapor densities.Figure 4 describes the water vapor density changes at different heights.It shows that the water vapor density significantly changes below 3.75 km and decreases in the latest few hours.In order to evaluate our optimized method, the tomography results are compared with those derived from the traditional tomography method in which exponential model is used as the vertical constraint and cubic spline is used for vertical interpolation.In addition, the ZHD and the humidity conversion factor before optimization are used in the traditional method.Tomography solutions are compared with external results derived from the radiosonde and COSMIC products.As the water vapor density remains stable above DLVT (3.75 km), whereas it changes significantly below DLVT (3.75 km), the statistics of tomography-derived water vapor density are made above 3.75 km and below 3.75 km.

Comparison of tomography-derived results and radiosonde-derived results
For the sake of convenient comparison, the tomography results in grams per cubic meter are transformed into the PWV in millimeters using Eq. ( 14) (Esteban et al., 2013): Applying Eq. ( 14), the water vapor densities have been transferred to PWV. Figure 5 shows the difference of PWV obtained from the radiosonde product and tomography in the same region (latitude 22.31 • and longitude 114.16 • ).
As can be seen from Fig. 5, the tomography solutions have a good agreement with radiosonde results.Using our optimized methods, the rms of the difference between tomography-derived and radiosonde-derived PWV is 1.8 mm, whereas it is 2.3 mm for the traditional method.
Radiosonde products provide 3-D distribution of atmospheric elements, such as temperature, pressure, dew point, mixing ratio and relative humidity.We can obtain the "wet" pressure according to the pressure and mixing ratio and further use the "wet" pressure to compute water vapor density (Song, 2004).Figure 6 shows the comparisons of water vapor densities obtained from tomography and radiosonde products on 16-17 February 2014.
The changing trends of water vapor density (WVD) between tomography-derived and radiosonde-derived results have a very good agreement.Despite a significant difference between tomography and radiosonde WVD in the lower  height, the result from our optimized method is better than the results from the traditional method since the blue curve is closer to the green curve.In addition, we also provide the statistical results between tomography-derived and radiosonde-derived water vapor densities above 3.75 km and below 3.75 km, respectively, using 28-day datasets in February 2014 in Hong Kong (Table 4).Table 4 provides both the mean and the rms of the differences between tomography-derived and radiosonde-derived water vapor density.In terms of the statistical results, the accuracy of tomography based on our optimized method is significantly better than the traditional method at heights above 3.75 km as well as below 3.75 km.Compared with radiosonde products, the rms statistical results indicate that the water vapor density accuracy from our optimized method is improved by 15 and 12 % compared to traditional method below 3.75 km and above 3.75 km, respectively.
Table 5.The rms and mean of differences between tomographyderived and COSMIC-derived water vapor density above 3.75 km and blow 3.75 km, respectively (g m −3 )."RD" represents the relative differences.COSMIC "wetPrf" profile provides "wet" pressure, temperature and pressure.Water vapor density (WVD) can be estimated using "wet" pressure (Song, 2004).We selected the COSMIC radio occultation events which occurred near or inside Hong Kong.Figure 7 depicts the comparisons of water vapor densities obtained from tomography and COSMIC products on 5 and 13 February 2014.
It can be seen that the changing trends of water vapor density have a good consistency between tomography-derived and COSMIC-derived results in Fig. 7.The WVD derived from our optimized method is better than the WVD from the traditional method since the pink curve is closer to the grey curve.Total eight COSMIC occultation events occurred near or inside Hong Kong in February 2014.These results are used to compare with the results of tomography-derived water vapor above and below the height of 3.75 km, respectively, as seen in Table 5.
Table 5 provides both the mean and the rms of the differences between tomography-derived and COSMIC-derived water vapor density.In terms of the statistical results, the accuracy of tomography based on our optimized method is significantly better than the traditional method.Compared with COSMIC products, the rms statistical results confirm that the water vapor density accuracy from our optimized method is improved by 15 and 19 % compared to the traditional method below 3.75 km and above 3.75 km, respectively.

Conclusions
The tomography technique is optimized in a few aspects under the help of radiosonde and COSMIC historical data.Firstly, in order to improve the zenith wet delay (ZWD) accuracy, the regional zenith hydrostatic delay (ZHD) models are optimized by compensating the estimates from the Saastamoinen model.Secondly, the regional conversion factor of converting the ZWD to the perceptible water vapor (PWV) is refined by improving the quality of the tropospheric weighted mean temperature.Next, we develop a method for dividing the tomography grid with an uneven voxel height and a var-  "COSMIC" represents the water vapor density derived from wet-Prf product."Impr" represents the tomography-derived water vapor density using our optimized method, and "Trad" represents the tomography-derived water vapor density using the traditional method.
ied water vapor layer top.Finally, we propose a Gaussian exponential vertical interpolation method which can better reflect the vertical variation characteristic of water vapor.
GPS datasets collected in Hong Kong in February 2014 are used to assess the optimized tomographic method with a comparison to the conventional method.The radiosondederived and COSMIC-derived water vapor density values are used as references to analyze the tomographic results.Using radiosonde products as references, the test results obtained from our optimized method indicate that the water vapor density rms is 2.52 and 0.86 g m −3 below and above 3.75 km, respectively, which is improved by 15 and 12 % compared to those of conventional method.Using the COSMIC products as references, the results indicate that the water vapor density rms is 1.24 and 0.72 g m −3 from our optimized method below and above 3.75 km, respectively, which is improved by 15 and 19 % compared to the conventional method.

Data availability
The Survey and Mapping Office of the Lands Department (2016), Hong Kong Special Administrative Region, provides GPS data and meteorological data.These datasets are from the Hong Kong Satellite Positioning Reference Station (SatRef) and can be made freely available for public access (https://www.geodetic.gov.hk/).The 45004 Kings Park radiosonde observations (UWYO, 2016) can be downloaded from the following website: http://weather.uwyo.edu/upperair/sounding.html.In addition, the Constellation Observing System for Meteorology, Ionosphere and Climate-1 program (CDAAC, 2016) offers the COSMIC radio occultation data, which can be freely obtained from http://cosmic-io.cosmic.ucar.edu/cdaac/index.html.

Figure 1 .
Figure 1.The GPS station distribution and horizontal grid division in Hong Kong.

Figure 3 .
Figure 3.The top of the water vapor layer and the top of the water vapor dense layer obtained by radiosonde and COSMIC products.
Figure 4 shows the 3-D distribution of atmospheric water vapor density with a WVLT height of 7.5 km on 19 February 2014.

Figure 5 .
Figure 5. Comparisons of PWV obtained from tomography and radiosonde product in February 2014."RadP-ImpP" represents the difference of radiosonde-derived and optimized tomography-derived PWV."RadP-TradP" represents the difference of radiosonde-derived and traditional tomography-derived PWV.

Figure 6 .
Figure6.Comparisons of water vapor densities (WVDs) obtained from tomography and radiosonde on 16-17 February 2014."Rad" represents the radiosonde-derived water vapor density."Impr" represents the tomography-derived water vapor density using our optimized method, and "Trad" represents the tomography-derived water vapor density using the traditional method.

Figure 7 .
Figure 7. Comparisons of water vapor densities obtained from tomography and COSMIC products on 5 and 13 February 2014."COSMIC"represents the water vapor density derived from wet-Prf product."Impr" represents the tomography-derived water vapor density using our optimized method, and "Trad" represents the tomography-derived water vapor density using the traditional method.

Table 1 .
Correlation coefficients between the ZHD offset and temperature and pressure.
Figure 2. ZHD offsets before and after calibrations with respect to radiosonde-derived ZHD using datasets in Hong Kong in February 2014.

Table 2 .
Derived humidity conversion factor for consecutive 7 days using improved method.Date format is YYYY/MM/DD.

Table 3 .
PWV comparison in Hong Kong in February 2014 using traditional and improved conversion factors (mm).

Table 4 .
The statistical results between tomography-derived and radiosonde-derived water vapor density above 3.75 km and below 3.75 km, respectively (g m −3 )."RD" represents the relative differences.