An Improved Tropospheric Tomographic Model Based on Artificial Neural Network

Global navigation satellite systems (GNSS) tropospheric tomography can be used to build a three-dimensional water vapor field. In traditional tomography, the signals crossing from the four sides of the tomographic region are not utilized. To make the best use of these valuable side-crossing signals, an improved tomographic model based on back propagation artificial neural network (BP-ANN) is proposed. In the new tomographic model, the inside part of the slant wet delay (SWD) of the side-crossing signal is divided into two sections: the isotropic and anisotropic components. The former is estimated by the zenith wet delay multiplied by the mapping function multiplied by an isotropic scale factor using a BP-ANN model, and the latter is estimated by horizontal gradients of the SWD multiplied by an anisotropic scale factor using an empirical model. The new tomographic model is experimentally evaluated using the HK CORS network measurements for the period of 21 days from 1 to 21 August 2019. Statistical results show that the root mean square error (RMSE) of slant water vapor reconstructed from the improved model is reduced to 1.35 from 2.85 mm of the traditional model. Compared with the traditional/height factor models, the percentages of the reduction in the RMSE of the tomographic result derived from the new model are 16%/9% and 22%/16%, respectively, using radiosonde and ERA5 data as references. These results suggest a good performance of the new model for GNSS tropospheric tomography.

high spatial-temporal resolution, all-weather conditions, low cost, nearly real-time and high accuracy of the global navigation satellite systems (GNSS) technique, water vapor retrieved from GNSS technique has been widely used in numerical weather prediction models and the research of severe weather phenomena [1], [2], [3], [4], [5]. To obtain four-dimensional (4-D) water vapor information, the tropospheric tomographic method is proposed by Flores et al. [6]. The tropospheric tomographic method has been rapidly developed during the last 20 years.
The traditional tropospheric tomography divides the tomographic region into a certain number of voxels in both horizontal and vertical domains under the assumption that the water vapor density (WVD) or wet refractivity in each voxel during the period of the GNSS tomography is constant. According to the distance of each GNSS signal travelling inside the voxels, the tomographic equation can be formed for the GNSS tomography. The tomographic observation equation matrix can be established using the slant water vapor (SWV) or slant wet delay (SWD) along the ray path of all GNSS signals during the period of the tomography. However, the shape formed from the GNSS signal ray paths during the tomographic period is similar to an inverted cone which is inconsistent with the box-shaped tomographic region. The tomographic equation matrix becomes an ill-conditioned or ill-posed problem due to the geometry of the signal distribution. To solve this problem, some iterative and non-iterative methods have been developed to resolve a tropospheric tomographic system. The singular value decomposition (SVD) is used to solve the inverse problem of a tropospheric tomographic system by Flores et al. [6], and the SVD is commonly used in noniterative methods. The damped least squares method is adopted to obtain a better tomographic result [7], and a wet refractivity Kalman filtering approach is proposed by Gradinarsky and Jarlemark [8]. Bender et al. [9] use several algebraic reconstruction techniques including algebraic reconstruction techniques, multiplicative algebraic reconstruction techniques (MART) and simultaneous iterative reconstruction technique to reconstruct a three-dimensional (3-D) wet refractivity field in Germany and the best tomographic results are obtained by using the MART. In addition, the compressive sensing method [10], and the genetic algorithm [11], are used for tropospheric tomography. For these tropospheric tomographic methods, the following three approaches including adding empirical constraint equations, introducing supplementary observations, and increasing the number of voxels crossed by GNSS signals are proposed for the regularization of the ill-conditioned tomographic system. The details of these methods are as follows.
Flores et al. [6] adopt empirical constraints including horizontal smoothing, vertical smoothing, and boundary conditions. Song et al. [12] propose a Gaussian weight function to define the weights of the horizontal constraint for all the voxels at the same height layer. In general, an exponential function is used for the vertical constraint [13]. However, it is difficult to know the real atmospheric condition in the period of the tomography, a priori information is introduced into the tomographic system. Flores et al. [6] adopt radiosonde profiles for the background field. Xia et al. [14] select radio occultation data as the a priori background field. Benevides et al. [15] analyze the impact of priori fields on the tomographic result using the prior wet refractivity profiles of atmosphere infrared sounder (AIRS) and radiosonde data.
For the second method, supplementary observations from other technologies are introduced into the tomographic system. Benevides et al. [16] and Heublein et al. [10], [17] add observations from interferometric synthetic aperture radar (InSAR) data into the tropospheric tomographic system. Benevides et al. [18] and Zhang et al. [19], [20] use water vapor products from moderate-resolution imaging spectroradiometer as supplementary observations for the GNSS tropospheric tomographic system. Zhang et al. [21] add water vapor products of the Fengyun-4A satellite into the GNSS tropospheric tomographic system. Chen and Liu [22] introduce water vapor observations from radiometers and sun photometers into the tropospheric tomographic system. Since the temporal/spatial resolutions of these observations are much lower than GNSS data, the number of the observations which can be added to the GNSS tomographic system is limited. In the traditional tomographic model, the GNSS signals crossing from the top boundary of the tomographic region are used, whilst the signals crossing from the side faces of the tomographic region are discarded. Therefore, recent studies tend to make the best use of the side-crossing signals and multi-GNSS signals in tropospheric tomography. Benevides et al. [23] and Dong et al. [24] utilize multi-GNSS observations to reduce the number of empty voxels. Rohm and Bosy [25] propose a method that used the University of New Brunswick version 3 modified model to obtain the SWD of the side-crossing signals inside the tomographic region. Benevides et al. [18] adopt an exponential law to represent the diminishing of the wet refractivity along the vertical direction and estimate the SWD inside the tomographic region for the side-crossing signals. Yao and Zhao [26] develop a water vapor unit index model based on both radiosonde data and the GNSS signals crossing from the top boundary of the tomographic region to estimate the initial WVD constraint for every height layer. Zhao et al. [27] develop an exponential height factor model to account for the side-crossing signals by expanding the tomographic region in the horizontal domain. Zhao et al. [28] improve this exponential model based on the Global Pressure and Temperature 2 wet model. By dividing the SWD into two parts-isotropic and anisotropic components, Zhang et al. [29] construct a height factor model (HFM) based on 30-year radiosonde data. In addition, Chen and Liu [30] and Yao and Zhao [31] use the optimal voxel distribution method to increase the number of voxels crossed by signals.
Water vapor varies highly in the spatial and temporal domains. However, in all the approaches using side-crossing signals mentioned above, the spatiotemporal characteristics of water vapor are not considered. In this study, the spatiotemporal characteristics are considered to develop the new tomographic model using the top-crossing and side-crossing signals. For the improved tomographic model, ERA5 hourly (0.125°× 0.125°) grid data are used to construct an improved isotropic scale factor model. The isotropic components of the SWD of side-crossing signals inside the tomographic region are estimated using the improved isotropic scale factor model. More specifically, the parameters of the isotropic scale factor model include not only the height, which is the only factor considered in the HFM but also the information of the spatiotemporal characteristics of the water vapor and the zenith wet delay (ZWD) at the GNSS station. In addition, the anisotropic component of the SWD of side-crossing signals inside the tomographic region is estimated by the experiential model. For the improved tomographic model, the SWD of the side-crossing signals inside the tomographic region is added into the tomographic system. Moreover, three tomographic schemes are designed for the performance evaluation of the improved tomographic model, and both radiosonde data and ERA5 grid data are selected as the two references.
The rest of this article is organized as follows. In Section II, the data and methods are introduced. Section III describes the new tomographic model using top-crossed and side-crossed GNSS signals. The validation of the new tomographic model is presented and analyzed in Section IV. Finally, Section V concludes this artice.

II. DATA AND METHODS
In this study, the following three types of data from the Hong Kong region are used.
1) GNSS data from the Hong Kong Satellite Positioning Reference Station Network. 2) Radiosonde data collected at the King's Park meteorological station. 3) ERA5 reanalysis grid data from the European Centre for Medium-Range Weather Forecast (ECMWF). The GNSS data are used to establish the tomographic equations. The radiosonde data are used to evaluate the results of the new tomographic model. The ERA5 reanalysis grid data are used to develop the isotropic scale factor model and validate the results of the new tomographic model. The methodology of GNSS tropospheric tomography and the artificial neural network used in this study are elaborated.
A. Data 1) GNSS Data: In this study, GNSS data with a 30 s sampling rate from 16 CORS stations in Hong Kong and three International GNSS Service (IGS) stations (BJFS, JGFN, and PIMO) during the 21 days from 1 to 21 August 2019 are selected for the tomographic experiment. The main reason for this selection is that this period contained sunny days, rainy days, and the value of water vapor is larger. In the GNSS data processing, the SWD is expressed as follows: Anisotropic component (1) where ZWD is the zenith wet delay estimated by the GNSS at MIT (GAMIT) / Global Kalman filter (GLOBK) software package [32], mf w (ε) and mf g (ε) are the mapping functions for the wet delay and gradients components, respectively, G w NS and G w EW are the horizontal gradients in the north-south and east-west directions, respectively, and ε, α are the elevation and azimuth angles of the signal path from the GNSS data processing, respectively. Since the SWV is used for tropospheric tomography in this study, the SWD retrieved from the GNSS signal needs to be transformed into the SWV by the following formula: where k 2 , k 3 , and R V are the specific gas constants [33]. T m is the weighted mean temperature which can be calculated by an empirical model from Liu et al. [34].
2) Radiosonde Data: In this article, the 21years (2000-2020) sounding data (twice per day), which are acquired from the Integrated Global Radiosonde Archive [35], are used at King's Park meteorological station. The used parameters of radiosonde data include the observations of temperature, pressure, water vapor pressure, and geopotential height at different pressure layers. The data are used to estimate the ZWD and WVD at the radiosonde station by the formulations [36], [37] where h s and h t are the heights of the bottom and top pressure layers, respectively, P w and T are the water vapor pressure (hPa) and temperature (K), respectively, and k 2 and k 3 are the ideal gas constants from Thayer [33] where ρ 0 is the WVD, P w and T are the water vapor pressure (hPa) and temperature (K), respectively, and R v is the specific gas constant for water vapor.
3) ERA5 Data: ERA5 is the latest reanalysis data product from ECMWF and provides meteorological data at 37 pressure layers from 1000 to 1 hPa, such as pressure, temperature, geopotential height, and specific humidity. The ERA5 hourly 0.125°× 0.125°grid data in the region of Hong Kong are used to calculate the ZWD in (3) and WVD in (4). The 9-years (2011-2019) ERA5 grid data [38] are used to construct the isotropic scale factor model contained in the improved tomographic model. The ERA5 grid data are also used to evaluate the isotropic scale factor model and the result of the tomography.

B. GNSS Tropospheric Tomography
The tropospheric tomographic region is divided into a certain number of voxels. The SWV of a GNSS signal ray inside the tomographic region can be calculated by the WVD multiplied by the distance of the signal ray path inside the voxels crossed by the signal [6] where SWV (mm) is the slant water vapor along the signal, i, j, k are the indices of the voxel in the longitude, latitude, and vertical domains, respectively, ρ i,j,k (g/m 3 ) and d i,j,k (km) denote the WVD and distance in the voxel (i, j, k) passed by the ray, respectively, and n is the number of all the voxels crossed by the ray. The observation equation system for all SWVs contained in the tomographic region can be expressed in the matrix form as follows: where SW V top includes the SWVs of the GNSS signals crossing from the top boundary of the tomographic region, A top is the coefficient matrix including the length of the top-crossing signals crossing through each voxel, and X is the column vector including the unknown WVD parameters in all voxels.
To obtain a better tomographic result, horizontal and vertical constraint equations are applied. The horizontal constraint for each height layer is a weighted mean WVD of all the voxels at the layer [6] w 1,1,k ρ 1,1,k + . . .
where w 1,1,k , w i−1,j−1,k . . . denote the weight coefficient of each voxel at the k height layer. ρ 1,1,k, , ρ i−1,j−1,k . . . denote the WVD of each voxel. i, j, and k are the indices of the voxel in the longitude, latitude, and vertical domain. The weight of a voxel is estimated using the Gauss weighting function [12]: Where i, j, and k are the indices of the target voxel in the three spatial dimensions, il, jl, and k are that of the other voxels at the same height layer, d il,jl,k is the distance between the target voxel and the other voxels, σ is the smoothing factor, and ne and nn are the numbers of voxels in the longitude and latitude domains, respectively. According to previous research [13], WVD usually decreases with height, and the variation trend is close to the exponential function. As a result, the exponential function is selected as the vertical constraint for the relationship between the WVDs of two vertical adjacent voxels with indices i, j, k and i, j, k + 1 where ρ (g/m 3 ) denotes the WVD, h (km) is the height of the vertical layer, SH is the scale height of the water vapor, for which a 2 km value is selected in this study, and i, j, k and i, j, k + 1 are the indices of two adjacent voxels in the vertical domain.
As we all know, the error in the SWV increases with the reduction in the elevation angle of the signal ray path. Hence, a function that expresses such a relationship is adopted for the weight of the observation equation matrix. The identity matrices are adopted for the weight of the horizontal and vertical constraints [17], [29].
The tomographic observation equation matrix is as follows: ⎛ where SW V top denotes the observation vector from the topcrossing signals, A top is the observation matrix of top-crossing signals, H and V are the horizontal and vertical constraint matrices, respectively, P top , P hr , P vr denote the weight matrices of observation equations for top-crossing signals, horizontal and vertical constraints, respectively, and X denotes the vector of the unknown WVD parameters of all voxels in the tomographic region.

C. BP-ANN
In previous studies, the feed-forward BP neural network (BP-NN) technique has been popularly applied in various geoscientific research [39], [40], [41], since it can solve non-linear problems with different structures and precisions. A BP-NN model is composed of one input layer, one or more hidden layers, and one output layer, and each layer contains at least one neuron. The training of the BP-NN has two processes including the forward propagation of input parameters and the backpropagation of the "loss" value. The forward propagation of input parameters is used to obtain predicted values. The backpropagation of the "loss" value which is the error of desired outputs [42] is used to update the weight matrix and offset vector. The forward propagation of input parameters can be simply expressed as follows: where X bp and Y bp are the input and predicted vectors, respectively, w bp 1 and a bp 1 are the weight matrix and offset vector between the input layer neurons and hidden layer neurons, respectively, w bp 2 and a bp 2 are the weight matrix and offset vector between the hidden layer neurons and output layer neurons, respectively, and f 1 and f 2 are two activation functions.
The purpose of the backpropagation is to update the weight matrix and offset vector for the forward propagation using the error between the desired value and the predicted result where w bp m+1 and a bp m+1 are the weight matrix and offset vector of the m+1 th iterative result, w bp m and a bp m are the weight matrix

D. Validation Method
To evaluate the accuracy of the isotropic scale factor from the new model, three error evaluation indicators, including bias, root mean square error (RMSE), and standard deviation (STD) are used. The accuracy of WVD from the new tomographic model is evaluated using three error evaluation indicators including bias, RMSE, and mean absolute error (Mae) where N is the number of samples, y and y are the isotropic scale factor / WVD from the tomographic model and validation data from radiosonde and ERA5 data, respectively, and μ is the mean value of the differences between two datasets.

III. NEW TOMOGRAPHIC MODEL BASED ON BP NEURAL NETWORK
In the traditional tropospheric tomographic model, only those signals that pass from the top boundary of the tomographic region are utilized. The signals that cross from the side faces of the tomographic region are ignored, which leads to a "waste" of these signals and a decrease in the number of voxels penetrated by these signals. In this study, a new method which improves the use of side-crossing signals is developed for tropospheric tomography. In Fig. 2(a), the 3-D view of the top-crossing and side-crossing signals is represented.
In the observation equation of GNSS data processing, the slant tropospheric delay can be expressed as a function of the zenith total delay (ZTD), and the ZTD is estimated from GNSS data processing. The ZWD can be calculated by the ZTD minus the zenith hydrostatic delay (ZHD) which can be obtained using the Saastamoinen model. The SWD can be obtained by the ZWD multiplied by a mapping function (the isotropic component) and horizontal gradients (the anisotropic component). In Fig. 2(b), the SWD of a side-crossing signal (PN), which is a part of the signal (PR) inside the tomographic region, can be expressed by its wet delay in the zenith direction (ZWD from P to N1) multiplied by a mapping function (the isotropic component) and its horizontal gradients (the anisotropic component). In this study, an isotropic scale factor model and an anisotropic scale factor model are developed to estimate the isotropic and anisotropic scale factors, respectively. The two scale factors are the ratio of the part of the signal inside the tomographic region to the whole signal for the isotropic and anisotropic components, respectively. Therefore, the SWD of a side-crossing signal inside the tomographic region can be expressed by the isotropic scale factor multiplied by the isotropic component plus the anisotropic scale factor multiplied by the anisotropic component. The procedure of the scale factor models is detailed in Section III-A and B. The tomographic equations of the new tomographic model are shown in Section III-C.

A. Isotropic Component of SWD From Side-Crossing Signals
The isotropic component of an SWD is obtained by the ZWD multiplied by a mapping function. The isotropic component of the SWD of a side-crossing signal inside the tomographic region can be obtained by the isotropic scale factor multiplied by the ZWD multiplied by the mapping function. The main idea of the isotropic scale factor model is as follows. The side-crossing signal (RP) in Fig. 2(b) is used as an example to show the process of developing the isotropic component model. 1) Obtaining the isotropic component of the SWD for the whole signal. In Fig. 2(b), the RP is a whole signal path and h P M1 (green vertical line) is the height of the tropopause at the GNSS station P. The isotropic part of a SWD and the zenith wet delay for the whole signal path (RP) can be expressed as follows: where ZWD total is the whole zenith wet delay at GNSS station P, RP is the signal path, h P and h M 1 are the height of the GNSS site and the top of the tomographic model, respectively, N w is the wet refractivity, SWD(RP ) iso is the isotropic part of the SWD of the signal (RP), mf w (ε) is the wet mapping function, and ε is the elevation angle of the signal (RP).
2) Obtaining the isotropic component of the SWD for the signal inside the tomographic region. The signal (PN) is a segment of the path of the signal (RP), and the wet delay of the signal (PN) can be expressed by the ZWD which is the integral of the wet refractivity from h P to h N 1 (h N 1 is the height of signal (RP) crossing from the side face of the tomographic model) multiplied by a mapping function. The two equations are as follows: where ZWD h P N1 is the partial ZWD, h P and h N 1 are the height of the GNSS site and the signal crossing from the tomographic model, respectively, SWD(PN) iso is the isotropic part of the SWD of signal (RP) inside the tomographic model, i.e., the SWD of signal path (PN). 3) Defining isotropic scale factor. According to (18) and (20), the isotropic scale factor of the side-crossing signal can be expressed as follows: where λ iso (h PN1 ) is the isotropic scale factor of the station P at the height h PN1 . Therefore, SWD(PN) iso can be expressed as follows: where SWD(PN) iso is the isotropic part of the SWD of the signal (RP) inside the tomographic model, i.e., the SWD of signal path (PN), λ iso is the isotropic scale factor at the height h PN1 and mf w (ε) is the mapping function, ε is the elevation angle of the signal (RP), and ZWD total is the whole zenith wet delay at GNSS station P. According to (22), the isotropic scale factor general form is as follows: where λ iso (h) is the isotropic scale factor at the height h, ZWD h and ZWD total are the partial ZWD inside the tomographic region and the whole ZWD at the station, respectively, and h is the height of the signal crossing from the tomographic model. 4) Calculating isotropic scale factor using ERA5 data.
In this research, ERA5 hourly grid data (with the horizontal resolution of 0.125 × 0.125°) in Hong Kong during the period from 2011 to 2019 are used to construct the regional isotropic scale factor model. The grid points used to develop the model are shown in Fig. 3.
To develop the isotropic scale factor model, ZWD h and ZWD total in (24) need to be calculated using ERA5 data with (19), (21), and (3) (26) where ZWD total and ZWD h l are the whole ZWD and partial ZWD from h S to h l at the grid point. h S , h TRP , and h l are the height of the lowest pressure layer of the ERA5 grid point, the height of the tropopause at the grid point, and the height of the other pressure layers of the ERA5 grid point. Then, the isotropic scale factors (λ iso (h)) of each grid point in the height of each pressure layer (under the height of the tropopause) can be calculated. These isotropic scale factors are used to develop the isotropic scale factor model. 5) Developing the isotropic scale factor model using BP-ANN. In the previous study, the exponential function is adopted to fit the scale factor model. However, the exponential model, like any fitting model, only reflects the general variation trend of λ iso (h) with low temporal and spatial resolutions. To construct a refined isotropic scale factor model, the BP-ANN is applied in this study. Fig. 4 shows the structure of the input layer, hidden layers, and output layer for the new isotropic scale factor model. For this model, the isotropic scale factor (λ iso (h)) is set as output variables, and seven variables including year, day of year (doy), hour of day (HoD), latitude, longitude, height (the height of the pressure layer), and the whole ZWD are set as the input variables. In this study, the isotropic scale factors from step 4) are used as the output data and the seven variables of each scale factor make up the input data. The input and output data are used to train the BP-ANN model. In the training process of the BP-ANN model, the Levenberg Marquardt algorithm is used to update the weight matrices. The ratios of the training, validation, and test datasets are 70%, 15%, and 15%, respectively. The threshold for error precision of the training process is set to 0.001. At present, some experimental equations are used to calculate the number of nodes of the hidden layer, but this value only is a reference value. To construct a better model, many schemes are tested for identifying the optimal number of nodes in the hidden layer, the optimal number of the hidden layers, and the optimal active functions. Finally, the structure with 7, 4, and 2, combined with the active functions of the logsig for the three hidden layers of the BP-ANN are determined. 6) Obtaining λ iso of GNSS site using the new model In a tropospheric tomography, information of time (year, doy, HoD), location of the GNSS station (longitude, latitude, and height of side-crossing signal crossing from the tomographic model) and the whole ZWD of each GNSS site can be obtained as the input value of the new isotropic scale factor model. Then, the output value of the new scale factor model, i.e., the λ iso of each side-crossing signal can be obtained. According to (23), the isotropic component of the SWD of the side-crossing signal can be expressed where SWD(h) iso is the SWD of the side-crossing signal path under the height h, λ iso (h) is the isotropic scale factor of the GNSS site at the height h, mf w (ε) is the mapping function of wet delay, ε is the elevation angle of the signal path, and h is the height of the side-crossing signal crossing from the tomographic region.

B. Anisotropic Component of SWD From Side-Crossing Signals
The horizontal gradients of the wet delay are another part of a SWD, and numerous studies have investigated the characteristic of the horizontal gradients using very long baseline interferometry [43], [44]. The horizontal gradient of the β direction can be expressed as the integral of the wet refractivity gradient with height where G β is the azimuthal asymmetry of the atmosphere, i.e., the horizontal gradient in the β direction, N β (h) is the wet refractivity gradient at height h, and β is the azimuth direction including north-south and east-west directions.
Since it is difficult to directly obtain N β (h) in a tomographic modelling process, it is assumed that the refractive gradient varies exponentially with height in this study, as the work of Zhang et al. [29] where N β (h) is the wet refractive gradient of the β direction at the height h, N S β is the wet refractive gradient of the β direction at the surface of the site, β is the azimuth directions including north-south and east-west directions, SH is the scale height of water vapor (1-2 km).
Substituting (29) into (28) leads to the following: where G β is the azimuthal asymmetry of the atmosphere, i.e., horizontal gradient, N S β is the wet refractive gradient of the β direction at the surface of the site, β is the azimuth directions including north-south and east-west directions, and SH is the scale height of water vapor.
Equation (30) can be expressed by partial integral calculus Therefore, horizontal gradient (anisotropic component) can be obtained Replacing h top with h leads to the following: where G h top β is the whole horizontal gradient at the β direction, G h β is the partial horizontal gradient of the β direction at the height h, β is the azimuth directions including north-south and east-west directions, and h top and h are the height of the tropopause and side-signal crossing from the tomographic model, respectively.
The main idea of the anisotropic scale factor model is as follows. The side-crossing signal (PR) in Fig. 2(b) is used as an example to show the process of developing the anisotropic component model. According to the anisotropic component term, i.e., the last term in (1), the anisotropic component of the signal (PR) can be expressed where SWD(PR) aniso is the anisotropic component of the SWD for the signal (PR), G h P M NS and G h P M EW are whole horizontal gradient components in north-south and east-west directions which are calculated by (33), h P M is the height of the tropopause, mf g (ε) is the mapping function of the horizontal gradient components, and ε, α are the elevation and azimuth angle of the signal (PR), respectively.
Similarly, for the part of the side-crossing signal (PR) that is within the tomographic region, i.e., PN in Fig. 2(b), its anisotropic component can be expressed as follows: where G h PN NS and G h PN EW are the north-south and east-west horizontal gradient components of the signal (PR) within the tomographic region, i.e., the north-south and east-west horizontal gradient components from the ground up to the height h PN , and they can be obtained from (34), h PN is the height of the side-crossing point of the signal (PR), mf g (ε) is the mapping function of the horizontal gradient components, and ε, α are the elevation and azimuth angle of the signal (PN), respectively.
Substituting (33) and (34) into (35) and (36), respectively, the anisotropic scale factor becomes as follows: The general form of (37) can be expressed by replacing h P N with h λ aniso (h) = where λ aniso (h) is the anisotropic scale factor, h is the height of the point from which the side-crossing signal crosses the tomographic region, h top is the height of the tropopause, which is also the same as the top of the tomographic region, and SH is the scale height of water vapor.
Considering (35) and (37), the anisotropic component of the signal (PN) shown in Fig. 2(b) can be expressed as follows: where SWD(PN) aniso and SWD(PR) aniso are the anisotropic components of SWDs from the signal (PN) and signal (PR).  Therefore, for any side-crossing signal, the anisotropic component of the SWD of the side-crossing signal inside the tomographic region can be obtained where SWD(h) aniso is the anisotropic component of the SWD of the side-crossing signal inside the tomographic region, SWD(h top ) aniso is the whole anisotropic component of the SWD at the site location, λ aniso (h) is the anisotropic scale factor at the height h, G h top NS and G h top EW are whole horizontal gradient components in north-south and east-west directions, h top is the height of the tropopause, which is also the same as the top of the tomographic region. h is the height of the side-crossing signal crossing from the tomographic region, and ε and α are the elevation and azimuth angle of the signal (PR), respectively.

C. New Tropospheric Tomographic Equation Matrix
Fig . 5 shows the flowchart of the new tropospheric tomographic model using top-crossed and side-crossed signals.
A whole SWD can be divided into the isotropic and anisotropic components in (1). Therefore, for a side-crossing signal, the SWD inside the tomographic model (SWD(h)) can be estimated using the isotropic and anisotropic scale factors at the height at which the side-crossing signal crosses from the tomographic model. The isotropic and anisotropic scale factors can be obtained by the isotropic scale factor model (BP-ANN model) and (38), respectively. Then, this SWD(h) can be expressed by merging (27) and (40) SWD(h) can be transformed into SWV(h) using (2). The observation equations for the side-crossing signals in the matrix form are as follows: where SW V side is the observation vectors from the sidecrossing signals, A side is the observation matrix of sidecrossing signals, and X is the same as that in (6). Therefore, the new tropospheric tomographic equation system can be expressed as follows: where SW V top and SW V side denote the observation vectors from the top and side-crossing signals, respectively, A top and A side are the observation matrices of the top and side-crossing signals, respectively, P top , P side P hr , P vr denote the weight matrices of observation equations for top-crossing signals, sidecrossing signals, horizontal and vertical constraints, respectively, and X denotes the vector of the unknown WVD parameters of all voxels in the tomographic region.

A. Evaluation of the Isotropic Scale Factor Model
In this section, the result of λ iso (h) from the new model are evaluated using radiosonde and ERA5 grid data in the year 2020 as the references. To exhibit the effect of the new model, the isotropic scale factor of the HFM is compared. In Zhang et al. [29], the exponential law defined in (44) is adopted to fit the isotropic scale factor model based on 30-year radiosonde data, and its formula is as follows: where λ exp iso (h) is the isotropic scale factor from the HFM, a 1 , a 2 , b 1 and b 2 are the fitting coefficients estimated by the leastsquares method, and h is the height of the signal crossing from the tomographic model.
The exponential model and the new model are called ISFEXP and ISFBP, respectively, in this study. In Table I

1) Validation With Isotropic Scale Factor Derived From Radiosonde Data:
The λ iso (h)−ISFEXP (the exponential model) and λ iso (h)−ISFBP (the new model) at the position of the radiosonde station for the year 2020 are calculated and the predicted results from the two models are compared with the monthly radiosonde result called λ iso (h) − RS. Their (a) bias, (b) RMSE, and (c) STD of the ISFBP and ISFEXP models are shown in Fig. 6. Fig. 6(a) indicates that the 8 monthly biases of λ iso (h)− −ISFEXP (blue) are smaller than that of λ iso (h) − ISFBP. The mean of bias values of λ iso (h) − ISFEXP is closer to zero, meaning that the ISFEXP result reflects the general variation trend of λ iso (h). The biases of nine months from the ISFBP are negative, meaning that the ISFBP underestimates λ iso (h) when radiosonde data are used as the reference. The means of all the monthly results shown in Fig. 4 are listed in Table II. It is noted that the mean bias of λ iso (h)−ISFEXP is smaller than that of λ iso (h)−ISFBP using the radiosonde data as the reference. This can be because the ISFEXP model is fitted from 20-year radiosonde data and reflects the general variation trend of the scale factor. However, λ iso (h)−ISFBP is based on the sample of ERA5 hourly grid data during the 9 years of 2011-2019.
In Fig. 6(b) and (c), all the monthly RMSEs and STDs of the ISFBP results, except for January and April, are smaller than IS-FEXP, and the statistical results are shown in Table II. The mean of the RMSEs from ISFEXP and ISFBP are 0.050 and 0.044, respectively. The percentages of the reduction in the RMSE made by ISFBP over ISFEXP is 14%. This implies that, in terms of RMSE, λ iso (h)−ISFBP outperforms λ iso (h)−ISFEXP. In addition, the two models have smaller RMSEs in the summer than in the winter, probably due to the larger concentration of water vapor in the summer.

2) Validation With Isotropic Scale Factor Derived From ERA5 Data:
Since only one radiosonde station is deployed in Hong Kong, the isotropic scale factor only in the location near the radiosonde station (the so-called co-location) can be validated when radiosonde data are used as the reference. If ERA5 grid data are used as the reference, all grid points in the region can be used for the validation. In this section, data from 12 ERA5 grid points in the Hong Kong region in the year 2020 are adopted as the reference to evaluate the above two selected models. Fig. 7 shows the monthly statistical results of all 12 grid points. For two models, nine monthly biases of isotropic scale factors are positive. Therefore, the two isotropic scale factor models overestimate the isotropic scale factor when ERA5 data are used as the reference. The nine monthly biases of ISFBP are closer to zero than that of ISFEXP. In Table III, which is the statistics of   FIG. 7 the results shown in Fig. 7, the mean of the 12 monthly biases from ISFBP and ISFEXP are 0.006 and 0.012, respectively, i.e., ISFBP is better in terms of the bias.
From Fig. 7(b) and (c), except for April, all the rest months' RMSEs and STDs from ISFBP are smaller than that of ISF-EXP. The maximum percentage of improvement in the RMSE made by ISFBP over ISFEXP is 46% in September, and the minimum value is 3% in January. The number of months with the improvement percentage above 20% is six, and the number of months with the improvement percentage above 10% is 10. However, the improvement percentage is −7.30% in April. The reason may be that the fluctuation of the scale factor in April is larger than that of the other months. In Table III, the mean of all the 12 monthly RMSEs of the ISFBP and ISFEXP are 0.045 and 0.055, respectively, and the mean of the improvement percentages is 22% made by ISFBP over ISFEXP.
To investigate the difference in the RMSEs and biases from the same model but from different grid points, or from the same grid point but from different models, the means of the 12 monthly RMSEs (the left column) and biases (the right column) at 12 grid points from ISFEXP (the upper row) and ISFBP (the bottom row) are shown in Fig. 8.
In the two subfigures (a) and (c) for the RMSE results, the former (ISFEXP) increases with the increase in latitude. The latter (ISFBP) varies little at all the 12 grid points, and the results of the 12 grid points are all considerably smaller than that of the former. Hence, in terms of accuracy, ISFBP outperforms ISFEXP. Similar results are also shown in the two subfigures in the right column for the bias results. Therefore, compared with the model that is based on the sample data from only one radiosonde station, the main advantage of using ERA5 grid data to develop the λ iso (h) model is a better spatial resolution.

B. Tomographic Experiment and Data-Processing Strategy
In this section, experimental results for the improved tropospheric tomographic model, the HFM, and the traditional model are compared. The radiosonde data (only the two epochs at 00:00 and 12:00 UTC each day are available) from the radiosonde GNSS data (with a 30s sampling rate) from 16 GNSS stations in the Hong Kong Satellite Positioning Reference Station Network and three IGS stations (BJFS, JFNG, and PIMO) during the 21 days from 1 to 21 August 2019 are selected as the sample data of the tomographic modelling for Hong Kong. Fig. 9 illustrates the spatial distribution of the 16 GNSS stations and the radiosonde station. To adapt to the above mentioned two epochs of radiosonde data, the ERA5 grid data at the same epochs on Since the water vapor content over a site varies slowly within a certain period [45], a 30-min is selected for the tomographic time window in this study, i.e., the temporal resolution of the tomographic reconstruction is 30 min.
The GAMIT/GLOBK (v.10.7) software package and the following strategies are used for the GNSS data processing. The ZTD and delay gradient are estimated with intervals of 0.5 and 2 h, and the ZTD is interpolated to 5 min by the piecewise linear interpolation. The ZHD is calculated by the Saastamoinen model [46], with the input of surface pressure measurements. Equation (1) is used to calculate the SWD. The VMF1 model is selected to calculate the mf w (ε). The mf g (ε) used in [47] is adopted. The SWD is transformed into the SWV by the conversion factor from (2).

C. Evaluation of Improved Tomographic Model
To evaluate the tomographic results, the number of the signals utilized and the residual of the reconstructed SWVs during the 21 days studied are used. For a comparative analysis, the traditional model which utilizes only those signals crossing from the top boundary of the tomographic region is also evaluated. Fig. 10 shows the 3-D distributions of the mean number of signals crossing the tomographic voxels for the traditional tropospheric tomographic model (a) and the improved tropospheric tomographic model (b), where white voxels indicate none-signal-crossing voxels. (c) shows the daily mean of the numbers of the signals used in the tomographic experiment and (d) shows the daily mean of the percentages of the voxels that are crossed by at least one signal ray. In the two subfigures (a) and (b), the number of empty voxels from the improved model is less than that of the traditional model at most vertical height layers, but the number of the empty voxels at the lowest height layer from two models is the same. This can be because the signals with the elevation angle under 15°are discarded. Therefore, it is difficult for the voxels in the lowest height layer to be crossed by the signal from the station which is in the other voxels at the same height layer. It is noted that only the data that are within the two 30-min periods of 00:00-00:30 UTC and 12:00-12:30 UTC are used to obtain these daily means for (c) and (d). The daily mean of the numbers of the signals and the percentage of the voxels crossed by signals from the improved model (red) are both larger than that of the traditional model (blue). The means of the 21 daily percentages of the voxels crossed by signals shown in (d) are 87% and 76% from the improved model and the traditional model. The percentage of the voxels crossed by signals from the improved tomographic model represents an improvement of 11% over the traditional tomographic model, which may be beneficial to the rank deficiency of the tomographic equations. Fig. 11 shows the variation in the residual of the improved model (red) and traditional model (blue) with the elevation angle of the signal. The residual is the difference between the SWV reconstructed from the models and the SWV obtained from the GNSS data processing. The residuals of reconstructed SWVs are in the range between −30 and 25 mm, and the absolute values of the residuals decrease with the increase in the elevation angle.
The residuals of the improved model are smaller than that of the traditional model, especially from low elevation angles, which indicates a better performance of the improved model.
The bias/RMSE of the reconstructed SWVs from the improved and traditional tomographic models are also calculated, and the SWVs at 00:00 and 12:00 UTC estimated from the GAMIT/GLOBK (v.10.7) software are used as the reference. The bias/RMSE of the improved and the traditional models are 0.009/1.345 and −0.573/2.852 mm, respectively, which is equivalent to a 53% reduction in the RMSE made by the improved model over the traditional model.

D. Evaluation of Improved Tomographic Model Using Radiosonde Data as Reference
For comparative analyses, three tomographic models including the improved model proposed in this study, the HFM [29], and the traditional model are tested, and their descriptions are shown in Table IV. The radiosonde data from radiosonde station 45004 (located at King's Park, Hong Kong) (see Fig. 9) during the 21 days studied are used as the reference for the three models' results. Fig. 12 shows the tomographic results, i.e., WVD profiles derived from the three tomographic models and the radiosonde data at 12:00 UTC on DOY 230 (rainy day, left pane) and DOY 220 (rainless day, right pane) in 2019. This figure only shows the results of two selected days for an illustration of WVD profiles. In Fig. 12(a)-(c) and (e)-(g), the 3-D water vapor fields at 12:00 UTC on DOY 230 (rainy day) and DOY 220 (rainless day) in 2019 are presented. The values of WVD from the improved model are smaller than that of the traditional and HFM models in the low layers from 1 to 3 layers. In the two subfigures (d) and (h), the WVD profiles from the three models are all close to the reference profiles calculated by (4) (blue), but the profiles from the improved model (red) are the closest, especially at the low layers (under 1 km) and high layers (above 4 km), roughly. In addition, the RMSE of the WVD profiles from the improved, HFM and traditional models on the selected two days are also calculated and the results are 1.50/0.89, 1.66/1.91, and 2.08/2.10 g/m 3 , respectively, which gives the numerical improvement made by the improved model over the other two models. Fig. 13 shows the RMSE of the WVDs from three tomographic models at 42 tomographic epochs using radiosonde data as the reference. For 35 tomographic epochs, the RMESs of the tomographic results made by the improved tomographic model are better than that of the traditional tomographic model. The RMSEs of tomographic results from the improved tomographic model are better than that of the HMF at 29 tomographic epochs.
The mean RMESs of 42 tomographic epochs are 1.64, 1.51, and 1.38g/m 3 for the traditional, HFM, and improved tomographic models, respectively.  This means that the improved model outperforms the other two models, and the HFM outperforms the traditional model using the radiosonde data as the reference.
Similar to the RMSE shown in Fig. 13, the other two statistical results including Mae and bias are also calculated in Table V.
The mean percentages of the reduction in the RMSE made by the improved model over the traditional and HFM models are 16% and 9%, respectively. In addition, the mean Mae of the WVD from the improved model is better than the other two models. However, the mean bias of the HFM is less than that of the improved model. It may be because the HFM is constructed based on the radiosonde data, which is also the reference of the model results. To further study the stability of the new tomographic model under different weather conditions, the mean RMSE of the WVD under nonrainy (9 days) and rainy (12 days) conditions are also shown in Table V. The statistical results show that the new tomographic model has obvious improvement in RMSE than the traditional and HMF models under both nonrainy (9 days) and rainy (12 days) conditions, and the improvement in the rainy days is better than that of the nonrainy days.

E. Evaluation of Improved Tomographic Model Using ERA5 Data as Reference
In this section, the WVDs calculated by (4) from ERA5 grid data (horizontal resolution of 0.125 × 0.125) at 0:00 and 12:00 UTC on the 21 days studied are used as the reference for the evaluation of the three tomographic models. The mean RMSEs of the WVD profiles at all tomographic epochs at 12 grid points are shown in Fig. 14, where the three subfigures (a), (b), and (c) are the results of the traditional model, the HFM, and the improved model, respectively. The RMSEs in each subfigure increase with the increase of longitude. For the comparison among the three subfigures, the RMSEs at all the 12 grid points in (c), i.e., the improved tomographic model, are significantly better than the other two models. This result is similar to the case when radiosonde data are used as the reference. Fig. 15 shows the mean RMSEs of the 12 grid points on each tomographic epoch from three models. For 34 tomographic epochs, the mean RMESs of the 12 grid points made by the improved tomographic model are better than that of the traditional tomographic model. The mean RMSEs of 12 grid points from the improved tomographic model are better than that of the HMF at 28 tomographic epochs. The mean RMESs of 42 tomographic epochs are 1.31, 1.22, and 1.03 g/m 3 for the traditional, HFM, and improved tomographic models, respectively. The three values indicate that the overall accuracy of the improved model significantly outperforms the other two models, and the HFM slightly outperforms the traditional model.
Similar to the RMSE shown in Fig. 15, the other two statistical results including Mae and bias are also calculated in Table VI. The mean percentages of the reduction in the RMSE made by the improved model over the traditional and HFM models are 22% and 16%, respectively. In addition, the mean Mae of the WVDs from the improved model is better (smaller) than the other two models. However, the mean bias of the HFM (−0.08) is less than that of the improved model (−0.18), meaning that both the models underestimate the WVD, and  the improved model is worse. The reasons are discussed in the previous section (radiosonde data as the reference) When ERA5 data are used as the reference, the new tomographic model still outperforms the HFM and traditional models under both rainy and nonrainy conditions, and the new model has greater improvement in the rainy days than that of the non-rainy days. The results are similar to the results using radiosonde data as the reference. These results mean that the new tomographic model is more suitable for the tomography of large water vapor variation.
To investigate the relationship between the RMSEs of each model results with altitude, the mean RMSEs of the WVDs at each height layer are calculated and the results are shown in Fig. 16. The results of the improved tomographic model (red) at most height layers are smaller than that of the other two models, especially at the layer that is closest to the surface. This means the improved model performs best at the lowest height layer. Moreover, at low-to-middle layers, i.e., from 1.8 to 3.4 km, roughly, the HFM and the improved model results are similar, and they are much better than the traditional model. At high altitude layers, i.e., above 3.4 km, roughly, the improved model slightly outperforms the other two models.

V. CONCLUSION
The traditional tropospheric tomographic model only uses the GNSS signals crossing from the top boundary of the tomographic region of interest, while the GNSS signals crossing from the four side faces of the tomographic region are not considered. This leads to a waste of side-crossing signals and a decrease in the number of voxels that are penetrated by these signals. In this study, an improved tropospheric tomographic method based on an artificial neural network is developed using the signals not only crossing from the top boundary of the tomographic region but also crossing from the four side faces of the tomographic region.
The accuracies of the isotropic scale factor from the improved tomographic model and that of the HFM model (both models use side-crossing signals) are compared using two references from radiosonde and ERA5 data. The statistical results show that the mean percentages of the reduction in the RMSE made by the improved model over the HFM are 14% and 22% using radiosonde and ERA5 data as the references, respectively. This suggests that the isotropic scale factor from the improved tropospheric tomographic model is better than that of the HFM. Moreover, the scale factor of the anisotropic component is estimated under the assumption that the variation in the refractivity gradient in the vertical direction is an exponential function of altitude.
To evaluate the improved tropospheric tomographic model, GNSS data from 16 GNSS stations in Hong Kong during the 21 days from 1 to 21 August 2019 are adopted to obtain the tomographic results, and the radiosonde and ERA5 data in the same region and time are used as the references. For the improved tropospheric tomographic model, the mean of the utilization percentages of voxels increases from 76% of the traditional model to 87% of the improved model. The percentage of the reduction in the RMSE of the reconstructed SWVs made by the improved model over the traditional model is 53%. The statistical results of the WVD profiles show that the mean percentages of the reduction in the RMSE made by the improved model over the traditional/HFM models are 16%/9% and 22%/16% using radiosonde and ERA5 data as the references, respectively. Similarly, the comparison of the RMSEs of the WVD profiles at 16 vertical height layers between the three models also shows the superiority of the improved model. Especially, the accuracy of the improved tomographic model has obvious improvement over the HMF and the traditional tomographic model in the lower height layer. These results suggest that the tomographic results from the improved tomographic model are more accurate than that from the other two models.
Our future work will be mainly on utilizing more water vapor data from various technologies, e.g., multi-GNSS, GNSS-R, InSAR and virtual signals from numerical weather prediction models into the tropospheric tomographic system for more refined or better tomographic results.
Moufeng Wan received the Ph.D. degree in geodesy and surveying engineering from the China University of Mining and Technology, Xuzhou, China, in 2023.
His research interests include robust and adaptive nonlinear Kalman filtering and its application in GNSS navigation. His current research interests include GNSS tomography, four-dimensional atmospheric water vapor modelling, meteorological applications of GNSS, and spatial distribution and variation mechanism of small/medium scale tropospheric water vapor. His main research interests include radio occultation and GNSS meteorology. Shangyi Liu received the B.Sc. degree in surveying engineering from the School of Geodesy and Geomatics, Wuhan University, Wuhan, China, in 2020. He is currently working toward the M.S. degree in surveying and mapping engineering in the School of Environment Science and Spatial Informatics, China University of Mining and Technology, Xuzhou, China.

Peng
His main research interests include global navigation satellite system meteorology.
Andong Hu received the Ph.D. degree in geospatial science from Royal Melbourne Institute of Technology University, Melbourne, Vic, Australia, in 2020.
He is currently working as a Postdoc with CU Boulder, Boulder, CO, USA, and NOAA Space Weather Prediction Center, Boulder. His current research interests include machine learning algorithm development and innovative applications of space weather modeling for solar wind prediction, solar wind and magnetosphere coupling, and ionosphere analysis.