A Combined Model to Predict GNSS Precipitable Water Vapor Based on Deep Learning

The precipitable water vapor (PWV) is a key parameter to reflect atmospheric water vapor, which can be derived by the global navigation satellite system (GNSS) technique with high accuracy and temporal resolution. PWV is an important parameter for weather forecasts and climate research. To develop a highly accurate PWV prediction model, we first combine the wavelet analysis (Wa), long short-term memory (LSTM) neural network, and autoregressive integrated moving average (ARIMA) algorithms as WLA model for the GNSS PWV prediction. Wa, LSTM, and ARIMA in WLA separate the random noise and predict the nonlinear and linear trends in PWV, respectively. Afterward, the WLA model is compared with LSTM, ARIMA, wavelet neural networks, and the multivariable linear regression (MLR) method. The WLA model shows the best result in the five prediction models in terms of the root-mean-square error (RMSE, 0.19–0.82 mm) and mean absolute error (0.01–0.07), which are 55.48% and 55.32% lower than other models, and Nash–Sutcliffe efficiency coefficient (NSE, 76.53%–99.7%) is 9.42% greater than other models. For further analysis, we also study the WLA performance in different months using one-month's data as training length. The result shows that WLA has good effects in predicting PWV in different months and the average NSE of WLA is 95%. In addition, the predicted PWV of the WLA model within 3 h is found to be accurate and reliable (RMSE < 2 mm, relative error < 0.1, NSE > 60%). This study demonstrates the good performance of WLA to predict GNSS PWV.

is very small, it is closely related to climate change [1], [2]. The amount of water vapor is difficult to measure due to strong variability in space and time. Recently, with the development of the satellite technique, GNSS becomes one of the most potent methods for water vapor detection due to its high accuracy and temporal resolution [3], [4]. Many previous studies use GNSS water vapor data as one of the important parameters to forecast short-term precipitation and the warning of extreme weather, such as torrential rains and hurricanes [5], [6], [7], [8], [9], [10]. Therefore, accurate prediction of GNSS water vapor has an important impact on climate research and weather forecast [11].
Much research in recent years has focused on predicting the precipitable water vapor (PWV) which discovered that there is potential for enhancing the accuracy of PWV prediction [12], [13], [14], [15], [16]. One way to improve the accuracy is to adopt a deep neural network that can improve the capability of the model by providing complex nonlinear systems and a higher level of abstraction for the model [17]. The deep neural network has been extensively used in meteorology, hydrology, transportation, and other fields in recent decades [18], [19], [20]. For example, Akbari Asanjan et al. [21], Khaniani et al. [22], and Liu et al. [1] completed precipitable predictions based on the long short-term memory (LSTM), support vector machine (SVM), and artificial neural network, respectively. In terms of atmospheric water vapor prediction, Pozo et al. [12] predicted PWV at ALMA site based on the YSU-Noah configuration in which root-mean-square error (RMSE) is 0.7 mm; Sharifi and Souri [13] used a hybrid LS-HE and LS-SVM model to predict PWV and the bias and standard deviations between the observed and predicted values are about 0.37 mm and 3 mm; Ge [23] and Ge et al. [24] used a wavelet neural network (WNN) to predict PWV, in which the RMSE is 0.2 mm on average; Xie et al. [25] showed that the RMSE of the genetic WNN prediction method used PWV sampled at a frequency of 5 min on August 12 at Wuhan station is 0.124 mm, which is smaller than the BP neural network and WNN; Jain et al. [14] predicted PWV using the least square estimation method, which achieves an RMSE of 0.1 mm for the forecasting internal of 5 min in the future; Huang et al. [26] conducted modeling and PWV prediction based on the improved BP neural network and the average RMSE at HKKP is about 5.142 mm. Turchi et al. [27] applied the autoregressive (AR) technique to predict PWV and RMSE < 1 mm for PWV ≤ 15 mm. An alternative approach to improve the accuracy is using the combined model to predict PWV. For instance, Liu et al. [28] used empirical mode decomposition combined with a neural network to predict PWV and Wang et al. [15] combined the improved adaptive Kalman filter and radial basis neural network to predict PWV. Compared with traditional models, the prediction accuracy of these combined models is obviously improved.
Previous experiments have established that the deep neural network and its combination models based on one or more stations have good precision in predicting PWV. However, the PWV variability is complex in different regions and seasons including linear and nonlinear terms. These models do not consider separating the random noise and the different components of linear or nonlinear trends in the PWV data, which may have an impact on the results of the forecast.
A combined model of wavelet analysis (Wa), LSTM neural network, and autoregressive integrated moving average (ARIMA) was proposed to predict GNSS PWV, hereinafter referred to as the WLA model. Wa in WLA is used to separate the random noise in PWV data. LSTM and ARIMA, respectively, predict nonlinear and linear parts of the data.

II. DATA AND METHODS
To get GNSS PWV data, observation files of ground GNSS stations need to be processed first. Tropospheric delays are estimated by GNSS data processing and are converted into PWV. After highly accurate GNSS PWV time series are derived, five models are used to predict the PWV. This part briefly introduces PWV derivation methods from ground GNSS stations and algorithms of prediction models.

A. Data Processing
We selected 22 ground GNSS stations at different latitudes in China for experiments. The latitude and longitude of 22 experimental stations are given in Table I. Among them, seven  stations (BJFS, CHAN In advance, GNSS observation files of 22 stations in 2018 are processed into hourly GNSS PWV data by GAMIT [29], [30]. GNSS signals are affected by the atmosphere, which can result in signal delay including ionospheric delay and tropospheric delay. Tropospheric delay mainly refers to the zenith tropospheric delay (ZTD). ZTD can be divided into zenith wet delay (ZWD) and zenith hydrostatic delay (ZHD), and accurate ZHD is necessary for data processing [31]. ZHD is calculated by the following equations [32], [33]: (1) where P s is the surface pressure of the station (unit: hPa); θ is the latitude of the station; H is the height of the station (unit: km); ZHD is the tropospheric zenith hydrostatic delay (unit: mm).
Finally, temperature (T ), local weighted mean temperature (T m ), the conversion factor (Π), and other parameters are used to obtain the PWV value [34], [35], [36], as shown in the following equations: where ρ w represents the density of liquid water, and the value is 1 × 10 6 kg/m 3 ; R v represents the vapor gas constant and is 461.495 J · kg −1 · k −1 ; k 2 and k 3 represent atmospheric physical constants. k 2 = 22.13 ± 2.20 k/hPa, k 3 = (3.739 ± 0.012) × 10 5 k 2 /hPa; T m represents the weighted mean temperature of the troposphere, which can be expressed as (5), where T and e are the absolute temperature and water vapor pressure, respectively, in the zenith direction.

B. Prediction Models
We use five prediction models in our studies. First, LSTM is a specific implementation of recurrent neural network (RNN) and is a powerful method for deep learning of time series data. LSTM introduces a memory cell to solve the problems of gradient disappearance and explosion that occur in RNN invariably. The memory cell contains a memory block, and each memory block has three gate structures, including the forgetting gate, the input gate, and the output gate. These three gate structures can read, write, and reset data [37]. Because the output value of the Sigmoid function is between 0 and 1 and it can Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
let the information flow through the door or not, the activation functions of the three gates are all S-type functions [38].
Second, ARIMA can identify complex patterns in data and generate predictions, which can be used to analyze and predict univariate time series data. The function of ARIMA is represented by (p, d, q), where p represents the number of autoregression items, d represents the number of nonseasonal differences, and q represents the number of lag prediction errors in the prediction equation. Three steps of establishing ARIMA are identification, estimation, and prediction [39].
Third, a WNN is formed by the combination of wavelet and neural network. Its structure is similar to the radial basis function (RBF) network and retains most of the advantages of the RBF network. Parameters of WNN such as the number of hidden nodes and weights are easier to determine. In addition, WNN has the advantages of fast convergence, high accuracy, and small network scale [40]. WNN adopts a three-layer network structure, including an input layer, hidden layer, and output layer. The number of neurons in the hidden layer is determined by the trying out method, and the activation function is the Morlet wavelet function [41].
Fourth, multivariable linear regression (MLR) refers to the method of establishing a prediction model through correlation analysis of two or more independent variables and one dependent variable. Assuming y is the dependent variable, x 1 , x 2 ,..., x n are the independent variables, and there is a linear relationship between the independent variables and the dependent variable, the MLR model is as follows: where b 0 , b 1 ,..., b n are regression coefficients of MLR model. Fifth, we combined the WLA model with Wa, LSTM, and ARIMA. Wa is an effective tool for studying nonstationary time series than Fourier Transform. Its main advantage is that it can obtain the time, position, and frequency information of the signal at the same time [39]. In the WLA model, Wa is used to separate, denoise, and reconstruct the original data to get the denoised PWV. Furthermore, the change of PWV sequence is not simply linear or nonlinear but a combination of the two. It is found in previous studies that LSTM is beneficial in predicting nonlinear sequences and ARIMA is good at predicting linear sequences, hence LSTM and ARIMA models are used to predict PWV. The nonlinear part of PWV can be largely predicted by LSTM, and the linear part of PWV can be largely predicted by ARIMA. Finally, the standard deviation weighting method is used to obtain the final PWV predicted value. WLA model not only separates the random noise existing in the original PWV sequence but also considers the linear and nonlinear trends. The main process of WLA is shown in Fig. 1.
The WLA prediction is divided into three steps.
Step 1: Wa is carried out to separate the random noise in PWV. The steps of Wa include wavelet decomposition, noise reduction, and reconstruction. "db2" is chosen as the basis function for wavelet decomposition and the number of wavelet decomposition layers is one. Noise reduction is applied to the detail coefficient after wavelet decomposition. Finally, the detail coefficient after noise reduction and the approximate coefficient  are reconstructed to obtain the denoised PWV. It showed original PWV and the sequences decomposed by Wa at CHAN in Fig. 2 as an example.
Step 2: The denoised PWV is input into the LSTM and ARIMA model to predict PWV. Enter the training and prediction set of the same time into both models, WL and WA are the predicted data of LSTM and ARIMA, respectively.
Step 3: The standard deviation weighting method, a method to determine the weight coefficient according to the standard deviation of data, is used to assign weights to LSTM and ARIMA. It is obtained by applying the following equations: where σ 2 L and σ 2 A are the variances of the predicted results of LSTM and ARIMA, w L and w A are the weights assigned to LSTM and ARIMA, respectively.
The final predicted value (Ỹ ) is obtained by the following equation:Ỹ

C. Evaluation Indicators
Evaluation indexes used in this experiment are RMSE, relative error (RE), and Nash-Sutcliffe efficiency (NSE) coefficient.
RMSE is the square root of the ratio of the sum of squares of the deviation between the predicted value and the true value (unit: mm). RE is the ratio of absolute error to the true value (no unit). The calculation formulas of RMSE and RE are shown as follows: where N is the total number, Y is the original value, and f (x i ) is the predicted value. NSE is generally used to verify the quality of hydrological model simulation results and here is used to verify the quality of PWV prediction models (−∞, 1). The closer NSE is to 1, the better the quality and reliability of the model. If NSE is much less than 0, the model is not trusted [42]. NSE (expressed as a percentage) is obtained by the following equation: where N is the total number, Y is the original value, f (x i ) is the predicted value, andȲ is the total average value of the original values.

III. EXPERIMENTS AND RESULTS
In this section, we compare the results of WLA with those of the other four models and further analyze the results of the five models in different months. In addition, experiments with different training lengths and different prediction steps based on the WLA model are carried out and the results are analyzed.

A. Model Evaluation
Five models are used to predict PWV in this part, including LSTM, ARIMA, WNN, MLR, and WLA. The prediction results of each model will be compared. These models are trained and evaluated based on one-month's data of PWV. The experiment uses 28 days' data as training data and 3 days' data (72 h) as testing data, and the previous 12 h of PWV are used to predict the next 1-h PWV. Fig. 3 shows the comparison of evaluation indicators at 22 stations of experimental results in January, and Fig. 4 shows the differences in evaluation indicators between WLA and other models. As can be seen from Fig. 3   Take BJFS station as an example, Fig. 5 shows the predicted value and error of the PWV prediction by five models at BJFS in January. It is found that the predicted value of PWV from the WLA model is closest to the original value and the error is relatively small compared with other models. Moreover, the peak of PWV value predicted by LSTM, ARIMA, WNN, and MLR in Fig. 5 has a time offset from the original one, while the value predicted by the WLA model has better agreement with the original time series by optimizing lag problems. Therefore, WLA has the best result of the PWV prediction among the five models.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  In addition, we conducted experiments on WLA and LA without adding Wa. This experiment used the previous 12 h of PWV to predict the next 1-h PWV in a year. Fig. 6 shows the variation of evaluation indicators between WLA and LA at 22 stations of one-year's data in 2018. It indicates that the RMSE and RE of WLA are smaller than LA and NSE of WLA is higher than that of LA. It can be concluded that Wa can reduce the error in predicting PWV. Fig. 7 is the predicted value and error of the PWV predicted by the WLA and LA at JFNG as an example.
PWV has a significant seasonal change, which is higher in summer and lower in winter. To further check the performance of prediction models in different months, WLA model with the same training lengths is used to predict PWV in each month in 2018. For one month, the first few days are the training period and the last three days (72 h) are the forecast period. Table II lists the range of the evaluation indicators of the different models in different months, from which we can see that the RMSE and RE of WLA are both smaller than those of other models, and NSE is larger than those of other models. Therefore, it can be concluded that the high prediction accuracy of WLA makes it more suitable for PWV prediction.    Fig. 8(a), the RMSE in January is relatively small, which is between 0.19 and 0.48 mm; the RMSE in April and October are between 0.24 and 1.06 mm and 0.20 and 0.88 mm, respectively; the RMSE in July is between 0.42 and 1.22 mm, which is relatively large. Fig. 8(b) shows that the RE in four different months is all small. The RE in January, April, July, and October are between 0.01 and 0.04, 0.01 and 0.07, 0.01 and 0.03, and 0.01 and 0.06. And RE of most stations in July is smaller than that of January, April, and October. In Fig. 8(c), we can see that the NSE of each station in four different months is around 95%. It is indicated that WLA has good effects in predicting PWV in different months. Fig. 8 shows that the RMSE of July is higher than that of other months. It also shows the situation of slightly large RE and low NSE at some stations, such as CHAN and LHAZ. For further analysis, we take CHAN and LHAZ as examples. Compared with other stations, the RMSE in July and RE in April of CHAN is significantly larger and the NSE of LHAZ in January is especially low. PWV predictions at different months    of CHAN and LHAZ are shown in Figs. 9 and 10. As can be seen from Fig. 9, the large value of PWV in July and the small value of PWV in April at CHAN make the RMSE and RE larger than in other months. What is more, PWV data at CHAN in April fluctuate less than in January, July, and October, which results in the NSE in April being the best. In Fig. 10, we can see that the PWV value of the four different months does not differ much, and therefore, the RMSE and RE do not differ much at LHAZ. The PWV of July is relatively larger and the PWV of January is smaller than other months, so the RMSE in July and RE in January are larger than those of others. On the other hand, the PWV fluctuates greatly in January, so the NSE is smaller than  in other months; the fluctuation is small in April, so the NSE is larger than in other months.
In general, the RMSE of PWV prediction in different months are all less than 1.22 mm, the RE in different months are all less than 0.1 and NSE are greater than 75.86%. It indicated that WLA can be applied to the PWV prediction of the different months with different PWV values.

B. WLA Based on Different Training Lengths
WLA model is used to predict PWV with different training lengths in this section. It predicts the last 3 days (72 h) PWV in December with the data of one month, half-year, and one year. Fig. 11 shows the evaluation indicators with different training lengths in December. For predicting the last 3 days (72 h) PWV by one-month's data, the range of RMSE, RE, and NSE is 0.28-1.01 mm, 0.02-0.1, and 83.98%-98.59%, respectively. For predicting the PWV by half-year's data, the range of RMSE, RE, and NSE is 0.26-1.02 mm, 0.01-0.1, 83.98%-98.59%, Fig. 12. Predicted PWV of different training lengths in December at AHAQ, where black represents the original PWV, blue represents the predicted PWV by one-month's data, green represents the predicted PWV by half-year's data, and red represents the predicted PWV by one-year's data. respectively. For predicting the PWV by one-year's data, the range of RMSE, RE, and NSE is 0.24-0.95 mm, 0.02-0.11, and 85.35%-98.25%, respectively.
We found that the differences in evaluation indicators among one month, half-year, and one year are very small, which indicates that increasing the training length has little effect on the accuracy of PWV prediction. It is observed that the predicted PWV by the data of one month, half-year, and one-year fit the original PWV almost as well in Fig. 12. All results show that PWV predictions of the WLA model can be made by using one month as the length of the training data.

C. WLA Based on Different Prediction Steps
PWV is of vital importance to the study of precipitation in meteorology and space geodesy. Precipitation forecast in meteorology is not simply predicting the next time but time series in the future invariably. Accordingly, this part increases the prediction step from 1 to 2, 3, 6, and 12 h.  Tables III and IV. It can be observed in Figs. 13 and 14 and Tables III and IV that the RMSE and RE of these 22 stations increase as the prediction step increases gradually. Appreciably, RMSE in June of all stations are less than 3 mm at 1 h (average 1.24 mm),    In the same way, as the prediction step increases, the NSE of the WLA is getting worse. In June, 14 of 22 stations have NSE at the prediction step in 3 h greater than 60%. NSE in December at the prediction step in 3 h is greater than 60% except for BJFS and JSLY, which illustrated that the WLA model has high reliability of PWV prediction within 3 h. When the prediction step is 6 h, the average NSE in June is 33.56%, and 12 of 22 stations have NSE greater than 40%, and 3 of 22 stations have negative NSE. And the average NSE in December is 40.66%, and 14 of 22 stations have NSE greater than 40%, and 2 of 22 stations have negative NSE. It means that the prediction results for the 6 h prediction step of the WLA model are close to the average level of original PWV. It also means that the overall results for the 6 h prediction step can be used as a reference when the accuracy requirement is not high. When the prediction step is 12 h, negative NSE is present at 16 of 22 stations in June and 17 of 22 stations in December. It has a large error and the results are not accurate, which means WLA model cannot predict PWV with 12 h of prediction steps accurately. Fig. 15 shows the results of the HKOH station using WLA to predict PWV under different prediction steps and the differences between prediction PWV and original PWV. It shows that the predicted PWV with the prediction step of 1, 2, and 3 h is verging on the original PWV while the error of 6 and 12 h is larger. In conclusion, the WLA model is reliable and accurate in predicting PWV within 3 h and is feasible to predict PWV in 6 h while the accuracy is not high.
IV. CONCLUSION WLA model that combined Wa, LSTM, and ARIMA was constructed to predict GNSS PWV. In total, 22 stations in different regions of high, low, and middle latitude in China are used to compare the WLA and other models, including LSTM, ARIMA, WNN, and MLR, and it is observed that WLA has the highest accuracy. The RMSE and RE ranges of the WLA model are 0.19-0.82 mm and 0.01-0.09, which is about 43.2% and 46.58% lower than other models. The NSE, in the range of 87.3%-99.3%, is about 17.62% greater than other models, which means that the WLA model is more stable and credible.
In addition, we analyzed the characteristics of WLA for different months, different training lengths, and different prediction steps. When comparing the results of different months, the average NSE of WLA is around 95%. And we concluded that the RMSE in January, April, and October are less than that in July, while the RE in July is smaller than that of January, April, and October. The results show that the RMSE and RE are closely related to the value of the PWV, and the WLA model has equally good prediction accuracy in different months. For experiments of different training lengths, there is little difference between the result of one month and that of half-year and one year, which can be concluded that one-month training data are enough to predict hourly PWV accurately. Furthermore, we found that the WLA model has high accuracy in predicting PWV with the prediction step of 1-3 h, and the average RMSE is less than 2 mm, RE is less than 0.1, and NSE is around 79.27%. The average RMSE and RE are 2.92 mm and 0.1 when the prediction step is 6 h. And 13 of 22 stations on average have NSE greater than 40%, indicating that the result of 6 h is credible, although the model is not stable enough.
Compared with the predicted results of other models, we found that WLA model has better accuracy in PWV prediction. The Wa in WLA can improve the prediction accuracy by decomposing the random noise from PWV time series. The LSTM and ARIMA models are combined with the help of the standard deviation weighting method. We conducted PWV prediction experiments for different months, different training lengths, and different prediction steps. The results of the PWV prediction experiment in each month show that WLA can be applied to PWV prediction in different months. What is more, we identified that the error of the PWV predicted by WLA model increased with the increase of the prediction step. And WLA can accurately predict the PWV in 3 h.
Overall, we demonstrate that the WLA model has high accuracy in the prediction of PWV and is feasible to predict PWV in the regions we selected and similar regions. It can also predict PWV in different months and different prediction steps, which will be an important promotion for the research meteorology. In future work, we can study how to improve the prediction accuracy and conduct experiments in a larger range by adding data (more stations and a longer period) to confirm the applicability of the model, which can provide powerful data support for further meteorological research.
Data Availability: The RINEX data of IGS stations are supplied by the IGS Data Center of Wuhan University (gnsswhu.cn). The RINEX data of CORS stations in Hong Kong are provided by Hong Kong Geodetic Survey Services (geodetic.gov.hk). It can provide the other data to readers by contacting the corresponding author.
Ming Shangguan received the Ph.D. degree in geodesy and geographic information science from TU Berlin, Berlin, Germany, in 2014.
She was a Researcher with GEOMAR, Kiel, Germany. From 2016 to 2020, she was a Lecturer with Southeast University, Nanjing, China. She is an Associate Professor with the School of Geography and Information, China University of Geosciences, Wuhan, China. Her research interests include GNSS meteorology, atmospheric sounding, and remote sensing.
Meng Dang received the bachelor's degree in 2020 from the China University of Geosciences, Wuhan, China, where she is currently working toward the postgraduate degree.
Her research focuses on GNSS meteorology and its applications. Yingchun Yue received the Ph.D. degree in geodesy from Wuhan University, Wuhan, China, in 2008.
She was a Post-doc with the China University of Geosciences, Wuhan, China, in 2011. She is an Associate Professor with the School of Geography and Information, China University of Geosciences. Her research interests include GNSS technology application, data processing, and GNSS meteorology in Geodesy. She is an Associate Professor with the Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan. Her research interests include geodetic crustal deformation analysis, Earth reference frame, and GNSS data processing.