A Novel Global Grid Model for Atmospheric Weighted Mean Temperature in Real-Time GNSS Precipitable Water Vapor Sounding

The atmospheric weighted mean temperature (Tm) is an important parameter in calculating the precipitable water vapor from Global Navigation Satellite System (GNSS) signals. As both GNSS positioning and GNSS precipitable water vapor detection require high spatial and temporal resolutions for calculating Tm, high-precision modeling of Tm has gained widespread attention in recent years. The previous models for calculating Tm have the limitation of too many model parameters or single-grid data. Therefore, this study presents a global high-precision Tm model (GGTm-H model) developed from the latest Modern-Era Retrospective Analysis for Research and Applications, version-2 (MERRA-2) atmospheric reanalysis data provided by the United States National Aeronautics and Space Administration. The accuracy of the GGTm-H model was verified by combining the MERRA-2 surface Tm data and 319 radiosonde data. The results highlighted that 1) When the MERRA-2 Tm data were used as a reference value, the mean annual RMSE of the GGTm-H model was observed to be 2.72 K. When compared with the Bevis model, GPT2w-5 model, and GPT2w-1 model, the GGTm-H model showed an improvement of 1.5, 0.33, and 0.21 K, respectively. 2) When the radiosonde data were used as a reference value, the mean bias and RMSE of the GGTm-H model were −0.41 K and 3.82 K, respectively. Compared with the other models, the GGTm-H model had the lowest mean annual bias and RMSE. The developed model does not consider any meteorological parameters while calculating Tm. Therefore, it has important applications in the real-time and high-precision monitoring of precipitable water vapor from GNSS signals.


I. INTRODUCTION
W ATER vapor makes up a small part of the atmosphere, but it plays a key role in a range of weather phenomena, from heavy rains to droughts [1], [2]. Precipitable water vapor (PWV) refers to the precipitation generated by the total amount of PWV contained in the large air column of a unit section from the ground to the top of the atmosphere, all of which condenses to the ground, which is an important parameter reflecting the change of atmospheric PWV [3], [4]. With the development of Global Navigation Satellite System (GNSS) technology, ground-based GNSS inversion can obtain PWV with higher temporal and spatial resolution, which has become one of the important means for PWV monitoring [5], [6], [7], [8], [9], [10], [11], [12]. Recently, high-precision GNSS PWV have played an important role in the monitoring and analysis of extreme weather such as heavy rain and typhoons [13], [14], [15].
The atmospheric weighted mean temperature (Tm) is a key parameter for calculating GNSS vapor conversion coefficient [16], [17], [18], [19], [20]. High-precision and real-time Tm values can be obtained by integrating the atmospheric reanalysis data; however, the data update is slow with limited horizontal resolution. This significantly restricts the development and application of GNSS meteorology in global regions [21], [22], [23], [24]. In recent years, several researchers have developed many regional or global Tm models to satisfy the needs of PWV inversion. The well-known Bevis model was originally developed for a specific area (27°N-65°N) in the northern hemisphere [25]. These models show good accuracy in a specific area but are not suitable for GNSS-based PWV estimation over a large area. Zhang et al. [26] observed that lapse rate is an important parameter for Tm elevation correction and developed an enhanced Tm model for China (GH-Tm). Studies by Huang et al. [27], [28] indicated that large topographic fluctuations in China lead to significant systematic error in the GPT2w-1 model [29] valuation. Therefore, the IGPT2W model was developed based on the lapse rate. Compared with the GPT2w model, the improved model has enhanced accuracy in high-altitude areas. Yao et al. [30], [31], [32] used spherical harmonics to develop a globally applicable empirical Tm model, GWMT GTm-III. This model considers the height and the periodicity of Tm. Yao et al. [33] considered the approximately linear relationship between the atmospheric weighted temperature and temperature and used the globally applicable ECMWF to establish the Tm model. The large-scale data involved in modeling and the very large selected grid range affect the stability of the model to some extent at a global level. Based on the Bevis model, Jiang et al. [34] considered the temporal variations in Ts-Tm for each Ts-Tm grid data and developed a global grid Ts-Tm model. Similarly, Sun et al. [35] developed a global Tm model considering nonmeteorological parameters (GTrop model). Although this model has good accuracy, too many model parameters make the model complicated and limit its practical applications to some extent. Yao et al. [36] used the ECMWF products to analyze the Tm distribution characteristics in the vertical direction and developed a new global Tm model. Tested with the ECMWF and the radiosonde (RS) data, the RMS of the model is 3.84 and 4.36 K, respectively, achieving accuracy improvements of 27% and 20% compared to the existing model. Sun et al. [37] developed a new global grid-based Tm empirical model named GGNTm. A three-order polynomial function was utilized to fit the vertical nonlinear variation in Tm at the grid points. The model made significant improvements at high-altitude pressure levels. Cao et al. [38] proposed to correct the bias of the GRAPES_MESO forecasting data using a linear model and a spherical cap harmonic model. Compared to the existing empirical models that only capture the tidal variations, the CTropGrid products capture well the nontidal variations of Tm.
Although the above-mentioned Tm empirical models show their advantages, some limitations still exist, such as the adoption of only single gridded data for modeling, and the model parameters need to be further optimized. Therefore, it is of great importance to develop a novel global Tm grid model for real-time GNSS PWV sounding. A sliding window algorithm is introduced to divide the global into regular windows of the same size, and the Modern-Era Retrospective Analysis for Research and Applications, version-2 (MERRA-2) data are adopted in the calculation of Tm information.
RS data are important meteorological observation data, with observations typically made only twice a day (UTC 0:00 and UTC 12:00). At present, there are more than 1500 comprehensive sounding stations. Each sounding station provides layered data on atmospheric pressure, temperature, relative humidity, and other meteorological parameters and surface parameters such as atmospheric PWV from the surface to the near-earth space extending about 30 km above [39]. Since the sounding data are obtained through actual measurements, they have reliable accuracy and, hence, are widely used in the model construction and accuracy verification analysis of key troposphere parameters, such as Tm.

B. Methods
Since the MERRA-2 reanalysis data with the horizontal resolution of 0.5°× 0.625°(latitude × longitude) are used for modeling, the global region was first divided to obtain a grid consistent with the horizontal resolution of the MERRA-2 reanalysis data. The sliding window algorithm [40], [41], [42] is applied to obtain the tropospheric model parameters in the window.
The integral method is used to calculate the Tm value within the stratified height-range MERRA-2 reanalysis data. The basic formula can be expressed as follows: where e denotes the water vapor pressure (hPa), T is the temperature (K), and e i is the average water vapor pressure of layer i of the atmosphere. T i is the average temperature of layer i of the atmosphere. ΔH i is the thickness of the i layer of the atmosphere. h top and h L are the top and bottom heights of the hierarchical data integration calculation, respectively.
As shown in Fig. 1, Tm data from 545 RS stations in 2017 were used to evaluate the accuracy of Tm calculated from the MERRA-2 reanalysis data based on the annual mean bias and RMSE.
As depicted in Fig. 1, the overall bias of Tm calculated from the MERRA-2 reanalysis data was close to 0, and the absolute bias values of most measurement stations were within 1 K. The high-latitude stations had the maximum positive bias and are centrally distributed in Asia. A small negative bias was observed in the equatorial and low-latitude regions, whereas Antarctic stations had the largest negative bias. Overall, the Tm calculated from the MERRA-2 reanalysis data showed a narrow global distribution, indicating that the Tm obtained through MERRA-2 data integration has reliable accuracy. The Tm calculated from the MERRA-2 reanalysis data showed the largest RMSE value in Greenland and Asia. This may be due to the relatively skewed spatial and temporal distribution of Tm in the high latitudes. The Tm obtained from MERRA-2 data integration showed a small RMSE in other regions, implying that the Tm calculated from the MERRA-2 reanalysis data has high accuracy and can be used for the construction of the Tm lapse rate model and the weighted average temperature model. Due to the difference between the user's location and the elevation of the grids, large errors will be caused when used directly. Therefore, it is important to assess the spatiotemporal characteristics of the Tm lapse rate on a global scale. The Tm stratified data of the global region from 2014 to 2017 were calculated from the MERRA-2 reanalysis data at 0.5°× 0.625°a nd then analyzed by the changes in Tm values with height. The reanalysis data at three resolutions (89.5°N × 179.375°W, 89°N × 178.75°W, and 88°N × 177.5°W) dated January 1 and  the window center grid data were selected, and the changes in Tm with height were observed. As shown in Fig. 2, the distribution in MERRA-2 layered Tm and elevation were almost linear. Tm showed a declining trend with increased elevation, and its slope is the lapse rate of Tm, which can be expressed as follows: where γ denotes the lapse rate of Tm, h is the height, and k is the constant. To obtain the lapse rate of Tm, the Tm profile information at 9/25/81 lattice outlets and the corresponding bit potential height data were fitted to each window at the three global horizontal resolutions. Four representative lattice outlets were selected to assess the temporal variations in the lapse rate coefficient. The Tm lapse rates from 2014 to 2017 were analyzed to obtain the periodic characteristics of the Tm lapse rate.
As observed in Fig. 3, the lapse rate of Tm showed a significant periodic change, which was consistent with the vertically decreasing rate of Tm, as observed in previous studies [43], [44].
The Tm stratified data were calculated by the integration of the hierarchical meteorological parameters at 0.5°× 0.625°obtained from the MERRA-2 reanalysis data. Furthermore, the Tm lapse rate was subsequently obtained by analyzing the changes in the Tm values with height. As shown in Fig. 4, the color bar shows that the lapse rate of Tm had the largest absolute value near the equator, indicating drastic changes in Tm near the equator with height. The absolute value of Tm at high latitudes was small, indicating a relatively small change in Tm with height.

A. Construction of the Global Tm Lapse Rate Model
The lapse rate of Tm can effectively improve the vertical accuracy of Tm at different heights. Therefore, evaluating the lapse rate of T m can provide valuable reference for developing a high-precision T m model. As the sliding window algorithm effectively improves the application efficiency of the reanalysis data, the limitations of using a single-lattice-dot data modeling are avoided. The optimization of the T m lapse rate model application in spatial interpolation and the improvement in the height-corrected efficiency of the T m lapse rate model were obtained by linearly fitting the T m stratified data and bit-high potential data of MRERRA-2 from nine grid outlets in each window to obtain the T m lapse rate in that window. Furthermore, the analysis of spatial and temporal changes in the Tm lapse rate highlighted that the lapse rate of T m was mainly manifested in the annual and semiannual cycle characteristics. The model can  be expressed as follows: where γ denotes the lapse rate of Tm; h a and h b denote the target height and reference height, respectively; T ma and T mb denote the Tm value of the target height and reference height, respectively. The lapse rate of T m under each window can be expressed separately using the following equation: where doy represents the day of the year, a 0 represents the annual mean of the Tm lapse rate, and (a 1 , a 2 ) and (a 3 , a 4 ) indicate the annual and semiannual cycle coefficients of the Tm lapse rate, respectively. In each window globally, the coefficients are estimated based on the Tm layered profile information at the nine MERRA-2 grid points contained in the window at a resolution of 6 h through least-squares adjustment. The five vertical reduction rates of Tm at the center point of each window can be calculated globally, and the coefficient factors can be stored in the form of lattice points in the geometric center of each window. The effect of the Tm lapse rate model on Tm vertical correction was evaluated using Tm stratified data at three resolutions. The corresponding bit-high potential data were then selected to fit the Tm lapse rate model. According to the sliding window algorithm, the global region was divided into three horizontal resolutions of 1°× 1.25°, 2°× 2.5°, and 4°× 5°, and the Tm data and the corresponding bit-high potential data in the global region from 2014 to 2017 were used to develop the Tm lapse rate models, namely, the Tm-H, Tm-H2, and Tm-H4 models, respectively. Considering that Tm vertical correction is mainly applied to the lower layers of the atmosphere, the height range from the surface to 10 km was chosen for the stratified data while fitting the Tm lapse rate. Fig. 5 represents the annual average of the Tm lapse rate coefficient, the annual periodic amplitude, and the semiannual amplitude distributions at the three horizontal resolutions. As shown in Fig. 5, the differences between the annual average, annual periodic amplitude, and semiannual periodic amplitude distributions of the Tm-H and Tm-H2 model coefficients were not evident. The Tm lapse rate was calculated from the reanalysis data at 2°× 2.5°horizontal resolution to reduce the usage parameters of the model. Moreover, the Tm-H4 model was evaluated at a horizontal resolution of 4°× 5°due to the large window size. The annual periodic signal and semiannual signal intensities were weakened while using the Tm profile information from 81 lattice outlets in the window, resulting in relatively small changes in the annual periodic amplitude and semiannual amplitude. The absolute annual averages of Tm-H and Tm-H2 models were large in low and middle latitudes compared to those in high latitudes and bipolar regions. The lowest annual mean was observed in Greenland and Antarctica, indicating relatively drastic changes in Tm in the low and middle latitudes and relatively moderate changes in the high latitudes. Large annual and semiannual periodic amplitudes were observed in some parts of Europe, Eastern China, Central Africa, Antarctica, and North America due to the stronger seasonal variation in Tm at the higher latitudes. Therefore, the model coefficient amplitude showed large values. The annual periodic amplitude and semiannual periodic amplitude were relatively large, probably due to the polar phenomenon. However, the variations in the average annual periodic amplitude and semiannual periodic amplitude of the Tm-H4 model were relatively small, which may affect the application of the Tm lapse rate for elevation correction.
1) Accuracy Test of the Tm Lapse Rate Model: To evaluate the application of the Tm lapse rate models for elevation correction at three horizontal resolutions, the surface Tm grid products of the 2017 MERRA-2 reanalysis data were initially selected. Furthermore, the Tm lapse rate model output corresponding to the hierarchical Tm grid at three horizontal resolutions was corrected, and the MERRA-2 stratified Tm products were used to verify the Tm accuracy after elevation correction, as highlighted by the values of bias and RMSE in Table I and Fig. 6. As shown in Table I, the bias using the lapse rate function model at three resolutions is −0.7 K after applying surface MERRA-2 grid Tm products to correct the lapse rate model elevation corresponding to the hierarchical MERRA-2 grid Tm products. The RMSE of the Tm-H model was 0.01 K higher than that of the Tm-H2 and Tm-H4 models. The maximum and minimum distributions of bias and RMSE indicated that lower resolution led to better stability of the model with minimum error because of the accurate data of the Tm profile at the center point of the grid window. The lower resolution also led to increased profile information of the window at 9/25/81. The Tm value is the value of the window center point; hence, the resolution was observed to be low, and the Tm lapse rate coefficient was found to be more stable. As shown in Fig. 6, large overall differences between the bias and the RMSE obtained by the model elevation correction results of different resolution lapse rates were not observed; however, sea-land distribution and regional distribution differences were observed. In terms of bias, it showed a significant negative bias at the equator and near the middle and low latitudes, and a large positive bias was observed in northeast China, the northern United States, Greenland, and parts of Antarctica. This indicated that the Tm values corrected by the lapse rate model were slightly larger than the MERRA-2 stratified Tm values due to the strong seasonal correction of Tm. In terms of RMSE at three resolutions, the maximum value was observed at the equator, high-latitude areas, and sea-land junction because the Tm is affected by terrain and tropical monsoon, resulting in the ideal correction effect through extrapolation of the simple linear relationship. Overall, the decreasing rate function model at all three resolutions can provide high-precision Tm values. The accuracy of bias and RMSE in specific regions was significantly enhanced with the improved resolution of the model, implying that the improvement in the horizontal resolution of the model parameters can lead to improvement in the model accuracy. Therefore, the selection of the Tm lapse rate model can provide a reference value for the spatial interpolation and the modeling of Tm.
The application of the Tm lapse rate model in elevation correction was further validated by selecting the Tm data in 2017 as the reference value to test the model's accuracy. The surface grid Tm product of MERRA-2 was corrected to the probe height for greater accuracy. As shown in Table II, the surface MERRA-2 grid Tm data were corrected to the minimum bias and RMSE of the elevation of the sounding station, namely −0.09 K and 1.47 K, respectively, using the Tm-H model. The decreasing resolution of the Tm lapse rate model led to the calculation of a larger negative bias and a gradual decrease in stability. RMSE increased by 0.21 K (14%), i.e., from 1.47 K for Tm-H to 1.68 K for Tm-H2. Similarly, RMSE increased by 0.44 K (26%), from 1.68 K for Tm-H2 to 2.12 K for Tm-H4. However, with the further expansion of the horizontal resolution, comparing the model output with the measured Tm provided by the sounding station leads to greater errors. This study temporarily considered the application of the three resolved Tm lapse rate models for elevation correction to provide users with more stable Tm values. As shown in Fig. 7, the test results of the three resolution models corrected to the sounding height were similar. Overall, most of the selected sounding sites are on land, and the distribution of these sites in the ocean is relatively sparse. Moreover, the RMSE in the marine area was low, which may be due to the distribution of sea and terrain affecting the Tm lapse rate and further affecting the vertical correction effect. The values of bias and RMSE were larger in high latitudes of the northern hemisphere, which is due to the strong seasonal variation of Tm in high latitudes. These variations are difficult to model using only the lapse rate coefficients, thus leading to large errors. As shown in Fig. 8, considering the sounding station as the reference value and 1°× 1.25°as an example, the largest elevation of the sounding station led to a bias of about −0.5 K and an RMSE of about 2 K, indicating that the lapse rate coefficient led to a positive effect in the elevation correction. The bias and RMSE present a normal distribution form, indicating that the model lapse rate can be used for elevation correction and has good stability.

2) Precision Analysis of the Spatial Interpolation of GGOS Atmospheric Grid Tm Products:
As GGOS can provide highaccuracy grid Tm products, this study utilized the GGOS grid Tm product correction to test the accuracy of the vertical correction model. Since the GGOS grid Tm product has only a resolution of 2°× 2.5°, the impact of two horizontal resolutions, namely, 2°× 2.5°and 4°× 5°, on the model was considered to combine the grid resolution and lapse rate model resolution.
As shown in Tables II and III and Fig. 8, the test results of the GGOS grid Tm product and surface MERRA-2 grid Tm data were found to be similar. Overall, the bias of the GGOS grid Tm product was high relative to that of the surface MERRA-2 grid products, indicating the better stability of the surface MERRA-2  grid products. The 2°× 2.5°lapse rate model performed relatively better in terms of RMSE. Therefore, 2°× 2.5°was selected as the plane resolution of the global grid window. As shown in Fig. 9, bias and RMSE were higher at high altitudes, similar to the variations in these parameters observed for the MERRA-2 grid product data from the sounding stations. When the GGOS grid Tm elevation correction test was used at exploration sites near Antarctica, the RMSE was relatively high, indicating that the polar day and polar night may affect the lapse rate model. This is consistent with the observation that the lapse rate coefficient had a large annual/semiannual amplitude in the polar region.

B. Global Tm Grid Model Construction
Most of the existing Tm models use data modeling with single-grid data, which leads to several model parameters and low model efficiency. It is of utmost necessity to establish an efficient Tm model to carry out real-time PWV inversion. The model's efficiency should also be improved as much as possible while ensuring its accuracy. The sliding window algorithm greatly enhances the modeling data utilization by harnessing the Tm data within the joint window range. Moreover, it can effectively reduce the number of modeling parameters and maintain model accuracy. In this study, a global Tm refinement model was developed by adopting a sliding window algorithm based on the directly modeled in-window Tm data. The steps involved in the development of the global Tm refinement model are as follows. Initially, the world was regularly divided based on the sliding window technique to obtain each regular grid window. Since the modeling requires Tm data at the elevation of the window, the MERRA-2 surface Tm grid data in the window should be corrected to the average elevation of each window using the Tm-H1 model. The MERRA-2 surface Tm grid data from 2015 to 2017 for nine grid outlets of each window were corrected to the average elevation of the window. The average elevation of the window is the average value of the elevation of the MERRA-2 grid data in the window, as represented by the following:  IV  VALIDATION ACCURACY OF THE GGTM-H, GPT2W-5, GPT2W-1, AND BEVIS  MODELS BASED ON THE GLOBAL GRID VALUES OF MERRA-2 TM FOR  2015 (UNIT, K) where H 0 represents the average elevation of the window, H i M represents the ith MERRA-2 lattice elevation in the window, and n represents the number of MERRA-2 grids. Thus, the Tm refinement model was developed and represented as( follows: where Tm H 0 represents the Tm value of H 0 at the average elevation, and ΔTm represents the Tm value at the surface elevation corrected to the average elevation. Tm mean represents the annual mean at Tm H 0 . amp 1 and amp 2 represent the Tm annual and semiannual cycle amplitudes at the mean window elevation, respectively. ϕ 1 and ϕ 2 represent annual and semiannual initial phases, respectively. Each parameter of the formed global grid data in the Tm refinement model is stored separately as a file according to the grid plane resolution of 1°× 1.25°. Thus, the Tm model parameters of each window are globally defined, and the global Tm refinement model that can be applied to any height and can be used at high spatial and temporal resolution was finally developed. The global Tm refinement model is abbreviated as the GGTm-H model for convenience in this study.  Table IV. Table IV shows the distribution of bias and RMSE for different models, with an average bias of 0.21 K for the GGTm-H model, i.e., 1.34 K more than the average bias of the Bevis model. The average bias of GPT2w-5 and GPT2w-1 was 0.07 K and 0.06 K, respectively, possibly because the GPT2w models have high spatiotemporal resolution and a complex grid layout. In addition, all four models had a negative average bias, indicating that the Tm values calculated by the model were low compared to those obtained from MERRA-2 integration. The distribution of bias indicated that the GGTm-H model had the lowest bias values, and its stability was better than that of the other models. The distribution of RMSE indicated that the average RMSE of the GGTm-H model was 2.72 K, i.e., 1.5, 0.33, and 0.21 K more than the average bias of the Bevis, GPT2w-5, and GPT2w-1 models, respectively. This indicated that the developed Tm refinement model using the sliding window technique could estimate Tm with greater accuracy. The Bevis model had the highest RMSE value of 4.22 K, whereas the RMSE values of the GPT2w-5 and GPT2w-1 models were relatively low and stable, i.e., around 3 K.

1) GGTm-H Model Accuracy Test Using
As shown in Fig. 10, the Bevis model showed a large bias in parts of Antarctica, Greenland, western China, and Africa. As the Bevis model is a regional Tm model developed based on exploration data from North America, it shows a large bias when used to calculate the Tm values on a global scale. The GPT2w-5 model showed a large bias in some regions of western China, Africa, North and South America, and Antarctica. The distribution of bias of the GPT2w-1 model was more uniform than that of the GPT2w-5 model, indicating that the accuracy of Tm estimation can be improved by improving the model's resolution. However, the GPT2w-5 and GPT2w-1 models do not consider the elevation correction when calculating Tm. Therefore, the models' stability is not good in areas with large regional fluctuations. Compared with the other models, the GGTm-H model showed a uniform bias distribution when applied on a global scale. As the model considers the effect of the Tm lapse rate on the elevation correction, it can better estimate the Tm data on a global scale. To conclude, the GGTm-H model significantly outperforms other models on a global scale, with better stability and wide applicability. As shown in Fig. 11, the Bevis model showed the highest RMSE in parts of Antarctica, Greenland, western China, and eastern North America, indicating that a simple linear relationship does not accurately estimate the change in Tm on a global scale. The GPT2w-5 model showed high RMSE in parts of southwest China and western Africa, and the accuracy of the GPT2w-1 model was better than that of the GPT2w-5 model in those regions, indicating that the high-resolution Tm model has certain advantages while calculating Tm. As the GPT2w series  models do not consider Tm vertical correction for Tm estimation, a high RMSE was observed in these models compared to that in the GGTm-H model. The established GGTm-H model and GPT2w-1 model have comparable accuracy when applied on a global scale; however, the GGTm-H model demonstrated significant improvement in accuracy at high altitude areas, such as southwest China, southern North America, Africa, and parts of Antarctica, indicating that elevation correction can effectively improve the accuracy of Tm estimation in undulating areas.
For evaluating the performance of different models at different latitudes, the global regions were divided into 90°N-180°N and 180°S-90°S, and the variations in bias and RMSE of different models at different latitudes were calculated, as depicted in Fig. 12. The Bevis model showed a significant positive bias at high latitudes and a significant negative bias at 25°N-35°S. The Bevis model showed greater accuracy in the middle latitudes of the northern hemisphere because it mainly relies on the sounding data of this region for validation. Both the GPT2w-5 and GPT2w-1 models showed relatively low average bias at each latitude on a global scale. It is mainly because the GPT2w model has more computational parameters compared with the Bevis model, indicating that the model's performance is relatively stable. The GPT2w model had the highest calculation bias at high latitudes of the southern hemisphere because of the relatively complex spatial and temporal variations in Tm at high latitudes. When compared with the other models, the developed GGTm-H model    model had the lowest annual bias and RMSE, and the variation in error was small, indicating optimal accuracy and stability of the GGTm-H model. When compared with the GPT2w-5 model, the GPT2w-1 model had a low average bias and RMSE and small variations in bias and RMSE, indicating that the GPT2w-1 model performs better than the GPT2w-5 model. Moreover, the high spatial and temporal resolution in the GPT2w-1 model can estimate more accurate Tm values. Both the models have a negative average bias, indicating that the Tm values calculated by the models were less than those obtained from the sounding stations. Compared with the other models, the Bevis model had the highest RMSE. The model formula does not incorporate the periodic changes in Tm, which inevitably leads to uncertainty while calculating the Tm values. As shown in Fig. 13, the Bevis model had a more observable positive bias at higher latitudes and a more significant negative bias near the equator. The Bevis model was built based on the sounding data from North America, and therefore, it may have a large computational bias on a global scale. The GPT2w-5 and GPT2w-1 models showed significantly better accuracy than the Bevis model on a global scale. They showed a significant positive bias at partial GNSS stations near the equator, a significant negative bias near north latitude 30°N, and the highest negative bias in Antarctica. Compared with the other models, GGTm-H had a more uniform bias distribution on a global scale. In addition, the Tm values calculated by the model were found to be more accurate and more consistent with the Tm values obtained from the sounding data. Thus, this signifies the global applicability of the GGTm-H model. Fig. 14 shows that the Bevis model had the highest RMSE in Greenland and Antarctica, with significant RMSE at high latitudes and lower RMSE at the equator and low latitudes. The RMSE of the Bevis model in parts of North America was smaller than that of the GPT2w-5 model because the Bevis model is constructed from widely distributed sounding data in this region, indicating its applicability to the region. Similar to the Bevis model, the GPT2w-5, GPT2w-1, and GGTm-H models exhibited lower RMSE near the equator and relatively high RMSE at mid-to high-latitudes, implying significant spatiotemporal changes in Tm at high latitudes and thus affecting the model's accuracy. The global RMSE distribution map indicated that the GPT2w-1 model performs better at higher latitudes  Table V reveal that the GGTm-H model is a better model than the other global models and can calculate the Tm value at any position on a global scale with greater accuracy.
To analyze the accuracy of each model at different heights, 319 global sounding sites were arranged in order of their heights, and the bias and RMSE values at different heights were calculated. As shown in Table VI, the average bias of the Bevis, GPT2w-5, GPT2w-1, and GGTm-H models was 0.27, −0.85, −0.59, and −0.36K, respectively, and the average RMSE was 4.59, 4.07, 3.92, and 3.76, respectively, for heights less than 500 m. The variations in bias and RMSE indicate that the GGTm-H model outperformed all the other models and showed the best computational performance and stability among them. When compared with the Bevis model, the bias and RMSE values of the GGZTD-H model were more uniformly distributed. When compared with the GPT2w-5 and GPT2w-1 models, the GGTM-H model showed small variations in annual bias and low RMSE values. The accuracy of the GPT2w-5 and GPT2w-1 models improved to 0.31 and 0.16 K, respectively. The GPT2w-5 and GPT2w-1 models showed small bias and RMSE values at these heights with negative average annual bias, indicating that the Tm values calculated by the model are lower than those calculated by the sounding integral. The accuracy of GPT2w-1 was slightly Due to a small number of statistical sounding stations and contingency errors at heights less than 500 m, the Bevis model performed better than the GPT2w-5 and GPT2w-1 models. However, overall, the GPT2w-1 model performed better than the GPT2w-5 model, and the GGTm-H model showed the most stable model performance. These findings are consistent with the statistical results from the previous studies. Six RS sites with reliable data were selected to analyze the change in Tm values with time compared to the Tm values obtained by the sounding integral. As shown in Fig. 15, the Bevis model exhibited a relatively large deviation when compared with the other models. Both the GPT2w and GGTm-H models showed a good fit. The Bevis model showed poor performance at low latitudes, and the linear relationship between Tm and Ts was affected by latitude. The GPT2w model incorporates the seasonal characteristics of the troposphere parameters through a function of the annual cycle and semiannual cycle. Therefore, the GPT2w model output depicts noticeable periodic changes. However, the GPT2w model does not consider the elevation correction, resulting in a high bias. In addition, Tm shows seasonal variations, with a maximum in summer of nearly 290 K and a minimum in winter of nearly 260 K in low-and mediumdimensional areas. Compared with RS data, the deviation of the GGTm model is larger in winter and smaller in summer.

IV. CONCLUSION
The Tm values computed from the MERRA-2 atmospheric reanalysis data of 294 IGS sites and 545 RS stations have high accuracy and stability and can be used to construct a global multidimensional Tm refinement model. This study used the sliding window algorithm to evaluate the Tm lapse rate model by calculating global Tm lapse rates from the MERRA-2 reanalysis data at different resolutions and elevations. The result shows that the 2°× 2.5°Tm lapse rate model has good stability. On this basis, we developed a global Tm refinement model, GGTm-H, which can be corrected to any height using the Tm lapse rate model based on the sliding window technique. When compared with other Tm models, the GGTm-H model shows smaller deviation and higher stability over the globe. The validation of the GGTm-H model using MERRA-2 Tm data shows that the RMSE remains stable within 4 K. The validation of the GGTm-H model using RS data shows that the bias is less than −0.72 K and the RMSE remains stable within 4.5 K at the middle and high altitudes, which is far better than the Bevis and GPT models. For any three-dimensional coordinates and time, the Tm value at the average height of the window and the Tm value at the user position can be obtained using the elevation calculation of the Tm lapse rate model. The model is easy to use and can be easily compared to the values obtained by combining the MERRA-2 data and the radio-sounding data.
The sliding window algorithm is adopted in the construction of the GGTm model. Although the data usage efficiency of a single window is increased, the horizontal resolution of model parameters is reduced. How to combine multisource data to determine the global average elevation of each window and refine the model parameters of the tropospheric delay model in each window remains to be further studied.