SHAtropE—A Regional Gridded ZTD Model for China and the Surrounding Areas

Abstract

A regional zenith tropospheric delay (ZTD) empirical model, referred to as SHAtropE (SHanghai Astronomical observatory tropospheric delay model—Extended), is developed and provides tropospheric propagation delay corrections for users in China and the surrounding areas with improved accuracy. The SHAtropE model was developed based on the ZTD time series of the continuous GNSS sites from the Crustal Movement Observation Network of China (CMONOC) and GNSS sites of surrounding areas. It combines the exponential and periodical functions and is provided as regional grids with a resolution of 2.5° × 2.0° in longitude and latitude. At each grid point, the exponential function converts the ZTD from the site height to the ellipsoid, and the periodical terms, including both annual and semi-annual periods, describe ZTD’s temporal variation. Moreover, SHAtropE also provides the predicted ZTD uncertainty, which is valuable in Precise Point Positioning (PPP) with ZTD being constrained for faster convergence. The data of 310 GNSS sites over 7 years were used to validate the new model. Results show that the SHAtropE ZTD has an accuracy of 3.5 cm in root mean square (RMS) quantity, which has a mean improvement of 35.2% and 5.4% over the UNB3m (5.4 cm) and GPT3 (3.7 cm) models, respectively. The predicted uncertainty of SHAtropE ZTD shows seasonal variations, where the values are larger in summer than in winter. By applying the SHAtropE model in the static PPP, the convergence time of GPS-only and BDS-only solutions are reduced by 8.1% and 14.5% respectively compared to the UNB3m model, and the reductions are 6.9% and 11.2% respectively for the GPT3 model. As no meteorological data are required for the implementation of the model, the SHAtropE could thus be a refined tropospheric model for GNSS users in mainland China and the surrounding areas. The method of modeling the ZTD uncertainty can also be used in further global tropospheric delay modeling.

1. Introduction

Through the path that radio signals traverse the troposphere, the signals are delayed to a magnitude of meters with respect to free-space propagation [1]. This delay, referred to as tropospheric delay, is a function of the tropospheric refractive index, and is dependent on the local temperature, pressure, and relative humidity. Tropospheric delay is one of the main error sources in space geodesy techniques, i.e., Global Navigation Satellite System (GNSS), Very Long Baseline Interferometry [2], satellite altimetry [3,4], etc. Tropospheric delay is usually modeled as the sum of a hydrostatic and a non-hydrostatic (wet) part, where the hydrostatic delay (referred as ZHD at zenith) accounts for 90% of the zenith total delay (ZTD), while the wet delays (referred as ZWD at zenith) contributes to the precipitable water vapor are hard to model since the water vapor varies a lot [5]. Consequently, wet delays are usually estimated as unknown parameters in GNSS data analysis. On the other hand, using the precise a priori tropospheric delay in GNSS data processing with proper constraint can improve the convergence time and the short-term precision [6,7,8,9,10,11].

The discrete tropospheric delays can be derived from radiosonde observations and Numerical Weather Model (NWM), e.g., European Centre for Medium-Range Weather Forecasts (ECMWF). The latter one has a global coverage and a ZTD accuracy of 1–2 cm can be achieved [12,13]. On the other hand, the empirical tropospheric delay models are widely used, as they provide real-time data and do not need any external information. The empirical models could be determined using discrete ZTD data from radiosonde, NWM, and GNSS products. The radiosonde observations are used for the early-stage models, e.g., the Hopfield model [14], the Saastamoinen [15], the UNB series models (UNB1 through UNB4) [16,17] and EGNOS model [18], etc. The Global Pressure and Temperature (GPT) [19,20] model is determined using the ECMWF, and provides the pressure and temperature at any location of the earth surface using the spherical harmonic function. The improved version of GPT, i.e., GPT2 [21], GPT2w [22], and GPT3 [23], all use a global grid with better accuracy and more tropospheric information. All these models determine the global meteorological parameters, e.g., pressure, temperature, water vapor, and then they are used as inputs of other models, e.g., the Saasamoinen and the Askne and Nordius [24] models to calculate tropospheric delays. Recently, some new models have provided tropospheric delay directly, such as the IGGtrop [25,26] and GZTD [27,28], where the ZTD time series from NWM are directly modeled. It has been demonstrated that by appropriately constraining the tropospheric parameter using the empirical models, the PPP convergence time could be accelerated [7,9].

With the expanding of GNSS ground tracking network, precise ZTD products have been derived with improved temporal and spatial resolutions. Covering the whole area of mainland China, the GNSS network of the Crustal Movement Observation Network of China (CMONOC) consists of around 260 well-distributed continuous sites. The GNSS analysis center at Shanghai Astronomical Observatory (abbreviated as SHA, [29]) has been providing the precise ZTD time series since 2011. By using the precise ZTD products of SHA in mainland China and the precise ZTD products of surrounding GNSS sites from the Nevada Geodetic Laboratory [30], we develop a regional tropospheric delay empirical model, namely SHAtropE (SHanghai Astronomical observatory tropospheric delay model—Extended), for China and surrounding areas. More importantly, for the first time we provide the a priori predicted uncertainty of the empirical ZTD at any time and any location, which can be used as a priori constraint in GNSS PPP to accelerate its convergence. The accuracy of other empirical models is available in several studies about the model assessment, but none of these models provide the epoch-wise location-specified ZTD uncertainty, which is achieved in this study.

This paper presents the details of the empirical SHAtropE ZTD model and its application in the PPP. Section 2 introduces the input ZTD data, where 7 years of data from January 2012 to December 2018 of 310 sites are used and their spatial and temporal characteristics are analyzed. Section 3 describes the method for the determination of SHAtropE model. Section 4 validates the SHAtropE model, where ZTD values of the 6 years are used for inner-consistency validation and the ZTDs in the year of 2018 are used for the validation of the predicted period [31], and one-month data are used for PPP assessment. Section 5 summaries the main conclusions.

2. Input Data

The CMONOC GNSS network starts full operation since 2011, and the SHA GNSS Analysis Center includes this regional network into its routine global GNSS network data analysis [29,32], providing site-wise ZTD products together with other precise products, e.g., site coordinates, satellite orbits and clocks, etc. The left panel in Figure 1 shows the distribution of these sites, including 223 sites from CMONOC and 87 sites from the local GNSS networks and IGS (http://www.igs.org/network). The sites have a good coverage over this region except for the Tibetan Plateau. The precise ZTD estimates from January 2012 to December 2018 of all CMONOC GNSS sites and other surrounding GNSS sites are used, where data of the first 6 years are used for the SHAtropE model determination and the data of the whole 7 years are used for the model assessment. The right panel in Figure 1 shows the ZTD usable days at each site during January 2012–December 2017, of which most sites have a period of valid data of more than 1500 days, and a minimum time of more than 2 years.

The SHA GNSS analysis center follows the IGS data analysis strategies using a zero-difference approach based on 5-min sampling data of selected global GNSS network [29,32]. SHA provides a precise ZTD product with a resolution of 1 h. The ZTD product has an agreement of around 4 mm compared to the IGS ZTD products [33], thus it is reliable to be used in the empirical ZTD modeling. The ZTD products of the Non-CMONOC sites are from the Nevada Geodetic Laboratory (NGL, http://geodesy.unr.edu/). For details regarding these sites we refer to Blewitt et al. [30]. All precise ZTD products from SHA and NGL are referred to as RAW ZTD in the following text.

Temporal Variation of RAW ZTD Time Series

The ZTD of a specific site at a specific day could be expressed as the combination of an empirical model and the residuals at this day [34]:

$$\begin{array}{l}ZTD\left(doy\right)=ZT{D}_{\mathrm{mod}}\left(doy\right)+res\left(doy\right)\\ ZT{D}_{\mathrm{mod}}\left(doy\right)=A_0+A_1\cos \left(\frac{2\mathsf{\pi }}{365.25}\right(doy-d_1\left)\right)+A_2\cos \left(\frac{4\mathsf{\pi }}{365.25}\left(doy-d_2\right)\right)\end{array}$$

(1)

where $ZTD\left(doy\right)$ is the RAW ZTD, $ZT{D}_{\mathrm{mod}}\left(doy\right)$ is the modeled ZTD and $res\left(doy\right)$ is the fitting residual; $A_0$ is the average annual ZTD; A₁ and d₁ are the amplitude and initial phase of the annual component; A₂ and d₂ are the amplitude and the initial phase of the semi-annual component.

The residual cannot be precisely modeled due to the high variable water vapor component and the atmosphere pressure fluctuation. It is for this reason that the empirical tropospheric delay suffers from accuracy loss. The residual variation, however, is highly related to the weather condition, e.g., the fluctuation in summer is much larger than that in winter, thus, the empirical model suffers a larger error in summer than in winter. The absolute value of the residuals can be modeled as:

$$\begin{array}{l}\left|res\left(doy\right)\right|=re{s}_{\mathrm{mod}}\left(doy\right)+dres\left(doy\right)\\ re{s}_{\mathrm{mod}}\left(doy\right)=B_0+B_1\cos \left(\frac{2\mathsf{\pi }}{365.25}\left(doy-f_1\right)\right)+B_2\cos \left(\frac{4\mathsf{\pi }}{365.25}\left(doy-f_2\right)\right)\end{array}$$

(2)

where $\left|res\left(doy\right)\right|$ is the absolute value of the RAW ZTD residuals, $re{s}_{\mathrm{mod}}\left(doy\right)$ is the modeled absolute value of the ZTD residuals and $dres\left(doy\right)$ is the fitting residual; B₀ is the average annual absolute value of empirical ZTD residual; B₁ and f₁ are the amplitude and initial phase of annual component; B₂ and f₂ are amplitude and initial phase of semi-annual component. By definition, the term $re{s}_{\mathrm{mod}}\left(doy\right)$ presents the uncertainty of the above-defined empirical tropospheric delay, and it can be used to predict the uncertainty of the empirical model, thus provides better a priori constraint of troposphere delays when applying the model in positioning.

Using 6 years of data, we first derive the empirical temporal ZTD variations model for each site. Figure 2 shows the model parameters and the fitting RMS of each site according to their geographical distribution. The fitting RMS values (panel (a)) show large differences among regions, where the largest RMS value is 7 cm and the smallest one is 1 cm. The fitting errors are generally smaller in western parts of the studied area and larger in east parts, where the content of water vapor and its variation is much larger. For the fitted parameters, each parameter shows spatial similarity, where the parameters are quite close for nearby sites. The spatial variation of the constant term A₀ (panel (b)) for all sites is from 1.4 to 2.5 m, where nearby sites show large similarity, and a rough linear regression between the value and the site height could be observed. The annual amplitude A₁ (panel (c)) ranges from −2 to −14 cm, where the maximum (absolute value) is in the vicinity of (25°N, 80°E) and the minimum is near (40°N, 65°E). The semi-annual amplitude A₂ (panel (e)) ranges from 1 to −5 cm, where the maximum is in the vicinity of (25°N, 80°E).

The fitted parameters of the absolute ZTD residuals $\left|res\left(doy\right)\right|$ using Equation (2) are presented in Figure 3. The annual average value of the absolute biases B₀ (panel (b)) ranges from 1 to 6 cm, and the maximum value is in the vicinity of (30°N, 120°E). This is consistent with Figure 2a, where the largest value of fit RMS is in the same area. The annual amplitude B₁ (panel (c)) could be as large as 2 cm, while the semi-annual amplitude B₂ (panel (e)) could be as large as 1.5 cm. The initial phases of annual and semi-annual terms (panel (d) and (f)) both show a similar pattern for nearby sites.

It should be noted that the RAW and modeled ZTD residuals, i.e., $res\left(doy\right)$ and $re{s}_{\mathrm{mod}}\left(doy\right)$ are different in terms of RMS statistics. The ratio between the RMS of RAW ZTD fit residuals $RMS\left(res\right)$ and the RMS of modeled residuals $RMS\left(re{s}_{\mathrm{mod}}\right)$ could be defined as: $ratio=\frac{RMS\left(res\right)}{RMS\left(re{s}_{\mathrm{mod}}\right)}$, which is shown in Figure 3a. The modeling of residuals introduces a down-scale of the RMS statistics in a magnitude of about 1.2, showing the modeled residuals $re{s}_{\mathrm{mod}}$ account for about 80% of the real residuals. Thus, to keep the consistency of modeled residuals and the RAW residuals, this ratio value is used for further model determination.

3. Model Determination of SHAtropE

It is widely acknowledged that ZTD is highly correlated with site height; thus, it should be carefully considered in the empirical ZTD model determination. There are two main approaches to account for the influence of the height difference on ZTD modeling. The first is to set up fitting models with site heights as input [26], and the other one is to correct the ZTD difference according to the site height difference with respect to a reference surface [35]. In the modeling of the SHAtropE, we choose the ellipsoid as the reference plane and the empirical temporal ZTD variations model of all sites are covert to this plane. The convert model functional parameters of each site are then interpolated to derive the functional parameters at pre-defined grids.

3.1. ZTD with the Ellipsoid as Reference Surface

Considering the site height shown in Figure 1a and the constant term A₀ in Figure 2b, the annual average ZTD is reversely proportional to site height. Figure 4 shows the relation of the annual average ZTD and ellipsoidal height of each site, where sites are separated into five groups according to site latitudes and each group is represented by different color, and all the sites are defined as the sixth group. For each group, the relation between annual average ZTD at the site ellipsoidal height $h$ ($ZT{D}_h$), and the annual average ZTD on the ellipsoid surface ($ZT{D}_e$), can be expressed as an exponential formula [36]:

$$ZT{D}_h=ZT{D}_e\times e^{-\frac{h}{H_0}}$$

(3)

where $e$ is Euler’s number, $H_0$ is the tropospheric scale height.

The fitted exponential curves for each group are shown in Figure 4, where the tropospheric scale height for each group is also presented (as in

Table 1. The fitted tropospheric scale height for different latitude groups.

Latitude	(m)
15°N–25°N	7065
25°N–30°N	7459
30°N–35°N	7519
35°N–40°N	7696
40°N–55°N	7747

). As shown in Figure 4, the RMS values of fit residuals are (0.6, 1.4, 1.6, 1.3, 1.1, 5.5) cm for the six groups, respectively. The fit RMS value (5.5 cm) of all the sites together is much larger than that of all other five separated groups. Therefore, by dividing the sites into different groups according to the latitude, the fitting result shows a significant improvement.

3.2. ZTD Temporal Variations on the Ellipsoid

Equation (3) coverts the ZTD from site ellipsoidal height $h$ to the ellipsoid, and the ZTD variation due to the site height differences is thus removed. Combining Equations (1) and (3), the empirical ZTD model for each site can thus be written as:

$$ZT{D}_h=\left[A_{0,e}+A_{1,e}\cos \left(\frac{2\pi }{365.25}\left(doy+d_1\right)\right)+A_{2,e}\cos \left(\frac{4\pi }{365.25}\left(doy+d_2\right)\right)\right]\times e^{-\frac{h}{H_0}}$$

(4)

where A_0,e, A_1,e and A_2,e are the mean ZTD, ZTD annual amplitude and ZTD semi-annual amplitude of the site on ellipsoid. With the tropospheric scale height $H_0$ from Table 1 (also Figure 4), A_0,e, A_1,e, and A_2,e can be directly calculated using the fitted parameters A₀, A₁, and A₂ in Section 2, through the following equations:

$$A_{0,e}=A_0÷e^{-\frac{h}{H_0}};\begin{array}{ccc}& A_{1,e}=A_1÷e^{-\frac{h}{H_0}};& \end{array}{A}_{2,e}=A_2÷e^{-\frac{h}{H_0}}$$

(5)

The parameters of A_0,e, A_1,e, and A_2,e of all the sites are plotted in Figure 5. The spatial variation of the term A_0,e (panel (a)) for all sites is from 2.4 to 2.6 m, where its variation range is much smaller than that of A₀ in Figure 2. The annual amplitude A_1,e (panel (b)) ranges from 2 to 14 cm, and the semi-annual amplitude A_2,e ranges from 1 to 5 cm (panel (c)).

3.3. Gridded ZTD Modeling of SHAtropE

As shown in the previous section, the parameters of empirical temporal ZTD variations model are rather close for nearby sites. Thus, we set up a spatial-correlated grid model, i.e., SHAtropE, for the study areas. The modeling steps are as follows:

• For each site, RAW ZTD from SHA (CMONOC sites) and NGL (Non-COMONC sites) are used to derive the five functional parameters A₀, A₁, d₁, A₂, d₂ and the five ZTD uncertainty functional parameters B₀, B₁, f₁, B₂, f₂, based on Equations (1) and (2), respectively.
• For each site, the five ZTD functional parameters of empirical model are converted to the ellipsoid using the exponential function and the constants in in Table 1, and the related parameters on the ellipsoid, i.e., A_0,e, A_1,e, and A_2,e are derived.
• Divide the study areas [70°E–135°E, 18°N–54°N] as grids with a resolution of 2.5° and 2.0° on longitude and latitude, respectively. The 5 ZTD functional parameters on the ellipsoid of each grid point is derived by the inverse distance weighted (IDW) function [37] using the parameters of nearby sites. And the ZTD uncertainty functional parameters B₀, B₁, f₁, B₂, f₂ for each grid could also be derived by the IDW approach.

The exponential constants of latitude groups, the ten ZTD/uncertainty functional parameters of each grid point are the basic components of the SHAtropE. The characteristics of the data sets used in this study are summarized in

Table 2. Characteristics of the data used in this study.

Raw Observation	RAW ZTD Spans	Spatial Resolution	Spatial Distribution
CMONOC, NGL	January 2012–December 2018	2.5° (Longitude) × 2.0° (Latitude)	70°E–135°E, 18°N–54°N

.

4. Assessment of SHAtropE

SHAtropE provides users with ten ZTD/uncertainty functional parameters of each grid point in an ASCII-file (30 KB with the resolution of 2.5° (Longitude) × 2.0° (Latitude). On the user side, the calculation of the ZTD and its uncertainty using the SHAtropE model requires the input of site location (longitude, latitude, ellipsoidal height) and the day of the year. In the first step, the user calculates the five ZTD functional parameters (A_0,e, A_1,e, d₁, A_2,e, and d₂) and uncertainty functional parameters (B₀, B₁, f₁, B₂, f₂) at its location using the nearby grid points with the bilinear interpolation scheme. Then the ZTD at the user location is calculated according to Equation (4) in the second step. As for the uncertainty of SHAtropE ZTD, user first calculate the modeled value according to Equation (2), and this value is multiplied by 1.2 to derive the uncertainty value.

In the following, we use three methods to assess the performance of the SHAtropE model. The ZTD from SHAtropE model is called SHAtropE ZTD in the following text. Three methods are: (1) Comparisons between SHAtropE ZTD and RAW ZTD using data from January 2012 to December 2017 of all the used sites. This comparison shows inner-accuracy, as the ZTDs of each site during this period are used for the SHAtropE modeling. (2) Comparisons between SHAtropE ZTD and RAW ZTD using the whole year (2018) data of all sites. This comparison shows external-accuracy, because the predicted SHAtropE ZTD values are used and validated. (3) Precise Point Positioning using SHAtropE ZTD as the a priori input, where the ZTD is constrained using the predicted uncertainty of the model.

4.1. ZTD Accuracy of SHAtropE

For the ZTD comparisons, RAW ZTD is chosen as reference and the mean bias and RMS are calculated as follows:

$$\begin{array}{l}RMS=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(ZT{D}_i^{RAW}-ZT{D}_i^{Model}\right)}^2}}\\ BIAS=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left(ZT{D}_i^{RAW}-ZT{D}_i^{Model}\right)}\end{array}$$

(6)

where $ZT{D}^{RAW}$ and $ZT{D}^{Model}$ are RAW ZTD and SHAtropE ZTD.

For comparison, the other two widely used models, e.g., UNB3m and GPT3, are also evaluated using the RAW ZTD products. UNB3m is a representative model based on radiosonde observations, and it uses a lookup table, where the different parameters are divided by latitude. GPT3 is the latest version of the tropospheric delay models developed by TU Wien. It is based on the numerical weather model, and the product is provided with a global grid. GPT3 model has two available resolutions, i.e., 1° × 1° and 5° × 5°. In this study the 1° × 1° grid product is used, as it has a better accuracy than the other one.

Figure 6 shows the RAW ZTD, SHAtropE ZTD and its uncertainty, and the differences between RAW/SHAtropE ZTD from January 2012 to December 2018 for sites NMER (at north-east China with low site altitude) and YNLJ (at south-west China with high site altitude). The RMS values of the ZTD differences between RAW ZTD and SHAtropE ZTD are 2.0 cm and 1.9 cm for the two sites, and that of the predicted uncertainty of SHAtropE (black line) are 1.8 cm and 2.1 cm, indicating a good agreement of the predicted uncertainty. Relative larger variations are observed in summer seasons for both sites. For site NMER with lower ellipsoidal height, the ZTD variation in summer season is so big that the residuals may even beyond ±2 times of the predicted uncertainty; while for site YNLJ with higher ellipsoidal height, the ZTD variation is almost all within ±1 times of the predicted uncertainty. As for the modeled uncertainty, apparent seasonal variation is observed. Comparing the results of the modeling and prediction periods, i.e., the 2012–2017 and 2018, respectively, no apparent divergence is observed, showing that SHAtropE is very reliable in both predicted ZTD and its uncertainty values.

Figure 7 shows the RMS and BIAS of the whole period from January 2012 to December 2018 for all sites. For comparisons, the same ZTD quantities are derived for UNB3m/GPT3 models. As shown in Figure 7 that SHAtropE has no significant bias for all sites, with the maximum bias of 2.9 mm and minimum one of −3.2 mm. It achieves the best accuracy with the RMS less than 3.5 cm for all sites, and the maximum RMS appears in the vicinity of (30°N, 120°E). UNB3m shows the worst accuracy with mean RMS of 5.4 cm, and the BIAS and RMS are apparently separated into two groups, where sites in lower latitude have much worse performance than those in higher latitude. GPT3 has a similar performance compared to SHAtropE, with a slightly larger RMS value of 3.7 cm.

Figure 8 shows the mean BIAS and RMS of all sites, arranged into the 2012–2017 and 2018 groups, where the first group represents the inner-accuracy and the latter for external-accuracy. For the mean BIAS, SHAtropE is of (0.0, 0.2) cm for inner-accuracy and external-accuracy, while they are of (0.4, 0.5) cm and (−0.2, 0.1) cm for UNB3m and GPT3, respectively. For the RMS, they are of (3.5, 3.6) cm, (5.4, 5.6) cm, and (3.7, 3.8) cm for SHAtropE, UNB3m and GPT3, respectively.

4.2. The Predicted ZTD Uncertainty of SHAtropE

Unlike other empirical ZTD models, for the first time SHAtropE provides not only the ZTD model but also the ZTD uncertainty. In this section we evaluate the predicted ZTD uncertainty of SHAtropE by comparing it to the RMS statistics of the differences between SHAtropE ZTD and RAW ZTD.

As shown in Figure 6, the predicted SHAtropE ZTD uncertainty is time variable. Figure 9 shows the predicted uncertainty of SHAtropE ZTD (in left panels) and the differences between the uncertainty and the RMS value of ZTD biases (the value is Section 4.1). The left subplots of Figure 9 show the seasonal mean uncertainty quantities, where four seasons are defined according to activities of the ZTD variation. The predicted ZTD uncertainty of SHAtropE shows clear seasonal variation, where the largest values are demonstrated in Jun-Jul-Aug, especially for sites with longitude greater than 105°E, while it is less than 2 cm in the season of Dec-Jan-Feb for most sites. The right panels in Figure 9 show the differences between the predicted ZTD uncertainty and the actual RMS of SHAtropE (shown in Figure 7). In general, the predicted ZTD uncertainty of SHAtropE has a good agreement with the RMS value of SHAtropE ZTD, with mean difference of less than 5 mm, and they show no seasonal variation for all sites.

Figure 10 shows the statistics of the RMS and predicted uncertainty of SHAtropE of all sites during different seasons, separating into the 2012–2017 and 2018 groups. Again, the statistics show the overall divergence of these two quantities is less than 2 mm for both modeling and prediction.

4.3. Precise Point Positioning Performance Improvement Using SHAtropE Model

To evaluate the effect of three types of tropospheric model (SHAtropE, UNB3m and GPT3) on the GNSS PPP performances, GPS and BDS observation datasets collected from 29 CMONOC sites for 1–28 February 2019 (DoY from 32 to 59) are selected and utilized for statistical analysis. All CMONOC sites are equipped with Trimble NetR9 receivers. Details of the models and strategies related to data processing for GPS/BDS static PPP are shown in

Table 3. Adopted models and strategies for multi-Global Navigation Satellite System (GNSS) static Precise Point Positioning (PPP).

Item	Models/Strategies
Frequency selection	GPS: L1/L2; BDS: B1/B2
Estimator	Kalman filter
Sampling rate	30 s
Elevation cutoff angle	10°
Satellite orbit and clock	Fixed to GFZ final orbit and clock offset products
Satellite differential code bias (DCB)	Correct using MGEX DCB products
Receiver and Satellite antenna	GPS PCO (phase center offset)/PCV (phase center variation) corrected with igs14.atx, BDS PCO corrected with the value released by ESA and PCV is not considered
Tropospheric delay	Modeled for the dry part and estimated for wet part as random-walk noise process; GMF [38] mapping function applied
Ionospheric delay	Eliminated by Ionosphere-free combinations
Tidal effects	Corrected by IERS Convention 2010, including solid tide and ocean tide loading [39]
Relativistic effects	Corrected by model
Phase windup	Corrected by model [40]
Weighing strategy	A priori precision of 0.003m and 0.3m for GPS phase and code; A priori precision of 0.003m and 0.6m for BDS phase and code; Elevation-dependent weighing (1 for $E>30°$ otherwise $2\cdot \sin (E)$) is used
Site coordinates	Estimated as constants
Receiver clock	Estimated as white noise process
Phase ambiguities	Estimated as float constants for each arc

. For the assessment of the PPP solution, the precise coordinate solution with an uncertainty of few millimeters from SHA analysis center is used as the reference.

Constraining the ZTD in PPP can speed up the ambiguity convergence process, especially for the vertical component. In the PPP tests, we define the convergence epoch for GPS PPP as the epoch, where the coordinate error is less than 10 cm in all the east (E), north (N) and up (U) components, while the convergence criterion of BDS PPP is set to 20 cm. This setting is according to the following two facts. One is that the accuracy of BDS precise orbit and clock offset is lower than that of GPS, the other point is that BDS-2 constellation consists mostly of GEO and IGSO satellites, resulting in a worse geometry of BDS observations compared to GPS [41]. Figure 11 shows the improvement of average convergence time of static PPP with SHAtropE model for GPS-only and BDS-only solutions compared to UNB3m and GPT3 models. The average convergence time of static GPS-only/BDS-only PPP with three tropospheric models for all test sites over one month is presented in Figure 12. For the GPS-only PPP with SHAtropE model, the convergence time of half of the test sites is shortened by more than 2 min and over 6 min for 5–6 sites. The average convergence time is reduced by 8.1% compared to the UNB3m model, and the reduction is 6.9% compared to the GPT3 model. In terms of BDS-only PPP, the improvement of convergence performance by applying the SHAtropE model is more obvious than the GPS-only solution, and the convergence time of most sites is reduced by 6–10 min. Compared with the UNB3m and GPT3 models, the time reduction of the SHAtropE model are up to 14.5% and 11.2%, respectively.

5. Conclusions

This study presents a regional ZTD empirical model, SHAtropE, for the regions of [70°E–135°E, 18°N–54°N]. The SHAtropE is a gridded model with a resolution of 2.5° × 2.0°, where at each grid SHAtropE provides users with five periodical ZTD functional parameters (A_0,e, A_1,e, d₁, A_2,e, d₂) together with five uncertainty functional parameters parameters (B₀, B₁, f₁, B₂, f₂). Users could derive their local functional parameters by bilinear interpolation and calculate the ZTD over ellipsoid. By applying the derived tropospheric scale height H₀ in this study, the ZTD over ellipsoid can be coverted to the true ZTD over site height.

Using 7 years’ data of 310 GNSS sites, the new model is derived and comprehensively validated. Results show that the ZTD precision of SHAtropE is of 3.5 cm in RMS, which has a mean improvement of 35.2% and 5.4% over the UNB3m and GPT3 models, respectively. SHAtropE also provides the predicted ZTD uncertainty, which is valuable when constraining the ZTD in PPP for faster convergence. By applying the SHAtropE model in the static PPP, the convergence time of GPS-only and BDS-only solutions is reduced by 8.1%/14.5% compared to the UNB3m model, and the reduction is 6.9%/11.2% for the GPT3 model.

The input of the SHAtropE is site location (longitude, latitude, ellipsoidal height) and day of the year, and neither temperature and pressure nor humidity measurements are required for the implementation of the model. Thus, its implementation is rather convenient. It is provided in an ASCII-file (30 KB with resolution of 2.5° (Longitude) × 2.0° (Latitude)) containing 10 parameters for each grid and is available at the SHA webpage (www.shao.ac.cn/sha_gnss_ac).