**1. Introduction**

Since the concept of GNSS meteorology was first proposed by Bevis et al., the water vapor information derived from GNSS has drawn increasing attention in the meteorological and GNSS communities [1]. The precipitable water vapor (PWV), which refers to the height of an equivalent column of water vapor [2], has been widely validated to achieve mm-level accuracy based on the conversion of GNSS zenith tropospheric delay [3]. Further, the three-dimensional water vapor information can also be inversed by using GPS signals as scanning rays in the research area, which is called water vapor tomography.

Braun et al. first proposed the concept of GPS water vapor tomography [4] and Flores et al. first realized it using the data from the Kilauea network in Hawaii [5]. The research region, covered by ground GPS receivers, is discretized into finite voxels according to its latitude, longitude, and altitude, and the unknown estimated parameter of the voxels are assumed to be constant during a given period. The GPS-derived slant water vapor is regarded as the observations for water vapor tomography.

In modeling the GPS water vapor tomography, it is found that the geometric distribution of the observed signals is an inverted cone due to the fixed structure of GPS sites and

**Citation:** Yang, F.; Wang, J.; Wang, H.; Gong, X.; Wang, L.; Huang, B. Assessment of the Water Vapor Tomography Based on Four Navigation Satellite Systems and Their Various Combinations. *Remote Sens.* **2022**, *14*, 3552. https://doi.org/ 10.3390/rs14153552

Academic Editors: Serdjo Kos, José Fernández and Juan F. Prieto

Received: 21 June 2022 Accepted: 22 July 2022 Published: 24 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

satellites [6]. The direct effect caused by this phenomenon is the presence of tomographic voxels without signal rays passing through, especially at the lower and edge layer of the area of interest. It also makes many voxels be penetrated by only a very small number of signals. From the perspective of the water vapor tomography model, it often leads to a large number of zero elements appearing in the coefficient matrix, which becomes a sparse matrix [7]. This is the fundamental cause of the ill-posed problem in GPS water vapor tomography, which seriously restricts its stable and high-accuracy solution. Obviously, enriching the observation equation of the GPS water vapor tomography is an effective way to overcome the above problem by introducing various observation information, and related research has been carried out.

Based on voxel horizontal boundary selection and non-uniform symmetrical division, Chen et al. and Yao et al. proposed an optimized approach of horizontal voxel division to introduce more signal rays penetrated from the top layer into the observation equation [8,9]. The similar effect can be obtained by the method of constructing the tomographic buffer area carried out by Trzcina et al. and Sa et al. [10,11]. These methods are limited to specific tomographic regions and certain experimental periods. Adavi et al. explored how to use the constructed virtual reference sites to augment location-specific GPS observations [12]. The virtual signals were also introduced to the tomography model using the calculated mapping function and ZWD/PWV of corresponding site and the elevation and azimuth of specified virtual satellite [6,13]. Studies have shown that it is a feasible method to incorporate the GPS signal rays passed through form the side face into the tomography model; for example, Zhao et al. constructed the unit scale factor for these signals using the radiosonde and reanalysis data [14], Zhang et al. and Hu et al. established the height factor models adapted to these signals from side face [15,16]. In addition, Zhao et al. tried to extend the observations of GPS sites outside the tomographic region into the tomography model based on the GPT2w and TMF models [17,18]. The above methods all rely on external data or models and tend to introduce new error for the observation information. On the other hand, some have attempted to add multi-source observation information from various sensors into the GPS tomography model, such as the COSMIC occultation data [19], the GNSS-R data [20], the InSAR data [21,22], WRF output data [23], LEO constellationaugmented data [24], PWV data derived from FY-3, and MODIS [25,26]. However, the spatiotemporal resolution, availability, and consistency with the tomographic region are the factors that seriously restrict the fusion of the above observations into the tomography model.

It is more reasonable and convenient to construct the GNSS water vapor tomography model together with the observations from GPS and the other three satellite navigation systems. Bender et al. simulated GPS, GLONASS, and Galileo data and introduced the method for obtaining three-dimensional water vapor information by tomography technique in multi-satellite systems [27]. Wang et al. compared the tomographic accuracy of BDS and GPS based on simulated data, and showed that the result using 9 satellites of BDS is basically comparable to that of GPS [28]. Xia et al. and Benevides et al. carried out the water vapor tomography experiments of GPS combined with GLONASS and GPS combined with Galileo in Hongkong and Lisbon regions, respectively [29,30]. Dong et al. and Zhao et al. utilized the measured data derived from different numbers of BDS2 satellites and combined it with GPS and GLONASS data to construct the tomography model in Wuhan and Guiyang, respectively [31,32]. With the gradual improvement of Galileo and the opening of BDS-3 services, the above experiments based on simulated data or incomplete satellite data cannot fully reflect the current status of water vapor tomography based on multi-GNSS. Therefore, this paper aims to explore the differences between the four satellite navigation systems and their combination in water vapor tomography, including the modeling process and the reconstructed results.

#### **2. Materials and Methods**

The observations in GNSS water vapor tomography are the slant water vapor (SWV) which can be converted from slant wet delay (SWD) as follows [33]:

$$\text{SWV} = \frac{10^6}{\rho\_{\text{uv}} \times \frac{R}{m\_{\text{w}}} \left(\frac{k\_3}{T\_{\text{w}}} + k\_2 - \frac{m\_{\text{r}}}{m\_d} \times k\_1\right)} \times \text{SWD} \tag{1}$$

where *<sup>ρ</sup><sup>w</sup>* refers to the liquid water density with the unit of <sup>g</sup>·m<sup>−</sup>3; *<sup>R</sup>* <sup>=</sup> <sup>8314</sup> Pa·m3·K−1·kmol−<sup>1</sup> denotes the universal gas constant; *mw* and *md* represent the molar mass of water and the dry atmosphere and their values are 18.02 kg·kmol−<sup>1</sup> and 28.96 kg·kmol−<sup>1</sup> , respectively; *Tm* is the weighted mean temperature, which can be calculated by using surface temperature [34,35]; *k*1, *k*2, and *k*<sup>3</sup> are the empirical physical constants, which are equal to 77.60 <sup>K</sup>·hPa−<sup>1</sup> , 70.4 <sup>K</sup>·hPa<sup>−</sup>1, and 3.739 <sup>×</sup> 105 <sup>K</sup>·hPa−<sup>1</sup> , respectively [36]. After mapping the zenith wet delay (ZWD) and the wet delay gradients into the elevation direction, the SWD can be obtained as follows [37]:

$$\begin{array}{lcl} \text{SWD} &=& m\_{\text{w}}(\text{el}\varepsilon) \times \text{ZWD} \\ &+& m\_{\text{w}}(\text{el}\varepsilon) \times \text{cot}(\text{el}\varepsilon) \times \left(\text{G}\_{NS}^{w} \times \cos(azi) + \text{G}\_{WE}^{w} \times \sin(azi)\right) \end{array} \tag{2}$$

where *mw* indicates the wet mapping function and the global mapping function (GMF) was used in this paper; *ele* and *azi* denote the satellite elevation and azimuth angles, respectively. *G<sup>w</sup> WE* and *<sup>G</sup><sup>w</sup> NS* represent the wet delay gradient parameters in the east-west and north-south directions, respectively. Affected by water vapor along the signal ray, ZWD is the wet component of zenith total delay (ZTD) which is the primary parameter retrieved from GNSS observation. To obtain ZWD, the zenith hydrostatic delay (ZHD) should be subtracted from ZTD [38]. In this paper, the Saastamoinen model is used to calculate the accurate ZHD using the pressure measurements as follows [39]:

$$\text{ZHD} = \frac{0.002277 \times P\_s}{1 - 0.00266 \times \cos(2\varphi) - 0.00028 \times H} \tag{3}$$

where *ϕ* and *H* denote the latitude and geodetic height of the GNSS site, respectively. *Ps* is the measured surface pressure.

In the water vapor tomography, the SWV value is also an integral expression of water vapor along the slant path from the ground receiver and GNSS satellite, given by the following:

$$\text{SWV} = 10^{-6} \cdot \int \rho(s) ds\tag{4}$$

where *<sup>ρ</sup>*(*s*) in <sup>g</sup>·m−<sup>3</sup> denotes the water vapor density and *ds* refers to the path traveled by a satellite signal ray. After discretizing the tomographic region into finite voxels, the observation equation of GNSS water vapor tomography can be established based on the distances of GNSS signal rays crossing the divided voxel and the unknown estimated water vapor density with each voxel. It can be expressed as follows:

$$\text{SWV} = \sum\_{i=1}^{n} d\_i \cdot x\_i \tag{5}$$

where *n* represents the total number of divided voxels in the research region. *di* denotes the distance of signal rays inside voxel *i*, which can be calculated by using the coordinates of GNSS sites and satellites. *xi* refers to the water vapor density of voxel *i*, which is the unknown estimated parameter.

In water vapor tomography, two types of constraints are widely used in tomographic modeling along with the observation equation, one is the horizontal constraint and the other is vertical constraint, since a spatial relation exists between water vapor in a specific voxel and its surrounding ones. For the horizontal constraint, it assumes that the distribution of water vapor density is relatively stable in the horizontal direction within a small region, and represents the relationship between the water vapor density of a certain voxel and those of its adjacent voxels in the same layer. For the vertical constraint, it refers to the exponential relationship between the water vapor density of voxels for two consecutive layers. These two constraints can be expressed as follows:

$$w\_1^h \mathbf{x}\_1 + w\_2^h \mathbf{x}\_2 + \dots + w\_{i-1}^h \mathbf{x}\_{i-1} - \mathbf{x}\_i + w\_{i+1}^h \mathbf{x}\_{i+1} + \dots \cdot w\_m^h \mathbf{x}\_m = \mathbf{0} \tag{6}$$

$$\mathbf{x}\_i - w\_{i+m}^v \mathbf{x}\_{i+m} = \mathbf{0} \tag{7}$$

where *m* is the total number of voxels in the same layer. *w<sup>h</sup>* and *w<sup>v</sup>* denote the horizontal weighted coefficient and the vertical weighted coefficient, respectively. The horizontal weighted coefficient is constructed based on the Gaussian inverse distance weighted function as the following equation:

$$w\_{i-1}^h = -\frac{e^{-\frac{d\_{i,i-1}^2}{2\sigma^2}}}{\sum\_{j=1}^m e^{-\frac{d\_{j,i}^2}{2\sigma^2}}}\tag{8}$$

where *d* is the distance between the two voxels. *j* is a number from 1 to *m*, represented the voxel ordering of the same horizontal layer. *σ* denotes the smoothing factor. The vertical weighted coefficient is constructed based on exponential function as follows:

$$w\_{i+m}^{\upsilon} = e^{(h\_{i+m} - h\_i)/H\_s} \tag{9}$$

where *h* represents the height of the corresponding voxel and *Hs* refers to the water vapor scale height with an empirical value of 1.5 km [40]. Note that each voxel has corresponding equations for horizontal and vertical constraints.

Thus, the tomography model for water vapor reconstruction can be established by combining the observation equation of multi-GNSS and the two types of constraint equations.

$$
\begin{pmatrix} y\_{swv}^G \\ y\_{swv}^C \\ y\_{swv}^R \\ y\_{swv}^F \\ 0 \\ 0 \end{pmatrix} = \begin{pmatrix} A\_G \\ A\_C \\ A\_R \\ A\_E \\ H \\ V \end{pmatrix} \cdot x \tag{10}
$$

where *yswv* denotes the vector with SWV values derived from these four satellite systems; *A* represents the coefficient matrices of the observation equation for different types of satellite systems; *H* and *V* are the coefficient matrices of the horizontal and vertical constraints, respectively. The tomography solution of the unknown water vapor density vector *x* can be obtained as follows:

$$\begin{array}{l} \mathfrak{X} = \left( A\_G^T P\_G A\_G + A\_C^T P\_C A\_C + A\_R^T P\_R A\_R + A\_E^T P\_E A\_E + H^T P\_T H + V^T P\_V V \right)^{-1} \\ \cdot \left( A\_G^T P\_G y\_{\text{surv}}^G + A\_C^T P\_C y\_{\text{surv}}^C + A\_R^T P\_R y\_{\text{surv}}^R + A\_E^T P\_E y\_{\text{surv}}^E \right) \end{array} \tag{11}$$

where *P* represents the weighting matrices of different equations, which are determined by an optimal weighting method using the variance components estimation and homogeneity test [41]. Note that the number of satellite systems in Equations (8) and (9) can be adjusted in the experiment.

#### **3. Results**

#### *3.1. Experimental Description*

In this paper, the Hong Kong satellite Positioning Reference Station Network (SatRef) was selected to conduct the water vapor tomography experiment. We divided this research

region into 560 voxels ranging from 113.87◦ to 114.35◦, from 22.19◦ to 22.54◦, and from 0 to 8 km for longitude, latitude, and altitude, respectively; that is, a voxel of 8 × 7 in the horizontal direction and 10 layers in the vertical direction. As shown in Figure 1, thirteen GNSS sites (T430, HKKT, HKLT, HKSL, HKNP, HKMW, HKPC, HKLM, HKOH, HKSC, HKST, HKSS, HKWS) in this region were used to provide SWV in the tomography modeling, and one GNSS (HKQT) and one radiosonde (45004) site were selected to validate the results of the water vapor tomography.

**Figure 1.** Distribution of GNSS, radiosonde sites, and the horizontal structure of the voxels.

We utilized the GAMIT 10.71 software to estimate the tropospheric parameters including ZTD and gradient parameters using the four GNSS systems. In this process, the elevation cutoff angle was set to 15◦, the IGS precise ephemeris was adopted. Three MEGX stations (JFNG, URUM, and LHAZ) were incorporated into the solution model to reduce the strong correlation of tropospheric parameters caused by the short baseline. The processing strategies were set to LC\_AUTCLN and BASELINE modes, meaning that the ionosphere-free linear combination was selected and the orbital parameters were fixed, respectively. The tropospheric parameters, including troposphere delay gradients and ZTD at 4 and 2 h intervals, are estimated and interpolated to a 30 s sampling rate in the GAMIT software. After calculating the ZHD using the measured pressure recorded by an automatic meteorological device, the SWV values of each satellite system were obtained by using Equations (1) and (2).

In this experiment, the GNSS observation data of one month from DOY 121 to 151, 2021 were selected to conduct the modeling and solution of water vapor tomography. For each tomographic solution, the period covered is 0.5 h. To assess the performance of water vapor tomography based on different satellite systems and different combinations of satellite system, each tomographic solution has 15 results, including those of a signal satellite system, the combination of two satellite systems, the combination of three satellite systems, and the combination of four satellite systems, namely GPS (G), BDS (C), GLONASS (R), Galileo (E), GPS+BDS (GC), GPS+GLONASS (GR), GPS+Galileo (GE), BDS+GLONASS (CR), BDS+Galileo (CE), GLONASS+Galileo (RE), GPS+BDS+GLONASS (GCR), GPS+BDS+Galileo (GCE), GPS+GLONASS+Galileo (GRE), BDS+GLONASS+Galileo (CRE), and GPS+BDS+GLONASS+Galileo (GCRE). Note that both BDS-2 and BDS-3 were included in the experiment.
