Parameter Estimation for Two-Dimensional Incoherently Distributed Source with Double Cross Arrays

Abstract

A direction of arrival (DOA) estimator for two-dimensional (2D) incoherently distributed (ID) sources is presented under proposed double cross arrays, satisfying both the small interval of parallel linear arrays and the aperture equalization in the elevation and azimuth dimensions. First, by virtue of a first-order Taylor expansion for array manifold vectors of parallel linear arrays, the received signal of arrays can be reconstructed by the products of generalized manifold matrices and extended signal vectors. Then, the rotating invariant relations concerning the nominal elevation and azimuth are derived. According to the rotating invariant relationships, the rotating operators are obtained through the subspace of the covariance matrix of the received vectors. Last, the angle matching approach and angular spreads are explored based on the Capon principle. The proposed method for estimating the DOA of 2D ID sources does not require a spectral search and prior knowledge of the angular power density function. The proposed DOA estimation has a significant advantage in terms of computational cost. Investigating the influence of experimental conditions and angular spreads on estimation, numerical simulations are carried out to validate the effectiveness of the proposed method. The experimental results show that the algorithm proposed in this paper has advantages in terms of estimation accuracy, with a similar number of sensors and the same experimental conditions when compared with existing methods, and that it shows a robustness in cases of model mismatch.

1. Introduction

In view of the problem of localization, traditional parameter estimation algorithms are based on a point source model. When spatial scattering characteristics of targets can be ignored and propagations are supposed to be a single straight path between targets and receive arrays, the point target model can simplify the calculation. However, in practical applications, due to scattering, reflection, diffraction and refraction in a complex environment, a large number of multipath phenomena result in the signal source expanding at a certain angle in space, which often has more complex spatial distribution characteristics than the point signal source, and distributed source models are presented in this context [1]. The distribution source exists in two situations. First, the size of the target cannot be ignored when compared with the distance. At this time, the reflection intensity of different parts of the target is random, and the direction is also variable. For instance, in an underwater sonar background, with the reduction of distance, the spatial geometry of a target cannot be ignored, many parts of the target reflecting signals; the point target model cannot describe the characteristics of the source effectively, and distributed sources models are suitable in this context. Secondly, due to the complex environment near the target, local strong scattering occurs. For example, with wireless communications in an environment with high-rise buildings, the intensity of the reflected signal changes with the orientation [2]. There are also similar scenarios, such as the area target or clutter involved in radar signal processing, whose array observation signal usually has a certain diffusion in the spatial and temporal domain [3,4]. At this time, the parameters estimation of targets based on a point source model assumption is seriously deteriorated, while the algorithms based on a distributed source model can obtain more accurate information.

The distributed signal source model has a broad application prospect in radar, sonar and mobile communication because it adopts the signal model, which is more in line with the actual signal characteristics. A distributed source can be thought of as an aggregation of point sources, which are also called scatterers, in a spatial range with a certain density distribution. According to the dependence of scatterers on the time domain, sources are divided into incoherently distributed (ID) sources and coherently distributed (CD) sources. The spatial distribution characteristics of ID and CD sources are described by the angular power density function (ADPF) and angular signal distribution function (ASDF), respectively. According to the spatial distribution characteristics of the distributed sources, APDF or ASDF can be expressed as a variety of distribution functions, such as Gaussian, uniform and asymmetric distributions [5]. Sources can be one-dimensional (1D) and two-dimensional (2D), depending on the spatial dimension of the sources. Whatever distribution function a source has, the APDF or ASDF of a 1D distributed source is described by two parameters: nominal angles, which can be called direction of arrival (DOA), and angular spreads, which reflect the spatial extension of a source. The APDF or ASDF of a 2D distributed source is depicted by four parameters including the nominal azimuth and nominal elevation, which can be collectively called DOA, and the elevation spread and azimuth spread, which can be called angular spreads reflecting the 2D spatial extension of a source. Compared with the 1D source, which has more parameters, the complexity of estimation for 2D sources significantly increases. In this paper, we are concerned with 2D ID sources.

Regarding point sources, the DOA estimation has recently made some important achievements in the field of mixed signal and multiple input multiple output (MIMO) systems, as well as sparse arrays. Under a uniform rectangular array, the 2D DOA estimation of circular and non-circular mixed signals is proposed in [6]. The authors of [7] have proposed an adaptive beamforming algorithm based on the minimum variance distortionless response principle. With coprime arrays, the authors of [8] have presented a virtual array interpolation estimator; the authors of [9] have proposed a sparse reconstruction using a proposed sliding window scheme for optimization; the authors of [10] have proposed an ESPRIT framework estimator; the authors of [11] have proposed an estimator utilizing Toeplitz matrix functions to resolve off-grid sources. An estimator was proposed in [12] for 2D DOA, based on the least squares technique in the electromagnetic vector sensor MIMO system. The authors of [13] have proposed a DOA estimator based on off-grid sparse Bayesian learning. In [14], the authors have proposed an estimator where transmitted and received DOA are first estimated via trilinear decomposition and refined by 2D searches.

Considering CD sources, scholars have presented estimators in [15,16,17,18,19,20,21,22,23,24,25]. For 1D ID distributed sources, scholars have proposed various parameter estimation algorithms, mainly including subspace algorithms such as DSPE [1] and DSPARE [4], covariance matching estimation techniques [26,27,28], maximum likelihood estimation algorithms [29,30,31,32] and beamforming algorithms [33,34,35]. For 2D ID sources, the authors of [36] have proposed a 2D covariance matching estimation technique, which uses four-dimensional nonlinear optimization with a large amount of computation. Scholars have proposed low-complexity estimators for 2D ID sources. The authors of [37] have proposed an estimator for DOAs where nominal elevations were first obtained through rotating invariant relations, after which nominal azimuths were solved by a Capon spectral search, which increased computational complexity.

In this paper, aiming at the parameters of 2D ID sources, an estimator is proposed under double cross arrays. Based on the rotating invariant relations derived through the first order of the Taylor expansion of array manifold vectors, received vectors of arrays are reconstructed via products of extended signal vectors and generalized manifold matrices. Nominal angles can be obtained through signal subspaces of covariance matrices of received vectors. According to the solved nominal angles, angular spreads can be searched via the 2D Capon principle. The main contributions of this article are listed below:

• Generally, the parameters of ID sources constitute an approximate conclusion under the condition of smaller sensor spacing. We propose a double cross array that satisfies both the small space between sensors and the aperture equalization in the elevation and azimuth dimensions.
• Based on the reconstructed received signal vectors and deduced rotating invariant relations, we propose an approach to a solution for nominal angles according to the ESPRIT framework, and we propose an angle matching as well as angular spreads solutions using the Capon estimation.
• The proposed method has an advantage in the DOA estimation with respect to the computational cost, shows advantages in terms of estimation accuracy with a similar number of sensors, and shows robustness in the case of a model mismatch.

The structure of the paper is as follows. Section 2 elaborates the source model and the received vectors of proposed arrays. In Section 3, based on rotating invariant relations and the reconstructed received signal vectors, approaches to a solution for DOA and angular spreads are detailed. Section 4 illustrates the simulation results, which are discussed. Section 5 draws the conclusions.

Notations: Scalar variables are denoted by italic letters, and vectors and matrices are denoted by bold letters. (•)⁻¹, (•)*, (•)^T and (•)^H means inverse, complex conjugate, transpose and Hermitian transposition of a matrix. E[•], (./) and (•)+ denote expectation, element-wise division and pseudo-inverse operations. [•]_k is the kth element of a vector. angle(•) is a phase of a complex number.

2. 2D ID Source Model and Double Cross Arrays

The proposed double cross arrays structure is shown in Figure 1. The array is composed of two pairs of parallel linear arrays X1 and X2, and Z1 and Z2. All subarrays are located on the xoz plane. Subarrays X1 and Z1 are located on the x axis and z axis, respectively. Subarrays X1 and Z1 share the sensor located on the origin. Both subarrays X1 and X2 have M sensors on both sides of the z axis, while both subarrays Z1 and Z2 have M sensors on both sides of the x axis. All linear arrays have 2M + 1 sensors with a spacing of d meters. The intervals between each pair of parallel linear arrays are all δ.

Assuming that k narrow-band 2D ID sources in the far field with nominal angles (θ_i, φ_i) (i = 1, 2, …, k) are impinging on the array. θ_i and φ_i denote the nominal elevation and nominal azimuth of the ith source, respectively. θ_i ∈ [0, π], φ_i ∈ [0, π]. Then, the received vectors of the arrays can be written as:

$$\{\begin{array}{l}{x}_m(t)={\displaystyle \sum _{i=1}^k{\displaystyle \iint {\mathsf{\alpha }}_m(\theta ,\phi )s_i(\theta ,\phi ,t)d\theta }d\phi }+n_{xm}(t)\\ z_m(t)={\displaystyle \sum _{i=1}^k{\displaystyle \iint {\mathsf{\beta }}_m(\theta ,\phi )s_i(\theta ,\phi ,t)d\theta }d\phi }+n_{zm}(t)\end{array},$$

(1)

where m = 1 or 2. x₁(t) and x₂(t) represent the received vectors of arrays X1 and X2, respectively; z₁(t) and z₂(t) are the received vectors of arrays Z1 and Z2. n_xm(t) and n_zm(t) are additive white Gaussian noises which are unrelated to the signals. The noise power is ${\sigma }_n^2$. α₁ (θ, φ) and α₂ (θ, φ) represent array manifold vectors of the point source with respect to X1 and X2. β₁ (θ, φ) and β₂ (θ, φ) denote array manifold vectors of the point source with respect to Z1 and Z2.

$$\{\begin{array}{l}{\mathsf{\alpha }}_1(\theta ,\phi )={[e^{j2\pi d(1-M)\cos \theta \sin \phi /\lambda },\cdots 1,\cdots ,e^{j2\pi (M-1)d\cos \theta \sin \phi /\lambda }]}^T\\ {\mathsf{\alpha }}_2(\theta ,\phi )={\mathsf{\alpha }}_1(\theta ,\phi )e^{j2\pi \delta \cos \phi /\lambda }\end{array},$$

(2)

$$\{\begin{array}{l}{\mathsf{\beta }}_1(\theta ,\phi )={[e^{j2\pi d(1-M)\cos \phi /\lambda },\cdots ,1,\cdots ,e^{j2\pi (M-1)d\cos \phi /\lambda }]}^T\\ {\mathsf{\beta }}_2(\theta ,\phi )={\mathsf{\beta }}_1(\theta ,\phi )e^{j2\pi \delta \cos \theta \sin \phi /\lambda }\end{array}.$$

(3)

s_i(θ,φ,t) is the complex random angular signal density of the ith distributed source, representing the scatterer intensity of the source from the direction (θ, φ) at the time t. Unlike a point source, the signal of a distributed source exists not only in a single direction (θ_i, φ_i) but in a spatial distribution around (θ_i, φ_i). An ID source means that different scatterers from one target generate uncorrelated signals. Therefore, one direction from s_i(θ,φ,t) is uncorrelated with other directions from s_i(θ,φ,t), which means that the following relationship exists:

$$E[s_i(\theta ,\phi ,t)s_i^{*}({\theta }^{\prime },{\phi }^{\prime },t)]={\sigma }_i^{2}f(\theta ,\phi ;u_i)\delta (\theta -{\theta }^{\prime })\delta (\phi -{\phi }^{\prime }),$$

(4)

where ${\sigma }_i^2$ is the power of the ith source, δ(•) is the Kronecker delta function and $p(\theta ,\phi ;u_i)$ is the normalized angular power density function (APDF). APDF have a parameter set u_i= [${\theta }_i$, ${\varphi }_i$, ${\sigma }_{\theta i}$, ${\sigma }_{\varphi i}$] denoting the nominal azimuth, nominal elevation, azimuth spread and elevation spread, respectively. $p(\theta ,\phi ;u_i)$ satisfies the following relationship:

$${\displaystyle \iint p(\theta ,\phi ;u_i)}d\theta d\phi =1.$$

(5)

3. Proposed Method

This part is composed of five sections. First, generalized manifold vectors and the rotating invariant relations of generalized manifold vectors are derived. Next, the received signal vectors can be reconstructed as products of generalized signal vectors and generalized manifold vectors. Based on the reconstructed form of the received signal vectors and rotating invariant relations, nominal angles can be resolved separately according to an ESPRIT-like framework. Then, in order to pair the nominal angles, an angle matching approach is proposed under the Capon principle. Afterwards, angular spreads can be obtained by two-dimensional spectral searching using the Capon estimation. Last, the computational procedure is introduced, and the complexity of the proposed approach in comparison with several existing methods is analyzed.

3.1. Rotating Invariant Relations

The array manifold vectors of arrays X1 and X2 are expanded with a first-order Taylor series at the point (θ_i, φ_i) as follows:

$${\alpha }_1(\theta ,\phi )\approx {\alpha }_1({\theta }_i,{\phi }_i)+[{\alpha }_1{({\theta }_i,{\phi }_i]}_{\theta }^{\prime }(\theta -{\theta }_i)+[{\alpha }_1{({\theta }_i,{\phi }_i]}_{\phi }^{\prime }(\phi -{\phi }_i),$$

(6)

$${\alpha }_2(\theta ,\phi )\approx {\alpha }_2({\theta }_i,{\phi }_i)+[{\alpha }_2{({\theta }_i,{\phi }_i]}_{\theta }^{\prime }(\theta -{\theta }_i)+[{\alpha }_2{({\theta }_i,{\phi }_i]}_{\phi }^{\prime }(\phi -{\phi }_i),$$

(7)

where ${[•]}_{\theta }^{\prime }$ and ${[•]}_{\phi }^{\prime }$ represent the first partial derivative of the function with respect to ${\theta }_i$ and ${\phi }_i$, respectively. Thus, we have the following relationships:

$$\{\begin{array}{l}{\alpha }_2({\theta }_i,{\phi }_i)={\alpha }_1({\theta }_i,{\phi }_i)e^{j2\pi \delta \cos {\phi }_i/\lambda }\\ {[{\alpha }_2({\theta }_i,{\phi }_i]}_{\theta }^{\prime }={[{\alpha }_1({\theta }_i,{\phi }_i)]}_{\theta }^{\prime }{e}^{j2\pi \delta \cos {\phi }_i/\lambda }\\ [{\alpha }_2{({\theta }_i,{\phi }_i]}_{\phi }^{\prime }={[{\alpha }_1({\theta }_i,{\phi }_i)]}_{\phi }^{\prime }{e}^{j2\pi \delta \cos {\phi }_i/\lambda }-(j2\pi \delta \sin {\phi }_i/\lambda ){\alpha }_1({\theta }_i,{\phi }_i)e^{j2\pi \delta \cos {\phi }_i/\lambda }\end{array}.$$

(8)

If $\delta /\lambda \ll 1$, the second item on the right side of $[{\alpha }_2{({\theta }_i,{\phi }_i]}_{\phi }^{\prime }$ is negligible, and thus the following approximation relationship exists:

$$[{\alpha }_2{({\theta }_i,{\phi }_i]}_{\phi }^{\prime }\approx {[a_1({\theta }_i,{\phi }_i)]}_{\phi }^{\prime }{e}^{j2\pi \delta \cos {\phi }_i/\lambda }.$$

(9)

Define the extended signal vectors $\overline{s}=[{\overline{s}}_1,{\overline{s}}_2,{\overline{s}}_3]$^H and $\tilde{s}=[{\overline{s}}_1,{\overline{s}}_3]$^H. ${\overline{s}}_1,{\overline{s}}_2$ and ${\overline{s}}_3$ can be written as follows:

$$\{\begin{array}{l}{\overline{s}}_1=[{\rho }_{10},{\rho }_{20},\cdots ,{\rho }_{k0}]\\ {\overline{s}}_2=[{\rho }_{1\theta },{\rho }_{2\theta },\cdots ,{\rho }_{k\theta }]\\ {\overline{s}}_3=[{\rho }_{1\phi },{\rho }_{2\phi },\cdots ,{\rho }_{k\phi }]\end{array},$$

(10)

where ${\rho }_{i0}$, ${\rho }_{i\theta }$ and ${\rho }_{i\phi }$ can be expressed as follows:

$$\{\begin{array}{l}{\rho }_{i0}={\displaystyle \iint s_i(\theta ,\phi ,t)d\theta }d\phi \\ {\rho }_{i\theta }={\displaystyle \iint (\theta -{\theta }_i)s_i(\theta ,\phi ,t)d\theta }d\phi \\ {\rho }_{i\phi }={\displaystyle \iint (\phi -{\phi }_i)s_i(\theta ,\phi ,t)d\theta }d\phi \end{array}.$$

(11)

$\overline{s}$ is a 3k × 1 dimensional vector, and $\tilde{s}$ is a 2k × 1 dimensional vector. It can be proven that the following relationship exists [37]:

$$E[{\rho }_{il}{\rho }_{in}]=\{\begin{array}{ll}{\sigma }_i^2& l=n=0\\ {\sigma }_i^2{M}_{\theta i}& l=n=\theta \\ {\sigma }_i^2{M}_{\phi i}& l=n=\phi \\ 0& l\ne n\end{array},$$

(12)

where $M_{\theta i}$ and $M_{\phi i}$ can be expressed as follows:

$$M_{\theta i}={\displaystyle \iint {(\theta -{\theta }_i)}^2{p}_i(\theta ,\phi ,t)d\theta }d\phi ,$$

(13)

$$M_{\phi i}={\displaystyle \iint {(\phi -{\phi }_i)}^2{p}_i(\theta ,\phi ,t)d\theta }d\phi .$$

(14)

Since the sources are unrelated, we have:

$$E[\overline{s}{\overline{s}}^H]=diag(\Lambda ,M_{\theta }\Lambda ,M_{\phi }\Lambda ),$$

(15)

$$E[\tilde{s}{\tilde{s}}^H]=diag(\Lambda ,M_{\phi }\Lambda ),$$

(16)

where:

$$\{\begin{array}{l}\Lambda =diag({\sigma }_1^2,{\sigma }_2^2,\cdots ,{\sigma }_k^2)\\ M_{\theta }=diag(M_{\theta 1},M_{\theta 2},\cdots ,M_{\theta k})\\ M_{\phi }=diag(M_{\phi 1},M_{\phi 2},\cdots ,M_{\phi k})\end{array}.$$

(17)

The generalized manifold matrix of the subarray X1 is defined as $[A_{11},A_{12},A_{13}]$. Blocks of the matrix can be expressed as follows:

$$\{\begin{array}{l}{A}_{11}=[{\alpha }_1({\theta }_1,{\phi }_1),{\alpha }_1({\theta }_2,{\phi }_2),\cdots ,{\alpha }_1({\theta }_k,{\phi }_k)]\\ A_{12}=[{[{\alpha }_1({\theta }_1,{\phi }_1)]}_{\theta }^{\prime },{[{\alpha }_1({\theta }_2,{\phi }_2)]}_{\theta }^{\prime },\cdots ,{[{\alpha }_1({\theta }_k,{\phi }_k)]}_{\theta }^{\prime }]\\ A_{13}=[{[{\alpha }_1({\theta }_1,{\phi }_1)]}_{\phi }^{\prime },{[{\alpha }_1({\theta }_2,{\phi }_2)]}_{\phi }^{\prime },\cdots ,{[{\alpha }_1({\theta }_k,{\phi }_k)]}_{\phi }^{\prime }]\end{array}.$$

(18)

Apparently, $[A_{11},A_{12},A_{13}]$ is a (2M + 1) × 3k dimensional generalized manifold matrix. The received vector of the X1 array can be written as:

$$x_1(t)=[A_{11},A_{12},A_{13}]\overline{s}+n_{x1}(t).$$

(19)

The generalized manifold matrix of the subarray X2 is defined as $[A_{21},A_{22},A_{23}]$, which is (2M + 1) × 3k dimensional. Blocks of the matrix can be written as follows:

$$\{\begin{array}{l}{A}_{21}=[{\alpha }_2({\theta }_1,{\phi }_1),{\alpha }_2({\theta }_2,{\phi }_2),\cdots ,{\alpha }_2({\theta }_k,{\phi }_k)]\\ A_{22}=[{[{\alpha }_2({\theta }_1,{\phi }_1)]}_{\theta }^{\prime },{[{\alpha }_2({\theta }_2,{\phi }_2)]}_{\theta }^{\prime },\cdots ,{[{\alpha }_2({\theta }_k,{\phi }_k)]}_{\theta }^{\prime }]\\ A_{23}=[{[{\alpha }_2({\theta }_1,{\phi }_1)]}_{\phi }^{\prime },{[{\alpha }_2({\theta }_2,{\phi }_2)]}_{\phi }^{\prime },\cdots ,{[{\alpha }_2({\theta }_k,{\phi }_k)]}_{\phi }^{\prime }]\end{array}.$$

(20)

The received vector of the X2 array can be written as:

$$x_2(t)=[A_{21},A_{22},A_{23}]\overline{s}+n_{x2}(t).$$

(21)

According to Equations (8) and (9), the following rotating invariant relation can be obtained:

$$[A_{21},A_{22},A_{23}]\approx [A_{11}\Phi ,A_{12}\Phi ,A_{13}\Phi ],$$

(22)

where the rotating operator can be written as:

$$\Phi =diag(e^{j2\pi \delta \cos {\phi }_1/\lambda },e^{j2\pi \delta \cos {\phi }_2/\lambda },\cdots ,e^{j2\pi \delta \cos {\phi }_k/\lambda }).$$

(23)

The array manifold vectors of arrays Z1 and Z2 are expanded with the first-order Taylor series at the point (θ_i, φ_i) as follows:

$${\mathsf{\beta }}_1(\theta ,\phi )\approx {\mathsf{\beta }}_1({\theta }_i,{\phi }_i)+{[{\mathsf{\beta }}_1({\theta }_i,{\phi }_i)]}_{\phi }^{\prime }(\phi -{\phi }_i),$$

(24)

$${\mathsf{\beta }}_2(\theta ,\phi )\approx {\mathsf{\beta }}_2({\theta }_i,{\phi }_i)+{[{\mathsf{\beta }}_2({\theta }_i,{\phi }_i)]}_{\theta }^{\prime }(\theta -{\theta }_i)+{[{\mathsf{\beta }}_2({\theta }_i,{\phi }_i)]}_{\phi }^{\prime }(\phi -{\phi }_i),$$

(25)

where:

$${[{\mathsf{\beta }}_2({\theta }_i,{\phi }_i)]}_{\theta }^{\prime }={\mathsf{\beta }}_1({\theta }_i,{\phi }_i)(-j2\pi \delta \sin {\theta }_i\sin {\phi }_i/\lambda )e^{j2\pi \delta \cos {\theta }_i\sin {\phi }_i/\lambda },$$

(26)

$${[{\mathsf{\beta }}_2({\theta }_i,{\phi }_i)]}_{\phi }^{\prime }={[{\mathsf{\beta }}_1({\theta }_i,{\phi }_i)]}_{\phi }^{\prime }{e}^{j2\pi \delta \cos {\theta }_i\sin {\phi }_i/\lambda }+(j2\pi \delta \cos {\theta }_i\cos {\phi }_i/\lambda )e^{j2\pi \delta \cos {\theta }_i\sin {\phi }_i/\lambda }{\mathsf{\beta }}_1({\theta }_i,{\phi }_i).$$

(27)

If $\delta /\lambda \ll 1$, the right side of Equation (26) and the second item on the right side of Equation (27) can be ignored, then we have the following relationships:

$$\{\begin{array}{l}{\mathsf{\beta }}_2({\theta }_i,{\phi }_i)={\mathsf{\beta }}_1({\theta }_i,{\phi }_i)e^{j2\pi \delta \cos {\theta }_i\sin {\phi }_i/\lambda }\\ {[{\mathsf{\beta }}_2({\theta }_i,{\phi }_i)]}_{\phi }^{\prime }\approx {[{\mathsf{\beta }}_1({\theta }_i,{\phi }_i)]}_{\phi }^{\prime }{e}^{j2\pi \delta \cos {\theta }_i\sin {\phi }_i/\lambda }\end{array}.$$

(28)

Define the generalized manifold matrix of array Z1 as $[B_{11},B_{12}]$, which is (2M + 1) × 2k dimensional. Blocks in the matrix can be written as follows:

$$\{\begin{array}{l}{B}_{11}=[{\mathsf{\beta }}_1({\theta }_1,{\phi }_1),{\mathsf{\beta }}_1({\theta }_2,{\phi }_2),\cdots ,{\mathsf{\beta }}_1({\theta }_k,{\phi }_k)]\\ B_{12}=[{[{\mathsf{\beta }}_1({\theta }_1,{\phi }_1)]}_{\phi }^{\prime },{[{\mathsf{\beta }}_1({\theta }_2,{\phi }_2)]}_{\phi }^{\prime },\cdots ,{[{\mathsf{\beta }}_1({\theta }_k,{\phi }_k)]}_{\phi }^{\prime }]\end{array}.$$

(29)

Then, the received vector of the Z1 array can be written as:

$$z_1(t)=[B_{11},B_{12}]\tilde{s}+n_{z1}(t).$$

(30)

Define the generalized manifold matrix of array Z2 as $[B_{21},B_{22}]$, which is (2M + 1) × 2k dimensional. Blocks in the matrix can be written as follows:

$$\{\begin{array}{l}{B}_{21}=[{\mathsf{\beta }}_2({\theta }_1,{\phi }_1),{\mathsf{\beta }}_2({\theta }_2,{\phi }_2),\cdots ,{\mathsf{\beta }}_2({\theta }_k,{\phi }_k)]\\ B_{22}=[{[{\mathsf{\beta }}_2({\theta }_1,{\phi }_1)]}_{\phi }^{\prime },{[{\mathsf{\beta }}_2({\theta }_2,{\phi }_2)]}_{\phi }^{\prime },\cdots ,{[{\mathsf{\beta }}_2({\theta }_k,{\phi }_k)]}_{\phi }^{\prime }]\end{array}.$$

(31)

The received vector of subarray Z2 can be written as:

$$z_2(t)=[B_{21},B_{22}]\tilde{s}+n_{z2}(t).$$

(32)

According to (28), the following rotating invariant relation is obtained:

$$[B_{21},B_{22}]\approx [B_{11}\Psi ,B_{12}\Psi ],$$

(33)

where the rotating operator can be written as:

$$\Psi =diag(e^{j2\pi \delta \cos {\theta }_1\sin {\phi }_1/\lambda },e^{j2\pi \delta \cos {\theta }_2\sin {\phi }_2/\lambda },\cdots ,e^{j2\pi \delta \cos {\theta }_k\sin {\phi }_k/\lambda }).$$

(34)

3.2. Estimation of Nominal Angles

Combine the received vectors of subarrays X1 and X2 as follows:

$$x_{12}(t)=\left[\begin{array}{l}{x}_1(t)\\ x_2(t)\end{array}\right]=[A_{x1},A_{x2},A_{x3}]\overline{s}+n_x(t),$$

(35)

where $[A_{x1},A_{x2},A_{x3}]$ is the generalized manifold matrix of the combined received vectors x₁₂(t), which can be written as:

$$\left[A_{x1},A_{x2},A_{x3}\right]=\left[\begin{array}{l}{A}_{11},A_{12},A_{13}\\ A_{21},A_{22},A_{23}\end{array}\right].$$

(36)

The combined noise vector of subarrays X1 and X2 can be written as:

$$n_x(t)=\left[\begin{array}{l}{n}_{x1}(t)\\ n_{x2}(t)\end{array}\right].$$

(37)

The covariance matrix of the combined receive vectors x₁₂(t) can be written as:

$$R_x^{12}=E[x_{12}(t)x_{12}^H(t)].$$

(38)

The sample covariance matrix ${\widehat{R}}_x^{12}$ can be substituted for $R_x^{12}$ in calculation, having the following expression:

$${\widehat{R}}_x^{12}=\frac{1}{N}{\displaystyle \sum _{t=1}^N{x}_{12}(t)x_{12}^H(t)}.$$

(39)

For ID sources, the rank of the covariance matrix is theoretically larger than the source number. Define $E_x$ as a subspace formed by the eigenvectors corresponding to the maximum k eigenvalues of the covariance matrix $R_x^{12}$. As both M_θi and M_φi are numbers far below 1 with the condition of a small angular spread, the subspace $E_x$ is the same as the subspace spanned by $A_{x1}$. Then, there exists a k × k nonsingular matrix T that satisfies the following relationship:

$$E_x=\left[\begin{array}{l}{A}_{11}\\ A_{11}\Phi \end{array}\right]T,$$

(40)

where $\Phi$ is the rotating operator described by Equation (23), A₁₁ is the first block of the generalized manifold matrix of the subarray X1 described by Equation (18).

Define $E_{x1}$ and $E_{x2}$ as the upper and lower kth rows of $E_x$, then we have:

$$E_{x1}=A_{11}T.$$

(41)

$$E_{x2}=A_{21}T.$$

(42)

Thus, the following relationship can be obtained:

$$E_{x2}=E_{x1}{T}^{-1}\Phi T.$$

(43)

Define ${\Omega }_x=T^{-1}\Phi T$, we have:

$${\Omega }_x=E_{x1}^{+}{E}_{x2}.$$

(44)

Therefore, the nominal elevations of the ID sources can be solved as follows:

$${\phi }_i=\arccos \frac{angle({\eta }_i)}{2\pi \delta /\lambda }i=1,2,\cdots ,k,$$

(45)

where ${\eta }_i$ is the ith eigenvalue of ${\Omega }_x$.

Combine the received vectors of subarrays Z1 and Z2 as follows:

$$z_{12}(t)=\left[\begin{array}{l}{z}_1(t)\\ z_2(t)\end{array}\right]=\left[B_{z1},B_{z2}\right]\tilde{s}+n_z(t),$$

(46)

where the generalized manifold matrix of the combined received vectors z₁₂(t) can be written as:

$$[B_{z1},B_{z2}]=\left[\begin{array}{l}{B}_{11},B_{12}\\ B_{21},B_{22}\end{array}\right].$$

(47)

The combined noise vector of subarrays Z1 and Z2 can be written as:

$$n_z(t)=\left[\begin{array}{l}{n}_{z1}(t)\\ n_{z2}(t)\end{array}\right].$$

(48)

The covariance matrix of the combined received vectors z₁₂(t) can be expressed as:

$$R_z^{12}=E[z_{12}(t)z_{12}^H(t)].$$

(49)

The sample covariance matrix ${\widehat{R}}_z^{12}$ can be substituted for $R_z^{12}$ in calculation, having the following expression:

$${\widehat{R}}_z^{12}=\frac{1}{N}{\displaystyle \sum _{t=1}^N{z}_{12}(t)z_{12}^H(t)}.$$

(50)

As M_φi is a number below 1 under the condition of a small angular spread, the subspace $E_z$ constructed by eigenvectors corresponding to the maximum k eigenvalues of the covariance matrix $R_z^{12}$ is the same as the subspace spanned by $B_{z1}$. Then, there exists a k × k nonsingular matrix Q that satisfies the following relationship:

$$E_z=\left[\begin{array}{l}{B}_{11}\\ B_{11}\Psi \end{array}\right]Q,$$

(51)

where $\Psi$ is the rotating operator described by Equation (34), B₁₁ is the first block of the generalized manifold matrix of the subarray Z1 described by Equation (29).

Define $E_{z1}$ and $E_{z2}$ as the upper and lower kth rows of $E_z$:

$$E_{z1}=B_{11}Q,$$

(52)

$$E_{z2}=B_{21}Q.$$

(53)

Thus, we can obtain:

$$E_{z2}=E_{z1}{Q}^{-1}\Psi Q.$$

(54)

Define ${\Omega }_z=Q^{-1}\Psi Q$, we have:

$${\Omega }_z=E_{z1}^{+}{E}_{z2}.$$

(55)

Then, the nominal azimuths of the ID sources can be obtained as follows:

$${\theta }_i=\arccos \frac{angle({\mu }_i)}{2\pi \delta /\lambda \sin {\phi }_i}i=1,2,\cdots ,k,$$

(56)

where ${\mu }_i$ is the ith eigenvalue of ${\Omega }_z$.

3.3. Angle Matching

Considering the combined generalized manifold matrix $[A_{x1},A_{x2},A_{x3}]$ of parallel subarrays X1 and X2, the Capon estimation criterion has the following expression:

$$\begin{array}{ll}\min & w^H{R}_x^{12}w\\ s.t.& w^H[A_{x1},A_{x2},A_{x3}]=1\end{array}$$

(57)

where $R_x^{12}$ is the covariance matrix of the combined received vectors x₁₂(t) described by Equation (38), and w is an unknown coefficient vector.

The above equation can be solved by the Lagrange function. The Capon cost function of subarrays X1 and X2 can be obtained as:

$$L(\theta ,\phi )=\frac{1}{{\Vert {[A_{x1},A_{x2},A_{x3}]}^H{({\widehat{R}}_x^{12})}^{-1}[A_{x1},A_{x2},A_{x3}]\Vert }_F}.$$

(58)

The steps for angle matching can be summarized as follows:

• Choose ${\widehat{\phi }}_i={\phi }_i$ as the estimated nominal elevation from the set $\{{\phi }_1,{\phi }_2,\cdots ,{\phi }_k\}$, traverse the set $\{angle({\mu }_1),angle({\mu }_2),\cdots ,angle({\mu }_k)\}$ and solve the paired nominal angles $({\theta }_j,{\widehat{\phi }}_i)$ $(j=1,2,\cdots ,k)$ according to Equation (56).
• Substitute $({\theta }_j,{\widehat{\phi }}_i)$ $(j=1,2,\cdots ,k)$ into (58) to select the pair that makes the maximum value of $L(\theta ,\phi )$ and denote the pair as $({\widehat{\theta }}_j,{\widehat{\phi }}_i)$.
• Remove ${\widehat{\phi }}_i$ from the set $\{{\phi }_1,{\phi }_2,\cdots ,{\phi }_k\}$ and remove $2\pi \delta /\lambda \cos {\widehat{\theta }}_i\sin {\widehat{\phi }}_i$ from the set $\{angle({\mu }_1),angle({\mu }_2),\cdots ,$ $angle({\mu }_k)\}$.
• Repeat steps 1–3.

All nominal angles can be paired by (k + 2)(k − 1)/2 calculations.

3.4. Estimation of Angular Spread

Considering subarray X1, the covariance matrix of the received vector has the following expression:

$$R(\theta ,\phi ;u_i)={\displaystyle \iint p(\theta ,\phi ;u_i){\alpha }_1(\theta ,\phi )}{\alpha }_1^H(\theta ,\phi )d\theta d\phi .$$

(59)

For Gaussian 2D ID sources, APDF has the following expression:

$$p(\theta ,\phi ;u_i)=\frac{1}{2\pi {\sigma }_{\theta i}{\sigma }_{\phi i}}\exp \left\{-0.5\left[{\left(\frac{\theta -{\theta }_i}{{\sigma }_{\theta i}}\right)}^2+{\left(\frac{\phi -{\phi }_i}{{\sigma }_{\phi i}}\right)}^2\right]\right\}.$$

(60)

Assuming $d/\lambda =0.5$, and within the premise of small angular spreads, the elements of $R(\theta ,\phi ;u_i)$ can be expressed as follows:

$${[R(\theta ,\phi ;u_i)]}_{lh}=e^{j\pi (l-h)\cos {\theta }_i\sin {\phi }_i}{e}^{-0.5\left\{{[\pi {\sigma }_{\theta i}(l-h)\sin {\theta }_k\sin {\phi }_i]}^2+{[\pi {\sigma }_{\phi i}(l-h)\cos {\theta }_i\cos {\phi }_i]}^2\right\}},$$

(61)

where [•]_lh represents the element of the lth row and hth column of a matrix.

For uniform 2D ID sources, APDF has the following expression:

$$p(\theta ,\phi ;u_i)=\{\begin{array}{ll}\frac{1}{4{\sigma }_{\theta i}{\sigma }_{\phi i}}& \left|\theta -{\theta }_i\right|\le {\sigma }_{\theta i}\mathrm{and}\left|\phi -{\phi }_i\right|\le {\sigma }_{\phi i}\\ 0& \left|\theta -{\theta }_i\right|\ge {\sigma }_{\theta i}\mathrm{or}\left|\phi -{\phi }_i\right|\ge {\sigma }_{\phi i}\end{array}.$$

(62)

$R(\theta ,\phi ;u_i)$ has the following expression:

$${[R(\theta ,\phi ;u_i)]}_{lh}=e^{j\pi (l-h)\cos {\theta }_i\sin {\phi }_i}\frac{\sin [\pi {\sigma }_{\phi i}(l-h)\cos {\theta }_i\cos {\phi }_i]·\sin [\pi {\sigma }_{\theta i}(l-h)\sin {\theta }_i\sin {\phi }_i]}{\pi {\sigma }_{\phi i}(l-h)\cos {\theta }_i\cos {\phi }_i·\pi {\sigma }_{\theta i}(l-h)\sin {\theta }_i\sin {\phi }_i}.$$

(63)

Azimuth spreads as well as elevation spreads can be obtained by two-dimensional spectral searching using the Capon estimation:

$$({\widehat{\sigma }}_{\theta },{\widehat{\sigma }}_{\phi })=\arg \max \frac{1}{{\sigma }_{\max }\{{({\widehat{R}}_x^1)}^{-1}R(\theta ,\phi ;u_i)\}},$$

(64)

where ${\sigma }_{\max }\{•\}$ represents the maximum eigenvalue of the matrix. ${\widehat{R}}_x^1$ is the sample covariance matrix of arrays X1, having the following expression:

$${\widehat{R}}_x^1=\frac{1}{N}{\displaystyle \sum _{t=1}^N{x}_1(t)x_1^H(t)}.$$

(65)

3.5. Algorithm Procedure and Complexity Analysis

According to the above analysis, the algorithm steps are summarized as follows:

• Calculate the sample covariance matrix ${\widehat{R}}_x^{12}$ and ${\widehat{R}}_z^{12}$.
• Obtain $E_x$ and $E_z$ by decomposition of ${\widehat{R}}_x^{12}$ and ${\widehat{R}}_z^{12}$. Calculate ${\Omega }_x$ and ${\Omega }_z$ according to (44) and (55).
• Conduct the decomposition of ${\Omega }_x$ and ${\Omega }_z$ to obtain the eigenvalues ${\eta }_i$ and ${\mu }_i$, and then calculate the nominal elevations ${\phi }_i$ and nominal azimuths ${\theta }_i$ from (45) and (56).
• Determine the nominal angles $({\widehat{\theta }}_i,{\widehat{\phi }}_i)$ $(j=1,2,\cdots ,k)$ according to the angle matching method.
• Substitute $({\widehat{\theta }}_i,{\widehat{\phi }}_i)$ $(j=1,2,\cdots ,k)$ into $R(\theta ,\phi ;u_i)$ according to Equation (61) or (63). Solve the angular spreads parameters ${\widehat{\sigma }}_{\theta i}$ and ${\widehat{\sigma }}_{\phi i}$ from (64).

The flow of the algorithm is shown in Figure 2.

The complexity of the algorithm proposed in this paper with respect to nominal angles contains four parts. First, calculating the covariance matrix ${\widehat{R}}_x^{12}$ and ${\widehat{R}}_z^{12}$ is O(NK²), where N is the snapshots number and K is the total sensor number. The second part is the eigendecomposition of ${\widehat{R}}_x^{12}$ and ${\widehat{R}}_z^{12}$, which is O(K³). The third is the eigendecomposition of ${\Omega }_x$ and ${\Omega }_z$, which is O(k³). The fourth is angle matching, which is O(k²). The computational cost of reference [36] mainly contains calculating sample covariance matrices O(NK²) and the alternating projection technique, which require four-dimensional optimizations: O(K⁴). The complexity of the estimation method of double parallel arrays proposed in reference [37] includes three parts. First, calculating the covariance matrix is O(NK²). The second part is the eigendecomposition of the covariance matrix, which is O(K³). Third, the 1D search is O(kL_φK²), where L_φ constitutes the search points of the nominal azimuth. The computational complexity of the angular spread parameters is mainly divided into two parts: the inversion of the covariance matrix ${\widehat{R}}_x^1$, which is O[(K/4)³], and the 2D search to obtain the angular spread parameters, which is O[(K/4)³kLσ_θLσ_φ], where Lσ_θ and Lσ_φ are the search points of the elevation spread and azimuth spread, respectively. As reference [37], which deals with the problem of DOA estimation, does not involve the estimation of angular spreads,

Table 1. Computational complexity of different methods for DOA estimation.

Method	Sample Covariance Matrix	Eigendecomposition/Alternating Projection Technique	Searching	Total
Proposed	O(NK2)	O(K3) + O(k3)		O(NK2) + O(K3) + O(k3)
Reference [36]	O(NK2)	O(K4)		O(NK2) + O(K4)
Reference [37]	O(NK2)	O(K3)	O(kLφK2)	O(NK2) + O(K3) + O(kLφK2)

shows the main computational costs of three methods for DOA estimation. From Table 1, we can clearly obtain that the computational cost for the DOA estimation of the proposed algorithm is lower than for the algorithms in refs. [36,37].

It should be noted that reference [36] also deals with the estimation of angular spreads, where the computational cost for angular spreads is involved in a four-dimensional optimization, O(K⁴). From a previous analysis, the computational cost of the proposed method for angular spreads is O[(K/4)³kLσ_θLσ_φ]. Generally speaking, angular spreads are within a few degrees for 2D ID sources. Supposing both Lσ_θ and Lσ_φ are set with 20 search points, which means that the search range is 2° with a step size 0.1°, when K > 6.25k the computational cost for angular spreads in reference [36] is larger than the proposed one. To summarize, the proposed DOA estimation has a significant advantage in terms of computational cost, whereas the angular spreads are not obvious.

4. Results and Discussions

The arrays structure is set with each linear array containing five sensors, which means that the total sensor number of the double cross arrays is 19. The array element spacing is set as $d=\lambda /2$ and the parallel array interval is $\delta =\lambda /10$. The root mean square error (RMSE) of the DOA is defined as RMSE_a with following expression:

$$RMS{E}_a=\sqrt{\frac{1}{Mc}{\displaystyle \sum _{\varsigma }^{Mc}{({\widehat{\theta }}^{\varsigma }-\theta )}^2+\frac{1}{Mc}{\displaystyle \sum _{\varsigma }^{Mc}{({\widehat{\phi }}^{\varsigma }-\phi )}^2}}}.$$

(66)

MC represents the Monte Carlo simulation times, while the superscript ^ and variable $\varsigma$ of ${\widehat{\theta }}^{\varsigma }$ and ${\widehat{\phi }}^{\varsigma }$ represent the estimated value in the $\varsigma$th Monte Carlo simulation. RMSE of the angular spread is defined as RMSE_s with following expression:

$$RMS{E}_s=\sqrt{\frac{1}{Mc}{\displaystyle \sum _{\varsigma }^{Mc}{({\widehat{\sigma }}_{\theta }^{\varsigma }-{\sigma }_{\theta })}^2+\frac{1}{Mc}{\displaystyle \sum _{\varsigma }^{Mc}{({\widehat{\sigma }}_{\phi }^{\varsigma }-{\sigma }_{\phi })}^2}}},$$

(67)

where ${\widehat{\sigma }}_{\theta }^{\varsigma }$ and ${\widehat{\sigma }}_{\phi }^{\varsigma }$ represent the estimated azimuth spread and elevation spread in the $\varsigma$th Monte Carlo simulation. RMSE_a denotes the estimation errors of DOA, while RMSE_s denotes the estimation errors of the angular spreads. The signal-to-noise ratio (SNR) is defined as 10log(1/${\sigma }_n^2$), where the noise is assumed to be the Gaussian white zero mean with a variance of ${\sigma }_n^2$.

In the first experiment, we investigate the influence of SNR and the snapshots number on the estimation. A 2D ID Gaussian source with parameters of [30°, 45°, 2°, 3°] is set to be estimated. The experiment is completed by 100 independent Monte Carlo simulations. Figure 3 shows the estimated RMSE_a and RMSE_s with SNR from 0 to 30 db when the snapshots number is equal to 200. Figure 4 shows the estimated RMSE_a and RMSE_s when changing the snapshots number when SNR is equal to 20 db. Figure 3 and Figure 4 also show the estimation results by the method used in refs. [36,37], using 20 sensors based on a uniform rectangular array (URA) and double parallel linear arrays (DPLA), respectively, as well as the Cramer–Rao bound (CRB) [36]. It can be seen that with the improvement of SNR or the snapshots number, the estimation accuracy of the proposed algorithm in this paper is better than that of the algorithm proposed in refs. [36,37] under the same experimental conditions. The method used in reference [37] unitizes double parallel arrays X1and X2 with sensor spacing in the elevation dimension of less than one tenth of the impinging signal wavelength and sensor spacing in the azimuth dimension of less than half of the impinging signal wavelength; consequently, the physical aperture of the array in the elevation dimension is far less than that of the azimuth dimension. Therefore, during the iterative estimation, the estimation of the nominal elevation will produce larger errors, which will also be transmitted to the solution of the nominal azimuth. Though the URA of the method used in reference [36] has a structural balance in the elevation and azimuth dimensions, the spacing of sensors of the double cross array proposed in both the elevation and azimuth dimensions is less than half of the impinging wavelength when the number of elements is constant; the aperture of the array proposed in this paper is larger than the URA in the two dimensions, exhibiting a higher estimation accuracy. To summarize, the advantage of the proposed double cross arrays is that they satisfy both the premise of a small interval of parallel linear arrays and the aperture equalization in the elevation and azimuth dimensions when compared with URA and DPLA.

In the second experiment, we examine the influence of the angular spread on estimation. For a 2D ID source with nominal angles (60°, 35°), two kinds of APDF (Gaussian and uniform) are considered, respectively. The azimuth spreads and elevation spreads are supposed to be equivalent in both kinds of sources. The estimated RMSE_a with an angular spread from 0° to 10° is investigated. The experiment is completed by 100 independent Monte Carlo simulations with the SNR at 20 db and snapshots number at 200. Figure 5 shows that the estimation accuracy decreases with an increase of the angular spreads. For the Gaussian source, when the angular spreads are 5°, RMSE_a is 0.03. As the angular spreads reach 10°, RMSE_a is 0.4. For the uniform source, RMSE_a is 0.04 and 0.7 when the angular spreads are 5° and 10°, respectively. These experimental results show that the estimation performance decreases with the increase of the angular spreads, but that it is still satisfactorily within 10°. The proposed algorithm is robust in the case of small angular spreads.

As can be seen from Figure 5, the DOA estimation accuracy is highest when the angular spreads are 0°. As a matter of fact, whether it is a Gaussian or uniform distributed source or any other distributed sources, when the angular spreads are 0°, the parameters describing the ID source only include the DOA parameters, which means that the distributed source can be equal to the point source. This experiment also shows that the method proposed in this paper is not only suitable for the DOA estimation of the ID source but also suitable for the DOA estimation of the point source. In the case of a mismatch between the point source and ID source model, it also has good robustness.

In the third experiment, we investigate the ability of the algorithm to estimate multiple distributed sources. Three ID sources, A, B and C, are set to be estimated. A and B are ID sources with a Gaussian APDF, C is an ID source with a uniform APDF. The three sources have the following parameters: [55°, 80°, 1.5°, 2°], [25°, 70°, 2°, 3°] and [35°, 40°, 3.5°, 4.5°]. SNR is set at 20 db, and the snapshots number is 200. The experiment is completed by 100 independent Monte Carlo simulations. Figure 6 shows the 100 results of the estimated DOAs of the three sources. Figure 5 illustrates that the proposed method can estimate nominal angles of multiple 2D ID source simultaneously without prior knowledge of APDF. Figure 7 shows the 100 results of the estimated angular spreads of the three sources, respectively. From the derivation process in Parts 3.2 and 3.3, as well as from the experimental results, we can see that the method proposed in this paper does not require knowledge of the specific form of APDF in the DOA estimation. As the angular spreads are parameters of the APDF function, it is impossible to discuss the angular spreads parameter without the APDF function. The method proposed in this paper does not need to know the APDF in the DOA estimation, but solving the angular spreads parameters requires prior information on the APDF form first, followed by searching and optimizing.

It should be noted that the rotating invariant relations in this paper are based on the assumption that the interval of the parallel linear array satisfies a small spacing. In practice, if the spacing of sensors is small, on the one hand, when taking into account the characteristics of the sensors, the mutual coupling effect between the sensors cannot be ignored. On the other hand, the installation error of sensors will also affect the estimation effect. Therefore, this is difficult to apply in some high frequency signal fields such as mobile communication, as the frequency of mobile communication rises from 800 MHz to 2.5 GHz. However, in the low-frequency sonic detection of the underwater background, the sonar frequency can be reduced to 100 HZ, and the signal wavelength can reach approximately 14.5 m. Similarly, with respect to a very high frequency (VHF) radar, the general working wavelength is selected to be within 1–3 m. In these cases, the interval of the parallel linear array can actually reach a value that satisfies the hypothesis.

5. Conclusions

In this paper, an estimation of the parameters of a 2D ID source is presented based on a double cross array that satisfies both the small interval of parallel linear arrays and the aperture equalization in the elevation and azimuth dimensions. The rotating invariant relations of the generalized manifold vectors with respect to nominal angles are derived under the proposed array. Received vectors are reconstructed as products of generalized manifold matrices and extended signal vectors by taking the first-order Taylor expansion of the array manifold vector. DOAs are calculated on the basis of the covariance matrix of the received vectors according to a proposed ESPRIT-like framework. Angular spreads are obtained by a two-dimensional search. Numerical simulations are conducted to verify the proposed method while considering different experimental conditions, angular spreads and multiple sources. The outcomes show that the proposed method shows a better estimation performance under the same experimental condition and sensor number when compared with existing methods.