# Multiradar Joint Tracking of Cluster Targets Based on Graph-LSTMs

Authors: Xirui Xue; Shucai Huang; Daozhi Wei; Jiahao Xie
Journal: Journal of Sensors
 (2022)

---

## Abstract

The cluster target brings a serious challenge to the traditional multisensor multitarget tracking algorithm because of its large number of members and the cooperative interaction between members. Using multiradar joint tracking cluster target is an alternative method to solve the problem of cluster target tracking, but it inevitably brings the problem of radar-target assignment and tracking information fusion. Aiming at the problem of radar-target assignment and tracking information fusion, a joint tracking method based on graph-long short-term memory neural nets (Graph-LSTMs) is proposed. Firstly, we use multivariable stochastic differential equations (SDE) to model the cooperative interaction of cluster members and transform the derived state space model of cluster members into the same form as the constant velocity (CV) motion model, and the target state equation of cluster which can be used for Bayesian filtering iteration is established. Secondly, based on the detection relationship between radars and cluster members, we introduce the detection confirmation matrix and propose a radar-target assignment method to achieve multiple measurements of single member and detection coverage of all cluster members. Then, each radar usesδ-GLMB filter to estimate the motion state of the assigned targets. Finally, on the basis of spatial discretization, the labels of multiple estimates of cluster member states are obtained. We use the designed Graph-LSTMs to learn the cooperative relationship between target states to fuse the labels and obtain better tracking effect. The experimental results show that the proposed method effectively simulates the cluster motion and realizes the joint estimation of cluster target motion state by multiradar. Our method makes up for the defect that a single radar cannot stably track adjacent multiple targets and achieves better estimation fusion effect than the expectation-maximization (EM) algorithm and mean method.


---

## Body

## 1. Introduction

The cluster target represented by UAV cluster is the research hotspot of target tracking. It has different tracking characteristics from multitarget and extended target. Compared with multitarget, there is cooperative interaction among cluster members, and the movement of members is consistent. However, due to the large number of members and the small distance between members, it is easy to lead to serious crossing of detection gates and difficult to carry out measurement association. Compared with the extended target, the cluster target is “nonrigid,” and the size and shape of the target are easy to change.At present, according to whether sensors can resolve cluster members, cluster target tracking can be divided into two research routes: cluster member tracking and overall cluster tracking. When the cluster is far away from the sensor and the members cannot be completely resolved, the common method is to model the cluster extension shape as an ellipse that is represented by a random matrix obeying the inverse Wishart distribution [1]. In this method, there is no need to consider the collaborative interaction between members but to track the cluster as a whole.When the members become resolvable to the sensor, the problem to be solved in cluster target tracking is the accurate state estimation of each cluster member. When solving this problem, the comprehensive application of various measurement methods [2] can ensure the resolution of the target to the greatest extent. For an example, multiradar networking, partition detection, and data fusion in the detection process are common methods in the tracking of UAV cluster. However, when using these methods, we need to solve not only the problem of accurate tracking of individual member but also the data fusion problem inevitably brought by multisensor joint tracking [3].This paper would focus on the tracking of resolvable cluster targets. All clusters in this paper refer to resolvable clusters. The first step to achieve accurate cluster target tracking is to obtain the prior information of cluster motion. Therefore, it is necessary to accurately establish the cluster target motion model. In terms of the general law of cluster movement, Reynolds summarizes that the cooperative rule is three different actions of separation, alignment, and cohesion between cluster members based on the distance from each other [4]. Viscek et al. study the speed consistency of cluster members and put forward the Viscek model [5]. In the model, they believe that the velocity consistency is achieved because each member adjusts the velocity direction to the average value of the velocity direction of other members in its neighborhood. Couzin et al. study the phenomenon of effective guidance and group decision-making and discuss the splitting phenomenon of clusters [6]. However, clusters described only by the above rules are prone to fracture due to the dispersion of member positions. Under the above rules, Olfati-Saber [7] adds the virtual leader and proposes the Olfati-Saber model to solve the problem of cluster fracture. The typical target we aim at in this paper is UAV cluster. The main description methods of UAV cluster motion include pilot following method [8], virtual structure method [9], behavior-based method [10], graph-based method [11], and consistency theory [12]. At present, most of the methods to describe UAV cluster motion are from the perspective of cluster control, which is difficult to adapt to the Bayesian iterative estimation process. In order to solve this problem, we try to deduce the motion model of a single-cluster member described in the state space.Radar-target assignment is another problem to be considered in the joint tracking of cluster members using multiradar [13, 14]. Under the limited radar resources, it is a contradiction to improve the detection density and ensure the full coverage of radar system to all cluster members. In this paper, the concept of detection confirmation matrix is introduced, and radar-target assignment problem is briefly discussed based on the concept.Aiming at the multitarget tracking (MTT) problem in the clutter environment, the multitarget tracking algorithm represented by the probabilistic data association [15] uses tracking gate to distinguish measurement and clutter and preliminarily solves the correlation problem between measurement and target state. Reference [16] proposes a multitarget tracking algorithm based on maximum entropy fuzzy C-means clustering joint probabilistic data association, which avoids confirmation matrix splitting and reduces the computational load of the JPDA algorithm. However, due to the large number of members in the cluster, the data association is very difficult, and only relying on the confirmation matrix for measuring data association is difficult to meet the requirements of tracking accuracy. At this time, due to the avoidance of data association, tracking methods based on random finite sets (RFS) gradually rise and are widely used [17]. These methods include probability hypothesis density (PHD) filter [18, 19], cardinalized probability hypothesis density (CPHD) filter [20, 21], multitarget multi-Bernoulli (MeMBer) filter [22, 23], and generalized label multi-Bernoulli (GLMB) filter [24, 25]. Reference [26] uses the MCMC particle filter algorithm to track the cooperative cluster with few members and solves the problems of cluster structure inference and joint estimation of cluster and member state. Reference [27] uses δ-GLMB filter to track cluster targets with splitting, merging, and reorganization behavior. Usually, RFS-based tracking methods set the detection probability as a constant a priori, but this does not conform to the actual sensor detection process. Reference [28] uses active sonar equation to model the detection probability and proposes the Pd-GM-PHD and Pd-GM-CPHD algorithms to achieve more accurate tracking effect. The existing filter based on RFS avoids the problem of data association and realizes the joint estimation of the number of targets and states. Compared with the probabilistic data association algorithm, the number of targets that can be tracked is increased. However, these methods cannot meet the requirements of the number of tracking cluster members only by using a single sensor, and they do not consider the problem of tracking information fusion caused by using multiple sensors in the actual situation.The traditional method of multiradar tracking information fusion is to approximate the posterior distribution of the target by Gaussian mixture. A fusion method of Gaussian mixture distribution using the Chebyshev information is proposed in [29]. A consensus CPHD filter based on Gaussian mixture distribution is designed in [30]. And a likelihood function distribution approximation method based on polynomial weighting is designed in [31]. However, these methods do not further consider the useful information that the relationship between targets may provide. In recent years, machine learning algorithms based on neural networks provide a new way for mining information in target tracking [32, 33].In this paper, in order to achieve better tracking effect, we try to use Graph-LSTMs [34–36] to realize the fusion of multiradar estimation information and mine the collaborative information among cluster members at the same time. The key contributions of this research are as follows.
(1)
We establish a stochastic differential equation model to describe the cooperative interaction relationship within the cluster and unify it under the framework of CV model, which would be used as the target state equation of the cluster in the iterative filtering process(2)
Considering the problem of multiradar detection of cluster targets, we propose a radar-target assignment method. Benefiting from the introduction of detection confirmation matrix, this method can not only ensure that multiple radar systems detect all cluster members but also realize multiple detection of the same member, so as to obtain multiple target state estimates(3)
Based on the cluster motion model proposed in the paper, the cluster motion trajectory dataset is simulated. We useδ-GLMB filter to realize the state estimation of radar to assigned targets(4)
The bipartite graph of radar detection relation is established and expanded into a tree structure. The space around the predicted value of the target state is discretized to obtain the label of the target state estimated value. The Graph-LSTM model is used to fuse the label of the estimated value, and a more accurate fusion value is obtained. Experimental results show that our method has more accurate fusion labels and smaller OSPA distance than the mean method and EM algorithmThe content of this paper is arranged as follows: Section2 introduces the construction method of the cooperative interaction model based on multivariable stochastic differential equation. Section 3 introduces the multiradar tracking method of cluster targets and puts forward a radar-target assignment method. On this basis, the state estimation of the assigned targets is realized by using the δ-GLMB filter. Section 4 introduces an estimation fusion method based on Graph-LSTMs. In Section 5, the cluster motion simulation is carried out, and the estimation fusion experiment is carried out using the simulation data set to verify the effectiveness of the proposed method.


## 2. Cooperative Interaction Motion Model

In this section, on the basis of cluster cooperation rules, the self-organizing cooperative interaction motion model of cluster is established, and its discrete form in state space is deduced. In the model, the motion of cluster members is directly related to other members in the neighborhood, and cluster members act according to the motion state of other members. Finally, the velocity of cluster members is gradually consistent, and there is no collision between members. Multivariable stochastic differential equation (SDE) would be an effective tool for us to establish cluster member cooperative interaction motion model [37].


### 2.1. Continuous-Time Collaborative Interaction Model

Considering that the motion state of cluster members is jointly determined by their own state and the states of other members in the neighborhood, the continuous-time collaborative interaction motion model of cluster memberi at time t is constructed as follows:
(1)dṡt,i=‐αṡt,i‐fṡt+βrist,idt+dWt,is+dBts,where st,i and ṡt,i, respectively, represent the position and velocity of the target i at time t, α represents the velocity control parameter, β represents the potential force control parameter, Wt,is represents the motion noise of the members, Bts represents the cluster motion noise, and fṡt represents the average velocity of the cluster where target i is located. Let the number of cluster members be N, and the expression of fṡt is as follows:
(2)fṡt=1N∑j=1Nṡt,j.The introduction offṡt makes the average velocity information of the cluster pass through the whole cluster, makes the topology of the cluster structure have weak connectivity, and realizes the gradual consistency of the velocity of members.rist,i represents the total potential force of member i by other members in the neighborhood, which can avoid collision between cluster members. Its expression is as follows:
(3)rist=∑∀j,j≠irst,i,st,j,where rst,i,st,j represents the potential force of member j to member i in the cluster, which is generated by the potential function Uj⋅. Uj⋅ needs to meet the following two properties. First, Uj⋅ should be a continuous, differentiable, and nonnegative function; second, Uj⋅ needs to obtain a unique minimum value at a certain required distance. In this paper, a definition of Uj⋅ is given as follows:
(4)Ujd=R11lndb+R12d+R12+R11db+R12d−db,d≤db,R21lndm−db+R22−d+dm+R22+R21dm−db+R22db−d,db<d<dm.Here,d=st,i−st,j and dr are the balance distance between two members. When the distance is less than dr, the potential force is repulsive force. Otherwise, it is attractive force. dm is the maximum distance that can generate potential force. The cluster potential field is shown in Figure 1. R11, R12, R21, and R22 are the parameters controlling the intensity of potential field. Under the above definition of potential function, as the negative gradient of potential function, the potential force rst,i,st,j is calculated as follows:
(5)rst,i,st,j=−gradUjd=R11d+R12−R11db+R12,d≤db,−R21dm−d+R22+R21dm−db+R22,db<d<dm.Figure 1
Schematic diagram of potential field.


### 2.2. Continuous-Time Joint Motion Model

Because the cooperative interaction model of cluster member contains the motion information of other members in cluster, it cannot be solved separately. The models of all members in the cluster need to be combined to obtain the analytical solution together.In 2-dimensional space, member positionst,i=xt,i,yt,i. Use Xt to represent the joint state of all members in the cluster; then, we can get Xt=xt,1,ẋt,1,yt,1,ẏt,1,⋯,xt,N,ẋt,N,yt,N,ẏt,NT. From (1), the linear differential equation of the joint state in continuous time can be obtained as follows:
(6)dXt=AXtdt+Rdt+CdWt+DdBt.MatrixA∈R4N×4N is defined as follows:
(7)A=A1A2⋯A2A2A1⋯A2⋮⋱⋮A2A2⋯A1,A1=010001−NαN0000010001−NαN,A2=00000αN000000000αN.MatrixR∈R4N×1 is defined as follows:
(8)R=0,βri,xst,0,βri,yst,⋯,0,βrN,xst,0,βrN,ystT.MatrixC∈R4N×2N is defined as follows:
(9)C=C10⋯00C1⋯0⋮⋱⋮00⋯C1,C1=00100001.MatrixD∈R4N×2 is defined as follows:
(10)D=C1T,C1T,⋯,C1TT.Motion noiseWt=Wt,1x,Wt,1y,⋯,Wt,Nx,Wt,NyT, and Bt=Btx,BtyT. According to the general practice, motion noise is usually modeled as white noise with Gaussian distribution, and their covariance matrices are QW=diagσx2,σy2,⋯,σx2,σy2 and QB=diagσg2,σg2, respectively. Because they have the same form, in order to simplify the model, the two noises can be combined into cluster joint motion noise Mt and its covariance matrix QM=diagσx2,σy2,⋯,σx2,σy2,σg2,σg2. Accordingly, after combining the two coefficient matrices, there is H=C,D. Equation (6) has the following form:
(11)dXt=AXtdt+Rdt+HdMt.The following analytical solution is obtained by solving (11).
(12)Xt+τ=eAτXt+∫h=th=t+τeAt+τ−hdh×R+∫h=th=t+τeAt+τ−hHdMt,where τ is the sampling time. Define the member joint state transition matrix Fcτ=eAτ. Define the potential force joint gain matrix Ecτ=∫h=th=t+τeAt+τ−hdh, and the noise joint covariance matrix is calculated as follows:
(13)Qcτ=∫h=th=t+τeAt+τ−hHQMHTeAt+τ−hTdh.It can be seen from (12) that Ecτ and Qcτ are only related to the sampling time τ. Therefore, the time discretization of (12) can be realized by making τ=1, and the following equation is obtained:
(14)Xk+1=FcXk+EcR+Γcwk.Equation (14) is the discrete-time joint state model of cluster, where Fc=Fc1 and Ec=Ec1; wk is the Gaussian white noise with covariance matrix Qc=Qc1; and Γc is the state noise factor matrix, which is used to convert the noise dimension and add the noise to the joint state Xk+1.


### 2.3. Discrete-Time State Space Model

After obtaining discrete-time joint state model of cluster, it is necessary to separate the state space model of a single member from the joint state model and further simplify it to a general form suitable for iterative tracking algorithm. By analyzing the coefficient matrix and covariance matrix, we can see that they have the following symmetric forms, which would be the basis for us to realize the above separation.(15)Fc=FcsFcc⋯FccFccFcs⋯Fcc⋮⋱⋮FccFcc⋯Fcs,Ec=EcsEcc⋯EccEccEcs⋯Ecc⋮⋱⋮EccEcc⋯Ecs,Qc=QcsQcc⋯QccQccQcs⋯Qcc⋮⋱⋮QccQcc⋯Qcs.Among them, the diagonal matrices with superscripts are related to the motion state of the members themselves, and other matrices with superscript c are related to the states of other members in the cluster.Further solving matrixFc shows that for any member i, its own state transition matrix Fcs and the cooperation matrix Fcc are as follows:
(16)Fcs=1f1000f200001f1000f2,Fcc=01−f1N−10001−f2N−1000001−f1N−10001−f2N−1.f1 and f2 are constants, which are related to the member number N and parameter α. Assuming that the influence of potential force on the motion state of member i at a certain time is only related to the potential force on member i at the current time, there are
(17)Ecs=112000100001120001.And Ectc=O4×4. O4×4 represents the fourth-order zero matrix.In the cluster movement dominated by collaborative behavior, the motion randomness of other members has little impact on the memberi and would offset each other due to synergy. The motion randomness of members is more expressed as endogenous motion uncertainty. Therefore, it can be considered that the noise covariance matrix of member i can be directly simplified to Qcts.Based on the above assumptions, the discrete-time state space model of cluster memberi can be given as follows:
(18)xk+1,i=Fcsxk,i+Fcc∑∀j,j≠ixk,j+ck,i+Γcswk,i,where
(19)ck,i=∑∀j,j≠ick,i,j=βri,x2,βri,x,βri,y2,βri,yT,ck,i,j=βri,xxk,i,xk,j2,βri,xxk,i,xk,j,βri,yxk,i,xk,j2,βri,yxk,i,xk,jT.Γcs is the state noise factor matrix, and wk,i is the Gaussian noise with covariance matrix Qcts.Although the above equation separates the motion states of memberi from other members in the cluster, the number of members still affects the values of f1 and f2.In order to completely separate the influence of collaborative interaction, simple transformation is carried out under the framework of constant velocity (CV) motion model, and the discrete-time state space model of cluster members is proposed as follows:(20)xk+1,i=Fkxk,i+wk,i0.Here,(21)Fk=1100010000110001.wk,i0 is a new collaborative noise, which includes all the effects of cooperative interaction and motion noise on a cluster member motion state. Its calculation method is shown in (22).
(22)wk,i0=Δck,i+Γkwk,i.Here,Γk=Γcs and Δck,i are calculated as follows:
(23)Δck,i=∑∀j,j≠iωkck,i,jωk−Fkẋk,i−ẋk,j,(24)ωk=1N−1I4×4−diagf1,f2,f1,f2.In (23), ẋk,i is the velocity vector and ẋk,i=0,ẋk,i,0,ẏk,iT.It is assumed that the interaction is a random variable obeying to Gaussian distribution, that is,ck,i,j∼Nc¯k,i,j,Sk,i,j. Then, the collaborative noise wk,i0 can also be modeled as Gaussian distribution, that is, wk,i0∼Nμk,i0,Qk,i0.Here,(25)μk,i0=∑∀j,j≠iωkc¯k,i,jωk−Fkẋ^k,i−ẋ^k,j,Qk,i0=∑∀j,j≠iSk,i,j+ωkFk∑∀j,j≠iPk,j−Pk,iFkTωkT+ΓkQcsΓkT,where Pk,i is the state covariance matrix of ẋk,i. The whole process of model construction is shown in Figure 2. The model construction has experienced the transformation from member to cluster and then to member and realizes the discretization of continuous-time model.Figure 2
Model construction process.However, even if model (20) is obtained, it is difficult to directly apply the traditional tracking algorithm. There are two reasons. On the one hand, the noise wk,i0 is not Gaussian white noise, but the noise with time-varying mean and covariance. On the other hand, due to the large number of members in the cluster and the limited number of radar channels, it is often difficult for a single radar to track the whole cluster, and it is prone to the problem of measurement-track association. In order to solve the above problems, multiradar joint tracking becomes an alternative method.


## 2.1. Continuous-Time Collaborative Interaction Model

Considering that the motion state of cluster members is jointly determined by their own state and the states of other members in the neighborhood, the continuous-time collaborative interaction motion model of cluster memberi at time t is constructed as follows:
(1)dṡt,i=‐αṡt,i‐fṡt+βrist,idt+dWt,is+dBts,where st,i and ṡt,i, respectively, represent the position and velocity of the target i at time t, α represents the velocity control parameter, β represents the potential force control parameter, Wt,is represents the motion noise of the members, Bts represents the cluster motion noise, and fṡt represents the average velocity of the cluster where target i is located. Let the number of cluster members be N, and the expression of fṡt is as follows:
(2)fṡt=1N∑j=1Nṡt,j.The introduction offṡt makes the average velocity information of the cluster pass through the whole cluster, makes the topology of the cluster structure have weak connectivity, and realizes the gradual consistency of the velocity of members.rist,i represents the total potential force of member i by other members in the neighborhood, which can avoid collision between cluster members. Its expression is as follows:
(3)rist=∑∀j,j≠irst,i,st,j,where rst,i,st,j represents the potential force of member j to member i in the cluster, which is generated by the potential function Uj⋅. Uj⋅ needs to meet the following two properties. First, Uj⋅ should be a continuous, differentiable, and nonnegative function; second, Uj⋅ needs to obtain a unique minimum value at a certain required distance. In this paper, a definition of Uj⋅ is given as follows:
(4)Ujd=R11lndb+R12d+R12+R11db+R12d−db,d≤db,R21lndm−db+R22−d+dm+R22+R21dm−db+R22db−d,db<d<dm.Here,d=st,i−st,j and dr are the balance distance between two members. When the distance is less than dr, the potential force is repulsive force. Otherwise, it is attractive force. dm is the maximum distance that can generate potential force. The cluster potential field is shown in Figure 1. R11, R12, R21, and R22 are the parameters controlling the intensity of potential field. Under the above definition of potential function, as the negative gradient of potential function, the potential force rst,i,st,j is calculated as follows:
(5)rst,i,st,j=−gradUjd=R11d+R12−R11db+R12,d≤db,−R21dm−d+R22+R21dm−db+R22,db<d<dm.Figure 1
Schematic diagram of potential field.


## 2.2. Continuous-Time Joint Motion Model

Because the cooperative interaction model of cluster member contains the motion information of other members in cluster, it cannot be solved separately. The models of all members in the cluster need to be combined to obtain the analytical solution together.In 2-dimensional space, member positionst,i=xt,i,yt,i. Use Xt to represent the joint state of all members in the cluster; then, we can get Xt=xt,1,ẋt,1,yt,1,ẏt,1,⋯,xt,N,ẋt,N,yt,N,ẏt,NT. From (1), the linear differential equation of the joint state in continuous time can be obtained as follows:
(6)dXt=AXtdt+Rdt+CdWt+DdBt.MatrixA∈R4N×4N is defined as follows:
(7)A=A1A2⋯A2A2A1⋯A2⋮⋱⋮A2A2⋯A1,A1=010001−NαN0000010001−NαN,A2=00000αN000000000αN.MatrixR∈R4N×1 is defined as follows:
(8)R=0,βri,xst,0,βri,yst,⋯,0,βrN,xst,0,βrN,ystT.MatrixC∈R4N×2N is defined as follows:
(9)C=C10⋯00C1⋯0⋮⋱⋮00⋯C1,C1=00100001.MatrixD∈R4N×2 is defined as follows:
(10)D=C1T,C1T,⋯,C1TT.Motion noiseWt=Wt,1x,Wt,1y,⋯,Wt,Nx,Wt,NyT, and Bt=Btx,BtyT. According to the general practice, motion noise is usually modeled as white noise with Gaussian distribution, and their covariance matrices are QW=diagσx2,σy2,⋯,σx2,σy2 and QB=diagσg2,σg2, respectively. Because they have the same form, in order to simplify the model, the two noises can be combined into cluster joint motion noise Mt and its covariance matrix QM=diagσx2,σy2,⋯,σx2,σy2,σg2,σg2. Accordingly, after combining the two coefficient matrices, there is H=C,D. Equation (6) has the following form:
(11)dXt=AXtdt+Rdt+HdMt.The following analytical solution is obtained by solving (11).
(12)Xt+τ=eAτXt+∫h=th=t+τeAt+τ−hdh×R+∫h=th=t+τeAt+τ−hHdMt,where τ is the sampling time. Define the member joint state transition matrix Fcτ=eAτ. Define the potential force joint gain matrix Ecτ=∫h=th=t+τeAt+τ−hdh, and the noise joint covariance matrix is calculated as follows:
(13)Qcτ=∫h=th=t+τeAt+τ−hHQMHTeAt+τ−hTdh.It can be seen from (12) that Ecτ and Qcτ are only related to the sampling time τ. Therefore, the time discretization of (12) can be realized by making τ=1, and the following equation is obtained:
(14)Xk+1=FcXk+EcR+Γcwk.Equation (14) is the discrete-time joint state model of cluster, where Fc=Fc1 and Ec=Ec1; wk is the Gaussian white noise with covariance matrix Qc=Qc1; and Γc is the state noise factor matrix, which is used to convert the noise dimension and add the noise to the joint state Xk+1.


## 2.3. Discrete-Time State Space Model

After obtaining discrete-time joint state model of cluster, it is necessary to separate the state space model of a single member from the joint state model and further simplify it to a general form suitable for iterative tracking algorithm. By analyzing the coefficient matrix and covariance matrix, we can see that they have the following symmetric forms, which would be the basis for us to realize the above separation.(15)Fc=FcsFcc⋯FccFccFcs⋯Fcc⋮⋱⋮FccFcc⋯Fcs,Ec=EcsEcc⋯EccEccEcs⋯Ecc⋮⋱⋮EccEcc⋯Ecs,Qc=QcsQcc⋯QccQccQcs⋯Qcc⋮⋱⋮QccQcc⋯Qcs.Among them, the diagonal matrices with superscripts are related to the motion state of the members themselves, and other matrices with superscript c are related to the states of other members in the cluster.Further solving matrixFc shows that for any member i, its own state transition matrix Fcs and the cooperation matrix Fcc are as follows:
(16)Fcs=1f1000f200001f1000f2,Fcc=01−f1N−10001−f2N−1000001−f1N−10001−f2N−1.f1 and f2 are constants, which are related to the member number N and parameter α. Assuming that the influence of potential force on the motion state of member i at a certain time is only related to the potential force on member i at the current time, there are
(17)Ecs=112000100001120001.And Ectc=O4×4. O4×4 represents the fourth-order zero matrix.In the cluster movement dominated by collaborative behavior, the motion randomness of other members has little impact on the memberi and would offset each other due to synergy. The motion randomness of members is more expressed as endogenous motion uncertainty. Therefore, it can be considered that the noise covariance matrix of member i can be directly simplified to Qcts.Based on the above assumptions, the discrete-time state space model of cluster memberi can be given as follows:
(18)xk+1,i=Fcsxk,i+Fcc∑∀j,j≠ixk,j+ck,i+Γcswk,i,where
(19)ck,i=∑∀j,j≠ick,i,j=βri,x2,βri,x,βri,y2,βri,yT,ck,i,j=βri,xxk,i,xk,j2,βri,xxk,i,xk,j,βri,yxk,i,xk,j2,βri,yxk,i,xk,jT.Γcs is the state noise factor matrix, and wk,i is the Gaussian noise with covariance matrix Qcts.Although the above equation separates the motion states of memberi from other members in the cluster, the number of members still affects the values of f1 and f2.In order to completely separate the influence of collaborative interaction, simple transformation is carried out under the framework of constant velocity (CV) motion model, and the discrete-time state space model of cluster members is proposed as follows:(20)xk+1,i=Fkxk,i+wk,i0.Here,(21)Fk=1100010000110001.wk,i0 is a new collaborative noise, which includes all the effects of cooperative interaction and motion noise on a cluster member motion state. Its calculation method is shown in (22).
(22)wk,i0=Δck,i+Γkwk,i.Here,Γk=Γcs and Δck,i are calculated as follows:
(23)Δck,i=∑∀j,j≠iωkck,i,jωk−Fkẋk,i−ẋk,j,(24)ωk=1N−1I4×4−diagf1,f2,f1,f2.In (23), ẋk,i is the velocity vector and ẋk,i=0,ẋk,i,0,ẏk,iT.It is assumed that the interaction is a random variable obeying to Gaussian distribution, that is,ck,i,j∼Nc¯k,i,j,Sk,i,j. Then, the collaborative noise wk,i0 can also be modeled as Gaussian distribution, that is, wk,i0∼Nμk,i0,Qk,i0.Here,(25)μk,i0=∑∀j,j≠iωkc¯k,i,jωk−Fkẋ^k,i−ẋ^k,j,Qk,i0=∑∀j,j≠iSk,i,j+ωkFk∑∀j,j≠iPk,j−Pk,iFkTωkT+ΓkQcsΓkT,where Pk,i is the state covariance matrix of ẋk,i. The whole process of model construction is shown in Figure 2. The model construction has experienced the transformation from member to cluster and then to member and realizes the discretization of continuous-time model.Figure 2
Model construction process.However, even if model (20) is obtained, it is difficult to directly apply the traditional tracking algorithm. There are two reasons. On the one hand, the noise wk,i0 is not Gaussian white noise, but the noise with time-varying mean and covariance. On the other hand, due to the large number of members in the cluster and the limited number of radar channels, it is often difficult for a single radar to track the whole cluster, and it is prone to the problem of measurement-track association. In order to solve the above problems, multiradar joint tracking becomes an alternative method.


## 3. Multiradar Joint Tracking Method

In order to solve the problem of obtaining measurement information in the process of cluster target detection and reduce the impact of collaborative noise on cluster tracking, the following solutions are proposed in this paper. Firstly, as shown in Figure3, multiple radars are used to detect a cluster at the same time, and each radar obtains some measurement information of cluster, respectively. The sum of all measurement information needs to ensure full coverage of all cluster members. Secondly, the states of the members are estimated by using the multitarget filter. Finally, the state estimation of each member is obtained by information fusion. In order to ensure the effect of information fusion, each member needs to estimate as much as possible.Figure 3
Multiradar detection.


### 3.1. Acquisition of Detection Information

#### 3.1.1. Detection Relationship Construction

In order to reflect the detection relationship between multiradar and cluster members, we introduce the concept of detection confirmation matrixΩ, which is defined as
(26)Ω=ωst,s=1,2,⋯,M,t=1,2,⋯,T,where ωst is a binary variable. ωst=1 indicates that cluster member t is detected by radar s. And ωst=0 indicates that cluster member t is not detected by radar s. The average detection density Ad is defined as the average detection times of all members in the cluster. Ω needs to meet the following three assumptions:
(1)
Each member is detected by at least one radar, i.e.,(27)∑s=1Mωst≥1,t=1,2,⋯,N(2)
The number of members detected by each radar is less than the number of radar channelsNc, i.e.,(28)∑t=1Nωst≤Nc,s=1,2,⋯,M(3)
The average detection densityAd should be greater than the given value A¯d, i.e.,(29)Ad=∑t=1N∑s=1MωstN≥A¯d


#### 3.1.2. Radar Measurement Model

When the radar detection relationship is determined, the measurement information model of the target follows the general measurement model. For memberi, there are the following measurement models:
(30)zk+1,i=Hk+1xk+1,i+vk+1,i,where
(31)Hk+1=10000010.The measurement noisevk+1,i is a Gaussian white noise with covariance matrix Rk+1.Equations (20) and (30) together constitute the system equation of the cluster target tracking algorithm. Considering the problem of radar-target assignment and tracking information fusion, the filter needs to have the ability to propagate multiple hypotheses in the label set of different trajectories and output the target label. We choose the δ-generalized label multi-Bernoulli (δ-GLMB) filter [38, 39] for target tracking.


### 3.2. Application ofδ-GLMB Filter

In the process of cluster target tracking, an accurate state estimation algorithm should include the estimation of cooperative interaction. However, we usually cannot obtain the priori information of cooperative interaction. Therefore, we use the method of multiradar estimation fusion to solve this problem, and each radar runs aδ-GLMB filter.Theδ-GLMB filter recursively transmits the δ-GLMB filtering density forward over time through Bayesian prediction and update equation. The difference from GLMB is that when predicting, δ-GLMB is summed on the label space I+,ϑ∈FL+×Ξ at the next time. Here, L+=L∪B. From the perspective of multitarget tracking, this is more intuitive because it shows how the prediction introduces a new target label. δ-GLMB is represented by parameter set wI,ϑ,pϑ:I+,ϑ∈FL+×Ξ. The parameter set of δ-GLMB can be regarded as an enumeration of all assumptions and their related weights and track density Ih,ϑh,wh,phh=1H.


#### 3.2.1. Prediction

At timek, the distribution of the random finite set of δ-GLMB is as follows:
(32)πX=ΔX∑I,ϑ∈FL×ΞwI,ϑδILXpϑX.Then, the multitarget prediction density at timek+1 is
(33)π+X+=ΔX+∑I+,ϑ∈FL+×Ξw+I+,ϑδI+LX+p+ϑX+,where
(34)w+I+,ϑ=wγI+∩BwSϑI+∩L,wsϑL=ηSϑL∑I⊇L1−ηSϑI−LwI,ϑ,ηSϑl=pS⋅,l,pϑ⋅,l,p+ϑx,l=1Llp+,Sϑx,l+1Blpγx,l,p+,Sϑx,l=pS⋅,lϕx⋅,l,pϑ⋅,lηSϑl.Here,wγI+∩B is the weight of newborn label I+∩B, and wSϑI+∩L is the weight of survival label. ϕx⋅,l is the single-target transfer density of trajectory l, and ηSϑl is the survival probability of trajectory l. pγx,l is the probability density of newborn state, and p+,Sϑx,l is the probability density of survival state.


#### 3.2.2. Update

After obtaining the multitarget prediction density at timek+1, it can be proved that the updated multitarget filter density is also δ-GLMB.The update process is shown as follows:(35)πXZ=ΔX∑I,ϑ∈FL×Ξ∑θ∈ΘIwI,ϑ,θZδILXpϑ,θ⋅ZX,where
(36)wI,ϑ,θZ=δθ−10:ZIηZϑ,θIwI,ϑ∑I,ϑ∈FL×Ξ∑θ∈ΘIδθ−10:ZIηZϑ,θIwI,ϑ∝ηZϑ,θIwI,ϑ,ηZϑ,θl=pϑ⋅,l,φZ⋅,l;θ,pϑ,θx,lZ=pϑx,lφZx,l;θηZϑ,θl.Here,ΘI represents a subset of the current association mapping with domain I. φZx,l;θ is defined as follows:
(37)φZx,l;θ=δ0θl1−pDx,l+1−δ0θlpDx,lgzθlx,lκzθl=pDx,lgzθlx,lκzθl,θl>0,1−pDx,l,θl=0.Here,pDx,l is the detection probability on x,l. gzθlx,l is the single-target likelihood function of zθl to x,l. κ⋅ is the Poisson clutter density.


#### 3.2.3. Realization ofδ-GLMB Filter

For cluster member motion model,pSx,l=pS, and ϕx+x,l=Nx+;Fkx,Qcs. Since the newborn density matrix is in the form of Gaussian mixture, if the single-target density is in the form of Gaussian mixture, there is
(38)pϑ⋅,l=∑i=1JϑlωiϑlNx;miϑl,Piθl.Then,(39)ηSϑl=pS,p+ϑx,l=1Ll∑i=1JϑlωiϑlNx;mS,iϑl,PS,iϑl+1Blpγl,where
(40)mS.iϑl=Fkmiϑl,PS,iϑl=FkPiϑlFkT+Qcs.When the motion model parameters are related to the label, we just need to substitutepS=pSl, Fk=Fkl, and Qcs=Qcsl into the above equation.In the update process,pDx,l=pD. For the measurement model shown in (30), the measurement likelihood is gzx,l=Nz;Hk+1x,Rk+1. Gaussian mixture expression provides the most general configuration for linear Gaussian model. It is assumed that each single-target density pϑ⋅,l is in the form of Gaussian mixture.
(41)pϑ⋅,l=∑n=1JϑlωnϑlNx;mnϑl,Pnϑl.Then,(42)ci,j=−lnpD∑n=1Jϑliωnϑliqnϑzj;li1−pDκzj.In addition, for the updated association processϑ,θ, there are
(43)ηZϑ,θl=∑n=1JϑlwZ,nϑ,θl,pϑ,θx,lZ=∑n=1JϑlwZ,nϑ,θlηZϑ,θlNx;mZ,nϑ,θl,Pnϑ,θl.Here,(44)wZ,nϑ,θl=ωnϑl×pDqnϑzθl;lκzθl,θl>0,1−pD,θl=0,qnϑz;l=Nz;Hk+1mnϑl,Hk+1PnϑlHk+1T+Rk+1,mZ,nϑ,θl=mnϑl+Gnϑ,θlzθl−Hk+1mnϑl,θl>0,mnϑl,θl=0,Pnϑ,θl=I−Gnϑ,θlHk+1Pnϑl,Gnϑ,θl=PnϑlHk+1THk+1PnϑlHk+1T+Rk+1−1,θl>0,0,θl=0.When the parameters of the measurement model are dependent on the labell, it is only necessary to substitute pD=pDl, Hk+1=Hk+1l, and Rk+1=Rk+1l into the above equation.


## 3.1. Acquisition of Detection Information

### 3.1.1. Detection Relationship Construction

In order to reflect the detection relationship between multiradar and cluster members, we introduce the concept of detection confirmation matrixΩ, which is defined as
(26)Ω=ωst,s=1,2,⋯,M,t=1,2,⋯,T,where ωst is a binary variable. ωst=1 indicates that cluster member t is detected by radar s. And ωst=0 indicates that cluster member t is not detected by radar s. The average detection density Ad is defined as the average detection times of all members in the cluster. Ω needs to meet the following three assumptions:
(1)
Each member is detected by at least one radar, i.e.,(27)∑s=1Mωst≥1,t=1,2,⋯,N(2)
The number of members detected by each radar is less than the number of radar channelsNc, i.e.,(28)∑t=1Nωst≤Nc,s=1,2,⋯,M(3)
The average detection densityAd should be greater than the given value A¯d, i.e.,(29)Ad=∑t=1N∑s=1MωstN≥A¯d


### 3.1.2. Radar Measurement Model

When the radar detection relationship is determined, the measurement information model of the target follows the general measurement model. For memberi, there are the following measurement models:
(30)zk+1,i=Hk+1xk+1,i+vk+1,i,where
(31)Hk+1=10000010.The measurement noisevk+1,i is a Gaussian white noise with covariance matrix Rk+1.Equations (20) and (30) together constitute the system equation of the cluster target tracking algorithm. Considering the problem of radar-target assignment and tracking information fusion, the filter needs to have the ability to propagate multiple hypotheses in the label set of different trajectories and output the target label. We choose the δ-generalized label multi-Bernoulli (δ-GLMB) filter [38, 39] for target tracking.


## 3.1.1. Detection Relationship Construction

In order to reflect the detection relationship between multiradar and cluster members, we introduce the concept of detection confirmation matrixΩ, which is defined as
(26)Ω=ωst,s=1,2,⋯,M,t=1,2,⋯,T,where ωst is a binary variable. ωst=1 indicates that cluster member t is detected by radar s. And ωst=0 indicates that cluster member t is not detected by radar s. The average detection density Ad is defined as the average detection times of all members in the cluster. Ω needs to meet the following three assumptions:
(1)
Each member is detected by at least one radar, i.e.,(27)∑s=1Mωst≥1,t=1,2,⋯,N(2)
The number of members detected by each radar is less than the number of radar channelsNc, i.e.,(28)∑t=1Nωst≤Nc,s=1,2,⋯,M(3)
The average detection densityAd should be greater than the given value A¯d, i.e.,(29)Ad=∑t=1N∑s=1MωstN≥A¯d


## 3.1.2. Radar Measurement Model

When the radar detection relationship is determined, the measurement information model of the target follows the general measurement model. For memberi, there are the following measurement models:
(30)zk+1,i=Hk+1xk+1,i+vk+1,i,where
(31)Hk+1=10000010.The measurement noisevk+1,i is a Gaussian white noise with covariance matrix Rk+1.Equations (20) and (30) together constitute the system equation of the cluster target tracking algorithm. Considering the problem of radar-target assignment and tracking information fusion, the filter needs to have the ability to propagate multiple hypotheses in the label set of different trajectories and output the target label. We choose the δ-generalized label multi-Bernoulli (δ-GLMB) filter [38, 39] for target tracking.


## 3.2. Application ofδ-GLMB Filter

In the process of cluster target tracking, an accurate state estimation algorithm should include the estimation of cooperative interaction. However, we usually cannot obtain the priori information of cooperative interaction. Therefore, we use the method of multiradar estimation fusion to solve this problem, and each radar runs aδ-GLMB filter.Theδ-GLMB filter recursively transmits the δ-GLMB filtering density forward over time through Bayesian prediction and update equation. The difference from GLMB is that when predicting, δ-GLMB is summed on the label space I+,ϑ∈FL+×Ξ at the next time. Here, L+=L∪B. From the perspective of multitarget tracking, this is more intuitive because it shows how the prediction introduces a new target label. δ-GLMB is represented by parameter set wI,ϑ,pϑ:I+,ϑ∈FL+×Ξ. The parameter set of δ-GLMB can be regarded as an enumeration of all assumptions and their related weights and track density Ih,ϑh,wh,phh=1H.


### 3.2.1. Prediction

At timek, the distribution of the random finite set of δ-GLMB is as follows:
(32)πX=ΔX∑I,ϑ∈FL×ΞwI,ϑδILXpϑX.Then, the multitarget prediction density at timek+1 is
(33)π+X+=ΔX+∑I+,ϑ∈FL+×Ξw+I+,ϑδI+LX+p+ϑX+,where
(34)w+I+,ϑ=wγI+∩BwSϑI+∩L,wsϑL=ηSϑL∑I⊇L1−ηSϑI−LwI,ϑ,ηSϑl=pS⋅,l,pϑ⋅,l,p+ϑx,l=1Llp+,Sϑx,l+1Blpγx,l,p+,Sϑx,l=pS⋅,lϕx⋅,l,pϑ⋅,lηSϑl.Here,wγI+∩B is the weight of newborn label I+∩B, and wSϑI+∩L is the weight of survival label. ϕx⋅,l is the single-target transfer density of trajectory l, and ηSϑl is the survival probability of trajectory l. pγx,l is the probability density of newborn state, and p+,Sϑx,l is the probability density of survival state.


### 3.2.2. Update

After obtaining the multitarget prediction density at timek+1, it can be proved that the updated multitarget filter density is also δ-GLMB.The update process is shown as follows:(35)πXZ=ΔX∑I,ϑ∈FL×Ξ∑θ∈ΘIwI,ϑ,θZδILXpϑ,θ⋅ZX,where
(36)wI,ϑ,θZ=δθ−10:ZIηZϑ,θIwI,ϑ∑I,ϑ∈FL×Ξ∑θ∈ΘIδθ−10:ZIηZϑ,θIwI,ϑ∝ηZϑ,θIwI,ϑ,ηZϑ,θl=pϑ⋅,l,φZ⋅,l;θ,pϑ,θx,lZ=pϑx,lφZx,l;θηZϑ,θl.Here,ΘI represents a subset of the current association mapping with domain I. φZx,l;θ is defined as follows:
(37)φZx,l;θ=δ0θl1−pDx,l+1−δ0θlpDx,lgzθlx,lκzθl=pDx,lgzθlx,lκzθl,θl>0,1−pDx,l,θl=0.Here,pDx,l is the detection probability on x,l. gzθlx,l is the single-target likelihood function of zθl to x,l. κ⋅ is the Poisson clutter density.


### 3.2.3. Realization ofδ-GLMB Filter

For cluster member motion model,pSx,l=pS, and ϕx+x,l=Nx+;Fkx,Qcs. Since the newborn density matrix is in the form of Gaussian mixture, if the single-target density is in the form of Gaussian mixture, there is
(38)pϑ⋅,l=∑i=1JϑlωiϑlNx;miϑl,Piθl.Then,(39)ηSϑl=pS,p+ϑx,l=1Ll∑i=1JϑlωiϑlNx;mS,iϑl,PS,iϑl+1Blpγl,where
(40)mS.iϑl=Fkmiϑl,PS,iϑl=FkPiϑlFkT+Qcs.When the motion model parameters are related to the label, we just need to substitutepS=pSl, Fk=Fkl, and Qcs=Qcsl into the above equation.In the update process,pDx,l=pD. For the measurement model shown in (30), the measurement likelihood is gzx,l=Nz;Hk+1x,Rk+1. Gaussian mixture expression provides the most general configuration for linear Gaussian model. It is assumed that each single-target density pϑ⋅,l is in the form of Gaussian mixture.
(41)pϑ⋅,l=∑n=1JϑlωnϑlNx;mnϑl,Pnϑl.Then,(42)ci,j=−lnpD∑n=1Jϑliωnϑliqnϑzj;li1−pDκzj.In addition, for the updated association processϑ,θ, there are
(43)ηZϑ,θl=∑n=1JϑlwZ,nϑ,θl,pϑ,θx,lZ=∑n=1JϑlwZ,nϑ,θlηZϑ,θlNx;mZ,nϑ,θl,Pnϑ,θl.Here,(44)wZ,nϑ,θl=ωnϑl×pDqnϑzθl;lκzθl,θl>0,1−pD,θl=0,qnϑz;l=Nz;Hk+1mnϑl,Hk+1PnϑlHk+1T+Rk+1,mZ,nϑ,θl=mnϑl+Gnϑ,θlzθl−Hk+1mnϑl,θl>0,mnϑl,θl=0,Pnϑ,θl=I−Gnϑ,θlHk+1Pnϑl,Gnϑ,θl=PnϑlHk+1THk+1PnϑlHk+1T+Rk+1−1,θl>0,0,θl=0.When the parameters of the measurement model are dependent on the labell, it is only necessary to substitute pD=pDl, Hk+1=Hk+1l, and Rk+1=Rk+1l into the above equation.


## 3.2.1. Prediction

At timek, the distribution of the random finite set of δ-GLMB is as follows:
(32)πX=ΔX∑I,ϑ∈FL×ΞwI,ϑδILXpϑX.Then, the multitarget prediction density at timek+1 is
(33)π+X+=ΔX+∑I+,ϑ∈FL+×Ξw+I+,ϑδI+LX+p+ϑX+,where
(34)w+I+,ϑ=wγI+∩BwSϑI+∩L,wsϑL=ηSϑL∑I⊇L1−ηSϑI−LwI,ϑ,ηSϑl=pS⋅,l,pϑ⋅,l,p+ϑx,l=1Llp+,Sϑx,l+1Blpγx,l,p+,Sϑx,l=pS⋅,lϕx⋅,l,pϑ⋅,lηSϑl.Here,wγI+∩B is the weight of newborn label I+∩B, and wSϑI+∩L is the weight of survival label. ϕx⋅,l is the single-target transfer density of trajectory l, and ηSϑl is the survival probability of trajectory l. pγx,l is the probability density of newborn state, and p+,Sϑx,l is the probability density of survival state.


## 3.2.2. Update

After obtaining the multitarget prediction density at timek+1, it can be proved that the updated multitarget filter density is also δ-GLMB.The update process is shown as follows:(35)πXZ=ΔX∑I,ϑ∈FL×Ξ∑θ∈ΘIwI,ϑ,θZδILXpϑ,θ⋅ZX,where
(36)wI,ϑ,θZ=δθ−10:ZIηZϑ,θIwI,ϑ∑I,ϑ∈FL×Ξ∑θ∈ΘIδθ−10:ZIηZϑ,θIwI,ϑ∝ηZϑ,θIwI,ϑ,ηZϑ,θl=pϑ⋅,l,φZ⋅,l;θ,pϑ,θx,lZ=pϑx,lφZx,l;θηZϑ,θl.Here,ΘI represents a subset of the current association mapping with domain I. φZx,l;θ is defined as follows:
(37)φZx,l;θ=δ0θl1−pDx,l+1−δ0θlpDx,lgzθlx,lκzθl=pDx,lgzθlx,lκzθl,θl>0,1−pDx,l,θl=0.Here,pDx,l is the detection probability on x,l. gzθlx,l is the single-target likelihood function of zθl to x,l. κ⋅ is the Poisson clutter density.


## 3.2.3. Realization ofδ-GLMB Filter

For cluster member motion model,pSx,l=pS, and ϕx+x,l=Nx+;Fkx,Qcs. Since the newborn density matrix is in the form of Gaussian mixture, if the single-target density is in the form of Gaussian mixture, there is
(38)pϑ⋅,l=∑i=1JϑlωiϑlNx;miϑl,Piθl.Then,(39)ηSϑl=pS,p+ϑx,l=1Ll∑i=1JϑlωiϑlNx;mS,iϑl,PS,iϑl+1Blpγl,where
(40)mS.iϑl=Fkmiϑl,PS,iϑl=FkPiϑlFkT+Qcs.When the motion model parameters are related to the label, we just need to substitutepS=pSl, Fk=Fkl, and Qcs=Qcsl into the above equation.In the update process,pDx,l=pD. For the measurement model shown in (30), the measurement likelihood is gzx,l=Nz;Hk+1x,Rk+1. Gaussian mixture expression provides the most general configuration for linear Gaussian model. It is assumed that each single-target density pϑ⋅,l is in the form of Gaussian mixture.
(41)pϑ⋅,l=∑n=1JϑlωnϑlNx;mnϑl,Pnϑl.Then,(42)ci,j=−lnpD∑n=1Jϑliωnϑliqnϑzj;li1−pDκzj.In addition, for the updated association processϑ,θ, there are
(43)ηZϑ,θl=∑n=1JϑlwZ,nϑ,θl,pϑ,θx,lZ=∑n=1JϑlwZ,nϑ,θlηZϑ,θlNx;mZ,nϑ,θl,Pnϑ,θl.Here,(44)wZ,nϑ,θl=ωnϑl×pDqnϑzθl;lκzθl,θl>0,1−pD,θl=0,qnϑz;l=Nz;Hk+1mnϑl,Hk+1PnϑlHk+1T+Rk+1,mZ,nϑ,θl=mnϑl+Gnϑ,θlzθl−Hk+1mnϑl,θl>0,mnϑl,θl=0,Pnϑ,θl=I−Gnϑ,θlHk+1Pnϑl,Gnϑ,θl=PnϑlHk+1THk+1PnϑlHk+1T+Rk+1−1,θl>0,0,θl=0.When the parameters of the measurement model are dependent on the labell, it is only necessary to substitute pD=pDl, Hk+1=Hk+1l, and Rk+1=Rk+1l into the above equation.


## 4. Multiradar Estimation Fusion Based on Graph-LSTMs

Through the above process, we can obtain multiple estimations of the same member state. An obvious fact is that in the tracking process, we only consider the explicit relationship between radar and target in detection and tracking but hardly consider the implicit relationship of motion consistency among cluster members which may play a positive role in improving the tracking accuracy. In this section, the purpose of applying the Graph-LSTM algorithm is to capture this implicit relationship.


### 4.1. Construction Tracking Relationship Graph

#### 4.1.1. Bipartite Graph Model

In order to reflect the detection and tracking relationship between radar and cluster members, we establish the graph of radar and cluster. Radar system and cluster can be divided into two subsetsS,T, and we connect the radar si that has detection relationship and obtains the target estimation with the corresponding target tj by edge E⊆S×T, and Figure 4 is as follows. According to the definition of bipartite graph, this graph is a bipartite graph, expressed as G=S,T,E.Figure 4
Bipartite graph of radar-target.


#### 4.1.2. Build Node Unfolding Tree

In the bipartite graph of detection and tracking relationship between radar and target, the size and topology of the graph neighborhood are variable. However, the common machine learning algorithms can only deal with feature vectors with fixed length. In order to apply the machine learning algorithm to mining information, an effective solution is to unfold the bipartite graph into a tree at a given depth and the learning algorithm is used to learn the sequence from the child node to the root node.For a bipartite graph, given the depthD=2, the tree Ttj is constructed with the target tj as the root node, and the unfolding results are shown in Figure 5. In the unfolding tree, the root node depth is 0. The radars with detection and tracking relationship with target tj are taken as child nodes, and the depth is 1. Other detection and tracking targets of these radars are the child nodes of the radar, with a depth of 2. In the process of building Ttj, when the depth d>1, the same target may be repeatedly detected by other radars. At this time, for the same node k, we would add superscript to represent its copy with k′ and k″. From the depth d≥2, there is an upward transmission of other target motion information in the tree, which makes the acquisition of fusion apply the implicit information between targets.Figure 5
An example of a graph and its unfolding tree at nodetj for depth 2.


### 4.2. State Space Labeling

In this paper, we would apply graph-long short-term memory neural nets (Graph-LSTMs) to sequence learning. The core of Graph-LSTMs is multilevel LSTMs, which is called multilevel sequence learners (MLSL). However, MLSL can only output the label of the target sate, and the continuous state of the target is difficult to adapt to this output form. In order to apply MLSL, it is necessary to label the state space and convert the true state and multiple estimations of the target into label space.The labeling of state space needs to be carried out in position space and velocity space, respectively. In position space, given the number of labelsNL and taking the predicted position as the center, we can generate an orthogonal grid of NL−1×NL−1, as shown in Figure 6. The space around the predicted position is discretized into NL−12 small spaces. We call the small space as label space. We use Lp to represent its label, and Lp=Lx,LyT. Similarly, we can get label Lv in velocity space and Lv=Lvx,LvyT. Then, we can obtain the label of the target estimation L=Lx,Lvx,Ly,LvyT.Figure 6
Position space labeling.At timek, the estimated state of target tj is xkk,j. At time k+1, the label Lk+1,ji of radar si estimation xk+1k+1,ji for target tj is calculated as follows:
(45)Lk+1,ji=xk+1k+1,ji−xk+1k,j⊙dΔL−1+Lk+1k,j,where xk+1k,j is the state prediction at time k+1; Lk+1,ji represents the label of the estimation xk+1k+1,ji of radar si; Lk+1k,j represents the label of xk+1k,j, which is usually used as the central label of label space; ⊙ represents the product of the corresponding element; and dΔL represents the side length of label space, which can reflect the fineness of state space labeling. In the paper, dΔL=dΔLx,dΔLvx,dΔLy,dΔLvyT. The elements of dΔL are preset by us. ⋅ means rounding down. Therefore, the smallest label is 0. After the label space is given, if the calculated target estimation label is outside the label space, the estimation would be discarded.


### 4.3. Sequence Learning Using MLSL

#### 4.3.1. Learning Strategy

Next, we would traverse the treeTtj from bottom to top; use the MLSL to calculate the state label ytj of root node. In the process of network training, ytj is compared with the expected state label y¯tj, the loss is calculated, and backpropagation is carried out to train the network parameters.


#### 4.3.2. MLSL Model

LSTM is the basis for efficient sequence learning in MLSL, which can alleviate the problem of gradient disappearance with the increase of time step in recurrent neural network. The basis of sequence learning ability of LSTM is the introduction of three gates in its basic unit. The structure and operation flow of LSTM unit are shown in Figure7.Figure 7
LSTMs cell.In the picture,id is the input gate, which controls the information of the input vector xd and updates the current information flow Cd. fd is the forget gate, which determines whether the value information is retained in Cd. od is the output gate, which includes the input information of level d+1 and determines the value of hidden state hd of level d. σ and tanh are the activation functions. The operation update mechanism of each gate is as follows:
(46)fd=σWf×hd+1,xd+bf,id=σWi×hd+1,xd+bi,C^d=tanhWC×hd+1,xd+bC,Cd=fd×Cd+1+id×C^d,od=σWo×hd+1,xd+bo,hd=od×tanhCd,where Wf, Wi, WC, and Wo are the weight parameters of each gate obtained by network training and bf, bi, bC, and bo are the bias obtained by network training.LSTMs composed of basic cells in series can accept any length of eigenvector as input and generate an output vector. If an LSTM can accept a sequence of vectors containingN vectors and generate an output vector of size K, we call its shape N,K. After obtaining the unfolding tree, D kinds of LSTM cell are applied to the tree with depth D, and the number of each kind of LSTM cell is the same as the number of nodes in its depth. At its depth, each LSTM cell processes the relevant node information to obtain the output vector.For LSTMs with depthd0<d<D, it is expressed as Ld and the shape is M+Kd+1,Kd. M indicates that there are M edges connected to the child nodes. Kd+1 indicates the output from the lower level LSTMs, and M+Kd+1 indicates the series connection of the two eigenvectors. Taking Figure 5 as an example, the output ftj of the target tj is calculated as follows:
(47)ftj=fsa∩gsa,tj+fsb∩gsb,tj+fsc∩gsc,tj,fsa=gsa,tk+gsa,tl,fsb=gsb,tk,fsc=gsc,tk+gsc,tl.Here,f⋅ and g⋅ are the input and output vectors of LSTMs, respectively. gsi,tj is the state estimation label of target tj obtained by radar si. The symbol ∩ represents the concatenation of feature vectors, and the symbol + represents the collection of information.


#### 4.3.3. Parameter Update Process

Ld with the same depth have the same parameters, which are expressed as parameter vector Wd. When updating, nodes in the same level independently forward and backward propagation to minimize the loss function. In the training process, repeatedly select the same root node tj∗, unfold graph Gtj∗ into tree Ttj∗, and carry out forward and backward propagation to update parameters.Forward propagation is carried out from the lowest level child node to the root node. For nodeu with arbitrary depth d, if it has d+1 level LSTMs, it obtains the output fu from d+1 level. The information that node u propagates to the parent node v is gv,u∩fu.In backward propagation, after obtaining the output eigenvectorytj=ftj, the loss Lytj is calculated through the loss function.
(48)Lytj=ytj−y¯tj.Here,y¯tj is the expected output eigenvector. In a training, after the loss is obtained, the derivative ∂L/∂ytj of the loss to the output vector can be obtained. Then, the derivative of the loss to any input vector is obtained. When the input vector is like gv,u∩fu, since gv,u does not need to propagate, the backward propagation of loss only needs to calculate ∂L/∂fu. The derivative of the loss propagates from the root node to the lowest level child node to stop when the depth is D and the input vector is only gv,u.Since there are manyLd at the same depth and have the same parameters, the Ld at the same depth need to calculate and update the parameters in parallel and then fuse the parameters for each Ld. Suppose that there are mLd in one level, which are labeled as Ld1, Ld2,…,Ldm. If the updated parameter ΔiWd can be obtained for Ldi, the updated parameter ΔWd is calculated as follows:
(49)ΔWd=Δ1Wd+Δ2Wd+⋯+ΔmWdm.


#### 4.3.4. Overall Network Structure and Output

Taking depthD=2 as an example, the estimation fusion network is shown in Figure 8. Firstly, a bipartite graph is constructed according to the detection and tracking relationship between radar and target, and we unfold the graph into a tree with the tracked target as the root node. Secondly, the target state is labeled, and the labels of the estimation are obtained from the velocity and position space, respectively. Then, we input the estimation label into the trained MLSL. Finally, we input the network output vector into the Softmax full connection layer to obtain the probability of each label in the label space. The state label fusion L^k+1,j of the target is the label with the highest probability in the label space.Figure 8
Overall network structure.After obtaining the target label fusionL^k+1,j, we need to convert the label fusion of the target into the state fusion. The transformation method is as follows:
(50)xk+1k+1,j=L^k+1,j−Lk+1k,j⊙dΔL+xk+1k,j.


## 4.1. Construction Tracking Relationship Graph

### 4.1.1. Bipartite Graph Model

In order to reflect the detection and tracking relationship between radar and cluster members, we establish the graph of radar and cluster. Radar system and cluster can be divided into two subsetsS,T, and we connect the radar si that has detection relationship and obtains the target estimation with the corresponding target tj by edge E⊆S×T, and Figure 4 is as follows. According to the definition of bipartite graph, this graph is a bipartite graph, expressed as G=S,T,E.Figure 4
Bipartite graph of radar-target.


### 4.1.2. Build Node Unfolding Tree

In the bipartite graph of detection and tracking relationship between radar and target, the size and topology of the graph neighborhood are variable. However, the common machine learning algorithms can only deal with feature vectors with fixed length. In order to apply the machine learning algorithm to mining information, an effective solution is to unfold the bipartite graph into a tree at a given depth and the learning algorithm is used to learn the sequence from the child node to the root node.For a bipartite graph, given the depthD=2, the tree Ttj is constructed with the target tj as the root node, and the unfolding results are shown in Figure 5. In the unfolding tree, the root node depth is 0. The radars with detection and tracking relationship with target tj are taken as child nodes, and the depth is 1. Other detection and tracking targets of these radars are the child nodes of the radar, with a depth of 2. In the process of building Ttj, when the depth d>1, the same target may be repeatedly detected by other radars. At this time, for the same node k, we would add superscript to represent its copy with k′ and k″. From the depth d≥2, there is an upward transmission of other target motion information in the tree, which makes the acquisition of fusion apply the implicit information between targets.Figure 5
An example of a graph and its unfolding tree at nodetj for depth 2.


## 4.1.1. Bipartite Graph Model

In order to reflect the detection and tracking relationship between radar and cluster members, we establish the graph of radar and cluster. Radar system and cluster can be divided into two subsetsS,T, and we connect the radar si that has detection relationship and obtains the target estimation with the corresponding target tj by edge E⊆S×T, and Figure 4 is as follows. According to the definition of bipartite graph, this graph is a bipartite graph, expressed as G=S,T,E.Figure 4
Bipartite graph of radar-target.


## 4.1.2. Build Node Unfolding Tree

In the bipartite graph of detection and tracking relationship between radar and target, the size and topology of the graph neighborhood are variable. However, the common machine learning algorithms can only deal with feature vectors with fixed length. In order to apply the machine learning algorithm to mining information, an effective solution is to unfold the bipartite graph into a tree at a given depth and the learning algorithm is used to learn the sequence from the child node to the root node.For a bipartite graph, given the depthD=2, the tree Ttj is constructed with the target tj as the root node, and the unfolding results are shown in Figure 5. In the unfolding tree, the root node depth is 0. The radars with detection and tracking relationship with target tj are taken as child nodes, and the depth is 1. Other detection and tracking targets of these radars are the child nodes of the radar, with a depth of 2. In the process of building Ttj, when the depth d>1, the same target may be repeatedly detected by other radars. At this time, for the same node k, we would add superscript to represent its copy with k′ and k″. From the depth d≥2, there is an upward transmission of other target motion information in the tree, which makes the acquisition of fusion apply the implicit information between targets.Figure 5
An example of a graph and its unfolding tree at nodetj for depth 2.


## 4.2. State Space Labeling

In this paper, we would apply graph-long short-term memory neural nets (Graph-LSTMs) to sequence learning. The core of Graph-LSTMs is multilevel LSTMs, which is called multilevel sequence learners (MLSL). However, MLSL can only output the label of the target sate, and the continuous state of the target is difficult to adapt to this output form. In order to apply MLSL, it is necessary to label the state space and convert the true state and multiple estimations of the target into label space.The labeling of state space needs to be carried out in position space and velocity space, respectively. In position space, given the number of labelsNL and taking the predicted position as the center, we can generate an orthogonal grid of NL−1×NL−1, as shown in Figure 6. The space around the predicted position is discretized into NL−12 small spaces. We call the small space as label space. We use Lp to represent its label, and Lp=Lx,LyT. Similarly, we can get label Lv in velocity space and Lv=Lvx,LvyT. Then, we can obtain the label of the target estimation L=Lx,Lvx,Ly,LvyT.Figure 6
Position space labeling.At timek, the estimated state of target tj is xkk,j. At time k+1, the label Lk+1,ji of radar si estimation xk+1k+1,ji for target tj is calculated as follows:
(45)Lk+1,ji=xk+1k+1,ji−xk+1k,j⊙dΔL−1+Lk+1k,j,where xk+1k,j is the state prediction at time k+1; Lk+1,ji represents the label of the estimation xk+1k+1,ji of radar si; Lk+1k,j represents the label of xk+1k,j, which is usually used as the central label of label space; ⊙ represents the product of the corresponding element; and dΔL represents the side length of label space, which can reflect the fineness of state space labeling. In the paper, dΔL=dΔLx,dΔLvx,dΔLy,dΔLvyT. The elements of dΔL are preset by us. ⋅ means rounding down. Therefore, the smallest label is 0. After the label space is given, if the calculated target estimation label is outside the label space, the estimation would be discarded.


## 4.3. Sequence Learning Using MLSL

### 4.3.1. Learning Strategy

Next, we would traverse the treeTtj from bottom to top; use the MLSL to calculate the state label ytj of root node. In the process of network training, ytj is compared with the expected state label y¯tj, the loss is calculated, and backpropagation is carried out to train the network parameters.


### 4.3.2. MLSL Model

LSTM is the basis for efficient sequence learning in MLSL, which can alleviate the problem of gradient disappearance with the increase of time step in recurrent neural network. The basis of sequence learning ability of LSTM is the introduction of three gates in its basic unit. The structure and operation flow of LSTM unit are shown in Figure7.Figure 7
LSTMs cell.In the picture,id is the input gate, which controls the information of the input vector xd and updates the current information flow Cd. fd is the forget gate, which determines whether the value information is retained in Cd. od is the output gate, which includes the input information of level d+1 and determines the value of hidden state hd of level d. σ and tanh are the activation functions. The operation update mechanism of each gate is as follows:
(46)fd=σWf×hd+1,xd+bf,id=σWi×hd+1,xd+bi,C^d=tanhWC×hd+1,xd+bC,Cd=fd×Cd+1+id×C^d,od=σWo×hd+1,xd+bo,hd=od×tanhCd,where Wf, Wi, WC, and Wo are the weight parameters of each gate obtained by network training and bf, bi, bC, and bo are the bias obtained by network training.LSTMs composed of basic cells in series can accept any length of eigenvector as input and generate an output vector. If an LSTM can accept a sequence of vectors containingN vectors and generate an output vector of size K, we call its shape N,K. After obtaining the unfolding tree, D kinds of LSTM cell are applied to the tree with depth D, and the number of each kind of LSTM cell is the same as the number of nodes in its depth. At its depth, each LSTM cell processes the relevant node information to obtain the output vector.For LSTMs with depthd0<d<D, it is expressed as Ld and the shape is M+Kd+1,Kd. M indicates that there are M edges connected to the child nodes. Kd+1 indicates the output from the lower level LSTMs, and M+Kd+1 indicates the series connection of the two eigenvectors. Taking Figure 5 as an example, the output ftj of the target tj is calculated as follows:
(47)ftj=fsa∩gsa,tj+fsb∩gsb,tj+fsc∩gsc,tj,fsa=gsa,tk+gsa,tl,fsb=gsb,tk,fsc=gsc,tk+gsc,tl.Here,f⋅ and g⋅ are the input and output vectors of LSTMs, respectively. gsi,tj is the state estimation label of target tj obtained by radar si. The symbol ∩ represents the concatenation of feature vectors, and the symbol + represents the collection of information.


### 4.3.3. Parameter Update Process

Ld with the same depth have the same parameters, which are expressed as parameter vector Wd. When updating, nodes in the same level independently forward and backward propagation to minimize the loss function. In the training process, repeatedly select the same root node tj∗, unfold graph Gtj∗ into tree Ttj∗, and carry out forward and backward propagation to update parameters.Forward propagation is carried out from the lowest level child node to the root node. For nodeu with arbitrary depth d, if it has d+1 level LSTMs, it obtains the output fu from d+1 level. The information that node u propagates to the parent node v is gv,u∩fu.In backward propagation, after obtaining the output eigenvectorytj=ftj, the loss Lytj is calculated through the loss function.
(48)Lytj=ytj−y¯tj.Here,y¯tj is the expected output eigenvector. In a training, after the loss is obtained, the derivative ∂L/∂ytj of the loss to the output vector can be obtained. Then, the derivative of the loss to any input vector is obtained. When the input vector is like gv,u∩fu, since gv,u does not need to propagate, the backward propagation of loss only needs to calculate ∂L/∂fu. The derivative of the loss propagates from the root node to the lowest level child node to stop when the depth is D and the input vector is only gv,u.Since there are manyLd at the same depth and have the same parameters, the Ld at the same depth need to calculate and update the parameters in parallel and then fuse the parameters for each Ld. Suppose that there are mLd in one level, which are labeled as Ld1, Ld2,…,Ldm. If the updated parameter ΔiWd can be obtained for Ldi, the updated parameter ΔWd is calculated as follows:
(49)ΔWd=Δ1Wd+Δ2Wd+⋯+ΔmWdm.


### 4.3.4. Overall Network Structure and Output

Taking depthD=2 as an example, the estimation fusion network is shown in Figure 8. Firstly, a bipartite graph is constructed according to the detection and tracking relationship between radar and target, and we unfold the graph into a tree with the tracked target as the root node. Secondly, the target state is labeled, and the labels of the estimation are obtained from the velocity and position space, respectively. Then, we input the estimation label into the trained MLSL. Finally, we input the network output vector into the Softmax full connection layer to obtain the probability of each label in the label space. The state label fusion L^k+1,j of the target is the label with the highest probability in the label space.Figure 8
Overall network structure.After obtaining the target label fusionL^k+1,j, we need to convert the label fusion of the target into the state fusion. The transformation method is as follows:
(50)xk+1k+1,j=L^k+1,j−Lk+1k,j⊙dΔL+xk+1k,j.


## 4.3.1. Learning Strategy

Next, we would traverse the treeTtj from bottom to top; use the MLSL to calculate the state label ytj of root node. In the process of network training, ytj is compared with the expected state label y¯tj, the loss is calculated, and backpropagation is carried out to train the network parameters.


## 4.3.2. MLSL Model

LSTM is the basis for efficient sequence learning in MLSL, which can alleviate the problem of gradient disappearance with the increase of time step in recurrent neural network. The basis of sequence learning ability of LSTM is the introduction of three gates in its basic unit. The structure and operation flow of LSTM unit are shown in Figure7.Figure 7
LSTMs cell.In the picture,id is the input gate, which controls the information of the input vector xd and updates the current information flow Cd. fd is the forget gate, which determines whether the value information is retained in Cd. od is the output gate, which includes the input information of level d+1 and determines the value of hidden state hd of level d. σ and tanh are the activation functions. The operation update mechanism of each gate is as follows:
(46)fd=σWf×hd+1,xd+bf,id=σWi×hd+1,xd+bi,C^d=tanhWC×hd+1,xd+bC,Cd=fd×Cd+1+id×C^d,od=σWo×hd+1,xd+bo,hd=od×tanhCd,where Wf, Wi, WC, and Wo are the weight parameters of each gate obtained by network training and bf, bi, bC, and bo are the bias obtained by network training.LSTMs composed of basic cells in series can accept any length of eigenvector as input and generate an output vector. If an LSTM can accept a sequence of vectors containingN vectors and generate an output vector of size K, we call its shape N,K. After obtaining the unfolding tree, D kinds of LSTM cell are applied to the tree with depth D, and the number of each kind of LSTM cell is the same as the number of nodes in its depth. At its depth, each LSTM cell processes the relevant node information to obtain the output vector.For LSTMs with depthd0<d<D, it is expressed as Ld and the shape is M+Kd+1,Kd. M indicates that there are M edges connected to the child nodes. Kd+1 indicates the output from the lower level LSTMs, and M+Kd+1 indicates the series connection of the two eigenvectors. Taking Figure 5 as an example, the output ftj of the target tj is calculated as follows:
(47)ftj=fsa∩gsa,tj+fsb∩gsb,tj+fsc∩gsc,tj,fsa=gsa,tk+gsa,tl,fsb=gsb,tk,fsc=gsc,tk+gsc,tl.Here,f⋅ and g⋅ are the input and output vectors of LSTMs, respectively. gsi,tj is the state estimation label of target tj obtained by radar si. The symbol ∩ represents the concatenation of feature vectors, and the symbol + represents the collection of information.


## 4.3.3. Parameter Update Process

Ld with the same depth have the same parameters, which are expressed as parameter vector Wd. When updating, nodes in the same level independently forward and backward propagation to minimize the loss function. In the training process, repeatedly select the same root node tj∗, unfold graph Gtj∗ into tree Ttj∗, and carry out forward and backward propagation to update parameters.Forward propagation is carried out from the lowest level child node to the root node. For nodeu with arbitrary depth d, if it has d+1 level LSTMs, it obtains the output fu from d+1 level. The information that node u propagates to the parent node v is gv,u∩fu.In backward propagation, after obtaining the output eigenvectorytj=ftj, the loss Lytj is calculated through the loss function.
(48)Lytj=ytj−y¯tj.Here,y¯tj is the expected output eigenvector. In a training, after the loss is obtained, the derivative ∂L/∂ytj of the loss to the output vector can be obtained. Then, the derivative of the loss to any input vector is obtained. When the input vector is like gv,u∩fu, since gv,u does not need to propagate, the backward propagation of loss only needs to calculate ∂L/∂fu. The derivative of the loss propagates from the root node to the lowest level child node to stop when the depth is D and the input vector is only gv,u.Since there are manyLd at the same depth and have the same parameters, the Ld at the same depth need to calculate and update the parameters in parallel and then fuse the parameters for each Ld. Suppose that there are mLd in one level, which are labeled as Ld1, Ld2,…,Ldm. If the updated parameter ΔiWd can be obtained for Ldi, the updated parameter ΔWd is calculated as follows:
(49)ΔWd=Δ1Wd+Δ2Wd+⋯+ΔmWdm.


## 4.3.4. Overall Network Structure and Output

Taking depthD=2 as an example, the estimation fusion network is shown in Figure 8. Firstly, a bipartite graph is constructed according to the detection and tracking relationship between radar and target, and we unfold the graph into a tree with the tracked target as the root node. Secondly, the target state is labeled, and the labels of the estimation are obtained from the velocity and position space, respectively. Then, we input the estimation label into the trained MLSL. Finally, we input the network output vector into the Softmax full connection layer to obtain the probability of each label in the label space. The state label fusion L^k+1,j of the target is the label with the highest probability in the label space.Figure 8
Overall network structure.After obtaining the target label fusionL^k+1,j, we need to convert the label fusion of the target into the state fusion. The transformation method is as follows:
(50)xk+1k+1,j=L^k+1,j−Lk+1k,j⊙dΔL+xk+1k,j.


## 5. Experiments and Results

### 5.1. Simulation Scene Setting

In order to evaluate the tracking performance of the algorithm proposed in this paper, we set up a simulation experiment scene in 2-dimensional space. In this simulation scenario, a total of 3000 simulation tracks are simulated. 30 targets form a cluster and move at the same time. Then, there are 100 clusters in total. Total time of cluster movementT=100s. 80% of the simulation data set is used for training and 20% for testing.Set the cluster motion model parameters as shown in Table1. 30 cluster members are randomly dropped in area 0,100×0,100 (m2), and the initial velocity in each direction is a random number in the range of [0,140] m/s. After the cluster is dropped, the members adjust the actual velocity through self-organizing cooperative interaction. Due to the cooperative interaction, the actual velocity of the member is not ẋ,ẏ in the state vector. If the sampling interval is τ, the calculation method of actual velocity v is shown in (51). The movement process of a cluster is shown in Figure 9, and the actual velocity of each member is shown in Figure 10.
(51)v=vxvy=1τxt+τyt+τ−xtyt.Table 1
Cluster motion model parameters.
ParameterValueParameterValuedb20dm50α0.1β0.2σx1σy1σg1R1140R123R2140R221Figure 9
A cluster movement process.Figure 10
Actual velocity of each member.As can be seen from Figures9 and 10, after the cluster falls, the members begin to diffuse under the action of repulsive force due to the small spacing between each member. Affected by the cooperative interaction, the actual velocity of each member begins to converge. vx converges to about 60 m/s, and vy converges to about 75 m/s in about 30 s. After the velocity of cluster converges, the cluster forms a stable structure and moves at almost the same velocity. At this time, the dispersion range of the cluster is about 1600×1600m2.Figure9 also shows that after the cluster forms a cooperative motion with almost the same velocity, the members in the cluster still have dynamic adjustment due to the action of motion noise and potential force. This adjustment action is very common, which would have a negative impact on the stable tracking of δ-GLMB filter.


### 5.2.δ-GLMB Filter Tracking Results

When the number of cluster membersN=30, given the minimum value of the average detection density A¯d=2.5, select the number of radars m=10 and set the maximum number of radar tracking targets Nc=8. In this case, the average detection density Ad≈2.7>A¯d. According to the above parameters, a detection confirmation matrix Ω is determined. The graphical representation of the matrix Ω is shown in Figure 11, in which 30 columns of the grid represent 30 targets and 10 rows represent 10 radars. Those marked with color in the grid are regarded as having detection relationship.Figure 11
Graphical representation of matrixΩ.Set the initial state covariance of cluster targetsP0=diag1×104,1×103,1×104,1×103T, the measurement noise covariance Rk=diag4×102,4×102T, and the potential force noise Sk=diag5,2.5,2.5,5T×10−4. The clutter obeys the Poisson distribution of the expected value λc=8 in time.Theδ-GLMB algorithm parameters are set as follows: target detection probability pD=0.98 and target survival probability ps=0.99. The newborn target model is determined according to the measurement set at time k, i.e., πγ=rγi,pγii=1Zk, where the newborn probability rγi=0.02 and state distribution of newborn targets pγix=Nx;mgi,Pg. Here, mγi=zk,i,0T and Pγ=diag1×104,5×103,1×104,5×103T.In the process of tracking 100 clusters,δ-GLMB filter runs 1000 times in total. The main filtering conditions are shown in Figure 12, and the proportion of each condition is shown in Table 2.Figure 12
δ-GLMB filtering condition.
(a)(b)(c)(d)Table 2
Proportion of filtering condition.
ConditionProportionConditionProportionFigure12(a)57.6%Figure12(b)28.3%Figure12(c)5.7%Figure12(d)7.1%Others1.3%As shown in Figure12, under the set parameters, the δ-GLMB algorithm can track all 8 assigned cluster members with a proportion of 91.6%. Figure 12(a) shows that the filter can accurately track all targets in a short time and achieve stable tracking in the whole process of target movement. This condition is the most ideal, accounting for 57.6% of all major conditions. Figure 12(b) shows that the filter can achieve stable tracking of 8 targets after 20 s. This is due to the slow acquisition of the target by the filter because of the wide distribution of the target, which is unfavorable to the accurate tracking of the target in the early stage of motion. Figure 12(c) reflects that in the later stage of tracking, the target filtering results of some filters have the problems of inaccurate track association and filtering tracks crossing, which is related to the close spacing of the assigned targets. The tracking condition shown in Figure 12(d) is not ideal. The filter has been unable to track the 2 assigned targets, accounting for 7.1%. However, we need not worry too much about this, because other filters may obtain the motion estimation of these targets and make up for these defects in the fusion stage. Other tracking conditions include discontinuous tracks, multiple track initiation, and other problems, accounting for a small proportion, which is caused by the defects of the random finite set tracking algorithm itself.Accumulate the number of targets tracked by 10 radars to obtain the total number of targets that can be tracked by the radar system, and calculate the average value of the total number of 100 experiments. The results are shown in Figure13.Figure 13
Total number of targets tracked.It can be seen from Figure13 that the radar system can only track about 20 targets at the initial time when Ad≈2.7. Then, the total number of targets tracked increases sharply to about 28. After 15 s, almost all targets can be tracked stably.


### 5.3. Estimation Fusion Results Using Graph-LSTMs

#### 5.3.1. Network Training Results

In order to verify the effectiveness of Graph-LSTMs for multiradar estimation fusion, set the label space lengthdΔLx=dΔLy=2 m and dΔLvx=dΔLvy=0.5 m/s. Set the number of labels NL=11 on each dimension, that is, 0,1,2,⋯,10. The training set is used to train the network parameters, and the test set is used to test the fusion effect. The results are shown in Figures 14 and 15.Figure 14
Training effect of 1-MLSL.
(a)(b)Figure 15
Training effect of different depth networks.Figure14 shows the label prediction accuracy of 1-MLSL with different training times. 1-MLSL indicates MLSL with depth D=1. Figure 14(b) shows the mean and standard deviation of the probability of correct label output from Softmax layer. It can be seen from Figure 14(a) that the more training times, the better the network fusion effect. When the training times are less than 10000, the prediction accuracy of target labels increases rapidly with the increase of training times. When the number of training times is more than 10000, the prediction accuracy increases slowly, and the network parameter training is completed. Figure 14(b) shows that after the accuracy of label prediction increases to a certain level, the prediction probability of correct label can only increase to about 0.38. However, with the increase of training times, the standard deviation decreases, and the network prediction is more stable.The curves shown in Figure15 are the fitting result of the prediction accuracy of the target label under different depths and different training times. Figure 14 shows that when the training times are less, the network prediction effect with depth D=1 is better. When the number of training times increases, the prediction effect of the network with depth D=1 improves significantly and exceeds the effect of D=1. The effect of depth D=1 is always slightly worse than that of D=1,2. This shows that when the number of training is less, the simple network that only considers the detection and tracking relationship between radar and target can obtain better fusion effect, and with the increase of training times, the complex network that can mine the cooperative relationship between targets can play a better prediction effect. When the network depth D=3, the relationship between radars is considered, which may have a better effect on mining deeper information, but it would take a lot of training time.


#### 5.3.2. Estimation Fusion Effect

In this experiment, we compare the prediction probability of the target true label by using the Graph-LSTM algorithm proposed in this paper, the commonly used EM algorithm [40, 41] and the mean method, and give the corresponding OSPA distance in target tracking. The Graph-LSTM algorithm sets the depth D=2. The OSPA distance is used to evaluate the tracking performance of the algorithm in the experiment. The OSPA distance is defined as follows [42].At timek, the true state set of cluster members is Xk=xk,1,xk,2,⋯,xk,m, and the target estimation state set given by the estimation algorithm is Xkk=xkk,1,xkk,2,⋯,xkk,n, where m and n represent the actual number of targets and the estimated number of targets in the cluster, respectively. xk,i∈Xk and xkk,i∈Xkk represent the actual value and estimated value of the target i state vector, respectively. The OSPA distance between the true state set and the estimated state set is calculated as follows:
(52)OSPAp,cXk,Xkk=1nminπ∈Πn∑i=1mdcxk,i,xkk,πip+n−m⋅cp1/p,m≤n,OSPAp,cXkk,Xk,m>n,where Πn=pnm represents all permutations of m elements in the state set.
(53)dcxk,i,xkk,πi=minc,xk,i−xkk,πi.minπ∈Πn represents the group with the smallest distance difference between the true point trace and the estimated point trace of all targets. c and p are distance sensitivity parameters and correlation sensitivity parameters, respectively.Setc=100 and p=1 in the experiment; the network fusion effect and OSPA distance are shown in Figure 16.Figure 16
Network fusion effect and OSPA distance.Figure16 shows the prediction accuracy of target state label and the average OSPA distance of 20 test sets in tracking. It can be seen from Figure 16 that at the initial time, the prediction accuracy of target state label is low and the corresponding OSPA distance is large because the radar system is difficult to track all targets. As the filtering process proceeds, all targets are tracked and the OSPA distance decreases. At about 80 s, the OSPA distance increases due to correlation errors in some filters.In the whole tracking process, the Graph-LSTM algorithm used in this paper has better prediction accuracy of correct labels than the comparison methods, which is related to the fact that the Graph-LSTM algorithm with depthD=2 can consider the cooperative interaction between cluster targets. Accordingly, the method of using Graph-LSTM algorithm to calculate the fusion has less tracking error, and it can better alleviate the problem of inaccurate target tracking caused by correlation error in the later stage of tracking.


### 5.4. More Simulation Samples

In order to reflect the improving effect of our method on the accuracy of obtaining multiradar fusion estimates, we set more cluster motion trajectory data and obtain more fusion results of target estimators. There are four starting points of cluster motion in the sensor’s field of view, namely,±8×103,±8×103, and the parameters of the cluster motion model are consistent with those above. The starting point and direction of cluster motion are random. We obtain 20 cluster motion trajectories. Cluster center trajectories are shown in Figure 17.Figure 17
Cluster center trajectories.Six detection allocation schemes are randomly generated for each cluster, and a set of detection schemes is shown in Figure18. We can obtain 120 sets of multiradar estimates. The proposed method, mean method, and EM algorithm are, respectively, used to fuse multiple estimates of the target. The average OSPA distance of the fusions is shown in Figure 19.Figure 18
A set of cluster detection schemes.Figure 19
The average OSPA distance of the fusions.According to Figure19, the average OSPA distance medians of the proposed method, EM algorithm, and mean method are 22.4, 27.8, and 34.3, respectively. Our method is superior to the mean method and EM algorithm in terms of fusion accuracy of estimates, which reflects the effectiveness and superiority of our method for multiple estimate fusion. However, the effect of the EM algorithm is more stable. This is because there are large differences between some test data and training data in the whole estimation process, and the Graph-LSTM algorithm parameters cannot be accurately obtained by learning, which makes our method not effective. In order to further improve the robustness of the proposed algorithm, the generalization ability of the algorithm can be improved by increasing the number of training samples or pretraining.


## 5.1. Simulation Scene Setting

In order to evaluate the tracking performance of the algorithm proposed in this paper, we set up a simulation experiment scene in 2-dimensional space. In this simulation scenario, a total of 3000 simulation tracks are simulated. 30 targets form a cluster and move at the same time. Then, there are 100 clusters in total. Total time of cluster movementT=100s. 80% of the simulation data set is used for training and 20% for testing.Set the cluster motion model parameters as shown in Table1. 30 cluster members are randomly dropped in area 0,100×0,100 (m2), and the initial velocity in each direction is a random number in the range of [0,140] m/s. After the cluster is dropped, the members adjust the actual velocity through self-organizing cooperative interaction. Due to the cooperative interaction, the actual velocity of the member is not ẋ,ẏ in the state vector. If the sampling interval is τ, the calculation method of actual velocity v is shown in (51). The movement process of a cluster is shown in Figure 9, and the actual velocity of each member is shown in Figure 10.
(51)v=vxvy=1τxt+τyt+τ−xtyt.Table 1
Cluster motion model parameters.
ParameterValueParameterValuedb20dm50α0.1β0.2σx1σy1σg1R1140R123R2140R221Figure 9
A cluster movement process.Figure 10
Actual velocity of each member.As can be seen from Figures9 and 10, after the cluster falls, the members begin to diffuse under the action of repulsive force due to the small spacing between each member. Affected by the cooperative interaction, the actual velocity of each member begins to converge. vx converges to about 60 m/s, and vy converges to about 75 m/s in about 30 s. After the velocity of cluster converges, the cluster forms a stable structure and moves at almost the same velocity. At this time, the dispersion range of the cluster is about 1600×1600m2.Figure9 also shows that after the cluster forms a cooperative motion with almost the same velocity, the members in the cluster still have dynamic adjustment due to the action of motion noise and potential force. This adjustment action is very common, which would have a negative impact on the stable tracking of δ-GLMB filter.


## 5.2.δ-GLMB Filter Tracking Results

When the number of cluster membersN=30, given the minimum value of the average detection density A¯d=2.5, select the number of radars m=10 and set the maximum number of radar tracking targets Nc=8. In this case, the average detection density Ad≈2.7>A¯d. According to the above parameters, a detection confirmation matrix Ω is determined. The graphical representation of the matrix Ω is shown in Figure 11, in which 30 columns of the grid represent 30 targets and 10 rows represent 10 radars. Those marked with color in the grid are regarded as having detection relationship.Figure 11
Graphical representation of matrixΩ.Set the initial state covariance of cluster targetsP0=diag1×104,1×103,1×104,1×103T, the measurement noise covariance Rk=diag4×102,4×102T, and the potential force noise Sk=diag5,2.5,2.5,5T×10−4. The clutter obeys the Poisson distribution of the expected value λc=8 in time.Theδ-GLMB algorithm parameters are set as follows: target detection probability pD=0.98 and target survival probability ps=0.99. The newborn target model is determined according to the measurement set at time k, i.e., πγ=rγi,pγii=1Zk, where the newborn probability rγi=0.02 and state distribution of newborn targets pγix=Nx;mgi,Pg. Here, mγi=zk,i,0T and Pγ=diag1×104,5×103,1×104,5×103T.In the process of tracking 100 clusters,δ-GLMB filter runs 1000 times in total. The main filtering conditions are shown in Figure 12, and the proportion of each condition is shown in Table 2.Figure 12
δ-GLMB filtering condition.
(a)(b)(c)(d)Table 2
Proportion of filtering condition.
ConditionProportionConditionProportionFigure12(a)57.6%Figure12(b)28.3%Figure12(c)5.7%Figure12(d)7.1%Others1.3%As shown in Figure12, under the set parameters, the δ-GLMB algorithm can track all 8 assigned cluster members with a proportion of 91.6%. Figure 12(a) shows that the filter can accurately track all targets in a short time and achieve stable tracking in the whole process of target movement. This condition is the most ideal, accounting for 57.6% of all major conditions. Figure 12(b) shows that the filter can achieve stable tracking of 8 targets after 20 s. This is due to the slow acquisition of the target by the filter because of the wide distribution of the target, which is unfavorable to the accurate tracking of the target in the early stage of motion. Figure 12(c) reflects that in the later stage of tracking, the target filtering results of some filters have the problems of inaccurate track association and filtering tracks crossing, which is related to the close spacing of the assigned targets. The tracking condition shown in Figure 12(d) is not ideal. The filter has been unable to track the 2 assigned targets, accounting for 7.1%. However, we need not worry too much about this, because other filters may obtain the motion estimation of these targets and make up for these defects in the fusion stage. Other tracking conditions include discontinuous tracks, multiple track initiation, and other problems, accounting for a small proportion, which is caused by the defects of the random finite set tracking algorithm itself.Accumulate the number of targets tracked by 10 radars to obtain the total number of targets that can be tracked by the radar system, and calculate the average value of the total number of 100 experiments. The results are shown in Figure13.Figure 13
Total number of targets tracked.It can be seen from Figure13 that the radar system can only track about 20 targets at the initial time when Ad≈2.7. Then, the total number of targets tracked increases sharply to about 28. After 15 s, almost all targets can be tracked stably.


## 5.3. Estimation Fusion Results Using Graph-LSTMs

### 5.3.1. Network Training Results

In order to verify the effectiveness of Graph-LSTMs for multiradar estimation fusion, set the label space lengthdΔLx=dΔLy=2 m and dΔLvx=dΔLvy=0.5 m/s. Set the number of labels NL=11 on each dimension, that is, 0,1,2,⋯,10. The training set is used to train the network parameters, and the test set is used to test the fusion effect. The results are shown in Figures 14 and 15.Figure 14
Training effect of 1-MLSL.
(a)(b)Figure 15
Training effect of different depth networks.Figure14 shows the label prediction accuracy of 1-MLSL with different training times. 1-MLSL indicates MLSL with depth D=1. Figure 14(b) shows the mean and standard deviation of the probability of correct label output from Softmax layer. It can be seen from Figure 14(a) that the more training times, the better the network fusion effect. When the training times are less than 10000, the prediction accuracy of target labels increases rapidly with the increase of training times. When the number of training times is more than 10000, the prediction accuracy increases slowly, and the network parameter training is completed. Figure 14(b) shows that after the accuracy of label prediction increases to a certain level, the prediction probability of correct label can only increase to about 0.38. However, with the increase of training times, the standard deviation decreases, and the network prediction is more stable.The curves shown in Figure15 are the fitting result of the prediction accuracy of the target label under different depths and different training times. Figure 14 shows that when the training times are less, the network prediction effect with depth D=1 is better. When the number of training times increases, the prediction effect of the network with depth D=1 improves significantly and exceeds the effect of D=1. The effect of depth D=1 is always slightly worse than that of D=1,2. This shows that when the number of training is less, the simple network that only considers the detection and tracking relationship between radar and target can obtain better fusion effect, and with the increase of training times, the complex network that can mine the cooperative relationship between targets can play a better prediction effect. When the network depth D=3, the relationship between radars is considered, which may have a better effect on mining deeper information, but it would take a lot of training time.


### 5.3.2. Estimation Fusion Effect

In this experiment, we compare the prediction probability of the target true label by using the Graph-LSTM algorithm proposed in this paper, the commonly used EM algorithm [40, 41] and the mean method, and give the corresponding OSPA distance in target tracking. The Graph-LSTM algorithm sets the depth D=2. The OSPA distance is used to evaluate the tracking performance of the algorithm in the experiment. The OSPA distance is defined as follows [42].At timek, the true state set of cluster members is Xk=xk,1,xk,2,⋯,xk,m, and the target estimation state set given by the estimation algorithm is Xkk=xkk,1,xkk,2,⋯,xkk,n, where m and n represent the actual number of targets and the estimated number of targets in the cluster, respectively. xk,i∈Xk and xkk,i∈Xkk represent the actual value and estimated value of the target i state vector, respectively. The OSPA distance between the true state set and the estimated state set is calculated as follows:
(52)OSPAp,cXk,Xkk=1nminπ∈Πn∑i=1mdcxk,i,xkk,πip+n−m⋅cp1/p,m≤n,OSPAp,cXkk,Xk,m>n,where Πn=pnm represents all permutations of m elements in the state set.
(53)dcxk,i,xkk,πi=minc,xk,i−xkk,πi.minπ∈Πn represents the group with the smallest distance difference between the true point trace and the estimated point trace of all targets. c and p are distance sensitivity parameters and correlation sensitivity parameters, respectively.Setc=100 and p=1 in the experiment; the network fusion effect and OSPA distance are shown in Figure 16.Figure 16
Network fusion effect and OSPA distance.Figure16 shows the prediction accuracy of target state label and the average OSPA distance of 20 test sets in tracking. It can be seen from Figure 16 that at the initial time, the prediction accuracy of target state label is low and the corresponding OSPA distance is large because the radar system is difficult to track all targets. As the filtering process proceeds, all targets are tracked and the OSPA distance decreases. At about 80 s, the OSPA distance increases due to correlation errors in some filters.In the whole tracking process, the Graph-LSTM algorithm used in this paper has better prediction accuracy of correct labels than the comparison methods, which is related to the fact that the Graph-LSTM algorithm with depthD=2 can consider the cooperative interaction between cluster targets. Accordingly, the method of using Graph-LSTM algorithm to calculate the fusion has less tracking error, and it can better alleviate the problem of inaccurate target tracking caused by correlation error in the later stage of tracking.


## 5.3.1. Network Training Results

In order to verify the effectiveness of Graph-LSTMs for multiradar estimation fusion, set the label space lengthdΔLx=dΔLy=2 m and dΔLvx=dΔLvy=0.5 m/s. Set the number of labels NL=11 on each dimension, that is, 0,1,2,⋯,10. The training set is used to train the network parameters, and the test set is used to test the fusion effect. The results are shown in Figures 14 and 15.Figure 14
Training effect of 1-MLSL.
(a)(b)Figure 15
Training effect of different depth networks.Figure14 shows the label prediction accuracy of 1-MLSL with different training times. 1-MLSL indicates MLSL with depth D=1. Figure 14(b) shows the mean and standard deviation of the probability of correct label output from Softmax layer. It can be seen from Figure 14(a) that the more training times, the better the network fusion effect. When the training times are less than 10000, the prediction accuracy of target labels increases rapidly with the increase of training times. When the number of training times is more than 10000, the prediction accuracy increases slowly, and the network parameter training is completed. Figure 14(b) shows that after the accuracy of label prediction increases to a certain level, the prediction probability of correct label can only increase to about 0.38. However, with the increase of training times, the standard deviation decreases, and the network prediction is more stable.The curves shown in Figure15 are the fitting result of the prediction accuracy of the target label under different depths and different training times. Figure 14 shows that when the training times are less, the network prediction effect with depth D=1 is better. When the number of training times increases, the prediction effect of the network with depth D=1 improves significantly and exceeds the effect of D=1. The effect of depth D=1 is always slightly worse than that of D=1,2. This shows that when the number of training is less, the simple network that only considers the detection and tracking relationship between radar and target can obtain better fusion effect, and with the increase of training times, the complex network that can mine the cooperative relationship between targets can play a better prediction effect. When the network depth D=3, the relationship between radars is considered, which may have a better effect on mining deeper information, but it would take a lot of training time.


## 5.3.2. Estimation Fusion Effect

In this experiment, we compare the prediction probability of the target true label by using the Graph-LSTM algorithm proposed in this paper, the commonly used EM algorithm [40, 41] and the mean method, and give the corresponding OSPA distance in target tracking. The Graph-LSTM algorithm sets the depth D=2. The OSPA distance is used to evaluate the tracking performance of the algorithm in the experiment. The OSPA distance is defined as follows [42].At timek, the true state set of cluster members is Xk=xk,1,xk,2,⋯,xk,m, and the target estimation state set given by the estimation algorithm is Xkk=xkk,1,xkk,2,⋯,xkk,n, where m and n represent the actual number of targets and the estimated number of targets in the cluster, respectively. xk,i∈Xk and xkk,i∈Xkk represent the actual value and estimated value of the target i state vector, respectively. The OSPA distance between the true state set and the estimated state set is calculated as follows:
(52)OSPAp,cXk,Xkk=1nminπ∈Πn∑i=1mdcxk,i,xkk,πip+n−m⋅cp1/p,m≤n,OSPAp,cXkk,Xk,m>n,where Πn=pnm represents all permutations of m elements in the state set.
(53)dcxk,i,xkk,πi=minc,xk,i−xkk,πi.minπ∈Πn represents the group with the smallest distance difference between the true point trace and the estimated point trace of all targets. c and p are distance sensitivity parameters and correlation sensitivity parameters, respectively.Setc=100 and p=1 in the experiment; the network fusion effect and OSPA distance are shown in Figure 16.Figure 16
Network fusion effect and OSPA distance.Figure16 shows the prediction accuracy of target state label and the average OSPA distance of 20 test sets in tracking. It can be seen from Figure 16 that at the initial time, the prediction accuracy of target state label is low and the corresponding OSPA distance is large because the radar system is difficult to track all targets. As the filtering process proceeds, all targets are tracked and the OSPA distance decreases. At about 80 s, the OSPA distance increases due to correlation errors in some filters.In the whole tracking process, the Graph-LSTM algorithm used in this paper has better prediction accuracy of correct labels than the comparison methods, which is related to the fact that the Graph-LSTM algorithm with depthD=2 can consider the cooperative interaction between cluster targets. Accordingly, the method of using Graph-LSTM algorithm to calculate the fusion has less tracking error, and it can better alleviate the problem of inaccurate target tracking caused by correlation error in the later stage of tracking.


## 5.4. More Simulation Samples

In order to reflect the improving effect of our method on the accuracy of obtaining multiradar fusion estimates, we set more cluster motion trajectory data and obtain more fusion results of target estimators. There are four starting points of cluster motion in the sensor’s field of view, namely,±8×103,±8×103, and the parameters of the cluster motion model are consistent with those above. The starting point and direction of cluster motion are random. We obtain 20 cluster motion trajectories. Cluster center trajectories are shown in Figure 17.Figure 17
Cluster center trajectories.Six detection allocation schemes are randomly generated for each cluster, and a set of detection schemes is shown in Figure18. We can obtain 120 sets of multiradar estimates. The proposed method, mean method, and EM algorithm are, respectively, used to fuse multiple estimates of the target. The average OSPA distance of the fusions is shown in Figure 19.Figure 18
A set of cluster detection schemes.Figure 19
The average OSPA distance of the fusions.According to Figure19, the average OSPA distance medians of the proposed method, EM algorithm, and mean method are 22.4, 27.8, and 34.3, respectively. Our method is superior to the mean method and EM algorithm in terms of fusion accuracy of estimates, which reflects the effectiveness and superiority of our method for multiple estimate fusion. However, the effect of the EM algorithm is more stable. This is because there are large differences between some test data and training data in the whole estimation process, and the Graph-LSTM algorithm parameters cannot be accurately obtained by learning, which makes our method not effective. In order to further improve the robustness of the proposed algorithm, the generalization ability of the algorithm can be improved by increasing the number of training samples or pretraining.


## 6. Discussion

### 6.1. Selection ofd△L

In the estimation fusion method using Graph-LSTMs proposed in this paper,dΔL is the key parameter to determine the fusion effect. Compared with dΔLv, dΔLp play a more major role. In order to compare the effects of different dΔLx and dΔLy on the estimation fusion, set dΔLx=dΔLy=dΔLp=1,2,5,10 m for experiments under the above experimental conditions. The experimental results are shown in Figure 20.Figure 20
Tracking effect at differentd△Lp.It can be seen from Figure20 that different dΔLp have different label prediction accuracy and the OSPA distance. Generally, the larger the dΔLp, the higher the label prediction accuracy. However, when dΔLp is too large, the OSPA distance cannot be reduced accordingly. The main reason for the above phenomenon is that the larger the dΔLp, the larger the actual space occupied by each label. Radar estimate and target true state are more likely to have the same label. Therefore, the higher the label prediction accuracy can be achieved. However, the essence of state space labeling is to use the state at the label to represent all the states in the state space represented by the label. According to (23), the fusion state is related to the actual state space occupied by each label. If it is larger, the OSPA distance between the fusion state and the true state of the target may be larger, and the performance of the algorithm would be reduced. For this experiment, when dΔLp=1, the OSPA distance is large due to the low accuracy of label prediction, while when dΔLp=10, the OSPA distance is large due to the excessive setting of dΔLp. dΔLp=2 is a better choice that can balance the prediction accuracy and OSPA distance under the experimental conditions of this paper. The determination of dΔLp is mainly related to the distance between the estimations of multiradar and the true state of target. In general, we can obtain some prior information to determine dΔLp from the posterior probability of δ-GLMB filter.


### 6.2. Selection ofAd

When the number of radars is constant in multiradar system, the number of members in cluster and the number of targets tracked by single radar are both important factors affecting the fusion accuracy of cluster member estimation, and the average detection densityAd reflects the relative quantity relationship between them. In order to measure the influence of Ad, the number of targets estimated by a single radar is gradually increased under the experimental conditions of 10 radars and 30 cluster members. The OSPA distance of the fusion at different average detection densities is shown in Figure 21.Figure 21
The OSPA distance of the fusion at different average detection densities.It can be seen from Figure21 that at the initial stage, with the increase of average detection density, the fusion error of cluster target state estimation decreases gradually. OSPA=17.9 is the minimum fusion error when Nc=11. This is because the fusion becomes more accurate as the number of target estimates increases when the tracking load is lower than the maximum tracking capability of a single radar. However, with the further increase of Nc, the tracking load exceeds the maximum tracking capability of radar. The inaccurate association increases first, and the estimation error starts to increase slowly. Subsequently, the proportion of missed detection gradually increases, and the actual correct tracking density does not increase with the detection density but decreases, and the fusion error of estimates increases rapidly.


### 6.3. Limitations and Improvements of Our Method

Our proposed radar-target assignment method preliminarily solves the problem of full coverage of cluster members, but the whole assignment process is given randomly without considering the optimal allocation scheme. In the actual tracking process, different members of the cluster have different values, and the members who are leaders should be tracked accurately first. At the same time, different radars also have different detection capabilities. The factors that affect the detection ability of radar include radar resolution and the position of the target in the beam. Higher-resolution radars should be used to detect higher-value targets. Therefore, we regard the problem of radar-target allocation as a nonlinear integer programming problem [43], and the mathematical model is as follows:
(54)Maxf=∑t∈TVt1−∏s∈S1−pstωst(55)s.t.∑t=1Nωst≤Nc,∀s∈S(56)∑s=1Mωst≥1,∀t∈T(57)∑t=1N∑s=1Mωst≥N⋅Ad,∀s∈S,t∈T(58)ωst∈z+,∀s∈S,t∈T,where T is the cluster member set, S is the radar set, Vt represents the value of target t, and pst represents the detection and tracking capability of radar s to target t. The definitions of other parameters are consistent with those in the paper. Equation (54) is the objective function to maximize the target value of the detection target cluster. Equations (55)–(58) are constraint conditions, which represent the full coverage of radar to members, the limit of the number of radar channels, and the reasonable detection density, respectively.To solve equation (54), we can use traditional optimization algorithms, such as branch and bound method and cut-plane method. We also can use intelligent optimization methods, such as genetic algorithm, particle swarm optimization algorithm, and ant colony algorithm. The optimization results of multiradar detection allocation can effectively support the generation of better tracking results.Our algorithm has achieved good results in tracking a stable cluster target with a large number of members, but more advanced methods need to be proposed to deal with the possible splitting or merging of clusters. For large-scale clusters, splitting and merging are important moving characteristics that reflect the intelligence of clusters, which brings greater difficulty to the application of tracking algorithms. The state estimation results of the merging cluster usingδ-GLMB filtering are shown in Figure 22.Figure 22
Tracking results during cluster merging.Figure22 shows the tracking results of δ-GLMB filter when two clusters are merging. It can be seen that the filter can still achieve effective tracking when the cluster is merging with a small maneuver. However, with cluster merging, the distance between the two cluster members decreases, and the filter loses track of one member at 60 s. When using multiradar to track clusters, the exchange of information between radars can make up for the problem of losing track, but it would still lead to the decline of effective tracking density and tracking accuracy.The reason for the above phenomenon is that when the cluster merges, the original cooperative relationship is destroyed, and the newly established cooperative relationship leads to the failure of the target state equation describing the original motion. Similarly, this phenomenon occurs when a cluster is split. In order to solve the above problems, we should establish the cluster equation of motion in a multimodel way [44]. During tracking, the algorithm automatically detects whether cluster splitting or merging occurs and switches the motion model accordingly.


## 6.1. Selection ofd△L

In the estimation fusion method using Graph-LSTMs proposed in this paper,dΔL is the key parameter to determine the fusion effect. Compared with dΔLv, dΔLp play a more major role. In order to compare the effects of different dΔLx and dΔLy on the estimation fusion, set dΔLx=dΔLy=dΔLp=1,2,5,10 m for experiments under the above experimental conditions. The experimental results are shown in Figure 20.Figure 20
Tracking effect at differentd△Lp.It can be seen from Figure20 that different dΔLp have different label prediction accuracy and the OSPA distance. Generally, the larger the dΔLp, the higher the label prediction accuracy. However, when dΔLp is too large, the OSPA distance cannot be reduced accordingly. The main reason for the above phenomenon is that the larger the dΔLp, the larger the actual space occupied by each label. Radar estimate and target true state are more likely to have the same label. Therefore, the higher the label prediction accuracy can be achieved. However, the essence of state space labeling is to use the state at the label to represent all the states in the state space represented by the label. According to (23), the fusion state is related to the actual state space occupied by each label. If it is larger, the OSPA distance between the fusion state and the true state of the target may be larger, and the performance of the algorithm would be reduced. For this experiment, when dΔLp=1, the OSPA distance is large due to the low accuracy of label prediction, while when dΔLp=10, the OSPA distance is large due to the excessive setting of dΔLp. dΔLp=2 is a better choice that can balance the prediction accuracy and OSPA distance under the experimental conditions of this paper. The determination of dΔLp is mainly related to the distance between the estimations of multiradar and the true state of target. In general, we can obtain some prior information to determine dΔLp from the posterior probability of δ-GLMB filter.


## 6.2. Selection ofAd

When the number of radars is constant in multiradar system, the number of members in cluster and the number of targets tracked by single radar are both important factors affecting the fusion accuracy of cluster member estimation, and the average detection densityAd reflects the relative quantity relationship between them. In order to measure the influence of Ad, the number of targets estimated by a single radar is gradually increased under the experimental conditions of 10 radars and 30 cluster members. The OSPA distance of the fusion at different average detection densities is shown in Figure 21.Figure 21
The OSPA distance of the fusion at different average detection densities.It can be seen from Figure21 that at the initial stage, with the increase of average detection density, the fusion error of cluster target state estimation decreases gradually. OSPA=17.9 is the minimum fusion error when Nc=11. This is because the fusion becomes more accurate as the number of target estimates increases when the tracking load is lower than the maximum tracking capability of a single radar. However, with the further increase of Nc, the tracking load exceeds the maximum tracking capability of radar. The inaccurate association increases first, and the estimation error starts to increase slowly. Subsequently, the proportion of missed detection gradually increases, and the actual correct tracking density does not increase with the detection density but decreases, and the fusion error of estimates increases rapidly.


## 6.3. Limitations and Improvements of Our Method

Our proposed radar-target assignment method preliminarily solves the problem of full coverage of cluster members, but the whole assignment process is given randomly without considering the optimal allocation scheme. In the actual tracking process, different members of the cluster have different values, and the members who are leaders should be tracked accurately first. At the same time, different radars also have different detection capabilities. The factors that affect the detection ability of radar include radar resolution and the position of the target in the beam. Higher-resolution radars should be used to detect higher-value targets. Therefore, we regard the problem of radar-target allocation as a nonlinear integer programming problem [43], and the mathematical model is as follows:
(54)Maxf=∑t∈TVt1−∏s∈S1−pstωst(55)s.t.∑t=1Nωst≤Nc,∀s∈S(56)∑s=1Mωst≥1,∀t∈T(57)∑t=1N∑s=1Mωst≥N⋅Ad,∀s∈S,t∈T(58)ωst∈z+,∀s∈S,t∈T,where T is the cluster member set, S is the radar set, Vt represents the value of target t, and pst represents the detection and tracking capability of radar s to target t. The definitions of other parameters are consistent with those in the paper. Equation (54) is the objective function to maximize the target value of the detection target cluster. Equations (55)–(58) are constraint conditions, which represent the full coverage of radar to members, the limit of the number of radar channels, and the reasonable detection density, respectively.To solve equation (54), we can use traditional optimization algorithms, such as branch and bound method and cut-plane method. We also can use intelligent optimization methods, such as genetic algorithm, particle swarm optimization algorithm, and ant colony algorithm. The optimization results of multiradar detection allocation can effectively support the generation of better tracking results.Our algorithm has achieved good results in tracking a stable cluster target with a large number of members, but more advanced methods need to be proposed to deal with the possible splitting or merging of clusters. For large-scale clusters, splitting and merging are important moving characteristics that reflect the intelligence of clusters, which brings greater difficulty to the application of tracking algorithms. The state estimation results of the merging cluster usingδ-GLMB filtering are shown in Figure 22.Figure 22
Tracking results during cluster merging.Figure22 shows the tracking results of δ-GLMB filter when two clusters are merging. It can be seen that the filter can still achieve effective tracking when the cluster is merging with a small maneuver. However, with cluster merging, the distance between the two cluster members decreases, and the filter loses track of one member at 60 s. When using multiradar to track clusters, the exchange of information between radars can make up for the problem of losing track, but it would still lead to the decline of effective tracking density and tracking accuracy.The reason for the above phenomenon is that when the cluster merges, the original cooperative relationship is destroyed, and the newly established cooperative relationship leads to the failure of the target state equation describing the original motion. Similarly, this phenomenon occurs when a cluster is split. In order to solve the above problems, we should establish the cluster equation of motion in a multimodel way [44]. During tracking, the algorithm automatically detects whether cluster splitting or merging occurs and switches the motion model accordingly.


## 7. Conclusion

In this paper, we propose a method of using multiradar to track cluster targets. According to the characteristics of cluster target motion, we simulate the cooperative interaction behavior with multivariable stochastic differential equations and obtain the state space motion model of cluster members that can be used for iterative tracking algorithm. In the construction of the detection relationship between multiradar and cluster members, this paper proposes a radar-target assignment method. After a single radar obtains the assigned target, the target state is estimated byδ-GLMB filter. After obtaining the estimate, the label of the estimate can be obtained by using the state space labeling method proposed in this paper. On this basis, multiple estimation labels of the same target are input into the Graph-LSTMs applied to obtain the fusion label, and then, the joint estimate of target states by multiradar is solved. In the simulation experiment, we simulate the cluster motion and analyze the estimation of the filtering algorithm. Then, we compare the fusion effect of the Graph-LSTM algorithm applied in this paper with the EM algorithm and the mean method and discuss the influence of parameters dΔL and Ad on the fusion effect. The results show that multiradar can effectively track cluster targets and obtain more accurate estimation fusion in the method proposed in this paper. However, the radar-target assignment method proposed in this paper is not the optimal detection situation, and the method does not consider the impact of different value members on detection. In the discussion, we explore the improvement of the assignment method. In addition, the tracking algorithm adopted in this paper does not take into account the merging or splitting of cluster targets, which would be the focus of our next research.



---

*Source: 8556477-2022-11-14.xml*