# Joint DOD and DOA Estimation of Bistatic MIMO Radar for Coprime Array Based on Array Elements Interpolation

Authors: Zhiyuan You; Guoping Hu; Guimei Zheng; Hao Zhou
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This paper proposes a method to address the problem of the joint direction of departure (DOD) and direction of arrival (DOA) estimation with augmented coprime array (CPA) bistatic multiple input multiple output (MIMO) radar using interpolating sensors. At first, we deduce the regular pattern of hole positions in the virtual array and interpolate a small number of sensors to augmented CPA to form a partially contiguous virtual array. Then, we use the diversity smoothing algorithm to reconstruct the Toeplitz matrix to obtain a spatial smoothing matrix. Finally, we combine the RD-MUSIC algorithm with the spatial smoothing matrix to estimate the spatial spectrum and achieve automatic matching of DODs and DOAs for the targets. Simulation results clearly illustrate the superiority of the method.


---

## Body

## 1. Introduction

Direction of arrival (DOA) estimation and Kalman filter algorithm for target tracking [1] are significant areas of research in array signal processing. It can estimate the angular position of different signals in a certain airspace and plays a key role in radar, sonar, and other target detection fields. Since the last century, DOA estimation has gone through three stages, which are beamforming [2], subspace fitting [3–6], and compressed sensing [7], and made a great contribution to the development of the DOA estimation algorithm. Bistatic multiple input multiple output (MIMO) radar has the advantages of waveform diversity [8] and spatial separation; it estimates target angle including DOA and direction of departure (DOD), so the joint estimation of DOD and DOA [9, 10] has also become the focus of research. Although the uniform linear array (ULA) can also solve the related problems of DOA estimation, it still has certain limitations. For example, the number of physical array sensors usually cannot be less than the number of signal sources; otherwise, the estimation accuracy will be badly affected, and thus the degree of freedom (DOF) is limited by physical array elements and other factors [11]. In order to solve these problems, a sparse array is introduced by scholars. The difference between the sparse array [12, 13] and the traditional uniform array mainly includes that the sparse array is formed by setting different interelement spacings in the array to form a sparse array structure, and part of the sparse array has a larger element spacing, so it forms longer virtual array to increase DOF. Coprime array (CPA) [14, 15] has been widely studied by scholars because of its larger element spacing. Thus, the CPA has a low mutual coupling effect and high accuracy of target angle estimation. Moreover, the CPA can construct a longer virtual array through its own differential coarray [16] and use its equivalent vector to estimate DOA. However, although the virtual array formed by CPA has a long virtual array aperture [17], there are array holes that make the virtual array discontinuous, which limits the expansion of DOF. Scholars have conducted a series of studies to settle the question.A method was advanced to fill the virtual array hole by moving the sparse array in [18], and it used the differential coarray of the original array and the differential coarray of the moved array to form a composite array without holes to estimate DOA effectively. But this method has a small estimated error that originated from array movement distance and signal source movement distance, and its theoretical feasibility is insufficient. Interpolated virtual array element was considered in [19] and solved the problem of virtual array discontinuity caused by virtual array holes, but the method will be affected when there are continuous virtual holes, and its DOF will reduce. Some scholars utilize the virtual array element interpolation theory [19] to the joint estimation of DOD and DOA of bistatic radar targets; this method reduces redundancy of the virtual array and improved DOF and estimated target resolution, but this method required more calculations. Combined with the traditional ESPRIT algorithm [20], Li et al. proposed a joint estimation of target DOD and DOA with bistatic CPA MIMO radar based on virtual aperture expansion in [21]. This method achieved better estimation performance than the traditional method, but it did not use discrete virtual array elements. Recently, a method which interpolated array elements in CPA was proposed to solve the virtual array hole problem in [22], but augmented CPA which the method used cannot nicely reflected the advantage of the interpolated array elements to expand the virtual array aperture. We propose a joint DOD and DOA estimation of bistatic MIMO radar for coprime array based on array elements interpolation. The method uses augmented CPAs as the transmit array and receive array and interpolates a small number of sensors to the particular holes in the virtual array to expand the aperture. Then, the method uses the selection matrix to reconstruct the Toeplitz matrix based on diversity smoothing to estimate the DODs and DOAs of the sources.The remaining sections are as follows: in Section2, we reduce the math model of the bistatic CPA MIMO radar. We deduce the law of hole position and propose a diversity smoothing algorithm for reconstructing the Toeplitz matrix to estimate the DODs and DOAs of the targets in Section 3. We perform the same simulation experiments with the proposed method and other methods, which clearly illustrates the superiority of the method in Section 4. Finally, we present conclusions for this paper in Section 5.Notations: we use italicized boldface characters to represent vectors and matrices in this paper. Superscripts.T and .H represent transpose and conjugate-transpose, respectively, diag [·] denotes diagonal matrix, and ⊗ and ∘ denote the Kronecker product and the Hadamard product, respectively.


## 2. Mathematical Model

The conventional bistatic CPA MIMO radar model is presented in Figure1. Transmit array and receive array in the model are augmented CPA which consists of two uniform linear arrays (ULA). One ULA has 2M1 − 1 sensors in transmit array, and red circles represent sensors of subarray 1; another ULA of transmit array has N1 sensors, and the black circle represents sensors of subarray 2. Similarly, receive array contains two ULAs, which, respectively, have 2M2 − 1 sensors and N2 sensors, where M1 and N1 are two coprime integers and M2 and N2 are two coprime integers. The number of sensors of arrays are, respectively, M = 2M1 + N1 − 1 and N = 2M2 + N2 − 1. The unit interelement spacing of the array is d, which is the half wavelength (λ/2). The sensor positions are given by the following equation:(1)St=−mN1d|0≤m≤2M1−1∪2nM1d|0≤n≤N1−1,Sr=−mN2d|0≤m≤2M2−1∪2nM2d|0≤n≤N2−1.Figure 1
Conventional bistatic CPA MIMO radar model.Suppose there areK uncorrelated signals in the space, the DOD and DOA of signals are given by φ=φ1,φ2,…,φK and θ=θ1,θ2,…,θK. pt=pt1,pt2,…,ptM represents the position of the sensors in the transmit array, and pr=pr1,pr2,…,prN denotes the position of the sensors in the receive array. Set the reflection coefficient of signals as st=s1t,s2t,…,sKtT, t=1,2,…,L, and L represents the number of snapshots. The received signal after matched filtering is given by the following equation:(2)xt=atφ1⊗arθ1,…,atφK⊗arθKst+nt=At∘Arst+nt=Ast+nt,where n(t) represents a matrix composed of Gaussian white noise and it follows the Gaussian distribution nt∼0,σ2. At and Ar are also given by the following equation:(3)Ar=arθ1,arθ2,…,arθk,…,arθK,At=atφ1,atφ2,…,atφk,…,atφK,where atφk and arθk, respectively, denote manifold matrices, and they are given by the following equation:(4)atφk=e−j2π/λpt1sinφk,…,e−j2π/λptMsinφkT,arθk=e−j2π/λpr1sinθk,…,e−j2π/λprNsinθkT.Therefore, the covariance matrix of the received signalR is given by the following equation:(5)R=ExtxHt=ARsAH+σn2IMN,where IMN represents an MN × MN dimensional identity matrix and Rs denotes the covariance matrix of the received target.(6)Rs=EstsHt=diagσ12,σ22,…,σk2,…,σK2,where σk2 represents the power of the kth target signal.


## 3. Joint Diversity Smoothing DOD and DOA Estimation Algorithm Based on Interpolated Sensors

### 3.1. Expansion of Virtual Array by Interpolating Sensors to CPA Holes

In this section, we will interpolate sensors to CPA holes to expand a longer virtual array aperture and then we can obtain a new equivalent vector to execute transmit-receive diversity smoothing. Assume transmit array and receive array are identical augmented CPA, and taking transmit array as an example, sensors distribution is illustrated in Figure2.Figure 2
Sensors distribution of augmented CPA.The virtual array is formed by the sum-difference array of CPA, and the position of virtual sensorsSc is given by the following equation:(7)Sc=Sc|Sc=Sct+Scr−S˜ct+S˜cr,Sct,S˜ct∈St,Scr,S˜cr∈Sr=Sc|Sc=Sct−S˜ct+Scr−S˜cr,Sct,S˜ct∈St,Scr,S˜cr∈Sr.According to (8), the virtual sensors are only generated by the difference coarray of transmit array and the difference coarray of receive array when DOD and DOA are estimated, respectively. Therefore, Sct−S˜ct corresponds to parameter φk and Scr−S˜cr corresponds to parameter θk.The position of the virtual sensors that are formed by the difference coarray of transmit array is given by the following equation:(8)Std=Std|Std=±M1nd+N1md,0≤m≤2M1−1,0≤n≤N1−1,where the value range of Std is [−(3M1N1 − M1 − N1), 3M1N1 − M1 − N1].Similarly, the position of the virtual sensors that are formed by the difference coarray of receive array is given by the following equation:(9)Srd=Srd|Srd=±M2nd+N2md,0≤m≤2M2−1,0≤n≤N2−1,where the value range of Srd is [−(3M2N2 − M2 − N2), 3M2N2 − M2 − N2].For transmit array,M1 and M2 can determine the distribution of the virtual sensors and the length of the virtual array aperture, and the distribution of virtual sensors that are formed by augmented CPA is shown in Figure 3.Figure 3
Distribution of virtual sensors and virtual holes.In Figure3, black filled circles represent virtual sensors and dotted circles represent virtual holes. Although augmented CPA can form more virtual sensors than traditional CPA, some virtual holes are located next to the center point. Therefore, other continuous virtual sensors cannot form a virtual array with a larger aperture.With the purpose of settling the question, we present a method that utilizes a small number of sensors to interpolate a part of virtual holes. The originally discontinuous virtual array will become continuous to achieve more DOF.We can know a special kind of holes next to the center point in Figure3 and call them central virtual holes. Central virtual holes are continuous virtual holes which locate next to the center point. When M1 = 2 and M2 = 3, there are two central virtual holes in the virtual array; when M1 = 3 and M2 = 4, there are four central virtual holes in the virtual array; and when M1 = 3 and M2 = 5, there are four central virtual holes in the virtual array.Assume we use sensors to interpolate the virtual central holes where next to the center point completely. Then, a new distribution of virtual sensors emerges in Figure4.Figure 4
The new distribution of virtual array after interpolating sensors.From Figure4, black filled circles represent virtual sensors, dotted circles represent virtual holes, red filled circles represent sensors, and red circles represent new formed virtual sensors. We found that we only need to interpolate the virtual holes located next to the center point to make the virtual array continuous. For example, when M1 = 3 and M2 = 4, we use two sensors to interpolate the virtual central holes which are located in {1d,2d}. Then, the array can form additional new virtual sensors which are located in {−1d,−2d,±5d}. Most of the original virtual holes are interpolated. Comparing Figure 3 with Figure 4, it can be found that the interpolated virtual array is a continuous virtual array that all the central virtual holes are replaced with sensors or new virtual sensors, and the DOF is 47. Therefore, we deduce the number of central virtual holes and DOF for different augmented CPAs from continuous experiments and find that the number of virtual holes to be interpolated is only related to M1. The derived regular pattern and formulas are shown in Table 1.Table 1
Number of elements for augmented CPA.
M1N1MHsHfM˜DOFN1 = M1 + 12362123233494247474512637979…………………N1 = M1 + 235104359595716841391397922126251251…………………N1 = M1 + 347146311111158178415915971023126279279…………………N1 = M1 + LM1N12M1 + N1 − 12(M1 − 1)M1 − 14M1N1 − 14M1N1 − 1M denotes the number of array sensors, Hs denotes the number of central virtual holes located next to the center point, Hf denotes the number of interpolated sensors, M˜ denotes the number of continuous virtual sensors, and DOF denotes the DOF after interpolating the holes.In general, we interpolate (M1 − 1) sensors into the virtual holes which located next to the center point, the original virtual array becomes continuous, virtual array apertures become larger to (4M1N1 − 1), and DOF also becomes larger to (4M1N1 − 1).


### 3.2. Reconstructing Toeplitz Matrix Algorithm Based on Diversity Smoothing

In recent years, scholars have proposed some spatial smoothing algorithms for bistatic radars in [23–25]. This paper chooses an algorithm that reconstructs the Toeplitz matrix based on diversity smoothing.After the sensors shave been interpolated, there areM˜ continuous virtual sensors in the virtual array of transmit and N˜ continuous virtual sensors in the virtual array of receive, and pt˜=pt1,pt2,…,ptM˜ denotes the sensors’ position in the virtual transmit array and pr˜=pr1,pr2,…,prN˜ denotes the sensors’ position in the virtual receive array. Therefore, the steering matrices of transmit and receive virtual array are given by the following equation:(10)A˜t˜=a˜t˜φ1,a˜t˜φ2,…,a˜t˜φk,A˜r˜=a˜r˜θ1,a˜r˜θ2,…,a˜r˜θk,where a˜t˜φk and a˜r˜θk, respectively, denote manifold matrices, and they are given by the following equation:(11)a˜t˜φk=e−j2π/λpt1sinφk,…,e−j2π/λptM˜sinφkT,a˜r˜θk=e−j2π/λpr1sinθk,…,e−j2π/λprN˜sinθkT.A new covariance matrixR˜ of the received signal can be formed by processing continuous virtual sensors. However, the ideal R˜ is hard to get in practice, so R^ is usually estimated by using L available snapshots and it is given by the following equation:(12)R^=1L∑i=1Lx^tx^Ht,where x^t is the new received signal of the virtual array.The new covariance matrix is now vectorized to obtain a new equivalent vector as follows:(13)r˜=A˜p+σn2i˜=A˜t∘A˜rp+σn2i˜,where p=σ12,σ22,…,σK2 and i˜ is a M˜N˜×1 vector whose elements are all zeros except for the M˜N˜+1/2th row.In fact, we regard vectorr˜ as the received signal of a single snapshot in bistatic MIMO radar. According to the principles of spatial smoothing algorithms, the number of smoothing times should be no less than the number of overlapping subarrays to form a full rank covariance matrix. Since the number of virtual sensors is always an odd number, the best smoothing results are achieved when the number of smoothing times is equal to the number of overlapping subarrays. Therefore, we divide the virtual array of transmit into Ms overlapping subarrays to carry out Ms smoothing and divide the virtual array of receive into Ns overlapping subarrays to carry out Ns smoothing. Thus, M˜=2Ms−1 and N˜=2Ns−1.Next, suppose there are two selection matricesSmt and Snr corresponding to the virtual array of transmit and the virtual array of receive, which are given by the following equations:(14)Smt=OMs×Ms−mIMs×MsOMs×m−1,(15)Snr=ONs×Ns−nINs×NsONs×n−1,where m=1,…,Ms and n=1,…,Ns. OMs×Ms−m denote the Ms×Ms−m zero matrix, IMs×Ms denote Ms×Ms identity matrix, and OMs×m−1 denote the Ms×m−1 zero matrix. Snr is similar.Utilize selection matricesSmt and Snr to form the transmission smoothing and receiving smoothing on the equivalent vector r˜. Then, we can change the values of m and n to achieve the smoothing of overlapping subarrays in the virtual array, m and n represent smoothing times. So subvector r˜m,n is given by the following equation:(16)r˜m,n=Smt⊗Snrr˜=SmtA˜t∘SnrA˜rp+σn2i˜m,n=A¯t∘A¯rΨtm−1Ψrn−1p+σn2i˜m,n,where A¯t=S1tA˜t, Ψt=diage−jπsinφ1,…,e−jπsinφk, A¯r=S1rA˜r, Ψr=diage−jπsinθ1,…,e−jπsinθk, and i˜m,n is the corresponding noise vector. The reconstructed Toeplitz matrix is given by the following equation:(17)RT=r˜1,1,r˜1,2,…,r˜1,Ns,r˜2,1,…,r˜Ms,Ns=A¯t∘A¯rRsA¯t∘A¯rH+σn2IMsNs=AsRsAsH+σn2IMsNs,where As represents the new steering matrix.Therefore, the covariance matrix of each sub-array is given by the following equation:(18)R¯m,n=r˜m,nr˜m,nH=AsΨtm−1Ψrn−1ppHΨt1−mΨr1−nAsH+σn4i˜m,ni˜m,nH+σn2ΑsΨtm−1Ψrn−1pi˜m,nH+σn2i˜m,npHΨtm−1Ψrn−1AsH.The values ofm and n are changed to achieve the effect of spatial smoothing in the virtual array and the spatial smoothing algorithm is used to solve the single snapshot problem caused by the vectorized covariance matrix. The spatial smoothing matrix is given by the following equation:(19)R¯=1MsNs∑m=1Ms∑n=1NsR¯m,n=1MsNsAsRsAsHAsRsAsH+σn4IMsNs+2σn2AsRsAsH.Through comparison, we find thatR¯ and Toeplitz matrix RT have the same form, so the spatial smoothing matrix R¯ is also given by the following equation:(20)R¯=1MsNsRT2.When performing spatial spectrum estimation, we choose the RD-MUSIC algorithm instead of the ESPRIT algorithm. Although we use the ESPRIT algorithm for spatial spectrum estimation to effectively reduce the complexity of the algorithm, the estimated DODs and DOAs of multiple targets need to be matched manually, whereas the RD-MUSIC algorithm can automatically match the DODs and DOAs of targets.The RD-MUSIC algorithm is combined withR¯ to estimate the spatial spectrum, and the spatial spectrum function is given by the following equation:(21)fmusic=1a^tφ⊗a^rθHEnEnHa^tφ⊗a^rθ,where En denotes signal subspace of the covariance matrix. a^tθ⊗a^rθ is also given by the following equation [25]:(22)a^tφ⊗a^rθ=a^tφ⊗INa^rθ.Therefore, the spatial spectrum function is also given by the following equation:(23)fmusic=1a^rθHa^tφ⊗INHEnEnHa^tφ⊗INa^rθ=1a^rθHVφa^rθ,where Vφ=a^tφ⊗INHEnEnHa^tφ⊗IN. The DOA of the kth signal is given by the following equation:(24)θ^k=argmin1e1TV−1φke1=argmaxe1TV−1φke1,where e1=1,0,…,0T.Similarly, The DOD of the kth signal is given by the following equation:(25)φ^k=argmin1e1TV−1θke1=argmaxe1TV−1θke1.The steps of the proposed Algorithm 1 are as follows.Algorithm 1: Array elements interpolation algorithm.
Input: receive signal:x(t) = As(t)+n(t), t = 1,2, …, L;Output: DODs and DOAs:φ^k,θ^k,k=1,2,…,K;Step:(1)
estimate covariance matrixR^ based on continuous virtual array as in (12);(2)
build the selection matrix of the transmit arraySmt as in (14) and the selection matrix of receive array Snr as in (15);(3)
rebuild a new Toeplitz matrixRT as in (17);(4)
form a spatial smoothing matrixR¯ as in (18);(5)
utilize the RD-MUSIC algorithm, to estimate DODs and DOAsφ^k,θ^k;


## 3.1. Expansion of Virtual Array by Interpolating Sensors to CPA Holes

In this section, we will interpolate sensors to CPA holes to expand a longer virtual array aperture and then we can obtain a new equivalent vector to execute transmit-receive diversity smoothing. Assume transmit array and receive array are identical augmented CPA, and taking transmit array as an example, sensors distribution is illustrated in Figure2.Figure 2
Sensors distribution of augmented CPA.The virtual array is formed by the sum-difference array of CPA, and the position of virtual sensorsSc is given by the following equation:(7)Sc=Sc|Sc=Sct+Scr−S˜ct+S˜cr,Sct,S˜ct∈St,Scr,S˜cr∈Sr=Sc|Sc=Sct−S˜ct+Scr−S˜cr,Sct,S˜ct∈St,Scr,S˜cr∈Sr.According to (8), the virtual sensors are only generated by the difference coarray of transmit array and the difference coarray of receive array when DOD and DOA are estimated, respectively. Therefore, Sct−S˜ct corresponds to parameter φk and Scr−S˜cr corresponds to parameter θk.The position of the virtual sensors that are formed by the difference coarray of transmit array is given by the following equation:(8)Std=Std|Std=±M1nd+N1md,0≤m≤2M1−1,0≤n≤N1−1,where the value range of Std is [−(3M1N1 − M1 − N1), 3M1N1 − M1 − N1].Similarly, the position of the virtual sensors that are formed by the difference coarray of receive array is given by the following equation:(9)Srd=Srd|Srd=±M2nd+N2md,0≤m≤2M2−1,0≤n≤N2−1,where the value range of Srd is [−(3M2N2 − M2 − N2), 3M2N2 − M2 − N2].For transmit array,M1 and M2 can determine the distribution of the virtual sensors and the length of the virtual array aperture, and the distribution of virtual sensors that are formed by augmented CPA is shown in Figure 3.Figure 3
Distribution of virtual sensors and virtual holes.In Figure3, black filled circles represent virtual sensors and dotted circles represent virtual holes. Although augmented CPA can form more virtual sensors than traditional CPA, some virtual holes are located next to the center point. Therefore, other continuous virtual sensors cannot form a virtual array with a larger aperture.With the purpose of settling the question, we present a method that utilizes a small number of sensors to interpolate a part of virtual holes. The originally discontinuous virtual array will become continuous to achieve more DOF.We can know a special kind of holes next to the center point in Figure3 and call them central virtual holes. Central virtual holes are continuous virtual holes which locate next to the center point. When M1 = 2 and M2 = 3, there are two central virtual holes in the virtual array; when M1 = 3 and M2 = 4, there are four central virtual holes in the virtual array; and when M1 = 3 and M2 = 5, there are four central virtual holes in the virtual array.Assume we use sensors to interpolate the virtual central holes where next to the center point completely. Then, a new distribution of virtual sensors emerges in Figure4.Figure 4
The new distribution of virtual array after interpolating sensors.From Figure4, black filled circles represent virtual sensors, dotted circles represent virtual holes, red filled circles represent sensors, and red circles represent new formed virtual sensors. We found that we only need to interpolate the virtual holes located next to the center point to make the virtual array continuous. For example, when M1 = 3 and M2 = 4, we use two sensors to interpolate the virtual central holes which are located in {1d,2d}. Then, the array can form additional new virtual sensors which are located in {−1d,−2d,±5d}. Most of the original virtual holes are interpolated. Comparing Figure 3 with Figure 4, it can be found that the interpolated virtual array is a continuous virtual array that all the central virtual holes are replaced with sensors or new virtual sensors, and the DOF is 47. Therefore, we deduce the number of central virtual holes and DOF for different augmented CPAs from continuous experiments and find that the number of virtual holes to be interpolated is only related to M1. The derived regular pattern and formulas are shown in Table 1.Table 1
Number of elements for augmented CPA.
M1N1MHsHfM˜DOFN1 = M1 + 12362123233494247474512637979…………………N1 = M1 + 235104359595716841391397922126251251…………………N1 = M1 + 347146311111158178415915971023126279279…………………N1 = M1 + LM1N12M1 + N1 − 12(M1 − 1)M1 − 14M1N1 − 14M1N1 − 1M denotes the number of array sensors, Hs denotes the number of central virtual holes located next to the center point, Hf denotes the number of interpolated sensors, M˜ denotes the number of continuous virtual sensors, and DOF denotes the DOF after interpolating the holes.In general, we interpolate (M1 − 1) sensors into the virtual holes which located next to the center point, the original virtual array becomes continuous, virtual array apertures become larger to (4M1N1 − 1), and DOF also becomes larger to (4M1N1 − 1).


## 3.2. Reconstructing Toeplitz Matrix Algorithm Based on Diversity Smoothing

In recent years, scholars have proposed some spatial smoothing algorithms for bistatic radars in [23–25]. This paper chooses an algorithm that reconstructs the Toeplitz matrix based on diversity smoothing.After the sensors shave been interpolated, there areM˜ continuous virtual sensors in the virtual array of transmit and N˜ continuous virtual sensors in the virtual array of receive, and pt˜=pt1,pt2,…,ptM˜ denotes the sensors’ position in the virtual transmit array and pr˜=pr1,pr2,…,prN˜ denotes the sensors’ position in the virtual receive array. Therefore, the steering matrices of transmit and receive virtual array are given by the following equation:(10)A˜t˜=a˜t˜φ1,a˜t˜φ2,…,a˜t˜φk,A˜r˜=a˜r˜θ1,a˜r˜θ2,…,a˜r˜θk,where a˜t˜φk and a˜r˜θk, respectively, denote manifold matrices, and they are given by the following equation:(11)a˜t˜φk=e−j2π/λpt1sinφk,…,e−j2π/λptM˜sinφkT,a˜r˜θk=e−j2π/λpr1sinθk,…,e−j2π/λprN˜sinθkT.A new covariance matrixR˜ of the received signal can be formed by processing continuous virtual sensors. However, the ideal R˜ is hard to get in practice, so R^ is usually estimated by using L available snapshots and it is given by the following equation:(12)R^=1L∑i=1Lx^tx^Ht,where x^t is the new received signal of the virtual array.The new covariance matrix is now vectorized to obtain a new equivalent vector as follows:(13)r˜=A˜p+σn2i˜=A˜t∘A˜rp+σn2i˜,where p=σ12,σ22,…,σK2 and i˜ is a M˜N˜×1 vector whose elements are all zeros except for the M˜N˜+1/2th row.In fact, we regard vectorr˜ as the received signal of a single snapshot in bistatic MIMO radar. According to the principles of spatial smoothing algorithms, the number of smoothing times should be no less than the number of overlapping subarrays to form a full rank covariance matrix. Since the number of virtual sensors is always an odd number, the best smoothing results are achieved when the number of smoothing times is equal to the number of overlapping subarrays. Therefore, we divide the virtual array of transmit into Ms overlapping subarrays to carry out Ms smoothing and divide the virtual array of receive into Ns overlapping subarrays to carry out Ns smoothing. Thus, M˜=2Ms−1 and N˜=2Ns−1.Next, suppose there are two selection matricesSmt and Snr corresponding to the virtual array of transmit and the virtual array of receive, which are given by the following equations:(14)Smt=OMs×Ms−mIMs×MsOMs×m−1,(15)Snr=ONs×Ns−nINs×NsONs×n−1,where m=1,…,Ms and n=1,…,Ns. OMs×Ms−m denote the Ms×Ms−m zero matrix, IMs×Ms denote Ms×Ms identity matrix, and OMs×m−1 denote the Ms×m−1 zero matrix. Snr is similar.Utilize selection matricesSmt and Snr to form the transmission smoothing and receiving smoothing on the equivalent vector r˜. Then, we can change the values of m and n to achieve the smoothing of overlapping subarrays in the virtual array, m and n represent smoothing times. So subvector r˜m,n is given by the following equation:(16)r˜m,n=Smt⊗Snrr˜=SmtA˜t∘SnrA˜rp+σn2i˜m,n=A¯t∘A¯rΨtm−1Ψrn−1p+σn2i˜m,n,where A¯t=S1tA˜t, Ψt=diage−jπsinφ1,…,e−jπsinφk, A¯r=S1rA˜r, Ψr=diage−jπsinθ1,…,e−jπsinθk, and i˜m,n is the corresponding noise vector. The reconstructed Toeplitz matrix is given by the following equation:(17)RT=r˜1,1,r˜1,2,…,r˜1,Ns,r˜2,1,…,r˜Ms,Ns=A¯t∘A¯rRsA¯t∘A¯rH+σn2IMsNs=AsRsAsH+σn2IMsNs,where As represents the new steering matrix.Therefore, the covariance matrix of each sub-array is given by the following equation:(18)R¯m,n=r˜m,nr˜m,nH=AsΨtm−1Ψrn−1ppHΨt1−mΨr1−nAsH+σn4i˜m,ni˜m,nH+σn2ΑsΨtm−1Ψrn−1pi˜m,nH+σn2i˜m,npHΨtm−1Ψrn−1AsH.The values ofm and n are changed to achieve the effect of spatial smoothing in the virtual array and the spatial smoothing algorithm is used to solve the single snapshot problem caused by the vectorized covariance matrix. The spatial smoothing matrix is given by the following equation:(19)R¯=1MsNs∑m=1Ms∑n=1NsR¯m,n=1MsNsAsRsAsHAsRsAsH+σn4IMsNs+2σn2AsRsAsH.Through comparison, we find thatR¯ and Toeplitz matrix RT have the same form, so the spatial smoothing matrix R¯ is also given by the following equation:(20)R¯=1MsNsRT2.When performing spatial spectrum estimation, we choose the RD-MUSIC algorithm instead of the ESPRIT algorithm. Although we use the ESPRIT algorithm for spatial spectrum estimation to effectively reduce the complexity of the algorithm, the estimated DODs and DOAs of multiple targets need to be matched manually, whereas the RD-MUSIC algorithm can automatically match the DODs and DOAs of targets.The RD-MUSIC algorithm is combined withR¯ to estimate the spatial spectrum, and the spatial spectrum function is given by the following equation:(21)fmusic=1a^tφ⊗a^rθHEnEnHa^tφ⊗a^rθ,where En denotes signal subspace of the covariance matrix. a^tθ⊗a^rθ is also given by the following equation [25]:(22)a^tφ⊗a^rθ=a^tφ⊗INa^rθ.Therefore, the spatial spectrum function is also given by the following equation:(23)fmusic=1a^rθHa^tφ⊗INHEnEnHa^tφ⊗INa^rθ=1a^rθHVφa^rθ,where Vφ=a^tφ⊗INHEnEnHa^tφ⊗IN. The DOA of the kth signal is given by the following equation:(24)θ^k=argmin1e1TV−1φke1=argmaxe1TV−1φke1,where e1=1,0,…,0T.Similarly, The DOD of the kth signal is given by the following equation:(25)φ^k=argmin1e1TV−1θke1=argmaxe1TV−1θke1.The steps of the proposed Algorithm 1 are as follows.Algorithm 1: Array elements interpolation algorithm.
Input: receive signal:x(t) = As(t)+n(t), t = 1,2, …, L;Output: DODs and DOAs:φ^k,θ^k,k=1,2,…,K;Step:(1)
estimate covariance matrixR^ based on continuous virtual array as in (12);(2)
build the selection matrix of the transmit arraySmt as in (14) and the selection matrix of receive array Snr as in (15);(3)
rebuild a new Toeplitz matrixRT as in (17);(4)
form a spatial smoothing matrixR¯ as in (18);(5)
utilize the RD-MUSIC algorithm, to estimate DODs and DOAsφ^k,θ^k;


## 4. Simulation Results

We assume the distance of array sensorsM1 = M2 = 2, N1 = N2 = 3, the number of transmit and receive array sensors M = N = 6, and the original positions of transmit and receive sensors are [−9d, −6d, −3d, 0d, 2d, 4d]. The number of central virtual holes to be interpolated is Hf = 1, and the positions of transmit and receive sensors after interpolating the holes are [−9d, −6d, −3d, 0d, 1d, 2d, 4d]. There are continuous virtual sensors in transmit array and receive array, overlapping subarrays of transmit and receive array, and smooth times are Ms = Ns = 12. The signals in the simulation are incoherent.Compare the performance in target estimation of four methods which include diversity smoothing algorithm based on interpolating sensors (DSIAS) in this article, interpolation virtual array element algorithm (IVAE) in [19], conventional CPA diversity smoothing algorithm (CDS) in [21], and conventional ULA RD-MUSIC algorithm in [26]. It includes three specific simulation experiments, namely, target number detection and angular resolution. In the simulation of root mean square error (RMSE), we not only compared the above-given methods, we also added a unitary dual-resolution ESPRIT (U-ESPRIT) method [27] to the comparative simulation.


### 4.1. Number of Detectable Targets

In this part, the signal-to-noise ratio (SNR) is set as 10 dB and the number of snapshots is set as 200. There are 15 signal targets distributed over the range [−70°, 70°], where located atφ= [−70°， −60°, −50°, −40°, −30°, −20°, −10°, 0°, 10°, 20°, 30°, 40°, 50°, 60°,70°] and θ= [−70°, −60°, −50°, −40°, −30°, −20°, −10°, 0°, 10°, 20°, 30°, 40°, 50°, 60°,70°]. Figure 5 shows the results of DOD and DOA joint estimation of targets. The four pictures in Figure 5 represent the algorithms in this article, [18, 20], and [25], respectively. As shown in Figure 5, it shows the contour map of the spatial spectrum peak, red circles represent the true direction of 15 targets and the spectral peak contour represents the estimated direction of 15 targets. The spectral peak contour completely overlaps in the red circle in Figure 5(a), so the DSIAS algorithm can accurately estimate 15 targets. Other methods’ spectral peak contours do not overlap completely in the red circle. The DSIAS algorithm is better than other methods for a number of detectable targets.Figure 5
Results of detecting 15 targets by different methods: (a) CPA with DSIAS algorithm, (b) CPA with IVAE algorithm in [19], (c) CPA with CDS algorithm in [21], and (d) ULA with RD-MUSIC algorithm in [26].
(a)(b)(c)(d)


### 4.2. Angular Resolution

In this part, with the purpose of comparing the performance of the four algorithms in angular resolution, assume that the SNR is 10 dB and the number of snapshots is 200. Figure6 shows situations of the angular resolution comparison of different algorithms. The position of two adjacent targets are φ1,θ1=11∘,11∘ and φ2,θ2=13∘,13∘, and the red line denotes the real target's direction. From the simulation results, we can also see that the MIMO radar which uses the DSIAS algorithm can estimate the targets well under the condition of two targets close to each other, and the other three algorithms cannot accurately distinguish two similar targets.Figure 6
Results of the angular resolution of different methods: (a) CPA with DSIAS algorithm, (b) CPA with IVAE algorithm in [19], (c) CPA with CDS algorithm in [21], and (d) ULA with RD-MUSIC algorithm in [26].
(a)(b)(c)(d)


### 4.3. Root Mean Square Error (RMSE)

We compare the RMSEs of different algorithms. RMSE is a common standard that reflects the accuracy of angle estimation, and the average RMSE is defined by the following equation:(26)RMSE=12×QK∑i=1Q∑k=1Kφ^ki−φk2+θ^ki−θk2,where Q denotes Monte Carlo simulation times, K denotes the number of targets, and φ^ki,θ^ki denotes the joint estimated DOD and DOA of the kth target for the ith Monte Carlo simulation, i = 1,2, …, Q. Estimation number of targets will affect the estimation accuracy. In order to reflect the comprehensiveness of the simulation results, we carried out simulation experiments for two targets and four targets, respectively, and the results of the simulation experiments are shown in Figures 7 and 8.Figure 7
RMSE versus SNR and number of snapshots for different methods (two targets): (a) RMSE versus SNR for different methods and (b) RMSE versus a number of snapshots for different methods.
(a)(b)Figure 8
RMSE versus SNR and number of snapshots for different methods (four targets): (a) RMSE versus SNR for different methods and (b) RMSE versus a number of snapshots for different methods.
(a)(b)Figure7 shows the relationship of RMSE with SNR and the number of snapshots under the condition of detecting two targets, and the targets locate at φ1,θ1=10∘,15∘, φ2,θ2=20∘,25∘. Figure 7(a) depicts the variation of the RMSE curve with SNR, where the number of snapshots is set as 200. The DSIAS algorithm has a higher estimation accuracy than other methods at low SNR but has similar performance with the other two methods [19, 21] at high SNR. Figure 7(b) depicts the variation of the RMSE curve with a number of snapshots, where the SNR is set as 10. The DSIAS algorithm has high estimation accuracy at a different number of snapshots.As shown in Figure7, the accuracy of all four algorithms is high when estimating two target angles, and the estimation accuracy of the DSIAS algorithm proposed in this paper is only slightly higher than the others. Figure 8 presents the relationship of RMSE with SNR and the number of snapshots under the condition of detecting four targets, and the targets locate at φ1,θ1=10∘,15∘, φ2,θ2=20∘,25∘, φ3,θ3=30∘,35∘, and φ4,θ4=40∘,45∘. Figure 8(a) depicts the variation of the RMSE curve with SNR, where the number of snapshots is set as 200. In contrast to the previous experiments, the DSIAS algorithm clearly performs better than other methods in estimation accuracy. Figure 8(b) depicts the variation of the RMSE curve with a number of snapshots, where the SNR is 10. It can be found that the DSIAS algorithm estimates more targets with greater accuracy by comparing two experiments which detects a different number of targets.


## 4.1. Number of Detectable Targets

In this part, the signal-to-noise ratio (SNR) is set as 10 dB and the number of snapshots is set as 200. There are 15 signal targets distributed over the range [−70°, 70°], where located atφ= [−70°， −60°, −50°, −40°, −30°, −20°, −10°, 0°, 10°, 20°, 30°, 40°, 50°, 60°,70°] and θ= [−70°, −60°, −50°, −40°, −30°, −20°, −10°, 0°, 10°, 20°, 30°, 40°, 50°, 60°,70°]. Figure 5 shows the results of DOD and DOA joint estimation of targets. The four pictures in Figure 5 represent the algorithms in this article, [18, 20], and [25], respectively. As shown in Figure 5, it shows the contour map of the spatial spectrum peak, red circles represent the true direction of 15 targets and the spectral peak contour represents the estimated direction of 15 targets. The spectral peak contour completely overlaps in the red circle in Figure 5(a), so the DSIAS algorithm can accurately estimate 15 targets. Other methods’ spectral peak contours do not overlap completely in the red circle. The DSIAS algorithm is better than other methods for a number of detectable targets.Figure 5
Results of detecting 15 targets by different methods: (a) CPA with DSIAS algorithm, (b) CPA with IVAE algorithm in [19], (c) CPA with CDS algorithm in [21], and (d) ULA with RD-MUSIC algorithm in [26].
(a)(b)(c)(d)


## 4.2. Angular Resolution

In this part, with the purpose of comparing the performance of the four algorithms in angular resolution, assume that the SNR is 10 dB and the number of snapshots is 200. Figure6 shows situations of the angular resolution comparison of different algorithms. The position of two adjacent targets are φ1,θ1=11∘,11∘ and φ2,θ2=13∘,13∘, and the red line denotes the real target's direction. From the simulation results, we can also see that the MIMO radar which uses the DSIAS algorithm can estimate the targets well under the condition of two targets close to each other, and the other three algorithms cannot accurately distinguish two similar targets.Figure 6
Results of the angular resolution of different methods: (a) CPA with DSIAS algorithm, (b) CPA with IVAE algorithm in [19], (c) CPA with CDS algorithm in [21], and (d) ULA with RD-MUSIC algorithm in [26].
(a)(b)(c)(d)


## 4.3. Root Mean Square Error (RMSE)

We compare the RMSEs of different algorithms. RMSE is a common standard that reflects the accuracy of angle estimation, and the average RMSE is defined by the following equation:(26)RMSE=12×QK∑i=1Q∑k=1Kφ^ki−φk2+θ^ki−θk2,where Q denotes Monte Carlo simulation times, K denotes the number of targets, and φ^ki,θ^ki denotes the joint estimated DOD and DOA of the kth target for the ith Monte Carlo simulation, i = 1,2, …, Q. Estimation number of targets will affect the estimation accuracy. In order to reflect the comprehensiveness of the simulation results, we carried out simulation experiments for two targets and four targets, respectively, and the results of the simulation experiments are shown in Figures 7 and 8.Figure 7
RMSE versus SNR and number of snapshots for different methods (two targets): (a) RMSE versus SNR for different methods and (b) RMSE versus a number of snapshots for different methods.
(a)(b)Figure 8
RMSE versus SNR and number of snapshots for different methods (four targets): (a) RMSE versus SNR for different methods and (b) RMSE versus a number of snapshots for different methods.
(a)(b)Figure7 shows the relationship of RMSE with SNR and the number of snapshots under the condition of detecting two targets, and the targets locate at φ1,θ1=10∘,15∘, φ2,θ2=20∘,25∘. Figure 7(a) depicts the variation of the RMSE curve with SNR, where the number of snapshots is set as 200. The DSIAS algorithm has a higher estimation accuracy than other methods at low SNR but has similar performance with the other two methods [19, 21] at high SNR. Figure 7(b) depicts the variation of the RMSE curve with a number of snapshots, where the SNR is set as 10. The DSIAS algorithm has high estimation accuracy at a different number of snapshots.As shown in Figure7, the accuracy of all four algorithms is high when estimating two target angles, and the estimation accuracy of the DSIAS algorithm proposed in this paper is only slightly higher than the others. Figure 8 presents the relationship of RMSE with SNR and the number of snapshots under the condition of detecting four targets, and the targets locate at φ1,θ1=10∘,15∘, φ2,θ2=20∘,25∘, φ3,θ3=30∘,35∘, and φ4,θ4=40∘,45∘. Figure 8(a) depicts the variation of the RMSE curve with SNR, where the number of snapshots is set as 200. In contrast to the previous experiments, the DSIAS algorithm clearly performs better than other methods in estimation accuracy. Figure 8(b) depicts the variation of the RMSE curve with a number of snapshots, where the SNR is 10. It can be found that the DSIAS algorithm estimates more targets with greater accuracy by comparing two experiments which detects a different number of targets.


## 5. Conclusions

We propose a joint estimation of the DOD and DOA method that interpolate a small number of sensors to a specific location in an augmented CPA virtual array. The method can expand the aperture of the virtual array and maintains the maximum DOF of augmented CPA despite the virtual holes which cannot be exploited in the virtual array. Meanwhile, we reconstructed the Toeplitz matrix based on diversity smoothing to obtain a spatial smoothing matrix. Finally, we combined the RD-MUSIC algorithm with the spatial smoothing matrix to estimate the spatial spectrum and accurately estimate the DOD and DOA of the targets. Simulation results illustrate the proposed method has better performance than other methods.



---

*Source: 3483778-2022-05-21.xml*