# Angle Estimation and Mutual Coupling Self-Calibration in Bistatic MIMO System with Arbitrary Geometry

Authors: Tingping Zhang; Di Wan; Xinhai Wang; Fangqing Wen
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

Ideal array responses are often desirable to a multiple-input multiple-output (MIMO) system. Unfortunately, it may not be guaranteed in practice as the mutual coupling (MC) effects always exist. Current works concerning MC in the MIMO system only account for the uniform array geometry scenario. In this paper, we generalize the issue of angle estimation and MC self-calibration in a bistatic MIMO system in the case of arbitrary sensor geometry. The MC effects corresponding to the transmit array and the receive array are modeled by two MC matrices with several distinct entities. Angle estimation is then recast to a linear constrained quadratic problem. Inspired by the MC transformation property, a multiple signal classification- (MUSIC-) like strategy is proposed, which can estimate the direction-of-departure (DOD) and direction-of-arrival (DOA) via two individual spectrum searches. Thereafter, the MC coefficients are obtained by exploiting the orthogonality between the signal subspace and the noise subspace. The proposed method is suitable for arbitrary sensor geometry. Detailed analyses with respect to computational complexity, identifiability, and Cramer-Rao bounds (CRBs) are provided. Simulation results validate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Multiple-input multiple-output (MIMO) is the technique with the most potential for the next-generation array radar system [1–3]. A MIMO system is characterized by multiple transmitting antennas and multiple receiving antennas. Unlike the traditional phased array radar, the MIMO system emits mutual orthogonal waveforms. The spatial diversity and waveform diversity enables the MIMO system to achieve a virtual aperture, which is much larger than the physical aperture of the radar system. Benefiting from the virtual aperture, the MIMO system outperforms the phased array radar with respect to clutter suppression, direction finding, and target tracking. Usually, the MIMO system can be divided into two categories in terms of its antenna configuration, namely, a distributed MIMO system and a colocated MIMO system [4, 5]. A distributed MIMO system is equipped with widely displaced transmitting and receiving Tx/Rx antennas, and it can illuminate the area of interest from various perspectives. Therefore, it is capable of extinguishing the radar cross-section (RCS) fluctuation issue. A colocated MIMO system consists of closely spaced Tx/Rx antennas, and it is able to provide super-resolution direction estimation. Herein, we focus on the bistatic MIMO system, which is a typical representative of the colocated MIMO system.As a canonical issue in bistatic MIMO systems, joint direction-of-departure (DOD) and direction-of-arrival (DOA) estimation has been studied extensively in the past decades. Various methodologies have been proposed, e.g., estimation method of signal parameters via rotational invariance technique (ESPRIT) [6], multiple signal classification (MUSIC) [7], maximum likelihood [8], tensor approaches [9, 10], and sparsity-aware methods [11, 12]. A common assumption in [6–11] is that the sensors are well-calibrated. Nevertheless, owing to the radiation effects [13], sensors may be susceptible to unknown mutual coupling (MC), which would bring perturbation to the signal model and thus lead to degraded estimation performance. In general, the MC effect between sensors can be described by a MC matrix, in which the entities reveal the strength between two sensors. Although the MC matrix can be directly obtained offline via the measurement techniques such as the electromotive force [14], it may be invalid in the scenario where MC varies with time. To pursue online MC calibration, the usage of instrumental sensors or auxiliary sources have been reported [15, 16], which require additional resource. As is known to us, the MC intensity of two sensors is inversely proportional to their distance. Therefore, the MC effects between sensors that are far apart can be neglected, i.e., the associated MC coefficients can be approximated by zeros. As a result, sparse arrays, e.g., coprime arrays [17, 18], may be free from MC, but it is impossible for applications with limited space.As mentioned earlier, the MC intensity between two sensors is inversely proportional to their interelement distance. Accordingly, the MC matrix associated with a uniform linear array (ULA) exhibits a special structure, which can be captured by a banded symmetric Toeplitz matrix. For a ULA-based MIMO system, several strategies have been proposed for MC self-calibration. In [19], a MUSIC-like estimator was presented, in which the MC coefficients are extracted into a single vector, and the DODs and DOAs are obtained via two one-dimensional (1D) peak searches. By interpreting the sensors at both ends of the arrays as instrumental sensors, the ESPRIT-like algorithms were introduced in [20, 21], where the MC effect is suppressed by using a selection matrix, and then the ESPRIT technique is followed to obtain close-form solutions for DOD and DOA estimation. In order to exploit the tensor structure of the array data, the tensor approaches were subsequently discussed in [22–24], which offer more accurate angle estimation performance than the matrix-based ones. Besides, some frameworks have investigated into the two-dimensional (2D) angle estimation issue using a polarization-sensitive array [25–28].It should be emphasized that the above algorithms are only effective for MIMO systems with ULA configurations. In practice, the ULA-based MIMO system may occasionally encounter sensor failure [29, 30]. Besides, a more complex array manifold would be adopted in the MIMO system [31]. Therefore, the MIMO system with arbitrary Tx/Rx geometries may suffer from unknown MC, and current estimators in [19–24] will fail to work. To the best of our knowledge, only a few efforts have been devoted to robust DOA estimation for arbitrary geometry sensor array with unknown MC. In such a case, the MC matrix exhibits no special structure except Hermitian. In [32], an optimization-based framework was investigated. It formulates the DOA estimation as a block sparse inverse problem, which is solved via a joint-sparse recovery algorithm. Like the traditional counterpart, it has high complexity due to the high-dimension nonconvex problem, making it unacceptable for the MIMO system. In [33], an iterative approach was considered, in which DOA and MC coefficients are updated alternately. More recently, the MUSIC-like method was applied to robust DOA estimation in [34]. Although the approaches in [33, 34] can be directly extended to angle estimation for a bistatic MIMO system, they are inefficient due to an exhaustive grid search.In this paper, we generalize the MC issue in a bistatic MIMO system. We consider a general scenario where the Tx/Rx arrays are in arbitrary geometries and MC effects occur in both the Tx/Rx arrays. In such a scenario, the MC matrices exhibit no special structure except symmetry. To tackle the DOD and DOA estimation problem, the eigendecomposition is firstly performed to obtain the noise subspace. Then, by exploiting the orthogonality between the noise subspace and the virtual steering vector, angle estimation is recast to a general quadratic problem. Therefore, the MC transformation property is adopted to transform the optimization problem into two different forms, in which only DOD and DOA are included. Consequently, DOD and DOA can be estimated via two individual spectrum searches. Finally, the MC coefficients can be easily obtained with the estimated DOD and DOA. The proposed method is analyzed in terms of computational complexity, identifiability, and Cramer-Rao bounds (CRBs). Numerical simulations are designed to show the effectiveness of the proposed method.Notation, bold capital letters, e.g.,X, and bold lowercase letters, e.g., x, denote matrices and vectors, respectively. The superscripts ·T, ·H, and ·† account for transpose, Hermitian transpose, and inverse, respectively. The identity matrix is denoted by I, the M×N full one matrix is denoted by 1M×N, the full zero matrix is denoted by 0. ⊗, ⊙, and ⊕ represent, respectively, the Kronecker product, the Khatri-Rao product, and the Hadamard product; diag· denotes the diagonalization operation; rank· denotes rank operator; and Re· and Im· return the real part and the image part of a vector, respectively.


## 2. Problem Formulation

We consider a bistatic MIMO system withM transmit sensors and N receive sensors, and there are K far-field targets appearing in the same range bins of the radar system. The 2D-DOA pair and the 2D-DOA of the kth k=1,2,⋯,K target are denoted by Θt,k=θt,k,ϕt,kT and Θr,k=θt,k,ϕt,kT, where θt,k and θr,k are elevation angles and ϕt,k and ϕr,k are azimuth angles. Without loss of generality, we assume that both the transmit sensors and the receive sensors are distributed in arbitrary geometries. The coordinate of the mth m=1,2,⋯,M transmit sensor is pt,m=xt,m,yt,m,zt,mT, and the coordinate of the nth n=1,2,⋯,M receive sensor is pr,n=xr,n,yr,n,zr,nT. The response vectors corresponding to the transmit array and the receive array are given by(1a)aΘt,k=exp−j2πτ1,kλ,exp−j2πτ2,kλ,⋯,exp−j2πτM,KλT,(1b)aΘr,k=exp−j2π∆1,Kλ,exp−j2π∆2,Kλ,⋯,exp−j2π∆N,KλT,where λ is the carrier wavelength, τm,k and∆n,k are the associated path differences with τm,k=pt,mTft,k, ∆n,k=pr,nTfr,k, ft,k=cosϕt,ksinθt,k,sinϕt,ksinθt,k,cosθt,kT, and fr,k=cosϕr,ksinθr,k,sinϕr,ksinθr,k,cosθr,kT, respectively. DefineAt=aΘt,1,aΘt,2,⋯,aΘt,k∈CM×K and Ar=aΘr,1,aΘr,2,⋯,aΘr,k∈CN×K are the direction matrices corresponding to the transmit array and receive array, respectively. In the absence of the MC effect, the output from the matched filters can be written as(2)xt=∑k=1kaΘr,k⊗aΘt,kskt+Nt=Ar⊙Atst+nt,where t is the slow time index; skt is the RCS coefficient of the kth targets; nt is the array noise, which is assumed to be Gaussian white with variance σ2; and st=s1t,s2t,⋯,sktT. Now, we consider a scenario where both the transit sensors and the receive sensors suffer from unknown MC and the MC matrices are denoted by Ct∈CM×M and Cr∈CN×N, respectively, which are Hermitian. In the presence of an unknown MC, the matched output in (2) becomes(3)yt=∑K=1KCraΘr,k⊗CraΘt,kskt+nt=CrAr⊙CtAt⏟Ast+nt.Suppose the RCS coefficients are uncorrelated, andst is irrelevant with nt, and then, the covariance matrix of yt is(4)R=ARsAH+σ2I,where Rs=diagγ1,γ2,⋯,γk is the covariance matrix of st, γk is the power of skt. With L available snapshots, i.e., t=1,2,⋯,L, R can be estimated via(5)R^=1L∑t=1LytyHt.


## 3. The Proposed Algorithm

### 3.1. DOD and DOA Estimation

Performing eigendecomposition onR, one can obtain the signal subspace and the noise subspace as(6)R=Us,UnΛsΛnUsHUnH,where Λs∈CK×Kand Λn∈CMN−K×MN−K are two diagonal matrices that contain the K dominate eigenvalues and the remind eigenvalues, respectively. Us∈CMN×K and Un∈CMN×MN−K consist of the corresponding eigenvectors, which are called the signal subspace and the noise subspace, respectively. It is known that Us is orthogonal to Un, i.e.,(7)UnHUs=0.SinceUs span the same subspace with the columns of A, we get(8)CraΘr,k⊗CtaΘt,kHUnUnHCraΘr,k⊗CtaΘt,k=0.The above property is similar to that in the 2D-MUSIC algorithm. However, the traditional MUSIC estimator is invalid as bothCt and Cr are unknown. To further process, the following conclusion in [33] is introduced.Theorem 1.
For a matrixC∈CM×M and a vector a∈CM×1, if there are only Q (Q<M) distinct entities c=c1,c2,⋯,cQT in C, then(9)Ca=Tc,where T∈CM×Q with the qth q=1,2,⋯,Q column given by(10)T:,q=Jqa,and Jq is defined as(11)Jqm,n=1,ifCm,n=cq,0,otherwise.In this paper, we build on top of Theorem1 to address the MC self-calibration challenge in the arbitrary geometry-based MIMO system. Suppose that there are P distinct entities ct=ct,1,ct,2,⋯,ct,pT in Ct, and Q distinct entities cr=cr,1,cr,2,⋯,cr,QT in Cr. From Theorem 1, we have(12a)CtaΘt,k=TΘt,kct,(12b)CraΘr,k=TΘr,kcr,where TΘt∈CM×P and TΘr∈CN×Q are constructed according to Theorem 1. In combination with the property A⊗BC⊗D=AC⊗BD, we can obtain(13a)CraΘr,k⊗CtaΘt,k=TrΘr,k⏟Acr⏟B⊗IM⏟CCtaΘt,k⏟D=TΘr,k⊗IMcr⊗CtaΘt,k,(13b)CraΘr,k⊗CtaΘt,k=TrΘr,k⏟Acr⏟B⊗IM⏟CCtaΘt,k⏟D=IN⊗TΘt,kCraΘr,k⊗ct.Inserting (13a) and (13b) into (9) establishes(14a)eHΘt,kTΘr,k⊗IMHUnUnHTΘr,k⊗IM⏟PΘr,keΘt,k=0,(14b)fHΘr,kIN⊗TΘt,kHUnUnHIN⊗TΘt,k⏟QΘt,kfΘr,k=0,where eΘt,k=ct⊗CtaΘt,k and fΘr,k=CraΘr,k⊗cr. To obtain the DOA, one needs to optimize(15)mineHPΘes.t.,deH=ρ.Notably, (15) is a quadratic optimization issue. To avoid the trivial solution e=0MP×1, the above optimization problem is constrained with dHe=ρ, i.e.,(16)mineHPΘes.t.,dHe=ρ,where ρ is a constant and d is a vector with the first entity being one and zeros elsewhere. To solve the above issue, we can construct the following Lagrange function:(17)LΘ=eHpΘe−TdHeρ−1,where T is a Lagrange multiplier. Enforcing ∂LΘ/c to zeros yields(18)e=ρ2P−1ΘddHP−1Θd.Plug (18) into (15) and remove the constant item gives(19)maxdHP−1Θd.By setting a grid to contain all possible DODs, the indexes corresponding to the maximum values in (19) reveal the DODs. Similarly, we can obtain the DOA estimation via computing(20)maxgHQ−1Θg,where g is a NQ×1 vector with the first entity one and zeros elsewhere.


### 3.2. MC Self-Calibration

Although DODΘt,k and DOA Θr,k can be separately estimated from (19) and (20), they require further pairing, since there are KK−1/2 possible DOD-DOA pairs. Recall the properties in (12) and we can rewrite (8) as(21)ctrHTΘr,k⊗TΘt,kHUnUnHTΘr,k⊗TΘt,k⏟TΘt,k,Θr,kctr=0,where ctr=cr⊗ct. Similar to (15), to find the correct DOD-DOA pairs, we need to calculate(22)minctrHTΘt,Θrctr,which can be accomplished via(23)maxhHT−1Θt,Θrh,where h is a PQ×1 vector with the first entity one and zeros elsewhere. List all the possible DOD-DOA pairs, and the real DOD-DOA pairs can be picked out by finding the K maximum values in (23).It should be pointed out that the optimal solution of (22) is achieved if(24)ctr=α2T−1Θt,ΘrhhHT−1Θt,Θrh,where α is a constant. After the DODs and DOAs have been paired, ctr can be estimated via(25)c^t,r=1K∑k=1Kα2T−1Θ^t,k,Θ^r,khhHT−1Θ^t,k,Θ^r,kh,where Θ^t,k and Θ^r,k denote the estimated Θt,k and Θr,k, respectively. Usually, we have ct,1=1 and cr,1=1. After normalizing c^t,r, the MC vectors can be estimated via(26a)c^t=JMc^t,r,(26b)c^r=JNc^t,r,where JM=1,0⊗I∈ℝM×MN,JN=I⊗1,0M−1×1T∈ℝN×MN.


## 3.1. DOD and DOA Estimation

Performing eigendecomposition onR, one can obtain the signal subspace and the noise subspace as(6)R=Us,UnΛsΛnUsHUnH,where Λs∈CK×Kand Λn∈CMN−K×MN−K are two diagonal matrices that contain the K dominate eigenvalues and the remind eigenvalues, respectively. Us∈CMN×K and Un∈CMN×MN−K consist of the corresponding eigenvectors, which are called the signal subspace and the noise subspace, respectively. It is known that Us is orthogonal to Un, i.e.,(7)UnHUs=0.SinceUs span the same subspace with the columns of A, we get(8)CraΘr,k⊗CtaΘt,kHUnUnHCraΘr,k⊗CtaΘt,k=0.The above property is similar to that in the 2D-MUSIC algorithm. However, the traditional MUSIC estimator is invalid as bothCt and Cr are unknown. To further process, the following conclusion in [33] is introduced.Theorem 1.
For a matrixC∈CM×M and a vector a∈CM×1, if there are only Q (Q<M) distinct entities c=c1,c2,⋯,cQT in C, then(9)Ca=Tc,where T∈CM×Q with the qth q=1,2,⋯,Q column given by(10)T:,q=Jqa,and Jq is defined as(11)Jqm,n=1,ifCm,n=cq,0,otherwise.In this paper, we build on top of Theorem1 to address the MC self-calibration challenge in the arbitrary geometry-based MIMO system. Suppose that there are P distinct entities ct=ct,1,ct,2,⋯,ct,pT in Ct, and Q distinct entities cr=cr,1,cr,2,⋯,cr,QT in Cr. From Theorem 1, we have(12a)CtaΘt,k=TΘt,kct,(12b)CraΘr,k=TΘr,kcr,where TΘt∈CM×P and TΘr∈CN×Q are constructed according to Theorem 1. In combination with the property A⊗BC⊗D=AC⊗BD, we can obtain(13a)CraΘr,k⊗CtaΘt,k=TrΘr,k⏟Acr⏟B⊗IM⏟CCtaΘt,k⏟D=TΘr,k⊗IMcr⊗CtaΘt,k,(13b)CraΘr,k⊗CtaΘt,k=TrΘr,k⏟Acr⏟B⊗IM⏟CCtaΘt,k⏟D=IN⊗TΘt,kCraΘr,k⊗ct.Inserting (13a) and (13b) into (9) establishes(14a)eHΘt,kTΘr,k⊗IMHUnUnHTΘr,k⊗IM⏟PΘr,keΘt,k=0,(14b)fHΘr,kIN⊗TΘt,kHUnUnHIN⊗TΘt,k⏟QΘt,kfΘr,k=0,where eΘt,k=ct⊗CtaΘt,k and fΘr,k=CraΘr,k⊗cr. To obtain the DOA, one needs to optimize(15)mineHPΘes.t.,deH=ρ.Notably, (15) is a quadratic optimization issue. To avoid the trivial solution e=0MP×1, the above optimization problem is constrained with dHe=ρ, i.e.,(16)mineHPΘes.t.,dHe=ρ,where ρ is a constant and d is a vector with the first entity being one and zeros elsewhere. To solve the above issue, we can construct the following Lagrange function:(17)LΘ=eHpΘe−TdHeρ−1,where T is a Lagrange multiplier. Enforcing ∂LΘ/c to zeros yields(18)e=ρ2P−1ΘddHP−1Θd.Plug (18) into (15) and remove the constant item gives(19)maxdHP−1Θd.By setting a grid to contain all possible DODs, the indexes corresponding to the maximum values in (19) reveal the DODs. Similarly, we can obtain the DOA estimation via computing(20)maxgHQ−1Θg,where g is a NQ×1 vector with the first entity one and zeros elsewhere.


## 3.2. MC Self-Calibration

Although DODΘt,k and DOA Θr,k can be separately estimated from (19) and (20), they require further pairing, since there are KK−1/2 possible DOD-DOA pairs. Recall the properties in (12) and we can rewrite (8) as(21)ctrHTΘr,k⊗TΘt,kHUnUnHTΘr,k⊗TΘt,k⏟TΘt,k,Θr,kctr=0,where ctr=cr⊗ct. Similar to (15), to find the correct DOD-DOA pairs, we need to calculate(22)minctrHTΘt,Θrctr,which can be accomplished via(23)maxhHT−1Θt,Θrh,where h is a PQ×1 vector with the first entity one and zeros elsewhere. List all the possible DOD-DOA pairs, and the real DOD-DOA pairs can be picked out by finding the K maximum values in (23).It should be pointed out that the optimal solution of (22) is achieved if(24)ctr=α2T−1Θt,ΘrhhHT−1Θt,Θrh,where α is a constant. After the DODs and DOAs have been paired, ctr can be estimated via(25)c^t,r=1K∑k=1Kα2T−1Θ^t,k,Θ^r,khhHT−1Θ^t,k,Θ^r,kh,where Θ^t,k and Θ^r,k denote the estimated Θt,k and Θr,k, respectively. Usually, we have ct,1=1 and cr,1=1. After normalizing c^t,r, the MC vectors can be estimated via(26a)c^t=JMc^t,r,(26b)c^r=JNc^t,r,where JM=1,0⊗I∈ℝM×MN,JN=I⊗1,0M−1×1T∈ℝN×MN.


## 4. Algorithm Analysis

### 4.1. Computational Complexity

The computational burden (the number of complex multiplication) of the proposed is summarized as follows:Step 1.
Estimate the covariance matrixR^⋯⋯⋯M2N2L.Step 2.
Perform eigendecomposition onR^ to obtain the noise subspace Un⋯⋯⋯OM3N3.Step 3.
Calculate the spectrum function in (19) to get the DOD estimation ⋯⋯⋯OM3N2Q.Step 4.
Search the maximums of (20) to achieve the DOA estimation ⋯⋯⋯OM2N3P.


### 4.2. Identifiability

The proposed algorithm is effective ifP−1Θ, Q−1Θ, and T−1Θt,Θr are valid. As rank UnUnH=MN−K, once MQ>MN−K, PΘ is rank deficit. Similarity, if NQ>MN−K, QΘ is nonfull rank, and if PQ>MN−K, rank deficiency would occur in TΘt,Θr. As a result, a sufficient condition of the proposed estimator is that K≤minMN−P,NM−Q,MN−PQ. Since M>P and N>Q, the above condition can be simplified by K≤minMN−P,NM−Q. In other words, the proposed estimator can identify minMN−P,NM−Q targets at most.


### 4.3. Cramer-Rao Bound (CRB)

For the sake of simplicity, (3) is rewritten as(27)yt=Ast+nt,t=1,2,⋯,L,where yt≜yt,st=st,andnt=nt. Then, all the measurements are arranged into a “column” vector as y=ytT,y2T⋯,yLTT∈CMNL×1. Suppose that the st is deterministic but unknown. Thereafter, the mean μ∈CMNL×1 of y is given by(28)μ=As1As2⋮AsL=HF,where H≜I⊗A∈ℂMNL×LK and F≜s1T,s2T⋯sLTT∈ℂLK×1.Define the unknown parameter vectorsθt≜θt,1,θt,2,⋯,θt,kT, ϕt≜ϕt,1,ϕt,1⋯,ϕt,1T, θr, θr,1,θr,1,⋯,θr,kT, ϕr≜ϕr,1,ϕr,2⋯,ϕr,kT, α≜θtT,θrT,ϕtT,ϕrTT∈R4K×1, β≜RecrT,RectT,ImcrT,ImctTT∈R2P+Q×1, and γ≜ReFT,ImFTT∈R2LK×1. The whole estimation parameter vector is formulated as ζ≜αT,βT,γTT∈R2LK+2K+P+Q×1. It is well known that the CRB matrix for ζ is given by(29)CRB=σ22ReΨHΨ−1,where Ψ≜∂μ/∂αT,∂μ/∂βT,∂μ/∂γT. Next, we will concentrate on each part of Ψ. Before the detailed derivations, the following variables are defined:(30a)A~t,1≜∂aΘt,1∂θt,1,∂aΘt,2∂θt,2,⋯,∂aΘt,K∂θt,K,(30b)A~t,2≜∂aΘt,1∂ϕt,1,∂aΘt,2∂ϕt,2,⋯,∂aΘt,K∂ϕt,K,(30c)A~r,1≜∂aΘr,1∂θr,1,∂aΘr,2∂θr,2,⋯,∂aΘr,K∂θr,K,(30d)A~r,2≜∂aΘr,1∂ϕr,1,∂aΘr,2∂ϕr,2,⋯,∂aΘr,K∂ϕr,K,(30e)C~t,p≜∂Ct∂Rect,1,∂Ct∂Rect,2,⋯,∂Ct∂Rect,P,(30f)C~r,q≜∂Cr∂Recr,1,∂Cr∂Recr,2,⋯,∂Cr∂Recr,Q,(30g)A~1≜CrAr⊙CtA~t,1,(30h)A~2≜CrAr⊙CtA~t,2,(30i)A~3≜CrA~r,1⊙CtAt(30j)A~4≜CrA~r,2⊙CtAt,(30k)A~t,p≜CrA⊙C~t,pAt,(30l)A~r,q≜C~r,qAr⊙CtAt,(30m)S≜s1,s2,⋯,sLT,(30n)A~=A~1,A~2,A~3,A~4,(30o)B~=A~t,1,⋯,A~t,P,A~r,1,⋯,A~r,Q⊗1,j.It is easy to find that∂μ/∂θTT=S⊙A~1. Stepping further, we can get(31)Δ≜∂μ∂αT∈ℂMNL×4K=11×4⊗S⊙A~.As∂μ/∂Rect,p=S⊙A~t,p, and ∂μ/∂Imct,p=jS⊙A~t,p, thus, we have(32)∇≜∂μ∂βT∈ℂMNL×2P+Q=11×2P+Q⊗S⊙B~.Besides, it is straightforward to get(33)∂μ∂γT=H,jH∈ℂMNL×2LK.Consequently, it is established that(34)J≜ReΨHΨ=ReΔH∇HHH−jHHΔ,∇,H,jH.Since we are interested in CRBs on angle and MC estimation only, we will remove the uninterested parts fromJ. Define(35a)PΔ≜HHH−1HHΔ∈ℂLK×4K,(35b)P∇≜HHH−1HH∇∈ℂLK×2P+Q.BothP∇−1 and P∇−1 are valid as HHH is nonsingular. Furthermore, the following transform matrix is defined:(36)V≜I0000I00−RePΔ−ReP∇I0−ImPΔ−ImP∇0I.Thereafter, we can obtainΔ,∇,H,jHV=Δ−HPΔ,∇−HP∇,H,jH. Let ΠH⊥ denote the orthogonal projection of HH onto null space, i.e.,(37)ΠH⊥≜I−HHHH−1HH.Notably,HH∏1/H=0. Then, the following result can be obtained:(38)VHJV=RcΔHΠH⊥∇HΠH⊥HH−jHHΠH⊥Δ,ΠH⊥∇,H,jH=ReΔHΠH⊥ΔΔHΠH⊥∇00∇HΠH⊥Δ∇HΠH⊥∇0000HHHjHHH00−jHHHHHH.Applying the property of a partitioned diagonal matrix, we get(39)J−1=VVHJV−1VH=ReΔHΠH⊥ΔReΔHΠH⊥∇⊘Re∇HΠH⊥ΔRe∇HΠH⊥∇⊘⊘⊘⊘−1,where ⊘ denotes the uninteresting term. As a result, the CRB matrix concerning angle and MC estimation is given by(40)CRBa,c=σ22ReΔHΠH⊥ΔReΔHΠH⊥∇Re∇HΠH⊥ΔRe∇HΠH⊥∇−1.According to the inverse property of a partitioned matrix, the CRB with respect to angle estimationCRBa and MC estimation CRBc are separately given by(41a)CRBa=ReΔHΠH⊥Δ−ReΔHΠH⊥∇Re−1∇HΠH⊥∇Re∇HΠH⊥Δ−1,(41b)CRBc=Re∇HΠH⊥∇−Re∇HΠH⊥ΔRe−1ΔHΠH⊥ΔReΔHΠH⊥∇−1.However, it may be impractical to calculateCRBa and CRBc via (41a) and (41b) due to the high-dimensional characteristics of the matrices ∆, ∇, and H. Taking ΔHΠH⊥∇ as an example, we will show a “faster” way to calculate CRBa and CRBc. Note that H=I⊗A; hence,(42)ΠH⊥=I⊗ΠA⊥.On the other hand, both∆ and ∇ can be divided into several blocks, as illustrated in (31) and (32), respectively. Let a~u and b~v be the u-th u∈1,2,⋯,4K and v-th v∈1,2,⋯,2P+Q columns of A~ and B~, respectively. Similarly, let su, ∆u and sv, ∇v denote the u-th column and the v-th column of 11×4⊗S, 11×2P+Q⊗S and ∆ and ∇, respectively. Recall the property A⊗BC⊗D=AC⊗BD yielding(43)ΔuHΠH⊥∇v=suH⊗a~uHI⊗ΠA⊥sv⊗b~v=suHsv⊗a~uHΠA⊥b~u=Ru,v⋅a~uHΠA⊥b~u,where Ru,v is the u,v-th entity of the extended covariance matrix 14×2P+Q⊗Rs. Based on (43), we have(44)M2=ΔHΠH⊥∇=A~HΠA⊥B~⊕14×2P+Q⊗Rs.Likewise, we define(45a)M1=ΔHΠH⊥Δ=A~HΠA⊥A~⊕14×4⊗Rs,(45b)M3=∇HΠH⊥Δ=B~HΠA⊥A~⊕12P+Q×4⊗Rs,(45c)M4=∇HΠH⊥∇=A~HΠB⊥B~⊕12P+Q×2P+Q⊗Rs.Inserting (44), (45a), (45b), and (45c) into (41a) and (41b) yields(46a)CRBa=M1−M2M4−1M3−1,(46b)CRBcM4−M3M1−1M2−1.


## 4.1. Computational Complexity

The computational burden (the number of complex multiplication) of the proposed is summarized as follows:Step 1.
Estimate the covariance matrixR^⋯⋯⋯M2N2L.Step 2.
Perform eigendecomposition onR^ to obtain the noise subspace Un⋯⋯⋯OM3N3.Step 3.
Calculate the spectrum function in (19) to get the DOD estimation ⋯⋯⋯OM3N2Q.Step 4.
Search the maximums of (20) to achieve the DOA estimation ⋯⋯⋯OM2N3P.


## 4.2. Identifiability

The proposed algorithm is effective ifP−1Θ, Q−1Θ, and T−1Θt,Θr are valid. As rank UnUnH=MN−K, once MQ>MN−K, PΘ is rank deficit. Similarity, if NQ>MN−K, QΘ is nonfull rank, and if PQ>MN−K, rank deficiency would occur in TΘt,Θr. As a result, a sufficient condition of the proposed estimator is that K≤minMN−P,NM−Q,MN−PQ. Since M>P and N>Q, the above condition can be simplified by K≤minMN−P,NM−Q. In other words, the proposed estimator can identify minMN−P,NM−Q targets at most.


## 4.3. Cramer-Rao Bound (CRB)

For the sake of simplicity, (3) is rewritten as(27)yt=Ast+nt,t=1,2,⋯,L,where yt≜yt,st=st,andnt=nt. Then, all the measurements are arranged into a “column” vector as y=ytT,y2T⋯,yLTT∈CMNL×1. Suppose that the st is deterministic but unknown. Thereafter, the mean μ∈CMNL×1 of y is given by(28)μ=As1As2⋮AsL=HF,where H≜I⊗A∈ℂMNL×LK and F≜s1T,s2T⋯sLTT∈ℂLK×1.Define the unknown parameter vectorsθt≜θt,1,θt,2,⋯,θt,kT, ϕt≜ϕt,1,ϕt,1⋯,ϕt,1T, θr, θr,1,θr,1,⋯,θr,kT, ϕr≜ϕr,1,ϕr,2⋯,ϕr,kT, α≜θtT,θrT,ϕtT,ϕrTT∈R4K×1, β≜RecrT,RectT,ImcrT,ImctTT∈R2P+Q×1, and γ≜ReFT,ImFTT∈R2LK×1. The whole estimation parameter vector is formulated as ζ≜αT,βT,γTT∈R2LK+2K+P+Q×1. It is well known that the CRB matrix for ζ is given by(29)CRB=σ22ReΨHΨ−1,where Ψ≜∂μ/∂αT,∂μ/∂βT,∂μ/∂γT. Next, we will concentrate on each part of Ψ. Before the detailed derivations, the following variables are defined:(30a)A~t,1≜∂aΘt,1∂θt,1,∂aΘt,2∂θt,2,⋯,∂aΘt,K∂θt,K,(30b)A~t,2≜∂aΘt,1∂ϕt,1,∂aΘt,2∂ϕt,2,⋯,∂aΘt,K∂ϕt,K,(30c)A~r,1≜∂aΘr,1∂θr,1,∂aΘr,2∂θr,2,⋯,∂aΘr,K∂θr,K,(30d)A~r,2≜∂aΘr,1∂ϕr,1,∂aΘr,2∂ϕr,2,⋯,∂aΘr,K∂ϕr,K,(30e)C~t,p≜∂Ct∂Rect,1,∂Ct∂Rect,2,⋯,∂Ct∂Rect,P,(30f)C~r,q≜∂Cr∂Recr,1,∂Cr∂Recr,2,⋯,∂Cr∂Recr,Q,(30g)A~1≜CrAr⊙CtA~t,1,(30h)A~2≜CrAr⊙CtA~t,2,(30i)A~3≜CrA~r,1⊙CtAt(30j)A~4≜CrA~r,2⊙CtAt,(30k)A~t,p≜CrA⊙C~t,pAt,(30l)A~r,q≜C~r,qAr⊙CtAt,(30m)S≜s1,s2,⋯,sLT,(30n)A~=A~1,A~2,A~3,A~4,(30o)B~=A~t,1,⋯,A~t,P,A~r,1,⋯,A~r,Q⊗1,j.It is easy to find that∂μ/∂θTT=S⊙A~1. Stepping further, we can get(31)Δ≜∂μ∂αT∈ℂMNL×4K=11×4⊗S⊙A~.As∂μ/∂Rect,p=S⊙A~t,p, and ∂μ/∂Imct,p=jS⊙A~t,p, thus, we have(32)∇≜∂μ∂βT∈ℂMNL×2P+Q=11×2P+Q⊗S⊙B~.Besides, it is straightforward to get(33)∂μ∂γT=H,jH∈ℂMNL×2LK.Consequently, it is established that(34)J≜ReΨHΨ=ReΔH∇HHH−jHHΔ,∇,H,jH.Since we are interested in CRBs on angle and MC estimation only, we will remove the uninterested parts fromJ. Define(35a)PΔ≜HHH−1HHΔ∈ℂLK×4K,(35b)P∇≜HHH−1HH∇∈ℂLK×2P+Q.BothP∇−1 and P∇−1 are valid as HHH is nonsingular. Furthermore, the following transform matrix is defined:(36)V≜I0000I00−RePΔ−ReP∇I0−ImPΔ−ImP∇0I.Thereafter, we can obtainΔ,∇,H,jHV=Δ−HPΔ,∇−HP∇,H,jH. Let ΠH⊥ denote the orthogonal projection of HH onto null space, i.e.,(37)ΠH⊥≜I−HHHH−1HH.Notably,HH∏1/H=0. Then, the following result can be obtained:(38)VHJV=RcΔHΠH⊥∇HΠH⊥HH−jHHΠH⊥Δ,ΠH⊥∇,H,jH=ReΔHΠH⊥ΔΔHΠH⊥∇00∇HΠH⊥Δ∇HΠH⊥∇0000HHHjHHH00−jHHHHHH.Applying the property of a partitioned diagonal matrix, we get(39)J−1=VVHJV−1VH=ReΔHΠH⊥ΔReΔHΠH⊥∇⊘Re∇HΠH⊥ΔRe∇HΠH⊥∇⊘⊘⊘⊘−1,where ⊘ denotes the uninteresting term. As a result, the CRB matrix concerning angle and MC estimation is given by(40)CRBa,c=σ22ReΔHΠH⊥ΔReΔHΠH⊥∇Re∇HΠH⊥ΔRe∇HΠH⊥∇−1.According to the inverse property of a partitioned matrix, the CRB with respect to angle estimationCRBa and MC estimation CRBc are separately given by(41a)CRBa=ReΔHΠH⊥Δ−ReΔHΠH⊥∇Re−1∇HΠH⊥∇Re∇HΠH⊥Δ−1,(41b)CRBc=Re∇HΠH⊥∇−Re∇HΠH⊥ΔRe−1ΔHΠH⊥ΔReΔHΠH⊥∇−1.However, it may be impractical to calculateCRBa and CRBc via (41a) and (41b) due to the high-dimensional characteristics of the matrices ∆, ∇, and H. Taking ΔHΠH⊥∇ as an example, we will show a “faster” way to calculate CRBa and CRBc. Note that H=I⊗A; hence,(42)ΠH⊥=I⊗ΠA⊥.On the other hand, both∆ and ∇ can be divided into several blocks, as illustrated in (31) and (32), respectively. Let a~u and b~v be the u-th u∈1,2,⋯,4K and v-th v∈1,2,⋯,2P+Q columns of A~ and B~, respectively. Similarly, let su, ∆u and sv, ∇v denote the u-th column and the v-th column of 11×4⊗S, 11×2P+Q⊗S and ∆ and ∇, respectively. Recall the property A⊗BC⊗D=AC⊗BD yielding(43)ΔuHΠH⊥∇v=suH⊗a~uHI⊗ΠA⊥sv⊗b~v=suHsv⊗a~uHΠA⊥b~u=Ru,v⋅a~uHΠA⊥b~u,where Ru,v is the u,v-th entity of the extended covariance matrix 14×2P+Q⊗Rs. Based on (43), we have(44)M2=ΔHΠH⊥∇=A~HΠA⊥B~⊕14×2P+Q⊗Rs.Likewise, we define(45a)M1=ΔHΠH⊥Δ=A~HΠA⊥A~⊕14×4⊗Rs,(45b)M3=∇HΠH⊥Δ=B~HΠA⊥A~⊕12P+Q×4⊗Rs,(45c)M4=∇HΠH⊥∇=A~HΠB⊥B~⊕12P+Q×2P+Q⊗Rs.Inserting (44), (45a), (45b), and (45c) into (41a) and (41b) yields(46a)CRBa=M1−M2M4−1M3−1,(46b)CRBcM4−M3M1−1M2−1.


## 5. Simulation Results

In this section, numerical simulations are carried out to show the effectiveness of the proposed estimator. Specifically, we assume that the MIMO system is configured withM transmit sensors and N receive sensors. Suppose that there are K targets, and each RCS fulfills the Swerling II model and L=500 snapshots are collected. Search interval in the simulation is set to ∆. The following two scenarios are considered.Scenario I.
A ULA-based MIMO system with sensor failure, whereM=12 and N=12, is where both the transmitters and the receivers are ULAs with the adjacent sensor distance λ/2. There are two damaged transmit sensors with locations 7,10, and two damaged receive sensors with locations 6,11. The MC coefficient between two adjacent sensors is 0.8+0.5j, the MC coefficient of sensors with distance λ is 0.2+0.1j, and there is no MC effect when the distance of sensors is larger than λ, i.e.,ct=cr=1,0.8+0.5j,0.2+0.1jT.Suppose that K=3 targets with DOD-DOA pairs are 25∘,15∘ and 60∘,70∘.Scenario II.
A 3D-ULA-based MIMO system, whereM=N=12, wherein the adjacent sensor distance is λ/2, as illustrated in Figure 1. The MC coefficient between two adjacent sensors is 0.8+0.5j, the MC coefficient of sensors with distance λ is 0.2+0.1j, the MC coefficient between two cross-adjacent sensors is 0.017+0.035j, and there is no MC effect when the distance of sensors is larger than λ, e.g., ct=cr=1,0.8+0.5j,0.2+0.1j,0.017+0.035jT. There are K=2 targets with transmit DOD-DOA pairs 60∘,70∘ and 25∘,15∘ and receive DOD-DOA pairs 65∘,20∘ and 25∘,20∘.Scenario III.
Arbitrary Tx/Rx geometries. We assumeM=8, N=12, and the Tx/Rx sensors are randomly placed in a three-dimensional space with xt,m,yt,m,zt,m∼unif−λ,λandxr,n,yr,n,zr,n∼unif−2λ,2λ, L=500, SNR=20dB, and search interval is fixed at 0.5°. There are K=2 far-field targets, the DOD pairs are (20°, 40°) and (35°, 25°), and the DOA pairs are (50°, 10°) and (30°, 25°). Moreover, ct=cr=1,0.8+0.5j,0.017+0.035j,0.2+0.1jT.Figure 1
Illustration of 3D array.In the first example, we test the spatial spectrum of the proposed estimator in ScenarioI, where ∆=0.01°, signal-to-noise ratio (SNR) is set to 15 dB, and 10 independent trials are recorded. The results are presented in Figure 2, from which we observe that the proposed estimator can accurately identify all the targets. For comparison purposes, the spectrum results of the traditional reduced-dimension MUSIC (RD-MUSIC) algorithm in [7] are added, which display great disparities for some angles. This is caused by the fact that the traditional RD-MUSIC cannot eliminate the MC effect.Figure 2
Spatial spectrum of the proposed estimator in ScenarioI: (a) spatial spectrums of DOD estimation; (b) spectrum results of DOA estimation.
(a)(b)In the second example, we give the 2D spectrum results of the proposed estimator in ScenarioII, where ∆=0.01° and SNR=15dB are considered. Figure 3 illustrates the 2D spectrum results of the proposed estimator. Clearly, the proposed estimator can correctly recover the 2D-DODs and the 2D-DOAs, which is evident that the proposed estimator is also suitable for the MIMO system with MC and 3D sensor geometry.Figure 3
Spatial spectrum of the proposed estimator in ScenarioII: (a) spatial spectrum of 2D-DOD estimation; (b) spectrum result of 2D-DOA estimation.
(a)(b)In the third example, we plot the 2D spectrum results of the proposed estimator in ScenarioIII. Figure 4 gives the results. Similar to the previous simulation, the proposed estimator can correctly estimate the 2D-DODs and the 2D-DOAs in such a condition.Figure 4
Spatial spectrum of the proposed estimator in ScenarioII: (a) spatial spectrum of 2D-DOD estimation; (b) spectrum result of 2D-DOA estimation.
(a)(b)In the fourth example, we examine the average running time of the proposed estimator against different search interval∆ in Scenario I, where SNR is fixed at 20 dB. The result is displayed in Figure 5. As expected, a refined search interval would result in more computational burden. For comparison purposes, the performances associated with the reduced-dimension MUSIC in [32] (marked with “RD-MUSIC”) and the iterative method in [33] (marked with “Iterative method”) as well as the MUSIC-like method in [34] (marked with “MUSIC-like”) are added. In contrast with RD-MUSIC, the proposed estimator requires a slightly more running time. However, it is shown that the running time of the proposed estimator is two orders of magnitude lower than that of the Iterative method and MUSIC-like. This is because both RD-MUSIC and the proposed estimator require two 1D spectrum searches, while the 2D spectrum search is essential in both the Iterative method and MUSIC-like.Figure 5
Average running time versus∆ in Scenario I.Finally, we plot the average root mean square error (RMSE) curves of the proposed estimator versus SNR in ScenarioI, where ∆=0.01° and 200 Monte Carlo trials are carried out. As shown in Figure 6, the proposed estimator deteriorates in the low SNR regime (e.g., SNR<4dB). With the increasing SNR, the proposed estimator would achieve improved RMSE performance, and it may attain the CRB at high SNR regions. An interesting observation is that the proposed estimator offers better performance than MUSIC-like when SNR is smaller than 6 dB, and it provides very close RMSE performance to the latter when SNR is larger than 6 dB, especially in MC estimation. Nevertheless, both RD-MUSIC and Iterative method provide the unacceptable RMSE performance over the entire SNR regions.Figure 6
RMSE performance versusL in Scenario I: (a) RMSE on angle estimation; (b) RMSE on MC coefficient estimation.
(a)(b)


## 6. Conclusion

In this paper, we presented a MUSIC-like estimator for bistatic MIMO radar with arbitrary geometries and MC. For the 1D sensor manifold scenario, the proposed estimator only needs two one-dimensional spectrum peak searches to obtain 1D-DODs and 1D-DOAs, while for 2D or 3D sensor geometries, the proposed estimator allows two 2D spectrum searches to achieve 2D-DODs and 2D-DOAs. The proposed estimator is efficient from the viewpoint of complexity. It should have a brighter prospect in numerous applications, e.g., military explorations and internet-of-vehicles.



---

*Source: 9984589-2021-07-26.xml*