A Tropospheric Zenith Delay Forecasting Model Based on a Long Short-Term Memory Neural Network and Its Impact on Precise Point Positioning

: Global navigation satellite system (GNSS) signals are affected by refraction when traveling through the troposphere, which result in tropospheric delay. Generally, the tropospheric delay is estimated as an unknown parameter in GNSS data processing. With the increasing demand for GNSS real-time applications, high-precision tropospheric delay augmentation information is vital to speed up the convergence of PPP. In this research, we estimate the zenith tropospheric delay (ZTD) from 2018 to 2019 by static precise point positioning (PPP) using the ﬁxed position mode; GNSS observations were obtained from the National Geomatics Center of China (NGCC). Firstly, ZTD outliers were detected, and data gaps were interpolated using the K-nearest neighbor algorithm (KNN). Secondly, The ZTD differences between the KNN and periodic model were employed as input datasets to train the long short-term memory (LSTM) neural network. Finally, LSTM forecasted ZTD differences and the ZTD periodic signals were combined to recover the ﬁnal forecasted ZTD results. In addition, the forecasted ZTD results were applied in static PPP as a prior constraint to reduce PPP convergence time. Numerical results show that the average root-mean-square error (RMSE) of predicting ZTD is about 1 cm. The convergence time of the PPP which was corrected by the LSTM-ZTD predictions is reduced by 13.9, 22.6, and 30.7% in the summer, autumn, and winter, respectively, over GPT2-ZTD corrected PPP and unconstrained conventional PPP for different seasons.


Introduction
Global navigation satellite system (GNSS) real-time kinematic positioning (RTK) and precise point positioning (PPP) are two popular precise positioning approaches.On the one hand, RTK eliminates common spatial-related atmospheric errors and clock errors to realize centimeter-level positioning in real-time [1].However, RTK requires a neighboring reference station, and the positioning accuracy is reduced when the distance between the rover and the reference station increases.On the other hand, PPP utilizes a standalone dual-frequency GNSS receiver with pseudo-range and carrier phase observations.By using post-mission precise satellite orbit and clock products, PPP can achieve millimeter-level absolute positioning accuracy [2][3][4][5].
In April 2013, the IGS began releasing real-time services (RTS), which consequently promote the rapid development of real-time PPP.However, the PPP applications are limited by a long convergence time, and the estimation of the tropospheric delay in solutions is one of the reasons for the long convergence time.The first-order ionospheric delay on its accuracy.For forward prediction, their model accuracy can reach 8.52 mm, given 140 epochs [33].Yang et al. (2017) constructed a ZTD forecasting model using GNSS data from Hong Kong using the UNB3m, a genetic algorithm, and a BP neural network.Although the accuracy is about 1.1 cm, the model needs more meteorological parameters, and the prediction accuracy is restricted by meteorological parameters [34].Zheng et al. (2015) constructed the ZTD forecast model in Jiangsu, China using meteorological parameters by associating the Hopfield model with the BP model.The RMSE was outstanding, at 2 mm, but the backward forecast was only 6 days [35].However, all of the above models need meteorological parameters, which brings inconvenience to the application of ZTD.Expanding the forecasting time can effectively reduce the update interval of the ZTD product and improve the prediction accuracy, which is essential for rapid PPP convergence.This research aims to expand the ZTD forecasting time using an LSTM neural network; training data involves only ZTD, without the help of additional meteorological parameters.
In this paper, the outliers in GNSS-derived ZTD were eliminated, and the ZTD gap was filled to obtain the preprocessed GNSS-ZTD data.Then, ZTD differences between the KNN ZTD and the periodic model ZTD were extracted to train the LSTM neural network and the forecast time was inputted to the trained neural network to obtain the output results.Finally, the LSTM-ZTD forecast products were recovered by combining the forecast results of the LSTM neural network with the ZTD periodic signals.The LSTM-ZTD accuracy was validated using the BPNN ZTD model, the BP-Hopfield ZTD model [35], and PPP perspectives.The convergence time of PPP between LSTM-ZTD predictions helped to correct PPP, and the GPT2 ZTD-corrected PPP was analyzed to validate the reliability of LSTM-ZTD data.Figure 1 shows the main procedures flow chart of the LSTM-ZTD forecast model.
than the Saastamoinen model, only at some sites [32].Xiao et al., (2018) constructed a tropospheric delay prediction model with a time resolution of 2 h in Japan through an improved BP neural network.They investigated the influence of the number of stations on the accuracy of the model but did not consider the effect of the duration of prediction on its accuracy.For forward prediction, their model accuracy can reach 8.52 mm, given 140 epochs [33].Yang et al. (2017) constructed a ZTD forecasting model using GNSS data from Hong Kong using the UNB3m, a genetic algorithm, and a BP neural network.Although the accuracy is about 1.1 cm, the model needs more meteorological parameters, and the prediction accuracy is restricted by meteorological parameters [34].Zheng et al. (2015) constructed the ZTD forecast model in Jiangsu, China using meteorological parameters by associating the Hopfield model with the BP model.The RMSE was outstanding, at 2 mm, but the backward forecast was only 6 days [35].However, all of the above models need meteorological parameters, which brings inconvenience to the application of ZTD.Expanding the forecasting time can effectively reduce the update interval of the ZTD product and improve the prediction accuracy, which is essential for rapid PPP convergence.This research aims to expand the ZTD forecasting time using an LSTM neural network; training data involves only ZTD, without the help of additional meteorological parameters.
In this paper, the outliers in GNSS-derived ZTD were eliminated, and the ZTD gap was filled to obtain the preprocessed GNSS-ZTD data.Then, ZTD differences between the KNN ZTD and the periodic model ZTD were extracted to train the LSTM neural network and the forecast time was inputted to the trained neural network to obtain the output results.Finally, the LSTM-ZTD forecast products were recovered by combining the forecast results of the LSTM neural network with the ZTD periodic signals.The LSTM-ZTD accuracy was validated using the BPNN ZTD model, the BP-Hopfield ZTD model [35], and PPP perspectives.The convergence time of PPP between LSTM-ZTD predictions helped to correct PPP, and the GPT2 ZTD-corrected PPP was analyzed to validate the reliability of LSTM-ZTD data.Figure 1 shows the main procedures flow chart of the LSTM-ZTD forecast model.

Materials and Methods
This study used GNSS observation data from stations of the National Geomatics Center of China (NGCC), their distribution is shown in Figure 2. The parameters and attributes of PPP are listed in Table 1 [24,36].The zenith hydrostatic delay (ZHD) is derived from the Saastamoinen model, and the zenith wet delay (ZWD) is estimated by static PPP using IGS's precise orbit and clock products, where the ZTD is obtained by adding the ZWD to the ZHD [37,38].

Materials and Methods
This study used GNSS observation data from stations of the National Geomatics Center of China (NGCC), their distribution is shown in Figure 2. The parameters and attributes of PPP are listed in Table 1 [24,36].The zenith hydrostatic delay (ZHD) is derived from the Saastamoinen model, and the zenith wet delay (ZWD) is estimated by static PPP using IGS's precise orbit and clock products, where the ZTD is obtained by adding the ZWD to the ZHD [37,38].

Data Preprocessing
Although the PPP-estimated ZTD is sampled with a 30-s interval, the ZTD interval was resampled to a one-hour resolution by retaining only the first ZTD values of each hour.To facilitate data processing, the day of the year (DOY) increases numerically from 1 January 2018, to 31 December 2019.As a result, the accumulated DOY ranges from 1.00 to 730.00, and the corresponding ZTD sample amount is 17,520 (730 × 24).

Data Preprocessing
Although the PPP-estimated ZTD is sampled with a 30-s interval, the ZTD interval was resampled to a one-hour resolution by retaining only the first ZTD values of each hour.To facilitate data processing, the day of the year (DOY) increases numerically from 1 January 2018, to 31 December 2019.As a result, the accumulated DOY ranges from 1.00 to 730.00, and the corresponding ZTD sample amount is 17,520 (730 × 24).

Outlier Elimination
The ZTD outliers are classified in the stochastic model, and all ZTDs are regarded as equal-weight observation samples [39][40][41].A fixed-length rolling window was employed to scan the long ZTD series, and if one of the ZTD samples within the window was greater than 3σ error, it was regarded as a suspicious outlier.Then, the length of the window was expanded to continue iterating until all the ZTD samples met the requirements.Figure 3 shows the ZTD scatter plot before and after the outliers were eliminated.Eliminated outliers and originally observed data gaps result in ZTD time series with discontinuities that form ZTD data gaps.

Outlier Elimination
The ZTD outliers are classified in the stochastic model, and all ZTDs are regarded as equal-weight observation samples [39][40][41].A fixed-length rolling window was employed to scan the long ZTD series, and if one of the ZTD samples within the window was greater than 3 error, it was regarded as a suspicious outlier.Then, the length of the window was expanded to continue iterating until all the ZTD samples met the requirements.Figure 3 shows the ZTD scatter plot before and after the outliers were eliminated.Eliminated outliers and originally observed data gaps result in ZTD time series with discontinuities that form ZTD data gaps.

Data Gap Interpolation
In this study, the ZTD data gaps were interpolated with the k-nearest neighbor algorithm (KNN).The KNN interpolates the missing data according to the mean value of the k-nearest samples that neighbor the gap.The principle of the KNN is to obtain ZTD samples based on the nearest Euclidean distance and interpolate the missing data by the inverse distance weighting.

Data Gap Interpolation
In this study, the ZTD data gaps were interpolated with the k-nearest neighbor algorithm (KNN).The KNN interpolates the missing data according to the mean value of the k-nearest samples that neighbor the gap.The principle of the KNN is to obtain ZTD samples based on the nearest Euclidean distance and interpolate the missing data by the inverse distance weighting.
Input parameters for the KNN algorithm are the sampling time and the corresponding ZTD.However, the ZTD series are discontinuous and unevenly spaced.A complete time series with hour-resolution were fed into the KNN algorithm, with the output parameter being interpolated ZTDs.The KNN was implemented with the 'knn-regression' module in the scikit-learn package [42].The k value was set to 720 for model training.The KNN only fills ZTD gaps, hence it retains the original ZTD series.
To validate the proposed data gap interpolation method, we use the periodic model as a reference, which is constructed by Equation (1).

ZTD(T) =
where T denotes the accumulated DOY from 2018 to 2019, A 0 is the mean value of the ZTD, and (A 1 , B 1 ) and (A 2 , B 2 ) are the coefficients of the periodic model, which were calculated with the least squares method.
Figure 4 shows the error comparison of 16-h ZTD interpolated data between the KNN and the periodic model at the CQPS station.
Input parameters for the KNN algorithm are the sampling time and the corresponding ZTD.However, the ZTD series are discontinuous and unevenly spaced.A complete time series with hour-resolution were fed into the KNN algorithm, with the output parameter being interpolated ZTDs.The KNN was implemented with the 'knnregression' module in the scikit-learn package [42].The k value was set to 720 for model training.The KNN only fills ZTD gaps, hence it retains the original ZTD series.
To validate the proposed data gap interpolation method, we use the periodic model as a reference, which is constructed by Equation (1).
where T denotes the accumulated DOY from 2018 to 2019,  0 is the mean value of the ZTD, and ( 1 ,  1 ) and ( 2 ,  2 ) are the coefficients of the periodic model, which were calculated with the least squares method.
Figure 4 shows the error comparison of 16-h ZTD interpolated data between the KNN and the periodic model at the CQPS station.Since the data gap length at each station is not uniform, we take the station CQPS as an example to demonstrate how the KNN fills ZTD gaps.The interpolated ZTDs via the periodic model and the k-nearest neighbor algorithm are smaller than 2 cm from 10:00 to 16:00.However, from 17:00 to 01:00, the errors increase to 6 cm, which implies the accuracies of both models decrease as the gap length increases.Furthermore, after 22:00, the k-nearest neighbor algorithm error begins to decrease from 5.3 to 3.8 cm, whereas the periodic model error remains at about 5.6 cm.We investigated the model accuracy indicators (bias, STD, RMSE, MAE, and R), where the equations for the indicators are as follows: (2)-( 5 (2) Since the data gap length at each station is not uniform, we take the station CQPS as an example to demonstrate how the KNN fills ZTD gaps.The interpolated ZTDs via the periodic model and the k-nearest neighbor algorithm are smaller than 2 cm from 10:00 to 16:00.However, from 17:00 to 01:00, the errors increase to 6 cm, which implies the accuracies of both models decrease as the gap length increases.Furthermore, after 22:00, the k-nearest neighbor algorithm error begins to decrease from 5.3 to 3.8 cm, whereas the periodic model error remains at about 5.6 cm.We investigated the model accuracy indicators (bias, STD, RMSE, MAE, and R), where the equations for the indicators are as follows: (2)-( 5), where N refers to the number of samples, ZTD  2. It demonstrates that the KNN is superior to the periodic model.Therefore, we utilize the KNN to interpolate the ZTD data gaps in this study.The global geodetic observing system (GGOS) provided VMF ZTDs with a time resolution of 6 h and a spatial resolution of 1 • × 1 • .We downloaded VMF3 ZTD products from 1 January 2018 to 31 December 2019 at https://vmf.geo.tuwien.ac.at/trop_products/GRID/1x1/ (accessed on 5 April 2020).The distribution of GNSS stations is sparse.We determine the VMF grid nearest to the GNSS station as a reference, where the criteria are that the longitude and the latitude between approximate GNSS positions and the VMF grid are both less than 0.5 • and the altitude difference is less than 100 m.The details are shown in Table 3.We resampled the above preprocessed ZTDs (after eliminating the outliers and filling the ZTD gaps) sampling interval to 6 h by retaining the first ZTD values of every 6 h to validate GNSS ZTD testing accuracy, which were evaluated using Bias, STD, MAE, RMSE, and R. The statistical accuracy indicators are shown in Figure 5.The mean values for Bias, STD, MAE, RMSE and R, which are the testing accuracy indicators, for GNSS ZTDs at GNSS stations are −0.04,2.44, 1.78, 2.51 and 0.87, respectively.

ZTD Differences Extraction
After data preprocessing, we obtained the ZTD sequences without outliers and data gaps.ZTD sequences were divided into two parts: training sets (T: 1-510) and testing sets (T: 510-730).The training set was modeled by Equation (1). Figure 6 shows the KNN, the

LSTM-ZTD Model Construction 2.2.1. ZTD Differences Extraction
After data preprocessing, we obtained the ZTD sequences without outliers and data gaps.ZTD sequences were divided into two parts: training sets (T: 1-510) and testing sets (T: 510-730).The training set was modeled by Equation (1). Figure 6 shows the KNN, the periodic model derived ZTDs, and their differences for the entire dataset.It demonstrates that the periodic model can reflect the annual and semi-annual variation trends.However, their differences vary within 20 cm, without obvious patterns.All seven GNSS stations show negligible bias.The maximum, minimum, and mean RMS values were 6.0, 4.1, and 4.8 cm, respectively.

LSTM-ZTD Model Construction
LSTM is one implementation of a recurrent neural network, which can effectively represent time series with long-term patterns [43], and it has been utilized in many fields, including meteorology, network traffic, and air pollution forecasting [44][45][46][47][48].In the previous section, we extracted differences between the k-nearest neighbor and the periodic model.In this section, we use these differences to train an LSTM neural network.Then, the trained network was used to forecast the ZTD difference during the coming epoch, and finally, the ZTD, by adding the periodic model ZTD.The proposed forecasting model is referred to as LSTM-ZTD in this study.
Figure 7 shows that the LSTM neural network is constructed with multiple layers.The data for the input layer contain ZTD differences (∆  ) between the KNN and the periodic model, with corresponding T. The hidden layer is composed of three LSTM layers (LSTM1, LSTM2, and LSTM3) and their corresponding dropout operation (Dropout1, Dropout2, and Dropout3).These representative features can be extracted as hidden layers from historical input data and connected to a dense layer, which is the fully connected layer that synthesizes features to obtain composite information,   .The output layer processes   to output the forecasting results, ∆  .

LSTM-ZTD Model Construction
LSTM is one implementation of a recurrent neural network, which can effectively represent time series with long-term patterns [43], and it has been utilized in many fields, including meteorology, network traffic, and air pollution forecasting [44][45][46][47][48].In the previous section, we extracted differences between the k-nearest neighbor and the periodic model.In this section, we use these differences to train an LSTM neural network.Then, the trained network was used to forecast the ZTD difference during the coming epoch, and finally, the ZTD, by adding the periodic model ZTD.The proposed forecasting model is referred to as LSTM-ZTD in this study.
Figure 7 shows that the LSTM neural network is constructed with multiple layers.The data for the input layer contain ZTD differences (∆ZTD t ) between the KNN and the periodic model, with corresponding T. The hidden layer is composed of three LSTM layers (LSTM1, LSTM2, and LSTM3) and their corresponding dropout operation (Dropout1, Dropout2, and Dropout3).These representative features can be extracted as hidden layers from historical input data and connected to a dense layer, which is the fully connected layer that synthesizes features to obtain composite information, y t .The output layer processes y t to output the forecasting results, ∆ZTD f .Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 16 A total of 70% (510 days) of the ZTD dataset was used for model training and 30% (220 days) was used for model validation.We chose to use adaptive moment estimation (ADAM) as the optimizer and RMS as the loss function indicator.
Figure 8 demonstrates the network training and validation flow chart.We use the ZTD differences between the KNN and the periodic model over 29 days in a sliding window to predict the next day.One day was treated as a group and the input parameters were divided into 510 groups.The activation function for fully connected layers is tanh, the dropout rate for each LSTM layer is 0.01, the batch size is 40, and the total epochs for training is 10.To obtain the optimum network hyperparameters, the cross-validation strategy was utilized in the training set.More specifically, the training set was divided into 10 equal parts, where 9 parts were used for training, and the remaining part was used for validation [47,49,50].It takes 6.64 h to train the model on a supercomputing system in the Supercomputing Center of Wuhan University.The supercomputing system is equipped with 2580.175MHz, 20 cores CPU, and 128 G RAM.A total of 70% (510 days) of the ZTD dataset was used for model training and 30% (220 days) was used for model validation.We chose to use adaptive moment estimation (ADAM) as the optimizer and RMS as the loss function indicator.
Figure 8 demonstrates the network training and validation flow chart.We use the ZTD differences between the KNN and the periodic model over 29 days in a sliding window to predict the next day.One day was treated as a group and the input parameters were divided into 510 groups.The activation function for fully connected layers is tanh, the dropout rate for each LSTM layer is 0.01, the batch size is 40, and the total epochs for training is 10.To obtain the optimum network hyperparameters, the cross-validation strategy was utilized in the training set.More specifically, the training set was divided into 10 equal parts, where 9 parts were used for training, and the remaining part was used for validation [47,49,50].It takes 6.64 h to train the model on a supercomputing system in the Supercomputing Center of Wuhan University.The supercomputing system is equipped with 2580.175MHz, 20 cores CPU, and 128 G RAM.A total of 70% (510 days) of the ZTD dataset was used for model training and 30% (220 days) was used for model validation.We chose to use adaptive moment estimation (ADAM) as the optimizer and RMS as the loss function indicator.
Figure 8 demonstrates the network training and validation flow chart.We use the ZTD differences between the KNN and the periodic model over 29 days in a sliding window to predict the next day.One day was treated as a group and the input parameters were divided into 510 groups.The activation function for fully connected layers is tanh, the dropout rate for each LSTM layer is 0.01, the batch size is 40, and the total epochs for training is 10.To obtain the optimum network hyperparameters, the cross-validation strategy was utilized in the training set.More specifically, the training set was divided into 10 equal parts, where 9 parts were used for training, and the remaining part was used for validation [47,49,50].It takes 6.64 h to train the model on a supercomputing system in the Supercomputing Center of Wuhan University.The supercomputing system is equipped with 2580.175MHz, 20 cores CPU, and 128 G RAM.

Results and Analysis
Firstly, ZTDs from T 510 to 730 are used to validate the LSTM-ZTD forecasting model.Secondly, the forecasted ZTDs are used for tropospheric delay correction for PPP.

Validation from ZTD Models Perspective
Figure 9 shows the ZTD series comparison for 220-day forecasts using different models (i.e., the BPNN, the BP-Hopfield, the KNN, and the LSTM-ZTD model).Among these models, the KNN derived ZTD series was selected as the reference.The LSTM-ZTD results are close to those of the KNN derived ZTD series, indicating that the proposed model can forecast complex ZTD changes.The BPNN model and the BP-Hopfield model products can only just reflect the approximate trend.

Results and Analysis
Firstly, ZTDs from T 510 to 730 are used to validate the LSTM-ZTD forecasting model.Secondly, the forecasted ZTDs are used for tropospheric delay correction for PPP.

Validation from ZTD Models Perspective
Figure 9 shows the ZTD series comparison for 220-day forecasts using different models (i.e., the BPNN, the BP-Hopfield, the KNN, and the LSTM-ZTD model).Among these models, the KNN derived ZTD series was as the reference.The LSTM-ZTD results are close to those of the KNN derived ZTD series, indicating that the proposed model can forecast complex ZTD changes.The BPNN model and the BP-Hopfield model products can only just reflect the approximate trend.

Validation from PPP Perspective
The tropospheric wet delay is generally estimated as an unknown parameter in PPP.The tropospheric delay can be expressed by the following model: where   0 modeled tropospheric delay,  ℎ () and   () are the projection functions of ZHD and ZWD, respectively, ∆  , is the unmodelled ZWD,   ,ℎ is the modeled ZHD, and   ,0 is the modelled ZWD [22,49,51]. ℎ () and   () are the NMF (Niell 1996), which do not require any meteorological data, and E is the elevation angle [24,50].More details of the corresponding processing solutions are shown in Tables 1 and 4. Tropospheric augmentation information can reduce PPP convergence time [22,52,53].Mendez Astudillo et al., (2018) found that the convergence time of PPP is affected by seasonal variations [49], and the forecasted ZTD series from T 510 to 730 involve three seasons, namely summer, autumn, and winter.The positioning experiment was carried out using observational data and ZTD series were forecasted for three GNSS stations (JSHZ, SXQY, SXYQ) in the summer (T 514, May 29), autumn (T 608, August 31), and winter (T 693, November 24).We investigated the convergence time to verify the feasibility of LSTM-ZTD corrected PPP (L-PPP), as well as GPT2 ZTD corrected PPP (G-PPP).Unlike conventional ionospheric-free combined PPP, LSTM-ZTD corrected PPP and GPT2 ZTD corrected PPP add a virtual observation for ZTD and its corresponding constraint on the observation equation [21].In this study, we define the convergence time as a reverse search from the last epoch to the first epoch, until the error during one epoch in one of the E, N, or U directions exceeds 10 cm.
Table 5 shows the convergence times of LSTM-ZTD corrected PPP compared with GPT2 ZTD corrected PPP and conventional PPP.It shows that the LSTM-ZTD corrected PPP convergence time is shorter than GPT2 ZTD corrected PPP and conventional PPP.The convergence time of GPT2 ZTD corrected PPP was reduced by 4.8%, and the

Validation from PPP Perspective
The tropospheric wet delay is generally estimated as an unknown parameter in PPP.The tropospheric delay can be expressed by the following model: where T r 0 is the modeled tropospheric delay, M h (E) and M w (E) are the projection functions of ZHD and ZWD, respectively, ∆T r z,w is the unmodelled ZWD, T r z,h is the modeled ZHD, and T r z,w0 is the modelled ZWD [22,49,51].M h (E) and M w (E) are the NMF (Niell 1996), which do not require any meteorological data, and E is the elevation angle [24,50].More details of the corresponding processing solutions are shown in Tables 1 and 4.  [49], and the forecasted ZTD series from T 510 to 730 involve three seasons, namely summer, autumn, and winter.The positioning experiment was carried out using observational data and ZTD series were forecasted for three GNSS stations (JSHZ, SXQY, SXYQ) in the summer (T 514, May 29), autumn (T 608, August 31), and winter (T 693, November 24).
We investigated the convergence time to verify the feasibility of LSTM-ZTD corrected PPP (L-PPP), as well as GPT2 ZTD corrected PPP (G-PPP).Unlike conventional ionosphericfree combined PPP, LSTM-ZTD corrected PPP and GPT2 ZTD corrected PPP add a virtual observation for ZTD and its corresponding constraint on the observation equation [21].In this study, we define the convergence time as a reverse search from the last epoch to the first epoch, until the error during one epoch in one of the E, N, or U directions exceeds 10 cm.
Table 5 shows the convergence times of LSTM-ZTD corrected PPP compared with GPT2 ZTD corrected PPP and conventional PPP.It shows that the LSTM-ZTD corrected PPP convergence time is shorter than GPT2 ZTD corrected PPP and conventional PPP.The convergence time of GPT2 ZTD corrected PPP was reduced by 4.8%, and the convergence time of LSTM-ZTD corrected PPP was reduced by 13.9% in the summer.Those two solutions were reduced by 13.7% and 22.6% in the autumn season and reduced by 13.8% and 30.7% in the winter season, respectively, compared to conventional PPP.The C-PPP denotes conventional PPP, G-PPP denotes GPT2 ZTD corrected PPP and L-PPP denotes LSTM-ZTD corrected PPP.

Conclusions
In this study, we proposed an LSTM-ZTD forecasting model.GNSS observations provided by the NGCC during 2018-2019 were used to establish and validate the proposed model.Firstly, the ZTDs were estimated using static PPP and the outliers were eliminated.Secondly, the k-nearest neighbor algorithm was applied to interpolate ZTD data gaps.Numerical results show that the KNN is superior to the periodical model in interpolating ZTD data gaps.Thirdly, ZTD differences between the k-nearest neighbor and periodic model were used to train an LSTM neural network.LSTM forecasted ZTD differences and periodic model derived ZTDs were combined to recover the forecasted ZTDs.Finally, a static PPP experiment was performed to validate and develop the LSTM-ZTD model.
Conclusions are summarized as follows: (1) The ZTD forecast error of the proposed LSTM model is about 1.0 cm, superior to that of the Hopfield (8.4 cm), BPNN (6.3 cm), and BP-Hopfield (11.1 cm), and periodic model (5.0 cm); (2) The GPT2 ZTD corrected PPP convergence time was reduced by 4.8, 13.7, and 13.8% in the summer, autumn, and winter seasons, respectively.The LSTM-ZTD corrected PPP convergence time was reduced by 13.9, 22.6, and 30.7% in the summer, autumn, and winter seasons, respectively, which are superior to conventional PPP.

Figure 1 .
Figure 1.Main procedures flow chart to establish the LSTM-ZTD forecast model.

Figure 1 .
Figure 1.Main procedures flow chart to establish the LSTM-ZTD forecast model.

Figure 2 .
Figure 2. Distribution of GNSS stations in this study.

Figure 2 .
Figure 2. Distribution of GNSS stations in this study.

Figure 3 .
Figure 3. ZTD time series of seven GNSS stations before (left) and after (right) the outliers were eliminated.A red box marks obvious outliers in the left panel.N = 730, representing two years, collected from 1 January 2018 to 31 December 2019.

Figure 3
Figure 3 illustrates that ZTDs are mostly continuous, and only parts of the ZTD data are missing.The statistical results show that the maximum data gap length is no more than 7 days, where most of the ZTD data gap lengths are several hours.

Figure 3 .
Figure 3. ZTD time series of seven GNSS stations before (left) and after (right) the outliers were eliminated.A red box marks obvious outliers in the left panel.N = 730, representing two years, collected from 1 January 2018 to 31 December 2019.

Figure 3
Figure 3 illustrates that ZTDs are mostly continuous, and only parts of the ZTD data are missing.The statistical results show that the maximum data gap length is no more than 7 days, where most of the ZTD data gap lengths are several hours.

Figure 4 .
Figure 4.The KNN and periodic model interpolated ZTD error comparison from 10:00 on 3 December 2018 to 1:00 on December 4 at the CQPS station.

Figure 4 .
Figure 4.The KNN and periodic model interpolated ZTD error comparison from 10:00 on 3 December 2018 to 1:00 on December 4 at the CQPS station.
pre i is the ZTD value that is output from the models, ZTD obs i is the original ZTD value as a reference, and ZTD obs i and ZTD pre i are the mean values of ZTD obs i and ZTD pre i , respectively.The periodic model and the KNN-filled ZTD gap within 16 h are shown in Table

Figure 6 .
Figure 6.Comparison between the KNN and the periodic model.N = 730, the blue and orange curves show the KNN and periodic model derived ZTDs, respectively.The green curve shows ZTD series differences between the two models.

Figure 6 .
Figure 6.Comparison between the KNN and the periodic model.N = 730, the blue and orange curves show the KNN and periodic model derived ZTDs, respectively.The green curve shows ZTD series differences between the two models.

Figure 7 .
Figure 7. LSTM neural network structure designed in this study.

Figure 8 .
Figure 8. Network training and validation flow chart.The red box shows the input parameter, which consists of T and the ZTD differences (∆  ) between the KNN and the periodic model.The blue box expresses the training sets, the orange box shows the validation sets, and the red circle represents the original differences between KNN and the periodic model.The green box shows the output forecasting results, ∆  .

Figure 7 .
Figure 7. LSTM neural network structure designed in this study.

Figure 7 .
Figure 7. LSTM neural network structure designed in this study.

Figure 8 .
Figure 8. Network training and validation flow chart.The red box shows the input parameter, which consists of T and the ZTD differences (∆  ) between the KNN and the periodic model.The blue box expresses the training sets, the orange box shows the validation sets, and the red circle represents the original differences between KNN and the periodic model.The green box shows the output forecasting results, ∆  .

Figure 8 .
Figure 8. Network training and validation flow chart.The red box shows the input parameter, which consists of T and the ZTD differences (∆ZTD t ) between the KNN and the periodic model.The blue box expresses the training sets, the orange box shows the validation sets, and the red circle represents the original differences between KNN and the periodic model.The green box shows the output forecasting results, ∆ZTD f .

Figure 9 .
Figure 9.The ZTD series and their errors.Comparisons among different models.Each panel represents a station, where the blue, purple, cyan, and red curves show the BPNN, BP-Hopfield, KNN, and LSTM-ZTD models, respectively.The dark green, orange, and light green curves show the errors in the BP-Hopfield, BPNN, and LSTM-ZTD models, respectively.

Figure 10
Figure 10 shows the RMSE comparison among the BPNN, BP-Hopfield, period model, and LSTM-ZTD errors.It demonstrates that the LSTM-ZTD model is superior to the other models, by more than about 1 cm.The mean RMSE error of the LSTM-ZTD model for seven stations is 1.0 cm.

Figure 9 .
Figure 9.The ZTD series and their errors.Comparisons among different models.Each panel represents a station, where the blue, purple, cyan, and red curves show the BPNN, BP-Hopfield, KNN, and LSTM-ZTD models, respectively.The dark green, orange, and light green curves show the errors in the BP-Hopfield, BPNN, and LSTM-ZTD models, respectively.

Figure 10
Figure 10 shows the RMSE comparison among the BPNN, BP-Hopfield, period model, and LSTM-ZTD errors.It demonstrates that the LSTM-ZTD model is superior to the other models, by more than about 1 cm.The mean RMSE error of the LSTM-ZTD model for seven stations is 1.0 cm.

Table 1 .
Strategies used in the PPP program.

Table 1 .
Strategies used in the PPP program.

Table 2 .
It demonstrates that the KNN is superior to the periodic model.Therefore, we utilize the KNN to interpolate the ZTD data gaps in this study.
), where N refers to the number of samples,    is the ZTD value that is output from the models,    is the original ZTD value as a reference, and    ̅̅̅̅̅̅̅̅̅̅ and    ̅̅̅̅̅̅̅̅̅̅ are the mean values of    and    , respectively.The periodic model and the KNN-filled ZTD gap within 16 h are shown in

Table 2 .
Accuracy comparison between the KNN and the periodical model for data gap interpolation.

Table 3 .
Spatial information of approximate GNSS station positions and the VMF grid.

Table 4 .
Processing solutions of PPP.

Table 4 .
Processing solutions of PPP.