Efficient Two-Dimensional Direction Finding Algorithm for Rectilinear Sources Under Unknown Mutual Coupling.

Digital communication signals in wireless systems may possess noncircularity, which can be used to enhance the degrees of freedom for direction-of-arrival (DOA) estimation in sensor array signal processing. On the other hand, the electromagnetic characteristics between sensors in uniform rectangular arrays (URAs), such as mutual coupling, may significantly deteriorate the estimation performance. To deal with this problem, a robust real-valued estimator for rectilinear sources was developed to alleviate unknown mutual coupling in URAs. An augmented covariance matrix was built up by extracting the real and imaginary parts of observations containing the circularity and noncircularity of signals. Then, the actual steering vector considering mutual coupling was reparameterized to make the rank reduction (RARE) property available. To reduce the computational complexity of two-dimensional (2D) spectral search, we individually estimated y-axis and x-axis direction-cosines in two stages following the principle of RARE. Finally, azimuth and elevation angle estimates were determined from the corresponding direction-cosines respectively. Compared with existing solutions, the proposed method is more computationally efficient, involving real-valued operations and decoupled 2D spectral searches into twice those of one-dimensional searches. Simulation results verified that the proposed method provides satisfactory estimation performance that is robust to unknown mutual coupling and close to the counterparts based on 2D spectral searches, but at the cost of much fewer calculations.


Introduction
Uniform rectangular array (URA), due to its versatility in providing both azimuthal and elevational coverage simultaneously, has attracted much attention in many applications, such as radar, sonar, and wireless communication [1]. Two-dimensional (2D) direction finding with URA is thus of practical use in many situations. Based on the subspace principle, the classical 1D super resolution algorithms such as multiple signal classification (MUSIC) [2] and compressive sensing [3] can be readily generalized to URA. Moreover, to further reduce the computational complexity or increase the accuracy of the direction finding results, some effective schemes have been proposed [4][5][6][7][8]. In [4], the authors proposed a 2D unitary ESPRIT algorithm for URA, which can be realized in both element space and beamspace. Also, its main procedure is formulated in terms of real-valued operation and is therefore computationally efficient. In [5], a dual-size spatial invariance array structure was proposed. Based on the ESPRIT theory, it utilizes unambiguous but low-accuracy estimates to disambiguate high-accuracy but periodically ambiguous direction cosine estimates. In [6], a 2D matrix pencil method is developed to convert the complex data matrix to a real one, which can reduce the complexity of the computation significantly.
symmetric Toeplitz property of the mutual coupling matrix and transformed this matrix into gain and phase perturbations of the array [21], where the perturbations of the middle subarray elements are identical. Consequently, the observation data from all the array elements can be fully utilized for DOA estimation via simple 1D RARE method. Furthermore, a fourth order cumulant (FOC) method was also proposed [22], which performs better especially for strong mutual coupling effects. In [23][24][25], the authors utilized sparse reconstruction methods to improve the DOA estimation performance under unknown mutual coupling. However, these methods are computationally demanding due to the l 1 -norm minimization problem. In [26], the authors presented a method for the localization of mixed sources using the noncircular information of the impinging signals under unknown mutual coupling. In [27], an ESPRIT-based algorithm is proposed for the ODA estimation of strictly noncircular sources with unknown mutual coupling, which has large virtual array aperture and low computational load. According to the latest work by B. Friedlander [28], the mutual coupling matrix is a direction dependent matrix, and not a constant matrix as is commonly assumed, which may impact some of the work presented in the DOA estimation literature. However, few of the works on DOA estimation have taken into account the direction-dependent mutual coupling effect. In [29], a decoupled 2D DOA estimation algorithm is proposed for uniform circular arrays. Yet, it assumes that the mutual coupling effect is only elevation-dependent. Also, it is an offline technique, for which is costly to build a calibration setup. In [30], a unified transformation method is proposed for the DOA estimation under direction-dependent mutual coupling. This method can be used in many array structures, but it requires convex optimization and iteration processes, which is computationally inefficient. In [31], the authors proposed a RARE-based low complexity DOA estimation method to solve the direction-dependent mutual coupling problem. Moreover, coprime sparse arrays are proposed to reduce the mutual coupling effect between the elements [32][33][34][35].
For the 2D array calibration problem, several effective methods have been developed [36][37][38][39][40]. In [36], the authors proposed a blind calibration scheme for uniform circular array based on azimuthal symmetric structures in the mutual coupling matrix. However, it only provides the azimuth angle estimates of the radiating sources, while the elevation angles are assumed to be zeros. In [37], a robust method was proposed to calibrate the mutual coupling and channel gain/phase inconsistency of uniform circular array. Yet, it requires calibration sources and 2D spectral searches. In [38], the auxiliary sensor-based algorithm is generalized for URA and uniform hexagon arrays, where the elements at the boundary of the array are set as the auxiliary ones. In order to mitigate the array aperture loss problem in the auxiliary sensor-based method, Wu et al. proposes a RARE-based auto-calibration method for URA [40], which can fully utilize the array observations. However, all these planar array calibration methods require exhaustive 2D spectral searches for the DOA estimation, which is computationally demanding.
To the extent of our knowledge, few works have been devoted to the 2D direction finding algorithms for rectilinear sources in the presence of unknown mutual coupling effects. Despite the various advantages that were manifested in the approaches cited above, there are some fundamentally inherent problems to be overcome in view of the previous analyses: (1) all of the above-mentioned 2D DOA estimation works only use the covariance matrix for direction finding and mutual coupling calibration, the rectilinearity of the sources are not fully utilized. (2) These works require 2D spectral searches for the azimuth and elevation angle estimations, which is computationally demanding.
(3) Generally, subspace-based direction finding algorithms require complex-valued operations for the calculation of covariance matrix, eigenvalue decomposition (EVD), and spectral searching. Many works have been devoted to converting the complex-valued covariance into a real-valued one. The computational complexity of EVD and spectral searches for this kind of transformation can be decreased by a factor of four. Yet, the transformation itself and the computation of the covariance still require additional complex-valued operations.
Therefore, aiming at addressing these problems, we developed a computationally efficient 2D direction finding algorithm for rectilinear signals under unknown mutual coupling. Based on the rectilinearity of the signal, the received data is transformed into the real-valued domain. Then the structure of the virtual steering vector is utilized to decouple the DOA parameter from other nuisance parameters. According to the RARE principle, the y-axis direction-cosine estimates are first obtained via a 1D spectral search. Consequently, the x-axis direction-cosines are then generated through another RARE estimator, which also only requires a 1D spectral search with the y-axis direction-cosine estimates.
The main contributions of this paper are as follows: (1) The proposed DOA estimation method avoids the two-dimensional spectral search and effectively reduces the computational complexity; (2) Only real-valued operations are needed, and its computation load is reduced to 1/4 of the complex value operations; (3) The estimation accuracy is improved by the efficient use of the source rectilinearity.
The rest of this paper is organized as follows. In Section 2, we present the signal model for the rectilinear sources of a URA under unknown mutual coupling. In Section 3, derivation of the proposed real-valued two stage 2D DOA estimation algorithm is detailed. Then we discuss the computational complexity and the maximum number of resolvable sources in Section 4. In Section 5, several numerical simulations are carried out to validate the performance of the proposed algorithm. In Section 6, the paper ends with conclusions. Finally, in the Appendix A, we have derived the Cramér-Rao Bound (CRB) for the problem.
Throughout the paper, we make use of the following notations as listed in Table 1.

Signal Model
As shown in Figure 1, consider a URA composed of M = M x × M y omnidirectional elements, with the inter element spacing along the x-axis and the y-axis being d x and d y , respectively. Suppose that there are K far-field narrowband sources radiating onto the URA from directions (θ 1 , φ 1 ), . . . ,(θ K , φ K ), where θ k and φ k are the azimuth and elevation angles of the k-th source. Under the ideal scenario, the received data of each element only depends on the incident wavefield. Accordingly, the ideal array received data vector without mutual coupling can be modeled as where , , However, in practical direction finding systems, a mismatch may be introduced by the electromagnetic interaction between the array elements. Therefore, taking into account the mutual coupling effect between the array elements, we can modify the received data model as CAs n (5) where C is the mutual coupling matrix, with its elements being mutual coupling coefficients (MCCs). As illustrated in [19], the mutual coupling effect between the elements is inversely proportional to their distance. Thus, the MCCs between the neighboring sensors are approximately equal. When two elements are far apart, the mutual coupling effect is negligible, and the coefficients approach zero. As in [30], the general mutual coupling model for URA is developed by the following block matrix Under the ideal scenario, the received data of each element only depends on the incident wavefield. Accordingly, the ideal array received data vector without mutual coupling can be modeled as T denotes the array noise vector, and A = [a(θ 1 , φ 1 ), · · · , a(θ K , φ K )] represents the ideal steering matrix, with its steering vectors being a(θ k , φ k ) = a y (ψ y,k ) ⊗ a x (ψ x,k ), k = 1, · · · , K where a x (ψ x,k ) = 1, e jψ x,k , · · · , e j(M x −1)ψ x,k T (3) a y (ψ y,k ) = 1, e jψ y,k , · · · , e j(M y −1)ψ y,k T Herein, ψ x,k and ψ y,k are the electrical angles of the kth source with respect to the x-axis and the y-axis, with ψ x,k = 2πd x cos θ k sin φ k /λ and ψ y,k = 2πd y sin θ k sin φ k /λ where λ is the carrier wavelength.
However, in practical direction finding systems, a mismatch may be introduced by the electromagnetic interaction between the array elements. Therefore, taking into account the mutual coupling effect between the array elements, we can modify the received data model as where C is the M x M y × M x M y mutual coupling matrix, with its elements being mutual coupling coefficients (MCCs). As illustrated in [19], the mutual coupling effect between the elements is inversely proportional to their distance. Thus, the MCCs between the neighboring sensors are approximately equal. When two elements are far apart, the mutual coupling effect is negligible, and the coefficients approach zero. As in [30], the general mutual coupling model for URA is developed by the following block matrix Sensors 2020, 20, 1914 and the M x × M x submatrix C p (p = 1, . . . ,P) is expressed as [27] The complex mutual coupling coefficient c m,n = ρ m,n e jξ m,n symbolizes the electromagnetic interaction between a reference sensor defined as the origin of a Cartesian coordinate system and the sensor located at (±md x , ±nd y ) under the corresponding coordinate [30]. Note that c 0,0 is especially defined as 1. Therefore, the mutual coupling matrix C can be further expressed as and J p is a M y × M y selection matrix, whose elements are Moreover, for rectilinear sources, the signal vector is of the following form [10] s(t) = Φs 0 (t) (10) where s 0 (t) = [s 0,1 (t), · · · , s 0,K (t)] T , with s 0,k (t) being the real-valued zero-phase version of the k-th source signal, Φ = diag{e jϕ 1 , · · · , e jϕ K } contains the noncircular phase shifts of the sources [9]. Owing to the above analysis, the received data vector of the planar array can be formulated as With the received data vector x(t), the main problem to be addressed in this paper is to obtain the 2D DOA estimates of rectilinear signals under unknown mutual coupling effects. Unlike the traditional methods which require multidimensional spectral searches and omit the rectilinearity of the signals, a novel algorithm is investigated in the next section for reducing the computational complexity and promoting the estimation performance.

Proposed Real-Valued 2D DOA Estimation Algorithm
In order to generate accurate 2D DOA estimates with low computational complexity, we creatively utilize the rectilinearity of signals and factorize the steering vector with respect to the DOA and other nuisance parameters. Our technique proceeds through the following stages: (1) based on the rectilinearity of the signal, the received data is transformed into the real-valued domain; (2) the structure of the virtual steering vector is utilized to decouple the DOA parameter from other nuisance parameters; (3) according to the RARE principle, the y-axis direction-cosine estimates can be obtained via a 1D spectral search; (4) with the y-axis direction-cosine estimates, the x-axis direction-cosine can be found through another RARE estimator, which also only requires a 1D spectral search. The following subsections give the details of the proposed algorithm.

Extended Real-Valued Signal Model
It is well known that by transforming complex-valued data into the real-valued domain, the computational cost can be effectively reduced by a factor of four. Therefore, it inspired us to look for a more efficient algorithm to reduce the complex-valued computational loads.
Based on the sources' rectilinearity, the transformation can be realized by extracting the real and imaginary parts of x(t), respectively.
From Equation (11), x R (t) and x I (t) can be further expressed as where are the real-valued manifold matrices of x R (t) and x I (t). n R (t) and n I (t) symbolize the real and imaginary components of sensor noise vector n(t), respectively. Accordingly, the real-valued virtual steering vectors in A R is of the following form Based on Equations (2) and (8), we can obtain the following equation from the Kronecker product property (A ⊗ B)(C ⊗ D) = (AC) ⊗ (BD) [41] Therefore, the real-valued virtual steering vector v R (θ k , φ k , ϕ k , C) can be further written as where a y,c (ψ y,k ) = 1, cos(ψ y,k ), · · · , cos (M y − 1)ψ y,k T Sensors 2020, 20, 1914 are two scalar parameters determined by the x-axis direction-cosine, rectilinearity phases, and MCCs. For p = 1, 2, · · · , P − 1, we have Similarly, the real-valued virtual steering vectors of A I can be defined as Based on the above real-valued signal model, we can concatenate the real part and the imaginary part to generate an augmented data vector in the real-valued domain. where is the augmented manifold matrix, and the corresponding real-valued steering vector of the k-th source is defined as Sensors 2020, 20, 1914 9 of 22 Therefore, the covariance matrix of the augmented real-valued data matrix y(t) is of the following form where R S 0 = E{s 0 (t)s T 0 (t)}, I 2M denotes a 2M × 2M identity matrix. In practical applications, R y can be approximated from the snapshots generated at L distinct time instantŝ Through the eigenvalue decomposition of R y , we can obtain that Herein, Λ s is a diagonal matrix comprised of the largest K eigenvalues of R y . E s is the signal subspace matrix which contains the eigenvectors corresponding to the largest K eigenvalues. E n is the noise subspace matrix, which is composed of the remaining 2M-K eigenvectors of R y .
Based on the subspace principle that the noise subspace is orthogonal to the signal subspace, we can construct the following MUSIC spectrum function for estimating the azimuth and elevation angles of the sources.
Unfortunately, this multidimensional spectrum function presented above requires a (P + 2)dimension joint spectral search for 2D DOAs, MCCs, and noncircular phases, which is computationally prohibitive. Therefore, to decrease the complexity of the algorithm, we first factorize the augmented steering vector a E (θ k , φ k , ϕ k , C) with respect to the y-axis direction-cosine and present a novel real-valued RARE estimator in the next subsection, which only requires a simple 1D spectral search procedure.

First Stage Rank Reduction DOA Estimation
In this stage, the real-valued steering vector a E in Equation (32) is first factorized with respect to the y-axis direction-cosine and other nuisance parameters. From Equations (17) and (30), the mutual coupling effect in the two virtual steering vectors can be deemed as angle-dependent gain and phase perturbations. However, it is clear that the central elements on the main diagonal of both Γ c,p (ψ x ) and Γ s,p (ψ x ) comprise a vector of all ones, which implies that the middle subarray can be handled as an unperturbed one. Therefore, we can obtain the following equations Similarly, for the y-axis direction-cosine related steering vector, we have that where T y,c (ψ y ) and T y,s (ψ y ) are M y × (2P − 1) real matrices, and α y,s are (2P − 1) × 1 real vectors, which can be expressed as For simplicity, we also denote T y,c (ψ y ) as T y,c in the following derivation. Based on the relationships derived above, we can factorize the virtual steering vectors v R and v I in Equations (17) and (30) as where α x,c , α x,s , α y,c , and α y,s are defined as For simplicity, we also denote , ( ) in the following derivation. Based on the relationships derived above, we can factorize the virtual steering vectors R v and I v in Equations (17) and (30) as where    T  I I  T  T  I I  T  T  I I  T  T  I I  T  α  α  α  α  α  α  α  α  T  I I  T  T  I I  T  T  I I  T  T  I I  T   (53) where Ix is a Mx × Mx identity matrix and Iy is a (2P-1) × (2P-1) identity matrix. By factorizing the first matrix, we can reorganize E a in the following form  T  I  T  I I  T  I  T  α  α  α  α  T  I  T  I  I  T  I  T  α where I x is a M x × M x identity matrix and I y is a (2P − 1) × (2P − 1) identity matrix. By factorizing the first matrix, we can reorganize a E in the following form  (17) and (30) as where    T  I I  T  T  I I  T  T  I I  T  T  I I  T  α  α  α  α  α  α  α  α  T  I I  T  T  I I  T  T  I I  T  T  I I  T   (53) where Ix is a Mx × Mx identity matrix and Iy is a (2P-1) × (2P-1) identity matrix. By factorizing the first matrix, we can reorganize E a in the following form , , ,  T  I  T  I I  T  I  T  α  α  α  α  T  I  T  I  I  T  I  T  α where T 1 (ψ y ) is a 2M x M y × 2M x (2P − 1) matrix which only depends on the y-axis direction-cosine, and α 1 (ψ x , ψ y , ϕ, C) is a 2M x (2P − 1) × 1 vector which depends on the x-axis and y-axis direction-cosines, MCCs, and noncircular phases.

Substituting the above equation into Equation (36), we have
It is obvious that , Q 1 (ψ y ) will generally be of full column rank. Note that α 1 is presumed to be a non-zero vector, thus Equation (55) holds only when Q 1 (ψ y ) rank drops, which implies that Q 1 (ψ y ) becomes a rank-deficient matrix when ψ y = ψ y,k , k = 1, · · · , K. Therefore, based on the principle of RARE, the y-axis direction-cosines can be estimated from the highest peaks of the following 1D spectral search function: Since the y-axis direction-cosine ψ y has been decoupled with other nuisance parameters in this stage, we can simplify the previous multidimensional optimization problem in Equation (36) into a 1D spectral search problem, which is computationally more efficient. Additionally, the operations involved in this stage are all real-valued, which means the computational load can be reduced by a factor of four compared with its complex-valued counterparts.
Our next objective is to estimate the values of x-axis direction-cosines and correctly pair them with the corresponding y-axis direction-cosines. Therefore, in the following subsection, we will present another RARE function which can generate the estimation of x-axis direction-cosines ψ x,k , k = 1, · · · , K without any additional matching procedures.

Second Stage Rank Reduction DOA Estimation
In the second stage, the real-valued steering vector a E in Equation (32) is factorized in another form with respect to the x-axis direction-cosine and other nuisance parameters. Since the y-axis direction-cosine ψ y has already been obtained in the previous stage, we can calculate the values of P p=1 J p a y,c (ψ y,k ) and P p=1 J p a y,s (ψ y,k ) . Also we denote that β y,c = P p=1 J p a y,c (ψ y,k ) , β y,s = P p=1 J p a y,s (ψ y,k ) Therefore, the virtual steering vectors v R and v I in Equations (17) and (30) can be also factorized as Utilizing the Kronecker product property β ⊗ (T · α) = (β ⊗ T) · α, we have the following expression v R (θ k , φ k , ϕ k , C) = β y,c ⊗ T x,c α x,c − β y,c ⊗ T x,s α x,s − β y,s ⊗ T x,s α x,c − β y,s ⊗ T x,c α x,s v I (θ k , φ k , ϕ k , C) = β y,c ⊗ T x,c α x,s + β y,c ⊗ T x,s α x,c − β y,s ⊗ T x,s α x,s + β y,s ⊗ T x,c α x,c where I β is a M y × M y identity matrix. Consequently, the extended real-valued steering vector a E can be written in the following form Utilizing the Kronecker product property ( ) ( ) ⊗ ⋅ = ⊗ ⋅ β T α β T α , we have the following matrix which only depends on the x-axis direction-cosine, and 2 ( , , ) x ψ ϕ α C is a 2(2 1) P − vector which depends on the x-axis direction-cosines, MCCs, and noncircular phases.
where T 2 (ψ x ) is a 2M x M y × 2(2P − 1) matrix which only depends on the x-axis direction-cosine, and α 2 (ψ x , ϕ, C) is a 2(2P − 1) vector which depends on the x-axis direction-cosines, MCCs, and noncircular phases. Again, according to the RARE principle, we substitute the above equation into Equation (36) is of full column rank in general. As both α x,c and α x,s are supposed to be non-zero vectors, α 2 must be non-zero as well. One can deduce that (64) holds only if the rank of Q 2 (ψ x ) drops, that is, rank deficiency occurs in Q 2 (ψ x ) provided that ψ x = ψ x,k , k = 1, · · · , K. By the principle of RARE, the x-axis direction-cosines can be obtained by searching for the highest peaks of the following spatial spectrum function: Since the x-axis direction-cosine ψ x has been decoupled with other nuisance parameters in the second stage, the multidimensional optimization problem can be reduced to a computationally efficient 1D spectral searching problem. Moreover, the computations involved are all real-valued, which means the computational load can be reduced by a factor of four compared with its complex-valued counterparts. Also, due to the orthogonality between the signal subspace and noise subspace, the 2D DOA parameters are automatically paired without any additional matching procedures. Therefore, the 2D DOA estimates can be obtained from the following equations:

Procedures of the Proposed Algorithm
The steps of the proposed algorithm can be summarized as follows.

Discussion
In this section, we compare the performance of the proposed algorithm with respect to two state-of-art 2D DOA estimation algorithms considering mutual coupling effects: the auxiliary sensor-based method (AUX) in [39] and the 2D rank reduction-based method (2D-RARE) in [40].

Computational Complexity
In this subsection, the computational complexity of the proposed algorithm is evaluated with respect to AUX and 2D-RARE. Note that the computational complexity is investigated in terms of the major multiplications in constructing the statistical matrices, EVDs, and spectral searches. Suppose that an M x × M y URA is utilized and the grid size for spectral searching is ∆θ. From [29], it can be observed that AUX constructs a (M x − 2P + 2)(M y − 2P + 2) × (M x − 2P + 2)(M y − 2P + 2) covariance matrix, takes the EVD of this matrix, and performs a 2D spectral search for DOA estimation, which require a computational complexity of O(4(M respectively. 2D-RARE constructs a M x M y × M x M y covariance matrix, performs the EVD, calculates the rank reduction testing matrix, and performs a 2D RARE spectral search, which require a computational [30].

Maximum Number of Resolvable Sources
In this subsection, we compare the maximum number of resolvable sources for the three algorithms. Note that AUX only utilizes the (M x − 2P + 2)(M y − 2P + 2) sensors in the middle of the URA. As a result of the subspace identification technique, at least one eigenvector is required for spanning the noise subspace. Therefore, AUX is capable of resolving no more than (M x − 2P + 2)(M y − 2P + 2) − 1 sources. Similarly, 2D-RARE requires that (2P − 1) 2 ≤ M x M y -K. Thus, it can also resolve M x M y − (2P − 1) 2 sources at most. For the proposed method, it requires 2M x M y − K ≥ 2M x (2P − 1) and 2M x M y − K ≥ 2(2P − 1) in the two consequent RARE estimators, which means the maximum number of resolvable sources for the proposed method is 2M x M y − 2M x (2P − 1). Therefore, it can be concluded that the alleviation of computational loads is at the cost of DOF reduction, since the signal subspace has been decoupled.

Dependence of Mutual Coupling on Directions
In this subsection, we discuss the dependence of antenna mutual coupling effects on the angles of signals. According to the latest work by B. Friedlander [28], the mutual coupling matrix is a direction dependent matrix, and not a constant matrix as is commonly assumed. In this work, the mutual coupling matrix of the URA is assumed to be direction-independent with a block banded Toeplitz structure, which has been simplified compared with the analytic manifold [28]. Therefore, the DOA estimation performance may deteriorate due to the mismatch between the assumed model and the physical characteristics of the array antennas.

Simulation Results
In this section, we present several numerical simulations to analyze the DOA estimation performance of the proposed algorithm with respect to 2D-MUSIC, AUX, and 2D-RARE. In the following simulations, a 10 × 10-elements URA is employed with the inter-element spacing being d x = d y =λ/2. The impinging sources are BPSK modulated with equal power and uncorrelated to each other, while the noise is amenable to the spatially white complex Gaussian distribution. The mutual coupling length is assumed to be P = 2, and MCCs are c 0,1 = c 1,0 = 0.7527 + j0.4854, c 1,1 = 0.2825 + j0.2801. The estimation performance is assessed in terms of spatial spectrum, root mean-square error (RMSE), and computational complexity. The RMSE is measured by 500 independent Monte Carlo runs: where θ k is the DOA of the k-th source, andθ n,k stands for the estimate of in the nth trial.
In the first case, we investigate the DOA spectra of the three algorithms. We consider that three incident rectilinear signals radiate from (θ 1 = −20 • , φ 1 = 40 • ), (θ 2 = 40 • , φ 2 = 30 • ), and (θ 3 = 80 • , φ 3 = 50 • ), respectively. The corresponding x-axis direction-cosines are (0.604, 0.383, 0.133) in radians, and y-axis direction-cosines are (-0.220, 0.321, 0.754) in radians, respectively. The signal-to-noise ratio (SNR) is 10 dB, and the number of snapshots is fixed at 1000. In Figure 2d, the proposed algorithm generates three sharp peaks corresponding to the y-axis direction-cosines of the three BPSK sources in the first DOA estimation stage. In Figure 2e, the three spectra produce three sharp peaks at the desired locations of the x-axis direction-cosines, respectively. Note that the DOA estimates of in the two stages are automatically paired without any postprocessing. From Figure 2a, it is can be seen that 2D-MUSIC fails to distinguish all three sources, and its two resolvable spectra peaks are biased due to the mismatch caused by the mutual coupling effects. Although AUX and 2D-RARE can generate three peaks corresponding to the sources in the 2D DOA spectrum plain, the former sacrifices eight sensors (i.e., the auxiliary sensors) to make the remaining part of the URA coupling-free while the latter requires exhaustive 2D spectral searches to obtain the joint azimuth and elevation DOA estimates. . The estimation performance is assessed in terms of spatial spectrum, root mean-square error (RMSE), and computational complexity. The RMSE is measured by 500 independent Monte Carlo runs: where k θ is the DOA of the k-th source, and ,n k θ stands for the estimate of k θ in the nth trial.
In the first case, we investigate the DOA spectra of the three algorithms. We consider that three incident rectilinear signals radiate from ( radians, and y-axis direction-cosines are (-0.220, 0.321, 0.754) in radians, respectively. The signal-tonoise ratio (SNR) is 10 dB, and the number of snapshots is fixed at 1000. In Figure 2d, the proposed algorithm generates three sharp peaks corresponding to the y-axis direction-cosines of the three BPSK sources in the first DOA estimation stage. In Figure 2e, the three spectra produce three sharp peaks at the desired locations of the x-axis direction-cosines, respectively. Note that the DOA estimates of in the two stages are automatically paired without any postprocessing. From Figure 2a, it is can be seen that 2D-MUSIC fails to distinguish all three sources, and its two resolvable spectra peaks are biased due to the mismatch caused by the mutual coupling effects. Although AUX and 2D-RARE can generate three peaks corresponding to the sources in the 2D DOA spectrum plain, the former sacrifices eight sensors (i.e., the auxiliary sensors) to make the remaining part of the URA coupling-free while the latter requires exhaustive 2D spectral searches to obtain the joint azimuth and elevation DOA estimates.
(a) 2D-MUSIC without calibration (b) AUX method  In the second experiment, we studied the RMSE of the DOA estimates as a function of the SNR. The number of snapshots remained at 1000. Other simulation settings were the same as the first experiment. Figure 3 depicts the RMSE of the proposed method while the RMSEs of the other three solutions, 2D-MUSIC, AUX, 2D-RARE, and the corresponding CRB for noncircular source DOA estimation in the presence of mutual coupling is also calculated for comparison. It is clear that 2D-RARE performs the best among the four algorithms, overall AUX takes second place, our method provides better accuracies than it at 10 dB and above but is inferior to AUX at low SNRs, and 2D-MUSIC has the worst performance due to the low robustness to the unknown mutual coupling. Besides, the performance of the proposed method, 2D-RARE, and AUX bears down on the CRB asymptotically at high SNRs, while the RMSE of 2D-MUSIC is stabilized at approximately 4  when the SNR is larger than 5dB. In the second experiment, we studied the RMSE of the DOA estimates as a function of the SNR. The number of snapshots remained at 1000. Other simulation settings were the same as the first experiment. Figure 3 depicts the RMSE of the proposed method while the RMSEs of the other three solutions, 2D-MUSIC, AUX, 2D-RARE, and the corresponding CRB for noncircular source DOA estimation in the presence of mutual coupling is also calculated for comparison. It is clear that 2D-RARE performs the best among the four algorithms, overall AUX takes second place, our method provides better accuracies than it at 10 dB and above but is inferior to AUX at low SNRs, and 2D-MUSIC has the worst performance due to the low robustness to the unknown mutual coupling. Besides, the performance of the proposed method, 2D-RARE, and AUX bears down on the CRB asymptotically at high SNRs, while the RMSE of 2D-MUSIC is stabilized at approximately 4 • when the SNR is larger than 5dB.
In the third experiment, we examine the variation of RMSE with the increase of the number of snapshots from 10 to 1000, where the SNR is set as 10dB. Figure 4 demonstrates that RMSEs of the proposed method, 2D-RARE, and AUX decline steadily as the snapshot size raises. This is because that the covariance and elliptic covariance matrices can be estimated more accurately from a larger number of observations. Additionally, the performance of our developed approach is almost the same as 2D-RARE and strictly superior to AUX, asymptotically approaching the corresponding CRB, while the RMSE of 2D-MUSIC saturates at 4 • as similar to the second scenario, which is the worst performance. In the third experiment, we examine the variation of RMSE with the increase of the number of snapshots from 10 to 1000, where the SNR is set as 10dB. Figure 4 demonstrates that RMSEs of the proposed method, 2D-RARE, and AUX decline steadily as the snapshot size raises. This is because that the covariance and elliptic covariance matrices can be estimated more accurately from a larger number of observations. Additionally, the performance of our developed approach is almost the same as 2D-RARE and strictly superior to AUX, asymptotically approaching the corresponding CRB, while the RMSE of 2D-MUSIC saturates at 4 as similar to the second scenario, which is the worst performance. In the last experiment, the computational burdens of the three algorithms versus the variation of the number of sensors were compared. The searching step size for DOA estimation was set as 0.1°,  In the third experiment, we examine the variation of RMSE with the increase of the number of snapshots from 10 to 1000, where the SNR is set as 10dB. Figure 4 demonstrates that RMSEs of the proposed method, 2D-RARE, and AUX decline steadily as the snapshot size raises. This is because that the covariance and elliptic covariance matrices can be estimated more accurately from a larger number of observations. Additionally, the performance of our developed approach is almost the same as 2D-RARE and strictly superior to AUX, asymptotically approaching the corresponding CRB, while the RMSE of 2D-MUSIC saturates at 4 o as similar to the second scenario, which is the worst performance.  In the last experiment, the computational burdens of the three algorithms versus the variation of the number of sensors were compared. The searching step size for DOA estimation was set as 0.1 • , and the number of snapshots was fixed at 200. Suppose the URA is square, i.e., M x = M y , and the number of elements ranges from 10 to 40. It can be observed from Figure 5 that our method is the most efficient compared with other two as it is only reliant upon a 1D spectral search, followed by AUX, and 2D-RARE requires the highest complexity.
In the last experiment, the computational burdens of the three algorithms versus the variation of the number of sensors were compared. The searching step size for DOA estimation was set as 0.1°, and the number of snapshots was fixed at 200. Suppose the URA is square, i.e., , and the number of elements ranges from 10 to 40. It can be observed from Figure 5 that our method is the most efficient compared with other two as it is only reliant upon a 1D spectral search, followed by AUX, and 2D-RARE requires the highest complexity.

Conclusions
In this paper, a 2D real-valued DOA estimator that is robust to unknown mutual coupling was developed for rectilinear sources. Taking advantage of the noncircularity in signals and the structured mutual coupling matrix revealed as block Toeplitz, a RARE estimator in real-value field was proposed to extract the angular information from the coupling coefficients. Compared with stateof-the-art approaches, our solution has satisfactory estimation performance with more efficient computations. A series of numerical simulations showed that the proposed method has its own advantages on estimation accuracy and computational complexity over previous algorithms with a favorable robustness to unknown mutual coupling. To apply the proposed algorithm to further applications, we may extend the proposed algorithm to considering the physical characteristics of the array antennas, such as direction-dependent mutual coupling and unstructured mutual coupling matrix in further studies.
where R¯x is the extended covariance matrix corresponding to the augmented signalx(t) = T , which is formed as: The FIM for the unknown parameters and their cross terms can be readily derived following the steps in [43]. Since F is a symmetric matrix, we only need to find the upper triangular parts of F, and the block matrices in F can be obtained as:   where the partial derivative of A with respect to the parameters θ, φ, ϕ, c R , c I in the above equations can be calculated as: Since the noise power is not coupled with other parameters, the cross terms with respect to σ 2 n are all zeros. Consequently, the FIM can be formulated as the following form: F θθ F θφ F θϕ F θc r F θc I F θρ 0 F φθ F φφ F φϕ F φc r F φc I F φρ 0 F ϕθ F ϕφ F ϕϕ F ϕc r F ϕc I F ϕρ 0 F c r θ F c r φ F c r ϕ F c r c r F c r c I F c r ρ 0 F c I θ F c I φ F c I ϕ F c I c r F c I c I F c I ρ 0 F ρθ F ρφ F ρϕ F ρc r F ρc I F ρρ 0 0 Sensors 2020, 20, 1914 21 of 22 Subsequently, the CRB for DOA parameters is defined as: