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

Authors: Tao Wu; Zhenghong Deng; Qingyue Gu; Jiwei Xu
Journal: International Journal of Antennas and Propagation
 (2018)

---

## 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.


---

## Body

## 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 [11–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–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.


## 2. Signal Model

Figure1 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.Figure 1
The L-shaped array configuration.TheK×1 output vectors of the arrays X and Y can be expressed as follows:
(1)xt=∬aθ,ϕsθ,ϕ,tdθdϕ+nxt,yt=∬bθ,ϕsθ,ϕ,tdθdϕ+nyt,where the steering vectors of the arrays X and Y can be expressed as follows:
(2)aθ,ϕ=1,ej2πdcosθsinϕ/λ,…,ej2πK−1dcosθsinϕ/λT,bθ,ϕ=1,ej2πdsinθsinϕ/λ,…,ej2πK−1dsinθsinϕ/λT.sθ,ϕ,t is the angular signal density function of the distributed source. According to the CD source assumption, sθ,ϕ,t can be expressed as
(3)sθ,ϕ,t=stfθ,ϕ;u,where st is a random process, P=Est2 is the power of the CD source, fθ,ϕ;u is the deterministic angular signal distribution function of the CD source, and u is the parameter set of the function.nxt and nyt are the K×1 additive noise vectors which can be combined into
(4)nt=nxtnyt.The noise is assumed to be uncorrelated between sensors and Gaussian white with zero mean(5)EntnHt′=ρI2Kδt−t′,EntnTt′=0∀t,t′,where ρ is the noise power, I2K 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(6)fθ,ϕ;u=∑i=1qwigiθ,ϕ;ui,where q is the number of Gaussian components. The ith Gaussian component is defined as
(7)giθ,ϕ;ui=12πσθiσϕiexp−0.5θ−θiσθi2+ϕ−ϕiσφi2,where ui=θ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 normalizefθ,ϕ;u, the weighting coefficient wi satisfies the following constraint:
(8)∑i=1qwi=1.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:
(9)xt=st∑i=1qwi∬aθ,ϕgiθ,ϕ;uidθdϕ+nxt,yt=st∑i=1qwi∬bθ,ϕgiθ,ϕ;uidθdϕ+nyt.


## 3. 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:
(10)xt=∑i=1qxit,yt=∑i=1qyit.Define the generalized steering vectors of the complete data for arraysX and Y as follows:
(11)αiui=∬aθ,ϕgiθ,ϕ;uidθdϕ,βiui=∬bθ,ϕgiθ,ϕ;uidθdϕ.Combine the incomplete dataxt and yt into
(12)zt=xtyt=st∑i=1qwiαiui∑i=1qwiβiui+nxtnyt.The noise is assumed to be distributed in the complete data equally. Combine the complete dataxit and yit into
(13)zit=xityit=stwiαiuiwiβiui+1qnxtnyt.Define two data selection matrices as follows:(14)Α=IK−1×K−1∣0K−1×1,B=0K−1×1∣IK−1×K−1,where IK−1×K−1 is the K−1×K−1 identity matrix and 0K−1×1 is the K−1×1 zero matrix. Divide X into two shifted subarrays XA and XB, and divide Y into two shifted subarrays YA and YB, which are depicted in Figure 1.The output vectors of the subarraysXA, XB, YA, and YB can be expressed as
(15)xAt,xBt,yAt,yBt=Axt,Bxt,Ayt,Byt.The complete data of output vectors of the subarraysXA, XB, YA, and YB is obtained as
(16)xAit,xBit,yAit,yBit=Axit,Bxit,Ayit,Byit.The generalized steering vectors of the complete data forXA, XB, YA, and YB can be written as
(17)αAiui,αBiui,βAiui,βBiui=Aαiui,Bαiui,Aβiui,Bβiui.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):
(18)αAiui≈ejπcosθisinϕiαBiui,βAiui≈ejπsinθisinϕiβBiui.


### 3.2. EM Algorithm

The sample covariance matrices of the complete dataR¯xi and R¯yi with N snapshots are defined as follows:
(19)R¯xi=∑t=1NxitxitH,R¯yi=∑t=1NyityitH,while the sample covariance matrix R¯zi with N snapshots is defined as
(20)R¯zi=∑t=1NxityitxityitH.The negative logarithm of the likelihood function can be simplified as(21)Lwi,ui=logRzi+trRzi−1R¯zi,where the covariance matrix of the complete data zit can be expressed as
(22)Rzi=EzitzitH=PiαiuiβiuiαiuiβiuiH+ρi.In formula (22), Pi is the power of the ith complete data. ρi=ρ/q2 is the noise power in the ith complete data. From formulas (12) and (13), zit can be denoted as
(23)zit=xityit=stwiαiuiwiβiui+1qxt−∑i=1qstwiαiuiyt−∑i=1qstwiβiui.AsR¯zi is a sufficient statistic of wi, θi, ϕi, σi, and Pi, the mth expectation step of the EM algorithm serves to calculate the expected value of the sufficient statistic as
(24)R¯zim=wim2PmαiuimβiuimαiuimβiuimH+1q2R¯x−Pm∑i=1qwimαiuim∑i=1qwimαiHuimR¯xy−Pm∑i=1qwimαiuim∑i=1qwimβiHuimR¯Hxy−Pm∑i=1qwimβiuim∑i=1qwimαiHuimR¯y−Pm∑i=1qwimβiuim∑i=1qwimβiHuim,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
(25)R¯x=1N∑t=1NxtxtH,R¯y=1N∑t=1NytytH,R¯xy=1N∑t=1NxtytH.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+1th iteration. The function (21) can be expressed as
(26)wim+1,θim+1,ϕim+1,σim+1,Pim+1=argminwi,uilogRzi+trRzi−1R¯zim,which means exploring the best m+1th parameters of the ith Gaussian component based on the sample covariance matrix of the complete data R¯zim 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 dataxit and yi̲t can be obtained as follows:
(27)R¯xim=IK×K0K×K0K×K0K×KR¯zim,R¯yim=0K×K0K×K0K×KIK×KR¯zim,where IK×K is the K×K identity matrix and 0K×K is the K×K zero matrix. By means of eigendecomposition of the covariance matrices R¯xim and R¯yim, we obtain
(28)R¯xim=cEximEximH+ρiNximNximH,R¯yim=cEyimEyimH+ρiNyimNyimH,where c is a constant, Exim and Eyim are the K-dimensional eigenvectors of R¯xim and R¯yim corresponding to the largest eigenvalues, respectively, while Nxim and Nyim are the K×K−1 matrices representing noise subspace.LetExAim and ExBim be the upper K−1 and the lower K−1 elements of Exim, respectively, corresponding to the subarrays XA and XB; EyAim and EyBim are the upper K−1 and the lower K−1 elements of Eyim, respectively, corresponding to the subarrays YA and YB. There exist relationships as follows:
(29)Exim=thefirstrowExBmi=ExAimthelastrow=thefirstrowαBiuic′=αAiuithelastrowc′,Eyim=thefirstrowEyBim=EyAimthelastrow=thefirstrowβBiuic″=βAiuithelastrowc″,where c′ and c″ are constants. Let hi=cosθisinϕi andvi=sinθisinϕi.From formula (29), we have
(30)him+1=ln11×K−1×ExBmi./ExAmijπK−1,vim+1=ln11×K−1×EyBmi./EyAmijπK−1,where I1×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+1th maximization step can be obtained from
(31)ϕim+1=arcsinhim+12+vim+12,θim+1=arctanvim+1him+1.According to the orthogonality of the subspace,σi can be obtained from one-dimensional searching
(32)σim+1=argminαiθim+1,ϕim+1,σiH,βiθim+1,ϕim+1,σiHNximNyimNximNyimHαiθim+1,ϕim+1,σiH,βiθim+1,ϕim+1,σiHH.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):
(33)αiuik=ejπk−1cosθisinϕie−0.5πσik−1sinθisinϕi2+πσik−1cosθicosϕi2,βiuik=ejπk−1sinθisinϕie−0.5πσik−1cosθisinϕi2+πσik−1sinθicosϕi2.The least-squares fits of the theoretical and sample covariance can be expressed as(34)U=R¯zi−Piriui−ρiI2KF2,where
(35)riui=αiuiβiuiαiuiβiuiH.Differentiating (34) with respect to Pi and ϱi setting the results to zero yield the following equation:
(36)ρi=trR¯zitrriuiriui∗−trriuiR¯zi∗trriui2Ktrriuiriui∗−trriuitrriui,(37)Pi=2KtrriuiR¯zi∗−trriuitrR¯zi2Ktrriuiriui∗−trriuitrriui,where superscript ·∗ denotes the conjugate operation and tr· is the trace of matrix.Let(38)γτm+1=Pim+1Pi+1m,τ=iwheni<q.FromPi=P/wi2 and formula (8), we can get
(39)wim+1=∏τ=iq−1γτm+1∏τ=1q−1γτm+1+∏τ=2q−1γτm+1+⋯+γq−1m+1+11≤i<q,1∏τ=1q−1γim+1+∏τ=2q−1γim+1+⋯+γq−1m+1+1i=q.At last, the power of the CD source can be obtained from(40)Pm+1=Pim+1wim+12.


### 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 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
(41)δ=14qεμm+1−εμmεμm,μ=1,2,…,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.Table 1
Computational complexity of different methods.
Method
Calculation of the sample covariance matrix
Eigendecomposition
Searching
Total
DSPE
o4NK2
o8K3
o24K3
o4NK2+o8K3+o24K3
TLS-ESPRIT
o4NK2
o8K3
o4NK2+o8K3
Proposed
o4MNK2
o8MK3
oM8K2+2K−12K−1
o4MNK2+o8MK3+oM8K2+2K−12K−1Now, our algorithm can be summarized as follows.Step1.
Determine the number of componentsq, and initialize P, wi, θi, ϕi, and σii=1,2,…,q.Step2.
Repeat the following substeps forM times which is a given sufficiently large number or until the iterative variation of all parameters reaches the condition δ≤0.01.
(1)
Compute the sample covariance matrices of complete dataR¯zim, R¯xim, and R¯yim using (24) and (27).
(2)
Find the eigenvectorsExim and Eyim corresponding to the largest eigenvalues through the eigendecomposition of R¯xim and R¯yim.
(3)
Calculate the nominal azimuthθim+1, nominal elevation ϕim+1, and angular spread σim+1 from (30), (31), and (32).
(4)
Estimate the power of each componentPim+1 using (37), calculate the weight of each component wim+1 using (39), and get the power of the distributed source from (40).
(5)
Repeat substeps 1 to 4 fori=1,2,…,q.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 componentsq 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.


## 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:
(10)xt=∑i=1qxit,yt=∑i=1qyit.Define the generalized steering vectors of the complete data for arraysX and Y as follows:
(11)αiui=∬aθ,ϕgiθ,ϕ;uidθdϕ,βiui=∬bθ,ϕgiθ,ϕ;uidθdϕ.Combine the incomplete dataxt and yt into
(12)zt=xtyt=st∑i=1qwiαiui∑i=1qwiβiui+nxtnyt.The noise is assumed to be distributed in the complete data equally. Combine the complete dataxit and yit into
(13)zit=xityit=stwiαiuiwiβiui+1qnxtnyt.Define two data selection matrices as follows:(14)Α=IK−1×K−1∣0K−1×1,B=0K−1×1∣IK−1×K−1,where IK−1×K−1 is the K−1×K−1 identity matrix and 0K−1×1 is the K−1×1 zero matrix. Divide X into two shifted subarrays XA and XB, and divide Y into two shifted subarrays YA and YB, which are depicted in Figure 1.The output vectors of the subarraysXA, XB, YA, and YB can be expressed as
(15)xAt,xBt,yAt,yBt=Axt,Bxt,Ayt,Byt.The complete data of output vectors of the subarraysXA, XB, YA, and YB is obtained as
(16)xAit,xBit,yAit,yBit=Axit,Bxit,Ayit,Byit.The generalized steering vectors of the complete data forXA, XB, YA, and YB can be written as
(17)αAiui,αBiui,βAiui,βBiui=Aαiui,Bαiui,Aβiui,Bβiui.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):
(18)αAiui≈ejπcosθisinϕiαBiui,βAiui≈ejπsinθisinϕiβBiui.


## 3.2. EM Algorithm

The sample covariance matrices of the complete dataR¯xi and R¯yi with N snapshots are defined as follows:
(19)R¯xi=∑t=1NxitxitH,R¯yi=∑t=1NyityitH,while the sample covariance matrix R¯zi with N snapshots is defined as
(20)R¯zi=∑t=1NxityitxityitH.The negative logarithm of the likelihood function can be simplified as(21)Lwi,ui=logRzi+trRzi−1R¯zi,where the covariance matrix of the complete data zit can be expressed as
(22)Rzi=EzitzitH=PiαiuiβiuiαiuiβiuiH+ρi.In formula (22), Pi is the power of the ith complete data. ρi=ρ/q2 is the noise power in the ith complete data. From formulas (12) and (13), zit can be denoted as
(23)zit=xityit=stwiαiuiwiβiui+1qxt−∑i=1qstwiαiuiyt−∑i=1qstwiβiui.AsR¯zi is a sufficient statistic of wi, θi, ϕi, σi, and Pi, the mth expectation step of the EM algorithm serves to calculate the expected value of the sufficient statistic as
(24)R¯zim=wim2PmαiuimβiuimαiuimβiuimH+1q2R¯x−Pm∑i=1qwimαiuim∑i=1qwimαiHuimR¯xy−Pm∑i=1qwimαiuim∑i=1qwimβiHuimR¯Hxy−Pm∑i=1qwimβiuim∑i=1qwimαiHuimR¯y−Pm∑i=1qwimβiuim∑i=1qwimβiHuim,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
(25)R¯x=1N∑t=1NxtxtH,R¯y=1N∑t=1NytytH,R¯xy=1N∑t=1NxtytH.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+1th iteration. The function (21) can be expressed as
(26)wim+1,θim+1,ϕim+1,σim+1,Pim+1=argminwi,uilogRzi+trRzi−1R¯zim,which means exploring the best m+1th parameters of the ith Gaussian component based on the sample covariance matrix of the complete data R¯zim 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 dataxit and yi̲t can be obtained as follows:
(27)R¯xim=IK×K0K×K0K×K0K×KR¯zim,R¯yim=0K×K0K×K0K×KIK×KR¯zim,where IK×K is the K×K identity matrix and 0K×K is the K×K zero matrix. By means of eigendecomposition of the covariance matrices R¯xim and R¯yim, we obtain
(28)R¯xim=cEximEximH+ρiNximNximH,R¯yim=cEyimEyimH+ρiNyimNyimH,where c is a constant, Exim and Eyim are the K-dimensional eigenvectors of R¯xim and R¯yim corresponding to the largest eigenvalues, respectively, while Nxim and Nyim are the K×K−1 matrices representing noise subspace.LetExAim and ExBim be the upper K−1 and the lower K−1 elements of Exim, respectively, corresponding to the subarrays XA and XB; EyAim and EyBim are the upper K−1 and the lower K−1 elements of Eyim, respectively, corresponding to the subarrays YA and YB. There exist relationships as follows:
(29)Exim=thefirstrowExBmi=ExAimthelastrow=thefirstrowαBiuic′=αAiuithelastrowc′,Eyim=thefirstrowEyBim=EyAimthelastrow=thefirstrowβBiuic″=βAiuithelastrowc″,where c′ and c″ are constants. Let hi=cosθisinϕi andvi=sinθisinϕi.From formula (29), we have
(30)him+1=ln11×K−1×ExBmi./ExAmijπK−1,vim+1=ln11×K−1×EyBmi./EyAmijπK−1,where I1×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+1th maximization step can be obtained from
(31)ϕim+1=arcsinhim+12+vim+12,θim+1=arctanvim+1him+1.According to the orthogonality of the subspace,σi can be obtained from one-dimensional searching
(32)σim+1=argminαiθim+1,ϕim+1,σiH,βiθim+1,ϕim+1,σiHNximNyimNximNyimHαiθim+1,ϕim+1,σiH,βiθim+1,ϕim+1,σiHH.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):
(33)αiuik=ejπk−1cosθisinϕie−0.5πσik−1sinθisinϕi2+πσik−1cosθicosϕi2,βiuik=ejπk−1sinθisinϕie−0.5πσik−1cosθisinϕi2+πσik−1sinθicosϕi2.The least-squares fits of the theoretical and sample covariance can be expressed as(34)U=R¯zi−Piriui−ρiI2KF2,where
(35)riui=αiuiβiuiαiuiβiuiH.Differentiating (34) with respect to Pi and ϱi setting the results to zero yield the following equation:
(36)ρi=trR¯zitrriuiriui∗−trriuiR¯zi∗trriui2Ktrriuiriui∗−trriuitrriui,(37)Pi=2KtrriuiR¯zi∗−trriuitrR¯zi2Ktrriuiriui∗−trriuitrriui,where superscript ·∗ denotes the conjugate operation and tr· is the trace of matrix.Let(38)γτm+1=Pim+1Pi+1m,τ=iwheni<q.FromPi=P/wi2 and formula (8), we can get
(39)wim+1=∏τ=iq−1γτm+1∏τ=1q−1γτm+1+∏τ=2q−1γτm+1+⋯+γq−1m+1+11≤i<q,1∏τ=1q−1γim+1+∏τ=2q−1γim+1+⋯+γq−1m+1+1i=q.At last, the power of the CD source can be obtained from(40)Pm+1=Pim+1wim+12.


## 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 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
(41)δ=14qεμm+1−εμmεμm,μ=1,2,…,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.Table 1
Computational complexity of different methods.
Method
Calculation of the sample covariance matrix
Eigendecomposition
Searching
Total
DSPE
o4NK2
o8K3
o24K3
o4NK2+o8K3+o24K3
TLS-ESPRIT
o4NK2
o8K3
o4NK2+o8K3
Proposed
o4MNK2
o8MK3
oM8K2+2K−12K−1
o4MNK2+o8MK3+oM8K2+2K−12K−1Now, our algorithm can be summarized as follows.Step1.
Determine the number of componentsq, and initialize P, wi, θi, ϕi, and σii=1,2,…,q.Step2.
Repeat the following substeps forM times which is a given sufficiently large number or until the iterative variation of all parameters reaches the condition δ≤0.01.
(1)
Compute the sample covariance matrices of complete dataR¯zim, R¯xim, and R¯yim using (24) and (27).
(2)
Find the eigenvectorsExim and Eyim corresponding to the largest eigenvalues through the eigendecomposition of R¯xim and R¯yim.
(3)
Calculate the nominal azimuthθim+1, nominal elevation ϕim+1, and angular spread σim+1 from (30), (31), and (32).
(4)
Estimate the power of each componentPim+1 using (37), calculate the weight of each component wim+1 using (39), and get the power of the distributed source from (40).
(5)
Repeat substeps 1 to 4 fori=1,2,…,q.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 componentsq 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.


## 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 Figure1 with K=4 sensors in both the x-axis and y-axis. d/λ is set at 1/2. SNR is defined as
(42)SNR=10logPρ.We use the root-mean-squared error (RMSE) to evaluate the estimation performance. The RMSE of the nominal angles is defined as(43)RMSEa=1Mc∑ςMcθ^ς−θ2+1Mc∑ςMcφ^ς−ϕ2,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(44)RMSEf=1Mc∑ςMc∬fθ,ϕ;u^ς−fθ,ϕ;u2dθdϕ,where fθ,ϕ;u^ς 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
(45)fθ,ϕ;u=0.2g40,40,2.5+0.2g45,40,2.5+0.2g50,40,2.5+0.3g40,45,2.5+0.1g40,50,3,where ga1,a2,a3 denotes a Gaussian component
(46)ga1,a2,a3=12πa32exp0.5θ−a1a32+ϕ−a2a32.Figure2 shows the constructed nonsymmetric angular signal distribution function which needs to be estimated.Figure 2
Probability density function (PDF) of the constructed nonsymmetric angular signal distribution.In the first example, we investigate the performance of the proposed method as well as a comparison with DSPE and TLS-ESPRIT. Figures3(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 N=200 while Figure 3(b) and 4(b) show the influence of the number of snapshots with SNR=15 dB. The EM iteration number M 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.Figure 3
(a) RMSEa estimated by the three methods versus SNR; (b) RMSEa estimated by the three methods versus number of snapshots.
(a)
(b)Figure 4
(a) RMSEf estimated by the three methods versus SNR; (b) RMSEf estimated by the three methods versus number of snapshots.
(a)
(b)Figure 5
The estimation of nonsymmetric angular signal distribution function by the proposed method.In the second example, we examine the changing processes of Gaussian components during the EM iterations. The initial nominal azimuth and nominal elevation of componentsA, 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 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 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 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.
(a)
(b)
(c)
(d)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. Figure7 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.Figure 7
(a) RMSEa estimated by different numbers of Gaussian components; (b) RMSEf estimated by different numbers of Gaussian components.
(a)
(b)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 componentA, B, C, and D 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 δ≤0.01.Figure 8
Squares of taking initial positions for each component.Figure9(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 RMSEfs change at a low level and in a small range, which means estimations maintain good accuracy with initial positions in the square.Figure 9
(a) Number of iteration versus initial position for componentA; (b) weight versus initial position for component A; (c) RMSEf versus initial position for component A.
(a)
(b)
(c)Figure10(a) shows that initial positions of component B 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 B with initial positions in this part has little contribution to estimation result; consequently, RMSEfs increase remarkably in Figure 10(c).Figure 10
(a) Number of iteration versus initial position for componentB; (b) weight versus initial position for component B; (c) RMSEf versus initial position for component B.
(a)
(b)
(c)From Figure11, 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 RMSEfs.Figure 11
(a) Number of iteration versus initial position for componentC; (b) weight versus initial position for component C; (c) RMSEf versus initial position for component C.
(a)
(b)
(c)Figure12 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 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.Figure 12
(a) Number of iteration versus initial position for componentD; (b) weight versus initial position for component D; (c) RMSEf versus initial position for component D.
(a)
(b)
(c)


## 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.



---

*Source: 5247919-2018-09-03.xml*