# Computationally Efficient 2D DOA Estimation for L-Shaped Array with Unknown Mutual Coupling

Authors: Yang-Yang Dong; Xin Chang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Although L-shaped array can provide good angle estimation performance and is easy to implement, its two-dimensional (2D) direction-of-arrival (DOA) performance degrades greatly in the presence of mutual coupling. To deal with the mutual coupling effect, a novel 2D DOA estimation method for L-shaped array with low computational complexity is developed in this paper. First, we generalize the conventional mutual coupling model for L-shaped array and compensate the mutual coupling blindly via sacrificing a few sensors as auxiliary elements. Then we apply the propagator method twice to mitigate the effect of strong source signal correlation effect. Finally, the estimations of azimuth and elevation angles are achieved simultaneously without pair matching via the complex eigenvalue technique. Compared with the existing methods, the proposed method is computationally efficient without spectrum search or polynomial rooting and also has fine angle estimation performance for highly correlated source signals. Theoretical analysis and simulation results have demonstrated the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Two-dimensional (2D) direction-of-arrival (DOA) estimation is an important area of array signal processing and has wide applications in wireless communication, radar, sonar, electronic warfare, and so far [1, 2]. In order to obtain unambiguous angle estimation, many researchers often restrict the minimum sensor spacing to not more than half-wavelength, which results in the mutual coupling effect [3, 4]. Various researches have indicated that mutual coupling changes the array manifold matrix and causes severe angle estimation performance degradation [5, 6]. To mitigate the mutual coupling effect, researchers have developed many 2D DOA estimation methods with unknown mutual coupling [3–14]. These methods can mainly be classified into three types: electromagnetic simulation [3, 4, 6], active calibration [7, 8], and blind calibration [5, 9–14]. The electromagnetic simulation methods in [3, 4] use the electromagnetic theory to calculate the mutual coupling and are high in accuracy. What is more, the mutual coupling coefficients are often easy to be affected by the environment and cannot be as stable as a constant, which causes the fact that this kind of methods cannot be used for the application fields with drastic change on environment, such as antiradiation missile (ARM) [2]. Meanwhile, for active calibration methods [7, 8], they require at least one calibration source to estimate the mutual coupling, where the additional active calibration sources are often not available in practice [9].While for blind calibration algorithms, they are very promising since there is no requirement for stable environment or calibration sources and 2D DOAs in the presence of unknown mutual coupling can be estimated only via signal processing. There are three primary techniques to deal with the unknown mutual coupling. The first one is realized by using iterative procedure to compensate the mutual coupling and estimate 2D DOAs [5], which can be applied for many array geometries. However, it sometimes converges slowly and may cause the increase in runtime and even wrong DOA estimations [9]. The second one is called RAnk REduction (RARE) technique [9, 12, 14], which exploits the special structure of coupled array manifold matrix to construct a cost function similar to that of MUSIC method and estimate the angles via minimizing the cost function. However, the angle estimation via RARE technique often needs multidimensional search, which is still high in computational complexity [14]. The third technique is named auxiliary sensor technique, which blindly compensates the mutual coupling at the cost of a few auxiliary sensors and generates new received array data [10, 13]. Since the data often possesses special structure and is much simpler compared to the RARE technique, many computationally efficient 2D DOA estimation methods available can be applied via a few simple modifications. Therefore, compared to the other two techniques, the auxiliary sensor technique is more computationally efficient.The blind calibration methods above are developed for various 2D arrays. Although Rübsamen and Gershman have developed novel sparse 2D arrays which can provide nonambiguous DOA estimates for the full360∘ azimuth field-of-view [15], the structures of the sparse 2D arrays are too complicated for hardware manufacturers to implement. Among the common arrays for 2D DOA estimation available and in consideration of the angle estimation performance and implementability, L-shaped array is of the optimal implementability [16]. Hence, 2D DOA estimation for L-shaped array with unknown mutual coupling is essential. Wu et al. [12] have proposed a RARE technique based method, which requires two-dimensional search and is still too high in computational complexity for real-time applications [2]. Luckily, Liang et al. [13] have developed an auxiliary sensor based method, which can estimate 2D angles in the presence of mutual coupling without spectrum search or pair matching, which is much lower than Wu et al.’s method [12] in terms of computational complexity. However, Liang et al.’s method has three main drawbacks: (1) it can only hold for the case that the maximum mutual coupling degree is equal to 2, which restricts its application in more general cases; (2) it cannot deal with highly correlated source signals, which is an unavoidable problem for many real applications, such as the electronic reconnaissance with deception jamming [2]; (3) it divides the 2D DOA estimation into two sequential steps, which causes the error accumulation effect and has a bad effect on angle estimations.To handle the above problems, in this paper, we first generalize the data model for any mutual coupling degree. Inspired by Liang et al.’s method [13], the proposed method obtains the new received array data via auxiliary sensor technique mentioned above. Then, we utilize the propagator method twice to decrease the source signal correlation and form two matrices associated with 2D DOAs which have the same eigenvectors and the different eigenvalues. Finally, the 2D DOAs are estimated simultaneously via the complex eigenvalue technique [17]. Simulation experiments and theoretical analysis have demonstrated the validity of the proposed method and proved that it can handle the three problems encountered in Liang et al.’s method [13]. What is more, compared with Wu et al.’s method [12] and Liang et al.’s method [13], the proposed method can handle the angle estimation for highly correlated source signals well with lower computational complexity.Notations. Matrices and vectors are denoted by boldfaced capital letters and lowercase letters, respectively. The superscripts (·)∗,(·)T,(·)H,(·)-1, and (·)† stand for conjugate, transpose, conjugate transpose, inverse, and pseudoinverse. The notations E{·}, ⊗, var⁡(·), ID, JD, 0m×n, diag⁡{·}, Re⁡{·}, and Im⁡{·} denote the statistical expectation, Kronecker product, the variance of a random variable, a D×D identity matrix, a D×D exchange matrix with the ones on its antidiagonal and zeros elsewhere, an m×n zero matrix, diagonalization, and the real and imaginary part of a complex number, respectively.


## 2. Signal Model

Two uniform linear orthogonal arrays with intersensor spacingd form an x-z plane L-shaped array configuration as shown in Figure 1 and the total number of sensors is 2M+1. We assume K noncoherent narrowband far-field signals {sk(n)}k=1K(n=1,2,…,N, where N denotes the number of snapshots) of wavelength λ from distinct directions with azimuth and elevation angles {(θk,ϕk)}k=1K (0∘≤θk≤180∘, 0∘≤ϕk≤180∘) impinging on the L-shaped array.Figure 1
L-shaped array configuration for 2D DOA estimation.Letαk=-2πdcos⁡θksin⁡ϕk/λ and βk=-2πdcos⁡ϕk/λ; the received array data vector at the nth snapshot can be expressed as [13](1)yn=CAsn+wn=∑k=1KCaαk,βkskn+wn,where C, A, a(αk,βk), s(n), and w(n) represent the unknown mutual coupling matrix, ideal manifold matrix, ideal steering vector, source signal vector, and noise vector, respectively. In detail, y(n)=[y1,0(n)⋯yM,0(n)y0,0(n)y0,1(n)⋯y0,M(n)]T, A=[a(α1,β1),a(α2,β2),…,a(αK,βK)], a(αk,βk)=[ejαk,…,ejMαk,1,ejβk,…,ejMβk]T, s(n)=[s1(n),s2(n),…,sK(n)]T, w(n)=[w1,0(n)⋯wM,0(n)w0,0(n)w0,1(n)⋯w0,M(n)]T.In [13], C is constructed via assuming that each sensor is only affected by its adjacent sensors within the distance of 5d, which is not general. To construct a more general mutual coupling model for L-shaped array, we firstly define the maximum coupling range as pd (p is a positive integer). Then, with the fact that the mutual coupling matrix of ULA can be modeled as a banded symmetric Toeplitz matrix, C can be expressed as a block symmetric matrix (2)C=FgHgT1gTHgF,where the mutual coupling matrix within each subarray is F=Toeplitz{[1,c1,…,cp,01×(M-p-1)]}, the mutual coupling matrix corresponding to the element at the origin is g=[c1,…,cp,01×(M-p)]T, the mutual coupling matrix between two subarrays is {Hl,m=Hm,l,m2+l2≤p;Hl,m=0,m2+l2>p}. (The subscript (m,l) stands for the sensor with coordinates (md,0,ld).)


## 3. Proposed Method

### 3.1. Blind Mutual Coupling Effect Compensation

Inspired by the analysis done in [13], let us define two selection matrices firstly (3)T1=0M+1-2p×p-1IM+1-2p0M+1-2p×M+1+p,T2=0M+1-2p×M+pIM+1-2p0M+1-2p×p,where T1 and T2 satisfy T1T1H=IM+1-2p and T1T2H=0(M+1-2p)×(M+1-2p).ApplyingT1 to (1), we have(4)xn=T1yn=AxΛxsn+wxn,where Ax=[ax(α1),ax(α2),…,ax(αK)] and ax(αk)=[1,ejαk,…,ej(M-2p)αk]T. Λx=diag⁡{γ1x,γ2x,…,γKx} and γkx=ejpαk+∑l=1p[cp-l+1ej(l-1)αk+clej(p+l)αk].Similarly, applyingT2 to (1), we can obtain(5)zn=T2yn=AzΛzsn+wzn,where Az=[az(β1),az(β2),…,az(βK)] and az(βk)=[1,ejβk,…,ej(M-2p)βk]T. Λz=diag⁡{γ1z,γ2z,…,γKz} and γkz=ejpβk+∑l=1p[cp-l+1ej(l-1)βk+clej(p+l)βk].Hence, the mutual coupling effect is compensated blindly at the cost of a few auxiliary sensors.


### 3.2. Joint Azimuth and Elevation Angle Estimation

With (4) and (5), Liang et al.’s method in [13] can be applied with a few modifications. However, Liang et al.’s method cannot handle source signals with moderate or high correlation and has the error accumulation effect for 2D DOA estimations.To solve the problems encountered in [13], we firstly apply the propagator method (PM) [18] to the new received array data y~(n)=[xT(n),zT(n)]T and obtain its corresponding propagator matrix P as follows:(6)P=R~1HR~1-1R~1HR~2,where R~1 and R~2 are two parts of the covariance matrix R~=E{y~(n)y~H(n)}=A~RssA~H+σw2I2M-4p+2=[R~1,R~2] and R~1 consists of first K columns of R~. Additionally, Rss=E{s(n)sH(n)} and A~=[(AxΛx)T,(AzΛz)T]T.Then, to exploit all array information available fromy~(n), we define a modified propagator matrix Pm=[IK,P]H satisfying(7)PmA~1=A~,where A~=[A~1T,A~2T]T. A~1 and A~2 are the first K and the last (2M-4p+2-K) rows of A~.To make full use of the potential rotational invariance property inPm, carry out the following calculations: (8)Q1=G3PmG1Pm†=A~z1ΛzΛx-1A~x1†,Q2=G3PmG2Pm†=A~z1Φ∗αΛzΛx-1A~x1†,Q3=G4PmG1Pm†=A~z1ΦβΛzΛx-1A~x1†,Q4=G4PmG2Pm†=A~z1Φ∗αΦβΛzΛx-1A~x1†,where G1=[IM-2p,0(M-2p)×(M-2p+2)], G2=[0(M-2p)×1,IM-2p,0(M-2p)×(M-2p+1)], G3=[0(M-2p)×(M-2p+1),IM-2p,0(M-2p)×1], and G4=[0(M-2p)×(M-2p+2),IM-2p]. Φ(α)=diag⁡{ejα1,…,ejαK}, Φ(β)=diag⁡{ejβ1,…,ejβK}, A~x1Λx=G1A~, and A~z1Λz=G3A~. From (8), it is noticed that they can always hold if Pm exists (i.e., the source signals are noncoherent), while (14)–(17) in [13] cannot hold for highly correlated source signals, because Rss cannot be approximated to a diagonal matrix for the signal correlation factor ρ>0.3. That is to say, (8) decreases the correlation effect.Next, to take full advantage of the information provided by (8), construct the following augmented matrix: (9)Raug=Q1,JM-2pQ4∗Q2,JM-2pQ3∗Q3,JM-2pQ2∗Q4,JM-2pQ1∗=A~z1A~z1Φ∗αA~z1ΦβA~z1Φ∗αΦβA~x1†,ΦαΦM-2p-1βΛzΛx-1A~x1†∗=AaugSaug.It is noticed that the construction of Rnew is different from the authors’ earlier work [19], where [19] just increases the effective array aperture, while the approach utilized here increases the effective array aperture and virtual snapshots simultaneously.Similar to (6), calculate the propagator matrix Paug corresponding to Raug as follows:(10)Paug=R1augR1augH-1R1augR2augH,where R1aug and R2aug are the first K rows and the last (4M-8p-K) rows of Raug. Define the modified propagator matrix Pmaug=[IK,Paug]H and do the calculations as follows: (11)Ψα=H1Pmaug†H2Pmaug=A1augΦαA1aug-1,Ψβ=H3Pmaug†H4Pmaug=A1augΦβA1aug-1,where H1=I2⊗[0(M-2p)×(M-2p),IM-2p], H2=I2⊗[IM-2p,0(M-2p)×(M-2p)], H3=[I2M-4p,0(2M-4p)×(2M-4p)], and H4=[0(2M-4p)×(2M-4p),I2M-4p]. A1aug is the first K rows of Aaug. It is noticed that Ψ(α) and Ψ(β) have the same eigenvectors and the different eigenvalues, which is a general problem in multidimensional spectral and array signal processing [17].To estimate azimuth and elevation simultaneously, constructΨ~ using the complex eigenvalue technique [17], (12)Ψ~=Ψα+IK-1Ψα-IK-jΨβ+IK-1Ψβ-IK=A1augΦ~A1aug-1,where Φ~=diag⁡{tan⁡(β1/2)+jtan⁡(α1/2),…,tan⁡(βK/2)+jtan⁡(αK/2)}. Then applying eigenvalue decomposition to Ψ~, we can obtain Φ~. Therefore, the azimuth and elevation angles can be estimated, (13)ϕ^k=arccos⁡-β^kλ2πd,(14)θ^k=arccos⁡-α^kλ2πdsin⁡ϕ^k,where α^k=2arctan⁡[Im(Φ~(k,k))], and β^K=2arctan⁡[Re⁡(Φ~(k,k))]. Φ~(k,k) denotes the kth row and kth column element of Φ~.From (11) to (14), we know that the proposed method is able to estimate the angles simultaneously, which avoids the error accumulation effect that occurs in [13].


## 3.1. Blind Mutual Coupling Effect Compensation

Inspired by the analysis done in [13], let us define two selection matrices firstly (3)T1=0M+1-2p×p-1IM+1-2p0M+1-2p×M+1+p,T2=0M+1-2p×M+pIM+1-2p0M+1-2p×p,where T1 and T2 satisfy T1T1H=IM+1-2p and T1T2H=0(M+1-2p)×(M+1-2p).ApplyingT1 to (1), we have(4)xn=T1yn=AxΛxsn+wxn,where Ax=[ax(α1),ax(α2),…,ax(αK)] and ax(αk)=[1,ejαk,…,ej(M-2p)αk]T. Λx=diag⁡{γ1x,γ2x,…,γKx} and γkx=ejpαk+∑l=1p[cp-l+1ej(l-1)αk+clej(p+l)αk].Similarly, applyingT2 to (1), we can obtain(5)zn=T2yn=AzΛzsn+wzn,where Az=[az(β1),az(β2),…,az(βK)] and az(βk)=[1,ejβk,…,ej(M-2p)βk]T. Λz=diag⁡{γ1z,γ2z,…,γKz} and γkz=ejpβk+∑l=1p[cp-l+1ej(l-1)βk+clej(p+l)βk].Hence, the mutual coupling effect is compensated blindly at the cost of a few auxiliary sensors.


## 3.2. Joint Azimuth and Elevation Angle Estimation

With (4) and (5), Liang et al.’s method in [13] can be applied with a few modifications. However, Liang et al.’s method cannot handle source signals with moderate or high correlation and has the error accumulation effect for 2D DOA estimations.To solve the problems encountered in [13], we firstly apply the propagator method (PM) [18] to the new received array data y~(n)=[xT(n),zT(n)]T and obtain its corresponding propagator matrix P as follows:(6)P=R~1HR~1-1R~1HR~2,where R~1 and R~2 are two parts of the covariance matrix R~=E{y~(n)y~H(n)}=A~RssA~H+σw2I2M-4p+2=[R~1,R~2] and R~1 consists of first K columns of R~. Additionally, Rss=E{s(n)sH(n)} and A~=[(AxΛx)T,(AzΛz)T]T.Then, to exploit all array information available fromy~(n), we define a modified propagator matrix Pm=[IK,P]H satisfying(7)PmA~1=A~,where A~=[A~1T,A~2T]T. A~1 and A~2 are the first K and the last (2M-4p+2-K) rows of A~.To make full use of the potential rotational invariance property inPm, carry out the following calculations: (8)Q1=G3PmG1Pm†=A~z1ΛzΛx-1A~x1†,Q2=G3PmG2Pm†=A~z1Φ∗αΛzΛx-1A~x1†,Q3=G4PmG1Pm†=A~z1ΦβΛzΛx-1A~x1†,Q4=G4PmG2Pm†=A~z1Φ∗αΦβΛzΛx-1A~x1†,where G1=[IM-2p,0(M-2p)×(M-2p+2)], G2=[0(M-2p)×1,IM-2p,0(M-2p)×(M-2p+1)], G3=[0(M-2p)×(M-2p+1),IM-2p,0(M-2p)×1], and G4=[0(M-2p)×(M-2p+2),IM-2p]. Φ(α)=diag⁡{ejα1,…,ejαK}, Φ(β)=diag⁡{ejβ1,…,ejβK}, A~x1Λx=G1A~, and A~z1Λz=G3A~. From (8), it is noticed that they can always hold if Pm exists (i.e., the source signals are noncoherent), while (14)–(17) in [13] cannot hold for highly correlated source signals, because Rss cannot be approximated to a diagonal matrix for the signal correlation factor ρ>0.3. That is to say, (8) decreases the correlation effect.Next, to take full advantage of the information provided by (8), construct the following augmented matrix: (9)Raug=Q1,JM-2pQ4∗Q2,JM-2pQ3∗Q3,JM-2pQ2∗Q4,JM-2pQ1∗=A~z1A~z1Φ∗αA~z1ΦβA~z1Φ∗αΦβA~x1†,ΦαΦM-2p-1βΛzΛx-1A~x1†∗=AaugSaug.It is noticed that the construction of Rnew is different from the authors’ earlier work [19], where [19] just increases the effective array aperture, while the approach utilized here increases the effective array aperture and virtual snapshots simultaneously.Similar to (6), calculate the propagator matrix Paug corresponding to Raug as follows:(10)Paug=R1augR1augH-1R1augR2augH,where R1aug and R2aug are the first K rows and the last (4M-8p-K) rows of Raug. Define the modified propagator matrix Pmaug=[IK,Paug]H and do the calculations as follows: (11)Ψα=H1Pmaug†H2Pmaug=A1augΦαA1aug-1,Ψβ=H3Pmaug†H4Pmaug=A1augΦβA1aug-1,where H1=I2⊗[0(M-2p)×(M-2p),IM-2p], H2=I2⊗[IM-2p,0(M-2p)×(M-2p)], H3=[I2M-4p,0(2M-4p)×(2M-4p)], and H4=[0(2M-4p)×(2M-4p),I2M-4p]. A1aug is the first K rows of Aaug. It is noticed that Ψ(α) and Ψ(β) have the same eigenvectors and the different eigenvalues, which is a general problem in multidimensional spectral and array signal processing [17].To estimate azimuth and elevation simultaneously, constructΨ~ using the complex eigenvalue technique [17], (12)Ψ~=Ψα+IK-1Ψα-IK-jΨβ+IK-1Ψβ-IK=A1augΦ~A1aug-1,where Φ~=diag⁡{tan⁡(β1/2)+jtan⁡(α1/2),…,tan⁡(βK/2)+jtan⁡(αK/2)}. Then applying eigenvalue decomposition to Ψ~, we can obtain Φ~. Therefore, the azimuth and elevation angles can be estimated, (13)ϕ^k=arccos⁡-β^kλ2πd,(14)θ^k=arccos⁡-α^kλ2πdsin⁡ϕ^k,where α^k=2arctan⁡[Im(Φ~(k,k))], and β^K=2arctan⁡[Re⁡(Φ~(k,k))]. Φ~(k,k) denotes the kth row and kth column element of Φ~.From (11) to (14), we know that the proposed method is able to estimate the angles simultaneously, which avoids the error accumulation effect that occurs in [13].


## 4. Algorithm Analysis

### 4.1. Limiting Performance Analysis

As shown in Section3.1, we compensate the mutual coupling effect blindly at the cost of a few auxiliary sensors. In detail, we choose the sensors with coordinates (md,0,ld) (where (l,m)∈{(1,0),…,(p-1,0),(M-p+1,0),…,(M,0),(0,0),(0,1),…,(0,p-1),(0,M-p+1),…,(0,M)}) as the auxiliary sensors and the others as the new received array data. That is to say, the received data of 2M-4p+2 sensors is used to estimate the 2D DOAs. What is more, in Section 3.2, to make (8) and (11) hold, the maximum number of sources Kmax≤M-2p and the source signal correlation factor ρ<1: that is, any two source signals cannot be coherent.Besides, the choice ofp affects the angle estimation performance greatly. Since the smaller p is, the larger array aperture and degrees of freedom can be utilized for 2D DOA estimation, which results in higher accuracy in angle estimation. However, p originates from the approximated model of the real antenna arrays with mutual coupling [5]. The larger p, the smaller error for the approximation of model in (1) to real model. Hence, for real antenna arrays, there exists an optimum p for angle estimation. However, the choice of optimum p is a difficult task and requires complicated electromagnetic simulations and real experiments, which is beyond the scope of this paper and needs further research. In particular, according to [20–22], the real mutual coupling can be reduced via various techniques and designs. Therefore, for simplicity, we will focus on the array with fixed p for the simulations in Section 5.


### 4.2. Computational Complexity Analysis

Before the algorithm computational complexity analysis, let us make some appointments as follows:(i)
As it is known that the complexity of complex multiplication is larger than that of complex add, we only calculate the complexity of the complex multiplication for analysis.(ii)
Since the left-multiplications of the selection matrices are just written to make the proposed method expressed well for paper writing and can be easily implemented without any multiplication operations, the computational complexity calculation does not include them, for instance, (4)-(5) and (9).According to Section3.2, the computational complexity of the proposed method includes construction of the propagator matrix P in (6): O{(2M-4p+2)2N+2(2M-4p+2)K2+K3+(2M-4p+2)(2M-4p+2-K)K}; construction of Q1, Q2, Q3, and Q4 in (8): O{8(M-2p)2K}; calculation of the propagator matrix Paug in (10): O{2(2M-4p)K2+K3+(4M-8p-K)(2M-4p)K}; calculation of Ψ(α), Ψ(β), and Ψ~ in (11) and (12): O{2(2M-4p)2K+2(2M-4p)K2+4K3}; and calculation of Φ~ via applying eigenvalue decomposition to Ψ~: O{K3}.For conventional 2D DOA estimation, asM≫K and N≫K, the total computational complexity of the proposed method is approximated as O{4(M-2p+1)2N+28(M-2p)2K}. While the total computational complexity of Wu et al.’s method [12] is about O{(2M+1)2N+2NθNϕ(p+ξ+1)(2M+1)2} (where Nθ and Nϕ denote the total search times within the search range of azimuth and elevation), that of Liang et al.’s method [13] is about O{4(M-4)2N+8(M-4)2K+(M-4-K)(M-4)η} (where η≫2 is a relatively large integer for known polynomial rooting algorithms [23]). Since Nθ,Nϕ≫M and η≫2, apparently, for p=2, the computational complexity of the proposed method is far less than that of the two methods.


### 4.3. Theoretical Performance Analysis

The theoretical analysis of the proposed method is based on the matrix perturbation theory [24, 25] and the perturbed data model is given as follows:(15)yn=y0n+Δyn=CAsn+wn.Lety~(n)=y~0(n)+Δy~(n), where y~0(n)=A~s(n) and Δy~(n)=[wxT(n),wzT(n)]T. Define R~=R~0+ΔR~0; then the perturbation of the covariance matrix can be expressed as (16)ΔR~0=EΔy~ny~0Hn+Ey~0nΔy~Hn+EΔy~nΔy~Hn.From (6), P=P0+ΔP, R~1=R~1,0+ΔR~1, and R~2=R~2,0+ΔR~2; we have R~2,0+ΔR~2=R~1,0+ΔR~1P0+ΔP. Neglect the second-order term ΔR~1ΔP to obtain(17)ΔP=R~1,0†ΔR~2-ΔR~1P0.Hence, Pm=Pm,0+ΔPm=Pm,0+[0K×K,ΔP]H.ForQ1, Q2, Q3, and Q4 in (8), and similar to (17), we have (18)ΔQ1=G3-Q1,0G1ΔPmG1Pm,0†,ΔQ2=G3-Q2,0G2ΔPmG2Pm,0†,ΔQ3=G4-Q1,0G1ΔPmG1Pm,0†,ΔQ4=G4-Q2,0G2ΔPmG2Pm,0†.LetRaug=AaugSaug+ΔRaug; then(19)ΔRaug=ΔQ1,JM-2pΔQ4∗ΔQ2,JM-2pΔQ3∗ΔQ3,JM-2pΔQ2∗ΔQ4,JM-2pΔQ1∗.Using the similar method in (17) to make the perturbation analysis of Paug in (10), we have(20)ΔPaug=R1,0augH†ΔR2augH-ΔR1augHPaug,0, where Paug=Paug,0+ΔPaug, R1aug=R1,0aug+ΔR1aug, and R2aug=R2,0aug+ΔR2aug.Hence,Pmaug=Pm,0aug+ΔPmaug=Pm,0aug+[0K×K,ΔPaug]H. Further, let Ψ(α)=Ψ0(α)+ΔΨ(α), Ψ(β)=Ψ0(β)+ΔΨ(β); we can obtain (21)ΔΨα=H1Pm,0aug†H2ΔPmaug-H1ΔPmaugΨ0α,ΔΨβ=H3Pm,0aug†H4ΔPmaug-H3ΔPmaugΨ0β.SubstitutingΔΨ(α) and ΔΨ(β) into (12), we use the Sherman-Morrison-Woodbury Formula [26] to obtain(22)ΔΨ~=Ψ0±αΔΨαIK-Ψ0±αΨ0-α-jΨ0±βΔΨβIK-Ψ0±βΨ0-β,where Ψ~=Ψ~0+ΔΨ~, Ψ0±(α)=Ψ0(α)+IK-1, Ψ0-(α)=Ψ0(α)-IK, Ψ0±(β)=Ψ0(β)+IK-1, and Ψ0-(β)=Ψ0(β)-IK.Calculate the first-order perturbation of the eigenvalues ofΨ~0 due to ΔΨ~,(23)Δλ~k=vkΔΨ~uk, where vk and uk are the left and right orthonormal eigenvectors of Ψ~0 associated with λ~k; that is, Ψ~0=[u1,…,uK]·diag⁡{λ~1,…,λ~K}·[v1T,…,vKT]T.Further, using the relationship thatαk=2arctan⁡(Im⁡(λ~k)), βk=2arctan⁡(Re⁡(λ~k)), we can obtain (24)Δαk=2cos2⁡αk2Im⁡Δλ~k,Δβk=2cos2⁡βk2Re⁡Δλ~k.Givenϕk=arccos⁡(-βkλ/(2πd)), θk=arccos⁡(-αkλ/(2πdsin⁡ϕk)), the perturbations of azimuth and elevation angles are (25)Δϕk=∂ϕk∂βkΔβk=CϕβΔβk,Δθk=∂θk∂αk,∂θk∂ϕkΔαkΔϕk=CθαΔαk+CθϕΔϕk,where Cϕβ=λ/(2πd)/1-βkλ/(2πd)2, Cθα=λ/2πdsin⁡ϕk/1-αkλ/2πdsin⁡ϕk2, Cθϕ=-αkλcos⁡ϕk/2πdsin2⁡ϕk/1-αkλ/2πdsin2⁡ϕk2.Therefore, the mean-squared error of estimated DOAs is(26)var⁡Δϕk=Cϕβ2var⁡Δβk=2cos2⁡βk2Cϕβ2var⁡Re⁡Δλ~k,var⁡Δθk=Cθα2var⁡Δαk+Cθϕ2var⁡Δϕk=2cos2⁡αk2Cθα2var⁡Im⁡Δλ~k+2cos2⁡βk2Cϕβ2Cθϕ2var⁡Re⁡Δλ~k.


## 4.1. Limiting Performance Analysis

As shown in Section3.1, we compensate the mutual coupling effect blindly at the cost of a few auxiliary sensors. In detail, we choose the sensors with coordinates (md,0,ld) (where (l,m)∈{(1,0),…,(p-1,0),(M-p+1,0),…,(M,0),(0,0),(0,1),…,(0,p-1),(0,M-p+1),…,(0,M)}) as the auxiliary sensors and the others as the new received array data. That is to say, the received data of 2M-4p+2 sensors is used to estimate the 2D DOAs. What is more, in Section 3.2, to make (8) and (11) hold, the maximum number of sources Kmax≤M-2p and the source signal correlation factor ρ<1: that is, any two source signals cannot be coherent.Besides, the choice ofp affects the angle estimation performance greatly. Since the smaller p is, the larger array aperture and degrees of freedom can be utilized for 2D DOA estimation, which results in higher accuracy in angle estimation. However, p originates from the approximated model of the real antenna arrays with mutual coupling [5]. The larger p, the smaller error for the approximation of model in (1) to real model. Hence, for real antenna arrays, there exists an optimum p for angle estimation. However, the choice of optimum p is a difficult task and requires complicated electromagnetic simulations and real experiments, which is beyond the scope of this paper and needs further research. In particular, according to [20–22], the real mutual coupling can be reduced via various techniques and designs. Therefore, for simplicity, we will focus on the array with fixed p for the simulations in Section 5.


## 4.2. Computational Complexity Analysis

Before the algorithm computational complexity analysis, let us make some appointments as follows:(i)
As it is known that the complexity of complex multiplication is larger than that of complex add, we only calculate the complexity of the complex multiplication for analysis.(ii)
Since the left-multiplications of the selection matrices are just written to make the proposed method expressed well for paper writing and can be easily implemented without any multiplication operations, the computational complexity calculation does not include them, for instance, (4)-(5) and (9).According to Section3.2, the computational complexity of the proposed method includes construction of the propagator matrix P in (6): O{(2M-4p+2)2N+2(2M-4p+2)K2+K3+(2M-4p+2)(2M-4p+2-K)K}; construction of Q1, Q2, Q3, and Q4 in (8): O{8(M-2p)2K}; calculation of the propagator matrix Paug in (10): O{2(2M-4p)K2+K3+(4M-8p-K)(2M-4p)K}; calculation of Ψ(α), Ψ(β), and Ψ~ in (11) and (12): O{2(2M-4p)2K+2(2M-4p)K2+4K3}; and calculation of Φ~ via applying eigenvalue decomposition to Ψ~: O{K3}.For conventional 2D DOA estimation, asM≫K and N≫K, the total computational complexity of the proposed method is approximated as O{4(M-2p+1)2N+28(M-2p)2K}. While the total computational complexity of Wu et al.’s method [12] is about O{(2M+1)2N+2NθNϕ(p+ξ+1)(2M+1)2} (where Nθ and Nϕ denote the total search times within the search range of azimuth and elevation), that of Liang et al.’s method [13] is about O{4(M-4)2N+8(M-4)2K+(M-4-K)(M-4)η} (where η≫2 is a relatively large integer for known polynomial rooting algorithms [23]). Since Nθ,Nϕ≫M and η≫2, apparently, for p=2, the computational complexity of the proposed method is far less than that of the two methods.


## 4.3. Theoretical Performance Analysis

The theoretical analysis of the proposed method is based on the matrix perturbation theory [24, 25] and the perturbed data model is given as follows:(15)yn=y0n+Δyn=CAsn+wn.Lety~(n)=y~0(n)+Δy~(n), where y~0(n)=A~s(n) and Δy~(n)=[wxT(n),wzT(n)]T. Define R~=R~0+ΔR~0; then the perturbation of the covariance matrix can be expressed as (16)ΔR~0=EΔy~ny~0Hn+Ey~0nΔy~Hn+EΔy~nΔy~Hn.From (6), P=P0+ΔP, R~1=R~1,0+ΔR~1, and R~2=R~2,0+ΔR~2; we have R~2,0+ΔR~2=R~1,0+ΔR~1P0+ΔP. Neglect the second-order term ΔR~1ΔP to obtain(17)ΔP=R~1,0†ΔR~2-ΔR~1P0.Hence, Pm=Pm,0+ΔPm=Pm,0+[0K×K,ΔP]H.ForQ1, Q2, Q3, and Q4 in (8), and similar to (17), we have (18)ΔQ1=G3-Q1,0G1ΔPmG1Pm,0†,ΔQ2=G3-Q2,0G2ΔPmG2Pm,0†,ΔQ3=G4-Q1,0G1ΔPmG1Pm,0†,ΔQ4=G4-Q2,0G2ΔPmG2Pm,0†.LetRaug=AaugSaug+ΔRaug; then(19)ΔRaug=ΔQ1,JM-2pΔQ4∗ΔQ2,JM-2pΔQ3∗ΔQ3,JM-2pΔQ2∗ΔQ4,JM-2pΔQ1∗.Using the similar method in (17) to make the perturbation analysis of Paug in (10), we have(20)ΔPaug=R1,0augH†ΔR2augH-ΔR1augHPaug,0, where Paug=Paug,0+ΔPaug, R1aug=R1,0aug+ΔR1aug, and R2aug=R2,0aug+ΔR2aug.Hence,Pmaug=Pm,0aug+ΔPmaug=Pm,0aug+[0K×K,ΔPaug]H. Further, let Ψ(α)=Ψ0(α)+ΔΨ(α), Ψ(β)=Ψ0(β)+ΔΨ(β); we can obtain (21)ΔΨα=H1Pm,0aug†H2ΔPmaug-H1ΔPmaugΨ0α,ΔΨβ=H3Pm,0aug†H4ΔPmaug-H3ΔPmaugΨ0β.SubstitutingΔΨ(α) and ΔΨ(β) into (12), we use the Sherman-Morrison-Woodbury Formula [26] to obtain(22)ΔΨ~=Ψ0±αΔΨαIK-Ψ0±αΨ0-α-jΨ0±βΔΨβIK-Ψ0±βΨ0-β,where Ψ~=Ψ~0+ΔΨ~, Ψ0±(α)=Ψ0(α)+IK-1, Ψ0-(α)=Ψ0(α)-IK, Ψ0±(β)=Ψ0(β)+IK-1, and Ψ0-(β)=Ψ0(β)-IK.Calculate the first-order perturbation of the eigenvalues ofΨ~0 due to ΔΨ~,(23)Δλ~k=vkΔΨ~uk, where vk and uk are the left and right orthonormal eigenvectors of Ψ~0 associated with λ~k; that is, Ψ~0=[u1,…,uK]·diag⁡{λ~1,…,λ~K}·[v1T,…,vKT]T.Further, using the relationship thatαk=2arctan⁡(Im⁡(λ~k)), βk=2arctan⁡(Re⁡(λ~k)), we can obtain (24)Δαk=2cos2⁡αk2Im⁡Δλ~k,Δβk=2cos2⁡βk2Re⁡Δλ~k.Givenϕk=arccos⁡(-βkλ/(2πd)), θk=arccos⁡(-αkλ/(2πdsin⁡ϕk)), the perturbations of azimuth and elevation angles are (25)Δϕk=∂ϕk∂βkΔβk=CϕβΔβk,Δθk=∂θk∂αk,∂θk∂ϕkΔαkΔϕk=CθαΔαk+CθϕΔϕk,where Cϕβ=λ/(2πd)/1-βkλ/(2πd)2, Cθα=λ/2πdsin⁡ϕk/1-αkλ/2πdsin⁡ϕk2, Cθϕ=-αkλcos⁡ϕk/2πdsin2⁡ϕk/1-αkλ/2πdsin2⁡ϕk2.Therefore, the mean-squared error of estimated DOAs is(26)var⁡Δϕk=Cϕβ2var⁡Δβk=2cos2⁡βk2Cϕβ2var⁡Re⁡Δλ~k,var⁡Δθk=Cθα2var⁡Δαk+Cθϕ2var⁡Δϕk=2cos2⁡αk2Cθα2var⁡Im⁡Δλ~k+2cos2⁡βk2Cϕβ2Cθϕ2var⁡Re⁡Δλ~k.


## 5. Simulation Results

In this section, the 2D DOA estimation performance of the proposed method is verified through simulation experiments in comparison with its theoretical analysis, Wu et al.’s method [12], Liang et al.’s method [13], and the Cramer-Rao bound (CRB) [27]. Without loss of generality, we assume that d=λ/2, all source signals have the same power σs2, and the signal to noise ratio (SNR) is defined as 10log10(σs2/σw2). The search ranges of azimuth and elevation angles for Wu et al.’s method are [0∘,180∘] with interval 0.01∘. In order to measure the overall performance of estimating the azimuth and elevation angles, define the root-mean-square-error (RMSE) of the 2D DOA estimations from V Monte Carlo trials as(27)RMSE=1KV∑v=1V∑k=1Kθ^kv-θk2+ϕ^kv-ϕk2.In addition, to evaluate the performance of the proposed method under different correlation factorρ between the incident signals, we generate two source signals s1(n) and s2(n) with determined ρ as follows:(28)s2n=s1n+1-ρ2s3n, where s1(n) is uncorrelated with s3(n). The magnitude of ρ varies from 0 to 1. Specially, for ρ=0, we usually call two source signals uncorrelated or independent. While for ρ=1, two source signals are called coherent.Since Wu et al.’s method and Liang et al.’s method can also be used forp≥2 with a few modifications, for simplicity, we choose p=2 for the following simulation examples. The mutual coupling coefficients are c1=0.4528-0.3764j, c2=0.1085+0.0987j, and c3=0.1828-0.3076j associated with distances d, 2d, and 2d.Example 1 (estimation performance versus SNR).
In the first example, the performance of the proposed method with respect to the SNR under different source signal correlations is investigated. Two sources are placed at(115∘,65∘) and (95∘,80∘). The number of subarray elements and the number of snapshots are fixed at M=8 and N=500. SNR varies from −10 dB to 30 dB with interval −2.5 dB. 500 Monte Carlo trials have been conducted for every fixed SNR. Figures 2 and 3 show the RMSEs of 2D DOA estimations for ρ=0 and ρ=0.6, that is, the uncorrelated and high correlated cases.
As we can see, the proposed method yields better 2D DOA estimation performance than Liang et al.’s method for SNR≥ −5 dB and its experimental performance coincides with the theoretical analysis well for SNR ≥ 0 dB, which verifies the correctness of the proposed method and its theoretical analysis. Specially, for ρ=0.6, the increase of angle estimation accuracy with the proposed method to Liang et al.’s method is greater than that of ρ=0, which proves that the proposed method can decrease the correlation between highly correlated source signals. Besides, from Figure 2, since the proposed method estimates angles simultaneously, it has better angle estimation performance than Liang et al.’s method which divides the angle estimations into two sequential steps and results in the error accumulation effect.
However, to estimate 2D DOAs with less time, the proposed method utilizes some array elements as auxiliary sensors and some of the array aperture is lost; hence, the angle estimation performance of proposed method is worse than that of Wu et al.’s method for both uncorrelated and correlated cases.Figure 2
RMSE versus SNR,M = 8, N = 500, ρ = 0.Figure 3
RMSE versus SNR,M = 8, N = 500, ρ = 0.6.Example 2 (estimation performance versus the number of snapshots).
In this example, we concentrate on the 2D DOA estimation performance of the proposed method versus the number of snapshots with differentρ’s. All simulation conditions are the same as Example 1 except that SNR is fixed at 10 dB and N ranges from 100 to 10000. The angle estimation results are shown in Figures 4 and 5.
We can know that the angle estimation performance of the proposed method outperforms that of Liang et al.’s method for differentN’s under SNR = 10 dB case. What is more, for ρ=0.6, the proposed method works well, while Liang et al.’s method does not, which verifies the decorrelation performance of the proposed method. The theoretical performance is very close to the experimental results under two different ρ’s. For ρ=0.6, the experimental results are always slightly higher than the theoretical results because the proposed method can only decrease the correlation effect not eliminate it.Figure 4
RMSE versusN, M = 8, SNR = 10 dB, ρ = 0.Figure 5
RMSE versusN, M = 8, SNR = 10 dB, ρ = 0.6.Similar to Example1, from Figures 4 and 5, we know that the 2D DOA estimation performance of the proposed method is still worse than that of Wu et al.’s method owing to the loss in array aperture.Example 3 (estimation performance versus the correlation factor).
In this example, we consider the angle estimation performance against the correlation factorρ and the simulation settings are similar to that of Example 1 except that SNR is equal to 10 dB, and the correlation factor ρ varies from 0 (uncorrelated) to 1 (coherent). The results are shown in Figure 6.
It is noticed that the 2D DOA estimation performance of the proposed method and Wu et al.’s method are both more robust to the source signal correlation than that of Liang et al.’s method for0≤ρ≤0.9. Importantly, similar to Wu et al.’s method, the proposed method can work well with extremely high correlation factor: that is, ρ=0.9. In addition, both of the three methods cannot work at coherent case: that is, ρ=1, which is an open problem for research in the future.Figure 6
RMSE versus correlation factor,M = 8, SNR = 10 dB, N = 500.Example 4 (runtime versus the number of subarray elements).
In this example, we investigate the runtime of the proposed method with respect toM. All simulation conditions are the same as Example 1 except that SNR = 10 dB and M ranges from 8 to 96 with interval 8. The implementation of Liang et al.’s method utilizes the Symbolic Toolbox,sym2poly, androots function in MATLAB. The measurements of cpu runtime for these three methods are shown in Figure 7, and Figure 8 gives the ratios of the runtime of the proposed method to those of the other two methods.Figure 7
Runtime versus the number of subarray elements, SNR = 10 dB,N = 500, ρ = 0.6.Figure 8
Ratio versus the number of subarray elements, SNR = 10 dB,N = 500, ρ = 0.6.Apparently, the proposed method is much more computationally efficient than Liang et al.’s method and Wu et al.’s method. What is more, the runtime of the proposed method is less than10-6 of Wu et al.’s and 10-2 of Liang et al.’s, which indicates that the proposed method is a significant improvement of Wu et al.’s method and Liang et al.’s method in computation time. Besides, as shown in Figure 8, with the increase in the number of subarray elements, the ratio of the cpu runtime of the proposed method to that of Liang et al.’s method decreases greatly. The reason is that their computational complexity expressions, neglecting the lower order terms and the constant coefficients including N, can be expressed approximately as O{M2} and O{Mη+1} (η≫2). Therefore, when M increases, the ratio decreases greatly as shown in Figure 8.Overall, from Figures2–8, the proposed method can estimate 2D DOAs against mutual coupling under highly correlated cases, which is similar to the ability of Wu et al.’s method. Importantly, the proposed method has the lowest computational complexity among the three methods. Besides, the angle estimation performance of the proposed method is outperformed by Wu et al.’s method in all cases, which results from the utilization of RARE technique guaranteeing less loss in array aperture and degrees of freedom than the auxiliary sensor method used by the proposed method.


## 6. Conclusions

In this paper, we have developed a computationally efficient 2D DOA estimation method for L-shaped array with unknown mutual coupling and its explicit expressions of theoretical performance for angle estimations have also been derived. Simulation results and theoretical analysis have verified the effectiveness of the proposed method and shown that the proposed method can handle highly correlated source signals with lower computational complexity and also has fine 2D DOA estimation performance against coupling effect. However, the proposed method has assumed that some elements are auxiliaries which is not consistent with the cases in real scenarios and results in larger estimation errors for deployment cases. Hence, further research on real application test is needed.



---

*Source: 5454719-2018-01-31.xml*