Distributed Space Debris Tracking with Consensus Labeled Random Finite Set Filtering †

Space debris tracking is a challenge for spacecraft operation because of the increasing number of both satellites and the amount of space debris. This paper investigates space debris tracking using marginalized δ-generalized labeled multi-Bernoulli filtering on a network of nodes consisting of a collection of sensors with different observation volumes. A consensus algorithm is used to achieve the global average by iterative regional averages. The sensor network can have unknown or time-varying topology. The proposed space debris tracking algorithm provides an efficient solution to the key challenges (e.g., detection uncertainty, data association uncertainty, clutter, etc.) for space situational awareness. The performance of the proposed algorithm is verified by simulation results.


Introduction
This paper presents a distributed space objects tracking approach in the context of space situational awareness (SSA) [1]. SSA is the ability to track and predict the velocity and location of space objects in orbit around the Earth. SSA has become a significant concern for both commercial and military systems. A U.S. Iridium 33 communications satellite was struck by a defunct Russian Cosmos 2251 military communications satellite in February 2009. Collisions like this generate space debris, which then increases the likelihood of further collisions [2]. Collisions may degrade the performance of a spacecraft or even fragment one if they involve enough energy [3]. Multi-object tracking algorithms play a significant role in space object tracking [4,5].
The catalog size which future space objects tracking systems should be able to handle is more than 100,000 resident space objects [6]. A possible approach to tracking space debris is by the use of a network with regional nodes. Multi-nodal systems which consist of multiple agents with sensing, processing, and communication capabilities are widely used in multiple object tracking systems. Information fusion methods used in multi-node systems fall into three categories: distributed, centralized, and hierarchical. The Joint Space Operations Center (JSpOC) currently employs the centralized information processing approach because it can provide the most accurate estimation [6]. As the space debris tracking system has more nodes, the computation load of centralized systems becomes very challenging. Distributed information fusion can be used as an alternative solution [7]. The main advantages of distributed systems include: the system is reliable (i.e., resilient to failures), framework, the software Turboprop [43] and UKF are used to approximate the transition density of space objects. Turboprop has a package of functions which can be used to calculate the trajectories of space debris.
The main contribution of this paper is the distributed fusion algorithm for sensors with different observation volumes. Most fusion methods assume that all sensors have the same observation volume or that targets stay in the combined observation volume all the time. When targets leave the observation volume, normal algorithms cannot keep the estimation of targets because there are no measurements available. In this paper, the detection probability outside of the observation volume is set to zero. The expectation maximization (EM) algorithm is used to approximate the densities across the observation volumes. Since there are no measurements for targets outside the observation volume, the estimation is essentially the prediction. Then, the fusion among sensors with different observation volumes can be performed in the same way as fusion among sensors with the same observation volume.
The paper is organized as follows. Background on labeled RFS and the network model are described in Section 2. The Mδ-GLMB filter recursion is presented in Section 3. The space debris dynamic model is provided in Section 4. Section 5 details the estimation for objects when they are outside the observation volume. Section 6 details the consensus with Mδ-GLMB filtering. Numerical results are presented in Section 7, and concluding remarks are given in Section 8.

Notation
In this paper, the Kronecker delta that takes arbitrary arguments is denoted by: The inclusion function is denoted by: The standard inner product notation is denoted by f , g f (x)g(x)dx, and the exponential notation h X Π x∈X h(x), where h is a real-valued function. The weighting operator for a given probability density function (PDF) p and a scalar α is defined as: [p(x)] α p α , 1 . (3)

Network Model
The nodal network can be denoted by a directed graph G = (M, A). M denotes the set of nodes and A = M × M denotes the connections among nodes. (i, j) ∈ A if node j can receive data from node i. M (j) {i ∈ M : (i, j) ∈ A} represents the set of neighbors for each node j ∈ M (including j itself).
Every node in the network has the same importance and performs the same actions: gathers measurements, carries out local computation, and exchanges information with neighbors. The fusion is performed among its own information and information from its neighbors.

Labeled RFS and Bayesian Multi-Object Filtering
An RFS is a finite set valued random variable. A unique label ∈ L = {α i : i ∈ N} is added to the dynamic state x ∈ X to incorporate the identity of objects, where N is the set of positive integers and the α i 's are distinct. Labels for objects are ordered pairs of integers = (k, i), where k is the time of object birth and i ∈ N is the index to distinguish objects born at the same time. L k = {k} × N is used to denote the label space for objects born at time k. Then, the new object born at time k has the state x ∈ X × L k . The label space can be constructed recursively by L 0:k = L 0:k−1 ∪ L k . The abbreviation L L 0:k is used for compactness.
Let L : X × L → L denote the projection L((x, )) = , then the function ∆(X) δ |X| (|L(X)|) is called the distinct label indicator, which means that X has the same cardinality as its labels L(X) = {L(x) : x ∈ X}. The labeled RFS can be unlabeled by discarding the labels. Let the labeled RFS distribution be denoted by π({(x 1 , 1 ), ..., (x n , n )}). Then, the unlabeled RFS is distributed according to: The cardinality distribution (the distribution of the number of objects) of a labeled RFS is the same as its unlabeled version. More information about labeled RFSs can be found in [28,29].
Assume that there are N(k) object states x k,1 , ..., x k,N(k) with state space X × L at time k, and M(k) measurements z k,1 , ..., z k,M(k) with observation space Z. Then, the set of objects and observations are treated as the multi-object state and multi-object observation, respectively: Let π k|k−1 denote the multi-object prediction density and π k (·|Z k ) denote the multi-object filtering density at time k. Then, the object density is propagated by the multi-object Bayes filter in time according to: where f k|k−1 (·|·) is the multi-object transition density to time k and g k (·|·) is the multi-object likelihood function at time k. The underlying models of object births, deaths, and motions are encapsulated in the multi-object transition density, while the underlying models for detections and false alarms are encapsulated in the multi-object likelihood function. For convenience, we omit the references to the time index k and denote g g k , f f k|k−1 , π + π k|k−1 , π π k , B L k , L + L ∪ B.

δ-Generalized Labeled Multi-Bernoulli
The δ-GLMB density has the form where each I represents a set of track labels, each ξ represents a history of association maps, and Ξ is a discrete space. The pair (I, ξ) ∈ F (L) × Ξ is called a hypothesis and the associated weight ω (I,ξ) means the probability of the hypothesis. p (ξ) is the density of the kinematic state [28,29].

The Mδ-GLMB Filter Recursion
Mδ-GLMB filter is an approximation to the δ-GLMB filter. It can be interpreted as performing a marginalization over the association histories. The Mδ-GLMB filter has the same cardinality distribution and PHD as the δ-GLMB filter [31].
An Mδ-GLMB density corresponding to the δ-GLMB density in (6) is of the form where A Mδ-GLMB is a special case of the δ-GLMB, and it is completely characterized by the parameter set π = {(ω (I) , p (I) ) : I ∈ F (L)}. The Mδ-GLMB prediction and update are outlined below.

Mδ-GLMB Prediction
For the current time and a given multi-object state X, each state (x, ) ∈ X either survives with probability p S (x, ) and evolves to the next step with new state (x + , + ), or dies with probability 1 − p S (x, ). The birth density of the new objects at the next time step is where p B and ω B are given parameters of the birth density. The multi-object state at the next time step X + includes two parts: the new-born objects and the surviving objects. Supposing that the births are independent of surviving objects and that objects evolve independently of each other, the multi-object transition density is given by [28,29]: where Let k ) : I ∈ F (L)} denote the Mδ-GLMB multi-object posterior density at time k. The multi-object prediction density is the Mδ-GLMB where ω (I)

Mδ-GLMB Update
Each state (x, ) ∈ X has the probability p D (x, ) to be detected, generates a measurement with likelihood function g(z|x, ), and has a probability 1 − p D (x, ) to be mis-detected. The measurement set Z is a superposition of detected objects and Poisson clutter with intensity function κ.
An association map is denoted by θ: L → {0, 1, ..., |Z|} such that θ(i) = θ(i )> 0 implies i = i . The set Θ of all such association maps is used to denote the association space, and the subset of association maps with domain I is denoted by Θ(I).
The ranked assignment and K-shortest paths algorithms are used in [29] to truncate the multi-object posterior and prediction densities, respectively. The algorithm has a cubic complexity in the number of measurements.
Given the posterior Mδ-GLMB density π = {(ω (I) , p (I) ) : I ∈ F (L)}, a tractable suboptimal multi-object estimation can be derived as follows: (1) Determine the maximum a posteriori cardinality estimation N * from (2) Among the label set with cardinality N * , determine the label set I * with highest weight ω(I): (3) Finally, determine the expected values of the kinematic states X * from p(x, ; I * ):

Space Debris Dynamic Model
The tracking algorithm needs a dynamic model of the debris and an observation model to track objects. The dynamic model is represented by the transition density function f k|k−1 (X k |X k−1 ). It is very difficult to calculate the Markov transition density function of space objects because the space dynamic model is much more complicated than a constant velocity or a constant turn model. However, we can approximate the transition density function of space debris with the help of the software Turboprop and the unscented transform [44]. Turboprop has a package of functions to calculate the trajectories of space objects, and can be called as a function in MATLAB. The elements in Turboprop include: an Earth orientation and atmospheric drag model; the lunar gravity models LP100K, GLGM-2, and LP150Q; the Earth gravity models JGM-3, GGM02C and WGS-84; JPL planetary ephemerides DE403 and DE405; and a solar radiation pressure model.
The Jet Propulsion Laboratory Development Ephemeris (JPL DE) are generally created to support spacecraft missions to the planets. JPL DE designates a series of models consisting of representations of accelerations, velocities, and positions of major Solar System bodies. The acceleration caused by the solar radiation pressure isr where p SR is the pressure of solar radiation in Pa. c R is the solar radiation coefficient and taken as 1.5 m is the mass of the space debris in kilograms. A is the cross-sectional area of the space debris facing the Sun in square meters. m and A are 0.05 kg and 0.01 m 2 , respectively.r represents the vector from the center of the Sun to the space debris. The mass concentration model is a gravity field, with the total acceleration determined by a series of point masses. This paper takes the Sun, the Earth, Venus, Jupiter, and the Moon into consideration. The unscented transform is used in the nonlinear projection of mean and covariance estimations. The unscented transform approximates the probability density function (PDF) by a bunch of sigma points. Suppose that each single object density p(x) is a Gaussian mixture of the form ∑ N i=1 ω i N (x; m i , P i ) and the Gaussian item is n-dimensional. Then, the state of space debris is propagated as shown in Table 1.  The sigma-points are chosen as follows: , j = 1, ..., n , where n represents the dimension of the state vector (6 in this paper), ∑ 2n j=0 W (j) = 1, (n + κ)P j is the jth row or column of the square root (n + κ)P, κ is the scaling parameter (2 in this paper).

Estimation of Objects Outside the Observation Volume
Most of the literature in object tracking assumes that objects do not go outside of the observation volume. This is not the case for space debris tracking. Filters with no special consideration of detection probability setting cannot provide estimation for objects outside the observation volume because there are no measurements available. However, since the observation model only affects the update and not the prediction, the estimation of objects can still be provided because the prediction is always available. Aiming to provide estimation for the whole area, we set the detection probability to zero for areas outside of the observation volume.
We assume that the single object density is represented by a Gaussian mixture (GM). Since UKF is used in this paper for the space debris transition function, each Gaussian item is represented by a bunch of sigma-points. If all sigma-points are inside or outside the observation volume, the detection probability for this Gaussian item is set to p Din and p Dout , respectively. p Din is the detection probability inside the observation volume and p Dout is the detection probability outside the observation volume. Otherwise, the EM algorithm is used to build a new GM to approximate the current Gaussian item.
An illustration is shown in Figure 1. The predicted probability density function is shown by a GM. The blue circle is one Gaussian item. V1 represents the part inside the observation volume, and V2 represents the outside part. * represents the sigma points obtained by the unscented transform. p Din = 0.98 is an example of the detection probability inside the observation volume. Since there is no measurement from outside the observation volume, the detection probability is set to p Dout = 0. The EM algorithm is used to build a new GM to represent V1 and V2, respectively.
The EM algorithm is an efficient method to find the maximum-likelihood estimate of the unknown parameters of an underlying distribution from a given data set. Usually, the data set has missing values or is incomplete. Each iteration of the EM algorithm has two procedures: The expectation step (E-step) and the maximization step (M-step). The E-step creates a function for the expectation of the log-likelihood evaluated using the current estimate for the parameters. In the M-step, the likelihood function is maximized under the assumption that the missing data are known.
If the sigma-points of one Gaussian item are on both sides of the observation edge, then N = 1000 points are sampled from this Gaussian distribution. Our experience shows that 1000 points are enough to represent the probability distribution in the experiment. Depending on whether the points are inside the observation volume or not, these points are divided into two groups: Points in and Points out . Points in are the sampled points inside the observation volume and Points out are the sampled points outside the observation volume. The EM algorithm is used afterwards for a Gaussian mixture parameter estimation. Each Gaussian item is represented by two new Gaussian mixtures. That is, The first GM is to represent V1 and the second GM is to represent V2. An illustration of the method is shown in Figure 2. The estimation of objects is available when they are outside the observation volume. The proposed consensus Mδ-GLMB can be used for sensors with different observation volumes.
The parameters in the density approximation are the number of Gaussian items, the weights for the items, the covariance matrices, and the means. These parameters are chosen as follows.
(1) The number of Gaussian components: The number of new Gaussian mixture components is important in the density approximation with the EM algorithm. However, it is very difficult to derive the number of Gaussian components. In practice, we choose a reasonably large number (J max = 10 in this paper) as the maximum number of new Gaussian mixture items in the approximation. The EM algorithm is used for situations with a different number of Gaussian components (J min = 1, ..., J max = 10). Then, the Bayesian information criterion (BIC) is used to determine the best number of Gaussian components. The BIC is a criterion for model selection among a finite set of models [45]. The model with the lowest BIC is preferred. The BIC is defined as follows: where J is the number Gaussian items set for the approximation, θ is the estimated parameter values that maximize the likelihood function, P is the sampled data, L(·) is the maximized value of the likelihood function, M is the number of parameters estimated (6 in this paper), and N is the number of sampled points (1000 in this paper).
(2) The weights for the Gaussian items: The weights can be the same and set as W j = 1 J .
(3) The covariance and the means: The means of the Gaussian position components can be set by uniformly sampling in the observation volume. The means of the Gaussian velocity components can be set to zero. The initial covariance is chosen by experience, which is the the same covariance matrix of the Gaussian item that was sampled from.

Distributed Information Fusion with Consensus for Mδ-GLMB
After each node finishes the state estimation of space debris, a distributed information fusion can be performed to achieve better results. The consensus algorithm is used in this paper to achieve distributed information fusion. This section shows that the Mδ-GLMB is algebraically closed under Kullback-Leibler averaging (KLA) (i.e., the KLA of Mδ-GLMB is also Mδ-GLMB) [18]. The closed form expression for KLA of the Mδ-GLMB is derived and used for consensus fusion of the Mδ-GLMB posterior density [18].

Consensus for Mδ-GLMB Filtering
The assumption that the estimation from each node in a multi-sensor network is independent or that the correlations among them are known is wrong. The correlation is commonly caused by data incest. Data incest happens when the same information takes several paths from the other sensors to the fusion center and the raw measurements are inadvertently used multiple times. Centralized information fusion is carried out by the Bayes solution: The result can be carried out as follows if the information is treated as being independent: Since there is a common measurement history and common process noise, the information from each node are not independent from each other. Data incest effects can be effectively counteracted by the consensus algorithm.
Given labeled multi-object densities π (i) on F (X × L), i ∈ I, and the normalized non-negative weights w (i) , i ∈ I, the weighted KLA is defined bȳ where is the Kullback-Leibler Divergence (KLD) of π (i) from π [20,22]. The set integral is used in the integral for any function f on F (X × L): Given normalized non-negative weights w (i) , i ∈ I and labeled multi-object densities π (i) , i ∈ I, then the weighted KLA in (35) is the normalized weighted geometric mean [18]. It was shown in [13] that independent, identically distributed cluster and Poisson RFSs are algebraically closed under KL averaging. The Mδ-GLMB family is also algebraically closed under KL averaging.
The following result shows that the KLA of the Mδ-GLMB densities is also a Mδ-GLMB density [18].
Given the Mδ-GLMB densities π (i) = {(ω (I) i , p (I) i ) : I ∈ F (L)}, i ∈ I and normalized non-negative weights w (i) , i ∈ I, the KLA, and hence the normalized weighted geometric mean, is an Mδ-GLMB given by [18] whereω The component (ω (L) ,p (L) (·)) of the KLA Mδ-GLMB can be rewritten as where (44) is the Chernoff fusion rule for the single-object PDFs. It can be seen from (43) and (44) that each fused Mδ-GLMB component (ω (L) ,p (L) (·)) can be independently calculated, which makes the fusion procedure fully parallelizable. Given a node network N with multi-object probability densities π (i) from each node i, and normalized non-negative weights w (i,j) for node i to nodes j ∈ N (i) , with ∑ j∈N (i) w (i,j) = 1, the KLA over the whole network can be calculated in a scalable and distributed way using the consensus algorithm. Suppose that each node starts with initial PDF π (i) 0 = π. The nth consensus can be calculated by where w (i,j) n is the (i, j)-th entry of the matrix Ω n . The (i, j)-th entry of the consensus matrix Ω is given by w (i,j) 1 N (i) (j). The consensus algorithm enjoys some good convergence properties. It was shown that if the consensus matrix is doubly stochastic (all rows and columns sum up to 1) and primitive (there exists an integer m that Ω m has all positive entries), then So, with a doubly stochastic and primitive consensus matrix, the global unweighted KLA of the densities over the whole network can be calculated by the consensus iterative of each node as the number of consensus steps tends to infinity [13].
A necessary condition for the consensus matrix to be primitive is that the network is strongly connected, which means that for any pair of nodes (i, j) ∈ N there exists a directed path from node i to node j and vice versa. This condition can also be satisfied if w (i,j) > 0 ∀i ∈ N , j ∈ N (i) . The consensus matrix is primitive and doubly stochastic for a undirected network (i.e., whenever node i sends information to node j, it also receives information from node j) if [46]: The global unweighted KLA of the multi-object posterior densities can be achieved by iterative local KLA averaging as the consensus step tends to infinity. In practice, the iteration is stopped at some finite number.
A Gaussian mixture (GM) is a typical choice to represent single-object densities. The fusion rule involves multiplication and exponentiation of GMs which generally do not provide a GM. A suitable approximation of the GM exponentiation has to be devised in order to preserve the GM form. For the sake of simplicity, let us consider only two agents (labeled a and b) with Gaussian mixture location densities The fused location PDF is also a Gaussian mixture, and can be approximated as follows: where The validity of the equations above depends on the cross-products of the different terms in the Gaussian mixture being negligible, which requires that the Gaussian components are well separated. This condition is well-met since a suitable merging step is used to fuse Gaussian components with Mahalanobis distance [47] below a given threshold. The fusion can be easily extended to more nodes by sequentially applying the pairwise fusion (50) and (51).
The other common way to represent a single-object density is by particles. The local filtering steps are more resource-demanding than a GM implementation. Moreover, the in-node computation burden is increased compared to GM representation because the information fusion requires additional techniques (e.g., least square estimation, kernel density estimation, or parametric model approaches).

Discussions on Consensus
Each node in the network performs local prediction and update. They then exchange information with their neighbors. Information fusion with the consensus algorithm is carried out afterwards. Depending on the requirements of the system, the procedure can be repeated N times, which is called N-steps of consensus. An illustration of the system design is shown in Figure 3.

Node Cluster
Perform GLMB filtering in parallel N1 N2 N3  The network used in this paper is a fully distributed one. Every node performs the same operations: updates its own information, exchanges information with its neighbors, and performs information fusion using the consensus algorithm.
The consensus algorithm achieves global averaging over the whole network with iterative averaging among neighboring nodes for each node. The main disadvantage is the computation complexity. With an increase of network size, the computation load for the whole network and the time it takes to reach global averaging significantly increases. The performance degradation is mainly due to the time it takes for the information from each node to be transferred to the rest of the nodes. In practice, the number of nodes in the whole network is chosen according to the computation power of each node so that the time required for the consensus algorithm is kept under a certain threshold.
It was proved in [48] that the complexity of the centralized Mδ-GLMB filter is linearly related to the number of sensors. Each node in the consensus Mδ-GLMB filtering system has to carry out local prediction and local update. The computation complexity is independent of the number of nodes. The computation complexity of the consensus algorithm for each node only depends on the number of neighboring nodes, and has nothing to do with the total number of nodes in the whole network.
As the number of targets increases, tracking becomes more challenging. Even though the ranked assignment algorithm and K-shortest paths algorithm are used to truncate the posterior densities and prediction densities, respectively, the filtering still has cubic complexity in the number of measurements. An efficient implementation of the GLMB filter by combining the prediction and update into a single step was proposed in [30]. The earlier implementation involves separate truncations in the prediction and update steps. The proposed implementation requires only one truncation procedure for each iteration. Furthermore, an efficient algorithm based on Gibbs sampling for GLMB filtering density truncation was also proposed in [30]. The implementation has a linear complexity in the number of measurements and a quadratic one in the number of objects. The Mδ-GLMB filter used in this paper was proposed in [31]. It is a tractable multi-object density approximation that can capture statistical dependence between objects. It matches the cardinality distribution and the first moment of the labeled multi-object distribution of interest. The performance of Mδ-GLMB filtering is sufficient for the situation in this paper.
The computation load of fusion for each node only depends on the consensus steps and the number of neighbors. For example, if the diameter of a network is 3, then it takes three steps at most for information in any node to be transferred to any other node in this network. As the number of nodes increases, it takes more time to reach the global average. However, the computation for each node at one time is the same, and only depends on the number of neighbors. Large networks will only be limited by the computation capability of each node. Therefore, the system is scalable and the processing load for each node is scalable with respect to the size of the network.
Fully decentralized systems suffer from the problem of synchronization and inconsistency. As the number of nodes increases, it becomes more difficult to achieve synchronization in the system. The effect of communication delays becomes more severe in larger networks. Even if the time delays are small, they can deteriorate the system's performance or even destabilize it. The output synchronization of nonlinear systems with communication time delays was discussed in [49]. A new framework was proposed in [50] to address the consensus with multi-agent systems and the synchronization of complex networks. Coordination can be achieved for the discrete-time delayed systems with linear dynamics [51] and switching topologies [52].
A space debris tracking approach with GM-CPHD and consensus algorithm was presented in [42]. The performance of the approach was shown with an example of 45 clusters of sensors to achieve global tracking of space debris. While the performance was very good, the approach used in this paper has much better results. This is because the Mδ-GLMB filter has a more accurate propagation of the posterior density than the CPHD filter, which makes the Mδ-GLMB filter better able to locate targets. Another disadvantage of the CPHD filter is the "spooky" effect which causes the CPHD filter to temporarily drop tracks which are subjected to missed detections and to declare multiple estimates for existing tracks in place of the dropped tracks [29]. Further, the CPHD filter cannot provide the trajectories of the targets. The estimates of the target from the CPHD filter are indistinguishable. The Mδ-GLMB filter used in this paper significantly outperforms the CPHD filter used in [42]. This paper presents a scalable solution to the problem of tracking space debris. The space tracking approach employed in this paper is an improvement over existing methods in several aspects. Firstly, the core elements of data association, detection, and tracking are solved with an integrated approach based on labeled RFS in this paper, while most SSA literature treats them separately. Secondly, most SSA methods use a single sensor (e.g., [4,35,41,53]) or a centralized tracking system as in [6] to track space objects. The consensus algorithm used in this paper has better performance than single-sensor methods and has a smaller computation load than a centralized tracking system. Thirdly, the consensus approach used in this paper can provide estimation for targets outside of the observation volume and can be used for information fusion for sensors with different observation volumes, unlike most consensus algorithms (e.g., [13,18,54]) or other distributed information fusion algorithms (e.g., [55,56]) in the current literature which can only be used for sensors with the same observation volume.
The consensus iterate of each node converges to the global unweighted KLA of the multi-object posterior densities as n tends to infinity. In practice, consensus steps N is chosen based on the network size. A methodology to achieve distributed average consensus in finite time was proposed in [57]. The proposed algorithm requires less memory at the cost of a slight increase in the number of steps required for termination. To relax the assumption that the nodes are aware of the upper bound on the network diameter, the authors provide an upper bound on the network diameter which, in the worst case, is twice the actual diameter. Ref. [58] reaches consensus in finite time using only linear iterations. The authors show that finite-time average consensus can always be achieved for undirected networks.
This paper treats the consensus step N as prior information. N is broadcasted to every node before the information fusion. This is certainly possible, but it becomes difficult when the size of the network grows. There are some approaches that do not require continuous data exchange nor knowledge on a global parameter. Ref. [59] reviews some of the main discrete and finite time average consensus implementations, from theoretical and empirical points of view. Three main aspects, namely computational analysis, packet loss resilience analysis, and stealth attacks resilience analysis, are analyzed. The authors of [59] first review synchronous and uniformly sampled approaches which are based on the standard average consensus iteration [60], on the flooding approach [57], or on a repetition of the max-consensus procedure [57]. They then present event-based finite-time average consensus implementations with point-to-point communication capability based on the Converge Cast [61] and on the token-passing approaches [62].

Numerical Results
The goal of this paper is to show the efficacy of consensus Mδ-GLMB in tracking space debris. We demonstrate this through three experiments. The first experiment was designed to show that the Mδ-GLMB filter proposed in this paper with special consideration of detection probability setting can provide estimation for objects outside the observation volume. The second experiment was designed to show that consensus fusion can be performed among sensors with different observation volumes. The third experiment was designed to show that the consensus algorithm can significantly improve the tracking performance when all objects are inside the observation volumes of all sensors.
The object state is a vector of position and velocity, x k = [p x p y p z v x v y v z ] T . The state dynamic model is: The first step is to build the space debris trajectories with the help of Turboprop. f (·) is the transition density function and w k ∼ N(·; 0, Q) is the process noise with Q = diag([σ 2 px , σ 2 py , σ 2 pz , σ 2 vx , σ 2 vy , σ 2 vz ]), σ px =σ py =σ pz =1 km, σ vx =σ vy =σ vz =0.01 km/s. The birth model to generate the simulation scenario is a labeled multi-Bernoulli RFS with parameters π B = {r  B is the birth intensity, which is the expected number of new objects born at each time step, at each object birth place i. All targets are initially born inside the observation volume of a sensor (not necessarily the same one). The situation where a target is born outside of the observation volume of all sensors and goes from outside to inside of the observation volume can be handled by a measurement-based birth model, which is beyond the scope of this paper. The space object trajectories generated by the birth model in this experiment are very common in real space object tracking scenarios.
Clutter is modeled as a Poisson RFS with intensity κ k (z) = λ c Vu(z). λ c is the average of clutter returns per unit volume, which was 100 in the simulation. u(·) is the uniform density over the observation region, V is the "volume" of the observation region.
All nodes used in this paper were located on the surface of the Earth and rotated with the Earth. The latitudes and longitudes of nodes were respectively. Node1 was initially on the x-z plane.
The sensor's observation is a noisy bearing and range vector given by where ε k ∼ N(·; 0, R k ) with ∆ = 40 s was the sampling period. The detection probability inside the observation volume p Din = 0.98. The δ-GLMB filter was capped to 10,000 components. In order to show statistical significance, results were shown over 100 Monte Carlo trials for the same object trajectories but different, independently generated, clutter and measurement noise realizations. The first experiment showed the efficacy of the algorithm with one realization. The second and third experiments showed with the mean value over 1000 Monte Carlo trials. Typically, the optimal sub-pattern assignment (OSPA) performance in the RFS area is shown with mean values. The OSPA distance is the sum of two parts: location error and cardinality error. The cardinality error is shown with not just the mean value, but also the standard deviation.
Performance evaluation of the multi-target algorithm is of great practical importance in the design and comparison of tracking systems. The difference between the estimation and true trajectories is very small compared with the motion range of space debris. Therefore, the OSPA metric was used to evaluate the performance with Euclidean distance p = 2, and cutoff parameter c = 20 km for the second experiment and c = 1 km for the third experiment. The OSPA is a mathematically and intuitively consistent metric for multi-object filtering performance evaluation [63]. OSPA distance has two components: localization and cardinality. The total OSPA distance is the sum of the two components. The cut-off parameter c determines the penalty weighting for cardinality errors as opposed to localization errors. The second experiment was mainly designed to show the estimation performance for targets when they go outside the observation volume. As there were no measurements available, the estimation error became larger and larger, and so did the OSPA value. Therefore, c = 20 km was used in this experiment. The third experiment was mainly intended to show the efficacy of the consensus algorithm when targets were all inside the combined observation volumes.
Since measurements were available all the time, the estimation error was smaller than in the second experiment. So, the cut-off value c = 2 km was used in the third experiment.
The consensus algorithm stops at N steps. The information of N should be broadcasted to every node beforehand and treated as prior information. In practice, N is chosen based on the size of the network. This paper shows the estimation performance with 1-step consensus, 2-step consensus, and 3-step consensus.

Estimation of Objects Outside the Observation Volume
This section demonstrates that estimation can be provided for objects outside the observation volume by setting the detection probability outside the observation volume to zero. Only node1 is used in this section. A fixed number of objects are presented here for simplicity. All objects are inside the observation volume initially, but as time goes on, some objects go outside. The field of view for node1 is This yields a sensor volume of 150 deg 2 . Some objects go outside of the observation volume with this setting.
The estimation from the Mδ-GLMB filter with uniform detection probability [13,18,48,54] is shown in Figure 4a. The filtering with uniform detection probability cannot provide estimation for targets when they are outside of the observation volumes. The estimation from Mδ-GLMB with detection probability outside of the observation volume set to zero is shown in Figure 4b.
Since this section is intended to show the estimation of target positions, estimation results from one simulation are shown here, instead of the mean values from all simulations. Note that all the simulations had similar results. It can be seen from Figure 4b that even when the measurements for objects outside the observation volume were not available, estimation could still be provided. In this case, the estimation for objects outside was essentially the prediction. The figure also demonstrates that the estimation became worse for objects outside the observation volume since there were no measurements to update.

Consensus among Nodes with Different Observation Volumes
This subsection demonstrates the effectiveness of information fusion with a consensus algorithm for nodes with different observation volumes. It is complicated to fuse information from nodes with no estimation for objects outside the observation volumes. The fusion algorithm must consider the observation volume overlap, which makes the fusion more difficult when more nodes are involved. With the estimation available for objects outside the observation volume, the fusion algorithm can be performed in the same way as for nodes with the same observation volume. This ensures that all objects were inside the combined observation volume of all sensors. The OSPA performance of a single-node and 5-node consensus is shown in Figure 5. The cardinality performance of a single-node and 5-node consensus is shown in Figure 6.  It can be seen that the single-node Mδ-GLMB filter [48] could provide good estimation for objects inside the observation volume. Since there was no update for objects outside the observation volume and only prediction was available, the estimation became worse with time. Two objects disappeared at the 140th time step, at which time they were outside the observation volume. A single node could not provide accurate cardinality estimation for objects outside the observation volume. Since all objects were in the combined observation volume of five nodes all the time, information fusion from five nodes could provide better estimation than a single node for objects outside the observation volume.
Since there was no measurement outside the observation volume, the best we could do is use the prediction model as the final estimation. The estimation performance was subject to the dynamic model which was characterized by the transition density noise. As time progressed, the estimation error became larger and larger without bound until new measurements were captured. The estimation confidence can be illustrated by the covariance matrix.

Consensus for Nodes with Similar Observation Volumes
This subsection is to show the performance of the consensus algorithm for nodes with similar observation volumes. The angle range of all nodes was 50 • , which was large enough for all objects to stay in the observation volume of all nodes all the time. The connection among the nodes is shown in Figure 7. A network with only five nodes was used in this paper to show the efficacy of the proposed method in distributed space debris tracking. A larger-scale network would involve more computation and processing time. We only show the estimation from node3. The numerical experimental results from a single sensor [48], 1-step consensus, 5-node consensus, and centralized Mδ-GLMB filtering [48] are shown in Figure 8. Five-node consensus was the fusion procedure for all the information from five nodes. It is shown that the performance of 1-step consensus was better than single-node estimation. Since 5-node consensus uses information from all the nodes, it had better performance than 1-step consensus. Centralized Mδ-GLMB filtering makes use of raw measurements from all sensors, and had the best performance. OSPA performance with different consensus steps is shown in Figure 9. It is shown in the figure that the OSPA performance was better with more consensus steps. With the designed network topology, the 3-step consensus had similar performance to the 5-node consensus. This is because it took 3 steps at most to transfer information from all of the other nodes to node3. After 3 steps of consensus, every node in the network had similar information, which is the average of the information from all nodes. There was no central node and the processing for each node was only carried among its neighbors, which makes the fusion scalable with respect to the size of the network.  OSPA performance with different consensus steps. More consensus steps yielded better performance.
The mean and standard deviation of the estimated cardinality of single-node Mδ-GLMB, 1-step consensus, 2-step consensus, 3-step consensus, 5-node consensus, and centralized Mδ-GLMB filtering versus time are shown in Figures 10 and 11, respectively. All filters estimated the cardinality accurately. The centralized Mδ-GLMB filtering had the best cardinality estimation performance. Consensus with more steps showed better cardinality estimation.    Figure 11. The standard deviation of the estimated cardinality of single-node Mδ-GLMB, 1-step consensus, 2-step consensus, 3-step consensus, 5-node consensus, and centralized Mδ-GLMB filtering versus time.
We can see that there are peaks in the plots in Figures 8, 9 and 11. The peaks happen at the time when there was a change in the number of objects. The change in the number of objects also contributed to the peaks in the cardinality standard deviation estimation. Cardinality estimation had a larger error than when the number of objects was not changing. The peaks in the total OSPA distance are because of the peaks in the cardinality OSPA distance.
The network can also be time-varying because the processing is performed only among neighbors. The time-varying topology is shown in Figure 12, and the performance of the consensus algorithm is shown in Figure 13. Node3 and node4 were in the network from the start to the 60th time step. Node3, node4, and node1 were in the network from the 61st time step to the 130th time step. All five nodes were in the network the rest of the time. Carrying out consensus only among a node's neighbors makes it possible for the method to handle a network topology that is time-varying.

Concluding Remarks
This paper presented a consensus algorithm for space debris tracking with labeled RFS. The key innovation lies in the fusion among sensors with different observation volumes, and most importantly, the processing load for each node is scalable with respect to the size of the network. The system is robust because of the avoidance of data incest. A salient feature is that the network topology can be time-varying. The scalability makes the approach appealing for future space object tracking systems as more nodes are connected to the network. Experiments confirm that the proposed algorithm has the potential to be used in distributed space debris tracking systems.
In recent years, there is more need for finite-time average consensus algorithms. Future work will focus on distributed space debris tracking systems with finite-time average consensus methods. Strengths and shortcomings of different consensus algorithms will be presented for different tracking scenarios in future work.