Tropospheric Refractivity Proﬁle Estimation by GNSS Measurement at China Big-Triangle Points

: Atmospheric radio refractivity has an obvious inﬂuence on the signal transmission path and communication group delay effect. The uncertainty of water vapor distribution is the main reason for the large error of tropospheric refractive index modeling. According to the distribution and characteristics of water vapor pressure, temperature, and pressure, which are the basic components of the refractive index, a method for retrieving atmospheric refractivity proﬁle based on GNSS (Global Navigation Satellite System) and meteorological sensor measurement is introduced and investigated in this study. The variation of the correlation between zenith wet delay and water vapor pressure is investigated and analyzed in detail. The partial pressure proﬁles of water vapor are retrieved with relevance vector machine method based on tropospheric zenith wet delay calculated by single ground-based GPS (Global Positioning System) receiver. The atmospheric temperature and pressure is calculated with the least square method, which is used to ﬁt the coefﬁcients of the polynomial model based on a large number of historical meteorological radiosonde data of local stations. By combining the water vapor pressure proﬁle retrieving from single ground-based GPS and temperature and pressure proﬁle from reference model, the refractivity proﬁle can be obtained, which is compared to radiosonde measurements. The comparison results show that results of the proposed method are consistent with the results of radiosonde. By using over ten years’ (through 2008 to 2017) historical radiosonde meteorological data of different months at China Big-Triangle Points, i.e., Qingdao, Sanya, Kashi, and Jiamusi radiosonde stations, tropospheric radio refractivity proﬁles are retrieved and modeled. The comparison results present that the accuracies of refractivity proﬁle of the proposed method at Qingdao, Sanya, Kashi, and Jiamusi are about 5.48, 5.63, 3.58, and 3.78 N-unit, respectively, and the annual average relative RMSE of refractivity at these stations are about 1.66, 1.53, 1.49, and 1.23%, respectively.


Introduction
It is well-known that refraction and time delay occur when electromagnetic wave propagates in the atmosphere as the changing properties of the atmosphere, the speed of propagation is altered, and the propagation path is also changed, which lead to a refraction error that severely disrupts the accuracy of radar detection, space tracking, satellite surveying and mapping, navigation, etc. The tasks of determining the profiles of the tropospheric refractive index and its characteristics are relevant for the both optical and radio ranges [1][2][3][4]. Tropospheric refractivity is a key parameter to evaluate the refraction error. As its variation in height is much greater than the horizontal direction, the accounting of atmospheric refractivity with different height is required so as to improve target tracking and navigation accuracy [5][6][7][8]. For this reason, some different atmospheric profile models were investigated based on a large number of data and theory. For example, Hopfield model [9] and segmented model [10], now as relatively accurate atmospheric prediction models, were developed by the atmospheric statics equation and ground meteorological data. While the actual atmosphere is not the same as the statics equation, the wet refractivity of those models lacks theoretical basis, and as a result, the results of the predicted models have a significant error with the actual atmosphere. Water vapor is a passive admixture and its dynamics depend on air temperature and wind speed fields, although water vapor is a small proportion in the whole atmosphere and has drastic changes with space and time, which makes it harder to monitor its distribution accurately [11][12][13].
There are also other methods to obtain atmospheric profile, such as radiosonde, Global Navigation Satellite System Meteorology (GNSS-MET), and water vapor radiometer. Radiosonde and water vapor radiometer are seen as high precision methods, while they suffer a higher cost and a lower temporal and spatial resolution. GNSS-MET is a potential method that has less cost than radiosonde, but higher precision than a traditional empirical model, moreover, it can be used in real time and has higher temporal and spatial resolution [14][15][16][17][18][19].
For ground-based GNSS remote sensing of the tropospheric refractivity profile, forward modeling is a process of calculating the refractivity profile based on the ray tracing algorithm, while inversion is a process of retrieving the refractivity profile by suitable algorithm based on tropospheric delay obtained by GNSS solution. Those inversion algorithms can be divided into two categories: one is to establish a distribution model of atmospheric parameters as a function of altitude, then solve model parameters using polynomial fitting, which is suitable for those atmospheric parameters changes regularly in profile distribution. The other is to build linear or non-linear regression statistics models of atmospheric parameters using a large amount of historical data, which is suitable for those atmospheric parameter changes in a random way and amplitude of fluctuations is more complicated and intense.
In this investigation, firstly, we retrieved the profile of water vapor pressure by nonlinear regression statistics models using GPS measurement based on Relevance Vector Machine [20]. Secondly, atmospheric temperature profile was calculated based on least square curve fitting principle. Thirdly, atmospheric pressure profile was computed according to the gradient of temperature profile. At last, atmospheric refractivity profile was computed by the profiles of atmospheric temperature, pressure, and water vapor pressure. A good consistency was shown between the retrieval refractivity and radiosonde refractivity. Furthermore, we can obtain the profile of atmospheric temperature, pressure, and water vapor pressure at the same time, which makes it flexible to use.

Data
The datasets used in this study include the historical meteorological data and GPS observation data, which were separated into two parts. The first part is the radiosonde observations from 2008 to 2016 to construct the model, and the second part is the radiosonde and GPS observations from the year 2017 to validate the proposed model. Our Figure 1. The triangle is composed of these stations representing some typical climates of China. The space measurement and control stations at the Big-Triangle points have been playing a vital role in aeronautics and astronautics, therefore investigation of the accurate tropospheric refractivity estimation is very important for improving the accuracy of space measurement. Radiosonde meteorological data including pressure, temperature, humidity, and wind were measured from the surface to the lower stratosphere by balloon-borne radiosonde, which was released twice a day at 8:00(LT) and 20:00(LT). The data were sorted by month and the sampling interval R Δ in the height of these data is shown as follows: The location of radiosonde stations and GPS observations sites of Qingdao, Sanya, Kashi, and Jiamusi are given in Table 1.

Method
In this section, we provide the method used to retrieve profile of water vapor pressure based on GPS measurement, as well as a method to fit profile of temperature, pressure, and some useful parameters. The radiosonde data of 10 years through 2008 to 2017 was conducted to build the inversion model, and the radiosonde data and GPS observation data of 2017 were used to validate the accuracy of these model. Radiosonde meteorological data including pressure, temperature, humidity, and wind were measured from the surface to the lower stratosphere by balloon-borne radiosonde, which was released twice a day at 8:00(LT) and 20:00(LT). The data were sorted by month and the sampling interval ∆R in the height of these data is shown as follows: The location of radiosonde stations and GPS observations sites of Qingdao, Sanya, Kashi, and Jiamusi are given in Table 1.

Method
In this section, we provide the method used to retrieve profile of water vapor pressure based on GPS measurement, as well as a method to fit profile of temperature, pressure, and some useful parameters. The radiosonde data of 10 years through 2008 to 2017 was conducted to build the inversion model, and the radiosonde data and GPS observation data of 2017 were used to validate the accuracy of these model.

Retrieving the Profile of Water Vapor Pressure Based on GPS Measurement
Tropospheric atmosphere is normally considered as a non-dispersive medium, and its influence on radiowave propagation can be expressed by refractive index n or refractivity N. In the neutral atmosphere of the Earth, n is very close to one, thus it is more convenient to use the so-called refractivity instead in engineering. The refractivity N (in "N-unit") is related to the refractive index by: The N-unit is a dimensionless unit, in terms of which, refractivity is expressed as atmosphere exceeds unity. Refractivity is a function of atmospheric temperature, pressure, and humidity, which has significant large-scale variations, and can be expressed by Equation (3) [22]: where and: where k 1 , k 2 , k 3 are empirical constants depending on frequency, the changes of ki brought by the frequency of 0-60 GHz are almost negligible [23], R is the universal gas constant, ρ is the density of dry air, M d and M w are the molar mass for dry air and water vapor, respectively, Z w is compressibility factors for water vapor, T is temperature, and p w is water vapor pressure. N h denotes hydrostatic term of refractivity only on the total density of air, N w denotes wet term of refractivity, which depends only on water vapor pressure and the temperature. While the hydrostatic part is larger than the wet part, the wet refractivity is much more variable and difficult to model. The uncertainty of refractivity is mainly caused by wet term, which accounts for about 10-30% of the entire refractive index.
The zenith tropospheric delay (ZTD) is the integrated refractivity along a vertical path through the neutral atmosphere: where c is the speed of light in a vacuum, τ is the delay measured in units of time, and N is the neutral atmospheric refractivity. ZTD is the sum of the hydrostatic or 'dry' delay (ZHD) and non-hydrostatic or zenith 'wet' delay (ZWD), due to the effects of dry gases and water vapor, respectively. In addition, the ZHD is proportional to the atmospheric pressure while the pressure is mainly related to height, and therefore, the ZHD is almost constant and the seasonal variations of zenith tropospheric delay (ZTD) are due primarily to the ZWD. Previous research showed that the accuracy of Saastamoinen model of ZHD caused by hydrostatic term of refractivity reach sub-millimeter level, while that of ZWD caused by wet term of refractivity is about 3 cm [24].

Correlation Analysis between Zenith Wet Delay (ZWD) and Water Vapor Pressure
As ZWD is highly variable due possibly to varying climate, relating to the temperature and water vapor, we analyzed the correlation between ZWD and partial pressure of water vapor of Qingdao radiosonde meteorological data at different heights. The wet component ZWD was obtained by subtracting ZHD from ZTD. The ZTD was computed by Equation (6) and the ZHD can be well calculated from surface meteorological data given by [25]: where Ps is the surface pressure in hPa, which can be measured by pressure sensor, f (λ, H) is a factor for correcting the local gravity as f (λ, H) = 1 − 0.00266 cos(2λ) − 0.28H, where λ is the latitude and H is the height above the ellipsoid in meters. Water vapor pressure E w can be calculated as a function of the relative humidity as follows: where Rh is the relative humidity and T is the temperature. Figure 2 indicates the statistical relationships between ZWD and water vapor pressure of Qingdao station at the height of 500, 2000, 5000, and 10,000 m, respectively. A well linear correlation was found between ZWD and water vapor pressure according to the figure, a stronger correlation was found at the lower height, and the correlation decreases with increase of height, this is due to the fact that most water vapor exists in the lower atmosphere and the water vapor density as well as water vapor pressure decrease with height.
As ZWD is highly variable due possibly to varying climate, relating to the temperature and water vapor, we analyzed the correlation between ZWD and partial pressure of water vapor of Qingdao radiosonde meteorological data at different heights. The wet component ZWD was obtained by subtracting ZHD from ZTD. The ZTD was computed by Equation (6) and the ZHD can be well calculated from surface meteorological data given by [25]: where Ps is the surface pressure in hPa, which can be measured by pressure sensor, where Rh is the relative humidity and T is the temperature. Figure 2 indicates the statistical relationships between ZWD and water vapor pressure of Qingdao station at the height of 500, 2000, 5000, and 10,000 m, respectively. A well linear correlation was found between ZWD and water vapor pressure according to the figure, a stronger correlation was found at the lower height, and the correlation decreases with increase of height, this is due to the fact that most water vapor exists in the lower atmosphere and the water vapor density as well as water vapor pressure decrease with height.  To investigate the variation characteristic of the correlation between zenith wet delay and water vapor pressure further, over 10 years' (through 2008 to 2017) radiosonde meteorological data of China Big-Triangle Points, Qingdao, Sanya, Kashi, and Jiamusi stations with the different months were processed and analyzed, respectively. Figure 3 depicts the variations of the correlation coefficients with height between zenith wet delay and water vapor pressure at different months of Qingdao, Sanya, Kashi, and Jiamusi stations, respectively. Tables A1-A4 in the Appendix A present the correlation coefficients at typical heights from January to December of those stations. We can find Atmosphere 2021, 12, 1468 6 of 18 that the similar variation trend was shown in Figure 3, and there is strong correlation for heights less than 5000 km and almost all the correlation coefficients are greater than 0.6. lay and water vapor pressure further, over 10 years' (through 2008 to 2017) radiosonde meteorological data of China Big-Triangle Points, Qingdao, Sanya, Kashi, and Jiamusi stations with the different months were processed and analyzed, respectively. Figure 3 depicts the variations of the correlation coefficients with height between zenith wet delay and water vapor pressure at different months of Qingdao, Sanya, Kashi, and Jiamusi stations, respectively. Tables A1-A4 in the Appendix A present the correlation coefficients at typical heights from January to December of those stations. We can find that the similar variation trend was shown in Figure 3, and there is strong correlation for heights less than 5000 km and almost all the correlation coefficients are greater than 0.6.

Estimation of Zenith Wet Delay (ZWD) by GNSS Precise Point Positioning
Due to the atmospheric refraction, GNSS signals propagating through the Earth's atmosphere are slightly lengthened with retarded speeds, and the delays of GNSS signals are called the neural atmospheric delay and ionospheric delay in the ionosphere, which were considered as a nuisance and error source for a long time, while nowadays these delays could be extracted and transformed into the useful atmospheric parameters. After being passed through the atmosphere into the ground, GNSS signals carry a large amount of atmospheric information. Therefore, GNSS precise point positioning (PPP) [26] can be used to detect tropospheric information as well as positioning. Here, a ground-based GNSS method is used to calculate the ZWD as follows: Firstly, ZTD is computed by GPS zero-differenced observation combinations, given by Equation (9):

Estimation of Zenith Wet Delay (ZWD) by GNSS Precise Point Positioning
Due to the atmospheric refraction, GNSS signals propagating through the Earth's atmosphere are slightly lengthened with retarded speeds, and the delays of GNSS signals are called the neural atmospheric delay and ionospheric delay in the ionosphere, which were considered as a nuisance and error source for a long time, while nowadays these delays could be extracted and transformed into the useful atmospheric parameters. After being passed through the atmosphere into the ground, GNSS signals carry a large amount of atmospheric information. Therefore, GNSS precise point positioning (PPP) [26] can be used to detect tropospheric information as well as positioning. Here, a ground-based GNSS method is used to calculate the ZWD as follows: Firstly, ZTD is computed by GPS zero-differenced observation combinations, given by Equation (9): where: P: code pseudorange measurements of two frequencies (m); ϕ: phase pseudorange measurements of two frequencies (L1, L2), respectively (m); ρ: the distance between satellite and antenna of GPS receiver (m); c: the speed of light (m/s); dt: the receiver clock error (s); dT: the satellite clock error (s); T trop : the signal delay caused by atmosphere (m); T iono : the signal delay caused by ionosphere (m); N: the integer ambiguity (cycle); λ: the wavelength of carrier phase (m); ε P : is the code error containing remaining errors, such as multipath and noise (m); ε ϕ : is the phase error containing remaining errors, such as multipath and noise (m). As T trop is the slant tropospheric delay along the signal path from the satellite to the receiver, each satellite with a different elevation angle corresponds to a different slant tropospheric delay (STD), a mapping relationship given by Equation (10): where M is the mapping function of tropospheric delays, the common mapping functions include Herring Mapping Function [27], Niell Mapping Function [28], Global Mapping Function [29], Vienna Mapping Function [30], etc. Herein, Global Mapping Function is recommended for use as it provides better precision than the Niell Mapping Function and smaller height biases with respect to Vienna Mapping Function, and it can be easy to implement because it only needs input parameters of station coordinates and day of year [29]. Secondly, ZHD is calculated by an empirical model based on ground surface pressure given by Equation (7).
Lastly, the ZWD can be obtained by subtracting ZHD from ZTD.

Building Inversion Model by Intelligent Optimization Algorithm
For the inversion process of the profile of water vapor pressure, first of all, forward calculation of input and output parameters should be carried out, that is, the tropospheric path delay ZWD should be calculated by ray-tracing algorithm [31,32] based on the known atmospheric parameter profile: ZWD = r T r0 (n w − 1)cscθdr = 10 −6 r T r0 N w cscθdr = 10 −6 r T r0 (3.73 × 10 5 e w T 2 )cscθdr (11) where θ is elevation of the signal path, it denotes zenith wet delay when its value equals to 90 degree. n w and N w are the wet term of atmospheric refractive index and refractivity respectively, e w and T are water vapor pressure and temperature, respectively. Then the inversion network can be built based on a large number of historical meteorological data. As there is a well linear correlation between water vapor pressure and zenith wet delay, one can build the relationship between the profile of water vapor pressure and ZWD. The inversion network should be trained by using a large number of historical meteorological data based on an intelligent optimization algorithm such as Neural Network [33], Support Vector Machine [34], Relevance Vector Machine (RVM) [21], etc. Here, we chose RVM to train the inversion network, because it has a number of advantages, which include the benefits of probabilistic predictions, automatic estimation of parameters, and the facility to utilize arbitrary basis functions, which are not necessary to be 'Mercer' kernels. The key feature of RVM is that as well as offering good generalization performance, the inferred predictors are exceedingly sparse in that they contain relatively few non-zero w i parameters. The majority of parameters are automatically set to zero during the learning process, giving a procedure that is extremely effective at discerning those basis functions that are 'relevant' for making good predictions.
The training of RVM is conducted under the Bayesian framework, which can be used for regression estimation and prediction to obtain the distribution of predicted values, given by: where K(x, x i ) is the function of kernels, w i is the weight coefficients of the model, the output is a linearly-weighted sum of N, generally nonlinear and fixed, basis functions K(x, x i ). With a wide domain of convergence and its applications to a variety of situations, Gaussian kernel function is usually used in RVM algorithm, which can be expressed by the following equation: where σ is width parameter, · denotes the norm of vector. Here, σ is taken as the standard deviation of the input training data vector, the input vector X is given by: where T, P, E w are the ground surface atmospheric temperature, pressure and water vapor pressure of historical meteorological data, respectively, ZWD is calculated by Equation (11). Normalize all data according to size range, and take standard deviation of normalized vector inner product: The inversion network training procedure is shown in Figure 4 and can be built by these steps as follows: 1.
The radiosonde data including atmospheric temperature, relative humidity, pressure, and height should be extracted from years of original radiosonde data. Then, we can obtain the water vapor pressure E w by Equation (8) and calculate ZWD by Equation (11).

2.
Divide all the data mentioned in step 1 into two parts: the train-data sets and the test-data sets, of which the proportion is 90 and 10 percent, respectively. 3.
Train the inversion network of RVM based on the train-data sets by adjusting the key parameters such as function of kernels, weight coefficients, and normalization variables.

4.
By putting surface temperature, water vapor pressure, pressure, and ZWD of testdata sets into inversion network built in step 3, we can get the retrieved profile of water vapor pressure. The accuracy of the model can be evaluated by comparing the retrieved profile of water vapor pressure with radiosonde profile of water vapor pressure. An accuracy threshold is set as 7 N-unit here, if the accuracy is enough, the modeling procedure is over, otherwise step 3 will be repeated until the error expectation is reached.

Retrieving the Profile of Water Vapor Pressure Based on GNSS Measurement
Once the inversion model has been built, water vapor pressure can be retrieved in real time based on original measurement of single ground-based GNSS and meteorological sensor. The flow diagram for retrieving tropospheric refractivity is shown in Figure 5. Firstly, the original GNSS observed data and broadcast ephemeris are collected by GNSS receiver, and the surface meteorological data is collected by meteorological sensor. Secondly, measurement data based on precise point positioning algorithm introduced in Section 3.1.2. are processed. As a result, the ZWD will be obtained. Thirdly, by inputting the ZWD, T, E w , and P at the surface to the inversion model, we can get the profile of water vapor pressure.

Data Preparation
Train-data sets (T, P, Ew at different height and ZWD) Test-data sets (T, P, Ew at different height and ZWD)

Retrieving the Profile of Water Vapor Pressure Based on GNSS Measurement
Once the inversion model has been built, water vapor pressure can be retrieved in real time based on original measurement of single ground-based GNSS and meteorological sensor. The flow diagram for retrieving tropospheric refractivity is shown in Figure 5. Firstly, the original GNSS observed data and broadcast ephemeris are collected by GNSS receiver, and the surface meteorological data is collected by meteorological sensor. Secondly, measurement data based on precise point positioning algorithm introduced in Section 3.1.2. are processed. As a result, the ZWD will be obtained. Thirdly, by inputting the ZWD, T, Ew, and P at the surface to the inversion model, we can get the profile of water vapor pressure.

Inversion Examples and Analysis
The inversion network of Sanya was built based on the method introduced in Section 3.1.3 and 10 years' historical meteorological data from 2008 to 2017 of Sanya radiosonde station was conducted as modeling data. The GPS observation data of 2017 at

Inversion Examples and Analysis
The inversion network of Sanya was built based on the method introduced in Section 3.1.3 and 10 years' historical meteorological data from 2008 to 2017 of Sanya radiosonde station was conducted as modeling data. The GPS observation data of 2017 at Sanya station was used to retrieve the pressure of water vapor and compare to radiosonde results as well as ITU model. A typical comparison result of the profile of water vapor pressure between retrieved, ITU model, and radiosonde is shown in Figure 6a.

Inversion Examples and Analysis
The inversion network of Sanya was built based on the method introduced in Section 3.1.3 and 10 years' historical meteorological data from 2008 to 2017 of Sanya radiosonde station was conducted as modeling data. The GPS observation data of 2017 at Sanya station was used to retrieve the pressure of water vapor and compare to radiosonde results as well as ITU model. A typical comparison result of the profile of water vapor pressure between retrieved, ITU model, and radiosonde is shown in Figure 6a.  It can be seen that the radiosonde water vapor pressure varies irregularly along the height. Although the retrieved results have a significant error at some heights, the variation trend in the whole height layer is in good agreement with the radiosonde data. Statistics showed that the average water vapor content of Sanya station in the layers up to 3 and 5 km were 91.61 and 97.58%, respectively. Figure 6b shows the root mean square error between the radiosonde and the retrieved profile of water vapor pressure of those input data and the maximum error of water vapor pressure of retrieved and ITU model can reach about 2.5 and 3.7 hPa respectively, and the mean of RMSE of water vapor pressure of retrieved and ITU model over the whole height are about 1.1 and 1.9 hPa, respectively.

Fitted Temperature and Pressure Profile of Reference Model
The distribution of atmospheric temperature with altitude is more regular and stable than that of water vapor pressure. Different models of atmospheric temperature distribution are selected according to different geographical locations and seasons, the least square method is used to fit the coefficients of the polynomial model based on a large numbers of historical meteorological radiosonde data of local stations. The temperature does not change obviously with the season in low-latitude area, the distribution model of temperature is given by [35]: Atmosphere 2021, 12, 1468

of 18
In summer of mid-latitude and high-latitude area, the atmospheric temperature profile model of different heights is expressed by: In winter of mid-latitude area, the atmospheric temperature profile model of different heights is expressed by: In winter of high-latitude area, the atmospheric temperature profile model of different heights is expressed by: where T 0 is the temperature of ground surface (K); T 1 , T 2 , T 3 are temperature of height h 1 , h 2 , h 3 , respectively; h 1 , h 2 , h 3 , h 4 are top height of troposphere, isothermal layer, thermal inversion layer, and second isothermal layer, respectively; h is height (km); a, b, c, d, e are fitting coefficients of historical radiosonde data of different areas and months. According to the definition of the least square method, the polynomial coefficient was calculated by historical temperature radiosonde data of years so that the sum of the squares of deviations at each point was minimized, then the temperature model was built. A typical comparison result of temperature profile between fitted and radiosonde is shown in Figure 7.  It can be seen that the variation of temperature profile is more regular and stable than that of water vapor pressure, the variation trend of the fitted temperature profile shows a good agreement with the radiosonde data. In order to validate the accuracy of the fitted method, 50 test-data sets were input to the network, Figure 7b shows the root mean square error of fitted and ITU model temperature compared to radiosonde temperature, the biggest temperature error of fitted and ITU model can reach over 2.4 and 7.4 K respectively, the mean RMSE of temperature of fitted and ITU model over the It can be seen that the variation of temperature profile is more regular and stable than that of water vapor pressure, the variation trend of the fitted temperature profile shows a good agreement with the radiosonde data. In order to validate the accuracy of the fitted method, 50 test-data sets were input to the network, Figure 7b shows the root mean square error of fitted and ITU model temperature compared to radiosonde temperature, the biggest temperature error of fitted and ITU model can reach over 2.4 and 7.4 K respectively, the mean RMSE of temperature of fitted and ITU model over the whole height are about 1.7 and 4.3 K, respectively.
The variation of pressure P(h) with height is calculated by the following formula based on the temperature profile: where P i , T i , L i are the pressure (hPa), temperature (K), and gradient of temperature (K/km) of the ith height layer, respectively. Figure 8 shows the comparison between fitted, ITU model, and radiosonde atmospheric pressure profile and the root mean square error of fitted pressure profile and ITU model pressure profile. The distribution of pressure with height is most stable compared to water vapor pressure and temperature profile, the mean RMSE of pressure calculated by fitted and ITU model over the whole height are about 2.2 and 5.6 hPa, respectively.

Results and Discussion
The retrieval refractivity profile can be obtained by Equation (3) after the profile of water vapor pressure, temperature profile, and pressure profile is calculated. Figure 9 shows the comparison between fitted, ITU model, and radiosonde refractivity profile and the RMSE of retrieved refractivity profile and ITU model refractivity profile. The mean RMSE of refractivity calculated by retrieved and ITU model over the whole height are about 5.3 and 8.5 N-unit, respectively.

Results and Discussion
The retrieval refractivity profile can be obtained by Equation (3) after the profile of water vapor pressure, temperature profile, and pressure profile is calculated. Figure 9 shows the comparison between fitted, ITU model, and radiosonde refractivity profile and the RMSE of retrieved refractivity profile and ITU model refractivity profile. The mean RMSE of refractivity calculated by retrieved and ITU model over the whole height are about 5.3 and 8.5 N-unit, respectively.
To verify the effectiveness of the proposed method further, ten years' historical meteorological data from 2008 to 2017 of Qingdao, Sanya, Kashi, and Jiamusi radiosonde stations for different months were used to build the inversion model. GPS observations of 2017 at these four stations were conducted to retrieve the refractivity and validate the accuracy, as shown in Figure 10, the annual average error of the proposed method at Qingdao, Sanya, Kashi, and Jiamusi are about 5.48, 5.63, 3.58, and 3.78 N-unit, respectively.
The retrieval refractivity profile can be obtained by Equation (3) after the profile of water vapor pressure, temperature profile, and pressure profile is calculated. Figure 9 shows the comparison between fitted, ITU model, and radiosonde refractivity profile and the RMSE of retrieved refractivity profile and ITU model refractivity profile. The mean RMSE of refractivity calculated by retrieved and ITU model over the whole height are about 5.3 and 8.5 N-unit, respectively. To verify the effectiveness of the proposed method further, ten years' historical meteorological data from 2008 to 2017 of Qingdao, Sanya, Kashi, and Jiamusi radiosonde stations for different months were used to build the inversion model. GPS observations of 2017 at these four stations were conducted to retrieve the refractivity and validate the accuracy, as shown in Figure 10, the annual average error of the proposed method at  The RMSE of Qingdao and Sanya station are relatively larger than that of Kashi and Jiamusi stations, this is because they have different climatic peculiarities. Qingdao and Sanya stations are located in the north temperate monsoon coastal region and the low latitude coastal region with a tropical ocean monsoon climate, respectively, while Kashi and Jiamusi stations are located in a warm temperate continental arid climate zone and the middle temperate humid climate zone with a continental monsoon climate, respectively. As is shown by Figures 11 and 12, the ground surface refractivity is generally higher in the summertime than in the wintertime, the ground surface refractivity of Qingdao and Sanya station are relatively larger than that of Kashi and Jiamusi station. Figure 13 shows the relative RMSE of refractivity profile of these four stations mentioned above in different months, which ranges from 1.1 to 2%, and the annual average relative RMSE of refractivity of the proposed method at Qingdao, Sanya, Kashi, and Jiamusi are about 1.66, 1.53, 1.49, and 1.23%, respectively. The RMSE of Qingdao and Sanya station are relatively larger than that of Kashi and Jiamusi stations, this is because they have different climatic peculiarities. Qingdao and Sanya stations are located in the north temperate monsoon coastal region and the low latitude coastal region with a tropical ocean monsoon climate, respectively, while Kashi and Jiamusi stations are located in a warm temperate continental arid climate zone and the middle temperate humid climate zone with a continental monsoon climate, respectively. As is shown by Figures 11 and 12, the ground surface refractivity is generally higher in the summertime than in the wintertime, the ground surface refractivity of Qingdao and Sanya station are relatively larger than that of Kashi and Jiamusi station.    Figure 11. Boxplot of the ground refractivity in different months of 10 years' radiosonde data, the red marks represent outlier data of refractivity. Figure 12. The mean ground refractivity in different months of 10 years' radiosonde data. Figure 12. The mean ground refractivity in different months of 10 years' radiosonde data. Figure 13 shows the relative RMSE of refractivity profile of these four stations mentioned above in different months, which ranges from 1.1 to 2%, and the annual average relative RMSE of refractivity of the proposed method at Qingdao, Sanya, Kashi, and Jiamusi are about 1.66, 1.53, 1.49, and 1.23%, respectively.

Conclusions
In this study, we proposed a method for retrieving low tropospheric refractivity profile from GPS measurement. While the hydrostatic part is larger than the wet part, the wet refractivity is much more variable and difficult to estimate, the variability of atmospheric refractivity is mainly attributed to the wet part, which accounts for about 10-30% of the entire refractive index. In previous studies and recommendations, this term is estimated from an empirical formula, which induced large errors. By analyzing the variations of the correlation between zenith wet delay and water vapor pressure based on over ten years' radiosonde meteorological data at Qingdao, Sanya, Kashi, and Jiamusi stations with different months, we found that there is a good linear correlation between zenith wet delay and water vapor pressure at different heights. The correlation decreases with increase of height due to the fact that water vapor mainly exists in the lower troposphere and the water vapor density as well as water vapor pressure rapidly decreasing with height. So we can build the relationship between the profile of water vapor pressure and ZWD, estimation of zenith tropospheric delay from GPS measurement is capable of providing accurate ZWD. We can retrieve water vapor pressure profile based on ZWD and the inversion model. By integrating the temperature and pressure profile adopted by reference atmosphere, we compared our results of atmospheric refractivity with the radiosonde measurement at China Big-Triangle Points, Qingdao, Sanya, Kashi, and Jiamusi station with a ten-year dataset. The RMSE of refractivity of the proposed method at Qingdao, Sanya, Kashi, and Jiamusi are about 5.48, 5.63, 3.58, and 3.78 N-Unit, respectively, and the annual average relative RMSE of refractivity of the proposed method at these stations are about 1.66, 1.53, 1.49, and 1.23%, respectively. This proposed method is very promising for radio meteorology research and atmospheric radio wave propagation correction for aerospace measurement and control.

Conclusions
In this study, we proposed a method for retrieving low tropospheric refractivity profile from GPS measurement. While the hydrostatic part is larger than the wet part, the wet refractivity is much more variable and difficult to estimate, the variability of atmospheric refractivity is mainly attributed to the wet part, which accounts for about 10-30% of the entire refractive index. In previous studies and recommendations, this term is estimated from an empirical formula, which induced large errors. By analyzing the variations of the correlation between zenith wet delay and water vapor pressure based on over ten years' radiosonde meteorological data at Qingdao, Sanya, Kashi, and Jiamusi stations with different months, we found that there is a good linear correlation between zenith wet delay and water vapor pressure at different heights. The correlation decreases with increase of height due to the fact that water vapor mainly exists in the lower troposphere and the water vapor density as well as water vapor pressure rapidly decreasing with height. So we can build the relationship between the profile of water vapor pressure and ZWD, estimation of zenith tropospheric delay from GPS measurement is capable of providing accurate ZWD. We can retrieve water vapor pressure profile based on ZWD and the inversion model. By integrating the temperature and pressure profile adopted by reference atmosphere, we compared our results of atmospheric refractivity with the radiosonde measurement at China Big-Triangle Points, Qingdao, Sanya, Kashi, and Jiamusi station with a ten-year dataset. The RMSE of refractivity of the proposed method at Qingdao, Sanya, Kashi, and Jiamusi are about 5.48, 5.63, 3.58, and 3.78 N-Unit, respectively, and the annual average relative RMSE of refractivity of the proposed method at these stations are about 1.66, 1.53, 1.49, and 1.23%, respectively. This proposed method is very promising for radio meteorology research and atmospheric radio wave propagation correction for aerospace measurement and control.