Three-Dimensional Multi-Target Tracking Using Dual-Orthogonal Baseline Interferometric Radar

Abstract

Multi-target tracking (MTT) generally needs either a Doppler radar network with spatially separated receivers or a single radar equipped with costly phased array antennas. However, Doppler radar networks have high computational complexity, attributed to the multiple receivers in the network. Moreover, array signal processing techniques for phased array radar also increase the computational burden on the processing unit. To resolve this issue, this paper investigates the problem of the detection and tracking of multiple targets in a three-dimensional (3D) Cartesian space based on range and 3D velocity measurements extracted from dual-orthogonal baseline interferometric radar. The contribution of this paper is twofold. First, a nonlinear 3D velocity measurement function, defining the relationship between the state of the target and 3D velocity measurements, is derived. Based on this measurement function, the design of the proposed algorithm includes the global nearest neighbor (GNN) technique for data association, an interacting multiple model estimator with a square-root cubature Kalman filter (IMM-SCKF) for state estimation, and a rule-based M/N logic for track management. Second, Monte Carlo simulation results for different multi-target scenarios are presented to demonstrate the performance of the algorithm in terms of track accuracy, computational complexity, and IMM mean model probabilities.

1. Introduction

Multi-target tracking (MTT) has received increased attention in real-time systems with a broad spectrum of applications, including aircraft tracking [1], surveillance [2], remote sensing [3,4], adaptive cruise control [5], robotics [6], biomedical engineering [7], image processing [8], and oceanography [9]. The main purpose of MTT algorithms is to identify the number of potential targets in the radar’s field of view (FOV) and to estimate their kinematic states from noisy radar measurements. Various algorithms to address the problem of MTT have been proposed in the literature [10,11,12]. MTT systems generally need either a Doppler radar network with spatially separated receivers or a single radar equipped with costly phased array antennas. The Doppler radar has the capability of determining the Doppler frequency shift of the target in the radar’s FoV, which is directly related to its radial velocity [13,14,15,16]. However, Doppler radar can provide the information of range and Doppler frequency shift only. Since it is not capable of extracting the direction-of-arrival (DOA), i.e., azimuth and elevation angles, information of the target, it is not possible to localize the target in the 3D Cartesian space using the range-Doppler data from single Doppler radar receiver only [17]. To resolve this issue, a distributed Doppler radar network with multiple sensors is typically employed, which can monitor the potential target of interest from different angles spatially. Various algorithms for MTT based on Doppler radar networks exist in the literature. Multistatic radar networks have become very popular in various real-time tracking applications over time [18,19,20]. An algorithm for target tracking, by monitoring the Doppler signals at multiple spatial points employing a network of four Doppler radars, was presented in [20]. The MTT algorithm using Doppler measurements from multistatic Doppler radar with an unknown probability of detection was proposed in [21]. A technique for tracking multiple targets by exploiting the range and Doppler information from multiple radar sensors was introduced in [22]. This technique, however, demands a multitude of iterations, which may cause the system to halt and generate latency. The MTT techniques presented in [17,23] are also based on multilateration of range-Doppler data from at least three radar sensors for localizing and tracking the targets without ambiguity. Multistatic radar systems for MTT by utilizing bistatic range and range-rate information from multiple radars were described in [24,25], which send data to a central station for estimating the spatial locations and velocities of the targets in the FoV. However, the integrated hardware due to multiple radar sensors in the network requires a high data transfer rate, making it difficult to realize real-time applications. The detection and tracking system using a single radar sensor demands a phased array antenna. However, a small-scale array of antennae offers poor angular resolution. Therefore, a costly large phased array is needed to provide angular information with better resolution. Furthermore, advanced array signal processing techniques such as the ESPRIT and MUSIC algorithms essentially increase the computational burden on the data processing unit caused by the complex matrix operations [26,27]. Furthermore, the MTT algorithm using a phased array antenna requires the number of receiving elements to be greater than the number of targets in the sensor’s FoV. To resolve this issue, in our previous work [28,29], we presented an algorithm for multiple targets’ detection and tracking in the 2D Cartesian space by measuring their range and 2D velocity (radial velocity and angular velocity) measurements using a dual-frequency frequency-modulated continuous wave (DF-FMCW) interferometric radar. Now, we propose to extend this concept to MTT in the 3D Cartesian space based on 3D velocity (radial velocity, azimuth angular velocity, and elevation angular velocity) measurements using a dual-orthogonal baseline DF-FMCW interferometric radar. The geometry of the dual-orthogonal baseline interferometric radar with a point source as the target is represented in Figure 1. The reference receiving antenna Rx1 is used to extract the initial ranges and radial velocities of the targets, and two orthogonal baselines with lengths $D_{12}$ and $D_{13}$ are used to measure the azimuth angular velocities and elevation angular velocities, respectively, of the targets in the radar’s FOV. $\theta$ and $\phi$ represent the azimuth and elevation angles of the target with the y-axis and xy-plane, respectively. $\rho$ is the range of the target relative to the observing radar in the 3D Cartesian space. The contribution of this paper is twofold:

• First, we present a mathematical model of the 3D velocity of multiple moving point sources and derive the nonlinear 3D velocity measurement function, which defines the relationship between the state of the target in terms of the 3D Cartesian space and the 3D velocity measurements extracted from the interferometric radar. Based on this measurement function, the design and implementation of the tracking algorithm are presented, which include (i) the GNN technique combined with the auction algorithm for data association, (ii) the IMM-SCKF estimator for state estimation, and (iii) the rule-based M/N logic for track management.
• Second, Monte Carlo simulation results with different multiple target scenarios are presented to validate the performance of the proposed algorithm in terms of track accuracy, computational complexity, and IMM mean model probabilities.

The layout of this paper is as follows. The mathematical formulation of the proposed MTT algorithm is presented in Section 2. Section 3 introduces the detailed design of the proposed MTT algorithm. The three main stages of the MTT algorithm including the GNN for data association, the IMM-SCKF estimator for state estimation, and the rule-based M/N logic for track management are described in Section 4, Section 5, Section 6 and Section 7. Section 8 presents the performance evaluation simulations. Finally, Section 9 provides the concluding remarks.

2. Mathematical Formulation of the Problem

The linear velocity of an arbitrary moving point source is a vector quantity. To completely localize a moving target in the 3D Cartesian space, it is necessary to measure the instantaneous 3D velocity vector ${\mathbf{v}}_i$ of the target, which is composed of radial velocity ${\mathbf{v}}_{r_{3D}}$, azimuth cross-radial velocity ${\mathbf{v}}_{{cr}_{\theta }}$, and elevation cross-radial velocity, i.e., ${\mathbf{v}}_{{cr}_{\phi }}$, ${\mathbf{v}}_i=\langle {\mathbf{v}}_{r_{3D}},{\mathbf{v}}_{{cr}_{\theta }},{\mathbf{v}}_{{cr}_{\phi }}\rangle$. The relationship between angular velocity and cross-radial velocity is defined as $\omega ={\mathbf{v}}_{cr}/\rho$.

2.1. The 3D Velocity of Point Sources

Assume an FMCW radar signal with carrier frequency $f_c$, bandwidth B, and sweep time T. The transmitted signal can be written as

$$S_T(t)=exp\{-j2\pi [f_c{t}_s+\frac{K}{2}{t}_s^2]\}$$

(1)

where $t_s=t-nT$ represents the time at the start of the nth sweep period and $K=B⁄T$ is the chirp rate. For an object initially present at a range of ${\rho }_0$ and moving with radial velocity $v_{r_{3D}}$ relative to the radar in the 3D Cartesian space, the signal reflected off of the object is the same as the transmitted signal, but delayed with a round trip time ${\tau }_{3D}$. Therefore, the received signal at antenna Rx1 for the FMCW radar can be written as

$$S_{R_1}(t)=exp\{-j2\pi [f_c(t_s-{\tau }_{3D})+\frac{K}{2}{(t_s-{\tau }_{3D})}^2]\}$$

(2)

where ${\tau }_{3D}=2\rho (t)/c$. The velocity of the object is considered to be slow enough that, during each pulse repetition interval, the object is generally expected to reside in the same range bin. Hence, the relative range of the target can be defined as

$$\rho (t)\approx \rho (nt)={\rho }_0+v_{r_{3D}}nT$$

(3)

According to the FMCW radar operating principle, the transmitted and received signals are mixed to obtain the beat frequency signal, which can be expressed as

$$\begin{array}{cc}\hfill S_{B_1}(t)& =S_T(t)S_{R_1}^{*}(t)\hfill \\ \hfill & \approx exp\{-j4\pi [\frac{K({\rho }_0+v_{r_{3D}}nT)t_s}{c}+\frac{{\rho }_0}{\lambda }+\frac{v_{r_{3D}}}{\lambda }nT]\}\hfill \end{array}$$

(4)

In the case of M number of point sources moving in the radar’s FOV, the beat frequency signal $S_{B_1}(t)$ can be expressed as

$$\begin{array}{c}\hfill {\displaystyle S_{B_1}(t)=\sum _{m=1}^{M}exp(-j4\pi [\frac{K({\rho }_{0m}+v_{{rm}_{3D}}nT)t_s}{c}+\frac{{\rho }_{0m}}{\lambda }+\frac{v_{{rm}_{3D}}}{\lambda }nT])}\end{array}$$

(5)

where ${\tau }_{m_{3D}}=2{\rho }_m⁄c$ for $m=1,2,\cdots ,M$. The relative radial motion of the object induces a Doppler frequency shift in the radar’s received signal. The radial velocity in terms of the Doppler frequency shift can be defined as

$$v_{{rm}_{3D}}=[\frac{f_{{dm}_{3D}}\lambda }{2}=\frac{c{f}_{{dm}_{3D}}}{2f_c},m=1,2,\cdots ,M]$$

(6)

The first baseline, constituting the receiving antennas Rx1 and Rx2 along the x-axis, is used to measure the azimuth angular velocity of the targets in the radar’s FOV. Considering the signal received by antenna Rx1 as the reference, the signal received by antenna Rx2, placed at a geometrical distance of $D_{12}$, can be represented as

$$\begin{array}{c}\hfill S_{R_2}(t)=exp\{-j2\pi [f_c(t_s-{\tau }_{3D}-{\tau }_{0_{\theta }})+\frac{K}{2}{(t_s-{\tau }_{3D}-{\tau }_{0_{\theta }})}^2]\}\end{array}$$

(7)

where ${\tau }_{0_{\theta }}=D_{12}sin\theta /c$. The second beat frequency signal $S_{B_2}(t)$ at receiving antenna Rx2 can be written as

$$\begin{array}{cc}\hfill S_{B_2}(t)& =S_T(t)S_{R_2}^{*}(t)\hfill \\ \hfill & \approx exp\{-j2\pi [\frac{2K({\rho }_0+v_{r_{3D}}nT)t_s}{c}+\frac{2v_{r_{3D}}nT}{\lambda }+\frac{2{\rho }_0}{\lambda }+(f_c+K{t}_s)\frac{D_{12}sin\theta }{c}]\}\hfill \end{array}$$

(8)

Following the interferometric radar principle, the beat frequency signals $S_{B_1}(t)$ and $S_{B_2}$ at the two receiving antennas are correlated to generate the interferometric output [30]:

$$\begin{array}{cc}\hfill y_{c_1}(t)& =S_{B_1}(t)S_{B_2}^{*}(t)\hfill \\ \hfill & =exp\{j2\pi (f_c+K{t}_s)\frac{D_{12}sin\theta }{c}\}\hfill \end{array}$$

(9)

In the case of M number of moving point sources, the interferometric output $y_{c_1}(t)$ represented by Equation (9) can be re-written as

$$\begin{array}{cc}\hfill y_{c_1}(t)& {\displaystyle =\sum _{m,b=1}^{M}exp(-j4\pi [\frac{K({\rho }_{0m}+v_{{rm}_{3D}}nT)t_s}{c}+\frac{{\rho }_{0m}}{\lambda }+\frac{v_{{rm}_{3D}}nT}{\lambda }])}\hfill \\ & exp(j2\pi [\frac{2K({\rho }_{0b}+v_{{rb}_{3D}}nT)t_s}{c}+\frac{2v_{{rb}_{3D}}nT}{\lambda }+\frac{2{\rho }_{0b}}{\lambda }+(f_c+K{t}_s)\frac{D_{12}sin{\theta }_b}{c}])\hfill \\ & {\displaystyle =\sum _{m=1}^{M}exp(j2\pi [f_c+K{t}_s]\frac{D_{12}sin{\theta }_m}{c})}\hfill \\ & {\displaystyle \sum _{m=1}^M\sum _{b=1,m\ne b}^{M}exp(-j4\pi [\frac{({\rho }_{0m}+v_{r{m}_{3D}}nT)-({\rho }_{0b}+v_{{rb}_{3D}}nT)}{c}K{t}_s}\hfill \\ & +\frac{(v_{{rm}_{3D}}nT-v_{{rb}_{3D}}nT)}{\lambda }+\frac{({\rho }_{0m}-{\rho }_{0b})}{\lambda }-(f_c+K{t}_s)\frac{D_{12}sin{\theta }_b}{2c}])\hfill \end{array}$$

(10)

As can be seen from Equation (10), the interferometric output is composed of two parts. The first part consists of M intra-correlation terms, generated by the angular velocities of moving objects, whereas the second part consists of $M(M-1)$ nuisance inter-correlation terms in which the radial and angular velocities of different moving objects are coupled together. In order to extract the angular velocities of the objects, these intermodulation terms need to be suppressed. Once these intermodulation terms are suppressed, azimuth angular velocity ${\omega }_{\theta }$ in terms of azimuth interferometric frequency shift $f_{\theta }$ is defined as

$${\omega }_{\theta m}=\frac{tial{\theta }_m}{tialt}=[\frac{f_{\theta m}{\lambda }_{t_s}}{D_{12}},m=1,2,\cdots ,M]$$

(11)

The second orthogonal baseline, constituting the receiving antennas Rx1 and Rx3 along the z-axis separated by geometrical distance $D_{13}$, is used to measure the elevation angular velocity of the targets in the radar’s FOV. Now, considering the signal received by the first receiving antenna Rx1 as the reference, the signal received by receiving antenna Rx3 can be written as

$$\begin{array}{cc}\hfill S_{R_3}(t)& =exp\{-j2\pi [f_c(t_s-{\tau }_{3D}-{\tau }_{0_{\phi }})+\frac{K}{2}{(t_s-{\tau }_{3D}-{\tau }_{0_{\phi }})}^2]\}\hfill \end{array}$$

(12)

where time delay ${\tau }_{0_{\phi }}=D_{13}sin\phi /c$. The transmitted and received signals represented by Equations (1) and (12), respectively, are mixed to generate the beat frequency signal. That is,

$$\begin{array}{cc}\hfill S_{B_3}(t)& =S_T(t)S_{R_3}^{*}(t)\hfill \\ \hfill & \approx exp\{-j2\pi [\frac{2K({\rho }_0+v_{r_{3D}}nT)t_s}{c}+\frac{2v_{r_{3D}}nT}{\lambda }+\frac{2{\rho }_0}{\lambda }+(f_c+K{t}_s)\frac{D_{13}sin\phi }{c}]\}\hfill \end{array}$$

(13)

The interferometric output, $y_{c_2}(t)$, of the second baseline along the z-axis is

$$\begin{array}{cc}\hfill y_{c_2}(t)& =S_{B_1}(t)S_{B_3}^{*}(t)\hfill \\ \hfill & =exp\{j2\pi (f_c+K{t}_s)\frac{D_{13}sin\phi }{c}\}\hfill \end{array}$$

(14)

In the case of M number of point sources moving in the radar’s FOV, the interferometric output $y_{c_2}(t)$ takes the following form.

$$\begin{array}{cc}\hfill y_{c_2}(t)& {\displaystyle =\sum _{m,b=1}^{M}exp(-j4\pi [\frac{K({\rho }_{0m}+v_{{rm}_{3D}}nT)t_s}{c}+\frac{{\rho }_{0m}}{\lambda }+\frac{v_{{rm}_{3D}}nT}{\lambda }])}\hfill \\ & exp(j2\pi [\frac{2K({\rho }_{0b}+v_{{rb}_{3D}}nT)t_s}{c}+\frac{2v_{{rb}_{3D}}nT}{\lambda }+\frac{2{\rho }_{0b}}{\lambda }+(f_c+K{t}_s)\frac{D_{13}sin{\phi }_b}{c}])\hfill \\ & {\displaystyle =\sum _{m=1}^{M}exp(j2\pi [f_c+K{t}_s]\frac{D_{13}sin{\phi }_m}{c})}\hfill \\ & {\displaystyle \sum _{m=1}^M\sum _{b=1,m\ne b}^{M}exp(-j4\pi [\frac{({\rho }_{0m}+v_{{rm}_{3D}}nT)-({\rho }_{0b}+v_{{rb}_{3D}}nT)}{c}K{t}_s}\hfill \\ & +\frac{(v_{{rm}_{3D}}nT-v_{{rb}_{3D}}nT)}{\lambda }+\frac{({\rho }_{0m}-{\rho }_{0b})}{\lambda }-(f_c+K{t}_s)\frac{D_{13}sin{\phi }_b}{2c}])\hfill \end{array}$$

(15)

Similar to the case of azimuth angular velocity measurement, the interferometric output $y_{c_2}(t)$ also consists of two parts, including M intra-correlation terms and $M(M-1)$ inter-correlation terms. The inter-correlation terms interfere with the extraction of the elevation angular velocities of the objects. The elevation angular velocity ${\omega }_{\phi }$, in terms of elevation interferometric frequency shift $f_{\phi }$, is defined as

$${\omega }_{\phi m}=\frac{tial{\phi }_m}{tialt}=[\frac{f_{\phi m}{\lambda }_{t_s}}{D_{13}},m=1,2,\cdots ,M]$$

(16)

2.2. Process Model

The process model describes the state transition between two consecutive time instants. Consider modeling the motion of the target by one of the i hypothesis models. These models can be represented by a set as $M_r:=\{1,2,\cdots ,r\}$. $M_{k-1}^j$ represents the event of model j being effective during time period $(t_{k-1},t_k]$. In this case, the target dynamic/motion for the jth hypothesis model can be described as

$${\mathbf{x}}_k={\mathbf{F}}_{k-1}^j{\mathbf{x}}_{k-1}+{\mathbf{v}}_{k-1}^j$$

(17)

where ${\mathbf{F}}_{k-1}^j$ represents the state transition matrix for motion model j being effective at time $k-1$, ${\mathbf{x}}_k$ represents the target’s state vector at time k, and ${\mathbf{v}}_{k-1}^j$ represents the process noise for the jth dynamic model, which is assumed to be independent and identically distributed (i.i.d.) zero-mean Gaussian noise with covariance ${\mathbf{Q}}_{k-1}^j$, such that ${\mathbf{v}}_{k-1}^j\sim \mathcal{N}(\mathbf{0},{\mathbf{Q}}_{k-1}^j)$. The target motion in the 3D Cartesian space is modeled with the nearly constant velocity (NCV) and nearly coordinated turn (NCT) process models. For uniform motion, the discrete-time NCV state dynamics combined with the DWNA model are represented by

$${\mathbf{x}}_k={\mathbf{F}}_{k-1}^{\mathrm{NCV}}{\mathbf{x}}_{k-1}+{\Gamma }_1{\mathbf{v}}_{1_{k-1}}$$

(18)

The state vector of the target with $n_x=6$ is defined as

$${\mathbf{x}}_k={[x_k,v_{x_k},y_k,v_{y_k},z_k,v_{z_k}]}^{\top }$$

(19)

where $x_k,y_k,z_k$ and $v_{x_k},v_{y_k},v_{z_k}$ denote the Cartesian coordinates and velocities of the object, respectively, at time instant k. $n_x$ represents the state vector dimension. The state transition matrix ${\mathbf{F}}^{\mathrm{NCV}}$ for the NCV model is

$$\begin{array}{c}\hfill {\mathbf{F}}_{k-1}^{\mathrm{NCV}}=\left[\begin{array}{cccccc}1& \Delta T& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& \Delta T& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& \Delta T\\ 0& 0& 0& 0& 0& 1\end{array}\right]\end{array}$$

(20)

The noise gain ${\Gamma }_1$ can be written as

$$\begin{array}{c}\hfill {\Gamma }_1=\left[\begin{array}{ccc}\frac{\Delta T^2}{2}& 0& 0\\ \Delta T& 0& 0\\ 0& \frac{\Delta T^2}{2}& 0\\ 0& \Delta T& 0\\ 0& 0& \frac{\Delta T^2}{2}\\ 0& 0& \Delta T\end{array}\right]\end{array}$$

(21)

The covariance of the process noise multiplied by gain ${\Gamma }_1$ is

$$\begin{array}{c}\hfill {\mathbf{Q}}_{1_{k-1}}=\left[\begin{array}{cccccc}\frac{\Delta T^4}{4}& \frac{\Delta T^3}{2}& 0& 0& 0& 0\\ \frac{\Delta T^3}{2}& \Delta T^2& 0& 0& 0& 0\\ 0& 0& \frac{\Delta T^4}{4}& \frac{\Delta T^3}{2}& 0& 0\\ 0& 0& \frac{\Delta T^3}{2}& \Delta T^2& 0& 0\\ 0& 0& 0& 0& \frac{\Delta T^4}{4}& \frac{\Delta T^3}{2}\\ 0& 0& 0& 0& \frac{\Delta T^3}{2}& \Delta T^2\end{array}\right]{\sigma }_{{\mathbf{v}}_1}^2\end{array}$$

(22)

where ${\sigma }_{{\mathbf{v}}_1}^2$ represents the variance of process noise ${\mathbf{v}}_1$. Here, ${\mathbf{v}}_1$ is a zero-mean Gaussian white noise used to model small accelerations, with an appropriate covariance ${\mathbf{Q}}_{1_{k-1}}$, which is a design parameter. For modeling the target maneuver, the discrete-time NCT state dynamics combined with the DWNA model in the 3D Cartesian coordinate system are represented by

$${\mathbf{x}}_k={\mathbf{F}}_{k-1}^{\mathrm{NCT}}{\mathbf{x}}_{k-1}+{\Gamma }_2{\mathbf{v}}_{2_{k-1}}$$

(23)

The state vector of the target augmented by turn rate $\Omega$, i.e., $n_x=7$, is defined as

$${\mathbf{x}}_k={[x_k,v_{x_k},y_k,v_{y_k},z_k,v_{z_k},{\Omega }_k]}^{\top }$$

(24)

The state transition matrix ${\mathbf{F}}^{\mathrm{NCT}}$ for the NCT model is

$$\begin{array}{c}\hfill {\mathbf{F}}_{k-1}^{\mathrm{NCT}}=\left[\begin{array}{ccccccc}1& \frac{sin({\Omega }_{k-1}\Delta T)}{{\Omega }_{k-1}}& 0& \frac{cos({\Omega }_{k-1}\Delta T)-1}{{\Omega }_{k-1}}& 0& 0& 0\\ 0& cos({\Omega }_{k-1}\Delta T)& 0& -sin({\Omega }_{k-1}\Delta T)& 0& 0& 0\\ 0& \frac{1-cos({\Omega }_k\Delta T)}{{\Omega }_{k-1}}& 1& \frac{sin({\Omega }_{k-1}\Delta T)}{{\Omega }_{k-1}}& 0& 0& 0\\ 0& sin({\Omega }_{k-1}\Delta T)& 0& cos({\Omega }_{k-1}\Delta T)& 0& 0& 0\\ 0& 0& 1& 0& 0& \Delta T& 0\\ 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]\end{array}$$

(25)

The noise gain ${\Gamma }_2$ for the DWNA model is defined as

$$\begin{array}{c}\hfill {\Gamma }_1=\left[\begin{array}{cccc}\frac{\Delta T^2}{2}& 0& 0& 0\\ \Delta T& 0& 0& 0\\ 0& \frac{\Delta T^2}{2}& 0& 0\\ 0& \Delta T& 0& 0\\ 0& 0& \frac{\Delta T^2}{2}& 0\\ 0& 0& \Delta T& 0\\ 0& 0& 0& \Delta T\end{array}\right]\end{array}$$

(26)

The covariance of the process noise multiplied by gain ${\Gamma }_2$ is

$$\begin{array}{c}\hfill {\mathbf{Q}}_{2_{k-1}}=\left[\begin{array}{ccccccc}\frac{\Delta T^4}{4}& \frac{\Delta T^3}{2}& 0& 0& 0& 0& 0\\ \frac{\Delta T^3}{2}& \Delta T^2& 0& 0& 0& 0& 0\\ 0& 0& \frac{\Delta T^4}{4}& \frac{T^3}{2}& 0& 0& 0\\ 0& 0& \frac{\Delta T^3}{2}& \Delta T^2& 0& 0& 0\\ 0& 0& 0& 0& \frac{\Delta T^4}{4}& \frac{\Delta T^3}{2}& 0\\ 0& 0& 0& 0& \frac{\Delta T^3}{2}& \Delta T^2& 0\\ 0& 0& 0& 0& 0& 0& \Delta T^2\end{array}\right]{\sigma }_{{\mathbf{v}}_2}^2\end{array}$$

(27)

where ${\sigma }_{{\mathbf{v}}_2}^2$ represents the variance of process noise ${\mathbf{v}}_2$.

2.3. Nonlinear Measurement Model

The nonlinear measurement model defining the relationship between the 3D velocity measurements received from the interferometric radar and the state of the target is defined as

$${\mathbf{z}}_k=\mathbf{h}({\mathbf{x}}_k)+{\mathbf{w}}_k$$

(28)

where $\mathbf{h}({\mathbf{x}}_k)$ is the nonlinear 3D velocity measurement function and ${\mathbf{z}}_k={[v_{r_{3D,k}},{\omega }_{{\theta }_k},{\omega }_{{\phi }_k}]}^{\top }$ with $n_z=3$. Here, $n_z$ denotes the dimension of the measurement vector. ${\mathbf{w}}_k$ is assumed to be an i.i.d. zero-mean Gaussian measurement noise with covariance ${\mathbf{R}}_k$, i.e., ${\mathbf{w}}_k\sim \mathcal{N}(\mathbf{0},{\mathbf{R}}_k)$ and $E\{{\mathbf{v}}_k{\mathbf{w}}_k^T\}=\mathbf{0}$.

Derivation of 3D Velocity Measurement Function

The range of the target relative to the observing radar in the 3D Cartesian space is expressed as

$${\rho }_k=\sqrt{{(x_k-x_s)}^2+{(y_k-y_s)}^2+{(z_k-z_s)}^2}$$

(29)

By taking the time derivative of Equation (29), the radial velocity of the target in the 3D Cartesian space can be determined. That is,

$${\dot{\rho }}_k=\frac{d\rho }{dt}=\frac{(x_k-x_s){\dot{x}}_k+(y_k-y_s){\dot{y}}_k+(z_k-z_s){\dot{z}}_k}{\sqrt{{(x_k-x_s)}^2+{(y_k-y_s)}^2+{(z_k-z_s)}^2}}$$

(30)

As we know that $v_{r_{3D,k}}={\dot{\rho }}_k$, $v_{x_k}={\dot{x}}_k$, $v_{y_k}={\dot{y}}_k$, and $v_{z_k}={\dot{z}}_k$, the expression for the radial velocity can be written as

$$v_{r_{3D,k}}={\dot{\rho }}_k=\frac{(x_k-x_s)v_{x_k}+(y_k-y_s)v_{y_k}+(z_k-z_s)v_{z_k}}{\sqrt{{(x_k-x_s)}^2+{(y_k-y_s)}^2+{(z_k-z_s)}^2}}$$

(31)

According to [29], the azimuth angular velocity is determined by

$${\omega }_{{\theta }_k}={\dot{\theta }}_k=\left[\frac{(y_k-y_s)v_{x_k}-(x_k-x_s)v_{y_k}}{{(x_k-x_s)}^2+{(y_k-y_s)}^2}\right]$$

(32)

The elevation angle ${\phi }_k$ relative to the observing radar is defined as

$$\phi ={tan}^{-1}\frac{(z_k-z_s)}{R_k}$$

(33)

Now, by re-arranging and taking the time derivative of Equation (33),

$$({sec}^2{\phi }_k)\dot{\phi }=\frac{R_k{\dot{z}}_k-(z_k-z_s){\dot{R}}_k}{{R_k}^2}$$

(34)

$${\omega }_{{\phi }_k}=\left[\frac{R_k^2{v}_{z_k}-(z_k-z_s)\{(x_k-x_s)v_{x_k}+(y_k-y_s)v_{y_k}\}}{R_k^3}\right]({cos}^2{\phi }_k)$$

(35)

where ${\omega }_{{\phi }_k}={\dot{\phi }}_k$. Since we know that

$$cos{\phi }_k=\frac{R_k}{{\rho }_k}$$

(36)

Equation (35) can be re-written as

$${\omega }_{{\phi }_k}={\dot{\phi }}_k=\left[\frac{R_k^2{v}_{z_k}-(z_k-z_s)\{(x_k-x_s)v_{x_k}+(y_k-y_s)v_{y_k}\}}{R_k^3}\right]\left[\frac{R^2}{{\rho }_k^2}\right]$$

(37)

$${\omega }_{{\phi }_k}=\left[\frac{R_k^2{v}_{z_k}-(x_k-x_s)(z_k-z_s)v_{x_k}-(y_k-y_s)(z_k-z_s)v_{y_k}}{{\rho }_k^2{R}_k}\right]$$

(38)

Proceeding from Equations (31), (32) and (38), the nonlinear 3D velocity measurement function can be defined as

$$\begin{array}{c}\hfill \mathbf{h}(x_k)=\left[\begin{array}{c}\frac{(x_k-x_s)v_{x_k}+(y_k-y_s)v_{y_k}+(z_k-z_s)v_{z_k}}{\sqrt{{(x_k-x_s)}^2+{(y_k-y_s)}^2+{(z_k-z_s)}^2}}\\ \\ \frac{(y_k-y_s)v_{x_k}-(x_k-x_s)v_{y_k}}{{(x_k-x_s)}^2+{(y_k-y_s)}^2}\\ \\ \frac{R_k^2{v}_{z_k}-(x_k-x_s)(z_k-z_s)v_{x_k}-(y_k-y_s)(z_k-z_s)v_{y_k}}{{\rho }_k^2{R}_k}\end{array}\right]\end{array}$$

(39)

Hence, the nonlinear measurement model defined by Equation (28) takes on the following form.

$$\begin{array}{c}\hfill {\mathbf{z}}_k=\left[\begin{array}{c}\frac{(x_k-x_s)v_{x_k}+(y_k-y_s)v_{y_k}+(z_k-z_s)v_{z_k}}{\sqrt{{(x_k-x_s)}^2+{(y_k-y_s)}^2+{(z_k-z_s)}^2}}\\ \\ \frac{(y_k-y_s)v_{x_k}-(x_k-x_s)v_{y_k}}{{(x_k-x_s)}^2+{(y_k-y_s)}^2}\\ \\ \frac{R_k^2{v}_{z_k}-(x_k-x_s)(z_k-z_s)v_{x_k}-(y_k-y_s)(z_k-z_s)v_{y_k}}{{\rho }_k^2{R}_k}\end{array}\right]+{\mathbf{w}}_k\end{array}$$

(40)

The objective of the proposed MTT algorithm is to estimate the state ${\widehat{\mathbf{x}}}_{k|k}=E\{{\mathbf{x}}_k|{\mathbf{z}}_k\}$ and error covariance ${\mathbf{P}}_{k|k}=E\{[{\mathbf{x}}_k-{\widehat{\mathbf{x}}}_{k|k}]{[{\mathbf{x}}_k-{\widehat{\mathbf{x}}}_{k|k}]}^{\top }|{\mathbf{z}}_k\}$ of each target in the radar’s FOV by exploiting the derived 3D velocity measurement function.

3. Design of the Proposed 3D MTT Algorithm

The detailed design of the proposed MTT algorithm based on 3D velocity measurements is delineated in Figure 2. The dual-orthogonal baseline DF-FMCW interferometric radar used to measure the 3D velocities of multiple moving objects has the capability to transmit FMCW signals with two different carrier frequencies, such that $f_{c_1}=6$ GHz and $f_{c_2}=24$ GHz. It consists of one transmitting antenna Tx and three receiving antennas Rx1, Rx2, and Rx3. The transmitting antenna Tx and receiving antenna Rx1 can switch between two operating frequencies $f_{c_1}$ and $f_{c_2}$ via RF switches, depending on the mode of operation. Based on the carrier frequency $f_{c_1}$, the lengths of the two orthogonal baselines were set as $D_{12}=D_{13}=60{\lambda }_{c_1}=3$ m. The detection and extraction of radial velocity measurements were performed at a higher carrier frequency $f_{c_2}$, whereas the azimuth and elevation angular velocity measurements were obtained at a lower carrier frequency $f_{c_1}$, thereby suppressing the intermodulation terms in the interferometric response. The interferometric frequencies’ information was preserved by increasing the baseline lengths $D_{12}$ and $D_{13}$ simultaneously [28]. First of all, the radar observes the region of interest (ROI) operating at carrier frequency $f_{c_2}$ and performs 2D-FFT on the received data for the detection of the potential targets by obtaining the range–radial velocity map. Once the targets are identified, STFT is applied to the received data, which provides the time-varying Doppler spectrogram of the targets in the FOV. Following Equation (6), the radial velocities of the targets are extracted from the Doppler spectrogram [29]. Then, for the azimuth and elevation angular velocity measurements at carrier frequency $f_{c_1}$, STFT is performed on the two interferometric outputs $y_{c_1}(t)$ and $y_{c_2}(t)$, respectively. The time-varying azimuth and elevation interferometric spectrograms thus obtained are used to calculate the azimuth and elevation angular velocities of the targets following Equations (11) and (16), respectively. In order to combine the 3D velocity measurements of each object, 2D FFT is performed on the three beat frequency signals at antennas Rx1, Rx2, and Rx3 in the interferometric mode corresponding to the carrier frequency $f_{c_1}$. The following relationships hold true in the interferometric mode.

$${\rho }_2-{\rho }_1=D_{12}sin\theta$$

(41)

$${\rho }_3-{\rho }_1=D_{13}sin\phi$$

(42)

where ${\rho }_1$, ${\rho }_2$, and ${\rho }_3$ represent the ranges of the object relative to antennas Rx1, Rx2, and Rx3, respectively. By taking the time derivative of Equations (41) and (42), the relationship between the radial and angular velocities of the object can be written as

$$v_{r_{3D,2}}-v_{r_{3D,1}}=D_{12}{\omega }_{\theta }cos\theta$$

(43)

$$v_{r_{3D,3}}-v_{r_{3D,1}}=D_{13}{\omega }_{\phi }cos\phi$$

(44)

where $v_{r_{3D,1}}$, $v_{r_{3D,2}}$, and $v_{r_{3D,2}}$ represent the radial velocity of the object at antennas Rx1, Rx2, and Rx3, respectively. The initial range measurements obtained by performing 2D FFT on the beat frequency signals are used to calculate the initial azimuth and elevation angles of the objects following Equations (41) and (42). Then, based on Equations (43) and (44), the initial values of the radial velocities along with the azimuth and elevation angle measurements are used to combine the 3D velocities of each target distinctly. Measurement-to-track data association and target state estimation are the two fundamental elements of the MTT algorithm. Moreover, the track management is also a crucial element of the MTT algorithm, which maintains the two major data structures, namely tentative and confirmed track lists. The track management unit handles the initiation of new target tracks, the confirmation of tentative tracks, and the deletion of tentative and confirmed tracks based on some predefined criteria. After the formation of the tentative and confirmed track lists from previous scans of the data, the new 3D velocity measurements received from the interferometric radar are assessed for association with existing target tracks or for initializing new tentative tracks. Based on the 3D velocity measurement function defined by Equation (39), the global nearest neighbor (GNN) method combined with the auction algorithm is used for measurement-to-track data association. The unassociated measurements are then tested for association with already existing tentative tracks. The measurements still unassociated with any existing tracks are then used to initiate new tracks. The 3D velocity measurements associated with the existing tracks are utilized in the filter measurement update step for the target’s state estimation. Since Equation (39) clearly demonstrates the nonlinear relation between the state of the target and the radar measurements, the proposed MTT algorithm utilizes the nonlinear IMM-SCKF estimator. Based on the 3D velocity measurement function, IMM-SCKF estimator is used to estimate the states of non-maneuvering and maneuvering targets. Different blocks of the MTT algorithm are thoroughly discussed in the following sections for better understanding of the algorithm.

4. Global Nearest Neighbor Data Association Method

Various data association techniques exist in the literature for the MTT algorithms, which include the NN, GNN, JPDA, and MHT algorithms. However, the GNN is the most commonly used optimal and robust technique for practical tracking systems. The GNN algorithm propagates a single global hypothesis and addresses the measurement-to-track association in the form of an optimal assignment problem [29,31,32]. Based on the 3D velocity measurement function expressed by Equation (39), the data association process follows three steps:

• Ellipsoid gating;
• GNN-algorithm-based assignment matrix formation;
• Auction algorithm for the optimal solution of the assignment matrix.

Gating is a technique for eliminating unlikely observation-to-track pairings. For a measurement to satisfy the ellipsoid gating relationship of a track, the norm of the measurement residual vector $d_k^2$ must fulfill the following criterion:

$$d_k^2={\tilde{\mathbf{v}}}_k^{\top }{\mathbf{S}}_{Z,k}^{-1}{\tilde{\mathbf{v}}}_k\le G$$

(45)

where $\tilde{\mathbf{v}}$ is the measurement residual vector, ${\mathbf{S}}_{Z,k}$ is the measurement residual covariance matrix at time instant k, and G represents the gate size. Based on the 3D velocity measurement function, the measurement residual vector can be expressed as

$${\tilde{\mathbf{v}}}_k={\mathbf{z}}_k-\left[\begin{array}{c}\frac{(x_k-x_s)v_{x_k}+(y_k-y_s)v_{y_k}+(z_k-z_s)v_{z_k}}{\sqrt{{(x_k-x_s)}^2+{(y_k-y_s)}^2+{(z_k-z_s)}^2}}\\ \\ \frac{(y_k-y_s)v_{x_k}-(x_k-x_s)v_{y_k}}{{(x_k-x_s)}^2+{(y_k-y_s)}^2}\\ \\ \frac{R_k^2{v}_{z_k}-(x_k-x_s)(z_k-z_s)v_{x_k}-(y_k-y_s)(z_k-z_s)v_{y_k}}{{\rho }_k^2{R}_k}\end{array}\right]{\widehat{\mathbf{x}}}_{k|k-1}$$

(46)

Generally, $d^2$ is assumed to have a chi distribution (${\chi }_{n_z}^2$) for correct measurement-to-track association, where $n_z$ represents the dimension of the measurement vector. The relationship between the gate size G and the probability of the 3D measurement vector falling inside the gate $P_G(n_z)$ is expressed as [11]

$$P_G(3)=2gc(\sqrt{G^{3D}})-\sqrt{\frac{2G^{3D}}{\pi }}exp(\frac{-G^{3D}}{2})$$

(47)

where $gc$ denotes standard Gaussian probability integral. Following the ellipsoid gating, the GNN algorithm is used to resolve measurement-to-track association conflicts by forming and solving a 2D assignment matrix [11]. For the formation of the assignment matrix, the rows are occupied by 3D velocity measurements, and the first columns of the assignment matrix are filled with the existing tracks. The rest of the columns are filled with new track hypotheses. The entries of the assignment matrix are expressed in the form of the score gain. That is,

$${\alpha }_{ij}=G-d_{ij}^2$$

(48)

where ${\alpha }_{ij}$ represents the margin by which the statistical distance passed the gate. The entries corresponding to the assignment of a measurement to a new track were set as zero. The non-allowed measurements are represented by −∞. The optimal solution is the set of assignments that produces the maximum score gain [11]. In the literature, various algorithms exist to solve the assignment problem, such as the Hungarian algorithm, the Munkres algorithm [33], the Jonker–Volgenant–Castanon (JVC) algorithm, and the auction algorithm [11]. However, the auction algorithm is presently the most efficient algorithm for the assignment problem, which has replaced the Munkres algorithm [12]. Therefore, in the current work, the GNN method with the auction algorithm was used for measurement-to-track data association.

5. IMM-SCKF Estimator for State Estimation

The main function of the estimation filter in the MTT algorithm is to estimate the state of the targets detected from noise-corrupted radar measurements [34]. The measurements allocated to the tracks present in the tentative or confirmed track lists during the data association step are utilized in the measurement update stage of the filtering procedure [24]. Referring to Equation (39), it is clear that there exists nonlinear relation between the 3D velocity measurements received from the interferometric radar and the state of the target. Therefore, a nonlinear filter is required for estimating the state of the target [31,32,35].

5.1. Square-Root Cubature Kalman Filter

Based on the properties of symmetry and positive semi-definiteness, the nonlinear SCKF is used for target state estimation. It is numerically more accurate, stable, and computationally efficient compared to the extended Kalman filter (EKF) [36], unscented Kalman filter (UKF) [37], and CKF [38]. The SCKF algorithm is outlined below [39];

5.1.1. Time Update

Starting at time instant k, assume that the posterior density with the square-root of error covariance matrix is given.

$$p({\mathbf{x}}_{k-1}|{\mathbf{Z}}_1^{k-1})=\mathcal{N}({\widehat{\mathbf{x}}}_{k-1|k-1}|{\mathbf{S}}_{k-1|k-1})$$

(49)

$${\mathbf{P}}_{k-1|k-1}={\mathbf{S}}_{k-1|k-1}{\mathbf{S}}_{k-1|k-1}^{\top }$$

(50)

where ${\mathbf{P}}_{k-1|k-1}$ denotes the error covariance matrix, ${\mathbf{S}}_{k-1|k-1}$ is the square-root of the error covariance matrix ${\mathbf{P}}_{k-1|k-1}$, and ${\mathbf{Z}}_1^{k-1}$ represents the set of measurements from the start to time instant $k-1$:

• Compute the cubature points ${\mathbf{X}}_i$:$(i=1,2,\cdots ,m$, where m = 2$n_x$):

  $${\mathbf{X}}_{i,k-1|k-1}={\mathbf{S}}_{k-1|k-1}{\xi }_i+{\widehat{\mathbf{x}}}_{k-1|k-1}$$

  (51)

  $${\xi }_i=\sqrt{\frac{m}{2}}{[1]}_i$$

  (52)

  where ${\xi }_i$ denotes the cubature points weights.
• Compute the cubature points propagated $(i=1,2,\cdots ,m)$.

  $${\mathbf{X}}_{i,k|k-1}^{*}=\mathbf{f}({\mathbf{X}}_{i,k-1|k-1})$$

  (53)

  where $\mathbf{f}(·)$ denotes the nonlinear function of the dynamic state equation.
• Calculate the predicted state.

  $${\widehat{\mathbf{x}}}_{k|k-1}=\frac{1}{m}\sum _{m=1}^M{\mathbf{X}}_{i,k|k-1}^{*}$$

  (54)
• Calculate the square-root predicted error covariance.

  $${\mathbf{S}}_{k|k-1}=\mathbf{Tria}([{\mathcal{X}}_{k|k-1}^{*}{\mathbf{S}}_{Q,k-1}])$$

  (55)

  $${\mathbf{Q}}_{k-1}={\mathbf{S}}_{Q,k-1}{\mathbf{S}}_{Q,k-1}^{\top }$$

  (56)

  $$\begin{array}{cc}\hfill {\mathcal{X}}_{k|k-1}^{*}& =\frac{1}{\sqrt{m}}[{\mathbf{X}}_{1,k|k-1}^{*}-{\widehat{\mathbf{x}}}_{k|k-1}\hfill \\ & {\mathbf{X}}_{2,k|k-1}^{*}-{\widehat{\mathbf{x}}}_{k|k-1}\cdots {\mathbf{X}}_{m,k|k-1}^{*}-{\widehat{\mathbf{x}}}_{k|k-1}]\hfill \end{array}$$

  (57)

  where ${\mathbf{S}}_{Q,k-1}$ is the square-root of process noise covariance ${\mathbf{Q}}_{k-1}$ and ${\mathcal{X}}_{k|k-1}^{*}$ is the weighted-centered matrix. $\mathbf{Tria}$ represents the QR-decomposition triangularization algorithm.

5.1.2. Measurement Update

• Compute the cubature points $(i=1,2,\cdots ,m)$.

  $${\mathbf{X}}_{i,k|k-1}={\mathbf{S}}_{k|k-1}{\xi }_i+{\widehat{\mathbf{x}}}_{k|k-1}$$

  (58)
• Compute the cubature points propagated $(i=1,2,...,m)$.

  $${\mathbf{Z}}_{i,k|k-1}^{*}=\mathbf{h}({\mathbf{X}}_{i,k|k-1})$$

  (59)

  where $\mathbf{h}(·)$ denotes the nonlinear 3D velocity measurement function defined by Equation (39).
• Calculate the predicted measurement.

  $${\widehat{\mathbf{z}}}_{k|k-1}=\frac{1}{m}\sum _{m=1}^M{\mathbf{Z}}_{i,k|k-1}^{*}$$

  (60)
• Compute the square-root innovation covariance matrix.

  $${\mathbf{S}}_{zz,k|k-1}=\mathbf{Tria}([{\mathfrak{J}}_{k|k-1}{\mathbf{S}}_{R,k}])$$

  (61)

  $${\mathbf{R}}_k={\mathbf{S}}_{R,k}{\mathbf{S}}_{R,k}^{\top }$$

  (62)

  $$\begin{array}{cc}\hfill {\mathfrak{J}}_{k|k-1}& =\frac{1}{\sqrt{m}}[{\mathbf{Z}}_{1,k|k-1}^{*}-{\widehat{\mathbf{z}}}_{k|k-1}\hfill \\ & {\mathbf{Z}}_{2,k|k-1}^{*}-{\widehat{\mathbf{x}}}_{k|k-1}\cdots {\mathbf{Z}}_{m,k|k-1}^{*}-{\widehat{\mathbf{z}}}_{k|k-1}]\hfill \end{array}$$

  (63)

  where ${\mathbf{S}}_{R,k}$ is the square-root of measurement noise covariance ${\mathbf{R}}_{k-1}$ and ${\mathfrak{J}}_{k|k-1}$ is the weighted-centered matrix.
• Compute the cross-covariance matrix.

  $${\mathbf{P}}_{xz,k|k-1}={\mathcal{X}}_{k|k-1}{\mathfrak{J}}_{k|k-1}^{\top }$$

  (64)

  $$\begin{array}{cc}\hfill {\mathcal{X}}_{k|k-1}& =\frac{1}{\sqrt{m}}[{\mathbf{X}}_{1,k|k-1}-{\widehat{\mathbf{x}}}_{k|k-1}\hfill \\ & {\mathbf{X}}_{2,k|k-1}-{\widehat{\mathbf{x}}}_{k|k-1}\cdots {\mathbf{X}}_{m,k|k-1}-{\widehat{\mathbf{x}}}_{k|k-1}]\hfill \end{array}$$

  (65)
• Calculate the square-root cubature Kalman gain.

  $${\mathbf{W}}_k=({\mathbf{P}}_{xz,k|k-1}/{\mathbf{S}}_{zz,k|k-1}^{\top })/{\mathbf{S}}_{zz,k|k-1}$$

  (66)
• Compute the updated state estimate.

  $${\widehat{\mathbf{x}}}_{k|k}={\widehat{\mathbf{x}}}_{k|k-1}+{\mathbf{W}}_k({\mathbf{z}}_k-{\widehat{\mathbf{z}}}_{k|k-1})$$

  (67)
• Calculate the square-root estimated error covariance.

  $${\mathbf{S}}_{k|k}=\mathbf{Tria}([{\mathcal{X}}_{k|k-1}-{\mathbf{W}}_k{\mathfrak{J}}_{k|k-1}{\mathbf{W}}_k{\mathbf{S}}_{R,k}])$$

  (68)

5.2. Interacting Multiple Model Estimator

Since the moving target does not necessarily follow a single dynamic model, modern tracking systems typically use the IMM filter [40] for maneuvering target tracking [29]. If the target dynamics are characterized by multiple switching models $(r>1)$, the posterior density of the target state vector is a mixture density [11,12,25,31]. The objective of the IMM filter is to combine all the mixture components into a single Gaussian distribution in a manner that the first and second moments are matched [24,25,29]. At the current time instant k, the IMM algorithm begins with r model-conditioned state estimates ${\widehat{\mathbf{x}}}_{k-1|k-1}^i$, square-root error covariance matrix ${\mathbf{S}}_{k-1|k-1}^i$, and model probabilities ${\mu }_{k-1|k-1}^i$ at instant $k-1$[24]. Upon receiving the current measurement ${\mathbf{z}}_k$, the IMM algorithm follows the steps outlined below [29,31,41]:

• Computing mixing probabilities:

  $${\overline{c}}^j=\sum _{i=1}^r{p}_{ij}{\mu }_{k-1|k-1}^i$$

  (69)

  $${\mu }_{k-1|k-1}^{(i,j)}=P\{M_{k-1}^i|M_k^j,{\mathbf{Z}}_1^{k-1}\}=\frac{1}{{\overline{c}}^j}{p}_{ij}{\mu }_{k|k-1}^i$$

  (70)

  where ${\mu }_{k-1|k-1}^{(i,j)}$ represents the model probability that model $M^i$ was effective at time instant $k-1$ provided that $M^j$ is effective at time instant k conditioned on ${\mathbf{Z}}_1^{k-1}$. Here, ${\mathbf{Z}}_1^{k-1}$ represents the measurement history from start to instant $k-1$. $p_{ij}$ denotes the elements of model transition probabilities matrix, ${\mu }_{k-1|k-1}^i$ is the model probability for ith SCKF, and ${\mu }_{k-1|k-1}^{(i,j)}$ is mixture probability for $i,j=1,\cdots ,r$.
• Interaction/mixing of state estimates:

  $${\widehat{\mathbf{x}}}_{k-1|k-1}^{0j}=\sum _{i=1}^r{\widehat{\mathbf{x}}}_{k-1|k-1}^i{\mu }_{k-1|k-1}^{(i,j)}$$

  (71)

  $${\mathbf{P}}_{k-1|k-1}^{0j}=\sum _{i=1}^r{\mu }_{k-1|k-1}^{(i,j)}\{{\mathbf{S}}_{k-1|k-1}^i{({\mathbf{S}}_{k-1|k-1}^i)}^{\top }+[{\widehat{\mathbf{x}}}_{k-1|k-1}^i-{\widehat{\mathbf{x}}}_{k-1|k-1}^{0j}]{[{\widehat{\mathbf{x}}}_{k-1|k-1}^i-{\widehat{\mathbf{x}}}_{k-1|k-1}^{0j}]}^{\top }\}$$

  (72)

  where ${\widehat{\mathbf{x}}}_{k-1|k-1}^{0j}$ represents the mixed initial state estimate for the jth SCKF, ${\mathbf{P}}_{k-1|k-1}^{0j}$ represents the error covariance corresponding to ${\widehat{\mathbf{x}}}_{k-1|k-1}^{0j}$, ${\widehat{\mathbf{x}}}_{k-1|k-1}^i$ represents the ith SCKF state estimate, and ${\mathbf{S}}_{k-1|k-1}^i$ is the square-root covariance matrix for ${\widehat{\mathbf{x}}}_{k-1|k-1}^i$i [41].
• Updating state estimates: The initial condition state estimate ${\widehat{\mathbf{x}}}_{k-1|k-1}^{0j}$ and error covariance matrix ${\mathbf{P}}_{k-1|k-1}^{0j}$ for each filter model are used to compute the updated state estimate ${\widehat{\mathbf{x}}}_{k|k}^j$ and square-root error covariance ${\mathbf{S}}_{k|k}^j$ for the SCKF model.
• Computing the model likelihood:

  $$\begin{array}{cc}{\tilde{\mathbf{v}}}_k^j\hfill & ={\mathbf{z}}_k-\mathbf{h}({\widehat{\mathbf{x}}}_{k|k-1}^j)\hfill \\ \hfill & ={\mathbf{z}}_k-\left[\begin{array}{c}\frac{(x_k-x_s)v_{x_k}+(y_k-y_s)v_{y_k}+(z_k-z_s)v_{z_k}}{\sqrt{{(x_k-x_s)}^2+{(y_k-y_s)}^2+{(z_k-z_s)}^2}}\\ \frac{(y_k-y_s)v_{x_k}-(x_k-x_s)v_{y_k}}{{(x_k-x_s)}^2+{(y_k-y_s)}^2}\\ \frac{R_k^2{v}_{z_k}-(x_k-x_s)(z_k-z_s)v_{x_k}-(y_k-y_s)(z_k-z_s)v_{y_k}}{{\rho }_k^2{R}_k}\end{array}\right]{\widehat{\mathbf{x}}}_{k|k-1}^j\hfill \end{array}$$

  (73)

  $${\wedge }_k^j=\frac{1}{\sqrt{|2\pi {\mathbf{S}}_{Z,k}^j|}}exp\{-0.5{[{\tilde{\mathbf{v}}}_k^j]}^T{[{\mathbf{S}}_{Z,k}^j]}^{-1}[{\tilde{\mathbf{v}}}_k^j]\}$$

  (74)

  $${\mathbf{S}}_{Z,k}^j={\mathbf{A}}_k^j{({\mathbf{A}}_k^j)}^{\top }$$

  (75)

  where ${\wedge }_k^j$ is the likelihood function and ${\mathbf{A}}_k^j$ represents the square-root measurement residual covariance matrix for the jth SCKF model.
• Updating the model probability:

  $$c=\sum _{j=1}^r{\wedge }_k^j{\overline{c}}^j$$

  (76)

  $${\mu }_{k|k}^j=P\{M_k^j|{\mathbf{Z}}_1^k\}=\frac{1}{c}{\wedge }_k^j{\overline{c}}^j$$

  (77)
• Combining the state mean and covariance estimates for the output:

  $${\widehat{\mathbf{x}}}_{k|k}=\sum _{j=1}^r{\widehat{\mathbf{x}}}_{k|k}^j{\mu }_{k|k}^j$$

  (78)

  $$\begin{array}{cc}\hfill {\mathbf{P}}_{k|k}& =\sum _{j=1}^r{\mu }_{k|k}^j\{{\mathbf{S}}_{k|k}^j{({\mathbf{S}}_{k|k}^j)}^{\top }\hfill \\ & +[{\widehat{\mathbf{x}}}_{k|k}-{\widehat{\mathbf{x}}}_{k|k}^j]{[{\widehat{\mathbf{x}}}_{k|k}-{\widehat{\mathbf{x}}}_{k|k}^j]}^{\top }\}\hfill \end{array}$$

  (79)

The proposed MTT algorithm utilizes the IMM-SCKF estimator with the 3D velocity measurement function for the confirmed tracks’ state estimation. For the tracks not receiving measurements during the data association process, the predicted state and covariance estimates become the measurement-updated state and covariance estimates. Since only a few tentative tracks fulfill the criteria to be inserted into the confirmed tracks’ list, a simple NCV-SCKF estimator is used for tentative track state estimation. Employing the IMM-SCKF estimator for the targets in the tentative tracks’ list will increase the computational load and complexity for the data processing unit.

6. Filter Initialization

Filter initialization is an extremely crucial problem in tracking applications, specifically for nonlinear systems, as nonlinear estimation filters generally depend on approximation theory [24,25]. The nonlinear measurement function in Equation (39), defining the relationship between 3D velocity measurements and the state of the target, makes finding the analytical solution impossible for filter initialization. Therefore, the initial range measurement and initial azimuth and elevation angles determined from Equations (41) and (42), respectively, are used to extract the 3D Cartesian position $(x_0,y_0,z_0)$ of the target. That is,

$$x_0=\rho cos\phi sin\theta$$

(80)

$$y_0=\rho cos\phi cos\theta$$

(81)

$$z_0=\rho sin\phi$$

(82)

The velocity of the target is calculated by integrating the position estimate at two consecutive time instants, i.e., $v_{x_0}=(x_1-x_0)⁄\Delta T$, $v_{y_0}=(y_1-y_0)⁄\Delta T$, and $v_{z_0}=(z_1-z_0)⁄\Delta T$. Since the azimuth and elevation angles calculated are highly sensitive to the range measurement error, $(R,\theta ,\phi )$ measurements are not reliable for target tracking in the 3D Cartesian space. These measurements are only utilized for filter initialization, and tracking is performed based on the 3D velocity information obtained from the interferometric radar.

7. Rule-Based M/N Logic for Track Management

Since the GNN algorithm is incapable of addressing the track management issue automatically, the rule-based M/N logic is used to handle the appearance and disappearance of the targets in the radar’s FOV. The track management involves new track initiation, tentative track confirmation, and tentative/confirmed track deletion. The measurement-to-track data association procedure is followed by the track management step, which is responsible for managing both the tentative tracks’ list and confirmed tracks’ list on the basis of predefined rules. Each track list contains the information of state estimate vector ${\widehat{\mathbf{x}}}_{k|k}$, square-root error covariance estimate matrix ${\mathbf{S}}_{k|k}$, measurement residual covariance matrix ${\mathbf{S}}_{Z,k}$, “Hit” counter, and “Miss” counter. Upon receiving a new set of 3D velocity measurements from the interferometric radar, these measurements are first tested for association with the confirmed tracks. If a confirmed track does not get associated with a measurement during this step, its “Hit” counter is decremented and its “Miss” counter is incremented. The measurements not associated with the confirmed tracks are tested for data association with existing tentative tracks. The “Hit” counters of the tentative tracks receiving measurements during the data association procedure get incremented, whereas their “Miss” counters remain unchanged. On the other hand, if a tentative track does not receive any measurement during the data association step, its “Hit” counter remains unchanged, whereas its “Miss” counter gets incremented. The still unassociated measurements are utilized by the track management unit to initiate new tentative tracks. The “Hit” and “Miss” counters of new tentative tracks are set to 1 and 0, respectively. The tentative tracks satisfying the 2/2 and 2/3 rules are inserted in the confirmed tracks’ list. This means that a tentative track receiving measurements during the first two consecutive scans of the data and then receiving measurements at least twice during the next three consecutive scans is eligible to be inserted in the confirmed tracks’ list. At this stage, the ”Hit” and “Miss” counters of the confirmed tracks are set as ”Hit = 5” and “Miss = 0”. The tentative track not satisfying the 2/2 and 2/3 criteria is deleted from the tentative track list. If a confirmed track does not receive any measurements during five consecutive scans of data, i.e., “Hit = 0” and “Miss = 5”, the delete flag for this track is set to 1 and the track is deleted from the confirmed tracks’ list [29].

8. Performance Evaluation Simulations

This section presents the performance evaluation of the proposed 3D MTT algorithm on the basis of a number of performance evaluation metrics for non-maneuvering and maneuvering multi-target scenarios in the presence of process noise, measurement noise, and clutter due to false alarms.

8.1. Performance Metrics

The performance evaluation metrics considered in the proposed research include the following.

8.1.1. Root-Mean-Squared Error in Position

Assume $[x_k,y_k,z_k]$ and $[{\widehat{x}}_k,{\widehat{y}}_k,{\widehat{z}}_k]$) represent the true and estimated positions, respectively, of a target at time instant k in 3D Cartesian coordinates. For M number of Monte Carlo runs, the RMSE in the position of the target is defined as [24,29]

$${\mathrm{RMSE}}_k^{pos}=\sqrt{\frac{1}{M}\sum _{i=1}^M{({\widehat{x}}_k^i-x_k^i)}^2+{({\widehat{y}}_k^i-y_k^i)}^2+{({\widehat{z}}_k^i-z_k^i)}^2}$$

(83)

8.1.2. Root-Mean-Squared Error in Velocity

Similarly, if $[v_{x_k},v_{y_k},v_{z_k}]$ and $[{\widehat{v}}_{x_k},{\widehat{v}}_{y_k},{\widehat{v}}_{z_k}]$ represent the true and estimated velocities, respectively, of a target, the RMSE in velocity for M number of Monte Carlo runs at time instant k can be written as

$${\mathrm{RMSE}}_k^{vel}=\sqrt{\frac{1}{M}\sum _{i=1}^M{({\widehat{v}}_{x_k}^i-v_{x_k}^i)}^2+{({\widehat{v}}_{y_k}^i-v_{y_k}^i)}^2+{({\widehat{v}}_{z_k}^i-v_{z_k}^i)}^2}$$

(84)

8.1.3. Mean Execution Time for One Data Scan

The mean execution time $T_E$ of the algorithm for one data scan was computed on a laptop computer. The specifications of the computer were 1.9 GHz processor, 4GB RAM, and Windows 10 for Matlab2018. The total execution time was composed of the time for data association $T_{DA}$, filtering for state estimation $T_F$, and track management $T_{TM}$ functions [24,29].

$$T_E=T_{DA}+T_F+T_{TM}$$

(85)

8.1.4. IMM Mean Model Probabilities

The IMM mean model probabilities for maneuvering targets reflect how efficiently the IMM algorithm can recognize the different dynamic motion models of the targets and switch between different filter models accordingly [29].

8.2. Parameter Selection and Simulated Data Generation

To evaluate the performance of the proposed MTT algorithm, the targets were considered as point sources in the far-field relative to the observing radar. A DF-FMCW interferometric radar with $f_{c_1}=6$ GHz and $f_{c_2}=24$ GHz was simulated with a bandwidth $B=500$ MHz and a sweep time $T=1$ ms. The sampling frequency was set to be $f_s=128$ kHz. The receiving antennas Rx1, Rx2, and Rx3, corresponding to $f_{c_1}$, were located at (0,0,0), ($-D_{12}$,0,0), and (0,0,$D_{13}$), where $D_{12}=D_{13}=3$ m [29]. The target trajectories were modeled using the NCV and NCT dynamic models with discrete white Gaussian process noise. The straight line motion of the targets followed the NCV model, whereas the maneuvering motion of the targets followed the NCT model. The standard deviations associated with the linear and circular segments of target motion were assumed to be ${\sigma }_{{\mathbf{v}}_1}=0.20$ and ${\sigma }_{{\mathbf{v}}_2}=0.10$. Further, the radial velocity, azimuth angular velocity, and elevation angular velocity measurement error standard deviations were set as ${\sigma }_{v_r}=0.14$ m/s, ${\sigma }_{{\omega }_{\theta }}=0.10$ rad/s, and ${\sigma }_{{\omega }_{\phi }}=0.10$ rad/s, respectively. The probability of gate detection was set as $P_G=0.99999$. The sampling interval was $\Delta T=0.02$ s [29]. The number of Monte Carlo simulation runs was $M=20$. The clutter points assumed a Poisson distribution and were uniformly distributed in the measurement region [42]. Each clutter point was composed of (i) a radial velocity measurement distributed in the range of [$v_{{\mathrm{r}}_{min}},v_{{\mathrm{r}}_{max}}$], (ii) an azimuth angular velocity measurement distributed in [${\omega }_{{\theta }_{\min }},{\omega }_{{\theta }_{\max }}$], and (iii) an elevation angular velocity measurement distributed in [${\omega }_{{\phi }_{\min }},{\omega }_{{\phi }_{\max }}$]. The average number of clutter points was set as 5 for each scan [29]. For the IMM filter with Model 1 as the NCV model and Model 2 as the NCT model, the model transition probabilities matrix is

$$\begin{array}{c}\hfill \pi =\left[\begin{array}{cc}0.99& 0.01\\ 0.01& 0.99\end{array}\right]\end{array}$$

(86)

The initial model probability for each filter model was chosen to be ${\mu }_j=0.5$, where $j=1,2$ [29].

8.3. Scenario 1: Two Non-Maneuvering Closely Spaced Targets

Scenario 1 for 3D tracking with two non-maneuvering targets and radar antennas is shown in Figure 3.

Table 1. Scenario 1.

Targets				Initial States of Targets
	$\mathit{x}\mathbf{\left(}\mathbf{m}\mathbf{\right)}$	${\mathit{v}}_{\mathit{x}}1\mathbf{(}\mathbf{m}\mathbf{/}\mathbf{s}\mathbf{)}$	$\mathit{y}\mathbf{\left(}\mathbf{m}\mathbf{\right)}$	${\mathit{v}}_{\mathit{y}}\mathbf{(}\mathbf{m}\mathbf{/}\mathbf{s}\mathbf{)}$	$\mathit{z}\mathbf{\left(}\mathbf{m}\mathbf{\right)}$	${\mathit{v}}_{\mathit{z}}\mathbf{(}\mathbf{m}\mathbf{/}\mathbf{s}\mathbf{)}$	$\mathbf{\Omega }\mathbf{(}\mathbf{rad}\mathbf{/}\mathbf{s}\mathbf{)}$
1	0	3.5	3	1	1.5	0	4.83
2	0	−2	2.6	1	1.8	0	−4.83

presents the initial states of the targets. The trajectories of the targets were modeled using the NCT dynamic model defined by Equation (23). However, the circular motion was along the xz-plane, rather than the xy-plane, as described by state transition matrix ${\mathbf{F}}_{k-1}^{\mathrm{NCT},3\mathrm{D}}$ in Equation (25).

The range–radial velocity map ($f_{c_2}=24$ GHz) represented by Figure 4a shows two targets with initial ranges of $R_{01}=3.47$ m and $R_{02}=3.11$ m and initial radial velocities of $v_{r_1}=1.31$ m/s and $v_{r_2}=1.1$ m/s, respectively.

The time-varying Doppler spectrogram ($f_{c_2}=24$ GHz), azimuth interferometric spectrogram ($f_{c_1}=6$ GHz), and elevation interferometric spectrogram ($f_{c_1}=6$ GHz) are shown in Figure 4b–d, respectively. The instantaneous Doppler frequency shift caused by the radial velocity of the target increases for the target moving away from the radar and decreases for the target moving towards the radar. Moreover, the sinusoidal patterns of the radial velocities represent the circular motions of the targets moving towards and away from the radar. The instantaneous azimuth angular frequency caused by the azimuth angular velocity is positive for the target moving along the positive x direction and negative for the target moving along the negative x direction. Similarly, the instantaneous elevation angular frequency caused by the elevation angular velocity is positive for the target moving in the positive $\phi$ direction and negative for the target moving in the negative $\phi$ direction.

Figure 5a–c represent the ideal and extracted radial, azimuth angular, and elevation angular velocities of the targets, providing complete 3D velocity measurement information. The extracted 3D velocity measurements are fed to the proposed GNN-IMM-SCKF algorithm for 3D tracking, which includes data association, state estimation, and track management functions.

The real target tracks following Equation (23) and the output of the proposed algorithm as estimated target tracks in the 3D Cartesian space are shown in Figure 6a. Moreover, Figure 6b–d represent the target tracks in the 2D xy, xz, and yz Cartesian spaces, respectively.

To access the performance of the proposed algorithm, the RMSEs in the positions and velocities of the targets are plotted in Figure 7a,b, respectively, which prove that the proposed algorithm based on 3D velocity measurements from the interferometric radar can be used for MTT in the 3D Cartesian space.

8.4. Scenario 2: Two Maneuvering Closely Spaced Targets

Scenario 2 for 3D target tracking with two closely spaced maneuvering targets and radar location is shown in Figure 8. The trajectories of the targets were modeled by the NCV (Equation (18)) and NCT (Equation (23)) dynamic models. The initial states of the targets are summarized in

Table 2. Scenario 2.

Targets				Initial States of Targets
	$\mathit{x}\mathbf{\left(}\mathbf{m}\mathbf{\right)}$	${\mathit{v}}_{\mathit{x}}\mathbf{1}\mathbf{(}\mathbf{m}\mathbf{/}\mathbf{s}\mathbf{)}$	$\mathit{y}\mathbf{\left(}\mathbf{m}\mathbf{\right)}$	${\mathit{v}}_{\mathit{y}}\mathbf{(}\mathbf{m}\mathbf{/}\mathbf{s}\mathbf{)}$	$\mathit{z}\mathbf{\left(}\mathbf{m}\mathbf{\right)}$	${\mathit{v}}_{\mathit{z}}\mathbf{(}\mathbf{m}\mathbf{/}\mathbf{s}\mathbf{)}$	$\mathbf{\Omega }\mathbf{(}\mathbf{rad}\mathbf{/}\mathbf{s}\mathbf{)}$
1	−2.25	3	3.5	0	−2	1	4.83
2	1.5	−2	3.5	0	2	−1	−4.83

.

The range–radial velocity map plotted in Figure 9a clearly shows two different targets in the radar’s FOV with initial ranges of $R_{01}=4.68$ m and $R_{02}=4.35$ m and initial radial velocities of $v_{r_1}=-1.91$ m/s and $v_{r_2}=-1.11$ m/s, respectively.

Time-varying Doppler, azimuth interferometric, and elevation interferometric spectrograms for Scenario 2 are shown in Figure 9b–d, respectively. The radial, azimuth angular, and elevation angular velocities extracted from the time-varying spectrograms are presented in Figure 10a–c, respectively, together with the ideal 3D velocity measurements. For targets moving towards the radar, the radial velocities decrease, and for targets moving away from the radar, the radial velocities increase. Similarly, the azimuth and elevation angular velocities of the targets moving in positive $\theta$ and $\phi$ are positive, respectively, and vice versa.

The real target tracks and the estimated target tracks obtained as the output of the GNN-IMM-SCKF algorithm in the 3D Cartesian space are plotted in Figure 11a. The real and estimated 2D tracks in the xy, xz, and yz Cartesian spaces are shown in Figure 11b–d, respectively.

For the performance evaluation, the RMSEs in the positions and velocities of the targets are presented in Figure 12a,b. Figure 12a,b assert the fact that the RMSEs for the GNN-IMM-SCKF algorithm are less compared to the RMSEs for the GNN-NCV-SCKF algorithm during the targets’ maneuvers. The IMM mean model probabilities for the NCV and NCT dynamic models presented in Figure 12c validate the fact that the proposed GNN-IMM-SCKF algorithm based on 3D velocity measurements obtained from the interferometric radar is applicable to 3D tracking systems for state estimation of maneuvering targets.

The execution times for different 3D target tracking scenarios are summarized in

Table 3. Execution times.

Scenario	Execution Time
	${\mathit{T}}_{\mathit{E}}\mathbf{\left(}\mathbf{ms}\mathbf{\right)}$	${\mathit{T}}_{\mathit{DA}}\mathbf{\left(}\mathbf{ms}\mathbf{\right)}$	${\mathit{T}}_{\mathit{F}}\mathbf{\left(}\mathbf{ms}\mathbf{\right)}$	${\mathit{T}}_{\mathit{TM}}\mathbf{\left(}\mathbf{ms}\mathbf{\right)}$
1 (GNN-IMM-SCKF)	17.7	10 (56%)	5.94 (34%)	1.80 (10%)
2 (GNN-NCV-SCKF)	15.1	7.25 (48%)	4.41 (29%)	3.44 (23%)
2 (GNN-IMM-SCKF)	19.2	9.18 (48%)	7.53 (39%)	2.48 (13%)

. It is evident that the data association block of the tracking algorithm consumes more time than the filtering and track management blocks, which is approximately 48–56% of the total execution time. Furthermore, the execution time for the GNN-IMM-SCKF algorithm is approximately 27% greater than the execution time for the GNN-NCV-SCKF algorithm. The total execution time $T_E$ being less than the measurement sampling interval $\Delta T$ proves the fact that the proposed MTT algorithm can be implemented for real-time applications.

9. Conclusions

MTT typically requires either a network of Doppler radar receivers at different locations or a single phased array radar. However, Doppler radar networks have high computational complexity and data throughput attributed to multiple receivers in the network. Moreover, array signal processing techniques for phased array radar are computationally expensive. To resolve the issue, this paper presented an algorithm for detection and tracking of multiple moving targets based on 3D velocity measurements obtained from a dual-orthogonal baseline interferometric radar. First, we presented the mathematical model of the 3D velocities of multiple moving point sources as targets. The radial, azimuth angular, and elevation angular velocities were extracted using reference receiving antenna, a baseline along the x-axis and a second orthogonal baseline along the z-axis, respectively. Then, we derived the nonlinear 3D velocity measurement function, which defines the relationship between the 3D velocity measurements and the state of the target. Based on the 3D velocity measurement function, we introduced the design and implementation of an MTT algorithm, which included the GNN for data association, the IMM-SCKF estimator for state estimation, and the rule-based M/N logic for track management. Then, we performed Monte Carlo simulations for different multi-target scenarios and evaluated the performance of the algorithm in terms of the RMSEs in position and velocity, the mean execution time, and the IMM mean model probability. In order to simulate a multi-target scenario close to a real environment, process noise in the dynamic motion models of the targets was modeled as the DWNA, taking into consideration small accelerations as noise. Moreover, zero-mean Gaussian measurement noise and clutter due to false alarms following a Poisson distribution were added to the received data. Consequently, the simulation results proved the fact that the proposed algorithm is robust and can be applied to practical 3D tracking systems.