Joint Detection and Localization in Distributed MIMO Radars Employing Waveforms With Imperfect Auto- and Cross-Correlation

This article deals with the problem of joint detection and localization of multiple targets in a distributed multiple-input multiple-output radar system where the emitted waveforms present imperfect auto- and cross-correlation functions. In this context, the sidelobe masking effects can considerably degrade the detection and localization performance of correlation-based detectors. To address this issue, we first derive a convenient discrete-time signal model, wherein the echoes generated by a target towards each receive antenna are regarded as a subspace signal. Then, we formulate the joint detection and localization problem as a composite multiple hypothesis test and derive the optimal solution according to a generalized information criterion, whose complexity scales exponentially with the number of prospective targets. To reduce the computational burden, we derive two iterative approximate solutions which detect and localize one target at a time from the noisy data, upon mitigating the sidelobe masking caused by the previously-detected ones. Finally, we provide an extensive numerical analysis to demonstrate the effectiveness of the proposed solutions, even in the challenging case where an unknown number of strong and weak targets is present in the covered area.


I. INTRODUCTION
M ULTIPLE-INPUT multiple-output (MIMO) radars transmitting linearly-independent waveforms have recently attracted a large interest for both civil and military applications [1], [2], [3], [4].In a co-located MIMO radar, the antennas are closely-spaced at both the transmit and receive sides and see the same target aspect angle [5], [6], [7], [8]; accordingly, waveform diversity and/or frequency diverse arrays can be exploited for adaptive beam-pattern shaping and signal processing.A distributed MIMO radar, instead, has widely-spaced antennas, which can see different aspect angles of the targets: the resulting spatial and geometric diversity can be exploited to reduce the scintillation of the radar cross section (RCS) [9], [10], [11], and a multidimensional measurement space can be obtained by separating all transmit-receive paths [3].Distributed MIMO radars have been effectively used to accomplish target detection [12], [13], [14], localization [15], [16], [17], and tracking [18]; in particular, the system geometry [19], antenna placement [20], the power allocation [21], and the waveform/code design [22] can be in principle optimized to run the desired task.
There are two main approaches to perform target detection and/or localization in the distributed MIMO radars [23].In the first approach [15], [16], [17], each receiver first obtains some position-related parameters of the detected targets, such as for example the angles and times of arrival of the corresponding echoes; then, a fusion center collects these estimated parameters and identifies the targets via triangulation.This two-step process usually requires a low transmission bandwidth and computational burden, but its performance may be degraded by the presence of missed detections and/or estimation errors at some receivers.In the second approach, the raw data measurements of all receivers are directly sent to the fusion center and jointly processed [24], [25], [26], [27], [28], [29], [30], [31]: since no preliminary elaboration is locally performed by each receiver, this solution can provide better detection and/or localization performance, especially when the signal-to-noise ratio (SNR) is low, at the price of a larger signaling overhead and computational cost.In this study, we consider this second approach, and the next paragraph analyzes the related works in more detail.
The studies in [24], [25], [26] rely on maximum likelihood (ML) estimates of the target parameters to develop localization algorithms, while [27], [28] employ grid-search methods to obtain an estimate of the target position; however, these works assume the existence of only one target and omit the detection phase.On the other hand, surveillance radars are often faced with the problem of detecting an unknown number of prospective targets at unknown locations.The joint detection and localization (JDL) of multiple targets poses two major challenges.One challenge is that the reflections of multiple targets may interfere with each other, resulting in the creation of shadow and ghost targets; the other one is that expanding the search space to account for the presence of multiple targets much increases the computational burden.In this context, [29] has employed data and geometric matching algorithms to address the first issue.Also, [31] has employed the energy modeling of the multiple transmit-receive paths and the compressive sensing theory to address the second issue.Finally, [30] has attempted to consider both challenges simultaneously; however, the proposed correlation-based solutions rely on the assumption that the radiated waveforms have thumb-tack auto-correlation functions and all-zero cross-correlation functions, whereby the echoes produced by different waveforms or possessing different delays can be perfectly separated by a matched filter (MF).In practice, no waveform of duration T can have a thumb-tack auto-correlation function: indeed, the auto-correlation is zero only for lags exceeding T , while, for lags smaller than T , a convenient waveform design can mitigate to some extent but not eliminate the sidelobes.Also, N distinct waveforms have all-zero cross-correlation functions for all lags of interest only if they are sufficiently well separated in time or if they occupy non-overlapping frequency intervals; however, this requires an N -fold increment of the scan time or bandwidth for a given range resolution.When the transmit antennas employ waveforms that have a long duration to reduce the peak transmit power and that overlap both in time and frequency to reduce the scan time and the system bandwidth, then the resulting sidelobe masking effects can much degrade the performance of a correlation-based detector/estimator [32], [33], [34].
Recently, [35], [36] have attempted to model the cross-terms at the output of the MF and have studied their impact on the detection performance in a single-target scenario; instead, a delay compensation approach has been proposed in [37] to improve the detection capability.However, to the best of our knowledge, no previous work has still considered and solved the problem of JDL in a distributed MIMO radar when the employed waveforms present imperfect cross-and auto-correlation functions, even though somehow similar problems have been recently considered in different application areas.For example, [38] has considered a co-located MIMO radar and has proposed a sparse recovery algorithm to enhance the detection of weak targets overwhelmed by the sidelobes of strong ones, while an iterative interference cancellation algorithm has been developed in [39] to recover weak targets in a co-located MIMO-OFDM dual-function radar-communication (DFRC) system.Interference cancellation algorithms have been also developed in the context of passive radar systems [40], [41]: indeed, here the direct signal (from the illuminator of opportunity to the passive radar receiver) and the indirect signal (from the illuminator of opportunity to target to the passive radar receiver) may overlap both in time and frequency, whereby the sidelobes of the auto-correlation function may impair the target detection and localization performance when a MF detector is employed.Finally, [42], [43], [44], [45] have considered a DFRC system where a mmWave communication signal is employed as a source of opportunity and, upon resorting to a generalized information criterion (GIC) [46], have derived practical subspace-based detectors.
This article considers the multi-target JDL problem in distributed MIMO radars where the transmit antennas employ waveforms with imperfect auto-and cross-correlation functions.Starting from the preliminary results in [47], we derive two novel subspace procedures that are inherently robust to the range sidelobe masking problem.More specifically, the main contribution can be summarized as follows.
r We present a general signal model, where the superposition of the echoes produced by each target in response to the emitted radar waveforms is regarded as a subspace signal, whose structure is uniquely specified by the target position, once the radar geometry is given.
r Upon formulating the multi-target JDL problem as a com- posite multiple hypothesis test directly applied on the raw data, we first derive a GIC-based detector.While providing some interesting insights, this solution is however not practical, as its complexity scales exponentially in the (unknown) number of prospective targets.
r We then derive two approximate solutions to the afore- mentioned testing problem, which are referred to as the matched subspace detector with iterative estimation of the interference subspace (MSD-IS) and the matched subspace detector with iterative estimation of the interference covariance matrix (MSD-ICM).Following the same logic adopted in [45], these solutions extract one target at a time from the noisy data, upon mitigation the cross-and auto-terms caused by the previously-detected targets by resorting to a zero-forcing or an interference-plus-noise whitening filter, respectively.r We provide a numerical study to show the merits of the proposed detectors, also in comparison with the previous solutions in [30] and some genie-aided benchmarks.The results confirm that our procedures effectively overcome the sidelobe masking produced by stronger targets on weaker ones, even in the presence of small localization errors, and are robust with respect to the increase in the number of interfering targets.The closest related works are [30], [45].While adopting a similar design approach, [45] considers a different problem, as the radar receiver is collocated with a communication transmitter in a monostatic DFRC system; our new detectors are instead able to handle a multi-static configuration where targets may even belong to dependent linear subspaces.On the other hand, while addressing the same problem tackled here, [30] makes several simplifying assumptions at the design stage that, as shown later, severely degrade the performance of the resulting detectors in the considered operating scenario.
The remainder of this work is organized as follows. 1  Section II presents the system description.Section III formulates the problem of joint target detection and localization and derives a decision rule based on the GIC.Section IV illustrates the proposed solutions.Section V contains the numerical analysis.Finally, Section VI gives the conclusions.

II. SYSTEM DESCRIPTION
Consider a distributed MIMO radar with N transmit antennas located at {x tx n } N n=1 and P receive antennas located at {x rx p } P p=1 , with 2 x tx n , x rx p ∈ R 2 ; all antennas share a common time-reference and operate on the same spectrum.We denote by TX n and RX p the n-th transmit and the p-th receive antenna, respectively, for n = 1, . . ., N and p = 1, . . ., P , and by sn (t) the baseband waveform emitted by the n-th transmit antenna.Hereafter, the positions of all radar antennas and the radiated waveforms are assumed known.

A. Radar Waveforms
At the design stage, we only assume that the waveforms {s n (t)} N n=1 are time-limited to [0, T ] and essentially frequency limited to [−W/2, W/2] [48] and that they possess the following properties: 1, where f max ≥ 0 is the maximum magnitude of the Doppler shift in any target echo; P2 the time-shifted replicas {s n (t − δ n )} N n=1 are linearly independent for any δ 1 , . . ., δ N .
The first property implies that the Doppler shift in the received echoes can be neglected; this is the typical case within each pulse repetition interval in pulse radars.The second property is required to maintain the delayed waveforms separable and implies W T ≥ N ; P2 much relaxes the requirement that the waveforms have all-zero cross-correlation functions, which would imply that {s n (t − δ n )} N n=1 be orthogonal for any δ 1 , . . ., δ N .Notice that no specific assumption is made here on the auto-correlation functions.
Several studies have investigated the problem of designing time-and/or frequency-coded waveforms possessing P2, which are commonly employed in radar and communication applications: see [49], [50], [51], [52], [53], [54] and references therein.For example, if time-coded waveforms are used, up to a scaling factor accounting for the transmitted energy, we have 1 Column vectors and matrices are denoted by lowercase and uppercase boldface letters, respectively.( • ) * , ( • ) , and ( • ) H denotes conjugate, transpose, and conjugate-transpose, respectively.I M is the M × M identity matrix.a is the Euclidean norm of the vector a, while cat{a 1 , . . ., a N } is the MN-dimensional vector obtained by concatenating the M -dimensional vectors {a n } N n=1 .diag{a 1 , . . ., a N } is the N × N diagonal matrix with the elements {a n } N n=1 on the main diagonal.rank{A} and A + are the rank and the pseudoinverse of the matrix A. det{A} is the determinant of the square matrix A. Finally, i and are the imaginary unit and the convolution operator, respectively, while E{•} denotes the statistical expectation. 2For simplicity, we consider a two-dimensional model.
where ψ(t) is the chip waveform with support in [0, T ψ ], T c is the chip interval, b n = (b n,1 , . . ., b n,L ) is the code employed by the n-th antenna, and T = (L − 1)T c + T ψ .The sequences {b n } N n=1 are usually chosen to have an impulse-like aperiodic auto-correlation function and a small aperiodic cross-correlation function.A common choice is to divide T into L equal chip intervals with constant amplitude and different phases to improve the efficiency of the power amplifiers: in this case, T c = T ψ = T /L and ψ(t) = Π(t/T c ), where Π(t) = 1 if t ∈ [0, 1] and zero otherwise.

B. Continuous-Time Received Signal
In the presence of K ≥ 0 point-like targets, 3 the baseband signal at the p-th receive antenna is [24] for p = 1, . . ., P , where the noise term wp (t) is a complex circularly-symmetric Gaussian process [55], independent across the receive antennas.The first term on the right-hand side is present only if K ≥ 1; in this case, a p,n,k sn (t − τ p,n (x k )) is the echo produced by the k-th target in response to the probing sn (t) emitted by the n-th transmit antenna, for k = 1, . . ., K and n = 1, . . ., N, where are the amplitude and the delay of the echo, respectively; also, τ min and τ max are the minimum and maximum delays in the inspected region, while where c is the speed of light.The number of targets and, for K ≥ 1, their positions and amplitudes are unknown.The signal rp (t) is sent to a linear time-invariant filter, whose impulse response φ(t) has bandwidth W and support in [0, T φ ].The purpose of this filter is to remove the out-of-bandwidth noise.Let s n (t) = sn (t) φ(t) be the filtered waveform corresponding to the n-th transmit antenna, which has bandwidth W and support in [0, T + T φ ], for n = 1, . . ., N; we assume that the filter φ(t) is chosen so that {s n (t − δ n )} N n=1 remain linearly independent for any δ 1 , . . ., δ N .Also, let w p (t) = wp (t) φ(t) be the filtered noise at the p-th receive antenna, for p = 1, . . ., P .Then, r p (t) = rp (t) φ(t) can be expanded as Hereafter we make the common assumption that two replicas of s n (t) with delay δ 1 and δ 2 are resolvable if |δ 1 − δ 2 | > 1/W .Also, we provide the following definitions.
Definition 1: K ≥ 2 targets located at {x k } K k=1 are separable if, for any target pair (k, q) with k = q, there exists at least one receive/transmit antenna pair (p, n) such that s n (t − τ p,n (x k )) and s n (t − τ p,n (x q )) are resolvable, i.e., min are resolvable for any k = q and (p, n), i.e., min Otherwise, they are said partially-separable. Fig. 1 reports two examples of separable target configurations.In the top plot, the two targets are partially-separable, as their echoes along the path originated from TX 1 and ending in RX 1 are resolvable, while those along the path originated from TX 1 and ending in RX 2 are not resolvable; in the bottom plot, instead, the two targets are fully-separable as their echoes are resolvable for each receive/transmit antenna pair.

C. Discrete-Time Received Signal
In the following, we process the signal received in the interval I = [τ min , T + T φ + τ max ); as shown in Fig. 2, we assume that T + T φ + τ max ≤ T pri , where T pri is the pulse repetition interval.r p (t) is sampled in I at rate 1/T s ≥ W , thus obtaining M = (T + T φ + τ max − τ min )/T s samples; notice that, under the considered assumptions, we have M > N.After collecting these observations into r p ∈ C M , we have for p = 1, . . ., P .In (7), and represents the signature of the echo generated by the k-th target towards the p-th receive antenna when it is illuminated by the n-th transmit antenna.Also, the full-rank mode matrix contains the signatures of the N echos of the k-th target at the p-th receive antenna, 4 while a p,k = cat{a p,1,k , . . ., a p,N,k } ∈ C N is the corresponding gain vector.Finally, w p is a complex circularly-symmetric Gaussian vector with covariance C p .
For any point x k , the mode matrix S p (x k ) can be readily obtained, as it only depends upon the waveforms {s n (t)} M n=1 and the positions of all transmit antennas {x tx n } N n=1 and of the p-th receive antenna x rx p , which are known.Notice also that the covariance matrices C 1 , . . ., C P are full-rank, as the disturbance always includes the white thermal noise of the front-end of the receiver.These matrices are assumed known in the subsequent developments; in practice an estimate could be obtained by resorting to secondary data and/or by exploiting prior knowledge of the environment.

III. JOINT TARGET DETECTION AND LOCALIZATION
In this section, we formulate the multi-target JDL problem in a distributed MIMO radar employing waveforms with imperfect auto-and cross-correlation functions; then, we derive and discuss a GIC-based solution.

A. Problem Formulation
We aim to jointly detect and localize an unknown number of targets in the region under inspection, given only the knowledge of the radar geometry, the transmitted waveforms, and the available measurements.From (7), it is verified by inspection that r p is the noisy superposition of an unknown number of subspace signals originated from as many targets; in particular, the k-th subspace signal belongs to the column span of S p (x k ), which only depends upon the position x k of the corresponding target.We are faced here with a composite testing problem, namely, where H 0 is the null hypothesis, H K is the hypothesis that K targets with unknown gain vectors {a p,1 , . . ., a p,K } P p=1 are present at the unknown positions x 1 , . . ., x K , for K = 1, . . ., K max , and K max ≥ 1 is an upper bound to the number of prospective targets. 5As usual, we start by assuming at the design stage that the targets are located on a finite uniform grid G with inter-element spacing Δ g ≤ c/W .Furthermore, we only search for targets which are separable.
The negative log-likelihood function at the p-th receive antenna is readily obtained as shown next.Under H 0 , we have Under H K , for K = 1, . . ., K max , we have where is the augmented vector specifying the positions of the targets, is the set specifying all possible positions of the targets, is the augmented mode matrix containing the signatures of the targets, and, finally, is the augmented gain vector. 5K max is a design parameter which depends upon the considered application.Based on its value, the other system parameters need to be properly chosen.To be more specific, notice that the received vector in (7) belongs to an M -dimensional space, while the echoes of the k-th target lays into an N -dimensional subspace.When K = K max , a necessary condition for correctly localizing all targets is that K max ≤ (M − N + 1).

B. GIC-Based Solution
We solve here the JDL problem via a GIC-based approach.The likelihood function under each hypothesis is penalized to account for the model order, and an estimate of the number of targets, say K, is computed as [46], [56] where In the previous equation, rank{C Σ p (x 1:K )} specifies the model order under H K , while αp,1:K (x 1:K ) is the ML estimate of the augmented gain vector at the p-th receive antenna when the prospective targets are located in x 1:K ; also, η is a penalty factor which, for example, can be set to have a given probability of false alarm P fa = Pr(reject H 0 under H 0 ).Notice that αp,1:K (x 1:K ) is available in closed-form, i.e., αp,1:K (x 1:K ) = arg min Hence, upon plugging ( 17) into (15) and neglecting irrelevant constant terms, the GIC-based decision rule can be recast as where E(0) = 0, for K ≥ 1, and is the orthogonal projector onto the column span of C −1/2 p Σ p (x 1:K ).The following remarks are now in order.Σ p (x 1:K ), while η rank{Π p (x 1:K )} is a penalty term that accounts for the dimension of this subspace.Accordingly, the scoring metric in (19) integrates the penalized Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
energies computed by each receive antenna in the candidate target subspace indexed by x 1:K and compares the maximum over G 1:K with a threshold; the argument of the maximum, say x1:K = cat{x 1 , . . ., xK }, provides an estimate of the target positions under H K .Finally, the rule in (18) compares the scoring metrics under all hypotheses and makes a decision on the number of targets.
Remark 2: When two or more inspected locations present a delay offset less than 1/W between the same transmitter and the p-th receiver, then the corresponding echoes may not be resolvable, and, consequently, C −1/2 p Σ p (x 1:K ) may be ill-conditioned.In this case, we can compute the orthogonal projector Π p (x 1:K ) in ( 20) from a truncated singular value decomposition of C −1/2 p Σ p (x 1:K ), wherein only the significant eigenvalues are retained, so as to avoid an unnecessary increment of the penalty term in (19).
Remark 3: For a large value of K max , the implementation of the GIC-based decision rule is not affordable in practice, mainly because the computation of E(K) in (19)

IV. PROPOSED DETECTORS
Leveraging the design methodology in [45], we derive next two low-complexity approximations of the GIC-based solution, which avoid the joint search of multiple targets: they are iterative subspace-based algorithms, referred to as the MSD-IS and the MSD-ICM, which extract one target at a time from the measurements.We also illustrate how these solutions can effectively handle small-scale localization errors.

A. Solution Based on the MSD-IS
We propose here to solve a sequence of composite binary hypothesis tests.In each problem, we aim to detect a subspace signal generated by a prospective target in the presence of the subspace interference caused by the previously-detected targets and of the noise.
To proceed, we denote by x(k) the estimated position of the target detected in the k-th iteration and by G (k) the search set in the k-th iteration, with G (1) = G and for k ≥ 2. The k-th testing problem is formulated as p,N } ∈ C N are the mode matrix of a prospective target located in x (k) ∈ G (k) and the corresponding gain vector, while S p (x (i) ) ∈ C M ×N and a (i) p ∈ C N are the mode matrix of the target detected in the i-th iteration and the corresponding gain vector, for i = 1, . . ., k − 1.In the following, we treat {a as unknown deterministic parameters.
The GIC-based rule solving ( 22) is max x (1) ∈G (1)   J (1) (x (1) ) where and When H (1) 1 is chosen, then the estimated position x(1) of the first detected target is the maximizer of (24), while the ML estimate of its gain vector is We now tackle the testing problem in (23).For k ≥ 2, the negative log-likelihood functions under H and Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
respectively, where p,0 = cat a (1)  p , . . ., are the augmented mode matrix and gain vector of the interference under H (k) 0 , while are the augmented mode matrix and gain vector accounting for the target and the interference under 1 .Accordingly, the GICbased rule is max where and are the ML estimate of the augmented gain vectors at the p-th receive antenna under 1 and H (k) 0 , respectively.We now elaborate on (33) to obtain some insights into the decision rule in (32).After plugging (34) into the log-likelihood ratio at the p-th receive antenna, we obtain [57] ln where and Notice here that Ξ (k)   p is the orthogonal projector on the interference subspace spanned by the columns of ) is the orthogonal projector on the part of the column r p contained therein.The penalty term at the p-th receive antenna accounts for the dimension of the subspace along with the signal is projected; indeed, we have Upon plugging ( 35) and ( 38) in (33), we obtain that The matrix ) may be ill-conditioned if the inspected position x (k) and those of some previouslydetected targets present a delay offset less than 1/W between the same transmit and the p-th receive antenna; in this case, we can again proceed as in Remark 2 to compute Π (k)  p (x (k) ).Finally, (32) compares the maximum value of the scoring metric in (39) over all points in G (k) with a threshold.When H (k) 1 is accepted, the estimated position x(k) of the k-th detected target is the argument of the maximum, while the ML estimate of the corresponding gain vector is Problems ( 24) and ( 32) are sequentially solved for k = 1, 2, . . .until no additional target is found or there is no remaining grid point to be inspected at the next iteration or k = K max .If the procedure ends at iteration K, there are K − 1 detected targets if the threshold has not been crossed in the last test, and K otherwise.The penalty factor η can be set to maintain a given probability of false alarm, defined here as P fa = Pr max x(1) ∈G (1)   J (1) (x (1) ) > 0 under The overall procedure is referred to as the MSD-IS and is summarized in Algorithm 1.

B. Solution Based on the MSD-ICM
We discuss here a slight modification of the MSD-IS algorithm, which avoids zero-forcing the subspace signals of the previously-detected targets at each iteration k ≥ 2; indeed, the zero-forcing operation may cause a detrimental noise enhancement.
To proceed, consider the testing problem in (23) when the gain vectors {a if k = 1 then 4: Compute J (1) (x (1) ) from ( 25) 5: else 6: Compute J (k) (x (k) ) from ( 39) or ( 46) for the MSD-IS or MSD-ICM, respectively 7: Compute G (k+1) from ( 21 where γ p,n ≥ 0 is a design parameter accounting for the target strength with respect to the receive/transmit antenna pair (p, n).Accordingly, the negative log-likelihood functions under H respectively, where is the interference-plus-noise covariance matrix.Under the above assumptions, the GIC-based rule is again expressed as in (32), once the scoring metric is computed as where Notice here that (C (k) p ) −1/2 r p is the measured data at the p-th receive antenna after an interference-plus-noise whitening transformation, while ).If this latter matrix is ill-conditioned, we can again proceed as in Remark 2 to compute 1 is accepted, the estimated position x(k) of the k-th detected target is obtained from the maximizer of (32), while the ML estimate of its gain vector is also an estimate of the target strength with respect to the receive/transmit antenna pair (p, n) is The resulting procedure is referred to as the MSD-ICM and is summarized in Algorithm 1.We underline that the MSD-IS and MSD-ICM algorithms share the same decision logic and also provide the same outcome at the first iteration; they only differ in how the interference from the previously-detected targets is handled in subsequent iterations.

C. Mitigation of Small-Scale Localization Errors
The proposed procedures may be sensitive to off-grid target placements; indeed, the estimated interference subspace at iteration k ≥ 2 may not fully contain the echoes of the previouslydetected targets, thus causing an interference spillover.Next, we give a practical fix involving a modified construction of the interference subspace, which extends the heuristics in [45] to the case where the multiple echoes of each target are produced by different probing signals.
To proceed, denote by G ⊃ G a finer uniform grid with element spacing Δg Δ g , by p,n the matrix with columns {s p,n (g)} g∈B (i) .The MSD-IS algorithm is modified as follows.After noise whitening and interference rejection, at the p-th receive antenna, we assume that the n-th echo of the target detected in the i-th iteration belongs to the augmented subspace spanned by the columns of with the understanding that Ξ (i)  p is an all-zero M × M matrix for i = 1.Since E (i) p,n may contain signatures corresponding to non-resolvable delays, it may be ill-conditioned and some care is required to avoid an unnecessary subspace enlargement.Let λ p,n,B (i) be the squared singular values of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
p,n arranged in a decreasing order and let u (i) p,n,1 , . . ., u (i) p,n,B (i) be the corresponding left singular vectors; then, we only retain the left singular vectors corresponding to the U (i) p,n largest singular values, where and is a design parameter.At this point, the matrix Ξ (k) p used at the k-th iteration in (36), (39), and ( 40) is replaced by the orthogonal projector onto the subspace spanned by Similarly, at the k-th iteration of the MSD-ICM algorithm, the following augmented covariance matrix is used in ( 46), ( .

D. Implementation Complexity
Both the MSD-IS and MSD-ICM solutions involve a search over a subset of the grid G at each iteration, whereby their implementation complexity only scales linearly with the size of G, which is a major cost saving with respect to the GIC-based solution in Section III-B.At the first iteration, both algorithms are equivalent; for a given grid point, the computational burden required to obtain the scoring metric in (25) is dominated by the computation of the orthogonal projector Π (1)  p (x (1) ) for p = 1, . . ., P , which costs at most O((M 2 N + MN 2 )P ) flops.Consider now the iteration k ≥ 2. For the MSD-IS algorithm, the computational burden required to obtain the scoring metric in (39) is dominated by the computation of the orthogonal projector Ξ (k)  p for p = 1, . . ., P , which costs at most O((M 2 (k − 1)N + M (k − 1) 2 N 2 )P ) flops; instead, for the MSD-ICM algorithm, the computational burden required to obtain the scoring metric in (46) is dominated by the computation of the square root of the inverse of the matrix C (k) p , for p = 1, . . ., P , which costs at most O(M 3 P )flops.

V. PERFORMANCE ANALYSIS
We consider here a radar employing 2 transmit and 2 receive antennas, with the system parameters in Table I and the geometry in Fig. 3.The targets are located in the inspected region Y = [5700, 6200] × [5700, 6200] m 2 .The time-coded waveforms in (1) are used, with randomly-generated four-phase code sequences, ψ(t) = Π(t/T c ), and φ(t) = 1/T c Π(t/T c ).We denote the k-th target by Q k and model its response along where φ p,n,k is a uniform random phase and ζ p,n,k accounts for the antenna gain, the target radar cross-section, and the propagation effects; also, we define its SNR as Finally, we denote the ratio between the SNR of the targets Q k and Q q as μ 2 k,q = SNR k /SNR q .Three configurations are studied in the following, as shown in Fig. 3. Targets Q 1 and Q 2 are always kept at x 1 = (6000, 6100) m and x 2 = (5800, 5900) m, respectively, so that their echoes are resolvable for each receive/transmit antenna pair, and have μ 2,1 = 0.5.Instead, targets Q 3 , . . ., Q K have a different position and strength in the considered scenarios (more on this infra).Target Q 1 is the target of interest, and its probability of detection (P d ) and the root mean square error (RMSE) in the estimation of its position are used as relevant performance metrics.In particular, under H0 , we say target Q 1 is detected if the event E 1 = {∃ k ∈ {1, . . ., K} : x(k) − x 1 ≤ c/W } is true; accordingly, we have Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Three benchmarks are considered for comparison.The first one is the JDL procedure with successive interference cancellation (JDL-SIC) in [30]; the second one is the GLRT detector operating in the absence of the interfering targets (namely, Q 2 , . . ., Q K ), referred to as GLRT with cleaned data (GLRT-CD); the third one is the GLRT with known locations of the interfering targets (GLRT-KL).All results are obtained by averaging over 2 × 10 4 independent instances and have been obtained by using the software MATLAB 2021a and a machine with two 2.5 GHz Intel(R) Xeon(R) E5-2678V3 CPUs.

A. Scenario 1: Impact of Target Strength
Consider here the Scenario 1 in Fig. 3, where target Q 3 is located at x 3 = (6020, 5740) m.Under this setting, the targets are partially-separable, as their echoes are resolvable for each Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.receive/transmit antenna pair, except for those two along the paths TX 2 − Q 1 − RX 2 and TX 2 − Q 3 − RX 2 .To illustrate the sidelobe masking resulting from the imperfect auto-and crosscorrelation functions of the adopted waveforms, we analyze in Fig. 4 (top), the normalized noise-free output of a standard correlator matched to s 2 (t), when the signal r 2 (t) observed by RX 2 is at the input and μ 3,1 = 4.75.In particular, the left plot shows the output auto-correlation for each target echo produced by s 2 (t) (auto-terms) and the output cross-correlation for each target echo produced by s 1 (t) (cross-terms); also, the middle plot reports the superposition of the auto-and cross-terms for each target; finally, the right plot reports the overall output obtained by superimposing all contributions.It is seen by inspection that the echos of targets Q 1 and Q 3 are not resolvable through the receive/transmit antenna pair (2, 2) and that target Q 2 is masked by the range sidelobes of the target Q 3 .
Next, we illustrate the evolution of the proposed solutions in a single snapshot.The output data plane (i.e., the value of the scoring metric over Y) is reported in Fig. 5 for the MSD-IS, MSD-ICM, and JDL-SIC solutions over four iterations, when SNR 1 = 14 dB and μ 3,1 = 4.75.The following remarks are now in order.
r At the 1-st iteration, all scoring metrics peak at the true position of target Q 3 , so that all algorithms can detect and localize strongest target; however, target Q 3 is also responsible for the presence here of large sidelobes.
r At the 2-nd iteration, an erroneous decision is made by the JDL-SIC algorithm, as there are many areas in Y where the value of its scoring metric is larger than that at the true positions of targets Q 1 and Q 2 ; the reason is that this algorithm neglects at the design stage the sidelobe interference from the first detected target; however, such interference is here strong enough to mask the other two targets.A similar issue is present at the 3-rd iteration of the JDL-SIC algorithm.
r The MSD-IS and MSD-ICM algorithms can instead detect and localize all targets; the reason is that they are explicitly designed to mitigate the interference from the previously detected targets at each iteration, thus overcoming the aforementioned masking effects.Finally, Figs. 6 and 7 report P d and RMSE in the position estimate of target Q 1 versus SNR 1 and μ 3,1 , respectively.The proposed algorithms perform only slightly worse than the single-target benchmark and are not significantly affected by the strength of target Q 3 : this confirms their capability of detecting and localizing multiple targets even in the presence of a power imbalance and partially overlapped subspaces.The MSD-ICM slightly outperforms the MSD-IS here, as the latter suffers from a larger noise enhancement when canceling the interference caused by the previously-detected targets.On the other hand, the performance of the JDL-SIC rapidly degrades as μ 3,1 increases, as this solution neglects the cross-and auto-terms at the output of the MF.

B. Scenario 2: Impact of Target Separation
We now consider the Scenario 2 in Fig. 3, where target Q 3 is located at x 3 = (6100, 5800) m.Under this setting, the target echoes are resolvable for each receive/transmit antenna pair Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.(i.e., we have fully-separable targets).In order to verify this latter point, Fig. 4 (bottom) reports the normalized noise-free output of a standard correlator matched to the waveform s 2 (t) emitted by TX 2 , when the signal r 2 (t) observed by RX 2 is at the input, μ 2,1 = 0.5, and μ 3,1 = 4.75.Unlike the Scenario 1, targets Q 1 and Q 3 are now seen in different range bins through the receive/transmit antenna pair (2, 2); however, the two weaker targets are again masked by the strongest one.
Fig. 8 compares the detection and localization performance for partially-separable (Scenario 1) and fully-separable (Scenario 2) targets when μ 2,1 = 0.5 and μ 3,1 = 2.All algorithms perform better in the presence of fully-separable targets: indeed, when Q 3 is detected before Q 1 , the subsequent interference mitigation step no longer cancels the echoes of Q 1 along the path TX 2 − RX 2 .Remarkably, for the fully-separable case, the proposed solutions only present here an SNR loss of 0.45 dB with respect to the single-target GLRT-CD benchmark when P d = 0.8; instead, the performance of the JDL-SIC still remains much worse.

C. Scenario 3: Impact of Target Number
Finally, we consider the Scenario 3 in Fig. 3, where targets Q 3 , . . ., Q K are randomly displaced in Y with SNR 1 = 13 dB, μ 2,1 = 0.5, and SNR k ∈ [25,50] dB for k = 3, . . ., K. Fig. 9 reports P d and RMSE in the position estimate of target Q 1 versus K.It can be seen that the performance of the proposed algorithms gracefully degrades with an increasing number of interfering targets (as opposed to the JDL-SIC), thus confirming their robustness to the sidelobe masking problem.The performance of the genie-aided solution GLRT-KL also degrades when K is increased; indeed, even when the interference from the other targets belongs to a known subspace, its elimination still may take away a part of the signal of interest.Interestingly, in the considered scenario, the proposed solutions present a relevant performance loss, as compared to the GLRT-KL benchmark, only when the number of overall targets K is above 10.

VI. CONCLUSION
In this article, we have addressed the problem of multi-target JDL in a distributed MIMO radar employing waveforms with imperfect auto-and cross-correlation functions, which can be recast as a composite multiple hypothesis testing problem.After modeling the signals at each receive antenna as the superposition of an unknown number of subspace signals and noise, we have shown that the JDL of multiple targets is tantamount to a joint model order selection and subspace matching problem, which can be tackled by resorting to a GIC.To limit the computational cost, we have then derived two approximate solutions, which convert the joint multi-target search into a sequence of single-target searches.A numerical study has been provided to demonstrate the effectiveness of our solutions under diverse instances of the target strength, separation, and number, also in comparison with another recently deployed algorithm, the single-target benchmark, and an ideal solution that has perfect knowledge of the numbers of interfering targets and their locations.
Future works may consider distributed radar networks with partial/no knowledge of the noise covariance matrices; the development of distributed algorithms for reducing the information exchange among the receive antennas may also be considered.Finally, the derivation of the Cramér-Rao bounds for the considered multi-target joint detection and localization problem still remains an open problem.

Fig. 1 .
Fig. 1.Examples of separable target configurations when K = 2, N = 1, and P = 2.In the top plot, the targets are partially-separable, while in the bottom plot they are fully-separable.

Remark 1 :
C −1/2 p r p is the data measured by the pth receive antenna after a noise whitening transformation, Π p (x 1:K )C −1/2 p r p 2 is the energy of C −1/2 p r p contained into the column span of C −1/2 p

2 )
involves a search over the set G 1:K whose cardinality scales exponentially with K.More specifically, assume that the whitened mode matrix C −1/2 p S p (g) and the whitened measurement C −1/2 p r p are precomputed and stored into a dedicated memory, for any grid point g ∈ G and p = 1, . . ., P .Then, the computation of the orthogonal projector Π p (x 1:K ) requires a pseudoinverse, which costs at most O(MK 2 N 2 ) floating point operations (flops), and a matrix multiplication, which costs O(M 2 KN ) flops; also, the computation of Π p (x 1:K )C flops.Hence, the total cost for computing(19) is at most O((M 2 KN + MK 2 N 2 )P |G 1:K |) flops.

S
p (x(k) ) not contained in the interference subspace, and Π(k)  p (x(k)

( 1 )Algorithm 1 :
p , . . ., a (k−1) p } N p=1 of the previously-detected targets are modeled as independent circularly-symmetric Gaussian vectors; in particular, assume that a (i) p has a diagonal covariance Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Implementation of the MSD-IS/ICM.1: provide G(1) = G, K = 0 and η ≥ 0 2: for k = 1, . . ., K max do 3: p-th receive antenna are

Fig. 4 .
Fig.4.Normalized noise-free output of a correlator matched to the waveform emitted by TX 2 when the signal received by RX 2 is at the input.The left plot shows the output auto-correlation for each target echo produced by s 2 (t) (auto-terms) and the output cross-correlation for each target echo produced by s 1 (t) (cross-terms); the middle plots report the superposition of the auto-and cross-terms for each target; the right plots report the overall output obtained by superimposing all contributions.The top and bottom plots refer to Scenarios 1 and 2 in Fig.3, respectively, when μ 2,1 = 0.5 and μ 3,1 = 4.75.

Fig. 5 .
Fig. 5. Output data plane of the MSD-IS (top), MSD-ICM (middle), and JDL-SIC solutions (bottom) in four consecutive iterations (from left to right).The triangle markers are the true target positions; the circle marker is the estimated position of the target detected (if any) in each iteration.Scenario 1 in Fig. 3 is considered with SNR 1 = 14 dB, μ 2,1 = 0.5, and μ 3,1 = 4.75.