Gauss process regression for real-time ionospheric delay estimation from GNSS observations

The number of devices equipped with global satellite positioning has exceeded seven billion recently. There are a wide variety of receivers regarding their accuracy and reliability. Low cost, multi-frequency units have been released on the market latterly; however, the number of single-frequency receivers is still significant. Since their measurements are influenced by ionospheric delay, accurate ionosphere models are of utmost importance to reduce the effect. This paper summarizes how Gauss process regression (GPR) can be applied to derive near real-time regional ionosphere models using raw Global Navigation Satellite System (GNSS) observations of permanent stations. While Gauss process is widely used in machine learning, GPR is a nonparametric, Bayesian approach to regression. GPR has several benefits for ionosphere monitoring since it is quite robust and efficient to derive a grid model from data available in irregular set of ionospheric pierce points. The corresponding instrumental delays are estimated by a parallel Kalman filter. The presented algorithm can be applied near real-time, however the results are offline calculated and are compared to two high quality TEC map products. Based on the analysis, the accuracy of the GPR modell is in 2 TECu range. The developed methods could be efficiently applied in the field of autonomous vehicle navigation with meeting both accuracy and integrity requirements.


Introduction
The ionosphere is the upper part of the Earth's atmosphere with an altitude between about 40 and 2000 km. The ionosphere is ionized by solar radiation which means that there are free electrons and positive ions in this region Schaer (1999). In addition to playing and important role in the Earth's magnetic field, the ionosphere influences radio signal propagation from satellites to the Earth.
One of the major parts of the error budget in the GNSS (Global Navigation Satellite System) measurements is the ionospheric delay. The effect depends on the number of charged particles, the Total Electron Content (TEC) along the path of the signal through the ionosphere. With multi-frequency devices, the delay can be removed by the so-called ionospheric-free linear combination Hofmann-Wellenhof et al. (2008). However, still major part of the devices operates on a single-frequency; thus the ionosphere modelling is crucial to mitigate its effect.
One of the earliest and simplest models is the Klubochar model Klobuchar (1987). The corresponding parameters are broadcast by the GPS satellites among navigational messages. In general, Klobuchar model can describe just a limited part of the whole effect (Feess and Stephens 1987;Filjar et al. 2009), however in extreme conditions, during ionospheric storms for instance, its accuracy can be degraded even more (Bergeot et al. 2010;Gordienko et al. 2005;Hu et al. 2005). The most recent satellite system, the European Galielo, broadcasts its own ionosphere model, the NeQuick one Angrisano et al. (2013). According to a throughout investigation and comparison of Klobuchar and NeQuick model Farah (2008), it can be concluded in general that NeQuick model offers a better behaviour when the ionosphere is stable but slightly poorer behaviour with higher variability of ionosphere or close to maximum TEC values.
High quality TEC maps can be derived from the measurements of permanent station with post-processing Schaer (1999). The spherical harmonics and B-spline based TEC models are adequate for global and regional representation, while the polynomial based models are limited for regional modelling. These models describe the TEC of the ionosphere as a function of time and space Hernández-Pajares et al. (2009). The commonly accepted exchange format is the IONEX Schaer et al. (1998), that represents the 2-and 3-dimensional TEC maps given in a geographic grid. Thanks to the rapid development of computer processing, providing global and regional tec maps can be done in near-real-time (Coster et al. 1992;Bergeot et al. 2014;Erdogan et al. 2017;Chang et al. 2019;Renga et al. 2018).
The Gauss process regression (GPR) method is well-known in the machine learning community and proved its utility in various domain. It is also worth to mention that the mathematical background of the GPR and popular Kriging in geostatistics are the same Stein (1999). Data from a single station in India was used to forecast ionospheric delay using GPR technique with almost 100% accuracy Lakshmi et al. (2020). Another GPR approach was used to predict daily TEC values based upon the TEC values recorded at various permanent GNSS station in Turkey Inyurt et al. (2020). A comparison of Gaussian process (GP) and neural network model was performed by Ackermann et al. (2011) over a test area in South Africa. The GP framework presented many advantages over competing modeling strategies, such as providing powerful and convenient ways of incorporating prior knowledge and requiring less training data than neural networks. Another recent study shows the GP ability to enhance the positioning performance by improved TEC estimation in real-time Lin et al. (2019).
The safety critical applications demand for high reliability besides accuracy. Auxiliary services are preferred to achieve this high standard of requirements. For the European region, the EGNOS (European Geostationary Navigation Overlay Service) grants real-time corrections from the measurements of RIMS (Ranging and Integrity Monitoring Stations) Ventura-Travest and Flament (2007). The stations are located across the European region in a sparse network (Fig. 1). The very first initiative to monitor the achievable accuracy and 1 3 integrity with EGNOS was performed by the EGNOS Data Collection Network (EDCN) project Soley et al. (2004). EGNOS performance can be slightly degraded at the boundary of coverage area due to lack of reliable ionosphere modelling over the periphery Grunwald et al. (2016) or Lupsic and Takács (2019). This paper presents a GPR based novel ionospheric model for real-time TEC map generation from RIMS data available in the European region. The accuracy of the developed model is investigated by comparing to widely used and referred global and regional models developed using raw GNSS observations.

Algorithm outline
A brief overview about the developed algorithm is presented in Fig. 2 and the mathematical details of each block will be featured in the following sections. A regional receiver network, in our case the EGNOS system provides code and phase measurements in double frequencies. Combined measurements cancel out the geometry dependency and the result directly reflects the state of the ionosphere and the combined instrumental delays of receiver and the tracked satellite. This is the so called L 4 combination and the rigorious mathematical description of its calculation is presented in Sect. 2.1. In case of known instrumental delays, the slant ionospheric delays from the L 4 observation can be calculated. A single layer model determines the conversion from the slant ionospheric delay to vertical ionospheric delay alongside the location of the ionospheric pierce point (IPP) [Sec. 2.1]. The vertical ionospheric delay with their corresponding locations are the input to the Gauss process regression model which is described in Sect. 2.2. The GPR model creates the Total Electron Content map and updates it in real time. The presented algorithm in the paper is designed to be capable of run in real time, therefore the instrumental delays shall be estimated, and monitored continuously in parallel with the GPR. For that purpose, a polynomial model was used to describe the vertical total electron content (VTEC) and the data was processed iterative with a Kalman filter. The Sect. 2.4 is dedicated to present this part of the algorithm. It is worth to mention any other near-real time estimation could be used which provides an adequate accuracy for instrumental delays.

Ionospheric TEC measurements from GNSS
The ionospheric refraction is frequency-dependent, therefore with a dual-frequency receiver this delay can be eliminated in first order by the ionospheric-free combination. On the other hand, with the geometry-free combination, one could extract information about the state of the ionosphere Ciraolo et al. (2007). The combination can be created from the pseudorange measurements and the phase measurements. The noise of phase observation is typically two or three orders smaller than the code but suffers from the initial ambiguity. This ambiguity can be resolved by quasi ionosphere-free (QIF) (Teunissen 1995;Zhang et al. 2016).
The ionospheric delay depends on the frequency and the total electron content along the signal path between satellite and receiver.
where N e is the electron density in el∕m 2 dimension, which varies along the ray path. The vertical total electron content (VTEC) definition follows the form of Eq. (1) but the integration path is vertical.
(1) STEC = ∫  Converting the ionospheric delay from one frequency to another can be easily done with the Eq. (2), where j is the target and i is the reference frequency. I s r,i is the delay on the reference signal. The code and phase observation can be divided to frequency dependent and independent parts.
where j indicates the carrier frequency, R s r,j , s r,j code and phase measurement respectively, is the real geometry distance between the r receiver and s satellite, dt r is the receiver clock bias, dt s stands for the satellite clock bias, T s r is the tropospheric delay, I s r,j the ionospheric delay in j frequency, d r,j , r,j are the hardware bias of the receiver in j frequency of the code and phase observation respectively, d s j , s j are for the hardware bias of the satellite in j frequency of the code and phase measurement respectively, s r,j , s r,j denote stochastic, Gaussian type noises, N r,j is the integer phase ambiguity, j the wavelength of j frequency, s r,j , s r,j denote the multipath. Subtracting the code and phase measurement in frequnecy 1, 2 yields a geometry-free combination.
where is a frequency dependent constant, b r , b s are the differential code biases, s r,P , s r,P are the combined noise, and multipath delay respectively. The same combination can be created from the phase measurement as well, where C s arc,r is the ambiguity bias, B r , B s are the receiver and satellite interfrequency bias (ISB), M s r, , s r, are the combined multipath delay and noise. The significantly lower noise level of the phase observation can be exploited in the L 4 combination. The CPB s r (Carrier Phase Bias) is the offset between the s r,I and R s r,I combination and it can be estimated averaging out the geomatry-free code and phase differences in a cycle slip free continuous arc.
where N is the number of observation in the continuously observed arc. The leveled geometry-free combination, L s r,4 for a given arc is To achieve a near real-time ionosphere estimation the levelling method is not quite suitable Xiang et al. (2017), instead of that a Hatch-filter smoothing was applied. The smoothed geometry-free code R s r,I (n) can be calculated as: where, n = k when k < N and n = N when k ≥ N . The N value was choosed to 50 in this paper, but to avoid the nonconvergent smoothed code R s r,I (n) create a significant bias to the estimation, the first 10 samples in each new arc were not used for TEC modelling. A cycle slip detection is substantial to avoid abrupt biases in L 4 smoothing and in case of this event the Hatch-filter algorithm resets.
Once the L4 is calculated and smoothed, one can create a direct link with corresponding STEC value and the satellite and receiver biases.
where the , b r , b s represents the same terms from Eq. (5) and L 4 is the noise of the smoothed L s r,4 measurement. The vertical total electron content is one of the main properties of the ionosphere. One of the most used, and simplest method to establish connection between the slant and vertical TEC is a single layer model (SLM).
where STEC is the total electron content between a receiver and a satellite, m(z) is the elevation dependent mapping function. Fig. 3. depicts the concept of SLM, where the the vertical dimension of ionosphere is reduced to a single layer. The slant total electron content is considered at the ionospheric point, where the line of sight crosses the single layer. The m(z) mapping function converts the STEC to VTEC value therefore the different line of sight measurements are comparable. The latitude and longitude of the sub-ionospheric point are the characteristic properties of the IPPs alongside with its VTEC value. Substituting Eq. (11) to Eq. (10) yields the following: In case of known b r and b s the VTEC value can be calculated by rearranging the Eq. (12).
The latitude and longitude of the IPP can be calculated from the single layer model and from the known position of the receiver and satellite. These are the input set to the GPR modell, and in case of RIMS stations 200-400 separate IPPs with their VTEC values can be calculated in each epoch (Fig. 4). The next section will present the details how one can create a consistent TEC map from the fore-mentioned snapshot data with Gauss-process regression.

Gauss process regression
The central understanding of Gaussian process lays in the multivariate normal distribution. A random vector y = (y 1 , … , y l ) has the multivariate normal distirbution if it is derived from the following transformation of z = (z 1 , … , z M ) , where each z m ∼ N(0, 1) Blitzstein and Hwang (2019): where A is a K × M matrix and ∈ ℝ K . This represented as y ∼ N( , ) , where = Cov(y, y) = A T . The mean vector represented by , and the covariance matrix is the . Therefore the probabilistic functions of y is Consider splitting a multivariate normal variable y into (y 1 , y 2 ) two sub-vectors. Each of it has multivariate normal distribution hence we can split the A matrix above accordingly. Rewriting the y vector in terms of its sub-vectors gives The conditional distribution of y 2 given y 1 is P(y 1 |y 2 ) = N( 1|2 , 1|2 ) , where Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution. Formally, an f function is a Gaussian process if any finite set of values f (x 1 ), ...f (x N ) has a multivariate normal distribution on any S domain, x 1 , ...x N ∈ S . GP is specifed by a mean function m(x) and a covariance function k(x, x ′ ) , which is denoted as f ∼ GP(m(x), k(x, x ′ )) . The mean and covariance of f GP for any . The goal is to estimate function values of GP conditioned on some training data. Denote the set of inputs as X ∈ ℝ N×D , and the corresponding function values as f ∈ ℝ N , where D is the domain dimension. In the simplest case, the mean function is assumed as constant zero, the f values are noise-free. The k(x, x ′ ) covariance function is driven by its hyperparameters, or kernel parameters. One of the most popular kernel function is the squared exponential given by the following formula: where 2 f and l are the hyperparameters, denoted by above. To estimate the function values f * for a new set of inputs X * similarly to Eq. (16) assume .
The K XX , K XX * , K X * X , K X * X * covariance matrices are constructed from the kernel function.
The conditional distribution of f * can be calculated from the conditional properties of multivariate normal.
This derived formula is called the predictive distribution and it models the distribution of an unobserved set of inputs. The function values can be estimated by taking the mean, and furthermore the covariance structure of the estimated values are generally known. In practice the function values are not accessed directly but can be observed with some noise.
where n ∼ N(0, 2 ) . The observation noise can be incorporated into the kernel function by adding 2 to every diagonal term in each covariance matrix. The new kernel takes the form of where I(⋅) represents the indicator function. The predicted distribution of the unobserved y * values of the X * set of input can be derived similarly as Eq. (21).

TEC map estimation with Gauss process regression
In the previous section the model was considered with a zero mean function, however, in many cases a trend function defines the relation between the input and output data. Consider a training sata set {(x i , y i ; i = 1, 2, … , n)} , with a following linear connection: where x i ∈ ℝ d and y i ∈ ℝ , ∼ N(0, 2 ) . The h(x) is a chosen base function and the coefficients and the standard deviation are estimated from the training data set. In case of ionosphere modelling, the where f (x) is a zero mean GP with covariance function k(x, x � ) , f (x) ∼ GP(0, k(x, x � )) . The basis functions h(x) transforms the original x vector to a ℝ p feature space. The represents the coefficients of the basis vector with its p-by-1 dimension. The GPR model for an instance of response y is the following: The GPR model is nonparametric, hence it introduces an f (x i ) latent variable for each x i observation. The vector form of the model is where The basis matrix H has a dimension of n × p , where n is the number of observations, or training data, and p is the dimension of feature space. In case of the missing H matrix, the form converts back to the Eq. (22). The simplest choice for basis function is constant, so the feature space p would be one dimension, so the H = 1 basis matrix is an n × 1 vector of ones. In this paper the constant base matrix was used. Worth to mention the linear and quadratic form, where H lin = [1, X] and H quad = [1, X, X 2 ] respectively. Extracted form of the matrices in details in case of ionosphere modelling are: The f latent variable in the Eqs. (29) and (28) has a joint distribution by definition, and it follows the form of where the K(X, X) has the form of The used k(x, x � ) covariance function in this paper is the Matern52 kernel defined as where r is the Euclidean distance between x i and x j . The hyperparameters in the above mentioned kernel are f and l . The k(x i , x j | ) form highlights the dependency of the k(x i , x j ) kernel from the hyperparameters. To create a TEC map, first the basis function coefficients, the 2 noise variance and the kernel's hyperparamaters shall be estimated from the given (X, y) training data set. The used approach is the P(y|X) likelihood maximization as a function of , and 2 . The estimation is conducted by maximization of the marginal log likelihood function: The marginal log likelihood function has a form of: where First maximize the log likelihood function with respect to for given and 2 . The ̂ ( , 2 ) estimation is Substituting Eq. (36) to Eq. (34) yields a -profiled log likelihood.
Equation (37) is maximized over and 2 to find their estimates, then back substitute the optimized and 2 to the Eq. (36) to fix the coefficients. The last step is to make prediction to the new y g TEC values in the given X g data set, to create a TEC map. The observation noise free prediction follows as where and y g , X g , H g have similar from as in Eq. (29) but with the new X g = (x

Differential code bias estimation with Kalman filter
Kalman filter is a common choice for real-time estimation. Consider extending the Eq. (12) where the following m(z) mapping function of the used single layer model Dach et al. (2007): (33) ,̂ ,̂2 = arg max , , 2 log P(y|X, , , 2 ).
(38) P(y g |y, X, X g ) = N(y g |H X g + , ), (39) = K(X g , X g ) − K(X g , X)(K(X, X) + 2 I) −1 K(X, X g ), (40) E(y g |y, X, X g ), , , 2 ) = (H X g + K(X, X g | ) ,  where The a j n , n, j ∈ 0, … 5 are the coefficients estimated by the Kalman filter are the same declared as in Eq. (25). The j 0 , j 0 , j ∈ 0, … 5 are the center of the polynomials. The ( , ) is a weight factor. The inverse distance of the IPP and the center of polynomial is a common choice for the weight factor but in this paper the closest polynomial is the only contributor.
is the Euclidean distance of the ( , ) IPP coordinates and the origo of the j-th polynomial.
The first four polynomials cover the European region. The number and position of them are arbitrary and summarized in Table 1. During the model alignment, it was found that even one polynomial is able to follow the main trend of the ionosphere which is suffice for DCB estimation, however usage of four polynomials has gained better and more consistent results. The number of six, and nine were also tested, but due to the non-uniform spatial distribution of IPPs, some polynomials suffered the lack of nearby measurement points and struggled to converge, causing biases in DCBs. The 1 3 4th and 5th polynomial are dedicated to HBK and MON station respectively. They are located in North-America and South-Africa, therefore their observations cannot contribute to the TEC calculation of the European region but served well for DCB estimation of the satellites. The Kalman filter has a prediction step and update step. The estimated x state vector contains the a i n polynomial coefficients and the DCBs of satellites and permanent stations. The equations follow as where F k is the transition matrix, which in this model corresponds to a simple identity matrix. The y k is the observation vector, which contains the L 4 combinations. The D k design matrix transforms the x k unknown parameters from the state space to the observation space according to Eqs. (12) and (43).
The k is the process noise vector and e k is the measurement noise vector. Both of them are assumed as a white noise with zero expected values and the two vector mutually independent of each other.
where k,l is the Kronecker delta. The covariance matrix of the process noise consist three different sub-blocks, poly for the polynomial coefficients, sv for satllites' DCBs and stat for stations' DCBs.
where sv and stat have 2 sv and 2 stat in diagonal and zeros in off-diagonal. The poly has a form of: where poly i ∈ i = [0, … , 5] have only diagonal elements but differ in respect of coefficients. The GNSS observation was considered An elevation dependent wighting scheme was applied on the y observation covariance matrix. The non-diagonal elements were set to zero, and the i y diagonal variances were calculated by the following formula (adapted from Wang et al. 2015): where z i is the zenith angle corresponding to the ith IPP, and the 2 comes from the likelihood maximization presented in Eq. (33).

Concept validation with simulation
This section presents the GPR's capabilities for TEC map generation from the available VTEC values at the corresponding coordinates of IPPs. The assessment consists a simulation of the VTEC values in the IPPs for a given epoch from a known reference map. The IPPs position have been calculated from the postion of RIMS and GPS satellites with an ionospheric layer height 350 km. The used reference map is the product of Royal Observatory of Belgium (ROB) Bergeot et al. (2014) at date of 2020-01-23. The published ROB map has a 15 minutes time update and 0.5 • × 0.5 • spatial resolution in lattidue and longitude. The VTEC values come from the ROB product, then a zero mean Gaussion noise scaled by the obliquity factor is added to the simulated data. As first step, the L s 4 combination was calculated from the reference ROB map.
where L s 4 is the simulated measurement, s The VTEC s values with the corresponding IPP coordinates are collected for a given epoch and feeded to the GPR model to estimate a TEC map at matching grid points of the ROB map. A day long data set was used, therefore 24 × 4 TEC map were estimated with 2.5 • × 2.5 • spatial resolution, altogether more then 20000 points. After differentiation with the reference VTEC values, the discrepancies have been evaulated. The standard deviation of the added noise varied from 0 to 10 TECu. In noise free case, an almost perfect match is expected with the reference map. Besides the GPR, two additional methods were used for interpolation to give more perspective to the results. The assisting methods are the following: radial base function (RBF) with multiquadratic function, and third degree polynomial fit (POLY), adapted from (Yilmaz et al. 2009;Yu et al. 2015). The 1 3 method with zero and s = 6[TECu] added noise. The Fig. 5 contains altogether six subplots. The Fig. 5a, b for the GPR, Fig. 5c, d for the RBF and Fig. 5e, f for the POLY method in case of zero and 6 TECu simulation noise. The Gaussian process regression model demonstrates its capabilaties to effectively estimates TEC map from a snapshot measurement dataset. In error free case scenario, the GPR has better characteristics than the RBF or polynomial fit method (Fig. 5.a). In case of large noise the polynomial fit shows inferior quality than the other two examined methods (Fig. 5f). In Fig. 6, the global mean of |d k n | showed in respect of the added noise. The error bars stand for the global standard deviation of d k n . In both quality indicators, the GPR shows better performance in overall domain of the increasing measurement noise Hhowever However the performance of the polynomial and the radial base function could be definitely improved with better parameter adjustments strategies. This simulation has no goal to state a comprehensive qaulity assessment of the three chosen estimation method, only to show GPR's ability to create regional TEC maps from epoch-wise observation data of a sparse monitoring network.

Real-time ionosphere map creation from RIMS observations
The performance of the GPR model was investigated in nine days from 2020, grouped by three consecutive days. The first three days are Jan 23-25, when the sunspot numbers reached the lowest in the current 11 year long cycle McIntosh et al. (2020). The next three days are April 19-21. They were chosen because of the observed high ionospheric activity on April 20. And the last three days are July 19-21. The EGNOS RIMS observation data were used to run the Kalman filter aided GPR model. The created TEC maps are compared to the ROB and CODE products, and they were regarded as reference. To give perspective of the derived quality indicators, the same comparison was also performed to the two reference ionospheric map. Three type of VTEC differences are derived based Eq. (55).
where m i , m j ∈ [ROB, CODE, GPR] indicates the type of TEC map.
The global mean, the absolute mean and the standard deviation of d k n values are derived from Eq. (55) and showed in Fig. 7a, b, c respectively for each day and combination. The 6 Global mean of absolute differences in respect of simulation noise mean absolute values of ROB-CODE combination are close to 0.5 TECu for each day, meanwhile the GPR combinations with the references have 1 TECu absolute mean. These discrepancies are consistent in the analyzed time intervals. The standard deviation of ROB-CODE are varies around 0.6-0.7 TECu. The standard deviations of GPR-ROB and GPR-CODE are less than twice the ROB-CODE and approximately 0.9-1.2 TECu. The global means of ROB-CODE are also closer to the expected zero than the GPR counterparts. Spatial characteristics of mean and standard deviations are depicted in Figs. 8, 9 respectively. There is a visible trend in the mean values. The GPR method tends to overestimate the ionospheric delay in the northern region, hence the values are negative. On the other hand, the GPR underestimates the TEC values in the Mediterranean region. This skew does not appear when the two reference maps are analyzed in Fig. 8c. The mean absolute discrepancies of GPR are less than 2 TECu in all grid points. The gridwise standard deviations are less in the terrestrial regions, and shows correlation with the means both in case of ROB-GPR and CODE-GPR.
The two-dimensional histograms help to visualize the distribution of TEC differences in respect of the estimated values (Fig. 10). The plotted data are collected from all 9 investigated day. One can observe again a skewness in the GPR data in the lower TEC region (Fig. 10a, b). The outliers are contained in range of ±5 TECu. This distortion is not visible when the differencies of reference maps are depicted (Fig. 10c).

Fig. 7
Daily statistics between the ROB-GPR, CODE-GPR and ROB-CODE products depicted by blue, red and yellow lines respectively. The upper subfigures refer to data collected from January, the middle subfigures' data are from April, and the lower subfigures' data are from July To gain better insight, six grid points were chosen to show the daily variation of the VTEC values (Fig. 11). The shaded area represents the 2 standard deviation. There is no shading in case of ROB product because of the lack of available variance information. It is noticeable that GPR estimates correctly the TEC variation in each case. Mark the bias at lower and higher latitude (Fig. 11a, b and e, f). This flaw in the estimation was already visible in Fig. 8, and it is present during the whole day, and not the result of a temporally localized disturbance. The authors assume that the main source of this systematic error is the insufficient DCB estimation of the RIMS receivers.

Conclusion
The presented Gaussian process regression approach is a novel and promising method for ionospheric model derivation from multi-frequency GNSS measurements. It is capable to accurately estimate regional Total Electron Content (TEC) maps from snapshot measurements of a relatively sparse monitoring station. Accuracy of the GPR-based models in this paper is about ±2 TECu, which is comparable to the ±1 TECu accuracy of the corresponding models in the literature. In addition to the TEC modelling, the hardware delays were also estimated by a Kalman filter continuously and independently of the GPR. In further research, this loosely coupled setup could be tightened with direct DCBs estimation by the GPR. The hyperparameters of the GPR are calculated   Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.