A Grey Wolf Optimization-Based Track-Before-Detect Method for Maneuvering Extended Target Detection and Tracking

A grey wolf optimization-based track-before-detect (GWO-TBD) method is developed for extended target detection and tracking. The aim of the GWO-TBD is tracking weak and maneuvering extended targets in a cluttered environment using the measurement points of an air surveillance radar. The optimal solution is the trajectory constituted by the points of an extended target. At the beginning of the GWO-TBD, the measurements of each scan are clustered into alternative sets. Secondly, closely sets are associated for tracklets. Each tracklet equals a candidate solution. Thirdly, the tracklets are further associated iteratively to find a better solution. An improved GWO algorithm is developed in the iteration for removal of unappreciated solution and acceleration of convergence. After the iteration of several generations, the optimal solution can be achieved, i.e. trajectory of an extended target. Both the real data and synthetic data are performed with the GWO-TBD and several existing algorithms in this work. Result infers that the GWO-TBD is superior to the others in detecting and tracking maneuvering targets. Meanwhile, much less prior information is necessary in the GWO-TBD. It makes the approach is engineering friendly.


Introduction
Maneuvering weak target detection and tracking is always a challenging problem in modern radar systems. Its purpose is to detect, track and identify targets from sequences of measurements and clutter. For the increased resolution of modern radar, radars are able to receive more than one measurement per time step from different corner reflectors of a single target. Various algorithms have been developed for extended target detection and tracking. The algorithms are mainly fall into two categories: extended target probability hypothesis density (ET-PHD) filters [1][2][3][4][5] and track-before-detect algorithms [6][7][8][9][10].
ET-PHD-based algorithms [1][2][3][4][5] are capable of estimating the target extent and measurement rates as well as the kinematic state of the target. Correct partitions are significant to achieve good tracking performance. Therefore, various partitioning methods have been also developed. Reference [1] shows the application of distance partitioning algorithms for the partitions of the measurement set in PHD filters. Distance thresholds are insufficient to generate enough partitions for the correct partition. Increases of unappreciated partitions make the extended target tracking process computationally intractable. Therefore, a novel fast partitioning algorithm with fuzzy adaptive resonance theory (ART) model for the ET-PHD filter is proposed in [11]. Then, affinity propagation clustering is introduced Firstly, unlike current ET-PHD-based filters that use only the data present in the current scan, our GWO-TBD uses multiple scans (including the current scan and some past scans). A higher detection rate can thus be achieved. Secondly, the GWO-TBD has fewer parameters to adjust. It makes the approach more flexibility and engineering-friendly. Thirdly, population-based metaheuristics generally have greater exploration compared to single solution-based algorithms. Fourthly, multiple candidate solutions assist each other in GWO algorithm to avoid locally optimal solutions.
The remainder of the work is organized as follows: in Section 2, models for extended target tracking are presented. Section 3 embeds the extended tracking problem into the GWO algorithm. Also, in this section, the detailed description and implementation of the GWO-TBD are presented. To evaluate the performance of the proposed method, four real scenarios and synthetic data are tested under various conditions in Section 4. Section 5 draws simple conclusions.

Target Model
The extended target state ξ k at k-th scan is defined as the triplet ξ k = (γ, x k , X k ) in [10]. Firstly, the random variable γ > 0 is the measurement rate that describes how many measurements the target, on average, generates per time scan. Secondly, x k = (p k , v k , α k ) T ∈R 4 is the kinematic state. p k describes the target's position where p k = (x k , y k ). v k denotes the velocity and α k represents the course of the targe. Finally, X k is the extension of the target and it describes the target's size and shape. The target shape is assumed to be an ellipse because it is a good combination of an informative shape model and low computational complexity. The size of the target is denoted by the major axis l and minor axis w of the ellipse. The dynamic models and sensor measurement processes related to the state of target ξ k at kth scan are given by Equations (1) and (2), respectively: where F (•) is the state propagation function and H (•) is the measurement function. Process noise σ and measurement noise ω are zero mean, white and uncorrelated Gaussian noise sequences. In (2), {z} k T denotes the measurements of extended target at kth scan. |{z} k T | is the number of elements in {z} k T . Then, according to the Poisson distribution [21] it has: The set of measurements generated by clutter is denoted by {z} k C . The set of measurements Z k obtained at time k is the collection of measurements generated by targets and clutters. Each measurement z k usually consists of a kinematic (position) measurement component (x k i , y k i ) and a time stamp records the received time t k i : The set of all the measurements in a time series is denoted by Z K , Z K = { Z k } k = 1 K and the input of the GWO-TBD is just the points set Z k .

Problem Statement
The TBD algorithm is a method to improve the detection of weak targets by integrating their signal returns over multiple consecutive scans, i.e. estimating the state of targets at each scan ξ by measurement Z k . The optimal estimation which has a maximum likelihood is: (5) For the enormity of the solution space of ξ, the estimation of the optimal solution ξ T is divided into two stages. Firstly, the measurements originated from the extended target are abstracted. According to the random matrix approach [22], the measurements of one scan should be clustered into sets. The measurements in a set are potentially generated by the extended target. Although two improved measurement partition methods are developed in [11,12], distance partitioning [1] is utilized here for its simplification and robustness. The sets of measurements partitioned by a single distance are represented by S k 1 , S k 1 , . . . , S k Mk where M k denotes the number of sets in k-th scan. One measurement must belong to only one set, i.e. (7): where z k i,j means j-th point in i-th set at the k-th scan. Similarly, alternative distance partitioning can be obtained by multiple distances. The partition result of the i-th distance partitioning can be represented The quantity of alternative partitions is represented by M k . Then, it has: where z k i,j,n means the n-th point in the j-th set under i-th partitioning distance in the k-th scan.
M k i denotes the number of sets if the measurements are partitioned by the i-th distance in the k-th scan. Meanwhile, the quantity of partition in this scan is represented by M k , it also means the quantity of partition distance. The alternative distance partitioning can be illustrated by Figure 1. The measurements of kth scan (green points) can be clustered into three sets with a small distance and two sets with a larger distance. The diagram in Figure 1 also infers that the quantity of alternative partitions M k is alternative and determined by the spatial distribution of measurements.  (5) For the enormity of the solution space of ξ, the estimation of the optimal solution ξT is divided into two stages. Firstly, the measurements originated from the extended target are abstracted. According to the random matrix approach [22], the measurements of one scan should be clustered into sets. The measurements in a set are potentially generated by the extended target. Although two improved measurement partition methods are developed in [11,12], distance partitioning [1] is utilized here for its simplification and robustness. The sets of measurements partitioned by a single distance are represented by Sk 1 , Sk 1 , …, Sk Mk where Mk denotes the number of sets in k-th scan. One measurement must belong to only one set, i.e. (7): where zk i,j means j-th point in i-th set at the k-th scan.
. , , , where zk i,j,n means the n-th point in the j-th set under i-th partitioning distance in the k-th scan. Mk i denotes the number of sets if the measurements are partitioned by the i-th distance in the k-th scan. Meanwhile, the quantity of partition in this scan is represented by Mk, it also means the quantity of partition distance. The alternative distance partitioning can be illustrated by Figure 1. The measurements of kth scan (green points) can be clustered into three sets with a small distance and two sets with a larger distance. The diagram in Figure 1 also infers that the quantity of alternative partitions Mk is alternative and determined by the spatial distribution of measurements. Since, the measurements of an extended target can be described by a series of measurement sets, the association of some measurement sets can be regarded as a potential trajectory. In this context, the association equals a solution in SI algorithms and s th candidate solution can be represented by Cs: Equation (10) assumes that the target is detected and lost at the ks-t h and ke-th scan, respectively. Figure 2 is presented to further explain the structure of the candidate solution. Figure 2a and Figure  Since, the measurements of an extended target can be described by a series of measurement sets, the association of some measurement sets can be regarded as a potential trajectory. In this context, the association equals a solution in SI algorithms and sth candidate solution can be represented by C s : Equation (10) assumes that the target is detected and lost at the k s -th and k e -th scan, respectively. Figure 2 is presented to further explain the structure of the candidate solution. Figure 2a,b showcase that a single distance and multiple distances are applied to partition the measurements. The measurement sets are associated by the black lines for forming the candidate solution 1. The dotted lines are corresponding to the association for the candidate solution 2. The dashed box in Figure 2a denotes the measurements of one scan. The dashed box in Figure 2b means a partition of the measurements in one scan. In Figure 2, Ø means an empty set. The probability of no measurements of the target is received by the sensor Pm is then given as: This assumes that the target is detected with probability Pd with a sensor. The weak target usually has a large Pm. For example, Pm is approximately equal to 0.374 when γ and Pd equal 1 and 0.99 respectively. Naturally, empty set is reasonable to be selected as a measurement set in a solution. Then, the aim of this stage becomes finding the optimal solution CT by measurement set Z k : For instance, the optimal solution CT in Figure 1 can be the association of sets {Sk-2 2,1 , Sk-1 1,2 , Sk 2,3 }. Secondly, the optimal solution ξT can be estimated by measurements of the extended target CT and a smoothing filter. In this work, an orthogonal least squares fit is applied for the ξT.

GWO Algorithm
The GWO algorithm is an adaptive metaheuristic search algorithm inspired by the hunting and searching behavior if wolves. In GWO, a complete wolf pack consists of alpha (α), beta (β), delta (δ), and omega (ω) wolves. The best wolves should be treated as α, β, and δ that assist other wolves (ω) in In Figure 2, Ø means an empty set. The probability of no measurements of the target is received by the sensor P m is then given as: This assumes that the target is detected with probability P d with a sensor. The weak target usually has a large P m . For example, P m is approximately equal to 0.374 when γ and P d equal 1 and 0.99 respectively. Naturally, empty set is reasonable to be selected as a measurement set in a solution. Then, the aim of this stage becomes finding the optimal solution C T by measurement set Z k : For instance, the optimal solution C T in Figure 1 can be the association of sets {S k-2 2,1 , S k-1 1,2 , S k 2,3 }.
Secondly, the optimal solution ξ T can be estimated by measurements of the extended target C T and a smoothing filter. In this work, an orthogonal least squares fit is applied for the ξ T .

GWO Algorithm
The GWO algorithm is an adaptive metaheuristic search algorithm inspired by the hunting and searching behavior if wolves. In GWO, a complete wolf pack consists of alpha (α), beta (β), delta (δ), and omega (ω) wolves. The best wolves should be treated as α, β, and δ that assist other wolves (ω) in exploring more favorable regions of solution space. The alphas are the leaders of the pack, responsible for making decisions. The alphas' decisions are dictated to the pack. The betas are subordinate wolves that help the alpha in decision making or other activities. The best candidates to be the alpha is mostly likely the betas. The omega wolves are the scapegoats of the pack, they have to submit to all the other dominant wolves. The deltas have to submit to alphas and betas, but they dominate the omegas [23]. The rank of the wolves equals the fitness of the solutions. According to the differences in the rank of the wolves, in order to have better knowledge about the potential location of prey, the alpha, beta and delta wolves are assumed to be the best, the second best and the third best candidate solution, respectively.
where t is iteration, → A and → B are random vectors, C indicates the vector of a grey wolf, and C P is location of the prey. The random → A and → C vectors are calculated as [19]: where the components of A smaller fluctuation range of → A means a smaller step towards the optimal solution. More iteration is necessary but more likely to achieve the optimal solution. The first three best candidate solutions obtained can lead other hunters (including the omegas) to update their positions according to the position of the best search agents [24], so the states of the updated solutions of wolves are determined by Equation (17) [11]: where t shows recent iteration and C 1 , C 2 , C 3 denote the final state of the updated solutions, they are defined as in Equations (18)- (20), respectively: where C α , C β , C γ denote the locations of alpha, beta, and delta wolves, respectively in the swarm at a given iteration t, where → B 3 are defined as representative random vectors. The updating of the parameter controls the tradeoff between exploration and exploitation in the grey wolf optimizer (GWO). The parameter a is linearly decreased in each iteration to range from 2 to 0 according to Equation (24): where MaxIter is the total number of iteration allowed for the optimization and t is the iteration number.

Introduction of the GWO-TBD
In this section, the basic idea of the GWO-TBD method is presented. The input of the GWO-TBD are the three dimensional points during a period of time, such as the measurements of several successive frames. As presented in Equation (4), the three dimensional points include a two dimensional positional information and its measuring time. The output is a 3D-line that consists of 3D-points regarding the estimated location and time of the target. Two stages exist in the method, finding the measurements of the target and smoothing the trajectory. In stage one, a modified GWO is introduced to find the measurements and it is the most fundamental constituent of this GWO-TBD. The candidate solution here is an association which selects one measurement set in each scan. The formation of the candidate solution has been presented in Figure 2b. The pool of solutions is a K (number of scans) dimensional space. The quantity of candidate solution in theory equals: An exhaustive search for the optimal association of measurement set in so huge the space of solutions is unpractical. Therefore, a modified GWO is exploited in measurement selection for the optimal trajectory. The overview of the GWO-TBD is presented by following steps: Step 1, measurements are clustered into sets, each set regards the measurements as potentially originated from the extended target, by the alternative partitioning approach in [1].
Step 2, the initial population (multiple solutions) is generated with the spatial relationship of measurements.
Step 3, the fitness function of each candidate solution is calculated. The value of fitness function, by definition, is the probability that the candidate solution is the trajectory of an extended target. The global-best solution (α wolf), the second-best individual (β wolf) and the third-best individual (γ wolf) are estimated.
Step 4, the position of each individual is updated by GWO. Some unappreciated solutions would be removed.
Step 5, tracklet fusion is conducted. Two tracklets may be combined to form a better trajectory.
Step 6, the parameters in GWO-TBD is updated for better exploration and exploitation of candidate individuals.
Step 7, if stopping criterion is met, stop and output the best solution achieved so far, otherwise, go to Step 3.
Step 8, an orthogonal least squares fit is applied on the best solution for a more smooth and accurate trajectory. This step is the stage two in the GWO-TBD.
The individuals will evolve through the course of selection, fusion and updating iteratively. This is exactly what the GWO is introduced for. The accurate implementation of the GWO-TBD is presented in the following section.

Initial Population
It is desirable that the initial population be scattered uniformly over the feasible solution space, so that the algorithm can explore the whole solution space evenly. Meanwhile, some impossible solutions should be avoided. If two sets are generated by the same target, the spatial distance between the two sets should be smaller than the product of the time interval and the maximum velocity V max . It means that, two sets should not put in one candidate solution when Equation (26) holds: Then the strategy to generate the initial population can be concluded as: associate the sets at different scans randomly and uniformly following the criterion in Equation (26). N p candidate solutions generated and the i-th candidate solution is denoted by parameter C i here. Selecting a larger N p means more calculation but taking the advantages of good stability and strong search ability: It is worth noting that the initial solution is merely a tracklet which starts at scan k s and ends at k e . The current candidate solution equals a wolf in GWO.

The Fitness Function
The fitness function F C (C a ) is designed to estimate the probability that the candidate solution is the trajectory of a target. The definition of the F C (C a ) is: P s and P e denote the probability of a target is emerge and disappear, respectively. F S (S k i,j , C a ) means the fitness of a set S k i,j . It is designed inversely proportional to the difference between S k i,j and the estimated set of kth scan by the other sets in this candidate solution. Meanwhile, it has three parts corresponding to the triple state of an extended target (γ, x k , X k ): The three functions F γ (•), F x (•) and F X (•) represent the probability that the set S k i,j is generated by the target using the information of measurement rate, target position, target extension respectively.
The measurement rate γ of extended target is assumed invariant. The estimated measurement rateγ equals the average number of the measurements in these sets: Then the F γ (C a , S k i,j ) in Equation (30) can be estimated by: (30) can be estimated by: N(d;0, σ) denotes a Gaussian distribution defined over the variable d with mean 0 and covariance σ.
As to the target extension F X (C a , S k i,j ), extension state X k is described an inverse Wishart probability distribution function (PDF). IW d (X; v, V) in Equation (35) denotes an inverse Wishart pdf defined over the matrix X with degrees of freedom v and parameter matrix V: where Γ d (•) is the multivariate gamma function, d X means the dimension of the matrix X and Tr(•) denotes the trace of a matrix. Degrees of freedom v and parameter matrix V can be estimated by the measurement sets S i,j t t = k s , . . . , k e ; t = k by [26]. The lower and upper limit of F S (S k i,j , C a ), the fitness of set S k i,j , are 0 and 1 respectively. Then, the fitness of all N p candidate solutions is calculated by Equation (29). The α, β, δ solutions based on fitness Cα, C β , C δ can be found.

Selection, Abruption, and Fusion
Selection operates in a way such that each and every member of the population has a chance to be selected, but the better the fitness value of a candidate solution F C (C i ) the more chance of being selected it will have. Some unappreciated solutions can be removed.
Then abruption and fusion are performed to find the updated grey wolf positions. Tracklet abruption and fusion mimic the reproduction in the nature. Given parent wolves, selected from the survival, child wolves which inherit the merits of parents are generated. Here, some better solutions maybe generated by existing candidate solutions. In abruption, unappreciated association of measurement sets should be broken up. The fitness function of sets F S (S k i,j , C a ) of each candidate solution has been calculated in the last step. A candidate solution can be broken up at the measurement set which has the lowest fitness function if Equation (36) holds: Equation (36) assumes that the measurement set S i,j k b has the lowest fitness function and the S i,j k b divides the candidate solution C a into two shorter tracklets C a1 and C a2 : Then, solution fusion is performed. Two shorter tracklets C a1 and C a2 can be combined together to form a longer one if Equation (38) holds: The fusion can be accelerated and the abruption can be restrained by setting a larger value of p e and p s . Two examples on abruption and fusion are given in Figure 3. The α, β, δ solutions are free from the selection, abruption, and fusion.

Exploration
In this step, the position of the omega wolves is updated towards the alpha, beta and delta wolves. The omega solutions would be further optimized towards the best three solutions. Different with the updating equation of continuous GWO in Equation (17) where Crossover (C1, C2, C3, Cb) is a suitable crossover between solutions C1, C2, C3 and C1, C2, C3 are discrete vectors representing the effect of wolf move towards the alpha, beta and delta grey wolves in order. Cb is the candidate solution which is updating its selection of measurement sets. C1, C2, C3 are calculated using Equations (40), (42), and (43), respectively: where Cα d is the selection of the alpha wolf in the d-th scan. Parameter cstepα d is the continuous valued step size for d-th scan and can be calculated using sigmoidal function as in Equation (41): The α, β, δ solutions are free from the selection, abruption, and fusion.

Exploration
In this step, the position of the omega wolves is updated towards the alpha, beta and delta wolves. The omega solutions would be further optimized towards the best three solutions. Different with the updating equation of continuous GWO in Equation (17), the main updating equation here can be formulated in Equation (39): where Crossover(C 1 , C 2 , C 3 , C b ) is a suitable crossover between solutions C 1 , C 2 , C 3 and C 1 , C 2 , C 3 are discrete vectors representing the effect of wolf move towards the alpha, beta and delta grey wolves in order. C b is the candidate solution which is updating its selection of measurement sets. C 1 , C 2 , C 3 are calculated using Equations (40), (42), and (43), respectively: where C α d is the selection of the alpha wolf in the d-th scan. Parameter cstep α d is the continuous valued step size for d-th scan and can be calculated using sigmoidal function as in Equation (41): where A α d , D α d are calculated using Equations (15) and (21) in original GWO algorithm. Similarly, it has: A simple stochastic crossover strategy is applied each scan to crossover C 1 , C 2 , C 3 solutions as shown in Equation (44): where rand is a random number drawn from uniform distribution in the range [0, 1]. The roadmap of the GWO-TBD is presented in Figure 4 for better description. Meanwhile, the diagram of utilizing the GWO-TBD in a real scenario is presented in Figure 5. In Figure 5a,c, the measurements are represented by the blue points. The blue points are showcased in a 3-dimensional Cartesian coordinate system. Values of points on thhe x-y axes represent the location of the target, while the value of the third axis denotes the measuring time. In Figure 5a, the initial candidate solutions are represented by the green arrows, each of which is associating two measurement sets in different scans. The red point on the tail of the green arrows denotes the centroid of the measurement set which is used to form the initial tracklet.
where rand is a random number drawn from uniform distribution in the range [0, 1]. The roadmap of the GWO-TBD is presented in Figure 4 for better description. Meanwhile, the diagram of utilizing the GWO-TBD in a real scenario is presented in Figure 5. In Figure 5a,c, the measurements are represented by the blue points. The blue points are showcased in a 3-dimensional Cartesian coordinate system. Values of points on thhe x-y axes represent the location of the target, while the value of the third axis denotes the measuring time. In Figure 5a, the initial candidate solutions are represented by the green arrows, each of which is associating two measurement sets in different scans. The red point on the tail of the green arrows denotes the centroid of the measurement set which is used to form the initial tracklet.

Synthetic Result
In order to evaluate the performance of the GWO-TBD algorithm, 200 Monte Carlo numerical simulations are performed on an Intel Core I7-4790 3.6 GHz CPU, equipped with 4 GB RAM in the MatLab R2016a environment. The specific trajectory of an airplane travelling at a constant speed is presented in Figure 6a. The airplane flies on a straight line at the beginning and maneuvers during the 21-st scan to the 40-th scan and during the 51-st scan to the 70-th scan. The whole trajectory can be divided into four parts, stage 1 (straight line, 1-st-20-th scans), stage 2 (maneuvering, 21-st-40-th scans), stage 3 (straight line, 31 -st 50-th scans) and stage 4 (highly maneuvering, 51-st-70-th scans). It is worth noting that the proposed approach is designed for detecting 2-dimensional trajectories of the target, the third axis in Figure 6a is the measuring time, not the altitude of the target.
In this work, six scenarios are considered to validate both the accuracy and robustness of the algorithms, the measurement rate of targets, the measurement noise and the false alarm rate are varied in each scenario. The detailed parameters of the scenarios are presented in Table 1. It is worth noting that the probability of no measurements generated by a target in a scan equals 36.78%, 13.53% and 1.83% when the measurement rate equals 1, 2 and 4, respectively. It is hard to detect or track an extended target if no measurements are generated by it. Meanwhile, a larger measurement noise means a larger localisation error. Larger measurement noise and dense clutter would significantly deteriorate the tracking performance. The parameters of the radar in this work is patched in Table A1.

Synthetic Result
In order to evaluate the performance of the GWO-TBD algorithm, 200 Monte Carlo numerical simulations are performed on an Intel Core I7-4790 3.6 GHz CPU, equipped with 4 GB RAM in the MatLab R2016a environment. The specific trajectory of an airplane travelling at a constant speed is presented in Figure 6a. The airplane flies on a straight line at the beginning and maneuvers during the 21-st scan to the 40-th scan and during the 51-st scan to the 70-th scan. The whole trajectory can be divided into four parts, stage 1 (straight line, 1-st-20-th scans), stage 2 (maneuvering, 21-st-40-th scans), stage 3 (straight line, 31-st50-th scans) and stage 4 (highly maneuvering, 51-st-70-th scans). It is worth noting that the proposed approach is designed for detecting 2-dimensional trajectories of the target, the third axis in Figure 6a is the measuring time, not the altitude of the target.
In this work, six scenarios are considered to validate both the accuracy and robustness of the algorithms, the measurement rate of targets, the measurement noise and the false alarm rate are varied in each scenario. The detailed parameters of the scenarios are presented in Table 1. It is worth noting that the probability of no measurements generated by a target in a scan equals 36.78%, 13.53% and 1.83% when the measurement rate equals 1, 2 and 4, respectively. It is hard to detect or track an extended target if no measurements are generated by it. Meanwhile, a larger measurement noise means a larger localisation error. Larger measurement noise and dense clutter would significantly deteriorate the tracking performance. The parameters of the radar in this work is patched in Table A1. Figure 6b-d show the synthetic data in scenario 2, scenario 3 and scenario 4. It is difficult to detect the targets in these scenarios with the naked eyes. The parameters of the radar in this work is presented in Table A1 of the Appendix A.  Figure 6b, Figure 6c and Figure 6d show the synthetic data in scenario 2, scenario 3 and scenario 4. It is difficult to detect the targets in these scenarios with the naked eyes. The parameters In this work, both the ET-PHD filter-based approaches and track-before-detect methods are compared with the GWO-TBD. In the category of PHD approaches, the distance partition method [1], the ART partition method [11] and the AP partition method [12] are combined with the ET-PHD [1] respectively. In the category of TBD approaches, the 3DHT-TBD [10] and the 4DHT-TBD [9] are used. Parameters, such as the false alarm rate, measurement rate, initial state of the extended target, In this work, both the ET-PHD filter-based approaches and track-before-detect methods are compared with the GWO-TBD. In the category of PHD approaches, the distance partition method [1], the ART partition method [11] and the AP partition method [12] are combined with the ET-PHD [1] respectively. In the category of TBD approaches, the 3DHT-TBD [10] and the 4DHT-TBD [9] are used. Parameters, such as the false alarm rate, measurement rate, initial state of the extended target, are fed to the PHD filter before its iteration. The optimal sub-pattern assignment (OSPA) distance [27] is used for evaluating the performance of the algorithms. The OSPA distance between the positions of n targets T = { T 1 , T 2 , . . . , T n } and the estimated positions p = { p 1 , p 2 , . . . , p n } in each scan can be calculated by: The results of the six algorithms at each scan are presented in Figure 7. Figure 7a corresponds to scenario 1, and so on. A smaller OSPA distance means a better tracking performance. In stage 1, the performance of the GWO-TBD is similar to the others. In stage 1 of scenario 4, the GWO-TBD is worse than the others. This is mainly because the initial state of the extended target is given in the ET-PHD filters. The Hough Transformation (HT)-based methods are intentionally designed for straight line detection. Meanwhile, the measurement rate of scenario 4 equals 1. With GWO-TBD it is hard to find the optimal trajectory because merely a few points are generated by the target. In the maneuvering stages, (stage 2 and stage 4), the performance of the ET-PHD filters and HT-based methods is greatly deteriorated. Especially in the stage 4, the HT-based methods could barely detect the target. However, there was almost no effect of maneuvering on the tracking of GWO-TBD. In stages 2 and 4 of all scenarios, the GWO-TBD performance is superior to the others. In stage 3, the ET-PHD filters are inferior to the GWO-TBD because the target is lost in stage 2 and detecting the trajectory in such scenarios is difficult. The results of the six algorithms at each scan are presented in Figure 7. Figure 7a corresponds to scenario 1, and so on. A smaller OSPA distance means a better tracking performance. In stage 1, the performance of the GWO-TBD is similar to the others. In stage 1 of scenario 4, the GWO-TBD is worse than the others. This is mainly because the initial state of the extended target is given in the ET-PHD filters. The Hough Transformation (HT)-based methods are intentionally designed for straight line detection. Meanwhile, the measurement rate of scenario 4 equals 1. With GWO-TBD it is hard to find the optimal trajectory because merely a few points are generated by the target. In the maneuvering stages, (stage 2 and stage 4), the performance of the ET-PHD filters and HT-based methods is greatly deteriorated. Especially in the stage 4, the HT-based methods could barely detect the target. However, there was almost no effect of maneuvering on the tracking of GWO-TBD. In stages 2 and 4 of all scenarios, the GWO-TBD performance is superior to the others. In stage 3, the ET-PHD filters are inferior to the GWO-TBD because the target is lost in stage 2 and detecting the The HT-based methods can obtain a better performance only in scenario 4 because of the deterioration caused by the low measurement rate in the GWO-TBD. In general, the average OSPA distance of the GWO-TBD is less than those of the others in all scenarios. The detailed values of the OSPA distance are listed in Table A2. The lowest OSPA distance in each scenario is emphasized in boldface. The measurement noise of scenario 2 is smaller than that of scenario 3. The comparison between the two scenarios infers that a lower OSPA distance can be achieved under a low measurement noise because the points are more centralized. Similarly, a comparison between scenario 2 and scenario 6 showcases that the higher the false alarm rate, the lower the tracking performance.
The performance of ET-PHD filters is related to the parameters of the scenario and movement of targets. With ET-PHD filters it is hard to achieve a satisfying result when the measurement noise or the false alarm rate is high, and when the target is weak or maneuvering. The measurement noise and the false alarm rate have little influence on the two HT-based methods and the performance would be greatly deteriorated when the target is maneuvering. The GWO-TBD can cope with the difficulties and is superior to the others in almost all the stages and scenarios (except stage 1 of scenario 4).
Meanwhile, the parameters in the ET-PHD filters, such as the measurement rate and false alarm rate, have been set to fit the simulated data of each scenario. Table A3 in the Appendix The HT-based methods can obtain a better performance only in scenario 4 because of the deterioration caused by the low measurement rate in the GWO-TBD. In general, the average OSPA distance of the GWO-TBD is less than those of the others in all scenarios. The detailed values of the OSPA distance are listed in Table A2. The lowest OSPA distance in each scenario is emphasized in boldface. The measurement noise of scenario 2 is smaller than that of scenario 3. The comparison between the two scenarios infers that a lower OSPA distance can be achieved under a low measurement noise because the points are more centralized. Similarly, a comparison between scenario 2 and scenario 6 showcases that the higher the false alarm rate, the lower the tracking performance.
The performance of ET-PHD filters is related to the parameters of the scenario and movement of targets. With ET-PHD filters it is hard to achieve a satisfying result when the measurement noise or the false alarm rate is high, and when the target is weak or maneuvering. The measurement noise and the false alarm rate have little influence on the two HT-based methods and the performance would be greatly deteriorated when the target is maneuvering. The GWO-TBD can cope with the difficulties and is superior to the others in almost all the stages and scenarios (except stage 1 of scenario 4).
Meanwhile, the parameters in the ET-PHD filters, such as the measurement rate and false alarm rate, have been set to fit the simulated data of each scenario. Table A3 in the Appendix A showcases the parameter values in several scenarios. It infers that the values of parameters in the ET-PHD are various in different scenarios. However, in the GWO-TBD, similar to the 3DHT-TBD, much fewer parameters are necessary to be given an appropriate value before the iteration. Fewer parameters allow the GWO-TBD more flexibility in use. Parameter values of the 3DHT-ET-TBD and the 4DHT are also presented in Table A4 of the Appendix A. The result infers that the GWO-TBD outperforms the others, especially when the target is maneuvering and little prior information is necessary.

Results with Real Data
To evaluate the performance of the proposed algorithm further, we conduct an experiment using an air surveillance radar located in a general airport of Pucheng City, ShannXi Provience, China. Acquisition of the radar data was performed in January, 2016. The real tracks of the targets are obtained by the Global Positioning System (GPS) in the airplane. The four real trajectories obtained by GPS are presented in Figure 8a,d,g,j. showcases the parameter values in several scenarios. It infers that the values of parameters in the ET-PHD are various in different scenarios. However, in the GWO-TBD, similar to the 3DHT-TBD, much fewer parameters are necessary to be given an appropriate value before the iteration. Fewer parameters allow the GWO-TBD more flexibility in use. Parameter values of the 3DHT-ET-TBD and the 4DHT are also presented in Table A4 of the Appendix. The result infers that the GWO-TBD outperforms the others, especially when the target is maneuvering and little prior information is necessary.

Results with Real Data
To evaluate the performance of the proposed algorithm further, we conduct an experiment using an air surveillance radar located in a general airport of Pucheng City, ShannXi Provience, China. Acquisition of the radar data was performed in January, 2016. The real tracks of the targets are obtained by the Global Positioning System (GPS) in the airplane. The four real trajectories obtained by GPS are presented in Figure 8a, 8d, 8g, and 8j. The colored curves represent the movement of the target in a Cartesian coordinate system. The measuring time of the target is represented by color from red to blue. Red and blue denote the starts and ends of the trajectories respectively. The measurements of the four scenarios are presented in Figure 8b, 8e, 8h and 8k. The measurement rate of the airplane is time varying and no The colored curves represent the movement of the target in a Cartesian coordinate system. The measuring time of the target is represented by color from red to blue. Red and blue denote the starts and ends of the trajectories respectively. The measurements of the four scenarios are presented in Figure 8b,e,h,k. The measurement rate of the airplane is time varying and no measurements are generated by the airplane in some scans. Some clutter arise randomly in the surveillance area. Then, the results of the GWO-TBD are also shown in Figure 8c,f,i,l. The other methods are also applied using the real data. The initial state of the target in the ET-PHD filters is set to the correct values obtained by GPS. Actually, accurate values of the measurement rate and clutter rate are unknown. To achieve a better performance, the parameters of the ET-TBD method are different in different real scenarios.
The specific values of the parameters can be found in Table A3. The OSPA distance of the four real scenarios is presented in Figure 9. Figure 9 infers that the OSPA distance of the GWO-TBD is much lower than the others, especially when the target is maneuvering. The ET-PHD filters substantially deteriorated when no points are generated by the target or the target is maneuvering. It is worth noting that the two HT-based methods are superior to the others in scenario 3, mainly due to the fact that the target is moving in a straight line and no points are generated by the target in eight successive scans (14-th-21-st scan). In the other three real scenarios where the target is maneuvering in some scans, the GWO-TBD are significantly outperformed the other methods. Especially in the 20-th-60-th scans of scenario 4, almost only the GWO-TBD works well in tracking such a weak maneuvering extended target. Comparison between scenario 4 and scenario 1 showcases that a lower measurement rate deteriorates the performance sharply because the target is hard to be detected when few measurements is originated by it. The OSPA distance of the four scenarios is also patched in Table A5 of the Appendix A. measurements are generated by the airplane in some scans. Some clutter arise randomly in the surveillance area. Then, the results of the GWO-TBD are also shown in Figure 8c, 8f, 8i and 8l. The other methods are also applied using the real data. The initial state of the target in the ET-PHD filters is set to the correct values obtained by GPS. Actually, accurate values of the measurement rate and clutter rate are unknown. To achieve a better performance, the parameters of the ET-TBD method are different in different real scenarios. The specific values of the parameters can be found in Table A3. The OSPA distance of the four real scenarios is presented in Figure 9. Figure 9 infers that the OSPA distance of the GWO-TBD is much lower than the others, especially when the target is maneuvering. The ET-PHD filters substantially deteriorated when no points are generated by the target or the target is maneuvering. It is worth noting that the two HT-based methods are superior to the others in scenario 3, mainly due to the fact that the target is moving in a straight line and no points are generated by the target in eight successive scans (14-th-21-st scan). In the other three real scenarios where the target is maneuvering in some scans, the GWO-TBD are significantly outperformed the other methods. Especially in the 20-th-60-th scans of scenario 4, almost only the GWO-TBD works well in tracking such a weak maneuvering extended target. Comparison between scenario 4 and scenario 1 showcases that a lower measurement rate deteriorates the performance sharply because the target is hard to be detected when few measurements is originated by it. The OSPA distance of the four scenarios is also patched in Table A5 of the Appendix. Based on the experiment and analysis above, we can safely say that the GWO-TBD is more engineering friendly and better in detection and tracking performance.

Conclusions
In this article, the GWO was implemented to track and detect an extended target in a radar system. The algorithm was able to find the optimal association of measurement sets among the multiple scans. Targets can be well detected and tracked with the GWO-TBD. It is superior to the existing methods, especially when the extended target signal is weak or the target is maneuvering. Based on the experiment and analysis above, we can safely say that the GWO-TBD is more engineering friendly and better in detection and tracking performance.

Conclusions
In this article, the GWO was implemented to track and detect an extended target in a radar system. The algorithm was able to find the optimal association of measurement sets among the multiple scans. Targets can be well detected and tracked with the GWO-TBD. It is superior to the existing methods, especially when the extended target signal is weak or the target is maneuvering. Meanwhile, far less prior information is necessary before the iteration of the GWO-TBD, such as clutter rate of the surveillance area, extension and initial position of targets. Experiment infers that the GWO-TBD is better in performance and more practical in the real world. However, some limitations still exist. The GWO-TBD only copes with one target at a time. In multiple target tracking scenarios, the targets can be well detected one by one when the targets are far away from each other. The performance would be deteriorated if several maneuvering extended targets are closely distributed because several optimal solutions will exist in this scenario simultaneously. We would like to develop more approaches which can be used to detect multiple closely maneuvering extended targets from strong clutter in our later work.

Conflicts of Interest:
The authors declare no conflict of interest.