A Tropospheric Tomography Method with a Novel Height Factor Model Including Two Parts: Isotropic and Anisotropic Height Factors

: Global Navigation Satellite System (GNSS) tomography has developed into an e ﬃ cient tool for sensing the high spatiotemporal variability of atmospheric water vapor. The integration of GNSS top signals and side rays for tropospheric tomography systems using a novel height factor model (HFM) is proposed and discussed in this paper. Within the HFM, the sectional slant wet delay (SWD) of inside signals (the part of the side signal inside the tomography area), which is considered a key factor for modeling side rays, is separated into isotropic and anisotropic components. Correspondingly, two height factors are deﬁned to calculate the isotropic and anisotropic part of tropospheric delays in the HFM. In addition, the dynamic tomography top boundary is ﬁrst analyzed and determined based on 30-year radiosonde data to reasonably divide signals into top and side rays. Four special experimental schemes based on di ﬀ erent tomography regions of Hong Kong are performed to assess the proposed HFM method, the results of which show increases of 33.42% in the mean utilization of rays, as well as decreases of 0.46 g / m 3 in the average root mean square error (RMSE), compared to the traditional approach, revealing the improvement of tomography solutions when the side signals are included in the modeling. Furthermore, compared with the existing correction model for modeling side rays, the water vapor proﬁles retrieved from the proposed improved model are closer to the radiosonde data, which highlights the advantages of the proposed HFM for optimizing the GNSS tomography model. gradient delay, is determined by formula deduction. The SWD of inside signals is calculated using the HFM from two components and added into the tomography system. Four experimental schemes are performed separately to test and assess the improved method. The water vapor proﬁles obtained by the traditional approach and proposed model are compared with those provided by radiosonde data, and the statistical results demonstrate the beneﬁt of including the side signals.


Introduction
Global Navigation Satellite System (GNSS) technology has become an efficient atmospheric water vapor detection tool for studying rainstorm events and the El Niño-Southern Oscillation (ENSO) phenomenon [1][2][3][4]. Highly accurate integrated water vapor (IWV) can be inverted through GNSS observation signals, which has been used for global meteorological assessment and analysis [5,6]. Furthermore, to obtain water vapor with a high spatial resolution, medical chromatography has been introduced to retrieve the three-dimensional water vapor fields over a study area [7]. Due to several significant advantages, including observation data accessibility, all-weather availability, wide-coverage, In addition, a further comparison of the proposed optimized model and the existing correction model is performed to illustrate the advantages of the proposed HFM in tomography modeling.

Height Factor Model for GNSS Tomography
The HFM, explicitly described in this section, provides the possibility to use side signals in the tomography model. The dynamic characteristic of the top boundary of the tomography area is studied for the first time in the HFM. Moreover, two height factors are defined to show the isotropic and anisotropic peculiarities of tropospheric delay, and the SWD of inside signals is estimated by calculating the two components.

Tropospheric Transmission of the GNSS Signal
Troposphere-induced path delays, having an evident effect on the analysis of GNSS, InSAR, Very Long Baseline Interferometry (VLBI), and other space sensors [26,28,29], are generally split into two parts: a dry component and wet component. Just as its name implies, the dry delay determined from surface pressure is caused by dry gases. On the contrary, water vapor is the primary cause of the wet delay, which can be estimated from the GNSS processing, for instance, GAMIT/GLOBK or Bernese software [16].
The wet delay along the GNSS ray path S, known as SWD, is derived from the following equation [30]: where s represents the signal path from the satellite to the receiver, and N w denotes the wet refractivity along the path. In general, SWD consists of isotropic and anisotropic parts [31]. The former, related to the ZWD of GNSS stations, is calculated with the wet mapping function m f w (ε) (e.g., VFM1 employed in this paper [32]), depending on the elevation angle ε. The latter, caused by the flow of atmospheric water vapor, can be estimated by the north-south wet gradient G w NS and the east-west wet gradient G w EW , as well as the horizontal gradient function m f g (ε) [33].

SWD = m f w (ε)·ZWD
isotropic + m f g (ε)· G w NS · cos(α) + G w EW · sin(α) where ZWD can be extracted from the zenith total delay (ZTD) estimated in the GNSS data processing [34]. α represents the azimuth angle of the GNSS signals. G w NS and G w EW , which represent north and east wet gradients, determine the variation of the delay with various azimuth angles. The SWV is practically proportional to the SWD, as follows [30]: where Π = 10 5 Tm +k 2 denotes the conversion factor. k 2 and k 3 are the constant factors with values of 16.48 K/hPa and 3.75 × 10 5 K 2 /hPa from [35], respectively, and R V = 461.53 J/(kg·K) represents the specific gas constants for water vapor. The weighted mean temperature T m is obtained using the following equation [30]: where T 0 is the surface temperature, which is available from the local ground meteorological station.

Dynamicity of the Tomography Top Boundary
The top boundary serves as a key factor for the GNSS water vapor tomography model. On the one hand, signals are classified into two categories (top and side signals) according to the positional Remote Sens. 2020, 12,1848 4 of 23 relationship between themselves and the top boundary. It is worth noting that a higher top boundary leads to a lower number of top rays in the traditional tomography model, according to Chen and Liu [36]. On the other hand, redundant voxels where the WVD is close to zero will be introduced into the tomography system with an increasing top boundary. As a consequence, it is essential to identify an appropriate tomography top boundary to utilize more signals and avoid unnecessary voxels.
Water vapor, having an approximately exponential decreasing trend with altitude, is considered to be a criterion for determining the top boundary. In Yao and Zhao [14], the principle that a certain height with a WVD of less than 0.2 g/m 3 can be considered as the top boundary of the tomography model is proposed.
In this study, the idea of determining the optimal tomography top boundary based on the water vapor vertical distribution derived from 30 years (1989-2018) of radiosonde data at Hong Kong King Park station (114.17 • E, 22.31 • N) is introduced. Accordingly, the water vapor profiles of different months in Hong Kong are shown in Figure 1. It can be seen that the distribution trend of water vapor in the atmosphere varies greatly in different seasons. Evidently, the optimal top boundary in summer and autumn is higher than that in spring and winter. The mean height values of the tomography top boundary for the four seasons (spring corresponding to March to May, summer corresponding to June to August, autumn corresponding to September to November, and winter corresponding to December to February) are approximately 10.2, 10.9, 10.6, and 9.6 km, respectively. As a consequence, the tomography top boundary changes in different seasons, even for the same region, which can be deduced from this research. What is more, the dynamic top boundary of the tomography area, particularly in regions with complex weather conditions, should be determined from the local water vapor profiles in different months.

Construction of the Height Factor Model
In the three-dimensional voxel tomography model, the criterion for judging the availability of GNSS signals is where they puncture. In other words, only GNSS signals passing from the top boundary are selected as effective rays for the tomography system, and a number of side signals penetrating from the side of the tomography area are eliminated. Evidently, this behavior not only reduces the utilization of GNSS signals, but also generates more empty voxels in the three-dimensional model, which leads to a decline in the accuracy and stability of the tomographic results. Consequently, the HFM introduced in this paper will show the process of side signal modeling. The key for modeling side rays is the estimation of the SWD of inside signals, which is discussed in the HFM.
For the integral rays penetrating from the top boundary, their SWDs are calculated by using Equation (2). For instance, the SWD of signal OM shown in Figure 2 is obtained by the ray path mapping of the ZWD of M , plus the gradient delay caused by the entire troposphere. Accordingly, the SWD of inside signals (e.g., signal OQ) can be estimated by the mapping of the ZWD of ray OQ and the gradient delay as a result of the partial tropospheric impact. Therefore, we present the estimation of the SWD of inside signals in terms of two aspects-isotropic and anisotropic components-which are carefully discussed in the following two subsections.
Generally, Equation (5) is discretized into the following form in GNSS tomography [26]: where (ℎ ) represents the wet refractivity of the ℎ layer and ℎ denotes the thickness of the ℎ layer. The ZWD from the station to any height ℎ above the site can be correctly estimated using Equation (6). Therefore, we can establish the HFM to analyze the relationship between the height ℎ Generally, Equation (5) is discretized into the following form in GNSS tomography [26]: where N w (h i ) represents the wet refractivity of the ith layer and h i denotes the thickness of the ith layer. The ZWD from the station to any height h above the site can be correctly estimated using Equation (6). Therefore, we can establish the HFM to analyze the relationship between the height h and the sectional ZWD, which may contribute to estimating the isotropic component of the SWD of inside signals. Radiosonde data, preserving a high accuracy and optimal vertical resolution of meteorological parameters along the tropospheric column, are used to construct the HFM in this paper. The detailed steps are given below.

1.
The initial radiosonde data is encrypted to produce uniform data with a vertical resolution of 100 m and the water vapor density of the interpolation points is calculated using the cubic spline interpolation algorithm [37]. In the processing of radiosonde observations, the vertical resolution of radiosonde data at different remote sensing time is variable, i.e., the resolution of data in the lower layer is higher than that in the upper layer. It should be noted that observations with an even distribution can generate a more accurate fitting model. Therefore, the initial radiosonde data should first be encrypted to produce uniform data.

2.
The partial ZWD is computed from the surface to a certain height using Equation (6).

3.
The isotropic height factor λ iso is defined, which represents the ratio between the partial content of ZWD and the total ZWD, as shown in Equation (7).
where ZWD h denotes the partial ZWD from the surface to the height h, and ZWD total derived from the GNSS data is the total ZWD of the stations. Figure 3 shows the relationship between the isotropic height factor and the height based on the 30-year radiosonde data from different months in Hong Kong.

4.
The function relationship between the λ iso (h) and h is analyzed. In view of the variety of atmospheric water vapor, owing to the fact that the radiosonde balloon in Hong Kong is launched twice a day, at 00:00 UTC and 12:00 UTC, the radiosonde data cannot provide sufficient observations for daily-scale analysis. In addition, according to the periodic characteristics of ZTD and precipitable water vapor (PWV) [1,38], the relationship between the isotropic height factor and the height using 30 years of radiosonde data from 12 months is presented in Figure 3. Besides, the long-term trends of the isotropic height factor will not have an important impact on the tomography results. On the other hand, the wet refractivity water vapor decreases exponentially with height in the troposphere, and it is found from Equation (7) that the isotropic height factor λ iso is related to the wet refractivity in the zenith direction. In the work of Yao and Zhao [14], Zhao et al. [39], and Zhao et al. [27], a similar scale factor or truncation factor λ is introduced, the trend of which changes with the height, as analyzed based on the exponential relationship in all three studies.

5.
Accordingly, the HFM is established as follows: where a 1 , b 1 , a 2 , and b 2 represent the coefficients of the HFM, which are determined by least-squares fitting [27]. It should be noted that the different weather conditions during the 30 years are not classified, which may have an impact on obtaining the SWV of side signals and needs to be analyzed in future work.
where , , , and represent the coefficients of the HFM, which are determined by leastsquares fitting [27]. It should be noted that the different weather conditions during the 30 years are not classified, which may have an impact on obtaining the SWV of side signals and needs to be analyzed in future work.  Table 1 shows the optimal coefficient for each month. Furthermore, the conclusion that the root mean square error (RMSE) of the 12 months of fitting results is less than 0.05 and the R-square representing the quality of the fitting results is greater than 0.98 can be inferred from   Table 1 shows the optimal coefficient for each month. Furthermore, the conclusion that the root mean square error (RMSE) of the 12 months of fitting results is less than 0.05 and the R-square representing the quality of the fitting results is greater than 0.98 can be inferred from Figure 4. It should be noted that the closer R-square gets to 1, the better the fitting results are [14]. Consequently, the isotropic component of SWD for inside signals could be accurately calculated using the HFM, as shown in Equation (9).

Anisotropic Component of SWD for Inside Signals
Gradient delays, characterizing the anisotropic water vapor distribution, have been intensively studied by a large number of scholars in VLBI analysis [28,32,33]. In the deduction of Boehm and Schuh [32], the two gradients mentioned above are obtained from the vertical integral of the water vapor refractivity gradient with height: where denotes the azimuth direction (e.g., or ), (ℎ) is the water vapor refractivity gradient at height ℎ, and represents the top of the troposphere in the integration.
Besides, the prevailing assumption that the refractivity gradient decreases exponentially with height was proposed in the work of Chen and Herring [33].
where is the refractivity gradient at the surface and denotes the scale height of the neutral atmosphere, with a value of approximately 2 km, according to the local ground meteorological data [40]. The integral of Equation (10) can be written with the following expression:

Anisotropic Component of SWD for Inside Signals
Gradient delays, characterizing the anisotropic water vapor distribution, have been intensively studied by a large number of scholars in VLBI analysis [28,32,33]. In the deduction of Boehm and Schuh [32], the two gradients mentioned above are obtained from the vertical integral of the water vapor refractivity gradient with height: where α denotes the azimuth direction (e.g., NS or EW), N α (h) is the water vapor refractivity gradient at height h, and H top represents the top of the troposphere in the integration. Besides, the prevailing assumption that the refractivity gradient decreases exponentially with height was proposed in the work of Chen and Herring [33].
where N S is the refractivity gradient at the surface and H scale denotes the scale height of the neutral atmosphere, with a value of approximately 2 km, according to the local ground meteorological data [40]. The integral of Equation (10) can be written with the following expression: In this section, an analysis of the gradient delay between the integral signal (e.g., OP in Figure 2) and inside signal (e.g., OQ in Figure 2) is performed. The anisotropic gradient delay of signal OP and OQ can be calculated with Equation (2).
where ∆L OP G and ∆L OQ G represent the anisotropic component of SWD of integral signal OP and inside signal OQ, respectively, and m f g (ε) = 1 sin(ε)· tan(ε)+0.003 denotes the horizontal gradient mapping function [33]. G OP α and G OQ α , the wet gradients of two rays, can be computed by using Equation (12).
Remote Sens. 2020, 12, 1848 9 of 23 Therefore, we can define the anisotropic height factor λ aniso to estimate the relationship between the ∆L OP G and ∆L OQ G .
where h OQ and H top represent the height of the zenith projection of the slant path OQ and OP along the zenith, respectively. During the tomography periods in Hong Kong in August 2017, the tomography top boundary H top was determined to be 11 km from Figure 1, and H scale was set to 2 km, according to the ground meteorological data in Hong Kong. Therefore, the anisotropic height factor λ aniso h OQ only has a functional relationship with height h OQ . For instance, h OQ increases equidistantly from 0 to 11 km with the interval of 0.1 km, and the corresponding λ iso h OQ could be calculated by Equation (15).
With the sample data of h OQ and λ iso h OQ , the red curve in Figure 5 could be drawn.
inside signal (e.g., in Figure 2) is performed. The anisotropic gradient delay of signal and can be calculated with Equation (2).
where ∆ and ∆ represent the anisotropic component of SWD of integral signal and inside signal , respectively, and ( ) = ( )• ( ) . denotes the horizontal gradient mapping function [33]. and , the wet gradients of two rays, can be computed by using Equation (12). Therefore, we can define the anisotropic height factor to estimate the relationship between the ∆ and ∆ .
where ℎ and represent the height of the zenith projection of the slant path and along the zenith, respectively. During the tomography periods in Hong Kong in August 2017, the tomography top boundary was determined to be 11 km from Figure 1, and was set to 2 km, according to the ground meteorological data in Hong Kong. Therefore, the anisotropic height factor ℎ only has a functional relationship with height ℎ . For instance, ℎ increases equidistantly from 0 to 11 km with the interval of 0.1 km, and the corresponding ℎ could be calculated by Equation (15). With the sample data of ℎ and ℎ , the red curve in Figure 5 could be drawn. Similarly, the anisotropic gradient delay of inside signals can be written in the following equation: where is the anisotropic delay component of inside signals on the same ray path with the observed signal, and ∆ represents the anisotropic delay component of the observed signal.
As a consequence, the SWD of inside signals can be computed with Similarly, the anisotropic gradient delay of inside signals can be written in the following equation: where SWD anisotropic is the anisotropic delay component of inside signals on the same ray path with the observed signal, and ∆L integral G represents the anisotropic delay component of the observed signal. As a consequence, the SWD of inside signals can be computed with In this work, ZWD total and ∆L integral G (ε, α). were calculated with the ZTD, G NS , and G EW , respectively, which were estimated by GAMIT/GBLOK(v.10.7) [41]. The SWD of inside signals could be converted into SWV using Equation (3) for the improved tomography model.

Modeling the GNSS Observations with the Proposed HFM
In common GNSS tomographic models, only GNSS signals passed from the top boundary are considered as effective observations for the tomographic model [10]. In this section, an improved GNSS tomography algorithm based on HFM is introduced, as shown in Figure 6. The process of GNSS tomography equations using integral signals and inside signals is discussed, respectively. The horizontal and vertical constraints are also described in detail.
horizontal and vertical constraints are also described in detail.
Additionally, in most tomography models, to obtain as many effective signals as possible, the influence of the elevation angle is neglected. However, signals which enter the neutral atmosphere at low elevation angles (e < 15°) are significantly affected by atmospheric bending, according to Möller and Landskron [38]. As the bending effect is not modeled by the tomographic algorithms, rays with an elevation angle below 15° are first eliminated.

Constructing Tomography Observation Equations Using Integral Signals
To reconstruct the three-dimensional water vapor field, the tomographic area is divided into a number of uniformly distributed voxels over a GNSS network. Generally, the WVD within a voxel, assumed to be constant for a tomographic epoch, is estimated as one parameter. GNSS signals, from the GNSS sites to the top or side boundary of the tomographic region, provide the distance traveled in each voxel. Accordingly, Equation (1) can be discretized to a linear equation with the following form: Additionally, in most tomography models, to obtain as many effective signals as possible, the influence of the elevation angle is neglected. However, signals which enter the neutral atmosphere at low elevation angles (e < 15) are significantly affected by atmospheric bending, according to Möller and Landskron [38]. As the bending effect is not modeled by the tomographic algorithms, rays with an elevation angle below 15 • are first eliminated.

Constructing Tomography Observation Equations Using Integral Signals
To reconstruct the three-dimensional water vapor field, the tomographic area is divided into a number of uniformly distributed voxels over a GNSS network. Generally, the WVD within a voxel, assumed to be constant for a tomographic epoch, is estimated as one parameter. GNSS signals, from the GNSS sites to the top or side boundary of the tomographic region, provide the distance traveled in each voxel. Accordingly, Equation (1) can be discretized to a linear equation with the following form: Note that the SWD of GNSS rays was converted into SWV using Equation (3) for the tomographic system. d ij denotes the distance traveled by the jth ray in the ith voxel, whose WVD is represented by x i . n is the total number of voxels intersected by the ray j.
For integral signals, passing from the top boundary of the tomography area, the observation equation can be written in matrix notation: where SWV integral is the observation vector of integral signals and A integral . represents the observation matrix, including the distance traveled in each voxel. X denotes the unknown vector for the three-dimensional tomographic model.

Constructing Tomography Observation Equations Using Inside Signals
For inside signals, penetrating from the side of the tomography area, the SWD should first be estimated by HFM. The specific steps of computing the SWD inside are as follows: 1.
The ZTD and horizontal gradients of effective side signals (e.g., OP in Figure 2) are estimated using GAMIT/BLOCK. The ZWD, north, and east wet gradients can be extracted from the estimated values [7]; 2.
The height h inter of the intersection point between the effective side signals and the side face of the tomographic area is calculated. The λ iso (h inter ) and λ aniso (h inter ) are estimated by Equations (8) and (15), respectively. Then, the SWD of inside signals can be calculated using Equation (17); 3.
The distance information of inside signals is obtained with the ray-tracing method [42], and the unknowns are related to the inside observations by the following equation: where SWV inside is the observation vector of inside signals and A inside represents the observation matrix derived from inside signals.

Constructing Tomography Constraint Equations
Although side signals with an elevation angle greater than 15 • are utilized, there are still some empty voxels in the three-dimensional tomographic model [10]. Usually, spatial constraints are introduced in the tomography system to relieve the ill-conditioned problem.
The following horizontal constraint, assuming that the WVD of a voxel is equal to the weighted mean of the WVD of the surrounding voxels in the same layer [7], is applied for the tomography system: where n h represents the number of voxels in a certain layer and x i is the WVD of the ith voxel. w ij denotes the horizontal weighted coefficient of the jth surrounding voxel in this layer corresponding to the ith voxel. The Gauss weighting function proposed in Song et al. [21] is introduced to calculate the coefficients: where d ij is the distance between the ith voxel and the jth surrounding voxel. σ represents the smoothing factor and is 1.5 times longer than the horizontal length of voxels [21]. In addition, the WVD can be described by exponential decay with an increasing height [43]: where x k and h k denote the WVD and height of the kth layer, and H scale is the scale height mentioned in Equation (11). As a result, the final tomographic system equations can be written as follows: where Y integral and Y inside are the observation vector of the integral signals and side signals, respectively. Similarly, both A integral and A inside represent the observation matrix, including the distance traveled in each voxel. The weighting matrix was also imposed on the tomographic observation equations. P integral and P inside , representing the weighting matrix of integral and inside rays, have a functional relationship with the elevation angle of the signal path [26]. However, for the weighting matrix of horizontal and vertical constraints, P hz and P vert , having a similar impact on voxels, are by default the identity matrix [26].
Due to the morbidity of the coefficient matrix, the Algebraic Reconstruction Technique (ART), with the advantages of simple iteration and fast convergence [10,44], is introduced to solve the tomography observation equations in this work.

Experiments and Results
In general, since voxels located on the side of the tomographic area do not travel by any ray, the water vapor density of these voxels can be estimated based on only constraint equations, which reduces the accuracy of the tomography solutions. The primary contribution of the side signals is to ensure that voxels not passed in the traditional model can be crossed, which allows the accuracy and stability of tomographic results to be improved. To test the HFM and evaluate the efficiency of side signals, four special experimental schemes were performed separately. Section 4.1 describes the implementation of the four experimental schemes; the contribution of the side signals to the tomography solutions is analyzed in Section 4.2. On the other hand, the existing improved model combining the side observations was reestablished using the GNSS data [14], and is compared to the proposed HFM in Section 4.3.
It should be noted that when assessing the accuracy of tomography solutions, the vertical profiles obtained from radiosonde were used as reference values. Although the radiosonde data has been made use of to calculate the isotropic height factor, its main contribution is to provide the vertical distribution characteristics of atmospheric water vapor, not to participate in the tomographic processing. For instance, the tomographic time was set to August 2017, according to the HFM, and the optimal tomography top boundary was set to 11 km. Additionally, the height factor model λ iso (h) = 1.121·e −0.051h − 0.191·e −1.019h was applied to calculate the isotropic height factor. Consequently, only GNSS observations went sent to the tomography system as input information for reconstructing the three-dimensional water vapor field, which shows that the proposed HFM does not interfere with the accuracy assessment of tomography results.

Experimental Schemes
The GNSS data with the sampling rate of 30 s from 17 stations (red dot in Figure 7) provided by the Hong Kong Satellite Positioning Reference Station Network (SatRef) were processed using GAMIT/GLOBK (v.10.7). In the processing, the four International GNSS Service (IGS) stations (BJFS station, HKWS station, LHAZ station, and SHAO station) were introduced. The ZTD and the gradient delay were estimated with the time resolution of 5 min, which was sufficient to reveal the variation characteristics of water vapor. The Saastamoinen model was set as the expression for the ZHD [45], and Vienna Mapping Function 1 (VMF1) was used for calculating the SWDs [46]. In addition, the post-fit residuals calculated by GAMIT/GLOBK exceeding 2.5 times the standard deviation were removed [47]. The 30-year radiosonde observation data was obtained from the King's Park Meteorological Station (HKKP, the blue dot in Figure 7a). The research area of the tomography model spanned from 113.82 • E to 114.36 • E and 22.16 • N to 22.56 • N. Moreover, the horizontal resolution was 0.09 • in longitude and 0.08 • in latitude, as shown in Figure 7a. In the current tomography studies based on the Hong Kong SatRef, about a month of GNSS data and the radiosonde observations derived from the HKKP station were processed for modeling the tomography area [14,20,36,41]. Therefore, the experimentation time in this work was set as August 2017, i.e., day of year (DOY) 213-243, when Hong Kong was in the summer and there were more rainstorms. The tomography region was segmented into 15 nonuniform layers from 0 to 11 km (August top boundary determined in Figure 2) [14]. 0.09° in longitude and 0.08° in latitude, as shown in Figure 7a. In the current tomography studies based on the Hong Kong SatRef, about a month of GNSS data and the radiosonde observations derived from the HKKP station were processed for modeling the tomography area [14,20,36,41]. Therefore, the experimentation time in this work was set as August 2017, i.e., day of year (DOY) 213-243, when Hong Kong was in the summer and there were more rainstorms. The tomography region was segmented into 15 nonuniform layers from 0 to 11 km (August top boundary determined in Figure 2) [14]. In the most common voxel-based tomography model, to assess the accuracy of tomographic results, water vapor profiles of the voxel columns (small blue rectangles in Figure 7b) where the radiosonde is located are compared with radiosonde-derived vertical profiles. However, the empty voxels are mainly distributed on the side of the three-dimensional tomographic model. Consequently, as shown in Figure 8, four special test schemes based on four different tomography areas (red boxes of Figure 8a to 8d) were implemented. The RS-voxels (the blue voxels, shown in Figure 7b, where the radiosonde is located) were placed on four boundaries of the tomographic region in the four schemes, for instance, RS-voxels on the east side (Figure 8a), RS-voxels on the west side (Figure 8b), RS-voxels on the south side (Figure 8c), and RS-voxels on the north side (Figure 8d). Table 2 shows the accurate range of the tomographic area and the number of voxels in these four schemes. In the most common voxel-based tomography model, to assess the accuracy of tomographic results, water vapor profiles of the voxel columns (small blue rectangles in Figure 7b) where the radiosonde is located are compared with radiosonde-derived vertical profiles. However, the empty voxels are mainly distributed on the side of the three-dimensional tomographic model. Consequently, as shown in Figure 8, four special test schemes based on four different tomography areas (red boxes of Figure 8a-d) were implemented. The RS-voxels (the blue voxels, shown in Figure 7b, where the radiosonde is located) were placed on four boundaries of the tomographic region in the four schemes, for instance, RS-voxels on the east side (Figure 8a), RS-voxels on the west side (Figure 8b), RS-voxels on the south side (Figure 8c), and RS-voxels on the north side (Figure 8d). Table 2 shows the accurate range of the tomographic area and the number of voxels in these four schemes.    In addition, to verify the accuracy and reliability of the HFM, two sub-groups of solutions were generated for the four schemes; one only used the rays travelling from the top boundary, while the others combined both rays passing from the top and side boundaries. It should be noted that if the effective observations from the HKSC station are included, the data is not consistent with the signal distribution of the general tomographic model. In addition, due to the RS-voxels always being crossed by rays in the four schemes, the benefit of the side signal for improving the accuracy of the water vapor density of the un-punctured voxels cannot be evaluated accurately. Therefore, all observations from HKSC station were excluded in the four schemes.

Contribution Analysis of the Side Signals
The isotropic and anisotropic SWV values of inside signals were first estimated by the HFM, as listed in Figure 9. A slight consistency between the total SWV and the isotropic one can be noticed from the four sub-figures. In contrast, the anisotropic parts were inconsistent with the totals and only accounted for an average of 2.31%, 2.56%, 2.61%, and 2.25% of the corresponding total component in the four schemes, which reveals that the isotropic portion has a more marked significance than the anisotropic portion in estimating the SWV value of the inside signals. However, the latter, with a value of approximately 20 mm for certain rays, is also important. Furthermore, we developed a detailed analysis of the number of rays crossing from the RS-voxel columns. Figure 10 demonstrates the benefit of absorbing inside signals into the traditional tomography model. Evidently, in the East, West, and North schemes, no signals passed from the voxel from layer 1 to layer 8 in the general tomography model, while only the first four layers of voxels were not penetrated by rays in the improved model. As far as the South scheme is concerned, it is noticeable that the RS-voxel columns were all passed by GNSS rays in the optimized model, except for the voxel on the first layer. The number of penetrated signals tripled for voxels from 2 to Furthermore, we developed a detailed analysis of the number of rays crossing from the RS-voxel columns. Figure 10 demonstrates the benefit of absorbing inside signals into the traditional tomography model. Evidently, in the East, West, and North schemes, no signals passed from the voxel from layer 1 to layer 8 in the general tomography model, while only the first four layers of voxels were not penetrated by rays in the improved model. As far as the South scheme is concerned, it is noticeable that the RS-voxel columns were all passed by GNSS rays in the optimized model, except for the voxel on the first layer. The number of penetrated signals tripled for voxels from 2 to 11 layers when adding the inside rays.  The three-dimensional WVD was reconstructed using the traditional method and the proposed approach, respectively. To assess the retrieval capabilities of the optimized tomographic model, the SWV differences between the observed SWV and reconstructed SWV from the three-dimensional tomographic water vapor field were calculated, and the SWV differences between the traditional method and the proposed approach were obtained, respectively. Figure 12 illustrates a comparison of the SWV differences between the two methods in the four schemes. The average number of effective rays and the mean signal utilization rate for the traditional model and optimized model during the tomographic period are also compared in Figure 11, where the histogram shows the number of signals and the line chart represents the signal utilization rate. The former refers to the 31-day average number of GNSS rays selected in the tomography model during the two epochs. Similarly, the latter represents the 31-day mean utilization rate of the number of rays used in the equations out of the total number of GNSS signals in each tomography epoch. Because the GNSS rays with an elevation angle below 15 • were eliminated, the mean utilization rate of integral signals and inside ones was about 80% in the four schemes. With side signals absorbed into the tomography system, the mean effective rays increased by 32.33%, whereas the average utilization rate of GNSS signals was enhanced by 33.42%, from 47.12% to 80.54%. In particular, the number of effective signals combining inside ones remained an approximative constant for any tomographic scheme, which enhanced the stability of the improved tomography model.  The three-dimensional WVD was reconstructed using the traditional method and the proposed approach, respectively. To assess the retrieval capabilities of the optimized tomographic model, the SWV differences between the observed SWV and reconstructed SWV from the three-dimensional tomographic water vapor field were calculated, and the SWV differences between the traditional method and the proposed approach were obtained, respectively. Figure 12 illustrates a comparison of the SWV differences between the two methods in the four schemes. The three-dimensional WVD was reconstructed using the traditional method and the proposed approach, respectively. To assess the retrieval capabilities of the optimized tomographic model, the SWV differences between the observed SWV and reconstructed SWV from the three-dimensional tomographic water vapor field were calculated, and the SWV differences between the traditional method and the proposed approach were obtained, respectively. Figure 12 illustrates a comparison of the SWV differences between the two methods in the four schemes. There is a noticeably similar tendency, in that the SWV differences of GNSS signals with an elevation angle larger than 45° are in the range of -5 to 5 mm, while those of rays with an elevation angle lower than 45° sharply increase with a decreasing elevation angle in all four schemes. As far as the comparison of SWV differences between the general and optimal tomography model is concerned, a remarkable improvement can be identified in the East scheme, West scheme, and South scheme, whereas the proposed method is slightly better than the traditional one in the North scheme. The same conclusion can also be deduced from Table 3, where the statistics, including the RMSE, standard deviation (STD), and bias, for the SWV differences from the four schemes are compared. The improvements of the combined observations with respect to only the integral one are 3.11/3.34, 2.55/2.19, 2.60/2.62, and 1.27/1.35 mm for the RMSE/STD in the four schemes, respectively. In the bias comparison, there is a visible advancement of the tomography solutions when inside signals are included.  There is a noticeably similar tendency, in that the SWV differences of GNSS signals with an elevation angle larger than 45 • are in the range of −5 to 5 mm, while those of rays with an elevation angle lower than 45 • sharply increase with a decreasing elevation angle in all four schemes. As far as the comparison of SWV differences between the general and optimal tomography model is concerned, a remarkable improvement can be identified in the East scheme, West scheme, and South scheme, whereas the proposed method is slightly better than the traditional one in the North scheme. The same conclusion can also be deduced from Table 3, where the statistics, including the RMSE, standard deviation (STD), and bias, for the SWV differences from the four schemes are compared. The improvements of the combined observations with respect to only the integral one are 3.11/3.34, 2.55/2.19, 2.60/2.62, and 1.27/1.35 mm for the RMSE/STD in the four schemes, respectively. In the bias comparison, there is a visible advancement of the tomography solutions when inside signals are included.
In addition, the tomographic results derived from the general model and improved model were compared with radiosonde data at 00:00 UTC and 12:00 UTC, when the radiosonde data provided accurate water vapor profiles at different altitudes. Table 4 lists the mean RMSE and STD of the tomography results retrieved from the integral signals and combined signals for August 2017. It can be observed that the tomography results are closer to the radiosonde data when the inside rays are absorbed. The improved tomography model displays a smaller RMSE than the conventional one in the four test schemes, with decreases of 33.12%, 30.86%, 19.27%, and 30.67% in the mean RMSE, and the approximate consistent decrease rate of the STD for the four schemes can also be computed from Table 4.

Comparison with the Existing Correction Model
In this section, the HFM in the current study is compared with the existing correction model in Yao and Zhao [14], which is named the Yao model in this work for convenience. In the Yao model, a similar exponential relationship between the scale factor and the height is established, as follows [14]: where α is the scale factor, and a and b represent the coefficients of the scale factor, respectively. h denotes the height of the intersection point between the effective side signals and the side face of the tomography area.
To compare the performance of the HFM and Yao model in retrieving the three-dimensional water vapor field, three further schemes were designed for tomography modeling based on all GNSS station observations from the Hong Kong SatRef, in addition to the HKSC station. The tomographic region for the three schemes covers from 113.82 • E to 114.36 • E in longitude and 22.16 • N to 22.56 • N in latitude. Other tomographic model configurations, such as a horizontal and vertical resolution, are consistent with the four special experimental schemes. The three schemes in this section are as follows.
• Scheme I: The traditional tomography method that only considers the GNSS signals crossing the top boundary is adopted to construct the observation equations; • Scheme II: The Yao model is used to estimate the SWV of side rays and build the equations with both rays passing from the top and side boundary; • Scheme III: The HFM is employed to calculate the SWV side and structure the tomography system with top and side observations.
Based on GNSS data at two tomography epochs of 00:00-00:30 UTC and 12:00-12:30 UTC during the period of DOY 213-243, 2017, the scale factor in the Yao model was calculated and is illustrated in Figure 13, where a similar exponential relationship can be observed. It is noticeable that owing to the neglection of GNSS signals with an elevation angle of less than 15 • , all of the heights of the intersection points are higher than approximately 1 km. Accordingly, the coefficients a and b were estimated by the least-squares method, and the results are given in Table 5 [27].  0.0145 With the scale factor in the Yao model and the HFM proposed in this paper, the SWV of the side signals could be estimated based on the height of the intersection points. Consequently, Figure 14 presents a comparison of the SWV value of effective side rays between the Yao model and HFM at two tomography epochs during each of the test days. It is clear that the HFM can estimate slightly higher SWV values than the Yao model for most signals, and the latter provides a rather small SWV estimation, with a value of approximately 20 mm, for several observation rays. To assess the water vapor reconstruction capabilities of different methods, the tomographic results derived from the three schemes, described in Figure 15, were employed to illustrate the water vapor vertical profile at 00:00 UTC and 12:00 UTC on DOY 220 (sunny day) and DOY 234 (rainy day).  With the scale factor in the Yao model and the HFM proposed in this paper, the SWV of the side signals could be estimated based on the height of the intersection points. Consequently, Figure 14 presents a comparison of the SWV value of effective side rays between the Yao model and HFM at two tomography epochs during each of the test days. It is clear that the HFM can estimate slightly higher SWV values than the Yao model for most signals, and the latter provides a rather small SWV estimation, with a value of approximately 20 mm, for several observation rays.  0.0145 With the scale factor in the Yao model and the HFM proposed in this paper, the SWV of the side signals could be estimated based on the height of the intersection points. Consequently, Figure 14 presents a comparison of the SWV value of effective side rays between the Yao model and HFM at two tomography epochs during each of the test days. It is clear that the HFM can estimate slightly higher SWV values than the Yao model for most signals, and the latter provides a rather small SWV estimation, with a value of approximately 20 mm, for several observation rays. To assess the water vapor reconstruction capabilities of different methods, the tomographic results derived from the three schemes, described in Figure 15, were employed to illustrate the water vapor vertical profile at 00:00 UTC and 12:00 UTC on DOY 220 (sunny day) and DOY 234 (rainy day). To assess the water vapor reconstruction capabilities of different methods, the tomographic results derived from the three schemes, described in Figure 15, were employed to illustrate the water vapor vertical profile at 00:00 UTC and 12:00 UTC on DOY 220 (sunny day) and DOY 234 (rainy day). It can be observed that the water vapor profiles derived from the three considered schemes agree with the reference profile obtained by the radiosonde data. Furthermore, a higher coincidence of the vertical profile can be seen on majority layers when combining the inside signals. As far as the quantitative comparison of WVD is concerned, the accuracy of the tomographic solutions from Scheme III (mean RMSE are 0.69 and 0.95 g/m 3 on both days) is superior to those from Scheme I (mean RMSE are 1.40 and 1.30 g/m 3 on both days) and Scheme II (mean RMSE are 0.73 and 1.20 g/m 3 on both days), which reveals that the proposed model has better reconstruction capabilities than the Yao model and traditional model in both weather conditions.
In addition, to further show the superiority of the proposed method in the WVD comparison, the RMSE of tomographic results from the three schemes was calculated, and the results are shown in Figure 16, including the 31-day period from DOY 213 to DOY 243, 2017, at 00:00 and 12:00 UTC daily. It is clear that, among the three schemes, Scheme III has the smallest RMSE compared to Scheme I and Scheme II in the majority of time periods. The average improvement of Scheme III with respect to Scheme I and Scheme II was 0.51 and 0.20 g/m 3 in terms of the RMSE, respectively. As a further assessment, Table 2 lists the statistics of the three tomography solutions, including the maximum, minimum, and mean values of RMSE, STD, and bias.   It can be observed that the water vapor profiles derived from the three considered schemes agree with the reference profile obtained by the radiosonde data. Furthermore, a higher coincidence of the vertical profile can be seen on majority layers when combining the inside signals. As far as the quantitative comparison of WVD is concerned, the accuracy of the tomographic solutions from Scheme III (mean RMSE are 0.69 and 0.95 g/m 3 on both days) is superior to those from Scheme I (mean RMSE are 1.40 and 1.30 g/m 3 on both days) and Scheme II (mean RMSE are 0.73 and 1.20 g/m 3 on both days), which reveals that the proposed model has better reconstruction capabilities than the Yao model and traditional model in both weather conditions.

RMSE STD Bias
In addition, to further show the superiority of the proposed method in the WVD comparison, the RMSE of tomographic results from the three schemes was calculated, and the results are shown in Figure 16, including the 31-day period from DOY 213 to DOY 243, 2017, at 00:00 and 12:00 UTC daily. It is clear that, among the three schemes, Scheme III has the smallest RMSE compared to Scheme I and Scheme II in the majority of time periods. The average improvement of Scheme III with respect to Scheme I and Scheme II was 0.51 and 0.20 g/m 3 in terms of the RMSE, respectively. As a further assessment, Table 2 lists the statistics of the three tomography solutions, including the maximum, minimum, and mean values of RMSE, STD, and bias. It can be observed that the water vapor profiles derived from the three considered schemes agree with the reference profile obtained by the radiosonde data. Furthermore, a higher coincidence of the vertical profile can be seen on majority layers when combining the inside signals. As far as the quantitative comparison of WVD is concerned, the accuracy of the tomographic solutions from Scheme III (mean RMSE are 0.69 and 0.95 g/m 3 on both days) is superior to those from Scheme I (mean RMSE are 1.40 and 1.30 g/m 3 on both days) and Scheme II (mean RMSE are 0.73 and 1.20 g/m 3 on both days), which reveals that the proposed model has better reconstruction capabilities than the Yao model and traditional model in both weather conditions.
In addition, to further show the superiority of the proposed method in the WVD comparison, the RMSE of tomographic results from the three schemes was calculated, and the results are shown in Figure 16, including the 31-day period from DOY 213 to DOY 243, 2017, at 00:00 and 12:00 UTC daily. It is clear that, among the three schemes, Scheme III has the smallest RMSE compared to Scheme I and Scheme II in the majority of time periods. The average improvement of Scheme III with respect to Scheme I and Scheme II was 0.51 and 0.20 g/m 3 in terms of the RMSE, respectively. As a further assessment, Table 2 lists the statistics of the three tomography solutions, including the maximum, minimum, and mean values of RMSE, STD, and bias.   It should be noted from Table 6 that all of the means of the three statistics in Scheme II and III are superior to those of Scheme I, which suggests that there is an improvement of the GNSS tomography solution when side observations are included in the tomography modeling. As far as Scheme II and III are concerned, the average RMSE and STD for the former are 1.28 and 1.25 g/m 3 and the latter are 1.08 and 1.07 g/m 3 , respectiely, and the mean accuracy improvement rates are 15.62% and 14.40% for both statistics, respectively, which demonstrates that the proposed HFM provides more accurate side SWV observations and three-dimensional water vapor fields. Following this, to quantitatively compare the accuracy of each layer's tomographic results, the reference value of WVD, along the vertical columns, was interpolated using the radiosonde data, and the mean difference between the tomographic results obtained from the three schemes and the radiosonde data during the tomographic period was calculated and is reported in Figure 17. It is evident that the absolute value of the average differences dramatically decreases with increasing height layers. Besides, it should be noted that in the earth surface layer from 0 to 2 km, the difference of most voxels is less than 0 in the three schemes, i.e., the WVD retrieved from the GNSS tomography solution is lower than that obtained from the radiosonde. However, it is noticeable that there is an impressive improvement of the WVD in these layers in Scheme II and Scheme III, although the voxels in this range are not penetrated by side signals according to Figure 10. This can be explained by the fact that when combining both integral signals and inside ones for GNSS tomography, the number of effective observations is obviously increased and the spatial geometric defect of the three-dimensional tomographic model and the ill condition of the tomographic observation equations are gradually remedied, which is beneficial for improving the accuracy and stability of the tomography results.  Table 6 that all of the means of the three statistics in Scheme II and III are superior to those of Scheme I, which suggests that there is an improvement of the GNSS tomography solution when side observations are included in the tomography modeling. As far as Scheme II and III are concerned, the average RMSE and STD for the former are 1.28 and 1.25 g/m 3 and the latter are 1.08 and 1.07 g/m 3 , respectiely, and the mean accuracy improvement rates are 15.62% and 14.40% for both statistics, respectively, which demonstrates that the proposed HFM provides more accurate side SWV observations and three-dimensional water vapor fields.

RMSE STD Bias
Following this, to quantitatively compare the accuracy of each layer's tomographic results, the reference value of WVD, along the vertical columns, was interpolated using the radiosonde data, and the mean difference between the tomographic results obtained from the three schemes and the radiosonde data during the tomographic period was calculated and is reported in Figure 17. It is evident that the absolute value of the average differences dramatically decreases with increasing height layers. Besides, it should be noted that in the earth surface layer from 0 to 2 km, the difference of most voxels is less than 0 in the three schemes, i.e., the WVD retrieved from the GNSS tomography solution is lower than that obtained from the radiosonde. However, it is noticeable that there is an impressive improvement of the WVD in these layers in Scheme II and Scheme III, although the voxels in this range are not penetrated by side signals according to Figure 10. This can be explained by the fact that when combining both integral signals and inside ones for GNSS tomography, the number of effective observations is obviously increased and the spatial geometric defect of the threedimensional tomographic model and the ill condition of the tomographic observation equations are gradually remedied, which is beneficial for improving the accuracy and stability of the tomography results.
Furthermore, the conclusion that the differences of Scheme III are smaller than those of Scheme II in each layer can be drawn from Figure 17, which highlights the advantages of the proposed HFM. This could be because of the higher side SWV observations estimated from the HFM in Figure 14, which results in the tomography solutions being closer to the radiosonde measurements.

Conclusions
A novel tropospheric tomography methodology based on the HFM, combining integral and side signals, has been introduced and demonstrated in this paper. The SWD of GNSS rays, which is a crucial observation for the GNSS tomography system, has been separated into two parts-isotropic and anisotropic components-corresponding to the isotropic and anisotropic height factors in the HFM. With these two factors, the SWD of the side rays can be estimated accurately. It should be noted that the signals with an elevation angle of less than 15° were eliminated due to bending of the ray Furthermore, the conclusion that the differences of Scheme III are smaller than those of Scheme II in each layer can be drawn from Figure 17, which highlights the advantages of the proposed HFM. This could be because of the higher side SWV observations estimated from the HFM in Figure 14, which results in the tomography solutions being closer to the radiosonde measurements.

Conclusions
A novel tropospheric tomography methodology based on the HFM, combining integral and side signals, has been introduced and demonstrated in this paper. The SWD of GNSS rays, which is a crucial observation for the GNSS tomography system, has been separated into two parts-isotropic and anisotropic components-corresponding to the isotropic and anisotropic height factors in the HFM. With these two factors, the SWD of the side rays can be estimated accurately. It should be noted that the signals with an elevation angle of less than 15 • were eliminated due to bending of the ray path. In addition, the dynamic tomography top boundary was determined based on 30-year radiosonde data, which preserved highly accurate water vapor profiles at different altitudes.
The proposed method was validated and assessed by four special schemes, in which the tomography area was adjusted to force the RS-voxels located in different tomographic boundaries. When the input of the tomography system was only integral signals, the average utilization of GNSS signals was 47.12%, and the mean RMSE and STD were 1.63 and 1.61 g/m 3 , respectively. When side rays were combined, the former was increased to 80.54% and was stable for the four schemes. The latter was improved to 1.16 and 1.15 g/m 3 , respectively. In the four experimental comparisons, whether in terms of RMSE or STD, it was observed that there was a significant improvement of the tomography results with the side signals. Furthermore, a special comparison of the proposed method and one existing correction approach was performed based on the three other schemes. Compared to the currently optimized model, the tomography results of the proposed model showed mean accuracy improvement rates of 15.62% and 14.40% in terms of the RMSE and STD, separately, which highlights the advantages of the proposed HFM in estimating the side SWV observations and retrieving the three-dimensional water vapor field.
In addition, the proposed HFM method provides a new way to address the tropospheric delays. From the perspective of isotropic and anisotropic components, some issues resulting from atmospheric water vapor may be skillfully handled, e.g., tropospheric delay correction in InSAR and VLBI, which may inspire new concepts for relevant researchers. Future work will focus on assimilating multi-source remote sensing data into tropospheric tomography, e.g., MODIS, MERIS, and InSAR, providing high-resolution PWV maps, which may be useful for GNSS tomography.