MIMO Radar: An H-Infinity Approach for Robust Multitarget Tracking in Unknown Cluttered Environment

In a multiple-input multiple-output (MIMO) radar, adaptive multitarget tracking can be achieved using subspace tracking algorithms in conjunction with super-resolution subspace localization algorithms. However, the presence of unknown radar clutter deteriorates or breaks algorithms that rely on the white noise assumption. In this article, an <inline-formula><tex-math notation="LaTeX">$H^{\infty }$</tex-math></inline-formula> approach is proposed for robust tracking of each target's range, direction-of-arrival and velocity with unknown clutter. Specifically, two “manifold extenders” are first proposed by combining the slow-time and fast-time dimensions of a pulse MIMO radar's received signal. Then, in the “extended” space, an <inline-formula><tex-math notation="LaTeX">$H^{\infty }$</tex-math></inline-formula> adaptive algorithm is proposed to track an equivalent noise subspace, which exists regardless of the noise assumption. Finally, the target parameters are extracted from the adaptively tracked noise subspace. Based on computer simulation studies, the performance of the proposed <inline-formula><tex-math notation="LaTeX">$H^{\infty }$</tex-math></inline-formula> tracking approach is evaluated using challenging tracking scenarios and compared against several existing subspace tracking algorithms that have been modified to operate on the “extended” space.

Trace of matrix A. eig min (A) Minimum eigenvalue of matrix A. col k (A)  The k th column of matrix A.

R, C
Set of real and complex numbers.

I. INTRODUCTION
A monostatic multiple-input multiple-output (MIMO) radar that uses antenna arrays at both the transmitter (Tx) and receiver (Rx) will be considered in this article for multitarget tracking in the presence of unknown clutter.In general, the target environment can be classified as stationary and nonstationary.In nonstationary environments, the locations of the targets are fast-varying and may need to be tracked on a snapshot-by-snapshot basis, e.g.[1], [2].However, in stationary environments, the locations of the targets remain unchanged for the observation period.In this case, the target locations can be estimated using various localization MIMO algorithms.Popular MIMO algorithms include the Capon and amplitude and phase estimation (APES) algorithms, as well as the combined Capon and APES (CAPES), and the combined Capon and approximate maximum likelihood (CAML) algorithms (see [3]).However, the performance of these algorithms decreases significantly for closely spaced targets due to high mutual target interference.In [4], a constrained multivariable optimization problem is formulated to suppress both the mutual target interference and noise, thus achieving superior target localization performance.Note that super-resolution subspace-based localization algorithms, which may exploit the concept of the array manifold extender [5], [6] for improved parameter identifiability and estimation accuracy, can also be used.
In addition to the stationary and nonstationary target cases, there is a third radar case where the target locations remain unchanged only over a small observation interval known as the coherent processing interval (CPI) [7], [8].Target tracking in this environment can be achieved through target localization within each CPI, by using each CPI data independently.Note that this is known as tracking based on repetitive localization from CPI to CPI.However, from CPI to CPI, the "target association" rule is crucial for ensuring that each track is actually caused by the same target across time [9].Information, such as a priori estimation of the target location, can be used for the association rule [10].
Instead of repetitive localization tracking approaches, a number of adaptive algorithms have been developed in [11], [12], [13], [14], [15], and [16] for tracking the signal or noise subspaces.These adaptive algorithms are known as subspace tracking.For instance, in [11], the data projection method (DPM) algorithm is developed for solving a constrained optimization problem.This requires constantly orthonormalizing the subspace estimate using the Gram-Schmidt orthogonalization procedure to force the algorithm to converge.A more efficient very well-known algorithm is the projection approximation subspace tracking (PAST) algorithm [12].Moreover, the yet another subspace tracker (YAST) algorithm, proposed in [13], outperforms many other approaches but is numerically unstable.In order to reduce the complexity and improve the numerical stability, various extensions of DPM, PAST, and YAST algorithms have been proposed over the years [14], [15], [16].These subspace trackers can be used to facilitate the solution to the adaptive target localization problem, such as "adaptive MUSIC" [17].
However, most of these subspace trackers are based on the white noise assumption.As in radar the target echo is corrupted by noise plus clutter, this provides an effective noise which is spatially colored, i.e. nonwhite noise.Thus, the covariance matrix of this noise is not diagonal anymore.In this case, prior color noise whitening method [18] can be used to diagonalize the noise-plus-clutter covariance matrix, e.g.[5], [19].Alternatively, oblique subspace tracking can be used [20], [21].For instance, the oblique PAST (ob-PAST) algorithm proposed in [20] embeds this prewhitening into the update procedure and does not need matrix decomposition.In [21], the YAST algorithm is extended to handle oblique subspace tracking in a Riemannian framework but only applies to real signals.Inevitably, both the prewhitening procedure and the oblique subspace tracking approach need a priori knowledge of the noise-plus-clutter covariance matrix.In [22], the tracking algorithm effectively tracks the target manifold vectors (signal subspace) by solving a constrained least-square (LS) problem.It involves the recursion and regularization phases, which implies that the parameter estimation is essential in every CPI in order to regularize the estimated subspace.
In addition, solutions to multitarget tracking can also be probabilistic.Probabilistic tracking approaches recursively estimate the unknown multitarget posterior probability density based on set-derivatives and set-integrals.Then an estimator, such as the maximum a posteriori (MAP) estimator, is applied to determine the target locations [23].For instance, the optimal Bayes multitarget filter, probability hypothesis density (PHD), cardinalised PHD (CPHD), and multi-Bernoulli filters [23].Using the covariance matrix estimated from multiple snapshots in each CPI, a probabilistic tracking algorithm based on the CPHD filter is proposed in [24] for multitarget direction-of-arrival (DOA) tracking.However, this type of approaches is out of the scope of this article.
In this article, a new tracking algorithm is proposed which belongs to the family of subspace tracking algorithms.The proposed algorithm does not rely on the white noise assumption and it is based on the H ∞ theoretical framework.This H ∞ algorithm is suitable for DOA, range and velocity multitarget tracking with the radar operating in the presence of colored clutter in an environment that is stationary over a CPI but nonstationary across different CPIs.It is important to point out that the H ∞ framework initially is developed from the minimax estimation in robust control theory with the aim of minimizing (or in the suboptimal case, bounding) the maximum energy gain1 from the disturbances to the estimation errors, e.g.[25], [26], [27], [28].The main advantages of H ∞ approaches include being more robust and less sensitive to model uncertainties, parameter variations, and the lack of statistical information on the noise signals [29], [30].In other words, the H ∞ estimator guarantees that if the disturbances are bounded (in energy), then the estimation errors will be as small as possible, no matter the disturbances and their statistics (and/or distributions).In contrast, for a nonrobust algorithm, such as LS (or H 2 ) based recursive least-squares (RLS) or Kalman filter, small disturbances may lead to large estimation errors if the assumptions on the disturbances are violated [31].The H ∞ tracking approach has been proposed in the past for applications, such as passive source tracking [32], communications [33], [34], and radar [8].Under the H ∞ framework, a propagator is recursively estimated, which can be used to form a subspace that is orthogonal to the signal subspace.In [8], the H ∞ tracking approach is first proposed for radar applications but does not examine the performance in a cluttered environment.Moreover, it considers only the tracking of range and DOA by employing an array manifold extender that exploits the fast-time dimension of the received signal of a pulse MIMO radar.As an extension of [8] and the Doppler-spatial array manifold extender used in MIMO radar [5], [35], [36], this article introduces another two manifold extenders which combine the slow-time and fast-time dimensions, thus capturing both the range and Doppler information for each target.The use of manifold extenders significantly increases the system degrees of freedom (DOFs) without adding more hardware.Higher DOFs in radar implies: • handling more targets; • increased detection and resolution capabilities; • smallest estimation error.This also facilitates tracking of all the parameters of interest in a unified way as the signal subspace is now defined collectively by the extended manifold vectors of all the targets.The rest of this article is organized as follows.In Section II, the transmit and receive signal models are presented.In Section III, two manifold extenders are introduced, termed as the Doppler-STAR and virtual Doppler-STAR manifold extenders.In Section IV, an equivalent noise subspace is introduced that is orthogonal to the signal subspace formed by the Doppler-STAR and virtual Doppler-STAR manifold vectors of all the targets.The proposed H ∞ tracking approach that is based on recursively estimating the "noise-subspace equivalent" is presented in Section V. Computer simulation studies are provided in Section VI to demonstrate the performance of the proposed H ∞ tracking approach in challenging scenarios and compare it with some of the existing subspace tracking algorithms.Finally, Section VII concludes this article.

II. MIMO RADAR SIGNAL MODELING
Consider a monostatic pulse radar system whose Tx and Rx employ antenna arrays of N and N antennas, respectively.Both the Tx and Rx arrays are assumed to be fully calibrated with known geometries.The Cartesian coordinates of the elements of the two arrays are represented in units of half-wavelength by the columns of the matrices r and r as follows2 : Consider the radar operates in the presence of K moving (relative to the radar) targets that are located in the array far-field, where K is assumed3 to be known or preestimated.

A. Transmit
With reference to Fig. 1, consider all the targets remain stationary over a CPI of D PRIs.The Tx-array transmits, during each pulse repetition interval (PRI) of duration T r , N s i.e. the elements of the vector symbol a [n] are simultaneously transmitted via N transmit antennas.The transmit symbol matrix is defined as which satisfies the orthonormality property In order to further increase the range resolution, each vector symbol is multiplied with a symbol compression code of , where T c is the chip period.Based on the above discussion, the data sequence transmitted during each pulse duration can be written as ( With reference to Point A in Fig. 2, this data sequence is driven to a digital-to-analog converter (DAC) to produce the baseband transmit vector signal m(t ) ∈ C N×1 .

B. Receive
The parameters including the range, DOA, and velocity of all the targets are assumed to remain approximately unchanged within each CPI.However, these parameters are varying across different CPIs.For a particular PRI of any CPI, with reference to Point B in Fig. 2, the time-continuous received vector signal x(t ) ∈ C N×1 can be modeled as ) where j = √ −1, and for the k th target, β k is the complexvalued path gain, R k is the range, and the Doppler frequency f k is related to the target radial velocity v k as In (7), λ = c/F c is the wavelength with c and F c denoting the wave propagation speed and the carrier frequency, respectively.Moreover, S k and S k represent, respectively, the k th target's Tx and Rx manifold vectors [39], which can be expressed as functions of the target DOA specified by azimuth θ k and elevation φ k as where is a unit-norm vector directed toward (θ k , φ k ).Finally, in (6), the vector n cn (t ) ∈ C N×1 is defined as n cn (t ) x c (t ) + n(t ) (11) where the vectors x c (t ) and n(t ) represent, respectively, the clutter and additive white Gaussian noise (AWGN) of zero mean and covariance σ 2 n I N with σ 2 n being the noise power.Assuming uncorrelated clutter and noise terms, the covariance matrix of n cn (t ) is given as which is no longer diagonal.

III. MANIFOLD EXTENDERS: COMBINING THE SLOW-TIME AND FAST-TIME DIMENSIONS
The parameters to be tracked for each target are the range R k , DOA 4 θ k and the radial velocity v k .In order to achieve a joint tracking of all these parameters, the concept of the array manifold extender [6], which allows additional system and channel parameters to be incorporated into the conventional (or spatial) array manifolds, will be exploited.Note that various types of manifold extenders have been proposed in radar systems for target localization and tracking [5], [8], [35], [40].Based on the Doppler-spatial array manifold extender, another two "extenders" are proposed in this section by combining the slow-time and fast-time dimensions.
In the Doppler-spatial array manifold extender, same range bin snapshots received from different PRIs are grouped and processed together to increase the dimensionality of the observation space [35], [36].This exploits the slow-time dimension of the received signal and represents the simplest manifold extender in MIMO radar.With reference to Point C in Fig. 2, consider the signal modeled by ( 6) is passed through an analog-to-digital converter (ADC) with a sampling period T c and, for each of D PRIs, N s N c snapshots containing the targets are collected.Then, a 3-D CPI data cube can be formed, as shown in Fig. 3, where the N × D slice associated with the l th snapshot of the n th symbol is denoted by X l [n].With reference to Point C1 in Fig. 4, the ND × 1 space-time snapshot x l [n] is obtained by vectorizing X l [n] as Then, consider in total 2N c snapshots associated with the n th and (n + 1)th symbol periods, an ND × 2N c space-time data matrix X[n] can be obtained and expressed as where which contains all the N c snapshots of the current symbol a[n] as well as parts of the previous symbol a[n − 1] and the next symbol a[n + 1].In (16), the vector c ∈ R 2N c ×1 is defined as and the target range bin l k ∈ [0, N c − 1] is calculated as where R = cT c /2 is the range resolution.It should be highlighted that at the Rx we consider only N s N c snapshots containing the targets for every PRI, which are associated with the transmitted data sequence in (5).This implies l k cannot be directly translated to the target range R k due to the modulo operation.Moreover, the matrix J is a 2N c × 2N c lower shift matrix defined as which has the property that its l k th power J l k or (J T ) l k , applied on the vector c, downshifts or upshifts c by l k elements.Finally, N cn [n] is the so-obtained ND × 2N c noise-plus-clutter data matrix.
In ( 14), the Doppler shifts within each pulse are neglected, assuming the product between the pulse duration N s N c T c and the Doppler frequency f k ∀k is small compared to unity.Note that the Doppler-spatial manifold vector of the k th target is given as which does not include the target range bin l k .Conventionally, a detection algorithm is required to run through all the range bins to detect the targets and extract the corresponding snapshots for localization and tracking.This implies the target range tracking is independent from DOA and velocity tracking and failed range estimation can lead to false tracks.
A. Proposed Manifold Extender 1: Doppler-STAR In order to further incorporate the target range bin l k into the Doppler-spatial manifold vector in (20), the data matrix X[n] modeled by ( 14) is vectorized to a 2N c ND × 1 column vector y[n].We will call this vector at Point C2 in Fig. 4 a "spatiotemporal snapshot", which can be modeled as where h k represents the 2N c ND × 1 dimensional Doppler-STAR manifold vector of the k th target associated with the current symbol a[n] and is defined as Note that h prev,k and h next,k in (21) are Doppler-STAR manifold vectors associated with the previous and the next symbols, respectively.These two are related to h k as Finally, the spatiotemporal noise-plus-clutter snapshot is obtained as ).Based on ( 23) and (24), Eqn ( 21) can be written compactly as where The transmit array can also be exploited at the Rx.This can be done by using the concept of the "virtual array" to further increase the DOFs [41].In particular, by taking the outer product between y[n] in (25) and the transmit symbol a[n], i.e. y[n]a H [n] ∀n ∈ [1, N s ], the 2N c NND × 1 "virtual spatiotemporal snapshot" [5] at Point C3 in Fig. 4 is obtained.That is, where z v is the further extended noise-plus-clutter snapshot, and the matrix T v is the following known transformation matrix: Furthermore, in (26), the k th column of H v is defined as and represents the 2N c NND × 1 dimensional extended manifold vector of the k th target, termed the "virtual Doppler-STAR."

IV. NOISE SUBSPACE IN THE "EXTENDED" SPACE
In this article, the snapshots at the output of two manifold extenders proposed in Sections III-A and III-B will be used (Point D in Fig. 2).To simplify the notation and to use the H ∞ framework proposed in this article for both manifold extenders 1 and 2, let us redefine the following parameters for the th CPI as follows: Moreover, the path gain vector β[ ] varying from CPI to CPI is assumed to contain zero-mean independent random variables and has covariance5 diag(σ 2 β ).Note that the two manifold extenders extend the N-dimensional observation space to an NN ext -dimensional complex space, where The target subspace spanned by the columns of H[ ] ∈ C NN ext ×K contains information along both the slow-time and fast-time dimensions, enabling various subspace-based methods to estimate or track all the parameters in a unified way.
Consider the extended snapshot y[ ] received in the th CPI of length NN ext , it can be y[n] in (25) for any fixed n for every CPI (e.g.y[N s ]) or y v in (26).This long vector can be partitioned into two subvectors y a [ ] and y b [ ] (the upper and lower subvectors) as follows: This implies that, the matrix H[ ] containing all the target extended manifold vectors and the noise-plus-clutter snapshot z[ ] (i.e.z[n] for a fixed value of n or z v ) are also partitioned accordingly as follows: H[ ], and z[ ], respectively.It is assumed that the partition is performed such that the K × K matrix H a [ ] is nonsingular, i.e. full rank.Note that this also requires any two targets not to share simultaneously the same DOA, Doppler frequency and range bin.In practice, we may choose more than K rows for y a [ ] to guarantee that the rank of H a [ ] is at least equal to K [42].
In any case, there exists a propagator (linear operator) From (35), it can be seen that where Equation (36) implies that the columns of U[ ] belong to the complement subspace spanned by the columns of H[ ], therefore it can be regarded as a "noise-subspace equivalent".Most importantly, this orthogonality holds regardless of the second-order statistics of z[ ].
The propagator L[ ] can be estimated from the data, if its covariance matrix R[ ] ∈ C NN ext ×NN ext is known or can be estimated.For manifold extender 1, the covariance matrix of y[n] in (25) can be calculated using all the available N s snapshots and written as where R t [ ] and R z [ ] represent the covariance matrices of the target return and noise-plus-clutter, respectively.Note that in (38), the covariance matrix R mm [ ] is It is clear from R t [ ] [see (38)] that the targets lie in a subspace spanned by the columns of H[ ] having dimensionality 3 K.Consider partitioning R[ ] in (38) as follows: where From (35), it can be seen that the following relation holds: Therefore, the propagator L[ ] can be estimated as

V. MULTITARGET H-INFINITY TRACKING
In this section, an H ∞ tracking approach that is applicable to the extended snapshot6 y[ ] obtained using both manifold extenders is proposed.Based on the discussion in Section IV, each target's DOA, Doppler frequency, and range can be traced by tracking the propagator L[ ] adaptively across different CPIs, then in each CPI forming the "noise-subspace equivalent" U[ ] and extracting the parameters of interest.It is assumed that any two targets do not have the same DOA, Doppler frequency, and range bin at the same time within any CPI.This implies that the matrices H[ ] ∀ have full column rank, i.e. the Doppler-STAR or the virtual Doppler-STAR manifold vectors of all the targets are linearly independent during the whole tracking period.
Noting that the target range bin l k [ ] only takes discrete integers between 0 and N c − 1, let us define the k th target's state vector in the th CPI p k [ ], containing the other two continuous parameters.That is, Then, the following state-space model can be formulated: where • the vector ε[ ] models the approximation error, i.e.
• the propagator L[ ] is the state variable in the th CPI; • the matrix L[ ] represents the unknown variation in the state variable from the ( − 1) th to the th CPI; • the matrix B is a user-defined bound using prior knowledge of the maximum rate of change in the state variable; • the vector p k [ ] represents the unknown variation of the k th target in the state vector p k [ ] from the ( − 1) th to the th CPI.
Under the H ∞ framework [29], the state variable L[ ] can be updated from L[ − 1] based on the following equation: (48) and with G[ ] following the Riccati recursion: + BB H (50) where γ represents the robustness of the H ∞ estimator with smaller values indicating greater robustness.
Once L[ ] is obtained, the target DOAs, Doppler frequencies, and range bins can be estimated by solving the following unconstrained optimization problem: where In (52), P U[ ] is the projection operator onto the subspace spanned by the columns of U[ ] given by (37).Furthermore, the parameters T( ) and h(θ, f) for the two manifold extenders are defined 7 as follows: for Manifold Extender 2. ( 54) However, (52) can be seen as the generalized eigenvalue decomposition (GEVD) formulation; thus, its minimization can be reduced to minimizing the following 1-D cost function: Once the target range bin is found, by assuming the rate of change of the target state vector is small over two adjacent CPIs, the Newton iteration formula can be employed to compute the variation in the state vector p k [ ] as follows: for the previous symbol; where ∇ p h(p) represents the gradient of h(p) with respect to the state vector.Moreover, in order to alleviate the computational burden in the tracking process, it should be noted that the calculation of the orthogonal projection operator P U[ ] , i.e.
can be simplified using the matrix inversion lemma8 as follows: This can be further simplified by performing the QR decomposition using the Householder transformation as where Q[ ] and R[ ] are the K × K unitary and upper triangular matrices, respectively.Therefore, (59) can be written as (61) Based on the previous discussion, the proposed H ∞ approach for robust multitarget tracking is summarized in Table I.It should be noted that no prior knowledge about the clutter statistics is required.

VI. COMPUTER SIMULATION STUDIES
Computer simulation studies are provided in this section to evaluate the performance of the proposed multitarget H ∞ tracking approach.Consider a monostatic MIMO radar that employs a Tx-array of N = 4 antennas and an Rx-array of N = 6 antennas.The Tx-array has a uniform linear array (ULA) geometry with interantenna spacing 3λ/2 while the Rx-array has a uniform circular array (UCA) geometry of λ/2 interantenna spacing.Fig. 5 shows both the Tx and Rx array geometries with the array centroids as their reference points (i.e.origin of their Cartesian coordinates).Note that the virtual array geometry, which is related to the manifold extender 2 is also shown in Fig. 6 and is given as the spatial convolution between the Tx and Rx arrays, i.e.The symbol compression code, that is the first N c chips of the vector c ∈ R 2N c ×1 used in the following simulations, is an m-sequence described by the primitive irreducible polynomial D 4 + D + 1 in Galois field 2 (GF-2).For the th CPI, the signal-to-noise-plus-clutter ratio (SNCR) of the k th target and the clutter-to-noise ratio (CNR) are defined, respectively, as follows: where P c [ ] represents the clutter power in the th CPI.
In the simulations, the clutter is assumed dominant with CNR fixed at 20 dB for all the CPIs.The clutter snapshots, assumed following the complex Weibull distribution, are generated according to a prescribed covariance matrix [43].
The path coefficients of all targets β k [ ] ∀k are modeled using the Swerling III model with equal mean and phase being uniformly distributed in [0, 2π ].The scalar parameters for H ∞ tracking are chosen to be μ = 0.1, ζ = 0.01, and γ = 5; however, these parameters can be further optimized given a particular tracking environment.Moreover, when partitioning the data, as in (33), the first K = 2N c rows are selected for y a [ ] and the rest for y b [ ] as the range bins containing the targets are unknown.Note that other adopted simulation parameters are specified in Table II.The proposed H ∞ tracking approach requires any two targets not to have the same range bin, DOA, and velocity at the same time for any CPI.In practice, targets may have their range, DOA, and velocity trajectories crossing each other at different CPIs.In this subsection, such crossing situations are considered.The MIMO radar is assumed to be operating in the presence of K = 2 targets moving toward the radar, with the tracking process being carried out over a time interval of 200 CPIs.Under 20 dB SNCR, Fig. 7 presents the range tracking of the two targets, while Figs. 8 and 9 present the DOA and velocity tracking results, respectively.It is clear that two targets' range, DOA, and velocity trajectories have been successfully tracked by using the extended snapshots obtained at the output of both the Doppler-STAR and virtual Doppler-STAR manifold extenders.

B. Appearance and Disappearance of Targets
In radar, the complex path gain of a slow-fluctuating target can be considered to be a random variable that follows the Swerling III model.This implies β k [ ] may have a very small magnitude due to target's RCS fluctuations.This makes the target to "disappear" in some CPIs and then to "appear" again.Figs.10-12 illustrate the capability of the proposed H ∞ tracking approach for handling the appearance/disappearance of K = 2 targets moving away from the radar.It can be observed that the first target disappears during the CPI interval 55-70 while the second target disappears from CPI 90 up to CPI 125.

C. Effects of Radar Clutter
In order to further examine the performance of the proposed multitarget H ∞ tracking algorithm, the Monte-Carlo simulation studies have been carried out based on the tracking environment of Figs.7-9.The number of Monte-Carlo trials is set at 100.The performance of the proposed H ∞ tracking is compared against   1) two popular subspace algorithms: the PAST [12], obPAST [20], and 2) two LS algorithms [22]: unconstrained and constrained LS.
In order to track all the parameters of interest for each target, these tracking algorithms are applied at Point D in   Fig. 2 using the extended snapshots.The PAST and obPAST have been modified/enhanced so that the signal subspace is defined collectively by the extended manifold vectors of all the targets.It should be emphasized that the original algorithms proposed in [12], [20], and [22] only support tracking some of the parameters, e.g.DOA.It should be  noted that the PAST and obPAST algorithms do not guarantee the orthonormality of the estimated subspaces, therefore a reorthonormalization procedure is performed after each update [14].Note also that, for the unconstrained LS algorithm, the regularization phase required in the constrained LS algorithm is omitted.The results presented in this subsection are based on manifold extender 1 with snapshots modeled by (25).However, similar results can be obtained by using manifold extender 2 with snapshots modeled by (26).The forgetting factor for these algorithms is set to 0.95.As the proposed H ∞ tracking algorithm initially estimates the target range bin, the probability of detection is examined with the results shown in Fig. 13.In the same figure the results of PAST, obPAST, and LS are also shown.It is important to point out that the obPAST fails due to the assumption that the noise-plus-clutter covariance matrix should be perfectly known, while for the proposed H ∞ algorithm this is unknown.Thus, in Fig. 13 (but also in Figs. 14 and 15), the obPAST has been implemented with full knowledge of the noise-plus-clutter covariance matrix.It can be clearly seen that the probability of detection of the PAST algorithm is poor and does not improve with the increase of the SNCR, which is due to the fact that the white noise assumption is violated.For the unconstrained LS algorithm, a low probability of detection implies that some other detection techniques are necessary to find the correct target range bins.Note that this also applies to the constrained LS algorithm as it differs from the unconstrained version only in that it has a post regularization phase.
Next, by assuming a successful detection of the two targets, the DOA and velocity tracking performances evaluated using the root-mean-square error (RMSE) criterion are plotted in Figs. 14 and 15 as a function of the SNCR.The RMSEs are defined as where θ k [ ] and v k [ ] represent the estimated DOA and velocity of the k th target in the th CPI, respectively, and the expectation is taken over 100 Monte-Carlo trials.Considering all the available 200 CPIs, from Figs. 14 and 15, it can be observed that the PAST algorithm has failed.At the same time, the proposed H ∞ tracking offers similar DOA tracking performance and more accurate velocity tracking compared to the obPAST and constrained LS algorithms.However, it should be emphasized that the obPAST algorithm relies on full knowledge of the noise-plus-clutter covariance matrix.
As for the constrained LS algorithm, it requires constantly regularizing the signal subspace at each CPI using the estimated parameters.This implies a higher computational complexity as the parameter estimation is essential during the whole tracking period, otherwise significant deviation from the true subspace would occur and this is clearly indicated by the unconstrained LS curves.

D. Tracking Robustness
The robustness of the proposed H ∞ tracking approach is expressed by the scalar parameter γ with smaller values indicating greater robustness.In this subsection, the H ∞ tracking robustness with different values of γ is investigated using snapshots collected at the outputs of both manifold extenders.Based on the environment used in Section VI-A, consider the tracking of Target 1 two ULAs of N = N = 3 antennas for both the Tx and Rx.The interantenna spacing for the Rx ULA is λ/2 and the Tx ULA is 3λ/2.In this case, the resultant virtual array given by (62) is a ULA with NN = 9 half-wavelength spaced elements.The results shown in Figs.16 and 17 imply that, for both manifold extenders at low SNCR regimes, smaller values of γ can reduce the sensitivity to random noise variations, thus leading to potentially more accurate tracking results.However, at high SNCR regimes, a larger value of γ is preferred such that it is more sensitive to the target DOA and velocity variations.Moreover, better DOA tracking accuracy can be achieved by manifold extender 2 due to increased spatial resolution of the virtual array; but the velocity tracking accuracy becomes worse when the SNCR is low.

VII. CONCLUSION
The problem of multitarget range, DOA, and velocity tracking with unknown colored radar clutter is studied in this article.In particular, two array manifold extenders are proposed that exploit both the slow-time and fast-time dimensions of a pulse MIMO radar received signal.Based on the extended snapshots, an H ∞ subspace tracking approach is proposed for handling the multitarget tracking problem in unknown cluttered environment.Unlike most of the existing subspace tracking algorithms, the proposed H ∞ tracking approach offers the following advantages: • Does not rely on the white noise assumption.
• Does not require prior knowledge of the noise-plusclutter covariance matrix.• Does not require any forms of reorthonormalization and/or regularization.
Computer simulation results demonstrate the performance of the proposed H ∞ framework.

Fig. 1 .
Fig. 1.Illustration of the i th Tx waveform.

Fig. 2 .
Fig. 2. Baseband representation of the transmitter, the receiver and the MIMO radar channel consisting of K targets.

Fig. 4 .
Fig. 4. Illustration of MIMO radar manifold extenders, where the snapshots at Point C2 or C3 can be taken as the output, i.e.Point D in Fig. 2.
) where y a [ ], H a [ ], and z a [ ] contain K rows while y b [ ], H b [ ], and z b [ ] contain the rest NN ext − K rows of y[ ],

TABLE I Proposed
H ∞ Approach for Robust Multitarget Tracking