Collaborative emitter tracking using Rao-Blackwellized random exchange diffusion particle filtering

We introduce in this paper the fully distributed, random exchange diffusion particle filter (ReDif-PF) to track a moving emitter using multiple received signal strength (RSS) sensors. We consider scenarios with both known and unknown sensor model parameters. In the unknown parameter case, a Rao-Blackwellized (RB) version of the random exchange diffusion particle filter, referred to as the RB ReDif-PF, is introduced. In a simulated scenario with a partially connected network, the proposed ReDif-PF outperformed a PF tracker that assimilates local neighboring measurements only and also outperformed a linearized random exchange distributed extended Kalman filter (ReDif-EKF). Furthermore, the novel ReDif-PF matched the tracking error performance of alternative suboptimal distributed PFs based respectively on iterative Markov chain move steps and selective average gossiping with an inter-node communication cost that is roughly two orders of magnitude lower than the corresponding cost for the Markov chain and selective gossip filters. Compared to a broadcast-based filter which exactly mimics the optimal centralized tracker or its equivalent (exact) consensus-based implementations, ReDif-PF showed a degradation in steady-state error performance. However, compared to the optimal consensus-based trackers, ReDif-PF is better suited for real-time applications since it does not require iterative inter-node communication between measurement arrivals.


Introduction
In several engineering applications, e.g., target tracking or fault detection, multiple agents [1] that are physically dispersed over remote nodes on a network cooperate to execute a global task, e.g., estimating a hidden signal or parameter, without relying on a global data fusion center. Each network node is normally equipped with one or more sensors that generate local measurements and can process those measurements independently of the rest of the network. At the same time, however, the network nodes are also able to communicate with each other in order to build in a collaborative fashion a joint estimate of the hidden signals or parameters of interest that depends both on local and remote measurements. Ideally, that joint estimate should be equal to or, at least approximate the optimal global estimate that would be generated by a centralized processor with access to all network measurements.
Most of the previous literature in distributed signal processing on networks is based on linear estimation methods. Specifically, distributed versions of the Kalman filter were proposed e.g., in [2][3][4] to track unknown state vectors in linear, Gaussian state-space models. In situations, however, where the state dynamic model or the sensor observation models are nonlinear, the posterior distribution of the states conditioned on the network measurements becomes non-Gaussian (even with Gaussian sensor noise) and, therefore, the linear minimum mean square error (LMMSE) estimate of the states provided, e.g., by an extended Kalman filter (EKF) may differ from the true minimum mean square error (MMSE) estimate given by the expected value of the state vector conditioned on the measurements. In this paper in particular, we focus specifically on an application where multiple passive received-signal-strength (RSS) sensors jointly track http://asp.eurasipjournals.com/content/2014/1/19 a moving emitter assuming, at each network node, nonlinear observation models with possibly unknown static parameters.

Distributed particle filtering
In nonlinear scenarios, an alternative to approximate the true MMSE estimate is to use a sequential Monte Carlo method like particle filters [5,6]. Several distributed particle filters have been proposed recently, see a comprehensive review in [7], to handle nonlinear distributed estimation tasks. An important constraint in the design of a distributed estimation algorithm is, however, that most networks of practical interest are only partially connected, i.e., each node can only directly access neighboring nodes in its immediate vicinity according to the network topology. In particular, assuming conditional independence of the different sensor measurements given the state vector, a distributed particle filter (PF) normally requires the computation of a product of likelihood functions that depend on local data only [8]. To compute that product over the network in a fully distributed fashion and with local neighborhood inter-node communication only, previous references suggest using iterative average consensus [8], iterative Markov chain Monte Carlo move steps [9], or selective gossip algorithms [10]. Alternatively, we proposed in [11] to compute the likelihood product exactly in a finite number of iterations using either iterative minimum consensus [12] or flooding techniques [13]. However, both consensus or flooding-based solutions are very costly in terms of bandwidth requirements as they require multiple iterative inter-node communication between two consecutive sensor measurements. Previous works, e.g. [8,14,15], propose approximations aimed at reducing the communication cost, but, in all aforementioned schemes, processing and sensing at different time scales are still required.

Diffusion particle filtering
An alternative to circumvent the high communication cost of consensus algorithms is to use diffusion algorithms [16] which, contrary to the former, do not require multiple iterative inter-node communication between consecutive measurements. Diffusion algorithms are, however, suboptimal in the sense that they do not simulate at each time step the behavior of the optimal global estimator, but rather, at best, approximate the optimal global solution asymptotically over time.
In the distributed linear estimation literature, most diffusion schemes are based on convex combinations of Kalman filters, see e.g., [3]. Kar et al. proposed in [2] a different approach based on random information dissemination. In a previous conference paper [17], we introduced the random exchange diffusion particle filter (ReDif-PF), which generalizes and extends the methodology in [2] to a PF framework by basically using random information dissemination to build at each network node different Monte Carlo representations of the posterior distribution of the states conditioned on random sets of measurements coming from the entire network. Reference [17] assumed, however, that the parameters of the sensor observation model were perfectly known. In this paper, we extend the algorithm to a scenario with unknown parameters and derive in detail a Rao-Blackwellized [18] version of the ReDif-PF. In the specific application under consideration in the paper, the unknown parameters are the sensor variances, but most of the methodology in the derivation of the RB ReDif-PF is general and could be easily adapted to other signal models and applications provided that, in a fully Bayesian framework, the dynamic posterior probability distribution of the unknown parameters conditioned on the observations and on the simulated particles is a conjugate prior [19] for the likelihood function of the measurements.
An abbreviated description of the RB ReDif-PF may be found also in the short paper [20]. This paper consolidates and extends both [17] and [20] including detailed derivations and additional simulation results and comparisons. We also detail approximate versions of the RB ReDif-PF where we use Gaussian mixture models (GMM) [21] and moment-matching techniques inspired by [22] to reduce communication requirements.

Paper outline
The paper is divided into six sections and three appendices. Section 1 is the introduction. Section 2 describes the state and sensor models. Section 3 describes the centralized PF and also briefly reviews the equivalent broadcast, consensus, and flooding implementations introduced in [11]. Section 4 derives the ReDif-PF algorithm considering alternate scenarios with both known and unknown parameters. In the unknown parameter case, we derive in detail the Rao-Blackwellized version of the ReDif-PF and introduce approximate versions thereof that enable significant reductions in communication cost. The performance of the proposed algorithms is evaluated with simulated data in a realistic scenario with 25 sensors in Section 5. We compare the ReDif-PF algorithm in the unknown parameter scenario to the optimal centralized PF and its equivalent consensus implementations. In the known parameter case, we also compare the proposed ReDif-PF tracker to the Markov chain Monte Carlo distributed particle filter (MCDPF) in [9], to a linearized random exchange distributed EKF, which is a variation of the algorithm proposed in [2], and to a distributed bootstrap particle filter based on selective gossip as proposed in [23]. Finally, we present our conclusions in Section 6. http://asp.eurasipjournals.com/content/2014/1/19 Appendices 1 and 2 show the proof of some key results in the paper, and Appendix 3 describes the ReDif-EKF algorithm used for comparison purposes in Section 5.

Problem setup
For simplicity of notation, we use lowercase letters in this paper to denote both random variables/vectors and real-valued samples of random variables/vectors with the proper interpretation implicit in context.
Without loss of generality, we assume that the emitter trajectory is described by the white noise acceleration model [24] x n+1 = Fx n + u n (1) where x n x nẋn y nẏn T is the hidden state vector at time step n consisting of the positions and velocities of the target's centroid respectively in dimensions x and y; F is the state transition matrix; and {u n } is a sequence of independent, identically distributed (i.i.d.) zero-mean Gaussian vectors with covariance matrix Q. Matrices F and Q, parameterized by the sampling period T and the acceleration noise σ 2 accel , are detailed in [11,24].

Observation model
Let N (m, σ 2 ) denote the Gaussian probability distribution with mean m and variance σ 2 and denote by IG(a, b) the inverse-gamma probability distribution with parameters a and b. The measurements z r,0:n = z r,0 , . . . , z r,n in decibels relative to one milliwatt (dBm) at the rth node of a network of R RSS sensors are modeled as where v r,n ∼ N (0, 1), σ 2 r ∼ IG(α, β), ∀r ∈ R {1, . . . , R}, and {x 0 , {u n } , v r,n , σ 2 r are mutually independent for all n ≥ 0 and for all r ∈ R. The nonlinear function g r (·) in (2) is in turn given by [25] where x r represents the rth sensor position, ||.|| is the Euclidean norm, (P 0 , d 0 , ζ r ) are known model parameters (see [25] for details), and H is a 2 × 4 projection matrix such that H(1, 1) = H(2, 3) = 1 and H(i, j) = 0 otherwise. We also denote by N r the set of nodes in the neighborhood of node r. The real-valued constants {α, β} are the model's hyper-parameters. Note that in (2), we take a fully Bayesian approach and model the unknown sensor noise variances σ 2 r , r ∈ R, as random variables that are mutually independent for s = r and identically distributed a priori with an inverse-gamma distribution.

Problem statement and goals
Let z 1:R,0:n denote the set z r,t for all network nodes r = 1, . . . , R and all time instants t = 0, . . . , n. Given z 1:R,0:n , we want to compute the MMSE estimatê x n|n = E x n |z 1:R,0:n (4) at each instant n ≥ 0, where E{x n |z 1:R,0:n } denotes the conditional expectation of x n given z 1:R,0:n .
In the sequel, we first describe in Section 3 a recursive, centralized PF algorithm that approximates the desired global MMSE in (4) at each time step n in a scenario with unknown sensor variance scales σ 2 r . Next, we review in Section 3.1 two fully distributed algorithms that operate on a partially connected network and allow exact innetwork computation of the state estimate in (4) without a global data fusion center and with inter-node communication limited to a node's immediate neighborhood according to the network topology. The network connectivity is described by a graph G = (R, E) where R = {1, . . . , R} is the set of nodes and the graph has an edge (u, v) ∈ E, (u, v) ∈ R × R if and only if nodes u and v can communicate directly with each other. The particular network graph used in the simulation scenarios in this paper is described in detail in Section 5.
Finally, we introduce in Section 4 a novel diffusionbased algorithm, which is also fully distributed and relies on local inter-node communication only specified as before by the network graph G but, rather than yielding an identical estimate (4) at each node, obtains at each node r a suboptimal estimatê where Z r,0:n is a random subset of z 1:R,0:n , which is different at each node r and includes measurements coming from random locations in the entire network, as opposed to measurements coming only from node r and its neighborhood. Compared to the exact distributed implementations of the optimal global estimate in Section 3.1, the diffusion solution in Section 4, although suboptimal, is designed to have a much lower inter-node communication cost and, therefore, is better suited for real-time applications.

Centralized particle filter
In a centralized architecture, all nodes in the network transmit their local measurements to a data fusion center which then runs a particle filter that approximates the MMSE estimate of the unknown state vector at each time instant n as E x n |z 1: is a properly weighted Monte Carlo set [5,6] that represents the posterior probability density function (PDF) p(x n |z 1:R,0:n ) in the sense that the sum on the right-hand side of (6) converges, according to some statistical criterion, to the expectation on the lefthand side when Q → ∞. The random samples x (q) n , also called particles, are sequentially generated according to a proposal probability distribution specified by a socalled importance PDF π(x n |x (q) 0:n−1 , z 1:R,0:n ). If the blind importance function [5] π(x n |x where ∝ denotes 'proportional to, ' z 1:R,n is an alternative notation for the set z r,n , r ∈ R, and the proportionality constant on the right-hand side of (7) is chosen such that From (8) and (9), it can be shown then that (see the proof in [11]) p(z 1:R,n |x . (10) Substituting now (10) into (7), the centralized weight update rule reduces to

Equivalent distributed implementation of the centralized particle filter
Note that each factor λ (q) r,n (x (q) n ) in the product on the right-hand side of (11) depends only on local observations. In a fully connected network, assuming that all nodes r ∈ R start out at instant n − 1 with the same particles x (q) n−1 , they can all synchronously draw [26] new particles x (q) n according to p(x n |x (q) n−1 ), locally compute their own local likelihood functions λ (q) r,n (x (q) n ), and then broadcast them to the entire network until all nodes have all the remote likelihood functions and can compute the product on the right-hand side of (11). Synchronous multinomial resampling according to the global weights followed by regularization may follow (see [11]) to mitigate particle degeneracy and impoverishment [5,6]. The algorithm described in this paragraph is referred to as the decentralized particle filter (DcPF) in [11] and [27].
As mentioned, however, in Section 1, real-world networks are only partially connected and fully distributed computations of the product in (11) are needed. One possibility is to approximate the product using iterative average consensus [28] as proposed, e.g., in [8] and [29]. Alternatively, we introduced in [11] a fully distributed computation of the global weights in (11) using either iterative minimum consensus [12] or flooding [13]. Both algorithms assume only local communication between nodes in immediate neighborhoods and, to achieve an exact computation of the global weights, require only a finite number of iterative message exchanges between nodes in the time interval between two consecutive sensor measurements.
Let D denote the diameter of the network graph, i.e., the maximum number of hops between any two nodes and, as before, denote by R the number of nodes in the network. By running R × D consecutive minimum consensus iterations [12] for each particle q, it is possible (see details in [11]) to build an identical ordered list of likelihood functions λ (q) r,n (x (q) n ) , r ∈ R, at all nodes. Each node can then locally compute the product of the likelihoods as in (10) and obtain identical, optimal global importance weights w (q) n . We refer to that (communicationintensive) minimum-consensus-based distributed tracking algorithm as CbPFa.
A more efficient way, however, to compute the exact optimal global weights at each node is to flood [13] the local node likelihoods over the network. Flooding protocols allow one to (iteratively) broadcast values over a network relying on local neighborhood inter-node communication only. Given a partially connected sensor network, one can simultaneously flood the R distinct likelihoods over the network as follows. First of all, each node r maintains an ordered list of distinct likelihoods. A likelihood in turn is flagged to indicate that it has not been sent to node r neighbors yet. Initially, the node r stores its local flagged likelihood in its list. At a given iteration, node r sends its lowest flagged likelihood to all neighbors and then unflags it. Conversely, it receives remote likelihoods from nodes s ∈ N r . If a received remote likelihood http://asp.eurasipjournals.com/content/2014/1 /19 is not included in node r's list yet, it is inserted with a flag in its list. This procedure is guaranteed to converge in a finite number of iterations as soon as each node has R distinct values in its ordered list of likelihoods. We refer to the flooding-based iterative tracker in this paper as the CbPFb algorithm. Figure 1 illustrates how the proposed flooding protocol iteratively creates at each node r an ordered list comprising all likelihoods across the network in a toy example with three nodes where node 1 is connected to node 2, node 2 is connected to nodes 1 and 3, and node 3 is connected to node 2 only. A star symbol is employed to indicate which likelihoods are flagged in the ordered list maintained by each node r at a given iteration j.
Although optimal in the sense of reproducing the centralized solution, the minimum consensus and flooding algorithms in [11] are still communication-intensive due to the requirement of iterative inter-node communication between sensor measurement arrivals. In the next sections, we describe an alternative fully distributed diffusion-based solution that drops this requirement and is the main topic of this paper.

Random exchange diffusion particle filter
In this section, we derive an alternative distributed PF based on random information dissemination that extends the methodology in [2] to a Monte Carlo framework. We also present a Rao-Blackwellized version of the proposed distributed PF in a scenario with unknown sensor parameters.
Let Z s,0:n−1 denote the set of all network measurements assimilated by node s up to instant n − 1. Next, let x (q) s,0:n−1 with associated weights w (q) s,n−1 , q ∈ Q, be a properly weighted set that represents the posterior PDF p(x 0:n−1 |Z s,0:n−1 ) at node s. Assume now that at instant n − 1, node s sends its particles and weights to a neighboring node r that can assimilate at instant n the measurements Z r,n = z i,n , i ∈ {r} ∪ N r . At instant n, the new particle set at node r, x (q) is now a properly weighted set to represent the updated posterior p(x 0:n |Z r,n , Z s,0:n−1 ), where {Z r,n , Z s,0:n−1 } is redefined as Z r,0:n . Resampling from the particle weights followed by regularization may be added to combat particle degeneracy and restore particle diversity, i.e., for q ∈ Q (see also [11]): r,n ) , and h > 0 is an empirically adjusted parameter.
• Reset the particle weights w locations in the entire network, it suffices to implement a protocol where each node r, starting from instant zero, exchanges its particles and weights with a randomly chosen neighboring node s, propagates the received particles using the blind importance function as in (12), and then updates their weights as in (13). Figure 2 illustrates the evolution of the marginal posterior at each node -in a linear network containing three nodes running the random exchange protocol -over four time instants. Initially, each node r ∈ {1, 2, 3} has a posterior at instant zero conditioned on the measurements Z r,0 = {z i,0 }, i ∈ {r} ∪ N r , in its vicinity only. At each time instant n ∈ {1, 2, 3}, network nodes perform the sequence of random exchanges as indicated in the rightmost column of Figure 2 and, then, update the received posterior by assimilating measurements in their respective neighborhoods.
Note that in the linear network topology shown in Figure 2, node 2 always performs two random exchanges at each time instant n. Generally speaking, however, at a given instant n, a node r exchanges its parameters at least one time with a randomly chosen neighbor s and, in the worst case, performs d(r) random exchanges between two measurement arrivals with nodes in its vicinity, where d(r) is the degree of node r, i.e., the number of neighbor nodes.
Unlike randomized gossip algorithms [30], this procedure diffuses information by randomly propagating posterior statistics across the network. More specifically, as the initial posterior statistics provided by a given node r 0 at time 0 follows a path P {r 0 , r 1 , . . . , r n } along the network, it assimilates the available measurements Z r,n in the neighborhood of each visited node r ∈ P. Since, as illustrated in Figure 2, the initial posteriors at each node follow different paths, the posterior available at node r n at time n will be different from those in the remaining nodes. Thus, network nodes will provide different estimates conditioned on distinct sets of measurements.

ReDif-PF with known sensor variances
If the parameters of the sensor observation model at each node r are deterministic and perfectly known, then At instant n, then, upon receiving (w s,n−1 ), q ∈ Q, from node s, the particle filter at node r samples as before and updates its weights as where Inter-node transmission requirements From the previous discussion, it follows that in the scenario with known variances at each instant n, it suffices for each node s to transmit to the chosen neighbor r the set of particles {x (q) s,n−1 } (4Q real numbers for a four-dimensional state space) and the respective set of importance weights {w (q) s,n−1 } (Q real numbers). In addition, node s also sends its scalar observation z s,n and the known observation model parameters ( ζ s , x s , σ 2 s ) (see (3)) to all nodes i in the neighborhood of s.

Rao-Blackwellized ReDif-PF with unknown sensor variances
Let IG(σ 2 |α, β) denote the PDF of a continuous random variable σ 2 with an inverse-gamma distribution specified by the parameters α and β, i.e. [19], and zero otherwise. In (16), (.) denotes the gamma function Similarly, let also N(x|m, ) denote the PDF of a Gaussian random vector taking values in L and with mean m and positive definite covariance matrix , i.e., where | | denotes the determinant of the matrix and the superscript T denotes the transpose of a vector.
In the scenario with unknown sensor variances, it can be shown (see Appendix 2) that if at instant n − 1, then where i ∈ {r}∪N r , and each factorλ r,n ) in the product on the right-hand side of (18) is computed by solving the integral where (·), as before, denotes the gamma function with g i (·) calculated as in (3). Furthermore, at node r and instant n, the updated parameter posterior PDF where α r,i,n and β (q) r,i,n are updated as in (20) and (21) if i ∈ {r} ∪ N r or, otherwise, are kept equal respectively to α s,i,n−1 and β (q) s,i,n−1 . If regularization is used to combat particle degeneracy, the posterior parameters {β r,n }. We follow, however, a different suboptimal strategy described in Section 4.3, which also allows a significant reduction in inter-node communication cost.

Inter-node transmission requirements
In the unknown variance scenario, based on the previous discussion, at each instant n , a node s has to transmit to its (randomly chosen) neighboring node r its particle set x (q) s,n−1 (4Q real numbers) plus the respective importance weights w (q) s,n−1 (Q real numbers) and the set of hyper-parameters (α s,i,n−1 , β (q) s,i,n−1 ) , i ∈ R, q ∈ Q (another R×(Q+1) real numbers), which specify the posterior PDF p(σ 2 1:R |x (q) s,0:n−1 , Z s,0:n−1 ). In addition, as before, node s also sends its scalar observation z s,n and the observation model parameters ( ζ s , x s ) to all nodes i in the neighborhood of s.

Approximate RB ReDif-PF
Although the exact ReDif-PF algorithms in Sections 4.1 and 4.2 converge asymptotically to the state estimate in (5) as the number of particles Q goes to infinity, their inter-node communication cost is still relatively high. To reduce the communication burden, we propose two suboptimal approximations which are described in detail in the sequel. s,n−1 , q ∈ Q, at node s using the Expectation-Maximization (EM) [31] algorithm. Node s now transmits to node r only the parameters that specify the GMM model, i.e., 15K real numbers for a four-dimensional state vector, as opposed to 5Q real numbers, where typically Q >> K (in the simulations in Section 5 for example, K is either 1 or 2, whereas Q is 500). Node r then locally resamples Q new particles x (q) s,n−1 according to the received GMM PDF and resets its importance weights w (q) s,n−1 to 1/Q. Since resampling from the GMM approximation is used, we omit the regularization step mentioned in Section 4.

Approximation of the posterior distribution of the sensor variances
In the particular situation where the sensor variances are unknown, in theory we should also locally resample the previous particle trajectories x . . , R for the resampled particle paths. To eliminate that curse of dimensionality, it is desirable to introduce a parametric approximation to p(σ 2 i | x (q) s,0:n−1 , Z s,0:n−1 ) that eliminates the dependence of that function on the particle label q and the simulated sequence x (q) s,0:n−1 . Specifically, we follow the lead in [11,22,32], and, for each i ∈ R, approximate the marginal posteriors p(σ 2 i |x s,0:n−1 )}, q ∈ Q, is a properly weight set available at node s at instant n − 1 to represent p(x 0:n−1 |Z s,0:n−1 ), we make the Monte Carlo approximation On the other hand, from the assumption that p(σ 2 1:R | x (q) s,0:n−1 , Z s,0:n−1 ) is a separable function factored as in (17), it follows that and, therefore, In the sequel, recall that if σ 2 ∼ IG(α, β), then the respective mean and variance of σ 2 are given by [19] Var Therefore, the parametersα s,i,n−1 andβ s,i,n−1 such that IG(σ 2 i |α s,i,n−1 ,β s,i,n−1 ) matches the mean and variance associated with the PDF on the right-hand side of (26) are found, following the procedure in [11,22,32] by making where s,0:n−1 , we get, at node r at instant n, new factorsλ i,n (.) such that where again for all q ∈ Q. The modified importance weight update rule at node r at instant n now becomes Inter-node communication cost By combining the GMM approximation and the moment-matching approximation described before, node s now transmits to its (randomly chosen) neighbor r only the GMM model parameters (15K real numbers as previously explained) plus 2R hyper-parameters (α s,i,n−1 ,β s,i,n−1 ), i ∈ R, as opposed to R × (Q + 1) hyper-parameters as before in the exact RB ReDif-PF algorithm.
Summary of the approximate RB-ReDif-PF Algorithm 1 summarizes the approximate RB ReDif-PF tracker at node r at instant n. In Algorithm 1, the symbol r,n denotes the set η (k) r,n , α r,i,nβr,i,n for i ∈ R and k ∈ K. Algorithm 1 Approximate Rao-Blackwellized random exchange diffusion particle filter 1: procedure REDIF-PF(z r,n , s,n−1 ) 2: Send z r,n to neighbors i ∈ N r 3: Block until receive all z i,n from nodes i ∈ N r 4: s,n−1 from s,n−1 , k ∈ K and resample x (q) s,n−1 , for all q ∈ Q from the GMM approximation defined by those parameters.

Differences between ReDif-PF and the Markov chain distributed particle filter
An alternative and different approach to distributed particle filtering is the MCDPF algorithm introduced in [9]. MCDPF, like other previous work in the distributed PF literature, assumes conditional independence of the sensor observations given the target state and, therefore, should be compared to the proposed ReDif-PF algorithm in this paper in the known sensor parameters scenario of Section 4.1 as opposed to the more general Rao-Blackwellized version of ReDif-PF proposed for unknown sensor parameters in Section 4.2.
The main idea in MCDPF is to move each particle and its associated weight multiple times between nodes in the time interval between instants n and n + 1, according to a Markov chain with transition probabilities defined by the normalized adjacency matrix A of the graph G that defines the network topology. Each time a given particle x n visits a network node r, its weight is multiplied by the pseudo-likelihood p(z r,n |x n ) 1/J φ(r) where φ(r) is the long-term stationary probability of the state of the Markov http://asp.eurasipjournals.com/content/2014/1/19 chain specified by A being equal to r, r = 1, . . . , R, and J is total number of Markov chain move steps between consecutive sensor measurements, which is set by the user. Since the number of visits to the node r divided by J converges to φ(r) [9] as J → ∞, it follows that if J is large enough so that particle x n not only visits all network nodes but also visits each node multiple times, then the aggregate update factor for its corresponding weight at the end of the random walk will approach R r=1 p(z r,n |x n ), which, under the assumption of conditional independence of the sensor measurements given the target state, is the exact update factor for the optimal global weight associated with particle x n . For a finite and especially low number of move steps, MCDPF is no longer optimal, meaning that the choice of the parameter J involves a tradeoff between inter-node communication cost and state estimation error. Contrary to MCDPF, the proposed ReDif-PF does not attempt to compute the exact optimal global posterior PDF p(x 0:n |z 1:R,0:n ) at all nodes r = 1, . . . , R at each instant n. Instead, as explained in previous sections, ReDif-PF builds at each node r and at each instant n a Monte Carlo representation of the posterior p(x 0:n |Z r,0:n ), where Z r,0:n is a random subset of z 1:R,0:n that changes from node to node. Such Monte Carlo representation is built in a way that between instants n and n + 1, each node makes only one request to exchange particles/weights (or equivalent parametric approximations of posterior distributions) with a randomly chosen neighbor, thus eliminating the need for multiple iterative inter-node communication between consecutive sensor measurements and resulting in a communication cost that is much lower than that of the MCDPF algorithm for a similar mean square state estimation error (see the numerical simulation results in Section 5.2).
Finally, we also note that compared to the non-iterative ReDif-PF, MCDPF is also computationally more intensive since each node r has to compute the local likelihoods p(z r,n |x (q) n ) for all its particles x (q) n multiple (namely J) times between instants n and n + 1. We also illustrate that point in the numerical simulations of Section 5.2.

Simulation results
We assessed the performance of the proposed algorithms using 100 Monte Carlo runs with simulated data in three distinct scenarios assuming both unknown and known sensor variances. In all scenarios, we used R = 25 RSS sensors with parameters P 0 = 1 dBm, d 0 = 1 m, ζ r = 3, ∀r ∈ R, and σ 2 r independently sampled at each node according to an IG distribution with mean 16. The nodes were deployed on a jittered grid within a square of size 100 m × 100 m. In the fully distributed algorithms, each node communicates with other nodes within a range of 40 m. All particle filters used Q = 500 particles. Figure 3 shows It also depicts the available network connections. The diameter of the sensor network is D = 5 hops and the minimum number of neighbors for any possible node is 3.

Scenario I: ReDif-PF vs. CbPF
In the first scenario, we assumed unknown sensor variances and evaluated the performance of the Rao-Blackwellized ReDif-PF and two consensus-based PF trackers using respectively iterative minimum consensus (CbPFa) and flooding (CbPFb) (see also [11]). The aforementioned algorithms were compared to the equivalent broadcast implementation of the optimal centralized PF tracker, referred to as DcPF in [11] and [27] and in Section 3.1 of this paper. We also assumed Gaussian priors with mean aforementioned Gaussian prior and, then, converted from polar to Cartesian coordinates. Figure 4 shows the evolution of the root mean square (RMS) error norm -averaged over all network nodes and Monte Carlo runs -of the emitter position estimates for the RB ReDif-PF and the CbPFa and CbPFb algorithms superimposed to the benchmark RMS error curve for the optimal DcPF algorithm. Furthermore, we also show in Figure 4 the average RMS error norm for the non-cooperative (isolated node) trackers and for a local cooperation scheme. In the former, each node runs a regularized PF tracker (see [11]) which assimilates local measurements only, while in the latter, a node r incorporates all measurements Z r,n in its vicinity in the same way as in the ReDif-PF tracker, but it does not exchange its updated posterior with its neighbors. The bars shown in Figure 4 represent the standard deviation of the error norm across all nodes in the network. There are no bars for the DcPF and CbPF algorithms since they provide the same state estimate at all nodes. The RMS error norm at time step 0 for all algorithms was calculated after the measurements z 1:R,0 were assimilated. We implemented the RB ReDif-PF in this scenario with the parametric approximations in Section 4.3 using only one Gaussian mode to represent p(x n−1 |Z s,0:n−1 ).
As expected, CbPFa and CbPFb match the performance of the DcPF tracker since both algorithms reproduce the optimal centralized PF tracker exactly, albeit with different communication and computational costs. On the other hand, as shown in Figure 4, the RB ReDif-PF tracker has a performance degradation compared to DcPF. This result is again theoretically expected since, in the RB ReDif-PF algorithm, the posterior at each node assimilates just a subset of the available measurements z 1:R,n in the whole network at each time step n. However, ReDif-PF offers an improvement in error performance compared to the local cooperation scheme by better diffusing the information across the network. We also note from Figure 4 that the standard deviation of the state estimate across the different network nodes is much lower in the ReDif-PF algorithm than in the local cooperation scheme. Note also that, as shown in Figure 4, isolated nodes were not able to properly track the emitter in the evaluated scenario.
As expected, as σ accel increases, there is a deterioration in the RMS error performance. However, the ratio between the RMS error performance of the suboptimal ReDif-PF tracker and the benchmark optimal DcPF/CbPFb algorithms remains approximately constant (close to a factor of two) along the simulation period for all three different values of σ accel employed.

Communication and computation cost
Considering a four-byte and a one-byte network representation respectively for real and Boolean values, the total amount of bytes transmitted and received by all nodes over the network was recorded while running each tracker in Figure 4. Table 1 summarizes the communication cost for each algorithm in the first scenario (unknown sensor variances) in terms of average transmission (TX) and average reception (RX) rates per node and also quantifies the processing cost for each algorithm in terms of average duty cycle per node, measured in a Intel Core i5 machine with 4GB RAM. The duty cycle of a given node is defined as the  ratio between the total node processing time and the simulation period 100 s. Finally, values in Table 1 are averaged over all Monte Carlo simulations. As shown in Table 1, the RB ReDif-PF tracker with the parametric approximations in Section 4.3 using only one Gaussian mode has a communication cost based on TX rate that is approximately one order of magnitude lower than the flooding-based CbPFb's communication requirements. Compared to the iterative minimum consensus solution (CbPFa), the average communication cost is reduced by two orders of magnitude.

Scenario II: ReDif-PF vs. ReDif-EKF
In the second scenario, the sensor variances are perfectly known and the ReDif-PF tracker is compared both to the optimal centralized PF and to a linearized random exchange extended Kalman filter (ReDif-EKF), which is summarized in Appendix 3. In the simulations, we assumed a non-informative prior for the sensor's initial position that is uniform in the entire surveillance space. The actual initial position of the emitter was, however, sampled from a Gaussian distribution centered at (5 m, 5 m) with standard deviation of 3 m in both dimensions. Figure 6 shows a normalized contour map for the posterior PDF p(x 0 , y 0 |z 1:R,0 ) at instant 0 as a function of x 0 and y 0 assuming the aforementioned non-informative prior. As seen from Figure 6, the initial posterior distribution of the target's position is non-Gaussian. Figure 7 shows the evolution of the RMS error norm assuming known sensor variances respectively for the ReDif-PF algorithm in Section 4.1 with a two-Gaussian GMM parametric approximation and the ReDif-EKF algorithm in Appendix 3. We also show the RMS curve for the optimal centralized PF tracker as a benchmark. The plots in Figure 7 show that, especially in the initial time steps, when the posterior distribution of the states is strongly non-Gaussian as suggested by Figure 6, the fully distributed ReDif-PF outperforms its linearized counterpart, the ReDif-EKF. As the emitter moves away from the near field of the initial dominant sensor, the performance of the ReDif-EKF slowly improves and approaches that of the ReDif-PF, albeit still with a slight degradation towards the end of the simulation.  Table 2 summarizes the communication and processing cost per node for each algorithm in the second scenario. As expected, the DcPF algorithm assuming known sensor variances has the same communication requirements as in the scenario with unknown variances since DcPF locally computes the likelihood functions and then broadcasts them to the entire network. However, as shown in Table 2, DcPF has a slightly lower processing cost when the sensor variances are known. The ReDif-PF tracker on the other hand outperformed the ReDif-EKF tracker in terms of the position RMS error at the expense of a greater communication and computational cost. However, as indicated in Table 2, the communication requirements  of the ReDif-PF and ReDif-EKF trackers still have the same order of magnitude.

Scenario III: ReDif-PF vs. MCDPF/selective gossip
In the third scenario, the ReDif-PF tracker is compared to two iterative algorithms from the literature -the MCDPF and the selective gossip from [9] and [23], respectivelyassuming perfectly known sensor variances as in the second scenario and the same Gaussian priors for the emitter's initial position and velocity used in the first scenario. Figure 8 shows the evolution of the RMS error norm assuming known sensor variances for the ReDif-PF algorithm in Section 4.1 with a single-mode GMM parametric approximation and the MCDPF algorithm in [9] for J ∈ {10, 30, 50, 100} iterations. Figure 9 shows the evolution of the RMS error norm for the ReDif-PF algorithm in Section 4.1 with a singlemode GMM parametric approximation and the selective gossip algorithm in [23] using respectively J ∈ {1, 000; 2, 000; 4, 000} iterations. More specifically, we first run J average gossip iterations considering only the particles in the top 10% bracket in terms of log-likelihood for each randomly selected pair of nodes at each iteration and, subsequently, we run J standard max gossip iterations for the averaged log-likelihood of the selected particle as proposed in [23] to ensure that all nodes have exactly the same weight update factors. Note that, since only one pair of nodes is active at each average gossip iteration and only 10% of the particles are being transmitted between the active nodes, the Selective Gossip algorithm has a lower inter-node communication cost than MCDPF even when a much larger number of iterations is used between consecutive sensor measurements. Table 3 summarizes the communication and processing cost per node for each algorithm in the third scenario. The MCDPF and the selective gossip algorithms have a RMS error performance similar to the ReDif-PF algorithm for J = 30 and J = 4, 000 iterations, respectively, at the expense of a communication cost approximately two orders of magnitude larger than that of the ReDif-PF tracker. Moreover, for a comparable RMS error, the measured ReDif-PF duty cycle is also approximately five and seven times lower than the duty cycle of the MCDPF and the selective gossip algorithms respectively. Note, however, that the selective gossip tracker converges to the same estimate at all nodes and the estimates at each node provided by the MCDPF tracker have a lower standard deviation than those provided by the ReDif-PF algorithm. We also note from Table 3 that with J = 100 Markov chain move steps between sensor measurements, the MCDPF RMS error approaches the error curve of the optimal flooding-based CbPFb tracker with a inter-node communication cost that is, however, roughly four times greater than that of the CbPFb algorithm.

Conclusions
We introduced in this paper a Rao-Blackwellized version of the random exchange diffusion particle filter which enables fully distributed tracking of hidden state vectors in cooperative sensor networks with unknown sensor parameters. Although the general structure of the algorithm can be generalized to arbitrary signal models, we specified the algorithm in this particular paper in an application where we track a moving emitter using multiple RSS sensors with unknown noise variances. The ReDif-PF tracker, introduced originally in a simpler version in [17], is based on random information dissemination and is well suited for real-time applications since, unlike consensusbased approaches, it does not require iterative inter-node communication between measurement arrivals.
The new Rao-Blackwellized version of the ReDif-PF was compared to an exact broadcast implementation of the optimal centralized PF solution, referred to as the DcPF algorithm, and to two equivalent, fully distributed PFs using respectively iterative minimum consensus (CbPFa) and flooding (CbPFb). As expected, due to its suboptimality, the ReDif-PF tracker showed a degradation in RMS error performance compared to both DcPF and the equivalent consensus implementations in our simulations, but required much lower communication bandwidth with savings of one order of magnitude compared to DcPF and CbPFb in terms of transmission rate, and two orders of magnitude compared to CbPFa. The communication cost savings in the RB ReDif-PF algorithm were possible due to suitable parametric approximations introduced in Section 4.3.
The RB ReDif-PF algorithm RMS error performance was also compared in the unknown variance scenario to a local cooperation scheme in which each node assimilates all available measurements in its neighborhood but does not exchange its posterior statistics with other nodes. By diffusing information over the network, the RB ReDif-PF tracker showed better error performance than the local cooperation scheme that uses local information only. Additionally, the standard deviation of the error norm considering all nodes in the network was much lower for RB ReDif-PF than in the local cooperation scheme, suggesting possible weak consensus.
Next, in a second scenario with perfectly known variances, we also compared a non-RB ReDif-PF tracker to its distributed linear filtering counterpart, the ReDif-EKF described in Appendix 3. Due to the non-Gaussianity of the posterior distribution of the states, the distributed PF solution outperformed the distributed EKF solution, albeit, as expected, at a greater computational and communication cost.
Finally, in a third scenario also with perfectly known variances, we compared the non-RB ReDif-PF tracker to two alternative distributed particle filters based respectively on iterative Markov chain move steps between sensor measurements as proposed in [9] and on iterative selective average gossiping as proposed in [23]. In our simulations, the novel ReDif-PF matched the RMS error performance with both the Markov chain and the selective gossip filters with an inter-node communication cost approximately two orders of magnitude lower and a required duty cycle that is reduced by a factor of 5 when compared to MCDPF and a factor of 7 when compared to the selective gossip scheme.
As future work, we plan to extend the ReDif-PF algorithm to perform joint detection and tracking-considering scenarios with probability of detection less than 1 and probability of false alarm greater than 0 as in [33]. We also plan to analyze the diffusion properties of ReDif-PF by investigating the long-term statistical properties of the sequence of visited nodes {r n } , n > 0, defined by the random exchange protocol starting from a random node r 0 .

Appendix 1
In this appendix, we use an importance sampling methodology (see [5,6]) to show that the augmented particle set x (q) r,n } obtained according to (12) and (13) in Section 4 is a properly weighted set to represent the posterior PDF p(x 0:n |Z r,n , Z s,0:n−1 ) in the sense that for any measurable function h(·), s,n−1 , q ∈ Q, be a properly weighted set that represents the posterior PDF p(x 0:n−1 |Z s,0:n−1 ) at node s. Assuming that the particle set x (q) s,0:n−1 was sampled according to some proposal importance function http://asp.eurasipjournals.com/content/2014/1/19 π(x 0:n−1 |Z s,0:n−1 ), the proper weights w (q) s,n−1 may be written as [5,6] w (q) where w(x . Assume next that node s sends its particle set and weights to a neighboring node r that can access at instant n the measurements Z r,n = z r,n ∪ z i,n i∈N r . For any measurable function h(·), we note that and w(x 0:n ) = p(x 0:n |Z r,n , Z s,0:n−1 ) p(x n |x n−1 )π(x 0:n−1 |Z s,0:n−1 ) = p(Z r,n |x 0:n , Z s,0:n−1 ) p(x n |x 0:n−1 , Z s,0:n−1 ) p(x n |x n−1 )p(Z r,n |Z s,0:n−1 ) × w(x 0:n−1 ).

Appendix 3
In a scenario with perfectly known sensor model parameters, assume that at instant n − 1, node s has a linear