Estimation for Two-Dimensional Nonsymmetric Coherently Distributed Source in L-Shaped Arrays

We explore the estimation of a two-dimensional (2D) nonsymmetric coherently distributed (CD) source using L-shaped arrays. Compared with a symmetric source, the modeling and estimation of a nonsymmetric source are more practical. A nonsymmetric CD source is established through modeling the deterministic angular signal distribution function as a summation of Gaussian probability density functions. Parameter estimation of the nonsymmetric distributed source is proposed under an expectation maximization (EM) framework. The proposed EM iterative calculation contains three steps in each cycle. Firstly, the nominal azimuth angles and nominal elevation angles of Gaussian components in the nonsymmetric source are obtained from the relationship of rotational invariance matrices. Then, angular spreads can be solved through one-dimensional (1D) searching based on nominal angles. Finally, the powers of Gaussian components are obtained by solving least-squares estimators. Simulations are conducted to verify the effectiveness of the nonsymmetric CD model and estimation technique.


Introduction
Although point source models are commonly used in applications such as wireless communications, radar and sonar systems, distributed source models which have considered multipath propagation and the surface features of targets tend to perform better.Position localization estimators based on distributed source models have proved to be more precise in multipath scenarios compared with point source models.With information of surface features, the distributed source models have potential application for high-resolution underwater acoustical imaging.
Sources may be classified into coherently and incoherently distributed sources [1].A coherently distributed (CD) source is defined as the signal components arriving from different angles within the range of extension that are coherent.For an incoherently distributed (ID) source, these components are uncorrelated.In this paper, the modeling and estimation of a CD source are considered.Some classical point source estimation techniques have been extended to CD sources.DSPE [1], DISPARE [2], and vec-MUSIC [3] are developed from MUSIC for distributed sources, where parameters are estimated by twodimensional (2D) spectral searching.In [4], ESPRIT is extended for distributed sources, where the TLS-ESPRIT algorithm is used to estimate nominal angles of sources firstly, and then angular spreads are obtained by onedimensional (1D) spectral searching.The authors of [5] have developed an efficient DSPE algorithm and proposed generalized beamforming estimators in [6].Generally, the parameters of CD sources are approximate solutions under the assumption of small angular spreads; some comparative studies on different methods have been proposed.The performance of two beamforming methods is analyzed in [7].
In the presence of model errors, the performance of DSPE algorithm is analyzed in [8], and the performance of MUSIC is analyzed in [9].Considering mismodeling of the spatial distribution of distributed sources, a robust estimator is proposed in [10] by means of exploiting some properties of the generalized steering vector in the case of CD sources and the covariance matrix in the case of ID sources.All these methods are based on the 1D distributed source models, which assume that signal sources are in the same plane, where distributed sources are characterized by the nominal direction of arrival (DOA) and angular spread.However, impinging signals are not generally in the same plane but in a three-dimensional space practically.The 2D distributed source models are usually described by the nominal azimuth DOA, nominal elevation DOA, azimuth angular spread, and elevation angular spread.Involving more parameters, 2D CD source estimation problem is more complicated.
Based on DSPE, several spectral searching methods for 2D CD sources have been proposed.In [11][12][13][14][15], algorithms for exponential CD sources are presented under L-shaped arrays, uniform circular arrays, and nested arrays, which transform a four-dimensional parameter searching problem into a 2D parameter searching problem.The authors of [16] have proposed an estimator for Gaussian CD sources under L-shaped arrays.
Using an array with sensors oriented along three axes of the Cartesian coordinates, the authors of [17] have extended the matrix pencil method for 2D Gaussian and Laplacian CD sources.The method needs not spectral searching, but its array structure is not beneficial to engineering realization and its concern is DOAs without angular spreads.
In [18], a sequential one-dimensional searching (SOS) algorithm is proposed for Gaussian CD sources using a pair of uniform circular arrays, which first estimates the nominal elevation by the TLS-ESPRIT method, followed by sequential one-dimensional searching to seek the nominal azimuth.As the method presented in [17], the SOS algorithm only deals with the DOAs.
Several low-complexity algorithms have been presented in [19][20][21][22][23][24] utilizing two closely spaced parallel ULAs, treble parallel ULAs, and conformal arrays.It is a common characteristic that though preliminary Taylor approximation to generalized steering vectors, the nominal elevation and azimuth are solved under ESPRIT or modified propagator.The authors of [25] have proposed a method using a centrosymmetric crossed array, where DOAs are obtained through the symmetric properties of the special array.These methods need not any spectral searching and it can deal with unknown angular distribution functions.The disadvantage of these algorithms is that angular spreads are not within consideration.
Utilizing two closely spaced parallel ULAs and L-shaped arrays, two estimators for CD noncircular signals are proposed [26,27] which exhibit better accuracy and resolution compared with circular signals.
Estimation techniques for distributed sources mentioned above are performed based on the assumption that the spatial distributions of sources are symmetric.Nevertheless, scatters are distributed irregularly or nonuniformly around targets; the distributions of scatters in practice are generally nonsymmetric.Considering complexity, people have paid less attention to nonsymmetric distributed sources.
Based on the principle that a nonsymmetric probability distribution can be composed of several symmetric distributions, some methods for 1D nonsymmetric ID sources have been presented.In [28], a nonsymmetric distribution is modeled by two Gaussian distributions.The shape of the nonsymmetric distribution would be figured via the variation of the ratio between two Gaussian distributions.In [29], the Gaussian mixture model (GMM) is employed to characterize the 1D nonsymmetric ID sources, and the expectation maximization (EM) algorithm is applied to solve the problem.
To the best of our knowledge, there are no algorithms for a 2D nonsymmetric distributed source.In this paper, we are concerned on the estimation of a 2D nonsymmetric CD source.Compared with a symmetric source, the modeling and estimation of a nonsymmetric source are more complex.Through modeling the nonsymmetric deterministic angular signal distribution function by Gaussian mixture, we have presented a parameter estimation method under an EM framework based on L-shaped arrays.In general EM frameworks, the maximization step is to maximize the likelihood function to get the best parameters which are hidden in expectation step results.In our method, the DOA parameters of Gaussian components of the 2D CD source are obtained through two approximate rotational invariance relations.By solving least-squares estimators, powers of Gaussian components are obtained.

Signal Model
Figure 1 shows the L-shaped array configuration, which uses the xoy plane.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 CD source with a nominal azimuth angle θ and a nominal elevation angle ϕ, θ ∈ −π/2, π/2 and ϕ ∈ 0, π/2 .λ is the wavelength of the impinging signal.The CD source is considered as stationary narrowband stochastic process.
The K × 1 output vectors of the arrays X and Y can be expressed as follows: x t = ∬ a θ, ϕ s θ, ϕ, t dθdϕ + n x t , y t = ∬ b θ, ϕ s θ, ϕ, t dθdϕ + n y t , 1 where the steering vectors of the arrays X and Y can be expressed as follows: The noise is assumed to be uncorrelated between sensors and Gaussian white with zero mean where ρ is the noise power, I 2K is the 2K × 2K identity matrix; superscript • T denotes the transpose and • H denotes the Hermitian transpose.δ • is the Kronecker delta function.The signal is assumed to be uncorrelated with the noise.
In this study, the deterministic angular signal distribution function of the CD source can be modeled as a summation of Gaussian distributions in order to express nonsymmetry, so the deterministic angular signal distribution function is obtained as where q is the number of Gaussian components.The ith Gaussian component is defined as

7
where u i = θ i , ϕ i , σ θi , σ ϕi is the parameter set with four elements denoting the nominal azimuth, nominal elevation, azimuth angular spread, and elevation angular spread, respectively.
To normalize f θ, ϕ ; u , the weighting coefficient w i satisfies the following constraint: The summation of Gaussian distributions in (6) can be substituted for f θ, ϕ ; u in (3) and (1).The output vectors of the arrays X and Y can be expressed as follows: w i∬ b θ, ϕ g i θ, ϕ ; u i dθdϕ + n y t 9

Proposed Method
The deterministic angular signal distribution function characterized by (6) has 5q unknown parameters.Compared with a symmetric CD source, there are much more parameters to be estimated.Due to nonlinear and high-dimensional properties, traditional methods such as maximum likelihood and spectral searching are hard to be used.The EM algorithm [30] is traditionally used to estimate Gaussian mixture parameters.Here, according to the iterative EM optimization framework, an iterative algorithm is presented for estimating the deterministic angular signal distribution function in (6).
3.1.Latent Variable Model.On the basis of the latent variable model described by the authors of [31], the incomplete data is simply the set of observations themselves.The many-to-one function, which is linear and maps the complete data to incomplete data, can be expressed as follows: Define the generalized steering vectors of the complete data for arrays X and Y as follows: Combine the incomplete data x t and y t into The L-shaped array configuration.

International Journal of Antennas and Propagation
The noise is assumed to be distributed in the complete data equally.Combine the complete data x i t and y i t into Define two data selection matrices as follows: where Divide X into two shifted subarrays X A and X B , and divide Y into two shifted subarrays Y A and Y B , which are depicted in Figure 1.
The output vectors of the subarrays X A , X B , Y A , and Y B can be expressed as The complete data of output vectors of the subarrays X A , X B , Y A , and Y B is obtained as The generalized steering vectors of the complete data for X A , X B , Y A , and Y B can be written as Assume that the angular spreads of Gaussian angular signal distribution in the complete data are small.The authors of [19] have proved that for d/λ = 1/2, there exists an approximate rotational invariance relation as follows (see Appendix A): The sample covariance matrices of the complete data R xi and R yi with N snapshots are defined as follows: while the sample covariance matrix R zi with N snapshots is defined as The negative logarithm of the likelihood function can be simplified as where the covariance matrix of the complete data z i t can be expressed as In formula (22), P i is the power of the ith complete data.ρ i = ρ/q 2 is the noise power in the ith complete data.From formulas (12) and ( 13), z i t can be denoted as As R zi is a sufficient statistic of w i , θ i , ϕ i , σ i , and P i , the m th expectation step of the EM algorithm serves to calculate the expected value of the sufficient statistic as International Journal of Antennas and Propagation where the superscript m denotes the value of variables at the mth iteration.The sample covariance matrices R x , R y , and R xy with N snapshots are defined as To simplify the estimation, we assume that the azimuth angular spread and elevation angular spread are the same; σ i is used to replace σ θi and σ ϕi for convenience.The maximization step will minimize the objective function (21) to find the optimal parameters in the m + 1 th iteration.The function ( 21) can be expressed as

26
which means exploring the best m + 1 th parameters of the ith Gaussian component based on the sample covariance matrix of the complete data R m zi obtained in the mth iteration.Minimizing (26) is too computationally complex mainly because it is nonlinear and involves a high dimensional matrix inversion.According to the rotational invariance relation, we can firstly estimate the nominal azimuth and nominal elevation of the complete data and then solve other parameters successively based on parameters that have been solved.
The sample covariance matrix of the complete data x i t and y i ̲ t can be obtained as follows: where I K×K is the K × K identity matrix and 0 K×K is the K × K zero matrix.By means of eigendecomposition of the covariance matrices R where c ′ and c ″ are constants.Let h i = cos θ i sin ϕ i and v i = sin θ i sin ϕ i .From formula (29), we have where I 1× K−1 is the 1 × K − 1 identity matrix and (./) denotes the element-wise division operation.The nominal azimuth and nominal elevation of the complete data in the m + 1 th maximization step can be obtained from According to the orthogonality of the subspace, σ i can be obtained from one-dimensional searching

International Journal of Antennas and Propagation
Under the condition that the parameters nominal azimuth θ i , nominal elevation ϕ i , angular spread σ i are solved, the generalized steering vectors of complete data can be expressed as follows (see Appendix B): The least-squares fits of the theoretical and sample covariance can be expressed as where Differentiating (34) with respect to P i and ϱ i setting the results to zero yield the following equation: where superscript • * denotes the conjugate operation and tr • is the trace of matrix.Let , τ = i when i < q 38 From P i = P/w 2 i and formula (8), we can get At last, the power of the CD source can be obtained from

Complexity Analysis and
Comparison.Now, we analyze the computational complexity of the proposed method in comparison with DSPE [16] and TLS-ESPRIT [19] which are two estimators for 2D symmetric CD sources under Lshaped arrays.The DSPE method estimates DOAs and angular spreads through two 2D spectral searches; its computation cost is mostly made of three parts: the calculation of the sample covariance matrix, the eigendecomposition of the matrix, and 2D searching.TLS-ESPRIT [19] is a method of DOA estimation; its computational cost mostly consists of the estimation and eigendecomposition of a 2K × 2K sample covariance matrix.
The stopping criterion of the EM algorithm is all parameters no longer changing [30,31].We define the iterative variation of all parameters as δ = 1 4q where ε m μ represents a parameter value in the mth iteration.When δ reaches a sufficiently small value, all parameters can be considered keeping stable.The computation cost of the proposed method in one EM circle mainly consists of three parts: the calculation of the 2K × 2K sample covariance matrix, the eigendecomposition of the matrix, and 1D searching.Assume that M is the EM iteration number.Table 1 shows the main computational cost of three methods.From Table 1, we can clearly obtain that the computational cost of the estimation for a nonsymmetric distributed source is significantly higher than that of a symmetric distributed source.Now, our algorithm can be summarized as follows.
Step1.Determine the number of components q, and initialize P, w i , θ i , ϕ i , and Step2.Repeat the following substeps for M times which is a given sufficiently large number or until the iterative variation of all parameters reaches the condition δ ≤ 0 01.(3) Calculate the nominal azimuth θ m+1 i , nominal elevation ϕ m+1 i , and angular spread σ m+1 i from ( 30), (31), and (32).α i u i k = e jπ k−1 cos θ i sin ϕ i e −0 5 πσ i k−1 sin θ i sin ϕ i 2 + πσ i k−1 cos θ i cos ϕ i 2 , International Journal of Antennas and Propagation (4) Estimate the power of each component P m+1 i using (37), calculate the weight of each component w m+1 i using (39), and get the power of the distributed source from (40).
It is noteworthy that the initial positions of the Gaussian components can be obtained by existing algorithms for a symmetric CD source and set uniformly around the guessed values.The distribution of a 2D nonsymmetric CD source is unknown so the true number of Gaussian components q is unknown; q needs to be set at a reasonable value initially.In the next section, the influence of initial parameters on estimation performance will be discussed.

Numerical Results
In this section, four simulation experiments are conducted to verify the effectiveness of the proposed technique.We consider the array geometry as depicted in Figure 1 with K = 4 sensors in both the x-axis and y-axis.d/λ is set at 1/2.SNR is defined as We use the root-mean-squared error (RMSE) to evaluate the estimation performance.The RMSE of the nominal angles is defined as where θ ς and ϕ ς are the estimated nominal azimuth and estimated nominal elevation of the CD source, respectively.The superscript ς denotes the estimated parameters from the ςth Monte Carlo run.Mc is the number of Monte Carlo simulations which is set at 500.Additionally, θ and ϕ are the true nominal azimuth and nominal elevation, respectively.
The nominal angle in the proposed algorithm is defined as the value corresponding to the maximum point of the deterministic angular signal distribution function, which can be obtained by the partial derivative of the function.
To examine the difference between the estimated and true nonsymmetric distributed source, the estimation of nominal angles is only part of the problem; in addition, the spatial distribution of the source should be emphasized compared with the estimation of a symmetric distributed source.To evaluate the estimation of the spatial distribution, the RMSE of the distributed function is defined as where f θ, ϕ ; ûς is the estimated deterministic angular signal distribution function from the ςth Monte Carlo run.In this section, we compare the performance of the proposed algorithm with DSPE [16] and TLS-ESPRIT [19] for a constructed 2D nonsymmetric CD source with the deterministic angular signal distribution function as   Clearly, the proposed algorithm provides better performance than other estimators with regard to RMSE a and especially to RMSE f .From Figures 4(a) and 4(b), we can find that there exists a big error of RMSE f by DSPE and TLS-ESPRIT.It can be concluded that utilizing traditional symmetric estimators for a nonsymmetric distributed source is invalid.Figure 5 displays the estimation of the deterministic angular     In the third example, we examine the relationship between the number of Gaussian components and estimation performance.The number of snapshots is set at 200 and SNR is 15 dB. Figure 7 shows that the error is large when there is only one Gaussian component.It can be found that RMSE f and RMSE a decrease as Gaussian components increase.Both the final RMSE f and final RMSE a change slightly as the number of Gaussian components changes from 3 to 6, and the convergence is markedly slow, which shows that increasing the number of Gaussian components will improve accuracy of estimation.Nevertheless, the effect will not be significantly improved as the number exceeds a certain extent.
The initial position parameters may be set around the guessed values; there will be initial Gaussian components outside the distributed source inevitably.To increase the robustness of the algorithm, the number of Gaussian      From Figure 11, we can find that the influence of initial positions of component C on the estimation results is similar to that of component B in respective regions.The lower left part which is close to the distributed source and the upper right part which is far from the distributed source have fewer iterations.As the similar region of component B, the upper right part of component C has small weights and large RMSE f s.
Figure 12 shows that numbers of iterations, weights of component D, and the distribution errors vary indistinctively in the whole range.No matter where the initial position is taken for component D, it is far away from the distributed source.From the example of components A, B, C, and D, it can be drawn that although the initial positions closer to distributed source have faster convergence rates, the weights and RMSE f s tend to be stable until the distances from the initial positions to the distributed source exceed a certain degree.When the distances exceed a robust range, the convergence rates increase, the weights decrease, and RMSE f s increase remarkably.The estimation of a 2D nonsymmetric CD source by the proposed algorithm shows robustness with regard to initial position of the Gaussian component.

Conclusions
In this study, we have considered the problem of estimating a nonsymmetric 2D CD source under L-shaped arrays.The  International Journal of Antennas and Propagation method we proposed is developed by modeling the nonsymmetric deterministic angular signal distribution function as Gaussian mixture.Parameters of a nonsymmetric CD source are more than those of a symmetric CD source.Parameter estimation based on an iterative EM framework has been introduced in detail.The computational cost of estimation for a nonsymmetric CD source is much higher compared with that of a symmetric CD source.For the sake of evaluating the estimation, two indicators are defined, and one reflects nominal angles and another for spatial distribution.Simulation results indicated that the proposed method is effective and robust for the estimation of a nonsymmetric CD source.e jπ k−1 cos θ i sin ϕ i +cos θ i cos ϕ i ϕ−sin θ i sin ϕ i θ dθdϕ = e jπ k−1 cos θ i sin ϕ i ∬ = e jπ k−1 cos θ i sin ϕ i e −0 5 πσ i k−1 sin θ i sin ϕ i 2 + πσ i k−1 cos θ i cos ϕ i 2 , where c is a constant, E m xi and E m yi are the K-dimensional eigenvectors of R m xi and R m yi corresponding to the largest eigenvalues, respectively, while N m xi and N m yi are the K × K − 1 matrices representing noise subspace.Let E m xAi and E m xBi be the upper K − 1 and the lower K − 1 elements of E m xi , respectively, corresponding to the subarrays X A and X B ; E m yAi and E m yBi are the upper K − 1 and the lower K − 1 elements of E m yi , respectively, corresponding to the subarrays Y A and Y B .There exist relationships as follows:

( 1 )( 2 )
Compute the sample covariance matrices of complete data R Find the eigenvectors E m xi and E m yi corresponding to the largest eigenvalues through the eigendecomposition of R m xi and R m yi .

where g a 1 ,θ − a 1 a 3 2 + ϕ − a 2 a 3 2 46Figure 2
Figure2shows the constructed nonsymmetric angular signal distribution function which needs to be estimated.

Figure 2 :
Figure 2: Probability density function (PDF) of the constructed nonsymmetric angular signal distribution.

Figure 3 :Figure 4 :Figure 5 :
Figure 3: (a) RMSE a estimated by the three methods versus SNR; (b) RMSE a estimated by the three methods versus number of snapshots.

Figure 6 :
Figure 6: (a) The changing process of central azimuths in EM iterations; (b) the changing process of central elevations in EM iterations; (c) the changing process of angular spreads in EM iterations; (d) the changing process of weights in EM iterations.

Figure 7 :Figure 8 :
Figure 7: (a) RMSE a estimated by different numbers of Gaussian components; (b) RMSE f estimated by different numbers of Gaussian components.

Figure 9 :
Figure 9: (a) Number of iteration versus initial position for component A; (b) weight versus initial position for component A; (c) RMSE f versus initial position for component A.

10
International Journal of Antennas and Propagation signal distribution function by the proposed algorithm.The result of the proposed method reflects spatial nonsymmetric distribution of the source and is closer to the given distributed source.In the second example, we examine the changing processes of Gaussian components during the EM iterations.The initial nominal azimuth and nominal elevation of components A, B, C, and D are (43 °, 43 °), (43 °, 53 °), (53 °, 43 °), and (53 °, 53 °), respectively.The initial angular spreads are set at 1 °.The number of snapshots is set at 200 and SNR is 15 dB.The changing processes of parameters in each component are shown in Figure6.All parameters converge to certain values after sufficient iterations.As shown in Figure6(d), a small weight 0.040 is developed by component D, which indicated that the nominal angles of the Gaussian component outside the scope of the distributed source make only a small contribution to the final result.

Figure 10 :
Figure 10: (a) Number of iteration versus initial position for component B; (b) weight versus initial position for component B; (c) RMSE f versus initial position for component B.

Figure 8
shows the investigated region of initial positions, where the purple region is the projection of the constructed nonsymmetric angular signal distribution function (45) on the θ-ϕ plane, and the color bar represents measurement of PDF.We investigate the influence of the initial position of a component on the number of iterations, final weight of the component, and RMSE f while initial positions of other components are set at centers.It should be noted that weights and positions of components would be changing during EM iterations, which have been shown in the second example.We focus on the sensitivity of initial positions to the final results; the weight and RMSE f can be persuasive.The stopping criterion of the EM algorithm is the iterative variation δ ≤ 0 01.

Figure 9 (
a) shows that initial positions of component A in the lower left part of the square have fewer iterations.The initial positions closer to the distributed source tend to have faster convergence rates.As shown in Figure 9(b), weights change slightly in the whole region.It can be concluded that the contribution of component A to the final PDF is stable when the initial positions are set in the investigated range.
Figure 9(c) shows RMSE f s change at a low level and in a small range, which means estimations maintain good accuracy with initial positions in the square.

Figure 10 (Figure 11 :
Figure 11: (a) Number of iteration versus initial position for component C; (b) weight versus initial position for component C; (c) RMSE f versus initial position for component C.
iterations.The lower left part has faster convergence rate due to closer distances from initial positons to the distributed source.The upper right part shows that the convergence rates increase when the distances from the initial positions to the distributed source exceed a certain degree.The weights in the upper right part decrease significantly in Figure10(b), which implies that component B with initial positions in this part has little contribution to estimation result; consequently, RMSE f s increase remarkably in Figure10(c).

Figure 12 :
Figure 12: (a) Number of iteration versus initial position for component D; (b) weight versus initial position for component D; (c) RMSE f versus initial position for component D.

Table 1 :
Computational complexity of different methods.