2D Direction of Arrival Estimation for Cross Array in the Presence of Mutual Coupling

Abstract

This paper proposes a new method for cross array to estimate two-dimensional direction of arrival (2-D DOA) in the presence of mutual coupling. In this method, the array elements which are affected by the same mutual coupling are chosen on 𝑥-axis and 𝑧-axis, respectively. Then a new matrix is constructed with the proper entries of cross covariance matrix of the chosen elements outputs on 𝑥-axis and 𝑧-axis. Propagation method (PM) and rotational invariance techniques for uniform linear array (ULA) are utilized in the constructed matrix to obtain two parameters correlated with elevations and azimuths. While calculating and pairing the two parameters, only once eigendecomposing and several division operations are required with the relationship among the matrix, its eigenvalues, and corresponding eigenvectors. Simulations are presented to validate the performance of the proposed method.

1. Introduction

Mutual coupling between elements has a serious effect on precision and resolution of the direction of arrival estimation algorithms [1, 2]. Therefore, there are lots of works addressing this issue of how to restrain the effects of mutual coupling. Many calibration algorithms have been proposed in the past decades. Friedlander and Weiss presented an iterative procedure to compensate the mutual coupling and perturbation of gain and phase [3]. In [4], Sellone and Serra constructed an objective function based on covariance matrix matching, dealt with mutual coupling matrix and its conjugate transpose as unrelated ones, and searched the optimal solution with iterative procedure. However, the multidimensional search needed for the associated nonlinear optimal search is computationally complicated in [3, 4]. Wang et al. made use of banded complex symmetric Toeplitz characteristic of the mutual coupling matrix of ULA and rank reduction (RARE) algorithm to obtain azimuths and mutual coupling coefficients [5], while Lin and Yang utilized banded complex symmetric circular Toeplitz characteristic of the mutual coupling matrix of uniform circle array (UCA) [6]. These two methods are only suitable for ULA and UCA, respectively. The idea from [5] was applied to cross array, but the transformed matrix related to the ideal steering vector has no universal form for reasons of the array configuration and the variable mutual coupling degree of freedom [7]. The literature [8] presented a decoupling algorithm for ULA. It chooses the proper elements whose ideal steering vector can be separated from the mutual coupling coefficients and then employs the outputs of the chosen elements to construct MUSIC spatial spectral estimator. Nevertheless this method can only estimate one-dimensional direction and needs the spectral peak search. In [9, 10], this method [8] was applied to uniform rectangular array (URA) and uniform hexagon array (UHA), respectively, to estimate two-dimensional directions, but the configurations of planar arrays are complicated and lots of auxiliary elements are needed.

In order to estimate 2-D DOA in the presence of mutual coupling, a method for cross array which has a simpler configuration than those for planar arrays is presented in this paper. The method is based on PM algorithm utilizing signal subspace eigenvectors algorithm for the array consisting of the chosen elements which are affected by the same mutual coupling. It is not an iterative algorithm, and the spectral peak search is not required. It can estimate the 2-D DOA accurately without estimating the mutual coupling coefficients.

2. Problem Formulation

Consider a cross array consisting of two ULAs located on 𝑥-axis and 𝑧-axis, respectively. The ULAs share the middle element as reference at origin of coordinates. Each ULA has 2𝑀+1sensors with element spacing of 𝑑, which means that there exist 4𝑀+1ones in the whole array because of the sharing element. Assume that𝐾narrowband, far-field, noncoherent signals with the same wave length 𝜆 (𝜆=2𝑑) impinge on this array, and the 𝑘th source’s elevation and azimuth are 𝛼_𝑘 and 𝛽_𝑘, respectively. The array configuration is shown in Figure 1.

Figure 1 

The cross array configuration.

Let 𝐂 and 𝑡 denote the mutual coupling matrix (MCM) and the time index, respectively. The output of the array can be expressed as 𝐘(𝑡)=𝐂𝐀𝐒(𝑡)+𝐍(𝑡), (1) where 𝐘, 𝐀, 𝐒, and 𝐍 represent the received signal vector, the ideal steering matrix, the source signal vector, and the stochastic noise (additive white Gauss noise whose mean is zero and variance is 𝜎²) vector, respectively. They can be described as (the time index𝑡is omitted for brevity) 𝐘=𝑥₁,…,𝑥_𝑀,…,𝑥_2𝑀+1,𝑧₁,…,𝑧_𝑀,…,𝑧_2𝑀+1^𝑇,𝐀=𝐚𝛼₁,𝛽₁,…,𝐚𝛼_𝑘,𝛽_𝑘,…,𝐚𝛼_𝐾,𝛽_𝐾,𝐚𝛼_𝑘,𝛽_𝑘=𝑒^{𝑗2𝜋𝑑𝑀cos𝛼𝑘cos𝛽𝑘/𝜆}𝑒,…,1,…,_{−𝑗2𝜋𝑑𝑀cos𝛼𝑘cos𝛽𝑘/𝜆},𝑒^{𝑗2𝜋𝑑𝑀sin𝛼𝑘/𝜆},…,1,…,𝑒^{−𝑗2𝜋𝑑𝑀sin𝛼𝑘/𝜆𝑇},𝐒=𝑠₁,…,𝑠_𝑘,…,𝑠_𝐾^𝑇,𝐍=𝑛_𝑥1,…,𝑛_𝑥𝑀,…,𝑛_𝑥2𝑀+1,𝑛_𝑧1,…,𝑛_𝑧𝑀,…,𝑛_𝑧2𝑀+1^𝑇, (2) where superscript 𝑇 denotes transpose.

As shown in [3], the MCM of ULA can be expressed as a banded complex symmetric Toeplitz matrix𝐁𝐁=Toeplitz1,𝑏₁,…,𝑏_𝑝, (3) where integer 𝑝 denotes the mutual coupling degree of freedom, which means the coefficient of mutual coupling is assumed to be zero if the distance between the two sensors exceeds𝑝𝑑.

Based on the MCM model of ULA by (3), the MCM of cross array has the following form: 𝐂=𝐁𝐃𝐃𝐁, (4) where 𝐁 is the MCM of ULA on 𝑥- or 𝑧-axis, and 𝐃 describes the mutual coupling between the two ULAs. The matrix 𝐃 can be presented by 𝐃=𝐉𝐃₁𝐉𝐡₁^𝑇𝐉𝐃₁𝐡₁1𝐡₁𝐉𝐃₁𝐉𝐡₁𝐉^𝑇𝐃₁, (5) where 𝐡₁=𝟎_{(𝑀−𝑝)×1}𝑏_𝑝⋯𝑏₁,𝐉=00⋯100⋰⋮⋮1⋱01⋯00^𝑀×𝑀. (6)𝐃₁=(𝑑_𝑖𝑗)_𝑀×𝑀,𝑑_𝑖𝑗denotes the coefficient of mutual coupling between the two elements spacing 𝑖²+𝑗²𝑑, and 𝑑_𝑖𝑗 satisfies 𝑑_𝑖𝑗=𝑑_𝑗𝑖,𝑑_𝑖𝑗=≠0𝑖²+𝑗²<𝑝=0𝑖²+𝑗²>𝑝. (7)

3. Proposed Method

The ideal steering vector of ULA elements affected by the same mutual coupling can be separated from the mutual coupling coefficients, which means they can be decoupled easily. Let 𝐘^′ and 𝐘^′′ be the outputs of elements affected by the same mutual coupling on𝑥-axis and𝑧-axis, respectively.

𝐘^′ can be written as 𝐘^′==𝐇𝐘_𝐾𝑘=1𝑏_𝑝𝑢_𝑘^𝑀+⋯+𝑢_𝑘^𝑀−𝑝+⋯𝑏_𝑝𝑢_𝑘^𝑀−2𝑝⋮𝑏_𝑝𝑢_𝑘^2𝑝+⋯+𝑢_𝑘^𝑝+⋯𝑏_𝑝𝑏_𝑝+⋯+𝑢_𝑘^−𝑝+⋯𝑏_𝑝𝑢_𝑘^−2𝑝⋮𝑏_𝑝𝑢_𝑘^2𝑝−𝑀+⋯+𝑢_𝑘^𝑝−𝑀+⋯𝑏_𝑝𝑢_𝑘^−𝑀𝑠_𝑘+𝐍₁=^𝐾_𝑘=1𝑓_𝑘𝑠_𝑘𝐚_𝑘+𝐍₁=𝐀₁𝐓₁𝐒+𝐍₁, (8) where 𝑓_𝑘=𝑏_𝑝𝑢_𝑘^𝑀+⋯+𝑢_𝑘^𝑀−𝑝+⋯𝑏_𝑝𝑢_𝑘^𝑀−2𝑝, 𝐇=𝐇₁𝟎_{}2(𝑀−2𝑝+1)×(2𝑀+1)},𝐇₁=^}_{(𝑀−2𝑝}𝟎_+1)×𝑝𝐈𝟎^}_{(𝑀−2𝑝+1)×2𝑝−1}𝟎𝟎^}_{(𝑀−2𝑝+1)×𝑝},𝐚𝟎𝟎𝟎𝐈𝟎_𝑘=1,…,𝑢_𝑘^{−(𝑀−2𝑝)},𝑢_𝑘^−𝑀,…,𝑢_𝑘^{−2𝑀+2𝑝𝑇},𝑢_𝑘=exp𝑗2𝜋𝑑cos𝛼_𝑘cos𝛽_𝑘𝜆,𝐀₁=𝐚₁,…,𝐚_𝑘,…,𝐚_𝐾,𝐓₁=diag𝑓₁,…,𝑓_𝑘,…,𝑓_𝐾. (9)

There are no entries of 𝐃 in (8) as the chosen elements on 𝑥-axis are not affected by the sensors on𝑧-axis, because the distances between them are out of 𝑝𝑑.

Similarly, 𝐘^′′ can be expressed as 𝐘^′′=𝐀₂𝐓₂𝐒+𝐍₂, (10) where 𝐀₂=[𝐛₁,…,𝐛_𝑘,…,𝐛_𝐾], 𝐓₂=diag𝑔₁,…,𝑔_𝑘,…,𝑔_𝐾,𝐛_𝑘=1,…,𝑣_𝑘^{−(𝑀−2𝑝)},𝑣_𝑘^−𝑀,…,𝑣_𝑘^{−2𝑀+2𝑝𝑇},𝑔_𝑘=𝑏_𝑝𝑣_𝑘^𝑀+⋯+𝑣_𝑘^𝑀−𝑝+⋯𝑏_𝑝𝑣_𝑘^𝑀−2𝑝,𝑣_𝑘=exp𝑗2𝜋𝑑sin𝛼_𝑘𝜆. (11)(12)(13)(14)

In order to introduce rotational invariance techniques for ULA to this array, 𝐘^′ is partitioned into two parts, so is 𝐘^′′. They are denoted by 𝐘₁₁, 𝐘₁₂, 𝐘₂₁, and 𝐘₂₂, respectively: 𝐘₁₁=𝐘^′_{(1∶𝑀−2𝑝)}𝐘^′_{(𝑀−2𝑝+2∶2𝑀−4𝑝+1)},𝐘₁₂=𝐘^′_{(2∶𝑀−2𝑝+1)}𝐘^′_{(𝑀−2𝑝+3∶2𝑀−4𝑝+2)},𝐘₂₁=𝐘^′′_{(1∶𝑀−2𝑝)}𝐘^′′_{(𝑀−2𝑝+2∶2𝑀−4𝑝+1)},𝐘₂₂=𝐘^′′_{(2∶𝑀−2𝑝+1)}𝐘^′′_{(𝑀−2𝑝+3∶2𝑀−4𝑝+2)}. (15)

𝐘₁₁ and 𝐘₂₁ are the submatrixes of 𝐘^′ and 𝐘^′′, respectively, so the cross covariance matrix of 𝐘₁₁ and 𝐘₂₁ can be built with the proper entries of cross covariance matrix of 𝐘^′ and 𝐘^′′: 𝐑₁=𝐀₁₁𝐑_𝑠𝑠𝐀^𝐻₂₁=𝐑₁₁𝐑₁₂𝐑₁₃𝐑₁₄, (16) where 𝐀₁₁ and 𝐀₂₁ are the ideal steering matrixes of arrays corresponding to 𝐘₁₁ and 𝐘₂₁, respectively, 𝐑=𝐸𝐘^′𝐘^′′𝐻𝐑₁₁,=𝐑_{(1∶𝑀−2𝑝,1∶𝑀−2𝑝),}𝐑₁₂=𝐑_{(1∶𝑀−2𝑝,𝑀−2𝑝+2∶2𝑀−4𝑝+1)},𝐑₁₃=𝐑_{(𝑀−2𝑝+2∶2𝑀−4𝑝+1,1∶𝑀−2𝑝),}𝐑₁₄=𝐑_{(𝑀−2𝑝+2∶2𝑀−4𝑝+1,𝑀−2𝑝+2∶2𝑀−4𝑝+1)},𝐑_𝑠𝑠=𝐓₁𝐸𝐒𝐒^𝐻𝐓₂^𝐻. (17)(18)(19)(20)(21)(22)

Similarly, let 𝐑₂, 𝐑₃, and 𝐑₄, which can be constructed with the proper entries chosen from 𝐑 as 𝐑₁ in (16), be the cross covariance matrixes of 𝐘₁₂ and 𝐘₂₁, 𝐘₁₁, and 𝐘₂₂, and 𝐘₁₂ and 𝐘₂₂, respectively, 𝐑₂=𝐀₁₂𝐑_𝑠𝑠𝐀^𝐻₂₁=𝐀₁₁𝚽₁𝐑_𝑠𝑠𝐀^𝐻₂₁𝐑₃,=𝐀₁₁𝐑_𝑠𝑠𝐀^𝐻₂₂=𝐀₁₁𝐑_𝑠𝑠𝐀₂₁𝚽₂^𝐻=𝐀₁₁𝚽₂^𝐻𝐑_𝑠𝑠𝐀^𝐻₂₁,𝐑₄=𝐀₁₂𝐑_𝑠𝑠𝐀^𝐻₂₂=𝐀₁₁𝚽₁𝚽₂^𝐻𝐑_𝑠𝑠𝐀^𝐻₂₁, (23)(24)(25) where 𝐀₁₂ and 𝐀₂₂ are the ideal steering matrixes of arrays corresponding to 𝐘₁₂ and 𝐘₂₂, respectively, 𝚽₁=diag𝑢₁⁻¹,…,𝑢_𝑘⁻¹,…,𝑢_𝐾⁻¹,𝚽₂=diag𝑣₁⁻¹,…,𝑣_𝑘⁻¹,…,𝑣_𝐾⁻¹. (26)

The new matrix 𝐑^′ consisting of the four cross covariance matrixes includes the rotational invariance factors, 𝐑^′=𝐑₁𝐑₂𝐑₃𝐑₄=𝐀₁₁𝐀₁₁𝚽₁𝐀₁₁𝚽₂^𝐻𝐀₁₁𝚽₁𝚽₂^𝐻𝐑_𝑠𝑠𝐀^𝐻₂₁. (27)

Similar to the conventional propagator method, 𝐑^′ can be partitioned as 𝐑^′=𝐆^𝑇𝐐^𝑇𝑇, (28) where 𝐆 and 𝐐 are the first 𝐾 and the last 8(𝑀−2𝑝)−𝐾 rows, respectively.

A matrix 𝐏 (the propagator parameter) is defined to satisfy 𝐏=min_𝐏‖𝐐−𝐏𝐆‖². (29)

The solution of (29) can be obtained as follows: 𝐏=𝐐𝐆⁺=𝐐𝐆^𝐻𝐆𝐆^𝐻−1. (30)

With the relationship between the propagator parameter and ideal steering matrix, we obtain an equation as 𝐏^′=𝐈𝐏=𝐏₁𝐏₂𝐏₃𝐏₄=𝐀₁₁𝐀₁₁𝚽₁𝐀₁₁𝚽₂^𝐻𝐀₁₁𝚽₁𝚽₂^𝐻𝐓₃, (31) where 𝐏_𝑖 (𝑖=1,2,3,4) is the (2𝑀−4𝑝)×𝐾 matrix and 𝐓₃is the 𝐾×𝐾 matrix with rank 𝐾.

Similar to ESPRIT, from (31) we can have the equations correlated with azimuths and elevations as 𝐏₁⁺𝐏₂=𝐓₃⁻¹𝚽₁𝐓₃,𝐏₃⁺𝐏₄=𝐓₃⁻¹𝚽₁𝐓₃,𝐏₁⁺𝐏₃=𝐓₃⁻¹𝚽₂^𝐻𝐓₃,𝐏₂⁺𝐏₄=𝐓₃⁻¹𝚽₂^𝐻𝐓₃. (32)

To solve the above equation sets and pair the azimuths and elevations, four times eigen-decomposing and 3[𝐾(𝐾+1)/2−1] times calculating inner product of vectors are necessary in a general way. Now a method, which needs only once eigen-decomposing and some division operations to obtain and pair the azimuths and the elevations, is introduced with the relationship among the matrix, its eigenvalues, and corresponding eigenvectors.

First eigendecompose 𝐏₁⁺𝐏₂, and the eigenvalues and the eigenvectors are denoted by 𝜆_1𝑘 and 𝝆_𝑘 (𝑘=1,2,…,𝐾), respectively Then substitute 𝝆_𝑘 into other equations, as the first one of the right equation sets in (32); the corresponding eigenvalue 𝜆_3𝑘 can be obtained from 𝐏₁⁺𝐏₃𝝆_𝑘=𝜆_3𝑘𝝆_𝑘. (33)

Let 𝐏₁⁺𝐏₃𝝆_𝑘=𝝎_𝑘, and 𝜆_3𝑘 can be expressed by 𝜆_3𝑘=1𝐾^𝐾_𝑖=1𝝎_𝑘𝑖𝝆_𝑘𝑖, (34) where 𝝎_𝑘𝑖 and 𝝆_𝑘𝑖 are the 𝑖th element of 𝝎_𝑘 and 𝝆_𝑘, respectively.

Assume that 𝜆_2𝑘 and 𝜆_4𝑘 denote the corresponding eigenvalues of 𝐏₃⁺𝐏₄, and 𝐏₂⁺𝐏₄ respectively, and they can be calculated by (34).

Using the relationship among 𝜆_𝑖𝑘 (𝑖=1,2,3,4), 𝑢_𝑘, 𝑣_𝑘, 𝛼_𝑘, and 𝛽_𝑘, the elevation 𝛼_𝑘 and the azimuth 𝛽_𝑘 of the 𝑘th source can be represented as 𝛼_𝑘=arcsinarg𝜆_3𝑘/𝜋+arcsinarg𝜆_4𝑘/𝜋2,𝛽_𝑘=arccos−arg𝜆_1𝑘/𝜋cos𝛼_𝑘+arccos−arg𝜆_2𝑘/𝜋cos𝛼_𝑘2. (35)

Now give the computation load on the proposed method and PM algorithm with known mutual coupling. The computation load on PM algorithm with known mutual coupling is focused on the inverse of MCM, cross matrix of output of elements on 𝑥-axis and 𝑧-axis, eigen-decomposing a matrix with rank 𝐾, and pairing the parameters, so the computation load is 𝑂((4𝑀+2)³)+𝑁(2𝑀+1)²+𝑂(4𝐾³)+3𝐾²[𝐾(𝐾+1)/2−1](𝑁is the snapshots number); the computation load on the proposed method is focused on cross matrix of output of elements on 𝑥-axis and 𝑧-axis, once eigen-decomposing a matrix with rank 𝐾, and pairing the parameters, so the computation load is 𝑁(2𝑀−4𝑝+2)²+𝑂(𝐾³)+3𝐾(𝐾²+𝐾).

Finally, the necessary conditions for the proposed method are given as follows:(1)for unique estimation, 2(𝑀−2𝑝)>𝐾 is required;(2)the sources can not be located at the blind angles which make the determinant of 𝐓₁ or 𝐓₂ zero.

4. Simulation Experiments

In this section, we carry out some representative simulation experiments to demonstrate the validity of the proposed calibration method. In all the experiments, we assume there are two uncorrelated equivalent-power sources at (20°, 10°) and (80°, 60°) that impinge on the cross array which consists of 53 elements (𝑀=13). The mutual coupling degree of freedom 𝑝 is two and the coefficients of mutual coupling are denoted by 𝑏₁=0.3065+0.3951𝑗, 𝑏₂=0.1018−0.1721𝑗, and 𝑑₁₁=0.2154−0.2738𝑗.

In the first experiment, the 2-D spatial spectrums of the PM algorithm with unknown mutual coupling, the PM algorithm with known mutual coupling, and the proposed method (estimate the 2-D DOAs to get the coefficients of mutual coupling with the method in [7]) are shown in Figures 2 to 4 in the condition that the signal-to-noise ration (SNR) is 5 dB and the number of the snapshots is 300.

Figure 2 

2-D spatial spectrum of the PM algorithm in the presence of unknown mutual coupling.

From Figure 2 we can see that the PM algorithm without calibration has very poor resolution in the presence of mutual coupling. Though the peaks appear at the signal directions, there exist lots of pseudopeaks especially at (80°, 10°) whose spatial spectrum is greater than the sources’, which may result in the inaccurate estimation.

Comparing with Figure 4 and Figure 3, the spatial spectrum of proposed method has two sharp peaks at the signal directions, which is also shown in Figure 3. In other words, the proposed method can restrain the effects of mutual coupling and has a high resolution. Though the performance of proposed method is worse than that in Figure 3 for reasons of the array elements consumption, it eliminates pseudo peaks caused by mutual coupling, which can distinguish the two sources visibly.

Figure 3 

2-D spatial spectrum of the PM algorithm in the presence of known mutual coupling.

Figure 4 

2-D spatial spectrum of the proposed method.

In the second experiment, the effect of SNR on the performance of proposed method is explored. The number of snapshots is 300 and the SNR varies from −10 to 10 dB with the interval 1 dB. The root mean square error (RMSE) which is defined as (36) is shown in Figure 5 with 100 independent experiments: 1RMSE=2²_𝑘=1𝐸𝛼_𝑘−𝛼_𝑘²+𝛽_𝑘−𝛽_𝑘². (36)

Figure 5 

RMSE versus SNR.

The simulation results show the proposed method has an obvious performance improvement over the PM algorithm with unknown coupling, though it is a little inferior to the PM algorithm with known coupling as we think. As the SNR increases, the difference of the RMSE between the proposed method and the PM algorithm with known coupling gets closer; however, the accuracy of DOA estimates of the PM algorithm with unknown coupling stays at a low level as the main error is caused by the manifold error between the ideal and the actual one.

In the last experiment, the influence of the number of the snapshots on the performance of proposed method is investigated. SNR is 0 dB and the number of the snapshots varies from 100 to 800 with the interval 100. 100 independent experiments are carried out, and the RMSE curves are shown in Figure 6.

Figure 6 

RMSE versus number of snapshots.

From Figure 6 we can see that the result is similar to that of the second experiment: the capability of proposed method is close to the one of the PM algorithm with known coupling, and the performance of the PM algorithm with unknown coupling is still poor when the number of the snapshots increases.

From the above three experiments, it is seen that the proposed method has compensated the effect caused by mutual coupling. It provides a significant improvement of DOA estimation accuracy and can achieve almost the same performance as that of PM algorithm with known mutual coupling.

5. Conclusion

A method compensating the effect of the mutual coupling is proposed for cross array to estimate the 2-D DOA. This method is based on PM algorithm and rotational invariance techniques for ULA. It can compensate the effect of the mutual coupling without any calibration sources, and it does not require any iterative procedure or spectral peak search. Simulations have shown that the proposed method can improve the DOA estimation accuracy well.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to thank the editor and reviewers for their helpful comments and suggestions. In addition, this work is supported by the National High Technology Research and Development Program of China (863 Program) under Grant 2012AA0601101, the Key Laboratory of Intelligent Perception and Systems for High-Dimensional Information (Nanjing University of Science and Technology), Ministry of Education under Grant 3092013012205, and Research Fund for the Doctoral Program of Higher Education of China under Grant 20133219110027.