Robust all-source positioning of UAVs based on belief propagation

For unmanned air vehicles (UAVs) to survive hostile operational environments, it is always preferable to utilize all wireless positioning sources available to fuse a robust position. While belief propagation is a well-established method for all source data fusion, it is not an easy job to handle all the mathematics therein. In this work, a comprehensive mathematical framework for belief propagation-based all-source positioning of UAVs is developed, taking wireless sources including Global Navigation Satellite Systems (GNSS) space vehicles, peer UAVs, ground control stations, and signal of opportunities. Based on the mathematical framework, a positioning algorithm named Belief propagation-based Opportunistic Positioning of UAVs (BOPU) is proposed, with an unscented particle filter for Bayesian approximation. The robustness of the proposed BOPU is evaluated by a fictitious scenario that a group of formation flying UAVs encounter GNSS countermeasures en route. Four different configurations of measurements availability are simulated. The results show that the performance of BOPU varies only slightly with different measurements availability.


Introduction
Unmanned air vehicles (UAVs) require an accurate estimate of their positions, velocities, and attitudes in order to control themselves, navigate, and reason about their environment. The way this is achieved varies greatly from systems to systems. While most current UAV navigation systems rely on a combination of the Global Navigation Satellite Systems (GNSS) and an inertial measurement unit (IMU), there is a trend towards the use of all navigation sources available to meet the endless pursuit of navigation robustness in the face of new threats and mission challenges [1][2][3].
For UAVs, besides GNSS signals, ranging with ground control stations and peer UAVs is readily achievable. Recently, signals of opportunities (SoOPs) from existing RF background infrastructure, such as digital terrestrial wireless TV signal, have also aroused much interests in the research community [4][5][6]. Under such a context, an innovative positioning algorithm is needed to fuse a right position utilizing all these measurements.
Existing positioning algorithms in the literature are enormous. Early pseudorange-only positioning algorithms for *Correspondence: chenxiee@tsinghua.edu.cn Space Center, Tsinghua University, Haidian District, Beijing, 100084, China GNSS were based on iterative least square and Kalman filtering [7,8]. Extended Kalman filtering is a natural extension to Kalman filtering for solving nonlinear problems using one-order linearization [9]. A big step forward for Kalman-based filtering is the invention of unscented Kalman filtering [10][11][12]. In cooperative positioning [13][14][15][16][17][18][19][20], nodes have not only pseudoranges from navigation satellites but also ranging information with wireless peers. Existing algorithms such as iterative least square and Kalman filters can be extended to cooperative positioning, which leads to cooperative least square and cooperative Kalman filtering algorithm [21]. In recent years, convex optimizations, including semidefinite programming (SDP) [22][23][24][25], sum of squares (SOS) relaxation [26,27], and distributed gradient algorithm, also found their place in cooperative positioning. Another category of positioning algorithms are Bayesian approaches with belief propagation as an outstanding representative. Belief propagation-based cooperative positioning was studied for wireless sensor networks, mobile ad hoc networks, vehicle ad hoc networks, etc. [28][29][30][31]. By jointly using ranges with peer nodes and pseudoranges from satellites, cooperative positioning dramatically improves the availability and accuracy of positioning, especially in GNSS-challenged environments. http://asp.eurasipjournals.com/content/2013/ 1/150 Despite that positioning has been treated in various settings, there lack a study for positioning of UAVs with all wireless sources available, which is of fatal importance for UAVs to survive hostile operational environments. In fact, the US Defense Advanced Research Projects Agency (DARPA) had already initiated its All Source Positioning and Navigation (ASPN) program as emerging threat scenarios become more sophisticated and widespread. Under all-source positioning context, flexibility and robustness, which allow for run-time join and leave of measurements, are of the first concern. Accuracy, however, is of the second concern. In existing work, algorithms based on belief propagations have proved to be the best candidate to meet the above requirements. While belief propagation is a well-established algorithm, it is not an easy job to handle the mathematics therein when all wireless positioning sources are considered.
For the above motivations, all-source positioning of UAVs based on belief propagations is studied in this work. The positioning sources include (1) pseudoranges and carrier phases from GNSS satellites, (2) ranges and closed-loop Doppler from peer UAVs, (3) ranging information and closed-loop Doppler with ground control stations [32], and (4) time difference of arrival (TDoA) from the signal of opportunities of background wireless infrastructure. The contributions are as follows: (1) A unified mathematical framework for position and velocity estimation is developed, taking all the above measurements and their statistical characteristics. (2) Based on the mathematical framework, a positioning algorithm, named Belief propagation-based Opportunistic Positioning of UAVs (BOPU), is proposed. (3) The factor products, which are mathematically intractable, are solved by an unscented particle filtering. For the accuracy performance of belief propagation with particle filters over Kalman filters and cooperative least square algorithms has already been proved by existing work such as [28,29], we focused on evaluating the robustness of BOPU. Simulations are conducted with a fictitious scenario that a group of formation flying UAVs are under the support of ground control stations but encounter GNSS countermeasures en route. Four different configurations of measurements are simulated and compared. The results show that the performance of BOPU varies only slightly with measurements availability.
The rest of the paper is organized as follows: Section 2 formulates the problem, Section 3 gives the details of the proposed positioning algorithm, Section 4 presents simulation results and discussions, and Section 5 concludes the paper.

Problem formulation
Consider a group of formation flying UAVs that are carrying out a mission. All UAVs, together with their ground control stations, form a wireless network. The set of UAVs is defined by a wireless node set M of cardinality |M|. Without loss of generosity, it is assumed that only one GNSS constellation is available with a set of satellites S of cardinality |S|. There is also a set of SoOP signal sources G of cardinality |G|. The epoch sequence is denoted by t (0) , t (1) , . . ., t (k) . For a selected wireless node m ∈ M, denote by M (k) m the nodes that node m ranges with, by S (k) m the subset of satellites m is in view, and by G (g,k) m the nodes that m shares TDoA about SoOP signal source g ∈ G. The location of node m at epoch k is denoted by The clock bias of node m with reference to the selected GNSS constellation is denoted as b (k) m in meters, which notionally equals the product of the speed of light multiplying the clock bias of node m in seconds. Define the state of node m as the vector Node m performs the following measurements: s represents the sum of correlated errors which generally include tropospheric error, ionospheric error, and ephemeris error, while ε (k) ρs represents all uncorrelated errors following a Gaussian distribution. 2. Carrier phase elapsed φ (k) s→m during epoch k and k − 1 from satellite s ∈ S, which is in the form where |φ| m is the true value and ε (k) φs is the carrier phase observation Gaussian error. 3. Ranges r where n is the index of the neighbor and ε (k) r is the ranging error following a Gaussian distribution. A neighbor can be a peer UAV or a ground control station. The position of a UAV is unknown, but that of a ground control station is assumed to be known. Note r n↔m with neighbors n ∈ M (k) m . By closed-loop Doppler measurement [32], the clock differential of wireless neighbors is eliminated and the resulting f (k) n↔m contains the Doppler that is only related to the http://asp.eurasipjournals.com/content/2013/1/150 relative movement of the neighbor nodes involved, specifically, where λ is the wavelength of the carrier used in Doppler measurement, and From Equation 4, we have f n→m in meters with neighbors n ∈ G (g,k) m referring to g, which is in the form In (6), c is the speed of light and ε (k) d is the error in meters, following a Gaussian distribution. Attaining TDoA measurement requires synchronization of nodes n and m. Nowadays, one way to achieve this is to use high-quality atomic clock re-synchronized at every departure. In such case, ε (k) d increases slowly as clock drift accumulates with time.
For brevity, we further define the following notations: The goal of the positioning is to find the a posteriori distribution of x (k) m at each epoch k, given all the available observationsÕ where (1 : k) denotes the epochs from 1 to k. At epoch

The Bayesian inference
It is reasonable to assume that ranges with peer UAVs and the control stations are independent, and it is also often the case that the nodes in M move independently. The pseudoranges are independent when ignoring b (k) s .
The error induced by ignoring b (k) s will be discussed later. While the movement of a node can be measured readily by IMUs in many cases, it is not the case in all-source positioning because IMUs and wireless measurement are loosely coupled for flexibility. In this work, the movement of node m ∈ M is modeled as a second-order Markov process. With these assumptions, (7) can be rewritten as where M\m denotes all variables inX To determine the marginal (7) recursively at each epoch k, the integrand of (9) can be further decomposed as Now for each node m ∈ M, we define the following: , denoting the range measurement model of node m with node n at epoch k.
denoting the TDoA measurement model of node m to n with reference to SoOP source g at epoch k.
, denoting the peer-to-peer Doppler measurements of nodes m and n at epoch k.
With the above definitions, we have With Equation 17, we have the factor subgraph of node m as given in Figure 1. The factor subgraphs of all nodes m ∈ M, when interconnected, make up the complete factor graph. Figure 2 illustrates the complete factor graph of nodes of the simulation scenario given by Figure 3.

The sum product update rule
A belief propagation algorithm defines the sum product messages and their update rules over the factor graph. In our case, there are six classes of messages: Compute: 4: h Δ m →x m using Equations 19 and 20, ∀m ∈ M 5: h x m →x n using Equation 26, ∀m ∈ M 6: Broadcast h x m →x n to all neighbors ∀m ∈ M 7: for iteration i = 1 to I do 8: for nodes m ∈ M in parallel do 9: Receive h x n →Δ n,m from all neighbors. 10: Compute: 11: h Γ s,m →x m using Equations 21 and 22; 12: h Θ s,m →x m using Equation 24; 13: 14: Communicate h x m →x n to n, ∀n ∈ M (k) m .

The messages in parametric form
A compact parametric form of the messages involved in the proposed algorithm is needed to make it permissible to transmit over a wireless network with limited bandwidth, which are given below: 1. Dead-reckon message h Δ m →x m is associated to the state of node m from epoch k − 2 : k − 1 to k. From Figure 1, we have  Figure 3 The simulation scenario. A group of six UAVs are flying from a start point under the control of ground control station N 0 to a destination with ground control station N 1 , but experience a GPS countermeasure in the midway.
From its definition, the dead-reckon message is a Gaussian probability density function with mean , respectively. Among the many ways for dead reckoning, we set Following (19), it is trivial to derive 2. Satellite factor messages h Υ s,m →x m are associated to the pseudoranges of node m. The pseudorange of node m from satellite s is assumed to be bias free except the receiver clock bias; thus, we have where σ 2 s→m is the pseudorange error power of satellite s at node m andρ (k) s is the carrier phase smoothed pseudorange, which can be expressed in a recursive form as [33] We note here that h Υ s,m →x m is only contributable to l m and b m . 3. Satellite carrier phase factor messages h Θ s,m →x m are associated to satellite carrier phase measurements. For where 1 m→s is the unit vector directed from node m to satellite s, δf (k) s is the clock drift rate of satellite s. So we have We note here that h Θ s,m →x m is only contributable to v m . 4. Messages to peers h x m →x n are messages that node m sends to all neighbors, in the following form: It is hard to find the exact expression of Equations 25 and 26, so we approximate the result of message multiplication as a Gaussian distribution: where z is the normalization factor. In practice, the values of μ x (k) m→n and Σ x (k) m→n are approximated by an unscented particle filtering as Algorithm 2, which will be disposed of later. 5. Messages from range neighbors h Γ n,m →x m represent the contribution of ranging information from wireless neighbors. For the position of ground control stations are known, h Γ n,m →x m associated to the ranging information of node m to a control station n can be expressed as where σ 2 n→m is the ranging error power andr (k) n is the Doppler smoothed range, which can be expressed in a recursive form aŝ Ground control stations' Doppler messages h Ψ n,m →x m are associated to the Doppler information of node m to ground control stations. Similar to h Θ s,m →x m , we have Messages from other UAVs whose positions are not known can be expressed as To find out the parametric form of this distribution, we follow a divide-and-conquer approach. First, we can see that the mean of position l (k) m in this distribution is a ball with radius r (k) n↔m and center μ l (k) m→n , and its covariance is Σ x (k) n→m + σ 2 n→m I. Any valid point on the surface of the ball is restricted by Similarly, 6. Messages from SoOP neighbors h Ω g n,m →x m represent the contribution of TDoA measurements referencing SoOP source g, which can be expressed as The mean position l (k) m in this distribution is a ball with center l g and radius ||l g − μ n→m || − d (k) n→m , and its covariance is Σ x (k) n→m +σ 2 d I. We now define a vector Ensure: Updatedμ x ,Σ x after product 1: for i from 1 to N do 2: Update the particle with UKF: 3: -Calculate 2n a + 1 sigma points: 4: -Propagate the particle into future: 6: -Incorporate new observations: 12: Evaluate the importance weights of the new particle: 21: end for 22: Normalize the importance weights of all particles so that Σ i w (k) i = 1 23: Estimate new mean and covariance using weighted particles:

Calculating the factor products
The products of factors, i.e., Equations 25, 26, and 41, are mathematically intractable. This is a problem that permeates in many disciplines of sciences. The most widely known methods are importance sampling [31], bootstrap particle filter [34], and unscented particle filters [35,36]. While importance sampling is convenient and attractive, it suffers from the sample degeneracy problem. Bootstrap filter and unscented particle filter try to avoid this degeneracy by context-aware resampling, which eliminates the particles having low importance weights and proliferates particles having high importance weights. A bootstrap particle filter uses state update information for resampling, while an unscented particle filter further improves the bootstrap particle filter by estimating the first-and second-order moments of the new state incorporating new observations using unscented transform before resampling.
We present the variation of unscented particle filter for calculating the factor products of this work in Algorithm 2. In Algorithm 2, n a , η, W c0 j , and W c1 j are parameters related to unscented transform. n a = 7 is the number of states, which is 7 in our case. η = 3α 2 − n a , W c0 0 = η/(n a + η), and W c1

Setup
For belief propagation combined with varied linear and nonlinear filters is widely available in the literature, we focused on evaluating the robustness of BOPU with some discussions on the appropriateness of the approximations in BOPU. Simulations are conducted by MATLAB with a fictitious scenario as given in Figure 3. In Figure 3, a group of six UAVs indexed by 0 to 5 are flying in a formation from a start point under the control of ground station N 0 to a destination with ground station N 1 and a SoOP source G 0 in view, but experience a GPS countermeasure en route.
The positions of satellites in view are given in Table 1. All UAVs follow the same velocity but with a different point of departure as also given in Table 1. The velocity of the formation flying is given in Figure 4a, and the route of node 0 is given in Figure 4b. The simulated cases include the following: • Case 0 : Whenever a node has at least four satellites in view, it will not participate in the belief propagation but will offer the statistics of its own position. In addition, the measurements with the two ground control stations N 0 and N 1 and the SoOP source G 0 signal are not utilized. Case 0 reduces the positioning traffic over the wireless network to a minimum but is expected to have the poorest positioning performance. • Case 1 : All nodes take part in the belief propagation process as stated in Algorithm 1 but do not utilize the measurements with the two ground control stations N 0 and N 1 and the SoOP source G 0 signal. • Case 2 : All nodes take part in the belief propagation process as stated in Algorithm 1, and the measurements with the two ground control stations N 0 and N 1 and the SoOP source G 0 signal are used. • Case 3 : All nodes take part in the belief propagation process. Besides, the control stations also provide pseudorange differential corrections at each epoch and broadcast to other nodes. The differential corrections help effectively remove b In the simulations, ground control stations 0 and 1 are placed at ENU(0, 0, 10) and ENU(50,000, 50,000, 10), respectively. The raw pseudorange observation error power   G 0 is placed at midway ENU(25,000, 25,000, 10). Nodes exchange time of arrival with wireless neighbors; thus, TDoA is also available between neighbors and D Figure 4 gives the simulation results. As can be seen, from case 0 to case 3, the positioning performance improves almost steadily. In case 0, node 0 and node 1 have at least four satellites, and they determine their position using a traditional weighted least square positioning algorithm in order to save communication bandwidth. Without iterations with node 0 and node 1, the positioning performance of other nodes is being restrained. In case 1, all nodes take part in the belief propagation process, which is helpful in improving the positioning performance of all nodes (especially nodes 2 to 5), so the overall performance of case 1 is better than that of case 0. The outperformance of case 2 over case 1 comes from the full usage of all available observations, especially ranges and closed-loop Doppler with ground control stations, and TDoA observations from SoOP source G 0 . The measurements with the ground control stations N 0 , N 1 and SoOP source G 0 help improve geometric dilution of positioning in a big way. In case 3, ground control stations also generate corrections for pseudoranges, which directly improves the quality of pseudoranges, thus the positioning precision.

Results and discussions
The positioning performance of all nodes is given in Figure 4e,f in terms of root mean square error (RMSE). It follows the fact that the more observables, the better precision. In case 0, node 1 uses only the observation from four satellites in view, so its horizonal RMSE is even less than those of node 2 and node 3. Nodes 2 to 5, which have less than four satellites in view under some given GNSS interference, can still achieve positioning by utilizing the peer to peer measurements and the measurements with control stations. Node 5 experienced the strongest interference. Without any satellite in view but with a better geometric position, node 5 achieved even better positioning performance than nodes 3 and 4 that have satellites in view.
The results show that the performance of BOPU varies only slightly with different measurements availability. The main approximations made in the proposed BOPU are as follows: (1) The correlated errors of pseudoranges b (k) s are ignored in factorization. The simulation results of case 3 and case 2 show that such an approximation is acceptable. (2) The dead-reckon message (19) and (20) actually ignored nonzero off-diagonal values, and Equation 11 ignores the smoothing ofÕ

Algorithmic complexity
Given a node m at epoch k, we have one dead-reckon message, |S m | messages from SoOP neighbors. It uses N particles in UPF and I iterations in the product estimate. The core of UPF is UKF whose complexity is O(n 3 a ) [37], where n a is the number of states as stated before. The complexities of main computations in BOPU are listed in Table 2. As was shown in Table 2, the complexity of BOPU is dominated by message multiplications needed by messages to peers, which scales as O( m |)n 3 a ).

Conclusions
Nowadays, UAVs are playing an increasingly important role both in the military and in civil affairs. Worries on the robustness of GNSS have also been increasing with the maturity of GNSS countermeasures and proliferation of wireless devices. With the aim of providing a better navigator for UAVs, we investigated the positioning of UAVs with all wireless sources via belief propagation with unscented particle filtering. By jointly using the measurements from GNSS satellites, peer UAVs, ground control stations, and signal of opportunities, the proposed algorithm, which is named BOPU, provides an improved positioning robustness and algorithmically proven high precision. By being opportunistic, BOPU allows for not only flexible variations of measurements availability but also agile compromise between wireless bandwidth consumption and positioning performance when put into practice.