Joint Tracking and Classification of Multiple Targets with Scattering Center Model and CBMeMBer Filter

This paper deals with joint tracking and classification (JTC) of multiple targets based on scattering center model (SCM) and wideband radar observations. We first introduce an SCM-based JTC method, where the SCM is used to generate the predicted high range resolution profile (HRRP) with the information of the target aspect angle, and target classification is implemented through the data correlation of observed HRRP with predicted HRRPs. To solve the problem of multi-target JTC in the presence of clutter and detection uncertainty, we then integrate the SCM-based JTC method into the CBMeMBer filter framework, and derive a novel SCM-JTC-CBMeMBer filter with Bayesian theory. To further tackle the complex integrals’ calculation involved in targets state and class estimation, we finally provide the sequential Monte Carlo (SMC) implementation of the proposed SCM-JTC-CBMeMBer filter. The effectiveness of the presented multi-target JTC method is validated by simulation results under the application scenario of maritime ship surveillance.


Introduction
Traditionally, target tracking and target classification are treated as two independent problems, and they are usually solved separately. However, these two problems are closely related. For example, tracking affects classification by providing flight envelope information for different air target classes, while classification affects tracking via selecting appropriate class-dependent kinematic models. Therefore, a good classification may benefit tracking and vice versa. For this reason, the joint tracking and classification (JTC) method is receiving more and more attention.
By now, many JTC methods have been proposed [1][2][3][4][5][6][7][8][9][10], and these methods can roughly be divided into three categories. The first category is the most popular one and is dedicated to point targets. In this case, the resolution of the tracking sensor is very limited, and realization of target classification has to exploit attribute/identity sensor (such as electronic support measure) information or target dynamics (such as class-dependent maneuverability) [1][2][3][4][5].
The second category treats the target as an extended target, and the measurement of the target is modeled as the extent in down-range direction (length of the target) or both in down-range and lateral-range directions (i.e., size of the observed target contour), based on the assumption that the target has an ellipsoidal shape [6][7][8]. Targets are classified with the feature information of different length or size. y k ] T . The subscript k denotes sampling time and T is the sign for vector transpose. F is the state transition matrix. w k ∼ N (w; 0, Q) represents the multi-dimensional Gaussian process noise vector, where N (ζ; µ, Σ) denotes Gaussian function with variable ζ, mean µ and covariance matrix Σ. F and Q can be further written as where t is the sampling interval and q is the acceleration variance.

Sensor Observation Model
The radar can provide range and bearing measurements of the target's centroid, and the observation model of kinematic (position) measurement can be described by where h(·) is the kinematic observation function and r k and β k represent noisy measurements of target range and bearing at time k, respectively. v k = v 1,k v 2,k T denotes the corresponding zero-mean observation noise with covariance matrix R = E[v k v T k ] = diag[σ 2 r , σ 2 β ]. As shown in Figure 1, for the low-speed maritime target, the heading is almost aligned with the axial direction of the target body because of its limited maneuverability. Under this condition, the aspect angle φ k of the target can be obtained as where θ k = tan −1 ( x k ) is the heading angle.

Target Motion Model
Without loss of generality, the considered target moves on the 2D plane with a nearly constant velocity, and the evolvement process of the target state is given as x y x y =   x represents target state including position component where t is the sampling interval and q is the acceleration variance.

Sensor Observation Model
The radar can provide range and bearing measurements of the target's centroid, and the observation model of kinematic (position) measurement can be described by where ( ) h ⋅ is the kinematic observation function and k r and k β represent noisy measurements of target range and bearing at time k , respectively.
denotes the corresponding zero-mean observation noise with covariance matrix As shown in Figure 1, for the low-speed maritime target, the heading is almost aligned with the axial direction of the target body because of its limited maneuverability. Under this condition, the aspect angle k φ of the target can be obtained as  An equivalent expression of aspect angle is in the form of An equivalent expression of aspect angle is in the form of The 3D-SCM is an equivalent of the target in geometry space to the radar response in the electromagnetic field. It provides a concise and physically relevant description of the target's scattering through a set of representative scattering parts, and thus a more effective way to characterize the target's electromagnetic scattering behavior. The 3D-SCM consists of a set of scattering center with a specific position, amplitude and type parameters, and it can be represented by S = a n , α n , x n , y n , z n N n=1 where a n is the amplitude of nth scattering center, α n is a frequency-dependent factor, (x n , y n , z n ) is the corresponding 3D spatial position in the target body coordinates and N is the number of scatters involved in the model. For a specific target class c, the associated 3D-SCM can be denoted as S c .
The whole target's backscattering with respect to radar instantaneous frequency f , viewing angles (i.e., azimuth angle φ and elevation angle γ) can be expressed as [13] (j f / f c ) α n a n (φ, γ) · exp(−j4π(x n cos γ cos φ + y n cos γ sin φ + z n sin γ)/λ) (7) where λ is the wavelength, f c is the central frequency of signal, j = √ −1 is the imaginary unit and a n (φ, γ) represents the amplitude of nth scattering center which may change with φ and γ.
At a specific viewing angle of (φ, γ), the corresponding projection position of the nth scattering center at the down-range direction is r n (φ, γ) = x n cos γ cos φ + y n cos γ sin φ + z n sin γ (8) Assuming that the bandwidth of the radar signal is B and the ith discrete frequency point is where ∆F is the frequency interval and I + 1 is the total number of frequency points. Then, the frequency response of the ith frequency point can be written as E i = E( f i , φ, γ, S). After a direct operation of inverse discrete Fourier transform (IDFT) on frequency response sequence E = [E 0 , E 1 , · · · , E I ], the desired HRRP can be immediately obtained.
When the target's motion is restricted on the 2D plane, the observation model of HRRP is represented as where d denotes the (I+1)-dimensional measurement of the HRRP and each component of d corresponds to one range resolution cell, g(φ, S) IDFT(E) denotes the compact form of observation function, f i and S are known parameters, γ ≈ 0, φ can be obtained through Equations (4) or (5) and n is the observation noise vector.

CBMeMBer Filter
The CBMeMBer filter is outlined as follows and the details can be found in [18,21]. It approximates the posterior multi-target density by a multi-Bernoulli RFS. The multi-Bernoulli RFS consists of M independent Bernoulli RFSs X (i) , that is, where r (i) ∈ [0, 1] is the target existence probability and p (i) (·) is a spatial distribution. Therefore, the probability density of multi-Bernoulli RFS X is given by where n is the number of targets. The multi-Bernoulli RFS X is completely described by the multi-Bernoulli parameter set , and the probability density of the multi-Bernoulli RFSs X can be abbreviated . The CBMeMBer filter consists of a prediction step and an update step.

Prediction Step
If the posterior probability density at time k − 1 is , then the predicted multi-target density is also a multi-Bernoulli formed by the union of the multi-Bernoulli for the surviving targets and target births is the predicted multi-Bernoulli for the target births and it is usually assumed to be known. (r (i) is the predicted multi-Bernoulli for the surviving targets and it is given by where f , g = f (x)g(x)dx denotes the inner product operation, p S,k (x k−1 ) is the survival probability of the surviving targets, f k|k−1 (x k |x k−1 ) is the single-target transition density. There are M k|k−1 = M k−1 + M Γ,k−1 predicted hypothesized tracks.

Update Step
Assuming that n k,z measurements are collected as Z k = z k,1 , · · · , z k,n k,z and the predicted probability density is π k|k−1 (X) = {(r , then the posterior multi-target density at time k can be approximated by a multi-Bernoulli as Sensors 2020, 20, 1679 6 of 25 The first term in Equation (16) corresponds to the multi-Bernoulli density for the legacy tracks and it can be given by where p D,k (x k ) is the detection probability.
The second term in Equation (16) corresponds to the multi-Bernoulli density for measurement-corrected tracks and it can be given as where g k (z|x k ) is the likelihood function, κ(z) is the clutter intensity function. There are M k = M k|k−1 + n k,z updated hypothesized tracks.

JTC Method Based on SCM and CBMeMBer Filter
In this section, the SCM-based JTC method is first presented by using the HRRP as the feature for target classification. Then, the SCM-JTC-CBMeMBer filter is derived for multi-target JTC.

SCM-Based JTC Method: Single-Target Case
The joint target state can be modeled as ξ k−1 (x k−1 , c), where x k−1 is the kinematic state and c is the class label that can be taken from the set of the target classes C = {c 1 , c 2 , · · · , c n c }. n c and c m represent the total number of the target class and the mth target class, respectively. In the SCM-based JTC method, the available measurement at time k consists of kinematic (position) measurement z p k = [r, θ] T and signature (HRRP) measurement z c k = d, and the joint measurement is denoted as . The purpose of Bayesian JTC is to estimate the target state and class simultaneously at time k, under the condition that the distribution p(x k−1 , c| Z k−1 ) at time k − 1 and the measurement z k at time k are available. That is, to obtain the posterior probability-mass distribution For target tracking, the class-dependent probability density function (PDF) for a specific target class c m can be represented as Sensors 2020, 20, 1679 Accordingly, for target classification, the probability function can be obtained by To obtain the recursive equations of the SCM-based JTC method, two assumptions should be followed. Assumption 1. All the targets have the same motion model, i.e., the single state transition function is the kinematic state transition function and is decided by the system model and f c k|k−1 (c j |c m ) is the class state transition function and can be represented by the Dirac function δ(·) as Assumption 2. The kinematic measurement and HRRP measurement are independent of each other, and the kinematic measurement is independent of the target class, so the measurement likelihood can be written as is the likelihood function of HRRP measurement, and is defined as normalized correlation coefficient of observed HRRP with model-predicted HRRP. S c m is the SCM corresponding to target class c m .
Therefore, the SCM-based JTC method can be constructed through the following two steps. The prediction steps of the target state and class are given by Similarly, the update steps of target state and class are with Sensors 2020, 20, 1679 8 of 25 Because of the complexity and high nonlinearity of the observation model, there is no analytic form to obtain the recursive estimation of target state and class, so we have to resort to the SMC technique (also known as particle filter, PF). Given the particle set w j k−1 , x j k−1 , l j n k−1,p j=1 , where the superscript j denotes the index of the particles, l j ∈ C represents the class label corresponding to the jth particle and n k−1,p is the number of particles at time k − 1. The posterior target state and associated class probability can be represented by Then, a complete recursive procedure from time k − 1 to k can be summarized as Algorithm 1. Step 1. Model prediction 1) Target state prediction: where S l j represents the 3D-SCM corresponding to target class l j . Step 2. Likelihood evaluation 1) Kinematic observation likelihood: Step 3. Particle weight evaluation 1) Joint weight calculation: The posterior target state estimationx k and class probabilities p(c m | Z k ) at time k can be obtained bŷ Sensors 2020, 20, 1679 To reduce the effect of particle degeneracy, the resampling operation [33] should be considered in method implementation. Specifically, the class-dependent resampling strategy is adopted to avoid particle degeneracy caused by an incorrect target classification. In this strategy, the maximum number of particles for the SCM-based JTC method is set as L max , while the minimum number of particles for each target class is set as L min . In this paper, the standard resampling operation is used to resample each class, that is, for the particle with l j = c m ,

SCM-JTC-CBMeMBer Filter: Multi-Target Case
For the JTC of multi-target, the available measurement set at time k is denoted as Z k = z k,l n k,z l=0 and the measurement set up to time k is Z k = {Z l } k l=1 . The posterior multi-target density at time k − 1 is modeled as a multi-Bernoulli In addition to Assumption 1 and Assumption 2, the following assumptions should also be followed to obtain the SCM-JTC-CBMeMBer filter.

Assumption 3. Each target evolves motion and generates measurements independently.
Assumption 4. The clutter is modeled as Poisson RFS with Poisson average rate λ c , and it is independent of target-originated measurements. The spatial distribution of the clutter is a uniform distribution, denoted by C( z). The clutter intensity function is κ( z) = λ c C( z).

Assumption 5.
The survival and detection probabilities are state-independent, i.e., p S,k (x, c) = p S,k , p D,k (x, c) = p D,k . Assumption 6. The PDF of birth targets at time k − 1 is also a multi-Bernoulli, namely (43) Proposition 1. If the posterior multi-target density at time k − 1 is a multi-Bernoulli, as shown in Equation (42), then the predicted multi-target density is also a multi-Bernoulli and is given by Sensors 2020, 20, 1679 10 of 25 with r (i) The proof of Proposition 1 is given in Appendix A.

Proposition 2.
If the predicted multi-target density at time k is a multi-Bernoulli Then, the posterior multi-target density can be approximated by a multi-Bernoulli as The proof of Proposition 2 is shown in Appendix B. The state extraction step is similar to the CBMeMBer filter, and the details can be found in [22].

SMC Implementation of the SCM-JTC-CBMeMBer Filter
In what follows, SMC implementation of the SCM-JTC-CBMeMBer filter recursion will be presented.
Supposing that the posterior multi-target density π k−1 = r is given, and each component p Then the class probability p

Proposition 3.
Given the importance density q , then the predicted multi-target density is also multi-Bernoulli and the SMC implementation is calculated as The proof of Proposition 3 is given in Appendix C.

Proposition 4. If the predicted multi-target density is the multi-
, then the updated multi-target density is also a multi-Bernoulli, and the SMC implementation can be computed by The proof of Proposition 4 is detailed in Appendix D.
Particle resampling is also needed for SMC implementation of the SCM-JTC-CBMeMBer filter, and the same resampling strategy as introduced in Section 3.1 is used. Similar to the CBMeMBer filter, in the proposed SCM-JTC-CBMeMBer filter, the pruning and merge strategy (refer to [22] for the details) should be adopted to reduce the computational burden.

Simulation Results
The effectiveness of the proposed SCM-based JTC method and SMC-JTC-CBMeMBer filter is evaluated by simulations. The observation precisions of range and bearing are σ r = 1 m and σ β = 0.3 • , respectively. The radar is located at the origin of the coordinates with center frequency f c = 35 GHz, bandwidth B = 150 MHz, frequency interval ∆F = 1 MHz. In this paper, three (n c = 3) different ship target classes are considered. The maximum and minimum number of particles used in all the simulations are L max = 2000 and L min = 300, respectively. The CAD models and the corresponding 3D-SCMs (denoted by a red asterisk "*") are shown in Figure 2.
The proof of Proposition 4 is detailed in Appendix D.
Particle resampling is also needed for SMC implementation of the SCM-JTC-CBMeMBer filter, and the same resampling strategy as introduced in Section 3.1 is used. Similar to the CBMeMBer filter, in the proposed SCM-JTC-CBMeMBer filter, the pruning and merge strategy (refer to [22] for the details) should be adopted to reduce the computational burden.

Simulation Results
The

SCM-Based JTC Method
In this simulation, the SCM-based JTC method will be evaluated under the scenario of single-target measurement without clutter. The process noise of the target motion state is characterized by

SCM-Based JTC Method
In this simulation, the SCM-based JTC method will be evaluated under the scenario of single-target measurement without clutter. The process noise of the target motion state is characterized by q = 0.5 m/s 2 . The three ship targets have the same motion model, and the initial state is x 0 = [1.2 km 1.5 km − 7.5 m/s 5.0 m/s] T . The sampling interval is t = 1s and the total duration is 40 s. At a certain time, only one ship target is present, and the JTC performance of different target is analyzed separately.
The trajectory tracking (for a single run) and target classification (averaged by 20 Monte Carlo runs) results are shown in Figures 3-5. As can be seen from the figures, the proposed method can not only accurately estimate the state of the target but also correctly classify the target simultaneously. It is also seen that the classification probability curve, which matches with the true target class, increases rapidly, and can approach one (100%) within 30 s under all the tested conditions. Conversely, the classification probability curves mismatching the true target class decrease gradually and reach 0 after about 10-30 estimate cycles, meaning that a high confidence classification result is obtained.
However, it does not mean that we have to wait for a period of 10-30 estimate cycles to obtain a reliable decision on target class, since the classification probability matching the true target class is obviously higher than others within only a few cycles.
The trajectory tracking (for a single run) and target classification (averaged by 20 Monte Carlo runs) results are shown in Figure 3-5. As can be seen from the figures, the proposed method can not only accurately estimate the state of the target but also correctly classify the target simultaneously. It is also seen that the classification probability curve, which matches with the true target class, increases rapidly, and can approach one (100%) within 30 s under all the tested conditions. Conversely, the classification probability curves mismatching the true target class decrease gradually and reach 0 after about 10-30 estimate cycles, meaning that a high confidence classification result is obtained. However, it does not mean that we have to wait for a period of 10-30 estimate cycles to obtain a reliable decision on target class, since the classification probability matching the true target class is obviously higher than others within only a few cycles.   The trajectory tracking (for a single run) and target classification (averaged by 20 Monte Carlo runs) results are shown in Figure 3-5. As can be seen from the figures, the proposed method can not only accurately estimate the state of the target but also correctly classify the target simultaneously. It is also seen that the classification probability curve, which matches with the true target class, increases rapidly, and can approach one (100%) within 30 s under all the tested conditions. Conversely, the classification probability curves mismatching the true target class decrease gradually and reach 0 after about 10-30 estimate cycles, meaning that a high confidence classification result is obtained. However, it does not mean that we have to wait for a period of 10-30 estimate cycles to obtain a reliable decision on target class, since the classification probability matching the true target class is obviously higher than others within only a few cycles.   The trajectory tracking (for a single run) and target classification (averaged by 20 Monte Carlo runs) results are shown in Figure 3-5. As can be seen from the figures, the proposed method can not only accurately estimate the state of the target but also correctly classify the target simultaneously. It is also seen that the classification probability curve, which matches with the true target class, increases rapidly, and can approach one (100%) within 30 s under all the tested conditions. Conversely, the classification probability curves mismatching the true target class decrease gradually and reach 0 after about 10-30 estimate cycles, meaning that a high confidence classification result is obtained. However, it does not mean that we have to wait for a period of 10-30 estimate cycles to obtain a reliable decision on target class, since the classification probability matching the true target class is obviously higher than others within only a few cycles.   As a comparison, the simulation also considers the target classification result directly obtained from the HRRP correlation method, where the HRRP templates (training data) are generated from the CAD model with electromagnetic simulation tool, and the test data are predicted from the 3D-SCM. In the simulation, for each target, 360 HRRP training samples (which cover 0-360 • in azimuth angle space) are generated with 1 • interval. Accordingly, 1440 test samples are generated from each 3D-SCM at the same viewing angle with azimuth angle interval 0.2 • . The classification results are shown in (93) where N total (m) denotes the total test samples for the mth target class, N correct (m) represents the correctly classified test samples for the mth target class. Compared the classification results in Table 1 (no tracking process involved) with those shown in Figures 3-5 (with JTC processing), it is seen that the SCM-based JTC method achieves a performance improvement of more than 0.1 (10%) in PCC after the tracking filter is stable, indicating the advantage of the proposed method in classification accuracy.

SCM-JTC-CBMeMBer Filter
In this simulation, all three classes of ship targets will appear in the surveillance area.
Target A appears at time k = 5 and disappears at time k = 55 with initial state x  The trajectory tracking results for a single run are shown in Figure 6. As can be clearly seen from Figure 6, under the clutter environments, the proposed SCM-JTC-CBMeMBer filter can estimate target number and state correctly, and it can also obtain correct target classification, which will be further validated by the repeated Monte Carlo trials. The trajectory tracking results for a single run are shown in Figure 6. As can be clearly seen from Figure 6, under the clutter environments, the proposed SCM-JTC-CBMeMBer filter can estimate target number and state correctly, and it can also obtain correct target classification, which will be further validated by the repeated Monte Carlo trials. To further test the performance of the proposed SCM-JTC-CBMeBer filter, 50 Monte Carlo runs are carried out under the same scenario as above. Specifically, the metric of Optimal Subpattern Assignment (OSPA) distance [34] is used to evaluate the multi-target tracking results, and the CBMeMBer filter is also considered as a comparison.
The OSPA distance and cardinality estimation are shown in Figure 7, and the target classification results are plotted in Figure 8.
From Figure 7, we can see that the SCM-JTC-CBMeMBer filter can effectively estimate the target state and target number. At the instant when the target appears and disappears, a slight degradation in estimation performance is observed, which is the normal phenomenon confronted in multi-target tracking. Compared with the conventional CBMeMBer filter (which can only be used for multi-target tracking purposes rather than targets classification), the SCM-JTC-CBMeMBer filter has almost the same performance in target tracking. As can be seen from Figure 8, the SCM-JTC-CBMeMBer filter can also correctly classify multiple targets, and the classification probability of each target is very high (almost reaches one). To further test the performance of the proposed SCM-JTC-CBMeBer filter, 50 Monte Carlo runs are carried out under the same scenario as above. Specifically, the metric of Optimal Subpattern Assignment (OSPA) distance [34] is used to evaluate the multi-target tracking results, and the CBMeMBer filter is also considered as a comparison.
The OSPA distance and cardinality estimation are shown in Figure 7, and the target classification results are plotted in Figure 8.
From Figure 7, we can see that the SCM-JTC-CBMeMBer filter can effectively estimate the target state and target number. At the instant when the target appears and disappears, a slight degradation in estimation performance is observed, which is the normal phenomenon confronted in multi-target tracking. Compared with the conventional CBMeMBer filter (which can only be used for multi-target tracking purposes rather than targets classification), the SCM-JTC-CBMeMBer filter has almost the same performance in target tracking.
As can be seen from Figure 8, the SCM-JTC-CBMeMBer filter can also correctly classify multiple targets, and the classification probability of each target is very high (almost reaches one). The trajectory tracking results for a single run are shown in Figure 6. As can be clearly seen from Figure 6, under the clutter environments, the proposed SCM-JTC-CBMeMBer filter can estimate target number and state correctly, and it can also obtain correct target classification, which will be further validated by the repeated Monte Carlo trials. To further test the performance of the proposed SCM-JTC-CBMeBer filter, 50 Monte Carlo runs are carried out under the same scenario as above. Specifically, the metric of Optimal Subpattern Assignment (OSPA) distance [34] is used to evaluate the multi-target tracking results, and the CBMeMBer filter is also considered as a comparison.
The OSPA distance and cardinality estimation are shown in Figure 7, and the target classification results are plotted in Figure 8.
From Figure 7, we can see that the SCM-JTC-CBMeMBer filter can effectively estimate the target state and target number. At the instant when the target appears and disappears, a slight degradation in estimation performance is observed, which is the normal phenomenon confronted in multi-target tracking. Compared with the conventional CBMeMBer filter (which can only be used for multi-target tracking purposes rather than targets classification), the SCM-JTC-CBMeMBer filter has almost the same performance in target tracking. As can be seen from Figure 8, the SCM-JTC-CBMeMBer filter can also correctly classify multiple targets, and the classification probability of each target is very high (almost reaches one).

Conclusions
In this paper, an SCM-based JTC method is first introduced. The presented method can implement target tracking and classification simultaneously by using a model based on HRRP prediction, target kinematic and HRRP measurements, and thus to alleviate the dependence of target classification on the requirements of target maneuvers or other support information (such as target attribute/identity) in the conventional methods. The SCM-based JTC is then integrated into the framework of the CBMeMBer filter, and the resulting SCM-JTC-CBMeMBer filter for multi-target JTC is derived under the condition with detection uncertainty. Finally, the SMC technique is adopted to implement the proposed filter in view of the complex calculation in multi-target state recursion. Simulations are carried out under the typical scenario with three different ship targets, and the results show that the developed method can not only effectively estimate the target state, but also obtain reliable target classification decision. Additionally, the proposed joint processing method can achieve better performance than separate HRRP classification without involving target tracking.

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

Conclusions
In this paper, an SCM-based JTC method is first introduced. The presented method can implement target tracking and classification simultaneously by using a model based on HRRP prediction, target kinematic and HRRP measurements, and thus to alleviate the dependence of target classification on the requirements of target maneuvers or other support information (such as target attribute/identity) in the conventional methods. The SCM-based JTC is then integrated into the framework of the CBMeMBer filter, and the resulting SCM-JTC-CBMeMBer filter for multi-target JTC is derived under the condition with detection uncertainty. Finally, the SMC technique is adopted to implement the proposed filter in view of the complex calculation in multi-target state recursion. Simulations are carried out under the typical scenario with three different ship targets, and the results show that the developed method can not only effectively estimate the target state, but also obtain reliable target classification decision. Additionally, the proposed joint processing method can achieve better performance than separate HRRP classification without involving target tracking.

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