Sensor Selection for Decentralized Large-Scale Multi-Target Tracking Network

A new optimization algorithm of sensor selection is proposed in this paper for decentralized large-scale multi-target tracking (MTT) network within a labeled random finite set (RFS) framework. The method is performed based on a marginalized δ-generalized labeled multi-Bernoulli RFS. The rule of weighted Kullback-Leibler average (KLA) is used to fuse local multi-target densities. A new metric, named as the label assignment (LA) metric, is proposed to measure the distance for two labeled sets. The lower bound of LA metric based mean square error between the labeled multi-target state set and its estimate is taken as the optimized objective function of sensor selection. The proposed bound is obtained by the information inequality to RFS measurement. Then, we present the sequential Monte Carlo and Gaussian mixture implementations for the bound. Another advantage of the bound is that it provides a basis for setting the weights of KLA. The coordinate descent method is proposed to compromise the computational cost of sensor selection and the accuracy of MTT. Simulations verify the effectiveness of our method under different signal-to- noise ratio scenarios.


Labeled RFS and Mδ-GLMB
In this paper, the unlabeled and labeled variables are, respectively, represented by the italics and bold. For example, the unlabeled state, measurement and their sets are noted as x, z, X and Z; the labeled state and its set are noted as x = (x, ) and X, where is the discrete label of x. Let L(X), |X| and X × L denote the label set, cardinality and space of X, where X and L are the spaces of the unlabeled state and label.
The state estimates of single-target and multi-targets derived from a measurement set Z are both the functions of Z. To make this clearer, they are, respectively, noted asx(Z) andX(Z).x(Z) andX(Z) are their labeled versions.
Let δ Y (X), 1 Y (X) and p X denote the functions of generalized Kronecker, inclusion indicator and multi-object exponential, For any real-valued function b(X) of X, its set integral b(X)δX is defined as b(X)δX = ∞ ∑ n=0 1 n! ∑ 1:n ∈L n X n b(X n )dx 1:n (4) where x 1:n = x 1 , . . . , x n and 1:n = 1 , . . . , n , X n = {x 1:n } is a n-element labeled set, X n and L n are the spaces of X n and 1:n . If X is a Mδ-GLMB RFS, then its density is described as [22] π(X) = ∆(X) ∑ I∈F (L) δ I (L(X))ω I p X I (5) where ∆(X) = δ |X| (|L(X)|) is a distinct indicator for the labels of X, I ∈ F (L) is a label set in the collection F (L) of finite subsets of L, the weight ω I is the existing probability of the label set I, p I (x) is the density of x involved in I. The Mδ-GLMB density is abbreviated as π = {(ω I , p I )} I∈F (L) and its cardinality distribution is where F n (L) is the collection of n-element subsets of L.

Information Inequality to RFS Measurement
Letx(Z m ) be an unbiased estimate of x derived from an m-element measurement set Z m and f (x, Z m ) be a joint density over the space X 1 × Z m . Assuming that regularity conditions hold and ∂ 2 log f (x, Z m )/∂x i ∂x j exists, the information inequality to RFS measurement is [35] Z m X 1 f (x, Z m ) x l −x l (Z m ) 2 dxdz 1:m ≥ J −1 m l,l , l = 1, . . . , L where z 1:m = z 1 , . . . , z m , L is the dimension of x, x l andx l (Z m ) are the lth components of the vectors x andx(Z m ), J m is the L × L Fisher information matrix (FIM) given |Z|= m , (7) holds with equality if and only if f (x, Z m ) satisfies the distribution of exponential family.

A New Metric for Labeled RFS
It is well known that the optimal sub-pattern assignment (OSPA) metric [39][40][41] has been extensively used to measure the distance for two unlabeled sets. Although the OSPA metric could also measure differences in the set labels, it may not be very appropriate for the labeled RFS in some application scenarios. As a result, a new metric between two labeled sets X and Y of order 1 ≤ p ≤ ∞ with cut-off c > 0 is proposed as follows. where L(X) and L(Y) are the label sets of X and Y, x ∈ X and y ∈ Y are the unlabeled elements corresponding to the label , |L(X) ∩ L(Y)| and |L(X) ∪ L(Y)| − |L(X) ∩ L(Y)| respectively indicates the number of elements in X and Y which have the same and different labels, denotes the 2-norm of x and y cut off at c > 0.
See Appendix A for the proof that the d (c) p (X, Y) is a metric. The proposed metric, which is named as the LA metric by us, is a different metric than the OSPA. A physical meaning about the LA metric between two labeled sets X and Y is explained as follows. For all x ∈ X and y ∈ Y, if x has the same label as y, then x is paired with y and a 'location' error d (c) (x, y) between them is involved in the metric; otherwise, x is unpaired with y and a 'penalty' error c for them is involved in the metric.
The most significant difference between the OSPA metric and the LA metric is: In the OSPA metric, the elements of X are paired with the elements of Y depending on the optimal assignment distance of their unlabeled versions. In contrast, in the LA metric, the elements of X are paired with the elements of Y completely depending on their labels.
Take the MTT with the labeled RFS state for example. The calculation of the OSPA error may pair an estimate with a state due to the rule of optimal assignment even if they have different labels. Instead, the calculation of LA error prohibits this kind of pairing even if their unlabeled states are sufficiently close to each other. As a result, the LA metric involves not only the estimation error arising from the target number and individual states as the OSPA metric but also the additional estimation error arising from the labels. In this sense, the LA metric is more demanding than the OSPA metric for measuring the error between the labeled state set and its estimate.

Problem Formulation
To simplify the formulas, the time index is omitted and the subscript '+' is used to indicate the predicted density.
Multiple targets independently move in region A with random birth and death. The multi-target states are modeled as a labeled RFS X. The dynamic of a single state x = (x, ) ∈ X is described by the survival probability p s (x) and transition density f ( x|x )δ ( ). The dynamic of multi-target states is described by the transition density f ( X|X ). Here x = (x , ) and X are, respectively, the state and state set at the last time.
Targets are observed by the decentralized sensor network shown in Figure 1. The network is composed of sensor nodes (SNs) and local fusion centers (LFCs). Each SN receives measurements and communicates with its superior LFC. Each LFC receives the measurements, conducts data processing and storage, communicates with the other LFCs connected to it and manages its subordinate SNs.
The network structure is completely described by a topological graph with parameter {N, C, A}, where N and C are the label sets of SN and LFC, A ⊆ C × C is the set of directed connections between LFCs. If the LFC j can receive data of the LFC i, then (i, j) ∈ A. Let C j = {i ∈ C : (i, j) ∈ A} be the label set of the LFCs connected to the LFC j (including itself) and N j be the label set of the SNs belonging to the LFC j. Each SN only belongs to one LFC, which indicates ∩ j∈C N j = ∅ and ∪ j∈C N j = N.
The most significant difference between the decentralized network and the centralized or hierarchical network is that the former has no global fusion center connected to all SNs or all LFCs. The network structure remains unchanged and all measurements are synchronized during the monitoring period. The network structure is completely described by a topological graph with parameter { } , , where N and C are the label sets of SN and LFC, A C C ⊆ × is the set of directed connections between LFCs. If the LFC j can receive data of the LFC i , then ( ) be the label set of the LFCs connected to the LFC j (including itself) and j N be the label set of the SNs belonging to the LFC j . Each SN only belongs to one LFC, which indicates = j j C N The most significant difference between the decentralized network and the centralized or hierarchical network is that the former has no global fusion center connected to all SNs or all LFCs. The network structure remains unchanged and all measurements are synchronized during the monitoring period.
For the SN s N ∈ , it may receive clutter and target measurement or miss the detection. The measurements are modeled as a RFS s Z over the space s  . s Figure 1. Diagram for decentralized sensor network. and denote SN and LFC, solid lines with arrows denote directed connections between LFCs, the dotted lines denote connections between LFC and its subordinate SNs.
For the SN s ∈ N, it may receive clutter and target measurement or miss the detection. The measurements are modeled as a RFS Z s over the space Z s . z s ∈ Z s is a single measurement. Clutter is modeled as a Poisson RFS with intensity κ s (z s ), and is the clutter rate. The multi-target likelihood of the SN s is obtained from [9] as where where p s d (x) and g s ( z s |x) are the single-target detection probability and likelihood, Θ s is a collection of association mapping θ s : L(X) → {0, 1, . . . , |Z s |} . θ s ( ) > 0 or θ s ( ) = 0 indicates that the track ∈ L(X) generates a measurement z s θ s ( ) ∈ Z s or be missed. Each track at most generates one measurement and each measurement is at most generated by one track, which indicates that = if θ s ( ) = θ s ( ) > 0. The number of association hypotheses is [42] Due to some restrictions, only part of SNs involved in the sub-network of each LFC can be activated to observe targets at each scan. Assume that the multi-target likelihoods of all the SNs in the network are independent of each other given the labeled state set X. Algorithm 1 presents the steps of sensor selection and MTT for the LFC j ∈ C under a Bayesian framework. Note that the sequence of measurement sets up to the last time is omitted here for simplifying the formulas of Bayesian recursion. Algorithm 1. Sensor selection and MTT for the LFC j.

Prediction:
Calculate the current predicted density by π j + (X) = f ( X|X )π j (X )δX , where π j (X ) is the fused density at the last time; 2. SN selection: Select the SN subset S j ⊆ N j and receive a collection of their measurement sets Z S j ; 3. Update: Calculate the current posterior density by π j X|Z S j = ∏ s∈S j g s ( Z s |X)π j + (X) ∏ s ∈S j g s ( Z s |X)π j + (X)δX and then transmit the density to the LFCs connected to the LFC j; 4. Fusion: Receive the posterior densities from the LFC set C j and then calculate the fused density by the weighted KLA rule π j (X) = is the preset normalized weight;

State extraction:
Extract the current state estimate from π j (X) as the output. Go to Step 1.

Remark 1.
It can be seen from Steps 3 and 4 of Algorithm 1 that in this network, each LFC only communicates once with the LFCs connected to it per recursion. Because of this, the consensus iteration [43,44] of each LFC cannot be carried out since each LFC has to communicate with other LFCs more than once in the iterative process. Finally, the consensus fusion [23] over the entire decentralized sensor network cannot be achieved. As a result, the multi-target estimates outputted by different LFCs may be different. In fact, this character is consistent with most practical application systems.
The SN selection in Step 2 of Algorithm 1 is described as the following constrained optimization problems: where ϑ j S j ; π j + , γ j i S j ; π j + ≥ 0 and ν j k S j ; π j + = 0 are the objective function, inequality and equality constraints of S j given the predicted density π j + (X). The ultimate goal for a MTT network is to optimize the tracking precision. As is known to all, the MSE between target state and its estimate is currently the most widely-used evaluation indicator for tracking accuracy. Given the selected SN subset S j for the LFC j ∈ C, the MSE σ j S j 2 between X and its Bayesian estimateX j Z S j is where Z S j is the joint measurement space of the SN set S j , f X, Z S j is the joint density of X, Z S j , e X,X j Z S j denotes the error distance between X andX j Z S j . In this paper, e X,X j Z S j is measured by the 2nd-order LA metric d , which provides an online indication on the limit of MTT accuracy within the labeled RFS framework. Treating π j + (X) as a default condition, (15) is finally rewritten as

Derivation of LA Bound
In Section 4.1, Section 4.2, Section 4.3 and Appendix B, the superscript 'j' for the index of LFC is omitted. For example, σ j S j 2 is abbreviated as σ 2 S . In order to derive σ 2 S , it is assumed that A1: Multi-target Bayesian recursion is a Mδ-GLMB RFS [22]. As a result, the predicted density π + (X) and posterior density π X|Z S can be described as π + = {(ω I,+ , p I,+ )} I∈F (L) and π ·|Z S = ω I Z S , p I ·|Z S I∈F (L) . A2: Although the optimal estimate of X can be extracted from the fused density π(X) by using the joint or marginal multi-target estimator [9], both the methods are very difficult to be implemented. Alternatively, the target number is firstly estimated according to the maximum a posterior (MAP) criterion and then the individual states are estimated according to the unbiased criterion under the obtained target number. In fact, the suboptimal method is applied in almost all multi-target Bayesian filters.
Let Z S m S = Z s 1 m s 1 , . . . , Z s |S| m s |S| be the collection of measurement sets from the SN set S and where m s i is the number of measurements received by the SN s i . Let q X n , Z S m S be the joint density over the space (X × L) n × Z S m S . According to Bayesian formula, q X n , Z S m S is written as where Ω n,m S is a normalization factor, (19) shows that Ω n,m S / m S ! · n! is actually the probability P X = n, Z S = m S , where m S ! = m s 1 ! · · · m s |S| ! and Z S = m S denotes |Z s 1 | = m s 1 , . . . , Z s |S| = m s |S| . Substituting A1 and (12) into (19) as well as using Lemma 12 in [17], Ω n,m S is obtained as where θ S ∈ Θ S denotes θ s 1 ∈ Θ s 1 , . . . , θ s |S| ∈ Θ s |S| , . Then, (21) can be rewritten as (22) where p I, denotes the probability that only the SN subset Y ⊆ S receives the measurement from state x while the others miss the measurement. (22) shows that ϕ I θ S is independent of the association mapping θ S . From (14), (20) and (22), Ω n,m S is finally obtained as Since q X n , Z S m S is permutation invariant over x 1:n , its marginal density over any of x 1:n is the same and is obtained by Substituting (18) into (25) as well as using A1 and the identical equation δ n (|{ , 2:n }|) = Substituting (12) into (26) and then simplifying the result, we get where q I x, Z S m S ; θ S is the marginal density of q X, Z S m S over any of x 1:n given the label set I and association mapping θ S , where θ s − θ s ( ) denotes the remaining association mapping in θ s except for θ s ( ). where Z Ŝ n,m S = Z s 1 n,m s 1 × · · · × Z s |S| n,m s |S| is the joint measurement subspace of the SN set S where the target number is estimated asn; Z S 0,m S , Z S 1,m S , . . . , Z S ∞,m S is a partition of Z S m S ; P |X|= n |Z S m S is the posterior probability of |X|= n given Z S m S . From (6), P |X|= n |Z S m S is written as where ω I Z S m S is the existing probability of the label set I given Z S m S . According to the update step of Mδ-GLMB [22], ω I Z S m S is obtained as (29) and (33).
Let Ψn ,n,m S be the integral of q X n , Z S m S over the space (X × L) n × Z Ŝ n,m S . From (18), Ψn ,n,m S can be written as shows that Ω n,m S Ψn ,n,m S / m S ! · n! is actually the probability =n, X = n, Z S = m S . Substituting A1 and (12) into (34) as well as using Lemma 12 in [17], Ψn ,n,m S is obtained as where Similar with ϕ I θ S , φ I,n θ S can be rewritten as (38) shows that φ I,n θ S is independent of the association mapping θ S . From (14), (35) and (38), Ψn ,n,m S is finally obtained as Theorem 1: Given A1, A2 and the SN set S, the lower bound for the 2nd-order LA metric based MSE of (16) is where c is the cut-off of the LA metric, L is the dimension of x, Ω n,m S and Ψn ,n,m S are given in (24) and (39), is actually the possible ratio of the number of common labels to the number of all distinct labels for the multi-target states and their estimates given where the integral region Z Ŝ n,m S for measurement is given by (30), where q n x, , Z S m S is given in (27). See Appendix B for proof of Theorem 1.

Remark 2.
The number of estimated targets is assumed to be unknown in the derivation of the proposed bound. Only the MAP rule, rather than the specific (or exact) number of estimated targets, is required to obtain the bound. The symboln used for calculating the bound in Theorem 1 is just an index for all possible (or unknown) number of estimated targets. This is similar with the symbols , k, n and m S in Theorem 1, which are just the indices for the labels, the number of common labels in true targets and their estimates, the number of true targets and the number of sensor measurements, respectively.
Furthermore, the reason for imposing the MAP rule has been explained in A2. (30) has also shown that the measurement space Z S m S can be divided into Z S 0,m S , Z S 1,m S , . . . , Z S ∞,m S according to the MAP rule. It is very helpful for the proof of Theorem 1.

Remark 3.
In general, the maximum number of targets and measurements can be presented by prior knowledge. Moreover, the label space L 1 = L 1 ∪ B 1 , where L 1 and B 1 are the label spaces for the last time and new-born targets. In general, B 1 can also be preseted by prior knowledge. According to these presets, the sum of infinite terms in (40) becomes the sum of finite terms.

Remark 4.
Once the specific forms of p I,+ (x), p d (x) and g s ( z s |x) are given, ∂q n x, , Z S m S /∂x i and ∂ 2 q n x, , Z S m S /∂x i ∂x j in Φ n x, , Z S m S can be obtained from (27) and (28).

Remark 5.
The formulas of Ψn ,n,m S and Jn ,n,m S ( ) contain the integral over the measurement subspace Z Ŝ n,m S , which is calculated via MC integration [45]. To improve computational efficiency, the samples of MC integration are selected as predicted ideal measurement sets (PIMS) [31]. The calculation steps are shown in Algorithm 2.

SMC and GM Implementations for the Bound
In order to derive the SMC implementation of the bound in Theorem 1, it is assumed that A3: Substituting (44) into (22), (28), (29), (33) and (38), , θ S , β I,Z S m S θ S and φ I,n (S) are rewritten as Finally, the SMC forms of Ω n,m S , q n x, Z S m S , Ψn ,n,m S and Φ n x, , Z S m S are respectively obtained by substituting (45)-(49) into (24), (27), (39) and (43), where ∂q n x, , Z S m S /∂x i and In order to derive the GM implementation of the bound in Theorem 1, it is assumed that A4: Each p I,+ (x) involved in the predicted Mδ-GLMB density π + = {(ω I,+ , p I,+ )} I∈F (L) is described by the GM form of where N ·; µ i I,+ ( ), Σ i I,+ ( ) denotes the Gaussian density with mean µ i I,+ ( ) and covariance matrix Σ i I,+ ( ), υ i I,+ ( ) and G I,+ ( ) are the weights and number of GM terms. A5: The detection probability p s d (x) is independent of x and the likelihood function g s ( z s |x) is linear Gaussian, where H s and R s are the observation function and covariance matrix for measurement noise. From A5 and (23), we have Substituting (53) into (22), ϕ I (S) is rewritten as Substituting (51), (53) and (54) into (28) as well as using Lemma 2 in [11], q I x, Z S m S ; θ S is rewritten as Similarly, substituting (51), (53) and (54) into (29), (33) and (38) ; Finally, the GM forms of Ω n,m S , q n x, Z S m S , Ψn ,n,m S and Φ n x, , Z S m S are respectively obtained by substituting (56)-(61) into (24), (27), (39) and (43).
Obviously, they no longer contain integrals of state x and have analytic forms except for Ψn ,n,m S . Here ; θ Y , Σ i I ( ; Y) /∂x i ∂x j , which can be obtained by the following formulas: If the observation model is non-linear, that is, where h s (x) is the nonlinear observation function of state x. In this case, the extended Kalman (EK) or unscented Kalman (UK) filter can be used to calculate the mean and covariance matrix of each GM term [46,47].

Sub-Optimization Based on Coordinate Descent
The computational cost of this method is composed of three parts: sensor selection, Mδ-GLMB filtering and weighted KLA fusion. The latter two have been analyzed in Reference [22,23]. This paper only studies the computational cost of sensor selection and its approximate algorithm. When the number of SNs is large, it has a much more significant effect on the total amount of computation than the last two.
As shown in (15) and (17), sensor selection is actually a constrained combinatorial optimization problem. To find the optimal solution by the exhaustive search method, the objective function needs to be repeatedly calculated for C |S| |N| = |N|! /(|S|!(|N|−|S|)! ) times. Obviously, its computational cost increases with the SN number |N| in the form of combination explosion. In order to reduce the computational cost, some heuristic optimization algorithms, such as genetic algorithm [39] and so forth, is used to tackle this problem. However, the convergence speed of the heuristic algorithms will become rather show when the objective function is relatively complex. As a result, to further improve the computational speed, coordinate descent method [37] is proposed to find a sub-optimal solution of (17). Its computational cost increases with the SN number |N| in an approximate polynomial form.
Set a binary switch variable ς s ∈ {0, 1}, s = 1, . . . , |N|, for each SN. ς s = 1 indicates s ∈ S while ς s = 0 indicates s / ∈ S. The vector ς = ς 1 , . . . , ς |N| is composed of the switch variables of all SNs belonging to the same LFC. Clearly, the set S is completely determined by ς. Then, (17) is relaxed to an unconstrained optimization problem by the augmented objective function where > 0 is a barrier factor, ∑ l i=1 γ −1 i (ς) and ∑ m j=1 ν 2 j (ς) are inequality and equality penalty terms. Algorithm 3 presents the iteration steps for handling the relaxed problem by using the coordinate descent method.

Algorithm 3. Coordinate descent method.
Step 1: Set initial iteration number i = 0, initial SN switch vector ς (0) , initial barrier factor and its reduction coefficient 0 < C < 1; Step 2: From s = 1 to s =|N|, calculate ς s (i+1) = argmin (i) are treated as constants; Step 3: If ς (i+1) = ς (i) , then go to Step 4; Otherwise, set i = i + 1, go to Step 2; Step 4: If ς (i+1) = ς (0) , then output ς (i+1) as the solution of (17); Otherwise, set = C , ς (0) = ς (i+1) , i = 0 and then go to Step 2. In order to improve the probability to converge to the global optimum and accelerate the convergence speed for the coordinate descent method, the initial barrier factor and its reduction coefficient C can be appropriately selected by the methods in Reference [48], the initial SN switch vector is set as ς (0) = ς * , where ς * is the outputted switch vector at the last time.

A1 indicates that the posterior density of each LFC is a Mδ-GLMB form of π i X|Z
, i ∈ C. Then, for the LFC j ∈ C, given the LFC subset C j connected to it and the normalized nonnegative fusion weights α j,i (i ∈ C j ), its fused density obtained by the weighted KLA rule is still a Mδ-GLMB form of π j (X) = ω j L , p j L (x, ) L∈F (L) [23], where the weight α j,i reflects the effect of local posterior density π i X|Z S i on the fusion of the LFC j. The larger the weight is, the greater the impact it has on the KLA fusion. The bound in Theorem 1 reflects the optimal MTT accuracy that is potentially achieved by a LFC after sensor selection. The larger the proposed bound of a LFC is, the worse the precision limit that it can achieve is. Therefore, the normalized weight α j,i in the KLA fusion of the LFC j should be set inversely proportional to the proposed bound σ i which indicates that the larger the proposed bound of the LFC is, the smaller the proportion of its posterior density in the KLA fusion should be; and vice versa.

Simulations
The main goal of the simulations is to verify the following two points under different SNR conditions. First, our method conducts the sensor selection more effectively than the CS divergence based methods for the decentralized large-scale MTT network. This case is much more obvious when the sensors have different observation performance. Second, the coordinate descent method significantly shortens the calculation time of genetic algorithm at the expense of slight loss in tracking accuracy. To highlight these, the specific scenarios, including the multi-target dynamic model, sensor network architecture, observation model for SNs and so forth, are designed as follows.

Multiple targets move in a constant velocity (CV) model [49] over a two-dimensional region
where F CV and Q are the transition matrix and process noise covariance matrix for unlabeled state, where ∆ is the sampling interval, q Q is the process noise standard deviation. In this example, ∆ = 10 s, q Q = 0.002 km/s 2 and survival probability p s (x) = 0.95.
Each SN of LFC2 and LFC6 only receives the distance measurement of target. Its h s (x) and R s (x) are h s (x) = x − u s , s ∈ N 2 or s ∈ N 6 (73) Each SN of LFC3 and LFC7 only receives the angle measurement of target. Its h s (x) and R s (x) are Each SN of LFC4 and LFC8 receives distance and Doppler measurements of target. Its h s (x) and .
In this example, there are three constraints for the sensor selection optimization of (17), which are C1: Due to the limitations of communication bandwidth, energy consumption, computation capacity and storage space, the LFC j can only select K j SNs at each scan (K j ≤ N j ), In this example, if j = 1, 4, 5, 8, then K j = 5; otherwise K j = 8. C2: The field of view (FoV) of the SN s is modeled as a circular area with the center u s and radius The FoVs of different SNs can be overlapped.
To ensure that the FoVs of the SN set S j can totally cover the region A, it is required that In this example, if j = 1, 4, 5, 8, then ρ s = 30 km; otherwise ρ s = 20 km. C3: To avoid mutual interference between the homogeneous SNs belonging to the same LFC, the distance between any two SNs in S j must be not smaller than the threshold D j , In this example, D j = 5 km for j = 1, . . . , 8. According to the objective function and calculation method used for sensor selection, our algorithm is abbreviated as LA bound with coordinate descent. It is firstly implemented by the SMC technique and coded by MATLAB R2018a. Each p I (·, ) involved in the Mδ-GLMB density is approximated by 500 particles on average. In this example, the maximum numbers of targets and measurements of each SN per scan are set as 25 and 200, the cut-off is set as c = 1000 m. The algorithm is testing on a desktop with the CPU of AMD Ryzen 7 2700X and 64 GB RAM. We conduct 500 MC simulations, each of which includes T = 25 scans (a total of 250 s). In these simulations, the target tracks (including the instants of birth and death), clutter and measurements originating from targets are generated independently according to the aforementioned models.
We firstly present the result of sensor selection obtained by the algorithm in one simulation. Figure 3 shows the target trajectories in the simulation, where a total of 15 targets are generated at different instants and locations. For easy description, the targets are named as T1~T15 in turn. The name and survival period of each target are marked at the start point of the target. During the surveillance period, the number of targets at the initial time is the least (3 targets). The number of targets at the 15th~18th scans is the most (13 targets). The target T1 intersects with the target T10 at the 15th scan.
tracks (including the instants of birth and death), clutter and measurements originating from targets are generated independently according to the aforementioned models.
We firstly present the result of sensor selection obtained by the algorithm in one simulation. Figure 3 shows the target trajectories in the simulation, where a total of 15 targets are generated at different instants and locations. For easy description, the targets are named as T1~T15 in turn. The name and survival period of each target are marked at the start point of the target. During the surveillance period, the number of targets at the initial time is the least (3 targets). The number of targets at the 15th~18th scans is the most (13 targets). The target T1 intersects with the target T10 at the 15th scan.  Figure 4a~f show the results of sensor selection at the 1st, 5th, 10th, 15th, 20th and 25th scans. It can be seen that the SNs selected by each LFC will change adaptively with the multi-target movement versus time. Specifically, in order to minimize the bound in Theorem 1, most of the selected SNs locate in the regions closer to the survival targets at each scan. A few SNs which are far from the survival targets are selected to satisfy Constraint C2. Moreover, due to Constraint C3, the homogeneous SNs of each LFC cannot be excessively concentrated in a small region.  Figures 4a~4f show the results of sensor selection at the 1st, 5th, 10th, 15th, 20th and 25th scans. It can be seen that the SNs selected by each LFC will change adaptively with the multi-target movement versus time. Specifically, in order to minimize the bound in Theorem 1, most of the selected SNs locate in the regions closer to the survival targets at each scan. A few SNs which are far from the survival targets are selected to satisfy Constraint C2. Moreover, due to Constraint C3, the homogeneous SNs of each LFC cannot be excessively concentrated in a small region.     In order to further verify the performance of this method in tracking accuracy and computational time, it is compared with the methods of LA bound with genetic algorithm, CS divergence with genetic algorithm and Random selection under the same test platform. In the genetic algorithm, the population size is 50, crossing rate is 0.9, mutation rate is 0.001, elite rate is 0.04 and the maximum number of iterations is 500. In the CS divergence method, the objective function of the LFC j is is unknown before sensor selection, the expected value rather than real value of the CS divergence is applied in (82). In order to calculate the expected value, the MC integration based on PIMS also needs to be used here. For comparison, both the OSPA and LA metrics are used to measure the error of multi-target position estimates. Figures 5a and 5b present the 500 MC averages of the OSPA and LA errors for the four methods. Note that the error here is selected as the average of all LFCs because of Remark 1.  In order to further verify the performance of this method in tracking accuracy and computational time, it is compared with the methods of LA bound with genetic algorithm, CS divergence with genetic algorithm and Random selection under the same test platform. In the genetic algorithm, the population size is 50, crossing rate is 0.9, mutation rate is 0.001, elite rate is 0.04 and the maximum number of iterations is 500. In the CS divergence method, the objective function of the LFC j is where D CS (φ, ϕ) denotes the CS divergence between the densities φ and ϕ. [31] presents the specific form of the CS divergence when φ and ϕ are both GLMB densities. Since the posterior density π j X|Z S j is unknown before sensor selection, the expected value rather than real value of the CS divergence is applied in (82). In order to calculate the expected value, the MC integration based on PIMS also needs to be used here. For comparison, both the OSPA and LA metrics are used to measure the error of multi-target position estimates. is unknown before sensor selection, the expected value rather than real value of the CS divergence is applied in (82). In order to calculate the expected value, the MC integration based on PIMS also needs to be used here. For comparison, both the OSPA and LA metrics are used to measure the error of multi-target position estimates. Figures 5a and 5b present the 500 MC averages of the OSPA and LA errors for the four methods. Note that the error here is selected as the average of all LFCs because of Remark 1.  Figure 5 shows that both the averaged OSPA and LA errors from all the four methods decrease with time. Furthermore, the LA errors are always larger than the relevant OSPA errors. The reason for this has been explained in Subsection 2.3. In both of Figures 5(a) and 5(b), the errors from Random selection are always the largest. The next is the errors from CS divergence with genetic algorithm. The errors from LA bound with genetic algorithm are always the smallest. The errors from LA bound with coordinate descent are slightly larger than those of LA bound with genetic algorithm. Compared with Random selection, the errors of the other three algorithms are approximately reduced by 40%, 60% and 55%, respectively. Obviously, the two LA bound based methods outperform the CS divergence based method in MTT accuracy. There are three reasons for this: 1) The LA bound has a clearer physical meaning than the CS divergence. This is because that the former indicates the achievable optimal MTT accuracy with labeled RFS state. In contrast, the latter is not directly related to the MTT accuracy since maximizing the CS divergence in (82) cannot guarantee to minimize the OSPA or LA error.
2) The CS divergence cannot provide a basis for setting the weights of KLA rule as the LA bound. Therefore, in the CS divergence based method, the KLA weights can only be set to the same by convention [32]. However, in this example the MTT accuracy of different LFCs is probably not the same because of the distinguishing observation performance of the diverse SNs. If the KLA weights are set to the same without discrimination, the fusion efficiency will decline dramatically in this case.
3) In the step of sensor selection optimization, the coordinate descent method may trap in one of local optimums. In contrast, the genetic algorithm may jump out of local optimums with certain   Figure 5 shows that both the averaged OSPA and LA errors from all the four methods decrease with time. Furthermore, the LA errors are always larger than the relevant OSPA errors. The reason for this has been explained in Section 2.3. In both of Figure 5a,b, the errors from Random selection are always the largest. The next is the errors from CS divergence with genetic algorithm. The errors from LA bound with genetic algorithm are always the smallest. The errors from LA bound with coordinate descent are slightly larger than those of LA bound with genetic algorithm. Compared with Random selection, the errors of the other three algorithms are approximately reduced by 40%, 60% and 55%, respectively. Obviously, the two LA bound based methods outperform the CS divergence based method in MTT accuracy. There are three reasons for this: 1) The LA bound has a clearer physical meaning than the CS divergence. This is because that the former indicates the achievable optimal MTT accuracy with labeled RFS state. In contrast, the latter is not directly related to the MTT accuracy since maximizing the CS divergence in (82) cannot guarantee to minimize the OSPA or LA error.
2) The CS divergence cannot provide a basis for setting the weights of KLA rule as the LA bound. Therefore, in the CS divergence based method, the KLA weights can only be set to the same by convention [32]. However, in this example the MTT accuracy of different LFCs is probably not the same because of the distinguishing observation performance of the diverse SNs. If the KLA weights are set to the same without discrimination, the fusion efficiency will decline dramatically in this case.
3) In the step of sensor selection optimization, the coordinate descent method may trap in one of local optimums. In contrast, the genetic algorithm may jump out of local optimums with certain probability by its randomness. This leads that the MTT accuracy of the former is a litter worse than the latter.
On the other hand, the computational cost of a method is in general measured by its CPU run time. Then, the averaged CPU run time per scan for LA bound with coordinate descent, LA bound with genetic algorithm and CS divergence with genetic algorithm are 0.62s, 6.18s and 5.69s, respectively. Note that the run time here is also selected as the average of all LFCs because of Remark 1. Random selection obviously does not need the optimization for sensor selection, so it has no time consumption of this step. It can be seen from this that although the coordinate descent method is slightly worse than the genetic algorithm in tracking accuracy, it significantly shortens the calculation time of sensor selection optimization. In addition, the time consumption of CS divergence with genetic algorithm is slightly less than that of LA bound with genetic algorithm. This case indicates that the computational cost of the proposed bound is a little larger than that of the CS divergence expectation in (82).
In order to show the influence of different SNR on our method, the clutter rate and detection probability of each SN are changed into (λ = 40, p d = 0.85), (λ = 60, p d = 0.75), (λ = 80, p d = 0.65) and (λ = 100, p d = 0.55). Tables 1 and 2 present the final values of the OSPA and LA errors for the four methods in each scenario after 500 MC run average. Moreover, the CPU times consumed by the sensor selection optimization for the first three methods are nearly the same for all the SNR scenarios. This is because that the sensor selection is irrelevant to the specific measurement realizations since it must be completed before the measurements are received.
It can be seen from Algorithm 1 and 2 that as clutter density increases and detection probability decreases, 1) The OSPA and LA errors of all the four methods increase in different degrees but the size order of them is always the same as that of Scenario 1; 2) Taking the error of Random selection as the benchmark, the improvement ratio of CS divergence with genetic algorithm is gradually reduced from about 40% to about 20%. Meanwhile, the improvement ratios of the two LA bound based methods are, respectively, maintained at about 60% and 55%.
The two points reflect that the lower the SNR is, the worse the sensor selection efficiency and tracking accuracy of the CS divergence based methods become. By contrast, the tracking accuracy of the LA bound based methods always maintains a good improvement ratio for all the SNR scenarios.
Assuming that the simulation scenarios remain unchanged, the above four methods are re-implemented by the GM technique where the non-linear likelihood of each SN is approximated by the EK filter,Ĥ s ≈ ∂h s (x) ∂x x=x + ;R s ≈ R s (x + ) The pruning and merging technology [11] are used to manage the GM terms. Let the number of the GM terms approximating to each p I (·, ) be no more than 30. The thresholds for merging, pruning and state extraction are set to 4, 10 −4 and 0.5, respectively. The simulation results of the GM implementation are basically consistent with those of the SMC implementation, except that the OSPA or LA error increases by about 8%. But the calculation time for sensor selection decreases by about 70%.

Conclusions and Future Work
A sensor selection optimization algorithm is proposed for the decentralized large-scale MTT network under the labeled RFS framework. The LA metric defined in this paper is used to measure the error between the labeled multi-target states and their estimates. The lower bound of the LA metric based MSE is taken as the cost function of sensor selection. The bound is derived by the information inequality and then, implemented by the SMC or GM technique. Then, the coordinate descent method is used to reduce the computational cost of sensor selection. Simulation results show that when the sensors of the decentralized network have different observation performance, our method outperforms the CS divergence based sensor selection algorithm in MTT accuracy.
Our future work will focus on the following two aspects: 1) Extend the proposed method to the cases of asynchronous measurement or correlated measurement noise; 2) Reference [50] has presented a very efficient implementation of the GLMB filter with linear complexity in the number of measurements and this filter has been demonstrated to handle over 1 million tracks simultaneously [51]. Therefore, it would be very helpful to improve our current study by the use of the methods proposed in Reference [50,51].