DOA estimation of incoherently distributed sources using importance sampling maximum likelihood

WU Tao ,DENG Zhenghong ,HU Xiaoxiang,LI Ao,and XU Jiwei

1.School of Automation,Northwestern Polytechnical University,Xi’an 710072,China;2.Equipment Management and UAV College,Air Force Engineering University,Xi’an 710051,China;3.School of Information Engineering,Xi’an University of Posts and Telecommunications,Xi’an 710061,China

Abstract: In this paper,an importance sampling maximum likelihood (ISML) estimator for direction-of-arrival (DOA) of incoherently distributed (ID) sources is proposed.Starting from the maximum likelihood estimation description of the uniform linear array(ULA),a decoupled concentrated likelihood function (CLF) is presented.A new objective function based on CLF which can obtain a closed-form solution of global maximum is constructed according to Pincus theorem.To obtain the optimal value of the objective function which is a complex high-dimensional integral,we propose an importance sampling approach based on Monte Carlo random calculation.Next,an importance function is derived,which can simplify the problem of generating random vector from a high-dimensional probability density function (PDF) to generate random variable from a one-dimensional PDF.Compared with the existing maximum likelihood (ML) algorithms for DOA estimation of ID sources,the proposed algorithm does not require initial estimates,and its performance is closer to Cramer-Rao lower bound (CRLB).The proposed algorithm performs better than the existing methods when the interval between sources to be estimated is small and in low signal to noise ratio (SNR)scenarios.

Keywords: direction-of-arrival (DOA) estimation,incoherently distributed (ID) sources,importance sampling maximum likelihood (ISML),Monte Carlo random calculation.

1.Introduction

In the field of array signal processing,such as radar,sonar,and communication,one of the most important practical needs is to estimate the direction-of-arrival (DOA) of targets.There are two ways to model targets,one being point source model,the other being the distributed source model.DOA estimation for point source targets has been relatively mature.At present,research directions mainly include estimations in special cases [1,2],applications of advanced array structure [3-5] and machine learning in estimation [6].As the distance between a receiving array and a target is shortened,spatial expansion of the target is observed.The distributed source model is more suitable to describe spatial characteristics of the target.Distributed sources can be divided into coherently distributed(CD) sources and incoherently distributed (ID) sources according to whether different scatterers of a target are correlated [7].By extending the classical methods of point sources,DOA estimation of CD sources has been well solved [8-10].However,for ID sources,the energy of signal is distributed throughout the observation space,and the covariance matrix of the signal receiver vector is full rank.If traditional algorithms of point source model are used to estimate the parameters of ID sources,the performance will decline [7,11].

DOA estimation of ID sources is one of the research hotspots in recent years.So far,DOA estimators of ID sources can be divided into five categories.One is subspace search algorithms.For example,distributed signal parameter estimator (DSPE) [7] and dispersed signal parameter estimator (DISPARE) [11] algorithms are extended from the classic point source estimator multiple signal classification (MUSIC).These technologies need spectral peak search and effective measurement with respect to dimensions of pseudo signal subspace.The subspace rotation algorithms which are widely studied at present [12-20] show great advantages in low complexity estimation.These algorithms mostly use the first order Taylor series approximation of the array manifold vectors to construct the generalized array manifold (GAM)vectors and the rotation invariant relationship between GAM vectors to estimate parameters.Since such methods are based on the estimate signal parameters via rotational invariance techniques (ESPRIT) framework,we call such methods ESPRIT-like algorithms.In such methods,[12-14] considered one dimensional ID source while[15,16] dealt with two-dimensional ID sources based on L-shaped array and circular array respectively.References [17,18] considered two-dimensional ID sources under multiple input multiple output (MIMO) systems background.Another class of low complexity algorithms is based on beamspace shift invariance structure via designing appropriate beamforming matrices such as[19,20] which have considered the noncircular ID sources.There are a class of nonparametric methods,such as beamforming and Capon spectral methods [21,22],which require high-dimensional searches and are ordinary in estimation accuracy.Covariance matching estimation techniques (COMET) methods [23-26] have also achieved good performance.However,[23,24] can only estimate a single source.In addition,initial estimates are required for iterative calculation in [25,26].The maximum likelihood (ML) estimation has advantages in estimation accuracy,whereas the global optimal solution of the ML function needs multi-dimensional grid search,the computational complexity increases geometrically with the number of sources;instead,several iterative methods have been proposed for ML approach.In theory,iterative methods cannot guarantee the convergence to the global optimal solution,and it needs to provide initial iterative values.Trump et al.[27] presented the ML estimation of ID source for the first time,which requires complex nonlinear optimization and initial estimates.An approximate ML method for joint estimation of DOA and angular spread considering multiple sources was proposed in[28],which also needs to provide initial estimates.Ghogho et al.[29] proposed an approximate ML estimation method for a simplified ID source model.Based on the same source model in [29],two approximate ML methods were proposed in [30],one of which needs initial estimates,and the other has a mediocre performance in the case of low signal to noise ratio (SNR);however,both algorithms handle single-source estimation.As it has significant advantages in estimation accuracy,ML method is widely used in DOA estimation for point sources.Aiming at joint estimation of DOA and other parameters as well as DOA estimation under special array,[31-35]achieved good performance by ML methods.References[36-38] proposed to use importance sampling to conduct Monte Carlo estimation of parameters,such as DOA and Doppler frequency.The estimation accuracy of these methods is significantly improved,and it does not depend on the initial estimation;nevertheless,the above methods can only deal with point source.

This paper presents a non-iterative ML estimation for DOA of ID sources.For ID sources,ML functions depend on DOAs and covariance matrices of receive vectors.In order to reduce parameters in process of likelihood function maximization,a concentrated likelihood function (CLF) is introduced to decouple the unnecessary parameters.To compute the optimal value of the CLF,we propose to apply the Pincus theorem [39] which can converge to the global optimal value but requires a multidimensional integral operation.Elaborating relationship between means of random vectors and the estimated values,the Monte Carlo random calculation technology[40,41] can be used to realize the multidimensional integral.The key solution then is designing an appropriate importance function as a probability density function(PDF) from which random vectors are realized.The ML method proposed in this paper is not limited to estimation for DOA of ID sources under uniform linear arrays(ULAs),and it can also be applied to different arrays by modifying the likelihood function.To show the contributions of this paper clearly,the main differences between the state-of-the-art methods and our work are listed as follows:

(i) Reference [12-20] used low complexity calculation methods,while this paper uses the ML method based on importance sampling,random calculations is utilized to find the optimal estimation.

(ii) Although [36-38] used the importance sampling method to calculate DOA,their research objects are all point sources,while this paper focuses on ID source.

(iii) Refences [27-28,30] used ML methods to estimate DOA of ID sources,they are all iterative methods and need to set initial values;the method proposed in this paper does not require initial values and can be applied to a wider range of array geometry types.The proposed method performs well in low SNR estimation and small snapshots size,especially for the signal source with relatively close distance.

This paper is organized as follows: In Section 2,we introduce the model of signal and array.In Section 3,the importance sampling maximum likelihood (ISML) estimator for DOA of ID sources is elaborated.Section 4 discusses the simulation results.Section 5 presents conclusions.

2.Signal and array model

As shown in Fig.1,assume that the number of elements in a ULA isMand the spacing isd.There areKtargets with nominal DOAsθ=[ θ1,θ2,···,θk] incident into the array,θk∈[-π/2,π/2].Different from the assumption of point source,ID source model assumes that thekth target is composed ofLkscatterers.Then the signal vector received by the array at timetcan be written as

Fig.1 Array structure of ULA with two ID sources

wheresk(t) is signal of thekth ID source,θk,l(t) is direction of thelth scatterer inkth ID source at timet,γk,l(t) is the random complex gain of thelth scatterer,n(t) is the additive white Gaussian noise with variance σ2,anda(θk,l(t)) denotes the array manifold of the ULA with respect to the point source.The ID source model assume that signal from different scatterers are uncorrelated,which means that the random complex gain oflth scatterer fromkth ID sources andl'th scatterer fromk'th ID sources [14,16,18] satisfies the following relationship:

whereδ(·) is the Kronecker delta function,and E[·] is the expectation operation.represents the total power gain of thekth source.The equation manifests that the power gain of a scatterer is 1/Lkof the total power gain of thekth source.

Letμ=2πd/λ,whereλis wavelength of signal,then themth element ofa(θk,l(t)) can be expressed as follows:

The direction of thelth scatterer within thekth ID source can be expressed as a sum of a deviation angle from nominal DOA of thekth source and nominal DOA,which can be expressed as follows:

wherea′(θk) is the partial derivative ofa(θk) atθk.According to the ID source model [7],scatterers of each source are assumed to be enormous,which may obey Gaussian,uniform,or other form distributions according to spatial geometry characteristics of targets.

According to (5),(1) can be expressed as follows:

Consequently,the receiver vector of ULA can be written as follows:

Equation (1) can be overwritten as a linear form approximately

where (·)Tmeans the transpose of a matrix,B(θ) isM×2Kdimensional and is regarded as the GAM of the ULA;g(t) is 2K×1 dimensional and is generalized signal vector.The covariance ofg(t) can be expressed as follows:

Therefore,the covariance of the received signalr(t)given by (12) can be derived as follows:

3.Proposed algorithm

3.1 ML estimation

The receiver vector of the ULA obeys a mean of 0 and a variance ofB(θ)ΛB(θ)H+σ2IMcomplex Gaussian distribution,which can be expressed as

Then the logarithmic likelihood function of the observed data can be written as

It can be seen that the unknown variables in above likelihood function includeΛand noise variance σ2in addition toθ.Reference [42] decouples the likelihood function fromΛand derives a CLF.The CLF with respect toθcan be rewritten as

Noise variance σ2can be estimated by eigendecomposition of R,taking an average of the 2K+1 toMeigenvalues of R [13].

3.2 Global maximization of the likelihood function

In order to obtain the ML estimate ofθ,we need to maximize (22) with respect toθ.Actually,(22) is high-dimensional and nonlinear,the iterative optimization approach cannot guarantee the convergence to the global optimal solution theoretically.To obtain the optimal value,the direct realization requires multi-dimensional grid search,and the computational complexity increases exponentially with the number of sources.For multivariable functions with multiple local maxima,Pincus [39] proposed a closed expression of the global maximum.

According to Pincus theorem [39],it is assumed thatxis ann-dimensional vector,andf(x) is a continuous function on a bounded domain closure inn-dimensional Euclidean space,and the global maximum off(x) is obtained at [x1,x2,···,xn].

Then,as the weight parameterρreaches infinity,the estimated value of the variablexican be obtained as

The degree of integration in (23) is the dimension of vectorx.A normalized function can be defined as follows:

In this case,Pincus theorem can be regarded as finding the expectation ofxwhich satisfies the PDF of(x).

Therefore,we can obtain the estimated DOA of thekth source by Pincus theorem as follows:

Equation (26) can be considered as the solving expectation of the random vectorθwhose pseudo PDF is(θ).Unfortunately,(θ) is a nonlinear high-dimensional function.The Monte Carlo random calculation technology[40,41] can be applied effectively.Its basic principle can be considered asis the unbiased estimate of (26),which has the following expression:

whereθqis theqth sample generated by the pseudo PDF.The sample mean converges to the optimal estimated value whenQgoes to infinity.Thereafter,the key question is how to get samples of random vectorθobeying the high-dimensional nonlinear distribution.

3.3 Importance sampling method

The importance sampling method is a random calculation approach to realize a complex integral by a designed and easily implemented pseudo PDF.

where F(x) andp(x) both satisfy the requirement of PDF,that is,they are non-negative,and the integral in domain is equal to 1.Equation (29) can be regarded as expectation ofh(x)F(x)/p(x) as the PDF ofxisp(x).p(x) is named as importance function.When the number of samples is large enough,according to the Monte Carlo random calculation principle,the expected value can be obtained by sample mean,which can be written as follows:

Now we consider designing the importance functionp(x).On one hand,the idealp(x) should be as similar to F(x)in shape as possible to reduce the variance of estimation by (30);on the other hand,the structure ofp(x)should be simple to facilitate the generation of random vectorx.Therefore,there is a trade-off between the two aspects when designingp(x).We consider the first term of CLF shown in (22) as follows:

In (31),entries inΦ=[B(θ)HB(θ)] can be written as follows:

whereA=M(M-1)/2,B=-(M-1)M(2M-1)/6.Ifk≠k',the following relationships can be proven:

ThenΦcan be approximately equal to the following expression:

Therefore,if normalization is not considered,the importance function is selected as

whereρ1is another weight parameter.

According to (34),

Considering its normalization,the importance function can be written as

Thus,eachθkcan be considered independent and has the same PDF as follows:

It can be seen that the importance sampling method can transform the problem of generating random vectors from a complex high-dimensional PDF into generating random variables from a number of one-dimensional PDFs in the context of calculating the estimated value through a large sample of random numbers.

3.4 DOA estimation

Circular mean can effectively reduce bias and computational complexity [43].Sinceθkis a cyclic random variable,we introduce the circular mean to calculate the estimated DOA.The circular mean of a random variableαwith PDF expressed by (44) and period [-π,π] is defined as

where ∠ represents the angle of a complex variable.Define variableωk=sinθk.The range of the variableωkis[-1,1].Then the circular mean form ofh(x) in (29) can be expressed as follows:

In the same way,eachωkcan be considered independent and has the same PDF as follows:

From the definition of the circular mean,according to importance sampling approach described by (29),we get the estimatedωkas follows:

ωqis theqth random vector realized byp(ω).[ωq]kis thekth variable inωq.Since the circular mean is the mean of angles of complex variable,the estimate ofωkcan be equivalent to the following expression:

y(ω) is given by dividing (46) and (47) but omitting the constantsin their denominators.This is because a complex angle after summation operation of complex variables does not change when the magnitude of each complex variable cancels out a constant term together.In this way,using the circular mean can reduce the complexity of the estimation calculation.It is worth noting that both the numerator and denominator in (51) are of an exponential form,which is prone to overflow during calculation.Hence,considering

Due to the nature of the circular mean,y′(ωq) only changes the amplitude of a complex variable and does not affect the complex angle.

3.5 Computational procedure and complexity analysis

The ML algorithm proposed in this paper can be summarized as follows:

Step 1CalculateFω(α),which is the approximate probability distribution function ofp(ω).Take a sufficient number of discrete points to approximate the probability distribution function.Set 2 000 discrete points.ωz=-1+z/1 000 (z=1,2,···,2 000).Then the probability distribution function ofp(ω) can be expressed as follows:

Step 2Generate random vectorsωfromp(ω).

Substep 2.1With the previous analysis,realize a single vectorωfrom PDFp(ω) which is equivalent to generateKrandom variables from PDFp(ω).To do this,according to the inverse probability integration method in [41],first operation is generating random numberuk(k=1,2,···,K) from [0,1] uniform distribution,then calculate

Substep 2.2Repeat Substep 2.1Qtimes to generate random vectorsωq(q=1,2,···,Q).

Step 3Estimate the value of DOA.

Substep 3.1Calculate an estimate of the circular mean(k=1,2,···,K) from (54).

Substep 3.2Estimated DOA of thekth ID source asθk=arcsinωk.The flow of the algorithm is shown in Fig.2.

Fig.2 Flow of the proposed algorithm

We analyze the computational complexity of the method proposed in comparison with representative ML methods [27,28] and a representative ESPRIT-like algorithm [14].Reference [27] solves the optimal value of likelihood function through Newton search,which is a highdimensional optimization process.Complexity of the algorithm are as follows: complexity of calculating the covariance matrix isO(MN2) and complexity of the optimization procedure isO(80LNKM3) whereLNis number of iterations of Newton search.Reference [28] separates the source power and noise power from the likelihood function,and seeks the optimal value of likelihood function through a 2Kdimensional search.The main calculation processes include calculating the covariance matrix,separating the power and noise power from the likelihood function and searching operation.Complexity of the three parts isrespectively whereLθandLσare the search points of DOA and angular spread.Reference [14] is a low complexity estimation method using the rotational invariant relationship of signal subspaces.The main calculation processes consist of calculating the covariance matrix,obtaining the signal subspace from the covariance matrix,conducting the polynomial root operation,and solving DOA from a linear equation.Complexity of the four parts isO(MN2),O(M3),O(8K3M3+8K3) andO(4MK2+2KM2) respectively.The computational cost of the proposed method mostly lies in calculating the approximate probability distribution functionFω(α),generating a large number of random vectorsωand calculating DOA.Complexity of the three parts isO(6LzNM) whereLzis number of discrete points ofFω(α),O(QK) andO(KM2+6KM) respectively.In light of above analysis,compared with [14],the computational complexity of ML algorithms is relatively high on the whole.The computational complexity of [27] depends on the number of iterationsLN,and that of [28] depends on the grid density of 2Kdimensional searchLθandLσ.Nevertheless,both methods need to provide initial estimates,consequently,the scale of iterationsLNand grid densityLθandLσdepend on the quality of the initial value.If the initial value is estimated well,the computational complexity is low,on the contrary,the computational complexity is high.The method proposed in this paper mainly depends on the numberQof random vectors and the fineness of the constructed distribution functionFω(α).

4.Numerical simulation

In this part,the performance of the proposed method is investigated by numerical simulations.Firstly,the influence of parameters of the ISML algorithm proposed in this paper is discussed.Secondly,two ML algorithms in[27,28],an ESPRIT-like algorithm [14] are compared with the proposed algorithm.The array structure of ULA is with nine elements,although the method in this paper is not limited to ULA.The SNR is defined as 10lg(1/σ2).In order to evaluate the estimation effect of different algorithms,we compare the root mean square error(RMSE) of the estimated results with the Cramer-Rao lower bound (CRLB) value in [27].

For the proposed method,ρandρ1need to be manually selected.In the following experiment,the influence of the two parameters will be examined: Two ID sources to be estimated are set as -4° and 9°,angular spreads are all 2°,Q=2 000,SNR=5 dB,the number of snapshots is set as 200.RMSE takes average of the two sources by 100 Monte Carlo simulations.The estimation results under different parameters are shown in Fig.3.As can be seen from Fig.3(a),whenρexceeds 8,the estimated result does not change significantly.However,as can be seen from Fig.3(b),with the increase ofρ1,the RMSE changes slowly at the beginning,then increases,and tends to remain unchanged when it achieves a certain extent.As parameterρtends to infinity,(27) takes the shape of aK-dimensional Dirac-delta function at the global optimal value point.There is a significant difference between other extreme points and global extreme points.In fact,when the value ofρexceeds a certain range,the global maximum points can already be separated from other extreme points.Forρ1,the requirement is satisfied as long as some similarity between the importance function and the CLF is ensured.Since the importance function is a simplification of CLF,its extreme value point is different from that of the CLF.Ifρ1is too large,the random samples will be generated around some points which deviate from the real extreme value,consequently,there will be errors in estimation.

Fig.3 RMSE versus parameter ρ and ρ1

Theoretically,the largerQis,the closer it is to the true value,but more computation is required.In this experiment,the parameters of sources to be estimated are the same as the first,SNR=5 dB,the number of snapshots is equal to 200.Fig.4 shows the RMSE estimated by 100 Monte Carlo simulations.WhenQis insufficient,the estimation is invalid,whereas the error will decrease withQincreasing.In practice,when the difference between several successive estimates is less than a given value,generating of random vectors will stop.WhenQis insufficient,there is a risk that not all sources to be estimated will be identified.Therefore,the number of unknown sources should be determined first,and then the decision conditions can be used to determine whether the random vector continues to be generated.

Fig.4 RMSE versus parameter Q with ρ=10 and ρ1=0.7

We discuss the performance of different estimators varies with SNR and number of snapshots when the sources are all distinguishable.The DOA of the two sources to be estimated are -8° and 9°.Angular spreads are all set as 2°.For ULA with nine array elements,the 3 dB beam width is approximately 12.8°,so the two sources are considered to be well apart.We chooseρ=10,ρ1=0.7,and we use 100 Monte Carlo simulations to calculate RMSE.ML algorithms [27,28] take the initial estimates as -11° and 13°.Fig.5 shows the estimated results when the number of snapshots is 100,Qis selected as 2 000,the SNR varies from -10 dB to 10 dB.Fig.5 also shows that estimation results of the proposed method are closer to CRLB compared with other methods.Especially,the method presented in this paper shows superiority under the condition of low SNR.Fig.6 shows the RMSE estimated with SNR of 0 dB,where experiment conditions remains the same as before and snapshots changes from 20 to 200.The results indicate that the proposed method shows advantages in the case of small snapshots.

Fig.5 RMSE versus SNR with number of snapshots equal to 100

Fig.6 RMSE versus number of snapshots with SNR equal to 0 dB

In next experiment,we investigate the computational complexity of different algorithms simultaneously.All the codes are written in Matlab,and the processor is the 11th Gen Inter(R) Core (TM) i7-1165G7.Table 1 shows the average estimated time for the different algorithms.It can be seen that ESPRIT-like algorithm [14] has higher computational efficiency than all other algorithms.Among the ML estimation methods,the calculation time of the proposed method is longer than algorithms [27,28].Therefore,the application scenarios of calculation need to be considered with respect to these estimators for ID sources.If the requirement for accuracy is high,the method proposed in this paper has advantages;if the requirement for real-time computing is relatively high,the algorithm proposed needs to be further studied to improve the computational efficiency.It should be noted that in estimation implementation,the algorithms [27,28] needs initial estimates,calculation of initial estimates also requires the calculation time;however,the process is not included in the statistics in this experiment.

Table 1 Estimation time of different algorithms

Since the distances of sources to be estimated have a great influence on estimation,we investigate the performance of estimators when the distance between sources is different.The DOA of the first sourceθ1is -8°.The DOA of the second source is (θ1+),whereis the interval of the two sources.Angular spreads are all 2°.SNR=5 dB,snapshots number=100,Q=4 000,ρ=10,andρ1=0.7.The initial estimate of first source of ML algorithms [27,28] is-11° and the second is 4° apart from the true value.Fig.6 shows RMSE with angles interval varying from 1° to 20°.As can be seen from Fig.7,the different estimators perform well beyond 9°.However,as the interval falls within 8°,the ESPRIT-like method [14] begins to deteriorate.Comparatively,ML algorithms [27,28] are robust to the variation of intervals between sources.It can be seen that the algorithm proposed in this paper can distinguish two sources even when the interval between two sources is 5°.Therefore,we can conclude that the method proposed has significant advantages compared with other methods in the context of near distance between two ID sources.

Fig.7 RMSE versus angles interval

In above experiments,the ML methods [27,28] require initial estimates.In order to investigate the effect of initial estimates on estimation,we discuss the relationship between initial estimates and RMSE.To get the best performance,we set SNR as 5 dB and the number of snapshots as 200.The DOAs of the two sources are -8° and 9°.Angular spreads are all 2°.Assume that the initial estimate for the first source is -11°.Fig.8 shows the RMSE curve as the initial estimate of the second source changes from 1° to 20°.It can be seen from Fig.8 that for algorithm in [27,28],if the initial estimate is significantly different from the real value,the estimated performance will be poor,whereas the algorithm in [28] is slightly better than the algorithm in [27].The ESPRIT-like algorithm does not need the initial estimation.Among ML algorithms,only the proposed algorithm can guarantee the robustness without relying on the initial estimation.

Fig.8 RMSE versus initial estimation value

5.Conclusions

In this paper,an importance sampling ML method for DOA estimation of ID sources has been developed,an objective function is constructed from the decoupled CLF which can achieve global optimal points according to the Pincus theorem.In the light of importance sampling,a selected important function is designed to realize random vectors in Monte Carlo random calculation while the circular mean concept is used to reduce the complexity.The method proposed does not require setting initial values and can be applied to a wider range of array geometry types.The proposed method performs well in low SNR estimation and small snapshots size,especially for the sources with relatively close distance.