# A Propagator Method for Bistatic Coprime EMVS-MIMO Radar

Authors: Yongqiang Yang; Ningjun Ruan; Guanjun Huang; Junpeng Shi; Fangqing Wen
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, a novel two-dimensional (2D) direction-of-departure (DOD) and 2D direction-of-arrival (DOA) estimate algorithm is proposed for bistatic multiple-input multiple-output (MIMO) radar system equipped with coprime electromagnetic vector sensors (EMVS) arrays. Firstly, we construct the propagator to obtain the signal subspace. Then, the ambiguous angles are estimated by using rotation invariant technique. Based on the characteristic of coprime array, the unambiguous angles estimation is achieved. Finally, all azimuth angles estimation is followed via vector cross product. Compared to the existing uniform linear array, coprime MIMO radar is occupying large array aperture, and the proposed algorithm does not need to obtain signal subspace by eigendecomposition. In contrast to the state-of-the-art algorithms, the proposed algorithm shows better estimation performance and simpler computation performance. The proposed algorithm’s effectiveness is proved by simulation results.


---

## Body

## 1. Introduction

In recent ten years, multiple-input multiple-output (MIMO) radar has aroused extensive attention. MIMO radar has many important branches, for example, detection, parameter estimation, waveform design, and synchronization. Among them, the estimation of DOD and DOA is one of the important research fields of bistatic MIMO radar. Up to now, many algorithms have been proposed to solve the above issue, for example, signal parameters estimation method based on rotational invariance technique (ESPRIT) [1], maximum-likelihood estimator [2], spectrum peak search method [3], tensor approach [4, 5], sparsity-aware strategy [6], and propagator method (PM) [7]. However, most of the existing algorithms can only solve the one-dimensional angle estimation problem. There are just a small number of algorithms focusing on the 2D-DOD and 2D-DOA estimation in MIMO radar [8–10]. For the 2D estimation algorithms in [8–10], the key of them is to utilize nonlinear geometry of the scalar sensor array [8–10]. On the other hand, some works propose the use of electromagnetic vector sensor (EMVS) to solve the 2D-DOD and 2D-DOA estimation in bistatic MIMO radar. Different from the scalar sensor array, a signal EMVS can not only perform the task of two-dimensional direction angle estimation but also provide the polarization status of the incoming signal.The EMVS is firstly proposed to be applied to bistatic MIMO radar in [11], where an EMVS array is mounted at transmitting and an EMVS array is at receiving. In [12], a new method for 2D-DOD and 2D-DOA estimation utilizing transmitting and receiving EMVS arrays is proposed. The framework used the ESPRIT-based method to estimate the 2D-DOD and 2D-DOA. Then, vector cross-product strategy is used to estimate the polarization status and azimuth angles. In [13], a method similar to PM is proposed to avoid the step of eigendecomposition in [12]. However, none of the above methods can make full use of the MIMO radar’s virtual aperture when we apply the selection matrix to the array measurement. In addition, the mentioned methods in [12–14] require extensive pair calculation and offer limited identifiability. For avoiding existing shortcomings in above algorithms, a parallel factor (PARAFAC) scheme was proposed in [15]. However, low computational efficiency and sensitivity to initialization are the drawbacks of PARAFAC decomposition. In [16], an improved PM-like method was introduced, which can realize the rotation invariance of the whole virtual steering vector. In addition, it can automatically pair the elevation angles associated with the receiving and transmitting arrays. More recently, the ESPRIT-like approach was investigated in [17], which is suitable for distributing EMVS sensor.In order to eliminate the effect of phase ambiguity, less than half wavelength is the limited requirement of distance between sensors, which limits the receive and transmit array’s virtual aperture. In addition, the mutual coupling problem may occur when sensors are closely displaced [18], which may reduce the estimation performance. The coprime array has a promising prospect to settle the above problems, so it has attracted extensive attention [19–21]. Also, the EMVS with nested array geometry has also been considered in [22, 23]. Two sparse and uniform linear arrays constitute a coprime linear array (ULA), in which the internal element spacings are coprime numbers. It is ambiguous to estimate the angle of two subarrays separately, but the coprime property can uniquely determine it. Due to the increase of array aperture, coprime array performs better than ULA in angle estimation. The results show that the coprime EMVS array can effectively improve the performance of parameter estimation [24]. However, there is little research on EMVS-MIMO radar’s parameter estimation.This paper presents a bistatic coprime EMVS-MIMO radar framework different from ULA manifold [12–16], since both receivers and transmitters are composed of coprime EMVS. On this basis, the method of using propagator to construct signal subspace is proposed. Then the transmitting elevation and receiving elevation are estimated and paired one-to-one, but the estimation result is ambiguous. Fortunately, the coprime characteristic allows a unique elevation. Finally, the cross-product technique is used to estimate the azimuth angle, and the results are automatically paired to estimate the elevation angle. Based on coprime array, this method increases the aperture, so it has better estimation performance. Numerical examples show that the estimation effect is improved.Throughout the paper, lowercase letters represent vectors and uppercase letters denote matrices, respectively. AnM×M identity matrix is represented by IM, 0M stands for the M×M full-zeros matrix, and 1n represents a row vector whose n−th entity is one and others are zero; ⊗ represents the Kronecker product and ⊙ is the Khatri-Rao product. The operations of transpose are indicated by superscript AT, and the superscript AH represents Hermitian transpose; meanwhile, A† and A−1 denote pseudo-inverse and inverse, respectively; Bp:q,: returns the rows from p−th to q−th of B; BnB denotes a diagonal matrix whose diagonal entities consist of the n−th row of B; ⋅F represents the Frobenius norm; angle⋅ means to gain the phase. real⋅ returns the real part in the bracket, and ⊕ is Hadamard product. The vector cross product between two column vectors e1=e1,e2,e3T and e2=e4,e5,e6T is defined as e1∗e2=0−e3e2e30−e1−e2e10e4e5e6.


## 2. The Proposed Algorithm

### 2.1. Problem Formulation

As shown in Figure1, suppose that there is a bistatic EMVS-MIMO radar system, whose receiving and transmitting are distributed in the coprime geometries. In terms of the transmit array, let us assume that Subarray 1 is equipped with Mt EMVS and Subarray 2 is equipped with Nt EMVS, where Mt and Nt are coprime integers. Subarray 1’s adjacent distance is Ntλ/2 and Subarray 2’s adjacent distance is Mtλ/2, where λ represents the carrier wavelength. Similar to transmitting array, two EMVS subarrays constitute the receiving array, which have Mr and Nr elements, respectively. What is more, it is assumed that there are K far-field targets appearing in same range bin, and k−th target’s 2D-DOA pair is θr,k,ϕr,k, and 2D-DOD pair is θt,k,ϕt,k, where the elevation angles are denoted by θt,k and θr,k, and the azimuth angles are denoted by ϕt,k and ϕr,k. The matched outputs can be expressed as(1)yτ=∑k=1Kat,k⊗bt,k⊗ar,k⊗br,kskτ+nτ=At⊙Bt⊙Ar⊙Brsτ+nτ,where τ represents the slow-time index (pulse index), and sτ=s1τ,s2τ,…,sKτT represents the target reflect coefficient vector; at,k=ejNtMt−1πsinθt,k,…,ejNtπsinθt,k,1,e−Mtπsinθt,k,…,e−jMtNt−1πsinθt,kT is the k−th transmit steering vector, and ar,k=ejNrMr−1πsinθr,k,…,ejNrπsinθr,k,1,e−Mrπsinθr,k,…,e−jMrNr−1πsinθr,kT is the k−th receive steering vector; bt,k∈ℂ6×K represents the k−th polarization response vectors of transmit array, and br,k∈ℂ6×K represents the k−th polarization response vectors of the receive array; At=at,1,at,2,…,at,K∈ℂM×K, Ar=ar,1,ar,2,…,ar,K∈ℂN×K, where M=Mt+Nt−1 and N=Mr+Nr−1; Br=br,1,br,2,…,br,K∈ℂ6×K and Bt=bt,1,bt,2,…,bt,K∈ℂ6×K stand for the associate polarization response matrices, respectively. The Gaussian noise vector with zero mean and variance of σ2 is represented by nτ. Moreover,(2a)bt,k=Dt,kvt,k,(2b)br,k=Dr,kvr,k,with(3a)Dt,k=cosϕt,kcosθt,k−sinϕt,ksinϕt,kcosθt,kcosϕt,k−sinθt,k0−sinϕt,k−cosϕt,kcosθt,kcosϕt,k−sinϕt,kcosθt,k0sinθt,k,(3b)Dr,k=cosϕr,kcosθr,k−sinϕr,ksinϕr,kcosθt/r,kcosϕr,k−sinθr,k0−sinϕr,k−cosϕr,kcosθr,kcosϕr,k−sinϕr,kcosθr,k0sinθr,k,(4a)vt,k=sinγt,kejηt,kcosγt,k,(4b)vr,k=sinγr,kejηr,kcosγr,k,where γt,k and γr,k represent the k−th transmit and receive auxiliary polarization angles, respectively. The corresponding polarization phase differences are indicated, respectively, by ηt,k and ηr,k. Let btp,k and brq,k be the p-th and q-th (p,q∈1,2,…,6) entities of bt,k and br,k; that is, bt,k=bt1,k,bt2,k,…,bt6,kT and br,k=br1,k,br2,k,…,br6,kT. Define the following vectors:(5a)bt1,k=bt1,k,bt2,k,bt3,kT,(5b)bt2,k=bt4,k,bt5,k,bt6,kT,(5c)br1,k=br1,k,br2,k,br3,kT,(5d)br2,k=br4,k,br5,k,br6,kT.Figure 1
Bistatic coprime EMVS-MIMO radar model.It should be emphasized thatbt1,k and br1,k represent the electric field information vectors, while bt2,k and br2,k account for the magnetic field information vectors. In addition,(6a)bt1,k∗bt2,k=sinθt,kcosϕt,ksinθt,ksinϕt,kcosθt,k,(6b)br1,k∗br2,k=sinθr,kcosϕr,ksinθr,ksinϕr,kcosθr,k.When the targets are uncorrelated, theyτ’s covariance matrix is(7)R=CRsCH+σ2I36MN,where C=Ct⊙Cr, Ct=At⊙Bt, and Cr=Ar⊙Br; Rs=diagr1,r2,…,rK is the sτ’s covariance matrix, and rk represents the power of skτ. In practical application, the estimation of R can be obtained via L available snapshots; that is,(8)R^=1L∑τ=1LyτyHτ.Herein, estimating the 2D-DOD and 2D-DOA fromR^ is our goal.


### 2.2. Propagator

Let the firstK rows of C be C1∈ℂK×K, and let the rest be C2∈ℂ36MN−K×K; that is,(9)C=C1C2.SinceC1 is a nonsingular matrix, taking a unique linear transformation of C1 can get C2; that is,(10)C2=PcC1,where Pc∈ℂ36MN−K×K represents the propagator. Let P=IKPc; we can get(11)PC1=C1C2=C.Denote the noiseless item ofR as R¯; that is, R¯=CRsCH. Divide R¯ into two submatrices as(12)R¯=G,H,where G∈ℂ36MN×K and H∈ℂ36MN×36MN−K. Obviously,(13a)G=CRsC1H,(13b)H=CRsC2H.According to the relationship in (10), we have(14)GPcH=H.Consequently,Pc can be estimated by fitting(15)minH^−G^PcHF,where G^ is the estimation of G and H^ is the estimation of H. According to equation (15), we can obtain the Pc’s least-squares solution by(16)P^cH=G^†H^.Thereafter, we can obtain theP’s estimation via P^=IKP^c. Based on equation (11), there is a full-rank matrix T so that(17)P^=CT.Obviously,P^ plays a similar role as the signal subspace that is obtained from eigendecomposition.


### 2.3. Ambiguous Elevation Angle Estimation

In what follows, we use the rotational invariance property to obtain the estimation of the elevation anglesθt,k and θr,k. Let us define Jt,1=IMt−1,0Mt−1×Nt∈ℂMt−1×M, Jt,2=0Mt−1×1,IMt−1,0Mt−1×Nt−1∈ℂMt−1×M, Jt,3=0Nt−1×Mt−1,INt−1,0Nt−1×1∈ℂNt−1×M, Jt,4=0Nt−1×Mt,INt−1∈ℂNt−1×M, Jr,1=IMr−1,0Mr−1×Nr∈ℂMr−1×N, Jr,2=0Mr−1×1,IMr−1,0Mr−1×Nr−1∈ℂMr−1×N, Jr,3=0Nr−1×Mr−1,INr−1,0Nr−1×1∈ℂNr−1×N, and Jr,4=0Nr−1×Mr,INr−1∈ℂNr−1×N. It is easy to draw a conclusion that(18a)Jt,1At=Jt,2AtΘt1,(18b)Jt,3At=Jt,4AtΘt2,(18c)Jr,1Ar=Jr,2ArΘr1,(18d)Jr,3Ar=Jr,4ArΘr2,where Θt1=diagejNtπsinθt,1,ejNtπsinθt,2,…,ejNtπsinθt,K, Θt2=diagejMtπsinθt,1,ejMtπsinθt,2,…,ejMtπsinθt,K, Θr1=diagejNrπsinθr,1,ejNrπsinθr,2,…,ejNrπsinθr,K, and Θr2=diagejMrπsinθr,1,ejMrπsinθr,2,…,ejMrπsinθr,K. Meanwhile, we define(19a)JN1=Jt,1⊗I36N,(19b)JN2=Jt,2⊗I36N,(19c)JN3=Jt,3⊗I36N,(19d)JN4=Jt,4⊗I36N,(19e)JM1=I6M⊗Jr,1⊗I6,(19f)JM2=I6M⊗Jr,2⊗I6,(19g)JM3=I6M⊗Jr,3⊗I6,(19h)JM4=I6M⊗Jr,4⊗I6.Because there is a property ofA⊗BC⊗D=AC⊗BD, we have(20a)JN1C=JN2CΘt1,(20b)JN3C=JN4CΘt2,(20c)JM1C=JM2CΘr1,(20d)JM3C=JM4CΘr2.Substitute (20a)–(20d) to (17), and we can get(21a)JN1P^=JN2P^T−1Θt1T,(21b)JN3P^=JN4P^T−1Θt2T.Equivalently,(22a)JN2P^†JN1P^=T−1Θt1T,(22b)JN4P^†JN3P^=T−1Θt2T.Letψt1=JN2P^†JN1P^ and ψt2=JN4P^†JN3P^. Obviously, Θt1 and Θt2 consist of the eigenvalues of ψt1 and ψt2, respectively. Moreover, the T’s k−th column is an eigenvector that corresponds to the k−th eigenvalue of both ψt1 and ψt2. Next, we have the eigendecomposition on ψt1 so that we get the eigenvector matrix T^ as well as the associate eigenvalue matrix Θ^t1=diagλt1,1,λt1,2,…,λt1,K; that is,(23)ψt1=T^−1Θ^T^t1.Multiplyψt2 left by T^ and right by T^−1 to get the Θt2’s estimation, which is denoted by Θ^t2. Let the k−th diagonal entity of Θ^t2 be λt2,k, and then the transmit elevation angle can be obtained from Subarray 1 and Subarray 2 as(24a)θt1,kest=arcsinangleλt1,kMtπ,(24b)θt2,kest=arcsinangleλt2,kNt/π.Due to the fact that the mapping functionangle⋅ is wrapped within the range −π,π, but Ntπsinθt,k and Mtπsinθt,k are within −Ntπ,Ntπ and −Mtπ,Mtπ, the transmit elevation angles and receive elevation angles estimated from (24a24b) are ambiguous.Similarly, we have(25a)JM2P^†JM1P^=T−1Θr1T,(25b)JM4P^†JM3P^=T−1Θr2T.Define(26a)ψr1=JM2P^†JM1P^,(26b)ψr2=JM4P^†JM3P^.Obviously,(27a)Tψr1T−1=Θ,r1(27b)Tψr2T−1=Θ.r2Since the eigenvector matrixT has been estimated previously, we can get the estimation of Θr1 and Θr2 (denoted by Θ^r1 and Θ^r2, respectively) by left and right multiplying operation on ψr1 and ψr2, respectively. Let λr1,k and λr2,k be the k-th diagonal elements of Θ^r1 and Θ^r2, respectively. The receive elevation angles can be estimated via(28a)θr1,kest=arcsinangleλr1,kMrπ,(28b)θr2,kest=arcsinangleλr2,kNrπ.Similar to the estimated transmit elevation angles, due to the coprime characteristic of the two subarrays, the receiver elevation angles estimated from (28a28b) are ambiguous.


### 2.4. Unique Elevation Angle Determination

It is worth noting thatθt1,kest and θr1,kest are obtained via the rotational invariance properties of Subarray 1, while θt2,kest and θr2,kest are achieved via Subarray 2’s rotational invariance characteristics. Because of the coprime between Mt, Nt and MrNr, θt1,kest, θr1,kest and θt2,kest, θr2,kest can be uniquely determined even though they are ambiguous.Since the transmit Subarray 1’s interelement spacing isNtλ/2, there are Nt possible answers, which include one from (24a). The connection between the nt-th (nt=1,2,…,Nt) possible solution θt1,knt and the estimation θt1,kest can be expressed as(29)sinθt1,knt=sinθt1,kest−2nt−1Nt.For the transmit Subarray 2, thek-th estimated transmit elevation angle should have Mt possible solutions, including one derived from (24b); the mt−thmt=1,2,…,Mt possible solution θt2,kmt should be(30)sinθt2,kmt=sinθt2,kest−2mt−1Mt.Theoretically, the unique solution can be determined by finding the coincident solutions in (29) and (30). As displayed in Figure 2, the results of (29) and (30) are presented, where θ=30°, Mt=3, and Nt=4 are considered. In practice, the recovered angles may be closely distributed rather than coincide. Therefore, it can be estimated by the average of the closest solutions of sinθt1,knt and sinθt2,kmt; that is,(31)θ^t,k=θ^t1,k+θ^t2,k2,where θ^t1,k and θ^t2,k denote the two nearest solutions’ associate angles. Similar to the transmit, the unique receive elevation angle estimation θ^r,k can be obtained by us.Figure 2
Illustration of the uniqueness of the coprime arrayDOA=30°,Mt=3,Nt=4.


### 2.5. Azimuth Angle Estimation

To estimate the azimuth angle, we need to estimateBt and Br first. According to (17), we can estimate C via(32)C^=P^T^−1.DefineJ1=1p⊗I6N∈ℂ6M×36MN, p∈1,2,…,6M, and we can find that(33)J1C=CrDpCt.We have that(34)Br=1N∑n=1NCr6n−5:6n,:Dn−1Ar.Once the elevation angles have been estimated, the direction matricesAt and Ar can be reconstructed (denoted as A^t and A^r, respectively). Then Br can be estimated via(35)B^r=1N∑n=1NC¯r6n−5:6n,:Dn−1A^r,where C¯r=J1C^=J1P^T^−1. Thereafter, we can compute(36)fr,k=b^r1,k∗b^r2,k,where b^r1,k represents a vector that is made up of the first three entities of B^r’s k−th column and b^r2,k is a vector that consists of the last three entities of the k−th column of B^r, respectively. Finally, we can estimate ϕr,k by(37)ϕ^r,k=arctanfr,k2fr,k1,where fr,k1 is the fr,k’s first entity and fr,k2 represents the second entity of fr,k.Similarly, we can defineJ2=I6M⊗1q∈ℂ6M×36MN, q∈1,2,…,6N, and construct C¯t=J2P^. According to (35)–(37), the receive azimuth angle ϕt,k can be obtained. Since the azimuth angle is one-to-one mapping with the elevation angle and the perturbation has been compensated, therefore, the estimated angles can be matched automatically.


### 2.6. 2D-TPA and 2D-RPA Estimation

As long as we complete the 2D-DOD and 2D-DOA estimates, we can get the estimates ofDtk and Drk, represented by D^tk and D^rk, respectively. Based on equations (3a) and (3b), we can get the estimations of vtk and vrk via(38a)v^t,k=D^t,k†b^t,k,(38b)v^r,k=D^r,k†b^r,k.Finally, 2D-TPA and 2D-RPA estimation is finished via(39a)γ^t,k=arctanv^t,k2v^t,k1,(39b)η^t,k=anglev^t,k2v^t,k1,(39c)γ^r,k=arctanv^r,k2v^r,k1,(39d)η^r,k=anglev^r,k2v^r,k1.Because the premise of estimating 2D-TPA and 2D-RPA is to get 2D-DOD and 2D-DOA estimations, anglesθ^t,k,ϕ^t,k,γ^t,k,η^t,k and θ^r,k,ϕ^r,k,γ^r,k,η^r,k are automatically paired.In order to understand the algorithm proposed in this paper more efficiently, the following are the algorithm steps of the estimator proposed in this paper:Step 1: obtain the estimation of the covarianceR firstly;Step 2: calculate the propagatorP^c by (16) and construct matrix P^;Step 3: build the selection matricesJt,n and Jr,nn=1,2,3,4, and get ψt via (21a), (21b) and (22a), (22b);Step 4: obtainT^ and corresponding eigenvalues by having eigendecomposition on ψt, and get θt1/2,kest via (24a) and (24b);Step 5: obtainψr via (26a) and (26b), and get the ambiguous elevation angles in receiver through (28a) and (28b);Step 6: recover all the estimations based on (29) and (30), and determine the unique elevation from (30);Step 7: constructJ1, J2, A^r, and A^t, and then calculate C¯t and C¯r. Thereafter, recover B^r and B^t. Finally, achieve the azimuth angles via (36) and (37).Step 8: finally, obtain 2D-TPA and 2D-RPA from (38a), (38b) and (39a)–(39d).


## 2.1. Problem Formulation

As shown in Figure1, suppose that there is a bistatic EMVS-MIMO radar system, whose receiving and transmitting are distributed in the coprime geometries. In terms of the transmit array, let us assume that Subarray 1 is equipped with Mt EMVS and Subarray 2 is equipped with Nt EMVS, where Mt and Nt are coprime integers. Subarray 1’s adjacent distance is Ntλ/2 and Subarray 2’s adjacent distance is Mtλ/2, where λ represents the carrier wavelength. Similar to transmitting array, two EMVS subarrays constitute the receiving array, which have Mr and Nr elements, respectively. What is more, it is assumed that there are K far-field targets appearing in same range bin, and k−th target’s 2D-DOA pair is θr,k,ϕr,k, and 2D-DOD pair is θt,k,ϕt,k, where the elevation angles are denoted by θt,k and θr,k, and the azimuth angles are denoted by ϕt,k and ϕr,k. The matched outputs can be expressed as(1)yτ=∑k=1Kat,k⊗bt,k⊗ar,k⊗br,kskτ+nτ=At⊙Bt⊙Ar⊙Brsτ+nτ,where τ represents the slow-time index (pulse index), and sτ=s1τ,s2τ,…,sKτT represents the target reflect coefficient vector; at,k=ejNtMt−1πsinθt,k,…,ejNtπsinθt,k,1,e−Mtπsinθt,k,…,e−jMtNt−1πsinθt,kT is the k−th transmit steering vector, and ar,k=ejNrMr−1πsinθr,k,…,ejNrπsinθr,k,1,e−Mrπsinθr,k,…,e−jMrNr−1πsinθr,kT is the k−th receive steering vector; bt,k∈ℂ6×K represents the k−th polarization response vectors of transmit array, and br,k∈ℂ6×K represents the k−th polarization response vectors of the receive array; At=at,1,at,2,…,at,K∈ℂM×K, Ar=ar,1,ar,2,…,ar,K∈ℂN×K, where M=Mt+Nt−1 and N=Mr+Nr−1; Br=br,1,br,2,…,br,K∈ℂ6×K and Bt=bt,1,bt,2,…,bt,K∈ℂ6×K stand for the associate polarization response matrices, respectively. The Gaussian noise vector with zero mean and variance of σ2 is represented by nτ. Moreover,(2a)bt,k=Dt,kvt,k,(2b)br,k=Dr,kvr,k,with(3a)Dt,k=cosϕt,kcosθt,k−sinϕt,ksinϕt,kcosθt,kcosϕt,k−sinθt,k0−sinϕt,k−cosϕt,kcosθt,kcosϕt,k−sinϕt,kcosθt,k0sinθt,k,(3b)Dr,k=cosϕr,kcosθr,k−sinϕr,ksinϕr,kcosθt/r,kcosϕr,k−sinθr,k0−sinϕr,k−cosϕr,kcosθr,kcosϕr,k−sinϕr,kcosθr,k0sinθr,k,(4a)vt,k=sinγt,kejηt,kcosγt,k,(4b)vr,k=sinγr,kejηr,kcosγr,k,where γt,k and γr,k represent the k−th transmit and receive auxiliary polarization angles, respectively. The corresponding polarization phase differences are indicated, respectively, by ηt,k and ηr,k. Let btp,k and brq,k be the p-th and q-th (p,q∈1,2,…,6) entities of bt,k and br,k; that is, bt,k=bt1,k,bt2,k,…,bt6,kT and br,k=br1,k,br2,k,…,br6,kT. Define the following vectors:(5a)bt1,k=bt1,k,bt2,k,bt3,kT,(5b)bt2,k=bt4,k,bt5,k,bt6,kT,(5c)br1,k=br1,k,br2,k,br3,kT,(5d)br2,k=br4,k,br5,k,br6,kT.Figure 1
Bistatic coprime EMVS-MIMO radar model.It should be emphasized thatbt1,k and br1,k represent the electric field information vectors, while bt2,k and br2,k account for the magnetic field information vectors. In addition,(6a)bt1,k∗bt2,k=sinθt,kcosϕt,ksinθt,ksinϕt,kcosθt,k,(6b)br1,k∗br2,k=sinθr,kcosϕr,ksinθr,ksinϕr,kcosθr,k.When the targets are uncorrelated, theyτ’s covariance matrix is(7)R=CRsCH+σ2I36MN,where C=Ct⊙Cr, Ct=At⊙Bt, and Cr=Ar⊙Br; Rs=diagr1,r2,…,rK is the sτ’s covariance matrix, and rk represents the power of skτ. In practical application, the estimation of R can be obtained via L available snapshots; that is,(8)R^=1L∑τ=1LyτyHτ.Herein, estimating the 2D-DOD and 2D-DOA fromR^ is our goal.


## 2.2. Propagator

Let the firstK rows of C be C1∈ℂK×K, and let the rest be C2∈ℂ36MN−K×K; that is,(9)C=C1C2.SinceC1 is a nonsingular matrix, taking a unique linear transformation of C1 can get C2; that is,(10)C2=PcC1,where Pc∈ℂ36MN−K×K represents the propagator. Let P=IKPc; we can get(11)PC1=C1C2=C.Denote the noiseless item ofR as R¯; that is, R¯=CRsCH. Divide R¯ into two submatrices as(12)R¯=G,H,where G∈ℂ36MN×K and H∈ℂ36MN×36MN−K. Obviously,(13a)G=CRsC1H,(13b)H=CRsC2H.According to the relationship in (10), we have(14)GPcH=H.Consequently,Pc can be estimated by fitting(15)minH^−G^PcHF,where G^ is the estimation of G and H^ is the estimation of H. According to equation (15), we can obtain the Pc’s least-squares solution by(16)P^cH=G^†H^.Thereafter, we can obtain theP’s estimation via P^=IKP^c. Based on equation (11), there is a full-rank matrix T so that(17)P^=CT.Obviously,P^ plays a similar role as the signal subspace that is obtained from eigendecomposition.


## 2.3. Ambiguous Elevation Angle Estimation

In what follows, we use the rotational invariance property to obtain the estimation of the elevation anglesθt,k and θr,k. Let us define Jt,1=IMt−1,0Mt−1×Nt∈ℂMt−1×M, Jt,2=0Mt−1×1,IMt−1,0Mt−1×Nt−1∈ℂMt−1×M, Jt,3=0Nt−1×Mt−1,INt−1,0Nt−1×1∈ℂNt−1×M, Jt,4=0Nt−1×Mt,INt−1∈ℂNt−1×M, Jr,1=IMr−1,0Mr−1×Nr∈ℂMr−1×N, Jr,2=0Mr−1×1,IMr−1,0Mr−1×Nr−1∈ℂMr−1×N, Jr,3=0Nr−1×Mr−1,INr−1,0Nr−1×1∈ℂNr−1×N, and Jr,4=0Nr−1×Mr,INr−1∈ℂNr−1×N. It is easy to draw a conclusion that(18a)Jt,1At=Jt,2AtΘt1,(18b)Jt,3At=Jt,4AtΘt2,(18c)Jr,1Ar=Jr,2ArΘr1,(18d)Jr,3Ar=Jr,4ArΘr2,where Θt1=diagejNtπsinθt,1,ejNtπsinθt,2,…,ejNtπsinθt,K, Θt2=diagejMtπsinθt,1,ejMtπsinθt,2,…,ejMtπsinθt,K, Θr1=diagejNrπsinθr,1,ejNrπsinθr,2,…,ejNrπsinθr,K, and Θr2=diagejMrπsinθr,1,ejMrπsinθr,2,…,ejMrπsinθr,K. Meanwhile, we define(19a)JN1=Jt,1⊗I36N,(19b)JN2=Jt,2⊗I36N,(19c)JN3=Jt,3⊗I36N,(19d)JN4=Jt,4⊗I36N,(19e)JM1=I6M⊗Jr,1⊗I6,(19f)JM2=I6M⊗Jr,2⊗I6,(19g)JM3=I6M⊗Jr,3⊗I6,(19h)JM4=I6M⊗Jr,4⊗I6.Because there is a property ofA⊗BC⊗D=AC⊗BD, we have(20a)JN1C=JN2CΘt1,(20b)JN3C=JN4CΘt2,(20c)JM1C=JM2CΘr1,(20d)JM3C=JM4CΘr2.Substitute (20a)–(20d) to (17), and we can get(21a)JN1P^=JN2P^T−1Θt1T,(21b)JN3P^=JN4P^T−1Θt2T.Equivalently,(22a)JN2P^†JN1P^=T−1Θt1T,(22b)JN4P^†JN3P^=T−1Θt2T.Letψt1=JN2P^†JN1P^ and ψt2=JN4P^†JN3P^. Obviously, Θt1 and Θt2 consist of the eigenvalues of ψt1 and ψt2, respectively. Moreover, the T’s k−th column is an eigenvector that corresponds to the k−th eigenvalue of both ψt1 and ψt2. Next, we have the eigendecomposition on ψt1 so that we get the eigenvector matrix T^ as well as the associate eigenvalue matrix Θ^t1=diagλt1,1,λt1,2,…,λt1,K; that is,(23)ψt1=T^−1Θ^T^t1.Multiplyψt2 left by T^ and right by T^−1 to get the Θt2’s estimation, which is denoted by Θ^t2. Let the k−th diagonal entity of Θ^t2 be λt2,k, and then the transmit elevation angle can be obtained from Subarray 1 and Subarray 2 as(24a)θt1,kest=arcsinangleλt1,kMtπ,(24b)θt2,kest=arcsinangleλt2,kNt/π.Due to the fact that the mapping functionangle⋅ is wrapped within the range −π,π, but Ntπsinθt,k and Mtπsinθt,k are within −Ntπ,Ntπ and −Mtπ,Mtπ, the transmit elevation angles and receive elevation angles estimated from (24a24b) are ambiguous.Similarly, we have(25a)JM2P^†JM1P^=T−1Θr1T,(25b)JM4P^†JM3P^=T−1Θr2T.Define(26a)ψr1=JM2P^†JM1P^,(26b)ψr2=JM4P^†JM3P^.Obviously,(27a)Tψr1T−1=Θ,r1(27b)Tψr2T−1=Θ.r2Since the eigenvector matrixT has been estimated previously, we can get the estimation of Θr1 and Θr2 (denoted by Θ^r1 and Θ^r2, respectively) by left and right multiplying operation on ψr1 and ψr2, respectively. Let λr1,k and λr2,k be the k-th diagonal elements of Θ^r1 and Θ^r2, respectively. The receive elevation angles can be estimated via(28a)θr1,kest=arcsinangleλr1,kMrπ,(28b)θr2,kest=arcsinangleλr2,kNrπ.Similar to the estimated transmit elevation angles, due to the coprime characteristic of the two subarrays, the receiver elevation angles estimated from (28a28b) are ambiguous.


## 2.4. Unique Elevation Angle Determination

It is worth noting thatθt1,kest and θr1,kest are obtained via the rotational invariance properties of Subarray 1, while θt2,kest and θr2,kest are achieved via Subarray 2’s rotational invariance characteristics. Because of the coprime between Mt, Nt and MrNr, θt1,kest, θr1,kest and θt2,kest, θr2,kest can be uniquely determined even though they are ambiguous.Since the transmit Subarray 1’s interelement spacing isNtλ/2, there are Nt possible answers, which include one from (24a). The connection between the nt-th (nt=1,2,…,Nt) possible solution θt1,knt and the estimation θt1,kest can be expressed as(29)sinθt1,knt=sinθt1,kest−2nt−1Nt.For the transmit Subarray 2, thek-th estimated transmit elevation angle should have Mt possible solutions, including one derived from (24b); the mt−thmt=1,2,…,Mt possible solution θt2,kmt should be(30)sinθt2,kmt=sinθt2,kest−2mt−1Mt.Theoretically, the unique solution can be determined by finding the coincident solutions in (29) and (30). As displayed in Figure 2, the results of (29) and (30) are presented, where θ=30°, Mt=3, and Nt=4 are considered. In practice, the recovered angles may be closely distributed rather than coincide. Therefore, it can be estimated by the average of the closest solutions of sinθt1,knt and sinθt2,kmt; that is,(31)θ^t,k=θ^t1,k+θ^t2,k2,where θ^t1,k and θ^t2,k denote the two nearest solutions’ associate angles. Similar to the transmit, the unique receive elevation angle estimation θ^r,k can be obtained by us.Figure 2
Illustration of the uniqueness of the coprime arrayDOA=30°,Mt=3,Nt=4.


## 2.5. Azimuth Angle Estimation

To estimate the azimuth angle, we need to estimateBt and Br first. According to (17), we can estimate C via(32)C^=P^T^−1.DefineJ1=1p⊗I6N∈ℂ6M×36MN, p∈1,2,…,6M, and we can find that(33)J1C=CrDpCt.We have that(34)Br=1N∑n=1NCr6n−5:6n,:Dn−1Ar.Once the elevation angles have been estimated, the direction matricesAt and Ar can be reconstructed (denoted as A^t and A^r, respectively). Then Br can be estimated via(35)B^r=1N∑n=1NC¯r6n−5:6n,:Dn−1A^r,where C¯r=J1C^=J1P^T^−1. Thereafter, we can compute(36)fr,k=b^r1,k∗b^r2,k,where b^r1,k represents a vector that is made up of the first three entities of B^r’s k−th column and b^r2,k is a vector that consists of the last three entities of the k−th column of B^r, respectively. Finally, we can estimate ϕr,k by(37)ϕ^r,k=arctanfr,k2fr,k1,where fr,k1 is the fr,k’s first entity and fr,k2 represents the second entity of fr,k.Similarly, we can defineJ2=I6M⊗1q∈ℂ6M×36MN, q∈1,2,…,6N, and construct C¯t=J2P^. According to (35)–(37), the receive azimuth angle ϕt,k can be obtained. Since the azimuth angle is one-to-one mapping with the elevation angle and the perturbation has been compensated, therefore, the estimated angles can be matched automatically.


## 2.6. 2D-TPA and 2D-RPA Estimation

As long as we complete the 2D-DOD and 2D-DOA estimates, we can get the estimates ofDtk and Drk, represented by D^tk and D^rk, respectively. Based on equations (3a) and (3b), we can get the estimations of vtk and vrk via(38a)v^t,k=D^t,k†b^t,k,(38b)v^r,k=D^r,k†b^r,k.Finally, 2D-TPA and 2D-RPA estimation is finished via(39a)γ^t,k=arctanv^t,k2v^t,k1,(39b)η^t,k=anglev^t,k2v^t,k1,(39c)γ^r,k=arctanv^r,k2v^r,k1,(39d)η^r,k=anglev^r,k2v^r,k1.Because the premise of estimating 2D-TPA and 2D-RPA is to get 2D-DOD and 2D-DOA estimations, anglesθ^t,k,ϕ^t,k,γ^t,k,η^t,k and θ^r,k,ϕ^r,k,γ^r,k,η^r,k are automatically paired.In order to understand the algorithm proposed in this paper more efficiently, the following are the algorithm steps of the estimator proposed in this paper:Step 1: obtain the estimation of the covarianceR firstly;Step 2: calculate the propagatorP^c by (16) and construct matrix P^;Step 3: build the selection matricesJt,n and Jr,nn=1,2,3,4, and get ψt via (21a), (21b) and (22a), (22b);Step 4: obtainT^ and corresponding eigenvalues by having eigendecomposition on ψt, and get θt1/2,kest via (24a) and (24b);Step 5: obtainψr via (26a) and (26b), and get the ambiguous elevation angles in receiver through (28a) and (28b);Step 6: recover all the estimations based on (29) and (30), and determine the unique elevation from (30);Step 7: constructJ1, J2, A^r, and A^t, and then calculate C¯t and C¯r. Thereafter, recover B^r and B^t. Finally, achieve the azimuth angles via (36) and (37).Step 8: finally, obtain 2D-TPA and 2D-RPA from (38a), (38b) and (39a)–(39d).


## 3. Algorithm Analysis

### 3.1. Complexity

In terms of the proposed algorithm, computingP^c results in major complexity. Since the proposed algorithm can pair angle estimation automatically, the final complexity is 72MNK2. Next, the complexity of PM [16], PARAFAC [15], and ESPRIT [12] algorithms is listed in Table 1 to facilitate comparison with the proposed algorithm. As it can be seen from Table 1, compared with the ESPRIT, PM, and PARAFAC algorithms, the proposed algorithm is more efficient.Table 1
The complexity comparison.
MethodComplexityPM72MNK2+OK3ESPRITO363M3N3+362M2N2K2PARAFACO6MK2+6NK2+LK2Proposed72MNK2


### 3.2. CRB

CRB is usually used to evaluate the accuracy of parameter estimation because it gives a lower bound for estimating unbiased parameter. Define a4K×1 vector containing 2D-DODs and 2D-DOAs as Θ=θt,1,…,θt,K,θr,1,…,θr,K,ϕt,1,…,ϕt,KT. The Θ’s CRB is obtained via(40)CRB=σ22LrealC˜HΘC⊥C˜⊕RsT⊗14×4−1,where C˜=∂c1/∂θt,1,…,∂cK/∂θt,K,∂c1/∂ϕt,1,…,∂cK/∂ϕr,K, ΘC⊥=I36MN−CCHC−1CH, and ck represents the C’s k−th column vector.


## 3.1. Complexity

In terms of the proposed algorithm, computingP^c results in major complexity. Since the proposed algorithm can pair angle estimation automatically, the final complexity is 72MNK2. Next, the complexity of PM [16], PARAFAC [15], and ESPRIT [12] algorithms is listed in Table 1 to facilitate comparison with the proposed algorithm. As it can be seen from Table 1, compared with the ESPRIT, PM, and PARAFAC algorithms, the proposed algorithm is more efficient.Table 1
The complexity comparison.
MethodComplexityPM72MNK2+OK3ESPRITO363M3N3+362M2N2K2PARAFACO6MK2+6NK2+LK2Proposed72MNK2


## 3.2. CRB

CRB is usually used to evaluate the accuracy of parameter estimation because it gives a lower bound for estimating unbiased parameter. Define a4K×1 vector containing 2D-DODs and 2D-DOAs as Θ=θt,1,…,θt,K,θr,1,…,θr,K,ϕt,1,…,ϕt,KT. The Θ’s CRB is obtained via(40)CRB=σ22LrealC˜HΘC⊥C˜⊕RsT⊗14×4−1,where C˜=∂c1/∂θt,1,…,∂cK/∂θt,K,∂c1/∂ϕt,1,…,∂cK/∂ϕr,K, ΘC⊥=I36MN−CCHC−1CH, and ck represents the C’s k−th column vector.


## 4. Simulation Results

We use the Monte Carlo method to verify that the proposed method is effective in this section. In the simulation, we equip the bistatic coprime EMVS-MIMO radar withM transmit sensors and N receive sensors, respectively. The receive array and transmit array are both distributed in coprime geometries, respectively, and the subarrays’ sensor numbers are Mt, Nt (M=Mt+Nt−1) and Mr, Nr (N=Mr+Nr−1), respectively. Assume that there are K=3 and θt=40°,20°,30°, ϕt=15°,25°,35°, γt=10°,22°,35°, ηt=36°,48°,56°, θr=24°,38°,16°, ϕr=21°,32°,55°, and γr=42°,33°,60° with ηr=17°,27°,39°, respectively. In the simulation, we set p and q to be 1 and 2, respectively. The reflected coefficients are modeled with a complex Gaussian process, and there are L snapshots. Each simulation curve is based on 200 independent trials. At the same time, the root mean square error (RMSE) and average running time (ART) of elevation estimation are used to evaluate the performance. Herein, RMSE is defined as(41)RMSE=12K∑k=1K1200∑t=1200θ^t,k−θt,k2+∑k=1K1200∑t=1200θ^r,k−θt,k2,where θ^t/r,k (transmit or elevation) is the θt/r,k’s estimation of the k−th trial. For the purpose of comparison, the performances of PM [16], PARAFAC [15], and ESPRIT [12] (all these algorithms are based on ULA geometry with N receivers and M transmitters) as well as CRB are added.The RMSE performance of different methods versus Signal-to-Noise Ratio (SNR) is presented in Figure3, where we assume that Mt=3, Nt=4, Mr=4, Nr=5, and L=200 are considered. It is shown that all algorithms achieve lower RMSE with higher SNR. When the SNR is lower than −4 dB, the proposed algorithm’s performance is worse than those of the other algorithms, because it is impossible to accurately estimate the propagator in the low SNR region. However, the estimation result of the proposed algorithm is better than those of the other algorithms when the SNR is larger than −4 dB, which implies that the coprime manifold is helpful to improve estimation accuracy. Figure 4 depicts the ART performance comparison versus SNR. The results show that the proposed algorithm runs less than ESPRIT algorithm in [12], which benefits from the fact that the proposed algorithm can avoid the eigendecomposition of the high-dimensional covariance matrix. Figure 5 shows the RMSE of polarization angles estimation versus SNR, from which we can see that, with the increase of SNR, the RMSE of polarization angles estimation is decreasing.Figure 3
RMSE comparison versus SNR.Figure 4
ART comparison versus SNR.Figure 5
RMSE of polarization angles estimation versus SNR.Figure6 shows the different algorithms’ estimation performance versus diverse snapshots, where the value range of snapshots L varies from 50 to 1000 with interval of 50. Meanwhile, we suppose that Mt=3, Nt=4 and Mr=4, Nr=5 with SNR=10dB. It is displayed that the performance of angle estimation which belongs to all algorithms improves with increase of L. Notably, it is obvious that the algorithm performs better than ESPRIT, PM, and PARAFAC. Figure 7 illustrates the ART performance. The result shows that the proposed algorithm’s ART as a whole is much smaller than that of ESPRIT. Besides, the ART performance of the proposed algorithm is quite close to that of PARAFAC and PM. Moreover, with the increased number of L, the ART curve of the proposed algorithm performs basically the same as that of PM. The RMSE performance of polarization angles estimation versus L is shown in Figure 8, and it can be clearly seen that the higher L is, the lower RMSE is.Figure 6
RMSE comparison versus L.Figure 7
ART comparison versus L.Figure 8
RMSE of polarization angles estimation versus L.Figure9 presents the RMSE performance versus N, where SNR=10dB, L=200, Mt=3, Nt=4, and Mr=1, respectively. Nr is set to be varying from 3 to 27 with interval of 4, and the reason for operating above is to ensure Mr and Nr as coprime numbers. As expected, all algorithms provide better RMSE when N increases. Interestingly, when N>7, the proposed algorithm provides obviously better RMSE performance compared to PM, ESPRIT, and PARAFAC in angle estimation. However, the RMSE performance of the proposed algorithm is worse than the remaining three when N<5, which is owing to the fact that the propagator cannot be accurately estimated in the presence of smaller sensor number. Figure 10 gives the comparison results of ART with different N. Obviously, the number of N has a significant impact on ART performance. It is clearly manifested that the ART required by all algorithms increases significantly with N getting more. In addition, the ART required by the proposed algorithm is much smaller compared with those of ESPRIT, PM, and PARAFAC. The RMSE performance of polarization angles estimation is shown in Figure 11. It can be seen that RMSE decreases with the increase of N.Figure 9
RMSE comparison versus N.Figure 10
ART comparison versus N.Figure 11
RMSE of polarization angles estimation versus N.The precision performance of 2D-TPA and 2D-RPA estimation is presented in Figure12, where we assume that γt=10°,22°,35°, γr=42°,33°,60°, and ηt=36°,48°,56° with ηr=17°,27°,39°, and Mt=3, Nt=4, Mr=4, and Nr=5. In addition, L=200 and SNR=10dB are considered. It can be observed from Figure 12 that the polarization angles, namely, by 2D-TPA and 2D-RPA, are accurately identified with high precision and can be correctly paired.Figure 12
Precision performance of polarization angles estimation.


## 5. Conclusion

This paper proposes a new angle estimation algorithm for bistatic coprime MIMO radar. The core of the algorithm is to complete 2D-DOA and 2D-DOD estimation by constructing propagator. Because it does not need to involve high eigendecomposition and use coprime geometry with larger aperture, it is shown to be more computationally efficient than the state-of-the-art subspace algorithm. Compared with the existing algorithms, the proposed algorithm can achieve better estimation performance due to the larger virtual aperture. The effectiveness and improvement of the proposed algorithm are verified by numerical simulation.



---

*Source: 9954573-2021-05-19.xml*