DOA Estimation for Transmit Beamspace MIMO Radar via Tensor Decomposition with Vandermonde Factor Matrix
11 DOA Estimation for Transmit Beamspace MIMORadar via Tensor Decomposition with VandermondeFactor Matrix
Feng Xu,
Student Member, IEEE , Matthew W. Morency,
Student Member, IEEE , and Sergiy A. Vorobyov,
Fellow,IEEE
Abstract —We address the problem of tensor decomposition inapplication to direction-of-arrival (DOA) estimation for transmitbeamspace (TB) multiple-input multiple-output (MIMO) radar.A general 4-order tensor model that enables computationallyefficient DOA estimation is designed. Whereas other tensordecomposition-based methods treat all factor matrices as arbi-trary, the essence of the proposed DOA estimation method is tofully exploit the Vandermonde structure of the factor matrices totake advantage of the shift-invariance between and within differ-ent subarrays. Specifically, the received signal of TB MIMO radaris expressed as a 4-order tensor. Depending on the target Dopplershifts, the constructed tensor is reshaped into two distinct 3-ordertensors. A computationally efficient tensor decomposition methodis proposed to decompose the Vandermonde factor matrices. Thegenerators of the Vandermonde factor matrices are computedto estimate the phase rotations between subarrays, which canbe utilized as a look-up table for finding target DOA. It isfurther shown that our proposed method can be used in a moregeneral scenario where the subarray structures can be arbitrarybut identical. The proposed DOA estimation method requiresno prior information about the tensor rank and is guaranteed toachieve precise decomposition result. Simulation results illustratethe performance improvement of the proposed DOA estimationmethod as compared to conventional DOA estimation techniquesfor TB MIMO Radar.
Index Terms —DOA estimation, Shift-invariance, TB MIMOradar, Tensor decomposition, Vandermonde factor matrix
I. I
NTRODUCTION T HE development of multiple-input multiple-output(MIMO) radar has been the focus of intensive research[1]–[5] over the last decade, and has opened new opportunitiesin target detection and parameter estimation. Many works havebeen reported in the literature showing the applications ofMIMO radar with widely separated antennas [1] or collocatedantennas [2]. Among these applications, direction-of-arrival
This work was supported in part by the Academy of Finland under Grant319822, in part by Huawei, and in part by the China Scholarship Council.This work was conducted while Feng Xu was a visiting doctoral studentwith the Department of Signal Processing and Acoustics, Aalto University.(
Corresponding author: Sergiy A. Vorobyov. )Feng Xu is with the School of Information and Electronics, Beijing Instituteof Technology, Beijing 100081, China, and also with the Department of SignalProcessing and Acoustics, Aalto University, Espoo 02150, Finland. (e-mail:[email protected], feng.xu@aalto.fi).Matthew W. Morency is with the Dept. Microelectronics, School ofElectrical Engineering, Mathematics, and Computer Science, Delft Universityof Technology, Mekelweg 4, 2628 CD Delft, The Netherlands. (e-mail:[email protected]).Sergiy A. Vorobyov is with the Department of Signal Processing andAcoustics, Aalto University, Espoo 02150, Finland. (e-mail: [email protected]). (DOA) estimation [3], [6]–[13] is one of the most fundamentalresearch topics. In this paper, we mainly focus on the DOAestimation problem for MIMO radar with collocated antennas.By ensuring that the transmitted waveforms are orthogonal[14], MIMO radar enables increasing the system’s degree offreedom (DoF), improving the spatial resolution and enhanc-ing the parameter identifiability. The essence behinds theseadvantages is the construction of a virtual array (VA), whichcan be regarded as a new array with larger aperture andmore elements [4], [5]. However, the omnidirectional transmitbeampattern in MIMO radar, resulting from the orthogonalwaveforms, deteriorates the parameter estimation performancesince most of the emitted energy is wasted as comparedto its phased-array counterpart. To tackle this problem, thetransmit beamspace (TB) technique has been introduced [3],[6], [15]. In TB MIMO radar, the transmitted energy canbe focused on a fixed region [3], [6] by using a number oflinear combinations of the transmitted waveforms via a TBmatrix. This benefit becomes more evident when the numberof elements in MIMO radar is large [15]. Specifically, at somenumber of waveforms, the gain from using more waveformsbegins to degrade the estimation performance. The trade-offbetween waveform diversity and spatial diversity implies thatthe performance of DOA estimation in TB MIMO radar canbe further improved with a carefully designed TB matrix.Meanwhile, many algorithms for DOA estimation in MIMOradar have been proposed. These algorithms can be sum-marized in two categories, signal covariance matrix-basedalgorithms [3], [6]–[10] and signal tensor decomposition-basedalgorithms [11]–[13], [16]–[21]. For example, the estimationof target spatial angles can be conducted by multiple signalclassification (MUSIC). The generalization of MUSIC to aplanar array requires a 2-dimension (2-D) spectrum searching[7], and thus suffers from high computational complexity.By exploiting the rotational invariance property (RIP) of thesignal subspace, estimation of signal parameters via rotationalinvariance technique (ESPRIT) [3], [6], [8] can be applied toestimate the target angles without a spectrum searching. TheRIP can be enforced in many ways, e.g., uniformly spacedantennas [8] and the design of TB matrix [3], [6]. To furtherreduce the computational complexity and increase the numberof snapshots, unitary-ESPRIT (U-ESPRIT) has been proposed[10]. Some algorithms like propagator method (PM) have beenstudied [9] to avoid the singular value decomposition (SVD)of the signal covariance matrix. The aforementioned DOA a r X i v : . [ c s . I T ] J a n estimation algorithms are mostly conducted on a per-pulsebasis to update the result from pulse to pulse. They ignorethe multi-linear structure of the received signal in MIMO radarand, therefore, lead to poor performance in low signal-to-noiseratio (SNR) region.The second category, signal tensor decomposition-basedalgorithms, has been proposed to address the problem ofpoor performance in low SNR. In particular, a 3-order tensoris introduced to store the whole received signal for MIMOradar in a single coherent processing interval (CPI). Methodslike high-order SVD (HOSVD) [13], [22] and parallel factor(PARAFAC) analysis [11], [12] can be applied to decomposethe factor matrices. The DOA estimation can be conductedby exploiting the factor matrix with the target angular infor-mation. For example, the widely used alternating least square(ALS) algorithm is a common way of computing the approx-imate low-rank factors of a tensor. These factor matrices canbe used to locate multiple targets simultaneously [12], [16].Although the application of the conventional ALS algorithmimproves the DOA estimation performance for MIMO radar,it usually requires the tensor rank as prior information, andthe computational complexity can be extremely high as theconvergence is unstable.Nevertheless, conventional tensor decomposition methodsare developed for tensors with arbitrary factor matrices. Inarray signal processing, special matrix structure like Toeplitz,Hankel, Vandermonde and columnwise orthonormal [17], [23]may exist in factor matrix when tensor model is applied tocollect the received signal. The Vandermonde structure, as themost common one, can be generated from the applicationof carrier frequency offset, e.g., frequency diversity array(FDA) [24] and orthogonal frequency-division multiplexing(OFDM) waveform [25], or uniformly spaced antennas, e.g.,uniform linear array (ULA) and uniform rectangular array(URA). While conventional tensor decomposition methods areusually designed for tensors with arbitrary factor matrices,the decomposition of a tensor with structured factor matricesdeserves further study as the structured factor matrix maypoint to a novel decomposition method and better uniquenessconditions. This is called constrained tensor decomposition [17], [23]. Moreover, transmit array interpolation is introducedfor MIMO radar with arbitrary array structure [13]. By solvingthe minimax optimization problem regarding interpolationmatrix design, the original transmit array is mapped to avirtual array with desired structure. The DOA estimationbias caused by interpolation errors has also been analyzedin [13]. However, the interpolation technique deteriorates theparameter identifiability, which makes it inappropriate for TBMIMO radar with arbitrary but identical subarrays.In this paper, we consider the problem of tensor decompo-sition in application to DOA estimation for TB MIMO radarwith multiple transmit subarrays. A general 4-order tensormodel that enables computationally efficient DOA estimationis designed. Whereas other tensor decomposition-based meth-ods treat all factor matrices as arbitrary, the proposed DOA Some preliminary ideas that have been extended and developed to thispaper we published in [26], [27]. estimation method fully exploits the Vandermonde structureof the factor matrix to take advantage of the shift-invariancebetween and within different subarrays. In particular, thereceived signal of TB MIMO radar is expressed as a 4-ordertensor. Depending on the target Doppler shifts, the constructedtensor is reshaped into two distinct 3-order tensors. A compu-tationally efficient tensor decomposition method, which can beconducted via linear algebra with no iterations, is proposed todecompose the factor matrices of the reshaped tensors. Then,the Vandermonde structure of the factor matrices is utilized toestimate the phase rotations between transmit subarrays, whichcan be applied as a look-up table for finding target DOA. Itis further shown that our proposed method can be used in amore general scenario where the subarray configurations arearbitrary but identical. By exploiting the shift-invariance, theproposed method improves the DOA estimation performanceover conventional methods, and it has no requirement of priorinformation about the tensor rank. Simulation results verifythat the proposed DOA estimation method has better accuracyand higher resolution.The rest of this paper is organized as follows. Some algebrapreliminaries about tensors and matrices are introduced at theend of Section I. A 4-order tensor model for TB MIMO radarwith uniformly spaced subarrays is designed in Section II. InSection III, the proposed tensor model is reshaped properly toachieve the uniqueness condition of tensor decomposition. TheDOA estimation is conducted by exploiting the shift-invariancebetween and within different subarrays. The parameter iden-tifiability is also analysed. Section IV generalizes the pro-posed DOA estimation method to TB MIMO radar with non-uniformly spaced subarrays, where multiple scales of shift-invariances can be found. Section V performs the simulationexamples while the conclusions are drawn in Section VI.
Notation : Scalars, vectors, matrices and tensors are de-noted by lower-case, boldface lower-case, boldface uppercase,and calligraphic letters, e.g., y , y , Y , and Y , respectively.The transposition, Hermitian transposition, inversion, pseudo-inversion, Hadamard product, outer product, Kronecker prod-uct and Khatri-Rao (KR) product operations are denotedby ( · ) T , ( · ) H , ( · ) − , ( · ) † , ∗ , ◦ , ⊗ , and (cid:12) , respectively, while vec ( · ) stands for the operator which stacks the elements ofa matrix/tensor one by one to a column vector. The notation diag ( y ) represents a diagonal matrix with its elements beingthe elements of y , while (cid:107) Y (cid:107) F and (cid:107) Y (cid:107) are the Frobeniusnorm and Euclidean norm of Y , respectively. Moreover, M × N and M × N denote an all-one matrix of dimension M × N and an all-zero matrix of size M × N , respectively,and I M stands for the identity matrix of size M × M . For B ∈ C M × N , the n -th column vector and ( m, n ) -th element aredenoted by b n and B mn , respectively, while the m -th elementof b ∈ C M × is given by b ( m ) . The estimates of B and b are given by ˆB and ˆb , while the rank and Kruskal-rank of B are denoted by r ( B ) and k B , respectively. To express twosubmatrices of B without the first and last row vector, B and B are applied. If the general form, B can be written as B (cid:44) [ β , β , · · · , β N ] , where β n (cid:44) [1 , z n , z n , · · · , z M − n ] T , i.e., B is a Vandermonde matrix, and z (cid:44) [ z , z , · · · , z N ] T ∈ C N × is the vector of generators. When each element is unique, z is considered to be distinct. A. Algebra Preliminaries for Tensors and Matrices
For an N -th order tensor Y ∈ C I × I ×···× I N , the followingfacts are introduced [16], [28]. Fact 1 (PARAFAC decomposition): The PARAFAC decom-position of an N -th order tensor is a linear combination ofthe minimum number of rank-one tensors, given by Y = L (cid:88) l =1 α (1) l ◦ α (2) l ◦ · · · ◦ α ( N ) l (cid:44) [[ A (1) , A (2) , · · · , A ( N ) ]] (1)where α ( n ) l is the l -th column of A ( n ) with A ( n ) being the n -th factor matrix of size I n × L , and L is the tensor rank. Fact 2 (Uniqueness of PARAFAC decomposition): ThePARAFAC decomposition is unique if all potential factormatrices satisfying (1) also match with ˜A ( n ) = A ( n ) Π ( n ) ∆ ( n ) (2)where Π ( n ) is a permutation matrix and ∆ ( n ) is a diagonalmatrix. The product of ∆ ( n ) , n = 1 , , · · · , N is an L × L identity matrix. Usually, the generic uniqueness condition isgiven by [28]: N (cid:80) n =1 k A ( n ) ≥ L + ( N − . Fact 3 (Mode- n unfolding of tensor): The mode- n unfoldingof a tensor Y ∈ C I × I ×···× I N is denoted by Y ( n ) , which isa matrix of size I · · · I n − I n +1 · · · I N × I n Y ( n ) = (cid:16) A (1) · · · (cid:12) A ( n − (cid:12) A ( n +1) · · · (cid:12) A ( N ) (cid:17) (cid:16) A ( n ) (cid:17) T . (3) Fact 4 (Tensor reshape): The reshape operator for an N -thorder tensor Y ∈ C I × I ×···× I N returns a new M -th ordertensor X ∈ C J × J ×···× J M ( M ≤ N ) with N (cid:81) n =1 I n = M (cid:81) m =1 J m and vec ( Y ) = vec ( X ) , e.g., if J m = I m , m = 1 , , · · · , M − and J M = N (cid:81) n = M I n , the mode- M unfolding of reshaped X is X ( m ) = (cid:16) A (1) (cid:12) · · · (cid:12) A ( M − (cid:17) (cid:16) A ( M ) (cid:12) · · · (cid:12) A ( N ) (cid:17) T . (4) Lemma 1 : For a 3-order tensor Y (cid:44) [[ A (1) , A (2) , A (3) ]] ,where A (1) is a Vandermonde matrix or the KR productof two Vandermonde matrices. The decomposition of Y isgenerically unique if the generators of A (1) are distinct and A (3) is column full rank. Proof
It is purely technical and is given in supplementalmaterial as Appendix A.
Lemma 2 : The following equalities hold true AB = A (cid:12) b T = b T (cid:12) AA (cid:12) b T (cid:12) C = b T (cid:12) A (cid:12) C = A (cid:12) C (cid:12) b T ( A (cid:12) B ) (cid:12) C = A (cid:12) ( B (cid:12) C ) vec ( ABD ) = (cid:0) D T (cid:12) A (cid:1) b ( A ⊗ C ) ( D ⊗ E ) = ( AD ) ⊗ ( CE ) (5) where A ∈ C M × N , C ∈ C Q × N , D ∈ C N × P , E ∈ C N × L and B = diag ( b ) ∈ C N × N .II. TB MIMO R ADAR T ENSOR MODEL
A. TB MIMO Radar with Linear Array
Consider a collocated MIMO radar with M transmit ele-ments and N receive elements. The transmit array is a ULAwith its elements spaced at half the working wavelength awayfrom each other. The receive elements are randomly placedwithin a fixed aperture. Assuming S subarrays are uniformlyspaced at the transmit side and also assuming that each subar-ray contains M elements, the indices of first elements in thosesubarrays are denoted by m s , s = 1 , , · · · , S . Without lossof generality, m s rises uniformly. The steering vectors of theentire transmit array and the first transmit subarray at direction θ can be given by α ( θ ) (cid:44) (cid:2) , e − jπ sin θ , · · · , e − j ( M − π sin θ (cid:3) T and α ( θ ) (cid:44) (cid:2) , e − jπ sin θ , · · · , e − j ( M − π sin θ (cid:3) T , respec-tively. The steering vector of the receive array can be writ-ten as β ( θ ) (cid:44) (cid:104) , e − j πλ x sin θ , · · · , e − j πλ x N sin θ (cid:105) T , where { x n | < x n ≤ D, n = 1 , · · · , N } and D is the aperture ofthe receive array.Accordingly, the transmit and receive steeringmatrices for L targets in { θ l } Ll =1 can be denotedby A (cid:44) [ α ( θ ) , α ( θ ) , · · · , α ( θ L )] and B (cid:44) [ β ( θ ) , β ( θ ) , · · · , β ( θ L )] , respectively, while the steeringmatrix for the first transmit subarray can be given as A (cid:44) [ α ( θ ) , α ( θ ) , · · · , α ( θ L )] . Note that A can alsobe regarded as the submatrix of A with first M rows.In conventional MIMO radar, the received signal at theoutput of the receive array after matched-filtering in matrixform can be modelled as [12]: Y = BΣA T + N , where Σ = diag ( σ ) , σ (cid:44) (cid:2) σ , σ , · · · , σ L (cid:3) T represents the vectorof target radar cross section (RCS) fading coefficients obeyingSwerling I model, and N is the noise residue of size N × M .When the TB technique is introduced [3], [6], the received sig-nal model after matched-filtering of K orthogonal waveforms( K ≤ M ) can be generalized as Y = BΣ ( W H A ) T + N ,where W (cid:44) [ w , w , · · · , w K ] M × K denotes the TB matrix.Hence, the received signal for the s -th transmit subarray, s = 1 , , · · · , S , and the whole receive array can be writtenas Y s = BΣ ( W Hs A s ) T + N s (6)where W s and A s represent the TB matrix and steer-ing matrix for the s -th transmit subarray, respectively, and N s is the noise residue of size N × M . Assume thatthe TB matrix for each subarray is identical, denoted by W ∈ C M × K . Note that since m s rises uniformly, thesteering matrix for the s -th transmit subarray can also beexpressed as A s = A Γ s , where Γ s = diag ( k s ) and k s (cid:44) (cid:2) e − jπ ( m s −
1) sin θ , · · · , e − jπ ( m s −
1) sin θ L (cid:3) T . Substitut-ing this relationship into (6) and vectorizing it, we have y s = (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3) Γ s σ + n s , where n s is the vectorizednoise residue.Considering the Doppler effect, the received signal during q -th pulse in a single CPI, q = 1 , , · · · , Q , can be written as y ( q ) s = (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3) Γ s c q + n ( q ) s (7) Reference SubarrayReference Subarray x - a x i s y-axis c d r d x - a x i s y-axis c d r d subarray x M y M J I ( ) , i j th − ( ) , i j m m Fig. 1. Transmit array configuration for TB MIMO radar with planar array. where c q = σ ∗ ¯c q , ¯c q (cid:44) (cid:2) e i πf qT , e i πf qT , · · · , e i πf L qT (cid:3) T , f l denotes the Doppler shift, T is the radar pulse duration,and n ( q ) s is the vectorized noise residue. Concatenate thereceived signal of S subarrays in q -th pulse, i.e., Y ( q ) (cid:44) (cid:104) y ( q )1 , y ( q )2 , · · · , y ( q ) S (cid:105) KN × S . The compact form can be writtenas Y ( q ) = (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3) (cid:2) c Tq (cid:12) K (cid:3) T + N ( q ) (8)where K (cid:44) [ k , k , · · · , k S ] T S × L and N ( q ) (cid:44) (cid:104) n ( q )1 , n ( q )2 , · · · , n ( q ) S (cid:105) . Note that the l -th column of K represents the phase rotations of l -th target for S subarrays. K can be named as transmit subarray steering matrix .Vectorizing (8), the KN S × vector can be given as z q = (cid:8)(cid:2) c Tq (cid:12) K (cid:3) (cid:12) (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3)(cid:9) L × + r q = (cid:2) K (cid:12) (cid:0) W H A (cid:1) (cid:12) B (cid:3) c q + r q (9)where r q is the vectorized noise residue of N ( q ) . Then,concatenate the received signal of Q pulses, i.e., Z (cid:44) [ z , z , · · · , z Q ] KNS × Q . The compact form can be formulatedas Z = (cid:2) K (cid:12) (cid:0) W H A (cid:1) (cid:12) B (cid:3) C T + R (10)where C (cid:44) [ c , c , · · · , c Q ] T Q × L and R (cid:44) [ r , r , · · · , r Q ] .Similarly, C can be named as Doppler steering matrix sinceeach column denotes the Doppler steering vector for one target(with additional RCS information). According to
Fact 3 , a 4-order tensor
Z ∈ C S × K × N × Q whose matricized version is Z in (10) can be constructed. Denote X (cid:44) W H A , then thistensor can be written as Z = L (cid:88) l =1 κ l ◦ χ l ◦ β l ◦ γ l + R (cid:44) [[ K , X , B , C ]] + R (11)where κ l , χ l , β l , γ l are the l -th columns of K , X , B , C ,respectively, L is the tensor rank, and R is the noise tensorof the same size. B. TB MIMO Radar with Planar Array
Consider the planar array case, as shown in Fig. 1.A URA with M = M x · M y elements spaced at halfthe working wavelength in both directions and a planar array with N randomly spaced elements are applied asthe transmit and receive array, respectively. The transmitsteering vector can be given as α ( θ, ϕ ) = u ( θ, ϕ ) ⊗ v ( θ, ϕ ) , where u ( θ, ϕ ) (cid:44) (cid:2) , e − jπu , · · · , e − j ( M y − πu (cid:3) T , v ( θ, ϕ ) (cid:44) (cid:2) , e − jπv , · · · , e − j ( M x − πv (cid:3) T , u (cid:44) sin ϕ sin θ , v (cid:44) sin ϕ cos θ , and ( θ, ϕ ) is the pair of azimuth and elevationof a target. The steering vector of the receive array can be writ-ten as β ( θ, ϕ ) (cid:44) (cid:104) , e − j πλ ( x v + y u ) , · · · , e − j πλ ( x N v + y N u ) (cid:105) T ,where { ( x n , y n ) | < x n ≤ D x , < y n ≤ D y } are the coordi-nates of the receive elements, and D x , D y denote the aperturesin two directions, respectively.Accordingly, assume S = I · J transmit subarrays are uni-formly spaced at the transmit side, which can be overlapped ornot. Each of them contains M = M x · M y elements. The firstsubarray is selected as the reference subarray. For L targets in { ( θ l , ϕ l ) } Ll =1 , the transmit and receive steering matrices canbe generalized as A (cid:44) [ α ( θ , ϕ ) , α ( θ , ϕ ) , · · · , α ( θ L , ϕ L )] and B (cid:44) [ β ( θ , ϕ ) , β ( θ , ϕ ) , · · · , β ( θ L , ϕ L )] , respectively.Note that the transmit array is a URA, thus, we have A = U (cid:12) V , where U (cid:44) [ u ( θ , ϕ ) , u ( θ , ϕ ) , · · · , u ( θ L , ϕ L )] and V (cid:44) [ v ( θ , ϕ ) , v ( θ , ϕ ) , · · · , v ( θ L , ϕ L )] . Similarly,the steering vector of the reference transmit subarray can bewritten as α ( θ, ϕ ) = u ( θ, ϕ ) ⊗ v ( θ, ϕ ) , where u ( θ, ϕ ) and v ( θ, ϕ ) contain the first M y and M x elements in u ( θ, ϕ ) and v ( θ, ϕ ) , respectively. The steering matrix for the referencetransmit subarray can be denoted by A = U (cid:12) V , where U and V are the submatrices of U and V that consist ofthe first M y and M x rows, respectively.For the ( i, j ) -th subarray (or equivalently, for the s -th trans-mit subarray where s = ( j − I + i ), the index of first elementis denoted by ( m i , m j ) , i = 1 , , · · · , I, j = 1 , , · · · , J .Both m i and m j rise uniformly. The steering matrix forthe ( i, j ) -th subarray can be given as A ij = U j (cid:12) V i ,where U j = U Γ j , V i = V Γ i , Γ j = diag ( h j ) , Γ i = diag ( d i ) , vectors h j (cid:44) (cid:2) e − jπ ( m j − u , · · · , e − jπ ( m j − u L (cid:3) T and d i (cid:44) (cid:2) e − jπ ( m i − v , · · · , e − jπ ( m i − v L (cid:3) T indicate thephase rotations for L targets in two directions, respectively.Generalizing (7), the received signal after matched-filteringfor the s -th transmit subarray and the whole receive array in q -th pulse can be written as y ( q ) s = (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3) Γ i Γ j c q + n ( q ) s . (12)Similarly, the concatenation of the received signal y ( q ) s forall S subarrays in q -th pulse can be expressed as Y ( q ) = (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3) (cid:0) c Tq (cid:12) H (cid:12) ∆ (cid:1) T + N ( q ) (13)where H (cid:44) [ h , h , · · · , h J ] T J × L and ∆ (cid:44) [ d , d , · · · , d I ] T I × L . Proof of (13) is purely technicaland is given in supplemental material as Appendix B. Then z q = vec ( Y ( q ) ) can be formulated as z q = (cid:2) H (cid:12) ∆ (cid:12) (cid:0) W H A (cid:1) (cid:12) B (cid:3) c q + r q . (14)After concatenating z q in the same way as (10), the receivedsignal of Q pulses in the URA case can be written as Z = (cid:2) H (cid:12) ∆ (cid:12) (cid:0) W H A (cid:1) (cid:12) B (cid:3) C T + R . (15) It is interesting that, (15) can be directly obtained from (10)by replacing K with H (cid:12) ∆ . Hence, H (cid:12) ∆ can be regarded asthe transmit subarray steering matrix for URA. Using Fact 3 ,a 5-order tensor Z whose matricized version is Z in (15) canbe constructed as Z = L (cid:88) l =1 η l ◦ δ l ◦ χ l ◦ β l ◦ γ l + R (cid:44) [[ H , ∆ , X , B , C ]] + R (16)where η l and δ l are the l -th columns of H and ∆ , respectively.Note that since all subarrays are uniformly spaced, K , ∆ and H are Vandermonde matrices and their vectors ofgenerators can be respectively denoted by ω (cid:44) (cid:2) e − jπ ∆ m sin θ , · · · , e − jπ ∆ m sin θ L (cid:3) T ω x (cid:44) (cid:2) e − jπ ∆ mx v , · · · , e − jπ ∆ mx v (cid:3) T ω y (cid:44) (cid:2) e − jπ ∆ my u , · · · , e − jπ ∆ my u L (cid:3) T (17)where the step sizes ∆ m = m s +1 − m s , ∆ m x = m i +1 − m i ,and ∆ m y = m j +1 − m j . We assume ω , ω x and ω y are distinct,which means that multiple targets are spatially distinct.III. DOA E STIMATION VIA T ENSOR D ECOMPOSITIONWITH V ANDERMONDE F ACTOR M ATRIX
We have shown that the received signal of TB MIMO radarwith transmit subarrays can be formulated as a high-ordertensor. It is useful to point out that (11) and (16) are identicalif the idea of tensor reshape is applied and K is replaced by H (cid:12) ∆ . Hence, a general 4-order tensor model can be used toexpress the received signal for TB MIMO radar with uniformlyspaced subarrays, given by Z (cid:44) [[ G , X , B , C ]] + R (18)where G ∈ C S × L is the transmit subarray steering matrix .Essentially, G can be interpreted as the result of element-wise spatial smoothing between the transmit elements. A newdimension is extended to express the phase rotations betweentransmit subarrays in the tensor model, which matches withthe derivations in (11) and (16).The tensor decomposition of Z can be regarded as the con-strained tensor decomposition, since one of the factor matricesis structured by the regular array configuration. Generally, theALS algorithm can be applied to decompose such a tensor.However, the convergence of the ALS algorithm heavily relieson the determination of tensor rank, which is an NP-hardproblem. The Vandermonde structure of the factor matrix isignored. The number of iterations in ALS algorithm is alsouncertain, which may lead to high computational complexity.In the literature [17], [18], [23], the uniqueness conditionof the tensor decomposition with special-structured factormatrices, e.g., Toeplitz, Hankel, Vandermonde and column-wise orthonormal, has been investigated. The structured factormatrix may change the uniqueness condition and, therefore,point to some new tensor decomposition methods.In this section, we mainly focus on the tensor decompositionwith Vandermonde factor matrix in application to DOA estima-tion for TB MIMO radar with uniformly spaced subarrays. Acomputationally efficient DOA estimation method is proposed and we discuss the application of the proposed method forboth linear and planar arrays.To begin with, a 3-order tensor F (cid:44) [[ G , ( X (cid:12) B ) , C )]] canbe reshaped from (18) (see Fact 4 ), whose mode-3 unfoldingis F (3) = ( G (cid:12) X (cid:12) B ) C T . Considering only the q -thpulse, F (3) is generically identical to that in (9) for lineararray or (14) for planar array. In other words, the signal co-variance matrix-based DOA estimation methods like MUSICand ESPRIT can be conducted by using R = 1 /Q F (3) F H (3) as the signal covariance matrix. Meanwhile, note that G iseither a Vandermonde matrix or the KR product of a pair ofVandermonde matrices. Thus, Lemma 1 can be applied toconduct a tensor decomposition-based DOA estimation if thesecond precondition is satisfied, i.e., r ( C ) = L .Take a ULA for example, let G = K . The SVD of F (3) is denoted by F (3) = UΛV H , where U ∈ C SKN × L , Λ ∈ C L × L , and V ∈ C Q × L . According to Lemma 1 , there mustexist a nonsingular matrix M of size L × L such that UM = K (cid:12) X (cid:12) B (19)or equivalently, U M = K (cid:12) X (cid:12) B , U M = K (cid:12) X (cid:12) B (20)where submatrices U = (cid:2) I KN ( S − , KN ( S − × KN (cid:3) U and U = (cid:2) KN ( S − × KN , I KN ( S − (cid:3) U are truncated from rowsof U . Since K is a Vandermonde matrix, K = KΩ , where Ω = diag ( ω ) . Substitute it into (20) to obtain U M = U MΩ . (21)Note that M and Ω are both full rank, U = U (cid:0) MΩM − (cid:1) . After the eigenvalue decomposition (EVD)of the matrix U † U , ω can be estimated as the vectorof eigenvalues and M is the matrix of the correspondingeigenvectors. Then, the target DOA can be computed by ˆ ω ( l ) = e − jπ ∆ m sin ¯ θ l ∆ m sin ¯ θ l − ∆ m sin ¯ θ (cid:48) l = ± k (22)where k ∈ (cid:0) − ∆ m , ∆ m (cid:1) is an integer, ¯ θ l is the true direction,and ¯ θ (cid:48) l denotes the potential grating lobes when ∆ m ≥ .The estimation of ω x and ω y for planar array is straight-forward, which can also be found in the proof of Lemma 1 .Consequently, ˆ u l and ˆ v l can be determined by (cid:26) ˆ ω y ( l ) = e − jπ ∆ my ¯ u l ∆ m y ¯ u l − ∆ m y ¯ u (cid:48) l = ± k y , (cid:26) ˆ ω x ( l ) = e − jπ ∆ mx ¯ v l ∆ m x ¯ v l − ∆ m x ¯ v (cid:48) l = ± k x (23)where k y ∈ (cid:16) − ∆ my , ∆ my (cid:17) and k x ∈ (cid:16) − ∆ mx , ∆ mx (cid:17) areintegers, respectively, ¯ u l and ¯ v l indicate the DOA informationof l -th target, while ¯ u (cid:48) l and ¯ v (cid:48) l correspond to the potentialgrating lobes. Since u l (cid:44) sin ϕ l sin θ l and v l (cid:44) sin ϕ l cos θ l ,the pair of (ˆ θ l , ˆ ϕ l ) can be denoted by ˆ θ l = arctan (cid:18) ¯ u l ¯ v l (cid:19) , ¯ ϕ l = arcsin (cid:18)(cid:113) ¯ u l + ¯ v l (cid:19) . (24)The process in (19)-(21) can be regarded as the gen-eralized ESPRIT method [29]. Compared to other tensor-decomposition based methods like PARAFAC, the Vander- monde structure of the factor matrix is exploited and the com-putational complexity is reduced significantly. No iterationsare required and the convergence is guaranteed.However, the precondition r ( C ) = L must be satisfied. Insome applications regarding target detection, it may happenthat two targets with similar Doppler shifts exist. Under thiscircumstance, two column vectors in C are considered to belinearly dependent. The rank deficiency problem limits theapplication of this computationally efficient DOA estimationmethod. Besides, the spatial ambiguity problem further re-stricts the placement of transmit elements. The distance ofphase centers between two adjacent subarrays should be nomore than half the working wavelength. The array aperture islimited. To tackle the problem of rank deficiency and obtain ahigher spatial resolution, (18) is reshaped by squeezing B and C into one dimension. The third factor matrix B (cid:12) C , as theKR product of a Vandermonde matrix and an arbitrary matrix,generically has rank min( QN, L ) . Two targets with identicalDoppler shift can be resolved, while the grating lobes can beeliminated by comparing the estimation result originated from G to the distinct target angular information obtained by X [21], [27]. A. Proposed Computationally Efficient DOA EstimationMethod for TB MIMO Radar with Uniformly Spaced TransmitSubarrays
Consider the noise-free version of (18), a 3-order tensor T (cid:44) [[ G , X , ( B (cid:12) C )]] can be reshaped. The mode-3 unfold-ing of T is given by T (3) = ( G (cid:12) X ) ( B (cid:12) C ) T (25)where G , X , and B (cid:12) C are the three factor matrices, respec-tively. The receive steering matrix and Doppler steering matrixare squeezed into one dimension. Note that the generators of G are distinct, the directions of all targets are unique withor without the existence of grating lobes. Hence, the thirdfactor matrix B (cid:12) C is column full rank [17]. Lemma 1 holds for tensor T . In the following, we develop methods forDOA estimation in TB MIMO radar with uniformly spacedsubarrays via the decomposition of T for linear and planararray sequentially.
1) ULA:
Let G = K . According to Lemma 1 , thedecomposition of T is unique. To obtain the factor matriceswith target DOA information, denote the SVD of T (3) as T (3) = UΛV H , where U ∈ C SK × L , Λ ∈ C L × L , and V ∈ C NQ × L . A nonsingular matrix E of size L × L satisfies UE = K (cid:12) X . (26)Owing to the operator of the KR product, we can write U E = K (cid:12) X , U E = K (cid:12) X (27)where U = (cid:2) I K ( S − , K ( S − × K (cid:3) U and U = (cid:2) K ( S − × K , I K ( S − (cid:3) U are truncated from rows of U ,respectively. Substitute K = KΩ into (27) to obtain U E = U EΩ . (28) Although there exists no deterministic formula for the rank of the KRproduct of a Vandermonde matrix and an arbitrary matrix, it is genericallyfull rank. See Appendix A.
Algorithm 1
DOA Estimation for 1-D TB MIMO Radar withUniformly Spaced Transmit Subarrays
Require:
Signal tensor
Z ∈ C S × K × N × Q from (11) Ensure:
Targets DOA information { θ l } Ll =1 Reshape Z into a 3-order tensor T ∈ C S × K × NQ , wherethe mode-3 unfolding of T is given by (25); Compute the SVD of the matrix T (3) = UΛV H ; Formulate two submatrices U , U satisfying (27); Calculate the EVD of the matrix U † U ; Estimate ˆ θ l via (22), which contains grating lobes; Construct χ l via (29); Define ˜W (cid:44) W − W (cid:48) , W (cid:48) (cid:44) (cid:2) χ l , K × (M − (cid:3) T ; Build a polynomial via F ( z l ) (cid:44) p H ( z l ) ˜W ˜W H p ( z l ) ,where p ( z l ) (cid:44) (cid:2) , z l , · · · , z lM − (cid:3) T and z l (cid:44) e − jπ sin θ l ; Compute the roots of the polynomial F ( z l ) and select theone closest to the unit circle as ˆ z l ; Estimate θ l via ˆ θ l = arcsin (cid:16) j ln(ˆ z l ) π (cid:17) ; Compare the results in step 5 and step 10; return { θ l } Ll =1 .Since E and Ω are both full rank, U = U (cid:0) EΩE − (cid:1) .The connection between U † U and EΩE − is revealed. Fromthe EVD of U † U , the generators ω can be estimated with E being the matrix of the corresponding eigenvectors. Then, (cid:110) ˆ θ l (cid:111) Ll =1 can be computed by (22).Note that (cid:16) κ Hl κ Hl κ l ⊗ I K (cid:17) ( κ l ⊗ χ l ) = χ l and κ Hl κ l = S ,the compact form of χ l is given as χ l = 1 /S (cid:0) κ Hl ⊗ I K (cid:1) Ue l . (29)Equation (29) provides an estimation of each column vectorof X . Given W as prior information, a polynomial rootingmethod [21] can be applied to estimate the unambiguous { θ l } Ll =1 in A independently. Instead of exploiting the signalsubspace shift-invariance of the transmit subarray steeringmatrix, the method in [21] focuses on the Vandermondestructure of A within a single subarray and reveals the re-lationship between the TB MIMO radar transmit beampatternand the generalized sidelobe canceller (GSC). Consequently,the estimation results originated from K and X both providethe target angular information. The grating lobes can beeliminated by comparing the results to each other. Note that(29) is conducted column by column, the angles are pairedautomatically before comparison.An outline of the proposed method for the DOA estimationin TB MIMO radar with linear array is given as Algorithm 1 .
2) URA:
First, substitute G = H (cid:12) ∆ into (25). Similarto (26), the SVD of T (3) is T (3) = UΛV H , and there is anonsingular matrix E ∈ C L × L such that UE = H (cid:12) ∆ (cid:12) X . (30) Considering the KR product, the Vandermonde structure ofboth H and ∆ is exploited via U E = H (cid:12) ∆ (cid:12) X = (cid:0) H (cid:12) ∆ (cid:12) X (cid:1) Ω y = U EΩ y U E = H (cid:12) ∆ (cid:12) X = (cid:0) H (cid:12) ∆ (cid:12) X (cid:1) Ω x = U EΩ x (31)where Ω y = diag ( ω y ) , Ω x = diag ( ω x ) , U , U , U and U are the submatrices truncated from rows of U , i.e., U = (cid:2) I IK ( J − , IK ( J − × IK (cid:3) UU = (cid:2) IK ( J − × IK , I IK ( J − (cid:3) UU = (cid:0) I J ⊗ (cid:2) I K ( I − , K ( I − × K (cid:3)(cid:1) UU = (cid:0) I J ⊗ (cid:2) K ( I − × K , I K ( I − (cid:3)(cid:1) U . (32)Like (28), the vectors ω y and ω x can be estimated as thecollections of eigenvalues of U † U and U † U , respectively,and E is the matrix of the corresponding eigenvectors. Then,the possible pairs of ( θ l , ϕ l ) can be computed by (23)-(24).To eliminate the grating lobes, the relationship between theTB MIMO radar transmit beampattern and the GSC is appliedagain to estimate the target DOA in 2-D case [27]. Specifically, (cid:16) κ Hl κ Hl κ l ⊗ I K (cid:17) ( κ l ⊗ χ l ) = χ l and κ Hl κ l = S still hold byreplacing κ l with h l (cid:12) δ l . Hence, each column of X can berestored by χ l = 1 /S (cid:2) ( h l (cid:12) δ l ) H ⊗ I K (cid:3) Ue l . (33)Note that the TB matrix W is given as a prior information, χ l = W H α ( θ l , ϕ l ) can be rewritten as K different linearequations χ l ( k ) = w Hk α ( θ l , ϕ l ) , k = 1 , , · · · , K (34)or equivalently, p Hk α ( θ l , ϕ l ) = 0 , where p k (cid:44) w k − [ χ l ( k ) , × ( M − ] . It can be seen that the linear equationsin (34) hold if and only if (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) P H a (ˆ θ l , ˆ φ l ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0 , where P = W − W , W (cid:44) [ χ l , K × ( M − ] T . Therefore, theestimation of a pair ( θ l , ϕ l ) can be found by solving thefollowing convex optimization problem [27] min (ˆ θ l , ˆ φ l ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) P H ( u (ˆ θ l , ˆ ϕ l ) ⊗ v (ˆ θ l , ˆ ϕ l )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (35)whose structure is similar with the TB MIMO radar transmitbeampattern. After obtaining the pair (ˆ u l , ˆ v l ) , the distincttarget DOA can be computed by (24). The computation isconducted via (33) column by column, hence the independentestimates of target DOA from G and X are paired automati-cally. By comparing them, the grating lobes can be mitigated.The primary procedures for the DOA estimation in TBMIMO radar with planar array is summarized as Algorithm 2 . B. Parameter Identifiability
As mentioned in
Fact 2 , the generical uniqueness conditionof tensor decomposition for a high-order tensor is given as N (cid:80) n =1 k A ( n ) ≥ L + ( N − , where N is the tensor order and L is the tensor rank. The upper bound of the tensor rank, i.e., themaximum number of targets that can be resolved, rises with Algorithm 2
DOA Estimation for 2-D TB MIMO radar withUniformly Spaced Transmit Subarrays
Require:
Signal Tensor
Z ∈ C J × I × K × N × Q from (16) Ensure:
Targets DOA information { θ l } Ll =1 and { ϕ l } Ll =1 Reshape Z into a 3-order tensor T ∈ C S × K × NQ , wherethe mode-3 unfolding of T is given by (25); Compute the SVD of the matrix T (3) = UΛV H ; Formulate four submatrices U , U , U , U via (32); Calculate the EVD of the matrices U † U and U † U ; Estimate the pair (ˆ u l , ˆ v l ) via (23) and compute the pair (ˆ θ l , ˆ φ l ) via (24); Construct χ l via (33), where E is the matrix of thecorresponding eigenvectors in step 4; Build the matrix P and solve the minimization problem(35) to obtain the pair (ˆ u l , ˆ v l ) ; Compute the unambiguous pair (ˆ θ l , ˆ ϕ l ) via (24); Compare the results in step 5 and step 8; return { θ l } Ll =1 and { ϕ l } Ll =1 .the increase of tensor order at the level of Kruskal-rank ofa matrix. However, if the factor matrix has special structure,the uniqueness condition is changed. An example of tensordecomposition with Vandermonde factor matrix is describedin Lemma 1 . It can be observed that the maximum number oftargets that can be resolved by (18) is determined by the pre-conditions. The first precondition is that the first factor matrixof the tensor, as a Vandermonde matrix or the KR productof two Vandermonde matrices, must have distinct vector ofgenerators. In this paper, we assume this condition holds sinceit means that each of the targets posses a unique direction,which is reasonable regarding DOA estimation problem.The second precondition requires that the third factor matrixhas rank L . Note that the Doppler steering vectors of any twotargets with the same Doppler shift are linearly dependent witha scale difference determined by the target RCS. Two scenariosare discussed next.When the target Doppler shifts are unique, Lemma 1 canbe applied directly and the tensor F can be used. To ensurethe uniqueness decomposition, it is required that [17] min(( S − KN, Q ) ≥ L. (36)In MIMO radar, the number of pulses during a single CPI isusually large. Thus, the maximum number of targets that canbe resolved is generically ( S − KN , which is better thanthat in Fact 2 [17]. However, the size of the transmit arrayis confined since the distance between phase centers of twoadjacent subarrays must be no more than half the workingwavelength to avoid the spatial ambiguity. This restrictiondegrades the spatial resolution and also raises the difficultyof the physical implementation of the array.When there are at least two targets that have identicalDoppler shift, tensor T is used to ensure the second precondi-tion. The receive steering matrix is squeezed together with theDoppler steering matrix to distinguish targets with identicalvelocity. Although the rank of a specific tensor remains the same when it is reshaped, it was proved that different reshapemanners are not equivalent from the performance point of view[30]. In our case, it means that the identifiability is changedand, therefore, the uniqueness condition of decomposition for T requires min (( S − K, N Q ) ≥ L. (37)The usage of T is more appropriate for the general case,since the rank deficiency problem caused by identical targetDoppler shift is solved. Additionally, by reshaping tensor Z into tensor T , the angular information can be estimatedindependently from G and X . The unambiguous estimationresult from X provides a second estimation of the target DOAand can be used to eliminate the grating lobes. Thus, it canbe concluded that the use of the tensor model T has at leastthe above two advantages. The maximum number of targetsthat can be resolved is reduced to ( S − · K . To improve theparameter identifiability, the increase of number of transmitsubarrays or transmit waveforms is worth considering.IV. A RBITRARY BUT I DENTICAL S UBARRAYS WITH M ULTIPLE S CALES OF S HIFT - INVARIANCES
In previous section, we have assumed that the transmitsubarrays are uniformly spaced to obtain a Vandermondestructure in the factor matrix of designed tensor. However, suchconstraint on subarray structure can be relaxed. The placementof all subarrays needs not be uniform, while the configurationwithin a single subarray can be arbitrary. The tensor modelin (18) is applicable for TB MIMO radar with any arbitrarybut identical subarrays, since the extended factor matrix thatrepresents the phase rotations between transmit subarrays ismerely determined by the coordinates of the transmit subarrayphase centers. The difference is that the array configurationvaries the structure of the factor matrix, which may causeextra steps to recover the target DOAs. A typical example hasbeen given earlier where the unambiguous spatial informationin X is exploited to eliminate the cyclic ambiguity in G .In the following, we discuss two general cases that the trans-mit array with multiple scales of shift-invariances is placed ona lattice and explain the use of the proposed computationallyefficient DOA estimation method in both scenarios. A. Generalized Vandermonde Matrix
Note that the Vandermonde structure of the steering matrix G is linked to the phase rotations between the transmitsubarrays, and it is exploited in a look up table for findingtarget DOAs. The Vandermonde structure is only a specialcase leading to phase rotation. Indeed, take, for example alinear array with its elements placed on a lattice, where alllattice cells are enumerated sequentially. The indices of theelements form a counted set of increasing positive integers. Itcan be shown that m s must be a subset of this set, since thefirst element of each subarray corresponds to a unique latticecell. Hence, m s may increase uniformly or non-uniformly.From (17), it can be observed that m s determines K .When it rises uniformly, K should be a Vandermonde matrix. This has been derived in Section II, and we can find that the shift-invariance between different subarrays is related to the step size of m s , ormore specifically, the coordinates of the transmit subarray phase centers. Determined by the step size of m s , i.e., ∆ m , the configurationof adjacent subarrays can be partly overlapped ( ∆ m < M ) ornon-overlapped ( ∆ m ≥ M ). The proposed DOA estimationmethod in last section can be used directly.In the case when m s rises non-uniformly, let us consideras an example the tensor model (11) with S = 7 subarraysand m s = { , , , , , , } . Each subarray contains threeelements, therefore, the original transmit array is a ULA with M = 11 elements. Then, K is a generalized Vandermondematrix [31], which can be written as K (cid:44) [ z , · · · , z L ] , where z l (cid:44) [1 , z l , z l , z l , z l , z l , z l ] T and z l (cid:44) e − jπ sin θ l .The idea of multiple invariance ESPRIT [32] is introducedto conduct the DOA estimation with non-uniformly spacedtransmit subarrays. Consequently, K can be interpreted as thecombination of a set of submatrices K ( sub ) denoting differentsub-ULAs associated with various shift-invariances, i.e., K ( sub ) (cid:44) (cid:20)(cid:16) K (1 , (cid:17) T , (cid:16) K (2 , (cid:17) T , (cid:16) K (1 , (cid:17) T (cid:21) T K (1 , (cid:44) (cid:104) z (1 , , · · · , z (1 , L (cid:105) K (2 , (cid:44) (cid:104) z (2 , , · · · , z (2 , L (cid:105) K (1 , (cid:44) (cid:104) z (1 , , · · · , z (1 , L (cid:105) (38)where z (1 , l is selected from z l with m s = { , , } , z (2 , l isselected from z l with m s = { , , } , and z (1 , l is selectedfrom z l with m s = { , , , , } . In other words, K (1 , isa submatrix of K that consists of first three rows with shift-invariance ∆ m = 1 . The other two submatrices are analogous.Note that (38) is not the only subarray construction method,but it contains all transmit subarrays with a minimal distinctshift-invariance set ∆ = { ∆ m | ∆ m = 1 , } .Substituting (38) to (25), we can write T ( sub )(3) = (cid:16) K ( sub ) (cid:12) X (cid:17) ( B (cid:12) C ) T . (39)Its SVD is given as T ( sub )(3) = U ( sub ) Λ ( sub ) (cid:0) V ( sub ) (cid:1) H . It canbe observed that Lemma 1 holds for (39). By constructing K ( sub ) , a new transmit subarray steering matrix that consistsof several Vandermonde submatrices can be introduced. Toexploit the Vandermonde structure, an extra row selection mustbe applied. Taking K (1 , , for example, we can generalize (26)to obtain K (1 , (cid:12) X = U (1 , E (40)where U (1 , is truncated from U ( sub ) in the same way as K (1 , from K ( sub ) . Thus, each column of K (1 , can beestimated. The estimates of K (1 , and K (2 , can be obtainedsimilarly. It is worth noting that if ∆ m > , the problemof grating lobes may still occur when recovering θ l from U ( N ∆ m , ∆ m ) , where N ∆ m represents the number of subarrayswhose shift-invariance is determined by ∆ m . Usually, theunambiguous spatial information in X can be exploited toeliminate the potential grating lobes. However, it requires eachsubarray to be dense ULA, which restricts the aperture of thetransmit subarray, and therefore, the spatial resolution.Note that the generators of the Vandermonde submatricesin K ( sub ) provide the target DOA information at different exponential levels, i.e., z ∆ m l . Based on this, a polynomialfunction is designed to estimate the target DOA without usingthe second factor matrix X .For every possible shift-invariance ∆ m , denote a (∆ m ) l (cid:44) (cid:104) κ (1 , ∆ m ) Tl , · · · , κ ( N ∆ m , ∆ m ) Tl (cid:105) T b (∆ m ) l (cid:44) (cid:104) κ (1 , ∆ m ) Tl , · · · , κ ( N ∆ m , ∆ m ) Tl (cid:105) T . (41)To illustrate (41), consider the array structure in (38). When ∆ m = 1 , there are two different submatrices/sub-ULAs,i.e., N = 2 , a (1) l = (cid:104) κ (1 , Tl , κ (2 , Tl (cid:105) T = (cid:2) , z l , z l , z l (cid:3) T ,and b (1) l = (cid:104) κ (1 , Tl , κ (2 , Tl (cid:105) T = (cid:2) z l , z l , z l , z l (cid:3) T . When ∆ m = 2 , only one submatrix/sub-ULA exists, i.e., N = 1 , a (2) l = κ (1 , l = (cid:2) , z l , z l , z l (cid:3) T and b (2) l = κ (1 , l = (cid:2) z l , z l , z l , z l (cid:3) T . Also, the following constraint should besatisfied a (∆ m ) l z ∆ m l = b (∆ m ) l . (42)It is proved in [31] that (42) can be achieved by rooting thepolynomial function f ( z l ) (cid:44) (cid:88) ∆ m ∈ ∆ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) a (∆ m ) l z ∆ m l − b (∆ m ) l (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) F (43)as long as two coprime numbers can be found in the shift-invariance set ∆ . By definition of z l , the root nearest to the unitcircle should be chosen as ˆ z l , which finally estimates the targetDOA as ˆ θ l = arcsin (cid:16) j ln(ˆ z l ) π (cid:17) . The construction of K ( sub ) enables the use of T in a more general scenario. The transmitsubarrays can be organized in a non-uniform way. If the shift-invariance set ∆ contains a pair of coprime integers, theproblem of spatial ambiguity can be solved with no limitationon the transmit subarray structure. Hence, the structures of thetransmit subarrays can be arbitrary but identical.An outline of the proposed DOA estimation method for TBMIMO radar with non-uniformly spaced arbitrary but identicalsubarrays is summarized in Algorithm 3 . Remarks : Multiple scales of shift-invariances can also befound in a Vandermonde matrix K [32]. A simple way tobuild K ( sub ) is to concatenate the submatrices of K , which,respectively, consist of the odd rows, even rows and all rows.Hence, Algorithm 3 is also applicable for TB MIMO radarwith uniformly spaced transmit subarrays. Since the manifoldof the subarray does need not to be dense ULA, it is possibleto place the transmit array on a larger lattice to obtain a higherspatial resolution. If some elements in a transmit subarrayare broken, a useful solution is to disable the elements inother subarrays accordingly to keep the manifolds identical.Moreover, we can select part of the elements in all subarraysto fulfill other purposes like communication in joint radar-communication system for example [33]. These remarks canbe extended to the case of planar array.
B. Multiscale Sensor Array
Another case when multiple scales of shift-invariances canbe found is called multiscale sensor array [34]–[36]. Generally,
Algorithm 3
DOA Estimation for TB MIMO radar with Non-Uniformly Spaced Arbitrary but Identical Subarrays
Require:
Signal Tensor
Z ∈ C S × K × N × Q from (18) Ensure:
Targets DOA information { θ l } Ll =1 Construct a new matrix K ( sub ) in (38), which can be di-vided into several Vandermonde submatrices and containsall transmit subarrays; Update the transmit subarray steering matrix and reshape Z into a 3-order tensor T ∈ C S (cid:48) × K × NQ ; Compute the SVD of the matrix T ( sub )(3) = U ( sub ) Λ ( sub ) V ( sub ) H ; Estimate each Vandermonde submatrix in K ( sub ) sequen-tially via (26)-(28); Build two vectors a (∆ m ) l and b (∆ m ) l from (41) for eachcolumn of estimated K ( sub ) ; Compute the roots of the polynomial function in (43) andselect the root nearest to the unit circle as ˆ z l ; Estimate ˆ θ l = arcsin (cid:16) j ln(ˆ z l ) π (cid:17) ; return { θ l } Ll =1 .a URA can be regarded as a 2-level multiscale sensor arraywith different scales of shift-invariances. As shown in Fig. 1,the generation process of such a URA contains two steps.First, consider a single subarray composed of M elementsas the reference subarray. Let I replica subarrays be placeduniformly across the x -axis, which form a larger subarray ata higher level. Then, J copies of this higher level subarrayare organized uniformly across the y -axis. Combining themtogether, a URA is constructed. Note that in this specificcase, I subarrays at level-1 are non-overlapped, while their J counterparts at level-2 are partly overlapped.From (16), it is clear that the transmit subarray steeringmatrices for subarrays at level-1 and level-2 are ∆ and H ,respectively. If the URA itself is repeated and spatially movedto other arbitrary but known locations, a new array thatyields a larger spatial aperture is created. In this way, an R -level multiscale sensor array can be constituted. For any S r subarrays at level- r , r = 1 , , · · · , R , define G ( r ) ∈ C S r × L as the transmit subarray steering matrix. The overall transmitsubarray steering matrix is given by G = G ( R ) (cid:12) G ( R − (cid:12) · · · (cid:12) G (1) (cid:44) R (cid:12) r =1 G ( r ) . (44)Substituting (44) into (18), the tensor model for TB MIMOradar with an R -level miltiscale sensor array at transmit sideis given. Note that this is an ( R + 3) -order tensor. The reshapeof it and the use of Lemma 1 in this case are quite flexible.The DOA estimation can be conducted via
Algorithm 2 or Algorithm 3 analogously.Take a cubic transmit array, for example. It is a 3-level multiscale sensor array, and we can immediatelywrite the three transmit subarray steering matrices as Repeat the URA in Fig. 1 D times across the z -axis with coordinates (0 , , m d ) , d = 1 , · · · , D . G (1) = ∆ , G (2) = H and G (3) = Ξ , respec-tively, where Ξ (cid:44) [ τ , τ , · · · , τ D ] T V × L and τ d (cid:44) (cid:2) e − jπ ( m d − cosϕ , e − jπ ( m d − cosϕ , · · · , e − jπ ( m d − cosϕ L (cid:3) T .The parameter identifiability for different reshaped 3-ordertensors T ∈ C I × I × I varies, which is determined by theuniqueness condition of tensor decomposition, i.e., min(( I − I , I ) ≥ L (45)where I , I and I can be regarded as permutations of theset { S , S , · · · , S R , K, N, Q } , and I I I = KN Q R (cid:81) r =1 S r .From (45), it is possible to estimate the target DOAs usingonly one single pulse. The transmit array with multiple scalesof shift-invariances is exploited via the tensor reshape to makeup for the lack of number of snapshots. This property canalso be used to distinguish coherent sources. An example oftwo targets with identical Doppler shift has been discussed inSection III-A. See [37] for more discussions about the partialidentifiability of the tensor decomposition, where specificconditions for coherent or collocated sources are investigated.V. S IMULATION R ESULTS
In this section, we investigate the DOA estimation per-formance of the proposed method in terms of the rootmean square error (RMSE) and probability of resolution ofclosely spaced targets for TB MIMO radar. Throughout thesimulations, there are Q = 50 pulses in a single CPI.We assume that L = 3 targets lie within a given spatialsteering vector determined by { θ l } Ll =1 in linear array and { ( θ l , ϕ l ) } Ll =1 in planar array, the normalized Doppler shiftsare f = − . , f = 0 . and f = 0 . . The number ofMonte Carlo trials is P = 200 . The RCS of every targetis drawn from a standard Gaussian distribution, and obeysthe Swerling I model. Note that the last two targets shareidentical Doppler shift, which cause C to drop rank. Thenoise signals are assumed to be Gaussian, zero-mean and whiteboth temporally and spatially. The K orthogonal waveformsare S k ( t ) = (cid:113) T e j π kT t , k = 1 , · · · , K . For both linear andplanar array, the tensor model in (18) is used and the TBmatrix is pre-designed [6], [26].For linear array, we assume a transmit ULA with S = 8 subarrays. Each transmit subarray has M = 10 elementsspaced at half the wavelength. The placement of transmitsubarrays can vary from totally overlapped case to non-overlapped case. The number of transmit elements is computedby M = M + ∆ m ( S − . The receive array has N = 12 elements, which are randomly selected from the transmit array.For planar array, the reference transmit subarray is a × URA. The number of subarrays is S = 6 , where J = 2 and I = 3 . The distances between subarrays in both directions arefixed as the working wavelength, which means that ∆ m x = 2 and ∆ m y = 2 . A number of N = 12 elements in the transmitarray are randomly chosen as the receive array.For comparisons, ESPRIT-based algorithm [8] that exploitsthe phase rotations between transmit subarrays and U-ESPRITalgorithm [10] that utilizes the conjugate symmetric propertyof array manifold are used as signal covariance matrix-based DOA estimation methods, while conventional ALS algorithm[12], [16] that decomposes the factor matrices iteratively isutilized as signal tensor decomposition-based DOA estimationmethod. The Cramer-Rao lower bound (CRLB) for MIMOradar is also provided. For target DOAs estimated by the factormatrix X , if applicable, we use a postfix to distinguish it,e.g., ALS-sub (Proposed-sub) refers to the estimation resultcomputed by X , while ALS (Proposed) denotes the estimationresult originated from G after tensor decomposition. A. Example 1: RMSE and Probability of Resolution for LinearArray with Non-overlapped Subarrays
Three targets are placed at θ l = [ − ◦ , ◦ , ◦ ] . Considerthe matricized form of Z in (18). The goal is to estimate θ l from Z = T (3) + τ R , where T (3) is given by (25)and G = K , the SNR is measured as: SN R [ dB ] =10 log (cid:16)(cid:13)(cid:13) T (3) (cid:13)(cid:13) F (cid:46) (cid:107) τ R (cid:107) F (cid:17) . The RMSE is computed by RM SE = (cid:118)(cid:117)(cid:117)(cid:116) P L L (cid:88) l =1 P (cid:88) p =1 (cid:16) ˆ θ l ( p ) − θ l ( p ) (cid:17) . As shown in Fig. 2a, the RMSE results decline graduallywith the rise of SNR for all methods. The ESPRIT-basedalgorithm merely exploits the phase rotations between transmitsubarrays and therefore the performance is quite poor. U-ESPRIT algorithm performs better since the number of snap-shots is doubled. For conventional ALS algorithm and ourproposed method, target angular information can be obtainedfrom both factor matrices K and X , which are used to compareto each other to eliminate the potential grating lobes. Theproposed method approaches the CRLB with a lower thresholdas compared to the ALS method, since the Vandermondestructure of the factor matrix is exploited. Therefore, theproposed method performs better at low SNR. Note that thecomplexity of our proposed method is reduced significantly, itrequires approximately the same number of flops as comparedto that of the ALS method in a single iteration. Also, thecomparison of the estimation results between G and X showsa reasonable difference. This is mainly caused by the differentapertures of the subarray and the whole transmit array.For the probability of resolution, we assume only twoclosely spaced targets located at θ l = [ − ◦ , − ◦ ] . Thesetwo targets are considered to be resolved when (cid:13)(cid:13)(cid:13) ˆ θ l − θ l (cid:13)(cid:13)(cid:13) ≤(cid:107) θ − θ (cid:107) / , l = 1 , . The Doppler shifts are both f = 0 . and the other parameters are the same as before.In Fig. 2b, the probability of resolution results for allmethods tested are shown and they are consistent with thosein Fig. 2a. All methods achieve absolute resolution in highSNR region, and resolution declines with the decrease of SNR.The ESPRIT method presents the worst performance whileperformance of the U-ESPRIT improves slightly. The resultsof the ALS-sub method and the Proposed-sub method arealmost the same. A gap of approximately 3 dB SNR can beobserved between the proposed method and the ALS method,which means that our proposed method enables the lowestSNR threshold. The performance of both accuracy and reso-lution for our proposed method surpasses the other methods -30 -20 -10 0 10 20 30 SNR[dB] -4 -3 -2 -1 R M SE [ D e g ] ESPRITU-ESPRITALS-subALSProposed-subProposedCRLB (a) RMSE versus SNR -30 -20 -10 0 10 20 30
SNR[dB] P r ob a b ili t y o f R es o l u t i on ESPRITU-ESPRITALS-subALSProposed-subProposed (b) Resolution versus SNR -30 -20 -10 0 10 20 30
SNR[dB] -3 -2 -1 R M SE [ D e g ] ESPRITU-ESPRITALS-2ALSSVD-2SVDCRLB -30 -20 -10 0 10 20 30
SNR[dB] -3 -2 -1 R M SE [ D e g ] ESPRITU-ESPRITALS-2ALSSVD-2SVDCRLB -30 -20 -10 0 10 20 30
SNR[dB] -3 -2 -1 R M SE [ D e g ] ESPRITU-ESPRITALS-2ALSSVD-2SVDCRLB -30 -20 -10 0 10 20 30
SNR[dB] -3 -2 -1 R M SE [ D e g ] ESPRITU-ESPRITALS-subALSProposed-subProposedCRLB (c) RMSE versus SNR -30 -20 -10 0 10 20 30
SNR[dB] P r ob a b ili t y o f R es o l u t i on ESPRITU-ESPRITALS-subALSProposed-subProposed (d) Resolution versus SNRFig. 2. DOA estimation performance for TB MIMO radar with uniformlyspaced subarrays for linear array (a)-(b) and planar array (c)-(d), 200 trials. since the shift-invariance between and within different transmitsubarrays are fully exploited.
B. Example 2: RMSE and Probability of Resolution for PlanarArray with Partly Overlapped Subarrays
In this example, three targets are placed at θ l =[ − ◦ , − ◦ , − ◦ ] and ϕ l = [25 ◦ , ◦ , ◦ ] . The signalmodel Z = T (3) + τ R is applied, where T (3) is given by(25) with G = H (cid:12) ∆ . The SNR is measured in the sameway as that in linear array. The RMSE for planar array iscompute by RMSE = (cid:118)(cid:117)(cid:117)(cid:116) P L L (cid:88) l =1 P (cid:88) p =1 (cid:20)(cid:16) ˆ θ l ( p ) − θ l ( p ) (cid:17) + ( ˆ ϕ l ( p ) − ϕ l ( p )) (cid:21) . In Fig. 2c, the RMSEs of ESPRIT, U-ESPRIT, ALS andthe proposed method are given. The CRLB is also provided.The performance of the ESPRIT method and the U-ESPRITmethod are relatively poor. It is because of the ignorance ofthe received signal shift-invariance within a single subarray.It can be observed that the Proposed-sub and the ALS-subsuccessfully estimate the target DOAs via X in the case ofplanar array, which proves the validity of (35). The resultscan be used to mitigate the spatial ambiguity in the followingestimations. Like their counterparts in linear array, the RMSEsof the proposed method and the ALS method are almost thesame for above 0 dB SNR while the performance of theproposed method in low SNR is better. The ALS methodignores the Vandermonde structure during tensor decomposi-tion. Compared to (35), the DOA estimation result in G takesadvantage of a larger aperture and therefore achieves a betterRMSE performance.To evaluate the resolution performance, only two targetsare reserved and the spatial directions are θ l = [ − ◦ , − ◦ ] and ϕ l = [15 ◦ , ◦ ] . The resolution is considered successfulif (cid:13)(cid:13)(cid:13) ˆ θ l − θ l (cid:13)(cid:13)(cid:13) ≤ (cid:107) θ − θ (cid:107) / , (cid:107) ˆ ϕ l − ϕ l (cid:107) ≤ (cid:107) ϕ − ϕ (cid:107) / , l =1 , . The target Doppler shifts are the same, given as f = 0 . .The other parameters are unchanged.Fig. 2d shows the results for all methods with respect to theprobability of resolution. The proposed method achieves thelowest SNR threshold, which benefits from the fully exploita-tion of the shift-invariance and the Vandermonde structureduring tensor decomposition. Note that the convergence of theALS method is unstable and can be influenced by the tensorsize. It can be observed that the resolution performance of theALS method is deteriorated as compared to its counterpart inFig. 2b. This conclusion implies that the robustness of ourproposed method is better regarding 2-D DOA estimation,since no iterations are required. C. Example 3: RMSE Performance for Linear Array withDifferent ∆ m In this example, we mainly consider the RMSE performancewhen ∆ m changes from one to at most M . The aperture isincreased gradually. The SNR is assumed to be 10 dB. Allother parameters are the same as those in Example 1.Given the number of subarrays and the structure of a singlesubarray, the aperture of the overall transmit ULA rises withthe increase of ∆ m while the number of elements shared bytwo adjacent subarrays declines. When ∆ m = 0 , this modelis identical to that for conventional ESPRIT method [6], [26],and there is no transmit subarray. When ∆ m rises, the distancebetween phase centers for two adjacent subarrays becomeslarger than half the working wavelength and grating lobes aregenerated. The locations of these grating lobes are determinedby (22), and can be eliminated. Meanwhile, the transmit arrayaperture is increased and the DOA estimation performanceshould be improved.To investigate the improvement, the RMSEs of three targetsare computed versus the rise of ∆ m . It can be seen in Fig. 3athat the RMSE results decrease steadily with the increaseof ∆ m . The ESPRIT method and U-ESPRIT method sufferfrom grating lobes and the received signal within a singlesubarray is not fully exploited, hence, they perform poorly.The RMSEs of the Proposed-sub and the ALS-sub are almostunchanged since the estimation is only based on the subarray,which is fixed during the simulation. Meanwhile, the proposedmethod and the ALS method achieve better accuracy thantheir counterparts originated from X when ∆ m > . It canbe noted in Fig. 2a that the convergence is satisfied for theALS method and our proposed method when SNR is above10 dB. Consequently, the RMSEs of the proposed method andthe ALS method are nearly coincident.To evaluate the RMSE performance versus ∆ m x or ∆ m y fora planar array, it is necessary to separately add a new subarrayin one direction while keeping the array structure in the otherdirection unchanged. This can be fulfilled by constructing anL-shaped transmit array, where each element is replaced by aURA subarray. However, this analysis would be beyond thescope of this paper. In general, it can be concluded that theproposed method can estimate the target DOAs via the phase m -3 -2 -1 R M SE [ D e g ] ESPRITU-ESPRITALS-subALSProposed-subProposedCRLB (a) RMSE versus ∆ m -30 -20 -10 0 10 20 30 SNR[dB] -3 -2 -1 R M SE [ D e g ] G-ESPRITALSProposedUniform(S = 7)Uniform(S = 9)CRLB(S = 7)CRLB(Non-uniform)CRLB(S = 9) (b) Generalized Vandermonde matrix -30 -20 -10 0 10 20 30
SNR[dB] -2 -1 E l eva t i on R M SE [ D e g ] ESPRITU-ESPRITInter-TEVProposedCRLB (c) Elevation RMSE versus SNR -30 -20 -10 0 10 20 30
SNR[dB] -2 -1 A z i m u t h R M SE [ D e g ] ESPRITU-ESPRITInter-TEVProposedCRLB (d) Azimuth RMSE versus SNRFig. 3. RMSE results for TB MIMO radar with different subarray configu-rations. rotations between transmit subarrays. If the placement of twoadjacent subarrays satisfies some conditions, e.g., ∆ m > for linear array, the RMSE performance is better than thatcomputed by a single subarray. Note that the received signal oftwo adjacent subarrays can be obtained by spatial smoothing[31], a proper spatial smoothing of the received signal canimprove the DOA estimation performance. D. Example 4: Generalized Vandermonde Factor Matrix forLinear Array with m s = { , , , , , } Here we evaluate the proposed DOA estimation methodfor TB MIMO radar with non-uniformly spaced transmitsubarrays. The transmit linear array has S = 7 subarrays with m s = { , , , , , , } . Each subarray contains M = 10 elements. The N = 12 elements are randomly chosen fromthe transmit array to form the receive array. Three targetsare placed at θ l = [ − ◦ , ◦ , ◦ ] with normalized Dopplershifts f l = [0 . , − . , − . . To simplify the signal model,each subarray is a ULA, which is not used during the DOAestimation in this example. Equations (38)-(43) can be applieddirectly, since the subarray structure stays identical. Twodifferent transmit arrays are introduced for comparison toillustrate the improved performance provided by constructing K ( sub ) . Both of them can be regarded as a linear array withuniformly spaced subarrays ( ∆ m = 1 ). The first one has S = 7 subarrays, while the second one has S = 9 subarrays toachieve the same aperture. The DOA estimation for these twotransmit arrays can be conducted by Algorithm 1 . Meanwhile,conventional ALS method can be applied to decompose thefactor matrix K ( sub ) , which will be used to estimate the targetDOAs by solving (43). The generalized-ESPRIT (G-ESPRIT)method in [36] is also used for comparison. The CRLBs ofthree different transmit arrays are also shown. From Fig 3b, it can be observed that the formulation of K ( sub ) exploits the multiple scales of shift-invariances ingeneralized Vandermonde matrix. By solving (43), the gratinglobes are eliminated efficiently. Hence, the structurs of transmitsubarrays can be arbitrary but identical, which provide moreflexibility for array design. The RMSE of the proposed methodsurpasses those of G-ESPRIT and ALS methods. Also, theperformance of the non-uniformly spaced transmit subarraysis better than that of the uniformly spaced transmit subarrays(S = 7). This is expected since the aperture is increased dueto sparsity. Compared to the fully spaced transmit subarraycase (S = 9), the performance of the proposed method isdeteriorated slightly. However, the fully spaced array can beextremely high-cost if the array aperture is further increased.By using the generalized Vandermonde matrix, the proposedmethod enables the sparsity in transmit array, which achieveshigher resolution with less elements. E. Example 5: Multiscale Sensor Array with Arbitrary butIdentical Subarrays
In the final example, we illustrate the performance of theproposed DOA estimation method for TB MIMO radar witharbitrary but identical subarrays. Specifically, a planar arraywith S = 4 × subarrays is considered, whose phase centersform a uniform rectangular grid with a distance of half theworking wavelength. For each subarray, M = 4 elementsare randomly placed in a circle centered on the phase centerwith a radius of a quarter of wavelength. The structure ofall subarrays are identical, hence, the transmit array can beregarded as an 3-level multiscale sensor array and the transmitsubarray steering matrix can be obtained from (44). The N =12 receive elements are randomly selected from the transmitarray. Three targets are placed at θ l = [ − ◦ , − ◦ , − ◦ ] and ϕ l = [11 ◦ , ◦ , ◦ ] . The other parameters are the sameas those in Example 2.Note that the subarray is arbitrary, the DOA information canonly be estimated by the phase rotations between the transmitsubarrays. Alternatively, the transmit array interpolation tech-nique [13] is introduced to map the original transmit arrayinto a × URA to enable the ESPRIT-like DOA estimation,which is referred to as Inter-TEV in Fig. 3c and Fig. 3d. It canbe observed that by carefully designing the mapping matrix,the RMSEs of the Inter-TEV method are better than thoseof ESPRIT and U-ESPRIT methods for both elevation andazimuth estimation. However, the proposed method surpassesthe other methods with a lower RMSE. This is because of thefull usage of the shift-invariance between and within differenttransmit subarrays. VI. C
ONCLUSION
The problem of tensor decomposition with Vandermondefactor matrix in application to DOA estimation for TB MIMOradar with arbitrary but identical transmit subarrays has beenconsidered. A general 4-order tensor that can be used toexpress the TB MIMO radar received signal in a variety ofscenarios, e.g., linear and planar arrays, uniformly and non-uniformly spaced subarrays, regular and irregular subarrays, has been designed. The shift-invariance of the received signalbetween and within different transmit subarrays have beenused to conduct DOA estimation. Specifically, a computation-ally efficient tensor decomposition method has been proposedto estimate the generators of the Vandermonde factor matrices,which can be used as a look-up table for finding target DOA.The proposed method fully exploits the shift-invariance ofthe received signal between and within different subarrays,which can be regarded as a generalized ESPRIT method.Comparing with conventional signal tensor decomposition-based techniques, our proposed method take advantage ofthe Vandermonde structure of factor matrices, and it requiresno iterations or any prior information about the tensor rank.The parameter identifiability of our tensor model has alsobeen studied via the discussion of the uniqueness conditionof tensor decomposition. Simulation results have verified thatthe proposed DOA estimation method has better accuracy andhigher resolution as compared to existing techniques.A PPENDIX AP ROOF OF L EMMA First, let A (1) ∈ C I × L be a Vandermonde matrix withdistinct generators, we have r (cid:0) A (1) (cid:12) A (2) (cid:1) = min( I I , L ) ,since it is the KR product of a Vandermonde matrix andan arbitrary matrix [17]. Assuming I ≥ L, I I ≥ L , thefollowing results r (cid:0) A (1) (cid:12) A (2) (cid:1) = L and r ( A (3) ) = L hold,since A (3) has full column rank. The proof of Lemma 1 inthis case is identical to that of Proposition III.2 in [17].Next, let A (1) = B (cid:12) C , where B ∈ C J × L and C ∈ C I × L are both Vandermonde matrices with distinctgenerators. Consider the rank of matrix A (1) (cid:12) A (2) = B (cid:12) C (cid:12) A (2) = Π (cid:0) B (cid:12) A (2) (cid:12) C (cid:1) , where Π is an exchangematrix. Again, the rank of B (cid:12) A (2) is min( JI , L ) while r (cid:0) B (cid:12) A (2) (cid:12) C (cid:1) = min( IJI , L ) . Since Π is nonsingular, r (cid:0) A (1) (cid:12) A (2) (cid:1) = L .The mode-3 unfolding of Y is Y (3) = (cid:0) A (1) (cid:12) A (2) (cid:1) A (3) T . The SVD of this matrix representationis denoted by Y (3) = UΛV H , where U ∈ C I I × L , Λ ∈ C L × L , and V ∈ C I × L . Since r (cid:0) A (1) (cid:12) A (2) (cid:1) = L and r (cid:0) A (3) (cid:1) = L , it can be derived that a nonsingular matrix E ∈ C L × L satisfies UE = A (1) (cid:12) A (2) . (46)The Vandermonde structure of both B and C can beexploited via U E = B (cid:12) C (cid:12) A (2) = (cid:16) B (cid:12) C (cid:12) A (2) (cid:17) Ω b = U EΩ b U E = B (cid:12) C (cid:12) A (2) = (cid:16) B (cid:12) C (cid:12) A (2) (cid:17) Ω c = U EΩ c (47)where Ω b = diag ( ω b ) , Ω c = diag ( ω c ) with ω b and ω c denoting the vectors of generators of B and C , respectively.The submatrices U , U , U and U are truncated from rowsof U according to the operator of the KR product, i.e., U = (cid:2) I IK ( J − , IK ( J − × IK (cid:3) UU = (cid:2) IK ( J − × IK , I IK ( J − (cid:3) UU = (cid:0) I J ⊗ (cid:2) I K ( I − , K ( I − × K (cid:3)(cid:1) UU = (cid:0) I J ⊗ (cid:2) K ( I − × K , I K ( I − (cid:3)(cid:1) U . (48) Note that E , Ω b and Ω c are full rank. We have U † U = EΩ b E − and U † U = EΩ c E − . Hence, the vectors ω b and ω c can be computed as the collections of eigenvalues of U † U and U † U , respectively, while E is the matrix of collectionof the corresponding eigenvectors. From the generators of B and C , the first factor matrix A (1) can be reconstructed.Meanwhile, it can be observed that (cid:32) α (1) Hl α (1) Hl α (1) l ⊗ I I (cid:33) (cid:16) α (1) l ⊗ α (2) l (cid:17) = α (2) l . (49)Assume that the column vectors of A (1) have unit form,then α (2) l can be written as α (2) l = ( α (1) Hl ⊗ I I ) Ue l , l = 1 , , · · · , L (50)where e l is the l -th column of E . Given A (1) and A (2) , thethird factor matrix can be computed by A (3) T = (cid:16)(cid:16) A (1) H A (1) (cid:17) ∗ (cid:16) A (2) H A (2) (cid:17)(cid:17) − × (cid:16) A (1) (cid:12) A (2) (cid:17) H Y (3) . (51)Therefore, the tensor decomposition of Y is genericallyunique, where A (1) can be a Vandermonde matrix or the KRproduct of two Vandermonde matrices with distinct generators,and A (3) is column full rank.A PPENDIX BP ROOF OF (13)Given y ( q ) s in (12), s = 1 , , · · · , I, I + 1 , · · · , IJ , concate-nate every I vectors together to form totally J matrices ofidentical dimension KN × I . These matrices are denoted by ¯Y ( q ) j , and are given as ¯Y ( q ) j = (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3) (cid:0) c Tq (cid:12) h Tj (cid:12) ∆ (cid:1) T + ¯N ( q ) j (52)where ¯N ( q ) j (cid:44) (cid:104) n ( q )( j − I +1 , n ( q )( j − I +2 , · · · , n ( q )( jI ) (cid:105) . The noise-free version of (52) can be rewritten as ¯Y ( q ) j = (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3) Γ q (cid:0) h Tj (cid:12) ∆ (cid:1) T (53)where Γ q = diag ( c q ) . Since (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3) Γ q is fixed, theconcatenation of J matrices merely depends on the concate-nation of (cid:0) h Tj (cid:12) ∆ (cid:1) T . Define a matrix ΘΘ (cid:44) h T (cid:12) ∆h T (cid:12) ∆ ... h TJ (cid:12) ∆ S × L (54)such that Θ performs the concatenation. From the definitionof the KR product, it can be observed that Θ = H (cid:12) ∆ .Therefore, the concatenation of ¯Y ( q ) j is given by ¯Y ( q ) = (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3) Γ q Θ T = (cid:2)(cid:0) W H A (cid:1) (cid:12) B (cid:3) (cid:0) c Tq (cid:12) H (cid:12) ∆ (cid:1) T . (55)Considering the noise term, equation (13) can be given. R EFERENCES[1] A. M. Haimovich, R. S. Blum, and L. J. Cimini, “MIMO radar withwidely separated antennas,”
IEEE Signal Process. Mag. , vol. 25, no. 1,pp. 116–129, Jan. 2008.[2] J. Li and P. Stoica, “MIMO radar with colocated antennas,”
IEEE SignalProcess. Mag. , vol. 24, no. 5, pp. 106–114, Sep. 2007.[3] A. Hassanien and S. A. Vorobyov, “Transmit energy focusing for DOAestimation in MIMO radar with colocated antennas,”
IEEE Trans. SignalProcess. , vol. 59, no. 6, pp. 2669–2682, Jun. 2011.[4] A. Hassanien and S. A. Vorobyov, “Phased-MIMO radar: A tradeoffbetween phased-array and MIMO radars,”
IEEE Trans. Signal Process. ,vol. 58, no. 6, pp. 3137–3151, Jun. 2010.[5] D. R. Fuhrmann, J. P. Browning, and M. Rangaswamy, “Signalingstrategies for the hybrid MIMO phased-array radar,”
IEEE J. Sel. TopicsSignal Process. , vol. 4, no. 1, pp. 66–78, Feb. 2010.[6] A. Khabbazibasmenj, A. Hassanien, S. A. Vorobyov, and M. W.Morency, “Efficient transmit beamspace design for search-free basedDOA estimation in MIMO radar,”
IEEE Trans. Signal Process. , vol. 62,no. 6, pp. 1490–1500, Mar. 2014.[7] Z. Guo, X. Wang, and W. Heng, “Millimeter-wave channel estimationbased on 2-D beamspace MUSIC method,”
IEEE Trans. Wireless Com-mun. , vol. 16, no. 8, pp. 5384–5394, Aug. 2017.[8] A. Hu, T. Lv, H. Gao, Z. Zhang, and S. Yang, “An ESPRIT-basedapproach for 2-D localization of incoherently distributed sources inmassive MIMO systems,”
IEEE J. Sel. Topics Signal Process. , vol. 8,no. 5, pp. 996–1011, Oct. 2014.[9] N. Tayem and H. M. Kwon, “L-shape 2-dimensional arrival angleestimation with propagator method,”
IEEE Trans. Antennas Propag. ,vol. 53, no. 5, pp. 1622–1630, May 2005.[10] M. D. Zoltowski, M. Haardt, and C. P. Mathews, “Closed-form 2-Dangle estimation with rectangular arrays in element space or beamspacevia unitary ESPRIT,”
IEEE Trans. Signal Process. , vol. 44, no. 2, pp.316–328, Feb. 1996.[11] B. Xu, Y. Zhao, Z. Cheng, and H. Li, “A novel unitary PARAFACmethod for DOD and DOA estimation in bistatic MIMO radar,”
SignalProcessing , vol. 138, pp. 273 – 279, Sep. 2017.[12] D. Nion and N. D. Sidiropoulos, “Tensor algebra and multidimensionalharmonic retrieval in signal processing for MIMO radar,”
IEEE Trans.Signal Process. , vol. 58, no. 11, pp. 5693–5705, Nov. 2010.[13] M. Cao, S. A. Vorobyov, and A. Hassanien, “Transmit array interpolationfor DOA estimation via tensor decomposition in 2-D MIMO radar,”
IEEE Trans. Signal Process. , vol. 65, no. 19, pp. 5225–5239, Oct. 2017.[14] S. D. Blunt and E. L. Mokole, “Overview of radar waveform diversity,”
IEEE Trans. Aerosp. Electron. Syst. , vol. 31, no. 11, pp. 2–42, Nov.2016.[15] A. Hassanien, M. W. Morency, A. Khabbazibasmenj, S. A. Vorobyov,J. Park, and S. Kim, “Two-dimensional transmit beamforming for MIMOradar with sparse symmetric arrays,” in
Proc. IEEE Radar Conf. , Ottawa,ON, Canada, Apr. 2013, pp. 1–6.[16] N. D. Sidiropoulos, L. De Lathauwer et al. , “Tensor decomposition forsignal processing and machine learning,”
IEEE Trans. Signal Process. ,vol. 65, no. 13, pp. 3551–3582, Jul. 2017.[17] M. Sørensen and L. De Lathauwer, “Blind signal separation via tensordecomposition with vandermonde factor: Canonical polyadic decompo-sition,”
IEEE Trans. Signal Process. , vol. 61, no. 22, pp. 5507–5519,Nov. 2013.[18] T. Jiang, N. D. Sidiropoulos, and J. M. F. ten Berge, “Almost-sureidentifiability of multidimensional harmonic retrieval,”
IEEE Trans.Signal Process. , vol. 49, no. 9, pp. 1849–1859, Sep 2001.[23] J. H. de M. Goulart, M. Boizard, R. Boyer, G. Favier, and P. Comon,“Tensor CP decomposition with structured factor matrices: Algorithmsand performance,”
IEEE J. Sel. Topics Signal Process. , vol. 10, no. 4,pp. 757–769, Jun. 2016. [19] F. Xu, S. A. Vorobyov, and X. Yang, “Joint DOD and DOA estimationin slow-time MIMO radar via PARAFAC decomposition,”
IEEE SignalProcess. Lett. , vol. 27, pp. 1495–1499, Aug. 2020.[20] A. Cichocki, D. Mandic, L. De Lathauwer, G. Zhou, Q. Zhao, C. Caiafa,and H. A. PHAN, “Tensor decompositions for signal processing appli-cations: From two-way to multiway component analysis,”
IEEE SignalProcess. Mag. , vol. 32, no. 2, pp. 145–163, Mar. 2015.[21] F. Xu, X. Yang, and T. Lan, “Search-free DOA estimation method basedon tensor decomposition and polynomial rooting for transmit beamspaceMIMO radar,” arXiv: 2010.03296 , Oct. 2020.[22] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinearsingular value decomposition,”
SIAM J. Matrix Anal. Appl. , vol. 21,no. 4, pp. 1253–1278, 2000.[24] W. Wang, H. C. So, and A. Farina, “An overview on time/frequencymodulated array processing,”
IEEE J. Sel. Topics Signal Process. ,vol. 11, no. 2, pp. 228–246, Mar. 2017.[25] L. Lu, G. Y. Li, A. L. Swindlehurst, A. Ashikhmin, and R. Zhang,“An overview of massive MIMO: Benefits and challenges,”
IEEE J. Sel.Topics Signal Process. , vol. 8, no. 5, pp. 742–758, Oct 2014.[26] M. W. Morency and S. A. Vorobyov, “Partially adaptive transmitbeamforming for search free 2D DOA estimation in MIMO radar,” in
Proc. 23rd Eur. Signal Process. Conf. , Nice, France, Aug. 2015, pp.2631–2635.[27] F. Xu and S. A. Vorobyov, “Constrained tensor decomposition for 2DDOA estimation in transmit beamspace MIMO radar with subarrays,”in
Submit to Proc. 46th Int. Conf. Acoust., Speech, Signal Process(ICASSP), 2021 , Jun. 2021.[28] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,”
SIAM Rev. , vol. 51, no. 3, pp. 455–500, 2009.[29] R. Roy and T. Kailath, “ESPRIT-estimation of signal parameters viarotational invariance techniques,”
IEEE Trans. Acoust., Speech, SignalProcess. , vol. 37, no. 7, pp. 984–995, Jul. 1989.[30] A. Phan, P. Tichavsk´y, and A. Cichocki, “CANDECOMP/PARAFACdecomposition of high-order tensors through tensor reshaping,”
IEEETrans. Signal Process. , vol. 61, no. 19, pp. 4847–4860, Oct. 2013.[31] M. Sørensen and L. De Lathauwer, “Multiple invariance ESPRIT fornonuniform linear arrays: A coupled canonical polyadic decompositionapproach,”
IEEE Trans. Signal Process. , vol. 64, no. 14, pp. 3693–3704,Jul. 2016.[32] A. L. Swindlehurst, B. Ottersten, R. Roy, and T. Kailath, “Multipleinvariance esprit,”
IEEE Trans. Signal Process. , vol. 40, no. 4, pp. 867–881, Apr. 1992.[33] K. V. Mishra, M. R. Bhavani Shankar, V. Koivunen, B. Ottersten, andS. A. Vorobyov, “Toward millimeter-wave joint radar communications:A signal processing perspective,”
IEEE Signal Process. Mag. , vol. 36,no. 5, pp. 100–114, Sep. 2019.[34] S. Miron, Y. Song, D. Brie, and K. T. Wong, “Multilinear directionfinding for sensor-array with multiple scales of invariance,”
IEEE Trans.Aerosp. Electron. Syst. , vol. 51, no. 3, pp. 2057–2070, Jul. 2015.[35] M. D. Zoltowski and K. T. Wong, “Closed-form eigenstructure-baseddirection finding using arbitrary but identical subarrays on a sparseuniform cartesian array grid,”
IEEE Trans. Signal Process. , vol. 48,no. 8, pp. 2205–2210, Aug. 2000.[36] B. Liao and S. Chan, “Direction-of-arrival estimation in subarrays-basedlinear sparse arrays with gain/phase uncertainties,”
IEEE Trans. Aerosp.Electron. Syst. , vol. 49, no. 4, pp. 2268–2280, Oct. 2013.[37] X. Guo, S. Miron, D. Brie, and A. Stegeman, “Uni-mode and partialuniqueness conditions for CANDECOMP/PARAFAC of three-way ar-rays with linearly dependent loadings,”