Robust-Extended Kalman Filter and Long Short-Term Memory Combination to Enhance the Quality of Single Point Positioning

In the recent years, multi-constellation and multi-frequency have improved the positioning precision in GNSS applications and significantly expanded the range of applications to new areas and services. However, the use of multiple signals presents advantages as well as disadvantages, since they may contain poor quality signals that negatively impact the position precision. The objective of this study is to improve the Single Point Positioning (SPP) accuracy using multi-GNSS data fusion. We propose the use of robust-Extended Kalman Filter (referred to as robust-EKF hereafter) to eliminate outliers. The robust-EKF used in the present work combines the Extended Kalman Filter with the Iterative ReWeighted Least Squares (IRWLS) and the Receiver Autonomous Integrity Monitoring (RAIM). The weight matrix in IRWLS is defined by the MM Estimation method which is a robust statistics approach for more efficient statistical data analysis with high breaking point. The RAIM algorithm is used to check the accuracy of the protection zone of the user. We apply the robust-EKF method along with the robust combination of GPS, Galileo and GLONASS data from ABMF base station, which significantly improves the position accuracy by about 84% compared to the non-robust data combination. ABMF station is a GNSS reception station managed by Météo-France in Guadeloupe . Thereafter, ABMF will refer to the acronym used to designate this station. Although robust-EKF demonstrates improvement in the position accuracy, its outputs might contain errors that are difficult to estimate. Therefore, an algorithm that can predetermine the error produced by robust-EKF is needed. For this purpose, the long short-term memory (LSTM) method is proposed as an adapted Deep Learning-Based approach. In this paper, LSTM is considered as a de-noising filter and the new method is proposed as a hybrid combination of robust-EKF and LSTM which is denoted rEKF-LSTM. The position precision greatly improves by about 95% compared to the non-robust combination of data from ABMF base station. In order to assess the rEKF-LSTM method, data from other base stations are tested. The position precision is enhanced by about 87%, 77% and 93% using the rEKF-LSTM compared to the non-robust combination of data from three other base stations AJAC, GRAC and LMMF in France, respectively.


Introduction
Nowadays, positioning services can greatly improve life quality by covering a wide range of applications such as automatic driving, intelligent transportation, agriculture and so on. In this This paper hence proposes to use LSTM as a denoising filter. Figure 1 summarizes the global flowchart for the position determination using the proposed approach. The remainder of this paper is organized as follows: Section 2 introduces the single point positioning technology; Section 3 presents the Robust-Extended Kalman Filter; Section 4 shows the experimental results of applying robust-EKF; Section 5 presents the de-noising filter method and experimental results; Section 6 presents the conclusions and perspectives.

Single Point Positioning Technology
The Single Point Positioning technology can provide meter-level positioning accuracy for emerging multi-GNSS integration and opens various new prospects. GNSS brings more signals to improve the user position and velocity accuracy. In this paper, observation models of three satellite systems (GPS, Galileo and GLONASS) are combined. For combining these satellites data, two important things are paid attention to: coordinate reference and time reference.
With respect to the coordinate reference, the coordinate systems of GPS, Galileo and GLONASS satellites adopt the broadcast orbits of WGS-84, GTRF and PZ90.11, respectively. Even though their coordinates are different, the disparity is around several centimeters as stated in [29]. As a result, their coordinate systems are considered the same.
As for the time reference, GPS, Galileo and GLONASS systems define their own time scale. These time systems are presented in the RINEX report version 3.03 [30]. GPS time (GPST) runs parallel to UTC (Universal Time Coordinated) (1 microsecond difference), but it is a continuous time scale that does not insert any leap seconds. Galileo runs on Galileo System Time (GST) which is nearly identical to GPS time (tens of nanoseconds difference). GLONASS is basically running on UTC or, more precisely, GLONASS time (GLON ASST) linked to UTC(SU). It is not a continuous time, i.e., it introduces the same leap seconds as UTC. The reported GLONASS time has the same hours as UTC and not UTC+3h as the original GLONASS System Time. Apart from the mild errors in the realization of the different time systems, the relation between them is given by: GLON ASST = UTC and GPST = GST = UTC + ∆tLS + τGPS (1) where ∆tLS is the time difference between GPST and UTC due to leap seconds (1 January 2019, ∆tLS = 18), τ GPS is the fractional time offset between the two systems (GPS and GLONASS). After synchronization of the coordinate and time systems, the position and velocity estimation models are presented in the following.
Pseudo-range measurements are considered as the sum of the distance to the satellite, the receiver clock error, the ionospheric and tropospheric errors and multi-paths. They are expressed as follows: ρ G,E,R = (x G,E,R − x u ) 2 + (y G,E,R − y u ) 2 + (z G,E,R − z u ) 2 + c (δt G,E,R − δTS G,E,R ) + Iono G,E,R + Tropo G,E,R + Mul pa G,E,R + ε ρ,G,E,R (2) where the superscripts/subscripts G, E and R refer to GPS, Galileo and GLONASS satellites respectively; ρ is the measured pseudo-range; c is the speed of light in vacuum; [x u , y u , z u ] is the user position; [x G,E,R , y G,E,R , z G,E,R ] are positions of GPS, Galileo and GLONASS satellites respectively; δt is the user clock offset; δTS is the satellite clock offset; Iono is the ionospheric delay; Tropo is the tropospheric delay; Mul pa are the multi-paths and ε ρ is measured noise.
Before the determination of user position and velocity, it is required to correct the GNSS errors that are the satelllites clocks bias, the ionospheric errors, the tropospheric errors, and the multi-paths. The corrections of the satellite clock for GPS and GLONASS are presented in [31] (pp. 88-89), whereas the correction of the satellite clock of GLONASS is discussed in [30] (p. 34). We use the Saastamoinen model to correct the tropospheric errors as in [32] (pp. 135-137). We correct the ionospheric errors using the dual-frequency method in [23] (p. 286). While the multi-paths are estimated by the combination of the code pseudo-ranges and the phase pseudo-ranges in [23] (pp. 290-291). After compensating for satellite clock bias, ionospheric errors, tropospheric errors, and multi-paths from Equation (2), the corrected pseudo-ranges can be written as: where ρ c is the corrected pseudo-range; r is the true range from the satellite to the user: Given the Doppler effect, the pseudo-range rateρ can be computed as [31]: where D is the Doppler effect; c is the speed of light and f is the satellite transmitted frequency. The velocity can be estimated from the pseudo-range rate, starting by differentiating Equation (3) to obtain:ρ where δṫ is the user's clock drift (sec/sec); ερ is the error in observation (meters/sec); andṙ is the true range rate expressed as: where v x,G,E,R , v y,G,E,R , v z,G,E,R is the satellite's velocity; v x,u , v y,u , v z,u is the true user's velocity; and 1 x , 1 y , 1 z is the true line of sight unit vector from the satellite to the user: Thus, the combined equations for the corrected pseudo-ranges (Equation (3)) and corrected pseudo-range rates (Equation (6)) can be as expressed as follows: To simultaneously estimate the user position and velocity, the state vector is defined from Equations (7)-(9) as follows: This vector will be used in robust extended Kalman filter algorithm detailed in Section 3.

Robust Extended Kalman Filter
Extended Kalman Filter (EKF) is an extended application of Kalman filter in solving nonlinear optimal filtering problems. It is used in a wide range of technological fields (radar, electronic vision, communication). However, in hybrid positioning, the large outliers can negatively impact its application to position estimation. According to [17,18], the resulting Kalman filter outliers can be sorted into three types: observation, innovation and structural outliers. In this paper, robust extended Kalman filter is used to handle observation and innovation outliers. Robust-EKF is presented as the combination of the Extended Kalman Filter with the Iterative Reweighted Least Squares algorithm (Section 3.3) with a weight matrix determined by MM-estimation (Section 3.2). Moreover, a combination of Robust-EKF with the RAIM algorithm (Section 3.4) is used to check the accuracy of user zone.

Robust Extended Kalman Filter Model
The state space representation for the extended Kalman filter model is defined by two equations as follows: State Equation : Observation Equation : where X k is the state vector and Z k is the measurement vector at time t k ; k is the vector that conveys the system error sources; e k is the vector that represents the measurement error sources; f () is the function for state transition and h() is the function for the measurement. IRWLS works on linear regression to handle the problem of outliers in the data. Therefore, converting a non-linear model into a linear model is necessary to apply the robust-EKF that consists in three steps: Linearizing equations, reforming filter and updating.
Step 1: Linearizing equations The linearized state in Equation (11) is described as: where A k is the Jacobian of the process model. The predicted state is determined by: The relation between the true state X K and its prediction X − k can be written as: where δ k is the error between the true state and its prediction. Linearizing Equation (12) using the first order Taylor series expansion with the predicted state vector X − k : By noting: ∂h (16) can be written as: Or Step 2: Reforming the filter This step will reform the filter to a regression equation. To perform this, combining Equations (15) and (18) together in a matrix form yields the following: where I is the identity matrix. Equation (19) is expressed in a compact form as: The covariance matrix ofȇ k is given by: where R k = E e k e T k is the measurement covariance matrix and P − k = E δ k δ T k is the covariance matrix of the predicted error. P − k is determined by: where Q k is the covariance matrix of the process. L k in Equation (21) is calculated using Cholesky decomposition. To perform pre-whitening, Equation (20) is multiplyed by L −1 k and by defining: (20) becomes the regression equation written as: Step 3: Updating This step consists in updating the estimated state and checking the position accuracy using RAIM algorithm presented in section "RAIM algorithm" (Section 3.4). The estimated stateX k in Equation (23) is optimized by the IRWLS algorithm presented in section "Iterative Reweighted Least Squares Algorithm" (Section 3.3). The final state estimateX k is given by: where W k is the diagonal weight matrix at time t k . This matrix is defined in section "MM Estimation theory" (Section 3.2) and is calculated in section "Iterative Reweighted Least Squares Algorithm" (Section 3.3).
The error covariance matrix is updated as follows: where K the gain of the filter given by:

MM Estimation Theory
In regression analysis, using the least squares method would not be enough to handle outliers or extreme observations. We therefore resort to the MM-Estimation which is known as a robust statistics method providing a more efficient analysis of statistical data with high breaking point as demonstrated in [20,22].
The linear equation in this paper is written as follows: where y is the new observation vector after combining system and observation data; G is the combined system and observation matrices corresponding to y; X is the unknown state vector and e is the noise vector corresponding to y. The robust method is used to bound the influence of outliers on the X state estimates. The goal of the robust estimator is to seek the best fitting criterion to the real observations in order to reduce the effect of abnormal data. Generally, a robust estimator minimizes the cost function J(r) given by: where y i is the ith observation, G i is the transformation vector for the ith observation, and y i − G i X is the residual error for the ith observation. µ(r) is the objective function that gives the contribution of each residual to the cost function J (r). The scaled estimation is used to minimize the weights out of the restrained observations dynamically: where s is the scale factor. Let ϕ be the derivative of the objective function for the unknown X, ϕ = µ | X . Differentiating the objective function in Equation (29) and setting the derivative to zero: Which is equivalent to: By introducing a weight w(ȓ s,i ) to consider the effects of each scaled residualȓ s,i = y i −G i X s , Equation (31) can be modified as: The solution of Equation (32) is a weighted least squares estimation. Nevertheless, the weights are determined by the scaled residuals, which rely upon the estimated coefficients, which themselves rely upon the weights. For this reason, the Iterative Re-Weighted Least Squares method is used and presented in section "Iterative Reweighted Least Squares Algorithm" (Section 3.3) to optimize the solution. The weight function plays an important role and some variants are presented in Table 1. Table 1. Weighting functions [22].

Name Weight Function
Huber In Table 1, r i is the ith observation residual. α is a tuning constant to optimize the estimation with a balance between a high breakdown point and efficiency of the estimator [22].

Iterative Reweighted Least Squares Algorithm
An Iterative Reweighted Least Squares Algorithm (IRWLS) is performed to eliminate the outliers. The steps for the robust estimation procedure are defined as follow: 1. Find an initial estimateX 0X 2. Estimate the vector for initial residuals of the observations r 0 : 3. Define the initial scale value s 0 [33] s 0 = 1.4826 median (|r 0 |) where "median" stands for the median function computed on the residuals vector r 0 .
4. Estimate the initial diagonal weight matrix by MM-Estimation.
where i = 1, . . . , n; n being the number of observations.

While (j is jth iteration)
(a) Update the value of the matrix G j at the jth iteration withX j−1 (b) Solve the estimated stateX j using the weighted least-squared method (c) Calculate the HPL (in Section 3.4) (h) Recalculate the diagonal weighted matrix using MM-Estimation: If not, they are rejected.
In this paper, we use the bi-Tukey weight function in Equations (36) and (40). The HPL is calculated using the RAIM algorithm presented in Section 3.4.

RAIM Algorithm
In this paper, RAIM algorithm is used to detect positioning errors exceeding the alert limit. One of RAIM's outputs is the HPL, that is defined as the radius of a circle in the horizontal plane, centered at the true position. The Horizontal Alert Limit (HAL) is the maximum allowed HPL, which means that the estimated positions are accepted when HPL ≤ H AL. Given the linear Equation (27), the estimated state vectorX is determined by IRWLS as explained in Section 3.3: where W is the W matrix last obtained matrix in Equation (40). The state error is calculated as the difference between the estimated state vectorX and the true state vector X : The estimated residualê is determined by computing the difference between the measured observation y and the estimated observation (GX): where: Note that S is a projector matrix and thus idempotent: S 2 = S. The a priori standard deviation of the pseudo-range is calculated as follows: where n is the number of observations and p is the number of parameters in the unknown state X.
The position error in the state error (δX) of Equation (42) is defined in the (x, y, z) coordinate system. In practice, a local (e, n, u: east, north, up) coordinate system is more suitable. To transform from (x, y, z) coordinates to (e, n, u) coordinates, the F T orthogonal transformation matrix is required [23]: where ϕ is the longitude and λ the latitude.
Here, only the parameters of coordinates are considered. For example, the state vector X is defined following Equation (10) with the first, third and fifth positions in relation to the three (x, y, z) coordinates. Therefore, the newG 0 matrix is obtained asG 0 = G(:, 1 : 2 : 5) and the position errors as δX 0 = δX(1 : 2 : 5). The M matrix related to the errors in easting (dE), northing (dN) and upping (dU) is calculated as follows: Since the horizontal error related to easting and northing is considered, a matrixM 0 can be defined asM 0 = M(1 : 2, :). In reality, none of the elements in e is equal to zero. However, when serious observation measure error occurs, some elements in e may be extremely larger than the others. Therefore, we only consider the ith observation that has a failure of magnitude β, while we suppose that the other elements are equal to zero. The error vector can thus be re-written as: The norm squared of δX EN (errors in easting and northing) assuming this specific choice of e is computed as: From Equation (43),ê = Se is recalled. Then S T S = S and the diagonal entry (i, i) of S is denoted s ii and we have: Combining Equations (49) and (50) gives: Or The obtained equation is a linear function with a straight line through the origin and a slope α i . The slope α i of the failure mode axis related to observation i is computed for all i = 1, . . . , n as: The mode axis with the largest α i is denoted α max and the HPL is defined as: where σ 0 is calculated following Equation (45).

Experimental Results of Applying Robust-EKF
Numerous investigations on robust-EKF for GNSS data exist in the literature [13][14][15]. In [13], robust-EKF is used for GEO/IGSO/GPS Raw-PPP/INS data fusion. In [14], it is applied for GPS/Galileo/GLONASS data combination. In [15], it is used for GPS/GLONASS data combination. All these studies assert the significant improvement of the position precision brought by the application of robust-EKF. In this paper, we combine the robust-EKF with the RAIM algorithm to ensure the position precision. The details of this robust-EKF and the experimental results for GNSS data are presented in the sequel in Sections 4.1 and 4.2 respectively.

Applying Robust-EKF
From Equation (10), we can write the compact user state of the three systems GPS, Galileo and GLONASS as follows: where [x u , y u , z u ] is the user position; [v x , v y , v z ] is the user velocity; b G , b E and b R are the errors in range due to the user's clock bias with GPS, Galileo and GLONASS time, respectively; d G , d E and d R are the user's clock drifts with GPS, Galileo and GLONASS time respectively.
We then proceed to the steps of the robust Extended Kalman Filter: Step 1: Linearizing equations By linearizing the function f (X k−1 ), the state transition matrix A k is obtained as a 6 × 6 block diagonal matrix of the form: where , and ∆T is the time interval between two epochs. The predicted state is then determined by: The relation between the true state (X k ) and its prediction (X − k ) can be written as: The observation equation is defined as: where Z ρ,G and Zρ ,G are the pseudo range (ρ G ) and the pseudo range rate (ρ G ) matrices of N gps GPS satellites tracked (dimension (N gps × 1) ); Z ρ,E and Zρ ,E are the pseudo range (ρ E ) and the pseudo range rate (ρ E ) matrices of N gal Galileo satellites tracked (dimension (N gal × 1)); Z ρ,R and Zρ ,R are the pseudo range (ρ R ) and the pseudo range rate (ρ R ) matrices of N glo GLONASS satellites tracked (dimension (N glo × 1) ). The ρ G ,ρ G , ρ E ,ρ E , ρ R , andρ R are calculated in Equation (9). Following step 1 in Section 3.1, the linearized "Observation equation" using the first order Taylor series expansion with the predicted state vector X − k is written as: where the matrix H k is determined as the derivative of the function h(X − k ) : H ρ,G , H ρ,E and H ρ,R are the Jacobians of the pseudo-range measurement model for GPS, Galileo and GLONASS data with the dimensions (N gps × 6), (N gal × 6) and (N glo × 6), respectively; Hρ ,G , Hρ ,E and Hρ ,R are the Jacobians of the pseudo-range rate measurement model for GPS, Galileo and GLONASS data with the dimensions (N gps × 6), (N gal × 6) and (N glo × 6), respectively. H ρ,G is computed as : And Hρ ,G is computed as: . . , N gps . H ρ,E , H ρ,R and Hρ ,E , Hρ ,R are computed similarly to H ρ,G , and Hρ ,G respectively.
Step 2: Reforming filter Following step 1 and step 2 in Section 3.1, from Equations (58) and (60) we have the equation as Equation (19). We compact it as Equation (20), the compact regression equation is as follows: The covariance matrix ofȇ is the measurement covariance matrix and P − k = E δ k δ k T is the covariance matrix of the predicted error, and both matrices are needed to determine L k . The determination of R k is not straightforward and in this work the Exponential model is proposed to estimate it. The Exponential model formulation [23] depends on the elevation angle and is given by: where l refers to the lth satellite; S l is the elevation uncertainty (m); x 0 and x 1 are constant; El(l) is the elevation angle of the lth satellite ( in degrees) and x 2 is a scaled value of the elevation error (in degrees).
x 0 , x 1 and x 2 must be estimated empirically. After performing many tests, we chose the values x 0 = 0.8, x 1 = 0.4 and x 2 = 12. The normalization elevation uncertainty σ 2 l is given by [7]: n s being the number of the tracked satellites. The covariance of the measurement noise of GPS satellite is determined as follows: where N gps is the number of the tracked GPS satellites. The covariance of the measurement noise of Galileo (R k,E ) and GLONASS (R k,R ) satellites are determined the same way as GPS satellite (R k,G ). Therefore, the covariance of the measurement noise of GPS, Galileo and GLONASS satellites is written as: where γR k,G , γR k,E and γR k,R are the covariance matrices of the pseudo-range rates of GPS, Galileo and GLONASS, respectively. γ is the covariance pseudo-range rate to the covariance pseudo-range ratio and is empirically chosen; After performing several tests, the value γ = 0.01 was chosen in this paper. The covariance prediction P − k is calculated by: where Q k is the process covariance matrix. It plays an important role in the optimal recursive process and is written as: The matrix ∑ is the spectral density matrix given by: where S x , S vx ; S y , S vy & S z , S vz are the power spectral densities of position noise and speed noise in the x-, y-and z-coordinates, respectively; S b is the power spectral density of the clock bias noise and S d is the power spectral density of the frequency drift noise. To simply, we consider to the ∑ which is unchangeable. The covariance matrix Q k in Equation (70) is rewritten as: With Q y and Q z are determined similarly to Q x ; and Q b,G is expressed as: Q b,E and Q b,R are determined similarly to Q b,G . After the calculation of the covariance matrix R k and the covariance prediction P − k , it is easy to determine the matrix L k using the Cholesky decomposition. Multiplying Equation (64) by L −1 k and defining: y k = L −1 kZ k , G k = L −1 kH k , ζ k = L −1 kȇ k , Equation (64) becomes the regression equation written as: Step 3: Updating The estimated stateX K in Equation (75) is optimized by the IRWLS algorithm previously detailed in Section 3.3. After testing all weight functions in Section 3.2, we chose the bi-Tukey weight function to eliminate the outliers data and α is set to 1.756 (empirically chosen after many tests) for the IRWLS algorithm. The final estimated state is calculated as Equation (24):X k = G T k W k G k −1 G T k W k y k . The error covariance matrix is updated as Equation (25): where K is the filter gain determined in Equation (26) : After the state estimation, RAIM algorithm in Section 3.4 is applied to check the position accuracy.

Experimental Results
To demonstrate the above algorithm, this study uses the GPS, Galileo and GLONASS data from ABMF station (RINEX version 2.11 (GPS and GLONASS) and version 3.02 (Galileo) format). This reference station is located in Guadeloupe. RINEX data and broadcast ephemeris data used for calculation were downloaded from RGP network and CDDIS for day of 1 January 2019. In order to assess the algorithm performance, five scenarios are investigated under MATLAB environment:  In Figure 3, the red curve pertains to scenario #1 with GPS data only. The black curve represents scenario #2 with Galileo data only. The cyan curve is related to scenario #3 with GLONASS data only. It can be seen that the position errors of GPS and Galileo data are small and stable whereas those from GLONASS data are large and unstable. The blue curve represents for scenario #4 with combination of GPS, Galileo, and GLONASS data. The position accuracy from the combination improves compared to the use of GLONASS data only. However, the errors still remain unsteady and large in the East axis. In particular, the positions are not determined around the times 22 h, 23 h, and 24 h, because of the large HPL (in Figure 4). Consequently, the robust-EKF method is applied to enhance the accuracy of user position in scenario #5 and the obtained position errors are represented in the green curve. They are more stable and more accurate than the other scenarios.     For more insights about the position accuracy obtained for the five scenarios, Table 2 presents a comparison of the corresponding position errors.  According to the RMS (Root Mean Square) error values listed in Table 2, the robust combination of GPS, Galileo, and GLONASS data improves the position accuracy by about 53 %, 41%, and 95% compared to the use of GPS data only, Galileo data only, and GLONASS data only, respectively. It is also improved by about 84% compared to the non-robust combination of GPS, Galileo, and GLONASS data.

RMS Error (m) GPS GAL GLO GPS/GAL/GLO Robust-EKF GPS/GAL/GLO
In the introduction section, it has been stated that the robust-EKF can be better than the Weighted Least-Squares Estimation (WLSE) to improve the position accuracy. To give tangible evidence on that, we compared the two methods and the results are illustrated in Figure 6. The Figure 6 presents the position errors obtained from the two methods using the combined GPS, Galileo and GLONASS data. This WLSE method was developed in our previous work in [8], where the weight matrix is defined using the "Exponential model". The position errors of WLSE are wavering; In particular, the positions around 23h-24h are unavailable because of high HPLs. On the other hand, the position errors of robust-EKF are stable and small.
Although the robust-EKF can greatly improve the position accuracy, its outputs might contain errors that cannot be estimated. Therefore, an algorithm that can predetermine the errors of robust-EKF behaviour is needed. For this purpose, a de-noising filter method is presented in Section 5.

De-Noising Filter Method and Experimental Results
The recent developments in the "Deep learning" field have allowed the use of its algorithms for GNSS applications [24,25,27,34]. In [34], an intelligent hybrid scheme consisting of an Artificial Neural Network (ANN) and Kalman Filter (KF) has been proposed to improve the accuracies of positional components as well as orientation components in real time. This scheme is able to overcome the limitations of KF. In [25], a "deep learning" algorithm was proposed to denoise the MEMS IMU output signals. The used algorithm is a Long Short Term Memory (LSTM) that was employed to filter the MEMS gyroscope outputs, in which the signals were treated as time series. LSTM is an artificial recurrent neural network (RNN) architecture that is useful for classifying, processing and making predictions based on time series data. Therefore, with a combination of the works presented in [25,34], LSTM is implemented in this paper to learn and compensate for the residual errors of robust-EKF in order to improve the position accuracy. Hence, the rEKF-LSTM hybrid schemes are proposed as a method capable of learning how the state vector behaves based on the dynamic of the filter. The details of rEKF-LSTM hybrid schemes and the experimental results for GNSS data are detailed in the upcoming Sections 5.1 and 5.2. Figure 7 presents the LSTM layered architecture. This design demonstrates the flow of a time series X (inputs) with N features of length S through an LSTM layer. The outputs are a time series "h" with H hidden units of length S. In the design, h k and c k denote the output (hidden state) and the "cell state" at time step t k , respectively. The first LSTM block uses the initial state of the network to calculate the first output (h 1 ) and the "updated cell state" (c 1 ) at the first time. At the t k , the block uses the previous state ( h k−1 , c k−1 ) to compute the output (h k ) and the "updated cell state" (c k ). A LSTM unit is composed of a cell, an "input gate", an "output gate" and a "forget gate". The basic structure of an LSTM unit is shown in Figure 8 illustrating the flow of data at instant t k . This diagram highlights the three parts: forget gate, input gate and output gate.  A LSTM layer has the learnable weights and bias: the input weights W lstm , the recurrent weights R lstm and the bias b lstm . The matrices W lstm , R lstm , and b lstm are the union of the input weights, the recurrent weights, and the bias of each part, respectively. These matrices are grouped as follows:

De-Noising Filter Model
where f , i, c, and o denote the forget gate, input gate, cell state, and output gate, respectively.
In each LSTM block, the objective is the calculation of the output (hidden state) and the cell state through three parts: forget gate, input gate and output gate. As shown in Figure 8, the "forget gate" is presented in the first part of the LSTM, which is used to decide what information will be kept in the cell state. The decision is made by a sigmoid layer called "forget gate layer". h k−1 and x k are inputs to the sigmoid function, and the output is a value ranging from 0 to 1 for each number in the cell state C k−1 . If the output is "1", the information is "completely kept" in the cell state. When the output is "0", the information is "completely cleared". The forget gate's activation vector f k is written as: where σ(.) is the sigmoid function, h k−1 is the hidden state at time t k−1 , and x k is the input vector at time t k . The second part is the "input gate", which is employed to decide which new information should enter the previous cell state. This gate is composed of two parts: (1) a "sigmoid" layer to decide what values will be updated. The output values i k for this layer ranging from 0 to 1. "0" means "not important" and "1" represents "important"; (2) another part is "tanh" layer which creates a vector of new candidate values C k between −1 and 1 to help regulate the network. The input gate's activation vector i k and the new candidate vector C k are calculated as follows: where "tanh" is the hyperbolic tangent function.
The new cell state C k is updated from the old cell state C k−1 as follows: In the last part, the "output gate" will decide the output. First, a sigmoid layer is used to determine what parts of the cell state will be output. After that, the cell state is passed through a "tanh" function and the ranging of the cell state values is from −1 to 1. Finally, the results are multiplied by the output of the sigmoid gate, and the output parts are decided. The output gate's activation vector and hidden state vector are given by: The above Equations (77)-(82) are the steps to compute the output and the cell state of a LSTM unit. LSTM network can be applied for classification or regression tasks with sequence and time series data. In this paper, we make use of it for the regression task. From Figure 7, we obtain outputs "h" with the dimensions (H × S), while we want have predicted results "Y" with dimesions (M × S); therefore we need a weight matrix with dimension (M × H) to tranforme "h" to "Y". As a result, we have the full architecture of the LSTM network with one LSTM layer for regression as the Figure 9. Applying this LSTM network, the inputs data (X) are the positions estimated by robust-EKF in Section 3, and the outputs data (Y) are the errors on the estimated positions.
For training and testing our model, we should divide our data into three distinct sets: training, validation and prediction. For the "training part", Figure 10 presents the rEKF-LSTM training architecture. We determine the estimated position by robust EKF using the combination of GPS, Galileo and GLONASS data. The ENU rEKF is the East, north, and up positions by robust EKF. Before the estimated positions enter the LSTM model, we need to standardize them to have zero mean and unit variance. As a result, the standardized estimated positions are now the training inputs of the LSTM model. The target outputs are the errors of the estimated positions by robust-EKF (δ ENU ). The errors in the estimated positions are the difference between the true and the estimated (robust-EFK) positions. Similarly, we standardize the target outputs to have zero mean and unit variance and the standardized target outputs are now the training target outputs of the LSTM model. The model is trained on the training dataset using supervised learning with optimization methods such as: stochastic gradient descent, Adam or RMSProp. Each time, the LSTM model is run with the training inputs and returns the training outputs, which are then compared with the targeted ones. Based on the results of the comparison, the parameters of the model (weights and biases) are adjusted. The training process works until the training errors reach the error threshold and the inferred LSTM model is referred to as the "fitting LSTM model". The subsequent part is the "validation part" for which the architecture is illustrated in Figure 11. Firstly, we standardize the estimated east, north, and up positions. After, we use the "fitting LSTM model" from the "training part" to predict the standardized ENU errors. Then, we unstandardize them to obtain the predicted ENU errors. Following that, we compare the predicted ENU errors with the true ENU errors. Based on the comparison, we can see how well the LSTM model is generalizing in the training part. If the model is overfitting or underfitting, it is necessary to go back to the training part to tune the model's hyperparameters (the number of layers, the number of hidden units in each layer etc.). We can return to the training part many times until we get the best LSTM model to ensure confidence on our model performance. After obtaining the best LSTM model, this model is used in the "prediction part". Figure 12 presents the prediction architecture of rEKF-LSTM. First, we standardize the estimated east, north, and up positions similarly to the validation part. After, we use the "best LSTM model" previously validated in the "validation part", to predict the standardization ENU errors. Then, we unstandardize them to obtain the predicted ENU errors. Finally, we obtain the estimated positions by rEKF-LSTM, which are computed as the difference between the estimated positions by robust-EKF and the predicted errors following Equation (83). The estimated positions by rEKF-LSTM are written as follows:

Experimental Results
In this section we aim at comparing the position estimation performance of the two methods : the robust-EKF and the rEKF-LSTM.
A total of 20 hidden layers were set up for the LSTM in this paper. The LSTM set up to H = 20 hidden layers, N = 3 inputs features (positions), S = 1440 number of timesteps and M = 3 outputs features (errors positions). The data used are still from ABMF station recorded on the 1 January 2019 (as in Section 4.2). Figure 13 represents the position errors with respect to time for the two methods: robust-EKF and rEKF-LSTM. The period of data is from 12 h to 24 h, since data from 0 h to 12 h is used to train the model.
In Figure 13, the green curve represents the position errors for the robust-EKF method, whereas the purple curve represents the errors produced by the rEKF-LSTM method. It can be seen that the position accuracy of rEKF-LSTM method largely improves and remains stable. To have more insights, Table 3 summarizes the RMS position errors of the two methods.
As shown in Table 3, the position accuracy improves by about 74.0 % using rEKF-LSTM compared to the robust EKF. To have a more general view, the five scenarios in Section 4.2 and two scenarios in this section are combined. As a result, six scenarios are considered since the scenario pertaining to robust-EKF is in common. The underlying results are summarized in Table 4.   Table 5 gives the RMS errors in position for three other base stations: AJAC, GRAC and LMMF. The position precision is enhanced by about 87%, 77% and 93% using the rEKF-LSTM compared to non-robust data combination from the three base stations AJAC, GRAC and LMMF in France on 1 January 2019, respectively. This corroborates the conclusion that the rEKF-LSTM method can significantly improve the position precision.

Conclusions and Perspectives
In this paper, the SPP technology is applied to determine the user position calculated by GPS, Galileo and GLONASS data fusion. Data fusion may be containing outliers which negatively impact the final results. To address this issue, the MM-Estimation method is used to eliminate the outliers and IRWLS method is used to optimize the estimated results. Although the MM-Estimation and the IRWLS are effective, a method to detect the positioning errors is needed. The RAIM algorithm was used for this purpose to detect positioning errors exceeding the alert limit. As a result, the proposed robust Extended Kalman Filter is built by combining the Extended Kalman Filter, the MM-Estimation, the IRWLS and the RAIM method to improve the position accuracy. Robust GPS, Galileo and GLONASS data fusion using the robust-EKF increases the position precision by about 84.0% compared to non-robust estimation using data of ABMF base station. Furthermore, the position precision is significantly enhanced using the rEKF-LSTM method by about 95% compared to non-robust data combination. To prove the power of rEKF-LSTM, the approach is used for data from three other base stations: AJAC, GRAC and LMMF. The position accuracy is significantly improved by about 87%, 77% and 93% compared to the non-robust data combination, respectively.
In this work, the emphasis is laid on the improvement of the positioning accuracy for a slowly moving receiver. Only data from base stations are used, for which velocities are almost equal to zero and helps enhancing the performance of the LSTM model. A perspective of our work could be assessing the performance and limitations of the method on data from a user in motion and adapting or potentially developing a new method taking into account additional velocity information. While considering data from moving users, sensor noise levels are not constant and the challenge is to estimate an adapted LSTM modeling. Moreover, an additional challenge concerns the GNSS signal blockages. A suggested solution is to combine GNSS and INS (Inertial Navigation System) data. Along these lines, the perspective for this work is to use an adapted rEKF-LSTM model for combining GNSS and INS data from users in motion.