Estimation of Two-Dimensional Non-Symmetric Incoherently Distributed Source with L-Shape Arrays

Abstract

In the field of array signal processing, distributed sources can be regarded as an assembly of point sources within a spatial distribution. In this study, a two-dimensional (2D) non-symmetric incoherently distributed (ID) source model is proposed; we explore the estimation of a 2D non-symmetric ID source using L-shape arrays. The 2D non-symmetric ID source is established by modeling the angular power density function (APDF) as a Gaussian mixture model. Estimation of the non-symmetric distributed source is proposed based on the expectation maximization (EM) framework. The proposed EM iterative framework contains three steps in the process of each circle. Firstly, the nominal azimuth and nominal elevation of each Gaussian component are obtained from the phase parts of elements in sample covariance matrices. Then the angular spreads can be solved through a one-dimensional (1D) search by the original generalized Capon estimator. Finally, weights of each Gaussian component are obtained by solving the least-squares estimator. Simulations are conducted to verify the effectiveness of the estimation technique.

1. Introduction

In array signal processing, applications such as wireless communications, radar and underwater acoustics, point source models (assuming that signals impinging into receive arrays are from point sources) are commonly used, which can simplify calculations. In the real surroundings of radar and sonar systems, because of multipath propagation between receive arrays and targets, especially when the distances of targets and receive arrays are short, the spatial scatterers of targets cannot be ignored, assumptive condition of point source is no longer valid and point source models cannot characterize sources effectively, which should be described by distributed source models [1]. Distributed sources can be regarded as an assembly of point sources within a spatial distribution. The shape of spatial distribution is related to geometry and surface property of a target, for instance, in underwater detection. Considering multipath propagation and the surface feature of targets, distributed source models are more appropriate in near field of radar or sonar detection.

Distributed sources may be classified into coherently distributed (CD) sources or incoherently distributed (ID) sources [1]. Characterized by deterministic angular signal distribution function (ASDF) representing spatial distribution, CD sources assume that scatterers within a source are coherent. On the contrary, scatterers of an ID source are supposed to be uncorrelated. Spatial distributions of ID sources are characterized by angular power density function (APDF).

Spatial distribution of a distributed source, both ASDF and APDF, can be generally modeled as Gaussian, uniform or any other distribution. Parameters of ASDF and APDF contain nominal angles and angular spreads. Nominal angles can also be expressed as nominal direction-of-arrival (DOA) representing the center of targets. Nominal spreads denote the spatial extension of targets. Traditional estimation of a distributed source mainly involves estimation of DOA and angular spreads, which are implied in ASDF or APDF.

For CD sources, representative achievements of parameters estimation are deriving the rotation invariance relation with respect to different array configuration by virtue of Taylor series expansions [2,3,4,5,6,7,8,9]. Performance of distributed signal parameter estimator (DSPE) algorithm is analyzed in [10] and the performance of multiple signal classification (MUSIC) is analyzed in [11]. In this paper, the modeling and estimation of an ID source are considered.

Investigators for ID sources have been developed from classical point sources estimation techniques. DSPE [1] and dispersed signal parametric estimation (DISPARE) [12] have been developed from MUSIC, both of which involve a two-dimensional (2D) spectral search. In [13], an estimator extended from estimation of signal parameters via rotational invariance techniques (ESPRIT) has been proposed for ID sources, which use the total least square-ESPRIT (TLS-ESPRIT) algorithm to estimate the nominal DOAs of sources firstly, and then the angular spreads are estimated using the central moments of the distribution. Generalizations of Capon’s methods have been proposed in [14,15,16], which involve a two-dimensional (2D) spectral search besides a high-order matrix inversion. The maximum likelihood (ML) approach [17,18,19] has better accuracy but leads to a multidimensional nonlinear optimization requiring high computational complexity. Developed from least squares estimators, covariance matching estimation techniques (COMET) [20,21,22,23,24] have lower computational complexity than ML but with the same large sample behavior. Applying sparse representation to first-order Taylor expansion of steering vectors, in the case of small angular spreads, the authors of [25] proposed an estimator via block sparse Bayesian learning for multiple incoherently distributed sources, which has presented better accuracy under fewer snapshots.

The aforementioned methods consider sources as one-dimensional (1D) ID models, which have two parameters: the nominal DOA and angular spread. However, impinging signals and arrays sensors are not in the same plane practically. 2D distributed source models characterized by 2D DOA and 2D angular spread should be more reasonable, which usually contain four parameters: nominal azimuth, nominal elevation, azimuth angular spread and elevation angular spread. Including more parameters, there have been relatively few studies on estimation of 2D ID sources. The authors of [26] have proposed an extension of COMET for 2D ID sources, which employs alternating projection principle [27] and estimates 2D DOAs and 2D angular spreads separating source powers and noise variances to reduce complexity. Based on double parallel uniform linear arrays, the authors of [28] have proposed a TLS-ESPRIT like approach for 2D DOAs, where estimation of nominal elevations is firstly implemented via rotational invariance relations based on the diffused steering vectors derived from the first-order Taylor series expansions, then the nominal azimuth is estimated by 1D searching. Authors of [29] have proposed a (ESPRIT) like algorithm for 2D ID sources based on uniform rectangular arrays. Receive vectors of arrays in [29] are described by generalized steering vector, which is composed of nominal steering vector and its first-order partial derivatives along the direction of elevation and azimuth. Nominal angles are estimated through rotational invariance relations of generalized steering vectors, which are derived under the assumption of small angular spreads and small distance between sensors. Generally, DOA estimation of ID sources by virtue of diffused steering vectors or generalized steering vectors are derived from the first-order Taylor series expansions under the assumption of small angular spreads [26] or the assumption of both small angular spreads and small distance between sensors [28,29]. Though different methods based on different arrays have analyzed the influence of angular spreads of Gaussian or uniform ID sources on DOA estimation, what is a satisfactory accuracy respect to a real specific target has not been carried out until now. The assumption that distance between sensors is far less than wavelength of the signal is applicable to few specific detection areas, such as underwater low frequency detection.

All the distributed sources estimation techniques mentioned above assume that the shape of sources is symmetric. Nevertheless, scatterers are distributed irregularly within targets and the surface of targets is generally non-symmetric. Models of 1D non-symmetric ID sources have been proposed according to the principle that a non-symmetric distribution can be constructed by symmetric distributions. The authors of [30] have proposed a non-symmetric ID source model, where the shape of APDF can be figured via the variation of the ratio between two Gaussian distributions. The authors of [31] have employed the Gaussian mixture model to characterize the 1D non-symmetric ID sources and proposed a COMET algorithm embedded in expectation maximization (EM) framework [32,33] to estimate the 1D non-symmetric APDF. Containing more parameters, estimation for a non-symmetric ID source mainly lies in estimation of the non-symmetric APDF, which is different from estimation of a symmetric ID source.

To the best of our knowledge, there are no algorithms for 2D non-symmetric distributed sources. In this paper, we are concerned with modeling and estimation of a 2D non-symmetric ID source. As the principle that symmetric distributions can form a non-symmetric distribution is also true in the case of 2D, we present a 2D non-symmetric ID source model by constituting APDF with the 2D Gaussian mixture model. The authors of [34] have explored Cramer-Rao bound of L-shape arrays composed of two orthogonal ULAs and shown that such arrays have better accuracy potential than traditional cross arrays using same number of sensors with respect to point sources. Utilizing the propagator method, the authors of [35] have said that L-shape arrays composed of two orthogonal ULAs with less number of elements can estimate better than parallel shape arrays using the same method. Several DOA estimation methods [36,37,38,39] have been developed based on a point source model with L-shaped array placed in xoz or xoy plane. Estimation for a 2D non-symmetric ID source is proposed under EM framework based on L-shape arrays. The general EM framework is to get the best parameters by maximization of the likelihood function during the iteration process. As a 2D non-symmetric ID source has more parameters than those of an ID symmetric source; maximization of the likelihood function makes it difficult to apply for a 2D non-symmetric ID source because of nonlinear and high dimensional property. In our method, parameters are obtained successively in each EM cycle. Deducing the covariance matrices via the first-order Taylor series expansion of the steering vectors in each Gaussian component, we find that the nominal DOA parameters are related to the phase parts in elements of covariance matrices. Accordingly, nominal azimuth and nominal elevation of each Gaussian component can be obtained by the sample covariance matrices. Then angular spreads can be searched through the original Generalization of Capon’s estimator base on DOA parameters. Weights are estimated by a least squares fit of theoretical covariance and sample covariance.

The rest of this paper is organized as follows. Section 2 formulates a 2D non-symmetric ID source model and describes the received signal vectors and covariance matrix under L-shape array. Section 3 details the proposed algorithm of estimation for parameters of the 2D non-symmetric ID source. Section 4 shows numerical simulations to validate the proposed estimation method. Section 5 draws the conclusion of this paper.

2. Distributed Source Model

Figure 1 shows the L-shape arrays configuration, which uses the xoy plane. Array X is composed of sensors on x axis, while sensors on y axis constitute array Y. Each linear array consists of K sensor elements separated by d meters, and the two linear arrays share an origin sensor. Suppose that there is a stationary narrow-band 2D ID source with azimuth angle θ and elevation angle φ distributed in a spatial distribution. θ ∈ [0, π/2], φ ∈ [0, π/2]. λ is the wavelength of the signal impinging into the array sensors. The K × 1 dimensional received signal vectors of the arrays X and Y can be expressed as follows:

$$x(t)={\left[x_1(t),x_2(t),\cdots ,x_k(t),\cdots ,x_K(t)\right]}^T,$$

(1)

$$y(t)={\left[y_1(t),y_2(t),\cdots ,y_k(t),\cdots ,y_K(t)\right]}^T,$$

(2)

where (•)^T denotes the transpose. x_k(t) and y_k(t) are the signal received by kth sensors in arrays X and Y, which can be expressed as follows:

$$x_k(t)=s(t){\displaystyle \sum _{l=1}^L{\alpha }_l(t)}{e}^{j2\pi d(k-1)\cos {\theta }_l\sin {\phi }_l/\lambda }+n_{xk}(t),$$

(3)

$$y_k(t)=s(t){\displaystyle \sum _{l=1}^L{\alpha }_l(t)}{e}^{j2\pi d(k-1)\sin {\theta }_l\sin {\phi }_l/\lambda }+n_{yk}(t).$$

(4)

An ID source means that different scatterers from one target generate uncorrelated multipath signals. s(t) in Equations (3) and (4) is the impinging signal into the target. L is number of scatterers. n_xk(t) and n_yk(t) are noises received. α_l(t) is random complex ray gain of the lth scatterer. The random complex ray gain α_l(t) is supposed to be white and independent from snapshot to snapshot as well as from scatterer to scatterer, which has the following relationship:

$$\{\begin{array}{l}E[{\alpha }_l(t){\alpha }_l(t^{\prime })]=0\\ E[{\alpha }_l(t){\alpha }_{l^{\prime }}^{\ast }(t^{\prime })]=\frac{{\sigma }_{\alpha }}{L}\delta (l-l^{\prime })\delta (t-t^{\prime })\end{array},$$

(5)

where |α_l(t)|² = σ_α/L, δ(•) is the Kronecker delta function, (•)^* denotes the conjugate operator. Define n_x(t) and n_y(t), which are the K × 1 dimensional additive noise vectors of arrays X and Y; they can be written as:

$$n_x(t)=\left[n_{x1}(t),n_{x2}(t),\cdots ,n_{xk}(t),\cdots ,n_{xK}(t)\right],$$

(6)

$$n_y(t)=\left[n_{y1}(t),n_{y2}(t),\cdots ,n_{yk}(t),\cdots ,n_{yK}(t)\right].$$

(7)

n_x(t) and n_y(t) can be combined into

$$n(t)=\left[\begin{array}{l}{n}_x(t)\\ n_y(t)\end{array}\right].$$

(8)

The noise is assumed to be uncorrelated with signal and uncorrelated between sensors as well as Gaussian white with zero mean:

$$E[n(t)n^H(t\prime )]=\rho I_{2K}\delta (t-t\prime ),$$

(9)

$$E[n(t)n^T(t\prime )]=0\forall t,t\prime ,$$

(10)

where ρ is the power of noise, I_2K is the 2K × 2K identity matrix, (•)^H denotes the Hermitian transpose.

Combine x(t) and y(t) into:

$$z(t)=\left[\begin{array}{l}x(t)\\ y(t)\end{array}\right].$$

(11)

The signal is assumed to be uncorrelated with the noise; the covariance matrix of z(t) can be written as follows:

$$R_z=E[z(t)z^H(t)]=\left[\begin{array}{cc}{R}_x& R_{xy}\\ R_{yx}& R_x\end{array}\right]=\left[\begin{array}{cc}E[x(t)x^H(t)]& E[x(t)y^H(t)]\\ E[y(t)x^H(t)]& E[y(t)y^H(t)]\end{array}\right].$$

(12)

The concept of APDF reflecting the distribution of scatterers of a source can be traced back to references [1] and [19], which can be approximately modeled as 2D Gaussian and uniform or any other distribution function according to the different characteristic of ID sources. Denote f(θ,φ) as APDF of an ID source. The (k,h)th noise free element of covariance matrix R_x is given by (proof can be seen in Appendix A):

$$\begin{array}{c}{\left[R_x\right]}_{kh}=E[x_k(t)x_h^{\ast }(t)]={\left|s(t)\right|}^{2}E\left[{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{l^{\prime }=1}^L{\alpha }_l(t)e^{j2\pi d(k-1)\cos {\theta }_l\sin {\phi }_l/\lambda }{\alpha }_{l^{\prime }}^{\ast }(t)e^{-j2\pi d(h-1)\cos {\theta }_{l^{\prime }}\sin {\phi }_{l^{\prime }}/\lambda }}}\right]\\ =P{\displaystyle \iint e^{j2\pi d(k-h)\cos \theta \sin \phi /\lambda }f(\theta ,\phi )}d\theta d\phi \end{array}.$$

(13)

where P = σ_α|s(t)|² is the received power of the target. Define a(θ,φ) and b(θ,φ) are K × 1 dimensional steering vectors of arrays X and Y, which can be written as follows:

$$a(\theta ,\phi )={\left[1,e^{j2\pi d\cos \theta \sin \phi /\lambda },\cdots e^{j2\pi (K-1)d\cos \theta \sin \phi /\lambda }\right]}^T,$$

(14)

$$b(\theta ,\phi )={\left[1,e^{j2\pi d\sin \theta \sin \phi /\lambda },\cdots e^{j2\pi (K-1)d\sin \theta \sin \phi /\lambda }\right]}^T.$$

(15)

Thus, the covariance matrices R_x, R_y, R_xy and R_yx can be expressed as follows (proof is in Appendix B):

$$\{\begin{array}{l}{R}_x=P{\displaystyle \iint f(\theta ,\phi )}a(\theta ,\phi )a^H(\theta ,\phi )d\theta d\phi +\rho I_K\\ R_y=P{\displaystyle \iint f(\theta ,\phi )}b(\theta ,\phi )b^H(\theta ,\phi )d\theta d\phi +\rho I_K\\ R_{xy}=P{\displaystyle \iint f(\theta ,\phi )}a(\theta ,\phi )b^H(\theta ,\phi )d\theta d\phi \\ R_{yx}=P{\displaystyle \iint f(\theta ,\phi )}b(\theta ,\phi )a^H(\theta ,\phi )d\theta d\phi \end{array}.$$

(16)

In this study, the Gaussian mixture model is used for angular power density of the distributed source in order to express the non-symmetric distribution, so the APDF of a non-symmetric ID source can be expressed as follows:

$$f(\theta ,\phi )\approx {\displaystyle \sum _{i=1}^q{w}_{i}g(\theta ,\phi ;u_i)}={\displaystyle \sum _{i=1}^q{w}_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\},$$

(17)

where the APDF consists of q Gaussian components g(θ,φ; u_i) i= 1,2,…,q. u_i = [${\theta }_i$,${\varphi }_i$,${\sigma }_{\theta i}$,${\sigma }_{\varphi i}$] is the parameter set of the ith Gaussian component denoting the nominal azimuth, nominal elevation, azimuth angular spread, and elevation angular spread, respectively.

To normalize f(θ,φ), the weighting coefficient w_i satisfies the following constraint:

$${\displaystyle \sum _{i=1}^q{w}_i}=1.$$

(18)

APDF in (17) can be shaped asymmetrically via the variation of the weighting coefficient w_i and parameter set of each Gaussian component.

Combine a(θ,φ) and b(θ,φ) into:

$$c(\theta ,\phi )=\left[\begin{array}{l}a(\theta ,\phi )\\ b(\theta ,\phi )\end{array}\right].$$

(19)

The summation of Gaussian mixture distributions in Equation (17) can be substituted for f(θ,φ). R_z can be expressed respectively as follows:

$$R_z\approx {\displaystyle \sum _{i=1}^q{P}_i}{\displaystyle \iint g(\theta ,\phi ;u_i)c(\theta ,\phi )c^H(\theta ,\phi )d\theta }d\phi +\rho I_{2K},$$

(20)

Equation (16) can be expressed as follows:

$$\{\begin{array}{l}{R}_x\approx {\displaystyle \sum _{i=1}^q{P}_i{\displaystyle \iint g(\theta ,\phi ;u_i)a(\theta ,\phi )a^H(\theta ,\phi )d\theta }d\phi }+\rho I_K\\ R_y\approx {\displaystyle \sum _{i=1}^q{P}_i{\displaystyle \iint g(\theta ,\phi ;u_i)b(\theta ,\phi )b^H(\theta ,\phi )d\theta }d\phi }+\rho I_K\\ R_{xy}\approx {\displaystyle \sum _{i=1}^q{P}_i{\displaystyle \iint g(\theta ,\phi ;u_i)a(\theta ,\phi )b^H(\theta ,\phi )d\theta }d\phi }\\ R_{yx}\approx {\displaystyle \sum _{i=1}^q{P}_i{\displaystyle \iint g(\theta ,\phi ;u_i)b(\theta ,\phi )a^H(\theta ,\phi )d\theta }d\phi }\end{array},$$

(21)

where I_K denotes the K × K identity matrix,$P_i=w_{i}P$ denotes power of the ith Gaussian component.

3. Proposed Method

For a 2D non-symmetric ID source, there are 5q unknown parameters in APDF. Compared with a 2D symmetric source, there are many more parameters to be estimated. Traditional methods such as the maximum likelihood function, COMET, and the subspace-based algorithms are difficult applications, on account of the high dimensional property of the 2D non-symmetric ID source. In this section, based on EM framework and the alternating projection principle, an iterative algorithm is proposed, which contains three steps in each iterative cycle: first, the nominal elevation angle and nominal azimuth angle of each Gaussian component are obtained utilizing feature of covariance matrices; then angular spreads are estimated by 1D searching; finally, estimation of weights are implemented though a least squares fit of theoretical and sample covariance.

3.1. Latent Variable Model

On the basis of the latent variable model [32,33], the observed received vectors x(t) and y(t) caused by source and noise can be considered as a combination of unobserved received vectors caused by each Gaussian component implied in APDF of (17) as well as noise accompanying each Gaussian component. x_i(t) and y_i(t) are supposed to be the K × 1 dimensional implied received vectors of arrays X and Y, reflecting received signal impinged by the ith Gaussian component and the noise accompanying the ith Gaussian component. According to the latent variable model, x(t) and y(t) are the observed data, which can be assumed as incomplete data; x_i(t) and y_i(t) can be assumed as complete data. The many-to-one function representing the relationship of incomplete data and complete data can be expressed as follows:

$$\{\begin{array}{l}x(t)={\displaystyle \sum _{i=1}^q{x}_i(t)}\\ y(t)={\displaystyle \sum _{i=1}^q{y}_i(t)}\end{array}.$$

(22)

Combine complete data x_i(t) and y_i(t) into:

$$z_i(t)=\left[\begin{array}{c}{x}_i(t)\\ y_i(t)\end{array}\right].$$

(23)

Thus, we have a relationship of incomplete data z(t) and complete data z_i(t) as follows:

$$z(t)={\displaystyle \sum _{i=1}^q{z}_i(t)}.$$

(24)

According to the incoherently distributed source assumption, signals from different Gaussian component are uncorrelated; so different implied received vectors z_i(t) are uncorrelated. Then covariance matrix of the incomplete data z(t), R_z is accordingly a summation of R_zi which is covariance matrix of the complete data z_i(t)._, R_z can be expressed as:

$$R_z={\displaystyle \sum _{i=1}^q{R}_{zi}}={\displaystyle \sum _{i=1}^{q}E[z_i(t)z_i^H(t)]}=\left[\begin{array}{cc}{\displaystyle \sum _{i=1}^q{R}_{xi}}& {\displaystyle \sum _{i=1}^q{R}_{xyi}}\\ {\displaystyle \sum _{i=1}^q{R}_{yxi}}& {\displaystyle \sum _{i=1}^q{R}_{yi}}\end{array}\right],$$

(25)

As z_i(t) is caused by the ith Gaussian component and its accompanied noise, covariance matrix of the complete data R_zi, can be a written as:

$$R_{zi}=E[z_i(t)z_i^H(t)]=P_i{\displaystyle \iint g(\theta ,\phi ;u_i)c(\theta ,\phi )c^H(\theta ,\phi )d\theta }d\phi +{\rho }_i{I}_{2K},$$

(26)

where ρ_i is noise power of complete data z_i(t). From Equations (22) and (23), we find that R_zi is composed of four parts—covariance matrix of the complete data x_i(t) R_xi, covariance matrix of the complete data y_i(t) R_yi, covariance matrix of the complete data x_i(t) and y_i(t) R_xyi, covariance matrix of the complete data y_i(t) and x_i(t) R_yxi, which means that:

$$R_{zi}=\left[\begin{array}{cc}{R}_{xi}& R_{xyi}\\ R_{yxi}& R_{xi}\end{array}\right]=\left[\begin{array}{cc}E[x_i(t)x_i{}^H(t)]& E[x_i(t)y_i{}^H(t)]\\ E[y_i(t)x_i{}^H(t)]& E[y_i(t)y_i{}^H(t)]\end{array}\right].$$

(27)

Denote the normalized noise-free covariance matrix of complete data z_i(t) as follows:

$$r_z(u_i)={\displaystyle \iint g(\theta ,\phi ;u_i)c(\theta ,\phi )c^H(\theta ,\phi )d\theta }d\phi =\left[\begin{array}{cc}{r}_x(u_i)& r_{xy}(u_i)\\ r_{yx}(u_i)& r_y(u_i)\end{array}\right].$$

(28)

The noise is assumed to be distributed in complete data equally so covariance matrix of z_i(t), R_zi expressed by Equation (26) can also be expressed as follows:

$$R_{zi}=P_i{r}_z(u_i)+\frac{\rho }{q}{I}_{2K}.$$

(29)

R_xi, R_yi, R_yxi and R_yxi can be expressed as follows:

$$\left[\begin{array}{cc}{R}_{xi}& R_{xyi}\\ R_{yxi}& R_{yi}\end{array}\right]=P_i\left[\begin{array}{cc}{r}_x(u_i)& r_{xy}(u_i)\\ r_{yx}(u_i)& r_y(u_i)\end{array}\right]+\frac{\rho }{q}{I}_{2K}.$$

(30)

Assuming that the angular spread of both azimuth and elevation in each Gaussian component is small, d/λ is set at 1/2; sinθsinφ and cosθsinφ can be approximated by the first term in the Taylor series expansions. Thus, elements in the normalized noise-free signal covariance matrix of Equation (30) can be expressed as follows (proof can be seen in Appendix C):

$$\{\begin{array}{l}{[r_x(u_i)]}_{kh}\approx e^{j\pi (k-h)\cos {\theta }_i\sin {\phi }_i}{e}^{-0.5\left\{{[\pi {\sigma }_{\theta i}(k-h)\sin {\theta }_i\sin {\phi }_i]}^2+{[\pi {\sigma }_{\phi i}(k-h)\cos {\theta }_i\cos {\phi }_i]}^2\right\}}\\ {[r_{xy}(u_i)]}_{kh}\approx e^{j\pi \left[(k-1)\cos {\theta }_i\sin {\phi }_i-(h-1)\sin {\theta }_i\sin {\phi }_i\right]}·{[C(u_i)]}_{kh}\\ [r_{yx}(u_i){]}_{kh}\approx e^{j\pi \left[(k-1)\sin {\theta }_i\sin {\phi }_i-(h-1)\cos {\theta }_i\sin {\phi }_i\right]}·{[D(u_i)]}_{kh}\\ [r_y(u_i){]}_{kh}\approx e^{j\pi (k-h)\sin {\theta }_i\sin {\phi }_i}{e}^{-0.5\left\{{[\pi {\sigma }_{\theta i}(k-h)\cos {\theta }_i\sin {\phi }_i]}^2+{[\pi {\sigma }_{\phi i}(k-h)\sin {\theta }_i\cos {\phi }_i]}^2\right\}}\\ {[C(u_i)]}_{kh}=e^{-0.5\left\{{\pi }^2{\sigma }^2{}_{\theta i}{[(k-1)\sin {\theta }_i\sin {\phi }_i+(h-1)\cos {\theta }_i\sin {\phi }_i]}^2+{\pi }^2{\sigma }^2{}_{\phi i}{[(k-1)\cos {\theta }_i\cos {\phi }_i-(h-1)\sin {\theta }_i\cos {\phi }_i]}^2\right\}}\\ [D(u_i){]}_{kh}=e^{-0.5\left\{{\pi }^2{\sigma }^2{}_{\theta i}{[(k-1)\cos {\theta }_i\sin {\phi }_i+(h-1)\sin {\theta }_i\sin {\phi }_i]}^2+{\pi }^2{\sigma }^2{}_{\phi i}{[(k-1)\sin {\theta }_i\cos {\phi }_i-(h-1)\cos {\theta }_i\cos {\phi }_i]}^2\right\}}\end{array},$$

(31)

where [•]_kh denotes the element of the kth row and the hth column in a matrix.

3.2. EM Algorithm

The EM algorithm is an iteration process containing two steps in turn: an expectation step (E-step) and a maximization step (M-step). The E-step serves to obtain parameters that are implied in each Gaussian component under the condition of incomplete data and the parameters values of E-step in the last EM circle. The M-step serves to update parameters based on the data from parameters obtained in the E-step, which is usually performed by maximizing the logarithm of the likelihood function.

The sample covariance matrix of incomplete data z(t), R_z can be replaced by ${\widehat{R}}_{\mathrm{z}}$ with N snapshots, which is defined as follows:

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

(32)

${\widehat{R}}_{zi}$,${\widehat{R}}_{xi}$ and ${\widehat{R}}_{yi}$ are the estimated sample covariance matrix of complete data z_i(t), x_i(t) and y_i(t). As ${\widehat{R}}_{zi}$ is a sufficient statistic of unknown parameters w_i, 𝜃_i, 𝜙_i, 𝜎_𝜃i, 𝜎_𝜙i and P, the (m+1)th E-step of the EM algorithm serves to calculate the expected value of sufficient statistics as follows (the proof can be seen in Appendix D):

$${\widehat{R}}_{zi}^{m+1}=E[{\widehat{R}}_{zi}|{\widehat{R}}_z]=R_{zi}^m{(R_z^m)}^{-1}{\widehat{R}}_z{(R_z^m)}^{-1}{R}_{zi}^m+R_{zi}^m-R_{zi}^m{(R_z^m)}^{-1}{R}_{zi}^m,$$

(33)

where the superscript m indicates the value at the mth iteration.

To simplify the calculation, we assume the same angular spread for both azimuth and elevation ${\sigma }_i={\sigma }_{\theta i}$=${\sigma }_{\varphi i}$. The M-step will minimize the negative logarithm of the likelihood function to find the optimal parameters in the (m + 1)th iteration, which can be expressed as follows:

$$w_i^{m+1},{\theta }_i^{m+1},{\phi }_i^{m+1},{\sigma }_i^{m+1},P_i^{m+1}=\underset{w_i,u_i}{\mathrm{argmin}}\left[L(u_i)=\log \left|R_{zi}\right|+\mathrm{Tr}(R_{zi}^{-1}{\widehat{R}}_{zi}^{m+1})\right](i=1,\cdots ,q),$$

(34)

where parameters to be estimated are implied in R_zi which can be expressed by Equations (29)–(31). Equation (34) means exploring the best (m + 1)th parameters of the ith Gaussian component based on the sample covariance matrix of the complete data ${\widehat{R}}_{zi}^{m+1}$. Minimizing the Equation (34) is computationally complicated because of nonlinearity and inversion of the high-dimensional matrix; therefore, solving all the unknown parameters simultaneously is impossible. According to the principle of alternating projection, parameters can be estimated successively based on other parameters that have already been solved in M-step.

According to Equation (27), the sample covariance matrices of the complete data x_i(t) and y_i(t) can be obtained as follows:

$${\widehat{\mathrm{R}}}_{\mathrm{x}\mathrm{i}}^{\mathrm{m}+1}=\left[\begin{array}{cc}{\mathrm{I}}_{\mathrm{K}}& 0_{\mathrm{K}}\\ 0_{\mathrm{K}}& 0_{\mathrm{K}}\end{array}\right]{\widehat{\mathrm{R}}}_{\mathrm{z}\mathrm{i}}^{\mathrm{m}+1},$$

(35)

$${\widehat{\mathrm{R}}}_{\mathrm{y}\mathrm{i}}^{\mathrm{m}+1}=\left[\begin{array}{cc}{0}_{\mathrm{K}}& 0_{\mathrm{K}}\\ 0_{\mathrm{K}}& {\mathrm{I}}_{\mathrm{K}}\end{array}\right]{\widehat{\mathrm{R}}}_{\mathrm{z}\mathrm{i}}^{\mathrm{m}+1},$$

(36)

where I_K is the K × K identity matrix and 0_K is K × K zero matrix.

From Equation (31), we show that the phase parts of the elements of r_x(u_i) and r_y(u_i) contain information of cosθ_isinφ_i and sinθ_isinφ_i, which means that we can obtain cosθ_isinφ_i and sinθ_isinφ_i from phase parts of R_xi and R_yi. Average phase parts of all elements in matrix ${\widehat{R}}_{xi}$ and ${\widehat{R}}_{yi}$ except the real diagonal elements, denote η_i as cosθ_isinφ_i and υ_i as sinθ_isinφ_i. Thus, we obtain:

$${\eta }_i^{m+1}=\frac{2}{K^2-K}{\displaystyle \sum _{h=1}^{k-1}{\displaystyle \sum _{k=2}^K\frac{\mathrm{angle}({[{\widehat{R}}_{xi}^{m+1}]}_{kh})}{\pi (k-h)}}},$$

(37)

$${\upsilon }_i^{m+1}=\frac{2}{K^2-K}{\displaystyle \sum _{h=1}^{k-1}{\displaystyle \sum _{k=2}^K\frac{\mathrm{angle}({[{\widehat{R}}_{yi}^{m+1}]}_{kh})}{\pi (k-h)}}},$$

(38)

where angle(•) denotes the phase of a complex number. Thus, we obtain θ_i and φ_i of the ith Gaussian component as:

$${\theta }_i^{m+1}=\arctan ({\upsilon }_i^{m+1}/{\eta }_i^{m+1}),$$

(39)

$${\phi }_i^{m+1}=\arcsin (\sqrt{{({\eta }_i^{m+1})}^2+{({\upsilon }_i^{m+1})}^2}).$$

(40)

After the nominal azimuth ${\theta }_i$ and nominal azimuth angle ${\varphi }_i$ are solved, the angular spread of the ith Gaussian components can be estimated by using the original generalized Capon estimator:

$${\sigma }_i^{m+1}=\underset{{\sigma }_i}{\mathrm{argmin}}{\mu }_{\max }[{({\widehat{R}}_{zi}^{m+1})}^{-1}{r}_z({\theta }_i^{m+1},{\phi }_i^{m+1},{\sigma }_i)],$$

(41)

where μ_max(•) represents the maximal eigenvalue of a matrix.

The least-squares fit of the theoretical and sample covariance can be expressed as:

$$U={{\Vert }{\widehat{R}}_{zi}-P_i{r}_z(u_i)-\frac{\rho }{q}{I}_{2K}{\Vert }}_{\mathrm{F}}^2.$$

(42)

Differentiating Equation (42) with respect to P_i and ρ setting the results to zero yields the following equation:

$$P_i^{m+1}=\frac{tr[{\widehat{R}}_{zi}^{m+1}{r}_z(u_i^{m+1})]-tr({\widehat{R}}_{zi}^{m+1})}{tr[r_z^2(u_i^{m+1})]-2K}.$$

(43)

At last, the weight of the ith Gaussian component can be obtained as follows:

$$w_i^{m+1}=\frac{P_i^{m+1}}{{\displaystyle \sum _i^q{P}_i^{m+1}}}.$$

(44)

After all parameters of the Gaussian component obtained in the (m + 1)th M-step, $R_{zi}^{m+1}$ can be expressed by Equations (29)–(31). $R_z^{m+1}$ is the summation of $R_{zi}^{m+1}$ according to Equation (25). Then the (m + 2)th iteration can be started.

3.3. Complexity Analysis and Comparison

Now, we analyze the computational complexity of the proposed method in comparison with COMET [26] and DISPARE [12], which are estimators for symmetric ID sources. It is noteworthy that the original algorithm in DISPARE [12] is for 1D ID sources and can be extended for 2D ID sources. COMET [26] is a method of estimation under alternating the projection algorithm framework; its computational cost mostly consists of calculation of the sample covariance matrix and the alternating projection technique with respect to cost functions. The DISPARE method estimates DOA and angular spreads through a three-dimensional (3D) spectrum-peak searching and its computation cost is mostly made of three parts: the calculation of the sample covariance matrix, the eigendecomposition of the matrix, and a 3D searching. The computation cost of the proposed method in one EM circle mainly consists of three operations: calculation of the 2K × 2K sample covariance matrix, 1D searching for angular spread, and calculation of weights. Assume that M is the EM iteration number.

Table 1. Computational complexity of different methods.

Method	Calculation of the Sample Covariance Matrix	Searching/AP Algorithm	Eigendecomposition/Calculation of Weights	Total
Proposed	o (4 NK2) + o(4 MK3)	o(8 MK3)	o(8 MK3)	o (64 MK3) + o (4 NK2)
COMET	o (4 NK2)	o (24 K3 + 12 K2)		o (24 K3 + 12 K2) + o (4 NK2)
DISPARE	o (4 NK2)	o (24 K2+6 K)	o (8K3)	o (8 K3) + o (4 NK2) + o (24 K2 + 6 K)

shows the main computational costs of three methods. From Table 1, we can clearly see if computational cost of estimation for non-symmetric distributed source is significantly higher than that of symmetric distributed source.

The stopping criterion of the EM algorithm is when all the parameters are no longer changing [30,31]. We define the iterative variation of all parameters as:

$$\mathsf{\Delta }=\frac{1}{4q}\left|\frac{{\epsilon }_{\gamma }^m-{\epsilon }_{\gamma }^m}{{\epsilon }_{\gamma }^m}\right|\gamma =1,2,\cdots ,4q,$$

(45)

where ${\epsilon }_{\gamma }^m$ represents a parameter value in the mth iteration. When $\mathsf{\Delta }$ reaches a sufficiently small value, all parameters can be considered keeping stable.

Now, our algorithm can be summarized as follows

Step 1: Determine the number q of Gaussian components. Initialize P, w_i, θ_i, φ_i, σ_i ( i = 1,2,…,q).

Step 2: Compute the sample covariance matrix of incomplete data ${\widehat{R}}_z$ using Equation (32).

Step 3: Repeat the following sub-steps for M times, which is a sufficiently large number or until the iterative variation of all parameters reaches the condition $\mathsf{\Delta }\le 0.001$.

• Compute the sample covariance matrices of complete data ${\widehat{R}}_{zi}^m$, ${\widehat{R}}_{xi}^m$ and ${\widehat{R}}_{yi}^m$ using Equations (33), (35) and (36).
• Calculate the nominal azimuth ${\theta }_i^{m+1}$ and nominal elevation ${\phi }_i^{m+1}$ from Equations (37)–(40).
• Calculate $r_z(u_i^{m+1})$ using Equation (27) and search angular spread ${\sigma }_i^{m+1}$ rough 1D search by Equation (41).
• Estimate the power of each component $P_i^{m+1}$ using Equation (39) and calculate weight of each component $w_i^{m+1}$ by Equation (44).
• Repeat sub-steps 1 to 4 for i = 1, 2, …, q.
• Superscript m = m + 1.

It is noteworthy that the distribution of a 2D non-symmetric ID source is unknown, so the true number of Gaussian components is unknown. Estimation is performed as q is an initial parameter, which needs to be set artificially. Step 1 can be considered as 0th iteration of the EM cycle. P⁰ is supposed to equal to ${\rho }^o$, which can be set at a unit power. ${\sigma }_i^0$ is set at a small value initially. $w_i^0$ can be set at 1/q. Then, ${\widehat{R}}_{zi}^0$ can be obtained according to Equations (29)–(31). ${\widehat{R}}_z^0$ can be obtained from Equation (25). As the incomplete data z(t) is observed, ${\widehat{R}}_{zi}^1$ can consequently be obtained from Equation (33).

4. Results and Discussion

In this section, we investigate the effectiveness of the proposed method though four simulation experiments. Assume that L-shape arrays have a configuration as in Figure 1 with sensors numbers K = 4 in both X and Y axis, d/λ is set at 1/2. SNR is defined as follows:

$$SNR=10\log \frac{P}{\rho }.$$

(46)

Root mean squared error (RMSE) is used to evaluate estimation performance. The RMSE of the nominal angle is defined as:

$$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}},$$

(47)

where ${\widehat{\theta }}^{\varsigma }$ and ${\widehat{\phi }}^{\varsigma }$ are the estimated nominal azimuth and estimated nominal elevation of the ID source, respectively. The superscript $\varsigma$ denotes the estimated parameters in $\varsigma$th Monte Carlo run. Mc is the number of Monte Carlo simulations, which is 500. θ and φ are the true nominal azimuth and nominal elevation, respectively.

We define the value corresponding to the maximum point of the APDF as the nominal angle of the non-symmetric ID source. Nominal angle can be obtained through partial derivative of the estimated APDF, according to the property of continuous distribution.

In addition to investigation of nominal angles, estimation of the spatial distribution should be emphasized with respect to a non-symmetric distributed source. To evaluate the estimation of spatial distribution, the RMSE of distributed function is defined as follows:

$$RMS{E}_f=\frac{1}{Mc}{\displaystyle \sum _{\varsigma }^{Mc}{\displaystyle \iint \sqrt{{\left[f(\theta ,\phi ;{\widehat{u}}^{\varsigma })-f(\theta ,\phi ;u)\right]}^2}}d\theta d\phi },$$

(48)

where $f(\theta ,\phi ;{\widehat{u}}^{\varsigma })$ is the estimated APDF in $\varsigma$th Monte Carlo run. In this section, we compare the performances of the proposed algorithm with two traditional estimators for ID sources, the COMET [26], which can be applied for any 2D arrays, and DISPARE [12], which can be extended for a 2D case. A 2D non-symmetric ID source with APDF is constructed as follows:

$$f(\theta ,\phi ;u)=0.2g(40,40,2.5)+0.2g(45,40,2.5)+0.2g(50,40,2.5)+0.3g(40,45,2.5)+0.1g(40,50,3),$$

(49)

where g(a₁,a₂,a₃) denotes a Gaussian component

$$g(a_1,a_2,a_3)=\frac{1}{2\pi a_3^2}\exp \left\{-0.5\left[{\left(\frac{\theta -a_1}{a_3}\right)}^2+{\left(\frac{\phi -a_2}{a_3}\right)}^2\right]\right\}.$$

(50)

The nominal angle of the APDF can be calculated as (40°, 44.4°). Figure 2a,b shows the constructed non-symmetric APDF. The purple region of Figure 2b is projection of the constructed non-symmetric APDF on the θ-φ plane and the color bar represents measurement of probability density function (PDF). Mark + represents the central position of each Gaussian component in Equation (49).

When estimating a 2D non-symmetric ID source, we do not know details of the non-symmetric APDF. The proposed algorithm is performed by setting initial parameters of the Gaussian mixture and iterates until the parameters no longer change. The following experiments examine the parameter variety process of each Gaussian component, different experimental conditions, number of Gaussian components, and initial positions of Gaussian components. Effectiveness of the estimated APDF can be measured though RMSE_a and RMSE_f .

In the first example, the variety processes of Gaussian components in EM iterations are investigated. We choose four Gaussian components to estimate the source and set number of snapshots N = 200 and SNR = 15dB. As the shape of APDF of a 2D non-symmetric ID source to be estimated is unknown, we can firstly get DOA estimated by the traditional method, which is defined as an assumptive nominal angle of the 2D non-symmetric ID source. To be more representative, the initial positions of Gaussian components are set uniformly around the assumptive value with same distance to the assumptive nominal angle. Thus, the nominal azimuth and nominal elevation of four Gaussian components—A, B, C and D—are set uniformly around the value (47°,48.5°), estimated by DISPARE and set at (43°, 43°), (41.5°, 52.5°), (52.5°, 44.5°) and (51°, 54°) respectively, where the distances from initial Gaussian components to the assumptive nominal angle are all 6.8°. The initial angular spreads are set at 1°. Figure 3 shows the variety processes of parameters of each component. Figure 4 shows the initial and final values of each Gaussian component. The ordered array in parentheses (${\theta }_i$,${\varphi }_i$,${\sigma }_{\theta i}$,${\sigma }_{\varphi i}$, $w_i$) of Figure 4 is the parameters of the ith Gaussian component. The beginning of the arrow represents the initial value, while the end of the arrow is the final value. The final estimated APDF is:

$$f(\theta ,\phi ;\widehat{u})=0.24g(40.1,40.3,2.7)+0.5g(39.9,46.4,3)+0.25g(46.4,40,2.8)+0.01g(39.4,45,0.9)$$

(51)

The nominal angle of the APDF is (39.5°, 45.9°), which is near the nominal angle of the given sources. RMSE_a is 1.58° and RMSE_f is 0.29. The APDF is displayed in Figure 5, which reflects the spatial non-symmetric feature of the source and is similar to the given source. It is seen that parameters will converge to certain values as a result of sufficient EM iterations. A noticeable phenomenon wherein a small weight, 0.01, developed from component D whose central position is originally far from the given source, indicates that the initial Gaussian component outside the scope of the distributed source makes hardly any contribution to the final results.

In the second example, we investigate the influence of SNR and the number of snapshots N on the performance of the proposed algorithm in comparison with COMET and DISPARE. The RMSE_a of the proposed algorithm and COMET, as well as DISPARE at different SNR and different number of snapshots N are shown, respectively, in Figure 6a,b, while RMSE_f of three methods at different SNR and N are shown, respectively, in Figure 7a,b. The EM iteration number M is set at 500. We also choose four Gaussian components to estimate the unknown source. The distances from initial Gaussian components to the assumptive nominal angle of the source are all set at 6.8°. The initial angular spreads are set at 1°. Then we randomly choose four positions with equal interval around the assumptive nominal angle (47°, 48.5°) in each Monte Carlo run. Five hundred independent Monte Carlo simulations are run to obtain the results. The number of snapshots N is set at 200 in experiments shown in Figure 6a and Figure 7a, while the SNR is set at 15dB in experiments shown in Figure 6b and Figure 7b. As the number of snapshots N or SNR increases, all algorithms provide better performance. Apparently, the proposed algorithm provides higher estimation accuracy than COMET and DISPARE algorithm with regard to RMSE_a and RMSE_f. As can be shown in Figure 6a,b, RMSE_a of COMET and DISPARE reach 4.3°, 4.3°, 7.1°, 7.3°. In Figure 7a,b, we have found that the RMSE_f of COMET and DISPARE reach 1.08, 1.1, 2.1, 2.3. As to RMSE_f, supposing in $\varsigma$th Monte Carlo trail, if the estimated APDF $f(\theta ,\phi ;{\widehat{u}}^{\varsigma })$ = 0, ${\displaystyle \iint \sqrt{{\left[f(\theta ,\phi ;{\widehat{u}}^{\varsigma })-f(\theta ,\phi ;u)\right]}^2}}d\theta d\phi$ = 1. RMSE_fs estimated by COMET and DISPARE, are big errors considering distribution of function. Therefore, Figure 7a,b show that COMET and DISPARE is invalid as to the spatial distribution estimation of the given non-symmetric distributed source on account of the large errors of RMSE_f.

In the third example, the influence of number of Gaussian components on performance is examined. The number of snapshots N = 200 and SNR = 15 dB. The initial angular spreads are set at 1°. When the performances of two or more Gaussian components are investigated, central positions of Gaussian components with equal interval around the assumptive nominal angle (47°, 48.5°) are randomly chosen in each Monte Carlo run; the distances from initial central positions to the assumptive nominal angle of the source are all set at 6.8°; 500 independent Monte Carlo simulations are run to obtain the result. As can be seen in Figure 8, the utilization of one Gaussian component provides a large error, which is an estimation considering sources as symmetric in essence. As the number of Gaussian components increase, RMSE_f and RMSE_a decrease. However, the final results of both RMSE_f and RMSE_a have little difference as the number changes from 3 to 6. Meanwhile, the convergence is markedly slower. It can be concluded that an increasing number of Gaussian components will provide a higher estimation accuracy, but the performance will not be notably improved as the number increase to a certain extent.

On the premise of not considering the cost of calculation, we can theoretically use any number of Gaussian components to fit a 2D non-symmetric ID source. In the third example, we examine the influence of different number of Gaussian components on estimation of a given source where the true number of Gaussian components is 5. It is found that the accuracy of estimation will not be significantly improved when q exceeds a certain extent. This phenomenon may be related to the shape of the given source. If the degree of asymmetry of the given 2D ID source is low, though it is composed of many Gaussian components, fewer Gaussian components can also fit the given 2D ID source well.

As the initial parameters of Gaussian components are set around assumptive value estimated by DISPARE or other methods for 2D symmetric sources, there will inevitably be Gaussian components with central positons outside the scope of the distributed source. To improve the robustness and accuracy of the algorithm, the number of Gaussian components should be set at a high level, and computing cost is also a matter of balance.

In the fourth example, we investigate the influence of initial positions to the final results. The number of snapshots N = 200 and SNR = 15dB. The initial angular spreads are set at 1°. The EM iteration is stopped under the condition Δ≤0.001. Three different assumptive nominal angles are considered. The first assumptive nominal angle is (47°,48.5°). Figure 9a shows that a circle is defined around the assumptive nominal angle (47°,48.5°). r is radius of the circle. We randomly select initial central positions of the Gaussian components, which are four points with equal interval on circle, such as the red dots in Figure 9a. 500 independent Monte Carlo simulations are run in a same circle, and then RMSE_a and RMSE_f are obtained with regard to each circle. We examine the radius of circle r changing from 0° to 20°. As shown in Figure 9b, both RMSE_a and RMSE_f change when the radius of the circle changes. The trends of the two curves are roughly the same. In general, as the radius increases, the estimation error decreases and then increases. When r is 8.5°, both RMSE_a and RMSE_f touch the bottom. The circle whose r equals 8.5° is drawn in dotted black line in Figure 9a. The second assumptive nominal angle is (45°, 44°), which is shown in Figure 10a. Setting process of initial positions is same as the first one. We examine the influence of radius of circle r changing on estimation. As shown in Figure 10b, both RMSE_a and RMSE_f decrease then increase with the radius of the circle increasing. The trends of the two curves are also roughly the same. When r is 6.8°, both RMSE_a and RMSE_f touch the bottom. The circle r equals 6.8° is drawn in dotted black line in Figure 10a. The third assumptive nominal angle is (50°, 42°), which is shown in Figure 11a. RMSE_a and RMSE_f changing with r are shown in Figure 11b. When r is 13.4°, both RMSE_a and RMSE_f touch the bottom. The circle r equals 13.4° is drawn in dotted black line in Figure 11a. The circles r, which equals 8.5°, 6.8° and 13.4° in the first, second and third trails, have common characteristics, such as initial positions of Gaussian components in these circles are within the given distributed sources with greater probability than any other circles. It is probable that estimation with initial positions of Gaussian components near the positions of Gaussian components in the given source has better accuracy than other parts, where the initial positions are far from the given source.

5. Conclusions

In this paper, we described the problem of estimating a 2D non-symmetric ID source based on L-shape arrays. The method we thus propose is developed by modeling the 2D non-symmetric APDF as a Gaussian mixture model. The estimation algorithm of a 2D non-symmetric ID source on the basis of iterative EM framework has been introduced in detail. The computational cost of a 2D non-symmetric ID source is much higher, when compared to the estimation of a 2D symmetric ID source. To evaluate the performance of estimation, we defined two indexes to reflect nominal angles and indicate spatial distribution. We investigated the parameter variety process of each Gaussian component, different SNR, number of snapshots, number of Gaussian components and initial positions of Gaussian components; the results of the numerical simulations show that the proposed method is effective for estimation of a non-symmetric ID source.