Cooperative Localization under Ionospheric Scintillation Events

Ionospheric scintillation causes major impairments to Global Navigation Satellite System (GNSS) in low-latitude regions. In severe scenarios, this event can lead to complete loss of lock, thus making GNSS measurements unusable for navigation. In this paper, we derive a cooperative localization algorithm where a set of partially connected aircraft exchange messages with neighboring nodes on the network to improve their own position estimates. We consider the scintillation events as abrupt changes in the measurement variance, which are modeled by a discrete-valued Markov process at the nodes which have access to GNSS measurements. Simulation results show that Markovian modeling and cooperation via factor graph message passing reduce the average 3D root mean square localization error and yield an average vertical position error that meets civil aviation standards for approach and landing.


Introduction
Global Navigation Satellite System (GNSS) is currently used in a wide range of applications, including transportation, surveillance, agriculture and meteorology [1][2][3][4][5][6][7][8][9].The system is composed of multiple satellite constellations which provide accurate and reliable positioning, navigation and timing information worldwide [10].Despite its benefits, in low-latitude regions the GNSS solution is frequently degraded by ionospheric scintillation [11][12][13][14].This phenomenon causes rapid fluctuations in the amplitude and phase of GNSS signals, which in turn can make the measurements unusable for positioning in critical applications such as civil aviation [15].
In the last decades, substantial effort has been made to deal with GNSS degradation caused by scintillation events.The use of multiple-frequency receivers is currently the most common approach to improve the positioning accuracy of faded signals [16][17][18].However, these sensors are nowadays undergoing certification in order to meet the strict civil aviation operational standards [19].Machine learning (ML) algorithms constitute another prevalent solution to predict scintillation events and improve the positioning performance [20][21][22][23].However, this approach often relies on huge datasets composed of past data for training the ML models.Advanced integration schemes are also complimentary solutions for improving positioning accuracy based on sensor fusion techniques [24][25][26].However, these schemes require the use of additional hardware devices, such as pseudolites or high-resolution cameras, further increasing the costs.
In this paper, we aim at mitigating the scintillation effects and improving positioning and navigation of a set of partially connected aircraft within a wireless communication network by means of cooperative localization [27][28][29][30][31][32][33]56].To that end, as in [34,35], we model this distributed estimation problem in a Bayesian fashion following the factor graph framework to apply the sum-product algorithm (SPA) (introduced in [36,37]) and solve for the marginal posterior distributions of the unknown state variables at each network aircraft, also referred in this paper as a network node.
In addition, we introduce some a priori information and make use of discrete-time Markov processes to model the scintillation occurrences at each node.Following this approach, it is possible to predict when a scintillation event is likely to occur and adjust the cooperative localization algorithm accordingly.
We follow the lead in [31] and present an efficient implementation of the message passing scheme based on Gaussian approximations that reduce the internode communication overhead.We also present a sequential Monte Carlo (SMC) filtering procedure [38,39] that runs the SPA in an efficient way at each network node and computes the approximate marginal posteriors of interest.By following this approach, the information is diffused over the network by building Monte Carlo representations of the marginal posteriors at each node conditioned on the available measurements over the entire network up to the current time step.
This paper is organized as follows.Section 1 is this Introduction.In Sect. 2 we review the main concepts from Factor graph theory.Section 3 presents the state model of the nodes and the measurement models.Section 4 presents the derived message passing scheme for the proposed factor graph and sum-product algorithm messages.In Sect.5, we introduce the Gaussian-SMC implementation of the messages.Section 6 shows the numerical simulation setup and results.Finally, Sect.7 summarizes the present work with conclusions.
Notation In this paper, we do not distinguish the notations for random variables/ vectors and real variables/ vectors, including samples of random variables or vectors.The proper interpretation is implied in context or otherwise explicitly stated, when required.We make use, however, of boldfaced fonts to denote vectors and matrices, and italic fonts to denote scalar objects.The notation p(x) denotes both prob- ability density functions (p.d.f.) and mixed probability density functions, which are defined as probability density functions in the continuous variables and probability mass functions in the discrete variables, see, for example, [40].The notation P(x) is reserved for probability mass functions (p.m.f.).

Factor graphs and sum-product algorithm
Let p(x 1 , x 2 , . . ., x n ) represent a joint p.d.f. of n continuous random vectors, which can be factorized as where {ψ 1 , ψ 2 , . . ., ψ N f } is the set of factors, X i ⊂ X is a subset of variables that are argu- ments of the factor ψ i , and X {x 1 , x 2 , . . ., x n }.
A factor graph [36,37] provides a visual representation of the factorization in (1), streamlining the computation of the marginal p.d.f.'s p(x j ) , where j ∈ {1, 2, . . ., n} .This computa- tional efficiency is achieved through message passing operations over the graph, extracting insights from the joint p.d.f.p(x 1 , x 2 , . . ., x n ) .A factor graph is composed of variable nodes representing arguments of the function on the left-hand side of (1) and factor nodes corresponding to the real-valued functions ψ i on the right-hand side.The interconnection between variable node x j and factor node ψ i occurs if and only if x j ∈ X i .As an example, Fig. 1 presents a simple factor graph with the factorization of the joint p.d.f.where factor and variable nodes are represented by squares and circles, respectively.The solid arrows are related to the message passing strategy presented next.
Sum-Product Algorithm The sum-product algorithm (SPA) is a message passing scheme that, in its broad configuration, calculates two messages for every edge in the graph, with one message propagated in each direction [37].
Let η ψ i →x j and η x j →ψ k denote a message from a factor node ψ i to a variable node x j and a message from a variable node x j to a factor node ψ k , respectively.Let also X i denote the set of all variable nodes that are connected to the factor node ψ i and, similarly, let j denote the set of all factor nodes that are connected to variable node x j .The SPA computes the mes- sages η ψ i →x j and η x j →ψ k over the respective edges of the graph as follows: where A\B {a|a ∈ A ∧ a / ∈ B} denotes the set difference operation between sets A and B.
(1) In particular, the message from a factor with no incoming connections, i.e., X i = ∅ , is the corresponding factor function ψ i (•) , and a message from a variable node x j that has no incoming messages, i.e., j = ∅ , is just a constant.Finally, the outgoing message of a variable node of degree 2 is just the copy of the incoming message.
Specifically, the message from a factor with no incoming connections (i.e., X i = ∅ ) consists of the associated factor function ψ i (•) .Similarly, a message from a variable node x j with no incoming messages (i.e., j = ∅ ) is just a constant Lastly, for a variable node with a degree of 2, the outgoing message is simply a duplication of the incoming message.
Following this framework, the exact expression for the marginal density function of x j in a tree-like graph is given by For a cyclic graph, the SPA is executed iteratively [37].The expression in (5) serves as an approximation for the marginal p.d.f. of interest.It is important to note that the SPA messages in (3), (4), and (5) must be normalized, such that η ψ i →x j (x j )dx j = 1 , η x j →ψ k (x j )dx j = 1 , and p(x j )dx j = 1.

Problem formulation
Consider a collection of N a aircraft, each equipped with onboard inertial naviga- tion system (INS), altimeter and GNSS receivers, which collaborate through a message passing algorithm over a partially connected network to improve their local position estimates.We model the network topology as an undirected graph, G = (V, E) , with vertices V {1, 2, . . ., N a } representing the aircraft and unconstrained edges E {(r, s) ∈ V × V ∧ s ∈ R(r)} denoting the communication links between nodes r and s, where R(r) represents the set of nodes directly connected to node r.We also model the ionospheric scintillation events at each GNSS receiver following a multimodal formulation, where we can have different observation models according to the absence or presence of scintillation.Following this approach, we extend the state vector with a discrete switching vector, which represents scintillation occurrences at the receiver.For this, denote as T k,r ] T the state vector containing the unknown states detailed in the sequel.

Position states/GNSS biases
Let x k,r [p T k,r b k,r ḃk,r ] T represent the unobserved continuous state vector composed of the tridimensional (3D) position state vector of aircraft r at the discrete-time instant k, the GNSS receiver clock bias and its derivative, respectively.The unobserved state vector is given at instant k + 1 by (5) p(x j ) ∝ m∈� j η ψ m →x j (x j ).(6) x k+1,r = F k,r x k,r + u k,r + w k,r where I n and 0 n denote, respectively, an identity matrix and a zero matrix of dimension n × n , u k,r [ x k,r y k,r z k,r 0 0] T represents an augmented input vector containing the dis- placement-estimated by an INS-within the 3D space between instants k and k + 1 , T s is the sampling period, and {w k,r } , k ≥ 0 , is an independent and identically distrib- uted (i.i.d.) sequence of Gaussian random vectors with zero mean and covariance matrix Q r = blkdiag{Q 1 , Q 2 } , and with σ 2 f , σ 2 b , σ 2 d representing, respectively, the 3D position, bias and bias drift noise variances.

Scintillation switching variable
Let k,r denote the unobserved scintillation switching vector of node r at instant k.The switching vector is modeled in terms of a discrete Markov chain with finite set of states M = {0, 1} N s , where N s denotes the number of in-view satellites at node r.The finite set M consists then of 2 N s = |M| elements denoted by {M i } i=1,...,|M| .The state transition matrix is composed of the probabilities of transition between different modes in M , e.g., p ij = P({� k,r = M i |� k−1,r = M j }), ∀i, j ∈ {1, . . ., |M|} , which vary according to the ionosphere dynamics.In this paper, however, we assume, without loss of generality, that the transition probabilities p ij are set a priori and fixed over the simulation window.

Pseudorange Model
Assuming that at instant k, the receiver r is connected to N s = |S(r)| satellites, where S(r) is the set of in-view satellites for receiver r, we write the pseudorange between r and s ∈ S(r)1 as where s k,s is the satellite s ephemeris for instant k received in the navigation message, || • || denotes the Euclidean norm in R 3 , m I k,s is the ionospheric delay error, σ 2 k,s stands for the pseudorange variance, which is a function of the discrete-valued mode switching variable k,r←s , and {ε k,r←s } , k ≥ 1 , is an independent, identically distributed (i.i.d.) sequence of zero-mean Gaussian variables with unit variance, i.e., ε k,r←s ∼ N (0, 1).
We assume for simplicity that the ionospheric delay m I k,s is known and can be subtracted from the receiver measurements since corrections to that effect are often received in the navigation messages or estimated off-line using current existing models, such the Autoregressive-Moving-Average (ARMA) model [41,42] and Global Ionospheric Maps (GIM/ CODE) [43], and broadcast by ground stations to the aircraft.However, ionospheric scintillations due to randomly occurring Equatorial Plasma Bubbles (EPBs), particularly in lowlatitude regions, lead to rapid fluctuations in the amplitude and phase of GPS signals [44], causing additional random positioning errors.Any positioning system must account then for both ionosphere delay estimation and detection of scintillation.
In this paper, we follow the approach in [45] and model the effect of ionosphere scintillation as random changes in the variance of the pseudorange measurement error.Specifically, the discrete-valued switching random variable k,r←s is defined such that k,r←s = 0 means that no scintillation occurs at instant k in the channel between aircraft r and satellite s.Conversely, k,r←s = 1 denotes presence of scintillation at instant k in the same channel.
In the nominal condition of absence of scintillation, σ 2 k,s ( k,r←s = 0) becomes the stand- ard GNSS positioning error variance, whereas, in scintillation mode, σ 2 k,s ( k,r←s = 1) becomes the scintillation error variance.The values of σ 2 k,s ( k,r←s ) were chosen according to an α-µ distribution, which is a simple and flexible fading model, introduced in [46], that describes the distribution of signal amplitudes, during scintillation events, in terms of fading parameters α and µ .The p.d.f. of the normalized amplitude envelope of the received signal, assuming that the resulting average signal power is equal to 1, is given by where Ŵ(•) is the gamma function.The use of this particular fading model presents better results when compared to other distributions, see [46][47][48][49] for more details.Following [48], we define σ 2 k,s ( k,r←s ) as where B n is the phase-locked loop (PLL) single-sided noise equivalent bandwidth, c/n 0 is the nominal carrier to noise density ratio for the coarse/acquisition (C/A) L1 carrier and η is the pre-detection integration period.

Peer-distance Model
We model the measured distance between nodes r and s ∈ R(r) at instant k as where {ζ k,r←s } , k ≥ 1 , is an i.i.d.sequence of zero-mean Gaussian random variables with variance σ 2 z . (9)

Problem Statement
We formulate the problem as a recursive estimation of the augmented state vector X k,r , given all available measurements Y k and Z k , where Y k is the collection of all available GNSS pseudorange measurements over the network in the set {y 1:k,r←s |r ∈ V ∧ s ∈ S(r)} and Z k is the collection of all peer-distance measurements over the network in the set {z 1:k,r←s |(r, s) ∈ V × V ∧ s ∈ R(r)} .This recursive estimation can be achieved with the marginalization of the mixed posterior density function, p(X k |Y k , Z k ) .However, as direct marginalization is computationally intensive, we make use of factor graphs and the sum-product algorithm (SPA) to obtain an approximation of p(X k,r |Y k , Z k ).

Message passing scheme
In this Section, we present the message passing scheme that approximates the marginal posteriors p(X k,r |Y k , Z k ) of each node r in the network.To accomplish that, we restrict the propagation of messages in time to be forward-only.Such assumption is in accordance with real-world applications, where messages flow only in one direction.This also reduces internode communication cost and CPU load [34].
From the conditional independence assumptions in the signal model, it follows that p(X k |Y k , Z k ) is factored as Figure 2 shows the factor graph representation of the factorization in (13) for a fixed time instant k, and the corresponding message passing scheme over the graph.In the illustration, we drop the variable nodes and message arguments for simplicity, thus any (13) r s Fig. 2 Factor graph and message passing scheme for the proposed cooperative localization scheme, where f represents the dynamic model, g the GNSS likelihood model and h communication likelihood.Also, red arrows refer to local messages, whereas blue arrows represent messages sent over the network function f k,r (X k,r ) becomes f k,r and any message η f k,r →X k,r (X k,r ) becomes η f k,r →X k,r .In addition, in the following derivation, we drop for simplicity of notation the conditional dependence on the observations, such that the posterior p(X k,r |Y k , Z k ) is denoted sim- ply as p k|k,r (X k,r ).Consider node r and time step k, using the dynamic model in (6), the factor f k,r is given by where with ′ k,r representing the value assumed by the scintillation switching vector at instant k and node r, and * k−1,r is the value assumed by the switching vector at instant k − 1 and node r.Also, where m ∈ R N x , represents the mean vector and S represents the covariance matrix.
Let Y k,r = {y k,r←s | s ∈ S(r)} be the collection of all pseudorange measurements available at node r at instant k.Factor g k,r is associated with the GNSS observation model in (8) and is given by Note that the satellite ephemeris s k,s does not appear as a variable node on the factor graph since it is received in the GNSS data and is, therefore, known at each node r.Likewise, as discussed in Sec.3.2, we assume that the ionospheric delay m I k,s is estimated off-line and broadcast to the GNSS receiver, so this term is also treated in our model as a deterministic and known parameter rather than a random state variable.
Let now Xk,r be a vector that collects both X k,r and {X k,s | s ∈ R(r)} , i.e., all state vectors in the closed neighborhood of node r at instant k.In turn, let Z k,r = {z k,r←s | s ∈ R(r)} be the set of all (unconstrained) internode distance meas- urements available at node r at instant k.Based on (12), the cooperative factor h k,r is given by ( 14)

Messages
We now introduce the messages outlined in Fig. 2. In the proposed factor graph, there are five different messages: temporal messages, GNSS messages, messages to neighbors, messages from neighbors and beliefs.The messages can be divided into local messages, where no network exchange is required, and cooperative messages.In this context, temporal and GNSS messages can be seen as local messages, whereas messages to/from neighbors are cooperative messages, which require internode communication.Beliefs are hybrid messages since, depending upon the cooperation status, they could be both local or cooperative messages.

Local Messages
The message η f k,r →X k,r is computed based on factor f k,r in ( 14), and the previous belief ) .Assuming that the previous belief has the form we write the message η f k,r →X k,r as ( 16)

GNSS messages
The GNSS messages η g k,r →X k,r simply propagate the GNSS factor g k,r in (15), since this factor has no incoming connections, i.e.,

Cooperative messages
As the factor graph in Fig. 2 presents a cyclic nature, cooperation can only be achieved in a loop with N L iterations.Therefore, in the following derivations, we include the superscript l in the messages to/from neighbors, denoting the iteration inside the cooperative loop.
Messages to neighbors At each iteration step of the cooperative loop, each node broadcasts the message η (l) X k,r →h k,s , representing the last calculated belief of X k,r .Thus, each node s in the neighborhood R(r) receives Messages from neighbors Messages coming from neighboring nodes at instant k use the likelihood of all peer-distance measurements z k,r←s , ∀s ∈ R(r) , received by node r, given the states X k,r and X k,s .This message is computed as follows: (19) Beliefs The belief at node r at instant k is an approximation of the marginal posterior p k|k,r (X k,r ) , based on all information available at node r.Following the SPA scheme, beliefs are computed as Note that, when no cooperation is performed we have N L = 0 , thus the last term of ( 22) is omitted, that is, node r is flying using GNSS-INSS fusion only.After N L iterations the belief p (N L ) k|k (X k,r ) is propagated to the next time step and the whole process is repeated for step k + 1 .Finally, it is noteworthy that when node r is flying in a GNSS-denied sce- nario, the GNSS messages are not calculated.

Gaussian-SMC implementation
In order to reduce computational and communication load, we compute the aforementioned messages following a parametric approximation.Furthermore, to assimilate the incoming nonlinear measurements at each step k, we employ a Monte Carlo representation of the beliefs and the message η (l) h k,r →X k,r (X k,r ) at each iteration loop l.

Temporal Messages
We assume that the belief p k−1|k−1,r (X k−1,r ) at instant k − 1 and at node r is repre- sented by a weighted sample set {(w (i)  k−1,r , X (i) k−1,r )} , i ∈ {1, . . ., N p } .Using a linear mini- mum mean square error (LMMSE) approach, we approximate p k−1|k−1,r (x k−1,r |� (i)  k−1,r ) by a multivariate Gaussian function with parameters defined as where the posterior means xk−1,r and ¯ k−1,r , and the posterior covariance matrix (22) are estimated using the sample means and sample covariance or cross-covariance matrices associated with the weighted sample set {(w Following this approach, we can approximate the temporal message η f k,r →X k,r in (18) as where We then sample new particles X k,r ∼ ηf k,r →X k,r (X k,r ) , i ∈ {1, . . ., N p } and update the particle weights as w (i,0) k,r ∝ η g k,r →X k,r (X (i,0) k,r ) , followed by a normalization, such that The weight update in this case corresponds to the assimilation of local GNSS pseudorange measurements.

Neighboring Messages
Cooperation is achieved by sharing the belief p k|k,r (X k,r ) with neighboring nodes s ∈ R(r) .However, since the belief follows a Monte Carlo representation this process results in communication burden.To overcome this, we build a Gaussian approximation of the last marginal belief p (l−1) k|k,r (x k,r ) of the continuous state vector x k,r at each coopera- tive iteration, using the Monte Carlo representation {(w )} .The parameters of this approximation are then sent to neighboring nodes of each node r.Upon receiving the Gaussian parameters from node s ∈ R(r) , node r resamples N p particles from that received parametric approximation.Therefore, at instant k and cooperative iteration l, each node r has an approximate Monte Carlo representation {(w 21) can then be approximated as Each node r then builds a new updated Monte Carlo representation {(w ) , at iteration l, by propagating the particles and updating the weights as After N L iterations of the cooperative loop, the approximated marginal beliefs p k|k,r (x k,r ) and P k|k,r (� k,r ) at each node r are computed based on the particle set {(w (i,N L ) k,r , X (i,N L ) k,r )} as follows: where δ in ( 28) is the Dirac delta and I in (29) is the indicator function.The weighted particle set is then propagated to the next time step.(25) η(l)

Complexity analysis
The complexity of the algorithm is dominated by the approximated message from neighbors in (25), which scales to O(CN L N 2 p ) , where C = |R(r)| is the number of neighboring nodes of node r.Also, for node r at instant k, the complexity of Algorithm 1 scales as O(N p (1 + N s ) + CN L N 2 p ) .Note that, in case of no cooperation, the complexity reduces to O(N p (1 + N s )).

Numerical simulations
In this Section, we first introduce the simulation setup and the performance metrics, and then, we present and discuss the numerical simulation results.

Simulation setup
We simulated N a = 9 aircraft flying on the surrounding area of a Brazilian airport for a period of 500 time steps, using flight trajectories obtained in [50].In the simulations, we considered a scenario where all aircraft have access to GNSS observations of N s = 6 satellites in view and cooperate in order to improve their own local estimates.The static network topology is depicted in Fig. 3 for time instants (150,175,200,225) seconds.
We then simulated our cooperative algorithm considering the Markov transition matrices depicted in Fig. 4.Each matrix, depicted in the panels of Fig. 4, was chosen in a manner that favors at least one channel in scintillation mode.At each simulation, we select for simplicity a subset of satellites in the full set of in-view satellites for the nodes, and assume that only the selected ones are subject to scintillation events.In addition, since EPBs are usually spread for hundreds of kilometers [51,52], the chosen subset of satellites is the same for the N a = 9 aircraft.In the transition matrix in Panel (a), we assume that only one GNSS channel can be faded by ionospheric scintillation events, whereas, in Panels (b) and (c), we assume, respectively, two and three affected channels.For the channels that are assumed to be unaffected by possible scintillation, we make the switching variable equal to zero at all time instants.The ground truth for the position/GNSS biases of each aircraft was generated using the process model in (6) with σ f = 3 m and GNSS clock noise parameters as in [53].The displacement of the aircraft from instant k to instant k + 1 was generated using the flight trajectories, in such way that the process noise w k,r also accounts for INS uncertainties as in [54].In addition, we make σ z = 5 m and set, for simplicity, the pseudorange stand- ard deviation σ k,s ( k,r←s ) to 4 m and 20 m for the cases in the absence and presence of scintillation, respectively.Also, the ground truth for the scintillation switching vector of each node was generated according to the steady-state distribution of the Markov chain, computed according to the transition probability matrix.
In order to properly initialize the filter, the prior belief of x 0,r was assumed to be a Gaussian p.d.f.centered at x0,r with covariance matrix P0,r = diag{25 2 , 25 2 , 25 2 , 1 2 , 0.1 2 } .The initial state estimate x0,r was then sampled from a Gaussian distribution centered at the true state with the same aforementioned covariance matrix.In addition, the prior belief of 0,r was initialized with equal probability for the different modes.
For the experiments, we fixed the maximum number of iterations in the cooperative loop to N L = 1 , since we did not empirically verify a sizeable performance gain when the number of iterations was increased.We performed M = 300 independent Monte Carlo runs and used N p = 500 particles.

Performance metrics
In order to measure the performance of the proposed algorithm, we define the following metrics.Let p (m)  k,r and p(m) k,r denote the ground truth and the estimated 3D position of aircraft r at instant k in the m-th Monte Carlo run, respectively.The root mean square error (RMSE) of the position estimate at node r at instant k is defined as where M is the total number of Monte Carlo runs.
We also compute the network RMSE considering the position estimates at all network nodes r ∈ {1, . . ., N a } at instant k as Finally, we define the overall RMSE considering the network position estimates at all time instants k ≥ 0 as where N is the total number of steps in the simulation.(30) (

Simulation results
We started our analysis by setting the number of channels subject to fading events to 2. We then simulated both non-cooperative and cooperative versions of the proposed localization algorithm, labeled 'GMarkov-SMC non-coop' and 'GMarkov-SMC coop, ' respectively, and compared them to two standard extended Kalman filters (EKF)-an optimistic version, referred as 'EKF-opt, ' which assumes complete absence of scintillation and a pessimistic version, 'EKF-pes, ' which assumes the presence of scintillation at every instant k and receiver r.For both EKF simulations, each aircraft performs localization with no cooperation.Next, we computed the network position RMSE using the metrics defined in Sec.6.2.The 3D positioning errors for each tested algorithm are shown in Fig. 5.As noted in the Figure, the aircraft can navigate safely, even without cooperation, as the network position RMSEs do not diverge.However, both EKF algorithms ('EKF-opt' and 'EKF-pes') present higher position estimation error compared to our proposed algorithms ('GMarkov-SMC non-coop' and 'GMarkov-SMC coop').
The overall position RMSE for the different algorithms is shown in Fig. 6.In this Figure , we see that the optimistic EKF yielded an overall position error of 10.33 meters, while the pessimistic version overall error was 13.28 m.The non-cooperative version of the proposed GMarkov-SMC method ended up with an average RMS error of 6.13 meters, thus resulting in more than 4 meters of difference when compared to the optimistic EKF.Considering the scenario where the aircraft cooperate with neighboring devices, we end up with an overall error of 4.89 m.
We also computed the cumulative probability distribution functions (c.d.f.s) of the global vertical position errors over all nodes in the network, shown in Fig. 7.One can readily note that for the cooperative version of the proposed filter ('GMarkov-SMC coop'), the vertical error remains below 4 m (approximately 3.3 m) for 95% of time, which meets some standards of navigation for aircraft in the approach and landing phases [55].
Next, we compared the proposed algorithm in this paper with the Gaussian-SMC algorithm as in [34,56], denoted here 'GSMC, ' whose equations do not handle scintillation effects.In the comparison, we simulated both non-cooperative and cooperative versions of the GMarkov-SMC and GSMC algorithms, and simulated different configurations regarding the number of channels that can be affected by scintillation. Figure 8 illustrates the overall position RMSE for each simulated algorithm.It can Fig. 6 Overall 3D position estimation error of the entire network at all time instants and all Monte Carlo simulation-assuming that only one channel can be affected by ionospheric scintillation events Fig. 7 CDF of vertical position error for cooperative, non-cooperative and EKF simulations be noted in Fig. 8 that our proposed GMarkov-SMC algorithm delivers the smallest positioning errors, both in the non-cooperative and cooperative versions, in all three scenarios where 1, 2, 3 and 4 channels can be affected by scintillation.Also, in the scenario where 3 channels are affected, GMarkov-SMC coop (hard blue bars) achieves an overall position error reduction of more than 1.5 meters when compared to GMarkov-SMC non-coop, and more than 5 meters when compared to GSMC  non-coop.For the scenario where 4 channels are affected, the 'GSMC non-coop' simulation delivered an overall position RMSE above 13 meters, while the proposed algorithm delivers 9.21 m.Also, when enabling cooperation for this scenario, the proposed approach reduces the overall position RMSE by 3 meters.Although the occurrence of scintillation events in several channels is not a likely event [17], the filter was able to deliver small stable RMSE values for these scenarios.
To test filter robustness, we then simulated a mismatched version of the algorithm.In this simulation setup, while the ground truth for the pseudorange measurements was generated using σ k,s ( k,r←s ) = 20 m, the filter equations considered σ k,s ( k,r←s ) = 12 m. Figure 9 presents the simulation results for the aforementioned setup.In the illustration, it can be readily seen in Fig. 9 that the filter can handle mismatched values of pseudorange standard deviation.
Figure 10 represents a simulation following the same setup described in Sect.6.1, but with no mismatch between the filter and the ground models, σ k,s ( k,r←s = 0) = 4 m and σ k,s ( k,r←s = 1) = 12 m.In Fig. 10, we noticed the same previous observed pat- tern in the filter responses.However, as the gap between the pseudorange standard deviations used in the simulations was diminished, the 'GSMC coop' algorithm delivered an overall position RMSE which is lower than that of the 'GMarkov-SMC noncoop' scheme.Once again, the 'GMarkov-SMC coop' algorithm delivered the lowest overall RMSE, despite the decreased gap between pseudorange standard deviations, respectively, in the absence and presence of scintillation.
Communication Cost At each iteration loop, each node r needs to broadcast its parametric representation of p (l−1) k|k,r (x k,r ) , where x k,r is a 5 × 1 vector that collects the 3D position p k,r of the aircraft r, the receiver bias b k,r , and the bias drift ḃk,r .As we use a Gaussian parametric approximation for internode communication, the broadcast messages over each internode communication link consist of 20 real values, i.e., 5 for the mean vector and 15 for the symmetric covariance matrix.For the network topology, depicted in Fig. 3, there are 9 nodes and 60 internode links, thus resulting in an average of 133.34 real numbers exchanged per node at each iteration loop l.Nodes with maximum cardinality (in this case, equal to 8) transmit 160 real numbers per loop, whereas nodes with minimum cardinality (equal to 4 in the assumed topology) send 80 real numbers per loop.As remarked before, only one loop iteration was sufficient in our simulations to achieve accurate positioning.

Conclusions
We presented in this paper a cooperative localization method following the factor graph framework and the SPA algorithm.We simulated a set of aircraft flying on the surrounding area of a Brazilian airport, exchanging messages with neighboring nodes in order to improve the self-localization problem and decrease the state estimation error.
We developed the factor graph formulation for the problem as a distributed message passing algorithm, using a Gaussian parametric approximation that reduces the internode communication load.Since each aircraft may be subject to scintillation events, commonly present in low-latitude regions, we also modeled the occurrence of these events with a discrete-time Markov process at each node.We then evaluated both noncooperative and cooperative versions of the algorithm according to different simulation setups.For each simulation, we evaluated the position estimation performance for both versions based on pre-defined RMSE metrics.
Simulation results showed improvements in the 3D position estimation error when compared to standard EKFs, where each node performs self-localization without communication, and Gaussian-SMC approaches that ignore scintillation events.In addition, cooperation significantly reduced the overall position RMSE even in scenarios where multiple channels can be affected by scintillation events.
In summary, our results suggest that the proposed cooperative localization algorithm is potentially a suitable candidate for air navigation in scenarios where the aircraft may be subject to ionospheric scintillation events.

2 Fig. 1
Fig. 1 Simple factor graph.Factor and variable nodes are represented by squares and circles, respectively.Bold lines represent the dependency of factors and variables

Fig. 3 Fig. 4
Fig. 3 Static network topology used in our simulation with aircraft (red circles) and communication links (blue lines)

Fig. 5
Fig. 5 Network 3D position estimation error for the EKF and non-cooperative and cooperative version of the proposed Gaussian-Markov filter

Fig. 9
Fig.9 Overall position estimation error of the non-cooperative and cooperative versions of the proposed filter and a Gaussian-SMC algorithm, for different configurations where several channels might be affected by scintillation, considering a mismatch of pseudorange standard deviation

Fig. 10
Fig. 10 Overall position estimation error of the non-cooperative and cooperative versions of the proposed filter and a Gaussian-SMC algorithm, for different configurations where several channels might be affected by scintillation, considering σ k,s ( k,r←s = 0) = 4 m and σ k,s ( k,r←s = 1) = 12 m Overall position estimation error of the non-cooperative and cooperative versions of the proposed filter and a Gaussian-SMC algorithm, for different configurations where 1, 2, 3 and 4 channels can be affected by scintillation