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

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.


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 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. 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 Sensors 2019, 19, 1226 4 of 21 array sensors. The K × 1 dimensional received signal vectors of the arrays X and Y can be expressed as follows:

Distributed Source Model
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: Sensors 2019, 19, x FOR PEER REVIEW 4 of 22 where (•) T denotes the transpose. xk(t) and yk(t) are the signal received by k th sensors in arrays X and Y, which can be expressed as follows: 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. nxk(t) and nyk(t) are noises received. αl(t) is random complex ray gain of the l th 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: where |αl(t)| 2 = σα/L, δ(•) is the Kronecker delta function, (•) * denotes the conjugate operator. Define nx(t) and ny(t), which are the K × 1 dimensional additive noise vectors of arrays X and Y; they can be written as: nx(t) and ny(t) can be combined into 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: where |α l (t)| 2 = σ α /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) = [n x1 (t), n x2 (t), · · · , n xk (t), · · · , n xK (t)], (6) n y (t) = n y1 (t), n y2 (t), · · · , n yk (t), · · · , n yK (t) .
n x (t) and n y (t) can be combined into The noise is assumed to be uncorrelated with signal and uncorrelated between sensors as well as Gaussian white with zero mean: 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: The signal is assumed to be uncorrelated with the noise; the covariance matrix of z(t) can be written as follows: 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): α l (t)e j2πd(k−1) cos θ l sin ϕ l /λ α * l (t)e −j2πd(h−1) cos θ l sin ϕ l /λ where P = σ α |s(t)| 2 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(θ, ϕ) = 1, e j2πd cos θ sin ϕ/λ , · · · e j2π(K−1)d cos θ sin ϕ/λ T , (14) b(θ, ϕ) = 1, e j2πd sin θ sin ϕ/λ , · · · e j2π(K−1)d sin θ sin ϕ/λ T .
Thus, the covariance matrices R x , R y , R xy and R yx can be expressed as follows (proof is in Appendix B): 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: where the APDF consists of q Gaussian components g(θ,ϕ; u i ) i= 1,2, . . . ,q. u i = [θ i ,φ i ,σ θi , σ φ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: 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: The summation of Gaussian mixture distributions in Equation (17) can be substituted for f (θ,ϕ). R z can be expressed respectively as follows: Equation (16) can be expressed as follows: where I K denotes the K × K identity matrix,P i = w i P denotes power of the ith Gaussian component.

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.

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: Combine complete data x i (t) and y i (t) into: Thus, we have a relationship of incomplete data z(t) and complete data z i (t) as follows: 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: 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: 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: Denote the normalized noise-free covariance matrix of complete data z i (t) as follows: 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 xi , R yi , R yxi and R yxi can be expressed as follows: 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): where [•] kh denotes the element of the kth row and the hth column in a matrix.

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 byR z with N snapshots, which is defined as follows: R zi ,R xi andR yi are the estimated sample covariance matrix of complete data z i (t), x i (t) and y i (t). AsR 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): 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 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: where parameters to be estimated are implied in R zi which can be expressed by Equations (29) According to Equation (27), the sample covariance matrices of the complete data x i (t) and y i (t) can be obtained as follows:R m+1 xi 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θ i sinϕ i and sinθ i sinϕ i , which means that we can obtain cosθ i sinϕ i and sinθ i sinϕ i from phase parts of R xi and R yi . Average phase parts of all elements in matrixR xi andR yi except the real diagonal elements, denote η i as cosθ i sinϕ i and υ i as sinθ i sinϕ i . Thus, we obtain: where angle(•) denotes the phase of a complex number. Thus, we obtain θ i and ϕ i of the ith Gaussian component as: After the nominal azimuth θ i and nominal azimuth angle φ i are solved, the angular spread of the ith Gaussian components can be estimated by using the original generalized Capon estimator: where µ max (•) represents the maximal eigenvalue of a matrix. The least-squares fit of the theoretical and sample covariance can be expressed as: Differentiating Equation (42) with respect to P i and ρ setting the results to zero yields the following equation: At last, the weight of the ith Gaussian component can be obtained as follows: After all parameters of the Gaussian component obtained in the (m + 1)th M-step, R m+1 zi can be expressed by Equations (29)- (31). R m+1 z is the summation of R m+1 zi according to Equation (25). Then the (m + 2)th iteration can be started.

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 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: where ε m γ represents a parameter value in the mth iteration. When ∆ 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 dataR 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 ∆ ≤ 0.001.

1.
Compute the sample covariance matrices of complete dataR
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 0 is supposed to equal to ρ o , which can be set at a unit power. σ 0 i is set at a small value initially. w 0 i can be set at 1/q. Then,R 0 zi can be obtained according to Equations (29)- (31).R 0 z can be obtained from Equation (25). As the incomplete data z(t) is observed,R 1 zi can consequently be obtained from Equation (33).

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: Root mean squared error (RMSE) is used to evaluate estimation performance. The RMSE of the nominal angle is defined as: whereθ ς andφ ς are the estimated nominal azimuth and estimated nominal elevation of the ID source, respectively. The superscript ς denotes the estimated parameters in ς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: where f (θ, ϕ;û ς ) is the estimated APDF in ς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: where g(a 1 ,a 2 ,a 3 ) denotes a Gaussian component 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 (θ i ,φ i ,σ θi , σ φ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: 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 RMSEa 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 RMSEf of three methods at different SNR and N are shown, respectively, in Figure 7a  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 Figures 6a and 7a, while the SNR is set at 15dB in experiments shown in Figures 6b and 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 ςth Monte Carlo trail, if the estimated APDF f (θ, ϕ;û ς ) = 0, [ f (θ, ϕ;û ς ) − f (θ, ϕ; u)] 2 dθdϕ = 1. RMSE f s 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 . and DISPARE reach 4.3°, 4.3°, 7.1°, 7.3°. In Figure 7a,b, we have found that the RMSEf of COMET and DISPARE reach 1.08, 1.1, 2.1, 2. 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 RMSEf. 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, RMSEf and RMSEa decrease. However, the final results 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 RMSEf. 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, RMSEf and RMSEa decrease. However, the final results  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.
of both RMSEf and RMSEa 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 RMSEa and RMSEf 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 RMSEa and RMSEf 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 RMSEa and RMSEf 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 RMSEa and RMSEf 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 RMSEa and RMSEf 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 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. Figure 11a. RMSEa and RMSEf changing with r are shown in Figure 11b. When r is 13.4°, both RMSEa and RMSEf 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.    Figure 11b. When r is 13.4°, both RMSEa and RMSEf 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.

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. Then [R x ] kh can be expressed as: [R x ] kh = |s(t)| 2 E L ∑ l=1 α l (t)e j2πd(k−1) cos θ l sin ϕ l /λ α * l (t)e −j2πd(h−1) cos θ l sin ϕ l /λ = σ α L |s(t)| 2 E L ∑ l=1 e j2πd(k−h) cos θ l sin ϕ l /λ (A3) As f (θ,ϕ) is the density function of scatterers of a source, we have: E e j2πd(k−h) cos θ l sin ϕ l /λ = e j2πd(k−h) cos θ sin ϕ/λ f (θ, ϕ)dθdϕ.
Considering the noise vector described by Equation (9), R x can be expressed as follows: The derivation processes of R y , R xy and R yx is similar to R x.

Appendix C
Check [r x (u i )] kh of Equation (27). Change variables (θ i + θ) for θ and (ϕ i + ϕ) for ϕ, where θ and ϕ are the small deviation of θ i and ϕ i . Thus, sinθsinϕ and cosθsinϕ can be approximated by the first term in the Taylor series expansions: [r x (u i )] kh ≈

Appendix D
Assume that complete data z i (t) and incomplete data z(t) are zero-mean complex Gaussian random vectors and their covariance matrices are R zi and R z , respectively. According to the incoherently