Abstract

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.

1. 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 two-dimensional (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 one-dimensional (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 [1115], 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 [1924] 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.

2. Signal Model

Figure 1 shows the L-shaped array configuration, which uses the xoy plane. Each linear array consists of sensor elements separated by 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 ϕ, and . is the wavelength of the impinging signal. The CD source is considered as stationary narrowband stochastic process.

The output vectors of the arrays and can be expressed as follows: where the steering vectors of the arrays and can be expressed as follows:

is the angular signal density function of the distributed source. According to the CD source assumption, can be expressed as where is a random process, is the power of the CD source, is the deterministic angular signal distribution function of the CD source, and is the parameter set of the function.

and are the additive noise vectors which can be combined into

The noise is assumed to be uncorrelated between sensors and Gaussian white with zero mean where is the noise power, is the identity matrix; superscript denotes the transpose and 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 is the number of Gaussian components. The th Gaussian component is defined as where is the parameter set with four elements denoting the nominal azimuth, nominal elevation, azimuth angular spread, and elevation angular spread, respectively.

To normalize , the weighting coefficient satisfies the following constraint:

The summation of Gaussian distributions in (6) can be substituted for in (3) and (1). The output vectors of the arrays and can be expressed as follows:

3. Proposed Method

The deterministic angular signal distribution function characterized by (6) has 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 and as follows:

Combine the incomplete data and into

The noise is assumed to be distributed in the complete data equally. Combine the complete data and into

Define two data selection matrices as follows: where is the identity matrix and is the zero matrix. Divide into two shifted subarrays and , and divide into two shifted subarrays and , which are depicted in Figure 1.

The output vectors of the subarrays , , , and can be expressed as

The complete data of output vectors of the subarrays , , , and is obtained as

The generalized steering vectors of the complete data for , , , and 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 , there exists an approximate rotational invariance relation as follows (see Appendix A):

3.2. EM Algorithm

The sample covariance matrices of the complete data and with snapshots are defined as follows: while the sample covariance matrix with snapshots is defined as

The negative logarithm of the likelihood function can be simplified as where the covariance matrix of the complete data can be expressed as

In formula (22), is the power of the th complete data. is the noise power in the th complete data. From formulas (12) and (13), can be denoted as

As is a sufficient statistic of , , , , and , the th expectation step of the EM algorithm serves to calculate the expected value of the sufficient statistic as where the superscript denotes the value of variables at the th iteration. The sample covariance matrices , , and with snapshots are defined as

To simplify the estimation, we assume that the azimuth angular spread and elevation angular spread are the same; is used to replace and for convenience. The maximization step will minimize the objective function (21) to find the optimal parameters in the th iteration. The function (21) can be expressed as which means exploring the best th parameters of the th Gaussian component based on the sample covariance matrix of the complete data obtained in the th 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 and can be obtained as follows: where is the identity matrix and is the zero matrix. By means of eigendecomposition of the covariance matrices and , we obtain where is a constant, and are the -dimensional eigenvectors of and corresponding to the largest eigenvalues, respectively, while and are the matrices representing noise subspace.

Let and be the upper and the lower elements of , respectively, corresponding to the subarrays and ; and are the upper and the lower elements of , respectively, corresponding to the subarrays and . There exist relationships as follows: where and are constants. Let and.

From formula (29), we have where is the identity matrix and (./) denotes the element-wise division operation. The nominal azimuth and nominal elevation of the complete data in the th maximization step can be obtained from

According to the orthogonality of the subspace, can be obtained from one-dimensional searching

Under the condition that the parameters nominal azimuth , nominal elevation , angular spread 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 and setting the results to zero yield the following equation: where superscript denotes the conjugate operation and is the trace of matrix.

Let

From and formula (8), we can get

At last, the power of the CD source can be obtained from

3.3. 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 L-shaped 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 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 where represents a parameter value in the th 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 sample covariance matrix, the eigendecomposition of the matrix, and 1D searching. Assume that 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 , and initialize , , , , and .

Step2. Repeat the following substeps for times which is a given sufficiently large number or until the iterative variation of all parameters reaches the condition . (1)Compute the sample covariance matrices of complete data , , and using (24) and (27).(2)Find the eigenvectors and corresponding to the largest eigenvalues through the eigendecomposition of and .(3)Calculate the nominal azimuth , nominal elevation , and angular spread from (30), (31), and (32).(4)Estimate the power of each component using (37), calculate the weight of each component using (39), and get the power of the distributed source from (40).(5)Repeat substeps 1 to 4 for .

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 is unknown; needs to be set at a reasonable value initially. In the next section, the influence of initial parameters on estimation performance will be discussed.

4. 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 sensors in both the -axis and -axis. 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 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 where denotes a Gaussian component

Figure 2 shows the constructed nonsymmetric angular signal distribution function which needs to be estimated.

In the first example, we investigate the performance of the proposed method as well as a comparison with DSPE and TLS-ESPRIT. Figures 3(a) and 3(b) display the estimation of RMSEa by the proposed algorithm, DSPE, as well as TLS-ESPRIT, while Figures 4(a) and 4(b) show the estimation of RMSEf. Figure 3(a) and 4(a) investigate the influence of SNR on performance with the number of snapshots while Figure 3(b) and 4(b) show the influence of the number of snapshots with  dB. The EM iteration number is set at 200. It can be observed that all algorithms perform better as the number of snapshots or SNR increases. Clearly, the proposed algorithm provides better performance than other estimators with regard to RMSEa and especially to RMSEf. From Figures 4(a) and 4(b), we can find that there exists a big error of RMSEf 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 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 , , , and 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 Figure 6. All parameters converge to certain values after sufficient iterations. As shown in Figure 6(d), a small weight 0.040 is developed by component , 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.

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 RMSEf and RMSEa decrease as Gaussian components increase. Both the final RMSEf and final RMSEa 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 components should be selected from a larger range; meanwhile, computing cost is a matter of balance.

In the fourth example, we examine the relationship between the initial positions of Gaussian components and estimation results. The initial positions of component , , , and are considered taking values within squares where the centers are set at (43°, 43°), (43°, 53°), (53°, 43°), and (53°, 53°). Initial weights are all set at 0.25. Each square takes 400 values uniformly. 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 RMSEf 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 RMSEf can be persuasive. The stopping criterion of the EM algorithm is the iterative variation .

Figure 9(a) shows that initial positions of component 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 to the final PDF is stable when the initial positions are set in the investigated range. Figure 9(c) shows RMSEfs change at a low level and in a small range, which means estimations maintain good accuracy with initial positions in the square.

Figure 10(a) shows that initial positions of component in the lower left part and upper right part have fewer 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 Figure 10(b), which implies that component with initial positions in this part has little contribution to estimation result; consequently, RMSEfs increase remarkably in Figure 10(c).

From Figure 11, we can find that the influence of initial positions of component on the estimation results is similar to that of component 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 , the upper right part of component has small weights and large RMSEfs.

Figure 12 shows that numbers of iterations, weights of component , and the distribution errors vary indistinctively in the whole range. No matter where the initial position is taken for component , it is far away from the distributed source. From the example of components , , , and , it can be drawn that although the initial positions closer to distributed source have faster convergence rates, the weights and RMSEfs 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 RMSEfs 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.

5. Conclusions

In this study, we have considered the problem of estimating a nonsymmetric 2D CD source under L-shaped arrays. The 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.

Appendix

A. Appendix A

Change the variables () for and for , where and are the small deviations of and . Thus, and can be approximated by the first term in the Taylor series expansions

Considering and are small variables, so we have

Thus, we can obtain following relationship

Similarly,

Consider the following relationship:

So we have

B. Appendix B

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the Chinese National Science Foundation under Grant no. 61471299.