Modeling of Residual GNSS Station Motions through Meteorological Data in a Machine Learning Approach

Abstract

Long-term Global Navigation Satellite System (GNSS) height residual time series contain signals that are related to environmental influences. A big part of the residuals can be explained by environmental surface loadings, expressed through physical models. This work aims to find a model that connects raw meteorological parameters with the GNSS residuals. The approach is to train a Temporal Convolutional Network (TCN) on 206 GNSS stations in central Europe, after which the resulting model is applied to 68 test stations in the same area. When comparing the Root Mean Square (RMS) error reduction of the time series reduced by physical models, and, by the TCN model, the latter reduction rate is, on average, 0.8% lower. In a second experiment, the TCN is utilized to further reduce the RMS of the time series, of which the loading models were already subtracted. This yields additional 2.7% of RMS reduction on average, resulting in a mean RMS reduction of 28.6% overall. The results suggests that a TCN, using meteorological features as input data, is able to reconstruct the reductions almost on the same level as physical models. Trained on the residuals, reduced by environmental loadings, the TCN is still able to slightly increase the overall reduction of variations in the GNSS station position time series.

1. Introduction

The first static Global Navigation Satellite System (GNSS) stations have been established almost three decades ago. Today, we profit from very long observation time series, where the height component allows us to resolve important information about the vertical movement of the Earth’s crust. This movement is affected by long-term trends, linear drifts, seasonal motions, and offsets. The seasonal variations are dominated by annual and semi-annual periodicities, that can partly be explained by the so-called environmental surface loadings. These can be categorized into hydrological (HYDL), non-tidal atmospheric (NTAL), and non-tidal oceanic loading (NTOL).

Already, in 1994, van Dam et al. [1] started analyzing the impact of atmospheric pressure on the GNSS height variance and (also van Dam et al. [2]) correlate atmospheric and oceanic leading with geoid deformation and further with GNSS height deviations. Based on global GNSS solutions, Dam et al. [3] also observe a correlation of continental water storage and vertical crust movement. Gegout et al. [4] introduce environmental loadings in the processing of GNSS height solutions and are able to improve the height residuals, but mainly in the northern hemisphere. In the southern hemisphere, most of their test stations were located closer to coastal regions, and the assumption was that the missing improvement can be traced back to mis-modeled tropospheric zenith delay. Recent studies look into the comparison of different environmental surface loading models [5,6,7,8], the comparison of geophysical models and information drawn from the Gravity Recovery and Climate Experiment (GRACE) [9], and the assimilation of both approaches [10,11]. Up to today, there is no standardized procedure of modeling and reducing the station position by environmental surface loadings [12]. Most studies achieve a Root Mean Square (RMS) improvement of around 20–40%, and some suggest mis-modeled tropospheric delays as a source of the remaining height variability [4,12].

Machine learning has become more important in recent years and is also emerging in the field of GNSS time series analysis, including a wide range of applications [13]. Several studies combine parameters derived from GNSS observations, such as precipitable or integrated water vapor (PVW, IVW) and meteorological data, in a deep learning approach to forecast heavy rainfall [14,15,16] or to perform storm nowcasting [17]. Another cluster of applications of using of machine learning and GNSS can be found in the modeling and prediction of tropospheric zenith delay [18,19,20,21]. In a recent study, Özdemir [22] take the approach of machine learning, but, instead of modeling the tropospheric delay first, they directly correct the GNSS position by meteorological parameters with the help of an artificial neural network (ANN) in a real-time application.

In this study, we are working on the post-processing of the dense network of GNSS stations in central Europe. We aim to improve GNSS residual height time series, by taking into account environmental loadings and then further reduce the remaining signal. In a first step, we take the GNSS height residuals and train the TCN directly on the meteorological parameters, and then compare it with the subtraction of the physical model. In a second step, we first reduce the GNSS residuals by the environmental loading models and then train the TCN on the already reduced time series. The idea is to find a model to partly explain the remaining signal and use it to further reduce the height residuals.

2. Data and Pre-Processing

2.1. GNSS Height Residuals

The Nevada Geodetic Laboratory (NGL) collects data from all available geodetic GNSS stations worldwide and makes their station position solutions publicly available [23]. For this work, the 24-h final solutions in the area of central Europe are used, where a network of stations with a nearby meteorological station is available (see Section 2.3). The dataset has daily sampling, ranging from January 1994 to October 2020. Only the stations that have a time span of minimum 3.5 years and not more than 20% of missing observations are used. Figure 1 shows the distribution and data availability of the selected NGL stations. The raw data files were downloaded from NGL (http://geodesy.unr.edu/index.php, accessed on 13 November 2020) in tenv3 format and further processed in the software package Hector [24]. Outliers are removed, and discontinuities are detected and corrected for. Then, a linear trend is fitted and subtracted, resulting in the detrended and cleaned signal.

2.2. Environmental Surface Loadings

The vertical displacements, caused by environmental surface loadings (HYDL, NTAL, NTOL), are available from the Earth System Modeling Group Repository of Deutsche Geoforschungszentrum Potsdam (ESMGFZ) [25]. All loadings are stored in a regular 0.5${}^{\circ }$× 0.5${}^{\circ }$ grid, and, while HYDL is provided in 24 h sampling, NTAL and NTOL have a 3 h sampling rate. First, NTAL and NTOL are downsampled to 24-h observation spacing by taking the daily average. In a next step, to reduce the GNSS residuals by the loading data, their values on the exact station positions are needed. These are derived from the grid by conducting a bilinear interpolation. Additionally, all resulting HYDL time series are reduced by a second-order polynomial fit, in order to correct for a second-order long-term trend, appearing at some stations. In Figure 2, the final individual loadings, as well as the sum of them, are depicted for the example station MAR3. More details about the preprocessing of the environmental surface loadings, as well as of the GNSS height residuals described in the previous section, can be found in Reference [26].

2.3. Meteorological Data

Meteorological data is provided by the European Climate Assessment and Dataset (ECAD) [27]. The dataset is downloadable as a predefined subset of daily observations, sorted by stations and parameters. The used dataset is blended, meaning that nearby stations (12.5 km distance and 25 m height difference) are used to fill any datagaps (https://www.ecad.eu/helptext/blending_help.html, accessed on 11 December 2021). This work includes time series of the daily mean temperature (degree celcius), precipitation amount (millimeters), sea level pressure (hectopascal), humidity (percent), and radiation (Watt per square meter). In Figure 3, all parameters are shown for the example station MAR3. The only stations used are those where all listed parameters are available, and not more than 20% of observations are missing. In order to match GNSS and meteo stations, the closest meteorological station to GNSS station was selected, with a maximum distance of 50 km. If no suitable meteorological station could be allocated, the GNSS station had to be removed from the dataset, which led to the final distribution of used GNSS stations. Figure 4 shows all selected meteorological and GNSS stations, as well as the distances from each GNSS station to the corresponding meteorological station.

2.4. Reduction of GNSS Residuals by Environmental Loadings

In Figure 5, an overview of the processing steps for the reduction of GNSS residuals by environmental surface loadings is shown. The original time series (GNSS and environmental loadings) are processed individually, according to the steps highlighted in the orange boxes, to $GNS{S}_p$ and ${\left[Loading\right]}_p$. The sum of the processed loading time series $SU{M}_p$ is then subtracted from $GNS{S}_p$, resulting the the reduced GNSS residuals $GNS{S}_{RED}$. From $GNS{S}_p$ and $GNS{S}_{RED}$, the Root Mean Square reduction $RM{S}_{RED}$ is computed. In Figure 6, the different processing steps are shown on one example station (MAR3) from the raw data to the loading reduced residuals. Figure 7 shows the spatial distribution and histogram of the achieved RMS reductions for all GNSS stations that were later used as test stations for the machine learning algorithm.

The average RMS reduction for these stations is 25.9%, with 98.5% of stations being positively reduced. The comparison to the state of the art in literature is not directly possible, as every study uses a different set of stations and loading models. Similar studies, which also focus on the area of Europe, are Bian [5], achieving an average RMS reduction of 16–25%, and Springer et al. [28], obtaining 20–30% with NTAL and NTOL and adding another 7% with HYDL.

3. Methodology

3.1. Temporal Convolutional Network

A Temporal Convolutional Network (TCN) is based on a sequence modeling problem, where a mapping function $f:{\mathcal{X}}^{T+1}\stackrel{}{\to }{\mathcal{Y}}^{T+1}$ maps an input sequence of length T, $x_0,\dots ,x_T$ to an output $y_0,\dots ,y_T$. The sequence modeling network tries to find the function f that minimizes the expected loss L between the true and predicted values: $L(y_0,\dots ,y_T,f(x_0,\dots ,x_T))$ [29]. A popular method for the analysis of a long time series is the Long Short-Term Memory (LSTM) Network, which is based on a Recurrent Neural Network (RNN). In recent research, it was shown that TCN often outperforms LSTM as it is using both the advantages of an RNN, by processing high-level features, as well as the low-level feature computation from Convolutional Neural Networks (CNN) [29,30,31,32]. By using a radically different network architecture compared to RNNs, TCN is also overcoming any gradient vanishing problems, that can occur in deep networks with very long input sequences [33]. The two main features of a TCN are:

(1)

the network can take an input sequence of any length and return the same length as output by using a 1D fully convolutional architecture, and

(2)

all convolutions are causal; therefore, $y_T$ can only depend on $x_0,\dots ,x_T$, and not any future inputs $x_{t+1},\dots ,x_T$ [29,31]. In Figure 8, the principle of dilated causal convolutions is depicted.

3.2. Implementation Details

After the previously described pre-processing steps (selection of stations and reduction by environmental surface loadings), all remaining datagaps are filled by linear interpolation over the whole time series. Overall, 342 stations are selected, from which 60% (206 stations) are used for training, 20% (68 stations) for validation, and 20% (68 stations) for testing. This ratio was adopted due to its wide use in literature and additional tests, that indicated better performance when increasing the number training stations. All stations were allocated randomly to the subsets. Figure 9 shows the spatial distribution of the station categories. The input features per station are all meteorological parameters (temperature, precipitation, sea level pressure, humidity, and radiation), as well as location information (latitude, longitude, and height) in the form of constants over the whole time series. Both the input feature matrix and the target vector are scaled between 0 and 1 for the training of the network, and the results are scaled back to the original range again for analysis. Each training sample is represented by a sequence of F features, which can be represented as a matrix of size (FxT), with F being the number of Features, and T the time. During training, samples are grouped in batches of size B, leading to a 3-dimensional object of size (BxFxT), where B is the number of elements in the batch. In

Table 1. Parameter values of the TCN architecture.

Parameter	Value(s)
batch size B	2
number of features F	8
sequence length (epochs) T	7
filter size k	8
dilations d	1, 2, 4, 8, 16, 32, 64, 128, 256
number of filters	24

below, the utilized parameter values can be found, which were determined empirically. The TCN is implemented in Python, based on Bai et al. [29], and a GitHub repository (https://github.com/philipperemy/keras-tcn, accessed on 11 December 2020).

Figure 10 gives a schematic overview of the training, validation, and testing of the TCN. After dividing the dataset into its subsets, the TCN is trained on the meteorological parameters as input features, and the GNSS residuals as target features. In the validation, the meteorological parameters from the validation dataset are put through the TCN and then compared to the validation GNSS residuals to get a measure of how well the TCN is performing, as well as given information on the ability to generalize. When the model is sufficiently trained, the meteorological time series of the test dataset are used to model the GNSS time series GNSS_MOD. These are subtracted from the corresponding GNSS residuals from the test dataset, which results in the reduced GNSS time series GNSS_RED.

4. Results and Discussion

4.1. TCN Modeling and GNSS Reduction and Comparison to Physical Loading Model Reduction

This section shows the results of modeling the GNSS residuals through the TCN using time series of meteorological parameters as input parameters. The modeled signal is then subtracted from the GNSS time series and the RMS reduction is computed. These results are compared to the reduction of GNSS residuals by physical loading models from GFZ. The results of a TCN-modeled signal in the time and frequency domain is shown on the example station MAR3 in Figure 11. There is a clear seasonal pattern in the modeled signal, with the strongest peak in the annual amplitude. However, the semi-annual and other inter-annual amplitudes also follow the pattern of original GNSS residuals. Figure 12 shows the RMS reduction rates for all stations when subtracting the TCN-modeled signal from the GNSS residuals. All stations are positively reduced, and the mean reduction rate is at 25.1%. The average reduction rate is slightly smaller than when using the physical models for reduction (see Figure 7), but, overall, 36 out of 68 stations (52.9%) have a higher reduction rate when modeling the signal with meteorological parameters through a TCN. In Figure 13, the differences in RMS reductions between the TCN and physical model reduction are depicted. A positive difference means that the TCN-modeled signal gave a higher reduction. The results demonstrate that it is possible to reproduce comparable reductions through a deep learning model using only meteorological parameters as input features. The physical models are, on the one hand, complex to compute, but, on the other hand, once they are computed, the complete dataset can be introduced into the GNSS residual computation at once and on all stations. The approach via a deep learning network works with minimally processed data, but the training of the network is time consuming and needs a lot of training data. In the evaluation of the trained model, not all stations can be reduced, only those that are dedicated test stations.

4.2. TCN Modeling Based on GNSS Residuals Reduced by Loading Models

The results of the direct subtraction of environmental loadings and of a meteorological-modeled time series from the GNSS residuals yield a similar level of RMS reduction. Therefore, in a second approach, the TCN model is trained on the GNSS residuals already reduced by environmental loadings. The idea is to increase the RMS reduction, by learning the potentially mis-modeled or leftover parts related to meteorological features. Once the model is trained, the correction signal is modeled for each test station and then subtracted from the reduced GNSS time series.

Figure 14 shows the results of the TCN-modeled time series, at MAR3. Similar to the results of the previous section (Figure 11), the seasonal patterns are very strong in the modeled signal. The annual peak is here even more distinct. The patterns of the underlying GNSS residuals are also different, as the environmental loadings were already subtracted. An interesting observation concerns the amplitude of the annual signal. Despite a smaller RMS variation overall, the annual peak is higher compared to the original GNSS residuals. This is due to a stronger annual signal in the environmental loading, mainly introduced by HYDL. HYDL has a clear annual pattern, which is, at most stations, not entirely in phase with the GNSS residuals. Phase differences between HYDL and GNSS have been quantified and discussed in several studies [34,35,36]. The reasons found are diverse and mostly trace back to local phase-inhomogeneities of surface water and the absence of a model for groundwater storage, but no correction models are proposed. Only recently, Michel et al. [37] made an attempt to model and correct for a misinterpretation of horizontal fluxes, which leads to an improvement of the phase advance over the GNSS seasonal signal. In our work, phase deviations are not corrected; therefore, the annual signal is inflated for some stations. These seasonal peaks get reduced by the TCN-modeled signal to a lower level than the original GNSS residuals. The different stages of reduced GNSS time series are shown in Figure 15.

The additional RMS reduction of the GNSS residuals is computed, as well, and is depicted in Figure 16. The achieved additional RMS reduction is, on average, 2.7%. Overall, the RMS of 84% of the stations could be further reduced.

In a final step, all reductions are applied to the GNSS residuals, and, again, the RMS reduction rate is computed and shown in Figure 17. The average RMS reduction is now at 28.6%, with a maximum reduction of 44% at station KALL. Only one station is negatively affected, meaning that the RMS of the variations in the height time series increased after all reductions.

5. Conclusions and Outlook

The aim of this work was to decrease the RMS of GNSS height residual time series by correcting them for environmental influences. First, the residuals were corrected for environmental surface loadings, which reduced the RMS of the time series, on average, by 25.9%. These results were compared to the RMS reduction, achieved by reducing the residuals by a time series modeled by a temporal convolutional neural network (TCN). The TCN was thereby trained on 206 stations to learn the relationship between meteorological parameters (temperature, precipitation, sea level pressure, humidity, and radiation) and the GNSS height residuals. The mean RMS reduction, achieved by reducing the GNSS residuals by the TCN-modeled signal, reached 25.1%. Although the mean RMS reduction through the TCN is slightly smaller, the majority of stations (52.9%) have a higher RMS reduction when using the signal modeled by the TCN. This suggests that the physical models can be recreated to great extend through a machine learning approach.

In a second experiment, the TCN was trained on the time series already reduced by environmental loading models. The mean additional RMS reduction of the residual time series was 2.7%, suggesting a relationship of the remaining residuals and meteorological features. Adding up both reductions, by environmental loading models and the TCN-modeled signal, the overall average RMS reduction was 28.6%, where 98.5% of the test stations were positively reduced, demonstrating the strong potential of applying machine learning models for GNSS time series modeling.

This study can be extended by using gridded meteorological data from numerical weather models to include more available GNSS stations. The current model could be improved by further hyperparameter tuning in an extensive grid search and complemented by a comparison to different machine learning algorithms. Furthermore, the approach could be extended by testing other study regions, and, in a following step, the derivation of global models could be evaluated.