A Cooperative Approach to Sensor Localisation in Distributed Fusion Networks

We consider self-localisation of networked sensor platforms, which are located disparately and collect cluttered and noisy measurements from an unknown number of objects (or, targets). These nodes perform local filtering of their measurements and exchange posterior densities of object states over the network to improve upon their myopic performance. Sensor locations need to be known, however, in order to register the incoming information in a common coordinate frame for fusion. In this work, we are interested in scenarios in which these locations need to be estimated solely based on the multi-object scene. We propose a cooperative scheme which features nodes using only the information they already receive for distributed fusion: we first introduce node-wise separable parameter likelihoods for sensor pairs, which are recursively updated using the incoming multi-object information and the local measurements. Second, we establish a network coordinate system through a pairwise Markov random field model which has the introduced likelihoods as its edge potentials. The resulting algorithm consists of consecutive edge potential updates and Belief Propagation message passing operations. These potentials are capable of incorporating multi-object information without the need to find explicit object-measurement associations and updated in linear complexity with the number of measurements. We demonstrate the efficacy of our algorithm through simulations with multiple objects and complex measurement models.


I. INTRODUCTION
F USION networks comprised of geographically dispersed and networked sensor platforms are one of the key enablers of wide area surveillance applications. These networks have the potential to enhance situational awareness in a number of aspects including coverage, accuracy and ease of deployment by providing ad-hoc deployability through a reconfigurable and scalable system structure. The sensor platforms have moderate sensing, computation and communication capabilities and energy resources (as opposed to the rather stringent resource constraints of the commonly considered wireless sensor networks). In order to monitor an unknown number of objects (or, targets), they use sensor signals with an inherently imperfect detection process and obtain noisy measurements from a subset of the objects in their coverage and false detections from the surroundings, all in their sensor centric coordinate system (SCCS). Given the limited bandwidth of the links and the high energy cost of communications, it is often not feasible to forward network wide collected measurements to a designated centre which in turn would have a high computational load [1].
Decentralised paradigms have more desirable properties in fusion networks such as better resource utilisation and flexibility [2]. Typically, the nodes locally filter their measurements to estimate the object trajectories. Then, they exchange the filtered distributions with other nodes over the network to improve upon the accuracy they achieve myopically based on only their local measurements(e.g., [3]). These information messages, however, can be combined only after they are registered in a common coordinate system [4], e.g., the local coordinate frame. Respective sensor locations constitute a fundamental component of sensor registration parameters that specify these coordinate transforms. Geographical routing algorithms underpinning the communication network also rely on a reasonably accurate knowledge of these locations [5].
We are interested in locating the sensors based solely on measurements from the multi-object scene. Such a constraint arise in a range of applications: For example, underwater fusion networks cannot exploit global navigation space systems (GNSS) due to signal propagation constraints of their environment [6]. In terresterial settings, GNSS might fail to perform reliably considering their vulnerabilities to deliberate interferences such as jamming [7]. The use of reference (or, cooperative) vehicles [8] does not match well the flexibility requirements. Another alternative which has been investigated intensely is to use communication front-end and/or network statistics such as received signal strength (RSS) and time of arrival (TOA) [9]. Localisation based on RSS and TOA type noisy distance measurements, however, is often not sufficiently accurate for networks with the degree of connectivity typical in fusion applications [10,Chp.6].
Our problem can be treated as a particular type of sensor registration (or, calibration [11]) using targets of opportunity. The latter topic has been investigated in the context of multitarget tracking, however, mostly for mitigating biases of model parameters [12]. The work which implicitly or explicitly address locating sensors include solutions in centralised settings based on conventional multi-target trackers [13], [14] and random finite sets (RFS) based multi-target filtering [15], [16]. These work involve selection of either a maximum likelihood (ML) or a Bayesian estimation paradigm and specifying the parameter likelihood of the problem (see, e.g., [17,Sec.IV]) in accordance with the multi-object and measurement models. This likelihood, however, requires all the target measurements collected across the network to be filtered together, and, in turn, centralised processing. Distributed alternatives often resort to joint filtering which embodies all the drawbacks of centralised fusion, both in the case of ML [18] and Bayesian [19] paradigms.
In this article, we propose a distibuted online selflocalisation scheme that avoids any form of centralisation and operates alongside distribution fusion. The computational structure is composed of iterative local message passings yielding full decentralisation. The algorithm is developed in two steps: First, we consider a pair of fusion sensors and develop a calibration likelihood which can be computed based solely on the multi-object distributions exchanged and local target measurements. Specifically, we approximate the parameter likelihood with a node-wise separable form which removes the need for local centralisations.
Second, we merge these likelihoods in a pairwise Markov random field (MRF) model which captures the communication structure of the network. Such models are equipped with estimation algorithms which operate in a message passing fashion and have proved useful for distributed estimation [20] and target tracking [21] applications in wireless sensor network as well as self-localisation based on RSS measurements [22]. In our model, we use the proposed node-wise separable likelihoods as edge potentials [23]. These potentials have a timerecursive structure as the respective parameters are indirectly observed through measurements from a hidden process.
In the resulting scheme, nodes simultaneously perform distributed fusion and update node-wise separable likelihoods -equivalently, edge potentials-with their neighbours using the incoming multi-object posteriors and local target measurements. At the end of a fixed length time window, the updated potentials are used in Belief Propagation (BP) [24] message passing iterations. We accommodate and benefit from the rich information from multiple targets of opportunity using Poisson multi-object models [25] in the update recursions. These models are propagated by local multi-object filters from which our scheme inherits the capability of operating with complex measurements involving false alarms and noisy measurements of the objects with imperfect detection rates. The potential updates feature linear complexity with the number of measurements without any need to explicitly find targetmeasurement associations.
The article is structured as follows: We provide the problem statement in Section II. Then, we introduce a pairwise MRF model for calibration and overview distributed estimation of calibration marginals in Section III. In Section IV, we develop node-wise separable likelihoods for calibration problems. We combine these likelihoods with the MRF model introduced and describe our collaborative scheme in Section V. We detail a Monte Carlo (MC) algorithm for self-localisation in Section VI and demonstrate its efficacy in Section VII. Finally, we conclude in Section VIII.
II. PROBLEM STATEMENT We consider networked sensor platforms listed as V " t1, ..., N u. The communication links available between pairs of sensors pi, jq specify the edge set E Ă VˆV with respect to the relation E " tpi, jq|i and j share a communication linku. We assume bidirectional communication links which is captured by using an undirected graph G for which if pi, jq P E, then it holds that pi, jq P E ô pj, iq P E. The neighbours of node i in G constitute the set nepiq fi tj|pi, jq P Eu. The network topology in G is connected and might contain cycles. We employ widely used models for capturing complex and uncertain interactions between sensors and an unknown number of manoeuvring objects.
An object in the network surveillance region is described by a state vector x in the state space X . Typically, x contains the Cartesian coordinates of the object in an Euclidean plane x l and its velocity x v , i.e., x " rx l , x v s.
The evolution of a moving object's state and the measurement it induces are described by a hidden Markov model (HMM). The state at time k is distributed according to the density πpx k |x k´1 q where x k´1 is the realisation of the previous state with initial density π 0 px 0 q 1 . An observation z k,j is generated at sensor j through a likelihood of the form lpz k,j |x k q independently from other sensors.
Sensor likelihoods are specified by models that characterise the detectors processing sensor signals. In multi-sensor settings, it is useful to explicitly condition the likelihoods further on calibration parameters which relate a given target state in a desired reference frame to sensor readings. These parameters might involve respective quantities such as the location and the orientation of the sensors with respect to the reference frame, as well as scale factors.
For example, consider a commonly used model involving a nonlinear mapping of the target state with additive uncertainties: where n j is a measurement noise realisation with probability density g j pn j q, and, h j is a known mapping. Here, rx k s j denotes x k in the SCCS of the jth sensor. For a sensor localisation problem, where θ j is the sensor location and x l k is the position component of the state vector x k , all represented in the reference frame. The likelihood density for this model is then found as l j pz k,j |x k ; θ j q " g j`zk,j´hj prx k s j q˘.
For more general calibration problems, it is useful to introduce a transform from the reference frame to the SCCS of sensor j parameterised by θ j , i.e., rx k s j " τ j px k ; θ j q.
For example, if the respective orientation angle α j in the plane is unknown in addition to the location of sensor j, then it is useful to consider the transform given by where Rpα j q is the rotation matrix for α j and θ l j is the position component of θ j in the reference frame. The likelihood density for this model is then l j pz k,j |x k ; θ j q " g j´zk,j´hj`τj px k , θ j q˘¯.
In practice, it is reasonable to assume that θ j takes values from a bounded set B Ă R d of dimensionality d. For example, d " 2 when localisation in a plane is considered. Our goal is to find these parameters θ fi pθ 1 , ..., θ N q using local message passings on G and based on measurements from the multi-object scene (or, targets of opportunity), models for which are discussed next 2 .
Consider the set of realisations of M k hidden processes X k fi tx 1 k , ..., x M k k u at time k. A point x P X k induces a measurement at sensor j with probability P D,j pxq, independently. Let us denote the set of measurements from the objects at sensor j byZ j k . In addition toZ j k , sensor j collects false detections from the surroundings, or clutter. We use a Poisson process to model the clutter points, i.e., C j is a Poisson realisation denoted by C j " P oisp.; λ C,j , s C,j pzqq where λ C,j is the average number of (Poisson distributed) clutter points and s C,j pzq is their spatial density 3 . Therefore, the set of measurements sensor j receives at time k is given by Sensor calibration, hence, involves finding the joint parameter likelihood relating the measurement histories tZ j 1:k u jPV and θ, i.e., l`Z 1 1:k , ..., Z N 1:k |θ˘. We consider a random θ and use this likelihood to update a prior distribution. Because θ is bounded, a uniform prior p 0 pθq over B N can be used, and, the posterior density is given by ppθ|Z 1 1:k , ..., Z N 1:k q 9 p 0 pθq l`Z 1 1:k , ..., Z N 1:k |θ˘.
We are interested in finding the minimum mean squared error (MMSE) estimate of θ based on this posterior 4 .

A. Centralised sensor calibration
In the centralised estimation of θ, a designated centre has access to all measurement histories tZ j 1:k u jPV to evaluate the parameter likelihood l`Z 1 1:k , ..., Z N 1:k |θ˘" where the factorisation follows from the chain rule of probabilities [17,Sec.IV]. The factors in the right hand side (RHS) are independent contributions of the measurement sets collected at each time step. The Markov property of the hidden object processes admit that the current measurements are conditionally independent of the measurement histories, given the current object states and the sensor locations. Using this relation together with the (conditional) mutual independence of measurements across the sensors, the instantaneous contribution of measurements collected at t`1 to the likelihood at k is found as [17] p` where the product term inside the integral is the multi-sensor likelihood. The second term is a prediction density for the multi-object scene at time t`1 based on the network history until t, and, can be found by the Bayesian filtering recursions, i.e., it is output at the prediction step of a "centralised" filter.
In addition to the difficulties in collecting Z j 1:k s at a designated fusion centre under resource constraints, the estimator obtained by substituting from (8) to (7) and (6) inherits the computational issues related to both multi-sensor filtering and parameter estimation in state space models: The centralised filter needed for the realisation of (8) faces the same computational complexity issues which are encountered in multi-object filtering and which are exacerbated in the case of multiple sensors 5 . The usual non-Gaussian/non-linear dynamics involved in such problem settings suggest Sequential Monte Carlo (SMC) methods to be used as a computational paradigm, and, in the case of parameter estimation problems including calibration, iterative sampling strategies with relatively high computational and memory complexities should be used to achieve reasonably accurate estimates [28]. Thus, the described centralised solution suffers from poor scalability with the number of sensors and objects.
We circumvent these problems through a number of modelling approximations introduced in the rest of this article, and provide a distributed scheme with efficient local computations.

B. Distributed fusion architecture and information exchange
Each sensor platform locally filters its measurement history and maintains a local representation of its environment including the multi-object state X k . There is no restriction on the multi-object models and inference algorithms used by these platforms. The neighbouring nodes in G exchange their representations to improve upon their myopic accuracy. As far as our self-localisation (calibration) scheme is concerned, the only requirement is that a Poisson RFS approximation of X k is available from the information received. In other words, the neighbours provide -directly or indirectly-a model in which X k is a realisation of a Poisson RFS conditioned on the sensor history, i.e., X k " P oisp.; λ k|k , s k|k pxqq, which first draws M k from a Poisson distribution with parameter λ k|k and then generates M k points from s k|k pxq. Hence, X k evaluates its density as where the product selects each element only once. Poisson multi-object models are provided directly by the Probability Hypothesis Density (PHD) filter [29]. In the case that another RFS model is provided (from filtering algorithms such as the Cardinalised PHD [30], labelled [31] and variational [32] multi-Bernoulli filters), a best Poisson approximation can be found using the first order statistical moment of this distribution [25]. More conventional algorithms such as multiple hypothesis tracking (MHT) consider hypotheses on measurement object correspondances, on the other hand, and find the likelihood of these association hypotheses together with single object posteriors. This form of information is convertible to RFS distributions [33] from which Poisson multi-object models can readily be obtained.
Besides that they can be obtained from a wide range of trackers, there are additional appealing features of Poisson multi-object models: because they are RFS models, the corresponding sensor likelihoods evaluate without the need to explicitly find measurement-object associations [25,Eq.(12.41)]. Moreover, marginalisations over the state variable as in the RHS of (8) have closed form expressions and can be computed using Monte Carlo methods with linear complexity in the number of measurements for a single sensor. These features of multi-object scalability are later combined with the MRF model we introduce in Sections III and IV for multi-sensor scalability.

III. A DYNAMIC PAIRWISE MARKOV RANDOM FIELD MODEL FOR CALIBRATION
In this section, we introduce an approximation to the centralised estimator (Eq.s(8), (7) and (6)) which enables the estimation of θ in a distributed fashion, through local message passing operations.
We consider a parsimonious representation for the calibration posterior in (6). Specifically, we assume that θ is Markov with respect to the communication topology G, i.e., if the sets of nodes A and B are separated by C on G, then the random variables associated with A, i.e., θ A " tθ i |i P Au, and θ B are conditionally independent given θ C . This conditional independence relation is often denoted by θ A K K θ B |θ C [23]. All such relations admitted by G factorise (6) to positive functions (or, potential functions) over the cliques of G (connected subsets of V) [23].
The distributed fusion architecture discussed in Section II-B involves nearest neighbour multi-object posterior exchanges. Consequently, the local information available at the nodes are related to the respective locations of pairs of (neighbouring) sensors. Therefore, we take into account only the singleton and two-node cliques and use a pairwise MRF modelp where the node potential functions ψ i s are arbitrary priors for θ i (e.g., uniform distributions over B) and the edge potentials ψ k ij s are predictive parameter likelihoods for the pairs pi, jqs based on sensor histories up to time k. These edge potentials have the time-recursive structure in (7), i.e., (11) and render a dynamical MRF.
One challenge in calibration using this model in distributed fusion networks is the computation of ψ k ij s in a distributed fashion, without communicating Z i k s. Another issue is the computational load even if the measurements could be exchanged among the neighbouring nodes. We introduce nodewise separable likelihoods to tackle these challenges later in Section IV. Next, we give a brief outline of decentralised estimation in sensor networks based on pairwise MRFs.

A. Decentralised Estimation Using Belief Propagation
Consider the MMSE parameter estimator based on the joint posterior in (10). This can be decomposed as a concatenation of the MMSE estimates of θ i s based on the marginal distributions. The pairwise MRF model in (10) allows the computation of the marginal densities through iterative local message passings such as Belief Propagation (BP) [24], when G contains no cycles. Specifically, the nodes maintain distributions over their local variables and update them based on messages from their neighbours which summarise the information neighbours have gained on these variables. This is described by the set of equations for all i P V, where k i s are scale factors ensuring thatp i pθ i qs integrate to unity. The fixed points of the set of equations above are the marginals of (10) in the case of a cycle-free G. In BP iterations, nodes simultaneously send messages to their neighbours using (12) (often starting with constants as the previously received messages) and update their local "belief" using (13). If G contains no cycles,p i s are guaranteed to converge to the fixed points, i.e., the marginals of (10), in a finite number of steps [24].
For the case in which G contains cycles, message and update equations are still well defined. BP iterations on loopy graphs, however, are not guaranteed to convergence to an equilibrium, in general [23, Sec.4]. The fixed points of Eq.s (13) and (12), if they exist, are often not equal to the marginals of the global distribution, but reside in vicinities of them. Additional details on the convergence of BP on graphs with cycles can be found in [23] and the references therein. Loopy versions of BP has been very successful in finding approximate marginal distributions in many applications including fusion problems in sensor network (see, for example [20], and the references therein).
As a result, we distribute the calibration task over the network by using BP iterations to compute the marginals of (10) (exactly or approximately, in the cases of tree or loppy structured Gs, respectively). Since the neighbours in G share a communication link, the messages of BP map directly onto the communication network and the ith sensor localises itself using, e.g., the MMSE rule, with its local marginal.
The next step involves the challenge of designing efficient computational procedures for approximating the edge potentials of the pairwise MRF model, which is discussed next.

IV. NODE-WISE SEPARABLE EDGE POTENTIALS
The MRF model in (10) and (11) does not readily allow for distributed calibration because the potential functions are the centralised likelihoods for sensor pairs (i.e., Eq.s (7) and (8) for sensors i and j). We overcome this problem by introducing an approximation for the instantaneous likelihood, or, update, term in the RHS of (11) which can be computed in a distributed fashion. This approximation has a node-wise separable form as a product of terms based on locally available information such as multi-object models received from the neighbours (for example, posterior multi-object Poisson distributions) and locally collected measurements of the multiobject scene. Specifically, we approximate the edge potential with a product of the form where l k ij is computed at node i in a recursive fashion and based on only Z i 1:k and the incoming posteriors from sensor j and vice versa.
Let us consider the update term in the RHS of (11) and denote the pair pθ i , θ j q by θ i,j . This term can be further factorised in alternative ways as follows: In the first and second lines above, the chain rule is used and the third equality follows from taking the geometric mean of the first two expressions. The conditioning of these factors to the measurement histories of both sensors prevents decentralisation. An alternative is to approximate the four factors in Eq.(15) by leaving out the history of sensor i (sensor j) in conditioning of the first and last (second and third) terms of (15) together with any current measurements, i.e., where κ above is the normalisation constant and it can easily be shown that κ " 1.
The appeal of this approximation is that its factors depend on single sensor histories allowing us to avoid centralisation. The computation of these factors in terms of local information and incoming multi-object posteriors from the neighbours and the utilisation of the update in (17) in a message passing scheme are detailed in Section V. Now, we consider the approximation quality in terms of the Kullback-Leibler (KL) divergence [34] between the centralised update term in the left hand side (LHS) of (17) with respect to its approximation on the RHS: Proposition 4.1: The KL divergence between the centralised update in (8) and the node-wise separable approximation in (16) equals to a sum of Mutual Information (MI) [34] terms given in (18).
The first two MI terms in Proposition 4.1 measure the departure of the current measurement and the sensor history from conditional independence when they are conditioned on the other sensor history instead of the current state X k for which Z i k K K Z i 1:k´1 |X k , θ i,j and Z j k K K Z j 1:k´1 |X k , θ i,j hold leading to The third term is a measure of departure from conditional independence for the current measurements given the measurement histories. Therefore, the RHS of (18) is zero if Z i k , K KZ i 1:k´1 |Z j 1:k´1 , θ i,j and Z j k K K Z j 1:k´1 |Z i 1:k´1 , θ i,j hold together with Z i k K K Z j k |Z i 1:k´1 , Z j 1:k´1 , θ i,j simultaneously. This condition is satisfied, for example, in the case that either of the measurement histories Z i 1:k´1 and Z j 1:k´1 are sufficient statistics for X k and the true state can be predicted by both sensors with probability one. One should not expect this level of prediction accuracy in realistic tracking scenarios, therefore, it is instructive to relate the KL divergence in (18) further to the uncertainty on X k given the sensor histories.

Corollary 4.2:
The KL divergence term in (18) is upper bounded by the difference between the total local state prediction entropies and the entropies of the joint prediction and its most uncertain single sensor update given in (19) with H denoting the Shannon Entropy [34]. Proof. See Appendix B.
Corollary 4.2 relates the approximation quality of the nodewise separable updates to the uncertainties in the object state prediction and estimation. The first two terms in the RHS of (19) measure the uncertainties in the locally predicted object states. Subtracted from their sum are the uncertainty in the joint prediction based on the histories of both of the sensors and the entropy of the single sensor update upon the joint prediction that has the highest value. Therefore, a better quality of approximation should be expected as the local prediction densities become more concentrated around a single point in the state space.
Consider, for example, range-bearing sensors providing noisy measurements of object positions in polar coordinates. Tracking filters can provide fairly accurate predictions and estimates of object locations and velocities using typical measurements and with increasing k. This in turn results with smaller values for the bound in (19). An alternative in which these conditions cannot be satisfied with guarantees is a bearing-only sensor scenario consisting of a single object and two identical sensors. Suppose that the object moves along the line-of-sight (LOS) of one of the sensors. The local target prediction as well as the posterior distribution at time k will typically have a probability mass spread around the line-ofsight. The amount of this spread, which can be measured by the entropy H, will drop significantly when state distributions are conditioned also on the other sensor's history [35] yielding a large value for (19).
The use of the node-wise separable term in (16) to update the dynamic MRF edge potentials given by (11) leads to the following recursive formulae: where the node-wise terms in (14) are the products of individual node-wise separable update factors over time defined in a recursive fashion: l k ji pθ i , θ j q fi l k´1 ji pθ i , θ j qppZ j k |Z i 1:k´1 , θ i,j q.
Here, subscript ij indicates that the associated term is computed at sensor i and will be transmitted to sensor j and vice versa.
V. THE COOPERATIVE CALIBRATION SCHEME In this section, we describe the cooperative calibration scheme based on the node-wise separable edge potentials. First, in order to facilitate online processing within the Bayesian paradigm, the measurement histories are partitioned into time windows of length T , and, θ is evolved with artificial dynamics between two consecutive windows leading to prediction-update cycles [28]. Let us denote the nth window of sensor j by Z j n fi Z j pn´1qT`1:nT . A recursive rule is, then, given by p n pθ n |Z 1 0:n , ..., Z N 0:n q9 lpZ 1 n , ..., Z N n |θ n qp n|n´1 pθ n |Z 1 0:n´1 , ..., Z N 0:n´1 q, p n|n´1,i pθ n,i |Z 1 0:n´1 , ..., Z N 0:n´1 q " ż β n pθ n,i |θ n´1,i qp n´1,i pθ n´1,i |Z 1 0:n´1 , ..., Z N 0:n´1 qdθ n´1,i , for i " 1, ..., N where subscript i denotes the ith marginal.
In the second line we apply dynamics to the single sensor parameters θ n´1,i independently which is specified by the conditional density β n pθ n,i |θ n´1,i q. β n models (small) Brownian motion steps, i.e., β n pθ n,i |θ n´1,i q " N pθ n,i´θn´1,i ; 0, D n q, where N is a multi-dimensional Gaussian with an all zero mean vector. D n " diagpσ 2 n,1 , ..., σ 2 n,d q is a dˆd diagonal covariance matrix and σ n,1 , ..., σ n,d are the Brownian motion parameters specifying the step sizes in each direction. Consequently, the prior distibution in (23) becomes the product of these marginals, i.e., p n|n´1 pθ n |Z 1 0:n´1 , ..., Z N 0:n´1 q " ź iPV p n|n´1,i pθ n,i |Z 1 0:n´1 , ..., Z N 0:n´1 q.
The prediction step in (24) can already be performed locally. The discussion in Sections III and IV relates to the update step in (23). In other words, we replace (23) with p n pθ n |Z 1 0:n , ..., Z N 0:n q9 ź pi,jqPEψ k i,j pθ n,i , θ n,j q ź iPV p n|n´1,i pθ n,i q. (26) We describe the cooperative update of the marginal distributions for a time-window of length T . Then, these steps are repeated in consecutive time windows.
The conceptual steps followed at platform i is given as pseudo-code in Algorithm 1. Initially, platform i specifies an a priori distribution over its calibration parameters and the multiobject scene. Here, k becomes a time index for traversing the current window n. In an infinite loop, platform i proceeds with distributed fusion by first collecting measurements from the multi-object scene and filtering them to find the local multiobject density f i pX k |Z i 1:k q. This density is exchanged with the neighbouring platforms and then fused with the incoming state densities from the neighbours.
For cooperative calibration, the local node-wise update term ppZ i k |Z j 1:k´1 , θ i,j q is found for all neighbours j P nepjq. The Filter Z i k and find the Poisson multi-object density f i pX k |Z i 1:k q Ź Local filtering

7:
Exchange f i pX k |Z i 1:k q s with j P nepiq Ź Distributed fusion 1

8:
Fuse f i pX k |Z i 1:k q and f j pX k |Z j 1:k qs usingθ i andθ j s Ź Distributed fusion 2

9:
Find ppZ i k |Z j 1:k´1 , θ i,j q using Z i k and f j pX k´1 |Z j 1:k´1 q for all j P nepiq Ź Eq.s (27) and (29) 10: Update l k ij with ppZ i k |Z j 1:k´1 , θ i,j q for all j P nepiq Ź Eq. (21) 11: if k " T then

12:
Exchange l k ij s with j P nepiq 13: Findψ k ij for j P nepiq Ź Eq. (14) 14: Perform loopy BP for S steps to findp n,i pθ i q Ź Eq.s (12) and (13) to estimate the ith marginal of (26).

21: end while
computation is carried out in two steps: First, based on the Poisson density received at k´1, i.e., f j prX k´1 s j |Z j 1:k´1 q " P oisprX k´1 s j ; λ k´1,j , s k´1,j prxs j qq, the corresponding prediction density is found. Here, we explicitly indicate that the argument of the incoming density is in the SCCS of sensor j. Poisson predictions often incorporate a probability for an object with state x to continue to exist at the next time step denoted by P S pxq. The predicted Poisson model is characterised by [25] λ k|k´1,j pθ i,j q " λ k´1,j ż P S pτ i˝τ´1 j px; θ i,j qqŝ k´1,j pxqdx, s k|k´1,j prxs i ; θ i,j q 9 ż πprxs i |τ i˝τ´1 j px 1 ; θ i,j qqP S pτ i˝τ´1 j px 1 ; θ i,j qqs k´1,j px 1 qdx 1 , (27) where π is the Markov transition modelling object motion (Section II). Here, the integral variables are in the jth SCCS and the composite function τ i˝τ´1 j transforms its argument in the jth SCCS to its counterpart in the ith SCCS given the calibration parameters θ i,j .
The predictive Poisson model f j pX k |Z j 1:k´1 ; θ i,j q specified by (27) has its X k argument in the ith SCCS. As a second step, the node-wise update term is computed based on this predictive density and the current measurements. Using the conditional independence relation Z i k K K Z j 1:k´1 |X k , it can easily be shown that where the integral variable is set valued and (28) is a set integral (see Appendix C for an overview of such integrals).
In addition, the RHS of (28) takes a simple form for a Poisson f j and the sensor measurement model described in Section II. In particular, ppZ i k |Z j 1:k´1 , θ i,j q " exp´´λ C,i´λk|k´1,j pθ i,j q ż P D,i pxq s k|k´1,j px; θ i,j qdx¯ź zPZ i k˜λ C,i s C,i pzq`λ k|k´1,j pθ i,j qż P D,i pxq l i pz|xq s k|k´1,j px; θ i,j qdx¸, (29) where subscript i denotes that the quantity belongs to the measurement model of node i (see, e.g., [36,Proposition 11.3], [29,Eq.(116)], or [16,Eq.(8)]). This term is then used to update the likelihood factor l k ij in (21). At the end of the time window of length T , the nodewise terms are exchanged with the neighbouring nodes and the edge potentials are constructed using (14). Then, S steps of loopy BP is performed using (12) and (13) to estimate the ith local marginal density of the MRF calibration model in (26). In the next step, the expected value of this marginal is found which becomes the current estimate of the parameters, and is exchanged with the neighbours. Before proceeding with the next time window, the local parameter distribution is convolved with a zero mean Gaussian realising the Brownian motion step given by (24) and (25).
It is worth noting that (28) is valid for arbitrary sensor coverages: If different sensor pairs observe different targets in common due to partially overlapping coverages, it suffices to select the probability of detection P D,i pxq embedded in the likelihood term as zero outside the coverage. The explicit evaluation given by (29), however, is valid when there are no targets observed by sensor i but not by sensor j. For arbitrary coverages, (29) would have the same form with the spurious measurement model (i.e., the terms related to clutter) adapting to include observations from those targets in the coverage of sensor i but not j. Further elaboration on this case remains as future work. For the rest of this article, we restrict our attention to the setting in which all sensors observe the same targets.

VI. COOPERATIVE SENSOR SELF-LOCALISATION USING MONTE CARLO METHODS
This section describes a realisation of the cooperative calibration scheme introduced in the previous section for sensor self-localisation using Monte Carlo methods. Central to the realisation of the conceptual procedure in Algorithm 1 are particle representations of the continuous densities involved and approximate computations which in turn require slight modifications to the original steps.
For sensor self-localisation, we consider a co-planar network and the respective calibration parameters as the sensor locations in a network coordinate system. An arbitrary node is selected as the centre of the network coordinate frame. For example, the SCCS of node 1 can be selected as the reference by setting θ 1 " r0, 0s T , without loss of generality. Thus, the transform in (2) maps the points in the 1st SCCS to their jth SCCS counterparts.
Consider the pair pi, jq P E. The composite function τ i˝τ´1 j appearing in (27) is found as Local filtering can be performed using, for example, Sequential Monte Carlo methods [37] within the filtering approaches discussed in Section II-B such as Sequential Monte Carlo realisations of RFS filters [38]. Using its local filter, node j will have transmitted an approximation to P oisp.; λ k´1,j , s k´1,j pxqq specified by an estimateλ k´1,j for the expected number of objects and a weighted set of independent, identically distributed (i.i.d) samples (or, particles) tx pmq k´1,j , ζ pmq k´1,j u M m"1 that encode an empirical distribution given byŜ where δ x is the Dirac measure concentrated at x. The predictive Poisson model in (27) is then approximated by substituting the empirical distribution in (31) together withλ k´1,j into (27), and using the Monte Carlo integration principle [39,Chp.3]. This is akin to the prediction stage of the Sequential Monte Carlo PHD filter [40] and the resulting approximation (also considering (30)) is given bŷ λ k|k´1,j pθ i,j q " λ k´1,j ÿ m ζ pmq k´1,j P S px pmq k´1,j´θ j`θi q, x pmq k|k´1,j " πpx|x pmq k´1,j´θ j`θi q, ζ pmq k|k´1,j " ζ pmq k´1,j P S px pmq k´1,j´θ j`θi q ř m 1 ζ The update term in (29) is computed using these quantities which yieldŝ The likelihood update ppZ i k |Z j 1:k´1 , θ i,j q local to node i can only be estimated for a finite set of θ i,j values which in turn will be used to estimate the node-wise likelihood term l k ij . It is beneficial to use the same set of parameter values at node j in order to estimate l k ji as the edge potentialψ k ij , then, can be found by simply taking the product of these estimates. For this reason, we generate L equally weighted samples from p n|n´1,i pθ i q and p n|n´1,j pθ j q, and, obtain tθ plq n|n´1,i u L l"1 and tθ plq n|n´1,j u L l"1 , respectively 6 . We use the concatenation of these values, i.e., in computing (33) at both node i and j.
In order to construct (34), the nodes exchange their local particle sets representing p n|n´1,i pθ i qs as an additional communication step to Algorithm 1, before starting to process a new time window. Thus, at each time step k within the window, the update term (33) is computed for all tθ plq i,j u L l"1 and for all j P nepiq in order to update tl k i,j pθ plq i , θ plq j qu L l"1 . At the last step of the time window k " T , the estimated node wise terms are exchanged and the edge potentials tψ k i,j pθ plq i , θ plq j qu L l"1 are found by simply taking the element-wise product of the nodewise separable terms (Eq.(14)).
The loopy BP realisation with these edge potentials follows the non-parametric BP (NBP) approach [41] and represents both the messages in (12) and the local marginals in (13) by particle sets.
First, we describe the message computations: Consider (12) for the MRF model in (26) and suppose that i.i.d samples from the (scaled) product of the jth local prior and the incoming messages from all neighbours except i are given, i.e., Based on these samples, the message from node j to i (scaled to one) is approximated by a smoothed density estimate in NBP [41]. We use Gaussian kernels leading to the approximation given bŷ where the kernel weights are the normalised edge potentials. Λ ji is related to a bandwidth parameter that can be found using Kernel Density Estimation (KDE) techniques. In particular, we use the rule-of-thumb method in [42] and find wherem ji andĈ ji are the empirical mean and covariance of the samples, respectively.
Second, we consider sampling from the updated marginal in (13) for the MRF in (26). We use the weighted bootstrap, or, sampling/importance resampling, approach [43] with samples generated from the (scaled) product of Gaussian densities with mean and covariance equal to the empirical mean and covariance of the particle sets, respectively. In other words, givenm ji andĈ ji as above, we generate θ plq i,β " βpθ i q, l " 1,¨¨¨, L, βpθ i q 9 N pθ i ;m n|n´1,i ,Ĉ n|n´1,i q ź jPnepiq N pθ i ;m ji ,Ĉ ji q.
The particle weights for these samples to represent the updated marginal is given by wherep n|n´1,i is a KDE found using tθ plq n|n´1,i u L l"1 with the rule-of-thumb method described above 7 . Thus, the local calibration marginal is estimated bŷ As the final step of the bootstrap, tθ plq i,β ,ω plq i,β u M l"1 is resampled (with replacement) leading to equally weighted particles from p n,i pθ i q, i.e., tθ plq n,i u L l"1 . We follow similar bootstrap steps in order to generate the samples in (35).
After nodes iterate the BP computations described above for S times, each node estimates its location by finding the 7 If p n|n´1,i can be evaluated exactly, it is replaced withp n|n´1,i . empirical mean of tθ plq n,i u L l"1 . Before proceeding with the next time window, Brownian motion is applied to tθ plq n,i u L l"1 yielding i.i.d. samples from p n`1|n,i pθ i q, i.e., θ plq n`1|n,i " θ plq n,i`ǫ plq n , ǫ plq n " N p.; 0, σ 2 n Iq, (38) where I is the 2ˆ2 identity matrix. The communication messages in this implementation strategy for the conceptual procedure in Algorithm 1 consist of particle sets which are equivalently (unordered) numerical arrays. For distributed fusion, node i broadcasts an array of average length OpM k M q to its neighbours in line 7 of Algorithm 1, where M k is the number of targets and M is the number of particles used per target.
The exchanges for cooperative localisation have a period of T . First, nodes broadcast L particles, i.e., an OpLq array, to construct (34). Because these particles are generated from the a priori localisation distributions in lines 3 and 16, it is convenient to have these broadcasts take place in these lines. Hence, the expectation in line 3 can be performed at the neighbours using the empirical mean of the messages. After the node-wise separable terms are computed, they are exchanged in line 12 which corresponds to the transmission of an OpLD i q array where D i is the degree (or, the number of neighbours) of node i. The loopy BP messages in line 14 are encoded by the L particles given in (35). The mixture representation in (36) is recovered from these particles by using the weigths already known both at the transmitting and the receiving nodes, and, the rule-of-thumb KDE method. As a result, S iterations of loopy BP corresponds to the transmission of an OpSLD i q array. Note that, the exchange of location estimates from the current posterior marginals in line 15 leads to a comparably negligible, constant communication load. As a result, the communication complexity of the cooperative scheme for node i is OpLppS`1qD i`1 qq.

VII. EXAMPLE
We demonstrate the cooperative self-calibration algorithm introduced in Section V for self-localisation using the Monte Carlo realisation described in Section VI. The example scenario is depicted in Figure 1 and consists of four objects moving in a surveillance region being observed by nine sensors for 200 steps 8 . The nodes networked through communication links which are the edges of the communication graph G (blue lines).
The object states are described by their planar position and velocity, i.e., x n " rpx l n q T , px v n q T s T . The trajectories are obtained using a linear constant velocity motion model with (slight) additive process noise: πpx n |x n´1 q " N px n ; F x n´1 , Qq Results using this scheme on a problem with a smaller scale can be found in [44]. where I and 0 are the 2ˆ2 identity and zero matrices, respectively.
The sensors measure the bearing angle and the range of the objects in their SCCS with standard deviations σ φ " 1˝and σ R " 5m, respectively. The likelihood model in (5) is then found as Each object is detected with P D " 0.98. The number of false alarms at time k is Poisson distributed with mean λ C,i " 2 for i P V and the associated spatial distribution is uniform in the sensor field-of-views. A typical realisation of measurements at sensor 7 over time can be seen in Figure 2.
The nodes filter their local measurements using the SMC PHD algorithm in [45]. The probability that an object with state x remains to exist in the next step is select P S pxq " 0.9. The resulting Poisson multi-object models are exchanged with the neighbours in G (for distributed fusion).
Node 1 is selected as the origin of the network coordinate system. For the other nodes, the localisation prior for the first time window, i.e., p 1|0,i pθ i q, is selected to be a uniform distribution over the sensing region. A window length of T " 10 is used for computing the node-wise separable edge potentials as detailed in Section VI. The first time window starts at k " 40 to avoid using the posteriors from the early (transient) stages of filtering. Because uniform localisation priors are used, the initial set of location values in (34)at which the node-wise terms are evaluated-are selected The field-of-view is the circular region around the sensor with a radius of 4000m. as random permutations of an L " 900 node uniform grid over the sensing region. Thus, at each time step k within the window, the local likelihood updates for these values are computed using (32)-(34) and the node-wise terms (21), (22) are updated. At the end of the time window, these terms are exchanged to find the edge potentials (20). This stage is followed by S " 8 steps of nonparametric BP (Eq.s (35)-(37)).
Following the last NBP step, each node estimates its location as the empirical mean of the local marginal represented by L " 900 particles. A total of 16 iterations of the proposed algorithm is performed starting from k " 50 to 200 in every T " 10 steps. The Brownian motion parameter σ n in (38) varies from σ 1 " 600 to σ 16 " 1 in accordance with a square law.
We consider the maximum localisation error in the network for 200 Monte Carlo runs. In Fig. 3 we present the box-plot of these errors for all runs with respect to the time window index n. The average error at the final iteration n " 16 is 5.95m which is comparably close to the uncertainties in sensor measurements. The highest error, or, the error margin is 12.09m. Next, we consider the error margins for all steps and normalise them with the minimum distance between the nodes in this deployment, i.e., 1000m. The semi-log plot of the normalised error margins can be seen in Fig. 4. Note that the decay speed is close to a log-linear regime and reaches a minimum corresponding to approximately 1.2% of the normalising distance at n " 16. These results demonstrate that the proposed scheme is capable of providing self-localisation with good accuracy and small error margins.

VIII. CONCLUSION
In this article, we propose a cooperative self-calibration scheme for sensor localisation in distribution fusion networks. The nodes in such networks collect measurements from a multi-object scene which involve uncertainties due to noise, probability of detection being less than one, and, the presence of false detections from the surroundings. Multi-sensor scalability is achieved by local filtering of these measurements, which in turn provides multi-object posteriors as the information entities to communicate and be fused for multi-sensor exploitation, as opposed to transmitting raw detections across the network.
We introduce a probabilistic model to estimate sensor calibration parameters without violating the locality structure of computations and by using nearest neighbour message exchanges. In particular, we develop node-wise separable parameter likelihoods which can be computed based on the locally available information for distributed fusion. The interactions between pairs of sensors sharing a communication link are used to build up a pairwise MRF model with these separable likelihoods as its edge potentials. Distributed estimation is, then, carried out using (loopy) BP over the communication graph.
Our approach inherets the capability of handling highly uncertain measurements from the local multi-object filters without posing any significant constraint on the system configuration. It can also be used in solving general calibration problems and at a fusion centre to avoid the computational burden of centralised multi-sensor filtering in parameter estimation. In the case of mobile sensor platforms, our framework can accommodate dynamical parameter models, as well.
There is a variety of possible extensions of this work. Different node wise separable parameter likelihoods can be found by using different strategies to approximate (15). More sophisticated sampling schemes can be adopted for the realisation of the algorithm. It is also possible to use more accurate multi-object models, in this framework, with the cost of increasing communication and/or computational demand.

A. Proof of Proposition 4.1
We begin the proof by substituting the distributions of concern in the definition of (conditional) KL divergence. For the sake of a shorter notation, we denote the sensor histories from time 1 to k´1, i.e., Z i 1:k´1 and Z j 1:k´1 by H i and H j , respectively. We also drop the time indices in the current measurements Z i k and Z j k . The KL divergence in (18) is given by D`ppZ i , Z j |H i , H j , θ i,j q||qpZ i , Z j |H i , H j , θ i,j q" ż δZ i δZ j δH i δH j dθ i,j ppZ i , Z j , H i , H j , θ i,j ql og ppZ i , Z j |H i , H j , θ i,j q qpZ i , Z j |H i , H j , θ i,j q " ż δZ i δZ j δH i δH j dθ i,j ppZ i , Z j , H i , H j , θ i,j ql ogˆp pZ i , Z j |H i , H j , θ i,j q ppZ i |H j , θ i,j qppZ j |H i , θ i,j q˙.
In the integration above, the variables regarding sensor measurements and histories are sets (as opposed to vectors), and, their integration is defined in Appendix C. As far as the proof is concerned, it suffices to have well-defined integration rules for the variables involved. The results hold true for any such variables including vector valued sensor measurements.
Next, we consider (39) and multiply the numerator and the denominator inside the log term by the product of ppZ i , Z j |H i , H j , θ i,j q, ppZ i |H i , H j , θ i,j q,ppZ j |H i , H j , θ i,j q, ppH i |H j , θ i,j q, and, ppH j |H i , θ i,j q. Thus, The second equality above follows from the definition of mutual information (MI) [34], and, in the penultimate line, the chain rule for information [34] is used with the positive terms.

B. Proof of Corollary 4.2
We use the simplified notation introduced in Appendix A for the proof of Proposition 4.1. The data processing inequality (DPI) [34] applied simultaneously to the first two terms in the RHS of (18) (or, the RHS of (40) in the simplified notation) in accordance with the (conditional) chains Z j Ø X k Ø H j |H i , θ i,j and Z i Ø X k Ø H i |H j , θ i,j , respectively, yields IpZ j ; H j |H i , θ i,j q`IpZ i ; H i |H j , θ i,j q ď IpX k ; H j |H i , θ i,j q`IpX k ; H i |H j , θ i,j q " HpX k |H i , θ i,j q´HpX k |H j , H i , θ i,j q HpX k |H j , θ i,j q´HpX k |H i , H j , θ i,j q. (41) For the third MI term in (18) (equivalently in (40)), there are two alternative ways to use the DPI in accordance with Z j k Ø X k Ø Z i k |H i , H j , θ i,j which hold true simultaneously: (42) and, IpZ j ; Z i |H j , H i , θ i,j q ď IpX k ; Z j |H j , H i , θ i,j q " HpX k |H j , H i , θ i,j q´HpX k |Z j , H j , H i , θ i,j q. (43) Therefore, combining (41), (42) and (43) leads to the upper bound for (18) given by (19).

C. Marginalisation of RFS variables
Throughout the article, RFS variables are often marginalised by integrating densities over all possible values. Let X be a RFS with elements in X . X takes values in the space of all finite sets denoted by F pX q. Consider a RFS probability density f : F pX q Ñ r0, 8q. The probability that X is contained in a closed subset S of X is given by the set integral [25,Chp.11] P rtX Ď Su " (44) Marginalisation of X involves integrating a probability density over F pX q which can equivalently be carried out through the set integral defined above: where µ is an appropriate measure, for example, the unnormalised distribution of a Poisson point process with a uniform rate (further details can be found in Section II.B in [40], and the references therein). This identity holds for any real valued measurable function f : F pX q Ñ R (Appendix B in [40]), thus, along with marginalisation of joint densities, computations regarding information measures and divergences of RFS distributions can be carried out using set integrals as well [46].