Accuracy Evaluation and Analysis of GNSS Tropospheric Delay Inversion from Meteorological Reanalysis Data

Abstract

Accurate estimation of tropospheric delay is significant for global navigation satellite system’s (GNSS) high-precision navigation and positioning. However, due to the random and contingent changes in weather conditions and water vapor factors, the classical tropospheric delay model cannot accurately reflect changes in tropospheric delay. In recent years, with the development of meteorological observation/detection and numerical weather prediction (NWP) technology, the accuracy and resolution of meteorological reanalysis data have been effectively improved, providing a new solution for the inversion and modeling of regional or global tropospheric delays. Here, we evaluate the consistency and accuracy of three different types of reanalysis data (i.e., ERA5, MERRA2, and CRA40) used to invert the zenith tropospheric delay (ZTD) from 436 international GNSS service (IGS) stations in 2020, based on the integral method. The results show that the ZTD inversion of the three types of reanalysis data was consistent with the IGS ZTD, even in heavy rain conditions. Furthermore, the average precision of the ZTD inversion of the ERA5 reanalysis data was higher, where the mean deviation (bias), mean absolute error (MAE), and root mean square (RMS) were –3.39, 9.69, and 12.55 mm, respectively. The ZTD average precisions of the MERRA2 and CRA40 inversions were comparable, showing slightly worse performance than the ERA5. In addition, we further analyzed the global distribution characteristics of the ZTD errors inverted from the reanalysis data. The results show that ZTD errors inverted from the reanalysis data were highly correlated with station latitude and climate type, and they were mainly concentrated in the tropical climate zone at low latitudes. Compared to dividing error areas by latitude, dividing error areas by climatic category could better reflect the global distribution of errors and would also provide a data reference for the establishment of tropospheric delay models considering climate type.

1. Introduction

GNSS tropospheric delay refers to the signal delay and path bending generated by satellite electromagnetic wave signals passing through the troposphere [1]. As one of the major error sources in GNSS navigation and positioning [2], tropospheric delay increases with decreases in the satellite altitude angle, such that the slant tropospheric delay (STD) is usually projected to the zenith direction by the mapping function (MF) to simplify data processing [3,4]. Previous research has shown that the ZTD is approximately 2.3 m, and approximately 90% of it is the zenith hydrostatic delay (ZHD). The ZHD can be estimated by measuring the surface pressure and temperature. Zenith wet delay (ZWD) accounts for approximately 10% and is mainly related to the water vapor content. However, due to the fact of its uneven spatial distribution and strong random variation, it is usually difficult to accurately calculate the ZWD [5,6,7].

A common approach to solving tropospheric delay is to establish a model that estimates or predicts the tropospheric delay by fitting long-term meteorological and tropospheric delay information [8,9,10]. The accuracy of the model method is generally at the centimeter level, which could satisfy ordinary positioning requirements such as the standard point positioning (SPP). In high-precision GNSS positioning data processing, the ZTD (or ZWD) is generally regarded as an unknown parameter to estimate and it can obtain a millimetric accuracy [11]. High-precision tropospheric delay information is the basis for establishing a tropospheric delay model. However, due to the limitations in the distribution of meteorological and GNSS stations, it is difficult to acquire long-term and high-resolution tropospheric delay information on a global scale. The emergence of meteorological reanalysis technology could effectively overcome this problem. Reanalysis data realize the restoration of historical meteorological records and make up for the shortcomings of uneven spatial and temporal distributions of atmospheric observation data by integrating meteorological observation and NWP products through data assimilation technology [12]. In recent years, with the development of meteorological observation technology and the realization of more accurate mid/long-term NWP models, the quality of atmospheric reanalysis data has been further strengthened, and research on tropospheric delay inversion and modeling using reanalysis data has also increased [13,14,15].

Several national and international meteorological organizations publish their meteorological reanalysis data. The ERA5 (ECMWF Reanalysis, Version 5) is the fifth-generation global reanalysis product released by the European Centre for Medium-Range Weather Forecasts (ECMWF). It has notable improvements in data assimilation technology and resolution compared to the previous generation product, ERA-Interim [16]. Previous studies have pointed out that the global IGS stations’ tropospheric delay retrieved from ERA5 reanalysis data has higher accuracy [17]. In addition, the substantial improvement in the temporal resolution is of great significance to the study of small-scale tropospheric delay series periodic changes [18,19]. ERA5 has also been used in GNSS meteorology, playing an important role in the modeling of meteorological parameters, atmospheric weighted-average temperature (Tm), and precipitable water vapor (PWV) [20,21,22,23]. The MERRA2 (Modern-Era Retrospective Analysis for Research and Applications, Version 2) is the second-generation global reanalysis dataset released by the National Aeronautics and Space Administration (NASA). Based on the MERRA dataset, MERRA2 supplements hyperspectral, microwave, and ozone profile observations and adds a radio occultation dataset, further strengthening the system in data assimilation [24]. Some scholars have studied the applicability of MERRA2 data to tropospheric delay inversion in China, confirming that the results have good accuracy and stability [25]. The MERRA2 dataset has also been widely used in precipitation prediction [26], pollutant assessment [27], and atmospheric anomaly research [28]. In May 2021, the China Meteorological Administration (CMA) released a Chinese, first-generation global atmospheric/land surface reanalysis product, CRA40 (CMA Reanalysis, 40 years). Official data show that the quality of CRA40 products is generally comparable to that of third-generation global reanalysis products. At present, the evaluation and application of CRA40 has been gradually carried out. Some scholars have argued that the CRA40 dataset is better than other datasets for calculating temperature, relative humidity, and precipitation [29,30,31] and have proven that CRA40 performs better with respect to reproducing the climate and its changing characteristics for Mainland China [32].

The development of global atmospheric reanalysis data has given rise to new opportunities for GNSS tropospheric delay inversion and modeling. Carrying out tropospheric delay inversion and evaluation based on reanalysis data could improve the application of reanalysis data in GNSS high-accuracy positioning [33]. However, few studies have been conducted on the unified inversion and accuracy evaluation of the tropospheric delay inversion of three common types of global reanalysis data (i.e., ERA5, MERRA2, and CRA40). Moreover, the consistency between the tropospheric delay retrieved from meteorological reanalysis data and that estimated by GNSS requires further exploration. Air temperature and water vapor, which could affect the accuracy of tropospheric delay inversion, are also used as indicators for judging the type of regional climate. Therefore, whether different climate types could affect the tropospheric delay inversion is worth studying.

This paper is organized as follows: The data and methodology for the ZTD invert from meteorological reanalysis data are introduced in Section 2. The tropospheric delay consistency between the meteorological reanalysis data and GNSS, the ZTD error distribution, and its relationship with climate types are analyzed in Section 3. A summary and conclusions are provided in Section 4.

2. Data and Methodology

2.1. Data Description

The data used in this paper were ERA5 reanalysis data from ECMWF, MERRA2 reanalysis data from NASA, CRA40 reanalysis data from CMA, and the ZTD final product from IGS, ranging from 1 January to 31 December in 2020. All reanalysis data included pressure levels (PLs) and land patterns. The PL data were used to calculate the major ZTD values, while the land data were used for compensation. The accuracy of the ZTD inverted from the three types of reanalysis data was evaluated by comparing it to the IGS ZTD on a global scale.

2.1.1. Reanalysis Products

The ERA5 pressure levels product, taking pressure (P), temperature (T), specific humidity (Q), and geopotential (G) into consideration, had the highest horizontal resolution of 0.25° × 0.25° (longitude × latitude) and a vertical resolution of 37 levels from 1000 to 1 hPa. The ERA5-land hourly data, with the highest horizontal resolution of 0.1° × 0.1°, 2 m surface pressure (SP), 2 m temperature (T2m), and 2 m dewpoint temperature (dT2m), were regarded as compensation. Then, the MERRA-2 6 hourly (0, 6, 12, and 18 UTC) pressure level data had a horizontal resolution of 0.625° × 0.5° and a vertical resolution of 42 levels from 1000 to 0.1 hPa, including pressure, temperature, specific humidity, and geopotential height (GH). The MERRA2 land hourly data had the highest horizontal resolution of 0.625° × 0.5°, containing surface pressure, surface temperature (ST), and surface specific humidity (SQ). Finally, the CRA40 6 hourly pressure level data had the highest horizontal resolution of 0.3125° × 0.3125° (approximately 34 km) and a vertical resolution of 47 levels including pressure, temperature, geopotential height, and specific humidity (37 levels). The CRA40/land 3 hourly data also had the highest horizontal resolution of 34 km, involving surface pressure, 2 m temperature, and 2 m specific humidity (SQ2m). The resolution and parameter statistics of the three types of reanalysis data are shown in

Table 1. Statistical table of the resolutions and parameters of three reanalysis data types.

Data Type		Temporal Resolution	Horizontal Resolution	Vertical Resolution	Parameters	Format
ERA5	PLs	1 h	0.25° × 0.25°	37 levels	P, T, Q, G	GRIB
	land	1 h	0.1° × 0.1°	1	SP, T2m, dT2m
MERRA2	PLs	6 h	0.625° × 0.5°	42 levels	P, T, Q, GH	NC4
	land	1 h	0.625° × 0.5°	1	SP, ST, SQ
CRA40	PLs	6 h	0.3125° × 0.3125°	47 levels	P, T, Q, GH	GRIB2
	land	3 h	0.3125° × 0.3125°	1	SP, T2m, SQ2m

. In order to obtain objective evaluation results, the horizontal resolution of all reanalysis data was chosen to be 0.5° × 0.5° uniformly, except for MERRA2 (0.625° × 0.5°), and the temporal resolution was 6 h.

2.1.2. IGS ZTD

In this paper, the final ZTD products from the IGS stations were used as a reference to evaluate the performance of the ZTD inverted by reanalysis data. The IGS ZTD products have a temporal resolution of 5 min and a high accuracy of 4 mm [34]. Considering missing data (stations with more than half of the ZTD data missing were removed) and the uniform distribution, data from 436 IGS stations in 2020 were collected. The distribution of IGS stations and the corresponding mean standard deviation (STD) of the IGS ZTD are shown in Figure 1 and Figure 2. It can be seen from Figure 2 that the STDs of the selected IGS stations were almost lower than 4 mm, which meant they could be used as high-precision ZTD reference values for error evaluation and analysis.

2.2. ZTD Inversion Method with Reanalysis Data

Compared to the ZTD model, which uses meteorological parameters only in each station (e.g., Saastamoinen [8]), the integral method is more accurate for inverting tropospheric delay, and uses all meteorological parameters above the station [35]. In fact, the reanalysis data provided limited meteorological parameters, because the data above and below the pressure levels were not considered. In order to refine the current integration method, the ZTD was divided into three parts in this study as shown in Figure 3 and Equation (1).

$$\begin{array}{c}\mathrm{ZTD}=ZT{D}_{PLs}+ZT{D}_{top}+ZT{D}_{bot}\end{array}$$

(1)

where $ZT{D}_{PLs}$ is the ZTD inverted by pressure levels data, and $ZT{D}_{top}$ and $ZT{D}_{bot}$ represent the ZTD above and below the pressure level, respectively. Generally, the ZTD is derived by integrating the atmospheric refraction, which can be calculated as follows [36]:

$$\begin{array}{c}N=N_h+N_w=k_1{R}_d\rho +\left(\left(k_2-k_1\frac{R_d}{R_w}\right)\frac{e}{T}+k_3\frac{e}{T^2}\right)\end{array}$$

(2)

where $N$ is the total atmospheric refraction, and $N_h$ and $N_w$ are the hydrostatic and wet refraction values. $k_1=77.6890 \mathrm{K}/\mathrm{hPa}$, $k_2=71.2952 \mathrm{K}/\mathrm{hPa}$, and $k_3=375,463 {\mathrm{K}}^2/\mathrm{hPa}$ represent the air refractivity parameters [37]. $R_d$ and $R_w$ are constants representing dry and wet gas. The air density ($\rho$) and water vapor pressure ($e$) can be derived using the following formulas [38]:

$$\begin{array}{c}\rho =\frac{P-e}{T{R}_d}+\frac{e}{T{R}_w}\end{array}$$

(3)

$$\begin{array}{c}e\approx \frac{QP}{0.622+0.378Q}\end{array}$$

(4)

After calculating the atmospheric refraction, the $ZT{D}_{PLs}$ was derived by integrating the refractivity first:

$$\begin{array}{c}ZT{D}_{PLs}={10}^{-6}\underset{S}{\overset{ }{{\displaystyle \int }}}N·dS={10}^{-6}{\displaystyle {\displaystyle \sum }_{i=1}^n}\Delta N_{hi}·\Delta S_i+{10}^{-6}{\displaystyle {\displaystyle \sum }_{i=1}^n}\Delta N_{wi}·\Delta S_i\end{array}$$

(5)

where $S$ is the height from the bottom level to the top. $\Delta N_{hi}$ and $\Delta N_{wi}$ denote the hydrostatic and wet atmospheric refraction in the $i$th integral area (see Figure 3), and $\Delta S_i$ is the height of the $i$th integral area. According to the declining characteristic of the atmospheric refraction, exponentially or linearly with height, $\Delta N_{hi}$ and $\Delta N_{wi}$ use the following two representations [39]:

$$\begin{array}{c}\Delta N_{hi}=\left(N_{i+1}-N_i\right)/\left(\ln N_i-\ln N_{i+1}\right)\end{array}$$

(6)

$$\begin{array}{c}\Delta N_{wi}=\left(N_{i+1}-N_i\right)/2\end{array}$$

(7)

where $N_i$ and are the top and bottom levels of the $i$th integral area. Secondly, the Saastamoinen model was used to calculate $ZT{D}_{top}$ by inputting the meteorological data for the top level, and the formulas are as follows:

$$\begin{array}{c}ZT{D}_{top}=0.002277\times \frac{\left[P+\left(0.05+1255/T\right)\right]}{f\left(\phi , H\right)}\end{array}$$

(8)

$$\begin{array}{c}f\left(\phi ,H\right)=1-0.00266cos2\phi -0.00028H\end{array}$$

(9)

where $\phi$ represents the latitude, and $H$ represents the height of the top level. Finally, the land meteorological reanalysis data and the first level of the pressure-level data above the station were used to calculate $ZT{D}_{bot}$ (calculate the bottom refractive first and then multiply by the bottom distance).

As the reanalysis data are provided by grid, it was necessary to interpolate the grid ZTD to the stations. The detailed process for inversion of station ZTD using meteorological reanalysis data is shown in Figure 4. Based on Figure 4 and the above equations, the specific steps for inversion of the station ZTD using reanalysis data are summarized as follows:

(1) Original data acquisition: searching for meteorological pressure levels and land grid data for the four grid points nearest to the measurement site based on station coordinates (i.e., longitude and latitude).

(2) Unified elevation system: converting the geopotential height of the meteorological data to GNSS geodetic height plus geoid height, which can be derived from GPT3 [20].

(3) Obtaining meteorological parameters for the station height grid: the two types of meteorological data acquired in step 1 are interpolated to the station height based on 1D interpolation.

(4) Obtaining the grid ZTD of station height: the ZTD of each grid point is calculated according to Equations (1)–(9). The surface data can be brought into Equations (8) and (9) to give direct access to the grid ZTD, but this was not covered in this paper.

(5) Gaining the station ZTD: Based on the coordinates and step 4, calculating the station ZTD using 2D interpolation. The interpolation method used in this paper was inverse distance weighted (IDW) with the following equations:

$$\begin{array}{c}{D}_i=\sqrt{{\left(L-L_i\right)}^2+{\left(\phi -{\phi }_i\right)}^2}\end{array}$$

(10)

$$\begin{array}{c}{\eta }_i={\left(1/D_i\right)}^2/{\displaystyle {\displaystyle \sum }_{i=1}^4}{\left(1/D_i\right)}^2\end{array}$$

(11)

$$\begin{array}{c}f\left(P\right)={\displaystyle {\displaystyle \sum }_{i=1}^4}{\eta }_{i}f\left(Q_i\right)\end{array}$$

(12)

where $D_i$ is the distance from the unknown point $P\left(L,\phi \right)$ to the grid point $Q_i\left(L_i,{\phi }_i\right)$, ${\eta }_i$ is the weight of the $i$th grid point, and $f\left(P\right)$ is the station ZTD.

The ZTD time series derived from the above steps were compared with the IGS ZTD to test its performance. The concordance correlation coefficient (CCC) was selected to evaluate the consistency between inverted ZTD and IGS ZTD; then, the mean deviation (bias), mean absolute error (MAE), and root mean square (RMS) error were calculated to assess the ZTD accuracy for each station. The formulas for calculating CCC, bias, MAE, and RMS are expressed as follows [40]:

$$\begin{array}{c}{\rho }_c=\raisebox{1ex}{$2{\rho }_p·{\sigma }_{ra}·{\sigma }_{igs}$}\!\left/ \!\raisebox{-1ex}{$\left[{\sigma }_{ra}^2+{\sigma }_{igs}^2+{\left({\mu }_{ra}-{\mu }_{igs}\right)}^2\right]$}\right.\end{array}$$

(13)

$$\begin{array}{c}bias={\displaystyle \sum }_{i=1}^N\left(ZT{D}_{i\_ra}-ZT{D}_{i\_igs}\right)/N\end{array}$$

(14)

$$\begin{array}{c}MAE={\displaystyle \sum }_{i=1}^N\left|ZT{D}_{i\_ra}-ZT{D}_{i\_igs}\right|/N\end{array}$$

(15)

$$\begin{array}{c}RMS={\left[{\displaystyle \sum }_{i=1}^N{\left(ZT{D}_{i\_ra}-ZT{D}_{i\_igs}\right)}^2/N\right]}^{1/2}\end{array}$$

(16)

where ${\rho }_c$ is the CCC, ${\rho }_p$ is Pearson’s correlation coefficient, ${\sigma }_{ra}$ and ${\sigma }_{igs}$ represent the standard deviation of the inverted ZTD and IGS ZTD, and ${\mu }_{ra}$ and ${\mu }_{igs}$ represent the annual mean of the inverted ZTD and IGS ZTD. $N$ is the total number of samples, $ZT{D}_{i\_ra}$ is the $i$th sample value calculated by the reanalysis data, and $ZT{D}_{i\_igs}$ is the reference value of the $i$th sample from IGS.

3. Results

The three types of reanalysis data were used to invert the ZTD of the 436 IGS stations from 2020 and to assess its performance including evaluating the consistency and accuracy of the ZTD inversion from the reanalysis data and analyzing the error distribution and climate correlation.

3.1. Consistency Evaluation

The concordance correlation coefficient is a common indicator used to evaluate the consistency of two measured results. This coefficient not only reflects the correlation between the two results, but it also takes into account the actual differences, thus achieving a comprehensive representation of correlation and consistency [40]. The indicators of consistency are shown in

Table 2. The indicators of consistency (by McBride).

Consistency	|ρc|
Almost perfect	>0.90
Substantial	0.80~0.90
Moderate	0.65~0.80
Poor	<0.65

[41]. We calculated the coefficient of consistency between the ZTD inversion from each type of reanalysis data and the IGS ZTD. The global distribution of consistency is shown in Figure 5, and the statistical results are shown in

Table 3. Statistics on the number of stations for each type of consistency from the three ZTD inversions.

Reanalysis Data	Consistency (Percentage)				Average CCC
	Almost Perfect	Substantial	Moderate	Poor
ERA5	409 (93.8%)	24 (5.5%)	3 (0.7%)	0 (0)	0.960
MERRA2	358 (82.1%)	66 (15.1%)	10 (2.3%)	2 (0.5%)	0.935
CRA40	338 (77.5%)	75 (17.2%)	20 (4.6%)	3 (0.7%)	0.927

.

The average station concordance correlation coefficient between the ZTD inversion from the three types of reanalysis data and the IGS ZTD was greater than 0.9, indicating that all inverted ZTDs had high levels of consistency with the IGS ZTD. ERA5 had the highest consistency among the three reanalysis data types, with more than 90% of stations showing almost perfect concordance with the IGS ZTD, but its concordance was reduced for some stations at lower latitudes and in Antarctica. The ZTD inversion from the MERR2 and CRA40 reanalysis data showed similar consistency results to the IGS ZTD, with approximately 80% of stations globally in good agreement with the IGS ZTD, while areas of decreased consistency were also concentrated at lower latitudes, in addition to some stations in the Middle East and western North America. The current uncertainty of ZTD inversion and modeling is an irregular periodic and stochastic variation in tropospheric delay. Previous studies have shown that the stochastic variation in ZTD is related to the atmospheric turbulence in the area where the station is located [42], and the large increase or decrease in ZTD over a short period of time is often the result of the influence of strong convective weather; it is also difficult to describe such drastic changes with traditional models [43,44]. Figure 6 shows the ZTD time series for the BJFS and WUH2 stations in 2020 using the three types of reanalysis data and the ZTD change under storm conditions.

As can be seen from Figure 6, the inversion results of all three types of reanalysis data were able to better reflect the trend in periodic and small amplitude changes in ZTD as well as showed good consistency with IGS ZTD, even under rainstorms, within a short time period. Therefore, it is of great importance to carry out analysis and forecasting of tropospheric delay inversion and short-term precipitation based on reanalysis data.

3.2. Accuracy Evaluation

We evaluated the accuracy of the ZTD inversion from the ERA5 (0.5° × 0.5°), MERRA2 (0.625° × 0.5°), and CRA40 (0.5° × 0.5°) reanalysis data. The global distribution of ZTD bias, MAE, and RMS errors are shown in Figure 7, and

Table 4. Statistical results of three inverted ZTD errors (Brackets indicate the error range).

Reanalysis Data	Resolution	Bias/mm	MAE/mm	RMS/mm
ERA5	0.5° × 0.5°	−3.39 [−19.35, 11.56]	9.69 [4.49, 25.59]	12.55 [6.07, 34.57]
MERRA2	0.625° × 0.5°	0.27 [−19.19, 17.44]	12.41 [4.73, 30.03]	16.69 [5.99, 39.53]
CRA40	0.5° × 0.5°	−3.69 [−19.50, 15.99]	12.76 [4.62, 28.97]	16.96 [5.91, 38.29]

shows the ZTD error statistics for the three reanalysis data inversions.

For ERA5, the inverted ZTD accuracy was higher than for the other data, as evidenced by the global average station bias, MAE, and RMS of −3.39, 9.69, and 12.55 mm, respectively. As can be seen from Figure 7, the inverted bias of the ERA5 data was negative and mainly concentrated at low latitudes, indicating that the tropospheric delay was underestimated in these regions, while the distribution of the bias was more consistent at middle and high latitudes. The MAE of the ZTD inversion using ERA5 was mostly concentrated within 10 mm in Europe, North America, and the Southern Hemisphere, while the values for Southeast Asia, Central Africa, and Central America were between 10 and 20 mm. Finally, the majority of RMS values for the global stations were less than 20 mm for the ERA5 inversion ZTD. The spatial distribution of the RMS correlated with the latitude, showing the characteristics of a low north–south, high middle, and symmetrical distribution.

In comparison to the ERA5 data, the ZTD inversion using MERRA2 data had low accuracy, with averages of −0.27, 12.41, and 16.69 mm for the bias, MAE, and RMS, respectively. The mean bias of the MERRA2 inversion of the ZTD was close to zero, whereas most of the bias values were positive at mid and high latitudes, and the overall distribution was uneven, indicating that the MERRA2 data inversion of the ZTD had unstable systematic errors. This also showed the probability of obtaining a falsely high accuracy by using the bias to represent error. The global average MAE error was approximately 3 mm higher than ERA5, with most stations having MAE values between 10 and 25 mm, but the MAE accuracy in Antarctica and Greenland was lower than 10 mm, indicating its advantage in the inversion of polar station ZTDs. As can be seen from Table 4, the RMS of the tropospheric delay in the MERRA2 data inversion was approximately 33% lower than that of the ERA5 data, and its accuracy at low latitudes needs to be improved.

The global station ZTD accuracy, which was inverted using CRA40 data, was comparable to MERRA2, with a mean bias, MAE, and RMS of −3.69, 12.76, and 16.96 mm, respectively. In terms of the bias error distribution, the negative bias region at low latitudes was within 30° of the north and south latitudes, which had an expansion of the negative ERA5 region, and the lack of ZTD accuracy at low latitudes still occurred in CRA40. The tropospheric delay with the CRA40 data had similar MAE and RMS error distributions. Compared to the MERRA2 data inversion results, the CRA40 inversion accuracy was worse at low latitudes, with MAE and RMS errors mainly distributed between 20 and 25 mm. However, in polar and mid- to high-latitude regions, the CRA40 inversion accuracy was higher, resulting in a small difference in average accuracy between the two types of reanalysis data on a global scale. To illustrate the distribution of station ZTD errors for the three reanalysis data inversions, the numbers of stations in different error zones were counted, and the results are shown in Figure 8.

Figure 8 shows that the ZTD error using ERA5 data was less biased from the IGS ZTD than the other two types of data, and all three types of data showed an underestimation of tropospheric delay at low latitudes. The MAE errors for the reanalysis data were mainly within 20 mm with peaks at approximately 9 mm for the ERA5 and CRA data, and a slightly worse performance for the MERRA2 data. The zones with errors larger than 20 mm were mainly stations using the CRA40 data, indicating that there is the probability of large errors in the CRA40 data inversion. In general, the ZTD accuracy of the ERA5 data inversion was higher, with the majority of the RMS values concentrated between 8 and 20 mm, while the MERRA2 and CRA40 inversions were mostly concentrated between 12 and 28 mm. In addition, CRA40 outperformed MERRA2 for lower accuracies than 10 mm and had a relatively high number of stations in the lower accuracy intervals.

3.3. Climate Correlation Analysis of Errors

As can be seen from the previous section, all three types of reanalysis data showed an excellent performance at mid to high latitudes. However, this also underlines the deficiency of accuracy in inversion ZTDs at low latitudes. Research has shown that drastic changes in temperature and water vapor could affect the accuracy of ZTD inversions using tropospheric delay models at low latitudes [45]. Temperature, precipitation, and light are important indicators of the climatic conditions of a region, indicating that different climate types may affect the inversion accuracy of tropospheric delay. In order to confirm this, we analyzed the climatic distribution of the ZTD errors in the inversion from the reanalysis data in this section.

Currently, the Köppen–Geiger climate classification system, as one of the most widely used standards on a global scale, classifies the climate into five main classes (A: tropical; B: arid; C: temperate; D: cold; E: polar) and 30 subtypes based on threshold values and seasonality of monthly air temperatures and precipitation [46]. We first classified the tropospheric delay of the RMS accuracy of global stations using ERA5 data (which had the highest accuracy of the three types of data) based on the latest Köppen–Geiger climate classification. Then, variations in bias, RMS with climate, longitude, latitude, and height were analyzed. The comprehensive statistical results are shown in Figure 9.

Figure 9 shows that the tropospheric delay bias and RMS errors in the inversion of the ERA5 reanalysis data varied significantly with latitude and elevation and less obviously with the longitude. The bias and RMS values showed a symmetrical distribution along the equator in the latitude interval, with RMS values within approximately 10 mm in the middle to high latitudes (north of 30°N and south of 30°S) and between 10 and 20 mm in the low latitudes (30°S to 30°N). In addition, the bias and RMS errors decreased with increasing elevation, which may be explained by the relatively small ZTD values and variability at higher elevations as well as the stabilization of climatic conditions and the clear separation of dry and rainy periods, creating favorable conditions for reanalysis data to estimate ZTD. It should be noted that the effect of the station elevation on the inverted ZTD accuracy needs more research, since few IGS stations (approximately 7%) have been constructed at high altitudes (1500–3500 m), and higher bias and RMS accuracy can be obtained in low-altitude areas.

The tropical zone is located at low latitudes near the equator, and the temperate zone is mainly located in the southeast of the continental shelf, in the central latitude coastal or island regions, while the arid, cold, and polar climate zones show a gradual spread from low to high latitudes, thus indicating that the global climate type is related to latitude and ocean distribution. In terms of the impact of different climate types on the accuracy of the inverse ZTD of the reanalysis data, areas of large error were mainly distributed in the tropical zone, while the errors in the cold and polar climate zones were relatively insignificant. In order to obtain the accuracy variation in different climatic zones and take into account the latitude distribution of the errors, the number of stations and the distribution of errors classified according to latitude and climate were counted in this paper as shown in

Table 5. ERA5 error statistics for different classification models.

Class	Latitude					Climate					Total
Subclass	NH	NM	L	SM	SH	A	B	C	D	E
Station Number	33	215	134	41	13	91	90	165	69	21	436
Percentage/%	7.57	49.31	30.73	9.40	2.98	20.87	20.64	37.84	15.83	4.82	100
Bias/mm	–2.63	–2.06	–5.52	–2.32	–8.82	–5.92	–2.62	–1.95	–3.57	–6.50	–3.39
MAE/mm	6.43	8.61	12.68	8.30	9.32	13.88	9.18	8.59	7.89	8.20	9.69
RMS/mm	8.39	11.23	16.39	10.93	10.59	17.87	12.09	11.22	10.10	10.05	12.55

and Figure 10. (We have divided the global region into five zones by latitude: Northern Hemisphere high latitude (NH): 60°N to 90°N, Northern Hemisphere mid latitude (NM): 30°N to 60°N, low latitude (L): 30°S to 30°N, Southern Hemisphere mid latitude (SM): 30°S to 60°S, and Southern Hemisphere high latitude (SH): 60°S to 90°S.)

As the latitude decreased, all errors in the inverted ZTD of the ERA5 reanalysis data increased significantly. The mean bias, MAE, and RMS in the low-latitude region increased by approximately 38.5%, 23.6%, and 23.4%, respectively, compared with the global average, with approximately 61.9%, 32.5%, and 31.8% increasing, respectively, in the mid-latitude region, and stations in the high-latitude region increasing by approximately 20.6%, 42.9%, and 45.0%. The mean error for the middle and high latitudes was lower than the global mean, except in the southern high-latitude region, indicating that station errors at low latitudes dominated the global error distribution. For the climatic distribution of errors, the stations in the tropical zone had significantly larger errors, with the mean bias, MAE, and RMS errors increasing by approximately 42.8%, 30.2%, and 29.8%, respectively, compared to the global average, and even increasing by approximately 6.9%, 8.6%, and 8.3% compared to the stations in the lower latitudes, indicating that the main errors in the inversion of ZTD using reanalysis data originated from the tropical zone at low latitudes. Improving the accuracy of ZTD inversions in the tropical zone will be a key point of future research into ZTD inversion using reanalysis data.

By comparing the error distribution of the two different classification models in Figure 10, it can be seen that the error distribution in the climate classification model was more concentrated than in the latitude classification model, especially in those regions with larger errors, and also the error variation in each climate zone was more regular compared to the latitude region. However, the climate-based classification not only reflects the gradual decrease in errors from low latitudes to the poles but may also compensate for the disadvantage of sparse stations at high latitudes, avoiding the overall accuracy being affected by large errors caused by individual stations. Therefore, a global-scale tropospheric delay inversion and modeling study that takes climate types into account will contribute to the establishment of high-precision GNSS tropospheric delay models.

4. Conclusions

High-accuracy ZTD is essential for GNSS positioning. Based on the meteorological reanalysis data inversion tropospheric delay integration method, in this paper, we first evaluated the consistency of the inversion ZTD of three types of reanalysis data (i.e., ERA5, MERRA2, and CRA40) with reference to the IGS ZTD of the 436 stations in 2020. Secondly, the accuracy of the inversion of the global IGS stations ZTD for the three types of reanalysis data was evaluated, and the differences in the regional accuracy between the three data types were analyzed. Finally, the relationship between the ZTD errors were analyzed using the reanalysis data and climate, longitude, latitude, and height with the following main conclusions being reached:

(1) The ZTD inversion from the three types of reanalysis data had high consistency with the IGS ZTD. At the same time, the results of the ZTD inversion using the reanalysis data effectively reflected the trend and regularity of the real tropospheric delay, even under rainstorm conditions;

(2) The ZTD inversion from the ERA5 reanalysis data at the global stations had higher accuracy with bias, MAE, and RMS errors of –3.39, 9.69, and 12.55 mm, respectively. The ZTD accuracy of MERRA2 and CRA40 inversion was comparable, with bias, MAE, and RMS errors of 0.27, 12.41, and 16.69 mm for MERRA2 and –3.69, 12.76, and 16.96 mm for CRA40, respectively, which were poorer than the ERA5 reanalysis data inversions at the global scale but had higher ZTD accuracy in the polar region inversion;

(3) The accuracy of the inversion of the ZTD from the reanalysis data was closely related to climate type and latitude, and for the tropical climate zone at low latitudes, the accuracy of using the reanalysis data was lower than in other latitudes. The global errors for all three types of reanalysis data showed a higher error in the lower latitudes and a lower error in the middle to high latitudes, with a symmetrical distribution along the equator. The distribution of global errors was better represented by climate type than by latitude. Furthermore, we suggest that future research on ZTD using reanalysis data should focus on improving the accuracy in the tropical climate zone at low latitudes and tropospheric delay modeling taking climate differences into consideration.