DOA Estimation under Unknown Mutual Coupling and Multipath with Improved Effective Array Aperture

Subspace-based high-resolution direction of arrival (DOA) estimation significantly deteriorates under array manifold perturbation and rank deficiency of the covariance matrix due to mutual coupling and multipath propagation, respectively. In this correspondence, the unknown mutual coupling can be circumvented by the proposed method without any passive or active calibration process, and the DOA of the coherent signals can be accurately estimated accordingly. With a newly constructed matrix, the deficient rank can be restored, and the effective array aperture can be extended compared with conventional spatial smoothing. The proposed method achieves a good robustness and DOA estimation accuracy with unknown mutual coupling. The simulation results demonstrate the validity and efficiency of the proposed method.


Introduction
In array signal processing, direction of arrival (DOA) estimation of incident signals is an important research issue. Many high-resolution methods based on the orthogonality of subspaces have been studied extensively over the years, such as MUSIC [1] and ESPRIT [2]. Most of these algorithms assume that the array manifold is known and the signals are uncorrelated. However, these two prerequisites are not always guaranteed in practice as mutual coupling between the elements is an intrinsic characteristic of antenna arrays and multipath propagation caused by reflective surfaces are very common in urban areas, both resulting in a failure of subspace-based methods [3,4].
In electromagnetic engineering, the mutual coupling matrix can be calculated using techniques such as finite element or method of moment techniques [5][6][7]. However, these calculation depends heavily on the physical shape of the antennas and are specific for each antenna array. It is thus often preferable to carry out a joint estimation of the mutual coupling matrix and signal DOA from a signal processing perspective. Weiss and Friedlander [3] proposed an iterative algorithm that compensates for mutual coupling and gain and phase perturbations, in uniform linear and circular arrays. Subsequently, Sellone et al. proposed another iterative method, which alternatingly minimises a cost function with respect to two complex valued symmetric Toeplitz matrices and a complex Hermitian Toeplitz matrix [8]. Compared with the method proposed in [3], this method is less sensitive to the array perturbations and can get an effective DOA estimate without a preliminary estimate in many cases. However, the computational complexities of both are high because the non-linear multidimensional search process is random in nature.
More recently, several DOA estimation algorithms which are not sensitive to mutual coupling have been developed [9][10][11]. By setting a group of sensors as auxiliary sensors, the DOAs can be directly estimated without compensating for mutual coupling. Moreover, the result can be refined iteratively by estimating the mutual coupling coefficients for compensation. The main advantage of these methods is that no iteration is required but the aperture is reduced by making some sensors auxiliary.
All the algorithms reviewed above are only effective in a scenario where the incident signals are uncorrelated to each other and do not undergo multipath propagation. In a typical wireless communication system where multipath propagation is unavoidable, these algorithms may not be effective.
To tackle these two thorny issues, a class of auto-calibration methods, which are able to estimate the mutual coupling coefficients and DOAs of coherent signals jointly, has been addressed in [12,13]. In [12], the mutual coupling coefficients and coherent DOAs can be estimated jointly in an alternating manner, but this results in a large computational complexity for the nonlinear multidimensional search and is not applicable in some systems due to the requirement of pilot symbols. A noniterative optimisation method is given in [13] to estimate the DOAs of mixed signals which include coherent signals, but the mutual coupling coefficients have to be estimated from uncorrelated signals first. Hence, it may not be effective when there are only coherent signals. Recently, a spatial smoothing scheme studied in [14] estimates the DOAs of coherent signals and is insensitive to unknown mutual coupling, but only the forward smoothing was considered and the conjugate information contained within the signals was not utilised. More recently, Toeplitz matrix reconstruction directly in data domain [15] showed good performance for real-time applications in decorrelating coherent signals without mutual coupling compensation, but at the cost of halving the effective array aperture which is very limiting in practice.
In this correspondence, we propose a new matrix reconstruction method which is insensitive to mutual coupling and has an improved aperture after rank restoration. The main differences between our method and the methods in [14,15] lie in: 1.
Our method exploits the conjugate of the received data for the purpose of rank restoration, which gives rise to more degrees of freedom (DOF) as well as an extended aperture, while both the algorithms in [14,15] ignore such information, which significantly restricts the identifiability of DOA estimates.

2.
Our method makes use of all entries of the covariance matrix of the middle subarray, and the performance in terms of estimation resolution will thus be improved, whereas the spatial smoothing technique adopted in [14], only averages the covariance matrices of the subarrays but does not take advantage of the cross correlations between them, and thus the estimation performance will be compromised to some extent.

3.
Our approach offers a clear advantage in how it adapts to changes in the number of signals, whereas the method in [15] has a fixed number of DOF regardless of the number of signals.
These differences combine together so that our method can offer better estimation performance than existing techniques that incorporate forward-only averaging.
The rest of the paper is organised as follows. In Section 2, the signal model when coherent signals and mutual coupling coexist is presented. In Section 3, we develop a preprocessing strategy to estimate the DOAs of the coherent signals which is blind to the mutual coupling effects and does not require array calibration. In Section 4, simulation results show the validity and efficiency of our proposed method. Finally, some concluding remarks are given in Section 5.
Throughout this paper, the notations that will be used are listed as follows. The operators (·) T , (·) * , (·) H , (·) −1 , E[·], · F , · , and · denote the operation of transpose, conjugate, conjugate transpose, inverse, expectation, the Frobenius norm, ceiling, and flooring of a decimal number, respectively. The symbol diag{z 1 , z 2 } represents a diagonal matrix with diagonal entries z 1 , z 2 and blkdiag{Z 1 , Z 2 } represents a block diagonal matrix with diagonal entries Z 1 , Z 2 . The symbol rank(Z) denotes the rank of a matrix Z. The symbol Z(a : b, c : d) denotes a constructed submatrix by the entries from rows a to b and columns c to d of Z, and the symbol Z(a, b) denotes the entry in the a-th row and b-th column of Z.

Problem Formulation
Consider a number of N narrowband far-field coherent signals impinging on a uniform linear array (ULA) with M identical omnidirectional sensors. Assume that these signals are classified into K groups, which result from K statistically uncorrelated far-field sources s k (t), k = 1, 2, · · · , K with P k multipath signals for each source. In the k-th coherent group, the signal coming from direction θ kp , p = 1, 2, · · · , P k corresponds to the p-th multipath propagation of the source s k (t), and the complex fading coefficient is α kp . It is readily seen that the total number of coherent signals satisfies N = ∑ K k=1 P k . Considering the effect of mutual coupling between the array elements, the corresponding M × 1 array output vector is then given by where a(θ) = 1, e j 2πd λ sin θ , · · · , e j 2π(M−1)d λ sin θ T ∈ C M is the steering vector with λ and d being the wavelength of carrier signal and the spacing between adjacent elements, respectively, C denotes the mutual coupling matrix (MCM), A = A 1 , · · · , A K with A k = a(θ k1 ), · · · , a(θ kP k ) , containing attenuation information of the k-th coherent group, s(t) = s 1 (t), · · · , s K (t) T , and n(t) is white Gaussian noise with the power σ 2 n for each entry. Besides, we assume that the array manifold A is unambiguous, i.e., the steering vectors {a(θ i )} N i=1 are linearly independent for any set of distinct {θ i } N i=1 . As described in [3,[12][13][14][15], it is often sufficient to consider the ULA coupling model that has just finite non-zero coefficients, and a banded symmetric Toeplitz matrix can be used to model mutual coupling. To elaborate on the structure of the MCM, we consider the coupling effects to the i-th element shown in Figure 1 as a instance, where S i,u = c |u−i| A(u, :)Γs(t), 1 ≤ i, u ≤ M, is the coupling contribution from the u-th sensor to the i-th sensor with c |u−i| being the mutual coupling coefficient. It is known that the effect of mutual coupling between two sensor elements is inversely related to their distance and can be ignored when the separation is more than a few wavelengths. More precisely, we assume that when the distance between two sensors is more than Pd where d is the inter-element spacings in a uniform linear array, the mutual coupling coefficient can be approximated as zero. We also assume that the mutual coupling coefficients only depend on the distance between two sensor elements. Based on these assumptions, the observed output at the i-th sensor can be expressed as where n i (t) is the additive noise at the i-th sensor. Considering mutual coupling for all the M sensors and stacking x 1 (t), x 2 (t), · · · , x M (t) in a column, one obtains Equation (1) where with 0 < |c 1 |, |c 2 |, · · · , |c P−1 | < c 0 = 1. From Equation (1), the array covariance matrix is given by 2 1 , · · · , σ 2 K } is the signal covariance matrix, and I M represents an M × M identity matrix. In the case of finite snapshots, the array covariance matrix can be calculated asR

Mutual Coupling Circumvention
Referring to [9][10][11], we know that the middle subarray, defined as the middle M − 2P + 2 elements, is insensitive to mutual coupling, and its actual steering vector is equivalent to an ideal one (i.e., no coupling effect) scaled by a scalar. In order to combat signal coherency as well as effect of unknown mutual coupling in the ULA, we only make use of the output of the middle subarray by virtue of their steering vector. All other elements are defined to belong to auxiliary subarrays as shown in Figure 2. Inspired by [9], a selection matrix J is defined as and the covariance matrix of the the middle subarray is then given bỹ auxiliary subarray middle subarray auxiliary subarray Figure 2. Allocation of the middle subarray and auxiliary subarray.

DOA Estimation of Coherent Signals
In this section, we will propose a rank restoration algorithm using an improved matrix reconstruction technique for DOA estimation of coherent signals in the presence of unknown mutual coupling.
Here a new matrix from the i-th row ofR x is constructed as where (1) and (9), the received signalx i (t) at the i-th element of the middle array is where β kp = e j 2πd λ sin θ kp and n i (t) are the white Gaussian noise at the i-th sensor. Then, the (i, j)-th entry ofR x can be expressed as where δ i,j is the Kronecker delta. Substituting Equation (12) into the matrix in Equation (10) allows B(i) to be written as where and their entries are defined as It is known that the effect of having only a finite number of snapshots, especially when the number is small, may cause biased estimates of subspaces. To further restore the rank and make the matrix reconstruction method more robust against the effects of noise and finite snapshots, all the M − 2P + 2 rows ofR x are then exploited to construct the following square matrix It can be readily verified that Defining Next we examine whether the rank of R d has been restored sufficiently to resolve all N incident coherent signals.
Proof: Since A c is unambiguous and m > N, one has rank(A c ) = N. It is easy to identify that rank Referring to the Lemma 1 in [16], we know that given H ∈ C I×J , the diagonal matrix Q has no identical entries on the diagonal, and rank (H) = r < I, then rank H, QH = r + 1. Substituting back into our problem, we find that the diagonal entries of Φ are not the same as each other, and rank (b(i)) = 1, and thus rank b(i), Φb(i) = 2.
If Lemma 1 is applying successively to the submatrices in G(i), one has a) if 2q ≥ N, then rank (G(i)) = N; b) rank (G(i)) = 2q, otherwise. If the latter holds, according to the assumption that the K coherent groups are uncorrelated to each other,R x has K linearly independent rows.

Remark 1.
The parameter m plays a significant role in rank restoration since it determines the array's DOF and effective aperture. We consider two extreme cases to discuss the choices of m. If m = M − 2P + 3 − P max 2 , i.e., m achieves its upper bound, then one barely restores the rank deficiency with the minimum number of columns of B(i), but this may result in some signals being undetected, especially at low SNRs or for few snapshots. On the other hand, If m = N + 1, i.e., m achieves its lower bound, then one restores the rank deficiency with redundant columns of B(i), but this may cause biased estimates due to the noise subspace being one-dimensional. Based on our simulation results, the proposed method performs well when m lies between these two bounds, away from the extremes. However, the optimum choice of m is still an open problem and is beyond the scope of this paper.
When m > N and 2q ≥ P K , the non-singularity of the smoothed covariance matrixR can be recovered, then the MUSIC algorithm can be employed to estimate the DOAs of coherent signals. More specifically, one can perform the eigen-decomposition ofR as where Σ s is a diagonal matrix consisting of the N largest eigenvalues, and Σ n is a diagonal matrix consisting of the m − N smallest eigenvalues. The columns of U s are the eigenvectors corresponding to the N largest eigenvalues, while the columns of U n are the eigenvectors corresponding to the m − N smallest eigenvalues. According to the subspace principle, the columns of U s span the signal subspace, which coincides with the range space of A, and the signal subspace is orthogonal to the noise subspace spanned by the columns of U n . Thus, This indicates that the DOAs of the coherent signals can be obtained by finding the N peaks from the spatial spectrum function So far, we have described the proposed algorithm for coherent signal DOA estimation in the presence of unknown mutual coupling. The major steps of the proposed algorithm are summarised as follows: Step 1. Obtain L snapshots of the received signal x(t) at t = t 1 , t 2 , · · · , t L , and form the following matrix as: Step 2. Calculate the covariance matrix using the above data matrix by Step 3. Select the middle array and calculate corresponding covariance matrixR x according to Equation (6); Step 4. Utilise rows ofR x to construct the matrix B(i), i = 1, 2, · · · , M − 2P + 2, according to Equation (10); Step 5. Construct the overall rank-restored matrixR as Equation (14); Step 6. Perform the eigen-decomposition ofR, and obtain the noise subspace U n ; Step 7. Scan the direction over [−90 • , 90 • ] with a step size of 0.1 • . Calculate the spacial spectrum using Equation (18), and obtain the DOA estimates θ 1 , θ 2 , · · · , θ N .
Although the mutual coupling information for a sensor array is a cumbersome calculation using electromagnetic simulators and measurement in practice, for some circumstances, such as high precision applications in satellite navigation and airborne early warning radar, it has to be obtained and compensated for a-priori. In the presence of known mutual coupling, we should first eliminate the effects of mutual coupling before applying the matrix reconstruction for rank restoration. To be specific, the following steps have to be carried out: Step 1. Perform eigen-decomposition of R x , one has SinceŨ s and CAΓ span the same signal subspace, there holdsŨ s = CAΓT where T is a nonsingular matrix.
Step 2. Reconstruct a covariance matrix by compensating for mutual coupling effects as Then, following the similar Steps 4-7 in the unknown mutual coupling case one can resolve the coherent signals.

Remark 2.
Once the mutual coupling information is known, one can compensate the coupling effects by Equation (22), and then the additional information of the auxiliary sensors, which had to be abandoned in the case of unknown mutual coupling, can now be exploited. Compared with Equation (14) there are more matrices {B(i)} M i=1 available after mutual coupling compensation, which means that more DOF and a larger effective aperture can be achieved, and improved estimation performance can thus be expected. In fact, our method can handle up to M − P max 2 coherent signals in the presence of known mutual coupling. The identifiability of DOA estimation for the proposed and comparative methods with unknown mutual coupling will be discussed in Section 3.3.

Separable Signal Number
The coherent signals are resolved with an increased effective array aperture compared with previous techniques. This makes better use of the DOF of the original ULA and allows more signals to be resolved. This motivates us to discuss the number of separable signals of the proposed method. Compared with the standard spatial smoothing (SS) technique, which estimates the coherent signals by averaging subarrays, under unknown mutual coupling [14], the maximum number of signals estimated by our method can be relaxed. If m > N and 2q ≥ P max , the proposed method can estimate a maximum number of M − 2P + 2 − P max 2 coherent signals since M − 2P + 2 ≥ N + P max 2 , while the standard spatial smoothing can estimate at most M − 2P + 2 − P max coherent signals since P max subarrays are required and M − 2P + 3 − P max > N. As we can see, the key factor influencing the number of resolvable coherent signals is the maximum number of signals in one coherent group, in other words, the larger P max is, the better our approach performs. Our approach also has better DOA estimation accuracy than SS as will be demonstrated by the simulation results in the next section. The method by Mao et al. [15], referred to as matrix reconstruction in data (MRD), can also achieve better DOA estimation accuracy than SS, provided that the effective array aperture is the same, but this comes at a cost of the number of signals which can be resolved. It can resolve at most M−2P+1 2 coherent signals and the array aperture is fixed regardless of the number of coherent sources. The spatial smoothing techniques are more flexible in this aspect as the number of subarrays and effective array aperture can be varied. Both the algorithms in [14,15] also ignore the conjugate of the received data in rank restoration, which significantly restricts the number of resolvable signals in practice. Table 1 lists the minimum number of array elements required to resolve a given number of signals by the three methods. For simplicity, we assume that each group has the same number of coherent signals. We can see that our method can use less array elements than the other two algorithms, to estimate the same number of signals. Table 1. Minimum number of array elements required. 2  1  2  2  6  7  5  2  2  2  4  8  11  7  3  1  4  4  12  13  10  3  2  3  6  13  17  12  4  1  6  6  18  19  15  4  2  4  8  18  23  16  4  2  6  12  24  31  21  3  3  6  18 28 41 23

Simulation Results and Discussion
In this section, a series of numerical experiments under different conditions are conducted to examine the performance of the proposed method. Simulations are carried out for a 14-element ULA with half-wavelength spacing between adjacent elements. For simplicity, we assume that all coherent signals have identical power σ 2 s , and the input SNR is defined as 10 log 10 (σ 2 s /σ 2 n ). SS, MRD, and the proposed method with known mutual coupling dealt with in Equation (22) are chosen for comparison. The mutual coupling is assumed can be negligible at a distance larger than 1.5λ. Hence, P = 3 and the mutual coupling coefficients are assumed to be c 0 = 1, c 1 = −0.1545 + 0.4755j, c 2 = 0.1618 − 0.1176j. The accuracy of the DOA estimate is measured from 1000 Monte Carlo runs in terms of the root mean square error (RMSE) which is defined as is the estimate of θ i for the n-th trial, and N is the number of coherent signals. Additionally, to assess the overall reliability of all the algorithms, the probability of resolution is defined as where F is the number of trials, and F r is the number of successful estimations for which the absolute DOA estimation errors are within 1 • . Here, for fair comparison it is first assumed that the three algorithms have an identical aperture after rank restoration, i.e., m = 5. From this figure, we can see that as the SNR and the number of snapshots increase, the RMSE of DOA estimation decreases gradually for all of the methods. The proposed method is significantly superior to MRD and SS, especially at low SNRs and few snapshots, since more column vectors have been utilised for rank restoration with the same array aperture, and is only inferior to the one when mutual coupling is known. Although MRD slightly outperforms SS, both have poor DOA estimation accuracy at low SNRs and few snapshots, and there is a clear discrepancy between these two and the proposed method even for snapshot sizes up to 1000 and beyond at 0 dB as shown in Figure 3b. Figure 4 depicts the probability of resolution of successful estimation versus input SNR and the number of snapshots for the same scenario as in Figure 3. It can be observed that all methods attain a 100% successful estimation above 11 dB. As the SNR decreases, the probability of successful estimation starts dropping for each method at a certain point, known as SNR threshold, until it eventually becomes zero. With unknown mutual coupling, the proposed method has the lowest SNR threshold followed by MRD and then SS which has the highest SNR threshold. The performance of the proposed method when the mutual coupling is known is also shown for comparison. Evidently the proposed method is superior to the counterparts. When the SNR is fixed at 0 dB, our approach can achieve 100% probability with unknown mutual coupling above 1800 snapshots. Even if the number of snapshots is greater than 2000, the MRD and SS algorithms can barely reach 50% success probability. Consistent with the illustration in Figure 3, MRD has a slightly better performance than SS in terms of probability of resolution versus SNR and the number of snapshots.

DOA Estimation of Coherent Signals From One Group
The simulations above demonstrate the performance of all algorithms with an identical aperture m = 5. As discussed in Section 3.3, the effective array aperture after rank restoration by MRD is always halved, in other words, the effective array aperture is fixed and cannot adapt to changes in the number of signals. In contrast, both the proposed and SS algorithm are adaptive to complex situations, such as multiple coherent groups, allowing a trade off between array aperture and rank restoration. To evaluate this advantage, we set m = 6 for the proposed and SS algorithm and compare them with the fixed array aperture MRD algorithm, i.e., m = 5, while keeping the other configurations the same as the first scenario.    Figures 5 and 6 show that compared with Figures 3 and 4 the performance of both our approach and SS meliorates, and our approach is still superior to SS and MRD and approaches to the case with known mutual coupling, while the SS method has the largest improvement in the accuracy and reliability. This is mainly because the effective array aperture is improved while keeping the rank completely restored. This group of simulations fully demonstrates that the proposed and SS algorithm have a clear advantage over the MRD algorithm, and the fixed aperture imposed by MRD is not optimal for rank restoration.  ], respectively. The total number of signals is equal to half the elements in the middle array, and thus MRD fails to work. To enable the proposed algorithm and SS work, we select m = 6, 7 for rank restoration.

DOA Estimation of Coherent Signals From Multiple Groups
Based on the simulation settings, the results of the RMSE versus SNR and the number of snapshots are depicted in Figure 7. As shown in this figure, the proposed approach performs the best over the whole range of SNR and the number of snapshots values for coherent signal estimation. When m = 6, the proposed algorithm with known or unknown mutual coupling is inferior to the first scenario. One possible explanation is that the fading coefficients in the two scenarios do not have same level of impact for rank restoration in general. SS has larger estimation errors than our method as it only exploits of 4 subarrays for rank restoration while our method effectively utilises twice as many. When m = 7, both the proposed and SS algorithm achieve a better performance since the limitation of effective array aperture is relaxed. Although the SS algorithm improves significantly for m = 7, it is still notably worse than our proposed method.
Next, in Figure 8 we plot the probability of resolution of the two methods by varying the SNR and total number of snapshots. It can be seen that the proposed approach still outperforms the SS algorithms for both configuration m = 6, 7. The performance of the SS method improves considerably when m = 7 but is strictly worse than our method.

Conclusions
This paper addresses the problem of DOA estimation when unknown mutual coupling and multipath coexist. An efficient approach has been presented that first circumvents the unknown mutual coupling effect using the output of the selected middle array, and then reconstructs the covariance matrix of the middle array for rank restoration. The proposed method has a better DOA estimation accuracy than the existing methods and can resolve more coherent signals. This is achieved due to two main factors. Firstly replacing standard spatial smoothing in the rank restoration step by a new matrix reconstruction algorithm takes advantage of the information in the off the main diagonal entries of the spatial covariance matrix. Secondly the new matrix reconstruction algorithm for restoring the rank of the coherent signal subspace exploits the conjugate information of the entries of the covariance matrix. These are ignored by standard spatial smoothing techniques and can effectively improve the array aperture or equivalently improve the rank restoration. Extensive simulation results demonstrate the validity and efficiency of the proposed method.