Development of a New Vertical Water Vapor Model for GNSS Water Vapor Tomography

Abstract

One of the main challenges of Global Navigation Satellite System (GNSS) tomography is in solving ill-conditioned system equations. Vertical constraint models are typically used in the solution procedure and play an important role in the quality of the GNSS tomography, in addition to helping resolve ill-posed problems in system equations. In this study, based on a water vapor (WV) parameter, namely IRPWV, a new vertical constraint model with six sets of coefficients for six different WV states was developed and tested throughout 2019 in the Hong Kong region with four tomographic schemes, which were carried out with the model and the traditional vertical constraint model using three different types of water vapor scale height parameters. Experimental results were numerically compared against their corresponding radiosonde-derived WV values. Compared with the tests that used the traditional model, our results showed that, first, for the daily relative error of WV density (WVD) less than 30%, the new model can lead to at least 10% and 49% improvement on average at the lower layers (below 3 km, except for the ground surface) and the upper layers (about 5–10 km), respectively. Second, the skill score of the monthly root-mean-square error (RMSE) of layered WVD above 10 accounted for about 83%, 87%, and 64%. Third, for the annual biases of layered WVD, the new model significantly decreased by 1.1–1.5 g/m³ at layers 2–3 (about 1 km), where all schemes showed the maximal bias value. Finally, for the annual RMSE of layered WVD, the new model at the lower (about 0.6–3 km) and upper layers improved by 13–42% and 5–47%, respectively. Overall, the new model performed better on GNSS tomography and significantly improved the accuracy of GNSS tomographic results, compared to the traditional model.

1. Introduction

Atmospheric water vapor (WV) is one of the most important components of the troposphere, and it plays a key role in atmospheric processes and affects almost all atmospheric processes [1], e.g., it is strongly associated with climate change, atmospheric radiation, weather patterns, and the hydrologic cycle [2,3]. Information on the distribution and variation of WV is essential for meteorological applications. Traditional WV detection techniques, such as radiosondes and automatic weather stations, can obtain WV profiles with high accuracy. However, the spatial and temporal resolution of WV cannot be described by demand in the related applications [4].

Bevis et al. [5] proposed a tropospheric tomography using meteorological Global Positioning System (GPS) networks to reconstruct three-dimensional WV fields based on GPS-derived slant delays. It made up for the shortcomings of the traditional techniques, with the advantages of low-cost, high accuracy, all weather conditions, and high temporal resolution. Hence, it has become one of the most promising techniques to reconstruct the four-dimensional distribution of atmospheric WV [6]. Plenty of Global Navigation Satellite System (GNSS) tropospheric WV tomography (GNSS tomography) studies and experiments have been carried out in different regions [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. A GNSS tomographic domain covering the ground-based GNSS receiver network is usually cuboid-shaped or box-shaped within the troposphere, and it is first needed to discretize into many tomographic voxels in the process of tomographic solving. On account of the geographical distribution of the receivers being sparse and uneven and the geometric characteristics of the satellite constellation of GNSS, the number of satellite signals received by each receiver at a tomographic epoch is limited. In addition, all observed GNSS signal ray paths from a receiver to satellites present an inverted triangular cone shape [10,26,27,28]. Therefore, at a tomographic epoch, the number of intersecting GNSS signal rays in a voxel is insufficient and a large number of voxels are not even crossed by any rays (i.e., non-intercepted voxels), especially the bottom voxels where the WV is most active and abundant. The situations described above lead to an ill-conditioned GNSS tomographic model, i.e., the singularity condition, in which an ill-posed problem appears in the inversion of the tomographic equations, so the inverse of the design tomographic matrix cannot be directly obtained [6,7,9,29,30]. Several methods have been developed for solving the ill-posed problem, such as (1) adding inter-spatial constraints [6,8,12,13,31,32,33,34,35]; (2) adding prior information of WV [31,36,37,38]; (3) expanding the tomographic area to increase the number of available observations [19,37,39,40,41,42]; (4) using multi-GNSS observations [41,43,44,45,46,47]; (5) utilizing other observations [18,48,49,50].

In the preliminary stage of GNSS tomography, the inter-spatial constraint is a smoothing function, which uses a weighted-mean method [6,8], or the inverse-distance-weighted interpolation algorithm [19], or a Gauss distance weighting function [12]. These approaches can perform well for neighbors in the horizontal direction. However, they cannot perform effectively in the vertical direction since simple smoothing functions cannot sufficiently reflect the rapid vertical variation in the WV. Some researchers employed an exponential relationship to express the vertical relationship between the WV in nearby vertical voxels, which is under the approximation that the WV follows the ideal-gas law distribution along the vertical direction [38,51]. The exponential function performs much better than the simple weighted-mean method on constraining the vertical distribution of WV [52]. As a result, it has become the most commonly used vertical constraint model in GNSS tomography, while the exponential function is determined by the WV scale height (H), which is a key parameter of the local troposphere and is empirically determined. The vertical constraint model typically used in the solution process of GNSS tomography is not only to overcome the ill-posed problem in the inversion of the tomographic equations but also to help retrieve the vertical structure of tropospheric WV [53]. Several investigations [12,32,33,35] have shown that the performance of GNSS tomography largely depends on the accuracy of the vertical constraint conditions used in the model, especially in the case that a ground-based GNSS network is relatively flat in the height distribution of all the GNSS receivers, i.e., the geographical distribution is not optimal because of the difficulty in inverting the vertical structure of WV [8,54]. The WV varies rapidly in time and space within the troposphere, and its vertical distribution also correlates with the WV state caused by some meteorological factors, e.g., temperature and WV pressure [55]. Its vertical distribution may not always strictly follow an exponential law under different WV states, especially in the lower troposphere where WV is most active and abundant. Therefore, using the exponential relation to estimate the WV content of an empty voxel by its nearby vertical voxels under different WV states may lead to large errors.

Since the information in the variation of WV cannot be accurately detected at a high spatial and temporal resolution with current techniques, the accuracy of its fitting model may not be as good as desired under various weather conditions. One of the main issues in current GNSS tomography is the low accuracy of the vertical constraint model [13]. In this study, based on the spatio-temporal characteristics of the WV parameter, i.e., IRPWV, under different WV states, a temporal model with six sets of coefficients has been developed for the six classified WV states over the Hong Kong region. The model of IRPWV is applied as a new vertical constraint condition in the GNSS tomographic system model, and the selection of its coefficients is based on the WV state at the tomographic epoch over the tomographic region.

This paper is organized as follows. Section 2 introduces the methodology of voxel-based GNSS tropospheric WV tomography; Section 3 describes the data and processing for the GNSS tomographic experiments performed in this study, followed by the newly proposed vertical constraint model; Section 4 presents the results of the GNSS tomographic experiments and Section 5 presents concluding remarks.

2. Voxel-Based GNSS Tropospheric WV Tomography

After the retrieving of GNSS-derived slant water vapor (SWV) of each GNSS signal ray path, the continuous tomographic domain is needed to uniformly discrete into a set of three-dimensional voxels over the research area, and the continuous and varying distribution of WV over the research area is discretized by the voxels under the assumption that the WV density (WVD) in each discrete voxel is invariable, i.e., a constant, during a tomographic period, and the small influence of geometric delay in signal propagation is neglected. By discretizing the tomographic domain into a finite number of voxels, an integral signal ray path is divided into several line segments, and the full GNSS-derived SWV carried by the integral signal ray path is divided into several segments by the voxels along the signal ray path. Each segment represents the WV content of its voxel, the partial SWV, which is obtained by combining the WVD of the voxel with its corresponding signal intercept. Each partial SWV along the path of the signal ray forms a basic tomographic observation equation in discrete form. The three-dimensional WV content in the tomographic domain can be reconstructed from a large number of such signal rays during the tomographic epoch.

From the above procedure, the GNSS-derived SWV along the path of the GNSS signal in the tomographic domain during the tomographic epoch can be discretized and expressed in terms of the tomographic observation equation as

$${\mathrm{SWV}}^{sta-sat}={\displaystyle \sum }_{i,j,k}^{m,n,l}{d}_{i,j,k}^{sta-sat}·{\mathrm{WVD}}_{i,j,k}+v$$

(1)

where ${\mathrm{SWV}}^{sta-sat}$ is the GNSS-derived SWV along the signal ray paths of the station and the satellite ($sta-sat$); $d_{i,j,k}^{sta-sat}$ (unit: km) and ${\mathrm{WVD}}_{i,j,k}$ (unit: g/m³) represent the distance traveled by the signal ray and the WVD within the (ith, jth, kth) voxel, respectively; the number of the voxels is m $\times$ n $\times$ l; v is the residual term.

During the tomographic period, according to the effective signals (i.e., the signal passes through the top of the tomographic domain), the matrix form of the tomographic observation equations, which consist of all effective GNSS-derived SWV observations from all satellites observed at the GNSS station, can be written

$$A·X=b$$

(2)

where A and b are the coefficient matrix composed of the signal intercept and the column vector made up of the effective GNSS-derived SWV observations, respectively; X is the column vector of the unknown parameters that are the WVD values in each voxel.

In general, the WVD parameter in a voxel can be theoretically obtained by solving the above linear equation using least-squares estimation. However, due to the ill-posed problem of the tomographic observation equations, i.e., Equation (2) is not uniquely solvable, the horizontal and vertical constraint conditions are usually combined with Equation (2) to establish the entire GNSS tomographic equations in the following matrix form.

$$\left[\begin{array}{c}A\\ A_H\\ A_V\end{array}\right]·X=\left[\begin{array}{c}b\\ 0\\ 0\end{array}\right]$$

(3)

where $A_H$ is the coefficient matrix of horizontal constraint equations (e.g., the Gauss-weighted function [12] was adopted in this study); $A_V$ is the coefficient matric for vertical constraint equations, which will be introduced in Section 3.

To analyze the impact of using different vertical constraints on the GNSS tomographic result, for all tomographic schemes performed in this study, the established model of GNSS tomography was resolved under the same configurations: utilizing the singular value decomposition (SVD) method, the same horizontal constraint, and the same weighting strategy, i.e., the three types of the observations in Equation (3) were set to equal weights (the identity matrix). The only difference between different tomographic schemes is the use of different vertical constraints.

3. GNSS Tomographic Experiment

3.1. Data Description and Data Processing

Hong Kong, China, was chosen as the study area. The distribution of the Satellite Positioning Reference Station Network (SatRef) and the King’s Park radiosonde station in the tomographic region are shown in Figure 1. SatRef provides the GNSS observations for the tomographic experiments, and the radiosonde station provides meteorological profiles that are used for modeling and also as references for tomographic results.

3.1.1. GNSS Data

The Survey and Mapping Office (SMO) of Lands Department in Hong Kong, China, has applied the Global Navigation Satellite System (GNSS) to develop the local satellite positioning system, i.e., SatRef, the GNSS reference station network, which consists of 18 Continuously Operating Reference Stations (CORS), including 16 reference stations and 2 Integrity Monitoring stations evenly distributed in Hong Kong. The stations receive GNSS satellite data round-the-clock and send it to a data center for further processing and analysis. Some stations are equipped with an automatic meteorological device to record temperature, pressure, and relative humidity. The largest altitude difference among the 17 GNSS stations selected for this study is about 330 m, implying that the GNSS stations constitute a relatively flat network of GNSS reference stations.

The GNSS observations for the whole year 2019 from 17 GNSS stations were preprocessed by the Bernese software (version 5.2), and the two epochs, 00 and 12 UTC on each day, were selected as the tomographic epochs for the testing of this study, since the observing time of the radiosonde data (used as the reference) are available only at these two epochs. The following processing strategies were adopted in the usage of the software: the precise point positioning (PPP) approach [56]; a 10° cut-off elevation angle; the global mapping function (GMF) [57]; a 30-s sampling rate of GNSS data; a 1-h temporal resolution of zenith total delay (ZTD) and horizontal gradient estimates. The wet delay gradient was obtained by subtracting the dry delay gradient (the mean of 12-h tropospheric gradients) from the tropospheric gradient [58]. The ZTD and the wet delay gradient were mapped with the mapping function and then transformed into SWV by the conversion factor [5] and expressed as

$$\mathrm{SWV}=\mathsf{\Pi }·\left({\mathrm{MF}}_w\left(e\right)·\left(\mathrm{ZTD}-\mathrm{ZHD}\right)+{\mathrm{MF}}_g\left(e\right)\left({\mathrm{G}}_w^{NS}\cos \left(\alpha \right)+{\mathrm{G}}_w^{EW}\sin \left(\alpha \right)\right)\right)+\mathsf{\epsilon }$$

(4)

where $\mathsf{\Pi }$ is the conversion factor; e and $\alpha$ are the elevation and azimuth angles of the ray path of the GNSS signal, respectively; ${\mathrm{MF}}_w$ and ${\mathrm{MF}}_g$ are the wet mapping function (GMF) and horizontal gradient mapping function [59], respectively; ${\mathrm{G}}_w^{NS}$ and ${\mathrm{G}}_w^{EW}$ are the wet gradients in the direction of north–south (NS) and east–west (EW), respectively; $\mathsf{\epsilon }$ is the residual unmodeled delay along the ray path [60]; ZHD is the zenith hydrostatic delay, which can be accurately estimated by the empirical Saastamoinen model [60] and using the surface pressure at the station as the input of the model.

3.1.2. Radiosonde Data

The sounding data over the 12 years from 2008 to 2019 at the King’s Park radiosonde station were downloaded from the Integrated Global Radiosonde Archive (IGRA). The data with a 12-h temporal resolution (observed at UTC 00:00 and 12:00) can provide high-accuracy various meteorological profiles, including pressure, temperature, WV pressure, relative humidity, etc.

According to the equation of state for WV, WVD (${\rho }_w$, unit: g/m³) can be obtained by a function of the WV partial pressure (e, unit: hPa) and temperature (T, unit: K) [61], and then PWV (unit: mm) can be calculated by the integral of WVD, expressed as

$${\rho }_w=\frac{e}{RT}$$

(5)

$$\mathrm{PWV}=\frac{1}{{\rho }_v}{{\displaystyle \int }}_{h_s}^{h_{trop}}{\rho }_{w}dh$$

(6)

where R is the specific gas constant of the WV (R = 0.4615 $\mathrm{J}/\left(\mathrm{K}·\mathrm{g}\right)$); ${\rho }_v$ is the water density (${\rho }_v=1$ g/cm³); $dh$ is the increment step (unit: km); $h_s$ and $h_{trop}$ are the heights of the ground surface and the tropopause (unit: km), respectively.

Based on the above two equations, the radiosonde-derived WVD (RS-derived WVD) profiles and PWV on each day of the year (doy) from 2008 to 2019 were calculated. Data from the first 11 years was used as the sample data to develop the vertical constraint model in Section 3.3, and data from the last year was used as a reference to validate the tomographic results.

3.2. Division Strategy for Tomographic Voxel

The tomographic region is from 113.81°E to 114.42°E and from 22.13°N to 22.57°N. A 0.1° interval in both longitude and latitude directions was adopted for the horizontal grid division, which yielded a 6 $\times$ 5-floor plan grid (the longitudinal and latitudinal axes were 6 and 5, respectively, as shown in the left pane of Figure 2); the troposphere over the research area is divided into 10 non-uniform layers from 0 km to the tropopause (about 10 km, where the WVD is close to 0 g/m³) along the vertical direction. Specifically, the three different height intervals were 0.6 km (below 3 km), 1 km (between 3 km and 6 km), and 2 km (above 6 km), as shown in the right pane of Figure 2. By combining the horizontal grid division and vertical layering, the tomographic domain was discretized into 6 $\times$ 5 $\times$ 10 voxels.

3.3. Vertical Constraint Conditions

3.3.1. Traditional Vertical Constraint Condition

The traditional vertical constraint condition used in GNSS tomography is an exponential function expressed as

$${{\rho }_w}_j={{\rho }_w}_i{e}^{\frac{-\left(h_j-h_i\right)}{H}}$$

(7)

where ${{\rho }_w}_i$ and ${{\rho }_w}_j$ are the WVD of the ith and jth layers, respectively; $h_i$ and $h_j$ are the height of the ith and jth layers, respectively; H is the WV scale height (unit: km), which can be calculated from PWV and the ground surface WVD (SWVD, ${{\rho }_w}_s$) by

$$H=\frac{\mathrm{PWV}}{{{\rho }_w}_s}$$

(8)

The H of a local troposphere is empirically determined and usually obtained using the historical sounding data of the local radiosonde station [62]. In the application of the traditional vertical constraint, three types of H parameters are generally taken, namely, a constant, a periodic function value, and a near real-time value, which are calculated as follows: (1) Based on the sample data in Section 3.1.2, the 11-year time series of H was calculated according to Equation (8), and the mean value of H ($H_{mean}$) was about 2.4 km. (2) According to the periodic nature of H [63,64], the second-order trigonometric periodic function ($H_f$) was used to fit the time series, and its fitted coefficients were $a_{H0}$ = 2.6149, $a_{H1}$ = 0.0487, $b_{H1}$ = 0.0607, $a_{H2}$ = −0.0253, and $b_{H2}$ = −0.0330. (3) Some of the 17 GNSS stations were equipped with an automatic meteorological device for recording temperature, pressure, and relative humidity. Therefore, a near real-time H ($H_{real-time}$) was obtained by a near real-time SWVD and GNSS-derived PWV [53]. Note that if there were more than one near real-time H at a tomography epoch, then their mean would be taken as the final result.

3.3.2. Newly Proposed Vertical Constraint Model

IRPWV (unit: 1/km), derived from our previous work [65], is a WV parameter of which the vertical distribution represents the variation in WV along the vertical direction over a site in the troposphere, and it is expressed as

$${\mathrm{IRPWV}}_h=\frac{{{\rho }_w}_h}{\mathrm{PWV}}$$

(9)

where ${\rho }_{w_h}$ is the WVD (unit: g/m³) at the height $h$ (unit: km).

The relationship between the WVDs at any two heights over the site is expressed in terms of their corresponding IRPWV as

$${\rho }_{w_j}={\rho }_{w_i}\frac{{\mathrm{IRPWV}}_j}{{\mathrm{IRPWV}}_i}$$

(10)

where ${\rho }_{w_i}$, ${\rho }_{w_j}$, ${\mathrm{IRPWV}}_i$ and ${\mathrm{IRPWV}}_j$ are the WVDs and IRPWVs at $h_i$ and $h_j$, respectively.

Based on the analysis conclusions that there is a strong correlation between the vertical structures of IRPWV and their corresponding WV states and the periodic variations in IRPWV, first, the vertical structures of IRPWV were classified into six vertical structures according to their WV states relative to the six WV states. Note that the six WV states are the states of maximal disturbance, sub-disturbance, normal, normal, sub-saturated, and saturated. Then, the spatiotemporal IRPWV model was developed at a specific height (i.e., the tomographic height layers of this study) for each of the six WV states. The spatio-temporal IRPWV models are represented in a unified form as

$${\mathrm{IRPWV}}_{h,r}=a_{0,h,r}+a_{1,h,r}\cos (\mathrm{doy}·w)+b_{1,h,r}\sin (\mathrm{doy}·w)+a_{2,h,r}\cos (\mathrm{doy}·2w)+b_{2,h,r}\sin (\mathrm{doy}·2w)$$

(11)

where $a_{0,h,r}$, $a_{1,h,r}$, $b_{1,h,r}$, $a_{2,h,r}$ and $b_{2,h,r}$ are the periodic coefficients; the subscript h and r are the indexes of the height layers and six WV states, respectively; $w=2\mathsf{\pi }/365.25$.

The spatiotemporal IRPWV model was employed as a new vertical constraint model in the GNSS tomographic system equation, and its coefficients were estimated and are shown in Appendix A. The procedure for applying the new vertical constraint model to GNSS tomography is as follows: Step 1—according to the WV state at the tomographic epoch (i.e., the doy parameter), the coefficients of the new vertical constraint model are selected from the corresponding WV states (

Table A2. Six sets of coefficients of the new vertical constraint model for the six classified WV states at the 10 tomographic layers corresponding to six PWV ranges (from 1–6) at the radiosonde station in Hong Kong, China.

Coefficient of PWV Range		Tomographic Layer
		1	2	3	4	5	6	7	8	9	10
1	$a_0$	0.478	0.363	0.249	0.183	0.130	0.088	0.047	0.026	0.013	0.005
	$a_1$	0.030	0.011	−0.018	0.000	0.007	−0.005	−0.011	−0.005	0.000	0.000
	$b_1$	−0.002	−0.017	−0.006	0.001	0.012	0.002	0.000	0.001	0.002	0.001
	$a_2$	0.020	0.003	−0.016	0.020	0.003	−0.012	−0.007	−0.002	0.000	0.000
	$b_2$	0.008	0.001	0.011	0.005	−0.001	−0.004	−0.005	0.000	−0.001	0.000
	$a_3$	−0.014	−0.017	−0.005	0.016	0.008	−0.002	−0.002	0.003	0.002	0.000
	$b_3$	0.012	0.012	0.013	−0.006	−0.010	−0.011	−0.003	−0.001	−0.001	0.000
2	$a_0$	0.409	0.343	0.274	0.215	0.147	0.094	0.049	0.027	0.011	0.003
	$a_1$	0.012	0.019	0.030	0.016	0.009	−0.019	−0.021	−0.010	−0.003	−0.001
	$b_1$	−0.012	−0.006	0.010	0.011	0.010	0.002	−0.007	−0.003	−0.001	0.000
	$a_2$	0.016	0.008	0.007	0.009	−0.005	−0.014	−0.008	0.001	0.001	0.000
	$b_2$	0.006	0.007	0.010	0.011	0.003	−0.006	−0.008	−0.003	−0.002	0.000
	$a_3$	−0.003	−0.005	−0.001	0.005	0.000	−0.004	0.002	0.001	0.001	0.000
	$b_3$	0.001	0.003	0.005	−0.001	−0.004	−0.003	0.001	0.000	−0.001	0.000
3	$a_0$	0.387	0.338	0.270	0.216	0.153	0.099	0.052	0.028	0.012	0.004
	$a_1$	0.008	0.027	0.032	0.021	0.011	−0.018	−0.020	−0.014	−0.005	−0.001
	$b_1$	−0.013	0.002	0.012	0.011	0.008	−0.001	−0.007	−0.004	−0.001	0.000
	$a_2$	0.011	0.006	0.007	0.007	−0.001	−0.011	−0.005	−0.001	0.001	0.000
	$b_2$	−0.002	0.002	0.012	0.015	0.000	−0.007	−0.007	−0.002	−0.001	0.000
	$a_3$	0.003	0.001	0.002	0.004	0.000	−0.005	−0.001	0.000	0.000	0.000
	$b_3$	−0.001	−0.001	0.004	0.001	−0.003	−0.004	0.002	0.002	0.000	0.000
4	$a_0$	0.365	0.325	0.265	0.218	0.157	0.105	0.056	0.031	0.014	0.004
	$a_1$	0.003	0.025	0.037	0.025	0.017	−0.014	−0.022	−0.016	−0.007	−0.002
	$b_1$	−0.001	0.012	0.012	0.007	0.004	−0.005	−0.010	−0.005	−0.002	0.000
	$a_2$	0.006	0.003	0.006	0.005	−0.001	−0.009	−0.006	0.000	0.001	0.001
	$b_2$	0.007	0.013	0.010	0.007	0.000	−0.008	−0.007	−0.003	−0.001	0.000
	$a_3$	−0.001	−0.002	0.000	0.003	−0.002	−0.002	0.001	0.001	0.000	0.000
	$b_3$	0.000	0.003	0.002	−0.001	0.000	−0.007	0.000	0.002	0.001	0.000
5	$a_0$	0.347	0.310	0.257	0.215	0.157	0.112	0.064	0.037	0.016	0.005
	$a_1$	0.000	0.020	0.031	0.024	0.017	−0.006	−0.016	−0.015	−0.009	−0.003
	$b_1$	0.001	0.016	0.015	0.010	0.006	−0.004	−0.010	−0.006	−0.003	−0.001
	$a_2$	−0.002	−0.004	0.003	0.002	0.000	−0.003	−0.003	0.000	0.000	0.000
	$b_2$	0.006	0.010	0.009	0.006	0.004	−0.007	−0.008	−0.004	−0.001	0.000
	$a_3$	−0.003	−0.005	−0.002	0.000	−0.001	0.003	0.001	0.002	0.000	0.000
	$b_3$	−0.001	−0.001	0.001	0.002	−0.001	−0.003	0.001	0.001	0.001	0.000
6	$a_0$	0.321	0.293	0.244	0.209	0.155	0.121	0.074	0.044	0.020	0.006
	$a_1$	−0.010	0.013	0.022	0.018	0.013	0.004	−0.007	−0.012	−0.009	−0.004
	$b_1$	0.005	0.017	0.013	0.010	0.005	−0.004	−0.008	−0.008	−0.004	−0.001
	$a_2$	−0.007	−0.004	0.001	0.002	0.000	0.004	−0.002	0.000	0.000	0.000
	$b_2$	0.005	0.009	0.006	0.004	0.002	−0.003	−0.006	−0.004	−0.001	0.000
	$a_3$	−0.006	−0.009	−0.005	−0.001	−0.002	0.003	0.004	0.003	0.001	0.000
	$b_3$	0.000	−0.001	−0.002	−0.002	−0.002	−0.002	−0.001	0.002	0.001	0.000

); Step 2—the coefficients and tomographic epoch are input into Equation (11) to calculate the IRPWV value at each tomographic layer; Step 3—the WVDs at any two tomographic layers along the vertical direction are expressed by Equation (10), and then the new vertical constraint equation is obtained.

4. Test Results of Using Different Vertical Constraints

In this section, WVDs resulting from the GNSS tomographic technique (hereinafter referred to as GNSS-Tom) using different vertical constraints are compared against the reference values. The GNSS-Tom experiments, with four sets of vertical constraints, are defined as different tomographic schemes, hereafter referred to as SCH1, SCH2, SCH3, and SCH4, and they are listed in

Table 1. Four tomographic schemes employing two different vertical constraints with different parameters, i.e., SCH1, SCH2, and SCH3 with the traditional vertical constraint using different H parameter values and SCH4 with the new vertical constraint.

Scheme	Vertical Constraint	Type of Parameter
SCH1	${{\rho }_w}_j={{\rho }_w}_i{e}^{\frac{-\left(h_j-h_i\right)}{H}}$	$H_{mean}$
SCH2		$H_f$1
SCH3		$H_{real-time}$
SCH4	${\rho }_{w_j}$$={\rho }_{w_i}\frac{{\mathrm{IRPWV}}_{h_j,r}}{{\mathrm{IRPWV}}_{h_i,r}}$	${\mathrm{IRPWV}}_{h,r}$

¹$H_f=a_{H0}+a_{H1}\cos (\mathrm{doy}·w)+b_{H1}\sin (\mathrm{doy}·w)+a_{H2}\cos (\mathrm{doy}·2w)+b_{H2}\sin (\mathrm{doy}·2w)$.

.

Before comparing to the RS-derived WVD, the GNSS-Tom-derived WVD in each tomographic voxel is the mean of the WVDs in that voxel, and therefore the RS-derived WVD profile needs to be interpolated to the central position of the voxel. The interpolation method used divides the intercepted RS-derived PWV within a voxel by the height of the voxel. From the left pane of Figure 2, it can be seen that the radiosonde station is located close to the center of the grid which is in row 4 from the top to the bottom along the latitudinal direction and in column 4 from the left to the right along the longitudinal direction. Hence, the GNSS-Tom-derived WVDs in the voxels located in row 4 and column 4 of the tomographic layers 1–10 were directly used to compare against the interpolated RS-derived WVDs in the same voxel.

The GNSS-Tom-derived WVD resulting from the four tomographic schemes for each doy of 2019 is compared against the RS-derived WVD. Four statistical indicators, including relative error (RE), bias, root mean square error (RMSE), and skill score (SS), were used to evaluate the performance of the four schemes. The four metrics are defined as follows.

$$\mathrm{RE}=\frac{\left|valu{e}_{model}^i-valu{e}_{reference}^i\right|}{valu{e}_{reference}^i}·100\%$$

(12)

$$\mathrm{bias}=\frac{1}{N}{\displaystyle \sum }_{i=1}^N\left(valu{e}_{model}^i-valu{e}_{reference}^i\right)$$

(13)

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\left(valu{e}_{model}^i-valu{e}_{reference}^i\right)}^2}$$

(14)

$$\mathrm{SS}=100·\left(1-\frac{{\mathrm{RMSE}}_1}{{\mathrm{RMSE}}_2}\right)$$

(15)

where the superscript i is the index of the sample data; the two subscripts denote the model and reference; N is the number of samples; ${\mathrm{RMSE}}_1$ and ${\mathrm{RMSE}}_2$ are the RMSE of experiments 1 and 2, respectively.

RE is defined as the percentage of the absolute difference between the model value and the reference value to the reference value, and in general, the smaller the RE value, the better the result. A positive SS indicates an improvement of experiment 1 over experiment 2, and the larger the SS value, the better the performance. When SS = 1, the best performance is achieved. A negative SS indicates that the application of the method (or data, or parameter, etc.) in experiment 1 had a detrimental impact. The daily RE, monthly and annual bias, and RMSE of the GNSS-Tom-derived WVD for each of the 10 tomographic layers during the 2019 tomographic experiment period, as well as the SS values for the monthly RMSEs, are all presented in the next two sections.

4.1. RE of Layered WVD

The resulting daily RE of the WVD for each of the four schemes at each of the 10 tomographic layers during 2019 is shown in Figure 3, and the information on the statistical proportion of the daily RE is given in

Table 2. The statistical proportion of daily RE (%) (RE < 30% and RE > 50%) of WVD resulted from the four tomographic schemes at each of the 10 tomographic layers during 2019. Note that for the RE > 50% case, because the values of the statistical proportions at tomographic layers 1 and 2 are small, the values with one decimal place were adopted, as shown in the blue underline.

Tomographic Layer (km)	RE < 30%				RE > 50%
	SCH1	SCH2	SCH3	SCH4	SCH1	SCH2	SCH3	SCH4
1 (0–0.6)	86	85	99	94	2.3	1.9	0.0	1.6
2 (0.6–1.2)	85	77	94	96	0.8	1.5	0.1	0.6
3 (1.2–1.8)	75	75	78	92	2	2	2	2
4 (1.8–2.4)	75	75	76	86	5	5	5	5
5 (2.4–3)	70	71	70	74	12	12	11	14
6 (3–4)	63	63	66	64	26	28	24	25
7 (4–5)	52	51	56	55	38	41	36	33
8 (5–6)	41	39	44	45	48	52	45	38
9 (6–8)	23	19	21	32	68	73	69	47
10 (8–10)	14	9	9	21	79	86	82	57

.

It can be seen from Figure 3 that, for all the schemes, the large daily RE values (e.g., RE > 50%) are mainly distributed at the upper layers and vary with the month. The distribution of the four schemes is all the sparest in the summer and the opposite in the winter. The RE distributions of SCH1, SCH2, and SCH3 are similar, while compared to the three schemes, the large daily REs of SCH4 are more sparsely distributed in the troposphere throughout 2019.

For a more detailed analysis, the daily REs were classified empirically into two categories: RE < 30% and RE > 50%, which can distinguish the differences between the four schemes more clearly than the other categories. The statistical proportions of RE < 30% and RE > 50% of the four schemes are shown in Table 2. For RE < 30%, the higher the proportion, the better the result; whilst for RE > 50%, the higher the proportion, the worse the result.

It can be seen from Table 2 that, for all schemes, the statistical proportion of RE < 30% decreases with the increase in height layer from 1–10, which is completely different from the statistical proportion of RE > 50%. For the RE < 30% case, at layers 1–5, on average, SCH4 was 13% and 15% better than SCH1 and SCH2, respectively; it was about 10% higher than SCH3, except for layer 1, where it was 5% lower. At layers 8–10, SCH4 outperformed all the other three schemes by at least 49%. At layers 6–8, the four schemes showed similar performance. SCH4 performed slightly better than SCH1 and SCH2, and it performed slightly better than SCH3 at layer 8, but not at layers 6–7. For the RE > 50% case, the results of SCH1, SCH2, and SCH4 decreased at first and then increased with the layers from 1–3, with the returning point between layers 2 and 4 (about 1.2–1.8 km), while the results of SCH3 only showed an upward trend. In more detail, SCH4 performed better than both SCH1 and SCH2, with SCH2 slightly worse than SCH1. The four schemes performed similarly at layers 3–6. At the middle and upper layers (4–10 km), SCH4 was smaller than the other three schemes, and its reduction of statistical proportion was decreased by at least 3–22% with the increase in height layer.

Overall, the above results indicated that the results of the four tomographic schemes all performed better in summer, especially at lower altitudes than at higher altitudes. The daily RE was greater at the upper layer than at the lower layer and in winter than in summer, mainly because a small difference between a GNSS-Tom-derived WVD and its reference may result in a large absolute error relative to the reference, which also leads to a large RE value in the height/season where/when the WVD is inherently tiny, especially at the upper layers in winter.

4.2. Bias and RMSE of Layered WVD

The monthly biases and RMSEs of the WVDs at each of the 10 tomographic layers resulting from each of the four schemes, as well as the skill scores (SSs) of the monthly RMSE, are shown in Figure 4 and Figure 5, respectively, during the experiment period of 2019. The results for the annual biases and RMSEs are shown in Figure 6.

4.2.1. Monthly Bias and RMSE

As can be seen from Figure 4a, the monthly biases of SCH1, SCH2, and SCH3 from 1–10 layers all presented a consistent wave shape during the 12 months of 2019, and their maximum biases occurred at layers 2–4 (about 1–2 km). All the 12 monthly biases of SCH3 at layer 1 were closer to zero than the other two schemes. Compared to the first three schemes, the monthly biases of SCH4 were around zero and slightly fluctuated with height at layers 1–10 in the 12 months, slightly varying with the variation in height. It can be seen from Figure 4b that, at layer 1, the monthly RMSEs of SCH3 were all significantly smaller than those of the other three schemes during the 12 months, wherein those of SCH4 were smaller than the remaining two schemes. At other layers, the results of SCH1, SCH2, and SCH3 were similar, while the monthly RMSEs of SCH4 were almost smaller than those of the other three schemes, except for layers 5–6, where all the schemes performed similarly. Compared to SCH1, SCH2, and SCH3, the monthly RMSEs of SCH4 were improved by at least 37%, 28%, 23%, and 21% at the lower layers (0.6–3 km), and 32%, 30%, 15%, and 18% at the upper layers (4–10 km) in the winter, spring, summer, and autumn, respectively.

4.2.2. Skill Score of Monthly RMSE

The monthly RMSEs of SCH4 and the other three schemes at each of the 10 layers were compared, and the performance of the vertical constraint model used in a tomographic scheme can be evaluated by the SS. According to Equation (15), the SS values of SCH4 relative to SCH1, SCH2, and SCH3 were obtained, namely SCH14, SCH24, and SCH34, respectively, and the results are shown in Figure 5. As we know, the only difference between the four tomographic schemes is in the vertical constraint strategy. Hence, a positive SS value indicates that the use of the new vertical constraint model proposed in this study in SCH4 is an improvement over the other three schemes. Conversely, a negative SS indicates that the assimilation of the new vertical constraint model in SCH4 has a detrimental effect. In our tomographic experiments, if the SS value was in the range of −10 to 10, then the performance of the two associated results was considered the same.

The SS values of SCH14 (left) and SCH24 (middle) in Figure 5 indicated that all the results of SCH4 were superior to those of SCH1 and SCH2. While the SSs of SCH34 revealed that SCH3 performed best at the first layer, its results at the other layers were similar to those of SCH1 and SCH2. At layers 2–5 and 8–10, the performances of SCH4 were better than those of SCH3. According to the statistics of the results from Figure 5, the numbers of SS values above 10 for SCH14, SCH24, and SCH34 were about 83%, 87%, and 64%, respectively. More specifically, the SS values at layers 1–4 of SCH14 and SCH24 were almost all above 10 with the mean of 30 and 32, respectively. Their values at layers 7–10 and 5–6 were approximately 25 and 33, and 11 and 10, respectively. The means of the SS values at layers 2–4 and 9–10 of SCH34 were about 28 and 40, respectively, and the values at layers 5, 7, and 8 were all about 10 on average.

4.2.3. Annual Bias and RMSE

From Figure 6a, compared with the other three schemes, the annual bias of SCH4 was closest to zero from the first to the last layers, and its largest absolute bias was at layers 2–3, with a value of about −0.5 g/m³. The other three schemes all showed a notable systematic deviation, and the signs of systematic deviations at layers 2–5 were the opposite of those of layers 7–9. Although all four schemes showed a large annual bias between layers 2 and 3 (about 1 km), the annual bias of SCH4 significantly decreased by 1.1–1.5 g/m³ compared to those of the other three schemes. In Figure 6b, the annual RMSEs of SCH4 at layers 1–6 and layers 7–10 have improvements of 23–34% and 13–42%, and 9–37% and 10–47%, respectively, over those of SCH1 and SCH2, while they outperformed those of SCH3 at layers 2–6 and 7–10 by 12–35% and 5–45%, respectively.

The monthly and annual layered results indicated that, first, the tomographic results of SCH4, which used the new model, were more stable than the other three schemes, which used the traditional vertical constraint model. The new model is a linear relationship between WVDs at different heights along the vertical direction. Unlike the exponential relationship expressed by the traditional vertical model, the new model can well reflect the relationship of WVs at different heights, even in the case that the WVD varies rapidly at a specific height in the troposphere, without the occurrence of very unreliable WV values. Secondly, the accuracy of the new model result significantly outperformed the other three schemes at both the lower layers (0.6–3 km) and upper layers (4–10 km) in all 12 months. The selection of the coefficient for the new model is based on the WV state at the tomographic epoch. As a result, its tomographic result agrees with the water vapor distribution better than the traditional model. Thirdly, at the middle layers (3–4 km), the new model performs similarly to the traditional models, which is mainly because the signs of the systematic errors in the tomographic results of SCH1, SCH2, and SCH3 at the lower layers are completely different from those at the upper layers, i.e., if the systematic errors are negative values at the lower layer, then they will be positive values at the upper layer. Thus, the systematic errors at the middle layers are most likely to be close to zero, which leads to the same performance of SCH1, SCH2, and SCH3 as SHC4 at these layers. It should be noted that the good result of SCH3 at layer 1 (ground surface) is due to its usage of near real-time H parameter calculated from real-time WV. At the other layers, however, the results of SCH3 are comparable to SCH1 and SCH2. Overall, the significantly better performance of SCH4 than the other three schemes suggests the effectiveness of the new model.

5. Conclusions

With the development of GNSS tomography for nearly 20 years, one of the main challenges of GNSS tomography has been to solve ill-conditioned model equations. Vertical constraint models are typically used in the solution process and play an important role in the quality of the GNSS tomography, in addition to resolving ill-posed problems in the system equations. In this study, based on the IRPWV parameter, a new vertical constraint model with six sets of coefficients was developed for the six classified WV states over the tomographic region for improving the accuracy of GNSS tomography, since the new model can provide a good representation of the vertical variation in the WV of the local troposphere under various WV states.

To test the performance of the new vertical constraint model on GNSS tomography, tomographic experiments with four tomographic schemes were carried out using GNSS-derived WV from 17 GNSS stations in Hong Kong, China throughout 2019. The first tomographic scheme used the new model (scheme 4, namely SCH4), and the other three schemes employed the traditional vertical constraint with three different types of water vapor scale height (H), i.e., a constant (scheme 1, namely SCH1), a periodic function (scheme 2, namely SCH2), and a near real-time value (scheme 3, namely SCH3). The results of the four schemes were compared against their corresponding radiosonde-derived WV density (WVD), and the results of the four schemes were compared on three different time scales: daily, monthly, and annual.

The results showed that the daily relative errors (RE) of all four schemes at the lower layer were smaller than those at the upper layer, especially in summer. This is mainly because the RE is sensitive to small values; even a tiny deviation will result in a large RE. For the case of the daily RE of WVD of less than 30%, the results of SCH4 were 13% and 15% on average better than those of SCH1 and SCH2 at the lower layers (below 3 km), respectively; they were about 10% higher than those of SCH3 at the same layers, except for the ground surface, where they were 5% lower. The four schemes showed similar performance at the middle layers (4–5), while at the upper layers (5–10 km), SCH4 outperformed all the other three schemes by at least 49%. It should be noted that the usage of the real-time WV value is the only reason for the higher accuracy of SCH3 at layer 1; hence, this will be taken into account in the development of an improved vertical constraint model in future.

The monthly biases of SCH4 were almost zero at all layers (1–10) within all 12 months, while that of SCH1, SCH2, and SCH3 were smaller than the reference at the lower layers and larger at the upper layers, which is mainly due to the systematic deviation of the traditional models. In contrast, the SCH4 result had no obvious deviation from the reference because the selection of the new model is based on the characteristics of water vapor at different heights under different water vapor states. The monthly RMSEs of SCH4 at the lower layers (0.6–3 km) were reduced at least by 37%, 28%, 23%, and 21% and at the middle and upper layers (4–10 km) by 32%, 30%, 15%, and 18% in the winter, spring, summer, and autumn, respectively. Additionally, the monthly skill scores of SCH4 can reflect the performance of the new model. By comparing the monthly layered RMSEs of SCH4 relative to those of SCH1, SCH2, and SCH3, the statistics of skill scores above 10 in the new model accounted for 83%, 87%, and 64%, respectively. Similar results were observed in the annual results as well as the monthly results. Although all four schemes showed the maximal annual absolute bias was at layers 2 and 3 (about 1 km), the annual bias of SCH4 was about −0.5 g/m³, and it significantly decreased by 1.1–1.5 g/m³ in comparison to those of the other schemes. The annual RMSEs of SCH4 compared with SCH1 and SCH2 at layers 1–6 (the ground surface to 3 km) decreased by 23–34% and 13–42%, respectively, and at the other layers (3–10 km) by 9–37% and 10–47%, respectively. Compared with SCH3, the annual RMSEs of SCH4 at layers 2–6 (0.6–3 km) and 7–10 (4–10 km) decreased by 12–35% and 5–45%, respectively.

All evaluation results demonstrate the better performance of the new model than the traditional ones. Based on the spatiotemporal distribution of the WV at different heights for different WV states, the new model obtained from the classification statistical analysis is more refined. In its application, the selection of an appropriate new vertical constraint model is based on the state of the WV at the current tomographic epoch over the tomographic region. The better performance of the new model compared to the traditional model, which is roughly a vertical constraint condition with a simple empirical value for the H parameter, suggests that the new model better reflects the vertical distribution of the WV over the tomographic region at the tomographic epoch. Our future work will focus on the construction and testing of new models for different weather conditions or regions.