**4. Numerical Study**

The advantages of multi-static Doppler radar such as lightweight, wide range of surveillance with high accuracy, and low power consumption lead to its broad applications in both civilian and military applications [43–45]. However,the number of the sensors in conjunction with the non-linear nature and and low observability of the Doppler type measurement leads to many numerical difficulties [44,45]. The use of multistatic Doppler-only measurements in a scenario of 10 receivers and one cooperative transmitter has been proposed in [46] and its extended version [47] for joint detection and tracking of one target.

This numerical study based on the model mentioned in [40] with 10 marine ships. Each ship at time *<sup>k</sup>* is represented by a 5 <sup>−</sup> *<sup>D</sup>* state vector *xk* in the surveillance of interest *xk* <sup>=</sup> *pT <sup>k</sup>* , *<sup>ν</sup><sup>T</sup> <sup>k</sup>* , *α<sup>k</sup>* .*T* , where *pk* = [*μk*, *λk*] *<sup>T</sup>* and *<sup>ν</sup><sup>k</sup>* = *μ*˙ *<sup>k</sup>*, *λ*˙ *<sup>k</sup>* .*<sup>T</sup>* denote the position and velocity in the longitude and latitude, respectively; *α<sup>k</sup>* is the course of the target, and *T* denotes the transpose operation. The target dynamic model can be given as follows:

$$\mathbf{x}\_{k} = \mathbf{F}\_{k|k-1} \left( \mathbf{x}\_{k-1} \right) + \mathbf{G} \mathbf{u}\_{k} \tag{51}$$

where

$$F(\mathbf{x}\_{k-1}) = \begin{bmatrix} 1 & \frac{\sin(at)}{a} & 0 & \frac{(\cos(at)-1)}{a} & 0\\ 0 & \cos(at) & 0 & -\sin(at) & 0\\ 0 & -\frac{(\cos(at)-1)}{a} & 1 & -\frac{\sin(at)}{a} & 0\\ 0 & \sin(at) & 0 & \cos(at) & 0\\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} \mathbf{x}\_{k-1}; \qquad \mathbf{G} = \begin{bmatrix} \frac{t^2}{2} & 0 & 0\\ t & 0 & 0\\ 0 & \frac{t^2}{2} & 0\\ 0 & t & 0\\ 0 & 0 & t \end{bmatrix},\tag{52}$$

and *t* is sample period, *nk* is a Gaussian noise vector of velocity and course noise components with zero-mean. Note that latitudinal and longitudinal measurements are in degrees (◦), the distance, speed and time are given in nautical miles (*M*), knots (*kn*), and hours (*h*), respectively.

**Remark 1.** *Equation* (52) *is resulted from the assumption that the surveillance region is not very far from the Equator.*

The new births are assumed to be distributed with labeled multi-Bernoulli RFS distributions of parameters *fB* (*x*) = *r* (*i*) *<sup>B</sup>* , *p* (*i*) *B* 4 *i*=1 where *r* (*i*) *<sup>B</sup>* is the i*th* common existence probability , and *<sup>p</sup>* (*i*) *<sup>B</sup>* (*x*) = N *x*, *x*ˆ (*i*) *<sup>B</sup>* , *PB* with

$$\begin{aligned} \hat{\mathfrak{x}}\_{B}^{(1)} &= \left[15.6^{\circ}N, 0, 113^{\circ}E, 0, 0\right]^{T}; \\ \hat{\mathfrak{x}}\_{B}^{(2)} &= \left[13.2^{\circ}N, 0, 107.5^{\circ}E, 0, 0\right]^{T} \\ \hat{\mathfrak{x}}\_{B}^{(3)} &= \left[18.2^{\circ}N, 0, 110.7^{\circ}E, 0, 0\right]^{T}; \\ \hat{\mathfrak{x}}\_{B}^{(4)} &= \left[22.3^{\circ}N, 0, 118.8^{\circ}E, 0, 0\right]^{T}; \\ P\_{B} &= \text{diag}\left(\left[2^{\prime}N, 30\left(kn\right), 2^{\prime}E, 30\left(kn\right), 6\pi/180\left(rads^{-1}\right)\right]\right) \end{aligned}$$

Table 1 lists out the initial state of ten targets with random time of appearance and disappearance, and the average course is *α*¯ = 2*π*/180(*rad*/*s*).

The parameters of the dynamic model are given in Table 2.


**Table 1.** Target initial states.

**Table 2.** Parameters of the Dynamic model.


Consider the configuration of multiple Doppler sensors system including two spatially distributed receivers and one cooperative transmitter located as in Figure 2. Based on Doppler effect, this system can measure the speed of a target at a distance by calculating the altered frequency of the returned signals which originate from the emitting pulses of radio signals and being reflected to radar after reaching target [48].

**Figure 2.** Configuration of multitarget tracking using MRS.

The observation of a target state *xk* at the *sth* receiver using Doppler measurement is given by

$$z\_k^{(s)} = -\nu\_k^T \left( \frac{p\_k - p\_r^{(s)}}{\left\| \left| p\_k - p\_r^{(s)} \right| \right\|} + \frac{p\_k - p\_t}{\left\| |p\_k - p\_t| \right\|} \right) \frac{f\_t}{c} + w\_{k'} \tag{53}$$

in which *pk* and *ν<sup>k</sup>* have been defined above the Equation (51), *pt* = [*μt*, *λt*] *<sup>T</sup>* is the position of the transmitter, and *p* (*s*) *<sup>r</sup>* = [*μ*(*s*) *<sup>r</sup>* , *<sup>λ</sup>*(*s*) *<sup>r</sup>* ] *<sup>T</sup>* is the location of the *<sup>s</sup>th*−receiver ; *wk* is zero-mean Gaussian noise, *wk* ∼ N (0, **Q***k*), with covariance **Q***<sup>k</sup>* = *diag* ' [1*Hz*2] ( ; and *ft* is the signal frequency emitted from the transmitter, and *c* is the light speed.

Since the targets are dynamic in different directions, the value of observation *z* (*s*) *<sup>k</sup>* in (53) can be negative or positive in the known interval [−*f*0, +*f*0] of the Doppler sensor. In this work, the measurement space for two receivers have the same measurement space of [−200*Hz*, 200*Hz*] . The parameters of the observation model are given in Table 3. It can be seen that, not only the state equation but also the measurement one are highly nonlinear.


**Table 3.** Parameters of observation model.

By using the proposed *B*−GLMB filter, the configuration of multiple marine ships tracking using multiple Doppler radars with ground truths and their tracking results are illustrated in Figure 2. For better visualization of multiple targets, each target is assigned to a distinct color. The results of longitudinal-latitudinal co-ordinate target trajectories are demonstrated in Figure 3.

**Figure 3.** Tracking in longitudinal-latitudinal coordinates.

For evaluating the effectiveness of the proposed method comparing to the fixed-GLMB filter and the Joint-CPHD, 100 Monte - Carlo run has been used, and the distance, location and cardinality errors are calculated via Optimal Sub-Pattern Assignment, OSPA, [49] and shown in Figure 4a. By using this metric, the distance between the set of true multitarget states and that of estimated target states is calculated at each time step. For measuring the error between two set of tracks, the use of OSPA is insufficient, and OSPA(2) is needed. The OSPA(2)errors [50] of the B-GLMB and fixed-GLMB filters are compared and plotted against time in Figure 4b, respectively.

**Figure 4.** Evaluation of tracking errors using (**a**) OSPA, and (**b**) OSPA(2).

For the *B*−GLMB filter and joint-CPHD filter, the clutter rate fluctuates in the range of *λ<sup>c</sup>* = [28, 70], and the detection probability changes from 0.75 to 0.98, i.e., *pD* = [0.75, 0.98]. The fixed-GLMB filter is used with fixed *pD* of 0.75 and 0.98 and fixed *λ* of 28 and 70, respectively. The window length used in OSPA(2) to obtain the differences between the true and estimate sets of trajectories in Figure 4b is set at *wl* = 10. Both the OSPA and OSPA(2) are used with cut-off parameter *c*<sup>0</sup> = 0 and *p* = 1.

Obviously from Figure 4a, the errors in distance and location between the set of true targets and the estimated ones using *B*−GLMB filter is the smallest values comparing to those of the fixed-GLMB and joint-CPHD filters. In addition, the errors in cardinality statistics using fixed-GLMB and *B*−GLMB are almost identical and better than error measured by joint-CPHD filter. The results of measuring errors between the set of true target tracks and that of the estimated tracks are given in Figure 4b. Once again, the effectiveness of the proposed method in reducing the errors in distances and locations of the target tracks is validated. The cardinality statistics for the *B*−GLMB filter, fixed-GLMB filter and joint-CPHD over 100 Monte Carlo run are given in Figure 5.

**Figure 5.** Cardinality tracking results.
